196 matrix_solver_tridiagonal
207 real(RP),
intent(out) :: qflx_sgs_momz(KA,IA,JA,3)
208 real(RP),
intent(out) :: qflx_sgs_momx(KA,IA,JA,3)
209 real(RP),
intent(out) :: qflx_sgs_momy(KA,IA,JA,3)
210 real(RP),
intent(out) :: qflx_sgs_rhot(KA,IA,JA,3)
211 real(RP),
intent(out) :: qflx_sgs_rhoq(KA,IA,JA,3,QA)
213 real(RP),
intent(inout) :: RHOQ_t (KA,IA,JA,QA)
215 real(RP),
intent(out) :: Km (KA,IA,JA)
216 real(RP),
intent(out) :: Ri (KA,IA,JA)
217 real(RP),
intent(out) :: Pr (KA,IA,JA)
219 real(RP),
intent(in) :: MOMZ (KA,IA,JA)
220 real(RP),
intent(in) :: MOMX (KA,IA,JA)
221 real(RP),
intent(in) :: MOMY (KA,IA,JA)
222 real(RP),
intent(in) :: RHOT (KA,IA,JA)
223 real(RP),
intent(in) :: DENS (KA,IA,JA)
224 real(RP),
intent(in) :: QTRC (KA,IA,JA,QA)
225 real(RP),
intent(in) :: N2 (KA,IA,JA)
227 real(RP),
intent(in) :: SFLX_MW (IA,JA)
228 real(RP),
intent(in) :: SFLX_MU (IA,JA)
229 real(RP),
intent(in) :: SFLX_MV (IA,JA)
230 real(RP),
intent(in) :: SFLX_SH (IA,JA)
231 real(RP),
intent(in) :: SFLX_Q (IA,JA,QA)
233 real(RP),
intent(in) :: GSQRT (KA,IA,JA,7)
234 real(RP),
intent(in) :: J13G (KA,IA,JA,7)
235 real(RP),
intent(in) :: J23G (KA,IA,JA,7)
236 real(RP),
intent(in) :: J33G
237 real(RP),
intent(in) :: MAPF (IA,JA,2,4)
238 real(DP),
intent(in) :: dt
241 real(RP) :: POTT (KA,IA,JA)
244 real(RP) :: S33_C (KA,IA,JA)
245 real(RP) :: S11_C (KA,IA,JA)
246 real(RP) :: S22_C (KA,IA,JA)
247 real(RP) :: S31_C (KA,IA,JA)
248 real(RP) :: S12_C (KA,IA,JA)
249 real(RP) :: S23_C (KA,IA,JA)
250 real(RP) :: S12_Z (KA,IA,JA)
251 real(RP) :: S23_X (KA,IA,JA)
252 real(RP) :: S31_Y (KA,IA,JA)
253 real(RP) :: S2 (KA,IA,JA)
255 real(RP) :: Kh (KA,IA,JA)
256 real(RP) :: l (KA,IA,JA)
258 real(RP) :: qflx_tke(KA,IA,JA,3)
260 real(RP) :: TEND(KA,IA,JA)
261 real(RP) :: TEND2(KA,IA,JA)
262 real(RP) :: a(KA,IA,JA)
263 real(RP) :: b(KA,IA,JA)
264 real(RP) :: c(KA,IA,JA)
270 integer :: k, i, j, iq
273 log_progress(*)
'atmosphere / physics / turbulence / D1980'
276 qflx_sgs_momz(:,:,:,:) = undef
277 qflx_sgs_momx(:,:,:,:) = undef
278 qflx_sgs_momy(:,:,:,:) = undef
279 qflx_sgs_rhot(:,:,:,:) = undef
280 qflx_sgs_rhoq(:,:,:,:,:) = undef
296 call check( __line__, rhot(k,i,j) )
297 call check( __line__, dens(k,i,j) )
299 pott(k,i,j) = rhot(k,i,j) / dens(k,i,j)
304 i = iundef; j = iundef; k = iundef
309 call calc_strain_tensor( &
310 s33_c, s11_c, s22_c, &
311 s31_c, s12_c, s23_c, &
312 s12_z, s23_x, s31_y, &
314 dens, momz, momx, momy, &
315 gsqrt, j13g, j23g, j33g, mapf )
328 call check( __line__, pott(k+1,i,j) )
329 call check( __line__, pott(k,i,j) )
330 call check( __line__, pott(k-1,i,j) )
331 call check( __line__, fdz(k) )
332 call check( __line__, fdz(k-1) )
333 call check( __line__, s2(k,i,j) )
335 ri(k,i,j) = n2(k,i,j) / s2(k,i,j)
340 i = iundef; j = iundef; k = iundef
348 if ( n2(k,i,j) > 1e-10_rp )
then
349 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 )
351 l(k,i,j) = delta(k,i,j)
357 i = iundef; j = iundef; k = iundef
364 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 )
365 pr(k,i,j) = 1.0_rp / ( 1.0_rp + 2.0_rp * l(k,i,j) / delta(k,i,j) )
366 kh(k,i,j) = km(k,i,j) / pr(k,i,j)
371 i = iundef; j = iundef; k = iundef
381 call check( __line__, dens(k,i,j) )
382 call check( __line__, km(k,i,j) )
383 call check( __line__, s33_c(k,i,j) )
385 qflx_sgs_momz(k,i,j,
zdir) = dens(k,i,j) * ( - 2.0_rp * km(k,i,j) * s33_c(k,i,j) )
390 i = iundef; j = iundef; k = iundef
395 qflx_sgs_momz(ks,i,j,
zdir) = 0.0_rp
396 qflx_sgs_momz(ke,i,j,
zdir) = 0.0_rp
400 i = iundef; j = iundef; k = iundef
408 call check( __line__, dens(k,i,j) )
409 call check( __line__, dens(k,i+1,j) )
410 call check( __line__, dens(k+1,i,j) )
411 call check( __line__, dens(k+1,i+1,j) )
412 call check( __line__, km(k,i,j) )
413 call check( __line__, km(k,i+1,j) )
414 call check( __line__, km(k+1,i,j) )
415 call check( __line__, km(k+1,i+1,j) )
416 call check( __line__, s31_y(k,i,j) )
418 qflx_sgs_momz(k,i,j,
xdir) = - 0.125_rp &
419 * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i+1,j)+dens(k+1,i+1,j) ) &
420 * ( km(k,i,j)+km(k+1,i,j)+km(k,i+1,j)+km(k+1,i+1,j)) &
426 i = iundef; j = iundef; k = iundef
434 call check( __line__, dens(k,i,j) )
435 call check( __line__, dens(k,i,j+1) )
436 call check( __line__, dens(k+1,i,j) )
437 call check( __line__, dens(k+1,i,j+1) )
438 call check( __line__, km(k,i,j) )
439 call check( __line__, km(k,i,j+1) )
440 call check( __line__, km(k+1,i,j) )
441 call check( __line__, km(k+1,i,j+1) )
442 call check( __line__, s23_x(k,i,j) )
444 qflx_sgs_momz(k,i,j,
ydir) = - 0.125_rp &
445 * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i,j+1)+dens(k+1,i,j+1) ) &
446 * ( km(k,i,j)+km(k+1,i,j)+km(k,i,j+1)+km(k+1,i,j+1) ) &
452 i = iundef; j = iundef; k = iundef
455 if ( atmos_phy_tb_d1980_implicit )
then
456 call calc_tend_momz( tend, &
458 gsqrt, j13g, j23g, j33g, mapf, &
466 ap = - fouroverthree * dt &
467 * dens(ks+1,i,j)*km(ks+1,i,j) &
468 * rcdz(ks+1) / gsqrt(ks+1,i,j,
i_xyz)
469 a(ks,i,j) = ap * rfdz(ks) / gsqrt(ks,i,j,
i_xyw)
471 b(ks,i,j) = - a(ks,i,j) + 0.5_rp * ( dens(ks,i,j)+dens(ks+1,i,j) )
473 c(k,i,j) = ap * rfdz(k+1) / gsqrt(k+1,i,j,
i_xyw)
474 ap = - fouroverthree * dt &
475 * dens(k+1,i,j)*km(k+1,i,j) &
476 * rcdz(k+1) / gsqrt(k+1,i,j,
i_xyz)
477 a(k,i,j) = ap * rfdz(k) / gsqrt(k,i,j,
i_xyw)
478 b(k,i,j) = - a(k,i,j) - c(k,i,j) + 0.5_rp * ( dens(k,i,j)+dens(k+1,i,j) )
481 c(ke-1,i,j) = ap * rfdz(ke) / gsqrt(ke,i,j,
i_xyw)
482 b(ke-1,i,j) = - c(ke-1,i,j) + 0.5_rp * ( dens(ke-1,i,j)+dens(ke,i,j) )
486 call matrix_solver_tridiagonal( ka, ks, ke-1, ia, iis, iie, ja, jjs, jje, &
487 a(:,:,:), b(:,:,:), c(:,:,:), tend(:,:,:), &
494 qflx_sgs_momz(k,i,j,
zdir) = qflx_sgs_momz(k,i,j,
zdir) &
495 - fouroverthree * dens(k,i,j) * km(k,i,j) * dt &
496 * ( tend2(k,i,j) - tend2(k-1,i,j) ) * rcdz(k) / gsqrt(k,i,j,
i_xyz)
511 call check( __line__, dens(k,i,j) )
512 call check( __line__, dens(k,i+1,j) )
513 call check( __line__, dens(k+1,i,j) )
514 call check( __line__, dens(k+1,i+1,j) )
515 call check( __line__, km(k,i,j) )
516 call check( __line__, km(k,i+1,j) )
517 call check( __line__, km(k+1,i,j) )
518 call check( __line__, km(k+1,i+1,j) )
519 call check( __line__, s31_y(k,i,j) )
521 qflx_sgs_momx(k,i,j,
zdir) = - 0.125_rp &
522 * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i+1,j)+dens(k+1,i+1,j) ) &
523 * ( km(k,i,j)+km(k+1,i,j)+km(k,i+1,j)+km(k+1,i+1,j) ) &
529 i = iundef; j = iundef; k = iundef
534 qflx_sgs_momx(ks-1,i,j,
zdir) = 0.0_rp
535 qflx_sgs_momx(ke ,i,j,
zdir) = 0.0_rp
539 i = iundef; j = iundef; k = iundef
547 call check( __line__, dens(k,i,j) )
548 call check( __line__, km(k,i,j) )
549 call check( __line__, s11_c(k,i,j) )
551 qflx_sgs_momx(k,i,j,
xdir) = dens(k,i,j) * ( - 2.0_rp * km(k,i,j) * s11_c(k,i,j) )
556 i = iundef; j = iundef; k = iundef
564 call check( __line__, dens(k,i,j) )
565 call check( __line__, dens(k,i+1,j) )
566 call check( __line__, dens(k,i,j+1) )
567 call check( __line__, dens(k,i+1,j+1) )
568 call check( __line__, km(k,i,j) )
569 call check( __line__, km(k,i+1,j) )
570 call check( __line__, km(k,i,j+1) )
571 call check( __line__, km(k,i+1,j+1) )
572 call check( __line__, s12_z(k,i,j) )
574 qflx_sgs_momx(k,i,j,
ydir) = - 0.125_rp &
575 * ( dens(k,i,j)+dens(k,i+1,j)+dens(k,i,j+1)+dens(k,i+1,j+1) ) &
576 * ( km(k,i,j)+km(k,i+1,j)+km(k,i,j+1)+km(k,i+1,j+1) ) &
582 i = iundef; j = iundef; k = iundef
585 if ( atmos_phy_tb_d1980_implicit )
then
586 call calc_tend_momx( tend, &
588 gsqrt, j13g, j23g, j33g, mapf, &
596 ap = - dt * 0.25_rp * ( dens(ks ,i ,j)*km(ks ,i ,j) &
597 + dens(ks+1,i ,j)*km(ks+1,i ,j) &
598 + dens(ks ,i+1,j)*km(ks ,i+1,j) &
599 + dens(ks+1,i+1,j)*km(ks+1,i+1,j) ) &
600 * rfdz(ks) / gsqrt(ks,i,j,
i_uyw)
601 a(ks,i,j) = ap * rcdz(ks) / gsqrt(ks,i,j,
i_uyz)
603 b(ks,i,j) = - a(ks,i,j) + 0.5_rp * ( dens(ks,i,j)+dens(ks,i+1,j) )
605 c(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,
i_uyz)
606 ap = - dt * 0.25_rp * ( dens(k ,i ,j)*km(k ,i ,j) &
607 + dens(k+1,i ,j)*km(k+1,i ,j) &
608 + dens(k ,i+1,j)*km(k ,i+1,j) &
609 + dens(k+1,i+1,j)*km(k+1,i+1,j) ) &
610 * rfdz(k) / gsqrt(k,i,j,
i_uyw)
611 a(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,
i_uyz)
612 b(k,i,j) = - a(k,i,j) - c(k,i,j) + 0.5_rp * ( dens(k,i,j)+dens(k,i+1,j) )
615 c(ke,i,j) = ap * rcdz(ke) / gsqrt(ke,i,j,
i_uyz)
616 b(ke,i,j) = - c(ke,i,j) + 0.5_rp * ( dens(ke,i,j)+dens(ke,i+1,j) )
620 call matrix_solver_tridiagonal( ka, ks, ke, ia, iis, iie, ja, jjs, jje, &
621 a(:,:,:), b(:,:,:), c(:,:,:), tend(:,:,:), &
628 qflx_sgs_momx(k,i,j,
zdir) = qflx_sgs_momx(k,i,j,
zdir) &
629 - 0.25_rp * ( dens(k ,i ,j)*km(k ,i ,j) &
630 + dens(k+1,i ,j)*km(k+1,i ,j) &
631 + dens(k ,i+1,j)*km(k ,i+1,j) &
632 + dens(k+1,i+1,j)*km(k+1,i+1,j) ) &
633 * dt * ( tend2(k+1,i,j) - tend2(k,i,j) ) * rfdz(k) / gsqrt(k,i,j,
i_uyw)
648 call check( __line__, dens(k,i,j) )
649 call check( __line__, dens(k,i,j+1) )
650 call check( __line__, dens(k+1,i,j) )
651 call check( __line__, dens(k+1,i,j+1) )
652 call check( __line__, km(k,i,j) )
653 call check( __line__, km(k,i,j+1) )
654 call check( __line__, km(k+1,i,j) )
655 call check( __line__, km(k+1,i,j+1) )
656 call check( __line__, s23_x(k,i,j) )
658 qflx_sgs_momy(k,i,j,
zdir) = - 0.125_rp &
659 * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i,j+1)+dens(k+1,i,j+1) ) &
660 * ( km(k,i,j)+km(k+1,i,j)+km(k,i,j+1)+km(k+1,i,j+1) ) &
666 i = iundef; j = iundef; k = iundef
671 qflx_sgs_momy(ks-1,i,j,
zdir) = 0.0_rp
672 qflx_sgs_momy(ke ,i,j,
zdir) = 0.0_rp
676 i = iundef; j = iundef; k = iundef
685 call check( __line__, dens(k,i,j) )
686 call check( __line__, dens(k,i+1,j) )
687 call check( __line__, dens(k,i,j+1) )
688 call check( __line__, dens(k,i+1,j+1) )
689 call check( __line__, km(k,i,j) )
690 call check( __line__, km(k,i+1,j) )
691 call check( __line__, km(k,i,j+1) )
692 call check( __line__, km(k,i+1,j+1) )
693 call check( __line__, s12_z(k,i,j) )
695 qflx_sgs_momy(k,i,j,
xdir) = - 0.125_rp &
696 * ( dens(k,i,j)+dens(k,i+1,j)+dens(k,i,j+1)+dens(k,i+1,j+1) ) &
697 * ( km(k,i,j)+km(k,i+1,j)+km(k,i,j+1)+km(k,i+1,j+1) ) &
703 i = iundef; j = iundef; k = iundef
712 call check( __line__, dens(k,i,j) )
713 call check( __line__, km(k,i,j) )
714 call check( __line__, s22_c(k,i,j) )
716 qflx_sgs_momy(k,i,j,
ydir) = dens(k,i,j) * ( - 2.0_rp * km(k,i,j) * s22_c(k,i,j) )
721 i = iundef; j = iundef; k = iundef
724 if ( atmos_phy_tb_d1980_implicit )
then
725 call calc_tend_momy( tend, &
727 gsqrt, j13g, j23g, j33g, mapf, &
735 ap = - dt * 0.25_rp * ( dens(ks ,i,j )*km(ks ,i,j ) &
736 + dens(ks+1,i,j )*km(ks+1,i,j ) &
737 + dens(ks ,i,j+1)*km(ks ,i,j+1) &
738 + dens(ks+1,i,j+1)*km(ks+1,i,j+1) ) &
739 * rfdz(ks) / gsqrt(ks,i,j,
i_xvw)
740 a(ks,i,j) = ap * rcdz(ks) / gsqrt(ks,i,j,
i_xvz)
742 b(ks,i,j) = - a(ks,i,j) + 0.5_rp * ( dens(ks,i,j)+dens(ks,i,j+1) )
744 c(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,
i_xvz)
745 ap = - dt * 0.25_rp * ( dens(k ,i,j )*km(k ,i,j ) &
746 + dens(k+1,i,j )*km(k+1,i,j ) &
747 + dens(k ,i,j+1)*km(k ,i,j+1) &
748 + dens(k+1,i,j+1)*km(k+1,i,j+1) ) &
749 * rfdz(k) / gsqrt(k,i,j,
i_xvw)
750 a(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,
i_xvz)
751 b(k,i,j) = - a(k,i,j) - c(k,i,j) + 0.5_rp * ( dens(k,i,j)+dens(k,i,j+1) )
754 c(ke,i,j) = ap * rcdz(ke) / gsqrt(ke,i,j,
i_xvz)
755 b(ke,i,j) = - c(ke,i,j) + 0.5_rp * ( dens(ke,i,j)+dens(ke,i,j+1) )
759 call matrix_solver_tridiagonal( ka, ks, ke, ia, iis, iie, ja, jjs, jje, &
760 a(:,:,:), b(:,:,:), c(:,:,:), tend(:,:,:), &
767 qflx_sgs_momy(k,i,j,
zdir) = qflx_sgs_momy(k,i,j,
zdir) &
768 - 0.25_rp * ( dens(k ,i,j )*km(k ,i,j ) &
769 + dens(k+1,i,j )*km(k+1,i,j ) &
770 + dens(k ,i,j+1)*km(k ,i,j+1) &
771 + dens(k+1,i,j+1)*km(k+1,i,j+1) ) &
772 * dt * ( tend2(k+1,i,j) - tend2(k,i,j) ) * rfdz(k) / gsqrt(k,i,j,
i_xvw)
782 if ( atmos_phy_tb_d1980_implicit )
then
789 ap = - dt * 0.25_rp * ( dens(ks,i,j)+dens(ks+1,i,j) ) &
790 * ( kh(ks,i,j)+kh(ks+1,i,j) ) &
791 * rfdz(ks) / gsqrt(ks,i,j,
i_xyw)
792 a(ks,i,j) = ap * rcdz(ks) / gsqrt(ks,i,j,
i_xyz)
794 b(ks,i,j) = - a(ks,i,j) + dens(ks,i,j)
796 c(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,
i_xyz)
797 ap = - dt * 0.25_rp * ( dens(k,i,j)+dens(k+1,i,j) ) &
798 * ( kh(k,i,j)+kh(k+1,i,j) ) &
799 * rfdz(k) / gsqrt(k,i,j,
i_xyw)
800 a(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,
i_xyz)
801 b(k,i,j) = - a(k,i,j) - c(k,i,j) + dens(k,i,j)
804 c(ke,i,j) = ap * rcdz(ke) / gsqrt(ke,i,j,
i_xyz)
805 b(ke,i,j) = - c(ke,i,j) + dens(ke,i,j)
812 call calc_flux_phi( &
814 dens, pott, kh, 1.0_rp, &
815 gsqrt, j13g, j23g, j33g, mapf, &
817 atmos_phy_tb_d1980_implicit, &
828 if ( iq ==
i_tke )
then
829 qflx_sgs_rhoq(:,:,:,:,iq) = 0.0_rp
839 call calc_flux_phi( &
840 qflx_sgs_rhoq(:,:,:,:,iq), &
841 dens, qtrc(:,:,:,iq), kh, 1.0_rp, &
842 gsqrt, j13g, j23g, j33g, mapf, &
844 atmos_phy_tb_d1980_implicit, &
851 iis = iundef; iie = iundef; jjs = iundef; jje = iundef
865 if ( atmos_phy_tb_d1980_implicit )
then
872 ap = - dt * 0.25_rp * ( dens(ks,i,j)+dens(ks+1,i,j) ) &
873 * 2.0_rp * ( km(ks,i,j)+km(ks+1,i,j) ) &
874 * rfdz(ks) / gsqrt(ks,i,j,
i_xyw)
875 a(ks,i,j) = ap * rcdz(ks) / gsqrt(ks,i,j,
i_xyz)
877 b(ks,i,j) = - a(ks,i,j) + dens(ks,i,j)
879 c(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,
i_xyz)
880 ap = - dt * 0.25_rp * ( dens(k,i,j)+dens(k+1,i,j) ) &
881 * 2.0_rp * ( km(k,i,j)+km(k+1,i,j) ) &
882 * rfdz(k) / gsqrt(k,i,j,
i_xyw)
883 a(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,
i_xyz)
884 b(k,i,j) = - a(k,i,j) - c(k,i,j) + dens(k,i,j)
887 c(ke,i,j) = ap * rcdz(ke) / gsqrt(ke,i,j,
i_xyz)
888 b(ke,i,j) = - c(ke,i,j) + dens(ke,i,j)
895 call calc_flux_phi( &
897 dens, qtrc(:,:,:,
i_tke), km, 2.0_rp, &
898 gsqrt, j13g, j23g, j33g, mapf, &
900 atmos_phy_tb_d1980_implicit, &
904 call calc_tend_phi( rhoq_t(:,:,:,
i_tke), &
906 gsqrt, j13g, j23g, j33g, mapf, &
912 rhoq_t(k,i,j,
i_tke) = max( &
913 rhoq_t(k,i,j,
i_tke) &
914 + ( km(k,i,j) * s2(k,i,j) &
915 - kh(k,i,j) * n2(k,i,j) &
916 - ( 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), &
917 - qtrc(k,i,j,
i_tke) * dens(k,i,j) / real(dt,kind=rp) )
922 i = iundef; j = iundef; k = iundef