Build a SLight coordinate column.
183 type(slight_CS),
intent(in) :: CS
184 type(EOS_type),
pointer :: eqn_of_state
185 real,
intent(in) :: H_to_Pa
186 real,
intent(in) :: H_subroundoff
187 integer,
intent(in) :: nz
188 real,
intent(in) :: depth
189 real,
dimension(nz),
intent(in) :: T_col
190 real,
dimension(nz),
intent(in) :: S_col
191 real,
dimension(nz),
intent(in) :: h_col
192 real,
dimension(nz),
intent(in) :: p_col
193 real,
dimension(nz+1),
intent(in) :: z_col
194 real,
dimension(nz+1),
intent(inout) :: z_col_new
195 real,
optional,
intent(in) :: h_neglect
197 real,
optional,
intent(in) :: h_neglect_edge
200 real,
dimension(nz) :: rho_col
201 real,
dimension(nz) :: T_f, S_f
202 logical,
dimension(nz+1) :: reliable
203 real,
dimension(nz+1) :: T_int, S_int
204 real,
dimension(nz+1) :: rho_tmp, drho_dp, p_IS, p_R
205 real,
dimension(nz+1) :: drhoIS_dT, drhoIS_dS
206 real,
dimension(nz+1) :: drhoR_dT, drhoR_dS
207 real,
dimension(nz+1) :: strat_rat
209 real :: drIS, drR, Fn_now, I_HStol, Fn_zero_val
225 logical :: maximum_depths_set
226 logical :: maximum_h_set
227 real :: k2_used, k2here, dz_sum, z_max
229 real :: h_tr, b_denom_1, b1, d1
230 real,
dimension(nz) :: c1
231 integer :: kur1, kur2
233 integer :: i, j, k, nkml
235 maximum_depths_set =
allocated(cs%max_interface_depths)
236 maximum_h_set =
allocated(cs%max_layer_thickness)
238 if (z_col(nz+1) - z_col(1) < nz*cs%min_thickness)
then
240 dz = (z_col(nz+1) - z_col(1)) / real(nz)
241 do k=2,nz ; z_col_new(k) = z_col(1) + dz*real(k-1) ;
enddo
243 call calculate_density(t_col, s_col, p_col, rho_col, 1, nz, &
248 call rho_interfaces_col(rho_col, h_col, z_col, cs%target_density, nz, &
249 z_col_new, cs, reliable, debug=.true., &
250 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
253 if (cs%min_thickness > 0.0)
then
255 do k=2,nz ;
if (z_col_new(k) < z_col_new(k-1) + cs%min_thickness)
then
256 z_col_new(k) = z_col_new(k-1) + cs%min_thickness
259 do k=nz,2,-1 ;
if (z_col_new(k) > z_col_new(k+1) - cs%min_thickness)
then
260 z_col_new(k) = z_col_new(k+1) - cs%min_thickness
271 do k=kur_ss,nz ;
if (.not.reliable(k))
then
277 do k=kur1+1,nz+1 ;
if (reliable(k))
then
278 kur2 = k-1 ; kur_ss = k ;
exit
280 if (kur2 < kur1)
call mom_error(fatal,
"Bad unreliable range.")
282 dz_ur = z_col_new(kur2+1) - z_col_new(kur1-1)
286 wgt = 1.0 ; cowgt = 0.0
288 z_col_new(k) = cowgt*z_col_new(k) + &
289 wgt * (z_col_new(kur1-1) + dz_ur*(k - (kur1-1)) / ((kur2 - kur1) + 2))
295 z_wt = 0.0 ; rho_x_z = 0.0
296 h_ml_av = cs%Rho_ml_avg_depth
298 if (z_wt + h_col(k) >= h_ml_av)
then
299 rho_x_z = rho_x_z + rho_col(k) * (h_ml_av - z_wt)
303 rho_x_z = rho_x_z + rho_col(k) * h_col(k)
304 z_wt = z_wt + h_col(k)
307 if (z_wt > 0.0) rho_ml_av = rho_x_z / z_wt
309 nkml = cs%nz_fixed_surface
311 if (rho_ml_av <= cs%target_density(nkml))
then
312 k_interior = cs%nlay_ml_offset + real(nkml)
313 elseif (rho_ml_av > cs%target_density(nz+1))
then
314 k_interior = real(nz+1)
316 if ((rho_ml_av >= cs%target_density(k)) .and. &
317 (rho_ml_av < cs%target_density(k+1)))
then
318 k_interior = (cs%nlay_ml_offset + k) + &
319 (rho_ml_av - cs%target_density(k)) / &
320 (cs%target_density(k+1) - cs%target_density(k))
324 if (k_interior > real(nz+1)) k_interior = real(nz+1)
327 k = int(ceiling(k_interior))
328 z_interior = (k-k_interior)*z_col_new(k-1) + (1.0+(k_interior-k))*z_col_new(k)
330 if (cs%fix_haloclines)
then
336 if (cs%halocline_filter_length > 0.0)
then
337 lfilt = cs%halocline_filter_length
340 h_tr = h_col(1) + h_subroundoff
341 b1 = 1.0 / (h_tr + lfilt) ; d1 = h_tr * b1
342 t_f(1) = (b1*h_tr)*t_col(1) ; s_f(1) = (b1*h_tr)*s_col(1)
345 h_tr = h_col(k) + h_subroundoff ; b_denom_1 = h_tr + d1*lfilt
346 b1 = 1.0 / (b_denom_1 + lfilt) ; d1 = b_denom_1 * b1
347 t_f(k) = b1 * (h_tr*t_col(k) + lfilt*t_f(k-1))
348 s_f(k) = b1 * (h_tr*s_col(k) + lfilt*s_f(k-1))
351 t_f(k) = t_f(k) + c1(k+1)*t_f(k+1) ; s_f(k) = s_f(k) + c1(k+1)*s_f(k+1)
354 do k=1,nz ; t_f(k) = t_col(k) ; s_f(k) = s_col(k) ;
enddo
357 t_int(1) = t_f(1) ; s_int(1) = s_f(1)
359 t_int(k) = 0.5*(t_f(k-1) + t_f(k)) ; s_int(k) = 0.5*(s_f(k-1) + s_f(k))
360 p_is(k) = z_col(k) * h_to_pa
361 p_r(k) = cs%ref_pressure + cs%compressibility_fraction * ( p_is(k) - cs%ref_pressure )
363 t_int(nz+1) = t_f(nz) ; s_int(nz+1) = s_f(nz)
364 p_is(nz+1) = z_col(nz+1) * h_to_pa
365 call calculate_density_derivs(t_int, s_int, p_is, drhois_dt, drhois_ds, 2, nz-1, &
367 call calculate_density_derivs(t_int, s_int, p_r, drhor_dt, drhor_ds, 2, nz-1, &
369 if (cs%compressibility_fraction > 0.0)
then
370 call calculate_compress(t_int, s_int, p_r, rho_tmp, drho_dp, 2, nz-1, &
373 do k=2,nz ; drho_dp(k) = 0.0 ;
enddo
376 h_to_cpa = cs%compressibility_fraction*h_to_pa
379 dris = drhois_dt(k) * (t_f(k) - t_f(k-1)) + &
380 drhois_ds(k) * (s_f(k) - s_f(k-1))
381 drr = (drhor_dt(k) * (t_f(k) - t_f(k-1)) + &
382 drhor_ds(k) * (s_f(k) - s_f(k-1))) + &
383 drho_dp(k) * (h_to_cpa*0.5*(h_col(k) + h_col(k-1)))
385 if (dris <= 0.0)
then
388 strat_rat(k) = 2.0*max(drr,0.0) / (dris + abs(drr))
391 strat_rat(nz+1) = 1.0
393 z_int_unst = 0.0 ; fn_now = 0.0
394 fn_zero_val = min(2.0*cs%halocline_strat_tol, &
395 0.5*(1.0 + cs%halocline_strat_tol))
396 if (cs%halocline_strat_tol > 0.0)
then
398 i_hstol = 0.0 ;
if (fn_zero_val - cs%halocline_strat_tol > 0.0) &
399 i_hstol = 1.0 / (fn_zero_val - cs%halocline_strat_tol)
400 do k=nz,1,-1 ;
if (cs%ref_pressure > p_is(k+1))
then
401 z_int_unst = z_int_unst + fn_now * h_col(k)
402 if (strat_rat(k) <= fn_zero_val)
then
403 if (strat_rat(k) <= cs%halocline_strat_tol)
then ; fn_now = 1.0
405 fn_now = max(fn_now, (fn_zero_val - strat_rat(k)) * i_hstol)
410 do k=nz,1,-1 ;
if (cs%ref_pressure > p_is(k+1))
then
411 z_int_unst = z_int_unst + fn_now * h_col(k)
412 if (strat_rat(k) <= cs%halocline_strat_tol) fn_now = 1.0
416 if (z_interior < z_int_unst)
then
418 kur1 = max(int(ceiling(k_interior)),2)
419 if (z_col_new(kur1-1) < z_interior)
then
421 do k = kur1,nz+1 ;
if (z_col_new(k) >= z_int_unst)
then
423 if (z_col_new(k-1) >= z_int_unst) &
424 call mom_error(fatal,
"build_grid_SLight, bad halocline structure.")
425 k_int2 = real(k-1) + (z_int_unst - z_col_new(k-1)) / &
426 (z_col_new(k) - z_col_new(k-1))
429 if (z_col_new(nz+1) < z_int_unst)
then
431 z_int_unst = z_col_new(nz+1) ; k_int2 = real(nz+1)
435 if (k_int2 > k_interior)
then
436 k_interior = k_int2 ; z_interior = z_int_unst
444 z_col_new(k) = min((k-1)*cs%dz_ml_min, &
445 z_col_new(nz+1) - cs%min_thickness*(nz+1-k))
447 z_ml_fix = z_col_new(nkml+1)
448 if (z_interior > z_ml_fix)
then
449 dz_dk = (z_interior - z_ml_fix) / (k_interior - (nkml+1))
450 do k=nkml+2,int(floor(k_interior))
451 z_col_new(k) = z_ml_fix + dz_dk * (k - (nkml+1))
455 if (z_col_new(k) <= z_col_new(cs%nz_fixed_surface+1))
then
456 z_col_new(k) = z_col_new(cs%nz_fixed_surface+1)
461 if (maximum_depths_set .and. maximum_h_set)
then ;
do k=2,nz
464 z_col_new(k) = min(z_col_new(k), cs%max_interface_depths(k), &
465 z_col_new(k-1) + cs%max_layer_thickness(k-1))
466 enddo ;
elseif (maximum_depths_set)
then ;
do k=2,nz
467 z_col_new(k) = min(z_col_new(k), cs%max_interface_depths(k))
468 enddo ;
elseif (maximum_h_set)
then ;
do k=2,nz
469 z_col_new(k) = min(z_col_new(k), z_col_new(k-1) + cs%max_layer_thickness(k-1))