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

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

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


Revision 1.45 - (show annotations) (download)
Sun Jul 8 20:24:02 2007 UTC (16 years, 11 months ago) by molod
Branch: MAIN
CVS Tags: checkpoint60, checkpoint61, checkpoint59q, checkpoint59p, checkpoint59r, checkpoint59e, checkpoint59g, checkpoint59f, checkpoint59m, checkpoint59l, checkpoint59o, checkpoint59n, checkpoint59i, checkpoint59h, checkpoint59k, checkpoint59j, checkpoint61f, checkpoint61g, checkpoint61d, checkpoint61e, checkpoint61b, checkpoint61c, checkpoint61a, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i
Changes since 1.44: +2 -3 lines
Some stuff for the strat version

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

  ViewVC Help
Powered by ViewVC 1.1.22