solid-body.cs-32x32x1/code/ini_psurf.F

aef96ed361 Jean*0001 #include "PACKAGES_CONFIG.h"
905f660335 Alis*0002 #include "CPP_OPTIONS.h"
                0003
aef96ed361 Jean*0004 CBOP
                0005 C     !ROUTINE: INI_PSURF
                0006 C     !INTERFACE:
905f660335 Alis*0007       SUBROUTINE INI_PSURF( myThid )
aef96ed361 Jean*0008
                0009 C     !DESCRIPTION: \bv
                0010 C     *==========================================================*
905f660335 Alis*0011 C     | SUBROUTINE INI_PSURF                                     |
                0012 C     | o Set model initial free-surface height/pressure.        |
aef96ed361 Jean*0013 C     *==========================================================*
905f660335 Alis*0014 C     | There are several options for setting the initial        |
                0015 C     | surface displacement (r unit) field.                     |
                0016 C     |  1. Inline code                                          |
                0017 C     |  2. Two-dimensional data from a file.                    |
aef96ed361 Jean*0018 C     *==========================================================*
                0019 C     \ev
905f660335 Alis*0020
aef96ed361 Jean*0021 C     !USES:
                0022       IMPLICIT NONE
905f660335 Alis*0023 C     === Global variables ===
                0024 #include "SIZE.h"
                0025 #include "EEPARAMS.h"
                0026 #include "PARAMS.h"
                0027 #include "GRID.h"
                0028 #include "DYNVARS.h"
aef96ed361 Jean*0029 #include "SURFACE.h"
                0030 #ifdef ALLOW_CD_CODE
                0031 # include "CD_CODE_VARS.h"
                0032 #endif
                0033 #ifdef ALLOW_SHELFICE
                0034 # include "SHELFICE.h"
                0035 #endif
905f660335 Alis*0036
aef96ed361 Jean*0037 C     !INPUT/OUTPUT PARAMETERS:
                0038 C     myThid :: my Thread Id number
905f660335 Alis*0039       INTEGER myThid
                0040
aef96ed361 Jean*0041 C     !LOCAL VARIABLES:
                0042 C     bi,bj  :: tiles indices
                0043 C     i, j   :: Loop counters
905f660335 Alis*0044       INTEGER bi, bj
aef96ed361 Jean*0045       INTEGER i, j
                0046       _RL omegaPrime, snFac, psFac
                0047 CEOP
905f660335 Alis*0048
                0049 C--   Initialise surface position anomaly to zero
aef96ed361 Jean*0050       omegaPrime = 80. _d 0 / rSphere
                0051 c     psFac = -(rSphere**2)*omegaPrime*(Omega+omegaPrime)
                0052 C     previous expression above was missing one "half" factor:
                0053       psfac = -rhoConst*(rSphere*rSphere)
                0054      &                 *omegaPrime*( Omega + omegaPrime*0.5 _d 0 )
                0055       snFac = 1. _d 0 / (4. _d 0*Omega*Omega)
905f660335 Alis*0056       DO bj = myByLo(myThid), myByHi(myThid)
                0057        DO bi = myBxLo(myThid), myBxHi(myThid)
aef96ed361 Jean*0058         DO j=1-OLy,sNy+OLy
                0059          DO i=1-OLx,sNx+OLx
                0060           etaN(i,j,bi,bj) = 0. _d 0
                0061      &     + psFac*( snFac*fCori(i,j,bi,bj)*fCori(i,j,bi,bj)
                0062      &             - 1. _d 0 / 3. _d 0 )
905f660335 Alis*0063          ENDDO
                0064         ENDDO
                0065        ENDDO
                0066       ENDDO
                0067 C     Read an initial state
aef96ed361 Jean*0068       IF ( pSurfInitFile .NE. ' ' ) THEN
905f660335 Alis*0069        CALL READ_FLD_XY_RL( pSurfInitFile, ' ', etaN, 0, myThid )
aef96ed361 Jean*0070 C      fill the overlap (+ BARRIER)
                0071        _EXCH_XY_RL(etaN, myThid)
905f660335 Alis*0072       ENDIF
                0073
138482fdf6 Ed H*0074 #ifdef ALLOW_CD_CODE
aef96ed361 Jean*0075 C--   By default, initialize etaNm1 with etaN :
905f660335 Alis*0076       DO bj=myByLo(myThid),myByHi(myThid)
                0077        DO bi=myBxLo(myThid),myBxHi(myThid)
aef96ed361 Jean*0078         DO j=1-OLy,sNy+OLy
                0079          DO i=1-OLx,sNx+OLx
                0080           etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
905f660335 Alis*0081          ENDDO
                0082         ENDDO
                0083        ENDDO
                0084       ENDDO
                0085 #endif
                0086
aef96ed361 Jean*0087 #ifdef EXACT_CONSERV
                0088 C--   By default, initialize etaH with etaN :
                0089       DO bj=myByLo(myThid),myByHi(myThid)
                0090        DO bi=myBxLo(myThid),myBxHi(myThid)
                0091         DO j=1-OLy,sNy+OLy
                0092          DO i=1-OLx,sNx+OLx
                0093           etaH(i,j,bi,bj) = etaN(i,j,bi,bj)
                0094           etaHnm1(i,j,bi,bj) = etaN(i,j,bi,bj)
                0095           dEtaHdt(i,j,bi,bj) = 0. _d 0
                0096          ENDDO
                0097         ENDDO
                0098        ENDDO
                0099       ENDDO
                0100 #endif /* EXACT_CONSERV */
                0101
                0102 #ifdef ALLOW_SHELFICE
                0103       IF ( useShelfIce .AND. usingZCoords ) THEN
                0104         DO bj=myByLo(myThid),myByHi(myThid)
                0105          DO bi=myBxLo(myThid),myBxHi(myThid)
                0106            DO j=1-OLy,sNy+OLy
                0107             DO i=1-OLx,sNx+OLx
                0108               phi0surf(i,j,bi,bj) = phi0surf(i,j,bi,bj)
                0109      &          + shelficeLoadAnomaly(i,j,bi,bj)*recip_rhoConst
                0110             ENDDO
                0111            ENDDO
                0112          ENDDO
                0113         ENDDO
                0114       ENDIF
                0115 #endif /* ALLOW_SHELFICE */
                0116
905f660335 Alis*0117       RETURN
                0118       END