Subversion Repositories pvinversion.ecmwf

      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