/[MITgcm]/MITgcm/pkg/fizhi/fizhi_turb.F
ViewVC logotype

Annotation of /MITgcm/pkg/fizhi/fizhi_turb.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.31 - (hide annotations) (download)
Thu Feb 24 21:05:23 2005 UTC (19 years, 4 months ago) by molod
Branch: MAIN
Changes since 1.30: +8 -1 lines
Bug fix -- initialize bottom edge velocities (to zero) before trbflx call

1 molod 1.31 C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_turb.F,v 1.30 2005/02/22 22:51:42 molod Exp $
2 molod 1.1 C $Name: $
3 molod 1.2
4 molod 1.14 #include "FIZHI_OPTIONS.h"
5 molod 1.6 subroutine turbio (im,jm,nlay,istrip,nymd,nhms,bi,bj
6 molod 1.10 1 ,ndturb,ptop, pz, uz, vz, tz, qz, ntracers,ptracers
7 molod 1.1 2 ,plz,plze,dpres,pkht,pkz,ctmt,xxmt,yymt,zetamt,xlmt,khmt,tke
8 molod 1.5 3 ,tgz,fracland,landtype
9 molod 1.1 4 ,tcanopy,ecanopy,tdeep,swetshal,swetroot,swetdeep,snodep,capac
10 molod 1.13 5 ,nchp,nchptot,nchplnd,chfr,chlt,chlon,igrd,ityp,alai,agrn,thkz
11     6 ,tprof
12 molod 1.1 8 ,duturb, dvturb, dtturb,dqturb,radlwg,st4,dst4,radswg,radswt
13     9 ,fdifpar,fdirpar,rainlsp,rainconv,snowfall,tempref
14     1 ,imstturblw,imstturbsw,qliqavelw,qliqavesw,fccavelw,fccavesw
15 molod 1.7 2 ,qqgrid,myid)
16 molod 1.1 c-----------------------------------------------------------------------
17     c subroutine turbio - model interface routine to trbflx, the turbulence
18     c parameterization, and tile, the land surface parameterization
19     c
20     c input:
21     c im - number of points in the longitude direction
22     c jm - number of points in the latitude direction
23     c nlay - number of vertical levels
24     c istrip - number of horizontal points to be handled at a time on
25     c nymd - year and date integer in YYMMDD format (ie, 790212)
26     c nhms - date and time integer in HHMMSS format (ie, 123000)
27     c ndturb - turbulence time step integer in HHMMSS format
28     c ptop - model top pressure - rigid lid assumed - real
29     c pz - surface pressure minus ptop in mb - real[lon,lat]
30     c uz - zonal wind in m/sec - real[lon,lat,level]
31     c vz - meridional wind in m/sec - real[lon,lat,level]
32     c tz - model theta (theta [deg K]/p0**k) - real[lon,lat,level]
33     c qz - specific humidity in kg/kg - real[lon,lat,level]
34     c ntracers- total number of tracers - integer
35     c ptracers- number of permanent tracers - integer
36     c pkht - pressure[mb]**k at bottom edges of levels - real[lon,lat,level]
37     c fracland- not being used - real[lon,lat]
38     c landtype- not being used - integer[lon,lat]
39     c nchp - nchplnd<nchp - total no chips (ocean too) - integer
40     c nchplnd - <=nchp - number of land chips - integer
41     c chfr - chip fraction - real[nchp]
42     c chlt - tile space latitude array - real[nchp]
43     c chlon - tile space longitude array - real[nchp]
44     c igrd - tile space grid number - integer[nchp]
45     c ityp - tile space vegetation type - integer[nchp]
46     c alai - leaf area index - real[nchp]
47     c agrn - greenness fraction - real[nchp]
48     c thkz - sea ice thickness in m (0. for no ice) - real[lon,lat]
49     c tprof - logical flag for point by point diagnostic output
50     c ndiagsiz- number of diagnostic 2-D arrays allocated
51     c ndlsm - number of tile diagnostic
52     c radlwg - net longwave flux at ground (up-down) in w/m**2 - real[lon,lat]
53     c st4 - upward longwave flux at ground in w/m**2 - real[lon,lat]
54     c dst4 - delta-sigma-T**4, ie, derivative of upward lw flux at
55     c ground with respect to ground Temperature - real[lon,lat]
56     c radswg - net shortwave flux at ground (down-up) NON-DIM - real[lon,lat]
57     c {NOTE: this field is divided by the incident shortwave
58     c at the top of the atmosphere to non-dimensionalize]
59     c radswt - incident shortwave at top of atmos in W/m**2 - real[lon,lat]
60     c fdifpar - incident diffuse-beam PAR at surface in W/m**2 - real[lon,lat]
61     c fdirpar - incident direct-beam PAR at surface in W/m**2 - real[lon,lat]
62     c rainlsp - large-scale (frontal,supersat) rainfall in mm/sec - real[lon,lat]
63     c rainconv- convective rainfall rate in mm/sec - real[lon,lat]
64     c snowfall- total snowfall rate in mm/sec - real[lon,lat]
65     c updated:
66     c tke - turbulent k.e. in m**2/s**2 - real[tiles,levels]
67     c tgz - surface skin temperature in deg K - real[lon,lat]
68     c tcanopy - canopy temperature in deg K real[tiles]
69     c (sea surface temp over the ocean tiles)
70     c ecanopy - canopy vapor pressure in mb real[tiles]
71     c (qstar at tground over the sea ice and ocean tiles)
72     c tdeep - deep soil temp in deg K real[tiles]
73     c swetshal- shallow level moisture field capacity fraction real[tiles]
74     c swetroot- root level moisture field capacity fraction real[tiles]
75     c swetdeep- deep soil level moisture field capacity fraction real[tiles]
76     c snodep - depth of snow pack in cm liquid water equiv real[tiles]
77     c capac - leaf canopy water reservoir in cm real[tiles]
78     c output:
79     c duturb - change in zonal wind component due to turbulent processes
80     c per unit time in m/sec**2 - real[lon,lat,levels]
81     c dvturb - change in meridional wind component due to turbulent processes
82     c per unit time in m/sec**2 - real[lon,lat,levels]
83     c dtturb - change in (model theta*pi) due to turbulent processes
84     c per unit time - real[lon,lat,levels] !! pi is pressure-ptop
85     c dqturb - change in (specific humidity*pi) due to turbulent processes
86     c per unit time - real[lon,lat,levels] !! pi is pressure-ptop
87     c qliqavelw - Moist Turbulence Liquid Water for Longwave - real[lon,lat,levels]
88     c qliqavesw - Moist Turbulence Liquid Water for Shortwave - real[lon,lat,levels]
89     c fccavelw - Moist Turbulence Cloud Fraction for Longwave - real[lon,lat,levels]
90     c fccavesw - Moist Turbulence Cloud Fraction for Shortwave - real[lon,lat,levels]
91     c qqgrid - Gridded Turbulent Kinetic Energy - real[lon,lat,levels]
92     c-----------------------------------------------------------------------
93     implicit none
94    
95 molod 1.7 #ifdef ALLOW_USE_MPI
96 molod 1.10 #include "mpif.h"
97 molod 1.7 #endif
98    
99 molod 1.2 #ifdef ALLOW_DIAGNOSTICS
100 molod 1.9 #include "SIZE.h"
101 molod 1.28 #include "DIAGNOSTICS_SIZE.h"
102     #include "DIAGNOSTICS.h"
103 molod 1.2 #endif
104 molod 1.1
105 molod 1.7 integer im,jm,nlay,istrip,nymd,nhms,bi,bj,ndturb
106 molod 1.10 integer ntracers, ptracers
107 molod 1.13 integer nchp,nchptot,nchplnd
108 molod 1.14 _RL ptop
109     _RL pz(im,jm),uz(im,jm,nlay),vz(im,jm,nlay),tz(im,jm,nlay)
110     _RL qz(im,jm,nlay,ntracers)
111     _RL plz(im,jm,nlay),plze(im,jm,nlay+1),dpres(im,jm,nlay)
112     _RL pkht(im,jm,nlay+1),pkz(im,jm,nlay)
113     _RL ctmt(nchp),xxmt(nchp),yymt(nchp),zetamt(nchp)
114     _RL xlmt(nchp,nlay),khmt(nchp,nlay),tke(nchp,nlay)
115     _RL tgz(im, jm),fracland(im,jm)
116 molod 1.7 integer landtype(im,jm)
117 molod 1.14 _RL tcanopy(nchp),tdeep(nchp),ecanopy(nchp),swetshal(nchp)
118     _RL swetroot(nchp),swetdeep(nchp),snodep(nchp),capac(nchp)
119     _RL chfr(nchp),chlt(nchp),chlon(nchp)
120 molod 1.7 integer igrd(nchp),ityp(nchp)
121 molod 1.14 _RL alai(nchp),agrn(nchp),thkz(im,jm)
122 molod 1.1 logical tprof
123 molod 1.14 _RL duturb(im,jm,nlay),dvturb(im,jm,nlay)
124     _RL dtturb(im,jm,nlay),dqturb(im,jm,nlay,ntracers)
125     _RL st4(im,jm),dst4(im,jm)
126     _RL radswg(im,jm),radswt(im,jm),radlwg(im,jm)
127     _RL fdifpar(im,jm),fdirpar(im,jm)
128     _RL rainlsp(im,jm),rainconv(im,jm),snowfall(im,jm)
129     _RL tempref (im,jm)
130 molod 1.7 integer imstturblw, imstturbsw
131 molod 1.14 _RL qliqavesw(im,jm,nlay),qliqavelw(im,jm,nlay)
132     _RL fccavelw (im,jm,nlay),fccavesw (im,jm,nlay)
133     _RL qqgrid (im,jm,nlay)
134 molod 1.7 integer myid
135    
136     C Local Variables
137 molod 1.1
138     integer numstrips
139     integer ijall
140 molod 1.14 _RL fmu,hice,tref,pref,cti,ed
141 molod 1.1 C Set fmu and ed to zero for no background diffusion
142     parameter ( fmu = 0.00000 )
143     parameter ( hice = 300. )
144     parameter ( tref = 258. )
145     parameter ( pref = 500. )
146     parameter ( cti = 0.0052 )
147     parameter ( ed = 0.0 )
148    
149 molod 1.14 _RL qliqtmp(im,jm,nlay)
150     _RL fcctmp(im,jm,nlay)
151     _RL tmpdiag(im,jm)
152     _RL thtgz(im*jm)
153 molod 1.29 logical diagnostics_is_on
154     external diagnostics_is_on
155 molod 1.1
156     integer nland
157 molod 1.14 _RL alwcoeff(nchp),blwcoeff(nchp)
158     _RL netsw(nchp)
159     _RL cnvprec(nchp),lsprec(nchp)
160     _RL snowprec(nchp)
161     _RL pardiff(nchp),pardirct(nchp)
162     _RL pmsc(nchp)
163     _RL netlw(nchp)
164     _RL sqscat(nchp), rsoil1(nchp)
165     _RL rsoil2(nchp)
166     _RL rdc(nchp),u2fac(nchp)
167     _RL z2ch(nchp)
168     _RL zoch(nchp),cdrc(nchp)
169     _RL cdsc(nchp)
170     _RL dqsdt(nchp)
171     _RL tground(nchp),qground(nchp)
172     _RL utility(nchp)
173     _RL qice(nchp)
174     _RL dqice(nchp)
175 molod 1.1
176 molod 1.14 _RL dumsc(nchp,nlay),dvmsc(nchp,nlay)
177     _RL dtmsc(nchp,nlay),dqmsc(nchp,nlay,ntracers)
178 molod 1.1
179 molod 1.14 _RL shg(nchp),z0(nchp),icethk(nchp)
180 molod 1.1 integer water(nchp)
181    
182 molod 1.14 _RL lats(istrip),lons(istrip),cosz(istrip),icest(istrip)
183     _RL rainls(istrip),raincon(istrip),newsnow(istrip)
184     _RL pardf(istrip),pardr(istrip),swnet(istrip)
185 molod 1.15 _RL hlwdwn(istrip),alwrad(istrip),blwrad(istrip)
186     _RL tmpnlay(istrip)
187 molod 1.14 _RL laistrip(istrip),grnstrip(istrip),z2str(istrip),cd(istrip)
188     _RL scatstr(istrip), rs1str(istrip), rs2str(istrip)
189     _RL rdcstr(istrip),u2fstr(istrip),dqsdtstr(istrip)
190     _RL eturb(istrip),dedqa(istrip),dedtc(istrip)
191     _RL hsturb(istrip),dhsdqa(istrip),dhsdtc(istrip)
192     _RL savetc(istrip),saveqa(istrip),lwstrip(istrip)
193     _RL chfrstr(istrip),psurf(istrip),shgstr(istrip)
194 molod 1.1 integer types(istrip),igrdstr(istrip)
195 molod 1.14 _RL evap(istrip),shflux(istrip),runoff(istrip),bomb(istrip)
196     _RL eint(istrip),esoi(istrip),eveg(istrip),esno(istrip)
197     _RL smelt(istrip),hlatn(istrip),hlwup(istrip),gdrain(istrip)
198     _RL runsrf(istrip),fwsoil(istrip),evpot(istrip)
199     _RL strdg1(istrip),strdg2(istrip),strdg3(istrip),strdg4(istrip)
200     _RL strdg5(istrip),strdg6(istrip),strdg7(istrip),strdg8(istrip)
201     _RL strdg9(istrip),tmpstrip(istrip),qicestr(istrip)
202     _RL dqicestr(istrip)
203    
204     _RL u(istrip,nlay+1), v(istrip,nlay+1), th(istrip,nlay+1)
205     _RL sh(istrip,nlay+1), thv(istrip,nlay+1), pe(istrip,nlay+1)
206     _RL tracers(istrip,nlay+1,ntracers)
207     _RL dpstr(istrip,nlay),pke(istrip,nlay+1)
208     _RL pk(istrip,nlay), qq(istrip,nlay), p(istrip,nlay)
209     _RL sri(istrip,nlay), skh(istrip,nlay), skm(istrip,nlay)
210     _RL stuflux(istrip,nlay), stvflux(istrip,nlay)
211     _RL sttflux(istrip,nlay), stqflux(istrip,nlay)
212     _RL frqtrb(istrip,nlay-1)
213     _RL dshdthg(istrip,nlay),dthdthg(istrip,nlay)
214     _RL dshdshg(istrip,nlay),dthdshg(istrip,nlay)
215     _RL transth(istrip,nlay), transsh(istrip,nlay)
216    
217     _RL tc(istrip),td(istrip),qa(istrip)
218     _RL swet1(istrip),swet2(istrip),swet3(istrip)
219     _RL capacity(istrip),snowdepth(istrip)
220     _RL stz0(istrip)
221     _RL stdiag(istrip)
222     _RL tends(istrip),sustar(istrip), sz0(istrip),pbldpth(istrip)
223     _RL sct(istrip), scu(istrip), swinds(istrip)
224     _RL stu2m(istrip),stv2m(istrip),stt2m(istrip),stq2m(istrip)
225     _RL stu10m(istrip),stv10m(istrip),stt10m(istrip),stq10m(istrip)
226 molod 1.1 integer stwatr(istrip)
227 molod 1.14 _RL wspeed(istrip)
228 molod 1.1
229 molod 1.15 _RL ctsave(istrip),xxsave(istrip),yysave(istrip)
230     _RL zetasave(istrip)
231 molod 1.14 _RL xlsave(istrip,nlay),khsave(istrip,nlay)
232     _RL qliq(istrip,nlay),turbfcc(istrip,nlay)
233     _RL qliqmsc(nchp,nlay),fccmsc(nchp,nlay)
234 molod 1.2
235     integer ndlsm
236 molod 1.21 parameter ( ndlsm = 1)
237 molod 1.14 _RL qdiaglsm(nchp,ndlsm)
238 molod 1.1
239 molod 1.14 _RL pi,secday,sdayopi2,rgas,akap,cp,alhl
240     _RL faceps,grav,caltoj,virtcon,getcon
241     _RL heatw,undef,timstp,delttrb,dttrb,ra
242     _RL edle,rmu,cltj10,atimstp,tice,const
243 molod 1.10 integer istnp1,istnlay,itrtrb,i,j,L,nn,nt
244 molod 1.1 integer nocean, nice
245     integer ndmoist,time_left,ndum
246     integer ntracedim
247 molod 1.25 _RL dtfac,timstp2,sum0
248 molod 1.1 C logical begin flag - set to true to indicate a cold start
249     logical qbeg
250    
251 molod 1.18 integer n,nsecf,nmonf,ndayf
252     nsecf(n) = n/10000*3600 + mod(n,10000)/100* 60 + mod(n,100)
253     nmonf(n) = mod(n,10000)/100
254     ndayf(n) = mod(n,100)
255    
256 molod 1.12 #ifdef CRAY
257     #ifdef f77
258 molod 1.1 cfpp$ expand (qsat)
259     #endif
260     #endif
261    
262     c compute variables that do not change
263     c
264 molod 1.17
265 molod 1.1 pi = 4.*atan(1.)
266     secday = getcon('SDAY')
267     sdayopi2 = getcon('SDAY') / (pi*2.)
268     rgas = getcon('RGAS')
269     akap = getcon('KAPPA')
270     cp = getcon('CP')
271     alhl = getcon('LATENT HEAT COND')
272     faceps = getcon('EPSFAC')
273     grav = getcon('GRAVITY')
274     caltoj = getcon('CALTOJ')
275     virtcon = getcon('VIRTCON')
276     heatw = getcon('HEATW')
277     undef = getcon('UNDEF')
278     ntracedim= max(ntracers-ptracers,1)
279    
280     call get_alarm ( 'moist',ndum,ndum,ndmoist,time_left )
281     timstp = nsecf(ndturb)
282     timstp2 = nsecf(ndmoist)
283     dtfac = min( 1.0, timstp/timstp2 )
284 molod 1.19
285 molod 1.1 c delttrb is the internal turbulence time step
286     c a value equal to ndturb means one internal iteration
287     delttrb = nsecf(ndturb)
288    
289     ijall = im * jm
290     istnp1 = istrip * (nlay+1)
291     istnlay = istrip * nlay
292     itrtrb = ( timstp / delttrb ) + 0.1
293     dttrb = timstp / float(itrtrb)
294     edle = ed * 0.2
295    
296     c coefficient of viscosity (background momentum diffusion)
297     c
298     rmu = fmu * tref * rgas / pref
299     cltj10 = 10. * caltoj
300     atimstp = 1. / timstp
301     tice = getcon('FREEZING-POINT')
302    
303     c **********************************************************************
304     c Check for Cold Start (if QQ is zero everywhere)
305     c **********************************************************************
306    
307     qbeg = .false.
308    
309 molod 1.25 sum0 = 0.0
310 molod 1.1 do L=1,nlay
311 molod 1.13 do n=1,nchptot
312 molod 1.25 sum0 = sum0 + tke(n,L)
313 molod 1.1 enddo
314     enddo
315 molod 1.10
316     #ifdef ALLOW_USE_MPI
317 molod 1.25 call mpi_allreduce(sum0,const,1,mpi_double_precision,mpi_sum,
318 molod 1.7 . mpi_comm_world,n)
319 molod 1.10 #else
320 molod 1.25 const = sum0
321 molod 1.10 #endif
322 molod 1.1
323     if( const.eq.0.0 ) then
324     qbeg = .true.
325 molod 1.27 if( myid.eq.1 .and. bi.eq.1 ) then
326 molod 1.1 print *
327     print *, 'Warning!'
328     print *, 'Turbulent Kinetic Energy has not been initialized.'
329     print *, 'Cold-Start will use Level 2.0 Turbulence.'
330     print *
331     endif
332     endif
333    
334 molod 1.29 call diagnostics_fill(tgz,'TGROUND ',0,1,3,bi,bj,myid)
335     c if(itground.gt.0) then
336     c do j =1,jm
337     c do i =1,im
338     c qdiag(i,j,itground,bi,bj) = qdiag(i,j,itground,bi,bj) + tgz(i,j)
339     c enddo
340     c enddo
341     c endif
342 molod 1.23
343 molod 1.1 c **********************************************************************
344     c Initialization
345     c **********************************************************************
346    
347     c Initialize diagnostic for ground temperature change
348     c ---------------------------------------------------
349     if(idtg.gt.0) then
350 molod 1.23 do j =1,jm
351     do i =1,im
352     qdiag(i,j,idtg,bi,bj) = qdiag(i,j,idtg,bi,bj) - tgz(i,j)
353     enddo
354 molod 1.1 enddo
355     endif
356    
357     c **********************************************************************
358     c entire turbulence and land surface package will run in 'tile space'
359     c do conversion of model state variables to tile space
360     c (ocean points appended to tile space land point arrays)
361     c **********************************************************************
362    
363 molod 1.13 numstrips = ((nchptot-1)/istrip) + 1
364 molod 1.1
365 molod 1.13 call grd2msc(pz(1,1),im,jm,igrd,pmsc,nchp,nchptot)
366 molod 1.16
367 molod 1.13 call grd2msc(tgz,im,jm,igrd,tground,nchp,nchptot)
368 molod 1.1 do i = 1,ijall
369     tmpdiag(i,1) = st4(i,1) + dst4(i,1)*(tgz(i,1)-tempref(i,1))
370     1 - dst4(i,1)* tgz(i,1)
371     enddo
372 molod 1.13 call grd2msc(tmpdiag,im,jm,igrd,alwcoeff,nchp,nchptot)
373 molod 1.1 do i = 1,ijall
374     tmpdiag(i,1) = dst4(i,1)
375     enddo
376 molod 1.13 call grd2msc(tmpdiag,im,jm,igrd,blwcoeff,nchp,nchptot)
377 molod 1.1 do i = 1,ijall
378     tmpdiag(i,1) = fdifpar(i,1) * radswt(i,1)
379     enddo
380 molod 1.13 call grd2msc(tmpdiag,im,jm,igrd,pardiff,nchp,nchptot)
381 molod 1.1 do i = 1,ijall
382     tmpdiag(i,1) = fdirpar(i,1) * radswt(i,1)
383     enddo
384 molod 1.13 call grd2msc(tmpdiag,im,jm,igrd,pardirct,nchp,nchptot)
385 molod 1.1 do i = 1,ijall
386     tmpdiag(i,1) = radswg(i,1) * radswt(i,1)
387     enddo
388 molod 1.13 call grd2msc(tmpdiag,im,jm,igrd,netsw,nchp,nchptot)
389 molod 1.1 do i = 1,ijall
390     tmpdiag(i,1) = radlwg(i,1) + dst4(i,1)*(tgz(i,1)-tempref(i,1))
391     enddo
392 molod 1.13 call grd2msc(tmpdiag,im,jm,igrd,netlw,nchp,nchptot)
393     call grd2msc(thkz,im,jm,igrd,icethk,nchp,nchptot)
394     call grd2msc(rainlsp,im,jm,igrd,lsprec,nchp,nchptot)
395     call grd2msc(rainconv,im,jm,igrd,cnvprec,nchp,nchptot)
396     call grd2msc(snowfall,im,jm,igrd,snowprec,nchp,nchptot)
397 molod 1.1
398     C Call chpprm to get non-varying vegetation and soil characteristics
399    
400     call chpprm(nymd,nhms,nchp,nchplnd,chlt,ityp,alai,
401     1 agrn,zoch,z2ch,cdrc,cdsc,sqscat,u2fac,rsoil1,rsoil2,rdc)
402    
403     c **********************************************************************
404     c **** surface specification ****
405     c **********************************************************************
406    
407     c set water
408    
409 molod 1.13 do i = 1,nchptot
410 molod 1.1 water(i) = 0
411     if((ityp(i).eq.100).and.(icethk(i).eq.0. ))water(i) = 1
412     enddo
413    
414     c roughness length z0
415     c
416 molod 1.13 do i =1,nchptot
417 molod 1.1 if (icethk(i).gt.0.) then
418     z0(i) = 1.e-4
419     else if (ityp(i).eq.100) then
420     z0(i) = 3.e-4
421     else
422     z0(i) = zoch(i)
423     endif
424     enddo
425    
426     c Fill Array Tground with canopy temperatures over land tiles
427     c (it has sst from the tgz array over the sea ice and ocean tiles)
428    
429     do i = 1,nchplnd
430     tground(i) = tcanopy(i)
431     enddo
432    
433     C value of sh at ground
434     C ---------------------
435 molod 1.13 do I =1,nchptot
436 molod 1.1 utility(I) = pmsc(i) + ptop
437     call qsat ( tground(i),utility(i),shg(i),dqsdt(i),.true. )
438     enddo
439    
440     c Fill Array Qground with canopy air specific humidity over land tiles
441     c (it has qstar at tground over the sea ice and ocean tiles)
442    
443     do i = 1,nchplnd
444     qground(i) = ecanopy(i)
445     enddo
446 molod 1.13 do i = nchplnd+1,nchptot
447 molod 1.1 qground(i) = shg(i)
448     enddo
449    
450     c Fill Array Swetshal with Value 1. over oceans and sea ice
451 molod 1.13 do i = nchplnd+1,nchptot
452 molod 1.1 swetshal(i) = 1.
453     enddo
454    
455     c compute heat conduction through ice
456     c -----------------------------------
457     const = ( cti / hice ) * cltj10
458 molod 1.13 do i =1,nchptot
459 molod 1.1 qice(i) = 0.0
460     dqice(i) = 0.0
461     if( icethk(i).gt.0.0 ) then
462     qice(i) = const*(tice-tground(i))
463     dqice(i) = -const
464     endif
465     enddo
466    
467     if( iqice.gt.0 ) then
468 molod 1.3 do i =1,ijall
469     tmpdiag(i,1) = 0.0
470     enddo
471 molod 1.13 call msc2grd (igrd,chfr,qice,nchp,nchptot,fracland,tmpdiag,im,jm)
472 molod 1.23 do j =1,jm
473     do i =1,im
474     qdiag(i,j,iqice,bi,bj) = qdiag(i,j,iqice,bi,bj) + tmpdiag(i,j)
475     enddo
476 molod 1.3 enddo
477 molod 1.1 nqice = nqice + 1
478     endif
479    
480     c**********************************************************************
481     c loop over regions
482     c**********************************************************************
483    
484     do 2000 nn = 1, numstrips
485    
486     call strip2tile(uz,igrd,u,nchp,ijall,istrip,nlay,nn)
487     call strip2tile(vz,igrd,v,nchp,ijall,istrip,nlay,nn)
488     call strip2tile(tz,igrd,th,nchp,ijall,istrip,nlay,nn)
489     call strip2tile(qz(1,1,1,1),igrd,sh,nchp,ijall,istrip,nlay,nn)
490     call strip2tile(dpres,igrd,dpstr,nchp,ijall,istrip,nlay,nn)
491     call strip2tile(plz,igrd,p,nchp,ijall,istrip,nlay,nn)
492     call strip2tile(plze,igrd,pe,nchp,ijall,istrip,nlay+1,nn)
493     call strip2tile(pkz,igrd,pk,nchp,ijall,istrip,nlay,nn)
494     call strip2tile(pkht,igrd,pke,nchp,ijall,istrip,nlay+1,nn)
495 molod 1.20 c do nt = 1,ntracers-ptracers
496     c call strip2tile(qz(1,1,1,ptracers+nt),igrd,tracers(1,1,nt),nchp,
497     c 1 ijall,istrip,nlay,nn)
498     c enddo
499 molod 1.1
500 molod 1.13 call stripit (z0,stz0,nchptot,nchp,istrip,1,nn)
501     call stripit (tground,th(1,nlay+1),nchptot,nchp,istrip,1,nn)
502     call stripit (pmsc,pe(1,nlay+1),nchptot,nchp,istrip,1,nn)
503     call stripit (tke,qq,nchptot,nchp,istrip,nlay-1,nn)
504     call stripit (ctmt,ctsave,nchptot,nchp,istrip,1,nn)
505     call stripit (xxmt,xxsave,nchptot,nchp,istrip,1,nn)
506     call stripit (yymt,yysave,nchptot,nchp,istrip,1,nn)
507     call stripit (zetamt,zetasave,nchptot,nchp,istrip,1,nn)
508     call stripit (xlmt,xlsave,nchptot,nchp,istrip,nlay,nn)
509     call stripit (khmt,khsave,nchptot,nchp,istrip,nlay,nn)
510     call stripitint (water,stwatr,nchptot,nchp,istrip,1,nn)
511    
512     call stripitint (igrd,igrdstr,nchptot,nchp,istrip,1,nn)
513     call stripit (chfr,chfrstr,nchptot,nchp,istrip,1,nn)
514     call stripit (icethk,icest,nchptot,nchp,istrip,1,nn)
515     call stripit (pardiff,pardf,nchptot,nchp,istrip,1,nn)
516     call stripit (pardirct,pardr,nchptot,nchp,istrip,1,nn)
517     call stripit (chlt,lats,nchptot,nchp,istrip,1,nn)
518     call stripit (chlon,lons,nchptot,nchp,istrip,1,nn)
519     call stripit (lsprec,rainls,nchptot,nchp,istrip,1,nn)
520     call stripit (cnvprec,raincon,nchptot,nchp,istrip,1,nn)
521     call stripit (snowprec,newsnow,nchptot,nchp,istrip,1,nn)
522     call stripit (netsw,swnet,nchptot,nchp,istrip,1,nn)
523     call stripit (netlw,lwstrip,nchptot,nchp,istrip,1,nn)
524     call stripit (alwcoeff,alwrad,nchptot,nchp,istrip,1,nn)
525     call stripit (blwcoeff,blwrad,nchptot,nchp,istrip,1,nn)
526     call stripit (alai,laistrip,nchptot,nchp,istrip,1,nn)
527     call stripit (agrn,grnstrip,nchptot,nchp,istrip,1,nn)
528     call stripit (z2ch,z2str,nchptot,nchp,istrip,1,nn)
529     call stripit (sqscat,scatstr,nchptot,nchp,istrip,1,nn)
530     call stripit (rsoil1,rs1str,nchptot,nchp,istrip,1,nn)
531     call stripit (rsoil2,rs2str,nchptot,nchp,istrip,1,nn)
532     call stripit (rdc,rdcstr,nchptot,nchp,istrip,1,nn)
533     call stripit (u2fac,u2fstr,nchptot,nchp,istrip,1,nn)
534     call stripit (shg,shgstr,nchptot,nchp,istrip,1,nn)
535     call stripit (dqsdt,dqsdtstr,nchptot,nchp,istrip,1,nn)
536     call stripit ( qice, qicestr,nchptot,nchp,istrip,1,nn)
537     call stripit (dqice,dqicestr,nchptot,nchp,istrip,1,nn)
538     call stripitint (ityp,types,nchptot,nchp,istrip,1,nn)
539    
540     call stripit (tground,tc,nchptot,nchp,istrip,1,nn)
541     call stripit (tdeep,td,nchptot,nchp,istrip,1,nn)
542     call stripit (qground,qa,nchptot,nchp,istrip,1,nn)
543     call stripit (swetshal,swet1,nchptot,nchp,istrip,1,nn)
544     call stripit (swetroot,swet2,nchptot,nchp,istrip,1,nn)
545     call stripit (swetdeep,swet3,nchptot,nchp,istrip,1,nn)
546     call stripit (snodep,snowdepth,nchptot,nchp,istrip,1,nn)
547     call stripit (capac,capacity,nchptot,nchp,istrip,1,nn)
548 molod 1.1
549     call astro ( nymd,nhms,lats,lons,istrip,cosz,ra )
550    
551     c we need to count up the land, sea ice and ocean points
552     nocean = 0
553     nland = 0
554     nice = 0
555     do i = 1,istrip
556     if( types(i).lt.100 ) nland = nland + 1
557     if( types(i).eq.100 ) nocean = nocean + 1
558     if( types(i).eq.100 .and. icest(i).gt.0.0 ) nice = nice + 1
559     enddo
560    
561     c Disable following ISTRIP check for MPI version
562     c ----------------------------------------------
563     c if( (nland+nocean).ne.istrip ) then
564     c print *
565     c print *,'Error!'
566     c print *,'Problem Stripping Land/Ocean/Ice points in Turbulence'
567     c print *
568     c stop
569     c endif
570    
571 molod 1.31 c Zero out velocities at the bottom edge of the model
572     c ---------------------------------------------------
573     do i =1,istrip
574     u(i,nlay+1) = 0.
575     v(i,nlay+1) = 0.
576     enddo
577    
578 molod 1.1 c convert temperature of level nlay+1 to theta & value of sh at ground
579     c --------------------------------------------------------------------
580     do i =1,istrip
581     th(i,nlay+1) = th(i,nlay+1) / pke(i,nlay+1)
582     sh(i,nlay+1) = qa(i)
583     enddo
584    
585 molod 1.23 c if(iqg.gt.0) then
586     c do i=1,istrip
587     c tmpstrip(i) = sh(i,nlay+1)*1000
588     c enddo
589     c call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,
590     c 1 qdiag(1,1,iqg,bi,bj),ijall,1,nn,.false.)
591     c endif
592 molod 1.1
593     c value of tracers at the ground
594     c ------------------------------
595 molod 1.20 c do nt = 1,ntracers-ptracers
596     C do i = 1,istrip
597     C tracers(i,nlay+1,nt) = 0.
598     C enddo
599     C enddo
600 molod 1.1
601     c compute virtual potential temperatures
602     c --------------------------------------
603     do L = 1,nlay+1
604     do i =1,istrip
605     thv(i,L) = 1. + virtcon * sh(i,L)
606     thv(i,L) = th(i,L) * thv(i,L)
607     enddo
608     enddo
609     do i =1,istrip
610     sh(i,nlay+1) = qa(i)
611     enddo
612    
613     c zero out arrays for output of qliq and fcc
614     do L =1,nlay
615     do i =1,istrip
616     qliq(i,L) = 0.
617     turbfcc(i,L) = 0.
618     enddo
619     enddo
620    
621     c zero out fluxes and derivatives
622     c -------------------------------
623     do i = 1,istrip
624     eturb(i) = 0.
625     scu(i) = 0.
626     dedqa(i) = 0.
627     dedtc(i) = 0.
628     hsturb(i) = 0.
629     dhsdqa(i) = 0.
630     dhsdtc(i) = 0.
631     enddo
632    
633     c increment diagnostic arrays for quantities calculated before trbfl
634     c ------------------------------------------------------------------
635 molod 1.23 c do i =1,istrip
636     c stdiag(i) = ( thv(i,nlay+1)-thv(i,nlay) ) / pke(i,nlay+1)
637     c enddo
638     c if(idtsrf.gt.0) then
639     c call paste2grd(stdiag,igrd,chfrstr,istrip,nchp,
640     c 1 qdiag(1,1,idtsrf,bi,bj),ijall,1,nn,.false.)
641     c endif
642 molod 1.1
643 molod 1.20
644 molod 1.24 if(2.eq.1)then
645 molod 1.21 print *,' In turb before trbflx - strip ',nn,' out of ',numstrips
646 molod 1.23 print *,' bi = ',bi
647 molod 1.20 print *,' ntracers ',ntracers,' ptracers ',ptracers
648     print *,'dttrb,itrtrb,rmu,edle ',dttrb,' ',itrtrb,' ',rmu,' ',edle
649 molod 1.21 print *,' nchp ',nchp,' nchptot ',nchptot,' nchplnd ',nchplnd
650 molod 1.20 print *,' qbeg, tprof ',qbeg,' ',tprof
651 molod 1.21 print *,'istrip,nlay,nymd,nhms ',istrip,' ',nlay,' ',nymd,' ',nhms
652     print *,' grav,cp,rgas,faceps,virtcon,undef ',
653     . grav,' ',cp,' ',rgas,' ',faceps,' ',virtcon,' ',undef
654     print *,' field: th ',th
655 molod 1.23 c print *,' field: thv ',thv
656     c print *,' field: sh ',sh
657     c print *,' field: u ',u
658     c print *,' field: v ',v
659 molod 1.21 print *,' field: p ',p
660 molod 1.23 c print *,' field: pe ',pe
661     c print *,' field: pk ',pk
662     c print *,' field: pke ',pke
663     c print *,' field: dpstr ',dpstr
664     c print *,' field: stwatr ',stwatr
665     c print *,' field: stz0 ',stz0
666 molod 1.20 endif
667    
668 molod 1.1 c call trbflx
669     c -----------
670     call trbflx(nn,th,thv,sh,u,v,qq,p,pe,pk,pke,dpstr,stwatr,stz0,
671     1 tracers,ntracers-ptracers,ntracedim,dttrb,itrtrb,rmu,edle,qbeg,
672     2 tprof,stuflux,stvflux,sri,skh,skm,swinds,sustar,sz0,frqtrb,
673     3 pbldpth,sct,scu,stu2m,stv2m,stt2m,stq2m,stu10m,stv10m,stt10m,
674     4 stq10m,istrip,nlay,nymd,nhms,grav,cp,rgas,faceps,virtcon,undef,
675     5 dshdthg,dshdshg,dthdthg,dthdshg,eturb,dedqa,dedtc,
676     6 hsturb,dhsdqa,dhsdtc,transth,transsh,
677 molod 1.11 7 ctsave,xxsave,yysave,zetasave,xlsave,khsave,qliq,turbfcc)
678 molod 1.1
679 molod 1.23
680 molod 1.20 if(2.eq.1)then
681 molod 1.23 print *,' In turbio, Just after trbflx for strip ',nn,' bi = ',bi
682 molod 1.21 print *,' field: stuflux ',stuflux
683     print *,' field: stvflux ',stvflux
684     print *,' field: dshdthg ',dshdthg
685     print *,' field: dshdshg ',dshdshg
686     print *,' field: dthdthg ',dthdthg
687     print *,' field: dthdshg ',dthdshg
688     print *,' field: scu ',scu
689     print *,' field: eturb ',eturb
690     print *,' field: dedqa ',dedqa
691     print *,' field: dedtc ',dedtc
692     print *,' field: hsturb ',hsturb
693     print *,' field: dhsdqa ',dhsdqa
694     print *,' field: dhsdtc ',dhsdtc
695     print *,' field: transth ',transth
696     print *,' field: transsh ',transsh
697 molod 1.20 endif
698    
699 molod 1.13 call pastit (qq,tke,istrip,nchp,nchptot,nlay,nn)
700     call pastit (ctsave,ctmt,istrip,nchp,nchptot,1,nn)
701     call pastit (xxsave,xxmt,istrip,nchp,nchptot,1,nn)
702     call pastit (yysave,yymt,istrip,nchp,nchptot,1,nn)
703     call pastit (zetasave,zetamt,istrip,nchp,nchptot,1,nn)
704     call pastit (xlsave,xlmt,istrip,nchp,nchptot,nlay,nn)
705     call pastit (khsave,khmt,istrip,nchp,nchptot,nlay,nn)
706 molod 1.1
707 molod 1.13 call pastit (qliq ,qliqmsc,istrip,nchp,nchptot,nlay,nn)
708     call pastit (turbfcc,fccmsc,istrip,nchp,nchptot,nlay,nn)
709 molod 1.1
710     c New diagnostic: potential evapotranspiration
711     do i = 1,istrip
712     evpot(i) = transsh(i,nlay) * (shgstr(i) - sh(i,nlay))
713     enddo
714    
715     C**********************************************************************
716     C Call Land Surface Module
717     C**********************************************************************
718    
719     do i = 1,istrip
720     savetc(i) = tc(i)
721     saveqa(i) = qa(i)
722     enddo
723     do i = 1,istrip
724     cosz(i) = max(cosz(i),0.0001)
725     cd(i) = scu(i)*scu(i)
726     tmpnlay(i) = th(i,nlay)*pk(i,nlay)
727     hlwdwn(i) = alwrad(i)+blwrad(i)*tc(i)-lwstrip(i)
728     psurf(i) = pe(i,nlay+1)
729     wspeed(i) = sqrt(u(i,nlay)*u(i,nlay) + v(i,nlay)*v(i,nlay))
730 molod 1.21 if(wspeed(i) .lt. 1.e-20) wspeed(i) = 1.e-20
731 molod 1.1 C Note: This LSM precip bug needs to be cleaned up
732     ccc newsnow(i) = newsnow(i)*dtfac
733     ccc raincon(i) = raincon(i)*dtfac
734     ccc rainls (i) = rainls (i)*dtfac
735     enddo
736    
737     do i = 1,istrip
738     eturb(i) = eturb(i) * pke(i,nlay+1)
739     dedqa(i) = dedqa(i) * pke(i,nlay+1)
740     hsturb(i) = hsturb(i) * pke(i,nlay+1)
741     enddo
742    
743     do i = 1,istrip
744     strdg1(i) = 0.
745     strdg2(i) = 0.
746     strdg3(i) = 0.
747     strdg4(i) = 0.
748     strdg5(i) = 0.
749     strdg6(i) = 0.
750     strdg7(i) = 0.
751     strdg8(i) = 0.
752     strdg9(i) = 0.
753     bomb(i) = 0.
754     runoff(i) = 0.
755     eint(i) = 0.
756     esoi(i) = 0.
757     eveg(i) = 0.
758     esno(i) = 0.
759     smelt(i) = 0.
760     hlatn(i) = 0.
761     hlwup(i) = 0.
762     gdrain(i) = 0.
763     runsrf(i) = 0.
764     fwsoil(i) = 0.
765     enddo
766    
767     c**********************************************************************
768     c diagnostics: fill arrays for lsm input fields
769     c**********************************************************************
770 molod 1.23 c if(isnowfall.gt.0) then
771     c do i = 1,istrip
772     c tmpstrip(i) = newsnow(i)*86400
773     c enddo
774     c call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,
775     c 1 qdiag(1,1,isnowfall,bi,bj),ijall,1,nn,.false.)
776     c endif
777 molod 1.29 if(diagnostics_is_on('RAINCON ',myid) ) then
778     do i = 1,istrip
779     tmpstrip(i) = raincon(i)*86400
780     enddo
781     call diag_vegtile_fill (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
782     . .false., 'RAINCON ', 1, 1, bi, bj, myid)
783     endif
784     if(diagnostics_is_on('RAINLSP ',myid) ) then
785     do i = 1,istrip
786     tmpstrip(i) = rainls(i)*86400
787     enddo
788     call diag_vegtile_fill (tmpstrip,igrd,chfrstr,istrip,nchp,nn,
789     . .false., 'RAINLSP ', 1, 1, bi, bj, myid)
790     endif
791 molod 1.23 c if(igreen.gt.0) then
792     c call paste2grd(grnstrip,igrd,chfrstr,istrip,nchp,
793     c 1 qdiag(1,1,igreen,bi,bj),ijall,1,nn,.false.)
794     c endif
795     c if(ilai.gt.0) then
796     c call paste2grd(laistrip,igrd,chfrstr,istrip,nchp,
797     c 1 qdiag(1,1,ilai,bi,bj),ijall,1,nn,.false.)
798     c endif
799 molod 1.29 if(diagnostics_is_on('PARDR ',myid) ) then
800     call diag_vegtile_fill (pardr,igrd,chfrstr,istrip,nchp,nn,
801     . .false., 'PARDR ', 1, 1, bi, bj, myid)
802     endif
803 molod 1.23 c if(ipardf.gt.0) then
804     c call paste2grd(pardf,igrd,chfrstr,istrip,nchp,
805     c 1 qdiag(1,1,ipardf,bi,bj),ijall,1,nn,.false.)
806     c endif
807     c if(idlwdtc.gt.0) then
808     c call paste2grd(blwrad,igrd,chfrstr,istrip,nchp,
809     c 1 qdiag(1,1,idlwdtc,bi,bj),ijall,1,nn,.false.)
810     c endif
811     c if(idhdtc.gt.0) then
812     c call paste2grd(dhsdtc,igrd,chfrstr,istrip,nchp,
813     c 1 qdiag(1,1,idhdtc,bi,bj),ijall,1,nn,.false.)
814     c endif
815     c if(idedtc.gt.0) then
816     c call paste2grd(dedtc,igrd,chfrstr,istrip,nchp,
817     c 1 qdiag(1,1,idedtc,bi,bj),ijall,1,nn,.false.)
818     c endif
819     c if(idhdqa.gt.0) then
820     c call paste2grd(dhsdqa,igrd,chfrstr,istrip,nchp,
821     c 1 qdiag(1,1,idhdqa,bi,bj),ijall,1,nn,.false.)
822     c endif
823     c if(idedqa.gt.0) then
824     c call paste2grd(dedqa,igrd,chfrstr,istrip,nchp,
825     c 1 qdiag(1,1,idedqa,bi,bj),ijall,1,nn,.false.)
826     c endif
827     c if(ilwgdown.gt.0) then
828     c call paste2grd(hlwdwn,igrd,chfrstr,istrip,nchp,
829     c 1 qdiag(1,1,ilwgdown,bi,bj),ijall,1,nn,.false.)
830     c endif
831 molod 1.1 c**********************************************************************
832    
833 molod 1.21 if(nland.gt.0)then
834    
835 molod 1.24 if(2.eq.1)then
836 molod 1.23 print *,' In turbio, Just before tile for strip ',nn,' bi = ',bi
837 molod 1.21 print *,' calling tile for ',nland,' land points '
838     print *,' field: types ',types
839 molod 1.23 print *,' field: chfrstr ',chfrstr
840     c print *,' field: rainls ',rainls
841     c print *,' field: newsnow ',newsnow
842     c print *,' field: wspeed ',wspeed
843 molod 1.26 print *,' field: eturb ',eturb
844     print *,' field: dedqa ',dedqa
845     print *,' field: dedtc ',dedtc
846     print *,' field: hsturb ',hsturb
847     print *,' field: dhsdqa ',dhsdqa
848     print *,' field: dhsdtc ',dhsdtc
849 molod 1.23 c print *,' field: tmpnlay ',tmpnlay
850     c print *,' field: sh(nlay) ',(sh(i,nlay),i=1,istrip)
851     c print *,' field: cd ',cd
852     c print *,' field: cosz ',cosz
853     c print *,' field: pardr ',pardr
854     c print *,' field: pardf ',pardf
855 molod 1.26 print *,' field: swnet ',swnet
856     print *,' field: hlwdwn ',hlwdwn
857 molod 1.23 c print *,' field: psurf ',psurf
858     c print *,' field: laistrip ',laistrip
859     c print *,' field: grnstrip ',grnstrip
860     c print *,' field: z2str ',z2str
861     c print *,' field: scatstr ',scatstr
862     c print *,' field: z2str ',z2str
863     c print *,' field: rs1str ',rs1str
864     c print *,' field: rs1str ',rs2str
865     c print *,' field: rdcstr ',rdcstr
866     c print *,' field: u2fstr ',u2fstr
867     c print *,' field: dqsdtstr ',dqsdtstr
868 molod 1.26 print *,' field: alwrad ',alwrad
869     print *,' field: blwrad ',blwrad
870 molod 1.21 print *,' field: tc ',tc
871 molod 1.26 print *,' field: td ',td
872 molod 1.23 c print *,' field: swet1 ',swet1
873     c print *,' field: swet2 ',swet2
874     c print *,' field: swet3 ',swet3
875     c print *,' field: capacity ',capacity
876     c print *,' field: snowdepth ',snowdepth
877     endif
878    
879 molod 1.21 call tile (
880     I nland, timstp, types, rainls, raincon, newsnow, wspeed,
881     I eturb, dedqa, dedtc, hsturb, dhsdqa, dhsdtc,
882     I tmpnlay, sh(1,nlay), cd, cosz, pardr, pardf,
883     I swnet, hlwdwn, psurf, laistrip, grnstrip, z2str,
884     I scatstr, rs1str, rs2str, rdcstr, u2fstr,
885     I shgstr, dqsdtstr, alwrad, blwrad,
886     U tc, td, qa, swet1, swet2, swet3, capacity, snowdepth,
887     O evap, shflux, runoff, bomb,
888     O eint, esoi, eveg, esno, smelt, hlatn,
889     O hlwup, gdrain, runsrf, fwsoil,
890     O strdg1, strdg2, strdg3, strdg4,
891     O strdg5, strdg6, strdg7, strdg8, strdg9)
892     endif
893 molod 1.20
894 molod 1.21 if(2.eq.1)then
895     print *,' In turbio, Just after tile for strip ',nn
896     print *,' calling tile for ',nland,' land points '
897     print *,' field: tc ',tc
898 molod 1.26 c print *,' field: td ',td
899 molod 1.21 print *,' field: strdg1 ',strdg1
900     print *,' field: strdg2 ',strdg2
901     print *,' field: strdg3 ',strdg3
902     print *,' field: strdg4 ',strdg4
903     print *,' field: strdg5 ',strdg5
904     print *,' field: strdg6 ',strdg6
905     print *,' field: strdg7 ',strdg7
906     print *,' field: strdg8 ',strdg8
907     print *,' field: strdg9 ',strdg9
908 molod 1.26 c print *,' field: swet1 ',swet1
909     c print *,' field: swet2 ',swet2
910     c print *,' field: swet3 ',swet3
911     c print *,' field: capacity ',capacity
912     c print *,' field: snowdepth ',snowdepth
913 molod 1.21 endif
914     if( nice.gt.0 ) then
915     print *,' Calling sea ice routine - SHOULD NOT BE HERE!'
916     call seaice ( nocean, timstp, hice,
917     . eturb(nland+1), dedtc(nland+1),
918     . hsturb(nland+1), dhsdtc(nland+1),
919     . qicestr(nland+1), dqicestr(nland+1),
920     . swnet(nland+1), lwstrip(nland+1), blwrad(nland+1),
921     . pke(nland+1,nlay+1), icest(nland+1),
922     . tc(nland+1), qa(nland+1) )
923     endif
924 molod 1.1
925     c***********************************************************************
926     c diagnostics: fill arrays for lsm output fields
927     c***********************************************************************
928    
929 molod 1.23 c if(irunoff.gt.0) then
930     c call paste2grd(runoff,igrd,chfrstr,istrip,nchp,
931     c 1 qdiag(1,1,irunoff,bi,bj),ijall,1,nn,.false.)
932     c endif
933     c if(ifwsoil.gt.0) then
934     c call paste2grd(fwsoil,igrd,chfrstr,istrip,nchp,
935     c 1 qdiag(1,1,ifwsoil,bi,bj),ijall,1,nn,.false.)
936     c endif
937     c if(igdrain.gt.0) then
938     c call paste2grd(gdrain,igrd,chfrstr,istrip,nchp,
939     c 1 qdiag(1,1,igdrain,bi,bj),ijall,1,nn,.false.)
940     c endif
941     c if(isnowmelt.gt.0) then
942     c call paste2grd(smelt,igrd,chfrstr,istrip,nchp,
943     c 1 qdiag(1,1,isnowmelt,bi,bj),ijall,1,nn,.false.)
944     c endif
945     c if(ieveg.gt.0) then
946     c call paste2grd(eveg,igrd,chfrstr,istrip,nchp,
947     c 1 qdiag(1,1,ieveg,bi,bj),ijall,1,nn,.false.)
948     c endif
949     c if(iesnow.gt.0) then
950     c call paste2grd(esno,igrd,chfrstr,istrip,nchp,
951     c 1 qdiag(1,1,iesnow,bi,bj),ijall,1,nn,.false.)
952     c endif
953     c if(iesoil.gt.0) then
954     c call paste2grd(esoi,igrd,chfrstr,istrip,nchp,
955     c 1 qdiag(1,1,iesoil,bi,bj),ijall,1,nn,.false.)
956     c endif
957     c if(ieresv.gt.0) then
958     c call paste2grd(eint,igrd,chfrstr,istrip,nchp,
959     c 1 qdiag(1,1,ieresv,bi,bj),ijall,1,nn,.false.)
960     c endif
961     c if(ievpot.gt.0) then
962     c call paste2grd(evpot,igrd,chfrstr,istrip,nchp,
963     c 1 qdiag(1,1,ievpot,bi,bj),ijall,1,nn,.false.)
964     c endif
965     c if(idtc.gt.0) then
966     c call paste2grd(strdg1,igrd,chfrstr,istrip,nchp,
967     c 1 qdiag(1,1,idtc,bi,bj),ijall,1,nn,.false.)
968     c endif
969     c if(idqc.gt.0) then
970     c call paste2grd(strdg2,igrd,chfrstr,istrip,nchp,
971     c 1 qdiag(1,1,idqc,bi,bj),ijall,1,nn,.false.)
972     c endif
973     c if(itcdtc.gt.0) then
974     c call paste2grd(strdg3,igrd,chfrstr,istrip,nchp,
975     c 1 qdiag(1,1,itcdtc,bi,bj),ijall,1,nn,.false.)
976     c endif
977     c if(iraddtc.gt.0) then
978     c call paste2grd(strdg4,igrd,chfrstr,istrip,nchp,
979     c 1 qdiag(1,1,iraddtc,bi,bj),ijall,1,nn,.false.)
980     c endif
981     c if(isensdtc.gt.0) then
982     c call paste2grd(strdg5,igrd,chfrstr,istrip,nchp,
983     c 1 qdiag(1,1,isensdtc,bi,bj),ijall,1,nn,.false.)
984     c endif
985     c if(ilatdtc.gt.0) then
986     c call paste2grd(strdg6,igrd,chfrstr,istrip,nchp,
987     c 1 qdiag(1,1,ilatdtc,bi,bj),ijall,1,nn,.false.)
988     c endif
989     c if(itddtc.gt.0) then
990     c call paste2grd(strdg7,igrd,chfrstr,istrip,nchp,
991     c 1 qdiag(1,1,itddtc,bi,bj),ijall,1,nn,.false.)
992     c endif
993     c if(iqcdtc.gt.0) then
994     c call paste2grd(strdg8,igrd,chfrstr,istrip,nchp,
995     c 1 qdiag(1,1,iqcdtc,bi,bj),ijall,1,nn,.false.)
996     c endif
997 molod 1.6 c***********************************************************************
998 molod 1.1
999     if( ndlsm.gt.1 ) then
1000 molod 1.13 call pstbitint(types,qdiaglsm(1,1),istrip,nchp,nchptot,1,nn)
1001     call pstbmpit(chfrstr,qdiaglsm(1,2),istrip,nchp,nchptot,1,nn)
1002     call pstbmpit(lats,qdiaglsm(1,3),istrip,nchp,nchptot,1,nn)
1003     call pstbmpit(lons,qdiaglsm(1,4),istrip,nchp,nchptot,1,nn)
1004     c call pstbmpit(igrdstr,qdiaglsm(1,5),istrip,nchp,nchptot,1,nn)
1005     call pstbmpit(tc,qdiaglsm(1,6),istrip,nchp,nchptot,1,nn)
1006     call pstbmpit(td,qdiaglsm(1,7),istrip,nchp,nchptot,1,nn)
1007     call pstbmpit(qa,qdiaglsm(1,8),istrip,nchp,nchptot,1,nn)
1008     call pstbmpit(swet1,qdiaglsm(1,9),istrip,nchp,nchptot,1,nn)
1009     call pstbmpit(swet2,qdiaglsm(1,10),istrip,nchp,nchptot,1,nn)
1010     call pstbmpit(swet3,qdiaglsm(1,11),istrip,nchp,nchptot,1,nn)
1011     call pstbmpit(capacity,qdiaglsm(1,12),istrip,nchp,nchptot,1,nn)
1012     call pstbmpit(snowdepth,qdiaglsm(1,13),istrip,nchp,nchptot,1,nn)
1013     call pstbmpit(eturb,qdiaglsm(1,14),istrip,nchp,nchptot,1,nn)
1014     call pstbmpit(hsturb,qdiaglsm(1,15),istrip,nchp,nchptot,1,nn)
1015     call pstbmpit(cd,qdiaglsm(1,16),istrip,nchp,nchptot,1,nn)
1016     call pstbmpit(laistrip,qdiaglsm(1,17),istrip,nchp,nchptot,1,nn)
1017     call pstbmpit(grnstrip,qdiaglsm(1,18),istrip,nchp,nchptot,1,nn)
1018     call pstbmpit(eint,qdiaglsm(1,19),istrip,nchp,nchptot,1,nn)
1019     call pstbmpit(esoi,qdiaglsm(1,20),istrip,nchp,nchptot,1,nn)
1020     call pstbmpit(eveg,qdiaglsm(1,21),istrip,nchp,nchptot,1,nn)
1021     call pstbmpit(esno,qdiaglsm(1,22),istrip,nchp,nchptot,1,nn)
1022     call pstbmpit(strdg1,qdiaglsm(1,23),istrip,nchp,nchptot,1,nn)
1023     call pstbmpit(strdg2,qdiaglsm(1,24),istrip,nchp,nchptot,1,nn)
1024     call pstbmpit(strdg3,qdiaglsm(1,25),istrip,nchp,nchptot,1,nn)
1025     call pstbmpit(strdg4,qdiaglsm(1,26),istrip,nchp,nchptot,1,nn)
1026     call pstbmpit(strdg5,qdiaglsm(1,27),istrip,nchp,nchptot,1,nn)
1027     call pstbmpit(strdg6,qdiaglsm(1,28),istrip,nchp,nchptot,1,nn)
1028     call pstbmpit(strdg7,qdiaglsm(1,29),istrip,nchp,nchptot,1,nn)
1029     call pstbmpit(strdg8,qdiaglsm(1,30),istrip,nchp,nchptot,1,nn)
1030     call pstbmpit(strdg9,qdiaglsm(1,31),istrip,nchp,nchptot,1,nn)
1031     call pstbmpit(smelt,qdiaglsm(1,32),istrip,nchp,nchptot,1,nn)
1032     call pstbmpit(gdrain,qdiaglsm(1,33),istrip,nchp,nchptot,1,nn)
1033     call pstbmpit(runsrf,qdiaglsm(1,34),istrip,nchp,nchptot,1,nn)
1034     call pstbmpit(fwsoil,qdiaglsm(1,35),istrip,nchp,nchptot,1,nn)
1035     call pstbmpit(evpot,qdiaglsm(1,36),istrip,nchp,nchptot,1,nn)
1036     call pstbmpit(stt2m,qdiaglsm(1,37),istrip,nchp,nchptot,1,nn)
1037     call pstbmpit(stq2m,qdiaglsm(1,38),istrip,nchp,nchptot,1,nn)
1038     endif
1039    
1040     call pastit (tc,tground,istrip,nchp,nchptot,1,nn)
1041     call pastit (td,tdeep,istrip,nchp,nchptot,1,nn)
1042     call pastit (qa,qground,istrip,nchp,nchptot,1,nn)
1043     call pastit (swet1,swetshal,istrip,nchp,nchptot,1,nn)
1044     call pastit (swet2,swetroot,istrip,nchp,nchptot,1,nn)
1045     call pastit (swet3,swetdeep,istrip,nchp,nchptot,1,nn)
1046     call pastit (capacity,capac,istrip,nchp,nchptot,1,nn)
1047     call pastit (snowdepth,snodep,istrip,nchp,nchptot,1,nn)
1048 molod 1.1
1049     c**********************************************************************
1050     c Now update the theta and sh profiles with the new ground temperature
1051     c**********************************************************************
1052    
1053     do i =1,istrip
1054     th(i,nlay+1) = tc(i) / pke(i,nlay+1)
1055     enddo
1056     do L = 1,nlay
1057     do i =1,istrip
1058     th(i,L) = th(i,L) + dthdthg(i,L)*(tc(i)-savetc(i))/pke(i,nlay+1)
1059     enddo
1060     enddo
1061    
1062     do i =1,istrip
1063     sh(i,nlay+1) = qa(i)
1064     enddo
1065     do L = 1,nlay
1066     do i =1,istrip
1067     sh(i,L) = sh(i,L) + dshdshg(i,L)*(qa(i)-saveqa(i))
1068     enddo
1069     enddo
1070    
1071     do L = 1,nlay
1072     do i =1,istrip
1073     sttflux(i,L) = transth(i,L) * (th(i,L+1)-th(i,L))
1074     stqflux(i,L) = transsh(i,L) * (sh(i,L+1)-sh(i,L))
1075     enddo
1076     enddo
1077    
1078 molod 1.21 if(2.eq.1)then
1079     print *,' In turbio, just after updating th and sh - strip ',nn
1080     print *,' field: th ',th
1081     print *,' field: sh ',sh
1082     endif
1083    
1084 molod 1.1 c tendency updates
1085     c ----------------
1086     do l=1,nlay
1087     call strip2tile(uz(1,1,l),igrd,tmpstrip,nchp,ijall,
1088     1 istrip,1,nn)
1089     do i =1,istrip
1090     tends(i) = ( u(i,l)-tmpstrip(i) )
1091     enddo
1092 molod 1.13 call pastit (tends,dumsc(1,l),istrip,nchp,nchptot,1,nn)
1093 molod 1.1
1094     call strip2tile(vz(1,1,l),igrd,tmpstrip,nchp,ijall,
1095     1 istrip,1,nn)
1096     do i =1,istrip
1097     tends(i) = ( v(i,l)-tmpstrip(i) )
1098     enddo
1099 molod 1.13 call pastit (tends,dvmsc(1,l),istrip,nchp,nchptot,1,nn)
1100 molod 1.1
1101     call strip2tile(tz(1,1,l),igrd,tmpstrip,nchp,ijall,
1102     1 istrip,1,nn)
1103     do i =1,istrip
1104     tends(i) = ( th(i,l)-tmpstrip(i) )
1105     enddo
1106 molod 1.21
1107 molod 1.22 if(2.eq.1)then
1108 molod 1.21 print *,' In turbio, tendencies for strip ',nn,' level ',l
1109     print *,' field: th ',tends
1110     endif
1111    
1112 molod 1.13 call pastit (tends,dtmsc(1,l),istrip,nchp,nchptot,1,nn)
1113 molod 1.1
1114     call strip2tile(qz(1,1,l,1),igrd,tmpstrip,nchp,ijall,
1115     1 istrip,1,nn)
1116     do i =1,istrip
1117     tends(i) = ( sh(i,l)-tmpstrip(i) )
1118     enddo
1119 molod 1.21
1120 molod 1.22 if(2.eq.1)then
1121 molod 1.21 print *,' In turbio, tendencies for strip ',nn,' level ',l
1122     print *,' field: sh ',tends
1123     endif
1124    
1125 molod 1.13 call pastit (tends,dqmsc(1,l,1),istrip,nchp,nchptot,1,nn)
1126 molod 1.1
1127 molod 1.20 c do nt = 1,ntracers-ptracers
1128     c call strip2tile(qz(1,1,L,ptracers+nt),igrd,tmpstrip,nchp,
1129     c 1 ijall,istrip,1,nn)
1130     c do i =1,istrip
1131     c tends(i) = ( tracers(i,L,nt)-tmpstrip(i) )
1132     c enddo
1133     c call pastit(tends,dqmsc(1,L,ptracers+nt),istrip,nchp,
1134     c . nchptot,1,nn)
1135     c enddo
1136 molod 1.1
1137     enddo
1138    
1139     c *********************************************************************
1140     c **** increment diagnostic arrays for quantities saved in trbflx
1141     c *********************************************************************
1142    
1143     c note: the order, logic, and scaling of the heat and moisture flux
1144     c diagnostics is critical!
1145     c ------------------------------
1146    
1147 molod 1.23 c if(ievap.gt.0) then
1148     c do i=1,istrip
1149     c tmpstrip(i) = stqflux(i,nlay) * 86400
1150     c enddo
1151     c call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,
1152     c 1 qdiag(1,1,ievap,bi,bj),ijall,1,nn,.false.)
1153     c endif
1154     c if(ieflux.gt.0) then
1155     c do i=1,istrip
1156     c tmpstrip(i) = stqflux(i,nlay) * alhl
1157     c enddo
1158     c call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,
1159     c 1 qdiag(1,1,ieflux,bi,bj),ijall,1,nn,.false.)
1160     c endif
1161     c if(ihflux.gt.0) then
1162     c do i=1,istrip
1163     c tmpstrip(i) = sttflux(i,nlay) * cp
1164     c enddo
1165     c call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,
1166     c 1 qdiag(1,1,ihflux,bi,bj),ijall,1,nn,.false.)
1167     c endif
1168 molod 1.30 if(diagnostics_is_on('TUFLUX ',myid) ) then
1169     call diag_vegtile_fill (stuflux,igrd,chfrstr,istrip,nchp,nn,
1170     . .false., 'TUFLUX ', 0, nlay, bi, bj, myid)
1171     endif
1172     if(diagnostics_is_on('TVFLUX ',myid) ) then
1173     call diag_vegtile_fill (stvflux,igrd,chfrstr,istrip,nchp,nn,
1174     . .false., 'TVFLUX ', 0, nlay, bi, bj, myid)
1175     endif
1176 molod 1.23 c if(ituflux.gt.0) then
1177     c call paste2grd(stuflux,igrd,chfrstr,istrip,nchp,
1178     c 1 qdiag(1,1,ituflux,bi,bj),ijall,nlay,nn,.false.)
1179     c endif
1180 molod 1.30 if(diagnostics_is_on('TTFLUX ',myid) ) then
1181     do l=1,nlay
1182     do i=1,istrip
1183     sttflux(i,l) = sttflux(i,l) * cp
1184     enddo
1185     enddo
1186     call diag_vegtile_fill (sttflux,igrd,chfrstr,istrip,nchp,nn,
1187     . .false., 'TTFLUX ', 0, nlay, bi, bj, myid)
1188     endif
1189     if(diagnostics_is_on('TQFLUX ',myid) ) then
1190     do l=1,nlay
1191     do i=1,istrip
1192     stqflux(i,l) = stqflux(i,l) * alhl
1193     enddo
1194     enddo
1195     call diag_vegtile_fill (stqflux,igrd,chfrstr,istrip,nchp,nn,
1196     . .false., 'TQFLUX ', 0, nlay, bi, bj, myid)
1197     endif
1198 molod 1.29 if(diagnostics_is_on('RI ',myid) ) then
1199     call diag_vegtile_fill (sri,igrd,chfrstr,istrip,nchp,nn,
1200     . .false., 'RI ', 0, nlay, bi, bj, myid)
1201     endif
1202     if(diagnostics_is_on('KH ',myid) ) then
1203     call diag_vegtile_fill (skh,igrd,chfrstr,istrip,nchp,nn,
1204     . .false., 'KH ', 0, nlay, bi, bj, myid)
1205     endif
1206     if(diagnostics_is_on('KM ',myid) ) then
1207     call diag_vegtile_fill (skm,igrd,chfrstr,istrip,nchp,nn,
1208     . .false., 'KM ', 0, nlay, bi, bj, myid)
1209     endif
1210 molod 1.23 c if(iri.gt.0) call paste2grd(sri,igrd,chfrstr,istrip,nchp,
1211     c 1 qdiag(1,1,iri,bi,bj),ijall,nlay,nn,.false.)
1212     c if(ikh.gt.0) call paste2grd(skh,igrd,chfrstr,istrip,nchp,
1213     c 1 qdiag(1,1,ikh,bi,bj),ijall,nlay,nn,.false.)
1214     c if(ikm.gt.0) call paste2grd(skm,igrd,chfrstr,istrip,nchp,
1215     c 1 qdiag(1,1,ikm,bi,bj),ijall,nlay,nn,.false.)
1216     c if(ict.gt.0) then
1217     c call paste2grd(sct,igrd,chfrstr,istrip,nchp,
1218     c 1 qdiag(1,1,ict,bi,bj),ijall,1,nn,.false.)
1219     c endif
1220     c if(icu.gt.0) then
1221     c call paste2grd(scu,igrd,chfrstr,istrip,nchp,
1222     c 1 qdiag(1,1,icu,bi,bj),ijall,1,nn,.false.)
1223     c endif
1224     c if(iwinds.gt.0) then
1225     c call paste2grd(swinds,igrd,chfrstr,istrip,nchp,
1226     c 1 qdiag(1,1,iwinds,bi,bj),ijall,1,nn,.false.)
1227     c endif
1228     c if(iuflux.gt.0) then
1229     c call paste2grd(stuflux(1,nlay),igrd,chfrstr,istrip,nchp,
1230     c 1 qdiag(1,1,iuflux,bi,bj),ijall,1,nn,.false.)
1231     c endif
1232     c if(ivflux.gt.0) then
1233     c call paste2grd(stvflux(1,nlay),igrd,chfrstr,istrip,nchp,
1234     c 1 qdiag(1,1,ivflux,bi,bj),ijall,1,nn,.false.)
1235     c endif
1236     c if(iustar.gt.0) then
1237     c call paste2grd(sustar,igrd,chfrstr,istrip,nchp,
1238     c 1 qdiag(1,1,iustar,bi,bj),ijall,1,nn,.false.)
1239     c endif
1240     c if(iz0.gt.0) then
1241     c call paste2grd(sz0,igrd,chfrstr,istrip,nchp,
1242     c 1 qdiag(1,1,iz0,bi,bj),ijall,1,nn,.false.)
1243     c endif
1244     c if(ifrqtrb.gt.0) then
1245     c call paste2grd(frqtrb,igrd,chfrstr,istrip,nchp,
1246     c 1 qdiag(1,1,ifrqtrb,bi,bj),ijall,1,nn,.false.)
1247     c endif
1248     c if(ipbl.gt.0) then
1249     c call paste2grd(pbldpth,igrd,chfrstr,istrip,nchp,
1250     c 1 qdiag(1,1,ipbl,bi,bj),ijall,1,nn,.false.)
1251     c endif
1252     c if(iu2m.gt.0) then
1253     c call paste2grd(stu2m,igrd,chfrstr,istrip,nchp,
1254     c 1 qdiag(1,1,iu2m,bi,bj),ijall,1,nn,.true.)
1255     c endif
1256     c if(iv2m.gt.0) then
1257     c call paste2grd(stv2m,igrd,chfrstr,istrip,nchp,
1258     c 1 qdiag(1,1,iv2m,bi,bj),ijall,1,nn,.true.)
1259     c endif
1260     c if(it2m.gt.0) then
1261     c call paste2grd(stt2m,igrd,chfrstr,istrip,nchp,
1262     c 1 qdiag(1,1,it2m,bi,bj),ijall,1,nn,.true.)
1263     c endif
1264     c if(iq2m.gt.0) then
1265     c do i=1,istrip
1266     c if( stq2m(i).ne.undef ) then
1267     c tmpstrip(i) = stq2m(i) * 1000
1268     c else
1269     c tmpstrip(i) = undef
1270     c endif
1271     c enddo
1272     c call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,
1273     c 1 qdiag(1,1,iq2m,bi,bj),ijall,1,nn,.true.)
1274     c endif
1275     c if(iu10m.gt.0) then
1276     c call paste2grd(stu10m,igrd,chfrstr,istrip,nchp,
1277     c 1 qdiag(1,1,iu10m,bi,bj),ijall,1,nn,.false.)
1278     c endif
1279     c if(iv10m.gt.0) then
1280     c call paste2grd(stv10m,igrd,chfrstr,istrip,nchp,
1281     c 1 qdiag(1,1,iv10m,bi,bj),ijall,1,nn,.false.)
1282     c endif
1283     c if(it10m.gt.0) then
1284     c call paste2grd(stt10m,igrd,chfrstr,istrip,nchp,
1285     c 1 qdiag(1,1,it10m,bi,bj),ijall,1,nn,.false.)
1286     c endif
1287     c if(iq10m.gt.0) then
1288     c do i=1,istrip
1289     c if( stq10m(i).ne.undef ) then
1290     c tmpstrip(i) = stq10m(i) * 1000
1291     c else
1292     c tmpstrip(i) = undef
1293     c endif
1294     c enddo
1295     c call paste2grd(tmpstrip,igrd,chfrstr,istrip,nchp,
1296     c 1 qdiag(1,1,iq10m,bi,bj),ijall,1,nn,.false.)
1297     c endif
1298 molod 1.1
1299     c**********************************************************************
1300     c more diagnostics: land surface model parameters
1301     c**********************************************************************
1302    
1303 molod 1.23 c if(itdeep.gt.0)call paste2grd(td,igrd,chfrstr,istrip,nchp,
1304     c . qdiag(1,1,itdeep,bi,bj),ijall,1,nn,.false.)
1305     c if(iqcanopy .gt.0)call paste2grd(qa,igrd,chfrstr,istrip,nchp,
1306     c . qdiag(1,1,iqcanopy,bi,bj) ,ijall,1,nn,.false.)
1307     c if(ismshal .gt.0)call paste2grd(swet1,igrd,chfrstr,istrip,nchp,
1308     c . qdiag(1,1,ismshal,bi,bj) ,ijall,1,nn,.false.)
1309     c if(ismroot .gt.0)call paste2grd(swet2,igrd,chfrstr,istrip,nchp,
1310     c . qdiag(1,1,ismroot,bi,bj) ,ijall,1,nn,.false.)
1311     c if(ismdeep .gt.0)call paste2grd(swet3,igrd,chfrstr,istrip,nchp,
1312     c . qdiag(1,1,ismdeep,bi,bj) ,ijall,1,nn,.false.)
1313     c if(icapacity.gt.0)call paste2grd(capacity,igrd,chfrstr,istrip,
1314     c . nchp,qdiag(1,1,icapacity,bi,bj),ijall,1,nn,.false.)
1315     c if(isnow.gt.0)call paste2grd(snowdepth,igrd,chfrstr,istrip,nchp,
1316     c . qdiag(1,1,isnow,bi,bj) ,ijall,1,nn,.false.)
1317 molod 1.1
1318     c**********************************************************************
1319     c end regions loop
1320    
1321     2000 continue
1322    
1323     c**********************************************************************
1324    
1325     c increment the counter for the accumulated fcc and qliq arrays
1326     c ---------------------------------------------------------------------
1327     imstturblw = imstturblw + 1
1328     imstturbsw = imstturbsw + 1
1329    
1330     c prevent ice or snow from melting
1331     c ---------------------------------------------------------------------
1332 molod 1.13 do i =1,nchptot
1333 molod 1.1 if( (icethk(i).gt.0.).and.(tground(i).gt.tice) ) tground(i)=tice
1334     enddo
1335    
1336     c Update tcanopy and ecanopy from the land points of the
1337     c tground and qground arrays
1338     c ---------------------------------------------------------------------
1339     do i =1,nchplnd
1340     tcanopy(i) = tground(i)
1341     ecanopy(i) = qground(i)
1342     enddo
1343    
1344     C Initialize Tendencies and Couplings
1345     c -----------------------------------
1346     do L = 1,nlay
1347     do i = 1,ijall
1348     duturb(i,1,L) = 0.
1349     dvturb(i,1,L) = 0.
1350     dtturb(i,1,L) = 0.
1351     qqgrid(i,1,L) = 0.
1352     qliqtmp(i,1,L) = 0.
1353     fcctmp(i,1,L) = 0.
1354     enddo
1355     do nt = 1,ntracers
1356     do i = 1,ijall
1357     dqturb(i,1,L,nt) = 0.
1358     enddo
1359     enddo
1360     enddo
1361    
1362     C Return Tendencies and Couplings to Grid Space
1363     c ---------------------------------------------
1364     do l = 1,nlay
1365 molod 1.13 call msc2grd(igrd,chfr,dumsc(1,L),nchp,nchptot,fracland,
1366 molod 1.1 . duturb(1,1,L),im,jm)
1367 molod 1.13 call msc2grd(igrd,chfr,dvmsc(1,L),nchp,nchptot,fracland,
1368 molod 1.1 . dvturb(1,1,L),im,jm)
1369 molod 1.13 call msc2grd(igrd,chfr,dtmsc(1,L),nchp,nchptot,fracland,
1370 molod 1.1 . dtturb(1,1,L),im,jm)
1371     do nt = 1,ntracers
1372 molod 1.13 call msc2grd(igrd,chfr,dqmsc(1,L,nt),nchp,nchptot,fracland,
1373 molod 1.1 . dqturb(1,1,L,nt),im,jm)
1374     enddo
1375 molod 1.13 call msc2grd(igrd,chfr, tke(1,L),nchp,nchptot,fracland,
1376 molod 1.1 . qqgrid(1,1,L),im,jm)
1377    
1378 molod 1.13 call msc2grd(igrd,chfr, fccmsc(1,L),nchp,nchptot,fracland,
1379 molod 1.1 . fcctmp(1,1,L),im,jm)
1380 molod 1.13 call msc2grd(igrd,chfr,qliqmsc(1,L),nchp,nchptot,fracland,
1381 molod 1.1 . qliqtmp(1,1,L),im,jm)
1382     enddo
1383    
1384     c Reduce clouds from conditionally unstable layer
1385     c -----------------------------------------------
1386     call ctei ( tz,qz,fcctmp,qliqtmp,plz,pkz,pkht,im*jm,nlay )
1387    
1388     C Bumb Total Cloud Liquid Water and Fraction by Instantanious Values
1389     c ------------------------------------------------------------------
1390     do l = 1,nlay
1391     do j=1,jm
1392     do i=1,im
1393     fccavesw(i,j,L) = fccavesw(i,j,L) + fcctmp(i,j,L)
1394     fccavelw(i,j,L) = fccavelw(i,j,L) + fcctmp(i,j,L)
1395     qliqavelw(i,j,L) = qliqavelw(i,j,L) + qliqtmp(i,j,L)
1396     qliqavesw(i,j,L) = qliqavesw(i,j,L) + qliqtmp(i,j,L)
1397     enddo
1398     enddo
1399    
1400     if (itrbfcc.gt.0) then
1401     do j=1,jm
1402     do i=1,im
1403 molod 1.6 qdiag(i,j,itrbfcc+L-1,bi,bj) = qdiag(i,j,itrbfcc+L-1,bi,bj) +
1404     . fcctmp(i,j,L)
1405 molod 1.1 enddo
1406     enddo
1407     endif
1408    
1409     if (itrbqliq.gt.0) then
1410     do j=1,jm
1411     do i=1,im
1412 molod 1.6 qdiag(i,j,itrbqliq+L-1,bi,bj)=qdiag(i,j,itrbqliq+L-1,bi,bj)+
1413 molod 1.1 . qliqtmp(i,j,L)*1.e6
1414     enddo
1415     enddo
1416     endif
1417     enddo
1418    
1419     C**********************************************************************
1420     C And some other variables to be transformed back to grid space:
1421     C Ground Temperature, snow depth and shallow layer ground wetness
1422     do j = 1,jm
1423     do i = 1,im
1424     tgz(i,j) = 0.
1425     enddo
1426     enddo
1427 molod 1.13 call msc2grd(igrd,chfr,tground ,nchp,nchptot,fracland,tgz ,im,jm)
1428 molod 1.1
1429     c *********************************************************************
1430     c **** increment diagnostic array for ground and surface temperatures,
1431     c *** ground temp tendency, and ground humidity
1432     c *********************************************************************
1433    
1434 molod 1.23 c if(itground.gt.0) then
1435     c do j =1,jm
1436     c do i =1,im
1437     c qdiag(i,j,itground,bi,bj) = qdiag(i,j,itground,bi,bj) + tgz(i,j)
1438     c enddo
1439     c enddo
1440     c endif
1441 molod 1.1
1442     if(itcanopy.gt.0) then
1443 molod 1.23 do j =1,jm
1444     do i =1,im
1445     qdiag(i,j,itcanopy,bi,bj) = qdiag(i,j,itcanopy,bi,bj) + tgz(i,j)
1446     enddo
1447 molod 1.1 enddo
1448     endif
1449    
1450     if(its.gt.0) then
1451 molod 1.23 do j =1,jm
1452     do i =1,im
1453     tmpdiag(i,j) = tz(i,j,nlay) * pkht(i,j,nlay)
1454     enddo
1455     enddo
1456     do j =1,jm
1457     do i =1,im
1458     qdiag(i,j,its,bi,bj) = qdiag(i,j,its,bi,bj) + tmpdiag(i,j)
1459 molod 1.1 enddo
1460     enddo
1461     endif
1462    
1463     if(idtg.gt.0) then
1464 molod 1.23 do j =1,jm
1465     do i =1,im
1466     qdiag(i,j,idtg,bi,bj) = qdiag(i,j,idtg,bi,bj) + tgz(i,j)
1467     enddo
1468 molod 1.1 enddo
1469     endif
1470    
1471     c *********************************************************************
1472     c **** increment diagnostic arrays for tendencies ****
1473     c *********************************************************************
1474     do L = 1,nlay
1475    
1476     if(iturbu.gt.0) then
1477 molod 1.23 do j =1,jm
1478     do i =1,im
1479     qdiag(i,j,iturbu+l-1,bi,bj) = qdiag(i,j,iturbu+l-1,bi,bj)
1480     . + duturb(i,j,l) * atimstp * secday
1481     enddo
1482 molod 1.1 enddo
1483     endif
1484    
1485     if(iturbv.gt.0) then
1486 molod 1.23 do j =1,jm
1487     do i =1,im
1488     qdiag(i,j,iturbv+l-1,bi,bj) = qdiag(i,j,iturbv+l-1,bi,bj)
1489     . + dvturb(i,j,l) * atimstp * secday
1490     enddo
1491 molod 1.1 enddo
1492     endif
1493    
1494     if(iturbq.gt.0) then
1495 molod 1.23 do j =1,jm
1496     do i =1,im
1497     qdiag(i,j,iturbq+l-1,bi,bj) = qdiag(i,j,iturbq+l-1,bi,bj)
1498     . + dqturb(i,j,l,1) * atimstp * secday * 1000
1499     enddo
1500 molod 1.1 enddo
1501     endif
1502    
1503     if(iturbt.gt.0) then
1504 molod 1.23 do j =1,jm
1505     do i =1,im
1506     qdiag(i,j,iturbt+l-1,bi,bj) = qdiag(i,j,iturbt+l-1,bi,bj)
1507 molod 1.25 . + dtturb(i,j,l) * pkz(i,j,l)*atimstp*secday
1508 molod 1.23 enddo
1509     enddo
1510     endif
1511    
1512 molod 1.1 enddo
1513    
1514     c pi-weight the theta and moisture tendencies
1515     c -------------------------------------------
1516     do i =1,ijall
1517     thtgz(i) = pz(i,1) * atimstp
1518     enddo
1519     do l =1,nlay
1520     do i =1,ijall
1521     duturb(i,1,l) = duturb(i,1,l)*atimstp
1522     dvturb(i,1,l) = dvturb(i,1,l)*atimstp
1523 molod 1.25 dtturb(i,1,l) = dtturb(i,1,l)*thtgz(i)
1524 molod 1.1 enddo
1525     do nt = 1,ntracers
1526     do i =1,ijall
1527 molod 1.25 dqturb(i,1,l,nt) = dqturb(i,1,l,nt)*thtgz(i)
1528 molod 1.1 enddo
1529     enddo
1530     enddo
1531    
1532     c *********************************************************************
1533     c **** zero out the accumulating rainfall and snowfall arrays ***
1534     c *********************************************************************
1535    
1536     if( time_left.lt.timstp ) then
1537     do j = 1,jm
1538     do i = 1,im
1539     rainlsp(i,j) = 0.
1540     rainconv(i,j) = 0.
1541     snowfall(i,j) = 0.
1542     enddo
1543     enddo
1544     endif
1545    
1546     c *********************************************************************
1547     c **** bump diagnostic counters ***
1548     c *********************************************************************
1549    
1550 molod 1.23 #ifdef ALLOW_DIAGNOSTICS
1551     if( (bi.eq.1) .and. (bj.eq.1) ) then
1552 molod 1.1 nturbt = nturbt + 1
1553     nturbq = nturbq + 1
1554     nturbu = nturbu + 1
1555     nturbv = nturbv + 1
1556     ntuflux = ntuflux + 1
1557     ntvflux = ntvflux + 1
1558     nttflux = nttflux + 1
1559     ntqflux = ntqflux + 1
1560     nwinds = nwinds + 1
1561     nkm = nkm + 1
1562     nkh = nkh + 1
1563     nri = nri + 1
1564     nct = nct + 1
1565     ncu = ncu + 1
1566     ntground = ntground + 1
1567     nts = nts + 1
1568     ndtg = ndtg + 1
1569     nqg = nqg + 1
1570     nqs = nqs + 1
1571     nhflux = nhflux + 1
1572     neflux = neflux + 1
1573     nevap = nevap + 1
1574     nuflux = nuflux + 1
1575     nvflux = nvflux + 1
1576     ndtsrf = ndtsrf + 1
1577     nustar = nustar + 1
1578     nz0 = nz0 + 1
1579     nfrqtrb = nfrqtrb + 1
1580     npbl = npbl + 1
1581     nu2m = nu2m + 1
1582     nv2m = nv2m + 1
1583     nt2m = nt2m + 1
1584     nq2m = nq2m + 1
1585     nu10m = nu10m + 1
1586     nv10m = nv10m + 1
1587     nt10m = nt10m + 1
1588     nq10m = nq10m + 1
1589     ntcanopy = ntcanopy + 1
1590     ntdeep = ntdeep + 1
1591     nqcanopy = nqcanopy + 1
1592     nsmshal = nsmshal + 1
1593     nsmroot = nsmroot + 1
1594     nsmdeep = nsmdeep + 1
1595     nsnow = nsnow + 1
1596     ncapacity = ncapacity + 1
1597     nraincon = nraincon + 1
1598     nrainlsp = nrainlsp + 1
1599     nsnowfall = nsnowfall + 1
1600     nrunoff = nrunoff + 1
1601     nfwsoil = nfwsoil + 1
1602     ngdrain = ngdrain + 1
1603     nsnowmelt = nsnowmelt + 1
1604     neresv = neresv + 1
1605     nesoil = nesoil + 1
1606     neveg = neveg + 1
1607     nesnow = nesnow + 1
1608     npardf = npardf + 1
1609     npardr = npardr + 1
1610     nlai = nlai + 1
1611     ngreen = ngreen + 1
1612     ndlwdtc = ndlwdtc + 1
1613     ndhdtc = ndhdtc + 1
1614     ndedtc = ndedtc + 1
1615     nevpot = nevpot + 1
1616     nlwgdown = nlwgdown + 1
1617     ndhdqa = ndhdqa + 1
1618     ndedqa = ndedqa + 1
1619     ndtc = ndtc + 1
1620     ndqc = ndqc + 1
1621     ntcdtc = ntcdtc + 1
1622     nraddtc = nraddtc + 1
1623     nsensdtc = nsensdtc + 1
1624     nlatdtc = nlatdtc + 1
1625     ntddtc = ntddtc + 1
1626     nqcdtc = nqcdtc + 1
1627     ntrbqliq = ntrbqliq + 1
1628     ntrbfcc = ntrbfcc + 1
1629    
1630 molod 1.23 nsdiag1 = nsdiag1 + 1
1631     nsdiag2 = nsdiag2 + 1
1632     endif
1633    
1634     #endif
1635    
1636 molod 1.1 return
1637     end
1638     SUBROUTINE TRBFLX (NN,TH,THV,SH,U,V,QQ,PL,PLE,PLK,PLKE,DPSTR,
1639     1 IWATER,Z0,tracers,ntrace,ntracedim,DTAU,ITRTRB,KMBG,KHBG,QBEG,
1640     2 TPROF,WU,WV,SRI,ET,EU,SWINDS,sustar,sz0,freqdg,pbldpth,
1641     3 sct,scu,stu2m,stv2m,stt2m,stq2m,stu10m,stv10m,stt10m,stq10m,
1642     4 irun,nlev,NYMD,NHMS,grav,cp,rgas,faceps,virtcon,undef,
1643     5 dshdthg,dshdshg,dthdthg,dthdshg,eturb,dedqa,dedtc,
1644     6 hsturb,dhsdqa,dhsdtc,transth,transsh,
1645     7 ctsave,xxsave,yysave,zetasave,xlsave,khsave,qliq,turbfcc)
1646     C**********************************************************************
1647     C SUBROUTINE TRBFLX - COMPUTES TURBULENT ADJUSTMENTS TO ATMOSPHERIC
1648     C PROFILE
1649     C - CALLED FROM PBL DRIVER
1650     C
1651     C ARGUMENTS ::
1652     C
1653     C INPUT:
1654     C ------
1655     C TH - POTENTIAL TEMPERATURE PROFILE
1656     C THV - VIRTUAL POTENTIAL TEMPERATURE PROFILE
1657     C SH - SPECIFIC HUMIDITY PROFILE
1658     C U - U - COMPONENT OF WIND PROFILE
1659     C V - V - COMPONENT OF WIND PROFILE
1660     C QQ - TURBULENT KINETIC ENERGY
1661     C PL - EVEN LEVEL PRESSURES
1662     C PLE - EDGE LEVEL PRESSURES
1663     C PLK - EVEN LEVEL PRESSURES ** KAPPA
1664     C PLKE - EDGE LEVEL PRESSURES ** KAPPA
1665     C DPSTR - PRESSURE INTERVALS
1666     C WATER - BIT ARRAY - '1' OVER OCEANS
1667     C Z0 - SURFACE ROUGHNESS
1668     C tracers - array of passive tracers
1669     C ntrace - number of tracers to be diffused
1670     C ntracedim - outer dimension of tracers array
1671 molod 1.20 C DTAU - TIME CHANGE PER ITERATION OF TRBFLX
1672     C ITRTRB - NUMBER OF ITERATIONS OF TRBFLX
1673 molod 1.1 C KMBG - BACKGROUND VALUE OF MOMENTUM TRANSFER COEF
1674     C KHBG - BACKGROUND VALUE OF HEAT TRANSFER COEF
1675     C NLEV - NUMBER OF ATMOSPHERIC LEVELS TO CALCULATE
1676     C QBEG - LOGICAL .TRUE. FOR INITIAL START OF GCM
1677     C TPROF - LOGICAL .TRUE. TO CALCULATE PT BY PT DIAGS
1678     C
1679     C OUTPUT:
1680     C -------
1681     C PROFILES RETURNED WITH UPDATED VALUES
1682     C
1683     C**********************************************************************
1684 molod 1.10 implicit none
1685    
1686     C Argument list declarations
1687 molod 1.11 integer nn,irun,nlev,ntrace,ntracedim,itrtrb,nhms,nymd
1688 molod 1.14 _RL TH(irun,NLEV+1),THV(irun,NLEV+1),SH(irun,NLEV+1)
1689     _RL U(irun,NLEV+1),V(irun,NLEV+1),QQ(irun,NLEV)
1690     _RL PL(irun,NLEV),PLE(irun,NLEV+1),PLK(irun,NLEV)
1691     _RL PLKE(irun,NLEV+1),DPSTR(irun,NLEV)
1692 molod 1.10 integer IWATER(irun)
1693 molod 1.14 _RL Z0(irun)
1694     _RL tracers(irun,nlev+1,ntracedim)
1695     _RL dtau,KMBG,KHBG
1696 molod 1.10 LOGICAL QBEG,TPROF
1697 molod 1.14 _RL SWINDS(irun)
1698     _RL SRI(irun,nlev), ET(irun,nlev)
1699     _RL EU (irun,nlev)
1700     _RL WU(irun,nlev)
1701     _RL WV (irun,nlev), pbldpth(irun)
1702     _RL sustar(irun), sz0(irun)
1703     _RL freqdg(irun,nlev-1)
1704     _RL sct(irun), scu(irun)
1705     _RL stu2m(irun),stv2m(irun),stt2m(irun),stq2m(irun)
1706     _RL stu10m(irun),stv10m(irun),stt10m(irun),stq10m(irun)
1707     _RL grav,cp,rgas,faceps,virtcon,undef
1708     _RL eturb(irun),dedqa(irun),dedtc(irun)
1709     _RL hsturb(irun),dhsdqa(irun),dhsdtc(irun)
1710     _RL dshdthg(irun,nlev),dthdthg(irun,nlev)
1711     _RL dshdshg(irun,nlev),dthdshg(irun,nlev)
1712     _RL transth(irun,nlev),transsh(irun,nlev)
1713     _RL ctsave(irun),xxsave(irun),yysave(irun)
1714     _RL zetasave(irun),xlsave(irun,nlev),khsave(irun,nlev)
1715     _RL qliq(irun,nlev),turbfcc(irun,nlev)
1716 molod 1.10
1717 molod 1.11 C Local Variables
1718 molod 1.14 _RL b1,b3,alpha,halpha,qqmin,qbustr
1719 molod 1.1 PARAMETER ( B1 = 16.6 )
1720     PARAMETER ( B3 = 1. / B1 )
1721     PARAMETER ( ALPHA = 0.1 )
1722     PARAMETER ( HALPHA = ALPHA * 0.5 )
1723     PARAMETER ( QQMIN = 0.005 )
1724     PARAMETER ( QBUSTR = 2.550952 )
1725 molod 1.14 _RL argmax, onethrd, z1pem25, b2, two
1726 molod 1.1 PARAMETER (ARGMAX = 30.)
1727     PARAMETER (ONETHRD = 1./3. )
1728     PARAMETER (Z1PEM25 = 1.E-25)
1729     PARAMETER ( B2 = 10.1 )
1730     PARAMETER ( two = 2.0 )
1731    
1732 molod 1.14 _RL AHS (irun), HS(irun)
1733     _RL XX (irun), YY(irun), CU(irun)
1734     _RL CT(irun), USTAR(irun)
1735     _RL RIB(irun), ZETA(irun), WS(irun)
1736     _RL DTHS(irun), DELTHS(irun)
1737     _RL DTHL(irun), DELTHL(irun)
1738     _RL RIBIN(irun),CUIN(irun)
1739     _RL CTIN(irun),ZETAIN(irun)
1740     _RL USTARIN(irun),RHOSIN(irun),Z0IN(irun)
1741     _RL qqcolmin(irun),qqcolmax(irun),levpbl(irun)
1742    
1743     _RL ADZ1(irun,nlev), DZ1TMP(irun,nlev)
1744     _RL DZ3(irun,nlev), TEMP(irun,nlev)
1745     _RL DV(irun,nlev), DTHV(irun,nlev)
1746     _RL DPK(irun,nlev), STRT(irun,nlev)
1747     _RL DW2(irun,nlev), RI(irun,nlev)
1748     _RL RHOZPK(irun,nlev), Q(irun,nlev)
1749     _RL RIINIT(irun,nlev), DU(irun,nlev)
1750     _RL QQINIT(irun,nlev), RHOKDZ(irun,nlev)
1751     _RL RHODZ2(irun,nlev)
1752     _RL KM(irun,nlev), KH(irun,nlev)
1753    
1754     _RL DELTH (irun,nlev+1), DELSH (irun,nlev+1)
1755     _RL FLXFAC (irun,nlev+1)
1756     _RL FLXFPK (irun,nlev+1)
1757    
1758     _RL ADZ2 (irun,nlev-1), RHODZ1(irun,nlev-1)
1759     _RL VKZE (irun,nlev-1), VKZM (irun,nlev-1)
1760     _RL XL (irun,nlev-1), QXLM (irun,nlev-1)
1761     _RL QQE (irun,nlev-1), QE (irun,nlev-1)
1762     _RL P3 (irun,nlev-1), XQ (irun,nlev-1)
1763     _RL XLDIAG (irun,nlev-1), FLXFCE(irun,nlev-1)
1764 molod 1.1
1765     LOGICAL FIRST,LAST
1766 molod 1.11 integer IBITSTB(irun,nlev),INTQ(irun,nlev)
1767 molod 1.1
1768     C arrays for use by moist bouyancy calculation
1769     C -----------------
1770 molod 1.14 _RL TL(irun,NLEV),DTH(irun,NLEV)
1771     _RL DSH(irun,NLEV)
1772     _RL SHL(irun,NLEV)
1773     _RL AA(irun,NLEV),BB(irun,NLEV),SSDEV(irun,NLEV)
1774     _RL ARG(irun,NLEV),XXZETA(irun),QBYU(irun)
1775     _RL SVAR(irun,NLEV),Q1M(irun,NLEV)
1776     _RL FCC(irun,NLEV)
1777     _RL BETAT(irun,NLEV),BETAW(irun,NLEV)
1778     _RL BETAL(irun,NLEV),BETAT1(irun,NLEV)
1779     _RL BETAW1(irun,NLEV),SBAR(irun,NLEV)
1780     _RL SHSAT(irun,NLEV)
1781 molod 1.1
1782     C Some space for variables to be used in called routines
1783     logical LWATER
1784     integer IVBITRIB(irun)
1785 molod 1.14 _RL VHZ(irun)
1786     _RL VH0(irun)
1787     _RL VPSIM(irun),VAPSIM(irun)
1788     _RL VPSIG(irun),VPSIHG(irun)
1789     _RL VTEMP(irun),VDZETA(irun)
1790     _RL VDZ0(irun),VDPSIM(irun)
1791     _RL VDPSIH(irun),VZH(irun)
1792     _RL VXX0(irun),VYY0(irun)
1793     _RL VAPSIHG(irun),VRIB1(irun),VWS1(irun)
1794     _RL VPSIH(irun),VZETAL(irun)
1795     _RL VZ0L(irun),VPSIH2(irun)
1796     _RL VX0PSIM(irun),VG(irun),VG0(irun),VR1MG0(irun)
1797     _RL VZ2(irun),VDZSEA(irun),VAZ0(irun),VXNUM1(irun)
1798     _RL VPSIGB2(irun),VDX(irun),VDXPSIM(irun),VDY(irun)
1799     _RL VXNUM2(irun),VDEN(irun),VAWS1(irun),VXNUM3(irun)
1800     _RL VXNUM(irun),VDZETA1(irun),VDZETA2(irun)
1801     _RL VZCOEF2(irun),VZCOEF1(irun),VTEMPLIN(irun)
1802     _RL VDPSIMC(irun),VDPSIHC(irun)
1803    
1804     _RL DZITRP(irun,nlev-1),STBFCN(irun,nlev)
1805     _RL XL0(irun,nlev),Q1(irun,nlev-1)
1806     _RL WRKIT1(irun,nlev-1)
1807     _RL WRKIT2(irun,nlev-1)
1808     _RL WRKIT3(irun,nlev-1)
1809     _RL WRKIT4(irun,nlev-1)
1810 molod 1.1 INTEGER INT1(irun,nlev), INT2(irun,nlev-1)
1811    
1812 molod 1.14 _RL vrt1con,pi,rsq2pi,p5sr,clh,vk,rvk,aitr,gbycp,fac1,fac2
1813     _RL getcon,dum,errf
1814 molod 1.11 integer istnlv,nlevm1,nlevm2,nlevml,nlevp1,istnm1,istnm2,istnp1
1815     integer istnml,istnmq,istlmq,nlevmq
1816     integer i,iter,init,n,nt,LL,L,Lp,Lp1,lmin,lminq,lminq1,ibit
1817 molod 1.1
1818     vk = getcon('VON KARMAN')
1819     rvk = 1./vk
1820     AITR = 1. / FLOAT(ITRTRB)
1821     ISTNLV = irun * NLEV
1822     NLEVM1 = NLEV - 1
1823     NLEVM2 = NLEV - 2
1824     NLEVP1 = NLEV + 1
1825     ISTNM1 = irun * NLEVM1
1826     ISTNM2 = irun * NLEVM2
1827     ISTNP1 = irun * NLEVP1
1828     GBYCP = GRAV / CP
1829    
1830     VRT1CON = 1. + VIRTCON
1831     PI = 4. * ATAN(1.)
1832     RSQ2PI = 1./ ((2.*PI)**0.5)
1833     P5SR = 0.5**0.5
1834     CLH = GETCON('LATENT HEAT COND') / CP
1835    
1836     C SET INITIAL NUMBER OF ITERATIONS OF SFCFLX
1837     C ------------------------------------------
1838     N = 6
1839     C DETERMINE IF INITIAL START
1840     C --------------------------
1841     INIT = 0
1842     IF(QBEG) INIT = 1
1843     C SET DIAGNOSTIC LOGICALS AND INITIALIZE DIAGNOSTIC ARRAYS
1844     C --------------------------------------------------------
1845     do I =1,istnlv
1846     wu(i,1) = 0.
1847     enddo
1848     do I =1,istnlv
1849     wv(i,1) = 0.
1850     enddo
1851     do I =1,istnlv
1852     eu(i,1) = 0.
1853     enddo
1854     do I =1,istnlv
1855     et(i,1) = 0.
1856     enddo
1857     if (tprof) then
1858     DO I =1,ISTNM1
1859     XLDIAG(I,1) = 0.
1860     enddo
1861     endif
1862     do I =1,irun
1863     wu(i,nlev) = 0.
1864     enddo
1865     do I =1,irun
1866     wv(i,nlev) = 0.
1867     enddo
1868     do I =1,irun
1869     scu(i) = 0.
1870     enddo
1871     do I =1,irun
1872     sct(i) = 0.
1873     enddo
1874     do I =1,irun
1875     pbldpth(i) = 0.
1876     enddo
1877     do I =1,irun
1878     sustar(i) = 0.
1879     enddo
1880     do I =1,irun
1881     sz0(i) = 0.
1882     enddo
1883     do I =1,ISTNM1
1884     FREQDG(I,1) = 0.
1885     enddo
1886     do I =1,irun
1887     stu2m(i) = 0.
1888     enddo
1889     do I =1,irun
1890     stv2m(i) = 0.
1891     enddo
1892     do I =1,irun
1893     stt2m(i) = 0.
1894     enddo
1895     do I =1,irun
1896     stq2m(i) = 0.
1897     enddo
1898     do I =1,irun
1899     stu10m(i) = 0.
1900     enddo
1901     do I =1,irun
1902     stv10m(i) = 0.
1903     enddo
1904     do I =1,irun
1905     stt10m(i) = 0.
1906     enddo
1907     do I =1,irun
1908     stq10m(i) = 0.
1909     enddo
1910    
1911     IF (INIT.EQ.1) THEN
1912     DO I = 1,ISTNM1
1913     XLSAVE(I,1) = 0.
1914     KHSAVE(I,1) = 0.
1915     ENDDO
1916     DO I = 1,irun
1917     CTSAVE(I) = 0.
1918     XXSAVE(I) = 0.
1919     YYSAVE(I) = 0.
1920     ZETASAVE(I) = 0.
1921     ENDDO
1922     ENDIF
1923    
1924     C COMPUTE VERTICAL GRID
1925     C ---------------------
1926     DO 9038 I =1,ISTNLV
1927     ADZ1(I,1) = (CP/GRAV)*(PLKE(I,2)-PLKE(I,1))
1928     ADZ1(I,1) = THV(I,1) * ADZ1(I,1)
1929     DZ1TMP(I,1) = ADZ1(I,1)
1930     9038 CONTINUE
1931     DO 9040 I =1,ISTNM1
1932     ADZ2(I,1) = 0.5 * (ADZ1(I,1)+ADZ1(I,2))
1933     9040 CONTINUE
1934     C DEPTH HS OF SURFACE LAYER
1935     C -------------------------
1936     DO 9042 I =1,irun
1937     HS(I) = 0.5 * ADZ1(I,NLEV)
1938     9042 CONTINUE
1939     C ALPHA * LAYER DEPTHS FOR TRBLEN
1940     C -------------------------------
1941     DO 9044 I =1,irun
1942     DZ3(I,1) = HALPHA * ADZ1(I,1)
1943     9044 CONTINUE
1944     DO 9046 I =1,ISTNM2
1945     DZ3(I,2) = ALPHA * ADZ1(I,2)
1946     9046 CONTINUE
1947     DO 9048 I =1,irun
1948     DZ3(I,NLEV) = ALPHA * HS(I)
1949     9048 CONTINUE
1950    
1951     C VK * HEIGHTS AT MID AND EDGE LEVELS
1952     C -----------------------------------
1953     DO 9050 I =1,ISTNM1
1954     TEMP(I,2) = VK * ADZ1(I,2)
1955     9050 CONTINUE
1956     DO 9052 I =1,irun
1957     VKZE(I,NLEVM1) = TEMP(I,NLEV)
1958     9052 CONTINUE
1959     DO 100 LL = 2,NLEVM1
1960     L = NLEV - LL
1961     LP1 = L + 1
1962     DO 9054 I =1,irun
1963     VKZE(I,L) = VKZE(I,LP1) + TEMP(I,LP1)
1964     9054 CONTINUE
1965     100 CONTINUE
1966     DO 9056 I =1,ISTNM1
1967     VKZM(I,1) = VKZE(I,1) - 0.5 * TEMP(I,2)
1968     9056 CONTINUE
1969     C COMPUTE RHO BY DZ AT MID AND EDGE LEVELS
1970     C ----------------------------------------
1971     DO 200 L = 1,NLEVM1
1972     LP1 = L + 1
1973     DO 9058 I =1,irun
1974     FAC1 = DPSTR(I,L) / ( DPSTR(I,L) + DPSTR(I,LP1) )
1975     FAC2 = 1. - FAC1
1976     RHODZ2(I,L) = FAC1 * THV(I,LP1)
1977     RHODZ2(I,L) = RHODZ2(I,L) + FAC2 * THV(I,L)
1978     9058 CONTINUE
1979     200 CONTINUE
1980     DO 9060 I =1,ISTNM1
1981     RHODZ2(I,1) = (RGAS*0.01) * RHODZ2(I,1)
1982     TEMP(I,1) = PLKE(I,2) * ADZ2(I,1)
1983     RHODZ2(I,1) = TEMP(I,1) * RHODZ2(I,1)
1984     RHODZ2(I,1) = PLE(I,2) / RHODZ2(I,1)
1985     RHOZPK(I,1) = RHODZ2(I,1) * PLKE(I,2)
1986     RHODZ1(I,1) = (RGAS*0.01) * THV(I,2)
1987     TEMP(I,1) = PLK(I,2) * ADZ1(I,2)
1988     RHODZ1(I,1) = TEMP(I,1) * RHODZ1(I,1)
1989     RHODZ1(I,1) = PL(I,2) / RHODZ1(I,1)
1990     9060 CONTINUE
1991     C COMPUTE FLXFAC FOR LAYERS AND EDGES
1992     C COMPUTE DTG / DT DUE TO RADIATION AND HEAT CONDUCTION THROUGH ICE
1993     C -----------------------------------------------------------------
1994     DO 9062 I =1,ISTNLV
1995     FLXFPK(I,1) = PLE(I,2) - PLE(I,1)
1996     FLXFPK(I,1) = FLXFPK(I,1) * PLK(I,1)
1997     FLXFPK(I,1) = (GRAV*DTAU*0.01) / FLXFPK(I,1)
1998     9062 CONTINUE
1999     DO 9064 I =1,irun
2000     FLXFPK(I,NLEVP1) = 0.
2001     9064 CONTINUE
2002     DO 9066 I =1,irun
2003     IF (IWATER(I).EQ.0 ) FLXFPK(I,NLEVP1) = 1. / PLKE(I,NLEVP1)
2004     9066 CONTINUE
2005     DO 9068 I =1,ISTNLV
2006     FLXFAC(I,1) = FLXFPK(I,1) * PLK(I,1)
2007     9068 CONTINUE
2008     DO 9070 I =1,irun
2009     FLXFAC(I,NLEVP1) = FLXFPK(I,NLEVP1)
2010     9070 CONTINUE
2011     DO 9074 I =1,irun
2012     FLXFPK(I,NLEVP1) = CP * FLXFPK(I,NLEVP1)
2013     9074 CONTINUE
2014     DO 9076 I =1,ISTNM1
2015     FLXFCE(I,1) = PL(I,2) - PL(I,1)
2016     9076 CONTINUE
2017     DO 9078 I =1,ISTNM1
2018     FLXFCE(I,1) = (GRAV*DTAU*0.01) / FLXFCE(I,1)
2019     9078 CONTINUE
2020     C COMPUTE RECIPROCALS OF DZ1, DZ2, HS
2021     C -----------------------------------
2022     DO 9084 I =1,ISTNLV
2023     ADZ1(I,1) = 1. / ADZ1(I,1)
2024     9084 CONTINUE
2025     DO 9086 I =1,ISTNM1
2026     ADZ2(I,1) = 1. / ADZ2(I,1)
2027     9086 CONTINUE
2028     DO 9088 I =1,irun
2029     AHS(I) = 1. / HS(I)
2030     9088 CONTINUE
2031     C COMPUTE GRADIENTS OF P**KAPPA
2032     C -----------------------------
2033     DO 9090 I =1,ISTNM1
2034     DPK(I,1) = ( PLK(I,2)-PLK(I,1) ) * ADZ2(I,1)
2035     9090 CONTINUE
2036     DO 9092 I =1,irun
2037     DPK(I,NLEV) = GBYCP / THV(I,NLEV)
2038     9092 CONTINUE
2039     C INITIALIZE Q ARRAY
2040     C ------------------
2041     DO 9094 I =1,ISTNM1
2042     Q(I,1) = 2. * QQ(I,1)
2043     Q(I,1) = SQRT( Q(I,1) )
2044     9094 CONTINUE
2045     FIRST = .TRUE.
2046     LAST = .FALSE.
2047     C**********************************************************************
2048     C**********************************************************************
2049     C MAIN LOOP
2050     C
2051     DO 2000 ITER = 1, ITRTRB
2052     C
2053     IF ( ITER .GE. ITRTRB ) LAST = .TRUE.
2054     C
2055     C CODE FOR MOIST BOUNDARY LAYER - NEW CALCULATION OF DTHV
2056     C
2057     IF(ITER.EQ.1) THEN
2058     DO I = 1,irun
2059     CT(I) = CTSAVE(I)
2060     XX(I) = XXSAVE(I)
2061     YY(I) = YYSAVE(I)
2062     ZETA(I) = ZETASAVE(I)
2063     ENDDO
2064     ENDIF
2065     C
2066     DO I = 1,irun
2067     TL(I,NLEV) = TH(I,NLEV)*PLK(I,NLEV)
2068     call qsat ( tl(i,nlev),pl(i,nlev),shsat(i,nlev),dum,.false. )
2069     ENDDO
2070    
2071     DO I = 1,irun
2072     BB(I,NLEV) = FACEPS*SHSAT(I,NLEV)/(TL(I,NLEV)*TL(I,NLEV))
2073     AA(I,NLEV) = 1. / (1. + CLH * BB(I,NLEV) )
2074     BB(I,NLEV) = BB(I,NLEV) * AA(I,NLEV) * plk(I,nlev)
2075     DTH(I,NLEV) = TH(I,NLEV)-TH(I,NLEVP1)
2076     DSH(I,NLEV) = SH(I,NLEV)-SH(I,NLEVP1)
2077     SBAR(I,NLEV) = AA(I,NLEV) * (SH(I,NLEV) - SHSAT(I,NLEV))
2078     SSDEV(I,NLEV)=CT(I)*(AA(I,NLEV)*DSH(I,NLEV)
2079     1 -BB(I,NLEV)*DTH(I,NLEV))
2080     XXZETA(I) = XX(I)-ZETA(I)
2081     IF(XXZETA(I).LT.0.1*XX(I)) XXZETA(I)=0.1*XX(I)
2082     IF(XXZETA(I).LE.0.) XXZETA(I)=0.1
2083     QBYU(I) =QBUSTR * XXZETA(I) ** ONETHRD
2084     SSDEV(I,NLEV) = B2*YY(I)*SSDEV(I,NLEV)*SSDEV(I,NLEV)/QBYU(I)
2085     SVAR(I,NLEV) = SQRT(SSDEV(I,NLEV))
2086     IF ( SVAR(I,NLEV).LT.Z1PEM25) SVAR(I,NLEV) = Z1PEM25
2087     Q1M(I,NLEV) = SBAR(I,NLEV) / SVAR(I,NLEV)
2088     FCC(I,NLEV) = (1./2.) * ( 1. + ERRF( P5SR*Q1M(I,NLEV) ) )
2089     SHL(I,NLEV) = FCC(I,NLEV) * SBAR(I,NLEV)
2090     ARG(I,NLEV) = (1./2.)*Q1M(I,NLEV)*Q1M(I,NLEV)
2091     IF(ARG(I,NLEV).LE.ARGMAX)
2092     1 SHL(I,NLEV) = SHL(I,NLEV)+RSQ2PI*SVAR(I,NLEV)*EXP(-ARG(I,NLEV))
2093     BETAT(I,NLEV) = 1. + VIRTCON*SH(I,NLEV) - VRT1CON*SHL(I,NLEV)
2094     BETAW(I,NLEV) = VIRTCON *
2095     1 ( TH(I,NLEV) + CLH * SHL(I,NLEV) * (1./plk(i,nlev)) )
2096     BETAL(I,NLEV) = (1.+VIRTCON*SH(I,NLEV)-TWO*VRT1CON*SHL(I,NLEV))
2097     1 * (1./plk(i,nlev)) * CLH - VRT1CON * TH(I,NLEV)
2098     BETAT1(I,NLEV) = BETAT(I,NLEV) - BB(I,NLEV)*FCC(I,NLEV)
2099     1 * BETAL(I,NLEV)
2100     BETAW1(I,NLEV) = BETAW(I,NLEV) + AA(I,NLEV) * FCC(I,NLEV)
2101     1 * BETAL(I,NLEV)
2102     DTHV(I,NLEV) = BETAT1(I,NLEV)*DTH(I,NLEV) +
2103     1 BETAW1(I,NLEV)*DSH(I,NLEV)
2104     THV(I,NLEVP1) = THV(I,NLEV) - DTHV(I,NLEV)
2105     ENDDO
2106    
2107     C SURFACE FLUX TRANSFER COEFFICIENTS
2108     C
2109     CALL SFCFLX(NN,U(1,NLEV),V(1,NLEV),
2110     1 THV(1,NLEV),
2111     2 THV(1,NLEVP1),TH(1,NLEV),TH(1,NLEVP1),
2112     3 SH(1,NLEV),SH(1,NLEVP1),PLK(1,NLEV),
2113     4 PLKE(1,NLEVP1),PLE(1,NLEVP1),Z0,
2114     5 IWATER,HS,AHS,
2115     6 FIRST,LAST,N,irun,aitr,RHODZ2(1,NLEV),RHOZPK(1,NLEV),
2116     7 KH(1,NLEV),KM(1,NLEV),USTAR,
2117     8 XX,YY,CU,
2118     9 CT,RIB,ZETA,WS,
2119     1 stu2m,stv2m,stt2m,stq2m,stu10m,stv10m,stt10m,stq10m,
2120     2 cp,rgas,undef,
2121     3 lwater, ivbitrib,
2122     4 VHZ,VPSIM,VAPSIM,VPSIG,VPSIHG,VTEMP,VDZETA,VDZ0,VDPSIM,
2123     5 VDPSIH,VZH,VXX0,VYY0,VAPSIHG,VRIB1,VWS1,VPSIH,
2124     9 VZETAL,VZ0L,VPSIH2,VH0,
2125     1 VX0PSIM,VG,VG0,VR1MG0,VZ2,VDZSEA,VAZ0,VXNUM1,VPSIGB2,VDX,
2126     2 VDXPSIM,VDY,VXNUM2,VDEN,VAWS1,VXNUM3,VXNUM,VDZETA1,VDZETA2,
2127     3 VZCOEF2,VZCOEF1,VTEMPLIN,VDPSIMC,VDPSIHC)
2128     CI
2129     C
2130     N = 1
2131     C
2132     C SET VALUES OF TURBULENT VELOCITY AND KINETIC ENERGY AT THE GROUND
2133     C
2134     CB
2135     DO 9098 I =1,irun
2136     Q(I,NLEV) = QBUSTR * USTAR(I)
2137     QQ(I,NLEV) = 0.5 * Q(I,NLEV) * Q(I,NLEV)
2138     9098 CONTINUE
2139     CE
2140     C
2141     C GRADIENTS
2142     C ---------
2143     DO 9100 I =1,ISTNM1
2144     DU(I,1) = ( U(I,1)- U(I,2) ) * ADZ2(I,1)
2145     DV(I,1) = ( V(I,1)- V(I,2) ) * ADZ2(I,1)
2146     9100 CONTINUE
2147    
2148    
2149     C NEW CODE FOR MOIST BOUNDARY LAYER - NEW CALCULATION OF DTHV
2150     C
2151     IF(ITER.EQ.1) THEN
2152     DO I = 1,ISTNM1
2153     XL(I,1) = XLSAVE(I,1)
2154     ENDDO
2155     ENDIF
2156     C
2157     DO I =1,ISTNM1
2158     DTH(I,1) = ( TH(I,1)-TH(I,2) ) * ADZ2(I,1)
2159     DSH(I,1) = ( SH(I,1)-SH(I,2) ) * ADZ2(I,1)
2160     TL(I,1) = TH(I,1)*PLK(I,1)
2161     ENDDO
2162     DO LL = 1,NLEVM1
2163     DO I = 1,irun
2164     call qsat ( tl(i,LL),pl(i,LL),shsat(i,LL),dum,.false. )
2165     ENDDO
2166     ENDDO
2167     DO I = 1,ISTNM1
2168     BB(I,1) = FACEPS*SHSAT(I,1)/(TL(I,1)*TL(I,1))
2169     AA(I,1) = 1. / (1. + CLH * BB(I,1) )
2170     COMMM BB(I,1) = BB(I,1) * AA(I,1) * plke(I,2)
2171     BB(I,1) = BB(I,1) * AA(I,1)
2172     SBAR(I,1) = AA(I,1) * (SH(I,1) - SHSAT(I,1))
2173     ENDDO
2174     DO I = 1,irun
2175     COMMM SSDEV(I,1) = XL(I,1)*(AA(I,1)*DSH(I,1)-BB(I,1)*DTH(I,1))
2176     SSDEV(I,1) = XL(I,1)*(AA(I,1)*DSH(I,1)-
2177     1 BB(I,1)*plke(I,2)*DTH(I,1))
2178     SSDEV(I,1) = B2 * KHSAVE(I,1) * SSDEV(I,1) * SSDEV(I,1)
2179     SVAR(I,1) = SQRT(SSDEV(I,1))
2180     IF ( SVAR(I,1).LT.Z1PEM25) SVAR(I,1) = Z1PEM25
2181     ENDDO
2182     DO I = 1,ISTNM2
2183     COMMM SSDEV(I,2) = XL(I,1)*(AA(I,2)*DSH(I,1)-BB(I,2)*DTH(I,1))
2184     SSDEV(I,2) = XL(I,1)*(AA(I,2)*DSH(I,1)-
2185     1 BB(I,2)*plke(I,2)*DTH(I,1))
2186     SSDEV(I,2) = B2 * KHSAVE(I,1) * SSDEV(I,1) * SSDEV(I,1)
2187     SVAR(I,2) = SQRT(SSDEV(I,2))
2188     COMMM SSDEV(I,2) = XL(I,2)*(AA(I,2)*DSH(I,2)-BB(I,2)*DTH(I,2))
2189     SSDEV(I,2) = XL(I,2)*(AA(I,2)*DSH(I,2)-
2190     1 BB(I,2)*plke(I,3)*DTH(I,2))
2191     SSDEV(I,2) = B2 * KHSAVE(I,2) * SSDEV(I,2) * SSDEV(I,2)
2192     TEMP(I,2) = SQRT(SSDEV(I,2))
2193     SVAR(I,2) = (1./2.) * (SVAR(I,2) + TEMP(I,2))
2194     IF ( SVAR(I,2).LT.Z1PEM25) SVAR(I,2) = Z1PEM25
2195     ENDDO
2196     DO I = 1,ISTNM1
2197     Q1M(I,1) = SBAR(I,1) / SVAR(I,1)
2198     FCC(I,1) = (1./2.) * ( 1. + ERRF( P5SR*Q1M(I,1) ) )
2199     SHL(I,1) = FCC(I,1) * SBAR(I,1)
2200     ARG(I,1) = (1./2.)*Q1M(I,1)*Q1M(I,1)
2201     IF(ARG(I,1).LE.ARGMAX)
2202     1 SHL(I,1) = SHL(I,1)+RSQ2PI*SVAR(I,1)*EXP(-ARG(I,1))
2203     BETAT(I,1) = 1. + VIRTCON * SH(I,1) - VRT1CON * SHL(I,1)
2204     BETAW(I,1) = VIRTCON *
2205     1 ( TH(I,1) + (CLH/plk(I,1)) * SHL(I,1) )
2206     BETAL(I,1) = ( 1. + VIRTCON*SH(I,1) - TWO*VRT1CON*SHL(I,1) )
2207     1 * (CLH/plke(I,2)) - VRT1CON * TH(I,1)
2208     COMMM BETAT1(I,1) = BETAT(I,1) - BB(I,1) * FCC(I,1) * BETAL(I,1)
2209     BETAT1(I,1) = BETAT(I,1) -
2210     1 BB(I,1)*plk(i,1) * FCC(I,1) * BETAL(I,1)
2211     BETAW1(I,1) = BETAW(I,1) + AA(I,1) * FCC(I,1) * BETAL(I,1)
2212     ENDDO
2213     DO I = 1,ISTNM1
2214     DTHV(I,1) = (1./2.)*((BETAT1(I,1)+BETAT1(I,2))*DTH(I,1)
2215     1 + (BETAW1(I,1)+BETAW1(I,2))*DSH(I,1))
2216     ENDDO
2217    
2218     C GRADIENTS AT THE TOP OF THE SURFACE LAYER
2219     C -----------------------------------------
2220     DO 9102 I =1,irun
2221     DU(I,NLEV) = CU(I)*XX(I)*AHS(I)*RVK
2222     DV(I,NLEV) = V(I,NLEV) * DU(I,NLEV)
2223     DU(I,NLEV) = U(I,NLEV) * DU(I,NLEV)
2224     DTHV(I,NLEV) = CT(I) * YY(I) *
2225     1 ((THV(I,NLEV)-THV(I,NLEVP1)) * RVK)* AHS(I)
2226     9102 CONTINUE
2227    
2228     C CALCULATE BRUNT-VAISALA FREQUENCIES, SHEARS, RICHARDSON NUMBERS
2229     C ---------------------------------------------------------------
2230     DO 9104 I =1,ISTNLV
2231     STRT(I,1) = CP * DTHV(I,1) * DPK(I,1)
2232     DW2(I,1) = DU(I,1) * DU(I,1) + DV(I,1) * DV(I,1)
2233     IF ( DW2(I,1) .LE. 1.e-4 ) DW2(I,1) = 1.e-4
2234     RI(I,1) = STRT(I,1) / DW2(I,1)
2235     9104 CONTINUE
2236     C FILL RICHARDSON NUMBER AND SURFACE WIND DIAGNOSTICS
2237     C (THOSE NEEDED FROM FIRST TRBFLX ITERATION)
2238     C ---------------------------------------------------
2239     DO 9106 I =1,ISTNM1
2240     SRI(I,1) = RI(I,1)
2241     9106 CONTINUE
2242     DO 9108 I =1,irun
2243     SRI(I,NLEV) = RIB(I)
2244     9108 CONTINUE
2245     DO 9110 I =1,irun
2246     SWINDS(I) = WS(I)
2247     9110 CONTINUE
2248     C INITIALIZE KH, KM, QE AND P3 AND ELIMINATE SMALL QQ
2249     C ---------------------------------------------------
2250     DO 9112 I =1,ISTNM1
2251     KH(I,1) = 0.
2252     KM(I,1) = 0.
2253     QQE(I,1) = 0.
2254     QE(I,1) = 0.
2255     P3(I,1) = 0.
2256     9112 CONTINUE
2257     DO 9414 I = 1,ISTNM1
2258     IBITSTB(I,1) = 0
2259     9414 CONTINUE
2260     DO 9314 I = 1,ISTNM1
2261     IF ( QQ(I,1) .GT. 1.e-8 ) THEN
2262     INTQ(I,1) = 1
2263     ELSE
2264     INTQ(I,1) = 0
2265     ENDIF
2266     9314 CONTINUE
2267     DO 9114 I = 1,ISTNM1
2268     IF ( QQ(I,1).LE.1.e-8 ) THEN
2269     QQ(I,1) = 0.
2270     Q(I,1) = 0.
2271     ENDIF
2272     9114 CONTINUE
2273     C
2274     DO 300 LMINQ = 1,NLEVM1
2275     IBIT = 0
2276     DO 9116 I = 1,irun
2277     IF ( QQ(I,LMINQ).GT.1.e-8 ) IBIT = IBIT + 1
2278     9116 CONTINUE
2279     IF(IBIT.GE.1)GO TO 310
2280     300 CONTINUE
2281     LMINQ = NLEV-1
2282     310 CONTINUE
2283     LMINQ = 1
2284     LMINQ1 = 1
2285     IF(LMINQ.GT.1)LMINQ1 = LMINQ - 1
2286     C LENGTH SCALE
2287     C ------------
2288     CALL TRBLEN(STRT,DW2,DZ3,Q,VKZE,VKZM,DTHV,DPK,DU,DV,XL,QXLM,
2289     1 NLEV,INIT,LMIN,LMINQ,LMINQ1,CP,INT1,INT2,
2290     2 DZITRP,STBFCN,XL0,Q1,WRKIT1,WRKIT2,WRKIT3,WRKIT4,irun)
2291     C QE AND DIMENSIONLESS COEFFS FROM LEVEL 2 MODEL
2292     C ----------------------------------------------
2293     IF( LMIN .LT. NLEV ) THEN
2294     NLEVML = NLEV - LMIN
2295     CALL TRBL20(RI(1,LMIN),STRT(1,LMIN),DW2(1,LMIN),XL(1,LMIN),
2296     1 KM(1,LMIN),KH(1,LMIN),QE(1,LMIN),QQE(1,LMIN),IBITSTB(1,LMIN),
2297     2 NLEVML,nlev,irun)
2298     ENDIF
2299     C FOR INITIAL START ONLY : USE EQUILIBRIUM MODEL
2300     C ----------------------------------------------
2301     IF ( INIT .EQ. 1 ) THEN
2302     DO 9180 I =1,ISTNM1
2303     QQ(I,1) = QQE(I,1)
2304     Q(I,1) = QE(I,1)
2305     9180 CONTINUE
2306     INIT = 2
2307     CALL TRBLEN(STRT,DW2,DZ3,Q,VKZE,VKZM,DTHV,DPK,DU,DV,XL,QXLM,
2308     1 NLEV,INIT,LMIN,LMINQ,LMINQ1,CP,INT1,INT2,
2309     2 DZITRP,STBFCN,XL0,Q1,WRKIT1,WRKIT2,WRKIT3,WRKIT4,irun)
2310     INIT = 0
2311     GO TO 550
2312     ENDIF
2313     C DIMENSIONLESS COEFFS AND P3 (Q LE QE)
2314     C -------------------------------------
2315     IF( LMIN .LT. NLEV ) THEN
2316     ISTNML = irun * NLEVML
2317     DO 9320 I = 1,ISTNML
2318     IF ( (IBITSTB(I,LMIN).EQ.1) .AND.
2319     1 ( Q(I,LMIN) .LE. QE(I,LMIN) ) ) THEN
2320     IBITSTB(I,LMIN) = 1
2321     ELSE
2322     IBITSTB(I,LMIN) = 0
2323     ENDIF
2324     9320 CONTINUE
2325     DO 9220 I = 1,ISTNML
2326     IF(IBITSTB(I,LMIN).EQ.1 ) THEN
2327     TEMP(I,LMIN) = Q(I,LMIN) / QE(I,LMIN)
2328     KH(I,LMIN) = TEMP(I,LMIN) * KH(I,LMIN)
2329     KM(I,LMIN) = TEMP(I,LMIN) * KM(I,LMIN)
2330     ENDIF
2331     TEMP(I,LMIN) = 0.01 * QQE(I,LMIN)
2332     IF((IBITSTB(I,LMIN).EQ.1) .AND.
2333     1 ( QQ(I,LMIN) .LE. TEMP(I,LMIN) )) THEN
2334     QQ(I,LMIN) = TEMP(I,LMIN)
2335     Q(I,LMIN) = 0.1 * QE(I,LMIN)
2336     ENDIF
2337     IF(IBITSTB(I,LMIN).EQ.1 ) P3(I,LMIN) = (2.*B3) *
2338     1 ( QE(I,LMIN) - Q(I,LMIN) )
2339     9220 CONTINUE
2340     ENDIF
2341     C DIMENSIONLESS COEFFS AND P3 (Q GT QE)
2342     C -------------------------------------
2343     NLEVML = NLEV - LMINQ
2344     CALL TRBL25(Q(1,LMINQ),XL(1,LMINQ),STRT(1,LMINQ),DW2(1,LMINQ),
2345     1 IBITSTB(1,LMINQ),INTQ(1,LMINQ),KM(1,LMINQ),KH(1,LMINQ),
2346     2 P3(1,LMINQ),NLEVML,nlev,irun)
2347     C CALCULATE SOURCE TERM P3
2348     C ------------------------
2349     IF ( LMINQ .LT. LMIN ) THEN
2350     LMIN = LMINQ
2351     ISTNML = irun * ( NLEV - LMIN )
2352     ENDIF
2353     IF( LMIN .LT. NLEV ) THEN
2354     DO 9122 I =1,ISTNML
2355     P3(I,LMIN) = P3(I,LMIN) * DTAU / XL(I,LMIN)
2356     TEMP(I,LMIN) = QQE(I,LMIN) * P3(I,LMIN)
2357     XQ(I,LMIN) = QQE(I,LMIN) - QQ(I,LMIN)
2358     9122 CONTINUE
2359     DO 9216 I = 1,ISTNML
2360     IF( ( (IBITSTB(I,LMIN).EQ.1) .AND.
2361     1 ( XQ(I,LMIN) .LT. TEMP(I,LMIN) ) )
2362     2 .OR.
2363     3 ( (IBITSTB(I,LMIN).EQ.0) .AND.
2364     4 ( XQ(I,LMIN) .GT. TEMP(I,LMIN) ) ) )
2365     5 P3(I,LMIN) = XQ(I,LMIN) / QQE(I,LMIN)
2366     9216 CONTINUE
2367     ENDIF
2368     550 CONTINUE
2369     C DIAGNOSTIC PROFILES : INITIAL RI AND QQ
2370     C ---------------------------------------
2371     IF ( TPROF .AND. FIRST ) THEN
2372     DO 9118 I =1,irun
2373     RIBIN(I) = RIB(I)
2374     CUIN(I) = CU(I)
2375     CTIN(I) = CT(I)
2376     USTARIN(I) = USTAR(I)
2377     RHOSIN(I) = RHODZ2(I,NLEV)
2378     Z0IN(I) = Z0(I)
2379     ZETAIN(I) = ZETA(I)
2380     9118 CONTINUE
2381     DO 9120 I =1,ISTNLV
2382     RIINIT(I,1) = RI(I,1)
2383     QQINIT(I,1) = QQ(I,1)
2384     9120 CONTINUE
2385     ENDIF
2386     C UPDATE TURBULENT KINETIC ENERGY QQ
2387     C ----------------------------------
2388     NLEVMQ = NLEV - LMINQ1
2389     ISTNMQ = irun * NLEVMQ
2390     DO 9306 I =1,ISTNMQ
2391     RHOKDZ(I,LMINQ1) = RHODZ1(I,LMINQ1)
2392     1 * QXLM(I,LMINQ1)
2393     9306 CONTINUE
2394     CALL TRBDIF(QQ(1,LMINQ1),P3(1,LMINQ1),RHOKDZ(1,LMINQ1),
2395 molod 1.15 1 FLXFCE(1,LMINQ1),DTHS,DELTHS,NLEVMQ,1,1.0 _d -20,irun)
2396 molod 1.1 C
2397     C SAVE KH BEFORE ADDING DIMENSIONS FOR USE BY MOIST BOUYANCY CALCULATION
2398     C
2399     DO I = 1,ISTNM1
2400     KHSAVE(I,1) = KH(I,1)
2401     ENDDO
2402     C
2403     C DIMENSIONAL DIFFUSION COEFFS INCLUDING BACKGROUND AMOUNTS
2404     C
2405     IF(LMINQ1.GT.1)THEN
2406     ISTLMQ = irun * (LMINQ1-1)
2407     CB
2408     DO 9124 I =1,ISTLMQ
2409     KM(I,1) = KMBG
2410     KH(I,1) = KHBG
2411     9124 CONTINUE
2412     CE
2413     ENDIF
2414     C
2415     CB
2416     DO 9126 I =1,ISTNMQ
2417     Q(I,LMINQ1) = 2. * QQ(I,LMINQ1)
2418     Q(I,LMINQ1) = SQRT(Q(I,LMINQ1))
2419     XQ(I,LMINQ1) = XL(I,LMINQ1) * Q(I,LMINQ1)
2420     KM(I,LMINQ1)=XQ(I,LMINQ1)*KM(I,LMINQ1)+KMBG
2421     KH(I,LMINQ1)=XQ(I,LMINQ1)*KH(I,LMINQ1)+KHBG
2422     9126 CONTINUE
2423     CE
2424     C
2425     C CALCULATE INTERNAL FLUXES AND UPDATE PROGNOSTIC VARIABLES: TH AND S
2426     C
2427     DO 9128 I =1,ISTNLV
2428     TEMP(I,1) = RHOZPK(I,1) * KH(I,1)
2429     9128 CONTINUE
2430     DO 9130 I =1,ISTNLV
2431     DELTH(I,1) = 0.
2432     9130 CONTINUE
2433     DO 9132 I =1,irun
2434     DELTH(I,NLEVP1) = 1.
2435     9132 CONTINUE
2436 molod 1.15 CALL TRBDIF(TH,DELTH,TEMP,FLXFPK,DTHS,DELTHS,NLEV,2,0. _d 0,irun)
2437 molod 1.1 do i = 1,irun
2438     hsturb(i) = -1.* dths(i)
2439     dhsdtc(i) = -1.* delths(i)
2440     enddo
2441     do L = 1,nlev
2442     do i = 1,irun
2443     dthdthg(i,L) = delth(i,L)
2444     enddo
2445     enddo
2446     do L = 1,nlev
2447     do i = 1,irun
2448     transth(i,L) = temp(i,L)
2449     enddo
2450     enddo
2451    
2452     DO 9134 I =1,ISTNLV
2453     RHOKDZ(I,1) = RHODZ2(I,1) * KH(I,1)
2454     9134 CONTINUE
2455     DO 9138 I =1,ISTNLV
2456     DELSH(I,1) = 0.
2457     9138 CONTINUE
2458     DO 9140 I =1,irun
2459     DELSH(I,NLEVP1) = 1.
2460     9140 CONTINUE
2461    
2462 molod 1.15 CALL TRBDIF(SH,DELSH,RHOKDZ,FLXFAC,DTHL,DELTHL,NLEV,
2463     . 2,0. _d 0,irun)
2464 molod 1.1 do i = 1,irun
2465     eturb(i) = -1.* dthl(i)
2466     dedqa(i) = -1.* delthl(i)
2467     enddo
2468     do L = 1,nlev
2469     do i = 1,irun
2470     dshdshg(i,L) = delsh(i,L)
2471     enddo
2472     enddo
2473     do L = 1,nlev
2474     do i = 1,irun
2475     transsh(i,L) = rhokdz(i,L)
2476     enddo
2477     enddo
2478    
2479     C
2480     C Update Tracers Due to Turbulent Diffusion
2481     C
2482     do i = 1,irun
2483     rhokdz(i,nlev) = 0.0
2484     enddo
2485    
2486 molod 1.20 c do nt = 1,ntrace
2487     c do i = 1,irun
2488     c tracers(i,nlev+1,nt) = tracers(i,nlev,nt)
2489     c enddo
2490     c CALL TRBDIF(tracers(1,1,nt),DELSH,RHOKDZ,FLXFAC,DTHL,DELTHL,
2491     c . NLEV,4,0. _d 0,irun)
2492     c enddo
2493 molod 1.1 C
2494     C CALCULATE INTERNAL FLUXES AND UPDATE PROGNOSTIC VARIABLES: U AND V
2495     C
2496     DO 9172 I =1,ISTNLV
2497     RHOKDZ(I,1) = RHODZ2(I,1) * KM(I,1)
2498     9172 CONTINUE
2499 molod 1.15 CALL TRBDIF(U,V,RHOKDZ,FLXFAC,DTHS,DELTHS,NLEV,3,0. _d 0,irun)
2500 molod 1.1 C ( FILL DIAGNOSTIC ARRAYS IF REQUIRED )
2501     DO 9174 I =1,ISTNLV
2502     WU(I,1) = WU(I,1) + RHOKDZ(I,1) * ( U(I,2) - U(I,1) )
2503     9174 CONTINUE
2504     DO 9176 I =1,ISTNLV
2505     WV(I,1) = WV(I,1) + RHOKDZ(I,1) * ( V(I,2) - V(I,1) )
2506     9176 CONTINUE
2507     DO 9300 I = 1,ISTNM1
2508     IF ( QQ(I,1) .GT. QQMIN ) THEN
2509     IBITSTB(I,1) = 1
2510     ELSE
2511     IBITSTB(I,1) = 0
2512     ENDIF
2513     IF( IBITSTB(I,1).EQ.1 ) FREQDG(I,1) = FREQDG(I,1) + aitr
2514     9300 CONTINUE
2515     do i = 1,irun
2516     qqcolmin(i) = qq(i,nlev)*0.1
2517     qqcolmax(i) = qq(i,nlev)
2518     levpbl(i) = nlev
2519     enddo
2520     DO L = nlev-1,1,-1
2521     DO I = 1,irun
2522     IF ( (qq(i,l).gt.qqcolmax(I)).and.(levpbl(i).eq.nlev))then
2523     qqcolmax(i) = qq(i,l)
2524     qqcolmin(i) = 0.1*qqcolmax(I)
2525     endif
2526     if((qq(i,l).lt.qqcolmin(i)).and.(levpbl(i).eq.nlev))
2527     1 levpbl(i)=l
2528     enddo
2529     enddo
2530     do i = 1,irun
2531     lp = levpbl(i)
2532     if(lp.lt.nlev)then
2533     pbldpth(I) = pbldpth(I) + ( (PLE(I,nlev+1)-PLE(I,Lp+2)) +
2534     1 ( (ple(i,lp+2)-ple(i,lp+1))*(qq(i,lp+1)-qqcolmin(i))
2535     2 / (qq(i,lp+1)-qq(i,lp)) ) ) * aitr
2536     else
2537     pbldpth(I) = pbldpth(I) + ( (PLE(I,nlev+1)-PLE(I,2)) +
2538     1 ( (ple(i,2)-ple(i,1))*(qq(i,1)-qqcolmin(i))
2539     2 / qq(i,1) ) ) * aitr
2540     endif
2541     enddo
2542     do i=1,irun
2543     sustar(i) = sustar(i) + aitr*ustar(i)
2544     enddo
2545     do i=1,irun
2546     sz0(i) = sz0(i) + aitr*z0(i)
2547     enddo
2548     DO I =1,ISTNLV
2549     EU(I,1) = EU(I,1) + AITR*KM(I,1)
2550     enddo
2551     DO I =1,ISTNLV
2552     ET(I,1) = ET(I,1) + AITR*KH(I,1)
2553     enddo
2554     DO I =1,irun
2555     scu(I) = scu(I) + AITR*cu(I)
2556     enddo
2557     DO I =1,irun
2558     sct(I) = sct(I) + AITR*ct(I)
2559     enddo
2560     IF(tprof) then
2561     do i=1,ISTNM1
2562     XLDIAG(I,1) = XLDIAG(I,1) + AITR*XL(I,1)
2563     enddo
2564     endif
2565     FIRST = .FALSE.
2566     C
2567     C SAVE XL,CT,XX,YY,ZETA FOR USE BY MOIST BOUYANCY CALCULATION
2568     C
2569     IF(ITER.EQ.ITRTRB)THEN
2570     DO I = 1,ISTNM1
2571     XLSAVE(I,1) = XL(I,1)
2572     ENDDO
2573     DO I = 1,irun
2574     CTSAVE(I) = CT(I)
2575     XXSAVE(I) = XX(I)
2576     YYSAVE(I) = YY(I)
2577     ZETASAVE(I) = ZETA(I)
2578     ENDDO
2579     ENDIF
2580    
2581     do i = 1,istnlv
2582     turbfcc(i,1) = turbfcc(i,1) + fcc(i,1) * aitr
2583     enddo
2584     do i = 1,irun*nlev
2585     qliq(i,1) = qliq(i,1) + shl(i,1) * aitr
2586     enddo
2587     C
2588     C END OF MAIN LOOP
2589     C
2590     2000 CONTINUE
2591     DO 9194 I =1,ISTNLV
2592     WU(I,1) = WU(I,1) * AITR
2593     WV(I,1) = WV(I,1) * AITR
2594     9194 CONTINUE
2595     C
2596     RETURN
2597     END
2598     SUBROUTINE SFCFLX(NN,VUS,VVS,VTHV1,VTHV2,VTH1,VTH2,VSH1,
2599     1 VSH2,VPK,VPKE,VPE,VZ0,IVWATER,VHS,
2600     2 VAHS,FIRST,LAST,N,IRUN,aitr,VRHO,VRHOZPK,VKH,VKM,
2601     3 VUSTAR,VXX,VYY,VCU,VCT,VRIB,VZETA,VWS,
2602     4 stu2m,stv2m,stt2m,stq2m,stu10m,stv10m,stt10m,stq10m,
2603     5 cp,rgas,undef,
2604     6 lwater, ivbitrib,
2605     7 VHZ,VPSIM,VAPSIM,VPSIG,VPSIHG,VTEMP,VDZETA,VDZ0,VDPSIM,
2606     8 VDPSIH,VZH,VXX0,VYY0,VAPSIHG,VRIB1,VWS1,VPSIH,VZETAL,
2607     9 VZ0L,VPSIH2,VH0,
2608     3 VX0PSIM,VG,VG0,VR1MG0,VZ2,VDZSEA,VAZ0,VXNUM1,VPSIGB2,VDX,
2609     4 VDXPSIM,VDY,VXNUM2,VDEN,VAWS1,VXNUM3,VXNUM,VDZETA1,VDZETA2,
2610     5 VZCOEF2,VZCOEF1,VTEMPLIN,VDPSIMC,VDPSIHC)
2611     C**********************************************************************
2612     C SUBROUTINE SFCFLX - COMPUTES SURFACE TRANSFER COEFFICIENTS
2613     C - CALLED FROM TRBFLX
2614     C
2615     C ARGUMENTS ::
2616     C
2617     C INPUT:
2618     C ------
2619     C US - U - COMPONENT OF SURFACE WIND
2620     C VS - V - COMPONENT OF SURFACE WIND
2621     C THV1 - VIRTUAL POTENTIAL TEMPERATURE AT NLAY
2622     C THV2 - VIRTUAL POTENTIAL TEMPERATURE AT GROUND
2623     C TH1 - POTENTIAL TEMPERATURE AT NLAY
2624     C TH2 - POTENTIAL TEMPERATURE AT GROUND
2625     C SH1 - SPECIFIC HUMIDITY AT NLAY
2626     C SH2 - SPECIFIC HUMIDITY AT GROUND
2627     C PK - EVEN LEVEL PRESSURE ** KAPPA AT LEVEL NLAY
2628     C PKE - EDGE LEVEL PRESSURE ** KAPPA AT GROUND
2629     C PE - SURFACE PRESSURE
2630     C Z0 - SURFACE ROUGHNESS
2631     C WATER - BIT ARRAY - '1' OVER OCEANS
2632     C HS - DEPTH OF SURFACE LAYER
2633     C AHS - ONE / HS
2634     C FIRST - LOGICAL .TRUE. FOR FIRST TRBFLX ITERATION
2635     C LAST - LOGICAL .TRUE. FOR LAST TRBFLX ITERATION
2636     C N - NUMBER OF SFCFLX ITERATIONS
2637     C OUTPUT:
2638     C -------
2639     C RHO - DENSITY AT 10M HEIGHT
2640     C RHOZPK - RHO * P**K AT THE SURFACE
2641     C KH - HEAT TRANSFER COEFFICIENT (CT*USTAR)
2642     C KM - MOMENTUM TRANSFER COEFFICIENT (CU*USTAR)
2643     C USTAR - FRICTION VELOCITY
2644     C XX - PHIM(ZETA) - DIMENSIONLESS WIND SHEAR
2645     C YY - PHIH(ZETA) - DIMENSIONLESS TEMP GRADIENT
2646     C CU - MOMENTUM TRANSPORT COEFFICIENT
2647     C CT - HEAT TRANSPORT COEFFICIENT
2648     C
2649     C**********************************************************************
2650 molod 1.11 implicit none
2651    
2652     C Argument List Declarations
2653     integer nn,n,irun
2654 molod 1.14 _RL aitr,cp,rgas,undef
2655     _RL VUS(IRUN),VVS(IRUN),VTHV1(IRUN),VTHV2(IRUN)
2656     _RL VTH1(IRUN),VTH2(IRUN),VSH1(IRUN),VSH2(IRUN)
2657     _RL VPK(IRUN),VPKE(IRUN),VPE(IRUN)
2658     _RL VZ0(IRUN),VHS(IRUN),VAHS(IRUN)
2659 molod 1.11 integer IVWATER(IRUN)
2660     LOGICAL FIRST,LAST
2661 molod 1.14 _RL VRHO(IRUN),VRHOZPK(IRUN)
2662     _RL VKM(IRUN),VKH(IRUN),VUSTAR(IRUN),VXX(IRUN)
2663     _RL VYY(IRUN),VCU(IRUN),VCT(IRUN),VRIB(IRUN)
2664     _RL VZETA(IRUN),VWS(IRUN)
2665     _RL stu2m(irun),stv2m(irun),stt2m(irun),stq2m(irun)
2666     _RL stu10m(irun),stv10m(irun),stt10m(irun),stq10m(irun)
2667 molod 1.11 LOGICAL LWATER
2668     integer IVBITRIB(irun)
2669 molod 1.14 _RL VHZ(irun),VPSIM(irun),VAPSIM(irun),VPSIG(irun),VPSIHG(irun)
2670     _RL VTEMP(irun),VDZETA(irun),VDZ0(irun),VDPSIM(irun)
2671     _RL VDPSIH(irun),VZH(irun),VXX0(irun),VYY0(irun)
2672     _RL VAPSIHG(irun),VRIB1(irun),VWS1(irun)
2673     _RL VPSIH(irun),VZETAL(irun),VZ0L(irun),VPSIH2(irun),VH0(irun)
2674     _RL VX0PSIM(irun),VG(irun),VG0(irun),VR1MG0(irun)
2675     _RL VZ2(irun),VDZSEA(irun),VAZ0(irun),VXNUM1(irun)
2676     _RL VPSIGB2(irun),VDX(irun),VDXPSIM(irun),VDY(irun)
2677     _RL VXNUM2(irun),VDEN(irun),VAWS1(irun),VXNUM3(irun)
2678     _RL VXNUM(irun),VDZETA1(irun),VDZETA2(irun)
2679     _RL VZCOEF2(irun),VZCOEF1(irun),VTEMPLIN(irun)
2680     _RL VDPSIMC(irun),VDPSIHC(irun)
2681 molod 1.1
2682 molod 1.11 C Local Variables
2683 molod 1.14 _RL USTMX3,USTZ0S,Z0MIN,H0BYZ0,USTH0S,H0VEG,Z0VEGM,PRFAC
2684     _RL XPFAC,DIFSQT
2685 molod 1.1 PARAMETER ( USTMX3 = 0.0632456)
2686     PARAMETER ( USTZ0S = 0.2030325E-5)
2687     PARAMETER ( Z0MIN = USTZ0S/USTMX3)
2688     PARAMETER ( H0BYZ0 = 30.0 )
2689     PARAMETER ( USTH0S = H0BYZ0*USTZ0S )
2690     PARAMETER ( H0VEG = 0.01 )
2691     PARAMETER ( Z0VEGM = 0.005 )
2692     PARAMETER ( PRFAC = 0.595864 )
2693     PARAMETER ( XPFAC = .55 )
2694     PARAMETER ( DIFSQT = 3.872983E-3)
2695    
2696 molod 1.14 _RL psihdiag(irun),psimdiag(irun)
2697     _RL getcon,vk,rvk,vk2,bmdl
2698 molod 1.11 integer iwater,itype
2699     integer i,iter
2700 molod 1.1 C
2701     vk = getcon('VON KARMAN')
2702     rvk = 1./vk
2703     vk2 = vk*vk
2704     BMDL = VK * XPFAC * PRFAC / DIFSQT
2705    
2706     C DETERMINE SURFACE WIND MAGNITUDE AND BULK RICHARDSON NUMBER
2707     C
2708     DO 9000 I = 1,IRUN
2709     VWS(I) = VUS(I) * VUS(I) + VVS(I) * VVS(I)
2710     IF ( VWS(I) .LE. 1.e-4) VWS(I) = 1.e-4
2711     VRIB(I) = ( CP * (VPKE(I)-VPK(I)) ) *
2712     1 (VTHV1(I) - VTHV2(I)) / VWS(I)
2713     VWS(I) = SQRT( VWS(I) )
2714     9000 CONTINUE
2715     C
2716     C INITIALIZATION (FIRST TRBFLX ITERATION)
2717     C INITIAL GUESS FOR ROUGHNESS LENGTH Z0 OVER WATER
2718     C
2719     IF (.NOT. FIRST) GO TO 100
2720     C
2721     IWATER = 0
2722     DO 9002 I = 1,IRUN
2723     IF (IVWATER(I).EQ.1) IWATER = IWATER + 1
2724     9002 CONTINUE
2725     LWATER = .FALSE.
2726     IF(IWATER.GE.1)LWATER = .TRUE.
2727     C
2728     IF(LWATER)THEN
2729     DO 9004 I = 1,IRUN
2730     IF (IVWATER(I).EQ.1) VZ0(I) = 0.0003
2731     9004 CONTINUE
2732     ENDIF
2733     do i = 1,irun
2734     vh0(i) = h0byz0 * vz0(i)
2735     if(vz0(i).ge.z0vegm)vh0(i) = h0veg
2736     enddo
2737    
2738     C CU AND PSIHG FOR NEUTRALLY STRATIFIED FLOW
2739     C
2740     DO 9006 I = 1,IRUN
2741     VHZ(I) = VHS(I) / VZ0(I)
2742     VPSIM(I) = LOG( VHZ(I) )
2743     VAPSIM(I) = 1. / VPSIM(I)
2744     VCU(I) = VK * VAPSIM(I)
2745     VUSTAR(I) = VCU(I) * VWS(I)
2746     C
2747     VPSIG(I) = VH0(I) * VUSTAR(I) - USTH0S
2748     if(VPSIG(I).lt.0.) VPSIG(I) = 0.
2749     VPSIG(I) = SQRT( VPSIG(I) )
2750     VPSIG(I) = BMDL * VPSIG(I)
2751     VPSIHG(I) = VPSIM(I) + VPSIG(I)
2752     9006 CONTINUE
2753     C
2754     C LINEAR CORRECTION FOR ERROR IN ROUGHNESS LENGTH Z0
2755     C
2756     IF(LWATER)THEN
2757     DO 9008 I = 1,IRUN
2758     VTEMP(I) = 0.
2759     9008 CONTINUE
2760     CALL LINADJ(NN,VRIB,VRIB,VWS,
2761     1 VWS,VZ0,VUSTAR,IVWATER,
2762     2 VAPSIM,VTEMP,VTEMP,
2763     3 VTEMP,VTEMP,VTEMP,
2764     4 VTEMP,VTEMP,1,.TRUE.,IRUN,VDZETA,
2765     5 VDZ0,VDPSIM,VDPSIH,
2766     6 IVBITRIB,
2767     3 VX0PSIM,VG,VG0,VR1MG0,VZ2,VDZSEA,VAZ0,VXNUM1,VPSIGB2,VDX,
2768     4 VDXPSIM,VDY,VXNUM2,VDEN,VAWS1,VXNUM3,VXNUM,VDZETA1,VDZETA2,
2769     5 VZCOEF2,VZCOEF1,VTEMPLIN,VDPSIMC,VDPSIHC)
2770     DO 9010 I = 1,IRUN
2771     IF ( IVWATER(I).EQ.1 ) THEN
2772     VCU(I) = VCU(I) * (1. - VDPSIM(I)*VAPSIM(I))
2773     VZ0(I) = VZ0(I) + VDZ0(I)
2774     ENDIF
2775     IF ( IVWATER(I).EQ.1) THEN
2776     IF ( VZ0(I) .LE. Z0MIN ) VZ0(I) = Z0MIN
2777     vh0(i) = h0byz0 * vz0(i)
2778     VPSIG(I) = VH0(I) * VCU(I) * VWS(I) - USTH0S
2779     if(VPSIG(I).lt.0.) VPSIG(I) = 0.
2780     VPSIG(I) = SQRT( VPSIG(I) )
2781     VPSIG(I) = BMDL * VPSIG(I)
2782     VPSIHG(I) = VPSIM(I) + VDPSIH(I) + VPSIG(I)
2783     ENDIF
2784     9010 CONTINUE
2785     ENDIF
2786     C
2787     C INITIAL GUESS FOR STABILITY PARAMETER ZETA
2788     C
2789     DO 9012 I = 1,IRUN
2790     VZETA(I) = VK2 * VRIB(I) / (VCU(I) * VCU(I) * VPSIHG(I))
2791     9012 CONTINUE
2792     C
2793     C RECOMPUTE CU, ESTIMATE PSIHG AND UPDATE ZETA AND Z0
2794     C
2795     DO 9014 I = 1,IRUN
2796     VZH(I) = VZ0(I) * VAHS(I)
2797     9014 CONTINUE
2798     CALL PSI (VZETA,VZH,VPSIM,
2799     1 VTEMP,IRUN,VXX,VXX0,VYY,
2800     2 VYY0,2)
2801     DO 9016 I = 1,IRUN
2802     VCU(I) = VK / VPSIM(I)
2803     VPSIG(I) = VH0(I) * VCU(I) * VWS(I) - USTH0S
2804     if(VPSIG(I).lt.0.) VPSIG(I) = 0.
2805     VPSIG(I) = SQRT(VPSIG(I))
2806     VPSIG(I) = BMDL * VPSIG(I)
2807     VPSIHG(I) = VPSIM(I) + VPSIG(I)
2808     VZETA(I) = VK2 * VRIB(I) / (VCU(I) * VCU(I) * VPSIHG(I))
2809     9016 CONTINUE
2810     C
2811     IF(LWATER)THEN
2812     CCCOOOMMMM ADDED 'WHERE WATER'
2813     DO 9018 I = 1,IRUN
2814     IF (IVWATER(I).EQ.1) VUSTAR(I) = VCU(I) * VWS(I)
2815     9018 CONTINUE
2816     CALL ZCSUB ( VUSTAR,VHZ,IVWATER,.FALSE.,IRUN,VTEMP)
2817     DO 9020 I = 1,IRUN
2818     IF (IVWATER(I).EQ.1 ) then
2819     VZ0(I) = VTEMP(I)
2820     IF ( VZ0(I) .LE. Z0MIN ) VZ0(I) = Z0MIN
2821     vh0(i) = h0byz0 * vz0(i)
2822     endif
2823     9020 CONTINUE
2824     ENDIF
2825     C
2826     GO TO 125
2827     C
2828     C LINEARLY UPDATE ZETA AND Z0 FOR SECOND OR GREATER TRBFLX ITERATION
2829     C
2830     100 CONTINUE
2831    
2832     CALL LINADJ(NN,VRIB1,VRIB,VWS1,
2833     1 VWS,VZ0,VUSTAR,IVWATER,
2834     2 VAPSIM,VAPSIHG,VPSIH,
2835     3 VPSIG,VXX,VXX0,
2836     4 VYY,VYY0,2,LWATER,IRUN,VDZETA,
2837     5 VDZ0,VDPSIM,VDPSIH,
2838     6 IVBITRIB,
2839     3 VX0PSIM,VG,VG0,VR1MG0,VZ2,VDZSEA,VAZ0,VXNUM1,VPSIGB2,VDX,
2840     4 VDXPSIM,VDY,VXNUM2,VDEN,VAWS1,VXNUM3,VXNUM,VDZETA1,VDZETA2,
2841     5 VZCOEF2,VZCOEF1,VTEMPLIN,VDPSIMC,VDPSIHC)
2842     C
2843     DO 9022 I = 1,IRUN
2844     VZETA(I) = VZETA(I) + VZETAL(I) * VDZETA(I)
2845     IF (IVBITRIB(I).EQ.1 )VZETA(I) =
2846     1 VPSIM(I) * VPSIM(I) * VRIB(I) * VCT(I) * RVK
2847     9022 CONTINUE
2848     C
2849     IF ( LWATER ) THEN
2850     DO 9024 I = 1,IRUN
2851     IF (IVWATER(I).EQ.1 ) then
2852     VZ0(I) = VZ0(I) + VZ0L(I) * VDZ0(I)
2853     IF (VZ0(I) .LE. Z0MIN ) VZ0(I) = Z0MIN
2854     vh0(i) = h0byz0 * vz0(i)
2855     endif
2856     9024 CONTINUE
2857     ENDIF
2858     C
2859     125 CONTINUE
2860     C
2861     C ITERATIVE LOOP - N ITERATIONS
2862     C COMPUTE CU AND CT
2863     C
2864     DO 200 ITER = 1,N
2865     DO 9026 I = 1,IRUN
2866     VZH(I) = VZ0(I) * VAHS(I)
2867     9026 CONTINUE
2868     CALL PSI (VZETA,VZH,VPSIM,
2869     1 VPSIH,IRUN,VXX,VXX0,VYY,
2870     2 VYY0,1)
2871     DO 9028 I = 1,IRUN
2872     VCU(I) = VK / VPSIM(I)
2873     VUSTAR(I) = VCU(I) * VWS(I)
2874     C
2875     VPSIG(I) = VH0(I) * VUSTAR(I) - USTH0S
2876     if(VPSIG(I).lt.0.) VPSIG(I) = 0.
2877     VPSIG(I) = SQRT(VPSIG(I))
2878     VPSIG(I) = BMDL * VPSIG(I)
2879     VPSIHG(I) = VPSIH(I) + VPSIG(I)
2880     C
2881     C LINEAR CORRECTIONS FOR CU, CT, ZETA, AND Z0
2882     C
2883     VAPSIM(I) = VCU(I) * RVK
2884     VAPSIHG(I) = 1. / VPSIHG(I)
2885     VRIB1(I) = VAPSIM(I) * VAPSIM(I) * VPSIHG(I) * VZETA(I)
2886     9028 CONTINUE
2887     C
2888     ITYPE = 3
2889     IF(ITER.EQ.N) ITYPE = 4
2890     IF( (ITYPE.EQ.4) .AND. (.NOT.LAST) ) ITYPE = 5
2891     C
2892     CALL LINADJ(NN,VRIB1,VRIB,VWS,
2893     1 VWS,VZ0,VUSTAR,IVWATER,
2894     2 VAPSIM,VAPSIHG,VPSIH,
2895     3 VPSIG,VXX,VXX0,
2896     4 VYY,VYY0,ITYPE,LWATER,IRUN,VDZETA,
2897     5 VDZ0,VDPSIM,VDPSIH,
2898     6 IVBITRIB,
2899     3 VX0PSIM,VG,VG0,VR1MG0,VZ2,VDZSEA,VAZ0,VXNUM1,VPSIGB2,VDX,
2900     4 VDXPSIM,VDY,VXNUM2,VDEN,VAWS1,VXNUM3,VXNUM,VDZETA1,VDZETA2,
2901     5 VZCOEF2,VZCOEF1,VTEMPLIN,VDPSIMC,VDPSIHC)
2902     C
2903     C UPDATES OF ZETA, Z0, CU AND CT
2904     C
2905     IF (ITYPE.EQ.5) THEN
2906     DO 9030 I = 1,IRUN
2907     VZETAL(I) = VZETA(I)
2908     VZ0L(I) = VZ0(I)
2909     9030 CONTINUE
2910     ENDIF
2911     C
2912     DO 9032 I = 1,IRUN
2913     VZETA(I) = VZETA(I) * ( 1. + VDZETA(I) )
2914     IF (IVBITRIB(I).EQ.1 ) VZETA(I) =
2915     1 VPSIM(I) * VPSIM(I) * VRIB(I) * VAPSIHG(I)
2916     9032 CONTINUE
2917     C
2918     IF ( LWATER ) THEN
2919     DO 9034 I = 1,IRUN
2920     IF (IVWATER(I).EQ.1 ) then
2921     VZ0(I) = VZ0(I) * ( 1. + VDZ0(I) )
2922     IF (VZ0(I) .LE. Z0MIN ) VZ0(I) = Z0MIN
2923     vh0(i) = h0byz0 * vz0(i)
2924     endif
2925     9034 CONTINUE
2926     ENDIF
2927     C
2928     IF ( ITER .EQ. N ) THEN
2929     DO 9036 I = 1,IRUN
2930     VPSIM(I) = VPSIM(I) + VDPSIM(I)
2931     VCU(I) = VK / VPSIM(I)
2932     VUSTAR(I) = VCU(I) * VWS(I)
2933     C
2934     VPSIG(I) = VH0(I) * VUSTAR(I) - USTH0S
2935     if(VPSIG(I).lt.0.) VPSIG(I) = 0.
2936     VPSIG(I) = SQRT(VPSIG(I))
2937     VPSIG(I) = BMDL * VPSIG(I)
2938     VPSIHG(I) = VPSIH(I) + VDPSIH(I) + VPSIG(I)
2939     VCT(I) = VK / VPSIHG(I)
2940     9036 CONTINUE
2941     ENDIF
2942     C
2943     C SAVE VALUES OF RIB AND WS FOR NEXT ITERATION OF TRBFLX
2944     C
2945     IF (ITYPE.EQ.5) THEN
2946     DO 9038 I = 1,IRUN
2947     VRIB1(I) = VRIB(I)
2948     VWS1(I) = VWS(I)
2949     9038 CONTINUE
2950     ENDIF
2951     C
2952     200 CONTINUE
2953     C
2954     C CALCULATE RHO-SURFACE ( KG / M**3 )
2955     C
2956     IF (FIRST) THEN
2957     DO I = 1,IRUN
2958     VTEMP(I) = 10. * VAHS(I) * VZETA(I)
2959     VZH(I) = VZ0(I) * 0.1
2960     ENDDO
2961     CALL PSI (VTEMP,VZH,VHZ,
2962     1 VPSIH2,IRUN,VHZ,VHZ,VHZ,
2963     2 VHZ,3)
2964     DO I = 1,IRUN
2965     VTEMP(I) = ( VPSIH2(I) + VPSIG(I) ) / VPSIHG(I)
2966     VRHO(I) = VPKE(I)*( VTH2(I) + VTEMP(I) * (VTH1(I)-VTH2(I)) )
2967     VRHO(I) = VPE(I)*100. / ( RGAS * VRHO(I) )
2968     ENDDO
2969     ENDIF
2970     C
2971     C interpolate uvtq to 2m and to 10 meters for diagnostic output
2972     C use psih and psim which represent non-dim change from ground
2973     C to specified level
2974     C and multiply theta by surface p**kappa to get temperatures
2975     C
2976     do i = 1,irun
2977     vtemp(i) = 2. * vahs(i) * vzeta(i)
2978     vzh(i) = vz0(i) * 0.5
2979     if(vz0(i).ge.2.)vzh(i) = 0.9
2980     enddo
2981     call psi(vtemp,vzh,psimdiag,psihdiag,irun,vhz,vhz,vhz,vhz,1)
2982     do i = 1,irun
2983     stu2m(i) = (psimdiag(i)/vpsim(i) * vus(i))
2984     stv2m(i) = (psimdiag(i)/vpsim(i) * vvs(i))
2985     stt2m(i) = ( (vth2(i) + ((psihdiag(i)+vpsig(i))/vpsihg(i))*
2986     1 (vth1(i)-vth2(i))) ) * vpke(i)
2987     stq2m(i) = (vsh2(i) + ((psihdiag(i)+vpsig(i))/vpsihg(i))*
2988     1 (vsh1(i)-vsh2(i)))
2989     if(vz0(i).ge.2.)then
2990     stu2m(i) = UNDEF
2991     stv2m(i) = UNDEF
2992     stt2m(i) = UNDEF
2993     stq2m(i) = UNDEF
2994     endif
2995     enddo
2996     do i = 1,irun
2997     vtemp(i) = 10. * vahs(i) * vzeta(i)
2998     vzh(i) = vz0(i) * 0.1
2999     enddo
3000     call psi(vtemp,vzh,psimdiag,psihdiag,irun,vhz,vhz,vhz,vhz,1)
3001     do i = 1,irun
3002     stu10m(i) = (psimdiag(i)/vpsim(i) * vus(i))
3003     stv10m(i) = (psimdiag(i)/vpsim(i) * vvs(i))
3004     stt10m(i) = ( (vth2(i) + ((psihdiag(i)+vpsig(i))/vpsihg(i))*
3005     1 (vth1(i)-vth2(i))) ) * vpke(i)
3006     stq10m(i) = (vsh2(i) + ((psihdiag(i)+vpsig(i))/vpsihg(i))*
3007     1 (vsh1(i)-vsh2(i)))
3008     enddo
3009     C
3010     C EVALUATE TURBULENT TRANSFER COEFFICIENTS
3011     C
3012     DO 9044 I = 1,IRUN
3013     VRHOZPK(I) = VRHO(I) * VPKE(I)
3014     VKH(I) = VUSTAR(I) * VCT(I)
3015     VKM(I) = VUSTAR(I) * VCU(I)
3016     9044 CONTINUE
3017     C
3018     RETURN
3019     END
3020     SUBROUTINE PHI(Z,PHIM,PHIH,IFLAG,N)
3021     C**********************************************************************
3022     C
3023     C FUNCTION PHI - SOLVES KEYPS EQUATIONS
3024     C - CALLED FROM PSI
3025     C
3026     C DESCRIPTION OF PARAMETERS
3027     C Z - INPUTED VALUE OF MONIN- OBUKHOV STABILITY PARAMETER ZETA
3028     C TIMES APPROPRIATE CONSTANT
3029     C PHIM - OUTPUTED SOLUTION OF KEYPS EQUATION FOR MOMENTUM
3030     C PHIH - OUTPUTED SOLUTION OF KEYPS EQUATION FOR SCALARS
3031     C IFLAG - FLAG TO DETERMINE IF X IS NEEDED (IFLAG=2), Y IS NEEDED
3032     C (IFLAG=3), OR BOTH (IFLAG=1)
3033     C N - LENGTH OF VECTOR TO BE SOLVED
3034     C
3035     C**********************************************************************
3036 molod 1.11 implicit none
3037    
3038     C Argument List Declarations
3039     integer n,iflag
3040 molod 1.14 _RL PHIM(N),PHIH(N),Z(N)
3041 molod 1.11
3042     C Local Variables
3043     integer I1(N),I2(N)
3044 molod 1.14 _RL ZSTAR(N),E1(N),E2(N),TEMP1(N)
3045 molod 1.1 C
3046 molod 1.14 _RL PHIM0(385),ZLINM1(75),ZLINM2(75),ZLINM3(36)
3047     _RL ZLOGM1(74),ZLOGM2(75),ZLOGM3(50)
3048     _RL PHIH0(385),ZLINH1(75),ZLINH2(75),ZLINH3(36)
3049     _RL ZLOGH1(74),ZLOGH2(75),ZLOGH3(50)
3050 molod 1.1 EQUIVALENCE (PHIM0(1),ZLINM1(1)),(PHIM0(76),ZLINM2(1))
3051     EQUIVALENCE (PHIM0(151),ZLINM3(1))
3052     EQUIVALENCE (PHIM0(187),ZLOGM1(1)),(PHIM0(261),ZLOGM2(1))
3053     EQUIVALENCE (PHIM0(336),ZLOGM3(1))
3054     EQUIVALENCE (PHIH0(1),ZLINH1(1)),(PHIH0(76),ZLINH2(1))
3055     EQUIVALENCE (PHIH0(151),ZLINH3(1))
3056     EQUIVALENCE (PHIH0(187),ZLOGH1(1)),(PHIH0(261),ZLOGH2(1))
3057     EQUIVALENCE (PHIH0(336),ZLOGH3(1))
3058     C
3059     DATA ZLOGM1/
3060     . 0.697894,0.678839,0.659598,0.640260,
3061     . 0.620910,0.601628,0.582486,0.563550,0.544877,
3062     . 0.526519,0.508516,0.490903,0.473708,0.456951,
3063     . 0.440649,0.424812,0.409446,0.394553,0.380133,
3064     . 0.366182,0.352695,0.339664,0.327082,0.314938,
3065     . 0.303222,0.291923,0.281029,0.270528,0.260409,
3066     . 0.250659,0.241267,0.232221,0.223509,0.215119,
3067     . 0.207041,0.199264,0.191776,0.184568,0.177628,
3068     . 0.170949,0.164519,0.158331,0.152374,0.146641,
3069     . 0.141123,0.135813,0.130702,0.125783,0.121048,
3070     . 0.116492,0.112107,0.107887,0.103826,0.0999177,
3071     . 0.0961563,0.0925364,0.0890528,0.0857003,0.0824739,
3072     . 0.0793690,0.0763810,0.0735054,0.0707380,0.0680749,
3073     . 0.0655120,0.0630455,0.0606720,0.0583877,0.0561895,
3074     . 0.0540740,0.0520382,0.0500790,0.0481936,0.0463791/
3075     DATA ZLOGM2/
3076     . 0.0446330,0.0429526,0.0413355,0.0397792,0.0382816,
3077     . 0.0368403,0.0354533,0.0341185,0.0328340,0.0315978,
3078     . 0.0304081,0.0292633,0.0281616,0.0271013,0.0260809,
3079     . 0.0250990,0.0241540,0.0232447,0.0223695,0.0215273,
3080     . 0.0207168,0.0199369,0.0191862,0.0184639,0.0177687,
3081     . 0.0170998,0.0164560,0.0158364,0.0152402,0.0146664,
3082     . 0.0141142,0.0135828,0.0130714,0.0125793,0.0121057,
3083     . 0.0116499,0.0112113,0.0107892,0.0103830,0.999210E-2,
3084     . 0.961590E-2,0.925387E-2,0.890547E-2,0.857018E-2,0.824752E-2,
3085     . 0.793701E-2,0.763818E-2,0.735061E-2,0.707386E-2,0.680754E-2,
3086     . 0.655124E-2,0.630459E-2,0.606722E-2,0.583880E-2,0.561897E-2,
3087     . 0.540742E-2,0.520383E-2,0.500791E-2,0.481937E-2,0.463792E-2,
3088     . 0.446331E-2,0.429527E-2,0.413355E-2,0.397793E-2,0.382816E-2,
3089     . 0.368403E-2,0.354533E-2,0.341185E-2,0.328340E-2,0.315978E-2,
3090     . 0.304082E-2,0.292633E-2,0.281616E-2,0.271013E-2,0.260809E-2/
3091     DATA ZLOGM3/
3092     . 0.250990E-2,0.241541E-2,0.232447E-2,0.223695E-2,0.215273E-2,
3093     . 0.207168E-2,0.199369E-2,0.191862E-2,0.184639E-2,0.177687E-2,
3094     . 0.170998E-2,0.164560E-2,0.158364E-2,0.152402E-2,0.146664E-2,
3095     . 0.141142E-2,0.135828E-2,0.130714E-2,0.125793E-2,0.121057E-2,
3096     . 0.116499E-2,0.112113E-2,0.107892E-2,0.103830E-2,0.999210E-3,
3097     . 0.961590E-3,0.925387E-3,0.890547E-3,0.857018E-3,0.824752E-3,
3098     . 0.793701E-3,0.763818E-3,0.735061E-3,0.707386E-3,0.680754E-3,
3099     . 0.655124E-3,0.630459E-3,0.606722E-3,0.583880E-3,0.561897E-3,
3100     . 0.540742E-3,0.520383E-3,0.500791E-3,0.481937E-3,0.463792E-3,
3101     . 0.446331E-3,0.429527E-3,0.413355E-3,0.397793E-3,0.382816E-3/
3102     DATA ZLOGH1/
3103     . 0.640529,0.623728,0.606937,0.590199,
3104     . 0.573552,0.557032,0.540672,0.524504,0.508553,
3105     . 0.492843,0.477397,0.462232,0.447365,0.432809,
3106     . 0.418574,0.404670,0.391103,0.377878,0.364999,
3107     . 0.352468,0.340284,0.328447,0.316954,0.305804,
3108     . 0.294992,0.284514,0.274364,0.264538,0.255028,
3109     . 0.245829,0.236933,0.228335,0.220026,0.211999,
3110     . 0.204247,0.196762,0.189537,0.182564,0.175837,
3111     . 0.169347,0.163088,0.157051,0.151231,0.145620,
3112     . 0.140211,0.134998,0.129974,0.125133,0.120469,
3113     . 0.115975,0.111645,0.107475,0.103458,0.995895E-1,
3114     . 0.958635E-1,0.922753E-1,0.888199E-1,0.854925E-1,0.822886E-1,
3115     . 0.792037E-1,0.762336E-1,0.733739E-1,0.706208E-1,0.679704E-1,
3116     . 0.654188E-1,0.629625E-1,0.605979E-1,0.583217E-1,0.561306E-1,
3117     . 0.540215E-1,0.519914E-1,0.500373E-1,0.481564E-1,0.463460E-1/
3118     DATA ZLOGH2/
3119     . 0.446034E-1,0.429263E-1,0.413120E-1,0.397583E-1,0.382629E-1,
3120     . 0.368237E-1,0.354385E-1,0.341053E-1,0.328222E-1,0.315873E-1,
3121     . 0.303988E-1,0.292550E-1,0.281541E-1,0.270947E-1,0.260750E-1,
3122     . 0.250937E-1,0.241494E-1,0.232405E-1,0.223658E-1,0.215240E-1,
3123     . 0.207139E-1,0.199342E-1,0.191839E-1,0.184618E-1,0.177669E-1,
3124     . 0.170981E-1,0.164545E-1,0.158351E-1,0.152390E-1,0.146653E-1,
3125     . 0.141133E-1,0.135820E-1,0.130707E-1,0.125786E-1,0.121051E-1,
3126     . 0.116494E-1,0.112108E-1,0.107888E-1,0.103826E-1,0.999177E-2,
3127     . 0.961561E-2,0.925360E-2,0.890523E-2,0.856997E-2,0.824733E-2,
3128     . 0.793684E-2,0.763803E-2,0.735048E-2,0.707375E-2,0.680743E-2,
3129     . 0.655114E-2,0.630450E-2,0.606715E-2,0.583873E-2,0.561891E-2,
3130     . 0.540737E-2,0.520379E-2,0.500787E-2,0.481933E-2,0.463789E-2,
3131     . 0.446328E-2,0.429524E-2,0.413353E-2,0.397790E-2,0.382814E-2,
3132     . 0.368401E-2,0.354532E-2,0.341184E-2,0.328338E-2,0.315977E-2,
3133     . 0.304081E-2,0.292632E-2,0.281615E-2,0.271012E-2,0.260809E-2/
3134     DATA ZLOGH3/
3135     . 0.250990E-2,0.241540E-2,0.232446E-2,0.223695E-2,0.215273E-2,
3136     . 0.207168E-2,0.199368E-2,0.191862E-2,0.184639E-2,0.177687E-2,
3137     . 0.170997E-2,0.164559E-2,0.158364E-2,0.152402E-2,0.146664E-2,
3138     . 0.141142E-2,0.135828E-2,0.130714E-2,0.125793E-2,0.121057E-2,
3139     . 0.116499E-2,0.112113E-2,0.107892E-2,0.103830E-2,0.999209E-3,
3140     . 0.961590E-3,0.925387E-3,0.890546E-3,0.857018E-3,0.824752E-3,
3141     . 0.793700E-3,0.763818E-3,0.735061E-3,0.707386E-3,0.680754E-3,
3142     . 0.655124E-3,0.630459E-3,0.606722E-3,0.583880E-3,0.561897E-3,
3143     . 0.540742E-3,0.520383E-3,0.500791E-3,0.481937E-3,0.463792E-3,
3144     . 0.446331E-3,0.429527E-3,0.413355E-3,0.397793E-3,0.382816E-3/
3145     C
3146     DATA ZLINM1/
3147     & 0.964508,0.962277,0.960062,0.957863,0.955680,
3148     & 0.953512,0.951359,0.949222,0.947100,0.944992,
3149     & 0.942899,0.940821,0.938758,0.936709,0.934673,
3150     & 0.932652,0.930645,0.928652,0.926672,0.924706,
3151     & 0.922753,0.920813,0.918886,0.916973,0.915072,
3152     & 0.913184,0.911308,0.909445,0.907594,0.905756,
3153     & 0.903930,0.902115,0.900313,0.898522,0.896743,
3154     & 0.894975,0.893219,0.891475,0.889741,0.888019,
3155     & 0.886307,0.884607,0.882917,0.881238,0.879569,
3156     & 0.877911,0.876264,0.874626,0.872999,0.871382,
3157     & 0.869775,0.868178,0.866591,0.865013,0.863445,
3158     & 0.861887,0.860338,0.858798,0.857268,0.855747,
3159     & 0.854235,0.852732,0.851238,0.849753,0.848277,
3160     & 0.846809,0.845350,0.843900,0.842458,0.841025,
3161     & 0.839599,0.838182,0.836774,0.835373,0.833980/
3162     DATA ZLINM2/
3163     & 0.832596,0.831219,0.829850,0.828489,0.827136,
3164     & 0.825790,0.824451,0.823121,0.821797,0.820481,
3165     & 0.819173,0.817871,0.816577,0.815289,0.814009,
3166     & 0.812736,0.811470,0.810210,0.808958,0.807712,
3167     & 0.806473,0.805240,0.804015,0.802795,0.801582,
3168     & 0.800376,0.799176,0.797982,0.796794,0.795613,
3169     & 0.794438,0.793269,0.792106,0.790949,0.789798,
3170     & 0.788652,0.787513,0.786380,0.785252,0.784130,
3171     & 0.783014,0.781903,0.780798,0.779698,0.778604,
3172     & 0.777516,0.776432,0.775354,0.774282,0.773215,
3173     & 0.772153,0.771096,0.770044,0.768998,0.767956,
3174     & 0.766920,0.765888,0.764862,0.763840,0.762824,
3175     & 0.761812,0.760805,0.759803,0.758805,0.757813,
3176     & 0.756824,0.755841,0.754862,0.753888,0.752918,
3177     & 0.751953,0.750992,0.750035,0.749083,0.748136/
3178     DATA ZLINM3/
3179     & 0.747192,0.746253,0.745318,0.744388,0.743462,
3180     & 0.742539,0.741621,0.740707,0.739798,0.738892,
3181     & 0.737990,0.737092,0.736198,0.735308,0.734423,
3182     & 0.733540,0.732662,0.731788,0.730917,0.730050,
3183     & 0.729187,0.728328,0.727472,0.726620,0.725772,
3184     & 0.724927,0.724086,0.723248,0.722414,0.721584,
3185     & 0.720757,0.719933,0.719113,0.718296,0.717483,
3186     & 0.716673/
3187     DATA ZLINH1/
3188     & 0.936397,0.932809,0.929287,0.925827,0.922429,
3189     & 0.919089,0.915806,0.912579,0.909405,0.906284,
3190     & 0.903212,0.900189,0.897214,0.894284,0.891399,
3191     & 0.888558,0.885759,0.883001,0.880283,0.877603,
3192     & 0.874962,0.872357,0.869788,0.867255,0.864755,
3193     & 0.862288,0.859854,0.857452,0.855081,0.852739,
3194     & 0.850427,0.848144,0.845889,0.843662,0.841461,
3195     & 0.839287,0.837138,0.835014,0.832915,0.830841,
3196     & 0.828789,0.826761,0.824755,0.822772,0.820810,
3197     & 0.818869,0.816949,0.815050,0.813170,0.811310,
3198     & 0.809470,0.807648,0.805845,0.804060,0.802293,
3199     & 0.800543,0.798811,0.797095,0.795396,0.793714,
3200     & 0.792047,0.790396,0.788761,0.787141,0.785535,
3201     & 0.783945,0.782369,0.780807,0.779259,0.777724,
3202     & 0.776204,0.774696,0.773202,0.771720,0.770251/
3203     DATA ZLINH2/
3204     & 0.768795,0.767351,0.765919,0.764499,0.763091,
3205     & 0.761694,0.760309,0.758935,0.757571,0.756219,
3206     & 0.754878,0.753547,0.752226,0.750916,0.749616,
3207     & 0.748326,0.747045,0.745775,0.744514,0.743262,
3208     & 0.742020,0.740787,0.739563,0.738348,0.737141,
3209     & 0.735944,0.734755,0.733574,0.732402,0.731238,
3210     & 0.730083,0.728935,0.727795,0.726664,0.725539,
3211     & 0.724423,0.723314,0.722213,0.721119,0.720032,
3212     & 0.718952,0.717880,0.716815,0.715756,0.714704,
3213     & 0.713660,0.712621,0.711590,0.710565,0.709547,
3214     & 0.708534,0.707529,0.706529,0.705536,0.704549,
3215     & 0.703567,0.702592,0.701623,0.700660,0.699702,
3216     & 0.698750,0.697804,0.696863,0.695928,0.694998,
3217     & 0.694074,0.693155,0.692241,0.691333,0.690430,
3218     & 0.689532,0.688639,0.687751,0.686868,0.685990/
3219     DATA ZLINH3/
3220     & 0.685117,0.684249,0.683386,0.682527,0.681673,
3221     & 0.680824,0.679979,0.679139,0.678303,0.677472,
3222     & 0.676645,0.675823,0.675005,0.674191,0.673381,
3223     & 0.672576,0.671775,0.670978,0.670185,0.669396,
3224     & 0.668611,0.667830,0.667054,0.666281,0.665512,
3225     & 0.664746,0.663985,0.663227,0.662473,0.661723,
3226     & 0.660977,0.660234,0.659495,0.658759,0.658027,
3227     & 0.657298/
3228 molod 1.11
3229     integer ibit1,ibit2,i
3230 molod 1.1 C
3231     IBIT1 = 0
3232     IBIT2 = 0
3233     C
3234     DO 9000 I = 1,N
3235     IF(Z(I).GE.0.15)IBIT1 = IBIT1 + 1
3236     IF(Z(I).GT.2. )IBIT2 = IBIT2 + 1
3237     9000 CONTINUE
3238     C
3239     IF( IBIT1 .LE. 0 ) GO TO 200
3240     C
3241     DO 9002 I = 1,N
3242     ZSTAR(I) = 100. * Z(I) - 14.
3243     9002 CONTINUE
3244     C
3245     IF( IBIT2 .LE. 0 ) GO TO 60
3246     DO 9004 I = 1,N
3247     TEMP1(I) = Z(I)*0.5
3248     IF( Z(I) .LE. 2. )TEMP1(I) = 1.
3249     TEMP1(I) = LOG10(TEMP1(I))
3250     TEMP1(I) = (TEMP1(I) + 9.3) * 20.
3251     IF( Z(I) .GT. 2. ) ZSTAR(I) = TEMP1(I)
3252     IF( Z(I).GT.1.78e10 ) ZSTAR(I) = 384.9999
3253     9004 CONTINUE
3254     C
3255     60 CONTINUE
3256     C
3257     DO 9006 I = 1,N
3258     I1(I) = ZSTAR(I)
3259     I2(I) = I1(I) + 1
3260     TEMP1(I) = ZSTAR(I) - I1(I)
3261     C
3262     9006 CONTINUE
3263     C
3264     IF( IFLAG .GT. 2 ) GO TO 100
3265     DO 9008 I = 1,N
3266     if( z(i).ge.0.15 ) then
3267     E1(I) = PHIM0( I1(I) )
3268     E2(I) = PHIM0( I2(I) )
3269     PHIM(I) = TEMP1(I) * ( E2(I)-E1(I) )
3270     PHIM(I) = PHIM(I) + E1(I)
3271     endif
3272     9008 CONTINUE
3273    
3274     100 CONTINUE
3275     C
3276     IF( IFLAG .EQ. 2 ) GO TO 200
3277     DO 9010 I = 1,N
3278     if( z(i).ge.0.15 ) then
3279     E1(I) = PHIH0( I1(I) )
3280     E2(I) = PHIH0( I2(I) )
3281     PHIH(I) = TEMP1(I) * ( E2(I)-E1(I) )
3282     PHIH(I) = PHIH(I) + E1(I)
3283     endif
3284     9010 CONTINUE
3285    
3286     200 CONTINUE
3287     IF( IBIT1 .GE. N ) GO TO 500
3288     C
3289     DO 9012 I = 1,N
3290     ZSTAR(I) = -Z(I)
3291     9012 CONTINUE
3292     C
3293     IF( IFLAG .GT. 2 ) GO TO 300
3294     DO 9014 I = 1,N
3295     IF( Z(I) .LT. 0.15 ) PHIM(I) = 1. + ZSTAR(I)
3296     2 *(0.25+ZSTAR(I)*(0.09375+ZSTAR(I)*
3297     3 (0.03125+0.00732422 * ZSTAR(I))))
3298     9014 CONTINUE
3299     C
3300     300 CONTINUE
3301     IF( IFLAG .EQ. 2 ) GO TO 500
3302     DO 9016 I = 1,N
3303     IF( Z(I) .LT. 0.15 ) THEN
3304     PHIH(I) =1.+ Z(I) * (0.5+ZSTAR(I)*(0.375+ZSTAR(I)*
3305     1 (0.5+ZSTAR(I)*(0.8203125+ZSTAR(I)*
3306     2 (1.5+2.93262*ZSTAR(I))))))
3307     PHIH(I) = 1. / PHIH(I)
3308     ENDIF
3309     9016 CONTINUE
3310     C
3311     500 CONTINUE
3312     RETURN
3313     END
3314     SUBROUTINE PSI(VZZ,VZH,VPSIM,VPSIH,IRUN,VX,VXS,VY,VYS,IFLAG)
3315     C**********************************************************************
3316     C
3317     C SUBROUTINE PSI - DETERMINES DIMENSIONLESS WIND AND
3318     C SCALAR PROFILES IN SURFACE LAYER
3319     C - CALLED FROM SFCFLX
3320     C
3321     C DESCRIPTION OF PARAMETERS
3322     C ZZ - INPUTED VALUE OF MONIN- OBUKHOV STABILITY PARAMETER ZETA
3323     C ZH - INPUTED VALUE OF PBL HEIGHT DIVIDED BY Z0
3324     C PSIM - OUTPUTED VALUE OF DIMENSIONLESS WIND
3325     C PSIH - OUTPUTED VALUE OF DIMENSIONLESS SCALAR
3326     C X - OUTPUTED VALUE OF PHIM(ZETA)
3327     C XS - OUTPUTED VALUE OF PHIM(ZETA0)
3328     C Y - OUTPUTED VALUE OF PHIH(ZETA)
3329     C YS - OUTPUTED VALUE OF PHIH(ZETA0)
3330     C IFLAG- FLAG TO DETERMINE IF CU IS NEEDED (IFLAG=2),
3331     C IF CT IS NEEDED (IFLAG=3), OR BOTH (IFLAG=1)
3332     C SUBPROGRAMS NEEDED
3333     C PHI - COMPUTES SIMILARITY FUNCTION FOR MOMENTUM AND SCALARS
3334     C
3335     C**********************************************************************
3336 molod 1.11 implicit none
3337    
3338     C Argument List Declarations
3339     integer irun,iflag
3340 molod 1.14 _RL VZZ(IRUN),VZH(IRUN),VPSIM(IRUN),VPSIH(IRUN),
3341 molod 1.11 1 VX(IRUN),VXS(IRUN),VY(IRUN),VYS(IRUN)
3342    
3343     C Local Variables
3344 molod 1.14 _RL ZWM,RZWM,Z0M,ZCM,RZCM,CM1,CM2,CM6,CM7,CM8ARG,YCM
3345 molod 1.1 PARAMETER ( ZWM = 1. )
3346     PARAMETER ( RZWM = 1./ZWM )
3347     PARAMETER ( Z0M = 0.2 )
3348     PARAMETER ( ZCM = 42. )
3349     PARAMETER ( RZCM = 1./ZCM )
3350     PARAMETER ( CM1 = 1./126. )
3351     PARAMETER ( CM2 = 1./(6.*CM1) )
3352     PARAMETER ( CM6 = 6. / ( 1. + 6.*CM1 ) )
3353     PARAMETER ( CM7 = CM2 + ZWM )
3354     PARAMETER ( CM8ARG = CM7*ZCM*RZWM / (CM2+ZCM) )
3355     PARAMETER ( YCM = 6. / ( 1. + 6.*CM1*ZCM ) )
3356    
3357 molod 1.11 integer INTSTB(irun),INTZ0(irun)
3358 molod 1.14 _RL ZZ0(irun),Z(irun),Z2(irun),Z1(irun),Z0(irun)
3359     _RL X0(irun),X1(irun),Y0(irun),Y1(irun)
3360     _RL PSI2(irun),TEMP(irun)
3361     _RL HZ(irun),ARG0(irun),ARG1(irun),DX(irun)
3362     _RL X0NUM(irun),X1NUM(irun),X0DEN(irun)
3363     _RL X1DEN(irun),Y1DEN(irun),Z2ZWM(irun)
3364     _RL cm3,cm4,cm5,cm8
3365 molod 1.23 integer ibit,indx
3366 molod 1.11 integer i
3367 molod 1.1 C
3368     CM3 = sqrt( 0.2/CM1-0.01 )
3369     CM4 = 1./CM3
3370     CM5 = (10.-CM1) / (10.*CM1*CM3)
3371     CM8 = 6. * LOG(CM8ARG)
3372     C
3373     DO 9000 I = 1,IRUN
3374     VPSIM(I) = 0.
3375     VPSIH(I) = 0.
3376     VX(I) = 0.
3377     VXS(I) = 0.
3378     VY(I) = 0.
3379     VYS(I) = 0.
3380     ZZ0(I) = VZH(I)*VZZ(I)
3381     9000 CONTINUE
3382     IBIT = 0
3383     DO 9122 I = 1,IRUN
3384     IF(VZZ(I).LE.-1.e-7)IBIT = IBIT + 1
3385     9122 CONTINUE
3386     DO 9022 I = 1,IRUN
3387     IF(VZZ(I).LE.-1.e-7)THEN
3388     INTSTB(I) = 1
3389     ELSE
3390     INTSTB(I) = 0
3391     ENDIF
3392     9022 CONTINUE
3393     C
3394     C ****************************************
3395     C ***** UNSTABLE SURFACE LAYER *****
3396     C ****************************************
3397     C
3398     IF(IBIT.LE.0) GO TO 100
3399     C
3400 molod 1.23 indx = 0
3401 molod 1.1 DO 9002 I = 1,IRUN
3402     IF (INTSTB(I).EQ.1)THEN
3403 molod 1.23 indx = indx + 1
3404     Z(indx) = VZZ(I)
3405     Z0(indx) = ZZ0(I)
3406 molod 1.1 ENDIF
3407     9002 CONTINUE
3408     C
3409     DO 9004 I = 1,IBIT
3410     Z(I) = -18. * Z(I)
3411     Z0(I) = -18. * Z0(I)
3412     9004 CONTINUE
3413    
3414     CALL PHI( Z,X1,Y1,IFLAG,IBIT )
3415     CALL PHI( Z0,X0,Y0,IFLAG,IBIT )
3416    
3417     C ****************************
3418     C ***** COMPUTE PSIM *****
3419     C ****************************
3420     C
3421     IF(IFLAG.GE.3) GO TO 75
3422     C
3423     DO 9006 I = 1,IBIT
3424     ARG1(I) = 1. - X1(I)
3425     IF ( Z(I) .LT. 0.013 ) ARG1(I) =
3426     1 Z(I) * ( 0.25 - 0.09375 * Z(I) )
3427     C
3428     ARG0(I) = 1. - X0(I)
3429     IF ( Z0(I) .LT. 0.013 ) ARG0(I) =
3430     1 Z0(I) * ( 0.25 - 0.09375 * Z0(I) )
3431     C
3432     ARG1(I) = ARG1(I) * ( 1.+X0(I) )
3433     ARG0(I) = ARG0(I) * ( 1.+X1(I) )
3434     DX(I) = X1(I) - X0(I)
3435     ARG1(I) = ARG1(I) / ARG0(I)
3436     ARG0(I) = -DX(I) / ( 1. + X1(I)*X0(I) )
3437     ARG0(I) = ATAN( ARG0(I) )
3438     ARG1(I) = LOG( ARG1(I) )
3439     PSI2(I) = 2. * ARG0(I) + ARG1(I)
3440     PSI2(I) = PSI2(I) + DX(I)
3441     9006 CONTINUE
3442     C
3443 molod 1.23 indx = 0
3444 molod 1.1 DO 9008 I = 1,IRUN
3445     IF( INTSTB(I).EQ.1 ) THEN
3446 molod 1.23 indx = indx + 1
3447     VPSIM(I) = PSI2(indx)
3448     VX(I) = X1(indx)
3449     VXS(I) = X0(indx)
3450 molod 1.1 ENDIF
3451     9008 CONTINUE
3452     C
3453     C ****************************
3454     C ***** COMPUTE PSIH *****
3455     C ****************************
3456     C
3457     IF(IFLAG.EQ.2) GO TO 100
3458     C
3459     75 CONTINUE
3460     DO 9010 I = 1,IBIT
3461     ARG1(I) = 1. - Y1(I)
3462     IF( Z(I) .LT. 0.0065 ) ARG1(I) =
3463     1 Z(I) * ( 0.5 - 0.625 * Z(I) )
3464     C
3465     ARG0(I) = 1. - Y0(I)
3466     IF( Z0(I) .LT. 0.0065 ) ARG0(I) =
3467     1 Z0(I) * ( 0.5 - 0.625 * Z0(I) )
3468     C
3469     ARG1(I) = ARG1(I) * ( 1. + Y0(I) )
3470     ARG0(I) = ARG0(I) * ( 1. + Y1(I) )
3471     ARG1(I) = ARG1(I) / ARG0(I)
3472     PSI2(I) = LOG( ARG1(I) )
3473     PSI2(I) = PSI2(I) - Y1(I) + Y0(I)
3474     9010 CONTINUE
3475     C
3476 molod 1.23 indx = 0
3477 molod 1.1 DO 9012 I = 1,IRUN
3478     IF( INTSTB(I).EQ.1 ) THEN
3479 molod 1.23 indx = indx + 1
3480     VPSIH(I) = PSI2(indx)
3481     VY(I) = Y1(indx)
3482     VYS(I) = Y0(indx)
3483 molod 1.1 ENDIF
3484     9012 CONTINUE
3485     C
3486     C **************************************
3487     C ***** STABLE SURFACE LAYER *****
3488     C **************************************
3489     C
3490     100 CONTINUE
3491     IBIT = 0
3492     DO 9114 I = 1,IRUN
3493     IF(VZZ(I).GT.-1.e-7)THEN
3494     IBIT = IBIT + 1
3495     ENDIF
3496     9114 CONTINUE
3497     DO 9014 I = 1,IRUN
3498     IF(VZZ(I).GT.-1.e-7)THEN
3499     INTSTB(I) = 1
3500     ELSE
3501     INTSTB(I) = 0
3502     ENDIF
3503     9014 CONTINUE
3504     IF(IBIT.LE.0) GO TO 300
3505 molod 1.23 indx = 0
3506 molod 1.12 #ifdef CRAY
3507 molod 1.1 CDIR$ NOVECTOR
3508     #endif
3509     DO 9016 I = 1,IRUN
3510     IF (INTSTB(I).EQ.1)THEN
3511 molod 1.23 indx = indx + 1
3512     Z(indx) = VZZ(I)
3513     Z0(indx) = ZZ0(I)
3514     ARG1(indx) = VZH(I)
3515 molod 1.1 ENDIF
3516     9016 CONTINUE
3517 molod 1.12 #ifdef CRAY
3518 molod 1.1 CDIR$ VECTOR
3519     #endif
3520    
3521     DO 9018 I = 1,IBIT
3522     HZ(I) = 1. / ARG1(I)
3523     Z1(I) = Z(I)
3524     Z2(I) = ZWM
3525     C
3526     IF ( Z(I) .GT. ZWM ) THEN
3527     Z1(I) = ZWM
3528     Z2(I) = Z(I)
3529     ENDIF
3530     C
3531     IF ( Z0(I) .GT. Z0M ) THEN
3532     Z0(I) = Z0M
3533     INTZ0(I) = 1
3534     ELSE
3535     INTZ0(I) = 0
3536     ENDIF
3537     C
3538     X1NUM(I) = 1. + 5. * Z1(I)
3539     X0NUM(I) = 1. + 5. * Z0(I)
3540     X1DEN(I) = 1. / (1. + CM1 * (X1NUM(I) * Z1(I)) )
3541     X0DEN(I) = 1. + CM1 * (X0NUM(I) * Z0(I))
3542     C
3543     IF ( (INTZ0(I).EQ.1) .OR. (Z(I).GT.ZWM) )
3544     1 HZ(I) = Z1(I) / Z0(I)
3545     ARG1(I) = HZ(I)*HZ(I)*X0DEN(I)*X1DEN(I)
3546     ARG1(I) = LOG( ARG1(I) )
3547     ARG1(I) = 0.5 * ARG1(I)
3548     ARG0(I) = (Z1(I) + 0.1) * (Z0(I) + 0.1)
3549     ARG0(I) = CM3 + ARG0(I) * CM4
3550     ARG0(I) = ( Z1(I) - Z0(I) ) / ARG0(I)
3551     ARG0(I) = ATAN( ARG0(I) )
3552     TEMP(I) = ARG1(I) + CM5 * ARG0(I)
3553     C
3554     X0(I) = X0NUM(I) / X0DEN(I)
3555     IF ( INTZ0(I).EQ.1 ) X0(I) = 0.
3556     Z2ZWM(I) = Z2(I) * RZWM
3557     9018 CONTINUE
3558     C
3559     C ****************************
3560     C ***** COMPUTE PSIM *****
3561     C ****************************
3562     C
3563     IF( IFLAG.GE.3 ) GO TO 225
3564     C
3565     DO 9020 I = 1,IBIT
3566     X1(I) = X1NUM(I) * X1DEN(I)
3567     ARG1(I) = LOG( Z2ZWM(I) )
3568     PSI2(I) = TEMP(I) + CM6 * ARG1(I)
3569     9020 CONTINUE
3570     C
3571 molod 1.23 indx = 0
3572 molod 1.1 DO 9030 I = 1,IRUN
3573     IF( INTSTB(I).EQ.1 ) THEN
3574 molod 1.23 indx = indx + 1
3575     VPSIM(I) = PSI2(indx)
3576     VX(I) = X1(indx)
3577     VXS(I) = X0(indx)
3578 molod 1.1 ENDIF
3579     9030 CONTINUE
3580     C
3581     C ****************************
3582     C ***** COMPUTE PSIH *****
3583     C ****************************
3584     C
3585     IF(IFLAG.EQ.2)GO TO 300
3586     C
3587     225 CONTINUE
3588     DO 9024 I = 1,IBIT
3589     Y1DEN(I) = 1. + CM1 * ( X1NUM(I) * Z(I) )
3590     Y1(I) = X1NUM(I) / Y1DEN(I)
3591     ARG1(I) = CM7 * Z2ZWM(I) / ( CM2 + Z2(I) )
3592     ARG0(I) = 6.
3593     IF ( Z2(I) .GT. ZCM ) THEN
3594     Y1(I) = YCM
3595     ARG1(I) = Z2(I) * RZCM
3596     ARG0(I) = YCM
3597     TEMP(I) = TEMP(I) + CM8
3598     ENDIF
3599     ARG1(I) = LOG( ARG1(I) )
3600     PSI2(I) = TEMP(I) + ARG0(I) * ARG1(I)
3601     9024 CONTINUE
3602     C
3603 molod 1.23 indx = 0
3604 molod 1.1 DO 9026 I = 1,IRUN
3605     IF( INTSTB(I).EQ.1 ) THEN
3606 molod 1.23 indx = indx + 1
3607     VPSIH(I) = PSI2(indx)
3608     VY(I) = Y1(indx)
3609     VYS(I) = X0(indx)
3610 molod 1.1 ENDIF
3611     9026 CONTINUE
3612     C
3613     300 CONTINUE
3614     C
3615     RETURN
3616     END
3617     SUBROUTINE TRBLEN (STRT,DW2,DZ3,Q,VKZE,VKZM,DTHV,DPK,DU,DV,XL,
3618     1 QXLM,NLEV,INIT,LMIN,LMINQ,LMINQ1,CP,INT1,INT2,
3619     2 DZITRP,STBFCN,XL0,Q1,WRKIT1,WRKIT2,WRKIT3,WRKIT4,irun)
3620     C**********************************************************************
3621     C
3622     C SUBROUTINE TRBLEN - COMPUTES TURBULENT LENGTH SCALE
3623     C - CALLED FROM TRBFLX
3624     C ARGUMENTS ::
3625     C
3626     C INPUT:
3627     C ------
3628     C STRT - BRUNT VAISALA FREQUENCY
3629     C DW2 - SQUARED SHEAR
3630     C DZ3 - LAYER THICKNESS FOR LENGTH SCALE CALC.
3631     C Q - TURBULENCE VELOCITY
3632     C VKZE - VK * Z AT LAYER EDGES
3633     C VKZM - VK * Z AT LAYER CENTERS
3634     C DTHV - VERTICAL GRADIENT OF THV
3635     C DPK - VERTICAL GRADIENT OF PK
3636     C DU - VERTICAL GRADIENT OF U
3637     C DV - VERTICAL GRADIENT OF V
3638     C NLEV - NUMBER OF ATMOSPHERIC LEVELS
3639     C INIT - INPUT FLAG : 1 = INITIAL START
3640     C 2 = 2ND CALL FOR INITIAL STAR
3641     C 0 = NON-INITIAL START
3642     C
3643     C OUTPUT:
3644     C -------
3645     C XL - TURBULENT LENGTH SCALE
3646     C QXLM - TURBULENT LENGTH SCALE * Q AT LAYER CENTER
3647     C LMIN - HIGHEST LAYER WHERE INSTABILITY OCCURS
3648     C LMINQ - HIGHEST LAYER WHERE TURBULENCE OCCURS
3649     C
3650     C SUBPROGRAMS NEEDED ::
3651     C
3652     C TRBITP - INTERPOLATES TO HEIGHT WHERE RI = RICR
3653     C
3654     C**********************************************************************
3655 molod 1.11 implicit none
3656    
3657     C Argument List Declarations
3658     integer irun,nlev,init,lmin,lminq,lminq1
3659 molod 1.14 _RL cp
3660     _RL STRT(irun,NLEV),DW2(irun,NLEV),DZ3(irun,NLEV)
3661     _RL Q(irun,NLEV),VKZM(irun,NLEV-1),VKZE(irun,NLEV-1)
3662     _RL DTHV(irun,NLEV),DPK(irun,NLEV),DU(irun,NLEV)
3663     _RL DV(irun,NLEV)
3664     _RL QXLM(irun,NLEV-1),XL(irun,NLEV-1)
3665     _RL DZITRP(irun,nlev-1),STBFCN(irun,nlev)
3666     _RL XL0(irun,nlev),Q1(irun,nlev-1)
3667     _RL WRKIT1(irun,nlev-1),WRKIT2(irun,nlev-1)
3668     _RL WRKIT3(irun,nlev-1)
3669     _RL WRKIT4(irun,nlev-1)
3670 molod 1.11 INTEGER INT1(irun,nlev), INT2(irun,nlev-1)
3671    
3672     C Local Variables
3673 molod 1.14 _RL rf1,rf2,e5,d4,d1,rfc,ricr,alpha,dzcnv,xl0cnv,xl0min
3674     _RL clmt,clmt53
3675 molod 1.1 PARAMETER ( RF1 = 0.2340678 )
3676     PARAMETER ( RF2 = 0.2231172 )
3677     PARAMETER ( E5 = 49.66 )
3678     PARAMETER ( D4 = 2.6532122E-2 )
3679     PARAMETER ( D1 = D4 * E5 )
3680     PARAMETER ( RFC = 0.1912323 )
3681     PARAMETER ( RICR = ( (RF1-RFC)*RFC ) / ( (RF2-RFC)*D1 ) )
3682     PARAMETER ( ALPHA = 0.1 )
3683     PARAMETER ( DZCNV = 100. )
3684     PARAMETER ( XL0CNV = DZCNV * ALPHA )
3685     PARAMETER ( XL0MIN = 1. )
3686     PARAMETER ( CLMT = 0.23 )
3687     PARAMETER ( CLMT53 = 5. * CLMT / 3. )
3688 molod 1.11
3689     integer ibit,nlevm1,nlevp1,istnlv,istnm1,nlevml,istnml,Lp1
3690     integer istnmq,istlmq,lminp,lm1,lmin1
3691     integer i,L,LL
3692 molod 1.1 C
3693     NLEVM1 = NLEV - 1
3694     NLEVP1 = NLEV + 1
3695     ISTNLV = irun * NLEV
3696     ISTNM1 = irun * NLEVM1
3697     C
3698     IF ( INIT.EQ.2 ) GO TO 1200
3699     C
3700     C COMPUTE DEPTHS OF UNSTABLE LAYERS
3701     C =================================
3702     DO 10 I=1,ISTNLV
3703     STBFCN(I,1) = STRT(I,1) - RICR * DW2(I,1)
3704     INT1(I,1) = 0
3705     IF( STBFCN(I,1).LE.0. ) INT1(I,1) = 1
3706     10 CONTINUE
3707     DO 20 I=1,ISTNM1
3708     INT2(I,1) = 0
3709     IF( INT1(I,1).EQ.1 .XOR. INT1(I,2).EQ.1 ) INT2(I,1) = 1
3710     20 CONTINUE
3711     C
3712     DO 40 LMIN = 1,NLEV
3713     IBIT = 0
3714     DO 30 I=1,irun
3715     IBIT = IBIT + INT1(I,LMIN)
3716     30 CONTINUE
3717     IF(IBIT.GE.1) GO TO 50
3718     40 CONTINUE
3719     LMIN = NLEVP1
3720     50 CONTINUE
3721     LMIN = 1
3722     C
3723     DO 60 I=1,ISTNM1
3724     XL0(I,1) = 0.
3725     60 CONTINUE
3726     DO 70 I=1,irun
3727     XL0(I,NLEV) = DZ3(I,NLEV)
3728     70 CONTINUE
3729     C
3730     IF(LMIN.GE.NLEVP1) GOTO 1100
3731     LMIN1 = LMIN - 1
3732     IF(LMIN1.EQ.0) LMIN1 = 1
3733     NLEVML = NLEV - LMIN1
3734     ISTNML = irun*NLEVML
3735     CALL TRBITP ( STBFCN(1,LMIN1),INT2(1,LMIN1),DTHV(1,LMIN1),
3736     . DPK(1,LMIN1), DU(1,LMIN1), DV(1,LMIN1),
3737     . DZITRP(1,LMIN1), NLEVML,
3738     . WRKIT1,WRKIT2,WRKIT3,WRKIT4,CP,irun )
3739     LP1 = LMIN1 + 1
3740     C
3741     DO 80 I=1,ISTNML
3742     INT2(I,LMIN1) = 0
3743     IF( INT1(I,LMIN1).EQ.1 .OR. INT1(I,LP1).EQ.1 ) INT2(I,LMIN1) = 1
3744     IF( INT2(I,LMIN1).EQ.1 )
3745     . XL0(I,LMIN1) = (0.5+DZITRP(I,LMIN1)) * DZ3(I,LP1)
3746     80 CONTINUE
3747     DO 90 I=1,irun
3748     INT2(I,NLEVM1) = INT1(I,NLEV)
3749     90 CONTINUE
3750     C
3751     DO 100 I=1,ISTNML
3752     IF( INT2(I,LMIN1).EQ.1 ) THEN
3753     XL0(I,LP1) = XL0(I,LP1) + ( (0.5-DZITRP(I,LMIN1)) * DZ3(I,LP1) )
3754     ENDIF
3755     100 CONTINUE
3756     IF (LMIN.GT.1) GOTO 400
3757     DO 110 I=1,irun
3758     IF( INT1(I,1).EQ.1 ) XL0(I,1) = XL0(I,1) + DZ3(I,1)
3759     110 CONTINUE
3760     400 CONTINUE
3761     C
3762     LMINP = LMIN + 1
3763     IF(LMINP.GT.NLEVM1) GOTO 550
3764     DO 500 L = LMINP,NLEVM1
3765     LM1 = L-1
3766     DO 120 I = 1,irun
3767     IF( INT1(I,LM1).EQ.1 ) XL0(I,L) = XL0(I,L) + XL0(I,LM1)
3768     120 CONTINUE
3769     500 CONTINUE
3770     550 CONTINUE
3771     IF(LMIN.GT.NLEVM1) GOTO 600
3772     DO 130 I = 1,irun
3773     IF( INT1(I,NLEVM1).EQ.1 .AND. INT1(I,NLEV).EQ.1 ) THEN
3774     XL0(I,NLEV) = XL0(I,NLEV) + XL0(I,NLEVM1)
3775     ENDIF
3776     130 CONTINUE
3777     IF(LMIN.GT.NLEV) GOTO 1100
3778     600 CONTINUE
3779     DO 1000 LL = LMIN,NLEV-1
3780     L = NLEVM1 + LMIN - LL
3781     LP1 = L+1
3782     DO 140 I = 1,irun
3783     IF( INT1(I,LP1).EQ.1 ) THEN
3784     IF( INT1(I,L) .EQ.1 ) THEN
3785     XL0(I,L) = XL0(I,LP1)
3786     ELSE
3787     XL0(I,L) = XL0(I,L) + XL0(I,LP1)
3788     ENDIF
3789     ENDIF
3790     140 CONTINUE
3791     1000 CONTINUE
3792     1100 CONTINUE
3793     C
3794     DO 150 I = 1,ISTNLV
3795     IF( XL0(I,1).LT.XL0CNV ) XL0(I,1) = XL0CNV
3796     150 CONTINUE
3797     C
3798     C *********************************************************************
3799     C **** DETERMINE MIXING LENGTHS FOR STABLE LAYERS ***
3800     C *********************************************************************
3801     C
3802     IF(INIT.EQ.1) GOTO 1400
3803     C
3804     IF(LMINQ.GT.1) THEN
3805     ISTLMQ = irun * LMINQ1
3806     DO 160 I = 1,ISTLMQ
3807     INT2(I,1) = 1 - INT1(I,1)
3808     160 CONTINUE
3809     ENDIF
3810     IF(LMINQ.LT.NLEV) THEN
3811     ISTNMQ = irun * (NLEV-LMINQ)
3812     DO 170 I = 1,ISTNMQ
3813     IF( INT1(I,LMINQ).EQ.0 ) THEN
3814     XL0(I,LMINQ) = Q(I,LMINQ) / XL0(I,LMINQ)
3815     XL0(I,LMINQ) = XL0(I,LMINQ) * XL0(I,LMINQ) + 1.0E-20
3816     XL0(I,LMINQ) = STBFCN(I,LMINQ) + XL0(I,LMINQ)
3817     XL0(I,LMINQ) = SQRT( XL0(I,LMINQ) )
3818     XL0(I,LMINQ) = Q(I,LMINQ) / XL0(I,LMINQ)
3819     ENDIF
3820     INT2(I,LMINQ) = 0
3821     IF( XL0(I,LMINQ).LT.XL0MIN ) INT2(I,LMINQ) = 1
3822     170 CONTINUE
3823     ENDIF
3824     C
3825     1200 CONTINUE
3826     C
3827     IF(INIT.EQ.2) THEN
3828     DO 180 I = 1,ISTNM1
3829     INT2(I,1) = 1 - INT1(I,1)
3830     180 CONTINUE
3831     ENDIF
3832     DO 190 I = 1,ISTNM1
3833     IF( INT2(I,1).EQ.1 ) XL0(I,1) = XL0MIN
3834     190 CONTINUE
3835     C
3836     C *********************************************************************
3837     C **** LENGTH SCALE XL FROM XL0 AND VKZE ****
3838     C *********************************************************************
3839     C
3840     1400 CONTINUE
3841     C
3842     DO 200 I = 1,ISTNM1
3843     XL(I,1) = XL0(I,1) * VKZE(I,1) / ( XL0(I,1)+VKZE(I,1) )
3844     200 CONTINUE
3845     C
3846     C *********************************************************************
3847     C **** CLMT53 TIMES Q TIMES LENGTH SCALE AT MID LEVELS ***
3848     C *********************************************************************
3849     C
3850     IF(INIT.EQ.1) GOTO 1700
3851     ISTNMQ = irun * (NLEV-LMINQ1)
3852     DO 210 I = 1,ISTNMQ
3853     Q1(I,LMINQ1) = Q(I,LMINQ1)
3854     INT1(I,LMINQ1) = 0
3855     IF( Q(I,LMINQ1).LE.Q(I,LMINQ1+1) ) INT1(I,LMINQ1) = 1
3856     IF( INT1(I,LMINQ1).EQ.1 ) THEN
3857     XL0(I,LMINQ1) = XL0(I,LMINQ1+1)
3858     Q1(I,LMINQ1) = Q(I,LMINQ1+1)
3859     ENDIF
3860     210 CONTINUE
3861     C
3862     DO 240 I = 1,ISTNMQ
3863     QXLM(I,LMINQ1) = XL0(I,LMINQ1)*VKZM(I,LMINQ1)
3864     . / ( XL0(I,LMINQ1)+VKZM(I,LMINQ1) )
3865     QXLM(I,LMINQ1) = CLMT53 * Q1(I,LMINQ1)*QXLM(I,LMINQ1)
3866     240 CONTINUE
3867     C
3868     1700 CONTINUE
3869     C
3870     RETURN
3871     END
3872     SUBROUTINE TRBITP ( STBFCN,INTCHG,DTHV,DPK,DU,DV,DZITRP,NLEV,
3873     . AAA,BBB,CCC,DDD,CP,irun )
3874     C**********************************************************************
3875     C
3876     C SUBROUTINE TRBITP - INTERPOLATES TO THE HEIGHT WHERE RI EQUALS RICR
3877     C - CALLED FROM TRBLEN
3878     C ARGUMENTS ::
3879     C
3880     C INPUT:
3881     C ------
3882     C STBFCN - DTHV * DPK - RICR*( DU*DU + DV*DV)
3883     C INTCHG - INT '1' AT LEVELS WHERE STBFCN CHANGES SIG
3884     C DTHV - VERTICAL GRADIENT OF THV
3885     C DPK - VERTICAL GRADIENT OF PK
3886     C DU - VERTICAL GRADIENT OF U
3887     C DV - VERTICAL GRADIENT OF V
3888     C NLEV - NUMBER OF LEVELS TO BE PROCESSED
3889     C
3890     C OUTPUT:
3891     C -------
3892     C DZITRP - INTERPOLATION COEFFICIENT
3893     C
3894     C**********************************************************************
3895 molod 1.11 implicit none
3896    
3897     C Argument List Declarations
3898     integer irun,nlev
3899 molod 1.14 _RL cp
3900     _RL STBFCN(irun,NLEV+1)
3901 molod 1.11 integer INTCHG(irun,NLEV)
3902 molod 1.14 _RL DTHV(irun,NLEV+1),DPK(irun,NLEV+1)
3903     _RL DU(irun,NLEV+1),DV(irun,NLEV+1)
3904     _RL DZITRP(irun,NLEV+1)
3905     _RL AAA(irun,NLEV),BBB(irun,NLEV)
3906     _RL CCC(irun,NLEV),DDD(irun,NLEV)
3907 molod 1.11
3908     C Local Variables
3909 molod 1.14 _RL rf1,rf2,e5,d4,d1,rfc,ricr
3910 molod 1.1 PARAMETER ( RF1 = 0.2340678 )
3911     PARAMETER ( RF2 = 0.2231172 )
3912     PARAMETER ( E5 = 49.66 )
3913     PARAMETER ( D4 = 2.6532122E-2 )
3914     PARAMETER ( D1 = D4 * E5 )
3915     PARAMETER ( RFC = 0.1912323 )
3916     PARAMETER ( RICR = ( (RF1-RFC)*RFC ) / ( (RF2-RFC)*D1 ) )
3917    
3918 molod 1.11 integer istnlv
3919     integer i
3920 molod 1.1 C
3921     C *********************************************************************
3922     C **** QUADRATIC INTERPOLATION OF RI TO RICR VIA ***
3923     C **** LINEAR INTERPOLATION OF DTHV, DPK, DU & DV ***
3924     C *********************************************************************
3925     C
3926     ISTNLV = irun*NLEV
3927     DO 10 I=1,ISTNLV
3928     DZITRP(I,1) = 0.
3929     10 CONTINUE
3930     DO 20 I=1,ISTNLV
3931     IF( INTCHG(I,1).EQ.1 ) THEN
3932     DDD(I,1) = ( CP *(DTHV(I,2)*DPK(I,1)
3933     . + DTHV(I,1)*DPK(I,2)) )
3934     . - ( (2.*RICR) * ( DU(I,2)* DU(I,1)
3935     . + DV(I,2)* DV(I,1)) )
3936     AAA(I,1) = STBFCN(I,1) + STBFCN(I,2)
3937     BBB(I,1) = STBFCN(I,1) - STBFCN(I,2)
3938     CCC(I,1) = 1. / BBB(I,1)
3939     DZITRP(I,1) = AAA(I,1) * CCC(I,1)
3940     AAA(I,1) = AAA(I,1) - DDD(I,1)
3941     DDD(I,1) = ( DDD(I,1) * DDD(I,1) )
3942     . - 4. * (STBFCN(I,2) * STBFCN(I,1) )
3943     DDD(I,1) = DDD(I,1)*CCC(I,1)*CCC(I,1)
3944     DDD(I,1) = SQRT( DDD(I,1) )
3945     ENDIF
3946     C
3947     IF( INTCHG(I,1).EQ.1 .AND. AAA(I,1).NE.0. ) THEN
3948     DZITRP(I,1) = ( BBB(I,1)*(1.-DDD(I,1)) ) / AAA(I,1)
3949     ENDIF
3950     C
3951     DZITRP(I,1) = 0.5 * DZITRP(I,1)
3952     20 CONTINUE
3953     C
3954     RETURN
3955     END
3956     SUBROUTINE TRBL20 (RI,STRT,DW2,XL,ZKM,ZKH,QE,QQE,INTSTB,NLEV,
3957     1 nlay,irun)
3958     C**********************************************************************
3959     C
3960     C SUBROUTINE TRBL20 - COMPUTES QE AND DIMLESS COEFS FROM
3961     C MELLOR-YAMADA LEVEL 2 MODEL
3962     C - CALLED FROM AND FROM TRBFLX
3963     C ARGUMENTS ::
3964     C
3965     C INPUT:
3966     C ------
3967     C RI - RICHARDSON NUMBER
3968     C STRT - BRUNT VAISALA FREQUENCY
3969     C DW2 - SQUARED SHEAR
3970     C XL - TURBULENT LENGTH SCALE
3971     C NLEV - NUMBER OF LEVELS TO BE PROCESSED
3972     C
3973     C OUTPUT:
3974     C -------
3975     C ZKM - MOMENTUM TRANSPORT COEFFICIENT
3976     C ZKH - HEAT TRANSPORT COEFFICIENT
3977     C QE - EQUILIBRIUM TURBULENT VELOCITY SCALE
3978     C QQE - EQUILIBRIUM TURBULENT KINETIC ENERGY
3979     C BITSTB - BIT '1' WHERE QE GREATER THAN ZERO
3980     C
3981     C**********************************************************************
3982 molod 1.11 implicit none
3983    
3984     C Argument List Declarations
3985     integer nlev,nlay,irun
3986 molod 1.14 _RL RI(irun,NLEV),STRT(irun,NLEV),DW2(irun,NLEV)
3987     _RL XL(irun,NLEV),ZKM(irun,NLEV),ZKH(irun,NLEV)
3988     _RL QE(irun,NLEV),QQE(irun,NLEV)
3989 molod 1.11 INTEGER INTSTB(irun,nlev)
3990 molod 1.14 _RL EE(irun,nlay-1),RF(irun,nlay-1)
3991 molod 1.11
3992     C Local Variables
3993 molod 1.14 _RL b1,b2,d3,rf1,rf2,d3b2,d2,e5,d4,d1,d1half,d2half
3994     _RL rfc,ricr,ch,cm
3995 molod 1.1 PARAMETER ( B1 = 16.6 )
3996     PARAMETER ( B2 = 10.1 )
3997     PARAMETER ( D3 = 0.29397643 )
3998     PARAMETER ( RF1 = 0.2340678 )
3999     PARAMETER ( RF2 = 0.2231172 )
4000     PARAMETER ( D3B2 = D3 / RF1 )
4001     PARAMETER ( D2 = RF1 )
4002     PARAMETER ( E5 = 49.66 )
4003     PARAMETER ( D4 = 2.6532122E-2 )
4004     PARAMETER ( D1 = D4 * E5 )
4005     PARAMETER ( D1HALF = 0.5 * D1 )
4006     PARAMETER ( D2HALF = 0.5 * D2 )
4007     PARAMETER ( RFC = 0.1912323 )
4008     PARAMETER ( RICR = ( (RF1-RFC)*RFC ) / ( (RF2-RFC)*D1 ) )
4009     PARAMETER ( CH = 2.5828674 )
4010     PARAMETER ( CM = CH / D1 )
4011    
4012 molod 1.11 integer istnlv
4013     integer i
4014    
4015 molod 1.1 ISTNLV = irun * NLEV
4016     C
4017     C *********************************************************************
4018     C **** COMPUTE FLUX RICHARDSON NUMBER ***
4019     C *********************************************************************
4020     C
4021     DO 10 I=1,ISTNLV
4022     EE(I,1) = D1HALF * RI(I,1) + D2HALF
4023     RF(I,1) = EE(I,1)* EE(I,1)
4024     RF(I,1) = RF(I,1)- D3*RI(I,1)
4025     RF(I,1) = SQRT( RF(I,1) )
4026     RF(I,1) = EE(I,1) - RF(I,1)
4027     C
4028     IF( RI(I,1).LE.1.e-4 .AND. RI(I,1).GE.-1.e-4 ) THEN
4029     RF(I,1) = D3B2*RI(I,1)
4030     ENDIF
4031     C
4032     C *********************************************************************
4033     C **** QE AND DIMENSIONLESS DIFFUSION COEFICIENTS ***
4034     C **** FROM LEVEL 2 CLOSURE MODEL ***
4035     C *********************************************************************
4036     C
4037     IF( RI(I,1).LT.RICR .AND. RF(I,1).LT.RFC ) THEN
4038     ZKH(I,1) = ( RFC-RF(I,1) ) / (1.-RF(I,1))
4039     ZKM(I,1) = CM * (RF1-RF(I,1))
4040     ZKM(I,1) = ZKH(I,1)*ZKM(I,1) / (RF2-RF(I,1))
4041     ZKH(I,1) = CH *ZKH(I,1)
4042     QE(I,1) = ZKM(I,1)*DW2(I,1) - ZKH(I,1)*STRT(I,1)
4043     ENDIF
4044     C
4045     IF( QE(I,1).LT.1.e-14 ) THEN
4046     INTSTB(I,1) = 0
4047     QE(I,1) = 0.
4048     ELSE
4049     INTSTB(I,1) = 1
4050     QE(I,1) = B1*QE(I,1)
4051     QE(I,1) = SQRT( QE(I,1) )
4052     QE(I,1) = XL(I,1)*QE(I,1)
4053     ENDIF
4054     QQE(I,1) = 0.5 * QE(I,1) * QE(I,1)
4055     10 CONTINUE
4056     C
4057     RETURN
4058     END
4059     SUBROUTINE TRBL25(Q,XL,STRT,DW2,INTSTB,INTQ,ZKM,ZKH,P3,NLEV,
4060     1 nlay,irun)
4061     C**********************************************************************
4062     C
4063     C SUBROUTINE TRBL25 - COMPUTES P3 AND DIMLESS COEFS FROM
4064     C MELLOR-YAMADA LEVEL 2.5 MODEL
4065     C - CALLED FROM TRBFLX
4066     C
4067     C ARGUMENTS ::
4068     C
4069     C INPUT:
4070     C ------
4071     C Q - TURBULENCE VELOCITY
4072     C XL - TURBULENT LENGTH SCALE
4073     C STRT - BRUNT VAISALA FREQUENCY
4074     C DW2 - SQUARED SHEAR
4075     C BITSTB - BIT '1' WHERE QE GREATER THAN ZERO
4076     C NLEV - NUMBER OF LEVELS TO BE PROCESSED
4077     C
4078     C OUTPUT:
4079     C -------
4080     C ZKM - MOMENTUM TRANSPORT COEFFICIENT
4081     C ZKH - HEAT TRANSPORT COEFFICIENT
4082     C P3 - PRODUCTION RATE OF TURBULENT KINETIC ENERG
4083     C
4084     C**********************************************************************
4085 molod 1.11 implicit none
4086    
4087     C Argument list Declarations
4088     integer nlev,nlay,irun
4089 molod 1.14 _RL Q(irun,NLEV),XL(irun,NLEV),STRT(irun,NLEV)
4090     _RL DW2(irun,NLEV)
4091 molod 1.11 INTEGER INTSTB(irun,nlay), INTQ(irun,nlay)
4092 molod 1.14 _RL ZKM(irun,NLEV),ZKH(irun,NLEV),P3(irun,NLEV)
4093 molod 1.11
4094     C Local Variables
4095 molod 1.14 _RL a1,a2,a4,c1,a5,a3,b1,b2,b3,ff2,ff3,ff4
4096 molod 1.1 PARAMETER ( A1 = 0.92 )
4097     PARAMETER ( A2 = 0.74 )
4098     PARAMETER ( A4 = 6. * A1 * A1)
4099     PARAMETER ( C1 = 0.08 )
4100     PARAMETER ( A5 = 3.*C1*(-1.) )
4101     PARAMETER ( A3 = A4 * A5*(-1.) )
4102     PARAMETER ( B1 = 16.6 )
4103     PARAMETER ( B2 = 10.1 )
4104     PARAMETER ( B3 = 1. / B1 )
4105     PARAMETER ( FF2 = 9. * A1 * A2 )
4106     PARAMETER ( FF3 = (3.*A2*B2) - (9.*A2*A2 ) )
4107     PARAMETER ( FF4 = (3.*A2*B2) + (12.*A1*A2 ) )
4108    
4109 molod 1.14 _RL F2(irun,nlay-1),F3(irun,nlay-1)
4110     _RL F4(irun,nlay-1),XQ(irun,nlay-1)
4111 molod 1.11
4112     integer istnlv
4113     integer i
4114 molod 1.1 C
4115     ISTNLV = irun * NLEV
4116     C
4117     C *********************************************************************
4118     C **** P3 AND DIMENSIONLESS DIFFUSION COEFICIENTS ***
4119     C **** FROM LEVEL 2.5 CLOSURE MODEL ***
4120     C *********************************************************************
4121     C
4122     DO 10 I=1,ISTNLV
4123     IF( INTQ(I,1).EQ.1 .AND. INTSTB(I,1).EQ.0 ) THEN
4124     XQ(I,1) = XL(I,1) / Q(I,1)
4125     XQ(I,1) = XQ(I,1) * XQ(I,1)
4126     STRT(I,1) = XQ(I,1) * STRT(I,1)
4127     DW2(I,1) = XQ(I,1) * DW2(I,1)
4128     F2(I,1) = 1.+FF2 * STRT(I,1)
4129     F3(I,1) = 1.+FF3 * STRT(I,1)
4130     F4(I,1) = 1.+FF4 * STRT(I,1)
4131     ZKH(I,1) = (F4(I,1) * F2(I,1))
4132     . + A4 * (F3(I,1) * DW2(I,1))
4133     ZKH(I,1) = (F2(I,1) + A3*DW2(I,1))
4134     . / ZKH(I,1)
4135     ZKM(I,1) = A1 * (F3(I,1)*ZKH(I,1)+A5)
4136     . / F2(I,1)
4137     ZKH(I,1) = A2 * ZKH(I,1)
4138     P3(I,1) = ZKH(I,1)*STRT(I,1) + B3
4139     P3(I,1) = 2. * ( ZKM(I,1)*DW2(I,1) - P3(I,1) )
4140     P3(I,1) = P3(I,1)*Q(I,1)
4141     C
4142     ENDIF
4143     10 CONTINUE
4144     C
4145     RETURN
4146     END
4147     SUBROUTINE TRBDIF ( XX1,XX2,RHOKDZ,FLXFAC,DXX1G,DXX2G,NLEV,
4148     . ITYPE,EPSL,irun )
4149     C
4150     C**********************************************************************
4151     C
4152     C ARGUMENTS ::
4153     C
4154     C INPUT:
4155     C ------
4156     C XX1 - FIRST PROPERTY TO BE DIFFUSED
4157     C (INPUT INCLUDES FORWARD PRODUCTION TERM)
4158     C XX2 - SECOND PROPERTY TO BE DIFFUSED (V-WIND)
4159     C (INPUT INCLUDES FORWARD PRODUCTION TERM)
4160     C -OR-
4161     C CHANGE IN XX1 DUE TO UNIT CHANGE IN THG
4162     C (TH OR SH PROFILES)
4163     C -OR-
4164     C BACKWARD PRODUCTION TERM (QQ)
4165     C RHOKDZ - RHO * K * WEIGHT / DZ AT INTERFACES
4166     C FLXFAC - G * DT / (DP*WEIGHT) AT EDGES
4167     C NLEV - NUMBER OF ATMOSPHERIC LEVELS
4168     C ITYPE - INTEGER FLAG FOR INPUT TYPE
4169     C 1 = QQ: COMPUTE BACKWARD PRODUCTION AND
4170     C USE UNDERFLOW CUTOFF
4171     C 2 = TH OR SH: COMPUTE TENDENCY DUE TO
4172     C SURFACE PERTURBATION
4173     C 3 = U AND V: COMPUTE BOTH FIELDS
4174     C EPSL - UNDERFLOW CUTOFF CRITERION (QQ ONLY)
4175     C
4176     C OUTPUT:
4177     C ------
4178     C XX1 - NEW VALUE RETURNED
4179     C XX2 - NEW VALUE RETURNED
4180     C DXX1G - SOURCE TERM FOR XX1 AT GROUND
4181     C DXX1G - SOURCE TERM FOR XX2 AT GROUND
4182     C
4183     C**********************************************************************
4184 molod 1.11 implicit none
4185    
4186     C Argument List Declarations
4187     integer nlev,itype,irun
4188 molod 1.14 _RL XX1(irun,NLEV+1),XX2(irun,NLEV+1)
4189     _RL RHOKDZ(irun,NLEV),FLXFAC(irun,NLEV+1)
4190     _RL DXX1G(irun),DXX2G(irun)
4191     _RL epsl
4192 molod 1.1 C
4193 molod 1.14 _RL AA(irun,nlev), BB(irun,nlev), CC(irun,nlev+1)
4194 molod 1.11 integer istnlv,istnm1,nlevp1,istnlx
4195     integer i
4196 molod 1.1 C
4197     ISTNLV = irun * NLEV
4198     ISTNM1 = ISTNLV - irun
4199     NLEVP1 = NLEV + 1
4200     ISTNLX = ISTNM1
4201     IF(ITYPE.EQ.2) ISTNLX = ISTNLV
4202     C
4203     C DEFINE MATRIX
4204     C
4205     DO 10 I=1,irun
4206     CC(I,1) = 0.
4207     10 CONTINUE
4208     DO 20 I=1,ISTNLX
4209     CC(I,2) = RHOKDZ(I,1) * FLXFAC(I,2)
4210     20 CONTINUE
4211     DO 30 I=1,ISTNLV
4212     BB(I,1) = RHOKDZ(I,1) * FLXFAC(I,1)
4213     AA(I,1) = 1. + CC(I,1) + BB(I,1)
4214     30 CONTINUE
4215     C
4216     C ADD IMPLICIT BACKWARD FORCING FOR QQ
4217     IF(ITYPE.EQ.1) THEN
4218     DO 40 I=1,ISTNLV
4219     AA(I,1) = AA(I,1) - XX2(I,1)
4220     40 CONTINUE
4221     ENDIF
4222     C
4223     C SOLVE MATRIX EQUATION FOR XX1
4224     CALL VTRI0(AA,BB,CC,XX1,XX1,NLEV,irun)
4225     C
4226     IF(ITYPE.EQ.2) THEN
4227     C COMPUTE CHANGE AT SURFACE
4228     C
4229     DO 50 I=1,irun
4230     DXX1G(I) = CC(I,NLEVP1) * ( XX1(I,NLEV)-XX1(I,NLEVP1) )
4231     50 CONTINUE
4232     C
4233     C SOLVE MATRIX FOR SURFACE PERTURBATION
4234     CALL VTRI1(AA,BB,XX2,NLEV,irun)
4235     DO 60 I=1,irun
4236     DXX2G(I) = CC(I,NLEVP1) * ( XX2(I,NLEV)-XX2(I,NLEVP1) )
4237     60 CONTINUE
4238     ENDIF
4239     C
4240     C SOLVE MATRIX EQUATION FOR XX2
4241     C
4242     IF(ITYPE.EQ.3) CALL VTRI2 (AA,BB,CC,XX2,XX2,NLEV,irun)
4243     C
4244     C ELIMINATE UNDERFLOW
4245     IF(ITYPE.EQ.1) THEN
4246     DO 70 I=1,ISTNLV
4247     IF( XX1(I,1).LT.EPSL ) XX1(I,1) = 0.
4248     70 CONTINUE
4249     ENDIF
4250     C
4251     RETURN
4252     END
4253     SUBROUTINE VTRI0 ( A,B,C,F,Y,K,irun)
4254 molod 1.11 implicit none
4255    
4256     integer k,irun
4257 molod 1.14 _RL A(irun,K),B(irun,K),C(irun,K),Y(irun,K+1)
4258     _RL F(irun,K)
4259 molod 1.11
4260     integer i,L,Lm1
4261 molod 1.1 C
4262     DO 9000 I = 1,irun
4263     A(I,1) = 1. / A(I,1)
4264     9000 CONTINUE
4265     C
4266     DO 100 L = 2,K
4267     LM1 = L - 1
4268     DO 9002 I = 1,irun
4269     C(I,L) = C(I,L) * A(I,LM1)
4270     A(I,L) = 1. / ( A(I,L) - B(I,LM1) * C(I,L) )
4271     F(I,L) = F(I,L) + F(I,LM1) * C(I,L)
4272     9002 CONTINUE
4273     100 CONTINUE
4274     C
4275     DO 200 L = K,1,-1
4276     DO 9004 I = 1,irun
4277     Y(I,L) = (F(I,L) + B(I,L) * Y(I,L+1)) * A(I,L)
4278     9004 CONTINUE
4279     200 CONTINUE
4280     C
4281     RETURN
4282     END
4283     C
4284     SUBROUTINE VTRI1 ( A,B,Y,K,irun)
4285 molod 1.11 implicit none
4286    
4287     integer k,irun
4288 molod 1.14 _RL A(irun,K),B(irun,K),Y(irun,K+1)
4289 molod 1.11
4290     integer i,L
4291 molod 1.1 C
4292     DO 200 L = K,1,-1
4293     DO 9000 I = 1,irun
4294     Y(I,L) = B(I,L) * Y(I,L+1) * A(I,L)
4295     9000 CONTINUE
4296     200 CONTINUE
4297     C
4298     RETURN
4299     END
4300     C
4301     SUBROUTINE VTRI2 ( A,B,C,F,Y,K,irun)
4302 molod 1.11 implicit none
4303    
4304     integer k,irun
4305 molod 1.14 _RL A(irun,K),B(irun,K),C(irun,K),F(irun,K)
4306     _RL Y(irun,K+1)
4307 molod 1.11
4308     integer i,L
4309 molod 1.1 C
4310     DO 100 L = 2,K
4311     DO 9000 I = 1,irun
4312     F(I,L) = F(I,L) + F(I,L-1) * C(I,L)
4313     9000 CONTINUE
4314     100 CONTINUE
4315     C
4316     DO 200 L = K,1,-1
4317     DO 9002 I = 1,irun
4318     Y(I,L) = (F(I,L) + B(I,L) * Y(I,L+1)) * A(I,L)
4319     9002 CONTINUE
4320     200 CONTINUE
4321     C
4322     RETURN
4323     END
4324     SUBROUTINE LINADJ ( NN,VRIB1,VRIB2,VWS1,VWS2,VZ1,VUSTAR,IWATER,
4325     1 VAPSIM, VAPSIHG,VPSIH,VPSIG,VX,VX0,VY,VY0,ITYPE,LWATER,IRUN,
4326     2 VDZETA,VDZ0,VDPSIM,VDPSIH,INTRIB,
4327     3 VX0PSIM,VG,VG0,VR1MG0,VZ2,VDZSEA,VAZ0,VXNUM1,VPSIGB2,VDX,
4328     4 VDXPSIM,VDY,VXNUM2,VDEN,VAWS1,VXNUM3,VXNUM,VDZETA1,VDZETA2,
4329     5 VZCOEF2,VZCOEF1,VTEMPLIN,VDPSIMC,VDPSIHC)
4330     C
4331     C**********************************************************************
4332     C
4333     C ARGUMENTS ::
4334     C
4335     C INPUT:
4336     C ------
4337     C RIB1 - BULK RICHARDSON NUMBER OF INPUT STATE
4338     C RIB2 - DESIRED BULK RICH NUMBER OF OUTPUT STATE
4339     C WS1 - SURFACE WIND SPEED OF INPUT STATE
4340     C WS2 - DESIRED SURFACE WIND SPEED OF OUTPUT STATE
4341     C Z1 - INPUT VALUE OF ROUGHNESS HEIGHT
4342     C USTAR - INPUT VALUE OF CU * WS
4343     C WATER - BIT ARRAY - '1' WHERE OCEAN
4344     C APSIM - (1/PSIM)
4345     C APSIHG - ( 1 / (PSIH+PSIG) )
4346     C PSIH - NON-DIM TEMP GRADIENT
4347     C PSIG - PSIH FOR THE MOLECULAR LAYER
4348     C X - PHIM(ZETA) - DERIVATIVE OF PSIM
4349     C X0 - PHIM(ZETA0)
4350     C Y - PHIH(ZETA) - DERIVATIVE OF PSIH
4351     C Y0 - PHIH(ZETA0)
4352     C ITYPE - INTEGER FLAG :
4353     C 1 = NEUTRAL ADJUSTMENT
4354     C 2 = ADJ FOR 2ND OR GREATER TRBFLX ITER
4355     C 3 - 5 = ADJUSTMENT INSIDE LOOP
4356     C 4 - 5 = ADJUST CU AND CT
4357     C 5 = PREPARATION FOR ITYPE = 2
4358     C LWATER - LOGICAL - .TRUE. IF THERE ARE WATER POINTS
4359     C
4360     C OUTPUT:
4361     C -------
4362     C DZETA - D LOG ZETA
4363     C DZ0 - D Z0 (ITYPE 1) OR D LOG Z0 (ITYPE 2-5)
4364     C DPSIM - D PSIM
4365     C DPSIH - D PSIH
4366     C BITRIB - BIT ARRAY - '1' WHERE RIB1 = 0
4367     C
4368     C**********************************************************************
4369 molod 1.11 implicit none
4370    
4371     C Argument List Declarations
4372     integer nn,irun,itype
4373 molod 1.14 _RL VRIB1(IRUN),VRIB2(IRUN)
4374     _RL VWS1(IRUN),VWS2(IRUN),VZ1(IRUN),VUSTAR(IRUN)
4375 molod 1.11 integer IWATER(IRUN)
4376 molod 1.14 _RL VAPSIM(IRUN),VAPSIHG(IRUN)
4377     _RL VPSIH(IRUN),VPSIG(IRUN),VX(IRUN)
4378     _RL VX0(IRUN),VY(IRUN),VY0(IRUN)
4379 molod 1.11 LOGICAL LWATER
4380 molod 1.14 _RL VDZETA(IRUN),VDZ0(IRUN),VDPSIM(IRUN)
4381     _RL VDPSIH(IRUN)
4382 molod 1.11 integer INTRIB(IRUN)
4383 molod 1.14 _RL VX0PSIM(irun),VG(irun),VG0(irun),VR1MG0(irun)
4384     _RL VZ2(irun),VDZSEA(irun),VAZ0(irun),VXNUM1(irun)
4385     _RL VPSIGB2(irun),VDX(irun),VDXPSIM(irun),VDY(irun)
4386     _RL VXNUM2(irun),VDEN(irun),VAWS1(irun),VXNUM3(irun)
4387     _RL VXNUM(irun),VDZETA1(irun),VDZETA2(irun)
4388     _RL VZCOEF2(irun),VZCOEF1(irun),VTEMPLIN(irun)
4389     _RL VDPSIMC(irun),VDPSIHC(irun)
4390 molod 1.11
4391     C Local Variables
4392 molod 1.14 _RL xx0max,prfac,xpfac,difsqt,ustz0s,h0byz0,usth0s
4393 molod 1.1 PARAMETER ( XX0MAX = 1.49821 )
4394     PARAMETER ( PRFAC = 0.595864 )
4395     PARAMETER ( XPFAC = .55 )
4396     PARAMETER ( DIFSQT = 3.872983E-3)
4397     PARAMETER ( USTZ0S = 0.2030325E-5)
4398     PARAMETER ( H0BYZ0 = 30.0 )
4399     PARAMETER ( USTH0S = H0BYZ0*USTZ0S )
4400    
4401 molod 1.11 integer VINT1(irun),VINT2(irun)
4402 molod 1.14 _RL getcon,vk,bmdl,b2uhs
4403 molod 1.11 integer i
4404 molod 1.1 C
4405     vk = getcon('VON KARMAN')
4406     BMDL = VK * XPFAC * PRFAC / DIFSQT
4407     B2UHS = BMDL * BMDL * USTH0S
4408    
4409     C COMPUTE X0/PSIM, 1/Z0, G, G0, 1/(1-G0),
4410     C DEL LOG Z0, D LOG ZO / D USTAR
4411     C
4412     CCCOOOMMMM ADDED 'WHERE WATER'
4413     IF ( (ITYPE.EQ.1) .AND. LWATER ) THEN
4414     DO 9000 I = 1,IRUN
4415     IF (IWATER(I).EQ.1) VX0PSIM(I) = VAPSIM(I)
4416     9000 CONTINUE
4417     ENDIF
4418     IF ( ITYPE .GE. 3 ) THEN
4419     DO 9002 I = 1,IRUN
4420     VX0PSIM(I) = VX0(I) * VAPSIM(I)
4421     9002 CONTINUE
4422     ENDIF
4423     IF ( ITYPE .NE. 2 ) THEN
4424     C
4425     DO 9004 I = 1,IRUN
4426     VDZ0(I) = 0.
4427     VG(I) = 0.
4428     VG0(I) = 0.
4429     VR1MG0(I) = 1.
4430     9004 CONTINUE
4431     C
4432     IF ( LWATER ) THEN
4433     CALL ZCSUB ( VUSTAR,VDZSEA,IWATER,.TRUE.,IRUN,VZ2)
4434     C
4435     DO 9006 I = 1,IRUN
4436     IF ( IWATER(I).EQ.1) THEN
4437     VAZ0(I) = 1. / VZ1(I)
4438     VG(I) = VDZSEA(I) * VAZ0(I)
4439     VG0(I) = VX0PSIM(I) * VG(I)
4440     VR1MG0(I) = 1. / ( 1. - VG0(I) )
4441     VDZ0(I) = ( VZ2(I) - VZ1(I) ) * VR1MG0(I)
4442     ENDIF
4443     9006 CONTINUE
4444     ENDIF
4445     ENDIF
4446     C
4447     IF ( LWATER .AND. (ITYPE.GE.3) ) THEN
4448     DO 9008 I = 1,IRUN
4449     IF (IWATER(I).EQ.1) VDZ0(I) = VDZ0(I) * VAZ0(I)
4450     9008 CONTINUE
4451     ENDIF
4452     C
4453     C COMPUTE NUM1,NUM2,NUM3, DEN
4454     C
4455     IF (ITYPE.GE.3) THEN
4456     DO 9010 I = 1,IRUN
4457     VXNUM1(I) = 0.
4458     IF (VRIB1(I).EQ.0.) THEN
4459     INTRIB(I) = 1
4460     ELSE
4461     INTRIB(I) = 0
4462     ENDIF
4463     IF ( INTRIB(I).EQ.0 ) VXNUM1(I) = 1. / VRIB1(I)
4464     VPSIGB2(I) = 0.
4465     if(vpsig(i).gt.0.)VPSIGB2(I) =
4466     1 0.5 * ( vpsig(i)*vpsig(i) + b2uhs ) / vpsig(i)
4467     VDX(I) = VX(I) - VX0(I)
4468     VDXPSIM(I) = VDX(I) * VAPSIM(I)
4469     VDY(I) = VY(I) - VY0(I)
4470     VXNUM3(I) = - VPSIGB2(I)
4471     C
4472     IF ( LWATER ) THEN
4473     CCCOOOMMMM ADDED 'WHERE WATER'
4474     IF (IWATER(I).EQ.1) THEN
4475     VDXPSIM(I) = VDXPSIM(I) * VR1MG0(I)
4476     VXNUM3(I) = VXNUM3(I) + VG(I) * ( VY0(I) - VPSIGB2(I) )
4477     VXNUM2(I) = VY0(I) - VPSIGB2(I) - VX0PSIM(I) * VPSIGB2(I)
4478     VXNUM2(I) = (VXNUM2(I) * VAPSIHG(I)) - 2. * VX0PSIM(I)
4479     VXNUM2(I) = VXNUM2(I) * VDZ0(I)
4480     ENDIF
4481     ENDIF
4482     C
4483     VDEN(I) = VDY(I) + VDXPSIM(I) * VXNUM3(I)
4484     VDEN(I) = ( 1. + VDEN(I) * VAPSIHG(I) ) - 2. * VDXPSIM(I)
4485     9010 CONTINUE
4486     ENDIF
4487     C
4488     IF (ITYPE.EQ.5) THEN
4489     DO 9012 I = 1,IRUN
4490     VAWS1(I) = VR1MG0(I) / VWS1(I)
4491     VXNUM3(I) = VXNUM3(I) * VAPSIHG(I)
4492     C
4493     IF ( LWATER ) THEN
4494     CCCOOOMMMM ADDED 'WHERE WATER'
4495     IF(IWATER(I).EQ.1) THEN
4496     VXNUM3(I) = VXNUM3(I) - 2. * VG0(I)
4497     VXNUM3(I) = VAWS1(I) * VXNUM3(I)
4498     ENDIF
4499     ENDIF
4500     9012 CONTINUE
4501     ENDIF
4502     C
4503     C COMPUTE D LOG ZETA
4504     C
4505     IF (ITYPE.GE.2) THEN
4506     DO 9014 I = 1,IRUN
4507     VXNUM(I) = VRIB2(I) - VRIB1(I)
4508     IF( (VX0(I).GT.XX0MAX).AND.(VXNUM(I).GE.0.) )VXNUM(I) = 0.
4509     VXNUM(I) = VXNUM1(I) * VXNUM(I)
4510     9014 CONTINUE
4511     ENDIF
4512     C
4513     IF ( ITYPE.EQ.2 )THEN
4514     DO 9016 I = 1,IRUN
4515     VDZETA1(I) = VDZETA(I)
4516     VXNUM(I) = VXNUM(I) + VXNUM3(I) * ( VWS2(I) - VWS1(I) )
4517     9016 CONTINUE
4518     ENDIF
4519     C
4520     IF (ITYPE.GE.3) THEN
4521     DO 9018 I = 1,IRUN
4522     VDZETA1(I) = VXNUM(I)
4523     IF(LWATER.AND.(IWATER(I).EQ.1)) VXNUM(I) = VXNUM(I) + VXNUM2(I)
4524     IF ( VDEN(I) .LT.0.1 ) VDEN(I) = 0.1
4525     9018 CONTINUE
4526     ENDIF
4527     C
4528     IF (ITYPE.GE.2) THEN
4529     DO 9020 I = 1,IRUN
4530     VDZETA(I) = VXNUM(I) / VDEN(I)
4531     9020 CONTINUE
4532     ENDIF
4533     IF (ITYPE.GE.3) THEN
4534     DO 9022 I = 1,IRUN
4535     IF ( (VRIB2(I).EQ.0.) .OR. (VDZETA(I).LE.-1.) ) THEN
4536     VINT1(I) = 1
4537     ELSE
4538     VINT1(I) = 0
4539     ENDIF
4540     IF ( VINT1(I).EQ.1 ) VDZETA(I) = VDZETA1(I)
4541     9022 CONTINUE
4542     ENDIF
4543     IF (ITYPE.EQ.2) THEN
4544     DO 9024 I = 1,IRUN
4545     VDZETA2(I) = VDZETA(I) + VDZETA1(I)
4546     IF ( (VRIB2(I).EQ.0.) .OR. (VDZETA2(I).LE.-1.) ) THEN
4547     VINT1(I) = 1
4548     ELSE
4549     VINT1(I) = 0
4550     ENDIF
4551     IF(VINT1(I).EQ.1)VDZETA(I)=VXNUM1(I)*VRIB2(I) - 1. - VDZETA1(I)
4552     9024 CONTINUE
4553     ENDIF
4554    
4555     C
4556     C COMPUTE D LOG Z0
4557     C
4558     IF ( LWATER .AND. (ITYPE.GE.3) )THEN
4559     DO 9026 I = 1,IRUN
4560     IF( IWATER(I).EQ.1 ) THEN
4561     VZCOEF2(I) = VG(I) * VDXPSIM(I)
4562     VDZ0(I) = VDZ0(I) - VZCOEF2(I) * VDZETA(I)
4563     ENDIF
4564     9026 CONTINUE
4565     ENDIF
4566     C
4567     IF ( LWATER .AND. (ITYPE.EQ.5) ) THEN
4568     DO 9028 I = 1,IRUN
4569     IF(IWATER(I).EQ.1) VZCOEF1(I) = VG(I) * VAWS1(I)
4570     9028 CONTINUE
4571     ENDIF
4572     C
4573     IF ( LWATER .AND. (ITYPE.EQ.2) ) THEN
4574     DO 9030 I = 1,IRUN
4575     IF (IWATER(I).EQ.1) VDZ0(I) =
4576     1 VZCOEF1(I) * ( VWS2(I) - VWS1(I) ) - VZCOEF2(I) * VDZETA(I)
4577     9030 CONTINUE
4578     ENDIF
4579     C
4580     C CALCULATE D PSIM AND D PSIH
4581     C
4582     IF ( (ITYPE.EQ.1) .AND. LWATER ) THEN
4583     DO 9032 I = 1,IRUN
4584     IF (IWATER(I).EQ.1) THEN
4585     VDPSIM(I) = - VDZ0(I) * VAZ0(I)
4586     VDPSIH(I) = VDPSIM(I)
4587     ENDIF
4588     9032 CONTINUE
4589     ENDIF
4590     C
4591     IF (ITYPE.GE.3) THEN
4592     DO 9034 I = 1,IRUN
4593     VDPSIM(I) = VDX(I) * VDZETA(I)
4594     VDPSIH(I) = VDY(I) * VDZETA(I)
4595     IF ( LWATER ) THEN
4596     IF (IWATER(I).EQ.1 ) THEN
4597     VDPSIM(I) = VDPSIM(I) - VX0(I) * VDZ0(I)
4598     VDPSIH(I) = VDPSIH(I) - VY0(I) * VDZ0(I)
4599     ENDIF
4600     ENDIF
4601     9034 CONTINUE
4602     ENDIF
4603     C
4604     C PREVENT OVERCORRECTION OF PSIM OR PSIH FOR UNSTABLE CASE
4605     C
4606     IF (ITYPE.GE.4) THEN
4607     DO 9036 I = 1,IRUN
4608     VDPSIMC(I) = -0.9 - VDPSIM(I) * VAPSIM(I)
4609     VDPSIHC(I) = -0.9 * VPSIH(I) - VDPSIH(I)
4610     IF ( VDPSIMC(I).GT.0. ) THEN
4611     VINT1(I) = 1
4612     ELSE
4613     VINT1(I) = 0
4614     ENDIF
4615     IF ( VDPSIHC(I).GT.0. ) THEN
4616     VINT2(I) = 1
4617     ELSE
4618     VINT2(I) = 0
4619     ENDIF
4620     VDZETA1(I) = 0.
4621     IF(VINT1(I).EQ.1) VDZETA1(I) = VDPSIMC(I) / VDXPSIM(I)
4622     IF((VINT1(I).EQ.1).OR.(VINT2(I).EQ.1)) VTEMPLIN(I) =
4623     1 VDY(I) + VY0(I) * VG(I) * VDXPSIM(I)
4624     IF (VINT2(I).EQ.1) then
4625     VDZETA2(I) = VDPSIHC(I) / VTEMPLIN(I)
4626     IF ( VDZETA2(I).LT.VDZETA1(I) ) VDZETA1(I) = VDZETA2(I)
4627     endif
4628     IF((VINT1(I).EQ.1).OR.(VINT2(I).EQ.1)) THEN
4629     VDZETA(I) = VDZETA1(I) + VDZETA(I)
4630     VDPSIM(I) = VDPSIM(I) + VDX(I) * VR1MG0(I) * VDZETA1(I)
4631     VDPSIH(I) = VDPSIH(I) + VTEMPLIN(I) * VDZETA1(I)
4632     IF ( IWATER(I).EQ.1 )
4633     1 VDZ0(I) = VDZ0(I) - VG(I) * VDXPSIM(I) * VDZETA1(I)
4634     ENDIF
4635     9036 CONTINUE
4636     ENDIF
4637     C
4638     RETURN
4639     END
4640     SUBROUTINE ZCSUB (VUSTAR,VDZSEA,IWATER,LDZSEA,IRUN,VZSEA)
4641     C**********************************************************************
4642     C FUNCTION ZSEA
4643     C PURPOSE
4644     C COMPUTES Z0 AS A FUNCTION OF USTAR OVER WATER SURFACES
4645     C USAGE
4646     C CALLED BY SFCFLX
4647     C DESCRIPTION OF PARAMETERS
4648     C USTAR - INPUTED VALUE OF SURFACE-STRESS VELOCITY
4649     C DZSEA - OUTPUTED VALUE OF DERIVATIVE D(ZSEA)/D(USTAR)
4650     C WATER - INPUTED BIT VECTOR TO DETERMINE WATER POINTS
4651     C LDZSEA- LOGICAL FLAG TO DETERMINE IF DZSEA SHOULD BE COMPUTED
4652     C ZSEA - OUTPUTED VALUE OF ROUGHNESS LENGTH
4653     C SUBPROGRAMS NEEDED
4654     C NONE
4655     C RECORD OF MODIFICATIONS
4656     C REMARKS:
4657     C COMPUTE ROUGHNESS LENGTH FOR OCEAN POINTS
4658     C BASED ON FUNCTIONS OF LARGE AND POND
4659     C AND OF KONDO --- DESIGNED FOR K = .4
4660     C *********************************************************************
4661 molod 1.11 implicit none
4662    
4663     C Argument List Delcarations
4664     integer irun
4665 molod 1.14 _RL VZSEA(IRUN),VUSTAR(IRUN),VDZSEA(IRUN)
4666 molod 1.11 integer IWATER(IRUN)
4667     LOGICAL LDZSEA
4668    
4669     C Local Variables
4670 molod 1.14 _RL USTMX1,USTMX2,USTMX3
4671 molod 1.1 PARAMETER ( USTMX1 = 1.14973 )
4672     PARAMETER ( USTMX2 = 0.381844 )
4673     PARAMETER ( USTMX3 = 0.0632456)
4674    
4675 molod 1.14 _RL AA(IRUN,5),TEMP(IRUN)
4676 molod 1.11 integer INT2(IRUN),INT3(IRUN),INT4(IRUN)
4677     integer i,k
4678    
4679 molod 1.14 _RL AA1(5),AA2(5),AA3(5),AA4(5)
4680 molod 1.1 DATA AA1/.2030325E-5,0.0,0.0,0.0,0.0/
4681     DATA AA2/-0.402451E-08,0.239597E-04,0.117484E-03,0.191918E-03,
4682     1 0.395649E-04/
4683     DATA AA3/-0.237910E-04,0.228221E-03,-0.860810E-03,0.176543E-02,
4684     1 0.784260E-04/
4685     DATA AA4/-0.343228E-04,0.552305E-03,-0.167541E-02,0.250208E-02,
4686     1 -0.153259E-03/
4687     C
4688     C**********************************************************************
4689     C***** LOWER CUTOFF CONDITION FOR USTAR ***
4690     C**********************************************************************
4691     C
4692     DO 9000 I = 1,IRUN
4693     IF(VUSTAR(I) .LT. 1.e-6)THEN
4694     INT3(I) = 1
4695     ELSE
4696     INT3(I) = 0
4697     ENDIF
4698     9000 CONTINUE
4699     DO 9002 I = 1,IRUN
4700     IF(INT3(I).EQ.1) VUSTAR(I) = 1.e-6
4701     9002 CONTINUE
4702     C
4703     C***********************************
4704     C***** LOAD THE ARRAY A(I,K) *****
4705     C***********************************
4706     C
4707     DO 9004 I = 1,IRUN
4708     IF( (VUSTAR(I) .GT. USTMX1) .AND. (IWATER(I).EQ.1) ) THEN
4709     INT4(I) = 1
4710     ELSE
4711     INT4(I) = 0
4712     ENDIF
4713     9004 CONTINUE
4714     DO 9006 I = 1,IRUN
4715     IF(VUSTAR(I) .GT. USTMX2) THEN
4716     INT3(I) = 1
4717     ELSE
4718     INT3(I) = 0
4719     ENDIF
4720     9006 CONTINUE
4721     DO 9008 I = 1,IRUN
4722     IF(VUSTAR(I) .GE. USTMX3) THEN
4723     INT2(I) = 1
4724     ELSE
4725     INT2(I) = 0
4726     ENDIF
4727     9008 CONTINUE
4728     C
4729     DO 100 K=1,5
4730     DO 9010 I = 1,IRUN
4731     AA(I,K) = AA1(K)
4732     IF( INT2(I).EQ.1 ) AA(I,K) = AA2(K)
4733     IF( INT3(I).EQ.1 ) AA(I,K) = AA3(K)
4734     IF( INT4(I).EQ.1 ) AA(I,K) = AA4(K)
4735     9010 CONTINUE
4736     100 CONTINUE
4737     C
4738     C********************************************************
4739     C***** EVALUATE THE ENHANCED POLYNOMIAL FOR ZSEA *****
4740     C********************************************************
4741     C
4742     DO 9012 I = 1,IRUN
4743     VDZSEA(I) = ( AA(I,4) + AA(I,5) * VUSTAR(I) ) * VUSTAR(I)
4744     VZSEA(I) = AA(I,2) + ( AA(I,3) + VDZSEA(I) ) * VUSTAR(I)
4745     TEMP(I) = AA(I,1) / VUSTAR(I)
4746     VZSEA(I) = VZSEA(I) + TEMP(I)
4747     9012 CONTINUE
4748     C
4749     C**********************************************************************
4750     C***** EVALUATE THE DERIVATIVE DZSEA IF LDZSEA IS TRUE ***
4751     C**********************************************************************
4752     C
4753     IF( LDZSEA ) THEN
4754     DO 9014 I = 1,IRUN
4755     VDZSEA(I) = 3. * VDZSEA(I) -(AA(I,4)*VUSTAR(I) - AA(I,3))
4756     VDZSEA(I) = VDZSEA(I) * VUSTAR(I) - TEMP(I)
4757     9014 CONTINUE
4758     ENDIF
4759     C
4760     RETURN
4761     END
4762    
4763     subroutine seaice ( nocean, timstp, hice,
4764     . eturb, dedtc,
4765     . hsturb, dhsdtc,
4766     . qice, dqice,
4767     . swnet, lwnet, dst4,
4768     . pke, seaic, tc, qa )
4769     implicit none
4770     integer nocean
4771 molod 1.14 _RL timstp
4772     _RL eturb(nocean),dedtc(nocean),hsturb(nocean),dhsdtc(nocean)
4773     _RL swnet(nocean),lwnet(nocean), dst4(nocean)
4774     _RL qice(nocean),dqice(nocean)
4775     _RL pke(nocean), tc(nocean), qa(nocean)
4776     _RL seaic(nocean)
4777 molod 1.1
4778     C rho*C = 1.93e6 J/(m**3 K) ; Peixoto & Oort
4779 molod 1.14 _RL rhoC
4780 molod 1.11 parameter (rhoC = 1.93e6)
4781 molod 1.1
4782 molod 1.14 _RL faceps,getcon,latent,codt,deltg,hice
4783 molod 1.1 integer i
4784    
4785     faceps = getcon('EPSFAC')
4786     latent = getcon('HEATI') * getcon('CALTOJ')
4787     C Note hice is in centimeters
4788     codt = rhoC * (hice/100) / timstp
4789    
4790     c Update TC and QA
4791     c ----------------
4792     do i =1,nocean
4793     if( seaic(i).gt.0.0 ) then
4794     deltg = ( swnet(i)-lwnet(i)-latent*eturb(i)-hsturb(i)+qice(i) )
4795     . / ( codt+dst4(i)+latent*dedtc(i)+dhsdtc(i)-dqice(i) )
4796     qa(i) = qa(i) + (faceps*qa(i)/(tc(i)*tc(i)))*deltg
4797     tc(i) = tc(i) + deltg
4798     endif
4799     enddo
4800    
4801     return
4802     end

  ViewVC Help
Powered by ViewVC 1.1.22