WebSVN – lagranto.ecmwf – Diff – /trunk/goodies/getmima.f


      PROGRAM getmima
 
c     ***********************************************************************
c     * Get minimum and maximum value of a varaiable on pressure            *
c     * or isentropic surface                                               *
c     * Michael Sprenger / Spring, summer 2016                              *
c     ***********************************************************************
 
      use netcdf
 
      implicit none
 
c     ----------------------------------------------------------------------
c     Declaration of parameters and variables
c     ----------------------------------------------------------------------
      
c     Interpolation method
      integer   ipom
      parameter (ipom=0)    
      
c     Flag for timecheck
      character*80 timecheck
      parameter    (timecheck = 'no' )        
      
c     netCDF fields      
      real,allocatable, dimension (:,:)   :: sp,varf
      real,allocatable, dimension (:,:,:) :: field,varl,tt,p3
      real,allocatable, dimension (:)     :: ak,bk
      integer                                stat
      character*80                           cdfnam,varnam
      character*80                           vnam(200)
      integer                                nvars
      integer	                             cdfid,ierr
      real	                                 mdv
      real                                   stagz
      
c     Grid description      
      real      xmin,xmax       ! Zonal grid extension
      real      ymin,ymax       ! Meridional grid extension
      integer   nx,ny,nz        ! Grid dimensions
      real      dx,dy           ! Horizontal grid resolution
      real	    pollon,pollat   ! Pole position
 
c     Output variables
      real      min,max                      ! Minimum & maximum value
      real      lonmin,lonmax,latmin,latmax  ! Position of min & max
 
c     Auxiliary variables      
      integer        iargc
      character*(80) arg
      real           rd
      integer        i,j
      integer        minindx,minindy,maxindx,maxindy
      character*80   tvar
      character*1    mode
      character*80   clev
      real	         xphys,yphys
      real	         level
      integer        flag
       
c     Externals      
      real      lmstolm,phstoph
      external  lmstolm,phstoph
 
c     ----------------------------------------------------------------------
c     Preparation - argument handling, grid paraemters, memor allocation
c     ----------------------------------------------------------------------
 
c     Check for sufficient requested arguments
      if (iargc().lt.2) then
        print*,'USAGE: getmima filename var ',
     >         '[lev in the form Pval or Tval]'
        call exit(1)
      endif
 
c     Read and transform input
      call getarg(1,arg)
      cdfnam=trim(arg)
      call getarg(2,arg)
      varnam=trim(arg)
 
      if (iargc().eq.3) then
        call getarg(3,arg)
        mode=arg(1:1)
        clev=arg(2:len(arg)) 
        call checkchar(clev,".",flag)
        if (flag.eq.0) clev=trim(clev)//"."
        read(clev,'(f10.2)') level
      else
        mode='X'
        level=0.
      endif
 
c     Init level type and minimum,maximum
      min= 1.e19
      max=-1.e19
 
c     Open netCDF file      
      call input_open (cdfid,cdfnam)
 
c     Get list of variables on file
      call input_getvars (cdfid,vnam,nvars)
 
C     Get infos about data domain - in particular: dimensions
      nx       = 1
      ny       = 1
      nz       = 1
      call input_grid (-cdfid,varnam,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
     >               0.,pollon,pollat,rd,rd,nz,rd,rd,stagz,timecheck)
 
C     Get memory for dynamic arrays
      allocate(sp(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array sp ***'
      allocate(varf(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array varf ***'
      allocate(varl(nx,ny,nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array varl ***'
      allocate(tt(nx,ny,nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tt ***'
      allocate(field(nx,ny,nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array field ***'
      allocate(ak(nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array ak ***'
      allocate(bk(nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array bk ***'
   
c     ----------------------------------------------------------------------
c     Load data
c     ----------------------------------------------------------------------
 
c     Get grid info for variable - non-constant part  (P, PS, AK, BK)   
      call input_grid                               
     >       (cdfid,varnam,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
     >        0.,pollon,pollat,varl,sp,nz,ak,bk,stagz,timecheck)
    
c     Load Variable    
      call input_wind (cdfid,varnam,field,0.,stagz,mdv,
     >                 xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
 
c     In case of isentropic level, read temperature
      if (mode.eq.'T') then
            tvar="T"
            do i=1,nvars
              if (vnam(i).eq."THETA") tvar="THETA"
              if (vnam(i).eq."TH")    tvar="TH"
            enddo
            if (tvar.eq.'T') then
                 call input_wind (cdfid,tvar,tt,0.,stagz,mdv,xmin,xmax,
     >                            ymin,ymax,dx,dy,nx,ny,nz,timecheck)
                 call pottemp(varl,tt,sp,nx,ny,nz,ak,bk)
            else
                 call input_wind (cdfid,tvar,tt,0.,stagz,mdv,xmin,xmax,
     >                            ymin,ymax,dx,dy,nx,ny,nz,timecheck)
            endif
      endif
 
c     ----------------------------------------------------------------------
c     Do interpolation and get minimum, maximum
c     ----------------------------------------------------------------------
 
C     If required do interpolation on pressure or theta level
      if (mode.eq.'P') then
        call vipo(field,varl,level,varf,nx,ny,nz,mdv,ipom)
      elseif (mode.eq.'T') then
        call thipo(field,varl,level,varf,nx,ny,nz,mdv,ipom)
      elseif (mode.eq.'L') then
        do i=1,nx
           do j=1,ny
               varf(i,j) = field(i,j,nint(level))
           enddo
        enddo
      elseif (mode.eq.'X') then
         do i=1,nx
            do j=1,ny
                varf(i,j) = field(i,j,1)
            enddo
         enddo
      endif
 
C     Determine minimum and maximum
      do i=1,nx
        do j=1,ny
 
            if ((varf(i,j).ne.mdv).and.(varf(i,j).lt.min)) then
               min     = varf(i,j)
               minindx = i
               minindy = j
            endif
            if ((varf(i,j).ne.mdv).and.(varf(i,j).gt.max)) then
               max     = varf(i,j)
               maxindx = i
               maxindy = j
          endif
        enddo
      enddo
      lonmin = xmin + real(minindx-1) * dx
      latmin = ymin + real(minindy-1) * dy
      lonmax = xmin + real(maxindx-1) * dx
      latmax = ymin + real(maxindy-1) * dy
 
c     Rotate position to true longitude/latitude 
      if ((pollon.ne.0.).or.(pollat.ne.90.)) then
        xphys=lmstolm(latmin,lonmin,pollat,pollon)
        yphys=phstoph(latmin,lonmin,pollat,pollon)
        lonmin=xphys
        latmin=yphys
        xphys=lmstolm(latmax,lonmax,pollat,pollon)
        yphys=phstoph(latmax,lonmax,pollat,pollon)
        lonmax=xphys
        latmax=yphys
      endif
 
c     Write output
      write(*,101) min,max,lonmin,latmin,nint(level),
     >                     lonmax,latmax,nint(level)
  101 format(2f10.3,2f8.2,i6,2f8.2,i6)
 
c     Close netCDF file
      call input_close(cdfid)
 
      end
      
c     ----------------------------------------------------------------------
C     Interpolates the 3d variable var3d on the pressure surface defined
C     by lev of the variable varl. The interpolated field is
C     returned as var.
C     ipom determines the way of vertical interpolation:
C       ipom=0 is for linear interpolation
C       ipom=1 is for logarithmic interpolation
c     ----------------------------------------------------------------------
 
      subroutine vipo(var3d,varl,lev,var,nx,ny,nz,mdv,ipom)
 
      integer   nx,ny,nz,ipom
      real      lev,mdv
      real      var3d(nx,ny,nz),varl(nx,ny,nz),var(nx,ny)
 
      integer   i,j,k
      real      kind
      real      int3dm
 
      do i=1,nx
      do j=1,ny
 
        do k=1,nz-1
          if ((varl(i,j,k).ge.lev).and.(varl(i,j,k+1).le.lev)) then
            kind=float(k)+(varl(i,j,k)-lev)/
     >                   (varl(i,j,k)-varl(i,j,k+1))
            goto 100
          endif
        enddo
 100    continue
 
        var(i,j)=int3dm(var3d,nx,ny,nz,float(i),float(j),kind,mdv)
 
      enddo
      enddo
 
      return
      end
      
c     ----------------------------------------------------------------------
C     Interpolates the 3d variable var3d on the isentropic surface defined
C     by lev of the variable varl. The interpolated field is
C     returned as var.
C     ipom determines the way of vertical interpolation:
C       ipom=0 is for linear interpolation
C       ipom=1 is for logarithmic interpolation
c     ----------------------------------------------------------------------
 
      subroutine thipo(var3d,th3d,lev,var,nx,ny,nz,mdv,mode)
 
      integer   nx,ny,nz,mode
      real      lev,mdv
      real      var3d(nx,ny,nz),th3d(nx,ny,nz),var(nx,ny)
 
      integer   i,j,k
      real      kind
      real      int3dm
 
      do i=1,nx
      do j=1,ny
 
        do k=1,nz-1
          if ((th3d(i,j,k).le.lev).and.(th3d(i,j,k+1).ge.lev)) then
            kind=float(k)+(th3d(i,j,k)-lev)/
     >                   (th3d(i,j,k)-th3d(i,j,k+1))
            goto 100
          endif
        enddo
 100    continue
 
        var(i,j)=int3dm(var3d,nx,ny,nz,float(i),float(j),kind,mdv)
 
      enddo
      enddo
 
      return
      end
      
c     ----------------------------------------------------------------------
c     Check whether character appears within string
c     ----------------------------------------------------------------------
 
      subroutine checkchar(string,char,flag)
 
      character*(*)     string
      character*(1)     char
      integer           n,flag
 
      flag=0
      do n=1,len(string)
        if (string(n:n).eq.char) then
          flag=n
          return
        endif
      enddo
      end
      
c     ----------------------------------------------------------------------
c     3D interpolation
c     ----------------------------------------------------------------------
      
      real function int3d(ar,n1,n2,n3,rid,rjd,rkd)
 
c     Purpose:
c        This subroutine interpolates a 3d-array to an arbitrary
c        location within the grid.
c     Arguments:
c        ar        real  input   surface pressure, define as ar(n1,n2,n3)
c        n1,n2,n3  int   input   dimensions of ar
c        ri,rj,rk  real  input   grid location to be interpolated to
c     History:
 
c     argument declarations
      integer   n1,n2,n3
      real      ar(n1,n2,n3), rid,rjd,rkd
 
c     local declarations
      integer   i,j,k,ip1,jp1,kp1,ih,jh,kh
      real      frac0i,frac0j,frac0k,frac1i,frac1j,frac1k,ri,rj,rk
 
c     do linear interpolation
      ri=amax1(1.,amin1(float(n1),rid))
      rj=amax1(1.,amin1(float(n2),rjd))
      rk=amax1(1.,amin1(float(n3),rkd))
      ih=nint(ri)
      jh=nint(rj)
      kh=nint(rk)
 
c     Check for interpolation in i
*     if (abs(float(ih)-ri).lt.1.e-3) then
*       i  =ih
*       ip1=ih
*     else
        i =min0(int(ri),n1-1)
        ip1=i+1
*     endif
 
c     Check for interpolation in j
      if (abs(float(jh)-rj).lt.1.e-3) then
        j  =jh
        jp1=jh
      else
        j =min0(int(rj),n2-1)
        jp1=j+1
      endif
 
c     Check for interpolation in k
*     if (abs(float(kh)-rk).lt.1.e-3) then
*       k  =kh
*       kp1=kh
*     else
        k =min0(int(rk),n3-1)
        kp1=k+1
*     endif
 
      if (k.eq.kp1) then
c       no interpolation in k
        if ((i.eq.ip1).and.(j.eq.jp1)) then
c          no interpolation at all
           int3d=ar(i,j,k)
c          print *,'int3d 00: ',rid,rjd,rkd,int3d
        else
c          horizontal interpolation only
           frac0i=ri-float(i)
           frac0j=rj-float(j)
           frac1i=1.-frac0i
           frac1j=1.-frac0j
           int3d = ar(i  ,j  ,k  ) * frac1i * frac1j
     &           + ar(i  ,jp1,k  ) * frac1i * frac0j
     &           + ar(ip1,j  ,k  ) * frac0i * frac1j
     &           + ar(ip1,jp1,k  ) * frac0i * frac0j
c          print *,'int3d 10: ',rid,rjd,rkd,int3d
        endif
      else 
        frac0k=rk-float(k)
        frac1k=1.-frac0k
        if ((i.eq.ip1).and.(j.eq.jp1)) then
c          vertical interpolation only
           int3d = ar(i  ,j  ,k  ) * frac1k
     &           + ar(i  ,j  ,kp1) * frac0k
c          print *,'int3d 01: ',rid,rjd,rkd,int3d
        else
c          full 3d interpolation
           frac0i=ri-float(i)
           frac0j=rj-float(j)
           frac1i=1.-frac0i
           frac1j=1.-frac0j
           int3d = ar(i  ,j  ,k  ) * frac1i * frac1j * frac1k
     &           + ar(i  ,jp1,k  ) * frac1i * frac0j * frac1k 
     &           + ar(ip1,j  ,k  ) * frac0i * frac1j * frac1k
     &           + ar(ip1,jp1,k  ) * frac0i * frac0j * frac1k
     &           + ar(i  ,j  ,kp1) * frac1i * frac1j * frac0k
     &           + ar(i  ,jp1,kp1) * frac1i * frac0j * frac0k 
     &           + ar(ip1,j  ,kp1) * frac0i * frac1j * frac0k
     &           + ar(ip1,jp1,kp1) * frac0i * frac0j * frac0k
c          print *,'int3d 11: ',rid,rjd,rkd,int3d
        endif
      endif
      end
      
c     ----------------------------------------------------------------------
c     3D interpolation including missing data check
c     ----------------------------------------------------------------------
  
      real function int3dm(ar,n1,n2,n3,rid,rjd,rkd,misdat)
 
c     Purpose:
c        This subroutine interpolates a 3d-array to an arbitrary
c        location within the grid. The interpolation includes the 
c        testing of the missing data flag 'misdat'. 
c     Arguments:
c        ar        real  input   surface pressure, define as ar(n1,n2,n3)
c        n1,n2,n3  int   input   dimensions of ar
c        ri,rj,rk  real  input   grid location to be interpolated to
c        misdat    real  input   missing data flag (on if misdat<>0)
c     Warning:
c        This routine has not yet been seriously tested
c     History:
 
c     argument declarations
      integer   n1,n2,n3
      real      ar(n1,n2,n3), rid,rjd,rkd, misdat
 
c     local declarations
      integer   i,j,k,ip1,jp1,kp1,ih,jh,kh
      real      frac0i,frac0j,frac0k,frac1i,frac1j,frac1k,ri,rj,rk,int3d
 
c     check if routine without missing data checking can be called instead
      if (misdat.eq.0.) then
        int3dm=int3d(ar,n1,n2,n3,rid,rjd,rkd)
        return
      endif
 
c     do linear interpolation
      ri=amax1(1.,amin1(float(n1),rid))
      rj=amax1(1.,amin1(float(n2),rjd))
      rk=amax1(1.,amin1(float(n3),rkd))
      ih=nint(ri)
      jh=nint(rj)
      kh=nint(rk)
 
c     Check for interpolation in i
*     if (abs(float(ih)-ri).lt.1.e-3) then
*       i  =ih
*       ip1=ih
*     else
        i =min0(int(ri),n1-1)
        ip1=i+1
*     endif
 
c     Check for interpolation in j
*     if (abs(float(jh)-rj).lt.1.e-3) then
*       j  =jh
*       jp1=jh
*     else
        j =min0(int(rj),n2-1)
        jp1=j+1
*     endif
 
c     Check for interpolation in k
*     if (abs(float(kh)-rk).lt.1.e-3) then
*       k  =kh
*       kp1=kh
*     else
        k =min0(int(rk),n3-1)
        kp1=k+1
*     endif
 
      if (k.eq.kp1) then
c       no interpolation in k
        if ((i.eq.ip1).and.(j.eq.jp1)) then
c          no interpolation at all
           if (misdat.eq.ar(i,j,k)) then
             int3dm=misdat
           else
             int3dm=ar(i,j,k)
           endif
c          print *,'int3dm 00: ',rid,rjd,rkd,int3dm
        else
c          horizontal interpolation only
           if ((misdat.eq.ar(i  ,j  ,k  )).or.
     &         (misdat.eq.ar(i  ,jp1,k  )).or.
     &         (misdat.eq.ar(ip1,j  ,k  )).or.
     &         (misdat.eq.ar(ip1,jp1,k  ))) then
             int3dm=misdat
           else
             frac0i=ri-float(i)
             frac0j=rj-float(j)
             frac1i=1.-frac0i
             frac1j=1.-frac0j
             int3dm = ar(i  ,j  ,k  ) * frac1i * frac1j
     &              + ar(i  ,jp1,k  ) * frac1i * frac0j
     &              + ar(ip1,j  ,k  ) * frac0i * frac1j
     &              + ar(ip1,jp1,k  ) * frac0i * frac0j
c            print *,'int3dm 10: ',rid,rjd,rkd,int3dm
           endif
        endif
      else 
        frac0k=rk-float(k)
        frac1k=1.-frac0k
        if ((i.eq.ip1).and.(j.eq.jp1)) then
c          vertical interpolation only
           if ((misdat.eq.ar(i  ,j  ,k  )).or.
     &         (misdat.eq.ar(i  ,j  ,kp1))) then
             int3dm=misdat
           else
             int3dm = ar(i  ,j  ,k  ) * frac1k
     &              + ar(i  ,j  ,kp1) * frac0k
c            print *,'int3dm 01: ',rid,rjd,rkd,int3dm
           endif
        else
c          full 3d interpolation
           if ((misdat.eq.ar(i  ,j  ,k  )).or.
     &         (misdat.eq.ar(i  ,jp1,k  )).or.
     &         (misdat.eq.ar(ip1,j  ,k  )).or.
     &         (misdat.eq.ar(ip1,jp1,k  )).or.
     &         (misdat.eq.ar(i  ,j  ,kp1)).or.
     &         (misdat.eq.ar(i  ,jp1,kp1)).or.
     &         (misdat.eq.ar(ip1,j  ,kp1)).or.
     &         (misdat.eq.ar(ip1,jp1,kp1))) then
             int3dm=misdat
           else
             frac0i=ri-float(i)
             frac0j=rj-float(j)
             frac1i=1.-frac0i
             frac1j=1.-frac0j
             int3dm = ar(i  ,j  ,k  ) * frac1i * frac1j * frac1k
     &              + ar(i  ,jp1,k  ) * frac1i * frac0j * frac1k
     &              + ar(ip1,j  ,k  ) * frac0i * frac1j * frac1k
     &              + ar(ip1,jp1,k  ) * frac0i * frac0j * frac1k
     &              + ar(i  ,j  ,kp1) * frac1i * frac1j * frac0k
     &              + ar(i  ,jp1,kp1) * frac1i * frac0j * frac0k
     &              + ar(ip1,j  ,kp1) * frac0i * frac1j * frac0k
     &              + ar(ip1,jp1,kp1) * frac0i * frac0j * frac0k
c            print *,'int3dm 11: ',rid,rjd,rkd,int3dm
           endif
        endif
      endif
      end
 
c     ----------------------------------------------------------------------
c     Calculate potential temperature
c     ----------------------------------------------------------------------
 
      subroutine pottemp(pt,t,sp,ie,je,ke,ak,bk)
 
c     argument declaration
      integer   ie,je,ke
      real      pt(ie,je,ke),t(ie,je,ke),sp(ie,je),
     >     	ak(ke),bk(ke)
 
c     variable declaration
      integer   i,j,k
      real      rdcp,tzero
      data      rdcp,tzero /0.286,273.15/
 
c     statement-functions for the computation of pressure
      real      pp,psrf
      integer   is
      pp(is)=ak(is)+bk(is)*psrf
 
c     computation of potential temperature
      do i=1,ie
      do j=1,je
        psrf=sp(i,j)
        do k=1,ke
c     distinction of temperature in K and deg C
          if (t(i,j,k).lt.100.) then
            pt(i,j,k)=(t(i,j,k)+tzero)*( (1000./pp(k))**rdcp )
          else
            pt(i,j,k)=t(i,j,k)*( (1000./pp(k))**rdcp )
          endif
        enddo
      enddo
      enddo
      end
      
c     ----------------------------------------------------------------------
c     Calculate 3D pressure
c     ----------------------------------------------------------------------
 
      subroutine pres(pr,sp,ie,je,ke,ak,bk)
 
c     argument declaration
      integer  ie,je,ke
      real,intent(OUT) :: pr(ie,je,ke)
      real,intent(IN)  :: sp(ie,je)
      real,intent(IN)  :: ak(ke),bk(ke)
 
c     variable declaration
      integer  i,j,k
 
c     computation pressure
      do i=1,ie
      do j=1,je
      do k=1,ke
        pr(i,j,k)=ak(k)+bk(k)*sp(i,j)
      enddo
      enddo
      enddo
      end
      
c     ----------------------------------------------------------------------
c     Coordinate rotations from COSMO
c     ----------------------------------------------------------------------      
      
      REAL FUNCTION PHTOPHS (PHI, LAM, POLPHI, POLLAM)
C
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
C
C**** PHTOPHS  -   FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI
C****                 AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   AUFRUF   :   PHI = PHTOPHS (PHI, LAM, POLPHI, POLLAM)
C**   ENTRIES  :   KEINE
C**   ZWECK    :   UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI AUF
C**                EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   VERSIONS-
C**   DATUM    :   03.05.90
C**
C**   EXTERNALS:   KEINE
C**   EINGABE-
C**   PARAMETER:   PHI    REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
C**                LAM    REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
C**                POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
C**                POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
C**   AUSGABE-
C**   PARAMETER:   ROTIERTE BREITE PHIS ALS WERT DER FUNKTION
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
C**
C**   COMMON-
C**   BLOECKE  :   KEINE
C**
C**   FEHLERBE-
C**   HANDLUNG :   KEINE
C**   VERFASSER:   G. DE MORSIER
 
      REAL        LAM,PHI,POLPHI,POLLAM
 
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
 
      ZSINPOL = SIN(ZPIR18*POLPHI)
      ZCOSPOL = COS(ZPIR18*POLPHI)
      ZLAMPOL = ZPIR18*POLLAM
      ZPHI    = ZPIR18*PHI
      ZLAM    = LAM
      IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
      ZLAM    = ZPIR18*ZLAM
      ZARG    = ZCOSPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL) + ZSINPOL*SIN(ZPHI)
 
      PHTOPHS = ZRPI18*ASIN(ZARG)
 
      RETURN
      END
      
      
      REAL FUNCTION PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)
C
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
C
C**** PHSTOPH  -   FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE FUER
C****                 EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
C****                 ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   AUFRUF   :   PHI = PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)
C**   ENTRIES  :   KEINE
C**   ZWECK    :   BERECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE FUER
C**                EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   VERSIONS-
C**   DATUM    :   03.05.90
C**
C**   EXTERNALS:   KEINE
C**   EINGABE-
C**   PARAMETER:   PHIS     REAL   GEOGR. BREITE DES PUNKTES IM ROT.SYS.
C**                LAMS     REAL   GEOGR. LAENGE DES PUNKTES IM ROT.SYS.
C**                POLPHI   REAL   WAHRE GEOGR. BREITE DES NORDPOLS
C**                POLLAM   REAL   WAHRE GEOGR. LAENGE DES NORDPOLS
C**   AUSGABE-
C**   PARAMETER:   WAHRE GEOGRAPHISCHE BREITE ALS WERT DER FUNKTION
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
C**
C**   COMMON-
C**   BLOECKE  :   KEINE
C**
C**   FEHLERBE-
C**   HANDLUNG :   KEINE
C**   VERFASSER:   D.MAJEWSKI
 
      REAL        LAMS,PHIS,POLPHI,POLLAM
 
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
 
      SINPOL = SIN(ZPIR18*POLPHI)
      COSPOL = COS(ZPIR18*POLPHI)
      ZPHIS  = ZPIR18*PHIS
      ZLAMS  = LAMS
      IF(ZLAMS.GT.180.0) ZLAMS = ZLAMS - 360.0
      ZLAMS  = ZPIR18*ZLAMS
      ARG     = COSPOL*COS(ZPHIS)*COS(ZLAMS) + SINPOL*SIN(ZPHIS)
 
      PHSTOPH = ZRPI18*ASIN(ARG)
 
      RETURN
      END
      REAL FUNCTION LMTOLMS (PHI, LAM, POLPHI, POLLAM)
C
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
C
C**** LMTOLMS  -   FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM
C****                 AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   AUFRUF   :   LAM = LMTOLMS (PHI, LAM, POLPHI, POLLAM)
C**   ENTRIES  :   KEINE
C**   ZWECK    :   UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM AUF
C**                EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   VERSIONS-
C**   DATUM    :   03.05.90
C**
C**   EXTERNALS:   KEINE
C**   EINGABE-
C**   PARAMETER:   PHI    REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
C**                LAM    REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
C**                POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
C**                POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
C**   AUSGABE-
C**   PARAMETER:   WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
C**
C**   COMMON-
C**   BLOECKE  :   KEINE
C**
C**   FEHLERBE-
C**   HANDLUNG :   KEINE
C**   VERFASSER:   G. DE MORSIER
 
      REAL        LAM,PHI,POLPHI,POLLAM
 
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
 
      ZSINPOL = SIN(ZPIR18*POLPHI)
      ZCOSPOL = COS(ZPIR18*POLPHI)
      ZLAMPOL =     ZPIR18*POLLAM
      ZPHI    =     ZPIR18*PHI
      ZLAM    = LAM
      IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
      ZLAM    = ZPIR18*ZLAM
 
      ZARG1   = - SIN(ZLAM-ZLAMPOL)*COS(ZPHI)
      ZARG2   = - ZSINPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL)+ZCOSPOL*SIN(ZPHI)
      IF (ABS(ZARG2).LT.1.E-30) THEN
        IF (ABS(ZARG1).LT.1.E-30) THEN
          LMTOLMS =   0.0
        ELSEIF (ZARG1.GT.0.) THEN
              LMTOLMS =  90.0
            ELSE
              LMTOLMS = -90.0
            ENDIF
      ELSE
        LMTOLMS = ZRPI18*ATAN2(ZARG1,ZARG2)
      ENDIF
 
      RETURN
      END
      
      
      REAL FUNCTION LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)
C
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
C
C**** LMSTOLM  -   FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE FUER
C****                 EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   AUFRUF   :   LAM = LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)
C**   ENTRIES  :   KEINE
C**   ZWECK    :   BERECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE FUER
C**                EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
C**                IM ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   VERSIONS-
C**   DATUM    :   03.05.90
C**
C**   EXTERNALS:   KEINE
C**   EINGABE-
C**   PARAMETER:   PHIS     REAL   GEOGR. BREITE DES PUNKTES IM ROT.SYS.
C**                LAMS     REAL   GEOGR. LAENGE DES PUNKTES IM ROT.SYS.
C**                POLPHI   REAL   WAHRE GEOGR. BREITE DES NORDPOLS
C**                POLLAM   REAL   WAHRE GEOGR. LAENGE DES NORDPOLS
C**   AUSGABE-
C**   PARAMETER:   WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
C**
C**   COMMON-
C**   BLOECKE  :   KEINE
C**
C**   FEHLERBE-
C**   HANDLUNG :   KEINE
C**   VERFASSER:   D.MAJEWSKI
 
      REAL        LAMS,PHIS,POLPHI,POLLAM
 
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
 
      ZSINPOL = SIN(ZPIR18*POLPHI)
      ZCOSPOL = COS(ZPIR18*POLPHI)
      ZLAMPOL = ZPIR18*POLLAM
      ZPHIS   = ZPIR18*PHIS
      ZLAMS   = LAMS
      IF(ZLAMS.GT.180.0) ZLAMS = ZLAMS - 360.0
      ZLAMS   = ZPIR18*ZLAMS
 
      ZARG1   = SIN(ZLAMPOL)*(- ZSINPOL*COS(ZLAMS)*COS(ZPHIS)  +
     1                          ZCOSPOL*           SIN(ZPHIS)) -
     2          COS(ZLAMPOL)*           SIN(ZLAMS)*COS(ZPHIS)
      ZARG2   = COS(ZLAMPOL)*(- ZSINPOL*COS(ZLAMS)*COS(ZPHIS)  +
     1                          ZCOSPOL*           SIN(ZPHIS)) +
     2          SIN(ZLAMPOL)*           SIN(ZLAMS)*COS(ZPHIS)
      IF (ABS(ZARG2).LT.1.E-30) THEN
        IF (ABS(ZARG1).LT.1.E-30) THEN
          LMSTOLM =   0.0
        ELSEIF (ZARG1.GT.0.) THEN
              LMSTOLAM =  90.0
            ELSE
              LMSTOLAM = -90.0
            ENDIF
      ELSE
        LMSTOLM = ZRPI18*ATAN2(ZARG1,ZARG2)
      ENDIF
 
      RETURN
      END

Rev 23	Rev 39
1	`PROGRAM getmima`	1	`PROGRAM getmima`
2		2
3	`c ***********************************************************************`	3	`c ***********************************************************************`
4	`c * Get minimum and maximum value of a varaiable on pressure *`	4	`c * Get minimum and maximum value of a varaiable on pressure *`
5	`c * or isentropic surface *`	5	`c * or isentropic surface *`
6	`c * Michael Sprenger / Spring, summer 2016 *`	6	`c * Michael Sprenger / Spring, summer 2016 *`
7	`c ***********************************************************************`	7	`c ***********************************************************************`
8		8
9	`use netcdf`	9	`use netcdf`
10		10
11	`implicit none`	11	`implicit none`
12		12
13	`c ----------------------------------------------------------------------`	13	`c ----------------------------------------------------------------------`
14	`c Declaration of parameters and variables`	14	`c Declaration of parameters and variables`
15	`c ----------------------------------------------------------------------`	15	`c ----------------------------------------------------------------------`
16		16
17	`c Interpolation method`	17	`c Interpolation method`
18	`integer ipom`	18	`integer ipom`
19	`parameter (ipom=0)`	19	`parameter (ipom=0)`
20		20
21	`c Flag for timecheck`	21	`c Flag for timecheck`
22	`character*80 timecheck`	22	`character*80 timecheck`
23	`parameter (timecheck = 'no' )`	23	`parameter (timecheck = 'no' )`
24		24
25	`c netCDF fields`	25	`c netCDF fields`
26	`real,allocatable, dimension (:,:) :: sp,varf`	26	`real,allocatable, dimension (:,:) :: sp,varf`
27	`real,allocatable, dimension (:,:,:) :: field,varl,tt,p3`	27	`real,allocatable, dimension (:,:,:) :: field,varl,tt,p3`
28	`real,allocatable, dimension (:) :: ak,bk`	28	`real,allocatable, dimension (:) :: ak,bk`
29	`integer stat`	29	`integer stat`
30	`character*80 cdfnam,varnam`	30	`character*80 cdfnam,varnam`
31	`character*80 vnam(200)`	31	`character*80 vnam(200)`
32	`integer nvars`	32	`integer nvars`
33	`integer cdfid,ierr`	33	`integer cdfid,ierr`
34	`real mdv`	34	`real mdv`
35	`real stagz`	35	`real stagz`
36		36
37	`c Grid description`	37	`c Grid description`
38	`real xmin,xmax ! Zonal grid extension`	38	`real xmin,xmax ! Zonal grid extension`
39	`real ymin,ymax ! Meridional grid extension`	39	`real ymin,ymax ! Meridional grid extension`
40	`integer nx,ny,nz ! Grid dimensions`	40	`integer nx,ny,nz ! Grid dimensions`
41	`real dx,dy ! Horizontal grid resolution`	41	`real dx,dy ! Horizontal grid resolution`
42	`real pollon,pollat ! Pole position`	42	`real pollon,pollat ! Pole position`
43		43
44	`c Output variables`	44	`c Output variables`
45	`real min,max ! Minimum & maximum value`	45	`real min,max ! Minimum & maximum value`
46	`real lonmin,lonmax,latmin,latmax ! Position of min & max`	46	`real lonmin,lonmax,latmin,latmax ! Position of min & max`
47		47
48	`c Auxiliary variables`	48	`c Auxiliary variables`
49	`integer iargc`	49	`integer iargc`
50	`character*(80) arg`	50	`character*(80) arg`
51	`real rd`	51	`real rd`
52	`integer i,j`	52	`integer i,j`
53	`integer minindx,minindy,maxindx,maxindy`	53	`integer minindx,minindy,maxindx,maxindy`
54	`character*80 tvar`	54	`character*80 tvar`
55	`character*1 mode`	55	`character*1 mode`
56	`character*80 clev`	56	`character*80 clev`
57	`real xphys,yphys`	57	`real xphys,yphys`
58	`real level`	58	`real level`
59	`integer flag`	59	`integer flag`
60		60
61	`c Externals`	61	`c Externals`
62	`real lmstolm,phstoph`	62	`real lmstolm,phstoph`
63	`external lmstolm,phstoph`	63	`external lmstolm,phstoph`
64		64
65	`c ----------------------------------------------------------------------`	65	`c ----------------------------------------------------------------------`
66	`c Preparation - argument handling, grid paraemters, memor allocation`	66	`c Preparation - argument handling, grid paraemters, memor allocation`
67	`c ----------------------------------------------------------------------`	67	`c ----------------------------------------------------------------------`
68		68
69	`c Check for sufficient requested arguments`	69	`c Check for sufficient requested arguments`
70	`if (iargc().lt.2) then`	70	`if (iargc().lt.2) then`
71	`print*,'USAGE: getmima filename var ',`	71	`print*,'USAGE: getmima filename var ',`
72	`> '[lev in the form Pval or Tval]'`	72	`> '[lev in the form Pval or Tval]'`
73	`call exit(1)`	73	`call exit(1)`
74	`endif`	74	`endif`
75		75
76	`c Read and transform input`	76	`c Read and transform input`
77	`call getarg(1,arg)`	77	`call getarg(1,arg)`
78	`cdfnam=trim(arg)`	78	`cdfnam=trim(arg)`
79	`call getarg(2,arg)`	79	`call getarg(2,arg)`
80	`varnam=trim(arg)`	80	`varnam=trim(arg)`
81		81
82	`if (iargc().eq.3) then`	82	`if (iargc().eq.3) then`
83	`call getarg(3,arg)`	83	`call getarg(3,arg)`
84	`mode=arg(1:1)`	84	`mode=arg(1:1)`
85	`clev=arg(2:len(arg))`	85	`clev=arg(2:len(arg))`
86	`call checkchar(clev,".",flag)`	86	`call checkchar(clev,".",flag)`
87	`if (flag.eq.0) clev=trim(clev)//"."`	87	`if (flag.eq.0) clev=trim(clev)//"."`
88	`read(clev,'(f10.2)') level`	88	`read(clev,'(f10.2)') level`
89	`else`	89	`else`
90	`mode='X'`	90	`mode='X'`
91	`level=0.`	91	`level=0.`
92	`endif`	92	`endif`
93		93
94	`c Init level type and minimum,maximum`	94	`c Init level type and minimum,maximum`
95	`min= 1.e19`	95	`min= 1.e19`
96	`max=-1.e19`	96	`max=-1.e19`
97		97
98	`c Open netCDF file`	98	`c Open netCDF file`
99	`call input_open (cdfid,cdfnam)`	99	`call input_open (cdfid,cdfnam)`
100		100
101	`c Get list of variables on file`	101	`c Get list of variables on file`
102	`call input_getvars (cdfid,vnam,nvars)`	102	`call input_getvars (cdfid,vnam,nvars)`
103		103
104	`C Get infos about data domain - in particular: dimensions`	104	`C Get infos about data domain - in particular: dimensions`
105	`nx = 1`	105	`nx = 1`
106	`ny = 1`	106	`ny = 1`
107	`nz = 1`	107	`nz = 1`
108	`call input_grid (-cdfid,varnam,xmin,xmax,ymin,ymax,dx,dy,nx,ny,`	108	`call input_grid (-cdfid,varnam,xmin,xmax,ymin,ymax,dx,dy,nx,ny,`
109	`> 0.,pollon,pollat,rd,rd,nz,rd,rd,rd,stagz,timecheck)`	109	`> 0.,pollon,pollat,rd,rd,nz,rd,rd,stagz,timecheck)`
110		110
111	`C Get memory for dynamic arrays`	111	`C Get memory for dynamic arrays`
112	`allocate(sp(nx,ny),stat=stat)`	112	`allocate(sp(nx,ny),stat=stat)`
113	`if (stat.ne.0) print,' error allocating array sp *'`	113	`if (stat.ne.0) print,' error allocating array sp *'`
114	`allocate(varf(nx,ny),stat=stat)`	114	`allocate(varf(nx,ny),stat=stat)`
115	`if (stat.ne.0) print,' error allocating array varf *'`	115	`if (stat.ne.0) print,' error allocating array varf *'`
116	`allocate(varl(nx,ny,nz),stat=stat)`	116	`allocate(varl(nx,ny,nz),stat=stat)`
117	`if (stat.ne.0) print,' error allocating array varl *'`	117	`if (stat.ne.0) print,' error allocating array varl *'`
118	`allocate(tt(nx,ny,nz),stat=stat)`	118	`allocate(tt(nx,ny,nz),stat=stat)`
119	`if (stat.ne.0) print,' error allocating array tt *'`	119	`if (stat.ne.0) print,' error allocating array tt *'`
120	`allocate(field(nx,ny,nz),stat=stat)`	120	`allocate(field(nx,ny,nz),stat=stat)`
121	`if (stat.ne.0) print,' error allocating array field *'`	121	`if (stat.ne.0) print,' error allocating array field *'`
122	`allocate(ak(nz),stat=stat)`	122	`allocate(ak(nz),stat=stat)`
123	`if (stat.ne.0) print,' error allocating array ak *'`	123	`if (stat.ne.0) print,' error allocating array ak *'`
124	`allocate(bk(nz),stat=stat)`	124	`allocate(bk(nz),stat=stat)`
125	`if (stat.ne.0) print,' error allocating array bk *'`	125	`if (stat.ne.0) print,' error allocating array bk *'`
126		126
127	`c ----------------------------------------------------------------------`	127	`c ----------------------------------------------------------------------`
128	`c Load data`	128	`c Load data`
129	`c ----------------------------------------------------------------------`	129	`c ----------------------------------------------------------------------`
130		130
131	`c Get grid info for variable - non-constant part (P, PS, AK, BK)`	131	`c Get grid info for variable - non-constant part (P, PS, AK, BK)`
132	`call input_grid`	132	`call input_grid`
133	`> (cdfid,varnam,xmin,xmax,ymin,ymax,dx,dy,nx,ny,`	133	`> (cdfid,varnam,xmin,xmax,ymin,ymax,dx,dy,nx,ny,`
134	`> 0.,pollon,pollat,varl,sp,nz,ak,bk,stagz,timecheck)`	134	`> 0.,pollon,pollat,varl,sp,nz,ak,bk,stagz,timecheck)`
135		135
136	`c Load Variable`	136	`c Load Variable`
137	`call input_wind (cdfid,varnam,field,0.,stagz,mdv,`	137	`call input_wind (cdfid,varnam,field,0.,stagz,mdv,`
138	`> xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)`	138	`> xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)`
139		139
140	`c In case of isentropic level, read temperature`	140	`c In case of isentropic level, read temperature`
141	`if (mode.eq.'T') then`	141	`if (mode.eq.'T') then`
142	`tvar="T"`	142	`tvar="T"`
143	`do i=1,nvars`	143	`do i=1,nvars`
144	`if (vnam(i).eq."THETA") tvar="THETA"`	144	`if (vnam(i).eq."THETA") tvar="THETA"`
145	`if (vnam(i).eq."TH") tvar="TH"`	145	`if (vnam(i).eq."TH") tvar="TH"`
146	`enddo`	146	`enddo`
147	`if (tvar.eq.'T') then`	147	`if (tvar.eq.'T') then`
148	`call input_wind (cdfid,tvar,tt,0.,stagz,mdv,xmin,xmax,`	148	`call input_wind (cdfid,tvar,tt,0.,stagz,mdv,xmin,xmax,`
149	`> ymin,ymax,dx,dy,nx,ny,nz,timecheck)`	149	`> ymin,ymax,dx,dy,nx,ny,nz,timecheck)`
150	`call pottemp(varl,tt,sp,nx,ny,nz,ak,bk)`	150	`call pottemp(varl,tt,sp,nx,ny,nz,ak,bk)`
151	`else`	151	`else`
152	`call input_wind (cdfid,tvar,tt,0.,stagz,mdv,xmin,xmax,`	152	`call input_wind (cdfid,tvar,tt,0.,stagz,mdv,xmin,xmax,`
153	`> ymin,ymax,dx,dy,nx,ny,nz,timecheck)`	153	`> ymin,ymax,dx,dy,nx,ny,nz,timecheck)`
154	`endif`	154	`endif`
155	`endif`	155	`endif`
156		156
157	`c ----------------------------------------------------------------------`	157	`c ----------------------------------------------------------------------`
158	`c Do interpolation and get minimum, maximum`	158	`c Do interpolation and get minimum, maximum`
159	`c ----------------------------------------------------------------------`	159	`c ----------------------------------------------------------------------`
160		160
161	`C If required do interpolation on pressure or theta level`	161	`C If required do interpolation on pressure or theta level`
162	`if (mode.eq.'P') then`	162	`if (mode.eq.'P') then`
163	`call vipo(field,varl,level,varf,nx,ny,nz,mdv,ipom)`	163	`call vipo(field,varl,level,varf,nx,ny,nz,mdv,ipom)`
164	`elseif (mode.eq.'T') then`	164	`elseif (mode.eq.'T') then`
165	`call thipo(field,varl,level,varf,nx,ny,nz,mdv,ipom)`	165	`call thipo(field,varl,level,varf,nx,ny,nz,mdv,ipom)`
166	`elseif (mode.eq.'L') then`	166	`elseif (mode.eq.'L') then`
167	`do i=1,nx`	167	`do i=1,nx`
168	`do j=1,ny`	168	`do j=1,ny`
169	`varf(i,j) = field(i,j,nint(level))`	169	`varf(i,j) = field(i,j,nint(level))`
170	`enddo`	170	`enddo`
171	`enddo`	171	`enddo`
172	`elseif (mode.eq.'X') then`	172	`elseif (mode.eq.'X') then`
173	`do i=1,nx`	173	`do i=1,nx`
174	`do j=1,ny`	174	`do j=1,ny`
175	`varf(i,j) = field(i,j,1)`	175	`varf(i,j) = field(i,j,1)`
176	`enddo`	176	`enddo`
177	`enddo`	177	`enddo`
178	`endif`	178	`endif`
179		179
180	`C Determine minimum and maximum`	180	`C Determine minimum and maximum`
181	`do i=1,nx`	181	`do i=1,nx`
182	`do j=1,ny`	182	`do j=1,ny`
183		183
184	`if ((varf(i,j).ne.mdv).and.(varf(i,j).lt.min)) then`	184	`if ((varf(i,j).ne.mdv).and.(varf(i,j).lt.min)) then`
185	`min = varf(i,j)`	185	`min = varf(i,j)`
186	`minindx = i`	186	`minindx = i`
187	`minindy = j`	187	`minindy = j`
188	`endif`	188	`endif`
189	`if ((varf(i,j).ne.mdv).and.(varf(i,j).gt.max)) then`	189	`if ((varf(i,j).ne.mdv).and.(varf(i,j).gt.max)) then`
190	`max = varf(i,j)`	190	`max = varf(i,j)`
191	`maxindx = i`	191	`maxindx = i`
192	`maxindy = j`	192	`maxindy = j`
193	`endif`	193	`endif`
194	`enddo`	194	`enddo`
195	`enddo`	195	`enddo`
196	`lonmin = xmin + real(minindx-1) * dx`	196	`lonmin = xmin + real(minindx-1) * dx`
197	`latmin = ymin + real(minindy-1) * dy`	197	`latmin = ymin + real(minindy-1) * dy`
198	`lonmax = xmin + real(maxindx-1) * dx`	198	`lonmax = xmin + real(maxindx-1) * dx`
199	`latmax = ymin + real(maxindy-1) * dy`	199	`latmax = ymin + real(maxindy-1) * dy`
200		200
201	`c Rotate position to true longitude/latitude`	201	`c Rotate position to true longitude/latitude`
202	`if ((pollon.ne.0.).or.(pollat.ne.90.)) then`	202	`if ((pollon.ne.0.).or.(pollat.ne.90.)) then`
203	`xphys=lmstolm(latmin,lonmin,pollat,pollon)`	203	`xphys=lmstolm(latmin,lonmin,pollat,pollon)`
204	`yphys=phstoph(latmin,lonmin,pollat,pollon)`	204	`yphys=phstoph(latmin,lonmin,pollat,pollon)`
205	`lonmin=xphys`	205	`lonmin=xphys`
206	`latmin=yphys`	206	`latmin=yphys`
207	`xphys=lmstolm(latmax,lonmax,pollat,pollon)`	207	`xphys=lmstolm(latmax,lonmax,pollat,pollon)`
208	`yphys=phstoph(latmax,lonmax,pollat,pollon)`	208	`yphys=phstoph(latmax,lonmax,pollat,pollon)`
209	`lonmax=xphys`	209	`lonmax=xphys`
210	`latmax=yphys`	210	`latmax=yphys`
211	`endif`	211	`endif`
212		212
213	`c Write output`	213	`c Write output`
214	`write(*,101) min,max,lonmin,latmin,nint(level),`	214	`write(*,101) min,max,lonmin,latmin,nint(level),`
215	`> lonmax,latmax,nint(level)`	215	`> lonmax,latmax,nint(level)`
216	`101 format(2f10.3,2f8.2,i6,2f8.2,i6)`	216	`101 format(2f10.3,2f8.2,i6,2f8.2,i6)`
217		217
218	`c Close netCDF file`	218	`c Close netCDF file`
219	`call input_close(cdfid)`	219	`call input_close(cdfid)`
220		220
221	`end`	221	`end`
222		222
223	`c ----------------------------------------------------------------------`	223	`c ----------------------------------------------------------------------`
224	`C Interpolates the 3d variable var3d on the pressure surface defined`	224	`C Interpolates the 3d variable var3d on the pressure surface defined`
225	`C by lev of the variable varl. The interpolated field is`	225	`C by lev of the variable varl. The interpolated field is`
226	`C returned as var.`	226	`C returned as var.`
227	`C ipom determines the way of vertical interpolation:`	227	`C ipom determines the way of vertical interpolation:`
228	`C ipom=0 is for linear interpolation`	228	`C ipom=0 is for linear interpolation`
229	`C ipom=1 is for logarithmic interpolation`	229	`C ipom=1 is for logarithmic interpolation`
230	`c ----------------------------------------------------------------------`	230	`c ----------------------------------------------------------------------`
231		231
232	`subroutine vipo(var3d,varl,lev,var,nx,ny,nz,mdv,ipom)`	232	`subroutine vipo(var3d,varl,lev,var,nx,ny,nz,mdv,ipom)`
233		233
234	`integer nx,ny,nz,ipom`	234	`integer nx,ny,nz,ipom`
235	`real lev,mdv`	235	`real lev,mdv`
236	`real var3d(nx,ny,nz),varl(nx,ny,nz),var(nx,ny)`	236	`real var3d(nx,ny,nz),varl(nx,ny,nz),var(nx,ny)`
237		237
238	`integer i,j,k`	238	`integer i,j,k`
239	`real kind`	239	`real kind`
240	`real int3dm`	240	`real int3dm`
241		241
242	`do i=1,nx`	242	`do i=1,nx`
243	`do j=1,ny`	243	`do j=1,ny`
244		244
245	`do k=1,nz-1`	245	`do k=1,nz-1`
246	`if ((varl(i,j,k).ge.lev).and.(varl(i,j,k+1).le.lev)) then`	246	`if ((varl(i,j,k).ge.lev).and.(varl(i,j,k+1).le.lev)) then`
247	`kind=float(k)+(varl(i,j,k)-lev)/`	247	`kind=float(k)+(varl(i,j,k)-lev)/`
248	`> (varl(i,j,k)-varl(i,j,k+1))`	248	`> (varl(i,j,k)-varl(i,j,k+1))`
249	`goto 100`	249	`goto 100`
250	`endif`	250	`endif`
251	`enddo`	251	`enddo`
252	`100 continue`	252	`100 continue`
253		253
254	`var(i,j)=int3dm(var3d,nx,ny,nz,float(i),float(j),kind,mdv)`	254	`var(i,j)=int3dm(var3d,nx,ny,nz,float(i),float(j),kind,mdv)`
255		255
256	`enddo`	256	`enddo`
257	`enddo`	257	`enddo`
258		258
259	`return`	259	`return`
260	`end`	260	`end`
261		261
262	`c ----------------------------------------------------------------------`	262	`c ----------------------------------------------------------------------`
263	`C Interpolates the 3d variable var3d on the isentropic surface defined`	263	`C Interpolates the 3d variable var3d on the isentropic surface defined`
264	`C by lev of the variable varl. The interpolated field is`	264	`C by lev of the variable varl. The interpolated field is`
265	`C returned as var.`	265	`C returned as var.`
266	`C ipom determines the way of vertical interpolation:`	266	`C ipom determines the way of vertical interpolation:`
267	`C ipom=0 is for linear interpolation`	267	`C ipom=0 is for linear interpolation`
268	`C ipom=1 is for logarithmic interpolation`	268	`C ipom=1 is for logarithmic interpolation`
269	`c ----------------------------------------------------------------------`	269	`c ----------------------------------------------------------------------`
270		270
271	`subroutine thipo(var3d,th3d,lev,var,nx,ny,nz,mdv,mode)`	271	`subroutine thipo(var3d,th3d,lev,var,nx,ny,nz,mdv,mode)`
272		272
273	`integer nx,ny,nz,mode`	273	`integer nx,ny,nz,mode`
274	`real lev,mdv`	274	`real lev,mdv`
275	`real var3d(nx,ny,nz),th3d(nx,ny,nz),var(nx,ny)`	275	`real var3d(nx,ny,nz),th3d(nx,ny,nz),var(nx,ny)`
276		276
277	`integer i,j,k`	277	`integer i,j,k`
278	`real kind`	278	`real kind`
279	`real int3dm`	279	`real int3dm`
280		280
281	`do i=1,nx`	281	`do i=1,nx`
282	`do j=1,ny`	282	`do j=1,ny`
283		283
284	`do k=1,nz-1`	284	`do k=1,nz-1`
285	`if ((th3d(i,j,k).le.lev).and.(th3d(i,j,k+1).ge.lev)) then`	285	`if ((th3d(i,j,k).le.lev).and.(th3d(i,j,k+1).ge.lev)) then`
286	`kind=float(k)+(th3d(i,j,k)-lev)/`	286	`kind=float(k)+(th3d(i,j,k)-lev)/`
287	`> (th3d(i,j,k)-th3d(i,j,k+1))`	287	`> (th3d(i,j,k)-th3d(i,j,k+1))`
288	`goto 100`	288	`goto 100`
289	`endif`	289	`endif`
290	`enddo`	290	`enddo`
291	`100 continue`	291	`100 continue`
292		292
293	`var(i,j)=int3dm(var3d,nx,ny,nz,float(i),float(j),kind,mdv)`	293	`var(i,j)=int3dm(var3d,nx,ny,nz,float(i),float(j),kind,mdv)`
294		294
295	`enddo`	295	`enddo`
296	`enddo`	296	`enddo`
297		297
298	`return`	298	`return`
299	`end`	299	`end`
300		300
301	`c ----------------------------------------------------------------------`	301	`c ----------------------------------------------------------------------`
302	`c Check whether character appears within string`	302	`c Check whether character appears within string`
303	`c ----------------------------------------------------------------------`	303	`c ----------------------------------------------------------------------`
304		304
305	`subroutine checkchar(string,char,flag)`	305	`subroutine checkchar(string,char,flag)`
306		306
307	`character() string`	307	`character() string`
308	`character*(1) char`	308	`character*(1) char`
309	`integer n,flag`	309	`integer n,flag`
310		310
311	`flag=0`	311	`flag=0`
312	`do n=1,len(string)`	312	`do n=1,len(string)`
313	`if (string(n:n).eq.char) then`	313	`if (string(n:n).eq.char) then`
314	`flag=n`	314	`flag=n`
315	`return`	315	`return`
316	`endif`	316	`endif`
317	`enddo`	317	`enddo`
318	`end`	318	`end`
319		319
320	`c ----------------------------------------------------------------------`	320	`c ----------------------------------------------------------------------`
321	`c 3D interpolation`	321	`c 3D interpolation`
322	`c ----------------------------------------------------------------------`	322	`c ----------------------------------------------------------------------`
323		323
324	`real function int3d(ar,n1,n2,n3,rid,rjd,rkd)`	324	`real function int3d(ar,n1,n2,n3,rid,rjd,rkd)`
325		325
326	`c Purpose:`	326	`c Purpose:`
327	`c This subroutine interpolates a 3d-array to an arbitrary`	327	`c This subroutine interpolates a 3d-array to an arbitrary`
328	`c location within the grid.`	328	`c location within the grid.`
329	`c Arguments:`	329	`c Arguments:`
330	`c ar real input surface pressure, define as ar(n1,n2,n3)`	330	`c ar real input surface pressure, define as ar(n1,n2,n3)`
331	`c n1,n2,n3 int input dimensions of ar`	331	`c n1,n2,n3 int input dimensions of ar`
332	`c ri,rj,rk real input grid location to be interpolated to`	332	`c ri,rj,rk real input grid location to be interpolated to`
333	`c History:`	333	`c History:`
334		334
335	`c argument declarations`	335	`c argument declarations`
336	`integer n1,n2,n3`	336	`integer n1,n2,n3`
337	`real ar(n1,n2,n3), rid,rjd,rkd`	337	`real ar(n1,n2,n3), rid,rjd,rkd`
338		338
339	`c local declarations`	339	`c local declarations`
340	`integer i,j,k,ip1,jp1,kp1,ih,jh,kh`	340	`integer i,j,k,ip1,jp1,kp1,ih,jh,kh`
341	`real frac0i,frac0j,frac0k,frac1i,frac1j,frac1k,ri,rj,rk`	341	`real frac0i,frac0j,frac0k,frac1i,frac1j,frac1k,ri,rj,rk`
342		342
343	`c do linear interpolation`	343	`c do linear interpolation`
344	`ri=amax1(1.,amin1(float(n1),rid))`	344	`ri=amax1(1.,amin1(float(n1),rid))`
345	`rj=amax1(1.,amin1(float(n2),rjd))`	345	`rj=amax1(1.,amin1(float(n2),rjd))`
346	`rk=amax1(1.,amin1(float(n3),rkd))`	346	`rk=amax1(1.,amin1(float(n3),rkd))`
347	`ih=nint(ri)`	347	`ih=nint(ri)`
348	`jh=nint(rj)`	348	`jh=nint(rj)`
349	`kh=nint(rk)`	349	`kh=nint(rk)`
350		350
351	`c Check for interpolation in i`	351	`c Check for interpolation in i`
352	`* if (abs(float(ih)-ri).lt.1.e-3) then`	352	`* if (abs(float(ih)-ri).lt.1.e-3) then`
353	`* i =ih`	353	`* i =ih`
354	`* ip1=ih`	354	`* ip1=ih`
355	`* else`	355	`* else`
356	`i =min0(int(ri),n1-1)`	356	`i =min0(int(ri),n1-1)`
357	`ip1=i+1`	357	`ip1=i+1`
358	`* endif`	358	`* endif`
359		359
360	`c Check for interpolation in j`	360	`c Check for interpolation in j`
361	`if (abs(float(jh)-rj).lt.1.e-3) then`	361	`if (abs(float(jh)-rj).lt.1.e-3) then`
362	`j =jh`	362	`j =jh`
363	`jp1=jh`	363	`jp1=jh`
364	`else`	364	`else`
365	`j =min0(int(rj),n2-1)`	365	`j =min0(int(rj),n2-1)`
366	`jp1=j+1`	366	`jp1=j+1`
367	`endif`	367	`endif`
368		368
369	`c Check for interpolation in k`	369	`c Check for interpolation in k`
370	`* if (abs(float(kh)-rk).lt.1.e-3) then`	370	`* if (abs(float(kh)-rk).lt.1.e-3) then`
371	`* k =kh`	371	`* k =kh`
372	`* kp1=kh`	372	`* kp1=kh`
373	`* else`	373	`* else`
374	`k =min0(int(rk),n3-1)`	374	`k =min0(int(rk),n3-1)`
375	`kp1=k+1`	375	`kp1=k+1`
376	`* endif`	376	`* endif`
377		377
378	`if (k.eq.kp1) then`	378	`if (k.eq.kp1) then`
379	`c no interpolation in k`	379	`c no interpolation in k`
380	`if ((i.eq.ip1).and.(j.eq.jp1)) then`	380	`if ((i.eq.ip1).and.(j.eq.jp1)) then`
381	`c no interpolation at all`	381	`c no interpolation at all`
382	`int3d=ar(i,j,k)`	382	`int3d=ar(i,j,k)`
383	`c print *,'int3d 00: ',rid,rjd,rkd,int3d`	383	`c print *,'int3d 00: ',rid,rjd,rkd,int3d`
384	`else`	384	`else`
385	`c horizontal interpolation only`	385	`c horizontal interpolation only`
386	`frac0i=ri-float(i)`	386	`frac0i=ri-float(i)`
387	`frac0j=rj-float(j)`	387	`frac0j=rj-float(j)`
388	`frac1i=1.-frac0i`	388	`frac1i=1.-frac0i`
389	`frac1j=1.-frac0j`	389	`frac1j=1.-frac0j`
390	`int3d = ar(i ,j ,k ) * frac1i * frac1j`	390	`int3d = ar(i ,j ,k ) * frac1i * frac1j`
391	`& + ar(i ,jp1,k ) * frac1i * frac0j`	391	`& + ar(i ,jp1,k ) * frac1i * frac0j`
392	`& + ar(ip1,j ,k ) * frac0i * frac1j`	392	`& + ar(ip1,j ,k ) * frac0i * frac1j`
393	`& + ar(ip1,jp1,k ) * frac0i * frac0j`	393	`& + ar(ip1,jp1,k ) * frac0i * frac0j`
394	`c print *,'int3d 10: ',rid,rjd,rkd,int3d`	394	`c print *,'int3d 10: ',rid,rjd,rkd,int3d`
395	`endif`	395	`endif`
396	`else`	396	`else`
397	`frac0k=rk-float(k)`	397	`frac0k=rk-float(k)`
398	`frac1k=1.-frac0k`	398	`frac1k=1.-frac0k`
399	`if ((i.eq.ip1).and.(j.eq.jp1)) then`	399	`if ((i.eq.ip1).and.(j.eq.jp1)) then`
400	`c vertical interpolation only`	400	`c vertical interpolation only`
401	`int3d = ar(i ,j ,k ) * frac1k`	401	`int3d = ar(i ,j ,k ) * frac1k`
402	`& + ar(i ,j ,kp1) * frac0k`	402	`& + ar(i ,j ,kp1) * frac0k`
403	`c print *,'int3d 01: ',rid,rjd,rkd,int3d`	403	`c print *,'int3d 01: ',rid,rjd,rkd,int3d`
404	`else`	404	`else`
405	`c full 3d interpolation`	405	`c full 3d interpolation`
406	`frac0i=ri-float(i)`	406	`frac0i=ri-float(i)`
407	`frac0j=rj-float(j)`	407	`frac0j=rj-float(j)`
408	`frac1i=1.-frac0i`	408	`frac1i=1.-frac0i`
409	`frac1j=1.-frac0j`	409	`frac1j=1.-frac0j`
410	`int3d = ar(i ,j ,k ) * frac1i * frac1j * frac1k`	410	`int3d = ar(i ,j ,k ) * frac1i * frac1j * frac1k`
411	`& + ar(i ,jp1,k ) * frac1i * frac0j * frac1k`	411	`& + ar(i ,jp1,k ) * frac1i * frac0j * frac1k`
412	`& + ar(ip1,j ,k ) * frac0i * frac1j * frac1k`	412	`& + ar(ip1,j ,k ) * frac0i * frac1j * frac1k`
413	`& + ar(ip1,jp1,k ) * frac0i * frac0j * frac1k`	413	`& + ar(ip1,jp1,k ) * frac0i * frac0j * frac1k`
414	`& + ar(i ,j ,kp1) * frac1i * frac1j * frac0k`	414	`& + ar(i ,j ,kp1) * frac1i * frac1j * frac0k`
415	`& + ar(i ,jp1,kp1) * frac1i * frac0j * frac0k`	415	`& + ar(i ,jp1,kp1) * frac1i * frac0j * frac0k`
416	`& + ar(ip1,j ,kp1) * frac0i * frac1j * frac0k`	416	`& + ar(ip1,j ,kp1) * frac0i * frac1j * frac0k`
417	`& + ar(ip1,jp1,kp1) * frac0i * frac0j * frac0k`	417	`& + ar(ip1,jp1,kp1) * frac0i * frac0j * frac0k`
418	`c print *,'int3d 11: ',rid,rjd,rkd,int3d`	418	`c print *,'int3d 11: ',rid,rjd,rkd,int3d`
419	`endif`	419	`endif`
420	`endif`	420	`endif`
421	`end`	421	`end`
422		422
423	`c ----------------------------------------------------------------------`	423	`c ----------------------------------------------------------------------`
424	`c 3D interpolation including missing data check`	424	`c 3D interpolation including missing data check`
425	`c ----------------------------------------------------------------------`	425	`c ----------------------------------------------------------------------`
426		426
427	`real function int3dm(ar,n1,n2,n3,rid,rjd,rkd,misdat)`	427	`real function int3dm(ar,n1,n2,n3,rid,rjd,rkd,misdat)`
428		428
429	`c Purpose:`	429	`c Purpose:`
430	`c This subroutine interpolates a 3d-array to an arbitrary`	430	`c This subroutine interpolates a 3d-array to an arbitrary`
431	`c location within the grid. The interpolation includes the`	431	`c location within the grid. The interpolation includes the`
432	`c testing of the missing data flag 'misdat'.`	432	`c testing of the missing data flag 'misdat'.`
433	`c Arguments:`	433	`c Arguments:`
434	`c ar real input surface pressure, define as ar(n1,n2,n3)`	434	`c ar real input surface pressure, define as ar(n1,n2,n3)`
435	`c n1,n2,n3 int input dimensions of ar`	435	`c n1,n2,n3 int input dimensions of ar`
436	`c ri,rj,rk real input grid location to be interpolated to`	436	`c ri,rj,rk real input grid location to be interpolated to`
437	`c misdat real input missing data flag (on if misdat<>0)`	437	`c misdat real input missing data flag (on if misdat<>0)`
438	`c Warning:`	438	`c Warning:`
439	`c This routine has not yet been seriously tested`	439	`c This routine has not yet been seriously tested`
440	`c History:`	440	`c History:`
441		441
442	`c argument declarations`	442	`c argument declarations`
443	`integer n1,n2,n3`	443	`integer n1,n2,n3`
444	`real ar(n1,n2,n3), rid,rjd,rkd, misdat`	444	`real ar(n1,n2,n3), rid,rjd,rkd, misdat`
445		445
446	`c local declarations`	446	`c local declarations`
447	`integer i,j,k,ip1,jp1,kp1,ih,jh,kh`	447	`integer i,j,k,ip1,jp1,kp1,ih,jh,kh`
448	`real frac0i,frac0j,frac0k,frac1i,frac1j,frac1k,ri,rj,rk,int3d`	448	`real frac0i,frac0j,frac0k,frac1i,frac1j,frac1k,ri,rj,rk,int3d`
449		449
450	`c check if routine without missing data checking can be called instead`	450	`c check if routine without missing data checking can be called instead`
451	`if (misdat.eq.0.) then`	451	`if (misdat.eq.0.) then`
452	`int3dm=int3d(ar,n1,n2,n3,rid,rjd,rkd)`	452	`int3dm=int3d(ar,n1,n2,n3,rid,rjd,rkd)`
453	`return`	453	`return`
454	`endif`	454	`endif`
455		455
456	`c do linear interpolation`	456	`c do linear interpolation`
457	`ri=amax1(1.,amin1(float(n1),rid))`	457	`ri=amax1(1.,amin1(float(n1),rid))`
458	`rj=amax1(1.,amin1(float(n2),rjd))`	458	`rj=amax1(1.,amin1(float(n2),rjd))`
459	`rk=amax1(1.,amin1(float(n3),rkd))`	459	`rk=amax1(1.,amin1(float(n3),rkd))`
460	`ih=nint(ri)`	460	`ih=nint(ri)`
461	`jh=nint(rj)`	461	`jh=nint(rj)`
462	`kh=nint(rk)`	462	`kh=nint(rk)`
463		463
464	`c Check for interpolation in i`	464	`c Check for interpolation in i`
465	`* if (abs(float(ih)-ri).lt.1.e-3) then`	465	`* if (abs(float(ih)-ri).lt.1.e-3) then`
466	`* i =ih`	466	`* i =ih`
467	`* ip1=ih`	467	`* ip1=ih`
468	`* else`	468	`* else`
469	`i =min0(int(ri),n1-1)`	469	`i =min0(int(ri),n1-1)`
470	`ip1=i+1`	470	`ip1=i+1`
471	`* endif`	471	`* endif`
472		472
473	`c Check for interpolation in j`	473	`c Check for interpolation in j`
474	`* if (abs(float(jh)-rj).lt.1.e-3) then`	474	`* if (abs(float(jh)-rj).lt.1.e-3) then`
475	`* j =jh`	475	`* j =jh`
476	`* jp1=jh`	476	`* jp1=jh`
477	`* else`	477	`* else`
478	`j =min0(int(rj),n2-1)`	478	`j =min0(int(rj),n2-1)`
479	`jp1=j+1`	479	`jp1=j+1`
480	`* endif`	480	`* endif`
481		481
482	`c Check for interpolation in k`	482	`c Check for interpolation in k`
483	`* if (abs(float(kh)-rk).lt.1.e-3) then`	483	`* if (abs(float(kh)-rk).lt.1.e-3) then`
484	`* k =kh`	484	`* k =kh`
485	`* kp1=kh`	485	`* kp1=kh`
486	`* else`	486	`* else`
487	`k =min0(int(rk),n3-1)`	487	`k =min0(int(rk),n3-1)`
488	`kp1=k+1`	488	`kp1=k+1`
489	`* endif`	489	`* endif`
490		490
491	`if (k.eq.kp1) then`	491	`if (k.eq.kp1) then`
492	`c no interpolation in k`	492	`c no interpolation in k`
493	`if ((i.eq.ip1).and.(j.eq.jp1)) then`	493	`if ((i.eq.ip1).and.(j.eq.jp1)) then`
494	`c no interpolation at all`	494	`c no interpolation at all`
495	`if (misdat.eq.ar(i,j,k)) then`	495	`if (misdat.eq.ar(i,j,k)) then`
496	`int3dm=misdat`	496	`int3dm=misdat`
497	`else`	497	`else`
498	`int3dm=ar(i,j,k)`	498	`int3dm=ar(i,j,k)`
499	`endif`	499	`endif`
500	`c print *,'int3dm 00: ',rid,rjd,rkd,int3dm`	500	`c print *,'int3dm 00: ',rid,rjd,rkd,int3dm`

Subversion Repositories lagranto.ecmwf

(root)/trunk/goodies/getmima.f – Rev 23 → 39