/[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.9 - (show annotations) (download)
Wed Oct 4 20:40:35 2017 UTC (6 years, 8 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, HEAD
Changes since 1.8: +25 -1 lines
- add new 2-d forcing field for time-dependent geopotential anomaly (e.g.,
  tidal forcing), in m^2/s^2 ;

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

  ViewVC Help
Powered by ViewVC 1.1.22