Blame | Last modification | View Log | Download | RSS feed
PROGRAM p2zc Calculate the geopotential and add it to the P filec Michael Sprenger / Spring 2006implicit nonec ---------------------------------------------------------------c Declaration of variablesc ---------------------------------------------------------------c Variables for input P file : model levelcharacter*80 ml_pfnreal ml_varmin(4),ml_varmax(4),ml_stag(4)integer ml_vardim(4)real ml_mdvinteger ml_ndiminteger ml_nx,ml_ny,ml_nzreal ml_xmin,ml_xmax,ml_ymin,ml_ymax,ml_dx,ml_dyinteger ml_ntimesreal ml_aklev(500),ml_bklev(500)real ml_aklay(500),ml_bklay(500)real ml_timereal ml_pollon,ml_pollatinteger ml_nvarscharacter*80 ml_vnam(100)integer ml_idate(5)real,allocatable, dimension (:,:) :: ml_ps,ml_zbreal,allocatable, dimension (:,:,:) :: ml_t3,ml_q3,ml_p3,ml_tv3real,allocatable, dimension (:,:,:) :: ml_z3c Variables for input Z file : pressure levelcharacter*80 pl_zfnreal pl_varmin(4),pl_varmax(4),pl_stag(4)integer pl_vardim(4)real pl_mdvinteger pl_ndiminteger pl_nx,pl_ny,pl_nzreal pl_xmin,pl_xmax,pl_ymin,pl_ymax,pl_dx,pl_dyinteger pl_ntimesreal pl_aklev(500),pl_bklev(500)real pl_aklay(500),pl_bklay(500)real pl_timereal pl_pollon,pl_pollatinteger pl_nvarscharacter*80 pl_vnam(100)integer pl_idate(5)real,allocatable, dimension (:,:,:) :: pl_z3,pl_p3c Variables for output Z file : height levelcharacter*80 zl_ofnreal zl_varmin(4),zl_varmax(4),zl_stag(4)integer zl_vardim(4)real zl_mdvinteger zl_ndiminteger zl_nx,zl_ny,zl_nzreal zl_xmin,zl_xmax,zl_ymin,zl_ymax,zl_dx,zl_dyinteger zl_ntimesreal zl_aklev(500),zl_bklev(500)real zl_aklay(500),zl_bklay(500)real zl_timereal zl_pollon,zl_pollatinteger zl_idate(5)real,allocatable, dimension (:,:,:) :: zl_fieldreal,allocatable, dimension (:,:,:) :: zl_pc Variables for input P,S file : model levelcharacter*80 in_ofnreal in_varmin(4),in_varmax(4),in_stag(4)integer in_vardim(4)real in_mdvinteger in_ndiminteger in_nx,in_ny,in_nzreal in_xmin,in_xmax,in_ymin,in_ymax,in_dx,in_dyinteger in_ntimesreal in_aklev(500),in_bklev(500)real in_aklay(500),in_bklay(500)real in_timereal in_pollon,in_pollatinteger in_nvarscharacter*80 in_vnam(100)integer in_idate(5)real,allocatable, dimension (:,:,:) :: in_fieldc Physical and numerical parametersreal gparameter (g=9.80665)real epsparameter (eps=0.01)real tzeroparameter (tzero=273.15)real kappaparameter (kappa=0.6078)real zerodivparameter (zerodiv=0.0000000001)real dpminparameter (dpmin=10.)real rdgparameter (rdg=29.271)c Flag for test modeinteger testparameter (test=0)character*80 testfileparameter (testfile='TEST')c Variables and levels for interpolation onto z levelscharacter*80 levelfile,varfileinteger nvar,nlevcharacter*80 varinp(100),varout(100),varsrc(100)real zlev(500)c Auxiliary variablesinteger ierrinteger cdfid,cstidcharacter*80 cfninteger statreal timereal tv1(1000),z1(1000),p1(1000),f1(1000)real spline_tv1(1000),spline_f1(1000),spline_z1(1000)real pu,po,zu,zo,p,z,dp,p0,tvu,tvo,ffinteger i,j,k,linteger lmin,ncharacter*80 varname,cdfnameinteger idate(5),stdate(5),datar(14)integer isokcharacter*80 namereal zmin,dzc -----------------------------------------------------------------c Read input fieldsc -----------------------------------------------------------------print*,'*********************************************************'print*,'* p2z *'print*,'*********************************************************'c Read in the parameter fileopen(10,file='fort.10')read(10,*) ml_pfnread(10,*) pl_zfnread(10,*) zl_ofnread(10,*) name,zminif ( trim(name).ne.'GEO_ZMIN') stopread(10,*) name,nlevif ( trim(name).ne.'GEO_NZ' ) stopread(10,*) name,dzif ( trim(name).ne.'GEO_DZ' ) stopdo i=1,nlevzlev(i)=zmin+real(i-1)*dzenddonvar=1102 continueread(10,*,end=103) varinp(nvar),varout(nvar),varsrc(nvar)nvar=nvar+1goto 102103 continuenvar=nvar-1print*do i=1,nvarwrite(*,'(a10,a10,a30)') trim(varinp(i)),trim(varout(i)),> trim(varsrc(i))enddoprint*close(10)c Get grid description for P file : model levelcall cdfopn(ml_pfn,cdfid,ierr)if (ierr.ne.0) goto 998call getcfn(cdfid,cfn,ierr)if (ierr.ne.0) goto 998call cdfopn(cfn,cstid,ierr)if (ierr.ne.0) goto 998call getvars(cdfid,ml_nvars,ml_vnam,ierr)varname='T'isok=0call check_varok (isok,varname,ml_vnam,ml_nvars)if (isok.ne.1) goto 998call getdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim,> ml_varmin,ml_varmax,ml_stag,ierr)if (ierr.ne.0) goto 998ml_nx =ml_vardim(1)ml_ny =ml_vardim(2)ml_nz =ml_vardim(3)ml_xmin=ml_varmin(1)ml_ymin=ml_varmin(2)call getlevs(cstid,ml_nz,ml_aklev,ml_bklev,ml_aklay,ml_bklay,ierr)call getgrid(cstid,ml_dx,ml_dy,ierr)ml_xmax=ml_xmin+real(ml_nx-1)*ml_dxml_ymax=ml_ymin+real(ml_ny-1)*ml_dycall gettimes(cdfid,ml_time,ml_ntimes,ierr)call getstart(cstid,ml_idate,ierr)call getpole(cstid,ml_pollon,ml_pollat,ierr)call clscdf(cstid,ierr)call clscdf(cdfid,ierr)c Get grid description for Z file : pressure levelcall cdfopn(pl_zfn,cdfid,ierr)if (ierr.ne.0) goto 998call getcfn(cdfid,cfn,ierr)if (ierr.ne.0) goto 998call cdfopn(cfn,cstid,ierr)if (ierr.ne.0) goto 998call getvars(cdfid,pl_nvars,pl_vnam,ierr)varname='Z'isok=0call check_varok (isok,varname,pl_vnam,pl_nvars)if (isok.ne.1) goto 998call getdef(cdfid,varname,pl_ndim,pl_mdv,pl_vardim,> pl_varmin,pl_varmax,pl_stag,ierr)if (ierr.ne.0) goto 998pl_nx =pl_vardim(1)pl_ny =pl_vardim(2)pl_nz =pl_vardim(3)pl_xmin=pl_varmin(1)pl_ymin=pl_varmin(2)call getlevs(cstid,pl_nz,pl_aklev,pl_bklev,pl_aklay,pl_bklay,ierr)call getgrid(cstid,pl_dx,pl_dy,ierr)pl_xmax=pl_xmin+real(pl_nx-1)*pl_dxpl_ymax=pl_ymin+real(pl_ny-1)*pl_dycall gettimes(cdfid,pl_time,pl_ntimes,ierr)call getstart(cstid,pl_idate,ierr)call getpole(cstid,pl_pollon,pl_pollat,ierr)call clscdf(cstid,ierr)call clscdf(cdfid,ierr)c Set grid description for output file : height levelzl_vardim(1) = pl_vardim(1)zl_vardim(2) = pl_vardim(2)zl_vardim(3) = nlevzl_varmin(1) = pl_varmin(1)zl_varmin(2) = pl_varmin(2)zl_varmin(3) = zlev(1)zl_varmax(1) = pl_varmax(1)zl_varmax(2) = pl_varmax(2)zl_varmax(3) = zlev(nlev)do i=1,nlevzl_aklay(i) = zlev(i)zl_bklay(i) = 0.zl_aklev(i) = zlev(i)zl_bklev(i) = 0.enddozl_dx = pl_dxzl_dy = pl_dyzl_time = pl_timezl_ntimes = pl_ntimeszl_ndim = 4zl_mdv = pl_mdvzl_nx = zl_vardim(1)zl_ny = zl_vardim(2)zl_nz = zl_vardim(3)zl_xmin = zl_varmin(1)zl_ymin = zl_varmin(2)zl_pollon = ml_pollonzl_pollat = ml_pollatdo i=1,5zl_idate(i) = ml_idate(i)enddoc Consitency check for the gridsif ( (ml_nx.ne.pl_nx).or.> (ml_ny.ne.pl_ny).or.> (abs(ml_xmin-pl_xmin ).gt.eps).or.> (abs(ml_ymin-pl_ymin ).gt.eps).or.> (abs(ml_xmax-pl_xmax ).gt.eps).or.> (abs(ml_ymax-pl_ymax ).gt.eps).or.> (abs(ml_dx -pl_dx ).gt.eps).or.> (abs(ml_dy -pl_dy ).gt.eps).or.> (abs(ml_time-pl_time ).gt.eps).or.> (abs(ml_pollon-pl_pollon).gt.eps).or.> (abs(ml_pollat-pl_pollat).gt.eps)) thenprint*,'Input P and Z grids are not consistent... Stop'print*print*,'Xmin : ',ml_xmin,pl_xminprint*,'Ymin : ',ml_ymin,pl_yminprint*,'Xmax : ',ml_xmax,pl_xmaxprint*,'Ymax : ',ml_ymax,pl_ymaxprint*,'Dx : ',ml_dx, pl_dxprint*,'Dy : ',ml_dy, pl_dyprint*,'Time : ',ml_time,pl_timeprint*,'Pollon : ',ml_pollon,pl_pollonprint*,'Pollat : ',ml_pollat,pl_pollatstopendifc Allocate memory for all fieldsallocate(ml_ps(ml_nx,ml_ny),stat=stat)if (stat.ne.0) print*,'*** error allocating array ml_ps ***'allocate(ml_zb(ml_nx,ml_ny),stat=stat)if (stat.ne.0) print*,'*** error allocating array ml_zb ***'allocate(ml_p3(ml_nx,ml_ny,ml_nz),stat=stat)if (stat.ne.0) print*,'*** error allocating array ml_p3 ***'allocate(ml_t3(ml_nx,ml_ny,ml_nz),stat=stat)if (stat.ne.0) print*,'*** error allocating array ml_t3 ***'allocate(ml_q3(ml_nx,ml_ny,ml_nz),stat=stat)if (stat.ne.0) print*,'*** error allocating array ml_q3 ***'allocate(ml_tv3(ml_nx,ml_ny,ml_nz),stat=stat)if (stat.ne.0) print*,'*** error allocating array ml_tv3 ***'allocate(ml_z3(ml_nx,ml_ny,ml_nz),stat=stat)if (stat.ne.0) print*,'*** error allocating array ml_z3 ***'allocate(in_field(ml_nx,ml_ny,ml_nz),stat=stat)if (stat.ne.0) print*,'*** error allocating array in_field ***'allocate(pl_z3(pl_nx,pl_ny,pl_nz),stat=stat)if (stat.ne.0) print*,'*** error allocating array pl_z3 ***'allocate(pl_p3(pl_nx,pl_ny,pl_nz),stat=stat)if (stat.ne.0) print*,'*** error allocating array pl_p3 ***'allocate(zl_field(zl_nx,zl_ny,zl_nz),stat=stat)if (stat.ne.0) print*,'*** error allocating array zl_field ***'allocate(zl_p(zl_nx,zl_ny,zl_nz),stat=stat)if (stat.ne.0) print*,'*** error allocating array zl_p ***'c Read T, Q, PS from P filecall cdfopn(ml_pfn,cdfid,ierr)if (ierr.ne.0) goto 998isok=0varname='T'call check_varok (isok,varname, ml_vnam,ml_nvars)varname='Q'call check_varok (isok,varname, ml_vnam,ml_nvars)varname='PS'call check_varok (isok,varname,ml_vnam,ml_nvars)if (isok.ne.3) goto 998print*,'R T ',trim(ml_pfn)call getdat(cdfid,'T',ml_time,0,ml_t3,ierr)print*,'R Q ',trim(ml_pfn)call getdat(cdfid,'Q',ml_time,0,ml_q3,ierr)print*,'R PS ',trim(ml_pfn)call getdat(cdfid,'PS',ml_time,1,ml_ps,ierr)call clscdf(cdfid,ierr)c Read Z from Z filecall cdfopn(pl_zfn,cdfid,ierr)if (ierr.ne.0) goto 998isok=0varname='Z'call check_varok (isok,varname,pl_vnam,pl_nvars)if (isok.ne.1) goto 998print*,'R Z ',trim(pl_zfn)call getdat(cdfid,varname,pl_time,0,pl_z3,ierr)call clscdf(cdfid,ierr)c Set the values for the pressure on the pressure levelsdo i=1,pl_nxdo j=1,pl_nydo k=1,pl_nzpl_p3(i,j,k)=pl_aklay(k)enddoenddoenddoc -----------------------------------------------------------------c Calculate geopotential on layersc -----------------------------------------------------------------c Calculate 3d pressure fieldprint*,'C P'do k=1,ml_nzdo i=1,ml_nxdo j=1,ml_nyml_p3(i,j,k)=ml_aklay(k)+ml_bklay(k)*ml_ps(i,j)enddoenddoenddoc Calculate 3d virtual temperatureprint*,'C TV'do k=1,ml_nzdo i=1,ml_nxdo j=1,ml_nyml_tv3(i,j,k) = (ml_t3(i,j,k)+tzero)*> (1.+kappa*ml_q3(i,j,k))enddoenddoenddoc Loop over all grid pointsprint*,'C HYDROSTATIC EQUATION'do i=1,ml_nxdo j=1,ml_nyc Make the virtual temperature profile availabledo k=1,ml_nzp1 (ml_nz-k+1)=ml_p3 (i,j,k)tv1(ml_nz-k+1)=ml_tv3(i,j,k)enddocall spline (p1,tv1,ml_nz,1.e30,1.e30,spline_tv1)c Loop over all model levelsdo k=1,ml_nzc Get pressure at the grid pointp = ml_p3(i,j,k)c Find nearest pressure level which is above topographylmin=pl_nzdo l=1,pl_nzif ((abs(p-pl_p3(i,j,l))).lt.(abs(p-pl_p3(i,j,lmin)))> .and.> (pl_p3(i,j,l).lt.ml_ps(i,j)) ) thenlmin=lendifenddoc Integrate hydrostatic equation from this level to the grid pointp0 = pl_p3(i,j,lmin)n = nint(abs(p-p0)/dpmin)if (n.lt.1) n=1dp = (p-p0)/real(n)pu = p0z = pl_z3(i,j,lmin)call splint(p1,tv1,spline_tv1,ml_nz,pu,tvu)do l=1,npo = pu+dpcall splint(p1,tv1,spline_tv1,ml_nz,po,tvo)z = z + rdg*0.5*(tvu+tvo)*alog(pu/po)tvu = tvopu = poenddoc Set the geopotential at the grid pointml_z3(i,j,k) = zenddoenddoenddoc -----------------------------------------------------------------c Calculate height of topographyc -----------------------------------------------------------------print*,'C TOPOGRAPHY'do i=1,ml_nxdo j=1,ml_nyc Make the z/p profile availabledo k=1,ml_nzp1(ml_nz-k+1)=ml_p3(i,j,k)z1(ml_nz-k+1)=ml_z3(i,j,k)enddoc Cubic spline interpolationcall spline (p1,z1,ml_nz,1.e30,1.e30,spline_z1)call splint (p1,z1,spline_z1,ml_nz,ml_ps(i,j),ml_zb(i,j))enddoenddoc -----------------------------------------------------------------c Interpolate pressure onto z levelsc -----------------------------------------------------------------do i=1,zl_nxdo j=1,zl_nyc Make 1d profile availabledo k=1,ml_nzz1(k)=ml_z3(i,j,k)f1(k)=ml_p3(i,j,k)enddocall spline (z1,f1,ml_nz,1.e30,1.e30,spline_f1)do k=1,zl_nzif (zl_aklay(k).gt.ml_zb(i,j)) thencall splint(z1,f1,spline_f1,ml_nz,zl_aklay(k),ff)zl_p(i,j,k)=ffelsezl_p(i,j,k)=zl_mdvendifenddoenddoenddoc -----------------------------------------------------------------c Write output to netcdf filec -----------------------------------------------------------------if (test.eq.1) thenc Create the output file (grid taken from temperature)cdfname=testfilecall cdfopn(ml_pfn,cdfid,ierr)if (ierr.ne.0) goto 998call getcfn(cdfid,cfn,ierr)if (ierr.ne.0) goto 998call clscdf(cdfid,ierr)call crecdf(cdfname,cdfid,> ml_varmin,ml_varmax,ml_ndim,cfn,ierr)if (ierr.ne.0) goto 997c Write the definitions of the output variablesvarname='PS'ml_vardim(3)=1call putdef(cdfid,varname,ml_ndim,ml_mdv,> ml_vardim,ml_varmin,ml_varmax,ml_stag,ierr)if (ierr.ne.0) goto 997varname='ORO'call putdef(cdfid,varname,ml_ndim,ml_mdv,> ml_vardim,ml_varmin,ml_varmax,ml_stag,ierr)if (ierr.ne.0) goto 997ml_vardim(3)=ml_nzvarname='Z'call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim,> ml_varmin,ml_varmax,ml_stag,ierr)if (ierr.ne.0) goto 997c Write output fieldsprint*,'W PS TEST'varname='PS'call putdat(cdfid,varname,ml_time,1,ml_ps,ierr)if (ierr.ne.0) goto 997print*,'W ZB TEST'varname='ORO'call putdat(cdfid,varname,ml_time,1,ml_zb,ierr)if (ierr.ne.0) goto 997print*,'W Z TEST'varname='Z'call putdat(cdfid,varname,ml_time,0,ml_z3,ierr)if (ierr.ne.0) goto 997c Close cdf filecall clscdf(cdfid,ierr)endifc -----------------------------------------------------------------c Interpolate fields onto a stack of z levels and write new filec -----------------------------------------------------------------c Create output filecfn=trim(zl_ofn)//'_cst'call crecdf(trim(zl_ofn),cdfid,zl_varmin,zl_varmax,> zl_ndim,trim(cfn),ierr)if (ierr.ne.0) goto 995c Write topographyprint*,'W ZB ',trim(zl_ofn)varname='ORO'zl_vardim(3)=1call putdef(cdfid,varname,zl_ndim,zl_mdv,zl_vardim,> zl_varmin,zl_varmax,zl_stag,ierr)zl_vardim(3)=zl_nzif (ierr.ne.0) goto 995call putdat(cdfid,varname,zl_time,1,ml_zb,ierr)if (ierr.ne.0) goto 995c Write pressureprint*,'W P ',trim(zl_ofn)varname='P'call putdef(cdfid,varname,4,zl_mdv,zl_vardim,> zl_varmin,zl_varmax,zl_stag,ierr)if (ierr.ne.0) goto 995call putdat(cdfid,varname,zl_time,0,zl_p,ierr)if (ierr.ne.0) goto 995c Close output filecall clscdf(cdfid,ierr)c Loop over all variablesdo l=1,nvarprint*,'I ',trim(varinp(l))c Get grid description for variablecall cdfopn(varsrc(l),cdfid,ierr)if (ierr.ne.0) goto 996call getcfn(cdfid,cfn,ierr)if (ierr.ne.0) goto 996call cdfopn(cfn,cstid,ierr)if (ierr.ne.0) goto 996call getvars(cdfid,in_nvars,in_vnam,ierr)isok=0varname=varinp(l)call check_varok (isok,varname, in_vnam,in_nvars)if (isok.ne.1) goto 996call getdef(cdfid,varinp(l),in_ndim,in_mdv,in_vardim,> in_varmin,in_varmax,in_stag,ierr)if (ierr.ne.0) goto 996in_nx =in_vardim(1)in_ny =in_vardim(2)in_nz =in_vardim(3)in_xmin=in_varmin(1)in_ymin=in_varmin(2)call getlevs(cstid,in_nz,in_aklev,in_bklev,> in_aklay,in_bklay,ierr)call getgrid(cstid,in_dx,in_dy,ierr)in_xmax=in_xmin+real(in_nx-1)*in_dxin_ymax=in_ymin+real(in_ny-1)*in_dycall gettimes(cdfid,in_time,in_ntimes,ierr)call getstart(cstid,in_idate,ierr)call getpole(cstid,in_pollon,in_pollat,ierr)call clscdf(cstid,ierr)call clscdf(cdfid,ierr)c Check whether this grid is consistent with P gridif ( (ml_nx.ne.in_nx).or.> (ml_ny.ne.in_ny).or.> (ml_nz.ne.in_nz).or.> (abs(ml_xmin-in_xmin ).gt.eps).or.> (abs(ml_ymin-in_ymin ).gt.eps).or.> (abs(ml_xmax-in_xmax ).gt.eps).or.> (abs(ml_ymax-in_ymax ).gt.eps).or.> (abs(ml_dx -in_dx ).gt.eps).or.> (abs(ml_dy -in_dy ).gt.eps).or.> (abs(ml_time-in_time ).gt.eps).or.> (abs(ml_stag(1)-in_stag(1)).gt.eps).or.> (abs(ml_stag(2)-in_stag(2)).gt.eps).or.> (abs(ml_stag(3)-in_stag(3)).gt.eps).or.> (abs(ml_pollon-in_pollon ).gt.eps).or.> (abs(ml_pollat-in_pollat ).gt.eps)) thenprint*,trim(varinp(l)),> ' :Input P and FIELD grids are not consistent... Stop'print*,'Stag: ',ml_stag, in_stagprint*,'Pol: ',ml_pollon,in_pollon,ml_pollat,in_pollatstopendifc Read variable from filecall cdfopn(varsrc(l),cdfid,ierr)if (ierr.ne.0) goto 996call getdat(cdfid,varinp(l),in_time,0,in_field,ierr)if (ierr.ne.0) goto 996call clscdf(cdfid,ierr)c Write the constants fileif (l.eq.1) thendatar(1)=zl_nxdatar(2)=zl_nydatar(3)=int(1000.*zl_varmax(2))datar(4)=int(1000.*zl_varmin(1))datar(5)=int(1000.*zl_varmin(2))datar(6)=int(1000.*zl_varmax(1))datar(7)=int(1000.*zl_dx)datar(8)=int(1000.*zl_dy)datar(9)=zl_nzdatar(10)=1datar(11)=1datar(12)=0datar(13)=int(1000.*zl_pollon)datar(14)=int(1000.*zl_pollat)cfn=trim(zl_ofn)//'_cst'call wricst(cfn,datar,zl_aklev,zl_bklev,> zl_aklay,zl_bklay,zl_idate)endifc Do the interpolationdo i=1,zl_nxdo j=1,zl_nyc Make 1d profile availabledo k=1,ml_nzz1(k)=ml_z3 (i,j,k)f1(k)=in_field(i,j,k)enddocall spline (z1,f1,ml_nz,1.e30,1.e30,spline_f1)do k=1,zl_nzif (zl_aklay(k).gt.ml_zb(i,j)) thencall splint(z1,f1,spline_f1,ml_nz,zl_aklay(k),ff)zl_field(i,j,k)=ffelsezl_field(i,j,k)=zl_mdvendifenddoenddoenddoc Write the output fieldprint*,'W ',trim(varout(l)),' ',trim(zl_ofn)call cdfwopn(trim(zl_ofn),cdfid,ierr)varname=trim(varout(l))call putdef(cdfid,varname,4,zl_mdv,zl_vardim,> zl_varmin,zl_varmax,zl_stag,ierr)if (ierr.ne.0) goto 995call putdat(cdfid,varname,zl_time,0,zl_field,ierr)if (ierr.ne.0) goto 995call clscdf(cdfid,ierr)enddoc -----------------------------------------------------------------c Write topography and geopotential also to the input P filec -----------------------------------------------------------------c Open the P filecall cdfwopn(trim(ml_pfn),cdfid,ierr)c Write topographyvarname='ORO'print*,'W ',trim(varname),' ',trim(ml_pfn)isok=0call check_varok (isok,varname,ml_vnam,ml_nvars)if (isok.eq.0) thenml_vardim(3)=1call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim,> ml_varmin,ml_varmax,ml_stag,ierr)ml_vardim(3)=ml_nzif (ierr.ne.0) goto 997endifcall putdat(cdfid,varname,ml_time,1,ml_zb,ierr)if (ierr.ne.0) goto 997c Write geopotential heightvarname='Z'print*,'W ',trim(varname),' ',trim(ml_pfn)isok=0call check_varok (isok,varname,ml_vnam,ml_nvars)if (isok.eq.0) thencall putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim,> ml_varmin,ml_varmax,ml_stag,ierr)if (ierr.ne.0) goto 997endifcall putdat(cdfid,varname,ml_time,0,ml_z3,ierr)if (ierr.ne.0) goto 997c Close P filecall clscdf(cdfid,ierr)c -----------------------------------------------------------------c Exception handling and format specsc -----------------------------------------------------------------stop998 print*,'Z: Problems with input from m level'stop997 print*,'Z: Problems with output on m level'stop996 print*,'F: Problems with input from m level'stop995 print*,'F: Problems with output on z level'stopendc -------------------------------------------------------------c Natural cubic splinec -------------------------------------------------------------SUBROUTINE spline(x,y,n,yp1,ypn,y2)INTEGER n,NMAXREAL yp1,ypn,x(n),y(n),y2(n)PARAMETER (NMAX=500)INTEGER i,kREAL p,qn,sig,un,u(NMAX)if (yp1.gt..99e30) theny2(1)=0.u(1)=0.elsey2(1)=-0.5u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)endifdo 11 i=2,n-1sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))p=sig*y2(i-1)+2.y2(i)=(sig-1.)/pu(i)=(6.*((y(i+1)-y(i))/(x(i+*1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig**u(i-1))/p11 continueif (ypn.gt..99e30) thenqn=0.un=0.elseqn=0.5un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))endify2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)do 12 k=n-1,1,-1y2(k)=y2(k)*y2(k+1)+u(k)12 continuereturnENDSUBROUTINE splint(xa,ya,y2a,n,x,y)INTEGER nREAL x,y,xa(n),y2a(n),ya(n)INTEGER k,khi,kloREAL a,b,hklo=1khi=n1 if (khi-klo.gt.1) thenk=(khi+klo)/2if(xa(k).gt.x)thenkhi=kelseklo=kendifgoto 1endifh=xa(khi)-xa(klo)if (h.eq.0.) pause 'bad xa input in splint'a=(xa(khi)-x)/hb=(x-xa(klo))/hy=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h***2)/6.returnENDc ----------------------------------------------------------------c Check whether variable is found on netcdf filec ----------------------------------------------------------------subroutine check_varok (isok,varname,varlist,nvars)c Check whether the variable <varname> is in the list <varlist(nvars)>.c If this is the case, <isok> is incremented by 1. Otherwise <isok>c keeps its value.implicit nonec Declaraion of subroutine parametersinteger isokinteger nvarscharacter*80 varnamecharacter*80 varlist(nvars)c Auxiliary variablesinteger ic Maindo i=1,nvarsif (trim(varname).eq.trim(varlist(i))) isok=isok+1enddoend