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

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

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


Revision 1.9 - (hide 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 jmc 1.9 C $Header: /u/gcmpack/MITgcm/model/src/apply_forcing.F,v 1.8 2016/02/14 00:08:49 jmc Exp $
2 jmc 1.1 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 jmc 1.2 _RL tmpVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65 jmc 1.1 #else
66     INTEGER kSurface
67 jmc 1.6 #endif /* USE_OLD_EXTERNAL_FORCING */
68 jmc 1.1 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 jmc 1.2 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 jmc 1.1 DO j=1-OLy,sNy+OLy
84     DO i=1-OLx,sNx+OLx
85 jmc 1.2 tmpVar(i,j) = gU(i,j,k,bi,bj) - locVar(i,j)
86 jmc 1.1 gU(i,j,k,bi,bj) = locVar(i,j)
87 jmc 1.2 ENDDO
88     ENDDO
89 jmc 1.4 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 jmc 1.2 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 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
113 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
120 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
127 jmc 1.1 #endif /* ALLOW_FIZHI */
128    
129 jmc 1.9 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 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
169 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
177 jmc 1.1
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 jmc 1.8 I myTime, myIter, myThid )
187 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
196 jmc 1.1 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 jmc 1.2 _RL tmpVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
255 jmc 1.1 #else
256     INTEGER kSurface
257 jmc 1.6 #endif /* USE_OLD_EXTERNAL_FORCING */
258 jmc 1.1 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 jmc 1.2 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 jmc 1.1 DO j=1-OLy,sNy+OLy
274     DO i=1-OLx,sNx+OLx
275 jmc 1.2 tmpVar(i,j) = gV(i,j,k,bi,bj) - locVar(i,j)
276 jmc 1.1 gV(i,j,k,bi,bj) = locVar(i,j)
277 jmc 1.2 ENDDO
278     ENDDO
279 jmc 1.4 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 jmc 1.2 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 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
303 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
310 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
317 jmc 1.1 #endif /* ALLOW_FIZHI */
318    
319 jmc 1.9 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 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
359 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
367 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
376 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
385 jmc 1.1 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 jmc 1.6 #ifndef USE_OLD_EXTERNAL_FORCING
443 jmc 1.1 INTEGER kSurface
444     INTEGER km, kc, kp
445 jmc 1.6 _RL tmpVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
446 jmc 1.1 _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 jmc 1.6 #endif /* USE_OLD_EXTERNAL_FORCING */
455     CEOP
456 jmc 1.1
457     #ifdef USE_OLD_EXTERNAL_FORCING
458    
459     DO j=1-OLy,sNy+OLy
460     DO i=1-OLx,sNx+OLx
461 jmc 1.6 gT(i,j,k,bi,bj) = 0. _d 0
462 jmc 1.1 ENDDO
463     ENDDO
464     CALL EXTERNAL_FORCING_T(
465     I iMin, iMax, jMin, jMax, bi, bj, k,
466     I myTime, myThid )
467 jmc 1.2 DO j=1-OLy,sNy+OLy
468     DO i=1-OLx,sNx+OLx
469 jmc 1.6 gT_arr(i,j) = gT_arr(i,j) + gT(i,j,k,bi,bj)
470 jmc 1.1 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 jmc 1.6 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 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
497 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
504 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
511 jmc 1.1 #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 jmc 1.6 DO j=0,sNy+1
518     DO i=0,sNx+1
519 jmc 1.1 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 jmc 1.6 DO j=0,sNy+1
529     DO i=0,sNx+1
530 jmc 1.1 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 jmc 1.6 DO j=0,sNy+1
554     DO i=0,sNx+1
555 jmc 1.1 gT_arr(i,j) = gT_arr(i,j)
556 jmc 1.7 & + frictionHeating(i,j,k,bi,bj)*tmpFac
557 jmc 1.1 & *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 jmc 1.6 DO j=0,sNy+1
580     DO i=0,sNx+1
581 jmc 1.1 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 jmc 1.6 DO j=0,sNy+1
588     DO i=0,sNx+1
589 jmc 1.1 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 jmc 1.6 DO j=0,sNy+1
600     DO i=0,sNx+1
601 jmc 1.1 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 jmc 1.6 DO j=0,sNy+1
610     DO i=0,sNx+1
611 jmc 1.1 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 jmc 1.6 & 'MoistCor', kc, 1, 2, bi,bj,myThid )
622 jmc 1.1 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 jmc 1.6 DO j=0,sNy+1
629     DO i=0,sNx+1
630 jmc 1.1 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 jmc 1.6 DO j=0,sNy+1
637     DO i=0,sNx+1
638 jmc 1.1 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 jmc 1.6 DO j=0,sNy+1
649     DO i=0,sNx+1
650 jmc 1.1 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 heimbach 1.5 #ifdef ALLOW_GEOTHERMAL_FLUX
659     IF ( usingZCoords ) THEN
660 jmc 1.6 DO j=0,sNy+1
661     DO i=0,sNx+1
662 heimbach 1.5 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 jmc 1.1 #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 jmc 1.6 DO j=0,sNy+1
688     DO i=0,sNx+1
689 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
705 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
713 jmc 1.1 #endif /* ALLOW_SHELFICE */
714    
715     #ifdef ALLOW_ICEFRONT
716     IF ( useICEFRONT )
717     & CALL ICEFRONT_TENDENCY_APPLY_T(
718     U gT_arr,
719 jmc 1.8 I k, bi, bj, myTime, myIter, myThid )
720 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
728 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
736 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
745 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
753 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
761 jmc 1.1 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 jmc 1.6 #ifndef USE_OLD_EXTERNAL_FORCING
819 jmc 1.1 INTEGER kSurface
820 jmc 1.6 #endif /* USE_OLD_EXTERNAL_FORCING */
821 jmc 1.1 CEOP
822    
823     #ifdef USE_OLD_EXTERNAL_FORCING
824    
825     DO j=1-OLy,sNy+OLy
826     DO i=1-OLx,sNx+OLx
827 jmc 1.6 gS(i,j,k,bi,bj) = 0. _d 0
828 jmc 1.1 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 jmc 1.6 gS_arr(i,j) = gS_arr(i,j) + gS(i,j,k,bi,bj)
836 jmc 1.1 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 jmc 1.6 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 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
862 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
869 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
876 jmc 1.1 #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 jmc 1.6 DO j=0,sNy+1
883     DO i=0,sNx+1
884 jmc 1.1 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 jmc 1.6 DO j=0,sNy+1
894     DO i=0,sNx+1
895 jmc 1.1 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 jmc 1.6 DO j=0,sNy+1
910     DO i=0,sNx+1
911 jmc 1.1 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 jmc 1.6 DO j=0,sNy+1
918     DO i=0,sNx+1
919 jmc 1.1 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 jmc 1.6 DO j=0,sNy+1
930     DO i=0,sNx+1
931 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
945 jmc 1.1 #endif /* ALLOW_SHELFICE */
946    
947     #ifdef ALLOW_ICEFRONT
948     IF ( useICEFRONT )
949     & CALL ICEFRONT_TENDENCY_APPLY_S(
950     U gS_arr,
951 jmc 1.8 I k, bi, bj, myTime, myIter, myThid )
952 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
960 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
968 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
977 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
985 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
993 jmc 1.1 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