/[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.18 - (hide annotations) (download)
Wed Jul 28 21:57:51 2004 UTC (19 years, 11 months ago) by molod
Branch: MAIN
Changes since 1.17: +6 -17 lines
Debugging - it runs!

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

  ViewVC Help
Powered by ViewVC 1.1.22