Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-02 18:36:29 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
6d57af55dc Jean*0001 #include "PACKAGES_CONFIG.h"
00b29feb62 Jean*0002 #include "CPP_OPTIONS.h"
f450521a99 Jean*0003 #ifdef ALLOW_EXCH2
                0004 # include "W2_OPTIONS.h"
                0005 #endif
00b29feb62 Jean*0006 
                0007 CBOP
                0008 C     !ROUTINE: CALC_R_STAR
                0009 C     !INTERFACE:
                0010       SUBROUTINE CALC_R_STAR( etaFld,
                0011      I                        myTime, myIter, myThid )
                0012 C     !DESCRIPTION: \bv
                0013 C     *==========================================================*
9dbedd0b0f Jean*0014 C     | SUBROUTINE CALC_R_STAR
00b29feb62 Jean*0015 C     | o Calculate new column thickness & scaling factor for r*
455a32c08d Jean*0016 C     |   according to the surface r-position (Non-Lin Free-Surf)
00b29feb62 Jean*0017 C     *==========================================================*
                0018 C     \ev
                0019 
                0020 C     !USES:
                0021       IMPLICIT NONE
                0022 C     == Global variables
                0023 #include "SIZE.h"
                0024 #include "EEPARAMS.h"
                0025 #include "PARAMS.h"
                0026 #include "GRID.h"
                0027 #include "SURFACE.h"
                0028 
                0029 C     !INPUT/OUTPUT PARAMETERS:
                0030 C     == Routine arguments ==
2c5c5e9c4a Jean*0031 C     etaFld    :: current eta field used to update the hFactor
                0032 C     myTime    :: current time in simulation
                0033 C     myIter    :: current iteration number in simulation
                0034 C     myThid    :: thread number for this instance of the routine.
ccbf39d1a4 Jean*0035       _RL etaFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
00b29feb62 Jean*0036       _RL myTime
                0037       INTEGER myIter
                0038       INTEGER myThid
                0039 
                0040 #ifdef NONLIN_FRSURF
                0041 
                0042 C     !LOCAL VARIABLES:
                0043 C     Local variables
7d2d5d27ed Jean*0044 C     rStarAreaWeight :: use area weighted average for rStarFac at U & V point
                0045 C     i,j, bi,bj      :: loop counter
                0046 C     numbWrite       :: count the Number of warning written on STD-ERR file
                0047 C     numbWrMax       ::  maximum  Number of warning written on STD-ERR file
                0048       LOGICAL rStarAreaWeight
2dd93b11e6 An T*0049       _RL maxhFacC
f04f2001af Jean*0050       INTEGER i,j,bi,bj
f450521a99 Jean*0051       INTEGER numbWrite, numbWrMax
f0b80b4d0a Mart*0052       INTEGER icntc1, icntc2, icntw, icnts
00b29feb62 Jean*0053       _RL tmpfldW, tmpfldS
                0054 CEOP
7d2d5d27ed Jean*0055 
7418e6b1e6 Jean*0056 #ifdef W2_FILL_NULL_REGIONS
                0057       INTEGER ii,jj
                0058 #endif
ced1fc8e72 Jean*0059       DATA numbWrite / 0 /
                0060       numbWrMax = Nx*Ny
00b29feb62 Jean*0061 
2dd93b11e6 An T*0062       maxhFacC    = 0. _d 0
                0063 
455a32c08d Jean*0064       rStarAreaWeight = .TRUE.
                0065 C-    Area-weighted average consistent with KE (& vert. advection):
                0066       IF ( vectorInvariantMomentum .AND.
                0067      &     (selectKEscheme.EQ.1 .OR. selectKEscheme.EQ.3)
                0068      &   ) rStarAreaWeight =.FALSE.
                0069 
00b29feb62 Jean*0070 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0071 
6d57af55dc Jean*0072 #ifdef ALLOW_DEBUG
                0073       IF (debugMode) CALL DEBUG_ENTER('CALC_R_STAR',myThid)
                0074 #endif
                0075 
00b29feb62 Jean*0076       DO bj=myByLo(myThid), myByHi(myThid)
                0077        DO bi=myBxLo(myThid), myBxHi(myThid)
455a32c08d Jean*0078 
72a058b866 Gael*0079 C--   before updating rStarFacC/S/W save current fields
ccbf39d1a4 Jean*0080          DO j=1-OLy,sNy+OLy
                0081            DO i=1-OLx,sNx+OLx
72a058b866 Gael*0082              rStarFacNm1C(i,j,bi,bj) = rStarFacC(i,j,bi,bj)
                0083              rStarFacNm1S(i,j,bi,bj) = rStarFacS(i,j,bi,bj)
                0084              rStarFacNm1W(i,j,bi,bj) = rStarFacW(i,j,bi,bj)
                0085            ENDDO
                0086          ENDDO
455a32c08d Jean*0087 
00b29feb62 Jean*0088 C-    1rst bi,bj loop :
                0089 
                0090 C-- copy rStarFacX -> rStarExpX
ccbf39d1a4 Jean*0091         DO j=1-OLy,sNy+OLy
                0092           DO i=1-OLx,sNx+OLx
00b29feb62 Jean*0093             rStarExpC(i,j,bi,bj) = rStarFacC(i,j,bi,bj)
                0094             rStarExpW(i,j,bi,bj) = rStarFacW(i,j,bi,bj)
                0095             rStarExpS(i,j,bi,bj) = rStarFacS(i,j,bi,bj)
                0096           ENDDO
f5ea98246f Jean*0097         ENDDO
00b29feb62 Jean*0098 
                0099 C-- Compute the new column thikness :
                0100         DO j=0,sNy+1
                0101          DO i=0,sNx+1
04a977827b Jean*0102           IF (kSurfC(i,j,bi,bj).LE.Nr ) THEN
9dbedd0b0f Jean*0103            rStarFacC(i,j,bi,bj) =
00b29feb62 Jean*0104      &      (etaFld(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj))
                0105      &      *recip_Rcol(i,j,bi,bj)
                0106           ELSE
                0107            rStarFacC(i,j,bi,bj) = 1.
                0108           ENDIF
                0109          ENDDO
                0110         ENDDO
7d2d5d27ed Jean*0111        IF ( rStarAreaWeight ) THEN
                0112 C-     Area weighted average
ced1fc8e72 Jean*0113         DO j=1,sNy
00b29feb62 Jean*0114          DO i=1,sNx+1
04a977827b Jean*0115           IF ( kSurfW(i,j,bi,bj).LE.Nr ) THEN
4e1e8e48f9 Jean*0116            tmpfldW = rSurfW(i,j,bi,bj) - rLowW(i,j,bi,bj)
9dbedd0b0f Jean*0117            rStarFacW(i,j,bi,bj) =
00b29feb62 Jean*0118      &       ( 0.5 _d 0 *( etaFld(i-1,j,bi,bj)*rA(i-1,j,bi,bj)
                0119      &                    +etaFld(i,j,bi,bj)*rA(i,j,bi,bj)
                0120      &                   )*recip_rAw(i,j,bi,bj)
                0121      &        +tmpfldW )/tmpfldW
                0122           ELSE
                0123            rStarFacW(i,j,bi,bj) = 1.
                0124           ENDIF
ced1fc8e72 Jean*0125          ENDDO
                0126         ENDDO
                0127         DO j=1,sNy+1
                0128          DO i=1,sNx
04a977827b Jean*0129           IF ( kSurfS(i,j,bi,bj).LE.Nr ) THEN
4e1e8e48f9 Jean*0130            tmpfldS = rSurfS(i,j,bi,bj) - rLowS(i,j,bi,bj)
9dbedd0b0f Jean*0131            rStarFacS(i,j,bi,bj) =
00b29feb62 Jean*0132      &       ( 0.5 _d 0 *( etaFld(i,j-1,bi,bj)*rA(i,j-1,bi,bj)
                0133      &                    +etaFld(i,j,bi,bj)*rA(i,j,bi,bj)
                0134      &                   )*recip_rAs(i,j,bi,bj)
                0135      &        +tmpfldS )/tmpfldS
                0136           ELSE
                0137            rStarFacS(i,j,bi,bj) = 1.
                0138           ENDIF
                0139          ENDDO
                0140         ENDDO
7d2d5d27ed Jean*0141        ELSE
                0142 C-     Simple average
                0143         DO j=1,sNy
                0144          DO i=1,sNx+1
                0145           IF ( kSurfW(i,j,bi,bj).LE.Nr ) THEN
                0146            tmpfldW = rSurfW(i,j,bi,bj) - rLowW(i,j,bi,bj)
                0147            rStarFacW(i,j,bi,bj) =
                0148      &       ( 0.5 _d 0 *( etaFld(i-1,j,bi,bj) + etaFld(i,j,bi,bj) )
                0149      &        +tmpfldW )/tmpfldW
                0150           ELSE
                0151            rStarFacW(i,j,bi,bj) = 1.
                0152           ENDIF
                0153          ENDDO
                0154         ENDDO
                0155         DO j=1,sNy+1
                0156          DO i=1,sNx
                0157           IF ( kSurfS(i,j,bi,bj).LE.Nr ) THEN
                0158            tmpfldS = rSurfS(i,j,bi,bj) - rLowS(i,j,bi,bj)
                0159            rStarFacS(i,j,bi,bj) =
                0160      &       ( 0.5 _d 0 *( etaFld(i,j-1,bi,bj) + etaFld(i,j,bi,bj) )
                0161      &        +tmpfldS )/tmpfldS
                0162           ELSE
                0163            rStarFacS(i,j,bi,bj) = 1.
                0164           ENDIF
                0165          ENDDO
                0166         ENDDO
                0167        ENDIF
2c5c5e9c4a Jean*0168 #ifdef ALLOW_OBCS
                0169        IF (useOBCS) THEN
                0170          CALL OBCS_APPLY_R_STAR(
                0171      I                    bi, bj, etaFld,
                0172      U                    rStarFacC, rStarFacW, rStarFacS,
                0173      I                    myTime, myIter, myThid )
                0174        ENDIF
                0175 #endif /* ALLOW_OBCS */
00b29feb62 Jean*0176 
ced1fc8e72 Jean*0177 C-    Needs to do something when r* ratio is too small ;
                0178 C     for now, just stop
f0b80b4d0a Mart*0179         icntc1 = 0
                0180         icntc2 = 0
                0181         icntw  = 0
                0182         icnts  = 0
ced1fc8e72 Jean*0183         DO j=1,sNy+1
                0184          DO i=1,sNx+1
9dbedd0b0f Jean*0185           IF ( rStarFacC(i,j,bi,bj).LT.hFacInf ) THEN
f0b80b4d0a Mart*0186             icntc1 = icntc1 + 1
ced1fc8e72 Jean*0187           ENDIF
9dbedd0b0f Jean*0188           IF ( rStarFacW(i,j,bi,bj).LT.hFacInf ) THEN
f0b80b4d0a Mart*0189             icntw = icntw + 1
ced1fc8e72 Jean*0190           ENDIF
9dbedd0b0f Jean*0191           IF ( rStarFacS(i,j,bi,bj).LT.hFacInf ) THEN
f0b80b4d0a Mart*0192             icnts = icnts + 1
ced1fc8e72 Jean*0193           ENDIF
f0b80b4d0a Mart*0194           IF ( rStarFacC(i,j,bi,bj).GT.hFacSup ) THEN
                0195             icntc2 = icntc2 + 1
2dd93b11e6 An T*0196             maxhFacC = max(rStarFacC(i,j,bi,bj),maxhFacC)
ced1fc8e72 Jean*0197           ENDIF
                0198          ENDDO
                0199         ENDDO
2dd93b11e6 An T*0200 
9dbedd0b0f Jean*0201         IF ( icntc1+icnts+icntw .GT. 0 ) THEN
                0202 C-    Print an error msg and then stop:
2dd93b11e6 An T*0203          DO j=1,sNy+1
                0204           DO i=1,sNx+1
                0205            IF ( rStarFacC(i,j,bi,bj).LT.hFacInf ) THEN
                0206             WRITE(errorMessageUnit,'(A,2I4,A,1F10.6,1P2E14.6)')
                0207      &       ' fail at i,j=',i,j,' ; rStarFacC,H,eta =',
                0208      &       rStarFacC(i,j,bi,bj),
                0209      &       Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj),
                0210      &       etaFld(i,j,bi,bj)
                0211            ENDIF
                0212            IF ( rStarFacW(i,j,bi,bj).LT.hFacInf ) THEN
                0213             tmpfldW = rSurfW(i,j,bi,bj) - rLowW(i,j,bi,bj)
                0214             WRITE(errorMessageUnit,'(A,2I4,A,1F10.6,1P3E14.6)')
                0215      &       ' fail at i,j=',i,j,' ; rStarFacW,H,eta =',
                0216      &        rStarFacW(i,j,bi,bj), tmpfldW,
                0217      &        etaFld(i-1,j,bi,bj), etaFld(i,j,bi,bj)
                0218            ENDIF
                0219            IF ( rStarFacS(i,j,bi,bj).LT.hFacInf ) THEN
                0220             tmpfldS = rSurfS(i,j,bi,bj) - rLowS(i,j,bi,bj)
                0221             WRITE(errorMessageUnit,'(A,2I4,A,1F10.6,1P3E14.6)')
                0222      &       ' fail at i,j=',i,j,' ; rStarFacS,H,eta =',
                0223      &        rStarFacS(i,j,bi,bj), tmpfldS,
                0224      &        etaFld(i,j-1,bi,bj), etaFld(i,j,bi,bj)
                0225            ENDIF
                0226           ENDDO
                0227          ENDDO
                0228          IF ( icntc1  .GT. 0 ) 
                0229      &    WRITE(errorMessageUnit,'(A,I8,A,3I4,I10)')
                0230      &     'WARNING: r*FacC < hFacInf at',icntc1,
                0231      &     ' pts : bi,bj,Thid,Iter=',bi,bj,myThid,myIter
                0232          IF ( icntw  .GT. 0 ) 
                0233      &    WRITE(errorMessageUnit,'(A,I8,A,3I4,I10)')
9dbedd0b0f Jean*0234      &     'WARNING: r*FacW < hFacInf at',icntw,
                0235      &     ' pts : bi,bj,Thid,Iter=',bi,bj,myThid,myIter
2dd93b11e6 An T*0236          IF ( icnts  .GT. 0 ) 
                0237      &    WRITE(errorMessageUnit,'(A,I8,A,3I4,I10)')
9dbedd0b0f Jean*0238      &     'WARNING: r*FacS < hFacInf at',icnts,
                0239      &     ' pts : bi,bj,Thid,Iter=',bi,bj,myThid,myIter
2dd93b11e6 An T*0240          WRITE(errorMessageUnit,'(A)')
                0241      &    'STOP in CALC_R_STAR : too SMALL rStarFac[C,W,S] !'
f0b80b4d0a Mart*0242          STOP 'ABNORMAL END: S/R CALC_R_STAR'
9dbedd0b0f Jean*0243         ENDIF
2dd93b11e6 An T*0244 
9dbedd0b0f Jean*0245 C-- Usefull warning when r*Fac becomes very large:
                0246         IF ( icntc2.GT.0 .AND. numbWrite.LE.numbWrMax ) THEN
                0247          numbWrite = numbWrite + 1
                0248          WRITE(errorMessageUnit,'(A,I8,A,3I4,I10)')
                0249      &    'WARNING: r*FacC > hFacSup at',icntc2,
                0250      &    ' pts : bi,bj,Thid,Iter=',bi,bj,myThid,myIter
2dd93b11e6 An T*0251          WRITE(errorMessageUnit,'(A,E14.6)')
                0252      &    'WARNING: max(hFacC) is ',maxhFacC
9dbedd0b0f Jean*0253         ENDIF
ced1fc8e72 Jean*0254 
00b29feb62 Jean*0255 C-    end 1rst bi,bj loop.
                0256        ENDDO
615c650f5e Jean*0257       ENDDO
00b29feb62 Jean*0258 
9dbedd0b0f Jean*0259        _EXCH_XY_RL( rStarFacC, myThid )
00b29feb62 Jean*0260       CALL EXCH_UV_XY_RL(rStarFacW,rStarFacS,.FALSE.,myThid)
                0261 
                0262 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0263 
                0264       DO bj=myByLo(myThid), myByHi(myThid)
                0265        DO bi=myBxLo(myThid), myBxHi(myThid)
                0266 C-    2nd bi,bj loop :
                0267 
6d57af55dc Jean*0268 #ifdef ALLOW_EXCH2
f450521a99 Jean*0269 #ifdef W2_FILL_NULL_REGIONS
6d57af55dc Jean*0270 C- Note: rStarFacC was non-zero EVERYWHERE before exch, but exch2 put zeros
ccbf39d1a4 Jean*0271 C        in the corner regions of the tile (e.g.:[1-OLx:0,1-OLy:0])
6d57af55dc Jean*0272 C       => need to add those lines (or to fix exch2):
ccbf39d1a4 Jean*0273         DO j=1,OLy
                0274          DO i=1,OLx
6d57af55dc Jean*0275           ii = sNx+i
                0276           jj = sNy+j
                0277 
04a977827b Jean*0278           IF (kSurfC(1-i,1-j,bi,bj).GT.Nr) rStarFacC(1-i,1-j,bi,bj)= 1.
                0279           IF (kSurfC(ii, 1-j,bi,bj).GT.Nr) rStarFacC(ii, 1-j,bi,bj)= 1.
                0280           IF (kSurfC(1-i,jj, bi,bj).GT.Nr) rStarFacC(1-i,jj, bi,bj)= 1.
                0281           IF (kSurfC(ii, jj, bi,bj).GT.Nr) rStarFacC(ii, jj, bi,bj)= 1.
f450521a99 Jean*0282 
04a977827b Jean*0283           IF (kSurfW(1-i,1-j,bi,bj).GT.Nr) rStarFacW(1-i,1-j,bi,bj)= 1.
                0284           IF (kSurfW(ii, 1-j,bi,bj).GT.Nr) rStarFacW(ii, 1-j,bi,bj)= 1.
                0285           IF (kSurfW(1-i,jj, bi,bj).GT.Nr) rStarFacW(1-i,jj, bi,bj)= 1.
                0286           IF (kSurfW(ii, jj, bi,bj).GT.Nr) rStarFacW(ii, jj, bi,bj)= 1.
f450521a99 Jean*0287 
04a977827b Jean*0288           IF (kSurfS(1-i,1-j,bi,bj).GT.Nr) rStarFacS(1-i,1-j,bi,bj)= 1.
                0289           IF (kSurfS(ii, 1-j,bi,bj).GT.Nr) rStarFacS(ii, 1-j,bi,bj)= 1.
                0290           IF (kSurfS(1-i,jj, bi,bj).GT.Nr) rStarFacS(1-i,jj, bi,bj)= 1.
                0291           IF (kSurfS(ii, jj, bi,bj).GT.Nr) rStarFacS(ii, jj, bi,bj)= 1.
f450521a99 Jean*0292 
6d57af55dc Jean*0293          ENDDO
                0294         ENDDO
f450521a99 Jean*0295 #endif /* W2_FILL_NULL_REGIONS */
6d57af55dc Jean*0296 #endif /* ALLOW_EXCH2 */
                0297 
ccbf39d1a4 Jean*0298         DO j=1-OLy,sNy+OLy
                0299          DO i=1-OLx,sNx+OLx
9dbedd0b0f Jean*0300            rStarDhCDt(i,j,bi,bj)=(rStarFacC(i,j,bi,bj)
ccbf39d1a4 Jean*0301      &                           -rStarExpC(i,j,bi,bj))/deltaTFreeSurf
00b29feb62 Jean*0302            rStarDhWDt(i,j,bi,bj)=(rStarFacW(i,j,bi,bj)
ccbf39d1a4 Jean*0303      &                           -rStarExpW(i,j,bi,bj))/deltaTFreeSurf
00b29feb62 Jean*0304            rStarDhSDt(i,j,bi,bj)=(rStarFacS(i,j,bi,bj)
ccbf39d1a4 Jean*0305      &                           -rStarExpS(i,j,bi,bj))/deltaTFreeSurf
00b29feb62 Jean*0306            rStarExpC(i,j,bi,bj) = rStarFacC(i,j,bi,bj)
                0307      &                          / rStarExpC(i,j,bi,bj)
                0308            rStarExpW(i,j,bi,bj) = rStarFacW(i,j,bi,bj)
                0309      &                          / rStarExpW(i,j,bi,bj)
                0310            rStarExpS(i,j,bi,bj) = rStarFacS(i,j,bi,bj)
                0311      &                          / rStarExpS(i,j,bi,bj)
                0312          ENDDO
                0313         ENDDO
                0314 
ccbf39d1a4 Jean*0315         IF ( fluidIsAir ) THEN
                0316          DO j=1-OLy,sNy+OLy
                0317           DO i=1-OLx,sNx+OLx
                0318            pStarFacK(i,j,bi,bj) = rStarFacC(i,j,bi,bj)**atm_kappa
                0319           ENDDO
                0320          ENDDO
                0321 #ifdef ALLOW_AUTODIFF
                0322         ELSE
                0323          DO j=1-OLy,sNy+OLy
                0324           DO i=1-OLx,sNx+OLx
                0325            pStarFacK(i,j,bi,bj) = 1. _d 0
                0326           ENDDO
                0327          ENDDO
                0328 #endif
                0329         ENDIF
                0330 
00b29feb62 Jean*0331 C-    end 2nd bi,bj loop.
                0332         ENDDO
                0333        ENDDO
                0334 
6d57af55dc Jean*0335 #ifdef ALLOW_DEBUG
                0336       IF (debugMode) CALL DEBUG_LEAVE('CALC_R_STAR',myThid)
                0337 #endif
                0338 
00b29feb62 Jean*0339 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0340 #endif /* NONLIN_FRSURF */
                0341 
                0342       RETURN
                0343       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!