/[MITgcm]/MITgcm/pkg/exf/exf_interp.F
ViewVC logotype

Diff of /MITgcm/pkg/exf/exf_interp.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.13 by heimbach, Sat Jan 14 18:09:45 2006 UTC revision 1.18 by dimitri, Sun Jul 2 13:02:22 2006 UTC
# Line 52  C     lon_0, lat_0 = lon and lat of sw c Line 52  C     lon_0, lat_0 = lon and lat of sw c
52  C     lon_inc      = scalar x-grid increment  C     lon_inc      = scalar x-grid increment
53  C     lat_inc      = vector y-grid increments  C     lat_inc      = vector y-grid increments
54  C     nx_in, ny_in = input x-grid and y-grid size  C     nx_in, ny_in = input x-grid and y-grid size
55  C     method       = 1 for bilinear 2 for bicubic  C     method       = 1,11,21 for bilinear; 2,12,22 for bicubic
56    C                    1,2 for tracer; 11,12 for U; 21,22 for V
57  C     mythid       = thread id  C     mythid       = thread id
58  C  C
59    
# Line 71  C subroutine variables Line 72  C subroutine variables
72        integer       method, mythid        integer       method, mythid
73    
74  C local variables  C local variables
       real*8   ne_fac,nw_fac,se_fac,sw_fac  
75        integer  e_ind(snx,sny),w_ind(snx,sny)        integer  e_ind(snx,sny),w_ind(snx,sny)
76        integer  n_ind(snx,sny),s_ind(snx,sny)        integer  n_ind(snx,sny),s_ind(snx,sny)
77        real*8   px_ind(4), py_ind(4), ew_val(4)        real*8   px_ind(4), py_ind(4), ew_val(4)
# Line 79  C local variables Line 79  C local variables
79        real*8   lagran        real*8   lagran
80        real*4   arrayin(-1:nx_in+2 ,      -1:ny_in+2)        real*4   arrayin(-1:nx_in+2 ,      -1:ny_in+2)
81        real*8   x_in   (-1:nx_in+2), y_in(-1:ny_in+2)        real*8   x_in   (-1:nx_in+2), y_in(-1:ny_in+2)
82          real*8   ninety
83          PARAMETER ( ninety = 90. )
84        integer  i, j, k, l, js, bi, bj, sp, interp_unit        integer  i, j, k, l, js, bi, bj, sp, interp_unit
85        _RS      xG(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RS      xG(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
86        _RS      threeSixtyRS        _RS      threeSixtyRS, NorthValue
       _RL      threeSixtyRL  
87        PARAMETER ( threeSixtyRS = 360. )        PARAMETER ( threeSixtyRS = 360. )
       PARAMETER ( threeSixtyRL = 360. )  
88    
89  !     nomalize xG  C     put xG in interval [ lon_0 , lon_0+360 [
90        do bj=myByLo(myThid),myByHi(myThid)        do bj=myByLo(myThid),myByHi(myThid)
91         do bi=myBxLo(myThid),myBxHi(myThid)         do bi=myBxLo(myThid),myBxHi(myThid)
92          do j=1-OLy,sNy+OLy          do j=1-OLy,sNy+OLy
93           do i=1-OLx,sNx+OLx           do i=1-OLx,sNx+OLx
94            xG(i,j,bi,bj) = xG_in(i,j,bi,bj)            xG(i,j,bi,bj) = xG_in(i,j,bi,bj)-lon_0
95            xG(i,j,bi,bj) = mod(xG(i,j,bi,bj),threeSixtyRS)       &                  + threeSixtyRS*2.
96            if ( xG(i,j,bi,bj) .lt. 0. )            xG(i,j,bi,bj) = lon_0+mod(xG(i,j,bi,bj),threeSixtyRS)
      &     xG(i,j,bi,bj) = threeSixtyRS+xG(i,j,bi,bj)  
97           enddo           enddo
98          enddo          enddo
99         enddo         enddo
100        enddo        enddo
101    
102         call exf_interp_read(         call exf_interp_read(
103       I   infile,       I   infile, filePrec,
      I   filePrec,  
104       O   arrayin,       O   arrayin,
105       I   irecord, xG, yG,       I   irecord, nx_in, ny_in, mythid)
      I   lon_0, lon_inc,  
      I   lat_0, lat_inc,  
      I   nx_in, ny_in, method, mythid)  
106        _BARRIER        _BARRIER
107    
108  C     _BEGIN_MASTER( myThid )  C     _BEGIN_MASTER( myThid )
109    
110  C setup input grid  C setup input longitude grid
111         do i=-1,nx_in+2        do i=-1,nx_in+2
112          x_in(i) = lon_0 + (i-1.)*lon_inc         x_in(i) = lon_0 + (i-1)*lon_inc
113          x_in(i) = mod(x_in(i),threeSixtyRL)        enddo
         if ( x_in(i) .lt. 0. )  
      &   x_in(i) = threeSixtyRL+x_in(i)  
        enddo  
114    
115         y_in(0) = lat_0 - lat_inc(1)  C setup input latitude grid
116         y_in(-1)= lat_0 - 2.*lat_inc(1)        y_in(0) = lat_0 - lat_inc(1)
117         y_in(1) = lat_0        y_in(-1)= lat_0 - 2.*lat_inc(1)
118         do j=2,ny_in        y_in(1) = lat_0
119          y_in(j) = y_in(j-1) + lat_inc(j-1)        do j=2,ny_in
120         enddo         y_in(j) = y_in(j-1) + lat_inc(j-1)
121         y_in(ny_in+1) = y_in(ny_in) + lat_inc(ny_in-1)        enddo
122         y_in(ny_in+2) = y_in(ny_in) + 2.*lat_inc(ny_in-1)        do j=ny_in+1,ny_in+2
123           if (y_in(j-1).eq.ninety) then
124            y_in(j) = 2 * ninety - y_in(j-2)
125           else
126            y_in(j) = min( y_in(j-1)+lat_inc(ny_in-1), ninety )
127           endif
128          enddo
129    
130  C enlarge boundary  C enlarge boundary
131         do j=1,ny_in        do j=1,ny_in
132          arrayin(0,j)       = arrayin(nx_in,j)         arrayin(0,j)       = arrayin(nx_in,j)
133          arrayin(-1,j)      = arrayin(nx_in-1,j)         arrayin(-1,j)      = arrayin(nx_in-1,j)
134          arrayin(nx_in+1,j) = arrayin(1,j)         arrayin(nx_in+1,j) = arrayin(1,j)
135          arrayin(nx_in+2,j) = arrayin(2,j)         arrayin(nx_in+2,j) = arrayin(2,j)
136         enddo        enddo
137         do i=-1,nx_in+2        do i=-1,nx_in+2
138          arrayin(i,0)       = arrayin(i,1)         arrayin(i,0)       = arrayin(i,1)
139          arrayin(i,-1)      = arrayin(i,1)         arrayin(i,-1)      = arrayin(i,1)
140          arrayin(i,ny_in+1) = arrayin(i,ny_in)         arrayin(i,ny_in+1) = arrayin(i,ny_in)
141          arrayin(i,ny_in+2) = arrayin(i,ny_in)         arrayin(i,ny_in+2) = arrayin(i,ny_in)
142         enddo        enddo
143    
144  C     _END_MASTER( myThid )  C     For tracer (method=1,2) set to northernmost zonal-mean value
145    C     at 90N to avoid sharp zonal gradients near the Pole.
146    C     For U (method=11,12) set to zero at 90N to minimize velocity
147    C     gradient at North Pole
148    C     For V (method=11,12) set to northernmost zonal value at 90N,
149    C     as is already done above in order to allow cross-PoleArctic flow
150          do j=ny_in,ny_in+2
151           if (y_in(j).eq.ninety) then
152            if (method.eq.1 .or. method.eq.2) then
153             NorthValue = 0
154             do i=1,nx_in
155              NorthValue = NorthValue + arrayin(i,j)
156             enddo
157             NorthValue = NorthValue / nx_in
158             do i=-1,nx_in+2
159              arrayin(i,j) = NorthValue
160             enddo
161            elseif (method.eq.11 .or. method.eq.12) then
162             do i=-1,nx_in+2
163              arrayin(i,j) = 0
164             enddo
165            endif
166           endif
167          enddo
168    
169    C     _END_MASTER( myThid )
170          
171        do bj = mybylo(mythid), mybyhi(mythid)        do bj = mybylo(mythid), mybyhi(mythid)
172         do bi = mybxlo(mythid), mybxhi(mythid)         do bi = mybxlo(mythid), mybxhi(mythid)
173    
# Line 192  C compute interpolation indices Line 215  C compute interpolation indices
215           enddo           enddo
216          enddo          enddo
217    
218          if (method .eq. 1) then          if (method.eq.1 .or. method.eq.11 .or. method.eq.21) then
219    
220  C bilinear interpolation  C bilinear interpolation
221           sp = 2           sp = 2
# Line 213  C bilinear interpolation Line 236  C bilinear interpolation
236             enddo             enddo
237            enddo            enddo
238           enddo           enddo
239          elseif (method .eq. 2) then          elseif (method .eq. 2 .or. method.eq.12 .or. method.eq.22) then
240    
241  C bicubic interpolation  C bicubic interpolation
242           sp = 4           sp = 4
Legend:
Removed from v.1.13  
changed lines
  Added in v.1.18
  ViewVC Help
Powered by ViewVC 1.1.22