program cloudsat use netcdf implicit none integer,parameter :: altmax2=125 character(len=200) :: filein, file2, modfileout character :: sds_varname*100, filetmp2*1024, filetmp3*1024, date*8 integer :: err, imonth, id, ifile,numfich, nc_id integer :: ii, it, ilon, ilat, ialt,ilid, nprof, iz, i, itemp integer,parameter :: altmax=22 integer,parameter :: lonmax=144,latmax=72 real,dimension(altmax) :: altmid real,dimension(altmax+1) :: altmod real,dimension(latmax+1) :: latmod real,dimension(lonmax+1) :: lonmod real,dimension(latmax) :: latmid real,dimension(lonmax) :: lonmid real,dimension(:),allocatable :: lat,lon,vartmp1 ! instant lat/lon real,dimension(2,altmax,latmax,lonmax) :: QR3d, QR3dind ! Cloud fraction !real,dimension(2,latmax,lonmax) :: RH3d, RH3dind ! Cloud fraction real,dimension(:,:),allocatable :: vartmp2, height !,RH real,dimension(:,:,:),allocatable :: QR,vartmp3 100 format(A200) ! Initialization QR3d(:,:,:,:)=0. !RH3d(:,:,:)=0. QR3dind(:,:,:,:)=0. !RH3dind(:,:,:)=0. altmod = (/1012.5, 987.5, 962.5, 937.5, 912.5, 875., 825., 775., 725., 675., 625., 575., 525., 475., 425., 375., 325., 275., 225., 175., 125., 75., 25./) altmid = (/1000., 975., 950., 925., 900., 850., 800., 750., 700., 650., 600., 550., 500., 450., 400., 350., 300., 250., 200., 150., 100., 50./) latmod(1)=-90 do ilat=1,latmax latmod(ilat+1)=latmod(ilat)+2.5 enddo latmid=latmod(1:latmax)+1 lonmod(1)=-180 do ilon=1,lonmax lonmod(ilon+1)=lonmod(ilon)+2.5 enddo lonmid=lonmod(1:lonmax)+1 !altmod(1)=0 !do ilat=1,altmax !altmod(ilat+1)=altmod(ilat)+0.48 !enddo !altmid=altmod(1:altmax)+0.24 !------------------- 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 ! Initialization of box and numfich numfich=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) !filetmp2='/bdd/CFMIP/OBS_LOCAL/CloudSat/FLXHR-LIDAR/ddddd/2007012003305_03767_CS_2B-FLXHR-LIDAR_GRANULE_P2_R04_E02.hdf' !filetmp2='/data1/jli/CLOUDSAT_data/CLOUDSAT_hdf_link/2009/2009001021530_14253_CS_2B-FLXHR-LIDAR_GRANULE_P2_R04_E02.hdf' ! ECMWF-AUX files for pressure !print *, trim(filetmp2(1:71)) !print *, ' espace' !print *, trim(filetmp2(86:95)) !print *, ' espace' !print *, trim(filetmp2(97:)) !print *, ' espace' filetmp3=filetmp2(1:73)//'ECMWF-AUX'//filetmp2(88:97)//filetmp2(99:) print *, trim(filetmp3) print *, '' print *, 'Read the CloudSat file data ...' !-------------------- Get the number of profiles ----------------! call sdsreadprof(filetmp2,it) !----------------------- Allocation -----------------------------! allocate(lat(it),lon(it),height(altmax2,it),QR(altmax2,it,2)) !,RH(it,2)) lat(:)=0. lon(:)=0. height(:,:)=0. !RH(:,:)=0. QR(:,:,:)=0. !------------------- Get thevariables of interest ---------------! 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) sds_varname='Pressure' print *, trim(sds_varname) call sdsreadeos2ECMWF(vartmp2,filetmp3,trim(sds_varname)) !call sdsreadeos2(vartmp2,filetmp2,trim(sds_varname)) height=vartmp2/100 deallocate(vartmp2) !sds_varname='RH' !print *, trim(sds_varname) !call sdsreadeos2(vartmp2,filetmp2,trim(sds_varname)) !RH=vartmp2 !where(RH.eq.-9.98) !RH=-9999. !endwhere !deallocate(vartmp2) sds_varname='QR' print *, trim(sds_varname) call sdsreadeos3(vartmp3,filetmp2,trim(sds_varname)) QR=vartmp3 where(QR.eq.-9.98) QR=-9999. endwhere deallocate(vartmp3) !allocate(vartmp2(altmax,it)) ! Flag night !vartmp2(:,:)=0. !vartmp2=QR(:,:,1); !where(vartmp2.eq.-9999) ! vartmp2=0. !endwhere !------------------------- AVERAGE ON 3D GRID ----------------------! do i=1,it ! profile ! if(sum(vartmp2(:,i),1).lt.0.01)then ! QR(:,i,1)=-9999. ! endif 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=altmax,1,-1 do ilid=1,altmax2 if ((height(ilid,i).le.altmod(iz)).and.(height(ilid,i).gt.altmod(iz+1)))then do ii=1,2 if(QR(ilid,i,ii).ne.-9999.)then QR3d(ii,iz,ilat,ilon)=QR3d(ii,iz,ilat,ilon)+QR(ilid,i,ii) QR3dind(ii,iz,ilat,ilon)=QR3dind(ii,iz,ilat,ilon)+1 endif !if(RH(i,ii).ne.-9999.)then ! RH3d(ii,ilat,ilon)=RH3d(ii,ilat,ilon)+RH(i,ii) ! RH3dind(ii,ilat,ilon)=RH3dind(ii,ilat,ilon)+1 !endif enddo endif enddo enddo endif enddo endif enddo enddo !------------------------- Deallocate ----------------------------! deallocate(lat,lon,height,QR)!,RH) !vartmp2) print *, 'next file' goto 888 ! next CloudSat file 999 continue !----------------------- Exclude nan values -----------------------! do ilon=1,lonmax do ilat=1,latmax do ii=1,2 do iz = 1, altmax ! loop by vertical bins if(QR3dind(ii,iz,ilat,ilon)>0.)then QR3d(ii,iz,ilat,ilon)=QR3d(ii,iz,ilat,ilon)/QR3dind(ii,iz,ilat,ilon) else QR3d(ii,iz,ilat,ilon)=-9999. endif enddo ! if(RH3dind(ii,ilat,ilon)>0.)then ! RH3d(ii,ilat,ilon)=RH3d(ii,ilat,ilon)/RH3dind(ii,ilat,ilon) ! else ! RH3d(ii,ilat,ilon)=-9999. ! endif enddo enddo enddo !---------------------- RECORD OUTPUT FILES ---------------------! print *, 'Record Output Files' print *, 'file1' modfileout='/home/gcesana/CloudSat/FLXHR-LIDAR_gridded/CloudSat_QR_'//trim(date)//'_ECMWF.nc' !modfileout='/home/gcesana/CloudSat/CloudSat_QR_'//trim(date)//'_ECMWF.nc' !modfileout='/bdd/CFMIP/OBS_LOCAL/CloudSat/CloudSat_QR_'//trim(date)//'_ECMWF.nc' print *, modfileout call create_3d(trim(modfileout),'QR','lon','lat','alt', & lonmax,latmax,altmax,QR3d,lonmid,latmid,altmid) !print *, 'file2' !modfileout='/home/gcesana/CloudSat/2B-FLXHR-LIDAR_gridded/CloudSat_RH_'//trim(date)//'.nc' !modfileout='/bdd/CFMIP/OBS_LOCAL/CloudSat/FLXHR-LIDAR_gridded/CloudSat_RH_'//trim(date)//'.nc' !call create_2d(trim(modfileout),'RH','lon','lat', & ! lonmax,latmax,RH3d,lonmid,latmid) !----------------------- END OF PROG ---------------------------! print *, 'End of prog' 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 !----------------------------------------------------------------------------! ! **** SDSREAD **** This routine allows to read a 2dim SDS-variable on ! hdf ! ! files ! ! !----------------------------------------------------------------------------! subroutine sdsreadeos(var,filename,varname) implicit none ! use hdf include "hdf.f90" include "dffunc.f90" real,dimension(:,:,:),allocatable :: var character :: varname*100 integer (kind=2), allocatable :: tempfld(:,:,:) character (len=128) :: filename,swathname,fieldname,attribute1,attribute2 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(:) ! print *,'Reading HDF file...' ! swathname = "2B-GEOPROF" ! fieldname = "Radar_Reflectivity" swathname = "2B-FLXHR-LIDAR" fieldname = "QR" fid = swopen(filename, DFACC_READ) if (fid .eq. -1) then print *,'CANNOT OPEN FILE',fid stop endif print *,fid !======= 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),dims(2),dims(3)),var(dims(1),dims(2),dims(3)),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 !deallocate(tempfld,start,stride) !print *, 'vartmp',tempfld(60,5000,1) end subroutine sdsreadeos !----------------------------------------------------------------------------! 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 !======= 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 var=tempfld(:) deallocate(tempfld,start,stride) end subroutine sdsreadeos1 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 sdsreadeos2(var,filename,fieldname) implicit none ! use hdf include "hdf.f90" include "dffunc.f90" real,allocatable :: var(:,:) character(len=*) :: fieldname integer (kind=2), allocatable :: tempfld(:,:) character (len=128) ::filename,swathname,attribute1,attribute2 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(:) real (kind=4) :: factor,offset ! field attributes (used for scaling) swathname = "2B-FLXHR-LIDAR" attribute1 = fieldname//".factor" attribute2 = fieldname//".offset" 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),dims(2)),var(dims(1),dims(2)),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 !======= read the factor attribute status = swrdattr(swid,attribute1,factor) !======= read the offset attribute status = swrdattr(swid,attribute2,offset) !======= scale the Radar_reflectivity field var = (tempfld - offset)/ factor !======= 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 deallocate(tempfld,start,stride) end subroutine sdsreadeos2 subroutine sdsreadeos2ECMWF(tempfld,filename,fieldname) implicit none ! use hdf include "hdf.f90" include "dffunc.f90" 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 = "ECMWF-AUX" 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),dims(2)),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 !======= 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 deallocate(start,stride) end subroutine sdsreadeos2ECMWF subroutine sdsreadeos3(var,filename,fieldname) implicit none ! use hdf include "hdf.f90" include "dffunc.f90" real,allocatable :: var(:,:,:) character(len=*) :: fieldname integer (kind=2), allocatable :: tempfld(:,:,:) character (len=128) ::filename,swathname,attribute1,attribute2 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(:) real (kind=4) :: factor,offset ! field attributes (used for scaling) swathname = "2B-FLXHR-LIDAR" attribute1 = fieldname//".factor" attribute2 = fieldname//".offset" 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 !print *,'-------------------------------------------------' !print *,'file: ',trim(filename) !print *,'field: ',trim(fieldname) !print *,'rank: ',rank !print *,'dimension list: ',trim(dimlist) !print *,'dimension sizes: ',dims(1:rank) ! print *,'datatype: ',gettype(datatype) !======= allocate the field to be read in allocate(tempfld(dims(1),dims(2),dims(3)),var(dims(1),dims(2),dims(3)),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 !======= read the factor attribute status = swrdattr(swid,attribute1,factor) !print *,'-------------------------------------------------' !print *,'attribute1: ',trim(attribute1) !print *,'value: ',factor !======= read the offset attribute status = swrdattr(swid,attribute2,offset) !print *,'-------------------------------------------------' !print *,'attribute2: ',trim(attribute2) !print *,'value: ',offset !======= scale the Radar_reflectivity field var = (tempfld - offset)/ factor !======= 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 deallocate(tempfld,start,stride) end subroutine sdsreadeos3 subroutine sdsread3d(var,filename,varname) implicit none include "hdf.f90" include "dffunc.f90" integer :: i, ret, start(2) = 0, stride(2) = 1 ,edges(3) integer :: sd_id, sds_id integer :: ivar, npts, nprofs, nb, OK_buffer real,dimension(:,:,:),allocatable :: var integer*4 :: index, rank, istat, attributes, num_type integer*4 :: dim_sizes (32) character(len=232) :: name, filename character :: varname*100 ! print *,'Reading HDF file...' sd_id = sfstart (filename, DFACC_READ) do i=2,100 nb=i ! Select and read the var sds_id = sfselect (sd_id, nb) istat = sfginfo (sds_id, name, rank, dim_sizes, num_type,attributes) if(name==trim(varname))exit enddo ! print *, name,'dimensions : ', dim_sizes(1:2) ivar = dim_sizes(1) npts = dim_sizes(2) nprofs = dim_sizes(3) allocate(var(ivar, npts, nprofs),stat = OK_buffer) if (OK_buffer/=0) print *,'--- buffer allocation of',trim(name),'error' edges = [ivar, npts, nprofs] ret = sfrdata (sds_id, start, stride, edges, var) ! print *,'Reading : ', name ! Do something with the data if (ret.eq.-1) then print *,'ERROR' end if ! Close the crap ret = sfendacc(sds_id) ret = sfend(sd_id) end subroutine sdsread3d ! 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=2 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,dimension(dim4,dim3,dim2,dim1) :: var5 real,dimension(dim1) :: lonmod real,dimension(dim2) :: latmod real,dimension(dim3) :: tempmod ! CREATION DU FICHIER nf90_create call check(nf90_create(trim(fname), NF90_CLOBBER, & nc_id)) call check(nf90_put_att(nc_id, NF90_GLOBAL, 'Description',& 'CloudSat Heating Rate gridded files')) call check(nf90_put_att(nc_id, NF90_GLOBAL,'Version','~/src/CloudSat.f90')) call check(nf90_put_att(nc_id, NF90_GLOBAL, 'Author',& 'Gregory CESANA; Cesana et al., 2015, Evaluation of Cloud and Heating rate & Profiles in Eight GCMs')) print *, 'attributes' ! 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)//'sw', NF90_FLOAT, dim,varid5)) call check(nf90_put_att(nc_id, varid5, '_FillValue',nan)) call check(nf90_def_var(nc_id, trim(vname)//'lw', NF90_FLOAT, dim,varid6)) call check(nf90_put_att(nc_id, varid6, '_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_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 integer,parameter :: dim4=2 integer :: i,dim(ndims),nc_id,dim1,dim2 integer :: varid1, varid2,varid5,varid6,varid7 integer :: dimid1,dimid2,dimid4, time_dimid real*4,dimension(dim4,dim2,dim1) :: var5 real,dimension(dim1) :: lonmod real,dimension(dim2) :: latmod ! 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 3CALIPSO & 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, 'time', NF90_UNLIMITED, time_dimid)) ! DIMENSION DES MATRICES dim = (/dimid2, dimid1, time_dimid/) ! VARIABLE 4D +time call check(nf90_def_var(nc_id, trim(vname)//'sw', NF90_FLOAT, dim,varid5)) call check(nf90_put_att(nc_id, varid5, '_FillValue',nan)) call check(nf90_def_var(nc_id, trim(vname)//'lw', NF90_FLOAT, dim,varid6)) call check(nf90_put_att(nc_id, varid6, '_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)) ! 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, varid5, var5(1,:,:))) call check(nf90_put_var(nc_id, varid6, var5(2,:,:))) call check(nf90_close(nc_id)) end subroutine create_2d end program