/[MITgcm]/MITgcm/model/src/apply_forcing.F
ViewVC logotype

Contents of /MITgcm/model/src/apply_forcing.F

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

Revision 1.6 - (show annotations) (download)
Fri Sep 5 21:04:34 2014 UTC (9 years, 9 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint65h, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e
Changes since 1.5: +58 -76 lines 
- change hard-coded loop-range in S/R APPLY_FORCING_T/S to set T & S
  forcing over 0:sNx+1, 0:sNy=1 (instead of just 1:sNx,1:sNy): Model
  dynamics requires valid tracer forcing only over tile interior but, for
  some diagnostics, may need to extend it over 1 point in tile halo region.
1 C $Header: /u/gcmpack/MITgcm/model/src/apply_forcing.F,v 1.5 2014/08/07 18:43:33 heimbach Exp $
2 C $Name: $
3
4 #include "PACKAGES_CONFIG.h"
5 #include "CPP_OPTIONS.h"
6
7 C-- File apply_forcing.F:
8 C-- Contents
9 C-- o APPLY_FORCING_U
10 C-- o APPLY_FORCING_V
11 C-- o APPLY_FORCING_T
12 C-- o APPLY_FORCING_S
13
14 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
15 CBOP
16 C !ROUTINE: APPLY_FORCING_U
17 C !INTERFACE:
18 SUBROUTINE APPLY_FORCING_U(
19 U gU_arr,
20 I iMin,iMax,jMin,jMax, k, bi, bj,
21 I myTime, myIter, myThid )
22 C !DESCRIPTION: \bv
23 C *==========================================================*
24 C | S/R APPLY_FORCING_U
25 C | o Contains problem specific forcing for zonal velocity.
26 C *==========================================================*
27 C | Adds terms to gU for forcing by external sources
28 C | e.g. wind stress, bottom friction etc ...
29 C *==========================================================*
30 C \ev
31
32 C !USES:
33 IMPLICIT NONE
34 C == Global data ==
35 #include "SIZE.h"
36 #include "EEPARAMS.h"
37 #include "PARAMS.h"
38 #include "GRID.h"
39 #include "DYNVARS.h"
40 #include "FFIELDS.h"
41
42 C !INPUT/OUTPUT PARAMETERS:
43 C gU_arr :: the tendency array
44 C iMin,iMax :: Working range of x-index for applying forcing.
45 C jMin,jMax :: Working range of y-index for applying forcing.
46 C k :: Current vertical level index
47 C bi,bj :: Current tile indices
48 C myTime :: Current time in simulation
49 C myIter :: Current iteration number
50 C myThid :: my Thread Id number
51 _RL gU_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
52 INTEGER iMin, iMax, jMin, jMax
53 INTEGER k, bi, bj
54 _RL myTime
55 INTEGER myIter
56 INTEGER myThid
57
58 C !LOCAL VARIABLES:
59 C i,j :: Loop counters
60 C kSurface :: index of surface level
61 INTEGER i, j
62 #ifdef USE_OLD_EXTERNAL_FORCING
63 _RL locVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
64 _RL tmpVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65 #else
66 INTEGER kSurface
67 #endif /* USE_OLD_EXTERNAL_FORCING */
68 CEOP
69
70 #ifdef USE_OLD_EXTERNAL_FORCING
71
72 DO j=1-OLy,sNy+OLy
73 DO i=1-OLx,sNx+OLx
74 locVar(i,j) = gU(i,j,k,bi,bj)
75 ENDDO
76 ENDDO
77 CALL EXTERNAL_FORCING_U(
78 I iMin, iMax, jMin, jMax, bi, bj, k,
79 I myTime, myThid )
80 C- Use 2-d local array tmpVar and split loop in 2 parts to avoid compiler
81 C to mess-up this part by re-arranging the order of instructions (wrong
82 C when gU and gU_arr are the same array, i.e., called with argument gU).
83 DO j=1-OLy,sNy+OLy
84 DO i=1-OLx,sNx+OLx
85 tmpVar(i,j) = gU(i,j,k,bi,bj) - locVar(i,j)
86 gU(i,j,k,bi,bj) = locVar(i,j)
87 ENDDO
88 ENDDO
89 C- not needed since APPLY_FORCING_U is no longer called with argument gU
90 c CALL FOOL_THE_COMPILER_RL( tmpVar(1,1) )
91 DO j=1-OLy,sNy+OLy
92 DO i=1-OLx,sNx+OLx
93 gU_arr(i,j) = gU_arr(i,j) + tmpVar(i,j)
94 ENDDO
95 ENDDO
96
97 #else /* USE_OLD_EXTERNAL_FORCING */
98
99 IF ( fluidIsAir ) THEN
100 kSurface = 0
101 ELSEIF ( usingPCoords ) THEN
102 kSurface = Nr
103 ELSE
104 kSurface = 1
105 ENDIF
106
107 C-- Forcing term
108 #ifdef ALLOW_AIM
109 IF ( useAIM ) CALL AIM_TENDENCY_APPLY_U(
110 U gU_arr,
111 I iMin,iMax,jMin,jMax, k, bi,bj,
112 I myTime, 0, myThid )
113 #endif /* ALLOW_AIM */
114
115 #ifdef ALLOW_ATM_PHYS
116 IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_U(
117 U gU_arr,
118 I iMin,iMax,jMin,jMax, k, bi,bj,
119 I myTime, 0, myThid )
120 #endif /* ALLOW_ATM_PHYS */
121
122 #ifdef ALLOW_FIZHI
123 IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_U(
124 U gU_arr,
125 I iMin,iMax,jMin,jMax, k, bi,bj,
126 I myTime, 0, myThid )
127 #endif /* ALLOW_FIZHI */
128
129 C Ocean: Add momentum surface forcing (e.g., wind-stress) in surface level
130 IF ( k .EQ. kSurface ) THEN
131 c DO j=1,sNy
132 C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNy+1]
133 DO j=0,sNy+1
134 DO i=1,sNx+1
135 gU_arr(i,j) = gU_arr(i,j)
136 & +foFacMom*surfaceForcingU(i,j,bi,bj)
137 & *recip_drF(k)*_recip_hFacW(i,j,k,bi,bj)
138 ENDDO
139 ENDDO
140 ELSEIF ( kSurface.EQ.-1 ) THEN
141 DO j=0,sNy+1
142 DO i=1,sNx+1
143 IF ( kSurfW(i,j,bi,bj).EQ.k ) THEN
144 gU_arr(i,j) = gU_arr(i,j)
145 & +foFacMom*surfaceForcingU(i,j,bi,bj)
146 & *recip_drF(k)*_recip_hFacW(i,j,k,bi,bj)
147 ENDIF
148 ENDDO
149 ENDDO
150 ENDIF
151
152 #ifdef ALLOW_EDDYPSI
153 CALL TAUEDDY_TENDENCY_APPLY_U(
154 U gU_arr,
155 I iMin,iMax,jMin,jMax, k, bi,bj,
156 I myTime, 0, myThid )
157 #endif
158
159 #ifdef ALLOW_RBCS
160 IF (useRBCS) THEN
161 CALL RBCS_ADD_TENDENCY(
162 U gU_arr,
163 I k, bi, bj, -1,
164 I myTime, 0, myThid )
165
166 ENDIF
167 #endif /* ALLOW_RBCS */
168
169 #ifdef ALLOW_OBCS
170 IF (useOBCS) THEN
171 CALL OBCS_SPONGE_U(
172 U gU_arr,
173 I iMin,iMax,jMin,jMax, k, bi,bj,
174 I myTime, 0, myThid )
175 ENDIF
176 #endif /* ALLOW_OBCS */
177
178 #ifdef ALLOW_MYPACKAGE
179 IF ( useMYPACKAGE ) THEN
180 CALL MYPACKAGE_TENDENCY_APPLY_U(
181 U gU_arr,
182 I iMin,iMax,jMin,jMax, k, bi,bj,
183 I myTime, 0, myThid )
184 ENDIF
185 #endif /* ALLOW_MYPACKAGE */
186
187 #endif /* USE_OLD_EXTERNAL_FORCING */
188
189 RETURN
190 END
191
192 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
193 CBOP
194 C !ROUTINE: APPLY_FORCING_V
195 C !INTERFACE:
196 SUBROUTINE APPLY_FORCING_V(
197 U gV_arr,
198 I iMin,iMax,jMin,jMax, k, bi, bj,
199 I myTime, myIter, myThid )
200 C !DESCRIPTION: \bv
201 C *==========================================================*
202 C | S/R APPLY_FORCING_V
203 C | o Contains problem specific forcing for merid velocity.
204 C *==========================================================*
205 C | Adds terms to gV for forcing by external sources
206 C | e.g. wind stress, bottom friction etc ...
207 C *==========================================================*
208 C \ev
209
210 C !USES:
211 IMPLICIT NONE
212 C == Global data ==
213 #include "SIZE.h"
214 #include "EEPARAMS.h"
215 #include "PARAMS.h"
216 #include "GRID.h"
217 #include "DYNVARS.h"
218 #include "FFIELDS.h"
219
220 C !INPUT/OUTPUT PARAMETERS:
221 C gV_arr :: the tendency array
222 C iMin,iMax :: Working range of x-index for applying forcing.
223 C jMin,jMax :: Working range of y-index for applying forcing.
224 C k :: Current vertical level index
225 C bi,bj :: Current tile indices
226 C myTime :: Current time in simulation
227 C myIter :: Current iteration number
228 C myThid :: my Thread Id number
229 _RL gV_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
230 INTEGER iMin, iMax, jMin, jMax
231 INTEGER k, bi, bj
232 _RL myTime
233 INTEGER myIter
234 INTEGER myThid
235
236 C !LOCAL VARIABLES:
237 C i,j :: Loop counters
238 C kSurface :: index of surface level
239 INTEGER i, j
240 #ifdef USE_OLD_EXTERNAL_FORCING
241 _RL locVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
242 _RL tmpVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
243 #else
244 INTEGER kSurface
245 #endif /* USE_OLD_EXTERNAL_FORCING */
246 CEOP
247
248 #ifdef USE_OLD_EXTERNAL_FORCING
249
250 DO j=1-OLy,sNy+OLy
251 DO i=1-OLx,sNx+OLx
252 locVar(i,j) = gV(i,j,k,bi,bj)
253 ENDDO
254 ENDDO
255 CALL EXTERNAL_FORCING_V(
256 I iMin, iMax, jMin, jMax, bi, bj, k,
257 I myTime, myThid )
258 C- Use 2-d local array tmpVar and split loop in 2 parts to avoid compiler
259 C to mess-up this part by re-arranging the order of instructions (wrong
260 C when gV and gV_arr are the same array, i.e., called with argument gV).
261 DO j=1-OLy,sNy+OLy
262 DO i=1-OLx,sNx+OLx
263 tmpVar(i,j) = gV(i,j,k,bi,bj) - locVar(i,j)
264 gV(i,j,k,bi,bj) = locVar(i,j)
265 ENDDO
266 ENDDO
267 C- not needed since APPLY_FORCING_V is no longer called with argument gV
268 c CALL FOOL_THE_COMPILER_RL( tmpVar(1,1) )
269 DO j=1-OLy,sNy+OLy
270 DO i=1-OLx,sNx+OLx
271 gV_arr(i,j) = gV_arr(i,j) + tmpVar(i,j)
272 ENDDO
273 ENDDO
274
275 #else /* USE_OLD_EXTERNAL_FORCING */
276
277 IF ( fluidIsAir ) THEN
278 kSurface = 0
279 ELSEIF ( usingPCoords ) THEN
280 kSurface = Nr
281 ELSE
282 kSurface = 1
283 ENDIF
284
285 C-- Forcing term
286 #ifdef ALLOW_AIM
287 IF ( useAIM ) CALL AIM_TENDENCY_APPLY_V(
288 U gV_arr,
289 I iMin,iMax,jMin,jMax, k, bi,bj,
290 I myTime, 0, myThid )
291 #endif /* ALLOW_AIM */
292
293 #ifdef ALLOW_ATM_PHYS
294 IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_V(
295 U gV_arr,
296 I iMin,iMax,jMin,jMax, k, bi,bj,
297 I myTime, 0, myThid )
298 #endif /* ALLOW_ATM_PHYS */
299
300 #ifdef ALLOW_FIZHI
301 IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_V(
302 U gV_arr,
303 I iMin,iMax,jMin,jMax, k, bi,bj,
304 I myTime, 0, myThid )
305 #endif /* ALLOW_FIZHI */
306
307 C Ocean: Add momentum surface forcing (e.g., wind-stress) in surface level
308 IF ( k .EQ. kSurface ) THEN
309 DO j=1,sNy+1
310 c DO i=1,sNx
311 C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNx+1]
312 DO i=0,sNx+1
313 gV_arr(i,j) = gV_arr(i,j)
314 & +foFacMom*surfaceForcingV(i,j,bi,bj)
315 & *recip_drF(k)*_recip_hFacS(i,j,k,bi,bj)
316 ENDDO
317 ENDDO
318 ELSEIF ( kSurface.EQ.-1 ) THEN
319 DO j=1,sNy+1
320 DO i=0,sNx+1
321 IF ( kSurfS(i,j,bi,bj).EQ.k ) THEN
322 gV_arr(i,j) = gV_arr(i,j)
323 & +foFacMom*surfaceForcingV(i,j,bi,bj)
324 & *recip_drF(k)*_recip_hFacS(i,j,k,bi,bj)
325 ENDIF
326 ENDDO
327 ENDDO
328 ENDIF
329
330 #ifdef ALLOW_EDDYPSI
331 CALL TAUEDDY_TENDENCY_APPLY_V(
332 U gV_arr,
333 I iMin,iMax,jMin,jMax, k, bi,bj,
334 I myTime, 0, myThid )
335 #endif
336
337 #ifdef ALLOW_RBCS
338 IF (useRBCS) THEN
339 CALL RBCS_ADD_TENDENCY(
340 U gV_arr,
341 I k, bi, bj, -2,
342 I myTime, 0, myThid )
343 ENDIF
344 #endif /* ALLOW_RBCS */
345
346 #ifdef ALLOW_OBCS
347 IF (useOBCS) THEN
348 CALL OBCS_SPONGE_V(
349 U gV_arr,
350 I iMin,iMax,jMin,jMax, k, bi,bj,
351 I myTime, 0, myThid )
352 ENDIF
353 #endif /* ALLOW_OBCS */
354
355 #ifdef ALLOW_MYPACKAGE
356 IF ( useMYPACKAGE ) THEN
357 CALL MYPACKAGE_TENDENCY_APPLY_V(
358 U gV_arr,
359 I iMin,iMax,jMin,jMax, k, bi,bj,
360 I myTime, 0, myThid )
361 ENDIF
362 #endif /* ALLOW_MYPACKAGE */
363
364 #endif /* USE_OLD_EXTERNAL_FORCING */
365
366 RETURN
367 END
368
369 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
370 CBOP
371 C !ROUTINE: APPLY_FORCING_T
372 C !INTERFACE:
373 SUBROUTINE APPLY_FORCING_T(
374 U gT_arr,
375 I iMin,iMax,jMin,jMax, k, bi, bj,
376 I myTime, myIter, myThid )
377 C !DESCRIPTION: \bv
378 C *==========================================================*
379 C | S/R APPLY_FORCING_T
380 C | o Contains problem specific forcing for temperature.
381 C *==========================================================*
382 C | Adds terms to gT for forcing by external sources
383 C | e.g. heat flux, climatalogical relaxation, etc ...
384 C *==========================================================*
385 C \ev
386
387 C !USES:
388 IMPLICIT NONE
389 C == Global data ==
390 #include "SIZE.h"
391 #include "EEPARAMS.h"
392 #include "PARAMS.h"
393 #include "GRID.h"
394 #include "DYNVARS.h"
395 #include "FFIELDS.h"
396 #include "SURFACE.h"
397
398 C !INPUT/OUTPUT PARAMETERS:
399 C gT_arr :: the tendency array
400 C iMin,iMax :: Working range of x-index for applying forcing.
401 C jMin,jMax :: Working range of y-index for applying forcing.
402 C k :: Current vertical level index
403 C bi,bj :: Current tile indices
404 C myTime :: Current time in simulation
405 C myIter :: Current iteration number
406 C myThid :: my Thread Id number
407 _RL gT_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
408 INTEGER iMin, iMax, jMin, jMax
409 INTEGER k, bi, bj
410 _RL myTime
411 INTEGER myIter
412 INTEGER myThid
413
414 C !LOCAL VARIABLES:
415 C i,j :: Loop counters
416 C kSurface :: index of surface level
417 INTEGER i, j
418 #ifndef USE_OLD_EXTERNAL_FORCING
419 INTEGER kSurface
420 INTEGER km, kc, kp
421 _RL tmpVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
422 _RL tmpFac, delPI
423 _RL recip_Cp
424 #ifdef SHORTWAVE_HEATING
425 _RL minusone
426 PARAMETER (minusOne=-1.)
427 _RL swfracb(2)
428 INTEGER kp1
429 #endif
430 #endif /* USE_OLD_EXTERNAL_FORCING */
431 CEOP
432
433 #ifdef USE_OLD_EXTERNAL_FORCING
434
435 DO j=1-OLy,sNy+OLy
436 DO i=1-OLx,sNx+OLx
437 gT(i,j,k,bi,bj) = 0. _d 0
438 ENDDO
439 ENDDO
440 CALL EXTERNAL_FORCING_T(
441 I iMin, iMax, jMin, jMax, bi, bj, k,
442 I myTime, myThid )
443 DO j=1-OLy,sNy+OLy
444 DO i=1-OLx,sNx+OLx
445 gT_arr(i,j) = gT_arr(i,j) + gT(i,j,k,bi,bj)
446 ENDDO
447 ENDDO
448
449 #else /* USE_OLD_EXTERNAL_FORCING */
450
451 IF ( fluidIsAir ) THEN
452 kSurface = 0
453 ELSEIF ( usingZCoords .AND. useShelfIce ) THEN
454 kSurface = -1
455 ELSEIF ( usingPCoords ) THEN
456 kSurface = Nr
457 ELSE
458 kSurface = 1
459 ENDIF
460 recip_Cp = 1. _d 0 / HeatCapacity_Cp
461
462 C-- Note on loop range: For model dynamics, only needs to get correct
463 C forcing (update gT_arr) in tile interior (1:sNx,1:sNy);
464 C However, for some diagnostics, we may want to get valid forcing
465 C extended over 1 point in tile halo region (0:sNx+1,0:sNy=1).
466
467 C-- Forcing term
468 #ifdef ALLOW_AIM
469 IF ( useAIM ) CALL AIM_TENDENCY_APPLY_T(
470 U gT_arr,
471 I iMin,iMax,jMin,jMax, k, bi,bj,
472 I myTime, 0, myThid )
473 #endif /* ALLOW_AIM */
474
475 #ifdef ALLOW_ATM_PHYS
476 IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_T(
477 U gT_arr,
478 I iMin,iMax,jMin,jMax, k, bi,bj,
479 I myTime, 0, myThid )
480 #endif /* ALLOW_ATM_PHYS */
481
482 #ifdef ALLOW_FIZHI
483 IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_T(
484 U gT_arr,
485 I iMin,iMax,jMin,jMax, k, bi,bj,
486 I myTime, 0, myThid )
487 #endif /* ALLOW_FIZHI */
488
489 #ifdef ALLOW_ADDFLUID
490 IF ( selectAddFluid.NE.0 .AND. temp_addMass.NE.UNSET_RL ) THEN
491 IF ( ( selectAddFluid.GE.1 .AND. nonlinFreeSurf.GT.0 )
492 & .OR. convertFW2Salt.EQ.-1. _d 0 ) THEN
493 DO j=0,sNy+1
494 DO i=0,sNx+1
495 gT_arr(i,j) = gT_arr(i,j)
496 & + addMass(i,j,k,bi,bj)*mass2rUnit
497 & *( temp_addMass - theta(i,j,k,bi,bj) )
498 & *recip_rA(i,j,bi,bj)
499 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
500 C & *recip_deepFac2C(k)*recip_rhoFacC(k)
501 ENDDO
502 ENDDO
503 ELSE
504 DO j=0,sNy+1
505 DO i=0,sNx+1
506 gT_arr(i,j) = gT_arr(i,j)
507 & + addMass(i,j,k,bi,bj)*mass2rUnit
508 & *( temp_addMass - tRef(k) )
509 & *recip_rA(i,j,bi,bj)
510 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
511 C & *recip_deepFac2C(k)*recip_rhoFacC(k)
512 ENDDO
513 ENDDO
514 ENDIF
515 ENDIF
516 #endif /* ALLOW_ADDFLUID */
517
518 #ifdef ALLOW_FRICTION_HEATING
519 IF ( addFrictionHeating ) THEN
520 IF ( fluidIsAir ) THEN
521 C conversion from in-situ Temp to Pot.Temp
522 tmpFac = (atm_Po/rC(k))**atm_kappa
523 C conversion from W/m^2/r_unit to K/s
524 tmpFac = (tmpFac/atm_Cp) * mass2rUnit
525 ELSE
526 C conversion from W/m^2/r_unit to K/s
527 tmpFac = recip_Cp * mass2rUnit
528 ENDIF
529 DO j=0,sNy+1
530 DO i=0,sNx+1
531 gT_arr(i,j) = gT_arr(i,j)
532 & + frictionHeating(i,j,k,bi,bj)
533 & *tmpFac*recip_rA(i,j,bi,bj)
534 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
535 ENDDO
536 ENDDO
537 ENDIF
538 #endif /* ALLOW_FRICTION_HEATING */
539
540 IF ( fluidIsAir .AND. atm_Rq.NE.zeroRL .AND. Nr.NE.1 ) THEN
541 C-- Compressible fluid: account for difference between moist and dry air
542 C specific volume in Enthalpy equation (+ V.dP term), since only the
543 C dry air part is accounted for in the (dry) Pot.Temp formulation.
544 C Used centered averaging from interface to center (consistent with
545 C conversion term in KE eq) and same discretisation ( [T*Q]_bar_k )
546 C as for Theta_v in CALC_PHI_HYD
547
548 C conversion from in-situ Temp to Pot.Temp
549 tmpFac = (atm_Po/rC(k))**atm_kappa
550 C conversion from W/kg to K/s
551 tmpFac = tmpFac/atm_Cp
552 km = k-1
553 kc = k
554 kp = k+1
555 IF ( k.EQ.1 ) THEN
556 DO j=0,sNy+1
557 DO i=0,sNx+1
558 tmpVar(i,j) = 0.
559 ENDDO
560 ENDDO
561 ELSE
562 delPI = atm_Cp*( (rC(km)/atm_Po)**atm_kappa
563 & - (rC(kc)/atm_Po)**atm_kappa )
564 DO j=0,sNy+1
565 DO i=0,sNx+1
566 tmpVar(i,j) = wVel(i,j,kc,bi,bj)*delPI*atm_Rq
567 & *( theta(i,j,km,bi,bj)*salt(i,j,km,bi,bj)
568 & + theta(i,j,kc,bi,bj)*salt(i,j,kc,bi,bj)
569 & )*maskC(i,j,km,bi,bj)*0.25 _d 0
570 ENDDO
571 ENDDO
572 ENDIF
573 IF ( k.LT.Nr ) THEN
574 delPI = atm_Cp*( (rC(kc)/atm_Po)**atm_kappa
575 & - (rC(kp)/atm_Po)**atm_kappa )
576 DO j=0,sNy+1
577 DO i=0,sNx+1
578 tmpVar(i,j) = tmpVar(i,j)
579 & + wVel(i,j,kp,bi,bj)*delPI*atm_Rq
580 & *( theta(i,j,kc,bi,bj)*salt(i,j,kc,bi,bj)
581 & + theta(i,j,kp,bi,bj)*salt(i,j,kp,bi,bj)
582 & )*maskC(i,j,kp,bi,bj)*0.25 _d 0
583 ENDDO
584 ENDDO
585 ENDIF
586 DO j=0,sNy+1
587 DO i=0,sNx+1
588 gT_arr(i,j) = gT_arr(i,j)
589 & + tmpVar(i,j)*tmpFac
590 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
591 ENDDO
592 ENDDO
593 #ifdef ALLOW_DIAGNOSTICS
594 IF ( useDiagnostics ) THEN
595 C conversion to W/m^2
596 tmpFac = rUnit2mass
597 CALL DIAGNOSTICS_SCALE_FILL( tmpVar, tmpFac, 1,
598 & 'MoistCor', kc, 1, 2, bi,bj,myThid )
599 ENDIF
600 #endif /* ALLOW_DIAGNOSTICS */
601 ENDIF
602
603 C Ocean: Add temperature surface forcing (e.g., heat-flux) in surface level
604 IF ( k .EQ. kSurface ) THEN
605 DO j=0,sNy+1
606 DO i=0,sNx+1
607 gT_arr(i,j) = gT_arr(i,j)
608 & +surfaceForcingT(i,j,bi,bj)
609 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
610 ENDDO
611 ENDDO
612 ELSEIF ( kSurface.EQ.-1 ) THEN
613 DO j=0,sNy+1
614 DO i=0,sNx+1
615 IF ( kSurfC(i,j,bi,bj).EQ.k ) THEN
616 gT_arr(i,j) = gT_arr(i,j)
617 & +surfaceForcingT(i,j,bi,bj)
618 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
619 ENDIF
620 ENDDO
621 ENDDO
622 ENDIF
623
624 IF (linFSConserveTr) THEN
625 DO j=0,sNy+1
626 DO i=0,sNx+1
627 IF (k .EQ. kSurfC(i,j,bi,bj)) THEN
628 gT_arr(i,j) = gT_arr(i,j)
629 & +TsurfCor*recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
630 ENDIF
631 ENDDO
632 ENDDO
633 ENDIF
634
635 #ifdef ALLOW_GEOTHERMAL_FLUX
636 IF ( usingZCoords ) THEN
637 DO j=0,sNy+1
638 DO i=0,sNx+1
639 IF ( k.EQ.kLowC(i,j,bi,bj) ) THEN
640 gT_arr(i,j)=gT_arr(i,j)
641 & + geothermalFlux(i,j,bi,bj)
642 & *recip_Cp*mass2rUnit
643 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
644 ENDIF
645 ENDDO
646 ENDDO
647 ENDIF
648 #endif /* ALLOW_GEOTHERMAL_FLUX */
649
650 #ifdef SHORTWAVE_HEATING
651 C Penetrating SW radiation
652 c IF ( usePenetratingSW ) THEN
653 swfracb(1)=abs(rF(k))
654 swfracb(2)=abs(rF(k+1))
655 CALL SWFRAC(
656 I 2, minusOne,
657 U swfracb,
658 I myTime, 1, myThid )
659 kp1 = k+1
660 IF (k.EQ.Nr) THEN
661 kp1 = k
662 swfracb(2)=0. _d 0
663 ENDIF
664 DO j=0,sNy+1
665 DO i=0,sNx+1
666 gT_arr(i,j) = gT_arr(i,j)
667 & -Qsw(i,j,bi,bj)*(swfracb(1)*maskC(i,j,k,bi,bj)
668 & -swfracb(2)*maskC(i,j,kp1, bi,bj))
669 & *recip_Cp*mass2rUnit
670 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
671 ENDDO
672 ENDDO
673 c ENDIF
674 #endif
675
676 #ifdef ALLOW_FRAZIL
677 IF ( useFRAZIL )
678 & CALL FRAZIL_TENDENCY_APPLY_T(
679 U gT_arr,
680 I iMin,iMax,jMin,jMax, k, bi,bj,
681 I myTime, 0, myThid )
682 #endif /* ALLOW_FRAZIL */
683
684 #ifdef ALLOW_SHELFICE
685 IF ( useShelfIce )
686 & CALL SHELFICE_FORCING_T(
687 U gT_arr,
688 I iMin,iMax,jMin,jMax, k, bi,bj,
689 I myTime, 0, myThid )
690 #endif /* ALLOW_SHELFICE */
691
692 #ifdef ALLOW_ICEFRONT
693 IF ( useICEFRONT )
694 & CALL ICEFRONT_TENDENCY_APPLY_T(
695 U gT_arr,
696 I k, bi, bj, myTime, 0, myThid )
697 #endif /* ALLOW_ICEFRONT */
698
699 #ifdef ALLOW_SALT_PLUME
700 IF ( useSALT_PLUME )
701 & CALL SALT_PLUME_TENDENCY_APPLY_T(
702 U gT_arr,
703 I iMin,iMax,jMin,jMax, k, bi,bj,
704 I myTime, 0, myThid )
705 #endif /* ALLOW_SALT_PLUME */
706
707 #ifdef ALLOW_RBCS
708 IF (useRBCS) THEN
709 CALL RBCS_ADD_TENDENCY(
710 U gT_arr,
711 I k, bi, bj, 1,
712 I myTime, 0, myThid )
713 ENDIF
714 #endif /* ALLOW_RBCS */
715
716 #ifdef ALLOW_OBCS
717 IF (useOBCS) THEN
718 CALL OBCS_SPONGE_T(
719 U gT_arr,
720 I iMin,iMax,jMin,jMax, k, bi,bj,
721 I myTime, 0, myThid )
722 ENDIF
723 #endif /* ALLOW_OBCS */
724
725 #ifdef ALLOW_BBL
726 IF ( useBBL ) CALL BBL_TENDENCY_APPLY_T(
727 U gT_arr,
728 I iMin,iMax,jMin,jMax, k, bi,bj,
729 I myTime, 0, myThid )
730 #endif /* ALLOW_BBL */
731
732 #ifdef ALLOW_MYPACKAGE
733 IF ( useMYPACKAGE ) THEN
734 CALL MYPACKAGE_TENDENCY_APPLY_T(
735 U gT_arr,
736 I iMin,iMax,jMin,jMax, k, bi,bj,
737 I myTime, 0, myThid )
738 ENDIF
739 #endif /* ALLOW_MYPACKAGE */
740
741 #endif /* USE_OLD_EXTERNAL_FORCING */
742
743 RETURN
744 END
745
746 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
747 CBOP
748 C !ROUTINE: APPLY_FORCING_S
749 C !INTERFACE:
750 SUBROUTINE APPLY_FORCING_S(
751 U gS_arr,
752 I iMin,iMax,jMin,jMax, k, bi, bj,
753 I myTime, myIter, myThid )
754 C !DESCRIPTION: \bv
755 C *==========================================================*
756 C | S/R APPLY_FORCING_S
757 C | o Contains problem specific forcing for merid velocity.
758 C *==========================================================*
759 C | Adds terms to gS for forcing by external sources
760 C | e.g. fresh-water flux, climatalogical relaxation, etc ...
761 C *==========================================================*
762 C \ev
763
764 C !USES:
765 IMPLICIT NONE
766 C == Global data ==
767 #include "SIZE.h"
768 #include "EEPARAMS.h"
769 #include "PARAMS.h"
770 #include "GRID.h"
771 #include "DYNVARS.h"
772 #include "FFIELDS.h"
773 #include "SURFACE.h"
774
775 C !INPUT/OUTPUT PARAMETERS:
776 C gS_arr :: the tendency array
777 C iMin,iMax :: Working range of x-index for applying forcing.
778 C jMin,jMax :: Working range of y-index for applying forcing.
779 C k :: Current vertical level index
780 C bi,bj :: Current tile indices
781 C myTime :: Current time in simulation
782 C myIter :: Current iteration number
783 C myThid :: my Thread Id number
784 _RL gS_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
785 INTEGER iMin, iMax, jMin, jMax
786 INTEGER k, bi, bj
787 _RL myTime
788 INTEGER myIter
789 INTEGER myThid
790
791 C !LOCAL VARIABLES:
792 C i,j :: Loop counters
793 C kSurface :: index of surface level
794 INTEGER i, j
795 #ifndef USE_OLD_EXTERNAL_FORCING
796 INTEGER kSurface
797 #endif /* USE_OLD_EXTERNAL_FORCING */
798 CEOP
799
800 #ifdef USE_OLD_EXTERNAL_FORCING
801
802 DO j=1-OLy,sNy+OLy
803 DO i=1-OLx,sNx+OLx
804 gS(i,j,k,bi,bj) = 0. _d 0
805 ENDDO
806 ENDDO
807 CALL EXTERNAL_FORCING_S(
808 I iMin, iMax, jMin, jMax, bi, bj, k,
809 I myTime, myThid )
810 DO j=1-OLy,sNy+OLy
811 DO i=1-OLx,sNx+OLx
812 gS_arr(i,j) = gS_arr(i,j) + gS(i,j,k,bi,bj)
813 ENDDO
814 ENDDO
815
816 #else /* USE_OLD_EXTERNAL_FORCING */
817
818 IF ( fluidIsAir ) THEN
819 kSurface = 0
820 ELSEIF ( usingZCoords .AND. useShelfIce ) THEN
821 kSurface = -1
822 ELSEIF ( usingPCoords ) THEN
823 kSurface = Nr
824 ELSE
825 kSurface = 1
826 ENDIF
827
828 C-- Note on loop range: For model dynamics, only needs to get correct
829 C forcing (update gS_arr) in tile interior (1:sNx,1:sNy);
830 C However, for some diagnostics, we may want to get valid forcing
831 C extended over 1 point in tile halo region (0:sNx+1,0:sNy=1).
832
833 C-- Forcing term
834 #ifdef ALLOW_AIM
835 IF ( useAIM ) CALL AIM_TENDENCY_APPLY_S(
836 U gS_arr,
837 I iMin,iMax,jMin,jMax, k, bi,bj,
838 I myTime, 0, myThid )
839 #endif /* ALLOW_AIM */
840
841 #ifdef ALLOW_ATM_PHYS
842 IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_S(
843 U gS_arr,
844 I iMin,iMax,jMin,jMax, k, bi,bj,
845 I myTime, 0, myThid )
846 #endif /* ALLOW_ATM_PHYS */
847
848 #ifdef ALLOW_FIZHI
849 IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_S(
850 U gS_arr,
851 I iMin,iMax,jMin,jMax, k, bi,bj,
852 I myTime, 0, myThid )
853 #endif /* ALLOW_FIZHI */
854
855 #ifdef ALLOW_ADDFLUID
856 IF ( selectAddFluid.NE.0 .AND. salt_addMass.NE.UNSET_RL ) THEN
857 IF ( ( selectAddFluid.GE.1 .AND. nonlinFreeSurf.GT.0 )
858 & .OR. convertFW2Salt.EQ.-1. _d 0 ) THEN
859 DO j=0,sNy+1
860 DO i=0,sNx+1
861 gS_arr(i,j) = gS_arr(i,j)
862 & + addMass(i,j,k,bi,bj)*mass2rUnit
863 & *( salt_addMass - salt(i,j,k,bi,bj) )
864 & *recip_rA(i,j,bi,bj)
865 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
866 C & *recip_deepFac2C(k)*recip_rhoFacC(k)
867 ENDDO
868 ENDDO
869 ELSE
870 DO j=0,sNy+1
871 DO i=0,sNx+1
872 gS_arr(i,j) = gS_arr(i,j)
873 & + addMass(i,j,k,bi,bj)*mass2rUnit
874 & *( salt_addMass - sRef(k) )
875 & *recip_rA(i,j,bi,bj)
876 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
877 C & *recip_deepFac2C(k)*recip_rhoFacC(k)
878 ENDDO
879 ENDDO
880 ENDIF
881 ENDIF
882 #endif /* ALLOW_ADDFLUID */
883
884 C Ocean: Add salinity surface forcing (e.g., fresh-water) in surface level
885 IF ( k .EQ. kSurface ) THEN
886 DO j=0,sNy+1
887 DO i=0,sNx+1
888 gS_arr(i,j) = gS_arr(i,j)
889 & +surfaceForcingS(i,j,bi,bj)
890 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
891 ENDDO
892 ENDDO
893 ELSEIF ( kSurface.EQ.-1 ) THEN
894 DO j=0,sNy+1
895 DO i=0,sNx+1
896 IF ( kSurfC(i,j,bi,bj).EQ.k ) THEN
897 gS_arr(i,j) = gS_arr(i,j)
898 & +surfaceForcingS(i,j,bi,bj)
899 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
900 ENDIF
901 ENDDO
902 ENDDO
903 ENDIF
904
905 IF (linFSConserveTr) THEN
906 DO j=0,sNy+1
907 DO i=0,sNx+1
908 IF (k .EQ. kSurfC(i,j,bi,bj)) THEN
909 gS_arr(i,j) = gS_arr(i,j)
910 & +SsurfCor*recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
911 ENDIF
912 ENDDO
913 ENDDO
914 ENDIF
915
916 #ifdef ALLOW_SHELFICE
917 IF ( useShelfIce )
918 & CALL SHELFICE_FORCING_S(
919 U gS_arr,
920 I iMin,iMax,jMin,jMax, k, bi,bj,
921 I myTime, 0, myThid )
922 #endif /* ALLOW_SHELFICE */
923
924 #ifdef ALLOW_ICEFRONT
925 IF ( useICEFRONT )
926 & CALL ICEFRONT_TENDENCY_APPLY_S(
927 U gS_arr,
928 I k, bi, bj, myTime, 0, myThid )
929 #endif /* ALLOW_ICEFRONT */
930
931 #ifdef ALLOW_SALT_PLUME
932 IF ( useSALT_PLUME )
933 & CALL SALT_PLUME_TENDENCY_APPLY_S(
934 U gS_arr,
935 I iMin,iMax,jMin,jMax, k, bi,bj,
936 I myTime, 0, myThid )
937 #endif /* ALLOW_SALT_PLUME */
938
939 #ifdef ALLOW_RBCS
940 IF (useRBCS) THEN
941 CALL RBCS_ADD_TENDENCY(
942 U gS_arr,
943 I k, bi, bj, 2,
944 I myTime, 0, myThid )
945 ENDIF
946 #endif /* ALLOW_RBCS */
947
948 #ifdef ALLOW_OBCS
949 IF (useOBCS) THEN
950 CALL OBCS_SPONGE_S(
951 U gS_arr,
952 I iMin,iMax,jMin,jMax, k, bi,bj,
953 I myTime, 0, myThid )
954 ENDIF
955 #endif /* ALLOW_OBCS */
956
957 #ifdef ALLOW_BBL
958 IF ( useBBL ) CALL BBL_TENDENCY_APPLY_S(
959 U gS_arr,
960 I iMin,iMax,jMin,jMax, k, bi,bj,
961 I myTime, 0, myThid )
962 #endif /* ALLOW_BBL */
963
964 #ifdef ALLOW_MYPACKAGE
965 IF ( useMYPACKAGE ) THEN
966 CALL MYPACKAGE_TENDENCY_APPLY_S(
967 U gS_arr,
968 I iMin,iMax,jMin,jMax, k, bi,bj,
969 I myTime, 0, myThid )
970 ENDIF
971 #endif /* ALLOW_MYPACKAGE */
972
973 #endif /* USE_OLD_EXTERNAL_FORCING */
974
975 RETURN
976 END
  ViewVC Help
Powered by ViewVC 1.1.22