This subroutine ensures that there is a degree of horizontal smoothness in the depths of the near-surface interfaces.
116 type(ocean_grid_type),
intent(inout) :: G
117 type(verticalGrid_type),
intent(in) :: GV
118 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
120 type(thermo_var_ptrs),
intent(inout) :: tv
123 real,
intent(in) :: dt
124 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
128 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
132 type(regularize_layers_CS),
pointer :: CS
135 real,
dimension(SZIB_(G),SZJ_(G)) :: &
137 real,
dimension(SZI_(G),SZJB_(G)) :: &
139 real,
dimension(SZI_(G),SZJ_(G)) :: &
141 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
145 real,
dimension(SZIB_(G),SZJ_(G)) :: &
146 def_rat_u_1b, def_rat_u_2, def_rat_u_2b, def_rat_u_3, def_rat_u_3b
147 real,
dimension(SZI_(G),SZJB_(G)) :: &
148 def_rat_v_1b, def_rat_v_2, def_rat_v_2b, def_rat_v_3, def_rat_v_3b
149 real,
dimension(SZI_(G),SZJB_(G)) :: &
150 def_rat_h2, def_rat_h3
151 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
155 real,
dimension(SZI_(G),SZK_(G)+1) :: &
157 real,
dimension(SZI_(G),SZK_(G)) :: &
171 real,
dimension(SZI_(G)) :: &
176 h_add_tgt, h_add_tot, &
177 h_tot1, th_tot1, sh_tot1, &
178 h_tot3, th_tot3, sh_tot3, &
179 h_tot2, th_tot2, sh_tot2
180 real,
dimension(SZK_(G)) :: &
185 real :: e_e, e_w, e_n, e_s
190 real,
dimension(SZK_(G)+1) :: &
191 int_flux, int_Tflux, int_Sflux, int_Rflux
198 real :: int_top, int_bot
203 logical :: cols_left, ent_any, more_ent_i(SZI_(G)), ent_i(SZI_(G))
204 logical :: det_any, det_i(SZI_(G))
205 logical :: do_j(SZJ_(G)), do_i(SZI_(G)), find_i(SZI_(G))
206 logical :: debug = .false.
207 logical :: fatal_error
208 character(len=256) :: mesg
209 integer :: i, j, k, is, ie, js, je, nz, nkmb, nkml, k1, k2, k3, ks, nz_filt
211 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
213 if (.not.
associated(cs))
call mom_error(fatal,
"MOM_regularize_layers: "//&
214 "Module must be initialized before it is used.")
216 if (gv%nkml<1)
return
217 nkmb = gv%nk_rho_varies ; nkml = gv%nkml
218 if (.not.
associated(tv%eqn_of_state))
call mom_error(fatal, &
219 "MOM_regularize_layers: This module now requires the use of temperature and "//&
220 "an equation of state.")
222 h_neglect = gv%H_subroundoff
223 debug = (debug .or. cs%debug)
226 if (cs%id_def_rat_2 > 0)
then
227 is = g%isc-1 ; ie = g%iec+1 ; js = g%jsc-1 ; je = g%jec+1
231 i_dtol = 1.0 / max(cs%h_def_tol2 - cs%h_def_tol1, 1e-40)
232 i_dtol34 = 1.0 / max(cs%h_def_tol4 - cs%h_def_tol3, 1e-40)
234 p_ref_cv(:) = tv%P_Ref
236 do j=js-1,je+1 ;
do i=is-1,ie+1
239 do k=1,nz ;
do j=js-1,je+1 ;
do i=is-1,ie+1
240 e(i,j,k+1) = e(i,j,k) - h(i,j,k)
241 enddo ;
enddo ;
enddo
244 call find_deficit_ratios(e, def_rat_u, def_rat_v, g, gv, cs, def_rat_u_1b, def_rat_v_1b, 1, h)
246 call find_deficit_ratios(e, def_rat_u, def_rat_v, g, gv, cs, h=h)
249 do j=js,je ; do_j(j) = .false. ;
enddo
250 do j=js,je ;
do i=is,ie
251 def_rat_h(i,j) = max(def_rat_u(i-1,j), def_rat_u(i,j), &
252 def_rat_v(i,j-1), def_rat_v(i,j))
253 if (def_rat_h(i,j) > cs%h_def_tol1) do_j(j) = .true.
257 if ((cs%id_def_rat_3 > 0) .or. (cs%id_e3 > 0) .or. &
258 (cs%id_def_rat_u_3 > 0) .or. (cs%id_def_rat_u_3b > 0) .or. &
259 (cs%id_def_rat_v_3 > 0) .or. (cs%id_def_rat_v_3b > 0) )
then
260 do j=js-1,je+1 ;
do i=is-1,ie+1
263 do k=2,nz+1 ;
do j=js,je ;
do i=is,ie
264 if (g%mask2dCu(i,j) <= 0.0)
then ; e_e = e(i,j,k) ;
else
265 e_e = max(e(i+1,j,k) + min(e(i,j,k) - e(i+1,j,nz+1), 0.0), e(i,j,nz+1))
267 if (g%mask2dCu(i-1,j) <= 0.0)
then ; e_w = e(i,j,k) ;
else
268 e_w = max(e(i-1,j,k) + min(e(i,j,k) - e(i-1,j,nz+1), 0.0), e(i,j,nz+1))
270 if (g%mask2dCv(i,j) <= 0.0)
then ; e_n = e(i,j,k) ;
else
271 e_n = max(e(i,j+1,k) + min(e(i,j,k) - e(i,j+1,nz+1), 0.0), e(i,j,nz+1))
273 if (g%mask2dCv(i,j-1) <= 0.0)
then ; e_s = e(i,j,k) ;
else
274 e_s = max(e(i,j-1,k) + min(e(i,j,k) - e(i,j-1,nz+1), 0.0), e(i,j,nz+1))
278 ef(i,j,k) = (1.0 - 0.5*wt) * e(i,j,k) + &
279 wt * 0.125 * ((e_e + e_w) + (e_n + e_s))
280 enddo ;
enddo ;
enddo
281 call find_deficit_ratios(ef, def_rat_u_3, def_rat_v_3, g, gv, cs, def_rat_u_3b, def_rat_v_3b)
284 do j=js,je ;
do i=is,ie
285 def_rat_h3(i,j) = max(def_rat_u_3(i-1,j), def_rat_u_3(i,j), &
286 def_rat_v_3(i,j-1), def_rat_v_3(i,j))
289 if (cs%id_e3 > 0)
call post_data(cs%id_e3, ef, cs%diag)
290 if (cs%id_def_rat_3 > 0)
call post_data(cs%id_def_rat_3, def_rat_h3, cs%diag)
291 if (cs%id_def_rat_u_3 > 0)
call post_data(cs%id_def_rat_u_3, def_rat_u_3, cs%diag)
292 if (cs%id_def_rat_u_3b > 0)
call post_data(cs%id_def_rat_u_3b, def_rat_u_3b, cs%diag)
293 if (cs%id_def_rat_v_3 > 0)
call post_data(cs%id_def_rat_v_3, def_rat_v_3, cs%diag)
294 if (cs%id_def_rat_v_3b > 0)
call post_data(cs%id_def_rat_v_3b, def_rat_v_3b, cs%diag)
314 do j=js,je ;
if (do_j(j))
then
321 do k=1,nz ;
do i=is,ie ; d_ea(i,k) = 0.0 ; d_eb(i,k) = 0.0 ;
enddo ;
enddo
325 do_i(i) = def_rat_h(i,j) > cs%h_def_tol1
326 if (def_rat_h(i,j) > max_def_rat) max_def_rat = def_rat_h(i,j)
328 nz_filt = nkmb+1 ;
if (max_def_rat > cs%h_def_tol3) nz_filt = nz+1
333 do k=1,nz_filt ;
do i=is,ie ;
if (do_i(i))
then
334 if (g%mask2dCu(i,j) <= 0.0)
then ; e_e = e(i,j,k) ;
else
335 e_e = max(e(i+1,j,k) + min(e(i,j,k) - e(i+1,j,nz+1), 0.0), &
336 e(i,j,nz+1) + (nz+1-k)*gv%Angstrom_H)
339 if (g%mask2dCu(i-1,j) <= 0.0)
then ; e_w = e(i,j,k) ;
else
340 e_w = max(e(i-1,j,k) + min(e(i,j,k) - e(i-1,j,nz+1), 0.0), &
341 e(i,j,nz+1) + (nz+1-k)*gv%Angstrom_H)
343 if (g%mask2dCv(i,j) <= 0.0)
then ; e_n = e(i,j,k) ;
else
344 e_n = max(e(i,j+1,k) + min(e(i,j,k) - e(i,j+1,nz+1), 0.0), &
345 e(i,j,nz+1) + (nz+1-k)*gv%Angstrom_H)
347 if (g%mask2dCv(i,j-1) <= 0.0)
then ; e_s = e(i,j,k) ;
else
348 e_s = max(e(i,j-1,k) + min(e(i,j,k) - e(i,j-1,nz+1), 0.0), &
349 e(i,j,nz+1) + (nz+1-k)*gv%Angstrom_H)
352 wt = max(0.0, min(1.0, i_dtol*(def_rat_h(i,j)-cs%h_def_tol1)))
354 e_filt(i,k) = (1.0 - 0.5*wt) * e(i,j,k) + &
355 wt * 0.125 * ((e_e + e_w) + (e_n + e_s))
357 endif ;
enddo ;
enddo
358 do k=1,nz ;
do i=is,ie
360 t_2d(i,k) = tv%T(i,j,k) ; s_2d(i,k) = tv%S(i,j,k)
364 do k=1,nz ;
do i=is,ie ;
if (do_i(i))
then
365 h_2d_init(i,k) = h(i,j,k)
366 t_2d_init(i,k) = tv%T(i,j,k) ; s_2d_init(i,k) = tv%S(i,j,k)
367 endif ;
enddo ;
enddo
373 more_ent_i(i) = .false. ; ent_i(i) = .false.
374 h_add_tgt(i) = 0.0 ; h_add_tot(i) = 0.0
375 if (do_i(i) .and. (e_2d(i,nkmb+1) > e_filt(i,nkmb+1)))
then
376 more_ent_i(i) = .true. ; ent_i(i) = .true. ; ent_any = .true.
377 h_add_tgt(i) = e_2d(i,nkmb+1) - e_filt(i,nkmb+1)
384 do i=is,ie ;
if (more_ent_i(i))
then
385 if (h_2d(i,k) - gv%Angstrom_H > h_neglect)
then
386 if (e_2d(i,nkmb+1)-e_filt(i,nkmb+1) > h_2d(i,k) - gv%Angstrom_H)
then
387 h_add = h_2d(i,k) - gv%Angstrom_H
388 h_2d(i,k) = gv%Angstrom_H
390 h_add = e_2d(i,nkmb+1)-e_filt(i,nkmb+1)
391 h_2d(i,k) = h_2d(i,k) - h_add
393 d_eb(i,k-1) = d_eb(i,k-1) + h_add
394 h_add_tot(i) = h_add_tot(i) + h_add
395 e_2d(i,nkmb+1) = e_2d(i,nkmb+1) - h_add
396 h_prev = h_2d(i,nkmb)
397 h_2d(i,nkmb) = h_2d(i,nkmb) + h_add
399 t_2d(i,nkmb) = (h_prev*t_2d(i,nkmb) + h_add*t_2d(i,k)) / h_2d(i,nkmb)
400 s_2d(i,nkmb) = (h_prev*s_2d(i,nkmb) + h_add*s_2d(i,k)) / h_2d(i,nkmb)
402 if ((e_2d(i,nkmb+1) <= e_filt(i,nkmb+1)) .or. &
403 (h_add_tot(i) > 0.6*h_add_tgt(i)))
then
404 more_ent_i(i) = .false.
412 if (.not.cols_left)
exit
416 do k=ks,nkmb,-1 ;
do i=is,ie ;
if (ent_i(i))
then
417 d_eb(i,k) = d_eb(i,k) + d_eb(i,k+1)
418 endif ;
enddo ;
enddo
430 if ((max_def_rat > cs%h_def_tol3) .and. (cs%reg_sfc_detrain))
then
432 det_i(i) = .false. ; rcv_tol(i) = 0.0
433 if (do_i(i) .and. (e_2d(i,nkmb+1) < e_filt(i,nkmb+1)) .and. &
434 (def_rat_h(i,j) > cs%h_def_tol3))
then
435 det_i(i) = .true. ; det_any = .true.
436 rcv_tol(i) = min((def_rat_h(i,j) - cs%h_def_tol3), 1.0)
441 call cpu_clock_begin(id_clock_eos)
443 call calculate_density(t_2d(:,k),s_2d(:,k),p_ref_cv,rcv(:,k), &
444 is,ie-is+1,tv%eqn_of_state)
446 call cpu_clock_end(id_clock_eos)
448 do i=is,ie ;
if (det_i(i))
then
455 rcv_min_det = gv%Rlay(k2) + 0.6*rcv_tol(i)*(gv%Rlay(k2-1)-gv%Rlay(k2))
457 rcv_max_det = gv%Rlay(k2) + 0.6*rcv_tol(i)*(gv%Rlay(k2+1)-gv%Rlay(k2))
459 rcv_max_det = gv%Rlay(nz) + 0.6*rcv_tol(i)*(gv%Rlay(nz)-gv%Rlay(nz-1))
461 if (rcv(i,k1) > rcv_max_det) &
464 h_deficit = (e_filt(i,k2)-e_filt(i,k2+1)) - h_2d(i,k2)
465 if ((e_filt(i,k2) > e_2d(i,k1+1)) .and. (h_deficit > 0.0) .and. &
466 (rcv(i,k1) < rcv_max_det) .and. (rcv(i,k1) > rcv_min_det))
then
468 h_add = min(e_filt(i,k2) - e_2d(i,k2), h_deficit )
469 if (h_add < h_2d(i,k1))
then
471 if (h_add > 0.0)
then
473 h_2d(i,k2) = h_2d(i,k2) + h_add
474 e_2d(i,k2) = e_2d(i,k2+1) + h_2d(i,k2)
475 d_ea(i,k2) = d_ea(i,k2) + h_add
477 t_2d(i,k2) = (h_prev*t_2d(i,k2) + h_add*t_2d(i,k1)) / h_2d(i,k2)
478 s_2d(i,k2) = (h_prev*s_2d(i,k2) + h_add*s_2d(i,k1)) / h_2d(i,k2)
479 h_det_tot = h_det_tot + h_add
481 h_2d(i,k1) = h_2d(i,k1) - h_add
482 do k3=k1,nkmb ; e_2d(i,k3+1) = e_2d(i,k3) - h_2d(i,k3) ;
enddo
483 do k3=k1+1,nkmb ; d_ea(i,k3) = d_ea(i,k3) + h_add ;
enddo
486 call mom_error(fatal,
"h_add is negative. Some logic is wrong.")
492 if (k2>nkmb+1) e_2d(i,k2) = e_2d(i,k2) + h_det_tot
496 h_2d(i,k2) = h_2d(i,k2) + h_add
497 e_2d(i,k2) = e_2d(i,k2+1) + h_2d(i,k2)
498 d_ea(i,k2) = d_ea(i,k2) + h_add
499 t_2d(i,k2) = (h_prev*t_2d(i,k2) + h_add*t_2d(i,k1)) / h_2d(i,k2)
500 s_2d(i,k2) = (h_prev*s_2d(i,k2) + h_add*s_2d(i,k1)) / h_2d(i,k2)
501 h_det_tot = h_det_tot + h_add
504 do k3=k1,nkmb ; e_2d(i,k3+1) = e_2d(i,k3) - h_2d(i,k3) ;
enddo
505 do k3=k1+1,nkmb ; d_ea(i,k3) = d_ea(i,k3) + h_add ;
enddo
514 if (k2>nkmb+1) e_2d(i,k2) = e_2d(i,k2) + h_det_tot
520 do k=nz-1,nkmb+1,-1 ;
do i=is,ie ;
if (det_i(i))
then
521 d_ea(i,k) = d_ea(i,k) + d_ea(i,k+1)
522 endif ;
enddo ;
enddo
525 do i=is,ie ; h_tot3(i) = 0.0 ; th_tot3(i) = 0.0 ; sh_tot3(i) = 0.0 ;
enddo
526 do k=1,nz ;
do i=is,ie ;
if (do_i(i))
then
527 h_tot3(i) = h_tot3(i) + h_2d(i,k)
528 th_tot3(i) = th_tot3(i) + h_2d(i,k) * t_2d(i,k)
529 sh_tot3(i) = sh_tot3(i) + h_2d(i,k) * s_2d(i,k)
530 endif ;
enddo ;
enddo
533 do i=is,ie ;
if (do_i(i))
then
536 scale = e_2d(i,nkmb+1) / e_filt(i,nkmb+1)
537 do k=2,nkmb+1 ; e_filt(i,k) = e_filt(i,k) * scale ;
enddo
541 if (e_filt(i,2) < e_2d(i,nkml))
then
542 scale = (e_2d(i,nkml) - e_filt(i,nkmb+1)) / &
543 ((e_filt(i,2) - e_filt(i,nkmb+1)) + h_neglect)
545 e_filt(i,k) = e_filt(i,nkmb+1) + scale * (e_filt(i,k) - e_filt(i,nkmb+1))
547 e_filt(i,2) = e_2d(i,nkml)
554 int_flux(k) = 0.0 ; int_rflux(k) = 0.0
555 int_tflux(k) = 0.0 ; int_sflux(k) = 0.0
558 int_bot = max(e_2d(i,k1+1),e_filt(i,k2+1))
559 h_add = int_top - int_bot
563 d_ea(i,k3) = d_ea(i,k3) + h_add
564 int_flux(k3) = int_flux(k3) + h_add
565 int_tflux(k3) = int_tflux(k3) + h_add*t_2d(i,k1)
566 int_sflux(k3) = int_sflux(k3) + h_add*s_2d(i,k1)
568 elseif (k1 > k2)
then
570 d_eb(i,k3) = d_eb(i,k3) + h_add
571 int_flux(k3+1) = int_flux(k3+1) - h_add
572 int_tflux(k3+1) = int_tflux(k3+1) - h_add*t_2d(i,k1)
573 int_sflux(k3+1) = int_sflux(k3+1) - h_add*s_2d(i,k1)
577 if (int_bot <= e_filt(i,k2+1))
then
580 elseif (int_bot <= e_2d(i,k1+1))
then
584 call mom_error(fatal, &
585 "Regularize_surface: Could not increment target or source.")
587 if ((k1 > nkmb) .or. (k2 > nkmb))
exit
591 call mom_error(fatal,
"Regularize_surface: Did not assign fluid to layer nkmb.")
595 do k=1,nkmb ; h_prev_1d(k) = h_2d(i,k) ;
enddo
596 h_2d(i,1) = h_2d(i,1) - int_flux(2)
598 h_2d(i,k) = h_2d(i,k) + (int_flux(k) - int_flux(k+1))
602 h_2d(i,nkmb) = h_2d(i,nkmb) + int_flux(nkmb)
604 t_2d(i,1) = (t_2d(i,1)*h_prev_1d(1) - int_tflux(2)) / h_2d(i,1)
605 s_2d(i,1) = (s_2d(i,1)*h_prev_1d(1) - int_sflux(2)) / h_2d(i,1)
607 t_2d(i,k) = (t_2d(i,k)*h_prev_1d(k) + (int_tflux(k) - int_tflux(k+1))) / h_2d(i,k)
608 s_2d(i,k) = (s_2d(i,k)*h_prev_1d(k) + (int_sflux(k) - int_sflux(k+1))) / h_2d(i,k)
610 t_2d(i,nkmb) = (t_2d(i,nkmb)*h_prev_1d(nkmb) + int_tflux(nkmb) ) / h_2d(i,nkmb)
611 s_2d(i,nkmb) = (s_2d(i,nkmb)*h_prev_1d(nkmb) + int_sflux(nkmb) ) / h_2d(i,nkmb)
616 do k=1,nz ;
do i=is,ie ;
if (do_i(i))
then
618 tv%T(i,j,k) = t_2d(i,k) ; tv%S(i,j,k) = s_2d(i,k)
619 ea(i,j,k) = ea(i,j,k) + d_ea(i,k)
620 eb(i,j,k) = eb(i,j,k) + d_eb(i,k)
621 endif ;
enddo ;
enddo
624 do i=is,ie ; h_tot1(i) = 0.0 ; th_tot1(i) = 0.0 ; sh_tot1(i) = 0.0 ;
enddo
625 do i=is,ie ; h_tot2(i) = 0.0 ; th_tot2(i) = 0.0 ; sh_tot2(i) = 0.0 ;
enddo
627 do k=1,nz ;
do i=is,ie ;
if (do_i(i))
then
628 h_tot1(i) = h_tot1(i) + h_2d_init(i,k)
629 h_tot2(i) = h_tot2(i) + h(i,j,k)
631 th_tot1(i) = th_tot1(i) + h_2d_init(i,k) * t_2d_init(i,k)
632 th_tot2(i) = th_tot2(i) + h(i,j,k) * tv%T(i,j,k)
633 sh_tot1(i) = sh_tot1(i) + h_2d_init(i,k) * s_2d_init(i,k)
634 sh_tot2(i) = sh_tot2(i) + h(i,j,k) * tv%S(i,j,k)
635 if (h(i,j,k) < 0.0) &
636 call mom_error(fatal,
"regularize_surface: Negative thicknesses.")
637 if (k==1)
then ; h_predicted = h_2d_init(i,k) + (d_eb(i,k) - d_ea(i,k+1))
638 elseif (k==nz)
then ; h_predicted = h_2d_init(i,k) + (d_ea(i,k) - d_eb(i,k-1))
640 h_predicted = h_2d_init(i,k) + ((d_ea(i,k) - d_eb(i,k-1)) + &
641 (d_eb(i,k) - d_ea(i,k+1)))
643 if (abs(h(i,j,k) - h_predicted) > max(1e-9*abs(h_predicted),gv%Angstrom_H)) &
644 call mom_error(fatal,
"regularize_surface: d_ea mismatch.")
645 endif ;
enddo ;
enddo
646 do i=is,ie ;
if (do_i(i))
then
647 fatal_error = .false.
648 if (abs(h_tot1(i) - h_tot2(i)) > 1e-12*h_tot1(i))
then
649 write(mesg,
'(ES11.4," became ",ES11.4," diff ",ES11.4)') &
650 h_tot1(i), h_tot2(i), (h_tot1(i) - h_tot2(i))
651 call mom_error(warning,
"regularize_surface: Mass non-conservation."//&
655 if (abs(th_tot1(i) - th_tot2(i)) > 1e-12*(th_tot1(i)+10.0*h_tot1(i)))
then
656 write(mesg,
'(ES11.4," became ",ES11.4," diff ",ES11.4," int diff ",ES11.4)') &
657 th_tot1(i), th_tot2(i), (th_tot1(i) - th_tot2(i)), (th_tot1(i) - th_tot3(i))
658 call mom_error(warning,
"regularize_surface: Heat non-conservation."//&
662 if (abs(sh_tot1(i) - sh_tot2(i)) > 1e-12*(sh_tot1(i)+10.0*h_tot1(i)))
then
663 write(mesg,
'(ES11.4," became ",ES11.4," diff ",ES11.4," int diff ",ES11.4)') &
664 sh_tot1(i), sh_tot2(i), (sh_tot1(i) - sh_tot2(i)), (sh_tot1(i) - sh_tot3(i))
665 call mom_error(warning,
"regularize_surface: Salinity non-conservation."//&
669 if (fatal_error)
then
670 write(mesg,
'("Error at lat/lon ",2(ES11.4))') g%geoLatT(i,j), g%geoLonT(i,j)
671 call mom_error(fatal,
"regularize_surface: Terminating with fatal error. "//&
679 if (cs%id_def_rat > 0)
call post_data(cs%id_def_rat, def_rat_h, cs%diag)
682 if (cs%id_e1 > 0)
call post_data(cs%id_e1, e, cs%diag)
683 if (cs%id_def_rat_u > 0)
call post_data(cs%id_def_rat_u, def_rat_u, cs%diag)
684 if (cs%id_def_rat_u_1b > 0)
call post_data(cs%id_def_rat_u_1b, def_rat_u_1b, cs%diag)
685 if (cs%id_def_rat_v > 0)
call post_data(cs%id_def_rat_v, def_rat_v, cs%diag)
686 if (cs%id_def_rat_v_1b > 0)
call post_data(cs%id_def_rat_v_1b, def_rat_v_1b, cs%diag)
688 if ((cs%id_def_rat_2 > 0) .or. &
689 (cs%id_def_rat_u_2 > 0) .or. (cs%id_def_rat_u_2b > 0) .or. &
690 (cs%id_def_rat_v_2 > 0) .or. (cs%id_def_rat_v_2b > 0) )
then
691 do j=js-1,je+1 ;
do i=is-1,ie+1
694 do k=1,nz ;
do j=js-1,je+1 ;
do i=is-1,ie+1
695 e(i,j,k+1) = e(i,j,k) - h(i,j,k)
696 enddo ;
enddo ;
enddo
698 call find_deficit_ratios(e, def_rat_u_2, def_rat_v_2, g, gv, cs, def_rat_u_2b, def_rat_v_2b, h=h)
701 do j=js,je ;
do i=is,ie
702 def_rat_h2(i,j) = max(def_rat_u_2(i-1,j), def_rat_u_2(i,j), &
703 def_rat_v_2(i,j-1), def_rat_v_2(i,j))
706 if (cs%id_def_rat_2 > 0)
call post_data(cs%id_def_rat_2, def_rat_h2, cs%diag)
707 if (cs%id_e2 > 0)
call post_data(cs%id_e2, e, cs%diag)
708 if (cs%id_def_rat_u_2 > 0)
call post_data(cs%id_def_rat_u_2, def_rat_u_2, cs%diag)
709 if (cs%id_def_rat_u_2b > 0)
call post_data(cs%id_def_rat_u_2b, def_rat_u_2b, cs%diag)
710 if (cs%id_def_rat_v_2 > 0)
call post_data(cs%id_def_rat_v_2, def_rat_v_2, cs%diag)
711 if (cs%id_def_rat_v_2b > 0)
call post_data(cs%id_def_rat_v_2b, def_rat_v_2b, cs%diag)