C $Header: /u/gcmpack/MITgcm/model/src/set_ref_state.F,v 1.10 2016/04/07 16:21:20 jmc Exp $
C $Name: $
#include "CPP_OPTIONS.h"
CBOP
C !ROUTINE: SET_REF_STATE
C !INTERFACE:
SUBROUTINE SET_REF_STATE(
I myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE SET_REF_STATE
C | o Set reference potential at level center and
C | level interface, using tRef,sRef profiles.
C | note: use same discretisation as in calc_phi_hyd
C | o Set also reference stratification here (for implicit
C | Internal Gravity Waves) and units conversion factor
C | for vertical velocity (for Non-Hydrostatic in p)
C | since both use also the same reference density.
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "EOS.h"
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
INTEGER myThid
C !LOCAL VARIABLES:
C == Local variables ==
C msgBuf :: Informational/error message buffer
CHARACTER*(MAX_LEN_MBUF) msgBuf
INTEGER k, ks, stdUnit
_RL rHalf(2*Nr+1), pRefLocF(Nr+1)
_RL rhoRef(Nr), tLoc(Nr)
_RL pLoc, rhoUp, rhoDw, rhoLoc
_RL ddPI, conv_theta2T, thetaLoc
C--
_RL maxResid
INTEGER maxIterNb, belowCritNb
CEOP
_BEGIN_MASTER( myThid )
C-- Initialise:
DO k=1,2*Nr
phiRef(k) = 0.
ENDDO
stdUnit = standardMessageUnit
DO k=1,Nr
rhoRef(k) = 0.
dBdrRef(k) = 0.
pRef4EOS(k) = 0.
rHalf(2*k-1) = rF(k)
rHalf(2*k) = rC(k)
ENDDO
rHalf(2*Nr+1) = rF(Nr+1)
DO k=1,Nr+1
rVel2wUnit(k) = 1. _d 0
wUnit2rVel(k) = 1. _d 0
ENDDO
C- Initialise density factor for anelastic formulation:
DO k=1,Nr
rhoFacC(k) = 1. _d 0
recip_rhoFacC(k) = 1. _d 0
ENDDO
DO k=1,Nr+1
rhoFacF(k) = 1. _d 0
recip_rhoFacF(k) = 1. _d 0
ENDDO
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C-- Oceanic: define reference pressure/geo-potential vertical scale:
IF ( buoyancyRelation.EQ.'OCEANIC' ) THEN
phiRef(1) = top_Pres*recip_rhoConst
pRefLocF(1) = top_Pres
IF ( gravityFile.EQ.' ' ) THEN
DO k=1,Nr
phiRef(2*k) = phiRef(1)
& + (rC(k) - rF(1))*gravity*gravitySign
phiRef(2*k+1) = phiRef(1)
& + (rF(k+1)-rF(1))*gravity*gravitySign
c pRef4EOS(k) = rhoConst*phiRef(2*k)
c pRefLocF(k+1) = rhoConst*phiRef(2*k+1)
C note: just to get the same previous machine truncation:
pRef4EOS(k) = pRefLocF(1)
& + rhoConst*(rC(k) - rF(1))*gravity*gravitySign
pRefLocF(k+1) = pRefLocF(1)
& + rhoConst*(rF(k+1)-rF(1))*gravity*gravitySign
ENDDO
ELSEIF ( integr_GeoPot.EQ.1 ) THEN
DO k=1,Nr
phiRef(2*k) = phiRef(2*k-1)
& + halfRL*drF(k)*gravity*gravFacC(k)
phiRef(2*k+1) = phiRef(2*k-1) + drF(k)*gravity*gravFacC(k)
pRef4EOS(k) = rhoConst*phiRef(2*k)
pRefLocF(k+1) = rhoConst*phiRef(2*k+1)
ENDDO
ELSE
phiRef(2) = phiRef(1) + drC(1)*gravity*gravFacF(1)
DO k=2,Nr
phiRef(2*k-1) = phiRef(2*k-2)
& + halfRL*drC(k)*gravity*gravFacF(k)
phiRef(2*k) = phiRef(2*k-2) + drC(k)*gravity*gravFacF(k)
ENDDO
k = Nr
phiRef(2*k+1) = phiRef(2*k) + drC(k+1)*gravity*gravFacF(k+1)
DO k=1,Nr
pRef4EOS(k) = rhoConst*phiRef(2*k)
pRefLocF(k+1) = rhoConst*phiRef(2*k+1)
ENDDO
ENDIF
ELSEIF ( buoyancyRelation.EQ.'OCEANICP' ) THEN
phiRef(2*Nr+1) = seaLev_Z*gravity
DO k=1,Nr
phiRef(2*k) = phiRef(2*Nr+1)
& - recip_rhoConst*( rC(k) - rF(Nr+1) )
phiRef(2*k-1) = phiRef(2*Nr+1)
& - recip_rhoConst*( rF(k) - rF(Nr+1) )
ENDDO
c ELSEIF (buoyancyRelation .EQ. 'ATMOSPHERIC') THEN
C phiRef is computed later (see below)
ENDIF
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
IF ( eosType.EQ.'POLY3' ) THEN
IF ( implicitIntGravWave ) THEN
WRITE(msgBuf,'(2A)') 'SET_REF_STATE:',
& ' need to compute reference density for Impl.IGW'
CALL PRINT_ERROR( msgBuf , myThid )
WRITE(msgBuf,'(2A)') 'SET_REF_STATE:',
& ' but FIND_RHO_SCALAR(EOS="POLY3") not (yet) implemented'
CALL PRINT_ERROR( msgBuf , myThid )
STOP 'ABNORMAL END: S/R SET_REF_STATE'
ELSEIF ( nonHydrostatic .AND.
& buoyancyRelation.EQ.'OCEANICP' ) THEN
WRITE(msgBuf,'(2A)') 'SET_REF_STATE:',
& ' need to compute reference density for Non-Hyd'
CALL PRINT_ERROR( msgBuf , myThid )
WRITE(msgBuf,'(2A)') 'SET_REF_STATE:',
& ' but FIND_RHO_SCALAR(EOS="POLY3") not (yet) implemented'
CALL PRINT_ERROR( msgBuf , myThid )
STOP 'ABNORMAL END: S/R SET_REF_STATE'
ELSE
WRITE(msgBuf,'(2A)') 'SET_REF_STATE:',
& ' Unable to compute reference stratification'
CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
& SQUEEZE_RIGHT , myThid )
WRITE(msgBuf,'(2A)') 'SET_REF_STATE:',
& ' with EOS="POLY3" ; set dBdrRef(1:Nr) to zeros'
CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
& SQUEEZE_RIGHT , myThid)
ENDIF
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
ELSEIF ( buoyancyRelation.EQ.'OCEANIC' ) THEN
C-- Compute reference density profile
DO k=1,Nr
pLoc = pRef4EOS(k)
CALL FIND_RHO_SCALAR(
I tRef(k), sRef(k), pLoc,
O rhoRef(k), myThid )
ENDDO
IF ( selectP_inEOS_Zc.GE.1 ) THEN
C-- Find reference pressure in hydrostatic balance with updated ref density
maxResid = rhoConst* 1. _d -14
belowCritNb = 5
maxIterNb = 10*Nr
CALL FIND_HYD_PRESS_1D(
O pRef4EOS, pRefLocF,
U rhoRef,
I tRef, sRef, maxResid,
I belowCritNb, maxIterNb, myThid )
ENDIF
C-- Compute reference stratification: N^2 = -(g/rho_c) * d.rho/dz @ const. p
dBdrRef(1) = 0. _d 0
DO k=2,Nr
pLoc = pRefLocF(k)
CALL FIND_RHO_SCALAR(
I tRef(k-1), sRef(k-1), pLoc,
O rhoUp, myThid )
CALL FIND_RHO_SCALAR(
I tRef(k), sRef(k), pLoc,
O rhoDw, myThid )
dBdrRef(k) = (rhoDw - rhoUp)*recip_drC(k)
& *recip_rhoConst*gravity*gravFacF(k)
IF ( eosType.EQ.'LINEAR' ) THEN
C- get more precise values (differences from above are due to machine round-off)
dBdrRef(k) = ( sBeta *(sRef(k)-sRef(k-1))
& -tAlpha*(tRef(k)-tRef(k-1))
& )*recip_drC(k)
& *rhoNil*recip_rhoConst*gravity*gravFacF(k)
ENDIF
ENDDO
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
ELSEIF ( buoyancyRelation.EQ.'OCEANICP' ) THEN
C-- Compute reference density profile
DO k=1,Nr
pLoc = rC(k)
CALL FIND_RHO_SCALAR(
I tRef(k), sRef(k), pLoc,
O rhoRef(k), myThid )
ENDDO
C-- Compute reference stratification: -d.alpha/dp @ constant p
dBdrRef(1) = 0. _d 0
DO k=1,Nr+1
pLoc = rF(k)
IF ( k.GE.2 ) CALL FIND_RHO_SCALAR(
I tRef(k-1), sRef(k-1), pLoc,
O rhoDw, myThid )
IF ( k.LE.Nr ) CALL FIND_RHO_SCALAR(
I tRef(k), sRef(k), pLoc,
O rhoUp, myThid )
IF ( k.GE.2 .AND. k.LE.Nr ) THEN
dBdrRef(k) = (rhoDw - rhoUp)*recip_drC(k)
& / (rhoDw*rhoUp)
rhoLoc = ( rhoDw + rhoUp )*0.5 _d 0
ELSEIF ( k.EQ.1 ) THEN
rhoLoc = rhoUp
ELSE
rhoLoc = rhoDw
ENDIF
C-- Units convertion factor for vertical velocity:
C wUnit2rVel = gravity*rhoRef : rVel [Pa/s] = wSpeed [m/s] * wUnit2rVel
C rVel2wUnit = 1/rVel2wUnit : wSpeed [m/s] = rVel [Pa/s] * rVel2wUnit
C note: wUnit2rVel & rVel2wUnit replace horiVertRatio & recip_horiVertRatio
wUnit2rVel(k) = gravity*rhoLoc
rVel2wUnit(k) = 1. _d 0 / wUnit2rVel(k)
ENDDO
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
ELSEIF ( buoyancyRelation.EQ.'ATMOSPHERIC' ) THEN
C-- Compute reference stratification: -d.alpha/dp @ constant p
dBdrRef(1) = 0. _d 0
DO k=2,Nr
conv_theta2T = (rF(k)/atm_Po)**atm_kappa
c dBdrRef(k) = (tRef(k) - tRef(k-1))*recip_drC(k)
c & * conv_theta2T*atm_Rd/rF(k)
ddPI=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
& -((rC( k )/atm_Po)**atm_kappa) )
dBdrRef(k) = (tRef(k) - tRef(k-1))*recip_drC(k)
& * ddPI*recip_drC(k)
ENDDO
C-- Units convertion factor for vertical velocity:
C wUnit2rVel = gravity/alpha : rVel [Pa/s] = wSpeed [m/s] * wUnit2rVel
C rVel2wUnit = alpha/gravity : wSpeed [m/s] = rVel [Pa/s] * rVel2wUnit
C with alpha = 1/rhoRef = (R.T/p) (ideal gas)
C note: wUnit2rVel & rVel2wUnit replace horiVertRatio & recip_horiVertRatio
DO k=1,Nr+1
IF ( k.EQ.1 ) THEN
thetaLoc = tRef(k)
ELSEIF ( k.GT.Nr ) THEN
thetaLoc = tRef(k-1)
ELSE
thetaLoc = (tRef(k) + tRef(k-1))*0.5 _d 0
ENDIF
IF ( thetaLoc.GT.0. _d 0 .AND. rF(k).GT.0. _d 0 ) THEN
conv_theta2T = (rF(k)/atm_Po)**atm_kappa
wUnit2rVel(k) = gravity
& * rF(k)/(atm_Rd*conv_theta2T*thetaLoc)
rVel2wUnit(k) = 1. _d 0 / wUnit2rVel(k)
ENDIF
ENDDO
C- Compute Reference Geopotential at Half levels :
C Tracer level k: phiRef(2k) ; Interface_W level k: phiRef(2k-1)
phiRef(1) = seaLev_Z*gravity
IF ( select_rStar.GE.1 .OR. selectSigmaCoord.GE.1 ) THEN
C- isothermal (theta=const) reference state
DO k=1,Nr
tLoc(k) = thetaConst
ENDDO
ELSE
C- horizontally uniform (tRef) reference state
DO k=1,Nr
tLoc(k) = tRef(k)
ENDDO
ENDIF
IF (integr_GeoPot.EQ.1) THEN
C- Finite Volume Form, linear by half level :
DO k=1,2*Nr
ks = (k+1)/2
ddPI=atm_Cp*( ((rHalf( k )/atm_Po)**atm_kappa)
& -((rHalf(k+1)/atm_Po)**atm_kappa) )
phiRef(k+1) = phiRef(k)+ddPI*tLoc(ks)
ENDDO
C------
ELSE
C- Finite Difference Form, linear between Tracer level :
C works with integr_GeoPot = 0, 2 or 3
k = 1
ddPI=atm_Cp*( ((rF(k)/atm_Po)**atm_kappa)
& -((rC(k)/atm_Po)**atm_kappa) )
phiRef(2*k) = phiRef(1) + ddPI*tLoc(k)
DO k=1,Nr-1
ddPI=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
& -((rC(k+1)/atm_Po)**atm_kappa) )
phiRef(2*k+1) = phiRef(2*k) + ddPI*0.5*tLoc(k)
phiRef(2*k+2) = phiRef(2*k)
& + ddPI*0.5*(tLoc(k)+tLoc(k+1))
ENDDO
k = Nr
ddPI=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
& -((rF(k+1)/atm_Po)**atm_kappa) )
phiRef(2*k+1) = phiRef(2*k) + ddPI*tLoc(k)
C------
ENDIF
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
ELSE
STOP 'SET_REF_STATE: Bad value of buoyancyRelation !'
C-- endif buoyancyRelation
ENDIF
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
IF ( usingZCoords .AND. rhoRefFile .NE. ' ' ) THEN
C-- anelastic formulation : set density factor from reference density profile
C surface-interface rho-factor has to be 1:
rhoFacF(1) = 1. _d 0
C rhoFac(k) = density ratio between layer k and top interface
DO k=1,Nr
rhoFacC(k) = rho1Ref(k)/rhoConst
c rhoFacC(k) = rho1Ref(k)*recip_rhoConst
ENDDO
DO k=2,Nr
C since rkSign=-1, recip_drC(k) = 1./(rC(k-1)-rC(k))
rhoFacF(k) = ( rhoFacC(k-1)*(rF(k)-rC(k))
& + rhoFacC(k)*(rC(k-1)-rF(k)) )*recip_drC(k)
ENDDO
C extrapolate down to the bottom:
rhoFacF(Nr+1) = ( rhoFacC(Nr)*(rF(Nr+1)-rF(Nr))
& + rhoFacF(Nr)*(rC(Nr)-rF(Nr+1))
& ) / (rC(Nr)-rF(Nr))
C- set reciprocal rho-factor:
DO k=1,Nr
recip_rhoFacC(k) = 1. _d 0/rhoFacC(k)
ENDDO
DO k=1,Nr+1
recip_rhoFacF(k) = 1. _d 0/rhoFacF(k)
ENDDO
ENDIF
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C-- Write to check :
IF ( gravityFile .NE. ' ' ) THEN
C- write gravity vertical profile factor to binary file:
CALL WRITE_GLVEC_RL( 'GravFacC',' ',gravFacC, Nr, -1, myThid )
CALL WRITE_GLVEC_RL( 'GravFacF',' ',gravFacF,Nr+1,-1, myThid )
ENDIF
IF ( selectP_inEOS_Zc.GE.1 ) THEN
C- write (to binary file) Hyd-Pressure used in EOS:
CALL WRITE_GLVEC_RL( 'PRef4EOS',' ',pRef4EOS, Nr, -1, myThid )
c CALL WRITE_GLVEC_RL( 'PRefLocF',' ',pRefLocF,Nr+1,-1, myThid )
ENDIF
IF ( buoyancyRelation.EQ.'ATMOSPHERIC' ) THEN
WRITE(msgBuf,'(A)') ' '
CALL PRINT_MESSAGE( msgBuf, stdUnit, SQUEEZE_RIGHT, myThid )
WRITE(msgBuf,'(A)')
& 'SET_REF_STATE: PhiRef/g [m] at level Center (integer)'
CALL PRINT_MESSAGE( msgBuf, stdUnit, SQUEEZE_RIGHT, myThid )
WRITE(msgBuf,'(A)')
& ' and at level Interface (half-int.) :'
CALL PRINT_MESSAGE( msgBuf, stdUnit, SQUEEZE_RIGHT, myThid )
DO k=1,2*Nr+1
WRITE(msgBuf,'(A,F5.1,A,F15.1,A,F13.3)')
& ' K=',k*0.5,' ; r=',rHalf(k),' ; phiRef/g=',
& phiRef(k)*recip_gravity
CALL PRINT_MESSAGE(msgBuf, stdUnit, SQUEEZE_RIGHT, myThid )
ENDDO
ELSE
C- Write reference density to binary file :
CALL WRITE_GLVEC_RL( 'RhoRef', ' ', rhoRef, Nr, -1, myThid )
ENDIF
IF ( usingZCoords .AND. rhoRefFile .NE. ' ' ) THEN
C- Write Anelastic density factor to binary file :
CALL WRITE_GLVEC_RL( 'RhoFacC',' ',rhoFacC, Nr , -1, myThid )
CALL WRITE_GLVEC_RL( 'RhoFacF',' ',rhoFacF, Nr+1, -1, myThid )
ENDIF
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
_END_MASTER( myThid )
_BARRIER
RETURN
END