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

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

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

revision 1.5 by heimbach, Mon Oct 11 16:41:01 2004 UTC revision 1.19 by dimitri, Tue Jan 29 11:25:53 2008 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
 #include "PACKAGES_CONFIG.h"  
4  #include "EXF_OPTIONS.h"  #include "EXF_OPTIONS.h"
5    
6        subroutine exf_set_uv(        subroutine exf_set_uv(
7       &     uvecfile, uvecstartdate, uvecperiod,       &     uvecfile, uvecstartdate, uvecperiod,
      &     uvecstartdate1, uvecstartdate2,  
8       &     exf_inscal_uvec, uvec, uvec0, uvec1, uvecmask,       &     exf_inscal_uvec, uvec, uvec0, uvec1, uvecmask,
9       &     uvec_lon0, uvec_lon_inc, uvec_lat0, uvec_lat_inc,       &     uvec_lon0, uvec_lon_inc, uvec_lat0, uvec_lat_inc,
10       &     uvec_nlon, uvec_nlat,       &     uvec_nlon, uvec_nlat,
11         &     uvec_remove_intercept, uvec_remove_slope,
12       &     vvecfile, vvecstartdate, vvecperiod,       &     vvecfile, vvecstartdate, vvecperiod,
      &     vvecstartdate1, vvecstartdate2,  
13       &     exf_inscal_vvec, vvec, vvec0, vvec1, vvecmask,       &     exf_inscal_vvec, vvec, vvec0, vvec1, vvecmask,
14       &     vvec_lon0, vvec_lon_inc, vvec_lat0, vvec_lat_inc,       &     vvec_lon0, vvec_lon_inc, vvec_lat0, vvec_lat_inc,
15       &     vvec_nlon, vvec_nlat,       &     vvec_nlon, vvec_nlat,
16         &     vvec_remove_intercept, vvec_remove_slope,
17       &     mycurrenttime, mycurrentiter, mythid )       &     mycurrenttime, mycurrentiter, mythid )
18    
19  c     ==================================================================  c     ==================================================================
# Line 40  c     == global variables == Line 39  c     == global variables ==
39  #include "DYNVARS.h"  #include "DYNVARS.h"
40  #include "GRID.h"  #include "GRID.h"
41    
42  #include "exf_param.h"  #include "EXF_PARAM.h"
43  #include "exf_fields.h"  #include "EXF_FIELDS.h"
44  #include "exf_constants.h"  #include "EXF_CONSTANTS.h"
45    
46  #ifdef ALLOW_AUTODIFF  #ifdef ALLOW_AUTODIFF
47  # include "ctrl.h"  # include "ctrl.h"
# Line 59  c     *vec_lat_inc         :: vector y-g Line 58  c     *vec_lat_inc         :: vector y-g
58    
59        character*(128) uvecfile, uvecfile0, uvecfile1        character*(128) uvecfile, uvecfile0, uvecfile1
60        _RL     uvecstartdate, uvecperiod        _RL     uvecstartdate, uvecperiod
       _RL     uvecstartdate1, uvecstartdate2  
61        _RL     exf_inscal_uvec        _RL     exf_inscal_uvec
62          _RL     uvec_remove_intercept, uvec_remove_slope
63        _RL     uvec  (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL     uvec  (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
64        _RL     uvec0 (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL     uvec0 (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
65        _RL     uvec1 (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL     uvec1 (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
# Line 70  c     *vec_lat_inc         :: vector y-g Line 69  c     *vec_lat_inc         :: vector y-g
69        INTEGER uvec_nlon, uvec_nlat        INTEGER uvec_nlon, uvec_nlat
70        character*(128) vvecfile, vvecfile0, vvecfile1        character*(128) vvecfile, vvecfile0, vvecfile1
71        _RL     vvecstartdate, vvecperiod        _RL     vvecstartdate, vvecperiod
       _RL     vvecstartdate1, vvecstartdate2  
72        _RL     exf_inscal_vvec        _RL     exf_inscal_vvec
73          _RL     vvec_remove_intercept, vvec_remove_slope
74        _RL     vvec  (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL     vvec  (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
75        _RL     vvec0 (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL     vvec0 (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
76        _RL     vvec1 (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)        _RL     vvec1 (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
# Line 93  c     == local variables == Line 92  c     == local variables ==
92        integer count0, count1        integer count0, count1
93        integer i, j, bi, bj        integer i, j, bi, bj
94        integer   il, interp_method        integer   il, interp_method
       parameter(interp_method=2)  
95        integer year0, year1        integer year0, year1
96    
97  c     == external ==  c     == external ==
# Line 103  c     == external == Line 101  c     == external ==
101    
102  c     == end of interface ==  c     == end of interface ==
103    
104        IF ( useCubedSphereExchange ) THEN        IF ( usingCurvilinearGrid ) THEN
105    
106           if ( uvecfile .NE. ' ' .and. vvecfile .NE. ' ' ) then          if ( uvecfile .NE. ' ' .and. vvecfile .NE. ' ' ) then
107    
108  c     get record numbers and interpolation factor  c     get record numbers and interpolation factor
109              call exf_GetFFieldRec(            call exf_GetFFieldRec(
110       I           uvecstartdate, uvecperiod       I        uvecstartdate, uvecperiod
111       I           , uvecstartdate1, uvecstartdate2       I        , useExfYearlyFields
112       I           , useExfYearlyFields       O        , fac, first, changed
113       O           , fac, first, changed       O        , count0, count1, year0, year1
114       O           , count0, count1, year0, year1       I        , mycurrenttime, mycurrentiter, mythid
115       I           , mycurrenttime, mycurrentiter, mythid       &        )
116       &           )  
117              if ( first ) then
118              if ( first ) then  
119                 if (useExfYearlyFields) then              call exf_GetYearlyFieldName(
120                    il = ilnblnk( uvecfile )       I         useExfYearlyFields, twoDigitYear, uvecperiod, year0,
121                    write(uvecfile0(1:128),'(2a,i4.4)')       I         uvecfile,
122       &                 uvecfile(1:il),'_',year0       O         uvecfile0,
123                    il = ilnblnk( vvecfile )       I         mycurrenttime, mycurrentiter, mythid )
124                    write(vvecfile0(1:128),'(2a,i4.4)')              call exf_GetYearlyFieldName(
125       &                 vvecfile(1:il),'_',year0       I         useExfYearlyFields, twoDigitYear, vvecperiod, year0,
126                 else       I         vvecfile,
127                    uvecfile0 = uvecfile       O         vvecfile0,
128                    vvecfile0 = vvecfile       I         mycurrenttime, mycurrentiter, mythid )
129                 endif  
130  c     scalar interpolation to (xC,yC) locations  c     scalar interpolation to (xC,yC) locations
131                 call exf_interp( uvecfile0, exf_iprec              interp_method=12
132       &              , tmp_u, count0, xC, yC              call exf_interp( uvecfile0, exf_iprec
133       &              , uvec_lon0,uvec_lon_inc       &          , tmp_u, count0, xC, yC
134       &              , uvec_lat0,uvec_lat_inc       &          , uvec_lon0,uvec_lon_inc
135       &              , uvec_nlon,uvec_nlat,interp_method,mythid       &          , uvec_lat0,uvec_lat_inc
136       &              )       &          , uvec_nlon,uvec_nlat,interp_method,mythid
137                 call exf_interp( vvecfile0, exf_iprec       &          )
138       &              , tmp_v, count0, xC, yC              interp_method=22
139       &              , vvec_lon0,vvec_lon_inc              call exf_interp( vvecfile0, exf_iprec
140       &              , vvec_lat0,vvec_lat_inc       &          , tmp_v, count0, xC, yC
141       &              , vvec_nlon,vvec_nlat,interp_method,mythid       &          , vvec_lon0,vvec_lon_inc
142       &              )       &          , vvec_lat0,vvec_lat_inc
143         &          , vvec_nlon,vvec_nlat,interp_method,mythid
144         &          )
145  c     vector rotation  c     vector rotation
146                 do bj = mybylo(mythid),mybyhi(mythid)              do bj = mybylo(mythid),mybyhi(mythid)
147                    do bi = mybxlo(mythid),mybxhi(mythid)                do bi = mybxlo(mythid),mybxhi(mythid)
148                       do j = 1,sny                  do j = 1,sny
149                          do i = 1,snx                    do i = 1,snx
150                             x1=xG(i,j,bi,bj)                      x1=xG(i,j,bi,bj)
151                             x2=xG(i+1,j,bi,bj)                      x2=xG(i+1,j,bi,bj)
152                             x3=xG(i,j+1,bi,bj)                      x3=xG(i,j+1,bi,bj)
153                             x4=xG(i+1,j+1,bi,bj)                      x4=xG(i+1,j+1,bi,bj)
154                             if ((x2-x1).gt.180) x2=x2-360                      if ((x2-x1).gt.180) x2=x2-360
155                             if ((x1-x2).gt.180) x2=x2+360                      if ((x1-x2).gt.180) x2=x2+360
156                             if ((x3-x1).gt.180) x3=x3-360                      if ((x3-x1).gt.180) x3=x3-360
157                             if ((x1-x3).gt.180) x3=x3+360                      if ((x1-x3).gt.180) x3=x3+360
158                             if ((x4-x1).gt.180) x4=x4-360                      if ((x4-x1).gt.180) x4=x4-360
159                             if ((x1-x4).gt.180) x4=x4+360                      if ((x1-x4).gt.180) x4=x4+360
160                             y1=yG(i,j,bi,bj)                      y1=yG(i,j,bi,bj)
161                             y2=yG(i+1,j,bi,bj)                      y2=yG(i+1,j,bi,bj)
162                             y3=yG(i,j+1,bi,bj)                      y3=yG(i,j+1,bi,bj)
163                             y4=yG(i+1,j+1,bi,bj)                      y4=yG(i+1,j+1,bi,bj)
164                             dx=0.5*(x3+x4-x1-x2)                      dx=0.5*(x3+x4-x1-x2)
165                             dx=dx*cos(deg2rad*yC(i,j,bi,bj))                      dx=dx*
166                             dy=0.5*(y3+y4-y1-y2)       &                  cos(deg2rad*yC(i,j,bi,bj))
167                             vvec1(i,j,bi,bj)=(tmp_u(i,j,bi,bj)*dx+                      dy=0.5*(y3+y4-y1-y2)
168       &                          tmp_v(i,j,bi,bj)*dy)/sqrt(dx*dx+dy*dy)                      vvec1(i,j,bi,bj)=
169                             dx=0.5*(x2+x4-x1-x3)       &                  (tmp_u(i,j,bi,bj)*dx+
170                             dx=dx*cos(deg2rad*yC(i,j,bi,bj))       &                  tmp_v(i,j,bi,bj)*dy)/
171                             dy=0.5*(y2+y4-y1-y3)       &                  sqrt(dx*dx+dy*dy)
172                             uvec1(i,j,bi,bj)=(tmp_u(i,j,bi,bj)*dx+                      dx=0.5*(x2+x4-x1-x3)
173       &                          tmp_v(i,j,bi,bj)*dy)/sqrt(dx*dx+dy*dy)                      dx=dx*
174                          enddo       &                  cos(deg2rad*yC(i,j,bi,bj))
175                       enddo                      dy=0.5*(y2+y4-y1-y3)
176                        uvec1(i,j,bi,bj)=
177         &                  (tmp_u(i,j,bi,bj)*dx+
178         &                  tmp_v(i,j,bi,bj)*dy)/
179         &                  sqrt(dx*dx+dy*dy)
180                    enddo                    enddo
181                 enddo                  enddo
182                  enddo
183                enddo
184  c     apply mask  c     apply mask
185                 if (exf_yftype .eq. 'RL') then              if (exf_yftype .eq. 'RL') then
186                    call exf_filter_rl( uvec1, uvecmask, mythid )                call exf_filter_rl( uvec1, uvecmask, mythid )
187                    call exf_filter_rl( vvec1, vvecmask, mythid )                call exf_filter_rl( vvec1, vvecmask, mythid )
188                 else              else
189                    call exf_filter_rs( uvec1, uvecmask, mythid )                call exf_filter_rs( uvec1, uvecmask, mythid )
190                    call exf_filter_rs( vvec1, vvecmask, mythid )                call exf_filter_rs( vvec1, vvecmask, mythid )
191                 end if              end if
192              endif            endif
193    
194              if (( first ) .or. ( changed )) then            if (( first ) .or. ( changed )) then
195                 call exf_SwapFFields( uvec0, uvec1, mythid )              call exf_SwapFFields( uvec0, uvec1, mythid )
196                 call exf_SwapFFields( vvec0, vvec1, mythid )              call exf_SwapFFields( vvec0, vvec1, mythid )
197                 if (useExfYearlyFields) then  
198                    il = ilnblnk( uvecfile )              call exf_GetYearlyFieldName(
199                    write(uvecfile1(1:128),'(2a,i4.4)')       I         useExfYearlyFields, twoDigitYear, uvecperiod, year1,
200       &                 uvecfile(1:il),'_',year0       I         uvecfile,
201                    il = ilnblnk( vvecfile )       O         uvecfile1,
202                    write(vvecfile1(1:128),'(2a,i4.4)')       I         mycurrenttime, mycurrentiter, mythid )
203       &                 vvecfile(1:il),'_',year0              call exf_GetYearlyFieldName(
204                 else       I         useExfYearlyFields, twoDigitYear, vvecperiod, year1,
205                    uvecfile1 = uvecfile       I         vvecfile,
206                    vvecfile1 = vvecfile       O         vvecfile1,
207                 endif       I         mycurrenttime, mycurrentiter, mythid )
208    
209  c     scalar interpolation to (xC,yC) locations  c     scalar interpolation to (xC,yC) locations
210                 call exf_interp( uvecfile1, exf_iprec              interp_method=12
211       &              , tmp_u, count1, xC, yC              call exf_interp( uvecfile1, exf_iprec
212       &              , uvec_lon0,uvec_lon_inc       &          , tmp_u, count1, xC, yC
213       &              , uvec_lat0,uvec_lat_inc       &          , uvec_lon0,uvec_lon_inc
214       &              , uvec_nlon,uvec_nlat,interp_method,mythid       &          , uvec_lat0,uvec_lat_inc
215       &              )       &          , uvec_nlon,uvec_nlat,interp_method,mythid
216                 call exf_interp( vvecfile1, exf_iprec       &          )
217       &              , tmp_v, count1, xC, yC              interp_method=22
218       &              , vvec_lon0,vvec_lon_inc              call exf_interp( vvecfile1, exf_iprec
219       &              , vvec_lat0,vvec_lat_inc       &          , tmp_v, count1, xC, yC
220       &              , vvec_nlon,vvec_nlat,interp_method,mythid       &          , vvec_lon0,vvec_lon_inc
221       &              )       &          , vvec_lat0,vvec_lat_inc
222         &          , vvec_nlon,vvec_nlat,interp_method,mythid
223         &          )
224  c     vector rotation  c     vector rotation
225                 do bj = mybylo(mythid),mybyhi(mythid)              do bj = mybylo(mythid),mybyhi(mythid)
226                    do bi = mybxlo(mythid),mybxhi(mythid)                do bi = mybxlo(mythid),mybxhi(mythid)
227                       do j = 1,sny                  do j = 1,sny
228                          do i = 1,snx                    do i = 1,snx
229                             x1=xG(i,j,bi,bj)                      x1=xG(i,j,bi,bj)
230                             x2=xG(i+1,j,bi,bj)                      x2=xG(i+1,j,bi,bj)
231                             x3=xG(i,j+1,bi,bj)                      x3=xG(i,j+1,bi,bj)
232                             x4=xG(i+1,j+1,bi,bj)                      x4=xG(i+1,j+1,bi,bj)
233                             if ((x2-x1).gt.180) x2=x2-360                      if ((x2-x1).gt.180) x2=x2-360
234                             if ((x1-x2).gt.180) x2=x2+360                      if ((x1-x2).gt.180) x2=x2+360
235                             if ((x3-x1).gt.180) x3=x3-360                      if ((x3-x1).gt.180) x3=x3-360
236                             if ((x1-x3).gt.180) x3=x3+360                      if ((x1-x3).gt.180) x3=x3+360
237                             if ((x4-x1).gt.180) x4=x4-360                      if ((x4-x1).gt.180) x4=x4-360
238                             if ((x1-x4).gt.180) x4=x4+360                      if ((x1-x4).gt.180) x4=x4+360
239                             y1=yG(i,j,bi,bj)                      y1=yG(i,j,bi,bj)
240                             y2=yG(i+1,j,bi,bj)                      y2=yG(i+1,j,bi,bj)
241                             y3=yG(i,j+1,bi,bj)                      y3=yG(i,j+1,bi,bj)
242                             y4=yG(i+1,j+1,bi,bj)                      y4=yG(i+1,j+1,bi,bj)
243                             dx=0.5*(x3+x4-x1-x2)                      dx=0.5*(x3+x4-x1-x2)
244                             dx=dx*cos(deg2rad*yC(i,j,bi,bj))                      dx=dx*
245                             dy=0.5*(y3+y4-y1-y2)       &                  cos(deg2rad*yC(i,j,bi,bj))
246                             vvec1(i,j,bi,bj)=(tmp_u(i,j,bi,bj)*dx+                      dy=0.5*(y3+y4-y1-y2)
247       &                          tmp_v(i,j,bi,bj)*dy)/sqrt(dx*dx+dy*dy)                      vvec1(i,j,bi,bj)=
248                             dx=0.5*(x2+x4-x1-x3)       &                  (tmp_u(i,j,bi,bj)*dx+
249                             dx=dx*cos(deg2rad*yC(i,j,bi,bj))       &                  tmp_v(i,j,bi,bj)*dy)/
250                             dy=0.5*(y2+y4-y1-y3)       &                  sqrt(dx*dx+dy*dy)
251                             uvec1(i,j,bi,bj)=(tmp_u(i,j,bi,bj)*dx+                      dx=0.5*(x2+x4-x1-x3)
252       &                          tmp_v(i,j,bi,bj)*dy)/sqrt(dx*dx+dy*dy)                      dx=dx*
253                          enddo       &                  cos(deg2rad*yC(i,j,bi,bj))
254                       enddo                      dy=0.5*(y2+y4-y1-y3)
255                        uvec1(i,j,bi,bj)=
256         &                  (tmp_u(i,j,bi,bj)*dx+
257         &                  tmp_v(i,j,bi,bj)*dy)/
258         &                  sqrt(dx*dx+dy*dy)
259                    enddo                    enddo
260                 enddo                  enddo
261                  enddo
262                enddo
263  c     apply mask  c     apply mask
264                 if (exf_yftype .eq. 'RL') then              if (exf_yftype .eq. 'RL') then
265                    call exf_filter_rl( uvec1, uvecmask, mythid )                call exf_filter_rl( uvec1, uvecmask, mythid )
266                    call exf_filter_rl( vvec1, vvecmask, mythid )                call exf_filter_rl( vvec1, vvecmask, mythid )
267                 else              else
268                    call exf_filter_rs( uvec1, uvecmask, mythid )                call exf_filter_rs( uvec1, uvecmask, mythid )
269                    call exf_filter_rs( vvec1, vvecmask, mythid )                call exf_filter_rs( vvec1, vvecmask, mythid )
270                 end if              end if
271              endif            endif
272    
273  c     Interpolate linearly onto the current time.  c     Interpolate linearly onto the current time.
274              do bj = mybylo(mythid),mybyhi(mythid)            do bj = mybylo(mythid),mybyhi(mythid)
275                 do bi = mybxlo(mythid),mybxhi(mythid)              do bi = mybxlo(mythid),mybxhi(mythid)
276                    do j = 1,sny                do j = 1,sny
277                       do i = 1,snx                  do i = 1,snx
278                          uvec(i,j,bi,bj) = exf_inscal_uvec * (                    uvec(i,j,bi,bj) = exf_inscal_uvec * (
279       &                       fac * uvec0(i,j,bi,bj) +       &                fac * uvec0(i,j,bi,bj) +
280       &                       (exf_one - fac) * uvec1(i,j,bi,bj) )       &                (exf_one - fac) * uvec1(i,j,bi,bj) )
281                          vvec(i,j,bi,bj) = exf_inscal_vvec * (                    vvec(i,j,bi,bj) = exf_inscal_vvec * (
282       &                       fac * vvec0(i,j,bi,bj) +       &                fac * vvec0(i,j,bi,bj) +
283       &                       (exf_one - fac) * vvec1(i,j,bi,bj) )       &                (exf_one - fac) * vvec1(i,j,bi,bj) )
284                       enddo                  enddo
285                    enddo                enddo
                enddo  
286              enddo              enddo
287              enddo
288           endif            
289            endif
290    
291        ELSE        ELSE
292  c     IF ( .NOT. useCubedSphereExchange )  c     IF ( .NOT. usingCurvilinearGrid )
   
             call exf_set_gen(  
      &           uvecfile, uvecstartdate, uvecperiod,  
      I           uvecstartdate1, uvecstartdate2,  
      &           exf_inscal_uvec,  
      &           uvec, uvec0, uvec1, uvecmask,  
      &           uvec_lon0, uvec_lon_inc, uvec_lat0, uvec_lat_inc,  
      &           uvec_nlon, uvec_nlat, xC, yC,  
      &           mycurrenttime, mycurrentiter, mythid )  
             call exf_set_gen(  
      &           vvecfile, vvecstartdate, vvecperiod,  
      I           vvecstartdate1, vvecstartdate2,  
      &           exf_inscal_vvec,  
      &           vvec, vvec0, vvec1, vvecmask,  
      &           vvec_lon0, vvec_lon_inc, vvec_lat0, vvec_lat_inc,  
      &           vvec_nlon, vvec_nlat, xC, yC,  
      &           mycurrenttime, mycurrentiter, mythid )  
293    
294           interp_method=12
295            call exf_set_gen(
296         &      uvecfile, uvecstartdate, uvecperiod,
297         &      exf_inscal_uvec,
298         &      uvec_remove_intercept, uvec_remove_slope,
299         &      uvec, uvec0, uvec1, uvecmask,
300         &      uvec_lon0, uvec_lon_inc, uvec_lat0, uvec_lat_inc,
301         &      uvec_nlon, uvec_nlat, xC, yC, interp_method,
302         &      mycurrenttime, mycurrentiter, mythid )
303            interp_method=22
304            call exf_set_gen(
305         &      vvecfile, vvecstartdate, vvecperiod,
306         &      exf_inscal_vvec,
307         &      vvec_remove_intercept, vvec_remove_slope,
308         &      vvec, vvec0, vvec1, vvecmask,
309         &      vvec_lon0, vvec_lon_inc, vvec_lat0, vvec_lat_inc,
310         &      vvec_nlon, vvec_nlat, xC, yC, interp_method,
311         &      mycurrenttime, mycurrentiter, mythid )
312            
313        ENDIF        ENDIF
314    
315  #endif /* USE_EXF_INTERPOLATION */  #endif /* USE_EXF_INTERPOLATION */
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.19
  ViewVC Help
Powered by ViewVC 1.1.22