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