| 0,0 → 1,846 |
|---|
| PROGRAM p2z |
| |
| c Calculate the geopotential and add it to the P file |
| c Michael Sprenger / Spring 2006 |
| |
| implicit none |
| |
| c --------------------------------------------------------------- |
| c Declaration of variables |
| c --------------------------------------------------------------- |
| |
| c Variables for input P file : model level |
| character*80 ml_pfn |
| real ml_varmin(4),ml_varmax(4),ml_stag(4) |
| integer ml_vardim(4) |
| real ml_mdv |
| integer ml_ndim |
| integer ml_nx,ml_ny,ml_nz |
| real ml_xmin,ml_xmax,ml_ymin,ml_ymax,ml_dx,ml_dy |
| integer ml_ntimes |
| real ml_aklev(500),ml_bklev(500) |
| real ml_aklay(500),ml_bklay(500) |
| real ml_time |
| real ml_pollon,ml_pollat |
| integer ml_nvars |
| character*80 ml_vnam(100) |
| integer ml_idate(5) |
| real,allocatable, dimension (:,:) :: ml_ps,ml_zb |
| real,allocatable, dimension (:,:,:) :: ml_t3,ml_q3,ml_p3,ml_tv3 |
| real,allocatable, dimension (:,:,:) :: ml_z3 |
| |
| c Variables for input Z file : pressure level |
| character*80 pl_zfn |
| real pl_varmin(4),pl_varmax(4),pl_stag(4) |
| integer pl_vardim(4) |
| real pl_mdv |
| integer pl_ndim |
| integer pl_nx,pl_ny,pl_nz |
| real pl_xmin,pl_xmax,pl_ymin,pl_ymax,pl_dx,pl_dy |
| integer pl_ntimes |
| real pl_aklev(500),pl_bklev(500) |
| real pl_aklay(500),pl_bklay(500) |
| real pl_time |
| real pl_pollon,pl_pollat |
| integer pl_nvars |
| character*80 pl_vnam(100) |
| integer pl_idate(5) |
| real,allocatable, dimension (:,:,:) :: pl_z3,pl_p3 |
| |
| c Variables for output Z file : height level |
| character*80 zl_ofn |
| real zl_varmin(4),zl_varmax(4),zl_stag(4) |
| integer zl_vardim(4) |
| real zl_mdv |
| integer zl_ndim |
| integer zl_nx,zl_ny,zl_nz |
| real zl_xmin,zl_xmax,zl_ymin,zl_ymax,zl_dx,zl_dy |
| integer zl_ntimes |
| real zl_aklev(500),zl_bklev(500) |
| real zl_aklay(500),zl_bklay(500) |
| real zl_time |
| real zl_pollon,zl_pollat |
| integer zl_idate(5) |
| real,allocatable, dimension (:,:,:) :: zl_field |
| real,allocatable, dimension (:,:,:) :: zl_p |
| |
| |
| c Variables for input P,S file : model level |
| character*80 in_ofn |
| real in_varmin(4),in_varmax(4),in_stag(4) |
| integer in_vardim(4) |
| real in_mdv |
| integer in_ndim |
| integer in_nx,in_ny,in_nz |
| real in_xmin,in_xmax,in_ymin,in_ymax,in_dx,in_dy |
| integer in_ntimes |
| real in_aklev(500),in_bklev(500) |
| real in_aklay(500),in_bklay(500) |
| real in_time |
| real in_pollon,in_pollat |
| integer in_nvars |
| character*80 in_vnam(100) |
| integer in_idate(5) |
| real,allocatable, dimension (:,:,:) :: in_field |
| |
| c Physical and numerical parameters |
| real g |
| parameter (g=9.80665) |
| real eps |
| parameter (eps=0.01) |
| real tzero |
| parameter (tzero=273.15) |
| real kappa |
| parameter (kappa=0.6078) |
| real zerodiv |
| parameter (zerodiv=0.0000000001) |
| real dpmin |
| parameter (dpmin=10.) |
| real rdg |
| parameter (rdg=29.271) |
| |
| c Flag for test mode |
| integer test |
| parameter (test=0) |
| character*80 testfile |
| parameter (testfile='TEST') |
| |
| c Variables and levels for interpolation onto z levels |
| character*80 levelfile,varfile |
| integer nvar,nlev |
| character*80 varinp(100),varout(100),varsrc(100) |
| real zlev(500) |
| |
| c Auxiliary variables |
| integer ierr |
| integer cdfid,cstid |
| character*80 cfn |
| integer stat |
| real time |
| real 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,ff |
| integer i,j,k,l |
| integer lmin,n |
| character*80 varname,cdfname |
| integer idate(5),stdate(5),datar(14) |
| integer isok |
| character*80 name |
| real zmin,dz |
| |
| c ----------------------------------------------------------------- |
| c Read input fields |
| c ----------------------------------------------------------------- |
| |
| print*,'*********************************************************' |
| print*,'* p2z *' |
| print*,'*********************************************************' |
| |
| c Read in the parameter file |
| open(10,file='fort.10') |
| |
| read(10,*) ml_pfn |
| read(10,*) pl_zfn |
| read(10,*) zl_ofn |
| |
| read(10,*) name,zmin |
| if ( trim(name).ne.'GEO_ZMIN') stop |
| read(10,*) name,nlev |
| if ( trim(name).ne.'GEO_NZ' ) stop |
| read(10,*) name,dz |
| if ( trim(name).ne.'GEO_DZ' ) stop |
| do i=1,nlev |
| zlev(i)=zmin+real(i-1)*dz |
| enddo |
| |
| nvar=1 |
| 102 continue |
| read(10,*,end=103) varinp(nvar),varout(nvar),varsrc(nvar) |
| nvar=nvar+1 |
| goto 102 |
| 103 continue |
| nvar=nvar-1 |
| print* |
| do i=1,nvar |
| write(*,'(a10,a10,a30)') trim(varinp(i)),trim(varout(i)), |
| > trim(varsrc(i)) |
| enddo |
| print* |
| |
| close(10) |
| |
| |
| c Get grid description for P file : model level |
| call cdfopn(ml_pfn,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getcfn(cdfid,cfn,ierr) |
| if (ierr.ne.0) goto 998 |
| call cdfopn(cfn,cstid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getvars(cdfid,ml_nvars,ml_vnam,ierr) |
| varname='T' |
| isok=0 |
| call check_varok (isok,varname,ml_vnam,ml_nvars) |
| if (isok.ne.1) goto 998 |
| call getdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim, |
| > ml_varmin,ml_varmax,ml_stag,ierr) |
| if (ierr.ne.0) goto 998 |
| ml_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_dx |
| ml_ymax=ml_ymin+real(ml_ny-1)*ml_dy |
| call 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 level |
| call cdfopn(pl_zfn,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getcfn(cdfid,cfn,ierr) |
| if (ierr.ne.0) goto 998 |
| call cdfopn(cfn,cstid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getvars(cdfid,pl_nvars,pl_vnam,ierr) |
| varname='Z' |
| isok=0 |
| call check_varok (isok,varname,pl_vnam,pl_nvars) |
| if (isok.ne.1) goto 998 |
| call getdef(cdfid,varname,pl_ndim,pl_mdv,pl_vardim, |
| > pl_varmin,pl_varmax,pl_stag,ierr) |
| if (ierr.ne.0) goto 998 |
| pl_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_dx |
| pl_ymax=pl_ymin+real(pl_ny-1)*pl_dy |
| call 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 level |
| zl_vardim(1) = pl_vardim(1) |
| zl_vardim(2) = pl_vardim(2) |
| zl_vardim(3) = nlev |
| zl_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,nlev |
| zl_aklay(i) = zlev(i) |
| zl_bklay(i) = 0. |
| zl_aklev(i) = zlev(i) |
| zl_bklev(i) = 0. |
| enddo |
| zl_dx = pl_dx |
| zl_dy = pl_dy |
| zl_time = pl_time |
| zl_ntimes = pl_ntimes |
| zl_ndim = 4 |
| zl_mdv = pl_mdv |
| zl_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_pollon |
| zl_pollat = ml_pollat |
| do i=1,5 |
| zl_idate(i) = ml_idate(i) |
| enddo |
| |
| c Consitency check for the grids |
| if ( (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)) then |
| print*,'Input P and Z grids are not consistent... Stop' |
| print* |
| print*,'Xmin : ',ml_xmin,pl_xmin |
| print*,'Ymin : ',ml_ymin,pl_ymin |
| print*,'Xmax : ',ml_xmax,pl_xmax |
| print*,'Ymax : ',ml_ymax,pl_ymax |
| print*,'Dx : ',ml_dx, pl_dx |
| print*,'Dy : ',ml_dy, pl_dy |
| print*,'Time : ',ml_time,pl_time |
| print*,'Pollon : ',ml_pollon,pl_pollon |
| print*,'Pollat : ',ml_pollat,pl_pollat |
| stop |
| endif |
| |
| c Allocate memory for all fields |
| allocate(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 file |
| call cdfopn(ml_pfn,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| isok=0 |
| varname='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 998 |
| print*,'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 file |
| call cdfopn(pl_zfn,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| isok=0 |
| varname='Z' |
| call check_varok (isok,varname,pl_vnam,pl_nvars) |
| if (isok.ne.1) goto 998 |
| print*,'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 levels |
| do i=1,pl_nx |
| do j=1,pl_ny |
| do k=1,pl_nz |
| pl_p3(i,j,k)=pl_aklay(k) |
| enddo |
| enddo |
| enddo |
| |
| c ----------------------------------------------------------------- |
| c Calculate geopotential on layers |
| c ----------------------------------------------------------------- |
| |
| c Calculate 3d pressure field |
| print*,'C P' |
| do k=1,ml_nz |
| do i=1,ml_nx |
| do j=1,ml_ny |
| ml_p3(i,j,k)=ml_aklay(k)+ml_bklay(k)*ml_ps(i,j) |
| enddo |
| enddo |
| enddo |
| |
| c Calculate 3d virtual temperature |
| print*,'C TV' |
| do k=1,ml_nz |
| do i=1,ml_nx |
| do j=1,ml_ny |
| ml_tv3(i,j,k) = (ml_t3(i,j,k)+tzero)* |
| > (1.+kappa*ml_q3(i,j,k)) |
| enddo |
| enddo |
| enddo |
| |
| c Loop over all grid points |
| print*,'C HYDROSTATIC EQUATION' |
| do i=1,ml_nx |
| do j=1,ml_ny |
| |
| c Make the virtual temperature profile available |
| do k=1,ml_nz |
| p1 (ml_nz-k+1)=ml_p3 (i,j,k) |
| tv1(ml_nz-k+1)=ml_tv3(i,j,k) |
| enddo |
| call spline (p1,tv1,ml_nz,1.e30,1.e30,spline_tv1) |
| |
| c Loop over all model levels |
| do k=1,ml_nz |
| |
| c Get pressure at the grid point |
| p = ml_p3(i,j,k) |
| |
| c Find nearest pressure level which is above topography |
| lmin=pl_nz |
| do l=1,pl_nz |
| if ((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)) ) then |
| lmin=l |
| endif |
| enddo |
| |
| c Integrate hydrostatic equation from this level to the grid point |
| p0 = pl_p3(i,j,lmin) |
| n = nint(abs(p-p0)/dpmin) |
| if (n.lt.1) n=1 |
| dp = (p-p0)/real(n) |
| |
| pu = p0 |
| z = pl_z3(i,j,lmin) |
| call splint(p1,tv1,spline_tv1,ml_nz,pu,tvu) |
| do l=1,n |
| po = pu+dp |
| call splint(p1,tv1,spline_tv1,ml_nz,po,tvo) |
| z = z + rdg*0.5*(tvu+tvo)*alog(pu/po) |
| tvu = tvo |
| pu = po |
| enddo |
| |
| c Set the geopotential at the grid point |
| ml_z3(i,j,k) = z |
| |
| enddo |
| |
| enddo |
| enddo |
| |
| c ----------------------------------------------------------------- |
| c Calculate height of topography |
| c ----------------------------------------------------------------- |
| |
| print*,'C TOPOGRAPHY' |
| do i=1,ml_nx |
| do j=1,ml_ny |
| |
| c Make the z/p profile available |
| do k=1,ml_nz |
| p1(ml_nz-k+1)=ml_p3(i,j,k) |
| z1(ml_nz-k+1)=ml_z3(i,j,k) |
| enddo |
| |
| c Cubic spline interpolation |
| call 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)) |
| |
| enddo |
| enddo |
| |
| c ----------------------------------------------------------------- |
| c Interpolate pressure onto z levels |
| c ----------------------------------------------------------------- |
| |
| do i=1,zl_nx |
| do j=1,zl_ny |
| |
| c Make 1d profile available |
| do k=1,ml_nz |
| z1(k)=ml_z3(i,j,k) |
| f1(k)=ml_p3(i,j,k) |
| enddo |
| call spline (z1,f1,ml_nz,1.e30,1.e30,spline_f1) |
| do k=1,zl_nz |
| if (zl_aklay(k).gt.ml_zb(i,j)) then |
| call splint(z1,f1,spline_f1,ml_nz,zl_aklay(k),ff) |
| zl_p(i,j,k)=ff |
| else |
| zl_p(i,j,k)=zl_mdv |
| endif |
| enddo |
| |
| enddo |
| enddo |
| |
| c ----------------------------------------------------------------- |
| c Write output to netcdf file |
| c ----------------------------------------------------------------- |
| |
| if (test.eq.1) then |
| |
| c Create the output file (grid taken from temperature) |
| cdfname=testfile |
| |
| call cdfopn(ml_pfn,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getcfn(cdfid,cfn,ierr) |
| if (ierr.ne.0) goto 998 |
| call clscdf(cdfid,ierr) |
| |
| call crecdf(cdfname,cdfid, |
| > ml_varmin,ml_varmax,ml_ndim,cfn,ierr) |
| if (ierr.ne.0) goto 997 |
| |
| c Write the definitions of the output variables |
| varname='PS' |
| ml_vardim(3)=1 |
| call putdef(cdfid,varname,ml_ndim,ml_mdv, |
| > ml_vardim,ml_varmin,ml_varmax,ml_stag,ierr) |
| if (ierr.ne.0) goto 997 |
| varname='ORO' |
| call putdef(cdfid,varname,ml_ndim,ml_mdv, |
| > ml_vardim,ml_varmin,ml_varmax,ml_stag,ierr) |
| if (ierr.ne.0) goto 997 |
| ml_vardim(3)=ml_nz |
| varname='Z' |
| call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim, |
| > ml_varmin,ml_varmax,ml_stag,ierr) |
| if (ierr.ne.0) goto 997 |
| |
| c Write output fields |
| print*,'W PS TEST' |
| varname='PS' |
| call putdat(cdfid,varname,ml_time,1,ml_ps,ierr) |
| if (ierr.ne.0) goto 997 |
| print*,'W ZB TEST' |
| varname='ORO' |
| call putdat(cdfid,varname,ml_time,1,ml_zb,ierr) |
| if (ierr.ne.0) goto 997 |
| print*,'W Z TEST' |
| varname='Z' |
| call putdat(cdfid,varname,ml_time,0,ml_z3,ierr) |
| if (ierr.ne.0) goto 997 |
| |
| c Close cdf file |
| call clscdf(cdfid,ierr) |
| |
| endif |
| |
| c ----------------------------------------------------------------- |
| c Interpolate fields onto a stack of z levels and write new file |
| c ----------------------------------------------------------------- |
| |
| c Create output file |
| cfn=trim(zl_ofn)//'_cst' |
| call crecdf(trim(zl_ofn),cdfid,zl_varmin,zl_varmax, |
| > zl_ndim,trim(cfn),ierr) |
| if (ierr.ne.0) goto 995 |
| |
| c Write topography |
| print*,'W ZB ',trim(zl_ofn) |
| varname='ORO' |
| zl_vardim(3)=1 |
| call putdef(cdfid,varname,zl_ndim,zl_mdv,zl_vardim, |
| > zl_varmin,zl_varmax,zl_stag,ierr) |
| zl_vardim(3)=zl_nz |
| if (ierr.ne.0) goto 995 |
| call putdat(cdfid,varname,zl_time,1,ml_zb,ierr) |
| if (ierr.ne.0) goto 995 |
| |
| c Write pressure |
| print*,'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 995 |
| call putdat(cdfid,varname,zl_time,0,zl_p,ierr) |
| if (ierr.ne.0) goto 995 |
| |
| c Close output file |
| call clscdf(cdfid,ierr) |
| |
| c Loop over all variables |
| do l=1,nvar |
| |
| print*,'I ',trim(varinp(l)) |
| |
| c Get grid description for variable |
| call cdfopn(varsrc(l),cdfid,ierr) |
| if (ierr.ne.0) goto 996 |
| call getcfn(cdfid,cfn,ierr) |
| if (ierr.ne.0) goto 996 |
| call cdfopn(cfn,cstid,ierr) |
| if (ierr.ne.0) goto 996 |
| call getvars(cdfid,in_nvars,in_vnam,ierr) |
| isok=0 |
| varname=varinp(l) |
| call check_varok (isok,varname, in_vnam,in_nvars) |
| if (isok.ne.1) goto 996 |
| call getdef(cdfid,varinp(l),in_ndim,in_mdv,in_vardim, |
| > in_varmin,in_varmax,in_stag,ierr) |
| if (ierr.ne.0) goto 996 |
| in_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_dx |
| in_ymax=in_ymin+real(in_ny-1)*in_dy |
| call 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 grid |
| if ( (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)) then |
| print*,trim(varinp(l)), |
| > ' :Input P and FIELD grids are not consistent... Stop' |
| print*,'Stag: ',ml_stag, in_stag |
| print*,'Pol: ',ml_pollon,in_pollon,ml_pollat,in_pollat |
| stop |
| endif |
| |
| c Read variable from file |
| call cdfopn(varsrc(l),cdfid,ierr) |
| if (ierr.ne.0) goto 996 |
| call getdat(cdfid,varinp(l),in_time,0,in_field,ierr) |
| if (ierr.ne.0) goto 996 |
| call clscdf(cdfid,ierr) |
| |
| c Write the constants file |
| if (l.eq.1) then |
| |
| datar(1)=zl_nx |
| datar(2)=zl_ny |
| datar(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_nz |
| datar(10)=1 |
| datar(11)=1 |
| datar(12)=0 |
| datar(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) |
| |
| endif |
| |
| c Do the interpolation |
| do i=1,zl_nx |
| do j=1,zl_ny |
| |
| c Make 1d profile available |
| do k=1,ml_nz |
| z1(k)=ml_z3 (i,j,k) |
| f1(k)=in_field(i,j,k) |
| enddo |
| call spline (z1,f1,ml_nz,1.e30,1.e30,spline_f1) |
| do k=1,zl_nz |
| if (zl_aklay(k).gt.ml_zb(i,j)) then |
| call splint(z1,f1,spline_f1,ml_nz,zl_aklay(k),ff) |
| zl_field(i,j,k)=ff |
| else |
| zl_field(i,j,k)=zl_mdv |
| endif |
| enddo |
| |
| enddo |
| enddo |
| |
| c Write the output field |
| print*,'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 995 |
| call putdat(cdfid,varname,zl_time,0,zl_field,ierr) |
| if (ierr.ne.0) goto 995 |
| call clscdf(cdfid,ierr) |
| |
| enddo |
| |
| c ----------------------------------------------------------------- |
| c Write topography and geopotential also to the input P file |
| c ----------------------------------------------------------------- |
| |
| c Open the P file |
| call cdfwopn(trim(ml_pfn),cdfid,ierr) |
| |
| c Write topography |
| varname='ORO' |
| print*,'W ',trim(varname),' ',trim(ml_pfn) |
| isok=0 |
| call check_varok (isok,varname,ml_vnam,ml_nvars) |
| if (isok.eq.0) then |
| ml_vardim(3)=1 |
| call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim, |
| > ml_varmin,ml_varmax,ml_stag,ierr) |
| ml_vardim(3)=ml_nz |
| if (ierr.ne.0) goto 997 |
| endif |
| call putdat(cdfid,varname,ml_time,1,ml_zb,ierr) |
| if (ierr.ne.0) goto 997 |
| |
| c Write geopotential height |
| varname='Z' |
| print*,'W ',trim(varname),' ',trim(ml_pfn) |
| isok=0 |
| call check_varok (isok,varname,ml_vnam,ml_nvars) |
| if (isok.eq.0) then |
| call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim, |
| > ml_varmin,ml_varmax,ml_stag,ierr) |
| if (ierr.ne.0) goto 997 |
| endif |
| call putdat(cdfid,varname,ml_time,0,ml_z3,ierr) |
| if (ierr.ne.0) goto 997 |
| |
| c Close P file |
| call clscdf(cdfid,ierr) |
| |
| c ----------------------------------------------------------------- |
| c Exception handling and format specs |
| c ----------------------------------------------------------------- |
| |
| stop |
| |
| 998 print*,'Z: Problems with input from m level' |
| stop |
| |
| 997 print*,'Z: Problems with output on m level' |
| stop |
| |
| 996 print*,'F: Problems with input from m level' |
| stop |
| |
| 995 print*,'F: Problems with output on z level' |
| stop |
| |
| |
| end |
| |
| c ------------------------------------------------------------- |
| c Natural cubic spline |
| c ------------------------------------------------------------- |
| |
| SUBROUTINE spline(x,y,n,yp1,ypn,y2) |
| INTEGER n,NMAX |
| REAL yp1,ypn,x(n),y(n),y2(n) |
| PARAMETER (NMAX=500) |
| INTEGER i,k |
| REAL p,qn,sig,un,u(NMAX) |
| if (yp1.gt..99e30) then |
| y2(1)=0. |
| u(1)=0. |
| else |
| y2(1)=-0.5 |
| u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1) |
| endif |
| do 11 i=2,n-1 |
| sig=(x(i)-x(i-1))/(x(i+1)-x(i-1)) |
| p=sig*y2(i-1)+2. |
| y2(i)=(sig-1.)/p |
| u(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))/p |
| 11 continue |
| if (ypn.gt..99e30) then |
| qn=0. |
| un=0. |
| else |
| qn=0.5 |
| un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1))) |
| endif |
| y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.) |
| do 12 k=n-1,1,-1 |
| y2(k)=y2(k)*y2(k+1)+u(k) |
| 12 continue |
| return |
| END |
| |
| |
| SUBROUTINE splint(xa,ya,y2a,n,x,y) |
| INTEGER n |
| REAL x,y,xa(n),y2a(n),ya(n) |
| INTEGER k,khi,klo |
| REAL a,b,h |
| klo=1 |
| khi=n |
| 1 if (khi-klo.gt.1) then |
| k=(khi+klo)/2 |
| if(xa(k).gt.x)then |
| khi=k |
| else |
| klo=k |
| endif |
| goto 1 |
| endif |
| h=xa(khi)-xa(klo) |
| if (h.eq.0.) pause 'bad xa input in splint' |
| a=(xa(khi)-x)/h |
| b=(x-xa(klo))/h |
| y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h** |
| *2)/6. |
| return |
| END |
| |
| c ---------------------------------------------------------------- |
| c Check whether variable is found on netcdf file |
| c ---------------------------------------------------------------- |
| |
| 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 none |
| |
| c Declaraion of subroutine parameters |
| integer isok |
| integer nvars |
| character*80 varname |
| character*80 varlist(nvars) |
| |
| c Auxiliary variables |
| integer i |
| |
| c Main |
| do i=1,nvars |
| if (trim(varname).eq.trim(varlist(i))) isok=isok+1 |
| enddo |
| |
| end |
| Property changes: |
|---|
| Added: svn:executable |
|---|