Subversion Repositories pvinversion.ecmwf

Rev

Go to most recent revision | Blame | Last modification | View Log | Download | RSS feed 

      PROGRAM prep_iteration

c     ************************************************************************
c     * Prepare the next step for the PV inversion                           *
c     * Michael Sprenger / Summer, Autumn 2006                               *
c     ************************************************************************

      implicit none

c     ------------------------------------------------------------------------
c     Declaration of variables and parameters
c     ------------------------------------------------------------------------
      
c     Input and output file
      character*80   anomafile
      character*80   iterafile

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

c     Numerical epsilon and other variables
      real           eps
      parameter      (eps=0.01)
      real           alpha

c     3d arrays
      real,allocatable,dimension (:,:,:) :: v_iter,v_anom
      real,allocatable,dimension (:,:,:) :: u_iter,u_anom
      real,allocatable,dimension (:,:,:) :: t_iter,t_anom
      real,allocatable,dimension (:,:,:) :: p_iter,p_anom

c     Auciliary variables 
      integer      i,j,k
      integer      stat
      character*80 varname
      character*80 name

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

      print*,'********************************************************'
      print*,'* PREP_ITERATION                                       *'
      print*,'********************************************************'

c     Read parameter file
      open(10,file='fort.10')
       read(10,*) iterafile
       read(10,*) anomafile
       read(10,*) name,alpha
       if ( trim(name).ne.'ALPHA' ) stop
      close(10) 
      print*
      print*,trim(anomafile)
      print*,trim(iterafile)
      print*

c     Get lat/lon gid parameters from input file
      call read_dim (nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv,
     >               anomafile)
      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*

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 3d arrays 
      allocate(u_anom (0:nx,0:ny,0:nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating u_anom'
      allocate(v_anom (0:nx,0:ny,0:nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating v_anom'
      allocate(t_anom  (0:nx,0:ny,0:nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating t_anom'
      allocate(p_anom  (0:nx,0:ny,0:nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating p_anom'
      
      allocate(u_iter (0:nx,0:ny,0:nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating u_iter'
      allocate(v_iter (0:nx,0:ny,0:nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating v_iter'
      allocate(t_iter  (0:nx,0:ny,0:nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating t_iter'
      allocate(p_iter  (0:nx,0:ny,0:nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating p_iter'
  
c     Read anomaly and iteration fields
      varname='U'
      call read_inp (u_anom,varname,anomafile,
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
      varname='V'
      call read_inp (v_anom,varname,anomafile,
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
      varname='T'
      call read_inp (t_anom,varname,anomafile,
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
      varname='P'
      call read_inp (p_anom,varname,anomafile,
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
      varname='U'
      call read_inp (u_iter,varname,iterafile,
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
      varname='V'
      call read_inp (v_iter,varname,iterafile,
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
      varname='T'
      call read_inp (t_iter,varname,iterafile,
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
      varname='P'
      call read_inp (p_iter,varname,iterafile,
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)

c     --------------------------------------------------------------------------------
c     Modify the iteration fields
c     --------------------------------------------------------------------------------

      do i=0,nx
         do j=0,ny
            do k=0,nz

c              Update zonal velocity
               if ((abs(u_anom(i,j,k)-mdv).gt.eps).and.
     >             (abs(u_iter(i,j,k)-mdv).gt.eps)) then
                  u_iter(i,j,k)=u_iter(i,j,k)-alpha*u_anom(i,j,k)
               else
                  u_iter(i,j,k)=mdv
               endif

c              Update meridional velocity
               if ((abs(v_anom(i,j,k)-mdv).gt.eps).and.
     >             (abs(v_iter(i,j,k)-mdv).gt.eps)) then
                  v_iter(i,j,k)=v_iter(i,j,k)-alpha*v_anom(i,j,k)
               else
                  v_iter(i,j,k)=mdv
               endif

c              Update temperature
               if ((abs(t_anom(i,j,k)-mdv).gt.eps).and.
     >             (abs(t_iter(i,j,k)-mdv).gt.eps)) then
                  t_iter(i,j,k)=t_iter(i,j,k)-alpha*t_anom(i,j,k)
               else
                  t_iter(i,j,k)=mdv
               endif

c              Update pressure
               if ((abs(p_anom(i,j,k)-mdv).gt.eps).and.
     >             (abs(p_iter(i,j,k)-mdv).gt.eps)) then
                  p_iter(i,j,k)=p_iter(i,j,k)-alpha*p_anom(i,j,k)
               else
                  p_iter(i,j,k)=mdv
               endif

            enddo
         enddo
      enddo


c     --------------------------------------------------------------------------------
c     Write output
c     --------------------------------------------------------------------------------

      varname='U'
      call write_inp (u_iter,varname,iterafile,nx,ny,nz)
      varname='V'
      call write_inp (v_iter,varname,iterafile,nx,ny,nz)
      varname='T'
      call write_inp (t_iter,varname,iterafile,nx,ny,nz)
      varname='P'
      call write_inp (p_iter,varname,iterafile,nx,ny,nz)


      end




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

c     --------------------------------------------------------------------------------
c     Write input field to netcdf
c     --------------------------------------------------------------------------------

      SUBROUTINE write_inp (field,fieldname,pvsrcfile,nx,ny,nz)

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. 

      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

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

c     Get grid parameters
      call cdfopn(pvsrcfile,cdfid,stat)
      if (stat.ne.0) goto 998
      call getvars(cdfid,nvars,vnam,stat)
      if (stat.ne.0) goto 998
      isok=0
      varname='TH'
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) goto 998
      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
      call clscdf(cdfid,stat)
      if (stat.ne.0) goto 998

c     Save variables (write definition, if necessary)
      call cdfwopn(pvsrcfile,cdfid,stat)
      if (stat.ne.0) goto 998
      isok=0
      varname=fieldname
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) then
         call putdef(cdfid,varname,ndimin,misdat,vardim,
     >               varmin,varmax,stag,stat)
         if (stat.ne.0) goto 998
      endif
      call putdat(cdfid,varname,time(1),0,field,stat)
      print*,'W ',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
      
      return

c     Exception handling
 998  print*,'Write_Inp: 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 from theta     
      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 <THETA> 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='TH'
      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