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

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

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


Revision 1.37 - (hide annotations) (download)
Thu May 5 14:41:53 2005 UTC (19 years, 1 month ago) by molod
Branch: MAIN
Changes since 1.36: +9 -1 lines
Check in fix for wierd behavior aloft

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

  ViewVC Help
Powered by ViewVC 1.1.22