C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/aim_v23/phy_vdifsc.F,v 1.1 2002/11/22 17:17:03 jmc Exp $
C $Name: checkpoint58b_post $

#include "AIM_OPTIONS.h"

      SUBROUTINE VDIFSC (dpFac,SE,RH,QA,QSAT,
     O                   TTENVD,QTENVD,
     I                   kGrd,bi,bj,myThid)
C--
C--   SUBROUTINE VDIFSC (UA,VA,SE,RH,QA,QSAT,PHI,
C--  &                   UTENVD,VTENVD,TTENVD,QTENVD)
C--
C--   Purpose: Compute tendencies of momentum, energy and moisture
C--            due to vertical diffusion and shallow convection
C--   Input:   UA     = u-wind                           (3-dim)
C--            VA     = v-wind                           (3-dim)
C              dpFac  = cell delta_P fraction            (3-dim)
C--            SE     = dry static energy                (3-dim)
C--            RH     = relative humidity [0-1]          (3-dim)
C--            QA     = specific humidity [g/kg]         (3-dim)
C--            QSAT   = saturation sp. humidity [g/kg]   (3-dim)
C--            PHI    = geopotential                     (3-dim)
C--   Output:  UTENVD = u-wind tendency                  (3-dim)
C--            VTENVD = v-wind tendency                  (3-dim)
C--            TTENVD = temperature tendency             (3-dim)
C--            QTENVD = sp. humidity tendency [g/(kg s)] (3-dim)
C    Input:    kGrd   = Ground level index               (2-dim)
C              bi,bj  = tile index
C           myThid    = Thread number for this instance of the routine
C-------
C  Note: a) dry static energy (SE,input) has been replaced by Pot.Temp.
C        b) UA,VA & U,V_TENVD not used => removed 
C-------
C In contrast to other Physics S/R, VDIFSC return a real tendency (dQ/dt,dT/dt)
C  nevertheless /dpFac is not applied here but later in AIM_AIM2DYN
C-------

      IMPLICIT NONE

C     Resolution parameters

C-- size for MITgcm & Physics package :
#include "AIM_SIZE.h"

#include "EEPARAMS.h"

C     Physical constants + functions of sigma and latitude
#include "com_physcon.h"

C     Vertical diffusion constants
#include "com_vdicon.h"

C-- Routine arguments:
c     _RL  UA(NGP,NLEV), VA(NGP,NLEV)
      _RL  dpFac(NGP,NLEV)
      _RL  SE(NGP,NLEV), RH(NGP,NLEV), QA(NGP,NLEV), QSAT(NGP,NLEV)
c     _RL  PHI(NGP,NLEV)

c     _RL  UTENVD(NGP,NLEV), VTENVD(NGP,NLEV)
      _RL  TTENVD(NGP,NLEV), QTENVD(NGP,NLEV)

      INTEGER  kGrd(NGP)
      INTEGER  bi,bj,myThid

#ifdef ALLOW_AIM

C-- Local variables:
      INTEGER J, K, Ktmp, NL1
      _RL  RSIG(NLEV)
      _RL  dSEdp(NGP,NLEV-1), DeltaPI(NLEV-1), factP

C- jmc: declare all local variables:
      _RL  CVDI(NGP), FSHCQ
      _RL  DRH0, DRH, DMSE, FLUXSE, FLUXQ
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C--   1. Initalization

C     N.B. In this routine, fluxes of dry static energy and humidity
C          are scaled in such a way that:
C          d_T/dt = d_F'(SE)/d_sigma,  d_Q/dt = d_F'(Q)/d_sigma

c_FM  NL1  = NLEV-1
c_FM  CSHC = DSIG(NLEV)/3600.
c_FM  CVDI = (SIGH(NL1)-SIGH(1))/((NL1-1)*3600.)

c_FM  FSHCQ  = CSHC/TRSHC
c_FM  FSHCSE = CSHC/(TRSHC*CP)

c_FM  FVDIQ  = CVDI/TRVDI
c_FM  FVDISE = CVDI/(TRVDS*CP)

      DO J=1,NGP
        NL1 = kGrd(J)-1
        CVDI(J) = 0.
        IF (NL1.GE.2) THEN
          CVDI(J) = (SIGH(NL1)-SIGH(1))/((NL1-1)*3600. _d 0)
        ENDIF
      ENDDO

      DO K=1,NLEV
        RSIG(K)=1./DSIG(K)
      ENDDO

      DO K=1,NLEV
        DO J=1,NGP
c         UTENVD(J,K) = 0.
c         VTENVD(J,K) = 0.
          TTENVD(J,K) = 0.
          QTENVD(J,K) = 0.
        ENDDO
      ENDDO

C ---------------------------------------------
C    Write Conditional stability based on Pot.Temp :
C    dSEdp(K) = Delta[Static-Energy] between 2 Plevels(k,k+1)
C    and corresponds to SE(K+1)-SE(K) in the original code
C -------
      DO K=1,NLEV-1
       factP = CP*SIGH(K)**(RD/CP)
       DO J=1,NGP
         dSEdp(J,K)=(SE(J,K+1)-SE(J,K))*factP
       ENDDO
       DeltaPI(K) = SIG(K+1)**(RD/CP) - SIG(K)**(RD/CP)
      ENDDO

C--   2. Shallow convection

      DO J=1,NGP
        Ktmp = kGrd(J)
        NL1 = Ktmp - 1
       IF (Ktmp.GE.2) THEN

        DRH0=RHGRAD*(SIG(Ktmp)-SIG(NL1))
        FSHCQ = DSIG(Ktmp)*dpFac(J,Ktmp)/(TRSHC*3600. _d 0)

c_FM    DMSE = (SE(J,NLEV)-SE(J,NL1))+ALHC*(QA(J,NLEV)-QSAT(J,NL1))
        DMSE = dSEdp(J,NL1)         + ALHC*(QA(J,Ktmp)-QSAT(J,NL1))
        DRH  = RH(J,Ktmp)-RH(J,NL1)

        IF (DMSE.GE.0.0) THEN

c_FM      FLUXSE         = FSHCSE*DMSE
          FLUXSE         = FSHCQ *DMSE/CP
          TTENVD(J,NL1)  = FLUXSE*RSIG(NL1)
          TTENVD(J,Ktmp) =-FLUXSE*RSIG(Ktmp)

          IF (DRH.GE.0.0) THEN
            FLUXQ          = FSHCQ*QSAT(J,Ktmp)*DRH
            QTENVD(J,NL1)  = FLUXQ*RSIG(NL1)
            QTENVD(J,Ktmp) =-FLUXQ*RSIG(Ktmp)
          ENDIF

        ELSE IF (DRH.GE.DRH0) THEN

c_FM      FLUXQ          = FVDIQ*QSAT(J,NL1)*DRH
          FLUXQ          = QSAT(J,NL1)*DRH*CVDI(J)/TRVDI
          QTENVD(J,NL1)  = FLUXQ*RSIG(NL1)
          QTENVD(J,Ktmp) =-FLUXQ*RSIG(Ktmp)

        ENDIF

       ENDIF
      ENDDO

C--   3. Vertical diffusion of moisture above the PBL

      DO J=1,NGP

        DO K=3,kGrd(J)-2

          DRH0=RHGRAD*(SIG(K+1)-SIG(K))

          DRH=RH(J,K+1)-RH(J,K)

          IF (DRH.GE.DRH0) THEN
c_FM        FLUXQ        = FVDIQ*QSAT(J,K)*DRH
            FLUXQ        = QSAT(J,K)*DRH*CVDI(J)/TRVDI
            QTENVD(J,K)  = QTENVD(J,K)  +FLUXQ*RSIG(K)
            QTENVD(J,K+1)= QTENVD(J,K+1)-FLUXQ*RSIG(K+1)
          ENDIF

        ENDDO

      ENDDO

C--   4. Damping of super-adiabatic lapse rate

      DO J=1,NGP
       DO K=1,kGrd(J)-1

c_FM     SE0 = SE(J,K+1)+SEGRAD*(PHI(J,K)-PHI(J,K+1))
         DMSE = dSEdp(J,K)
     &        +SEGRAD*CP*DeltaPI(K)*(SE(J,K+1)+SE(J,K))*0.5 _d 0

c_FM     IF (SE(J,K).LT.SE0) THEN
c_FM       FLUXSE        = FVDISE*(SE0-SE(J,K))
         IF (DMSE.GT.0.) THEN
           FLUXSE        = DMSE*CVDI(J)/(TRVDS*CP)
           TTENVD(J,K  ) = TTENVD(J,K  )+FLUXSE*RSIG(K)
           TTENVD(J,K+1) = TTENVD(J,K+1)-FLUXSE*RSIG(K+1)
         ENDIF

       ENDDO
      ENDDO
C--
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#endif /* ALLOW_AIM */

      RETURN
      END