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)
242 real(RP) :: S2 (KA,IA,JA)
244 real(RP) :: Kh (KA,IA,JA)
245 real(RP) :: l (KA,IA,JA)
247 real(RP) :: qflx_tke(KA,IA,JA,3)
249 real(RP) :: TEND(KA,IA,JA)
250 real(RP) :: a(KA,IA,JA)
251 real(RP) :: b(KA,IA,JA)
252 real(RP) :: c(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) &
457 * rcdz(ks+1) / gsqrt(ks+1,i,j,
i_xyz)
458 a(ks,i,j) = ap * rfdz(ks) / gsqrt(ks,i,j,
i_xyw)
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) )
470 c(ke-1,i,j) = ap * rfdz(ke) / gsqrt(ke,i,j,
i_xyw)
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) ) &
589 * rfdz(ks) / gsqrt(ks,i,j,
i_uyw)
590 a(ks,i,j) = ap * rcdz(ks) / gsqrt(ks,i,j,
i_uyz)
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) ) &
599 * rfdz(k) / gsqrt(k,i,j,
i_uyw)
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) )
604 c(ke,i,j) = ap * rcdz(ke) / gsqrt(ke,i,j,
i_uyz)
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) ) &
728 * rfdz(ks) / gsqrt(ks,i,j,
i_xvw)
729 a(ks,i,j) = ap * rcdz(ks) / gsqrt(ks,i,j,
i_xvz)
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) ) &
738 * rfdz(k) / gsqrt(k,i,j,
i_xvw)
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) )
743 c(ke,i,j) = ap * rcdz(ke) / gsqrt(ke,i,j,
i_xvz)
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) ) &
780 * rfdz(ks) / gsqrt(ks,i,j,
i_xyw)
781 a(ks,i,j) = ap * rcdz(ks) / gsqrt(ks,i,j,
i_xyz)
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) ) &
788 * rfdz(k) / gsqrt(k,i,j,
i_xyw)
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)
793 c(ke,i,j) = ap * rcdz(ke) / gsqrt(ke,i,j,
i_xyz)
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) ) &
863 * rfdz(ks) / gsqrt(ks,i,j,
i_xyw)
864 a(ks,i,j) = ap * rcdz(ks) / gsqrt(ks,i,j,
i_xyz)
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) ) &
871 * rfdz(k) / gsqrt(k,i,j,
i_xyw)
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)
876 c(ke,i,j) = ap * rcdz(ke) / gsqrt(ke,i,j,
i_xyz)
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