/[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.50 - (show annotations) (download)
Tue Mar 27 15:48:27 2012 UTC (12 years, 1 month ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint64, checkpoint65, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65o, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, HEAD
Changes since 1.49: +7 -42 lines
clean-up turbulence cold-start switch: decided in fizhi_init_vars.F, stored
in common bloc (fizhi_coms.h) and then passed as argument up to S/R TURBIO.

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

  ViewVC Help
Powered by ViewVC 1.1.22