/[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.12 - (show annotations) (download)
Fri Jul 16 20:08:08 2004 UTC (19 years, 11 months ago) by molod
Branch: MAIN
CVS Tags: checkpoint54c_post
Changes since 1.11: +5 -5 lines
Cleaning

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

  ViewVC Help
Powered by ViewVC 1.1.22