Calculate new temperatures and salinities that have been subject to full convective mixing.
22 type(ocean_grid_type),
intent(in) :: G
23 type(verticalGrid_type),
intent(in) :: GV
24 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
26 type(thermo_var_ptrs),
intent(in) :: tv
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
110 call smoothed_drdt_drds(h, tv, kddt_smooth, drho_dt, drho_ds, g, gv, j, p_surf, halo)
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