59 real(RP),
private,
parameter :: C1 = 0.19_rp
60 real(RP),
private,
parameter :: C2 = 0.51_rp
62 real(RP),
private,
parameter :: OneOverThree = 1.0_rp / 3.0_rp
63 real(RP),
private,
parameter :: TwoOverThree = 2.0_rp / 3.0_rp
64 real(RP),
private,
parameter :: FourOverThree = 4.0_rp / 3.0_rp
66 real(RP),
private,
allocatable :: delta(:,:,:)
68 integer,
private :: I_TKE
70 real(RP),
private :: ATMOS_PHY_TB_D1980_Km_MAX = 1000.0_rp
71 logical,
private :: ATMOS_PHY_TB_D1980_implicit = .false.
86 character(len=*),
intent(in) :: type_tb
87 integer,
intent(out) :: i_tke_out
91 log_info(
"ATMOS_PHY_TB_d1980_config",*)
'Setup'
92 log_info(
"ATMOS_PHY_TB_d1980_config",*)
'Tracers for Deardorff (1980) 1.5th TKE Model'
94 if ( type_tb /=
'D1980' )
then
95 log_error(
"ATMOS_PHY_TB_d1980_config",*)
'ATMOS_PHY_TB_TYPE is not D1980. Check!'
102 (/
'turbulent kinetic energy (D1980)' /), &
118 real(rp),
intent(in) :: cdz(
ka)
119 real(rp),
intent(in) :: cdx(
ia)
120 real(rp),
intent(in) :: cdy(
ja)
121 real(rp),
intent(in) :: cz (
ka,
ia,
ja)
123 namelist / param_atmos_phy_tb_d1980 / &
124 atmos_phy_tb_d1980_km_max, &
125 atmos_phy_tb_d1980_implicit
132 log_info(
"ATMOS_PHY_TB_d1980_setup",*)
'Setup'
133 log_info(
"ATMOS_PHY_TB_d1980_setup",*)
'Deardorff (1980) 1.5th TKE Model'
137 read(
io_fid_conf,nml=param_atmos_phy_tb_d1980,iostat=ierr)
139 log_info(
"ATMOS_PHY_TB_d1980_setup",*)
'Not found namelist. Default used.'
140 elseif( ierr > 0 )
then
141 log_error(
"ATMOS_PHY_TB_d1980_setup",*)
'Not appropriate names in namelist PARAM_ATMOS_PHY_TB_D1980. Check!'
144 log_nml(param_atmos_phy_tb_d1980)
146 allocate( delta(
ka,
ia,
ja) )
155 delta(
k,i,j) = ( cdz(
k) * cdx(i) * cdy(j) )**(1.0_rp/3.0_rp)
165 qflx_sgs_momz, qflx_sgs_momx, qflx_sgs_momy, &
166 qflx_sgs_rhot, qflx_sgs_rhoq, &
169 MOMZ, MOMX, MOMY, RHOT, DENS, QTRC, N2, &
170 SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_Q, &
171 GSQRT, J13G, J23G, J33G, MAPF, dt )
196 real(rp),
intent(out) :: qflx_sgs_momz(
ka,
ia,
ja,3)
197 real(rp),
intent(out) :: qflx_sgs_momx(
ka,
ia,
ja,3)
198 real(rp),
intent(out) :: qflx_sgs_momy(
ka,
ia,
ja,3)
199 real(rp),
intent(out) :: qflx_sgs_rhot(
ka,
ia,
ja,3)
200 real(rp),
intent(out) :: qflx_sgs_rhoq(
ka,
ia,
ja,3,
qa)
202 real(rp),
intent(inout) :: rhoq_t (
ka,
ia,
ja,
qa)
204 real(rp),
intent(out) :: km (
ka,
ia,
ja)
205 real(rp),
intent(out) :: ri (
ka,
ia,
ja)
206 real(rp),
intent(out) :: pr (
ka,
ia,
ja)
208 real(rp),
intent(in) :: momz (
ka,
ia,
ja)
209 real(rp),
intent(in) :: momx (
ka,
ia,
ja)
210 real(rp),
intent(in) :: momy (
ka,
ia,
ja)
211 real(rp),
intent(in) :: rhot (
ka,
ia,
ja)
212 real(rp),
intent(in) :: dens (
ka,
ia,
ja)
213 real(rp),
intent(in) :: qtrc (
ka,
ia,
ja,
qa)
214 real(rp),
intent(in) :: n2 (
ka,
ia,
ja)
216 real(rp),
intent(in) :: sflx_mw (
ia,
ja)
217 real(rp),
intent(in) :: sflx_mu (
ia,
ja)
218 real(rp),
intent(in) :: sflx_mv (
ia,
ja)
219 real(rp),
intent(in) :: sflx_sh (
ia,
ja)
220 real(rp),
intent(in) :: sflx_q (
ia,
ja,
qa)
222 real(rp),
intent(in) :: gsqrt (
ka,
ia,
ja,7)
223 real(rp),
intent(in) :: j13g (
ka,
ia,
ja,7)
224 real(rp),
intent(in) :: j23g (
ka,
ia,
ja,7)
225 real(rp),
intent(in) :: j33g
226 real(rp),
intent(in) :: mapf (
ia,
ja,2,4)
227 real(
dp),
intent(in) :: dt
230 real(rp) :: pott (
ka,
ia,
ja)
233 real(rp) :: s33_c (
ka,
ia,
ja)
234 real(rp) :: s11_c (
ka,
ia,
ja)
235 real(rp) :: s22_c (
ka,
ia,
ja)
236 real(rp) :: s31_c (
ka,
ia,
ja)
237 real(rp) :: s12_c (
ka,
ia,
ja)
238 real(rp) :: s23_c (
ka,
ia,
ja)
239 real(rp) :: s12_z (
ka,
ia,
ja)
240 real(rp) :: s23_x (
ka,
ia,
ja)
241 real(rp) :: s31_y (
ka,
ia,
ja)
247 real(rp) :: qflx_tke(
ka,
ia,
ja,3)
249 real(rp) :: tend(
ka,
ia,
ja)
259 integer ::
k, i, j, iq
262 log_progress(*)
'atmosphere / physics / turbulence / D1980'
265 qflx_sgs_momz(:,:,:,:) = undef
266 qflx_sgs_momx(:,:,:,:) = undef
267 qflx_sgs_momy(:,:,:,:) = undef
268 qflx_sgs_rhot(:,:,:,:) = undef
269 qflx_sgs_rhoq(:,:,:,:,:) = undef
285 call check( __line__, rhot(
k,i,j) )
286 call check( __line__, dens(
k,i,j) )
288 pott(
k,i,j) = rhot(
k,i,j) / dens(
k,i,j)
293 i = iundef; j = iundef;
k = iundef
298 call calc_strain_tensor( &
299 s33_c, s11_c, s22_c, &
300 s31_c, s12_c, s23_c, &
301 s12_z, s23_x, s31_y, &
303 dens, momz, momx, momy, &
304 gsqrt, j13g, j23g, j33g, mapf )
317 call check( __line__, pott(
k+1,i,j) )
318 call check( __line__, pott(
k,i,j) )
319 call check( __line__, pott(
k-1,i,j) )
320 call check( __line__, fdz(
k) )
321 call check( __line__, fdz(
k-1) )
322 call check( __line__, s2(
k,i,j) )
324 ri(
k,i,j) = n2(
k,i,j) / s2(
k,i,j)
329 i = iundef; j = iundef;
k = iundef
337 if ( n2(
k,i,j) > 1e-10_rp )
then
338 l(
k,i,j) = max( min( 0.76_rp * sqrt(qtrc(
k,i,j,i_tke)/n2(
k,i,j)), delta(
k,i,j) ), 1e-10_rp )
340 l(
k,i,j) = delta(
k,i,j)
346 i = iundef; j = iundef;
k = iundef
353 km(
k,i,j) = min( 0.10_rp * l(
k,i,j) * sqrt(qtrc(
k,i,j,i_tke)), atmos_phy_tb_d1980_km_max )
354 pr(
k,i,j) = 1.0_rp / ( 1.0_rp + 2.0_rp * l(
k,i,j) / delta(
k,i,j) )
355 kh(
k,i,j) = km(
k,i,j) / pr(
k,i,j)
360 i = iundef; j = iundef;
k = iundef
370 call check( __line__, dens(
k,i,j) )
371 call check( __line__, km(
k,i,j) )
372 call check( __line__, s33_c(
k,i,j) )
374 qflx_sgs_momz(
k,i,j,
zdir) = dens(
k,i,j) * ( - 2.0_rp * km(
k,i,j) * s33_c(
k,i,j) )
379 i = iundef; j = iundef;
k = iundef
384 qflx_sgs_momz(
ks,i,j,
zdir) = 0.0_rp
385 qflx_sgs_momz(
ke,i,j,
zdir) = 0.0_rp
389 i = iundef; j = iundef;
k = iundef
397 call check( __line__, dens(
k,i,j) )
398 call check( __line__, dens(
k,i+1,j) )
399 call check( __line__, dens(
k+1,i,j) )
400 call check( __line__, dens(
k+1,i+1,j) )
401 call check( __line__, km(
k,i,j) )
402 call check( __line__, km(
k,i+1,j) )
403 call check( __line__, km(
k+1,i,j) )
404 call check( __line__, km(
k+1,i+1,j) )
405 call check( __line__, s31_y(
k,i,j) )
407 qflx_sgs_momz(
k,i,j,
xdir) = - 0.125_rp &
408 * ( dens(
k,i,j)+dens(
k+1,i,j)+dens(
k,i+1,j)+dens(
k+1,i+1,j) ) &
409 * ( km(
k,i,j)+km(
k+1,i,j)+km(
k,i+1,j)+km(
k+1,i+1,j)) &
415 i = iundef; j = iundef;
k = iundef
423 call check( __line__, dens(
k,i,j) )
424 call check( __line__, dens(
k,i,j+1) )
425 call check( __line__, dens(
k+1,i,j) )
426 call check( __line__, dens(
k+1,i,j+1) )
427 call check( __line__, km(
k,i,j) )
428 call check( __line__, km(
k,i,j+1) )
429 call check( __line__, km(
k+1,i,j) )
430 call check( __line__, km(
k+1,i,j+1) )
431 call check( __line__, s23_x(
k,i,j) )
433 qflx_sgs_momz(
k,i,j,
ydir) = - 0.125_rp &
434 * ( dens(
k,i,j)+dens(
k+1,i,j)+dens(
k,i,j+1)+dens(
k+1,i,j+1) ) &
435 * ( km(
k,i,j)+km(
k+1,i,j)+km(
k,i,j+1)+km(
k+1,i,j+1) ) &
441 i = iundef; j = iundef;
k = iundef
444 if ( atmos_phy_tb_d1980_implicit )
then
445 call calc_tend_momz( tend, &
447 gsqrt, j13g, j23g, j33g, mapf, &
455 ap = - fouroverthree * dt &
456 * dens(
ks+1,i,j)*km(
ks+1,i,j) &
460 b(
ks,i,j) = - a(
ks,i,j) + 0.5_rp * ( dens(
ks,i,j)+dens(
ks+1,i,j) )
462 c(
k,i,j) = ap * rfdz(
k+1) / gsqrt(
k+1,i,j,
i_xyw)
463 ap = - fouroverthree * dt &
464 * dens(
k+1,i,j)*km(
k+1,i,j) &
465 * rcdz(
k+1) / gsqrt(
k+1,i,j,
i_xyz)
466 a(
k,i,j) = ap * rfdz(
k) / gsqrt(
k,i,j,
i_xyw)
467 b(
k,i,j) = - a(
k,i,j) - c(
k,i,j) + 0.5_rp * ( dens(
k,i,j)+dens(
k+1,i,j) )
471 b(
ke-1,i,j) = - c(
ke-1,i,j) + 0.5_rp * ( dens(
ke-1,i,j)+dens(
ke,i,j) )
477 call diffusion_solver( &
479 a(:,i,j), b(:,i,j), c(:,i,j), d, &
483 qflx_sgs_momz(
k,i,j,
zdir) = qflx_sgs_momz(
k,i,j,
zdir) &
484 - fouroverthree * dens(
k,i,j) * km(
k,i,j) * dt &
485 * ( tend(
k,i,j) - tend(
k-1,i,j) ) * rcdz(
k) / gsqrt(
k,i,j,
i_xyz)
500 call check( __line__, dens(
k,i,j) )
501 call check( __line__, dens(
k,i+1,j) )
502 call check( __line__, dens(
k+1,i,j) )
503 call check( __line__, dens(
k+1,i+1,j) )
504 call check( __line__, km(
k,i,j) )
505 call check( __line__, km(
k,i+1,j) )
506 call check( __line__, km(
k+1,i,j) )
507 call check( __line__, km(
k+1,i+1,j) )
508 call check( __line__, s31_y(
k,i,j) )
510 qflx_sgs_momx(
k,i,j,
zdir) = - 0.125_rp &
511 * ( dens(
k,i,j)+dens(
k+1,i,j)+dens(
k,i+1,j)+dens(
k+1,i+1,j) ) &
512 * ( km(
k,i,j)+km(
k+1,i,j)+km(
k,i+1,j)+km(
k+1,i+1,j) ) &
518 i = iundef; j = iundef;
k = iundef
523 qflx_sgs_momx(
ks-1,i,j,
zdir) = 0.0_rp
524 qflx_sgs_momx(
ke ,i,j,
zdir) = 0.0_rp
528 i = iundef; j = iundef;
k = iundef
536 call check( __line__, dens(
k,i,j) )
537 call check( __line__, km(
k,i,j) )
538 call check( __line__, s11_c(
k,i,j) )
540 qflx_sgs_momx(
k,i,j,
xdir) = dens(
k,i,j) * ( - 2.0_rp * km(
k,i,j) * s11_c(
k,i,j) )
545 i = iundef; j = iundef;
k = iundef
553 call check( __line__, dens(
k,i,j) )
554 call check( __line__, dens(
k,i+1,j) )
555 call check( __line__, dens(
k,i,j+1) )
556 call check( __line__, dens(
k,i+1,j+1) )
557 call check( __line__, km(
k,i,j) )
558 call check( __line__, km(
k,i+1,j) )
559 call check( __line__, km(
k,i,j+1) )
560 call check( __line__, km(
k,i+1,j+1) )
561 call check( __line__, s12_z(
k,i,j) )
563 qflx_sgs_momx(
k,i,j,
ydir) = - 0.125_rp &
564 * ( dens(
k,i,j)+dens(
k,i+1,j)+dens(
k,i,j+1)+dens(
k,i+1,j+1) ) &
565 * ( km(
k,i,j)+km(
k,i+1,j)+km(
k,i,j+1)+km(
k,i+1,j+1) ) &
571 i = iundef; j = iundef;
k = iundef
574 if ( atmos_phy_tb_d1980_implicit )
then
575 call calc_tend_momx( tend, &
577 gsqrt, j13g, j23g, j33g, mapf, &
585 ap = - dt * 0.25_rp * ( dens(
ks ,i ,j)*km(
ks ,i ,j) &
586 + dens(
ks+1,i ,j)*km(
ks+1,i ,j) &
587 + dens(
ks ,i+1,j)*km(
ks ,i+1,j) &
588 + dens(
ks+1,i+1,j)*km(
ks+1,i+1,j) ) &
592 b(
ks,i,j) = - a(
ks,i,j) + 0.5_rp * ( dens(
ks,i,j)+dens(
ks,i+1,j) )
594 c(
k,i,j) = ap * rcdz(
k) / gsqrt(
k,i,j,
i_uyz)
595 ap = - dt * 0.25_rp * ( dens(
k ,i ,j)*km(
k ,i ,j) &
596 + dens(
k+1,i ,j)*km(
k+1,i ,j) &
597 + dens(
k ,i+1,j)*km(
k ,i+1,j) &
598 + dens(
k+1,i+1,j)*km(
k+1,i+1,j) ) &
600 a(
k,i,j) = ap * rcdz(
k) / gsqrt(
k,i,j,
i_uyz)
601 b(
k,i,j) = - a(
k,i,j) - c(
k,i,j) + 0.5_rp * ( dens(
k,i,j)+dens(
k,i+1,j) )
605 b(
ke,i,j) = - c(
ke,i,j) + 0.5_rp * ( dens(
ke,i,j)+dens(
ke,i+1,j) )
611 call diffusion_solver( &
613 a(:,i,j), b(:,i,j), c(:,i,j), d, &
617 qflx_sgs_momx(
k,i,j,
zdir) = qflx_sgs_momx(
k,i,j,
zdir) &
618 - 0.25_rp * ( dens(
k ,i ,j)*km(
k ,i ,j) &
619 + dens(
k+1,i ,j)*km(
k+1,i ,j) &
620 + dens(
k ,i+1,j)*km(
k ,i+1,j) &
621 + dens(
k+1,i+1,j)*km(
k+1,i+1,j) ) &
622 * dt * ( tend(
k+1,i,j) - tend(
k,i,j) ) * rfdz(
k) / gsqrt(
k,i,j,
i_uyw)
637 call check( __line__, dens(
k,i,j) )
638 call check( __line__, dens(
k,i,j+1) )
639 call check( __line__, dens(
k+1,i,j) )
640 call check( __line__, dens(
k+1,i,j+1) )
641 call check( __line__, km(
k,i,j) )
642 call check( __line__, km(
k,i,j+1) )
643 call check( __line__, km(
k+1,i,j) )
644 call check( __line__, km(
k+1,i,j+1) )
645 call check( __line__, s23_x(
k,i,j) )
647 qflx_sgs_momy(
k,i,j,
zdir) = - 0.125_rp &
648 * ( dens(
k,i,j)+dens(
k+1,i,j)+dens(
k,i,j+1)+dens(
k+1,i,j+1) ) &
649 * ( km(
k,i,j)+km(
k+1,i,j)+km(
k,i,j+1)+km(
k+1,i,j+1) ) &
655 i = iundef; j = iundef;
k = iundef
660 qflx_sgs_momy(
ks-1,i,j,
zdir) = 0.0_rp
661 qflx_sgs_momy(
ke ,i,j,
zdir) = 0.0_rp
665 i = iundef; j = iundef;
k = iundef
674 call check( __line__, dens(
k,i,j) )
675 call check( __line__, dens(
k,i+1,j) )
676 call check( __line__, dens(
k,i,j+1) )
677 call check( __line__, dens(
k,i+1,j+1) )
678 call check( __line__, km(
k,i,j) )
679 call check( __line__, km(
k,i+1,j) )
680 call check( __line__, km(
k,i,j+1) )
681 call check( __line__, km(
k,i+1,j+1) )
682 call check( __line__, s12_z(
k,i,j) )
684 qflx_sgs_momy(
k,i,j,
xdir) = - 0.125_rp &
685 * ( dens(
k,i,j)+dens(
k,i+1,j)+dens(
k,i,j+1)+dens(
k,i+1,j+1) ) &
686 * ( km(
k,i,j)+km(
k,i+1,j)+km(
k,i,j+1)+km(
k,i+1,j+1) ) &
692 i = iundef; j = iundef;
k = iundef
701 call check( __line__, dens(
k,i,j) )
702 call check( __line__, km(
k,i,j) )
703 call check( __line__, s22_c(
k,i,j) )
705 qflx_sgs_momy(
k,i,j,
ydir) = dens(
k,i,j) * ( - 2.0_rp * km(
k,i,j) * s22_c(
k,i,j) )
710 i = iundef; j = iundef;
k = iundef
713 if ( atmos_phy_tb_d1980_implicit )
then
714 call calc_tend_momy( tend, &
716 gsqrt, j13g, j23g, j33g, mapf, &
724 ap = - dt * 0.25_rp * ( dens(
ks ,i,j )*km(
ks ,i,j ) &
725 + dens(
ks+1,i,j )*km(
ks+1,i,j ) &
726 + dens(
ks ,i,j+1)*km(
ks ,i,j+1) &
727 + dens(
ks+1,i,j+1)*km(
ks+1,i,j+1) ) &
731 b(
ks,i,j) = - a(
ks,i,j) + 0.5_rp * ( dens(
ks,i,j)+dens(
ks,i,j+1) )
733 c(
k,i,j) = ap * rcdz(
k) / gsqrt(
k,i,j,
i_xvz)
734 ap = - dt * 0.25_rp * ( dens(
k ,i,j )*km(
k ,i,j ) &
735 + dens(
k+1,i,j )*km(
k+1,i,j ) &
736 + dens(
k ,i,j+1)*km(
k ,i,j+1) &
737 + dens(
k+1,i,j+1)*km(
k+1,i,j+1) ) &
739 a(
k,i,j) = ap * rcdz(
k) / gsqrt(
k,i,j,
i_xvz)
740 b(
k,i,j) = - a(
k,i,j) - c(
k,i,j) + 0.5_rp * ( dens(
k,i,j)+dens(
k,i,j+1) )
744 b(
ke,i,j) = - c(
ke,i,j) + 0.5_rp * ( dens(
ke,i,j)+dens(
ke,i,j+1) )
750 call diffusion_solver( &
752 a(:,i,j), b(:,i,j), c(:,i,j), d, &
756 qflx_sgs_momy(
k,i,j,
zdir) = qflx_sgs_momy(
k,i,j,
zdir) &
757 - 0.25_rp * ( dens(
k ,i,j )*km(
k ,i,j ) &
758 + dens(
k+1,i,j )*km(
k+1,i,j ) &
759 + dens(
k ,i,j+1)*km(
k ,i,j+1) &
760 + dens(
k+1,i,j+1)*km(
k+1,i,j+1) ) &
761 * dt * ( tend(
k+1,i,j) - tend(
k,i,j) ) * rfdz(
k) / gsqrt(
k,i,j,
i_xvw)
771 if ( atmos_phy_tb_d1980_implicit )
then
778 ap = - dt * 0.25_rp * ( dens(
ks,i,j)+dens(
ks+1,i,j) ) &
779 * ( kh(
ks,i,j)+kh(
ks+1,i,j) ) &
783 b(
ks,i,j) = - a(
ks,i,j) + dens(
ks,i,j)
785 c(
k,i,j) = ap * rcdz(
k) / gsqrt(
k,i,j,
i_xyz)
786 ap = - dt * 0.25_rp * ( dens(
k,i,j)+dens(
k+1,i,j) ) &
787 * ( kh(
k,i,j)+kh(
k+1,i,j) ) &
789 a(
k,i,j) = ap * rcdz(
k) / gsqrt(
k,i,j,
i_xyz)
790 b(
k,i,j) = - a(
k,i,j) - c(
k,i,j) + dens(
k,i,j)
794 b(
ke,i,j) = - c(
ke,i,j) + dens(
ke,i,j)
801 call calc_flux_phi( &
803 dens, pott, kh, 1.0_rp, &
804 gsqrt, j13g, j23g, j33g, mapf, &
806 atmos_phy_tb_d1980_implicit, &
817 if ( iq == i_tke )
then
818 qflx_sgs_rhoq(:,:,:,:,iq) = 0.0_rp
828 call calc_flux_phi( &
829 qflx_sgs_rhoq(:,:,:,:,iq), &
830 dens, qtrc(:,:,:,iq), kh, 1.0_rp, &
831 gsqrt, j13g, j23g, j33g, mapf, &
833 atmos_phy_tb_d1980_implicit, &
840 iis = iundef; iie = iundef; jjs = iundef; jje = iundef
854 if ( atmos_phy_tb_d1980_implicit )
then
861 ap = - dt * 0.25_rp * ( dens(
ks,i,j)+dens(
ks+1,i,j) ) &
862 * 2.0_rp * ( km(
ks,i,j)+km(
ks+1,i,j) ) &
866 b(
ks,i,j) = - a(
ks,i,j) + dens(
ks,i,j)
868 c(
k,i,j) = ap * rcdz(
k) / gsqrt(
k,i,j,
i_xyz)
869 ap = - dt * 0.25_rp * ( dens(
k,i,j)+dens(
k+1,i,j) ) &
870 * 2.0_rp * ( km(
k,i,j)+km(
k+1,i,j) ) &
872 a(
k,i,j) = ap * rcdz(
k) / gsqrt(
k,i,j,
i_xyz)
873 b(
k,i,j) = - a(
k,i,j) - c(
k,i,j) + dens(
k,i,j)
877 b(
ke,i,j) = - c(
ke,i,j) + dens(
ke,i,j)
884 call calc_flux_phi( &
886 dens, qtrc(:,:,:,i_tke), km, 2.0_rp, &
887 gsqrt, j13g, j23g, j33g, mapf, &
889 atmos_phy_tb_d1980_implicit, &
893 call calc_tend_phi( rhoq_t(:,:,:,i_tke), &
895 gsqrt, j13g, j23g, j33g, mapf, &
901 rhoq_t(
k,i,j,i_tke) = max( &
902 rhoq_t(
k,i,j,i_tke) &
903 + ( km(
k,i,j) * s2(
k,i,j) &
904 - kh(
k,i,j) * n2(
k,i,j) &
905 - ( c1 + c2*l(
k,i,j)/delta(
k,i,j) ) * qtrc(
k,i,j,i_tke)**(1.5_rp) / l(
k,i,j) ) * dens(
k,i,j), &
906 - qtrc(
k,i,j,i_tke) * dens(
k,i,j) / real(dt,kind=rp) )
911 i = iundef; j = iundef;
k = iundef