C $Header: /u/gcmpack/MITgcm/model/src/ini_nh_fields.F,v 1.6 2013/08/11 20:23:39 jmc Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: INI_NH_FIELDS
C     !INTERFACE:
      SUBROUTINE INI_NH_FIELDS( myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE INI_NH_FIELDS
C     | o Set model initial non-hydrostatic fields.
C     *==========================================================*
C     | Note: If using NH form,
C     |  call this S/R whether starting or restarting simulation.
C     | This is different from other "true" ini_fields type S/R
C     |  (e.g., INI_VEL) which are called only when starting.
C     | Reason: no real physical field to initialise (since wVel
C     |  is diagnose from continuity) but needs to set few arrays
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "RESTART.h"
#include "NH_VARS.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     myThid  :: My Thread Id number
      INTEGER myThid

#ifdef ALLOW_NONHYDROSTATIC
C     !LOCAL VARIABLES:
C     == Local variables ==
      INTEGER bi,bj
      INTEGER i,j
      INTEGER ks
      CHARACTER*(MAX_LEN_MBUF) msgBuf
#ifdef NONLIN_FRSURF
      INTEGER k
#endif
CEOP

      IF ( startTime .EQ. baseTime .AND.  nIter0 .EQ. 0
     &     .AND. pickupSuff .EQ. ' ' ) THEN
C--   Case where starting from initial conditions

C--   Read an initial non-hydrostatic pressure field
c       IF (phiNHinitFile .NE. ' ') THEN
c        CALL READ_FLD_XYZ_RL( phiNHinitFile, ' ', phi_nh, 0, myThid )
c        _EXCH_XYZ_RL(phi_nh, myThid)
c       ENDIF

      ELSE
C--   Case where restarting from a pickup

       _BEGIN_MASTER(myThid)
       WRITE(msgBuf,'(A,I4)')
     &   'INI_NH_FIELDS: dPhiNHstatus=', dPhiNHstatus
       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                     SQUEEZE_RIGHT, myThid )
       _END_MASTER(myThid)
       IF ( exactConserv .AND. dPhiNHstatus.EQ.0 ) THEN
c      IF ( exactConserv ) THEN
C--   Separate the Hydrostatic Surface Pressure adjusment (=> put it in dPhiNH)
C     from the Non-hydrostatic pressure (since cg3d_x contains both contribution)
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          IF ( select_rStar.EQ.0 .AND. uniformFreeSurfLev ) THEN
C-       Z coordinate: assume surface @ level k=1
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
              dPhiNH(i,j,bi,bj) = phi_nh(i,j,1,bi,bj)
c             dPhiNH(i,j,bi,bj) = 0.
            ENDDO
           ENDDO
          ELSEIF ( select_rStar.EQ.0 ) THEN
C-       Other than Z coordinate: no assumption on surface level index
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             ks = kSurfC(i,j,bi,bj)
             IF ( ks.LE.Nr ) THEN
              dPhiNH(i,j,bi,bj) = phi_nh(i,j,ks,bi,bj)
             ELSE
              dPhiNH(i,j,bi,bj) = 0.
             ENDIF
            ENDDO
           ENDDO
#ifdef NONLIN_FRSURF
          ELSE
C        rStar : take vertical average of P_NH as Hyd.Surf.Press adjustment
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
              dPhiNH(i,j,bi,bj) = 0.
            ENDDO
           ENDDO
           DO k=1,Nr
            DO j=1-OLy,sNy+OLy
             DO i=1-OLx,sNx+OLx
              dPhiNH(i,j,bi,bj) = dPhiNH(i,j,bi,bj)
     &          + phi_nh(i,j,k,bi,bj)*drF(k)*hFacC(i,j,k,bi,bj)
             ENDDO
            ENDDO
           ENDDO
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
              dPhiNH(i,j,bi,bj) = dPhiNH(i,j,bi,bj)
     &                           *recip_Rcol(i,j,bi,bj)
            ENDDO
           ENDDO
#endif /* NONLIN_FRSURF */
          ENDIF
         ENDDO
        ENDDO
C-     end of if-block: dPhiNH_status
       ENDIF

      ENDIF

#endif /* ALLOW_NONHYDROSTATIC */
      RETURN
      END