C $Header: /u/gcmpack/MITgcm/pkg/ggl90/ggl90_init_varia.F,v 1.8 2015/02/21 20:10:08 jmc Exp $
C $Name:  $

#include "GGL90_OPTIONS.h"

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

C     !DESCRIPTION: \bv
C     *================================================================*
C     | SUBROUTINE GGL90_INIT_VARIA
C     | o initialize variables
C     *================================================================*
C     \ev

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

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     myThid -  Number of this instance of GGL90_INIT
      INTEGER myThid

#ifdef ALLOW_GGL90
C     !LOCAL VARIABLES:
C     == Local variables ==
C     bi,bj,I,J,K,Km1 - Loop counters
      INTEGER bi, bj
      INTEGER I, J, K
      INTEGER Km1
CEOP

C--   Over all tiles
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)

C-      3D arrays
        DO K=1,Nr
         Km1 = MAX(K-1,1)
         DO J=1-OLy,sNy+OLy
          DO I=1-OLx,sNx+OLx
           GGL90viscArU(I,J,K,bi,bj) = 0. _d 0
           GGL90viscArV(I,J,K,bi,bj) = 0. _d 0
           GGL90diffKr(I,J,K,bi,bj)  = 0. _d 0
           IF ( useIDEMIX) THEN
            GGL90TKE(I,J,K,bi,bj)=GGL90eps*maskC(I,J,K,bi,bj)
           ELSE
CML           GGL90TKE(I,J,K,bi,bj)=0. _d 0
            GGL90TKE(I,J,K,bi,bj)=GGL90TKEmin*maskC(I,J,K,bi,bj)
           ENDIF
          ENDDO
         ENDDO
        ENDDO

C-      2D arrays

       ENDDO
      ENDDO

#ifdef ALLOW_GGL90_IDEMIX
c-----------------------------------------------------------------------
c  IDEMIX
c-----------------------------------------------------------------------
C--   Over all tiles
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)

        DO K=1,Nr
         DO J=1-OLy,sNy+OLy
          DO I=1-OLx,sNx+OLx
           IDEMIX_E(I,J,K,bi,bj)     = 0. _d 0
           IDEMIX_V0(I,J,K,bi,bj)    = 0. _d 0
           IDEMIX_tau_d(I,J,K,bi,bj) = 0. _d 0
          ENDDO
         ENDDO
        ENDDO

        DO J=1-OLy,sNy+OLy
         DO I=1-OLx,sNx+OLx
          IDEMIX_F_B(I,J,bi,bj) = 0. _d 0
          IDEMIX_F_S(I,J,bi,bj) = 0. _d 0
         ENDDO
        ENDDO

       ENDDO
      ENDDO
c-----------------------------------------------------------------------
c      read in forcing fields for IDEMIX
c-----------------------------------------------------------------------
      IF ( useIDEMIX .AND. IDEMIX_tidal_file .NE. ' ' ) THEN
       CALL READ_REC_XY_RL( IDEMIX_tidal_file,IDEMIX_F_b,1,0,myThid )
       _EXCH_XY_RL(IDEMIX_F_b,myThid)

       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO J=1-OLy,sNy+OLy
          DO I=1-OLx,sNx+OLx
           IDEMIX_F_b(i,j,bi,bj) = -MAX(0. _d 0,
     &          MIN(1. _d 0,IDEMIX_F_b(i,j,bi,bj)))
C     fraction that enters IW field
           IDEMIX_F_b(i,j,bi,bj) = IDEMIX_frac_F_b*
CML     &          IDEMIX_F_b(i,j,bi,bj)*recip_rhoConst
     &          IDEMIX_F_b(i,j,bi,bj)/1024. _d 0
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDIF

      IF ( useIDEMIX .AND. IDEMIX_wind_file .NE. ' ' ) THEN
       CALL READ_REC_XY_RL( IDEMIX_wind_file ,IDEMIX_F_s,1,0,myThid )
       _EXCH_XY_RL(IDEMIX_F_s,myThid)
       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO J=1-OLy,sNy+OLy
          DO I=1-OLx,sNx+OLx
           IDEMIX_F_s(i,j,bi,bj) = MAX(0. _d 0,
     &          MIN(1. _d 0,IDEMIX_F_s(i,j,bi,bj)))
C     fraction that leaves mixed layer
           IDEMIX_F_s(i,j,bi,bj) = IDEMIX_frac_F_s*
CML     &          IDEMIX_F_s(i,j,bi,bj)*recip_rhoConst
     &          IDEMIX_F_s(i,j,bi,bj)/1024.
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDIF
#endif /* ALLOW_GGL90_IDEMIX */

      IF ( nIter0.NE.0 .OR. pickupSuff.NE.' ' ) THEN
       CALL GGL90_READ_PICKUP( nIter0, myThid )
      ELSE
C--   Read initial turbulent kinetic energy thickness from file if available.
       IF ( GGL90TKEFile .NE. ' ' ) THEN
        CALL READ_FLD_XYZ_RL( GGL90TKEFile, ' ', GGL90TKE, 0, myThid )
        _EXCH_XYZ_RL(GGL90TKE,myThid)
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO K=1,Nr
           DO J=1-OLy,sNy+OLy
            DO I=1-OLx,sNx+OLx
             GGL90TKE(I,J,K,bi,bj) = MAX(GGL90TKE(I,J,K,bi,bj),
     &            GGL90TKEmin)*maskC(I,J,K,bi,bj)
            ENDDO
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ENDIF
      ENDIF

#endif /*  ALLOW_GGL90  */

      RETURN
      END