/[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.30 - (show annotations) (download)
Tue Feb 22 22:51:42 2005 UTC (19 years, 4 months ago) by molod
Branch: MAIN
Changes since 1.29: +26 -22 lines
More diagnostic calls

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

  ViewVC Help
Powered by ViewVC 1.1.22