/[MITgcm]/MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/external_forcing.F
ViewVC logotype

Annotation of /MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/external_forcing.F

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

Revision 1.2 - (hide annotations) (download)
Tue Apr 22 04:32:26 2014 UTC (11 years, 3 months ago) by atn
Branch: MAIN
Changes since 1.1: +10 -7 lines 
in progress:
1. sort out bi,bj loop
2. skip remove saltplumeflux in external_forcing_surf, (thus) skip kpp
3. move saltplume outside k-loop in [salt,temp]_integrate.F
4. add diagnostics
1 atn 1.2 C $Header: /u/gcmpack/MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/external_forcing.F,v 1.1 2014/04/20 04:03:07 atn Exp $
2 atn 1.1 C $Name: $
3    
4     #include "PACKAGES_CONFIG.h"
5     #include "CPP_OPTIONS.h"
6    
7     C-- File external_forcing.F:
8     C-- Contents
9     C-- o EXTERNAL_FORCING_U
10     C-- o EXTERNAL_FORCING_V
11     C-- o EXTERNAL_FORCING_T
12     C-- o EXTERNAL_FORCING_S
13    
14     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
15     CBOP
16     C !ROUTINE: EXTERNAL_FORCING_U
17     C !INTERFACE:
18     SUBROUTINE EXTERNAL_FORCING_U(
19     I iMin,iMax, jMin,jMax, bi,bj, kLev,
20     I myTime, myThid )
21     C !DESCRIPTION: \bv
22     C *==========================================================*
23     C | S/R EXTERNAL_FORCING_U
24     C | o Contains problem specific forcing for zonal velocity.
25     C *==========================================================*
26     C | Adds terms to gU for forcing by external sources
27     C | e.g. wind stress, bottom friction etc ...
28     C *==========================================================*
29     C \ev
30    
31     C !USES:
32     IMPLICIT NONE
33     C == Global data ==
34     #include "SIZE.h"
35     #include "EEPARAMS.h"
36     #include "PARAMS.h"
37     #include "GRID.h"
38     #include "DYNVARS.h"
39     #include "FFIELDS.h"
40    
41     C !INPUT/OUTPUT PARAMETERS:
42     C == Routine arguments ==
43     C iMin,iMax :: Working range of x-index for applying forcing.
44     C jMin,jMax :: Working range of y-index for applying forcing.
45     C bi,bj :: Current tile indices
46     C kLev :: Current vertical level index
47     C myTime :: Current time in simulation
48     C myThid :: Thread Id number
49     INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
50     _RL myTime
51     INTEGER myThid
52    
53     C !LOCAL VARIABLES:
54     C == Local variables ==
55     C i,j :: Loop counters
56     C kSurface :: index of surface level
57     INTEGER i, j
58     INTEGER kSurface
59     CEOP
60    
61     IF ( fluidIsAir ) THEN
62     kSurface = 0
63     ELSEIF ( usingPCoords ) THEN
64     kSurface = Nr
65     ELSE
66     kSurface = 1
67     ENDIF
68    
69     C-- Forcing term
70     #ifdef ALLOW_AIM
71     IF ( useAIM ) CALL AIM_TENDENCY_APPLY_U(
72     & iMin,iMax, jMin,jMax, bi,bj, kLev,
73     & myTime, myThid )
74     #endif /* ALLOW_AIM */
75    
76     #ifdef ALLOW_ATM_PHYS
77     IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_U(
78     & iMin,iMax, jMin,jMax, bi,bj, kLev,
79     & myTime, myThid )
80     #endif /* ALLOW_ATM_PHYS */
81    
82     #ifdef ALLOW_FIZHI
83     IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_U(
84     & iMin,iMax, jMin,jMax, bi,bj, kLev,
85     & myTime, myThid )
86     #endif /* ALLOW_FIZHI */
87    
88     C Ocean: Add momentum surface forcing (e.g., wind-stress) in surface level
89     IF ( kLev .EQ. kSurface ) THEN
90     c DO j=1,sNy
91     C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNy+1]
92     DO j=0,sNy+1
93     DO i=1,sNx+1
94     gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)
95     & +foFacMom*surfaceForcingU(i,j,bi,bj)
96     & *recip_drF(kLev)*_recip_hFacW(i,j,kLev,bi,bj)
97     ENDDO
98     ENDDO
99     ELSEIF ( kSurface.EQ.-1 ) THEN
100     DO j=0,sNy+1
101     DO i=1,sNx+1
102     IF ( kSurfW(i,j,bi,bj).EQ.kLev ) THEN
103     gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)
104     & +foFacMom*surfaceForcingU(i,j,bi,bj)
105     & *recip_drF(kLev)*_recip_hFacW(i,j,kLev,bi,bj)
106     ENDIF
107     ENDDO
108     ENDDO
109     ENDIF
110    
111     #ifdef ALLOW_EDDYPSI
112     CALL TAUEDDY_EXTERNAL_FORCING_U(
113     I iMin,iMax, jMin,jMax, bi,bj, kLev,
114     I myTime, myThid )
115     #endif
116    
117     #ifdef ALLOW_RBCS
118     IF (useRBCS) THEN
119     CALL RBCS_ADD_TENDENCY( bi, bj, klev, -1,
120     & myTime, myThid )
121     ENDIF
122     #endif
123    
124     #ifdef ALLOW_OBCS
125     IF (useOBCS) THEN
126     CALL OBCS_SPONGE_U(
127     I iMin,iMax, jMin,jMax, bi,bj, kLev,
128     I myTime, myThid )
129     ENDIF
130     #endif
131    
132     #ifdef ALLOW_MYPACKAGE
133     IF ( useMYPACKAGE ) CALL MYPACKAGE_TENDENCY_APPLY_U(
134     & iMin,iMax, jMin,jMax, bi,bj, kLev,
135     & myTime, myThid )
136     #endif /* ALLOW_MYPACKAGE */
137    
138     RETURN
139     END
140    
141     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
142     CBOP
143     C !ROUTINE: EXTERNAL_FORCING_V
144     C !INTERFACE:
145     SUBROUTINE EXTERNAL_FORCING_V(
146     I iMin,iMax, jMin,jMax, bi,bj, kLev,
147     I myTime, myThid )
148     C !DESCRIPTION: \bv
149     C *==========================================================*
150     C | S/R EXTERNAL_FORCING_V
151     C | o Contains problem specific forcing for merid velocity.
152     C *==========================================================*
153     C | Adds terms to gV for forcing by external sources
154     C | e.g. wind stress, bottom friction etc ...
155     C *==========================================================*
156     C \ev
157    
158     C !USES:
159     IMPLICIT NONE
160     C == Global data ==
161     #include "SIZE.h"
162     #include "EEPARAMS.h"
163     #include "PARAMS.h"
164     #include "GRID.h"
165     #include "DYNVARS.h"
166     #include "FFIELDS.h"
167    
168     C !INPUT/OUTPUT PARAMETERS:
169     C == Routine arguments ==
170     C iMin,iMax :: Working range of x-index for applying forcing.
171     C jMin,jMax :: Working range of y-index for applying forcing.
172     C bi,bj :: Current tile indices
173     C kLev :: Current vertical level index
174     C myTime :: Current time in simulation
175     C myThid :: Thread Id number
176     INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
177     _RL myTime
178     INTEGER myThid
179    
180     C !LOCAL VARIABLES:
181     C == Local variables ==
182     C i,j :: Loop counters
183     C kSurface :: index of surface level
184     INTEGER i, j
185     INTEGER kSurface
186     CEOP
187    
188     IF ( fluidIsAir ) THEN
189     kSurface = 0
190     ELSEIF ( usingPCoords ) THEN
191     kSurface = Nr
192     ELSE
193     kSurface = 1
194     ENDIF
195    
196     C-- Forcing term
197     #ifdef ALLOW_AIM
198     IF ( useAIM ) CALL AIM_TENDENCY_APPLY_V(
199     & iMin,iMax, jMin,jMax, bi,bj, kLev,
200     & myTime, myThid )
201     #endif /* ALLOW_AIM */
202    
203     #ifdef ALLOW_ATM_PHYS
204     IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_V(
205     & iMin,iMax, jMin,jMax, bi,bj, kLev,
206     & myTime, myThid )
207     #endif /* ALLOW_ATM_PHYS */
208    
209     #ifdef ALLOW_FIZHI
210     IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_V(
211     & iMin,iMax, jMin,jMax, bi,bj, kLev,
212     & myTime, myThid )
213     #endif /* ALLOW_FIZHI */
214    
215     C Ocean: Add momentum surface forcing (e.g., wind-stress) in surface level
216     IF ( kLev .EQ. kSurface ) THEN
217     DO j=1,sNy+1
218     c DO i=1,sNx
219     C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNx+1]
220     DO i=0,sNx+1
221     gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)
222     & +foFacMom*surfaceForcingV(i,j,bi,bj)
223     & *recip_drF(kLev)*_recip_hFacS(i,j,kLev,bi,bj)
224     ENDDO
225     ENDDO
226     ELSEIF ( kSurface.EQ.-1 ) THEN
227     DO j=1,sNy+1
228     DO i=0,sNx+1
229     IF ( kSurfS(i,j,bi,bj).EQ.kLev ) THEN
230     gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)
231     & +foFacMom*surfaceForcingV(i,j,bi,bj)
232     & *recip_drF(kLev)*_recip_hFacS(i,j,kLev,bi,bj)
233     ENDIF
234     ENDDO
235     ENDDO
236     ENDIF
237    
238     #ifdef ALLOW_EDDYPSI
239     CALL TAUEDDY_EXTERNAL_FORCING_V(
240     I iMin,iMax, jMin,jMax, bi,bj, kLev,
241     I myTime, myThid )
242     #endif
243    
244     #ifdef ALLOW_RBCS
245     IF (useRBCS) THEN
246     CALL RBCS_ADD_TENDENCY( bi, bj, klev, -2,
247     & myTime, myThid )
248     ENDIF
249     #endif
250    
251     #ifdef ALLOW_OBCS
252     IF (useOBCS) THEN
253     CALL OBCS_SPONGE_V(
254     I iMin,iMax, jMin,jMax, bi,bj, kLev,
255     I myTime, myThid )
256     ENDIF
257     #endif
258    
259     #ifdef ALLOW_MYPACKAGE
260     IF ( useMYPACKAGE ) CALL MYPACKAGE_TENDENCY_APPLY_V(
261     & iMin,iMax, jMin,jMax, bi,bj, kLev,
262     & myTime, myThid )
263     #endif /* ALLOW_MYPACKAGE */
264    
265     RETURN
266     END
267    
268     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
269     CBOP
270     C !ROUTINE: EXTERNAL_FORCING_T
271     C !INTERFACE:
272     SUBROUTINE EXTERNAL_FORCING_T(
273     I iMin,iMax, jMin,jMax, bi,bj, kLev,
274     I myTime, myThid )
275     C !DESCRIPTION: \bv
276     C *==========================================================*
277     C | S/R EXTERNAL_FORCING_T
278     C | o Contains problem specific forcing for temperature.
279     C *==========================================================*
280     C | Adds terms to gT for forcing by external sources
281     C | e.g. heat flux, climatalogical relaxation, etc ...
282     C *==========================================================*
283     C \ev
284    
285     C !USES:
286     IMPLICIT NONE
287     C == Global data ==
288     #include "SIZE.h"
289     #include "EEPARAMS.h"
290     #include "PARAMS.h"
291     #include "GRID.h"
292     #include "DYNVARS.h"
293     #include "FFIELDS.h"
294     #include "SURFACE.h"
295    
296     C !INPUT/OUTPUT PARAMETERS:
297     C == Routine arguments ==
298     C iMin,iMax :: Working range of x-index for applying forcing.
299     C jMin,jMax :: Working range of y-index for applying forcing.
300     C bi,bj :: Current tile indices
301     C kLev :: Current vertical level index
302     C myTime :: Current time in simulation
303     C myThid :: Thread Id number
304     INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
305     _RL myTime
306     INTEGER myThid
307    
308     C !LOCAL VARIABLES:
309     C == Local variables ==
310     C i,j :: Loop counters
311     C kSurface :: index of surface level
312     INTEGER i, j
313     INTEGER kSurface
314     INTEGER km, kc, kp
315     _RL tmpVar(1:sNx,1:sNy)
316     _RL tmpFac, delPI
317     CEOP
318     c#ifdef ALLOW_FRICTION_HEATING
319     c _RL tmpFac
320     c#endif
321     #ifdef SHORTWAVE_HEATING
322     _RL minusone
323     PARAMETER (minusOne=-1.)
324     _RL swfracb(2)
325     INTEGER kp1
326     #endif
327    
328     IF ( fluidIsAir ) THEN
329     kSurface = 0
330     ELSEIF ( usingZCoords .AND. useShelfIce ) THEN
331     kSurface = -1
332     ELSEIF ( usingPCoords ) THEN
333     kSurface = Nr
334     ELSE
335     kSurface = 1
336     ENDIF
337    
338     C-- Forcing term
339     #ifdef ALLOW_AIM
340     IF ( useAIM ) CALL AIM_TENDENCY_APPLY_T(
341     & iMin,iMax, jMin,jMax, bi,bj, kLev,
342     & myTime, myThid )
343     #endif /* ALLOW_AIM */
344    
345     #ifdef ALLOW_ATM_PHYS
346     IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_T(
347     & iMin,iMax, jMin,jMax, bi,bj, kLev,
348     & myTime, myThid )
349     #endif /* ALLOW_ATM_PHYS */
350    
351     #ifdef ALLOW_FIZHI
352     IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_T(
353     & iMin,iMax, jMin,jMax, bi,bj, kLev,
354     & myTime, myThid )
355     #endif /* ALLOW_FIZHI */
356    
357     #ifdef ALLOW_ADDFLUID
358     IF ( selectAddFluid.NE.0 .AND. temp_addMass.NE.UNSET_RL ) THEN
359     IF ( ( selectAddFluid.GE.1 .AND. nonlinFreeSurf.GT.0 )
360     & .OR. convertFW2Salt.EQ.-1. _d 0 ) THEN
361     DO j=1,sNy
362     DO i=1,sNx
363     gT(i,j,kLev,bi,bj) = gT(i,j,kLev,bi,bj)
364     & + addMass(i,j,kLev,bi,bj)*mass2rUnit
365     & *( temp_addMass - theta(i,j,kLev,bi,bj) )
366     & *recip_rA(i,j,bi,bj)
367     & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
368     C & *recip_deepFac2C(kLev)*recip_rhoFacC(kLev)
369     ENDDO
370     ENDDO
371     ELSE
372     DO j=1,sNy
373     DO i=1,sNx
374     gT(i,j,kLev,bi,bj) = gT(i,j,kLev,bi,bj)
375     & + addMass(i,j,kLev,bi,bj)*mass2rUnit
376     & *( temp_addMass - tRef(kLev) )
377     & *recip_rA(i,j,bi,bj)
378     & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
379     C & *recip_deepFac2C(kLev)*recip_rhoFacC(kLev)
380     ENDDO
381     ENDDO
382     ENDIF
383     ENDIF
384     #endif /* ALLOW_ADDFLUID */
385    
386     #ifdef ALLOW_FRICTION_HEATING
387     IF ( addFrictionHeating ) THEN
388     IF ( fluidIsAir ) THEN
389     C conversion from in-situ Temp to Pot.Temp
390     tmpFac = (atm_Po/rC(kLev))**atm_kappa
391     C conversion from W/m^2/r_unit to K/s
392     tmpFac = (tmpFac/atm_Cp) * mass2rUnit
393     ELSE
394     C conversion from W/m^2/r_unit to K/s
395     tmpFac = recip_Cp * mass2rUnit
396     ENDIF
397     DO j=1,sNy
398     DO i=1,sNx
399     gT(i,j,kLev,bi,bj) = gT(i,j,kLev,bi,bj)
400     & + frictionHeating(i,j,kLev,bi,bj)
401     & *tmpFac*recip_rA(i,j,bi,bj)
402     & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
403     ENDDO
404     ENDDO
405     ENDIF
406     #endif /* ALLOW_FRICTION_HEATING */
407    
408     IF ( fluidIsAir .AND. atm_Rq.NE.zeroRL .AND. Nr.NE.1 ) THEN
409     C-- Compressible fluid: account for difference between moist and dry air
410     C specific volume in Enthalpy equation (+ V.dP term), since only the
411     C dry air part is accounted for in the (dry) Pot.Temp formulation.
412     C Used centered averaging from interface to center (consistent with
413     C conversion term in KE eq) and same discretisation ( [T*Q]_bar_k )
414     C as for Theta_v in CALC_PHI_HYD
415    
416     C conversion from in-situ Temp to Pot.Temp
417     tmpFac = (atm_Po/rC(kLev))**atm_kappa
418     C conversion from W/kg to K/s
419     tmpFac = tmpFac/atm_Cp
420     km = kLev-1
421     kc = kLev
422     kp = kLev+1
423     IF ( kLev.EQ.1 ) THEN
424     DO j=1,sNy
425     DO i=1,sNx
426     tmpVar(i,j) = 0.
427     ENDDO
428     ENDDO
429     ELSE
430     delPI = atm_Cp*( (rC(km)/atm_Po)**atm_kappa
431     & - (rC(kc)/atm_Po)**atm_kappa )
432     DO j=1,sNy
433     DO i=1,sNx
434     tmpVar(i,j) = wVel(i,j,kc,bi,bj)*delPI*atm_Rq
435     & *( theta(i,j,km,bi,bj)*salt(i,j,km,bi,bj)
436     & + theta(i,j,kc,bi,bj)*salt(i,j,kc,bi,bj)
437     & )*maskC(i,j,km,bi,bj)*0.25 _d 0
438     ENDDO
439     ENDDO
440     ENDIF
441     IF ( kLev.LT.Nr ) THEN
442     delPI = atm_Cp*( (rC(kc)/atm_Po)**atm_kappa
443     & - (rC(kp)/atm_Po)**atm_kappa )
444     DO j=1,sNy
445     DO i=1,sNx
446     tmpVar(i,j) = tmpVar(i,j)
447     & + wVel(i,j,kp,bi,bj)*delPI*atm_Rq
448     & *( theta(i,j,kc,bi,bj)*salt(i,j,kc,bi,bj)
449     & + theta(i,j,kp,bi,bj)*salt(i,j,kp,bi,bj)
450     & )*maskC(i,j,kp,bi,bj)*0.25 _d 0
451     ENDDO
452     ENDDO
453     ENDIF
454     DO j=1,sNy
455     DO i=1,sNx
456     gT(i,j,kLev,bi,bj) = gT(i,j,kLev,bi,bj)
457     & + tmpVar(i,j)*tmpFac
458     & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
459     ENDDO
460     ENDDO
461     #ifdef ALLOW_DIAGNOSTICS
462     IF ( useDiagnostics ) THEN
463     C conversion to W/m^2
464     tmpFac = rUnit2mass
465     CALL DIAGNOSTICS_SCALE_FILL( tmpVar, tmpFac, 1,
466     & 'MoistCor', kc, 1, 3, bi,bj,myThid )
467     ENDIF
468     #endif /* ALLOW_DIAGNOSTICS */
469     ENDIF
470    
471     C Ocean: Add temperature surface forcing (e.g., heat-flux) in surface level
472     IF ( kLev .EQ. kSurface ) THEN
473     DO j=1,sNy
474     DO i=1,sNx
475     gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
476     & +surfaceForcingT(i,j,bi,bj)
477     & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
478     ENDDO
479     ENDDO
480     ELSEIF ( kSurface.EQ.-1 ) THEN
481     DO j=1,sNy
482     DO i=1,sNx
483     IF ( kSurfC(i,j,bi,bj).EQ.kLev ) THEN
484     gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
485     & +surfaceForcingT(i,j,bi,bj)
486     & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
487     ENDIF
488     ENDDO
489     ENDDO
490     ENDIF
491    
492     IF (linFSConserveTr) THEN
493     DO j=1,sNy
494     DO i=1,sNx
495     IF (kLev .EQ. kSurfC(i,j,bi,bj)) THEN
496     gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
497     & +TsurfCor*recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
498     ENDIF
499     ENDDO
500     ENDDO
501     ENDIF
502    
503     #ifdef ALLOW_FRAZIL
504     IF ( useFRAZIL )
505     & CALL FRAZIL_TENDENCY_APPLY_T(
506     I iMin,iMax, jMin,jMax, bi,bj, kLev,
507     I myTime, myThid )
508     #endif /* ALLOW_FRAZIL */
509    
510     #ifdef ALLOW_SHELFICE
511     IF ( useShelfIce )
512     & CALL SHELFICE_FORCING_T(
513     I iMin,iMax, jMin,jMax, bi,bj, kLev,
514     I myTime, myThid )
515     #endif /* ALLOW_SHELFICE */
516    
517     #ifdef ALLOW_ICEFRONT
518     IF ( useICEFRONT )
519     & CALL ICEFRONT_TENDENCY_APPLY_T(
520     & bi,bj, kLev, myTime, myThid )
521     #endif /* ALLOW_ICEFRONT */
522    
523     #ifdef SHORTWAVE_HEATING
524     C Penetrating SW radiation
525     c IF ( usePenetratingSW ) THEN
526     swfracb(1)=abs(rF(klev))
527     swfracb(2)=abs(rF(klev+1))
528     CALL SWFRAC(
529     I 2, minusOne,
530     U swfracb,
531     I myTime, 1, myThid )
532     kp1 = klev+1
533     IF (klev.EQ.Nr) THEN
534     kp1 = klev
535     swfracb(2)=0. _d 0
536     ENDIF
537     DO j=1,sNy
538     DO i=1,sNx
539     gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj)
540     & -Qsw(i,j,bi,bj)*(swfracb(1)*maskC(i,j,klev,bi,bj)
541     & -swfracb(2)*maskC(i,j,kp1, bi,bj))
542     & *recip_Cp*mass2rUnit
543     & *recip_drF(klev)*_recip_hFacC(i,j,kLev,bi,bj)
544     ENDDO
545     ENDDO
546     c ENDIF
547     #endif
548    
549     #ifdef ALLOW_RBCS
550     IF (useRBCS) THEN
551     CALL RBCS_ADD_TENDENCY(bi,bj,klev, 1,
552     & myTime, myThid )
553     ENDIF
554     #endif
555    
556     #ifdef ALLOW_OBCS
557     IF (useOBCS) THEN
558     CALL OBCS_SPONGE_T(
559     I iMin,iMax, jMin,jMax, bi,bj, kLev,
560     I myTime, myThid )
561     ENDIF
562     #endif
563    
564     #ifdef ALLOW_BBL
565     IF ( useBBL ) CALL BBL_TENDENCY_APPLY_T(
566     & iMin,iMax, jMin,jMax, bi,bj, kLev,
567     & myTime, myThid )
568     #endif /* ALLOW_BBL */
569    
570     #ifdef ALLOW_MYPACKAGE
571     IF ( useMYPACKAGE ) CALL MYPACKAGE_TENDENCY_APPLY_T(
572     & iMin,iMax, jMin,jMax, bi,bj, kLev,
573     & myTime, myThid )
574     #endif /* ALLOW_MYPACKAGE */
575    
576     RETURN
577     END
578    
579     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
580     CBOP
581     C !ROUTINE: EXTERNAL_FORCING_S
582     C !INTERFACE:
583     SUBROUTINE EXTERNAL_FORCING_S(
584     I iMin,iMax, jMin,jMax, bi,bj, kLev,
585     I myTime, myThid )
586    
587     C !DESCRIPTION: \bv
588     C *==========================================================*
589     C | S/R EXTERNAL_FORCING_S
590     C | o Contains problem specific forcing for merid velocity.
591     C *==========================================================*
592     C | Adds terms to gS for forcing by external sources
593     C | e.g. fresh-water flux, climatalogical relaxation, etc ...
594     C *==========================================================*
595     C \ev
596    
597     C !USES:
598     IMPLICIT NONE
599     C == Global data ==
600     #include "SIZE.h"
601     #include "EEPARAMS.h"
602     #include "PARAMS.h"
603     #include "GRID.h"
604     #include "DYNVARS.h"
605     #include "FFIELDS.h"
606     #include "SURFACE.h"
607    
608     C !INPUT/OUTPUT PARAMETERS:
609     C == Routine arguments ==
610     C iMin,iMax :: Working range of x-index for applying forcing.
611     C jMin,jMax :: Working range of y-index for applying forcing.
612     C bi,bj :: Current tile indices
613     C kLev :: Current vertical level index
614     C myTime :: Current time in simulation
615     C myThid :: Thread Id number
616     INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
617     _RL myTime
618     INTEGER myThid
619    
620     C !LOCAL VARIABLES:
621     C == Local variables ==
622     C i,j :: Loop counters
623     C kSurface :: index of surface level
624     INTEGER i, j
625     INTEGER kSurface
626     CEOP
627    
628     IF ( fluidIsAir ) THEN
629     kSurface = 0
630     ELSEIF ( usingZCoords .AND. useShelfIce ) THEN
631     kSurface = -1
632     ELSEIF ( usingPCoords ) THEN
633     kSurface = Nr
634     ELSE
635     kSurface = 1
636     ENDIF
637    
638     C-- Forcing term
639     #ifdef ALLOW_AIM
640     IF ( useAIM ) CALL AIM_TENDENCY_APPLY_S(
641     & iMin,iMax, jMin,jMax, bi,bj, kLev,
642     & myTime, myThid )
643     #endif /* ALLOW_AIM */
644    
645     #ifdef ALLOW_ATM_PHYS
646     IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_S(
647     & iMin,iMax, jMin,jMax, bi,bj, kLev,
648     & myTime, myThid )
649     #endif /* ALLOW_ATM_PHYS */
650    
651     #ifdef ALLOW_FIZHI
652     IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_S(
653     & iMin,iMax, jMin,jMax, bi,bj, kLev,
654     & myTime, myThid )
655     #endif /* ALLOW_FIZHI */
656    
657     #ifdef ALLOW_ADDFLUID
658     IF ( selectAddFluid.NE.0 .AND. salt_addMass.NE.UNSET_RL ) THEN
659     IF ( ( selectAddFluid.GE.1 .AND. nonlinFreeSurf.GT.0 )
660     & .OR. convertFW2Salt.EQ.-1. _d 0 ) THEN
661     DO j=1,sNy
662     DO i=1,sNx
663     gS(i,j,kLev,bi,bj) = gS(i,j,kLev,bi,bj)
664     & + addMass(i,j,kLev,bi,bj)*mass2rUnit
665     & *( salt_addMass - salt(i,j,kLev,bi,bj) )
666     & *recip_rA(i,j,bi,bj)
667     & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
668     C & *recip_deepFac2C(kLev)*recip_rhoFacC(kLev)
669     ENDDO
670     ENDDO
671     ELSE
672     DO j=1,sNy
673     DO i=1,sNx
674     gS(i,j,kLev,bi,bj) = gS(i,j,kLev,bi,bj)
675     & + addMass(i,j,kLev,bi,bj)*mass2rUnit
676     & *( salt_addMass - sRef(kLev) )
677     & *recip_rA(i,j,bi,bj)
678     & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
679     C & *recip_deepFac2C(kLev)*recip_rhoFacC(kLev)
680     ENDDO
681     ENDDO
682     ENDIF
683     ENDIF
684     #endif /* ALLOW_ADDFLUID */
685    
686     C Ocean: Add salinity surface forcing (e.g., fresh-water) in surface level
687     IF ( kLev .EQ. kSurface ) THEN
688     DO j=1,sNy
689     DO i=1,sNx
690     gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
691     & +surfaceForcingS(i,j,bi,bj)
692     & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
693     ENDDO
694     ENDDO
695     ELSEIF ( kSurface.EQ.-1 ) THEN
696     DO j=1,sNy
697     DO i=1,sNx
698     IF ( kSurfC(i,j,bi,bj).EQ.kLev ) THEN
699     gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
700     & +surfaceForcingS(i,j,bi,bj)
701     & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
702     ENDIF
703     ENDDO
704     ENDDO
705     ENDIF
706    
707     IF (linFSConserveTr) THEN
708     DO j=1,sNy
709     DO i=1,sNx
710     IF (kLev .EQ. kSurfC(i,j,bi,bj)) THEN
711     gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
712     & +SsurfCor*recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
713     ENDIF
714     ENDDO
715     ENDDO
716     ENDIF
717    
718     #ifdef ALLOW_SHELFICE
719     IF ( useShelfIce )
720     & CALL SHELFICE_FORCING_S(
721     I iMin,iMax, jMin,jMax, bi,bj, kLev,
722     I myTime, myThid )
723     #endif /* ALLOW_SHELFICE */
724    
725     #ifdef ALLOW_ICEFRONT
726     IF ( useICEFRONT )
727     & CALL ICEFRONT_TENDENCY_APPLY_S(
728     & bi,bj, kLev, myTime, myThid )
729     #endif /* ALLOW_ICEFRONT */
730    
731 atn 1.2 Catn: org. version of SP: do within k-loop
732     Catn new version: outside k-loop; called from [temp,salt]_integrate.F
733 atn 1.1 #ifdef ALLOW_SALT_PLUME
734 atn 1.2 #ifndef SALT_PLUME_VOLUME
735 atn 1.1 IF ( useSALT_PLUME )
736     & CALL SALT_PLUME_TENDENCY_APPLY_S(
737     I iMin,iMax, jMin,jMax, bi,bj, kLev,
738     I myTime, myThid )
739 atn 1.2 C#ifdef SALT_PLUME_VOLUME
740     C IF ( useSALT_PLUME )
741     C & CALL SALT_PLUME_TENDENCY_APPLY_T(
742     C I iMin,iMax, jMin,jMax, bi,bj, kLev,
743     C I myTime, myThid )
744     #endif /* ndef SALT_PLUME_VOLUME */
745 atn 1.1 #endif /* ALLOW_SALT_PLUME */
746    
747     #ifdef ALLOW_RBCS
748     IF (useRBCS) THEN
749     CALL RBCS_ADD_TENDENCY(bi,bj,klev, 2,
750     & myTime, myThid )
751     ENDIF
752     #endif /* ALLOW_RBCS */
753    
754     #ifdef ALLOW_OBCS
755     IF (useOBCS) THEN
756     CALL OBCS_SPONGE_S(
757     I iMin,iMax, jMin,jMax, bi,bj, kLev,
758     I myTime, myThid )
759     ENDIF
760     #endif /* ALLOW_OBCS */
761    
762     #ifdef ALLOW_BBL
763     IF ( useBBL ) CALL BBL_TENDENCY_APPLY_S(
764     & iMin,iMax, jMin,jMax, bi,bj, kLev,
765     & myTime, myThid )
766     #endif /* ALLOW_BBL */
767    
768     #ifdef ALLOW_MYPACKAGE
769     IF ( useMYPACKAGE ) CALL MYPACKAGE_TENDENCY_APPLY_S(
770     & iMin,iMax, jMin,jMax, bi,bj, kLev,
771     & myTime, myThid )
772     #endif /* ALLOW_MYPACKAGE */
773    
774     RETURN
775     END
  ViewVC Help
Powered by ViewVC 1.1.22