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


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

Subversion Repositories lagranto.ecmwf

(root)/trunk/goodies/getmima.f – Rev 3 → 13