/[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.23 - (show annotations) (download)
Tue Aug 10 15:13:47 2004 UTC (19 years, 10 months ago) by molod
Branch: MAIN
Changes since 1.22: +540 -404 lines
Debugging

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

  ViewVC Help
Powered by ViewVC 1.1.22