/[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.8 - (hide annotations) (download)
Sun Feb 14 00:08:49 2016 UTC (8 years, 3 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65v, checkpoint65w, checkpoint65u
Changes since 1.7: +36 -36 lines
pass myIter (instead of just 0) to various specific forcing S/R

1 jmc 1.8 C $Header: /u/gcmpack/MITgcm/model/src/apply_forcing.F,v 1.7 2015/01/21 14:35:16 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     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 jmc 1.8 I myTime, myIter, myThid )
157 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
165 jmc 1.1
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 jmc 1.8 I myTime, myIter, myThid )
175 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
184 jmc 1.1 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 jmc 1.2 _RL tmpVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
243 jmc 1.1 #else
244     INTEGER kSurface
245 jmc 1.6 #endif /* USE_OLD_EXTERNAL_FORCING */
246 jmc 1.1 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 jmc 1.2 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 jmc 1.1 DO j=1-OLy,sNy+OLy
262     DO i=1-OLx,sNx+OLx
263 jmc 1.2 tmpVar(i,j) = gV(i,j,k,bi,bj) - locVar(i,j)
264 jmc 1.1 gV(i,j,k,bi,bj) = locVar(i,j)
265 jmc 1.2 ENDDO
266     ENDDO
267 jmc 1.4 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 jmc 1.2 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 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
291 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
298 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
305 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
335 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
343 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
352 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
361 jmc 1.1 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 jmc 1.6 #ifndef USE_OLD_EXTERNAL_FORCING
419 jmc 1.1 INTEGER kSurface
420     INTEGER km, kc, kp
421 jmc 1.6 _RL tmpVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
422 jmc 1.1 _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 jmc 1.6 #endif /* USE_OLD_EXTERNAL_FORCING */
431     CEOP
432 jmc 1.1
433     #ifdef USE_OLD_EXTERNAL_FORCING
434    
435     DO j=1-OLy,sNy+OLy
436     DO i=1-OLx,sNx+OLx
437 jmc 1.6 gT(i,j,k,bi,bj) = 0. _d 0
438 jmc 1.1 ENDDO
439     ENDDO
440     CALL EXTERNAL_FORCING_T(
441     I iMin, iMax, jMin, jMax, bi, bj, k,
442     I myTime, myThid )
443 jmc 1.2 DO j=1-OLy,sNy+OLy
444     DO i=1-OLx,sNx+OLx
445 jmc 1.6 gT_arr(i,j) = gT_arr(i,j) + gT(i,j,k,bi,bj)
446 jmc 1.1 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 jmc 1.6 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 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
473 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
480 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
487 jmc 1.1 #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 jmc 1.6 DO j=0,sNy+1
494     DO i=0,sNx+1
495 jmc 1.1 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 jmc 1.6 DO j=0,sNy+1
505     DO i=0,sNx+1
506 jmc 1.1 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 jmc 1.6 DO j=0,sNy+1
530     DO i=0,sNx+1
531 jmc 1.1 gT_arr(i,j) = gT_arr(i,j)
532 jmc 1.7 & + frictionHeating(i,j,k,bi,bj)*tmpFac
533 jmc 1.1 & *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 jmc 1.6 DO j=0,sNy+1
556     DO i=0,sNx+1
557 jmc 1.1 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 jmc 1.6 DO j=0,sNy+1
564     DO i=0,sNx+1
565 jmc 1.1 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 jmc 1.6 DO j=0,sNy+1
576     DO i=0,sNx+1
577 jmc 1.1 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 jmc 1.6 DO j=0,sNy+1
586     DO i=0,sNx+1
587 jmc 1.1 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 jmc 1.6 & 'MoistCor', kc, 1, 2, bi,bj,myThid )
598 jmc 1.1 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 jmc 1.6 DO j=0,sNy+1
605     DO i=0,sNx+1
606 jmc 1.1 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 jmc 1.6 DO j=0,sNy+1
613     DO i=0,sNx+1
614 jmc 1.1 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 jmc 1.6 DO j=0,sNy+1
625     DO i=0,sNx+1
626 jmc 1.1 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 heimbach 1.5 #ifdef ALLOW_GEOTHERMAL_FLUX
635     IF ( usingZCoords ) THEN
636 jmc 1.6 DO j=0,sNy+1
637     DO i=0,sNx+1
638 heimbach 1.5 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 jmc 1.1 #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 jmc 1.6 DO j=0,sNy+1
664     DO i=0,sNx+1
665 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
681 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
689 jmc 1.1 #endif /* ALLOW_SHELFICE */
690    
691     #ifdef ALLOW_ICEFRONT
692     IF ( useICEFRONT )
693     & CALL ICEFRONT_TENDENCY_APPLY_T(
694     U gT_arr,
695 jmc 1.8 I k, bi, bj, myTime, myIter, myThid )
696 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
704 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
712 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
721 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
729 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
737 jmc 1.1 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 jmc 1.6 #ifndef USE_OLD_EXTERNAL_FORCING
795 jmc 1.1 INTEGER kSurface
796 jmc 1.6 #endif /* USE_OLD_EXTERNAL_FORCING */
797 jmc 1.1 CEOP
798    
799     #ifdef USE_OLD_EXTERNAL_FORCING
800    
801     DO j=1-OLy,sNy+OLy
802     DO i=1-OLx,sNx+OLx
803 jmc 1.6 gS(i,j,k,bi,bj) = 0. _d 0
804 jmc 1.1 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 jmc 1.6 gS_arr(i,j) = gS_arr(i,j) + gS(i,j,k,bi,bj)
812 jmc 1.1 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 jmc 1.6 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 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
838 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
845 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
852 jmc 1.1 #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 jmc 1.6 DO j=0,sNy+1
859     DO i=0,sNx+1
860 jmc 1.1 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 jmc 1.6 DO j=0,sNy+1
870     DO i=0,sNx+1
871 jmc 1.1 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 jmc 1.6 DO j=0,sNy+1
886     DO i=0,sNx+1
887 jmc 1.1 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 jmc 1.6 DO j=0,sNy+1
894     DO i=0,sNx+1
895 jmc 1.1 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 jmc 1.6 DO j=0,sNy+1
906     DO i=0,sNx+1
907 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
921 jmc 1.1 #endif /* ALLOW_SHELFICE */
922    
923     #ifdef ALLOW_ICEFRONT
924     IF ( useICEFRONT )
925     & CALL ICEFRONT_TENDENCY_APPLY_S(
926     U gS_arr,
927 jmc 1.8 I k, bi, bj, myTime, myIter, myThid )
928 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
936 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
944 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
953 jmc 1.1 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 jmc 1.8 I myTime, myIter, myThid )
961 jmc 1.1 #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 jmc 1.8 I myTime, myIter, myThid )
969 jmc 1.1 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