/[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.44 by dimitri, Mon Jul 23 00:17:44 2007 UTC revision 1.66 by heimbach, Wed May 21 10:44:59 2014 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "PACKAGES_CONFIG.h"  #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    #ifdef ALLOW_SALT_PLUME
7    #include "SALT_PLUME_OPTIONS.h"
8    #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 45  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 65  C--   Forcing term Line 76  C--   Forcing term
76       &                      myTime, myThid )       &                      myTime, myThid )
77  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
78    
79  #ifdef ALLOW_FIZHI  #ifdef ALLOW_ATM_PHYS
80        IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_U(        IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_U(
81       &                      iMin,iMax, jMin,jMax, bi,bj, kLev,       &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
82       &                      myTime, myThid )       &                      myTime, myThid )
83  #endif /* ALLOW_FIZHI */  #endif /* ALLOW_ATM_PHYS */
84    
85  #ifdef ALLOW_MYPACKAGE  #ifdef ALLOW_FIZHI
86        IF ( useMYPACKAGE ) CALL MYPACKAGE_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_MYPACKAGE */  #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 104  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 149  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 169  C--   Forcing term Line 203  C--   Forcing term
203       &                      myTime, myThid )       &                      myTime, myThid )
204  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
205    
206  #ifdef ALLOW_FIZHI  #ifdef ALLOW_ATM_PHYS
207        IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_V(        IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_V(
208       &                      iMin,iMax, jMin,jMax, bi,bj, kLev,       &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
209       &                      myTime, myThid )       &                      myTime, myThid )
210  #endif /* ALLOW_FIZHI */  #endif /* ALLOW_ATM_PHYS */
211    
212  #ifdef ALLOW_MYPACKAGE  #ifdef ALLOW_FIZHI
213        IF ( useMYPACKAGE ) CALL MYPACKAGE_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_MYPACKAGE */  #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 208  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 254  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  CEOP  CEOP
321    c#ifdef ALLOW_FRICTION_HEATING
322    c     _RL tmpFac
323    c#endif
324  #ifdef SHORTWAVE_HEATING  #ifdef SHORTWAVE_HEATING
       integer two  
325        _RL minusone        _RL minusone
326        parameter (two=2,minusone=-1.)        PARAMETER (minusOne=-1.)
327        _RL swfracb(two)        _RL swfracb(2)
328        INTEGER kp1        INTEGER kp1
329  #endif  #endif
330    
331        IF ( fluidIsAir ) THEN        IF ( fluidIsAir ) THEN
332         kSurface = 0         kSurface = 0
333          ELSEIF ( usingZCoords .AND. useShelfIce ) THEN
334           kSurface = -1
335        ELSEIF ( usingPCoords ) THEN        ELSEIF ( usingPCoords ) THEN
336         kSurface = Nr         kSurface = Nr
337        ELSE        ELSE
# Line 281  C--   Forcing term Line 345  C--   Forcing term
345       &                      myTime, myThid )       &                      myTime, myThid )
346  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
347    
348    #ifdef ALLOW_ATM_PHYS
349          IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_T(
350         &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
351         &                      myTime, myThid )
352    #endif /* ALLOW_ATM_PHYS */
353    
354  #ifdef ALLOW_FIZHI  #ifdef ALLOW_FIZHI
355        IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_T(        IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_T(
356       &                      iMin,iMax, jMin,jMax, bi,bj, kLev,       &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
357       &                      myTime, myThid )       &                      myTime, myThid )
358  #endif /* ALLOW_FIZHI */  #endif /* ALLOW_FIZHI */
359    
360  #ifdef ALLOW_MYPACKAGE  #ifdef ALLOW_ADDFLUID
361        IF ( useMYPACKAGE ) CALL MYPACKAGE_TENDENCY_APPLY_T(        IF ( selectAddFluid.NE.0 .AND. temp_addMass.NE.UNSET_RL ) THEN
362       &                      iMin,iMax, jMin,jMax, bi,bj, kLev,         IF ( ( selectAddFluid.GE.1 .AND. nonlinFreeSurf.GT.0 )
363       &                      myTime, myThid )       &      .OR. convertFW2Salt.EQ.-1. _d 0 ) THEN
364  #endif /* ALLOW_MYPACKAGE */           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     Add heat in top-layer  C     Ocean: Add temperature surface forcing (e.g., heat-flux) in surface level
475        IF ( kLev .EQ. kSurface ) THEN        IF ( kLev .EQ. kSurface ) THEN
476         DO j=1,sNy         DO j=1,sNy
477          DO i=1,sNx          DO i=1,sNx
478           gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)            gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
479       &     +surfaceForcingT(i,j,bi,bj)       &      +surfaceForcingT(i,j,bi,bj)
480       &     *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)       &      *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
481            ENDDO
482           ENDDO
483          ELSEIF ( kSurface.EQ.-1 ) THEN
484           DO j=1,sNy
485            DO i=1,sNx
486             IF ( kSurfC(i,j,bi,bj).EQ.kLev ) THEN
487              gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
488         &      +surfaceForcingT(i,j,bi,bj)
489         &      *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
490             ENDIF
491          ENDDO          ENDDO
492         ENDDO         ENDDO
493        ENDIF        ENDIF
494    
 #ifndef ALLOW_AUTODIFF_TAMC  
495        IF (linFSConserveTr) THEN        IF (linFSConserveTr) THEN
496         DO j=1,sNy         DO j=1,sNy
497          DO i=1,sNx          DO i=1,sNx
498            IF (kLev .EQ. ksurfC(i,j,bi,bj)) THEN            IF (kLev .EQ. kSurfC(i,j,bi,bj)) THEN
499              gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)              gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
500       &        +TsurfCor*recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)       &        +TsurfCor*recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
501            ENDIF            ENDIF
502          ENDDO          ENDDO
503         ENDDO         ENDDO
504        ENDIF        ENDIF
505  #endif /* ndfef ALLOW_AUTODIFF_TAMC */  
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  #ifdef ALLOW_SHELFICE
514        IF ( useShelfIce )        IF ( useShelfIce )
# Line 324  C     Add heat in top-layer Line 517  C     Add heat in top-layer
517       I     myTime, myThid )       I     myTime, myThid )
518  #endif /* ALLOW_SHELFICE */  #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  c     IF ( usePenetratingSW ) THEN  c     IF ( usePenetratingSW ) THEN
529         swfracb(1)=abs(rF(klev))         swfracb(1)=abs(rF(klev))
530         swfracb(2)=abs(rF(klev+1))         swfracb(2)=abs(rF(klev+1))
531         CALL SWFRAC(         CALL SWFRAC(
532       I             two, minusone,       I             2, minusOne,
533       U             swfracb,       U             swfracb,
534       I             myTime, 1, myThid )       I             myTime, 1, myThid )
535         kp1 = klev+1         kp1 = klev+1
# Line 343  c     IF ( usePenetratingSW ) THEN Line 542  c     IF ( usePenetratingSW ) THEN
542           gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj)           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)       &   -Qsw(i,j,bi,bj)*(swfracb(1)*maskC(i,j,klev,bi,bj)
544       &                   -swfracb(2)*maskC(i,j,kp1, bi,bj))       &                   -swfracb(2)*maskC(i,j,kp1, bi,bj))
545       &    *recip_Cp*recip_rhoConst       &    *mass2rUnit/HeatCapacity_Cp
546       &    *recip_drF(klev)*_recip_hFacC(i,j,kLev,bi,bj)       &    *recip_drF(klev)*_recip_hFacC(i,j,kLev,bi,bj)
547          ENDDO          ENDDO
548         ENDDO         ENDDO
# Line 351  c     ENDIF Line 550  c     ENDIF
550  #endif  #endif
551    
552  #ifdef ALLOW_RBCS  #ifdef ALLOW_RBCS
553         if (useRBCS) then         IF (useRBCS) THEN
554            call RBCS_ADD_TENDENCY(bi,bj,klev, 1,            CALL RBCS_ADD_TENDENCY(bi,bj,klev, 1,
555       &                            myTime, myThid )       &                            myTime, myThid )
556         endif         ENDIF
557  #endif  #endif
558    
559  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
# Line 365  c     ENDIF Line 564  c     ENDIF
564        ENDIF        ENDIF
565  #endif  #endif
566    
567    #ifdef ALLOW_BBL
568          IF ( useBBL ) CALL BBL_TENDENCY_APPLY_T(
569         &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
570         &                      myTime, myThid )
571    #endif /* ALLOW_BBL */
572    
573    #ifdef ALLOW_MYPACKAGE
574          IF ( useMYPACKAGE ) CALL MYPACKAGE_TENDENCY_APPLY_T(
575         &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
576         &                      myTime, myThid )
577    #endif /* ALLOW_MYPACKAGE */
578    
579        RETURN        RETURN
580        END        END
581    
# Line 396  C     == Global data == Line 607  C     == Global data ==
607  #include "DYNVARS.h"  #include "DYNVARS.h"
608  #include "FFIELDS.h"  #include "FFIELDS.h"
609  #include "SURFACE.h"  #include "SURFACE.h"
 #ifdef ALLOW_SALT_PLUME  
 #ifdef ALLOW_SEAICE  
 #include "SEAICE_PARAMS.h"  
 #endif /* ALLOW_SEAICE */  
 #endif /* ALLOW_SALT_PLUME */  
610    
611  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
612  C     == Routine arguments ==  C     == Routine arguments ==
# Line 417  C     myThid    :: Thread Id number Line 623  C     myThid    :: Thread Id number
623  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
624  C     == Local variables ==  C     == Local variables ==
625  C     i,j       :: Loop counters  C     i,j       :: Loop counters
626  C     kSurface  :: index of surface layer  C     kSurface  :: index of surface level
627        INTEGER i, j        INTEGER i, j
628        INTEGER kSurface        INTEGER kSurface
629  CEOP  CEOP
 #ifdef ALLOW_SALT_PLUME  
       _RL saltPlume  
 #endif /* ALLOW_SALT_PLUME */  
630    
631        IF ( fluidIsAir ) THEN        IF ( fluidIsAir ) THEN
632         kSurface = 0         kSurface = 0
633          ELSEIF ( usingZCoords .AND. useShelfIce ) THEN
634           kSurface = -1
635        ELSEIF ( usingPCoords ) THEN        ELSEIF ( usingPCoords ) THEN
636         kSurface = Nr         kSurface = Nr
637        ELSE        ELSE
# Line 440  C--   Forcing term Line 645  C--   Forcing term
645       &                      myTime, myThid )       &                      myTime, myThid )
646  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
647    
648    #ifdef ALLOW_ATM_PHYS
649          IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_S(
650         &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
651         &                      myTime, myThid )
652    #endif /* ALLOW_ATM_PHYS */
653    
654  #ifdef ALLOW_FIZHI  #ifdef ALLOW_FIZHI
655        IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_S(        IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_S(
656       &                      iMin,iMax, jMin,jMax, bi,bj, kLev,       &                      iMin,iMax, jMin,jMax, bi,bj, kLev,
657       &                      myTime, myThid )       &                      myTime, myThid )
658  #endif /* ALLOW_FIZHI */  #endif /* ALLOW_FIZHI */
659    
660  #ifdef ALLOW_MYPACKAGE  #ifdef ALLOW_ADDFLUID
661        IF ( useMYPACKAGE ) CALL MYPACKAGE_TENDENCY_APPLY_S(        IF ( selectAddFluid.NE.0 .AND. salt_addMass.NE.UNSET_RL ) THEN
662       &                      iMin,iMax, jMin,jMax, bi,bj, kLev,         IF ( ( selectAddFluid.GE.1 .AND. nonlinFreeSurf.GT.0 )
663       &                      myTime, myThid )       &      .OR. convertFW2Salt.EQ.-1. _d 0 ) THEN
664  #endif /* ALLOW_MYPACKAGE */           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     Add fresh-water in top-layer  C     Ocean: Add salinity surface forcing (e.g., fresh-water) in surface level
690        IF ( kLev .EQ. kSurface ) THEN        IF ( kLev .EQ. kSurface ) THEN
691         DO j=1,sNy         DO j=1,sNy
692          DO i=1,sNx          DO i=1,sNx
693           gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)            gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
694       &     +surfaceForcingS(i,j,bi,bj)       &      +surfaceForcingS(i,j,bi,bj)
695       &     *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)       &      *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
696            ENDDO
697           ENDDO
698          ELSEIF ( kSurface.EQ.-1 ) THEN
699           DO j=1,sNy
700            DO i=1,sNx
701             IF ( kSurfC(i,j,bi,bj).EQ.kLev ) THEN
702              gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
703         &      +surfaceForcingS(i,j,bi,bj)
704         &      *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
705             ENDIF
706          ENDDO          ENDDO
707         ENDDO         ENDDO
708        ENDIF        ENDIF
709    
 #ifndef ALLOW_AUTODIFF_TAMC  
710        IF (linFSConserveTr) THEN        IF (linFSConserveTr) THEN
711         DO j=1,sNy         DO j=1,sNy
712          DO i=1,sNx          DO i=1,sNx
713            IF (kLev .EQ. ksurfC(i,j,bi,bj)) THEN            IF (kLev .EQ. kSurfC(i,j,bi,bj)) THEN
714              gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)              gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
715       &        +SsurfCor*recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)       &        +SsurfCor*recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
716            ENDIF            ENDIF
717          ENDDO          ENDDO
718         ENDDO         ENDDO
719        ENDIF        ENDIF
 #endif /* ndfef ALLOW_AUTODIFF_TAMC */  
720    
721  #ifdef ALLOW_SHELFICE  #ifdef ALLOW_SHELFICE
722        IF ( useShelfIce )        IF ( useShelfIce )
# Line 483  C     Add fresh-water in top-layer Line 725  C     Add fresh-water in top-layer
725       I     myTime, myThid )       I     myTime, myThid )
726  #endif /* ALLOW_SHELFICE */  #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  #ifdef ALLOW_SALT_PLUME
737  C saltPlume is the amount of salt rejected by ice while freezing;  CC#ifndef SALT_PLUME_VOLUME
738  C it is here redistributed to multiple vertical levels as per        IF ( useSALT_PLUME )
739  C Duffy et al. (GRL 1999)       &     CALL SALT_PLUME_TENDENCY_APPLY_S(
740         DO j=1,sNy       I     iMin,iMax, jMin,jMax, bi,bj, kLev,
741          DO i=1,sNx       I     myTime, myThid )
742            saltPlume = 0.  #ifdef SALT_PLUME_VOLUME
743  #ifdef ALLOW_SEAICE        IF ( useSALT_PLUME )
744            IF ( saltFlux(i,j,bi,bj) .GT. 0. .AND.       &     CALL SALT_PLUME_TENDENCY_APPLY_T(
745       &         salt(i,j,kSurface,bi,bj)  .GT. SEAICE_salinity ) THEN       I     iMin,iMax, jMin,jMax, bi,bj, kLev,
746             saltPlume = (salt(i,j,kSurface,bi,bj)-SEAICE_salinity) *       I     myTime, myThid )
747       &          saltFlux(i,j,bi,bj) / salt(i,j,kSurface,bi,bj)  #endif /* ndef SALT_PLUME_VOLUME */
           ENDIF  
 #endif /* ALLOW_SEAICE */  
           IF ( SaltPlumeDepth(i,j,bi,bj) .GT. -rF(kLev) ) THEN  
            gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)  
      &          +saltPlume*horiVertRatio*recip_rhoConst  
      &          *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)  
      &          *max(drF(kLev),SaltPlumeDepth(i,j,bi,bj)+rF(kLev))  
      &          /SaltPlumeDepth(i,j,bi,bj)  
           ENDIF  
         ENDDO  
        ENDDO  
748  #endif /* ALLOW_SALT_PLUME */  #endif /* ALLOW_SALT_PLUME */
749    
750  #ifdef ALLOW_RBCS  #ifdef ALLOW_RBCS
751         if (useRBCS) then         IF (useRBCS) THEN
752            call RBCS_ADD_TENDENCY(bi,bj,klev, 2,            CALL RBCS_ADD_TENDENCY(bi,bj,klev, 2,
753       &                            myTime, myThid )       &                            myTime, myThid )
754         endif         ENDIF
755  #endif /* ALLOW_RBCS */  #endif /* ALLOW_RBCS */
756    
757  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
# Line 523  C Duffy et al. (GRL 1999) Line 762  C Duffy et al. (GRL 1999)
762        ENDIF        ENDIF
763  #endif /* ALLOW_OBCS */  #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.44  
changed lines
  Added in v.1.66
  ViewVC Help
Powered by ViewVC 1.1.22