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

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

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


Revision 1.29 - (show annotations) (download)
Sat Jun 7 20:16:58 2014 UTC (10 years, 6 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64y, checkpoint64z, checkpoint65, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65o, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e
Changes since 1.28: +3 -3 lines
minor adjustment

1 C $Header: /u/gcmpack/MITgcm/pkg/exf/exf_set_uv.F,v 1.28 2014/06/07 17:52:18 jmc Exp $
2 C $Name: $
3
4 #include "EXF_OPTIONS.h"
5
6 SUBROUTINE EXF_SET_UV(
7 I uvecfile, uvecstartdate, uvecperiod,
8 I exf_inscal_uvec, uvec_remove_intercept, uvec_remove_slope,
9 U uvec, uvec0, uvec1, uvecmask,
10 I vvecfile, vvecstartdate, vvecperiod,
11 I exf_inscal_vvec, vvec_remove_intercept, vvec_remove_slope,
12 U vvec, vvec0, vvec1, vvecmask,
13 #ifdef USE_EXF_INTERPOLATION
14 I uvec_lon0, uvec_lon_inc, uvec_lat0, uvec_lat_inc,
15 I uvec_nlon, uvec_nlat, u_interp_method,
16 I vvec_lon0, vvec_lon_inc, vvec_lat0, vvec_lat_inc,
17 I vvec_nlon, vvec_nlat, v_interp_method, uvInterp,
18 #endif /* USE_EXF_INTERPOLATION */
19 I myTime, myIter, 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 #include "EEPARAMS.h"
38 #include "SIZE.h"
39 #include "PARAMS.h"
40 #include "GRID.h"
41 #include "EXF_PARAM.h"
42 #include "EXF_FIELDS.h"
43 #include "EXF_CONSTANTS.h"
44
45 C == routine arguments ==
46 C *vec_lon_0, :: longitude and latitude of SouthWest
47 C *vec_lat_0 corner of global input grid for *vec
48 C *vec_nlon, *vec_nlat :: input x-grid and y-grid size for *vec
49 C *vec_lon_inc :: scalar x-grid increment for *vec
50 C *vec_lat_inc :: vector y-grid increments for *vec
51
52 CHARACTER*(128) uvecfile
53 _RL uvecstartdate, uvecperiod
54 _RL exf_inscal_uvec
55 _RL uvec_remove_intercept, uvec_remove_slope
56 _RL uvec (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
57 _RL uvec0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
58 _RL uvec1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
59 CHARACTER*1 uvecmask
60 CHARACTER*(128) vvecfile
61 _RL vvecstartdate, vvecperiod
62 _RL exf_inscal_vvec
63 _RL vvec_remove_intercept, vvec_remove_slope
64 _RL vvec (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
65 _RL vvec0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
66 _RL vvec1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
67 CHARACTER*1 vvecmask
68 #ifdef USE_EXF_INTERPOLATION
69 _RL uvec_lon0, uvec_lon_inc
70 _RL uvec_lat0, uvec_lat_inc(MAX_LAT_INC)
71 INTEGER uvec_nlon, uvec_nlat, u_interp_method
72 _RL vvec_lon0, vvec_lon_inc
73 _RL vvec_lat0, vvec_lat_inc(MAX_LAT_INC)
74 INTEGER vvec_nlon, vvec_nlat, v_interp_method
75 LOGICAL uvInterp
76 #endif /* USE_EXF_INTERPOLATION */
77 _RL myTime
78 INTEGER myIter
79 INTEGER myThid
80
81 C == Functions ==
82 INTEGER ILNBLNK
83 EXTERNAL ILNBLNK
84
85 C == local variables ==
86 #ifdef USE_EXF_INTERPOLATION
87 C msgBuf :: Informational/error message buffer
88 CHARACTER*(MAX_LEN_MBUF) msgBuf
89 CHARACTER*(128) uvecfile0, uvecfile1
90 CHARACTER*(128) vvecfile0, vvecfile1
91 LOGICAL first, changed
92 _RL fac
93 #ifdef EXF_USE_OLD_VEC_ROTATION
94 _RL x1, x2, x3, x4, y1, y2, y3, y4, dx, dy
95 #endif
96 _RL tmp_u (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
97 _RL tmp_v (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
98 INTEGER count0, count1
99 INTEGER i, j, bi, bj
100 INTEGER year0, year1
101 #endif /* USE_EXF_INTERPOLATION */
102
103 C == end of interface ==
104
105 #ifdef USE_EXF_INTERPOLATION
106 IF ( u_interp_method.GE.1 .AND. v_interp_method.GE.1 .AND.
107 & uvecfile.NE.' ' .AND. vvecfile.NE.' ' .AND.
108 & (usingCurvilinearGrid .OR. rotateGrid .OR. uvInterp) ) THEN
109
110 IF ( uvecperiod .EQ. -12. ) THEN
111 C- genperiod=-12 means input file contains 12 monthly means
112 C records, corresponding to Jan. (rec=1) through Dec. (rec=12)
113 CALL cal_GetMonthsRec(
114 O fac, first, changed,
115 O count0, count1,
116 I myTime, myIter, myThid )
117
118 ELSEIF ( uvecperiod .LE. 0. ) THEN
119 j = ILNBLNK(uvecfile)
120 WRITE(msgBuf,'(A,1PE16.8,2A)')
121 & 'EXF_SET_UV: Invalid uvecperiod=', uvecperiod,
122 & ' for file: ', uvecfile(1:j)
123 CALL PRINT_ERROR( msgBuf, myThid )
124 STOP 'ABNORMAL END: S/R EXF_SET_UV'
125 ELSE
126 C- get record numbers and interpolation factor
127 CALL exf_GetFFieldRec(
128 I uvecstartdate, uvecperiod,
129 I useExfYearlyFields,
130 O fac, first, changed,
131 O count0, count1, year0, year1,
132 I myTime, myIter, myThid )
133 ENDIF
134 IF ( exf_debugLev.GE.debLevD ) THEN
135 _BEGIN_MASTER( myThid )
136 i = ILNBLNK(uvecfile)
137 j = ILNBLNK(vvecfile)
138 WRITE(msgBuf,'(5A)') ' EXF_SET_UV: ',
139 & 'processing: ', uvecfile(1:i), ' & ', vvecfile(1:j)
140 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
141 & SQUEEZE_RIGHT, myThid )
142 WRITE(msgBuf,'(2A,I10,2I7)') ' EXF_SET_UV: ',
143 & ' myIter, count0, count1:', myIter, count0, count1
144 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
145 & SQUEEZE_RIGHT, myThid )
146 WRITE(msgBuf,'(2A,2(L2,2X),E16.9)') ' EXF_SET_UV: ',
147 & ' first, changed, fac: ', first, changed, fac
148 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
149 & SQUEEZE_RIGHT, myThid )
150 _END_MASTER( myThid )
151 ENDIF
152
153 IF ( first ) THEN
154 C-- Load and interpolate a new reccord (= 1rst one of this run)
155
156 CALL exf_GetYearlyFieldName(
157 I useExfYearlyFields, twoDigitYear, uvecperiod, year0,
158 I uvecfile,
159 O uvecfile0,
160 I myTime, myIter, myThid )
161 CALL exf_GetYearlyFieldName(
162 I useExfYearlyFields, twoDigitYear, vvecperiod, year0,
163 I vvecfile,
164 O vvecfile0,
165 I myTime, myIter, myThid )
166 IF ( exf_debugLev.GE.debLevC ) THEN
167 _BEGIN_MASTER(myThid)
168 j = ILNBLNK(uvecfile0)
169 WRITE(msgBuf,'(A,I10,A,I6,2A)')
170 & ' EXF_SET_UV: it=', myIter, ' loading rec=', count0,
171 & ' from: ', uvecfile0(1:j)
172 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
173 & SQUEEZE_RIGHT, myThid )
174 j = ILNBLNK(vvecfile0)
175 WRITE(msgBuf,'(A,I10,A,I6,2A)')
176 & ' EXF_SET_UV: it=', myIter, ' loading rec=', count0,
177 & ' from: ', vvecfile0(1:j)
178 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
179 & SQUEEZE_RIGHT, myThid )
180 _END_MASTER(myThid)
181 ENDIF
182
183 IF ( uvInterp ) THEN
184 C- vector interpolation to (xC,yC) locations
185 CALL EXF_INTERP_UV(
186 I uvecfile0, vvecfile0, exf_iprec, count0,
187 I uvec_nlon, uvec_nlat,
188 I uvec_lon0, uvec_lon_inc, uvec_lat0, uvec_lat_inc,
189 O tmp_u, tmp_v,
190 I xC, yC,
191 I u_interp_method, v_interp_method, myIter, myThid )
192 ELSE
193 C- scalar interpolation to (xC,yC) locations
194 CALL EXF_INTERP(
195 I uvecfile0, exf_iprec,
196 O tmp_u,
197 I count0, xC, yC,
198 I uvec_lon0, uvec_lon_inc, uvec_lat0, uvec_lat_inc,
199 I uvec_nlon, uvec_nlat, u_interp_method,
200 I myIter, myThid )
201 CALL EXF_INTERP(
202 I vvecfile0, exf_iprec,
203 O tmp_v,
204 I count0, xC, yC,
205 I vvec_lon0, vvec_lon_inc, vvec_lat0, vvec_lat_inc,
206 I vvec_nlon, vvec_nlat, v_interp_method,
207 I myIter, myThid )
208 ENDIF
209
210 C- vector rotation
211 IF ( usingCurvilinearGrid .OR. rotateGrid ) THEN
212 DO bj = myByLo(myThid),myByHi(myThid)
213 DO bi = myBxLo(myThid),myBxHi(myThid)
214 DO j = 1,sNy
215 DO i = 1,sNx
216 #ifdef EXF_USE_OLD_VEC_ROTATION
217 x1=xG(i,j,bi,bj)
218 x2=xG(i+1,j,bi,bj)
219 x3=xG(i,j+1,bi,bj)
220 x4=xG(i+1,j+1,bi,bj)
221 IF ((x2-x1).GT.180) x2=x2-360
222 IF ((x1-x2).GT.180) x2=x2+360
223 IF ((x3-x1).GT.180) x3=x3-360
224 IF ((x1-x3).GT.180) x3=x3+360
225 IF ((x4-x1).GT.180) x4=x4-360
226 IF ((x1-x4).GT.180) x4=x4+360
227 y1=yG(i,j,bi,bj)
228 y2=yG(i+1,j,bi,bj)
229 y3=yG(i,j+1,bi,bj)
230 y4=yG(i+1,j+1,bi,bj)
231 dx=0.5*(x3+x4-x1-x2)
232 dx=dx*
233 & cos(deg2rad*yC(i,j,bi,bj))
234 dy=0.5*(y3+y4-y1-y2)
235 vvec1(i,j,bi,bj)=
236 & (tmp_u(i,j,bi,bj)*dx+
237 & tmp_v(i,j,bi,bj)*dy)/
238 & SQRT(dx*dx+dy*dy)
239 dx=0.5*(x2+x4-x1-x3)
240 dx=dx*
241 & cos(deg2rad*yC(i,j,bi,bj))
242 dy=0.5*(y2+y4-y1-y3)
243 uvec1(i,j,bi,bj)=
244 & (tmp_u(i,j,bi,bj)*dx+
245 & tmp_v(i,j,bi,bj)*dy)/
246 & SQRT(dx*dx+dy*dy)
247 #else /* EXF_USE_OLD_VEC_ROTATION */
248 uvec1(i,j,bi,bj) =
249 & angleCosC(i,j,bi,bj)*tmp_u(i,j,bi,bj)
250 & +angleSinC(i,j,bi,bj)*tmp_v(i,j,bi,bj)
251 vvec1(i,j,bi,bj) =
252 & -angleSinC(i,j,bi,bj)*tmp_u(i,j,bi,bj)
253 & +angleCosC(i,j,bi,bj)*tmp_v(i,j,bi,bj)
254 #endif /* EXF_USE_OLD_VEC_ROTATION */
255 ENDDO
256 ENDDO
257 ENDDO
258 ENDDO
259 ELSE
260 DO bj = myByLo(myThid),myByHi(myThid)
261 DO bi = myBxLo(myThid),myBxHi(myThid)
262 DO j = 1,sNy
263 DO i = 1,sNx
264 uvec1(i,j,bi,bj) = tmp_u(i,j,bi,bj)
265 vvec1(i,j,bi,bj) = tmp_v(i,j,bi,bj)
266 ENDDO
267 ENDDO
268 ENDDO
269 ENDDO
270 ENDIF
271 C- apply mask
272 CALL EXF_FILTER_RL( uvec1, uvecmask, myThid )
273 CALL EXF_FILTER_RL( vvec1, vvecmask, myThid )
274 ENDIF
275
276 IF ( first .OR. changed ) THEN
277 C-- Load and interpolate a new reccord
278
279 CALL exf_SwapFFields( uvec0, uvec1, myThid )
280 CALL exf_SwapFFields( vvec0, vvec1, myThid )
281
282 CALL exf_GetYearlyFieldName(
283 I useExfYearlyFields, twoDigitYear, uvecperiod, year1,
284 I uvecfile,
285 O uvecfile1,
286 I myTime, myIter, myThid )
287 CALL exf_GetYearlyFieldName(
288 I useExfYearlyFields, twoDigitYear, vvecperiod, year1,
289 I vvecfile,
290 O vvecfile1,
291 I myTime, myIter, myThid )
292 IF ( exf_debugLev.GE.debLevC ) THEN
293 _BEGIN_MASTER(myThid)
294 j = ILNBLNK(uvecfile1)
295 WRITE(msgBuf,'(A,I10,A,I6,2A)')
296 & ' EXF_SET_UV: it=', myIter, ' loading rec=', count1,
297 & ' from: ', uvecfile1(1:j)
298 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
299 & SQUEEZE_RIGHT, myThid )
300 j = ILNBLNK(vvecfile1)
301 WRITE(msgBuf,'(A,I10,A,I6,2A)')
302 & ' EXF_SET_UV: it=', myIter, ' loading rec=', count1,
303 & ' from: ', vvecfile1(1:j)
304 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
305 & SQUEEZE_RIGHT, myThid )
306 _END_MASTER(myThid)
307 ENDIF
308
309 IF ( uvInterp ) THEN
310 C- vector interpolation to (xC,yC) locations
311 CALL EXF_INTERP_UV(
312 I uvecfile1, vvecfile1, exf_iprec, count1,
313 I uvec_nlon, uvec_nlat,
314 I uvec_lon0, uvec_lon_inc, uvec_lat0, uvec_lat_inc,
315 O tmp_u, tmp_v,
316 I xC, yC,
317 I u_interp_method, v_interp_method, myIter, myThid )
318 ELSE
319 C- scalar interpolation to (xC,yC) locations
320 CALL EXF_INTERP(
321 I uvecfile1, exf_iprec,
322 O tmp_u,
323 I count1, xC, yC,
324 I uvec_lon0, uvec_lon_inc, uvec_lat0, uvec_lat_inc,
325 I uvec_nlon, uvec_nlat, u_interp_method,
326 I myIter, myThid )
327 CALL EXF_INTERP(
328 I vvecfile1, exf_iprec,
329 O tmp_v,
330 I count1, xC, yC,
331 I vvec_lon0, vvec_lon_inc, vvec_lat0, vvec_lat_inc,
332 I vvec_nlon, vvec_nlat, v_interp_method,
333 I myIter, myThid )
334 ENDIF
335
336 C- vector rotation
337 IF ( usingCurvilinearGrid .OR. rotateGrid ) THEN
338 DO bj = myByLo(myThid),myByHi(myThid)
339 DO bi = myBxLo(myThid),myBxHi(myThid)
340 DO j = 1,sNy
341 DO i = 1,sNx
342 #ifdef EXF_USE_OLD_VEC_ROTATION
343 x1=xG(i,j,bi,bj)
344 x2=xG(i+1,j,bi,bj)
345 x3=xG(i,j+1,bi,bj)
346 x4=xG(i+1,j+1,bi,bj)
347 IF ((x2-x1).GT.180) x2=x2-360
348 IF ((x1-x2).GT.180) x2=x2+360
349 IF ((x3-x1).GT.180) x3=x3-360
350 IF ((x1-x3).GT.180) x3=x3+360
351 IF ((x4-x1).GT.180) x4=x4-360
352 IF ((x1-x4).GT.180) x4=x4+360
353 y1=yG(i,j,bi,bj)
354 y2=yG(i+1,j,bi,bj)
355 y3=yG(i,j+1,bi,bj)
356 y4=yG(i+1,j+1,bi,bj)
357 dx=0.5*(x3+x4-x1-x2)
358 dx=dx*
359 & cos(deg2rad*yC(i,j,bi,bj))
360 dy=0.5*(y3+y4-y1-y2)
361 vvec1(i,j,bi,bj)=
362 & (tmp_u(i,j,bi,bj)*dx+
363 & tmp_v(i,j,bi,bj)*dy)/
364 & SQRT(dx*dx+dy*dy)
365 dx=0.5*(x2+x4-x1-x3)
366 dx=dx*
367 & cos(deg2rad*yC(i,j,bi,bj))
368 dy=0.5*(y2+y4-y1-y3)
369 uvec1(i,j,bi,bj)=
370 & (tmp_u(i,j,bi,bj)*dx+
371 & tmp_v(i,j,bi,bj)*dy)/
372 & SQRT(dx*dx+dy*dy)
373 #else /* EXF_USE_OLD_VEC_ROTATION */
374 uvec1(i,j,bi,bj) =
375 & angleCosC(i,j,bi,bj)*tmp_u(i,j,bi,bj)
376 & +angleSinC(i,j,bi,bj)*tmp_v(i,j,bi,bj)
377 vvec1(i,j,bi,bj) =
378 & -angleSinC(i,j,bi,bj)*tmp_u(i,j,bi,bj)
379 & +angleCosC(i,j,bi,bj)*tmp_v(i,j,bi,bj)
380 #endif /* EXF_USE_OLD_VEC_ROTATION */
381 ENDDO
382 ENDDO
383 ENDDO
384 ENDDO
385 ELSE
386 DO bj = myByLo(myThid),myByHi(myThid)
387 DO bi = myBxLo(myThid),myBxHi(myThid)
388 DO j = 1,sNy
389 DO i = 1,sNx
390 uvec1(i,j,bi,bj) = tmp_u(i,j,bi,bj)
391 vvec1(i,j,bi,bj) = tmp_v(i,j,bi,bj)
392 ENDDO
393 ENDDO
394 ENDDO
395 ENDDO
396 ENDIF
397 C- apply mask
398 CALL EXF_FILTER_RL( uvec1, uvecmask, myThid )
399 CALL EXF_FILTER_RL( vvec1, vvecmask, myThid )
400 ENDIF
401
402 C-- Interpolate linearly onto the current time.
403 DO bj = myByLo(myThid),myByHi(myThid)
404 DO bi = myBxLo(myThid),myBxHi(myThid)
405 DO j = 1,sNy
406 DO i = 1,sNx
407 uvec(i,j,bi,bj) = exf_inscal_uvec * (
408 & fac * uvec0(i,j,bi,bj) +
409 & (exf_one - fac) * uvec1(i,j,bi,bj) )
410 vvec(i,j,bi,bj) = exf_inscal_vvec * (
411 & fac * vvec0(i,j,bi,bj) +
412 & (exf_one - fac) * vvec1(i,j,bi,bj) )
413 ENDDO
414 ENDDO
415 ENDDO
416 ENDDO
417
418 ELSE
419 C case no-interpolation
420 C or ( .NOT.usingCurvilinearGrid & .NOT.rotateGrid & .NOT.uvInterp )
421 #else /* USE_EXF_INTERPOLATION */
422 IF ( .TRUE. ) THEN
423 #endif /* USE_EXF_INTERPOLATION */
424
425 CALL EXF_SET_GEN(
426 & uvecfile, uvecstartdate, uvecperiod,
427 & exf_inscal_uvec,
428 & uvec_remove_intercept, uvec_remove_slope,
429 & uvec, uvec0, uvec1, uvecmask,
430 #ifdef USE_EXF_INTERPOLATION
431 & uvec_lon0, uvec_lon_inc, uvec_lat0, uvec_lat_inc,
432 & uvec_nlon, uvec_nlat, xC, yC, u_interp_method,
433 #endif /* USE_EXF_INTERPOLATION */
434 & myTime, myIter, myThid )
435
436 CALL EXF_SET_GEN(
437 & vvecfile, vvecstartdate, vvecperiod,
438 & exf_inscal_vvec,
439 & vvec_remove_intercept, vvec_remove_slope,
440 & vvec, vvec0, vvec1, vvecmask,
441 #ifdef USE_EXF_INTERPOLATION
442 & vvec_lon0, vvec_lon_inc, vvec_lat0, vvec_lat_inc,
443 & vvec_nlon, vvec_nlat, xC, yC, v_interp_method,
444 #endif /* USE_EXF_INTERPOLATION */
445 & myTime, myIter, myThid )
446
447 ENDIF
448
449 RETURN
450 END

  ViewVC Help
Powered by ViewVC 1.1.22