CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
myjpbl_wrapper.F90
1
3
5
6 USE machine, only: kfpt => kind_phys, &
7 kind_phys
8
9 contains
10
11 subroutine myjpbl_wrapper_init (do_myjpbl,errmsg,errflg)
12
13 logical, intent(in) :: do_myjpbl
14 character(len=*), intent(out) :: errmsg
15 integer, intent(out) :: errflg
16
17 ! Initialize CCPP error handling variables
18 errmsg = ''
19 errflg = 0
20
21 ! Consistency checks
22 if (.not. do_myjpbl) then
23 write(errmsg,fmt='(*(a))') 'Logic error: do_myjpbl=.false.'
24 errflg = 1
25 return
26 end if
27 end subroutine myjpbl_wrapper_init
28
33!###===================================================================
34 SUBROUTINE myjpbl_wrapper_run( &
35 restart,do_myjsfc, &
36 im,levs,dt_phs, &
37 kdt,ntrac,ntke, &
38 ntcw,ntiw,ntrw,ntsw,ntgl, &
39 ugrs, vgrs, tgrs, qgrs, &
40 prsl, prsi, phii, hprime1, &
41 prsik_1, prslk_1, prslki, tsfc, qsfc, &
42 phy_myj_qsfc, phy_myj_thz0, phy_myj_qz0, &
43 phy_myj_uz0, phy_myj_vz0, phy_myj_z0base, &
44 phy_myj_akhs, phy_myj_akms, &
45 phy_myj_chkqlm, phy_myj_elflx, &
46 phy_myj_a1u, phy_myj_a1t, phy_myj_a1q, &
47 pblh, kpbl, kinver, slmsk, &
48 garea, ustar, cm, ch, wind, &
49 snowd, zorl, evap, hflx, &
50 dudt, dvdt, dtdt, dqdt, &
51 dusfc,dvsfc,dtsfc,dqsfc, &
52 dkt,xkzm_m, xkzm_h,xkzm_s, gamt,gamq, &
53 con_cp,con_g,con_rd, &
54 me, lprnt, gen_tend, ldiag3d, dtend, dtidx, &
55 index_of_temperature, index_of_x_wind, &
56 index_of_y_wind, index_of_process_pbl, &
57 ntqv, errmsg, errflg )
58
59!
60
61 use module_bl_myjpbl, only: myjpbl_init,myjpbl
62
63!-------------------------------------------------------------------
64 implicit none
65
66! integer,parameter:: &
67! klog=4 & ! logical variables
68! ,kint=4 & ! integer variables
69! !,kfpt=4 & ! floating point variables
70! ,kfpt=8 & ! floating point variables
71! ,kdbl=8 ! double precision
72
73!-------------------------------------------------------------------
74! --- constant parameters:
75!For reference
76! real , parameter :: karman = 0.4
77! real , parameter :: g = 9.81
78! real , parameter :: r_d = 287.
79! real , parameter :: cp = 7.*r_d/2.
80!
81! real, parameter :: g = 9.81, r_d=287., cp= 7.*r_d/2.
82! real, parameter :: rd=r_d, rk=cp/rd
83! real, parameter :: elwv=2.501e6, eliv=2.834e6
84! real, parameter :: reliw=eliv/elwv,
85 real, parameter :: xkgdx=25000.,xkzinv=0.15
86
87! real, parameter :: g_inv=1./con_g, cappa=con_rd/con_cp
88
89 character(len=*), intent(out) :: errmsg
90 integer, intent(out) :: errflg
91
92 real(kind=kind_phys), intent(inout), optional :: dtend(:,:,:)
93 integer, intent(in) :: dtidx(:,:)
94 integer, intent(in) :: index_of_temperature, index_of_x_wind, &
95 index_of_y_wind, index_of_process_pbl, ntqv
96
97!MYJ-1D
98 integer,intent(in) :: im, levs
99 integer,intent(in) :: kdt, me
100 integer,intent(in) :: ntrac,ntke,ntcw,ntiw,ntrw,ntsw,ntgl
101 logical,intent(in) :: restart,do_myjsfc,lprnt,ldiag3d,gen_tend
102 real(kind=kind_phys),intent(in) :: con_cp, con_g, con_rd
103 real(kind=kind_phys),intent(in) :: dt_phs, xkzm_m, xkzm_h, xkzm_s
104
105!MYJ-2D
106 real(kind=kind_phys),dimension(:),intent(in) :: &
107 prsik_1, prslk_1, prslki, slmsk, garea, &
108 snowd, evap, hflx, cm, ch, wind, hprime1
109 real(kind=kind_phys),dimension(:),intent(inout) :: &
110 zorl, ustar, tsfc, qsfc
111 real(kind=kind_phys),dimension(:),intent(inout),optional :: &
112 phy_myj_thz0, phy_myj_z0base, phy_myj_chkqlm, &
113 phy_myj_akhs, phy_myj_akms, phy_myj_qz0, &
114 phy_myj_qsfc, phy_myj_elflx, phy_myj_a1u, &
115 phy_myj_a1t, phy_myj_a1q, phy_myj_uz0, &
116 phy_myj_vz0
117 real(kind=kind_phys),dimension(:),intent(out) :: &
118 pblh,dusfc,dvsfc,dtsfc,dqsfc,gamt,gamq
119 integer,dimension(:),intent(out) :: kpbl
120 integer,dimension(:),intent(in) :: kinver
121
122!MYJ-3D
123 real(kind=kind_phys),dimension(:,:),intent(in) :: &
124 phii, prsi
125 real(kind=kind_phys),dimension(:,:),intent(in) :: &
126 ugrs, vgrs, tgrs, prsl
127! real(kind=kind_phys),dimension(:,:),intent(inout) :: &
128! dudt, dvdt, dtdt, dkt
129 real(kind=kind_phys),dimension(:,:),intent(inout) :: &
130 dudt, dvdt, dtdt
131 real(kind=kind_phys),dimension(:,:),intent(out) :: &
132 dkt
133
134!MYJ-4D
135 real(kind=kind_phys),dimension(:,:,:),intent(inout) :: &
136 qgrs,dqdt
137
138!LOCAL
139 integer :: ntsd, k, k1, i, kx1
140 integer :: i_min, i_max, k_min, k_max
141
142 logical :: lprnt1,lprnt2
143 integer :: ict, ide, lm, me1
144 real(kind=kfpt) :: dt_myj, tem, tem1, tem2, ptem
145 integer,dimension(im) :: kpbl_myj
146 real(kind=kfpt),dimension(1:levs-1):: epsl
147 real(kind=kfpt),dimension(1:levs):: epsq2
148 real(kind=kfpt),dimension(im) :: &
149 xland, sice, snowd1, ht, stdh, tsk, &
150 ustar1,z0,pblh_myj, &
151 elflx,mixht,ct
152 real(kind=kfpt), dimension(im,levs) :: &
153 u_myj, v_myj, t_myj, q_myj, th_myj, &
154 cw, dz_myj, pmid, q2, exner, del
155 real(kind=kfpt), dimension(im,levs+1) :: pint
156 real(kind=kfpt),dimension(im,levs) :: &
157 rublten,rvblten,rthblten,rqvblten,rqcblten
158 real(kind=kfpt),dimension(im,levs) :: el_myj
159 real(kind=kfpt),dimension(im) :: &
160 dusfc1,dvsfc1,dtsfc1,dqsfc1
161 real(kind=kfpt),dimension(im) :: thlm,qlm
162 real(kind=kfpt),dimension(im,13) :: phy_f2d_myj
163 real(kind=kfpt), dimension(im,levs) :: xcofh &
164 ,xkzo,xkzmo
165 real(kind=kind_phys) :: g, r_d, g_inv, cappa
166 real(kind=kind_phys) :: thz0, qz0, a1u, a1t, a1q
167 real(kind=kind_phys) :: z0m, aa1u, aa1t, z1uov, z1tox
168 real(kind=kind_phys) :: tmax,tmin,t_myj1
169 real(kind=kind_phys),dimension(im) :: &
170 thsfc,sfcz,tsfc1, &
171 sm,work3,wind1,work4 &
172 ,rho,qfc1,gdx,xkzm_hx,xkzm_mx,tx1, tx2
173! real(kind=kind_phys), dimension(im,levs,ntrac) :: &
174! & qgrs_myj
175 real(kind=kind_phys),dimension(im,levs) :: dkt2
176 integer :: uidx, vidx, tidx, qidx
177
178 ! Initialize CCPP error handling variables
179 errmsg = ''
180 errflg = 0
181
182
183! if (lprnt) then
184! write(0,*)"=============================================="
185! write(0,*)"in myj wrapper..."
186! endif
187
188 ntsd = kdt - 1
189
190 lprnt1=.false.
191 lprnt2=.false.
192
193 if (lprnt1) then
194 if(me.eq.0)print*,'ntsd=', ntsd
195 end if
196
197!prep MYJ-only variables
198
199 r_d = con_rd
200 g = con_g
201 g_inv = 1./con_g
202 cappa = con_rd/con_cp
203
204 do i=1,im
205 work3(i)=prsik_1(i) / prslk_1(i)
206 sice(i)=slmsk(i)*0.5
207 if(sice(i) < 0.7)sice(i)=0
208 sm(i)=1.; if(slmsk(i) > 0.5 ) sm(i)=0.
209 z0(i)=zorl(i)*0.01
210 xland(i)=sm(i)+1.
211 sfcz(i)=phii(i,1)*g_inv
212 work4(i)=(1.e5/prsi(i,1))**cappa
213 thsfc(i)=tsfc(i)*work4(i) ! thsfc
214 enddo
215
216 do k=1,levs
217 k1=levs+1-k
218 do i=1,im
219 u_myj(i,k)=ugrs(i,k1)
220 v_myj(i,k)=vgrs(i,k1)
221 t_myj(i,k)=tgrs(i,k1)
222 q_myj(i,k)=qgrs(i,k1,1)
223 cw(i,k) =qgrs(i,k1,ntcw)
224! if(ntrw.gt.0)cw(i,k) = cw(i,k) + qgrs(i,k1,ntrw)
225! if(ntiw.gt.0)cw(i,k) = cw(i,k) + qgrs(i,k1,ntiw)
226! if(ntsw.gt.0)cw(i,k) = cw(i,k) + qgrs(i,k1,ntsw)
227! if(ntgl.gt.0)cw(i,k) = cw(i,k) + qgrs(i,k1,ntgl)
228 if(ntke.gt.0)then
229 q2(i,k) =max(0.02,qgrs(i,k1,ntke)*2.)
230 else
231 q2(i,k) =0.02
232 end if
233! fmid(i,k) =phil(i,k1)
234 pmid(i,k) =prsl(i,k1)
235 exner(i,k)=(prsl(i,k1)*1.e-5)**cappa
236 th_myj(i,k)=tgrs(i,k1)/exner(i,k)
237 end do
238 end do
239 do k=1,levs+1
240 k1=levs+2-k
241 do i=1,im
242 pint(i,k) =prsi(i,k1)
243 end do
244 end do
245
246 do i=1,im
247 gdx(i) = sqrt(garea(i))
248 enddo
249
250 do i=1,im
251 kx1 = 1
252 tx1(i) = 1.0 / prsi(i,1)
253 tx2(i) = tx1(i)
254 if(gdx(i) >= xkgdx) then
255 xkzm_hx(i) = xkzm_h
256 xkzm_mx(i) = xkzm_m
257 else
258 tem = 1. / (xkgdx - 5.)
259 tem1 = (xkzm_h - 0.01) * tem
260 tem2 = (xkzm_m - 0.01) * tem
261 ptem = gdx(i) - 5.
262 xkzm_hx(i) = 0.01 + tem1 * ptem
263 xkzm_mx(i) = 0.01 + tem2 * ptem
264 endif
265 enddo
266 xkzo = 0.0
267 xkzmo = 0.0
268 do k = 1,levs-1
269 do i=1,im
270 if (k < kinver(i)) then
271! vertical background diffusivity
272 ptem = prsi(i,k+1) * tx1(i)
273 tem1 = 1.0 - ptem
274 tem1 = tem1 * tem1 * 10.0
275 xkzo(i,k) = xkzm_hx(i) * min(1.0, exp(-tem1))
276 xkzo(i,k) = min(xkzo(i,k),xkzinv)
277! vertical background diffusivity for momentum
278 if (ptem >= xkzm_s) then
279 xkzmo(i,k) = xkzm_mx(i)
280 kx1 = k + 1
281 else
282 if (k == kx1 .and. k > 1) tx2(i) = 1.0 / prsi(i,k)
283 tem1 = 1.0 - prsi(i,k+1) * tx2(i)
284 tem1 = tem1 * tem1 * 5.0
285 xkzmo(i,k) = xkzm_mx(i) * min(1.0, exp(-tem1))
286 xkzmo(i,k) = min(xkzmo(i,k),xkzinv)
287 endif
288 endif
289 enddo
290 enddo
291
292! change vertical coordinate
293 do k=1,levs
294 k1=levs+1-k
295 do i=1,im
296 xcofh(i,k1)=xkzo(i,k) ! temp use xcofh
297 el_myj(i,k1)=xkzmo(i,k) ! temp use EL_MYJ
298 end do
299 end do
300
301 do k=1,levs
302 do i=1,im
303 xkzo(i,k)=xcofh(i,k)
304 xkzmo(i,k)=el_myj(i,k)
305 end do
306 end do
307
308 do k=1,levs-1
309 epsq2(k)=0.02
310 epsl(k)=sqrt(epsq2(k)*0.5)
311! if (xkzo(i,k) .gt. 0.01) then
312! epsl(k)=1.0
313! end if
314 end do
315 epsq2(levs)=epsq2(levs-1)
316
317 do k = 1, levs
318 k1 = levs-k+1
319 do i = 1, im
320 del(i,k) = prsi(i,k1) - prsi(i,k1+1)
321 dz_myj(i,k) = (phii(i,k1+1)-phii(i,k1)) * g_inv
322 enddo
323 enddo
324
325 do i = 1, im
326 wind1(i)=max(wind(i),1.0)
327 end do
328
329 if(.not.do_myjsfc)then
330 do i=1,im
331 if(sm(i).gt.0.5.and.sice(i).le.0.5) then
332 z0m=max(0.018*g_inv*ustar(i)*ustar(i),1.59e-5)
333 z1uov=0.35*30.*sqrt(sqrt(z0m*ustar(i)/1.5e-5))/ustar(i)
334 aa1u=cm(i)*wind1(i)*z1uov
335 a1u=aa1u/(1.-aa1u)
336 z1tox=0.84*z1uov
337 aa1t=ch(i)*wind1(i)*z1tox
338 a1t=aa1t/(1.-aa1t)
339!
340! a1u=0.3
341! a1t=0.25
342!
343 a1q=a1t
344 else
345 z0m=zorl(i)*0.01
346 a1u=0.
347 a1t=0.
348 a1q=0.
349 end if
350 phy_myj_a1u(i) = a1u
351 phy_myj_a1t(i) = a1t
352 phy_myj_a1q(i) = a1q
353 phy_myj_akhs(i) = ch(i)*wind1(i)*(1.+a1t)
354 phy_myj_akms(i) = cm(i)*wind1(i)*(1.+a1u)
355 phy_myj_uz0(i) = u_myj(i,levs)*a1u/(1.+a1u)
356 phy_myj_vz0(i) = v_myj(i,levs)*a1u/(1.+a1u)
357 phy_myj_z0base(i)= z0m
358
359 if(ntsd.eq.0)then
360! if(sm(i).gt.0.5)then
361 qz0=max(evap(i)/phy_myj_akhs(i)+q_myj(i,levs),1.e-9)
362 thz0=hflx(i)/phy_myj_akhs(i)+th_myj(i,levs)
363! else
364! if(sice(i).gt.0.5)then
365! qsfc(i)=qss_ice(i)
366! else
367! qsfc(i)=qss_land(i)
368! end if
369! endif
370 phy_myj_thz0(i) = thz0
371 phy_myj_qz0(i) = qz0
372 end if
373 if(cw(i,levs).gt.1.e-9)then
374 phy_myj_chkqlm(i)= 0.
375 else
376 phy_myj_chkqlm(i)= 1.
377 end if
378 end do
379 end if
380
381 if(do_myjsfc)then
382 do i=1,im
383 phy_myj_akhs(i)=phy_myj_akhs(i)*wind1(i)/wind(i)
384 phy_myj_akms(i)=phy_myj_akms(i)*wind1(i)/wind(i)
385 end do
386 end if
387
388! update qsfc, thz0, qz0 and elflx after Land/Ocean model.
389 do i=1,im
390 phy_myj_elflx(i) = evap(i)
391 qsfc(i)=max(evap(i)/(ch(i)*wind1(i))+q_myj(i,levs),1.e-9)
392 tsfc1(i)=(hflx(i)/(ch(i)*wind1(i))+th_myj(i,levs))/work4(i)
393 phy_myj_qsfc(i) = qsfc(i)
394 thz0 = phy_myj_thz0(i)
395 thlm(i)=th_myj(i,levs)
396 qlm(i)=q_myj(i,levs)
397! a1t=phy_myj_a1t(i)
398! thsfc(i)=hflx(i)/phy_myj_akhs(i)+th_myj(i,levs)
399! phy_myj_thz0(i)=((a1t*thlm(i)+thsfc(i))/(a1t+1.)+thz0)*0.5 ! thz0
400 phy_myj_thz0(i)=0.5*(thz0+ &
401 hflx(i)/phy_myj_akhs(i)+th_myj(i,levs))
402! a1q=phy_myj_a1q(i)
403 qz0=phy_myj_qz0(i)
404! phy_myj_qz0(i) = ((a1q*q_myj(i,levs)+qsfc(i))/(a1q+1.)+qz0)*0.5
405 phy_myj_qz0(i) = 0.5*(qz0+ &
406 max(evap(i)/phy_myj_akhs(i)+q_myj(i,levs),1.e-9))
407 enddo
408
409 rthblten = 0.
410 rqvblten = 0.
411 rqcblten= 0.
412 rublten = 0.
413 rvblten = 0.
414! rtrblten= 0.
415 xcofh=0.
416 kpbl(:)=levs-1
417 ict=1 ! no longer used
418
419 if (lprnt1) then
420
421 if (me.eq.0.and.ntsd.lt.2)then
422 print*,'Qingfu test starts PBL'
423 print*,'ntsd,me,im,levs,ict=',ntsd,me,im,levs,ict
424 print*,'dt_phs,sfcz,dz_myj=',dt_phs,sfcz(1),dz_myj(1,5)
425 print*,'pmid,pint,th_myj=',pmid(1,5),pint(1,5),th_myj(1,5)
426 print*,'t_myj,exner,q_myj=',t_myj(1,5),exner(1,5),q_myj(1,5)
427 print*,'cw,u_myj,v_myj=',cw(1,5),u_myj(1,5),v_myj(1,5)
428 print*,'tsfc,xland,sice,snowd=',tsfc(1),xland(1),sice(1),snowd(1)
429 print*,'ustar,z0,pblh,kpbl=',ustar(1),z0(1),pblh(1),kpbl(1)
430 print*,'q2,xcofh=',q2(1,5),xcofh(1,5)
431! print*,'Tbd%phy_f2d_myj(1,1-5)=',(Tbd%phy_f2d_myj(1,i),i=1,5)
432! print*,'Tbd%phy_f2d_myj(1,6-10)=',(Tbd%phy_f2d_myj(1,i),i=6,10)
433! print*,'Tbd%phy_f2d_myj(1,11-13)=',(Tbd%phy_f2d_myj(1,i),i=11,13)
434 print*,'thlm,thsfc=',thlm(i),thsfc(i)
435 end if
436
437 do k=1,levs
438 do i=1,im
439 if(t_myj(i,k).gt.390..or.t_myj(i,k).lt.110.)then
440 print*,'Qingfu test starts PBL',i,k,t_myj(i,k)
441 print*,'ntsd,me,im,levs,ict=',ntsd,me,im,levs,ict
442 print*,'dt_phs,sfcz,dz_myj=',dt_phs,sfcz(i),dz_myj(i,k)
443 print*,'pmid,pint,th_myj=',pmid(i,k),pint(i,k),th_myj(i,k)
444 print*,'t_myj,exner,q_myj=',t_myj(i,k),exner(i,k),q_myj(i,k)
445 print*,'cw,u_myj,v_myj=',cw(i,k),u_myj(i,k),v_myj(i,k)
446 print*,'tsfc,xland,sice,snowd=',tsfc(i),xland(i),sice(i),snowd(i)
447 print*,'ustar,z0,pblh,kpbl=',ustar(i),z0(i),pblh(i),kpbl(i)
448 print*,'q2,xcofh=',q2(i,k),xcofh(i,k)
449 end if
450 end do
451 end do
452
453 tmax=-1.e-5
454 tmin=1.e5
455 do k=1,levs
456 k1=levs+1-k
457 do i=1,im
458 if(tmax.lt.t_myj(i,k1))then
459 tmax=t_myj(i,k1)
460 i_max=i
461 k_max=k
462 end if
463 if(tmin.gt.t_myj(i,k1))then
464 tmin=t_myj(i,k1)
465 i_min=i
466 k_min=k
467 end if
468 end do
469 end do
470! print*,'before i_min,k_min,i_max,k_max=',i_min,k_min,i_max,k_max
471! print*,'ntsd,me,tmin,tmax=',ntsd,me,tmin,tmax
472! if(me.eq.me1.and.tmin.lt.113.6.or.tmax.gt.350.)then
473! i=i_max
474! print*,'before bad bad tmin,tmax=',tmin,tmax,i_min,k_min,i_max,k_max
475! print*,'ntsd,me,tmin,tmax=',ntsd,me,tmin,tmax
476! end if
477
478 end if
479
480 ct=0.
481 ide=im
482 lm=levs
483 dt_myj=dt_phs
484 do i=1,im
485 ustar1(i)=ustar(i)
486 ht(i)=phii(i,1)*g_inv
487 stdh(i)=hprime1(i)
488 tsk(i)=tsfc(i)
489 snowd1(i)=snowd(i)
490 phy_f2d_myj(i,1) = phy_myj_qsfc(i)
491 phy_f2d_myj(i,2) = phy_myj_thz0(i)
492 phy_f2d_myj(i,3) = phy_myj_qz0(i)
493 phy_f2d_myj(i,4) = phy_myj_uz0(i)
494 phy_f2d_myj(i,5) = phy_myj_vz0(i)
495 phy_f2d_myj(i,6) = phy_myj_z0base(i)
496 phy_f2d_myj(i,7) = phy_myj_akhs(i)
497 phy_f2d_myj(i,8) = phy_myj_akms(i)
498 phy_f2d_myj(i,9) = phy_myj_chkqlm(i)
499 phy_f2d_myj(i,10) = phy_myj_elflx(i)
500 phy_f2d_myj(i,11) = phy_myj_a1u(i)
501 phy_f2d_myj(i,12) = phy_myj_a1t(i)
502 phy_f2d_myj(i,13) = phy_myj_a1q(i)
503! do k=1,13
504! phy_f2d_myj(i,k)=Tbd%phy_f2d_myj(i,k)
505! end do
506 end do
507
508! do i = 1, im
509! rho(i) = prsl(i,1)/(r_d*tgrs(i,1) &
510! *(0.608*qgrs(i,1,1)+1.-qgrs(i,1,ntcw)))
511! if(sm(i).lt.0.5)then
512! qfc1(i)=elwv*rho(i)
513! if(snowd(i).gt.0..or.sice(i).gt.0.5)then
514! qfc1(i)=qfc1(i)*reliw
515! end if
516! else
517! qfc1(i)=elwv*rho(i)
518! end if
519! phy_f2d_myj(i,10)=qfc1(i)*phy_f2d_myj(i,10) ! convert units
520! end do
521
522 if(ntsd.eq.0.or.restart)then
523 if(.not.restart) xcofh=0.
524 call myjpbl_init( &
525 1,ide,1,1,lm, &
526 1,ide,1,1, &
527 1,ide,1,1)
528 end if
529
530 call myjpbl(ntsd,me,dt_myj,epsl,epsq2,ht,stdh,dz_myj,del &
531 ,pmid,pint,th_myj,t_myj,exner,q_myj,cw,u_myj,v_myj &
532 ,tsk,phy_f2d_myj(1:im,1),phy_f2d_myj(1:im,9) &
533 ,phy_f2d_myj(1:im,2),phy_f2d_myj(1:im,3) &
534 ,phy_f2d_myj(1:im,4),phy_f2d_myj(1:im,5) &
535 ,xland,sice,snowd1 &
536 ,q2,xcofh,ustar1,z0,el_myj,pblh_myj,kpbl_myj,ct &
537 ,phy_f2d_myj(1:im,7),phy_f2d_myj(1:im,8) &
538 ,phy_f2d_myj(1:im,10),mixht,thlm,qlm &
539 ,rublten,rvblten,rthblten,rqvblten,rqcblten &
540 ,dusfc1,dvsfc1,dtsfc1,dqsfc1,xkzo,xkzmo,ict &
541 ,1,ide,1,1 &
542 ,1,ide,1,1 &
543 ,1,ide,1,1,lm)
544
545 do i=1,im
546 zorl(i)=z0(i)*100.
547 dusfc(i)=dusfc1(i)
548 dvsfc(i)=dvsfc1(i)
549 dtsfc(i)=dtsfc1(i)
550 dqsfc(i)=dqsfc1(i)
551 pblh(i)=pblh_myj(i)
552 kpbl(i)=levs-kpbl_myj(i)
553! ustar(i)=ustar1(i)
554 phy_myj_qsfc(i) = phy_f2d_myj(i,1)
555 phy_myj_thz0(i) = phy_f2d_myj(i,2)
556 phy_myj_qz0(i) = phy_f2d_myj(i,3)
557 phy_myj_uz0(i) = phy_f2d_myj(i,4)
558 phy_myj_vz0(i) = phy_f2d_myj(i,5)
559 phy_myj_z0base(i) = phy_f2d_myj(i,6)
560 phy_myj_akhs(i) = phy_f2d_myj(i,7)
561 phy_myj_akms(i) = phy_f2d_myj(i,8)
562 phy_myj_chkqlm(i) = phy_f2d_myj(i,9)
563 phy_myj_elflx(i) = phy_f2d_myj(i,10)
564 phy_myj_a1u(i) = phy_f2d_myj(i,11)
565 phy_myj_a1t(i) = phy_f2d_myj(i,12)
566 phy_myj_a1q(i) = phy_f2d_myj(i,13)
567! do k=1,13
568! Tbd%phy_f2d_myj(i,k)=phy_f2d_myj(i,k)
569! end do
570 end do
571
572 dkt=0.
573 do k=1,levs
574 k1=levs-k+1
575 do i=1,im
576! dkt(i,k)=max(xcofh(i,k1),xkzo(i,k))
577 dkt(i,k)=xcofh(i,k1)
578 end do
579 end do
580 if(ntke.gt.0)then
581 do k=1,levs
582 k1=levs+1-k
583 qgrs(:,k,ntke)=q2(:,k1)*0.5
584 end do
585 end if
586 gamt=0.
587 gamq=0.
588
589 do k=1,levs
590 k1=levs+1-k
591 do i=1,im
592 dudt(i,k)=dudt(i,k)+rublten(i,k1)
593 dvdt(i,k)=dvdt(i,k)+rvblten(i,k1)
594 dtdt(i,k)=dtdt(i,k)+rthblten(i,k1)*exner(i,k1)
595 dqdt(i,k,1)=dqdt(i,k,1)+rqvblten(i,k1)
596 dqdt(i,k,ntcw)=dqdt(i,k,ntcw)+rqcblten(i,k1)
597 end do
598 end do
599 if (ldiag3d .and. .not. gen_tend) then
600 uidx = dtidx(index_of_x_wind,index_of_process_pbl)
601 vidx = dtidx(index_of_y_wind,index_of_process_pbl)
602 tidx = dtidx(index_of_temperature,index_of_process_pbl)
603 qidx = dtidx(ntqv+100,index_of_process_pbl)
604 ! NOTE: The code that was here before was wrong. It replaced the
605 ! cumulative value with the instantaneous value.
606 do k=1,levs
607 k1=levs+1-k
608 if(uidx>=1) dtend(:,k,uidx)=dtend(:,k,uidx)+rublten(:,k1)*dt_phs
609 if(vidx>=1) dtend(:,k,vidx)=dtend(:,k,vidx)+rvblten(:,k1)*dt_phs
610 if(tidx>=1) dtend(:,k,tidx)=dtend(:,k,tidx)+rthblten(:,k1)*exner(:,k1)*dt_phs
611 if(qidx>=1) dtend(:,k,qidx)=dtend(:,k,qidx)+rqvblten(:,k1)*dt_phs
612 end do
613 end if
614
615 if (lprnt1) then
616
617 do i=1,im
618 if(tsfc(i).gt.350.)then
619 print*,'21tsfc,tsfc1,hflx=',tsfc(i),tsfc1(i),hflx(i)
620 print*,'21qsfc,evap=',qsfc(i),evap(i)
621 end if
622 end do
623 tmax=-1.e-5
624 tmin=1.e5
625 do k=1,levs
626 k1=levs+1-k
627 do i=1,im
628! t_myj1=t_myj(i,k1)+rthblten(i,k1)*exner(i,k1)*dt_phs
629 t_myj1=t_myj(i,k1)+dtdt(i,k)*dt_phs
630 if(tmax.lt.t_myj1)then
631 tmax=t_myj1
632 i_max=i
633 k_max=k
634 me1=me
635 end if
636 if(tmin.gt.t_myj1)then
637 tmin=t_myj1
638 i_min=i
639 k_min=k
640 end if
641 end do
642 end do
643! print*,'2after i_min,k_min,i_max,k_max=',i_min,k_min,i_max,k_max
644! print*,'ntsd,me,tmin,tmax=',ntsd,me,tmin,tmax
645
646 if(me.eq.me1.and.tmin.lt.113.6.or.tmax.gt.350.)then
647 i=i_max
648 print*,'bad bad tmin,tmax=',tmin,tmax,i_min,k_min,i_max,k_max
649
650 do k=1,levs
651 print*,'delt,t_myj=',k,dtdt(i,k)*dt_phs,tgrs(i,k)
652 end do
653
654 print*,'ide,levs,ntsd=',ide,lm,ntsd,dt_myj
655 print*,'epsl,epsq2,ht,stdh,xland,sice,snowd1=', &
656 epsl(i),epsq2(i),ht(i),stdh(i),xland(i),sice(i),snowd1(i)
657 print*,'phy_f2d_myj=', &
658 (phy_f2d_myj(i,k),k=1,13)
659 print*,'tsk(i),ustar1,z0,pblh_myj,kpbl_myj=', &
660 tsk(i),ustar1(i),z0(i),pblh_myj(i),kpbl_myj(i)
661 print*,'mixht=',mixht(i)
662 do k=1,levs
663 print*,'u,v,t=',k,u_myj(i,k),v_myj(i,k), &
664 t_myj(i,k)
665 end do
666 do k=1,levs
667 print*,'q,th,dz_myj=',k,q_myj(i,k),th_myj(i,k),dz_myj(i,k)
668 end do
669 do k=1,levs
670 print*,'del,pmid,pint,=',k,del(i,k), &
671 pmid(i,k),pint(i,k+1)
672 end do
673 do k=1,levs
674 print*,'exner,cw,q2=',k,exner(i,k),cw(i,k), &
675 q2(i,k)
676 end do
677 do k=1,levs
678 print*,'xcofh,el_myj,dkt=',k,xcofh(i,k),el_myj(i,k),dkt(i,k)
679 end do
680 end if
681
682 end if ! lprnt1
683
684 if (lprnt2) then
685
686 tmax=-1.e-5
687 tmin=1.e5
688 do k=1,levs
689 k1=levs+1-k
690 do i=1,im
691! t_myj1=t_myj(i,k1)+rthblten(i,k1)*exner(i,k1)*dt_phs
692 t_myj1=t_myj(i,k1)+dtdt(i,k)*dt_phs
693 if(tmax.lt.t_myj1)then
694 tmax=t_myj1
695 i_max=i
696 k_max=k
697 me1=me
698 end if
699 if(tmin.gt.t_myj1)then
700 tmin=t_myj1
701 i_min=i
702 k_min=k
703 end if
704 end do
705 end do
706 print*,'2after me i_min,k_min,i_max,k_max=',me,i_min,k_min,i_max,k_max
707 print*,'ntsd,tmin,tmax=',ntsd,tmin,tmax
708 print*,'dtdt(i,j)=',dtdt(i_max,k_max)*dt_phs,t_myj(i_max,k_max)
709
710 tmax=-1.e-5
711 tmin=1.e5
712 do k=1,levs
713 k1=levs+1-k
714 do i=1,im
715! t_myj1=t_myj(i,k1)+rthblten(i,k1)*exner(i,k1)*dt_phs
716 t_myj1=ugrs(i,k)+dudt(i,k)*dt_phs
717! t_myj1=dudt(i,k)*dt_phs
718 if(tmax.lt.t_myj1)then
719 tmax=t_myj1
720 i_max=i
721 k_max=k
722 end if
723 if(tmin.gt.t_myj1)then
724 tmin=t_myj1
725 i_min=i
726 k_min=k
727 end if
728 end do
729 end do
730 print*,'3after i_min,k_min,i_max,k_max=',i_min,k_min,i_max,k_max
731 print*,'ntsd,me,tmin,tmax=',ntsd,me,tmin,tmax
732 print*,'dudt(i,k)=',dudt(i_max,k_max)*dt_phs,ugrs(i_max,k_max)
733
734 if(tmax.gt.200.or.tmin.lt.-200)then
735 print*,'bad,bad,bad=',dudt(i_max,k_max)*dt_phs,ugrs(i_max,k_max)
736 do k=1,levs
737 print*,'k,dudt*dt_phs,ugrs=',k,dudt(i_max,k)*dt_phs,ugrs(i_max,k)
738 end do
739 end if
740
741 tmax=-1.e-5
742 tmin=1.e5
743 do k=1,levs
744 k1=levs+1-k
745 do i=1,im
746! t_myj1=t_myj(i,k1)+rthblten(i,k1)*exner(i,k1)*dt_phs
747 t_myj1=vgrs(i,k)+dvdt(i,k)*dt_phs
748! t_myj1=dvdt(i,k)*dt_phs
749 if(tmax.lt.t_myj1)then
750 tmax=t_myj1
751 i_max=i
752 k_max=k
753 end if
754 if(tmin.gt.t_myj1)then
755 tmin=t_myj1
756 i_min=i
757 k_min=k
758 end if
759 end do
760 end do
761 print*,'4after i_min,k_min,i_max,k_max=',i_min,k_min,i_max,k_max
762 print*,'ntsd,me,tmin,tmax=',ntsd,me,tmin,tmax
763 print*,'dvdt(i,k)=',dvdt(i_max,k_max)*dt_phs,vgrs(i_max,k_max)
764
765 tmax=-1.e-5
766 tmin=1.e5
767 do k=1,levs
768 k1=levs+1-k
769 do i=1,im
770! t_myj1=q_myj(i,k1)+rthblten(i,k1)*exner(i,k1)*dt_phs
771 t_myj1=q_myj(i,k1)+dqdt(i,k,1)*dt_phs
772! t_myj1=dqdt(i,k,1)*dt_phs
773 if(tmax.lt.t_myj1)then
774 tmax=t_myj1
775 i_max=i
776 k_max=k
777 end if
778 if(tmin.gt.t_myj1)then
779 tmin=t_myj1
780 i_min=i
781 k_min=k
782 end if
783 end do
784 end do
785 print*,'5after i_min,k_min,i_max,k_max=',i_min,k_min,i_max,k_max
786 print*,'ntsd,me,tmin,tmax=',ntsd,me,tmin,tmax
787 print*,'dqdt(i,k)=',dqdt(i_max,k_max,1)*dt_phs,qgrs(i_max,k_max,1)
788
789 tmax=-1.e-5
790 tmin=1.e5
791 do k=1,levs
792 k1=levs+1-k
793 do i=1,im
794! t_myj1=t_myj(i,k1)+rthblten(i,k1)*exner(i,k1)*dt_phs
795 t_myj1=cw(i,k1)+dqdt(i,k,ntcw)*dt_phs
796! t_myj1=dqdt(i,k,ntcw)*dt_phs
797 if(tmax.lt.t_myj1)then
798 tmax=t_myj1
799 i_max=i
800 k_max=k
801 end if
802 if(tmin.gt.t_myj1)then
803 tmin=t_myj1
804 i_min=i
805 k_min=k
806 end if
807 end do
808 end do
809 print*,'6after i_min,k_min,i_max,k_max=',i_min,k_min,i_max,k_max
810 print*,'ntsd,me,tmin,tmax=',ntsd,me,tmin,tmax
811 print*,'dqdt(i,k,ntcw)=',dqdt(i_max,k_max,ntcw)*dt_phs,qgrs(i_max,k_max,ntcw)
812
813 end if ! lprnt2
814
815! if (lprnt) then
816! print*
817! print*,"===Finished with myj_bl_driver; output:"
818! print*
819! endif
820
821 END SUBROUTINE myjpbl_wrapper_run
822
823!###=================================================================
824
825END MODULE myjpbl_wrapper