program convert_pres2height use netcdf implicit none integer :: vdim, fdim1, fdim2, fdim3, fdim4,dummy,dimidsp,dimidsp2(2),& & dimsp3(3),dimsp4(4),dimsp42(4), nc_id integer :: dimidlon, dimidlat, dimidlevel, dimidlevel2,dimidtime, dimidbnds integer :: dimbndslon(2), dimbndslat(2), dimlonlat(2),varidall(37) !integer,parameter :: lonmax=320, latmax=160, levmax=46,timemax=72 character*8 :: ndim1, ndim2, ndim3, ndim4 ! nvar=name variable, ndim= name dim character*13 :: nvar character(len=200) :: filein,modfile,modfile3,fileout,file2,modfileout,modfilego,modfileoutwap character(len=200) :: taname,zgname,cliname,clwname integer,parameter :: levmax2=40, levmaxp=37 integer :: lonmax, latmax, levmax,timemax integer :: j,l,k,t, i,ii , err, imonth, id, ifile, ilat, daymax integer :: varid1,varid2,varid3,varid4,varid5,varid6,varid7,varid8,varid9,varid10 real,parameter :: R=287.058, g=9.81, mu=29 real*4 :: p0 real*4,dimension(levmax2) :: varintz,alt_mid real*4,dimension(:),allocatable :: varintp real*4,dimension(:,:,:,:),allocatable :: taz, height,heighthalf,& &ta,clw,cli,zg,& &heighthalfint,preshalf real*4,dimension(:),allocatable :: b,ap real*4,dimension(:,:,:),allocatable :: ps !!!! NEW PHASE !!! character(len=200) :: clliqname,clicename integer,parameter :: tempmax=38 integer :: itemp, dimidtemp, dimsp43(4) real*4,parameter :: nan=1.e20 real*4,dimension(:,:,:,:),allocatable :: clice,clliq,phase,phasesim real*4,dimension(:,:,:,:),allocatable :: cltemp,cltemp2,ind2,indsim2,phase2,phasesim2 real*4,dimension(tempmax) :: tempmid real*4,dimension(tempmax+1) :: tempmod real*4,dimension(tempmax,2) :: tempmod_bound ! loading the grid of the DepolSR boxes value open(20,file='/users/gcesana/grilles_lmdz/tempmod39') do itemp=1,tempmax+1 read(20,*)tempmod(itemp) enddo close(20) do itemp=1,tempmax tempmid(itemp) = (tempmod(itemp)+tempmod(itemp+1))/2 tempmod_bound(itemp,1)=tempmod(itemp); tempmod_bound(itemp,2)=tempmod(itemp+1); enddo !!!! NEW PHASE !!! alt_mid=(/0.,480.,960.,1440.,1920.,2400.,2880.,3360.,3840.,4320.,4800.,5280.,5760.,6240.,6720.,7200.,& 7680.,8160.,8640.,9120.,9600.,10080.,10560.,11040.,11520.,12000.,12480.,12960.,13440.,13920.,14400.,& 14880.,15360.,15840.,16320.,16800.,17280.,17760.,18240.,18720./) alt_mid=alt_mid+240. print *, alt_mid print *, '' !!!! NEW PHASE !!! do write (*,'(a)',advance='no') 'clice' read *, clicename if (err==0) exit enddo do write (*,'(a)',advance='no') 'clliq' read *, clliqname if (err==0) exit enddo !!!! NEW PHASE !!! do write (*,'(a)',advance='no') 'ta' read *, taname if (err==0) exit enddo do write (*,'(a)',advance='no') 'zg' read *, zgname if (err==0) exit enddo ! output file name do write (*,'(a)',advance='no') 'fileout' read *, fileout if (err==0) exit enddo ! output file name do write (*,'(a)',advance='no') 'lon lat lev time' read *, lonmax,latmax,levmax,timemax if (err==0) exit enddo print *, ' test ' !!!! NEW PHASE !!! print *, trim(clicename) print *, trim(clliqname) !!!! NEW PHASE !!! print *, trim(zgname) print *, trim(taname) print *, trim(fileout) print *, lonmax,latmax,levmax,timemax print *, "" print *, "allocating" !!! ALLOCATE VARIABLES !!! allocate(taz(lonmax,latmax,levmax2,timemax)) allocate(ta(lonmax,latmax,levmaxp,timemax)) allocate(zg(lonmax,latmax,levmaxp,timemax)) !!!! NEW PHASE !!! allocate(clliq(lonmax,latmax,levmax2,timemax)) allocate(clice(lonmax,latmax,levmax2,timemax)) !!!! NEW PHASE !!! print *, "allocation done" print *, "" print *, "initialization" taz(:,:,:,:)=0. ta(:,:,:,:)=0. zg(:,:,:,:)=0. print *, "initialization done" print *, "" !!!! NEW PHASE !!! clliq(:,:,:,:)=0. clice(:,:,:,:)=0. !!!! NEW PHASE !!! print *, "allocation done" !!! READ VARIABLES !!! call rdnc4(trim(taname),timemax,levmaxp,latmax,lonmax,ta,'ta') call rdnc4(trim(zgname),timemax,levmaxp,latmax,lonmax,zg,'geopt') !!!! NEW PHASE !!! call rdnc4(trim(clliqname),timemax,levmax2,latmax,lonmax,clliq,'clcalipso_liq') call rdnc4(trim(clicename),timemax,levmax2,latmax,lonmax,clice,'clcalipso_ice') !!!! NEW PHASE !!! print *, "allocation done" print *, "" print *, "interp ta from levp to simulator height and model pressure" do t=1,timemax do k=1,latmax do j=1,lonmax varintz(:)=1.e20 call interp(ta(j,k,:,t),varintz,zg(j,k,:,t),alt_mid,levmaxp,levmax2) taz(j,k,:,t)=varintz enddo enddo enddo print *, "interp of ta done" print *, "" print *, "deallocation" !!! DEALLOCATE !!! deallocate(ta,zg) print *, "allocating phase and ind" print *, "" allocate(phasesim2(lonmax,latmax,tempmax,timemax)) allocate(phasesim(lonmax,latmax,levmax2,timemax)) allocate(indsim2(lonmax,latmax,tempmax,timemax)) allocate(cltemp(lonmax,latmax,levmax2,timemax)) allocate(cltemp2(lonmax,latmax,tempmax,timemax)) cltemp=0. cltemp2=0. phasesim=0. phasesim2=0. indsim2=0. print *, "building phase variables" print *, "" do t=1,timemax do k=1,latmax do j=1,lonmax do l=1,levmax2 if (clliq(j,k,l,t).lt.nan) then if (clice(j,k,l,t).lt.nan) then phasesim(j,k,l,t)=(clliq(j,k,l,t)+clice(j,k,l,t)) else phasesim(j,k,l,t)=clliq(j,k,l,t) endif elseif(clice(j,k,l,t).lt.nan) then phasesim(j,k,l,t)=clice(j,k,l,t) else phasesim(j,k,l,t)=nan endif enddo enddo enddo enddo cltemp=phasesim where(phasesim.lt.nan) phasesim=clice/phasesim elsewhere phasesim=nan endwhere deallocate(clliq,clice) print *, "phase done" print *, "" print *, "phase temp" print *, "" taz=taz-273.15 do t=1,timemax do k=1,latmax do j=1,lonmax do l=1,levmax2 do itemp=1,tempmax if ( (taz(j,k,l,t).gt.tempmod(itemp)).and. & (taz(j,k,l,t) .le.tempmod(itemp+1)) )then if(phasesim2(j,k,itemp,t).lt.nan)then phasesim2(j,k,itemp,t)=phasesim2(j,k,itemp,t)+phasesim(j,k,l,t) indsim2(j,k,itemp,t)=indsim2(j,k,itemp,t)+1 cltemp2(j,k,itemp,t)=cltemp2(j,k,itemp,t)+cltemp(j,k,l,t) endif endif enddo enddo enddo enddo enddo deallocate(phasesim,taz) where(indsim2.gt.0) phasesim2=phasesim2/indsim2 cltemp2=cltemp2/indsim2 elsewhere phasesim2=nan cltemp2=nan endwhere print *, "phase2 done" print *, "" print *, "writing in file2" ! CREATION DU FICHIER nf90_create call check(nf90_create(trim(fileout), NF90_CLOBBER,nc_id)) ! DEFINITION DES DIMENSIONS call check(nf90_def_dim(nc_id,'lon', lonmax, dimidlon)) call check(nf90_def_dim(nc_id,'lat', latmax, dimidlat )) call check(nf90_def_dim(nc_id,'time', timemax, dimidtime)) !!!! NEW PHASE !!! call check(nf90_def_dim(nc_id,'temp', tempmax, dimidtemp)) !!!! NEW PHASE !!! ! call check(nf90_def_dim(nc_id,ndim6, NF90_UNLIMITED, dimid6)) ! call check(nf90_def_var(nc_id, 'time', NF90_FLOAT, dimid6, varid6)) !!!! NEW PHASE !!! dimsp43 = (/dimidlon, dimidlat, dimidtemp, dimidtime/) call check(nf90_def_var(nc_id, 'phasecaltemp', NF90_FLOAT, dimsp43, varid8)) call check(nf90_put_att(nc_id, varid8, '_FillValue',1.e20)) call check(nf90_put_att(nc_id, varid8, 'missing_value',1.e20)) call check(nf90_def_var(nc_id, 'clcaltemp', NF90_FLOAT, dimsp43, varid9)) call check(nf90_put_att(nc_id, varid9, '_FillValue',1.e20)) call check(nf90_put_att(nc_id, varid9, 'missing_value',1.e20)) !!!! NEW PHASE !!! print *, "writing in file2 definition done" ! FIN DU MODE DEFINITION call check(nf90_enddef(nc_id)) print *, "writing in file2 variables" !!! WRITE THE HEIGHT VARIABLE !!! call check(nf90_put_var(nc_id, varid8,phasesim2 )) call check(nf90_put_var(nc_id, varid9,cltemp2 )) print *, "closing file" ! CLOSE NETCDF FILE call check(nf90_close(nc_id)) print *, "deallocation" !!! DEALLOCATE !!! deallocate(phasesim2,indsim2) deallocate(cltemp,cltemp2) file2=trim(fileout(12:200)) ! CREATION DU FICHIER nf90_create call check(nf90_create(trim(file2), NF90_CLOBBER,nc_id)) ! DEFINITION DES DIMENSIONS call check(nf90_def_dim(nc_id,'lon', lonmax, dimidlon)) call check(nf90_def_dim(nc_id,'lat', latmax, dimidlat )) call check(nf90_def_dim(nc_id,'time', timemax, dimidtime)) !!!! NEW PHASE !!! call check(nf90_def_dim(nc_id,'alt', levmax2, dimidtemp)) !!!! NEW PHASE !!! ! call check(nf90_def_dim(nc_id,ndim6, NF90_UNLIMITED, dimid6)) ! call check(nf90_def_var(nc_id, 'time', NF90_FLOAT, dimid6, varid6)) !!!! NEW PHASE !!! dimsp43 = (/dimidlon, dimidlat, dimidtemp, dimidtime/) call check(nf90_def_var(nc_id, 'taz', NF90_FLOAT, dimsp43, varid8)) call check(nf90_put_att(nc_id, varid8, '_FillValue',1.e20)) call check(nf90_put_att(nc_id, varid8, 'missing_value',1.e20)) !!!! NEW PHASE !!! print *, "writing in file2 definition done" ! FIN DU MODE DEFINITION call check(nf90_enddef(nc_id)) print *, "writing in file2 variables" !!! WRITE THE HEIGHT VARIABLE !!! call check(nf90_put_var(nc_id, varid8,taz )) print *, "closing file" ! CLOSE NETCDF FILE call check(nf90_close(nc_id)) !!! SUBROUTINES !!! contains subroutine check(status) implicit none include 'netcdf.inc' integer, intent ( in) :: status if(status /= nf90_noerr) then print *, trim(nf90_strerror(status)) stop "Stopped" end if end subroutine check ! Read 2D netcdf variable subroutine rdncvar(fname,var,nvar) use netcdf implicit none character(len=*) :: nvar, fname integer :: ncid,varid real*4 :: var print *, trim(fname) call check(NF90_OPEN(trim(fname),NF90_NOWRITE,ncid)) call check(NF90_inq_varid(ncid,nvar,varid)) call check(NF90_get_var(ncid,varid,var)) call check(NF90_CLOSE(ncid)) end subroutine rdncvar ! Read 2D netcdf variable subroutine rdnc(fname,alt,var,nvar) use netcdf implicit none character(len=*) :: nvar, fname integer :: alt integer :: ncid,varid real*4,dimension(alt) :: var print *, trim(fname) call check(NF90_OPEN(trim(fname),NF90_NOWRITE,ncid)) call check(NF90_inq_varid(ncid,nvar,varid)) call check(NF90_get_var(ncid,varid,var)) call check(NF90_CLOSE(ncid)) end subroutine rdnc ! Read 2D netcdf variable subroutine rdnc2(fname,prof,alt,var,nvar) use netcdf implicit none character(len=*) :: nvar, fname integer :: alt,prof integer :: ncid,varid real*4,dimension(alt,prof) :: var print *, trim(fname) call check(NF90_OPEN(trim(fname),NF90_NOWRITE,ncid)) call check(NF90_inq_varid(ncid,nvar,varid)) call check(NF90_get_var(ncid,varid,var)) call check(NF90_CLOSE(ncid)) end subroutine rdnc2 ! Read 3D netcdf variable subroutine rdnc3(fname,lon,lat,alt,var,nvar) use netcdf implicit none character(len=*) :: nvar, fname integer :: lon,lat,alt integer :: ncid,varid real*4,dimension(alt,lat,lon) :: var print *, trim(fname) call check(NF90_OPEN(trim(fname),NF90_NOWRITE,ncid)) call check(NF90_inq_varid(ncid,nvar,varid)) call check(NF90_get_var(ncid,varid,var)) call check(NF90_CLOSE(ncid)) end subroutine rdnc3 ! Read 3D netcdf variable subroutine rdnc4(fname,lon,lat,alt,hydro,var,nvar) !call rdnc4(trim(taname),timemax,levmax,latmax,lonmax,ta,'ta') use netcdf implicit none character(len=*) :: nvar, fname integer :: lon,lat,alt,hydro integer :: ncid,varid real*4,dimension(hydro,alt,lat,lon) :: var print *, trim(fname) call check(NF90_OPEN(trim(fname),NF90_NOWRITE,ncid)) call check(NF90_inq_varid(ncid,nvar,varid)) call check(NF90_get_var(ncid,varid,var)) call check(NF90_CLOSE(ncid)) end subroutine rdnc4 ! ex : call interp(pres,pres2,altm,altl,i,it) ! !----------------------------------------------------------------------------! subroutine interp(var,varint,lev,levint,ind,indint) implicit none ! Indexes & parameters real :: a,b integer :: z,zint integer :: ind integer :: indint ! META variables real*4,dimension(indint) :: levint real*4,dimension(ind) :: lev real*4,dimension(ind) :: var real*4,dimension(indint) :: varint do z=2,ind ! exclude the nan value if( ( (var(z).lt.0.).or. & (var(z).gt.500.) ).or.( (var(z-1).lt.0.) & .or.(var(z-1).gt.500.) ) )then do zint=1,indint if ((levint(zint).le.lev(z)).and.(levint(zint).gt.lev(z-1)))then varint(zint)=1.e20 endif enddo else ! calculation of coefficient a=(var(z)-var(z-1))/(lev(z)-lev(z-1)) b=var(z)-a*lev(z) ! calculation of new interpolated variable do zint=1,indint !print *, z,zint,levint(zint),lev(z) if ((levint(zint).le.lev(z)).and.(levint(zint).gt.lev(z-1)))then varint(zint)=a*levint(zint)+b endif enddo endif enddo end subroutine interp !----------------------------------------------------------------------------! ! ex : call interp(pres,pres2,altm,altl,i,it) ! !----------------------------------------------------------------------------! subroutine interppres(var,varint,lev,levint,ind,indint) implicit none ! Indexes & parameters real :: a,b integer :: z,zint integer :: ind integer :: indint ! META variables real*4,dimension(indint) :: levint real*4,dimension(ind) :: lev real*4,dimension(ind) :: var real*4,dimension(indint) :: varint do z=2,ind ! exclude the nan value if( ( (var(z).lt.0.).or. & (var(z).gt.500.) ).or.( (var(z-1).lt.0.) & .or.(var(z-1).gt.500.) ) )then do zint=1,indint if ((levint(zint).ge.lev(z)).and.(levint(zint).lt.lev(z-1)))then varint(zint)=1.e20 endif enddo else ! calculation of coefficient a=(var(z)-var(z-1))/(lev(z)-lev(z-1)) b=var(z)-a*lev(z) ! calculation of new interpolated variable do zint=1,indint !print *, z,zint,levint(zint),lev(z) if ((levint(zint).ge.lev(z)).and.(levint(zint).lt.lev(z-1)))then varint(zint)=a*levint(zint)+b endif enddo endif enddo end subroutine interppres !----------------------------------------------------------------------------! end program