WebSVN – lagranto.ecmwf – Diff – /trunk/lib/inter.f


c      *************************************************************
c      * This package provides a general interpolaton routine      *
c      *************************************************************	
 
c     The main interface routines are:
c         get_index3,4 : get the grid indices for interpolation
c         int_index3,4 : interpolate to the grid position
 
c     -------------------------------------------------------------
c     Get index in grid space for interpolation
c     -------------------------------------------------------------
 
      subroutine get_index4 (rid,rjd,rkd,xpo,ypo,ppo,rtp,
     >                       vert0,vert1,surf0,surf1,mode,
     >                       nx,ny,nz,lonw,lats,dlon,dlat,misdat)
 
c     Purpose:
c        This subroutine determines the indices (rid,rjd,rkd) in grid 
c        space for a point in physical space (xpo,ypo,ppo). The 
c        horizontal grid is specified by the south-west point (lats,lonw)
c        and the grid spacing (dlat,dlon). The vertical grid is given
c        by <vert(n1,n2,n3)>. The lower boundary (typicall surface 
c        pressure) is given by <surf(n1,n2)>.
c     Arguments:
c        rid,rjd,rkd  real  output   grid location to be interpolated to
c        xpo,ypo,ppo  real  input    physical coordinates
c        rtp          real  input    relative time position (0=beginning, 1=end)
c        n1,n2,n3     int   input    grid dimensions in x-, y- and p-direction
c        lats,lonw    real  input    south and west boundary of grid space
c        vert         real  input    vertical coordinate grid
c        surf         real  input    lower boundary (surface pressure)
c        mode         int   input    direction of vertical axis (p=1,th=-1)
c                                        1: linear, 1 -> nz (th)
c                                        2: linear, nz -> 1 (pv)
c                                        3: binary (p)
 
      implicit none
 
c     Declartion of function parameters
      integer   nx,ny,nz
      real      xpo,ypo,ppo,rtp
      real      vert0(nx*ny*nz),vert1(nx*ny*nz)
      real      surf0(nx*ny)   ,surf1(nx*ny*nz)
      real      rid,rjd,rkd
      real      dlat,dlon,lats,lonw
      real      misdat
      integer   mode
 
c     Set numerical parameters
      real      eps
      parameter (eps=1.e-8)
 
c     Auxiliary variables
      real      rid0,rjd0,rkd0,rid1,rjd1,rkd1
      integer   i
 
c     Externals 
      real      int_time
      external  int_time
 
c     Get the inidices
      if (abs(rtp).lt.eps) then
         call  get_index3 (rid,rjd,rkd,xpo,ypo,ppo,mode,
     >                     vert0,surf0,nx,ny,nz,lonw,lats,dlon,dlat)
      elseif (abs(rtp-1.).lt.eps) then
         call  get_index3 (rid,rjd,rkd,xpo,ypo,ppo,mode,
     >                     vert1,surf1,nx,ny,nz,lonw,lats,dlon,dlat)
      else
         call  get_index3 (rid0,rjd0,rkd0,xpo,ypo,ppo,mode,
     >                     vert0,surf0,nx,ny,nz,lonw,lats,dlon,dlat)
         call  get_index3 (rid1,rjd1,rkd1,xpo,ypo,ppo,mode,
     >                     vert1,surf1,nx,ny,nz,lonw,lats,dlon,dlat)
         rid = int_time (rid0,rid1,rtp,misdat)
         rjd = int_time (rjd0,rjd1,rtp,misdat)
         rkd = int_time (rkd0,rkd1,rtp,misdat)
 
      endif
 
      end
 
c     -------------------------------------------------------------
c     Interpolate to an arbitrary position in grid space and time
c     -------------------------------------------------------------
 
      real function int_index4 (ar0,ar1,n1,n2,n3,rid,rjd,rkd,rtp,misdat)
 
c     Purpose:
c        This subroutine interpolates a 3d-array to an arbitrary
c        location within the grid including a linear time-interpolation. 
c     Arguments:
c        rid,rjd,rkd  real  output   grid location to be interpolated to
c        xpo,ypo,ppo  real  input    physical coordinates
c        n1,n2,n3     int   input    grid dimensions in x-, y- and p-direction
c        lats,lonw    real  input    south and west boundary of grid space
c        vert         real  input    vertical coordinate grid
c        surf         real  input    lower boundary (surface pressure)
 
      implicit none
 
c     Declartion of function parameters
      integer   n1,n2,n3
      real      ar0(n1*n2*n3),ar1(n1*n2*n3)
      real      rid,rjd,rkd
      real      rtp
      real      misdat
 
c     Set numerical parameters
      real      eps
      parameter (eps=1.e-8)
 
c     Externals  
      real      int_index3,int_time
      external  int_index3,int_time
 
c     Auxiliary variables
      real      val0,val1,val
 
c     Do the 3d-interpolation
      if (abs(rtp).lt.eps) then
         val = int_index3 (ar0,n1,n2,n3,rid,rjd,rkd,misdat)
      elseif (abs(rtp-1.).lt.eps) then
         val = int_index3 (ar1,n1,n2,n3,rid,rjd,rkd,misdat)
      else
         val0 = int_index3 (ar0,n1,n2,n3,rid,rjd,rkd,misdat)
         val1 = int_index3 (ar1,n1,n2,n3,rid,rjd,rkd,misdat)
         val  = int_time (val0,val1,rtp,misdat)
      endif
 
c     Return value
      int_index4 = val
 
      return
      end
 
 
c     -------------------------------------------------------------
c     Interpolate to an arbitrary position in grid space
c     -------------------------------------------------------------
 
      real function int_index3 (ar,n1,n2,n3,rid,rjd,rkd,misdat)
 
c     Purpose:
c        This subroutine interpolates a 3d-array to an arbitrary
c        location within the grid. The interpolation includes the 
c        testing of the missing data flag 'misdat'. If one dimension
c        is 1, a 2d-interpolation is performed; if two dimensions
c        are 1, it is a 1d-interpolation; if all three dimensions are
c        1, no interpolation is performed and the input value is
c        returned.
c     Arguments:
c        ar        real  input   input data array
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)
 
      implicit none
 
c     Declartion of function parameters 
      integer   n1,n2,n3
      real      ar(n1*n2*n3)
      real      rid,rjd,rkd
      real      misdat
 
c     Set numerical parameters
      real      eps
      parameter (eps=1.e-8)
 
c     Local variables
      integer   i,j,k,ip1,jp1,kp1
      real      frac0i,frac0j,frac0k,frac1i,frac1j,frac1k
      real      ri,rj,rk
      real      val000,val001,val010,val011,val100,val101,val110,val111
      real      frc000,frc001,frc010,frc011,frc100,frc101,frc110,frc111
      real      frc
      real      mdv
      real      val
 
c     Elementary test for dimensions
      if ( (n1.lt.1).or.(n2.lt.1).or.(n3.lt.1) ) then
         print*,'Invalid grid dimensions ',n1,n2,n3
         stop
      endif
 
c     Activate or inactive the missing data check (quick and dirty)
      if (misdat.ne.0.) then
         mdv = misdat
      else
         mdv = 257.22725394015
      endif
 
c     Bring the indices into the grid space
      ri = amax1(1.,amin1(float(n1),rid))
      rj = amax1(1.,amin1(float(n2),rjd))
      rk = amax1(1.,amin1(float(n3),rkd))
 
c     Get the index of the west-south-bottom corner of the box
      i   = min0(int(ri),n1-1)
      ip1 = i+1
      j   = min0(int(rj),n2-1)
      jp1 = j+1
      k   = min0(int(rk),n3-1)
      kp1 = k+1
 
c     Special handling for 2d arrays
      if (n3.eq.1) then
         k=1
         kp1=1
      endif
 
c     Get location relative to grid box
      if ( i.ne.ip1 ) then
         frac0i = ri-float(i)
         frac1i = 1.-frac0i
      else
         frac0i = 0.
         frac1i = 1.
      endif
      if ( j.ne.jp1 ) then
         frac0j = rj-float(j)
         frac1j = 1.-frac0j
      else
         frac0j = 0.
         frac1j = 1.
      endif
      if ( k.ne.kp1 ) then
         frac0k = rk-float(k)
         frac1k = 1.-frac0k
      else
         frac0k = 0.
         frac1k = 1.
      endif
 
c     On a grid point - take the grid point value 
      if ( ( abs(frac0i).lt.eps ).and.
     >     ( abs(frac0j).lt.eps ).and.
     >     ( abs(frac0k).lt.eps ) ) then
         
         val = ar( i + n1*(j -1) + n1*n2*(k -1) )
         goto 100
         
      endif
 
c     Init the fractions
      frc000 = frac1i * frac1j * frac1k
      frc001 = frac0i * frac1j * frac1k
      frc010 = frac1i * frac0j * frac1k
      frc011 = frac0i * frac0j * frac1k
      frc100 = frac1i * frac1j * frac0k
      frc101 = frac0i * frac1j * frac0k
      frc110 = frac1i * frac0j * frac0k
      frc111 = frac0i * frac0j * frac0k
 
c     Init the values
      val000 = ar( i   + n1*(j  -1) + n1*n2*(k  -1) )
      val001 = ar( ip1 + n1*(j  -1) + n1*n2*(k  -1) )
      val010 = ar( i   + n1*(jp1-1) + n1*n2*(k  -1) )
      val011 = ar( ip1 + n1*(jp1-1) + n1*n2*(k  -1) )
      val100 = ar( i   + n1*(j  -1) + n1*n2*(kp1-1) )
      val101 = ar( ip1 + n1*(j  -1) + n1*n2*(kp1-1) )
      val110 = ar( i   + n1*(jp1-1) + n1*n2*(kp1-1) )
      val111 = ar( ip1 + n1*(jp1-1) + n1*n2*(kp1-1) )
 
c     Handle missing data
      if ( abs(val000-mdv).lt.eps ) frc000 = 0.
      if ( abs(val001-mdv).lt.eps ) frc001 = 0.
      if ( abs(val010-mdv).lt.eps ) frc010 = 0.
      if ( abs(val011-mdv).lt.eps ) frc011 = 0.
      if ( abs(val100-mdv).lt.eps ) frc100 = 0.
      if ( abs(val101-mdv).lt.eps ) frc101 = 0.
      if ( abs(val110-mdv).lt.eps ) frc110 = 0.
      if ( abs(val111-mdv).lt.eps ) frc111 = 0.
 
c     Build the final value
      frc = frc000 + frc001 + frc010 + frc011 + 
     >      frc100 + frc101 + frc110 + frc111   
      if ( frc.gt.0. ) then
         val = 1./frc * ( frc000 * val000 + frc001 * val001 +
     >                    frc010 * val010 + frc011 * val011 +
     >                    frc100 * val100 + frc101 * val101 +
     >                    frc110 * val110 + frc111 * val111 )
      else
         val = misdat
      endif
 
c     Return the value 
 100  continue
 
      int_index3 = val
 
      end
 
 
c     -------------------------------------------------------------
c     Time interpolation
c     -------------------------------------------------------------
 
      real function int_time (val0,val1,reltpos,misdat)
 
c     Purpose:
c        This subroutine interpolates linearly in time between two
c        values.
c     Arguments:
c        val0      real  input   value at time 0
c        val1      real  input   value at time 1
c        reltpos   real  input   relative time (between 0 and 1)
c        misdat    real  input   missing data flag (on if misdat<>0)
 
      implicit none
 
c     Declaration of parameters
      real      val0
      real      val1
      real      reltpos
      real      misdat
 
c     Numerical epsilon
      real      eps
      parameter (eps=1.e-8)
 
c     Local variables
      real      val
      real      mdv
 
c     Activate or inactive the missing data check (quick and dirty)
      if (misdat.ne.0.) then
         mdv = misdat
      else
         mdv = 257.22725394015
      endif
 
c     Do the linear interpolation
      if ( abs(reltpos).lt.eps ) then
         val = val0
      elseif ( abs(reltpos-1.).lt.eps ) then
         val = val1
      elseif ( (abs(val0-mdv).gt.eps).and.
     >         (abs(val1-mdv).gt.eps) ) then
         val = (1.-reltpos)*val0+reltpos*val1
      else
         val = mdv
      endif
 
c     Return value
      int_time = val
 
      end
 
 
c     -------------------------------------------------------------
c     Get the position of a physical point in grid space
c     -------------------------------------------------------------
 
      subroutine get_index3 (rid,rjd,rkd,xpo,ypo,ppo,mode,
     >                       vert,surf,nx,ny,nz,lonw,lats,dlon,dlat)
 
c     Purpose:
c        This subroutine determines the indices (rid,rjd,rkd) in grid 
c        space for a point in physical space (xpo,ypo,ppo). The 
c        horizontal grid is specified by the south-west point (lats,lonw)
c        and the grid spacing (dlat,dlon). The vertical grid is given
c        by <vert(n1,n2,n3)>. The lower boundary (typicall surface 
c        pressure) is given by <surf(n1,n2)>.
c     Arguments:
c        rid,rjd,rkd  real  output   grid location to be interpolated to
c        xpo,ypo,ppo  real  input    physical coordinates
c        n1,n2,n3     int   input    grid dimensions in x-, y- and p-direction
c        lats,lonw    real  input    south and west boundary of grid space
c        vert         real  input    vertical coordinate grid
c        surf         real  input    lower boundary (surface pressure)
c        mode         int   input    direction of vertical axis 
c                                        1: linear, 1 -> nz (th)
c                                        2: linear, nz -> 1 (pv)
c                                        3: binary (p)
 
 
      implicit none
 
c     Declartion of function parameters
      integer   nx,ny,nz
      real      vert(nx*ny*nz)
      real      surf(nx*ny)
      real      rid,rjd,rkd
      real      xpo,ypo,ppo
      real      dlat,dlon,lats,lonw
      integer   mode
 
c     Numerical epsilon
      real      eps
      parameter (eps=1.e-8)
 
c     Local variables
      integer   i,j,k
      real      ppo0,ppo1,ppom,psur
      integer   i0,im,i1
      
c     Externals 
      real      int_index3
      external  int_index3
 
c     Get the horizontal grid indices
      rid=(xpo-lonw)/dlon+1.
      rjd=(ypo-lats)/dlat+1.
 
c     Two-dimensional interpolation on horizontal plane: return
      if ( nz.eq.1 ) then
         rkd = 1.
         goto 100
      endif
         
c     Lowest-level interpolation: return
      if ( abs(ppo-1050.).lt.eps ) then
         rkd = 1.
         goto 100
      endif
 
c     Get the pressure at the lowest level and at the surface
      ppo0 = int_index3(vert,nx,ny,nz,rid,rjd,real(1),0.)
      psur = int_index3(surf,nx,ny, 1,rid,rjd,real(1),0.)
 
c     Decide whether the pressure is ascending or descending
      ppo1 = int_index3(vert,nx,ny,nz,rid,rjd,real(nz),0.)
      if ( ppo1.gt.ppo0 ) then
         ppo0 = ppo1
      endif
 
c     The point is between the surface and the lowest level: return
      if ( mode.eq.3 ) then
        if ( (ppo.ge.ppo0).and.(ppo.le.psur).or.
     >       (ppo.le.ppo0).and.(ppo.ge.psur) )
     >  then
           psur = int_index3(surf,nx,ny, 1,rid,rjd,real(1),0.)
           rkd  = (psur-ppo)/(psur-ppo0)
           goto 100
        endif
      endif
 
c     Full-level search (TH): linear ascending scanning through all levels
      if ( mode.eq.1 ) then
        
         ppo0 = int_index3(vert,nx,ny,nz,rid,rjd,real(1),0.)
         rkd=0
         do i=1,nz-1
            ppo1 = int_index3(vert,nx,ny,nz,rid,rjd,real(i+1),0.)
            if ( (ppo.gt.ppo0).and.(ppo.le.ppo1) ) then
               rkd=real(i)+(ppo-ppo0)/(ppo1-ppo0)
               goto 100
            endif
            ppo0 = ppo1
         enddo
 
c     Full-level search (PV): linear descending scanning through all levels
      elseif ( mode.eq.2 ) then
 
         ppo1 = int_index3(vert,nx,ny,nz,rid,rjd,real(nz),0.)
         rkd=0
         do i=nz-1,1,-1
            ppo0 = int_index3(vert,nx,ny,nz,rid,rjd,real(i),0.)
            if ( (ppo1.gt.ppo).and.(ppo0.le.ppo) ) then
               rkd=real(i)+(ppo0-ppo)/(ppo0-ppo1)
               goto 100
            endif
            ppo1 = ppo0
         enddo
 
c     Full-level search (P):  binary search 
      elseif ( mode.eq.3 ) then
 
         rkd  = 0
         i0   = 1
         i1   = nz
         ppo0 = int_index3(vert,nx,ny,nz,rid,rjd,real( 1),0.)
         ppo1 = int_index3(vert,nx,ny,nz,rid,rjd,real(nz),0.)
         
         do while ( i1.gt.(i0+1) )
            im   = (i0+i1)/2
            ppom = int_index3(vert,nx,ny,nz,rid,rjd,real(im),0.)
            if (ppom.lt.ppo) then
               i1   = im
               ppo1 = ppom
            else
               i0   = im
               ppo0 = ppom
            endif
         enddo
            
         rkd=real(i0)+(ppo0-ppo)/(ppo0-ppo1)
 
      endif
 
c     Exit point for subroutine
 100  continue
 
      end
 

Rev 13	Rev 15
1	`c *************************************************************`	1	`c *************************************************************`
2	`c * This package provides a general interpolaton routine *`	2	`c * This package provides a general interpolaton routine *`
3	`c *************************************************************`	3	`c *************************************************************`
4		4
5	`c The main interface routines are:`	5	`c The main interface routines are:`
6	`c get_index3,4 : get the grid indices for interpolation`	6	`c get_index3,4 : get the grid indices for interpolation`
7	`c int_index3,4 : interpolate to the grid position`	7	`c int_index3,4 : interpolate to the grid position`
8		8
9	`c -------------------------------------------------------------`	9	`c -------------------------------------------------------------`
10	`c Get index in grid space for interpolation`	10	`c Get index in grid space for interpolation`
11	`c -------------------------------------------------------------`	11	`c -------------------------------------------------------------`
12		12
13	`subroutine get_index4 (rid,rjd,rkd,xpo,ypo,ppo,rtp,`	13	`subroutine get_index4 (rid,rjd,rkd,xpo,ypo,ppo,rtp,`
14	`> vert0,vert1,surf0,surf1,mode,`	14	`> vert0,vert1,surf0,surf1,mode,`
15	`> nx,ny,nz,lonw,lats,dlon,dlat,misdat)`	15	`> nx,ny,nz,lonw,lats,dlon,dlat,misdat)`
16		16
17	`c Purpose:`	17	`c Purpose:`
18	`c This subroutine determines the indices (rid,rjd,rkd) in grid`	18	`c This subroutine determines the indices (rid,rjd,rkd) in grid`
19	`c space for a point in physical space (xpo,ypo,ppo). The`	19	`c space for a point in physical space (xpo,ypo,ppo). The`
20	`c horizontal grid is specified by the south-west point (lats,lonw)`	20	`c horizontal grid is specified by the south-west point (lats,lonw)`
21	`c and the grid spacing (dlat,dlon). The vertical grid is given`	21	`c and the grid spacing (dlat,dlon). The vertical grid is given`
22	`c by <vert(n1,n2,n3)>. The lower boundary (typicall surface`	22	`c by <vert(n1,n2,n3)>. The lower boundary (typicall surface`
23	`c pressure) is given by <surf(n1,n2)>.`	23	`c pressure) is given by <surf(n1,n2)>.`
24	`c Arguments:`	24	`c Arguments:`
25	`c rid,rjd,rkd real output grid location to be interpolated to`	25	`c rid,rjd,rkd real output grid location to be interpolated to`
26	`c xpo,ypo,ppo real input physical coordinates`	26	`c xpo,ypo,ppo real input physical coordinates`
27	`c rtp real input relative time position (0=beginning, 1=end)`	27	`c rtp real input relative time position (0=beginning, 1=end)`
28	`c n1,n2,n3 int input grid dimensions in x-, y- and p-direction`	28	`c n1,n2,n3 int input grid dimensions in x-, y- and p-direction`
29	`c lats,lonw real input south and west boundary of grid space`	29	`c lats,lonw real input south and west boundary of grid space`
30	`c vert real input vertical coordinate grid`	30	`c vert real input vertical coordinate grid`
31	`c surf real input lower boundary (surface pressure)`	31	`c surf real input lower boundary (surface pressure)`
32	`c mode int input direction of vertical axis (p=1,th=-1)`	32	`c mode int input direction of vertical axis (p=1,th=-1)`
33	`c 1: linear, 1 -> nz (th)`	33	`c 1: linear, 1 -> nz (th)`
34	`c 2: linear, nz -> 1 (pv)`	34	`c 2: linear, nz -> 1 (pv)`
35	`c 3: binary (p)`	35	`c 3: binary (p)`
36		36
37	`implicit none`	37	`implicit none`
38		38
39	`c Declartion of function parameters`	39	`c Declartion of function parameters`
40	`integer nx,ny,nz`	40	`integer nx,ny,nz`
41	`real xpo,ypo,ppo,rtp`	41	`real xpo,ypo,ppo,rtp`
42	`real vert0(nxnynz),vert1(nxnynz)`	42	`real vert0(nxnynz),vert1(nxnynz)`
43	`real surf0(nxny) ,surf1(nxny*nz)`	43	`real surf0(nxny) ,surf1(nxny*nz)`
44	`real rid,rjd,rkd`	44	`real rid,rjd,rkd`
45	`real dlat,dlon,lats,lonw`	45	`real dlat,dlon,lats,lonw`
46	`real misdat`	46	`real misdat`
47	`integer mode`	47	`integer mode`
48		48
49	`c Set numerical parameters`	49	`c Set numerical parameters`
50	`real eps`	50	`real eps`
51	`parameter (eps=1.e-8)`	51	`parameter (eps=1.e-8)`
52		52
53	`c Auxiliary variables`	53	`c Auxiliary variables`
54	`real rid0,rjd0,rkd0,rid1,rjd1,rkd1`	54	`real rid0,rjd0,rkd0,rid1,rjd1,rkd1`
55	`integer i`	55	`integer i`
56		56
57	`c Externals`	57	`c Externals`
58	`real int_time`	58	`real int_time`
59	`external int_time`	59	`external int_time`
60		60
61	`c Get the inidices`	61	`c Get the inidices`
62	`if (abs(rtp).lt.eps) then`	62	`if (abs(rtp).lt.eps) then`
63	`call get_index3 (rid,rjd,rkd,xpo,ypo,ppo,mode,`	63	`call get_index3 (rid,rjd,rkd,xpo,ypo,ppo,mode,`
64	`> vert0,surf0,nx,ny,nz,lonw,lats,dlon,dlat)`	64	`> vert0,surf0,nx,ny,nz,lonw,lats,dlon,dlat)`
65	`elseif (abs(rtp-1.).lt.eps) then`	65	`elseif (abs(rtp-1.).lt.eps) then`
66	`call get_index3 (rid,rjd,rkd,xpo,ypo,ppo,mode,`	66	`call get_index3 (rid,rjd,rkd,xpo,ypo,ppo,mode,`
67	`> vert1,surf1,nx,ny,nz,lonw,lats,dlon,dlat)`	67	`> vert1,surf1,nx,ny,nz,lonw,lats,dlon,dlat)`
68	`else`	68	`else`
69	`call get_index3 (rid0,rjd0,rkd0,xpo,ypo,ppo,mode,`	69	`call get_index3 (rid0,rjd0,rkd0,xpo,ypo,ppo,mode,`
70	`> vert0,surf0,nx,ny,nz,lonw,lats,dlon,dlat)`	70	`> vert0,surf0,nx,ny,nz,lonw,lats,dlon,dlat)`
71	`call get_index3 (rid1,rjd1,rkd1,xpo,ypo,ppo,mode,`	71	`call get_index3 (rid1,rjd1,rkd1,xpo,ypo,ppo,mode,`
72	`> vert1,surf1,nx,ny,nz,lonw,lats,dlon,dlat)`	72	`> vert1,surf1,nx,ny,nz,lonw,lats,dlon,dlat)`
73	`rid = int_time (rid0,rid1,rtp,misdat)`	73	`rid = int_time (rid0,rid1,rtp,misdat)`
74	`rjd = int_time (rjd0,rjd1,rtp,misdat)`	74	`rjd = int_time (rjd0,rjd1,rtp,misdat)`
75	`rkd = int_time (rkd0,rkd1,rtp,misdat)`	75	`rkd = int_time (rkd0,rkd1,rtp,misdat)`
76		76
77	`endif`	77	`endif`
78		78
79	`end`	79	`end`
80		80
81	`c -------------------------------------------------------------`	81	`c -------------------------------------------------------------`
82	`c Interpolate to an arbitrary position in grid space and time`	82	`c Interpolate to an arbitrary position in grid space and time`
83	`c -------------------------------------------------------------`	83	`c -------------------------------------------------------------`
84		84
85	`real function int_index4 (ar0,ar1,n1,n2,n3,rid,rjd,rkd,rtp,misdat)`	85	`real function int_index4 (ar0,ar1,n1,n2,n3,rid,rjd,rkd,rtp,misdat)`
86		86
87	`c Purpose:`	87	`c Purpose:`
88	`c This subroutine interpolates a 3d-array to an arbitrary`	88	`c This subroutine interpolates a 3d-array to an arbitrary`
89	`c location within the grid including a linear time-interpolation.`	89	`c location within the grid including a linear time-interpolation.`
90	`c Arguments:`	90	`c Arguments:`
91	`c rid,rjd,rkd real output grid location to be interpolated to`	91	`c rid,rjd,rkd real output grid location to be interpolated to`
92	`c xpo,ypo,ppo real input physical coordinates`	92	`c xpo,ypo,ppo real input physical coordinates`
93	`c n1,n2,n3 int input grid dimensions in x-, y- and p-direction`	93	`c n1,n2,n3 int input grid dimensions in x-, y- and p-direction`
94	`c lats,lonw real input south and west boundary of grid space`	94	`c lats,lonw real input south and west boundary of grid space`
95	`c vert real input vertical coordinate grid`	95	`c vert real input vertical coordinate grid`
96	`c surf real input lower boundary (surface pressure)`	96	`c surf real input lower boundary (surface pressure)`
97		97
98	`implicit none`	98	`implicit none`
99		99
100	`c Declartion of function parameters`	100	`c Declartion of function parameters`
101	`integer n1,n2,n3`	101	`integer n1,n2,n3`
102	`real ar0(n1n2n3),ar1(n1n2n3)`	102	`real ar0(n1n2n3),ar1(n1n2n3)`
103	`real rid,rjd,rkd`	103	`real rid,rjd,rkd`
104	`real rtp`	104	`real rtp`
105	`real misdat`	105	`real misdat`
106		106
107	`c Set numerical parameters`	107	`c Set numerical parameters`
108	`real eps`	108	`real eps`
109	`parameter (eps=1.e-8)`	109	`parameter (eps=1.e-8)`
110		110
111	`c Externals`	111	`c Externals`
112	`real int_index3,int_time`	112	`real int_index3,int_time`
113	`external int_index3,int_time`	113	`external int_index3,int_time`
114		114
115	`c Auxiliary variables`	115	`c Auxiliary variables`
116	`real val0,val1,val`	116	`real val0,val1,val`
117		117
118	`c Do the 3d-interpolation`	118	`c Do the 3d-interpolation`
119	`if (abs(rtp).lt.eps) then`	119	`if (abs(rtp).lt.eps) then`
120	`val = int_index3 (ar0,n1,n2,n3,rid,rjd,rkd,misdat)`	120	`val = int_index3 (ar0,n1,n2,n3,rid,rjd,rkd,misdat)`
121	`elseif (abs(rtp-1.).lt.eps) then`	121	`elseif (abs(rtp-1.).lt.eps) then`
122	`val = int_index3 (ar1,n1,n2,n3,rid,rjd,rkd,misdat)`	122	`val = int_index3 (ar1,n1,n2,n3,rid,rjd,rkd,misdat)`
123	`else`	123	`else`
124	`val0 = int_index3 (ar0,n1,n2,n3,rid,rjd,rkd,misdat)`	124	`val0 = int_index3 (ar0,n1,n2,n3,rid,rjd,rkd,misdat)`
125	`val1 = int_index3 (ar1,n1,n2,n3,rid,rjd,rkd,misdat)`	125	`val1 = int_index3 (ar1,n1,n2,n3,rid,rjd,rkd,misdat)`
126	`val = int_time (val0,val1,rtp,misdat)`	126	`val = int_time (val0,val1,rtp,misdat)`
127	`endif`	127	`endif`
128		128
129	`c Return value`	129	`c Return value`
130	`int_index4 = val`	130	`int_index4 = val`
131		131
132	`return`	132	`return`
133	`end`	133	`end`
134		134
135		135
136	`c -------------------------------------------------------------`	136	`c -------------------------------------------------------------`
137	`c Interpolate to an arbitrary position in grid space`	137	`c Interpolate to an arbitrary position in grid space`
138	`c -------------------------------------------------------------`	138	`c -------------------------------------------------------------`
139		139
140	`real function int_index3 (ar,n1,n2,n3,rid,rjd,rkd,misdat)`	140	`real function int_index3 (ar,n1,n2,n3,rid,rjd,rkd,misdat)`
141		141
142	`c Purpose:`	142	`c Purpose:`
143	`c This subroutine interpolates a 3d-array to an arbitrary`	143	`c This subroutine interpolates a 3d-array to an arbitrary`
144	`c location within the grid. The interpolation includes the`	144	`c location within the grid. The interpolation includes the`
145	`c testing of the missing data flag 'misdat'. If one dimension`	145	`c testing of the missing data flag 'misdat'. If one dimension`
146	`c is 1, a 2d-interpolation is performed; if two dimensions`	146	`c is 1, a 2d-interpolation is performed; if two dimensions`
147	`c are 1, it is a 1d-interpolation; if all three dimensions are`	147	`c are 1, it is a 1d-interpolation; if all three dimensions are`
148	`c 1, no interpolation is performed and the input value is`	148	`c 1, no interpolation is performed and the input value is`
149	`c returned.`	149	`c returned.`
150	`c Arguments:`	150	`c Arguments:`
151	`c ar real input input data array`	151	`c ar real input input data array`
152	`c n1,n2,n3 int input dimensions of ar`	152	`c n1,n2,n3 int input dimensions of ar`
153	`c ri,rj,rk real input grid location to be interpolated to`	153	`c ri,rj,rk real input grid location to be interpolated to`
154	`c misdat real input missing data flag (on if misdat<>0)`	154	`c misdat real input missing data flag (on if misdat<>0)`
155		155
156	`implicit none`	156	`implicit none`
157		157
158	`c Declartion of function parameters`	158	`c Declartion of function parameters`
159	`integer n1,n2,n3`	159	`integer n1,n2,n3`
160	`real ar(n1n2n3)`	160	`real ar(n1n2n3)`
161	`real rid,rjd,rkd`	161	`real rid,rjd,rkd`
162	`real misdat`	162	`real misdat`
163		163
164	`c Set numerical parameters`	164	`c Set numerical parameters`
165	`real eps`	165	`real eps`
166	`parameter (eps=1.e-8)`	166	`parameter (eps=1.e-8)`
167		167
168	`c Local variables`	168	`c Local variables`
169	`integer i,j,k,ip1,jp1,kp1`	169	`integer i,j,k,ip1,jp1,kp1`
170	`real frac0i,frac0j,frac0k,frac1i,frac1j,frac1k`	170	`real frac0i,frac0j,frac0k,frac1i,frac1j,frac1k`
171	`real ri,rj,rk`	171	`real ri,rj,rk`
172	`real val000,val001,val010,val011,val100,val101,val110,val111`	172	`real val000,val001,val010,val011,val100,val101,val110,val111`
173	`real frc000,frc001,frc010,frc011,frc100,frc101,frc110,frc111`	173	`real frc000,frc001,frc010,frc011,frc100,frc101,frc110,frc111`
174	`real frc`	174	`real frc`
175	`real mdv`	175	`real mdv`
176	`real val`	176	`real val`
177		177
178	`c Elementary test for dimensions`	178	`c Elementary test for dimensions`
179	`if ( (n1.lt.1).or.(n2.lt.1).or.(n3.lt.1) ) then`	179	`if ( (n1.lt.1).or.(n2.lt.1).or.(n3.lt.1) ) then`
180	`print*,'Invalid grid dimensions ',n1,n2,n3`	180	`print*,'Invalid grid dimensions ',n1,n2,n3`
181	`stop`	181	`stop`
182	`endif`	182	`endif`
183		183
184	`c Activate or inactive the missing data check (quick and dirty)`	184	`c Activate or inactive the missing data check (quick and dirty)`
185	`if (misdat.ne.0.) then`	185	`if (misdat.ne.0.) then`
186	`mdv = misdat`	186	`mdv = misdat`
187	`else`	187	`else`
188	`mdv = 257.22725394015`	188	`mdv = 257.22725394015`
189	`endif`	189	`endif`
190		190
191	`c Bring the indices into the grid space`	191	`c Bring the indices into the grid space`
192	`ri = amax1(1.,amin1(float(n1),rid))`	192	`ri = amax1(1.,amin1(float(n1),rid))`
193	`rj = amax1(1.,amin1(float(n2),rjd))`	193	`rj = amax1(1.,amin1(float(n2),rjd))`
194	`rk = amax1(1.,amin1(float(n3),rkd))`	194	`rk = amax1(1.,amin1(float(n3),rkd))`
195		195
196	`c Get the index of the west-south-bottom corner of the box`	196	`c Get the index of the west-south-bottom corner of the box`
197	`i = min0(int(ri),n1-1)`	197	`i = min0(int(ri),n1-1)`
198	`ip1 = i+1`	198	`ip1 = i+1`
199	`j = min0(int(rj),n2-1)`	199	`j = min0(int(rj),n2-1)`
200	`jp1 = j+1`	200	`jp1 = j+1`
201	`k = min0(int(rk),n3-1)`	201	`k = min0(int(rk),n3-1)`
202	`kp1 = k+1`	202	`kp1 = k+1`
203		203
204	`c Special handling for 2d arrays`	204	`c Special handling for 2d arrays`
205	`if (n3.eq.1) then`	205	`if (n3.eq.1) then`
206	`k=1`	206	`k=1`
207	`kp1=1`	207	`kp1=1`
208	`endif`	208	`endif`
209		209
210	`c Get location relative to grid box`	210	`c Get location relative to grid box`
211	`if ( i.ne.ip1 ) then`	211	`if ( i.ne.ip1 ) then`
212	`frac0i = ri-float(i)`	212	`frac0i = ri-float(i)`
213	`frac1i = 1.-frac0i`	213	`frac1i = 1.-frac0i`
214	`else`	214	`else`
215	`frac0i = 0.`	215	`frac0i = 0.`
216	`frac1i = 1.`	216	`frac1i = 1.`
217	`endif`	217	`endif`
218	`if ( j.ne.jp1 ) then`	218	`if ( j.ne.jp1 ) then`
219	`frac0j = rj-float(j)`	219	`frac0j = rj-float(j)`
220	`frac1j = 1.-frac0j`	220	`frac1j = 1.-frac0j`
221	`else`	221	`else`
222	`frac0j = 0.`	222	`frac0j = 0.`
223	`frac1j = 1.`	223	`frac1j = 1.`
224	`endif`	224	`endif`
225	`if ( k.ne.kp1 ) then`	225	`if ( k.ne.kp1 ) then`
226	`frac0k = rk-float(k)`	226	`frac0k = rk-float(k)`
227	`frac1k = 1.-frac0k`	227	`frac1k = 1.-frac0k`
228	`else`	228	`else`
229	`frac0k = 0.`	229	`frac0k = 0.`
230	`frac1k = 1.`	230	`frac1k = 1.`
231	`endif`	231	`endif`
232		232
233	`c On a grid point - take the grid point value`	233	`c On a grid point - take the grid point value`
234	`if ( ( abs(frac0i).lt.eps ).and.`	234	`if ( ( abs(frac0i).lt.eps ).and.`
235	`> ( abs(frac0j).lt.eps ).and.`	235	`> ( abs(frac0j).lt.eps ).and.`
236	`> ( abs(frac0k).lt.eps ) ) then`	236	`> ( abs(frac0k).lt.eps ) ) then`
237		237
238	`val = ar( i + n1(j -1) + n1n2*(k -1) )`	238	`val = ar( i + n1(j -1) + n1n2*(k -1) )`
239	`goto 100`	239	`goto 100`
240		240
241	`endif`	241	`endif`
242		242
243	`c Init the fractions`	243	`c Init the fractions`
244	`frc000 = frac1i * frac1j * frac1k`	244	`frc000 = frac1i * frac1j * frac1k`
245	`frc001 = frac0i * frac1j * frac1k`	245	`frc001 = frac0i * frac1j * frac1k`
246	`frc010 = frac1i * frac0j * frac1k`	246	`frc010 = frac1i * frac0j * frac1k`
247	`frc011 = frac0i * frac0j * frac1k`	247	`frc011 = frac0i * frac0j * frac1k`
248	`frc100 = frac1i * frac1j * frac0k`	248	`frc100 = frac1i * frac1j * frac0k`
249	`frc101 = frac0i * frac1j * frac0k`	249	`frc101 = frac0i * frac1j * frac0k`
250	`frc110 = frac1i * frac0j * frac0k`	250	`frc110 = frac1i * frac0j * frac0k`
251	`frc111 = frac0i * frac0j * frac0k`	251	`frc111 = frac0i * frac0j * frac0k`
252		252
253	`c Init the values`	253	`c Init the values`
254	`val000 = ar( i + n1(j -1) + n1n2*(k -1) )`	254	`val000 = ar( i + n1(j -1) + n1n2*(k -1) )`
255	`val001 = ar( ip1 + n1(j -1) + n1n2*(k -1) )`	255	`val001 = ar( ip1 + n1(j -1) + n1n2*(k -1) )`
256	`val010 = ar( i + n1(jp1-1) + n1n2*(k -1) )`	256	`val010 = ar( i + n1(jp1-1) + n1n2*(k -1) )`
257	`val011 = ar( ip1 + n1(jp1-1) + n1n2*(k -1) )`	257	`val011 = ar( ip1 + n1(jp1-1) + n1n2*(k -1) )`
258	`val100 = ar( i + n1(j -1) + n1n2*(kp1-1) )`	258	`val100 = ar( i + n1(j -1) + n1n2*(kp1-1) )`
259	`val101 = ar( ip1 + n1(j -1) + n1n2*(kp1-1) )`	259	`val101 = ar( ip1 + n1(j -1) + n1n2*(kp1-1) )`
260	`val110 = ar( i + n1(jp1-1) + n1n2*(kp1-1) )`	260	`val110 = ar( i + n1(jp1-1) + n1n2*(kp1-1) )`
261	`val111 = ar( ip1 + n1(jp1-1) + n1n2*(kp1-1) )`	261	`val111 = ar( ip1 + n1(jp1-1) + n1n2*(kp1-1) )`
262		262
263	`c Handle missing data`	263	`c Handle missing data`
264	`if ( abs(val000-mdv).lt.eps ) frc000 = 0.`	264	`if ( abs(val000-mdv).lt.eps ) frc000 = 0.`
265	`if ( abs(val001-mdv).lt.eps ) frc001 = 0.`	265	`if ( abs(val001-mdv).lt.eps ) frc001 = 0.`
266	`if ( abs(val010-mdv).lt.eps ) frc010 = 0.`	266	`if ( abs(val010-mdv).lt.eps ) frc010 = 0.`
267	`if ( abs(val011-mdv).lt.eps ) frc011 = 0.`	267	`if ( abs(val011-mdv).lt.eps ) frc011 = 0.`
268	`if ( abs(val100-mdv).lt.eps ) frc100 = 0.`	268	`if ( abs(val100-mdv).lt.eps ) frc100 = 0.`
269	`if ( abs(val101-mdv).lt.eps ) frc101 = 0.`	269	`if ( abs(val101-mdv).lt.eps ) frc101 = 0.`
270	`if ( abs(val110-mdv).lt.eps ) frc110 = 0.`	270	`if ( abs(val110-mdv).lt.eps ) frc110 = 0.`
271	`if ( abs(val111-mdv).lt.eps ) frc111 = 0.`	271	`if ( abs(val111-mdv).lt.eps ) frc111 = 0.`
272		272
273	`c Build the final value`	273	`c Build the final value`
274	`frc = frc000 + frc001 + frc010 + frc011 +`	274	`frc = frc000 + frc001 + frc010 + frc011 +`
275	`> frc100 + frc101 + frc110 + frc111`	275	`> frc100 + frc101 + frc110 + frc111`
276	`if ( frc.gt.0. ) then`	276	`if ( frc.gt.0. ) then`
277	`val = 1./frc * ( frc000 * val000 + frc001 * val001 +`	277	`val = 1./frc * ( frc000 * val000 + frc001 * val001 +`
278	`> frc010 * val010 + frc011 * val011 +`	278	`> frc010 * val010 + frc011 * val011 +`
279	`> frc100 * val100 + frc101 * val101 +`	279	`> frc100 * val100 + frc101 * val101 +`
280	`> frc110 * val110 + frc111 * val111 )`	280	`> frc110 * val110 + frc111 * val111 )`
281	`else`	281	`else`
282	`val = misdat`	282	`val = misdat`
283	`endif`	283	`endif`
284		284
285	`c Return the value`	285	`c Return the value`
286	`100 continue`	286	`100 continue`
287		287
288	`int_index3 = val`	288	`int_index3 = val`
289		289
290	`end`	290	`end`
291		291
292		292
293	`c -------------------------------------------------------------`	293	`c -------------------------------------------------------------`
294	`c Time interpolation`	294	`c Time interpolation`
295	`c -------------------------------------------------------------`	295	`c -------------------------------------------------------------`
296		296
297	`real function int_time (val0,val1,reltpos,misdat)`	297	`real function int_time (val0,val1,reltpos,misdat)`
298		298
299	`c Purpose:`	299	`c Purpose:`
300	`c This subroutine interpolates linearly in time between two`	300	`c This subroutine interpolates linearly in time between two`
301	`c values.`	301	`c values.`
302	`c Arguments:`	302	`c Arguments:`
303	`c val0 real input value at time 0`	303	`c val0 real input value at time 0`
304	`c val1 real input value at time 1`	304	`c val1 real input value at time 1`
305	`c reltpos real input relative time (between 0 and 1)`	305	`c reltpos real input relative time (between 0 and 1)`
306	`c misdat real input missing data flag (on if misdat<>0)`	306	`c misdat real input missing data flag (on if misdat<>0)`
307		307
308	`implicit none`	308	`implicit none`
309		309
310	`c Declaration of parameters`	310	`c Declaration of parameters`
311	`real val0`	311	`real val0`
312	`real val1`	312	`real val1`
313	`real reltpos`	313	`real reltpos`
314	`real misdat`	314	`real misdat`
315		315
316	`c Numerical epsilon`	316	`c Numerical epsilon`
317	`real eps`	317	`real eps`
318	`parameter (eps=1.e-8)`	318	`parameter (eps=1.e-8)`
319		319
320	`c Local variables`	320	`c Local variables`
321	`real val`	321	`real val`
322	`real mdv`	322	`real mdv`
323		323
324	`c Activate or inactive the missing data check (quick and dirty)`	324	`c Activate or inactive the missing data check (quick and dirty)`
325	`if (misdat.ne.0.) then`	325	`if (misdat.ne.0.) then`
326	`mdv = misdat`	326	`mdv = misdat`
327	`else`	327	`else`
328	`mdv = 257.22725394015`	328	`mdv = 257.22725394015`
329	`endif`	329	`endif`
330		330
331	`c Do the linear interpolation`	331	`c Do the linear interpolation`
332	`if ( abs(reltpos).lt.eps ) then`	332	`if ( abs(reltpos).lt.eps ) then`
333	`val = val0`	333	`val = val0`
334	`elseif ( abs(reltpos-1.).lt.eps ) then`	334	`elseif ( abs(reltpos-1.).lt.eps ) then`
335	`val = val1`	335	`val = val1`
336	`elseif ( (abs(val0-mdv).gt.eps).and.`	336	`elseif ( (abs(val0-mdv).gt.eps).and.`
337	`> (abs(val1-mdv).gt.eps) ) then`	337	`> (abs(val1-mdv).gt.eps) ) then`
338	`val = (1.-reltpos)val0+reltposval1`	338	`val = (1.-reltpos)val0+reltposval1`
339	`else`	339	`else`
340	`val = mdv`	340	`val = mdv`
341	`endif`	341	`endif`
342		342
343	`c Return value`	343	`c Return value`
344	`int_time = val`	344	`int_time = val`
345		345
346	`end`	346	`end`
347		347
348		348
349	`c -------------------------------------------------------------`	349	`c -------------------------------------------------------------`
350	`c Get the position of a physical point in grid space`	350	`c Get the position of a physical point in grid space`
351	`c -------------------------------------------------------------`	351	`c -------------------------------------------------------------`
352		352
353	`subroutine get_index3 (rid,rjd,rkd,xpo,ypo,ppo,mode,`	353	`subroutine get_index3 (rid,rjd,rkd,xpo,ypo,ppo,mode,`
354	`> vert,surf,nx,ny,nz,lonw,lats,dlon,dlat)`	354	`> vert,surf,nx,ny,nz,lonw,lats,dlon,dlat)`
355		355
356	`c Purpose:`	356	`c Purpose:`
357	`c This subroutine determines the indices (rid,rjd,rkd) in grid`	357	`c This subroutine determines the indices (rid,rjd,rkd) in grid`
358	`c space for a point in physical space (xpo,ypo,ppo). The`	358	`c space for a point in physical space (xpo,ypo,ppo). The`
359	`c horizontal grid is specified by the south-west point (lats,lonw)`	359	`c horizontal grid is specified by the south-west point (lats,lonw)`
360	`c and the grid spacing (dlat,dlon). The vertical grid is given`	360	`c and the grid spacing (dlat,dlon). The vertical grid is given`
361	`c by <vert(n1,n2,n3)>. The lower boundary (typicall surface`	361	`c by <vert(n1,n2,n3)>. The lower boundary (typicall surface`
362	`c pressure) is given by <surf(n1,n2)>.`	362	`c pressure) is given by <surf(n1,n2)>.`
363	`c Arguments:`	363	`c Arguments:`
364	`c rid,rjd,rkd real output grid location to be interpolated to`	364	`c rid,rjd,rkd real output grid location to be interpolated to`
365	`c xpo,ypo,ppo real input physical coordinates`	365	`c xpo,ypo,ppo real input physical coordinates`
366	`c n1,n2,n3 int input grid dimensions in x-, y- and p-direction`	366	`c n1,n2,n3 int input grid dimensions in x-, y- and p-direction`
367	`c lats,lonw real input south and west boundary of grid space`	367	`c lats,lonw real input south and west boundary of grid space`
368	`c vert real input vertical coordinate grid`	368	`c vert real input vertical coordinate grid`
369	`c surf real input lower boundary (surface pressure)`	369	`c surf real input lower boundary (surface pressure)`
370	`c mode int input direction of vertical axis`	370	`c mode int input direction of vertical axis`
371	`c 1: linear, 1 -> nz (th)`	371	`c 1: linear, 1 -> nz (th)`
372	`c 2: linear, nz -> 1 (pv)`	372	`c 2: linear, nz -> 1 (pv)`
373	`c 3: binary (p)`	373	`c 3: binary (p)`
374		374
375		375
376	`implicit none`	376	`implicit none`
377		377
378	`c Declartion of function parameters`	378	`c Declartion of function parameters`
379	`integer nx,ny,nz`	379	`integer nx,ny,nz`
380	`real vert(nxnynz)`	380	`real vert(nxnynz)`
381	`real surf(nx*ny)`	381	`real surf(nx*ny)`
382	`real rid,rjd,rkd`	382	`real rid,rjd,rkd`
383	`real xpo,ypo,ppo`	383	`real xpo,ypo,ppo`
384	`real dlat,dlon,lats,lonw`	384	`real dlat,dlon,lats,lonw`
385	`integer mode`	385	`integer mode`
386		386
387	`c Numerical epsilon`	387	`c Numerical epsilon`
388	`real eps`	388	`real eps`
389	`parameter (eps=1.e-8)`	389	`parameter (eps=1.e-8)`
390		390
391	`c Local variables`	391	`c Local variables`
392	`integer i,j,k`	392	`integer i,j,k`
393	`real ppo0,ppo1,ppom,psur`	393	`real ppo0,ppo1,ppom,psur`
394	`integer i0,im,i1`	394	`integer i0,im,i1`
395		395
396	`c Externals`	396	`c Externals`
397	`real int_index3`	397	`real int_index3`
398	`external int_index3`	398	`external int_index3`
399		399
400	`c Get the horizontal grid indices`	400	`c Get the horizontal grid indices`
401	`rid=(xpo-lonw)/dlon+1.`	401	`rid=(xpo-lonw)/dlon+1.`
402	`rjd=(ypo-lats)/dlat+1.`	402	`rjd=(ypo-lats)/dlat+1.`
403		403
404	`c Two-dimensional interpolation on horizontal plane: return`	404	`c Two-dimensional interpolation on horizontal plane: return`
405	`if ( nz.eq.1 ) then`	405	`if ( nz.eq.1 ) then`
406	`rkd = 1.`	406	`rkd = 1.`
407	`goto 100`	407	`goto 100`
408	`endif`	408	`endif`
409		409
410	`c Lowest-level interpolation: return`	410	`c Lowest-level interpolation: return`
411	`if ( abs(ppo-1050.).lt.eps ) then`	411	`if ( abs(ppo-1050.).lt.eps ) then`
412	`rkd = 1.`	412	`rkd = 1.`
413	`goto 100`	413	`goto 100`
414	`endif`	414	`endif`
415		415
416	`c Get the pressure at the lowest level and at the surface`	416	`c Get the pressure at the lowest level and at the surface`
417	`ppo0 = int_index3(vert,nx,ny,nz,rid,rjd,real(1),0.)`	417	`ppo0 = int_index3(vert,nx,ny,nz,rid,rjd,real(1),0.)`
418	`psur = int_index3(surf,nx,ny, 1,rid,rjd,real(1),0.)`	418	`psur = int_index3(surf,nx,ny, 1,rid,rjd,real(1),0.)`
419		419
420	`c Decide whether the pressure is ascending or descending`	420	`c Decide whether the pressure is ascending or descending`
421	`ppo1 = int_index3(vert,nx,ny,nz,rid,rjd,real(nz),0.)`	421	`ppo1 = int_index3(vert,nx,ny,nz,rid,rjd,real(nz),0.)`
422	`if ( ppo1.gt.ppo0 ) then`	422	`if ( ppo1.gt.ppo0 ) then`
423	`ppo0 = ppo1`	423	`ppo0 = ppo1`
424	`endif`	424	`endif`
425		425
426	`c The point is between the surface and the lowest level: return`	426	`c The point is between the surface and the lowest level: return`
427	`if ( mode.eq.3 ) then`	427	`if ( mode.eq.3 ) then`
428	`if ( (ppo.ge.ppo0).and.(ppo.le.psur).or.`	428	`if ( (ppo.ge.ppo0).and.(ppo.le.psur).or.`
429	`> (ppo.le.ppo0).and.(ppo.ge.psur) )`	429	`> (ppo.le.ppo0).and.(ppo.ge.psur) )`
430	`> then`	430	`> then`
431	`psur = int_index3(surf,nx,ny, 1,rid,rjd,real(1),0.)`	431	`psur = int_index3(surf,nx,ny, 1,rid,rjd,real(1),0.)`
432	`rkd = (psur-ppo)/(psur-ppo0)`	432	`rkd = (psur-ppo)/(psur-ppo0)`
433	`goto 100`	433	`goto 100`
434	`endif`	434	`endif`
435	`endif`	435	`endif`
436		436
437	`c Full-level search (TH): linear ascending scanning through all levels`	437	`c Full-level search (TH): linear ascending scanning through all levels`
438	`if ( mode.eq.1 ) then`	438	`if ( mode.eq.1 ) then`
439		439
440	`ppo0 = int_index3(vert,nx,ny,nz,rid,rjd,real(1),0.)`	440	`ppo0 = int_index3(vert,nx,ny,nz,rid,rjd,real(1),0.)`
441	`rkd=0`	441	`rkd=0`
442	`do i=1,nz-1`	442	`do i=1,nz-1`
443	`ppo1 = int_index3(vert,nx,ny,nz,rid,rjd,real(i+1),0.)`	443	`ppo1 = int_index3(vert,nx,ny,nz,rid,rjd,real(i+1),0.)`
444	`if ( (ppo0.lt.ppo).and.(ppo1.ge.ppo) ) then`	444	`if ( (ppo.gt.ppo0).and.(ppo.le.ppo1) ) then`
445	`rkd=real(i)+(ppo0-ppo)/(ppo0-ppo1)`	445	`rkd=real(i)+(ppo-ppo0)/(ppo1-ppo0)`
446	`goto 100`	446	`goto 100`
447	`endif`	447	`endif`
448	`ppo0 = ppo1`	448	`ppo0 = ppo1`
449	`enddo`	449	`enddo`
450		450
451	`c Full-level search (PV): linear descending scanning through all levels`	451	`c Full-level search (PV): linear descending scanning through all levels`
452	`elseif ( mode.eq.2 ) then`	452	`elseif ( mode.eq.2 ) then`
453		453
454	`ppo1 = int_index3(vert,nx,ny,nz,rid,rjd,real(nz),0.)`	454	`ppo1 = int_index3(vert,nx,ny,nz,rid,rjd,real(nz),0.)`
455	`rkd=0`	455	`rkd=0`
456	`do i=nz-1,1,-1`	456	`do i=nz-1,1,-1`
457	`ppo0 = int_index3(vert,nx,ny,nz,rid,rjd,real(i),0.)`	457	`ppo0 = int_index3(vert,nx,ny,nz,rid,rjd,real(i),0.)`
458	`if ( (ppo1.gt.ppo).and.(ppo0.le.ppo) ) then`	458	`if ( (ppo1.gt.ppo).and.(ppo0.le.ppo) ) then`
459	`rkd=real(i)+(ppo0-ppo)/(ppo0-ppo1)`	459	`rkd=real(i)+(ppo0-ppo)/(ppo0-ppo1)`
460	`goto 100`	460	`goto 100`
461	`endif`	461	`endif`
462	`ppo1 = ppo0`	462	`ppo1 = ppo0`
463	`enddo`	463	`enddo`
464		464
465	`c Full-level search (P): binary search`	465	`c Full-level search (P): binary search`
466	`elseif ( mode.eq.3 ) then`	466	`elseif ( mode.eq.3 ) then`
467		467
468	`rkd = 0`	468	`rkd = 0`
469	`i0 = 1`	469	`i0 = 1`
470	`i1 = nz`	470	`i1 = nz`
471	`ppo0 = int_index3(vert,nx,ny,nz,rid,rjd,real( 1),0.)`	471	`ppo0 = int_index3(vert,nx,ny,nz,rid,rjd,real( 1),0.)`
472	`ppo1 = int_index3(vert,nx,ny,nz,rid,rjd,real(nz),0.)`	472	`ppo1 = int_index3(vert,nx,ny,nz,rid,rjd,real(nz),0.)`
473		473
474	`do while ( i1.gt.(i0+1) )`	474	`do while ( i1.gt.(i0+1) )`
475	`im = (i0+i1)/2`	475	`im = (i0+i1)/2`
476	`ppom = int_index3(vert,nx,ny,nz,rid,rjd,real(im),0.)`	476	`ppom = int_index3(vert,nx,ny,nz,rid,rjd,real(im),0.)`
477	`if (ppom.lt.ppo) then`	477	`if (ppom.lt.ppo) then`
478	`i1 = im`	478	`i1 = im`
479	`ppo1 = ppom`	479	`ppo1 = ppom`
480	`else`	480	`else`
481	`i0 = im`	481	`i0 = im`
482	`ppo0 = ppom`	482	`ppo0 = ppom`
483	`endif`	483	`endif`
484	`enddo`	484	`enddo`
485		485
486	`rkd=real(i0)+(ppo0-ppo)/(ppo0-ppo1)`	486	`rkd=real(i0)+(ppo0-ppo)/(ppo0-ppo1)`
487		487
488	`endif`	488	`endif`
489		489
490	`c Exit point for subroutine`	490	`c Exit point for subroutine`
491	`100 continue`	491	`100 continue`
492		492
493	`end`	493	`end`
494		494

Subversion Repositories lagranto.ecmwf

(root)/trunk/lib/inter.f – Rev 13 → 15