Back to home page

MITgcm
 
 

    

File indexing completed on 2020-02-14 06:10:28 UTC

view on githubraw file Latest commit dff94812 on 2019-12-16 20:16:51 UTC
f15994caab Jean*0001 c#include "PACKAGES_CONFIG.h"
                0002 #include "CPP_OPTIONS.h"
                0003 
                0004 CBOP
                0005 C     !ROUTINE: INI_SIGMA_HFAC
                0006 C     !INTERFACE:
                0007       SUBROUTINE INI_SIGMA_HFAC( myThid )
                0008 C     !DESCRIPTION: \bv
                0009 C     *==========================================================*
                0010 C     | SUBROUTINE INI_SIGMA_HFAC
                0011 C     | o Initialise grid factors when using Sigma coordiante
                0012 C     *==========================================================*
                0013 C     | These arrays are used throughout the code and describe
                0014 C     | fractional height factors.
                0015 C     *==========================================================*
                0016 C     \ev
                0017 
                0018 C     !USES:
                0019       IMPLICIT NONE
                0020 C     === Global variables ===
                0021 #include "SIZE.h"
                0022 #include "EEPARAMS.h"
                0023 #include "PARAMS.h"
                0024 #include "GRID.h"
                0025 #include "SURFACE.h"
                0026 
                0027 C     !INPUT/OUTPUT PARAMETERS:
                0028 C     == Routine arguments ==
                0029 C     myThid  ::  Number of this instance of INI_SIGMA_HFAC
                0030       INTEGER myThid
                0031 
                0032 C     !LOCAL VARIABLES:
                0033 C     == Local variables ==
                0034 C     bi, bj     :: tile indices
                0035 C     i, j, k    :: Loop counters
                0036 C     rEmpty     :: empty column r-position
                0037 C     rFullDepth :: maximum depth of a full column
                0038 C     tmpFld     :: Temporary array used to compute & write Total Depth
                0039 C     min_hFac   :: actual minimum of cell-centered hFac
                0040 C     msgBuf     :: Informational/error message buffer
                0041       INTEGER bi, bj
                0042       INTEGER i, j, k
                0043       _RS rEmpty
                0044       _RL rFullDepth
                0045       _RL tmpFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0046       _RL min_hFac
                0047       _RL hFactmp
                0048       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0049 CEOP
                0050 
                0051 C     r(ij,k,t) = rLow(ij) + aHybSigm(k)*[rF(1)-rF(Nr+1)]
                0052 C               + bHybSigm(k)*[eta(ij,t)+Ro_surf(ij) - rLow(ij)]
                0053 
                0054       IF ( usingPCoords ) rEmpty = rF(Nr+1)
                0055       IF ( usingZCoords ) rEmpty = rF(1)
                0056       rFullDepth = rF(1)-rF(Nr+1)
                0057 
                0058 C---  Calculate partial-cell factor hFacC :
                0059       min_hFac = 1.
                0060       DO bj=myByLo(myThid), myByHi(myThid)
                0061        DO bi=myBxLo(myThid), myBxHi(myThid)
                0062 C-    Remove column (mask=0) thinner than hFacMin*rFullDepth
                0063 C       ensures hFac > hFacMin (assuming we use pure Sigma)
                0064 C Note: because of unfortunate hFacMin default value (=1) (would produce
                0065 C       unexpected empty column), for now, use hFacInf instead of hFacMin
dff94812d5 Jean*0066         DO j=1-OLy,sNy+OLy
                0067          DO i=1-OLx,sNx+OLx
f15994caab Jean*0068            tmpFld(i,j) = Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj)
                0069 c          IF ( tmpFld(i,j).LT.hFacMin*rFullDepth )
                0070            IF ( tmpFld(i,j).LT.hFacInf*rFullDepth )
                0071      &       tmpFld(i,j) = 0. _d 0
                0072          ENDDO
                0073         ENDDO
                0074 c#ifdef ALLOW_SHELFICE
                0075 C--   Would need a specific call here similar to SHELFICE_UPDATE_MASKS
                0076 c     IF ( useShelfIce ) THEN
                0077 c     ENDIF
                0078 c#endif /* ALLOW_SHELFICE */
                0079 C-    Set (or reset) other 2-D cell-centered fields
dff94812d5 Jean*0080         DO j=1-OLy,sNy+OLy
                0081          DO i=1-OLx,sNx+OLx
f15994caab Jean*0082            IF ( tmpFld(i,j).GT.0. _d 0 ) THEN
                0083              kSurfC (i,j,bi,bj) = 1
                0084              kLowC  (i,j,bi,bj) = Nr
                0085              maskInC(i,j,bi,bj) = 1.
                0086              recip_Rcol(i,j,bi,bj) = 1. _d 0 / tmpFld(i,j)
                0087            ELSE
                0088              kSurfC (i,j,bi,bj) = Nr+1
                0089              kLowC  (i,j,bi,bj) = 0
                0090              maskInC(i,j,bi,bj) = 0.
                0091              recip_Rcol(i,j,bi,bj) = 0. _d 0
                0092              Ro_surf(i,j,bi,bj) = rEmpty
                0093              R_low(i,j,bi,bj)   = rEmpty
                0094            ENDIF
                0095          ENDDO
                0096         ENDDO
                0097 C-    Set 3-D hFacC
                0098         DO k=1, Nr
dff94812d5 Jean*0099          DO j=1-OLy,sNy+OLy
                0100           DO i=1-OLx,sNx+OLx
f15994caab Jean*0101            IF ( maskInC(i,j,bi,bj).NE.0. _d 0 ) THEN
                0102              hFactmp = ( dAHybSigF(k)*rFullDepth
                0103      &                 + dBHybSigF(k)*tmpFld(i,j)
                0104      &                 )*recip_drF(k)
                0105              hFacC(i,j,k,bi,bj) = hFactmp
                0106              min_hFac = MIN( min_hFac, hFactmp )
                0107            ELSE
                0108              hFacC(i,j,k,bi,bj) = 0.
                0109            ENDIF
                0110           ENDDO
                0111          ENDDO
                0112         ENDDO
                0113 C-    end bi,bj loops.
                0114        ENDDO
                0115       ENDDO
                0116 
                0117       WRITE(msgBuf,'(A,1PE14.6)')
                0118      &     'S/R INI_SIGMA_HFAC: minimum hFacC=', min_hFac
                0119       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0120      &                    SQUEEZE_RIGHT, myThid )
                0121 
                0122 c     CALL PLOT_FIELD_XYRS(R_low,
                0123 c    &         'Model R_low (ini_masks_etc)', 1, myThid)
                0124 c     CALL PLOT_FIELD_XYRS(Ro_surf,
                0125 c    &         'Model Ro_surf (ini_masks_etc)', 1, myThid)
                0126 
                0127 C--   Set Western & Southern fields (at U and V points)
                0128       DO bj=myByLo(myThid), myByHi(myThid)
                0129        DO bi=myBxLo(myThid), myBxHi(myThid)
                0130 C-    set 2-D mask and rLow & reference rSurf at Western & Southern edges
                0131         i = 1-OlX
dff94812d5 Jean*0132         DO j=1-OLy,sNy+OLy
f15994caab Jean*0133            rSurfW(i,j,bi,bj) = rEmpty
                0134            rLowW (i,j,bi,bj) = rEmpty
                0135            maskInW(i,j,bi,bj)= 0.
                0136         ENDDO
dff94812d5 Jean*0137         j = 1-OLy
                0138         DO i=1-OLx,sNx+OLx
f15994caab Jean*0139            rSurfS(i,j,bi,bj) = rEmpty
                0140            rLowS (i,j,bi,bj) = rEmpty
                0141            maskInS(i,j,bi,bj)= 0.
                0142         ENDDO
dff94812d5 Jean*0143         DO j=1-OLy,sNy+OLy
                0144          DO i=2-OLx,sNx+OLx
f15994caab Jean*0145            maskInW(i,j,bi,bj)= maskInC(i-1,j,bi,bj)*maskInC(i,j,bi,bj)
                0146            rSurfW(i,j,bi,bj) =
                0147      &               ( Ro_surf(i-1,j,bi,bj)
                0148      &               + Ro_surf( i, j,bi,bj) )*0.5 _d 0
                0149            rLowW(i,j,bi,bj)  =
                0150      &                 ( R_low(i-1,j,bi,bj)
                0151      &                 + R_low( i, j,bi,bj) )*0.5 _d 0
                0152 c          rSurfW(i,j,bi,bj) =
                0153 c    &               ( Ro_surf(i-1,j,bi,bj)*rA(i-1,j,bi,bj)
                0154 c    &               + Ro_surf( i, j,bi,bj)*rA( i, j,bi,bj)
                0155 c    &               )*recip_rAw(i,j,bi,bj)*0.5 _d 0
                0156 c          rLowW(i,j,bi,bj)  =
                0157 c    &                 ( R_low(i-1,j,bi,bj)*rA(i-1,j,bi,bj)
                0158 c    &                 + R_low( i, j,bi,bj)*rA( i, j,bi,bj)
                0159 c    &                 )*recip_rAw(i,j,bi,bj)*0.5 _d 0
                0160            IF ( maskInW(i,j,bi,bj).EQ.0. ) THEN
                0161              rSurfW(i,j,bi,bj) = rEmpty
                0162              rLowW (i,j,bi,bj) = rEmpty
                0163            ENDIF
                0164          ENDDO
                0165         ENDDO
dff94812d5 Jean*0166         DO j=2-OLy,sNy+OLy
                0167          DO i=1-OLx,sNx+OLx
f15994caab Jean*0168            maskInS(i,j,bi,bj)= maskInC(i,j-1,bi,bj)*maskInC(i,j,bi,bj)
                0169            rSurfS(i,j,bi,bj) =
                0170      &               ( Ro_surf(i,j-1,bi,bj)
                0171      &               + Ro_surf(i, j, bi,bj) )*0.5 _d 0
                0172            rLowS(i,j,bi,bj)  =
                0173      &                 ( R_low(i,j-1,bi,bj)
                0174      &                 + R_low(i, j, bi,bj) )*0.5 _d 0
                0175 c          rSurfS(i,j,bi,bj) =
                0176 c    &               ( Ro_surf(i,j-1,bi,bj)*rA(i,j-1,bi,bj)
                0177 c    &               + Ro_surf(i, j, bi,bj)*rA(i, j, bi,bj)
                0178 c    &               )*recip_rAs(i,j,bi,bj)*0.5 _d 0
                0179 c          rLowS(i,j,bi,bj)  =
                0180 c    &                 ( R_low(i,j-1,bi,bj)*rA(i,j-1,bi,bj)
                0181 c    &                 + R_low(i, j, bi,bj)*rA(i, j, bi,bj)
                0182 c    &                 )*recip_rAs(i,j,bi,bj)*0.5 _d 0
                0183            IF ( maskInS(i,j,bi,bj).EQ.0. ) THEN
                0184              rSurfS(i,j,bi,bj) = rEmpty
                0185              rLowS (i,j,bi,bj) = rEmpty
                0186            ENDIF
                0187          ENDDO
                0188         ENDDO
                0189        ENDDO
                0190       ENDDO
                0191       CALL EXCH_UV_XY_RS( rSurfW,  rSurfS,  .FALSE., myThid )
                0192       CALL EXCH_UV_XY_RS( rLowW,   rLowS,   .FALSE., myThid )
                0193       CALL EXCH_UV_XY_RS( maskInW, maskInS, .FALSE., myThid )
                0194 
                0195 C-    Set hFacW and hFacS (at U and V points)
                0196       DO bj=myByLo(myThid), myByHi(myThid)
                0197        DO bi=myBxLo(myThid), myBxHi(myThid)
                0198         DO k=1, Nr
dff94812d5 Jean*0199          DO j=1-OLy,sNy+OLy
                0200           DO i=1-OLx,sNx+OLx
f15994caab Jean*0201             hFactmp =
                0202      &          ( dAHybSigF(k)*rFullDepth
                0203      &          + dBHybSigF(k)*( rSurfW(i,j,bi,bj)-rLowW(i,j,bi,bj) )
                0204      &          )*recip_drF(k)
                0205             hFacW(i,j,k,bi,bj) = hFactmp*maskInW(i,j,bi,bj)
                0206           ENDDO
                0207          ENDDO
                0208         ENDDO
                0209         DO k=1, Nr
dff94812d5 Jean*0210          DO j=1-OLy,sNy+OLy
                0211           DO i=1-OLx,sNx+OLx
f15994caab Jean*0212             hFactmp =
                0213      &          ( dAHybSigF(k)*rFullDepth
                0214      &          + dBHybSigF(k)*( rSurfS(i,j,bi,bj)-rLowS(i,j,bi,bj) )
                0215      &          )*recip_drF(k)
                0216             hFacS(i,j,k,bi,bj) = hFactmp*maskInS(i,j,bi,bj)
                0217           ENDDO
                0218          ENDDO
                0219         ENDDO
                0220 C-    Set surface k index for interface W & S (U & V points)
dff94812d5 Jean*0221         DO j=1-OLy,sNy+OLy
                0222          DO i=1-OLx,sNx+OLx
f15994caab Jean*0223            kSurfW(i,j,bi,bj) = Nr+1
                0224            kSurfS(i,j,bi,bj) = Nr+1
                0225            IF ( maskInW(i,j,bi,bj).NE.0. ) kSurfW(i,j,bi,bj) = 1
                0226            IF ( maskInS(i,j,bi,bj).NE.0. ) kSurfS(i,j,bi,bj) = 1
                0227          ENDDO
                0228         ENDDO
                0229 C-    end bi,bj loops.
                0230        ENDDO
                0231       ENDDO
                0232 
dff94812d5 Jean*0233 C--   Additional closing of Western and Southern grid-cell edges: for example,
                0234 C     a) might add some "thin walls" in specific location
                0235 C     b) close non-periodic N & S boundaries of lat-lon grid at the N/S poles.
                0236 C     new: location now reccorded as kSurfW/S = Nr+2
                0237       CALL ADD_WALLS2MASKS( rEmpty, myThid )
                0238 
f15994caab Jean*0239       RETURN
                0240       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!