/[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.35 by stephd, Mon Dec 19 19:09:35 2005 UTC revision 1.68 by atn, Thu May 22 22:00:36 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_EXF  #ifdef ALLOW_SALT_PLUME
7  # include "EXF_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:
# Line 48  C     myThid    :: Thread Id number Line 56  C     myThid    :: Thread Id number
56  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
57  C     == Local variables ==  C     == Local variables ==
58  C     i,j       :: Loop counters  C     i,j       :: Loop counters
59  C     kSurface  :: index of surface layer  C     kSurface  :: index of surface level
60        INTEGER i, j        INTEGER i, j
61        INTEGER kSurface        INTEGER kSurface
62  CEOP  CEOP
# Line 68  C--   Forcing term Line 76  C--   Forcing term
76       &                      myTime, myThid )       &                      myTime, myThid )
77  #endif /* ALLOW_AIM */  #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  #ifdef ALLOW_FIZHI
86        IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_U(        IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_U(
87       &                      iMin,iMax, jMin,jMax, bi,bj, kLev,       &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
88       &                      myTime, myThid )       &                      myTime, myThid )
89  #endif /* ALLOW_FIZHI */  #endif /* ALLOW_FIZHI */
90    
91  C     Add windstress momentum impulse into the top-layer  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  c      DO j=1,sNy  c      DO j=1,sNy
94  C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNy+1]  C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNy+1]
95         DO j=0,sNy+1         DO j=0,sNy+1
96          DO i=1,sNx+1          DO i=1,sNx+1
97           gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)            gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)
98       &   +foFacMom*surfaceForcingU(i,j,bi,bj)       &      +foFacMom*surfaceForcingU(i,j,bi,bj)
99       &   *recip_drF(kLev)*recip_hFacW(i,j,kLev,bi,bj)       &      *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_TAU_EDDY))  #ifdef ALLOW_EDDYPSI
115         CALL TAUEDDY_EXTERNAL_FORCING_U(         CALL TAUEDDY_EXTERNAL_FORCING_U(
116       I           iMin,iMax, jMin,jMax, bi,bj, kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
117       I           myTime, myThid )       I           myTime, myThid )
118  #endif  #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  #ifdef ALLOW_OBCS
128        IF (useOBCS) THEN        IF (useOBCS) THEN
129         CALL OBCS_SPONGE_U(         CALL OBCS_SPONGE_U(
# Line 101  C-jmc: Without CD-scheme, this is OK ; b Line 132  C-jmc: Without CD-scheme, this is OK ; b
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    
# Line 146  C     myThid    :: Thread Id number Line 183  C     myThid    :: Thread Id number
183  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
184  C     == Local variables ==  C     == Local variables ==
185  C     i,j       :: Loop counters  C     i,j       :: Loop counters
186  C     kSurface  :: index of surface layer  C     kSurface  :: index of surface level
187        INTEGER i, j        INTEGER i, j
188        INTEGER kSurface        INTEGER kSurface
189  CEOP  CEOP
# Line 166  C--   Forcing term Line 203  C--   Forcing term
203       &                      myTime, myThid )       &                      myTime, myThid )
204  #endif /* ALLOW_AIM */  #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  #ifdef ALLOW_FIZHI
213        IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_V(        IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_V(
214       &                      iMin,iMax, jMin,jMax, bi,bj, kLev,       &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
215       &                      myTime, myThid )       &                      myTime, myThid )
216  #endif /* ALLOW_FIZHI */  #endif /* ALLOW_FIZHI */
217    
218  C     Add windstress momentum impulse into the top-layer  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=1,sNy+1         DO j=1,sNy+1
221  c       DO i=1,sNx  c       DO i=1,sNx
222  C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNx+1]  C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNx+1]
223          DO i=0,sNx+1          DO i=0,sNx+1
224           gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)            gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)
225       &   +foFacMom*surfaceForcingV(i,j,bi,bj)       &      +foFacMom*surfaceForcingV(i,j,bi,bj)
226       &   *recip_drF(kLev)*recip_hFacS(i,j,kLev,bi,bj)       &      *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_TAU_EDDY))  #ifdef ALLOW_EDDYPSI
242         CALL TAUEDDY_EXTERNAL_FORCING_V(         CALL TAUEDDY_EXTERNAL_FORCING_V(
243       I           iMin,iMax, jMin,jMax, bi,bj, kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
244       I           myTime, myThid )       I           myTime, myThid )
245  #endif  #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  #ifdef ALLOW_OBCS
255        IF (useOBCS) THEN        IF (useOBCS) THEN
256         CALL OBCS_SPONGE_V(         CALL OBCS_SPONGE_V(
# Line 199  C-jmc: Without CD-scheme, this is OK ; b Line 259  C-jmc: Without CD-scheme, this is OK ; b
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    
# Line 228  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    #include "SURFACE.h"
298    
299  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
300  C     == Routine arguments ==  C     == Routine arguments ==
# Line 244  C     myThid    :: Thread Id number Line 311  C     myThid    :: Thread Id number
311  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
312  C     == Local variables ==  C     == Local variables ==
313  C     i,j       :: Loop counters  C     i,j       :: Loop counters
314  C     kSurface  :: index of surface layer  C     kSurface  :: index of surface level
315        INTEGER i, j        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          _RL recip_Cp
321  CEOP  CEOP
322    c#ifdef ALLOW_FRICTION_HEATING
323    c     _RL tmpFac
324    c#endif
325  #ifdef SHORTWAVE_HEATING  #ifdef SHORTWAVE_HEATING
       integer two  
326        _RL minusone        _RL minusone
327        parameter (two=2,minusone=-1.)        PARAMETER (minusOne=-1.)
328        _RL swfracb(two)        _RL swfracb(2)
329        INTEGER kp1        INTEGER kp1
330  #endif  #endif
331    
332        IF ( fluidIsAir ) THEN        IF ( fluidIsAir ) THEN
333         kSurface = 0         kSurface = 0
334          ELSEIF ( usingZCoords .AND. useShelfIce ) THEN
335           kSurface = -1
336        ELSEIF ( usingPCoords ) THEN        ELSEIF ( usingPCoords ) THEN
337         kSurface = Nr         kSurface = Nr
338        ELSE        ELSE
339         kSurface = 1         kSurface = 1
340        ENDIF        ENDIF
341          recip_Cp = 1. _d 0 / HeatCapacity_Cp
342    
343  C--   Forcing term  C--   Forcing term
344  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
# Line 271  C--   Forcing term Line 347  C--   Forcing term
347       &                      myTime, myThid )       &                      myTime, myThid )
348  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
349    
350    #ifdef ALLOW_ATM_PHYS
351          IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_T(
352         &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
353         &                      myTime, myThid )
354    #endif /* ALLOW_ATM_PHYS */
355    
356  #ifdef ALLOW_FIZHI  #ifdef ALLOW_FIZHI
357        IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_T(        IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_T(
358       &                      iMin,iMax, jMin,jMax, bi,bj, kLev,       &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
359       &                      myTime, myThid )       &                      myTime, myThid )
360  #endif /* ALLOW_FIZHI */  #endif /* ALLOW_FIZHI */
361    
362  C     Add heat in top-layer  #ifdef ALLOW_ADDFLUID
363          IF ( selectAddFluid.NE.0 .AND. temp_addMass.NE.UNSET_RL ) THEN
364           IF ( ( selectAddFluid.GE.1 .AND. nonlinFreeSurf.GT.0 )
365         &      .OR. convertFW2Salt.EQ.-1. _d 0 ) THEN
366             DO j=1,sNy
367              DO i=1,sNx
368                gT(i,j,kLev,bi,bj) = gT(i,j,kLev,bi,bj)
369         &        + addMass(i,j,kLev,bi,bj)*mass2rUnit
370         &          *( temp_addMass - theta(i,j,kLev,bi,bj) )
371         &          *recip_rA(i,j,bi,bj)
372         &          *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
373    C    &          *recip_deepFac2C(kLev)*recip_rhoFacC(kLev)
374              ENDDO
375             ENDDO
376           ELSE
377             DO j=1,sNy
378              DO i=1,sNx
379                gT(i,j,kLev,bi,bj) = gT(i,j,kLev,bi,bj)
380         &        + addMass(i,j,kLev,bi,bj)*mass2rUnit
381         &          *( temp_addMass - tRef(kLev) )
382         &          *recip_rA(i,j,bi,bj)
383         &          *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
384    C    &          *recip_deepFac2C(kLev)*recip_rhoFacC(kLev)
385              ENDDO
386             ENDDO
387           ENDIF
388          ENDIF
389    #endif /* ALLOW_ADDFLUID */
390    
391    #ifdef ALLOW_FRICTION_HEATING
392          IF ( addFrictionHeating ) THEN
393            IF ( fluidIsAir ) THEN
394    C         conversion from in-situ Temp to Pot.Temp
395              tmpFac = (atm_Po/rC(kLev))**atm_kappa
396    C         conversion from W/m^2/r_unit to K/s
397              tmpFac = (tmpFac/atm_Cp) * mass2rUnit
398            ELSE
399    C         conversion from W/m^2/r_unit to K/s
400              tmpFac = recip_Cp * mass2rUnit
401            ENDIF
402            DO j=1,sNy
403              DO i=1,sNx
404                gT(i,j,kLev,bi,bj) = gT(i,j,kLev,bi,bj)
405         &         + frictionHeating(i,j,kLev,bi,bj)
406         &          *tmpFac*recip_rA(i,j,bi,bj)
407         &          *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
408              ENDDO
409            ENDDO
410          ENDIF
411    #endif /* ALLOW_FRICTION_HEATING */
412    
413          IF ( fluidIsAir .AND. atm_Rq.NE.zeroRL .AND. Nr.NE.1 ) THEN
414    C--   Compressible fluid: account for difference between moist and dry air
415    C     specific volume in Enthalpy equation (+ V.dP term), since only the
416    C     dry air part is accounted for in the (dry) Pot.Temp formulation.
417    C     Used centered averaging from interface to center (consistent with
418    C     conversion term in KE eq) and same discretisation ( [T*Q]_bar_k )
419    C     as for Theta_v in CALC_PHI_HYD
420    
421    C     conversion from in-situ Temp to Pot.Temp
422            tmpFac = (atm_Po/rC(kLev))**atm_kappa
423    C     conversion from W/kg to K/s
424            tmpFac = tmpFac/atm_Cp
425            km = kLev-1
426            kc = kLev
427            kp = kLev+1
428            IF ( kLev.EQ.1 ) THEN
429              DO j=1,sNy
430               DO i=1,sNx
431                tmpVar(i,j) = 0.
432               ENDDO
433              ENDDO
434            ELSE
435              delPI = atm_Cp*( (rC(km)/atm_Po)**atm_kappa
436         &                   - (rC(kc)/atm_Po)**atm_kappa )
437              DO j=1,sNy
438               DO i=1,sNx
439                tmpVar(i,j) = wVel(i,j,kc,bi,bj)*delPI*atm_Rq
440         &                  *( theta(i,j,km,bi,bj)*salt(i,j,km,bi,bj)
441         &                   + theta(i,j,kc,bi,bj)*salt(i,j,kc,bi,bj)
442         &                   )*maskC(i,j,km,bi,bj)*0.25 _d 0
443               ENDDO
444              ENDDO
445            ENDIF
446            IF ( kLev.LT.Nr ) THEN
447              delPI = atm_Cp*( (rC(kc)/atm_Po)**atm_kappa
448         &                   - (rC(kp)/atm_Po)**atm_kappa )
449              DO j=1,sNy
450               DO i=1,sNx
451                tmpVar(i,j) = tmpVar(i,j)
452         &                  + wVel(i,j,kp,bi,bj)*delPI*atm_Rq
453         &                  *( theta(i,j,kc,bi,bj)*salt(i,j,kc,bi,bj)
454         &                   + theta(i,j,kp,bi,bj)*salt(i,j,kp,bi,bj)
455         &                   )*maskC(i,j,kp,bi,bj)*0.25 _d 0
456               ENDDO
457              ENDDO
458            ENDIF
459            DO j=1,sNy
460              DO i=1,sNx
461                gT(i,j,kLev,bi,bj) = gT(i,j,kLev,bi,bj)
462         &         + tmpVar(i,j)*tmpFac
463         &          *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
464              ENDDO
465            ENDDO
466    #ifdef ALLOW_DIAGNOSTICS
467            IF ( useDiagnostics ) THEN
468    C     conversion to W/m^2
469              tmpFac = rUnit2mass
470              CALL DIAGNOSTICS_SCALE_FILL( tmpVar, tmpFac, 1,
471         &                     'MoistCor', kc, 1, 3, bi,bj,myThid )
472            ENDIF
473    #endif /* ALLOW_DIAGNOSTICS */
474          ENDIF
475    
476    C     Ocean: Add temperature surface forcing (e.g., heat-flux) in surface level
477        IF ( kLev .EQ. kSurface ) THEN        IF ( kLev .EQ. kSurface ) THEN
478         DO j=1,sNy         DO j=1,sNy
479          DO i=1,sNx          DO i=1,sNx
480           gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)            gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
481       &     +surfaceForcingT(i,j,bi,bj)       &      +surfaceForcingT(i,j,bi,bj)
482       &     *recip_drF(kLev)*recip_hFacC(i,j,kLev,bi,bj)       &      *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
483            ENDDO
484           ENDDO
485          ELSEIF ( kSurface.EQ.-1 ) THEN
486           DO j=1,sNy
487            DO i=1,sNx
488             IF ( kSurfC(i,j,bi,bj).EQ.kLev ) THEN
489              gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
490         &      +surfaceForcingT(i,j,bi,bj)
491         &      *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
492             ENDIF
493          ENDDO          ENDDO
494         ENDDO         ENDDO
495        ENDIF        ENDIF
496    
497          IF (linFSConserveTr) THEN
498           DO j=1,sNy
499            DO i=1,sNx
500              IF (kLev .EQ. kSurfC(i,j,bi,bj)) THEN
501                gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
502         &        +TsurfCor*recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
503              ENDIF
504            ENDDO
505           ENDDO
506          ENDIF
507    
508    #ifdef ALLOW_FRAZIL
509          IF ( useFRAZIL )
510         &     CALL FRAZIL_TENDENCY_APPLY_T(
511         I     iMin,iMax, jMin,jMax, bi,bj, kLev,
512         I     myTime, myThid )
513    #endif /* ALLOW_FRAZIL */
514    
515    #ifdef ALLOW_SHELFICE
516          IF ( useShelfIce )
517         &     CALL SHELFICE_FORCING_T(
518         I     iMin,iMax, jMin,jMax, bi,bj, kLev,
519         I     myTime, myThid )
520    #endif /* ALLOW_SHELFICE */
521    
522    #ifdef ALLOW_ICEFRONT
523          IF ( useICEFRONT )
524         &     CALL ICEFRONT_TENDENCY_APPLY_T(
525         &     bi,bj, kLev, myTime, myThid )
526    #endif /* ALLOW_ICEFRONT */
527    
528  #ifdef SHORTWAVE_HEATING  #ifdef SHORTWAVE_HEATING
529  C Penetrating SW radiation  C Penetrating SW radiation
530  c     IF ( usePenetratingSW ) THEN  c     IF ( usePenetratingSW ) THEN
531         swfracb(1)=abs(rF(klev))         swfracb(1)=abs(rF(klev))
532         swfracb(2)=abs(rF(klev+1))         swfracb(2)=abs(rF(klev+1))
533         CALL SWFRAC(         CALL SWFRAC(
534       I     two,minusone,       I             2, minusOne,
535       I     myTime,myThid,       U             swfracb,
536       U     swfracb)       I             myTime, 1, myThid )
537         kp1 = klev+1         kp1 = klev+1
538         IF (klev.EQ.Nr) THEN         IF (klev.EQ.Nr) THEN
539          kp1 = klev          kp1 = klev
# Line 307  c     IF ( usePenetratingSW ) THEN Line 544  c     IF ( usePenetratingSW ) THEN
544           gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj)           gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj)
545       &   -Qsw(i,j,bi,bj)*(swfracb(1)*maskC(i,j,klev,bi,bj)       &   -Qsw(i,j,bi,bj)*(swfracb(1)*maskC(i,j,klev,bi,bj)
546       &                   -swfracb(2)*maskC(i,j,kp1, bi,bj))       &                   -swfracb(2)*maskC(i,j,kp1, bi,bj))
547       &    *recip_Cp*recip_rhoConst       &    *recip_Cp*mass2rUnit
548       &    *recip_drF(klev)*recip_hFacC(i,j,kLev,bi,bj)       &    *recip_drF(klev)*_recip_hFacC(i,j,kLev,bi,bj)
549          ENDDO          ENDDO
550         ENDDO         ENDDO
551  c     ENDIF  c     ENDIF
552  #endif  #endif
553    
554    #ifdef ALLOW_SALT_PLUME
555    #ifdef SALT_PLUME_VOLUME
556          IF ( useSALT_PLUME )
557         &     CALL SALT_PLUME_TENDENCY_APPLY_T(
558         I     iMin,iMax, jMin,jMax, bi,bj, kLev,
559         I     myTime, myThid )
560    #endif /* SALT_PLUME_VOLUME */
561    #endif /* ALLOW_SALT_PLUME */
562    
563  #ifdef ALLOW_RBCS  #ifdef ALLOW_RBCS
564         if (useRBCS) then         IF (useRBCS) THEN
565            call RBCS_ADD_TENDENCY(bi,bj,klev, 1,            CALL RBCS_ADD_TENDENCY(bi,bj,klev, 1,
566       &                            myTime, myThid )       &                            myTime, myThid )
        endif  
 #endif  
   
 #ifdef ALLOW_CLIMTEMP_RELAXATION  
        IF ( tauThetaClimRelax3Dim .NE. 0. ) THEN  
         DO j=1,sNy  
          DO i=1,sNx  
           gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj)  
      &     -1./tauThetaClimRelax3Dim  
      &         *(theta(i,j,klev,bi,bj)-thetaStar(i,j,klev,bi,bj))  
      &         *hFacC(i,j,klev,bi,bj)*recip_hFacC(i,j,kLev,bi,bj)  
          ENDDO  
         ENDDO  
567         ENDIF         ENDIF
568  #endif  #endif
569    
# Line 342  c     ENDIF Line 575  c     ENDIF
575        ENDIF        ENDIF
576  #endif  #endif
577    
578    #ifdef ALLOW_BBL
579          IF ( useBBL ) CALL BBL_TENDENCY_APPLY_T(
580         &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
581         &                      myTime, myThid )
582    #endif /* ALLOW_BBL */
583    
584    #ifdef ALLOW_MYPACKAGE
585          IF ( useMYPACKAGE ) CALL MYPACKAGE_TENDENCY_APPLY_T(
586         &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
587         &                      myTime, myThid )
588    #endif /* ALLOW_MYPACKAGE */
589    
590        RETURN        RETURN
591        END        END
592    
# Line 372  C     == Global data == Line 617  C     == Global data ==
617  #include "GRID.h"  #include "GRID.h"
618  #include "DYNVARS.h"  #include "DYNVARS.h"
619  #include "FFIELDS.h"  #include "FFIELDS.h"
620    #include "SURFACE.h"
621    
622  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
623  C     == Routine arguments ==  C     == Routine arguments ==
# Line 388  C     myThid    :: Thread Id number Line 634  C     myThid    :: Thread Id number
634  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
635  C     == Local variables ==  C     == Local variables ==
636  C     i,j       :: Loop counters  C     i,j       :: Loop counters
637  C     kSurface  :: index of surface layer  C     kSurface  :: index of surface level
638        INTEGER i, j        INTEGER i, j
639        INTEGER kSurface        INTEGER kSurface
640  CEOP  CEOP
641    
642        IF ( fluidIsAir ) THEN        IF ( fluidIsAir ) THEN
643         kSurface = 0         kSurface = 0
644          ELSEIF ( usingZCoords .AND. useShelfIce ) THEN
645           kSurface = -1
646        ELSEIF ( usingPCoords ) THEN        ELSEIF ( usingPCoords ) THEN
647         kSurface = Nr         kSurface = Nr
648        ELSE        ELSE
# Line 408  C--   Forcing term Line 656  C--   Forcing term
656       &                      myTime, myThid )       &                      myTime, myThid )
657  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
658    
659    #ifdef ALLOW_ATM_PHYS
660          IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_S(
661         &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
662         &                      myTime, myThid )
663    #endif /* ALLOW_ATM_PHYS */
664    
665  #ifdef ALLOW_FIZHI  #ifdef ALLOW_FIZHI
666        IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_S(        IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_S(
667       &                      iMin,iMax, jMin,jMax, bi,bj, kLev,       &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
668       &                      myTime, myThid )       &                      myTime, myThid )
669  #endif /* ALLOW_FIZHI */  #endif /* ALLOW_FIZHI */
670    
671  C     Add fresh-water in top-layer  #ifdef ALLOW_ADDFLUID
672          IF ( selectAddFluid.NE.0 .AND. salt_addMass.NE.UNSET_RL ) THEN
673           IF ( ( selectAddFluid.GE.1 .AND. nonlinFreeSurf.GT.0 )
674         &      .OR. convertFW2Salt.EQ.-1. _d 0 ) THEN
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 - salt(i,j,kLev,bi,bj) )
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           ELSE
686             DO j=1,sNy
687              DO i=1,sNx
688                gS(i,j,kLev,bi,bj) = gS(i,j,kLev,bi,bj)
689         &        + addMass(i,j,kLev,bi,bj)*mass2rUnit
690         &          *( salt_addMass - sRef(kLev) )
691         &          *recip_rA(i,j,bi,bj)
692         &          *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
693    C    &          *recip_deepFac2C(kLev)*recip_rhoFacC(kLev)
694              ENDDO
695             ENDDO
696           ENDIF
697          ENDIF
698    #endif /* ALLOW_ADDFLUID */
699    
700    C     Ocean: Add salinity surface forcing (e.g., fresh-water) in surface level
701        IF ( kLev .EQ. kSurface ) THEN        IF ( kLev .EQ. kSurface ) THEN
702         DO j=1,sNy         DO j=1,sNy
703          DO i=1,sNx          DO i=1,sNx
704           gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)            gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
705       &     +surfaceForcingS(i,j,bi,bj)       &      +surfaceForcingS(i,j,bi,bj)
706       &     *recip_drF(kLev)*recip_hFacC(i,j,kLev,bi,bj)       &      *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
707            ENDDO
708           ENDDO
709          ELSEIF ( kSurface.EQ.-1 ) THEN
710           DO j=1,sNy
711            DO i=1,sNx
712             IF ( kSurfC(i,j,bi,bj).EQ.kLev ) THEN
713              gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
714         &      +surfaceForcingS(i,j,bi,bj)
715         &      *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
716             ENDIF
717          ENDDO          ENDDO
718         ENDDO         ENDDO
719        ENDIF        ENDIF
720    
721          IF (linFSConserveTr) THEN
722           DO j=1,sNy
723            DO i=1,sNx
724              IF (kLev .EQ. kSurfC(i,j,bi,bj)) THEN
725                gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
726         &        +SsurfCor*recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
727              ENDIF
728            ENDDO
729           ENDDO
730          ENDIF
731    
732    #ifdef ALLOW_SHELFICE
733          IF ( useShelfIce )
734         &     CALL SHELFICE_FORCING_S(
735         I     iMin,iMax, jMin,jMax, bi,bj, kLev,
736         I     myTime, myThid )
737    #endif /* ALLOW_SHELFICE */
738    
739    #ifdef ALLOW_ICEFRONT
740          IF ( useICEFRONT )
741         &     CALL ICEFRONT_TENDENCY_APPLY_S(
742         &     bi,bj, kLev, myTime, myThid )
743    #endif /* ALLOW_ICEFRONT */
744    
745    #ifdef ALLOW_SALT_PLUME
746          IF ( useSALT_PLUME )
747         &     CALL SALT_PLUME_TENDENCY_APPLY_S(
748         I     iMin,iMax, jMin,jMax, bi,bj, kLev,
749         I     myTime, myThid )
750    #endif /* ALLOW_SALT_PLUME */
751    
752  #ifdef ALLOW_RBCS  #ifdef ALLOW_RBCS
753         if (useRBCS) then         IF (useRBCS) THEN
754            call RBCS_ADD_TENDENCY(bi,bj,klev, 2,            CALL RBCS_ADD_TENDENCY(bi,bj,klev, 2,
755       &                            myTime, myThid )       &                            myTime, myThid )
        endif  
 #endif  
   
 #ifdef ALLOW_CLIMSALT_RELAXATION  
        IF ( tauSaltClimRelax3Dim .NE. 0. ) THEN  
         DO j=1,sNy  
          DO i=1,sNx  
           gS(i,j,klev,bi,bj) = gS(i,j,klev,bi,bj)  
      &     -1./tauSaltClimRelax3Dim  
      &         *(salt(i,j,klev,bi,bj)-saltStar(i,j,klev,bi,bj))  
      &         *hFacC(i,j,klev,bi,bj)*recip_hFacC(i,j,kLev,bi,bj)  
          ENDDO  
         ENDDO  
756         ENDIF         ENDIF
757  #endif  #endif /* ALLOW_RBCS */
758    
759  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
760        IF (useOBCS) THEN        IF (useOBCS) THEN
# Line 451  C     Add fresh-water in top-layer Line 762  C     Add fresh-water in top-layer
762       I           iMin,iMax, jMin,jMax, bi,bj, kLev,       I           iMin,iMax, jMin,jMax, bi,bj, kLev,
763       I           myTime, myThid )       I           myTime, myThid )
764        ENDIF        ENDIF
765  #endif  #endif /* ALLOW_OBCS */
766    
767    #ifdef ALLOW_BBL
768          IF ( useBBL ) CALL BBL_TENDENCY_APPLY_S(
769         &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
770         &                      myTime, myThid )
771    #endif /* ALLOW_BBL */
772    
773    #ifdef ALLOW_MYPACKAGE
774          IF ( useMYPACKAGE ) CALL MYPACKAGE_TENDENCY_APPLY_S(
775         &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
776         &                      myTime, myThid )
777    #endif /* ALLOW_MYPACKAGE */
778    
779        RETURN        RETURN
780        END        END
Legend:
Removed from v.1.35  
changed lines
  Added in v.1.68
  ViewVC Help
Powered by ViewVC 1.1.22