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

Diff of /MITgcm/model/src/external_forcing.F

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

revision 1.19.2.1 by edhill, Thu Oct 2 18:10:45 2003 UTC revision 1.66 by heimbach, Wed May 21 10:44:59 2014 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "PACKAGES_CONFIG.h"  #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6  #ifdef ALLOW_OBCS  #ifdef ALLOW_SALT_PLUME
7  # include "OBCS_OPTIONS.h"  #include "SALT_PLUME_OPTIONS.h"
8  #endif  #endif
9    
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  CBOP
19  C     !ROUTINE: EXTERNAL_FORCING_U  C     !ROUTINE: EXTERNAL_FORCING_U
20  C     !INTERFACE:  C     !INTERFACE:
21        SUBROUTINE EXTERNAL_FORCING_U(        SUBROUTINE EXTERNAL_FORCING_U(
22       I           iMin, iMax, jMin, jMax,bi,bj,kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
23       I           myCurrentTime,myThid)       I           myTime, myThid )
24  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
25  C     *==========================================================*  C     *==========================================================*
26  C     | S/R EXTERNAL_FORCING_U                                      C     | S/R EXTERNAL_FORCING_U
27  C     | o Contains problem specific forcing for zonal velocity.    C     | o Contains problem specific forcing for zonal velocity.
28  C     *==========================================================*  C     *==========================================================*
29  C     | Adds terms to gU for forcing by external sources            C     | Adds terms to gU for forcing by external sources
30  C     | e.g. wind stress, bottom friction etc..................    C     | e.g. wind stress, bottom friction etc ...
31  C     *==========================================================*  C     *==========================================================*
32  C     \ev  C     \ev
33    
# Line 35  C     == Global data == Line 43  C     == Global data ==
43    
44  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
45  C     == Routine arguments ==  C     == Routine arguments ==
46  C     iMin - Working range of tile for applying forcing.  C     iMin,iMax :: Working range of x-index for applying forcing.
47  C     iMax  C     jMin,jMax :: Working range of y-index for applying forcing.
48  C     jMin  C     bi,bj     :: Current tile indices
49  C     jMax  C     kLev      :: Current vertical level index
50  C     kLev  C     myTime    :: Current time in simulation
51    C     myThid    :: Thread Id number
52        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
53        _RL myCurrentTime        _RL myTime
54        INTEGER myThid        INTEGER myThid
55    
56  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
57  C     == Local variables ==  C     == Local variables ==
58  C     Loop counters  C     i,j       :: Loop counters
59        INTEGER I, J  C     kSurface  :: index of surface level
60  C     number of surface interface layer        INTEGER i, j
61        INTEGER kSurface        INTEGER kSurface
62  CEOP  CEOP
63    
64        if ( buoyancyRelation .eq. 'OCEANICP' ) then        IF ( fluidIsAir ) THEN
65           kSurface = 0
66          ELSEIF ( usingPCoords ) THEN
67         kSurface = Nr         kSurface = Nr
68        else        ELSE
69         kSurface = 1         kSurface = 1
70        endif        ENDIF
71    
72  C--   Forcing term  C--   Forcing term
73  C     Add windstress momentum impulse into the top-layer  #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        IF ( kLev .EQ. kSurface ) THEN
93         DO j=jMin,jMax  c      DO j=1,sNy
94          DO i=iMin,iMax  C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNy+1]
95           gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)         DO j=0,sNy+1
96       &   +foFacMom*surfaceTendencyU(i,j,bi,bj)          DO i=1,sNx+1
97       &   *_maskW(i,j,kLev,bi,bj)            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          ENDDO
111         ENDDO         ENDDO
112        ENDIF        ENDIF
113    
114  #if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE))  #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        IF (useOBCS) THEN
129         CALL OBCS_SPONGE_U(         CALL OBCS_SPONGE_U(
130       I           iMin, iMax, jMin, jMax,bi,bj,kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
131       I           myCurrentTime,myThid)       I           myTime, myThid )
132        ENDIF        ENDIF
133  #endif  #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        RETURN
142        END        END
143    
144    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
145  CBOP  CBOP
146  C     !ROUTINE: EXTERNAL_FORCING_V  C     !ROUTINE: EXTERNAL_FORCING_V
147  C     !INTERFACE:  C     !INTERFACE:
148        SUBROUTINE EXTERNAL_FORCING_V(        SUBROUTINE EXTERNAL_FORCING_V(
149       I           iMin, iMax, jMin, jMax,bi,bj,kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
150       I           myCurrentTime,myThid)       I           myTime, myThid )
151  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
152  C     *==========================================================*  C     *==========================================================*
153  C     | S/R EXTERNAL_FORCING_V                                      C     | S/R EXTERNAL_FORCING_V
154  C     | o Contains problem specific forcing for merid velocity.    C     | o Contains problem specific forcing for merid velocity.
155  C     *==========================================================*  C     *==========================================================*
156  C     | Adds terms to gV for forcing by external sources            C     | Adds terms to gV for forcing by external sources
157  C     | e.g. wind stress, bottom friction etc..................    C     | e.g. wind stress, bottom friction etc ...
158  C     *==========================================================*  C     *==========================================================*
159  C     \ev  C     \ev
160    
# Line 108  C     == Global data == Line 170  C     == Global data ==
170    
171  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
172  C     == Routine arguments ==  C     == Routine arguments ==
173  C     iMin - Working range of tile for applying forcing.  C     iMin,iMax :: Working range of x-index for applying forcing.
174  C     iMax  C     jMin,jMax :: Working range of y-index for applying forcing.
175  C     jMin  C     bi,bj     :: Current tile indices
176  C     jMax  C     kLev      :: Current vertical level index
177  C     kLev  C     myTime    :: Current time in simulation
178    C     myThid    :: Thread Id number
179        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
180        _RL myCurrentTime        _RL myTime
181        INTEGER myThid        INTEGER myThid
182    
183  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
184  C     == Local variables ==  C     == Local variables ==
185  C     Loop counters  C     i,j       :: Loop counters
186        INTEGER I, J  C     kSurface  :: index of surface level
187  C     number of surface interface layer        INTEGER i, j
188        INTEGER kSurface        INTEGER kSurface
189  CEOP  CEOP
190    
191        if ( buoyancyRelation .eq. 'OCEANICP' ) then        IF ( fluidIsAir ) THEN
192           kSurface = 0
193          ELSEIF ( usingPCoords ) THEN
194         kSurface = Nr         kSurface = Nr
195        else        ELSE
196         kSurface = 1         kSurface = 1
197        endif        ENDIF
198    
199  C--   Forcing term  C--   Forcing term
200  C     Add windstress momentum impulse into the top-layer  #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        IF ( kLev .EQ. kSurface ) THEN
220         DO j=jMin,jMax         DO j=1,sNy+1
221          DO i=iMin,iMax  c       DO i=1,sNx
222           gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)  C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNx+1]
223       &   +foFacMom*surfaceTendencyV(i,j,bi,bj)          DO i=0,sNx+1
224       &   *_maskS(i,j,kLev,bi,bj)            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          ENDDO
238         ENDDO         ENDDO
239        ENDIF        ENDIF
240    
241  #if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE))  #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        IF (useOBCS) THEN
256         CALL OBCS_SPONGE_V(         CALL OBCS_SPONGE_V(
257       I           iMin, iMax, jMin, jMax,bi,bj,kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
258       I           myCurrentTime,myThid)       I           myTime, myThid )
259        ENDIF        ENDIF
260  #endif  #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        RETURN
269        END        END
270    
271    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
272  CBOP  CBOP
273  C     !ROUTINE: EXTERNAL_FORCING_T  C     !ROUTINE: EXTERNAL_FORCING_T
274  C     !INTERFACE:  C     !INTERFACE:
275        SUBROUTINE EXTERNAL_FORCING_T(        SUBROUTINE EXTERNAL_FORCING_T(
276       I           iMin, iMax, jMin, jMax,bi,bj,kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
277       I           myCurrentTime,myThid)       I           myTime, myThid )
278  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
279  C     *==========================================================*  C     *==========================================================*
280  C     | S/R EXTERNAL_FORCING_T                                      C     | S/R EXTERNAL_FORCING_T
281  C     | o Contains problem specific forcing for temperature.        C     | o Contains problem specific forcing for temperature.
282  C     *==========================================================*  C     *==========================================================*
283  C     | Adds terms to gT for forcing by external sources            C     | Adds terms to gT for forcing by external sources
284  C     | e.g. heat flux, climatalogical relaxation..............    C     | e.g. heat flux, climatalogical relaxation, etc ...
285  C     *==========================================================*  C     *==========================================================*
286  C     \ev  C     \ev
287    
# Line 178  C     == Global data == Line 294  C     == Global data ==
294  #include "GRID.h"  #include "GRID.h"
295  #include "DYNVARS.h"  #include "DYNVARS.h"
296  #include "FFIELDS.h"  #include "FFIELDS.h"
297  #ifdef SHORTWAVE_HEATING  #include "SURFACE.h"
       integer two  
       _RL minusone  
       parameter (two=2,minusone=-1.)  
       _RL swfracb(two)  
 #endif  
298    
299  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
300  C     == Routine arguments ==  C     == Routine arguments ==
301  C     iMin - Working range of tile for applying forcing.  C     iMin,iMax :: Working range of x-index for applying forcing.
302  C     iMax  C     jMin,jMax :: Working range of y-index for applying forcing.
303  C     jMin  C     bi,bj     :: Current tile indices
304  C     jMax  C     kLev      :: Current vertical level index
305  C     kLev  C     myTime    :: Current time in simulation
306    C     myThid    :: Thread Id number
307        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
308        _RL myCurrentTime        _RL myTime
309        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
310    
311  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
312  C     == Local variables ==  C     == Local variables ==
313  C     Loop counters  C     i,j       :: Loop counters
314        INTEGER I, J  C     kSurface  :: index of surface level
315  C     number of surface interface layer        INTEGER i, j
316        INTEGER kSurface        INTEGER kSurface
317          INTEGER km, kc, kp
318          _RL tmpVar(1:sNx,1:sNy)
319          _RL tmpFac, delPI
320  CEOP  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 ( buoyancyRelation .eq. 'OCEANICP' ) then        IF ( fluidIsAir ) THEN
332           kSurface = 0
333          ELSEIF ( usingZCoords .AND. useShelfIce ) THEN
334           kSurface = -1
335          ELSEIF ( usingPCoords ) THEN
336         kSurface = Nr         kSurface = Nr
337        else        ELSE
338         kSurface = 1         kSurface = 1
339        endif        ENDIF
340    
341  C--   Forcing term  C--   Forcing term
342  C     Add heat in top-layer  #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 = mass2rUnit/HeatCapacity_Cp
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        IF ( kLev .EQ. kSurface ) THEN
476         DO j=jMin,jMax         DO j=1,sNy
477          DO i=iMin,iMax          DO i=1,sNx
478           gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)            gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
479       &     +maskC(i,j,kLev,bi,bj)*surfaceTendencyT(i,j,bi,bj)       &      +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          ENDDO
492         ENDDO         ENDDO
493        ENDIF        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  #ifdef SHORTWAVE_HEATING
527  C Penetrating SW radiation  C Penetrating SW radiation
528        swfracb(1)=abs(rF(klev))  c     IF ( usePenetratingSW ) THEN
529        swfracb(2)=abs(rF(klev+1))         swfracb(1)=abs(rF(klev))
530        call SWFRAC(         swfracb(2)=abs(rF(klev+1))
531       I     two,minusone,         CALL SWFRAC(
532       I     myCurrentTime,myThid,       I             2, minusOne,
533       U     swfracb)       U             swfracb,
534        DO j=jMin,jMax       I             myTime, 1, myThid )
535         DO i=iMin,iMax         kp1 = klev+1
536          gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj)         IF (klev.EQ.Nr) THEN
537       &   -maskC(i,j,klev,bi,bj)*Qsw(i,j,bi,bj)*(swfracb(1)-swfracb(2))          kp1 = klev
538       &    *recip_Cp*recip_rhoConst*recip_drF(klev)          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         &    *mass2rUnit/HeatCapacity_Cp
546         &    *recip_drF(klev)*_recip_hFacC(i,j,kLev,bi,bj)
547            ENDDO
548         ENDDO         ENDDO
549        ENDDO  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  #endif
558    
559  #if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE))  #ifdef ALLOW_OBCS
560        IF (useOBCS) THEN        IF (useOBCS) THEN
561         CALL OBCS_SPONGE_T(         CALL OBCS_SPONGE_T(
562       I           iMin, iMax, jMin, jMax,bi,bj,kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
563       I           myCurrentTime,myThid)       I           myTime, myThid )
564        ENDIF        ENDIF
565  #endif  #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        RETURN
580        END        END
581    
582    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
583  CBOP  CBOP
584  C     !ROUTINE: EXTERNAL_FORCING_S  C     !ROUTINE: EXTERNAL_FORCING_S
585  C     !INTERFACE:  C     !INTERFACE:
586        SUBROUTINE EXTERNAL_FORCING_S(        SUBROUTINE EXTERNAL_FORCING_S(
587       I           iMin, iMax, jMin, jMax,bi,bj,kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
588       I           myCurrentTime,myThid)       I           myTime, myThid )
589    
590  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
591  C     *==========================================================*  C     *==========================================================*
592  C     | S/R EXTERNAL_FORCING_S                                      C     | S/R EXTERNAL_FORCING_S
593  C     | o Contains problem specific forcing for merid velocity.    C     | o Contains problem specific forcing for merid velocity.
594  C     *==========================================================*  C     *==========================================================*
595  C     | Adds terms to gS for forcing by external sources            C     | Adds terms to gS for forcing by external sources
596  C     | e.g. fresh-water flux, climatalogical relaxation.......    C     | e.g. fresh-water flux, climatalogical relaxation, etc ...
597  C     *==========================================================*  C     *==========================================================*
598  C     \ev  C     \ev
599    
# Line 275  C     == Global data == Line 606  C     == Global data ==
606  #include "GRID.h"  #include "GRID.h"
607  #include "DYNVARS.h"  #include "DYNVARS.h"
608  #include "FFIELDS.h"  #include "FFIELDS.h"
609    #include "SURFACE.h"
610    
611  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
612  C     == Routine arguments ==  C     == Routine arguments ==
613  C     iMin - Working range of tile for applying forcing.  C     iMin,iMax :: Working range of x-index for applying forcing.
614  C     iMax  C     jMin,jMax :: Working range of y-index for applying forcing.
615  C     jMin  C     bi,bj     :: Current tile indices
616  C     jMax  C     kLev      :: Current vertical level index
617  C     kLev  C     myTime    :: Current time in simulation
618    C     myThid    :: Thread Id number
619        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj        INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
620        _RL myCurrentTime        _RL myTime
621        INTEGER myThid        INTEGER myThid
622    
623  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
624  C     == Local variables ==  C     == Local variables ==
625  C     Loop counters  C     i,j       :: Loop counters
626        INTEGER I, J  C     kSurface  :: index of surface level
627  C     number of surface interface layer        INTEGER i, j
628        INTEGER kSurface        INTEGER kSurface
629  CEOP  CEOP
630    
631        if ( buoyancyRelation .eq. 'OCEANICP' ) then        IF ( fluidIsAir ) THEN
632           kSurface = 0
633          ELSEIF ( usingZCoords .AND. useShelfIce ) THEN
634           kSurface = -1
635          ELSEIF ( usingPCoords ) THEN
636         kSurface = Nr         kSurface = Nr
637        else        ELSE
638         kSurface = 1         kSurface = 1
639        endif        ENDIF
   
640    
641  C--   Forcing term  C--   Forcing term
642  C     Add fresh-water in top-layer  #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        IF ( kLev .EQ. kSurface ) THEN
691         DO j=jMin,jMax         DO j=1,sNy
692          DO i=iMin,iMax          DO i=1,sNx
693           gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)            gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
694       &   +maskC(i,j,kLev,bi,bj)*surfaceTendencyS(i,j,bi,bj)       &      +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          ENDDO
707         ENDDO         ENDDO
708        ENDIF        ENDIF
709    
710  #if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE))        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    Catn: org. version of SP: do within k-loop
735    Catn  new version: outside k-loop; called from [temp,salt]_integrate.F
736    #ifdef ALLOW_SALT_PLUME
737    CC#ifndef SALT_PLUME_VOLUME
738          IF ( useSALT_PLUME )
739         &     CALL SALT_PLUME_TENDENCY_APPLY_S(
740         I     iMin,iMax, jMin,jMax, bi,bj, kLev,
741         I     myTime, myThid )
742    #ifdef SALT_PLUME_VOLUME
743          IF ( useSALT_PLUME )
744         &     CALL SALT_PLUME_TENDENCY_APPLY_T(
745         I     iMin,iMax, jMin,jMax, bi,bj, kLev,
746         I     myTime, myThid )
747    #endif /* ndef SALT_PLUME_VOLUME */
748    #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        IF (useOBCS) THEN
759         CALL OBCS_SPONGE_S(         CALL OBCS_SPONGE_S(
760       I           iMin, iMax, jMin, jMax,bi,bj,kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
761       I           myCurrentTime,myThid)       I           myTime, myThid )
762        ENDIF        ENDIF
763  #endif  #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        RETURN
778        END        END
Legend:
Removed from v.1.19.2.1  
changed lines
  Added in v.1.66
  ViewVC Help
Powered by ViewVC 1.1.22