/[MITgcm]/MITgcm/verification/fizhi-cs-aqualev20/code/fizhi_turb.F
ViewVC logotype

Contents of /MITgcm/verification/fizhi-cs-aqualev20/code/fizhi_turb.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (show annotations) (download)
Mon Apr 3 20:55:14 2006 UTC (18 years, 1 month ago) by molod
Branch: MAIN
CVS Tags: checkpoint58r_post, checkpoint58l_post, checkpoint58k_post, checkpoint58e_post, checkpoint58n_post, checkpoint58v_post, checkpoint58h_post, checkpoint58s_post, checkpoint58p_post, checkpoint58f_post, checkpoint58t_post, checkpoint58m_post, checkpoint58q_post, checkpoint58j_post, checkpoint58d_post, checkpoint58i_post, checkpoint58u_post, checkpoint58g_post, checkpoint58o_post
New check-in for fizhi aquaplanet experiment to match the APE control (will replace fizhi-cs-aqualev10)

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

  ViewVC Help
Powered by ViewVC 1.1.22