Imposes the diapycnal mass fluxes and the accompanying diapycnal advection of momentum and tracers using the original MOM6 algorithms.
1913 type(ocean_grid_type),
intent(inout) :: G
1914 type(verticalGrid_type),
intent(in) :: GV
1915 type(unit_scale_type),
intent(in) :: US
1916 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: u
1917 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: v
1918 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: h
1919 type(thermo_var_ptrs),
intent(inout) :: tv
1921 real,
dimension(:,:),
pointer :: Hml
1922 type(forcing),
intent(inout) :: fluxes
1924 type(vertvisc_type),
intent(inout) :: visc
1925 type(accel_diag_ptrs),
intent(inout) :: ADp
1928 type(cont_diag_ptrs),
intent(inout) :: CDp
1929 real,
intent(in) :: dt
1930 type(time_type),
intent(in) :: Time_end
1931 type(diabatic_CS),
pointer :: CS
1932 type(Wave_parameters_CS),
optional,
pointer :: Waves
1934 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
1950 real,
dimension(SZI_(G),SZJ_(G)) :: &
1953 real,
dimension(SZI_(G),SZJ_(G),G%ke) :: h_diag
1954 real,
dimension(SZI_(G),SZJ_(G),G%ke) :: temp_diag
1955 real,
dimension(SZI_(G),SZJ_(G),G%ke) :: saln_diag
1956 real,
dimension(SZI_(G),SZJ_(G)) :: tendency_2d
1960 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
target :: &
1966 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
target :: &
1977 real,
pointer,
dimension(:,:,:) :: &
1983 integer :: kb(SZI_(G),SZJ_(G))
1986 real :: p_ref_cv(SZI_(G))
1990 logical :: in_boundary(SZI_(G))
2010 real :: htot(SZIB_(G))
2011 real :: b1(SZIB_(G)), d1(SZIB_(G))
2012 real :: c1(SZIB_(G),SZK_(G))
2020 logical :: showCallTree
2021 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, m, halo
2026 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
2027 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
2028 nkmb = gv%nk_rho_varies
2029 h_neglect = gv%H_subroundoff ; h_neglect2 = h_neglect*h_neglect
2030 kd_heat(:,:,:) = 0.0 ; kd_salt(:,:,:) = 0.0
2033 showcalltree = calltree_showquery()
2034 if (showcalltree)
call calltree_enter(
"layered_diabatic(), MOM_diabatic_driver.F90")
2037 eaml => eatr ; ebml => ebtr
2038 dt_in_t = dt * us%s_to_T
2041 call enable_averaging(dt, time_end, cs%diag)
2043 if ((cs%ML_mix_first > 0.0) .or. cs%use_geothermal)
then
2044 halo = cs%halo_TS_diff
2046 do k=1,nz ;
do j=js-halo,je+halo ;
do i=is-halo,ie+halo
2047 h_orig(i,j,k) = h(i,j,k) ; eaml(i,j,k) = 0.0 ; ebml(i,j,k) = 0.0
2048 enddo ;
enddo ;
enddo
2051 if (cs%use_geothermal)
then
2052 call cpu_clock_begin(id_clock_geothermal)
2053 call geothermal(h, tv, dt, eaml, ebml, g, gv, cs%geothermal_CSp, halo=cs%halo_TS_diff)
2054 call cpu_clock_end(id_clock_geothermal)
2055 if (showcalltree)
call calltree_waypoint(
"geothermal (diabatic)")
2056 if (cs%debugConservation)
call mom_state_stats(
'geothermal', u, v, h, tv%T, tv%S, g)
2061 call diag_update_remap_grids(cs%diag)
2066 if (
associated(cs%optics)) &
2067 call set_pen_shortwave(cs%optics, fluxes, g, gv, cs%diabatic_aux_CSp, cs%opacity_CSp, cs%tracer_flow_CSp)
2069 if (cs%bulkmixedlayer)
then
2070 if (cs%debug)
call mom_forcing_chksum(
"Before mixedlayer", fluxes, g, us, haloshift=0)
2072 if (cs%ML_mix_first > 0.0)
then
2080 call find_uv_at_h(u, v, h, u_h, v_h, g, gv)
2082 call cpu_clock_begin(id_clock_mixedlayer)
2083 if (cs%ML_mix_first < 1.0)
then
2085 call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt_in_t*cs%ML_mix_first, &
2086 eaml,ebml, g, gv, us, cs%bulkmixedlayer_CSp, cs%optics, &
2087 hml, cs%aggregate_FW_forcing, dt_in_t, last_call=.false.)
2088 if (cs%salt_reject_below_ML) &
2089 call insert_brine(h, tv, g, gv, fluxes, nkmb, cs%diabatic_aux_CSp, &
2090 dt*cs%ML_mix_first, cs%id_brine_lay)
2093 call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt_in_t, eaml, ebml, &
2094 g, gv, us, cs%bulkmixedlayer_CSp, cs%optics, &
2095 hml, cs%aggregate_FW_forcing, dt_in_t, last_call=.true.)
2104 if (
associated(tv%S) .and.
associated(tv%salt_deficit)) &
2105 call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)
2106 call cpu_clock_end(id_clock_mixedlayer)
2108 call mom_state_chksum(
"After mixedlayer ", u, v, h, g, gv, haloshift=0)
2109 call mom_forcing_chksum(
"After mixedlayer", fluxes, g, us, haloshift=0)
2111 if (showcalltree)
call calltree_waypoint(
"done with 1st bulkmixedlayer (diabatic)")
2112 if (cs%debugConservation)
call mom_state_stats(
'1st bulkmixedlayer', u, v, h, tv%T, tv%S, g)
2117 call mom_state_chksum(
"before find_uv_at_h", u, v, h, g, gv, haloshift=0)
2118 if (cs%use_kappa_shear .or. cs%use_CVMix_shear)
then
2119 if ((cs%ML_mix_first > 0.0) .or. cs%use_geothermal)
then
2120 call find_uv_at_h(u, v, h_orig, u_h, v_h, g, gv, eaml, ebml)
2122 call hchksum(eaml,
"after find_uv_at_h eaml",g%HI, scale=gv%H_to_m)
2123 call hchksum(ebml,
"after find_uv_at_h ebml",g%HI, scale=gv%H_to_m)
2126 call find_uv_at_h(u, v, h, u_h, v_h, g, gv)
2128 if (showcalltree)
call calltree_waypoint(
"done with find_uv_at_h (diabatic)")
2131 call cpu_clock_begin(id_clock_set_diffusivity)
2134 if ((cs%halo_TS_diff > 0) .and. (cs%ML_mix_first > 0.0))
then
2135 if (
associated(tv%T))
call pass_var(tv%T, g%Domain, halo=cs%halo_TS_diff, complete=.false.)
2136 if (
associated(tv%T))
call pass_var(tv%S, g%Domain, halo=cs%halo_TS_diff, complete=.false.)
2137 call pass_var(h, g%domain, halo=cs%halo_TS_diff, complete=.true.)
2139 call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt_in_t, g, gv, us, &
2140 cs%set_diff_CSp, kd_lay, kd_int)
2141 call cpu_clock_end(id_clock_set_diffusivity)
2142 if (showcalltree)
call calltree_waypoint(
"done with set_diffusivity (diabatic)")
2145 call mom_state_chksum(
"after set_diffusivity ", u, v, h, g, gv, haloshift=0)
2146 call mom_forcing_chksum(
"after set_diffusivity ", fluxes, g, us, haloshift=0)
2147 call mom_thermovar_chksum(
"after set_diffusivity ", tv, g)
2148 call hchksum(kd_lay,
"after set_diffusivity Kd_lay", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
2149 call hchksum(kd_int,
"after set_diffusivity Kd_Int", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
2154 call cpu_clock_begin(id_clock_kpp)
2160 call calculatebuoyancyflux2d(g, gv, us, fluxes, cs%optics, h, tv%T, tv%S, tv, &
2161 cs%KPP_buoy_flux, cs%KPP_temp_flux, cs%KPP_salt_flux)
2167 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
2168 kd_salt(i,j,k) = kd_int(i,j,k)
2169 kd_heat(i,j,k) = kd_int(i,j,k)
2170 enddo ;
enddo ;
enddo
2172 if (
associated(visc%Kd_extra_S))
then
2174 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
2175 kd_salt(i,j,k) = kd_salt(i,j,k) + visc%Kd_extra_S(i,j,k)
2176 enddo ;
enddo ;
enddo
2178 if (
associated(visc%Kd_extra_T))
then
2180 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
2181 kd_heat(i,j,k) = kd_heat(i,j,k) + visc%Kd_extra_T(i,j,k)
2182 enddo ;
enddo ;
enddo
2185 call kpp_compute_bld(cs%KPP_CSp, g, gv, us, h, tv%T, tv%S, u, v, tv%eqn_of_state, &
2186 fluxes%ustar, cs%KPP_buoy_flux, waves=waves)
2188 call kpp_calculate(cs%KPP_CSp, g, gv, us, h, fluxes%ustar, cs%KPP_buoy_flux, kd_heat, &
2189 kd_salt, visc%Kv_shear, cs%KPP_NLTheat, cs%KPP_NLTscalar, waves=waves)
2191 if (
associated(hml))
then
2192 call kpp_get_bld(cs%KPP_CSp, hml(:,:), g)
2193 call pass_var(hml, g%domain, halo=1)
2195 if (
associated(visc%MLD)) visc%MLD(:,:) = hml(:,:)
2198 if (.not. cs%KPPisPassive)
then
2200 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
2201 kd_int(i,j,k) = min( kd_salt(i,j,k), kd_heat(i,j,k) )
2202 enddo ;
enddo ;
enddo
2203 if (
associated(visc%Kd_extra_S))
then
2205 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
2206 visc%Kd_extra_S(i,j,k) = (kd_salt(i,j,k) - kd_int(i,j,k))
2207 enddo ;
enddo ;
enddo
2209 if (
associated(visc%Kd_extra_T))
then
2211 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
2212 visc%Kd_extra_T(i,j,k) = (kd_heat(i,j,k) - kd_int(i,j,k))
2213 enddo ;
enddo ;
enddo
2217 call cpu_clock_end(id_clock_kpp)
2218 if (showcalltree)
call calltree_waypoint(
"done with KPP_calculate (diabatic)")
2220 call mom_state_chksum(
"after KPP", u, v, h, g, gv, haloshift=0)
2221 call mom_forcing_chksum(
"after KPP", fluxes, g, us, haloshift=0)
2222 call mom_thermovar_chksum(
"after KPP", tv, g)
2223 call hchksum(kd_lay,
"after KPP Kd_lay", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
2224 call hchksum(kd_int,
"after KPP Kd_Int", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
2230 if (cs%use_CVMix_conv)
then
2231 call calculate_cvmix_conv(h, tv, g, gv, us, cs%CVMix_conv_csp, hml)
2233 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
2234 kd_int(i,j,k) = kd_int(i,j,k) + cs%CVMix_conv_csp%kd_conv(i,j,k)
2235 visc%Kv_slow(i,j,k) = visc%Kv_slow(i,j,k) + cs%CVMix_conv_csp%kv_conv(i,j,k)
2236 enddo ;
enddo ;
enddo
2240 call cpu_clock_begin(id_clock_kpp)
2242 call hchksum(cs%KPP_temp_flux,
"before KPP_applyNLT netHeat", g%HI, haloshift=0, scale=gv%H_to_m)
2243 call hchksum(cs%KPP_salt_flux,
"before KPP_applyNLT netSalt", g%HI, haloshift=0, scale=gv%H_to_m)
2244 call hchksum(cs%KPP_NLTheat,
"before KPP_applyNLT NLTheat", g%HI, haloshift=0)
2245 call hchksum(cs%KPP_NLTscalar,
"before KPP_applyNLT NLTscalar", g%HI, haloshift=0)
2249 call kpp_nonlocaltransport_temp(cs%KPP_CSp, g, gv, h, cs%KPP_NLTheat, cs%KPP_temp_flux, dt, tv%T, tv%C_p)
2250 call kpp_nonlocaltransport_saln(cs%KPP_CSp, g, gv, h, cs%KPP_NLTscalar, cs%KPP_salt_flux, dt, tv%S)
2251 call cpu_clock_end(id_clock_kpp)
2252 if (showcalltree)
call calltree_waypoint(
"done with KPP_applyNonLocalTransport (diabatic)")
2253 if (cs%debugConservation)
call mom_state_stats(
'KPP_applyNonLocalTransport', u, v, h, tv%T, tv%S, g)
2256 call mom_state_chksum(
"after KPP_applyNLT ", u, v, h, g, gv, haloshift=0)
2257 call mom_forcing_chksum(
"after KPP_applyNLT ", fluxes, g, us, haloshift=0)
2258 call mom_thermovar_chksum(
"after KPP_applyNLT ", tv, g)
2264 if (
associated(visc%Kd_extra_T) .and.
associated(visc%Kd_extra_S) .and.
associated(tv%T))
then
2266 call cpu_clock_begin(id_clock_differential_diff)
2267 call differential_diffuse_t_s(h, tv, visc, dt_in_t, g, gv)
2268 call cpu_clock_end(id_clock_differential_diff)
2269 if (showcalltree)
call calltree_waypoint(
"done with differential_diffuse_T_S (diabatic)")
2270 if (cs%debugConservation)
call mom_state_stats(
'differential_diffuse_T_S', u, v, h, tv%T, tv%S, g)
2274 if (.not. cs%useKPP)
then
2276 do k=2,nz ;
do j=js,je ;
do i=is,ie
2277 kd_heat(i,j,k) = kd_heat(i,j,k) + visc%Kd_extra_T(i,j,k)
2278 kd_salt(i,j,k) = kd_salt(i,j,k) + visc%Kd_extra_S(i,j,k)
2279 enddo ;
enddo ;
enddo
2289 call cpu_clock_begin(id_clock_entrain)
2292 call entrainment_diffusive(h, tv, fluxes, dt_in_t, g, gv, us, cs%entrain_diffusive_CSp, &
2293 ea, eb, kb, kd_lay=kd_lay, kd_int=kd_int)
2294 call cpu_clock_end(id_clock_entrain)
2295 if (showcalltree)
call calltree_waypoint(
"done with Entrainment_diffusive (diabatic)")
2298 call mom_forcing_chksum(
"after calc_entrain ", fluxes, g, us, haloshift=0)
2299 call mom_thermovar_chksum(
"after calc_entrain ", tv, g)
2300 call mom_state_chksum(
"after calc_entrain ", u, v, h, g, gv, haloshift=0)
2301 call hchksum(ea,
"after calc_entrain ea", g%HI, haloshift=0, scale=gv%H_to_m)
2302 call hchksum(eb,
"after calc_entrain eb", g%HI, haloshift=0, scale=gv%H_to_m)
2306 if (cs%boundary_forcing_tendency_diag)
then
2307 do k=1,nz ;
do j=js,je ;
do i=is,ie
2308 h_diag(i,j,k) = h(i,j,k)
2309 temp_diag(i,j,k) = tv%T(i,j,k)
2310 saln_diag(i,j,k) = tv%S(i,j,k)
2311 enddo ;
enddo ;
enddo
2325 hold(i,j,1) = h(i,j,1)
2326 h(i,j,1) = h(i,j,1) + (eb(i,j,1) - ea(i,j,2))
2327 hold(i,j,nz) = h(i,j,nz)
2328 h(i,j,nz) = h(i,j,nz) + (ea(i,j,nz) - eb(i,j,nz-1))
2329 if (h(i,j,1) <= 0.0)
then
2330 h(i,j,1) = gv%Angstrom_H
2332 if (h(i,j,nz) <= 0.0)
then
2333 h(i,j,nz) = gv%Angstrom_H
2336 do k=2,nz-1 ;
do i=is,ie
2337 hold(i,j,k) = h(i,j,k)
2338 h(i,j,k) = h(i,j,k) + ((ea(i,j,k) - eb(i,j,k-1)) + &
2339 (eb(i,j,k) - ea(i,j,k+1)))
2340 if (h(i,j,k) <= 0.0)
then
2341 h(i,j,k) = gv%Angstrom_H
2346 call diag_update_remap_grids(cs%diag)
2349 call mom_state_chksum(
"after negative check ", u, v, h, g, gv, haloshift=0)
2350 call mom_forcing_chksum(
"after negative check ", fluxes, g, us, haloshift=0)
2351 call mom_thermovar_chksum(
"after negative check ", tv, g)
2353 if (showcalltree)
call calltree_waypoint(
"done with h=ea-eb (diabatic)")
2354 if (cs%debugConservation)
call mom_state_stats(
'h=ea-eb', u, v, h, tv%T, tv%S, g)
2360 if (cs%bulkmixedlayer)
then
2362 if (
associated(tv%T))
then
2363 call cpu_clock_begin(id_clock_tridiag)
2368 if (cs%massless_match_targets)
then
2372 h_tr = hold(i,j,1) + h_neglect
2373 b1(i) = 1.0 / (h_tr + eb(i,j,1))
2374 d1(i) = h_tr * b1(i)
2375 tv%T(i,j,1) = b1(i) * (h_tr*tv%T(i,j,1))
2376 tv%S(i,j,1) = b1(i) * (h_tr*tv%S(i,j,1))
2378 do k=2,nkmb ;
do i=is,ie
2379 c1(i,k) = eb(i,j,k-1) * b1(i)
2380 h_tr = hold(i,j,k) + h_neglect
2381 b_denom_1 = h_tr + d1(i)*ea(i,j,k)
2382 b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
2383 if (k<nkmb) d1(i) = b_denom_1 * b1(i)
2384 tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,k-1))
2385 tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,k-1))
2388 do k=nkmb+1,nz ;
do i=is,ie
2389 if (k == kb(i,j))
then
2390 c1(i,k) = eb(i,j,k-1) * b1(i)
2391 d1(i) = (((eb(i,j,nkmb)-eb(i,j,k-1)) + hold(i,j,nkmb) + h_neglect) + &
2392 d1(i)*ea(i,j,nkmb)) * b1(i)
2393 h_tr = hold(i,j,k) + h_neglect
2394 b_denom_1 = h_tr + d1(i)*ea(i,j,k)
2395 b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
2396 d1(i) = b_denom_1 * b1(i)
2397 tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,nkmb))
2398 tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,nkmb))
2399 elseif (k > kb(i,j))
then
2400 c1(i,k) = eb(i,j,k-1) * b1(i)
2401 h_tr = hold(i,j,k) + h_neglect
2402 b_denom_1 = h_tr + d1(i)*ea(i,j,k)
2403 b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
2404 d1(i) = b_denom_1 * b1(i)
2405 tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,k-1))
2406 tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,k-1))
2407 elseif (eb(i,j,k) < eb(i,j,k-1))
then
2412 tv%T(i,j,nkmb) = tv%T(i,j,nkmb) + b1(i) * (eb(i,j,k-1) - eb(i,j,k)) * tv%T(i,j,k)
2413 tv%S(i,j,nkmb) = tv%S(i,j,nkmb) + b1(i) * (eb(i,j,k-1) - eb(i,j,k)) * tv%S(i,j,k)
2417 do k=nz-1,nkmb,-1 ;
do i=is,ie
2418 if (k >= kb(i,j))
then
2419 tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1)
2420 tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1)
2423 do i=is,ie ;
if (kb(i,j) <= nz)
then
2424 tv%T(i,j,nkmb) = tv%T(i,j,nkmb) + c1(i,kb(i,j))*tv%T(i,j,kb(i,j))
2425 tv%S(i,j,nkmb) = tv%S(i,j,nkmb) + c1(i,kb(i,j))*tv%S(i,j,kb(i,j))
2427 do k=nkmb-1,1,-1 ;
do i=is,ie
2428 tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1)
2429 tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1)
2436 if (cs%tracer_tridiag)
then
2437 call tracer_vertdiff(hold, ea, eb, dt, tv%T, g, gv)
2438 call tracer_vertdiff(hold, ea, eb, dt, tv%S, g, gv)
2440 call tridiagts(g, gv, is, ie, js, je, hold, ea, eb, tv%T, tv%S)
2443 call cpu_clock_end(id_clock_tridiag)
2446 if (cs%debugConservation)
call mom_state_stats(
'BML tridiag', u, v, h, tv%T, tv%S, g)
2448 if ((cs%ML_mix_first > 0.0) .or. cs%use_geothermal)
then
2452 do k=1,nz ;
do j=js,je ;
do i=is,ie
2453 hold(i,j,k) = h_orig(i,j,k)
2454 ea(i,j,k) = ea(i,j,k) + eaml(i,j,k)
2455 eb(i,j,k) = eb(i,j,k) + ebml(i,j,k)
2456 enddo ;
enddo ;
enddo
2458 call hchksum(ea,
"after ea = ea + eaml",g%HI,haloshift=0, scale=gv%H_to_m)
2459 call hchksum(eb,
"after eb = eb + ebml",g%HI,haloshift=0, scale=gv%H_to_m)
2463 if (cs%ML_mix_first < 1.0)
then
2472 call find_uv_at_h(u, v, hold, u_h, v_h, g, gv, ea, eb)
2473 if (cs%debug)
call mom_state_chksum(
"find_uv_at_h1 ", u, v, h, g, gv, haloshift=0)
2475 dt_mix = min(dt_in_t, dt_in_t*(1.0 - cs%ML_mix_first))
2476 call cpu_clock_begin(id_clock_mixedlayer)
2478 call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt_mix, ea, eb, &
2479 g, gv, us, cs%bulkmixedlayer_CSp, cs%optics, &
2480 hml, cs%aggregate_FW_forcing, dt_in_t, last_call=.true.)
2482 if (cs%salt_reject_below_ML) &
2483 call insert_brine(h, tv, g, gv, fluxes, nkmb, cs%diabatic_aux_CSp, us%T_to_s*dt_mix, &
2492 if (
associated(tv%S) .and.
associated(tv%salt_deficit)) &
2493 call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)
2495 call cpu_clock_end(id_clock_mixedlayer)
2496 if (showcalltree)
call calltree_waypoint(
"done with 2nd bulkmixedlayer (diabatic)")
2497 if (cs%debugConservation)
call mom_state_stats(
'2nd bulkmixedlayer', u, v, h, tv%T, tv%S, g)
2503 if (
associated(tv%T))
then
2506 call hchksum(ea,
"before triDiagTS ea ",g%HI,haloshift=0, scale=gv%H_to_m)
2507 call hchksum(eb,
"before triDiagTS eb ",g%HI,haloshift=0, scale=gv%H_to_m)
2509 call cpu_clock_begin(id_clock_tridiag)
2517 if (
associated(tv%S) .and.
associated(tv%salt_deficit)) &
2518 call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)
2520 if (cs%diabatic_diff_tendency_diag)
then
2521 do k=1,nz ;
do j=js,je ;
do i=is,ie
2522 temp_diag(i,j,k) = tv%T(i,j,k)
2523 saln_diag(i,j,k) = tv%S(i,j,k)
2524 enddo ;
enddo ;
enddo
2528 if (cs%tracer_tridiag)
then
2529 call tracer_vertdiff(hold, ea, eb, dt, tv%T, g, gv)
2530 call tracer_vertdiff(hold, ea, eb, dt, tv%S, g, gv)
2532 call tridiagts(g, gv, is, ie, js, je, hold, ea, eb, tv%T, tv%S)
2538 if (cs%diabatic_diff_tendency_diag)
then
2539 call diagnose_diabatic_diff_tendency(tv, hold, temp_diag, saln_diag, dt, g, gv, cs)
2540 if (cs%id_diabatic_diff_h > 0)
call post_data(cs%id_diabatic_diff_h, hold, cs%diag, alt_h = hold)
2543 call cpu_clock_end(id_clock_tridiag)
2544 if (showcalltree)
call calltree_waypoint(
"done with triDiagTS (diabatic)")
2547 if (cs%debugConservation)
call mom_state_stats(
'triDiagTS', u, v, h, tv%T, tv%S, g)
2552 call mom_state_chksum(
"after mixed layer ", u, v, h, g, gv, haloshift=0)
2553 call mom_thermovar_chksum(
"after mixed layer ", tv, g)
2554 call hchksum(ea,
"after mixed layer ea", g%HI, scale=gv%H_to_m)
2555 call hchksum(eb,
"after mixed layer eb", g%HI, scale=gv%H_to_m)
2558 call cpu_clock_begin(id_clock_remap)
2559 call regularize_layers(h, tv, dt, ea, eb, g, gv, cs%regularize_layers_CSp)
2560 call cpu_clock_end(id_clock_remap)
2561 if (showcalltree)
call calltree_waypoint(
"done with regularize_layers (diabatic)")
2562 if (cs%debugConservation)
call mom_state_stats(
'regularize_layers', u, v, h, tv%T, tv%S, g)
2566 if (cs%id_dudt_dia > 0 .or. cs%id_dvdt_dia > 0) &
2568 call pass_var(h, g%domain, to_west+to_south+omit_corners, halo=1)
2569 call diag_update_remap_grids(cs%diag)
2573 if ((cs%id_Tdif > 0) .or. (cs%id_Tadv > 0))
then
2574 do j=js,je ;
do i=is,ie
2575 tdif_flx(i,j,1) = 0.0 ; tdif_flx(i,j,nz+1) = 0.0
2576 tadv_flx(i,j,1) = 0.0 ; tadv_flx(i,j,nz+1) = 0.0
2579 do k=2,nz ;
do j=js,je ;
do i=is,ie
2580 tdif_flx(i,j,k) = (idt * 0.5*(ea(i,j,k) + eb(i,j,k-1))) * &
2581 (tv%T(i,j,k-1) - tv%T(i,j,k))
2582 tadv_flx(i,j,k) = (idt * (ea(i,j,k) - eb(i,j,k-1))) * &
2583 0.5*(tv%T(i,j,k-1) + tv%T(i,j,k))
2584 enddo ;
enddo ;
enddo
2586 if ((cs%id_Sdif > 0) .or. (cs%id_Sadv > 0))
then
2587 do j=js,je ;
do i=is,ie
2588 sdif_flx(i,j,1) = 0.0 ; sdif_flx(i,j,nz+1) = 0.0
2589 sadv_flx(i,j,1) = 0.0 ; sadv_flx(i,j,nz+1) = 0.0
2592 do k=2,nz ;
do j=js,je ;
do i=is,ie
2593 sdif_flx(i,j,k) = (idt * 0.5*(ea(i,j,k) + eb(i,j,k-1))) * &
2594 (tv%S(i,j,k-1) - tv%S(i,j,k))
2595 sadv_flx(i,j,k) = (idt * (ea(i,j,k) - eb(i,j,k-1))) * &
2596 0.5*(tv%S(i,j,k-1) + tv%S(i,j,k))
2597 enddo ;
enddo ;
enddo
2601 call cpu_clock_begin(id_clock_tracers)
2602 if (cs%mix_boundary_tracers)
then
2603 tr_ea_bbl = gv%Z_to_H * sqrt(dt_in_t*cs%Kd_BBL_tr)
2607 ebtr(i,j,nz) = eb(i,j,nz)
2609 in_boundary(i) = (g%mask2dT(i,j) > 0.0)
2611 do k=nz,2,-1 ;
do i=is,ie
2612 if (in_boundary(i))
then
2613 htot(i) = htot(i) + h(i,j,k)
2622 add_ent = ((dt_in_t * cs%Kd_min_tr) * gv%Z_to_H**2) * &
2623 ((h(i,j,k-1)+h(i,j,k)+h_neglect) / &
2624 (h(i,j,k-1)*h(i,j,k)+h_neglect2)) - &
2625 0.5*(ea(i,j,k) + eb(i,j,k-1))
2626 if (htot(i) < tr_ea_bbl)
then
2627 add_ent = max(0.0, add_ent, &
2628 (tr_ea_bbl - htot(i)) - min(ea(i,j,k),eb(i,j,k-1)))
2629 elseif (add_ent < 0.0)
then
2630 add_ent = 0.0 ; in_boundary(i) = .false.
2633 ebtr(i,j,k-1) = eb(i,j,k-1) + add_ent
2634 eatr(i,j,k) = ea(i,j,k) + add_ent
2636 ebtr(i,j,k-1) = eb(i,j,k-1) ; eatr(i,j,k) = ea(i,j,k)
2638 if (
associated(visc%Kd_extra_S))
then ;
if (visc%Kd_extra_S(i,j,k) > 0.0)
then
2639 add_ent = ((dt_in_t * visc%Kd_extra_S(i,j,k)) * gv%Z_to_H**2) / &
2640 (0.25 * ((h(i,j,k-1) + h(i,j,k)) + (hold(i,j,k-1) + hold(i,j,k))) + &
2642 ebtr(i,j,k-1) = ebtr(i,j,k-1) + add_ent
2643 eatr(i,j,k) = eatr(i,j,k) + add_ent
2646 do i=is,ie ; eatr(i,j,1) = ea(i,j,1) ;
enddo
2650 call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, hml, dt, g, gv, tv, &
2651 cs%optics, cs%tracer_flow_CSp, cs%debug)
2653 elseif (
associated(visc%Kd_extra_S))
then
2655 do j=js,je ;
do i=is,ie
2656 ebtr(i,j,nz) = eb(i,j,nz) ; eatr(i,j,1) = ea(i,j,1)
2659 do k=nz,2,-1 ;
do j=js,je ;
do i=is,ie
2660 if (visc%Kd_extra_S(i,j,k) > 0.0)
then
2661 add_ent = ((dt_in_t * visc%Kd_extra_S(i,j,k)) * gv%Z_to_H**2) / &
2662 (0.25 * ((h(i,j,k-1) + h(i,j,k)) + (hold(i,j,k-1) + hold(i,j,k))) + &
2667 ebtr(i,j,k-1) = eb(i,j,k-1) + add_ent
2668 eatr(i,j,k) = ea(i,j,k) + add_ent
2669 enddo ;
enddo ;
enddo
2671 call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, hml, dt, g, gv, tv, &
2672 cs%optics, cs%tracer_flow_CSp, cs%debug)
2675 call call_tracer_column_fns(hold, h, ea, eb, fluxes, hml, dt, g, gv, tv, &
2676 cs%optics, cs%tracer_flow_CSp, cs%debug)
2680 call cpu_clock_end(id_clock_tracers)
2683 if (cs%use_sponge)
then
2684 call cpu_clock_begin(id_clock_sponge)
2686 if (cs%bulkmixedlayer .and.
associated(tv%eqn_of_state))
then
2687 do i=is,ie ; p_ref_cv(i) = tv%P_Ref ;
enddo
2690 call calculate_density(tv%T(:,j,1), tv%S(:,j,1), p_ref_cv, rcv_ml(:,j), &
2691 is, ie-is+1, tv%eqn_of_state)
2693 call apply_sponge(h, dt, g, gv, ea, eb, cs%sponge_CSp, rcv_ml)
2695 call apply_sponge(h, dt, g, gv, ea, eb, cs%sponge_CSp)
2697 call cpu_clock_end(id_clock_sponge)
2699 call mom_state_chksum(
"apply_sponge ", u, v, h, g, gv, haloshift=0)
2700 call mom_thermovar_chksum(
"apply_sponge ", tv, g)
2705 if (
associated(cdp%diapyc_vel))
then
2708 do k=2,nz ;
do i=is,ie
2709 cdp%diapyc_vel(i,j,k) = idt * (gv%H_to_m * (ea(i,j,k) - eb(i,j,k-1)))
2712 cdp%diapyc_vel(i,j,1) = 0.0
2713 cdp%diapyc_vel(i,j,nz+1) = 0.0
2721 if (cs%bulkmixedlayer)
then
2723 call hchksum(ea,
"before net flux rearrangement ea",g%HI, scale=gv%H_to_m)
2724 call hchksum(eb,
"before net flux rearrangement eb",g%HI, scale=gv%H_to_m)
2728 do k=2,gv%nkml ;
do i=is,ie
2729 net_ent = ea(i,j,k) - eb(i,j,k-1)
2730 ea(i,j,k) = max(net_ent, 0.0)
2731 eb(i,j,k-1) = max(-net_ent, 0.0)
2735 call hchksum(ea,
"after net flux rearrangement ea",g%HI, scale=gv%H_to_m)
2736 call hchksum(eb,
"after net flux rearrangement eb",g%HI, scale=gv%H_to_m)
2744 hold(i,js-1,k) = gv%Angstrom_H ; ea(i,js-1,k) = 0.0 ; eb(i,js-1,k) = 0.0
2745 hold(i,je+1,k) = gv%Angstrom_H ; ea(i,je+1,k) = 0.0 ; eb(i,je+1,k) = 0.0
2748 hold(is-1,j,k) = gv%Angstrom_H ; ea(is-1,j,k) = 0.0 ; eb(is-1,j,k) = 0.0
2749 hold(ie+1,j,k) = gv%Angstrom_H ; ea(ie+1,j,k) = 0.0 ; eb(ie+1,j,k) = 0.0
2753 call cpu_clock_begin(id_clock_pass)
2754 if (g%symmetric)
then ; dir_flag = to_all+omit_corners
2755 else ; dir_flag = to_west+to_south+omit_corners ;
endif
2756 call create_group_pass(cs%pass_hold_eb_ea, hold, g%Domain, dir_flag, halo=1)
2757 call create_group_pass(cs%pass_hold_eb_ea, eb, g%Domain, dir_flag, halo=1)
2758 call create_group_pass(cs%pass_hold_eb_ea, ea, g%Domain, dir_flag, halo=1)
2759 call do_group_pass(cs%pass_hold_eb_ea, g%Domain)
2760 call cpu_clock_end(id_clock_pass)
2766 call mom_state_chksum(
"before u/v tridiag ", u, v, h, g, gv, haloshift=0)
2767 call hchksum(ea,
"before u/v tridiag ea",g%HI, scale=gv%H_to_m)
2768 call hchksum(eb,
"before u/v tridiag eb",g%HI, scale=gv%H_to_m)
2769 call hchksum(hold,
"before u/v tridiag hold",g%HI, scale=gv%H_to_m)
2771 call cpu_clock_begin(id_clock_tridiag)
2775 if (
associated(adp%du_dt_dia)) adp%du_dt_dia(i,j,1) = u(i,j,1)
2776 hval = (hold(i,j,1) + hold(i+1,j,1)) + (ea(i,j,1) + ea(i+1,j,1)) + h_neglect
2777 b1(i) = 1.0 / (hval + (eb(i,j,1) + eb(i+1,j,1)))
2778 d1(i) = hval * b1(i)
2779 u(i,j,1) = b1(i) * (hval * u(i,j,1))
2781 do k=2,nz ;
do i=isq,ieq
2782 if (
associated(adp%du_dt_dia)) adp%du_dt_dia(i,j,k) = u(i,j,k)
2783 c1(i,k) = (eb(i,j,k-1)+eb(i+1,j,k-1)) * b1(i)
2784 eaval = ea(i,j,k) + ea(i+1,j,k)
2785 hval = hold(i,j,k) + hold(i+1,j,k) + h_neglect
2786 b1(i) = 1.0 / ((eb(i,j,k) + eb(i+1,j,k)) + (hval + d1(i)*eaval))
2787 d1(i) = (hval + d1(i)*eaval) * b1(i)
2788 u(i,j,k) = (hval*u(i,j,k) + eaval*u(i,j,k-1))*b1(i)
2790 do k=nz-1,1,-1 ;
do i=isq,ieq
2791 u(i,j,k) = u(i,j,k) + c1(i,k+1)*u(i,j,k+1)
2792 if (
associated(adp%du_dt_dia)) &
2793 adp%du_dt_dia(i,j,k) = (u(i,j,k) - adp%du_dt_dia(i,j,k)) * idt
2795 if (
associated(adp%du_dt_dia))
then
2797 adp%du_dt_dia(i,j,nz) = (u(i,j,nz)-adp%du_dt_dia(i,j,nz)) * idt
2802 call mom_state_chksum(
"aft 1st loop tridiag ", u, v, h, g, gv, haloshift=0)
2807 if (
associated(adp%dv_dt_dia)) adp%dv_dt_dia(i,j,1) = v(i,j,1)
2808 hval = (hold(i,j,1) + hold(i,j+1,1)) + (ea(i,j,1) + ea(i,j+1,1)) + h_neglect
2809 b1(i) = 1.0 / (hval + (eb(i,j,1) + eb(i,j+1,1)))
2810 d1(i) = hval * b1(i)
2811 v(i,j,1) = b1(i) * (hval * v(i,j,1))
2813 do k=2,nz ;
do i=is,ie
2814 if (
associated(adp%dv_dt_dia)) adp%dv_dt_dia(i,j,k) = v(i,j,k)
2815 c1(i,k) = (eb(i,j,k-1)+eb(i,j+1,k-1)) * b1(i)
2816 eaval = ea(i,j,k) + ea(i,j+1,k)
2817 hval = hold(i,j,k) + hold(i,j+1,k) + h_neglect
2818 b1(i) = 1.0 / ((eb(i,j,k) + eb(i,j+1,k)) + (hval + d1(i)*eaval))
2819 d1(i) = (hval + d1(i)*eaval) * b1(i)
2820 v(i,j,k) = (hval*v(i,j,k) + eaval*v(i,j,k-1))*b1(i)
2822 do k=nz-1,1,-1 ;
do i=is,ie
2823 v(i,j,k) = v(i,j,k) + c1(i,k+1)*v(i,j,k+1)
2824 if (
associated(adp%dv_dt_dia)) &
2825 adp%dv_dt_dia(i,j,k) = (v(i,j,k) - adp%dv_dt_dia(i,j,k)) * idt
2827 if (
associated(adp%dv_dt_dia))
then
2829 adp%dv_dt_dia(i,j,nz) = (v(i,j,nz)-adp%dv_dt_dia(i,j,nz)) * idt
2833 call cpu_clock_end(id_clock_tridiag)
2835 call mom_state_chksum(
"after u/v tridiag ", u, v, h, g, gv, haloshift=0)
2838 call disable_averaging(cs%diag)
2840 call enable_averaging(dt, time_end, cs%diag)
2842 if (cs%id_Kd_interface > 0)
call post_data(cs%id_Kd_interface, kd_int, cs%diag)
2843 if (cs%id_Kd_heat > 0)
call post_data(cs%id_Kd_heat, kd_heat, cs%diag)
2844 if (cs%id_Kd_salt > 0)
call post_data(cs%id_Kd_salt, kd_salt, cs%diag)
2845 if (cs%id_Kd_ePBL > 0)
call post_data(cs%id_Kd_ePBL, kd_epbl, cs%diag)
2847 if (cs%id_ea > 0)
call post_data(cs%id_ea, ea, cs%diag)
2848 if (cs%id_eb > 0)
call post_data(cs%id_eb, eb, cs%diag)
2850 if (cs%id_dudt_dia > 0)
call post_data(cs%id_dudt_dia, adp%du_dt_dia, cs%diag)
2851 if (cs%id_dvdt_dia > 0)
call post_data(cs%id_dvdt_dia, adp%dv_dt_dia, cs%diag)
2852 if (cs%id_wd > 0)
call post_data(cs%id_wd, cdp%diapyc_vel, cs%diag)
2854 if (cs%id_Tdif > 0)
call post_data(cs%id_Tdif, tdif_flx, cs%diag)
2855 if (cs%id_Tadv > 0)
call post_data(cs%id_Tadv, tadv_flx, cs%diag)
2856 if (cs%id_Sdif > 0)
call post_data(cs%id_Sdif, sdif_flx, cs%diag)
2857 if (cs%id_Sadv > 0)
call post_data(cs%id_Sadv, sadv_flx, cs%diag)
2859 call disable_averaging(cs%diag)
2861 if (showcalltree)
call calltree_leave(
"layered_diabatic()")