This subroutine calculates shear-driven diffusivity and TKE in a single column.
664 type(verticalGrid_type),
intent(in) :: GV
665 type(unit_scale_type),
intent(in) :: US
666 real,
dimension(SZK_(GV)+1), &
667 intent(inout) :: kappa
668 real,
dimension(SZK_(GV)+1), &
671 integer,
intent(in) :: nzc
672 real,
intent(in) :: f2
673 real,
intent(in) :: surface_pres
674 real,
dimension(SZK_(GV)), &
676 real,
dimension(SZK_(GV)), &
678 real,
dimension(SZK_(GV)), &
680 real,
dimension(SZK_(GV)), &
682 real,
dimension(SZK_(GV)), &
684 real,
dimension(SZK_(GV)+1), &
685 intent(out) :: kappa_avg
686 real,
dimension(SZK_(GV)+1), &
687 intent(out) :: tke_avg
688 real,
intent(in) :: dt
689 type(thermo_var_ptrs),
intent(in) :: tv
692 type(Kappa_shear_CS),
pointer :: CS
694 real,
dimension(SZK_(GV)+1), &
695 optional,
intent(out) :: I_Ld2_1d
696 real,
dimension(SZK_(GV)+1), &
697 optional,
intent(out) :: dz_Int_1d
700 real,
dimension(nzc) :: &
709 real,
dimension(nzc+1) :: &
744 real :: dist_from_bot
752 real :: tol_dksrc, tol2
753 real :: tol_dksrc_low
768 logical :: use_temperature
770 integer :: ks_kappa, ke_kappa
771 integer :: dt_halvings
774 integer :: dt_refinements
777 integer :: k, itt, itt_dt
779 integer :: max_debug_itt ;
parameter(max_debug_itt=20)
780 real :: wt(SZK_(GV)+1), wt_tot, I_wt_tot, wt_itt
781 real,
dimension(SZK_(GV)+1) :: &
782 Ri_k, tke_prev, dtke, dkappa, dtke_norm, &
785 real,
dimension(SZK_(GV)+1,0:max_debug_itt) :: &
786 tke_it1, N2_it1, Sh2_it1, ksrc_it1, kappa_it1, kprev_it1
787 real,
dimension(SZK_(GV)+1,1:max_debug_itt) :: &
788 dkappa_it1, wt_it1, K_Q_it1, d_dkappa_it1, dkappa_norm
789 real,
dimension(SZK_(GV),0:max_debug_itt) :: &
790 u_it1, v_it1, rho_it1, T_it1, S_it1
791 real,
dimension(0:max_debug_itt) :: &
792 dk_wt_it1, dkpos_wt_it1, dkneg_wt_it1, k_mag
793 real,
dimension(max_debug_itt) :: dt_it1
796 ri_crit = cs%Rino_crit
797 gr0 = gv%z_to_H*gv%H_to_Pa
798 g_r0 = (us%L_to_Z**2 * gv%g_Earth) / gv%Rho0
802 tol_dksrc = 10.0 ; tol_dksrc_low = 0.95 ; tol2 = 2.0*cs%kappa_tol_err
804 use_temperature = .false. ;
if (
associated(tv%T)) use_temperature = .true.
808 do k=1,nzc ; idz(k) = 1.0 / dz(k) ;
enddo
811 i_dz_int(1) = 2.0 / dz(1)
812 dist_from_top(1) = 0.0
814 i_dz_int(k) = 2.0 / (dz(k-1) + dz(k))
815 dist_from_top(k) = dist_from_top(k-1) + dz(k-1)
817 i_dz_int(nzc+1) = 2.0 / dz(nzc)
822 a1(2) = k0dt*i_dz_int(2)
823 b1 = 1.0 / (dz(1) + a1(2))
824 u(1) = b1 * u0xdz(1) ; v(1) = b1 * v0xdz(1)
825 t(1) = b1 * t0xdz(1) ; sal(1) = b1 * s0xdz(1)
826 c1(2) = a1(2) * b1 ; d1 = dz(1) * b1
828 bd1 = dz(k) + d1*a1(k)
829 a1(k+1) = k0dt*i_dz_int(k+1)
830 b1 = 1.0 / (bd1 + a1(k+1))
831 u(k) = b1 * (u0xdz(k) + a1(k)*u(k-1))
832 v(k) = b1 * (v0xdz(k) + a1(k)*v(k-1))
833 t(k) = b1 * (t0xdz(k) + a1(k)*t(k-1))
834 sal(k) = b1 * (s0xdz(k) + a1(k)*sal(k-1))
835 c1(k+1) = a1(k+1) * b1 ; d1 = bd1 * b1
840 b1 = 1.0 / ((dz(nzc) + d1*a1(nzc)) + k0dt*i_dz_int(nzc+1))
841 u(nzc) = b1 * (u0xdz(nzc) + a1(nzc)*u(nzc-1))
842 v(nzc) = b1 * (v0xdz(nzc) + a1(nzc)*v(nzc-1))
844 b1 = 1.0 / (dz(nzc) + d1*a1(nzc))
845 t(nzc) = b1 * (t0xdz(nzc) + a1(nzc)*t(nzc-1))
846 sal(nzc) = b1 * (s0xdz(nzc) + a1(nzc)*sal(nzc-1))
848 u(k) = u(k) + c1(k+1)*u(k+1) ; v(k) = v(k) + c1(k+1)*v(k+1)
849 t(k) = t(k) + c1(k+1)*t(k+1) ; sal(k) = sal(k) + c1(k+1)*sal(k+1)
853 b1 = 1.0 / (dz(1) + k0dt*i_dz_int(2))
854 u(1) = b1 * u0xdz(1) ; v(1) = b1 * v0xdz(1)
856 t(1) = b1 * t0xdz(1) ; sal(1) = b1 * s0xdz(1)
865 dz_int(1) = 0.0 ; dz_int(2) = dz(1)
867 norm = 1.0 / (dz(k)*(dz(k-1)+dz(k+1)) + 2.0*dz(k-1)*dz(k+1))
868 dz_int(k) = dz_int(k) + dz(k) * ( ((dz(k)+dz(k+1)) * dz(k-1)) * norm)
869 dz_int(k+1) = dz(k) * ( ((dz(k-1)+dz(k)) * dz(k+1)) * norm)
871 dz_int(nzc) = dz_int(nzc) + dz(nzc) ; dz_int(nzc+1) = 0.0
875 dist_from_bot = dist_from_bot + dz(k)
876 i_l2_bdry(k) = (dist_from_top(k) + dist_from_bot)**2 / &
877 (dist_from_top(k) * dist_from_bot)**2
881 if (use_temperature)
then
882 pressure(1) = surface_pres
884 pressure(k) = pressure(k-1) + gr0*dz(k-1)
885 t_int(k) = 0.5*(t(k-1) + t(k))
886 sal_int(k) = 0.5*(sal(k-1) + sal(k))
888 call calculate_density_derivs(t_int, sal_int, pressure, dbuoy_dt, &
889 dbuoy_ds, 2, nzc-1, tv%eqn_of_state)
891 dbuoy_dt(k) = -g_r0*dbuoy_dt(k)
892 dbuoy_ds(k) = -g_r0*dbuoy_ds(k)
895 do k=1,nzc+1 ; dbuoy_dt(k) = -g_r0 ; dbuoy_ds(k) = 0.0 ;
enddo
899 n2_debug(1) = 0.0 ; n2_debug(nzc+1) = 0.0
901 n2_debug(k) = max((dbuoy_dt(k) * (t0xdz(k-1)*idz(k-1) - t0xdz(k)*idz(k)) + &
902 dbuoy_ds(k) * (s0xdz(k-1)*idz(k-1) - s0xdz(k)*idz(k))) * &
906 u_it1(k,0) = u0xdz(k)*idz(k) ; v_it1(k,0) = v0xdz(k)*idz(k)
907 t_it1(k,0) = t0xdz(k)*idz(k) ; s_it1(k,0) = s0xdz(k)*idz(k)
910 kprev_it1(k,0) = kappa(k) ; kappa_it1(k,0) = kappa(k)
912 n2_it1(k,0) = n2_debug(k) ; sh2_it1(k,0) = s2(k) ; ksrc_it1(k,0) = k_src(k)
915 u_it1(k,0) = 0.0 ; v_it1(k,0) = 0.0
916 t_it1(k,0) = 0.0 ; s_it1(k,0) = 0.0
917 kprev_it1(k+1,0) = 0.0 ; kappa_it1(k+1,0) = 0.0 ; tke_it1(k+1,0) = 0.0
918 n2_it1(k+1,0) = 0.0 ; sh2_it1(k+1,0) = 0.0 ; ksrc_it1(k+1,0) = 0.0
920 do itt=1,max_debug_itt
923 u_it1(k,itt) = 0.0 ; v_it1(k,itt) = 0.0
924 t_it1(k,itt) = 0.0 ; s_it1(k,itt) = 0.0
928 kprev_it1(k,itt) = 0.0 ; kappa_it1(k,itt) = 0.0 ; tke_it1(k,itt) = 0.0
929 n2_it1(k,itt) = 0.0 ; sh2_it1(k,itt) = 0.0
930 ksrc_it1(k,itt) = 0.0
931 dkappa_it1(k,itt) = 0.0 ; wt_it1(k,itt) = 0.0
932 k_q_it1(k,itt) = 0.0 ; d_dkappa_it1(k,itt) = 0.0
935 do k=1,gv%ke+1 ; ksrc_av(k) = 0.0 ;
enddo
939 call calculate_projected_state(kappa, u, v, t, sal, 0.0, nzc, dz, i_dz_int, &
940 dbuoy_dt, dbuoy_ds, u, v, t, sal, gv, us, &
941 n2=n2, s2=s2, vel_underflow=cs%vel_underflow)
948 kappa_avg(k) = 0.0 ; tke_avg(k) = 0.0
949 local_src_avg(k) = 0.0
951 if ( dz_int(k) > 0.0 ) &
952 local_src_avg(k) = 0.1 * k0dt * i_dz_int(k) / dz_int(k)
958 do itt=1,cs%max_KS_it
965 ri_k(k) = 1e3 ;
if (s2(k) > 1e-3*n2(k)) ri_k(k) = n2(k) / s2(k)
966 if (itt > 1)
then ; tke_prev(k) = tke(k) ;
else ; tke_prev(k) = 0.0 ;
endif
971 call find_kappa_tke(n2, s2, kappa, idz, dz_int, i_l2_bdry, f2, &
972 nzc, cs, gv, us, k_q, tke, kappa_out, kappa_src, local_src)
977 ks_kappa = gv%ke+1 ; ke_kappa = 0
978 do k=2,nzc ;
if (kappa_out(k) > 0.0)
then
981 do k=nzc,ks_kappa,-1 ;
if (kappa_out(k) > 0.0)
then
984 if (ke_kappa == nzc) kappa_out(nzc+1) = 0.0
990 if ((ke_kappa < ks_kappa) .or. (itt==cs%max_RiNo_it))
then
996 tol_max(k) = kappa_src(k) + tol_dksrc * local_src(k)
997 tol_min(k) = kappa_src(k) - tol_dksrc_low * local_src(k)
998 tol_chg(k) = tol2 * local_src_avg(k)
1001 do itt_dt=1,(cs%max_KS_it+1-itt)/2
1009 call calculate_projected_state(kappa_out, u, v, t, sal, 0.5*dt_test, nzc, dz, i_dz_int, &
1010 dbuoy_dt, dbuoy_ds, u_test, v_test, t_test, s_test, &
1011 gv, us, n2, s2, ks_int=ks_kappa, ke_int=ke_kappa, &
1012 vel_underflow=cs%vel_underflow)
1014 idtt = 1.0 / dt_test
1015 do k=max(ks_kappa-1,2),min(ke_kappa+1,nzc)
1016 if (n2(k) < ri_crit * s2(k))
then
1017 k_src(k) = (2.0 * cs%Shearmix_rate * sqrt(s2(k))) * &
1018 ((ri_crit*s2(k) - n2(k)) / (ri_crit*s2(k) + cs%FRi_curvature*n2(k)))
1019 if ((k_src(k) > max(tol_max(k), kappa_src(k) + idtt*tol_chg(k))) .or. &
1020 (k_src(k) < min(tol_min(k), kappa_src(k) - idtt*tol_chg(k))))
then
1021 valid_dt = .false. ;
exit
1024 if (0.0 < min(tol_min(k), kappa_src(k) - idtt*tol_chg(k)))
then
1025 valid_dt = .false. ; k_src(k) = 0.0 ;
exit
1031 dt_test = 0.5*dt_test
1033 if ((dt_test < dt_rem) .and. valid_dt)
then
1034 dt_inc = 0.5*dt_test
1035 do itt_dt=1,dt_refinements
1036 call calculate_projected_state(kappa_out, u, v, t, sal, 0.5*(dt_test+dt_inc), &
1037 nzc, dz, i_dz_int, dbuoy_dt, dbuoy_ds, u_test, v_test, t_test, s_test, &
1038 gv, us, n2, s2, ks_int=ks_kappa, ke_int=ke_kappa, vel_underflow=cs%vel_underflow)
1040 idtt = 1.0 / (dt_test+dt_inc)
1041 do k=max(ks_kappa-1,2),min(ke_kappa+1,nzc)
1042 if (n2(k) < ri_crit * s2(k))
then
1043 k_src(k) = (2.0 * cs%Shearmix_rate * sqrt(s2(k))) * &
1044 ((ri_crit*s2(k) - n2(k)) / (ri_crit*s2(k) + cs%FRi_curvature*n2(k)))
1045 if ((k_src(k) > max(tol_max(k), kappa_src(k) + idtt*tol_chg(k))) .or. &
1046 (k_src(k) < min(tol_min(k), kappa_src(k) - idtt*tol_chg(k))))
then
1047 valid_dt = .false. ;
exit
1050 if (0.0 < min(tol_min(k), kappa_src(k) - idtt*tol_chg(k)))
then
1051 valid_dt = .false. ; k_src(k) = 0.0 ;
exit
1056 if (valid_dt) dt_test = dt_test + dt_inc
1063 dt_now = min(dt_test*(1.0+cs%kappa_tol_err)+dt_inc, dt_rem)
1065 local_src_avg(k) = local_src_avg(k) + dt_now * local_src(k)
1072 if (ke_kappa < ks_kappa)
then
1075 dt_wt = dt_rem / dt ; dt_rem = 0.0
1080 tke_avg(k) = tke_avg(k) + dt_wt*tke(k)
1082 tke_pred(k) = tke(k) ; kappa_pred(k) = 0.0 ; kappa(k) = 0.0
1088 call calculate_projected_state(kappa_out, u, v, t, sal, dt_now, nzc, dz, i_dz_int, &
1089 dbuoy_dt, dbuoy_ds, u_test, v_test, t_test, s_test, &
1090 gv, us, n2=n2, s2=s2, ks_int=ks_kappa, ke_int=ke_kappa, &
1091 vel_underflow=cs%vel_underflow)
1095 do k=1,nzc+1 ; k_q_tmp(k) = k_q(k) ;
enddo
1096 call find_kappa_tke(n2, s2, kappa_out, idz, dz_int, i_l2_bdry, f2, &
1097 nzc, cs, gv, us, k_q_tmp, tke_pred, kappa_pred)
1100 ks_kappa = gv%ke+1 ; ke_kappa = 0
1102 kappa_mid(k) = 0.5*(kappa_out(k) + kappa_pred(k))
1103 if ((kappa_mid(k) > 0.0) .and. (k<ks_kappa)) ks_kappa = k
1104 if (kappa_mid(k) > 0.0) ke_kappa = k
1108 call calculate_projected_state(kappa_mid, u, v, t, sal, dt_now, nzc, dz, i_dz_int, &
1109 dbuoy_dt, dbuoy_ds, u_test, v_test, t_test, s_test, &
1110 gv, us, n2=n2, s2=s2, ks_int=ks_kappa, ke_int=ke_kappa, &
1111 vel_underflow=cs%vel_underflow)
1115 call find_kappa_tke(n2, s2, kappa_out, idz, dz_int, i_l2_bdry, f2, &
1116 nzc, cs, gv, us, k_q, tke_pred, kappa_pred)
1120 dt_wt = dt_now / dt ; dt_rem = dt_rem - dt_now
1122 kappa_mid(k) = 0.5*(kappa_out(k) + kappa_pred(k))
1123 kappa_avg(k) = kappa_avg(k) + kappa_mid(k)*dt_wt
1124 tke_avg(k) = tke_avg(k) + dt_wt*0.5*(tke_pred(k) + tke(k))
1125 kappa(k) = kappa_pred(k)
1130 if (dt_rem > 0.0)
then
1133 call calculate_projected_state(kappa_mid, u, v, t, sal, dt_now, nzc, &
1134 dz, i_dz_int, dbuoy_dt, dbuoy_ds, u, v, t, sal, &
1135 gv, us, n2, s2, vel_underflow=cs%vel_underflow)
1140 if (itt <= max_debug_itt)
then
1141 dt_it1(itt) = dt_now
1142 dk_wt_it1(itt) = 0.0 ; dkpos_wt_it1(itt) = 0.0 ; dkneg_wt_it1(itt) = 0.0
1144 wt_itt = 1.0/real(itt) ; wt_tot = 0.0
1146 ksrc_av(k) = (1.0-wt_itt)*ksrc_av(k) + wt_itt*k_src(k)
1147 wt_tot = wt_tot + dz_int(k) * ksrc_av(k)
1150 i_wt_tot = 0.0 ;
if (wt_tot > 0.0) i_wt_tot = 1.0/wt_tot
1153 wt(k) = (dz_int(k)*ksrc_av(k)) * i_wt_tot
1154 k_mag(itt) = k_mag(itt) + wt(k)*kappa_mid(k)
1155 dkappa_it1(k,itt) = kappa_pred(k) - kappa_out(k)
1156 dk_wt_it1(itt) = dk_wt_it1(itt) + wt(k)*dkappa_it1(k,itt)
1157 if (dkappa_it1(k,itt) > 0.0)
then
1158 dkpos_wt_it1(itt) = dkpos_wt_it1(itt) + wt(k)*dkappa_it1(k,itt)
1160 dkneg_wt_it1(itt) = dkneg_wt_it1(itt) + wt(k)*dkappa_it1(k,itt)
1162 wt_it1(k,itt) = wt(k)
1166 ri_k(k) = 1e3 ;
if (n2(k) < 1e3 * s2(k)) ri_k(k) = n2(k) / s2(k)
1167 dtke(k) = tke_pred(k) - tke(k)
1168 dtke_norm(k) = dtke(k) / (0.5*(tke(k) + tke_pred(k)))
1169 dkappa(k) = kappa_pred(k) - kappa_out(k)
1171 if (itt <= max_debug_itt)
then
1173 u_it1(k,itt) = u(k) ; v_it1(k,itt) = v(k)
1174 t_it1(k,itt) = t(k) ; s_it1(k,itt) = sal(k)
1177 kprev_it1(k,itt) = kappa_out(k)
1178 kappa_it1(k,itt) = kappa_mid(k) ; tke_it1(k,itt) = 0.5*(tke(k)+tke_pred(k))
1179 n2_it1(k,itt)=n2(k) ; sh2_it1(k,itt)=s2(k)
1180 ksrc_it1(k,itt) = kappa_src(k)
1181 k_q_it1(k,itt) = kappa_out(k) / (tke(k))
1183 if (abs(dkappa_it1(k,itt-1)) > 1e-20) &
1184 d_dkappa_it1(k,itt) = dkappa_it1(k,itt) / dkappa_it1(k,itt-1)
1186 dkappa_norm(k,itt) = dkappa(k) / max(0.5*(kappa_pred(k) + kappa_out(k)), us%m2_s_to_Z2_T*1e-100)
1191 if (dt_rem <= 0.0)
exit
1195 #ifdef ADD_DIAGNOSTICS
1196 if (
present(i_ld2_1d))
then
1197 do k=1,gv%ke+1 ; i_ld2_1d(k) = 0.0 ;
enddo
1198 do k=2,nzc ;
if (tke(k) > 0.0) &
1199 i_ld2_1d(k) = i_l2_bdry(k) + (n2(k) / cs%lambda**2 + f2) / tke(k)
1202 if (
present(dz_int_1d))
then
1203 do k=1,nzc+1 ; dz_int_1d(k) = dz_int(k) ;
enddo
1204 do k=nzc+2,gv%ke ; dz_int_1d(k) = 0.0 ;
enddo