C $Header: /u/gcmpack/MITgcm/model/src/ini_linear_phisurf.F,v 1.24 2014/05/12 01:27:47 jmc Exp $ C $Name: $ #include "PACKAGES_CONFIG.h" #include "CPP_OPTIONS.h" CBOP C !ROUTINE: INI_LINEAR_PHISURF C !INTERFACE: SUBROUTINE INI_LINEAR_PHISURF( myThid ) C !DESCRIPTION: \bv C *==========================================================* C | SUBROUTINE INI_LINEAR_PHISURF C | o Initialise the Linear Relation Phi_surf(eta) C *==========================================================* C | Initialise -Buoyancy at surface level (Bo_surf) C | to setup the Linear relation: Phi_surf(eta)=Bo_surf*eta C | Initialise phi0surf = starting point for integrating C | phiHyd (= phiHyd at r=RoSurf) C *==========================================================* C \ev C !USES: IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "SURFACE.h" C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C myThid :: my Thread Id number INTEGER myThid C == Local variables in common == C topoHloc had to be in common for multi threading but no longer C needed since MDSIO now allows (2009/06/07) to write local arrays C !LOCAL VARIABLES: C === Local variables === C topoHloc :: Temporary array used to write surface topography C bi,bj :: tile indices C i,j,k :: Loop counters _RS topoHloc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) INTEGER bi, bj INTEGER i, j, k _RL pLoc, rhoLoc _RL dPIdp CEOP C-- Initialisation #ifdef ALLOW_AUTODIFF DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx Bo_surf(i,j,bi,bj) = 0. _d 0 recip_Bo(i,j,bi,bj) = 0. _d 0 ENDDO ENDDO ENDDO ENDDO #endif /* ALLOW_AUTODIFF */ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C-- Initialise -Buoyancy at surface level : Bo_surf C Bo_surf is defined as d/dr(Phi_surf) and set to g/z2rUnit with C z2rUnit = conversion factor from z-unit to r-unit: [z] * z2rUnit = [r] C an accurate formulation includes P_surf and T,S_ref effects on rho_surf: C (setting uniformLin_PhiSurf=.FALSE.): C z-coord (z2rUnit=1): Bo_surf = - Buoyancy C = g * rho_surf(Tref,Sref,pSurf_0)/rho_0 C p-coord (z2rUnit=rho*g): Bo_surf = 1/rho(Tref(ksurf),pSurf_0) C Note: on Phi_surf splitting : Non-linear Time-dependent effects on B_surf C [through eta & (T-tRef)_surf] are included in PhiHyd rather than in Bo_surf C-- IF ( usingZCoords ) THEN C- gBaro = gravity (except for External mode test with reduced gravity) DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx Bo_surf(i,j,bi,bj) = gBaro recip_Bo(i,j,bi,bj) = 1. _d 0 / gBaro ENDDO ENDDO ENDDO ENDDO ELSEIF ( uniformLin_PhiSurf ) THEN C- use a linear (in ps) uniform relation : Phi'_surf = 1/rhoConst * ps'_surf DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx c Bo_surf(i,j,bi,bj) = rVel2wUnit(1)*gravity c recip_Bo(i,j,bi,bj) = wUnit2rVel(1)*recip_gravity Bo_surf(i,j,bi,bj) = recip_rhoConst recip_Bo(i,j,bi,bj) = rhoConst ENDDO ENDDO ENDDO ENDDO ELSEIF ( fluidIsWater ) THEN C-- More precise than uniformLin_PhiSurf case but inconsistent C with nonlinFreeSurf=4 in CALC_PHI_HYD (eta contribution to phiHyd) DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx IF ( Ro_surf(i,j,bi,bj).GT.0. _d 0 & .AND. kSurfC(i,j,bi,bj).LE.Nr ) THEN pLoc = Ro_surf(i,j,bi,bj) #ifdef ALLOW_OPENAD CALL FIND_RHO_SCALAR( I tRef(kSurfC(i,j,bi,bj)), I sRef(kSurfC(i,j,bi,bj)), I pLoc, O rhoLoc, myThid ) #else /* ALLOW_OPENAD */ k = kSurfC(i,j,bi,bj) CALL FIND_RHO_SCALAR( I tRef(k), sRef(k), pLoc, O rhoLoc, myThid ) #endif /* ALLOW_OPENAD */ IF ( rhoLoc .EQ. 0. _d 0 ) THEN Bo_surf(i,j,bi,bj) = 0. _d 0 ELSE Bo_surf(i,j,bi,bj) = 1. _d 0/rhoLoc ENDIF recip_Bo(i,j,bi,bj) = rhoLoc ELSE Bo_surf(i,j,bi,bj) = 0. _d 0 recip_Bo(i,j,bi,bj) = 0. _d 0 ENDIF ENDDO ENDDO ENDDO ENDDO ELSEIF ( fluidIsAir ) THEN C- use a linearized (in ps) Non-uniform relation : Bo_surf(Po_surf,tRef_surf) C Bo = d/d_p(Phi_surf) = tRef_surf*d/d_p(PI) ; PI = Cp*(p/Po)^kappa C and atm_Cp*atm_kappa = atm_Rd DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) IF ( select_rStar.GE.1 .OR. selectSigmaCoord.GE.1 ) THEN C- isothermal (theta=const) reference state DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx IF ( Ro_surf(i,j,bi,bj).GT.0. _d 0 & .AND. kSurfC(i,j,bi,bj).LE.Nr ) THEN dPIdp = (atm_Rd/atm_Po)* & (Ro_surf(i,j,bi,bj)/atm_Po)**(atm_kappa-1. _d 0) Bo_surf(i,j,bi,bj) = dPIdp*thetaConst recip_Bo(i,j,bi,bj) = 1. _d 0 / Bo_surf(i,j,bi,bj) ELSE Bo_surf(i,j,bi,bj) = 0. recip_Bo(i,j,bi,bj) = 0. ENDIF ENDDO ENDDO ELSE C- horizontally uniform (tRef) reference state DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx IF ( Ro_surf(i,j,bi,bj).GT.0. _d 0 & .AND. kSurfC(i,j,bi,bj).LE.Nr ) THEN dPIdp = (atm_Rd/atm_Po)* & (Ro_surf(i,j,bi,bj)/atm_Po)**(atm_kappa-1. _d 0) Bo_surf(i,j,bi,bj) = dPIdp*tRef(kSurfC(i,j,bi,bj)) recip_Bo(i,j,bi,bj) = 1. _d 0 / Bo_surf(i,j,bi,bj) ELSE Bo_surf(i,j,bi,bj) = 0. recip_Bo(i,j,bi,bj) = 0. ENDIF ENDDO ENDDO ENDIF ENDDO ENDDO ELSE STOP 'INI_LINEAR_PHISURF: We should never reach this point!' ENDIF C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C-- Update overlap regions (jmc: is it really needed ?) c _EXCH_XY_RL(Bo_surf, myThid) c _EXCH_XY_RL(recip_Bo, myThid) IF ( usingPCoords .AND. .NOT.uniformLin_PhiSurf ) THEN CALL WRITE_FLD_XY_RL( 'Bo_surf',' ',Bo_surf,0,myThid) ENDIF C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C-- Initialise phi0surf: used for atmos. surf. P-loading (ocean, z-coord) C or topographic geopotential anom. (p-coord) DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx phi0surf(i,j,bi,bj) = 0. ENDDO ENDDO ENDDO ENDDO IF ( fluidIsAir .AND. topoFile.NE.' ' ) THEN #ifdef ALLOW_AUTODIFF STOP 'CANNOT PRESENTLY USE THIS OPTION WITH ADJOINT' #else C-- Compute topoH = PhiRef(Po_surf)/g ; is different from original C topoZ(read from file) because of truncation of Po_surf. C NOTE: not clear for now which topoZ needs to be saved in common block C-- AND set phi0surf = starting point for integrating Geopotential; CALL INI_P_GROUND( -2, O topoHloc, I Ro_surf, myThid ) IF (selectFindRoSurf.NE.0) THEN _EXCH_XY_RS(phi0surf, myThid) CALL WRITE_FLD_XY_RS( 'phi0surf',' ',phi0surf,0,myThid) ENDIF CALL WRITE_FLD_XY_RS( 'topo_H',' ',topoHloc,0,myThid) #endif /* ALLOW_AUTODIFF */ ENDIF C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| RETURN END