program caliporbit use netcdf implicit none character(len=200) :: filein,month2, file2, filemod, modfileout, modfileout2,daytime character(len=200) :: varfile1, varfile2 character(len=20) :: gcm, modrep, modname integer :: err, imonth, id, ifile,day, month,numfich,counter,itcs integer :: numday,numdaytmp character :: charday*3, filetmp2*1024, date*6, sds_varname*100 character*2 :: dayc,monthc,charmin,charhour real :: t2,hour integer :: nummin,numhour,t integer :: dummy1, dummy2, dummy3, dummy4 integer :: ilon, ilat, ialt, nprof, iz, i, icat, itemp,n, lontmp, lattmp integer,parameter :: lonmax=143,latmax=72,altmax=22, daymax2=1460, daymax=30 real,dimension(latmax+1) :: latmod real,dimension(lonmax+1) :: lonmod real,dimension(latmax) :: latmid real,dimension(lonmax) :: lonmid real,dimension(altmax) :: altmid real,dimension(altmax+1) :: alt real,dimension(:),allocatable :: lat,lon,vartmp1 ! instant lat/lon real,dimension(lonmax+1,latmax+1,altmax) :: QRsw,QRswind, QRlw, QRlwind real,dimension(:,:,:,:),allocatable :: tntrsw, tntrlw 100 format(A200) !------------------- Open input data file ------------------------! 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 !------------------- Output data file name ------------------------! do write (*,'(a)',advance='no') 'DATE = ' read *, date if (err==0) exit enddo !------------------- GCM's name ------------------------! do write (*,'(a)',advance='no') 'GCM = ' read *, gcm if (err==0) exit enddo allocate(tntrsw(lonmax+1,latmax+1,altmax,daymax2),tntrlw(lonmax+1,latmax+1,altmax,daymax2)) open(5,file='/home/gcesana/'//trim(gcm)//'.p') read(5,*)modrep, modname,dummy1, dummy2, dummy3, dummy4 close(5) print *, 'parameters' print *, modrep, modname QRsw(:,:,:)=0. QRswind(:,:,:)=0. QRlw(:,:,:)=0. QRlwind(:,:,:)=0. do i=1,143 lon(i+1)=lon(i)+2.5 enddo do i=1,72 lat(i+1)=lat(i)+2.5 enddo varfile1='/dav/data/archive/'//trim(modrep)//'/expt1/tntrsw/'//trim(modname)//'.tntrsw.'//date(1:4)//'0101-'//date(1:4)//'1231.nc' varfile2='/dav/data/archive/'//trim(modrep)//'/expt1/tntrlw/'//trim(modname)//'.tntrlw.'//date(1:4)//'0101-'//date(1:4)//'1231.nc' print *, trim(varfile1) print *, trim(varfile2) call rdnc4(trim(varfile1),lonmax+1,latmax+1,altmax,daymax2,tntrsw,'tntrsw') call rdnc4(trim(varfile2),lonmax+1,latmax+1,altmax,daymax2,tntrlw,'tntrlw') print *, 'model files read' ! Initialization of box and numfich numfich=0 numdaytmp=0 counter=0 !------------------ Read the list of hdf files -------------------! 888 read(1,100,end=999)file2 print *, 'Processing with ',trim(file2) numfich=numfich+1 print *, 'lecture du fichier numero ',numfich filetmp2=trim(file2) charday=filetmp2(38:40) charhour=filetmp2(41:42) charmin=filetmp2(43:44) print *, charday,' ',charhour,charmin ! listfile a calculer write(numday,'(a3)')charday write(numhour,'(a2)')charhour write(nummin,'(a2)')charmin if(numday.gt.numdaytmp)then counter=counter+1 endif print *, 'Read the CloudSat file data ...' !-------------------- Get the number of profiles ----------------! call sdsreadprof(filetmp2,itcs) !----------------------- Allocation -----------------------------! allocate(lat(itcs),lon(itcs)) !,RH(it,2)) lat(:)=0. lon(:)=0. sds_varname='Latitude' print *, trim(sds_varname) call sdsreadeos1(vartmp1,filetmp2,trim(sds_varname)) lat=vartmp1 deallocate(vartmp1) sds_varname='Longitude' print *, trim(sds_varname) call sdsreadeos1(vartmp1,filetmp2,trim(sds_varname)) lon=vartmp1 deallocate(vartmp1) hour=1/60*nummin+numhour hour=nint(hour) print *, hour t2=(numday-1)*4+(hour/6+1) t2=nint(t2) t=t2 lontmp=1 lattmp=1 if(t.le.daymax2)then do i=1,itcs ! 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).le.latmod(ilat)) .and. (lat(i).gt.latmod(ilat+1)) )then if( (i.gt.1).and. ((lontmp.ne.ilon).or.(lattmp.ne.ilat)) )then lontmp=ilon lattmp=ilat !print *, ilon,ilat do iz=1,22 if (tntrsw(lontmp,lattmp,iz,t).lt.105) then QRsw(lontmp,lattmp,iz)=QRsw(lontmp,lattmp,iz)+tntrsw(lontmp,lattmp,iz,t) QRswind(lontmp,lattmp,iz)=QRswind(lontmp,lattmp,iz)+1 endif if (tntrlw(lontmp,lattmp,iz,t).lt.105) then QRlw(lontmp,lattmp,iz)=QRlw(lontmp,lattmp,iz)+tntrlw(lontmp,lattmp,iz,t) QRlwind(lontmp,lattmp,iz)=QRlwind(lontmp,lattmp,iz)+1 endif enddo endif endif enddo endif enddo enddo endif deallocate(lat,lon) !print *, 'Next file' goto 888 999 continue deallocate(tntrsw,tntrlw) where(QRswind.gt.0) QRsw=QRsw/QRswind elsewhere QRsw=-9999. endwhere where(QRlwind.gt.0) QRlw=QRlw/QRlwind elsewhere QRlw=-9999. endwhere !do day=1,30 !do ilon=1,lonmax+1 !do ilat=1,latmax+1 !do iz = 1, altmax ! loop by vertical bins ! if(cf3dind(ilon,ilat,iz,day)>0.)then ! cf3d(ilon,ilat,iz,day)=cf3d(ilon,ilat,iz,day)/cf3dind(ilon,ilat,iz,day) ! else ! cf3d(ilon,ilat,iz,day)=-9999. ! endif !enddo !enddo !enddo !enddo ! Recording daily files print *, 'Recording files' modfileout2='/home/gcesana/'//trim(modname)//'_QRlw_'//date//'_CloudSat.nc' modfileout='/home/gcesana/'//trim(modname)//'_QRsw_'//date//'_CloudSat.nc' altmid(:)=0. !varfile1='/dav/data/archive/'//trim(modrep)//'/expt1/tntrsw/'//trim(modname)//'.tntrsw.'//date(1:4)//'0101-'//date(1:4)//'1231.nc' call create_3d(trim(modfileout),'QRsw','lon','lat','alt', & lonmax+1,latmax+1,altmax,QRsw,lonmod,latmod,altmid) call create_3d(trim(modfileout2),'QRlw','lon','lat','alt', & lonmax+1,latmax+1,altmax,QRsw,lonmod,latmod,altmid) 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 subroutine getdim(fname,vname,fdim) use netcdf implicit none integer :: ncid, dimid, fdim character(len=*) :: fname,vname !print *, trim(fname) call check(NF90_OPEN(trim(fname),NF90_NOWRITE,ncid)) call check(NF90_inq_dimid(ncid, trim(vname), dimid)) call check(NF90_inquire_dimension(ncid, dimid, len = fdim)) end subroutine getdim subroutine rdnc4(fname,lon,lat,alt,day,var,nvar) use netcdf implicit none character(len=*) :: nvar, fname integer :: lon,lat,alt,day integer :: ncid,varid real,dimension(lon,lat,alt,day) :: 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 rdnc4 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 character*3 :: ndim4='day' integer,parameter :: ndims=4 real,parameter :: nan=-9999 integer,parameter :: dim4=1 integer :: i,dim(ndims),nc_id,dim1,dim2,dim3 integer :: varid1, varid2,tiime_varid,varid3,varid4,varid5 integer :: dimid1,dimid2,dimid3,dimid4, time_dimid real*4,dimension(dim1,dim2,dim3,dim4) :: var5 real,dimension(dim1) :: lonmod real,dimension(dim2) :: latmod real,dimension(dim3) :: altmod real,dimension(dim4) :: 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', & 'LMDZ CALIOP Orbit files')) call check(nf90_put_att(nc_id, NF90_GLOBAL, 'Version','caliporbit2')) 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,ndim4, NF90_UNLIMITED, dimid4)) ! call check(nf90_def_dim(nc_id, 'time', NF90_UNLIMITED, time_dimid)) ! DIMENSION DES MATRICES dim = (/dimid1, dimid2, dimid3, dimid4/) ! 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, 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)) call check(nf90_def_var(nc_id, ndim4, NF90_FLOAT, dimid4, varid4)) ! FIN DU MODE DEFINITION call check(nf90_enddef(nc_id)) ! ENREGISTREMENT DES VARIABLES/VECTEURS DIMENSIONS 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, altmod)) call check(nf90_put_var(nc_id, varid4, tempmod)) call check(nf90_put_var(nc_id, varid5, var5)) call check(nf90_close(nc_id)) end subroutine create_3d subroutine create_2d(fname,vname,ndim1,ndim2,dim1,dim2,var5,lonmod,latmod) use netcdf implicit none character(LEN=*) :: fname,ndim1,ndim2,vname integer,parameter :: ndims=3 real,parameter :: nan=-9999 character*3 :: ndim3='day' integer,parameter :: dim3=30,dim4=4 integer :: i,dim(ndims),nc_id,dim1,dim2 integer :: varid1, varid2,time_varid,varid3,varid4,varid5,varid6,varid7,varid8 integer :: dimid1,dimid2,dimid3,dimid4, time_dimid real*4,dimension(dim1,dim2,dim3,dim4) :: var5 real,dimension(dim1) :: lonmod real,dimension(dim2) :: latmod real,dimension(dim3) :: tempmod tempmod(:)=0. do i=1,30 tempmod(i)=i enddo ! CREATION DU FICHIER nf90_create call check(nf90_create(fname, NF90_CLOBBER, & nc_id)) call check(nf90_put_att(nc_id, NF90_GLOBAL, 'Description', & 'LMDZ CALIOP Orbit files')) call check(nf90_put_att(nc_id, NF90_GLOBAL, 'Version','caliporbit2')) 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, NF90_UNLIMITED, dimid3 )) ! call check(nf90_def_dim(nc_id, 'time', NF90_UNLIMITED, time_dimid)) ! DIMENSION DES MATRICES dim = (/dimid1,dimid2, dimid3/) ! VARIABLE 4D +time call check(nf90_def_var(nc_id, trim(vname)//'tot', NF90_FLOAT, dim, varid5)) call check(nf90_put_att(nc_id, varid5, '_FillValue',nan)) call check(nf90_def_var(nc_id, trim(vname)//'high', NF90_FLOAT, dim, varid6)) call check(nf90_put_att(nc_id, varid6, '_FillValue',nan)) call check(nf90_def_var(nc_id, trim(vname)//'mid', NF90_FLOAT, dim, varid7)) call check(nf90_put_att(nc_id, varid7, '_FillValue',nan)) call check(nf90_def_var(nc_id, trim(vname)//'low', NF90_FLOAT, dim, varid8)) call check(nf90_put_att(nc_id, varid8, '_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, 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_put_var(nc_id, varid8, var5(:,:,:,4))) call check(nf90_close(nc_id)) end subroutine create_2d subroutine sdsreadprof(filename,iprof) implicit none ! use hdf include "hdf.f90" include "dffunc.f90" character (len=128) ::filename,swathname, fieldname integer (kind=4) :: fid,swid,status,astat integer (kind=4) :: swopen,swattach,swdetach,swclose integer (kind=4) :: swfldinfo,swrdfld,swrdattr integer (kind=4) :: rank,dims(7),datatype,i,j character (len=255) :: dimlist integer :: iprof swathname = "2B-FLXHR-LIDAR" fieldname = "Latitude" fid = swopen(filename, DFACC_READ) if (fid .eq. -1) then print *,'CANNOT OPEN FILE',fid stop endif !======= attach to the swath swid = swattach(fid,swathname) if (swid .eq. -1) then print *,'CANNOT ATTACH TO SWATH',swid stop endif !======= get field information status = swfldinfo(swid,fieldname,rank,dims,datatype,dimlist) if (status .ne. 0) then print *,'CANNOT GET FIELD INFO',status stop endif !====== get the number of profiles iprof = dims(1) !======= detach from the swath status = swdetach(swid) if (status .ne. 0) then print *,'CANNOT DETACH FROM SWATH',status stop endif !======= close the file status = swclose(fid) if (status .ne. 0) then print *,'CANNOT CLOSE FILE',status stop endif end subroutine sdsreadprof subroutine sdsreadeos1(var,filename,fieldname) implicit none ! use hdf include "hdf.f90" include "dffunc.f90" real,allocatable :: var(:) character(len=*) :: fieldname real, allocatable :: tempfld(:) character (len=128) ::filename,swathname integer (kind=4) :: fid,swid,status,astat integer (kind=4) :: swopen,swattach,swdetach,swclose integer (kind=4) :: swfldinfo,swrdfld,swrdattr integer (kind=4) :: rank,dims(7),datatype,i,j character (len=255) :: dimlist integer (kind=4), allocatable :: start(:),stride(:) swathname = "2B-FLXHR-LIDAR" fid = swopen(filename, DFACC_READ) if (fid .eq. -1) then print *,'CANNOT OPEN FILE',fid stop endif !======= attach to the swath swid = swattach(fid,swathname) if (swid .eq. -1) then print *,'CANNOT ATTACH TO SWATH',swid stop endif !======= get field information status = swfldinfo(swid,fieldname,rank,dims,datatype,dimlist) if (status .ne. 0) then print *,'CANNOT GET FIELD INFO',status stop endif !======= allocate the field to be read in allocate(tempfld(dims(1)),var(dims(1)),stat=astat) if (astat .ne. 0) then print *,'CANNOT ALLOCATE tempfld',astat stop endif allocate(start(rank),stride(rank),stat=astat) if (astat .ne. 0) then print *,'CANNOT ALLOCATE start or stride',astat stop endif start(:) = 0 stride(:) = 1 !======= read the field status = swrdfld(swid,fieldname,start,stride,dims(1:rank),tempfld) if (status .ne. 0) then print *,'CANNOT READ FIELD',status stop endif !======= close the file status = swclose(fid) if (status .ne. 0) then print *,'CANNOT CLOSE FILE',status stop endif var=tempfld(:) deallocate(tempfld,start,stride) end subroutine sdsreadeos1 end program