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

Subversion Repositories lagranto.um

(root)/tags/1.0/goodies/getmima.f – Rev 3 → 4