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

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

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


Revision 1.7 - (show annotations) (download)
Mon Jul 12 21:33:37 2004 UTC (19 years, 11 months ago) by molod
Branch: MAIN
Changes since 1.6: +27 -40 lines
Corrections for argument list match and proper declarations

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

  ViewVC Help
Powered by ViewVC 1.1.22