model/src/calc_surf_dr.F

6c49d7f2f2 Jean*0001 #include "PACKAGES_CONFIG.h"
be47244872 Jean*0002 #include "CPP_OPTIONS.h"
                0003
9366854e02 Chri*0004 CBOP
                0005 C     !ROUTINE: CALC_SURF_DR
                0006 C     !INTERFACE:
6c49d7f2f2 Jean*0007       SUBROUTINE CALC_SURF_DR( etaFld,
                0008      I                         myTime, myIter, myThid )
9366854e02 Chri*0009 C     !DESCRIPTION: \bv
                0010 C     *==========================================================*
dff69a4ed4 Jean*0011 C     | SUBROUTINE CALC_SURF_DR
                0012 C     | o Calculate the new surface level thickness according to
                0013 C     |   the surface r-position  (Non-Linear Free-Surf)
                0014 C     | o take decision if grid box becomes too thin or too thick
9366854e02 Chri*0015 C     *==========================================================*
                0016 C     \ev
be47244872 Jean*0017
9366854e02 Chri*0018 C     !USES:
                0019       IMPLICIT NONE
be47244872 Jean*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
9366854e02 Chri*0027 C     !INPUT/OUTPUT PARAMETERS:
be47244872 Jean*0028 C     == Routine arguments ==
6c49d7f2f2 Jean*0029 C     etaFld    :: current eta field used to update the hFactor
                0030 C     myTime    :: current time in simulation
                0031 C     myIter    :: current iteration number in simulation
                0032 C     myThid    :: thread number for this instance of the routine.
31a65e9638 Jean*0033       _RL etaFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
be47244872 Jean*0034       _RL myTime
                0035       INTEGER myIter
                0036       INTEGER myThid
                0037
                0038 #ifdef NONLIN_FRSURF
                0039
9366854e02 Chri*0040 C     !LOCAL VARIABLES:
be47244872 Jean*0041 C     Local variables
9366854e02 Chri*0042 C     i,j,k,bi,bj  :: loop counter
                0043 C     rSurftmp     :: free surface r-position that is used to compute hFac_surf
                0044 C     adjust_nb_pt :: Nb of grid points where rSurf is adjusted (hFactInf)
                0045 C     adjust_volum :: adjustment effect on the volume (domain size)
da1ba8e1fa Jean*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
9e9a4cf401 Jean*0048       INTEGER i,j,bi,bj
2e7aec9951 dngo*0049       INTEGER ks, nTmp
                0050       INTEGER numbWrite, numbWrMax
dff69a4ed4 Jean*0051       _RL hFactmp, adjust_nb_pt, adjust_volum
23e0a84a0a Jean*0052       _RL rSurftmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
be47244872 Jean*0053       _RS hhm, hhp
55e9ea8a90 Jean*0054 c     CHARACTER*(MAX_LEN_MBUF) suff
9366854e02 Chri*0055 CEOP
da1ba8e1fa Jean*0056       DATA numbWrite / 0 /
                0057       numbWrMax = Nx*Ny
be47244872 Jean*0058
                0059 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0060
23e0a84a0a Jean*0061       adjust_nb_pt = 0.
                0062       adjust_volum = 0.
                0063
be47244872 Jean*0064       DO bj=myByLo(myThid), myByHi(myThid)
                0065        DO bi=myBxLo(myThid), myBxHi(myThid)
6c49d7f2f2 Jean*0066
72a058b866 Gael*0067 C--   before updating hFac_surfC/S/W save current fields
cb905f3199 Jean*0068         DO j=1-OLy,sNy+OLy
                0069          DO i=1-OLx,sNx+OLx
                0070            hFac_surfNm1C(i,j,bi,bj) = hFac_surfC(i,j,bi,bj)
                0071            hFac_surfNm1S(i,j,bi,bj) = hFac_surfS(i,j,bi,bj)
                0072            hFac_surfNm1W(i,j,bi,bj) = hFac_surfW(i,j,bi,bj)
72a058b866 Gael*0073          ENDDO
cb905f3199 Jean*0074         ENDDO
23e0a84a0a Jean*0075
                0076 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
6c49d7f2f2 Jean*0077 C-- Compute the new fractional thickness of surface level (kSurfC):
23e0a84a0a Jean*0078
                0079         DO j=0,sNy+1
                0080          DO i=0,sNx+1
                0081           rSurftmp(i,j) = Ro_surf(i,j,bi,bj)+etaFld(i,j,bi,bj)
6c49d7f2f2 Jean*0082           ks = kSurfC(i,j,bi,bj)
be47244872 Jean*0083           IF (ks.LE.Nr) THEN
cb905f3199 Jean*0084            IF ( rSurftmp(i,j).LT.Rmin_surf(i,j,bi,bj) ) THEN
23e0a84a0a Jean*0085 C-- Needs to do something :
872ca36aca Jean*0086             IF ( numbWrite.LE.numbWrMax .AND.
                0087      &           ( i.GE.1.AND.i.LE.sNx .AND.
                0088      &             j.GE.1.AND.j.LE.sNy ) ) THEN
55e9ea8a90 Jean*0089              numbWrite = numbWrite + 1
333081f7c0 Jean*0090              hFactmp = h0FacC(i,j,ks,bi,bj)
                0091      &         + ( rSurftmp(i,j) - Ro_surf(i,j,bi,bj) )*recip_drF(ks)
be47244872 Jean*0092              IF (hFactmp.LT.hFacInf) THEN
55e9ea8a90 Jean*0093               WRITE(errorMessageUnit,'(2A,6I4,I10)')
da1ba8e1fa Jean*0094      &         'WARNING: hFacC < hFacInf at:',
be47244872 Jean*0095      &         ' i,j,k,bi,bj,Thid,Iter=',i,j,ks,bi,bj,myThid,myIter
23e0a84a0a Jean*0096              ELSE
55e9ea8a90 Jean*0097               WRITE(errorMessageUnit,'(2A,6I4,I10)')
da1ba8e1fa Jean*0098      &         'WARNING: hFac < hFacInf near:',
be47244872 Jean*0099      &         ' i,j,k,bi,bj,Thid,Iter=',i,j,ks,bi,bj,myThid,myIter
                0100              ENDIF
cb905f3199 Jean*0101              WRITE(errorMessageUnit,'(A,2F10.6,1PE14.6)')
872ca36aca Jean*0102      &         ' hFac_n-1,hFac_n,eta =',
be47244872 Jean*0103      &          hfacC(i,j,ks,bi,bj), hFactmp, etaFld(i,j,bi,bj)
da1ba8e1fa Jean*0104             ENDIF
23e0a84a0a Jean*0105 C-- Decide to STOP :
cb905f3199 Jean*0106 c            WRITE(errorMessageUnit,'(A)')
da1ba8e1fa Jean*0107 c    &        'STOP in CALC_SURF_DR : too SMALL hFac !'
cb905f3199 Jean*0108 c            STOP 'ABNORMAL END: S/R CALC_SURF_DR'
                0109 C-- Or continue with Rmin_surf:
                0110             IF ( i.GE.1.AND.i.LE.sNx .AND.
                0111      &           j.GE.1.AND.j.LE.sNy ) THEN
                0112               adjust_nb_pt = adjust_nb_pt + 1.
                0113               adjust_volum = adjust_volum
23e0a84a0a Jean*0114      &          + rA(i,j,bi,bj)*(Rmin_surf(i,j,bi,bj)-rSurftmp(i,j))
cb905f3199 Jean*0115             ENDIF
                0116             rSurftmp(i,j) = Rmin_surf(i,j,bi,bj)
23e0a84a0a Jean*0117 C----------
                0118            ENDIF
                0119
                0120 C-- Set hFac_surfC :
333081f7c0 Jean*0121            hFac_surfC(i,j,bi,bj) = h0FacC(i,j,ks,bi,bj)
                0122      &            + ( rSurftmp(i,j) - Ro_surf(i,j,bi,bj)
                0123      &              )*recip_drF(ks)*maskC(i,j,ks,bi,bj)
23e0a84a0a Jean*0124
                0125 C-- Usefull warning when hFac becomes very large:
c69f1949e6 Jean*0126            IF ( numbWrite.LE.numbWrMax .AND.
                0127      &          hFac_surfC(i,j,bi,bj).GT.hFacSup ) THEN
55e9ea8a90 Jean*0128               numbWrite = numbWrite + 1
                0129               WRITE(errorMessageUnit,'(2A,6I4,I10)')
da1ba8e1fa Jean*0130      &         'WARNING: hFacC > hFacSup at:',
23e0a84a0a Jean*0131      &         ' i,j,k,bi,bj,Thid,Iter=',i,j,ks,bi,bj,myThid,myIter
55e9ea8a90 Jean*0132               WRITE(errorMessageUnit,'(A,2F10.6,1PE14.6)')
872ca36aca Jean*0133      &         ' hFac_n-1,hFac_n,eta =', hfacC(i,j,ks,bi,bj),
da1ba8e1fa Jean*0134      &          hFac_surfC(i,j,bi,bj), etaFld(i,j,bi,bj)
be47244872 Jean*0135            ENDIF
c69f1949e6 Jean*0136 C----------
be47244872 Jean*0137           ENDIF
23e0a84a0a Jean*0138
be47244872 Jean*0139          ENDDO
                0140         ENDDO
                0141
                0142 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0143 C-- Compute fractional thickness of surface level, for U & V point:
                0144
                0145         DO j=1,sNy
                0146          DO i=1,sNx+1
6c49d7f2f2 Jean*0147           ks = kSurfW(i,j,bi,bj)
be47244872 Jean*0148           IF (ks.LE.Nr) THEN
cb905f3199 Jean*0149 C-  allows hFacW to be larger than surrounding hFacC=1 @ edge of a step with
                0150 C   different kSurfC on either side (topo in p-coords, ice-shelf in z-coords)
                0151             hhm = rSurftmp(i-1,j)
                0152             hhp = rSurftmp(i,j)
                0153 C-  make sure hFacW is not larger than the 2 surrounding hFacC
                0154 c           hhm = rF(ks)
                0155 c           IF(ks.EQ.kSurfC(i-1,j,bi,bj)) hhm = rSurftmp(i-1,j)
                0156 c           hhp = rF(ks)
                0157 c           IF(ks.EQ.kSurfC(i,j,bi,bj))   hhp = rSurftmp(i,j)
31a65e9638 Jean*0158             hFac_surfW(i,j,bi,bj) = h0FacW(i,j,ks,bi,bj)
872ca36aca Jean*0159      &            + ( MIN(hhm,hhp) - rSurfW(i,j,bi,bj)
31a65e9638 Jean*0160      &              )*recip_drF(ks)*maskW(i,j,ks,bi,bj)
be47244872 Jean*0161           ENDIF
                0162          ENDDO
                0163         ENDDO
23e0a84a0a Jean*0164
be47244872 Jean*0165         DO j=1,sNy+1
                0166          DO i=1,sNx
6c49d7f2f2 Jean*0167           ks = kSurfS(i,j,bi,bj)
be47244872 Jean*0168           IF (ks.LE.Nr) THEN
cb905f3199 Jean*0169 C-  allows hFacS to be larger than surrounding hFacC=1 @ edge of a step with
                0170 C   different kSurfC on either side (topo in p-coords, ice-shelf in z-coords)
                0171             hhm = rSurftmp(i,j-1)
                0172             hhp = rSurftmp(i,j)
                0173 C-  make sure hFacS is not larger than the 2 surrounding hFacC
                0174 c           hhm = rF(ks)
                0175 c           IF(ks.EQ.kSurfC(i,j-1,bi,bj)) hhm = rSurftmp(i,j-1)
                0176 c           hhp = rF(ks)
                0177 c           IF(ks.EQ.kSurfC(i,j,bi,bj))   hhp = rSurftmp(i,j)
31a65e9638 Jean*0178             hFac_surfS(i,j,bi,bj) = h0FacS(i,j,ks,bi,bj)
872ca36aca Jean*0179      &            + ( MIN(hhm,hhp) - rSurfS(i,j,bi,bj)
31a65e9638 Jean*0180      &              )*recip_drF(ks)*maskS(i,j,ks,bi,bj)
be47244872 Jean*0181           ENDIF
                0182          ENDDO
                0183         ENDDO
23e0a84a0a Jean*0184
6c49d7f2f2 Jean*0185 #ifdef ALLOW_OBCS
                0186 C-- Apply OBC to hFac_surfW,S before the EXCH calls
                0187         IF ( useOBCS ) THEN
                0188           CALL OBCS_APPLY_SURF_DR(
                0189      I                    bi, bj, etaFld,
                0190      U                    hFac_surfC, hFac_surfW, hFac_surfS,
3814debfe6 Jean*0191      I                    myTime, myIter, myThid )
6c49d7f2f2 Jean*0192         ENDIF
                0193 #endif /* ALLOW_OBCS */
                0194
23e0a84a0a Jean*0195 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0196
                0197 C-    end bi,bj loop.
be47244872 Jean*0198        ENDDO
                0199       ENDDO
                0200
23e0a84a0a Jean*0201 C-- Global diagnostic :
12c8b75709 Jean*0202       _GLOBAL_SUM_RL( adjust_nb_pt , myThid )
31a65e9638 Jean*0203       IF ( adjust_nb_pt.GE.1. ) THEN
                0204         _GLOBAL_SUM_RL( adjust_volum , myThid )
55e9ea8a90 Jean*0205         _BEGIN_MASTER( myThid )
2e7aec9951 dngo*0206 C   just to avoid issue with OpenAD: copy to integer nTmp before printing
                0207         nTmp = NINT(adjust_nb_pt)
55e9ea8a90 Jean*0208         WRITE(standardMessageUnit,'(2(A,I10),1PE16.8)')
da1ba8e1fa Jean*0209      &    ' SURF_ADJUSTMENT: Iter=', myIter,
2e7aec9951 dngo*0210      &    ' Nb_pts,Vol=', nTmp, adjust_volum
da1ba8e1fa Jean*0211         _END_MASTER( myThid )
23e0a84a0a Jean*0212       ENDIF
                0213
12c8b75709 Jean*0214       _EXCH_XY_RS(hFac_surfC, myThid )
be47244872 Jean*0215       CALL EXCH_UV_XY_RS(hFac_surfW,hFac_surfS,.FALSE.,myThid)
                0216
23e0a84a0a Jean*0217 C-----
6c49d7f2f2 Jean*0218 C Note: testing kSurfW,S is equivalent to a full height mask
23e0a84a0a Jean*0219 C   ==> no need for applying the mask here.
                0220 C and with "partial thin wall" ==> mask could be applied in S/R UPDATE_SURF_DR
                0221 C-----
be47244872 Jean*0222
37d60d46f4 Jean*0223 c     IF ( myIter.GE.0 ) THEN
                0224 c       WRITE(suff,'(I10.10)') myIter
                0225 c       CALL WRITE_FLD_XY_RS( 'hFac_surfC.', suff, hFac_surfC,
55e9ea8a90 Jean*0226 c    &                         myIter, myThid )
37d60d46f4 Jean*0227 c     ENDIF
be47244872 Jean*0228
                0229 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0230 #endif /* NONLIN_FRSURF */
                0231
                0232       RETURN
                0233       END