C $Header: /u/gcmpack/MITgcm/pkg/aim/phy_suflux.F,v 1.7 2003/12/24 00:25:36 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(1), RDUMMY(1)
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