/[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.9 by adcroft, Wed Nov 29 22:29:23 2000 UTC revision 1.66 by heimbach, Wed May 21 10:44:59 2014 UTC
# Line 1  Line 1 
1  C     $Header$  C $Header$
2    C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    #ifdef ALLOW_SALT_PLUME
7    #include "SALT_PLUME_OPTIONS.h"
8    #endif
9    
10  CStartOfInterface  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(        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     /==========================================================\  C     !DESCRIPTION: \bv
25  C     | S/R EXTERNAL_FORCING_U                                   |  C     *==========================================================*
26  C     | o Contains problem specific forcing for zonal velocity.  |  C     | S/R EXTERNAL_FORCING_U
27  C     |==========================================================|  C     | o Contains problem specific forcing for zonal velocity.
28  C     | Adds terms to gU for forcing by external sources         |  C     *==========================================================*
29  C     | e.g. wind stress, bottom friction etc..................  |  C     | Adds terms to gU for forcing by external sources
30  C     \==========================================================/  C     | e.g. wind stress, bottom friction etc ...
31        IMPLICIT NONE  C     *==========================================================*
32    C     \ev
33    
34    C     !USES:
35          IMPLICIT NONE
36  C     == Global data ==  C     == Global data ==
37  #include "SIZE.h"  #include "SIZE.h"
38  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 22  C     == Global data == Line 40  C     == Global data ==
40  #include "GRID.h"  #include "GRID.h"
41  #include "DYNVARS.h"  #include "DYNVARS.h"
42  #include "FFIELDS.h"  #include "FFIELDS.h"
43    
44    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
 CEndOfInterface  
55    
56    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          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  C--   Forcing term
73  C     Add windstress momentum impulse into the top-layer  #ifdef ALLOW_AIM
74        IF ( kLev .EQ. 1 ) THEN        IF ( useAIM ) CALL AIM_TENDENCY_APPLY_U(
75         DO j=jMin,jMax       &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
76          DO i=iMin,iMax       &                      myTime, myThid )
77           gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)  #endif /* ALLOW_AIM */
78       &   +foFacMom*surfaceTendencyU(i,j,bi,bj)  
79       &   *_maskW(i,j,kLev,bi,bj)  #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          ENDDO
111         ENDDO         ENDDO
112        ENDIF        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        RETURN
142        END        END
143  CStartOfInterface  
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(        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     /==========================================================\  C     !DESCRIPTION: \bv
152  C     | S/R EXTERNAL_FORCING_V                                   |  C     *==========================================================*
153  C     | o Contains problem specific forcing for merid velocity.  |  C     | S/R EXTERNAL_FORCING_V
154  C     |==========================================================|  C     | o Contains problem specific forcing for merid velocity.
155  C     | Adds terms to gV for forcing by external sources         |  C     *==========================================================*
156  C     | e.g. wind stress, bottom friction etc..................  |  C     | Adds terms to gV for forcing by external sources
157  C     \==========================================================/  C     | e.g. wind stress, bottom friction etc ...
158        IMPLICIT NONE  C     *==========================================================*
159    C     \ev
160    
161    C     !USES:
162          IMPLICIT NONE
163  C     == Global data ==  C     == Global data ==
164  #include "SIZE.h"  #include "SIZE.h"
165  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 72  C     == Global data == Line 168  C     == Global data ==
168  #include "DYNVARS.h"  #include "DYNVARS.h"
169  #include "FFIELDS.h"  #include "FFIELDS.h"
170    
171    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  CEndOfInterface  
183    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          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  C--   Forcing term
200  C     Add windstress momentum impulse into the top-layer  #ifdef ALLOW_AIM
201        IF ( kLev .EQ. 1 ) THEN        IF ( useAIM ) CALL AIM_TENDENCY_APPLY_V(
202         DO j=jMin,jMax       &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
203          DO i=iMin,iMax       &                      myTime, myThid )
204           gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)  #endif /* ALLOW_AIM */
205       &   +foFacMom*surfaceTendencyV(i,j,bi,bj)  
206       &   *_maskS(i,j,kLev,bi,bj)  #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          ENDDO
238         ENDDO         ENDDO
239        ENDIF        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        RETURN
269        END        END
270  CStartOfInterface  
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(        SUBROUTINE EXTERNAL_FORCING_T(
276       I           iMin, iMax, jMin, jMax,bi,bj,kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
277       I           maskC,       I           myTime, myThid )
278       I           myCurrentTime,myThid)  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        IMPLICIT NONE  C     \ev
287    
288    C     !USES:
289          IMPLICIT NONE
290  C     == Global data ==  C     == Global data ==
291  #include "SIZE.h"  #include "SIZE.h"
292  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 122  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:
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        _RS maskC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  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:
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          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  C--   Forcing term
342  C     Add heat in top-layer  #ifdef ALLOW_AIM
343        IF ( kLev .EQ. 1 ) THEN        IF ( useAIM ) CALL AIM_TENDENCY_APPLY_T(
344         DO j=jMin,jMax       &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
345          DO i=iMin,iMax       &                      myTime, myThid )
346           gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)  #endif /* ALLOW_AIM */
347       &     +maskC(i,j)*surfaceTendencyT(i,j,bi,bj)  
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
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          ENDDO
503         ENDDO         ENDDO
504        ENDIF        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       O     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)*Qsw(i,j,bi,bj)*(swfracb(1)-swfracb(2))          kp1 = klev
538       &    *recip_Cp*recip_rhoNil*recip_dRf(1)          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    #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        RETURN
580        END        END
581  CStartOfInterface  
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(        SUBROUTINE EXTERNAL_FORCING_S(
587       I           iMin, iMax, jMin, jMax,bi,bj,kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
588       I           maskC,       I           myTime, myThid )
589       I           myCurrentTime,myThid)  
590  C     /==========================================================\  C     !DESCRIPTION: \bv
591  C     | S/R EXTERNAL_FORCING_S                                   |  C     *==========================================================*
592  C     | o Contains problem specific forcing for merid velocity.  |  C     | S/R EXTERNAL_FORCING_S
593  C     |==========================================================|  C     | o Contains problem specific forcing for merid velocity.
594  C     | Adds terms to gS for forcing by external sources         |  C     *==========================================================*
595  C     | e.g. fresh-water flux, climatalogical relaxation.......  |  C     | Adds terms to gS for forcing by external sources
596  C     \==========================================================/  C     | e.g. fresh-water flux, climatalogical relaxation, etc ...
597        IMPLICIT NONE  C     *==========================================================*
598    C     \ev
599    
600    C     !USES:
601          IMPLICIT NONE
602  C     == Global data ==  C     == Global data ==
603  #include "SIZE.h"  #include "SIZE.h"
604  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 196  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:
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        _RS maskC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  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
 CEndOfInterface  
622    
623    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          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  C--   Forcing term
642  C     Add fresh-water in top-layer  #ifdef ALLOW_AIM
643        IF ( kLev .EQ. 1 ) THEN        IF ( useAIM ) CALL AIM_TENDENCY_APPLY_S(
644         DO j=jMin,jMax       &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
645          DO i=iMin,iMax       &                      myTime, myThid )
646           gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)  #endif /* ALLOW_AIM */
647       &   +maskC(i,j)*surfaceTendencyS(i,j,bi,bj)  
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          ENDDO
718         ENDDO         ENDDO
719        ENDIF        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
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        RETURN
778        END        END
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.66
  ViewVC Help
Powered by ViewVC 1.1.22