/[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.25 - (show annotations) (download)
Thu Aug 26 14:41:56 2004 UTC (19 years, 9 months ago) by molod
Branch: MAIN
CVS Tags: checkpoint54e_post
Changes since 1.24: +9 -11 lines
Get rid of compiler warming messages about intrinsic internal functions

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

  ViewVC Help
Powered by ViewVC 1.1.22