73 real(RP),
private,
allocatable :: nu_fact (:,:,:)
75 real(RP),
private :: cs = 0.13_rp
76 real(RP),
private,
parameter :: ck = 0.1_rp
77 real(RP),
private,
parameter :: prn = 0.7_rp
78 real(RP),
private,
parameter :: ric = 0.25_rp
79 real(RP),
private,
parameter :: fmc = 16.0_rp
80 real(RP),
private,
parameter :: fhb = 40.0_rp
81 real(RP),
private :: rprn
82 real(RP),
private :: rric
83 real(RP),
private :: prnovric
85 real(RP),
private,
parameter :: oneoverthree = 1.0_rp / 3.0_rp
86 real(RP),
private,
parameter :: twooverthree = 2.0_rp / 3.0_rp
87 real(RP),
private,
parameter :: fouroverthree = 4.0_rp / 3.0_rp
89 real(RP),
private :: atmos_phy_tb_smg_nu_max = 10000.0_rp
90 logical,
private :: atmos_phy_tb_smg_implicit = .false.
91 logical,
private :: atmos_phy_tb_smg_bottom = .false.
93 real(RP),
private :: tke_fact
106 character(len=*),
intent(in) :: TYPE_TB
108 real(RP),
intent(in) :: CDZ(
ka)
109 real(RP),
intent(in) :: CDX(
ia)
110 real(RP),
intent(in) :: CDY(
ja)
111 real(RP),
intent(in) :: CZ (
ka,
ia,
ja)
113 real(RP) :: ATMOS_PHY_TB_SMG_Cs
114 real(RP) :: ATMOS_PHY_TB_SMG_filter_fact = 2.0_rp
115 logical :: ATMOS_PHY_TB_SMG_consistent_tke = .true.
117 namelist / param_atmos_phy_tb_smg / &
118 atmos_phy_tb_smg_cs, &
119 atmos_phy_tb_smg_nu_max, &
120 atmos_phy_tb_smg_filter_fact, &
121 atmos_phy_tb_smg_implicit, &
122 atmos_phy_tb_smg_consistent_tke, &
123 atmos_phy_tb_smg_bottom
131 if(
io_l )
write(
io_fid_log,*)
'++++++ Module[TURBULENCE] / Categ[ATMOS PHYSICS] / Origin[SCALElib]' 132 if(
io_l )
write(
io_fid_log,*)
'+++ Smagorinsky-type Eddy Viscocity Model' 134 if ( type_tb /=
'SMAGORINSKY' )
then 135 write(*,*)
'xxx ATMOS_PHY_TB_TYPE is not SMAGORINSKY. Check!' 139 atmos_phy_tb_smg_cs = cs
143 read(
io_fid_conf,nml=param_atmos_phy_tb_smg,iostat=ierr)
145 if(
io_l )
write(
io_fid_log,*)
'*** Not found namelist. Default used.' 146 elseif( ierr > 0 )
then 147 write(*,*)
'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_TB_SMG. Check!' 152 cs = atmos_phy_tb_smg_cs
156 prnovric = (1- prn) * rric
158 allocate( nu_fact(
ka,
ia,
ja) )
161 nu_fact(:,:,:) = undef
166 nu_fact(k,i,j) = ( cs *
mixlen(cdz(k),cdx(i),cdy(j),cz(k,i,j),atmos_phy_tb_smg_filter_fact) )**2
171 i = iundef; j = iundef; k = iundef
174 if ( atmos_phy_tb_smg_consistent_tke )
then 186 qflx_sgs_momz, qflx_sgs_momx, qflx_sgs_momy, &
187 qflx_sgs_rhot, qflx_sgs_rhoq, &
188 tke, tke_t, nu, Ri, Pr, N2, &
189 MOMZ, MOMX, MOMY, RHOT, DENS, QTRC, &
190 SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_QV, &
191 GSQRT, J13G, J23G, J33G, MAPF, dt )
230 real(RP),
intent(out) :: qflx_sgs_momz(
ka,
ia,
ja,3)
231 real(RP),
intent(out) :: qflx_sgs_momx(
ka,
ia,
ja,3)
232 real(RP),
intent(out) :: qflx_sgs_momy(
ka,
ia,
ja,3)
233 real(RP),
intent(out) :: qflx_sgs_rhot(
ka,
ia,
ja,3)
234 real(RP),
intent(out) :: qflx_sgs_rhoq(
ka,
ia,
ja,3,
qa)
236 real(RP),
intent(inout) :: tke(
ka,
ia,
ja)
237 real(RP),
intent(out) :: tke_t(
ka,
ia,
ja)
238 real(RP),
intent(out) :: nu (
ka,
ia,
ja)
239 real(RP),
intent(out) :: Ri (
ka,
ia,
ja)
240 real(RP),
intent(out) :: Pr (
ka,
ia,
ja)
241 real(RP),
intent(out) :: N2 (
ka,
ia,
ja)
243 real(RP),
intent(in) :: MOMZ(
ka,
ia,
ja)
244 real(RP),
intent(in) :: MOMX(
ka,
ia,
ja)
245 real(RP),
intent(in) :: MOMY(
ka,
ia,
ja)
246 real(RP),
intent(in) :: RHOT(
ka,
ia,
ja)
247 real(RP),
intent(in) :: DENS(
ka,
ia,
ja)
248 real(RP),
intent(in) :: QTRC(
ka,
ia,
ja,
qa)
250 real(RP),
intent(in) :: SFLX_MW(
ia,
ja)
251 real(RP),
intent(in) :: SFLX_MU(
ia,
ja)
252 real(RP),
intent(in) :: SFLX_MV(
ia,
ja)
253 real(RP),
intent(in) :: SFLX_SH(
ia,
ja)
254 real(RP),
intent(in) :: SFLX_QV(
ia,
ja)
256 real(RP),
intent(in) :: GSQRT (
ka,
ia,
ja,7)
257 real(RP),
intent(in) :: J13G (
ka,
ia,
ja,7)
258 real(RP),
intent(in) :: J23G (
ka,
ia,
ja,7)
259 real(RP),
intent(in) :: J33G
260 real(RP),
intent(in) :: MAPF(
ia,
ja,2,4)
261 real(DP),
intent(in) :: dt
264 real(RP) :: POTT (
ka,
ia,
ja)
267 real(RP) :: S33_C (
ka,
ia,
ja)
268 real(RP) :: S11_C (
ka,
ia,
ja)
269 real(RP) :: S22_C (
ka,
ia,
ja)
270 real(RP) :: S31_C (
ka,
ia,
ja)
271 real(RP) :: S12_C (
ka,
ia,
ja)
272 real(RP) :: S23_C (
ka,
ia,
ja)
273 real(RP) :: S12_Z (
ka,
ia,
ja)
274 real(RP) :: S23_X (
ka,
ia,
ja)
275 real(RP) :: S31_Y (
ka,
ia,
ja)
280 real(RP) :: TEND(
ka,
ia,
ja)
290 integer :: k, i, j, iq
293 if(
io_l )
write(
io_fid_log,*)
'*** Physics step: Turbulence(smagorinsky)' 296 qflx_sgs_momz(:,:,:,:) = undef
297 qflx_sgs_momx(:,:,:,:) = undef
298 qflx_sgs_momy(:,:,:,:) = undef
299 qflx_sgs_rhot(:,:,:,:) = undef
300 qflx_sgs_rhoq(:,:,:,:,:) = undef
319 call check( __line__, rhot(k,i,j) )
320 call check( __line__, dens(k,i,j) )
322 pott(k,i,j) = rhot(k,i,j) / dens(k,i,j)
327 i = iundef; j = iundef; k = iundef
332 call calc_strain_tensor( &
333 s33_c, s11_c, s22_c, &
334 s31_c, s12_c, s23_c, &
335 s12_z, s23_x, s31_y, &
337 dens, momz, momx, momy, &
338 gsqrt, j13g, j23g, j33g, mapf )
350 call check( __line__, pott(k+1,i,j) )
351 call check( __line__, pott(k,i,j) )
352 call check( __line__, pott(k-1,i,j) )
353 call check( __line__, fdz(k) )
354 call check( __line__, fdz(k-1) )
355 call check( __line__, s2(k,i,j) )
357 n2(k,i,j) = grav * ( pott(k+1,i,j) - pott(k-1,i,j) ) * j33g &
358 / ( ( fdz(k) + fdz(k-1) ) * gsqrt(k,i,j,
i_xyz) * pott(k,i,j) )
359 ri(k,i,j) = n2(k,i,j) / s2(k,i,j)
364 i = iundef; j = iundef; k = iundef
369 call check( __line__, pott(
ks+1,i,j) )
370 call check( __line__, pott(
ks,i,j) )
371 call check( __line__, rfdz(
ks) )
372 call check( __line__, s2(
ks,i,j) )
374 n2(
ks,i,j) = grav * ( pott(
ks+1,i,j) - pott(
ks,i,j) ) * j33g &
376 ri(
ks,i,j) = grav * ( pott(
ks+1,i,j) - pott(
ks,i,j) ) * j33g * rfdz(
ks) &
377 / ( gsqrt(
ks,i,j,
i_xyz) * pott(
ks,i,j) * s2(
ks,i,j) )
381 i = iundef; j = iundef; k = iundef
386 call check( __line__, pott(
ke,i,j) )
387 call check( __line__, pott(
ke-1,i,j) )
388 call check( __line__, rfdz(
ke-1) )
389 call check( __line__, s2(
ke,i,j) )
391 n2(
ke,i,j) = grav * ( pott(
ke,i,j) - pott(
ke-1,i,j) ) * j33g &
392 / ( fdz(
ke-1) * gsqrt(
ke,i,j,
i_xyz) * pott(
ke,i,j) )
393 ri(
ke,i,j) = grav * ( pott(
ke,i,j) - pott(
ke-1,i,j) ) * j33g * rfdz(
ke-1) &
394 / ( gsqrt(
ke,i,j,
i_xyz) * pott(
ke,i,j) * s2(
ke,i,j) )
398 i = iundef; j = iundef; k = iundef
404 call check( __line__, ri(k,i,j) )
405 call check( __line__, nu_fact(k,i,j) )
406 call check( __line__, s2(k,i,j) )
408 if ( ri(k,i,j) < 0.0_rp )
then 409 nu(k,i,j) = nu_fact(k,i,j) &
410 * sqrt( s2(k,i,j) * (1.0_rp - fmc*ri(k,i,j)) )
411 else if ( ri(k,i,j) < ric )
then 412 nu(k,i,j) = nu_fact(k,i,j) &
413 * sqrt( s2(k,i,j) ) * ( 1.0_rp - ri(k,i,j)*rric )**4
421 i = iundef; j = iundef; k = iundef
428 call check( __line__, ri(k,i,j) )
430 if ( ri(k,i,j) < 0.0_rp )
then 431 pr(k,i,j) = sqrt( ( 1.0_rp - fmc*ri(k,i,j) ) &
432 / ( 1.0_rp - fhb*ri(k,i,j) ) ) * prn
433 else if ( ri(k,i,j) < ric )
then 434 pr(k,i,j) = prn / ( 1.0_rp - prnovric * ri(k,i,j) )
438 kh(k,i,j) = min( nu(k,i,j)/pr(k,i,j), atmos_phy_tb_smg_nu_max )
439 nu(k,i,j) = min( nu(k,i,j), atmos_phy_tb_smg_nu_max )
440 pr(k,i,j) = nu(k,i,j) / ( kh(k,i,j) + ( 0.5_rp - sign(0.5_rp, kh(k,i,j)-eps) ) )
445 i = iundef; j = iundef; k = iundef
454 call check( __line__, nu(k,i,j) )
455 call check( __line__, nu_fact(k,i,j) )
457 tke(k,i,j) = ( nu(k,i,j) * cs / ck )**2 / nu_fact(k,i,j)
462 i = iundef; j = iundef; k = iundef
471 call check( __line__, dens(k,i,j) )
472 call check( __line__, nu(k,i,j) )
473 call check( __line__, s33_c(k,i,j) )
474 call check( __line__, s11_c(k,i,j) )
475 call check( __line__, s22_c(k,i,j) )
476 call check( __line__, tke(k,i,j) )
478 qflx_sgs_momz(k,i,j,
zdir) = dens(k,i,j) * ( &
479 - 2.0_rp * nu(k,i,j) &
480 * ( s33_c(k,i,j) - ( s11_c(k,i,j) + s22_c(k,i,j) + s33_c(k,i,j) ) * oneoverthree * tke_fact ) &
481 + twooverthree * tke(k,i,j) * tke_fact )
486 i = iundef; j = iundef; k = iundef
490 qflx_sgs_momz(
ks,i,j,
zdir) = 0.0_rp
491 qflx_sgs_momz(
ke,i,j,
zdir) = 0.0_rp
495 i = iundef; j = iundef; k = iundef
502 call check( __line__, dens(k,i,j) )
503 call check( __line__, dens(k,i+1,j) )
504 call check( __line__, dens(k+1,i,j) )
505 call check( __line__, dens(k+1,i+1,j) )
506 call check( __line__, nu(k,i,j) )
507 call check( __line__, nu(k,i+1,j) )
508 call check( __line__, nu(k+1,i,j) )
509 call check( __line__, nu(k+1,i+1,j) )
510 call check( __line__, s31_y(k,i,j) )
512 qflx_sgs_momz(k,i,j,
xdir) = - 0.125_rp &
513 * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i+1,j)+dens(k+1,i+1,j) ) &
514 * ( nu(k,i,j)+nu(k+1,i,j)+nu(k,i+1,j)+nu(k+1,i+1,j)) &
520 i = iundef; j = iundef; k = iundef
527 call check( __line__, dens(k,i,j) )
528 call check( __line__, dens(k,i,j+1) )
529 call check( __line__, dens(k+1,i,j) )
530 call check( __line__, dens(k+1,i,j+1) )
531 call check( __line__, nu(k,i,j) )
532 call check( __line__, nu(k,i,j+1) )
533 call check( __line__, nu(k+1,i,j) )
534 call check( __line__, nu(k+1,i,j+1) )
535 call check( __line__, s23_x(k,i,j) )
537 qflx_sgs_momz(k,i,j,
ydir) = - 0.125_rp &
538 * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i,j+1)+dens(k+1,i,j+1) ) &
539 * ( nu(k,i,j)+nu(k+1,i,j)+nu(k,i,j+1)+nu(k+1,i,j+1) ) &
545 i = iundef; j = iundef; k = iundef
548 if ( atmos_phy_tb_smg_implicit )
then 549 call calc_tend_momz( tend, &
551 gsqrt, j13g, j23g, j33g, mapf, &
557 ap = - fouroverthree * dt &
558 * dens(
ks+1,i,j)*nu(
ks+1,i,j) &
562 b(
ks,i,j) = - a(
ks,i,j) + 0.5_rp * ( dens(
ks,i,j)+dens(
ks+1,i,j) )
564 c(k,i,j) = ap * rfdz(k+1) / gsqrt(k+1,i,j,
i_xyw)
565 ap = - fouroverthree * dt &
566 * dens(k+1,i,j)*nu(k+1,i,j) &
567 * rcdz(k+1) / gsqrt(k+1,i,j,
i_xyz)
568 a(k,i,j) = ap * rfdz(k) / gsqrt(k,i,j,
i_xyw)
569 b(k,i,j) = - a(k,i,j) - c(k,i,j) + 0.5_rp * ( dens(k,i,j)+dens(k+1,i,j) )
573 b(
ke-1,i,j) = - c(
ke-1,i,j) + 0.5_rp * ( dens(
ke-1,i,j)+dens(
ke,i,j) )
579 call diffusion_solver( &
581 a(:,i,j), b(:,i,j), c(:,i,j), d, &
585 qflx_sgs_momz(k,i,j,
zdir) = qflx_sgs_momz(k,i,j,
zdir) &
586 - fouroverthree * dens(k,i,j) * nu(k,i,j) * dt &
587 * ( tend(k,i,j) - tend(k-1,i,j) ) * rcdz(k) / gsqrt(k,i,j,
i_xyz)
601 call check( __line__, dens(k,i,j) )
602 call check( __line__, dens(k,i+1,j) )
603 call check( __line__, dens(k+1,i,j) )
604 call check( __line__, dens(k+1,i+1,j) )
605 call check( __line__, nu(k,i,j) )
606 call check( __line__, nu(k,i+1,j) )
607 call check( __line__, nu(k+1,i,j) )
608 call check( __line__, nu(k+1,i+1,j) )
609 call check( __line__, s31_y(k,i,j) )
611 qflx_sgs_momx(k,i,j,
zdir) = - 0.125_rp &
612 * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i+1,j)+dens(k+1,i+1,j) ) &
613 * ( nu(k,i,j)+nu(k+1,i,j)+nu(k,i+1,j)+nu(k+1,i+1,j) ) &
619 i = iundef; j = iundef; k = iundef
623 qflx_sgs_momx(
ks-1,i,j,
zdir) = 0.0_rp
624 qflx_sgs_momx(
ke ,i,j,
zdir) = 0.0_rp
628 i = iundef; j = iundef; k = iundef
635 call check( __line__, dens(k,i,j) )
636 call check( __line__, nu(k,i,j) )
637 call check( __line__, s11_c(k,i,j) )
638 call check( __line__, s22_c(k,i,j) )
639 call check( __line__, s33_c(k,i,j) )
640 call check( __line__, tke(k,i,j) )
642 qflx_sgs_momx(k,i,j,
xdir) = dens(k,i,j) * ( &
643 - 2.0_rp * nu(k,i,j) &
644 * ( s11_c(k,i,j) - ( s11_c(k,i,j) + s22_c(k,i,j) + s33_c(k,i,j) ) * oneoverthree * tke_fact ) &
645 + twooverthree * tke(k,i,j) * tke_fact )
650 i = iundef; j = iundef; k = iundef
657 call check( __line__, dens(k,i,j) )
658 call check( __line__, dens(k,i+1,j) )
659 call check( __line__, dens(k,i,j+1) )
660 call check( __line__, dens(k,i+1,j+1) )
661 call check( __line__, nu(k,i,j) )
662 call check( __line__, nu(k,i+1,j) )
663 call check( __line__, nu(k,i,j+1) )
664 call check( __line__, nu(k,i+1,j+1) )
665 call check( __line__, s12_z(k,i,j) )
667 qflx_sgs_momx(k,i,j,
ydir) = - 0.125_rp &
668 * ( dens(k,i,j)+dens(k,i+1,j)+dens(k,i,j+1)+dens(k,i+1,j+1) ) &
669 * ( nu(k,i,j)+nu(k,i+1,j)+nu(k,i,j+1)+nu(k,i+1,j+1) ) &
675 i = iundef; j = iundef; k = iundef
678 if ( atmos_phy_tb_smg_implicit )
then 679 call calc_tend_momx( tend, &
681 gsqrt, j13g, j23g, j33g, mapf, &
687 ap = - dt * 0.25_rp * ( dens(
ks ,i ,j)*nu(
ks ,i ,j) &
688 + dens(
ks+1,i ,j)*nu(
ks+1,i ,j) &
689 + dens(
ks ,i+1,j)*nu(
ks ,i+1,j) &
690 + dens(
ks+1,i+1,j)*nu(
ks+1,i+1,j) ) &
694 b(
ks,i,j) = - a(
ks,i,j) + 0.5_rp * ( dens(
ks,i,j)+dens(
ks,i+1,j) )
696 c(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,
i_uyz)
697 ap = - dt * 0.25_rp * ( dens(k ,i ,j)*nu(k ,i ,j) &
698 + dens(k+1,i ,j)*nu(k+1,i ,j) &
699 + dens(k ,i+1,j)*nu(k ,i+1,j) &
700 + dens(k+1,i+1,j)*nu(k+1,i+1,j) ) &
701 * rfdz(k) / gsqrt(k,i,j,
i_uyw)
702 a(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,
i_uyz)
703 b(k,i,j) = - a(k,i,j) - c(k,i,j) + 0.5_rp * ( dens(k,i,j)+dens(k,i+1,j) )
707 b(
ke,i,j) = - c(
ke,i,j) + 0.5_rp * ( dens(
ke,i,j)+dens(
ke,i+1,j) )
713 call diffusion_solver( &
715 a(:,i,j), b(:,i,j), c(:,i,j), d, &
719 qflx_sgs_momx(k,i,j,
zdir) = qflx_sgs_momx(k,i,j,
zdir) &
720 - 0.25_rp * ( dens(k ,i ,j)*nu(k ,i ,j) &
721 + dens(k+1,i ,j)*nu(k+1,i ,j) &
722 + dens(k ,i+1,j)*nu(k ,i+1,j) &
723 + dens(k+1,i+1,j)*nu(k+1,i+1,j) ) &
724 * dt * ( tend(k+1,i,j) - tend(k,i,j) ) * rfdz(k) / gsqrt(k,i,j,
i_uyw)
738 call check( __line__, dens(k,i,j) )
739 call check( __line__, dens(k,i,j+1) )
740 call check( __line__, dens(k+1,i,j) )
741 call check( __line__, dens(k+1,i,j+1) )
742 call check( __line__, nu(k,i,j) )
743 call check( __line__, nu(k,i,j+1) )
744 call check( __line__, nu(k+1,i,j) )
745 call check( __line__, nu(k+1,i,j+1) )
746 call check( __line__, s23_x(k,i,j) )
748 qflx_sgs_momy(k,i,j,
zdir) = - 0.125_rp &
749 * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i,j+1)+dens(k+1,i,j+1) ) &
750 * ( nu(k,i,j)+nu(k+1,i,j)+nu(k,i,j+1)+nu(k+1,i,j+1) ) &
756 i = iundef; j = iundef; k = iundef
760 qflx_sgs_momy(
ks-1,i,j,
zdir) = 0.0_rp
761 qflx_sgs_momy(
ke ,i,j,
zdir) = 0.0_rp
765 i = iundef; j = iundef; k = iundef
773 call check( __line__, dens(k,i,j) )
774 call check( __line__, dens(k,i+1,j) )
775 call check( __line__, dens(k,i,j+1) )
776 call check( __line__, dens(k,i+1,j+1) )
777 call check( __line__, nu(k,i,j) )
778 call check( __line__, nu(k,i+1,j) )
779 call check( __line__, nu(k,i,j+1) )
780 call check( __line__, nu(k,i+1,j+1) )
781 call check( __line__, s12_z(k,i,j) )
783 qflx_sgs_momy(k,i,j,
xdir) = - 0.125_rp &
784 * ( dens(k,i,j)+dens(k,i+1,j)+dens(k,i,j+1)+dens(k,i+1,j+1) ) &
785 * ( nu(k,i,j)+nu(k,i+1,j)+nu(k,i,j+1)+nu(k,i+1,j+1) ) &
791 i = iundef; j = iundef; k = iundef
799 call check( __line__, dens(k,i,j) )
800 call check( __line__, nu(k,i,j) )
801 call check( __line__, s11_c(k,i,j) )
802 call check( __line__, s22_c(k,i,j) )
803 call check( __line__, s33_c(k,i,j) )
804 call check( __line__, tke(k,i,j) )
806 qflx_sgs_momy(k,i,j,
ydir) = dens(k,i,j) * ( &
807 - 2.0_rp * nu(k,i,j) &
808 * ( s22_c(k,i,j) - ( s11_c(k,i,j) + s22_c(k,i,j) + s33_c(k,i,j) ) * oneoverthree * tke_fact ) &
809 + twooverthree * tke(k,i,j) * tke_fact)
814 i = iundef; j = iundef; k = iundef
817 if ( atmos_phy_tb_smg_implicit )
then 818 call calc_tend_momy( tend, &
820 gsqrt, j13g, j23g, j33g, mapf, &
826 ap = - dt * 0.25_rp * ( dens(
ks ,i,j )*nu(
ks ,i,j ) &
827 + dens(
ks+1,i,j )*nu(
ks+1,i,j ) &
828 + dens(
ks ,i,j+1)*nu(
ks ,i,j+1) &
829 + dens(
ks+1,i,j+1)*nu(
ks+1,i,j+1) ) &
833 b(
ks,i,j) = - a(
ks,i,j) + 0.5_rp * ( dens(
ks,i,j)+dens(
ks,i,j+1) )
835 c(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,
i_xvz)
836 ap = - dt * 0.25_rp * ( dens(k ,i,j )*nu(k ,i,j ) &
837 + dens(k+1,i,j )*nu(k+1,i,j ) &
838 + dens(k ,i,j+1)*nu(k ,i,j+1) &
839 + dens(k+1,i,j+1)*nu(k+1,i,j+1) ) &
840 * rfdz(k) / gsqrt(k,i,j,
i_xvw)
841 a(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,
i_xvz)
842 b(k,i,j) = - a(k,i,j) - c(k,i,j) + 0.5_rp * ( dens(k,i,j)+dens(k,i,j+1) )
846 b(
ke,i,j) = - c(
ke,i,j) + 0.5_rp * ( dens(
ke,i,j)+dens(
ke,i,j+1) )
852 call diffusion_solver( &
854 a(:,i,j), b(:,i,j), c(:,i,j), d, &
858 qflx_sgs_momy(k,i,j,
zdir) = qflx_sgs_momy(k,i,j,
zdir) &
859 - 0.25_rp * ( dens(k ,i,j )*nu(k ,i,j ) &
860 + dens(k+1,i,j )*nu(k+1,i,j ) &
861 + dens(k ,i,j+1)*nu(k ,i,j+1) &
862 + dens(k+1,i,j+1)*nu(k+1,i,j+1) ) &
863 * dt * ( tend(k+1,i,j) - tend(k,i,j) ) * rfdz(k) / gsqrt(k,i,j,
i_xvw)
873 if ( atmos_phy_tb_smg_implicit )
then 878 ap = - dt * 0.25_rp * ( dens(
ks,i,j)+dens(
ks+1,i,j) ) &
879 * ( kh(
ks,i,j)+kh(
ks+1,i,j) ) &
883 b(
ks,i,j) = - a(
ks,i,j) + dens(
ks,i,j)
885 c(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,
i_xyz)
886 ap = - dt * 0.25_rp * ( dens(k,i,j)+dens(k+1,i,j) ) &
887 * ( kh(k,i,j)+kh(k+1,i,j) ) &
888 * rfdz(k) / gsqrt(k,i,j,
i_xyw)
889 a(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,
i_xyz)
890 b(k,i,j) = - a(k,i,j) - c(k,i,j) + dens(k,i,j)
894 b(
ke,i,j) = - c(
ke,i,j) + dens(
ke,i,j)
901 call calc_flux_phi( &
903 dens, pott, kh, 1.0_rp, &
904 gsqrt, j13g, j23g, j33g, mapf, &
906 atmos_phy_tb_smg_implicit, &
921 call calc_flux_phi( &
922 qflx_sgs_rhoq(:,:,:,:,iq), &
923 dens, qtrc(:,:,:,iq), kh, 1.0_rp, &
924 gsqrt, j13g, j23g, j33g, mapf, &
926 atmos_phy_tb_smg_implicit, &
933 iis = iundef; iie = iundef; jjs = iundef; jje = iundef
947 function mixlen(dz, dx, dy, z, filter_fact)
951 real(RP),
intent(in) :: dz
952 real(RP),
intent(in) :: dx
953 real(RP),
intent(in) :: dy
954 real(RP),
intent(in) :: z
955 real(RP),
intent(in) :: filter_fact
960 d0 =
fact(dz, dx, dy) * filter_fact * ( dz * dx * dy )**oneoverthree
961 if ( atmos_phy_tb_smg_bottom )
then 962 mixlen = sqrt( 1.0_rp / ( 1.0_rp/d0**2 + 1.0_rp/(karman*z)**2 ) )
970 function fact(dz, dx, dy)
971 real(RP),
intent(in) :: dz
972 real(RP),
intent(in) :: dx
973 real(RP),
intent(in) :: dy
976 real(RP),
parameter :: oot = -1.0_rp/3.0_rp
977 real(RP),
parameter :: fot = 5.0_rp/3.0_rp
978 real(RP),
parameter :: eot = 11.0_rp/3.0_rp
979 real(RP),
parameter :: tof = -3.0_rp/4.0_rp
980 real(RP) :: a1, a2, b1, b2, dmax
983 dmax = max(dz, dx, dy)
984 if ( dz == dmax )
then 987 else if ( dx == dmax )
then 997 fact = 1.736_rp * (a1*a2)**oot &
998 * ( 4.0_rp*p1(b1)*a1**oot + 0.222_rp*p2(b1)*a1**fot + 0.077*p3(b1)*a1**eot - 3.0_rp*b1 &
999 + 4.0_rp*p1(b2)*a2**oot + 0.222_rp*p2(b2)*a2**fot + 0.077*p3(b2)*a2**eot - 3.0_rp*b2 &
1004 real(RP),
intent(in) :: z
1007 p1 = 2.5_rp * p2(z) - 1.5_rp * sin(z) * cos(z)**twooverthree
1011 real(RP),
intent(in) :: z
1014 p2 = 0.986_rp * z + 0.073_rp * z**2 - 0.418_rp * z**3 + 0.120_rp * z**4
1018 real(RP),
intent(in) :: z
1021 p3 = 0.976_rp * z + 0.188_rp * z**2 - 1.169_rp * z**3 + 0.755_rp * z**4 - 0.151_rp * z**5
integer, public is
start point of inner domain: x, local
subroutine, public atmos_phy_tb_calc_tend_momz(MOMZ_t_TB, QFLX_MOMZ, GSQRT, J13G, J23G, J33G, MAPF, IIS, IIE, JJS, JJE)
integer, public je
end point of inner domain: y, local
real(rp), dimension(:), allocatable, public grid_rcdy
reciprocal of center-dy
subroutine, public prc_mpistop
Abort MPI.
subroutine, public atmos_phy_tb_diffusion_solver(phi, a, b, c, d, KE_TB)
subroutine, public atmos_phy_tb_calc_flux_phi(qflx_phi, DENS, PHI, Kh, FACT, GSQRT, J13G, J23G, J33G, MAPF, a, b, c, dt, implicit, IIS, IIE, JJS, JJE)
real(rp), dimension(:), allocatable, public grid_fdy
y-length of grid(j+1) to grid(j) [m]
integer, public iblock
block size for cache blocking: x
subroutine, public atmos_phy_tb_calc_tend_momy(MOMY_t_TB, QFLX_MOMY, GSQRT, J13G, J23G, J33G, MAPF, IIS, IIE, JJS, JJE)
logical, public io_l
output log or not? (this process)
real(rp), dimension(:), allocatable, public grid_rcdx
reciprocal of center-dx
integer, parameter, public zdir
real(rp) function fact(dz, dx, dy)
integer, parameter, public ydir
integer, public ke
end point of inner domain: z, local
integer, parameter, public xdir
real(rp), dimension(:), allocatable, public grid_rfdy
reciprocal of face-dy
subroutine, public check(current_line, v)
Undefined value checker.
real(rp), dimension(:), allocatable, public grid_rcdz
reciprocal of center-dz
real(rp), parameter, public const_karman
von Karman constant
real(rp), public const_undef
real(rp) function mixlen(dz, dx, dy, z, filter_fact)
subroutine, public atmos_phy_tb_calc_tend_phi(phi_t_TB, QFLX_phi, GSQRT, J13G, J23G, J33G, MAPF, IIS, IIE, JJS, JJE)
module ATMOSPHERE / Physics Turbulence
integer, public ia
of x whole cells (local, with HALO)
subroutine, public atmos_phy_tb_calc_tend_momx(MOMX_t_TB, QFLX_MOMX, GSQRT, J13G, J23G, J33G, MAPF, IIS, IIE, JJS, JJE)
real(rp), dimension(:), allocatable, public grid_fdz
z-length of grid(k+1) to grid(k) [m]
integer, public ka
of z whole cells (local, with HALO)
integer, public jblock
block size for cache blocking: y
subroutine, public atmos_phy_tb_smg(qflx_sgs_momz, qflx_sgs_momx, qflx_sgs_momy, qflx_sgs_rhot, qflx_sgs_rhoq, tke, tke_t, nu, Ri, Pr, N2, MOMZ, MOMX, MOMY, RHOT, DENS, QTRC, SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_QV, GSQRT, J13G, J23G, J33G, MAPF, dt)
real(rp), public const_grav
standard acceleration of gravity [m/s2]
integer, public js
start point of inner domain: y, local
subroutine, public atmos_phy_tb_calc_strain_tensor(S33_C, S11_C, S22_C, S31_C, S12_C, S23_C, S12_Z, S23_X, S31_Y, S2, DENS, MOMZ, MOMX, MOMY, GSQRT, J13G, J23G, J33G, MAPF)
integer, parameter, public const_undef2
undefined value (INT2)
subroutine, public atmos_phy_tb_smg_setup(TYPE_TB, CDZ, CDX, CDY, CZ)
integer, public ks
start point of inner domain: z, local
module ATMOSPHERE / Physics Turbulence
integer, public ie
end point of inner domain: x, local
real(rp), public const_eps
small number
logical, public io_lnml
output log or not? (for namelist, this process)
real(rp), dimension(:), allocatable, public grid_fdx
x-length of grid(i+1) to grid(i) [m]
real(rp), dimension(:), allocatable, public grid_rfdx
reciprocal of face-dx
real(rp), dimension(:), allocatable, public grid_rfdz
reciprocal of face-dz
integer, public io_fid_conf
Config file ID.
integer, public io_fid_log
Log file ID.
integer, public ja
of y whole cells (local, with HALO)