1 |
C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_set_uv.F,v 1.36 2017/02/21 01:04:26 jmc Exp $ |
2 |
C $Name: $ |
3 |
|
4 |
#include "EXF_OPTIONS.h" |
5 |
|
6 |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
7 |
|
8 |
C !ROUTINE: EXF_INIT_FLD |
9 |
C !INTERFACE: |
10 |
SUBROUTINE EXF_SET_UV( |
11 |
I uVecName, uVecFile, uVecMask, |
12 |
I uVecStartTime, uVecPeriod, uVecRepeatCycle, |
13 |
I uVec_inScale, uVec_remove_intercept, uVec_remove_slope, |
14 |
U uVec, uVec0, uVec1, |
15 |
I vVecName, vVecFile, vVecMask, |
16 |
I vVecStartTime, vVecPeriod, vVecRepeatCycle, |
17 |
I vVec_inScale, vVec_remove_intercept, vVec_remove_slope, |
18 |
U vVec, vVec0, vVec1, |
19 |
#ifdef USE_EXF_INTERPOLATION |
20 |
I uVec_lon0, uVec_lon_inc, uVec_lat0, uVec_lat_inc, |
21 |
I uVec_nlon, uVec_nlat, u_interp_method, |
22 |
I vVec_lon0, vVec_lon_inc, vVec_lat0, vVec_lat_inc, |
23 |
I vVec_nlon, vVec_nlat, v_interp_method, uvInterp, |
24 |
#endif /* USE_EXF_INTERPOLATION */ |
25 |
I myTime, myIter, myThid ) |
26 |
|
27 |
C !DESCRIPTION: \bv |
28 |
C *=================================================================* |
29 |
C | SUBROUTINE EXF_SET_UV |
30 |
C | o Read-in, interpolate, and rotate wind or wind stress vectors |
31 |
C | from a spherical-polar input grid to an arbitrary output grid. |
32 |
C *=================================================================* |
33 |
C | menemenlis@jpl.nasa.gov, 8-Dec-2003 |
34 |
C *=================================================================* |
35 |
C \ev |
36 |
|
37 |
C !USES: |
38 |
IMPLICIT NONE |
39 |
C == global variables == |
40 |
#include "EEPARAMS.h" |
41 |
#include "SIZE.h" |
42 |
#include "PARAMS.h" |
43 |
#include "GRID.h" |
44 |
#include "EXF_INTERP_SIZE.h" |
45 |
#include "EXF_PARAM.h" |
46 |
#include "EXF_CONSTANTS.h" |
47 |
#include "EXF_FIELDS.h" |
48 |
|
49 |
C !INPUT/OUTPUT PARAMETERS: |
50 |
C *VecName :: vector compon. short name (to print mesg) |
51 |
C *VecFile :: file-name for this vector compon. field |
52 |
C *VecStartTime :: corresponding starting time (in sec) for this vec |
53 |
C *VecPeriod :: time period (in sec) between 2 reccords |
54 |
C *VecRepeatCycle :: time duration of a repeating cycle |
55 |
C *Vec_inScale :: input field scaling factor |
56 |
C *VecRemove_intercept :: |
57 |
C *VecRemove_slope :: |
58 |
C *Vec :: field array containing current time values |
59 |
C *Vec0 :: field array holding previous reccord |
60 |
C *Vec1 :: field array holding next reccord |
61 |
#ifdef USE_EXF_INTERPOLATION |
62 |
C *vec_lon0, *vec_lat0 :: longitude and latitude of SouthWest |
63 |
C :: corner of global input grid for *vec |
64 |
C *vec_nlon, *vec_nlat :: input x-grid and y-grid size for *vec |
65 |
C *_interp_method :: select interpolation method for *vec |
66 |
C *vec_lon_inc :: scalar x-grid increment for *vec |
67 |
C *vec_lat_inc :: vector y-grid increments for *vec |
68 |
#endif /* USE_EXF_INTERPOLATION */ |
69 |
C myTime :: Current time (in sec) in simulation |
70 |
C myIter :: Current iteration number |
71 |
C myThid :: My Thread Id number |
72 |
CHARACTER*(*) uVecName |
73 |
CHARACTER*(128) uVecFile |
74 |
CHARACTER*1 uVecMask |
75 |
_RL uVecStartTime, uVecPeriod, uVecRepeatCycle |
76 |
_RL uVec_inScale |
77 |
_RL uVec_remove_intercept, uVec_remove_slope |
78 |
_RL uVec (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
79 |
_RL uVec0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
80 |
_RL uVec1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
81 |
CHARACTER*(*) vVecName |
82 |
CHARACTER*(128) vVecFile |
83 |
CHARACTER*1 vVecMask |
84 |
_RL vVecStartTime, vVecPeriod, vVecRepeatCycle |
85 |
_RL vVec_inScale |
86 |
_RL vVec_remove_intercept, vVec_remove_slope |
87 |
_RL vVec (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
88 |
_RL vVec0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
89 |
_RL vVec1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
90 |
#ifdef USE_EXF_INTERPOLATION |
91 |
_RL uVec_lon0, uVec_lon_inc |
92 |
_RL uVec_lat0, uVec_lat_inc(MAX_LAT_INC) |
93 |
INTEGER uVec_nlon, uVec_nlat, u_interp_method |
94 |
_RL vVec_lon0, vVec_lon_inc |
95 |
_RL vVec_lat0, vVec_lat_inc(MAX_LAT_INC) |
96 |
INTEGER vVec_nlon, vVec_nlat, v_interp_method |
97 |
LOGICAL uvInterp |
98 |
#endif /* USE_EXF_INTERPOLATION */ |
99 |
_RL myTime |
100 |
INTEGER myIter |
101 |
INTEGER myThid |
102 |
|
103 |
C !FUNCTIONS: |
104 |
INTEGER ILNBLNK |
105 |
EXTERNAL ILNBLNK |
106 |
|
107 |
C !LOCAL VARIABLES: |
108 |
INTEGER i, j, bi, bj |
109 |
_RL tmp_u (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
110 |
_RL tmp_v (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) |
111 |
#ifdef USE_EXF_INTERPOLATION |
112 |
C msgBuf :: Informational/error message buffer |
113 |
CHARACTER*(MAX_LEN_MBUF) msgBuf |
114 |
CHARACTER*(128) uVecFile0, uVecFile1 |
115 |
CHARACTER*(128) vVecFile0, vVecFile1 |
116 |
CHARACTER*(MAX_LEN_FNAM) out_uVecFile, out_vVecFile |
117 |
LOGICAL first, changed |
118 |
_RL fac |
119 |
#ifdef EXF_USE_OLD_VEC_ROTATION |
120 |
_RL x1, x2, x3, x4, y1, y2, y3, y4, dx, dy |
121 |
#endif |
122 |
INTEGER count0, count1 |
123 |
INTEGER year0, year1 |
124 |
# ifndef EXF_INTERP_USE_DYNALLOC |
125 |
_RL bufArr1( exf_interp_bufferSize ) |
126 |
_RL bufArr2( exf_interp_bufferSize ) |
127 |
# endif |
128 |
#endif /* USE_EXF_INTERPOLATION */ |
129 |
CEOP |
130 |
|
131 |
#ifdef USE_EXF_INTERPOLATION |
132 |
IF ( u_interp_method.GE.1 .AND. v_interp_method.GE.1 .AND. |
133 |
& uVecFile.NE.' ' .AND. vVecFile.NE.' ' .AND. |
134 |
& (usingCurvilinearGrid .OR. rotateGrid .OR. uvInterp) ) THEN |
135 |
|
136 |
IF ( exf_debugLev.GE.debLevD ) THEN |
137 |
_BEGIN_MASTER( myThid ) |
138 |
i = ILNBLNK(uVecFile) |
139 |
j = ILNBLNK(vVecFile) |
140 |
WRITE(msgBuf,'(6A)') 'EXF_SET_UV: ', |
141 |
& 'processing fields "', uVecName, '" & "', vVecName, '"' |
142 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
143 |
& SQUEEZE_RIGHT, myThid ) |
144 |
WRITE(msgBuf,'(6A)') 'EXF_SET_UV: ', |
145 |
& ' files: ', uVecFile(1:i), ' & ', vVecFile(1:j) |
146 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
147 |
& SQUEEZE_RIGHT, myThid ) |
148 |
_END_MASTER( myThid ) |
149 |
ENDIF |
150 |
IF ( useCAL .AND. uVecPeriod .EQ. -12. ) THEN |
151 |
#ifdef ALLOW_CAL |
152 |
C- genperiod=-12 means input file contains 12 monthly means |
153 |
C records, corresponding to Jan. (rec=1) through Dec. (rec=12) |
154 |
CALL cal_GetMonthsRec( |
155 |
O fac, first, changed, |
156 |
O count0, count1, |
157 |
I myTime, myIter, myThid ) |
158 |
#endif /* ALLOW_CAL */ |
159 |
ELSEIF ( uVecPeriod .LT. 0. ) THEN |
160 |
j = ILNBLNK(uVecFile) |
161 |
WRITE(msgBuf,'(4A,1PE16.8,2A)') 'EXF_SET_UV: ', |
162 |
& '"', uVecName, '", Invalid uVecPeriod=', uVecPeriod, |
163 |
& ' for file: ', uVecFile(1:j) |
164 |
CALL PRINT_ERROR( msgBuf, myThid ) |
165 |
STOP 'ABNORMAL END: S/R EXF_SET_UV' |
166 |
ELSE |
167 |
C- get record numbers and interpolation factor |
168 |
CALL EXF_GetFFieldRec( |
169 |
I uVecStartTime, uVecPeriod, uVecRepeatCycle, |
170 |
I uVecName, useExfYearlyFields, |
171 |
O fac, first, changed, |
172 |
O count0, count1, year0, year1, |
173 |
I myTime, myIter, myThid ) |
174 |
ENDIF |
175 |
IF ( exf_debugLev.GE.debLevD ) THEN |
176 |
_BEGIN_MASTER( myThid ) |
177 |
WRITE(msgBuf,'(2A,I10,2I7)') 'EXF_SET_UV: ', |
178 |
& ' myIter, count0, count1:', myIter, count0, count1 |
179 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
180 |
& SQUEEZE_RIGHT, myThid ) |
181 |
WRITE(msgBuf,'(2A,2(L2,2X),E16.9)') 'EXF_SET_UV: ', |
182 |
& ' first, changed, fac: ', first, changed, fac |
183 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
184 |
& SQUEEZE_RIGHT, myThid ) |
185 |
_END_MASTER( myThid ) |
186 |
ENDIF |
187 |
|
188 |
IF ( first ) THEN |
189 |
C-- Load and interpolate a new reccord (= 1rst one of this run) |
190 |
|
191 |
CALL exf_GetYearlyFieldName( |
192 |
I useExfYearlyFields, twoDigitYear, uVecPeriod, year0, |
193 |
I uVecFile, |
194 |
O uVecFile0, |
195 |
I myTime, myIter, myThid ) |
196 |
CALL exf_GetYearlyFieldName( |
197 |
I useExfYearlyFields, twoDigitYear, vVecPeriod, year0, |
198 |
I vVecFile, |
199 |
O vVecFile0, |
200 |
I myTime, myIter, myThid ) |
201 |
IF ( exf_debugLev.GE.debLevC ) THEN |
202 |
_BEGIN_MASTER(myThid) |
203 |
WRITE(msgBuf,'(6A,I10)') 'EXF_SET_UV: ', |
204 |
& 'fields "', uVecName, '" & "', vVecName, '", it=', myIter |
205 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
206 |
& SQUEEZE_RIGHT, myThid ) |
207 |
j = ILNBLNK(uVecFile0) |
208 |
WRITE(msgBuf,'(2A,I6,3A)') 'EXF_SET_UV: ', |
209 |
& 'loading rec=', count0, ' from file: "', uVecFile0(1:j), '"' |
210 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
211 |
& SQUEEZE_RIGHT, myThid ) |
212 |
j = ILNBLNK(vVecFile0) |
213 |
WRITE(msgBuf,'(2A,I6,3A)') 'EXF_SET_UV: ', |
214 |
& 'loading rec=', count0, ' from file: "', vVecFile0(1:j), '"' |
215 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
216 |
& SQUEEZE_RIGHT, myThid ) |
217 |
_END_MASTER(myThid) |
218 |
ENDIF |
219 |
|
220 |
IF ( uvInterp ) THEN |
221 |
C- vector interpolation to (xC,yC) locations |
222 |
CALL EXF_INTERP_UV( |
223 |
I uVecFile0, vVecFile0, exf_iprec, count0, |
224 |
I uVec_nlon, uVec_nlat, |
225 |
I uVec_lon0, uVec_lon_inc, uVec_lat0, uVec_lat_inc, |
226 |
#ifdef EXF_INTERP_USE_DYNALLOC |
227 |
O tmp_u, tmp_v, |
228 |
#else |
229 |
O tmp_u, tmp_v, bufArr1, bufArr2, |
230 |
#endif |
231 |
I xC, yC, |
232 |
I u_interp_method, v_interp_method, myIter, myThid ) |
233 |
ELSE |
234 |
C- scalar interpolation to (xC,yC) locations |
235 |
CALL EXF_INTERP( |
236 |
I uVecFile0, exf_iprec, |
237 |
#ifdef EXF_INTERP_USE_DYNALLOC |
238 |
O tmp_u, |
239 |
#else |
240 |
O tmp_u, bufArr1, |
241 |
#endif |
242 |
I count0, xC, yC, |
243 |
I uVec_lon0, uVec_lon_inc, uVec_lat0, uVec_lat_inc, |
244 |
I uVec_nlon, uVec_nlat, u_interp_method, |
245 |
I myIter, myThid ) |
246 |
CALL EXF_INTERP( |
247 |
I vVecFile0, exf_iprec, |
248 |
#ifdef EXF_INTERP_USE_DYNALLOC |
249 |
O tmp_v, |
250 |
#else |
251 |
O tmp_v, bufArr2, |
252 |
#endif |
253 |
I count0, xC, yC, |
254 |
I vVec_lon0, vVec_lon_inc, vVec_lat0, vVec_lat_inc, |
255 |
I vVec_nlon, vVec_nlat, v_interp_method, |
256 |
I myIter, myThid ) |
257 |
ENDIF |
258 |
|
259 |
C- apply mask: Note: done after applying scaling factor and rotation |
260 |
c CALL EXF_FILTER_RL( tmp_u, uVecMask, myThid ) |
261 |
c CALL EXF_FILTER_RL( tmp_v, vVecMask, myThid ) |
262 |
|
263 |
IF ( exf_output_interp ) THEN |
264 |
j = ILNBLNK(uVecFile0) |
265 |
WRITE(out_uVecFile,'(2A)') uVecFile0(1:j), '_out' |
266 |
IF ( count0.NE.1 ) |
267 |
& CALL WRITE_REC_XY_RL(out_uVecFile,tmp_u,1,myIter,myThid) |
268 |
CALL WRITE_REC_XY_RL(out_uVecFile,tmp_u,count0,myIter,myThid) |
269 |
j = ILNBLNK(vVecFile0) |
270 |
WRITE(out_vVecFile,'(2A)') vVecFile0(1:j), '_out' |
271 |
IF ( count0.NE.1 ) |
272 |
& CALL WRITE_REC_XY_RL(out_vVecFile,tmp_v,1,myIter,myThid) |
273 |
CALL WRITE_REC_XY_RL(out_vVecFile,tmp_v,count0,myIter,myThid) |
274 |
ENDIF |
275 |
|
276 |
C- scaling factor and vector rotation |
277 |
IF ( usingCurvilinearGrid .OR. rotateGrid ) THEN |
278 |
DO bj = myByLo(myThid),myByHi(myThid) |
279 |
DO bi = myBxLo(myThid),myBxHi(myThid) |
280 |
DO j = 1,sNy |
281 |
DO i = 1,sNx |
282 |
tmp_u(i,j,bi,bj) = uVec_inScale*tmp_u(i,j,bi,bj) |
283 |
tmp_v(i,j,bi,bj) = vVec_inScale*tmp_v(i,j,bi,bj) |
284 |
ENDDO |
285 |
ENDDO |
286 |
DO j = 1,sNy |
287 |
DO i = 1,sNx |
288 |
#ifdef EXF_USE_OLD_VEC_ROTATION |
289 |
x1=xG(i,j,bi,bj) |
290 |
x2=xG(i+1,j,bi,bj) |
291 |
x3=xG(i,j+1,bi,bj) |
292 |
x4=xG(i+1,j+1,bi,bj) |
293 |
IF ((x2-x1).GT.180) x2=x2-360 |
294 |
IF ((x1-x2).GT.180) x2=x2+360 |
295 |
IF ((x3-x1).GT.180) x3=x3-360 |
296 |
IF ((x1-x3).GT.180) x3=x3+360 |
297 |
IF ((x4-x1).GT.180) x4=x4-360 |
298 |
IF ((x1-x4).GT.180) x4=x4+360 |
299 |
y1=yG(i,j,bi,bj) |
300 |
y2=yG(i+1,j,bi,bj) |
301 |
y3=yG(i,j+1,bi,bj) |
302 |
y4=yG(i+1,j+1,bi,bj) |
303 |
dx=0.5*(x3+x4-x1-x2) |
304 |
dx=dx* |
305 |
& cos(deg2rad*yC(i,j,bi,bj)) |
306 |
dy=0.5*(y3+y4-y1-y2) |
307 |
vVec1(i,j,bi,bj)= |
308 |
& (tmp_u(i,j,bi,bj)*dx+ |
309 |
& tmp_v(i,j,bi,bj)*dy)/ |
310 |
& SQRT(dx*dx+dy*dy) |
311 |
dx=0.5*(x2+x4-x1-x3) |
312 |
dx=dx* |
313 |
& cos(deg2rad*yC(i,j,bi,bj)) |
314 |
dy=0.5*(y2+y4-y1-y3) |
315 |
uVec1(i,j,bi,bj)= |
316 |
& (tmp_u(i,j,bi,bj)*dx+ |
317 |
& tmp_v(i,j,bi,bj)*dy)/ |
318 |
& SQRT(dx*dx+dy*dy) |
319 |
#else /* EXF_USE_OLD_VEC_ROTATION */ |
320 |
uVec1(i,j,bi,bj) = |
321 |
& angleCosC(i,j,bi,bj)*tmp_u(i,j,bi,bj) |
322 |
& +angleSinC(i,j,bi,bj)*tmp_v(i,j,bi,bj) |
323 |
vVec1(i,j,bi,bj) = |
324 |
& -angleSinC(i,j,bi,bj)*tmp_u(i,j,bi,bj) |
325 |
& +angleCosC(i,j,bi,bj)*tmp_v(i,j,bi,bj) |
326 |
#endif /* EXF_USE_OLD_VEC_ROTATION */ |
327 |
ENDDO |
328 |
ENDDO |
329 |
ENDDO |
330 |
ENDDO |
331 |
ELSE |
332 |
DO bj = myByLo(myThid),myByHi(myThid) |
333 |
DO bi = myBxLo(myThid),myBxHi(myThid) |
334 |
DO j = 1,sNy |
335 |
DO i = 1,sNx |
336 |
uVec1(i,j,bi,bj) = uVec_inScale*tmp_u(i,j,bi,bj) |
337 |
vVec1(i,j,bi,bj) = vVec_inScale*tmp_v(i,j,bi,bj) |
338 |
ENDDO |
339 |
ENDDO |
340 |
ENDDO |
341 |
ENDDO |
342 |
ENDIF |
343 |
C- apply mask (after scaling factor and rotation) |
344 |
CALL EXF_FILTER_RL( uVec1, uVecMask, myThid ) |
345 |
CALL EXF_FILTER_RL( vVec1, vVecMask, myThid ) |
346 |
|
347 |
C- end if ( first ) block |
348 |
ENDIF |
349 |
|
350 |
IF ( first .OR. changed ) THEN |
351 |
C-- Load and interpolate a new reccord |
352 |
|
353 |
CALL exf_SwapFFields( uVec0, uVec1, myThid ) |
354 |
CALL exf_SwapFFields( vVec0, vVec1, myThid ) |
355 |
|
356 |
CALL exf_GetYearlyFieldName( |
357 |
I useExfYearlyFields, twoDigitYear, uVecPeriod, year1, |
358 |
I uVecFile, |
359 |
O uVecFile1, |
360 |
I myTime, myIter, myThid ) |
361 |
CALL exf_GetYearlyFieldName( |
362 |
I useExfYearlyFields, twoDigitYear, vVecPeriod, year1, |
363 |
I vVecFile, |
364 |
O vVecFile1, |
365 |
I myTime, myIter, myThid ) |
366 |
IF ( exf_debugLev.GE.debLevC ) THEN |
367 |
_BEGIN_MASTER(myThid) |
368 |
WRITE(msgBuf,'(6A,I10)') 'EXF_SET_UV: ', |
369 |
& 'fields "', uVecName, '" & "', vVecName, '", it=', myIter |
370 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
371 |
& SQUEEZE_RIGHT, myThid ) |
372 |
j = ILNBLNK(uVecFile1) |
373 |
WRITE(msgBuf,'(2A,I6,3A)') 'EXF_SET_UV: ', |
374 |
& 'loading rec=', count1, ' from file: "', uVecFile1(1:j), '"' |
375 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
376 |
& SQUEEZE_RIGHT, myThid ) |
377 |
j = ILNBLNK(vVecFile1) |
378 |
WRITE(msgBuf,'(2A,I6,3A)') 'EXF_SET_UV: ', |
379 |
& 'loading rec=', count1, ' from file: "', vVecFile1(1:j), '"' |
380 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
381 |
& SQUEEZE_RIGHT, myThid ) |
382 |
_END_MASTER(myThid) |
383 |
ENDIF |
384 |
|
385 |
IF ( uvInterp ) THEN |
386 |
C- vector interpolation to (xC,yC) locations |
387 |
CALL EXF_INTERP_UV( |
388 |
I uVecFile1, vVecFile1, exf_iprec, count1, |
389 |
I uVec_nlon, uVec_nlat, |
390 |
I uVec_lon0, uVec_lon_inc, uVec_lat0, uVec_lat_inc, |
391 |
#ifdef EXF_INTERP_USE_DYNALLOC |
392 |
O tmp_u, tmp_v, |
393 |
#else |
394 |
O tmp_u, tmp_v, bufArr1, bufArr2, |
395 |
#endif |
396 |
I xC, yC, |
397 |
I u_interp_method, v_interp_method, myIter, myThid ) |
398 |
ELSE |
399 |
C- scalar interpolation to (xC,yC) locations |
400 |
CALL EXF_INTERP( |
401 |
I uVecFile1, exf_iprec, |
402 |
#ifdef EXF_INTERP_USE_DYNALLOC |
403 |
O tmp_u, |
404 |
#else |
405 |
O tmp_u, bufArr1, |
406 |
#endif |
407 |
I count1, xC, yC, |
408 |
I uVec_lon0, uVec_lon_inc, uVec_lat0, uVec_lat_inc, |
409 |
I uVec_nlon, uVec_nlat, u_interp_method, |
410 |
I myIter, myThid ) |
411 |
CALL EXF_INTERP( |
412 |
I vVecFile1, exf_iprec, |
413 |
#ifdef EXF_INTERP_USE_DYNALLOC |
414 |
O tmp_v, |
415 |
#else |
416 |
O tmp_v, bufArr2, |
417 |
#endif |
418 |
I count1, xC, yC, |
419 |
I vVec_lon0, vVec_lon_inc, vVec_lat0, vVec_lat_inc, |
420 |
I vVec_nlon, vVec_nlat, v_interp_method, |
421 |
I myIter, myThid ) |
422 |
ENDIF |
423 |
|
424 |
C- apply mask: Note: done after applying scaling factor and rotation |
425 |
c CALL EXF_FILTER_RL( tmp_u, uVecMask, myThid ) |
426 |
c CALL EXF_FILTER_RL( tmp_v, vVecMask, myThid ) |
427 |
|
428 |
IF ( exf_output_interp ) THEN |
429 |
j = ILNBLNK(uVecFile1) |
430 |
WRITE(out_uVecFile,'(2A)') uVecFile1(1:j), '_out' |
431 |
CALL WRITE_REC_XY_RL(out_uVecFile,tmp_u,count1,myIter,myThid) |
432 |
j = ILNBLNK(vVecFile1) |
433 |
WRITE(out_vVecFile,'(2A)') vVecFile1(1:j), '_out' |
434 |
CALL WRITE_REC_XY_RL(out_vVecFile,tmp_v,count1,myIter,myThid) |
435 |
ENDIF |
436 |
|
437 |
C- scaling factor and vector rotation |
438 |
IF ( usingCurvilinearGrid .OR. rotateGrid ) THEN |
439 |
DO bj = myByLo(myThid),myByHi(myThid) |
440 |
DO bi = myBxLo(myThid),myBxHi(myThid) |
441 |
DO j = 1,sNy |
442 |
DO i = 1,sNx |
443 |
tmp_u(i,j,bi,bj) = uVec_inScale*tmp_u(i,j,bi,bj) |
444 |
tmp_v(i,j,bi,bj) = vVec_inScale*tmp_v(i,j,bi,bj) |
445 |
ENDDO |
446 |
ENDDO |
447 |
DO j = 1,sNy |
448 |
DO i = 1,sNx |
449 |
#ifdef EXF_USE_OLD_VEC_ROTATION |
450 |
x1=xG(i,j,bi,bj) |
451 |
x2=xG(i+1,j,bi,bj) |
452 |
x3=xG(i,j+1,bi,bj) |
453 |
x4=xG(i+1,j+1,bi,bj) |
454 |
IF ((x2-x1).GT.180) x2=x2-360 |
455 |
IF ((x1-x2).GT.180) x2=x2+360 |
456 |
IF ((x3-x1).GT.180) x3=x3-360 |
457 |
IF ((x1-x3).GT.180) x3=x3+360 |
458 |
IF ((x4-x1).GT.180) x4=x4-360 |
459 |
IF ((x1-x4).GT.180) x4=x4+360 |
460 |
y1=yG(i,j,bi,bj) |
461 |
y2=yG(i+1,j,bi,bj) |
462 |
y3=yG(i,j+1,bi,bj) |
463 |
y4=yG(i+1,j+1,bi,bj) |
464 |
dx=0.5*(x3+x4-x1-x2) |
465 |
dx=dx* |
466 |
& cos(deg2rad*yC(i,j,bi,bj)) |
467 |
dy=0.5*(y3+y4-y1-y2) |
468 |
vVec1(i,j,bi,bj)= |
469 |
& (tmp_u(i,j,bi,bj)*dx+ |
470 |
& tmp_v(i,j,bi,bj)*dy)/ |
471 |
& SQRT(dx*dx+dy*dy) |
472 |
dx=0.5*(x2+x4-x1-x3) |
473 |
dx=dx* |
474 |
& cos(deg2rad*yC(i,j,bi,bj)) |
475 |
dy=0.5*(y2+y4-y1-y3) |
476 |
uVec1(i,j,bi,bj)= |
477 |
& (tmp_u(i,j,bi,bj)*dx+ |
478 |
& tmp_v(i,j,bi,bj)*dy)/ |
479 |
& SQRT(dx*dx+dy*dy) |
480 |
#else /* EXF_USE_OLD_VEC_ROTATION */ |
481 |
uVec1(i,j,bi,bj) = |
482 |
& angleCosC(i,j,bi,bj)*tmp_u(i,j,bi,bj) |
483 |
& +angleSinC(i,j,bi,bj)*tmp_v(i,j,bi,bj) |
484 |
vVec1(i,j,bi,bj) = |
485 |
& -angleSinC(i,j,bi,bj)*tmp_u(i,j,bi,bj) |
486 |
& +angleCosC(i,j,bi,bj)*tmp_v(i,j,bi,bj) |
487 |
#endif /* EXF_USE_OLD_VEC_ROTATION */ |
488 |
ENDDO |
489 |
ENDDO |
490 |
ENDDO |
491 |
ENDDO |
492 |
ELSE |
493 |
DO bj = myByLo(myThid),myByHi(myThid) |
494 |
DO bi = myBxLo(myThid),myBxHi(myThid) |
495 |
DO j = 1,sNy |
496 |
DO i = 1,sNx |
497 |
uVec1(i,j,bi,bj) = uVec_inScale*tmp_u(i,j,bi,bj) |
498 |
vVec1(i,j,bi,bj) = vVec_inScale*tmp_v(i,j,bi,bj) |
499 |
ENDDO |
500 |
ENDDO |
501 |
ENDDO |
502 |
ENDDO |
503 |
ENDIF |
504 |
C- apply mask (after scaling factor and rotation) |
505 |
CALL EXF_FILTER_RL( uVec1, uVecMask, myThid ) |
506 |
CALL EXF_FILTER_RL( vVec1, vVecMask, myThid ) |
507 |
|
508 |
C- end if ( first or changed ) block |
509 |
ENDIF |
510 |
|
511 |
C-- Interpolate linearly onto the current time. |
512 |
DO bj = myByLo(myThid),myByHi(myThid) |
513 |
DO bi = myBxLo(myThid),myBxHi(myThid) |
514 |
DO j = 1,sNy |
515 |
DO i = 1,sNx |
516 |
uVec(i,j,bi,bj) = fac * uVec0(i,j,bi,bj) |
517 |
& + (exf_one - fac)* uVec1(i,j,bi,bj) |
518 |
vVec(i,j,bi,bj) = fac * vVec0(i,j,bi,bj) |
519 |
& + (exf_one - fac)* vVec1(i,j,bi,bj) |
520 |
ENDDO |
521 |
ENDDO |
522 |
ENDDO |
523 |
ENDDO |
524 |
|
525 |
ELSE |
526 |
C case no-interpolation |
527 |
C or ( .NOT.usingCurvilinearGrid & .NOT.rotateGrid & .NOT.uvInterp ) |
528 |
#else /* USE_EXF_INTERPOLATION */ |
529 |
IF ( .TRUE. ) THEN |
530 |
#endif /* USE_EXF_INTERPOLATION */ |
531 |
|
532 |
CALL EXF_SET_FLD( |
533 |
I uVecName, uVecFile, uVecMask, |
534 |
I uVecStartTime, uVecPeriod, uVecRepeatCycle, |
535 |
I uVec_inScale, uVec_remove_intercept, uVec_remove_slope, |
536 |
U uVec, uVec0, uVec1, |
537 |
#ifdef USE_EXF_INTERPOLATION |
538 |
I uVec_lon0, uVec_lon_inc, uVec_lat0, uVec_lat_inc, |
539 |
I uVec_nlon, uVec_nlat, xC, yC, u_interp_method, |
540 |
#endif /* USE_EXF_INTERPOLATION */ |
541 |
I myTime, myIter, myThid ) |
542 |
|
543 |
CALL EXF_SET_FLD( |
544 |
I vVecName, vVecFile, vVecMask, |
545 |
I vVecStartTime, vVecPeriod, vVecRepeatCycle, |
546 |
I vVec_inScale, vVec_remove_intercept, vVec_remove_slope, |
547 |
U vVec, vVec0, vVec1, |
548 |
#ifdef USE_EXF_INTERPOLATION |
549 |
I vVec_lon0, vVec_lon_inc, vVec_lat0, vVec_lat_inc, |
550 |
I vVec_nlon, vVec_nlat, xC, yC, v_interp_method, |
551 |
#endif /* USE_EXF_INTERPOLATION */ |
552 |
I myTime, myIter, myThid ) |
553 |
|
554 |
C- vector rotation |
555 |
IF ( rotateStressOnAgrid ) THEN |
556 |
DO bj = myByLo(myThid),myByHi(myThid) |
557 |
DO bi = myBxLo(myThid),myBxHi(myThid) |
558 |
DO j = 1,sNy |
559 |
DO i = 1,sNx |
560 |
tmp_u(i,j,bi,bj) = uVec(i,j,bi,bj) |
561 |
tmp_v(i,j,bi,bj) = vVec(i,j,bi,bj) |
562 |
ENDDO |
563 |
ENDDO |
564 |
ENDDO |
565 |
ENDDO |
566 |
DO bj = myByLo(myThid),myByHi(myThid) |
567 |
DO bi = myBxLo(myThid),myBxHi(myThid) |
568 |
DO j = 1,sNy |
569 |
DO i = 1,sNx |
570 |
uVec(i,j,bi,bj) = |
571 |
& angleCosC(i,j,bi,bj)*tmp_u(i,j,bi,bj) |
572 |
& +angleSinC(i,j,bi,bj)*tmp_v(i,j,bi,bj) |
573 |
vVec(i,j,bi,bj) = |
574 |
& -angleSinC(i,j,bi,bj)*tmp_u(i,j,bi,bj) |
575 |
& +angleCosC(i,j,bi,bj)*tmp_v(i,j,bi,bj) |
576 |
ENDDO |
577 |
ENDDO |
578 |
ENDDO |
579 |
ENDDO |
580 |
ENDIF |
581 |
|
582 |
ENDIF |
583 |
|
584 |
RETURN |
585 |
END |