/[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.32 - (hide annotations) (download)
Thu Feb 24 23:26:03 2005 UTC (19 years, 4 months ago) by molod
Branch: MAIN
Changes since 1.31: +320 -514 lines
Changes to get rid of diagnostics common blocks - almost there

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

  ViewVC Help
Powered by ViewVC 1.1.22