C $Header: /u/gcmpack/MITgcm/pkg/aim_v23/phy_suflux_prep.F,v 1.3 2004/06/25 18:24:17 jmc Exp $
C $Name: $
#include "AIM_OPTIONS.h"
CBOP
C !ROUTINE: SUFLUX_PREP
C !INTERFACE:
SUBROUTINE SUFLUX_PREP(
I PSA,TA,QA,RH,ThA,Vsurf2,WVS,CLAT,FOROG,
I FMASK,TLAND,TSEA,TSICE,SSR,
O SPEED0,DRAG,DENVV,dTskin,T1,T0,Q0,
I kGrd,bi,bj,myThid)
C !DESCRIPTION: \bv
C *==========================================================*
C | S/R SUFLUX_PREP
C | o prepare surface flux calculation
C *==========================================================*
C | o contain 1rst part of original S/R SUFLUX (Speedy code)
C *==========================================================*
C--
C-- SUBROUTINE SUFLUX (PSA,UA,VA,TA,QA,RH,PHI,
C-- & PHI0,FMASK,TLAND,TSEA,SWAV,SSR,SLRD,
C-- & USTR,VSTR,SHF,EVAP,SLRU,
C-- & TSFC,TSKIN,U0,V0,T0,Q0)
C--
C-- Purpose: Compute surface fluxes of momentum, energy and moisture,
C-- and define surface skin temperature from energy balance
C *==========================================================*
C \ev
C !USES:
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 Surface flux constants
#include "com_sflcon.h"
C !INPUT/OUTPUT PARAMETERS:
C == Routine Arguments ==
C-- Input:
C PSA :: norm. surface pressure [p/p0] (2-dim)
C TA :: temperature (3-dim)
C QA :: specific humidity [g/kg] (3-dim)
C RH :: relative humidity [0-1] (3-dim)
C ThA :: Pot.temperature [K] (3-dim)
C Vsurf2 :: square of surface wind speed (2-dim,input)
C ==> UA,VA are no longer used
C WVS :: weights for near surf interp (2-dim)
C CLAT :: cos(lat) (2-dim)
C FOROG :: orographic factor (surf. drag) (2-dim)
C FMASK :: fraction land - sea - sea-ice (2.5-dim)
C TLAND :: land-surface temperature (2-dim)
C TSEA :: sea-surface temperature (2-dim)
C TSICE :: sea-ice surface temperature (2-dim)
C SSR :: sfc sw radiation (net flux) (2-dim)
C-- Output:
C SPEED0 :: effective surface wind speed (2-dim)
C DRAG :: surface Drag term (= Cd*Rho*|V|)(2-dim)
C ==> USTR,VSTR are no longer used
C DENVV :: surface flux (sens,lat.) coeff. (=Rho*|V|) [kg/m2/s]
C dTskin :: temp. correction for daily-cycle heating [K]
C T1 :: near-surface air temperature (from Pot.Temp)
C T0 :: near-surface air temperature (2-dim)
C Q0 :: near-surface sp. humidity [g/kg](2-dim)
C-- Input:
C kGrd :: Ground level index (2-dim)
C bi,bj :: tile index
C myThid :: Thread number for this instance of the routine
C--
_RL PSA(NGP), TA(NGP,NLEV), QA(NGP,NLEV), RH(NGP,NLEV)
_RL ThA(NGP,NLEV)
_RL Vsurf2(NGP), WVS(NGP), CLAT(NGP), FOROG(NGP)
_RL FMASK(NGP,3), TLAND(NGP), TSEA(NGP), TSICE(NGP)
_RL SSR(NGP)
_RL SPEED0(NGP), DRAG(NGP,0:3), T1(NGP), DENVV(NGP)
_RL dTskin(NGP), T0(NGP), Q0(NGP)
INTEGER kGrd(NGP)
INTEGER bi,bj,myThid
CEOP
#ifdef ALLOW_AIM
C-- Local variables:
_RL QSAT0(NGP,2)
INTEGER J, Ktmp, NL1
_RL tmpRH(NGP)
_RL factWind2, kappa
C- jmc: declare all local variables:
_RL GTEMP0, GHUM0, RCP, PRD, VG2
c _RL RDTH, FSLAND, FSSEA, FSSICE
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C-- 1. Extrapolation of wind, temp, hum. and density to the surface
C 1.1 Wind components
c DO J=1,NGP
c U0(J) = 0.0
c V0(J) = 0.0
c Ktmp = kGrd(J)
c IF ( Ktmp.GT.0 ) THEN
c U0(J) = FWIND0*UA(J,Ktmp)
c V0(J) = FWIND0*VA(J,Ktmp)
c ENDIF
c ENDDO
C 1.2 Temperature
GTEMP0 = 1.-FTEMP0
RCP = 1. _d 0 /CP
kappa = RD/CP
C
DO J=1,NGP
Ktmp = kGrd(J)
NL1 = Ktmp-1
IF ( Ktmp.GT.1 ) THEN
c_FM T0(J) = TA(J,NLEV)+WVI(NLEV,2)*(TA(J,NLEV)-TA(J,NL1))
c_FM T1(J) = TA(J,NLEV)+RCP*(PHI(J,NLEV)-PHI0(J))
T0(J) = TA(J,Ktmp) + WVS(J)*(TA(J,Ktmp)-TA(J,NL1))
Cjmc: used previously but not valid with partial cell !
c T1(J) = TA(J,Ktmp)*(SIGH(Ktmp)/SIG(Ktmp))**kappa
T1(J) = ThA(J,Ktmp)*(PSA(J)**kappa)
tmpRH(J)=RH(J,Ktmp)
ELSE
T0(J) = 273.16 _d 0
T1(J) = 273.16 _d 0
tmpRH(J)= 0.
ENDIF
ENDDO
DO J=1,NGP
c T0(J) = FTEMP0*T0(J)+GTEMP0*T1(J)
T0(J) = FTEMP0*MIN(T0(J),T1(J))+GTEMP0*T1(J)
ENDDO
C 1.3 Spec. humidity
GHUM0 = 1.-FHUM0
CALL SHTORH (-1,NGP,T0, PSA, 1. _d 0, Q0, tmpRH, QSAT0, myThid)
DO J=1,NGP
IF ( kGrd(J) .GT. 0 ) THEN
Q0(J)=FHUM0*Q0(J)+GHUM0*QA(J,kGrd(J))
ENDIF
ENDDO
C 1.4 Density * wind speed (including gustiness factor)
PRD = P0/RD
VG2 = VGUST*VGUST
factWind2 = FWIND0*FWIND0
DO J=1,NGP
c_FM DENVV(J)=(PRD*PSA(J)/T0(J))*
c_FM & SQRT(U0(J)*U0(J)+V0(J)*V0(J)+VG2)
SPEED0(J)=SQRT(factWind2*Vsurf2(J)+VG2)
DENVV(J)=(PRD*PSA(J)/T0(J))*SPEED0(J)
ENDDO
C 1.5 Define effective skin temperature to compensate for
C non-linearity of heat/moisture fluxes during the daily cycle
C Tskin = Tland + dTskin
DO J=1,NGP
dTskin(J)=CTDAY*CLAT(J)*SSR(J)*PSA(J)
ENDDO
C-- 2. Computation of fluxes over land and sea
C 2.1 Wind stress
C Orographic correction
DO J=1,NGP
c CDENVV(J,1)=CDL*DENVV(J)*FOROG(J)
c CDENVV(J,2)=CDS*DENVV(J)
DRAG(J,1) = CDL*DENVV(J)*FOROG(J)
DRAG(J,2) = CDS*DENVV(J)
DRAG(J,3) = CDS*DENVV(J)
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 Here compute only surface drag term (= C_drag*Rho*|V| )
c DO J=1,NGP
c USTR(J,1) = -CDENVV(J,1)*UA(J,NLEV)
c VSTR(J,1) = -CDENVV(J,1)*VA(J,NLEV)
c USTR(J,2) = -CDENVV(J,2)*UA(J,NLEV)
c VSTR(J,2) = -CDENVV(J,2)*VA(J,NLEV)
c ENDDO
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#endif /* ALLOW_AIM */
RETURN
END
SUBROUTINE SFLSET (PHI0, FOROG, bi,bj,myThid)
C--
C-- SUBROUTINE SFLSET (PHI0)
C--
C-- Purpose: compute orographic factor for land surface drag
C-- Input: PHI0 = surface geopotential (2-dim)
C Output: FOROG = orographic factor (surf. drag) (2-dim)
C-- (originally in common blocks: SFLFIX)
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 Surface flux constants
#include "com_sflcon.h"
C-- Routine arguments:
INTEGER bi,bj,myThid
_RL PHI0(NGP)
_RL FOROG(NGP)
#ifdef ALLOW_AIM
C-- Local variables:
INTEGER J
_RL RHDRAG
RHDRAG = 1./(GG*HDRAG)
DO J=1,NGP
FOROG(J) = 1. _d 0
& + FHDRAG*(1. _d 0 - EXP(-MAX(PHI0(J),0. _d 0)*RHDRAG) )
ENDDO
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#endif /* ALLOW_AIM */
RETURN
END