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

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

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


Revision 1.11 - (hide annotations) (download)
Mon Jun 5 19:23:31 2006 UTC (18 years ago) by dimitri
Branch: MAIN
Changes since 1.10: +2 -2 lines
exf_set_uv.F vector rotation applied when usingCurvilinearGrid instead of
more restrictive condition useCubedSphereExchange

1 dimitri 1.11 C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_set_uv.F,v 1.10 2006/03/02 02:53:23 heimbach Exp $
2 dimitri 1.1 C $Name: $
3    
4     #include "EXF_OPTIONS.h"
5    
6     subroutine exf_set_uv(
7     & uvecfile, uvecstartdate, uvecperiod,
8 heimbach 1.5 & uvecstartdate1, uvecstartdate2,
9 dimitri 1.1 & exf_inscal_uvec, uvec, uvec0, uvec1, uvecmask,
10     & uvec_lon0, uvec_lon_inc, uvec_lat0, uvec_lat_inc,
11 heimbach 1.10 & uvec_nlon, uvec_nlat,
12     & uvec_remove_intercept, uvec_remove_slope,
13 dimitri 1.1 & vvecfile, vvecstartdate, vvecperiod,
14 heimbach 1.5 & vvecstartdate1, vvecstartdate2,
15 dimitri 1.1 & exf_inscal_vvec, vvec, vvec0, vvec1, vvecmask,
16     & vvec_lon0, vvec_lon_inc, vvec_lat0, vvec_lat_inc,
17 dimitri 1.3 & vvec_nlon, vvec_nlat,
18 heimbach 1.10 & vvec_remove_intercept, vvec_remove_slope,
19 dimitri 1.1 & mycurrenttime, mycurrentiter, mythid )
20    
21     c ==================================================================
22     c SUBROUTINE exf_set_uv
23     c ==================================================================
24     c
25     c o Read-in, interpolate, and rotate wind or wind stress vectors
26     c from a spherical-polar input grid to an arbitrary output grid.
27     c
28     c menemenlis@jpl.nasa.gov, 8-Dec-2003
29     c
30     c ==================================================================
31     c SUBROUTINE exf_set_uv
32     c ==================================================================
33    
34     implicit none
35    
36     c == global variables ==
37    
38     #include "EEPARAMS.h"
39     #include "SIZE.h"
40     #include "PARAMS.h"
41     #include "DYNVARS.h"
42     #include "GRID.h"
43    
44     #include "exf_param.h"
45     #include "exf_fields.h"
46     #include "exf_constants.h"
47    
48     #ifdef ALLOW_AUTODIFF
49     # include "ctrl.h"
50     # include "ctrl_dummy.h"
51     #endif
52    
53     c == routine arguments ==
54    
55     c *vec_lon_0, :: longitude and latitude of SouthWest
56     c *vec_lat_0 corner of global input grid for *vec
57     c *vec_nlon, *vec_nlat :: input x-grid and y-grid size for *vec
58     c *vec_lon_inc :: scalar x-grid increment for *vec
59     c *vec_lat_inc :: vector y-grid increments for *vec
60    
61 heimbach 1.5 character*(128) uvecfile, uvecfile0, uvecfile1
62 dimitri 1.1 _RL uvecstartdate, uvecperiod
63 jmc 1.6 integer uvecstartdate1, uvecstartdate2
64 dimitri 1.1 _RL exf_inscal_uvec
65 heimbach 1.10 _RL uvec_remove_intercept, uvec_remove_slope
66 dimitri 1.1 _RL uvec (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
67     _RL uvec0 (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
68     _RL uvec1 (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
69     character*1 uvecmask
70     _RL uvec_lon0, uvec_lon_inc
71     _RL uvec_lat0, uvec_lat_inc(MAX_LAT_INC)
72     INTEGER uvec_nlon, uvec_nlat
73 heimbach 1.5 character*(128) vvecfile, vvecfile0, vvecfile1
74 dimitri 1.1 _RL vvecstartdate, vvecperiod
75 jmc 1.6 integer vvecstartdate1, vvecstartdate2
76 dimitri 1.1 _RL exf_inscal_vvec
77 heimbach 1.10 _RL vvec_remove_intercept, vvec_remove_slope
78 dimitri 1.1 _RL vvec (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
79     _RL vvec0 (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
80     _RL vvec1 (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
81     character*1 vvecmask
82     _RL vvec_lon0, vvec_lon_inc
83     _RL vvec_lat0, vvec_lat_inc(MAX_LAT_INC)
84     INTEGER vvec_nlon, vvec_nlat
85     _RL mycurrenttime
86     integer mycurrentiter
87     integer mythid
88    
89 heimbach 1.4 #ifdef USE_EXF_INTERPOLATION
90 dimitri 1.1 c == local variables ==
91    
92     logical first, changed
93     _RL fac, x1, x2, x3, x4, y1, y2, y3, y4, dx, dy
94     _RL tmp_u (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
95     _RL tmp_v (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
96     integer count0, count1
97     integer i, j, bi, bj
98 heimbach 1.5 integer il, interp_method
99 dimitri 1.1 parameter(interp_method=2)
100 heimbach 1.5 integer year0, year1
101    
102     c == external ==
103    
104     integer ilnblnk
105     external ilnblnk
106 dimitri 1.1
107     c == end of interface ==
108    
109 dimitri 1.11 IF ( usingCurvilinearGrid ) THEN
110 dimitri 1.1
111     if ( uvecfile .NE. ' ' .and. vvecfile .NE. ' ' ) then
112    
113     c get record numbers and interpolation factor
114     call exf_GetFFieldRec(
115     I uvecstartdate, uvecperiod
116 heimbach 1.5 I , uvecstartdate1, uvecstartdate2
117     I , useExfYearlyFields
118 dimitri 1.1 O , fac, first, changed
119 heimbach 1.5 O , count0, count1, year0, year1
120 dimitri 1.1 I , mycurrenttime, mycurrentiter, mythid
121     & )
122    
123     if ( first ) then
124 heimbach 1.5 if (useExfYearlyFields) then
125 dimitri 1.8 C Complete filenames with YR or _YEAR extension
126 heimbach 1.5 il = ilnblnk( uvecfile )
127 dimitri 1.8 if (twoDigitYear) then
128     if (year0.ge.2000) then
129     write(uvecfile0(1:128),'(a,i2.2)')
130     & uvecfile(1:il),year0-2000
131     else
132     write(uvecfile0(1:128),'(a,i2.2)')
133     & uvecfile(1:il),year0-1900
134     endif
135     else
136     write(uvecfile0(1:128),'(2a,i4.4)')
137     & uvecfile(1:il),'_',year0
138     endif
139 heimbach 1.5 il = ilnblnk( vvecfile )
140 dimitri 1.8 if (twoDigitYear) then
141     if (year0.ge.2000) then
142     write(vvecfile0(1:128),'(a,i2.2)')
143     & vvecfile(1:il),year0-2000
144     else
145     write(vvecfile0(1:128),'(a,i2.2)')
146     & vvecfile(1:il),year0-1900
147     endif
148     else
149     write(vvecfile0(1:128),'(2a,i4.4)')
150     & vvecfile(1:il),'_',year0
151     endif
152 heimbach 1.5 else
153     uvecfile0 = uvecfile
154     vvecfile0 = vvecfile
155     endif
156 dimitri 1.1 c scalar interpolation to (xC,yC) locations
157 heimbach 1.5 call exf_interp( uvecfile0, exf_iprec
158 dimitri 1.1 & , tmp_u, count0, xC, yC
159     & , uvec_lon0,uvec_lon_inc
160     & , uvec_lat0,uvec_lat_inc
161     & , uvec_nlon,uvec_nlat,interp_method,mythid
162     & )
163 heimbach 1.5 call exf_interp( vvecfile0, exf_iprec
164 dimitri 1.1 & , tmp_v, count0, xC, yC
165     & , vvec_lon0,vvec_lon_inc
166     & , vvec_lat0,vvec_lat_inc
167     & , vvec_nlon,vvec_nlat,interp_method,mythid
168     & )
169     c vector rotation
170     do bj = mybylo(mythid),mybyhi(mythid)
171     do bi = mybxlo(mythid),mybxhi(mythid)
172     do j = 1,sny
173     do i = 1,snx
174     x1=xG(i,j,bi,bj)
175     x2=xG(i+1,j,bi,bj)
176     x3=xG(i,j+1,bi,bj)
177     x4=xG(i+1,j+1,bi,bj)
178     if ((x2-x1).gt.180) x2=x2-360
179     if ((x1-x2).gt.180) x2=x2+360
180     if ((x3-x1).gt.180) x3=x3-360
181     if ((x1-x3).gt.180) x3=x3+360
182     if ((x4-x1).gt.180) x4=x4-360
183     if ((x1-x4).gt.180) x4=x4+360
184     y1=yG(i,j,bi,bj)
185     y2=yG(i+1,j,bi,bj)
186     y3=yG(i,j+1,bi,bj)
187     y4=yG(i+1,j+1,bi,bj)
188     dx=0.5*(x3+x4-x1-x2)
189     dx=dx*cos(deg2rad*yC(i,j,bi,bj))
190     dy=0.5*(y3+y4-y1-y2)
191     vvec1(i,j,bi,bj)=(tmp_u(i,j,bi,bj)*dx+
192     & tmp_v(i,j,bi,bj)*dy)/sqrt(dx*dx+dy*dy)
193     dx=0.5*(x2+x4-x1-x3)
194     dx=dx*cos(deg2rad*yC(i,j,bi,bj))
195     dy=0.5*(y2+y4-y1-y3)
196     uvec1(i,j,bi,bj)=(tmp_u(i,j,bi,bj)*dx+
197     & tmp_v(i,j,bi,bj)*dy)/sqrt(dx*dx+dy*dy)
198     enddo
199     enddo
200     enddo
201     enddo
202     c apply mask
203     if (exf_yftype .eq. 'RL') then
204     call exf_filter_rl( uvec1, uvecmask, mythid )
205     call exf_filter_rl( vvec1, vvecmask, mythid )
206     else
207     call exf_filter_rs( uvec1, uvecmask, mythid )
208     call exf_filter_rs( vvec1, vvecmask, mythid )
209     end if
210     endif
211    
212     if (( first ) .or. ( changed )) then
213     call exf_SwapFFields( uvec0, uvec1, mythid )
214     call exf_SwapFFields( vvec0, vvec1, mythid )
215 heimbach 1.5 if (useExfYearlyFields) then
216 dimitri 1.8 C Complete filenames with YR or _YEAR extension
217 heimbach 1.5 il = ilnblnk( uvecfile )
218 dimitri 1.8 if (twoDigitYear) then
219     if (year1.ge.2000) then
220     write(uvecfile1(1:128),'(a,i2.2)')
221     & uvecfile(1:il),year1-2000
222     else
223     write(uvecfile1(1:128),'(a,i2.2)')
224     & uvecfile(1:il),year1-1900
225     endif
226     else
227     write(uvecfile1(1:128),'(2a,i4.4)')
228     & uvecfile(1:il),'_',year1
229     endif
230 heimbach 1.5 il = ilnblnk( vvecfile )
231 dimitri 1.8 if (twoDigitYear) then
232     if (year1.ge.2000) then
233     write(vvecfile1(1:128),'(a,i2.2)')
234     & vvecfile(1:il),year1-2000
235     else
236     write(vvecfile1(1:128),'(a,i2.2)')
237     & vvecfile(1:il),year1-1900
238     endif
239     else
240     write(vvecfile1(1:128),'(2a,i4.4)')
241     & vvecfile(1:il),'_',year1
242     endif
243 heimbach 1.5 else
244     uvecfile1 = uvecfile
245     vvecfile1 = vvecfile
246     endif
247 dimitri 1.1 c scalar interpolation to (xC,yC) locations
248 heimbach 1.5 call exf_interp( uvecfile1, exf_iprec
249 dimitri 1.1 & , tmp_u, count1, xC, yC
250     & , uvec_lon0,uvec_lon_inc
251     & , uvec_lat0,uvec_lat_inc
252     & , uvec_nlon,uvec_nlat,interp_method,mythid
253     & )
254 heimbach 1.5 call exf_interp( vvecfile1, exf_iprec
255 dimitri 1.1 & , tmp_v, count1, xC, yC
256     & , vvec_lon0,vvec_lon_inc
257     & , vvec_lat0,vvec_lat_inc
258     & , vvec_nlon,vvec_nlat,interp_method,mythid
259     & )
260     c vector rotation
261     do bj = mybylo(mythid),mybyhi(mythid)
262     do bi = mybxlo(mythid),mybxhi(mythid)
263     do j = 1,sny
264     do i = 1,snx
265     x1=xG(i,j,bi,bj)
266     x2=xG(i+1,j,bi,bj)
267     x3=xG(i,j+1,bi,bj)
268     x4=xG(i+1,j+1,bi,bj)
269     if ((x2-x1).gt.180) x2=x2-360
270     if ((x1-x2).gt.180) x2=x2+360
271     if ((x3-x1).gt.180) x3=x3-360
272     if ((x1-x3).gt.180) x3=x3+360
273     if ((x4-x1).gt.180) x4=x4-360
274     if ((x1-x4).gt.180) x4=x4+360
275     y1=yG(i,j,bi,bj)
276     y2=yG(i+1,j,bi,bj)
277     y3=yG(i,j+1,bi,bj)
278     y4=yG(i+1,j+1,bi,bj)
279     dx=0.5*(x3+x4-x1-x2)
280     dx=dx*cos(deg2rad*yC(i,j,bi,bj))
281     dy=0.5*(y3+y4-y1-y2)
282     vvec1(i,j,bi,bj)=(tmp_u(i,j,bi,bj)*dx+
283     & tmp_v(i,j,bi,bj)*dy)/sqrt(dx*dx+dy*dy)
284     dx=0.5*(x2+x4-x1-x3)
285     dx=dx*cos(deg2rad*yC(i,j,bi,bj))
286     dy=0.5*(y2+y4-y1-y3)
287     uvec1(i,j,bi,bj)=(tmp_u(i,j,bi,bj)*dx+
288     & tmp_v(i,j,bi,bj)*dy)/sqrt(dx*dx+dy*dy)
289     enddo
290     enddo
291     enddo
292     enddo
293     c apply mask
294     if (exf_yftype .eq. 'RL') then
295     call exf_filter_rl( uvec1, uvecmask, mythid )
296     call exf_filter_rl( vvec1, vvecmask, mythid )
297     else
298     call exf_filter_rs( uvec1, uvecmask, mythid )
299     call exf_filter_rs( vvec1, vvecmask, mythid )
300     end if
301     endif
302    
303     c Interpolate linearly onto the current time.
304     do bj = mybylo(mythid),mybyhi(mythid)
305     do bi = mybxlo(mythid),mybxhi(mythid)
306     do j = 1,sny
307     do i = 1,snx
308     uvec(i,j,bi,bj) = exf_inscal_uvec * (
309     & fac * uvec0(i,j,bi,bj) +
310     & (exf_one - fac) * uvec1(i,j,bi,bj) )
311     vvec(i,j,bi,bj) = exf_inscal_vvec * (
312     & fac * vvec0(i,j,bi,bj) +
313     & (exf_one - fac) * vvec1(i,j,bi,bj) )
314     enddo
315     enddo
316     enddo
317     enddo
318    
319     endif
320    
321     ELSE
322     c IF ( .NOT. useCubedSphereExchange )
323    
324     call exf_set_gen(
325     & uvecfile, uvecstartdate, uvecperiod,
326 heimbach 1.5 I uvecstartdate1, uvecstartdate2,
327 dimitri 1.1 & exf_inscal_uvec,
328 heimbach 1.10 & uvec_remove_intercept, uvec_remove_slope,
329 dimitri 1.1 & uvec, uvec0, uvec1, uvecmask,
330     & uvec_lon0, uvec_lon_inc, uvec_lat0, uvec_lat_inc,
331 dimitri 1.7 & uvec_nlon, uvec_nlat, xC, yC, interp_method,
332 dimitri 1.1 & mycurrenttime, mycurrentiter, mythid )
333     call exf_set_gen(
334     & vvecfile, vvecstartdate, vvecperiod,
335 heimbach 1.5 I vvecstartdate1, vvecstartdate2,
336 dimitri 1.1 & exf_inscal_vvec,
337 heimbach 1.10 & vvec_remove_intercept, vvec_remove_slope,
338 dimitri 1.1 & vvec, vvec0, vvec1, vvecmask,
339     & vvec_lon0, vvec_lon_inc, vvec_lat0, vvec_lat_inc,
340 dimitri 1.7 & vvec_nlon, vvec_nlat, xC, yC, interp_method,
341 dimitri 1.1 & mycurrenttime, mycurrentiter, mythid )
342    
343     ENDIF
344    
345 heimbach 1.4 #endif /* USE_EXF_INTERPOLATION */
346    
347 dimitri 1.1 return
348     end

  ViewVC Help
Powered by ViewVC 1.1.22