/[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.19 - (hide annotations) (download)
Thu Jul 29 15:21:50 2004 UTC (19 years, 11 months ago) by molod
Branch: MAIN
Changes since 1.18: +4 -1 lines
Again, integer versus real time step

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

  ViewVC Help
Powered by ViewVC 1.1.22