11 implicit none ;
private
20 real :: min_thickness = 0.
27 logical :: integrate_downward_for_e = .false.
30 real,
allocatable,
dimension(:) :: target_density
49 subroutine init_coord_rho(CS, nk, ref_pressure, target_density, interp_CS, rho_scale)
51 integer,
intent(in) :: nk
52 real,
intent(in) :: ref_pressure
53 real,
dimension(:),
intent(in) :: target_density
55 real,
optional,
intent(in) :: rho_scale
57 if (
associated(cs))
call mom_error(fatal,
"init_coord_rho: CS already associated!")
59 allocate(cs%target_density(nk+1))
62 cs%ref_pressure = ref_pressure
63 cs%target_density(:) = target_density(:)
64 cs%interp_CS = interp_cs
65 cs%kg_m3_to_R = 1.0 ;
if (
present(rho_scale)) cs%kg_m3_to_R = rho_scale
74 if (.not.
associated(cs))
return
75 deallocate(cs%target_density)
80 subroutine set_rho_params(CS, min_thickness, integrate_downward_for_e, interp_CS)
82 real,
optional,
intent(in) :: min_thickness
83 logical,
optional,
intent(in) :: integrate_downward_for_e
89 if (.not.
associated(cs))
call mom_error(fatal,
"set_rho_params: CS not associated")
91 if (
present(min_thickness)) cs%min_thickness = min_thickness
92 if (
present(integrate_downward_for_e)) cs%integrate_downward_for_e = integrate_downward_for_e
93 if (
present(interp_cs)) cs%interp_CS = interp_cs
100 subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, &
101 h_neglect, h_neglect_edge)
103 integer,
intent(in) :: nz
104 real,
intent(in) :: depth
105 real,
dimension(nz),
intent(in) :: h
106 real,
dimension(nz),
intent(in) :: t
107 real,
dimension(nz),
intent(in) :: s
108 type(
eos_type),
pointer :: eqn_of_state
109 real,
dimension(CS%nk+1), &
110 intent(inout) :: z_interface
111 real,
optional,
intent(in) :: h_neglect
114 real,
optional,
intent(in) :: h_neglect_edge
118 integer :: k, count_nonzero_layers
119 integer,
dimension(nz) :: mapping
120 real,
dimension(nz) :: p, h_nv
121 real,
dimension(nz) :: densities
122 real,
dimension(nz+1) :: xtmp
123 real,
dimension(CS%nk) :: h_new
124 real,
dimension(CS%nk+1) :: x1
129 if (count_nonzero_layers > 1)
then
131 do k = 1,count_nonzero_layers
132 xtmp(k+1) = xtmp(k) + h_nv(k)
136 p(:) = cs%ref_pressure
137 call calculate_density(t, s, p, densities, 1, nz, eqn_of_state, scale=cs%kg_m3_to_R)
138 do k = 1,count_nonzero_layers
139 densities(k) = densities(mapping(k))
144 h_nv, xtmp, cs%target_density, cs%nk, h_new, &
145 x1, h_neglect, h_neglect_edge)
151 x1(1) = 0.0 ;
do k = 1,cs%nk ; x1(k+1) = x1(k) + h_new(k) ;
enddo
153 h_new(k) = x1(k+1) - x1(k)
157 if (nz == cs%nk)
then
166 if (cs%integrate_downward_for_e)
then
170 z_interface(k+1) = z_interface(k) - h_new(k)
174 z_interface(cs%nk+1) = -depth
176 z_interface(k) = z_interface(k+1) + h_new(k)
195 zInterface, h_neglect, h_neglect_edge)
198 integer,
intent(in) :: nz
199 real,
intent(in) :: depth
200 real,
dimension(nz),
intent(in) :: h
201 real,
dimension(nz),
intent(in) :: T
202 real,
dimension(nz),
intent(in) :: S
203 type(
eos_type),
pointer :: eqn_of_state
204 real,
dimension(nz+1),
intent(inout) :: zInterface
205 real,
optional,
intent(in) :: h_neglect
208 real,
optional,
intent(in) :: h_neglect_edge
213 integer :: count_nonzero_layers
218 real,
dimension(nz) :: p, densities, T_tmp, S_tmp, Tmp
219 integer,
dimension(nz) :: mapping
220 real,
dimension(nz) :: h0, h1, hTmp
221 real,
dimension(nz+1) :: x0, x1, xTmp
223 threshold = cs%min_thickness
224 p(:) = cs%ref_pressure
237 if ( count_nonzero_layers <= 1 )
then
243 do k = 1,count_nonzero_layers
244 xtmp(k+1) = xtmp(k) + htmp(k)
249 1, nz, eqn_of_state, scale=cs%kg_m3_to_R)
251 do k = 1,count_nonzero_layers
252 densities(k) = densities(mapping(k))
258 htmp, xtmp, cs%target_density, nz, h1, x1, h_neglect, h_neglect_edge)
261 x1(1) = 0.0 ;
do k = 1,nz ; x1(k+1) = x1(k) + h1(k) ;
enddo
265 h1(k) = x1(k+1) - x1(k)
268 call remapping_core_h(remapcs, nz, h0, s, nz, h1, tmp, h_neglect, h_neglect_edge)
271 call remapping_core_h(remapcs, nz, h0, t, nz, h1, tmp, h_neglect, h_neglect_edge)
279 x0(k) = x0(k-1) + h0(k-1)
280 x1(k) = x1(k-1) + h1(k-1)
281 deviation = deviation + (x0(k)-x1(k))**2
283 deviation = sqrt( deviation / (nz-1) )
292 if (cs%integrate_downward_for_e)
then
295 zinterface(k+1) = zinterface(k) - h1(k)
301 zinterface(nz+1) = -depth
303 zinterface(k) = zinterface(k+1) + h1(k)
313 integer,
intent(in) :: nk
314 real,
dimension(nk),
intent(in) :: h_in
315 real,
intent(in) :: threshold
316 integer,
intent(out) :: nout
317 real,
dimension(nk),
intent(out) :: h_out
318 integer,
dimension(nk),
intent(out) :: mapping
320 integer :: k, k_thickest
321 real :: thickness_in_vanished, thickest_h_out
325 thickness_in_vanished = 0.0
326 thickest_h_out = h_in(1)
331 if (h_in(k) > threshold)
then
335 h_out(nout) = h_in(k)
336 if (h_out(nout) > thickest_h_out)
then
337 thickest_h_out = h_out(nout)
342 thickness_in_vanished = thickness_in_vanished + h_in(k)
347 if (nout <= 1)
return
350 h_out(k_thickest) = h_out(k_thickest) + thickness_in_vanished
359 real,
intent(in) :: min_thickness
360 integer,
intent(in) :: nk
361 real,
dimension(:),
intent(inout) :: h
366 integer :: count_nonzero_layers
372 count_nonzero_layers = 0
374 if ( h(k) > min_thickness )
then
375 count_nonzero_layers = count_nonzero_layers + 1
380 if ( count_nonzero_layers == nk )
return
383 if ( count_nonzero_layers == 0 )
then
393 if ( h(k) <= min_thickness )
then
394 delta = min_thickness - h(k)
395 correction = correction + delta
404 if ( h(k) > maxthickness )
then
410 h(k_found) = h(k_found) - correction