12 implicit none ;
private
28 real :: compressibility_fraction
32 real :: rho_ml_avg_depth
35 real :: nlay_ml_offset
38 integer :: nz_fixed_surface = 2
45 logical :: fix_haloclines = .false.
49 real :: halocline_filter_length
52 real :: halocline_strat_tol
55 real,
allocatable,
dimension(:) :: target_density
61 real,
allocatable,
dimension(:) :: max_interface_depths
64 real,
allocatable,
dimension(:) :: max_layer_thickness
75 subroutine init_coord_slight(CS, nk, ref_pressure, target_density, interp_CS, m_to_H, rho_scale)
77 integer,
intent(in) :: nk
78 real,
intent(in) :: ref_pressure
79 real,
dimension(:),
intent(in) :: target_density
81 real,
optional,
intent(in) :: m_to_h
82 real,
optional,
intent(in) :: rho_scale
84 real :: m_to_h_rescale
86 if (
associated(cs))
call mom_error(fatal,
"init_coord_slight: CS already associated!")
88 allocate(cs%target_density(nk+1))
90 m_to_h_rescale = 1.0 ;
if (
present(m_to_h)) m_to_h_rescale = m_to_h
93 cs%ref_pressure = ref_pressure
94 cs%target_density(:) = target_density(:)
95 cs%interp_CS = interp_cs
98 cs%compressibility_fraction = 0.
99 cs%Rho_ML_avg_depth = 1.0 * m_to_h_rescale
100 cs%nlay_ml_offset = 2.0
101 cs%dz_ml_min = 1.0 * m_to_h_rescale
102 cs%halocline_filter_length = 2.0 * m_to_h_rescale
103 cs%halocline_strat_tol = 0.25
104 cs%kg_m3_to_R = 1.0 ;
if (
present(rho_scale)) cs%kg_m3_to_R = rho_scale
113 if (.not.
associated(cs))
return
114 deallocate(cs%target_density)
120 min_thickness, compressibility_fraction, dz_ml_min, &
121 nz_fixed_surface, Rho_ML_avg_depth, nlay_ML_offset, fix_haloclines, &
122 halocline_filter_length, halocline_strat_tol, interp_CS)
124 real,
dimension(:), &
125 optional,
intent(in) :: max_interface_depths
126 real,
dimension(:), &
127 optional,
intent(in) :: max_layer_thickness
128 real,
optional,
intent(in) :: min_thickness
130 real,
optional,
intent(in) :: compressibility_fraction
133 real,
optional,
intent(in) :: dz_ml_min
135 integer,
optional,
intent(in) :: nz_fixed_surface
137 real,
optional,
intent(in) :: rho_ml_avg_depth
139 real,
optional,
intent(in) :: nlay_ml_offset
141 logical,
optional,
intent(in) :: fix_haloclines
143 real,
optional,
intent(in) :: halocline_filter_length
145 real,
optional,
intent(in) :: halocline_strat_tol
148 optional,
intent(in) :: interp_cs
150 if (.not.
associated(cs))
call mom_error(fatal,
"set_slight_params: CS not associated")
152 if (
present(max_interface_depths))
then
153 if (
size(max_interface_depths) /= cs%nk+1) &
154 call mom_error(fatal,
"set_slight_params: max_interface_depths inconsistent size")
155 allocate(cs%max_interface_depths(cs%nk+1))
156 cs%max_interface_depths(:) = max_interface_depths(:)
159 if (
present(max_layer_thickness))
then
160 if (
size(max_layer_thickness) /= cs%nk) &
161 call mom_error(fatal,
"set_slight_params: max_layer_thickness inconsistent size")
162 allocate(cs%max_layer_thickness(cs%nk))
163 cs%max_layer_thickness(:) = max_layer_thickness(:)
166 if (
present(min_thickness)) cs%min_thickness = min_thickness
167 if (
present(compressibility_fraction)) cs%compressibility_fraction = compressibility_fraction
169 if (
present(dz_ml_min)) cs%dz_ml_min = dz_ml_min
170 if (
present(nz_fixed_surface)) cs%nz_fixed_surface = nz_fixed_surface
171 if (
present(rho_ml_avg_depth)) cs%Rho_ML_avg_depth = rho_ml_avg_depth
172 if (
present(nlay_ml_offset)) cs%nlay_ML_offset = nlay_ml_offset
173 if (
present(fix_haloclines)) cs%fix_haloclines = fix_haloclines
174 if (
present(halocline_filter_length)) cs%halocline_filter_length = halocline_filter_length
175 if (
present(halocline_strat_tol))
then
176 if (halocline_strat_tol > 1.0)
call mom_error(fatal,
"set_slight_params: "//&
177 "HALOCLINE_STRAT_TOL must not exceed 1.0.")
178 cs%halocline_strat_tol = halocline_strat_tol
181 if (
present(interp_cs)) cs%interp_CS = interp_cs
186 nz, depth, h_col, T_col, S_col, p_col, z_col, z_col_new, &
187 h_neglect, h_neglect_edge)
189 type(
eos_type),
pointer :: eqn_of_state
190 real,
intent(in) :: h_to_pa
191 real,
intent(in) :: h_subroundoff
192 integer,
intent(in) :: nz
193 real,
intent(in) :: depth
194 real,
dimension(nz),
intent(in) :: t_col
195 real,
dimension(nz),
intent(in) :: s_col
196 real,
dimension(nz),
intent(in) :: h_col
197 real,
dimension(nz),
intent(in) :: p_col
198 real,
dimension(nz+1),
intent(in) :: z_col
199 real,
dimension(nz+1),
intent(inout) :: z_col_new
200 real,
optional,
intent(in) :: h_neglect
202 real,
optional,
intent(in) :: h_neglect_edge
205 real,
dimension(nz) :: rho_col
206 real,
dimension(nz) :: t_f, s_f
207 logical,
dimension(nz+1) :: reliable
208 real,
dimension(nz+1) :: t_int, s_int
209 real,
dimension(nz+1) :: rho_tmp
210 real,
dimension(nz+1) :: drho_dp
211 real,
dimension(nz+1) :: p_is, p_r
212 real,
dimension(nz+1) :: drhois_dt
214 real,
dimension(nz+1) :: drhois_ds
216 real,
dimension(nz+1) :: drhor_dt
218 real,
dimension(nz+1) :: drhor_ds
220 real,
dimension(nz+1) :: strat_rat
223 real :: fn_now, i_hstol, fn_zero_val
239 logical :: maximum_depths_set
240 logical :: maximum_h_set
241 real :: k2_used, k2here, dz_sum, z_max
243 real :: h_tr, b_denom_1, b1, d1
244 real,
dimension(nz) :: c1
245 integer :: kur1, kur2
247 integer :: i, j, k, nkml
249 maximum_depths_set =
allocated(cs%max_interface_depths)
250 maximum_h_set =
allocated(cs%max_layer_thickness)
252 if (z_col(nz+1) - z_col(1) < nz*cs%min_thickness)
then
254 dz = (z_col(nz+1) - z_col(1)) / real(nz)
255 do k=2,nz ; z_col_new(k) = z_col(1) + dz*real(k-1) ;
enddo
258 eqn_of_state, scale=cs%kg_m3_to_R)
263 z_col_new, cs, reliable, debug=.true., &
264 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
267 if (cs%min_thickness > 0.0)
then
269 do k=2,nz ;
if (z_col_new(k) < z_col_new(k-1) + cs%min_thickness)
then
270 z_col_new(k) = z_col_new(k-1) + cs%min_thickness
273 do k=nz,2,-1 ;
if (z_col_new(k) > z_col_new(k+1) - cs%min_thickness)
then
274 z_col_new(k) = z_col_new(k+1) - cs%min_thickness
285 do k=kur_ss,nz ;
if (.not.reliable(k))
then
291 do k=kur1+1,nz+1 ;
if (reliable(k))
then
292 kur2 = k-1 ; kur_ss = k ;
exit
294 if (kur2 < kur1)
call mom_error(fatal,
"Bad unreliable range.")
296 dz_ur = z_col_new(kur2+1) - z_col_new(kur1-1)
300 wgt = 1.0 ; cowgt = 0.0
302 z_col_new(k) = cowgt*z_col_new(k) + &
303 wgt * (z_col_new(kur1-1) + dz_ur*(k - (kur1-1)) / ((kur2 - kur1) + 2))
309 z_wt = 0.0 ; rho_x_z = 0.0
310 h_ml_av = cs%Rho_ml_avg_depth
312 if (z_wt + h_col(k) >= h_ml_av)
then
313 rho_x_z = rho_x_z + rho_col(k) * (h_ml_av - z_wt)
317 rho_x_z = rho_x_z + rho_col(k) * h_col(k)
318 z_wt = z_wt + h_col(k)
321 if (z_wt > 0.0) rho_ml_av = rho_x_z / z_wt
323 nkml = cs%nz_fixed_surface
325 if (rho_ml_av <= cs%target_density(nkml))
then
326 k_interior = cs%nlay_ml_offset + real(nkml)
327 elseif (rho_ml_av > cs%target_density(nz+1))
then
328 k_interior = real(nz+1)
330 if ((rho_ml_av >= cs%target_density(k)) .and. &
331 (rho_ml_av < cs%target_density(k+1)))
then
332 k_interior = (cs%nlay_ml_offset + k) + &
333 (rho_ml_av - cs%target_density(k)) / &
334 (cs%target_density(k+1) - cs%target_density(k))
338 if (k_interior > real(nz+1)) k_interior = real(nz+1)
341 k = int(ceiling(k_interior))
342 z_interior = (k-k_interior)*z_col_new(k-1) + (1.0+(k_interior-k))*z_col_new(k)
344 if (cs%fix_haloclines)
then
350 if (cs%halocline_filter_length > 0.0)
then
351 lfilt = cs%halocline_filter_length
354 h_tr = h_col(1) + h_subroundoff
355 b1 = 1.0 / (h_tr + lfilt) ; d1 = h_tr * b1
356 t_f(1) = (b1*h_tr)*t_col(1) ; s_f(1) = (b1*h_tr)*s_col(1)
359 h_tr = h_col(k) + h_subroundoff ; b_denom_1 = h_tr + d1*lfilt
360 b1 = 1.0 / (b_denom_1 + lfilt) ; d1 = b_denom_1 * b1
361 t_f(k) = b1 * (h_tr*t_col(k) + lfilt*t_f(k-1))
362 s_f(k) = b1 * (h_tr*s_col(k) + lfilt*s_f(k-1))
365 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)
368 do k=1,nz ; t_f(k) = t_col(k) ; s_f(k) = s_col(k) ;
enddo
371 t_int(1) = t_f(1) ; s_int(1) = s_f(1)
373 t_int(k) = 0.5*(t_f(k-1) + t_f(k)) ; s_int(k) = 0.5*(s_f(k-1) + s_f(k))
374 p_is(k) = z_col(k) * h_to_pa
375 p_r(k) = cs%ref_pressure + cs%compressibility_fraction * ( p_is(k) - cs%ref_pressure )
377 t_int(nz+1) = t_f(nz) ; s_int(nz+1) = s_f(nz)
378 p_is(nz+1) = z_col(nz+1) * h_to_pa
380 eqn_of_state, scale=cs%kg_m3_to_R)
382 eqn_of_state, scale=cs%kg_m3_to_R)
383 if (cs%compressibility_fraction > 0.0)
then
387 do k=2,nz ; drho_dp(k) = 0.0 ;
enddo
390 h_to_cpa = cs%compressibility_fraction*cs%kg_m3_to_R*h_to_pa
393 dris = drhois_dt(k) * (t_f(k) - t_f(k-1)) + &
394 drhois_ds(k) * (s_f(k) - s_f(k-1))
395 drr = (drhor_dt(k) * (t_f(k) - t_f(k-1)) + &
396 drhor_ds(k) * (s_f(k) - s_f(k-1))) + &
397 drho_dp(k) * (h_to_cpa*0.5*(h_col(k) + h_col(k-1)))
399 if (dris <= 0.0)
then
402 strat_rat(k) = 2.0*max(drr,0.0) / (dris + abs(drr))
405 strat_rat(nz+1) = 1.0
407 z_int_unst = 0.0 ; fn_now = 0.0
408 fn_zero_val = min(2.0*cs%halocline_strat_tol, &
409 0.5*(1.0 + cs%halocline_strat_tol))
410 if (cs%halocline_strat_tol > 0.0)
then
412 i_hstol = 0.0 ;
if (fn_zero_val - cs%halocline_strat_tol > 0.0) &
413 i_hstol = 1.0 / (fn_zero_val - cs%halocline_strat_tol)
414 do k=nz,1,-1 ;
if (cs%ref_pressure > p_is(k+1))
then
415 z_int_unst = z_int_unst + fn_now * h_col(k)
416 if (strat_rat(k) <= fn_zero_val)
then
417 if (strat_rat(k) <= cs%halocline_strat_tol)
then ; fn_now = 1.0
419 fn_now = max(fn_now, (fn_zero_val - strat_rat(k)) * i_hstol)
424 do k=nz,1,-1 ;
if (cs%ref_pressure > p_is(k+1))
then
425 z_int_unst = z_int_unst + fn_now * h_col(k)
426 if (strat_rat(k) <= cs%halocline_strat_tol) fn_now = 1.0
430 if (z_interior < z_int_unst)
then
432 kur1 = max(int(ceiling(k_interior)),2)
433 if (z_col_new(kur1-1) < z_interior)
then
435 do k = kur1,nz+1 ;
if (z_col_new(k) >= z_int_unst)
then
437 if (z_col_new(k-1) >= z_int_unst) &
438 call mom_error(fatal,
"build_grid_SLight, bad halocline structure.")
439 k_int2 = real(k-1) + (z_int_unst - z_col_new(k-1)) / &
440 (z_col_new(k) - z_col_new(k-1))
443 if (z_col_new(nz+1) < z_int_unst)
then
445 z_int_unst = z_col_new(nz+1) ; k_int2 = real(nz+1)
449 if (k_int2 > k_interior)
then
450 k_interior = k_int2 ; z_interior = z_int_unst
458 z_col_new(k) = min((k-1)*cs%dz_ml_min, &
459 z_col_new(nz+1) - cs%min_thickness*(nz+1-k))
461 z_ml_fix = z_col_new(nkml+1)
462 if (z_interior > z_ml_fix)
then
463 dz_dk = (z_interior - z_ml_fix) / (k_interior - (nkml+1))
464 do k=nkml+2,int(floor(k_interior))
465 z_col_new(k) = z_ml_fix + dz_dk * (k - (nkml+1))
469 if (z_col_new(k) <= z_col_new(cs%nz_fixed_surface+1))
then
470 z_col_new(k) = z_col_new(cs%nz_fixed_surface+1)
475 if (maximum_depths_set .and. maximum_h_set)
then ;
do k=2,nz
478 z_col_new(k) = min(z_col_new(k), cs%max_interface_depths(k), &
479 z_col_new(k-1) + cs%max_layer_thickness(k-1))
480 enddo ;
elseif (maximum_depths_set)
then ;
do k=2,nz
481 z_col_new(k) = min(z_col_new(k), cs%max_interface_depths(k))
482 enddo ;
elseif (maximum_h_set)
then ;
do k=2,nz
483 z_col_new(k) = min(z_col_new(k), z_col_new(k-1) + cs%max_layer_thickness(k-1))
493 CS, reliable, debug, h_neglect, h_neglect_edge)
494 integer,
intent(in) :: nz
495 real,
dimension(nz),
intent(in) :: rho_col
496 real,
dimension(nz),
intent(in) :: h_col
497 real,
dimension(nz+1),
intent(in) :: z_col
498 real,
dimension(nz+1),
intent(in) :: rho_tgt
499 real,
dimension(nz+1),
intent(inout) :: z_col_new
501 logical,
dimension(nz+1),
intent(inout) :: reliable
503 logical,
optional,
intent(in) :: debug
504 real,
optional,
intent(in) :: h_neglect
507 real,
optional,
intent(in) :: h_neglect_edge
511 real,
dimension(nz+1) :: ru_max_int
512 real,
dimension(nz+1) :: ru_min_int
513 real,
dimension(nz) :: ru_max_lay
514 real,
dimension(nz) :: ru_min_lay
515 real,
dimension(nz,2) :: ppoly_i_E
516 real,
dimension(nz,2) :: ppoly_i_S
517 real,
dimension(nz,DEGREE_MAX+1) :: ppoly_i_coefficients
518 logical,
dimension(nz) :: unstable_lay
519 logical,
dimension(nz+1) :: unstable_int
524 real :: zf1, zf2, rfn1, rfn2
525 real :: drfn_dzf, sgn, delta_zf, zf_prev
529 integer :: ppoly_degree
530 integer :: k, k1, k1_min, itt, max_itt, m
535 debugging = .false. ;
if (
present(debug)) debugging = debug
539 z_sgn = 1.0 ;
if ( z_col(1) > z_col(nz+1) ) z_sgn = -1.0
542 if (abs((z_col(k+1) - z_col(k)) - z_sgn*h_col(k)) > &
543 1.0e-14*(abs(z_col(k+1)) + abs(z_col(k)) + abs(h_col(k))) ) &
544 call mom_error(fatal,
"rho_interfaces_col: Inconsistent z_col and h_col")
548 if ( z_col(1) == z_col(nz+1) )
then
550 do k=1,nz+1 ; z_col_new(k) = z_col(1) ; reliable(k) = .true. ;
enddo
555 call regridding_set_ppolys(cs%interp_CS, rho_col, nz, h_col, ppoly_i_e, ppoly_i_s, &
556 ppoly_i_coefficients, ppoly_degree, h_neglect, h_neglect_edge)
562 unstable_int(1) = .false.
563 ru_max_int(1) = ppoly_i_e(1,1)
565 unstable_lay(1) = (ppoly_i_e(1,1) > ppoly_i_e(1,2))
566 ru_max_lay(1) = max(ppoly_i_e(1,1), ppoly_i_e(1,2))
569 unstable_int(k) = (ppoly_i_e(k-1,2) > ppoly_i_e(k,1))
570 ru_max_int(k) = max(ppoly_i_e(k-1,2), ppoly_i_e(k,1))
571 ru_min_int(k) = min(ppoly_i_e(k-1,2), ppoly_i_e(k,1))
572 if (unstable_int(k) .and. unstable_lay(k-1)) &
573 ru_max_int(k) = max(ru_max_lay(k-1), ru_max_int(k))
575 unstable_lay(k) = (ppoly_i_e(k,1) > ppoly_i_e(k,2))
576 ru_max_lay(k) = max(ppoly_i_e(k,1), ppoly_i_e(k,2))
577 ru_min_lay(k) = min(ppoly_i_e(k,1), ppoly_i_e(k,2))
578 if (unstable_lay(k) .and. unstable_int(k)) &
579 ru_max_lay(k) = max(ru_max_int(k), ru_max_lay(k))
581 unstable_int(nz+1) = .false.
582 ru_min_int(nz+1) = ppoly_i_e(nz,2)
585 if (unstable_lay(k) .and. unstable_int(k+1)) &
586 ru_min_lay(k) = min(ru_min_int(k+1), ru_min_lay(k))
588 if (unstable_int(k) .and. unstable_lay(k)) &
589 ru_min_int(k) = min(ru_min_lay(k), ru_min_int(k))
592 z_col_new(1) = z_col(1) ; reliable(1) = .true.
600 if (rt <= ppoly_i_e(k1_min,1))
then
601 z_col_new(k) = z_col(k1_min)
604 elseif (k1_min == nz+1)
then
605 z_col_new(k) = z_col(nz+1)
608 if (unstable_int(k) .and. (rt >= ru_min_int(k)) .and. (rt <= ru_max_int(k)))
then
610 z_col_new(k) = z_col(k) ; reliable(k) = .false.
611 k1_min = k ; k_found = .true.
612 elseif ((rt >= ppoly_i_e(k-1,2)) .and. (rt <= ppoly_i_e(k,1)))
then
614 z_col_new(k) = z_col(k) ; reliable(k) = .true.
615 k1_min = k ; k_found = .true.
616 elseif (rt < ppoly_i_e(k-1,2))
then
619 if ((rt < ppoly_i_e(k1,2)) .and. (rt > ppoly_i_e(k1,1)))
then
622 k1_min = k1 ; k_found = .true. ;
exit
623 elseif (unstable_lay(k1) .and. (rt >= ru_min_lay(k1)) .and. (rt <= ru_max_lay(k1)))
then
626 z_col_new(k) = z_col(k1+1) ; reliable(k) = .false.
627 k1_min = k1 ; k_found = .true. ;
exit
630 if (k1 > 1)
then ;
if ((rt <= ppoly_i_e(k1,1)) .and. (rt >= ppoly_i_e(k1-1,2)))
then
632 z_col_new(k) = z_col(k1) ; reliable(k) = .true.
633 k1_min = k1 ; k_found = .true. ;
exit
634 elseif (unstable_int(k1) .and. (rt >= ru_min_int(k1)) .and. (rt <= ru_max_int(k1)))
then
637 z_col_new(k) = z_col(k1) ; reliable(k) = .false.
638 k1_min = k1 ; k_found = .true. ;
exit
642 if (.not.k_found)
then
645 z_col_new(k) = z_col(k1_min)
647 z_col_new(k) = z_col(k1_min)
653 if ((rt < ppoly_i_e(k1,2)) .and. (rt > ppoly_i_e(k1,1)))
then
656 k1_min = k1 ; k_found = .true. ;
exit
657 elseif (unstable_lay(k1) .and. (rt >= ru_min_lay(k1)) .and. (rt <= ru_max_lay(k1)))
then
660 z_col_new(k) = z_col(k1)
661 reliable(k) = .false.
662 k1_min = k1 ; k_found = .true. ;
exit
664 if (k1 < nz)
then ;
if ((rt <= ppoly_i_e(k1+1,1)) .and. (rt >= ppoly_i_e(k1,2)))
then
667 z_col_new(k) = z_col(k1+1) ; reliable(k) = .true.
668 k1_min = k1+1 ; k_found = .true. ;
exit
669 elseif (unstable_int(k1+1) .and. (rt >= ru_min_int(k1+1)) .and. (rt <= ru_max_int(k1+1)))
then
672 z_col_new(k) = z_col(k1+1)
673 reliable(k) = .false.
674 k1_min = k1+1 ; k_found = .true. ;
exit
677 if (.not.k_found)
then
678 z_col_new(k) = z_col(nz+1)
679 if (rt >= ppoly_i_e(nz,2))
then
682 reliable(k) = .false.
687 if (k_layer > 0)
then
689 if (.not.(ppoly_i_e(k1,2) > ppoly_i_e(k1,1)))
call mom_error(fatal, &
690 "build_grid_SLight: Erroneously searching for an interface in an unstratified layer.")
693 zf = (rt - ppoly_i_e(k1,1)) / (ppoly_i_e(k1,2) - ppoly_i_e(k1,1))
695 if (ppoly_degree > 1)
then
696 a(:) = 0.0 ; a(1) = ppoly_i_coefficients(k_layer,1) - rt
697 do m=2,ppoly_degree+1 ; a(m) = ppoly_i_coefficients(k_layer,m) ;
enddo
699 zf1 = 0.0 ; rfn1 = a(1)
700 zf2 = 1.0 ; rfn2 = a(1) + (a(2) + (a(3) + (a(4) + a(5))))
701 if (rfn1 * rfn2 > 0.0)
call mom_error(fatal,
"build_grid_SLight: Bad bracketing.")
704 rfn = a(1) + zf*(a(2) + zf*(a(3) + zf*(a(4) + zf*a(5))))
706 if (rfn * rfn1 > 0.0)
then
707 zf1 = zf ; rfn1 = rfn
709 zf2 = zf ; rfn2 = rfn
711 if (rfn1 == rfn2)
exit
713 drfn_dzf = (a(2) + zf*(2.0*a(3) + zf*(3.0*a(4) + zf*4.0*a(5))))
714 sgn = 1.0 ;
if (drfn_dzf < 0.0) sgn = -1.0
716 if ((sgn*(zf - rfn) >= zf1 * abs(drfn_dzf)) .and. &
717 (sgn*(zf - rfn) <= zf2 * abs(drfn_dzf)))
then
718 delta_zf = -rfn / drfn_dzf
722 zf = ( rfn2 * zf1 - rfn1 * zf2 ) / (rfn2 - rfn1)
723 delta_zf = zf - zf_prev
726 if (abs(delta_zf) < tol)
exit
729 z_col_new(k) = z_col(k_layer) + zf * z_sgn * h_col(k_layer)
736 z_col_new(nz+1) = z_col(nz+1) ; reliable(nz+1) = .true.