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

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

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


Revision 1.2 - (show annotations) (download)
Tue Jun 15 14:47:24 2004 UTC (19 years, 11 months ago) by molod
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +0 -0 lines
FILE REMOVED
Add CVS header and name stuff (and rename)

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

  ViewVC Help
Powered by ViewVC 1.1.22