/[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.7 - (show annotations) (download)
Wed Jan 21 14:35:16 2015 UTC (9 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65t, checkpoint65j, checkpoint65k, checkpoint65i, checkpoint65n, checkpoint65o, checkpoint65l, checkpoint65m
Changes since 1.6: +2 -3 lines 
change units of frictionHeating field (from W to W/m^2)
1 C $Header: /u/gcmpack/MITgcm/model/src/apply_forcing.F,v 1.6 2014/09/05 21:04:34 jmc 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)*tmpFac
533 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
534 ENDDO
535 ENDDO
536 ENDIF
537 #endif /* ALLOW_FRICTION_HEATING */
538
539 IF ( fluidIsAir .AND. atm_Rq.NE.zeroRL .AND. Nr.NE.1 ) THEN
540 C-- Compressible fluid: account for difference between moist and dry air
541 C specific volume in Enthalpy equation (+ V.dP term), since only the
542 C dry air part is accounted for in the (dry) Pot.Temp formulation.
543 C Used centered averaging from interface to center (consistent with
544 C conversion term in KE eq) and same discretisation ( [T*Q]_bar_k )
545 C as for Theta_v in CALC_PHI_HYD
546
547 C conversion from in-situ Temp to Pot.Temp
548 tmpFac = (atm_Po/rC(k))**atm_kappa
549 C conversion from W/kg to K/s
550 tmpFac = tmpFac/atm_Cp
551 km = k-1
552 kc = k
553 kp = k+1
554 IF ( k.EQ.1 ) THEN
555 DO j=0,sNy+1
556 DO i=0,sNx+1
557 tmpVar(i,j) = 0.
558 ENDDO
559 ENDDO
560 ELSE
561 delPI = atm_Cp*( (rC(km)/atm_Po)**atm_kappa
562 & - (rC(kc)/atm_Po)**atm_kappa )
563 DO j=0,sNy+1
564 DO i=0,sNx+1
565 tmpVar(i,j) = wVel(i,j,kc,bi,bj)*delPI*atm_Rq
566 & *( theta(i,j,km,bi,bj)*salt(i,j,km,bi,bj)
567 & + theta(i,j,kc,bi,bj)*salt(i,j,kc,bi,bj)
568 & )*maskC(i,j,km,bi,bj)*0.25 _d 0
569 ENDDO
570 ENDDO
571 ENDIF
572 IF ( k.LT.Nr ) THEN
573 delPI = atm_Cp*( (rC(kc)/atm_Po)**atm_kappa
574 & - (rC(kp)/atm_Po)**atm_kappa )
575 DO j=0,sNy+1
576 DO i=0,sNx+1
577 tmpVar(i,j) = tmpVar(i,j)
578 & + wVel(i,j,kp,bi,bj)*delPI*atm_Rq
579 & *( theta(i,j,kc,bi,bj)*salt(i,j,kc,bi,bj)
580 & + theta(i,j,kp,bi,bj)*salt(i,j,kp,bi,bj)
581 & )*maskC(i,j,kp,bi,bj)*0.25 _d 0
582 ENDDO
583 ENDDO
584 ENDIF
585 DO j=0,sNy+1
586 DO i=0,sNx+1
587 gT_arr(i,j) = gT_arr(i,j)
588 & + tmpVar(i,j)*tmpFac
589 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
590 ENDDO
591 ENDDO
592 #ifdef ALLOW_DIAGNOSTICS
593 IF ( useDiagnostics ) THEN
594 C conversion to W/m^2
595 tmpFac = rUnit2mass
596 CALL DIAGNOSTICS_SCALE_FILL( tmpVar, tmpFac, 1,
597 & 'MoistCor', kc, 1, 2, bi,bj,myThid )
598 ENDIF
599 #endif /* ALLOW_DIAGNOSTICS */
600 ENDIF
601
602 C Ocean: Add temperature surface forcing (e.g., heat-flux) in surface level
603 IF ( k .EQ. kSurface ) THEN
604 DO j=0,sNy+1
605 DO i=0,sNx+1
606 gT_arr(i,j) = gT_arr(i,j)
607 & +surfaceForcingT(i,j,bi,bj)
608 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
609 ENDDO
610 ENDDO
611 ELSEIF ( kSurface.EQ.-1 ) THEN
612 DO j=0,sNy+1
613 DO i=0,sNx+1
614 IF ( kSurfC(i,j,bi,bj).EQ.k ) THEN
615 gT_arr(i,j) = gT_arr(i,j)
616 & +surfaceForcingT(i,j,bi,bj)
617 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
618 ENDIF
619 ENDDO
620 ENDDO
621 ENDIF
622
623 IF (linFSConserveTr) THEN
624 DO j=0,sNy+1
625 DO i=0,sNx+1
626 IF (k .EQ. kSurfC(i,j,bi,bj)) THEN
627 gT_arr(i,j) = gT_arr(i,j)
628 & +TsurfCor*recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
629 ENDIF
630 ENDDO
631 ENDDO
632 ENDIF
633
634 #ifdef ALLOW_GEOTHERMAL_FLUX
635 IF ( usingZCoords ) THEN
636 DO j=0,sNy+1
637 DO i=0,sNx+1
638 IF ( k.EQ.kLowC(i,j,bi,bj) ) THEN
639 gT_arr(i,j)=gT_arr(i,j)
640 & + geothermalFlux(i,j,bi,bj)
641 & *recip_Cp*mass2rUnit
642 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
643 ENDIF
644 ENDDO
645 ENDDO
646 ENDIF
647 #endif /* ALLOW_GEOTHERMAL_FLUX */
648
649 #ifdef SHORTWAVE_HEATING
650 C Penetrating SW radiation
651 c IF ( usePenetratingSW ) THEN
652 swfracb(1)=abs(rF(k))
653 swfracb(2)=abs(rF(k+1))
654 CALL SWFRAC(
655 I 2, minusOne,
656 U swfracb,
657 I myTime, 1, myThid )
658 kp1 = k+1
659 IF (k.EQ.Nr) THEN
660 kp1 = k
661 swfracb(2)=0. _d 0
662 ENDIF
663 DO j=0,sNy+1
664 DO i=0,sNx+1
665 gT_arr(i,j) = gT_arr(i,j)
666 & -Qsw(i,j,bi,bj)*(swfracb(1)*maskC(i,j,k,bi,bj)
667 & -swfracb(2)*maskC(i,j,kp1, bi,bj))
668 & *recip_Cp*mass2rUnit
669 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
670 ENDDO
671 ENDDO
672 c ENDIF
673 #endif
674
675 #ifdef ALLOW_FRAZIL
676 IF ( useFRAZIL )
677 & CALL FRAZIL_TENDENCY_APPLY_T(
678 U gT_arr,
679 I iMin,iMax,jMin,jMax, k, bi,bj,
680 I myTime, 0, myThid )
681 #endif /* ALLOW_FRAZIL */
682
683 #ifdef ALLOW_SHELFICE
684 IF ( useShelfIce )
685 & CALL SHELFICE_FORCING_T(
686 U gT_arr,
687 I iMin,iMax,jMin,jMax, k, bi,bj,
688 I myTime, 0, myThid )
689 #endif /* ALLOW_SHELFICE */
690
691 #ifdef ALLOW_ICEFRONT
692 IF ( useICEFRONT )
693 & CALL ICEFRONT_TENDENCY_APPLY_T(
694 U gT_arr,
695 I k, bi, bj, myTime, 0, myThid )
696 #endif /* ALLOW_ICEFRONT */
697
698 #ifdef ALLOW_SALT_PLUME
699 IF ( useSALT_PLUME )
700 & CALL SALT_PLUME_TENDENCY_APPLY_T(
701 U gT_arr,
702 I iMin,iMax,jMin,jMax, k, bi,bj,
703 I myTime, 0, myThid )
704 #endif /* ALLOW_SALT_PLUME */
705
706 #ifdef ALLOW_RBCS
707 IF (useRBCS) THEN
708 CALL RBCS_ADD_TENDENCY(
709 U gT_arr,
710 I k, bi, bj, 1,
711 I myTime, 0, myThid )
712 ENDIF
713 #endif /* ALLOW_RBCS */
714
715 #ifdef ALLOW_OBCS
716 IF (useOBCS) THEN
717 CALL OBCS_SPONGE_T(
718 U gT_arr,
719 I iMin,iMax,jMin,jMax, k, bi,bj,
720 I myTime, 0, myThid )
721 ENDIF
722 #endif /* ALLOW_OBCS */
723
724 #ifdef ALLOW_BBL
725 IF ( useBBL ) CALL BBL_TENDENCY_APPLY_T(
726 U gT_arr,
727 I iMin,iMax,jMin,jMax, k, bi,bj,
728 I myTime, 0, myThid )
729 #endif /* ALLOW_BBL */
730
731 #ifdef ALLOW_MYPACKAGE
732 IF ( useMYPACKAGE ) THEN
733 CALL MYPACKAGE_TENDENCY_APPLY_T(
734 U gT_arr,
735 I iMin,iMax,jMin,jMax, k, bi,bj,
736 I myTime, 0, myThid )
737 ENDIF
738 #endif /* ALLOW_MYPACKAGE */
739
740 #endif /* USE_OLD_EXTERNAL_FORCING */
741
742 RETURN
743 END
744
745 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
746 CBOP
747 C !ROUTINE: APPLY_FORCING_S
748 C !INTERFACE:
749 SUBROUTINE APPLY_FORCING_S(
750 U gS_arr,
751 I iMin,iMax,jMin,jMax, k, bi, bj,
752 I myTime, myIter, myThid )
753 C !DESCRIPTION: \bv
754 C *==========================================================*
755 C | S/R APPLY_FORCING_S
756 C | o Contains problem specific forcing for merid velocity.
757 C *==========================================================*
758 C | Adds terms to gS for forcing by external sources
759 C | e.g. fresh-water flux, climatalogical relaxation, etc ...
760 C *==========================================================*
761 C \ev
762
763 C !USES:
764 IMPLICIT NONE
765 C == Global data ==
766 #include "SIZE.h"
767 #include "EEPARAMS.h"
768 #include "PARAMS.h"
769 #include "GRID.h"
770 #include "DYNVARS.h"
771 #include "FFIELDS.h"
772 #include "SURFACE.h"
773
774 C !INPUT/OUTPUT PARAMETERS:
775 C gS_arr :: the tendency array
776 C iMin,iMax :: Working range of x-index for applying forcing.
777 C jMin,jMax :: Working range of y-index for applying forcing.
778 C k :: Current vertical level index
779 C bi,bj :: Current tile indices
780 C myTime :: Current time in simulation
781 C myIter :: Current iteration number
782 C myThid :: my Thread Id number
783 _RL gS_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
784 INTEGER iMin, iMax, jMin, jMax
785 INTEGER k, bi, bj
786 _RL myTime
787 INTEGER myIter
788 INTEGER myThid
789
790 C !LOCAL VARIABLES:
791 C i,j :: Loop counters
792 C kSurface :: index of surface level
793 INTEGER i, j
794 #ifndef USE_OLD_EXTERNAL_FORCING
795 INTEGER kSurface
796 #endif /* USE_OLD_EXTERNAL_FORCING */
797 CEOP
798
799 #ifdef USE_OLD_EXTERNAL_FORCING
800
801 DO j=1-OLy,sNy+OLy
802 DO i=1-OLx,sNx+OLx
803 gS(i,j,k,bi,bj) = 0. _d 0
804 ENDDO
805 ENDDO
806 CALL EXTERNAL_FORCING_S(
807 I iMin, iMax, jMin, jMax, bi, bj, k,
808 I myTime, myThid )
809 DO j=1-OLy,sNy+OLy
810 DO i=1-OLx,sNx+OLx
811 gS_arr(i,j) = gS_arr(i,j) + gS(i,j,k,bi,bj)
812 ENDDO
813 ENDDO
814
815 #else /* USE_OLD_EXTERNAL_FORCING */
816
817 IF ( fluidIsAir ) THEN
818 kSurface = 0
819 ELSEIF ( usingZCoords .AND. useShelfIce ) THEN
820 kSurface = -1
821 ELSEIF ( usingPCoords ) THEN
822 kSurface = Nr
823 ELSE
824 kSurface = 1
825 ENDIF
826
827 C-- Note on loop range: For model dynamics, only needs to get correct
828 C forcing (update gS_arr) in tile interior (1:sNx,1:sNy);
829 C However, for some diagnostics, we may want to get valid forcing
830 C extended over 1 point in tile halo region (0:sNx+1,0:sNy=1).
831
832 C-- Forcing term
833 #ifdef ALLOW_AIM
834 IF ( useAIM ) CALL AIM_TENDENCY_APPLY_S(
835 U gS_arr,
836 I iMin,iMax,jMin,jMax, k, bi,bj,
837 I myTime, 0, myThid )
838 #endif /* ALLOW_AIM */
839
840 #ifdef ALLOW_ATM_PHYS
841 IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_S(
842 U gS_arr,
843 I iMin,iMax,jMin,jMax, k, bi,bj,
844 I myTime, 0, myThid )
845 #endif /* ALLOW_ATM_PHYS */
846
847 #ifdef ALLOW_FIZHI
848 IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_S(
849 U gS_arr,
850 I iMin,iMax,jMin,jMax, k, bi,bj,
851 I myTime, 0, myThid )
852 #endif /* ALLOW_FIZHI */
853
854 #ifdef ALLOW_ADDFLUID
855 IF ( selectAddFluid.NE.0 .AND. salt_addMass.NE.UNSET_RL ) THEN
856 IF ( ( selectAddFluid.GE.1 .AND. nonlinFreeSurf.GT.0 )
857 & .OR. convertFW2Salt.EQ.-1. _d 0 ) THEN
858 DO j=0,sNy+1
859 DO i=0,sNx+1
860 gS_arr(i,j) = gS_arr(i,j)
861 & + addMass(i,j,k,bi,bj)*mass2rUnit
862 & *( salt_addMass - salt(i,j,k,bi,bj) )
863 & *recip_rA(i,j,bi,bj)
864 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
865 C & *recip_deepFac2C(k)*recip_rhoFacC(k)
866 ENDDO
867 ENDDO
868 ELSE
869 DO j=0,sNy+1
870 DO i=0,sNx+1
871 gS_arr(i,j) = gS_arr(i,j)
872 & + addMass(i,j,k,bi,bj)*mass2rUnit
873 & *( salt_addMass - sRef(k) )
874 & *recip_rA(i,j,bi,bj)
875 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
876 C & *recip_deepFac2C(k)*recip_rhoFacC(k)
877 ENDDO
878 ENDDO
879 ENDIF
880 ENDIF
881 #endif /* ALLOW_ADDFLUID */
882
883 C Ocean: Add salinity surface forcing (e.g., fresh-water) in surface level
884 IF ( k .EQ. kSurface ) THEN
885 DO j=0,sNy+1
886 DO i=0,sNx+1
887 gS_arr(i,j) = gS_arr(i,j)
888 & +surfaceForcingS(i,j,bi,bj)
889 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
890 ENDDO
891 ENDDO
892 ELSEIF ( kSurface.EQ.-1 ) THEN
893 DO j=0,sNy+1
894 DO i=0,sNx+1
895 IF ( kSurfC(i,j,bi,bj).EQ.k ) THEN
896 gS_arr(i,j) = gS_arr(i,j)
897 & +surfaceForcingS(i,j,bi,bj)
898 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
899 ENDIF
900 ENDDO
901 ENDDO
902 ENDIF
903
904 IF (linFSConserveTr) THEN
905 DO j=0,sNy+1
906 DO i=0,sNx+1
907 IF (k .EQ. kSurfC(i,j,bi,bj)) THEN
908 gS_arr(i,j) = gS_arr(i,j)
909 & +SsurfCor*recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
910 ENDIF
911 ENDDO
912 ENDDO
913 ENDIF
914
915 #ifdef ALLOW_SHELFICE
916 IF ( useShelfIce )
917 & CALL SHELFICE_FORCING_S(
918 U gS_arr,
919 I iMin,iMax,jMin,jMax, k, bi,bj,
920 I myTime, 0, myThid )
921 #endif /* ALLOW_SHELFICE */
922
923 #ifdef ALLOW_ICEFRONT
924 IF ( useICEFRONT )
925 & CALL ICEFRONT_TENDENCY_APPLY_S(
926 U gS_arr,
927 I k, bi, bj, myTime, 0, myThid )
928 #endif /* ALLOW_ICEFRONT */
929
930 #ifdef ALLOW_SALT_PLUME
931 IF ( useSALT_PLUME )
932 & CALL SALT_PLUME_TENDENCY_APPLY_S(
933 U gS_arr,
934 I iMin,iMax,jMin,jMax, k, bi,bj,
935 I myTime, 0, myThid )
936 #endif /* ALLOW_SALT_PLUME */
937
938 #ifdef ALLOW_RBCS
939 IF (useRBCS) THEN
940 CALL RBCS_ADD_TENDENCY(
941 U gS_arr,
942 I k, bi, bj, 2,
943 I myTime, 0, myThid )
944 ENDIF
945 #endif /* ALLOW_RBCS */
946
947 #ifdef ALLOW_OBCS
948 IF (useOBCS) THEN
949 CALL OBCS_SPONGE_S(
950 U gS_arr,
951 I iMin,iMax,jMin,jMax, k, bi,bj,
952 I myTime, 0, myThid )
953 ENDIF
954 #endif /* ALLOW_OBCS */
955
956 #ifdef ALLOW_BBL
957 IF ( useBBL ) CALL BBL_TENDENCY_APPLY_S(
958 U gS_arr,
959 I iMin,iMax,jMin,jMax, k, bi,bj,
960 I myTime, 0, myThid )
961 #endif /* ALLOW_BBL */
962
963 #ifdef ALLOW_MYPACKAGE
964 IF ( useMYPACKAGE ) THEN
965 CALL MYPACKAGE_TENDENCY_APPLY_S(
966 U gS_arr,
967 I iMin,iMax,jMin,jMax, k, bi,bj,
968 I myTime, 0, myThid )
969 ENDIF
970 #endif /* ALLOW_MYPACKAGE */
971
972 #endif /* USE_OLD_EXTERNAL_FORCING */
973
974 RETURN
975 END
  ViewVC Help
Powered by ViewVC 1.1.22