pkg/mom_common/mom_calc_hfacz.F

9c05b3873e Alis*0001 #include "MOM_COMMON_OPTIONS.h"
14bb46b4fa Jean*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
aecc8b0f47 Mart*0005 #ifdef ALLOW_CTRL
                0006 # include "CTRL_OPTIONS.h"
                0007 #endif
aea29c8517 Alis*0008
71207ba845 Alis*0009 CBOP
                0010 C !ROUTINE: MOM_CALC_HFACZ
                0011
                0012 C !INTERFACE: ==========================================================
aea29c8517 Alis*0013       SUBROUTINE MOM_CALC_HFACZ(
3678fae08f Jean*0014      I        bi, bj, k,
                0015      O        hFacZ, r_hFacZ,
                0016      I        myThid )
aea29c8517 Alis*0017
71207ba845 Alis*0018 C !DESCRIPTION:
                0019 C Calculates the fractional thickness at vorticity points
                0020
                0021 C !USES: ===============================================================
                0022       IMPLICIT NONE
aea29c8517 Alis*0023 #include "SIZE.h"
9354528577 Jean*0024 #include "EEPARAMS.h"
                0025 #include "PARAMS.h"
aea29c8517 Alis*0026 #include "GRID.h"
3678fae08f Jean*0027 #ifdef ALLOW_DEPTH_CONTROL
7c50f07931 Mart*0028 # ifdef ALLOW_AUTODIFF_TAMC
3678fae08f Jean*0029 #  include "tamc.h"
                0030 # endif
                0031 #else /* ALLOW_DEPTH_CONTROL */
                0032 # ifdef ALLOW_EXCH2
                0033 #  include "W2_EXCH2_SIZE.h"
                0034 #  include "W2_EXCH2_TOPOLOGY.h"
                0035 # endif
                0036 #endif /* ALLOW_DEPTH_CONTROL */
616600b8d2 Patr*0037
71207ba845 Alis*0038 C !INPUT PARAMETERS: ===================================================
3678fae08f Jean*0039 C  bi,bj      :: tile indices
                0040 C  k          :: vertical level
                0041 C  myThid     :: my Thread Id number
                0042       INTEGER bi, bj, k
71207ba845 Alis*0043       INTEGER myThid
                0044
                0045 C !OUTPUT PARAMETERS: ==================================================
3678fae08f Jean*0046 C  hFacZ      :: fractional thickness at vorticity points
                0047 C  r_hFacZ    :: reciprocal
aea29c8517 Alis*0048       _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0049       _RS r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0050
71207ba845 Alis*0051 C !LOCAL VARIABLES: ====================================================
3678fae08f Jean*0052 C  i,j        :: loop indices
9354528577 Jean*0053 C  hZoption   :: forward mode option to select the way hFacZ is computed:
                0054 C                0 : = minimum of 4 hFacW,hFacS arround (consistent with
                0055 C                    definition of partial cell & mask near topography)
                0056 C                1 : = minimum of 2 average (hFacW)_j,(hFacS)_i
                0057 C                2 : = average of 4 hFacW,hFacS arround (consistent with
                0058 C                    how free surface affects hFacW,hFacS it using r* and
                0059 C                    without topography)
3678fae08f Jean*0060       INTEGER i,j
616600b8d2 Patr*0061 #ifdef ALLOW_DEPTH_CONTROL
351ff83338 Jean*0062       _RS hFacZOpenI(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0063       _RS hFacZOpenJ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
616600b8d2 Patr*0064 # ifdef USE_SMOOTH_MIN
bcc0cf4625 Jean*0065       _RS      SMOOTHMIN_RS
                0066       EXTERNAL SMOOTHMIN_RS
616600b8d2 Patr*0067 # endif /* USE_SMOOTH_MIN */
7c50f07931 Mart*0068 # ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0069       INTEGER kkey
7c50f07931 Mart*0070 # endif
3678fae08f Jean*0071 #else /* ALLOW_DEPTH_CONTROL */
0e80da5578 Jean*0072       _RS     hFacZOpen
9354528577 Jean*0073       INTEGER hZoption
                0074       LOGICAL northWestCorner, northEastCorner,
                0075      &        southWestCorner, southEastCorner
                0076       INTEGER myFace
3678fae08f Jean*0077 # ifdef ALLOW_EXCH2
9354528577 Jean*0078       INTEGER myTile
3678fae08f Jean*0079 # endif /* ALLOW_EXCH2 */
9354528577 Jean*0080       PARAMETER ( hZoption = 0 )
                0081 #endif /* ALLOW_DEPTH_CONTROL */
3678fae08f Jean*0082 CEOP
                0083
                0084 C--   Calculate open water fraction at vorticity points
aea29c8517 Alis*0085
616600b8d2 Patr*0086 #ifdef ALLOW_DEPTH_CONTROL
3678fae08f Jean*0087
                0088 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0089       kkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
                0090       kkey = k + (kkey-1)*Nr
616600b8d2 Patr*0091 #endif /* ALLOW_AUTODIFF_TAMC */
                0092
bcc0cf4625 Jean*0093       DO j=1-OLy,sNy+OLy
                0094        DO i=1-OLx,sNx+OLx
616600b8d2 Patr*0095         hFacZ(i,j)     =0.
                0096         r_hFacZ(i,j)   =0.
3678fae08f Jean*0097         hFacZOpenI(i,j)=0.
616600b8d2 Patr*0098         hFacZOpenJ(i,j)=0.
                0099        ENDDO
                0100       ENDDO
                0101
                0102 #ifdef    ALLOW_AUTODIFF_TAMC
                0103 CADJ STORE      hFacZ(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
                0104 CADJ STORE    r_hFacZ(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
                0105 #endif    /* ALLOW_AUTODIFF_TAMC */
bcc0cf4625 Jean*0106       DO j=2-OLy,sNy+OLy
                0107        DO i=2-OLx,sNx+OLx
616600b8d2 Patr*0108         hFacZOpenJ(i,j)=
                0109 #ifdef USE_SMOOTH_MIN
bcc0cf4625 Jean*0110      &       SMOOTHMIN_RS(_hFacW(i  ,j  ,k,bi,bj),
616600b8d2 Patr*0111 #else
                0112      &                MIN(_hFacW(i  ,j  ,k,bi,bj),
                0113 #endif /* USE_SMOOTH_MIN */
                0114      &                    _hFacW(i  ,j-1,k,bi,bj))
                0115      &       *maskW(i,j,k,bi,bj)*maskW(i,j-1,k,bi,bj)
                0116         hFacZOpenI(i,j)=
                0117 #ifdef USE_SMOOTH_MIN
bcc0cf4625 Jean*0118      &       SMOOTHMIN_RS(_hFacS(i  ,j  ,k,bi,bj),
616600b8d2 Patr*0119 #else
                0120      &                MIN(_hFacS(i  ,j  ,k,bi,bj),
                0121 #endif /* USE_SMOOTH_MIN */
                0122      &                    _hFacS(i-1,j  ,k,bi,bj))
                0123      &         *maskS(i,j,k,bi,bj)*maskS(i-1,j,k,bi,bj)
                0124        ENDDO
                0125       ENDDO
3678fae08f Jean*0126 #ifdef ALLOW_AUTODIFF_TAMC
616600b8d2 Patr*0127 CADJ STORE hFacZOpenI(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
                0128 CADJ STORE hFacZOpenJ(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
3678fae08f Jean*0129 #endif /* ALLOW_AUTODIFF_TAMC */
bcc0cf4625 Jean*0130       DO j=2-OLy,sNy+OLy
                0131        DO i=2-OLx,sNx+OLx
616600b8d2 Patr*0132         hFacZ(i,j) =
                0133 #ifdef USE_SMOOTH_MIN
bcc0cf4625 Jean*0134      &       SMOOTHMIN_RS(hFacZOpenI(i,j),hFacZOpenJ(i,j))
616600b8d2 Patr*0135 #else
                0136      &                MIN(hFacZOpenI(i,j),hFacZOpenJ(i,j))
                0137 #endif /* USE_SMOOTH_MIN */
                0138      &         *maskW(i,j,k,bi,bj)*maskW(i,j-1,k,bi,bj)
                0139      &         *maskS(i,j,k,bi,bj)*maskS(i-1,j,k,bi,bj)
                0140        ENDDO
                0141       ENDDO
3678fae08f Jean*0142 #ifdef ALLOW_AUTODIFF_TAMC
616600b8d2 Patr*0143 CADJ STORE hFacZ(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
3678fae08f Jean*0144 #endif /* ALLOW_AUTODIFF_TAMC */
bcc0cf4625 Jean*0145       DO j=2-OLy,sNy+OLy
                0146        DO i=2-OLx,sNx+OLx
616600b8d2 Patr*0147         IF (hFacZ(i,j).EQ.0.) THEN
3678fae08f Jean*0148          r_hFacZ(i,j) = 0. _d 0
616600b8d2 Patr*0149         ELSE
3678fae08f Jean*0150          r_hFacZ(i,j) = 1. _d 0/hFacZ(i,j)
616600b8d2 Patr*0151         ENDIF
                0152        ENDDO
                0153       ENDDO
3678fae08f Jean*0154 #ifdef ALLOW_AUTODIFF_TAMC
616600b8d2 Patr*0155 CADJ STORE    r_hFacZ(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
3678fae08f Jean*0156 #endif /* ALLOW_AUTODIFF_TAMC */
616600b8d2 Patr*0157
                0158 #else /* not ALLOW_DEPTH_CONTROL */
                0159
9354528577 Jean*0160 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0161
                0162 C-    Initialize hFacZ:
bcc0cf4625 Jean*0163 c     DO j=1-OLy,sNy+OLy
                0164 c      DO i=1-OLx,sNx+OLx
9354528577 Jean*0165 c        hFacZ(i,j)=0.
                0166 c      ENDDO
                0167 c     ENDDO
                0168
                0169 C--   1rst row & column are not computed: fill with zero
bcc0cf4625 Jean*0170       DO i=1-OLx,sNx+OLx
                0171         hFacZ(i,1-OLy)=0.
aea29c8517 Alis*0172       ENDDO
bcc0cf4625 Jean*0173       DO j=2-OLy,sNy+OLy
                0174         hFacZ(1-OLx,j)=0.
9354528577 Jean*0175       ENDDO
                0176
                0177 C--   Calculate open water fraction at vorticity points
                0178
                0179       IF ( hZoption.EQ.2 ) THEN
bcc0cf4625 Jean*0180         DO j=2-OLy,sNy+OLy
                0181          DO i=2-OLx,sNx+OLx
9354528577 Jean*0182 c         hFacZOpen=
                0183 c    &               ( _hFacW(i, j ,k,bi,bj)*rAw(i, j ,bi,bj)
                0184 c    &                +_hFacW(i,j-1,k,bi,bj)*rAw(i,j-1,bi,bj) )
                0185 c    &             + ( _hFacS( i ,j,k,bi,bj)*rAs( i ,j,bi,bj)
                0186 c    &                +_hFacS(i-1,j,k,bi,bj)*rAs(i-1,j,bi,bj) )
                0187 c         hFacZ(i,j) = 0.25 _d 0 * hFacZOpen*recip_rAz(i,j,bi,bj)
                0188           hFacZOpen=
                0189      &               ( _hFacW(i, j ,k,bi,bj)
                0190      &                +_hFacW(i,j-1,k,bi,bj) )
                0191      &             + ( _hFacS( i ,j,k,bi,bj)
                0192      &                +_hFacS(i-1,j,k,bi,bj) )
                0193           hFacZ(i,j) = 0.25 _d 0 * hFacZOpen
                0194          ENDDO
                0195         ENDDO
                0196       ELSEIF ( hZoption.EQ.1 ) THEN
bcc0cf4625 Jean*0197         DO j=2-OLy,sNy+OLy
                0198          DO i=2-OLx,sNx+OLx
9354528577 Jean*0199 c         hFacZOpen=MIN(
                0200 c    &                  _hFacW(i, j ,k,bi,bj)*rAw(i, j ,bi,bj)
                0201 c    &                + _hFacW(i,j-1,k,bi,bj)*rAw(i,j-1,bi,bj)
                0202 c    &                , _hFacS( i ,j,k,bi,bj)*rAs( i ,j,bi,bj)
                0203 c    &                + _hFacS(i-1,j,k,bi,bj)*rAs(i-1,j,bi,bj)
                0204 c    &                 )
                0205 c         hFacZ(i,j) = 0.5 _d 0 * hFacZOpen*recip_rAz(i,j,bi,bj)
                0206           hFacZOpen=MIN(
                0207      &                  _hFacW(i, j ,k,bi,bj)
                0208      &                + _hFacW(i,j-1,k,bi,bj)
                0209      &                , _hFacS( i ,j,k,bi,bj)
                0210      &                + _hFacS(i-1,j,k,bi,bj)
                0211      &                 )
                0212           hFacZ(i,j) = 0.5 _d 0 * hFacZOpen
                0213          ENDDO
                0214         ENDDO
                0215       ELSE
bcc0cf4625 Jean*0216         DO j=2-OLy,sNy+OLy
                0217          DO i=2-OLx,sNx+OLx
9354528577 Jean*0218           hFacZOpen=MIN(_hFacW(i,j,k,bi,bj),
                0219      &                  _hFacW(i,j-1,k,bi,bj))
                0220           hFacZOpen=MIN(_hFacS(i,j,k,bi,bj),hFacZOpen)
                0221           hFacZOpen=MIN(_hFacS(i-1,j,k,bi,bj),hFacZOpen)
                0222           hFacZ(i,j)=hFacZOpen
                0223          ENDDO
                0224         ENDDO
                0225       ENDIF
                0226
3678fae08f Jean*0227 C-----------------------------------------
9354528577 Jean*0228 C     Special stuff for Cubed Sphere
                0229       IF ( useCubedSphereExchange .AND. hZoption.GE.1 ) THEN
                0230
                0231 #ifdef ALLOW_EXCH2
c424ee7cc7 Jean*0232         myTile = W2_myTileList(bi,bj)
9354528577 Jean*0233         myFace = exch2_myFace(myTile)
                0234         southWestCorner = exch2_isWedge(myTile).EQ.1
3678fae08f Jean*0235      &              .AND. exch2_isSedge(myTile).EQ.1
9354528577 Jean*0236         southEastCorner = exch2_isEedge(myTile).EQ.1
                0237      &              .AND. exch2_isSedge(myTile).EQ.1
                0238         northEastCorner = exch2_isEedge(myTile).EQ.1
                0239      &              .AND. exch2_isNedge(myTile).EQ.1
                0240         northWestCorner = exch2_isWedge(myTile).EQ.1
                0241      &              .AND. exch2_isNedge(myTile).EQ.1
                0242 #else
                0243         myFace = bi
                0244         southWestCorner = .TRUE.
                0245         southEastCorner = .TRUE.
                0246         northWestCorner = .TRUE.
                0247         northEastCorner = .TRUE.
                0248 #endif /* ALLOW_EXCH2 */
                0249
                0250         IF ( southWestCorner ) THEN
                0251          i=1
                0252          j=1
                0253          IF ( hZoption.EQ.1 ) THEN
                0254           hFacZOpen=MIN(_hFacW(i,j,k,bi,bj),
                0255      &                  _hFacW(i,j-1,k,bi,bj))
                0256           hFacZOpen=MIN(_hFacS(i,j,k,bi,bj),hFacZOpen)
                0257           hFacZ(i,j)=hFacZOpen
                0258          ELSE
                0259           IF ( MOD(myFace,2).EQ.1 ) THEN
                0260             hFacZOpen=
                0261      &               ( _hFacW(i,j-1,k,bi,bj)
                0262      &                +_hFacS( i ,j,k,bi,bj) )
                0263      &               + _hFacW(i, j ,k,bi,bj)
                0264           ELSE
                0265             hFacZOpen=
                0266      &               ( _hFacW(i, j ,k,bi,bj)
                0267      &                +_hFacW(i,j-1,k,bi,bj) )
                0268      &               + _hFacS( i ,j,k,bi,bj)
                0269           ENDIF
                0270           hFacZ(i,j) = hFacZOpen / 3. _d 0
                0271          ENDIF
                0272         ENDIF
                0273
                0274         IF ( southEastCorner ) THEN
3678fae08f Jean*0275          i=sNx+1
                0276          j=1
9354528577 Jean*0277 C-    to get the same truncation, independent from the face Nb:
                0278          IF ( hZoption.EQ.1 ) THEN
                0279           hFacZOpen=MIN(_hFacW(i,j,k,bi,bj),
                0280      &                  _hFacW(i,j-1,k,bi,bj))
                0281           hFacZOpen=MIN(_hFacS(i-1,j,k,bi,bj),hFacZOpen)
                0282           hFacZ(i,j)=hFacZOpen
                0283          ELSE
                0284           IF ( myFace.EQ.4 ) THEN
                0285             hFacZOpen=
                0286      &               ( _hFacS(i-1,j,k,bi,bj)
                0287      &                +_hFacW(i,j-1,k,bi,bj) )
                0288      &               + _hFacW(i, j ,k,bi,bj)
                0289           ELSEIF ( myFace.EQ.6 ) THEN
                0290             hFacZOpen=
                0291      &               ( _hFacW(i,j-1,k,bi,bj)
                0292      &                +_hFacW(i, j ,k,bi,bj) )
                0293      &               + _hFacS(i-1,j,k,bi,bj)
                0294           ELSE
                0295             hFacZOpen=
                0296      &               ( _hFacW(i, j ,k,bi,bj)
                0297      &                +_hFacS(i-1,j,k,bi,bj) )
                0298      &               + _hFacW(i,j-1,k,bi,bj)
                0299           ENDIF
                0300           hFacZ(i,j) = hFacZOpen / 3. _d 0
                0301          ENDIF
                0302         ENDIF
                0303
                0304         IF ( northWestCorner ) THEN
                0305          i=1
                0306          j=sNy+1
                0307 C-    to get the same truncation, independent from the face Nb:
                0308          IF ( hZoption.EQ.1 ) THEN
                0309           hFacZOpen=MIN(_hFacW(i,j,k,bi,bj),
                0310      &                  _hFacW(i,j-1,k,bi,bj))
                0311           hFacZOpen=MIN(_hFacS(i,j,k,bi,bj),hFacZOpen)
                0312           hFacZ(i,j)=hFacZOpen
                0313          ELSE
                0314           IF ( myFace.EQ.1 ) THEN
                0315             hFacZOpen=
                0316      &               ( _hFacS( i ,j,k,bi,bj)
                0317      &                +_hFacW(i, j ,k,bi,bj) )
                0318      &               + _hFacW(i,j-1,k,bi,bj)
                0319           ELSEIF ( myFace.EQ.5 ) THEN
                0320             hFacZOpen=
                0321      &               ( _hFacW(i, j ,k,bi,bj)
                0322      &                +_hFacW(i,j-1,k,bi,bj) )
                0323      &               + _hFacS( i ,j,k,bi,bj)
                0324           ELSE
                0325             hFacZOpen=
                0326      &               ( _hFacW(i,j-1,k,bi,bj)
                0327      &                +_hFacS( i ,j,k,bi,bj) )
                0328      &               + _hFacW(i, j ,k,bi,bj)
                0329           ENDIF
                0330           hFacZ(i,j) = hFacZOpen / 3. _d 0
                0331          ENDIF
                0332         ENDIF
                0333
                0334         IF ( northEastCorner ) THEN
                0335          i=sNx+1
                0336          j=sNy+1
                0337          IF ( hZoption.EQ.1 ) THEN
                0338           hFacZOpen=MIN(_hFacW(i,j,k,bi,bj),
                0339      &                  _hFacW(i,j-1,k,bi,bj))
                0340           hFacZOpen=MIN(_hFacS(i-1,j,k,bi,bj),hFacZOpen)
                0341           hFacZ(i,j)=hFacZOpen
                0342          ELSE
                0343           IF ( MOD(myFace,2).EQ.1 ) THEN
                0344             hFacZOpen=
                0345      &               ( _hFacW(i,j-1,k,bi,bj)
                0346      &                +_hFacW(i, j ,k,bi,bj) )
                0347      &               + _hFacS(i-1,j,k,bi,bj)
                0348           ELSE
                0349             hFacZOpen=
                0350      &               ( _hFacW(i, j ,k,bi,bj)
                0351      &                +_hFacS(i-1,j,k,bi,bj) )
                0352      &               + _hFacW(i,j-1,k,bi,bj)
                0353           ENDIF
                0354           hFacZ(i,j) = hFacZOpen / 3. _d 0
                0355          ENDIF
                0356         ENDIF
                0357
                0358       ENDIF
3678fae08f Jean*0359 C     Special Cubed Sphere block ends here
                0360 C-----------------------------------------
9354528577 Jean*0361
edb6656069 Mart*0362 C--   Calculate reciprocal:
bcc0cf4625 Jean*0363       DO j=1-OLy,sNy+OLy
                0364        DO i=1-OLx,sNx+OLx
aea29c8517 Alis*0365         IF (hFacZ(i,j).EQ.0.) THEN
9354528577 Jean*0366          r_hFacZ(i,j) = 0.
aea29c8517 Alis*0367         ELSE
9354528577 Jean*0368          r_hFacZ(i,j) = 1. _d 0/hFacZ(i,j)
aea29c8517 Alis*0369         ENDIF
                0370        ENDDO
                0371       ENDDO
9354528577 Jean*0372
                0373 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
616600b8d2 Patr*0374 #endif /* ALLOW_DEPTH_CONTROL */
aea29c8517 Alis*0375
                0376       RETURN
                0377       END