model/src/correction_step.F

aacd3965ce Jean*0001 #include "PACKAGES_CONFIG.h"
1dbaea09ee Chri*0002 #include "CPP_OPTIONS.h"
0bb99fb476 Alis*0003
9366854e02 Chri*0004 CBOP
                0005 C     !ROUTINE: CORRECTION_STEP
                0006 C     !INTERFACE:
0bb99fb476 Alis*0007       SUBROUTINE CORRECTION_STEP( bi, bj, iMin, iMax, jMin, jMax,
6c16fd0209 Jean*0008      I                      phiSurfX, phiSurfY,
                0009      I                      myTime, myIter, myThid )
9366854e02 Chri*0010 C     !DESCRIPTION: \bv
                0011 C     *==========================================================*
cb7fa97db9 Jean*0012 C     | S/R CORRECTION_STEP
                0013 C     | o Corrects the horizontal flow fields with the surface
                0014 C     |   pressure (and Non-Hydrostatic pressure).
9366854e02 Chri*0015 C     *==========================================================*
                0016 C     \ev
fba1188151 Patr*0017
9366854e02 Chri*0018 C     !USES:
                0019       IMPLICIT NONE
fba1188151 Patr*0020 C     == Global variables ==
0bb99fb476 Alis*0021 #include "SIZE.h"
81bc00c2f0 Chri*0022 #include "EEPARAMS.h"
0bb99fb476 Alis*0023 #include "PARAMS.h"
                0024 #include "GRID.h"
d1b81ea0bc Jean*0025 #include "DYNVARS.h"
88830be691 Alis*0026 #ifdef ALLOW_NONHYDROSTATIC
df999eca2c Jean*0027 # include "NH_VARS.h"
                0028 #endif
                0029 #ifdef ALLOW_DIAGNOSTICS
                0030 # include "FFIELDS.h"
88830be691 Alis*0031 #endif
9366854e02 Chri*0032
                0033 C     !INPUT/OUTPUT PARAMETERS:
0bb99fb476 Alis*0034 C     == Routine Arguments ==
16f5093311 Jean*0035 C     bi, bj              :: Tile indices
                0036 C     iMin,iMax,jMin,jMax :: Loop counters range
                0037 C     phiSurfX, phiSurfY  :: Surface Potential gradient
                0038 C     myTime              :: Current time in simulation
6c16fd0209 Jean*0039 C     myIter              :: Current iteration number in simulation
16f5093311 Jean*0040 C     myThid              :: my Thread Id number
                0041       _RL     phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0042       _RL     phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
6c16fd0209 Jean*0043       INTEGER bi, bj
                0044       INTEGER iMin, iMax, jMin, jMax
16f5093311 Jean*0045       _RL     myTime
6c16fd0209 Jean*0046       INTEGER myIter
0bb99fb476 Alis*0047       INTEGER myThid
88ea655fd6 Chri*0048
df999eca2c Jean*0049 C     !FUNCTIONS:
                0050 #ifdef ALLOW_DIAGNOSTICS
                0051       LOGICAL  DIAGNOSTICS_IS_ON
                0052       EXTERNAL DIAGNOSTICS_IS_ON
                0053 #endif
                0054
9366854e02 Chri*0055 C     !LOCAL VARIABLES:
0bb99fb476 Alis*0056 C     == Local variables ==
6c16fd0209 Jean*0057 C     i, j         :: Loop counters
                0058 C     k            :: Level index
16f5093311 Jean*0059 C     psFac, nhFac :: Scaling parameters for supressing gradients
6c16fd0209 Jean*0060 C     gU_dpx       :: implicit part of pressure gradient tendency
                0061 C     gV_dpy       :: implicit part of pressure gradient tendency
0bb99fb476 Alis*0062       INTEGER i,j
6c16fd0209 Jean*0063       INTEGER k
16f5093311 Jean*0064       _RL     psFac, nhFac
6c16fd0209 Jean*0065       _RL     gU_dpx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0066       _RL     gV_dpy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
df999eca2c Jean*0067 #ifdef ALLOW_DIAGNOSTICS
                0068       LOGICAL dPhiDiagIsOn, implDisDiagIsOn
622cf2fb44 Jean*0069 # ifdef ALLOW_SOLVE4_PS_AND_DRAG
df999eca2c Jean*0070       _RL tmpFac
622cf2fb44 Jean*0071 # endif
df999eca2c Jean*0072 #endif /* ALLOW_DIAGNOSTICS */
9366854e02 Chri*0073 CEOP
0bb99fb476 Alis*0074
df999eca2c Jean*0075 #ifdef ALLOW_DIAGNOSTICS
                0076       dPhiDiagIsOn = .FALSE.
                0077       IF ( useDiagnostics )
                0078      &  dPhiDiagIsOn = DIAGNOSTICS_IS_ON( 'Um_dPhiX', myThid )
                0079      &            .OR. DIAGNOSTICS_IS_ON( 'Vm_dPhiY', myThid )
                0080       implDisDiagIsOn = .FALSE.
                0081       IF ( useDiagnostics )
                0082      &  implDisDiagIsOn = DIAGNOSTICS_IS_ON( 'Um_ImplD', myThid )
                0083      &               .OR. DIAGNOSTICS_IS_ON( 'Vm_ImplD', myThid )
                0084 #endif
                0085
6c16fd0209 Jean*0086 C--   Loop over all layers, top to bottom
                0087       DO k=1,Nr
                0088
855d57fc61 Jean*0089 #ifdef ALLOW_SOLVE4_PS_AND_DRAG
6c16fd0209 Jean*0090        IF ( selectImplicitDrag.EQ.2 ) THEN
855d57fc61 Jean*0091
df999eca2c Jean*0092 #ifdef ALLOW_DIAGNOSTICS
                0093 C--    Momentum update: separate contribution from surf.pressure
                0094 C      and contribution from implicit viscosity with bottom drag
                0095         IF ( dPhiDiagIsOn ) THEN
                0096          psFac = pfFacMom*implicSurfPress
                0097      &          *recip_deepFacC(k)*recip_rhoFacC(k)
                0098          DO j=jMin,jMax
                0099           DO i=iMin,iMax
                0100            gU_dpx(i,j) = -psFac*phiSurfX(i,j)*_maskW(i,j,k,bi,bj)
                0101            gV_dpy(i,j) = -psFac*phiSurfY(i,j)*_maskS(i,j,k,bi,bj)
                0102           ENDDO
                0103          ENDDO
                0104          CALL DIAGNOSTICS_FILL( gU_dpx,
                0105      &                         'Um_dPhiX', k, 1, -2, bi, bj, myThid )
                0106          CALL DIAGNOSTICS_FILL( gV_dpy,
                0107      &                         'Vm_dPhiY', k, 1, -2, bi, bj, myThid )
                0108         ENDIF
                0109         IF ( implDisDiagIsOn ) THEN
                0110          psFac = pfFacMom*implicSurfPress
                0111          tmpFac = recip_deepFacC(k)*recip_rhoFacC(k)
                0112          DO j=jMin,jMax
                0113           DO i=iMin,iMax
                0114            gU_dpx(i,j) = -psFac*( dU_psFacX(i,j,k,bi,bj) - tmpFac )
                0115      &                   *phiSurfX(i,j)*_maskW(i,j,k,bi,bj)
                0116            gV_dpy(i,j) = -psFac*( dV_psFacY(i,j,k,bi,bj) - tmpFac )
                0117      &                   *phiSurfY(i,j)*_maskS(i,j,k,bi,bj)
                0118           ENDDO
                0119          ENDDO
                0120          CALL DIAGNOSTICS_FILL( gU_dpx,
                0121      &                         'Um_ImplD', k, 1, -2, bi, bj, myThid )
                0122          CALL DIAGNOSTICS_FILL( gV_dpy,
                0123      &                         'Vm_ImplD', k, 1, -2, bi, bj, myThid )
                0124         ENDIF
                0125 #endif /* ALLOW_DIAGNOSTICS */
                0126
6c16fd0209 Jean*0127 C     On/off scaling parameter
                0128         psFac = pfFacMom*implicSurfPress
855d57fc61 Jean*0129
6c16fd0209 Jean*0130 C     Pressure gradient tendency (zonal mom): Implicit part
                0131         DO j=jMin,jMax
                0132          DO i=iMin,iMax
                0133           gU_dpx(i,j) =
                0134      &         -psFac*dU_psFacX(i,j,k,bi,bj)*phiSurfX(i,j)
                0135 c    &                  *_maskW(i,j,k,bi,bj)
                0136          ENDDO
855d57fc61 Jean*0137         ENDDO
                0138
6c16fd0209 Jean*0139 C     Pressure gradient tendency (merid mom): Implicit part
                0140         DO j=jMin,jMax
                0141          DO i=iMin,iMax
                0142           gV_dpy(i,j) =
                0143      &         -psFac*dV_psFacY(i,j,k,bi,bj)*phiSurfY(i,j)
                0144 c    &                  *_maskS(i,j,k,bi,bj)
                0145          ENDDO
855d57fc61 Jean*0146         ENDDO
                0147
6c16fd0209 Jean*0148        ELSE
855d57fc61 Jean*0149 #endif /* ALLOW_SOLVE4_PS_AND_DRAG */
                0150
4606c28752 Jean*0151 C     On/off scaling parameters (including anelastic & deep-model factors)
6c16fd0209 Jean*0152         psFac = pfFacMom*implicSurfPress
16f5093311 Jean*0153      &         *recip_deepFacC(k)*recip_rhoFacC(k)
df999eca2c Jean*0154
                0155 C     Pressure gradient tendency (zonal mom): Implicit part
6c16fd0209 Jean*0156         IF ( use3Dsolver ) THEN
                0157          nhFac = pfFacMom*implicitNHPress
                0158      &          *recip_deepFacC(k)*recip_rhoFacC(k)
df999eca2c Jean*0159          DO j=jMin,jMax
                0160           DO i=iMin,iMax
                0161            gU_dpx(i,j) = -(
6c16fd0209 Jean*0162      &          psFac*phiSurfX(i,j)
                0163 #ifdef ALLOW_NONHYDROSTATIC
                0164      &        + nhFac*_recip_dxC(i,j,bi,bj)
                0165      &           *(phi_nh(i,j,k,bi,bj)-phi_nh(i-1,j,k,bi,bj))
                0166 #endif
                0167      &                   )*_maskW(i,j,k,bi,bj)
df999eca2c Jean*0168           ENDDO
6c16fd0209 Jean*0169          ENDDO
df999eca2c Jean*0170         ELSE
                0171          DO j=jMin,jMax
                0172           DO i=iMin,iMax
                0173            gU_dpx(i,j) = -psFac*phiSurfX(i,j)*_maskW(i,j,k,bi,bj)
                0174           ENDDO
                0175          ENDDO
                0176         ENDIF
6c16fd0209 Jean*0177
                0178 C     Pressure gradient tendency (merid mom): Implicit part
df999eca2c Jean*0179         IF ( use3Dsolver ) THEN
                0180          DO j=jMin,jMax
                0181           DO i=iMin,iMax
                0182            gV_dpy(i,j) = -(
6c16fd0209 Jean*0183      &          psFac*phiSurfY(i,j)
                0184 #ifdef ALLOW_NONHYDROSTATIC
                0185      &        + nhFac*_recip_dyC(i,j,bi,bj)
                0186      &           *(phi_nh(i,j,k,bi,bj)-phi_nh(i,j-1,k,bi,bj))
                0187 #endif
                0188      &                   )*_maskS(i,j,k,bi,bj)
df999eca2c Jean*0189           ENDDO
6c16fd0209 Jean*0190          ENDDO
df999eca2c Jean*0191         ELSE
                0192          DO j=jMin,jMax
                0193           DO i=iMin,iMax
                0194            gV_dpy(i,j) = -psFac*phiSurfY(i,j)*_maskS(i,j,k,bi,bj)
                0195           ENDDO
                0196          ENDDO
                0197         ENDIF
                0198
                0199 #ifdef ALLOW_DIAGNOSTICS
                0200         IF ( dPhiDiagIsOn ) THEN
                0201          CALL DIAGNOSTICS_FILL( gU_dpx,
                0202      &                         'Um_dPhiX', k, 1, -2, bi, bj, myThid )
                0203          CALL DIAGNOSTICS_FILL( gV_dpy,
                0204      &                         'Vm_dPhiY', k, 1, -2, bi, bj, myThid )
                0205         ENDIF
                0206 #endif /* ALLOW_DIAGNOSTICS */
6c16fd0209 Jean*0207
                0208 #ifdef ALLOW_SOLVE4_PS_AND_DRAG
855d57fc61 Jean*0209        ENDIF
6c16fd0209 Jean*0210 #endif /* ALLOW_SOLVE4_PS_AND_DRAG */
0bb99fb476 Alis*0211
6c16fd0209 Jean*0212 C     Update zonal velocity: add implicit pressure gradient tendency
855d57fc61 Jean*0213        DO j=jMin,jMax
                0214         DO i=iMin,iMax
6c16fd0209 Jean*0215           uVel(i,j,k,bi,bj)=( gU(i,j,k,bi,bj)
                0216      &                      + deltaTMom*gU_dpx(i,j)
                0217      &                      )*_maskW(i,j,k,bi,bj)
aacd3965ce Jean*0218 #ifdef ALLOW_OBCS
6c16fd0209 Jean*0219      &                       *maskInW(i,j,bi,bj)
aacd3965ce Jean*0220 #endif
855d57fc61 Jean*0221         ENDDO
0bb99fb476 Alis*0222        ENDDO
                0223
6c16fd0209 Jean*0224 C     Update merid. velocity: add implicit pressure gradient tendency
855d57fc61 Jean*0225        DO j=jMin,jMax
                0226         DO i=iMin,iMax
6c16fd0209 Jean*0227           vVel(i,j,k,bi,bj)=( gV(i,j,k,bi,bj)
                0228      &                      + deltaTMom*gV_dpy(i,j)
                0229      &                      )*_maskS(i,j,k,bi,bj)
aacd3965ce Jean*0230 #ifdef ALLOW_OBCS
6c16fd0209 Jean*0231      &                       *maskInS(i,j,bi,bj)
aacd3965ce Jean*0232 #endif
855d57fc61 Jean*0233         ENDDO
0bb99fb476 Alis*0234        ENDDO
855d57fc61 Jean*0235
6c16fd0209 Jean*0236 C-    end of k loop
                0237       ENDDO
0bb99fb476 Alis*0238
df999eca2c Jean*0239 #ifdef ALLOW_DIAGNOSTICS
                0240       IF ( useDiagnostics .AND. selectImplicitDrag.EQ.2 ) THEN
                0241        IF ( DIAGNOSTICS_IS_ON( 'botTauX ', myThid ) ) THEN
                0242          IF ( usingZCoords ) THEN
                0243 C         kLowC = 0 if dry column
                0244           DO j=jMin,jMax
                0245            DO i=iMin,iMax
                0246             k = MAX( 1, MIN( kLowC(i-1,j,bi,bj), kLowC(i,j,bi,bj) ) )
                0247             botDragU(i,j,bi,bj) = -botDragU(i,j,bi,bj)
                0248      &                          * uVel(i,j,k,bi,bj)
                0249            ENDDO
                0250           ENDDO
                0251          ELSE
                0252 C         kSurfW = Nr+1 if dry column
                0253           DO j=jMin,jMax
                0254            DO i=iMin,iMax
                0255             k = MIN( Nr, kSurfW(i,j,bi,bj) )
                0256             botDragU(i,j,bi,bj) = -botDragU(i,j,bi,bj)
                0257      &                          * uVel(i,j,k,bi,bj)
                0258            ENDDO
                0259           ENDDO
                0260          ENDIF
a54ff488da Jean*0261          CALL DIAGNOSTICS_FILL_RS( botDragU, 'botTauX ',
                0262      &                             0, 1, 1, bi, bj, myThid )
df999eca2c Jean*0263        ENDIF
                0264        IF ( DIAGNOSTICS_IS_ON( 'botTauY ', myThid ) ) THEN
                0265          IF ( usingZCoords ) THEN
                0266 C         kLowC = 0 if dry column
                0267           DO j=jMin,jMax
                0268            DO i=iMin,iMax
                0269             k = MAX( 1, MIN( kLowC(i,j-1,bi,bj), kLowC(i,j,bi,bj) ) )
                0270             botDragV(i,j,bi,bj) = -botDragV(i,j,bi,bj)
                0271      &                          * vVel(i,j,k,bi,bj)
                0272            ENDDO
                0273           ENDDO
                0274          ELSE
                0275 C         kSurfS = Nr+1 if dry column
                0276           DO j=jMin,jMax
                0277            DO i=iMin,iMax
                0278             k = MIN( Nr, kSurfS(i,j,bi,bj) )
                0279             botDragV(i,j,bi,bj) = -botDragV(i,j,bi,bj)
                0280      &                          * vVel(i,j,k,bi,bj)
                0281            ENDDO
                0282           ENDDO
                0283          ENDIF
a54ff488da Jean*0284          CALL DIAGNOSTICS_FILL_RS( botDragV, 'botTauY ',
                0285      &                             0, 1, 1, bi, bj, myThid )
df999eca2c Jean*0286        ENDIF
                0287       ENDIF
1fa1235c63 Jean*0288
                0289 # ifdef ALLOW_SHELFICE
                0290       IF ( useShelfIce .AND. useDiagnostics ) THEN
                0291         CALL SHELFICE_DIAGNOSTICS_DRAG(
                0292      I                uVel(1-OLx,1-OLy,1,bi,bj),
                0293      I                vVel(1-OLx,1-OLy,1,bi,bj),
                0294      I                bi, bj, myIter, myThid )
                0295       ENDIF
                0296 # endif /* ALLOW_SHELFICE */
df999eca2c Jean*0297 #endif /* ALLOW_DIAGNOSTICS */
                0298
0bb99fb476 Alis*0299       RETURN
                0300       END