7 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine
26 use mom_time_manager,
only : time_type, real_to_time,
operator(+),
operator(-)
32 implicit none ;
private 34 #include <MOM_memory.h> 39 # define WHALOI_ MAX(BTHALO_-NIHALO_,0) 40 # define WHALOJ_ MAX(BTHALO_-NJHALO_,0) 41 # define NIMEMW_ 1-WHALOI_:NIMEM_+WHALOI_ 42 # define NJMEMW_ 1-WHALOJ_:NJMEM_+WHALOJ_ 43 # define NIMEMBW_ -WHALOI_:NIMEM_+WHALOI_ 44 # define NJMEMBW_ -WHALOJ_:NJMEM_+WHALOJ_ 45 # define SZIW_(G) NIMEMW_ 46 # define SZJW_(G) NJMEMW_ 47 # define SZIBW_(G) NIMEMBW_ 48 # define SZJBW_(G) NJMEMBW_ 54 # define SZIW_(G) G%isdw:G%iedw 55 # define SZJW_(G) G%jsdw:G%jedw 56 # define SZIBW_(G) G%isdw-1:G%iedw 57 # define SZJBW_(G) G%jsdw-1:G%jedw 60 public btcalc, bt_mass_source, btstep, barotropic_init, barotropic_end
61 public register_barotropic_restarts, set_dtbt, barotropic_get_tav
70 real,
dimension(:,:),
pointer :: cg_u => null()
71 real,
dimension(:,:),
pointer :: cg_v => null()
72 real,
dimension(:,:),
pointer :: h_u => null()
73 real,
dimension(:,:),
pointer :: h_v => null()
74 real,
dimension(:,:),
pointer :: uhbt => null()
76 real,
dimension(:,:),
pointer :: vhbt => null()
78 real,
dimension(:,:),
pointer :: ubt_outer => null()
80 real,
dimension(:,:),
pointer :: vbt_outer => null()
82 real,
dimension(:,:),
pointer :: eta_outer_u => null()
84 real,
dimension(:,:),
pointer :: eta_outer_v => null()
86 logical :: apply_u_obcs
87 logical :: apply_v_obcs
89 integer :: is_u_obc, ie_u_obc, js_u_obc, je_u_obc
90 integer :: is_v_obc, ie_v_obc, js_v_obc, je_v_obc
92 logical :: is_alloced = .false.
94 type(group_pass_type) :: pass_uv
95 type(group_pass_type) :: pass_uhvh
96 type(group_pass_type) :: pass_h
97 type(group_pass_type) :: pass_cg
98 type(group_pass_type) :: pass_eta_outer
103 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: frhatu
105 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: frhatv
107 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_) :: idatu
109 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_) :: lin_drag_u
112 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_) :: ubt_ic
115 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_) :: ubtav
117 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_) :: idatv
119 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_) :: lin_drag_v
122 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_) :: vbt_ic
125 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_) :: vbtav
127 real allocable_,
dimension(NIMEM_,NJMEM_) :: eta_cor
131 real allocable_,
dimension(NIMEM_,NJMEM_) :: eta_cor_bound
134 real allocable_,
dimension(NIMEMW_,NJMEMW_) :: &
135 ua_polarity, & !< Test vector components for checking grid polarity.
136 va_polarity, & !< Test vector components for checking grid polarity.
138 real allocable_,
dimension(NIMEMW_,NJMEMW_) :: iareat
141 real allocable_,
dimension(NIMEMBW_,NJMEMW_) :: &
142 d_u_cor, & !< A simply averaged depth at u points [Z ~> m].
143 dy_cu, & !< A copy of G%dy_Cu with wide halos [L ~> m].
145 real allocable_,
dimension(NIMEMW_,NJMEMBW_) :: &
146 d_v_cor, & !< A simply averaged depth at v points [Z ~> m].
147 dx_cv, & !< A copy of G%dx_Cv with wide halos [L ~> m].
149 real allocable_,
dimension(NIMEMBW_,NJMEMBW_) :: &
152 real,
dimension(:,:,:),
pointer :: frhatu1 => null()
153 real,
dimension(:,:,:),
pointer :: frhatv1 => null()
159 real :: dtbt_fraction
166 integer :: nstep_last = 0
174 logical :: bound_bt_corr
179 logical :: gradual_bt_ics
190 logical :: integral_bt_cont
195 logical :: nonlinear_continuity
197 integer :: nonlin_cont_update_period
201 logical :: bt_project_velocity
208 logical :: nonlin_stress
211 real :: bt_coriolis_scale
213 logical :: answers_2018
217 logical :: dynamic_psurf
219 real :: dmin_dyn_psurf
221 real :: ice_strength_length
224 real :: const_dyn_psurf
229 integer :: hvel_scheme
233 logical :: strong_drag
235 logical :: linear_wave_drag
238 logical :: linearized_bt_pv
241 logical :: use_wide_halos
243 logical :: clip_velocity
249 real :: vel_underflow
256 real :: maxcfl_bt_cont
261 logical :: bt_cont_bounds
263 logical :: visc_rem_u_uh0
267 logical :: adjust_bt_cont
270 logical :: use_old_coriolis_bracket_bug
273 type(time_type),
pointer :: time => null()
279 logical :: module_is_initialized = .false.
286 type(group_pass_type) :: pass_q_dcor
287 type(group_pass_type) :: pass_gtot
288 type(group_pass_type) :: pass_tmp_uv
289 type(group_pass_type) :: pass_eta_bt_rem
290 type(group_pass_type) :: pass_force_hbt0_cor_ref
291 type(group_pass_type) :: pass_dat_uv
292 type(group_pass_type) :: pass_eta_ubt
293 type(group_pass_type) :: pass_etaav
294 type(group_pass_type) :: pass_ubt_cor
295 type(group_pass_type) :: pass_ubta_uhbta
296 type(group_pass_type) :: pass_e_anom
299 integer :: id_pfu_bt = -1, id_pfv_bt = -1, id_coru_bt = -1, id_corv_bt = -1
300 integer :: id_ubtforce = -1, id_vbtforce = -1, id_uaccel = -1, id_vaccel = -1
301 integer :: id_visc_rem_u = -1, id_visc_rem_v = -1, id_eta_cor = -1
302 integer :: id_ubt = -1, id_vbt = -1, id_eta_bt = -1, id_ubtav = -1, id_vbtav = -1
303 integer :: id_ubt_st = -1, id_vbt_st = -1, id_eta_st = -1
304 integer :: id_ubtdt = -1, id_vbtdt = -1
305 integer :: id_ubt_hifreq = -1, id_vbt_hifreq = -1, id_eta_hifreq = -1
306 integer :: id_uhbt_hifreq = -1, id_vhbt_hifreq = -1, id_eta_pred_hifreq = -1
307 integer :: id_gtotn = -1, id_gtots = -1, id_gtote = -1, id_gtotw = -1
308 integer :: id_uhbt = -1, id_frhatu = -1, id_vhbt = -1, id_frhatv = -1
309 integer :: id_frhatu1 = -1, id_frhatv1 = -1
311 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
312 integer :: id_btc_ubt_ee = -1, id_btc_ubt_ww = -1
313 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
314 integer :: id_btc_vbt_nn = -1, id_btc_vbt_ss = -1
315 integer :: id_btc_fa_u_rat0 = -1, id_btc_fa_v_rat0 = -1, id_btc_fa_h_rat0 = -1
316 integer :: id_uhbt0 = -1, id_vhbt0 = -1
376 integer :: isdw, iedw, jsdw, jedw
381 integer :: id_clock_sync=-1, id_clock_calc=-1
382 integer :: id_clock_calc_pre=-1, id_clock_calc_post=-1
383 integer :: id_clock_pass_step=-1, id_clock_pass_pre=-1, id_clock_pass_post=-1
387 integer,
parameter :: harmonic = 1
388 integer,
parameter :: arithmetic = 2
389 integer,
parameter :: hybrid = 3
390 integer,
parameter :: from_bt_cont = 4
391 integer,
parameter :: hybrid_bt_cont = 5
392 character*(20),
parameter :: hybrid_string =
"HYBRID" 393 character*(20),
parameter :: harmonic_string =
"HARMONIC" 394 character*(20),
parameter :: arithmetic_string =
"ARITHMETIC" 395 character*(20),
parameter :: bt_cont_string =
"FROM_BT_CONT" 406 subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, &
407 eta_PF_in, U_Cor, V_Cor, accel_layer_u, accel_layer_v, &
408 eta_out, uhbtav, vhbtav, G, GV, US, CS, &
409 visc_rem_u, visc_rem_v, etaav, ADp, OBC, BT_cont, eta_PF_start, &
410 taux_bot, tauy_bot, uh0, vh0, u_uh0, v_vh0)
414 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: u_in
416 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: v_in
418 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: eta_in
420 real,
intent(in) :: dt
421 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: bc_accel_u
423 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: bc_accel_v
426 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: pbce
429 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: eta_pf_in
435 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: u_cor
437 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: v_cor
439 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(out) :: accel_layer_u
441 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(out) :: accel_layer_v
443 real,
dimension(SZI_(G),SZJ_(G)),
intent(out) :: eta_out
445 real,
dimension(SZIB_(G),SZJ_(G)),
intent(out) :: uhbtav
448 real,
dimension(SZI_(G),SZJB_(G)),
intent(out) :: vhbtav
453 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: visc_rem_u
459 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: visc_rem_v
460 real,
dimension(SZI_(G),SZJ_(G)),
optional,
intent(out) :: etaav
467 real,
dimension(:,:),
optional,
pointer :: eta_pf_start
470 real,
dimension(:,:),
optional,
pointer :: taux_bot
472 real,
dimension(:,:),
optional,
pointer :: tauy_bot
474 real,
dimension(:,:,:),
optional,
pointer :: uh0
476 real,
dimension(:,:,:),
optional,
pointer :: u_uh0
478 real,
dimension(:,:,:),
optional,
pointer :: vh0
480 real,
dimension(:,:,:),
optional,
pointer :: v_vh0
484 real :: ubt_cor(szib_(g),szj_(g))
485 real :: vbt_cor(szi_(g),szjb_(g))
487 real :: wt_u(szib_(g),szj_(g),szk_(g))
488 real :: wt_v(szi_(g),szjb_(g),szk_(g))
491 real,
dimension(SZIB_(G),SZJ_(G)) :: &
496 real,
dimension(SZI_(G),SZJB_(G)) :: &
501 real,
dimension(SZI_(G),SZJ_(G)) :: &
508 real :: q(szibw_(cs),szjbw_(cs))
510 real,
dimension(SZIBW_(CS),SZJW_(CS)) :: &
545 real,
dimension(SZIW_(CS),SZJBW_(CS)) :: &
578 real,
target,
dimension(SZIW_(CS),SZJW_(CS)) :: &
582 real,
dimension(:,:),
pointer :: &
585 real,
dimension(SZIW_(CS),SZJW_(CS)) :: &
613 real,
dimension(SZIBW_(CS),SZJW_(CS)) :: &
614 ubt_prev, ubt_sum_prev, ubt_wtd_prev, &
615 uhbt_prev, uhbt_sum_prev, &
618 real,
dimension(SZIW_(CS),SZJBW_(CS)) :: &
619 vbt_prev, vbt_sum_prev, vbt_wtd_prev, &
620 vhbt_prev, vhbt_sum_prev, &
624 real :: mass_accel_to_z
650 logical :: do_hifreq_output
651 logical :: use_bt_cont, do_ave, find_etaav, find_pf, find_cor
652 logical :: integral_bt_cont
654 logical :: ice_is_rigid, nonblock_setup, interp_eta_pf
655 logical :: project_velocity, add_uh0
659 real :: ice_strength = 0.0
666 real :: u_max_cor, v_max_cor
667 real :: uint_cor, vint_cor
670 real :: accel_underflow
675 real,
allocatable,
dimension(:) :: wt_vel, wt_eta, wt_accel, wt_trans, wt_accel2
676 real :: sum_wt_vel, sum_wt_eta, sum_wt_accel, sum_wt_trans
677 real :: i_sum_wt_vel, i_sum_wt_eta, i_sum_wt_accel, i_sum_wt_trans
679 real :: trans_wt1, trans_wt2
682 logical :: apply_obcs, apply_obc_flather, apply_obc_open
684 character(len=200) :: mesg
685 integer :: isv, iev, jsv, jev
687 integer :: isvf, ievf, jsvf, jevf, num_cycles
688 integer :: i, j, k, n
689 integer :: is, ie, js, je, nz, isq, ieq, jsq, jeq
690 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
691 integer :: ioff, joff
694 if (.not.
associated(cs))
call mom_error(fatal, &
695 "btstep: Module MOM_barotropic must be initialized before it is used.")
696 if (.not.cs%split)
return 697 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
698 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
699 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
700 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
701 ms%isdw = cs%isdw ; ms%iedw = cs%iedw ; ms%jsdw = cs%jsdw ; ms%jedw = cs%jedw
702 h_neglect = gv%H_subroundoff
705 accel_underflow = cs%vel_underflow * idt
707 use_bt_cont = .false.
708 if (
present(bt_cont)) use_bt_cont = (
associated(bt_cont))
709 integral_bt_cont = use_bt_cont .and. cs%integral_BT_cont
711 interp_eta_pf = .false.
712 if (
present(eta_pf_start)) interp_eta_pf = (
associated(eta_pf_start))
714 project_velocity = cs%BT_project_velocity
718 if ((.not.use_bt_cont) .and. cs%Nonlinear_continuity .and. &
719 (cs%Nonlin_cont_update_period > 0)) stencil = 2
721 do_ave = query_averaging_enabled(cs%diag)
722 find_etaav =
present(etaav)
723 find_pf = (do_ave .and. ((cs%id_PFu_bt > 0) .or. (cs%id_PFv_bt > 0)))
724 find_cor = (do_ave .and. ((cs%id_Coru_bt > 0) .or. (cs%id_Corv_bt > 0)))
727 if (
present(uh0)) add_uh0 =
associated(uh0)
728 if (add_uh0 .and. .not.(
present(vh0) .and.
present(u_uh0) .and. &
729 present(v_vh0)))
call mom_error(fatal, &
730 "btstep: vh0, u_uh0, and v_vh0 must be present if uh0 is used.")
731 if (add_uh0 .and. .not.(
associated(vh0) .and.
associated(u_uh0) .and. &
732 associated(v_vh0)))
call mom_error(fatal, &
733 "btstep: vh0, u_uh0, and v_vh0 must be associated if uh0 is used.")
736 nonblock_setup = g%nonblocking_updates
738 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
740 apply_obcs = .false. ; cs%BT_OBC%apply_u_OBCs = .false. ; cs%BT_OBC%apply_v_OBCs = .false.
741 apply_obc_open = .false.
742 apply_obc_flather = .false.
743 if (
present(obc))
then ;
if (
associated(obc))
then 744 cs%BT_OBC%apply_u_OBCs = obc%open_u_BCs_exist_globally .or. obc%specified_u_BCs_exist_globally
745 cs%BT_OBC%apply_v_OBCs = obc%open_v_BCs_exist_globally .or. obc%specified_v_BCs_exist_globally
746 apply_obc_flather = open_boundary_query(obc, apply_flather_obc=.true.)
747 apply_obc_open = open_boundary_query(obc, apply_open_obc=.true.)
748 apply_obcs = open_boundary_query(obc, apply_specified_obc=.true.) .or. &
749 apply_obc_flather .or. apply_obc_open
751 if (apply_obc_flather .and. .not.gv%Boussinesq)
call mom_error(fatal, &
752 "btstep: Flather open boundary conditions have not yet been "// &
753 "implemented for a non-Boussinesq model.")
757 if (cs%use_wide_halos) &
758 num_cycles = min((is-cs%isdw) / stencil, (js-cs%jsdw) / stencil)
759 isvf = is - (num_cycles-1)*stencil ; ievf = ie + (num_cycles-1)*stencil
760 jsvf = js - (num_cycles-1)*stencil ; jevf = je + (num_cycles-1)*stencil
762 nstep = ceiling(dt/cs%dtbt - 0.0001)
763 if (is_root_pe() .and. (nstep /= cs%nstep_last))
then 764 write(mesg,
'("btstep is using a dynamic barotropic timestep of ", ES12.6, & 765 & " seconds, max ", ES12.6, ".")') (us%T_to_s*dt/nstep), us%T_to_s*cs%dtbt_max
766 call mom_mesg(mesg, 3)
768 cs%nstep_last = nstep
771 instep = 1.0 / real(nstep)
776 mass_accel_to_z = 1.0 / gv%Rho0
777 mass_to_z = us%m_to_Z / gv%Rho0
780 if (project_velocity)
then 781 trans_wt1 = (1.0 + be_proj); trans_wt2 = -be_proj
783 trans_wt1 = bebt ; trans_wt2 = (1.0-bebt)
786 do_hifreq_output = .false.
787 if ((cs%id_ubt_hifreq > 0) .or. (cs%id_vbt_hifreq > 0) .or. &
788 (cs%id_eta_hifreq > 0) .or. (cs%id_eta_pred_hifreq > 0) .or. &
789 (cs%id_uhbt_hifreq > 0) .or. (cs%id_vhbt_hifreq > 0))
then 790 do_hifreq_output = query_averaging_enabled(cs%diag, time_int_in, time_end_in)
791 if (do_hifreq_output) &
792 time_bt_start = time_end_in - real_to_time(us%T_to_s*dt)
796 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
797 if (.not. cs%linearized_BT_PV)
then 802 if ((isq > is-1) .or. (jsq > js-1)) &
805 to_all+scalar_pair, agrid)
807 to_all+scalar_pair, agrid)
809 if (cs%dynamic_psurf) &
811 if (interp_eta_pf)
then 817 if (integral_bt_cont) &
824 cs%BT_Domain, to_all+scalar_pair)
825 if (cs%linear_wave_drag) &
827 cs%BT_Domain, to_all+scalar_pair)
830 if (((g%isd > cs%isdw) .or. (g%jsd > cs%jsdw)) .or. (isq <= is-1) .or. (jsq <= js-1)) &
831 call create_group_pass(cs%pass_force_hbt0_Cor_ref, bt_force_u, bt_force_v, cs%BT_Domain)
832 if (add_uh0)
call create_group_pass(cs%pass_force_hbt0_Cor_ref, uhbt0, vhbt0, cs%BT_Domain)
833 call create_group_pass(cs%pass_force_hbt0_Cor_ref, cor_ref_u, cor_ref_v, cs%BT_Domain)
834 if (.not. use_bt_cont)
then 835 call create_group_pass(cs%pass_Dat_uv, datu, datv, cs%BT_Domain, to_all+scalar_pair)
839 if (integral_bt_cont)
then 842 if (apply_obc_open) &
856 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
862 if (cs%linearized_BT_PV)
then 864 do j=jsvf-2,jevf+1 ;
do i=isvf-2,ievf+1
868 do j=jsvf-1,jevf+1 ;
do i=isvf-2,ievf+1
869 dcor_u(i,j) = cs%D_u_Cor(i,j)
872 do j=jsvf-2,jevf+1 ;
do i=isvf-1,ievf+1
873 dcor_v(i,j) = cs%D_v_Cor(i,j)
876 q(:,:) = 0.0 ; dcor_u(:,:) = 0.0 ; dcor_v(:,:) = 0.0
877 if (gv%Boussinesq)
then 879 do j=js,je ;
do i=is-1,ie
880 dcor_u(i,j) = 0.5 * (max(gv%Z_to_H*g%bathyT(i+1,j) + eta_in(i+1,j), 0.0) + &
881 max(gv%Z_to_H*g%bathyT(i,j) + eta_in(i,j), 0.0) )
884 do j=js-1,je ;
do i=is,ie
885 dcor_v(i,j) = 0.5 * (max(gv%Z_to_H*g%bathyT(i,j+1) + eta_in(i+1,j), 0.0) + &
886 max(gv%Z_to_H*g%bathyT(i,j) + eta_in(i,j), 0.0) )
889 do j=js-1,je ;
do i=is-1,ie
890 q(i,j) = 0.25 * (cs%BT_Coriolis_scale * g%CoriolisBu(i,j)) * &
891 ((g%areaT(i,j) + g%areaT(i+1,j+1)) + (g%areaT(i+1,j) + g%areaT(i,j+1))) / &
892 (max((g%areaT(i,j) * max(gv%Z_to_H*g%bathyT(i,j) + eta_in(i,j), 0.0) + &
893 g%areaT(i+1,j+1) * max(gv%Z_to_H*g%bathyT(i+1,j+1) + eta_in(i+1,j+1), 0.0)) + &
894 (g%areaT(i+1,j) * max(gv%Z_to_H*g%bathyT(i+1,j) + eta_in(i+1,j), 0.0) + &
895 g%areaT(i,j+1) * max(gv%Z_to_H*g%bathyT(i,j+1) + eta_in(i,j+1), 0.0)), h_neglect) )
899 do j=js,je ;
do i=is-1,ie
900 dcor_u(i,j) = 0.5 * (eta_in(i+1,j) + eta_in(i,j))
903 do j=js-1,je ;
do i=is,ie
904 dcor_v(i,j) = 0.5 * (eta_in(i,j+1) + eta_in(i,j))
907 do j=js-1,je ;
do i=is-1,ie
908 q(i,j) = 0.25 * (cs%BT_Coriolis_scale * g%CoriolisBu(i,j)) * &
909 ((g%areaT(i,j) + g%areaT(i+1,j+1)) + (g%areaT(i+1,j) + g%areaT(i,j+1))) / &
910 (max((g%areaT(i,j) * eta_in(i,j) + g%areaT(i+1,j+1) * eta_in(i+1,j+1)) + &
911 (g%areaT(i+1,j) * eta_in(i+1,j) + g%areaT(i,j+1) * eta_in(i,j+1)), h_neglect) )
919 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
920 if (nonblock_setup)
then 921 call start_group_pass(cs%pass_q_DCor, cs%BT_Domain, clock=id_clock_pass_pre)
923 call do_group_pass(cs%pass_q_DCor, cs%BT_Domain, clock=id_clock_pass_pre)
925 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
930 do j=cs%jsdw,cs%jedw ;
do i=cs%isdw,cs%iedw
931 gtot_e(i,j) = 0.0 ; gtot_w(i,j) = 0.0
932 gtot_n(i,j) = 0.0 ; gtot_s(i,j) = 0.0
935 if (interp_eta_pf)
then 936 eta_pf_1(i,j) = 0.0 ; d_eta_pf(i,j) = 0.0
938 if (integral_bt_cont)
then 941 p_surf_dyn(i,j) = 0.0
942 if (cs%dynamic_psurf) dyn_coef_eta(i,j) = 0.0
948 do j=cs%jsdw,cs%jedw ;
do i=cs%isdw-1,cs%iedw
949 cor_ref_u(i,j) = 0.0 ; bt_force_u(i,j) = 0.0 ; ubt(i,j) = 0.0
950 datu(i,j) = 0.0 ; bt_rem_u(i,j) = 0.0 ; uhbt0(i,j) = 0.0
953 do j=cs%jsdw-1,cs%jedw ;
do i=cs%isdw,cs%iedw
954 cor_ref_v(i,j) = 0.0 ; bt_force_v(i,j) = 0.0 ; vbt(i,j) = 0.0
955 datv(i,j) = 0.0 ; bt_rem_v(i,j) = 0.0 ; vhbt0(i,j) = 0.0
959 if (interp_eta_pf)
then 961 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
962 eta(i,j) = eta_in(i,j)
963 eta_pf_1(i,j) = eta_pf_start(i,j)
964 d_eta_pf(i,j) = eta_pf_in(i,j) - eta_pf_start(i,j)
968 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
969 eta(i,j) = eta_in(i,j)
970 eta_pf(i,j) = eta_pf_in(i,j)
973 if (integral_bt_cont)
then 975 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
976 eta_ic(i,j) = eta_in(i,j)
981 do k=1,nz ;
do j=js,je ;
do i=is-1,ie
984 if (visc_rem_u(i,j,k) <= 0.0)
then ; visc_rem = 0.0
985 elseif (visc_rem_u(i,j,k) >= 1.0)
then ; visc_rem = 1.0
986 elseif (visc_rem_u(i,j,k)**2 > visc_rem_u(i,j,k) - 0.5*instep)
then 987 visc_rem = visc_rem_u(i,j,k)
988 else ; visc_rem = 1.0 - 0.5*instep/visc_rem_u(i,j,k) ;
endif 989 wt_u(i,j,k) = cs%frhatu(i,j,k) * visc_rem
990 enddo ;
enddo ;
enddo 992 do k=1,nz ;
do j=js-1,je ;
do i=is,ie
994 if (visc_rem_v(i,j,k) <= 0.0)
then ; visc_rem = 0.0
995 elseif (visc_rem_v(i,j,k) >= 1.0)
then ; visc_rem = 1.0
996 elseif (visc_rem_v(i,j,k)**2 > visc_rem_v(i,j,k) - 0.5*instep)
then 997 visc_rem = visc_rem_v(i,j,k)
998 else ; visc_rem = 1.0 - 0.5*instep/visc_rem_v(i,j,k) ;
endif 999 wt_v(i,j,k) = cs%frhatv(i,j,k) * visc_rem
1000 enddo ;
enddo ;
enddo 1005 do j=js-1,je+1 ;
do i=is-1,ie ; ubt_cor(i,j) = 0.0 ;
enddo ;
enddo 1007 do j=js-1,je ;
do i=is-1,ie+1 ; vbt_cor(i,j) = 0.0 ;
enddo ;
enddo 1009 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
1010 ubt_cor(i,j) = ubt_cor(i,j) + wt_u(i,j,k) * u_cor(i,j,k)
1011 enddo ;
enddo ;
enddo 1013 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
1014 vbt_cor(i,j) = vbt_cor(i,j) + wt_v(i,j,k) * v_cor(i,j,k)
1015 enddo ;
enddo ;
enddo 1023 do k=1,nz ;
do i=is-1,ie
1024 gtot_e(i,j) = gtot_e(i,j) + pbce(i,j,k) * wt_u(i,j,k)
1025 gtot_w(i+1,j) = gtot_w(i+1,j) + pbce(i+1,j,k) * wt_u(i,j,k)
1030 do k=1,nz ;
do i=is,ie
1031 gtot_n(i,j) = gtot_n(i,j) + pbce(i,j,k) * wt_v(i,j,k)
1032 gtot_s(i,j+1) = gtot_s(i,j+1) + pbce(i,j+1,k) * wt_v(i,j,k)
1037 call tidal_forcing_sensitivity(g, cs%tides_CSp, det_de)
1038 dgeo_de = 1.0 + det_de + cs%G_extra
1040 dgeo_de = 1.0 + cs%G_extra
1043 if (nonblock_setup .and. .not.cs%linearized_BT_PV)
then 1044 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1045 call complete_group_pass(cs%pass_q_DCor, cs%BT_Domain, clock=id_clock_pass_pre)
1046 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1051 if (integral_bt_cont)
then 1052 call set_local_bt_cont_types(bt_cont, btcl_u, btcl_v, g, us, ms, cs%BT_Domain, 1+ievf-ie, dt_baroclinic=dt)
1053 elseif (use_bt_cont)
then 1054 call set_local_bt_cont_types(bt_cont, btcl_u, btcl_v, g, us, ms, cs%BT_Domain, 1+ievf-ie)
1056 if (cs%Nonlinear_continuity)
then 1057 call find_face_areas(datu, datv, g, gv, us, cs, ms, eta, 1)
1059 call find_face_areas(datu, datv, g, gv, us, cs, ms, halo=1)
1064 if (apply_obcs)
then 1065 call set_up_bt_obc(obc, eta, cs%BT_OBC, cs%BT_Domain, g, gv, us, ms, ievf-ie, use_bt_cont, &
1066 integral_bt_cont, dt, datu, datv, btcl_u, btcl_v)
1073 do j=js,je ;
do i=is-1,ie ; uhbt(i,j) = 0.0 ; ubt(i,j) = 0.0 ;
enddo ;
enddo 1075 do j=js-1,je ;
do i=is,ie ; vhbt(i,j) = 0.0 ; vbt(i,j) = 0.0 ;
enddo ;
enddo 1076 if (cs%visc_rem_u_uh0)
then 1078 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
1079 uhbt(i,j) = uhbt(i,j) + uh0(i,j,k)
1080 ubt(i,j) = ubt(i,j) + wt_u(i,j,k) * u_uh0(i,j,k)
1081 enddo ;
enddo ;
enddo 1083 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
1084 vhbt(i,j) = vhbt(i,j) + vh0(i,j,k)
1085 vbt(i,j) = vbt(i,j) + wt_v(i,j,k) * v_vh0(i,j,k)
1086 enddo ;
enddo ;
enddo 1089 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
1090 uhbt(i,j) = uhbt(i,j) + uh0(i,j,k)
1091 ubt(i,j) = ubt(i,j) + cs%frhatu(i,j,k) * u_uh0(i,j,k)
1092 enddo ;
enddo ;
enddo 1094 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
1095 vhbt(i,j) = vhbt(i,j) + vh0(i,j,k)
1096 vbt(i,j) = vbt(i,j) + cs%frhatv(i,j,k) * v_vh0(i,j,k)
1097 enddo ;
enddo ;
enddo 1099 if ((use_bt_cont .or. integral_bt_cont) .and. cs%adjust_BT_cont)
then 1104 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1105 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
1106 call pass_vector(ubt, vbt, cs%BT_Domain, complete=.false., halo=1+ievf-ie)
1107 call pass_vector(uhbt, vhbt, cs%BT_Domain, complete=.true., halo=1+ievf-ie)
1108 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
1109 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1111 if (integral_bt_cont)
then 1112 call adjust_local_bt_cont_types(ubt, uhbt, vbt, vhbt, btcl_u, btcl_v, &
1113 g, us, ms, halo=1+ievf-ie, dt_baroclinic=dt)
1115 call adjust_local_bt_cont_types(ubt, uhbt, vbt, vhbt, btcl_u, btcl_v, &
1116 g, us, ms, halo=1+ievf-ie)
1119 if (integral_bt_cont)
then 1121 do j=js,je ;
do i=is-1,ie
1122 uhbt0(i,j) = uhbt(i,j) - find_uhbt(dt*ubt(i,j), btcl_u(i,j)) * idt
1125 do j=js-1,je ;
do i=is,ie
1126 vhbt0(i,j) = vhbt(i,j) - find_vhbt(dt*vbt(i,j), btcl_v(i,j)) * idt
1128 elseif (use_bt_cont)
then 1130 do j=js,je ;
do i=is-1,ie
1131 uhbt0(i,j) = uhbt(i,j) - find_uhbt(ubt(i,j), btcl_u(i,j))
1134 do j=js-1,je ;
do i=is,ie
1135 vhbt0(i,j) = vhbt(i,j) - find_vhbt(vbt(i,j), btcl_v(i,j))
1139 do j=js,je ;
do i=is-1,ie
1140 uhbt0(i,j) = uhbt(i,j) - datu(i,j)*ubt(i,j)
1143 do j=js-1,je ;
do i=is,ie
1144 vhbt0(i,j) = vhbt(i,j) - datv(i,j)*vbt(i,j)
1147 if (cs%BT_OBC%apply_u_OBCs)
then 1149 do j=js,je ;
do i=is-1,ie ;
if (obc%segnum_u(i,j) /= obc_none)
then 1151 endif ;
enddo ;
enddo 1153 if (cs%BT_OBC%apply_v_OBCs)
then 1155 do j=js-1,je ;
do i=is,ie ;
if (obc%segnum_v(i,j) /= obc_none)
then 1157 endif ;
enddo ;
enddo 1162 if (integral_bt_cont)
then 1164 do j=jsvf-1,jevf+1 ;
do i=isvf-2,ievf+1
1165 ubt(i,j) = 0.0 ; uhbt(i,j) = 0.0 ; u_accel_bt(i,j) = 0.0
1166 ubt_int(i,j) = 0.0 ; uhbt_int(i,j) = 0.0
1169 do j=jsvf-2,jevf+1 ;
do i=isvf-1,ievf+1
1170 vbt(i,j) = 0.0 ; vhbt(i,j) = 0.0 ; v_accel_bt(i,j) = 0.0
1171 vbt_int(i,j) = 0.0 ; vhbt_int(i,j) = 0.0
1175 do j=jsvf-1,jevf+1 ;
do i=isvf-2,ievf+1
1176 ubt(i,j) = 0.0 ; uhbt(i,j) = 0.0 ; u_accel_bt(i,j) = 0.0
1179 do j=jsvf-2,jevf+1 ;
do i=isvf-1,ievf+1
1180 vbt(i,j) = 0.0 ; vhbt(i,j) = 0.0 ; v_accel_bt(i,j) = 0.0
1184 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
1185 ubt(i,j) = ubt(i,j) + wt_u(i,j,k) * u_in(i,j,k)
1186 enddo ;
enddo ;
enddo 1188 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
1189 vbt(i,j) = vbt(i,j) + wt_v(i,j,k) * v_in(i,j,k)
1190 enddo ;
enddo ;
enddo 1192 do j=js,je ;
do i=is-1,ie
1193 if (abs(ubt(i,j)) < cs%vel_underflow) ubt(i,j) = 0.0
1196 do j=js-1,je ;
do i=is,ie
1197 if (abs(vbt(i,j)) < cs%vel_underflow) vbt(i,j) = 0.0
1200 if (apply_obcs)
then 1201 ubt_first(:,:) = ubt(:,:) ; vbt_first(:,:) = vbt(:,:)
1211 do j=js,je ;
do i=is-1,ie ;
if (g%mask2dCu(i,j) > 0.0)
then 1212 if (cs%nonlin_stress)
then 1213 if (gv%Boussinesq)
then 1214 htot_avg = 0.5*(max(cs%bathyT(i,j)*gv%Z_to_H + eta(i,j), 0.0) + &
1215 max(cs%bathyT(i+1,j)*gv%Z_to_H + eta(i+1,j), 0.0))
1217 htot_avg = 0.5*(eta(i,j) + eta(i+1,j))
1219 if (htot_avg*cs%dy_Cu(i,j) <= 0.0)
then 1221 elseif (integral_bt_cont)
then 1222 cs%IDatu(i,j) = cs%dy_Cu(i,j) / (max(find_duhbt_du(ubt(i,j)*dt, btcl_u(i,j)), &
1223 cs%dy_Cu(i,j)*htot_avg) )
1224 elseif (use_bt_cont)
then 1225 cs%IDatu(i,j) = cs%dy_Cu(i,j) / (max(find_duhbt_du(ubt(i,j), btcl_u(i,j)), &
1226 cs%dy_Cu(i,j)*htot_avg) )
1228 cs%IDatu(i,j) = 1.0 / htot_avg
1232 bt_force_u(i,j) = forces%taux(i,j) * mass_accel_to_z * cs%IDatu(i,j)*visc_rem_u(i,j,1)
1234 bt_force_u(i,j) = 0.0
1235 endif ;
enddo ;
enddo 1237 do j=js-1,je ;
do i=is,ie ;
if (g%mask2dCv(i,j) > 0.0)
then 1238 if (cs%nonlin_stress)
then 1239 if (gv%Boussinesq)
then 1240 htot_avg = 0.5*(max(cs%bathyT(i,j)*gv%Z_to_H + eta(i,j), 0.0) + &
1241 max(cs%bathyT(i,j+1)*gv%Z_to_H + eta(i,j+1), 0.0))
1243 htot_avg = 0.5*(eta(i,j) + eta(i,j+1))
1245 if (htot_avg*cs%dx_Cv(i,j) <= 0.0)
then 1247 elseif (integral_bt_cont)
then 1248 cs%IDatv(i,j) = cs%dx_Cv(i,j) / (max(find_dvhbt_dv(vbt(i,j)*dt, btcl_v(i,j)), &
1249 cs%dx_Cv(i,j)*htot_avg) )
1250 elseif (use_bt_cont)
then 1251 cs%IDatv(i,j) = cs%dx_Cv(i,j) / (max(find_dvhbt_dv(vbt(i,j), btcl_v(i,j)), &
1252 cs%dx_Cv(i,j)*htot_avg) )
1254 cs%IDatv(i,j) = 1.0 / htot_avg
1258 bt_force_v(i,j) = forces%tauy(i,j) * mass_accel_to_z * cs%IDatv(i,j)*visc_rem_v(i,j,1)
1260 bt_force_v(i,j) = 0.0
1261 endif ;
enddo ;
enddo 1262 if (
present(taux_bot) .and.
present(tauy_bot))
then 1263 if (
associated(taux_bot) .and.
associated(tauy_bot))
then 1265 do j=js,je ;
do i=is-1,ie ;
if (g%mask2dCu(i,j) > 0.0)
then 1266 bt_force_u(i,j) = bt_force_u(i,j) - taux_bot(i,j) * mass_to_z * cs%IDatu(i,j)
1267 endif ;
enddo ;
enddo 1269 do j=js-1,je ;
do i=is,ie ;
if (g%mask2dCv(i,j) > 0.0)
then 1270 bt_force_v(i,j) = bt_force_v(i,j) - tauy_bot(i,j) * mass_to_z * cs%IDatv(i,j)
1271 endif ;
enddo ;
enddo 1278 do j=js,je ;
do k=1,nz ;
do i=isq,ieq
1279 bt_force_u(i,j) = bt_force_u(i,j) + wt_u(i,j,k) * bc_accel_u(i,j,k)
1280 enddo ;
enddo ;
enddo 1282 do j=jsq,jeq ;
do k=1,nz ;
do i=is,ie
1283 bt_force_v(i,j) = bt_force_v(i,j) + wt_v(i,j,k) * bc_accel_v(i,j,k)
1284 enddo ;
enddo ;
enddo 1286 if (cs%gradual_BT_ICs)
then 1288 do j=js,je ;
do i=is-1,ie
1289 bt_force_u(i,j) = bt_force_u(i,j) + (ubt(i,j) - cs%ubt_IC(i,j)) * idt
1290 ubt(i,j) = cs%ubt_IC(i,j)
1291 if (abs(ubt(i,j)) < cs%vel_underflow) ubt(i,j) = 0.0
1294 do j=js-1,je ;
do i=is,ie
1295 bt_force_v(i,j) = bt_force_v(i,j) + (vbt(i,j) - cs%vbt_IC(i,j)) * idt
1296 vbt(i,j) = cs%vbt_IC(i,j)
1297 if (abs(vbt(i,j)) < cs%vel_underflow) vbt(i,j) = 0.0
1301 if ((isq > is-1) .or. (jsq > js-1))
then 1304 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1305 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
1306 tmp_u(:,:) = 0.0 ; tmp_v(:,:) = 0.0
1307 do j=js,je ;
do i=isq,ieq ; tmp_u(i,j) = bt_force_u(i,j) ;
enddo ;
enddo 1308 do j=jsq,jeq ;
do i=is,ie ; tmp_v(i,j) = bt_force_v(i,j) ;
enddo ;
enddo 1309 if (nonblock_setup)
then 1310 call start_group_pass(cs%pass_tmp_uv, g%Domain)
1312 call do_group_pass(cs%pass_tmp_uv, g%Domain)
1313 do j=jsd,jed ;
do i=isdb,iedb ; bt_force_u(i,j) = tmp_u(i,j) ;
enddo ;
enddo 1314 do j=jsdb,jedb ;
do i=isd,ied ; bt_force_v(i,j) = tmp_v(i,j) ;
enddo ;
enddo 1316 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
1317 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1320 if (nonblock_setup)
then 1321 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1322 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
1323 call start_group_pass(cs%pass_gtot, cs%BT_Domain)
1324 call start_group_pass(cs%pass_ubt_Cor, g%Domain)
1325 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
1326 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1331 do j=jsvf-1,jevf ;
do i=isvf-1,ievf+1
1332 if (cs%Sadourny)
then 1333 amer(i-1,j) = dcor_u(i-1,j) * q(i-1,j)
1334 bmer(i,j) = dcor_u(i,j) * q(i,j)
1335 cmer(i,j+1) = dcor_u(i,j+1) * q(i,j)
1336 dmer(i-1,j+1) = dcor_u(i-1,j+1) * q(i-1,j)
1338 amer(i-1,j) = dcor_u(i-1,j) * &
1339 ((q(i,j) + q(i-1,j-1)) + q(i-1,j)) / 3.0
1340 bmer(i,j) = dcor_u(i,j) * &
1341 (q(i,j) + (q(i-1,j) + q(i,j-1))) / 3.0
1342 cmer(i,j+1) = dcor_u(i,j+1) * &
1343 (q(i,j) + (q(i-1,j) + q(i,j+1))) / 3.0
1344 dmer(i-1,j+1) = dcor_u(i-1,j+1) * &
1345 ((q(i,j) + q(i-1,j+1)) + q(i-1,j)) / 3.0
1350 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf
1351 if (cs%Sadourny)
then 1352 azon(i,j) = dcor_v(i+1,j) * q(i,j)
1353 bzon(i,j) = dcor_v(i,j) * q(i,j)
1354 czon(i,j) = dcor_v(i,j-1) * q(i,j-1)
1355 dzon(i,j) = dcor_v(i+1,j-1) * q(i,j-1)
1357 azon(i,j) = dcor_v(i+1,j) * &
1358 (q(i,j) + (q(i+1,j) + q(i,j-1))) / 3.0
1359 bzon(i,j) = dcor_v(i,j) * &
1360 (q(i,j) + (q(i-1,j) + q(i,j-1))) / 3.0
1361 czon(i,j) = dcor_v(i,j-1) * &
1362 ((q(i,j) + q(i-1,j-1)) + q(i,j-1)) / 3.0
1363 dzon(i,j) = dcor_v(i+1,j-1) * &
1364 ((q(i,j) + q(i+1,j-1)) + q(i,j-1)) / 3.0
1369 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1370 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
1371 if (nonblock_setup)
then 1372 if ((isq > is-1) .or. (jsq > js-1))
then 1373 call complete_group_pass(cs%pass_tmp_uv, g%Domain)
1374 do j=jsd,jed ;
do i=isdb,iedb ; bt_force_u(i,j) = tmp_u(i,j) ;
enddo ;
enddo 1375 do j=jsdb,jedb ;
do i=isd,ied ; bt_force_v(i,j) = tmp_v(i,j) ;
enddo ;
enddo 1377 call complete_group_pass(cs%pass_gtot, cs%BT_Domain)
1378 call complete_group_pass(cs%pass_ubt_Cor, g%Domain)
1380 call do_group_pass(cs%pass_gtot, cs%BT_Domain)
1381 call do_group_pass(cs%pass_ubt_Cor, g%Domain)
1385 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf+1
1386 if (cs%ua_polarity(i,j) < 0.0)
call swap(gtot_e(i,j), gtot_w(i,j))
1387 if (cs%va_polarity(i,j) < 0.0)
call swap(gtot_n(i,j), gtot_s(i,j))
1391 do j=js,je ;
do i=is-1,ie
1393 ((azon(i,j) * vbt_cor(i+1,j) + czon(i,j) * vbt_cor(i ,j-1)) + &
1394 (bzon(i,j) * vbt_cor(i ,j) + dzon(i,j) * vbt_cor(i+1,j-1)))
1397 do j=js-1,je ;
do i=is,ie
1398 cor_ref_v(i,j) = -1.0 * &
1399 ((amer(i-1,j) * ubt_cor(i-1,j) + cmer(i ,j+1) * ubt_cor(i ,j+1)) + &
1400 (bmer(i ,j) * ubt_cor(i ,j) + dmer(i-1,j+1) * ubt_cor(i-1,j+1)))
1404 if (nonblock_setup)
then 1405 if (.not.use_bt_cont) &
1406 call start_group_pass(cs%pass_Dat_uv, cs%BT_Domain)
1409 call start_group_pass(cs%pass_force_hbt0_Cor_ref, cs%BT_Domain)
1411 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
1412 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1415 do j=js-1,je+1 ;
do i=is-1,ie ; av_rem_u(i,j) = 0.0 ;
enddo ;
enddo 1417 do j=js-1,je ;
do i=is-1,ie+1 ; av_rem_v(i,j) = 0.0 ;
enddo ;
enddo 1419 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
1420 av_rem_u(i,j) = av_rem_u(i,j) + cs%frhatu(i,j,k) * visc_rem_u(i,j,k)
1421 enddo ;
enddo ;
enddo 1423 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
1424 av_rem_v(i,j) = av_rem_v(i,j) + cs%frhatv(i,j,k) * visc_rem_v(i,j,k)
1425 enddo ;
enddo ;
enddo 1426 if (cs%strong_drag)
then 1428 do j=js,je ;
do i=is-1,ie
1429 bt_rem_u(i,j) = g%mask2dCu(i,j) * &
1430 ((nstep * av_rem_u(i,j)) / (1.0 + (nstep-1)*av_rem_u(i,j)))
1433 do j=js-1,je ;
do i=is,ie
1434 bt_rem_v(i,j) = g%mask2dCv(i,j) * &
1435 ((nstep * av_rem_v(i,j)) / (1.0 + (nstep-1)*av_rem_v(i,j)))
1439 do j=js,je ;
do i=is-1,ie
1441 if (g%mask2dCu(i,j) * av_rem_u(i,j) > 0.0) &
1442 bt_rem_u(i,j) = g%mask2dCu(i,j) * (av_rem_u(i,j)**instep)
1445 do j=js-1,je ;
do i=is,ie
1447 if (g%mask2dCv(i,j) * av_rem_v(i,j) > 0.0) &
1448 bt_rem_v(i,j) = g%mask2dCv(i,j) * (av_rem_v(i,j)**instep)
1451 if (cs%linear_wave_drag)
then 1453 do j=js,je ;
do i=is-1,ie ;
if (cs%lin_drag_u(i,j) > 0.0)
then 1454 htot = 0.5 * (eta(i,j) + eta(i+1,j))
1455 if (gv%Boussinesq) &
1456 htot = htot + 0.5*gv%Z_to_H * (cs%bathyT(i,j) + cs%bathyT(i+1,j))
1457 bt_rem_u(i,j) = bt_rem_u(i,j) * (htot / (htot + cs%lin_drag_u(i,j) * dtbt))
1459 rayleigh_u(i,j) = cs%lin_drag_u(i,j) / htot
1460 endif ;
enddo ;
enddo 1462 do j=js-1,je ;
do i=is,ie ;
if (cs%lin_drag_v(i,j) > 0.0)
then 1463 htot = 0.5 * (eta(i,j) + eta(i,j+1))
1464 if (gv%Boussinesq) &
1465 htot = htot + 0.5*gv%Z_to_H * (cs%bathyT(i,j) + cs%bathyT(i+1,j+1))
1466 bt_rem_v(i,j) = bt_rem_v(i,j) * (htot / (htot + cs%lin_drag_v(i,j) * dtbt))
1468 rayleigh_v(i,j) = cs%lin_drag_v(i,j) / htot
1469 endif ;
enddo ;
enddo 1473 if (find_etaav)
then 1475 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf+1
1476 eta_sum(i,j) = 0.0 ; eta_wtd(i,j) = 0.0
1480 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf+1
1485 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf
1486 ubt_sum(i,j) = 0.0 ; uhbt_sum(i,j) = 0.0
1487 pfu_bt_sum(i,j) = 0.0 ; coru_bt_sum(i,j) = 0.0
1488 ubt_wtd(i,j) = 0.0 ; ubt_trans(i,j) = 0.0
1491 do j=jsvf-1,jevf ;
do i=isvf-1,ievf+1
1492 vbt_sum(i,j) = 0.0 ; vhbt_sum(i,j) = 0.0
1493 pfv_bt_sum(i,j) = 0.0 ; corv_bt_sum(i,j) = 0.0
1494 vbt_wtd(i,j) = 0.0 ; vbt_trans(i,j) = 0.0
1499 do j=jsvf-1,jevf+1;
do i=isvf-1,ievf+1 ; eta_src(i,j) = 0.0 ;
enddo ;
enddo 1500 if (cs%bound_BT_corr)
then ;
if ((use_bt_cont.or.integral_bt_cont) .and. cs%BT_cont_bounds)
then 1501 do j=js,je ;
do i=is,ie ;
if (g%mask2dT(i,j) > 0.0)
then 1502 if (cs%eta_cor(i,j) > 0.0)
then 1505 if (integral_bt_cont)
then 1506 uint_cor = g%dxT(i,j) * cs%maxCFL_BT_cont
1507 vint_cor = g%dyT(i,j) * cs%maxCFL_BT_cont
1508 eta_cor_max = (cs%IareaT(i,j) * &
1509 (((find_uhbt(uint_cor, btcl_u(i,j)) + dt*uhbt0(i,j)) - &
1510 (find_uhbt(-uint_cor, btcl_u(i-1,j)) + dt*uhbt0(i-1,j))) + &
1511 ((find_vhbt(vint_cor, btcl_v(i,j)) + dt*vhbt0(i,j)) - &
1512 (find_vhbt(-vint_cor, btcl_v(i,j-1)) + dt*vhbt0(i,j-1))) ))
1514 u_max_cor = g%dxT(i,j) * (cs%maxCFL_BT_cont*idt)
1515 v_max_cor = g%dyT(i,j) * (cs%maxCFL_BT_cont*idt)
1516 eta_cor_max = dt * (cs%IareaT(i,j) * &
1517 (((find_uhbt(u_max_cor, btcl_u(i,j)) + uhbt0(i,j)) - &
1518 (find_uhbt(-u_max_cor, btcl_u(i-1,j)) + uhbt0(i-1,j))) + &
1519 ((find_vhbt(v_max_cor, btcl_v(i,j)) + vhbt0(i,j)) - &
1520 (find_vhbt(-v_max_cor, btcl_v(i,j-1)) + vhbt0(i,j-1))) ))
1522 cs%eta_cor(i,j) = min(cs%eta_cor(i,j), max(0.0, eta_cor_max))
1526 if (gv%Boussinesq) htot = cs%bathyT(i,j)*gv%Z_to_H + eta(i,j)
1528 cs%eta_cor(i,j) = max(cs%eta_cor(i,j), -max(0.0,htot))
1530 endif ;
enddo ;
enddo 1531 else ;
do j=js,je ;
do i=is,ie
1532 if (abs(cs%eta_cor(i,j)) > dt*cs%eta_cor_bound(i,j)) &
1533 cs%eta_cor(i,j) = sign(dt*cs%eta_cor_bound(i,j), cs%eta_cor(i,j))
1534 enddo ;
enddo ;
endif ;
endif 1536 do j=js,je ;
do i=is,ie
1537 eta_src(i,j) = g%mask2dT(i,j) * (instep * cs%eta_cor(i,j))
1541 if (cs%dynamic_psurf)
then 1542 ice_is_rigid = (
associated(forces%rigidity_ice_u) .and. &
1543 associated(forces%rigidity_ice_v))
1544 h_min_dyn = gv%Z_to_H * cs%Dmin_dyn_psurf
1545 if (ice_is_rigid .and. use_bt_cont) &
1546 call bt_cont_to_face_areas(bt_cont, datu, datv, g, us, ms, 0, .true.)
1547 if (ice_is_rigid)
then 1549 do j=js,je ;
do i=is,ie
1555 idt_max2 = 0.5 * (dgeo_de * (1.0 + 2.0*bebt)) * (g%IareaT(i,j) * &
1556 ((gtot_e(i,j) * (datu(i,j)*g%IdxCu(i,j)) + &
1557 gtot_w(i,j) * (datu(i-1,j)*g%IdxCu(i-1,j))) + &
1558 (gtot_n(i,j) * (datv(i,j)*g%IdyCv(i,j)) + &
1559 gtot_s(i,j) * (datv(i,j-1)*g%IdyCv(i,j-1)))) + &
1560 ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
1561 (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2)) * cs%BT_Coriolis_scale**2 )
1562 h_eff_dx2 = max(h_min_dyn * ((g%IdxT(i,j))**2 + (g%IdyT(i,j))**2), &
1564 ((datu(i,j)*g%IdxCu(i,j) + datu(i-1,j)*g%IdxCu(i-1,j)) + &
1565 (datv(i,j)*g%IdyCv(i,j) + datv(i,j-1)*g%IdyCv(i,j-1)) ) )
1566 dyn_coef_max = cs%const_dyn_psurf * max(0.0, 1.0 - dtbt**2 * idt_max2) / &
1567 (dtbt**2 * h_eff_dx2)
1570 ice_strength = ((forces%rigidity_ice_u(i,j) + forces%rigidity_ice_u(i-1,j)) + &
1571 (forces%rigidity_ice_v(i,j) + forces%rigidity_ice_v(i,j-1))) / &
1572 (cs%ice_strength_length**2 * dtbt)
1575 dyn_coef_eta(i,j) = min(dyn_coef_max, ice_strength * gv%H_to_Z)
1576 enddo ;
enddo ;
endif 1579 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1580 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
1581 if (nonblock_setup)
then 1582 call start_group_pass(cs%pass_eta_bt_rem, cs%BT_Domain)
1585 call do_group_pass(cs%pass_eta_bt_rem, cs%BT_Domain)
1586 if (.not.use_bt_cont) &
1587 call do_group_pass(cs%pass_Dat_uv, cs%BT_Domain)
1588 call do_group_pass(cs%pass_force_hbt0_Cor_ref, cs%BT_Domain)
1590 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
1591 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1594 if (nonblock_setup)
then 1595 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1596 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
1598 if (.not.use_bt_cont)
call complete_group_pass(cs%pass_Dat_uv, cs%BT_Domain)
1599 call complete_group_pass(cs%pass_force_hbt0_Cor_ref, cs%BT_Domain)
1600 call complete_group_pass(cs%pass_eta_bt_rem, cs%BT_Domain)
1602 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
1603 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1607 call uvchksum(
"BT [uv]hbt", uhbt, vhbt, cs%debug_BT_HI, haloshift=0, &
1608 scale=us%s_to_T*us%L_to_m**2*gv%H_to_m)
1609 call uvchksum(
"BT Initial [uv]bt", ubt, vbt, cs%debug_BT_HI, haloshift=0, scale=us%L_T_to_m_s)
1610 call hchksum(eta,
"BT Initial eta", cs%debug_BT_HI, haloshift=0, scale=gv%H_to_m)
1611 call uvchksum(
"BT BT_force_[uv]", bt_force_u, bt_force_v, &
1612 cs%debug_BT_HI, haloshift=0, scale=us%L_T2_to_m_s2)
1613 if (interp_eta_pf)
then 1614 call hchksum(eta_pf_1,
"BT eta_PF_1",cs%debug_BT_HI,haloshift=0, scale=gv%H_to_m)
1615 call hchksum(d_eta_pf,
"BT d_eta_PF",cs%debug_BT_HI,haloshift=0, scale=gv%H_to_m)
1617 call hchksum(eta_pf,
"BT eta_PF",cs%debug_BT_HI,haloshift=0, scale=gv%H_to_m)
1618 call hchksum(eta_pf_in,
"BT eta_PF_in",g%HI,haloshift=0, scale=gv%H_to_m)
1620 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)
1621 call uvchksum(
"BT [uv]hbt0", uhbt0, vhbt0, cs%debug_BT_HI, haloshift=0, &
1622 scale=us%L_to_m**2*us%s_to_T*gv%H_to_m)
1623 if (.not. use_bt_cont)
then 1624 call uvchksum(
"BT Dat[uv]", datu, datv, cs%debug_BT_HI, haloshift=1, scale=us%L_to_m*gv%H_to_m)
1626 call uvchksum(
"BT wt_[uv]", wt_u, wt_v, g%HI, haloshift=0, &
1627 symmetric=.true., omit_corners=.true., scalar_pair=.true.)
1628 call uvchksum(
"BT frhat[uv]", cs%frhatu, cs%frhatv, g%HI, haloshift=0, &
1629 symmetric=.true., omit_corners=.true., scalar_pair=.true.)
1630 call uvchksum(
"BT bc_accel_[uv]", bc_accel_u, bc_accel_v, g%HI, haloshift=0, scale=us%L_T2_to_m_s2)
1631 call uvchksum(
"BT IDat[uv]", cs%IDatu, cs%IDatv, g%HI, haloshift=0, &
1632 scale=us%m_to_Z, scalar_pair=.true.)
1633 call uvchksum(
"BT visc_rem_[uv]", visc_rem_u, visc_rem_v, g%HI, &
1634 haloshift=1, scalar_pair=.true.)
1637 if (cs%id_ubtdt > 0)
then 1638 do j=js-1,je+1 ;
do i=is-1,ie
1639 ubt_st(i,j) = ubt(i,j)
1642 if (cs%id_vbtdt > 0)
then 1643 do j=js-1,je ;
do i=is-1,ie+1
1644 vbt_st(i,j) = vbt(i,j)
1648 if (query_averaging_enabled(cs%diag))
then 1649 if (cs%id_eta_st > 0)
call post_data(cs%id_eta_st, eta(isd:ied,jsd:jed), cs%diag)
1650 if (cs%id_ubt_st > 0)
call post_data(cs%id_ubt_st, ubt(isdb:iedb,jsd:jed), cs%diag)
1651 if (cs%id_vbt_st > 0)
call post_data(cs%id_vbt_st, vbt(isd:ied,jsdb:jedb), cs%diag)
1654 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1655 if (id_clock_calc > 0)
call cpu_clock_begin(id_clock_calc)
1657 if (project_velocity)
then ; eta_pf_bt => eta ;
else ; eta_pf_bt => eta_pred ;
endif 1659 if (cs%dt_bt_filter >= 0.0)
then 1660 dt_filt = 0.5 * max(0.0, min(cs%dt_bt_filter, 2.0*dt))
1662 dt_filt = 0.5 * max(0.0, dt * min(-cs%dt_bt_filter, 2.0))
1664 nfilter = ceiling(dt_filt / dtbt)
1666 if (nstep+nfilter==0 )
call mom_error(fatal, &
1667 "btstep: number of barotropic step (nstep+nfilter) is 0")
1670 sum_wt_vel = 0.0 ; sum_wt_eta = 0.0 ; sum_wt_accel = 0.0 ; sum_wt_trans = 0.0
1671 allocate(wt_vel(nstep+nfilter)) ;
allocate(wt_eta(nstep+nfilter))
1672 allocate(wt_trans(nstep+nfilter+1)) ;
allocate(wt_accel(nstep+nfilter+1))
1673 allocate(wt_accel2(nstep+nfilter+1))
1674 do n=1,nstep+nfilter
1678 if ( (n==nstep) .or. (dt_filt - abs(n-nstep)*dtbt >= 0.0))
then 1679 wt_vel(n) = 1.0 ; wt_eta(n) = 1.0
1680 elseif (dtbt + dt_filt - abs(n-nstep)*dtbt > 0.0)
then 1681 wt_vel(n) = 1.0 + (dt_filt / dtbt) - abs(n-nstep) ; wt_eta(n) = wt_vel(n)
1683 wt_vel(n) = 0.0 ; wt_eta(n) = 0.0
1690 sum_wt_vel = sum_wt_vel + wt_vel(n) ; sum_wt_eta = sum_wt_eta + wt_eta(n)
1692 wt_trans(nstep+nfilter+1) = 0.0 ; wt_accel(nstep+nfilter+1) = 0.0
1693 do n=nstep+nfilter,1,-1
1694 wt_trans(n) = wt_trans(n+1) + wt_eta(n)
1695 wt_accel(n) = wt_accel(n+1) + wt_vel(n)
1696 sum_wt_accel = sum_wt_accel + wt_accel(n) ; sum_wt_trans = sum_wt_trans + wt_trans(n)
1699 i_sum_wt_vel = 1.0 / sum_wt_vel ; i_sum_wt_accel = 1.0 / sum_wt_accel
1700 i_sum_wt_eta = 1.0 / sum_wt_eta ; i_sum_wt_trans = 1.0 / sum_wt_trans
1701 do n=1,nstep+nfilter
1702 wt_vel(n) = wt_vel(n) * i_sum_wt_vel
1703 if (cs%answers_2018)
then 1704 wt_accel2(n) = wt_accel(n)
1707 wt_accel2(n) = wt_accel(n) * i_sum_wt_accel
1708 wt_trans(n) = wt_trans(n) * i_sum_wt_trans
1710 wt_accel(n) = wt_accel(n) * i_sum_wt_accel
1711 wt_eta(n) = wt_eta(n) * i_sum_wt_eta
1714 sum_wt_vel = 0.0 ; sum_wt_eta = 0.0 ; sum_wt_accel = 0.0 ; sum_wt_trans = 0.0
1717 isv=is ; iev=ie ; jsv=js ; jev=je
1718 do n=1,nstep+nfilter
1720 sum_wt_vel = sum_wt_vel + wt_vel(n)
1721 sum_wt_eta = sum_wt_eta + wt_eta(n)
1722 sum_wt_accel = sum_wt_accel + wt_accel2(n)
1723 sum_wt_trans = sum_wt_trans + wt_trans(n)
1725 if (cs%clip_velocity)
then 1726 do j=jsv,jev ;
do i=isv-1,iev
1727 if ((ubt(i,j) * (dt * g%dy_Cu(i,j))) * g%IareaT(i+1,j) < -cs%CFL_trunc)
then 1729 ubt(i,j) = (-0.95*cs%CFL_trunc) * (g%areaT(i+1,j) / (dt * g%dy_Cu(i,j)))
1730 elseif ((ubt(i,j) * (dt * g%dy_Cu(i,j))) * g%IareaT(i,j) > cs%CFL_trunc)
then 1732 ubt(i,j) = (0.95*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dy_Cu(i,j)))
1735 do j=jsv-1,jev ;
do i=isv,iev
1736 if ((vbt(i,j) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j+1) < -cs%CFL_trunc)
then 1738 vbt(i,j) = (-0.9*cs%CFL_trunc) * (g%areaT(i,j+1) / (dt * g%dx_Cv(i,j)))
1739 elseif ((vbt(i,j) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j) > cs%CFL_trunc)
then 1741 vbt(i,j) = (0.9*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dx_Cv(i,j)))
1746 if ((iev - stencil < ie) .or. (jev - stencil < je))
then 1747 if (id_clock_calc > 0)
call cpu_clock_end(id_clock_calc)
1748 call do_group_pass(cs%pass_eta_ubt, cs%BT_Domain, clock=id_clock_pass_step)
1749 isv = isvf ; iev = ievf ; jsv = jsvf ; jev = jevf
1750 if (id_clock_calc > 0)
call cpu_clock_begin(id_clock_calc)
1752 isv = isv+stencil ; iev = iev-stencil
1753 jsv = jsv+stencil ; jev = jev-stencil
1756 if ((.not.use_bt_cont) .and. cs%Nonlinear_continuity .and. &
1757 (cs%Nonlin_cont_update_period > 0))
then 1758 if ((n>1) .and. (mod(n-1,cs%Nonlin_cont_update_period) == 0)) &
1759 call find_face_areas(datu, datv, g, gv, us, cs, ms, eta, 1+iev-ie)
1762 if (integral_bt_cont)
then 1764 do j=jsv-1,jev+1 ;
do i=isv-2,iev+1
1765 ubt_int_prev(i,j) = ubt_int(i,j) ; uhbt_int_prev(i,j) = uhbt_int(i,j)
1768 do j=jsv-2,jev+1 ;
do i=isv-1,iev+1
1769 vbt_int_prev(i,j) = vbt_int(i,j) ; vhbt_int_prev(i,j) = vhbt_int(i,j)
1774 if (cs%dynamic_psurf .or. .not.project_velocity)
then 1775 if (integral_bt_cont)
then 1777 do j=jsv-1,jev+1 ;
do i=isv-2,iev+1
1778 uhbt_int(i,j) = find_uhbt(ubt_int(i,j) + dtbt*ubt(i,j), btcl_u(i,j)) + n*dtbt*uhbt0(i,j)
1782 do j=jsv-2,jev+1 ;
do i=isv-1,iev+1
1783 vhbt_int(i,j) = find_vhbt(vbt_int(i,j) + dtbt*vbt(i,j), btcl_v(i,j)) + n*dtbt*vhbt0(i,j)
1786 do j=jsv-1,jev+1 ;
do i=isv-1,iev+1
1787 eta_pred(i,j) = (eta_ic(i,j) + n*eta_src(i,j)) + cs%IareaT(i,j) * &
1788 ((uhbt_int(i-1,j) - uhbt_int(i,j)) + (vhbt_int(i,j-1) - vhbt_int(i,j)))
1790 elseif (use_bt_cont)
then 1792 do j=jsv-1,jev+1 ;
do i=isv-2,iev+1
1793 uhbt(i,j) = find_uhbt(ubt(i,j), btcl_u(i,j)) + uhbt0(i,j)
1796 do j=jsv-2,jev+1 ;
do i=isv-1,iev+1
1797 vhbt(i,j) = find_vhbt(vbt(i,j), btcl_v(i,j)) + vhbt0(i,j)
1800 do j=jsv-1,jev+1 ;
do i=isv-1,iev+1
1801 eta_pred(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * cs%IareaT(i,j)) * &
1802 ((uhbt(i-1,j) - uhbt(i,j)) + (vhbt(i,j-1) - vhbt(i,j)))
1806 do j=jsv-1,jev+1 ;
do i=isv-1,iev+1
1807 eta_pred(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * cs%IareaT(i,j)) * &
1808 (((datu(i-1,j)*ubt(i-1,j) + uhbt0(i-1,j)) - &
1809 (datu(i,j)*ubt(i,j) + uhbt0(i,j))) + &
1810 ((datv(i,j-1)*vbt(i,j-1) + vhbt0(i,j-1)) - &
1811 (datv(i,j)*vbt(i,j) + vhbt0(i,j))))
1815 if (cs%dynamic_psurf)
then 1817 do j=jsv-1,jev+1 ;
do i=isv-1,iev+1
1818 p_surf_dyn(i,j) = dyn_coef_eta(i,j) * (eta_pred(i,j) - eta(i,j))
1826 if (find_etaav)
then 1828 do j=js,je ;
do i=is,ie
1829 eta_sum(i,j) = eta_sum(i,j) + wt_accel2(n) * eta_pf_bt(i,j)
1834 if (interp_eta_pf)
then 1837 do j=jsv-1,jev+1 ;
do i=isv-1,iev+1
1838 eta_pf(i,j) = eta_pf_1(i,j) + wt_end*d_eta_pf(i,j)
1842 if (apply_obc_flather .or. apply_obc_open)
then 1844 do j=jsv,jev ;
do i=isv-2,iev+1
1845 ubt_old(i,j) = ubt(i,j)
1848 do j=jsv-2,jev+1 ;
do i=isv,iev
1849 vbt_old(i,j) = vbt(i,j)
1853 if (apply_obcs)
then 1854 if (mod(n+g%first_direction,2)==1)
then 1860 if (cs%BT_OBC%apply_u_OBCs)
then 1862 do j=jsv-joff,jev+joff ;
do i=isv-1,iev
1863 ubt_prev(i,j) = ubt(i,j) ; uhbt_prev(i,j) = uhbt(i,j)
1864 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)
1868 if (cs%BT_OBC%apply_v_OBCs)
then 1870 do j=jsv-1,jev ;
do i=isv-ioff,iev+ioff
1871 vbt_prev(i,j) = vbt(i,j) ; vhbt_prev(i,j) = vhbt(i,j)
1872 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)
1877 if (mod(n+g%first_direction,2)==1)
then 1880 do j=jsv-1,jev ;
do i=isv-1,iev+1
1881 cor_v(i,j) = -1.0*((amer(i-1,j) * ubt(i-1,j) + cmer(i,j+1) * ubt(i,j+1)) + &
1882 (bmer(i,j) * ubt(i,j) + dmer(i-1,j+1) * ubt(i-1,j+1))) - cor_ref_v(i,j)
1883 pfv(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_n(i,j) - &
1884 (eta_pf_bt(i,j+1)-eta_pf(i,j+1))*gtot_s(i,j+1)) * &
1885 dgeo_de * cs%IdyCv(i,j)
1888 if (cs%dynamic_psurf)
then 1890 do j=jsv-1,jev ;
do i=isv-1,iev+1
1891 pfv(i,j) = pfv(i,j) + (p_surf_dyn(i,j) - p_surf_dyn(i,j+1)) * cs%IdyCv(i,j)
1896 if (cs%BT_OBC%apply_v_OBCs)
then 1898 do j=jsv-1,jev ;
do i=isv-1,iev+1 ;
if (obc%segnum_v(i,j) /= obc_none)
then 1900 endif ;
enddo ;
enddo 1904 do j=jsv-1,jev ;
do i=isv-1,iev+1
1906 vbt(i,j) = bt_rem_v(i,j) * (vbt(i,j) + &
1907 dtbt * ((bt_force_v(i,j) + cor_v(i,j)) + pfv(i,j)))
1908 vbt_trans(i,j) = trans_wt1*vbt(i,j) + trans_wt2*vel_prev
1910 if (cs%linear_wave_drag)
then 1911 v_accel_bt(i,j) = v_accel_bt(i,j) + wt_accel(n) * &
1912 ((cor_v(i,j) + pfv(i,j)) - vbt(i,j)*rayleigh_v(i,j))
1914 v_accel_bt(i,j) = v_accel_bt(i,j) + wt_accel(n) * (cor_v(i,j) + pfv(i,j))
1918 if (integral_bt_cont)
then 1920 do j=jsv-1,jev ;
do i=isv-1,iev+1
1921 vbt_int(i,j) = vbt_int(i,j) + dtbt * vbt_trans(i,j)
1922 vhbt_int(i,j) = find_vhbt(vbt_int(i,j), btcl_v(i,j)) + n*dtbt*vhbt0(i,j)
1924 vhbt(i,j) = (vhbt_int(i,j) - vhbt_int_prev(i,j)) * idtbt
1926 elseif (use_bt_cont)
then 1928 do j=jsv-1,jev ;
do i=isv-1,iev+1
1929 vhbt(i,j) = find_vhbt(vbt_trans(i,j), btcl_v(i,j)) + vhbt0(i,j)
1934 do j=jsv-1,jev ;
do i=isv-1,iev+1
1935 vhbt(i,j) = datv(i,j)*vbt_trans(i,j) + vhbt0(i,j)
1939 if (cs%BT_OBC%apply_v_OBCs)
then 1941 do j=jsv-1,jev ;
do i=isv-1,iev+1 ;
if (obc%segnum_v(i,j) /= obc_none)
then 1942 vbt(i,j) = vbt_prev(i,j) ; vhbt(i,j) = vhbt_prev(i,j)
1943 endif ;
enddo ;
enddo 1947 do j=jsv,jev ;
do i=isv-1,iev
1948 cor_u(i,j) = ((azon(i,j) * vbt(i+1,j) + czon(i,j) * vbt(i,j-1)) + &
1949 (bzon(i,j) * vbt(i,j) + dzon(i,j) * vbt(i+1,j-1))) - &
1951 pfu(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_e(i,j) - &
1952 (eta_pf_bt(i+1,j)-eta_pf(i+1,j))*gtot_w(i+1,j)) * &
1953 dgeo_de * cs%IdxCu(i,j)
1957 if (cs%dynamic_psurf)
then 1959 do j=jsv,jev ;
do i=isv-1,iev
1960 pfu(i,j) = pfu(i,j) + (p_surf_dyn(i,j) - p_surf_dyn(i+1,j)) * cs%IdxCu(i,j)
1965 if (cs%BT_OBC%apply_u_OBCs)
then 1967 do j=jsv,jev ;
do i=isv-1,iev ;
if (obc%segnum_u(i,j) /= obc_none)
then 1969 endif ;
enddo ;
enddo 1973 do j=jsv,jev ;
do i=isv-1,iev
1975 ubt(i,j) = bt_rem_u(i,j) * (ubt(i,j) + &
1976 dtbt * ((bt_force_u(i,j) + cor_u(i,j)) + pfu(i,j)))
1977 if (abs(ubt(i,j)) < cs%vel_underflow) ubt(i,j) = 0.0
1978 ubt_trans(i,j) = trans_wt1*ubt(i,j) + trans_wt2*vel_prev
1980 if (cs%linear_wave_drag)
then 1981 u_accel_bt(i,j) = u_accel_bt(i,j) + wt_accel(n) * &
1982 ((cor_u(i,j) + pfu(i,j)) - ubt(i,j)*rayleigh_u(i,j))
1984 u_accel_bt(i,j) = u_accel_bt(i,j) + wt_accel(n) * (cor_u(i,j) + pfu(i,j))
1989 if (integral_bt_cont)
then 1991 do j=jsv,jev ;
do i=isv-1,iev
1992 ubt_int(i,j) = ubt_int(i,j) + dtbt * ubt_trans(i,j)
1993 uhbt_int(i,j) = find_uhbt(ubt_int(i,j), btcl_u(i,j)) + n*dtbt*uhbt0(i,j)
1995 uhbt(i,j) = (uhbt_int(i,j) - uhbt_int_prev(i,j)) * idtbt
1997 elseif (use_bt_cont)
then 1999 do j=jsv,jev ;
do i=isv-1,iev
2000 uhbt(i,j) = find_uhbt(ubt_trans(i,j), btcl_u(i,j)) + uhbt0(i,j)
2004 do j=jsv,jev ;
do i=isv-1,iev
2005 uhbt(i,j) = datu(i,j)*ubt_trans(i,j) + uhbt0(i,j)
2008 if (cs%BT_OBC%apply_u_OBCs)
then 2010 do j=jsv,jev ;
do i=isv-1,iev ;
if (obc%segnum_u(i,j) /= obc_none)
then 2011 ubt(i,j) = ubt_prev(i,j) ; uhbt(i,j) = uhbt_prev(i,j)
2012 endif ;
enddo ;
enddo 2017 do j=jsv-1,jev+1 ;
do i=isv-1,iev
2018 cor_u(i,j) = ((azon(i,j) * vbt(i+1,j) + czon(i,j) * vbt(i,j-1)) + &
2019 (bzon(i,j) * vbt(i,j) + dzon(i,j) * vbt(i+1,j-1))) - &
2021 pfu(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_e(i,j) - &
2022 (eta_pf_bt(i+1,j)-eta_pf(i+1,j))*gtot_w(i+1,j)) * &
2023 dgeo_de * cs%IdxCu(i,j)
2027 if (cs%dynamic_psurf)
then 2029 do j=jsv-1,jev+1 ;
do i=isv-1,iev
2030 pfu(i,j) = pfu(i,j) + (p_surf_dyn(i,j) - p_surf_dyn(i+1,j)) * cs%IdxCu(i,j)
2035 if (cs%BT_OBC%apply_u_OBCs)
then 2037 do j=jsv,jev ;
do i=isv-1,iev ;
if (obc%segnum_u(i,j) /= obc_none)
then 2039 endif ;
enddo ;
enddo 2043 do j=jsv-1,jev+1 ;
do i=isv-1,iev
2045 ubt(i,j) = bt_rem_u(i,j) * (ubt(i,j) + &
2046 dtbt * ((bt_force_u(i,j) + cor_u(i,j)) + pfu(i,j)))
2047 if (abs(ubt(i,j)) < cs%vel_underflow) ubt(i,j) = 0.0
2048 ubt_trans(i,j) = trans_wt1*ubt(i,j) + trans_wt2*vel_prev
2050 if (cs%linear_wave_drag)
then 2051 u_accel_bt(i,j) = u_accel_bt(i,j) + wt_accel(n) * &
2052 ((cor_u(i,j) + pfu(i,j)) - ubt(i,j)*rayleigh_u(i,j))
2054 u_accel_bt(i,j) = u_accel_bt(i,j) + wt_accel(n) * (cor_u(i,j) + pfu(i,j))
2058 if (integral_bt_cont)
then 2060 do j=jsv-1,jev+1 ;
do i=isv-1,iev
2061 ubt_int(i,j) = ubt_int(i,j) + dtbt * ubt_trans(i,j)
2062 uhbt_int(i,j) = find_uhbt(ubt_int(i,j), btcl_u(i,j)) + n*dtbt*uhbt0(i,j)
2064 uhbt(i,j) = (uhbt_int(i,j) - uhbt_int_prev(i,j)) * idtbt
2066 elseif (use_bt_cont)
then 2068 do j=jsv-1,jev+1 ;
do i=isv-1,iev
2069 uhbt(i,j) = find_uhbt(ubt_trans(i,j), btcl_u(i,j)) + uhbt0(i,j)
2074 do j=jsv-1,jev+1 ;
do i=isv-1,iev
2075 uhbt(i,j) = datu(i,j)*ubt_trans(i,j) + uhbt0(i,j)
2079 if (cs%BT_OBC%apply_u_OBCs)
then 2081 do j=jsv-1,jev+1 ;
do i=isv-1,iev ;
if (obc%segnum_u(i,j) /= obc_none)
then 2082 ubt(i,j) = ubt_prev(i,j) ; uhbt(i,j) = uhbt_prev(i,j)
2083 endif ;
enddo ;
enddo 2087 if (cs%use_old_coriolis_bracket_bug)
then 2089 do j=jsv-1,jev ;
do i=isv,iev
2090 cor_v(i,j) = -1.0*((amer(i-1,j) * ubt(i-1,j) + bmer(i,j) * ubt(i,j)) + &
2091 (cmer(i,j+1) * ubt(i,j+1) + dmer(i-1,j+1) * ubt(i-1,j+1))) - cor_ref_v(i,j)
2092 pfv(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_n(i,j) - &
2093 (eta_pf_bt(i,j+1)-eta_pf(i,j+1))*gtot_s(i,j+1)) * &
2094 dgeo_de * cs%IdyCv(i,j)
2099 do j=jsv-1,jev ;
do i=isv,iev
2100 cor_v(i,j) = -1.0*((amer(i-1,j) * ubt(i-1,j) + cmer(i,j+1) * ubt(i,j+1)) + &
2101 (bmer(i,j) * ubt(i,j) + dmer(i-1,j+1) * ubt(i-1,j+1))) - cor_ref_v(i,j)
2102 pfv(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_n(i,j) - &
2103 (eta_pf_bt(i,j+1)-eta_pf(i,j+1))*gtot_s(i,j+1)) * &
2104 dgeo_de * cs%IdyCv(i,j)
2109 if (cs%dynamic_psurf)
then 2111 do j=jsv-1,jev ;
do i=isv,iev
2112 pfv(i,j) = pfv(i,j) + (p_surf_dyn(i,j) - p_surf_dyn(i,j+1)) * cs%IdyCv(i,j)
2117 if (cs%BT_OBC%apply_v_OBCs)
then 2119 do j=jsv-1,jev ;
do i=isv-1,iev+1 ;
if (obc%segnum_v(i,j) /= obc_none)
then 2121 endif ;
enddo ;
enddo 2125 do j=jsv-1,jev ;
do i=isv,iev
2127 vbt(i,j) = bt_rem_v(i,j) * (vbt(i,j) + &
2128 dtbt * ((bt_force_v(i,j) + cor_v(i,j)) + pfv(i,j)))
2129 if (abs(vbt(i,j)) < cs%vel_underflow) vbt(i,j) = 0.0
2130 vbt_trans(i,j) = trans_wt1*vbt(i,j) + trans_wt2*vel_prev
2132 if (cs%linear_wave_drag)
then 2133 v_accel_bt(i,j) = v_accel_bt(i,j) + wt_accel(n) * &
2134 ((cor_v(i,j) + pfv(i,j)) - vbt(i,j)*rayleigh_v(i,j))
2136 v_accel_bt(i,j) = v_accel_bt(i,j) + wt_accel(n) * (cor_v(i,j) + pfv(i,j))
2140 if (integral_bt_cont)
then 2142 do j=jsv-1,jev ;
do i=isv,iev
2143 vbt_int(i,j) = vbt_int(i,j) + dtbt * vbt_trans(i,j)
2144 vhbt_int(i,j) = find_vhbt(vbt_int(i,j), btcl_v(i,j)) + n*dtbt*vhbt0(i,j)
2146 vhbt(i,j) = (vhbt_int(i,j) - vhbt_int_prev(i,j)) * idtbt
2148 elseif (use_bt_cont)
then 2150 do j=jsv-1,jev ;
do i=isv,iev
2151 vhbt(i,j) = find_vhbt(vbt_trans(i,j), btcl_v(i,j)) + vhbt0(i,j)
2155 do j=jsv-1,jev ;
do i=isv,iev
2156 vhbt(i,j) = datv(i,j)*vbt_trans(i,j) + vhbt0(i,j)
2159 if (cs%BT_OBC%apply_v_OBCs)
then 2161 do j=jsv-1,jev ;
do i=isv,iev ;
if (obc%segnum_v(i,j) /= obc_none)
then 2162 vbt(i,j) = vbt_prev(i,j); vhbt(i,j) = vhbt_prev(i,j)
2163 endif ;
enddo ;
enddo 2168 if (cs%debug_bt)
then 2169 write(mesg,
'("BT vel update ",I4)') n
2170 call uvchksum(trim(mesg)//
" PF[uv]", pfu, pfv, cs%debug_BT_HI, haloshift=iev-ie, &
2171 scale=us%L_T_to_m_s*us%s_to_T)
2172 call uvchksum(trim(mesg)//
" Cor_[uv]", cor_u, cor_v, cs%debug_BT_HI, haloshift=iev-ie, &
2173 scale=us%L_T_to_m_s*us%s_to_T)
2174 call uvchksum(trim(mesg)//
" BT_force_[uv]", bt_force_u, bt_force_v, cs%debug_BT_HI, haloshift=iev-ie, &
2175 scale=us%L_T_to_m_s*us%s_to_T)
2176 call uvchksum(trim(mesg)//
" BT_rem_[uv]", bt_rem_u, bt_rem_v, cs%debug_BT_HI, haloshift=iev-ie)
2177 call uvchksum(trim(mesg)//
" [uv]bt", ubt, vbt, cs%debug_BT_HI, haloshift=iev-ie, &
2178 scale=us%L_T_to_m_s)
2179 call uvchksum(trim(mesg)//
" [uv]bt_trans", ubt_trans, vbt_trans, cs%debug_BT_HI, haloshift=iev-ie, &
2180 scale=us%L_T_to_m_s)
2181 call uvchksum(trim(mesg)//
" [uv]hbt", uhbt, vhbt, cs%debug_BT_HI, haloshift=iev-ie, &
2182 scale=us%s_to_T*us%L_to_m**2*gv%H_to_m)
2183 if (integral_bt_cont) &
2184 call uvchksum(trim(mesg)//
" [uv]hbt_int", uhbt_int, vhbt_int, cs%debug_BT_HI, haloshift=iev-ie, &
2185 scale=us%L_to_m**2*gv%H_to_m)
2190 do j=js,je ;
do i=is-1,ie
2191 pfu_bt_sum(i,j) = pfu_bt_sum(i,j) + wt_accel2(n) * pfu(i,j)
2195 do j=js-1,je ;
do i=is,ie
2196 pfv_bt_sum(i,j) = pfv_bt_sum(i,j) + wt_accel2(n) * pfv(i,j)
2202 do j=js,je ;
do i=is-1,ie
2203 coru_bt_sum(i,j) = coru_bt_sum(i,j) + wt_accel2(n) * cor_u(i,j)
2207 do j=js-1,je ;
do i=is,ie
2208 corv_bt_sum(i,j) = corv_bt_sum(i,j) + wt_accel2(n) * cor_v(i,j)
2214 do j=js,je ;
do i=is-1,ie
2215 ubt_sum(i,j) = ubt_sum(i,j) + wt_trans(n) * ubt_trans(i,j)
2216 uhbt_sum(i,j) = uhbt_sum(i,j) + wt_trans(n) * uhbt(i,j)
2217 ubt_wtd(i,j) = ubt_wtd(i,j) + wt_vel(n) * ubt(i,j)
2221 do j=js-1,je ;
do i=is,ie
2222 vbt_sum(i,j) = vbt_sum(i,j) + wt_trans(n) * vbt_trans(i,j)
2223 vhbt_sum(i,j) = vhbt_sum(i,j) + wt_trans(n) * vhbt(i,j)
2224 vbt_wtd(i,j) = vbt_wtd(i,j) + wt_vel(n) * vbt(i,j)
2228 if (apply_obcs)
then 2231 call apply_velocity_obcs(obc, ubt, vbt, uhbt, vhbt, &
2232 ubt_trans, vbt_trans, eta, ubt_old, vbt_old, cs%BT_OBC, &
2233 g, ms, us, iev-ie, dtbt, bebt, use_bt_cont, integral_bt_cont, &
2234 n*dtbt, datu, datv, btcl_u, btcl_v, uhbt0, vhbt0, &
2235 ubt_int_prev, vbt_int_prev, uhbt_int_prev, vhbt_int_prev)
2238 if (cs%BT_OBC%apply_u_OBCs)
then 2240 do j=js,je ;
do i=is-1,ie
2241 if (obc%segnum_u(i,j) /= obc_none)
then 2243 ubt_sum(i,j) = ubt_sum_prev(i,j) + wt_trans(n) * ubt_trans(i,j)
2244 uhbt_sum(i,j) = uhbt_sum_prev(i,j) + wt_trans(n) * uhbt(i,j)
2245 ubt_wtd(i,j) = ubt_wtd_prev(i,j) + wt_vel(n) * ubt(i,j)
2246 if (integral_bt_cont)
then 2247 uhbt_int(i,j) = uhbt_int_prev(i,j) + dtbt * uhbt(i,j)
2248 ubt_int(i,j) = ubt_int_prev(i,j) + dtbt * ubt_trans(i,j)
2253 if (cs%BT_OBC%apply_v_OBCs)
then 2255 do j=js-1,je ;
do i=is,ie
2256 if (obc%segnum_v(i,j) /= obc_none)
then 2258 vbt_sum(i,j) = vbt_sum_prev(i,j) + wt_trans(n) * vbt_trans(i,j)
2259 vhbt_sum(i,j) = vhbt_sum_prev(i,j) + wt_trans(n) * vhbt(i,j)
2260 vbt_wtd(i,j) = vbt_wtd_prev(i,j) + wt_vel(n) * vbt(i,j)
2261 if (integral_bt_cont)
then 2262 vbt_int(i,j) = vbt_int_prev(i,j) + dtbt * vbt_trans(i,j)
2263 vhbt_int(i,j) = vhbt_int_prev(i,j) + dtbt * vhbt(i,j)
2270 if (cs%debug_bt)
then 2271 call uvchksum(
"BT [uv]hbt just after OBC", uhbt, vhbt, cs%debug_BT_HI, haloshift=iev-ie, &
2272 scale=us%s_to_T*us%L_to_m**2*gv%H_to_m)
2273 if (integral_bt_cont) &
2274 call uvchksum(
"BT [uv]hbt_int just after OBC", uhbt_int, vhbt_int, cs%debug_BT_HI, &
2275 haloshift=iev-ie, scale=us%L_to_m**2*gv%H_to_m)
2278 if (integral_bt_cont)
then 2280 do j=jsv,jev ;
do i=isv,iev
2281 eta(i,j) = (eta_ic(i,j) + n*eta_src(i,j)) + cs%IareaT(i,j) * &
2282 ((uhbt_int(i-1,j) - uhbt_int(i,j)) + (vhbt_int(i,j-1) - vhbt_int(i,j)))
2283 eta_wtd(i,j) = eta_wtd(i,j) + eta(i,j) * wt_eta(n)
2287 do j=jsv,jev ;
do i=isv,iev
2288 eta(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * cs%IareaT(i,j)) * &
2289 ((uhbt(i-1,j) - uhbt(i,j)) + (vhbt(i,j-1) - vhbt(i,j)))
2290 eta_wtd(i,j) = eta_wtd(i,j) + eta(i,j) * wt_eta(n)
2295 if (do_hifreq_output)
then 2296 time_step_end = time_bt_start + real_to_time(n*us%T_to_s*dtbt)
2297 call enable_averaging(us%T_to_s*dtbt, time_step_end, cs%diag)
2298 if (cs%id_ubt_hifreq > 0)
call post_data(cs%id_ubt_hifreq, ubt(isdb:iedb,jsd:jed), cs%diag)
2299 if (cs%id_vbt_hifreq > 0)
call post_data(cs%id_vbt_hifreq, vbt(isd:ied,jsdb:jedb), cs%diag)
2300 if (cs%id_eta_hifreq > 0)
call post_data(cs%id_eta_hifreq, eta(isd:ied,jsd:jed), cs%diag)
2301 if (cs%id_uhbt_hifreq > 0)
call post_data(cs%id_uhbt_hifreq, uhbt(isdb:iedb,jsd:jed), cs%diag)
2302 if (cs%id_vhbt_hifreq > 0)
call post_data(cs%id_vhbt_hifreq, vhbt(isd:ied,jsdb:jedb), cs%diag)
2303 if (cs%id_eta_pred_hifreq > 0)
call post_data(cs%id_eta_pred_hifreq, eta_pf_bt(isd:ied,jsd:jed), cs%diag)
2306 if (cs%debug_bt)
then 2307 write(mesg,
'("BT step ",I4)') n
2308 call uvchksum(trim(mesg)//
" [uv]bt", ubt, vbt, cs%debug_BT_HI, haloshift=iev-ie, &
2309 scale=us%L_T_to_m_s)
2310 call hchksum(eta, trim(mesg)//
" eta", cs%debug_BT_HI, haloshift=iev-ie, scale=gv%H_to_m)
2313 if (gv%Boussinesq)
then 2314 do j=js,je ;
do i=is,ie
2315 if (eta(i,j) < -gv%Z_to_H*g%bathyT(i,j)) &
2316 call mom_error(warning,
"btstep: eta has dropped below bathyT.")
2319 do j=js,je ;
do i=is,ie
2320 if (eta(i,j) < 0.0) &
2321 call mom_error(warning,
"btstep: negative eta in a non-Boussinesq barotropic solver.")
2326 if (id_clock_calc > 0)
call cpu_clock_end(id_clock_calc)
2327 if (id_clock_calc_post > 0)
call cpu_clock_begin(id_clock_calc_post)
2330 if (do_hifreq_output)
call enable_averaging(time_int_in, time_end_in, cs%diag)
2332 if (cs%answers_2018)
then 2333 i_sum_wt_vel = 1.0 / sum_wt_vel ; i_sum_wt_eta = 1.0 / sum_wt_eta
2334 i_sum_wt_accel = 1.0 / sum_wt_accel ; i_sum_wt_trans = 1.0 / sum_wt_trans
2336 i_sum_wt_vel = 1.0 ; i_sum_wt_eta = 1.0 ; i_sum_wt_accel = 1.0 ; i_sum_wt_trans = 1.0
2339 if (find_etaav)
then ;
do j=js,je ;
do i=is,ie
2340 etaav(i,j) = eta_sum(i,j) * i_sum_wt_accel
2341 enddo ;
enddo ;
endif 2342 do j=js-1,je+1 ;
do i=is-1,ie+1 ; e_anom(i,j) = 0.0 ;
enddo ;
enddo 2343 if (interp_eta_pf)
then 2344 do j=js,je ;
do i=is,ie
2345 e_anom(i,j) = dgeo_de * (0.5 * (eta(i,j) + eta_in(i,j)) - &
2346 (eta_pf_1(i,j) + 0.5*d_eta_pf(i,j)))
2349 do j=js,je ;
do i=is,ie
2350 e_anom(i,j) = dgeo_de * (0.5 * (eta(i,j) + eta_in(i,j)) - eta_pf(i,j))
2353 if (apply_obcs)
then 2355 if (cs%BT_OBC%apply_u_OBCs)
then 2357 do j=js,je ;
do i=is-1,ie
2358 l_seg = obc%segnum_u(i,j)
2359 if (l_seg == obc_none) cycle
2361 if (obc%segment(l_seg)%direction == obc_direction_e)
then 2362 e_anom(i+1,j) = e_anom(i,j)
2363 elseif (obc%segment(l_seg)%direction == obc_direction_w)
then 2364 e_anom(i,j) = e_anom(i+1,j)
2369 if (cs%BT_OBC%apply_v_OBCs)
then 2371 do j=js-1,je ;
do i=is,ie
2372 l_seg = obc%segnum_v(i,j)
2373 if (l_seg == obc_none) cycle
2375 if (obc%segment(l_seg)%direction == obc_direction_n)
then 2376 e_anom(i,j+1) = e_anom(i,j)
2377 elseif (obc%segment(l_seg)%direction == obc_direction_s)
then 2378 e_anom(i,j) = e_anom(i,j+1)
2385 do j=js,je ;
do i=is,ie
2386 eta_out(i,j) = eta_wtd(i,j) * i_sum_wt_eta
2389 if (id_clock_calc_post > 0)
call cpu_clock_end(id_clock_calc_post)
2390 if (id_clock_pass_post > 0)
call cpu_clock_begin(id_clock_pass_post)
2391 if (g%nonblocking_updates)
then 2392 call start_group_pass(cs%pass_e_anom, g%Domain)
2394 if (find_etaav)
call do_group_pass(cs%pass_etaav, g%Domain)
2395 call do_group_pass(cs%pass_e_anom, g%Domain)
2397 if (id_clock_pass_post > 0)
call cpu_clock_end(id_clock_pass_post)
2398 if (id_clock_calc_post > 0)
call cpu_clock_begin(id_clock_calc_post)
2400 if (cs%answers_2018)
then 2401 do j=js,je ;
do i=is-1,ie
2402 cs%ubtav(i,j) = ubt_sum(i,j) * i_sum_wt_trans
2403 uhbtav(i,j) = uhbt_sum(i,j) * i_sum_wt_trans
2404 ubt_wtd(i,j) = ubt_wtd(i,j) * i_sum_wt_vel
2407 do j=js-1,je ;
do i=is,ie
2408 cs%vbtav(i,j) = vbt_sum(i,j) * i_sum_wt_trans
2409 vhbtav(i,j) = vhbt_sum(i,j) * i_sum_wt_trans
2410 vbt_wtd(i,j) = vbt_wtd(i,j) * i_sum_wt_vel
2413 do j=js,je ;
do i=is-1,ie
2414 cs%ubtav(i,j) = ubt_sum(i,j)
2415 uhbtav(i,j) = uhbt_sum(i,j)
2418 do j=js-1,je ;
do i=is,ie
2419 cs%vbtav(i,j) = vbt_sum(i,j)
2420 vhbtav(i,j) = vhbt_sum(i,j)
2425 if (id_clock_calc_post > 0)
call cpu_clock_end(id_clock_calc_post)
2426 if (id_clock_pass_post > 0)
call cpu_clock_begin(id_clock_pass_post)
2427 if (g%nonblocking_updates)
then 2428 call complete_group_pass(cs%pass_e_anom, g%Domain)
2429 if (find_etaav)
call start_group_pass(cs%pass_etaav, g%Domain)
2430 call start_group_pass(cs%pass_ubta_uhbta, g%DoMain)
2432 call do_group_pass(cs%pass_ubta_uhbta, g%Domain)
2434 if (id_clock_pass_post > 0)
call cpu_clock_end(id_clock_pass_post)
2435 if (id_clock_calc_post > 0)
call cpu_clock_begin(id_clock_calc_post)
2440 do j=js,je ;
do i=is-1,ie
2441 accel_layer_u(i,j,k) = (u_accel_bt(i,j) - &
2442 ((pbce(i+1,j,k) - gtot_w(i+1,j)) * e_anom(i+1,j) - &
2443 (pbce(i,j,k) - gtot_e(i,j)) * e_anom(i,j)) * cs%IdxCu(i,j) )
2444 if (abs(accel_layer_u(i,j,k)) < accel_underflow) accel_layer_u(i,j,k) = 0.0
2446 do j=js-1,je ;
do i=is,ie
2447 accel_layer_v(i,j,k) = (v_accel_bt(i,j) - &
2448 ((pbce(i,j+1,k) - gtot_s(i,j+1)) * e_anom(i,j+1) - &
2449 (pbce(i,j,k) - gtot_n(i,j)) * e_anom(i,j)) * cs%IdyCv(i,j) )
2450 if (abs(accel_layer_v(i,j,k)) < accel_underflow) accel_layer_v(i,j,k) = 0.0
2454 if (apply_obcs)
then 2457 if (cs%BT_OBC%apply_u_OBCs)
then ;
do j=js,je ;
do i=is-1,ie
2458 if (obc%segnum_u(i,j) /= obc_none)
then 2459 u_accel_bt(i,j) = (ubt_wtd(i,j) - ubt_first(i,j)) / dt
2460 do k=1,nz ; accel_layer_u(i,j,k) = u_accel_bt(i,j) ;
enddo 2462 enddo ;
enddo ;
endif 2463 if (cs%BT_OBC%apply_v_OBCs)
then ;
do j=js-1,je ;
do i=is,ie
2464 if (obc%segnum_v(i,j) /= obc_none)
then 2465 v_accel_bt(i,j) = (vbt_wtd(i,j) - vbt_first(i,j)) / dt
2466 do k=1,nz ; accel_layer_v(i,j,k) = v_accel_bt(i,j) ;
enddo 2468 enddo ;
enddo ;
endif 2471 if (id_clock_calc_post > 0)
call cpu_clock_end(id_clock_calc_post)
2474 if (query_averaging_enabled(cs%diag))
then 2476 if (cs%gradual_BT_ICs)
then 2477 do j=js,je ;
do i=is-1,ie ; cs%ubt_IC(i,j) = ubt_wtd(i,j) ;
enddo ;
enddo 2478 do j=js-1,je ;
do i=is,ie ; cs%vbt_IC(i,j) = vbt_wtd(i,j) ;
enddo ;
enddo 2482 if (cs%id_PFu_bt > 0)
then 2483 do j=js,je ;
do i=is-1,ie
2484 pfu_bt_sum(i,j) = pfu_bt_sum(i,j) * i_sum_wt_accel
2486 call post_data(cs%id_PFu_bt, pfu_bt_sum(isdb:iedb,jsd:jed), cs%diag)
2488 if (cs%id_PFv_bt > 0)
then 2489 do j=js-1,je ;
do i=is,ie
2490 pfv_bt_sum(i,j) = pfv_bt_sum(i,j) * i_sum_wt_accel
2492 call post_data(cs%id_PFv_bt, pfv_bt_sum(isd:ied,jsdb:jedb), cs%diag)
2494 if (cs%id_Coru_bt > 0)
then 2495 do j=js,je ;
do i=is-1,ie
2496 coru_bt_sum(i,j) = coru_bt_sum(i,j) * i_sum_wt_accel
2498 call post_data(cs%id_Coru_bt, coru_bt_sum(isdb:iedb,jsd:jed), cs%diag)
2500 if (cs%id_Corv_bt > 0)
then 2501 do j=js-1,je ;
do i=is,ie
2502 corv_bt_sum(i,j) = corv_bt_sum(i,j) * i_sum_wt_accel
2504 call post_data(cs%id_Corv_bt, corv_bt_sum(isd:ied,jsdb:jedb), cs%diag)
2506 if (cs%id_ubtdt > 0)
then 2507 do j=js,je ;
do i=is-1,ie
2508 ubt_dt(i,j) = (ubt_wtd(i,j) - ubt_st(i,j))*idt
2510 call post_data(cs%id_ubtdt, ubt_dt(isdb:iedb,jsd:jed), cs%diag)
2512 if (cs%id_vbtdt > 0)
then 2513 do j=js-1,je ;
do i=is,ie
2514 vbt_dt(i,j) = (vbt_wtd(i,j) - vbt_st(i,j))*idt
2516 call post_data(cs%id_vbtdt, vbt_dt(isd:ied,jsdb:jedb), cs%diag)
2519 if (cs%id_ubtforce > 0)
call post_data(cs%id_ubtforce, bt_force_u(isdb:iedb,jsd:jed), cs%diag)
2520 if (cs%id_vbtforce > 0)
call post_data(cs%id_vbtforce, bt_force_v(isd:ied,jsdb:jedb), cs%diag)
2521 if (cs%id_uaccel > 0)
call post_data(cs%id_uaccel, u_accel_bt(isdb:iedb,jsd:jed), cs%diag)
2522 if (cs%id_vaccel > 0)
call post_data(cs%id_vaccel, v_accel_bt(isd:ied,jsdb:jedb), cs%diag)
2524 if (cs%id_eta_cor > 0)
call post_data(cs%id_eta_cor, cs%eta_cor, cs%diag)
2525 if (cs%id_eta_bt > 0)
call post_data(cs%id_eta_bt, eta_out, cs%diag)
2526 if (cs%id_gtotn > 0)
call post_data(cs%id_gtotn, gtot_n(isd:ied,jsd:jed), cs%diag)
2527 if (cs%id_gtots > 0)
call post_data(cs%id_gtots, gtot_s(isd:ied,jsd:jed), cs%diag)
2528 if (cs%id_gtote > 0)
call post_data(cs%id_gtote, gtot_e(isd:ied,jsd:jed), cs%diag)
2529 if (cs%id_gtotw > 0)
call post_data(cs%id_gtotw, gtot_w(isd:ied,jsd:jed), cs%diag)
2530 if (cs%id_ubt > 0)
call post_data(cs%id_ubt, ubt_wtd(isdb:iedb,jsd:jed), cs%diag)
2531 if (cs%id_vbt > 0)
call post_data(cs%id_vbt, vbt_wtd(isd:ied,jsdb:jedb), cs%diag)
2532 if (cs%id_ubtav > 0)
call post_data(cs%id_ubtav, cs%ubtav, cs%diag)
2533 if (cs%id_vbtav > 0)
call post_data(cs%id_vbtav, cs%vbtav, cs%diag)
2534 if (cs%id_visc_rem_u > 0)
call post_data(cs%id_visc_rem_u, visc_rem_u, cs%diag)
2535 if (cs%id_visc_rem_v > 0)
call post_data(cs%id_visc_rem_v, visc_rem_v, cs%diag)
2537 if (cs%id_frhatu > 0)
call post_data(cs%id_frhatu, cs%frhatu, cs%diag)
2538 if (cs%id_uhbt > 0)
call post_data(cs%id_uhbt, uhbtav, cs%diag)
2539 if (cs%id_frhatv > 0)
call post_data(cs%id_frhatv, cs%frhatv, cs%diag)
2540 if (cs%id_vhbt > 0)
call post_data(cs%id_vhbt, vhbtav, cs%diag)
2541 if (cs%id_uhbt0 > 0)
call post_data(cs%id_uhbt0, uhbt0(isdb:iedb,jsd:jed), cs%diag)
2542 if (cs%id_vhbt0 > 0)
call post_data(cs%id_vhbt0, vhbt0(isd:ied,jsdb:jedb), cs%diag)
2544 if (cs%id_frhatu1 > 0)
call post_data(cs%id_frhatu1, cs%frhatu1, cs%diag)
2545 if (cs%id_frhatv1 > 0)
call post_data(cs%id_frhatv1, cs%frhatv1, cs%diag)
2547 if (use_bt_cont)
then 2548 if (cs%id_BTC_FA_u_EE > 0)
call post_data(cs%id_BTC_FA_u_EE, bt_cont%FA_u_EE, cs%diag)
2549 if (cs%id_BTC_FA_u_E0 > 0)
call post_data(cs%id_BTC_FA_u_E0, bt_cont%FA_u_E0, cs%diag)
2550 if (cs%id_BTC_FA_u_W0 > 0)
call post_data(cs%id_BTC_FA_u_W0, bt_cont%FA_u_W0, cs%diag)
2551 if (cs%id_BTC_FA_u_WW > 0)
call post_data(cs%id_BTC_FA_u_WW, bt_cont%FA_u_WW, cs%diag)
2552 if (cs%id_BTC_uBT_EE > 0)
call post_data(cs%id_BTC_uBT_EE, bt_cont%uBT_EE, cs%diag)
2553 if (cs%id_BTC_uBT_WW > 0)
call post_data(cs%id_BTC_uBT_WW, bt_cont%uBT_WW, cs%diag)
2554 if (cs%id_BTC_FA_u_rat0 > 0)
then 2556 do j=js,je ;
do i=is-1,ie
2557 if ((g%mask2dCu(i,j) > 0.0) .and. (bt_cont%FA_u_W0(i,j) > 0.0))
then 2558 tmp_u(i,j) = (bt_cont%FA_u_E0(i,j)/ bt_cont%FA_u_W0(i,j))
2563 call post_data(cs%id_BTC_FA_u_rat0, tmp_u, cs%diag)
2565 if (cs%id_BTC_FA_v_NN > 0)
call post_data(cs%id_BTC_FA_v_NN, bt_cont%FA_v_NN, cs%diag)
2566 if (cs%id_BTC_FA_v_N0 > 0)
call post_data(cs%id_BTC_FA_v_N0, bt_cont%FA_v_N0, cs%diag)
2567 if (cs%id_BTC_FA_v_S0 > 0)
call post_data(cs%id_BTC_FA_v_S0, bt_cont%FA_v_S0, cs%diag)
2568 if (cs%id_BTC_FA_v_SS > 0)
call post_data(cs%id_BTC_FA_v_SS, bt_cont%FA_v_SS, cs%diag)
2569 if (cs%id_BTC_vBT_NN > 0)
call post_data(cs%id_BTC_vBT_NN, bt_cont%vBT_NN, cs%diag)
2570 if (cs%id_BTC_vBT_SS > 0)
call post_data(cs%id_BTC_vBT_SS, bt_cont%vBT_SS, cs%diag)
2571 if (cs%id_BTC_FA_v_rat0 > 0)
then 2573 do j=js-1,je ;
do i=is,ie
2574 if ((g%mask2dCv(i,j) > 0.0) .and. (bt_cont%FA_v_S0(i,j) > 0.0))
then 2575 tmp_v(i,j) = (bt_cont%FA_v_N0(i,j)/ bt_cont%FA_v_S0(i,j))
2580 call post_data(cs%id_BTC_FA_v_rat0, tmp_v, cs%diag)
2582 if (cs%id_BTC_FA_h_rat0 > 0)
then 2584 do j=js,je ;
do i=is,ie
2586 if ((g%mask2dCu(i,j) > 0.0) .and. (bt_cont%FA_u_W0(i,j) > 0.0) .and. (bt_cont%FA_u_E0(i,j) > 0.0))
then 2587 if (bt_cont%FA_u_W0(i,j) > bt_cont%FA_u_E0(i,j))
then 2588 tmp_h(i,j) = max(tmp_h(i,j), (bt_cont%FA_u_W0(i,j)/ bt_cont%FA_u_E0(i,j)))
2590 tmp_h(i,j) = max(tmp_h(i,j), (bt_cont%FA_u_E0(i,j)/ bt_cont%FA_u_W0(i,j)))
2593 if ((g%mask2dCu(i-1,j) > 0.0) .and. (bt_cont%FA_u_W0(i-1,j) > 0.0) .and. (bt_cont%FA_u_E0(i-1,j) > 0.0))
then 2594 if (bt_cont%FA_u_W0(i-1,j) > bt_cont%FA_u_E0(i-1,j))
then 2595 tmp_h(i,j) = max(tmp_h(i,j), (bt_cont%FA_u_W0(i-1,j)/ bt_cont%FA_u_E0(i-1,j)))
2597 tmp_h(i,j) = max(tmp_h(i,j), (bt_cont%FA_u_E0(i-1,j)/ bt_cont%FA_u_W0(i-1,j)))
2600 if ((g%mask2dCv(i,j) > 0.0) .and. (bt_cont%FA_v_S0(i,j) > 0.0) .and. (bt_cont%FA_v_N0(i,j) > 0.0))
then 2601 if (bt_cont%FA_v_S0(i,j) > bt_cont%FA_v_N0(i,j))
then 2602 tmp_h(i,j) = max(tmp_h(i,j), (bt_cont%FA_v_S0(i,j)/ bt_cont%FA_v_N0(i,j)))
2604 tmp_h(i,j) = max(tmp_h(i,j), (bt_cont%FA_v_N0(i,j)/ bt_cont%FA_v_S0(i,j)))
2607 if ((g%mask2dCv(i,j-1) > 0.0) .and. (bt_cont%FA_v_S0(i,j-1) > 0.0) .and. (bt_cont%FA_v_N0(i,j-1) > 0.0))
then 2608 if (bt_cont%FA_v_S0(i,j-1) > bt_cont%FA_v_N0(i,j-1))
then 2609 tmp_h(i,j) = max(tmp_h(i,j), (bt_cont%FA_v_S0(i,j-1)/ bt_cont%FA_v_N0(i,j-1)))
2611 tmp_h(i,j) = max(tmp_h(i,j), (bt_cont%FA_v_N0(i,j-1)/ bt_cont%FA_v_S0(i,j-1)))
2615 call post_data(cs%id_BTC_FA_h_rat0, tmp_h, cs%diag)
2619 if (cs%id_frhatu1 > 0) cs%frhatu1(:,:,:) = cs%frhatu(:,:,:)
2620 if (cs%id_frhatv1 > 0) cs%frhatv1(:,:,:) = cs%frhatv(:,:,:)
2623 if ((
present(adp)) .and. (
associated(adp%diag_hfrac_u)))
then 2624 do k=1,nz ;
do j=js,je ;
do i=is-1,ie
2625 adp%diag_hfrac_u(i,j,k) = cs%frhatu(i,j,k)
2626 enddo ;
enddo ;
enddo 2628 if ((
present(adp)) .and. (
associated(adp%diag_hfrac_v)))
then 2629 do k=1,nz ;
do j=js-1,je ;
do i=is,ie
2630 adp%diag_hfrac_v(i,j,k) = cs%frhatv(i,j,k)
2631 enddo ;
enddo ;
enddo 2634 if (g%nonblocking_updates)
then 2635 if (find_etaav)
call complete_group_pass(cs%pass_etaav, g%Domain)
2636 call complete_group_pass(cs%pass_ubta_uhbta, g%Domain)
2639 end subroutine btstep
2643 subroutine set_dtbt(G, GV, US, CS, eta, pbce, BT_cont, gtot_est, SSH_add)
2648 real,
dimension(SZI_(G),SZJ_(G)),
optional,
intent(in) :: eta
2650 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
optional,
intent(in) :: pbce
2656 real,
optional,
intent(in) :: gtot_est
2658 real,
optional,
intent(in) :: ssh_add
2663 real,
dimension(SZI_(G),SZJ_(G)) :: &
2670 real,
dimension(SZIBS_(G),SZJ_(G)) :: &
2673 real,
dimension(SZI_(G),SZJBS_(G)) :: &
2690 logical :: use_bt_cont
2693 character(len=200) :: mesg
2694 integer :: i, j, k, is, ie, js, je, nz
2696 if (.not.
associated(cs))
call mom_error(fatal, &
2697 "set_dtbt: Module MOM_barotropic must be initialized before it is used.")
2698 if (.not.cs%split)
return 2699 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
2700 ms%isdw = g%isd ; ms%iedw = g%ied ; ms%jsdw = g%jsd ; ms%jedw = g%jed
2702 if (.not.(
present(pbce) .or.
present(gtot_est)))
call mom_error(fatal, &
2703 "set_dtbt: Either pbce or gtot_est must be present.")
2705 add_ssh = 0.0 ;
if (
present(ssh_add)) add_ssh = ssh_add
2707 use_bt_cont = .false.
2708 if (
present(bt_cont)) use_bt_cont = (
associated(bt_cont))
2710 if (use_bt_cont)
then 2711 call bt_cont_to_face_areas(bt_cont, datu, datv, g, us, ms, 0, .true.)
2712 elseif (cs%Nonlinear_continuity .and.
present(eta))
then 2713 call find_face_areas(datu, datv, g, gv, us, cs, ms, eta=eta, halo=0)
2715 call find_face_areas(datu, datv, g, gv, us, cs, ms, halo=0, add_max=add_ssh)
2719 if (cs%tides)
call tidal_forcing_sensitivity(g, cs%tides_CSp, det_de)
2720 dgeo_de = 1.0 + max(0.0, det_de + cs%G_extra)
2721 if (
present(pbce))
then 2722 do j=js,je ;
do i=is,ie
2723 gtot_e(i,j) = 0.0 ; gtot_w(i,j) = 0.0
2724 gtot_n(i,j) = 0.0 ; gtot_s(i,j) = 0.0
2726 do k=1,nz ;
do j=js,je ;
do i=is,ie
2727 gtot_e(i,j) = gtot_e(i,j) + pbce(i,j,k) * cs%frhatu(i,j,k)
2728 gtot_w(i,j) = gtot_w(i,j) + pbce(i,j,k) * cs%frhatu(i-1,j,k)
2729 gtot_n(i,j) = gtot_n(i,j) + pbce(i,j,k) * cs%frhatv(i,j,k)
2730 gtot_s(i,j) = gtot_s(i,j) + pbce(i,j,k) * cs%frhatv(i,j-1,k)
2731 enddo ;
enddo ;
enddo 2733 do j=js,je ;
do i=is,ie
2734 gtot_e(i,j) = gtot_est * gv%H_to_Z ; gtot_w(i,j) = gtot_est * gv%H_to_Z
2735 gtot_n(i,j) = gtot_est * gv%H_to_Z ; gtot_s(i,j) = gtot_est * gv%H_to_Z
2739 min_max_dt2 = 1.0e38*us%s_to_T**2
2740 do j=js,je ;
do i=is,ie
2743 idt_max2 = 0.5 * (1.0 + 2.0*cs%bebt) * (g%IareaT(i,j) * &
2744 ((gtot_e(i,j)*datu(i,j)*g%IdxCu(i,j) + gtot_w(i,j)*datu(i-1,j)*g%IdxCu(i-1,j)) + &
2745 (gtot_n(i,j)*datv(i,j)*g%IdyCv(i,j) + gtot_s(i,j)*datv(i,j-1)*g%IdyCv(i,j-1))) + &
2746 ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
2747 (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2)) * cs%BT_Coriolis_scale**2 )
2748 if (idt_max2 * min_max_dt2 > 1.0) min_max_dt2 = 1.0 / idt_max2
2750 dtbt_max = sqrt(min_max_dt2 / dgeo_de)
2751 if (id_clock_sync > 0)
call cpu_clock_begin(id_clock_sync)
2752 call min_across_pes(dtbt_max)
2753 if (id_clock_sync > 0)
call cpu_clock_end(id_clock_sync)
2755 cs%dtbt = cs%dtbt_fraction * dtbt_max
2756 cs%dtbt_max = dtbt_max
2757 end subroutine set_dtbt
2762 subroutine apply_velocity_obcs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans, eta, &
2763 ubt_old, vbt_old, BT_OBC, G, MS, US, halo, dtbt, bebt, &
2764 use_BT_cont, integral_BT_cont, dt_elapsed, Datu, Datv, &
2765 BTCL_u, BTCL_v, uhbt0, vhbt0, ubt_int, vbt_int, uhbt_int, vhbt_int)
2770 real,
dimension(SZIBW_(MS),SZJW_(MS)),
intent(inout) :: ubt
2771 real,
dimension(SZIBW_(MS),SZJW_(MS)),
intent(inout) :: uhbt
2773 real,
dimension(SZIBW_(MS),SZJW_(MS)),
intent(inout) :: ubt_trans
2775 real,
dimension(SZIW_(MS),SZJBW_(MS)),
intent(inout) :: vbt
2777 real,
dimension(SZIW_(MS),SZJBW_(MS)),
intent(inout) :: vhbt
2779 real,
dimension(SZIW_(MS),SZJBW_(MS)),
intent(inout) :: vbt_trans
2781 real,
dimension(SZIW_(MS),SZJW_(MS)),
intent(in) :: eta
2783 real,
dimension(SZIBW_(MS),SZJW_(MS)),
intent(in) :: ubt_old
2785 real,
dimension(SZIW_(MS),SZJBW_(MS)),
intent(in) :: vbt_old
2791 integer,
intent(in) :: halo
2792 real,
intent(in) :: dtbt
2793 real,
intent(in) :: bebt
2795 logical,
intent(in) :: use_BT_cont
2797 logical,
intent(in) :: integral_BT_cont
2800 real,
intent(in) :: dt_elapsed
2802 real,
dimension(SZIBW_(MS),SZJW_(MS)),
intent(in) :: Datu
2804 real,
dimension(SZIW_(MS),SZJBW_(MS)),
intent(in) :: Datv
2812 real,
dimension(SZIBW_(MS),SZJW_(MS)),
intent(in) :: uhbt0
2816 real,
dimension(SZIW_(MS),SZJBW_(MS)),
intent(in) :: vhbt0
2820 real,
dimension(SZIBW_(MS),SZJW_(MS)),
intent(in) :: ubt_int
2822 real,
dimension(SZIBW_(MS),SZJW_(MS)),
intent(in) :: uhbt_int
2824 real,
dimension(SZIW_(MS),SZJBW_(MS)),
intent(in) :: vbt_int
2826 real,
dimension(SZIW_(MS),SZJBW_(MS)),
intent(in) :: vhbt_int
2838 real :: uhbt_int_new
2839 real :: vhbt_int_new
2841 real :: cff, Cx, Cy, tau
2842 real :: dhdt, dhdx, dhdy
2844 integer :: i, j, is, ie, js, je
2845 real,
dimension(SZIB_(G),SZJB_(G)) :: grad
2846 real,
parameter :: eps = 1.0e-20
2847 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
2849 if (.not.(bt_obc%apply_u_OBCs .or. bt_obc%apply_v_OBCs))
return 2853 if (bt_obc%apply_u_OBCs)
then 2854 do j=js,je ;
do i=is-1,ie ;
if (obc%segnum_u(i,j) /= obc_none)
then 2855 if (obc%segment(obc%segnum_u(i,j))%specified)
then 2856 uhbt(i,j) = bt_obc%uhbt(i,j)
2857 ubt(i,j) = bt_obc%ubt_outer(i,j)
2858 vel_trans = ubt(i,j)
2859 elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e)
then 2860 if (obc%segment(obc%segnum_u(i,j))%Flather)
then 2861 cfl = dtbt * bt_obc%Cg_u(i,j) * g%IdxCu(i,j)
2862 u_inlet = cfl*ubt_old(i-1,j) + (1.0-cfl)*ubt_old(i,j)
2863 h_in = eta(i,j) + (0.5-cfl)*(eta(i,j)-eta(i-1,j))
2864 h_u = bt_obc%H_u(i,j)
2866 ubt(i,j) = 0.5*((u_inlet + bt_obc%ubt_outer(i,j)) + &
2867 (bt_obc%Cg_u(i,j)/h_u) * (h_in-bt_obc%eta_outer_u(i,j)))
2868 vel_trans = (1.0-bebt)*vel_prev + bebt*ubt(i,j)
2869 elseif (obc%segment(obc%segnum_u(i,j))%gradient)
then 2870 ubt(i,j) = ubt(i-1,j)
2871 vel_trans = ubt(i,j)
2873 elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w)
then 2874 if (obc%segment(obc%segnum_u(i,j))%Flather)
then 2875 cfl = dtbt * bt_obc%Cg_u(i,j) * g%IdxCu(i,j)
2876 u_inlet = cfl*ubt_old(i+1,j) + (1.0-cfl)*ubt_old(i,j)
2877 h_in = eta(i+1,j) + (0.5-cfl)*(eta(i+1,j)-eta(i+2,j))
2879 h_u = bt_obc%H_u(i,j)
2881 ubt(i,j) = 0.5*((u_inlet + bt_obc%ubt_outer(i,j)) + &
2882 (bt_obc%Cg_u(i,j)/h_u) * (bt_obc%eta_outer_u(i,j)-h_in))
2884 vel_trans = (1.0-bebt)*vel_prev + bebt*ubt(i,j)
2885 elseif (obc%segment(obc%segnum_u(i,j))%gradient)
then 2886 ubt(i,j) = ubt(i+1,j)
2887 vel_trans = ubt(i,j)
2891 if (.not. obc%segment(obc%segnum_u(i,j))%specified)
then 2892 if (integral_bt_cont)
then 2893 uhbt_int_new = find_uhbt(ubt_int(i,j) + dtbt*vel_trans, btcl_u(i,j)) + &
2894 dt_elapsed*uhbt0(i,j)
2895 uhbt(i,j) = (uhbt_int_new - uhbt_int(i,j)) * idtbt
2896 elseif (use_bt_cont)
then 2897 uhbt(i,j) = find_uhbt(vel_trans, btcl_u(i,j)) + uhbt0(i,j)
2899 uhbt(i,j) = datu(i,j)*vel_trans + uhbt0(i,j)
2903 ubt_trans(i,j) = vel_trans
2904 endif ;
enddo ;
enddo 2907 if (bt_obc%apply_v_OBCs)
then 2908 do j=js-1,je ;
do i=is,ie ;
if (obc%segnum_v(i,j) /= obc_none)
then 2909 if (obc%segment(obc%segnum_v(i,j))%specified)
then 2910 vhbt(i,j) = bt_obc%vhbt(i,j)
2911 vbt(i,j) = bt_obc%vbt_outer(i,j)
2912 vel_trans = vbt(i,j)
2913 elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n)
then 2914 if (obc%segment(obc%segnum_v(i,j))%Flather)
then 2915 cfl = dtbt * bt_obc%Cg_v(i,j) * g%IdyCv(i,j)
2916 v_inlet = cfl*vbt_old(i,j-1) + (1.0-cfl)*vbt_old(i,j)
2917 h_in = eta(i,j) + (0.5-cfl)*(eta(i,j)-eta(i,j-1))
2919 h_v = bt_obc%H_v(i,j)
2921 vbt(i,j) = 0.5*((v_inlet + bt_obc%vbt_outer(i,j)) + &
2922 (bt_obc%Cg_v(i,j)/h_v) * (h_in-bt_obc%eta_outer_v(i,j)))
2924 vel_trans = (1.0-bebt)*vel_prev + bebt*vbt(i,j)
2925 elseif (obc%segment(obc%segnum_v(i,j))%gradient)
then 2926 vbt(i,j) = vbt(i,j-1)
2927 vel_trans = vbt(i,j)
2929 elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s)
then 2930 if (obc%segment(obc%segnum_v(i,j))%Flather)
then 2931 cfl = dtbt * bt_obc%Cg_v(i,j) * g%IdyCv(i,j)
2932 v_inlet = cfl*vbt_old(i,j+1) + (1.0-cfl)*vbt_old(i,j)
2933 h_in = eta(i,j+1) + (0.5-cfl)*(eta(i,j+1)-eta(i,j+2))
2935 h_v = bt_obc%H_v(i,j)
2937 vbt(i,j) = 0.5*((v_inlet + bt_obc%vbt_outer(i,j)) + &
2938 (bt_obc%Cg_v(i,j)/h_v) * (bt_obc%eta_outer_v(i,j)-h_in))
2940 vel_trans = (1.0-bebt)*vel_prev + bebt*vbt(i,j)
2941 elseif (obc%segment(obc%segnum_v(i,j))%gradient)
then 2942 vbt(i,j) = vbt(i,j+1)
2943 vel_trans = vbt(i,j)
2947 if (.not. obc%segment(obc%segnum_v(i,j))%specified)
then 2948 if (integral_bt_cont)
then 2949 vhbt_int_new = find_vhbt(vbt_int(i,j) + dtbt*vel_trans, btcl_v(i,j)) + &
2950 dt_elapsed*vhbt0(i,j)
2951 vhbt(i,j) = (vhbt_int_new - vhbt_int(i,j)) * idtbt
2952 elseif (use_bt_cont)
then 2953 vhbt(i,j) = find_vhbt(vel_trans, btcl_v(i,j)) + vhbt0(i,j)
2955 vhbt(i,j) = vel_trans*datv(i,j) + vhbt0(i,j)
2959 vbt_trans(i,j) = vel_trans
2960 endif ;
enddo ;
enddo 2963 end subroutine apply_velocity_obcs
2967 subroutine set_up_bt_obc(OBC, eta, BT_OBC, BT_Domain, G, GV, US, MS, halo, use_BT_cont, &
2968 integral_BT_cont, dt_baroclinic, Datu, Datv, BTCL_u, BTCL_v)
2972 real,
dimension(SZIW_(MS),SZJW_(MS)),
intent(in) :: eta
2981 integer,
intent(in) :: halo
2982 logical,
intent(in) :: use_BT_cont
2984 logical,
intent(in) :: integral_BT_cont
2987 real,
intent(in) :: dt_baroclinic
2989 real,
dimension(SZIBW_(MS),SZJW_(MS)),
intent(in) :: Datu
2991 real,
dimension(SZIW_(MS),SZJBW_(MS)),
intent(in) :: Datv
3002 integer :: i, j, k, is, ie, js, je, n, nz, Isq, Ieq, Jsq, Jeq
3003 integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
3004 integer :: isdw, iedw, jsdw, jedw
3008 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
3009 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
3010 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
3011 isdw = ms%isdw ; iedw = ms%iedw ; jsdw = ms%jsdw ; jedw = ms%jedw
3013 i_dt = 1.0 / dt_baroclinic
3015 if ((isdw < isd) .or. (jsdw < jsd))
then 3016 call mom_error(fatal,
"set_up_BT_OBC: Open boundary conditions are not "//&
3017 "yet fully implemented with wide barotropic halos.")
3020 if (.not. bt_obc%is_alloced)
then 3021 allocate(bt_obc%Cg_u(isdw-1:iedw,jsdw:jedw)) ; bt_obc%Cg_u(:,:) = 0.0
3022 allocate(bt_obc%H_u(isdw-1:iedw,jsdw:jedw)) ; bt_obc%H_u(:,:) = 0.0
3023 allocate(bt_obc%uhbt(isdw-1:iedw,jsdw:jedw)) ; bt_obc%uhbt(:,:) = 0.0
3024 allocate(bt_obc%ubt_outer(isdw-1:iedw,jsdw:jedw)) ; bt_obc%ubt_outer(:,:) = 0.0
3025 allocate(bt_obc%eta_outer_u(isdw-1:iedw,jsdw:jedw)) ; bt_obc%eta_outer_u(:,:) = 0.0
3027 allocate(bt_obc%Cg_v(isdw:iedw,jsdw-1:jedw)) ; bt_obc%Cg_v(:,:) = 0.0
3028 allocate(bt_obc%H_v(isdw:iedw,jsdw-1:jedw)) ; bt_obc%H_v(:,:) = 0.0
3029 allocate(bt_obc%vhbt(isdw:iedw,jsdw-1:jedw)) ; bt_obc%vhbt(:,:) = 0.0
3030 allocate(bt_obc%vbt_outer(isdw:iedw,jsdw-1:jedw)) ; bt_obc%vbt_outer(:,:) = 0.0
3031 allocate(bt_obc%eta_outer_v(isdw:iedw,jsdw-1:jedw)) ; bt_obc%eta_outer_v(:,:)=0.0
3032 bt_obc%is_alloced = .true.
3033 call create_group_pass(bt_obc%pass_uv, bt_obc%ubt_outer, bt_obc%vbt_outer, bt_domain)
3035 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)
3036 call create_group_pass(bt_obc%pass_h, bt_obc%H_u, bt_obc%H_v, bt_domain,to_all+scalar_pair)
3037 call create_group_pass(bt_obc%pass_cg, bt_obc%Cg_u, bt_obc%Cg_v, bt_domain,to_all+scalar_pair)
3040 if (bt_obc%apply_u_OBCs)
then 3041 if (obc%specified_u_BCs_exist_globally)
then 3042 do n = 1, obc%number_of_segments
3043 segment => obc%segment(n)
3044 if (segment%is_E_or_W .and. segment%specified)
then 3045 do j=segment%HI%jsd,segment%HI%jed ;
do i=segment%HI%IsdB,segment%HI%IedB
3046 bt_obc%uhbt(i,j) = 0.
3048 do k=1,nz ;
do j=segment%HI%jsd,segment%HI%jed ;
do i=segment%HI%IsdB,segment%HI%IedB
3049 bt_obc%uhbt(i,j) = bt_obc%uhbt(i,j) + segment%normal_trans(i,j,k)
3050 enddo ;
enddo ;
enddo 3054 do j=js,je ;
do i=is-1,ie ;
if (obc%segnum_u(i,j) /= obc_none)
then 3056 if (obc%segment(obc%segnum_u(i,j))%specified)
then 3057 if (integral_bt_cont)
then 3058 bt_obc%ubt_outer(i,j) = uhbt_to_ubt(bt_obc%uhbt(i,j)*dt_baroclinic, btcl_u(i,j)) * i_dt
3059 elseif (use_bt_cont)
then 3060 bt_obc%ubt_outer(i,j) = uhbt_to_ubt(bt_obc%uhbt(i,j), btcl_u(i,j))
3062 if (datu(i,j) > 0.0) bt_obc%ubt_outer(i,j) = bt_obc%uhbt(i,j) / datu(i,j)
3065 if (gv%Boussinesq)
then 3066 if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e)
then 3067 bt_obc%H_u(i,j) = g%bathyT(i,j)*gv%Z_to_H + eta(i,j)
3068 elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w)
then 3069 bt_obc%H_u(i,j) = g%bathyT(i+1,j)*gv%Z_to_H + eta(i+1,j)
3072 if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e)
then 3073 bt_obc%H_u(i,j) = eta(i,j)
3074 elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w)
then 3075 bt_obc%H_u(i,j) = eta(i+1,j)
3078 bt_obc%Cg_u(i,j) = sqrt(gv%g_prime(1) * gv%H_to_Z*bt_obc%H_u(i,j))
3080 endif ;
enddo ;
enddo 3081 if (obc%Flather_u_BCs_exist_globally)
then 3082 do n = 1, obc%number_of_segments
3083 segment => obc%segment(n)
3084 if (segment%is_E_or_W .and. segment%Flather)
then 3085 do j=segment%HI%jsd,segment%HI%jed ;
do i=segment%HI%IsdB,segment%HI%IedB
3086 bt_obc%ubt_outer(i,j) = segment%normal_vel_bt(i,j)
3087 bt_obc%eta_outer_u(i,j) = segment%eta(i,j)
3094 if (bt_obc%apply_v_OBCs)
then 3095 if (obc%specified_v_BCs_exist_globally)
then 3096 do n = 1, obc%number_of_segments
3097 segment => obc%segment(n)
3098 if (segment%is_N_or_S .and. segment%specified)
then 3099 do j=segment%HI%JsdB,segment%HI%JedB ;
do i=segment%HI%isd,segment%HI%ied
3100 bt_obc%vhbt(i,j) = 0.
3102 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB ;
do i=segment%HI%isd,segment%HI%ied
3103 bt_obc%vhbt(i,j) = bt_obc%vhbt(i,j) + segment%normal_trans(i,j,k)
3104 enddo ;
enddo ;
enddo 3108 do j=js-1,je ;
do i=is,ie ;
if (obc%segnum_v(i,j) /= obc_none)
then 3110 if (obc%segment(obc%segnum_v(i,j))%specified)
then 3111 if (integral_bt_cont)
then 3112 bt_obc%vbt_outer(i,j) = vhbt_to_vbt(bt_obc%vhbt(i,j)*dt_baroclinic, btcl_v(i,j)) * i_dt
3113 elseif (use_bt_cont)
then 3114 bt_obc%vbt_outer(i,j) = vhbt_to_vbt(bt_obc%vhbt(i,j), btcl_v(i,j))
3116 if (datv(i,j) > 0.0) bt_obc%vbt_outer(i,j) = bt_obc%vhbt(i,j) / datv(i,j)
3119 if (gv%Boussinesq)
then 3120 if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n)
then 3121 bt_obc%H_v(i,j) = g%bathyT(i,j)*gv%Z_to_H + eta(i,j)
3122 elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s)
then 3123 bt_obc%H_v(i,j) = g%bathyT(i,j+1)*gv%Z_to_H + eta(i,j+1)
3126 if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n)
then 3127 bt_obc%H_v(i,j) = eta(i,j)
3128 elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s)
then 3129 bt_obc%H_v(i,j) = eta(i,j+1)
3132 bt_obc%Cg_v(i,j) = sqrt(gv%g_prime(1) * gv%H_to_Z*bt_obc%H_v(i,j))
3134 endif ;
enddo ;
enddo 3135 if (obc%Flather_v_BCs_exist_globally)
then 3136 do n = 1, obc%number_of_segments
3137 segment => obc%segment(n)
3138 if (segment%is_N_or_S .and. segment%Flather)
then 3139 do j=segment%HI%JsdB,segment%HI%JedB ;
do i=segment%HI%isd,segment%HI%ied
3140 bt_obc%vbt_outer(i,j) = segment%normal_vel_bt(i,j)
3141 bt_obc%eta_outer_v(i,j) = segment%eta(i,j)
3148 call do_group_pass(bt_obc%pass_uv, bt_domain)
3149 call do_group_pass(bt_obc%pass_uhvh, bt_domain)
3150 call do_group_pass(bt_obc%pass_eta_outer, bt_domain)
3151 call do_group_pass(bt_obc%pass_h, bt_domain)
3152 call do_group_pass(bt_obc%pass_cg, bt_domain)
3154 end subroutine set_up_bt_obc
3157 subroutine destroy_bt_obc(BT_OBC)
3162 if (bt_obc%is_alloced)
then 3163 deallocate(bt_obc%Cg_u)
3164 deallocate(bt_obc%H_u)
3165 deallocate(bt_obc%uhbt)
3166 deallocate(bt_obc%ubt_outer)
3167 deallocate(bt_obc%eta_outer_u)
3169 deallocate(bt_obc%Cg_v)
3170 deallocate(bt_obc%H_v)
3171 deallocate(bt_obc%vhbt)
3172 deallocate(bt_obc%vbt_outer)
3173 deallocate(bt_obc%eta_outer_v)
3174 bt_obc%is_alloced = .false.
3176 end subroutine destroy_bt_obc
3183 subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC)
3186 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
3190 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
3191 optional,
intent(in) :: h_u
3192 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
3193 optional,
intent(in) :: h_v
3194 logical,
optional,
intent(in) :: may_use_default
3203 real :: hatutot(szib_(g))
3204 real :: hatvtot(szi_(g))
3205 real :: ihatutot(szib_(g))
3206 real :: ihatvtot(szi_(g))
3214 real :: e_u(szib_(g),szk_(g)+1)
3215 real :: e_v(szi_(g),szk_(g)+1)
3216 real :: d_shallow_u(szi_(g))
3217 real :: d_shallow_v(szib_(g))
3221 logical :: use_default, test_dflt, apply_obcs
3222 integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz, i, j, k
3223 integer :: iss, ies, n
3227 if (.not.
associated(cs))
call mom_error(fatal, &
3228 "btcalc: Module MOM_barotropic must be initialized before it is used.")
3229 if (.not.cs%split)
return 3231 use_default = .false.
3232 test_dflt = .false. ;
if (
present(may_use_default)) test_dflt = may_use_default
3235 if (.not.((
present(h_u) .and.
present(h_v)) .or. &
3236 (cs%hvel_scheme == harmonic) .or. (cs%hvel_scheme == hybrid) .or.&
3237 (cs%hvel_scheme == arithmetic))) use_default = .true.
3239 if (.not.((
present(h_u) .and.
present(h_v)) .or. &
3240 (cs%hvel_scheme == harmonic) .or. (cs%hvel_scheme == hybrid) .or.&
3241 (cs%hvel_scheme == arithmetic)))
call mom_error(fatal, &
3242 "btcalc: Inconsistent settings of optional arguments and hvel_scheme.")
3245 apply_obcs = .false.
3246 if (
present(obc))
then ;
if (
associated(obc))
then ;
if (obc%OBC_pe)
then 3249 apply_obcs = (obc%number_of_segments > 0)
3250 endif ;
endif ;
endif 3252 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
3253 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
3254 h_neglect = gv%H_subroundoff
3262 if (
present(h_u))
then 3263 do i=is-1,ie ; hatutot(i) = h_u(i,j,1) ;
enddo 3264 do k=2,nz ;
do i=is-1,ie
3265 hatutot(i) = hatutot(i) + h_u(i,j,k)
3267 do i=is-1,ie ; ihatutot(i) = g%mask2dCu(i,j) / (hatutot(i) + h_neglect) ;
enddo 3268 do k=1,nz ;
do i=is-1,ie
3269 cs%frhatu(i,j,k) = h_u(i,j,k) * ihatutot(i)
3272 if (cs%hvel_scheme == arithmetic)
then 3274 cs%frhatu(i,j,1) = 0.5 * (h(i+1,j,1) + h(i,j,1))
3275 hatutot(i) = cs%frhatu(i,j,1)
3277 do k=2,nz ;
do i=is-1,ie
3278 cs%frhatu(i,j,k) = 0.5 * (h(i+1,j,k) + h(i,j,k))
3279 hatutot(i) = hatutot(i) + cs%frhatu(i,j,k)
3281 elseif (cs%hvel_scheme == hybrid .or. use_default)
then 3283 e_u(i,nz+1) = -0.5 * gv%Z_to_H * (g%bathyT(i+1,j) + g%bathyT(i,j))
3284 d_shallow_u(i) = -gv%Z_to_H * min(g%bathyT(i+1,j), g%bathyT(i,j))
3287 do k=nz,1,-1 ;
do i=is-1,ie
3288 e_u(i,k) = e_u(i,k+1) + 0.5 * (h(i+1,j,k) + h(i,j,k))
3289 h_arith = 0.5 * (h(i+1,j,k) + h(i,j,k))
3290 if (e_u(i,k+1) >= d_shallow_u(i))
then 3291 cs%frhatu(i,j,k) = h_arith
3293 h_harm = (h(i+1,j,k) * h(i,j,k)) / (h_arith + h_neglect)
3294 if (e_u(i,k) <= d_shallow_u(i))
then 3295 cs%frhatu(i,j,k) = h_harm
3297 wt_arith = (e_u(i,k) - d_shallow_u(i)) / (h_arith + h_neglect)
3298 cs%frhatu(i,j,k) = wt_arith*h_arith + (1.0-wt_arith)*h_harm
3301 hatutot(i) = hatutot(i) + cs%frhatu(i,j,k)
3303 elseif (cs%hvel_scheme == harmonic)
then 3305 cs%frhatu(i,j,1) = 2.0*(h(i+1,j,1) * h(i,j,1)) / &
3306 ((h(i+1,j,1) + h(i,j,1)) + h_neglect)
3307 hatutot(i) = cs%frhatu(i,j,1)
3309 do k=2,nz ;
do i=is-1,ie
3310 cs%frhatu(i,j,k) = 2.0*(h(i+1,j,k) * h(i,j,k)) / &
3311 ((h(i+1,j,k) + h(i,j,k)) + h_neglect)
3312 hatutot(i) = hatutot(i) + cs%frhatu(i,j,k)
3315 do i=is-1,ie ; ihatutot(i) = g%mask2dCu(i,j) / (hatutot(i) + h_neglect) ;
enddo 3316 do k=1,nz ;
do i=is-1,ie
3317 cs%frhatu(i,j,k) = cs%frhatu(i,j,k) * ihatutot(i)
3325 if (
present(h_v))
then 3326 do i=is,ie ; hatvtot(i) = h_v(i,j,1) ;
enddo 3327 do k=2,nz ;
do i=is,ie
3328 hatvtot(i) = hatvtot(i) + h_v(i,j,k)
3330 do i=is,ie ; ihatvtot(i) = g%mask2dCv(i,j) / (hatvtot(i) + h_neglect) ;
enddo 3331 do k=1,nz ;
do i=is,ie
3332 cs%frhatv(i,j,k) = h_v(i,j,k) * ihatvtot(i)
3335 if (cs%hvel_scheme == arithmetic)
then 3337 cs%frhatv(i,j,1) = 0.5 * (h(i,j+1,1) + h(i,j,1))
3338 hatvtot(i) = cs%frhatv(i,j,1)
3340 do k=2,nz ;
do i=is,ie
3341 cs%frhatv(i,j,k) = 0.5 * (h(i,j+1,k) + h(i,j,k))
3342 hatvtot(i) = hatvtot(i) + cs%frhatv(i,j,k)
3344 elseif (cs%hvel_scheme == hybrid .or. use_default)
then 3346 e_v(i,nz+1) = -0.5 * gv%Z_to_H * (g%bathyT(i,j+1) + g%bathyT(i,j))
3347 d_shallow_v(i) = -gv%Z_to_H * min(g%bathyT(i,j+1), g%bathyT(i,j))
3350 do k=nz,1,-1 ;
do i=is,ie
3351 e_v(i,k) = e_v(i,k+1) + 0.5 * (h(i,j+1,k) + h(i,j,k))
3352 h_arith = 0.5 * (h(i,j+1,k) + h(i,j,k))
3353 if (e_v(i,k+1) >= d_shallow_v(i))
then 3354 cs%frhatv(i,j,k) = h_arith
3356 h_harm = (h(i,j+1,k) * h(i,j,k)) / (h_arith + h_neglect)
3357 if (e_v(i,k) <= d_shallow_v(i))
then 3358 cs%frhatv(i,j,k) = h_harm
3360 wt_arith = (e_v(i,k) - d_shallow_v(i)) / (h_arith + h_neglect)
3361 cs%frhatv(i,j,k) = wt_arith*h_arith + (1.0-wt_arith)*h_harm
3364 hatvtot(i) = hatvtot(i) + cs%frhatv(i,j,k)
3366 elseif (cs%hvel_scheme == harmonic)
then 3368 cs%frhatv(i,j,1) = 2.0*(h(i,j+1,1) * h(i,j,1)) / &
3369 ((h(i,j+1,1) + h(i,j,1)) + h_neglect)
3370 hatvtot(i) = cs%frhatv(i,j,1)
3372 do k=2,nz ;
do i=is,ie
3373 cs%frhatv(i,j,k) = 2.0*(h(i,j+1,k) * h(i,j,k)) / &
3374 ((h(i,j+1,k) + h(i,j,k)) + h_neglect)
3375 hatvtot(i) = hatvtot(i) + cs%frhatv(i,j,k)
3378 do i=is,ie ; ihatvtot(i) = g%mask2dCv(i,j) / (hatvtot(i) + h_neglect) ;
enddo 3379 do k=1,nz ;
do i=is,ie
3380 cs%frhatv(i,j,k) = cs%frhatv(i,j,k) * ihatvtot(i)
3385 if (apply_obcs)
then ;
do n=1,obc%number_of_segments
3386 if (.not. obc%segment(n)%on_pe) cycle
3387 if (obc%segment(n)%direction == obc_direction_n)
then 3388 j = obc%segment(n)%HI%JsdB
3389 if ((j >= js-1) .and. (j <= je))
then 3390 iss = max(is,obc%segment(n)%HI%isd) ; ies = min(ie,obc%segment(n)%HI%ied)
3391 do i=iss,ies ; hatvtot(i) = h(i,j,1) ;
enddo 3392 do k=2,nz ;
do i=iss,ies
3393 hatvtot(i) = hatvtot(i) + h(i,j,k)
3396 ihatvtot(i) = g%mask2dCv(i,j) / (hatvtot(i) + h_neglect)
3398 do k=1,nz ;
do i=iss,ies
3399 cs%frhatv(i,j,k) = h(i,j,k) * ihatvtot(i)
3402 elseif (obc%segment(n)%direction == obc_direction_s)
then 3403 j = obc%segment(n)%HI%JsdB
3404 if ((j >= js-1) .and. (j <= je))
then 3405 iss = max(is,obc%segment(n)%HI%isd) ; ies = min(ie,obc%segment(n)%HI%ied)
3406 do i=iss,ies ; hatvtot(i) = h(i,j+1,1) ;
enddo 3407 do k=2,nz ;
do i=iss,ies
3408 hatvtot(i) = hatvtot(i) + h(i,j+1,k)
3411 ihatvtot(i) = g%mask2dCv(i,j) / (hatvtot(i) + h_neglect)
3413 do k=1,nz ;
do i=iss,ies
3414 cs%frhatv(i,j,k) = h(i,j+1,k) * ihatvtot(i)
3417 elseif (obc%segment(n)%direction == obc_direction_e)
then 3418 i = obc%segment(n)%HI%IsdB
3419 if ((i >= is-1) .and. (i <= ie))
then 3420 do j = max(js,obc%segment(n)%HI%jsd), min(je,obc%segment(n)%HI%jed)
3422 do k=2,nz ; htot = htot + h(i,j,k) ;
enddo 3423 ihtot = g%mask2dCu(i,j) / (htot + h_neglect)
3424 do k=1,nz ; cs%frhatu(i,j,k) = h(i,j,k) * ihtot ;
enddo 3427 elseif (obc%segment(n)%direction == obc_direction_w)
then 3428 i = obc%segment(n)%HI%IsdB
3429 if ((i >= is-1) .and. (i <= ie))
then 3430 do j = max(js,obc%segment(n)%HI%jsd), min(je,obc%segment(n)%HI%jed)
3432 do k=2,nz ; htot = htot + h(i+1,j,k) ;
enddo 3433 ihtot = g%mask2dCu(i,j) / (htot + h_neglect)
3434 do k=1,nz ; cs%frhatu(i,j,k) = h(i+1,j,k) * ihtot ;
enddo 3438 call mom_error(fatal,
"btcalc encountered and OBC segment of indeterminate direction.")
3443 call uvchksum(
"btcalc frhat[uv]", cs%frhatu, cs%frhatv, g%HI, &
3444 haloshift=0, symmetric=.true., omit_corners=.true., &
3446 if (
present(h_u) .and.
present(h_v)) &
3447 call uvchksum(
"btcalc h_[uv]", h_u, h_v, g%HI, haloshift=0, &
3448 symmetric=.true., omit_corners=.true., scale=gv%H_to_m, &
3450 call hchksum(h,
"btcalc h",g%HI, haloshift=1, scale=gv%H_to_m)
3453 end subroutine btcalc
3458 function find_uhbt(u, BTC)
result(uhbt)
3459 real,
intent(in) :: u
3469 elseif (u < btc%uBT_EE)
then 3470 uhbt = (u - btc%uBT_EE) * btc%FA_u_EE + btc%uh_EE
3471 elseif (u < 0.0)
then 3472 uhbt = u * (btc%FA_u_E0 + btc%uh_crvE * u**2)
3473 elseif (u <= btc%uBT_WW)
then 3474 uhbt = u * (btc%FA_u_W0 + btc%uh_crvW * u**2)
3476 uhbt = (u - btc%uBT_WW) * btc%FA_u_WW + btc%uh_WW
3479 end function find_uhbt
3483 function find_duhbt_du(u, BTC)
result(duhbt_du)
3484 real,
intent(in) :: u
3492 duhbt_du = 0.5*(btc%FA_u_E0 + btc%FA_u_W0)
3493 elseif (u < btc%uBT_EE)
then 3494 duhbt_du = btc%FA_u_EE
3495 elseif (u < 0.0)
then 3496 duhbt_du = (btc%FA_u_E0 + 3.0*btc%uh_crvE * u**2)
3497 elseif (u <= btc%uBT_WW)
then 3498 duhbt_du = (btc%FA_u_W0 + 3.0*btc%uh_crvW * u**2)
3500 duhbt_du = btc%FA_u_WW
3503 end function find_duhbt_du
3508 function uhbt_to_ubt(uhbt, BTC, guess)
result(ubt)
3509 real,
intent(in) :: uhbt
3516 real,
optional,
intent(in) :: guess
3522 real :: ubt_min, ubt_max
3525 real :: uherr_min, uherr_max
3527 real,
parameter :: tol = 1.0e-10
3530 real,
parameter :: vs1 = 1.25
3531 real,
parameter :: vs2 = 2.0
3533 integer :: itt, max_itt = 20
3536 if (uhbt == 0.0)
then 3538 elseif (uhbt < btc%uh_EE)
then 3539 ubt = btc%uBT_EE + (uhbt - btc%uh_EE) / btc%FA_u_EE
3540 elseif (uhbt < 0.0)
then 3543 ubt_min = btc%uBT_EE ; uherr_min = btc%uh_EE - uhbt
3544 ubt_max = 0.0 ; uherr_max = -uhbt
3546 ubt = btc%uBT_EE * (uhbt / btc%uh_EE)
3548 uhbt_err = ubt * (btc%FA_u_E0 + btc%uh_crvE * ubt**2) - uhbt
3550 if (abs(uhbt_err) < tol*abs(uhbt))
exit 3551 if (uhbt_err > 0.0)
then ; ubt_max = ubt ; uherr_max = uhbt_err ;
endif 3552 if (uhbt_err < 0.0)
then ; ubt_min = ubt ; uherr_min = uhbt_err ;
endif 3554 derr_du = btc%FA_u_E0 + 3.0 * btc%uh_crvE * ubt**2
3555 if ((uhbt_err >= derr_du*(ubt - ubt_min)) .or. &
3556 (-uhbt_err >= derr_du*(ubt_max - ubt)) .or. (derr_du <= 0.0))
then 3558 ubt = ubt_max + (ubt_min-ubt_max) * (uherr_max / (uherr_max-uherr_min))
3560 ubt = ubt - uhbt_err / derr_du
3561 if (abs(uhbt_err) < (0.01*tol)*abs(ubt_min*derr_du))
exit 3564 elseif (uhbt <= btc%uh_WW)
then 3566 ubt_min = 0.0 ; uherr_min = -uhbt
3567 ubt_max = btc%uBT_WW ; uherr_max = btc%uh_WW - uhbt
3569 ubt = btc%uBT_WW * (uhbt / btc%uh_WW)
3571 uhbt_err = ubt * (btc%FA_u_W0 + btc%uh_crvW * ubt**2) - uhbt
3573 if (abs(uhbt_err) < tol*abs(uhbt))
exit 3574 if (uhbt_err > 0.0)
then ; ubt_max = ubt ; uherr_max = uhbt_err ;
endif 3575 if (uhbt_err < 0.0)
then ; ubt_min = ubt ; uherr_min = uhbt_err ;
endif 3577 derr_du = btc%FA_u_W0 + 3.0 * btc%uh_crvW * ubt**2
3578 if ((uhbt_err >= derr_du*(ubt - ubt_min)) .or. &
3579 (-uhbt_err >= derr_du*(ubt_max - ubt)) .or. (derr_du <= 0.0))
then 3581 ubt = ubt_min + (ubt_max-ubt_min) * (-uherr_min / (uherr_max-uherr_min))
3583 ubt = ubt - uhbt_err / derr_du
3584 if (abs(uhbt_err) < (0.01*tol)*(ubt_max*derr_du))
exit 3588 ubt = btc%uBT_WW + (uhbt - btc%uh_WW) / btc%FA_u_WW
3591 if (
present(guess))
then 3592 dvel = abs(ubt) - vs1*abs(guess)
3593 if (dvel > 0.0)
then 3594 if (dvel < 40.0 * (abs(guess)*(vs2-vs1)) )
then 3595 vsr = vs2 - (vs2-vs1) * exp(-dvel / (abs(guess)*(vs2-vs1)))
3599 ubt = sign(vsr * guess, ubt)
3603 end function uhbt_to_ubt
3608 function find_vhbt(v, BTC)
result(vhbt)
3609 real,
intent(in) :: v
3618 elseif (v < btc%vBT_NN)
then 3619 vhbt = (v - btc%vBT_NN) * btc%FA_v_NN + btc%vh_NN
3620 elseif (v < 0.0)
then 3621 vhbt = v * (btc%FA_v_N0 + btc%vh_crvN * v**2)
3622 elseif (v <= btc%vBT_SS)
then 3623 vhbt = v * (btc%FA_v_S0 + btc%vh_crvS * v**2)
3625 vhbt = (v - btc%vBT_SS) * btc%FA_v_SS + btc%vh_SS
3628 end function find_vhbt
3632 function find_dvhbt_dv(v, BTC)
result(dvhbt_dv)
3633 real,
intent(in) :: v
3641 dvhbt_dv = 0.5*(btc%FA_v_N0 + btc%FA_v_S0)
3642 elseif (v < btc%vBT_NN)
then 3643 dvhbt_dv = btc%FA_v_NN
3644 elseif (v < 0.0)
then 3645 dvhbt_dv = btc%FA_v_N0 + 3.0*btc%vh_crvN * v**2
3646 elseif (v <= btc%vBT_SS)
then 3647 dvhbt_dv = btc%FA_v_S0 + 3.0*btc%vh_crvS * v**2
3649 dvhbt_dv = btc%FA_v_SS
3652 end function find_dvhbt_dv
3657 function vhbt_to_vbt(vhbt, BTC, guess)
result(vbt)
3658 real,
intent(in) :: vhbt
3665 real,
optional,
intent(in) :: guess
3671 real :: vbt_min, vbt_max
3674 real :: vherr_min, vherr_max
3676 real,
parameter :: tol = 1.0e-10
3679 real,
parameter :: vs1 = 1.25
3680 real,
parameter :: vs2 = 2.0
3682 integer :: itt, max_itt = 20
3685 if (vhbt == 0.0)
then 3687 elseif (vhbt < btc%vh_NN)
then 3688 vbt = btc%vBT_NN + (vhbt - btc%vh_NN) / btc%FA_v_NN
3689 elseif (vhbt < 0.0)
then 3692 vbt_min = btc%vBT_NN ; vherr_min = btc%vh_NN - vhbt
3693 vbt_max = 0.0 ; vherr_max = -vhbt
3695 vbt = btc%vBT_NN * (vhbt / btc%vh_NN)
3697 vhbt_err = vbt * (btc%FA_v_N0 + btc%vh_crvN * vbt**2) - vhbt
3699 if (abs(vhbt_err) < tol*abs(vhbt))
exit 3700 if (vhbt_err > 0.0)
then ; vbt_max = vbt ; vherr_max = vhbt_err ;
endif 3701 if (vhbt_err < 0.0)
then ; vbt_min = vbt ; vherr_min = vhbt_err ;
endif 3703 derr_dv = btc%FA_v_N0 + 3.0 * btc%vh_crvN * vbt**2
3704 if ((vhbt_err >= derr_dv*(vbt - vbt_min)) .or. &
3705 (-vhbt_err >= derr_dv*(vbt_max - vbt)) .or. (derr_dv <= 0.0))
then 3707 vbt = vbt_max + (vbt_min-vbt_max) * (vherr_max / (vherr_max-vherr_min))
3709 vbt = vbt - vhbt_err / derr_dv
3710 if (abs(vhbt_err) < (0.01*tol)*abs(derr_dv*vbt_min))
exit 3713 elseif (vhbt <= btc%vh_SS)
then 3715 vbt_min = 0.0 ; vherr_min = -vhbt
3716 vbt_max = btc%vBT_SS ; vherr_max = btc%vh_SS - vhbt
3718 vbt = btc%vBT_SS * (vhbt / btc%vh_SS)
3720 vhbt_err = vbt * (btc%FA_v_S0 + btc%vh_crvS * vbt**2) - vhbt
3722 if (abs(vhbt_err) < tol*abs(vhbt))
exit 3723 if (vhbt_err > 0.0)
then ; vbt_max = vbt ; vherr_max = vhbt_err ;
endif 3724 if (vhbt_err < 0.0)
then ; vbt_min = vbt ; vherr_min = vhbt_err ;
endif 3726 derr_dv = btc%FA_v_S0 + 3.0 * btc%vh_crvS * vbt**2
3727 if ((vhbt_err >= derr_dv*(vbt - vbt_min)) .or. &
3728 (-vhbt_err >= derr_dv*(vbt_max - vbt)) .or. (derr_dv <= 0.0))
then 3730 vbt = vbt_min + (vbt_max-vbt_min) * (-vherr_min / (vherr_max-vherr_min))
3732 vbt = vbt - vhbt_err / derr_dv
3733 if (abs(vhbt_err) < (0.01*tol)*(vbt_max*derr_dv))
exit 3737 vbt = btc%vBT_SS + (vhbt - btc%vh_SS) / btc%FA_v_SS
3740 if (
present(guess))
then 3741 dvel = abs(vbt) - vs1*abs(guess)
3742 if (dvel > 0.0)
then 3743 if (dvel < 40.0 * (abs(guess)*(vs2-vs1)) )
then 3744 vsr = vs2 - (vs2-vs1) * exp(-dvel / (abs(guess)*(vs2-vs1)))
3748 vbt = sign(guess * vsr, vbt)
3752 end function vhbt_to_vbt
3756 subroutine set_local_bt_cont_types(BT_cont, BTCL_u, BTCL_v, G, US, MS, BT_Domain, halo, dt_baroclinic)
3770 integer,
optional,
intent(in) :: halo
3771 real,
optional,
intent(in) :: dt_baroclinic
3776 real,
dimension(SZIBW_(MS),SZJW_(MS)) :: &
3779 FA_u_EE, FA_u_E0, FA_u_W0, FA_u_WW
3780 real,
dimension(SZIW_(MS),SZJBW_(MS)) :: &
3783 FA_v_NN, FA_v_N0, FA_v_S0, FA_v_SS
3785 real,
parameter :: C1_3 = 1.0/3.0
3786 integer :: i, j, is, ie, js, je, hs
3788 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
3789 hs = 1 ;
if (
present(halo)) hs = max(halo,0)
3790 dt = 1.0 ;
if (
present(dt_baroclinic)) dt = dt_baroclinic
3797 do j=js-hs,je+hs ;
do i=is-hs-1,ie+hs
3798 u_polarity(i,j) = 1.0
3799 ubt_ee(i,j) = 0.0 ; ubt_ww(i,j) = 0.0
3800 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
3803 do j=js-hs-1,je+hs ;
do i=is-hs,ie+hs
3804 v_polarity(i,j) = 1.0
3805 vbt_nn(i,j) = 0.0 ; vbt_ss(i,j) = 0.0
3806 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
3809 do j=js,je;
do i=is-1,ie
3810 ubt_ee(i,j) = bt_cont%uBT_EE(i,j) ; ubt_ww(i,j) = bt_cont%uBT_WW(i,j)
3811 fa_u_ee(i,j) = bt_cont%FA_u_EE(i,j) ; fa_u_e0(i,j) = bt_cont%FA_u_E0(i,j)
3812 fa_u_w0(i,j) = bt_cont%FA_u_W0(i,j) ; fa_u_ww(i,j) = bt_cont%FA_u_WW(i,j)
3815 do j=js-1,je;
do i=is,ie
3816 vbt_nn(i,j) = bt_cont%vBT_NN(i,j) ; vbt_ss(i,j) = bt_cont%vBT_SS(i,j)
3817 fa_v_nn(i,j) = bt_cont%FA_v_NN(i,j) ; fa_v_n0(i,j) = bt_cont%FA_v_N0(i,j)
3818 fa_v_s0(i,j) = bt_cont%FA_v_S0(i,j) ; fa_v_ss(i,j) = bt_cont%FA_v_SS(i,j)
3822 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
3823 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
3825 call create_group_pass(bt_cont%pass_polarity_BT, u_polarity, v_polarity, bt_domain)
3829 call create_group_pass(bt_cont%pass_FA_uv, fa_u_ee, fa_v_nn, bt_domain, to_all+scalar_pair)
3830 call create_group_pass(bt_cont%pass_FA_uv, fa_u_e0, fa_v_n0, bt_domain, to_all+scalar_pair)
3831 call create_group_pass(bt_cont%pass_FA_uv, fa_u_w0, fa_v_s0, bt_domain, to_all+scalar_pair)
3832 call create_group_pass(bt_cont%pass_FA_uv, fa_u_ww, fa_v_ss, bt_domain, to_all+scalar_pair)
3835 call do_group_pass(bt_cont%pass_polarity_BT, bt_domain)
3836 call do_group_pass(bt_cont%pass_FA_uv, bt_domain)
3837 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
3838 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
3842 do j=js-hs,je+hs ;
do i=is-hs-1,ie+hs
3843 btcl_u(i,j)%FA_u_EE = fa_u_ee(i,j) ; btcl_u(i,j)%FA_u_E0 = fa_u_e0(i,j)
3844 btcl_u(i,j)%FA_u_W0 = fa_u_w0(i,j) ; btcl_u(i,j)%FA_u_WW = fa_u_ww(i,j)
3845 btcl_u(i,j)%uBT_EE = dt*ubt_ee(i,j) ; btcl_u(i,j)%uBT_WW = dt*ubt_ww(i,j)
3847 if (u_polarity(i,j) < 0.0)
then 3848 call swap(btcl_u(i,j)%FA_u_EE, btcl_u(i,j)%FA_u_WW)
3849 call swap(btcl_u(i,j)%FA_u_E0, btcl_u(i,j)%FA_u_W0)
3850 call swap(btcl_u(i,j)%uBT_EE, btcl_u(i,j)%uBT_WW)
3853 btcl_u(i,j)%uh_EE = btcl_u(i,j)%uBT_EE * &
3854 (c1_3 * (2.0*btcl_u(i,j)%FA_u_E0 + btcl_u(i,j)%FA_u_EE))
3855 btcl_u(i,j)%uh_WW = btcl_u(i,j)%uBT_WW * &
3856 (c1_3 * (2.0*btcl_u(i,j)%FA_u_W0 + btcl_u(i,j)%FA_u_WW))
3858 btcl_u(i,j)%uh_crvE = 0.0 ; btcl_u(i,j)%uh_crvW = 0.0
3859 if (abs(btcl_u(i,j)%uBT_WW) > 0.0) btcl_u(i,j)%uh_crvW = &
3860 (c1_3 * (btcl_u(i,j)%FA_u_WW - btcl_u(i,j)%FA_u_W0)) / btcl_u(i,j)%uBT_WW**2
3861 if (abs(btcl_u(i,j)%uBT_EE) > 0.0) btcl_u(i,j)%uh_crvE = &
3862 (c1_3 * (btcl_u(i,j)%FA_u_EE - btcl_u(i,j)%FA_u_E0)) / btcl_u(i,j)%uBT_EE**2
3865 do j=js-hs-1,je+hs ;
do i=is-hs,ie+hs
3866 btcl_v(i,j)%FA_v_NN = fa_v_nn(i,j) ; btcl_v(i,j)%FA_v_N0 = fa_v_n0(i,j)
3867 btcl_v(i,j)%FA_v_S0 = fa_v_s0(i,j) ; btcl_v(i,j)%FA_v_SS = fa_v_ss(i,j)
3868 btcl_v(i,j)%vBT_NN = dt*vbt_nn(i,j) ; btcl_v(i,j)%vBT_SS = dt*vbt_ss(i,j)
3870 if (v_polarity(i,j) < 0.0)
then 3871 call swap(btcl_v(i,j)%FA_v_NN, btcl_v(i,j)%FA_v_SS)
3872 call swap(btcl_v(i,j)%FA_v_N0, btcl_v(i,j)%FA_v_S0)
3873 call swap(btcl_v(i,j)%vBT_NN, btcl_v(i,j)%vBT_SS)
3876 btcl_v(i,j)%vh_NN = btcl_v(i,j)%vBT_NN * &
3877 (c1_3 * (2.0*btcl_v(i,j)%FA_v_N0 + btcl_v(i,j)%FA_v_NN))
3878 btcl_v(i,j)%vh_SS = btcl_v(i,j)%vBT_SS * &
3879 (c1_3 * (2.0*btcl_v(i,j)%FA_v_S0 + btcl_v(i,j)%FA_v_SS))
3881 btcl_v(i,j)%vh_crvN = 0.0 ; btcl_v(i,j)%vh_crvS = 0.0
3882 if (abs(btcl_v(i,j)%vBT_SS) > 0.0) btcl_v(i,j)%vh_crvS = &
3883 (c1_3 * (btcl_v(i,j)%FA_v_SS - btcl_v(i,j)%FA_v_S0)) / btcl_v(i,j)%vBT_SS**2
3884 if (abs(btcl_v(i,j)%vBT_NN) > 0.0) btcl_v(i,j)%vh_crvN = &
3885 (c1_3 * (btcl_v(i,j)%FA_v_NN - btcl_v(i,j)%FA_v_N0)) / btcl_v(i,j)%vBT_NN**2
3888 end subroutine set_local_bt_cont_types
3895 subroutine adjust_local_bt_cont_types(ubt, uhbt, vbt, vhbt, BTCL_u, BTCL_v, &
3896 G, US, MS, halo, dt_baroclinic)
3898 real,
dimension(SZIBW_(MS),SZJW_(MS)), &
3900 real,
dimension(SZIBW_(MS),SZJW_(MS)), &
3903 real,
dimension(SZIW_(MS),SZJBW_(MS)), &
3905 real,
dimension(SZIW_(MS),SZJBW_(MS)), &
3909 intent(out) :: BTCL_u
3911 intent(out) :: BTCL_v
3914 integer,
optional,
intent(in) :: halo
3915 real,
optional,
intent(in) :: dt_baroclinic
3919 real,
dimension(SZIBW_(MS),SZJW_(MS)) :: &
3920 u_polarity, uBT_EE, uBT_WW, FA_u_EE, FA_u_E0, FA_u_W0, FA_u_WW
3921 real,
dimension(SZIW_(MS),SZJBW_(MS)) :: &
3922 v_polarity, vBT_NN, vBT_SS, FA_v_NN, FA_v_N0, FA_v_S0, FA_v_SS
3924 real,
parameter :: C1_3 = 1.0/3.0
3925 integer :: i, j, is, ie, js, je, hs
3927 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
3928 hs = 1 ;
if (
present(halo)) hs = max(halo,0)
3929 dt = 1.0 ;
if (
present(dt_baroclinic)) dt = dt_baroclinic
3932 do j=js-hs,je+hs ;
do i=is-hs-1,ie+hs
3933 if ((dt*ubt(i,j) > btcl_u(i,j)%uBT_WW) .and. (dt*uhbt(i,j) > btcl_u(i,j)%uh_WW))
then 3935 btcl_u(i,j)%ubt_WW = dt * ubt(i,j)
3936 if (3.0*uhbt(i,j) < 2.0*ubt(i,j) * btcl_u(i,j)%FA_u_W0)
then 3938 btcl_u(i,j)%uh_crvW = (uhbt(i,j) - ubt(i,j) * btcl_u(i,j)%FA_u_W0) / (dt**2 * ubt(i,j)**3)
3940 btcl_u(i,j)%FA_u_W0 = 1.5*uhbt(i,j) / ubt(i,j)
3941 btcl_u(i,j)%uh_crvW = -0.5*uhbt(i,j) / (dt**2 * ubt(i,j)**3)
3943 btcl_u(i,j)%uh_WW = dt * uhbt(i,j)
3946 elseif ((dt*ubt(i,j) < btcl_u(i,j)%uBT_EE) .and. (dt*uhbt(i,j) < btcl_u(i,j)%uh_EE))
then 3948 btcl_u(i,j)%ubt_EE = dt * ubt(i,j)
3949 if (3.0*uhbt(i,j) < 2.0*ubt(i,j) * btcl_u(i,j)%FA_u_E0)
then 3951 btcl_u(i,j)%uh_crvE = (uhbt(i,j) - ubt(i,j) * btcl_u(i,j)%FA_u_E0) / (dt**2 * ubt(i,j)**3)
3953 btcl_u(i,j)%FA_u_E0 = 1.5*uhbt(i,j) / ubt(i,j)
3954 btcl_u(i,j)%uh_crvE = -0.5*uhbt(i,j) / (dt**2 * ubt(i,j)**3)
3956 btcl_u(i,j)%uh_EE = dt * uhbt(i,j)
3962 do j=js-hs-1,je+hs ;
do i=is-hs,ie+hs
3963 if ((dt*vbt(i,j) > btcl_v(i,j)%vBT_SS) .and. (dt*vhbt(i,j) > btcl_v(i,j)%vh_SS))
then 3965 btcl_v(i,j)%vbt_SS = dt * vbt(i,j)
3966 if (3.0*vhbt(i,j) < 2.0*vbt(i,j) * btcl_v(i,j)%FA_v_S0)
then 3968 btcl_v(i,j)%vh_crvS = (vhbt(i,j) - vbt(i,j) * btcl_v(i,j)%FA_v_S0) / (dt**2 * vbt(i,j)**3)
3970 btcl_v(i,j)%FA_v_S0 = 1.5*vhbt(i,j) / (vbt(i,j))
3971 btcl_v(i,j)%vh_crvS = -0.5*vhbt(i,j) / (dt**2 * vbt(i,j)**3)
3973 btcl_v(i,j)%vh_SS = dt * vhbt(i,j)
3976 elseif ((dt*vbt(i,j) < btcl_v(i,j)%vBT_NN) .and. (dt*vhbt(i,j) < btcl_v(i,j)%vh_NN))
then 3978 btcl_v(i,j)%vbt_NN = dt * vbt(i,j)
3979 if (3.0*vhbt(i,j) < 2.0*vbt(i,j) * btcl_v(i,j)%FA_v_N0)
then 3981 btcl_v(i,j)%vh_crvN = (vhbt(i,j) - vbt(i,j) * btcl_v(i,j)%FA_v_N0) / (dt**2 * vbt(i,j)**3)
3983 btcl_v(i,j)%FA_v_N0 = 1.5*vhbt(i,j) / (vbt(i,j))
3984 btcl_v(i,j)%vh_crvN = -0.5*vhbt(i,j) / (dt**2 * vbt(i,j)**3)
3986 btcl_v(i,j)%vh_NN = dt * vhbt(i,j)
3992 end subroutine adjust_local_bt_cont_types
3996 subroutine bt_cont_to_face_areas(BT_cont, Datu, Datv, G, US, MS, halo, maximize)
4001 real,
dimension(MS%isdw-1:MS%iedw,MS%jsdw:MS%jedw), &
4003 real,
dimension(MS%isdw:MS%iedw,MS%jsdw-1:MS%jedw), &
4007 integer,
optional,
intent(in) :: halo
4008 logical,
optional,
intent(in) :: maximize
4013 integer :: i, j, is, ie, js, je, hs
4014 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
4015 hs = 1 ;
if (
present(halo)) hs = max(halo,0)
4016 find_max = .false. ;
if (
present(maximize)) find_max = maximize
4019 do j=js-hs,je+hs ;
do i=is-1-hs,ie+hs
4020 datu(i,j) = max(bt_cont%FA_u_EE(i,j), bt_cont%FA_u_E0(i,j), &
4021 bt_cont%FA_u_W0(i,j), bt_cont%FA_u_WW(i,j))
4023 do j=js-1-hs,je+hs ;
do i=is-hs,ie+hs
4024 datv(i,j) = max(bt_cont%FA_v_NN(i,j), bt_cont%FA_v_N0(i,j), &
4025 bt_cont%FA_v_S0(i,j), bt_cont%FA_v_SS(i,j))
4028 do j=js-hs,je+hs ;
do i=is-1-hs,ie+hs
4029 datu(i,j) = 0.5 * (bt_cont%FA_u_E0(i,j) + bt_cont%FA_u_W0(i,j))
4031 do j=js-1-hs,je+hs ;
do i=is-hs,ie+hs
4032 datv(i,j) = 0.5 * (bt_cont%FA_v_N0(i,j) + bt_cont%FA_v_S0(i,j))
4036 end subroutine bt_cont_to_face_areas
4039 subroutine swap(a,b)
4040 real,
intent(inout) :: a
4041 real,
intent(inout) :: b
4043 tmp = a ; a = b ; b = tmp
4048 subroutine find_face_areas(Datu, Datv, G, GV, US, CS, MS, eta, halo, add_max)
4050 real,
dimension(MS%isdw-1:MS%iedw,MS%jsdw:MS%jedw), &
4052 real,
dimension(MS%isdw:MS%iedw,MS%jsdw-1:MS%jedw), &
4059 real,
dimension(MS%isdw:MS%iedw,MS%jsdw:MS%jedw), &
4060 optional,
intent(in) :: eta
4062 integer,
optional,
intent(in) :: halo
4063 real,
optional,
intent(in) :: add_max
4068 integer :: i, j, is, ie, js, je, hs
4069 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
4070 hs = 1 ;
if (
present(halo)) hs = max(halo,0)
4074 if (
present(eta))
then 4076 if (gv%Boussinesq)
then 4078 do j=js-hs,je+hs ;
do i=is-1-hs,ie+hs
4079 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)
4080 datu(i,j) = 0.0 ;
if ((h1 > 0.0) .and. (h2 > 0.0)) &
4081 datu(i,j) = cs%dy_Cu(i,j) * (2.0 * h1 * h2) / (h1 + h2)
4085 do j=js-1-hs,je+hs ;
do i=is-hs,ie+hs
4086 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)
4087 datv(i,j) = 0.0 ;
if ((h1 > 0.0) .and. (h2 > 0.0)) &
4088 datv(i,j) = cs%dx_Cv(i,j) * (2.0 * h1 * h2) / (h1 + h2)
4093 do j=js-hs,je+hs ;
do i=is-1-hs,ie+hs
4094 datu(i,j) = 0.0 ;
if ((eta(i,j) > 0.0) .and. (eta(i+1,j) > 0.0)) &
4095 datu(i,j) = cs%dy_Cu(i,j) * (2.0 * eta(i,j) * eta(i+1,j)) / &
4096 (eta(i,j) + eta(i+1,j))
4100 do j=js-1-hs,je+hs ;
do i=is-hs,ie+hs
4101 datv(i,j) = 0.0 ;
if ((eta(i,j) > 0.0) .and. (eta(i,j+1) > 0.0)) &
4102 datv(i,j) = cs%dx_Cv(i,j) * (2.0 * eta(i,j) * eta(i,j+1)) / &
4103 (eta(i,j) + eta(i,j+1))
4107 elseif (
present(add_max))
then 4109 do j=js-hs,je+hs ;
do i=is-1-hs,ie+hs
4110 datu(i,j) = cs%dy_Cu(i,j) * gv%Z_to_H * &
4111 (max(cs%bathyT(i+1,j), cs%bathyT(i,j)) + add_max)
4114 do j=js-1-hs,je+hs ;
do i=is-hs,ie+hs
4115 datv(i,j) = cs%dx_Cv(i,j) * gv%Z_to_H * &
4116 (max(cs%bathyT(i,j+1), cs%bathyT(i,j)) + add_max)
4120 do j=js-hs,je+hs ;
do i=is-1-hs,ie+hs
4123 if (cs%bathyT(i+1,j)+cs%bathyT(i,j)>0.) &
4124 datu(i,j) = 2.0*cs%dy_Cu(i,j) * gv%Z_to_H * &
4125 (cs%bathyT(i+1,j) * cs%bathyT(i,j)) / &
4126 (cs%bathyT(i+1,j) + cs%bathyT(i,j))
4129 do j=js-1-hs,je+hs ;
do i=is-hs,ie+hs
4132 if (cs%bathyT(i,j+1)+cs%bathyT(i,j)>0.) &
4133 datv(i,j) = 2.0*cs%dx_Cv(i,j) * gv%Z_to_H * &
4134 (cs%bathyT(i,j+1) * cs%bathyT(i,j)) / &
4135 (cs%bathyT(i,j+1) + cs%bathyT(i,j))
4140 end subroutine find_face_areas
4146 subroutine bt_mass_source(h, eta, set_cor, G, GV, CS)
4149 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
4150 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: eta
4152 logical,
intent(in) :: set_cor
4160 real :: h_tot(szi_(g))
4161 real :: eta_h(szi_(g))
4165 integer :: is, ie, js, je, nz, i, j, k
4167 if (.not.
associated(cs))
call mom_error(fatal,
"bt_mass_source: "// &
4168 "Module MOM_barotropic must be initialized before it is used.")
4169 if (.not.cs%split)
return 4171 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
4175 do i=is,ie ; h_tot(i) = h(i,j,1) ;
enddo 4176 if (gv%Boussinesq)
then 4177 do i=is,ie ; eta_h(i) = h(i,j,1) - g%bathyT(i,j)*gv%Z_to_H ;
enddo 4179 do i=is,ie ; eta_h(i) = h(i,j,1) ;
enddo 4181 do k=2,nz ;
do i=is,ie
4182 eta_h(i) = eta_h(i) + h(i,j,k)
4183 h_tot(i) = h_tot(i) + h(i,j,k)
4188 d_eta = eta_h(i) - eta(i,j)
4189 cs%eta_cor(i,j) = d_eta
4193 d_eta = eta_h(i) - eta(i,j)
4194 cs%eta_cor(i,j) = cs%eta_cor(i,j) + d_eta
4199 end subroutine bt_mass_source
4204 subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, &
4205 restart_CS, calc_dtbt, BT_cont, tides_CSp)
4209 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
4211 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
4213 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
4215 real,
dimension(SZI_(G),SZJ_(G)), &
4218 type(time_type),
target,
intent(in) :: time
4220 type(
diag_ctrl),
target,
intent(inout) :: diag
4225 logical,
intent(out) :: calc_dtbt
4232 pointer :: tides_csp
4236 #include "version_variable.h" 4238 character(len=40) :: mdl =
"MOM_barotropic" 4239 real :: datu(szibs_(g),szj_(g))
4240 real :: datv(szi_(g),szjbs_(g))
4241 real :: gtot_estimate
4246 real :: wave_drag_scale
4248 character(len=200) :: inputdir
4249 character(len=200) :: wave_drag_file
4251 character(len=80) :: wave_drag_var
4259 real,
allocatable,
dimension(:,:) :: lin_drag_h
4261 type(group_pass_type) :: pass_static_data, pass_q_d_cor
4262 type(group_pass_type) :: pass_bt_hbt_btav, pass_a_polarity
4263 logical :: default_2018_answers
4264 logical :: apply_bt_drag, use_bt_cont_type
4265 character(len=48) :: thickness_units, flux_units
4266 character*(40) :: hvel_str
4267 integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz
4268 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
4269 integer :: isdw, iedw, jsdw, jedw
4271 integer :: wd_halos(2), bt_halo_sz
4272 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
4273 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
4274 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
4275 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
4276 ms%isdw = g%isd ; ms%iedw = g%ied ; ms%jsdw = g%jsd ; ms%jedw = g%jed
4278 if (cs%module_is_initialized)
then 4279 call mom_error(warning,
"barotropic_init called with a control structure "// &
4280 "that has already been initialized.")
4283 cs%module_is_initialized = .true.
4285 cs%diag => diag ; cs%Time => time
4286 if (
present(tides_csp))
then 4287 if (
associated(tides_csp)) cs%tides_CSp => tides_csp
4291 call get_param(param_file, mdl,
"SPLIT", cs%split, default=.true., do_not_log=.true.)
4292 call log_version(param_file, mdl, version,
"", log_to_all=.true., layout=cs%split, &
4293 debugging=cs%split, all_default=.not.cs%split)
4294 call get_param(param_file, mdl,
"SPLIT", cs%split, &
4295 "Use the split time stepping if true.", default=.true.)
4296 if (.not.cs%split)
return 4298 call get_param(param_file, mdl,
"USE_BT_CONT_TYPE", use_bt_cont_type, &
4299 "If true, use a structure with elements that describe "//&
4300 "effective face areas from the summed continuity solver "//&
4301 "as a function the barotropic flow in coupling between "//&
4302 "the barotropic and baroclinic flow. This is only used "//&
4303 "if SPLIT is true.", default=.true.)
4304 call get_param(param_file, mdl,
"INTEGRAL_BT_CONTINUITY", cs%integral_bt_cont, &
4305 "If true, use the time-integrated velocity over the barotropic steps "//&
4306 "to determine the integrated transports used to update the continuity "//&
4307 "equation. Otherwise the transports are the sum of the transports based on "//&
4308 "a series of instantaneous velocities and the BT_CONT_TYPE for transports. "//&
4309 "This is only valid if USE_BT_CONT_TYPE = True.", &
4310 default=.false., do_not_log=.not.use_bt_cont_type)
4311 call get_param(param_file, mdl,
"BOUND_BT_CORRECTION", cs%bound_BT_corr, &
4312 "If true, the corrective pseudo mass-fluxes into the "//&
4313 "barotropic solver are limited to values that require "//&
4314 "less than maxCFL_BT_cont to be accommodated.",default=.false.)
4315 call get_param(param_file, mdl,
"BT_CONT_CORR_BOUNDS", cs%BT_cont_bounds, &
4316 "If true, and BOUND_BT_CORRECTION is true, use the "//&
4317 "BT_cont_type variables to set limits determined by "//&
4318 "MAXCFL_BT_CONT on the CFL number of the velocities "//&
4319 "that are likely to be driven by the corrective mass fluxes.", &
4320 default=.true., do_not_log=.not.cs%bound_BT_corr)
4321 call get_param(param_file, mdl,
"ADJUST_BT_CONT", cs%adjust_BT_cont, &
4322 "If true, adjust the curve fit to the BT_cont type "//&
4323 "that is used by the barotropic solver to match the "//&
4324 "transport about which the flow is being linearized.", &
4325 default=.false., do_not_log=.not.use_bt_cont_type)
4326 call get_param(param_file, mdl,
"GRADUAL_BT_ICS", cs%gradual_BT_ICs, &
4327 "If true, adjust the initial conditions for the "//&
4328 "barotropic solver to the values from the layered "//&
4329 "solution over a whole timestep instead of instantly. "//&
4330 "This is a decent approximation to the inclusion of "//&
4331 "sum(u dh_dt) while also correcting for truncation errors.", &
4333 call get_param(param_file, mdl,
"BT_USE_VISC_REM_U_UH0", cs%visc_rem_u_uh0, &
4334 "If true, use the viscous remnants when estimating the "//&
4335 "barotropic velocities that were used to calculate uh0 "//&
4336 "and vh0. False is probably the better choice.", default=.false.)
4337 call get_param(param_file, mdl,
"BT_USE_WIDE_HALOS", cs%use_wide_halos, &
4338 "If true, use wide halos and march in during the "//&
4339 "barotropic time stepping for efficiency.", default=.true., &
4341 call get_param(param_file, mdl,
"BTHALO", bt_halo_sz, &
4342 "The minimum halo size for the barotropic solver.", default=0, &
4344 #ifdef STATIC_MEMORY_ 4345 if ((bt_halo_sz > 0) .and. (bt_halo_sz /= bthalo_))
call mom_error(fatal, &
4346 "barotropic_init: Run-time values of BTHALO must agree with the "//&
4347 "macro BTHALO_ with STATIC_MEMORY_.")
4348 wd_halos(1) = whaloi_+nihalo_ ; wd_halos(2) = whaloj_+njhalo_
4350 wd_halos(1) = bt_halo_sz; wd_halos(2) = bt_halo_sz
4352 call log_param(param_file, mdl,
"!BT x-halo", wd_halos(1), &
4353 "The barotropic x-halo size that is actually used.", &
4355 call log_param(param_file, mdl,
"!BT y-halo", wd_halos(2), &
4356 "The barotropic y-halo size that is actually used.", &
4359 call get_param(param_file, mdl,
"NONLINEAR_BT_CONTINUITY", cs%Nonlinear_continuity, &
4360 "If true, use nonlinear transports in the barotropic "//&
4361 "continuity equation. This does not apply if "//&
4362 "USE_BT_CONT_TYPE is true.", default=.false., do_not_log=use_bt_cont_type)
4363 call get_param(param_file, mdl,
"NONLIN_BT_CONT_UPDATE_PERIOD", cs%Nonlin_cont_update_period, &
4364 "If NONLINEAR_BT_CONTINUITY is true, this is the number "//&
4365 "of barotropic time steps between updates to the face "//&
4366 "areas, or 0 to update only before the barotropic stepping.", &
4367 units=
"nondim", default=1, do_not_log=.not.cs%Nonlinear_continuity)
4369 call get_param(param_file, mdl,
"BT_PROJECT_VELOCITY", cs%BT_project_velocity,&
4370 "If true, step the barotropic velocity first and project "//&
4371 "out the velocity tendency by 1+BEBT when calculating the "//&
4372 "transport. The default (false) is to use a predictor "//&
4373 "continuity step to find the pressure field, and then "//&
4374 "to do a corrector continuity step using a weighted "//&
4375 "average of the old and new velocities, with weights "//&
4376 "of (1-BEBT) and BEBT.", default=.false.)
4377 call get_param(param_file, mdl,
"BT_NONLIN_STRESS", cs%nonlin_stress, &
4378 "If true, use the full depth of the ocean at the start of the barotropic "//&
4379 "step when calculating the surface stress contribution to the barotropic "//&
4380 "acclerations. Otherwise use the depth based on bathyT.", default=.false.)
4382 call get_param(param_file, mdl,
"DYNAMIC_SURFACE_PRESSURE", cs%dynamic_psurf, &
4383 "If true, add a dynamic pressure due to a viscous ice "//&
4384 "shelf, for instance.", default=.false.)
4385 call get_param(param_file, mdl,
"ICE_LENGTH_DYN_PSURF", cs%ice_strength_length, &
4386 "The length scale at which the Rayleigh damping rate due "//&
4387 "to the ice strength should be the same as if a Laplacian "//&
4388 "were applied, if DYNAMIC_SURFACE_PRESSURE is true.", &
4389 units=
"m", default=1.0e4, scale=us%m_to_L, do_not_log=.not.cs%dynamic_psurf)
4390 call get_param(param_file, mdl,
"DEPTH_MIN_DYN_PSURF", cs%Dmin_dyn_psurf, &
4391 "The minimum depth to use in limiting the size of the "//&
4392 "dynamic surface pressure for stability, if "//&
4393 "DYNAMIC_SURFACE_PRESSURE is true..", &
4394 units=
"m", default=1.0e-6, scale=us%m_to_Z, do_not_log=.not.cs%dynamic_psurf)
4395 call get_param(param_file, mdl,
"CONST_DYN_PSURF", cs%const_dyn_psurf, &
4396 "The constant that scales the dynamic surface pressure, "//&
4397 "if DYNAMIC_SURFACE_PRESSURE is true. Stable values "//&
4398 "are < ~1.0.", units=
"nondim", default=0.9, do_not_log=.not.cs%dynamic_psurf)
4400 call get_param(param_file, mdl,
"BT_CORIOLIS_SCALE", cs%BT_Coriolis_scale, &
4401 "A factor by which the barotropic Coriolis anomaly terms are scaled.", &
4402 units=
"nondim", default=1.0)
4403 call get_param(param_file, mdl,
"DEFAULT_2018_ANSWERS", default_2018_answers, &
4404 "This sets the default value for the various _2018_ANSWERS parameters.", &
4406 call get_param(param_file, mdl,
"BAROTROPIC_2018_ANSWERS", cs%answers_2018, &
4407 "If true, use expressions for the barotropic solver that recover the answers "//&
4408 "from the end of 2018. Otherwise, use more efficient or general expressions.", &
4409 default=default_2018_answers)
4411 call get_param(param_file, mdl,
"TIDES", cs%tides, &
4412 "If true, apply tidal momentum forcing.", default=.false.)
4413 call get_param(param_file, mdl,
"SADOURNY", cs%Sadourny, &
4414 "If true, the Coriolis terms are discretized with the "//&
4415 "Sadourny (1975) energy conserving scheme, otherwise "//&
4416 "the Arakawa & Hsu scheme is used. If the internal "//&
4417 "deformation radius is not resolved, the Sadourny scheme "//&
4418 "should probably be used.", default=.true.)
4420 call get_param(param_file, mdl,
"BT_THICK_SCHEME", hvel_str, &
4421 "A string describing the scheme that is used to set the "//&
4422 "open face areas used for barotropic transport and the "//&
4423 "relative weights of the accelerations. Valid values are:\n"//&
4424 "\t ARITHMETIC - arithmetic mean layer thicknesses \n"//&
4425 "\t HARMONIC - harmonic mean layer thicknesses \n"//&
4426 "\t HYBRID (the default) - use arithmetic means for \n"//&
4427 "\t layers above the shallowest bottom, the harmonic \n"//&
4428 "\t mean for layers below, and a weighted average for \n"//&
4429 "\t layers that straddle that depth \n"//&
4430 "\t FROM_BT_CONT - use the average thicknesses kept \n"//&
4431 "\t in the h_u and h_v fields of the BT_cont_type", &
4432 default=bt_cont_string)
4433 select case (hvel_str)
4434 case (hybrid_string) ; cs%hvel_scheme = hybrid
4435 case (harmonic_string) ; cs%hvel_scheme = harmonic
4436 case (arithmetic_string) ; cs%hvel_scheme = arithmetic
4437 case (bt_cont_string) ; cs%hvel_scheme = from_bt_cont
4439 call mom_mesg(
'barotropic_init: BT_THICK_SCHEME ="'//trim(hvel_str)//
'"', 0)
4440 call mom_error(fatal,
"barotropic_init: Unrecognized setting "// &
4441 "#define BT_THICK_SCHEME "//trim(hvel_str)//
" found in input file.")
4443 if ((cs%hvel_scheme == from_bt_cont) .and. .not.use_bt_cont_type) &
4444 call mom_error(fatal,
"barotropic_init: BT_THICK_SCHEME FROM_BT_CONT "//&
4445 "can only be used if USE_BT_CONT_TYPE is defined.")
4447 call get_param(param_file, mdl,
"BT_STRONG_DRAG", cs%strong_drag, &
4448 "If true, use a stronger estimate of the retarding "//&
4449 "effects of strong bottom drag, by making it implicit "//&
4450 "with the barotropic time-step instead of implicit with "//&
4451 "the baroclinic time-step and dividing by the number of "//&
4452 "barotropic steps.", default=.false.)
4453 call get_param(param_file, mdl,
"BT_LINEAR_WAVE_DRAG", cs%linear_wave_drag, &
4454 "If true, apply a linear drag to the barotropic velocities, "//&
4455 "using rates set by lin_drag_u & _v divided by the depth of "//&
4456 "the ocean. This was introduced to facilitate tide modeling.", &
4458 call get_param(param_file, mdl,
"BT_WAVE_DRAG_FILE", wave_drag_file, &
4459 "The name of the file with the barotropic linear wave drag "//&
4460 "piston velocities.", default=
"", do_not_log=.not.cs%linear_wave_drag)
4461 call get_param(param_file, mdl,
"BT_WAVE_DRAG_VAR", wave_drag_var, &
4462 "The name of the variable in BT_WAVE_DRAG_FILE with the "//&
4463 "barotropic linear wave drag piston velocities at h points.", &
4464 default=
"rH", do_not_log=.not.cs%linear_wave_drag)
4465 call get_param(param_file, mdl,
"BT_WAVE_DRAG_SCALE", wave_drag_scale, &
4466 "A scaling factor for the barotropic linear wave drag "//&
4467 "piston velocities.", default=1.0, units=
"nondim", &
4468 do_not_log=.not.cs%linear_wave_drag)
4470 call get_param(param_file, mdl,
"CLIP_BT_VELOCITY", cs%clip_velocity, &
4471 "If true, limit any velocity components that exceed "//&
4472 "CFL_TRUNCATE. This should only be used as a desperate "//&
4473 "debugging measure.", default=.false.)
4474 call get_param(param_file, mdl,
"CFL_TRUNCATE", cs%CFL_trunc, &
4475 "The value of the CFL number that will cause velocity "//&
4476 "components to be truncated; instability can occur past 0.5.", &
4477 units=
"nondim", default=0.5, do_not_log=.not.cs%clip_velocity)
4478 call get_param(param_file, mdl,
"MAXVEL", cs%maxvel, &
4479 "The maximum velocity allowed before the velocity "//&
4480 "components are truncated.", units=
"m s-1", default=3.0e8, scale=us%m_s_to_L_T, &
4481 do_not_log=.not.cs%clip_velocity)
4482 call get_param(param_file, mdl,
"MAXCFL_BT_CONT", cs%maxCFL_BT_cont, &
4483 "The maximum permitted CFL number associated with the "//&
4484 "barotropic accelerations from the summed velocities "//&
4485 "times the time-derivatives of thicknesses.", units=
"nondim", &
4487 call get_param(param_file, mdl,
"VEL_UNDERFLOW", cs%vel_underflow, &
4488 "A negligibly small velocity magnitude below which velocity "//&
4489 "components are set to 0. A reasonable value might be "//&
4490 "1e-30 m/s, which is less than an Angstrom divided by "//&
4491 "the age of the universe.", units=
"m s-1", default=0.0, scale=us%m_s_to_L_T)
4493 call get_param(param_file, mdl,
"DT_BT_FILTER", cs%dt_bt_filter, &
4494 "A time-scale over which the barotropic mode solutions "//&
4495 "are filtered, in seconds if positive, or as a fraction "//&
4496 "of DT if negative. When used this can never be taken to "//&
4497 "be longer than 2*dt. Set this to 0 to apply no filtering.", &
4498 units=
"sec or nondim", default=-0.25)
4499 if (cs%dt_bt_filter > 0.0) cs%dt_bt_filter = us%s_to_T*cs%dt_bt_filter
4500 call get_param(param_file, mdl,
"G_BT_EXTRA", cs%G_extra, &
4501 "A nondimensional factor by which gtot is enhanced.", &
4502 units=
"nondim", default=0.0)
4503 call get_param(param_file, mdl,
"SSH_EXTRA", ssh_extra, &
4504 "An estimate of how much higher SSH might get, for use "//&
4505 "in calculating the safe external wave speed. The "//&
4506 "default is the minimum of 10 m or 5% of MAXIMUM_DEPTH.", &
4507 units=
"m", default=min(10.0,0.05*g%max_depth*us%Z_to_m), scale=us%m_to_Z)
4509 call get_param(param_file, mdl,
"DEBUG", cs%debug, &
4510 "If true, write out verbose debugging data.", &
4511 default=.false., debuggingparam=.true.)
4512 call get_param(param_file, mdl,
"DEBUG_BT", cs%debug_bt, &
4513 "If true, write out verbose debugging data within the "//&
4514 "barotropic time-stepping loop. The data volume can be "//&
4515 "quite large if this is true.", default=cs%debug, &
4516 debuggingparam=.true.)
4518 call get_param(param_file, mdl,
"LINEARIZED_BT_CORIOLIS", cs%linearized_BT_PV, &
4519 "If true use the bottom depth instead of the total water column thickness "//&
4520 "in the barotropic Coriolis term calculations.", default=.true.)
4521 call get_param(param_file, mdl,
"BEBT", cs%bebt, &
4522 "BEBT determines whether the barotropic time stepping "//&
4523 "uses the forward-backward time-stepping scheme or a "//&
4524 "backward Euler scheme. BEBT is valid in the range from "//&
4525 "0 (for a forward-backward treatment of nonrotating "//&
4526 "gravity waves) to 1 (for a backward Euler treatment). "//&
4527 "In practice, BEBT must be greater than about 0.05.", &
4528 units=
"nondim", default=0.1)
4529 call get_param(param_file, mdl,
"DTBT", dtbt_input, &
4530 "The barotropic time step, in s. DTBT is only used with "//&
4531 "the split explicit time stepping. To set the time step "//&
4532 "automatically based the maximum stable value use 0, or "//&
4533 "a negative value gives the fraction of the stable value. "//&
4534 "Setting DTBT to 0 is the same as setting it to -0.98. "//&
4535 "The value of DTBT that will actually be used is an "//&
4536 "integer fraction of DT, rounding down.", units=
"s or nondim",&
4538 call get_param(param_file, mdl,
"BT_USE_OLD_CORIOLIS_BRACKET_BUG", &
4539 cs%use_old_coriolis_bracket_bug , &
4540 "If True, use an order of operations that is not bitwise "//&
4541 "rotationally symmetric in the meridional Coriolis term of "//&
4542 "the barotropic solver.", default=.false.)
4545 call clone_mom_domain(g%Domain, cs%BT_Domain, min_halo=wd_halos, symmetric=.true.)
4546 #ifdef STATIC_MEMORY_ 4547 if (wd_halos(1) /= whaloi_+nihalo_)
call mom_error(fatal,
"barotropic_init: "//&
4548 "Barotropic x-halo sizes are incorrectly resized with STATIC_MEMORY_.")
4549 if (wd_halos(2) /= whaloj_+njhalo_)
call mom_error(fatal,
"barotropic_init: "//&
4550 "Barotropic y-halo sizes are incorrectly resized with STATIC_MEMORY_.")
4552 if (bt_halo_sz > 0)
then 4553 if (wd_halos(1) > bt_halo_sz) &
4554 call mom_mesg(
"barotropic_init: barotropic x-halo size increased.", 3)
4555 if (wd_halos(2) > bt_halo_sz) &
4556 call mom_mesg(
"barotropic_init: barotropic y-halo size increased.", 3)
4560 cs%isdw = g%isc-wd_halos(1) ; cs%iedw = g%iec+wd_halos(1)
4561 cs%jsdw = g%jsc-wd_halos(2) ; cs%jedw = g%jec+wd_halos(2)
4562 isdw = cs%isdw ; iedw = cs%iedw ; jsdw = cs%jsdw ; jedw = cs%jedw
4564 alloc_(cs%frhatu(isdb:iedb,jsd:jed,nz)) ; alloc_(cs%frhatv(isd:ied,jsdb:jedb,nz))
4565 alloc_(cs%eta_cor(isd:ied,jsd:jed))
4566 if (cs%bound_BT_corr)
then 4567 alloc_(cs%eta_cor_bound(isd:ied,jsd:jed)) ; cs%eta_cor_bound(:,:) = 0.0
4569 alloc_(cs%IDatu(isdb:iedb,jsd:jed)) ; alloc_(cs%IDatv(isd:ied,jsdb:jedb))
4571 alloc_(cs%ua_polarity(isdw:iedw,jsdw:jedw))
4572 alloc_(cs%va_polarity(isdw:iedw,jsdw:jedw))
4574 cs%frhatu(:,:,:) = 0.0 ; cs%frhatv(:,:,:) = 0.0
4575 cs%eta_cor(:,:) = 0.0
4576 cs%IDatu(:,:) = 0.0 ; cs%IDatv(:,:) = 0.0
4578 cs%ua_polarity(:,:) = 1.0 ; cs%va_polarity(:,:) = 1.0
4579 call create_group_pass(pass_a_polarity, cs%ua_polarity, cs%va_polarity, cs%BT_domain, to_all, agrid)
4580 call do_group_pass(pass_a_polarity, cs%BT_domain)
4582 if (use_bt_cont_type) &
4583 call alloc_bt_cont_type(bt_cont, g, (cs%hvel_scheme == from_bt_cont))
4586 allocate(cs%debug_BT_HI)
4587 cs%debug_BT_HI%isc=g%isc
4588 cs%debug_BT_HI%iec=g%iec
4589 cs%debug_BT_HI%jsc=g%jsc
4590 cs%debug_BT_HI%jec=g%jec
4591 cs%debug_BT_HI%IscB=g%isc-1
4592 cs%debug_BT_HI%IecB=g%iec
4593 cs%debug_BT_HI%JscB=g%jsc-1
4594 cs%debug_BT_HI%JecB=g%jec
4595 cs%debug_BT_HI%isd=cs%isdw
4596 cs%debug_BT_HI%ied=cs%iedw
4597 cs%debug_BT_HI%jsd=cs%jsdw
4598 cs%debug_BT_HI%jed=cs%jedw
4599 cs%debug_BT_HI%IsdB=cs%isdw-1
4600 cs%debug_BT_HI%IedB=cs%iedw
4601 cs%debug_BT_HI%JsdB=cs%jsdw-1
4602 cs%debug_BT_HI%JedB=cs%jedw
4603 cs%debug_BT_HI%turns = g%HI%turns
4607 alloc_(cs%IareaT(cs%isdw:cs%iedw,cs%jsdw:cs%jedw)) ; cs%IareaT(:,:) = 0.0
4608 alloc_(cs%bathyT(cs%isdw:cs%iedw,cs%jsdw:cs%jedw)) ; cs%bathyT(:,:) = gv%Angstrom_m
4609 alloc_(cs%IdxCu(cs%isdw-1:cs%iedw,cs%jsdw:cs%jedw)) ; cs%IdxCu(:,:) = 0.0
4610 alloc_(cs%IdyCv(cs%isdw:cs%iedw,cs%jsdw-1:cs%jedw)) ; cs%IdyCv(:,:) = 0.0
4611 alloc_(cs%dy_Cu(cs%isdw-1:cs%iedw,cs%jsdw:cs%jedw)) ; cs%dy_Cu(:,:) = 0.0
4612 alloc_(cs%dx_Cv(cs%isdw:cs%iedw,cs%jsdw-1:cs%jedw)) ; cs%dx_Cv(:,:) = 0.0
4613 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
4614 cs%IareaT(i,j) = g%IareaT(i,j)
4615 cs%bathyT(i,j) = g%bathyT(i,j)
4620 do j=g%jsd,g%jed ;
do i=g%IsdB,g%IedB
4621 cs%IdxCu(i,j) = g%IdxCu(i,j) ; cs%dy_Cu(i,j) = g%dy_Cu(i,j)
4623 do j=g%JsdB,g%JedB ;
do i=g%isd,g%ied
4624 cs%IdyCv(i,j) = g%IdyCv(i,j) ; cs%dx_Cv(i,j) = g%dx_Cv(i,j)
4628 call create_group_pass(pass_static_data, cs%IdxCu, cs%IdyCv, cs%BT_domain, to_all+scalar_pair)
4629 call create_group_pass(pass_static_data, cs%dy_Cu, cs%dx_Cv, cs%BT_domain, to_all+scalar_pair)
4630 call do_group_pass(pass_static_data, cs%BT_domain)
4632 if (cs%linearized_BT_PV)
then 4633 alloc_(cs%q_D(cs%isdw-1:cs%iedw,cs%jsdw-1:cs%jedw))
4634 alloc_(cs%D_u_Cor(cs%isdw-1:cs%iedw,cs%jsdw:cs%jedw))
4635 alloc_(cs%D_v_Cor(cs%isdw:cs%iedw,cs%jsdw-1:cs%jedw))
4636 cs%q_D(:,:) = 0.0 ; cs%D_u_Cor(:,:) = 0.0 ; cs%D_v_Cor(:,:) = 0.0
4639 do j=js,je ;
do i=is-1,ie
4640 cs%D_u_Cor(i,j) = 0.5 * (max(mean_sl+g%bathyT(i+1,j),0.0) + max(mean_sl+g%bathyT(i,j),0.0))
4642 do j=js-1,je ;
do i=is,ie
4643 cs%D_v_Cor(i,j) = 0.5 * (max(mean_sl+g%bathyT(i,j+1),0.0) + max(mean_sl+g%bathyT(i,j),0.0))
4645 do j=js-1,je ;
do i=is-1,ie
4646 if (g%mask2dT(i,j)+g%mask2dT(i,j+1)+g%mask2dT(i+1,j)+g%mask2dT(i+1,j+1)>0.)
then 4647 cs%q_D(i,j) = 0.25 * (cs%BT_Coriolis_scale * g%CoriolisBu(i,j)) * &
4648 ((g%areaT(i,j) + g%areaT(i+1,j+1)) + (g%areaT(i+1,j) + g%areaT(i,j+1))) / &
4649 (max(((g%areaT(i,j) * max(mean_sl+g%bathyT(i,j),0.0) + &
4650 g%areaT(i+1,j+1) * max(mean_sl+g%bathyT(i+1,j+1),0.0)) + &
4651 (g%areaT(i+1,j) * max(mean_sl+g%bathyT(i+1,j),0.0) + &
4652 g%areaT(i,j+1) * max(mean_sl+g%bathyT(i,j+1),0.0))), gv%H_to_Z*gv%H_subroundoff) )
4659 call create_group_pass(pass_q_d_cor, cs%q_D, cs%BT_Domain, to_all, position=corner)
4662 call do_group_pass(pass_q_d_cor, cs%BT_Domain)
4665 if (cs%linear_wave_drag)
then 4666 alloc_(cs%lin_drag_u(isdb:iedb,jsd:jed)) ; cs%lin_drag_u(:,:) = 0.0
4667 alloc_(cs%lin_drag_v(isd:ied,jsdb:jedb)) ; cs%lin_drag_v(:,:) = 0.0
4669 if (len_trim(wave_drag_file) > 0)
then 4670 inputdir =
"." ;
call get_param(param_file, mdl,
"INPUTDIR", inputdir)
4671 wave_drag_file = trim(slasher(inputdir))//trim(wave_drag_file)
4672 call log_param(param_file, mdl,
"INPUTDIR/BT_WAVE_DRAG_FILE", wave_drag_file)
4674 allocate(lin_drag_h(isd:ied,jsd:jed)) ; lin_drag_h(:,:) = 0.0
4676 call mom_read_data(wave_drag_file, wave_drag_var, lin_drag_h, g%Domain, scale=us%m_to_Z*us%T_to_s)
4677 call pass_var(lin_drag_h, g%Domain)
4678 do j=js,je ;
do i=is-1,ie
4679 cs%lin_drag_u(i,j) = (gv%Z_to_H * wave_drag_scale) * &
4680 0.5 * (lin_drag_h(i,j) + lin_drag_h(i+1,j))
4682 do j=js-1,je ;
do i=is,ie
4683 cs%lin_drag_v(i,j) = (gv%Z_to_H * wave_drag_scale) * &
4684 0.5 * (lin_drag_h(i,j) + lin_drag_h(i,j+1))
4686 deallocate(lin_drag_h)
4690 cs%dtbt_fraction = 0.98 ;
if (dtbt_input < 0.0) cs%dtbt_fraction = -dtbt_input
4695 if ((us%s_to_T_restart /= 0.0) .and. (us%s_to_T_restart /= us%s_to_T)) &
4696 dtbt_tmp = (us%s_to_T / us%s_to_T_restart) * cs%dtbt
4701 do k=1,g%ke ; gtot_estimate = gtot_estimate + gv%g_prime(k) ;
enddo 4702 call set_dtbt(g, gv, us, cs, gtot_est=gtot_estimate, ssh_add=ssh_extra)
4704 if (dtbt_input > 0.0)
then 4705 cs%dtbt = us%s_to_T * dtbt_input
4706 elseif (dtbt_tmp > 0.0)
then 4709 if ((dtbt_tmp > 0.0) .and. (dtbt_input > 0.0)) calc_dtbt = .false.
4711 call log_param(param_file, mdl,
"DTBT as used", cs%dtbt*us%T_to_s)
4712 call log_param(param_file, mdl,
"estimated maximum DTBT", cs%dtbt_max*us%T_to_s)
4717 if (gv%Boussinesq)
then 4718 thickness_units =
"m" ; flux_units =
"m3 s-1" 4720 thickness_units =
"kg m-2" ; flux_units =
"kg s-1" 4723 cs%id_PFu_bt = register_diag_field(
'ocean_model',
'PFuBT', diag%axesCu1, time, &
4724 'Zonal Anomalous Barotropic Pressure Force Force Acceleration',
'm s-2', conversion=us%L_T2_to_m_s2)
4725 cs%id_PFv_bt = register_diag_field(
'ocean_model',
'PFvBT', diag%axesCv1, time, &
4726 'Meridional Anomalous Barotropic Pressure Force Acceleration',
'm s-2', conversion=us%L_T2_to_m_s2)
4727 cs%id_Coru_bt = register_diag_field(
'ocean_model',
'CoruBT', diag%axesCu1, time, &
4728 'Zonal Barotropic Coriolis Acceleration',
'm s-2', conversion=us%L_T2_to_m_s2)
4729 cs%id_Corv_bt = register_diag_field(
'ocean_model',
'CorvBT', diag%axesCv1, time, &
4730 'Meridional Barotropic Coriolis Acceleration',
'm s-2', conversion=us%L_T2_to_m_s2)
4731 cs%id_uaccel = register_diag_field(
'ocean_model',
'u_accel_bt', diag%axesCu1, time, &
4732 'Barotropic zonal acceleration',
'm s-2', conversion=us%L_T2_to_m_s2)
4733 cs%id_vaccel = register_diag_field(
'ocean_model',
'v_accel_bt', diag%axesCv1, time, &
4734 'Barotropic meridional acceleration',
'm s-2', conversion=us%L_T2_to_m_s2)
4735 cs%id_ubtforce = register_diag_field(
'ocean_model',
'ubtforce', diag%axesCu1, time, &
4736 'Barotropic zonal acceleration from baroclinic terms',
'm s-2', conversion=us%L_T2_to_m_s2)
4737 cs%id_vbtforce = register_diag_field(
'ocean_model',
'vbtforce', diag%axesCv1, time, &
4738 'Barotropic meridional acceleration from baroclinic terms',
'm s-2', conversion=us%L_T2_to_m_s2)
4739 cs%id_ubtdt = register_diag_field(
'ocean_model',
'ubt_dt', diag%axesCu1, time, &
4740 'Barotropic zonal acceleration',
'm s-2', conversion=us%L_T2_to_m_s2)
4741 cs%id_vbtdt = register_diag_field(
'ocean_model',
'vbt_dt', diag%axesCv1, time, &
4742 'Barotropic meridional acceleration',
'm s-2', conversion=us%L_T2_to_m_s2)
4744 cs%id_eta_bt = register_diag_field(
'ocean_model',
'eta_bt', diag%axesT1, time, &
4745 'Barotropic end SSH', thickness_units, conversion=gv%H_to_m)
4746 cs%id_ubt = register_diag_field(
'ocean_model',
'ubt', diag%axesCu1, time, &
4747 'Barotropic end zonal velocity',
'm s-1', conversion=us%L_T_to_m_s)
4748 cs%id_vbt = register_diag_field(
'ocean_model',
'vbt', diag%axesCv1, time, &
4749 'Barotropic end meridional velocity',
'm s-1', conversion=us%L_T_to_m_s)
4750 cs%id_eta_st = register_diag_field(
'ocean_model',
'eta_st', diag%axesT1, time, &
4751 'Barotropic start SSH', thickness_units, conversion=gv%H_to_m)
4752 cs%id_ubt_st = register_diag_field(
'ocean_model',
'ubt_st', diag%axesCu1, time, &
4753 'Barotropic start zonal velocity',
'm s-1', conversion=us%L_T_to_m_s)
4754 cs%id_vbt_st = register_diag_field(
'ocean_model',
'vbt_st', diag%axesCv1, time, &
4755 'Barotropic start meridional velocity',
'm s-1', conversion=us%L_T_to_m_s)
4756 cs%id_ubtav = register_diag_field(
'ocean_model',
'ubtav', diag%axesCu1, time, &
4757 'Barotropic time-average zonal velocity',
'm s-1', conversion=us%L_T_to_m_s)
4758 cs%id_vbtav = register_diag_field(
'ocean_model',
'vbtav', diag%axesCv1, time, &
4759 'Barotropic time-average meridional velocity',
'm s-1', conversion=us%L_T_to_m_s)
4760 cs%id_eta_cor = register_diag_field(
'ocean_model',
'eta_cor', diag%axesT1, time, &
4761 'Corrective mass flux',
'm s-1', conversion=gv%H_to_m)
4762 cs%id_visc_rem_u = register_diag_field(
'ocean_model',
'visc_rem_u', diag%axesCuL, time, &
4763 'Viscous remnant at u',
'nondim')
4764 cs%id_visc_rem_v = register_diag_field(
'ocean_model',
'visc_rem_v', diag%axesCvL, time, &
4765 'Viscous remnant at v',
'nondim')
4766 cs%id_gtotn = register_diag_field(
'ocean_model',
'gtot_n', diag%axesT1, time, &
4767 'gtot to North',
'm s-2', conversion=gv%m_to_H*(us%L_T_to_m_s**2))
4768 cs%id_gtots = register_diag_field(
'ocean_model',
'gtot_s', diag%axesT1, time, &
4769 'gtot to South',
'm s-2', conversion=gv%m_to_H*(us%L_T_to_m_s**2))
4770 cs%id_gtote = register_diag_field(
'ocean_model',
'gtot_e', diag%axesT1, time, &
4771 'gtot to East',
'm s-2', conversion=gv%m_to_H*(us%L_T_to_m_s**2))
4772 cs%id_gtotw = register_diag_field(
'ocean_model',
'gtot_w', diag%axesT1, time, &
4773 'gtot to West',
'm s-2', conversion=gv%m_to_H*(us%L_T_to_m_s**2))
4774 cs%id_eta_hifreq = register_diag_field(
'ocean_model',
'eta_hifreq', diag%axesT1, time, &
4775 'High Frequency Barotropic SSH', thickness_units, conversion=gv%H_to_m)
4776 cs%id_ubt_hifreq = register_diag_field(
'ocean_model',
'ubt_hifreq', diag%axesCu1, time, &
4777 'High Frequency Barotropic zonal velocity',
'm s-1', conversion=us%L_T_to_m_s)
4778 cs%id_vbt_hifreq = register_diag_field(
'ocean_model',
'vbt_hifreq', diag%axesCv1, time, &
4779 'High Frequency Barotropic meridional velocity',
'm s-1', conversion=us%L_T_to_m_s)
4780 cs%id_eta_pred_hifreq = register_diag_field(
'ocean_model',
'eta_pred_hifreq', diag%axesT1, time, &
4781 'High Frequency Predictor Barotropic SSH', thickness_units, &
4782 conversion=gv%H_to_m)
4783 cs%id_uhbt_hifreq = register_diag_field(
'ocean_model',
'uhbt_hifreq', diag%axesCu1, time, &
4784 'High Frequency Barotropic zonal transport',
'm3 s-1', &
4785 conversion=gv%H_to_m*us%L_to_m*us%L_T_to_m_s)
4786 cs%id_vhbt_hifreq = register_diag_field(
'ocean_model',
'vhbt_hifreq', diag%axesCv1, time, &
4787 'High Frequency Barotropic meridional transport',
'm3 s-1', &
4788 conversion=gv%H_to_m*us%L_to_m*us%L_T_to_m_s)
4789 cs%id_frhatu = register_diag_field(
'ocean_model',
'frhatu', diag%axesCuL, time, &
4790 'Fractional thickness of layers in u-columns',
'nondim')
4791 cs%id_frhatv = register_diag_field(
'ocean_model',
'frhatv', diag%axesCvL, time, &
4792 'Fractional thickness of layers in v-columns',
'nondim')
4793 cs%id_frhatu1 = register_diag_field(
'ocean_model',
'frhatu1', diag%axesCuL, time, &
4794 'Predictor Fractional thickness of layers in u-columns',
'nondim')
4795 cs%id_frhatv1 = register_diag_field(
'ocean_model',
'frhatv1', diag%axesCvL, time, &
4796 'Predictor Fractional thickness of layers in v-columns',
'nondim')
4797 cs%id_uhbt = register_diag_field(
'ocean_model',
'uhbt', diag%axesCu1, time, &
4798 'Barotropic zonal transport averaged over a baroclinic step',
'm3 s-1', &
4799 conversion=gv%H_to_m*us%L_to_m*us%L_T_to_m_s)
4800 cs%id_vhbt = register_diag_field(
'ocean_model',
'vhbt', diag%axesCv1, time, &
4801 'Barotropic meridional transport averaged over a baroclinic step',
'm3 s-1', &
4802 conversion=gv%H_to_m*us%L_to_m*us%L_T_to_m_s)
4804 if (use_bt_cont_type)
then 4805 cs%id_BTC_FA_u_EE = register_diag_field(
'ocean_model',
'BTC_FA_u_EE', diag%axesCu1, time, &
4806 'BTCont type far east face area',
'm2', conversion=us%L_to_m*gv%H_to_m)
4807 cs%id_BTC_FA_u_E0 = register_diag_field(
'ocean_model',
'BTC_FA_u_E0', diag%axesCu1, time, &
4808 'BTCont type near east face area',
'm2', conversion=us%L_to_m*gv%H_to_m)
4809 cs%id_BTC_FA_u_WW = register_diag_field(
'ocean_model',
'BTC_FA_u_WW', diag%axesCu1, time, &
4810 'BTCont type far west face area',
'm2', conversion=us%L_to_m*gv%H_to_m)
4811 cs%id_BTC_FA_u_W0 = register_diag_field(
'ocean_model',
'BTC_FA_u_W0', diag%axesCu1, time, &
4812 'BTCont type near west face area',
'm2', conversion=us%L_to_m*gv%H_to_m)
4813 cs%id_BTC_ubt_EE = register_diag_field(
'ocean_model',
'BTC_ubt_EE', diag%axesCu1, time, &
4814 'BTCont type far east velocity',
'm s-1', conversion=us%L_T_to_m_s)
4815 cs%id_BTC_ubt_WW = register_diag_field(
'ocean_model',
'BTC_ubt_WW', diag%axesCu1, time, &
4816 'BTCont type far west velocity',
'm s-1', conversion=us%L_T_to_m_s)
4820 cs%id_BTC_FA_v_NN = register_diag_field(
'ocean_model',
'BTC_FA_v_NN', diag%axesCv1, time, &
4821 'BTCont type far north face area',
'm2', conversion=us%L_to_m*gv%H_to_m)
4822 cs%id_BTC_FA_v_N0 = register_diag_field(
'ocean_model',
'BTC_FA_v_N0', diag%axesCv1, time, &
4823 'BTCont type near north face area',
'm2', conversion=us%L_to_m*gv%H_to_m)
4824 cs%id_BTC_FA_v_SS = register_diag_field(
'ocean_model',
'BTC_FA_v_SS', diag%axesCv1, time, &
4825 'BTCont type far south face area',
'm2', conversion=us%L_to_m*gv%H_to_m)
4826 cs%id_BTC_FA_v_S0 = register_diag_field(
'ocean_model',
'BTC_FA_v_S0', diag%axesCv1, time, &
4827 'BTCont type near south face area',
'm2', conversion=us%L_to_m*gv%H_to_m)
4828 cs%id_BTC_vbt_NN = register_diag_field(
'ocean_model',
'BTC_vbt_NN', diag%axesCv1, time, &
4829 'BTCont type far north velocity',
'm s-1', conversion=us%L_T_to_m_s)
4830 cs%id_BTC_vbt_SS = register_diag_field(
'ocean_model',
'BTC_vbt_SS', diag%axesCv1, time, &
4831 'BTCont type far south velocity',
'm s-1', conversion=us%L_T_to_m_s)
4838 cs%id_uhbt0 = register_diag_field(
'ocean_model',
'uhbt0', diag%axesCu1, time, &
4839 'Barotropic zonal transport difference',
'm3 s-1', conversion=gv%H_to_m*us%L_to_m**2*us%s_to_T)
4840 cs%id_vhbt0 = register_diag_field(
'ocean_model',
'vhbt0', diag%axesCv1, time, &
4841 'Barotropic meridional transport difference',
'm3 s-1', conversion=gv%H_to_m*us%L_to_m**2*us%s_to_T)
4843 if (cs%id_frhatu1 > 0)
call safe_alloc_ptr(cs%frhatu1, isdb,iedb,jsd,jed,nz)
4844 if (cs%id_frhatv1 > 0)
call safe_alloc_ptr(cs%frhatv1, isd,ied,jsdb,jedb,nz)
4848 call btcalc(h, g, gv, cs, may_use_default=.true.)
4849 cs%ubtav(:,:) = 0.0 ; cs%vbtav(:,:) = 0.0
4850 do k=1,nz ;
do j=js,je ;
do i=is-1,ie
4851 cs%ubtav(i,j) = cs%ubtav(i,j) + cs%frhatu(i,j,k) * u(i,j,k)
4852 enddo ;
enddo ;
enddo 4853 do k=1,nz ;
do j=js-1,je ;
do i=is,ie
4854 cs%vbtav(i,j) = cs%vbtav(i,j) + cs%frhatv(i,j,k) * v(i,j,k)
4855 enddo ;
enddo ;
enddo 4856 elseif ((us%s_to_T_restart*us%m_to_L_restart /= 0.0) .and. &
4857 (us%m_to_L*us%s_to_T_restart) /= (us%m_to_L_restart*us%s_to_T))
then 4858 vel_rescale = (us%m_to_L*us%s_to_T_restart) / (us%m_to_L_restart*us%s_to_T)
4859 do j=js,je ;
do i=is-1,ie ; cs%ubtav(i,j) = vel_rescale * cs%ubtav(i,j) ;
enddo ;
enddo 4860 do j=js-1,je ;
do i=is,ie ; cs%vbtav(i,j) = vel_rescale * cs%vbtav(i,j) ;
enddo ;
enddo 4863 if (cs%gradual_BT_ICs)
then 4866 do j=js,je ;
do i=is-1,ie ; cs%ubt_IC(i,j) = cs%ubtav(i,j) ;
enddo ;
enddo 4867 do j=js-1,je ;
do i=is,ie ; cs%vbt_IC(i,j) = cs%vbtav(i,j) ;
enddo ;
enddo 4868 elseif ((us%s_to_T_restart*us%m_to_L_restart /= 0.0) .and. &
4869 (us%m_to_L*us%s_to_T_restart) /= (us%m_to_L_restart*us%s_to_T))
then 4870 vel_rescale = (us%m_to_L*us%s_to_T_restart) / (us%m_to_L_restart*us%s_to_T)
4871 do j=js,je ;
do i=is-1,ie ; cs%ubt_IC(i,j) = vel_rescale * cs%ubt_IC(i,j) ;
enddo ;
enddo 4872 do j=js-1,je ;
do i=is,ie ; cs%vbt_IC(i,j) = vel_rescale * cs%vbt_IC(i,j) ;
enddo ;
enddo 4877 if (.not.cs%nonlin_stress)
then 4879 do j=js,je ;
do i=is-1,ie
4880 if (g%mask2dCu(i,j)>0.)
then 4881 cs%IDatu(i,j) = g%mask2dCu(i,j) * 2.0 / ((g%bathyT(i+1,j) + g%bathyT(i,j)) + 2.0*mean_sl)
4886 do j=js-1,je ;
do i=is,ie
4887 if (g%mask2dCv(i,j)>0.)
then 4888 cs%IDatv(i,j) = g%mask2dCv(i,j) * 2.0 / ((g%bathyT(i,j+1) + g%bathyT(i,j)) + 2.0*mean_sl)
4895 call find_face_areas(datu, datv, g, gv, us, cs, ms, halo=1)
4896 if ((cs%bound_BT_corr) .and. .not.(use_bt_cont_type .and. cs%BT_cont_bounds))
then 4899 do j=js,je ;
do i=is,ie
4900 cs%eta_cor_bound(i,j) = g%IareaT(i,j) * 0.1 * cs%maxvel * &
4901 ((datu(i-1,j) + datu(i,j)) + (datv(i,j) + datv(i,j-1)))
4905 if (cs%gradual_BT_ICs) &
4908 call do_group_pass(pass_bt_hbt_btav, g%Domain)
4911 id_clock_calc_pre = cpu_clock_id(
'(Ocean BT pre-calcs only)', grain=clock_routine)
4912 id_clock_pass_pre = cpu_clock_id(
'(Ocean BT pre-step halo updates)', grain=clock_routine)
4913 id_clock_calc = cpu_clock_id(
'(Ocean BT stepping calcs only)', grain=clock_routine)
4914 id_clock_pass_step = cpu_clock_id(
'(Ocean BT stepping halo updates)', grain=clock_routine)
4915 id_clock_calc_post = cpu_clock_id(
'(Ocean BT post-calcs only)', grain=clock_routine)
4916 id_clock_pass_post = cpu_clock_id(
'(Ocean BT post-step halo updates)', grain=clock_routine)
4917 if (dtbt_input <= 0.0) &
4918 id_clock_sync = cpu_clock_id(
'(Ocean BT global synch)', grain=clock_routine)
4920 end subroutine barotropic_init
4923 subroutine barotropic_get_tav(CS, ubtav, vbtav, G, US)
4926 real,
dimension(SZIB_(G),SZJ_(G)),
intent(inout) :: ubtav
4928 real,
dimension(SZI_(G),SZJB_(G)),
intent(inout) :: vbtav
4934 do j=g%jsc,g%jec ;
do i=g%isc-1,g%iec
4935 ubtav(i,j) = cs%ubtav(i,j)
4938 do j=g%jsc-1,g%jec ;
do i=g%isc,g%iec
4939 vbtav(i,j) = cs%vbtav(i,j)
4942 end subroutine barotropic_get_tav
4946 subroutine barotropic_end(CS)
4948 dealloc_(cs%frhatu) ; dealloc_(cs%frhatv)
4949 dealloc_(cs%IDatu) ; dealloc_(cs%IDatv)
4950 dealloc_(cs%ubtav) ; dealloc_(cs%vbtav)
4951 dealloc_(cs%eta_cor)
4952 dealloc_(cs%ua_polarity) ; dealloc_(cs%va_polarity)
4953 if (cs%bound_BT_corr)
then 4954 dealloc_(cs%eta_cor_bound)
4957 call destroy_bt_obc(cs%BT_OBC)
4960 end subroutine barotropic_end
4964 subroutine register_barotropic_restarts(HI, GV, param_file, CS, restart_CS)
4974 character(len=40) :: mdl =
"MOM_barotropic" 4975 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
4977 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed
4978 isdb = hi%IsdB ; iedb = hi%IedB ; jsdb = hi%JsdB ; jedb = hi%JedB
4980 if (
associated(cs))
then 4981 call mom_error(warning,
"register_barotropic_restarts called with an associated "// &
4982 "control structure.")
4987 call get_param(param_file, mdl,
"GRADUAL_BT_ICS", cs%gradual_BT_ICs, &
4988 "If true, adjust the initial conditions for the "//&
4989 "barotropic solver to the values from the layered "//&
4990 "solution over a whole timestep instead of instantly. "//&
4991 "This is a decent approximation to the inclusion of "//&
4992 "sum(u dh_dt) while also correcting for truncation errors.", &
4993 default=.false., do_not_log=.true.)
4995 alloc_(cs%ubtav(isdb:iedb,jsd:jed)) ; cs%ubtav(:,:) = 0.0
4996 alloc_(cs%vbtav(isd:ied,jsdb:jedb)) ; cs%vbtav(:,:) = 0.0
4997 if (cs%gradual_BT_ICs)
then 4998 alloc_(cs%ubt_IC(isdb:iedb,jsd:jed)) ; cs%ubt_IC(:,:) = 0.0
4999 alloc_(cs%vbt_IC(isd:ied,jsdb:jedb)) ; cs%vbt_IC(:,:) = 0.0
5002 vd(2) = var_desc(
"ubtav",
"m s-1",
"Time mean barotropic zonal velocity", &
5003 hor_grid=
'u', z_grid=
'1')
5004 vd(3) = var_desc(
"vbtav",
"m s-1",
"Time mean barotropic meridional velocity",&
5005 hor_grid=
'v', z_grid=
'1')
5008 if (cs%gradual_BT_ICs)
then 5009 vd(2) = var_desc(
"ubt_IC",
"m s-1", &
5010 longname=
"Next initial condition for the barotropic zonal velocity", &
5011 hor_grid=
'u', z_grid=
'1')
5012 vd(3) = var_desc(
"vbt_IC",
"m s-1", &
5013 longname=
"Next initial condition for the barotropic meridional velocity",&
5014 hor_grid=
'v', z_grid=
'1')
5020 longname=
"Barotropic timestep", units=
"seconds")
5022 end subroutine register_barotropic_restarts
Wraps the FMS time manager functions.
This module implements boundary forcing for MOM6.
Ocean grid type. See mom_grid for details.
A structure that can be parsed to read and document run-time parameters.
Provides the ocean grid type.
Open boundary segment data structure.
Wraps the MPP cpu clock functions.
Register fields for restarts.
This module contains I/O framework code.
The MOM6 facility to parse input files for runtime parameters.
Defines the horizontal index type (hor_index_type) used for providing index ranges.
Pointers to arrays with accelerations, which can later be used for derived diagnostics,...
Register a pair of restart fieilds whose rotations map onto each other.
An overloaded interface to log the values of various types of parameters.
A desciption of the functional dependence of transport at a u-point.
Do a halo update on a pair of arrays representing the two components of a vector.
Container for horizontal index ranges for data, computational and global domains.
Structure that contains pointers to the mechanical forcing at the surface used to drive the liquid oc...
A desciption of the functional dependence of transport at a v-point.
Make a diagnostic available for averaging or output.
A restart registry and the control structure for restarts.
Container for information about the summed layer transports and how they will vary as the barotropic ...
Describes various unit conversion factors.
Copy one MOM_domain_type into another.
A container for passing around active tracer point memory limits.
Tidal contributions to geopotential.
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Describes the decomposed MOM domain and has routines for communications across PEs.
Routines for error handling and I/O management.
The MOM6 facility for reading and writing restart files, and querying what has been read.
Type for describing a variable, typically a tracer.
An overloaded interface to log version information about modules.
Describes the vertical ocean grid, including unit conversion factors.
The MOM_domain_type contains information about the domain decompositoin.
Set up a group of halo updates.
The control structure for the MOM_tidal_forcing module.
Indicate whether a field has been read from a restart file.
The barotropic stepping control stucture.
Controls where open boundary conditions are applied.
Provides a transparent vertical ocean grid type and supporting routines.
Provides transparent structures with groups of MOM6 variables and supporting routines.
Do a halo update on an array.
Read a data field from a file.
The barotropic stepping open boundary condition type.
An overloaded interface to read and log the values of various types of parameters.
Provides checksumming functions for debugging.