c ------------------ c make_cdf_ARCS_02.f c ------------------ program converting_to_netcdf c ---------------------------- c This program provides a model for converting a data set to netCDF. c The basic strategy used in this program is to open an existing netCDF c file, query the file for the ID's of the variables it contains, and c then write the data to those variables. c The output netCDF file must be created **before** this program is run. c The simplest way to do this is to cd to your scratch directory and c % cp $FER_DIR/doc/converting_to_netcdf.basic converting_to_netcdf.cdl c and then edit converting_to_netcdf.cdl (an ASCII file) to describe YOUR data c set. If your data set requires unequally spaced axes, climatological time c axes, staggered grids, etc. then converting_to_netcdf.supplement may be c a better starting point then the "basic" file used above. c After you edit converting_to_netcdf.cdl then create the netCDF file with c the command c % ncgen -o converting_to_netcdf.cdf converting_to_netcdf.cdl c Now we will read in **your** data (gridded oceanic temperature and c salt in this example) and write it out into the netCDF file c converting_to_netcdf.cdf. Note that the axis coordinates can be written c out exactly the same methodology - including time step values (as below). ************************************************************************* c An alternative to modifying this program is to use the command: c ncgen -f converting_to_netcdf.cdl c This will create a large source code to which select lines can c be added to write out your data. ************************************************************************* c To compile and link converting_to_netcdf.f, use: c f77 -o converting_to_netcdf converting_to_netcdf.f -lnetcdf ************************************************************************* c include file necessary for netCDF c update for your location include '/usr/local/netcdf/netcdf-3.5.0/include/netcdf.inc' c may be found in $FER_DIR/fmt/cmn ************************************************************************* c parameters relevant to the data being read in c THESE NORMALLY MATCH THE DIMENSIONS IN THE CDL FILE c (except nt which may be "unlimited") integer imt, jmt, km, nt, lnew, inlun parameter (imt=144, jmt=25, km=1, nt=240) c imt is longitude, jmt latitude, km depth, and nt number of time steps ************************************************************************* c variable declaration real datain(imt) real st(imt,jmt) real precip(imt,jmt) real slp(imt,jmt) real taux(imt,jmt) real tauy(imt,jmt) real lhflx(imt,jmt) real shflx(imt,jmt) real sswradflx(imt,jmt) real tswradflx(imt,jmt) real slwradflx(imt,jmt) real tlwradflx(imt,jmt) real tswclradflx(imt,jmt) real tlwclradflx(imt,jmt) real ws(imt,jmt) real rh(imt,jmt) real hold1(imt,jmt) real hold2(imt,jmt) real sum_st(imt,jmt,12) real sum_precip(imt,jmt,12) real sum_slp(imt,jmt,12) real sum_taux(imt,jmt,12) real sum_tauy(imt,jmt,12) real sum_lhflx(imt,jmt,12) real sum_shflx(imt,jmt,12) real sum_sswradflx(imt,jmt,12) real sum_tswradflx(imt,jmt,12) real sum_slwradflx(imt,jmt,12) real sum_tlwradflx(imt,jmt,12) real sum_tswclradflx(imt,jmt,12) real sum_tlwclradflx(imt,jmt,12) real sum_ws(imt,jmt,12) real sum_rh(imt,jmt,12) real tmp1,tmp2,tmp3 real tkf,elhe,elhs,shcv,shcw,shci,rdry,rvap,ppvsf,ts real elhv,ppvs,ppv,rhs real miss integer cdfid, rcode, index c netcdf field ids integer timeaxid integer stid integer precipid integer slpid integer tauxid integer tauyid integer lhflxid integer shflxid integer sswradflxid integer tswradflxid integer slwradflxid integer tlwradflxid integer tswclradflxid integer tlwclradflxid integer wsid integer rhid integer itime,i,j,k,l ************************************************************************* c dimension corner and step for defining size of gridded data integer corner1(4) integer step1(4) integer corner2(3) integer step2(3) c corner and step are used to define the size of the gridded data c to be written out. Since temp and salt are four dimensional arrays, c corner and step must be four dimensions as well. In each output c to my_data.cdf within the do loop, the entire array of data (160 long. c pts, 100 lat. pts., 27 depth pts.) will be written for one time step. c Corner tells netCDF where to start, and step indicates how many steps c in each dimension to take. data corner1/1, 1, 1, -1/ data corner2/1, 1, -1/ data step1/imt, jmt, km, 1/ data step2/imt, jmt, 1/ real dayofdat(12) data dayofdat / 14, 45, 74, 105, 135, 166, & 196, 227, 258, 288, 319, 349 / c values for specific humidity to relative humidity tkf = 273.16 ekhe = 2500000.0 elhs = 2834000.0 shcv = 1911.0 shcw = 4185.0 shci = 2060.0 rdry = 287.0 rvap = 461.5 ppvsf = 610.571 c ***NOTE*** Since Fortran and C access data differently, the order of c the variables in the Fortran code must be opposite that in the cdl c file. In Fortran, the first variable varies fastest while in C, the c last variables varies fastest. ************************************************************************** c initialize cdfid by using ncopn cdfid = ncopn('make_cdf_25x144_grid_clim_01.cdf', ncwrite, rcode) if (rcode.ne.ncnoerr) stop 'error with ncopn' ************************************************************************** miss = -9.99e12 c get variable id's by using ncvid c THE VARIABLE NAMES MUST MATCH THE CDL FILE (case sensitive) stid = ncvid(cdfid, 'st', rcode) if (rcode.ne.ncnoerr) stop 'error with stid' precipid = ncvid(cdfid, 'precip', rcode) if (rcode.ne.ncnoerr) stop 'error with precipid' slpid = ncvid(cdfid, 'slp', rcode) if (rcode.ne.ncnoerr) stop 'error with slpid' tauxid = ncvid(cdfid, 'taux', rcode) if (rcode.ne.ncnoerr) stop 'error with tauxid' tauyid = ncvid(cdfid, 'tauy', rcode) if (rcode.ne.ncnoerr) stop 'error with tauyid' lhflxid = ncvid(cdfid, 'lhflx', rcode) if (rcode.ne.ncnoerr) stop 'error with lhflxid' shflxid = ncvid(cdfid, 'shflx', rcode) if (rcode.ne.ncnoerr) stop 'error with shflxid' sswradflxid = ncvid(cdfid, 'sswradflx', rcode) if (rcode.ne.ncnoerr) stop 'error with sswradflxid' tswradflxid = ncvid(cdfid, 'tswradflx', rcode) if (rcode.ne.ncnoerr) stop 'error with tswradflxid' slwradflxid = ncvid(cdfid, 'slwradflx', rcode) if (rcode.ne.ncnoerr) stop 'error with slwradflxid' tlwradflxid = ncvid(cdfid, 'tlwradflx', rcode) if (rcode.ne.ncnoerr) stop 'error with tlwradflxid' tswclradflxid = ncvid(cdfid, 'tswclradflx', rcode) if (rcode.ne.ncnoerr) stop 'error with tswclradflxid' tlwclradflxid = ncvid(cdfid, 'tlwclradflx', rcode) if (rcode.ne.ncnoerr) stop 'error with tlwclradflxid' wsid = ncvid(cdfid, 'ws', rcode) if (rcode.ne.ncnoerr) stop 'error with wsid' rhid = ncvid(cdfid, 'rh', rcode) if (rcode.ne.ncnoerr) stop 'error with rhid' timeaxid = ncvid(cdfid, 'time', rcode) if (rcode.ne.ncnoerr) stop 'error with timeaxid' ************************************************************************** c this is a good place to open your data file open(unit=30, * file= * 'make_cdf_25x144_ts.dat', * status='old',form='unformatted') open(unit=31, * file= * 'make_cdf_25x144_precc.dat', * status='old',form='unformatted') open(unit=32, * file= * 'make_cdf_25x144_precc.dat', * status='old',form='unformatted') open(unit=33, * file= * 'make_cdf_25x144_psl.dat', * status='old',form='unformatted') open(unit=34, * file= * 'make_cdf_25x144_taux.dat', * status='old',form='unformatted') open(unit=35, * file= * 'make_cdf_25x144_tauy.dat', * status='old',form='unformatted') open(unit=36, * file= * 'make_cdf_25x144_lhflx.dat', * status='old',form='unformatted') open(unit=37, * file= * 'make_cdf_25x144_shflx.dat', * status='old',form='unformatted') open(unit=38, * file= * 'make_cdf_25x144_surswflx.dat', * status='old',form='unformatted') open(unit=39, * file= * 'make_cdf_25x144_toaswflx.dat', * status='old',form='unformatted') open(unit=40, * file= * 'make_cdf_25x144_surlwflx.dat', * status='old',form='unformatted') open(unit=41, * file= * 'make_cdf_25x144_toalwflx.dat', * status='old',form='unformatted') open(unit=42, * file= * 'make_cdf_25x144_surclskswflx.dat', * status='old',form='unformatted') open(unit=43, * file= * 'make_cdf_25x144_toaclskswflx.dat', * status='old',form='unformatted') open(unit=44, * file= * 'make_cdf_25x144_surclsklwflx.dat', * status='old',form='unformatted') open(unit=45, * file= * 'make_cdf_25x144_toaclsklwflx.dat', * status='old',form='unformatted') open(unit=46, * file= * 'make_cdf_25x144_u.dat', * status='old',form='unformatted') open(unit=47, * file= * 'make_cdf_25x144_v.dat', * status='old',form='unformatted') open(unit=48, * file= * 'make_cdf_25x144_q.dat', * status='old',form='unformatted') ************************************************************************** c begin do loop. Each step will read in one time step of data c and then write it out to my_data.cdf. c zero the climatology fields do l=1,12 do i=1,imt do j=1,jmt sum_st(i,j,l) = 0.0 sum_precip(i,j,l) = 0.0 sum_slp(i,j,l) = 0.0 sum_taux(i,j,l) = 0.0 sum_tauy(i,j,l) = 0.0 sum_lhflx(i,j,l) = 0.0 sum_shflx(i,j,l) = 0.0 sum_sswradflx(i,j,l) = 0.0 sum_tswradflx(i,j,l) = 0.0 sum_slwradflx(i,j,l) = 0.0 sum_tlwradflx(i,j,l) = 0.0 sum_tswclradflx(i,j,l) = 0.0 sum_tlwclradflx(i,j,l) = 0.0 sum_ws(i,j,l) = 0.0 sum_rh(i,j,l) = 0.0 end do end do end do index = 0 do itime = 1, nt write(6,*) itime index = index +1 if(index .eq. 13)then index = 1 end if c read data do j=1, jmt read(30) datain do i=1, imt if(datain(i) .lt. -1.0e30)then sum_st(i,j,index) = miss else st(i,j) = datain(i) - tkf sum_st(i,j,index) = sum_st(i,j,index) & + datain(i) end if end do end do do j=1, jmt read(31) datain do i=1, imt if(datain(i) .lt. -1.0e30)then sum_precip(i,j,index) = miss else sum_precip(i,j,index) = sum_precip(i,j,index) & + datain(i) end if end do end do do j=1, jmt read(32) datain do i=1, imt if(datain(i) .lt. -1.0e30)then sum_precip(i,j,index) = miss else sum_precip(i,j,index) = sum_precip(i,j,index) & + datain(i) end if end do end do do j=1, jmt read(33) datain do i=1, imt if(datain(i) .lt. -1.0e30)then sum_slp(i,j,index) = miss else slp(i,j) = (datain(i)/1000.0) - 1013.25 sum_slp(i,j,index) = sum_slp(i,j,index) & + datain(i) end if end do end do do j=1, jmt read(34) datain do i=1, imt if(datain(i) .lt. -1.0e30)then sum_taux(i,j,index) = miss else sum_taux(i,j,index) = sum_taux(i,j,index) & + datain(i) end if end do end do do j=1, jmt read(35) datain do i=1, imt if(datain(i) .lt. -1.0e30)then sum_tauy(i,j,index) = miss else sum_tauy(i,j,index) = sum_tauy(i,j,index) & + datain(i) end if end do end do do j=1, jmt read(36) datain do i=1, imt if(datain(i) .lt. -1.0e30)then sum_lhflx(i,j,index) = miss else sum_lhflx(i,j,index) = sum_lhflx(i,j,index) & + datain(i) end if end do end do do j=1, jmt read(37) datain do i=1, imt if(datain(i) .lt. -1.0e30)then sum_shflx(i,j,index) = miss else sum_shflx(i,j,index) = sum_shflx(i,j,index) & + datain(i) end if end do end do do j=1, jmt read(38) datain do i=1, imt if(datain(i) .lt. -1.0e30)then sum_sswradflx(i,j,index) = miss else sum_sswradflx(i,j,index) = sum_sswradflx(i,j,index) & + datain(i) end if end do end do do j=1, jmt read(39) datain do i=1, imt if(datain(i) .lt. -1.0e30)then sum_tswradflx(i,j,index) = miss else sum_tswradflx(i,j,index) = sum_tswradflx(i,j,index) & + datain(i) end if end do end do do j=1, jmt read(40) datain do i=1, imt if(datain(i) .lt. -1.0e30)then sum_slwradflx(i,j,index) = miss else sum_slwradflx(i,j,index) = sum_slwradflx(i,j,index) & + datain(i) end if end do end do do j=1, jmt read(41) datain do i=1, imt if(datain(i) .lt. -1.0e30)then sum_tlwradflx(i,j,index) = miss else sum_tlwradflx(i,j,index) = sum_tlwradflx(i,j,index) & + datain(i) end if end do end do do j=1, jmt read(42) datain do i=1, imt if(datain(i) .lt. -1.0e30)then hold1(i,j) = miss else hold1(i,j) = datain(i) end if end do end do do j=1, jmt read(43) datain do i=1, imt if(datain(i) .lt. -1.0e30)then hold2(i,j) = miss else hold2(i,j) = datain(i) end if end do end do do j=1, jmt do i=1, imt if(hold1(i,j) .lt. -1.0e30)then sum_tswclradflx(i,j,index) = miss else sum_tswclradflx(i,j,index) = & sum_tswclradflx(i,j,index) + & hold2(i,j) - hold1(i,j) end if end do end do do j=1, jmt read(44) datain do i=1, imt if(datain(i) .lt. -1.0e30)then hold1(i,j) = miss else hold1(i,j) = datain(i) end if end do end do do j=1, jmt read(45) datain do i=1, imt if(datain(i) .lt. -1.0e30)then hold2(i,j) = miss else hold2(i,j) = datain(i) end if end do end do do j=1, jmt do i=1, imt if(hold1(i,j) .lt. -1.0e30)then sum_tlwclradflx(i,j,index) = miss else sum_tlwclradflx(i,j,index) = & sum_tlwclradflx(i,j,index) + & hold2(i,j) - hold1(i,j) end if end do end do do j=1, jmt read(46) datain do i=1, imt if(datain(i) .lt. -1.0e30)then hold1(i,j) = miss else hold1(i,j) = datain(i) end if end do end do do j=1, jmt read(47) datain do i=1, imt if(datain(i) .lt. -1.0e30)then hold2(i,j) = miss else hold2(i,j) = datain(i) end if end do end do do j=1, jmt do i=1, imt if(hold1(i,j) .lt. -1.0e30)then sum_ws(i,j,index) = miss else sum_ws(i,j,index) = & sum_ws(i,j,index) + & sqrt(hold2(i,j)**2 + hold1(i,j)**2) end if end do end do do j=1, jmt read(48) datain do i=1, imt if(datain(i) .lt. -1.0e30)then rh(i,j) = 0.0 else c convert Kg/Kg to g/Kg hold2(i,j) = datain(i)*1000.00 end if end do end do c compute relative humidity from specific humidity do j=1, jmt do i=1, imt if(st(i,j) .gt. 0.0)then elhv = (shcv-shcw)*st(i,j) else elhs = (shcv-shci)*st(i,j) end if ppvs = ppvsf * exp( & elhv*(1.0/tkf - 1.0/(st(i,j)+tkf))/rvap & ) ppv = rvap*(hold2(i,j)*0.001)* & (slp(i,j)*100.0+101325.0)/ & (rdry+(rvap-rdry)*(hold2(i,j)*0.001)) rhs = 100.0*ppv/ppvs sum_rh(i,j,index) = sum_rh(i,j,index) + rhs end do end do end do c compute number of years tmp1 = float(nt)/12.0 write(6,*) 'the number of years is ', tmp1 do l=1,12 do i=1,imt do j=1,jmt if(sum_st(i,j,l) .ne. miss)then sum_st(i,j,l) = sum_st(i,j,l)/tmp1 end if if(sum_precip(i,j,l) .ne. miss)then sum_precip(i,j,l) = sum_precip(i,j,l)/tmp1 end if if(sum_slp(i,j,l) .ne. miss)then sum_slp(i,j,l) = sum_slp(i,j,l)/tmp1 end if if(sum_taux(i,j,l) .ne. miss)then sum_taux(i,j,l) = sum_taux(i,j,l)/tmp1 end if if(sum_tauy(i,j,l) .ne. miss)then sum_tauy(i,j,l) = sum_tauy(i,j,l)/tmp1 end if if(sum_lhflx(i,j,l) .ne. miss)then sum_lhflx(i,j,l) = sum_lhflx(i,j,l)/tmp1 end if if(sum_shflx(i,j,l) .ne. miss)then sum_shflx(i,j,l) = sum_shflx(i,j,l)/tmp1 end if if(sum_sswradflx(i,j,l) .ne. miss)then sum_sswradflx(i,j,l) = sum_sswradflx(i,j,l)/tmp1 end if if(sum_tswradflx(i,j,l) .ne. miss)then sum_tswradflx(i,j,l) = sum_tswradflx(i,j,l)/tmp1 end if if(sum_slwradflx(i,j,l) .ne. miss)then sum_slwradflx(i,j,l) = sum_slwradflx(i,j,l)/tmp1 end if if(sum_tlwradflx(i,j,l) .ne. miss)then sum_tlwradflx(i,j,l) = sum_tlwradflx(i,j,l)/tmp1 end if if(sum_tswclradflx(i,j,l) .ne. miss)then sum_tswclradflx(i,j,l) = sum_tswclradflx(i,j,l)/tmp1 end if if(sum_tlwclradflx(i,j,l) .ne. miss)then sum_tlwclradflx(i,j,l) = sum_tlwclradflx(i,j,l)/tmp1 end if if(sum_ws(i,j,l) .ne. miss)then sum_ws(i,j,l) = sum_ws(i,j,l)/tmp1 end if if(sum_rh(i,j,l) .ne. miss)then sum_rh(i,j,l) = sum_rh(i,j,l)/tmp1 end if end do end do end do close(30) close(31) close(32) ************************************************************************** c begin do loop. Each step will read in one time step of data c and then write it out to my_data.cdf. write(6,*) 'Data writing routine' index = 0 do itime = 1, 12 time_step = dayofdat(itime) index = index +1 if(index .eq. 13)then index = 1 end if c initialize time step in corner corner2(3) = itime corner1(4) = itime do i=1,imt do j=1,jmt c convert from K to C st(i,j) = sum_st(i,j,itime) - 273.16 c convert form m/sec to mm/day precip(i,j) = sum_precip(i,j,itime)*8.64e7 slp(i,j) = sum_slp(i,j,itime) taux(i,j) = sum_taux(i,j,itime) tauy(i,j) = sum_tauy(i,j,itime) lhflx(i,j) = sum_lhflx(i,j,itime) shflx(i,j) = sum_shflx(i,j,itime) sswradflx(i,j) = sum_sswradflx(i,j,itime) tswradflx(i,j) = sum_tswradflx(i,j,itime) slwradflx(i,j) = sum_slwradflx(i,j,itime) tlwradflx(i,j) = sum_tlwradflx(i,j,itime) tswclradflx(i,j) = sum_tswclradflx(i,j,itime) tlwclradflx(i,j) = sum_tlwclradflx(i,j,itime) ws(i,j) = sum_ws(i,j,itime) rh(i,j) = sum_rh(i,j,itime) end do end do c writing the data to the netCDF file write (6,*) 'writing time step: ',itime, time_step write (6,*) 'writing st' call ncvpt(cdfid,stid,corner2,step2,st(1,1),rcode) if (rcode.ne.ncnoerr) stop 'error with st-put' write (6,*) 'writing precip' call ncvpt(cdfid,precipid,corner2,step2,precip(1,1),rcode) if (rcode.ne.ncnoerr) stop 'error with precip-put' write (6,*) 'writing slp' call ncvpt(cdfid,slpid,corner2,step2,slp(1,1),rcode) if (rcode.ne.ncnoerr) stop 'error with slp-put' write (6,*) 'writing taux' call ncvpt(cdfid,tauxid,corner2,step2,taux(1,1),rcode) if (rcode.ne.ncnoerr) stop 'error with taux-put' write (6,*) 'writing tauy' call ncvpt(cdfid,tauyid,corner2,step2,tauy(1,1),rcode) if (rcode.ne.ncnoerr) stop 'error with tauy-put' write (6,*) 'writing lhflx' call ncvpt(cdfid,lhflxid,corner2,step2,lhflx(1,1),rcode) if (rcode.ne.ncnoerr) stop 'error with lhflx-put' write (6,*) 'writing shflx' call ncvpt(cdfid,shflxid,corner2,step2,shflx(1,1),rcode) if (rcode.ne.ncnoerr) stop 'error with shflx-put' write (6,*) 'writing sswradflx' call ncvpt(cdfid,sswradflxid,corner2,step2,sswradflx(1,1),rcode) if (rcode.ne.ncnoerr) stop 'error with sswradflx-put' write (6,*) 'writing tswradflx' call ncvpt(cdfid,tswradflxid,corner2,step2,tswradflx(1,1),rcode) if (rcode.ne.ncnoerr) stop 'error with tswradflx-put' write (6,*) 'writing slwradflx' call ncvpt(cdfid,slwradflxid,corner2,step2,slwradflx(1,1),rcode) if (rcode.ne.ncnoerr) stop 'error with slwradflx-put' write (6,*) 'writing tlwradflx' call ncvpt(cdfid,tlwradflxid,corner2,step2,tlwradflx(1,1),rcode) if (rcode.ne.ncnoerr) stop 'error with tlwradflx-put' write (6,*) 'writing tswclradflx' call ncvpt(cdfid,tswclradflxid,corner2,step2, & tswclradflx(1,1),rcode) if (rcode.ne.ncnoerr) stop 'error with tswclradflx-put' write (6,*) 'writing tlwclradflx' call ncvpt(cdfid,tlwclradflxid,corner2,step2, & tlwclradflx(1,1),rcode) if (rcode.ne.ncnoerr) stop 'error with tlwclradflx-put' write (6,*) 'writing ws' call ncvpt(cdfid,wsid,corner2,step2,ws(1,1),rcode) if (rcode.ne.ncnoerr) stop 'error with ws-put' write (6,*) 'writing rh' call ncvpt(cdfid,rhid,corner2,step2,rh(1,1),rcode) if (rcode.ne.ncnoerr) stop 'error with rh-put' call ncvpt1(cdfid,timeaxid,itime,time_step,rcode) if (rcode.ne.ncnoerr) stop 'error with timax-put' end do close(32) close(33) close(34) close(36) close(37) close(38) close(39) close(40) close(41) close(42) close(43) close(44) close(46) close(47) close(48) ************************************************************************** c close my_data.cdf using ncclos call ncclos(cdfid, rcode) if (rcode.ne.ncnoerr) stop 'error with ncclos' ************************************************************************** stop end