CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
ysuvdif.F90
1
9!----------------------------------------------------------------------
10
11 module ysuvdif
12 contains
13
14 subroutine ysuvdif_init (do_ysu,errmsg,errflg)
15
16 logical, intent(in) :: do_ysu
17 character(len=*), intent(out) :: errmsg
18 integer, intent(out) :: errflg
19
20 ! Initialize CCPP error handling variables
21 errmsg = ''
22 errflg = 0
23
24 ! Consistency checks
25 if (.not. do_ysu) then
26 write(errmsg,fmt='(*(a))') 'Logic error: do_ysu = .false.'
27 errflg = 1
28 return
29 end if
30 end subroutine ysuvdif_init
31
39!-------------------------------------------------------------------------------
40 subroutine ysuvdif_run(im,km,ux,vx,tx,qx,p2d,p2di,pi2d,karman, &
41 utnp,vtnp,ttnp,qtnp, &
42 swh,hlw,xmu,ntrac,ndiff,ntcw,ntiw, &
43 phii,phil,psfcpa, &
44 zorl,stress,hpbl,psim,psih, &
45 landmask,heat,evap,wspd,br, &
46 g,rd,cp,rv,ep1,ep2,xlv, &
47 dusfc,dvsfc,dtsfc,dqsfc, &
48 dt,kpbl1d,u10,v10,lssav,ldiag3d,qdiag3d, &
49 flag_for_pbl_generic_tend,ntoz,ntqv,dtend,dtidx, &
50 index_of_temperature,index_of_x_wind,index_of_y_wind, &
51 index_of_process_pbl,errmsg,errflg )
52
53 use machine , only : kind_phys
54!
55!-------------------------------------------------------------------------------
56 implicit none
57!-------------------------------------------------------------------------------
58 real(kind=kind_phys),parameter :: xkzminm = 0.1,xkzminh = 0.01
59 real(kind=kind_phys),parameter :: xkzmin = 0.01,xkzmax = 1000.,rimin = -100.
60 real(kind=kind_phys),parameter :: rlam = 30.,prmin = 0.25,prmax = 4.
61 real(kind=kind_phys),parameter :: brcr_ub = 0.0,brcr_sb = 0.25,cori = 1.e-4
62 real(kind=kind_phys),parameter :: afac = 6.8,bfac = 6.8,pfac = 2.0,pfac_q = 2.0
63 real(kind=kind_phys),parameter :: phifac = 8.,sfcfrac = 0.1
64 real(kind=kind_phys),parameter :: d1 = 0.02, d2 = 0.05, d3 = 0.001
65 real(kind=kind_phys),parameter :: h1 = 0.33333335, h2 = 0.6666667
66 real(kind=kind_phys),parameter :: zfmin = 1.e-8,aphi5 = 5.,aphi16 = 16.
67 real(kind=kind_phys),parameter :: tmin=1.e-2
68 real(kind=kind_phys),parameter :: gamcrt = 3.,gamcrq = 2.e-3
69 real(kind=kind_phys),parameter :: xka = 2.4e-5
70 real(kind=kind_phys),parameter :: rcl = 1.0
71 real(kind=kind_phys),intent(in) :: karman
72 integer,parameter :: imvdif = 1
73 integer,parameter :: ysu_topdown_pblmix = 1
74!
75!-------------------------------------------------------------------------------------
76! input variables
77 integer, intent(in ) :: im,km,ntrac,ndiff,ntcw,ntiw,ntoz
78 real(kind=kind_phys), intent(in ) :: g,cp,rd,rv,ep1,ep2,xlv,dt
79
80 real(kind=kind_phys), dimension( :,: ), &
81 intent(in) :: pi2d,p2d,phil,ux,vx,swh,hlw,tx
82
83 real(kind=kind_phys), dimension( :,:,: ) , &
84 intent(in ) :: qx
85
86 real(kind=kind_phys), dimension( :,: ) , &
87 intent(in ) :: p2di,phii
88
89 real(kind=kind_phys), dimension( : ) , &
90 intent(in) :: stress,zorl,heat,evap,wspd,br,psim,psih,psfcpa, &
91 u10,v10,xmu
92 integer, dimension(:) ,&
93 intent(in ) :: landmask
94 logical, intent(in ) :: lssav, ldiag3d, qdiag3d, &
95 flag_for_pbl_generic_tend
96!
97!----------------------------------------------------------------------------------
98! input/output variables
99!
100 real(kind=kind_phys), dimension( :,: ) , &
101 intent(inout) :: utnp,vtnp,ttnp
102 real(kind=kind_phys), dimension( :,:,: ) , &
103 intent(inout) :: qtnp
104 real(kind=kind_phys), optional, intent(inout), optional :: dtend(:,:,:)
105 integer, intent(in) :: dtidx(:,:), ntqv, index_of_temperature, &
106 index_of_x_wind, index_of_y_wind, index_of_process_pbl
107!
108!---------------------------------------------------------------------------------
109! output variables
110 integer, dimension( : ), intent(out ) :: kpbl1d
111 real(kind=kind_phys), dimension( : ), &
112 intent(out) :: hpbl
113 real(kind=kind_phys), dimension( : ), &
114 intent(out) :: dusfc,dvsfc, dtsfc,dqsfc
115
116 ! error messages
117 character(len=*), intent(out) :: errmsg
118 integer, intent(out) :: errflg
119!
120!--------------------------------------------------------------------------------
121!
122! local vars
123!
124 real(kind=kind_phys), dimension( im ) :: hol
125 real(kind=kind_phys), dimension( im, km+1 ) :: zq
126!
127 real(kind=kind_phys), dimension( im, km ) :: &
128 thx,thvx,thlix, &
129 del, &
130 dza, &
131 dzq, &
132 xkzom, &
133 xkzoh, &
134 za
135!
136 real(kind=kind_phys), dimension( im ) :: &
137 rhox, &
138 govrth, &
139 zl1,thermal, &
140 wscale, &
141 hgamt,hgamq, &
142 brdn,brup, &
143 phim,phih, &
144 prpbl, &
145 wspd1,thermalli
146!
147 real(kind=kind_phys), dimension( im, km ) :: xkzm,xkzh, &
148 f1,f2, &
149 r1,r2, &
150 ad,au, &
151 cu, &
152 al, &
153 xkzq, &
154 zfac, &
155 rhox2, &
156 hgamt2
157!
158 real(kind=kind_phys), dimension( im ) :: &
159 brcr, &
160 sflux, &
161 zol1, &
162 brcr_sbro
163!
164 real(kind=kind_phys), dimension( im ) :: xland
165 real(kind=kind_phys), dimension( im ) :: ust
166 real(kind=kind_phys), dimension( im ) :: hfx
167 real(kind=kind_phys), dimension( im ) :: qfx
168 real(kind=kind_phys), dimension( im ) :: znt
169 real(kind=kind_phys), dimension( im ) :: uox
170 real(kind=kind_phys), dimension( im ) :: vox
171!
172 real(kind=kind_phys), dimension( im, km, ndiff) :: r3,f3
173 integer, dimension( im ) :: kpbl,kpblold
174!
175 logical, dimension( im ) :: pblflg, &
176 sfcflg, &
177 stable, &
178 cloudflg
179
180 logical :: definebrup
181!
182 integer :: n,i,k,l,ic,is,kk
183 integer :: klpbl, ktrace1, ktrace2, ktrace3
184!
185!
186 real(kind=kind_phys) :: dt2,rdt,spdk2,fm,fh,hol1,gamfac,vpert,prnum,prnum0
187 real(kind=kind_phys) :: ss,ri,qmean,tmean,alph,chi,zk,rl2,dk,sri
188 real(kind=kind_phys) :: brint,dtodsd,dtodsu,rdz,dsdzt,dsdzq,dsdz2,rlamdz
189 real(kind=kind_phys) :: utend,vtend,ttend,qtend
190 real(kind=kind_phys) :: dtstep,govrthv
191 real(kind=kind_phys) :: cont, conq, conw, conwrc, rovcp
192!
193
194 real(kind=kind_phys), dimension( im, km ) :: wscalek,wscalek2
195 real(kind=kind_phys), dimension( im ) :: wstar
196 real(kind=kind_phys), dimension( im ) :: delta
197 real(kind=kind_phys), dimension( im, km ) :: xkzml,xkzhl, &
198 zfacent,entfac
199 real(kind=kind_phys), dimension( im ) :: ust3, &
200 wstar3, &
201 wstar3_2, &
202 hgamu,hgamv, &
203 wm2, we, &
204 bfxpbl, &
205 hfxpbl,qfxpbl, &
206 ufxpbl,vfxpbl, &
207 dthvx
208 real(kind=kind_phys) :: prnumfac,bfx0,hfx0,qfx0,delb,dux,dvx, &
209 dsdzu,dsdzv,wm3,dthx,dqx,wspd10,ross,tem1,dsig,tvcon,conpr, &
210 prfac,prfac2,phim8z,radsum,tmp1,templ,rvls,temps,ent_eff, &
211 rcldb,bruptmp,radflux
212 integer :: idtend
213!
214!-------------------------------------------------------------------------------
215!
216! Initialize CCPP error handling variables
217 errmsg = ''
218 errflg = 0
219
220 klpbl = km
221!
222 rovcp=rd/cp
223 cont=cp/g
224 conq=xlv/g
225 conw=1./g
226 conwrc = conw*sqrt(rcl)
227 conpr = bfac*karman*sfcfrac
228!
229! change xland values
230 do i=1,im
231 if(landmask(i).eq.0) then !ocean
232 xland(i) = 2
233 else
234 xland(i) = 1 !land
235 end if
236 end do
237!
238 do k = 1,km
239 do i = 1,im
240 thx(i,k) = tx(i,k)/pi2d(i,k)
241 thlix(i,k) = (tx(i,k)-xlv*qx(i,k,ntcw)/cp-2.834e6*qx(i,k,ntiw)/cp)/pi2d(i,k)
242 enddo
243 enddo
244!
245 do k = 1,km
246 do i = 1,im
247 tvcon = (1.+ep1*qx(i,k,1))
248 thvx(i,k) = thx(i,k)*tvcon
249 enddo
250 enddo
251!
252 do i = 1,im
253 tvcon = (1.+ep1*qx(i,1,1))
254 rhox(i) = psfcpa(i)/(rd*tx(i,1)*tvcon)
255 govrth(i) = g/thx(i,1)
256 hfx(i) = heat(i)*rhox(i)*cp ! reset to the variable in WRF
257 qfx(i) = evap(i)*rhox(i) ! reset to the variable in WRF
258 ust(i) = sqrt(stress(i)) ! reset to the variable in WRF
259 znt(i) = 0.01*zorl(i) ! reset to the variable in WRF
260 uox(i) = 0.0
261 vox(i) = 0.0
262 enddo
263!
264!-----compute the height of full- and half-sigma levels above ground
265! level, and the layer thicknesses.
266!
267 do i = 1,im
268 zq(i,1) = 0.
269 enddo
270!
271 do k = 1,km
272 do i = 1,im
273 zq(i,k+1) = phii(i,k+1)*conw
274 tvcon = (1.+ep1*qx(i,k,1))
275 rhox2(i,k) = p2d(i,k)/(rd*tx(i,k)*tvcon)
276 enddo
277 enddo
278!
279 do k = 1,km
280 do i = 1,im
281 za(i,k) = phil(i,k)*conw
282 dzq(i,k) = zq(i,k+1)-zq(i,k)
283 del(i,k) = p2di(i,k)-p2di(i,k+1)
284 enddo
285 enddo
286!
287 do i = 1,im
288 dza(i,1) = za(i,1)
289 enddo
290!
291 do k = 2,km
292 do i = 1,im
293 dza(i,k) = za(i,k)-za(i,k-1)
294 enddo
295 enddo
296
297! write(0,*)"===CALLING ysu; input:"
298! print*,"t:",tx(1,1),tx(1,2),tx(1,km)
299! print*,"u:",ux(1,1),ux(1,2),ux(1,km)
300! print*,"v:",vx(1,1),vx(1,2),vx(1,km)
301! print*,"q:",qx(1,1,1),qx(1,2,1),qx(1,km,1)
302! print*,"exner:",pi2d(1,1),pi2d(1,2),pi2d(1,km)
303! print*,"phii:",zq(1,1),zq(1,2),zq(1,km+1)
304! print*,"phil:",za(1,1),za(1,2),za(1,km)
305! print*,"p2d:",p2d(1,1),p2d(1,2),p2d(1,km)
306! print*,"p2di:",p2di(1,1),p2di(1,2),p2di(1,km+1)
307! print *,"del:",del(1,1),del(1,2),del(1,km)
308! print*,"znt,ust,wspd:",znt(1),ust(1),wspd(1)
309! print*,"hfx,qfx,xland:",hfx(1),qfx(1),xland(1)
310! print*,"rd,rv,g:",rd,rv,g
311! print*,"ep1,ep2,xlv:",ep1,ep2,xlv
312! print*,"br,psim,psih:",br(1),psim(1),psih(1)
313! print*,"u10,v10:",u10(1),v10(1)
314! print*,"psfcpa,cp:",psfcpa(1),cp
315! print*,"ntrac,ndiff,ntcw,ntiw:",ntrac,ndiff,ntcw,ntiw
316!
317!
318!-----initialize vertical tendencies and
319!
320! utnp(:,:) = 0.
321! vtnp(:,:) = 0.
322! ttnp(:,:) = 0.
323! qtnp(:,:,:) = 0.
324!
325 do i = 1,im
326 wspd1(i) = sqrt( (ux(i,1)-uox(i))*(ux(i,1)-uox(i)) + (vx(i,1)-vox(i))*(vx(i,1)-vox(i)) )+1.e-9
327 enddo
328!
329!---- compute vertical diffusion
330!
331! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
332! compute preliminary variables
333!
334 dtstep = dt
335 dt2 = 2.*dtstep
336 rdt = 1./dt2
337!
338 do i = 1,im
339 bfxpbl(i) = 0.0
340 hfxpbl(i) = 0.0
341 qfxpbl(i) = 0.0
342 ufxpbl(i) = 0.0
343 vfxpbl(i) = 0.0
344 hgamu(i) = 0.0
345 hgamv(i) = 0.0
346 delta(i) = 0.0
347 wstar3_2(i) = 0.0
348 enddo
349!
350 do k = 1,klpbl
351 do i = 1,im
352 wscalek(i,k) = 0.0
353 wscalek2(i,k) = 0.0
354 enddo
355 enddo
356!
357 do k = 1,klpbl
358 do i = 1,im
359 zfac(i,k) = 0.0
360 enddo
361 enddo
362 do k = 1,klpbl-1
363 do i = 1,im
364 xkzom(i,k) = xkzminm
365 xkzoh(i,k) = xkzminh
366 enddo
367 enddo
368!
369 do i = 1,im
370 dusfc(i) = 0.
371 dvsfc(i) = 0.
372 dtsfc(i) = 0.
373 dqsfc(i) = 0.
374 enddo
375!
376 do i = 1,im
377 hgamt(i) = 0.
378 hgamq(i) = 0.
379 wscale(i) = 0.
380 kpbl(i) = 1
381 hpbl(i) = zq(i,1)
382 zl1(i) = za(i,1)
383 thermal(i)= thvx(i,1)
384 thermalli(i) = thlix(i,1)
385 pblflg(i) = .true.
386 sfcflg(i) = .true.
387 sflux(i) = hfx(i)/rhox(i)/cp + qfx(i)/rhox(i)*ep1*thx(i,1)
388 if(br(i).gt.0.0) sfcflg(i) = .false.
389 enddo
390!
391! compute the first guess of pbl height
392!
393 do i = 1,im
394 stable(i) = .false.
395 brup(i) = br(i)
396 brcr(i) = brcr_ub
397 enddo
398!
399 do k = 2,klpbl
400 do i = 1,im
401 if(.not.stable(i))then
402 brdn(i) = brup(i)
403 spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.)
404 brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2
405 kpbl(i) = k
406 stable(i) = brup(i).gt.brcr(i)
407 endif
408 enddo
409 enddo
410!
411 do i = 1,im
412 k = kpbl(i)
413 if(brdn(i).ge.brcr(i))then
414 brint = 0.
415 elseif(brup(i).le.brcr(i))then
416 brint = 1.
417 else
418 brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i))
419 endif
420 hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
421 if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
422 if(kpbl(i).le.1) pblflg(i) = .false.
423 enddo
424!
425 do i = 1,im
426 fm = psim(i)
427 fh = psih(i)
428 zol1(i) = max(br(i)*fm*fm/fh,rimin)
429 if(sfcflg(i))then
430 zol1(i) = min(zol1(i),-zfmin)
431 else
432 zol1(i) = max(zol1(i),zfmin)
433 endif
434 hol1 = zol1(i)*hpbl(i)/zl1(i)*sfcfrac
435 if(sfcflg(i))then
436 phim(i) = (1.-aphi16*hol1)**(-1./4.)
437 phih(i) = (1.-aphi16*hol1)**(-1./2.)
438 bfx0 = max(sflux(i),0.)
439 hfx0 = max(hfx(i)/rhox(i)/cp,0.)
440 qfx0 = max(ep1*thx(i,1)*qfx(i)/rhox(i),0.)
441 wstar3(i) = (govrth(i)*bfx0*hpbl(i))
442 wstar(i) = (wstar3(i))**h1
443 else
444 phim(i) = (1.+aphi5*hol1)
445 phih(i) = phim(i)
446 wstar(i) = 0.
447 wstar3(i) = 0.
448 endif
449 ust3(i) = ust(i)**3.
450 wscale(i) = (ust3(i)+phifac*karman*wstar3(i)*0.5)**h1
451 wscale(i) = min(wscale(i),ust(i)*aphi16)
452 wscale(i) = max(wscale(i),ust(i)/aphi5)
453 enddo
454!
455! compute the surface variables for pbl height estimation
456! under unstable conditions
457!
458 do i = 1,im
459 if(sfcflg(i).and.sflux(i).gt.0.0)then
460 gamfac = bfac/rhox(i)/wscale(i)
461 hgamt(i) = min(gamfac*hfx(i)/cp,gamcrt)
462 hgamq(i) = min(gamfac*qfx(i),gamcrq)
463 vpert = (hgamt(i)+ep1*thx(i,1)*hgamq(i))/bfac*afac
464 thermal(i) = thermal(i)+max(vpert,0.)*min(za(i,1)/(sfcfrac*hpbl(i)),1.0)
465 thermalli(i)= thermalli(i)+max(vpert,0.)*min(za(i,1)/(sfcfrac*hpbl(i)),1.0)
466 hgamt(i) = max(hgamt(i),0.0)
467 hgamq(i) = max(hgamq(i),0.0)
468 brint = -15.9*ust(i)*ust(i)/wspd(i)*wstar3(i)/(wscale(i)**4.)
469 hgamu(i) = brint*ux(i,1)
470 hgamv(i) = brint*vx(i,1)
471 else
472 pblflg(i) = .false.
473 endif
474 enddo
475!
476! enhance the pbl height by considering the thermal
477!
478 do i = 1,im
479 if(pblflg(i))then
480 kpbl(i) = 1
481 hpbl(i) = zq(i,1)
482 endif
483 enddo
484!
485 do i = 1,im
486 if(pblflg(i))then
487 stable(i) = .false.
488 brup(i) = br(i)
489 brcr(i) = brcr_ub
490 endif
491 enddo
492!
493 do k = 2,klpbl
494 do i = 1,im
495 if(.not.stable(i).and.pblflg(i))then
496 brdn(i) = brup(i)
497 spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.)
498 brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2
499 kpbl(i) = k
500 stable(i) = brup(i).gt.brcr(i)
501 endif
502 enddo
503 enddo
504!
505! enhance pbl by theta-li
506!
507 if (ysu_topdown_pblmix.eq.1)then
508 do i = 1,im
509 kpblold(i) = kpbl(i)
510 definebrup=.false.
511 do k = kpblold(i), km-1
512 spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.)
513 bruptmp = (thlix(i,k)-thermalli(i))*(g*za(i,k)/thlix(i,1))/spdk2
514 stable(i) = bruptmp.ge.brcr(i)
515 if (definebrup) then
516 kpbl(i) = k
517 brup(i) = bruptmp
518 definebrup=.false.
519 endif
520 if (.not.stable(i)) then !overwrite brup brdn values
521 brdn(i)=bruptmp
522 definebrup=.true.
523 pblflg(i)=.true.
524 endif
525 enddo
526 enddo
527 endif
528
529 do i = 1,im
530 if(pblflg(i)) then
531 k = kpbl(i)
532 if(brdn(i).ge.brcr(i))then
533 brint = 0.
534 elseif(brup(i).le.brcr(i))then
535 brint = 1.
536 else
537 brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i))
538 endif
539 hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
540 if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
541 if(kpbl(i).le.1) pblflg(i) = .false.
542 endif
543 enddo
544!
545! stable boundary layer
546!
547 do i = 1,im
548 if((.not.sfcflg(i)).and.hpbl(i).lt.zq(i,2)) then
549 brup(i) = br(i)
550 stable(i) = .false.
551 else
552 stable(i) = .true.
553 endif
554 enddo
555!
556 do i = 1,im
557 if((.not.stable(i)).and.((xland(i)-1.5).ge.0))then
558 wspd10 = u10(i)*u10(i) + v10(i)*v10(i)
559 wspd10 = sqrt(wspd10)
560 ross = wspd10 / (cori*znt(i))
561 brcr_sbro(i) = min(0.16*(1.e-7*ross)**(-0.18),.3)
562 endif
563 enddo
564!
565 do i = 1,im
566 if(.not.stable(i))then
567 if((xland(i)-1.5).ge.0)then
568 brcr(i) = brcr_sbro(i)
569 else
570 brcr(i) = brcr_sb
571 endif
572 endif
573 enddo
574!
575 do k = 2,klpbl
576 do i = 1,im
577 if(.not.stable(i))then
578 brdn(i) = brup(i)
579 spdk2 = max(ux(i,k)**2+vx(i,k)**2,1.)
580 brup(i) = (thvx(i,k)-thermal(i))*(g*za(i,k)/thvx(i,1))/spdk2
581 kpbl(i) = k
582 stable(i) = brup(i).gt.brcr(i)
583 endif
584 enddo
585 enddo
586!
587 do i = 1,im
588 if((.not.sfcflg(i)).and.hpbl(i).lt.zq(i,2)) then
589 k = kpbl(i)
590 if(brdn(i).ge.brcr(i))then
591 brint = 0.
592 elseif(brup(i).le.brcr(i))then
593 brint = 1.
594 else
595 brint = (brcr(i)-brdn(i))/(brup(i)-brdn(i))
596 endif
597 hpbl(i) = za(i,k-1)+brint*(za(i,k)-za(i,k-1))
598 if(hpbl(i).lt.zq(i,2)) kpbl(i) = 1
599 if(kpbl(i).le.1) pblflg(i) = .false.
600 endif
601 enddo
602!
603! estimate the entrainment parameters
604!
605 do i = 1,im
606 cloudflg(i)=.false.
607 if(pblflg(i)) then
608 k = kpbl(i) - 1
609 wm3 = wstar3(i) + 5. * ust3(i)
610 wm2(i) = wm3**h2
611 bfxpbl(i) = -0.15*thvx(i,1)/g*wm3/hpbl(i)
612 dthvx(i) = max(thvx(i,k+1)-thvx(i,k),tmin)
613 we(i) = max(bfxpbl(i)/dthvx(i),-sqrt(wm2(i)))
614 if((qx(i,k,ntcw)+qx(i,k,ntiw)).gt.0.01e-3.and.ysu_topdown_pblmix.eq.1)then
615 if ( kpbl(i) .ge. 2) then
616 cloudflg(i)=.true.
617 templ=thlix(i,k)*(p2di(i,k+1)/100000)**rovcp
618 !rvls is ws at full level
619 rvls=100.*6.112*exp(17.67*(templ-273.16)/(templ-29.65))*(ep2/p2di(i,k+1))
620 temps=templ + ((qx(i,k,1)+qx(i,k,ntcw))-rvls)/(cp/xlv + &
621 ep2*xlv*rvls/(rd*templ**2))
622 rvls=100.*6.112*exp(17.67*(temps-273.15)/(temps-29.65))*(ep2/p2di(i,k+1))
623 rcldb=max((qx(i,k,1)+qx(i,k,ntcw))-rvls,0.)
624 !entrainment efficiency
625 dthvx(i) = (thlix(i,k+2)+thx(i,k+2)*ep1*(qx(i,k+2,1)+qx(i,k+2,ntcw))) &
626 - (thlix(i,k) + thx(i,k) *ep1*(qx(i,k,1) +qx(i,k,ntcw)))
627 dthvx(i) = max(dthvx(i),0.1)
628 tmp1 = xlv/cp * rcldb/(pi2d(i,k)*dthvx(i))
629 ent_eff = 0.2 * 8. * tmp1 +0.2
630
631 radsum=0.
632 do kk = 1,kpbl(i)-1
633 radflux=swh(i,kk)*xmu(i)+hlw(i,kk) !radiative heating rate temp/s
634 radflux=radflux*cp/g*(p2di(i,kk)-p2di(i,kk+1)) ! converts temp/s to W/m^2
635 if (radflux < 0.0 ) radsum=abs(radflux)+radsum
636 enddo
637 radsum=max(radsum,0.0)
638
639 !recompute entrainment from sfc thermals
640 bfx0 = max(max(sflux(i),0.0)-radsum/rhox2(i,k)/cp,0.)
641 bfx0 = max(sflux(i),0.0)
642 wm3 = (govrth(i)*bfx0*hpbl(i))+5. * ust3(i)
643 wm2(i) = wm3**h2
644 bfxpbl(i) = -0.15*thvx(i,1)/g*wm3/hpbl(i)
645 dthvx(i) = max(thvx(i,k+1)-thvx(i,k),tmin)
646 we(i) = max(bfxpbl(i)/dthvx(i),-sqrt(wm2(i)))
647
648 !entrainment from PBL top thermals
649 bfx0 = max(radsum/rhox2(i,k)/cp-max(sflux(i),0.0),0.)
650 bfx0 = max(radsum/rhox2(i,k)/cp,0.)
651 wm3 = (g/thvx(i,k)*bfx0*hpbl(i)) ! this is wstar3(i)
652 wm2(i) = wm2(i)+wm3**h2
653 bfxpbl(i) = - ent_eff * bfx0
654 dthvx(i) = max(thvx(i,k+1)-thvx(i,k),0.1)
655 we(i) = we(i) + max(bfxpbl(i)/dthvx(i),-sqrt(wm3**h2))
656
657 !wstar3_2
658 bfx0 = max(radsum/rhox2(i,k)/cp,0.)
659 wstar3_2(i) = (g/thvx(i,k)*bfx0*hpbl(i))
660 !recompute hgamt
661 wscale(i) = (ust3(i)+phifac*karman*(wstar3(i)+wstar3_2(i))*0.5)**h1
662 wscale(i) = min(wscale(i),ust(i)*aphi16)
663 wscale(i) = max(wscale(i),ust(i)/aphi5)
664 gamfac = bfac/rhox(i)/wscale(i)
665 hgamt(i) = min(gamfac*hfx(i)/cp,gamcrt)
666 hgamq(i) = min(gamfac*qfx(i),gamcrq)
667 gamfac = bfac/rhox2(i,k)/wscale(i)
668 hgamt2(i,k) = min(gamfac*radsum/cp,gamcrt)
669 hgamt(i) = max(hgamt(i),0.0) + max(hgamt2(i,k),0.0)
670 brint = -15.9*ust(i)*ust(i)/wspd(i)*(wstar3(i)+wstar3_2(i))/(wscale(i)**4.)
671 hgamu(i) = brint*ux(i,1)
672 hgamv(i) = brint*vx(i,1)
673 endif
674 endif
675 prpbl(i) = 1.0
676 dthx = max(thx(i,k+1)-thx(i,k),tmin)
677 dqx = min(qx(i,k+1,1)-qx(i,k,1),0.0)
678 hfxpbl(i) = we(i)*dthx
679 qfxpbl(i) = we(i)*dqx
680!
681 dux = ux(i,k+1)-ux(i,k)
682 dvx = vx(i,k+1)-vx(i,k)
683 if(dux.gt.tmin) then
684 ufxpbl(i) = max(prpbl(i)*we(i)*dux,-ust(i)*ust(i))
685 elseif(dux.lt.-tmin) then
686 ufxpbl(i) = min(prpbl(i)*we(i)*dux,ust(i)*ust(i))
687 else
688 ufxpbl(i) = 0.0
689 endif
690 if(dvx.gt.tmin) then
691 vfxpbl(i) = max(prpbl(i)*we(i)*dvx,-ust(i)*ust(i))
692 elseif(dvx.lt.-tmin) then
693 vfxpbl(i) = min(prpbl(i)*we(i)*dvx,ust(i)*ust(i))
694 else
695 vfxpbl(i) = 0.0
696 endif
697 delb = govrth(i)*d3*hpbl(i)
698 delta(i) = min(d1*hpbl(i) + d2*wm2(i)/delb,100.)
699 endif
700 enddo
701!
702 do k = 1,klpbl
703 do i = 1,im
704 if(pblflg(i).and.k.ge.kpbl(i))then
705 entfac(i,k) = ((zq(i,k+1)-hpbl(i))/delta(i))**2.
706 else
707 entfac(i,k) = 1.e30
708 endif
709 enddo
710 enddo
711!
712! compute diffusion coefficients below pbl
713!
714 do k = 1,klpbl
715 do i = 1,im
716 if(k.lt.kpbl(i)) then
717 zfac(i,k) = min(max((1.-(zq(i,k+1)-zl1(i))/(hpbl(i)-zl1(i))),zfmin),1.)
718 zfacent(i,k) = (1.-zfac(i,k))**3.
719 wscalek(i,k) = (ust3(i)+phifac*karman*wstar3(i)*(1.-zfac(i,k)))**h1
720 wscalek2(i,k) = (phifac*karman*wstar3_2(i)*(zfac(i,k)))**h1
721 if(sfcflg(i)) then
722 prfac = conpr
723 prfac2 = 15.9*(wstar3(i)+wstar3_2(i))/ust3(i)/(1.+4.*karman*(wstar3(i)+wstar3_2(i))/ust3(i))
724 prnumfac = -3.*(max(zq(i,k+1)-sfcfrac*hpbl(i),0.))**2./hpbl(i)**2.
725 else
726 prfac = 0.
727 prfac2 = 0.
728 prnumfac = 0.
729 phim8z = 1.+aphi5*zol1(i)*zq(i,k+1)/zl1(i)
730 wscalek(i,k) = ust(i)/phim8z
731 wscalek(i,k) = max(wscalek(i,k),0.001)
732 endif
733 prnum0 = (phih(i)/phim(i)+prfac)
734 prnum0 = max(min(prnum0,prmax),prmin)
735 xkzm(i,k) = wscalek(i,k) *karman* zq(i,k+1) * zfac(i,k)**pfac+ &
736 wscalek2(i,k)*karman*(hpbl(i)-zq(i,k+1))*(1-zfac(i,k))**pfac
737 !Do not include xkzm at kpbl-1 since it changes entrainment
738 if (k.eq.kpbl(i)-1.and.cloudflg(i).and.we(i).lt.0.0) then
739 xkzm(i,k) = 0.0
740 endif
741 prnum = 1. + (prnum0-1.)*exp(prnumfac)
742 xkzq(i,k) = xkzm(i,k)/prnum*zfac(i,k)**(pfac_q-pfac)
743 prnum0 = prnum0/(1.+prfac2*karman*sfcfrac)
744 prnum = 1. + (prnum0-1.)*exp(prnumfac)
745 xkzh(i,k) = xkzm(i,k)/prnum
746 xkzm(i,k) = xkzm(i,k)+xkzom(i,k)
747 xkzh(i,k) = xkzh(i,k)+xkzoh(i,k)
748 xkzq(i,k) = xkzq(i,k)+xkzoh(i,k)
749 xkzm(i,k) = min(xkzm(i,k),xkzmax)
750 xkzh(i,k) = min(xkzh(i,k),xkzmax)
751 xkzq(i,k) = min(xkzq(i,k),xkzmax)
752 endif
753 enddo
754 enddo
755!
756! compute diffusion coefficients over pbl (free atmosphere)
757!
758 do k = 1,km-1
759 do i = 1,im
760 if(k.ge.kpbl(i)) then
761 ss = ((ux(i,k+1)-ux(i,k))*(ux(i,k+1)-ux(i,k)) &
762 +(vx(i,k+1)-vx(i,k))*(vx(i,k+1)-vx(i,k))) &
763 /(dza(i,k+1)*dza(i,k+1))+1.e-9
764 govrthv = g/(0.5*(thvx(i,k+1)+thvx(i,k)))
765 ri = govrthv*(thvx(i,k+1)-thvx(i,k))/(ss*dza(i,k+1))
766 if(imvdif.eq.1.and.ntcw.ge.2.and.ntiw.ge.2)then
767 if((qx(i,k,ntcw)+qx(i,k,ntiw)).gt.0.01e-3.and.(qx(i &
768 ,k+1,ntcw)+qx(i,k+1,ntiw)).gt.0.01e-3)then
769! in cloud
770 qmean = 0.5*(qx(i,k,1)+qx(i,k+1,1))
771 tmean = 0.5*(tx(i,k)+tx(i,k+1))
772 alph = xlv*qmean/rd/tmean
773 chi = xlv*xlv*qmean/cp/rv/tmean/tmean
774 ri = (1.+alph)*(ri-g*g/ss/tmean/cp*((chi-alph)/(1.+chi)))
775 endif
776 endif
777 zk = karman*zq(i,k+1)
778 rlamdz = min(max(0.1*dza(i,k+1),rlam),300.)
779 rlamdz = min(dza(i,k+1),rlamdz)
780 rl2 = (zk*rlamdz/(rlamdz+zk))**2
781 dk = rl2*sqrt(ss)
782 if(ri.lt.0.)then
783! unstable regime
784 ri = max(ri, rimin)
785 sri = sqrt(-ri)
786 xkzm(i,k) = dk*(1+8.*(-ri)/(1+1.746*sri))
787 xkzh(i,k) = dk*(1+8.*(-ri)/(1+1.286*sri))
788 else
789! stable regime
790 xkzh(i,k) = dk/(1+5.*ri)**2
791 prnum = 1.0+2.1*ri
792 prnum = min(prnum,prmax)
793 xkzm(i,k) = xkzh(i,k)*prnum
794 endif
795!
796 xkzm(i,k) = xkzm(i,k)+xkzom(i,k)
797 xkzh(i,k) = xkzh(i,k)+xkzoh(i,k)
798 xkzm(i,k) = min(xkzm(i,k),xkzmax)
799 xkzh(i,k) = min(xkzh(i,k),xkzmax)
800 xkzml(i,k) = xkzm(i,k)
801 xkzhl(i,k) = xkzh(i,k)
802 endif
803 enddo
804 enddo
805!
806! compute tridiagonal matrix elements for heat
807!
808 do k = 1,km
809 do i = 1,im
810 au(i,k) = 0.
811 al(i,k) = 0.
812 ad(i,k) = 0.
813 f1(i,k) = 0.
814 enddo
815 enddo
816!
817 do i = 1,im
818 ad(i,1) = 1.
819 f1(i,1) = thx(i,1)-300.+hfx(i)/cont/del(i,1)*dt2
820 enddo
821!
822 do k = 1,km-1
823 do i = 1,im
824 dtodsd = dt2/del(i,k)
825 dtodsu = dt2/del(i,k+1)
826 dsig = p2d(i,k)-p2d(i,k+1)
827 rdz = 1./dza(i,k+1)
828 tem1 = dsig*xkzh(i,k)*rdz
829 if(pblflg(i).and.k.lt.kpbl(i)) then
830 dsdzt = tem1*(-hgamt(i)/hpbl(i)-hfxpbl(i)*zfacent(i,k)/xkzh(i,k))
831 f1(i,k) = f1(i,k)+dtodsd*dsdzt
832 f1(i,k+1) = thx(i,k+1)-300.-dtodsu*dsdzt
833 elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then
834 xkzh(i,k) = -we(i)*dza(i,kpbl(i))*exp(-entfac(i,k))
835 xkzh(i,k) = sqrt(xkzh(i,k)*xkzhl(i,k))
836 xkzh(i,k) = max(xkzh(i,k),xkzoh(i,k))
837 xkzh(i,k) = min(xkzh(i,k),xkzmax)
838 f1(i,k+1) = thx(i,k+1)-300.
839 else
840 f1(i,k+1) = thx(i,k+1)-300.
841 endif
842 tem1 = dsig*xkzh(i,k)*rdz
843 dsdz2 = tem1*rdz
844 au(i,k) = -dtodsd*dsdz2
845 al(i,k) = -dtodsu*dsdz2
846 ad(i,k) = ad(i,k)-au(i,k)
847 ad(i,k+1) = 1.-al(i,k)
848 enddo
849 enddo
850!
851! copies here to avoid duplicate input args for tridin
852!
853 do k = 1,km
854 do i = 1,im
855 cu(i,k) = au(i,k)
856 r1(i,k) = f1(i,k)
857 enddo
858 enddo
859!
860 call tridin_ysu(al,ad,cu,r1,au,f1,im,km,1)
861!
862! recover tendencies of heat
863!
864 do k = km,1,-1
865 do i = 1,im
866 ttend = (f1(i,k)-thx(i,k)+300.)*rdt*pi2d(i,k)
867 ttnp(i,k) = ttnp(i,k)+ttend
868 dtsfc(i) = dtsfc(i)+ttend*cont*del(i,k)
869 enddo
870 enddo
871 if(lssav .and. ldiag3d .and. .not. flag_for_pbl_generic_tend) then
872 idtend = dtidx(index_of_temperature,index_of_process_pbl)
873 if(idtend>=1) then
874 dtend(:,:,idtend) = dtend(:,:,idtend) + dtstep*(f1-thx+300.)*rdt*pi2d
875 endif
876 endif
877!
878! compute tridiagonal matrix elements for moisture, clouds, and gases
879!
880 do k = 1,km
881 do i = 1,im
882 au(i,k) = 0.
883 al(i,k) = 0.
884 ad(i,k) = 0.
885 enddo
886 enddo
887!
888 do ic = 1,ndiff
889 do i = 1,im
890 do k = 1,km
891 f3(i,k,ic) = 0.
892 enddo
893 enddo
894 enddo
895!
896 do i = 1,im
897 ad(i,1) = 1.
898 f3(i,1,1) = qx(i,1,1)+qfx(i)*g/del(i,1)*dt2
899 enddo
900!
901 if(ndiff.ge.2) then
902 do ic = 2,ndiff
903 do i = 1,im
904 f3(i,1,ic) = qx(i,1,ic)
905 enddo
906 enddo
907 endif
908!
909 do k = 1,km-1
910 do i = 1,im
911 if(k.ge.kpbl(i)) then
912 xkzq(i,k) = xkzh(i,k)
913 endif
914 enddo
915 enddo
916!
917 do k = 1,km-1
918 do i = 1,im
919 dtodsd = dt2/del(i,k)
920 dtodsu = dt2/del(i,k+1)
921 dsig = p2d(i,k)-p2d(i,k+1)
922 rdz = 1./dza(i,k+1)
923 tem1 = dsig*xkzq(i,k)*rdz
924 if(pblflg(i).and.k.lt.kpbl(i)) then
925 dsdzq = tem1*(-qfxpbl(i)*zfacent(i,k)/xkzq(i,k))
926 f3(i,k,1) = f3(i,k,1)+dtodsd*dsdzq
927 f3(i,k+1,1) = qx(i,k+1,1)-dtodsu*dsdzq
928 elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then
929 xkzq(i,k) = -we(i)*dza(i,kpbl(i))*exp(-entfac(i,k))
930 xkzq(i,k) = sqrt(xkzq(i,k)*xkzhl(i,k))
931 xkzq(i,k) = max(xkzq(i,k),xkzoh(i,k))
932 xkzq(i,k) = min(xkzq(i,k),xkzmax)
933 f3(i,k+1,1) = qx(i,k+1,1)
934 else
935 f3(i,k+1,1) = qx(i,k+1,1)
936 endif
937 tem1 = dsig*xkzq(i,k)*rdz
938 dsdz2 = tem1*rdz
939 au(i,k) = -dtodsd*dsdz2
940 al(i,k) = -dtodsu*dsdz2
941 ad(i,k) = ad(i,k)-au(i,k)
942 ad(i,k+1) = 1.-al(i,k)
943 enddo
944 enddo
945!
946 if(ndiff.ge.2) then
947 do ic = 2,ndiff
948 do k = 1,km-1
949 do i = 1,im
950 f3(i,k+1,ic) = qx(i,k+1,ic)
951 enddo
952 enddo
953 enddo
954 endif
955!
956! copies here to avoid duplicate input args for tridin
957!
958 do k = 1,km
959 do i = 1,im
960 cu(i,k) = au(i,k)
961 enddo
962 enddo
963!
964 do ic = 1,ndiff
965 do k = 1,km
966 do i = 1,im
967 r3(i,k,ic) = f3(i,k,ic)
968 enddo
969 enddo
970 enddo
971!
972! solve tridiagonal problem for moisture, clouds, and gases
973!
974 call tridin_ysu(al,ad,cu,r3,au,f3,im,km,ndiff)
975!
976! recover tendencies of heat and moisture
977!
978 do k = km,1,-1
979 do i = 1,im
980 qtend = (f3(i,k,1)-qx(i,k,1))*rdt
981 qtnp(i,k,1) = qtnp(i,k,1)+qtend
982 dqsfc(i) = dqsfc(i)+qtend*conq*del(i,k)
983 enddo
984 enddo
985 if(lssav .and. ldiag3d .and. qdiag3d .and. .not. flag_for_pbl_generic_tend) then
986 idtend = dtidx(ntqv+100,index_of_process_pbl)
987 if(idtend>=1) then
988 dtend(:,:,idtend) = dtend(:,:,idtend) + dtstep*(f3(:,:,1)-qx(:,:,1))*rdt
989 endif
990 endif
991!
992 if(ndiff.ge.2) then
993 do ic = 2,ndiff
994 do k = km,1,-1
995 do i = 1,im
996 qtend = (f3(i,k,ic)-qx(i,k,ic))*rdt
997 qtnp(i,k,ic) = qtnp(i,k,ic)+qtend
998 enddo
999 enddo
1000 enddo
1001 if(lssav .and. ldiag3d .and. ntoz>0 .and. qdiag3d .and. &
1002 & .not. flag_for_pbl_generic_tend) then
1003 idtend = dtidx(100+ntoz,index_of_process_pbl)
1004 if(idtend>=1) then
1005 dtend(:,:,idtend) = dtend(:,:,idtend) + f3(:,:,ntoz)-qx(:,:,ntoz)
1006 endif
1007 endif
1008 endif
1009!
1010! compute tridiagonal matrix elements for momentum
1011!
1012 do i = 1,im
1013 do k = 1,km
1014 au(i,k) = 0.
1015 al(i,k) = 0.
1016 ad(i,k) = 0.
1017 f1(i,k) = 0.
1018 f2(i,k) = 0.
1019 enddo
1020 enddo
1021!
1022 do i = 1,im
1023 ad(i,1) = 1.+ust(i)**2/wspd1(i)*rhox(i)*g/del(i,1)*dt2 &
1024 *(wspd1(i)/wspd(i))**2
1025 f1(i,1) = ux(i,1)+uox(i)*ust(i)**2*g/del(i,1)*dt2/wspd1(i)
1026 f2(i,1) = vx(i,1)+vox(i)*ust(i)**2*g/del(i,1)*dt2/wspd1(i)
1027 enddo
1028!
1029 do k = 1,km-1
1030 do i = 1,im
1031 dtodsd = dt2/del(i,k)
1032 dtodsu = dt2/del(i,k+1)
1033 dsig = p2d(i,k)-p2d(i,k+1)
1034 rdz = 1./dza(i,k+1)
1035 tem1 = dsig*xkzm(i,k)*rdz
1036 if(pblflg(i).and.k.lt.kpbl(i))then
1037 dsdzu = tem1*(-hgamu(i)/hpbl(i)-ufxpbl(i)*zfacent(i,k)/xkzm(i,k))
1038 dsdzv = tem1*(-hgamv(i)/hpbl(i)-vfxpbl(i)*zfacent(i,k)/xkzm(i,k))
1039 f1(i,k) = f1(i,k)+dtodsd*dsdzu
1040 f1(i,k+1) = ux(i,k+1)-dtodsu*dsdzu
1041 f2(i,k) = f2(i,k)+dtodsd*dsdzv
1042 f2(i,k+1) = vx(i,k+1)-dtodsu*dsdzv
1043 elseif(pblflg(i).and.k.ge.kpbl(i).and.entfac(i,k).lt.4.6) then
1044 xkzm(i,k) = prpbl(i)*xkzh(i,k)
1045 xkzm(i,k) = sqrt(xkzm(i,k)*xkzml(i,k))
1046 xkzm(i,k) = max(xkzm(i,k),xkzom(i,k))
1047 xkzm(i,k) = min(xkzm(i,k),xkzmax)
1048 f1(i,k+1) = ux(i,k+1)
1049 f2(i,k+1) = vx(i,k+1)
1050 else
1051 f1(i,k+1) = ux(i,k+1)
1052 f2(i,k+1) = vx(i,k+1)
1053 endif
1054 tem1 = dsig*xkzm(i,k)*rdz
1055 dsdz2 = tem1*rdz
1056 au(i,k) = -dtodsd*dsdz2
1057 al(i,k) = -dtodsu*dsdz2
1058 ad(i,k) = ad(i,k)-au(i,k)
1059 ad(i,k+1) = 1.-al(i,k)
1060 enddo
1061 enddo
1062!
1063! copies here to avoid duplicate input args for tridin
1064!
1065 do k = 1,km
1066 do i = 1,im
1067 cu(i,k) = au(i,k)
1068 r1(i,k) = f1(i,k)
1069 r2(i,k) = f2(i,k)
1070 enddo
1071 enddo
1072!
1073! solve tridiagonal problem for momentum
1074!
1075 call tridi1n(al,ad,cu,r1,r2,au,f1,f2,im,km,1)
1076!
1077! recover tendencies of momentum
1078!
1079 do k = km,1,-1
1080 do i = 1,im
1081 utend = (f1(i,k)-ux(i,k))*rdt
1082 vtend = (f2(i,k)-vx(i,k))*rdt
1083 utnp(i,k) = utnp(i,k)+utend
1084 vtnp(i,k) = vtnp(i,k)+vtend
1085 dusfc(i) = dusfc(i) + utend*conwrc*del(i,k)
1086 dvsfc(i) = dvsfc(i) + vtend*conwrc*del(i,k)
1087 enddo
1088 enddo
1089 if(lssav .and. ldiag3d .and. .not. flag_for_pbl_generic_tend) then
1090 idtend = dtidx(index_of_x_wind,index_of_process_pbl)
1091 if(idtend>=1) then
1092 dtend(:,:,idtend) = dtend(:,:,idtend) + dtstep*(f1-ux)*rdt
1093 endif
1094
1095 idtend = dtidx(index_of_y_wind,index_of_process_pbl)
1096 if(idtend>=1) then
1097 dtend(:,:,idtend) = dtend(:,:,idtend) + dtstep*(f2-vx)*rdt
1098 endif
1099 endif
1100!
1101!---- end of vertical diffusion
1102!
1103 do i = 1,im
1104 kpbl1d(i) = kpbl(i)
1105 enddo
1106!
1107!
1108 end subroutine ysuvdif_run
1109!-------------------------------------------------------------------------------
1110!
1111!-------------------------------------------------------------------------------
1112 subroutine tridi1n(cl,cm,cu,r1,r2,au,f1,f2,im,km,nt)
1113 use machine , only : kind_phys
1114!-------------------------------------------------------------------------------
1115 implicit none
1116!-------------------------------------------------------------------------------
1117!
1118 integer, intent(in ) :: im, km, nt
1119!
1120 real(kind=kind_phys), dimension( im, 2:km+1 ) , &
1121 intent(in ) :: cl
1122!
1123 real(kind=kind_phys), dimension( im, km ) , &
1124 intent(in ) :: cm, &
1125 r1
1126 real(kind=kind_phys), dimension( im, km,nt ) , &
1127 intent(in ) :: r2
1128!
1129 real(kind=kind_phys), dimension( im, km ) , &
1130 intent(inout) :: au, &
1131 cu, &
1132 f1
1133 real(kind=kind_phys), dimension( im, km,nt ) , &
1134 intent(inout) :: f2
1135!
1136 real(kind=kind_phys) :: fk
1137 integer :: i,k,l,n,it
1138!
1139!-------------------------------------------------------------------------------
1140!
1141 l = im
1142 n = km
1143!
1144 do i = 1,l
1145 fk = 1./cm(i,1)
1146 au(i,1) = fk*cu(i,1)
1147 f1(i,1) = fk*r1(i,1)
1148 enddo
1149!
1150 do it = 1,nt
1151 do i = 1,l
1152 fk = 1./cm(i,1)
1153 f2(i,1,it) = fk*r2(i,1,it)
1154 enddo
1155 enddo
1156!
1157 do k = 2,n-1
1158 do i = 1,l
1159 fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1))
1160 au(i,k) = fk*cu(i,k)
1161 f1(i,k) = fk*(r1(i,k)-cl(i,k)*f1(i,k-1))
1162 enddo
1163 enddo
1164!
1165 do it = 1,nt
1166 do k = 2,n-1
1167 do i = 1,l
1168 fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1))
1169 f2(i,k,it) = fk*(r2(i,k,it)-cl(i,k)*f2(i,k-1,it))
1170 enddo
1171 enddo
1172 enddo
1173!
1174 do i = 1,l
1175 fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1))
1176 f1(i,n) = fk*(r1(i,n)-cl(i,n)*f1(i,n-1))
1177 enddo
1178!
1179 do it = 1,nt
1180 do i = 1,l
1181 fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1))
1182 f2(i,n,it) = fk*(r2(i,n,it)-cl(i,n)*f2(i,n-1,it))
1183 enddo
1184 enddo
1185!
1186 do k = n-1,1,-1
1187 do i = 1,l
1188 f1(i,k) = f1(i,k)-au(i,k)*f1(i,k+1)
1189 enddo
1190 enddo
1191!
1192 do it = 1,nt
1193 do k = n-1,1,-1
1194 do i = 1,l
1195 f2(i,k,it) = f2(i,k,it)-au(i,k)*f2(i,k+1,it)
1196 enddo
1197 enddo
1198 enddo
1199!
1200 end subroutine tridi1n
1201!-------------------------------------------------------------------------------
1202!
1203!-------------------------------------------------------------------------------
1204 subroutine tridin_ysu(cl,cm,cu,r2,au,f2,im,km,nt)
1205 use machine , only : kind_phys
1206!-------------------------------------------------------------------------------
1207 implicit none
1208!-------------------------------------------------------------------------------
1209!
1210 integer, intent(in ) :: im, km, nt
1211!
1212 real(kind=kind_phys), dimension( im, 2:km+1 ) , &
1213 intent(in ) :: cl
1214!
1215 real(kind=kind_phys), dimension( im, km ) , &
1216 intent(in ) :: cm
1217 real(kind=kind_phys), dimension( im, km,nt ) , &
1218 intent(in ) :: r2
1219!
1220 real(kind=kind_phys), dimension( im, km ) , &
1221 intent(inout) :: au, &
1222 cu
1223 real(kind=kind_phys), dimension( im, km,nt ) , &
1224 intent(inout) :: f2
1225!
1226 real(kind=kind_phys) :: fk
1227 integer :: i,k,l,n,it
1228!
1229!-------------------------------------------------------------------------------
1230!
1231 l = im
1232 n = km
1233!
1234 do it = 1,nt
1235 do i = 1,l
1236 fk = 1./cm(i,1)
1237 au(i,1) = fk*cu(i,1)
1238 f2(i,1,it) = fk*r2(i,1,it)
1239 enddo
1240 enddo
1241!
1242 do it = 1,nt
1243 do k = 2,n-1
1244 do i = 1,l
1245 fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1))
1246 au(i,k) = fk*cu(i,k)
1247 f2(i,k,it) = fk*(r2(i,k,it)-cl(i,k)*f2(i,k-1,it))
1248 enddo
1249 enddo
1250 enddo
1251!
1252 do it = 1,nt
1253 do i = 1,l
1254 fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1))
1255 f2(i,n,it) = fk*(r2(i,n,it)-cl(i,n)*f2(i,n-1,it))
1256 enddo
1257 enddo
1258!
1259 do it = 1,nt
1260 do k = n-1,1,-1
1261 do i = 1,l
1262 f2(i,k,it) = f2(i,k,it)-au(i,k)*f2(i,k+1,it)
1263 enddo
1264 enddo
1265 enddo
1266!
1267 end subroutine tridin_ysu
1268!-------------------------------------------------------------------------------
1269end module ysuvdif
1270!-------------------------------------------------------------------------------