Subversion Repositories pvinversion.ecmwf

      PROGRAM inv_prep

c     ********************************************************************************
c     * CALCULATE REFERENCE PROFILE, CORIOLIS PARAMETER AND GRID PARAMETERS          *
c     * Rene Fehlmann 1994 / Code re-organization: Michael Sprenger, 2006            *
c     ********************************************************************************

c     --------------------------------------------------------------------------------
c     Declaration of variables, parameters, externals and common blocks
c     --------------------------------------------------------------------------------

      implicit none

c     Input and output file
      character*80   pvsrcfile,referfile
      integer        mode
      real           radius

c     Grid parameters
      integer        nx,ny,nz
      real           xmin,ymin,zmin,xmax,ymax,zmax
      real           dx,dy,dz
      real           mdv

c     Reference state
      real,   allocatable,dimension (:) :: nsqref
      real,   allocatable,dimension (:) :: thetaref
      real,   allocatable,dimension (:) :: rhoref
      real,   allocatable,dimension (:) :: pressref
      real,   allocatable,dimension (:) :: zref
      real    pressn,thetan

c     3d fields for calculation of reference profile
      real,allocatable,dimension (:,:,:) :: th,rho,nsq,p,z

c     2d weight for mean
      real,allocatable,dimension (:,:) :: weight

c     Auxiliary variables
      integer      i,j,k
      integer      stat
      integer      jj,kk
      character*80 varname
      integer      istep
      real         sum,max,min
      integer      cnt,nmd
      integer      step
      integer      i0,j0
      real         lon0,lat0,lon1,lat1,dist
      integer      vardim(4),ndim,cdfid,ierr
      real         varmin(4),varmax(4),stag(4)
      real         mdv
      character*80 name

c     Externals
      real      sdis
      external  sdis

c     --------------------------------------------------------------------------------
c     Input
c     --------------------------------------------------------------------------------

      print*,'********************************************************'
      print*,'* ref_grid                                             *'
      print*,'********************************************************'

c     Read parameter file
      open(10,file='fort.10')
       read(10,*) pvsrcfile
       read(10,*) referfile
       read(10,*) name,radius
       if ( trim(name).ne.'REF_R') stop
      close(10) 
      print*
      print*,trim(pvsrcfile)
      print*,trim(referfile)
      print*,radius
      print*

c     Get lat/lon gid parameters from input file
      call read_dim (nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv,
     >               pvsrcfile)
      print*,'Read_Dim: nx,ny,nz         = ',nx,ny,nz
      print*,'          dx,dy,dz         = ',dx,dy,dz
      print*,'          xmin,ymin,zmin   = ',xmin,ymin,zmin
      print*,'          mdv              = ',mdv
      print*
      xmax = xmin + real(nx-1) * dx
      ymax = ymin + real(ny-1) * dy

c     Count from 0, not from 1: consistent with <inv_cart.f>.
      nx=nx-1
      ny=ny-1
      nz=nz-1

c     Allocate memory for reference profile
      allocate(rhoref  (0:2*nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating rhoref'
      allocate(pressref(0:2*nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating pressref'
      allocate(thetaref(0:2*nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating thetaref'
      allocate(nsqref  (0:2*nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating nsqref'
      allocate(zref    (0:2*nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating zref'
         
c     Allocate memory for calculatation of reference profile
      allocate(th (0:nx,0:ny,0:nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating th'
      allocate(rho(0:nx,0:ny,0:nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating rho'
      allocate(p  (0:nx,0:ny,0:nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating p'
      allocate(nsq(0:nx,0:ny,0:nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating nsq'
      allocate(z(0:nx,0:ny,0:nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating z'

c     Allocate memory for weight
      allocate(weight(0:nx,0:ny),STAT=stat)
      if (stat.ne.0) print*,'error allocating weight'

c     --------------------------------------------------------------------------------
c     Calculate the reference profile and put it onto netcdf file
c     --------------------------------------------------------------------------------

c     Read data from file
      varname='TH'
      call read_inp (th,varname,pvsrcfile,
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
      varname='RHO'
      call read_inp (rho,varname,pvsrcfile,
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
      varname='NSQ'
      call read_inp (nsq,varname,pvsrcfile,
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
      varname='P'
      call read_inp (p,varname,pvsrcfile,
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)

c     Init the height field (not really necessary, but code becomes more symmetrical)
      do i=0,nx
         do j=0,ny
            do k=0,nz
               z(i,j,k)=zmin+real(k)*dz
            enddo
         enddo
      enddo

c     Define the weight
      if ( radius.lt.0 ) then
        do i=0,nx
          do j=0,ny
              weight(i,j) = 1.
          enddo
        enddo
      else
        i0 = nx/2
        j0 = ny/2
        lon0 = xmin + real(i0-1) * dx
        lat0 = ymin + real(j0-1) * dy
        weight(i0,j0)=1.
        do i=0,nx
          do j=0,ny
             lon1 = xmin + real(i-1) * dx
             lat1 = ymin + real(j-1) * dy
             dist = sdis(lon0,lat0,lon1,lat1)
             if ( dist.lt.radius ) then
                weight(i,j) = 1.
             else
                weight(i,j) = 0.
             endif
          enddo
        enddo
      endif

c     Determine the reference profile (mean over domain, split levels and layers)
      call mean(zref,    z,  nx,ny,nz,mdv,weight)
      call mean(nsqref,  nsq,nx,ny,nz,mdv,weight)
      call mean(rhoref,  rho,nx,ny,nz,mdv,weight)
      call mean(thetaref,th, nx,ny,nz,mdv,weight)
      call mean(pressref,p,  nx,ny,nz,mdv,weight)

c     Write reference file to netcdf file
      call write_ref (nsqref,rhoref,thetaref,pressref,zref,
     >                nz,referfile)      

c     Write some info
      print*
      print*,'Ref:            z         p       rho       nsq     theta'
      step=2*nz/10
      if (step.lt.1) step=1
      do k=0,2*nz,step
         write(*,'(8x,f10.1,2f10.2,f10.6,f10.2)')
     >         zref(k),pressref(k),rhoref(k),nsqref(k),thetaref(k) 
      enddo
      print*

c     Write weighht to REF file
      call cdfwopn(trim(referfile),cdfid,ierr)
      varname='WEIGHT'
      vardim(1)=nx+1
      vardim(2)=ny+1
      vardim(3)=1
      vardim(4)=1
      varmin(1)=xmin
      varmin(2)=ymin
      varmin(3)=0.
      varmax(1)=xmax
      varmax(2)=ymax
      varmax(3)=0.
      ndim=3
      mdv=-999.
      call putdef(cdfid,varname,ndim,mdv,vardim,
     >            varmin,varmax,stag,ierr)
      call putdat(cdfid,varname,0.,1,weight,ierr)

c     --------------------------------------------------------------------------------
c     Format specifications
c     --------------------------------------------------------------------------------

 111  format (5f20.9)
 106  format (2f20.9)
          
      end

c     --------------------------------------------------------------------------
c     Spherical distance between lat/lon points
c     --------------------------------------------------------------------------

      real function sdis(xp,yp,xq,yq)
c
c     calculates spherical distance (in km) between two points given
c     by their spherical coordinates (xp,yp) and (xq,yq), respectively.
c
      real      re
      parameter (re=6370.)
      real      pi180
      parameter (pi180=3.14159/180.)
      real      xp,yp,xq,yq,arg

      arg=sin(pi180*yp)*sin(pi180*yq)+
     >    cos(pi180*yp)*cos(pi180*yq)*cos(pi180*(xp-xq))
      if (arg.lt.-1.) arg=-1.
      if (arg.gt.1.) arg=1.

      sdis=re*acos(arg)

      end

c     ********************************************************************************
c     * NETCDF INPUT AND OUTPUT                                                      *
c     ********************************************************************************


c     --------------------------------------------------------------------------------
c     Write reference file to netcdf
c     --------------------------------------------------------------------------------

      SUBROUTINE write_ref (nsqref,rhoref,thetaref,pressref,zref,
     >                      nz,pvsrcfile) 

c     Write the reference profile to file
c
c         thetaref             : Reference potential temperature (K)
c         pressref             : Reference pressure (Pa)
c         rhoref               : Reference density (kg/m^3)
c         nsqref               : Stratification (s^-1)
c         zref                 : Reference height (m)
c         nz                   : Grid dimension in z direction
c         pvsrcfile            : Output file

      implicit   none

c     Declaration of subroutine parameters
      integer              nx,ny,nz
      real                 nsqref  (0:2*nz)
      real                 thetaref(0:2*nz)
      real                 rhoref  (0:2*nz)
      real                 pressref(0:2*nz)
      real                 zref    (0:2*nz)
      character*80         pvsrcfile

c     Numerical and physical parameters
      real                 eps
      parameter            (eps=0.01)

c     Auxiliary variables
      integer              cdfid,stat
      integer              vardim(4)
      real                 misdat
      integer              ndimin
      real                 varmin(4),varmax(4),stag(4)
      integer              i,j,k,nf1
      integer              ntimes
      real                 time(200)  
      character*80         vnam(100),varname
      integer              nvars
      integer              isok,ierr

c     Get grid description from topography
      call cdfopn(pvsrcfile,cdfid,stat)
      if (stat.ne.0) goto 997
      call getvars(cdfid,nvars,vnam,stat)
      if (stat.ne.0) goto 997
      isok=0
      varname='ORO'
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) goto 997
      call getdef(cdfid,varname,ndimin,misdat,vardim,
     >            varmin,varmax,stag,stat)    
      if (stat.ne.0) goto 997
      time(1)=0.
      call gettimes(cdfid,time,ntimes,stat)  
      if (stat.ne.0) goto 997
      call clscdf(cdfid,stat)
      if (stat.ne.0) goto 997

c     Open output netcdf file
      call cdfwopn(pvsrcfile,cdfid,stat)
      if (stat.ne.0) goto 997
      
c     Create the variable if necessary
      call getvars(cdfid,nvars,vnam,stat)
      if (stat.ne.0) goto 997
      isok=0
      varname='NSQREF'
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) then
         vardim(1)=1
         vardim(2)=1
         vardim(3)=2*nz+1
         call putdef(cdfid,varname,ndimin,misdat,vardim,
     >               varmin,varmax,stag,stat)
         if (stat.ne.0) goto 997
      endif
      isok=0
      varname='RHOREF'
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) then
         vardim(1)=1
         vardim(2)=1
         vardim(3)=2*nz+1
         call putdef(cdfid,varname,ndimin,misdat,vardim,
     >               varmin,varmax,stag,stat)
         if (stat.ne.0) goto 997
      endif
      isok=0
      varname='PREREF'
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) then
         vardim(1)=1
         vardim(2)=1
         vardim(3)=2*nz+1
         call putdef(cdfid,varname,ndimin,misdat,vardim,
     >               varmin,varmax,stag,stat)
         if (stat.ne.0) goto 997
      endif
      isok=0
      varname='THETAREF'
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) then
         vardim(1)=1
         vardim(2)=1
         vardim(3)=2*nz+1
         call putdef(cdfid,varname,ndimin,misdat,vardim,
     >               varmin,varmax,stag,stat)
         if (stat.ne.0) goto 997
      endif
      isok=0
      varname='ZREF'
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) then
         vardim(1)=1
         vardim(2)=1
         vardim(3)=2*nz+1
         call putdef(cdfid,varname,ndimin,misdat,vardim,
     >               varmin,varmax,stag,stat)
         if (stat.ne.0) goto 997
      endif

c     Write data
      varname='NSQREF'
      print*,'W NSQREF     ',trim(pvsrcfile)
      call putdat(cdfid,varname,time(1),0,nsqref,stat)
      if (stat.ne.0) goto 997
      varname='RHOREF'
      print*,'W RHOREF     ',trim(pvsrcfile)
      call putdat(cdfid,varname,time(1),0,rhoref,stat)
      if (stat.ne.0) goto 997
      varname='THETAREF'
      print*,'W THETAREF   ',trim(pvsrcfile)
      call putdat(cdfid,varname,time(1),0,thetaref,stat)
      if (stat.ne.0) goto 997
      varname='PREREF'
      print*,'W PREREF     ',trim(pvsrcfile)
      call putdat(cdfid,varname,time(1),0,pressref,stat)
      if (stat.ne.0) goto 997
      varname='ZREF'
      print*,'W ZREF       ',trim(pvsrcfile)
      call putdat(cdfid,varname,time(1),0,zref,stat)
      if (stat.ne.0) goto 997

c     Close output netcdf file
      call clscdf(cdfid,stat)
      if (stat.ne.0) goto 997

      return

c     Exception handling
 997  print*,'Write_Ref: Problem with input netcdf file... Stop'
      stop

      end


c     --------------------------------------------------------------------------------
c     Read input fields for reference profile
c     --------------------------------------------------------------------------------

      SUBROUTINE read_inp (field,fieldname,pvsrcfile,
     >                     nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)

c     Read <fieldname> from netcdf file <pvsrcfile> into <field>. The grid is specified 
c     by <nx,ny,nz,dx,dy,dz,xmin,ymin,zmin>. A check is performed whether the input 
c     files are consitent with this grid. The missing data value is set to <mdv>.

      implicit   none

c     Declaration of subroutine parameters
      integer              nx,ny,nz
      real                 field(0:nx,0:ny,0:nz)
      character*80         fieldname
      character*80         pvsrcfile
      real                 dx,dy,dz,xmin,ymin,zmin
      real                 mdv

c     Numerical and physical parameters
      real                 eps
      parameter            (eps=0.01)

c     Auxiliary variables
      integer              cdfid,stat,cdfid99
      integer              vardim(4)
      real                 misdat
      real                 varmin(4),varmax(4),stag(4)
      integer              ndimin,outid,i,j,k
      real                 max_th
      real                 tmp(nx,ny,nz)
      integer              ntimes
      real                 time(200)       
      integer              nvars
      character*80         vnam(100),varname
      integer              isok

c     Open the input netcdf file
      call cdfopn(pvsrcfile,cdfid,stat)
      if (stat.ne.0) goto 998

c     Check whether needed variables are on file
      call getvars(cdfid,nvars,vnam,stat)
      if (stat.ne.0) goto 998
      isok=0
      varname=trim(fieldname)
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) goto 998

c     Get the grid parameters 
      call getdef(cdfid,varname,ndimin,misdat,vardim,
     >            varmin,varmax,stag,stat)    
      if (stat.ne.0) goto 998
      time(1)=0.
      call gettimes(cdfid,time,ntimes,stat)  
      if (stat.ne.0) goto 998

c     Check whether grid parameters are consistent
      if ( (vardim(1).ne.(nx+1)).or.
     >     (vardim(2).ne.(ny+1)).or.
     >     (vardim(3).ne.(nz+1)).or.
     >     (abs(varmin(1)-xmin).gt.eps).or.
     >     (abs(varmin(2)-ymin).gt.eps).or.
     >     (abs(varmin(3)-zmin).gt.eps).or.
     >     (abs((varmax(1)-varmin(1))/real(vardim(1)-1)-dx).gt.eps).or.
     >     (abs((varmax(2)-varmin(2))/real(vardim(2)-1)-dy).gt.eps).or.
     >     (abs((varmax(3)-varmin(3))/real(vardim(3)-1)-dz).gt.eps) ) 
     >then
         print*,'Input grid inconsitency...'
         print*,'  Nx      = ',vardim(1),nx+1
         print*,'  Ny      = ',vardim(2),ny+1
         print*,'  Nz      = ',vardim(3),nz+1
         print*,'  Varminx = ',varmin(1),xmin
         print*,'  Varminy = ',varmin(2),ymin
         print*,'  Varminz = ',varmin(3),zmin
         print*,'  Dx      = ',(varmax(1)-varmin(1))/real(nx-1),dx
         print*,'  Dy      = ',(varmax(2)-varmin(2))/real(ny-1),dy
         print*,'  Dz      = ',(varmax(3)-varmin(3))/real(nz-1),dz
         goto 998
      endif

c     Load variables
      call getdef(cdfid,varname,ndimin,misdat,vardim,
     >            varmin,varmax,stag,stat)
      if (stat.ne.0) goto 998
      call getdat(cdfid,varname,time(1),0,field,stat)
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
      if (stat.ne.0) goto 998

c     Close input netcdf file
      call clscdf(cdfid,stat)
      if (stat.ne.0) goto 998

c     Set missing data value to <mdv>
      do i=1,nx
         do j=1,ny
            do k=1,nz
               if (abs(field(i,j,k)-misdat).lt.eps) then
                  field(i,j,k)=mdv
               endif
            enddo
         enddo
      enddo

      return

c     Exception handling
 998  print*,'Read_Inp: Problem with input netcdf file... Stop'
      stop

      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)>. If this is 
C     the case, <isok> is incremented by 1. Otherwise <isok> 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

c     --------------------------------------------------------------------------------
c     Get grid parameters
c     --------------------------------------------------------------------------------

      subroutine read_dim (nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv,
     >                     pvsrcfile)

c     Get the grid parameters from the variable <TH> on the input file <pvsrcfile>.
c     The grid parameters are
c        nx,ny,nz                : Number of grid points in x, y and z direction
c        xmin,ymin,zmin          : Minimum domain coordinates in x, y and z direction
c        xmax,ymax,zmax          : Maximal domain coordinates in x, y and z direction
c        dx,dy,dz                : Horizontal and vertical resolution
c     Additionally, it is checked whether the vertical grid is equally spaced. If ok,
c     the grid paramters are transformed from lon/lat to distance (in meters)

      implicit none

c     Declaration of subroutine parameters
      character*80   pvsrcfile
      integer        nx,ny,nz
      real           dx,dy,dz
      real           xmin,ymin,zmin,xmax,ymax,zmax
      real           mdv

c     Numerical epsilon and other physical/geoemtrical parameters
      real           eps
      parameter      (eps=0.01)

c     Auxiliary variables
      integer        cdfid,cstid
      integer        ierr
      character*80   vnam(100),varname
      integer        nvars
      integer        isok
      integer        vardim(4)
      real           misdat
      real           varmin(4),varmax(4),stag(4)
      real           aklev(1000),bklev(1000),aklay(1000),bklay(1000)
      real           dh
      character*80   csn
      integer        ndim
      integer        i

c     Get all grid parameters
      call cdfopn(pvsrcfile,cdfid,ierr)
      if (ierr.ne.0) goto 998
      call getvars(cdfid,nvars,vnam,ierr)
      if (ierr.ne.0) goto 998
      isok=0
      varname='T'
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) goto 998
      call getcfn(cdfid,csn,ierr)
      if (ierr.ne.0) goto 998
      call cdfopn(csn,cstid,ierr)
      if (ierr.ne.0) goto 998
      call getdef(cdfid,varname,ndim,misdat,vardim,varmin,varmax,
     >            stag,ierr)
      if (ierr.ne.0) goto 998
      nx=vardim(1)
      ny=vardim(2)
      nz=vardim(3)
      xmin=varmin(1)
      ymin=varmin(2)
      zmin=varmin(3)
      call getlevs(cstid,nz,aklev,bklev,aklay,bklay,ierr)
      if (ierr.ne.0) goto 998
      call getgrid(cstid,dx,dy,ierr)
      if (ierr.ne.0) goto 998
      xmax=varmax(1)
      ymax=varmax(2)
      zmax=varmax(3)
      dz=(zmax-zmin)/real(nz-1)
      call clscdf(cstid,ierr)
      if (ierr.ne.0) goto 998
      call clscdf(cdfid,ierr)
      if (ierr.ne.0) goto 998

c     Check whether the grid is equallay spaced in the vertical
      do i=1,nz-1
         dh=aklev(i+1)-aklev(i)
         if (abs(dh-dz).gt.eps) then
            print*,'Aklev: Vertical grid must be equally spaced... Stop'
            print*,(aklev(i),i=1,nz)
            stop
         endif
         dh=aklay(i+1)-aklay(i)
         if (abs(dh-dz).gt.eps) then
            print*,'Aklay: Vertical grid must be equally spaced... Stop'
            print*,(aklay(i),i=1,nz)
            stop
         endif
      enddo

c     Set missing data value
      mdv=misdat

      return

c     Exception handling
 998  print*,'Read_Dim: Problem with input netcdf file... Stop'
      stop

      end


c     ********************************************************************************
c     * DEFINE REFERENCE PROFILE                                                     *
c     ********************************************************************************

c     --------------------------------------------------------------------------------
c     Calculate area mean
c     --------------------------------------------------------------------------------

      SUBROUTINE mean(a,m,nx,ny,nz,mdv,weight)

c     Calculate the area-mean of <m> and save it on <a>.

      implicit none

c     Declaration of subroutine parameters
      real       mdv
      integer    nx,ny,nz
      real       m(0:nx,0:ny,0:nz),a(0:2*nz)
      real       weight(0:nx,0:ny)

c     Numerical epsilon
      real       eps
      parameter  (eps=0.01)

c     Auxiliary varaibles
      real       mea(0:nz)
      integer    i,j,k,kk
      real       counter

c     Determine the mean over all levels (handle missing data)
      do k=0,nz
         mea(k)=0.
         counter=0.
         do i=0,nx
            do j=0,ny
               if (abs(m(i,j,k)-mdv).gt.eps) then
                  mea(k)=mea(k)+m(i,j,k)*weight(i,j)
                  counter=counter+weight(i,j)
               endif
               
            enddo
         enddo
         if (counter.gt.0) then
            mea(k)=mea(k)/counter
         else
            mea(k)=mdv
         endif
      enddo

c     Prepare the output array: split layers and levels
      do k=0,nz-1
         kk=2*k
         a(kk)=mea(k)
         a(kk+1)=0.5*(mea(k)+mea(k+1))
      enddo
      a(2*nz)=mea(nz)
      
      end