!copyright (C) 2001 MSC-RPN COMM %%%RPNPHY%%%
*** S/P DIASURF
*
#include "phy_macros_f.h"
SUBROUTINE DIASURF1(UZ,VZ,TZ,QZ,NI,U,V,TG,QG,Z0,Z0T,ILMO,ZA, 5,4
1 H,UE2,FTEMP,FVAP,ZU,ZT,LAT)
#include "impnone.cdk"
INTEGER NI
REAL ZT(NI),ZU(NI)
REAL UZ(NI),VZ(NI),TZ(NI),QZ(NI),ZA(NI),U(NI),V(NI)
REAL TG(NI),QG(NI),UE2(NI),FTEMP(NI),FVAP(NI)
REAL LAT(NI),ILMO(NI),Z0T(NI),Z0(NI),H(NI)
*Author
* Yves Delage (Aug1990)
*
*Revision
* 001 G. Pellerin(JUN94)
* Adaptation to new surface formulation
* 002 B. Bilodeau (Nov 95) - Replace VK by KARMAN
* 003 R. Sarrazin (Jan 96) - Prevent problems if zu < za
* 004 G. Pellerin (Feb 96) - Rewrite stable formulation
* 005 Y. Delage and B. Bilodeau (Jul 97) - Cleanup
* 006 Y. Delage (Feb 98) - Addition of HMIN
* 007 G. Pellerin (Mai 03) - Conversion IBM
* - calls to vslog routine (from massvp4 library)
*
*Object
* to calculate the diagnostic values of U, V, T, Q
* near the surface (ZU and ZT)
*
*Arguments
*
* - Output -
* UZ U component of the wind at Z=ZU
* VZ V component of the wind at Z=ZU
* TZ temperature in kelvins at Z=ZT
* QZ specific humidity at Z=ZT
*
* - Input -
* NI number of points to process
* U U component of wind at Z=ZA
* V V component of wind at Z=ZA
* TG temperature at the surface (Z=0) in Kelvins
* QG specific humidity
* PS surface pressure at the surface
* ILMO inverse of MONIN-OBUKHOV lenth
* H height of boundary layer
* Z0 roughness lenth for winds
* Z0T roughness lenth for temperature and moisture
* FTEMP temperature flux at surface
* FVAP vapor flux at surface
* ZA top of surface layer
* ZU heights for computation of wind components
* ZT heights for computation of temperature and moisture
* LAT LATITUDE
*
REAL ANG,ANGI,VITS,LZZ0
REAL CT,DANG,CM,ANGMAX
REAL X,X0,Y,Y0,FH,FM,UE
REAL RAC3
SAVE ANGMAX
INTEGER J
*
*
*******************************************************
* AUTOMATIC ARRAYS
*******************************************************
*
AUTOMATIC ( LZZ0T , REAL , (NI) )
AUTOMATIC ( LZZ0U , REAL , (NI) )
*
*******************************************************
*Implicites
#include "surfcon.cdk"
#include "consphy.cdk"
*
#include "dfsn.cdk"
#include "fsn.cdk"
*
*MODULES
DATA ANGMAX /0.85/
RAC3=SQRT(3.)
DO J=1,NI
LZZ0T(J)=ZT(J)/Z0T(J)+1
LZZ0U(J)=ZU(J)/Z0(J)+1
ENDDO
*
call vslog(LZZ0T,LZZ0T,NI)
call vslog(LZZ0U,LZZ0U,NI)
*
DO 10 J=1,NI
IF(ILMO(J).LE.0.) THEN
***********************************************************************
* BRANCHE INSTABLE
***********************************************************************
LZZ0=LZZ0T(J)
Y=(1-BETA*CI*(ZT(J)+Z0T(J))*ILMO(J))**(1./3)
Y0=(1-BETA*CI*Z0T(J)*ILMO(J))**(1./3)
FH=BETA*(LZZ0+1.5*ALOG((Y0**2+Y0+1)/(Y**2+Y+1))+RAC3*
1 ATAN(RAC3*2*(Y-Y0)/((2*Y0+1)*(2*Y+1)+3)))
CT=KARMAN/FH
LZZ0=LZZ0U(J)
X=(1-BETA*CI*(ZU(J)+Z0(J))*ILMO(J))**(1./6)
X0=(1-BETA*CI*Z0(J)*ILMO(J))**(1./6)
FM=LZZ0+ALOG((X0+1)**2*SQRT(X0**2-X0+1)*(X0**2+X0+1)**1.5
1 /((X+1)**2*SQRT(X**2-X+1)*(X**2+X+1)**1.5))
2 +RAC3*ATAN(RAC3*((X**2-1)*X0-(X0**2-1)*X)/
3 ((X0**2-1)*(X**2-1)+3*X*X0))
CM=KARMAN/FM
ELSE
***********************************************************************
* BRANCHE STABLE
***********************************************************************
unsl=ilmo(j)
hi=1/MAX(HMIN,H(J),(ZA(J)+2*Z0(J))*factn,factn/d(ILMO(J)))
LZZ0=LZZ0T(J)
fh=BETA*(LZZ0+psi(ZT(J)+Z0T(J),hi)-psi(Z0T(J),hi))
CT=KARMAN/fh
LZZ0=LZZ0U(J)
fm=LZZ0+psi(zu(J)+Z0(J),hi)-psi(Z0(J),hi)
CM=KARMAN/fm
ENDIF
***********************************************************************
* BRANCHE COMMUNE
***********************************************************************
UE=SQRT(UE2(J))
TZ(J)=TG(J)-FTEMP(J)/(CT*UE)-GRAV/CPD*ZT(J)
QZ(J)=QG(J)-FVAP(J)/(CT*UE)
VITS=UE/CM
* CALCUL DU CHANGEMENT D'ANGLE
DANG=SIGN(MIN(1.,max(0.,(ZA(J)-ZU(J)))/H(J)*.20*PI),LAT(J))
* DANG= (ZA(J)-ZU(J))/H(J)*ANGMAX*SIN(LAT(J))
ANGI=ATAN2(V(J),SIGN(ABS(U(J))+1.e-05,U(J)))
*
IF(ILMO(J).GE.0.) THEN
ANG=ANGI+DANG
ELSE
ANG=ANGI
ENDIF
UZ(J)=VITS*COS(ANG)
VZ(J)=VITS*SIN(ANG)
10 CONTINUE
RETURN
END