C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_turb.F,v 1.50 2012/03/27 15:48:27 jmc Exp $
C $Name: $
#include "FIZHI_OPTIONS.h"
C-- File fizhi_turb.F:
C-- Contents
C-- o TURBIO
C-- o TRBFLX
C-- o SFCFLX
C-- o PHI
C-- o PSI
C-- o TRBLEN
C-- o TRBITP
C-- o TRBL20
C-- o TRBL25
C-- o TRBDIF
C-- o VTRI0
C-- o VTRI1
C-- o VTRI2
C-- o LINADJ
C-- o ZCSUB
C-- o SEAICE
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
SUBROUTINE TURBIO (im,jm,nlay,istrip,nymd,nhms,bi,bj, qbeg
1 ,ndturb,nltop,ptop, pz, uz, vz, tz, qz, ntracers,ptracers
2 ,plz,plze,dpres,pkht,pkz,ctmt,xxmt,yymt,zetamt,xlmt,khmt,tke
3 ,tgz,fracland,landtype
4 ,tcanopy,ecanopy,tdeep,swetshal,swetroot,swetdeep,snodep,capac
5 ,nchp,nchptot,nchplnd,chfr,chlt,chlon,igrd,ityp,alai,agrn,thkz
6 ,tprof
8 ,duturb, dvturb, dtturb,dqturb,radlwg,st4,dst4,radswg,radswt
9 ,fdifpar,fdirpar,rainlsp,rainconv,snowfall,tempref
1 ,imstturblw,imstturbsw,qliqavelw,qliqavesw,fccavelw,fccavesw
2 ,qqgrid,myid)
C-----------------------------------------------------------------------
C subroutine turbio - model interface routine to trbflx, the turbulence
C parameterization, and tile, the land surface parameterization
C input:
C im - number of points in the longitude direction
C jm - number of points in the latitude direction
C nlay - number of vertical levels
C istrip - number of horizontal points to be handled at a time on
C nymd - year and date integer in YYMMDD format (ie, 790212)
C nhms - date and time integer in HHMMSS format (ie, 123000)
C qbeg - set to true to indicate a cold start for turbulence
C ndturb - turbulence time step integer in HHMMSS format
C nltop - Top level at which to allow turbulence
C ptop - model top pressure - rigid lid assumed - real
C pz - surface pressure minus ptop in mb - real[lon,lat]
C uz - zonal wind in m/sec - real[lon,lat,level]
C vz - meridional wind in m/sec - real[lon,lat,level]
C tz - model theta (theta [deg K]/p0**k) - real[lon,lat,level]
C qz - specific humidity in kg/kg - real[lon,lat,level]
C ntracers- total number of tracers - integer
C ptracers- number of permanent tracers - integer
C pkht - pressure[mb]**k at bottom edges of levels - real[lon,lat,level]
C fracland- not being used - real[lon,lat]
C landtype- not being used - integer[lon,lat]
C nchp - nchplnd
C nchplnd - <=nchp - number of land chips - integer
C chfr - chip fraction - real[nchp]
C chlt - tile space latitude array - real[nchp]
C chlon - tile space longitude array - real[nchp]
C igrd - tile space grid number - integer[nchp]
C ityp - tile space vegetation type - integer[nchp]
C alai - leaf area index - real[nchp]
C agrn - greenness fraction - real[nchp]
C thkz - sea ice thickness in m (0. for no ice) - real[lon,lat]
C tprof - logical flag for point by point diagnostic output
C radlwg - net longwave flux at ground (up-down) in w/m**2 - real[lon,lat]
C st4 - upward longwave flux at ground in w/m**2 - real[lon,lat]
C dst4 - delta-sigma-T**4, ie, derivative of upward lw flux at
C ground with respect to ground Temperature - real[lon,lat]
C radswg - net shortwave flux at ground (down-up) NON-DIM - real[lon,lat]
C {NOTE: this field is divided by the incident shortwave
C at the top of the atmosphere to non-dimensionalize]
C radswt - incident shortwave at top of atmos in W/m**2 - real[lon,lat]
C fdifpar - incident diffuse-beam PAR at surface in W/m**2 - real[lon,lat]
C fdirpar - incident direct-beam PAR at surface in W/m**2 - real[lon,lat]
C rainlsp - large-scale (frontal,supersat) rainfall in mm/sec - real[lon,lat]
C rainconv- convective rainfall rate in mm/sec - real[lon,lat]
C snowfall- total snowfall rate in mm/sec - real[lon,lat]
C updated:
C tke - turbulent k.e. in m**2/s**2 - real[tiles,levels]
C tgz - surface skin temperature in deg K - real[lon,lat]
C tcanopy - canopy temperature in deg K real[tiles]
C (sea surface temp over the ocean tiles)
C ecanopy - canopy vapor pressure in mb real[tiles]
C (qstar at tground over the sea ice and ocean tiles)
C tdeep - deep soil temp in deg K real[tiles]
C swetshal- shallow level moisture field capacity fraction real[tiles]
C swetroot- root level moisture field capacity fraction real[tiles]
C swetdeep- deep soil level moisture field capacity fraction real[tiles]
C snodep - depth of snow pack in cm liquid water equiv real[tiles]
C capac - leaf canopy water reservoir in cm real[tiles]
C output:
C duturb - change in zonal wind component due to turbulent processes
C per unit time in m/sec**2 - real[lon,lat,levels]
C dvturb - change in meridional wind component due to turbulent processes
C per unit time in m/sec**2 - real[lon,lat,levels]
C dtturb - change in (model theta*pi) due to turbulent processes
C per unit time - real[lon,lat,levels] !! pi is pressure-ptop
C dqturb - change in (specific humidity*pi) due to turbulent processes
C per unit time - real[lon,lat,levels] !! pi is pressure-ptop
C qliqavelw - Moist Turbulence Liquid Water for Longwave - real[lon,lat,levels]
C qliqavesw - Moist Turbulence Liquid Water for Shortwave - real[lon,lat,levels]
C fccavelw - Moist Turbulence Cloud Fraction for Longwave - real[lon,lat,levels]
C fccavesw - Moist Turbulence Cloud Fraction for Shortwave - real[lon,lat,levels]
C qqgrid - Gridded Turbulent Kinetic Energy - real[lon,lat,levels]
C-----------------------------------------------------------------------
implicit none
integer im,jm,nlay,istrip,nymd,nhms,bi,bj,ndturb,nltop
integer ntracers, ptracers
integer nchp,nchptot,nchplnd
logical qbeg
_RL ptop
_RL pz(im*jm,1),uz(im*jm,1,nlay),vz(im*jm,1,nlay)
_RL tz(im*jm,1,nlay),qz(im*jm,1,nlay,ntracers)
_RL plz(im*jm,1,nlay),plze(im*jm,1,nlay+1),dpres(im*jm,1,nlay)
_RL pkht(im*jm,1,nlay+1),pkz(im*jm,1,nlay)
_RL ctmt(nchp),xxmt(nchp),yymt(nchp),zetamt(nchp)
_RL xlmt(nchp,nlay),khmt(nchp,nlay),tke(nchp,nlay)
_RL tgz(im*jm,1),fracland(im*jm,1)
integer landtype(im*jm,1)
_RL tcanopy(nchp),tdeep(nchp),ecanopy(nchp),swetshal(nchp)
_RL swetroot(nchp),swetdeep(nchp),snodep(nchp),capac(nchp)
_RL chfr(nchp),chlt(nchp),chlon(nchp)
integer igrd(nchp),ityp(nchp)
_RL alai(nchp),agrn(nchp),thkz(im*jm,1)
logical tprof
_RL duturb(im*jm,1,nlay),dvturb(im*jm,1,nlay)
_RL dtturb(im*jm,1,nlay),dqturb(im*jm,1,nlay,ntracers)
_RL st4(im*jm,1),dst4(im*jm,1)
_RL radswg(im*jm,1),radswt(im*jm,1),radlwg(im*jm,1)
_RL fdifpar(im*jm,1),fdirpar(im*jm,1)
_RL rainlsp(im*jm,1),rainconv(im*jm,1),snowfall(im*jm,1)
_RL tempref (im*jm,1)
integer imstturblw, imstturbsw
_RL qliqavesw(im*jm,1,nlay),qliqavelw(im*jm,1,nlay)
_RL fccavelw (im*jm,1,nlay),fccavesw (im*jm,1,nlay)
_RL qqgrid (im*jm,1,nlay)
integer myid
C Local Variables
integer numstrips
integer ijall
_RL fmu,hice,tref,pref,cti,ed
C Set fmu and ed to zero for no background diffusion
parameter ( fmu = 0.00000 )
parameter ( hice = 300. )
parameter ( tref = 258. )
parameter ( pref = 500. )
parameter ( cti = 0.0052 )
parameter ( ed = 0.0 )
_RL qliqtmp(im*jm,1,nlay)
_RL fcctmp(im*jm,1,nlay)
_RL tmpdiag(im*jm,1), tmpFac
_RL thtgz(im*jm)
logical diagnostics_is_on
external
integer nland
_RL alwcoeff(nchp),blwcoeff(nchp)
_RL netsw(nchp)
_RL cnvprec(nchp),lsprec(nchp)
_RL snowprec(nchp)
_RL pardiff(nchp),pardirct(nchp)
_RL pmsc(nchp)
_RL netlw(nchp)
_RL sqscat(nchp), rsoil1(nchp)
_RL rsoil2(nchp)
_RL rdc(nchp),u2fac(nchp)
_RL z2ch(nchp)
_RL zoch(nchp),cdrc(nchp)
_RL cdsc(nchp)
_RL dqsdt(nchp)
_RL tground(nchp),qground(nchp)
_RL utility(nchp)
_RL qice(nchp)
_RL dqice(nchp)
_RL dumsc(nchp,nlay),dvmsc(nchp,nlay)
_RL dtmsc(nchp,nlay),dqmsc(nchp,nlay,ntracers)
_RL shg(nchp),z0(nchp),icethk(nchp)
integer water(nchp)
_RL lats(istrip),lons(istrip),cosz(istrip),icest(istrip)
_RL rainls(istrip),raincon(istrip),newsnow(istrip)
_RL pardf(istrip),pardr(istrip),swnet(istrip)
_RL hlwdwn(istrip),alwrad(istrip),blwrad(istrip)
_RL tmpnlay(istrip)
_RL laistrip(istrip),grnstrip(istrip),z2str(istrip),cd(istrip)
_RL scatstr(istrip), rs1str(istrip), rs2str(istrip)
_RL rdcstr(istrip),u2fstr(istrip),dqsdtstr(istrip)
_RL eturb(istrip),dedqa(istrip),dedtc(istrip)
_RL hsturb(istrip),dhsdqa(istrip),dhsdtc(istrip)
_RL savetc(istrip),saveqa(istrip),lwstrip(istrip)
_RL chfrstr(istrip),psurf(istrip),shgstr(istrip)
integer types(istrip),igrdstr(istrip)
_RL evap(istrip),shflux(istrip),runoff(istrip),bomb(istrip)
_RL eint(istrip),esoi(istrip),eveg(istrip),esno(istrip)
_RL smelt(istrip),hlatn(istrip),hlwup(istrip),gdrain(istrip)
_RL runsrf(istrip),fwsoil(istrip),evpot(istrip)
_RL strdg1(istrip),strdg2(istrip),strdg3(istrip),strdg4(istrip)
_RL strdg5(istrip),strdg6(istrip),strdg7(istrip),strdg8(istrip)
_RL strdg9(istrip),tmpstrip(istrip),qicestr(istrip)
_RL dqicestr(istrip)
_RL u(istrip,nlay+1), v(istrip,nlay+1), th(istrip,nlay+1)
_RL sh(istrip,nlay+1), thv(istrip,nlay+1), pe(istrip,nlay+1)
_RL tracers(istrip,nlay+1,ntracers)
_RL dpstr(istrip,nlay),pke(istrip,nlay+1)
_RL pk(istrip,nlay), qq(istrip,nlay), p(istrip,nlay)
_RL sri(istrip,nlay), skh(istrip,nlay), skm(istrip,nlay)
_RL stuflux(istrip,nlay), stvflux(istrip,nlay)
_RL sttflux(istrip,nlay), stqflux(istrip,nlay)
_RL frqtrb(istrip,nlay-1)
_RL dshdthg(istrip,nlay),dthdthg(istrip,nlay)
_RL dshdshg(istrip,nlay),dthdshg(istrip,nlay)
_RL transth(istrip,nlay), transsh(istrip,nlay)
_RL tc(istrip),td(istrip),qa(istrip)
_RL swet1(istrip),swet2(istrip),swet3(istrip)
_RL capacity(istrip),snowdepth(istrip)
_RL stz0(istrip)
_RL stdiag(istrip)
_RL tends(istrip),sustar(istrip), sz0(istrip),pbldpth(istrip)
_RL sct(istrip), scu(istrip), swinds(istrip)
_RL stu2m(istrip),stv2m(istrip),stt2m(istrip),stq2m(istrip)
_RL stu10m(istrip),stv10m(istrip),stt10m(istrip),stq10m(istrip)
integer stwatr(istrip)
_RL wspeed(istrip)
_RL ctsave(istrip),xxsave(istrip),yysave(istrip)
_RL zetasave(istrip)
_RL xlsave(istrip,nlay),khsave(istrip,nlay)
_RL qliq(istrip,nlay),turbfcc(istrip,nlay)
_RL qliqmsc(nchp,nlay),fccmsc(nchp,nlay)
_RL pi,secday,sdayopi2,rgas,akap,cp,alhl
_RL faceps,grav,caltoj,virtcon,getcon
_RL heatw,undef,timstp,delttrb,dttrb,ra
_RL edle,rmu,cltj10,atimstp,tice,const
integer istnp1,istnlay,itrtrb,i,L,nn,nt
integer nocean, nice
integer ndmoist,time_left, ndum0,ndum1
integer ntracedim
_RL dtfac,timstp2
integer n,nsecf,nmonf,ndayf
nsecf(n) = n/10000*3600 + mod(n,10000)/100* 60 + mod(n,100)
nmonf(n) = mod(n,10000)/100
ndayf(n) = mod(n,100)
#ifdef CRAY
#ifdef f77
cfpp$ expand (qsat)
#endif
#endif
C compute variables that do not change
pi = 4.*atan(1.)
secday = getcon('SDAY')
sdayopi2 = getcon('SDAY') / (pi*2.)
rgas = getcon('RGAS')
akap = getcon('KAPPA')
cp = getcon('CP')
alhl = getcon('LATENT HEAT COND')
faceps = getcon('EPSFAC')
grav = getcon('GRAVITY')
caltoj = getcon('CALTOJ')
virtcon = getcon('VIRTCON')
heatw = getcon('HEATW')
undef = getcon('UNDEF')
ntracedim= max(ntracers-ptracers,1)
call GET_ALARM ( 'moist',ndum0,ndum1,ndmoist,time_left )
timstp = nsecf(ndturb)
timstp2 = nsecf(ndmoist)
dtfac = min( 1.0 _d 0, timstp/timstp2 )
C delttrb is the internal turbulence time step
C a value equal to ndturb means one internal iteration
delttrb = nsecf(ndturb)
ijall = im * jm
istnp1 = istrip * (nlay+1)
istnlay = istrip * nlay
itrtrb = ( timstp / delttrb ) + 0.1
dttrb = timstp / float(itrtrb)
edle = ed * 0.2
C coefficient of viscosity (background momentum diffusion)
rmu = fmu * tref * rgas / pref
cltj10 = 10. * caltoj
atimstp = 1. / timstp
tice = getcon('FREEZING-POINT')
C **********************************************************************
C Initialization
C **********************************************************************
C Initialize diagnostic for ground temperature change
C ---------------------------------------------------
tmpFac = -1.
CALL DIAGNOSTICS_SCALE_FILL( tgz,tmpFac,1,
& 'DTG ',0,1,-3,bi,bj,myid)
C **********************************************************************
C entire turbulence and land surface package will run in 'tile space'
C do conversion of model state variables to tile space
C (ocean points appended to tile space land point arrays)
C **********************************************************************
numstrips = ((nchptot-1)/istrip) + 1
call GRD2MSC(pz(1,1),im,jm,igrd,pmsc,nchp,nchptot)
call GRD2MSC(tgz,im,jm,igrd,tground,nchp,nchptot)
do i = 1,ijall
tmpdiag(i,1) = st4(i,1) + dst4(i,1)*(tgz(i,1)-tempref(i,1))
1 - dst4(i,1)* tgz(i,1)
enddo
call GRD2MSC(tmpdiag,im,jm,igrd,alwcoeff,nchp,nchptot)
do i = 1,ijall
tmpdiag(i,1) = dst4(i,1)
enddo
call GRD2MSC(tmpdiag,im,jm,igrd,blwcoeff,nchp,nchptot)
do i = 1,ijall
tmpdiag(i,1) = fdifpar(i,1) * radswt(i,1)
enddo
call GRD2MSC(tmpdiag,im,jm,igrd,pardiff,nchp,nchptot)
do i = 1,ijall
tmpdiag(i,1) = fdirpar(i,1) * radswt(i,1)
enddo
call GRD2MSC(tmpdiag,im,jm,igrd,pardirct,nchp,nchptot)
do i = 1,ijall
tmpdiag(i,1) = radswg(i,1) * radswt(i,1)
enddo
call GRD2MSC(tmpdiag,im,jm,igrd,netsw,nchp,nchptot)
do i = 1,ijall
tmpdiag(i,1) = radlwg(i,1) + dst4(i,1)*(tgz(i,1)-tempref(i,1))
enddo
call GRD2MSC(tmpdiag,im,jm,igrd,netlw,nchp,nchptot)
call GRD2MSC(thkz,im,jm,igrd,icethk,nchp,nchptot)
call GRD2MSC(rainlsp,im,jm,igrd,lsprec,nchp,nchptot)
call GRD2MSC(rainconv,im,jm,igrd,cnvprec,nchp,nchptot)
call GRD2MSC(snowfall,im,jm,igrd,snowprec,nchp,nchptot)
C Call chpprm to get non-varying vegetation and soil characteristics
call CHPPRM(nymd,nhms,nchp,nchplnd,chlt,ityp,alai,
1 agrn,zoch,z2ch,cdrc,cdsc,sqscat,u2fac,rsoil1,rsoil2,rdc)
C **********************************************************************
C **** surface specification ****
C **********************************************************************
C set water
do i = 1,nchptot
water(i) = 0
if((ityp(i).eq.100).and.(icethk(i).eq.0. ))water(i) = 1
enddo
C roughness length z0
do i =1,nchptot
if (icethk(i).gt.0.) then
z0(i) = 1.e-4
else if (ityp(i).eq.100) then
z0(i) = 3.e-4
else
z0(i) = zoch(i)
endif
enddo
C Fill Array Tground with canopy temperatures over land tiles
C (it has sst from the tgz array over the sea ice and ocean tiles)
do i = 1,nchplnd
tground(i) = tcanopy(i)
enddo
C value of sh at ground
C ---------------------
do I =1,nchptot
utility(I) = pmsc(i) + ptop
call QSAT ( tground(i),utility(i),shg(i),dqsdt(i),.true. )
enddo
C Fill Array Qground with canopy air specific humidity over land tiles
C (it has qstar at tground over the sea ice and ocean tiles)
do i = 1,nchplnd
qground(i) = ecanopy(i)
enddo
do i = nchplnd+1,nchptot
qground(i) = shg(i)
enddo
C Fill Array Swetshal with Value 1. over oceans and sea ice
do i = nchplnd+1,nchptot
swetshal(i) = 1.
enddo
C compute heat conduction through ice
C -----------------------------------
const = ( cti / hice ) * cltj10
do i =1,nchptot
qice(i) = 0.0
dqice(i) = 0.0
if( icethk(i).gt.0.0 ) then
qice(i) = const*(tice-tground(i))
dqice(i) = -const
endif
enddo
if(diagnostics_is_on('QICE ',myid) ) then
do i =1,ijall
tmpdiag(i,1) = 0.0
enddo
call MSC2GRD (igrd,chfr,qice,nchp,nchptot,fracland,tmpdiag,im,jm)
call DIAGNOSTICS_FILL(tmpdiag,'QICE ',0,1,3,bi,bj,myid)
endif
C***********************************************************************
C loop over regions
C***********************************************************************
do 2000 nn = 1, numstrips
call STRIP2TILE(uz,igrd,u,nchp,ijall,istrip,nlay,nn)
call STRIP2TILE(vz,igrd,v,nchp,ijall,istrip,nlay,nn)
call STRIP2TILE(tz,igrd,th,nchp,ijall,istrip,nlay,nn)
call STRIP2TILE(qz(1,1,1,1),igrd,sh,nchp,ijall,istrip,nlay,nn)
call STRIP2TILE(dpres,igrd,dpstr,nchp,ijall,istrip,nlay,nn)
call STRIP2TILE(plz,igrd,p,nchp,ijall,istrip,nlay,nn)
call STRIP2TILE(plze,igrd,pe,nchp,ijall,istrip,nlay+1,nn)
call STRIP2TILE(pkz,igrd,pk,nchp,ijall,istrip,nlay,nn)
call STRIP2TILE(pkht,igrd,pke,nchp,ijall,istrip,nlay+1,nn)
c do nt = 1,ntracers-ptracers
c call strip2tile(qz(1,1,1,ptracers+nt),igrd,tracers(1,1,nt),nchp,
c 1 ijall,istrip,nlay,nn)
c enddo
call STRIPIT (z0,stz0,nchptot,nchp,istrip,1,nn)
call STRIPIT (tground,th(1,nlay+1),nchptot,nchp,istrip,1,nn)
call STRIPIT (pmsc,pe(1,nlay+1),nchptot,nchp,istrip,1,nn)
call STRIPIT (tke,qq,nchptot,nchp,istrip,nlay-1,nn)
call STRIPIT (ctmt,ctsave,nchptot,nchp,istrip,1,nn)
call STRIPIT (xxmt,xxsave,nchptot,nchp,istrip,1,nn)
call STRIPIT (yymt,yysave,nchptot,nchp,istrip,1,nn)
call STRIPIT (zetamt,zetasave,nchptot,nchp,istrip,1,nn)
call STRIPIT (xlmt,xlsave,nchptot,nchp,istrip,nlay,nn)
call STRIPIT (khmt,khsave,nchptot,nchp,istrip,nlay,nn)
call STRIPITINT (water,stwatr,nchptot,nchp,istrip,1,nn)
call STRIPITINT (igrd,igrdstr,nchptot,nchp,istrip,1,nn)
call STRIPIT (chfr,chfrstr,nchptot,nchp,istrip,1,nn)
call STRIPIT (icethk,icest,nchptot,nchp,istrip,1,nn)
call STRIPIT (pardiff,pardf,nchptot,nchp,istrip,1,nn)
call STRIPIT (pardirct,pardr,nchptot,nchp,istrip,1,nn)
call STRIPIT (chlt,lats,nchptot,nchp,istrip,1,nn)
call STRIPIT (chlon,lons,nchptot,nchp,istrip,1,nn)
call STRIPIT (lsprec,rainls,nchptot,nchp,istrip,1,nn)
call STRIPIT (cnvprec,raincon,nchptot,nchp,istrip,1,nn)
call STRIPIT (snowprec,newsnow,nchptot,nchp,istrip,1,nn)
call STRIPIT (netsw,swnet,nchptot,nchp,istrip,1,nn)
call STRIPIT (netlw,lwstrip,nchptot,nchp,istrip,1,nn)
call STRIPIT (alwcoeff,alwrad,nchptot,nchp,istrip,1,nn)
call STRIPIT (blwcoeff,blwrad,nchptot,nchp,istrip,1,nn)
call STRIPIT (alai,laistrip,nchptot,nchp,istrip,1,nn)
call STRIPIT (agrn,grnstrip,nchptot,nchp,istrip,1,nn)
call STRIPIT (z2ch,z2str,nchptot,nchp,istrip,1,nn)
call STRIPIT (sqscat,scatstr,nchptot,nchp,istrip,1,nn)
call STRIPIT (rsoil1,rs1str,nchptot,nchp,istrip,1,nn)
call STRIPIT (rsoil2,rs2str,nchptot,nchp,istrip,1,nn)
call STRIPIT (rdc,rdcstr,nchptot,nchp,istrip,1,nn)
call STRIPIT (u2fac,u2fstr,nchptot,nchp,istrip,1,nn)
call STRIPIT (shg,shgstr,nchptot,nchp,istrip,1,nn)
call STRIPIT (dqsdt,dqsdtstr,nchptot,nchp,istrip,1,nn)
call STRIPIT ( qice, qicestr,nchptot,nchp,istrip,1,nn)
call STRIPIT (dqice,dqicestr,nchptot,nchp,istrip,1,nn)
call STRIPITINT (ityp,types,nchptot,nchp,istrip,1,nn)
call STRIPIT (tground,tc,nchptot,nchp,istrip,1,nn)
call STRIPIT (tdeep,td,nchptot,nchp,istrip,1,nn)
call STRIPIT (qground,qa,nchptot,nchp,istrip,1,nn)
call STRIPIT (swetshal,swet1,nchptot,nchp,istrip,1,nn)
call STRIPIT (swetroot,swet2,nchptot,nchp,istrip,1,nn)
call STRIPIT (swetdeep,swet3,nchptot,nchp,istrip,1,nn)
call STRIPIT (snodep,snowdepth,nchptot,nchp,istrip,1,nn)
call STRIPIT (capac,capacity,nchptot,nchp,istrip,1,nn)
#ifdef FIZHI_USE_FIXED_DAY
call ASTRO ( 20040321,nhms,lats,lons,istrip,cosz,ra )
#else
call ASTRO ( nymd,nhms,lats,lons,istrip,cosz,ra )
#endif
C we need to count up the land, sea ice and ocean points
nocean = 0
nland = 0
nice = 0
do i = 1,istrip
if( types(i).lt.100 ) nland = nland + 1
if( types(i).eq.100 ) nocean = nocean + 1
if( types(i).eq.100 .and. icest(i).gt.0.0 ) nice = nice + 1
enddo
C Zero out velocities at the bottom edge of the model
C ---------------------------------------------------
do i =1,istrip
u(i,nlay+1) = 0.
v(i,nlay+1) = 0.
enddo
C convert temperature of level nlay+1 to theta & value of sh at ground
C --------------------------------------------------------------------
do i =1,istrip
th(i,nlay+1) = th(i,nlay+1) / pke(i,nlay+1)
sh(i,nlay+1) = qa(i)
enddo
if(diagnostics_is_on('QG ',myid) ) then
do i=1,istrip
tmpstrip(i) = sh(i,nlay+1)*1000
enddo
call DIAG_VEGTILE_FILL (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
& .false., 'QG ', 1, 1, bi, bj, myid)
endif
C value of tracers at the ground
C ------------------------------
c do nt = 1,ntracers-ptracers
c do i = 1,istrip
c tracers(i,nlay+1,nt) = 0.
c enddo
c enddo
C compute virtual potential temperatures
C --------------------------------------
do L = 1,nlay+1
do i =1,istrip
thv(i,L) = 1. + virtcon * sh(i,L)
thv(i,L) = th(i,L) * thv(i,L)
enddo
enddo
do i =1,istrip
sh(i,nlay+1) = qa(i)
enddo
C zero out arrays for output of qliq and fcc
do L =1,nlay
do i =1,istrip
qliq(i,L) = 0.
turbfcc(i,L) = 0.
enddo
enddo
C zero out fluxes and derivatives
C -------------------------------
do i = 1,istrip
eturb(i) = 0.
scu(i) = 0.
dedqa(i) = 0.
dedtc(i) = 0.
hsturb(i) = 0.
dhsdqa(i) = 0.
dhsdtc(i) = 0.
enddo
C increment diagnostic arrays for quantities calculated before trbfl
C ------------------------------------------------------------------
if(diagnostics_is_on('DTSRF ',myid) ) then
do i=1,istrip
stdiag(i) = ( thv(i,nlay+1)-thv(i,nlay) ) / pke(i,nlay+1)
enddo
call DIAG_VEGTILE_FILL (stdiag,igrd,chfrstr,istrip,nchp,nn,
& .false., 'DTSRF ', 1, 1, bi, bj, myid)
endif
C call trbflx
C -----------
call TRBFLX(nn,th,thv,sh,u,v,qq,p,pe,pk,pke,dpstr,stwatr,stz0,
1 tracers,ntracers-ptracers,ntracedim,dttrb,itrtrb,rmu,edle,qbeg,
2 tprof,stuflux,stvflux,sri,skh,skm,swinds,sustar,sz0,frqtrb,
3 pbldpth,sct,scu,stu2m,stv2m,stt2m,stq2m,stu10m,stv10m,stt10m,
4 stq10m,istrip,nlay,nltop,nymd,nhms,grav,cp,rgas,faceps,virtcon,
5 undef,dshdthg,dshdshg,dthdthg,dthdshg,eturb,dedqa,dedtc,
6 hsturb,dhsdqa,dhsdtc,transth,transsh,
7 ctsave,xxsave,yysave,zetasave,xlsave,khsave,qliq,turbfcc)
call PASTIT (qq,tke,istrip,nchp,nchptot,nlay,nn)
call PASTIT (ctsave,ctmt,istrip,nchp,nchptot,1,nn)
call PASTIT (xxsave,xxmt,istrip,nchp,nchptot,1,nn)
call PASTIT (yysave,yymt,istrip,nchp,nchptot,1,nn)
call PASTIT (zetasave,zetamt,istrip,nchp,nchptot,1,nn)
call PASTIT (xlsave,xlmt,istrip,nchp,nchptot,nlay,nn)
call PASTIT (khsave,khmt,istrip,nchp,nchptot,nlay,nn)
call PASTIT (qliq ,qliqmsc,istrip,nchp,nchptot,nlay,nn)
call PASTIT (turbfcc,fccmsc,istrip,nchp,nchptot,nlay,nn)
C New diagnostic: potential evapotranspiration
do i = 1,istrip
evpot(i) = transsh(i,nlay) * (shgstr(i) - sh(i,nlay))
enddo
C**********************************************************************
C Call Land Surface Module
C**********************************************************************
do i = 1,istrip
savetc(i) = tc(i)
saveqa(i) = qa(i)
enddo
do i = 1,istrip
cosz(i) = max(cosz(i),0.0001 _d 0)
cd(i) = scu(i)*scu(i)
tmpnlay(i) = th(i,nlay)*pk(i,nlay)
hlwdwn(i) = alwrad(i)+blwrad(i)*tc(i)-lwstrip(i)
psurf(i) = pe(i,nlay+1)
wspeed(i) = sqrt(u(i,nlay)*u(i,nlay) + v(i,nlay)*v(i,nlay))
if(wspeed(i) .lt. 1.e-20) wspeed(i) = 1.e-20
C Note: This LSM precip bug needs to be cleaned up
c newsnow(i) = newsnow(i)*dtfac
c raincon(i) = raincon(i)*dtfac
c rainls (i) = rainls (i)*dtfac
enddo
do i = 1,istrip
eturb(i) = eturb(i) * pke(i,nlay+1)
dedqa(i) = dedqa(i) * pke(i,nlay+1)
hsturb(i) = hsturb(i) * pke(i,nlay+1)
enddo
do i = 1,istrip
strdg1(i) = 0.
strdg2(i) = 0.
strdg3(i) = 0.
strdg4(i) = 0.
strdg5(i) = 0.
strdg6(i) = 0.
strdg7(i) = 0.
strdg8(i) = 0.
strdg9(i) = 0.
bomb(i) = 0.
runoff(i) = 0.
eint(i) = 0.
esoi(i) = 0.
eveg(i) = 0.
esno(i) = 0.
smelt(i) = 0.
hlatn(i) = 0.
hlwup(i) = 0.
gdrain(i) = 0.
runsrf(i) = 0.
fwsoil(i) = 0.
enddo
C**********************************************************************
C diagnostics: fill arrays for lsm input fields
C**********************************************************************
if(diagnostics_is_on('SNOWFALL',myid) ) then
do i = 1,istrip
tmpstrip(i) = newsnow(i)*86400
enddo
call DIAG_VEGTILE_FILL (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
& .false., 'SNOWFALL', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('RAINCON ',myid) ) then
do i = 1,istrip
tmpstrip(i) = raincon(i)*86400
enddo
call DIAG_VEGTILE_FILL (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
& .false., 'RAINCON ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('RAINLSP ',myid) ) then
do i = 1,istrip
tmpstrip(i) = rainls(i)*86400
enddo
call DIAG_VEGTILE_FILL (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
& .false., 'RAINLSP ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('GREEN ',myid) ) then
call DIAG_VEGTILE_FILL (grnstrip,igrd,chfrstr,istrip,nchp,nn,
& .false., 'GREEN ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('LAI ',myid) ) then
call DIAG_VEGTILE_FILL (laistrip,igrd,chfrstr,istrip,nchp,nn,
& .false., 'LAI ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('PARDR ',myid) ) then
call DIAG_VEGTILE_FILL (pardr,igrd,chfrstr,istrip,nchp,nn,
& .false., 'PARDR ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('PARDF ',myid) ) then
call DIAG_VEGTILE_FILL (pardf,igrd,chfrstr,istrip,nchp,nn,
& .false., 'PARDF ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('DLWDTC ',myid) ) then
call DIAG_VEGTILE_FILL (blwrad,igrd,chfrstr,istrip,nchp,nn,
& .false., 'DLWDTC ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('DHDTC ',myid) ) then
call DIAG_VEGTILE_FILL (dhsdtc,igrd,chfrstr,istrip,nchp,nn,
& .false., 'DHDTC ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('DEDTC ',myid) ) then
call DIAG_VEGTILE_FILL (dedtc,igrd,chfrstr,istrip,nchp,nn,
& .false., 'DEDTC ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('DHDQA ',myid) ) then
call DIAG_VEGTILE_FILL (dhsdqa,igrd,chfrstr,istrip,nchp,nn,
& .false., 'DHDQA ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('DEDQA ',myid) ) then
call DIAG_VEGTILE_FILL (dedqa,igrd,chfrstr,istrip,nchp,nn,
& .false., 'DEDQA ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('LWGDOWN ',myid) ) then
call DIAG_VEGTILE_FILL (hlwdwn,igrd,chfrstr,istrip,nchp,nn,
& .false., 'LWGDOWN ', 1, 1, bi, bj, myid)
endif
C**********************************************************************
if(nland.gt.0)then
call TILE (
I nland, timstp, types, rainls, raincon, newsnow, wspeed,
I eturb, dedqa, dedtc, hsturb, dhsdqa, dhsdtc,
I tmpnlay, sh(1,nlay), cd, cosz, pardr, pardf,
I swnet, hlwdwn, psurf, laistrip, grnstrip, z2str,
I scatstr, rs1str, rs2str, rdcstr, u2fstr,
I shgstr, dqsdtstr, alwrad, blwrad,
U tc, td, qa, swet1, swet2, swet3, capacity, snowdepth,
O evap, shflux, runoff, bomb,
O eint, esoi, eveg, esno, smelt, hlatn,
O hlwup, gdrain, runsrf, fwsoil,
O strdg1, strdg2, strdg3, strdg4,
O strdg5, strdg6, strdg7, strdg8, strdg9)
endif
if( nice.gt.0 ) then
call SEAICE ( nocean, timstp, hice,
& eturb(nland+1), dedtc(nland+1),
& hsturb(nland+1), dhsdtc(nland+1),
& qicestr(nland+1), dqicestr(nland+1),
& swnet(nland+1), lwstrip(nland+1), blwrad(nland+1),
& pke(nland+1,nlay+1), icest(nland+1),
& tc(nland+1), qa(nland+1) )
endif
C***********************************************************************
C diagnostics: fill arrays for lsm output fields
C***********************************************************************
if(diagnostics_is_on('RUNOFF ',myid) ) then
call DIAG_VEGTILE_FILL (runoff,igrd,chfrstr,istrip,nchp,nn,
& .false., 'RUNOFF ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('FWSOIL ',myid) ) then
call DIAG_VEGTILE_FILL (fwsoil,igrd,chfrstr,istrip,nchp,nn,
& .false., 'FWSOIL ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('GDRAIN ',myid) ) then
call DIAG_VEGTILE_FILL (gdrain,igrd,chfrstr,istrip,nchp,nn,
& .false., 'GDRAIN ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('SNOWMELT',myid) ) then
call DIAG_VEGTILE_FILL (smelt,igrd,chfrstr,istrip,nchp,nn,
& .false., 'SNOWMELT', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('EVEG ',myid) ) then
call DIAG_VEGTILE_FILL (eveg,igrd,chfrstr,istrip,nchp,nn,
& .false., 'EVEG ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('ESNOW ',myid) ) then
call DIAG_VEGTILE_FILL (esno,igrd,chfrstr,istrip,nchp,nn,
& .false., 'ESNOW ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('ESOIL ',myid) ) then
call DIAG_VEGTILE_FILL (esoi,igrd,chfrstr,istrip,nchp,nn,
& .false., 'ESOIL ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('ERESV ',myid) ) then
call DIAG_VEGTILE_FILL (eint,igrd,chfrstr,istrip,nchp,nn,
& .false., 'ERESV ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('EVPOT ',myid) ) then
call DIAG_VEGTILE_FILL (evpot,igrd,chfrstr,istrip,nchp,nn,
& .false., 'EVPOT ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('DTC ',myid) ) then
call DIAG_VEGTILE_FILL (strdg1,igrd,chfrstr,istrip,nchp,nn,
& .false., 'DTC ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('DQC ',myid) ) then
call DIAG_VEGTILE_FILL (strdg2,igrd,chfrstr,istrip,nchp,nn,
& .false., 'DQC ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('TCDTC ',myid) ) then
call DIAG_VEGTILE_FILL (strdg3,igrd,chfrstr,istrip,nchp,nn,
& .false., 'TCDTC ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('RADDTC ',myid) ) then
call DIAG_VEGTILE_FILL (strdg4,igrd,chfrstr,istrip,nchp,nn,
& .false., 'RADDTC ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('SENSDTC ',myid) ) then
call DIAG_VEGTILE_FILL (strdg5,igrd,chfrstr,istrip,nchp,nn,
& .false., 'SENSDTC ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('LATDTC ',myid) ) then
call DIAG_VEGTILE_FILL (strdg6,igrd,chfrstr,istrip,nchp,nn,
& .false., 'LATDTC ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('TDDTC ',myid) ) then
call DIAG_VEGTILE_FILL (strdg7,igrd,chfrstr,istrip,nchp,nn,
& .false., 'TDDTC ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('QCDTC ',myid) ) then
call DIAG_VEGTILE_FILL (strdg8,igrd,chfrstr,istrip,nchp,nn,
& .false., 'QCDTC ', 1, 1, bi, bj, myid)
endif
C***********************************************************************
call PASTIT (tc,tground,istrip,nchp,nchptot,1,nn)
call PASTIT (td,tdeep,istrip,nchp,nchptot,1,nn)
call PASTIT (qa,qground,istrip,nchp,nchptot,1,nn)
call PASTIT (swet1,swetshal,istrip,nchp,nchptot,1,nn)
call PASTIT (swet2,swetroot,istrip,nchp,nchptot,1,nn)
call PASTIT (swet3,swetdeep,istrip,nchp,nchptot,1,nn)
call PASTIT (capacity,capac,istrip,nchp,nchptot,1,nn)
call PASTIT (snowdepth,snodep,istrip,nchp,nchptot,1,nn)
C**********************************************************************
C Now update the theta and sh profiles with the new ground temperature
C**********************************************************************
do i =1,istrip
th(i,nlay+1) = tc(i) / pke(i,nlay+1)
enddo
do L = 1,nlay
do i =1,istrip
th(i,L) = th(i,L) + dthdthg(i,L)*(tc(i)-savetc(i))/pke(i,nlay+1)
enddo
enddo
do i =1,istrip
sh(i,nlay+1) = qa(i)
enddo
do L = 1,nlay
do i =1,istrip
sh(i,L) = sh(i,L) + dshdshg(i,L)*(qa(i)-saveqa(i))
enddo
enddo
do L = 1,nlay
do i =1,istrip
sttflux(i,L) = transth(i,L) * (th(i,L+1)-th(i,L))
stqflux(i,L) = transsh(i,L) * (sh(i,L+1)-sh(i,L))
enddo
enddo
C tendency updates
C ----------------
do l=1,nlay
call STRIP2TILE(uz(1,1,l),igrd,tmpstrip,nchp,ijall,
1 istrip,1,nn)
do i =1,istrip
tends(i) = ( u(i,l)-tmpstrip(i) )
enddo
call PASTIT (tends,dumsc(1,l),istrip,nchp,nchptot,1,nn)
call STRIP2TILE(vz(1,1,l),igrd,tmpstrip,nchp,ijall,
1 istrip,1,nn)
do i =1,istrip
tends(i) = ( v(i,l)-tmpstrip(i) )
enddo
call PASTIT (tends,dvmsc(1,l),istrip,nchp,nchptot,1,nn)
call STRIP2TILE(tz(1,1,l),igrd,tmpstrip,nchp,ijall,
1 istrip,1,nn)
do i =1,istrip
tends(i) = ( th(i,l)-tmpstrip(i) )
enddo
call PASTIT (tends,dtmsc(1,l),istrip,nchp,nchptot,1,nn)
call STRIP2TILE(qz(1,1,l,1),igrd,tmpstrip,nchp,ijall,
1 istrip,1,nn)
do i =1,istrip
tends(i) = ( sh(i,l)-tmpstrip(i) )
enddo
call PASTIT (tends,dqmsc(1,l,1),istrip,nchp,nchptot,1,nn)
c do nt = 1,ntracers-ptracers
c call strip2tile(qz(1,1,L,ptracers+nt),igrd,tmpstrip,nchp,
c 1 ijall,istrip,1,nn)
c do i =1,istrip
c tends(i) = ( tracers(i,L,nt)-tmpstrip(i) )
c enddo
c call pastit(tends,dqmsc(1,L,ptracers+nt),istrip,nchp,
c . nchptot,1,nn)
c enddo
enddo
C *********************************************************************
C **** increment diagnostic arrays for quantities saved in trbflx
C *********************************************************************
C note: the order, logic, and scaling of the heat and moisture flux
C diagnostics is critical!
C ------------------------------
if(diagnostics_is_on('EVAP ',myid) ) then
do i=1,istrip
tmpstrip(i) = stqflux(i,nlay) * 86400
enddo
call DIAG_VEGTILE_FILL (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
& .false., 'EVAP ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('EFLUX ',myid) ) then
do i=1,istrip
tmpstrip(i) = stqflux(i,nlay) * alhl
enddo
call DIAG_VEGTILE_FILL (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
& .false., 'EFLUX ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('HFLUX ',myid) ) then
do i=1,istrip
tmpstrip(i) = sttflux(i,nlay) * cp
enddo
call DIAG_VEGTILE_FILL (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
& .false., 'HFLUX ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('TUFLUX ',myid) ) then
call DIAG_VEGTILE_FILL (stuflux,igrd,chfrstr,istrip,nchp,nn,
& .false., 'TUFLUX ', 0, nlay, bi, bj, myid)
endif
if(diagnostics_is_on('TVFLUX ',myid) ) then
call DIAG_VEGTILE_FILL (stvflux,igrd,chfrstr,istrip,nchp,nn,
& .false., 'TVFLUX ', 0, nlay, bi, bj, myid)
endif
if(diagnostics_is_on('TTFLUX ',myid) ) then
do l=1,nlay
do i=1,istrip
sttflux(i,l) = sttflux(i,l) * cp
enddo
enddo
call DIAG_VEGTILE_FILL (sttflux,igrd,chfrstr,istrip,nchp,nn,
& .false., 'TTFLUX ', 0, nlay, bi, bj, myid)
endif
if(diagnostics_is_on('TQFLUX ',myid) ) then
do l=1,nlay
do i=1,istrip
stqflux(i,l) = stqflux(i,l) * alhl
enddo
enddo
call DIAG_VEGTILE_FILL (stqflux,igrd,chfrstr,istrip,nchp,nn,
& .false., 'TQFLUX ', 0, nlay, bi, bj, myid)
endif
if(diagnostics_is_on('RI ',myid) ) then
call DIAG_VEGTILE_FILL (sri,igrd,chfrstr,istrip,nchp,nn,
& .false., 'RI ', 0, nlay, bi, bj, myid)
endif
if(diagnostics_is_on('KH ',myid) ) then
call DIAG_VEGTILE_FILL (skh,igrd,chfrstr,istrip,nchp,nn,
& .false., 'KH ', 0, nlay, bi, bj, myid)
endif
if(diagnostics_is_on('KM ',myid) ) then
call DIAG_VEGTILE_FILL (skm,igrd,chfrstr,istrip,nchp,nn,
& .false., 'KM ', 0, nlay, bi, bj, myid)
endif
if(diagnostics_is_on('CT ',myid) ) then
call DIAG_VEGTILE_FILL (sct,igrd,chfrstr,istrip,nchp,nn,
& .false., 'CT ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('CU ',myid) ) then
call DIAG_VEGTILE_FILL (scu,igrd,chfrstr,istrip,nchp,nn,
& .false., 'CU ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('WINDS ',myid) ) then
call DIAG_VEGTILE_FILL (swinds,igrd,chfrstr,istrip,nchp,nn,
& .false., 'WINDS ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('UFLUX ',myid) ) then
call DIAG_VEGTILE_FILL (stuflux(1,nlay),igrd,chfrstr,istrip,nchp,
& nn,.false., 'UFLUX ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('VFLUX ',myid) ) then
call DIAG_VEGTILE_FILL (stvflux(1,nlay),igrd,chfrstr,istrip,nchp,
& nn,.false., 'VFLUX ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('USTAR ',myid) ) then
call DIAG_VEGTILE_FILL (sustar,igrd,chfrstr,istrip,nchp,nn,
& .false., 'USTAR ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('Z0 ',myid) ) then
call DIAG_VEGTILE_FILL (sz0,igrd,chfrstr,istrip,nchp,nn,
& .false., 'Z0 ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('FRQTRB ',myid) ) then
call DIAG_VEGTILE_FILL (frqtrb,igrd,chfrstr,istrip,nchp,nn,
& .false., 'FRQTRB ', 0, nlay-1, bi, bj, myid)
endif
if(diagnostics_is_on('PBL ',myid) ) then
call DIAG_VEGTILE_FILL (pbldpth,igrd,chfrstr,istrip,nchp,nn,
& .false., 'PBL ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('U2M ',myid) ) then
call DIAG_VEGTILE_FILL (stu2m,igrd,chfrstr,istrip,nchp,nn,
& .false., 'U2M ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('V2M ',myid) ) then
call DIAG_VEGTILE_FILL (stv2m,igrd,chfrstr,istrip,nchp,nn,
& .false., 'V2M ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('T2M ',myid) ) then
call DIAG_VEGTILE_FILL (stt2m,igrd,chfrstr,istrip,nchp,nn,
& .false., 'T2M ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('Q2M ',myid) ) then
do i=1,istrip
if( stq2m(i).ne.undef ) then
tmpstrip(i) = stq2m(i) * 1000
else
tmpstrip(i) = undef
endif
enddo
call DIAG_VEGTILE_FILL (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
& .false., 'Q2M ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('U10M ',myid) ) then
call DIAG_VEGTILE_FILL (stu10m,igrd,chfrstr,istrip,nchp,nn,
& .false., 'U10M ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('V10M ',myid) ) then
call DIAG_VEGTILE_FILL (stv10m,igrd,chfrstr,istrip,nchp,nn,
& .false., 'V10M ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('T10M ',myid) ) then
call DIAG_VEGTILE_FILL (stt10m,igrd,chfrstr,istrip,nchp,nn,
& .false., 'T10M ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('Q10M ',myid) ) then
do i=1,istrip
if( stq10m(i).ne.undef ) then
tmpstrip(i) = stq10m(i) * 1000
else
tmpstrip(i) = undef
endif
enddo
call DIAG_VEGTILE_FILL (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
& .false., 'Q10M ', 1, 1, bi, bj, myid)
endif
C**********************************************************************
C more diagnostics: land surface model parameters
C**********************************************************************
if(diagnostics_is_on('TDEEP ',myid) ) then
call DIAG_VEGTILE_FILL (td,igrd,chfrstr,istrip,nchp,nn,
& .false., 'TDEEP ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('QCANOPY ',myid) ) then
call DIAG_VEGTILE_FILL (qa,igrd,chfrstr,istrip,nchp,nn,
& .false., 'QCANOPY ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('SMSHAL ',myid) ) then
call DIAG_VEGTILE_FILL (swet1,igrd,chfrstr,istrip,nchp,nn,
& .false., 'SMSHAL ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('SMROOT ',myid) ) then
call DIAG_VEGTILE_FILL (swet2,igrd,chfrstr,istrip,nchp,nn,
& .false., 'SMROOT ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('SMDEEP ',myid) ) then
call DIAG_VEGTILE_FILL (swet3,igrd,chfrstr,istrip,nchp,nn,
& .false., 'SMDEEP ', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('CAPACITY',myid) ) then
call DIAG_VEGTILE_FILL (capacity,igrd,chfrstr,istrip,nchp,nn,
& .false., 'CAPACITY', 1, 1, bi, bj, myid)
endif
if(diagnostics_is_on('SNOW ',myid) ) then
call DIAG_VEGTILE_FILL (snowdepth,igrd,chfrstr,istrip,nchp,nn,
& .false., 'SNOW ', 1, 1, bi, bj, myid)
endif
C**********************************************************************
C end regions loop
2000 continue
C**********************************************************************
C increment the counter for the accumulated fcc and qliq arrays
C ---------------------------------------------------------------------
imstturblw = imstturblw + 1
imstturbsw = imstturbsw + 1
C prevent ice or snow from melting
C ---------------------------------------------------------------------
do i =1,nchptot
if( (icethk(i).gt.0.).and.(tground(i).gt.tice) ) tground(i)=tice
enddo
C Update tcanopy and ecanopy from the points of the
C tground and qground arrays
C ---------------------------------------------------------------------
do i =1,nchptot
tcanopy(i) = tground(i)
ecanopy(i) = qground(i)
enddo
C Initialize Tendencies and Couplings
C -----------------------------------
do L = 1,nlay
do i = 1,ijall
duturb(i,1,L) = 0.
dvturb(i,1,L) = 0.
dtturb(i,1,L) = 0.
qqgrid(i,1,L) = 0.
qliqtmp(i,1,L) = 0.
fcctmp(i,1,L) = 0.
enddo
do nt = 1,ntracers
do i = 1,ijall
dqturb(i,1,L,nt) = 0.
enddo
enddo
enddo
C Return Tendencies and Couplings to Grid Space
C ---------------------------------------------
do l = 1,nlay
call MSC2GRD(igrd,chfr,dumsc(1,L),nchp,nchptot,fracland,
& duturb(1,1,L),im,jm)
call MSC2GRD(igrd,chfr,dvmsc(1,L),nchp,nchptot,fracland,
& dvturb(1,1,L),im,jm)
call MSC2GRD(igrd,chfr,dtmsc(1,L),nchp,nchptot,fracland,
& dtturb(1,1,L),im,jm)
do nt = 1,ntracers
call MSC2GRD(igrd,chfr,dqmsc(1,L,nt),nchp,nchptot,fracland,
& dqturb(1,1,L,nt),im,jm)
enddo
call MSC2GRD(igrd,chfr, tke(1,L),nchp,nchptot,fracland,
& qqgrid(1,1,L),im,jm)
call MSC2GRD(igrd,chfr, fccmsc(1,L),nchp,nchptot,fracland,
& fcctmp(1,1,L),im,jm)
call MSC2GRD(igrd,chfr,qliqmsc(1,L),nchp,nchptot,fracland,
& qliqtmp(1,1,L),im,jm)
enddo
C Reduce clouds from conditionally unstable layer
C -----------------------------------------------
call CTEI ( tz,qz,fcctmp,qliqtmp,plz,pkz,pkht,im*jm,nlay )
C Bump Total Cloud Liquid Water and Fraction by Instantaneous Values
C ------------------------------------------------------------------
do l = 1,nlay
do i = 1,ijall
fccavesw(i,1,L) = fccavesw(i,1,L) + fcctmp(i,1,L)
fccavelw(i,1,L) = fccavelw(i,1,L) + fcctmp(i,1,L)
qliqavelw(i,1,L) = qliqavelw(i,1,L) + qliqtmp(i,1,L)
qliqavesw(i,1,L) = qliqavesw(i,1,L) + qliqtmp(i,1,L)
enddo
enddo
tmpFac = 1.e6
CALL DIAGNOSTICS_FILL(fcctmp,'TRBFCC ',0,nlay,3,bi,bj,myid)
CALL DIAGNOSTICS_SCALE_FILL( qliqtmp,tmpFac,1,
& 'TRBQLIQ ',0,nlay,3,bi,bj,myid)
C**********************************************************************
C And some other variables to be transformed back to grid space:
C Ground Temperature, snow depth and shallow layer ground wetness
do i = 1,ijall
tgz(i,1) = 0.
enddo
call MSC2GRD(igrd,chfr,tground ,nchp,nchptot,fracland,tgz ,im,jm)
C *********************************************************************
C **** increment diagnostic array for ground and surface temperatures,
C *** ground temp tendency, and ground humidity
C *********************************************************************
call DIAGNOSTICS_FILL(tgz,'TGROUND ',0,1,3,bi,bj,myid)
call DIAGNOSTICS_FILL(tgz,'TCANOPY ',0,1,3,bi,bj,myid)
if(diagnostics_is_on('TS ',myid) ) then
do i = 1,ijall
tmpdiag(i,1) = tz(i,1,nlay) * pkht(i,1,nlay)
enddo
call DIAGNOSTICS_FILL(tmpdiag,'TS ',0,1,3,bi,bj,myid)
endif
call DIAGNOSTICS_FILL(tgz,'DTG ',0,1,3,bi,bj,myid)
C *********************************************************************
C **** increment diagnostic arrays for tendencies ****
C *********************************************************************
tmpFac = atimstp * secday
CALL DIAGNOSTICS_SCALE_FILL(dvturb,tmpFac,1,
& 'TURBV ',0,nlay,3,bi,bj,myid)
CALL DIAGNOSTICS_SCALE_FILL(duturb,tmpFac,1,
& 'TURBU ',0,nlay,3,bi,bj,myid)
tmpFac = atimstp * secday * 1000.
CALL DIAGNOSTICS_SCALE_FILL(dqturb,tmpFac,1,
& 'TURBQ ',0,nlay,3,bi,bj,myid)
if(diagnostics_is_on('TURBT ',myid) ) then
do L = 1,nlay
do i = 1,ijall
tmpdiag(i,1) = dtturb(i,1,l) * pkz(i,1,l)*atimstp*secday
enddo
call DIAGNOSTICS_FILL(tmpdiag,'TURBT ',L,1,3,bi,bj,myid)
enddo
endif
C pi-weight the theta and moisture tendencies
C -------------------------------------------
do i =1,ijall
thtgz(i) = pz(i,1) * atimstp
enddo
do l =1,nlay
do i =1,ijall
duturb(i,1,l) = duturb(i,1,l)*atimstp
dvturb(i,1,l) = dvturb(i,1,l)*atimstp
dtturb(i,1,l) = dtturb(i,1,l)*thtgz(i)
enddo
do nt = 1,ntracers
do i =1,ijall
dqturb(i,1,l,nt) = dqturb(i,1,l,nt)*thtgz(i)
enddo
enddo
enddo
C *********************************************************************
C **** zero out the accumulating rainfall and snowfall arrays ***
C *********************************************************************
if( time_left.lt.timstp ) then
do i =1,ijall
rainlsp(i,1) = 0.
rainconv(i,1) = 0.
snowfall(i,1) = 0.
enddo
endif
return
end
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
SUBROUTINE TRBFLX (NN,TH,THV,SH,U,V,QQ,PL,PLE,PLK,PLKE,DPSTR,
1 IWATER,Z0,tracers,ntrace,ntracedim,DTAU,ITRTRB,KMBG,KHBG,QBEG,
2 TPROF,WU,WV,SRI,ET,EU,SWINDS,sustar,sz0,freqdg,pbldpth,
3 sct,scu,stu2m,stv2m,stt2m,stq2m,stu10m,stv10m,stt10m,stq10m,
4 irun,nlev,nltop,NYMD,NHMS,grav,cp,rgas,faceps,virtcon,undef,
5 dshdthg,dshdshg,dthdthg,dthdshg,eturb,dedqa,dedtc,
6 hsturb,dhsdqa,dhsdtc,transth,transsh,
7 ctsave,xxsave,yysave,zetasave,xlsave,khsave,qliq,turbfcc)
C**********************************************************************
C SUBROUTINE TRBFLX - COMPUTES TURBULENT ADJUSTMENTS TO ATMOSPHERIC
C PROFILE
C - CALLED FROM PBL DRIVER
C ARGUMENTS ::
C INPUT:
C ------
C TH - POTENTIAL TEMPERATURE PROFILE
C THV - VIRTUAL POTENTIAL TEMPERATURE PROFILE
C SH - SPECIFIC HUMIDITY PROFILE
C U - U - COMPONENT OF WIND PROFILE
C V - V - COMPONENT OF WIND PROFILE
C QQ - TURBULENT KINETIC ENERGY
C PL - EVEN LEVEL PRESSURES
C PLE - EDGE LEVEL PRESSURES
C PLK - EVEN LEVEL PRESSURES ** KAPPA
C PLKE - EDGE LEVEL PRESSURES ** KAPPA
C DPSTR - PRESSURE INTERVALS
C WATER - BIT ARRAY - '1' OVER OCEANS
C Z0 - SURFACE ROUGHNESS
C tracers - array of passive tracers
C ntrace - number of tracers to be diffused
C ntracedim - outer dimension of tracers array
C DTAU - TIME CHANGE PER ITERATION OF TRBFLX
C ITRTRB - NUMBER OF ITERATIONS OF TRBFLX
C KMBG - BACKGROUND VALUE OF MOMENTUM TRANSFER COEF
C KHBG - BACKGROUND VALUE OF HEAT TRANSFER COEF
C NLEV - NUMBER OF ATMOSPHERIC LEVELS TO CALCULATE
C nltop - Top level allowed for turbulence
C QBEG - LOGICAL .TRUE. FOR INITIAL START OF GCM
C TPROF - LOGICAL .TRUE. TO CALCULATE PT BY PT DIAGS
C
C OUTPUT:
C -------
C PROFILES RETURNED WITH UPDATED VALUES
C**********************************************************************
implicit none
C Argument list declarations
integer nn,irun,nlev,nltop,ntrace,ntracedim,itrtrb,nhms,nymd
_RL TH(irun,NLEV+1),THV(irun,NLEV+1),SH(irun,NLEV+1)
_RL U(irun,NLEV+1),V(irun,NLEV+1),QQ(irun,NLEV)
_RL PL(irun,NLEV),PLE(irun,NLEV+1),PLK(irun,NLEV)
_RL PLKE(irun,NLEV+1),DPSTR(irun,NLEV)
integer IWATER(irun)
_RL Z0(irun)
_RL tracers(irun,nlev+1,ntracedim)
_RL dtau,KMBG,KHBG
LOGICAL QBEG,TPROF
_RL SWINDS(irun)
_RL SRI(irun,nlev), ET(irun,nlev)
_RL EU (irun,nlev)
_RL WU(irun,nlev)
_RL WV (irun,nlev), pbldpth(irun)
_RL sustar(irun), sz0(irun)
_RL freqdg(irun,nlev-1)
_RL sct(irun), scu(irun)
_RL stu2m(irun),stv2m(irun),stt2m(irun),stq2m(irun)
_RL stu10m(irun),stv10m(irun),stt10m(irun),stq10m(irun)
_RL grav,cp,rgas,faceps,virtcon,undef
_RL eturb(irun),dedqa(irun),dedtc(irun)
_RL hsturb(irun),dhsdqa(irun),dhsdtc(irun)
_RL dshdthg(irun,nlev),dthdthg(irun,nlev)
_RL dshdshg(irun,nlev),dthdshg(irun,nlev)
_RL transth(irun,nlev),transsh(irun,nlev)
_RL ctsave(irun),xxsave(irun),yysave(irun)
_RL zetasave(irun),xlsave(irun,nlev),khsave(irun,nlev)
_RL qliq(irun,nlev),turbfcc(irun,nlev)
C Local Variables
_RL b1,b3,alpha,halpha,qqmin,qbustr
PARAMETER ( B1 = 16.6 )
PARAMETER ( B3 = 1. / B1 )
PARAMETER ( ALPHA = 0.1 )
PARAMETER ( HALPHA = ALPHA * 0.5 )
PARAMETER ( QQMIN = 0.005 )
PARAMETER ( QBUSTR = 2.550952 )
_RL argmax, onethrd, z1pem25, b2, two
PARAMETER (ARGMAX = 30.)
PARAMETER (ONETHRD = 1./3. )
PARAMETER (Z1PEM25 = 1.E-25)
PARAMETER ( B2 = 10.1 )
PARAMETER ( two = 2.0 )
_RL AHS (irun), HS(irun)
_RL XX (irun), YY(irun), CU(irun)
_RL CT(irun), USTAR(irun)
_RL RIB(irun), ZETA(irun), WS(irun)
_RL DTHS(irun), DELTHS(irun)
_RL DTHL(irun), DELTHL(irun)
_RL RIBIN(irun),CUIN(irun)
_RL CTIN(irun),ZETAIN(irun)
_RL USTARIN(irun),RHOSIN(irun),Z0IN(irun)
_RL qqcolmin(irun),qqcolmax(irun),levpbl(irun)
_RL ADZ1(irun,nlev), DZ1TMP(irun,nlev)
_RL DZ3(irun,nlev), TEMP(irun,nlev)
_RL DV(irun,nlev), DTHV(irun,nlev)
_RL DPK(irun,nlev), STRT(irun,nlev)
_RL DW2(irun,nlev), RI(irun,nlev)
_RL RHOZPK(irun,nlev), Q(irun,nlev)
_RL RIINIT(irun,nlev), DU(irun,nlev)
_RL QQINIT(irun,nlev), RHOKDZ(irun,nlev)
_RL RHODZ2(irun,nlev)
_RL KM(irun,nlev), KH(irun,nlev)
_RL DELTH (irun,nlev+1), DELSH (irun,nlev+1)
_RL FLXFAC (irun,nlev+1)
_RL FLXFPK (irun,nlev+1)
_RL ADZ2 (irun,nlev-1), RHODZ1(irun,nlev-1)
_RL VKZE (irun,nlev-1), VKZM (irun,nlev-1)
_RL XL (irun,nlev-1), QXLM (irun,nlev-1)
_RL QQE (irun,nlev-1), QE (irun,nlev-1)
_RL P3 (irun,nlev-1), XQ (irun,nlev-1)
_RL XLDIAG (irun,nlev-1), FLXFCE(irun,nlev-1)
LOGICAL FIRST,LAST
integer IBITSTB(irun,nlev),INTQ(irun,nlev)
C arrays for use by moist bouyancy calculation
C -----------------
_RL TL(irun,NLEV),DTH(irun,NLEV)
_RL DSH(irun,NLEV)
_RL SHL(irun,NLEV)
_RL AA(irun,NLEV),BB(irun,NLEV),SSDEV(irun,NLEV)
_RL ARG(irun,NLEV),XXZETA(irun),QBYU(irun)
_RL SVAR(irun,NLEV),Q1M(irun,NLEV)
_RL FCC(irun,NLEV)
_RL BETAT(irun,NLEV),BETAW(irun,NLEV)
_RL BETAL(irun,NLEV),BETAT1(irun,NLEV)
_RL BETAW1(irun,NLEV),SBAR(irun,NLEV)
_RL SHSAT(irun,NLEV)
C Some space for variables to be used in called routines
logical LWATER
integer IVBITRIB(irun)
_RL VHZ(irun)
_RL VH0(irun)
_RL VPSIM(irun),VAPSIM(irun)
_RL VPSIG(irun),VPSIHG(irun)
_RL VTEMP(irun),VDZETA(irun)
_RL VDZ0(irun),VDPSIM(irun)
_RL VDPSIH(irun),VZH(irun)
_RL VXX0(irun),VYY0(irun)
_RL VAPSIHG(irun),VRIB1(irun),VWS1(irun)
_RL VPSIH(irun),VZETAL(irun)
_RL VZ0L(irun),VPSIH2(irun)
_RL VX0PSIM(irun),VG(irun),VG0(irun),VR1MG0(irun)
_RL VZ2(irun),VDZSEA(irun),VAZ0(irun),VXNUM1(irun)
_RL VPSIGB2(irun),VDX(irun),VDXPSIM(irun),VDY(irun)
_RL VXNUM2(irun),VDEN(irun),VAWS1(irun),VXNUM3(irun)
_RL VXNUM(irun),VDZETA1(irun),VDZETA2(irun)
_RL VZCOEF2(irun),VZCOEF1(irun),VTEMPLIN(irun)
_RL VDPSIMC(irun),VDPSIHC(irun)
_RL DZITRP(irun,nlev-1),STBFCN(irun,nlev)
_RL XL0(irun,nlev),Q1(irun,nlev-1)
_RL WRKIT1(irun,nlev-1)
_RL WRKIT2(irun,nlev-1)
_RL WRKIT3(irun,nlev-1)
_RL WRKIT4(irun,nlev-1)
INTEGER INT1(irun,nlev), INT2(irun,nlev-1)
_RL vrt1con,pi,rsq2pi,p5sr,clh,vk,rvk,aitr,gbycp,fac1,fac2
_RL getcon,dum,errf
integer istnlv,nlevm1,nlevm2,nlevml,nlevp1,istnm1,istnm2,istnp1
integer istnml,istnmq,istlmq,nlevmq
integer i,iter,init,n,LL,L,Lp,Lp1,lmin,lminq,lminq1,ibit
vk = getcon('VON KARMAN')
rvk = 1./vk
AITR = 1. / FLOAT(ITRTRB)
ISTNLV = irun * NLEV
NLEVM1 = NLEV - 1
NLEVM2 = NLEV - 2
NLEVP1 = NLEV + 1
ISTNM1 = irun * NLEVM1
ISTNM2 = irun * NLEVM2
ISTNP1 = irun * NLEVP1
GBYCP = GRAV / CP
VRT1CON = 1. + VIRTCON
PI = 4. * ATAN(1.)
RSQ2PI = 1./ ((2.*PI)**0.5)
P5SR = 0.5**0.5
CLH = GETCON('LATENT HEAT COND') / CP
C SET INITIAL NUMBER OF ITERATIONS OF SFCFLX
C ------------------------------------------
N = 6
C DETERMINE IF INITIAL START
C --------------------------
INIT = 0
IF(QBEG) INIT = 1
C SET DIAGNOSTIC LOGICALS AND INITIALIZE DIAGNOSTIC ARRAYS
C --------------------------------------------------------
DO L =1,NLEV
DO I =1,irun
wu(I,L) = 0.
wv(I,L) = 0.
eu(I,L) = 0.
et(I,L) = 0.
ENDDO
ENDDO
if (tprof) then
DO L =1,NLEVM1
DO I =1,irun
XLDIAG(I,L) = 0.
ENDDO
ENDDO
endif
DO L =1,NLEVM1
DO I =1,irun
FREQDG(I,L) = 0.
ENDDO
ENDDO
do I =1,irun
scu(i) = 0.
enddo
do I =1,irun
sct(i) = 0.
enddo
do I =1,irun
pbldpth(i) = 0.
enddo
do I =1,irun
sustar(i) = 0.
enddo
do I =1,irun
sz0(i) = 0.
enddo
c do I =1,ISTNM1
c FREQDG(I,1) = 0.
c enddo
do I =1,irun
stu2m(i) = 0.
enddo
do I =1,irun
stv2m(i) = 0.
enddo
do I =1,irun
stt2m(i) = 0.
enddo
do I =1,irun
stq2m(i) = 0.
enddo
do I =1,irun
stu10m(i) = 0.
enddo
do I =1,irun
stv10m(i) = 0.
enddo
do I =1,irun
stt10m(i) = 0.
enddo
do I =1,irun
stq10m(i) = 0.
enddo
IF (INIT.EQ.1) THEN
DO L =1,NLEVM1
DO I =1,irun
XLSAVE(I,L) = 0.
KHSAVE(I,L) = 0.
ENDDO
ENDDO
DO I = 1,irun
CTSAVE(I) = 0.
XXSAVE(I) = 0.
YYSAVE(I) = 0.
ZETASAVE(I) = 0.
ENDDO
ENDIF
C COMPUTE VERTICAL GRID
C ---------------------
DO L =1,NLEV
DO I =1,irun
ADZ1(I,L) = (CP/GRAV)*(PLKE(I,L+1)-PLKE(I,L))
ADZ1(I,L) = THV(I,L) * ADZ1(I,L)
DZ1TMP(I,L) = ADZ1(I,L)
ENDDO
ENDDO
DO L =1,NLEVM1
DO I =1,irun
ADZ2(I,L) = 0.5 * (ADZ1(I,L)+ADZ1(I,L+1))
ENDDO
ENDDO
C DEPTH HS OF SURFACE LAYER
C -------------------------
DO 9042 I =1,irun
HS(I) = 0.5 * ADZ1(I,NLEV)
9042 CONTINUE
C ALPHA * LAYER DEPTHS FOR TRBLEN
C -------------------------------
DO 9044 I =1,irun
DZ3(I,1) = HALPHA * ADZ1(I,1)
9044 CONTINUE
DO L =2,NLEVM1
DO I =1,irun
DZ3(I,L) = ALPHA * ADZ1(I,L)
ENDDO
ENDDO
DO 9048 I =1,irun
DZ3(I,NLEV) = ALPHA * HS(I)
9048 CONTINUE
C VK * HEIGHTS AT MID AND EDGE LEVELS
C -----------------------------------
DO L =2,NLEV
DO I =1,irun
TEMP(I,L) = VK * ADZ1(I,L)
ENDDO
ENDDO
DO 9052 I =1,irun
VKZE(I,NLEVM1) = TEMP(I,NLEV)
9052 CONTINUE
DO 100 LL = 2,NLEVM1
L = NLEV - LL
LP1 = L + 1
DO 9054 I =1,irun
VKZE(I,L) = VKZE(I,LP1) + TEMP(I,LP1)
9054 CONTINUE
100 CONTINUE
DO L =1,NLEVM1
DO I =1,irun
VKZM(I,L) = VKZE(I,L) - 0.5 * TEMP(I,L+1)
ENDDO
ENDDO
C COMPUTE RHO BY DZ AT MID AND EDGE LEVELS
C ----------------------------------------
DO 200 L = 1,NLEVM1
LP1 = L + 1
DO 9058 I =1,irun
FAC1 = DPSTR(I,L) / ( DPSTR(I,L) + DPSTR(I,LP1) )
FAC2 = 1. - FAC1
RHODZ2(I,L) = FAC1 * THV(I,LP1)
RHODZ2(I,L) = RHODZ2(I,L) + FAC2 * THV(I,L)
9058 CONTINUE
200 CONTINUE
DO L =1,NLEVM1
DO I =1,irun
RHODZ2(I,L) = (RGAS*0.01) * RHODZ2(I,L)
TEMP(I,L) = PLKE(I,L+1) * ADZ2(I,L)
RHODZ2(I,L) = TEMP(I,L) * RHODZ2(I,L)
RHODZ2(I,L) = PLE(I,L+1) / RHODZ2(I,L)
RHOZPK(I,L) = RHODZ2(I,L) * PLKE(I,L+1)
RHODZ1(I,L) = (RGAS*0.01) * THV(I,L+1)
TEMP(I,L) = PLK(I,L+1) * ADZ1(I,L+1)
RHODZ1(I,L) = TEMP(I,L) * RHODZ1(I,L)
RHODZ1(I,L) = PL(I,L+1) / RHODZ1(I,L)
ENDDO
ENDDO
C COMPUTE FLXFAC FOR LAYERS AND EDGES
C COMPUTE DTG / DT DUE TO RADIATION AND HEAT CONDUCTION THROUGH ICE
C -----------------------------------------------------------------
DO L =1,NLEV
DO I =1,irun
FLXFPK(I,L) = PLE(I,L+1) - PLE(I,L)
FLXFPK(I,L) = FLXFPK(I,L) * PLK(I,L)
FLXFPK(I,L) = (GRAV*DTAU*0.01) / FLXFPK(I,L)
ENDDO
ENDDO
DO 9064 I =1,irun
FLXFPK(I,NLEVP1) = 0.
9064 CONTINUE
DO 9066 I =1,irun
IF (IWATER(I).EQ.0 ) FLXFPK(I,NLEVP1) = 1. / PLKE(I,NLEVP1)
9066 CONTINUE
DO L =1,NLEV
DO I =1,irun
FLXFAC(I,L) = FLXFPK(I,L) * PLK(I,L)
ENDDO
ENDDO
DO 9070 I =1,irun
FLXFAC(I,NLEVP1) = FLXFPK(I,NLEVP1)
9070 CONTINUE
DO 9074 I =1,irun
FLXFPK(I,NLEVP1) = CP * FLXFPK(I,NLEVP1)
9074 CONTINUE
DO L =1,NLEVM1
DO I =1,irun
FLXFCE(I,L) = PL(I,L+1) - PL(I,L)
FLXFCE(I,L) = (GRAV*DTAU*0.01) / FLXFCE(I,L)
ENDDO
ENDDO
C COMPUTE RECIPROCALS OF DZ1, DZ2, HS
C -----------------------------------
DO L =1,NLEV
DO I =1,irun
ADZ1(I,L) = 1. / ADZ1(I,L)
ENDDO
ENDDO
DO L =1,NLEVM1
DO I =1,irun
ADZ2(I,L) = 1. / ADZ2(I,L)
ENDDO
ENDDO
DO 9088 I =1,irun
AHS(I) = 1. / HS(I)
9088 CONTINUE
C COMPUTE GRADIENTS OF P**KAPPA
C -----------------------------
DO L =1,NLEVM1
DO I =1,irun
DPK(I,L) = ( PLK(I,L+1)-PLK(I,L) ) * ADZ2(I,L)
ENDDO
ENDDO
DO 9092 I =1,irun
DPK(I,NLEV) = GBYCP / THV(I,NLEV)
9092 CONTINUE
C INITIALIZE Q ARRAY
C ------------------
DO L =1,NLEVM1
DO I =1,irun
Q(I,L) = 2. * QQ(I,L)
Q(I,L) = SQRT( Q(I,L) )
ENDDO
ENDDO
FIRST = .TRUE.
LAST = .FALSE.
C**********************************************************************
C**********************************************************************
C MAIN LOOP
DO 2000 ITER = 1, ITRTRB
IF ( ITER .GE. ITRTRB ) LAST = .TRUE.
C CODE FOR MOIST BOUNDARY LAYER - NEW CALCULATION OF DTHV
IF(ITER.EQ.1) THEN
DO I = 1,irun
CT(I) = CTSAVE(I)
XX(I) = XXSAVE(I)
YY(I) = YYSAVE(I)
ZETA(I) = ZETASAVE(I)
ENDDO
ENDIF
DO I = 1,irun
TL(I,NLEV) = TH(I,NLEV)*PLK(I,NLEV)
call QSAT ( tl(i,nlev),pl(i,nlev),shsat(i,nlev),dum,.false. )
ENDDO
DO I = 1,irun
BB(I,NLEV) = FACEPS*SHSAT(I,NLEV)/(TL(I,NLEV)*TL(I,NLEV))
AA(I,NLEV) = 1. / (1. + CLH * BB(I,NLEV) )
BB(I,NLEV) = BB(I,NLEV) * AA(I,NLEV) * plk(I,nlev)
DTH(I,NLEV) = TH(I,NLEV)-TH(I,NLEVP1)
DSH(I,NLEV) = SH(I,NLEV)-SH(I,NLEVP1)
SBAR(I,NLEV) = AA(I,NLEV) * (SH(I,NLEV) - SHSAT(I,NLEV))
SSDEV(I,NLEV)=CT(I)*(AA(I,NLEV)*DSH(I,NLEV)
1 -BB(I,NLEV)*DTH(I,NLEV))
XXZETA(I) = XX(I)-ZETA(I)
IF(XXZETA(I).LT.0.1*XX(I)) XXZETA(I)=0.1*XX(I)
IF(XXZETA(I).LE.0.) XXZETA(I)=0.1
QBYU(I) =QBUSTR * XXZETA(I) ** ONETHRD
SSDEV(I,NLEV) = B2*YY(I)*SSDEV(I,NLEV)*SSDEV(I,NLEV)/QBYU(I)
SVAR(I,NLEV) = SQRT(SSDEV(I,NLEV))
IF ( SVAR(I,NLEV).LT.Z1PEM25) SVAR(I,NLEV) = Z1PEM25
Q1M(I,NLEV) = SBAR(I,NLEV) / SVAR(I,NLEV)
FCC(I,NLEV) = (1./2.) * ( 1. + ERRF( P5SR*Q1M(I,NLEV) ) )
SHL(I,NLEV) = FCC(I,NLEV) * SBAR(I,NLEV)
ARG(I,NLEV) = (1./2.)*Q1M(I,NLEV)*Q1M(I,NLEV)
IF(ARG(I,NLEV).LE.ARGMAX)
1 SHL(I,NLEV) = SHL(I,NLEV)+RSQ2PI*SVAR(I,NLEV)*EXP(-ARG(I,NLEV))
BETAT(I,NLEV) = 1. + VIRTCON*SH(I,NLEV) - VRT1CON*SHL(I,NLEV)
BETAW(I,NLEV) = VIRTCON *
1 ( TH(I,NLEV) + CLH * SHL(I,NLEV) * (1./plk(i,nlev)) )
BETAL(I,NLEV) = (1.+VIRTCON*SH(I,NLEV)-TWO*VRT1CON*SHL(I,NLEV))
1 * (1./plk(i,nlev)) * CLH - VRT1CON * TH(I,NLEV)
BETAT1(I,NLEV) = BETAT(I,NLEV) - BB(I,NLEV)*FCC(I,NLEV)
1 * BETAL(I,NLEV)
BETAW1(I,NLEV) = BETAW(I,NLEV) + AA(I,NLEV) * FCC(I,NLEV)
1 * BETAL(I,NLEV)
DTHV(I,NLEV) = BETAT1(I,NLEV)*DTH(I,NLEV) +
1 BETAW1(I,NLEV)*DSH(I,NLEV)
THV(I,NLEVP1) = THV(I,NLEV) - DTHV(I,NLEV)
ENDDO
C SURFACE FLUX TRANSFER COEFFICIENTS
CALL SFCFLX(NN,U(1,NLEV),V(1,NLEV),
1 THV(1,NLEV),
2 THV(1,NLEVP1),TH(1,NLEV),TH(1,NLEVP1),
3 SH(1,NLEV),SH(1,NLEVP1),PLK(1,NLEV),
4 PLKE(1,NLEVP1),PLE(1,NLEVP1),Z0,
5 IWATER,HS,AHS,
6 FIRST,LAST,N,irun,aitr,RHODZ2(1,NLEV),RHOZPK(1,NLEV),
7 KH(1,NLEV),KM(1,NLEV),USTAR,
8 XX,YY,CU,
9 CT,RIB,ZETA,WS,
1 stu2m,stv2m,stt2m,stq2m,stu10m,stv10m,stt10m,stq10m,
2 cp,rgas,undef,
3 lwater, ivbitrib,
4 VHZ,VPSIM,VAPSIM,VPSIG,VPSIHG,VTEMP,VDZETA,VDZ0,VDPSIM,
5 VDPSIH,VZH,VXX0,VYY0,VAPSIHG,VRIB1,VWS1,VPSIH,
9 VZETAL,VZ0L,VPSIH2,VH0,
1 VX0PSIM,VG,VG0,VR1MG0,VZ2,VDZSEA,VAZ0,VXNUM1,VPSIGB2,VDX,
2 VDXPSIM,VDY,VXNUM2,VDEN,VAWS1,VXNUM3,VXNUM,VDZETA1,VDZETA2,
3 VZCOEF2,VZCOEF1,VTEMPLIN,VDPSIMC,VDPSIHC)
CI
N = 1
C SET VALUES OF TURBULENT VELOCITY AND KINETIC ENERGY AT THE GROUND
CB
DO 9098 I =1,irun
Q(I,NLEV) = QBUSTR * USTAR(I)
QQ(I,NLEV) = 0.5 * Q(I,NLEV) * Q(I,NLEV)
9098 CONTINUE
CE
C GRADIENTS
C ---------
DO L =1,NLEVM1
DO I =1,irun
DU(I,L) = ( U(I,L)- U(I,L+1) ) * ADZ2(I,L)
DV(I,L) = ( V(I,L)- V(I,L+1) ) * ADZ2(I,L)
ENDDO
ENDDO
C NEW CODE FOR MOIST BOUNDARY LAYER - NEW CALCULATION OF DTHV
IF(ITER.EQ.1) THEN
DO L =1,NLEVM1
DO I =1,irun
XL(I,L) = XLSAVE(I,L)
ENDDO
ENDDO
ENDIF
DO L =1,NLEVM1
DO I =1,irun
DTH(I,L) = ( TH(I,L)-TH(I,L+1) ) * ADZ2(I,L)
DSH(I,L) = ( SH(I,L)-SH(I,L+1) ) * ADZ2(I,L)
TL(I,L) = TH(I,L)*PLK(I,L)
ENDDO
ENDDO
DO LL = 1,NLEVM1
DO I = 1,irun
call QSAT ( tl(i,LL),pl(i,LL),shsat(i,LL),dum,.false. )
ENDDO
ENDDO
DO L =1,NLEVM1
DO I =1,irun
BB(I,L) = FACEPS*SHSAT(I,L)/(TL(I,L)*TL(I,L))
AA(I,L) = 1. / (1. + CLH * BB(I,L) )
COMMM BB(I,L) = BB(I,L) * AA(I,L) * plke(I,L+1)
BB(I,L) = BB(I,L) * AA(I,L)
SBAR(I,L) = AA(I,L) * (SH(I,L) - SHSAT(I,L))
ENDDO
ENDDO
DO I = 1,irun
COMMM SSDEV(I,1) = XL(I,1)*(AA(I,1)*DSH(I,1)-BB(I,1)*DTH(I,1))
SSDEV(I,1) = XL(I,1)*(AA(I,1)*DSH(I,1)-
1 BB(I,1)*plke(I,2)*DTH(I,1))
SSDEV(I,1) = B2 * KHSAVE(I,1) * SSDEV(I,1) * SSDEV(I,1)
SVAR(I,1) = SQRT(SSDEV(I,1))
IF ( SVAR(I,1).LT.Z1PEM25) SVAR(I,1) = Z1PEM25
ENDDO
DO L =2,NLEVM1
DO I =1,irun
COMMM SSDEV(I,L) = XL(I,L-1)*(AA(I,L)*DSH(I,L-1)-BB(I,L)*DTH(I,L-1))
SSDEV(I,L) = XL(I,L-1)*(AA(I,L)*DSH(I,L-1)-
1 BB(I,L)*plke(I,L)*DTH(I,L-1))
SSDEV(I,L) = B2 * KHSAVE(I,L-1) * SSDEV(I,L-1) * SSDEV(I,L-1)
SVAR(I,L) = SQRT(SSDEV(I,L))
COMMM SSDEV(I,L) = XL(I,L)*(AA(I,L)*DSH(I,L)-BB(I,L)*DTH(I,L))
SSDEV(I,L) = XL(I,L)*(AA(I,L)*DSH(I,L)-
1 BB(I,L)*plke(I,L+1)*DTH(I,L))
SSDEV(I,L) = B2 * KHSAVE(I,L) * SSDEV(I,L) * SSDEV(I,L)
TEMP(I,L) = SQRT(SSDEV(I,L))
SVAR(I,L) = (1./2.) * (SVAR(I,L) + TEMP(I,L))
IF ( SVAR(I,L).LT.Z1PEM25) SVAR(I,L) = Z1PEM25
ENDDO
ENDDO
DO L =1,NLEVM1
DO I =1,irun
Q1M(I,L) = SBAR(I,L) / SVAR(I,L)
FCC(I,L) = (1./2.) * ( 1. + ERRF( P5SR*Q1M(I,L) ) )
SHL(I,L) = FCC(I,L) * SBAR(I,L)
ARG(I,L) = (1./2.)*Q1M(I,L)*Q1M(I,L)
IF(ARG(I,L).LE.ARGMAX)
1 SHL(I,L) = SHL(I,L)+RSQ2PI*SVAR(I,L)*EXP(-ARG(I,L))
BETAT(I,L) = 1. + VIRTCON * SH(I,L) - VRT1CON * SHL(I,L)
BETAW(I,L) = VIRTCON *
1 ( TH(I,L) + (CLH/plk(I,L)) * SHL(I,L) )
BETAL(I,L) = ( 1. + VIRTCON*SH(I,L) - TWO*VRT1CON*SHL(I,L) )
1 * (CLH/plke(I,L+1)) - VRT1CON * TH(I,L)
COMMM BETAT1(I,L) = BETAT(I,L) - BB(I,L) * FCC(I,L) * BETAL(I,L)
BETAT1(I,L) = BETAT(I,L) -
1 BB(I,L)*plk(i,L) * FCC(I,L) * BETAL(I,L)
BETAW1(I,L) = BETAW(I,L) + AA(I,L) * FCC(I,L) * BETAL(I,L)
ENDDO
ENDDO
DO L =1,NLEVM1
DO I =1,irun
DTHV(I,L) = (1./2.)*((BETAT1(I,L)+BETAT1(I,L+1))*DTH(I,L)
1 + (BETAW1(I,L)+BETAW1(I,L+1))*DSH(I,L))
ENDDO
ENDDO
C GRADIENTS AT THE TOP OF THE SURFACE LAYER
C -----------------------------------------
DO 9102 I =1,irun
DU(I,NLEV) = CU(I)*XX(I)*AHS(I)*RVK
DV(I,NLEV) = V(I,NLEV) * DU(I,NLEV)
DU(I,NLEV) = U(I,NLEV) * DU(I,NLEV)
DTHV(I,NLEV) = CT(I) * YY(I) *
1 ((THV(I,NLEV)-THV(I,NLEVP1)) * RVK)* AHS(I)
9102 CONTINUE
C CALCULATE BRUNT-VAISALA FREQUENCIES, SHEARS, RICHARDSON NUMBERS
C ---------------------------------------------------------------
DO L =1,NLEV
DO I =1,irun
STRT(I,L) = CP * DTHV(I,L) * DPK(I,L)
DW2(I,L) = DU(I,L) * DU(I,L) + DV(I,L) * DV(I,L)
IF ( DW2(I,L) .LE. 1.e-4 ) DW2(I,L) = 1.e-4
RI(I,L) = STRT(I,L) / DW2(I,L)
ENDDO
ENDDO
C FILL RICHARDSON NUMBER AND SURFACE WIND DIAGNOSTICS
C (THOSE NEEDED FROM FIRST TRBFLX ITERATION)
C ---------------------------------------------------
DO L =1,NLEVM1
DO I =1,irun
SRI(I,L) = RI(I,L)
ENDDO
ENDDO
DO 9108 I =1,irun
SRI(I,NLEV) = RIB(I)
9108 CONTINUE
DO 9110 I =1,irun
SWINDS(I) = WS(I)
9110 CONTINUE
C INITIALIZE KH, KM, QE AND P3 AND ELIMINATE SMALL QQ
C ---------------------------------------------------
DO L =1,NLEVM1
DO I =1,irun
KH(I,L) = 0.
KM(I,L) = 0.
QQE(I,L) = 0.
QE(I,L) = 0.
P3(I,L) = 0.
IBITSTB(I,L) = 0
IF ( QQ(I,L) .GT. 1.e-8 ) THEN
INTQ(I,L) = 1
ELSE
INTQ(I,L) = 0
ENDIF
IF ( QQ(I,L).LE.1.e-8 ) THEN
QQ(I,L) = 0.
Q(I,L) = 0.
ENDIF
ENDDO
ENDDO
DO 300 LMINQ = 1,NLEVM1
IBIT = 0
DO 9116 I = 1,irun
IF ( QQ(I,LMINQ).GT.1.e-8 ) IBIT = IBIT + 1
9116 CONTINUE
IF(IBIT.GE.1)GO TO 310
300 CONTINUE
LMINQ = NLEV-1
310 CONTINUE
LMINQ = 1
LMINQ1 = 1
IF(LMINQ.GT.1)LMINQ1 = LMINQ - 1
C LENGTH SCALE
C ------------
CALL TRBLEN(STRT,DW2,DZ3,Q,VKZE,VKZM,DTHV,DPK,DU,DV,XL,QXLM,
1 NLEV,INIT,LMIN,LMINQ,LMINQ1,CP,INT1,INT2,
2 DZITRP,STBFCN,XL0,Q1,WRKIT1,WRKIT2,WRKIT3,WRKIT4,irun)
C QE AND DIMENSIONLESS COEFFS FROM LEVEL 2 MODEL
C ----------------------------------------------
IF( LMIN .LT. NLEV ) THEN
NLEVML = NLEV - LMIN
CALL TRBL20(RI(1,LMIN),STRT(1,LMIN),DW2(1,LMIN),XL(1,LMIN),
1 KM(1,LMIN),KH(1,LMIN),QE(1,LMIN),QQE(1,LMIN),IBITSTB(1,LMIN),
2 NLEVML,nlev,irun)
ENDIF
C FOR INITIAL START ONLY : USE EQUILIBRIUM MODEL
C ----------------------------------------------
IF ( INIT .EQ. 1 ) THEN
DO L =1,NLEVM1
DO I =1,irun
QQ(I,L) = QQE(I,L)
Q(I,L) = QE(I,L)
ENDDO
ENDDO
INIT = 2
CALL TRBLEN(STRT,DW2,DZ3,Q,VKZE,VKZM,DTHV,DPK,DU,DV,XL,QXLM,
1 NLEV,INIT,LMIN,LMINQ,LMINQ1,CP,INT1,INT2,
2 DZITRP,STBFCN,XL0,Q1,WRKIT1,WRKIT2,WRKIT3,WRKIT4,irun)
INIT = 0
GO TO 550
ENDIF
C DIMENSIONLESS COEFFS AND P3 (Q LE QE)
C -------------------------------------
IF( LMIN .LT. NLEV ) THEN
ISTNML = irun * NLEVML
DO L =LMIN,NLEVM1
DO I =1,irun
IF ( (IBITSTB(I,L).EQ.1) .AND.
1 ( Q(I,L) .LE. QE(I,L) ) ) THEN
IBITSTB(I,L) = 1
ELSE
IBITSTB(I,L) = 0
ENDIF
ENDDO
ENDDO
DO L =LMIN,NLEVM1
DO I =1,irun
IF(IBITSTB(I,L).EQ.1 ) THEN
TEMP(I,L) = Q(I,L) / QE(I,L)
KH(I,L) = TEMP(I,L) * KH(I,L)
KM(I,L) = TEMP(I,L) * KM(I,L)
ENDIF
TEMP(I,L) = 0.01 * QQE(I,L)
IF((IBITSTB(I,L).EQ.1) .AND.
1 ( QQ(I,L) .LE. TEMP(I,L) )) THEN
QQ(I,L) = TEMP(I,L)
Q(I,L) = 0.1 * QE(I,L)
ENDIF
IF(IBITSTB(I,L).EQ.1 ) P3(I,L) = (2.*B3) *
1 ( QE(I,L) - Q(I,L) )
ENDDO
ENDDO
ENDIF
C DIMENSIONLESS COEFFS AND P3 (Q GT QE)
C -------------------------------------
NLEVML = NLEV - LMINQ
CALL TRBL25(Q(1,LMINQ),XL(1,LMINQ),STRT(1,LMINQ),DW2(1,LMINQ),
1 IBITSTB(1,LMINQ),INTQ(1,LMINQ),KM(1,LMINQ),KH(1,LMINQ),
2 P3(1,LMINQ),NLEVML,nlev,irun)
C CALCULATE SOURCE TERM P3
C ------------------------
IF ( LMINQ .LT. LMIN ) THEN
LMIN = LMINQ
ISTNML = irun * ( NLEV - LMIN )
ENDIF
IF( LMIN .LT. NLEV ) THEN
DO L =LMIN,NLEVM1
DO I =1,irun
P3(I,L) = P3(I,L) * DTAU / XL(I,L)
TEMP(I,L) = QQE(I,L) * P3(I,L)
XQ(I,L) = QQE(I,L) - QQ(I,L)
ENDDO
ENDDO
DO L =LMIN,NLEVM1
DO I =1,irun
IF( ( (IBITSTB(I,L).EQ.1) .AND.
1 ( XQ(I,L) .LT. TEMP(I,L) ) )
2 .OR.
3 ( (IBITSTB(I,L).EQ.0) .AND.
4 ( XQ(I,L) .GT. TEMP(I,L) ) ) )
5 P3(I,L) = XQ(I,L) / QQE(I,L)
ENDDO
ENDDO
ENDIF
550 CONTINUE
C DIAGNOSTIC PROFILES : INITIAL RI AND QQ
C ---------------------------------------
IF ( TPROF .AND. FIRST ) THEN
DO 9118 I =1,irun
RIBIN(I) = RIB(I)
CUIN(I) = CU(I)
CTIN(I) = CT(I)
USTARIN(I) = USTAR(I)
RHOSIN(I) = RHODZ2(I,NLEV)
Z0IN(I) = Z0(I)
ZETAIN(I) = ZETA(I)
9118 CONTINUE
DO L =1,NLEV
DO I =1,irun
RIINIT(I,L) = RI(I,L)
QQINIT(I,L) = QQ(I,L)
ENDDO
ENDDO
ENDIF
C Zero diffusion of TKE above 10 mb
do L = 1,nlev-1
do i = 1,irun
c if(pl(i,L).le.10.) then
c qxlm(i,L) = 0.
c endif
if(pl(i,L).le.150.) then
qxlm(i,L) = min(qxlm(i,L),5. _d 0)
endif
enddo
enddo
C UPDATE TURBULENT KINETIC ENERGY QQ
C ----------------------------------
NLEVMQ = NLEV - LMINQ1
ISTNMQ = irun * NLEVMQ
DO L =LMINQ1,NLEVM1
DO I =1,irun
RHOKDZ(I,L) = RHODZ1(I,L) * QXLM(I,L)
ENDDO
ENDDO
CALL TRBDIF(QQ(1,LMINQ1),P3(1,LMINQ1),RHOKDZ(1,LMINQ1),
1 FLXFCE(1,LMINQ1),DTHS,DELTHS,NLEVMQ,1,1.0 _d -20,irun)
C SAVE KH BEFORE ADDING DIMENSIONS FOR USE BY MOIST BOUYANCY CALCULATION
DO L =1,NLEVM1
DO I =1,irun
KHSAVE(I,L) = KH(I,L)
ENDDO
ENDDO
C DIMENSIONAL DIFFUSION COEFFS INCLUDING BACKGROUND AMOUNTS
IF(LMINQ1.GT.1)THEN
ISTLMQ = irun * (LMINQ1-1)
CB
DO L =1,LMINQ1-1
DO I =1,irun
KM(I,L) = KMBG
KH(I,L) = KHBG
ENDDO
ENDDO
CE
ENDIF
CB
DO L =LMINQ1,NLEVM1
DO I =1,irun
Q(I,L) = 2. * QQ(I,L)
Q(I,L) = SQRT(Q(I,L))
XQ(I,L)= XL(I,L) * Q(I,L)
KM(I,L)=XQ(I,L)*KM(I,L)+KMBG
KH(I,L)=XQ(I,L)*KH(I,L)+KHBG
ENDDO
ENDDO
CE
C Zero diffusion of u,v,t,q above 10 mb
do L = 1,nlev-1
do i = 1,irun
c if(pl(i,L).le.10.) then
c kh(i,L) = 0.
c km(i,L) = 0.
c endif
if(pl(i,L).le.150.) then
kh(i,L) = min(kh(i,L),5. _d 0)
km(i,L) = min(km(i,L),5. _d 0)
endif
enddo
enddo
C CALCULATE INTERNAL FLUXES AND UPDATE PROGNOSTIC VARIABLES: TH AND S
DO L =1,NLEV
DO I =1,irun
TEMP(I,L) = RHOZPK(I,L) * KH(I,L)
DELTH(I,L) = 0.
ENDDO
ENDDO
DO 9132 I =1,irun
DELTH(I,NLEVP1) = 1.
9132 CONTINUE
CALL TRBDIF(TH,DELTH,TEMP,FLXFPK,DTHS,DELTHS,NLEV,2,0. _d 0,irun)
do i = 1,irun
hsturb(i) = -1.* dths(i)
dhsdtc(i) = -1.* delths(i)
enddo
do L = 1,nlev
do i = 1,irun
dthdthg(i,L) = delth(i,L)
enddo
enddo
do L = 1,nlev
do i = 1,irun
transth(i,L) = temp(i,L)
enddo
enddo
DO L =1,NLEV
DO I =1,irun
RHOKDZ(I,L) = RHODZ2(I,L) * KH(I,L)
DELSH(I,L) = 0.
ENDDO
ENDDO
DO 9140 I =1,irun
DELSH(I,NLEVP1) = 1.
9140 CONTINUE
CALL TRBDIF(SH,DELSH,RHOKDZ,FLXFAC,DTHL,DELTHL,NLEV,
& 2,0. _d 0,irun)
do i = 1,irun
eturb(i) = -1.* dthl(i)
dedqa(i) = -1.* delthl(i)
enddo
do L = 1,nlev
do i = 1,irun
dshdshg(i,L) = delsh(i,L)
enddo
enddo
do L = 1,nlev
do i = 1,irun
transsh(i,L) = rhokdz(i,L)
enddo
enddo
C Update Tracers Due to Turbulent Diffusion
do i = 1,irun
rhokdz(i,nlev) = 0.0
enddo
c do nt = 1,ntrace
c do i = 1,irun
c tracers(i,nlev+1,nt) = tracers(i,nlev,nt)
c enddo
c CALL TRBDIF(tracers(1,1,nt),DELSH,RHOKDZ,FLXFAC,DTHL,DELTHL,
c . NLEV,4,0. _d 0,irun)
c enddo
C CALCULATE INTERNAL FLUXES AND UPDATE PROGNOSTIC VARIABLES: U AND V
DO L =1,NLEV
DO I =1,irun
RHOKDZ(I,L) = RHODZ2(I,L) * KM(I,L)
ENDDO
ENDDO
CALL TRBDIF(U,V,RHOKDZ,FLXFAC,DTHS,DELTHS,NLEV,3,0. _d 0,irun)
C ( FILL DIAGNOSTIC ARRAYS IF REQUIRED )
DO L =1,NLEV
DO I =1,irun
WU(I,L) = WU(I,L) + RHOKDZ(I,L) * ( U(I,L+1) - U(I,L) )
WV(I,L) = WV(I,L) + RHOKDZ(I,L) * ( V(I,L+1) - V(I,L) )
ENDDO
ENDDO
DO L =1,NLEVM1
DO I =1,irun
IF ( QQ(I,L) .GT. QQMIN ) THEN
IBITSTB(I,L) = 1
ELSE
IBITSTB(I,L) = 0
ENDIF
IF( IBITSTB(I,L).EQ.1 ) FREQDG(I,L) = FREQDG(I,L) + aitr
ENDDO
ENDDO
do i = 1,irun
qqcolmin(i) = qq(i,nlev)*0.1
qqcolmax(i) = qq(i,nlev)
levpbl(i) = nlev
enddo
DO L = nlev-1,1,-1
DO I = 1,irun
IF ( (qq(i,l).gt.qqcolmax(I)).and.(levpbl(i).eq.nlev))then
qqcolmax(i) = qq(i,l)
qqcolmin(i) = 0.1*qqcolmax(I)
endif
if((qq(i,l).lt.qqcolmin(i)).and.(levpbl(i).eq.nlev))
1 levpbl(i)=l
enddo
enddo
do i = 1,irun
lp = levpbl(i)
if(lp.lt.nlev)then
pbldpth(I) = pbldpth(I) + ( (PLE(I,nlev+1)-PLE(I,Lp+2)) +
1 ( (ple(i,lp+2)-ple(i,lp+1))*(qq(i,lp+1)-qqcolmin(i))
2 / (qq(i,lp+1)-qq(i,lp)) ) ) * aitr
else
pbldpth(I) = pbldpth(I) + ( (PLE(I,nlev+1)-PLE(I,2)) +
1 ( (ple(i,2)-ple(i,1))*(qq(i,1)-qqcolmin(i))
2 / qq(i,1) ) ) * aitr
endif
enddo
do i=1,irun
sustar(i) = sustar(i) + aitr*ustar(i)
enddo
do i=1,irun
sz0(i) = sz0(i) + aitr*z0(i)
enddo
DO L =1,NLEV
DO I =1,irun
EU(I,L) = EU(I,L) + AITR*KM(I,L)
ET(I,L) = ET(I,L) + AITR*KH(I,L)
ENDDO
ENDDO
DO I =1,irun
scu(I) = scu(I) + AITR*cu(I)
enddo
DO I =1,irun
sct(I) = sct(I) + AITR*ct(I)
enddo
IF(tprof) then
DO L =1,NLEVM1
DO I =1,irun
XLDIAG(I,L) = XLDIAG(I,L) + AITR*XL(I,L)
ENDDO
ENDDO
endif
FIRST = .FALSE.
C SAVE XL,CT,XX,YY,ZETA FOR USE BY MOIST BOUYANCY CALCULATION
IF(ITER.EQ.ITRTRB)THEN
DO L =1,NLEVM1
DO I =1,irun
XLSAVE(I,L) = XL(I,L)
ENDDO
ENDDO
DO I = 1,irun
CTSAVE(I) = CT(I)
XXSAVE(I) = XX(I)
YYSAVE(I) = YY(I)
ZETASAVE(I) = ZETA(I)
ENDDO
ENDIF
DO L =1,NLEV
DO I =1,irun
turbfcc(I,L) = turbfcc(I,L) + fcc(I,L) * aitr
qliq(I,L) = qliq(I,L) + shl(I,L) * aitr
ENDDO
ENDDO
C END OF MAIN LOOP
2000 CONTINUE
DO L =1,NLEV
DO I =1,irun
WU(I,L) = WU(I,L) * AITR
WV(I,L) = WV(I,L) * AITR
ENDDO
ENDDO
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
SUBROUTINE SFCFLX(NN,VUS,VVS,VTHV1,VTHV2,VTH1,VTH2,VSH1,
1 VSH2,VPK,VPKE,VPE,VZ0,IVWATER,VHS,
2 VAHS,FIRST,LAST,N,IRUN,aitr,VRHO,VRHOZPK,VKH,VKM,
3 VUSTAR,VXX,VYY,VCU,VCT,VRIB,VZETA,VWS,
4 stu2m,stv2m,stt2m,stq2m,stu10m,stv10m,stt10m,stq10m,
5 cp,rgas,undef,
6 lwater, ivbitrib,
7 VHZ,VPSIM,VAPSIM,VPSIG,VPSIHG,VTEMP,VDZETA,VDZ0,VDPSIM,
8 VDPSIH,VZH,VXX0,VYY0,VAPSIHG,VRIB1,VWS1,VPSIH,VZETAL,
9 VZ0L,VPSIH2,VH0,
3 VX0PSIM,VG,VG0,VR1MG0,VZ2,VDZSEA,VAZ0,VXNUM1,VPSIGB2,VDX,
4 VDXPSIM,VDY,VXNUM2,VDEN,VAWS1,VXNUM3,VXNUM,VDZETA1,VDZETA2,
5 VZCOEF2,VZCOEF1,VTEMPLIN,VDPSIMC,VDPSIHC)
C**********************************************************************
C SUBROUTINE SFCFLX - COMPUTES SURFACE TRANSFER COEFFICIENTS
C - CALLED FROM TRBFLX
C ARGUMENTS ::
C INPUT:
C ------
C US - U - COMPONENT OF SURFACE WIND
C VS - V - COMPONENT OF SURFACE WIND
C THV1 - VIRTUAL POTENTIAL TEMPERATURE AT NLAY
C THV2 - VIRTUAL POTENTIAL TEMPERATURE AT GROUND
C TH1 - POTENTIAL TEMPERATURE AT NLAY
C TH2 - POTENTIAL TEMPERATURE AT GROUND
C SH1 - SPECIFIC HUMIDITY AT NLAY
C SH2 - SPECIFIC HUMIDITY AT GROUND
C PK - EVEN LEVEL PRESSURE ** KAPPA AT LEVEL NLAY
C PKE - EDGE LEVEL PRESSURE ** KAPPA AT GROUND
C PE - SURFACE PRESSURE
C Z0 - SURFACE ROUGHNESS
C WATER - BIT ARRAY - '1' OVER OCEANS
C HS - DEPTH OF SURFACE LAYER
C AHS - ONE / HS
C FIRST - LOGICAL .TRUE. FOR FIRST TRBFLX ITERATION
C LAST - LOGICAL .TRUE. FOR LAST TRBFLX ITERATION
C N - NUMBER OF SFCFLX ITERATIONS
C OUTPUT:
C -------
C RHO - DENSITY AT 10M HEIGHT
C RHOZPK - RHO * P**K AT THE SURFACE
C KH - HEAT TRANSFER COEFFICIENT (CT*USTAR)
C KM - MOMENTUM TRANSFER COEFFICIENT (CU*USTAR)
C USTAR - FRICTION VELOCITY
C XX - PHIM(ZETA) - DIMENSIONLESS WIND SHEAR
C YY - PHIH(ZETA) - DIMENSIONLESS TEMP GRADIENT
C CU - MOMENTUM TRANSPORT COEFFICIENT
C CT - HEAT TRANSPORT COEFFICIENT
C**********************************************************************
implicit none
C Argument List Declarations
integer nn,n,irun
_RL aitr,cp,rgas,undef
_RL VUS(IRUN),VVS(IRUN),VTHV1(IRUN),VTHV2(IRUN)
_RL VTH1(IRUN),VTH2(IRUN),VSH1(IRUN),VSH2(IRUN)
_RL VPK(IRUN),VPKE(IRUN),VPE(IRUN)
_RL VZ0(IRUN),VHS(IRUN),VAHS(IRUN)
integer IVWATER(IRUN)
LOGICAL FIRST,LAST
_RL VRHO(IRUN),VRHOZPK(IRUN)
_RL VKM(IRUN),VKH(IRUN),VUSTAR(IRUN),VXX(IRUN)
_RL VYY(IRUN),VCU(IRUN),VCT(IRUN),VRIB(IRUN)
_RL VZETA(IRUN),VWS(IRUN)
_RL stu2m(irun),stv2m(irun),stt2m(irun),stq2m(irun)
_RL stu10m(irun),stv10m(irun),stt10m(irun),stq10m(irun)
LOGICAL LWATER
integer IVBITRIB(irun)
_RL VHZ(irun),VPSIM(irun),VAPSIM(irun),VPSIG(irun),VPSIHG(irun)
_RL VTEMP(irun),VDZETA(irun),VDZ0(irun),VDPSIM(irun)
_RL VDPSIH(irun),VZH(irun),VXX0(irun),VYY0(irun)
_RL VAPSIHG(irun),VRIB1(irun),VWS1(irun)
_RL VPSIH(irun),VZETAL(irun),VZ0L(irun),VPSIH2(irun),VH0(irun)
_RL VX0PSIM(irun),VG(irun),VG0(irun),VR1MG0(irun)
_RL VZ2(irun),VDZSEA(irun),VAZ0(irun),VXNUM1(irun)
_RL VPSIGB2(irun),VDX(irun),VDXPSIM(irun),VDY(irun)
_RL VXNUM2(irun),VDEN(irun),VAWS1(irun),VXNUM3(irun)
_RL VXNUM(irun),VDZETA1(irun),VDZETA2(irun)
_RL VZCOEF2(irun),VZCOEF1(irun),VTEMPLIN(irun)
_RL VDPSIMC(irun),VDPSIHC(irun)
C Local Variables
_RL USTMX3,USTZ0S,Z0MIN,H0BYZ0,USTH0S,H0VEG,Z0VEGM,PRFAC
_RL XPFAC,DIFSQT
PARAMETER ( USTMX3 = 0.0632456)
PARAMETER ( USTZ0S = 0.2030325E-5)
PARAMETER ( Z0MIN = USTZ0S/USTMX3)
PARAMETER ( H0BYZ0 = 30.0 )
PARAMETER ( USTH0S = H0BYZ0*USTZ0S )
PARAMETER ( H0VEG = 0.01 )
PARAMETER ( Z0VEGM = 0.005 )
PARAMETER ( PRFAC = 0.595864 )
PARAMETER ( XPFAC = .55 )
PARAMETER ( DIFSQT = 3.872983E-3)
_RL psihdiag(irun),psimdiag(irun)
_RL getcon,vk,rvk,vk2,bmdl
integer iwater,itype
integer i,iter
vk = getcon('VON KARMAN')
rvk = 1./vk
vk2 = vk*vk
BMDL = VK * XPFAC * PRFAC / DIFSQT
C DETERMINE SURFACE WIND MAGNITUDE AND BULK RICHARDSON NUMBER
DO 9000 I = 1,IRUN
VWS(I) = VUS(I) * VUS(I) + VVS(I) * VVS(I)
IF ( VWS(I) .LE. 1.e-4) VWS(I) = 1.e-4
VRIB(I) = ( CP * (VPKE(I)-VPK(I)) ) *
1 (VTHV1(I) - VTHV2(I)) / VWS(I)
VWS(I) = SQRT( VWS(I) )
9000 CONTINUE
C INITIALIZATION (FIRST TRBFLX ITERATION)
C INITIAL GUESS FOR ROUGHNESS LENGTH Z0 OVER WATER
IF (.NOT. FIRST) GO TO 100
IWATER = 0
DO 9002 I = 1,IRUN
IF (IVWATER(I).EQ.1) IWATER = IWATER + 1
9002 CONTINUE
LWATER = .FALSE.
IF(IWATER.GE.1)LWATER = .TRUE.
IF(LWATER)THEN
DO 9004 I = 1,IRUN
IF (IVWATER(I).EQ.1) VZ0(I) = 0.0003
9004 CONTINUE
ENDIF
do i = 1,irun
vh0(i) = h0byz0 * vz0(i)
if(vz0(i).ge.z0vegm)vh0(i) = h0veg
enddo
C CU AND PSIHG FOR NEUTRALLY STRATIFIED FLOW
DO 9006 I = 1,IRUN
VHZ(I) = VHS(I) / VZ0(I)
VPSIM(I) = LOG( VHZ(I) )
VAPSIM(I) = 1. / VPSIM(I)
VCU(I) = VK * VAPSIM(I)
VUSTAR(I) = VCU(I) * VWS(I)
VPSIG(I) = VH0(I) * VUSTAR(I) - USTH0S
if(VPSIG(I).lt.0.) VPSIG(I) = 0.
VPSIG(I) = SQRT( VPSIG(I) )
VPSIG(I) = BMDL * VPSIG(I)
VPSIHG(I) = VPSIM(I) + VPSIG(I)
9006 CONTINUE
C LINEAR CORRECTION FOR ERROR IN ROUGHNESS LENGTH Z0
IF(LWATER)THEN
DO 9008 I = 1,IRUN
VTEMP(I) = 0.
9008 CONTINUE
CALL LINADJ(NN,VRIB,VRIB,VWS,
1 VWS,VZ0,VUSTAR,IVWATER,
2 VAPSIM,VTEMP,VTEMP,
3 VTEMP,VTEMP,VTEMP,
4 VTEMP,VTEMP,1,.TRUE.,IRUN,VDZETA,
5 VDZ0,VDPSIM,VDPSIH,
6 IVBITRIB,
3 VX0PSIM,VG,VG0,VR1MG0,VZ2,VDZSEA,VAZ0,VXNUM1,VPSIGB2,VDX,
4 VDXPSIM,VDY,VXNUM2,VDEN,VAWS1,VXNUM3,VXNUM,VDZETA1,VDZETA2,
5 VZCOEF2,VZCOEF1,VTEMPLIN,VDPSIMC,VDPSIHC)
DO 9010 I = 1,IRUN
IF ( IVWATER(I).EQ.1 ) THEN
VCU(I) = VCU(I) * (1. - VDPSIM(I)*VAPSIM(I))
VZ0(I) = VZ0(I) + VDZ0(I)
ENDIF
IF ( IVWATER(I).EQ.1) THEN
IF ( VZ0(I) .LE. Z0MIN ) VZ0(I) = Z0MIN
vh0(i) = h0byz0 * vz0(i)
VPSIG(I) = VH0(I) * VCU(I) * VWS(I) - USTH0S
if(VPSIG(I).lt.0.) VPSIG(I) = 0.
VPSIG(I) = SQRT( VPSIG(I) )
VPSIG(I) = BMDL * VPSIG(I)
VPSIHG(I) = VPSIM(I) + VDPSIH(I) + VPSIG(I)
ENDIF
9010 CONTINUE
ENDIF
C INITIAL GUESS FOR STABILITY PARAMETER ZETA
DO 9012 I = 1,IRUN
VZETA(I) = VK2 * VRIB(I) / (VCU(I) * VCU(I) * VPSIHG(I))
9012 CONTINUE
C RECOMPUTE CU, ESTIMATE PSIHG AND UPDATE ZETA AND Z0
DO 9014 I = 1,IRUN
VZH(I) = VZ0(I) * VAHS(I)
9014 CONTINUE
CALL PSI (VZETA,VZH,VPSIM,
1 VTEMP,IRUN,VXX,VXX0,VYY,
2 VYY0,2)
DO 9016 I = 1,IRUN
VCU(I) = VK / VPSIM(I)
VPSIG(I) = VH0(I) * VCU(I) * VWS(I) - USTH0S
if(VPSIG(I).lt.0.) VPSIG(I) = 0.
VPSIG(I) = SQRT(VPSIG(I))
VPSIG(I) = BMDL * VPSIG(I)
VPSIHG(I) = VPSIM(I) + VPSIG(I)
VZETA(I) = VK2 * VRIB(I) / (VCU(I) * VCU(I) * VPSIHG(I))
9016 CONTINUE
IF(LWATER)THEN
CCCOOOMMMM ADDED 'WHERE WATER'
DO 9018 I = 1,IRUN
IF (IVWATER(I).EQ.1) VUSTAR(I) = VCU(I) * VWS(I)
9018 CONTINUE
CALL ZCSUB ( VUSTAR,VHZ,IVWATER,.FALSE.,IRUN,VTEMP)
DO 9020 I = 1,IRUN
IF (IVWATER(I).EQ.1 ) then
VZ0(I) = VTEMP(I)
IF ( VZ0(I) .LE. Z0MIN ) VZ0(I) = Z0MIN
vh0(i) = h0byz0 * vz0(i)
endif
9020 CONTINUE
ENDIF
GO TO 125
C LINEARLY UPDATE ZETA AND Z0 FOR SECOND OR GREATER TRBFLX ITERATION
100 CONTINUE
CALL LINADJ(NN,VRIB1,VRIB,VWS1,
1 VWS,VZ0,VUSTAR,IVWATER,
2 VAPSIM,VAPSIHG,VPSIH,
3 VPSIG,VXX,VXX0,
4 VYY,VYY0,2,LWATER,IRUN,VDZETA,
5 VDZ0,VDPSIM,VDPSIH,
6 IVBITRIB,
3 VX0PSIM,VG,VG0,VR1MG0,VZ2,VDZSEA,VAZ0,VXNUM1,VPSIGB2,VDX,
4 VDXPSIM,VDY,VXNUM2,VDEN,VAWS1,VXNUM3,VXNUM,VDZETA1,VDZETA2,
5 VZCOEF2,VZCOEF1,VTEMPLIN,VDPSIMC,VDPSIHC)
DO 9022 I = 1,IRUN
VZETA(I) = VZETA(I) + VZETAL(I) * VDZETA(I)
IF (IVBITRIB(I).EQ.1 )VZETA(I) =
1 VPSIM(I) * VPSIM(I) * VRIB(I) * VCT(I) * RVK
9022 CONTINUE
IF ( LWATER ) THEN
DO 9024 I = 1,IRUN
IF (IVWATER(I).EQ.1 ) then
VZ0(I) = VZ0(I) + VZ0L(I) * VDZ0(I)
IF (VZ0(I) .LE. Z0MIN ) VZ0(I) = Z0MIN
vh0(i) = h0byz0 * vz0(i)
endif
9024 CONTINUE
ENDIF
125 CONTINUE
C ITERATIVE LOOP - N ITERATIONS
C COMPUTE CU AND CT
DO 200 ITER = 1,N
DO 9026 I = 1,IRUN
VZH(I) = VZ0(I) * VAHS(I)
9026 CONTINUE
CALL PSI (VZETA,VZH,VPSIM,
1 VPSIH,IRUN,VXX,VXX0,VYY,
2 VYY0,1)
DO 9028 I = 1,IRUN
VCU(I) = VK / VPSIM(I)
VUSTAR(I) = VCU(I) * VWS(I)
VPSIG(I) = VH0(I) * VUSTAR(I) - USTH0S
if(VPSIG(I).lt.0.) VPSIG(I) = 0.
VPSIG(I) = SQRT(VPSIG(I))
VPSIG(I) = BMDL * VPSIG(I)
VPSIHG(I) = VPSIH(I) + VPSIG(I)
C LINEAR CORRECTIONS FOR CU, CT, ZETA, AND Z0
VAPSIM(I) = VCU(I) * RVK
VAPSIHG(I) = 1. / VPSIHG(I)
VRIB1(I) = VAPSIM(I) * VAPSIM(I) * VPSIHG(I) * VZETA(I)
9028 CONTINUE
ITYPE = 3
IF(ITER.EQ.N) ITYPE = 4
IF( (ITYPE.EQ.4) .AND. (.NOT.LAST) ) ITYPE = 5
CALL LINADJ(NN,VRIB1,VRIB,VWS,
1 VWS,VZ0,VUSTAR,IVWATER,
2 VAPSIM,VAPSIHG,VPSIH,
3 VPSIG,VXX,VXX0,
4 VYY,VYY0,ITYPE,LWATER,IRUN,VDZETA,
5 VDZ0,VDPSIM,VDPSIH,
6 IVBITRIB,
3 VX0PSIM,VG,VG0,VR1MG0,VZ2,VDZSEA,VAZ0,VXNUM1,VPSIGB2,VDX,
4 VDXPSIM,VDY,VXNUM2,VDEN,VAWS1,VXNUM3,VXNUM,VDZETA1,VDZETA2,
5 VZCOEF2,VZCOEF1,VTEMPLIN,VDPSIMC,VDPSIHC)
C UPDATES OF ZETA, Z0, CU AND CT
IF (ITYPE.EQ.5) THEN
DO 9030 I = 1,IRUN
VZETAL(I) = VZETA(I)
VZ0L(I) = VZ0(I)
9030 CONTINUE
ENDIF
DO 9032 I = 1,IRUN
VZETA(I) = VZETA(I) * ( 1. + VDZETA(I) )
IF (IVBITRIB(I).EQ.1 ) VZETA(I) =
1 VPSIM(I) * VPSIM(I) * VRIB(I) * VAPSIHG(I)
9032 CONTINUE
IF ( LWATER ) THEN
DO 9034 I = 1,IRUN
IF (IVWATER(I).EQ.1 ) then
VZ0(I) = VZ0(I) * ( 1. + VDZ0(I) )
IF (VZ0(I) .LE. Z0MIN ) VZ0(I) = Z0MIN
vh0(i) = h0byz0 * vz0(i)
endif
9034 CONTINUE
ENDIF
IF ( ITER .EQ. N ) THEN
DO 9036 I = 1,IRUN
VPSIM(I) = VPSIM(I) + VDPSIM(I)
VCU(I) = VK / VPSIM(I)
VUSTAR(I) = VCU(I) * VWS(I)
VPSIG(I) = VH0(I) * VUSTAR(I) - USTH0S
if(VPSIG(I).lt.0.) VPSIG(I) = 0.
VPSIG(I) = SQRT(VPSIG(I))
VPSIG(I) = BMDL * VPSIG(I)
VPSIHG(I) = VPSIH(I) + VDPSIH(I) + VPSIG(I)
VCT(I) = VK / VPSIHG(I)
9036 CONTINUE
ENDIF
C SAVE VALUES OF RIB AND WS FOR NEXT ITERATION OF TRBFLX
IF (ITYPE.EQ.5) THEN
DO 9038 I = 1,IRUN
VRIB1(I) = VRIB(I)
VWS1(I) = VWS(I)
9038 CONTINUE
ENDIF
200 CONTINUE
C CALCULATE RHO-SURFACE ( KG / M**3 )
IF (FIRST) THEN
DO I = 1,IRUN
VTEMP(I) = 10. * VAHS(I) * VZETA(I)
VZH(I) = VZ0(I) * 0.1
ENDDO
CALL PSI (VTEMP,VZH,VHZ,
1 VPSIH2,IRUN,VHZ,VHZ,VHZ,
2 VHZ,3)
DO I = 1,IRUN
VTEMP(I) = ( VPSIH2(I) + VPSIG(I) ) / VPSIHG(I)
VRHO(I) = VPKE(I)*( VTH2(I) + VTEMP(I) * (VTH1(I)-VTH2(I)) )
VRHO(I) = VPE(I)*100. / ( RGAS * VRHO(I) )
ENDDO
ENDIF
C interpolate uvtq to 2m and to 10 meters for diagnostic output
C use psih and psim which represent non-dim change from ground
C to specified level
C and multiply theta by surface p**kappa to get temperatures
do i = 1,irun
vtemp(i) = 2. * vahs(i) * vzeta(i)
vzh(i) = vz0(i) * 0.5
if(vz0(i).ge.2.)vzh(i) = 0.9
enddo
call PSI(vtemp,vzh,psimdiag,psihdiag,irun,vhz,vhz,vhz,vhz,1)
do i = 1,irun
stu2m(i) = (psimdiag(i)/vpsim(i) * vus(i))
stv2m(i) = (psimdiag(i)/vpsim(i) * vvs(i))
stt2m(i) = ( (vth2(i) + ((psihdiag(i)+vpsig(i))/vpsihg(i))*
1 (vth1(i)-vth2(i))) ) * vpke(i)
stq2m(i) = (vsh2(i) + ((psihdiag(i)+vpsig(i))/vpsihg(i))*
1 (vsh1(i)-vsh2(i)))
if(vz0(i).ge.2.)then
stu2m(i) = UNDEF
stv2m(i) = UNDEF
stt2m(i) = UNDEF
stq2m(i) = UNDEF
endif
enddo
do i = 1,irun
vtemp(i) = 10. * vahs(i) * vzeta(i)
vzh(i) = vz0(i) * 0.1
enddo
call PSI(vtemp,vzh,psimdiag,psihdiag,irun,vhz,vhz,vhz,vhz,1)
do i = 1,irun
stu10m(i) = (psimdiag(i)/vpsim(i) * vus(i))
stv10m(i) = (psimdiag(i)/vpsim(i) * vvs(i))
stt10m(i) = ( (vth2(i) + ((psihdiag(i)+vpsig(i))/vpsihg(i))*
1 (vth1(i)-vth2(i))) ) * vpke(i)
stq10m(i) = (vsh2(i) + ((psihdiag(i)+vpsig(i))/vpsihg(i))*
1 (vsh1(i)-vsh2(i)))
enddo
C EVALUATE TURBULENT TRANSFER COEFFICIENTS
DO 9044 I = 1,IRUN
VRHOZPK(I) = VRHO(I) * VPKE(I)
VKH(I) = VUSTAR(I) * VCT(I)
VKM(I) = VUSTAR(I) * VCU(I)
9044 CONTINUE
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
SUBROUTINE PHI(Z,PHIM,PHIH,IFLAG,N)
C**********************************************************************
C FUNCTION PHI - SOLVES KEYPS EQUATIONS
C - CALLED FROM PSI
C DESCRIPTION OF PARAMETERS
C Z - INPUTED VALUE OF MONIN- OBUKHOV STABILITY PARAMETER ZETA
C TIMES APPROPRIATE CONSTANT
C PHIM - OUTPUTED SOLUTION OF KEYPS EQUATION FOR MOMENTUM
C PHIH - OUTPUTED SOLUTION OF KEYPS EQUATION FOR SCALARS
C IFLAG - FLAG TO DETERMINE IF X IS NEEDED (IFLAG=2), Y IS NEEDED
C (IFLAG=3), OR BOTH (IFLAG=1)
C N - LENGTH OF VECTOR TO BE SOLVED
C**********************************************************************
implicit none
C Argument List Declarations
integer n,iflag
_RL PHIM(N),PHIH(N),Z(N)
C Local Variables
integer I1(N),I2(N)
_RL ZSTAR(N),E1(N),E2(N),TEMP1(N)
_RL PHIM0(385),ZLINM1(75),ZLINM2(75),ZLINM3(36)
_RL ZLOGM1(74),ZLOGM2(75),ZLOGM3(50)
_RL PHIH0(385),ZLINH1(75),ZLINH2(75),ZLINH3(36)
_RL ZLOGH1(74),ZLOGH2(75),ZLOGH3(50)
EQUIVALENCE (PHIM0(1),ZLINM1(1)),(PHIM0(76),ZLINM2(1))
EQUIVALENCE (PHIM0(151),ZLINM3(1))
EQUIVALENCE (PHIM0(187),ZLOGM1(1)),(PHIM0(261),ZLOGM2(1))
EQUIVALENCE (PHIM0(336),ZLOGM3(1))
EQUIVALENCE (PHIH0(1),ZLINH1(1)),(PHIH0(76),ZLINH2(1))
EQUIVALENCE (PHIH0(151),ZLINH3(1))
EQUIVALENCE (PHIH0(187),ZLOGH1(1)),(PHIH0(261),ZLOGH2(1))
EQUIVALENCE (PHIH0(336),ZLOGH3(1))
DATA ZLOGM1/
& 0.697894,0.678839,0.659598,0.640260,
& 0.620910,0.601628,0.582486,0.563550,0.544877,
& 0.526519,0.508516,0.490903,0.473708,0.456951,
& 0.440649,0.424812,0.409446,0.394553,0.380133,
& 0.366182,0.352695,0.339664,0.327082,0.314938,
& 0.303222,0.291923,0.281029,0.270528,0.260409,
& 0.250659,0.241267,0.232221,0.223509,0.215119,
& 0.207041,0.199264,0.191776,0.184568,0.177628,
& 0.170949,0.164519,0.158331,0.152374,0.146641,
& 0.141123,0.135813,0.130702,0.125783,0.121048,
& 0.116492,0.112107,0.107887,0.103826,0.0999177,
& 0.0961563,0.0925364,0.0890528,0.0857003,0.0824739,
& 0.0793690,0.0763810,0.0735054,0.0707380,0.0680749,
& 0.0655120,0.0630455,0.0606720,0.0583877,0.0561895,
& 0.0540740,0.0520382,0.0500790,0.0481936,0.0463791/
DATA ZLOGM2/
& 0.0446330,0.0429526,0.0413355,0.0397792,0.0382816,
& 0.0368403,0.0354533,0.0341185,0.0328340,0.0315978,
& 0.0304081,0.0292633,0.0281616,0.0271013,0.0260809,
& 0.0250990,0.0241540,0.0232447,0.0223695,0.0215273,
& 0.0207168,0.0199369,0.0191862,0.0184639,0.0177687,
& 0.0170998,0.0164560,0.0158364,0.0152402,0.0146664,
& 0.0141142,0.0135828,0.0130714,0.0125793,0.0121057,
& 0.0116499,0.0112113,0.0107892,0.0103830,0.999210E-2,
& 0.961590E-2,0.925387E-2,0.890547E-2,0.857018E-2,0.824752E-2,
& 0.793701E-2,0.763818E-2,0.735061E-2,0.707386E-2,0.680754E-2,
& 0.655124E-2,0.630459E-2,0.606722E-2,0.583880E-2,0.561897E-2,
& 0.540742E-2,0.520383E-2,0.500791E-2,0.481937E-2,0.463792E-2,
& 0.446331E-2,0.429527E-2,0.413355E-2,0.397793E-2,0.382816E-2,
& 0.368403E-2,0.354533E-2,0.341185E-2,0.328340E-2,0.315978E-2,
& 0.304082E-2,0.292633E-2,0.281616E-2,0.271013E-2,0.260809E-2/
DATA ZLOGM3/
& 0.250990E-2,0.241541E-2,0.232447E-2,0.223695E-2,0.215273E-2,
& 0.207168E-2,0.199369E-2,0.191862E-2,0.184639E-2,0.177687E-2,
& 0.170998E-2,0.164560E-2,0.158364E-2,0.152402E-2,0.146664E-2,
& 0.141142E-2,0.135828E-2,0.130714E-2,0.125793E-2,0.121057E-2,
& 0.116499E-2,0.112113E-2,0.107892E-2,0.103830E-2,0.999210E-3,
& 0.961590E-3,0.925387E-3,0.890547E-3,0.857018E-3,0.824752E-3,
& 0.793701E-3,0.763818E-3,0.735061E-3,0.707386E-3,0.680754E-3,
& 0.655124E-3,0.630459E-3,0.606722E-3,0.583880E-3,0.561897E-3,
& 0.540742E-3,0.520383E-3,0.500791E-3,0.481937E-3,0.463792E-3,
& 0.446331E-3,0.429527E-3,0.413355E-3,0.397793E-3,0.382816E-3/
DATA ZLOGH1/
& 0.640529,0.623728,0.606937,0.590199,
& 0.573552,0.557032,0.540672,0.524504,0.508553,
& 0.492843,0.477397,0.462232,0.447365,0.432809,
& 0.418574,0.404670,0.391103,0.377878,0.364999,
& 0.352468,0.340284,0.328447,0.316954,0.305804,
& 0.294992,0.284514,0.274364,0.264538,0.255028,
& 0.245829,0.236933,0.228335,0.220026,0.211999,
& 0.204247,0.196762,0.189537,0.182564,0.175837,
& 0.169347,0.163088,0.157051,0.151231,0.145620,
& 0.140211,0.134998,0.129974,0.125133,0.120469,
& 0.115975,0.111645,0.107475,0.103458,0.995895E-1,
& 0.958635E-1,0.922753E-1,0.888199E-1,0.854925E-1,0.822886E-1,
& 0.792037E-1,0.762336E-1,0.733739E-1,0.706208E-1,0.679704E-1,
& 0.654188E-1,0.629625E-1,0.605979E-1,0.583217E-1,0.561306E-1,
& 0.540215E-1,0.519914E-1,0.500373E-1,0.481564E-1,0.463460E-1/
DATA ZLOGH2/
& 0.446034E-1,0.429263E-1,0.413120E-1,0.397583E-1,0.382629E-1,
& 0.368237E-1,0.354385E-1,0.341053E-1,0.328222E-1,0.315873E-1,
& 0.303988E-1,0.292550E-1,0.281541E-1,0.270947E-1,0.260750E-1,
& 0.250937E-1,0.241494E-1,0.232405E-1,0.223658E-1,0.215240E-1,
& 0.207139E-1,0.199342E-1,0.191839E-1,0.184618E-1,0.177669E-1,
& 0.170981E-1,0.164545E-1,0.158351E-1,0.152390E-1,0.146653E-1,
& 0.141133E-1,0.135820E-1,0.130707E-1,0.125786E-1,0.121051E-1,
& 0.116494E-1,0.112108E-1,0.107888E-1,0.103826E-1,0.999177E-2,
& 0.961561E-2,0.925360E-2,0.890523E-2,0.856997E-2,0.824733E-2,
& 0.793684E-2,0.763803E-2,0.735048E-2,0.707375E-2,0.680743E-2,
& 0.655114E-2,0.630450E-2,0.606715E-2,0.583873E-2,0.561891E-2,
& 0.540737E-2,0.520379E-2,0.500787E-2,0.481933E-2,0.463789E-2,
& 0.446328E-2,0.429524E-2,0.413353E-2,0.397790E-2,0.382814E-2,
& 0.368401E-2,0.354532E-2,0.341184E-2,0.328338E-2,0.315977E-2,
& 0.304081E-2,0.292632E-2,0.281615E-2,0.271012E-2,0.260809E-2/
DATA ZLOGH3/
& 0.250990E-2,0.241540E-2,0.232446E-2,0.223695E-2,0.215273E-2,
& 0.207168E-2,0.199368E-2,0.191862E-2,0.184639E-2,0.177687E-2,
& 0.170997E-2,0.164559E-2,0.158364E-2,0.152402E-2,0.146664E-2,
& 0.141142E-2,0.135828E-2,0.130714E-2,0.125793E-2,0.121057E-2,
& 0.116499E-2,0.112113E-2,0.107892E-2,0.103830E-2,0.999209E-3,
& 0.961590E-3,0.925387E-3,0.890546E-3,0.857018E-3,0.824752E-3,
& 0.793700E-3,0.763818E-3,0.735061E-3,0.707386E-3,0.680754E-3,
& 0.655124E-3,0.630459E-3,0.606722E-3,0.583880E-3,0.561897E-3,
& 0.540742E-3,0.520383E-3,0.500791E-3,0.481937E-3,0.463792E-3,
& 0.446331E-3,0.429527E-3,0.413355E-3,0.397793E-3,0.382816E-3/
DATA ZLINM1/
& 0.964508,0.962277,0.960062,0.957863,0.955680,
& 0.953512,0.951359,0.949222,0.947100,0.944992,
& 0.942899,0.940821,0.938758,0.936709,0.934673,
& 0.932652,0.930645,0.928652,0.926672,0.924706,
& 0.922753,0.920813,0.918886,0.916973,0.915072,
& 0.913184,0.911308,0.909445,0.907594,0.905756,
& 0.903930,0.902115,0.900313,0.898522,0.896743,
& 0.894975,0.893219,0.891475,0.889741,0.888019,
& 0.886307,0.884607,0.882917,0.881238,0.879569,
& 0.877911,0.876264,0.874626,0.872999,0.871382,
& 0.869775,0.868178,0.866591,0.865013,0.863445,
& 0.861887,0.860338,0.858798,0.857268,0.855747,
& 0.854235,0.852732,0.851238,0.849753,0.848277,
& 0.846809,0.845350,0.843900,0.842458,0.841025,
& 0.839599,0.838182,0.836774,0.835373,0.833980/
DATA ZLINM2/
& 0.832596,0.831219,0.829850,0.828489,0.827136,
& 0.825790,0.824451,0.823121,0.821797,0.820481,
& 0.819173,0.817871,0.816577,0.815289,0.814009,
& 0.812736,0.811470,0.810210,0.808958,0.807712,
& 0.806473,0.805240,0.804015,0.802795,0.801582,
& 0.800376,0.799176,0.797982,0.796794,0.795613,
& 0.794438,0.793269,0.792106,0.790949,0.789798,
& 0.788652,0.787513,0.786380,0.785252,0.784130,
& 0.783014,0.781903,0.780798,0.779698,0.778604,
& 0.777516,0.776432,0.775354,0.774282,0.773215,
& 0.772153,0.771096,0.770044,0.768998,0.767956,
& 0.766920,0.765888,0.764862,0.763840,0.762824,
& 0.761812,0.760805,0.759803,0.758805,0.757813,
& 0.756824,0.755841,0.754862,0.753888,0.752918,
& 0.751953,0.750992,0.750035,0.749083,0.748136/
DATA ZLINM3/
& 0.747192,0.746253,0.745318,0.744388,0.743462,
& 0.742539,0.741621,0.740707,0.739798,0.738892,
& 0.737990,0.737092,0.736198,0.735308,0.734423,
& 0.733540,0.732662,0.731788,0.730917,0.730050,
& 0.729187,0.728328,0.727472,0.726620,0.725772,
& 0.724927,0.724086,0.723248,0.722414,0.721584,
& 0.720757,0.719933,0.719113,0.718296,0.717483,
& 0.716673/
DATA ZLINH1/
& 0.936397,0.932809,0.929287,0.925827,0.922429,
& 0.919089,0.915806,0.912579,0.909405,0.906284,
& 0.903212,0.900189,0.897214,0.894284,0.891399,
& 0.888558,0.885759,0.883001,0.880283,0.877603,
& 0.874962,0.872357,0.869788,0.867255,0.864755,
& 0.862288,0.859854,0.857452,0.855081,0.852739,
& 0.850427,0.848144,0.845889,0.843662,0.841461,
& 0.839287,0.837138,0.835014,0.832915,0.830841,
& 0.828789,0.826761,0.824755,0.822772,0.820810,
& 0.818869,0.816949,0.815050,0.813170,0.811310,
& 0.809470,0.807648,0.805845,0.804060,0.802293,
& 0.800543,0.798811,0.797095,0.795396,0.793714,
& 0.792047,0.790396,0.788761,0.787141,0.785535,
& 0.783945,0.782369,0.780807,0.779259,0.777724,
& 0.776204,0.774696,0.773202,0.771720,0.770251/
DATA ZLINH2/
& 0.768795,0.767351,0.765919,0.764499,0.763091,
& 0.761694,0.760309,0.758935,0.757571,0.756219,
& 0.754878,0.753547,0.752226,0.750916,0.749616,
& 0.748326,0.747045,0.745775,0.744514,0.743262,
& 0.742020,0.740787,0.739563,0.738348,0.737141,
& 0.735944,0.734755,0.733574,0.732402,0.731238,
& 0.730083,0.728935,0.727795,0.726664,0.725539,
& 0.724423,0.723314,0.722213,0.721119,0.720032,
& 0.718952,0.717880,0.716815,0.715756,0.714704,
& 0.713660,0.712621,0.711590,0.710565,0.709547,
& 0.708534,0.707529,0.706529,0.705536,0.704549,
& 0.703567,0.702592,0.701623,0.700660,0.699702,
& 0.698750,0.697804,0.696863,0.695928,0.694998,
& 0.694074,0.693155,0.692241,0.691333,0.690430,
& 0.689532,0.688639,0.687751,0.686868,0.685990/
DATA ZLINH3/
& 0.685117,0.684249,0.683386,0.682527,0.681673,
& 0.680824,0.679979,0.679139,0.678303,0.677472,
& 0.676645,0.675823,0.675005,0.674191,0.673381,
& 0.672576,0.671775,0.670978,0.670185,0.669396,
& 0.668611,0.667830,0.667054,0.666281,0.665512,
& 0.664746,0.663985,0.663227,0.662473,0.661723,
& 0.660977,0.660234,0.659495,0.658759,0.658027,
& 0.657298/
integer ibit1,ibit2,i
IBIT1 = 0
IBIT2 = 0
DO 9000 I = 1,N
IF(Z(I).GE.0.15)IBIT1 = IBIT1 + 1
IF(Z(I).GT.2. )IBIT2 = IBIT2 + 1
9000 CONTINUE
IF( IBIT1 .LE. 0 ) GO TO 200
DO 9002 I = 1,N
ZSTAR(I) = 100. * Z(I) - 14.
9002 CONTINUE
IF( IBIT2 .LE. 0 ) GO TO 60
DO 9004 I = 1,N
TEMP1(I) = Z(I)*0.5
IF( Z(I) .LE. 2. )TEMP1(I) = 1.
TEMP1(I) = LOG10(TEMP1(I))
TEMP1(I) = (TEMP1(I) + 9.3) * 20.
IF( Z(I) .GT. 2. ) ZSTAR(I) = TEMP1(I)
IF( Z(I).GT.1.78e10 ) ZSTAR(I) = 384.9999
9004 CONTINUE
60 CONTINUE
DO 9006 I = 1,N
I1(I) = ZSTAR(I)
I2(I) = I1(I) + 1
TEMP1(I) = ZSTAR(I) - I1(I)
9006 CONTINUE
IF( IFLAG .GT. 2 ) GO TO 100
DO 9008 I = 1,N
if( z(i).ge.0.15 ) then
E1(I) = PHIM0( I1(I) )
E2(I) = PHIM0( I2(I) )
PHIM(I) = TEMP1(I) * ( E2(I)-E1(I) )
PHIM(I) = PHIM(I) + E1(I)
endif
9008 CONTINUE
100 CONTINUE
IF( IFLAG .EQ. 2 ) GO TO 200
DO 9010 I = 1,N
if( z(i).ge.0.15 ) then
E1(I) = PHIH0( I1(I) )
E2(I) = PHIH0( I2(I) )
PHIH(I) = TEMP1(I) * ( E2(I)-E1(I) )
PHIH(I) = PHIH(I) + E1(I)
endif
9010 CONTINUE
200 CONTINUE
IF( IBIT1 .GE. N ) GO TO 500
DO 9012 I = 1,N
ZSTAR(I) = -Z(I)
9012 CONTINUE
IF( IFLAG .GT. 2 ) GO TO 300
DO 9014 I = 1,N
IF( Z(I) .LT. 0.15 ) PHIM(I) = 1. + ZSTAR(I)
2 *(0.25+ZSTAR(I)*(0.09375+ZSTAR(I)*
3 (0.03125+0.00732422 * ZSTAR(I))))
9014 CONTINUE
300 CONTINUE
IF( IFLAG .EQ. 2 ) GO TO 500
DO 9016 I = 1,N
IF( Z(I) .LT. 0.15 ) THEN
PHIH(I) =1.+ Z(I) * (0.5+ZSTAR(I)*(0.375+ZSTAR(I)*
1 (0.5+ZSTAR(I)*(0.8203125+ZSTAR(I)*
2 (1.5+2.93262*ZSTAR(I))))))
PHIH(I) = 1. / PHIH(I)
ENDIF
9016 CONTINUE
500 CONTINUE
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
SUBROUTINE PSI(VZZ,VZH,VPSIM,VPSIH,IRUN,VX,VXS,VY,VYS,IFLAG)
C**********************************************************************
C SUBROUTINE PSI - DETERMINES DIMENSIONLESS WIND AND
C SCALAR PROFILES IN SURFACE LAYER
C - CALLED FROM SFCFLX
C DESCRIPTION OF PARAMETERS
C ZZ - INPUTED VALUE OF MONIN- OBUKHOV STABILITY PARAMETER ZETA
C ZH - INPUTED VALUE OF PBL HEIGHT DIVIDED BY Z0
C PSIM - OUTPUTED VALUE OF DIMENSIONLESS WIND
C PSIH - OUTPUTED VALUE OF DIMENSIONLESS SCALAR
C X - OUTPUTED VALUE OF PHIM(ZETA)
C XS - OUTPUTED VALUE OF PHIM(ZETA0)
C Y - OUTPUTED VALUE OF PHIH(ZETA)
C YS - OUTPUTED VALUE OF PHIH(ZETA0)
C IFLAG- FLAG TO DETERMINE IF CU IS NEEDED (IFLAG=2),
C IF CT IS NEEDED (IFLAG=3), OR BOTH (IFLAG=1)
C SUBPROGRAMS NEEDED
C PHI - COMPUTES SIMILARITY FUNCTION FOR MOMENTUM AND SCALARS
C**********************************************************************
implicit none
C Argument List Declarations
integer irun,iflag
_RL VZZ(IRUN),VZH(IRUN),VPSIM(IRUN),VPSIH(IRUN),
1 VX(IRUN),VXS(IRUN),VY(IRUN),VYS(IRUN)
C Local Variables
_RL ZWM,RZWM,Z0M,ZCM,RZCM,CM1,CM2,CM6,CM7,CM8ARG,YCM
PARAMETER ( ZWM = 1. )
PARAMETER ( RZWM = 1./ZWM )
PARAMETER ( Z0M = 0.2 )
PARAMETER ( ZCM = 42. )
PARAMETER ( RZCM = 1./ZCM )
PARAMETER ( CM1 = 1./126. )
PARAMETER ( CM2 = 1./(6.*CM1) )
PARAMETER ( CM6 = 6. / ( 1. + 6.*CM1 ) )
PARAMETER ( CM7 = CM2 + ZWM )
PARAMETER ( CM8ARG = CM7*ZCM*RZWM / (CM2+ZCM) )
PARAMETER ( YCM = 6. / ( 1. + 6.*CM1*ZCM ) )
integer INTSTB(irun),INTZ0(irun)
_RL ZZ0(irun),Z(irun),Z2(irun),Z1(irun),Z0(irun)
_RL X0(irun),X1(irun),Y0(irun),Y1(irun)
_RL PSI2(irun),TEMP(irun)
_RL HZ(irun),ARG0(irun),ARG1(irun),DX(irun)
_RL X0NUM(irun),X1NUM(irun),X0DEN(irun)
_RL X1DEN(irun),Y1DEN(irun),Z2ZWM(irun)
_RL cm3,cm4,cm5,cm8
integer ibit,indx
integer i
CM3 = sqrt( 0.2/CM1-0.01 )
CM4 = 1./CM3
CM5 = (10.-CM1) / (10.*CM1*CM3)
CM8 = 6. * LOG(CM8ARG)
DO 9000 I = 1,IRUN
VPSIM(I) = 0.
VPSIH(I) = 0.
VX(I) = 0.
VXS(I) = 0.
VY(I) = 0.
VYS(I) = 0.
ZZ0(I) = VZH(I)*VZZ(I)
9000 CONTINUE
IBIT = 0
DO 9122 I = 1,IRUN
IF(VZZ(I).LE.-1.e-7)IBIT = IBIT + 1
9122 CONTINUE
DO 9022 I = 1,IRUN
IF(VZZ(I).LE.-1.e-7)THEN
INTSTB(I) = 1
ELSE
INTSTB(I) = 0
ENDIF
9022 CONTINUE
C ****************************************
C ***** UNSTABLE SURFACE LAYER *****
C ****************************************
IF(IBIT.LE.0) GO TO 100
indx = 0
DO 9002 I = 1,IRUN
IF (INTSTB(I).EQ.1)THEN
indx = indx + 1
Z(indx) = VZZ(I)
Z0(indx) = ZZ0(I)
ENDIF
9002 CONTINUE
DO 9004 I = 1,IBIT
Z(I) = -18. * Z(I)
Z0(I) = -18. * Z0(I)
9004 CONTINUE
CALL PHI( Z,X1,Y1,IFLAG,IBIT )
CALL PHI( Z0,X0,Y0,IFLAG,IBIT )
C ****************************
C ***** COMPUTE PSIM *****
C ****************************
IF(IFLAG.GE.3) GO TO 75
DO 9006 I = 1,IBIT
ARG1(I) = 1. - X1(I)
IF ( Z(I) .LT. 0.013 ) ARG1(I) =
1 Z(I) * ( 0.25 - 0.09375 * Z(I) )
ARG0(I) = 1. - X0(I)
IF ( Z0(I) .LT. 0.013 ) ARG0(I) =
1 Z0(I) * ( 0.25 - 0.09375 * Z0(I) )
ARG1(I) = ARG1(I) * ( 1.+X0(I) )
ARG0(I) = ARG0(I) * ( 1.+X1(I) )
DX(I) = X1(I) - X0(I)
ARG1(I) = ARG1(I) / ARG0(I)
ARG0(I) = -DX(I) / ( 1. + X1(I)*X0(I) )
ARG0(I) = ATAN( ARG0(I) )
ARG1(I) = LOG( ARG1(I) )
PSI2(I) = 2. * ARG0(I) + ARG1(I)
PSI2(I) = PSI2(I) + DX(I)
9006 CONTINUE
indx = 0
DO 9008 I = 1,IRUN
IF( INTSTB(I).EQ.1 ) THEN
indx = indx + 1
VPSIM(I) = PSI2(indx)
VX(I) = X1(indx)
VXS(I) = X0(indx)
ENDIF
9008 CONTINUE
C ****************************
C ***** COMPUTE PSIH *****
C ****************************
IF(IFLAG.EQ.2) GO TO 100
75 CONTINUE
DO 9010 I = 1,IBIT
ARG1(I) = 1. - Y1(I)
IF( Z(I) .LT. 0.0065 ) ARG1(I) =
1 Z(I) * ( 0.5 - 0.625 * Z(I) )
ARG0(I) = 1. - Y0(I)
IF( Z0(I) .LT. 0.0065 ) ARG0(I) =
1 Z0(I) * ( 0.5 - 0.625 * Z0(I) )
ARG1(I) = ARG1(I) * ( 1. + Y0(I) )
ARG0(I) = ARG0(I) * ( 1. + Y1(I) )
ARG1(I) = ARG1(I) / ARG0(I)
PSI2(I) = LOG( ARG1(I) )
PSI2(I) = PSI2(I) - Y1(I) + Y0(I)
9010 CONTINUE
indx = 0
DO 9012 I = 1,IRUN
IF( INTSTB(I).EQ.1 ) THEN
indx = indx + 1
VPSIH(I) = PSI2(indx)
VY(I) = Y1(indx)
VYS(I) = Y0(indx)
ENDIF
9012 CONTINUE
C **************************************
C ***** STABLE SURFACE LAYER *****
C **************************************
100 CONTINUE
IBIT = 0
DO 9114 I = 1,IRUN
IF(VZZ(I).GT.-1.e-7)THEN
IBIT = IBIT + 1
ENDIF
9114 CONTINUE
DO 9014 I = 1,IRUN
IF(VZZ(I).GT.-1.e-7)THEN
INTSTB(I) = 1
ELSE
INTSTB(I) = 0
ENDIF
9014 CONTINUE
IF(IBIT.LE.0) GO TO 300
indx = 0
#ifdef CRAY
CDIR$ NOVECTOR
#endif
DO 9016 I = 1,IRUN
IF (INTSTB(I).EQ.1)THEN
indx = indx + 1
Z(indx) = VZZ(I)
Z0(indx) = ZZ0(I)
ARG1(indx) = VZH(I)
ENDIF
9016 CONTINUE
#ifdef CRAY
CDIR$ VECTOR
#endif
DO 9018 I = 1,IBIT
HZ(I) = 1. / ARG1(I)
Z1(I) = Z(I)
Z2(I) = ZWM
IF ( Z(I) .GT. ZWM ) THEN
Z1(I) = ZWM
Z2(I) = Z(I)
ENDIF
IF ( Z0(I) .GT. Z0M ) THEN
Z0(I) = Z0M
INTZ0(I) = 1
ELSE
INTZ0(I) = 0
ENDIF
X1NUM(I) = 1. + 5. * Z1(I)
X0NUM(I) = 1. + 5. * Z0(I)
X1DEN(I) = 1. / (1. + CM1 * (X1NUM(I) * Z1(I)) )
X0DEN(I) = 1. + CM1 * (X0NUM(I) * Z0(I))
IF ( (INTZ0(I).EQ.1) .OR. (Z(I).GT.ZWM) )
1 HZ(I) = Z1(I) / Z0(I)
ARG1(I) = HZ(I)*HZ(I)*X0DEN(I)*X1DEN(I)
ARG1(I) = LOG( ARG1(I) )
ARG1(I) = 0.5 * ARG1(I)
ARG0(I) = (Z1(I) + 0.1) * (Z0(I) + 0.1)
ARG0(I) = CM3 + ARG0(I) * CM4
ARG0(I) = ( Z1(I) - Z0(I) ) / ARG0(I)
ARG0(I) = ATAN( ARG0(I) )
TEMP(I) = ARG1(I) + CM5 * ARG0(I)
X0(I) = X0NUM(I) / X0DEN(I)
IF ( INTZ0(I).EQ.1 ) X0(I) = 0.
Z2ZWM(I) = Z2(I) * RZWM
9018 CONTINUE
C ****************************
C ***** COMPUTE PSIM *****
C ****************************
IF( IFLAG.GE.3 ) GO TO 225
DO 9020 I = 1,IBIT
X1(I) = X1NUM(I) * X1DEN(I)
ARG1(I) = LOG( Z2ZWM(I) )
PSI2(I) = TEMP(I) + CM6 * ARG1(I)
9020 CONTINUE
indx = 0
DO 9030 I = 1,IRUN
IF( INTSTB(I).EQ.1 ) THEN
indx = indx + 1
VPSIM(I) = PSI2(indx)
VX(I) = X1(indx)
VXS(I) = X0(indx)
ENDIF
9030 CONTINUE
C ****************************
C ***** COMPUTE PSIH *****
C ****************************
IF(IFLAG.EQ.2)GO TO 300
225 CONTINUE
DO 9024 I = 1,IBIT
Y1DEN(I) = 1. + CM1 * ( X1NUM(I) * Z(I) )
Y1(I) = X1NUM(I) / Y1DEN(I)
ARG1(I) = CM7 * Z2ZWM(I) / ( CM2 + Z2(I) )
ARG0(I) = 6.
IF ( Z2(I) .GT. ZCM ) THEN
Y1(I) = YCM
ARG1(I) = Z2(I) * RZCM
ARG0(I) = YCM
TEMP(I) = TEMP(I) + CM8
ENDIF
ARG1(I) = LOG( ARG1(I) )
PSI2(I) = TEMP(I) + ARG0(I) * ARG1(I)
9024 CONTINUE
indx = 0
DO 9026 I = 1,IRUN
IF( INTSTB(I).EQ.1 ) THEN
indx = indx + 1
VPSIH(I) = PSI2(indx)
VY(I) = Y1(indx)
VYS(I) = X0(indx)
ENDIF
9026 CONTINUE
300 CONTINUE
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
SUBROUTINE TRBLEN (STRT,DW2,DZ3,Q,VKZE,VKZM,DTHV,DPK,DU,DV,XL,
1 QXLM,NLEV,INIT,LMIN,LMINQ,LMINQ1,CP,INT1,INT2,
2 DZITRP,STBFCN,XL0,Q1,WRKIT1,WRKIT2,WRKIT3,WRKIT4,irun)
C**********************************************************************
C SUBROUTINE TRBLEN - COMPUTES TURBULENT LENGTH SCALE
C - CALLED FROM TRBFLX
C ARGUMENTS ::
C INPUT:
C ------
C STRT - BRUNT VAISALA FREQUENCY
C DW2 - SQUARED SHEAR
C DZ3 - LAYER THICKNESS FOR LENGTH SCALE CALC.
C Q - TURBULENCE VELOCITY
C VKZE - VK * Z AT LAYER EDGES
C VKZM - VK * Z AT LAYER CENTERS
C DTHV - VERTICAL GRADIENT OF THV
C DPK - VERTICAL GRADIENT OF PK
C DU - VERTICAL GRADIENT OF U
C DV - VERTICAL GRADIENT OF V
C NLEV - NUMBER OF ATMOSPHERIC LEVELS
C INIT - INPUT FLAG : 1 = INITIAL START
C 2 = 2ND CALL FOR INITIAL STAR
C 0 = NON-INITIAL START
C OUTPUT:
C -------
C XL - TURBULENT LENGTH SCALE
C QXLM - TURBULENT LENGTH SCALE * Q AT LAYER CENTER
C LMIN - HIGHEST LAYER WHERE INSTABILITY OCCURS
C LMINQ - HIGHEST LAYER WHERE TURBULENCE OCCURS
C
C SUBPROGRAMS NEEDED ::
C TRBITP - INTERPOLATES TO HEIGHT WHERE RI = RICR
C**********************************************************************
implicit none
C Argument List Declarations
integer irun,nlev,init,lmin,lminq,lminq1
_RL cp
_RL STRT(irun,NLEV),DW2(irun,NLEV),DZ3(irun,NLEV)
_RL Q(irun,NLEV),VKZM(irun,NLEV-1),VKZE(irun,NLEV-1)
_RL DTHV(irun,NLEV),DPK(irun,NLEV),DU(irun,NLEV)
_RL DV(irun,NLEV)
_RL QXLM(irun,NLEV-1),XL(irun,NLEV-1)
_RL DZITRP(irun,nlev-1),STBFCN(irun,nlev)
_RL XL0(irun,nlev),Q1(irun,nlev-1)
_RL WRKIT1(irun,nlev-1),WRKIT2(irun,nlev-1)
_RL WRKIT3(irun,nlev-1)
_RL WRKIT4(irun,nlev-1)
INTEGER INT1(irun,nlev), INT2(irun,nlev-1)
C Local Variables
_RL rf1,rf2,e5,d4,d1,rfc,ricr,alpha,dzcnv,xl0cnv,xl0min
_RL clmt,clmt53
PARAMETER ( RF1 = 0.2340678 )
PARAMETER ( RF2 = 0.2231172 )
PARAMETER ( E5 = 49.66 )
PARAMETER ( D4 = 2.6532122E-2 )
PARAMETER ( D1 = D4 * E5 )
PARAMETER ( RFC = 0.1912323 )
PARAMETER ( RICR = ( (RF1-RFC)*RFC ) / ( (RF2-RFC)*D1 ) )
PARAMETER ( ALPHA = 0.1 )
PARAMETER ( DZCNV = 100. )
PARAMETER ( XL0CNV = DZCNV * ALPHA )
PARAMETER ( XL0MIN = 1. )
PARAMETER ( CLMT = 0.23 )
PARAMETER ( CLMT53 = 5. * CLMT / 3. )
integer ibit,nlevm1,nlevp1,istnlv,istnm1,nlevml,istnml,Lp1
integer istnmq,istlmq,lminp,lm1,lmin1
integer i,L,LL
NLEVM1 = NLEV - 1
NLEVP1 = NLEV + 1
ISTNLV = irun * NLEV
ISTNM1 = irun * NLEVM1
IF ( INIT.EQ.2 ) GO TO 1200
C COMPUTE DEPTHS OF UNSTABLE LAYERS
C =================================
DO L =1,NLEV
DO I =1,irun
STBFCN(I,L) = STRT(I,L) - RICR * DW2(I,L)
INT1(I,L) = 0
IF( STBFCN(I,L).LE.0. ) INT1(I,L) = 1
ENDDO
ENDDO
DO L =1,NLEVM1
DO I =1,irun
INT2(I,L) = 0
IF( (INT1(I,L).EQ.1) .NEQV. (INT1(I,L+1).EQ.1) ) INT2(I,L) = 1
ENDDO
ENDDO
DO 40 LMIN = 1,NLEV
IBIT = 0
DO 30 I=1,irun
IBIT = IBIT + INT1(I,LMIN)
30 CONTINUE
IF(IBIT.GE.1) GO TO 50
40 CONTINUE
LMIN = NLEVP1
50 CONTINUE
LMIN = 1
DO L =1,NLEVM1
DO I =1,irun
XL0(I,L) = 0.
ENDDO
ENDDO
DO 70 I=1,irun
XL0(I,NLEV) = DZ3(I,NLEV)
70 CONTINUE
IF(LMIN.GE.NLEVP1) GOTO 1100
LMIN1 = LMIN - 1
IF(LMIN1.EQ.0) LMIN1 = 1
NLEVML = NLEV - LMIN1
ISTNML = irun*NLEVML
CALL TRBITP ( STBFCN(1,LMIN1),INT2(1,LMIN1),DTHV(1,LMIN1),
& DPK(1,LMIN1), DU(1,LMIN1), DV(1,LMIN1),
& DZITRP(1,LMIN1), NLEVML,
& WRKIT1,WRKIT2,WRKIT3,WRKIT4,CP,irun )
LP1 = LMIN1 + 1
DO L =LMIN1,NLEVM1
DO I =1,irun
INT2(I,L) = 0
IF ( INT1(I,L).EQ.1 .OR. INT1(I,L+1).EQ.1 ) INT2(I,L) = 1
IF ( INT2(I,L).EQ.1 )
& XL0(I,L) = (0.5+DZITRP(I,L)) * DZ3(I,L+1)
ENDDO
ENDDO
DO 90 I=1,irun
INT2(I,NLEVM1) = INT1(I,NLEV)
90 CONTINUE
DO L =LMIN1,NLEVM1
DO I =1,irun
IF ( INT2(I,L).EQ.1 ) THEN
XL0(I,L+1) = XL0(I,L+1) + ( (0.5-DZITRP(I,L)) * DZ3(I,L+1) )
ENDIF
ENDDO
ENDDO
IF (LMIN.GT.1) GOTO 400
DO 110 I=1,irun
IF( INT1(I,1).EQ.1 ) XL0(I,1) = XL0(I,1) + DZ3(I,1)
110 CONTINUE
400 CONTINUE
LMINP = LMIN + 1
IF(LMINP.GT.NLEVM1) GOTO 550
DO 500 L = LMINP,NLEVM1
LM1 = L-1
DO 120 I = 1,irun
IF( INT1(I,LM1).EQ.1 ) XL0(I,L) = XL0(I,L) + XL0(I,LM1)
120 CONTINUE
500 CONTINUE
550 CONTINUE
IF(LMIN.GT.NLEVM1) GOTO 600
DO 130 I = 1,irun
IF( INT1(I,NLEVM1).EQ.1 .AND. INT1(I,NLEV).EQ.1 ) THEN
XL0(I,NLEV) = XL0(I,NLEV) + XL0(I,NLEVM1)
ENDIF
130 CONTINUE
IF(LMIN.GT.NLEV) GOTO 1100
600 CONTINUE
DO 1000 LL = LMIN,NLEV-1
L = NLEVM1 + LMIN - LL
LP1 = L+1
DO 140 I = 1,irun
IF( INT1(I,LP1).EQ.1 ) THEN
IF( INT1(I,L) .EQ.1 ) THEN
XL0(I,L) = XL0(I,LP1)
ELSE
XL0(I,L) = XL0(I,L) + XL0(I,LP1)
ENDIF
ENDIF
140 CONTINUE
1000 CONTINUE
1100 CONTINUE
DO L =1,NLEV
DO I =1,irun
IF ( XL0(I,L).LT.XL0CNV ) XL0(I,L) = XL0CNV
ENDDO
ENDDO
C *********************************************************************
C **** DETERMINE MIXING LENGTHS FOR STABLE LAYERS ***
C *********************************************************************
IF(INIT.EQ.1) GOTO 1400
IF(LMINQ.GT.1) THEN
ISTLMQ = irun * LMINQ1
DO L =1,LMINQ1
DO I =1,irun
INT2(I,L) = 1 - INT1(I,L)
ENDDO
ENDDO
ENDIF
IF(LMINQ.LT.NLEV) THEN
ISTNMQ = irun * (NLEV-LMINQ)
DO L =LMINQ,NLEVM1
DO I =1,irun
IF( INT1(I,L).EQ.0 ) THEN
XL0(I,L) = Q(I,L) / XL0(I,L)
XL0(I,L) = XL0(I,L) * XL0(I,L) + 1.0E-20
XL0(I,L) = STBFCN(I,L) + XL0(I,L)
XL0(I,L) = SQRT( XL0(I,L) )
XL0(I,L) = Q(I,L) / XL0(I,L)
ENDIF
INT2(I,L) = 0
IF( XL0(I,L).LT.XL0MIN ) INT2(I,L) = 1
ENDDO
ENDDO
ENDIF
1200 CONTINUE
IF(INIT.EQ.2) THEN
DO L =1,NLEVM1
DO I =1,irun
INT2(I,L) = 1 - INT1(I,L)
IF ( INT2(I,L).EQ.1 ) XL0(I,L) = XL0MIN
ENDDO
ENDDO
ENDIF
DO L =1,NLEVM1
DO I =1,irun
IF ( INT2(I,L).EQ.1 ) XL0(I,L) = XL0MIN
ENDDO
ENDDO
C *********************************************************************
C **** LENGTH SCALE XL FROM XL0 AND VKZE ****
C *********************************************************************
1400 CONTINUE
DO L =1,NLEVM1
DO I =1,irun
XL(I,L) = XL0(I,L) * VKZE(I,L) / ( XL0(I,L)+VKZE(I,L) )
ENDDO
ENDDO
C *********************************************************************
C **** CLMT53 TIMES Q TIMES LENGTH SCALE AT MID LEVELS ***
C *********************************************************************
IF(INIT.EQ.1) GOTO 1700
ISTNMQ = irun * (NLEV-LMINQ1)
DO L =LMINQ1,NLEVM1
DO I =1,irun
Q1(I,L) = Q(I,L)
INT1(I,L) = 0
IF ( Q(I,L).LE.Q(I,L+1) ) INT1(I,L) = 1
IF ( INT1(I,L).EQ.1 ) THEN
XL0(I,L) = XL0(I,L+1)
Q1(I,L) = Q(I,L+1)
ENDIF
ENDDO
ENDDO
DO L =LMINQ1,NLEVM1
DO I =1,irun
QXLM(I,L) = XL0(I,L)*VKZM(I,L)
& / ( XL0(I,L)+VKZM(I,L) )
QXLM(I,L) = CLMT53 * Q1(I,L)*QXLM(I,L)
ENDDO
ENDDO
1700 CONTINUE
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
SUBROUTINE TRBITP ( STBFCN,INTCHG,DTHV,DPK,DU,DV,DZITRP,NLEV,
& AAA,BBB,CCC,DDD,CP,irun )
C**********************************************************************
C SUBROUTINE TRBITP - INTERPOLATES TO THE HEIGHT WHERE RI EQUALS RICR
C - CALLED FROM TRBLEN
C ARGUMENTS ::
C INPUT:
C ------
C STBFCN - DTHV * DPK - RICR*( DU*DU + DV*DV)
C INTCHG - INT '1' AT LEVELS WHERE STBFCN CHANGES SIG
C DTHV - VERTICAL GRADIENT OF THV
C DPK - VERTICAL GRADIENT OF PK
C DU - VERTICAL GRADIENT OF U
C DV - VERTICAL GRADIENT OF V
C NLEV - NUMBER OF LEVELS TO BE PROCESSED
C OUTPUT:
C -------
C DZITRP - INTERPOLATION COEFFICIENT
C**********************************************************************
implicit none
C Argument List Declarations
integer irun,nlev
_RL cp
_RL STBFCN(irun,NLEV+1)
integer INTCHG(irun,NLEV)
_RL DTHV(irun,NLEV+1),DPK(irun,NLEV+1)
_RL DU(irun,NLEV+1),DV(irun,NLEV+1)
_RL DZITRP(irun,NLEV+1)
_RL AAA(irun,NLEV),BBB(irun,NLEV)
_RL CCC(irun,NLEV),DDD(irun,NLEV)
C Local Variables
_RL rf1,rf2,e5,d4,d1,rfc,ricr
PARAMETER ( RF1 = 0.2340678 )
PARAMETER ( RF2 = 0.2231172 )
PARAMETER ( E5 = 49.66 )
PARAMETER ( D4 = 2.6532122E-2 )
PARAMETER ( D1 = D4 * E5 )
PARAMETER ( RFC = 0.1912323 )
PARAMETER ( RICR = ( (RF1-RFC)*RFC ) / ( (RF2-RFC)*D1 ) )
integer istnlv
integer I,L
C *********************************************************************
C **** QUADRATIC INTERPOLATION OF RI TO RICR VIA ***
C **** LINEAR INTERPOLATION OF DTHV, DPK, DU & DV ***
C *********************************************************************
ISTNLV = irun*NLEV
DO L =1,NLEV
DO I =1,irun
DZITRP(I,L) = 0.
ENDDO
ENDDO
DO L =1,NLEV
DO I =1,irun
IF ( INTCHG(I,L).EQ.1 ) THEN
DDD(I,L) = ( CP *(DTHV(I,L+1)*DPK(I,L)
& + DTHV(I,L)*DPK(I,L+1)) )
& - ( (2.*RICR) * ( DU(I,L+1)* DU(I,L)
& + DV(I,L+1)* DV(I,L)) )
AAA(I,L) = STBFCN(I,L) + STBFCN(I,L+1)
BBB(I,L) = STBFCN(I,L) - STBFCN(I,L+1)
CCC(I,L) = 1. / BBB(I,L)
DZITRP(I,L) = AAA(I,L) * CCC(I,L)
AAA(I,L) = AAA(I,L) - DDD(I,L)
DDD(I,L) = ( DDD(I,L) * DDD(I,L) )
& - 4. * (STBFCN(I,L+1) * STBFCN(I,L) )
DDD(I,L) = DDD(I,L)*CCC(I,L)*CCC(I,L)
DDD(I,L) = SQRT( DDD(I,L) )
ENDIF
IF( INTCHG(I,L).EQ.1 .AND. AAA(I,L).NE.0. ) THEN
DZITRP(I,L) = ( BBB(I,L)*(1.-DDD(I,L)) ) / AAA(I,L)
ENDIF
DZITRP(I,L) = 0.5 * DZITRP(I,L)
ENDDO
ENDDO
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
SUBROUTINE TRBL20 (RI,STRT,DW2,XL,ZKM,ZKH,QE,QQE,INTSTB,NLEV,
1 nlay,irun)
C**********************************************************************
C SUBROUTINE TRBL20 - COMPUTES QE AND DIMLESS COEFS FROM
C MELLOR-YAMADA LEVEL 2 MODEL
C - CALLED FROM AND FROM TRBFLX
C ARGUMENTS ::
C INPUT:
C ------
C RI - RICHARDSON NUMBER
C STRT - BRUNT VAISALA FREQUENCY
C DW2 - SQUARED SHEAR
C XL - TURBULENT LENGTH SCALE
C NLEV - NUMBER OF LEVELS TO BE PROCESSED
C OUTPUT:
C -------
C ZKM - MOMENTUM TRANSPORT COEFFICIENT
C ZKH - HEAT TRANSPORT COEFFICIENT
C QE - EQUILIBRIUM TURBULENT VELOCITY SCALE
C QQE - EQUILIBRIUM TURBULENT KINETIC ENERGY
C BITSTB - BIT '1' WHERE QE GREATER THAN ZERO
C**********************************************************************
implicit none
C Argument List Declarations
integer nlev,nlay,irun
_RL RI(irun*NLEV,1),STRT(irun*NLEV,1),DW2(irun*NLEV,1)
_RL XL(irun*NLEV,1),ZKM(irun*NLEV,1), ZKH(irun*NLEV,1)
_RL QE(irun*NLEV,1),QQE(irun*NLEV,1)
INTEGER INTSTB(irun*NLEV,1)
_RL EE(irun*(nlay-1),1),RF(irun*(nlay-1),1)
C Local Variables
_RL b1,b2,d3,rf1,rf2,d3b2,d2,e5,d4,d1,d1half,d2half
_RL rfc,ricr,ch,cm
PARAMETER ( B1 = 16.6 )
PARAMETER ( B2 = 10.1 )
PARAMETER ( D3 = 0.29397643 )
PARAMETER ( RF1 = 0.2340678 )
PARAMETER ( RF2 = 0.2231172 )
PARAMETER ( D3B2 = D3 / RF1 )
PARAMETER ( D2 = RF1 )
PARAMETER ( E5 = 49.66 )
PARAMETER ( D4 = 2.6532122E-2 )
PARAMETER ( D1 = D4 * E5 )
PARAMETER ( D1HALF = 0.5 * D1 )
PARAMETER ( D2HALF = 0.5 * D2 )
PARAMETER ( RFC = 0.1912323 )
PARAMETER ( RICR = ( (RF1-RFC)*RFC ) / ( (RF2-RFC)*D1 ) )
PARAMETER ( CH = 2.5828674 )
PARAMETER ( CM = CH / D1 )
integer istnlv
integer i
ISTNLV = irun * NLEV
C *********************************************************************
C **** COMPUTE FLUX RICHARDSON NUMBER ***
C *********************************************************************
DO 10 I=1,ISTNLV
EE(I,1) = D1HALF * RI(I,1) + D2HALF
RF(I,1) = EE(I,1)* EE(I,1)
RF(I,1) = RF(I,1)- D3*RI(I,1)
RF(I,1) = SQRT( RF(I,1) )
RF(I,1) = EE(I,1) - RF(I,1)
IF( RI(I,1).LE.1.e-4 .AND. RI(I,1).GE.-1.e-4 ) THEN
RF(I,1) = D3B2*RI(I,1)
ENDIF
C *********************************************************************
C **** QE AND DIMENSIONLESS DIFFUSION COEFICIENTS ***
C **** FROM LEVEL 2 CLOSURE MODEL ***
C *********************************************************************
IF( RI(I,1).LT.RICR .AND. RF(I,1).LT.RFC ) THEN
ZKH(I,1) = ( RFC-RF(I,1) ) / (1.-RF(I,1))
ZKM(I,1) = CM * (RF1-RF(I,1))
ZKM(I,1) = ZKH(I,1)*ZKM(I,1) / (RF2-RF(I,1))
ZKH(I,1) = CH *ZKH(I,1)
QE(I,1) = ZKM(I,1)*DW2(I,1) - ZKH(I,1)*STRT(I,1)
ENDIF
IF( QE(I,1).LT.1.e-14 ) THEN
INTSTB(I,1) = 0
QE(I,1) = 0.
ELSE
INTSTB(I,1) = 1
QE(I,1) = B1*QE(I,1)
QE(I,1) = SQRT( QE(I,1) )
QE(I,1) = XL(I,1)*QE(I,1)
ENDIF
QQE(I,1) = 0.5 * QE(I,1) * QE(I,1)
10 CONTINUE
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
SUBROUTINE TRBL25(Q,XL,STRT,DW2,INTSTB,INTQ,ZKM,ZKH,P3,NLEV,
1 nlay,irun)
C**********************************************************************
C SUBROUTINE TRBL25 - COMPUTES P3 AND DIMLESS COEFS FROM
C MELLOR-YAMADA LEVEL 2.5 MODEL
C - CALLED FROM TRBFLX
C ARGUMENTS ::
C INPUT:
C ------
C Q - TURBULENCE VELOCITY
C XL - TURBULENT LENGTH SCALE
C STRT - BRUNT VAISALA FREQUENCY
C DW2 - SQUARED SHEAR
C BITSTB - BIT '1' WHERE QE GREATER THAN ZERO
C NLEV - NUMBER OF LEVELS TO BE PROCESSED
C OUTPUT:
C -------
C ZKM - MOMENTUM TRANSPORT COEFFICIENT
C ZKH - HEAT TRANSPORT COEFFICIENT
C P3 - PRODUCTION RATE OF TURBULENT KINETIC ENERG
C**********************************************************************
implicit none
C Argument list Declarations
integer nlev,nlay,irun
_RL Q(irun*NLEV,1),XL(irun*NLEV,1),STRT(irun*NLEV,1)
_RL DW2(irun*NLEV,1)
INTEGER INTSTB(irun*nlay,1), INTQ(irun*nlay,1)
_RL ZKM(irun*NLEV,1),ZKH(irun*NLEV,1),P3(irun*NLEV,1)
C Local Variables
_RL a1,a2,a4,c1,a5,a3,b1,b2,b3,ff2,ff3,ff4
PARAMETER ( A1 = 0.92 )
PARAMETER ( A2 = 0.74 )
PARAMETER ( A4 = 6. * A1 * A1)
PARAMETER ( C1 = 0.08 )
PARAMETER ( A5 = 3.*C1*(-1.) )
PARAMETER ( A3 = A4 * A5*(-1.) )
PARAMETER ( B1 = 16.6 )
PARAMETER ( B2 = 10.1 )
PARAMETER ( B3 = 1. / B1 )
PARAMETER ( FF2 = 9. * A1 * A2 )
PARAMETER ( FF3 = (3.*A2*B2) - (9.*A2*A2 ) )
PARAMETER ( FF4 = (3.*A2*B2) + (12.*A1*A2 ) )
_RL F2(irun*(nlay-1),1),F3(irun*(nlay-1),1)
_RL F4(irun*(nlay-1),1),XQ(irun*(nlay-1),1)
integer istnlv
integer i
ISTNLV = irun * NLEV
C *********************************************************************
C **** P3 AND DIMENSIONLESS DIFFUSION COEFICIENTS ***
C **** FROM LEVEL 2.5 CLOSURE MODEL ***
C *********************************************************************
DO 10 I=1,ISTNLV
IF( INTQ(I,1).EQ.1 .AND. INTSTB(I,1).EQ.0 ) THEN
XQ(I,1) = XL(I,1) / Q(I,1)
XQ(I,1) = XQ(I,1) * XQ(I,1)
STRT(I,1) = XQ(I,1) * STRT(I,1)
DW2(I,1) = XQ(I,1) * DW2(I,1)
F2(I,1) = 1.+FF2 * STRT(I,1)
F3(I,1) = 1.+FF3 * STRT(I,1)
F4(I,1) = 1.+FF4 * STRT(I,1)
ZKH(I,1) = (F4(I,1) * F2(I,1))
& + A4 * (F3(I,1) * DW2(I,1))
ZKH(I,1) = (F2(I,1) + A3*DW2(I,1))
& / ZKH(I,1)
ZKM(I,1) = A1 * (F3(I,1)*ZKH(I,1)+A5)
& / F2(I,1)
ZKH(I,1) = A2 * ZKH(I,1)
P3(I,1) = ZKH(I,1)*STRT(I,1) + B3
P3(I,1) = 2. * ( ZKM(I,1)*DW2(I,1) - P3(I,1) )
P3(I,1) = P3(I,1)*Q(I,1)
ENDIF
10 CONTINUE
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
SUBROUTINE TRBDIF ( XX1,XX2,RHOKDZ,FLXFAC,DXX1G,DXX2G,NLEV,
& ITYPE,EPSL,irun )
C**********************************************************************
C ARGUMENTS ::
C INPUT:
C ------
C XX1 - FIRST PROPERTY TO BE DIFFUSED
C (INPUT INCLUDES FORWARD PRODUCTION TERM)
C XX2 - SECOND PROPERTY TO BE DIFFUSED (V-WIND)
C (INPUT INCLUDES FORWARD PRODUCTION TERM)
C -OR-
C CHANGE IN XX1 DUE TO UNIT CHANGE IN THG
C (TH OR SH PROFILES)
C -OR-
C BACKWARD PRODUCTION TERM (QQ)
C RHOKDZ - RHO * K * WEIGHT / DZ AT INTERFACES
C FLXFAC - G * DT / (DP*WEIGHT) AT EDGES
C NLEV - NUMBER OF ATMOSPHERIC LEVELS
C ITYPE - INTEGER FLAG FOR INPUT TYPE
C 1 = QQ: COMPUTE BACKWARD PRODUCTION AND
C USE UNDERFLOW CUTOFF
C 2 = TH OR SH: COMPUTE TENDENCY DUE TO
C SURFACE PERTURBATION
C 3 = U AND V: COMPUTE BOTH FIELDS
C EPSL - UNDERFLOW CUTOFF CRITERION (QQ ONLY)
C OUTPUT:
C ------
C XX1 - NEW VALUE RETURNED
C XX2 - NEW VALUE RETURNED
C DXX1G - SOURCE TERM FOR XX1 AT GROUND
C DXX1G - SOURCE TERM FOR XX2 AT GROUND
C**********************************************************************
implicit none
C Argument List Declarations
integer nlev,itype,irun
_RL XX1(irun,NLEV+1),XX2(irun,NLEV+1)
_RL RHOKDZ(irun,NLEV),FLXFAC(irun,NLEV+1)
_RL DXX1G(irun),DXX2G(irun)
_RL epsl
_RL AA(irun,nlev), BB(irun,nlev), CC(irun,nlev+1)
integer NLEVP1,NLEVX
integer I,L
NLEVP1 = NLEV + 1
NLEVX = NLEV - 1
IF(ITYPE.EQ.2) NLEVX = NLEV
C DEFINE MATRIX
DO 10 I=1,irun
CC(I,1) = 0.
10 CONTINUE
DO L =1,NLEVX
DO I =1,irun
CC(I,L+1) = RHOKDZ(I,L) * FLXFAC(I,L+1)
ENDDO
ENDDO
DO L =1,NLEV
DO I =1,irun
BB(I,L) = RHOKDZ(I,L) * FLXFAC(I,L)
AA(I,L) = 1. + CC(I,L) + BB(I,L)
ENDDO
ENDDO
C ADD IMPLICIT BACKWARD FORCING FOR QQ
IF(ITYPE.EQ.1) THEN
DO L =1,NLEV
DO I =1,irun
AA(I,L) = AA(I,L) - XX2(I,L)
ENDDO
ENDDO
ENDIF
C SOLVE MATRIX EQUATION FOR XX1
CALL VTRI0(AA,BB,CC,XX1,XX1,NLEV,irun)
IF(ITYPE.EQ.2) THEN
C COMPUTE CHANGE AT SURFACE
DO 50 I=1,irun
DXX1G(I) = CC(I,NLEVP1) * ( XX1(I,NLEV)-XX1(I,NLEVP1) )
50 CONTINUE
C SOLVE MATRIX FOR SURFACE PERTURBATION
CALL VTRI1(AA,BB,XX2,NLEV,irun)
DO 60 I=1,irun
DXX2G(I) = CC(I,NLEVP1) * ( XX2(I,NLEV)-XX2(I,NLEVP1) )
60 CONTINUE
ENDIF
C SOLVE MATRIX EQUATION FOR XX2
IF(ITYPE.EQ.3) CALL VTRI2 (AA,BB,CC,XX2,XX2,NLEV,irun)
C ELIMINATE UNDERFLOW
IF(ITYPE.EQ.1) THEN
DO L =1,NLEV
DO I =1,irun
IF ( XX1(I,L).LT.EPSL ) XX1(I,L) = 0.
ENDDO
ENDDO
ENDIF
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
SUBROUTINE VTRI0 ( A,B,C,F,Y,K,irun)
implicit none
integer k,irun
_RL A(irun,K),B(irun,K),C(irun,K),Y(irun,K+1)
_RL F(irun,K)
integer i,L,Lm1
DO 9000 I = 1,irun
A(I,1) = 1. / A(I,1)
9000 CONTINUE
DO 100 L = 2,K
LM1 = L - 1
DO 9002 I = 1,irun
C(I,L) = C(I,L) * A(I,LM1)
A(I,L) = 1. / ( A(I,L) - B(I,LM1) * C(I,L) )
F(I,L) = F(I,L) + F(I,LM1) * C(I,L)
9002 CONTINUE
100 CONTINUE
DO 200 L = K,1,-1
DO 9004 I = 1,irun
Y(I,L) = (F(I,L) + B(I,L) * Y(I,L+1)) * A(I,L)
9004 CONTINUE
200 CONTINUE
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
SUBROUTINE VTRI1 ( A,B,Y,K,irun)
implicit none
integer k,irun
_RL A(irun,K),B(irun,K),Y(irun,K+1)
integer i,L
DO 200 L = K,1,-1
DO 9000 I = 1,irun
Y(I,L) = B(I,L) * Y(I,L+1) * A(I,L)
9000 CONTINUE
200 CONTINUE
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
SUBROUTINE VTRI2 ( A,B,C,F,Y,K,irun)
implicit none
integer k,irun
_RL A(irun,K),B(irun,K),C(irun,K),F(irun,K)
_RL Y(irun,K+1)
integer i,L
DO 100 L = 2,K
DO 9000 I = 1,irun
F(I,L) = F(I,L) + F(I,L-1) * C(I,L)
9000 CONTINUE
100 CONTINUE
DO 200 L = K,1,-1
DO 9002 I = 1,irun
Y(I,L) = (F(I,L) + B(I,L) * Y(I,L+1)) * A(I,L)
9002 CONTINUE
200 CONTINUE
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
SUBROUTINE LINADJ ( NN,VRIB1,VRIB2,VWS1,VWS2,VZ1,VUSTAR,IWATER,
1 VAPSIM, VAPSIHG,VPSIH,VPSIG,VX,VX0,VY,VY0,ITYPE,LWATER,IRUN,
2 VDZETA,VDZ0,VDPSIM,VDPSIH,INTRIB,
3 VX0PSIM,VG,VG0,VR1MG0,VZ2,VDZSEA,VAZ0,VXNUM1,VPSIGB2,VDX,
4 VDXPSIM,VDY,VXNUM2,VDEN,VAWS1,VXNUM3,VXNUM,VDZETA1,VDZETA2,
5 VZCOEF2,VZCOEF1,VTEMPLIN,VDPSIMC,VDPSIHC)
C**********************************************************************
C ARGUMENTS ::
C INPUT:
C ------
C RIB1 - BULK RICHARDSON NUMBER OF INPUT STATE
C RIB2 - DESIRED BULK RICH NUMBER OF OUTPUT STATE
C WS1 - SURFACE WIND SPEED OF INPUT STATE
C WS2 - DESIRED SURFACE WIND SPEED OF OUTPUT STATE
C Z1 - INPUT VALUE OF ROUGHNESS HEIGHT
C USTAR - INPUT VALUE OF CU * WS
C WATER - BIT ARRAY - '1' WHERE OCEAN
C APSIM - (1/PSIM)
C APSIHG - ( 1 / (PSIH+PSIG) )
C PSIH - NON-DIM TEMP GRADIENT
C PSIG - PSIH FOR THE MOLECULAR LAYER
C X - PHIM(ZETA) - DERIVATIVE OF PSIM
C X0 - PHIM(ZETA0)
C Y - PHIH(ZETA) - DERIVATIVE OF PSIH
C Y0 - PHIH(ZETA0)
C ITYPE - INTEGER FLAG :
C 1 = NEUTRAL ADJUSTMENT
C 2 = ADJ FOR 2ND OR GREATER TRBFLX ITER
C 3 - 5 = ADJUSTMENT INSIDE LOOP
C 4 - 5 = ADJUST CU AND CT
C 5 = PREPARATION FOR ITYPE = 2
C LWATER - LOGICAL - .TRUE. IF THERE ARE WATER POINTS
C OUTPUT:
C -------
C DZETA - D LOG ZETA
C DZ0 - D Z0 (ITYPE 1) OR D LOG Z0 (ITYPE 2-5)
C DPSIM - D PSIM
C DPSIH - D PSIH
C BITRIB - BIT ARRAY - '1' WHERE RIB1 = 0
C**********************************************************************
implicit none
C Argument List Declarations
integer nn,irun,itype
_RL VRIB1(IRUN),VRIB2(IRUN)
_RL VWS1(IRUN),VWS2(IRUN),VZ1(IRUN),VUSTAR(IRUN)
integer IWATER(IRUN)
_RL VAPSIM(IRUN),VAPSIHG(IRUN)
_RL VPSIH(IRUN),VPSIG(IRUN),VX(IRUN)
_RL VX0(IRUN),VY(IRUN),VY0(IRUN)
LOGICAL LWATER
_RL VDZETA(IRUN),VDZ0(IRUN),VDPSIM(IRUN)
_RL VDPSIH(IRUN)
integer INTRIB(IRUN)
_RL VX0PSIM(irun),VG(irun),VG0(irun),VR1MG0(irun)
_RL VZ2(irun),VDZSEA(irun),VAZ0(irun),VXNUM1(irun)
_RL VPSIGB2(irun),VDX(irun),VDXPSIM(irun),VDY(irun)
_RL VXNUM2(irun),VDEN(irun),VAWS1(irun),VXNUM3(irun)
_RL VXNUM(irun),VDZETA1(irun),VDZETA2(irun)
_RL VZCOEF2(irun),VZCOEF1(irun),VTEMPLIN(irun)
_RL VDPSIMC(irun),VDPSIHC(irun)
C Local Variables
_RL xx0max,prfac,xpfac,difsqt,ustz0s,h0byz0,usth0s
PARAMETER ( XX0MAX = 1.49821 )
PARAMETER ( PRFAC = 0.595864 )
PARAMETER ( XPFAC = .55 )
PARAMETER ( DIFSQT = 3.872983E-3)
PARAMETER ( USTZ0S = 0.2030325E-5)
PARAMETER ( H0BYZ0 = 30.0 )
PARAMETER ( USTH0S = H0BYZ0*USTZ0S )
integer VINT1(irun),VINT2(irun)
_RL getcon,vk,bmdl,b2uhs
integer i
vk = getcon('VON KARMAN')
BMDL = VK * XPFAC * PRFAC / DIFSQT
B2UHS = BMDL * BMDL * USTH0S
C COMPUTE X0/PSIM, 1/Z0, G, G0, 1/(1-G0),
C DEL LOG Z0, D LOG ZO / D USTAR
CCCOOOMMMM ADDED 'WHERE WATER'
IF ( (ITYPE.EQ.1) .AND. LWATER ) THEN
DO 9000 I = 1,IRUN
IF (IWATER(I).EQ.1) VX0PSIM(I) = VAPSIM(I)
9000 CONTINUE
ENDIF
IF ( ITYPE .GE. 3 ) THEN
DO 9002 I = 1,IRUN
VX0PSIM(I) = VX0(I) * VAPSIM(I)
9002 CONTINUE
ENDIF
IF ( ITYPE .NE. 2 ) THEN
DO 9004 I = 1,IRUN
VDZ0(I) = 0.
VG(I) = 0.
VG0(I) = 0.
VR1MG0(I) = 1.
9004 CONTINUE
IF ( LWATER ) THEN
CALL ZCSUB ( VUSTAR,VDZSEA,IWATER,.TRUE.,IRUN,VZ2)
DO 9006 I = 1,IRUN
IF ( IWATER(I).EQ.1) THEN
VAZ0(I) = 1. / VZ1(I)
VG(I) = VDZSEA(I) * VAZ0(I)
VG0(I) = VX0PSIM(I) * VG(I)
VR1MG0(I) = 1. / ( 1. - VG0(I) )
VDZ0(I) = ( VZ2(I) - VZ1(I) ) * VR1MG0(I)
ENDIF
9006 CONTINUE
ENDIF
ENDIF
IF ( LWATER .AND. (ITYPE.GE.3) ) THEN
DO 9008 I = 1,IRUN
IF (IWATER(I).EQ.1) VDZ0(I) = VDZ0(I) * VAZ0(I)
9008 CONTINUE
ENDIF
C COMPUTE NUM1,NUM2,NUM3, DEN
IF (ITYPE.GE.3) THEN
DO 9010 I = 1,IRUN
VXNUM1(I) = 0.
IF (VRIB1(I).EQ.0.) THEN
INTRIB(I) = 1
ELSE
INTRIB(I) = 0
ENDIF
IF ( INTRIB(I).EQ.0 ) VXNUM1(I) = 1. / VRIB1(I)
VPSIGB2(I) = 0.
if(vpsig(i).gt.0.)VPSIGB2(I) =
1 0.5 * ( vpsig(i)*vpsig(i) + b2uhs ) / vpsig(i)
VDX(I) = VX(I) - VX0(I)
VDXPSIM(I) = VDX(I) * VAPSIM(I)
VDY(I) = VY(I) - VY0(I)
VXNUM3(I) = - VPSIGB2(I)
IF ( LWATER ) THEN
CCCOOOMMMM ADDED 'WHERE WATER'
IF (IWATER(I).EQ.1) THEN
VDXPSIM(I) = VDXPSIM(I) * VR1MG0(I)
VXNUM3(I) = VXNUM3(I) + VG(I) * ( VY0(I) - VPSIGB2(I) )
VXNUM2(I) = VY0(I) - VPSIGB2(I) - VX0PSIM(I) * VPSIGB2(I)
VXNUM2(I) = (VXNUM2(I) * VAPSIHG(I)) - 2. * VX0PSIM(I)
VXNUM2(I) = VXNUM2(I) * VDZ0(I)
ENDIF
ENDIF
VDEN(I) = VDY(I) + VDXPSIM(I) * VXNUM3(I)
VDEN(I) = ( 1. + VDEN(I) * VAPSIHG(I) ) - 2. * VDXPSIM(I)
9010 CONTINUE
ENDIF
IF (ITYPE.EQ.5) THEN
DO 9012 I = 1,IRUN
VAWS1(I) = VR1MG0(I) / VWS1(I)
VXNUM3(I) = VXNUM3(I) * VAPSIHG(I)
IF ( LWATER ) THEN
CCCOOOMMMM ADDED 'WHERE WATER'
IF(IWATER(I).EQ.1) THEN
VXNUM3(I) = VXNUM3(I) - 2. * VG0(I)
VXNUM3(I) = VAWS1(I) * VXNUM3(I)
ENDIF
ENDIF
9012 CONTINUE
ENDIF
C COMPUTE D LOG ZETA
IF (ITYPE.GE.2) THEN
DO 9014 I = 1,IRUN
VXNUM(I) = VRIB2(I) - VRIB1(I)
IF( (VX0(I).GT.XX0MAX).AND.(VXNUM(I).GE.0.) )VXNUM(I) = 0.
VXNUM(I) = VXNUM1(I) * VXNUM(I)
9014 CONTINUE
ENDIF
IF ( ITYPE.EQ.2 )THEN
DO 9016 I = 1,IRUN
VDZETA1(I) = VDZETA(I)
VXNUM(I) = VXNUM(I) + VXNUM3(I) * ( VWS2(I) - VWS1(I) )
9016 CONTINUE
ENDIF
IF (ITYPE.GE.3) THEN
DO 9018 I = 1,IRUN
VDZETA1(I) = VXNUM(I)
IF(LWATER.AND.(IWATER(I).EQ.1)) VXNUM(I) = VXNUM(I) + VXNUM2(I)
IF ( VDEN(I) .LT.0.1 ) VDEN(I) = 0.1
9018 CONTINUE
ENDIF
IF (ITYPE.GE.2) THEN
DO 9020 I = 1,IRUN
VDZETA(I) = VXNUM(I) / VDEN(I)
9020 CONTINUE
ENDIF
IF (ITYPE.GE.3) THEN
DO 9022 I = 1,IRUN
IF ( (VRIB2(I).EQ.0.) .OR. (VDZETA(I).LE.-1.) ) THEN
VINT1(I) = 1
ELSE
VINT1(I) = 0
ENDIF
IF ( VINT1(I).EQ.1 ) VDZETA(I) = VDZETA1(I)
9022 CONTINUE
ENDIF
IF (ITYPE.EQ.2) THEN
DO 9024 I = 1,IRUN
VDZETA2(I) = VDZETA(I) + VDZETA1(I)
IF ( (VRIB2(I).EQ.0.) .OR. (VDZETA2(I).LE.-1.) ) THEN
VINT1(I) = 1
ELSE
VINT1(I) = 0
ENDIF
IF(VINT1(I).EQ.1)VDZETA(I)=VXNUM1(I)*VRIB2(I) - 1. - VDZETA1(I)
9024 CONTINUE
ENDIF
C COMPUTE D LOG Z0
IF ( LWATER .AND. (ITYPE.GE.3) )THEN
DO 9026 I = 1,IRUN
IF( IWATER(I).EQ.1 ) THEN
VZCOEF2(I) = VG(I) * VDXPSIM(I)
VDZ0(I) = VDZ0(I) - VZCOEF2(I) * VDZETA(I)
ENDIF
9026 CONTINUE
ENDIF
IF ( LWATER .AND. (ITYPE.EQ.5) ) THEN
DO 9028 I = 1,IRUN
IF(IWATER(I).EQ.1) VZCOEF1(I) = VG(I) * VAWS1(I)
9028 CONTINUE
ENDIF
IF ( LWATER .AND. (ITYPE.EQ.2) ) THEN
DO 9030 I = 1,IRUN
IF (IWATER(I).EQ.1) VDZ0(I) =
1 VZCOEF1(I) * ( VWS2(I) - VWS1(I) ) - VZCOEF2(I) * VDZETA(I)
9030 CONTINUE
ENDIF
C CALCULATE D PSIM AND D PSIH
IF ( (ITYPE.EQ.1) .AND. LWATER ) THEN
DO 9032 I = 1,IRUN
IF (IWATER(I).EQ.1) THEN
VDPSIM(I) = - VDZ0(I) * VAZ0(I)
VDPSIH(I) = VDPSIM(I)
ENDIF
9032 CONTINUE
ENDIF
IF (ITYPE.GE.3) THEN
DO 9034 I = 1,IRUN
VDPSIM(I) = VDX(I) * VDZETA(I)
VDPSIH(I) = VDY(I) * VDZETA(I)
IF ( LWATER ) THEN
IF (IWATER(I).EQ.1 ) THEN
VDPSIM(I) = VDPSIM(I) - VX0(I) * VDZ0(I)
VDPSIH(I) = VDPSIH(I) - VY0(I) * VDZ0(I)
ENDIF
ENDIF
9034 CONTINUE
ENDIF
C PREVENT OVERCORRECTION OF PSIM OR PSIH FOR UNSTABLE CASE
IF (ITYPE.GE.4) THEN
DO 9036 I = 1,IRUN
VDPSIMC(I) = -0.9 - VDPSIM(I) * VAPSIM(I)
VDPSIHC(I) = -0.9 * VPSIH(I) - VDPSIH(I)
IF ( VDPSIMC(I).GT.0. ) THEN
VINT1(I) = 1
ELSE
VINT1(I) = 0
ENDIF
IF ( VDPSIHC(I).GT.0. ) THEN
VINT2(I) = 1
ELSE
VINT2(I) = 0
ENDIF
VDZETA1(I) = 0.
IF(VINT1(I).EQ.1) VDZETA1(I) = VDPSIMC(I) / VDXPSIM(I)
IF((VINT1(I).EQ.1).OR.(VINT2(I).EQ.1)) VTEMPLIN(I) =
1 VDY(I) + VY0(I) * VG(I) * VDXPSIM(I)
IF (VINT2(I).EQ.1) then
VDZETA2(I) = VDPSIHC(I) / VTEMPLIN(I)
IF ( VDZETA2(I).LT.VDZETA1(I) ) VDZETA1(I) = VDZETA2(I)
endif
IF((VINT1(I).EQ.1).OR.(VINT2(I).EQ.1)) THEN
VDZETA(I) = VDZETA1(I) + VDZETA(I)
VDPSIM(I) = VDPSIM(I) + VDX(I) * VR1MG0(I) * VDZETA1(I)
VDPSIH(I) = VDPSIH(I) + VTEMPLIN(I) * VDZETA1(I)
IF ( IWATER(I).EQ.1 )
1 VDZ0(I) = VDZ0(I) - VG(I) * VDXPSIM(I) * VDZETA1(I)
ENDIF
9036 CONTINUE
ENDIF
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
SUBROUTINE ZCSUB (VUSTAR,VDZSEA,IWATER,LDZSEA,IRUN,VZSEA)
C**********************************************************************
C FUNCTION ZSEA
C PURPOSE
C COMPUTES Z0 AS A FUNCTION OF USTAR OVER WATER SURFACES
C USAGE
C CALLED BY SFCFLX
C DESCRIPTION OF PARAMETERS
C USTAR - INPUTED VALUE OF SURFACE-STRESS VELOCITY
C DZSEA - OUTPUTED VALUE OF DERIVATIVE D(ZSEA)/D(USTAR)
C WATER - INPUTED BIT VECTOR TO DETERMINE WATER POINTS
C LDZSEA- LOGICAL FLAG TO DETERMINE IF DZSEA SHOULD BE COMPUTED
C ZSEA - OUTPUTED VALUE OF ROUGHNESS LENGTH
C SUBPROGRAMS NEEDED
C NONE
C RECORD OF MODIFICATIONS
C REMARKS:
C COMPUTE ROUGHNESS LENGTH FOR OCEAN POINTS
C BASED ON FUNCTIONS OF LARGE AND POND
C AND OF KONDO --- DESIGNED FOR K = .4
C *********************************************************************
implicit none
C Argument List Delcarations
integer irun
_RL VZSEA(IRUN),VUSTAR(IRUN),VDZSEA(IRUN)
integer IWATER(IRUN)
LOGICAL LDZSEA
C Local Variables
_RL USTMX1,USTMX2,USTMX3
PARAMETER ( USTMX1 = 1.14973 )
PARAMETER ( USTMX2 = 0.381844 )
PARAMETER ( USTMX3 = 0.0632456)
_RL AA(IRUN,5),TEMP(IRUN)
integer INT2(IRUN),INT3(IRUN),INT4(IRUN)
integer i,k
_RL AA1(5),AA2(5),AA3(5),AA4(5)
DATA AA1/.2030325E-5,0.0,0.0,0.0,0.0/
DATA AA2/-0.402451E-08,0.239597E-04,0.117484E-03,0.191918E-03,
1 0.395649E-04/
DATA AA3/-0.237910E-04,0.228221E-03,-0.860810E-03,0.176543E-02,
1 0.784260E-04/
DATA AA4/-0.343228E-04,0.552305E-03,-0.167541E-02,0.250208E-02,
1 -0.153259E-03/
C**********************************************************************
C***** LOWER CUTOFF CONDITION FOR USTAR ***
C**********************************************************************
DO 9000 I = 1,IRUN
IF(VUSTAR(I) .LT. 1.e-6)THEN
INT3(I) = 1
ELSE
INT3(I) = 0
ENDIF
9000 CONTINUE
DO 9002 I = 1,IRUN
IF(INT3(I).EQ.1) VUSTAR(I) = 1.e-6
9002 CONTINUE
C***********************************
C***** LOAD THE ARRAY A(I,K) *****
C***********************************
DO 9004 I = 1,IRUN
IF( (VUSTAR(I) .GT. USTMX1) .AND. (IWATER(I).EQ.1) ) THEN
INT4(I) = 1
ELSE
INT4(I) = 0
ENDIF
9004 CONTINUE
DO 9006 I = 1,IRUN
IF(VUSTAR(I) .GT. USTMX2) THEN
INT3(I) = 1
ELSE
INT3(I) = 0
ENDIF
9006 CONTINUE
DO 9008 I = 1,IRUN
IF(VUSTAR(I) .GE. USTMX3) THEN
INT2(I) = 1
ELSE
INT2(I) = 0
ENDIF
9008 CONTINUE
DO 100 K=1,5
DO 9010 I = 1,IRUN
AA(I,K) = AA1(K)
IF( INT2(I).EQ.1 ) AA(I,K) = AA2(K)
IF( INT3(I).EQ.1 ) AA(I,K) = AA3(K)
IF( INT4(I).EQ.1 ) AA(I,K) = AA4(K)
9010 CONTINUE
100 CONTINUE
C********************************************************
C***** EVALUATE THE ENHANCED POLYNOMIAL FOR ZSEA *****
C********************************************************
DO 9012 I = 1,IRUN
VDZSEA(I) = ( AA(I,4) + AA(I,5) * VUSTAR(I) ) * VUSTAR(I)
VZSEA(I) = AA(I,2) + ( AA(I,3) + VDZSEA(I) ) * VUSTAR(I)
TEMP(I) = AA(I,1) / VUSTAR(I)
VZSEA(I) = VZSEA(I) + TEMP(I)
9012 CONTINUE
C**********************************************************************
C***** EVALUATE THE DERIVATIVE DZSEA IF LDZSEA IS TRUE ***
C**********************************************************************
IF( LDZSEA ) THEN
DO 9014 I = 1,IRUN
VDZSEA(I) = 3. * VDZSEA(I) -(AA(I,4)*VUSTAR(I) - AA(I,3))
VDZSEA(I) = VDZSEA(I) * VUSTAR(I) - TEMP(I)
9014 CONTINUE
ENDIF
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
SUBROUTINE SEAICE ( nocean, timstp, hice,
& eturb, dedtc,
& hsturb, dhsdtc,
& qice, dqice,
& swnet, lwnet, dst4,
& pke, seaic, tc, qa )
implicit none
integer nocean
_RL timstp
_RL eturb(nocean),dedtc(nocean),hsturb(nocean),dhsdtc(nocean)
_RL swnet(nocean),lwnet(nocean), dst4(nocean)
_RL qice(nocean),dqice(nocean)
_RL pke(nocean), tc(nocean), qa(nocean)
_RL seaic(nocean)
C rho*C = 1.93e6 J/(m**3 K) ; Peixoto & Oort
_RL rhoC
parameter (rhoC = 1.93e6)
_RL faceps,getcon,latent,codt,deltg,hice
integer i
faceps = getcon('EPSFAC')
latent = getcon('HEATI') * getcon('CALTOJ')
C Note hice is in centimeters
codt = rhoC * (hice/100) / timstp
C Update TC and QA
C ----------------
do i =1,nocean
if( seaic(i).gt.0.0 ) then
deltg = ( swnet(i)-lwnet(i)-latent*eturb(i)-hsturb(i)+qice(i) )
& / ( codt+dst4(i)+latent*dedtc(i)+dhsdtc(i)-dqice(i) )
qa(i) = qa(i) + (faceps*qa(i)/(tc(i)*tc(i)))*deltg
tc(i) = tc(i) + deltg
endif
enddo
return
end