C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_turb.F,v 1.44 2006/06/18 00:10:40 molod Exp $
C $Name: $
#include "FIZHI_OPTIONS.h"
subroutine TURBIO (im,jm,nlay,istrip,nymd,nhms,bi,bj
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
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 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 ndiagsiz- number of diagnostic 2-D arrays allocated
c ndlsm - number of tile diagnostic
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
#ifdef ALLOW_USE_MPI
#include "mpif.h"
#endif
integer im,jm,nlay,istrip,nymd,nhms,bi,bj,ndturb,nltop
integer ntracers, ptracers
integer nchp,nchptot,nchplnd
_RL ptop
_RL pz(im,jm),uz(im,jm,nlay),vz(im,jm,nlay),tz(im,jm,nlay)
_RL qz(im,jm,nlay,ntracers)
_RL plz(im,jm,nlay),plze(im,jm,nlay+1),dpres(im,jm,nlay)
_RL pkht(im,jm,nlay+1),pkz(im,jm,nlay)
_RL ctmt(nchp),xxmt(nchp),yymt(nchp),zetamt(nchp)
_RL xlmt(nchp,nlay),khmt(nchp,nlay),tke(nchp,nlay)
_RL tgz(im, jm),fracland(im,jm)
integer landtype(im,jm)
_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)
logical tprof
_RL duturb(im,jm,nlay),dvturb(im,jm,nlay)
_RL dtturb(im,jm,nlay),dqturb(im,jm,nlay,ntracers)
_RL st4(im,jm),dst4(im,jm)
_RL radswg(im,jm),radswt(im,jm),radlwg(im,jm)
_RL fdifpar(im,jm),fdirpar(im,jm)
_RL rainlsp(im,jm),rainconv(im,jm),snowfall(im,jm)
_RL tempref (im,jm)
integer imstturblw, imstturbsw
_RL qliqavesw(im,jm,nlay),qliqavelw(im,jm,nlay)
_RL fccavelw (im,jm,nlay),fccavesw (im,jm,nlay)
_RL qqgrid (im,jm,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,nlay)
_RL fcctmp(im,jm,nlay)
_RL tmpdiag(im,jm)
_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)
integer ndlsm
parameter ( ndlsm = 1)
_RL qdiaglsm(nchp,ndlsm)
_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,j,L,nn,nt
integer nocean, nice
integer ndmoist,time_left,ndum
integer ntracedim
_RL dtfac,timstp2,sum0
C logical begin flag - set to true to indicate a cold start
logical qbeg
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
c
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',ndum,ndum,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)
c
rmu = fmu * tref * rgas / pref
cltj10 = 10. * caltoj
atimstp = 1. / timstp
tice = getcon('FREEZING-POINT')
c **********************************************************************
c Check for Cold Start (if QQ is zero everywhere)
c **********************************************************************
qbeg = .false.
sum0 = 0.0
do L=1,nlay
do n=1,nchptot
sum0 = sum0 + tke(n,L)
enddo
enddo
#ifdef ALLOW_USE_MPI
call MPI_ALLREDUCE(sum0,const,1,mpi_double_precision,mpi_sum,
. mpi_comm_world,n)
#else
const = sum0
#endif
if( const.eq.0.0 ) then
qbeg = .true.
if( myid.eq.1 .and. bi.eq.1 ) then
print *
print *, 'Warning!'
print *, 'Turbulent Kinetic Energy has not been initialized.'
print *, 'Cold-Start will use Level 2.0 Turbulence.'
print *
endif
endif
c **********************************************************************
c Initialization
c **********************************************************************
c Initialize diagnostic for ground temperature change
c ---------------------------------------------------
if(diagnostics_is_on('DTG ',myid) ) then
do j =1,jm
do i =1,im
tmpdiag(i,j) = -tgz(i,j)
enddo
enddo
call DIAGNOSTICS_FILL(tmpdiag,'DTG ',0,1,-3,bi,bj,myid)
endif
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
c
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)
call ASTRO ( nymd,nhms,lats,lons,istrip,cosz,ra )
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
ccc newsnow(i) = newsnow(i)*dtfac
ccc raincon(i) = raincon(i)*dtfac
ccc 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