model/src/set_ref_state.F

769c28d1ef Jean*0001 #include "CPP_OPTIONS.h"
                0002
                0003 CBOP
d87b427b2d Jean*0004 C     !ROUTINE: SET_REF_STATE
769c28d1ef Jean*0005 C     !INTERFACE:
d87b427b2d Jean*0006       SUBROUTINE SET_REF_STATE(
                0007      I                          myThid )
769c28d1ef Jean*0008 C     !DESCRIPTION: \bv
                0009 C     *==========================================================*
d87b427b2d Jean*0010 C     | SUBROUTINE SET_REF_STATE
12a0dca2f3 Jean*0011 C     | o Set reference potential at level center and
769c28d1ef Jean*0012 C     |   level interface, using tRef,sRef profiles.
12a0dca2f3 Jean*0013 C     | note: use same discretisation as in calc_phi_hyd
                0014 C     | o Set also reference stratification here (for implicit
f056be794c Jean*0015 C     |   Internal Gravity Waves) and units conversion factor
                0016 C     |   for vertical velocity (for Non-Hydrostatic in p)
                0017 C     |   since both use also the same reference density.
769c28d1ef Jean*0018 C     *==========================================================*
                0019 C     \ev
                0020
                0021 C     !USES:
                0022       IMPLICIT NONE
                0023 C     == Global variables ==
                0024 #include "SIZE.h"
                0025 #include "EEPARAMS.h"
                0026 #include "PARAMS.h"
                0027 #include "GRID.h"
12a0dca2f3 Jean*0028 #include "EOS.h"
769c28d1ef Jean*0029
                0030 C     !INPUT/OUTPUT PARAMETERS:
5b172de0d2 Jean*0031 C     myThid   :: my Thread Id number
769c28d1ef Jean*0032       INTEGER myThid
                0033
                0034 C     !LOCAL VARIABLES:
5b172de0d2 Jean*0035 C     msgBuf   :: Informational/error message buffer
                0036 C     pRefIntF :: reference pressure at level interface
769c28d1ef Jean*0037       CHARACTER*(MAX_LEN_MBUF) msgBuf
12a0dca2f3 Jean*0038       INTEGER k, ks, stdUnit
5b172de0d2 Jean*0039       _RL rHalf(2*Nr+1), pRefIntF(Nr+1)
                0040       _RL tLoc(Nr)
f056be794c Jean*0041       _RL pLoc, rhoUp, rhoDw, rhoLoc
                0042       _RL ddPI, conv_theta2T, thetaLoc
759f945c40 Jean*0043 C--
                0044       _RL     maxResid
                0045       INTEGER maxIterNb, belowCritNb
769c28d1ef Jean*0046 CEOP
                0047
                0048       _BEGIN_MASTER( myThid )
                0049
d87b427b2d Jean*0050 C--   Initialise:
12a0dca2f3 Jean*0051       DO k=1,2*Nr
                0052         phiRef(k) = 0.
769c28d1ef Jean*0053       ENDDO
                0054       stdUnit = standardMessageUnit
                0055
12a0dca2f3 Jean*0056       DO k=1,Nr
5b172de0d2 Jean*0057         rUnit2z(k) = 1. _d 0
                0058         z2rUnit(k) = 1. _d 0
                0059         rhoRef(k)  = rhoConst
                0060         dBdrRef(k) = 0. _d 0
                0061         pRef4EOS(k)  = 0. _d 0
12a0dca2f3 Jean*0062         rHalf(2*k)   = rC(k)
769c28d1ef Jean*0063       ENDDO
                0064
f056be794c Jean*0065       DO k=1,Nr+1
5b172de0d2 Jean*0066         pRefIntF(k)  = 0. _d 0
                0067         rHalf(2*k-1) = rF(k)
f056be794c Jean*0068         rVel2wUnit(k) = 1. _d 0
                0069         wUnit2rVel(k) = 1. _d 0
                0070       ENDDO
                0071
759f945c40 Jean*0072 C-    Initialise density factor for anelastic formulation:
d87b427b2d Jean*0073       DO k=1,Nr
                0074         rhoFacC(k) = 1. _d 0
                0075         recip_rhoFacC(k) = 1. _d 0
                0076       ENDDO
                0077       DO k=1,Nr+1
                0078         rhoFacF(k) = 1. _d 0
                0079         recip_rhoFacF(k) = 1. _d 0
                0080       ENDDO
                0081
769c28d1ef Jean*0082 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
759f945c40 Jean*0083 C--   Oceanic: define reference pressure/geo-potential vertical scale:
                0084       IF ( buoyancyRelation.EQ.'OCEANIC' ) THEN
2ad79bdf32 Jean*0085         phiRef(1)   = top_Pres*recip_rhoConst
5b172de0d2 Jean*0086         pRefIntF(1) = top_Pres
759f945c40 Jean*0087         IF ( gravityFile.EQ.' ' ) THEN
                0088          DO k=1,Nr
2ad79bdf32 Jean*0089           phiRef(2*k)   = phiRef(1)
                0090      &                  + (rC(k) - rF(1))*gravity*gravitySign
                0091           phiRef(2*k+1) = phiRef(1)
                0092      &                  + (rF(k+1)-rF(1))*gravity*gravitySign
759f945c40 Jean*0093 c         pRef4EOS(k)   = rhoConst*phiRef(2*k)
5b172de0d2 Jean*0094 c         pRefIntF(k+1) = rhoConst*phiRef(2*k+1)
759f945c40 Jean*0095 C note: just to get the same previous machine truncation:
5b172de0d2 Jean*0096           pRef4EOS(k)   = pRefIntF(1)
759f945c40 Jean*0097      &                  + rhoConst*(rC(k) - rF(1))*gravity*gravitySign
5b172de0d2 Jean*0098           pRefIntF(k+1) = pRefIntF(1)
759f945c40 Jean*0099      &                  + rhoConst*(rF(k+1)-rF(1))*gravity*gravitySign
                0100          ENDDO
                0101         ELSEIF ( integr_GeoPot.EQ.1 ) THEN
                0102          DO k=1,Nr
                0103           phiRef(2*k)   = phiRef(2*k-1)
                0104      &                           + halfRL*drF(k)*gravity*gravFacC(k)
                0105           phiRef(2*k+1) = phiRef(2*k-1) + drF(k)*gravity*gravFacC(k)
                0106           pRef4EOS(k)   = rhoConst*phiRef(2*k)
5b172de0d2 Jean*0107           pRefIntF(k+1) = rhoConst*phiRef(2*k+1)
759f945c40 Jean*0108          ENDDO
                0109         ELSE
2ad79bdf32 Jean*0110          phiRef(2)   = phiRef(1) + drC(1)*gravity*gravFacF(1)
759f945c40 Jean*0111          DO k=2,Nr
                0112           phiRef(2*k-1) = phiRef(2*k-2)
                0113      &                           + halfRL*drC(k)*gravity*gravFacF(k)
                0114           phiRef(2*k)   = phiRef(2*k-2) + drC(k)*gravity*gravFacF(k)
                0115          ENDDO
                0116          k = Nr
                0117          phiRef(2*k+1) = phiRef(2*k) + drC(k+1)*gravity*gravFacF(k+1)
                0118          DO k=1,Nr
                0119           pRef4EOS(k)   = rhoConst*phiRef(2*k)
5b172de0d2 Jean*0120           pRefIntF(k+1) = rhoConst*phiRef(2*k+1)
759f945c40 Jean*0121          ENDDO
                0122         ENDIF
                0123       ELSEIF ( buoyancyRelation.EQ.'OCEANICP' ) THEN
2ad79bdf32 Jean*0124         phiRef(2*Nr+1) = seaLev_Z*gravity
759f945c40 Jean*0125         DO k=1,Nr
2ad79bdf32 Jean*0126           phiRef(2*k)   = phiRef(2*Nr+1)
                0127      &                  - recip_rhoConst*( rC(k) - rF(Nr+1) )
                0128           phiRef(2*k-1) = phiRef(2*Nr+1)
                0129      &                  - recip_rhoConst*( rF(k) - rF(Nr+1) )
5b172de0d2 Jean*0130           pRef4EOS(k)   = rC(k)
                0131           pRefIntF(k)   = rF(k)
759f945c40 Jean*0132         ENDDO
5b172de0d2 Jean*0133         pRefIntF(Nr+1) = rF(Nr+1)
759f945c40 Jean*0134 c     ELSEIF (buoyancyRelation .EQ. 'ATMOSPHERIC') THEN
                0135 C     phiRef is computed later (see below)
                0136       ENDIF
                0137
                0138 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
12a0dca2f3 Jean*0139       IF ( eosType.EQ.'POLY3' ) THEN
                0140        IF ( implicitIntGravWave ) THEN
d87b427b2d Jean*0141          WRITE(msgBuf,'(2A)') 'SET_REF_STATE:',
                0142      &    ' need to compute reference density for Impl.IGW'
f056be794c Jean*0143          CALL PRINT_ERROR( msgBuf , myThid )
d87b427b2d Jean*0144          WRITE(msgBuf,'(2A)') 'SET_REF_STATE:',
                0145      &    ' but FIND_RHO_SCALAR(EOS="POLY3") not (yet) implemented'
f056be794c Jean*0146          CALL PRINT_ERROR( msgBuf , myThid )
d87b427b2d Jean*0147          STOP 'ABNORMAL END: S/R SET_REF_STATE'
f056be794c Jean*0148        ELSEIF ( nonHydrostatic .AND.
759f945c40 Jean*0149      &          buoyancyRelation.EQ.'OCEANICP' ) THEN
d87b427b2d Jean*0150          WRITE(msgBuf,'(2A)') 'SET_REF_STATE:',
                0151      &    ' need to compute reference density for Non-Hyd'
f056be794c Jean*0152          CALL PRINT_ERROR( msgBuf , myThid )
d87b427b2d Jean*0153          WRITE(msgBuf,'(2A)') 'SET_REF_STATE:',
                0154      &    ' but FIND_RHO_SCALAR(EOS="POLY3") not (yet) implemented'
f056be794c Jean*0155          CALL PRINT_ERROR( msgBuf , myThid )
d87b427b2d Jean*0156          STOP 'ABNORMAL END: S/R SET_REF_STATE'
12a0dca2f3 Jean*0157        ELSE
d87b427b2d Jean*0158          WRITE(msgBuf,'(2A)') 'SET_REF_STATE:',
                0159      &    ' Unable to compute reference stratification'
12a0dca2f3 Jean*0160          CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
f056be794c Jean*0161      &                       SQUEEZE_RIGHT , myThid )
d87b427b2d Jean*0162          WRITE(msgBuf,'(2A)') 'SET_REF_STATE:',
                0163      &    '  with EOS="POLY3" ; set dBdrRef(1:Nr) to zeros'
12a0dca2f3 Jean*0164          CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0165      &                       SQUEEZE_RIGHT , myThid)
                0166        ENDIF
                0167 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
759f945c40 Jean*0168       ELSEIF ( buoyancyRelation.EQ.'OCEANIC' ) THEN
12a0dca2f3 Jean*0169
759f945c40 Jean*0170 C--   Compute reference density profile
12a0dca2f3 Jean*0171         DO k=1,Nr
759f945c40 Jean*0172           pLoc = pRef4EOS(k)
12a0dca2f3 Jean*0173           CALL FIND_RHO_SCALAR(
                0174      I                          tRef(k), sRef(k), pLoc,
                0175      O                          rhoRef(k), myThid )
                0176         ENDDO
                0177
d6720d7423 Jean*0178         IF ( selectP_inEOS_Zc.GE.1 ) THEN
759f945c40 Jean*0179 C--   Find reference pressure in hydrostatic balance with updated ref density
                0180           maxResid = rhoConst* 1. _d -14
                0181           belowCritNb = 5
                0182           maxIterNb = 10*Nr
                0183           CALL FIND_HYD_PRESS_1D(
5b172de0d2 Jean*0184      O                            pRef4EOS, pRefIntF,
759f945c40 Jean*0185      U                            rhoRef,
                0186      I                            tRef, sRef, maxResid,
                0187      I                            belowCritNb, maxIterNb, myThid )
                0188         ENDIF
                0189
12a0dca2f3 Jean*0190 C--   Compute reference stratification: N^2 = -(g/rho_c) * d.rho/dz @ const. p
                0191         dBdrRef(1) = 0. _d 0
                0192         DO k=2,Nr
5b172de0d2 Jean*0193           pLoc = pRefIntF(k)
12a0dca2f3 Jean*0194           CALL FIND_RHO_SCALAR(
                0195      I                          tRef(k-1), sRef(k-1), pLoc,
                0196      O                          rhoUp, myThid )
                0197           CALL FIND_RHO_SCALAR(
                0198      I                          tRef(k), sRef(k), pLoc,
                0199      O                          rhoDw, myThid )
                0200           dBdrRef(k) = (rhoDw - rhoUp)*recip_drC(k)
759f945c40 Jean*0201      &               *recip_rhoConst*gravity*gravFacF(k)
                0202           IF ( eosType.EQ.'LINEAR' ) THEN
12a0dca2f3 Jean*0203 C- get more precise values (differences from above are due to machine round-off)
                0204             dBdrRef(k) = ( sBeta *(sRef(k)-sRef(k-1))
                0205      &                    -tAlpha*(tRef(k)-tRef(k-1))
                0206      &                   )*recip_drC(k)
759f945c40 Jean*0207      &                 *rhoNil*recip_rhoConst*gravity*gravFacF(k)
12a0dca2f3 Jean*0208           ENDIF
                0209         ENDDO
                0210
                0211 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
759f945c40 Jean*0212       ELSEIF ( buoyancyRelation.EQ.'OCEANICP' ) THEN
12a0dca2f3 Jean*0213
5b172de0d2 Jean*0214 C--   Compute reference density profile (from tRef,sRef):
12a0dca2f3 Jean*0215         DO k=1,Nr
5b172de0d2 Jean*0216           pLoc = pRef4EOS(k)
12a0dca2f3 Jean*0217           CALL FIND_RHO_SCALAR(
                0218      I                          tRef(k), sRef(k), pLoc,
                0219      O                          rhoRef(k), myThid )
5b172de0d2 Jean*0220 C-    set units convertion factor @ level center:
                0221 C       z2rUnit = gravity*rhoRef : dr [Pa] = dz [m] * z2rUnit
                0222 C       rUnit2z = 1/z2rUnit      : dz [m] = dr [Pa] * rUnit2z
                0223           z2rUnit(k) = gravity*rhoRef(k)
                0224           rUnit2z(k) = 1. _d 0 / z2rUnit(k)
12a0dca2f3 Jean*0225         ENDDO
                0226
                0227 C--   Compute reference stratification: -d.alpha/dp @ constant p
                0228         dBdrRef(1) = 0. _d 0
f056be794c Jean*0229         DO k=1,Nr+1
5b172de0d2 Jean*0230           pLoc = pRefIntF(k)
f056be794c Jean*0231           IF ( k.GE.2 )  CALL FIND_RHO_SCALAR(
                0232      I                             tRef(k-1), sRef(k-1), pLoc,
                0233      O                             rhoDw, myThid )
                0234           IF ( k.LE.Nr ) CALL FIND_RHO_SCALAR(
                0235      I                             tRef(k), sRef(k), pLoc,
                0236      O                             rhoUp, myThid )
                0237           IF ( k.GE.2 .AND. k.LE.Nr ) THEN
                0238             dBdrRef(k) = (rhoDw - rhoUp)*recip_drC(k)
b7e4f9b618 Jean*0239      &                 / (rhoDw*rhoUp)
                0240             rhoLoc = ( rhoDw + rhoUp )*0.5 _d 0
f056be794c Jean*0241           ELSEIF ( k.EQ.1 ) THEN
b7e4f9b618 Jean*0242             rhoLoc = rhoUp
f056be794c Jean*0243           ELSE
b7e4f9b618 Jean*0244             rhoLoc = rhoDw
f056be794c Jean*0245           ENDIF
                0246 C--   Units convertion factor for vertical velocity:
                0247 C       wUnit2rVel = gravity*rhoRef : rVel  [Pa/s] = wSpeed [m/s] * wUnit2rVel
                0248 C       rVel2wUnit = 1/rVel2wUnit   : wSpeed [m/s] = rVel  [Pa/s] * rVel2wUnit
                0249 C     note: wUnit2rVel & rVel2wUnit replace horiVertRatio & recip_horiVertRatio
                0250           wUnit2rVel(k) = gravity*rhoLoc
                0251           rVel2wUnit(k) = 1. _d 0 / wUnit2rVel(k)
12a0dca2f3 Jean*0252         ENDDO
                0253
                0254 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
759f945c40 Jean*0255       ELSEIF ( buoyancyRelation.EQ.'ATMOSPHERIC' ) THEN
12a0dca2f3 Jean*0256
                0257 C--   Compute reference stratification: -d.alpha/dp @ constant p
                0258         dBdrRef(1) = 0. _d 0
                0259         DO k=2,Nr
                0260           conv_theta2T = (rF(k)/atm_Po)**atm_kappa
                0261 c         dBdrRef(k) = (tRef(k) - tRef(k-1))*recip_drC(k)
                0262 c    &               * conv_theta2T*atm_Rd/rF(k)
                0263           ddPI=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
                0264      &                 -((rC( k )/atm_Po)**atm_kappa) )
                0265           dBdrRef(k) = (tRef(k) - tRef(k-1))*recip_drC(k)
                0266      &               * ddPI*recip_drC(k)
                0267         ENDDO
                0268
5b172de0d2 Jean*0269 C--   Set units convertion factor @ level center:
                0270 C       z2rUnit = gravity/alpha : dr [Pa] = dz [m] * z2rUnit
                0271 C       rUnit2z = alpha/gravity : dz [m] = dr [Pa] * rUnit2z
                0272         DO k=1,Nr
                0273           pRef4EOS(k)  = rC(k)
                0274           thetaLoc = tRef(k)
                0275           IF ( thetaLoc.GT.0. _d 0 .AND. rC(k).GT.0. _d 0 ) THEN
                0276             conv_theta2T = (rC(k)/atm_Po)**atm_kappa
                0277             z2rUnit(k) = gravity
                0278      &                 * rC(k)/(atm_Rd*conv_theta2T*thetaLoc)
                0279             rUnit2z(k) = 1. _d 0 / z2rUnit(k)
                0280           ENDIF
                0281         ENDDO
                0282
f056be794c Jean*0283 C--   Units convertion factor for vertical velocity:
                0284 C       wUnit2rVel = gravity/alpha : rVel  [Pa/s] = wSpeed [m/s] * wUnit2rVel
                0285 C       rVel2wUnit = alpha/gravity : wSpeed [m/s] = rVel  [Pa/s] * rVel2wUnit
                0286 C       with alpha = 1/rhoRef = (R.T/p) (ideal gas)
                0287 C     note: wUnit2rVel & rVel2wUnit replace horiVertRatio & recip_horiVertRatio
                0288         DO k=1,Nr+1
5b172de0d2 Jean*0289           pRefIntF(k)  = rF(k)
f056be794c Jean*0290           IF ( k.EQ.1 ) THEN
                0291             thetaLoc = tRef(k)
                0292           ELSEIF ( k.GT.Nr ) THEN
                0293             thetaLoc = tRef(k-1)
                0294           ELSE
                0295             thetaLoc = (tRef(k) + tRef(k-1))*0.5 _d 0
                0296           ENDIF
                0297           IF ( thetaLoc.GT.0. _d 0 .AND. rF(k).GT.0. _d 0 ) THEN
                0298             conv_theta2T  = (rF(k)/atm_Po)**atm_kappa
                0299             wUnit2rVel(k) = gravity
                0300      &                    * rF(k)/(atm_Rd*conv_theta2T*thetaLoc)
                0301             rVel2wUnit(k) = 1. _d 0 / wUnit2rVel(k)
                0302           ENDIF
                0303         ENDDO
                0304
12a0dca2f3 Jean*0305 C-    Compute Reference Geopotential at Half levels :
d533c5b431 Jean*0306 C      Tracer level k: phiRef(2k)  ;  Interface_W level k: phiRef(2k-1)
769c28d1ef Jean*0307
2ad79bdf32 Jean*0308        phiRef(1) = seaLev_Z*gravity
d533c5b431 Jean*0309        IF ( select_rStar.GE.1 .OR. selectSigmaCoord.GE.1 ) THEN
                0310 C-      isothermal (theta=const) reference state
                0311         DO k=1,Nr
                0312          tLoc(k) = thetaConst
                0313         ENDDO
                0314        ELSE
                0315 C-      horizontally uniform (tRef) reference state
                0316         DO k=1,Nr
                0317          tLoc(k) = tRef(k)
                0318         ENDDO
                0319        ENDIF
769c28d1ef Jean*0320
                0321        IF (integr_GeoPot.EQ.1) THEN
12a0dca2f3 Jean*0322 C-    Finite Volume Form, linear by half level :
                0323         DO k=1,2*Nr
                0324           ks = (k+1)/2
                0325           ddPI=atm_Cp*( ((rHalf( k )/atm_Po)**atm_kappa)
                0326      &                 -((rHalf(k+1)/atm_Po)**atm_kappa) )
d533c5b431 Jean*0327           phiRef(k+1) = phiRef(k)+ddPI*tLoc(ks)
769c28d1ef Jean*0328         ENDDO
                0329 C------
                0330        ELSE
12a0dca2f3 Jean*0331 C-    Finite Difference Form, linear between Tracer level :
                0332 C      works with integr_GeoPot = 0, 2 or 3
                0333         k = 1
                0334           ddPI=atm_Cp*( ((rF(k)/atm_Po)**atm_kappa)
                0335      &                 -((rC(k)/atm_Po)**atm_kappa) )
d533c5b431 Jean*0336           phiRef(2*k)   = phiRef(1) + ddPI*tLoc(k)
12a0dca2f3 Jean*0337         DO k=1,Nr-1
                0338           ddPI=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
                0339      &                 -((rC(k+1)/atm_Po)**atm_kappa) )
d533c5b431 Jean*0340           phiRef(2*k+1) = phiRef(2*k) + ddPI*0.5*tLoc(k)
12a0dca2f3 Jean*0341           phiRef(2*k+2) = phiRef(2*k)
d533c5b431 Jean*0342      &                  + ddPI*0.5*(tLoc(k)+tLoc(k+1))
769c28d1ef Jean*0343         ENDDO
12a0dca2f3 Jean*0344         k = Nr
                0345           ddPI=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
                0346      &                 -((rF(k+1)/atm_Po)**atm_kappa) )
d533c5b431 Jean*0347           phiRef(2*k+1) = phiRef(2*k) + ddPI*tLoc(k)
769c28d1ef Jean*0348 C------
                0349        ENDIF
                0350
                0351 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
12a0dca2f3 Jean*0352       ELSE
d87b427b2d Jean*0353         STOP 'SET_REF_STATE: Bad value of buoyancyRelation !'
12a0dca2f3 Jean*0354 C--   endif buoyancyRelation
                0355       ENDIF
769c28d1ef Jean*0356
12a0dca2f3 Jean*0357 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0358
d87b427b2d Jean*0359       IF ( usingZCoords .AND. rhoRefFile .NE. ' ' ) THEN
                0360 C--   anelastic formulation : set density factor from reference density profile
                0361 C       surface-interface rho-factor has to be 1:
                0362         rhoFacF(1)   = 1. _d 0
                0363 C       rhoFac(k) = density ratio between layer k and top interface
                0364         DO k=1,Nr
                0365           rhoFacC(k) = rho1Ref(k)/rhoConst
759f945c40 Jean*0366 c         rhoFacC(k) = rho1Ref(k)*recip_rhoConst
d87b427b2d Jean*0367         ENDDO
                0368         DO k=2,Nr
fd7bdc3397 Jean*0369 C       since rkSign=-1, recip_drC(k) = 1./(rC(k-1)-rC(k))
                0370           rhoFacF(k) = ( rhoFacC(k-1)*(rF(k)-rC(k))
                0371      &                 + rhoFacC(k)*(rC(k-1)-rF(k)) )*recip_drC(k)
d87b427b2d Jean*0372         ENDDO
                0373 C       extrapolate down to the bottom:
fd7bdc3397 Jean*0374         rhoFacF(Nr+1) = ( rhoFacC(Nr)*(rF(Nr+1)-rF(Nr))
                0375      &                  + rhoFacF(Nr)*(rC(Nr)-rF(Nr+1))
                0376      &                  ) / (rC(Nr)-rF(Nr))
d87b427b2d Jean*0377 C-      set reciprocal rho-factor:
                0378         DO k=1,Nr
                0379           recip_rhoFacC(k) = 1. _d 0/rhoFacC(k)
                0380         ENDDO
                0381         DO k=1,Nr+1
                0382           recip_rhoFacF(k) = 1. _d 0/rhoFacF(k)
                0383         ENDDO
                0384       ENDIF
                0385
                0386 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0387
12a0dca2f3 Jean*0388 C--   Write to check :
759f945c40 Jean*0389       IF ( gravityFile .NE. ' ' ) THEN
                0390 C-    write gravity vertical profile factor to binary file:
                0391         CALL WRITE_GLVEC_RL( 'GravFacC',' ',gravFacC, Nr, -1, myThid )
                0392         CALL WRITE_GLVEC_RL( 'GravFacF',' ',gravFacF,Nr+1,-1, myThid )
                0393       ENDIF
d6720d7423 Jean*0394       IF ( selectP_inEOS_Zc.GE.1 ) THEN
759f945c40 Jean*0395 C-    write (to binary file) Hyd-Pressure used in EOS:
                0396         CALL WRITE_GLVEC_RL( 'PRef4EOS',' ',pRef4EOS, Nr, -1, myThid )
5b172de0d2 Jean*0397 c       CALL WRITE_GLVEC_RL( 'PRefIntF',' ',pRefIntF,Nr+1,-1, myThid )
759f945c40 Jean*0398       ENDIF
                0399       IF ( buoyancyRelation.EQ.'ATMOSPHERIC' ) THEN
b2b9a7d77d Jean*0400        WRITE(msgBuf,'(A)') ' '
769c28d1ef Jean*0401        CALL PRINT_MESSAGE( msgBuf, stdUnit, SQUEEZE_RIGHT, myThid )
12a0dca2f3 Jean*0402        WRITE(msgBuf,'(A)')
d87b427b2d Jean*0403      &  'SET_REF_STATE: PhiRef/g [m] at level Center (integer)'
769c28d1ef Jean*0404        CALL PRINT_MESSAGE( msgBuf, stdUnit, SQUEEZE_RIGHT, myThid )
12a0dca2f3 Jean*0405        WRITE(msgBuf,'(A)')
769c28d1ef Jean*0406      &  '                     and at level Interface (half-int.) :'
                0407        CALL PRINT_MESSAGE( msgBuf, stdUnit, SQUEEZE_RIGHT, myThid )
12a0dca2f3 Jean*0408        DO k=1,2*Nr+1
b7e4f9b618 Jean*0409         WRITE(msgBuf,'(A,F5.1,A,F15.1,A,F13.3)')
12a0dca2f3 Jean*0410      &    ' K=',k*0.5,'  ;  r=',rHalf(k),'  ;  phiRef/g=',
                0411      &    phiRef(k)*recip_gravity
769c28d1ef Jean*0412         CALL PRINT_MESSAGE(msgBuf, stdUnit, SQUEEZE_RIGHT, myThid )
                0413        ENDDO
d533c5b431 Jean*0414       ELSE
fd7bdc3397 Jean*0415 C-    Write reference density to binary file :
                0416         CALL WRITE_GLVEC_RL( 'RhoRef', ' ', rhoRef, Nr, -1, myThid )
                0417       ENDIF
                0418       IF ( usingZCoords .AND. rhoRefFile .NE. ' ' ) THEN
                0419 C-    Write Anelastic density factor to binary file :
                0420         CALL WRITE_GLVEC_RL( 'RhoFacC',' ',rhoFacC, Nr ,  -1, myThid )
                0421         CALL WRITE_GLVEC_RL( 'RhoFacF',' ',rhoFacF, Nr+1, -1, myThid )
d533c5b431 Jean*0422       ENDIF
769c28d1ef Jean*0423
                0424 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0425
                0426       _END_MASTER( myThid )
06bb0cec77 Jean*0427       _BARRIER
769c28d1ef Jean*0428
                0429       RETURN
                0430       END