!copyright (C) 2001 MSC-RPN COMM %%%RPNPHY%%%
***S/P FLXSURF2
*
SUBROUTINE FLXSURF2(CMU, CTU, RIB, FTEMP, FVAP, ILMO, 5,12
X UE, FCOR, TA , QA , ZA, VA,
Y TG , QG , H , Z0 , Z0T,
% LZZ0, LZZ0T, FM, FH, DFM, DFH,N )
#include "impnone.cdk"
INTEGER N
REAL CMU(N),CTU(N),RIB(N),FCOR(N),ILMO(N)
REAL FTEMP(N),FVAP(N),TA(N),QA(N),ZA(N),VA(N)
REAL TG(N),QG(N),H(N),Z0(N),UE(N)
REAL Z0T(N),LZZ0(N),LZZ0T(N)
REAL fm(N),fh(N),dfm(N),dfh(N)
*
*Author
* Y.Delage (Jul 1990)
*Revision
* 001 G. Pellerin (Jun 94) New function for unstable case
* 002 G. Pellerin (Jui 94) New formulation for stable case
* 003 B. Bilodeau (Nov 95) Replace VK by KARMAN
* 004 M. Desgagne (Dec 95) Add safety code in function ff
* and ensures that RIB is non zero
* 005 R. Sarrazin (Jan 96) Correction for H
* 006 C. Girard (Nov 95) - Diffuse T instead of Tv
* 007 G. Pellerin (Feb 96) Revised calculation for H (stable)
* 008 G. Pellerin (Feb 96) Remove corrective terms to CTU
* 009 Y. Delage and B. Bilodeau (Jul 97) - Cleanup
* 010 Y. Delage (Feb 98) - Addition of HMIN
* 011 D. Talbot and Y. Delage (Jan 02) -
* Correct bug of zero divide by dg in loop 35
*
*
*Object
* to calculate surface layer transfer coefficients and fluxes
*
*Arguments
*
* - Output -
* CMU transfer coefficient of momentum times UE
* CTU transfer coefficient of temperature times UE
* RIB bulk Richardson number
* FTEMP temperature flux
* FVAP vapor flux
* ILMO (1/length of Monin-Obukov)
* UE friction velocity (squared for output)
*
* - Input -
* FCOR Coriolis factor
* ZA Height of the surface layer
* TA potential temperature at last predictive level
* QA specific humidity " " " "
* VA wind speed " " " "
* TG surface temperature
* QG specific humidity at the surface
*
* - Input/Output -
* H height of the boundary layer
*
* - Input -
* Z0 roughness length for momentum flux calculations
* Z0T roughness length for heat/moisture flux calculations
* LZZ0 work field (N)
* LZZ0T work field (N)
*
* - Output -
* FM momentum function
* FH heat function
* DFM derivative of FM
* DFH derivative of FH
*
* - Input -
* N horizontal dimension
*
*Notes
* SEE DELAGE AND GIRARD BLM 58 (19-31)
* " BLM 82 (23-48)
*
*IMPLICITES
*
#include "surfcon.cdk"
#include "consphy.cdk"
*MODULES
*
EXTERNAL BORNES
*
**
*
INTEGER J
INTEGER IT,ITMAX
REAL HMAX,CORMIN
REAL RAC3,CM,CT
REAL F,G,DG
real*8 X,X0,Y,Y0
REAL UNSHN,UNSH,HE,HS
REAL UNSLN,DFUSL,FMN,FHN
REAL DTHV,TVA,TVS
REAL HL,AA,HP
SAVE HMAX,CORMIN,EPSLN,AA
SAVE ITMAX
SAVE N0RIB
*
#include "ribcom.cdk"
DATA CORMIN, HMAX /0.7E-4 , 1000.0/
DATA ITMAX / 4 /
DATA AA / 0.6 /
*
#include "dfsn.cdk"
#include "fsn.cdk"
*
RAC3=sqrt(3.0)
*
DO 1 J=1,N
*
* CALCUL DU NOMBRE DE RICHARDSON
*
tva=(1+DELTA*QA(J))*TA(J)
tvs=(1+DELTA*QG(J))*TG(J)
dthv=tva-tvs
RIB(J)=GRAV/(tvs+0.5*dthv)*ZA(J)*dthv/(VA(J)**2)
if (rib(j).ge.0.0) rib(j) = max(rib(j), n0rib)
if (rib(j).lt.0.0) rib(j) = min(rib(j),-n0rib)
*
* CALCULS AUXILIAIRES
*
LZZ0(J)=ALOG(1+ZA(J)/Z0(J))
LZZ0T(J)=ALOG(1+ZA(J)/Z0T(J))
1 CONTINUE
*
*********************************************************************
* BRANCHE STABLE
*********************************************************************
DO 5 J=1,N
IF(RIB(J).GT.0.) THEN
* RELATION ENTRE 1/L et RIB EN PREMIERE APPROXIMATION
ILMO(J)=RIB(J)*(LZZ0(J)+AS*RIB(J))**2/
% ((LZZ0T(J)+AS*RIB(J))*ZA(J))
H(J)=MAX(H(J),HMIN,(ZA(J)+2*Z0(J))*factn,factn/d(ILMO(J)))
UNSH=1/H(J)
unsl=ILMO(J)
fm(J)=LZZ0(J)+psi(ZA(J)+Z0(J),unsh)-psi(Z0(J),unsh)
fh(J)=BETA*(LZZ0T(J)+psi(ZA(J)+Z0T(J),unsh)-psi(Z0T(J),unsh))
dfm(J)=0.
dfh(J)=0.
ENDIF
5 CONTINUE
* DEBUT DES ITERATIONS POUR CALCULER LA VALEUR DE ILMO ET H
DO 35 IT=1,ITMAX
DO 35 J=1,N
IF(RIB(J).GT.0.) THEN
g=RIB(J)-ZA(J)*ILMO(J)*fh(J)/fm(J)**2
dg=-ZA(J)*fh(J)/fm(J)**2*(1+ILMO(J)*(dfh(J)/fh(J)-2*dfm(J)/fm(J)))
unsln=MAX(EPSLN,ILMO(J)- g/sign(max(abs(dg),1.e-9),dg))
H(J)=MAX(H(J),HMIN,(ZA(J)+2*Z0(J))*factn,factn/d(unsln))
unshn=1./H(J)
unsl=unsln
fmn=LZZ0(J)+psi(ZA(J)+Z0(J),unshn)-psi(Z0(J),unshn)
fhn=BETA*(LZZ0T(J)+psi(ZA(J)+Z0T(J),unshn)-psi(Z0T(J),unshn))
dfusl=unsln-ILMO(J)
dfm(J)=(fmn-fm(J))/sign(max(EPSLN,abs(dfusl)),dfusl)
dfh(J)=(fhn-fh(J))/sign(max(EPSLN,abs(dfusl)),dfusl)
* formulation pour h=BS*sqrt(L*UE/F)
ILMO(J)=unsln
F=MAX(ABS(FCOR(J)),CORMIN)
hs=BS*sqrt(KARMAN*VA(J)/(ILMO(J)*F*fmn))
hl=(ZA(J)+2*Z0(J))*factn
hp=aa*hs+(1-aa)*h(j)
H(J)=MAX(HMIN,hp,hl,factn/d(ILMO(J)))
unsh=1/H(J)
fm(J)=LZZ0(J)+psi(ZA(J)+Z0(J),unsh)-psi(Z0(J),unsh)
fh(J)=BETA*(LZZ0T(J)+psi(ZA(J)+Z0T(J),unsh)-psi(Z0T(J),unsh))
*
***********************************************************************
* BRANCHE INSTABLE
***********************************************************************
*
ELSE
ILMO(J)=MIN(0.,ILMO(J))
X=(1-CI*(ZA(J)+Z0(J))*BETA*ILMO(J))**(1./6)
X0=(1-CI*Z0(J)*BETA*ILMO(J))**(1./6.)
FM(J)=LZZ0(J)+DLOG((X0+1)**2*SQRT(X0**2-X0+1)*(X0**2+X0+1)**1.5
% /((X+1)**2*SQRT(X**2-X+1)*(X**2+X+1)**1.5))
% +RAC3*ATAN(RAC3*((X**2-1)*X0-(X0**2-1)*X)/
% ((X0**2-1)*(X**2-1)+3*X*X0))
Y=(1-CI*(ZA(J)+Z0T(J))*BETA*ILMO(J))**(1./3)
Y0=(1-CI*Z0T(J)*BETA*ILMO(J))**(1./3)
FH(J)=BETA*(LZZ0T(J)+1.5*DLOG((Y0**2+Y0+1)/(Y**2+Y+1))+RAC3*
% ATAN(RAC3*2*(Y-Y0)/((2*Y0+1)*(2*Y+1)+3)))
IF(IT.LT.ITMAX) THEN
G=RIB(J)-FH(J)/FM(J)**2*ZA(J)*ILMO(J)
DG=-ZA(J)*FH(J)/FM(J)**2*(1+1/FH(J)*(1/Y-1/Y0)-2/FM(J)*
% (1/X-1/X0))
ILMO(J)=ILMO(J)-G/DG
ENDIF
ENDIF
35 CONTINUE
*
***********************************************************************
* BRANCHE COMMUNE
***********************************************************************
DO 80 J=1,N
CM=KARMAN/FM(J)
CT=KARMAN/FH(J)
UE(J)=VA(J)*CM
CMU(J)=CM*UE(J)
CTU(J)=CT*UE(J)
* la couche limite neutre 0.3*ue/coriolis
F=MAX(ABS(FCOR(J)),CORMIN)
he=max(HMIN,0.3*UE(J)/F,ZA(J)+2*Z0(J))
if (rib(j).gt.0.0) then
* cas stable
H(J)=MIN(H(J),he,hmax)
else
* cas instable
H(J)=MIN(he,hmax)
endif
UE(J)=UE(J)*UE(J)
FTEMP(J)=-CTU(J)*(TA(J)-TG(J))
FVAP(J)=-CTU(J)*(QA(J)-QG(J))
80 CONTINUE
RETURN
END