C $Header: /u/gcmpack/MITgcm/pkg/aim_v23/aim_dyn2aim.F,v 1.6 2014/05/12 01:32:41 jmc Exp $
C $Name:  $

#include "AIM_OPTIONS.h"

CBOP
C     !ROUTINE: AIM_DYN2AIM
C     !INTERFACE:
      SUBROUTINE AIM_DYN2AIM(
     O           TA, QA, ThA, Vsurf2, PSA, dpFac,
     O           kGrd,
     I           bi,bj, myTime, myIter, myThid)
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R AIM_DYN2AIM
C     | o Map dynamics conforming arrays to AIM internal arrays.
C     *==========================================================*
C     | this routine transfers grid information
C     | and all dynamical variables (T,Q, ...) to AIM physics
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
C-- size for MITgcm & Physics package :
#include "AIM_SIZE.h"

#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"
#include "DYNVARS.h"

#include "AIM_GRID.h"
#include "com_physcon.h"

C     !INPUT/OUTPUT PARAMETERS:
C--   Input:
C     bi,bj  :: Tile index
C     myTime :: Current time of simulation ( s )
C     myIter :: Current iteration number in simulation
C     myThid :: Number of this instance of the routine
C--   Output:  TA     :: temperature  [K}                        (3-dim)
C              QA     :: specific humidity [g/kg]                (3-dim)
C              ThA    :: Pot.Temp. [K] (replace dry stat. energy)(3-dim)
C              Vsurf2 :: square of surface wind speed            (2-dim)
C              PSA    :: norm. surface pressure [p/p0]           (2-dim)
C              dpFac  :: cell delta_P fraction                   (3-dim)
C              kGrd   :: Ground level index                      (2-dim)
C--  Updated common blocks: AIM_GRID_R
C             WVSurf  :: weights for near surf interpolation     (2-dim)
C             fOrogr  :: orographic factor (for surface drag)    (2-dim)
C         snLat,csLat :: sin(Lat) & cos(Lat)                     (2-dim)
      _RL TA(NGP,NLEV), QA(NGP,NLEV), ThA(NGP,NLEV)
      _RL Vsurf2(NGP), PSA(NGP), dpFac(NGP,NLEV)
      INTEGER kGrd(NGP)
      INTEGER bi, bj, myIter, myThid
      _RL myTime
CEOP

#ifdef ALLOW_AIM
C     !LOCAL VARIABLES:
C     i, j, I2 :: Loop counters
C     k, Katm  :: Loop counters
C     msgBuf   :: Informational/error message buffer
      CHARACTER*(MAX_LEN_MBUF) msgBuf
      INTEGER i, j, I2, k, Katm
      _RL conv_theta2T, temp1, temp2

c     _RL hInitC(5), hInitF(5)
c     DATA    hInitC / 17338.0,10090.02,5296.88,2038.54,418.038/
c     DATA    hInitF / 15090.4, 8050.96, 4087.75, 1657.54, 0. /

C-    Compute Sin & Cos (Latitude) :
      DO j = 1,sNy
       DO i = 1,sNx
        I2 = i+(j-1)*sNx
        snLat(I2,myThid) = SIN(yC(i,j,bi,bj)*deg2rad)
        csLat(I2,myThid) = COS(yC(i,j,bi,bj)*deg2rad)
       ENDDO
      ENDDO

C-    Set surface level index :
      DO j = 1,sNy
       DO i = 1,sNx
        I2 = i+(j-1)*sNx
        kGrd(I2) = (Nr+1) - kSurfC(i,j,bi,bj)
       ENDDO
      ENDDO

C-    Set (normalized) surface pressure :
      DO j=1,sNy
       DO i=1,sNx
        I2 = i+(j-1)*sNx
        k = kSurfC(i,j,bi,bj)
        IF ( k.LE.Nr ) THEN
c        PSA(I2) = rF(k)/atm_Po
         PSA(I2) = Ro_surf(i,j,bi,bj)/atm_Po
        ELSE
         PSA(I2) = 1.
        ENDIF
       ENDDO
      ENDDO

C-    Set cell delta_P fraction (of the full delta.P = drF_k):
#ifdef NONLIN_FRSURF
      IF ( staggerTimeStep .AND. nonlinFreeSurf.GT.0 ) THEN
       IF ( select_rStar.GT.0 ) THEN
        DO k = 1,Nr
         Katm = _KD2KA( k )
         DO j = 1,sNy
          DO i = 1,sNx
           I2 = i+(j-1)*sNx
           dpFac(I2,Katm) = h0FacC(i,j,k,bi,bj)*rStarFacC(i,j,bi,bj)
c          dpFac(I2,Katm) = 1. _d 0
          ENDDO
         ENDDO
        ENDDO
       ELSE
        DO k = 1,Nr
         Katm = _KD2KA( k )
         DO j = 1,sNy
          DO i = 1,sNx
           I2 = i+(j-1)*sNx
           IF ( k.EQ.kSurfC(i,j,bi,bj) ) THEN
            dpFac(I2,Katm) = hFac_surfC(i,j,bi,bj)
           ELSE
            dpFac(I2,Katm) = hFacC(i,j,k,bi,bj)
           ENDIF
c          dpFac(I2,Katm) = 1. _d 0
          ENDDO
         ENDDO
        ENDDO
       ENDIF
      ELSE
#else /* ndef NONLIN_FRSURF */
      IF (.TRUE.) THEN
#endif /* NONLIN_FRSURF */
        DO k = 1,Nr
         Katm = _KD2KA( k )
         DO j = 1,sNy
          DO i = 1,sNx
           I2 = i+(j-1)*sNx
           dpFac(I2,Katm) = hFacC(i,j,k,bi,bj)
c          dpFac(I2,Katm) = 1. _d 0
          ENDDO
         ENDDO
        ENDDO
      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C     Physics package works with sub-domains 1:sNx,1:sNy,1:Nr.
C     Internal index mapping is linear in X and Y with a second
C     dimension for the vertical.

C-    Dynamical var --> AIM var :
C       note: UA & VA are not used  => removed
      temp1 = lwTemp1
      temp2 = lwTemp2
      DO k = 1,Nr
       conv_theta2T = (rC(k)/atm_Po)**atm_kappa
       Katm = _KD2KA( k )
       DO j = 1,sNy
        DO i = 1,sNx
         I2 = i+(j-1)*sNx
         IF ( maskC(i,j,k,bi,bj).EQ.oneRS ) THEN
c         UA(I2,Katm)  = uVel(i,j,k,bi,bj)
c         VA(I2,Katm)  = vVel(i,j,k,bi,bj)
C     Physics works with temperature - not potential temp.
          TA(I2,Katm)  = theta(i,j,k,bi,bj)*conv_theta2T
c         TA(I2,Katm)  = max(temp1,min(temp2,
c    &                       theta(i,j,k,bi,bj)*conv_theta2T ))
C     In atm.Phys, water vapor must be > 0 :
          QA(I2,Katm)  = MAX( salt(i,j,k,bi,bj), zeroRL )
C     Dry static energy replaced by Pot.Temp:
          ThA(I2,Katm) = theta(i,j,k,bi,bj)
         ELSE
          TA(I2,Katm)  = 300. _d 0
          QA(I2,Katm)  =   0. _d 0
          ThA(I2,Katm) = 300. _d 0
         ENDIF
        ENDDO
       ENDDO
#ifdef NONLIN_FRSURF
       IF ( select_rStar.GE.1 ) THEN
        DO j = 1,sNy
         DO i = 1,sNx
          I2 = i+(j-1)*sNx
          TA(I2,Katm)  = TA(I2,Katm)*pStarFacK(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDIF
#endif /* NONLIN_FRSURF */
      ENDDO

C_jmc: add square of surface wind speed (center of C grid) = 2 * KE_surf
      DO j = 1,sNy
       DO i = 1,sNx
        I2 = i+(j-1)*sNx
        k = kSurfC(i,j,bi,bj)
        IF (k.LE.Nr) THEN
         Vsurf2(I2) = 0.5 * (
     &                uVel(i,j,k,bi,bj)*uVel(i,j,k,bi,bj)
     &              + uVel(i+1,j,k,bi,bj)*uVel(i+1,j,k,bi,bj)
     &              + vVel(i,j,k,bi,bj)*vVel(i,j,k,bi,bj)
     &              + vVel(i,j+1,k,bi,bj)*vVel(i,j+1,k,bi,bj)
     &                        )
        ELSE
         Vsurf2(I2) = 0.
        ENDIF
       ENDDO
      ENDDO

C-    Check that Temp is OK for LW Radiation scheme :
      DO k = 1,Nr
       Katm = _KD2KA( k )
       DO I2=1,NGP
        IF (  TA(I2,Katm).LT.lwTemp1 .OR.
     &        TA(I2,Katm).GT.lwTemp2 ) THEN
         i = 1 + mod((I2-1),sNx)
         j = 1 + int((I2-1)/sNx)
         WRITE(msgBuf,'(A,1PE20.13,A,2I4)')
     &    'AIM_DYN2AIM: Temp=', TA(I2,Katm),
     &    ' out of range ',lwTemp1,lwTemp2
         CALL PRINT_ERROR( msgBuf , myThid)
         WRITE(msgBuf,'(A,3I4,3I3,I6,2F9.3)')
     &    'AIM_DYN2AIM: Pb in i,j,k,bi,bj,myThid,I2,X,Y=',
     &        i,j,k,bi,bj,myThid,I2,xC(i,j,bi,bj),yC(i,j,bi,bj)
         CALL PRINT_ERROR( msgBuf , myThid)
         STOP 'ABNORMAL END: S/R AIM_DYN2AIM'
        ENDIF
       ENDDO
      ENDDO

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C-    Set geopotential surfaces
c     DO Katm=1,NLEV
c      DO I2=1,NGP
c       PHIG1(I2,Katm) = gravity*HinitC(Katm)
c      ENDDO
c     ENDDO

C-    Weights for vertical interpolation down to the surface
C     Fsurf = Ffull(nlev)+WVS*(Ffull(nlev)-Ffull(nlev-1))
      DO j = 1,sNy
       DO i = 1,sNx
        I2 = i+(j-1)*sNx
        WVSurf(I2,myThid) = 0.
        k = kGrd(I2)
        IF (k.GT.1) THEN
C- full cell version of Franco Molteni formula:
c         WVSurf(I2,myThid) = (LOG(SIGH(k))-SIGL(k))*WVI(k-1,2)
C- partial cell version using true log-P extrapolation:
          WVSurf(I2,myThid) = (LOG(PSA(I2))-SIGL(k))*WVI(k-1,1)
C- like in the old code:
c         WVSurf(I2,myThid) = WVI(k,2)
        ENDIF
       ENDDO
      ENDDO
      IF (myIter.EQ.nIter0)
     &  CALL AIM_WRITE_PHYS( 'aim_WeightSurf', '', 1, WVSurf,
     &                       1, bi, bj, 1, myIter, myThid )

#endif /* ALLOW_AIM */

      RETURN
      END