11 implicit none ;
private
13 #include <MOM_memory.h>
20 subroutine full_convection(G, GV, h, tv, T_adj, S_adj, p_surf, Kddt_smooth, &
24 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
28 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
30 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
32 real,
dimension(:,:),
pointer :: p_surf
33 real,
intent(in) :: kddt_smooth
35 real,
optional,
intent(in) :: kddt_convect
37 integer,
optional,
intent(in) :: halo
40 real,
dimension(SZI_(G),SZK_(G)+1) :: &
46 real,
dimension(SZI_(G),SZK0_(G)) :: &
55 real,
dimension(SZI_(G),SZK_(G)+1) :: &
64 real,
dimension(SZI_(G),SZK_(G)+1) :: &
75 real,
dimension(SZI_(G),SZK_(G)+1) :: &
83 logical,
dimension(SZI_(G)) :: do_i
84 logical,
dimension(SZI_(G)) :: last_down
85 integer,
dimension(SZI_(G)) :: change_ct
87 integer :: changed_col
88 integer :: i, j, k, is, ie, js, je, nz, itt
90 if (
present(halo))
then
91 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
93 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
97 if (.not.
associated(tv%eqn_of_state))
return
99 h_neglect = gv%H_subroundoff
101 if (
present(kddt_convect)) kap_dt_x2 = 2.0*kddt_convect
102 mix_len = (1.0e20 * nz) * (g%max_depth * gv%Z_to_H)
103 h0 = 1.0e-16*sqrt(kddt_smooth) + h_neglect
106 mix(:,:) = 0.0 ; d_b(:,:) = 1.0
108 te_b(:,:) = 0.0 ; se_b(:,:) = 0.0
113 do_i(i) = (g%mask2dT(i,j) > 0.0)
116 last_down(i) = .true.
118 te_a(i,0) = 0.0 ; se_a(i,0) = 0.0
124 do i=is,ie ; change_ct(i) = 0 ;
enddo
127 do k=2,nz ;
do i=is,ie ;
if (do_i(i))
then
129 h_a = h(i,j,k-1) + h_neglect ; h_b = h(i,j,k) + h_neglect
130 if (mix(i,k) <= 0.0)
then
131 if (
is_unstable(drho_dt(i,k), drho_ds(i,k), h_a, h_b, mix(i,k-1), mix(i,k+1), &
132 tv%T(i,j,k-1), tv%T(i,j,k), tv%S(i,j,k-1), tv%S(i,j,k), &
133 te_a(i,k-2), te_b(i,k+1), se_a(i,k-2), se_b(i,k+1), &
134 d_a(i,k-1), d_b(i,k+1)))
then
136 if (kap_dt_x2 > 0.0) mix(i,k) = kap_dt_x2 / ((h(i,j,k-1)+h(i,j,k)) + h0)
137 change_ct(i) = change_ct(i) + 1
141 b_a = 1.0 / ((h_a + d_a(i,k-1)*mix(i,k-1)) + mix(i,k))
142 if (mix(i,k) <= 0.0)
then
143 c_a(i,k) = 0.0 ; d_a(i,k) = 1.0
145 d_a(i,k) = b_a * (h_a + d_a(i,k-1)*mix(i,k-1))
146 c_a(i,k) = 1.0 ;
if (d_a(i,k) > epsilon(b_a)) c_a(i,k) = b_a * mix(i,k)
150 te_a(i,k-1) = b_a * (h_a*tv%T(i,j,k-1) + mix(i,k-1)*te_a(i,k-2))
151 se_a(i,k-1) = b_a * (h_a*tv%S(i,j,k-1) + mix(i,k-1)*se_a(i,k-2))
153 te_a(i,k-1) = b_a * (h_a*tv%T(i,j,k-1))
154 se_a(i,k-1) = b_a * (h_a*tv%S(i,j,k-1))
156 endif ;
enddo ;
enddo
160 do i=is,ie ;
if (do_i(i))
then
161 if (change_ct(i) == 0)
then
162 last_down(i) = .true. ; do_i(i) = .false.
164 changed_col = changed_col + 1 ; change_ct(i) = 0
167 if (changed_col == 0)
exit
170 do k=nz,2,-1 ;
do i=is,ie ;
if (do_i(i))
then
172 h_a = h(i,j,k-1) + h_neglect ; h_b = h(i,j,k) + h_neglect
173 if (mix(i,k) <= 0.0)
then
174 if (
is_unstable(drho_dt(i,k), drho_ds(i,k), h_a, h_b, mix(i,k-1), mix(i,k+1), &
175 tv%T(i,j,k-1), tv%T(i,j,k), tv%S(i,j,k-1), tv%S(i,j,k), &
176 te_a(i,k-2), te_b(i,k+1), se_a(i,k-2), se_b(i,k+1), &
177 d_a(i,k-1), d_b(i,k+1)))
then
179 if (kap_dt_x2 > 0.0) mix(i,k) = kap_dt_x2 / ((h(i,j,k-1)+h(i,j,k)) + h0)
180 change_ct(i) = change_ct(i) + 1
184 b_b = 1.0 / ((h_b + d_b(i,k+1)*mix(i,k+1)) + mix(i,k))
185 if (mix(i,k) <= 0.0)
then
186 c_b(i,k) = 0.0 ; d_b(i,k) = 1.0
188 d_b(i,k) = b_b * (h_b + d_b(i,k+1)*mix(i,k+1))
189 c_b(i,k) = 1.0 ;
if (d_b(i,k) > epsilon(b_b)) c_b(i,k) = b_b * mix(i,k)
193 te_b(i,k) = b_b * (h_b*tv%T(i,j,k) + mix(i,k+1)*te_b(i,k+1))
194 se_b(i,k) = b_b * (h_b*tv%S(i,j,k) + mix(i,k+1)*se_b(i,k+1))
196 te_b(i,k) = b_b * (h_b*tv%T(i,j,k))
197 se_b(i,k) = b_b * (h_b*tv%S(i,j,k))
199 endif ;
enddo ;
enddo
203 do i=is,ie ;
if (do_i(i))
then
204 if (change_ct(i) == 0)
then
205 last_down(i) = .false. ; do_i(i) = .false.
207 changed_col = changed_col + 1 ; change_ct(i) = 0
210 if (changed_col == 0)
exit
215 do i=is,ie ; do_i(i) = ((g%mask2dT(i,j) > 0.0) .and. last_down(i)) ;
enddo
216 do i=is,ie ;
if (do_i(i))
then
217 h_a = h(i,j,nz) + h_neglect
218 b_a = 1.0 / (h_a + d_a(i,nz)*mix(i,nz))
219 t_adj(i,j,nz) = b_a * (h_a*tv%T(i,j,nz) + mix(i,nz)*te_a(i,nz-1))
220 s_adj(i,j,nz) = b_a * (h_a*tv%S(i,j,nz) + mix(i,nz)*se_a(i,nz-1))
222 do k=nz-1,1,-1 ;
do i=is,ie ;
if (do_i(i))
then
223 t_adj(i,j,k) = te_a(i,k) + c_a(i,k+1)*t_adj(i,j,k+1)
224 s_adj(i,j,k) = se_a(i,k) + c_a(i,k+1)*s_adj(i,j,k+1)
225 endif ;
enddo ;
enddo
227 do i=is,ie ;
if (do_i(i))
then
234 do_i(i) = ((g%mask2dT(i,j) > 0.0) .and. .not.last_down(i))
236 h_b = h(i,j,1) + h_neglect
237 b_b = 1.0 / (h_b + d_b(i,2)*mix(i,2))
238 t_adj(i,j,1) = b_b * (h_b*tv%T(i,j,1) + mix(i,2)*te_b(i,2))
239 s_adj(i,j,1) = b_b * (h_b*tv%S(i,j,1) + mix(i,2)*se_b(i,2))
240 elseif (g%mask2dT(i,j) <= 0.0)
then
241 t_adj(i,j,1) = tv%T(i,j,1) ; s_adj(i,j,1) = tv%S(i,j,1)
244 do k=2,nz ;
do i=is,ie
246 t_adj(i,j,k) = te_b(i,k) + c_b(i,k)*t_adj(i,j,k-1)
247 s_adj(i,j,k) = se_b(i,k) + c_b(i,k)*s_adj(i,j,k-1)
248 elseif (g%mask2dT(i,j) <= 0.0)
then
249 t_adj(i,j,k) = tv%T(i,j,k) ; s_adj(i,j,k) = tv%S(i,j,k)
253 do i=is,ie ;
if (do_i(i))
then
282 function is_unstable(dRho_dT, dRho_dS, h_a, h_b, mix_A, mix_B, T_a, T_b, S_a, S_b, &
283 Te_aa, Te_bb, Se_aa, Se_bb, d_A, d_B)
284 real,
intent(in) :: drho_dt
285 real,
intent(in) :: drho_ds
286 real,
intent(in) :: h_a
287 real,
intent(in) :: h_b
288 real,
intent(in) :: mix_a
289 real,
intent(in) :: mix_b
290 real,
intent(in) :: t_a
291 real,
intent(in) :: t_b
292 real,
intent(in) :: s_a
293 real,
intent(in) :: s_b
294 real,
intent(in) :: te_aa
295 real,
intent(in) :: te_bb
296 real,
intent(in) :: se_aa
297 real,
intent(in) :: se_bb
298 real,
intent(in) :: d_a
299 real,
intent(in) :: d_b
307 is_unstable = (drho_dt * ((h_a * h_b * (t_b - t_a) + &
308 mix_a*mix_b * (d_a*te_bb - d_b*te_aa)) + &
309 (h_a*mix_b * (te_bb - d_b*t_a) + &
310 h_b*mix_a * (d_a*t_b - te_aa)) ) + &
311 drho_ds * ((h_a * h_b * (s_b - s_a) + &
312 mix_a*mix_b * (d_a*se_bb - d_b*se_aa)) + &
313 (h_a*mix_b * (se_bb - d_b*s_a) + &
314 h_b*mix_a * (d_a*s_b - se_aa)) ) < 0.0)
323 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
327 real,
intent(in) :: Kddt
328 real,
dimension(SZI_(G),SZK_(G)+1), &
331 real,
dimension(SZI_(G),SZK_(G)+1), &
334 integer,
intent(in) :: j
335 real,
dimension(:,:),
pointer :: p_surf
336 integer,
optional,
intent(in) :: halo
339 real :: mix(SZI_(G),SZK_(G)+1)
341 real :: b1(SZI_(G)), d1(SZI_(G))
342 real :: c1(SZI_(G),SZK_(G))
343 real :: T_f(SZI_(G),SZK_(G))
344 real :: S_f(SZI_(G),SZK_(G))
345 real :: pres(SZI_(G))
346 real :: T_EOS(SZI_(G))
347 real :: S_EOS(SZI_(G))
349 real :: h_neglect, h0
352 integer :: i, k, is, ie, nz
354 if (
present(halo))
then
355 is = g%isc-halo ; ie = g%iec+halo
357 is = g%isc ; ie = g%iec
361 h_neglect = gv%H_subroundoff
364 if (kddt <= 0.0)
then
365 do k=1,nz ;
do i=is,ie
366 t_f(i,k) = tv%T(i,j,k) ; s_f(i,k) = tv%S(i,j,k)
369 h0 = 1.0e-16*sqrt(kddt) + h_neglect
371 mix(i,2) = kap_dt_x2 / ((h(i,j,1)+h(i,j,2)) + h0)
373 h_tr = h(i,j,1) + h_neglect
374 b1(i) = 1.0 / (h_tr + mix(i,2))
375 d1(i) = b1(i) * h(i,j,1)
376 t_f(i,1) = (b1(i)*h_tr)*tv%T(i,j,1)
377 s_f(i,1) = (b1(i)*h_tr)*tv%S(i,j,1)
379 do k=2,nz-1 ;
do i=is,ie
380 mix(i,k+1) = kap_dt_x2 / ((h(i,j,k)+h(i,j,k+1)) + h0)
382 c1(i,k) = mix(i,k) * b1(i)
383 h_tr = h(i,j,k) + h_neglect
384 b1(i) = 1.0 / ((h_tr + d1(i)*mix(i,k)) + mix(i,k+1))
385 d1(i) = b1(i) * (h_tr + d1(i)*mix(i,k))
386 t_f(i,k) = b1(i) * (h_tr*tv%T(i,j,k) + mix(i,k)*t_f(i,k-1))
387 s_f(i,k) = b1(i) * (h_tr*tv%S(i,j,k) + mix(i,k)*s_f(i,k-1))
390 c1(i,nz) = mix(i,nz) * b1(i)
391 h_tr = h(i,j,nz) + h_neglect
392 b1(i) = 1.0 / (h_tr + d1(i)*mix(i,nz))
393 t_f(i,nz) = b1(i) * (h_tr*tv%T(i,j,nz) + mix(i,nz)*t_f(i,nz-1))
394 s_f(i,nz) = b1(i) * (h_tr*tv%S(i,j,nz) + mix(i,nz)*s_f(i,nz-1))
396 do k=nz-1,1,-1 ;
do i=is,ie
397 t_f(i,k) = t_f(i,k) + c1(i,k+1)*t_f(i,k+1)
398 s_f(i,k) = s_f(i,k) + c1(i,k+1)*s_f(i,k+1)
402 if (
associated(p_surf))
then
403 do i=is,ie ; pres(i) = p_surf(i,j) ;
enddo
405 do i=is,ie ; pres(i) = 0.0 ;
enddo
408 is-g%isd+1, ie-is+1, tv%eqn_of_state)
409 do i=is,ie ; pres(i) = pres(i) + h(i,j,1)*gv%H_to_Pa ;
enddo
412 t_eos(i) = 0.5*(t_f(i,k-1) + t_f(i,k))
413 s_eos(i) = 0.5*(s_f(i,k-1) + s_f(i,k))
416 is-g%isd+1, ie-is+1, tv%eqn_of_state)
417 do i=is,ie ; pres(i) = pres(i) + h(i,j,k)*gv%H_to_Pa ;
enddo
420 is-g%isd+1, ie-is+1, tv%eqn_of_state)