model/src/ini_linear_phisurf.F

b6b8988e60 Jean*0001 #include "PACKAGES_CONFIG.h"
074ef64531 Jean*0002 #include "CPP_OPTIONS.h"
                0003
9366854e02 Chri*0004 CBOP
                0005 C     !ROUTINE: INI_LINEAR_PHISURF
                0006 C     !INTERFACE:
074ef64531 Jean*0007       SUBROUTINE INI_LINEAR_PHISURF( myThid )
                0008
9366854e02 Chri*0009 C     !DESCRIPTION: \bv
                0010 C     *==========================================================*
05f412be17 Jean*0011 C     | SUBROUTINE INI_LINEAR_PHISURF
                0012 C     | o Initialise the Linear Relation Phi_surf(eta)
9366854e02 Chri*0013 C     *==========================================================*
c8f29584cf Jean*0014 C     | Initialise -Buoyancy at surface level (Bo_surf)
05f412be17 Jean*0015 C     |  to setup the Linear relation: Phi_surf(eta)=Bo_surf*eta
                0016 C     | Initialise phi0surf = starting point for integrating
463053c692 Jean*0017 C     |                       phiHyd (= phiHyd at r=RoSurf)
9366854e02 Chri*0018 C     *==========================================================*
                0019 C     \ev
                0020
                0021 C     !USES:
                0022       IMPLICIT NONE
074ef64531 Jean*0023 C     === Global variables ===
                0024 #include "SIZE.h"
                0025 #include "EEPARAMS.h"
                0026 #include "PARAMS.h"
                0027 #include "GRID.h"
                0028 #include "SURFACE.h"
                0029
9366854e02 Chri*0030 C     !INPUT/OUTPUT PARAMETERS:
aa03c27196 Jean*0031 C     myThid :: my Thread Id number
074ef64531 Jean*0032       INTEGER myThid
                0033
46ac67bdf3 Jean*0034 C     == Local variables in common ==
3bffe370e9 Jean*0035 C     topoHloc had to be in common for multi threading but no longer
                0036 C     needed since MDSIO now allows (2009/06/07) to write local arrays
46ac67bdf3 Jean*0037
9366854e02 Chri*0038 C     !LOCAL VARIABLES:
3bffe370e9 Jean*0039 C     topoHloc  :: Temporary array used to write surface topography
e437317be2 Jean*0040 C     bi,bj  :: tile indices
065ea7c980 Jean*0041 C     i,j,k  :: Loop counters
9f5240b52a Jean*0042 #ifndef ALLOW_AUTODIFF
e437317be2 Jean*0043       _RS topoHloc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
9f5240b52a Jean*0044 #endif
074ef64531 Jean*0045       INTEGER bi, bj
065ea7c980 Jean*0046       INTEGER i, j, k
69d75ba75d Jean*0047       _RL     pLoc, rhoLoc
074ef64531 Jean*0048       _RL     dPIdp
9366854e02 Chri*0049 CEOP
074ef64531 Jean*0050
065ea7c980 Jean*0051 C--   Initialisation
b6b8988e60 Jean*0052 #ifdef ALLOW_AUTODIFF
5c39e917f7 Patr*0053       DO bj=myByLo(myThid),myByHi(myThid)
                0054        DO bi=myBxLo(myThid),myBxHi(myThid)
8885a61e2b Jean*0055         DO j=1-OLy,sNy+OLy
                0056          DO i=1-OLx,sNx+OLx
065ea7c980 Jean*0057            Bo_surf(i,j,bi,bj)  = 0. _d 0
                0058            recip_Bo(i,j,bi,bj) = 0. _d 0
                0059          ENDDO
                0060         ENDDO
8885a61e2b Jean*0061        ENDDO
                0062       ENDDO
b6b8988e60 Jean*0063 #endif /* ALLOW_AUTODIFF */
3bffe370e9 Jean*0064
                0065 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0066
c8f29584cf Jean*0067 C-- Initialise -Buoyancy at surface level : Bo_surf
05f412be17 Jean*0068 C   Bo_surf is defined as d/dr(Phi_surf) and set to g/z2rUnit with
                0069 C     z2rUnit = conversion factor from z-unit to r-unit: [z] * z2rUnit = [r]
                0070 C   an accurate formulation includes P_surf and T,S_ref effects on rho_surf:
074ef64531 Jean*0071 C    (setting uniformLin_PhiSurf=.FALSE.):
c8f29584cf Jean*0072 C    z-coord (z2rUnit=1): Bo_surf = - Buoyancy
05f412be17 Jean*0073 C                                 = g * rho_surf(Tref,Sref,pSurf_0)/rho_0
                0074 C    p-coord (z2rUnit=rho*g): Bo_surf = 1/rho(Tref(ksurf),pSurf_0)
                0075 C   Note: on Phi_surf splitting : Non-linear Time-dependent effects on B_surf
                0076 C   [through eta & (T-tRef)_surf] are included in PhiHyd rather than in Bo_surf
074ef64531 Jean*0077 C--
aa03c27196 Jean*0078       IF ( usingZCoords ) THEN
074ef64531 Jean*0079 C-  gBaro = gravity (except for External mode test with reduced gravity)
                0080         DO bj=myByLo(myThid),myByHi(myThid)
                0081          DO bi=myBxLo(myThid),myBxHi(myThid)
cfbe180681 Jean*0082           DO j=1-OLy,sNy+OLy
                0083            DO i=1-OLx,sNx+OLx
                0084              Bo_surf(i,j,bi,bj) = gBaro
                0085              recip_Bo(i,j,bi,bj) = 1. _d 0 / gBaro
074ef64531 Jean*0086            ENDDO
                0087           ENDDO
                0088          ENDDO
                0089         ENDDO
                0090       ELSEIF ( uniformLin_PhiSurf ) THEN
                0091 C-  use a linear (in ps) uniform relation : Phi'_surf = 1/rhoConst * ps'_surf
                0092         DO bj=myByLo(myThid),myByHi(myThid)
                0093          DO bi=myBxLo(myThid),myBxHi(myThid)
cfbe180681 Jean*0094           DO j=1-OLy,sNy+OLy
                0095            DO i=1-OLx,sNx+OLx
                0096 c            Bo_surf(i,j,bi,bj)  = rVel2wUnit(1)*gravity
                0097 c            recip_Bo(i,j,bi,bj) = wUnit2rVel(1)*recip_gravity
                0098              Bo_surf(i,j,bi,bj)  = recip_rhoConst
                0099              recip_Bo(i,j,bi,bj) = rhoConst
074ef64531 Jean*0100            ENDDO
                0101           ENDDO
                0102          ENDDO
                0103         ENDDO
aa03c27196 Jean*0104       ELSEIF ( fluidIsWater ) THEN
cfbe180681 Jean*0105 C--   More precise than uniformLin_PhiSurf case but inconsistent
                0106 C     with nonlinFreeSurf=4 in CALC_PHI_HYD (eta contribution to phiHyd)
60c223928f Mart*0107         DO bj=myByLo(myThid),myByHi(myThid)
                0108          DO bi=myBxLo(myThid),myBxHi(myThid)
cfbe180681 Jean*0109           DO j=1-OLy,sNy+OLy
                0110            DO i=1-OLx,sNx+OLx
                0111             IF ( Ro_surf(i,j,bi,bj).GT.0. _d 0
                0112      &          .AND. kSurfC(i,j,bi,bj).LE.Nr ) THEN
                0113              pLoc = Ro_surf(i,j,bi,bj)
1177955fb6 Jean*0114 #ifdef ALLOW_OPENAD
                0115              CALL FIND_RHO_SCALAR(
cfbe180681 Jean*0116      I            tRef(kSurfC(i,j,bi,bj)),
                0117      I            sRef(kSurfC(i,j,bi,bj)),
1177955fb6 Jean*0118      I            pLoc,
                0119      O            rhoLoc, myThid )
                0120 #else /* ALLOW_OPENAD */
cfbe180681 Jean*0121              k = kSurfC(i,j,bi,bj)
05f412be17 Jean*0122              CALL FIND_RHO_SCALAR(
69d75ba75d Jean*0123      I            tRef(k), sRef(k), pLoc,
                0124      O            rhoLoc, myThid )
1177955fb6 Jean*0125 #endif /* ALLOW_OPENAD */
05f412be17 Jean*0126              IF ( rhoLoc .EQ. 0. _d 0 ) THEN
cfbe180681 Jean*0127               Bo_surf(i,j,bi,bj) = 0. _d 0
05f412be17 Jean*0128              ELSE
cfbe180681 Jean*0129               Bo_surf(i,j,bi,bj) = 1. _d 0/rhoLoc
05f412be17 Jean*0130              ENDIF
cfbe180681 Jean*0131              recip_Bo(i,j,bi,bj) =  rhoLoc
60c223928f Mart*0132             ELSE
cfbe180681 Jean*0133               Bo_surf(i,j,bi,bj)  = 0. _d 0
                0134               recip_Bo(i,j,bi,bj) = 0. _d 0
60c223928f Mart*0135             ENDIF
                0136            ENDDO
                0137           ENDDO
                0138          ENDDO
                0139         ENDDO
aa03c27196 Jean*0140       ELSEIF ( fluidIsAir ) THEN
074ef64531 Jean*0141 C-  use a linearized (in ps) Non-uniform relation : Bo_surf(Po_surf,tRef_surf)
cfbe180681 Jean*0142 C    Bo = d/d_p(Phi_surf) = tRef_surf*d/d_p(PI) ; PI = Cp*(p/Po)^kappa
                0143 C   and atm_Cp*atm_kappa = atm_Rd
074ef64531 Jean*0144         DO bj=myByLo(myThid),myByHi(myThid)
                0145          DO bi=myBxLo(myThid),myBxHi(myThid)
cfbe180681 Jean*0146           IF ( select_rStar.GE.1 .OR. selectSigmaCoord.GE.1 ) THEN
                0147 C-      isothermal (theta=const) reference state
                0148            DO j=1-OLy,sNy+OLy
                0149             DO i=1-OLx,sNx+OLx
                0150              IF ( Ro_surf(i,j,bi,bj).GT.0. _d 0
                0151      &          .AND. kSurfC(i,j,bi,bj).LE.Nr ) THEN
                0152               dPIdp = (atm_Rd/atm_Po)*
                0153      &         (Ro_surf(i,j,bi,bj)/atm_Po)**(atm_kappa-1. _d 0)
                0154               Bo_surf(i,j,bi,bj) = dPIdp*thetaConst
                0155               recip_Bo(i,j,bi,bj) = 1. _d 0 / Bo_surf(i,j,bi,bj)
                0156              ELSE
                0157               Bo_surf(i,j,bi,bj) = 0.
                0158               recip_Bo(i,j,bi,bj) = 0.
                0159              ENDIF
                0160             ENDDO
074ef64531 Jean*0161            ENDDO
cfbe180681 Jean*0162           ELSE
                0163 C-      horizontally uniform (tRef) reference state
                0164            DO j=1-OLy,sNy+OLy
                0165             DO i=1-OLx,sNx+OLx
                0166              IF ( Ro_surf(i,j,bi,bj).GT.0. _d 0
                0167      &          .AND. kSurfC(i,j,bi,bj).LE.Nr ) THEN
                0168               dPIdp = (atm_Rd/atm_Po)*
                0169      &         (Ro_surf(i,j,bi,bj)/atm_Po)**(atm_kappa-1. _d 0)
                0170               Bo_surf(i,j,bi,bj) = dPIdp*tRef(kSurfC(i,j,bi,bj))
                0171               recip_Bo(i,j,bi,bj) = 1. _d 0 / Bo_surf(i,j,bi,bj)
                0172              ELSE
                0173               Bo_surf(i,j,bi,bj) = 0.
                0174               recip_Bo(i,j,bi,bj) = 0.
                0175              ENDIF
                0176             ENDDO
                0177            ENDDO
                0178           ENDIF
074ef64531 Jean*0179          ENDDO
                0180         ENDDO
60c223928f Mart*0181       ELSE
                0182         STOP 'INI_LINEAR_PHISURF: We should never reach this point!'
074ef64531 Jean*0183       ENDIF
                0184
                0185 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0186
05f412be17 Jean*0187 C--   Update overlap regions (jmc: is it really needed ?)
cfbe180681 Jean*0188 c     _EXCH_XY_RL(Bo_surf, myThid)
                0189 c     _EXCH_XY_RL(recip_Bo, myThid)
074ef64531 Jean*0190
aa03c27196 Jean*0191       IF ( usingPCoords .AND. .NOT.uniformLin_PhiSurf ) THEN
05f412be17 Jean*0192         CALL WRITE_FLD_XY_RL( 'Bo_surf',' ',Bo_surf,0,myThid)
074ef64531 Jean*0193       ENDIF
                0194
463053c692 Jean*0195 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0196
05f412be17 Jean*0197 C--   Initialise phi0surf: used for atmos. surf. P-loading (ocean, z-coord)
463053c692 Jean*0198 C                               or topographic geopotential anom. (p-coord)
                0199
0320e25227 Mart*0200       DO bj=myByLo(myThid),myByHi(myThid)
                0201        DO bi=myBxLo(myThid),myBxHi(myThid)
                0202         DO j=1-OLy,sNy+OLy
                0203          DO i=1-OLx,sNx+OLx
                0204           phi0surf(i,j,bi,bj) = 0.
463053c692 Jean*0205          ENDDO
                0206         ENDDO
0320e25227 Mart*0207        ENDDO
                0208       ENDDO
                0209
                0210       IF ( geoPotAnomFile .NE. ' '  ) THEN
                0211        CALL READ_FLD_XY_RS( geoPotAnomFile, ' ', phi0Surf, 0, myThid )
                0212       ENDIF
                0213       CALL EXCH_XY_RS( phi0Surf  , myThid )
463053c692 Jean*0214
aa03c27196 Jean*0215       IF ( fluidIsAir .AND. topoFile.NE.' ' ) THEN
463053c692 Jean*0216
b6b8988e60 Jean*0217 #ifdef ALLOW_AUTODIFF
ecaea33887 Patr*0218          STOP 'CANNOT PRESENTLY USE THIS OPTION WITH ADJOINT'
                0219 #else
                0220
05f412be17 Jean*0221 C--   Compute topoH = PhiRef(Po_surf)/g ; is different from original
46ac67bdf3 Jean*0222 C      topoZ(read from file) because of truncation of Po_surf.
                0223 C     NOTE: not clear for now which topoZ needs to be saved in common block
                0224 C--   AND set phi0surf = starting point for integrating Geopotential;
                0225
05f412be17 Jean*0226         CALL INI_P_GROUND( -2,
                0227      O                     topoHloc,
463053c692 Jean*0228      I                     Ro_surf, myThid )
                0229
46ac67bdf3 Jean*0230        IF (selectFindRoSurf.NE.0) THEN
463053c692 Jean*0231         _EXCH_XY_RS(phi0surf, myThid)
05f412be17 Jean*0232         CALL WRITE_FLD_XY_RS( 'phi0surf',' ',phi0surf,0,myThid)
46ac67bdf3 Jean*0233        ENDIF
463053c692 Jean*0234
06bb0cec77 Jean*0235         CALL WRITE_FLD_XY_RS( 'topo_H',' ',topoHloc,0,myThid)
                0236
b6b8988e60 Jean*0237 #endif /* ALLOW_AUTODIFF */
ecaea33887 Patr*0238
463053c692 Jean*0239       ENDIF
                0240
                0241 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
074ef64531 Jean*0242       RETURN
                0243       END