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