program cloudsat use netcdf implicit none integer,parameter :: altmax2=125 character(len=200) :: filein, file2, modfileout character :: sds_varname*100, filetmp2*1024, date*8 integer :: err, imonth, id, ifile,numfich integer :: ii, it, ilon, ilat, ialt,ilid, nprof, iz, i, itemp integer,parameter :: altmax=40 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 :: RH,height,vartmp2 real,dimension(:,:,:),allocatable :: QR,vartmp3 100 format(A200) ! Initialization QR3d(:,:,:,:)=0. RH3d(:,:,:)=0. QR3dind(:,:,:,:)=0. RH3d(:,:,:)=0. 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) print *, 'Read the CloudSat file data ...' !-------------------- Get the number of profiles ----------------! call sdsreadprof(filetmp2,it) !----------------------- Allocation -----------------------------! allocate(lat(it),lon(it),height(altmax2,it),RH(it,2),QR(altmax2,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='Height' print *, trim(sds_varname) call sdsreadeos2(vartmp2,filetmp2,trim(sds_varname)) height=vartmp2 height=height/1000. deallocate(vartmp2) sds_varname='RH' print *, trim(sds_varname) call sdsreadeos2(vartmp2,filetmp2,trim(sds_varname)) RH=vartmp2 where(RH.eq.-9.99) 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.99) QR=-9999. endwhere deallocate(vartmp3) !sds_varname='QR' !print *, trim(sds_varname),'varname' !call sdsreadeos3(QR,filetmp2,trim(sds_varname)) !------------------------- AVERAGE ON 3D GRID ----------------------! do i=1,it ! profile 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).ge.altmod(iz)).and.(height(ilid,i).lt.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) 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='/bdd/CFMIP/OBS_LOCAL/CloudSat/FLXHR_gridded/CloudSat_QR_'//trim(date)//'.nc' call create_3d(trim(modfileout),'QR','lon','lat','alt', & lonmax,latmax,altmax,QR3d,lonmid,latmid,altmid) print *, 'file2' modfileout='/bdd/CFMIP/OBS_LOCAL/CloudSat/FLXHR_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" 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" 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" 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" 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 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" 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*4,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(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')) ! 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