! Program phase_stats_ST3 ! last edition 03/05/2015 ! by G. Cesana ! ! This program uses the netcdf instantaneous ST files to compute ! cloud and cloud phase statistics onto a 2deg x 2deg grid as a ! function of the height (40lvl from 0 to 19.2km, every 480m) or ! the temperature (3degC bins from -93C to 21C) ! ! Input: list of netcdf instantaneous ST files (daily broken) ! Outputs: 2 daily statistics files (height and temperature) ! ! Outputs can be concatenate using NCO to form monthly files: ! ncrcat liste_of_netcdf_files outputs.nc ! ! Compilation: ! ifort phase_stats_ST3.f90 -I/usr/include/hdf -L/usr/lib64/hdf ! -lmfhdf -ldf -ljpeg -lz -I/opt/netcdf3/ifort/include -L/opt/netcdf3/ifort/lib ! -lnetcdf -o phase_stats_ST3.e program progcf3d use netcdf implicit none character(len=200) :: filein, file2, file3, modfileout, modfileout2 integer :: err, imonth, id, ifile integer :: ilon, ilat, ialt, nprof, iz, i, itemp integer,parameter :: tempmax=38 integer,parameter :: altmax=40 integer,parameter :: lonmax=180,latmax=90 integer,parameter :: catmax=3 real,dimension(tempmax+1) :: tempmod real,dimension(tempmax) :: tempmid real,dimension(altmax) :: altmid real,dimension(altmax+1) :: alt real,dimension(latmax+1) :: latmod real,dimension(lonmax+1) :: lonmod real,dimension(latmax) :: latmid real,dimension(lonmax) :: lonmid real,dimension(:,:,:),allocatable :: cf ! 1=all, 2=ice, 3=liq, 4=attenuation real,dimension(:,:),allocatable :: temp ! instant temperature real,dimension(:),allocatable :: lat,lon ! instant lat/lon real,dimension(altmax,latmax,lonmax,catmax) :: cf3d ! Cloud fraction 2x2xL40 grid; 1=all, 2=ice, 3=liq real,dimension(tempmax,latmax,lonmax,catmax) :: cf3dtemp ! CF=f(temperature) 1=all, 2=ice, 3=liq real,dimension(altmax,latmax,lonmax) :: cf3dind real,dimension(tempmax,latmax,lonmax) :: cf3dtempind 100 format(A200) ! Initialization cf3d(:,:,:,:)=0. cf3dtemp(:,:,:,:)=0. cf3dind(:,:,:)=0. cf3dtempind(:,:,:)=0. ! Reading the daily fime list do write (*,'(a)',advance='no') 'input file name = ' read *, filein open(unit=1,file=filein,iostat=err,status='OLD') if (err==0) exit print *,'--- input file not found' enddo print *, 'file read' ! Definition of dimension vectors latmod(1)=-90 do ilat=1,latmax latmod(ilat+1)=latmod(ilat)+2 enddo latmid=latmod(1:latmax)+1 lonmod(1)=-180 do ilon=1,lonmax lonmod(ilon+1)=lonmod(ilon)+2 enddo lonmid=lonmod(1:lonmax)+1 alt(1)=0 do ilat=1,altmax alt(ilat+1)=alt(ilat)+0.48 enddo altmid=alt(1:altmax)+0.24 tempmod=(/-120.,-90.,-87.,-84.,-81.,-78.,-75.,-72.,-69.,-66.,-63.,-60.,-57.,-54.,& -51.,-48.,-45.,-42.,-39.,-36.,-33.,-30.,-27.,-24.,-21.,-18.,-15.,-12.,& -9.,-6.,-3.,0.,3.,6.,9.,12.,15.,18.,60. /) tempmid=tempmod(1:tempmax)+1.5 tempmid(1)=-91.5 ! Reading the files 888 read(1,100,end=999)file2 print *, 'file :',trim(file2) call getdim('/bdd/CFMIP/GOCCP/upload/'//trim(file2),nprof) allocate(cf(altmax+2,nprof,catmax+1),temp(altmax,nprof)) allocate(lon(nprof),lat(nprof)) call rdnc2('/bdd/CFMIP/GOCCP/upload/'//trim(file2),nprof,altmax+2,cf(:,:,1),'cloud_fraction') call rdnc2('/bdd/CFMIP/GOCCP/upload/'//trim(file2),nprof,altmax+2,cf(:,:,2),'ice_cloud_fraction') call rdnc2('/bdd/CFMIP/GOCCP/upload/'//trim(file2),nprof,altmax+2,cf(:,:,3),'water_cloud_fraction') call rdnc2('/bdd/CFMIP/GOCCP/upload/'//trim(file2),nprof,altmax+2,cf(:,:,4),'signal_attenuation') !!!!!!! you must add the temperature here that should be the averaged !temprature of the 480m pixel computed for each profile !!!!!!! !call rdnc2(trim(file2),altmax,nprof,temp,'temperature') !!! Here I use the GOCCP temperature instead of the ST temperature !!! file3='/bdd/CFMIP/GOCCP/instant_SR_CR_DR/temp/instant_SR_CR_DR_'//file2(1:19)//'ZD_day_CFMIP2_2.90.nc' print *, trim(file3) call rdnc2(trim(file3),nprof,altmax,temp,'TEMP') call rdnc('/bdd/CFMIP/GOCCP/upload/'//trim(file2),nprof,lon,'longitude') call rdnc('/bdd/CFMIP/GOCCP/upload/'//trim(file2),nprof,lat,'latitude') ! Profile loop do i=1,nprof ! Cloud Fraction 3D do ilon=1,lonmax !longitude if( (lon(i).ge.lonmod(ilon)) .and. (lon(i).lt.lonmod(ilon+1)) )then do ilat=1,latmax !latitude if ( (lat(i).ge.latmod(ilat)) .and. (lat(i).lt.latmod(ilat+1)) )then do iz = 1,altmax ! loop by vertical bins if(cf(iz,i,4).ne.0.)then ! excluding the attenuated pixels cf3d(iz,ilat,ilon,:)=cf3d(iz,ilat,ilon,:)+cf(iz,i,1:3) cf3dind(iz,ilat,ilon)=cf3dind(iz,ilat,ilon)+1. endif ! Calculation of cloudfraction as a function of the temperature do itemp=1,tempmax if ( (temp(iz,i).gt.tempmod(itemp)).and. & (temp(iz,i) .le.tempmod(itemp+1)) )then if(cf(iz,i,4).ne.0.)then ! excluding the attenuated pixels cf3dtempind(itemp,ilat,ilon)=cf3dtempind(itemp,ilat,ilon)+1. cf3dtemp(itemp,ilat,ilon,:)=cf3dtemp(itemp,ilat,ilon,:)+cf(iz,i,1:3) endif endif enddo enddo ! Altitude loop endif ! Lat enddo endif ! Lon enddo enddo ! Profile loop deallocate(cf,temp,lat,lon) print *, 'Next file' goto 888 999 continue ! Calculation of the daily statistics and attribuation of nan values (-9999) do ilon=1,lonmax do ilat=1,latmax do iz = 1, altmax ! loop by vertical bins if(cf3dind(iz,ilat,ilon)>0.)then cf3d(iz,ilat,ilon,1)=cf3d(iz,ilat,ilon,1)/cf3dind(iz,ilat,ilon) cf3d(iz,ilat,ilon,2)=cf3d(iz,ilat,ilon,2)/cf3dind(iz,ilat,ilon) cf3d(iz,ilat,ilon,3)=cf3d(iz,ilat,ilon,3)/cf3dind(iz,ilat,ilon) else cf3d(iz,ilat,ilon,:)=-9999. endif enddo enddo enddo do ilon=1,lonmax do ilat=1,latmax do itemp = 1, tempmax ! loop by vertical bins if(cf3dtempind(itemp,ilat,ilon)>0.)then cf3dtemp(itemp,ilat,ilon,1)=cf3dtemp(itemp,ilat,ilon,1)/cf3dtempind(itemp,ilat,ilon) cf3dtemp(itemp,ilat,ilon,2)=cf3dtemp(itemp,ilat,ilon,2)/cf3dtempind(itemp,ilat,ilon) cf3dtemp(itemp,ilat,ilon,3)=cf3dtemp(itemp,ilat,ilon,3)/cf3dtempind(itemp,ilat,ilon) else cf3dtemp(itemp,ilat,ilon,:)=-9999. endif enddo enddo enddo ! Recording daily files print *, 'Recording files' modfileout='/bdd/CFMIP/workspace/gcesana/ST_outputs/CF3D330m_date_CFMIP2_ST.nc' modfileout2='/bdd/CFMIP/workspace/gcesana/ST_outputs/CF3Dtemp330m_date_CFMIP2_ST.nc' call create_3d(trim(modfileout),'clcalipso','lon','lat','alt', & lonmax,latmax,altmax,cf3d(1:altmax,:,:,:),lonmid,latmid,altmid) call create_3d(trim(modfileout2),'clcalipsotmp','lon','lat','temp', & lonmax,latmax,tempmax,cf3dtemp(1:tempmax,:,:,:),lonmid,latmid,tempmid) contains ! Check the status of the nf90 command 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 ! Get the dimension of instantaneous CALIPSO files subroutine getdim(fname,fdim) use netcdf implicit none integer :: ncid, dimid, fdim character(len=*) :: fname !print *, trim(fname) call check(NF90_OPEN(trim(fname),NF90_NOWRITE,ncid)) call check(NF90_inq_dimid(ncid, 'profile', dimid)) call check(NF90_inquire_dimension(ncid, dimid, len = fdim)) end subroutine getdim ! Read 1D netcdf variable subroutine rdnc(fname,fdim,var,nvar) use netcdf implicit none character(len=*) :: nvar, fname integer :: fdim integer :: ncid,varid real,dimension(fdim) :: var 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,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 ! Record 3D variable into a netcdf file ! It includes a 1D unilimited time value to allow operation with NCO subroutine create_3d(fname,vname,ndim1,ndim2,ndim3,dim1,dim2,dim3,var5,lonmod,latmod,tempmod) use netcdf implicit none character(LEN=*) :: fname,ndim1,ndim2,ndim3,vname integer,parameter :: ndims=4 real,parameter :: nan=-9999 integer,parameter :: dim4=3 integer :: i,dim(ndims),nc_id,dim1,dim2,dim3 integer :: varid1, varid2,time_varid,varid3,varid5,varid6,varid7 integer :: dimid1,dimid2,dimid3,dimid4, time_dimid real*4,dimension(dim3,dim2,dim1,dim4) :: var5 real,dimension(dim1) :: lonmod real,dimension(dim2) :: latmod real,dimension(dim3) :: tempmod ! CREATION DU FICHIER nf90_create call check(nf90_create(fname, NF90_CLOBBER, & nc_id)) call check(nf90_put_att(nc_id, NF90_GLOBAL, 'Description', & 'ST statistics files')) call check(nf90_put_att(nc_id, NF90_GLOBAL, 'Version','phase_stats_ST3')) call check(nf90_put_att(nc_id, NF90_GLOBAL, 'Author', & 'Gregory CESANA; Cesana et al., 2015, Evaluation of 3 CALIPSO Cloud Phase Products Using In-Situ Airborne Measurements')) ! DEFINITION DES DIMENSIONS call check(nf90_def_dim(nc_id,ndim1, dim1, dimid1)) call check(nf90_def_dim(nc_id,ndim2, dim2, dimid2 )) call check(nf90_def_dim(nc_id,ndim3, dim3, dimid3)) call check(nf90_def_dim(nc_id, 'time', NF90_UNLIMITED, time_dimid)) ! DIMENSION DES MATRICES dim = (/dimid3, dimid2, dimid1, time_dimid/) ! VARIABLE 4D +time call check(nf90_def_var(nc_id, trim(vname), NF90_FLOAT, dim, varid5)) call check(nf90_put_att(nc_id, varid5, '_FillValue',nan)) call check(nf90_def_var(nc_id, trim(vname)//'_ice', NF90_FLOAT, dim, varid6)) call check(nf90_put_att(nc_id, varid6, '_FillValue',nan)) call check(nf90_def_var(nc_id, trim(vname)//'_liq', NF90_FLOAT, dim, varid7)) call check(nf90_put_att(nc_id, varid7, '_FillValue',nan)) call check(nf90_def_var(nc_id, ndim1, NF90_FLOAT, dimid1, varid1)) call check(nf90_def_var(nc_id, ndim2, NF90_FLOAT, dimid2, varid2)) call check(nf90_def_var(nc_id, ndim3, NF90_FLOAT, dimid3, varid3)) ! FIN DU MODE DEFINITION call check(nf90_enddef(nc_id)) ! ENREGISTREMENT DES VARIABLES/VECTEURS DIMENSIONS ! call check(nf90_put_var(nc_id, varid, var1)) ! call check(nf90_put_var(nc_id, varid2, var2)) call check(nf90_put_var(nc_id, varid1, lonmod)) call check(nf90_put_var(nc_id, varid2, latmod)) call check(nf90_put_var(nc_id, varid3, tempmod)) call check(nf90_put_var(nc_id, varid5, var5(:,:,:,1))) call check(nf90_put_var(nc_id, varid6, var5(:,:,:,2))) call check(nf90_put_var(nc_id, varid7, var5(:,:,:,3))) call check(nf90_close(nc_id)) end subroutine create_3d end program