!copyright (C) 2001 MSC-RPN COMM %%%CCC%%%
***S/P SGOFLX3
*
subroutine sgoflx3(uu,vv,utend,vtend, 1,10
+ tth,ttf,ss,ssh,
+ ilev,ilg,il1,il2,
+ grav,rgas,rgocp,tau,taufac,
+ gc,height,slope,xcent,mtdir,
+ psurf,fcor,
+ gwdrag,blocking,orolift,leewave,
+ applytend)
*
#include "phy_macros_f.h"
#include "impnone.cdk"
*
*
logical gwdrag,blocking,orolift,leewave,applytend
*
integer ilev,ilg,il1,il2
real grav,rgas,rgocp,tau,taufac
real*8 uu(ilg,ilev), vv(ilg,ilev), utend(ilg,ilev),
+ vtend(ilg,ilev), tth(ilg,ilev), ttf(ilg,ilev),
+ ss(ilg,ilev), ssh(ilg,ilev), gc(ilg),
+ height(ilg), slope(ilg), xcent(ilg),
+ mtdir(ilg), psurf(ilg), fcor(ilg)
*
*Author
* A. Zadra - May 2002 - From lin_sgoflx1
*
*Revisions
* 001 B. Bilodeau and A. Zadra (Apr 2003) - Add smoothing of BV frequency.
* 002 J.-P. Toviessi (May 2003) - IBM conversion
* - calls to vsqrt routine (from massvp4 library)
* - calls to vdiv routine (from massvp4 library)
* - calls to vrec routine (from massvp4 library)
* - calls to vlog routine (from massvp4 library)
* 003 A. Zadra and B. Bilodeau (Oct 2003) - Modifications to the
* new code in order to validate with oper. code
*
*
*Object
* Simplified version of subgrid orographic drag (sgoflx2) scheme:
* - reduced, non-smoothed buoyancy frequency
* - shortened gravity-wave drag (McFarlane 87)
* - shortened low-level blocking (Lott & Miller 97)
* - orographic lift (not yet included)
* - lee-wave breaking (not yet included)
*
*
*Arguments
* - Input/Output -
* UU U component of wind as input
* U component of wind modified by the gravity wave
* drag as output (if applytend is .true.)
* VV V component of wind as input
* V component of wind modified by the gravity wave
* drag as output (if applytend is .true.)
*
* - Input -
* TTH virtual temperature level means
* TTF virtual temperature
* SS sigma levels
* SSH intermediate levels
*
* - Output -
* UTEND total tendency on U (gwd + blocking)
* wind
* VTEND total tendency on V (gwd + blocking)
*
* - Input -
* ILEV number of levels
* ILG total horizontal dimension
* IL1 1st dimension of U,V,T to start calculation
* IL2 index on horizontal dimension to stop calculation
* GRAV gravitational constant
* RGAS gas constant
* RGOCP CAPPA (see thermodynamic constants)
* TAU timestep times a factor: 1 for two time-level models
* 2 for three time-level models
* TAUFAC 1/(length scale).
* GC ground cover (land-sea mask) : between 0(sea) and 1(land)
* HEIGHT launching height (variance associated with unresolved orography)
* SLOPE mountain slope
* XCENT mountain eccentricity
* MTDIR mountain direction
* PSURF surface pressure
* FCOR Coriolis factor
* GWDRAG .true. for gravity wave drag
* BLOCKING .true. for blocking
* OROLIFT .true. for orographic lift
* LEEWAVE .true. for lee wave parameterization
* APPLYTEND.true. if tendencies are to be applied in this subroutine
*
*
**
*
************************************************************************
* AUTOMATIC ARRAYS
************************************************************************
*
AUTOMATIC ( VMOD , REAL*8 , (ILG ) )
AUTOMATIC ( UUB , REAL*8 , (ILG ) )
AUTOMATIC ( VVB , REAL*8 , (ILG ) )
AUTOMATIC ( DRAG , INTEGER , (ILG ) )
AUTOMATIC ( UB , REAL*8 , (ILG ) )
AUTOMATIC ( VB , REAL*8 , (ILG ) )
AUTOMATIC ( VMODB , REAL*8 , (ILG ) )
AUTOMATIC ( ENV , REAL*8 , (ILG ) )
AUTOMATIC ( SLP2 , REAL*8 , (ILG ) )
AUTOMATIC ( SLPF , REAL*8 , (ILG ) )
AUTOMATIC ( GAMMA , REAL*8 , (ILG ) )
AUTOMATIC ( THETA , REAL*8 , (ILG ) )
AUTOMATIC ( IZT1 , INTEGER , (ILG ) )
AUTOMATIC ( IZT2 , INTEGER , (ILG ) )
AUTOMATIC ( IZT3 , INTEGER , (ILG ) )
AUTOMATIC ( IZB , INTEGER , (ILG ) )
AUTOMATIC ( HBLK , REAL*8 , (ILG ) )
AUTOMATIC ( UAV , REAL*8 , (ILG ) )
AUTOMATIC ( VAV , REAL*8 , (ILG ) )
AUTOMATIC ( VELAV , REAL*8 , (ILG ) )
AUTOMATIC ( DELZ , REAL*8 , (ILG ) )
AUTOMATIC ( FDIR , REAL*8 , (ILG ) )
AUTOMATIC ( BLOFF , REAL*8 , (ILG ) )
c
AUTOMATIC ( U , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( V , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( TF , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( TH , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( S , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( SH , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( SEXPK , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( SHEXPK , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( BVFREQ , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( UTENDGWD , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( VTENDGWD , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( VELN , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( ASQ , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( ASQI , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( ASQS , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( DFAC , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( DEPFAC , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( GRAD , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( DENFAC , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( ZB , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( UTENDLLB , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( VTENDLLB , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( UTENDLFT , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( VTENDLFT , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( UTENDLWB , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( VTENDLWB , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( UTENDTOT , REAL*8 , (ILG,ILEV) )
AUTOMATIC ( VTENDTOT , REAL*8 , (ILG,ILEV) )
*
AUTOMATIC ( QDvtmp1 , REAL*8 , (ILG ) )
AUTOMATIC ( QDvtmp2 , REAL*8 , (ILG ) )
AUTOMATIC ( QDvtmp3 , REAL*8 , (ILG ) )
*
************************************************************************
*
integer i,l,ii,len,lref,lrefm,jyes,jno
real*8 aux,eta,dz,uparl,piotwo,psi,cpsi,spsi,ratio,
+ fvert,amp,ampmax,vmin,v0,zero,unit,cdblk,
+ hmin1,hmin2
real*8 QDrgas,QDratio, QDgrav,QDhmin2
real*8 QDtmp
*
vmin = 2.
v0 = 1.e-30
hmin1 = 10.
hmin2 = 3.
QDhmin2= 1.0/hmin2
QDrgas = 1.0/dble(rgas)
QDgrav = 1.0/dble(grav)
zero = 0.
unit = 1.
cdblk = 1.
*
len = il2 - il1 + 1
lref = ilev
lrefm = lref - 1
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c INITIAL STEPS
c
c--------------------------------------------------------------------
c Initialize total tendency
c
do l=1,ilev
do i=il1,il2
utend(i,l) = zero
vtend(i,l) = zero
enddo
enddo
c-------------------------------------------------------------------
c Find and gather active grid columns
c
c
jyes = 0
jno = len + 1
c
do i=il1,il2
if ( gc(i).eq.-1. ) then
jyes = jyes + 1
drag(jyes) = i
else
jno = jno - 1
drag(jno) = i
endif
enddo
c
c Check if there is AT LEAST ONE active column
c
if (jyes.le.0) then
goto 600
endif
c
do i=1,len
ii = drag(i) + il1 - 1
env(i) = height(ii)
slp2(i) = slope(ii)
gamma(i) = xcent(ii)
theta(i) = mtdir(ii)
enddo
c
do i=1,len
if (env(i) .lt. hmin2) then
slpf(i) = slp2(i)*Qdhmin2
else
slpf(i) = slp2(i)/env(i)
endif
enddo
c
do l=1,ilev
do i=1,len
ii = drag(i) + il1 - 1
u(i,l) = uu(ii,l)
v(i,l) = vv(ii,l)
tf(i,l) = ttf(ii,l)
th(i,l) = tth(ii,l)
s(i,l) = ss(ii,l)
sh(i,l) = ssh(ii,l)
enddo
enddo
c
call vpown1 (sexpk ,s ,dble(rgocp),len*ilev)
call vpown1 (shexpk,sh,dble(rgocp),len*ilev)
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c GRAVITY-WAVE DRAG
c
c--------------------------------------------------------------------
c Initialize tendencies
do l=1,ilev
do i=1,len
utendgwd(i,l) = zero
vtendgwd(i,l) = zero
enddo
enddo
c
if (gwdrag) then
c
c--------------------------------------------------------------------
c Wind and unit vector at reference level LREFM:
c
do i=1,jyes
vmodb(i) = u(i,lrefm)*u(i,lrefm) + v(i,lrefm)*v(i,lrefm)
enddo
call VSQRT(vmodb,vmodb,jyes)
do i=1,jyes
if (vmodb(i).le.unit) vmodb(i) = unit
QDtmp = 1.d0/vmodb(i)
ub(i) = u(i,lrefm)*QDtmp
vb(i) = v(i,lrefm)*QDtmp
enddo
c
c
c-------------------------------------------------------------------
c Calculate BF frequency:
c
do l=2,ilev
do i=1,jyes
aux = (th(i,l)/shexpk(i,l)-th(i,l-1)/shexpk(i,l-1))
+ /(sh(i,l)-sh(i,l-1))
if ( aux.ge.(-5./s(i,l)) ) aux = -5./s(i,l)
aux = -aux*s(i,l)*grav/(rgas*tf(i,l))
bvfreq(i,l) = grav*aux*sexpk(i,l)/tf(i,l)
enddo
call vsqrt(bvfreq(1,l),bvfreq(1,l),jyes)
enddo
do i=1,jyes
bvfreq(i,1) = bvfreq(i,2)
enddo
c
c * Smooth B V frequency
do l=2,ilev
* ratio = 5.*alog(s(l)/s(l-1))
call VLOG(QDvtmp2,s(1,L), len)
call VLOG(QDvtmp1,s(1,L-1),len)
do i=1,len
ratio = 5.*(QDvtmp2(i) - QDvtmp1(i))
bvfreq(i,l) = (bvfreq(i,l-1) + ratio*bvfreq(i,l))
1 /(1.+ratio)
end do
end do
c
c--------------------------------------------------------------------
c Project wind field on reference wind:
c
do l=1,ilev
do i=1,jyes
veln(i,l) = u(i,l)*ub(i)+v(i,l)*vb(i)
if (veln(i,l).le.v0) veln(i,l) = v0
enddo
enddo
c
c--------------------------------------------------------------------
c Stress field
c
c Compute stress at reference level:
c
do i=1,jyes
asq(i,lref) = env(i)*env(i)
asqs(i,lref) = asq(i,lref)
asqi(i,lref) = asq(i,lref)
dfac(i,lref) = bvfreq(i,lrefm)*sh(i,lrefm)*vmodb(i)
+ /tf(i,lrefm)
enddo
c
c Compute stress at other levels (bottom-up):
c
do l=lrefm,1,-1
do i=1,jyes
dfac(i,l) = bvfreq(i,l)*s(i,l)*veln(i,l)/tf(i,l)
asqi(i,l) = asq(i,l+1)*dfac(i,l+1)/dfac(i,l)
if (veln(i,l).ge.1.) then
QDtmp = veln(i,l)/bvfreq(i,l)
asqs(i,l) = 0.5*Qdtmp*Qdtmp
else
asqs(i,l) = zero
endif
if (asqi(i,l).le.asqs(i,l)) then
asq(i,l) = asqi(i,l)
else
asq(i,l) = asqs(i,l)
endif
depfac(i,l) = (taufac*grav*QDrgas)*asq(i,l)*dfac(i,l)
enddo
enddo
do i=1,jyes
depfac(i,lref) = depfac(i,lrefm)
enddo
c
c--------------------------------------------------------------------
c Compute gwd tendencies:
c
call VREC(QDvtmp1,veln(1,1),jyes)
do i=1,jyes
if ((depfac(i,2) - depfac(i,1)).gt.1.e-10) then
eta = 1.
else
eta = 0.
endif
c
QDtmp = 2.*sh(i,1) + eta*3.*tau*depfac(i,1)*QDvtmp1(i)
QDtmp = 1.0/QDtmp
grad(i,1) = 2.*eta*depfac(i,1)
+ *QDtmp
utendgwd(i,1) = -ub(i)*grad(i,1)
vtendgwd(i,1) = -vb(i)*grad(i,1)
denfac(i,1) = grad(i,1)*3.*tau*depfac(i,1)*QDvtmp1(i)
utendgwd(i,lref) = zero
vtendgwd(i,lref) = zero
enddo
c
do l=2,lrefm
call VREC(QDvtmp1,veln(1,l),jyes)
do i=1,jyes
if ((depfac(i,l) - depfac(i,l-1)).gt.1.e-10) then
eta = 1.
else
eta = 0.
endif
QDtmp = 2.*(sh(i,l)-sh(i,l-1)) +
+ eta*3.*tau*depfac(i,l)*QDvtmp1(i)
QDtmp = 1.0/QDtmp
grad(i,l) = ( 2.*depfac(i,l)-2.*depfac(i,l-1) +
+ eta*denfac(i,l-1) )*QDtmp
utendgwd(i,l) = -ub(i)*grad(i,l)
vtendgwd(i,l) = -vb(i)*grad(i,l)
denfac(i,l) = grad(i,l)*3.*tau*depfac(i,l)*QDvtmp1(i)
enddo
enddo
c
do l=1,ilev
do i=1,jyes
if ( vmodb(i).lt.vmin .or. env(i).lt.hmin1 ) then
utendgwd(i,l) = zero
vtendgwd(i,l) = zero
endif
enddo
enddo
c
endif
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c LOW-LEVEL BLOCKING
c
c--------------------------------------------------------------------
c Initialize arrays
do l=1,ilev
do i=1,len
zb(i,l) = zero
utendllb(i,l) = zero
vtendllb(i,l) = zero
enddo
enddo
c
do i=1,len
izt1(i) = ilev-1
izt2(i) = ilev
izt3(i) = ilev
hblk(i) = zero
izb(i) = ilev
uav(i) = v0
vav(i) = v0
velav(i) = v0
delz(i) = zero
fdir(i) = unit
bloff(i) = 0
enddo
c
if (blocking) then
c--------------------------------------------------------------------
C Recalculate non-smoothed buoyancy frequency
c
C
do l=2,ilev
do i=1,jyes
aux = (th(i,l)/shexpk(i,l)-th(i,l-1)/shexpk(i,l-1))
+ /(sh(i,l)-sh(i,l-1))
if ( aux.ge.(-5./s(i,l)) ) aux = -5./s(i,l)
aux = -aux*s(i,l)*grav/(rgas*tf(i,l))
bvfreq(i,l) = grav*aux*sexpk(i,l)/tf(i,l)
enddo
call vsqrt(bvfreq(1,l),bvfreq(1,l),jyes)
enddo
do i=1,jyes
bvfreq(i,1) = bvfreq(i,2)
enddo
C
c Build elevation field:
c
call VLOG(QDvtmp1,sh(1,ilev),jyes)
do i=1,jyes
zb(i,ilev) = -(rgas*QDgrav)*tf(i,ilev)*QDvtmp1(i)
enddo
do l=ilev-1,1,-1
call VLOG(QDvtmp1,sh(1,l+1),jyes)
call VLOG(QDvtmp2,sh(1,l ),jyes)
do i=1,jyes
zb(i,l) = zb(i,l+1) +
+ (rgas*QDgrav)*tf(i,l)*(QDvtmp1(i)-QDvtmp2(i))
enddo
enddo
c
c--------------------------------------------------------------------
c Blocking height
c
c Find maximum blocking level, upper level for averaging and
c lower level for averaging:
do l=ilev-2,1,-1
do i=1,jyes
if (zb(i,l).lt.(1.5*env(i))) izt3(i) = l
if (zb(i,l).lt. env(i) ) izt1(i) = l
enddo
enddo
c
do l=ilev-1,1,-1
do i=1,jyes
if (zb(i,l).lt.(0.5*env(i))) izt2(i) = l
enddo
enddo
c
c Compute averages:
do l=ilev,2,-1
do i=1,jyes
if (l.le.izt2(i) .and. l.ge.izt1(i)) then
dz = zb(i,l-1) - zb(i,l)
delz(i) = delz(i) + dz
uav(i) = uav(i) + dz*u(i,l)
vav(i) = vav(i) + dz*v(i,l)
endif
enddo
enddo
call VREC(QDvtmp1,delz(1),jyes)
do i=1,jyes
uav(i) = uav(i)*QDvtmp1(i)
vav(i) = vav(i)*QDvtmp1(i)
if (abs(vav(i)).lt.v0 .and. abs(uav(i)).lt.v0) then
velav(i) = v0
else
velav(i) = uav(i)*uav(i) + vav(i)*vav(i)
endif
enddo
c
call vsqrt(velav,velav,jyes)
c
c Compute blocking height and blocking level:
c
do l=2,ilev
do i=1,jyes
if (l.ge.izt3(i) .and. bloff(i).eq.0) then
dz = zb(i,l-1) - zb(i,l)
uparl = (u(i,l)*uav(i) + v(i,l)*vav(i))/velav(i)
if (uparl .lt. v0) then
izb(i) = l-1
bloff(i) = 1
else
hblk(i) = hblk(i) + dz*bvfreq(i,l)/uparl
if (hblk(i) .gt. 0.5) then
izb(i) = l-1
bloff(i) = 1
endif
endif
endif
enddo
enddo
c
c--------------------------------------------------------------------
c Compute directional factor:
c
piotwo = .5*acos(-1.)
do i=1,jyes
c
c Angle between mean wind and topography:
if ( abs(vav(i)) .lt. v0 .and. abs(uav(i)) .lt. v0) then
psi = zero
else
psi = theta(i) - atan2(vav(i),uav(i))
if (psi .gt. piotwo ) psi = psi - 2.*piotwo
if (psi .lt. (-piotwo)) psi = psi + 2.*piotwo
endif
cpsi = (cos(psi))*( cos(psi))
spsi = ( sin(psi) )*( sin(psi) )
c
c Directional factor:
amp = cpsi + gamma(i)*spsi
if (amp .lt. 1.e-10) then
amp = zero
else
QDtmp = cpsi + gamma(i)*spsi
QDtmp = 1.0/QDtmp
amp = 2. - ( gamma(i)*cpsi + spsi )*QDtmp
if (amp.lt.zero) amp = zero
endif
fdir(i) = amp*
+ ( (1.-.18*gamma(i) -.04*(gamma(i)*gamma(i)))*cpsi
+ +( .48*gamma(i) +.30*(gamma(i)*gamma(i)))*spsi )
c
enddo
c--------------------------------------------------------------------
c Compute llb tendencies:
c
do l=ilev,1,-1
do i=1,jyes
if ( velav(i).ge.vmin .and.
+ l.gt.izb(i) .and. zb(i,izb(i)).ge.hmin2 ) then
c
c Vertical factor:
fvert = sqrt( (zb(i,izb(i)) - zb(i,l))
+ /(0.5*env(i) + zb(i,l)) )
c
c Drag intensity:
amp = 0.5*cdblk*slpf(i)*fdir(i)*fvert*
+ (u(i,l)*u(i,l) + v(i,l)*v(i,l))
ampmax = sqrt( u(i,l)*u(i,l) + v(i,l)*v(i,l) )/tau
if (amp.gt.ampmax) then
utendllb(i,l) = -u(i,l)/tau
vtendllb(i,l) = -v(i,l)/tau
else
utendllb(i,l) = -0.5*cdblk*slpf(i)*fdir(i)*fvert*
+ u(i,l)*sqrt(u(i,l)*u(i,l) + v(i,l)*v(i,l))
vtendllb(i,l) = -0.5*cdblk*slpf(i)*fdir(i)*fvert*
+ v(i,l)*sqrt(u(i,l)*u(i,l) + v(i,l)*v(i,l))
endif
endif
enddo
enddo
c
endif
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c OROGRAPHIC LIFT
c
c--------------------------------------------------------------------
c Initialize arrays
do l=1,ilev
do i=1,len
utendlft(i,l) = zero
vtendlft(i,l) = zero
enddo
enddo
c
if (orolift) then
c
endif
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c LEE-WAVE BREAKING
c
c--------------------------------------------------------------------
c Initialize arrays
do l=1,ilev
do i=1,len
utendlwb(i,l) = zero
vtendlwb(i,l) = zero
enddo
enddo
c
if (leewave) then
c
endif
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c TOTAL DRAG AND RESULTING WIND FIELD
c
c--------------------------------------------------------------------
do l=1,ilev
do i=1,len
utendtot(i,l) = zero
vtendtot(i,l) = zero
enddo
enddo
c
c Add and scatter tendencies
do l=1,ilev
do i=1,len
utendtot(i,l) = utendgwd(i,l) +
+ utendllb(i,l) +
+ utendlft(i,l) +
+ utendlwb(i,l)
vtendtot(i,l) = vtendgwd(i,l) +
+ vtendllb(i,l) +
+ vtendlft(i,l) +
+ vtendlwb(i,l)
enddo
enddo
c
do l=1,ilev
do i=1,len
ii = drag(i) + il1 - 1
utend(ii,l) = utendtot(i,l)
vtend(ii,l) = vtendtot(i,l)
enddo
enddo
c
c--------------------------------------------------------------------
600 continue
c--------------------------------------------------------------------
c Apply sgo tendency onto U and V:
c
if (applytend) then
c
do l=1,ilev
do i=il1,il2
uu(i,l) = uu(i,l) + tau*utend(i,l)
vv(i,l) = vv(i,l) + tau*vtend(i,l)
enddo
enddo
c
endif
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
return
end