C $Header: /home/ubuntu/mnt/e9_copy/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