C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/aim/Attic/phy_suflux.F,v 1.6 2002/09/27 20:05:11 jmc Exp $
C $Name:  $

#include "AIM_OPTIONS.h"

      SUBROUTINE SUFLUX (PSA,UA,VA,TA,QA,RH,QSAT,Vsurfsq,PHI,
     *                   PHI0,FMASK,TLAND,TSEA,SWAV,SSR,SLR,
     *                   DRAG,USTR,VSTR,SHF,EVAP,T0,Q0,QSAT0,SPEED0,
     &                   myThid)
C--
C--   SUBROUTINE SUFLUX (PSA,UA,VA,TA,QA,RH,PHI,
C--  *                   PHI0,FMASK,TLAND,TSEA,SWAV,
C--  *                   USTR,VSTR,SHF,EVAP)
C--
C--   Purpose: Compute surface fluxes of momentum, energy and moisture
C--   Input:   PSA    = norm. surface pressure [p/p0] (2-dim)
C--            UA     = u-wind                        (3-dim)
C--            VA     = v-wind                        (3-dim)
C--            TA     = temperature                   (3-dim)
C--            QA     = specific humidity [g/kg]      (3-dim)
C--            RH     = relative humidity             (3-dim)
C--            PHI    = geopotential                  (3-dim)
C--            PHI0   = surface geopotential          (2-dim)
C--            FMASK  = fractional land-sea mask         (2-dim)
C--            TLAND  = land-surface temperature         (2-dim)
C--            TSEA   =  sea-surface temperature         (2-dim)
C--            SWET   = soil wetness availability [0-1]  (2-dim)
C--   Output:  USTR   = u stress                         (2-dim)
C--            VSTR   = v stress                         (2-dim)
C--            SHF    = sensible heat flux               (2-dim)
C--            EVAP   = evaporation [g/(m^2 s)]          (2-dim)
C - added (jmc) :
C--           Vsurfsq = square of surface wind speed (2-dim,input)
C--                       ==> UA,VA are no longer used
C--            DRAG   = surface Drag term (= Cd*Rho*|V|) (2-dim,output)
C--                       ==> USTR,VSTR are no longer used
C--         myThid    = Instance number of this instance of the
C--                     the routine.
C--

      IMPLICIT NONE

C     Resolution parameters

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

#include "EEPARAMS.h"

#include "AIM_GRID.h"

C     Physical constants + functions of sigma and latitude

#include "com_physcon.h"

C     Surface flux constants

#include "com_sflcon.h"

C-- Routine arguments:
      INTEGER myThid
      _RL PSA(NGP), UA(NGP,NLEV), VA(NGP,NLEV), TA(NGP,NLEV),
     *     QA(NGP,NLEV), RH(NGP,NLEV), QSAT(NGP,NLEV), PHI(NGP,NLEV),
     *     PHI0(NGP), FMASK(NGP), TLAND(NGP), TSEA(NGP), SWAV(NGP),
     *     SSR(NGP), SLR(NGP)
      _RL Vsurfsq(NGP), DRAG(NGP)
      _RL USTR(NGP,3), VSTR(NGP,3), SHF(NGP,3), EVAP(NGP,3)
      _RL T0(NGP,2),Q0(NGP),QSAT0(NGP,2), SPEED0(NGP)

#ifdef ALLOW_AIM

C-- Local variables: 
      _RL U0(NGP),V0(NGP),DENVV(NGP,3)
      _RL WORK(NGP)
      INTEGER NL1(NGP)
      _RL AUX(NGP)
C
      _RL    pSurfs(NLEV)
      DATA    pSurfs /75 _d 2,250 _d 2,500 _d 2,775 _d 2,950 _d 2 /
C
      _RL    pSurfw(NLEV)
      DATA    pSurfw /150 _d 2,350 _d 2,650 _d 2,900 _d 2,1000 _d 2/
Cchdbg
c     INTEGER NPAS
c     SAVE NPAS
c     LOGICAL Ifirst
c     SAVE Ifirst
c     DATA Ifirst /.TRUE./
c     _RL T0moy1(NGP)
c     _RL T0moy2(NGP)
c     SAVE T0moy2, T0moy1
c     _RL denmoy1(NGP)
c     _RL denmoy2(NGP)
c     SAVE denmoy1, denmoy2
c     _RL Q0moy(NGP)
c     SAVE Q0moy

      INTEGER J
      _RL factwind2

C- jmc: declare all local variables: 
      _RL GTEMP0, GHUM0, RCP, PRD, VG2, DTHETAF, FSTAB, RDTH
      _RL FSLAND, FSSEA, CHLCP, CHSCP, QDUMMY, RDUMMY
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C--   1. Extrapolation of wind, temp, hum. and density to the surface
C
C     1.1 Wind components
C   
      DO J=1,NGP
         U0(J) = 0.0
         V0(J) = 0.0
        IF ( NLEVxyU(J,myThid) .GT. 0 ) THEN
         U0(J) = FWIND0*UA(J,NLEVxyU(J,myThid))
        ENDIF
        IF ( NLEVxyV(J,myThid) .GT. 0 ) THEN
         V0(J) = FWIND0*VA(J,NLEVxyV(J,myThid))
        ENDIF
      ENDDO
C
C     1.2 Temperature
C
      GTEMP0 = 1.D0-FTEMP0
      RCP = 1.D0/CP
      DO J=1,NGP
        NL1(J)=NLEVxy(J,myThid)-1
      ENDDO
C
      DO J=1,NGP
       IF ( NLEVxy(J,myThid) .GT. 0 ) THEN
        T0(J,1) = TA(J,NLEVxy(J,myThid))+WVI(NLEVxy(J,myThid),2)*
     &                          (TA(J,NLEVxy(J,myThid))-TA(J,NL1(J)))
ccccc        T0(J,2) = TA(J,NLEVxy(J))+RCP*(PHI(J,NLEVxy(J))-PHI0(J))
        T0(J,2) = TA(J,NLEVxy(J,myThid))*
     &         ((Psurfw(NLEVxy(J,myThid))/
     &           Psurfs(Nlevxy(J,myThid)))**(RD/CP))
       ELSE
        T0(J,1) = 273.16 _d 0
        T0(J,2) = 273.16 _d 0
       ENDIF
      ENDDO
C
      DO J=1,NGP
        T0(J,1) = FTEMP0*T0(J,1)+GTEMP0*T0(J,2)
      ENDDO
C
C     1.3 Spec. humidity
C
      GHUM0 = 1. _d 0-FHUM0
C
      DO J=1,NGP
       IF ( NLEVxy(J,myThid) .GT. 0 ) THEN
        WORK(J)=RH(J,Nlevxy(J,myThid))
cchdbg
c        WORK(J) = RH(J,NLEVxy(J))+WVI(NLEVxy(J),2)*
c     &                          (RH(J,NLEVxy(J))-RH(J,NL1(J)))
cchdbg
       ENDIF
      ENDDO
C
      CALL SHTORH (-1,NGP,T0,PSA,1. _d 0,Q0,WORK,QSAT0,myThid)
C
      DO J=1,NGP
       IF ( NLEVxy(J,myThid) .GT. 0 ) THEN
        Q0(J)=FHUM0*Q0(J)+GHUM0*QA(J,NLEVxy(J,myThid))
       ENDIF
      ENDDO

C     1.4 Density * wind speed (including gustiness factor)

      PRD = P0/RD
      VG2 = VGUST*VGUST
      factwind2 = FWIND0*FWIND0
C
      DO J=1,NGP
c_jmc   SPEED0(J)=SQRT(U0(J)*U0(J)+V0(J)*V0(J)+VG2)
        SPEED0(J)=SQRT(factwind2*Vsurfsq(J)+VG2)
        DENVV(J,3)=(PRD*PSA(J)/T0(J,1))*SPEED0(J)
c_jmc   DENVV(J,3)=(PRD*PSA(J)/T0(J,1))*
c_jmc*           SQRT(U0(J)*U0(J)+V0(J)*V0(J)+VG2)
cchdbg        DENVV(J,3)=0.7*(PRD*PSA(J)/T0(J,1))*
cchdbg     *           SQRT(U0(J)*U0(J)+V0(J)*V0(J)+VG2)
      ENDDO
C
C     1.5 Stability correction for heat/moisture fluxes
C
      DTHETAF = 3. _d 0
      FSTAB = 0.67 _d 0
      RDTH = FSTAB/DTHETAF
C
      DO J=1,NGP
        FSLAND=1.+MIN(DTHETAF,MAX(-DTHETAF,TLAND(J)-T0(J,2)))*RDTH
        FSSEA =1.+MIN(DTHETAF,MAX(-DTHETAF, TSEA(J)-T0(J,2)))*RDTH
        aux(J)=FSLAND
        DENVV(J,1)=DENVV(J,3)*FSLAND
        DENVV(J,2)=DENVV(J,3)*FSSEA
cchdbg        DENVV(J,1)=DENVV(J,3)
cchdbg        DENVV(J,2)=DENVV(J,3)
      ENDDO
C
C--   2. Computation of fluxes over land and sea
C
C     2.1 Wind stress
C
c_jmc DO J=1,NGP
c_jmc  IF ( NLEVxyu(J) .GT. 0  ) THEN
c_jmc   USTR(J,1) = -CDL*DENVV(J,3)*UA(J,NLEVxyu(J))
c_jmc   USTR(J,2) = -CDS*DENVV(J,3)*UA(J,NLEVxyu(J))
c_jmc  ENDIF
c_jmc  IF ( NLEVxyv(J) .GT. 0  ) THEN
c_jmc   VSTR(J,1) = -CDL*DENVV(J,3)*VA(J,NLEVxyv(J))
c_jmc   VSTR(J,2) = -CDS*DENVV(J,3)*VA(J,NLEVxyv(j))
c_jmc  ENDIF
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C   Compute surface drag term (= C_drag*|V| ) to allow direct computation 
C     of surface stress on C-grid.
C     add Land + Sea contributions ; Convert to surface pressure level 
      DO J=1,NGP
       IF ( NLEVxy(J,myThid) .GT. 0  ) THEN
        DRAG(J) = ( CDS+FMASK(J)*(CDL-CDS) ) * DENVV(J,3)
       ELSE
        DRAG(J) = 0.
       ENDIF
      ENDDO
C - Notes: 
C   because of a different mapping between the Drag and the Wind (A/C-grid)
C   the surface stress is computed later, in "External Forcing", 
C  => USTR,VSTR is no longer used. only here for diagnostic of old version.
      DO J=1,NGP
        USTR(J,3) = 0.
        VSTR(J,3) = 0.
       IF ( NLEVxyU(J,myThid) .GT. 0  ) 
     &  USTR(J,3) = -DRAG(J)*UA(J,NLEVxyU(J,myThid))
       IF ( NLEVxyV(J,myThid) .GT. 0  ) 
     &  VSTR(J,3) = -DRAG(J)*VA(J,NLEVxyV(J,myThid))
      ENDDO
c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C
C     2.2 Sensible heat flux (from clim. TS over land)
C
      CHLCP = CHL*CP
      CHSCP = CHS*CP
C
      DO J=1,NGP
        SHF(J,1) = CHLCP*DENVV(J,1)*(TLAND(J)-T0(J,1))
        SHF(J,2) = CHSCP*DENVV(J,2)*(TSEA(J) -T0(J,1))
      ENDDO
C
C ****************************************************
cchdbg
c     IF (Ifirst) then 
c       npas=0.
c       do J=1,NGP
c         T0moy1(J)=0. 
c         T0moy2(J)=0. 
c         denmoy1(J)=0. 
c         denmoy2(J)=0. 
c         Q0moy(J)=0.
c       enddo
c       ifirst=.false.
c     ENDIF

c     npas=npas+1
c     DO J=1,NGP
c       T0moy1(J)=T0moy1(J) + T0(J,1)
c       T0moy2(J)=T0moy2(J) + T0(J,2)
c       denmoy1(J)=denmoy1(J) + DENVV(J,1)
c       denmoy2(J)=denmoy2(J) + DENVV(J,2)
c       Q0moy(J)=Q0moy(J) + Q0(J)
c     ENDDO

c     if(npas.eq.5760) then
c       DO J=1,NGP
c         T0moy1(J)=T0moy1(J)/float(npas)
c         T0moy2(J)=T0moy2(J)/float(npas)
c         denmoy1(J)=denmoy1(J)/float(npas)
c         denmoy2(J)=denmoy2(J)/float(npas)
c         Q0moy(J)=Q0moy(J)/float(npas)
c       ENDDO

c       open(73,file='Tmoy1',form='unformatted')
c       write(73) T0moy1
c       close(73)
c       open(74,file='Tmoy2',form='unformatted')
c       write(74) T0moy2
c       close(74)
c       open(73,file='denmoy1',form='unformatted')
c       write(73) denmoy1
c       close(73)
c       open(74,file='denmoy2',form='unformatted')
c       write(74) denmoy2
c       close(74)
c       open(74,file='Q0moy',form='unformatted')
c       write(74) Q0moy
c       close(74)
c     ENDIF
cchdbg
C ****************************************************
C
c        CALL DUMP_WRITE2D ( NGP,1,'T0.', nIter,T0 , iErr)
c        CALL DUMP_WRITE2D ( NGP,3,'DENVV.', nIter,DENVV , iErr)
c        CALL DUMP_WRITE2D ( NGP,1,'TSEA.', nIter,TSEA , iErr)
c        CALL DUMP_WRITE2D ( NGP,1,'TLAND.', nIter,TLAND , iErr)
c        CALL DUMP_WRITE2D ( NGP,1,'AUX.', nIter,aux , iErr)
C
C     2.3 Evaporation
C
      CALL SHTORH (0,NGP,TLAND,PSA,1. _d 0,QDUMMY,RDUMMY,QSAT0(1,1),myThid)
      CALL SHTORH (0,NGP,TSEA ,PSA,1. _d 0,QDUMMY,RDUMMY,QSAT0(1,2),myThid)

      DO J=1,NGP
Cdj    EVAP(J,1) = CHL*DENVV(J,1)*MAX(0. _d 0,SWAV(J)*QSAT0(J,1)-Q0(J))
       EVAP(J,1)=CHL*DENVV(J,1)*MIN(1. _d 0,SWAV(J))*(QSAT0(J,1)-Q0(J))
cdj try new scheme : assume that the net heat flux on land is zero
cdj                  at all time and at all points (this is equivalent
cdj                  to say that the land has a zero heat capacity)
cdj        EVAP(J,1) = SSR(J)-SLR(J)-SHF(J,1)
       EVAP(J,2) = CHS*DENVV(J,2)*(QSAT0(J,2)-Q0(J))
C - test(jmc): only positive evaporation :
c      EVAP(J,1)=CHL*DENVV(J,1)*MAX(0. _d 0,SWAV(J)*(QSAT0(J,1)-Q0(J)))
c      EVAP(J,2)=CHS*DENVV(J,2)*MAX(0. _d 0,QSAT0(J,2)-Q0(J))
      ENDDO


C--   3. Weighted average of fluxes according to land-sea mask

      DO J=1,NGP
c_jmc   USTR(J,3) = USTR(J,2)+FMASK(J)*(USTR(J,1)-USTR(J,2))
c_jmc   VSTR(J,3) = VSTR(J,2)+FMASK(J)*(VSTR(J,1)-VSTR(J,2))
         SHF(J,3) =  SHF(J,2)+FMASK(J)*( SHF(J,1)- SHF(J,2))
        EVAP(J,3) = EVAP(J,2)+FMASK(J)*(EVAP(J,1)-EVAP(J,2))
cdj
        QSAT0(J,1) = QSAT0(J,2)+FMASK(J)*(QSAT0(J,1)-QSAT0(J,2))
cdj
      ENDDO

C--
#endif /* ALLOW_AIM */ 

      RETURN
      END