/[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.61 by jmc, Fri Nov 9 22:49:17 2012 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_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
# Line 156  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 176  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
# Line 268  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  #ifdef ALLOW_FRICTION_HEATING  c#ifdef ALLOW_FRICTION_HEATING
322        _RL tmpFac  c     _RL tmpFac
323  #endif  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 298  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,
# Line 342  C         conversion from W/m^2/r_unit t Line 395  C         conversion from W/m^2/r_unit t
395            tmpFac = (tmpFac/atm_Cp) * mass2rUnit            tmpFac = (tmpFac/atm_Cp) * mass2rUnit
396          ELSE          ELSE
397  C         conversion from W/m^2/r_unit to K/s  C         conversion from W/m^2/r_unit to K/s
398            tmpFac = recip_Cp * mass2rUnit            tmpFac = mass2rUnit/HeatCapacity_Cp
399          ENDIF          ENDIF
400          DO j=1,sNy          DO j=1,sNy
401            DO i=1,sNx            DO i=1,sNx
# Line 355  C         conversion from W/m^2/r_unit t Line 408  C         conversion from W/m^2/r_unit t
408        ENDIF        ENDIF
409  #endif /* ALLOW_FRICTION_HEATING */  #endif /* ALLOW_FRICTION_HEATING */
410    
411  C     Add heat in top-layer        IF ( fluidIsAir .AND. atm_Rq.NE.zeroRL .AND. Nr.NE.1 ) THEN
412    C--   Compressible fluid: account for difference between moist and dry air
413    C     specific volume in Enthalpy equation (+ V.dP term), since only the
414    C     dry air part is accounted for in the (dry) Pot.Temp formulation.
415    C     Used centered averaging from interface to center (consistent with
416    C     conversion term in KE eq) and same discretisation ( [T*Q]_bar_k )
417    C     as for Theta_v in CALC_PHI_HYD
418    
419    C     conversion from in-situ Temp to Pot.Temp
420            tmpFac = (atm_Po/rC(kLev))**atm_kappa
421    C     conversion from W/kg to K/s
422            tmpFac = tmpFac/atm_Cp
423            km = kLev-1
424            kc = kLev
425            kp = kLev+1
426            IF ( kLev.EQ.1 ) THEN
427              DO j=1,sNy
428               DO i=1,sNx
429                tmpVar(i,j) = 0.
430               ENDDO
431              ENDDO
432            ELSE
433              delPI = atm_Cp*( (rC(km)/atm_Po)**atm_kappa
434         &                   - (rC(kc)/atm_Po)**atm_kappa )
435              DO j=1,sNy
436               DO i=1,sNx
437                tmpVar(i,j) = wVel(i,j,kc,bi,bj)*delPI*atm_Rq
438         &                  *( theta(i,j,km,bi,bj)*salt(i,j,km,bi,bj)
439         &                   + theta(i,j,kc,bi,bj)*salt(i,j,kc,bi,bj)
440         &                   )*maskC(i,j,km,bi,bj)*0.25 _d 0
441               ENDDO
442              ENDDO
443            ENDIF
444            IF ( kLev.LT.Nr ) THEN
445              delPI = atm_Cp*( (rC(kc)/atm_Po)**atm_kappa
446         &                   - (rC(kp)/atm_Po)**atm_kappa )
447              DO j=1,sNy
448               DO i=1,sNx
449                tmpVar(i,j) = tmpVar(i,j)
450         &                  + wVel(i,j,kp,bi,bj)*delPI*atm_Rq
451         &                  *( theta(i,j,kc,bi,bj)*salt(i,j,kc,bi,bj)
452         &                   + theta(i,j,kp,bi,bj)*salt(i,j,kp,bi,bj)
453         &                   )*maskC(i,j,kp,bi,bj)*0.25 _d 0
454               ENDDO
455              ENDDO
456            ENDIF
457            DO j=1,sNy
458              DO i=1,sNx
459                gT(i,j,kLev,bi,bj) = gT(i,j,kLev,bi,bj)
460         &         + tmpVar(i,j)*tmpFac
461         &          *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
462              ENDDO
463            ENDDO
464    #ifdef ALLOW_DIAGNOSTICS
465            IF ( useDiagnostics ) THEN
466    C     conversion to W/m^2
467              tmpFac = rUnit2mass
468              CALL DIAGNOSTICS_SCALE_FILL( tmpVar, tmpFac, 1,
469         &                     'MoistCor', kc, 1, 3, bi,bj,myThid )
470            ENDIF
471    #endif /* ALLOW_DIAGNOSTICS */
472          ENDIF
473    
474    C     Ocean: Add temperature surface forcing (e.g., heat-flux) in surface level
475        IF ( kLev .EQ. kSurface ) THEN        IF ( kLev .EQ. kSurface ) THEN
476         DO j=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
# Line 403  c     IF ( usePenetratingSW ) THEN Line 529  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 416  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*mass2rUnit       &    *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 497  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
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 517  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,
# Line 552  C    &          *recip_deepFac2C(kLev)*r Line 686  C    &          *recip_deepFac2C(kLev)*r
686        ENDIF        ENDIF
687  #endif /* ALLOW_ADDFLUID */  #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
# Line 587  C     Add fresh-water in top-layer Line 731  C     Add fresh-water in top-layer
731       &     bi,bj, kLev, myTime, myThid )       &     bi,bj, kLev, myTime, myThid )
732  #endif /* ALLOW_ICEFRONT */  #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    CC#ifndef SALT_PLUME_VOLUME
738        IF ( useSALT_PLUME )        IF ( useSALT_PLUME )
739       &     CALL SALT_PLUME_TENDENCY_APPLY_S(       &     CALL SALT_PLUME_TENDENCY_APPLY_S(
740       I     iMin,iMax, jMin,jMax, bi,bj, kLev,       I     iMin,iMax, jMin,jMax, bi,bj, kLev,
741       I     myTime, myThid )       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 */  #endif /* ALLOW_SALT_PLUME */
749    
750  #ifdef ALLOW_RBCS  #ifdef ALLOW_RBCS
Legend:
Removed from v.1.61  
changed lines
  Added in v.1.66
  ViewVC Help
Powered by ViewVC 1.1.22