7 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine
25 use mom_time_manager,
only : time_type, real_to_time,
operator(+),
operator(-)
30 implicit none ;
private
32 #include <MOM_memory.h>
37 # define WHALOI_ MAX(BTHALO_-NIHALO_,0)
38 # define WHALOJ_ MAX(BTHALO_-NJHALO_,0)
39 # define NIMEMW_ 1-WHALOI_:NIMEM_+WHALOI_
40 # define NJMEMW_ 1-WHALOJ_:NJMEM_+WHALOJ_
41 # define NIMEMBW_ -WHALOI_:NIMEM_+WHALOI_
42 # define NJMEMBW_ -WHALOJ_:NJMEM_+WHALOJ_
43 # define SZIW_(G) NIMEMW_
44 # define SZJW_(G) NJMEMW_
45 # define SZIBW_(G) NIMEMBW_
46 # define SZJBW_(G) NJMEMBW_
52 # define SZIW_(G) G%isdw:G%iedw
53 # define SZJW_(G) G%jsdw:G%jedw
54 # define SZIBW_(G) G%isdw-1:G%iedw
55 # define SZJBW_(G) G%jsdw-1:G%jedw
58 public btcalc, bt_mass_source, btstep, barotropic_init, barotropic_end
59 public register_barotropic_restarts, set_dtbt, barotropic_get_tav
68 real,
dimension(:,:),
pointer :: cg_u => null()
69 real,
dimension(:,:),
pointer :: cg_v => null()
70 real,
dimension(:,:),
pointer :: h_u => null()
71 real,
dimension(:,:),
pointer :: h_v => null()
72 real,
dimension(:,:),
pointer :: uhbt => null()
74 real,
dimension(:,:),
pointer :: vhbt => null()
76 real,
dimension(:,:),
pointer :: ubt_outer => null()
78 real,
dimension(:,:),
pointer :: vbt_outer => null()
80 real,
dimension(:,:),
pointer :: eta_outer_u => null()
82 real,
dimension(:,:),
pointer :: eta_outer_v => null()
84 logical :: apply_u_obcs
85 logical :: apply_v_obcs
87 integer :: is_u_obc, ie_u_obc, js_u_obc, je_u_obc
88 integer :: is_v_obc, ie_v_obc, js_v_obc, je_v_obc
90 logical :: is_alloced = .false.
92 type(group_pass_type) :: pass_uv
93 type(group_pass_type) :: pass_uhvh
94 type(group_pass_type) :: pass_h
95 type(group_pass_type) :: pass_cg
96 type(group_pass_type) :: pass_eta_outer
101 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: frhatu
103 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: frhatv
105 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_) :: idatu
107 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_) :: lin_drag_u
110 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_) :: uhbt_ic
113 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_) :: ubt_ic
116 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_) :: ubtav
118 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_) :: idatv
120 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_) :: lin_drag_v
123 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_) :: vhbt_ic
126 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_) :: vbt_ic
129 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_) :: vbtav
131 real allocable_,
dimension(NIMEM_,NJMEM_) :: eta_cor
135 real allocable_,
dimension(NIMEM_,NJMEM_) :: eta_cor_bound
138 real allocable_,
dimension(NIMEMW_,NJMEMW_) :: &
139 ua_polarity, & !< Test vector components for checking grid polarity.
140 va_polarity, & !< Test vector components for checking grid polarity.
142 real allocable_,
dimension(NIMEMW_,NJMEMW_) :: iareat
145 real allocable_,
dimension(NIMEMBW_,NJMEMW_) :: &
146 d_u_cor, & !< A simply averaged depth at u points [Z ~> m].
147 dy_cu, & !< A copy of G%dy_Cu with wide halos [L ~> m].
149 real allocable_,
dimension(NIMEMW_,NJMEMBW_) :: &
150 d_v_cor, & !< A simply averaged depth at v points [Z ~> m].
151 dx_cv, & !< A copy of G%dx_Cv with wide halos [L ~> m].
153 real allocable_,
dimension(NIMEMBW_,NJMEMBW_) :: &
156 real,
dimension(:,:,:),
pointer :: frhatu1 => null()
157 real,
dimension(:,:,:),
pointer :: frhatv1 => null()
164 real :: dtbt_fraction
171 integer :: nstep_last = 0
179 logical :: bound_bt_corr
184 logical :: gradual_bt_ics
195 logical :: nonlinear_continuity
197 integer :: nonlin_cont_update_period
201 logical :: bt_project_velocity
209 logical :: dynamic_psurf
211 real :: dmin_dyn_psurf
213 real :: ice_strength_length
216 real :: const_dyn_psurf
221 integer :: hvel_scheme
225 logical :: strong_drag
227 logical :: linear_wave_drag
230 logical :: linearized_bt_pv
233 logical :: use_wide_halos
235 logical :: clip_velocity
241 real :: vel_underflow
248 real :: maxcfl_bt_cont
253 logical :: bt_cont_bounds
256 logical :: visc_rem_u_uh0
260 logical :: adjust_bt_cont
263 logical :: use_old_coriolis_bracket_bug
266 type(time_type),
pointer :: time => null()
272 logical :: module_is_initialized = .false.
279 type(group_pass_type) :: pass_q_dcor
280 type(group_pass_type) :: pass_gtot
281 type(group_pass_type) :: pass_tmp_uv
282 type(group_pass_type) :: pass_eta_bt_rem
283 type(group_pass_type) :: pass_force_hbt0_cor_ref
284 type(group_pass_type) :: pass_dat_uv
285 type(group_pass_type) :: pass_eta_ubt
286 type(group_pass_type) :: pass_etaav
287 type(group_pass_type) :: pass_ubt_cor
288 type(group_pass_type) :: pass_ubta_uhbta
289 type(group_pass_type) :: pass_e_anom
292 integer :: id_pfu_bt = -1, id_pfv_bt = -1, id_coru_bt = -1, id_corv_bt = -1
293 integer :: id_ubtforce = -1, id_vbtforce = -1, id_uaccel = -1, id_vaccel = -1
294 integer :: id_visc_rem_u = -1, id_visc_rem_v = -1, id_eta_cor = -1
295 integer :: id_ubt = -1, id_vbt = -1, id_eta_bt = -1, id_ubtav = -1, id_vbtav = -1
296 integer :: id_ubt_st = -1, id_vbt_st = -1, id_eta_st = -1
297 integer :: id_ubt_hifreq = -1, id_vbt_hifreq = -1, id_eta_hifreq = -1
298 integer :: id_uhbt_hifreq = -1, id_vhbt_hifreq = -1, id_eta_pred_hifreq = -1
299 integer :: id_gtotn = -1, id_gtots = -1, id_gtote = -1, id_gtotw = -1
300 integer :: id_uhbt = -1, id_frhatu = -1, id_vhbt = -1, id_frhatv = -1
301 integer :: id_frhatu1 = -1, id_frhatv1 = -1
303 integer :: id_btc_fa_u_ee = -1, id_btc_fa_u_e0 = -1, id_btc_fa_u_w0 = -1, id_btc_fa_u_ww = -1
304 integer :: id_btc_ubt_ee = -1, id_btc_ubt_ww = -1
305 integer :: id_btc_fa_v_nn = -1, id_btc_fa_v_n0 = -1, id_btc_fa_v_s0 = -1, id_btc_fa_v_ss = -1
306 integer :: id_btc_vbt_nn = -1, id_btc_vbt_ss = -1
307 integer :: id_uhbt0 = -1, id_vhbt0 = -1
354 integer :: isdw, iedw, jsdw, jedw
359 integer :: id_clock_sync=-1, id_clock_calc=-1
360 integer :: id_clock_calc_pre=-1, id_clock_calc_post=-1
361 integer :: id_clock_pass_step=-1, id_clock_pass_pre=-1, id_clock_pass_post=-1
365 integer,
parameter :: harmonic = 1
366 integer,
parameter :: arithmetic = 2
367 integer,
parameter :: hybrid = 3
368 integer,
parameter :: from_bt_cont = 4
369 integer,
parameter :: hybrid_bt_cont = 5
370 character*(20),
parameter :: hybrid_string =
"HYBRID"
371 character*(20),
parameter :: harmonic_string =
"HARMONIC"
372 character*(20),
parameter :: arithmetic_string =
"ARITHMETIC"
373 character*(20),
parameter :: bt_cont_string =
"FROM_BT_CONT"
384 subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, &
385 eta_PF_in, U_Cor, V_Cor, accel_layer_u, accel_layer_v, &
386 eta_out, uhbtav, vhbtav, G, GV, US, CS, &
387 visc_rem_u, visc_rem_v, etaav, OBC, &
388 BT_cont, eta_PF_start, &
389 taux_bot, tauy_bot, uh0, vh0, u_uh0, v_vh0)
393 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: u_in
394 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: v_in
395 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: eta_in
397 real,
intent(in) :: dt
398 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: bc_accel_u
399 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: bc_accel_v
402 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: pbce
405 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: eta_pf_in
411 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: u_cor
413 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: v_cor
414 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(out) :: accel_layer_u
416 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(out) :: accel_layer_v
418 real,
dimension(SZI_(G),SZJ_(G)),
intent(out) :: eta_out
420 real,
dimension(SZIB_(G),SZJ_(G)),
intent(out) :: uhbtav
423 real,
dimension(SZI_(G),SZJB_(G)),
intent(out) :: vhbtav
428 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: visc_rem_u
434 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: visc_rem_v
435 real,
dimension(SZI_(G),SZJ_(G)),
optional,
intent(out) :: etaav
441 real,
dimension(:,:),
optional,
pointer :: eta_pf_start
444 real,
dimension(:,:),
optional,
pointer :: taux_bot
446 real,
dimension(:,:),
optional,
pointer :: tauy_bot
448 real,
dimension(:,:,:),
optional,
pointer :: uh0
450 real,
dimension(:,:,:),
optional,
pointer :: u_uh0
451 real,
dimension(:,:,:),
optional,
pointer :: vh0
453 real,
dimension(:,:,:),
optional,
pointer :: v_vh0
456 real :: ubt_cor(szib_(g),szj_(g))
457 real :: vbt_cor(szi_(g),szjb_(g))
459 real :: wt_u(szib_(g),szj_(g),szk_(g))
460 real :: wt_v(szi_(g),szjb_(g),szk_(g))
463 real,
dimension(SZIB_(G),SZJ_(G)) :: &
466 real,
dimension(SZI_(G),SZJB_(G)) :: &
469 real,
dimension(SZI_(G),SZJ_(G)) :: &
475 real :: q(szibw_(cs),szjbw_(cs))
476 real,
dimension(SZIBW_(CS),SZJW_(CS)) :: &
509 real,
dimension(SZIW_(CS),SZJBW_(CS)) :: &
540 real,
target,
dimension(SZIW_(CS),SZJW_(CS)) :: &
544 real,
dimension(:,:),
pointer :: &
547 real,
dimension(SZIW_(CS),SZJW_(CS)) :: &
574 real,
dimension(SZIBW_(CS),SZJW_(CS)) :: &
575 ubt_prev, uhbt_prev, ubt_sum_prev, uhbt_sum_prev, ubt_wtd_prev
576 real,
dimension(SZIW_(CS),SZJBW_(CS)) :: &
577 vbt_prev, vhbt_prev, vbt_sum_prev, vhbt_sum_prev, vbt_wtd_prev
605 logical :: do_hifreq_output
606 logical :: use_bt_cont, do_ave, find_etaav, find_pf, find_cor
607 logical :: ice_is_rigid, nonblock_setup, interp_eta_pf
608 logical :: project_velocity, add_uh0
612 real :: ice_strength = 0.0
620 real :: u_max_cor, v_max_cor
623 real :: accel_underflow
625 real,
allocatable,
dimension(:) :: wt_vel, wt_eta, wt_accel, wt_trans, wt_accel2
626 real :: sum_wt_vel, sum_wt_eta, sum_wt_accel, sum_wt_trans
627 real :: i_sum_wt_vel, i_sum_wt_eta, i_sum_wt_accel, i_sum_wt_trans
629 real :: trans_wt1, trans_wt2
632 logical :: apply_obcs, apply_obc_flather, apply_obc_open
634 character(len=200) :: mesg
635 integer :: isv, iev, jsv, jev
637 integer :: isvf, ievf, jsvf, jevf, num_cycles
638 integer :: i, j, k, n
639 integer :: is, ie, js, je, nz, isq, ieq, jsq, jeq
640 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
641 integer :: ioff, joff
643 if (.not.
associated(cs))
call mom_error(fatal, &
644 "btstep: Module MOM_barotropic must be initialized before it is used.")
645 if (.not.cs%split)
return
646 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
647 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
648 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
649 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
650 ms%isdw = cs%isdw ; ms%iedw = cs%iedw ; ms%jsdw = cs%jsdw ; ms%jedw = cs%jedw
651 dt_in_t = us%s_to_T*dt
653 accel_underflow = us%L_T_to_m_s*cs%vel_underflow * us%s_to_T*idt
655 use_bt_cont = .false.
656 if (
present(bt_cont)) use_bt_cont = (
associated(bt_cont))
658 interp_eta_pf = .false.
659 if (
present(eta_pf_start)) interp_eta_pf = (
associated(eta_pf_start))
661 project_velocity = cs%BT_project_velocity
665 if ((.not.use_bt_cont) .and. cs%Nonlinear_continuity .and. &
666 (cs%Nonlin_cont_update_period > 0)) stencil = 2
668 do_ave = query_averaging_enabled(cs%diag)
669 find_etaav =
present(etaav)
670 find_pf = (do_ave .and. ((cs%id_PFu_bt > 0) .or. (cs%id_PFv_bt > 0)))
671 find_cor = (do_ave .and. ((cs%id_Coru_bt > 0) .or. (cs%id_Corv_bt > 0)))
674 if (
present(uh0)) add_uh0 =
associated(uh0)
675 if (add_uh0 .and. .not.(
present(vh0) .and.
present(u_uh0) .and. &
676 present(v_vh0)))
call mom_error(fatal, &
677 "btstep: vh0, u_uh0, and v_vh0 must be present if uh0 is used.")
678 if (add_uh0 .and. .not.(
associated(vh0) .and.
associated(u_uh0) .and. &
679 associated(v_vh0)))
call mom_error(fatal, &
680 "btstep: vh0, u_uh0, and v_vh0 must be associated if uh0 is used.")
683 nonblock_setup = g%nonblocking_updates
685 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
687 apply_obcs = .false. ; cs%BT_OBC%apply_u_OBCs = .false. ; cs%BT_OBC%apply_v_OBCs = .false.
688 apply_obc_open = .false.
689 apply_obc_flather = .false.
690 if (
present(obc))
then ;
if (
associated(obc))
then
691 cs%BT_OBC%apply_u_OBCs = obc%open_u_BCs_exist_globally .or. obc%specified_u_BCs_exist_globally
692 cs%BT_OBC%apply_v_OBCs = obc%open_v_BCs_exist_globally .or. obc%specified_v_BCs_exist_globally
693 apply_obc_flather = open_boundary_query(obc, apply_flather_obc=.true.)
694 apply_obc_open = open_boundary_query(obc, apply_open_obc=.true.)
695 apply_obcs = open_boundary_query(obc, apply_specified_obc=.true.) .or. &
696 apply_obc_flather .or. apply_obc_open
698 if (apply_obc_flather .and. .not.gv%Boussinesq)
call mom_error(fatal, &
699 "btstep: Flather open boundary conditions have not yet been "// &
700 "implemented for a non-Boussinesq model.")
704 if (cs%use_wide_halos) &
705 num_cycles = min((is-cs%isdw) / stencil, (js-cs%jsdw) / stencil)
706 isvf = is - (num_cycles-1)*stencil ; ievf = ie + (num_cycles-1)*stencil
707 jsvf = js - (num_cycles-1)*stencil ; jevf = je + (num_cycles-1)*stencil
709 nstep = ceiling(dt/cs%dtbt - 0.0001)
710 if (is_root_pe() .and. (nstep /= cs%nstep_last))
then
711 write(mesg,
'("btstep is using a dynamic barotropic timestep of ", ES12.6, &
712 & " seconds, max ", ES12.6, ".")') (dt/nstep), cs%dtbt_max
713 call mom_mesg(mesg, 3)
715 cs%nstep_last = nstep
718 instep = 1.0 / real(nstep)
719 dtbt = dt_in_t * instep
722 mass_to_z = us%m_to_L*us%T_to_s**2 * us%m_to_Z / gv%Rho0
725 if (project_velocity)
then
726 trans_wt1 = (1.0 + be_proj); trans_wt2 = -be_proj
728 trans_wt1 = bebt ; trans_wt2 = (1.0-bebt)
731 do_hifreq_output = .false.
732 if ((cs%id_ubt_hifreq > 0) .or. (cs%id_vbt_hifreq > 0) .or. &
733 (cs%id_eta_hifreq > 0) .or. (cs%id_eta_pred_hifreq > 0) .or. &
734 (cs%id_uhbt_hifreq > 0) .or. (cs%id_vhbt_hifreq > 0))
then
735 do_hifreq_output = query_averaging_enabled(cs%diag, time_int_in, time_end_in)
736 if (do_hifreq_output) &
737 time_bt_start = time_end_in - real_to_time(dt)
741 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
742 if (.not. cs%linearized_BT_PV)
then
747 if ((isq > is-1) .or. (jsq > js-1)) &
750 to_all+scalar_pair, agrid)
752 to_all+scalar_pair, agrid)
754 if (cs%dynamic_psurf) &
756 if (interp_eta_pf)
then
767 cs%BT_Domain, to_all+scalar_pair)
768 if (cs%linear_wave_drag) &
770 cs%BT_Domain, to_all+scalar_pair)
773 if (((g%isd > cs%isdw) .or. (g%jsd > cs%jsdw)) .or. (isq <= is-1) .or. (jsq <= js-1)) &
774 call create_group_pass(cs%pass_force_hbt0_Cor_ref, bt_force_u, bt_force_v, cs%BT_Domain)
775 if (add_uh0)
call create_group_pass(cs%pass_force_hbt0_Cor_ref, uhbt0, vhbt0, cs%BT_Domain)
776 call create_group_pass(cs%pass_force_hbt0_Cor_ref, cor_ref_u, cor_ref_v, cs%BT_Domain)
777 if (.not. use_bt_cont)
then
778 call create_group_pass(cs%pass_Dat_uv, datu, datv, cs%BT_Domain, to_all+scalar_pair)
793 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
799 if (cs%linearized_BT_PV)
then
801 do j=jsvf-2,jevf+1 ;
do i=isvf-2,ievf+1
805 do j=jsvf-1,jevf+1 ;
do i=isvf-2,ievf+1
806 dcor_u(i,j) = cs%D_u_Cor(i,j)
809 do j=jsvf-2,jevf+1 ;
do i=isvf-1,ievf+1
810 dcor_v(i,j) = cs%D_v_Cor(i,j)
813 q(:,:) = 0.0 ; dcor_u(:,:) = 0.0 ; dcor_v(:,:) = 0.0
817 do j=js,je ;
do i=is-1,ie
818 dcor_u(i,j) = 0.5 * (g%bathyT(i+1,j) + g%bathyT(i,j))
821 do j=js-1,je ;
do i=is,ie
822 dcor_v(i,j) = 0.5 * (g%bathyT(i,j+1) + g%bathyT(i,j))
825 do j=js-1,je ;
do i=is-1,ie
826 q(i,j) = 0.25 * g%CoriolisBu(i,j) * &
827 ((g%areaT(i,j) + g%areaT(i+1,j+1)) + (g%areaT(i+1,j) + g%areaT(i,j+1))) / &
828 ((g%areaT(i,j) * g%bathyT(i,j) + g%areaT(i+1,j+1) * g%bathyT(i+1,j+1)) + &
829 (g%areaT(i+1,j) * g%bathyT(i+1,j) + g%areaT(i,j+1) * g%bathyT(i,j+1)) )
836 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
837 if (nonblock_setup)
then
838 call start_group_pass(cs%pass_q_DCor, cs%BT_Domain, clock=id_clock_pass_pre)
840 call do_group_pass(cs%pass_q_DCor, cs%BT_Domain, clock=id_clock_pass_pre)
842 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
847 do j=cs%jsdw,cs%jedw ;
do i=cs%isdw,cs%iedw
848 gtot_e(i,j) = 0.0 ; gtot_w(i,j) = 0.0
849 gtot_n(i,j) = 0.0 ; gtot_s(i,j) = 0.0
852 if (interp_eta_pf)
then
853 eta_pf_1(i,j) = 0.0 ; d_eta_pf(i,j) = 0.0
855 p_surf_dyn(i,j) = 0.0
856 if (cs%dynamic_psurf) dyn_coef_eta(i,j) = 0.0
862 do j=cs%jsdw,cs%jedw ;
do i=cs%isdw-1,cs%iedw
863 cor_ref_u(i,j) = 0.0 ; bt_force_u(i,j) = 0.0 ; ubt(i,j) = 0.0
864 datu(i,j) = 0.0 ; bt_rem_u(i,j) = 0.0 ; uhbt0(i,j) = 0.0
867 do j=cs%jsdw-1,cs%jedw ;
do i=cs%isdw,cs%iedw
868 cor_ref_v(i,j) = 0.0 ; bt_force_v(i,j) = 0.0 ; vbt(i,j) = 0.0
869 datv(i,j) = 0.0 ; bt_rem_v(i,j) = 0.0 ; vhbt0(i,j) = 0.0
873 if (interp_eta_pf)
then
875 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
876 eta(i,j) = eta_in(i,j)
877 eta_pf_1(i,j) = eta_pf_start(i,j)
878 d_eta_pf(i,j) = eta_pf_in(i,j) - eta_pf_start(i,j)
882 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
883 eta(i,j) = eta_in(i,j)
884 eta_pf(i,j) = eta_pf_in(i,j)
889 do k=1,nz ;
do j=js,je ;
do i=is-1,ie
892 if (visc_rem_u(i,j,k) <= 0.0)
then ; visc_rem = 0.0
893 elseif (visc_rem_u(i,j,k) >= 1.0)
then ; visc_rem = 1.0
894 elseif (visc_rem_u(i,j,k)**2 > visc_rem_u(i,j,k) - 0.5*instep)
then
895 visc_rem = visc_rem_u(i,j,k)
896 else ; visc_rem = 1.0 - 0.5*instep/visc_rem_u(i,j,k) ;
endif
897 wt_u(i,j,k) = cs%frhatu(i,j,k) * visc_rem
898 enddo ;
enddo ;
enddo
900 do k=1,nz ;
do j=js-1,je ;
do i=is,ie
902 if (visc_rem_v(i,j,k) <= 0.0)
then ; visc_rem = 0.0
903 elseif (visc_rem_v(i,j,k) >= 1.0)
then ; visc_rem = 1.0
904 elseif (visc_rem_v(i,j,k)**2 > visc_rem_v(i,j,k) - 0.5*instep)
then
905 visc_rem = visc_rem_v(i,j,k)
906 else ; visc_rem = 1.0 - 0.5*instep/visc_rem_v(i,j,k) ;
endif
907 wt_v(i,j,k) = cs%frhatv(i,j,k) * visc_rem
908 enddo ;
enddo ;
enddo
913 do j=js-1,je+1 ;
do i=is-1,ie ; ubt_cor(i,j) = 0.0 ;
enddo ;
enddo
915 do j=js-1,je ;
do i=is-1,ie+1 ; vbt_cor(i,j) = 0.0 ;
enddo ;
enddo
917 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
918 ubt_cor(i,j) = ubt_cor(i,j) + wt_u(i,j,k) * us%m_s_to_L_T*u_cor(i,j,k)
919 enddo ;
enddo ;
enddo
921 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
922 vbt_cor(i,j) = vbt_cor(i,j) + wt_v(i,j,k) * us%m_s_to_L_T*v_cor(i,j,k)
923 enddo ;
enddo ;
enddo
931 do k=1,nz ;
do i=is-1,ie
932 gtot_e(i,j) = gtot_e(i,j) + pbce(i,j,k) * wt_u(i,j,k)
933 gtot_w(i+1,j) = gtot_w(i+1,j) + pbce(i+1,j,k) * wt_u(i,j,k)
938 do k=1,nz ;
do i=is,ie
939 gtot_n(i,j) = gtot_n(i,j) + pbce(i,j,k) * wt_v(i,j,k)
940 gtot_s(i,j+1) = gtot_s(i,j+1) + pbce(i,j+1,k) * wt_v(i,j,k)
945 call tidal_forcing_sensitivity(g, cs%tides_CSp, det_de)
946 dgeo_de = 1.0 + det_de + cs%G_extra
948 dgeo_de = 1.0 + cs%G_extra
951 if (nonblock_setup .and. .not.cs%linearized_BT_PV)
then
952 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
953 call complete_group_pass(cs%pass_q_DCor, cs%BT_Domain, clock=id_clock_pass_pre)
954 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
959 if (use_bt_cont)
then
960 call set_local_bt_cont_types(bt_cont, btcl_u, btcl_v, g, us, ms, cs%BT_Domain, 1+ievf-ie)
962 if (cs%Nonlinear_continuity)
then
963 call find_face_areas(datu, datv, g, gv, us, cs, ms, eta, 1)
965 call find_face_areas(datu, datv, g, gv, us, cs, ms, halo=1)
971 call set_up_bt_obc(obc, eta, cs%BT_OBC, cs%BT_Domain, g, gv, us, ms, ievf-ie, use_bt_cont, &
972 datu, datv, btcl_u, btcl_v)
982 do j=js,je ;
do i=is-1,ie
986 bt_force_u(i,j) = forces%taux(i,j) * mass_to_z * cs%IDatu(i,j)*visc_rem_u(i,j,1)
989 do j=js-1,je ;
do i=is,ie
993 bt_force_v(i,j) = forces%tauy(i,j) * mass_to_z * cs%IDatv(i,j)*visc_rem_v(i,j,1)
995 if (
present(taux_bot) .and.
present(tauy_bot))
then
996 if (
associated(taux_bot) .and.
associated(tauy_bot))
then
998 do j=js,je ;
do i=is-1,ie
999 bt_force_u(i,j) = bt_force_u(i,j) - taux_bot(i,j) * mass_to_z * cs%IDatu(i,j)
1002 do j=js-1,je ;
do i=is,ie
1003 bt_force_v(i,j) = bt_force_v(i,j) - tauy_bot(i,j) * mass_to_z * cs%IDatv(i,j)
1011 do j=js,je ;
do k=1,nz ;
do i=isq,ieq
1012 bt_force_u(i,j) = bt_force_u(i,j) + wt_u(i,j,k) * us%m_to_L*us%T_to_s**2*bc_accel_u(i,j,k)
1013 enddo ;
enddo ;
enddo
1015 do j=jsq,jeq ;
do k=1,nz ;
do i=is,ie
1016 bt_force_v(i,j) = bt_force_v(i,j) + wt_v(i,j,k) * us%m_to_L*us%T_to_s**2*bc_accel_v(i,j,k)
1017 enddo ;
enddo ;
enddo
1023 do j=js,je ;
do i=is-1,ie ; uhbt(i,j) = 0.0 ; ubt(i,j) = 0.0 ;
enddo ;
enddo
1025 do j=js-1,je ;
do i=is,ie ; vhbt(i,j) = 0.0 ; vbt(i,j) = 0.0 ;
enddo ;
enddo
1026 if (cs%visc_rem_u_uh0)
then
1028 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
1029 uhbt(i,j) = uhbt(i,j) + us%T_to_s*us%m_to_L**2*uh0(i,j,k)
1030 ubt(i,j) = ubt(i,j) + wt_u(i,j,k) * us%m_s_to_L_T*u_uh0(i,j,k)
1031 enddo ;
enddo ;
enddo
1033 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
1034 vhbt(i,j) = vhbt(i,j) + us%T_to_s*us%m_to_L**2*vh0(i,j,k)
1035 vbt(i,j) = vbt(i,j) + wt_v(i,j,k) * us%m_s_to_L_T*v_vh0(i,j,k)
1036 enddo ;
enddo ;
enddo
1039 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
1040 uhbt(i,j) = uhbt(i,j) + us%T_to_s*us%m_to_L**2*uh0(i,j,k)
1041 ubt(i,j) = ubt(i,j) + cs%frhatu(i,j,k) * us%m_s_to_L_T*u_uh0(i,j,k)
1042 enddo ;
enddo ;
enddo
1044 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
1045 vhbt(i,j) = vhbt(i,j) + us%T_to_s*us%m_to_L**2*vh0(i,j,k)
1046 vbt(i,j) = vbt(i,j) + cs%frhatv(i,j,k) * us%m_s_to_L_T*v_vh0(i,j,k)
1047 enddo ;
enddo ;
enddo
1049 if (use_bt_cont)
then
1050 if (cs%adjust_BT_cont)
then
1055 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1056 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
1057 call pass_vector(ubt, vbt, cs%BT_Domain, complete=.false., halo=1+ievf-ie)
1058 call pass_vector(uhbt, vhbt, cs%BT_Domain, complete=.true., halo=1+ievf-ie)
1059 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
1060 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1062 call adjust_local_bt_cont_types(ubt, uhbt, vbt, vhbt, btcl_u, btcl_v, &
1063 g, us, ms, 1+ievf-ie)
1066 do j=js,je ;
do i=is-1,ie
1067 uhbt0(i,j) = uhbt(i,j) - find_uhbt(ubt(i,j), btcl_u(i,j), us)
1070 do j=js-1,je ;
do i=is,ie
1071 vhbt0(i,j) = vhbt(i,j) - find_vhbt(vbt(i,j), btcl_v(i,j), us)
1075 do j=js,je ;
do i=is-1,ie
1076 uhbt0(i,j) = uhbt(i,j) - datu(i,j)*ubt(i,j)
1079 do j=js-1,je ;
do i=is,ie
1080 vhbt0(i,j) = vhbt(i,j) - datv(i,j)*vbt(i,j)
1083 if (cs%BT_OBC%apply_u_OBCs)
then
1085 do j=js,je ;
do i=is-1,ie ;
if (obc%segnum_u(i,j) /= obc_none)
then
1087 endif ;
enddo ;
enddo
1089 if (cs%BT_OBC%apply_v_OBCs)
then
1091 do j=js-1,je ;
do i=is,ie ;
if (obc%segnum_v(i,j) /= obc_none)
then
1093 endif ;
enddo ;
enddo
1099 do j=jsvf-1,jevf+1 ;
do i=isvf-2,ievf+1
1100 ubt(i,j) = 0.0 ; uhbt(i,j) = 0.0 ; u_accel_bt(i,j) = 0.0
1103 do j=jsvf-2,jevf+1 ;
do i=isvf-1,ievf+1
1104 vbt(i,j) = 0.0 ; vhbt(i,j) = 0.0 ; v_accel_bt(i,j) = 0.0
1107 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
1108 ubt(i,j) = ubt(i,j) + wt_u(i,j,k) * us%m_s_to_L_T*u_in(i,j,k)
1109 enddo ;
enddo ;
enddo
1111 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
1112 vbt(i,j) = vbt(i,j) + wt_v(i,j,k) * us%m_s_to_L_T*v_in(i,j,k)
1113 enddo ;
enddo ;
enddo
1115 do j=js,je ;
do i=is-1,ie
1116 if (abs(ubt(i,j)) < cs%vel_underflow) ubt(i,j) = 0.0
1119 do j=js-1,je ;
do i=is,ie
1120 if (abs(vbt(i,j)) < cs%vel_underflow) vbt(i,j) = 0.0
1123 if (apply_obcs)
then
1124 ubt_first(:,:) = ubt(:,:) ; vbt_first(:,:) = vbt(:,:)
1127 if (cs%gradual_BT_ICs)
then
1129 do j=js,je ;
do i=is-1,ie
1130 bt_force_u(i,j) = bt_force_u(i,j) + (ubt(i,j) - cs%ubt_IC(i,j)) * idt
1131 ubt(i,j) = cs%ubt_IC(i,j)
1132 if (abs(ubt(i,j)) < cs%vel_underflow) ubt(i,j) = 0.0
1135 do j=js-1,je ;
do i=is,ie
1136 bt_force_v(i,j) = bt_force_v(i,j) + (vbt(i,j) - cs%vbt_IC(i,j)) * idt
1137 vbt(i,j) = cs%vbt_IC(i,j)
1138 if (abs(vbt(i,j)) < cs%vel_underflow) vbt(i,j) = 0.0
1142 if ((isq > is-1) .or. (jsq > js-1))
then
1145 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1146 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
1147 tmp_u(:,:) = 0.0 ; tmp_v(:,:) = 0.0
1148 do j=js,je ;
do i=isq,ieq ; tmp_u(i,j) = bt_force_u(i,j) ;
enddo ;
enddo
1149 do j=jsq,jeq ;
do i=is,ie ; tmp_v(i,j) = bt_force_v(i,j) ;
enddo ;
enddo
1150 if (nonblock_setup)
then
1151 call start_group_pass(cs%pass_tmp_uv, g%Domain)
1153 call do_group_pass(cs%pass_tmp_uv, g%Domain)
1154 do j=jsd,jed ;
do i=isdb,iedb ; bt_force_u(i,j) = tmp_u(i,j) ;
enddo ;
enddo
1155 do j=jsdb,jedb ;
do i=isd,ied ; bt_force_v(i,j) = tmp_v(i,j) ;
enddo ;
enddo
1157 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
1158 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1161 if (nonblock_setup)
then
1162 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1163 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
1164 call start_group_pass(cs%pass_gtot, cs%BT_Domain)
1165 call start_group_pass(cs%pass_ubt_Cor, g%Domain)
1166 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
1167 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1172 do j=jsvf-1,jevf ;
do i=isvf-1,ievf+1
1173 if (cs%Sadourny)
then
1174 amer(i-1,j) = dcor_u(i-1,j) * q(i-1,j)
1175 bmer(i,j) = dcor_u(i,j) * q(i,j)
1176 cmer(i,j+1) = dcor_u(i,j+1) * q(i,j)
1177 dmer(i-1,j+1) = dcor_u(i-1,j+1) * q(i-1,j)
1179 amer(i-1,j) = dcor_u(i-1,j) * &
1180 ((q(i,j) + q(i-1,j-1)) + q(i-1,j)) / 3.0
1181 bmer(i,j) = dcor_u(i,j) * &
1182 (q(i,j) + (q(i-1,j) + q(i,j-1))) / 3.0
1183 cmer(i,j+1) = dcor_u(i,j+1) * &
1184 (q(i,j) + (q(i-1,j) + q(i,j+1))) / 3.0
1185 dmer(i-1,j+1) = dcor_u(i-1,j+1) * &
1186 ((q(i,j) + q(i-1,j+1)) + q(i-1,j)) / 3.0
1191 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf
1192 if (cs%Sadourny)
then
1193 azon(i,j) = dcor_v(i+1,j) * q(i,j)
1194 bzon(i,j) = dcor_v(i,j) * q(i,j)
1195 czon(i,j) = dcor_v(i,j-1) * q(i,j-1)
1196 dzon(i,j) = dcor_v(i+1,j-1) * q(i,j-1)
1198 azon(i,j) = dcor_v(i+1,j) * &
1199 (q(i,j) + (q(i+1,j) + q(i,j-1))) / 3.0
1200 bzon(i,j) = dcor_v(i,j) * &
1201 (q(i,j) + (q(i-1,j) + q(i,j-1))) / 3.0
1202 czon(i,j) = dcor_v(i,j-1) * &
1203 ((q(i,j) + q(i-1,j-1)) + q(i,j-1)) / 3.0
1204 dzon(i,j) = dcor_v(i+1,j-1) * &
1205 ((q(i,j) + q(i+1,j-1)) + q(i,j-1)) / 3.0
1210 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1211 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
1212 if (nonblock_setup)
then
1213 if ((isq > is-1) .or. (jsq > js-1))
then
1214 call complete_group_pass(cs%pass_tmp_uv, g%Domain)
1215 do j=jsd,jed ;
do i=isdb,iedb ; bt_force_u(i,j) = tmp_u(i,j) ;
enddo ;
enddo
1216 do j=jsdb,jedb ;
do i=isd,ied ; bt_force_v(i,j) = tmp_v(i,j) ;
enddo ;
enddo
1218 call complete_group_pass(cs%pass_gtot, cs%BT_Domain)
1219 call complete_group_pass(cs%pass_ubt_Cor, g%Domain)
1221 call do_group_pass(cs%pass_gtot, cs%BT_Domain)
1222 call do_group_pass(cs%pass_ubt_Cor, g%Domain)
1226 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf+1
1227 if (cs%ua_polarity(i,j) < 0.0)
call swap(gtot_e(i,j), gtot_w(i,j))
1228 if (cs%va_polarity(i,j) < 0.0)
call swap(gtot_n(i,j), gtot_s(i,j))
1232 do j=js,je ;
do i=is-1,ie
1234 ((azon(i,j) * vbt_cor(i+1,j) + czon(i,j) * vbt_cor(i ,j-1)) + &
1235 (bzon(i,j) * vbt_cor(i ,j) + dzon(i,j) * vbt_cor(i+1,j-1)))
1238 do j=js-1,je ;
do i=is,ie
1239 cor_ref_v(i,j) = -1.0 * &
1240 ((amer(i-1,j) * ubt_cor(i-1,j) + cmer(i ,j+1) * ubt_cor(i ,j+1)) + &
1241 (bmer(i ,j) * ubt_cor(i ,j) + dmer(i-1,j+1) * ubt_cor(i-1,j+1)))
1245 if (nonblock_setup)
then
1246 if (.not.use_bt_cont) &
1247 call start_group_pass(cs%pass_Dat_uv, cs%BT_Domain)
1250 call start_group_pass(cs%pass_force_hbt0_Cor_ref, cs%BT_Domain)
1252 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
1253 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1264 do j=js-1,je+1 ;
do i=is-1,ie ; av_rem_u(i,j) = 0.0 ;
enddo ;
enddo
1266 do j=js-1,je ;
do i=is-1,ie+1 ; av_rem_v(i,j) = 0.0 ;
enddo ;
enddo
1268 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
1269 av_rem_u(i,j) = av_rem_u(i,j) + cs%frhatu(i,j,k) * visc_rem_u(i,j,k)
1270 enddo ;
enddo ;
enddo
1272 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
1273 av_rem_v(i,j) = av_rem_v(i,j) + cs%frhatv(i,j,k) * visc_rem_v(i,j,k)
1274 enddo ;
enddo ;
enddo
1275 if (cs%strong_drag)
then
1277 do j=js,je ;
do i=is-1,ie
1278 bt_rem_u(i,j) = g%mask2dCu(i,j) * &
1279 ((nstep * av_rem_u(i,j)) / (1.0 + (nstep-1)*av_rem_u(i,j)))
1282 do j=js-1,je ;
do i=is,ie
1283 bt_rem_v(i,j) = g%mask2dCv(i,j) * &
1284 ((nstep * av_rem_v(i,j)) / (1.0 + (nstep-1)*av_rem_v(i,j)))
1288 do j=js,je ;
do i=is-1,ie
1290 if (g%mask2dCu(i,j) * av_rem_u(i,j) > 0.0) &
1291 bt_rem_u(i,j) = g%mask2dCu(i,j) * (av_rem_u(i,j)**instep)
1294 do j=js-1,je ;
do i=is,ie
1296 if (g%mask2dCv(i,j) * av_rem_v(i,j) > 0.0) &
1297 bt_rem_v(i,j) = g%mask2dCv(i,j) * (av_rem_v(i,j)**instep)
1300 if (cs%linear_wave_drag)
then
1302 do j=js,je ;
do i=is-1,ie ;
if (cs%lin_drag_u(i,j) > 0.0)
then
1303 htot = 0.5 * (eta(i,j) + eta(i+1,j))
1304 if (gv%Boussinesq) &
1305 htot = htot + 0.5*gv%Z_to_H * (cs%bathyT(i,j) + cs%bathyT(i+1,j))
1306 bt_rem_u(i,j) = bt_rem_u(i,j) * (htot / (htot + cs%lin_drag_u(i,j) * dtbt))
1308 rayleigh_u(i,j) = cs%lin_drag_u(i,j) / htot
1309 endif ;
enddo ;
enddo
1311 do j=js-1,je ;
do i=is,ie ;
if (cs%lin_drag_v(i,j) > 0.0)
then
1312 htot = 0.5 * (eta(i,j) + eta(i,j+1))
1313 if (gv%Boussinesq) &
1314 htot = htot + 0.5*gv%Z_to_H * (cs%bathyT(i,j) + cs%bathyT(i+1,j+1))
1315 bt_rem_v(i,j) = bt_rem_v(i,j) * (htot / (htot + cs%lin_drag_v(i,j) * dtbt))
1317 rayleigh_v(i,j) = cs%lin_drag_v(i,j) / htot
1318 endif ;
enddo ;
enddo
1322 if (find_etaav)
then
1324 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf+1
1325 eta_sum(i,j) = 0.0 ; eta_wtd(i,j) = 0.0
1329 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf+1
1334 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf
1335 ubt_sum(i,j) = 0.0 ; uhbt_sum(i,j) = 0.0
1336 pfu_bt_sum(i,j) = 0.0 ; coru_bt_sum(i,j) = 0.0
1337 ubt_wtd(i,j) = 0.0 ; ubt_trans(i,j) = 0.0
1340 do j=jsvf-1,jevf ;
do i=isvf-1,ievf+1
1341 vbt_sum(i,j) = 0.0 ; vhbt_sum(i,j) = 0.0
1342 pfv_bt_sum(i,j) = 0.0 ; corv_bt_sum(i,j) = 0.0
1343 vbt_wtd(i,j) = 0.0 ; vbt_trans(i,j) = 0.0
1348 do j=jsvf-1,jevf+1;
do i=isvf-1,ievf+1 ; eta_src(i,j) = 0.0 ;
enddo ;
enddo
1349 if (cs%bound_BT_corr)
then ;
if (use_bt_cont .and. cs%BT_cont_bounds)
then
1350 do j=js,je ;
do i=is,ie ;
if (g%mask2dT(i,j) > 0.0)
then
1351 if (cs%eta_cor(i,j) > 0.0)
then
1355 u_max_cor = us%m_to_L*g%dxT(i,j) * (cs%maxCFL_BT_cont*idt)
1356 v_max_cor = us%m_to_L*g%dyT(i,j) * (cs%maxCFL_BT_cont*idt)
1357 eta_cor_max = dt_in_t * (cs%IareaT(i,j) * &
1358 (((find_uhbt(u_max_cor, btcl_u(i,j), us) + uhbt0(i,j)) - &
1359 (find_uhbt(-u_max_cor, btcl_u(i-1,j), us) + uhbt0(i-1,j))) + &
1360 ((find_vhbt(v_max_cor, btcl_v(i,j), us) + vhbt0(i,j)) - &
1361 (find_vhbt(-v_max_cor, btcl_v(i,j-1), us) + vhbt0(i,j-1))) ))
1362 cs%eta_cor(i,j) = min(cs%eta_cor(i,j), max(0.0, eta_cor_max))
1367 if (gv%Boussinesq) htot = cs%bathyT(i,j)*gv%Z_to_H + eta(i,j)
1369 cs%eta_cor(i,j) = max(cs%eta_cor(i,j), -max(0.0,htot))
1371 endif ;
enddo ;
enddo
1372 else ;
do j=js,je ;
do i=is,ie
1373 if (abs(cs%eta_cor(i,j)) > dt_in_t*cs%eta_cor_bound(i,j)) &
1374 cs%eta_cor(i,j) = sign(dt_in_t*cs%eta_cor_bound(i,j), cs%eta_cor(i,j))
1375 enddo ;
enddo ;
endif ;
endif
1377 do j=js,je ;
do i=is,ie
1378 eta_src(i,j) = g%mask2dT(i,j) * (instep * cs%eta_cor(i,j))
1382 if (cs%dynamic_psurf)
then
1383 ice_is_rigid = (
associated(forces%rigidity_ice_u) .and. &
1384 associated(forces%rigidity_ice_v))
1385 h_min_dyn = gv%Z_to_H * cs%Dmin_dyn_psurf
1386 if (ice_is_rigid .and. use_bt_cont) &
1387 call bt_cont_to_face_areas(bt_cont, datu, datv, g, us, ms, 0, .true.)
1388 if (ice_is_rigid)
then
1390 do j=js,je ;
do i=is,ie
1396 idt_max2 = 0.5 * (dgeo_de * (1.0 + 2.0*bebt)) * (us%L_to_m**2*g%IareaT(i,j) * &
1397 ((gtot_e(i,j) * (datu(i,j)*us%L_to_m*g%IdxCu(i,j)) + &
1398 gtot_w(i,j) * (datu(i-1,j)*us%L_to_m*g%IdxCu(i-1,j))) + &
1399 (gtot_n(i,j) * (datv(i,j)*us%L_to_m*g%IdyCv(i,j)) + &
1400 gtot_s(i,j) * (datv(i,j-1)*us%L_to_m*g%IdyCv(i,j-1)))) + &
1401 ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
1402 (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2)))
1403 h_eff_dx2 = max(h_min_dyn * ((us%L_to_m*g%IdxT(i,j))**2 + (us%L_to_m*g%IdyT(i,j))**2), &
1404 us%L_to_m**2*g%IareaT(i,j) * &
1405 ((datu(i,j)*us%L_to_m*g%IdxCu(i,j) + datu(i-1,j)*us%L_to_m*g%IdxCu(i-1,j)) + &
1406 (datv(i,j)*us%L_to_m*g%IdyCv(i,j) + datv(i,j-1)*us%L_to_m*g%IdyCv(i,j-1)) ) )
1407 dyn_coef_max = cs%const_dyn_psurf * max(0.0, 1.0 - dtbt**2 * idt_max2) / &
1408 (dtbt**2 * h_eff_dx2)
1411 ice_strength = us%m_to_L**4*us%Z_to_m*us%T_to_s* &
1412 ((forces%rigidity_ice_u(i,j) + forces%rigidity_ice_u(i-1,j)) + &
1413 (forces%rigidity_ice_v(i,j) + forces%rigidity_ice_v(i,j-1))) / &
1414 (cs%ice_strength_length**2 * dtbt)
1417 dyn_coef_eta(i,j) = min(dyn_coef_max, ice_strength * gv%H_to_Z)
1418 enddo ;
enddo ;
endif
1421 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1422 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
1423 if (nonblock_setup)
then
1424 call start_group_pass(cs%pass_eta_bt_rem, cs%BT_Domain)
1427 call do_group_pass(cs%pass_eta_bt_rem, cs%BT_Domain)
1428 if (.not.use_bt_cont) &
1429 call do_group_pass(cs%pass_Dat_uv, cs%BT_Domain)
1430 call do_group_pass(cs%pass_force_hbt0_Cor_ref, cs%BT_Domain)
1432 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
1433 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1436 if (nonblock_setup)
then
1437 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1438 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
1440 if (.not.use_bt_cont) &
1441 call complete_group_pass(cs%pass_Dat_uv, cs%BT_Domain)
1442 call complete_group_pass(cs%pass_force_hbt0_Cor_ref, cs%BT_Domain)
1443 call complete_group_pass(cs%pass_eta_bt_rem, cs%BT_Domain)
1445 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
1446 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1450 call uvchksum(
"BT [uv]hbt", uhbt, vhbt, cs%debug_BT_HI, haloshift=0, &
1451 scale=us%s_to_T*us%L_to_m**2*gv%H_to_m)
1452 call uvchksum(
"BT Initial [uv]bt", ubt, vbt, cs%debug_BT_HI, haloshift=0, scale=us%L_T_to_m_s)
1453 call hchksum(eta,
"BT Initial eta", cs%debug_BT_HI, haloshift=0, scale=gv%H_to_m)
1454 call uvchksum(
"BT BT_force_[uv]", bt_force_u, bt_force_v, &
1455 cs%debug_BT_HI, haloshift=0, scale=us%L_T2_to_m_s2)
1456 if (interp_eta_pf)
then
1457 call hchksum(eta_pf_1,
"BT eta_PF_1",cs%debug_BT_HI,haloshift=0, scale=gv%H_to_m)
1458 call hchksum(d_eta_pf,
"BT d_eta_PF",cs%debug_BT_HI,haloshift=0, scale=gv%H_to_m)
1460 call hchksum(eta_pf,
"BT eta_PF",cs%debug_BT_HI,haloshift=0, scale=gv%H_to_m)
1461 call hchksum(eta_pf_in,
"BT eta_PF_in",g%HI,haloshift=0, scale=gv%H_to_m)
1463 call uvchksum(
"BT Cor_ref_[uv]", cor_ref_u, cor_ref_v, cs%debug_BT_HI, haloshift=0, scale=us%L_T2_to_m_s2)
1464 call uvchksum(
"BT [uv]hbt0", uhbt0, vhbt0, cs%debug_BT_HI, haloshift=0, &
1465 scale=us%L_to_m**2*us%s_to_T*gv%H_to_m)
1466 if (.not. use_bt_cont)
then
1467 call uvchksum(
"BT Dat[uv]", datu, datv, cs%debug_BT_HI, haloshift=1, scale=us%L_to_m*gv%H_to_m)
1469 call uvchksum(
"BT wt_[uv]", wt_u, wt_v, g%HI, 0, .true., .true.)
1470 call uvchksum(
"BT frhat[uv]", cs%frhatu, cs%frhatv, g%HI, 0, .true., .true.)
1471 call uvchksum(
"BT bc_accel_[uv]", bc_accel_u, bc_accel_v, g%HI, haloshift=0)
1472 call uvchksum(
"BT IDat[uv]", cs%IDatu, cs%IDatv, g%HI, haloshift=0, scale=us%m_to_Z)
1473 call uvchksum(
"BT visc_rem_[uv]", visc_rem_u, visc_rem_v, g%HI, haloshift=1)
1476 if (query_averaging_enabled(cs%diag))
then
1477 if (cs%id_eta_st > 0)
call post_data(cs%id_eta_st, eta(isd:ied,jsd:jed), cs%diag)
1478 if (cs%id_ubt_st > 0)
call post_data(cs%id_ubt_st, ubt(isdb:iedb,jsd:jed), cs%diag)
1479 if (cs%id_vbt_st > 0)
call post_data(cs%id_vbt_st, vbt(isd:ied,jsdb:jedb), cs%diag)
1482 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1483 if (id_clock_calc > 0)
call cpu_clock_begin(id_clock_calc)
1485 if (project_velocity)
then ; eta_pf_bt => eta ;
else ; eta_pf_bt => eta_pred ;
endif
1487 if (cs%dt_bt_filter >= 0.0)
then
1488 dt_filt = 0.5 * max(0.0, min(cs%dt_bt_filter, 2.0*dt_in_t))
1490 dt_filt = 0.5 * max(0.0, dt_in_t * min(-cs%dt_bt_filter, 2.0))
1492 nfilter = ceiling(dt_filt / dtbt)
1494 if (nstep+nfilter==0 )
call mom_error(fatal, &
1495 "btstep: number of barotropic step (nstep+nfilter) is 0")
1498 sum_wt_vel = 0.0 ; sum_wt_eta = 0.0 ; sum_wt_accel = 0.0 ; sum_wt_trans = 0.0
1499 allocate(wt_vel(nstep+nfilter)) ;
allocate(wt_eta(nstep+nfilter))
1500 allocate(wt_trans(nstep+nfilter+1)) ;
allocate(wt_accel(nstep+nfilter+1))
1501 allocate(wt_accel2(nstep+nfilter+1))
1502 do n=1,nstep+nfilter
1504 if ( (n==nstep) .or. (dt_filt - abs(n-nstep)*dtbt >= 0.0))
then
1505 wt_vel(n) = 1.0 ; wt_eta(n) = 1.0
1506 elseif (dtbt + dt_filt - abs(n-nstep)*dtbt > 0.0)
then
1507 wt_vel(n) = 1.0 + (dt_filt / dtbt) - abs(n-nstep) ; wt_eta(n) = wt_vel(n)
1509 wt_vel(n) = 0.0 ; wt_eta(n) = 0.0
1515 sum_wt_vel = sum_wt_vel + wt_vel(n) ; sum_wt_eta = sum_wt_eta + wt_eta(n)
1517 wt_trans(nstep+nfilter+1) = 0.0 ; wt_accel(nstep+nfilter+1) = 0.0
1518 do n=nstep+nfilter,1,-1
1519 wt_trans(n) = wt_trans(n+1) + wt_eta(n)
1520 wt_accel(n) = wt_accel(n+1) + wt_vel(n)
1521 sum_wt_accel = sum_wt_accel + wt_accel(n) ; sum_wt_trans = sum_wt_trans + wt_trans(n)
1524 i_sum_wt_vel = 1.0 / sum_wt_vel ; i_sum_wt_accel = 1.0 / sum_wt_accel
1525 i_sum_wt_eta = 1.0 / sum_wt_eta ; i_sum_wt_trans = 1.0 / sum_wt_trans
1526 do n=1,nstep+nfilter
1527 wt_vel(n) = wt_vel(n) * i_sum_wt_vel
1528 wt_accel2(n) = wt_accel(n)
1529 wt_accel(n) = wt_accel(n) * i_sum_wt_accel
1530 wt_eta(n) = wt_eta(n) * i_sum_wt_eta
1535 sum_wt_vel = 0.0 ; sum_wt_eta = 0.0 ; sum_wt_accel = 0.0 ; sum_wt_trans = 0.0
1538 isv=is ; iev=ie ; jsv=js ; jev=je
1539 do n=1,nstep+nfilter
1541 sum_wt_vel = sum_wt_vel + wt_vel(n)
1542 sum_wt_eta = sum_wt_eta + wt_eta(n)
1543 sum_wt_accel = sum_wt_accel + wt_accel2(n)
1544 sum_wt_trans = sum_wt_trans + wt_trans(n)
1546 if (cs%clip_velocity)
then
1547 do j=jsv,jev ;
do i=isv-1,iev
1548 if ((ubt(i,j) * (dt_in_t * us%m_to_L*g%dy_Cu(i,j))) * us%L_to_m**2*g%IareaT(i+1,j) < -cs%CFL_trunc)
then
1550 ubt(i,j) = (-0.95*cs%CFL_trunc) * (us%m_to_L**2*g%areaT(i+1,j) / (dt_in_t * us%m_to_L*g%dy_Cu(i,j)))
1551 elseif ((ubt(i,j) * (dt_in_t * us%m_to_L*g%dy_Cu(i,j))) * us%L_to_m**2*g%IareaT(i,j) > cs%CFL_trunc)
then
1553 ubt(i,j) = (0.95*cs%CFL_trunc) * (us%m_to_L**2*g%areaT(i,j) / (dt_in_t * us%m_to_L*g%dy_Cu(i,j)))
1556 do j=jsv-1,jev ;
do i=isv,iev
1557 if ((vbt(i,j) * (dt_in_t * us%m_to_L*g%dx_Cv(i,j))) * us%L_to_m**2*g%IareaT(i,j+1) < -cs%CFL_trunc)
then
1559 vbt(i,j) = (-0.9*cs%CFL_trunc) * (us%m_to_L**2*g%areaT(i,j+1) / (dt_in_t * us%m_to_L*g%dx_Cv(i,j)))
1560 elseif ((vbt(i,j) * (dt_in_t * us%m_to_L*g%dx_Cv(i,j))) * us%L_to_m**2*g%IareaT(i,j) > cs%CFL_trunc)
then
1562 vbt(i,j) = (0.9*cs%CFL_trunc) * (us%m_to_L**2*g%areaT(i,j) / (dt_in_t * us%m_to_L*g%dx_Cv(i,j)))
1567 if ((iev - stencil < ie) .or. (jev - stencil < je))
then
1568 if (id_clock_calc > 0)
call cpu_clock_end(id_clock_calc)
1569 call do_group_pass(cs%pass_eta_ubt, cs%BT_Domain, clock=id_clock_pass_step)
1570 isv = isvf ; iev = ievf ; jsv = jsvf ; jev = jevf
1571 if (id_clock_calc > 0)
call cpu_clock_begin(id_clock_calc)
1573 isv = isv+stencil ; iev = iev-stencil
1574 jsv = jsv+stencil ; jev = jev-stencil
1577 if ((.not.use_bt_cont) .and. cs%Nonlinear_continuity .and. &
1578 (cs%Nonlin_cont_update_period > 0))
then
1579 if ((n>1) .and. (mod(n-1,cs%Nonlin_cont_update_period) == 0)) &
1580 call find_face_areas(datu, datv, g, gv, us, cs, ms, eta, 1+iev-ie)
1584 if (cs%dynamic_psurf .or. .not.project_velocity)
then
1585 if (use_bt_cont)
then
1587 do j=jsv-1,jev+1 ;
do i=isv-2,iev+1
1588 uhbt(i,j) = find_uhbt(ubt(i,j), btcl_u(i,j), us) + uhbt0(i,j)
1591 do j=jsv-2,jev+1 ;
do i=isv-1,iev+1
1592 vhbt(i,j) = find_vhbt(vbt(i,j), btcl_v(i,j), us) + vhbt0(i,j)
1595 do j=jsv-1,jev+1 ;
do i=isv-1,iev+1
1596 eta_pred(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * cs%IareaT(i,j)) * &
1597 ((uhbt(i-1,j) - uhbt(i,j)) + (vhbt(i,j-1) - vhbt(i,j)))
1601 do j=jsv-1,jev+1 ;
do i=isv-1,iev+1
1602 eta_pred(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * cs%IareaT(i,j)) * &
1603 (((datu(i-1,j)*ubt(i-1,j) + uhbt0(i-1,j)) - &
1604 (datu(i,j)*ubt(i,j) + uhbt0(i,j))) + &
1605 ((datv(i,j-1)*vbt(i,j-1) + vhbt0(i,j-1)) - &
1606 (datv(i,j)*vbt(i,j) + vhbt0(i,j))))
1610 if (cs%dynamic_psurf)
then
1612 do j=jsv-1,jev+1 ;
do i=isv-1,iev+1
1613 p_surf_dyn(i,j) = dyn_coef_eta(i,j) * (eta_pred(i,j) - eta(i,j))
1621 if (find_etaav)
then
1623 do j=js,je ;
do i=is,ie
1624 eta_sum(i,j) = eta_sum(i,j) + wt_accel2(n) * eta_pf_bt(i,j)
1628 if (interp_eta_pf)
then
1631 do j=jsv-1,jev+1 ;
do i=isv-1,iev+1
1632 eta_pf(i,j) = eta_pf_1(i,j) + wt_end*d_eta_pf(i,j)
1636 if (apply_obc_flather .or. apply_obc_open)
then
1638 do j=jsv,jev ;
do i=isv-2,iev+1
1639 ubt_old(i,j) = ubt(i,j)
1642 do j=jsv-2,jev+1 ;
do i=isv,iev
1643 vbt_old(i,j) = vbt(i,j)
1648 if (apply_obcs)
then
1649 if (mod(n+g%first_direction,2)==1)
then
1655 if (cs%BT_OBC%apply_u_OBCs)
then
1657 do j=jsv-joff,jev+joff ;
do i=isv-1,iev
1658 ubt_prev(i,j) = ubt(i,j) ; uhbt_prev(i,j) = uhbt(i,j)
1659 ubt_sum_prev(i,j) = ubt_sum(i,j) ; uhbt_sum_prev(i,j) = uhbt_sum(i,j) ; ubt_wtd_prev(i,j) = ubt_wtd(i,j)
1663 if (cs%BT_OBC%apply_v_OBCs)
then
1665 do j=jsv-1,jev ;
do i=isv-ioff,iev+ioff
1666 vbt_prev(i,j) = vbt(i,j) ; vhbt_prev(i,j) = vhbt(i,j)
1667 vbt_sum_prev(i,j) = vbt_sum(i,j) ; vhbt_sum_prev(i,j) = vhbt_sum(i,j) ; vbt_wtd_prev(i,j) = vbt_wtd(i,j)
1673 if (mod(n+g%first_direction,2)==1)
then
1676 do j=jsv-1,jev ;
do i=isv-1,iev+1
1677 cor_v(i,j) = -1.0*((amer(i-1,j) * ubt(i-1,j) + cmer(i,j+1) * ubt(i,j+1)) + &
1678 (bmer(i,j) * ubt(i,j) + dmer(i-1,j+1) * ubt(i-1,j+1))) - cor_ref_v(i,j)
1679 pfv(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_n(i,j) - &
1680 (eta_pf_bt(i,j+1)-eta_pf(i,j+1))*gtot_s(i,j+1)) * &
1681 dgeo_de * cs%IdyCv(i,j)
1683 if (cs%dynamic_psurf)
then
1685 do j=jsv-1,jev ;
do i=isv-1,iev+1
1686 pfv(i,j) = pfv(i,j) + (p_surf_dyn(i,j) - p_surf_dyn(i,j+1)) * cs%IdyCv(i,j)
1690 if (cs%BT_OBC%apply_v_OBCs)
then
1692 do j=jsv-1,jev ;
do i=isv-1,iev+1 ;
if (obc%segnum_v(i,j) /= obc_none)
then
1694 endif ;
enddo ;
enddo
1697 do j=jsv-1,jev ;
do i=isv-1,iev+1
1699 vbt(i,j) = bt_rem_v(i,j) * (vbt(i,j) + &
1700 dtbt * ((bt_force_v(i,j) + cor_v(i,j)) + pfv(i,j)))
1701 vbt_trans(i,j) = trans_wt1*vbt(i,j) + trans_wt2*vel_prev
1703 if (cs%linear_wave_drag)
then
1704 v_accel_bt(i,j) = v_accel_bt(i,j) + wt_accel(n) * &
1705 ((cor_v(i,j) + pfv(i,j)) - vbt(i,j)*rayleigh_v(i,j))
1707 v_accel_bt(i,j) = v_accel_bt(i,j) + wt_accel(n) * (cor_v(i,j) + pfv(i,j))
1711 if (use_bt_cont)
then
1713 do j=jsv-1,jev ;
do i=isv-1,iev+1
1714 vhbt(i,j) = find_vhbt(vbt_trans(i,j), btcl_v(i,j), us) + vhbt0(i,j)
1718 do j=jsv-1,jev ;
do i=isv-1,iev+1
1719 vhbt(i,j) = datv(i,j)*vbt_trans(i,j) + vhbt0(i,j)
1722 if (cs%BT_OBC%apply_v_OBCs)
then
1724 do j=jsv-1,jev ;
do i=isv-1,iev+1 ;
if (obc%segnum_v(i,j) /= obc_none)
then
1725 vbt(i,j) = vbt_prev(i,j) ; vhbt(i,j) = vhbt_prev(i,j)
1726 endif ;
enddo ;
enddo
1730 do j=jsv,jev ;
do i=isv-1,iev
1731 cor_u(i,j) = ((azon(i,j) * vbt(i+1,j) + czon(i,j) * vbt(i,j-1)) + &
1732 (bzon(i,j) * vbt(i,j) + dzon(i,j) * vbt(i+1,j-1))) - &
1734 pfu(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_e(i,j) - &
1735 (eta_pf_bt(i+1,j)-eta_pf(i+1,j))*gtot_w(i+1,j)) * &
1736 dgeo_de * cs%IdxCu(i,j)
1739 if (cs%dynamic_psurf)
then
1741 do j=jsv,jev ;
do i=isv-1,iev
1742 pfu(i,j) = pfu(i,j) + (p_surf_dyn(i,j) - p_surf_dyn(i+1,j)) * cs%IdxCu(i,j)
1746 if (cs%BT_OBC%apply_u_OBCs)
then
1748 do j=jsv,jev ;
do i=isv-1,iev ;
if (obc%segnum_u(i,j) /= obc_none)
then
1750 endif ;
enddo ;
enddo
1753 do j=jsv,jev ;
do i=isv-1,iev
1755 ubt(i,j) = bt_rem_u(i,j) * (ubt(i,j) + &
1756 dtbt * ((bt_force_u(i,j) + cor_u(i,j)) + pfu(i,j)))
1757 if (abs(ubt(i,j)) < cs%vel_underflow) ubt(i,j) = 0.0
1758 ubt_trans(i,j) = trans_wt1*ubt(i,j) + trans_wt2*vel_prev
1760 if (cs%linear_wave_drag)
then
1761 u_accel_bt(i,j) = u_accel_bt(i,j) + wt_accel(n) * &
1762 ((cor_u(i,j) + pfu(i,j)) - ubt(i,j)*rayleigh_u(i,j))
1764 u_accel_bt(i,j) = u_accel_bt(i,j) + wt_accel(n) * (cor_u(i,j) + pfu(i,j))
1768 if (use_bt_cont)
then
1770 do j=jsv,jev ;
do i=isv-1,iev
1771 uhbt(i,j) = find_uhbt(ubt_trans(i,j), btcl_u(i,j), us) + uhbt0(i,j)
1775 do j=jsv,jev ;
do i=isv-1,iev
1776 uhbt(i,j) = datu(i,j)*ubt_trans(i,j) + uhbt0(i,j)
1779 if (cs%BT_OBC%apply_u_OBCs)
then
1781 do j=jsv,jev ;
do i=isv-1,iev ;
if (obc%segnum_u(i,j) /= obc_none)
then
1782 ubt(i,j) = ubt_prev(i,j); uhbt(i,j) = uhbt_prev(i,j)
1783 endif ;
enddo ;
enddo
1788 do j=jsv-1,jev+1 ;
do i=isv-1,iev
1789 cor_u(i,j) = ((azon(i,j) * vbt(i+1,j) + czon(i,j) * vbt(i,j-1)) + &
1790 (bzon(i,j) * vbt(i,j) + dzon(i,j) * vbt(i+1,j-1))) - &
1792 pfu(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_e(i,j) - &
1793 (eta_pf_bt(i+1,j)-eta_pf(i+1,j))*gtot_w(i+1,j)) * &
1794 dgeo_de * cs%IdxCu(i,j)
1797 if (cs%dynamic_psurf)
then
1799 do j=jsv-1,jev+1 ;
do i=isv-1,iev
1800 pfu(i,j) = pfu(i,j) + (p_surf_dyn(i,j) - p_surf_dyn(i+1,j)) * cs%IdxCu(i,j)
1804 if (cs%BT_OBC%apply_u_OBCs)
then
1806 do j=jsv,jev ;
do i=isv-1,iev ;
if (obc%segnum_u(i,j) /= obc_none)
then
1808 endif ;
enddo ;
enddo
1812 do j=jsv-1,jev+1 ;
do i=isv-1,iev
1814 ubt(i,j) = bt_rem_u(i,j) * (ubt(i,j) + &
1815 dtbt * ((bt_force_u(i,j) + cor_u(i,j)) + pfu(i,j)))
1816 if (abs(ubt(i,j)) < cs%vel_underflow) ubt(i,j) = 0.0
1817 ubt_trans(i,j) = trans_wt1*ubt(i,j) + trans_wt2*vel_prev
1819 if (cs%linear_wave_drag)
then
1820 u_accel_bt(i,j) = u_accel_bt(i,j) + wt_accel(n) * &
1821 ((cor_u(i,j) + pfu(i,j)) - ubt(i,j)*rayleigh_u(i,j))
1823 u_accel_bt(i,j) = u_accel_bt(i,j) + wt_accel(n) * (cor_u(i,j) + pfu(i,j))
1827 if (use_bt_cont)
then
1829 do j=jsv-1,jev+1 ;
do i=isv-1,iev
1830 uhbt(i,j) = find_uhbt(ubt_trans(i,j), btcl_u(i,j), us) + uhbt0(i,j)
1834 do j=jsv-1,jev+1 ;
do i=isv-1,iev
1835 uhbt(i,j) = datu(i,j)*ubt_trans(i,j) + uhbt0(i,j)
1838 if (cs%BT_OBC%apply_u_OBCs)
then
1840 do j=jsv-1,jev+1 ;
do i=isv-1,iev ;
if (obc%segnum_u(i,j) /= obc_none)
then
1841 ubt(i,j) = ubt_prev(i,j); uhbt(i,j) = uhbt_prev(i,j)
1842 endif ;
enddo ;
enddo
1846 if (cs%use_old_coriolis_bracket_bug)
then
1848 do j=jsv-1,jev ;
do i=isv,iev
1849 cor_v(i,j) = -1.0*((amer(i-1,j) * ubt(i-1,j) + bmer(i,j) * ubt(i,j)) + &
1850 (cmer(i,j+1) * ubt(i,j+1) + dmer(i-1,j+1) * ubt(i-1,j+1))) - cor_ref_v(i,j)
1851 pfv(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_n(i,j) - &
1852 (eta_pf_bt(i,j+1)-eta_pf(i,j+1))*gtot_s(i,j+1)) * &
1853 dgeo_de * cs%IdyCv(i,j)
1857 do j=jsv-1,jev ;
do i=isv,iev
1858 cor_v(i,j) = -1.0*((amer(i-1,j) * ubt(i-1,j) + cmer(i,j+1) * ubt(i,j+1)) + &
1859 (bmer(i,j) * ubt(i,j) + dmer(i-1,j+1) * ubt(i-1,j+1))) - cor_ref_v(i,j)
1860 pfv(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_n(i,j) - &
1861 (eta_pf_bt(i,j+1)-eta_pf(i,j+1))*gtot_s(i,j+1)) * &
1862 dgeo_de * cs%IdyCv(i,j)
1866 if (cs%dynamic_psurf)
then
1868 do j=jsv-1,jev ;
do i=isv,iev
1869 pfv(i,j) = pfv(i,j) + (p_surf_dyn(i,j) - p_surf_dyn(i,j+1)) * cs%IdyCv(i,j)
1873 if (cs%BT_OBC%apply_v_OBCs)
then
1875 do j=jsv-1,jev ;
do i=isv-1,iev+1 ;
if (obc%segnum_v(i,j) /= obc_none)
then
1877 endif ;
enddo ;
enddo
1881 do j=jsv-1,jev ;
do i=isv,iev
1883 vbt(i,j) = bt_rem_v(i,j) * (vbt(i,j) + &
1884 dtbt * ((bt_force_v(i,j) + cor_v(i,j)) + pfv(i,j)))
1885 if (abs(vbt(i,j)) < cs%vel_underflow) vbt(i,j) = 0.0
1886 vbt_trans(i,j) = trans_wt1*vbt(i,j) + trans_wt2*vel_prev
1888 if (cs%linear_wave_drag)
then
1889 v_accel_bt(i,j) = v_accel_bt(i,j) + wt_accel(n) * &
1890 ((cor_v(i,j) + pfv(i,j)) - vbt(i,j)*rayleigh_v(i,j))
1892 v_accel_bt(i,j) = v_accel_bt(i,j) + wt_accel(n) * (cor_v(i,j) + pfv(i,j))
1895 if (use_bt_cont)
then
1897 do j=jsv-1,jev ;
do i=isv,iev
1898 vhbt(i,j) = find_vhbt(vbt_trans(i,j), btcl_v(i,j), us) + vhbt0(i,j)
1902 do j=jsv-1,jev ;
do i=isv,iev
1903 vhbt(i,j) = datv(i,j)*vbt_trans(i,j) + vhbt0(i,j)
1906 if (cs%BT_OBC%apply_v_OBCs)
then
1908 do j=jsv-1,jev ;
do i=isv,iev ;
if (obc%segnum_v(i,j) /= obc_none)
then
1909 vbt(i,j) = vbt_prev(i,j); vhbt(i,j) = vhbt_prev(i,j)
1910 endif ;
enddo ;
enddo
1918 do j=js,je ;
do i=is-1,ie
1919 pfu_bt_sum(i,j) = pfu_bt_sum(i,j) + wt_accel2(n) * pfu(i,j)
1922 do j=js-1,je ;
do i=is,ie
1923 pfv_bt_sum(i,j) = pfv_bt_sum(i,j) + wt_accel2(n) * pfv(i,j)
1928 do j=js,je ;
do i=is-1,ie
1929 coru_bt_sum(i,j) = coru_bt_sum(i,j) + wt_accel2(n) * cor_u(i,j)
1932 do j=js-1,je ;
do i=is,ie
1933 corv_bt_sum(i,j) = corv_bt_sum(i,j) + wt_accel2(n) * cor_v(i,j)
1938 do j=js,je ;
do i=is-1,ie
1939 ubt_sum(i,j) = ubt_sum(i,j) + wt_trans(n) * ubt_trans(i,j)
1940 uhbt_sum(i,j) = uhbt_sum(i,j) + wt_trans(n) * uhbt(i,j)
1941 ubt_wtd(i,j) = ubt_wtd(i,j) + wt_vel(n) * ubt(i,j)
1944 do j=js-1,je ;
do i=is,ie
1945 vbt_sum(i,j) = vbt_sum(i,j) + wt_trans(n) * vbt_trans(i,j)
1946 vhbt_sum(i,j) = vhbt_sum(i,j) + wt_trans(n) * vhbt(i,j)
1947 vbt_wtd(i,j) = vbt_wtd(i,j) + wt_vel(n) * vbt(i,j)
1951 if (apply_obcs)
then
1952 if (cs%BT_OBC%apply_u_OBCs)
then
1954 do j=js,je ;
do i=is-1,ie
1955 if (obc%segnum_u(i,j) /= obc_none)
then
1956 ubt_sum(i,j) = ubt_sum_prev(i,j) ; uhbt_sum(i,j) = uhbt_sum_prev(i,j)
1957 ubt_wtd(i,j) = ubt_wtd_prev(i,j)
1962 if (cs%BT_OBC%apply_v_OBCs)
then
1964 do j=js-1,je ;
do i=is,ie
1965 if (obc%segnum_v(i,j) /= obc_none)
then
1966 vbt_sum(i,j) = vbt_sum_prev(i,j) ; vhbt_sum(i,j) = vhbt_sum_prev(i,j)
1967 vbt_wtd(i,j) = vbt_wtd_prev(i,j)
1972 call apply_velocity_obcs(obc, ubt, vbt, uhbt, vhbt, &
1973 ubt_trans, vbt_trans, eta, ubt_old, vbt_old, cs%BT_OBC, &
1974 g, ms, us, iev-ie, dtbt, bebt, use_bt_cont, datu, datv, btcl_u, btcl_v, &
1976 if (cs%BT_OBC%apply_u_OBCs)
then ;
do j=js,je ;
do i=is-1,ie
1977 if (obc%segnum_u(i,j) /= obc_none)
then
1978 ubt_sum(i,j) = ubt_sum(i,j) + wt_trans(n) * ubt_trans(i,j)
1979 uhbt_sum(i,j) = uhbt_sum(i,j) + wt_trans(n) * uhbt(i,j)
1980 ubt_wtd(i,j) = ubt_wtd(i,j) + wt_vel(n) * ubt(i,j)
1982 enddo ;
enddo ;
endif
1983 if (cs%BT_OBC%apply_v_OBCs)
then ;
do j=js-1,je ;
do i=is,ie
1984 if (obc%segnum_v(i,j) /= obc_none)
then
1985 vbt_sum(i,j) = vbt_sum(i,j) + wt_trans(n) * vbt_trans(i,j)
1986 vhbt_sum(i,j) = vhbt_sum(i,j) + wt_trans(n) * vhbt(i,j)
1987 vbt_wtd(i,j) = vbt_wtd(i,j) + wt_vel(n) * vbt(i,j)
1989 enddo ;
enddo ;
endif
1992 if (cs%debug_bt)
then
1993 call uvchksum(
"BT [uv]hbt just after OBC", uhbt, vhbt, cs%debug_BT_HI, haloshift=iev-ie, &
1994 scale=us%s_to_T*us%L_to_m**2*gv%H_to_m)
1998 do j=jsv,jev ;
do i=isv,iev
1999 eta(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * cs%IareaT(i,j)) * &
2000 ((uhbt(i-1,j) - uhbt(i,j)) + (vhbt(i,j-1) - vhbt(i,j)))
2001 eta_wtd(i,j) = eta_wtd(i,j) + eta(i,j) * wt_eta(n)
2005 if (do_hifreq_output)
then
2006 time_step_end = time_bt_start + real_to_time(n*us%T_to_s*dtbt)
2007 call enable_averaging(us%T_to_s*dtbt, time_step_end, cs%diag)
2008 if (cs%id_ubt_hifreq > 0)
call post_data(cs%id_ubt_hifreq, ubt(isdb:iedb,jsd:jed), cs%diag)
2009 if (cs%id_vbt_hifreq > 0)
call post_data(cs%id_vbt_hifreq, vbt(isd:ied,jsdb:jedb), cs%diag)
2010 if (cs%id_eta_hifreq > 0)
call post_data(cs%id_eta_hifreq, eta(isd:ied,jsd:jed), cs%diag)
2011 if (cs%id_uhbt_hifreq > 0)
call post_data(cs%id_uhbt_hifreq, uhbt(isdb:iedb,jsd:jed), cs%diag)
2012 if (cs%id_vhbt_hifreq > 0)
call post_data(cs%id_vhbt_hifreq, vhbt(isd:ied,jsdb:jedb), cs%diag)
2013 if (cs%id_eta_pred_hifreq > 0)
call post_data(cs%id_eta_pred_hifreq, eta_pf_bt(isd:ied,jsd:jed), cs%diag)
2016 if (cs%debug_bt)
then
2017 write(mesg,
'("BT step ",I4)') n
2018 call uvchksum(trim(mesg)//
" [uv]bt", ubt, vbt, cs%debug_BT_HI, haloshift=iev-ie, &
2019 scale=us%L_T_to_m_s)
2020 call hchksum(eta, trim(mesg)//
" eta", cs%debug_BT_HI, haloshift=iev-ie, scale=gv%H_to_m)
2024 if (id_clock_calc > 0)
call cpu_clock_end(id_clock_calc)
2025 if (id_clock_calc_post > 0)
call cpu_clock_begin(id_clock_calc_post)
2028 if (do_hifreq_output)
call enable_averaging(time_int_in, time_end_in, cs%diag)
2030 i_sum_wt_vel = 1.0 / sum_wt_vel ; i_sum_wt_eta = 1.0 / sum_wt_eta
2031 i_sum_wt_accel = 1.0 / sum_wt_accel ; i_sum_wt_trans = 1.0 / sum_wt_trans
2033 if (find_etaav)
then ;
do j=js,je ;
do i=is,ie
2034 etaav(i,j) = eta_sum(i,j) * i_sum_wt_accel
2035 enddo ;
enddo ;
endif
2036 do j=js-1,je+1 ;
do i=is-1,ie+1 ; e_anom(i,j) = 0.0 ;
enddo ;
enddo
2037 if (interp_eta_pf)
then
2038 do j=js,je ;
do i=is,ie
2039 e_anom(i,j) = dgeo_de * (0.5 * (eta(i,j) + eta_in(i,j)) - &
2040 (eta_pf_1(i,j) + 0.5*d_eta_pf(i,j)))
2043 do j=js,je ;
do i=is,ie
2044 e_anom(i,j) = dgeo_de * (0.5 * (eta(i,j) + eta_in(i,j)) - eta_pf(i,j))
2047 if (apply_obcs)
then
2049 if (cs%BT_OBC%apply_u_OBCs)
then
2051 do j=js,je ;
do i=is-1,ie
2052 if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e)
then
2053 e_anom(i+1,j) = e_anom(i,j)
2054 elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w)
then
2055 e_anom(i,j) = e_anom(i+1,j)
2060 if (cs%BT_OBC%apply_v_OBCs)
then
2062 do j=js-1,je ;
do i=is,ie
2063 if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n)
then
2064 e_anom(i,j+1) = e_anom(i,j)
2065 elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s)
then
2066 e_anom(i,j) = e_anom(i,j+1)
2073 do j=js,je ;
do i=is,ie
2074 eta_out(i,j) = eta_wtd(i,j) * i_sum_wt_eta
2077 if (id_clock_calc_post > 0)
call cpu_clock_end(id_clock_calc_post)
2078 if (id_clock_pass_post > 0)
call cpu_clock_begin(id_clock_pass_post)
2079 if (g%nonblocking_updates)
then
2080 call start_group_pass(cs%pass_e_anom, g%Domain)
2082 if (find_etaav)
call do_group_pass(cs%pass_etaav, g%Domain)
2083 call do_group_pass(cs%pass_e_anom, g%Domain)
2085 if (id_clock_pass_post > 0)
call cpu_clock_end(id_clock_pass_post)
2086 if (id_clock_calc_post > 0)
call cpu_clock_begin(id_clock_calc_post)
2088 do j=js,je ;
do i=is-1,ie
2089 cs%ubtav(i,j) = ubt_sum(i,j) * i_sum_wt_trans
2090 uhbtav(i,j) = us%s_to_T*us%L_to_m**2*uhbt_sum(i,j) * i_sum_wt_trans
2093 ubt_wtd(i,j) = ubt_wtd(i,j) * i_sum_wt_vel
2096 do j=js-1,je ;
do i=is,ie
2097 cs%vbtav(i,j) = vbt_sum(i,j) * i_sum_wt_trans
2098 vhbtav(i,j) = us%s_to_T*us%L_to_m**2*vhbt_sum(i,j) * i_sum_wt_trans
2101 vbt_wtd(i,j) = vbt_wtd(i,j) * i_sum_wt_vel
2104 if (id_clock_calc_post > 0)
call cpu_clock_end(id_clock_calc_post)
2105 if (id_clock_pass_post > 0)
call cpu_clock_begin(id_clock_pass_post)
2106 if (g%nonblocking_updates)
then
2107 call complete_group_pass(cs%pass_e_anom, g%Domain)
2108 if (find_etaav)
call start_group_pass(cs%pass_etaav, g%Domain)
2109 call start_group_pass(cs%pass_ubta_uhbta, g%DoMain)
2111 call do_group_pass(cs%pass_ubta_uhbta, g%Domain)
2113 if (id_clock_pass_post > 0)
call cpu_clock_end(id_clock_pass_post)
2114 if (id_clock_calc_post > 0)
call cpu_clock_begin(id_clock_calc_post)
2119 do j=js,je ;
do i=is-1,ie
2120 accel_layer_u(i,j,k) = us%L_T2_to_m_s2 * (u_accel_bt(i,j) - &
2121 ((pbce(i+1,j,k) - gtot_w(i+1,j)) * e_anom(i+1,j) - &
2122 (pbce(i,j,k) - gtot_e(i,j)) * e_anom(i,j)) * cs%IdxCu(i,j) )
2123 if (abs(accel_layer_u(i,j,k)) < accel_underflow) accel_layer_u(i,j,k) = 0.0
2125 do j=js-1,je ;
do i=is,ie
2126 accel_layer_v(i,j,k) = us%L_T2_to_m_s2 * (v_accel_bt(i,j) - &
2127 ((pbce(i,j+1,k) - gtot_s(i,j+1)) * e_anom(i,j+1) - &
2128 (pbce(i,j,k) - gtot_n(i,j)) * e_anom(i,j)) * cs%IdyCv(i,j) )
2129 if (abs(accel_layer_v(i,j,k)) < accel_underflow) accel_layer_v(i,j,k) = 0.0
2133 if (apply_obcs)
then
2136 if (cs%BT_OBC%apply_u_OBCs)
then ;
do j=js,je ;
do i=is-1,ie
2137 if (obc%segnum_u(i,j) /= obc_none)
then
2138 u_accel_bt(i,j) = (ubt_wtd(i,j) - ubt_first(i,j)) / dt_in_t
2139 do k=1,nz ; accel_layer_u(i,j,k) = us%L_T2_to_m_s2*u_accel_bt(i,j) ;
enddo
2141 enddo ;
enddo ;
endif
2142 if (cs%BT_OBC%apply_v_OBCs)
then ;
do j=js-1,je ;
do i=is,ie
2143 if (obc%segnum_v(i,j) /= obc_none)
then
2144 v_accel_bt(i,j) = (vbt_wtd(i,j) - vbt_first(i,j)) / dt_in_t
2145 do k=1,nz ; accel_layer_v(i,j,k) = us%L_T2_to_m_s2*v_accel_bt(i,j) ;
enddo
2147 enddo ;
enddo ;
endif
2150 if (id_clock_calc_post > 0)
call cpu_clock_end(id_clock_calc_post)
2153 if (query_averaging_enabled(cs%diag))
then
2155 do j=js,je ;
do i=is-1,ie ; cs%ubt_IC(i,j) = ubt_wtd(i,j) ;
enddo ;
enddo
2156 do j=js-1,je ;
do i=is,ie ; cs%vbt_IC(i,j) = vbt_wtd(i,j) ;
enddo ;
enddo
2157 if (use_bt_cont)
then
2158 do j=js,je ;
do i=is-1,ie
2159 cs%uhbt_IC(i,j) = find_uhbt(ubt_wtd(i,j), btcl_u(i,j), us) + uhbt0(i,j)
2161 do j=js-1,je ;
do i=is,ie
2162 cs%vhbt_IC(i,j) = find_vhbt(vbt_wtd(i,j), btcl_v(i,j), us) + vhbt0(i,j)
2165 do j=js,je ;
do i=is-1,ie
2166 cs%uhbt_IC(i,j) = ubt_wtd(i,j) * datu(i,j) + uhbt0(i,j)
2168 do j=js-1,je ;
do i=is,ie
2169 cs%vhbt_IC(i,j) = vbt_wtd(i,j) * datv(i,j) + vhbt0(i,j)
2174 if (cs%id_PFu_bt > 0)
then
2175 do j=js,je ;
do i=is-1,ie
2176 pfu_bt_sum(i,j) = pfu_bt_sum(i,j) * i_sum_wt_accel
2178 call post_data(cs%id_PFu_bt, pfu_bt_sum(isdb:iedb,jsd:jed), cs%diag)
2180 if (cs%id_PFv_bt > 0)
then
2181 do j=js-1,je ;
do i=is,ie
2182 pfv_bt_sum(i,j) = pfv_bt_sum(i,j) * i_sum_wt_accel
2184 call post_data(cs%id_PFv_bt, pfv_bt_sum(isd:ied,jsdb:jedb), cs%diag)
2186 if (cs%id_Coru_bt > 0)
then
2187 do j=js,je ;
do i=is-1,ie
2188 coru_bt_sum(i,j) = coru_bt_sum(i,j) * i_sum_wt_accel
2190 call post_data(cs%id_Coru_bt, coru_bt_sum(isdb:iedb,jsd:jed), cs%diag)
2192 if (cs%id_Corv_bt > 0)
then
2193 do j=js-1,je ;
do i=is,ie
2194 corv_bt_sum(i,j) = corv_bt_sum(i,j) * i_sum_wt_accel
2196 call post_data(cs%id_Corv_bt, corv_bt_sum(isd:ied,jsdb:jedb), cs%diag)
2198 if (cs%id_ubtforce > 0)
call post_data(cs%id_ubtforce, bt_force_u(isdb:iedb,jsd:jed), cs%diag)
2199 if (cs%id_vbtforce > 0)
call post_data(cs%id_vbtforce, bt_force_v(isd:ied,jsdb:jedb), cs%diag)
2200 if (cs%id_uaccel > 0)
call post_data(cs%id_uaccel, u_accel_bt(isdb:iedb,jsd:jed), cs%diag)
2201 if (cs%id_vaccel > 0)
call post_data(cs%id_vaccel, v_accel_bt(isd:ied,jsdb:jedb), cs%diag)
2203 if (cs%id_eta_cor > 0)
call post_data(cs%id_eta_cor, cs%eta_cor, cs%diag)
2204 if (cs%id_eta_bt > 0)
call post_data(cs%id_eta_bt, eta_out, cs%diag)
2205 if (cs%id_gtotn > 0)
call post_data(cs%id_gtotn, gtot_n(isd:ied,jsd:jed), cs%diag)
2206 if (cs%id_gtots > 0)
call post_data(cs%id_gtots, gtot_s(isd:ied,jsd:jed), cs%diag)
2207 if (cs%id_gtote > 0)
call post_data(cs%id_gtote, gtot_e(isd:ied,jsd:jed), cs%diag)
2208 if (cs%id_gtotw > 0)
call post_data(cs%id_gtotw, gtot_w(isd:ied,jsd:jed), cs%diag)
2209 if (cs%id_ubt > 0)
call post_data(cs%id_ubt, ubt_wtd(isdb:iedb,jsd:jed), cs%diag)
2210 if (cs%id_vbt > 0)
call post_data(cs%id_vbt, vbt_wtd(isd:ied,jsdb:jedb), cs%diag)
2211 if (cs%id_ubtav > 0)
call post_data(cs%id_ubtav, cs%ubtav, cs%diag)
2212 if (cs%id_vbtav > 0)
call post_data(cs%id_vbtav, cs%vbtav, cs%diag)
2213 if (cs%id_visc_rem_u > 0)
call post_data(cs%id_visc_rem_u, visc_rem_u, cs%diag)
2214 if (cs%id_visc_rem_v > 0)
call post_data(cs%id_visc_rem_v, visc_rem_v, cs%diag)
2216 if (cs%id_frhatu > 0)
call post_data(cs%id_frhatu, cs%frhatu, cs%diag)
2217 if (cs%id_uhbt > 0)
call post_data(cs%id_uhbt, uhbtav, cs%diag)
2218 if (cs%id_frhatv > 0)
call post_data(cs%id_frhatv, cs%frhatv, cs%diag)
2219 if (cs%id_vhbt > 0)
call post_data(cs%id_vhbt, vhbtav, cs%diag)
2220 if (cs%id_uhbt0 > 0)
call post_data(cs%id_uhbt0, uhbt0(isdb:iedb,jsd:jed), cs%diag)
2221 if (cs%id_vhbt0 > 0)
call post_data(cs%id_vhbt0, vhbt0(isd:ied,jsdb:jedb), cs%diag)
2223 if (cs%id_frhatu1 > 0)
call post_data(cs%id_frhatu1, cs%frhatu1, cs%diag)
2224 if (cs%id_frhatv1 > 0)
call post_data(cs%id_frhatv1, cs%frhatv1, cs%diag)
2226 if (use_bt_cont)
then
2227 if (cs%id_BTC_FA_u_EE > 0)
call post_data(cs%id_BTC_FA_u_EE, bt_cont%FA_u_EE, cs%diag)
2228 if (cs%id_BTC_FA_u_E0 > 0)
call post_data(cs%id_BTC_FA_u_E0, bt_cont%FA_u_E0, cs%diag)
2229 if (cs%id_BTC_FA_u_W0 > 0)
call post_data(cs%id_BTC_FA_u_W0, bt_cont%FA_u_W0, cs%diag)
2230 if (cs%id_BTC_FA_u_WW > 0)
call post_data(cs%id_BTC_FA_u_WW, bt_cont%FA_u_WW, cs%diag)
2231 if (cs%id_BTC_uBT_EE > 0)
call post_data(cs%id_BTC_uBT_EE, bt_cont%uBT_EE, cs%diag)
2232 if (cs%id_BTC_uBT_WW > 0)
call post_data(cs%id_BTC_uBT_WW, bt_cont%uBT_WW, cs%diag)
2233 if (cs%id_BTC_FA_v_NN > 0)
call post_data(cs%id_BTC_FA_v_NN, bt_cont%FA_v_NN, cs%diag)
2234 if (cs%id_BTC_FA_v_N0 > 0)
call post_data(cs%id_BTC_FA_v_N0, bt_cont%FA_v_N0, cs%diag)
2235 if (cs%id_BTC_FA_v_S0 > 0)
call post_data(cs%id_BTC_FA_v_S0, bt_cont%FA_v_S0, cs%diag)
2236 if (cs%id_BTC_FA_v_SS > 0)
call post_data(cs%id_BTC_FA_v_SS, bt_cont%FA_v_SS, cs%diag)
2237 if (cs%id_BTC_vBT_NN > 0)
call post_data(cs%id_BTC_vBT_NN, bt_cont%vBT_NN, cs%diag)
2238 if (cs%id_BTC_vBT_SS > 0)
call post_data(cs%id_BTC_vBT_SS, bt_cont%vBT_SS, cs%diag)
2241 if (cs%id_frhatu1 > 0) cs%frhatu1(:,:,:) = cs%frhatu(:,:,:)
2242 if (cs%id_frhatv1 > 0) cs%frhatv1(:,:,:) = cs%frhatv(:,:,:)
2245 if (g%nonblocking_updates)
then
2246 if (find_etaav)
call complete_group_pass(cs%pass_etaav, g%Domain)
2247 call complete_group_pass(cs%pass_ubta_uhbta, g%Domain)
2250 end subroutine btstep
2254 subroutine set_dtbt(G, GV, US, CS, eta, pbce, BT_cont, gtot_est, SSH_add)
2259 real,
dimension(SZI_(G),SZJ_(G)),
optional,
intent(in) :: eta
2261 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
optional,
intent(in) :: pbce
2267 real,
optional,
intent(in) :: gtot_est
2269 real,
optional,
intent(in) :: ssh_add
2274 real,
dimension(SZI_(G),SZJ_(G)) :: &
2281 real,
dimension(SZIBS_(G),SZJ_(G)) :: &
2284 real,
dimension(SZI_(G),SZJBS_(G)) :: &
2301 logical :: use_bt_cont
2304 character(len=200) :: mesg
2305 integer :: i, j, k, is, ie, js, je, nz
2307 if (.not.
associated(cs))
call mom_error(fatal, &
2308 "set_dtbt: Module MOM_barotropic must be initialized before it is used.")
2309 if (.not.cs%split)
return
2310 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
2311 ms%isdw = g%isd ; ms%iedw = g%ied ; ms%jsdw = g%jsd ; ms%jedw = g%jed
2313 if (.not.(
present(pbce) .or.
present(gtot_est)))
call mom_error(fatal, &
2314 "set_dtbt: Either pbce or gtot_est must be present.")
2316 add_ssh = 0.0 ;
if (
present(ssh_add)) add_ssh = ssh_add
2318 use_bt_cont = .false.
2319 if (
present(bt_cont)) use_bt_cont = (
associated(bt_cont))
2321 if (use_bt_cont)
then
2322 call bt_cont_to_face_areas(bt_cont, datu, datv, g, us, ms, 0, .true.)
2323 elseif (cs%Nonlinear_continuity .and.
present(eta))
then
2324 call find_face_areas(datu, datv, g, gv, us, cs, ms, eta=eta, halo=0)
2326 call find_face_areas(datu, datv, g, gv, us, cs, ms, halo=0, add_max=add_ssh)
2330 if (cs%tides)
call tidal_forcing_sensitivity(g, cs%tides_CSp, det_de)
2331 dgeo_de = 1.0 + max(0.0, det_de + cs%G_extra)
2332 if (
present(pbce))
then
2333 do j=js,je ;
do i=is,ie
2334 gtot_e(i,j) = 0.0 ; gtot_w(i,j) = 0.0
2335 gtot_n(i,j) = 0.0 ; gtot_s(i,j) = 0.0
2337 do k=1,nz ;
do j=js,je ;
do i=is,ie
2338 gtot_e(i,j) = gtot_e(i,j) + pbce(i,j,k) * cs%frhatu(i,j,k)
2339 gtot_w(i,j) = gtot_w(i,j) + pbce(i,j,k) * cs%frhatu(i-1,j,k)
2340 gtot_n(i,j) = gtot_n(i,j) + pbce(i,j,k) * cs%frhatv(i,j,k)
2341 gtot_s(i,j) = gtot_s(i,j) + pbce(i,j,k) * cs%frhatv(i,j-1,k)
2342 enddo ;
enddo ;
enddo
2344 do j=js,je ;
do i=is,ie
2345 gtot_e(i,j) = gtot_est * gv%H_to_Z ; gtot_w(i,j) = gtot_est * gv%H_to_Z
2346 gtot_n(i,j) = gtot_est * gv%H_to_Z ; gtot_s(i,j) = gtot_est * gv%H_to_Z
2350 min_max_dt2 = 1.0e38*us%s_to_T**2
2351 do j=js,je ;
do i=is,ie
2354 idt_max2 = 0.5 * (1.0 + 2.0*cs%bebt) * (us%L_to_m**2*g%IareaT(i,j) * &
2355 ((gtot_e(i,j)*datu(i,j)*us%L_to_m*g%IdxCu(i,j) + gtot_w(i,j)*datu(i-1,j)*us%L_to_m*g%IdxCu(i-1,j)) + &
2356 (gtot_n(i,j)*datv(i,j)*us%L_to_m*g%IdyCv(i,j) + gtot_s(i,j)*datv(i,j-1)*us%L_to_m*g%IdyCv(i,j-1))) + &
2357 ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
2358 (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2)))
2359 if (idt_max2 * min_max_dt2 > 1.0) min_max_dt2 = 1.0 / idt_max2
2361 dtbt_max = sqrt(min_max_dt2 / dgeo_de)
2362 if (id_clock_sync > 0)
call cpu_clock_begin(id_clock_sync)
2363 call min_across_pes(dtbt_max)
2364 if (id_clock_sync > 0)
call cpu_clock_end(id_clock_sync)
2366 cs%dtbt = cs%dtbt_fraction * us%T_to_s * dtbt_max
2367 cs%dtbt_max = us%T_to_s * dtbt_max
2368 end subroutine set_dtbt
2373 subroutine apply_velocity_obcs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans, &
2374 eta, ubt_old, vbt_old, BT_OBC, &
2375 G, MS, US, halo, dtbt, bebt, use_BT_cont, Datu, Datv, &
2376 BTCL_u, BTCL_v, uhbt0, vhbt0)
2381 real,
dimension(SZIBW_(MS),SZJW_(MS)),
intent(inout) :: ubt
2382 real,
dimension(SZIBW_(MS),SZJW_(MS)),
intent(inout) :: uhbt
2384 real,
dimension(SZIBW_(MS),SZJW_(MS)),
intent(inout) :: ubt_trans
2386 real,
dimension(SZIW_(MS),SZJBW_(MS)),
intent(inout) :: vbt
2387 real,
dimension(SZIW_(MS),SZJBW_(MS)),
intent(inout) :: vhbt
2389 real,
dimension(SZIW_(MS),SZJBW_(MS)),
intent(inout) :: vbt_trans
2391 real,
dimension(SZIW_(MS),SZJW_(MS)),
intent(in) :: eta
2393 real,
dimension(SZIBW_(MS),SZJW_(MS)),
intent(in) :: ubt_old
2395 real,
dimension(SZIW_(MS),SZJBW_(MS)),
intent(in) :: vbt_old
2401 integer,
intent(in) :: halo
2402 real,
intent(in) :: dtbt
2403 real,
intent(in) :: bebt
2405 logical,
intent(in) :: use_BT_cont
2407 real,
dimension(SZIBW_(MS),SZJW_(MS)),
intent(in) :: Datu
2409 real,
dimension(SZIW_(MS),SZJBW_(MS)),
intent(in) :: Datv
2417 real,
dimension(SZIBW_(MS),SZJW_(MS)),
intent(in) :: uhbt0
2421 real,
dimension(SZIW_(MS),SZJBW_(MS)),
intent(in) :: vhbt0
2436 real :: cff, Cx, Cy, tau
2437 real :: dhdt, dhdx, dhdy
2438 integer :: i, j, is, ie, js, je
2439 real,
dimension(SZIB_(G),SZJB_(G)) :: grad
2440 real,
parameter :: eps = 1.0e-20
2441 real :: rx_max, ry_max
2442 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
2443 rx_max = obc%rx_max ; ry_max = obc%rx_max
2445 if (bt_obc%apply_u_OBCs)
then
2446 do j=js,je ;
do i=is-1,ie ;
if (obc%segnum_u(i,j) /= obc_none)
then
2447 if (obc%segment(obc%segnum_u(i,j))%specified)
then
2448 uhbt(i,j) = bt_obc%uhbt(i,j)
2449 ubt(i,j) = bt_obc%ubt_outer(i,j)
2450 vel_trans = ubt(i,j)
2451 elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e)
then
2452 if (obc%segment(obc%segnum_u(i,j))%Flather)
then
2453 cfl = dtbt * bt_obc%Cg_u(i,j) * us%L_to_m*g%IdxCu(i,j)
2454 u_inlet = cfl*ubt_old(i-1,j) + (1.0-cfl)*ubt_old(i,j)
2455 h_in = eta(i,j) + (0.5-cfl)*(eta(i,j)-eta(i-1,j))
2456 h_u = bt_obc%H_u(i,j)
2458 ubt(i,j) = 0.5*((u_inlet + bt_obc%ubt_outer(i,j)) + &
2459 (bt_obc%Cg_u(i,j)/h_u) * (h_in-bt_obc%eta_outer_u(i,j)))
2460 vel_trans = (1.0-bebt)*vel_prev + bebt*ubt(i,j)
2461 elseif (obc%segment(obc%segnum_u(i,j))%gradient)
then
2462 ubt(i,j) = ubt(i-1,j)
2463 vel_trans = ubt(i,j)
2465 elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w)
then
2466 if (obc%segment(obc%segnum_u(i,j))%Flather)
then
2467 cfl = dtbt * bt_obc%Cg_u(i,j) * us%L_to_m*g%IdxCu(i,j)
2468 u_inlet = cfl*ubt_old(i+1,j) + (1.0-cfl)*ubt_old(i,j)
2469 h_in = eta(i+1,j) + (0.5-cfl)*(eta(i+1,j)-eta(i+2,j))
2471 h_u = bt_obc%H_u(i,j)
2473 ubt(i,j) = 0.5*((u_inlet + bt_obc%ubt_outer(i,j)) + &
2474 (bt_obc%Cg_u(i,j)/h_u) * (bt_obc%eta_outer_u(i,j)-h_in))
2476 vel_trans = (1.0-bebt)*vel_prev + bebt*ubt(i,j)
2477 elseif (obc%segment(obc%segnum_u(i,j))%gradient)
then
2478 ubt(i,j) = ubt(i+1,j)
2479 vel_trans = ubt(i,j)
2483 if (.not. obc%segment(obc%segnum_u(i,j))%specified)
then
2484 if (use_bt_cont)
then
2485 uhbt(i,j) = find_uhbt(vel_trans, btcl_u(i,j), us) + uhbt0(i,j)
2487 uhbt(i,j) = datu(i,j)*vel_trans + uhbt0(i,j)
2491 ubt_trans(i,j) = vel_trans
2492 endif ;
enddo ;
enddo
2495 if (bt_obc%apply_v_OBCs)
then
2496 do j=js-1,je ;
do i=is,ie ;
if (obc%segnum_v(i,j) /= obc_none)
then
2497 if (obc%segment(obc%segnum_v(i,j))%specified)
then
2498 vhbt(i,j) = bt_obc%vhbt(i,j)
2499 vbt(i,j) = bt_obc%vbt_outer(i,j)
2500 vel_trans = vbt(i,j)
2501 elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n)
then
2502 if (obc%segment(obc%segnum_v(i,j))%Flather)
then
2503 cfl = dtbt * bt_obc%Cg_v(i,j) * us%L_to_m*g%IdyCv(i,j)
2504 v_inlet = cfl*vbt_old(i,j-1) + (1.0-cfl)*vbt_old(i,j)
2505 h_in = eta(i,j) + (0.5-cfl)*(eta(i,j)-eta(i,j-1))
2507 h_v = bt_obc%H_v(i,j)
2509 vbt(i,j) = 0.5*((v_inlet + bt_obc%vbt_outer(i,j)) + &
2510 (bt_obc%Cg_v(i,j)/h_v) * (h_in-bt_obc%eta_outer_v(i,j)))
2512 vel_trans = (1.0-bebt)*vel_prev + bebt*vbt(i,j)
2513 elseif (obc%segment(obc%segnum_v(i,j))%gradient)
then
2514 vbt(i,j) = vbt(i,j-1)
2515 vel_trans = vbt(i,j)
2517 elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s)
then
2518 if (obc%segment(obc%segnum_v(i,j))%Flather)
then
2519 cfl = dtbt * bt_obc%Cg_v(i,j) * us%L_to_m*g%IdyCv(i,j)
2520 v_inlet = cfl*vbt_old(i,j+1) + (1.0-cfl)*vbt_old(i,j)
2521 h_in = eta(i,j+1) + (0.5-cfl)*(eta(i,j+1)-eta(i,j+2))
2523 h_v = bt_obc%H_v(i,j)
2525 vbt(i,j) = 0.5*((v_inlet + bt_obc%vbt_outer(i,j)) + &
2526 (bt_obc%Cg_v(i,j)/h_v) * (bt_obc%eta_outer_v(i,j)-h_in))
2528 vel_trans = (1.0-bebt)*vel_prev + bebt*vbt(i,j)
2529 elseif (obc%segment(obc%segnum_v(i,j))%gradient)
then
2530 vbt(i,j) = vbt(i,j+1)
2531 vel_trans = vbt(i,j)
2535 if (.not. obc%segment(obc%segnum_v(i,j))%specified)
then
2536 if (use_bt_cont)
then
2537 vhbt(i,j) = find_vhbt(vel_trans, btcl_v(i,j), us) + vhbt0(i,j)
2539 vhbt(i,j) = vel_trans*datv(i,j) + vhbt0(i,j)
2543 vbt_trans(i,j) = vel_trans
2544 endif ;
enddo ;
enddo
2547 end subroutine apply_velocity_obcs
2551 subroutine set_up_bt_obc(OBC, eta, BT_OBC, BT_Domain, G, GV, US, MS, halo, use_BT_cont, Datu, Datv, BTCL_u, BTCL_v)
2555 real,
dimension(SZIW_(MS),SZJW_(MS)),
intent(in) :: eta
2564 integer,
intent(in) :: halo
2565 logical,
intent(in) :: use_BT_cont
2567 real,
dimension(SZIBW_(MS),SZJW_(MS)),
intent(in) :: Datu
2569 real,
dimension(SZIW_(MS),SZJBW_(MS)),
intent(in) :: Datv
2579 integer :: i, j, k, is, ie, js, je, n, nz, Isq, Ieq, Jsq, Jeq
2580 integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
2581 integer :: isdw, iedw, jsdw, jedw
2586 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
2587 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
2588 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
2589 isdw = ms%isdw ; iedw = ms%iedw ; jsdw = ms%jsdw ; jedw = ms%jedw
2592 if ((isdw < isd) .or. (jsdw < jsd))
then
2593 call mom_error(fatal,
"set_up_BT_OBC: Open boundary conditions are not "//&
2594 "yet fully implemented with wide barotropic halos.")
2597 if (.not. bt_obc%is_alloced)
then
2598 allocate(bt_obc%Cg_u(isdw-1:iedw,jsdw:jedw)) ; bt_obc%Cg_u(:,:) = 0.0
2599 allocate(bt_obc%H_u(isdw-1:iedw,jsdw:jedw)) ; bt_obc%H_u(:,:) = 0.0
2600 allocate(bt_obc%uhbt(isdw-1:iedw,jsdw:jedw)) ; bt_obc%uhbt(:,:) = 0.0
2601 allocate(bt_obc%ubt_outer(isdw-1:iedw,jsdw:jedw)) ; bt_obc%ubt_outer(:,:) = 0.0
2602 allocate(bt_obc%eta_outer_u(isdw-1:iedw,jsdw:jedw)) ; bt_obc%eta_outer_u(:,:) = 0.0
2604 allocate(bt_obc%Cg_v(isdw:iedw,jsdw-1:jedw)) ; bt_obc%Cg_v(:,:) = 0.0
2605 allocate(bt_obc%H_v(isdw:iedw,jsdw-1:jedw)) ; bt_obc%H_v(:,:) = 0.0
2606 allocate(bt_obc%vhbt(isdw:iedw,jsdw-1:jedw)) ; bt_obc%vhbt(:,:) = 0.0
2607 allocate(bt_obc%vbt_outer(isdw:iedw,jsdw-1:jedw)) ; bt_obc%vbt_outer(:,:) = 0.0
2608 allocate(bt_obc%eta_outer_v(isdw:iedw,jsdw-1:jedw)) ; bt_obc%eta_outer_v(:,:)=0.0
2609 bt_obc%is_alloced = .true.
2610 call create_group_pass(bt_obc%pass_uv, bt_obc%ubt_outer, bt_obc%vbt_outer, bt_domain)
2612 call create_group_pass(bt_obc%pass_eta_outer, bt_obc%eta_outer_u, bt_obc%eta_outer_v, bt_domain,to_all+scalar_pair)
2613 call create_group_pass(bt_obc%pass_h, bt_obc%H_u, bt_obc%H_v, bt_domain,to_all+scalar_pair)
2614 call create_group_pass(bt_obc%pass_cg, bt_obc%Cg_u, bt_obc%Cg_v, bt_domain,to_all+scalar_pair)
2617 if (bt_obc%apply_u_OBCs)
then
2618 if (obc%specified_u_BCs_exist_globally)
then
2619 do n = 1, obc%number_of_segments
2620 segment => obc%segment(n)
2621 if (segment%is_E_or_W .and. segment%specified)
then
2622 do j=segment%HI%jsd,segment%HI%jed ;
do i=segment%HI%IsdB,segment%HI%IedB
2623 bt_obc%uhbt(i,j) = 0.
2625 do k=1,nz ;
do j=segment%HI%jsd,segment%HI%jed ;
do i=segment%HI%IsdB,segment%HI%IedB
2626 bt_obc%uhbt(i,j) = bt_obc%uhbt(i,j) + us%T_to_s*us%m_to_L**2*segment%normal_trans(i,j,k)
2627 enddo ;
enddo ;
enddo
2631 do j=js,je ;
do i=is-1,ie ;
if (obc%segnum_u(i,j) /= obc_none)
then
2633 if (obc%segment(obc%segnum_u(i,j))%specified)
then
2634 if (use_bt_cont)
then
2635 bt_obc%ubt_outer(i,j) = uhbt_to_ubt(bt_obc%uhbt(i,j), btcl_u(i,j), us)
2637 if (datu(i,j) > 0.0) bt_obc%ubt_outer(i,j) = bt_obc%uhbt(i,j) / datu(i,j)
2640 if (gv%Boussinesq)
then
2641 if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e)
then
2642 bt_obc%H_u(i,j) = g%bathyT(i,j)*gv%Z_to_H + eta(i,j)
2643 elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w)
then
2644 bt_obc%H_u(i,j) = g%bathyT(i+1,j)*gv%Z_to_H + eta(i+1,j)
2647 if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e)
then
2648 bt_obc%H_u(i,j) = eta(i,j)
2649 elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w)
then
2650 bt_obc%H_u(i,j) = eta(i+1,j)
2653 bt_obc%Cg_u(i,j) = sqrt(gv%g_prime(1) * gv%H_to_Z*bt_obc%H_u(i,j))
2655 endif ;
enddo ;
enddo
2656 if (obc%Flather_u_BCs_exist_globally)
then
2657 do n = 1, obc%number_of_segments
2658 segment => obc%segment(n)
2659 if (segment%is_E_or_W .and. segment%Flather)
then
2660 do j=segment%HI%jsd,segment%HI%jed ;
do i=segment%HI%IsdB,segment%HI%IedB
2661 bt_obc%ubt_outer(i,j) = us%m_s_to_L_T*segment%normal_vel_bt(i,j)
2662 bt_obc%eta_outer_u(i,j) = segment%eta(i,j)
2669 if (bt_obc%apply_v_OBCs)
then
2670 if (obc%specified_v_BCs_exist_globally)
then
2671 do n = 1, obc%number_of_segments
2672 segment => obc%segment(n)
2673 if (segment%is_N_or_S .and. segment%specified)
then
2674 do j=segment%HI%JsdB,segment%HI%JedB ;
do i=segment%HI%isd,segment%HI%ied
2675 bt_obc%vhbt(i,j) = 0.
2677 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB ;
do i=segment%HI%isd,segment%HI%ied
2678 bt_obc%vhbt(i,j) = bt_obc%vhbt(i,j) + us%T_to_s*us%m_to_L**2*segment%normal_trans(i,j,k)
2679 enddo ;
enddo ;
enddo
2683 do j=js-1,je ;
do i=is,ie ;
if (obc%segnum_v(i,j) /= obc_none)
then
2685 if (obc%segment(obc%segnum_v(i,j))%specified)
then
2686 if (use_bt_cont)
then
2687 bt_obc%vbt_outer(i,j) = vhbt_to_vbt(bt_obc%vhbt(i,j), btcl_v(i,j), us)
2689 if (datv(i,j) > 0.0) bt_obc%vbt_outer(i,j) = bt_obc%vhbt(i,j) / datv(i,j)
2692 if (gv%Boussinesq)
then
2693 if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n)
then
2694 bt_obc%H_v(i,j) = g%bathyT(i,j)*gv%Z_to_H + eta(i,j)
2695 elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s)
then
2696 bt_obc%H_v(i,j) = g%bathyT(i,j+1)*gv%Z_to_H + eta(i,j+1)
2699 if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n)
then
2700 bt_obc%H_v(i,j) = eta(i,j)
2701 elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s)
then
2702 bt_obc%H_v(i,j) = eta(i,j+1)
2705 bt_obc%Cg_v(i,j) = sqrt(gv%g_prime(1) * gv%H_to_Z*bt_obc%H_v(i,j))
2707 endif ;
enddo ;
enddo
2708 if (obc%Flather_v_BCs_exist_globally)
then
2709 do n = 1, obc%number_of_segments
2710 segment => obc%segment(n)
2711 if (segment%is_N_or_S .and. segment%Flather)
then
2712 do j=segment%HI%JsdB,segment%HI%JedB ;
do i=segment%HI%isd,segment%HI%ied
2713 bt_obc%vbt_outer(i,j) = us%m_s_to_L_T*segment%normal_vel_bt(i,j)
2714 bt_obc%eta_outer_v(i,j) = segment%eta(i,j)
2721 call do_group_pass(bt_obc%pass_uv, bt_domain)
2722 call do_group_pass(bt_obc%pass_uhvh, bt_domain)
2723 call do_group_pass(bt_obc%pass_eta_outer, bt_domain)
2724 call do_group_pass(bt_obc%pass_h, bt_domain)
2725 call do_group_pass(bt_obc%pass_cg, bt_domain)
2727 end subroutine set_up_bt_obc
2730 subroutine destroy_bt_obc(BT_OBC)
2735 if (bt_obc%is_alloced)
then
2736 deallocate(bt_obc%Cg_u)
2737 deallocate(bt_obc%H_u)
2738 deallocate(bt_obc%uhbt)
2739 deallocate(bt_obc%ubt_outer)
2740 deallocate(bt_obc%eta_outer_u)
2742 deallocate(bt_obc%Cg_v)
2743 deallocate(bt_obc%H_v)
2744 deallocate(bt_obc%vhbt)
2745 deallocate(bt_obc%vbt_outer)
2746 deallocate(bt_obc%eta_outer_v)
2747 bt_obc%is_alloced = .false.
2749 end subroutine destroy_bt_obc
2756 subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC)
2759 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
2763 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
2764 optional,
intent(in) :: h_u
2765 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
2766 optional,
intent(in) :: h_v
2767 logical,
optional,
intent(in) :: may_use_default
2776 real :: hatutot(szib_(g))
2777 real :: hatvtot(szi_(g))
2778 real :: ihatutot(szib_(g))
2779 real :: ihatvtot(szi_(g))
2787 real :: e_u(szib_(g),szk_(g)+1)
2788 real :: e_v(szi_(g),szk_(g)+1)
2789 real :: d_shallow_u(szi_(g))
2790 real :: d_shallow_v(szib_(g))
2794 logical :: use_default, test_dflt, apply_obcs
2795 integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz, i, j, k
2796 integer :: iss, ies, n
2800 if (.not.
associated(cs))
call mom_error(fatal, &
2801 "btcalc: Module MOM_barotropic must be initialized before it is used.")
2802 if (.not.cs%split)
return
2804 use_default = .false.
2805 test_dflt = .false. ;
if (
present(may_use_default)) test_dflt = may_use_default
2808 if (.not.((
present(h_u) .and.
present(h_v)) .or. &
2809 (cs%hvel_scheme == harmonic) .or. (cs%hvel_scheme == hybrid) .or.&
2810 (cs%hvel_scheme == arithmetic))) use_default = .true.
2812 if (.not.((
present(h_u) .and.
present(h_v)) .or. &
2813 (cs%hvel_scheme == harmonic) .or. (cs%hvel_scheme == hybrid) .or.&
2814 (cs%hvel_scheme == arithmetic)))
call mom_error(fatal, &
2815 "btcalc: Inconsistent settings of optional arguments and hvel_scheme.")
2818 apply_obcs = .false.
2819 if (
present(obc))
then ;
if (
associated(obc))
then ;
if (obc%OBC_pe)
then
2822 apply_obcs = (obc%number_of_segments > 0)
2823 endif ;
endif ;
endif
2825 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
2826 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
2827 h_neglect = gv%H_subroundoff
2835 if (
present(h_u))
then
2836 do i=is-1,ie ; hatutot(i) = h_u(i,j,1) ;
enddo
2837 do k=2,nz ;
do i=is-1,ie
2838 hatutot(i) = hatutot(i) + h_u(i,j,k)
2840 do i=is-1,ie ; ihatutot(i) = g%mask2dCu(i,j) / (hatutot(i) + h_neglect) ;
enddo
2841 do k=1,nz ;
do i=is-1,ie
2842 cs%frhatu(i,j,k) = h_u(i,j,k) * ihatutot(i)
2845 if (cs%hvel_scheme == arithmetic)
then
2847 cs%frhatu(i,j,1) = 0.5 * (h(i+1,j,1) + h(i,j,1))
2848 hatutot(i) = cs%frhatu(i,j,1)
2850 do k=2,nz ;
do i=is-1,ie
2851 cs%frhatu(i,j,k) = 0.5 * (h(i+1,j,k) + h(i,j,k))
2852 hatutot(i) = hatutot(i) + cs%frhatu(i,j,k)
2854 elseif (cs%hvel_scheme == hybrid .or. use_default)
then
2856 e_u(i,nz+1) = -0.5 * gv%Z_to_H * (g%bathyT(i+1,j) + g%bathyT(i,j))
2857 d_shallow_u(i) = -gv%Z_to_H * min(g%bathyT(i+1,j), g%bathyT(i,j))
2860 do k=nz,1,-1 ;
do i=is-1,ie
2861 e_u(i,k) = e_u(i,k+1) + 0.5 * (h(i+1,j,k) + h(i,j,k))
2862 h_arith = 0.5 * (h(i+1,j,k) + h(i,j,k))
2863 if (e_u(i,k+1) >= d_shallow_u(i))
then
2864 cs%frhatu(i,j,k) = h_arith
2866 h_harm = (h(i+1,j,k) * h(i,j,k)) / (h_arith + h_neglect)
2867 if (e_u(i,k) <= d_shallow_u(i))
then
2868 cs%frhatu(i,j,k) = h_harm
2870 wt_arith = (e_u(i,k) - d_shallow_u(i)) / (h_arith + h_neglect)
2871 cs%frhatu(i,j,k) = wt_arith*h_arith + (1.0-wt_arith)*h_harm
2874 hatutot(i) = hatutot(i) + cs%frhatu(i,j,k)
2876 elseif (cs%hvel_scheme == harmonic)
then
2878 cs%frhatu(i,j,1) = 2.0*(h(i+1,j,1) * h(i,j,1)) / &
2879 ((h(i+1,j,1) + h(i,j,1)) + h_neglect)
2880 hatutot(i) = cs%frhatu(i,j,1)
2882 do k=2,nz ;
do i=is-1,ie
2883 cs%frhatu(i,j,k) = 2.0*(h(i+1,j,k) * h(i,j,k)) / &
2884 ((h(i+1,j,k) + h(i,j,k)) + h_neglect)
2885 hatutot(i) = hatutot(i) + cs%frhatu(i,j,k)
2888 do i=is-1,ie ; ihatutot(i) = g%mask2dCu(i,j) / (hatutot(i) + h_neglect) ;
enddo
2889 do k=1,nz ;
do i=is-1,ie
2890 cs%frhatu(i,j,k) = cs%frhatu(i,j,k) * ihatutot(i)
2898 if (
present(h_v))
then
2899 do i=is,ie ; hatvtot(i) = h_v(i,j,1) ;
enddo
2900 do k=2,nz ;
do i=is,ie
2901 hatvtot(i) = hatvtot(i) + h_v(i,j,k)
2903 do i=is,ie ; ihatvtot(i) = g%mask2dCv(i,j) / (hatvtot(i) + h_neglect) ;
enddo
2904 do k=1,nz ;
do i=is,ie
2905 cs%frhatv(i,j,k) = h_v(i,j,k) * ihatvtot(i)
2908 if (cs%hvel_scheme == arithmetic)
then
2910 cs%frhatv(i,j,1) = 0.5 * (h(i,j+1,1) + h(i,j,1))
2911 hatvtot(i) = cs%frhatv(i,j,1)
2913 do k=2,nz ;
do i=is,ie
2914 cs%frhatv(i,j,k) = 0.5 * (h(i,j+1,k) + h(i,j,k))
2915 hatvtot(i) = hatvtot(i) + cs%frhatv(i,j,k)
2917 elseif (cs%hvel_scheme == hybrid .or. use_default)
then
2919 e_v(i,nz+1) = -0.5 * gv%Z_to_H * (g%bathyT(i,j+1) + g%bathyT(i,j))
2920 d_shallow_v(i) = -gv%Z_to_H * min(g%bathyT(i,j+1), g%bathyT(i,j))
2923 do k=nz,1,-1 ;
do i=is,ie
2924 e_v(i,k) = e_v(i,k+1) + 0.5 * (h(i,j+1,k) + h(i,j,k))
2925 h_arith = 0.5 * (h(i,j+1,k) + h(i,j,k))
2926 if (e_v(i,k+1) >= d_shallow_v(i))
then
2927 cs%frhatv(i,j,k) = h_arith
2929 h_harm = (h(i,j+1,k) * h(i,j,k)) / (h_arith + h_neglect)
2930 if (e_v(i,k) <= d_shallow_v(i))
then
2931 cs%frhatv(i,j,k) = h_harm
2933 wt_arith = (e_v(i,k) - d_shallow_v(i)) / (h_arith + h_neglect)
2934 cs%frhatv(i,j,k) = wt_arith*h_arith + (1.0-wt_arith)*h_harm
2937 hatvtot(i) = hatvtot(i) + cs%frhatv(i,j,k)
2939 elseif (cs%hvel_scheme == harmonic)
then
2941 cs%frhatv(i,j,1) = 2.0*(h(i,j+1,1) * h(i,j,1)) / &
2942 ((h(i,j+1,1) + h(i,j,1)) + h_neglect)
2943 hatvtot(i) = cs%frhatv(i,j,1)
2945 do k=2,nz ;
do i=is,ie
2946 cs%frhatv(i,j,k) = 2.0*(h(i,j+1,k) * h(i,j,k)) / &
2947 ((h(i,j+1,k) + h(i,j,k)) + h_neglect)
2948 hatvtot(i) = hatvtot(i) + cs%frhatv(i,j,k)
2951 do i=is,ie ; ihatvtot(i) = g%mask2dCv(i,j) / (hatvtot(i) + h_neglect) ;
enddo
2952 do k=1,nz ;
do i=is,ie
2953 cs%frhatv(i,j,k) = cs%frhatv(i,j,k) * ihatvtot(i)
2958 if (apply_obcs)
then ;
do n=1,obc%number_of_segments
2959 if (.not. obc%segment(n)%on_pe) cycle
2960 if (obc%segment(n)%direction == obc_direction_n)
then
2961 j = obc%segment(n)%HI%JsdB
2962 if ((j >= js-1) .and. (j <= je))
then
2963 iss = max(is,obc%segment(n)%HI%isd) ; ies = min(ie,obc%segment(n)%HI%ied)
2964 do i=iss,ies ; hatvtot(i) = h(i,j,1) ;
enddo
2965 do k=2,nz ;
do i=iss,ies
2966 hatvtot(i) = hatvtot(i) + h(i,j,k)
2969 ihatvtot(i) = g%mask2dCv(i,j) / (hatvtot(i) + h_neglect)
2971 do k=1,nz ;
do i=iss,ies
2972 cs%frhatv(i,j,k) = h(i,j,k) * ihatvtot(i)
2975 elseif (obc%segment(n)%direction == obc_direction_s)
then
2976 j = obc%segment(n)%HI%JsdB
2977 if ((j >= js-1) .and. (j <= je))
then
2978 iss = max(is,obc%segment(n)%HI%isd) ; ies = min(ie,obc%segment(n)%HI%ied)
2979 do i=iss,ies ; hatvtot(i) = h(i,j+1,1) ;
enddo
2980 do k=2,nz ;
do i=iss,ies
2981 hatvtot(i) = hatvtot(i) + h(i,j+1,k)
2984 ihatvtot(i) = g%mask2dCv(i,j) / (hatvtot(i) + h_neglect)
2986 do k=1,nz ;
do i=iss,ies
2987 cs%frhatv(i,j,k) = h(i,j+1,k) * ihatvtot(i)
2990 elseif (obc%segment(n)%direction == obc_direction_e)
then
2991 i = obc%segment(n)%HI%IsdB
2992 if ((i >= is-1) .and. (i <= ie))
then
2993 do j = max(js,obc%segment(n)%HI%jsd), min(je,obc%segment(n)%HI%jed)
2995 do k=2,nz ; htot = htot + h(i,j,k) ;
enddo
2996 ihtot = g%mask2dCu(i,j) / (htot + h_neglect)
2997 do k=1,nz ; cs%frhatu(i,j,k) = h(i,j,k) * ihtot ;
enddo
3000 elseif (obc%segment(n)%direction == obc_direction_w)
then
3001 i = obc%segment(n)%HI%IsdB
3002 if ((i >= is-1) .and. (i <= ie))
then
3003 do j = max(js,obc%segment(n)%HI%jsd), min(je,obc%segment(n)%HI%jed)
3005 do k=2,nz ; htot = htot + h(i+1,j,k) ;
enddo
3006 ihtot = g%mask2dCu(i,j) / (htot + h_neglect)
3007 do k=1,nz ; cs%frhatu(i,j,k) = h(i+1,j,k) * ihtot ;
enddo
3011 call mom_error(fatal,
"btcalc encountered and OBC segment of indeterminate direction.")
3016 call uvchksum(
"btcalc frhat[uv]", cs%frhatu, cs%frhatv, g%HI, 0, .true., .true.)
3017 if (
present(h_u) .and.
present(h_v)) &
3018 call uvchksum(
"btcalc h_[uv]", h_u, h_v, g%HI, 0, .true., .true., scale=gv%H_to_m)
3019 call hchksum(h,
"btcalc h",g%HI, haloshift=1, scale=gv%H_to_m)
3022 end subroutine btcalc
3025 function find_uhbt(u, BTC, US)
result(uhbt)
3026 real,
intent(in) :: u
3036 elseif (u < btc%uBT_EE)
then
3037 uhbt = (u - btc%uBT_EE) * btc%FA_u_EE + btc%uh_EE
3038 elseif (u < 0.0)
then
3039 uhbt = u * (btc%FA_u_E0 + btc%uh_crvE * u**2)
3040 elseif (u <= btc%uBT_WW)
then
3041 uhbt = u * (btc%FA_u_W0 + btc%uh_crvW * u**2)
3043 uhbt = (u - btc%uBT_WW) * btc%FA_u_WW + btc%uh_WW
3046 end function find_uhbt
3050 function uhbt_to_ubt(uhbt, BTC, US, guess)
result(ubt)
3051 real,
intent(in) :: uhbt
3057 real,
optional,
intent(in) :: guess
3062 real :: ubt_min, ubt_max, uhbt_err, derr_du
3063 real :: uherr_min, uherr_max
3064 real,
parameter :: tol = 1.0e-10
3067 real,
parameter :: vs1 = 1.25
3068 real,
parameter :: vs2 = 2.0
3070 integer :: itt, max_itt = 20
3073 if (uhbt == 0.0)
then
3075 elseif (uhbt < btc%uh_EE)
then
3076 ubt = btc%uBT_EE + (uhbt - btc%uh_EE) / btc%FA_u_EE
3077 elseif (uhbt < 0.0)
then
3080 ubt_min = btc%uBT_EE ; uherr_min = btc%uh_EE - uhbt
3081 ubt_max = 0.0 ; uherr_max = -uhbt
3083 ubt = btc%uBT_EE * (uhbt / btc%uh_EE)
3085 uhbt_err = ubt * (btc%FA_u_E0 + btc%uh_crvE * ubt**2) - uhbt
3087 if (abs(uhbt_err) < tol*abs(uhbt))
exit
3088 if (uhbt_err > 0.0)
then ; ubt_max = ubt ; uherr_max = uhbt_err ;
endif
3089 if (uhbt_err < 0.0)
then ; ubt_min = ubt ; uherr_min = uhbt_err ;
endif
3091 derr_du = btc%FA_u_E0 + 3.0 * btc%uh_crvE * ubt**2
3092 if ((uhbt_err >= derr_du*(ubt - ubt_min)) .or. &
3093 (-uhbt_err >= derr_du*(ubt_max - ubt)) .or. (derr_du <= 0.0))
then
3095 ubt = ubt_max + (ubt_min-ubt_max) * (uherr_max / (uherr_max-uherr_min))
3097 ubt = ubt - uhbt_err / derr_du
3098 if (abs(uhbt_err) < (0.01*tol)*abs(ubt_min*derr_du))
exit
3101 elseif (uhbt <= btc%uh_WW)
then
3103 ubt_min = 0.0 ; uherr_min = -uhbt
3104 ubt_max = btc%uBT_WW ; uherr_max = btc%uh_WW - uhbt
3106 ubt = btc%uBT_WW * (uhbt / btc%uh_WW)
3108 uhbt_err = ubt * (btc%FA_u_W0 + btc%uh_crvW * ubt**2) - uhbt
3110 if (abs(uhbt_err) < tol*abs(uhbt))
exit
3111 if (uhbt_err > 0.0)
then ; ubt_max = ubt ; uherr_max = uhbt_err ;
endif
3112 if (uhbt_err < 0.0)
then ; ubt_min = ubt ; uherr_min = uhbt_err ;
endif
3114 derr_du = btc%FA_u_W0 + 3.0 * btc%uh_crvW * ubt**2
3115 if ((uhbt_err >= derr_du*(ubt - ubt_min)) .or. &
3116 (-uhbt_err >= derr_du*(ubt_max - ubt)) .or. (derr_du <= 0.0))
then
3118 ubt = ubt_min + (ubt_max-ubt_min) * (-uherr_min / (uherr_max-uherr_min))
3120 ubt = ubt - uhbt_err / derr_du
3121 if (abs(uhbt_err) < (0.01*tol)*(ubt_max*derr_du))
exit
3125 ubt = btc%uBT_WW + (uhbt - btc%uh_WW) / btc%FA_u_WW
3128 if (
present(guess))
then
3129 dvel = abs(ubt) - vs1*abs(guess)
3130 if (dvel > 0.0)
then
3131 if (dvel < 40.0 * (abs(guess)*(vs2-vs1)) )
then
3132 vsr = vs2 - (vs2-vs1) * exp(-dvel / (abs(guess)*(vs2-vs1)))
3136 ubt = sign(vsr * guess, ubt)
3140 end function uhbt_to_ubt
3143 function find_vhbt(v, BTC, US)
result(vhbt)
3144 real,
intent(in) :: v
3153 elseif (v < btc%vBT_NN)
then
3154 vhbt = (v - btc%vBT_NN) * btc%FA_v_NN + btc%vh_NN
3155 elseif (v < 0.0)
then
3156 vhbt = v * (btc%FA_v_N0 + btc%vh_crvN * v**2)
3157 elseif (v <= btc%vBT_SS)
then
3158 vhbt = v * (btc%FA_v_S0 + btc%vh_crvS * v**2)
3160 vhbt = (v - btc%vBT_SS) * btc%FA_v_SS + btc%vh_SS
3163 end function find_vhbt
3167 function vhbt_to_vbt(vhbt, BTC, US, guess)
result(vbt)
3168 real,
intent(in) :: vhbt
3174 real,
optional,
intent(in) :: guess
3179 real :: vbt_min, vbt_max, vhbt_err, derr_dv
3180 real :: vherr_min, vherr_max
3181 real,
parameter :: tol = 1.0e-10
3184 real,
parameter :: vs1 = 1.25
3185 real,
parameter :: vs2 = 2.0
3187 integer :: itt, max_itt = 20
3190 if (vhbt == 0.0)
then
3192 elseif (vhbt < btc%vh_NN)
then
3193 vbt = btc%vBT_NN + (vhbt - btc%vh_NN) / btc%FA_v_NN
3194 elseif (vhbt < 0.0)
then
3197 vbt_min = btc%vBT_NN ; vherr_min = btc%vh_NN - vhbt
3198 vbt_max = 0.0 ; vherr_max = -vhbt
3200 vbt = btc%vBT_NN * (vhbt / btc%vh_NN)
3202 vhbt_err = vbt * (btc%FA_v_N0 + btc%vh_crvN * vbt**2) - vhbt
3204 if (abs(vhbt_err) < tol*abs(vhbt))
exit
3205 if (vhbt_err > 0.0)
then ; vbt_max = vbt ; vherr_max = vhbt_err ;
endif
3206 if (vhbt_err < 0.0)
then ; vbt_min = vbt ; vherr_min = vhbt_err ;
endif
3208 derr_dv = btc%FA_v_N0 + 3.0 * btc%vh_crvN * vbt**2
3209 if ((vhbt_err >= derr_dv*(vbt - vbt_min)) .or. &
3210 (-vhbt_err >= derr_dv*(vbt_max - vbt)) .or. (derr_dv <= 0.0))
then
3212 vbt = vbt_max + (vbt_min-vbt_max) * (vherr_max / (vherr_max-vherr_min))
3214 vbt = vbt - vhbt_err / derr_dv
3215 if (abs(vhbt_err) < (0.01*tol)*abs(derr_dv*vbt_min))
exit
3218 elseif (vhbt <= btc%vh_SS)
then
3220 vbt_min = 0.0 ; vherr_min = -vhbt
3221 vbt_max = btc%vBT_SS ; vherr_max = btc%vh_SS - vhbt
3223 vbt = btc%vBT_SS * (vhbt / btc%vh_SS)
3225 vhbt_err = vbt * (btc%FA_v_S0 + btc%vh_crvS * vbt**2) - vhbt
3227 if (abs(vhbt_err) < tol*abs(vhbt))
exit
3228 if (vhbt_err > 0.0)
then ; vbt_max = vbt ; vherr_max = vhbt_err ;
endif
3229 if (vhbt_err < 0.0)
then ; vbt_min = vbt ; vherr_min = vhbt_err ;
endif
3231 derr_dv = btc%FA_v_S0 + 3.0 * btc%vh_crvS * vbt**2
3232 if ((vhbt_err >= derr_dv*(vbt - vbt_min)) .or. &
3233 (-vhbt_err >= derr_dv*(vbt_max - vbt)) .or. (derr_dv <= 0.0))
then
3235 vbt = vbt_min + (vbt_max-vbt_min) * (-vherr_min / (vherr_max-vherr_min))
3237 vbt = vbt - vhbt_err / derr_dv
3238 if (abs(vhbt_err) < (0.01*tol)*(vbt_max*derr_dv))
exit
3242 vbt = btc%vBT_SS + (vhbt - btc%vh_SS) / btc%FA_v_SS
3245 if (
present(guess))
then
3246 dvel = abs(vbt) - vs1*abs(guess)
3247 if (dvel > 0.0)
then
3248 if (dvel < 40.0 * (abs(guess)*(vs2-vs1)) )
then
3249 vsr = vs2 - (vs2-vs1) * exp(-dvel / (abs(guess)*(vs2-vs1)))
3253 vbt = sign(guess * vsr, vbt)
3257 end function vhbt_to_vbt
3261 subroutine set_local_bt_cont_types(BT_cont, BTCL_u, BTCL_v, G, US, MS, BT_Domain, halo)
3275 integer,
optional,
intent(in) :: halo
3278 real,
dimension(SZIBW_(MS),SZJW_(MS)) :: &
3279 u_polarity, uBT_EE, uBT_WW, FA_u_EE, FA_u_E0, FA_u_W0, FA_u_WW
3280 real,
dimension(SZIW_(MS),SZJBW_(MS)) :: &
3281 v_polarity, vBT_NN, vBT_SS, FA_v_NN, FA_v_N0, FA_v_S0, FA_v_SS
3282 real,
parameter :: C1_3 = 1.0/3.0
3283 integer :: i, j, is, ie, js, je, hs
3284 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
3285 hs = 1 ;
if (
present(halo)) hs = max(halo,0)
3292 do j=js-hs,je+hs ;
do i=is-hs-1,ie+hs
3293 u_polarity(i,j) = 1.0
3294 ubt_ee(i,j) = 0.0 ; ubt_ww(i,j) = 0.0
3295 fa_u_ee(i,j) = 0.0 ; fa_u_e0(i,j) = 0.0 ; fa_u_w0(i,j) = 0.0 ; fa_u_ww(i,j) = 0.0
3298 do j=js-hs-1,je+hs ;
do i=is-hs,ie+hs
3299 v_polarity(i,j) = 1.0
3300 vbt_nn(i,j) = 0.0 ; vbt_ss(i,j) = 0.0
3301 fa_v_nn(i,j) = 0.0 ; fa_v_n0(i,j) = 0.0 ; fa_v_s0(i,j) = 0.0 ; fa_v_ss(i,j) = 0.0
3304 do j=js,je;
do i=is-1,ie
3305 ubt_ee(i,j) = bt_cont%uBT_EE(i,j) ; ubt_ww(i,j) = bt_cont%uBT_WW(i,j)
3306 fa_u_ee(i,j) = bt_cont%FA_u_EE(i,j) ; fa_u_e0(i,j) = bt_cont%FA_u_E0(i,j)
3307 fa_u_w0(i,j) = bt_cont%FA_u_W0(i,j) ; fa_u_ww(i,j) = bt_cont%FA_u_WW(i,j)
3310 do j=js-1,je;
do i=is,ie
3311 vbt_nn(i,j) = bt_cont%vBT_NN(i,j) ; vbt_ss(i,j) = bt_cont%vBT_SS(i,j)
3312 fa_v_nn(i,j) = bt_cont%FA_v_NN(i,j) ; fa_v_n0(i,j) = bt_cont%FA_v_N0(i,j)
3313 fa_v_s0(i,j) = bt_cont%FA_v_S0(i,j) ; fa_v_ss(i,j) = bt_cont%FA_v_SS(i,j)
3317 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
3318 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
3320 call create_group_pass(bt_cont%pass_polarity_BT, u_polarity, v_polarity, bt_domain)
3324 call create_group_pass(bt_cont%pass_FA_uv, fa_u_ee, fa_v_nn, bt_domain, to_all+scalar_pair)
3325 call create_group_pass(bt_cont%pass_FA_uv, fa_u_e0, fa_v_n0, bt_domain, to_all+scalar_pair)
3326 call create_group_pass(bt_cont%pass_FA_uv, fa_u_w0, fa_v_s0, bt_domain, to_all+scalar_pair)
3327 call create_group_pass(bt_cont%pass_FA_uv, fa_u_ww, fa_v_ss, bt_domain, to_all+scalar_pair)
3330 call do_group_pass(bt_cont%pass_polarity_BT, bt_domain)
3331 call do_group_pass(bt_cont%pass_FA_uv, bt_domain)
3332 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
3333 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
3340 do j=js-hs,je+hs ;
do i=is-hs-1,ie+hs
3341 btcl_u(i,j)%FA_u_EE = fa_u_ee(i,j) ; btcl_u(i,j)%FA_u_E0 = fa_u_e0(i,j)
3342 btcl_u(i,j)%FA_u_W0 = fa_u_w0(i,j) ; btcl_u(i,j)%FA_u_WW = fa_u_ww(i,j)
3343 btcl_u(i,j)%uBT_EE = ubt_ee(i,j) ; btcl_u(i,j)%uBT_WW = ubt_ww(i,j)
3345 if (u_polarity(i,j) < 0.0)
then
3346 call swap(btcl_u(i,j)%FA_u_EE, btcl_u(i,j)%FA_u_WW)
3347 call swap(btcl_u(i,j)%FA_u_E0, btcl_u(i,j)%FA_u_W0)
3348 call swap(btcl_u(i,j)%uBT_EE, btcl_u(i,j)%uBT_WW)
3351 btcl_u(i,j)%uh_EE = btcl_u(i,j)%uBT_EE * &
3352 (c1_3 * (2.0*btcl_u(i,j)%FA_u_E0 + btcl_u(i,j)%FA_u_EE))
3353 btcl_u(i,j)%uh_WW = btcl_u(i,j)%uBT_WW * &
3354 (c1_3 * (2.0*btcl_u(i,j)%FA_u_W0 + btcl_u(i,j)%FA_u_WW))
3356 btcl_u(i,j)%uh_crvE = 0.0 ; btcl_u(i,j)%uh_crvW = 0.0
3357 if (abs(btcl_u(i,j)%uBT_WW) > 0.0) btcl_u(i,j)%uh_crvW = &
3358 (c1_3 * (btcl_u(i,j)%FA_u_WW - btcl_u(i,j)%FA_u_W0)) / btcl_u(i,j)%uBT_WW**2
3359 if (abs(btcl_u(i,j)%uBT_EE) > 0.0) btcl_u(i,j)%uh_crvE = &
3360 (c1_3 * (btcl_u(i,j)%FA_u_EE - btcl_u(i,j)%FA_u_E0)) / btcl_u(i,j)%uBT_EE**2
3363 do j=js-hs-1,je+hs ;
do i=is-hs,ie+hs
3364 btcl_v(i,j)%FA_v_NN = fa_v_nn(i,j) ; btcl_v(i,j)%FA_v_N0 = fa_v_n0(i,j)
3365 btcl_v(i,j)%FA_v_S0 = fa_v_s0(i,j) ; btcl_v(i,j)%FA_v_SS = fa_v_ss(i,j)
3366 btcl_v(i,j)%vBT_NN = vbt_nn(i,j) ; btcl_v(i,j)%vBT_SS = vbt_ss(i,j)
3368 if (v_polarity(i,j) < 0.0)
then
3369 call swap(btcl_v(i,j)%FA_v_NN, btcl_v(i,j)%FA_v_SS)
3370 call swap(btcl_v(i,j)%FA_v_N0, btcl_v(i,j)%FA_v_S0)
3371 call swap(btcl_v(i,j)%vBT_NN, btcl_v(i,j)%vBT_SS)
3374 btcl_v(i,j)%vh_NN = btcl_v(i,j)%vBT_NN * &
3375 (c1_3 * (2.0*btcl_v(i,j)%FA_v_N0 + btcl_v(i,j)%FA_v_NN))
3376 btcl_v(i,j)%vh_SS = btcl_v(i,j)%vBT_SS * &
3377 (c1_3 * (2.0*btcl_v(i,j)%FA_v_S0 + btcl_v(i,j)%FA_v_SS))
3379 btcl_v(i,j)%vh_crvN = 0.0 ; btcl_v(i,j)%vh_crvS = 0.0
3380 if (abs(btcl_v(i,j)%vBT_SS) > 0.0) btcl_v(i,j)%vh_crvS = &
3381 (c1_3 * (btcl_v(i,j)%FA_v_SS - btcl_v(i,j)%FA_v_S0)) / btcl_v(i,j)%vBT_SS**2
3382 if (abs(btcl_v(i,j)%vBT_NN) > 0.0) btcl_v(i,j)%vh_crvN = &
3383 (c1_3 * (btcl_v(i,j)%FA_v_NN - btcl_v(i,j)%FA_v_N0)) / btcl_v(i,j)%vBT_NN**2
3386 end subroutine set_local_bt_cont_types
3391 subroutine adjust_local_bt_cont_types(ubt, uhbt, vbt, vhbt, BTCL_u, BTCL_v, &
3394 real,
dimension(SZIBW_(MS),SZJW_(MS)), &
3396 real,
dimension(SZIBW_(MS),SZJW_(MS)), &
3399 real,
dimension(SZIW_(MS),SZJBW_(MS)), &
3401 real,
dimension(SZIW_(MS),SZJBW_(MS)), &
3405 intent(out) :: BTCL_u
3407 intent(out) :: BTCL_v
3410 integer,
optional,
intent(in) :: halo
3413 real,
dimension(SZIBW_(MS),SZJW_(MS)) :: &
3414 u_polarity, uBT_EE, uBT_WW, FA_u_EE, FA_u_E0, FA_u_W0, FA_u_WW
3415 real,
dimension(SZIW_(MS),SZJBW_(MS)) :: &
3416 v_polarity, vBT_NN, vBT_SS, FA_v_NN, FA_v_N0, FA_v_S0, FA_v_SS
3417 real,
parameter :: C1_3 = 1.0/3.0
3418 integer :: i, j, is, ie, js, je, hs
3419 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
3420 hs = 1 ;
if (
present(halo)) hs = max(halo,0)
3423 do j=js-hs,je+hs ;
do i=is-hs-1,ie+hs
3424 if ((ubt(i,j) > btcl_u(i,j)%uBT_WW) .and. (uhbt(i,j) > btcl_u(i,j)%uh_WW))
then
3426 btcl_u(i,j)%ubt_WW = ubt(i,j)
3427 if (3.0*uhbt(i,j) < 2.0*ubt(i,j) * btcl_u(i,j)%FA_u_W0)
then
3429 btcl_u(i,j)%uh_crvW = (uhbt(i,j) - ubt(i,j) * btcl_u(i,j)%FA_u_W0) / ubt(i,j)**3
3431 btcl_u(i,j)%FA_u_W0 = 1.5*uhbt(i,j) / ubt(i,j)
3432 btcl_u(i,j)%uh_crvW = -0.5*uhbt(i,j) / ubt(i,j)**3
3434 btcl_u(i,j)%uh_WW = uhbt(i,j)
3437 elseif ((ubt(i,j) < btcl_u(i,j)%uBT_EE) .and. (uhbt(i,j) < btcl_u(i,j)%uh_EE))
then
3439 btcl_u(i,j)%ubt_EE = ubt(i,j)
3440 if (3.0*uhbt(i,j) < 2.0*ubt(i,j) * btcl_u(i,j)%FA_u_E0)
then
3442 btcl_u(i,j)%uh_crvE = (uhbt(i,j) - ubt(i,j) * btcl_u(i,j)%FA_u_E0) / ubt(i,j)**3
3444 btcl_u(i,j)%FA_u_E0 = 1.5*uhbt(i,j) / ubt(i,j)
3445 btcl_u(i,j)%uh_crvE = -0.5*uhbt(i,j) / ubt(i,j)**3
3447 btcl_u(i,j)%uh_EE = uhbt(i,j)
3453 do j=js-hs-1,je+hs ;
do i=is-hs,ie+hs
3454 if ((vbt(i,j) > btcl_v(i,j)%vBT_SS) .and. (vhbt(i,j) > btcl_v(i,j)%vh_SS))
then
3456 btcl_v(i,j)%vbt_SS = vbt(i,j)
3457 if (3.0*vhbt(i,j) < 2.0*vbt(i,j) * btcl_v(i,j)%FA_v_S0)
then
3459 btcl_v(i,j)%vh_crvS = (vhbt(i,j) - vbt(i,j) * btcl_v(i,j)%FA_v_S0) / vbt(i,j)**3
3461 btcl_v(i,j)%FA_v_S0 = 1.5*vhbt(i,j) / (vbt(i,j))
3462 btcl_v(i,j)%vh_crvS = -0.5*vhbt(i,j) / vbt(i,j)**3
3464 btcl_v(i,j)%vh_SS = vhbt(i,j)
3467 elseif ((vbt(i,j) < btcl_v(i,j)%vBT_NN) .and. (vhbt(i,j) < btcl_v(i,j)%vh_NN))
then
3469 btcl_v(i,j)%vbt_NN = vbt(i,j)
3470 if (3.0*vhbt(i,j) < 2.0*vbt(i,j) * btcl_v(i,j)%FA_v_N0)
then
3472 btcl_v(i,j)%vh_crvN = (vhbt(i,j) - vbt(i,j) * btcl_v(i,j)%FA_v_N0) / vbt(i,j)**3
3474 btcl_v(i,j)%FA_v_N0 = 1.5*vhbt(i,j) / (vbt(i,j))
3475 btcl_v(i,j)%vh_crvN = -0.5*vhbt(i,j) / vbt(i,j)**3
3477 btcl_v(i,j)%vh_NN = vhbt(i,j)
3483 end subroutine adjust_local_bt_cont_types
3487 subroutine bt_cont_to_face_areas(BT_cont, Datu, Datv, G, US, MS, halo, maximize)
3492 real,
dimension(MS%isdw-1:MS%iedw,MS%jsdw:MS%jedw), &
3494 real,
dimension(MS%isdw:MS%iedw,MS%jsdw-1:MS%jedw), &
3498 integer,
optional,
intent(in) :: halo
3499 logical,
optional,
intent(in) :: maximize
3504 integer :: i, j, is, ie, js, je, hs
3505 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
3506 hs = 1 ;
if (
present(halo)) hs = max(halo,0)
3507 find_max = .false. ;
if (
present(maximize)) find_max = maximize
3510 do j=js-hs,je+hs ;
do i=is-1-hs,ie+hs
3511 datu(i,j) = max(bt_cont%FA_u_EE(i,j), bt_cont%FA_u_E0(i,j), &
3512 bt_cont%FA_u_W0(i,j), bt_cont%FA_u_WW(i,j))
3514 do j=js-1-hs,je+hs ;
do i=is-hs,ie+hs
3515 datv(i,j) = max(bt_cont%FA_v_NN(i,j), bt_cont%FA_v_N0(i,j), &
3516 bt_cont%FA_v_S0(i,j), bt_cont%FA_v_SS(i,j))
3519 do j=js-hs,je+hs ;
do i=is-1-hs,ie+hs
3520 datu(i,j) = 0.5 * (bt_cont%FA_u_E0(i,j) + bt_cont%FA_u_W0(i,j))
3522 do j=js-1-hs,je+hs ;
do i=is-hs,ie+hs
3523 datv(i,j) = 0.5 * (bt_cont%FA_v_N0(i,j) + bt_cont%FA_v_S0(i,j))
3527 end subroutine bt_cont_to_face_areas
3530 subroutine swap(a,b)
3531 real,
intent(inout) :: a
3532 real,
intent(inout) :: b
3534 tmp = a ; a = b ; b = tmp
3539 subroutine find_face_areas(Datu, Datv, G, GV, US, CS, MS, eta, halo, add_max)
3541 real,
dimension(MS%isdw-1:MS%iedw,MS%jsdw:MS%jedw), &
3543 real,
dimension(MS%isdw:MS%iedw,MS%jsdw-1:MS%jedw), &
3550 real,
dimension(MS%isdw:MS%iedw,MS%jsdw:MS%jedw), &
3551 optional,
intent(in) :: eta
3553 integer,
optional,
intent(in) :: halo
3554 real,
optional,
intent(in) :: add_max
3559 integer :: i, j, is, ie, js, je, hs
3560 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
3561 hs = 1 ;
if (
present(halo)) hs = max(halo,0)
3565 if (
present(eta))
then
3567 if (gv%Boussinesq)
then
3569 do j=js-hs,je+hs ;
do i=is-1-hs,ie+hs
3570 h1 = cs%bathyT(i,j)*gv%Z_to_H + eta(i,j) ; h2 = cs%bathyT(i+1,j)*gv%Z_to_H + eta(i+1,j)
3571 datu(i,j) = 0.0 ;
if ((h1 > 0.0) .and. (h2 > 0.0)) &
3572 datu(i,j) = cs%dy_Cu(i,j) * (2.0 * h1 * h2) / (h1 + h2)
3576 do j=js-1-hs,je+hs ;
do i=is-hs,ie+hs
3577 h1 = cs%bathyT(i,j)*gv%Z_to_H + eta(i,j) ; h2 = cs%bathyT(i,j+1)*gv%Z_to_H + eta(i,j+1)
3578 datv(i,j) = 0.0 ;
if ((h1 > 0.0) .and. (h2 > 0.0)) &
3579 datv(i,j) = cs%dx_Cv(i,j) * (2.0 * h1 * h2) / (h1 + h2)
3584 do j=js-hs,je+hs ;
do i=is-1-hs,ie+hs
3585 datu(i,j) = 0.0 ;
if ((eta(i,j) > 0.0) .and. (eta(i+1,j) > 0.0)) &
3586 datu(i,j) = cs%dy_Cu(i,j) * (2.0 * eta(i,j) * eta(i+1,j)) / &
3587 (eta(i,j) + eta(i+1,j))
3591 do j=js-1-hs,je+hs ;
do i=is-hs,ie+hs
3592 datv(i,j) = 0.0 ;
if ((eta(i,j) > 0.0) .and. (eta(i,j+1) > 0.0)) &
3593 datv(i,j) = cs%dx_Cv(i,j) * (2.0 * eta(i,j) * eta(i,j+1)) / &
3594 (eta(i,j) + eta(i,j+1))
3598 elseif (
present(add_max))
then
3600 do j=js-hs,je+hs ;
do i=is-1-hs,ie+hs
3601 datu(i,j) = cs%dy_Cu(i,j) * gv%Z_to_H * &
3602 (max(cs%bathyT(i+1,j), cs%bathyT(i,j)) + add_max)
3605 do j=js-1-hs,je+hs ;
do i=is-hs,ie+hs
3606 datv(i,j) = cs%dx_Cv(i,j) * gv%Z_to_H * &
3607 (max(cs%bathyT(i,j+1), cs%bathyT(i,j)) + add_max)
3611 do j=js-hs,je+hs ;
do i=is-1-hs,ie+hs
3614 if (cs%bathyT(i+1,j)+cs%bathyT(i,j)>0.) &
3615 datu(i,j) = 2.0*cs%dy_Cu(i,j) * gv%Z_to_H * &
3616 (cs%bathyT(i+1,j) * cs%bathyT(i,j)) / &
3617 (cs%bathyT(i+1,j) + cs%bathyT(i,j))
3620 do j=js-1-hs,je+hs ;
do i=is-hs,ie+hs
3623 if (cs%bathyT(i,j+1)+cs%bathyT(i,j)>0.) &
3624 datv(i,j) = 2.0*cs%dx_Cv(i,j) * gv%Z_to_H * &
3625 (cs%bathyT(i,j+1) * cs%bathyT(i,j)) / &
3626 (cs%bathyT(i,j+1) + cs%bathyT(i,j))
3631 end subroutine find_face_areas
3637 subroutine bt_mass_source(h, eta, set_cor, G, GV, CS)
3640 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
3641 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: eta
3643 logical,
intent(in) :: set_cor
3651 real :: h_tot(szi_(g))
3652 real :: eta_h(szi_(g))
3658 integer :: is, ie, js, je, nz, i, j, k
3659 real,
parameter :: frac_cor = 0.25
3660 real,
parameter :: slow_rate = 0.125
3662 if (.not.
associated(cs))
call mom_error(fatal,
"bt_mass_source: "// &
3663 "Module MOM_barotropic must be initialized before it is used.")
3664 if (.not.cs%split)
return
3666 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
3670 do i=is,ie ; h_tot(i) = h(i,j,1) ;
enddo
3671 if (gv%Boussinesq)
then
3672 do i=is,ie ; eta_h(i) = h(i,j,1) - g%bathyT(i,j)*gv%Z_to_H ;
enddo
3674 do i=is,ie ; eta_h(i) = h(i,j,1) ;
enddo
3676 do k=2,nz ;
do i=is,ie
3677 eta_h(i) = eta_h(i) + h(i,j,k)
3678 h_tot(i) = h_tot(i) + h(i,j,k)
3683 d_eta = eta_h(i) - eta(i,j)
3684 cs%eta_cor(i,j) = d_eta
3688 d_eta = eta_h(i) - eta(i,j)
3689 cs%eta_cor(i,j) = cs%eta_cor(i,j) + d_eta
3694 end subroutine bt_mass_source
3699 subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, &
3700 restart_CS, calc_dtbt, BT_cont, tides_CSp)
3704 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
3706 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
3708 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
3710 real,
dimension(SZI_(G),SZJ_(G)), &
3713 type(time_type),
target,
intent(in) :: time
3715 type(
diag_ctrl),
target,
intent(inout) :: diag
3720 logical,
intent(out) :: calc_dtbt
3727 pointer :: tides_csp
3731 #include "version_variable.h"
3733 character(len=40) :: mdl =
"MOM_barotropic"
3734 real :: datu(szibs_(g),szj_(g))
3735 real :: datv(szi_(g),szjbs_(g))
3736 real :: gtot_estimate
3739 real :: dtbt_input, dtbt_tmp
3740 real :: wave_drag_scale
3742 character(len=200) :: inputdir
3743 character(len=200) :: wave_drag_file
3745 character(len=80) :: wave_drag_var
3751 real,
allocatable,
dimension(:,:) :: lin_drag_h
3753 type(group_pass_type) :: pass_static_data, pass_q_d_cor
3754 type(group_pass_type) :: pass_bt_hbt_btav, pass_a_polarity
3755 logical :: apply_bt_drag, use_bt_cont_type
3756 character(len=48) :: thickness_units, flux_units
3757 character*(40) :: hvel_str
3758 integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz
3759 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
3760 integer :: isdw, iedw, jsdw, jedw
3762 integer :: wd_halos(2), bt_halo_sz
3763 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
3764 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
3765 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
3766 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
3767 ms%isdw = g%isd ; ms%iedw = g%ied ; ms%jsdw = g%jsd ; ms%jedw = g%jed
3769 if (cs%module_is_initialized)
then
3770 call mom_error(warning,
"barotropic_init called with a control structure "// &
3771 "that has already been initialized.")
3774 cs%module_is_initialized = .true.
3776 cs%diag => diag ; cs%Time => time
3777 if (
present(tides_csp))
then
3778 if (
associated(tides_csp)) cs%tides_CSp => tides_csp
3783 call get_param(param_file, mdl,
"SPLIT", cs%split, &
3784 "Use the split time stepping if true.", default=.true.)
3785 if (.not.cs%split)
return
3787 call get_param(param_file, mdl,
"BOUND_BT_CORRECTION", cs%bound_BT_corr, &
3788 "If true, the corrective pseudo mass-fluxes into the "//&
3789 "barotropic solver are limited to values that require "//&
3790 "less than maxCFL_BT_cont to be accommodated.",default=.false.)
3791 call get_param(param_file, mdl,
"BT_CONT_CORR_BOUNDS", cs%BT_cont_bounds, &
3792 "If true, and BOUND_BT_CORRECTION is true, use the "//&
3793 "BT_cont_type variables to set limits determined by "//&
3794 "MAXCFL_BT_CONT on the CFL number of the velocities "//&
3795 "that are likely to be driven by the corrective mass fluxes.", &
3797 call get_param(param_file, mdl,
"ADJUST_BT_CONT", cs%adjust_BT_cont, &
3798 "If true, adjust the curve fit to the BT_cont type "//&
3799 "that is used by the barotropic solver to match the "//&
3800 "transport about which the flow is being linearized.", default=.false.)
3801 call get_param(param_file, mdl,
"GRADUAL_BT_ICS", cs%gradual_BT_ICs, &
3802 "If true, adjust the initial conditions for the "//&
3803 "barotropic solver to the values from the layered "//&
3804 "solution over a whole timestep instead of instantly. "//&
3805 "This is a decent approximation to the inclusion of "//&
3806 "sum(u dh_dt) while also correcting for truncation errors.", &
3808 call get_param(param_file, mdl,
"BT_USE_VISC_REM_U_UH0", cs%visc_rem_u_uh0, &
3809 "If true, use the viscous remnants when estimating the "//&
3810 "barotropic velocities that were used to calculate uh0 "//&
3811 "and vh0. False is probably the better choice.", default=.false.)
3812 call get_param(param_file, mdl,
"BT_USE_WIDE_HALOS", cs%use_wide_halos, &
3813 "If true, use wide halos and march in during the "//&
3814 "barotropic time stepping for efficiency.", default=.true., &
3816 call get_param(param_file, mdl,
"BTHALO", bt_halo_sz, &
3817 "The minimum halo size for the barotropic solver.", default=0, &
3819 #ifdef STATIC_MEMORY_
3820 if ((bt_halo_sz > 0) .and. (bt_halo_sz /= bthalo_))
call mom_error(fatal, &
3821 "barotropic_init: Run-time values of BTHALO must agree with the "//&
3822 "macro BTHALO_ with STATIC_MEMORY_.")
3823 wd_halos(1) = whaloi_+nihalo_ ; wd_halos(2) = whaloj_+njhalo_
3825 wd_halos(1) = bt_halo_sz; wd_halos(2) = bt_halo_sz
3827 call log_param(param_file, mdl,
"!BT x-halo", wd_halos(1), &
3828 "The barotropic x-halo size that is actually used.", &
3830 call log_param(param_file, mdl,
"!BT y-halo", wd_halos(2), &
3831 "The barotropic y-halo size that is actually used.", &
3834 call get_param(param_file, mdl,
"USE_BT_CONT_TYPE", use_bt_cont_type, &
3835 "If true, use a structure with elements that describe "//&
3836 "effective face areas from the summed continuity solver "//&
3837 "as a function the barotropic flow in coupling between "//&
3838 "the barotropic and baroclinic flow. This is only used "//&
3839 "if SPLIT is true. \n", default=.true.)
3840 call get_param(param_file, mdl,
"NONLINEAR_BT_CONTINUITY", &
3841 cs%Nonlinear_continuity, &
3842 "If true, use nonlinear transports in the barotropic "//&
3843 "continuity equation. This does not apply if "//&
3844 "USE_BT_CONT_TYPE is true.", default=.false.)
3845 cs%Nonlin_cont_update_period = 1
3846 if (cs%Nonlinear_continuity) &
3847 call get_param(param_file, mdl,
"NONLIN_BT_CONT_UPDATE_PERIOD", &
3848 cs%Nonlin_cont_update_period, &
3849 "If NONLINEAR_BT_CONTINUITY is true, this is the number "//&
3850 "of barotropic time steps between updates to the face "//&
3851 "areas, or 0 to update only before the barotropic stepping.",&
3852 units=
"nondim", default=1)
3853 call get_param(param_file, mdl,
"BT_PROJECT_VELOCITY", cs%BT_project_velocity,&
3854 "If true, step the barotropic velocity first and project "//&
3855 "out the velocity tendency by 1+BEBT when calculating the "//&
3856 "transport. The default (false) is to use a predictor "//&
3857 "continuity step to find the pressure field, and then "//&
3858 "to do a corrector continuity step using a weighted "//&
3859 "average of the old and new velocities, with weights "//&
3860 "of (1-BEBT) and BEBT.", default=.false.)
3862 call get_param(param_file, mdl,
"DYNAMIC_SURFACE_PRESSURE", cs%dynamic_psurf, &
3863 "If true, add a dynamic pressure due to a viscous ice "//&
3864 "shelf, for instance.", default=.false.)
3865 if (cs%dynamic_psurf)
then
3866 call get_param(param_file, mdl,
"ICE_LENGTH_DYN_PSURF", cs%ice_strength_length, &
3867 "The length scale at which the Rayleigh damping rate due "//&
3868 "to the ice strength should be the same as if a Laplacian "//&
3869 "were applied, if DYNAMIC_SURFACE_PRESSURE is true.", &
3870 units=
"m", default=1.0e4, scale=us%m_to_L)
3871 call get_param(param_file, mdl,
"DEPTH_MIN_DYN_PSURF", cs%Dmin_dyn_psurf, &
3872 "The minimum depth to use in limiting the size of the "//&
3873 "dynamic surface pressure for stability, if "//&
3874 "DYNAMIC_SURFACE_PRESSURE is true..", &
3875 units=
"m", default=1.0e-6, scale=us%m_to_Z)
3876 call get_param(param_file, mdl,
"CONST_DYN_PSURF", cs%const_dyn_psurf, &
3877 "The constant that scales the dynamic surface pressure, "//&
3878 "if DYNAMIC_SURFACE_PRESSURE is true. Stable values "//&
3879 "are < ~1.0.", units=
"nondim", default=0.9)
3882 call get_param(param_file, mdl,
"TIDES", cs%tides, &
3883 "If true, apply tidal momentum forcing.", default=.false.)
3884 call get_param(param_file, mdl,
"SADOURNY", cs%Sadourny, &
3885 "If true, the Coriolis terms are discretized with the "//&
3886 "Sadourny (1975) energy conserving scheme, otherwise "//&
3887 "the Arakawa & Hsu scheme is used. If the internal "//&
3888 "deformation radius is not resolved, the Sadourny scheme "//&
3889 "should probably be used.", default=.true.)
3891 call get_param(param_file, mdl,
"BT_THICK_SCHEME", hvel_str, &
3892 "A string describing the scheme that is used to set the "//&
3893 "open face areas used for barotropic transport and the "//&
3894 "relative weights of the accelerations. Valid values are:\n"//&
3895 "\t ARITHMETIC - arithmetic mean layer thicknesses \n"//&
3896 "\t HARMONIC - harmonic mean layer thicknesses \n"//&
3897 "\t HYBRID (the default) - use arithmetic means for \n"//&
3898 "\t layers above the shallowest bottom, the harmonic \n"//&
3899 "\t mean for layers below, and a weighted average for \n"//&
3900 "\t layers that straddle that depth \n"//&
3901 "\t FROM_BT_CONT - use the average thicknesses kept \n"//&
3902 "\t in the h_u and h_v fields of the BT_cont_type", &
3903 default=bt_cont_string)
3904 select case (hvel_str)
3905 case (hybrid_string) ; cs%hvel_scheme = hybrid
3906 case (harmonic_string) ; cs%hvel_scheme = harmonic
3907 case (arithmetic_string) ; cs%hvel_scheme = arithmetic
3908 case (bt_cont_string) ; cs%hvel_scheme = from_bt_cont
3910 call mom_mesg(
'barotropic_init: BT_THICK_SCHEME ="'//trim(hvel_str)//
'"', 0)
3911 call mom_error(fatal,
"barotropic_init: Unrecognized setting "// &
3912 "#define BT_THICK_SCHEME "//trim(hvel_str)//
" found in input file.")
3914 if ((cs%hvel_scheme == from_bt_cont) .and. .not.use_bt_cont_type) &
3915 call mom_error(fatal,
"barotropic_init: BT_THICK_SCHEME FROM_BT_CONT "//&
3916 "can only be used if USE_BT_CONT_TYPE is defined.")
3918 call get_param(param_file, mdl,
"BT_STRONG_DRAG", cs%strong_drag, &
3919 "If true, use a stronger estimate of the retarding "//&
3920 "effects of strong bottom drag, by making it implicit "//&
3921 "with the barotropic time-step instead of implicit with "//&
3922 "the baroclinic time-step and dividing by the number of "//&
3923 "barotropic steps.", default=.false.)
3924 call get_param(param_file, mdl,
"BT_LINEAR_WAVE_DRAG", cs%linear_wave_drag, &
3925 "If true, apply a linear drag to the barotropic velocities, "//&
3926 "using rates set by lin_drag_u & _v divided by the depth of "//&
3927 "the ocean. This was introduced to facilitate tide modeling.", &
3929 call get_param(param_file, mdl,
"BT_WAVE_DRAG_FILE", wave_drag_file, &
3930 "The name of the file with the barotropic linear wave drag "//&
3931 "piston velocities.", default=
"", do_not_log=.not.cs%linear_wave_drag)
3932 call get_param(param_file, mdl,
"BT_WAVE_DRAG_VAR", wave_drag_var, &
3933 "The name of the variable in BT_WAVE_DRAG_FILE with the "//&
3934 "barotropic linear wave drag piston velocities at h points.", &
3935 default=
"rH", do_not_log=.not.cs%linear_wave_drag)
3936 call get_param(param_file, mdl,
"BT_WAVE_DRAG_SCALE", wave_drag_scale, &
3937 "A scaling factor for the barotropic linear wave drag "//&
3938 "piston velocities.", default=1.0, units=
"nondim", &
3939 do_not_log=.not.cs%linear_wave_drag)
3941 call get_param(param_file, mdl,
"CLIP_BT_VELOCITY", cs%clip_velocity, &
3942 "If true, limit any velocity components that exceed "//&
3943 "CFL_TRUNCATE. This should only be used as a desperate "//&
3944 "debugging measure.", default=.false.)
3945 call get_param(param_file, mdl,
"CFL_TRUNCATE", cs%CFL_trunc, &
3946 "The value of the CFL number that will cause velocity "//&
3947 "components to be truncated; instability can occur past 0.5.", &
3948 units=
"nondim", default=0.5, do_not_log=.not.cs%clip_velocity)
3949 call get_param(param_file, mdl,
"MAXVEL", cs%maxvel, &
3950 "The maximum velocity allowed before the velocity "//&
3951 "components are truncated.", units=
"m s-1", default=3.0e8, scale=us%m_s_to_L_T, &
3952 do_not_log=.not.cs%clip_velocity)
3953 call get_param(param_file, mdl,
"MAXCFL_BT_CONT", cs%maxCFL_BT_cont, &
3954 "The maximum permitted CFL number associated with the "//&
3955 "barotropic accelerations from the summed velocities "//&
3956 "times the time-derivatives of thicknesses.", units=
"nondim", &
3958 call get_param(param_file, mdl,
"VEL_UNDERFLOW", cs%vel_underflow, &
3959 "A negligibly small velocity magnitude below which velocity "//&
3960 "components are set to 0. A reasonable value might be "//&
3961 "1e-30 m/s, which is less than an Angstrom divided by "//&
3962 "the age of the universe.", units=
"m s-1", default=0.0, scale=us%m_s_to_L_T)
3964 call get_param(param_file, mdl,
"DT_BT_FILTER", cs%dt_bt_filter, &
3965 "A time-scale over which the barotropic mode solutions "//&
3966 "are filtered, in seconds if positive, or as a fraction "//&
3967 "of DT if negative. When used this can never be taken to "//&
3968 "be longer than 2*dt. Set this to 0 to apply no filtering.", &
3969 units=
"sec or nondim", default=-0.25)
3970 if (cs%dt_bt_filter > 0.0) cs%dt_bt_filter = us%s_to_T*cs%dt_bt_filter
3971 call get_param(param_file, mdl,
"G_BT_EXTRA", cs%G_extra, &
3972 "A nondimensional factor by which gtot is enhanced.", &
3973 units=
"nondim", default=0.0)
3974 call get_param(param_file, mdl,
"SSH_EXTRA", ssh_extra, &
3975 "An estimate of how much higher SSH might get, for use "//&
3976 "in calculating the safe external wave speed. The "//&
3977 "default is the minimum of 10 m or 5% of MAXIMUM_DEPTH.", &
3978 units=
"m", default=min(10.0,0.05*g%max_depth*us%Z_to_m), scale=us%m_to_Z)
3980 call get_param(param_file, mdl,
"DEBUG", cs%debug, &
3981 "If true, write out verbose debugging data.", &
3982 default=.false., debuggingparam=.true.)
3983 call get_param(param_file, mdl,
"DEBUG_BT", cs%debug_bt, &
3984 "If true, write out verbose debugging data within the "//&
3985 "barotropic time-stepping loop. The data volume can be "//&
3986 "quite large if this is true.", default=cs%debug, &
3987 debuggingparam=.true.)
3989 cs%linearized_BT_PV = .true.
3990 call get_param(param_file, mdl,
"BEBT", cs%bebt, &
3991 "BEBT determines whether the barotropic time stepping "//&
3992 "uses the forward-backward time-stepping scheme or a "//&
3993 "backward Euler scheme. BEBT is valid in the range from "//&
3994 "0 (for a forward-backward treatment of nonrotating "//&
3995 "gravity waves) to 1 (for a backward Euler treatment). "//&
3996 "In practice, BEBT must be greater than about 0.05.", &
3997 units=
"nondim", default=0.1)
3998 call get_param(param_file, mdl,
"DTBT", dtbt_input, &
3999 "The barotropic time step, in s. DTBT is only used with "//&
4000 "the split explicit time stepping. To set the time step "//&
4001 "automatically based the maximum stable value use 0, or "//&
4002 "a negative value gives the fraction of the stable value. "//&
4003 "Setting DTBT to 0 is the same as setting it to -0.98. "//&
4004 "The value of DTBT that will actually be used is an "//&
4005 "integer fraction of DT, rounding down.", units=
"s or nondim",&
4007 call get_param(param_file, mdl,
"BT_USE_OLD_CORIOLIS_BRACKET_BUG", &
4008 cs%use_old_coriolis_bracket_bug , &
4009 "If True, use an order of operations that is not bitwise "//&
4010 "rotationally symmetric in the meridional Coriolis term of "//&
4011 "the barotropic solver.", default=.false.)
4014 call clone_mom_domain(g%Domain, cs%BT_Domain, min_halo=wd_halos, symmetric=.true.)
4015 #ifdef STATIC_MEMORY_
4016 if (wd_halos(1) /= whaloi_+nihalo_)
call mom_error(fatal,
"barotropic_init: "//&
4017 "Barotropic x-halo sizes are incorrectly resized with STATIC_MEMORY_.")
4018 if (wd_halos(2) /= whaloj_+njhalo_)
call mom_error(fatal,
"barotropic_init: "//&
4019 "Barotropic y-halo sizes are incorrectly resized with STATIC_MEMORY_.")
4021 if (bt_halo_sz > 0)
then
4022 if (wd_halos(1) > bt_halo_sz) &
4023 call mom_mesg(
"barotropic_init: barotropic x-halo size increased.", 3)
4024 if (wd_halos(2) > bt_halo_sz) &
4025 call mom_mesg(
"barotropic_init: barotropic y-halo size increased.", 3)
4029 cs%isdw = g%isc-wd_halos(1) ; cs%iedw = g%iec+wd_halos(1)
4030 cs%jsdw = g%jsc-wd_halos(2) ; cs%jedw = g%jec+wd_halos(2)
4031 isdw = cs%isdw ; iedw = cs%iedw ; jsdw = cs%jsdw ; jedw = cs%jedw
4033 alloc_(cs%frhatu(isdb:iedb,jsd:jed,nz)) ; alloc_(cs%frhatv(isd:ied,jsdb:jedb,nz))
4034 alloc_(cs%eta_cor(isd:ied,jsd:jed))
4035 if (cs%bound_BT_corr)
then
4036 alloc_(cs%eta_cor_bound(isd:ied,jsd:jed)) ; cs%eta_cor_bound(:,:) = 0.0
4038 alloc_(cs%IDatu(isdb:iedb,jsd:jed)) ; alloc_(cs%IDatv(isd:ied,jsdb:jedb))
4040 alloc_(cs%ua_polarity(isdw:iedw,jsdw:jedw))
4041 alloc_(cs%va_polarity(isdw:iedw,jsdw:jedw))
4043 cs%frhatu(:,:,:) = 0.0 ; cs%frhatv(:,:,:) = 0.0
4044 cs%eta_cor(:,:) = 0.0
4045 cs%IDatu(:,:) = 0.0 ; cs%IDatv(:,:) = 0.0
4047 cs%ua_polarity(:,:) = 1.0 ; cs%va_polarity(:,:) = 1.0
4048 call create_group_pass(pass_a_polarity, cs%ua_polarity, cs%va_polarity, cs%BT_domain, to_all, agrid)
4049 call do_group_pass(pass_a_polarity, cs%BT_domain)
4051 if (use_bt_cont_type) &
4052 call alloc_bt_cont_type(bt_cont, g, (cs%hvel_scheme == from_bt_cont))
4055 allocate(cs%debug_BT_HI)
4056 cs%debug_BT_HI%isc=g%isc
4057 cs%debug_BT_HI%iec=g%iec
4058 cs%debug_BT_HI%jsc=g%jsc
4059 cs%debug_BT_HI%jec=g%jec
4060 cs%debug_BT_HI%IscB=g%isc-1
4061 cs%debug_BT_HI%IecB=g%iec
4062 cs%debug_BT_HI%JscB=g%jsc-1
4063 cs%debug_BT_HI%JecB=g%jec
4064 cs%debug_BT_HI%isd=cs%isdw
4065 cs%debug_BT_HI%ied=cs%iedw
4066 cs%debug_BT_HI%jsd=cs%jsdw
4067 cs%debug_BT_HI%jed=cs%jedw
4068 cs%debug_BT_HI%IsdB=cs%isdw-1
4069 cs%debug_BT_HI%IedB=cs%iedw
4070 cs%debug_BT_HI%JsdB=cs%jsdw-1
4071 cs%debug_BT_HI%JedB=cs%jedw
4075 alloc_(cs%IareaT(cs%isdw:cs%iedw,cs%jsdw:cs%jedw)) ; cs%IareaT(:,:) = 0.0
4076 alloc_(cs%bathyT(cs%isdw:cs%iedw,cs%jsdw:cs%jedw)) ; cs%bathyT(:,:) = gv%Angstrom_m
4077 alloc_(cs%IdxCu(cs%isdw-1:cs%iedw,cs%jsdw:cs%jedw)) ; cs%IdxCu(:,:) = 0.0
4078 alloc_(cs%IdyCv(cs%isdw:cs%iedw,cs%jsdw-1:cs%jedw)) ; cs%IdyCv(:,:) = 0.0
4079 alloc_(cs%dy_Cu(cs%isdw-1:cs%iedw,cs%jsdw:cs%jedw)) ; cs%dy_Cu(:,:) = 0.0
4080 alloc_(cs%dx_Cv(cs%isdw:cs%iedw,cs%jsdw-1:cs%jedw)) ; cs%dx_Cv(:,:) = 0.0
4081 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
4082 cs%IareaT(i,j) = us%L_to_m**2*g%IareaT(i,j)
4083 cs%bathyT(i,j) = g%bathyT(i,j)
4088 do j=g%jsd,g%jed ;
do i=g%IsdB,g%IedB
4089 cs%IdxCu(i,j) = us%L_to_m*g%IdxCu(i,j) ; cs%dy_Cu(i,j) = us%m_to_L*g%dy_Cu(i,j)
4091 do j=g%JsdB,g%JedB ;
do i=g%isd,g%ied
4092 cs%IdyCv(i,j) = us%L_to_m*g%IdyCv(i,j) ; cs%dx_Cv(i,j) = us%m_to_L*g%dx_Cv(i,j)
4096 call create_group_pass(pass_static_data, cs%IdxCu, cs%IdyCv, cs%BT_domain, to_all+scalar_pair)
4097 call create_group_pass(pass_static_data, cs%dy_Cu, cs%dx_Cv, cs%BT_domain, to_all+scalar_pair)
4098 call do_group_pass(pass_static_data, cs%BT_domain)
4100 if (cs%linearized_BT_PV)
then
4101 alloc_(cs%q_D(cs%isdw-1:cs%iedw,cs%jsdw-1:cs%jedw))
4102 alloc_(cs%D_u_Cor(cs%isdw-1:cs%iedw,cs%jsdw:cs%jedw))
4103 alloc_(cs%D_v_Cor(cs%isdw:cs%iedw,cs%jsdw-1:cs%jedw))
4104 cs%q_D(:,:) = 0.0 ; cs%D_u_Cor(:,:) = 0.0 ; cs%D_v_Cor(:,:) = 0.0
4105 do j=js,je ;
do i=is-1,ie
4106 cs%D_u_Cor(i,j) = 0.5 * (g%bathyT(i+1,j) + g%bathyT(i,j))
4108 do j=js-1,je ;
do i=is,ie
4109 cs%D_v_Cor(i,j) = 0.5 * (g%bathyT(i,j+1) + g%bathyT(i,j))
4111 do j=js-1,je ;
do i=is-1,ie
4112 if (g%mask2dT(i,j)+g%mask2dT(i,j+1)+g%mask2dT(i+1,j)+g%mask2dT(i+1,j+1)>0.)
then
4113 cs%q_D(i,j) = 0.25 * g%CoriolisBu(i,j) * &
4114 ((g%areaT(i,j) + g%areaT(i+1,j+1)) + (g%areaT(i+1,j) + g%areaT(i,j+1))) / &
4115 ((g%areaT(i,j) * g%bathyT(i,j) + g%areaT(i+1,j+1) * g%bathyT(i+1,j+1)) + &
4116 (g%areaT(i+1,j) * g%bathyT(i+1,j) + g%areaT(i,j+1) * g%bathyT(i,j+1)) )
4123 call create_group_pass(pass_q_d_cor, cs%q_D, cs%BT_Domain, to_all, position=corner)
4126 call do_group_pass(pass_q_d_cor, cs%BT_Domain)
4129 if (cs%linear_wave_drag)
then
4130 alloc_(cs%lin_drag_u(isdb:iedb,jsd:jed)) ; cs%lin_drag_u(:,:) = 0.0
4131 alloc_(cs%lin_drag_v(isd:ied,jsdb:jedb)) ; cs%lin_drag_v(:,:) = 0.0
4133 if (len_trim(wave_drag_file) > 0)
then
4134 inputdir =
"." ;
call get_param(param_file, mdl,
"INPUTDIR", inputdir)
4135 wave_drag_file = trim(slasher(inputdir))//trim(wave_drag_file)
4136 call log_param(param_file, mdl,
"INPUTDIR/BT_WAVE_DRAG_FILE", wave_drag_file)
4138 allocate(lin_drag_h(isd:ied,jsd:jed)) ; lin_drag_h(:,:) = 0.0
4140 call mom_read_data(wave_drag_file, wave_drag_var, lin_drag_h, g%Domain, scale=us%m_to_Z*us%T_to_s)
4141 call pass_var(lin_drag_h, g%Domain)
4142 do j=js,je ;
do i=is-1,ie
4143 cs%lin_drag_u(i,j) = (gv%Z_to_H * wave_drag_scale) * &
4144 0.5 * (lin_drag_h(i,j) + lin_drag_h(i+1,j))
4146 do j=js-1,je ;
do i=is,ie
4147 cs%lin_drag_v(i,j) = (gv%Z_to_H * wave_drag_scale) * &
4148 0.5 * (lin_drag_h(i,j) + lin_drag_h(i,j+1))
4150 deallocate(lin_drag_h)
4154 cs%dtbt_fraction = 0.98 ;
if (dtbt_input < 0.0) cs%dtbt_fraction = -dtbt_input
4161 do k=1,g%ke ; gtot_estimate = gtot_estimate + gv%g_prime(k) ;
enddo
4162 call set_dtbt(g, gv, us, cs, gtot_est=gtot_estimate, ssh_add=ssh_extra)
4164 if (dtbt_input > 0.0)
then
4165 cs%dtbt = dtbt_input
4166 elseif (dtbt_tmp > 0.0)
then
4169 if ((dtbt_tmp > 0.0) .and. (dtbt_input > 0.0)) calc_dtbt = .false.
4171 call log_param(param_file, mdl,
"DTBT as used", cs%dtbt)
4172 call log_param(param_file, mdl,
"estimated maximum DTBT", cs%dtbt_max)
4177 if (gv%Boussinesq)
then
4178 thickness_units =
"m" ; flux_units =
"m3 s-1"
4180 thickness_units =
"kg m-2" ; flux_units =
"kg s-1"
4183 cs%id_PFu_bt = register_diag_field(
'ocean_model',
'PFuBT', diag%axesCu1, time, &
4184 'Zonal Anomalous Barotropic Pressure Force Force Acceleration',
'm s-2', conversion=us%L_T2_to_m_s2)
4185 cs%id_PFv_bt = register_diag_field(
'ocean_model',
'PFvBT', diag%axesCv1, time, &
4186 'Meridional Anomalous Barotropic Pressure Force Acceleration',
'm s-2', conversion=us%L_T2_to_m_s2)
4187 cs%id_Coru_bt = register_diag_field(
'ocean_model',
'CoruBT', diag%axesCu1, time, &
4188 'Zonal Barotropic Coriolis Acceleration',
'm s-2', conversion=us%L_T2_to_m_s2)
4189 cs%id_Corv_bt = register_diag_field(
'ocean_model',
'CorvBT', diag%axesCv1, time, &
4190 'Meridional Barotropic Coriolis Acceleration',
'm s-2', conversion=us%L_T2_to_m_s2)
4191 cs%id_uaccel = register_diag_field(
'ocean_model',
'u_accel_bt', diag%axesCu1, time, &
4192 'Barotropic zonal acceleration',
'm s-2', conversion=us%L_T2_to_m_s2)
4193 cs%id_vaccel = register_diag_field(
'ocean_model',
'v_accel_bt', diag%axesCv1, time, &
4194 'Barotropic meridional acceleration',
'm s-2', conversion=us%L_T2_to_m_s2)
4195 cs%id_ubtforce = register_diag_field(
'ocean_model',
'ubtforce', diag%axesCu1, time, &
4196 'Barotropic zonal acceleration from baroclinic terms',
'm s-2', conversion=us%L_T2_to_m_s2)
4197 cs%id_vbtforce = register_diag_field(
'ocean_model',
'vbtforce', diag%axesCv1, time, &
4198 'Barotropic meridional acceleration from baroclinic terms',
'm s-2', conversion=us%L_T2_to_m_s2)
4200 cs%id_eta_bt = register_diag_field(
'ocean_model',
'eta_bt', diag%axesT1, time, &
4201 'Barotropic end SSH', thickness_units)
4202 cs%id_ubt = register_diag_field(
'ocean_model',
'ubt', diag%axesCu1, time, &
4203 'Barotropic end zonal velocity',
'm s-1', conversion=us%L_T_to_m_s)
4204 cs%id_vbt = register_diag_field(
'ocean_model',
'vbt', diag%axesCv1, time, &
4205 'Barotropic end meridional velocity',
'm s-1', conversion=us%L_T_to_m_s)
4206 cs%id_eta_st = register_diag_field(
'ocean_model',
'eta_st', diag%axesT1, time, &
4207 'Barotropic start SSH', thickness_units)
4208 cs%id_ubt_st = register_diag_field(
'ocean_model',
'ubt_st', diag%axesCu1, time, &
4209 'Barotropic start zonal velocity',
'm s-1', conversion=us%L_T_to_m_s)
4210 cs%id_vbt_st = register_diag_field(
'ocean_model',
'vbt_st', diag%axesCv1, time, &
4211 'Barotropic start meridional velocity',
'm s-1', conversion=us%L_T_to_m_s)
4212 cs%id_ubtav = register_diag_field(
'ocean_model',
'ubtav', diag%axesCu1, time, &
4213 'Barotropic time-average zonal velocity',
'm s-1', conversion=us%L_T_to_m_s)
4214 cs%id_vbtav = register_diag_field(
'ocean_model',
'vbtav', diag%axesCv1, time, &
4215 'Barotropic time-average meridional velocity',
'm s-1', conversion=us%L_T_to_m_s)
4216 cs%id_eta_cor = register_diag_field(
'ocean_model',
'eta_cor', diag%axesT1, time, &
4217 'Corrective mass flux',
'm s-1')
4218 cs%id_visc_rem_u = register_diag_field(
'ocean_model',
'visc_rem_u', diag%axesCuL, time, &
4219 'Viscous remnant at u',
'nondim')
4220 cs%id_visc_rem_v = register_diag_field(
'ocean_model',
'visc_rem_v', diag%axesCvL, time, &
4221 'Viscous remnant at v',
'nondim')
4222 cs%id_gtotn = register_diag_field(
'ocean_model',
'gtot_n', diag%axesT1, time, &
4223 'gtot to North',
'm s-2', conversion=us%L_T_to_m_s**2)
4224 cs%id_gtots = register_diag_field(
'ocean_model',
'gtot_s', diag%axesT1, time, &
4225 'gtot to South',
'm s-2', conversion=us%L_T_to_m_s**2)
4226 cs%id_gtote = register_diag_field(
'ocean_model',
'gtot_e', diag%axesT1, time, &
4227 'gtot to East',
'm s-2', conversion=us%L_T_to_m_s**2)
4228 cs%id_gtotw = register_diag_field(
'ocean_model',
'gtot_w', diag%axesT1, time, &
4229 'gtot to West',
'm s-2', conversion=us%L_T_to_m_s**2)
4230 cs%id_eta_hifreq = register_diag_field(
'ocean_model',
'eta_hifreq', diag%axesT1, time, &
4231 'High Frequency Barotropic SSH', thickness_units)
4232 cs%id_ubt_hifreq = register_diag_field(
'ocean_model',
'ubt_hifreq', diag%axesCu1, time, &
4233 'High Frequency Barotropic zonal velocity',
'm s-1', conversion=us%L_T_to_m_s)
4234 cs%id_vbt_hifreq = register_diag_field(
'ocean_model',
'vbt_hifreq', diag%axesCv1, time, &
4235 'High Frequency Barotropic meridional velocity',
'm s-1', conversion=us%L_T_to_m_s)
4236 cs%id_eta_pred_hifreq = register_diag_field(
'ocean_model',
'eta_pred_hifreq', diag%axesT1, time, &
4237 'High Frequency Predictor Barotropic SSH', thickness_units)
4238 cs%id_uhbt_hifreq = register_diag_field(
'ocean_model',
'uhbt_hifreq', diag%axesCu1, time, &
4239 'High Frequency Barotropic zonal transport',
'm3 s-1')
4240 cs%id_vhbt_hifreq = register_diag_field(
'ocean_model',
'vhbt_hifreq', diag%axesCv1, time, &
4241 'High Frequency Barotropic meridional transport',
'm3 s-1')
4242 cs%id_frhatu = register_diag_field(
'ocean_model',
'frhatu', diag%axesCuL, time, &
4243 'Fractional thickness of layers in u-columns',
'nondim')
4244 cs%id_frhatv = register_diag_field(
'ocean_model',
'frhatv', diag%axesCvL, time, &
4245 'Fractional thickness of layers in v-columns',
'nondim')
4246 cs%id_frhatu1 = register_diag_field(
'ocean_model',
'frhatu1', diag%axesCuL, time, &
4247 'Predictor Fractional thickness of layers in u-columns',
'nondim')
4248 cs%id_frhatv1 = register_diag_field(
'ocean_model',
'frhatv1', diag%axesCvL, time, &
4249 'Predictor Fractional thickness of layers in v-columns',
'nondim')
4250 cs%id_uhbt = register_diag_field(
'ocean_model',
'uhbt', diag%axesCu1, time, &
4251 'Barotropic zonal transport averaged over a baroclinic step',
'm3 s-1')
4252 cs%id_vhbt = register_diag_field(
'ocean_model',
'vhbt', diag%axesCv1, time, &
4253 'Barotropic meridional transport averaged over a baroclinic step',
'm3 s-1')
4255 if (use_bt_cont_type)
then
4256 cs%id_BTC_FA_u_EE = register_diag_field(
'ocean_model',
'BTC_FA_u_EE', diag%axesCu1, time, &
4257 'BTCont type far east face area',
'm2', conversion=us%L_to_m*gv%H_to_m)
4258 cs%id_BTC_FA_u_E0 = register_diag_field(
'ocean_model',
'BTC_FA_u_E0', diag%axesCu1, time, &
4259 'BTCont type near east face area',
'm2', conversion=us%L_to_m*gv%H_to_m)
4260 cs%id_BTC_FA_u_WW = register_diag_field(
'ocean_model',
'BTC_FA_u_WW', diag%axesCu1, time, &
4261 'BTCont type far west face area',
'm2', conversion=us%L_to_m*gv%H_to_m)
4262 cs%id_BTC_FA_u_W0 = register_diag_field(
'ocean_model',
'BTC_FA_u_W0', diag%axesCu1, time, &
4263 'BTCont type near west face area',
'm2', conversion=us%L_to_m*gv%H_to_m)
4264 cs%id_BTC_ubt_EE = register_diag_field(
'ocean_model',
'BTC_ubt_EE', diag%axesCu1, time, &
4265 'BTCont type far east velocity',
'm s-1')
4266 cs%id_BTC_ubt_WW = register_diag_field(
'ocean_model',
'BTC_ubt_WW', diag%axesCu1, time, &
4267 'BTCont type far west velocity',
'm s-1')
4268 cs%id_BTC_FA_v_NN = register_diag_field(
'ocean_model',
'BTC_FA_v_NN', diag%axesCv1, time, &
4269 'BTCont type far north face area',
'm2', conversion=us%L_to_m*gv%H_to_m)
4270 cs%id_BTC_FA_v_N0 = register_diag_field(
'ocean_model',
'BTC_FA_v_N0', diag%axesCv1, time, &
4271 'BTCont type near north face area',
'm2', conversion=us%L_to_m*gv%H_to_m)
4272 cs%id_BTC_FA_v_SS = register_diag_field(
'ocean_model',
'BTC_FA_v_SS', diag%axesCv1, time, &
4273 'BTCont type far south face area',
'm2', conversion=us%L_to_m*gv%H_to_m)
4274 cs%id_BTC_FA_v_S0 = register_diag_field(
'ocean_model',
'BTC_FA_v_S0', diag%axesCv1, time, &
4275 'BTCont type near south face area',
'm2', conversion=us%L_to_m*gv%H_to_m)
4276 cs%id_BTC_vbt_NN = register_diag_field(
'ocean_model',
'BTC_vbt_NN', diag%axesCv1, time, &
4277 'BTCont type far north velocity',
'm s-1', conversion=us%L_T_to_m_s)
4278 cs%id_BTC_vbt_SS = register_diag_field(
'ocean_model',
'BTC_vbt_SS', diag%axesCv1, time, &
4279 'BTCont type far south velocity',
'm s-1', conversion=us%L_T_to_m_s)
4281 cs%id_uhbt0 = register_diag_field(
'ocean_model',
'uhbt0', diag%axesCu1, time, &
4282 'Barotropic zonal transport difference',
'm3 s-1', conversion=gv%H_to_m*us%L_to_m**2*us%s_to_T)
4283 cs%id_vhbt0 = register_diag_field(
'ocean_model',
'vhbt0', diag%axesCv1, time, &
4284 'Barotropic meridional transport difference',
'm3 s-1', conversion=gv%H_to_m*us%L_to_m**2*us%s_to_T)
4286 if (cs%id_frhatu1 > 0)
call safe_alloc_ptr(cs%frhatu1, isdb,iedb,jsd,jed,nz)
4287 if (cs%id_frhatv1 > 0)
call safe_alloc_ptr(cs%frhatv1, isd,ied,jsdb,jedb,nz)
4291 call btcalc(h, g, gv, cs, may_use_default=.true.)
4292 cs%ubtav(:,:) = 0.0 ; cs%vbtav(:,:) = 0.0
4293 do k=1,nz ;
do j=js,je ;
do i=is-1,ie
4294 cs%ubtav(i,j) = cs%ubtav(i,j) + cs%frhatu(i,j,k) * us%m_s_to_L_T*u(i,j,k)
4295 enddo ;
enddo ;
enddo
4296 do k=1,nz ;
do j=js-1,je ;
do i=is,ie
4297 cs%vbtav(i,j) = cs%vbtav(i,j) + cs%frhatv(i,j,k) * us%m_s_to_L_T*v(i,j,k)
4298 enddo ;
enddo ;
enddo
4299 elseif ((us%s_to_T_restart*us%m_to_L_restart /= 0.0) .and. &
4300 (us%m_to_L*us%s_to_T_restart) /= (us%m_to_L_restart*us%s_to_T))
then
4301 do j=js,je ;
do i=is-1,ie ; cs%ubtav(i,j) = vel_rescale * cs%ubtav(i,j) ;
enddo ;
enddo
4302 do j=js-1,je ;
do i=is,ie ; cs%vbtav(i,j) = vel_rescale * cs%vbtav(i,j) ;
enddo ;
enddo
4307 do j=js,je ;
do i=is-1,ie ; cs%ubt_IC(i,j) = cs%ubtav(i,j) ;
enddo ;
enddo
4308 do j=js-1,je ;
do i=is,ie ; cs%vbt_IC(i,j) = cs%vbtav(i,j) ;
enddo ;
enddo
4309 elseif ((us%s_to_T_restart*us%m_to_L_restart /= 0.0) .and. &
4310 (us%m_to_L*us%s_to_T_restart) /= (us%m_to_L_restart*us%s_to_T))
then
4311 vel_rescale = (us%m_to_L*us%s_to_T_restart) / (us%m_to_L_restart*us%s_to_T)
4312 do j=js,je ;
do i=is-1,ie ; cs%ubt_IC(i,j) = vel_rescale * cs%ubt_IC(i,j) ;
enddo ;
enddo
4313 do j=js-1,je ;
do i=is,ie ; cs%vbt_IC(i,j) = vel_rescale * cs%vbt_IC(i,j) ;
enddo ;
enddo
4320 do j=js,je ;
do i=is-1,ie
4321 if (g%mask2dCu(i,j)>0.)
then
4322 cs%IDatu(i,j) = g%mask2dCu(i,j) * 2.0 / (g%bathyT(i+1,j) + g%bathyT(i,j))
4327 do j=js-1,je ;
do i=is,ie
4328 if (g%mask2dCv(i,j)>0.)
then
4329 cs%IDatv(i,j) = g%mask2dCv(i,j) * 2.0 / (g%bathyT(i,j+1) + g%bathyT(i,j))
4343 call find_face_areas(datu, datv, g, gv, us, cs, ms, halo=1)
4344 if (cs%bound_BT_corr)
then
4347 do j=js,je ;
do i=is,ie
4348 cs%eta_cor_bound(i,j) = gv%m_to_H * us%L_to_m**2*g%IareaT(i,j) * 0.1 * cs%maxvel * &
4349 ((datu(i-1,j) + datu(i,j)) + (datv(i,j) + datv(i,j-1)))
4355 do j=js,je ;
do i=is-1,ie ; cs%uhbt_IC(i,j) = cs%ubtav(i,j) * datu(i,j) ;
enddo ;
enddo
4356 do j=js-1,je ;
do i=is,ie ; cs%vhbt_IC(i,j) = cs%vbtav(i,j) * datv(i,j) ;
enddo ;
enddo
4357 elseif ((us%s_to_T_restart * us%m_to_L_restart * gv%m_to_H_restart /= 0.0) .and. &
4358 ((us%s_to_T_restart * us%m_to_L**2 * gv%m_to_H) /= &
4359 (us%s_to_T * us%m_to_L_restart**2 * gv%m_to_H_restart)))
then
4360 uh_rescale = (us%s_to_T_restart * us%m_to_L**2 * gv%m_to_H) / &
4361 (us%s_to_T * us%m_to_L_restart**2 * gv%m_to_H_restart)
4362 do j=js,je ;
do i=is-1,ie ; cs%uhbt_IC(i,j) = uh_rescale * cs%uhbt_IC(i,j) ;
enddo ;
enddo
4363 do j=js-1,je ;
do i=is,ie ; cs%vhbt_IC(i,j) = uh_rescale * cs%vhbt_IC(i,j) ;
enddo ;
enddo
4369 call do_group_pass(pass_bt_hbt_btav, g%Domain)
4372 id_clock_calc_pre = cpu_clock_id(
'(Ocean BT pre-calcs only)', grain=clock_routine)
4373 id_clock_pass_pre = cpu_clock_id(
'(Ocean BT pre-step halo updates)', grain=clock_routine)
4374 id_clock_calc = cpu_clock_id(
'(Ocean BT stepping calcs only)', grain=clock_routine)
4375 id_clock_pass_step = cpu_clock_id(
'(Ocean BT stepping halo updates)', grain=clock_routine)
4376 id_clock_calc_post = cpu_clock_id(
'(Ocean BT post-calcs only)', grain=clock_routine)
4377 id_clock_pass_post = cpu_clock_id(
'(Ocean BT post-step halo updates)', grain=clock_routine)
4378 if (dtbt_input <= 0.0) &
4379 id_clock_sync = cpu_clock_id(
'(Ocean BT global synch)', grain=clock_routine)
4381 end subroutine barotropic_init
4384 subroutine barotropic_get_tav(CS, ubtav, vbtav, G, US)
4387 real,
dimension(SZIB_(G),SZJ_(G)),
intent(inout) :: ubtav
4389 real,
dimension(SZI_(G),SZJB_(G)),
intent(inout) :: vbtav
4395 do j=g%jsc,g%jec ;
do i=g%isc-1,g%iec
4396 ubtav(i,j) = us%L_T_to_m_s*cs%ubtav(i,j)
4399 do j=g%jsc-1,g%jec ;
do i=g%isc,g%iec
4400 vbtav(i,j) = us%L_T_to_m_s*cs%vbtav(i,j)
4403 end subroutine barotropic_get_tav
4407 subroutine barotropic_end(CS)
4409 dealloc_(cs%frhatu) ; dealloc_(cs%frhatv)
4410 dealloc_(cs%IDatu) ; dealloc_(cs%IDatv)
4411 dealloc_(cs%ubtav) ; dealloc_(cs%vbtav)
4412 dealloc_(cs%eta_cor)
4413 dealloc_(cs%ua_polarity) ; dealloc_(cs%va_polarity)
4414 if (cs%bound_BT_corr)
then
4415 dealloc_(cs%eta_cor_bound)
4418 call destroy_bt_obc(cs%BT_OBC)
4421 end subroutine barotropic_end
4425 subroutine register_barotropic_restarts(HI, GV, param_file, CS, restart_CS)
4436 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
4437 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed
4438 isdb = hi%IsdB ; iedb = hi%IedB ; jsdb = hi%JsdB ; jedb = hi%JedB
4440 if (
associated(cs))
then
4441 call mom_error(warning,
"register_barotropic_restarts called with an associated "// &
4442 "control structure.")
4447 alloc_(cs%ubtav(isdb:iedb,jsd:jed)) ; cs%ubtav(:,:) = 0.0
4448 alloc_(cs%vbtav(isd:ied,jsdb:jedb)) ; cs%vbtav(:,:) = 0.0
4449 alloc_(cs%ubt_IC(isdb:iedb,jsd:jed)) ; cs%ubt_IC(:,:) = 0.0
4450 alloc_(cs%vbt_IC(isd:ied,jsdb:jedb)) ; cs%vbt_IC(:,:) = 0.0
4451 alloc_(cs%uhbt_IC(isdb:iedb,jsd:jed)) ; cs%uhbt_IC(:,:) = 0.0
4452 alloc_(cs%vhbt_IC(isd:ied,jsdb:jedb)) ; cs%vhbt_IC(:,:) = 0.0
4454 vd(2) = var_desc(
"ubtav",
"m s-1",
"Time mean barotropic zonal velocity", &
4455 hor_grid=
'u', z_grid=
'1')
4456 vd(3) = var_desc(
"vbtav",
"m s-1",
"Time mean barotropic meridional velocity",&
4457 hor_grid=
'v', z_grid=
'1')
4461 vd(2) = var_desc(
"ubt_IC",
"m s-1", &
4462 longname=
"Next initial condition for the barotropic zonal velocity", &
4463 hor_grid=
'u', z_grid=
'1')
4464 vd(3) = var_desc(
"vbt_IC",
"m s-1", &
4465 longname=
"Next initial condition for the barotropic meridional velocity",&
4466 hor_grid=
'v', z_grid=
'1')
4470 if (gv%Boussinesq)
then
4471 vd(2) = var_desc(
"uhbt_IC",
"m3 s-1", &
4472 longname=
"Next initial condition for the barotropic zonal transport", &
4473 hor_grid=
'u', z_grid=
'1')
4474 vd(3) = var_desc(
"vhbt_IC",
"m3 s-1", &
4475 longname=
"Next initial condition for the barotropic meridional transport",&
4476 hor_grid=
'v', z_grid=
'1')
4478 vd(2) = var_desc(
"uhbt_IC",
"kg s-1", &
4479 longname=
"Next initial condition for the barotropic zonal transport", &
4480 hor_grid=
'u', z_grid=
'1')
4481 vd(3) = var_desc(
"vhbt_IC",
"kg s-1", &
4482 longname=
"Next initial condition for the barotropic meridional transport",&
4483 hor_grid=
'v', z_grid=
'1')
4489 longname=
"Barotropic timestep", units=
"seconds")
4491 end subroutine register_barotropic_restarts