/[MITgcm]/MITgcm/pkg/fizhi/fizhi_turb.F
ViewVC logotype

Contents of /MITgcm/pkg/fizhi/fizhi_turb.F

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


Revision 1.48 - (show annotations) (download)
Tue Jan 19 14:24:36 2010 UTC (14 years, 5 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63g, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint63, checkpoint62c, checkpoint62b, checkpoint62g, checkpoint62f, checkpoint62e, checkpoint62d, checkpoint62k, checkpoint62j, checkpoint62i, checkpoint62h, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62p, checkpoint62w, checkpoint62v, checkpoint62u, checkpoint62t, checkpoint62z, checkpoint62y, checkpoint62x
Changes since 1.47: +4 -4 lines
add 3 "_d 0"s so that fizhi compiles with xlf on iblade

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

  ViewVC Help
Powered by ViewVC 1.1.22