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