7 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine
30 implicit none ;
private
32 #include <MOM_memory.h>
37 # define WHALOI_ MAX(BTHALO_-NIHALO_,0)
38 # define WHALOJ_ MAX(BTHALO_-NJHALO_,0)
39 # define NIMEMW_ 1-WHALOI_:NIMEM_+WHALOI_
40 # define NJMEMW_ 1-WHALOJ_:NJMEM_+WHALOJ_
41 # define NIMEMBW_ -WHALOI_:NIMEM_+WHALOI_
42 # define NJMEMBW_ -WHALOJ_:NJMEM_+WHALOJ_
43 # define SZIW_(G) NIMEMW_
44 # define SZJW_(G) NJMEMW_
45 # define SZIBW_(G) NIMEMBW_
46 # define SZJBW_(G) NJMEMBW_
52 # define SZIW_(G) G%isdw:G%iedw
53 # define SZJW_(G) G%jsdw:G%jedw
54 # define SZIBW_(G) G%isdw-1:G%iedw
55 # define SZJBW_(G) G%jsdw-1:G%jedw
68 real,
dimension(:,:),
pointer :: cg_u => null()
69 real,
dimension(:,:),
pointer :: cg_v => null()
70 real,
dimension(:,:),
pointer :: h_u => null()
71 real,
dimension(:,:),
pointer :: h_v => null()
72 real,
dimension(:,:),
pointer :: uhbt => null()
74 real,
dimension(:,:),
pointer :: vhbt => null()
76 real,
dimension(:,:),
pointer :: ubt_outer => null()
78 real,
dimension(:,:),
pointer :: vbt_outer => null()
80 real,
dimension(:,:),
pointer :: eta_outer_u => null()
82 real,
dimension(:,:),
pointer :: eta_outer_v => null()
84 logical :: apply_u_obcs
85 logical :: apply_v_obcs
87 integer :: is_u_obc, ie_u_obc, js_u_obc, je_u_obc
88 integer :: is_v_obc, ie_v_obc, js_v_obc, je_v_obc
90 logical :: is_alloced = .false.
92 type(group_pass_type) :: pass_uv
93 type(group_pass_type) :: pass_uhvh
94 type(group_pass_type) :: pass_h
95 type(group_pass_type) :: pass_cg
96 type(group_pass_type) :: pass_eta_outer
101 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: frhatu
103 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: frhatv
105 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_) :: idatu
107 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_) :: lin_drag_u
110 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_) :: uhbt_ic
113 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_) :: ubt_ic
116 real allocable_,
dimension(NIMEMB_PTR_,NJMEM_) :: ubtav
118 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_) :: idatv
120 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_) :: lin_drag_v
123 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_) :: vhbt_ic
126 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_) :: vbt_ic
129 real allocable_,
dimension(NIMEM_,NJMEMB_PTR_) :: vbtav
131 real allocable_,
dimension(NIMEM_,NJMEM_) :: eta_cor
135 real allocable_,
dimension(NIMEM_,NJMEM_) :: eta_cor_bound
138 real allocable_,
dimension(NIMEMW_,NJMEMW_) :: &
139 ua_polarity, & !< Test vector components for checking grid polarity.
140 va_polarity, & !< Test vector components for checking grid polarity.
142 real allocable_,
dimension(NIMEMW_,NJMEMW_) :: iareat
145 real allocable_,
dimension(NIMEMBW_,NJMEMW_) :: &
146 d_u_cor, & !< A simply averaged depth at u points [Z ~> m].
147 dy_cu, & !< A copy of G%dy_Cu with wide halos [L ~> m].
149 real allocable_,
dimension(NIMEMW_,NJMEMBW_) :: &
150 d_v_cor, & !< A simply averaged depth at v points [Z ~> m].
151 dx_cv, & !< A copy of G%dx_Cv with wide halos [L ~> m].
153 real allocable_,
dimension(NIMEMBW_,NJMEMBW_) :: &
156 real,
dimension(:,:,:),
pointer :: frhatu1 => null()
157 real,
dimension(:,:,:),
pointer :: frhatv1 => null()
163 real :: dtbt_fraction
170 integer :: nstep_last = 0
178 logical :: bound_bt_corr
183 logical :: gradual_bt_ics
194 logical :: nonlinear_continuity
196 integer :: nonlin_cont_update_period
200 logical :: bt_project_velocity
208 logical :: dynamic_psurf
210 real :: dmin_dyn_psurf
212 real :: ice_strength_length
215 real :: const_dyn_psurf
220 integer :: hvel_scheme
224 logical :: strong_drag
226 logical :: linear_wave_drag
229 logical :: linearized_bt_pv
232 logical :: use_wide_halos
234 logical :: clip_velocity
240 real :: vel_underflow
247 real :: maxcfl_bt_cont
252 logical :: bt_cont_bounds
255 logical :: visc_rem_u_uh0
259 logical :: adjust_bt_cont
262 logical :: use_old_coriolis_bracket_bug
265 type(time_type),
pointer :: time => null()
271 logical :: module_is_initialized = .false.
278 type(group_pass_type) :: pass_q_dcor
279 type(group_pass_type) :: pass_gtot
280 type(group_pass_type) :: pass_tmp_uv
281 type(group_pass_type) :: pass_eta_bt_rem
282 type(group_pass_type) :: pass_force_hbt0_cor_ref
283 type(group_pass_type) :: pass_dat_uv
284 type(group_pass_type) :: pass_eta_ubt
285 type(group_pass_type) :: pass_etaav
286 type(group_pass_type) :: pass_ubt_cor
287 type(group_pass_type) :: pass_ubta_uhbta
288 type(group_pass_type) :: pass_e_anom
291 integer :: id_pfu_bt = -1, id_pfv_bt = -1, id_coru_bt = -1, id_corv_bt = -1
292 integer :: id_ubtforce = -1, id_vbtforce = -1, id_uaccel = -1, id_vaccel = -1
293 integer :: id_visc_rem_u = -1, id_visc_rem_v = -1, id_eta_cor = -1
294 integer :: id_ubt = -1, id_vbt = -1, id_eta_bt = -1, id_ubtav = -1, id_vbtav = -1
295 integer :: id_ubt_st = -1, id_vbt_st = -1, id_eta_st = -1
296 integer :: id_ubt_hifreq = -1, id_vbt_hifreq = -1, id_eta_hifreq = -1
297 integer :: id_uhbt_hifreq = -1, id_vhbt_hifreq = -1, id_eta_pred_hifreq = -1
298 integer :: id_gtotn = -1, id_gtots = -1, id_gtote = -1, id_gtotw = -1
299 integer :: id_uhbt = -1, id_frhatu = -1, id_vhbt = -1, id_frhatv = -1
300 integer :: id_frhatu1 = -1, id_frhatv1 = -1
302 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
303 integer :: id_btc_ubt_ee = -1, id_btc_ubt_ww = -1
304 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
305 integer :: id_btc_vbt_nn = -1, id_btc_vbt_ss = -1
306 integer :: id_uhbt0 = -1, id_vhbt0 = -1
353 integer :: isdw, iedw, jsdw, jedw
383 subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, &
384 eta_PF_in, U_Cor, V_Cor, accel_layer_u, accel_layer_v, &
385 eta_out, uhbtav, vhbtav, G, GV, US, CS, &
386 visc_rem_u, visc_rem_v, etaav, OBC, BT_cont, eta_PF_start, &
387 taux_bot, tauy_bot, uh0, vh0, u_uh0, v_vh0)
391 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: u_in
393 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: v_in
395 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: eta_in
397 real,
intent(in) :: dt
398 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: bc_accel_u
400 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: bc_accel_v
403 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: pbce
406 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: eta_pf_in
412 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: u_cor
414 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: v_cor
416 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(out) :: accel_layer_u
418 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(out) :: accel_layer_v
420 real,
dimension(SZI_(G),SZJ_(G)),
intent(out) :: eta_out
422 real,
dimension(SZIB_(G),SZJ_(G)),
intent(out) :: uhbtav
425 real,
dimension(SZI_(G),SZJB_(G)),
intent(out) :: vhbtav
430 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: visc_rem_u
436 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: visc_rem_v
437 real,
dimension(SZI_(G),SZJ_(G)),
optional,
intent(out) :: etaav
443 real,
dimension(:,:),
optional,
pointer :: eta_pf_start
446 real,
dimension(:,:),
optional,
pointer :: taux_bot
448 real,
dimension(:,:),
optional,
pointer :: tauy_bot
450 real,
dimension(:,:,:),
optional,
pointer :: uh0
452 real,
dimension(:,:,:),
optional,
pointer :: u_uh0
454 real,
dimension(:,:,:),
optional,
pointer :: vh0
456 real,
dimension(:,:,:),
optional,
pointer :: v_vh0
460 real :: ubt_cor(szib_(g),szj_(g))
461 real :: vbt_cor(szi_(g),szjb_(g))
463 real :: wt_u(szib_(g),szj_(g),szk_(g))
464 real :: wt_v(szi_(g),szjb_(g),szk_(g))
467 real,
dimension(SZIB_(G),SZJ_(G)) :: &
470 real,
dimension(SZI_(G),SZJB_(G)) :: &
473 real,
dimension(SZI_(G),SZJ_(G)) :: &
479 real :: q(szibw_(cs),szjbw_(cs))
480 real,
dimension(SZIBW_(CS),SZJW_(CS)) :: &
513 real,
dimension(SZIW_(CS),SZJBW_(CS)) :: &
544 real,
target,
dimension(SZIW_(CS),SZJW_(CS)) :: &
548 real,
dimension(:,:),
pointer :: &
551 real,
dimension(SZIW_(CS),SZJW_(CS)) :: &
578 real,
dimension(SZIBW_(CS),SZJW_(CS)) :: &
579 ubt_prev, uhbt_prev, ubt_sum_prev, uhbt_sum_prev, ubt_wtd_prev
580 real,
dimension(SZIW_(CS),SZJBW_(CS)) :: &
581 vbt_prev, vhbt_prev, vbt_sum_prev, vhbt_sum_prev, vbt_wtd_prev
584 real :: mass_accel_to_z
609 logical :: do_hifreq_output
610 logical :: use_bt_cont, do_ave, find_etaav, find_pf, find_cor
611 logical :: ice_is_rigid, nonblock_setup, interp_eta_pf
612 logical :: project_velocity, add_uh0
616 real :: ice_strength = 0.0
623 real :: u_max_cor, v_max_cor
626 real :: accel_underflow
628 real,
allocatable,
dimension(:) :: wt_vel, wt_eta, wt_accel, wt_trans, wt_accel2
629 real :: sum_wt_vel, sum_wt_eta, sum_wt_accel, sum_wt_trans
630 real :: i_sum_wt_vel, i_sum_wt_eta, i_sum_wt_accel, i_sum_wt_trans
632 real :: trans_wt1, trans_wt2
635 logical :: apply_obcs, apply_obc_flather, apply_obc_open
637 character(len=200) :: mesg
638 integer :: isv, iev, jsv, jev
640 integer :: isvf, ievf, jsvf, jevf, num_cycles
641 integer :: i, j, k, n
642 integer :: is, ie, js, je, nz, isq, ieq, jsq, jeq
643 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
644 integer :: ioff, joff
646 if (.not.
associated(cs))
call mom_error(fatal, &
647 "btstep: Module MOM_barotropic must be initialized before it is used.")
648 if (.not.cs%split)
return
649 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
650 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
651 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
652 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
653 ms%isdw = cs%isdw ; ms%iedw = cs%iedw ; ms%jsdw = cs%jsdw ; ms%jedw = cs%jedw
656 accel_underflow = cs%vel_underflow * idt
658 use_bt_cont = .false.
659 if (
present(bt_cont)) use_bt_cont = (
associated(bt_cont))
661 interp_eta_pf = .false.
662 if (
present(eta_pf_start)) interp_eta_pf = (
associated(eta_pf_start))
664 project_velocity = cs%BT_project_velocity
668 if ((.not.use_bt_cont) .and. cs%Nonlinear_continuity .and. &
669 (cs%Nonlin_cont_update_period > 0)) stencil = 2
671 do_ave = query_averaging_enabled(cs%diag)
672 find_etaav =
present(etaav)
673 find_pf = (do_ave .and. ((cs%id_PFu_bt > 0) .or. (cs%id_PFv_bt > 0)))
674 find_cor = (do_ave .and. ((cs%id_Coru_bt > 0) .or. (cs%id_Corv_bt > 0)))
677 if (
present(uh0)) add_uh0 =
associated(uh0)
678 if (add_uh0 .and. .not.(
present(vh0) .and.
present(u_uh0) .and. &
680 "btstep: vh0, u_uh0, and v_vh0 must be present if uh0 is used.")
681 if (add_uh0 .and. .not.(
associated(vh0) .and.
associated(u_uh0) .and. &
682 associated(v_vh0)))
call mom_error(fatal, &
683 "btstep: vh0, u_uh0, and v_vh0 must be associated if uh0 is used.")
686 nonblock_setup = g%nonblocking_updates
690 apply_obcs = .false. ; cs%BT_OBC%apply_u_OBCs = .false. ; cs%BT_OBC%apply_v_OBCs = .false.
691 apply_obc_open = .false.
692 apply_obc_flather = .false.
693 if (
present(obc))
then ;
if (
associated(obc))
then
694 cs%BT_OBC%apply_u_OBCs = obc%open_u_BCs_exist_globally .or. obc%specified_u_BCs_exist_globally
695 cs%BT_OBC%apply_v_OBCs = obc%open_v_BCs_exist_globally .or. obc%specified_v_BCs_exist_globally
696 apply_obc_flather = open_boundary_query(obc, apply_flather_obc=.true.)
697 apply_obc_open = open_boundary_query(obc, apply_open_obc=.true.)
698 apply_obcs = open_boundary_query(obc, apply_specified_obc=.true.) .or. &
699 apply_obc_flather .or. apply_obc_open
701 if (apply_obc_flather .and. .not.gv%Boussinesq)
call mom_error(fatal, &
702 "btstep: Flather open boundary conditions have not yet been "// &
703 "implemented for a non-Boussinesq model.")
707 if (cs%use_wide_halos) &
708 num_cycles = min((is-cs%isdw) / stencil, (js-cs%jsdw) / stencil)
709 isvf = is - (num_cycles-1)*stencil ; ievf = ie + (num_cycles-1)*stencil
710 jsvf = js - (num_cycles-1)*stencil ; jevf = je + (num_cycles-1)*stencil
712 nstep = ceiling(dt/cs%dtbt - 0.0001)
713 if (
is_root_pe() .and. (nstep /= cs%nstep_last))
then
714 write(mesg,
'("btstep is using a dynamic barotropic timestep of ", ES12.6, &
715 & " seconds, max ", ES12.6, ".")') (us%T_to_s*dt/nstep), us%T_to_s*cs%dtbt_max
718 cs%nstep_last = nstep
721 instep = 1.0 / real(nstep)
725 mass_accel_to_z = 1.0 / gv%Rho0
726 mass_to_z = us%m_to_Z / gv%Rho0
729 if (project_velocity)
then
730 trans_wt1 = (1.0 + be_proj); trans_wt2 = -be_proj
732 trans_wt1 = bebt ; trans_wt2 = (1.0-bebt)
735 do_hifreq_output = .false.
736 if ((cs%id_ubt_hifreq > 0) .or. (cs%id_vbt_hifreq > 0) .or. &
737 (cs%id_eta_hifreq > 0) .or. (cs%id_eta_pred_hifreq > 0) .or. &
738 (cs%id_uhbt_hifreq > 0) .or. (cs%id_vhbt_hifreq > 0))
then
739 do_hifreq_output = query_averaging_enabled(cs%diag, time_int_in, time_end_in)
740 if (do_hifreq_output) &
741 time_bt_start = time_end_in -
real_to_time(us%T_to_s*dt)
746 if (.not. cs%linearized_BT_PV)
then
751 if ((isq > is-1) .or. (jsq > js-1)) &
754 to_all+scalar_pair, agrid)
756 to_all+scalar_pair, agrid)
758 if (cs%dynamic_psurf) &
760 if (interp_eta_pf)
then
771 cs%BT_Domain, to_all+scalar_pair)
772 if (cs%linear_wave_drag) &
774 cs%BT_Domain, to_all+scalar_pair)
777 if (((g%isd > cs%isdw) .or. (g%jsd > cs%jsdw)) .or. (isq <= is-1) .or. (jsq <= js-1)) &
778 call create_group_pass(cs%pass_force_hbt0_Cor_ref, bt_force_u, bt_force_v, cs%BT_Domain)
779 if (add_uh0)
call create_group_pass(cs%pass_force_hbt0_Cor_ref, uhbt0, vhbt0, cs%BT_Domain)
780 call create_group_pass(cs%pass_force_hbt0_Cor_ref, cor_ref_u, cor_ref_v, cs%BT_Domain)
781 if (.not. use_bt_cont)
then
782 call create_group_pass(cs%pass_Dat_uv, datu, datv, cs%BT_Domain, to_all+scalar_pair)
803 if (cs%linearized_BT_PV)
then
805 do j=jsvf-2,jevf+1 ;
do i=isvf-2,ievf+1
809 do j=jsvf-1,jevf+1 ;
do i=isvf-2,ievf+1
810 dcor_u(i,j) = cs%D_u_Cor(i,j)
813 do j=jsvf-2,jevf+1 ;
do i=isvf-1,ievf+1
814 dcor_v(i,j) = cs%D_v_Cor(i,j)
817 q(:,:) = 0.0 ; dcor_u(:,:) = 0.0 ; dcor_v(:,:) = 0.0
821 do j=js,je ;
do i=is-1,ie
822 dcor_u(i,j) = 0.5 * (g%bathyT(i+1,j) + g%bathyT(i,j))
825 do j=js-1,je ;
do i=is,ie
826 dcor_v(i,j) = 0.5 * (g%bathyT(i,j+1) + g%bathyT(i,j))
829 do j=js-1,je ;
do i=is-1,ie
830 q(i,j) = 0.25 * g%CoriolisBu(i,j) * &
831 ((g%areaT(i,j) + g%areaT(i+1,j+1)) + (g%areaT(i+1,j) + g%areaT(i,j+1))) / &
832 ((g%areaT(i,j) * g%bathyT(i,j) + g%areaT(i+1,j+1) * g%bathyT(i+1,j+1)) + &
833 (g%areaT(i+1,j) * g%bathyT(i+1,j) + g%areaT(i,j+1) * g%bathyT(i,j+1)) )
841 if (nonblock_setup)
then
851 do j=cs%jsdw,cs%jedw ;
do i=cs%isdw,cs%iedw
852 gtot_e(i,j) = 0.0 ; gtot_w(i,j) = 0.0
853 gtot_n(i,j) = 0.0 ; gtot_s(i,j) = 0.0
856 if (interp_eta_pf)
then
857 eta_pf_1(i,j) = 0.0 ; d_eta_pf(i,j) = 0.0
859 p_surf_dyn(i,j) = 0.0
860 if (cs%dynamic_psurf) dyn_coef_eta(i,j) = 0.0
866 do j=cs%jsdw,cs%jedw ;
do i=cs%isdw-1,cs%iedw
867 cor_ref_u(i,j) = 0.0 ; bt_force_u(i,j) = 0.0 ; ubt(i,j) = 0.0
868 datu(i,j) = 0.0 ; bt_rem_u(i,j) = 0.0 ; uhbt0(i,j) = 0.0
871 do j=cs%jsdw-1,cs%jedw ;
do i=cs%isdw,cs%iedw
872 cor_ref_v(i,j) = 0.0 ; bt_force_v(i,j) = 0.0 ; vbt(i,j) = 0.0
873 datv(i,j) = 0.0 ; bt_rem_v(i,j) = 0.0 ; vhbt0(i,j) = 0.0
877 if (interp_eta_pf)
then
879 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
880 eta(i,j) = eta_in(i,j)
881 eta_pf_1(i,j) = eta_pf_start(i,j)
882 d_eta_pf(i,j) = eta_pf_in(i,j) - eta_pf_start(i,j)
886 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
887 eta(i,j) = eta_in(i,j)
888 eta_pf(i,j) = eta_pf_in(i,j)
893 do k=1,nz ;
do j=js,je ;
do i=is-1,ie
896 if (visc_rem_u(i,j,k) <= 0.0)
then ; visc_rem = 0.0
897 elseif (visc_rem_u(i,j,k) >= 1.0)
then ; visc_rem = 1.0
898 elseif (visc_rem_u(i,j,k)**2 > visc_rem_u(i,j,k) - 0.5*instep)
then
899 visc_rem = visc_rem_u(i,j,k)
900 else ; visc_rem = 1.0 - 0.5*instep/visc_rem_u(i,j,k) ;
endif
901 wt_u(i,j,k) = cs%frhatu(i,j,k) * visc_rem
902 enddo ;
enddo ;
enddo
904 do k=1,nz ;
do j=js-1,je ;
do i=is,ie
906 if (visc_rem_v(i,j,k) <= 0.0)
then ; visc_rem = 0.0
907 elseif (visc_rem_v(i,j,k) >= 1.0)
then ; visc_rem = 1.0
908 elseif (visc_rem_v(i,j,k)**2 > visc_rem_v(i,j,k) - 0.5*instep)
then
909 visc_rem = visc_rem_v(i,j,k)
910 else ; visc_rem = 1.0 - 0.5*instep/visc_rem_v(i,j,k) ;
endif
911 wt_v(i,j,k) = cs%frhatv(i,j,k) * visc_rem
912 enddo ;
enddo ;
enddo
917 do j=js-1,je+1 ;
do i=is-1,ie ; ubt_cor(i,j) = 0.0 ;
enddo ;
enddo
919 do j=js-1,je ;
do i=is-1,ie+1 ; vbt_cor(i,j) = 0.0 ;
enddo ;
enddo
921 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
922 ubt_cor(i,j) = ubt_cor(i,j) + wt_u(i,j,k) * u_cor(i,j,k)
923 enddo ;
enddo ;
enddo
925 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
926 vbt_cor(i,j) = vbt_cor(i,j) + wt_v(i,j,k) * v_cor(i,j,k)
927 enddo ;
enddo ;
enddo
935 do k=1,nz ;
do i=is-1,ie
936 gtot_e(i,j) = gtot_e(i,j) + pbce(i,j,k) * wt_u(i,j,k)
937 gtot_w(i+1,j) = gtot_w(i+1,j) + pbce(i+1,j,k) * wt_u(i,j,k)
942 do k=1,nz ;
do i=is,ie
943 gtot_n(i,j) = gtot_n(i,j) + pbce(i,j,k) * wt_v(i,j,k)
944 gtot_s(i,j+1) = gtot_s(i,j+1) + pbce(i,j+1,k) * wt_v(i,j,k)
950 dgeo_de = 1.0 + det_de + cs%G_extra
952 dgeo_de = 1.0 + cs%G_extra
955 if (nonblock_setup .and. .not.cs%linearized_BT_PV)
then
963 if (use_bt_cont)
then
966 if (cs%Nonlinear_continuity)
then
975 call set_up_bt_obc(obc, eta, cs%BT_OBC, cs%BT_Domain, g, gv, us, ms, ievf-ie, use_bt_cont, &
976 datu, datv, btcl_u, btcl_v)
986 do j=js,je ;
do i=is-1,ie
990 bt_force_u(i,j) = forces%taux(i,j) * mass_accel_to_z * cs%IDatu(i,j)*visc_rem_u(i,j,1)
993 do j=js-1,je ;
do i=is,ie
997 bt_force_v(i,j) = forces%tauy(i,j) * mass_accel_to_z * cs%IDatv(i,j)*visc_rem_v(i,j,1)
999 if (
present(taux_bot) .and.
present(tauy_bot))
then
1000 if (
associated(taux_bot) .and.
associated(tauy_bot))
then
1002 do j=js,je ;
do i=is-1,ie
1003 bt_force_u(i,j) = bt_force_u(i,j) - taux_bot(i,j) * mass_to_z * cs%IDatu(i,j)
1006 do j=js-1,je ;
do i=is,ie
1007 bt_force_v(i,j) = bt_force_v(i,j) - tauy_bot(i,j) * mass_to_z * cs%IDatv(i,j)
1015 do j=js,je ;
do k=1,nz ;
do i=isq,ieq
1016 bt_force_u(i,j) = bt_force_u(i,j) + wt_u(i,j,k) * bc_accel_u(i,j,k)
1017 enddo ;
enddo ;
enddo
1019 do j=jsq,jeq ;
do k=1,nz ;
do i=is,ie
1020 bt_force_v(i,j) = bt_force_v(i,j) + wt_v(i,j,k) * bc_accel_v(i,j,k)
1021 enddo ;
enddo ;
enddo
1027 do j=js,je ;
do i=is-1,ie ; uhbt(i,j) = 0.0 ; ubt(i,j) = 0.0 ;
enddo ;
enddo
1029 do j=js-1,je ;
do i=is,ie ; vhbt(i,j) = 0.0 ; vbt(i,j) = 0.0 ;
enddo ;
enddo
1030 if (cs%visc_rem_u_uh0)
then
1032 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
1033 uhbt(i,j) = uhbt(i,j) + uh0(i,j,k)
1034 ubt(i,j) = ubt(i,j) + wt_u(i,j,k) * u_uh0(i,j,k)
1035 enddo ;
enddo ;
enddo
1037 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
1038 vhbt(i,j) = vhbt(i,j) + vh0(i,j,k)
1039 vbt(i,j) = vbt(i,j) + wt_v(i,j,k) * v_vh0(i,j,k)
1040 enddo ;
enddo ;
enddo
1043 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
1044 uhbt(i,j) = uhbt(i,j) + uh0(i,j,k)
1045 ubt(i,j) = ubt(i,j) + cs%frhatu(i,j,k) * u_uh0(i,j,k)
1046 enddo ;
enddo ;
enddo
1048 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
1049 vhbt(i,j) = vhbt(i,j) + vh0(i,j,k)
1050 vbt(i,j) = vbt(i,j) + cs%frhatv(i,j,k) * v_vh0(i,j,k)
1051 enddo ;
enddo ;
enddo
1053 if (use_bt_cont)
then
1054 if (cs%adjust_BT_cont)
then
1061 call pass_vector(ubt, vbt, cs%BT_Domain, complete=.false., halo=1+ievf-ie)
1062 call pass_vector(uhbt, vhbt, cs%BT_Domain, complete=.true., halo=1+ievf-ie)
1067 g, us, ms, 1+ievf-ie)
1070 do j=js,je ;
do i=is-1,ie
1071 uhbt0(i,j) = uhbt(i,j) -
find_uhbt(ubt(i,j), btcl_u(i,j), us)
1074 do j=js-1,je ;
do i=is,ie
1075 vhbt0(i,j) = vhbt(i,j) -
find_vhbt(vbt(i,j), btcl_v(i,j), us)
1079 do j=js,je ;
do i=is-1,ie
1080 uhbt0(i,j) = uhbt(i,j) - datu(i,j)*ubt(i,j)
1083 do j=js-1,je ;
do i=is,ie
1084 vhbt0(i,j) = vhbt(i,j) - datv(i,j)*vbt(i,j)
1087 if (cs%BT_OBC%apply_u_OBCs)
then
1089 do j=js,je ;
do i=is-1,ie ;
if (obc%segnum_u(i,j) /= obc_none)
then
1091 endif ;
enddo ;
enddo
1093 if (cs%BT_OBC%apply_v_OBCs)
then
1095 do j=js-1,je ;
do i=is,ie ;
if (obc%segnum_v(i,j) /= obc_none)
then
1097 endif ;
enddo ;
enddo
1103 do j=jsvf-1,jevf+1 ;
do i=isvf-2,ievf+1
1104 ubt(i,j) = 0.0 ; uhbt(i,j) = 0.0 ; u_accel_bt(i,j) = 0.0
1107 do j=jsvf-2,jevf+1 ;
do i=isvf-1,ievf+1
1108 vbt(i,j) = 0.0 ; vhbt(i,j) = 0.0 ; v_accel_bt(i,j) = 0.0
1111 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
1112 ubt(i,j) = ubt(i,j) + wt_u(i,j,k) * u_in(i,j,k)
1113 enddo ;
enddo ;
enddo
1115 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
1116 vbt(i,j) = vbt(i,j) + wt_v(i,j,k) * v_in(i,j,k)
1117 enddo ;
enddo ;
enddo
1119 do j=js,je ;
do i=is-1,ie
1120 if (abs(ubt(i,j)) < cs%vel_underflow) ubt(i,j) = 0.0
1123 do j=js-1,je ;
do i=is,ie
1124 if (abs(vbt(i,j)) < cs%vel_underflow) vbt(i,j) = 0.0
1127 if (apply_obcs)
then
1128 ubt_first(:,:) = ubt(:,:) ; vbt_first(:,:) = vbt(:,:)
1131 if (cs%gradual_BT_ICs)
then
1133 do j=js,je ;
do i=is-1,ie
1134 bt_force_u(i,j) = bt_force_u(i,j) + (ubt(i,j) - cs%ubt_IC(i,j)) * idt
1135 ubt(i,j) = cs%ubt_IC(i,j)
1136 if (abs(ubt(i,j)) < cs%vel_underflow) ubt(i,j) = 0.0
1139 do j=js-1,je ;
do i=is,ie
1140 bt_force_v(i,j) = bt_force_v(i,j) + (vbt(i,j) - cs%vbt_IC(i,j)) * idt
1141 vbt(i,j) = cs%vbt_IC(i,j)
1142 if (abs(vbt(i,j)) < cs%vel_underflow) vbt(i,j) = 0.0
1146 if ((isq > is-1) .or. (jsq > js-1))
then
1151 tmp_u(:,:) = 0.0 ; tmp_v(:,:) = 0.0
1152 do j=js,je ;
do i=isq,ieq ; tmp_u(i,j) = bt_force_u(i,j) ;
enddo ;
enddo
1153 do j=jsq,jeq ;
do i=is,ie ; tmp_v(i,j) = bt_force_v(i,j) ;
enddo ;
enddo
1154 if (nonblock_setup)
then
1157 call do_group_pass(cs%pass_tmp_uv, g%Domain)
1158 do j=jsd,jed ;
do i=isdb,iedb ; bt_force_u(i,j) = tmp_u(i,j) ;
enddo ;
enddo
1159 do j=jsdb,jedb ;
do i=isd,ied ; bt_force_v(i,j) = tmp_v(i,j) ;
enddo ;
enddo
1165 if (nonblock_setup)
then
1176 do j=jsvf-1,jevf ;
do i=isvf-1,ievf+1
1177 if (cs%Sadourny)
then
1178 amer(i-1,j) = dcor_u(i-1,j) * q(i-1,j)
1179 bmer(i,j) = dcor_u(i,j) * q(i,j)
1180 cmer(i,j+1) = dcor_u(i,j+1) * q(i,j)
1181 dmer(i-1,j+1) = dcor_u(i-1,j+1) * q(i-1,j)
1183 amer(i-1,j) = dcor_u(i-1,j) * &
1184 ((q(i,j) + q(i-1,j-1)) + q(i-1,j)) / 3.0
1185 bmer(i,j) = dcor_u(i,j) * &
1186 (q(i,j) + (q(i-1,j) + q(i,j-1))) / 3.0
1187 cmer(i,j+1) = dcor_u(i,j+1) * &
1188 (q(i,j) + (q(i-1,j) + q(i,j+1))) / 3.0
1189 dmer(i-1,j+1) = dcor_u(i-1,j+1) * &
1190 ((q(i,j) + q(i-1,j+1)) + q(i-1,j)) / 3.0
1195 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf
1196 if (cs%Sadourny)
then
1197 azon(i,j) = dcor_v(i+1,j) * q(i,j)
1198 bzon(i,j) = dcor_v(i,j) * q(i,j)
1199 czon(i,j) = dcor_v(i,j-1) * q(i,j-1)
1200 dzon(i,j) = dcor_v(i+1,j-1) * q(i,j-1)
1202 azon(i,j) = dcor_v(i+1,j) * &
1203 (q(i,j) + (q(i+1,j) + q(i,j-1))) / 3.0
1204 bzon(i,j) = dcor_v(i,j) * &
1205 (q(i,j) + (q(i-1,j) + q(i,j-1))) / 3.0
1206 czon(i,j) = dcor_v(i,j-1) * &
1207 ((q(i,j) + q(i-1,j-1)) + q(i,j-1)) / 3.0
1208 dzon(i,j) = dcor_v(i+1,j-1) * &
1209 ((q(i,j) + q(i+1,j-1)) + q(i,j-1)) / 3.0
1216 if (nonblock_setup)
then
1217 if ((isq > is-1) .or. (jsq > js-1))
then
1219 do j=jsd,jed ;
do i=isdb,iedb ; bt_force_u(i,j) = tmp_u(i,j) ;
enddo ;
enddo
1220 do j=jsdb,jedb ;
do i=isd,ied ; bt_force_v(i,j) = tmp_v(i,j) ;
enddo ;
enddo
1225 call do_group_pass(cs%pass_gtot, cs%BT_Domain)
1226 call do_group_pass(cs%pass_ubt_Cor, g%Domain)
1230 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf+1
1231 if (cs%ua_polarity(i,j) < 0.0)
call swap(gtot_e(i,j), gtot_w(i,j))
1232 if (cs%va_polarity(i,j) < 0.0)
call swap(gtot_n(i,j), gtot_s(i,j))
1236 do j=js,je ;
do i=is-1,ie
1238 ((azon(i,j) * vbt_cor(i+1,j) + czon(i,j) * vbt_cor(i ,j-1)) + &
1239 (bzon(i,j) * vbt_cor(i ,j) + dzon(i,j) * vbt_cor(i+1,j-1)))
1242 do j=js-1,je ;
do i=is,ie
1243 cor_ref_v(i,j) = -1.0 * &
1244 ((amer(i-1,j) * ubt_cor(i-1,j) + cmer(i ,j+1) * ubt_cor(i ,j+1)) + &
1245 (bmer(i ,j) * ubt_cor(i ,j) + dmer(i-1,j+1) * ubt_cor(i-1,j+1)))
1249 if (nonblock_setup)
then
1250 if (.not.use_bt_cont) &
1268 do j=js-1,je+1 ;
do i=is-1,ie ; av_rem_u(i,j) = 0.0 ;
enddo ;
enddo
1270 do j=js-1,je ;
do i=is-1,ie+1 ; av_rem_v(i,j) = 0.0 ;
enddo ;
enddo
1272 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
1273 av_rem_u(i,j) = av_rem_u(i,j) + cs%frhatu(i,j,k) * visc_rem_u(i,j,k)
1274 enddo ;
enddo ;
enddo
1276 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
1277 av_rem_v(i,j) = av_rem_v(i,j) + cs%frhatv(i,j,k) * visc_rem_v(i,j,k)
1278 enddo ;
enddo ;
enddo
1279 if (cs%strong_drag)
then
1281 do j=js,je ;
do i=is-1,ie
1282 bt_rem_u(i,j) = g%mask2dCu(i,j) * &
1283 ((nstep * av_rem_u(i,j)) / (1.0 + (nstep-1)*av_rem_u(i,j)))
1286 do j=js-1,je ;
do i=is,ie
1287 bt_rem_v(i,j) = g%mask2dCv(i,j) * &
1288 ((nstep * av_rem_v(i,j)) / (1.0 + (nstep-1)*av_rem_v(i,j)))
1292 do j=js,je ;
do i=is-1,ie
1294 if (g%mask2dCu(i,j) * av_rem_u(i,j) > 0.0) &
1295 bt_rem_u(i,j) = g%mask2dCu(i,j) * (av_rem_u(i,j)**instep)
1298 do j=js-1,je ;
do i=is,ie
1300 if (g%mask2dCv(i,j) * av_rem_v(i,j) > 0.0) &
1301 bt_rem_v(i,j) = g%mask2dCv(i,j) * (av_rem_v(i,j)**instep)
1304 if (cs%linear_wave_drag)
then
1306 do j=js,je ;
do i=is-1,ie ;
if (cs%lin_drag_u(i,j) > 0.0)
then
1307 htot = 0.5 * (eta(i,j) + eta(i+1,j))
1308 if (gv%Boussinesq) &
1309 htot = htot + 0.5*gv%Z_to_H * (cs%bathyT(i,j) + cs%bathyT(i+1,j))
1310 bt_rem_u(i,j) = bt_rem_u(i,j) * (htot / (htot + cs%lin_drag_u(i,j) * dtbt))
1312 rayleigh_u(i,j) = cs%lin_drag_u(i,j) / htot
1313 endif ;
enddo ;
enddo
1315 do j=js-1,je ;
do i=is,ie ;
if (cs%lin_drag_v(i,j) > 0.0)
then
1316 htot = 0.5 * (eta(i,j) + eta(i,j+1))
1317 if (gv%Boussinesq) &
1318 htot = htot + 0.5*gv%Z_to_H * (cs%bathyT(i,j) + cs%bathyT(i+1,j+1))
1319 bt_rem_v(i,j) = bt_rem_v(i,j) * (htot / (htot + cs%lin_drag_v(i,j) * dtbt))
1321 rayleigh_v(i,j) = cs%lin_drag_v(i,j) / htot
1322 endif ;
enddo ;
enddo
1326 if (find_etaav)
then
1328 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf+1
1329 eta_sum(i,j) = 0.0 ; eta_wtd(i,j) = 0.0
1333 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf+1
1338 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf
1339 ubt_sum(i,j) = 0.0 ; uhbt_sum(i,j) = 0.0
1340 pfu_bt_sum(i,j) = 0.0 ; coru_bt_sum(i,j) = 0.0
1341 ubt_wtd(i,j) = 0.0 ; ubt_trans(i,j) = 0.0
1344 do j=jsvf-1,jevf ;
do i=isvf-1,ievf+1
1345 vbt_sum(i,j) = 0.0 ; vhbt_sum(i,j) = 0.0
1346 pfv_bt_sum(i,j) = 0.0 ; corv_bt_sum(i,j) = 0.0
1347 vbt_wtd(i,j) = 0.0 ; vbt_trans(i,j) = 0.0
1352 do j=jsvf-1,jevf+1;
do i=isvf-1,ievf+1 ; eta_src(i,j) = 0.0 ;
enddo ;
enddo
1353 if (cs%bound_BT_corr)
then ;
if (use_bt_cont .and. cs%BT_cont_bounds)
then
1354 do j=js,je ;
do i=is,ie ;
if (g%mask2dT(i,j) > 0.0)
then
1355 if (cs%eta_cor(i,j) > 0.0)
then
1359 u_max_cor = g%dxT(i,j) * (cs%maxCFL_BT_cont*idt)
1360 v_max_cor = g%dyT(i,j) * (cs%maxCFL_BT_cont*idt)
1361 eta_cor_max = dt * (cs%IareaT(i,j) * &
1362 (((
find_uhbt(u_max_cor, btcl_u(i,j), us) + uhbt0(i,j)) - &
1363 (
find_uhbt(-u_max_cor, btcl_u(i-1,j), us) + uhbt0(i-1,j))) + &
1364 ((
find_vhbt(v_max_cor, btcl_v(i,j), us) + vhbt0(i,j)) - &
1365 (
find_vhbt(-v_max_cor, btcl_v(i,j-1), us) + vhbt0(i,j-1))) ))
1366 cs%eta_cor(i,j) = min(cs%eta_cor(i,j), max(0.0, eta_cor_max))
1371 if (gv%Boussinesq) htot = cs%bathyT(i,j)*gv%Z_to_H + eta(i,j)
1373 cs%eta_cor(i,j) = max(cs%eta_cor(i,j), -max(0.0,htot))
1375 endif ;
enddo ;
enddo
1376 else ;
do j=js,je ;
do i=is,ie
1377 if (abs(cs%eta_cor(i,j)) > dt*cs%eta_cor_bound(i,j)) &
1378 cs%eta_cor(i,j) = sign(dt*cs%eta_cor_bound(i,j), cs%eta_cor(i,j))
1379 enddo ;
enddo ;
endif ;
endif
1381 do j=js,je ;
do i=is,ie
1382 eta_src(i,j) = g%mask2dT(i,j) * (instep * cs%eta_cor(i,j))
1386 if (cs%dynamic_psurf)
then
1387 ice_is_rigid = (
associated(forces%rigidity_ice_u) .and. &
1388 associated(forces%rigidity_ice_v))
1389 h_min_dyn = gv%Z_to_H * cs%Dmin_dyn_psurf
1390 if (ice_is_rigid .and. use_bt_cont) &
1392 if (ice_is_rigid)
then
1394 do j=js,je ;
do i=is,ie
1400 idt_max2 = 0.5 * (dgeo_de * (1.0 + 2.0*bebt)) * (g%IareaT(i,j) * &
1401 ((gtot_e(i,j) * (datu(i,j)*g%IdxCu(i,j)) + &
1402 gtot_w(i,j) * (datu(i-1,j)*g%IdxCu(i-1,j))) + &
1403 (gtot_n(i,j) * (datv(i,j)*g%IdyCv(i,j)) + &
1404 gtot_s(i,j) * (datv(i,j-1)*g%IdyCv(i,j-1)))) + &
1405 ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
1406 (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2)))
1407 h_eff_dx2 = max(h_min_dyn * ((g%IdxT(i,j))**2 + (g%IdyT(i,j))**2), &
1409 ((datu(i,j)*g%IdxCu(i,j) + datu(i-1,j)*g%IdxCu(i-1,j)) + &
1410 (datv(i,j)*g%IdyCv(i,j) + datv(i,j-1)*g%IdyCv(i,j-1)) ) )
1411 dyn_coef_max = cs%const_dyn_psurf * max(0.0, 1.0 - dtbt**2 * idt_max2) / &
1412 (dtbt**2 * h_eff_dx2)
1415 ice_strength = us%m_to_L**4*us%Z_to_m*us%T_to_s* &
1416 ((forces%rigidity_ice_u(i,j) + forces%rigidity_ice_u(i-1,j)) + &
1417 (forces%rigidity_ice_v(i,j) + forces%rigidity_ice_v(i,j-1))) / &
1418 (cs%ice_strength_length**2 * dtbt)
1421 dyn_coef_eta(i,j) = min(dyn_coef_max, ice_strength * gv%H_to_Z)
1422 enddo ;
enddo ;
endif
1427 if (nonblock_setup)
then
1431 call do_group_pass(cs%pass_eta_bt_rem, cs%BT_Domain)
1432 if (.not.use_bt_cont) &
1433 call do_group_pass(cs%pass_Dat_uv, cs%BT_Domain)
1434 call do_group_pass(cs%pass_force_hbt0_Cor_ref, cs%BT_Domain)
1440 if (nonblock_setup)
then
1444 if (.not.use_bt_cont) &
1454 call uvchksum(
"BT [uv]hbt", uhbt, vhbt, cs%debug_BT_HI, haloshift=0, &
1455 scale=us%s_to_T*us%L_to_m**2*gv%H_to_m)
1456 call uvchksum(
"BT Initial [uv]bt", ubt, vbt, cs%debug_BT_HI, haloshift=0, scale=us%L_T_to_m_s)
1457 call hchksum(eta,
"BT Initial eta", cs%debug_BT_HI, haloshift=0, scale=gv%H_to_m)
1458 call uvchksum(
"BT BT_force_[uv]", bt_force_u, bt_force_v, &
1459 cs%debug_BT_HI, haloshift=0, scale=us%L_T2_to_m_s2)
1460 if (interp_eta_pf)
then
1461 call hchksum(eta_pf_1,
"BT eta_PF_1",cs%debug_BT_HI,haloshift=0, scale=gv%H_to_m)
1462 call hchksum(d_eta_pf,
"BT d_eta_PF",cs%debug_BT_HI,haloshift=0, scale=gv%H_to_m)
1464 call hchksum(eta_pf,
"BT eta_PF",cs%debug_BT_HI,haloshift=0, scale=gv%H_to_m)
1465 call hchksum(eta_pf_in,
"BT eta_PF_in",g%HI,haloshift=0, scale=gv%H_to_m)
1467 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)
1468 call uvchksum(
"BT [uv]hbt0", uhbt0, vhbt0, cs%debug_BT_HI, haloshift=0, &
1469 scale=us%L_to_m**2*us%s_to_T*gv%H_to_m)
1470 if (.not. use_bt_cont)
then
1471 call uvchksum(
"BT Dat[uv]", datu, datv, cs%debug_BT_HI, haloshift=1, scale=us%L_to_m*gv%H_to_m)
1473 call uvchksum(
"BT wt_[uv]", wt_u, wt_v, g%HI, 0, .true., .true.)
1474 call uvchksum(
"BT frhat[uv]", cs%frhatu, cs%frhatv, g%HI, 0, .true., .true.)
1475 call uvchksum(
"BT bc_accel_[uv]", bc_accel_u, bc_accel_v, g%HI, haloshift=0, scale=us%L_T2_to_m_s2)
1476 call uvchksum(
"BT IDat[uv]", cs%IDatu, cs%IDatv, g%HI, haloshift=0, scale=us%m_to_Z)
1477 call uvchksum(
"BT visc_rem_[uv]", visc_rem_u, visc_rem_v, g%HI, haloshift=1)
1480 if (query_averaging_enabled(cs%diag))
then
1481 if (cs%id_eta_st > 0)
call post_data(cs%id_eta_st, eta(isd:ied,jsd:jed), cs%diag)
1482 if (cs%id_ubt_st > 0)
call post_data(cs%id_ubt_st, ubt(isdb:iedb,jsd:jed), cs%diag)
1483 if (cs%id_vbt_st > 0)
call post_data(cs%id_vbt_st, vbt(isd:ied,jsdb:jedb), cs%diag)
1489 if (project_velocity)
then ; eta_pf_bt => eta ;
else ; eta_pf_bt => eta_pred ;
endif
1491 if (cs%dt_bt_filter >= 0.0)
then
1492 dt_filt = 0.5 * max(0.0, min(cs%dt_bt_filter, 2.0*dt))
1494 dt_filt = 0.5 * max(0.0, dt * min(-cs%dt_bt_filter, 2.0))
1496 nfilter = ceiling(dt_filt / dtbt)
1498 if (nstep+nfilter==0 )
call mom_error(fatal, &
1499 "btstep: number of barotropic step (nstep+nfilter) is 0")
1502 sum_wt_vel = 0.0 ; sum_wt_eta = 0.0 ; sum_wt_accel = 0.0 ; sum_wt_trans = 0.0
1503 allocate(wt_vel(nstep+nfilter)) ;
allocate(wt_eta(nstep+nfilter))
1504 allocate(wt_trans(nstep+nfilter+1)) ;
allocate(wt_accel(nstep+nfilter+1))
1505 allocate(wt_accel2(nstep+nfilter+1))
1506 do n=1,nstep+nfilter
1508 if ( (n==nstep) .or. (dt_filt - abs(n-nstep)*dtbt >= 0.0))
then
1509 wt_vel(n) = 1.0 ; wt_eta(n) = 1.0
1510 elseif (dtbt + dt_filt - abs(n-nstep)*dtbt > 0.0)
then
1511 wt_vel(n) = 1.0 + (dt_filt / dtbt) - abs(n-nstep) ; wt_eta(n) = wt_vel(n)
1513 wt_vel(n) = 0.0 ; wt_eta(n) = 0.0
1519 sum_wt_vel = sum_wt_vel + wt_vel(n) ; sum_wt_eta = sum_wt_eta + wt_eta(n)
1521 wt_trans(nstep+nfilter+1) = 0.0 ; wt_accel(nstep+nfilter+1) = 0.0
1522 do n=nstep+nfilter,1,-1
1523 wt_trans(n) = wt_trans(n+1) + wt_eta(n)
1524 wt_accel(n) = wt_accel(n+1) + wt_vel(n)
1525 sum_wt_accel = sum_wt_accel + wt_accel(n) ; sum_wt_trans = sum_wt_trans + wt_trans(n)
1528 i_sum_wt_vel = 1.0 / sum_wt_vel ; i_sum_wt_accel = 1.0 / sum_wt_accel
1529 i_sum_wt_eta = 1.0 / sum_wt_eta ; i_sum_wt_trans = 1.0 / sum_wt_trans
1530 do n=1,nstep+nfilter
1531 wt_vel(n) = wt_vel(n) * i_sum_wt_vel
1532 wt_accel2(n) = wt_accel(n)
1533 wt_accel(n) = wt_accel(n) * i_sum_wt_accel
1534 wt_eta(n) = wt_eta(n) * i_sum_wt_eta
1539 sum_wt_vel = 0.0 ; sum_wt_eta = 0.0 ; sum_wt_accel = 0.0 ; sum_wt_trans = 0.0
1542 isv=is ; iev=ie ; jsv=js ; jev=je
1543 do n=1,nstep+nfilter
1545 sum_wt_vel = sum_wt_vel + wt_vel(n)
1546 sum_wt_eta = sum_wt_eta + wt_eta(n)
1547 sum_wt_accel = sum_wt_accel + wt_accel2(n)
1548 sum_wt_trans = sum_wt_trans + wt_trans(n)
1550 if (cs%clip_velocity)
then
1551 do j=jsv,jev ;
do i=isv-1,iev
1552 if ((ubt(i,j) * (dt * g%dy_Cu(i,j))) * g%IareaT(i+1,j) < -cs%CFL_trunc)
then
1554 ubt(i,j) = (-0.95*cs%CFL_trunc) * (g%areaT(i+1,j) / (dt * g%dy_Cu(i,j)))
1555 elseif ((ubt(i,j) * (dt * g%dy_Cu(i,j))) * g%IareaT(i,j) > cs%CFL_trunc)
then
1557 ubt(i,j) = (0.95*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dy_Cu(i,j)))
1560 do j=jsv-1,jev ;
do i=isv,iev
1561 if ((vbt(i,j) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j+1) < -cs%CFL_trunc)
then
1563 vbt(i,j) = (-0.9*cs%CFL_trunc) * (g%areaT(i,j+1) / (dt * g%dx_Cv(i,j)))
1564 elseif ((vbt(i,j) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j) > cs%CFL_trunc)
then
1566 vbt(i,j) = (0.9*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dx_Cv(i,j)))
1571 if ((iev - stencil < ie) .or. (jev - stencil < je))
then
1574 isv = isvf ; iev = ievf ; jsv = jsvf ; jev = jevf
1577 isv = isv+stencil ; iev = iev-stencil
1578 jsv = jsv+stencil ; jev = jev-stencil
1581 if ((.not.use_bt_cont) .and. cs%Nonlinear_continuity .and. &
1582 (cs%Nonlin_cont_update_period > 0))
then
1583 if ((n>1) .and. (mod(n-1,cs%Nonlin_cont_update_period) == 0)) &
1588 if (cs%dynamic_psurf .or. .not.project_velocity)
then
1589 if (use_bt_cont)
then
1591 do j=jsv-1,jev+1 ;
do i=isv-2,iev+1
1592 uhbt(i,j) =
find_uhbt(ubt(i,j), btcl_u(i,j), us) + uhbt0(i,j)
1595 do j=jsv-2,jev+1 ;
do i=isv-1,iev+1
1596 vhbt(i,j) =
find_vhbt(vbt(i,j), btcl_v(i,j), us) + vhbt0(i,j)
1599 do j=jsv-1,jev+1 ;
do i=isv-1,iev+1
1600 eta_pred(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * cs%IareaT(i,j)) * &
1601 ((uhbt(i-1,j) - uhbt(i,j)) + (vhbt(i,j-1) - vhbt(i,j)))
1605 do j=jsv-1,jev+1 ;
do i=isv-1,iev+1
1606 eta_pred(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * cs%IareaT(i,j)) * &
1607 (((datu(i-1,j)*ubt(i-1,j) + uhbt0(i-1,j)) - &
1608 (datu(i,j)*ubt(i,j) + uhbt0(i,j))) + &
1609 ((datv(i,j-1)*vbt(i,j-1) + vhbt0(i,j-1)) - &
1610 (datv(i,j)*vbt(i,j) + vhbt0(i,j))))
1614 if (cs%dynamic_psurf)
then
1616 do j=jsv-1,jev+1 ;
do i=isv-1,iev+1
1617 p_surf_dyn(i,j) = dyn_coef_eta(i,j) * (eta_pred(i,j) - eta(i,j))
1625 if (find_etaav)
then
1627 do j=js,je ;
do i=is,ie
1628 eta_sum(i,j) = eta_sum(i,j) + wt_accel2(n) * eta_pf_bt(i,j)
1632 if (interp_eta_pf)
then
1635 do j=jsv-1,jev+1 ;
do i=isv-1,iev+1
1636 eta_pf(i,j) = eta_pf_1(i,j) + wt_end*d_eta_pf(i,j)
1640 if (apply_obc_flather .or. apply_obc_open)
then
1642 do j=jsv,jev ;
do i=isv-2,iev+1
1643 ubt_old(i,j) = ubt(i,j)
1646 do j=jsv-2,jev+1 ;
do i=isv,iev
1647 vbt_old(i,j) = vbt(i,j)
1652 if (apply_obcs)
then
1653 if (mod(n+g%first_direction,2)==1)
then
1659 if (cs%BT_OBC%apply_u_OBCs)
then
1661 do j=jsv-joff,jev+joff ;
do i=isv-1,iev
1662 ubt_prev(i,j) = ubt(i,j) ; uhbt_prev(i,j) = uhbt(i,j)
1663 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)
1667 if (cs%BT_OBC%apply_v_OBCs)
then
1669 do j=jsv-1,jev ;
do i=isv-ioff,iev+ioff
1670 vbt_prev(i,j) = vbt(i,j) ; vhbt_prev(i,j) = vhbt(i,j)
1671 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)
1677 if (mod(n+g%first_direction,2)==1)
then
1680 do j=jsv-1,jev ;
do i=isv-1,iev+1
1681 cor_v(i,j) = -1.0*((amer(i-1,j) * ubt(i-1,j) + cmer(i,j+1) * ubt(i,j+1)) + &
1682 (bmer(i,j) * ubt(i,j) + dmer(i-1,j+1) * ubt(i-1,j+1))) - cor_ref_v(i,j)
1683 pfv(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_n(i,j) - &
1684 (eta_pf_bt(i,j+1)-eta_pf(i,j+1))*gtot_s(i,j+1)) * &
1685 dgeo_de * cs%IdyCv(i,j)
1687 if (cs%dynamic_psurf)
then
1689 do j=jsv-1,jev ;
do i=isv-1,iev+1
1690 pfv(i,j) = pfv(i,j) + (p_surf_dyn(i,j) - p_surf_dyn(i,j+1)) * cs%IdyCv(i,j)
1694 if (cs%BT_OBC%apply_v_OBCs)
then
1696 do j=jsv-1,jev ;
do i=isv-1,iev+1 ;
if (obc%segnum_v(i,j) /= obc_none)
then
1698 endif ;
enddo ;
enddo
1701 do j=jsv-1,jev ;
do i=isv-1,iev+1
1703 vbt(i,j) = bt_rem_v(i,j) * (vbt(i,j) + &
1704 dtbt * ((bt_force_v(i,j) + cor_v(i,j)) + pfv(i,j)))
1705 vbt_trans(i,j) = trans_wt1*vbt(i,j) + trans_wt2*vel_prev
1707 if (cs%linear_wave_drag)
then
1708 v_accel_bt(i,j) = v_accel_bt(i,j) + wt_accel(n) * &
1709 ((cor_v(i,j) + pfv(i,j)) - vbt(i,j)*rayleigh_v(i,j))
1711 v_accel_bt(i,j) = v_accel_bt(i,j) + wt_accel(n) * (cor_v(i,j) + pfv(i,j))
1715 if (use_bt_cont)
then
1717 do j=jsv-1,jev ;
do i=isv-1,iev+1
1718 vhbt(i,j) =
find_vhbt(vbt_trans(i,j), btcl_v(i,j), us) + vhbt0(i,j)
1722 do j=jsv-1,jev ;
do i=isv-1,iev+1
1723 vhbt(i,j) = datv(i,j)*vbt_trans(i,j) + vhbt0(i,j)
1726 if (cs%BT_OBC%apply_v_OBCs)
then
1728 do j=jsv-1,jev ;
do i=isv-1,iev+1 ;
if (obc%segnum_v(i,j) /= obc_none)
then
1729 vbt(i,j) = vbt_prev(i,j) ; vhbt(i,j) = vhbt_prev(i,j)
1730 endif ;
enddo ;
enddo
1734 do j=jsv,jev ;
do i=isv-1,iev
1735 cor_u(i,j) = ((azon(i,j) * vbt(i+1,j) + czon(i,j) * vbt(i,j-1)) + &
1736 (bzon(i,j) * vbt(i,j) + dzon(i,j) * vbt(i+1,j-1))) - &
1738 pfu(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_e(i,j) - &
1739 (eta_pf_bt(i+1,j)-eta_pf(i+1,j))*gtot_w(i+1,j)) * &
1740 dgeo_de * cs%IdxCu(i,j)
1743 if (cs%dynamic_psurf)
then
1745 do j=jsv,jev ;
do i=isv-1,iev
1746 pfu(i,j) = pfu(i,j) + (p_surf_dyn(i,j) - p_surf_dyn(i+1,j)) * cs%IdxCu(i,j)
1750 if (cs%BT_OBC%apply_u_OBCs)
then
1752 do j=jsv,jev ;
do i=isv-1,iev ;
if (obc%segnum_u(i,j) /= obc_none)
then
1754 endif ;
enddo ;
enddo
1757 do j=jsv,jev ;
do i=isv-1,iev
1759 ubt(i,j) = bt_rem_u(i,j) * (ubt(i,j) + &
1760 dtbt * ((bt_force_u(i,j) + cor_u(i,j)) + pfu(i,j)))
1761 if (abs(ubt(i,j)) < cs%vel_underflow) ubt(i,j) = 0.0
1762 ubt_trans(i,j) = trans_wt1*ubt(i,j) + trans_wt2*vel_prev
1764 if (cs%linear_wave_drag)
then
1765 u_accel_bt(i,j) = u_accel_bt(i,j) + wt_accel(n) * &
1766 ((cor_u(i,j) + pfu(i,j)) - ubt(i,j)*rayleigh_u(i,j))
1768 u_accel_bt(i,j) = u_accel_bt(i,j) + wt_accel(n) * (cor_u(i,j) + pfu(i,j))
1772 if (use_bt_cont)
then
1774 do j=jsv,jev ;
do i=isv-1,iev
1775 uhbt(i,j) =
find_uhbt(ubt_trans(i,j), btcl_u(i,j), us) + uhbt0(i,j)
1779 do j=jsv,jev ;
do i=isv-1,iev
1780 uhbt(i,j) = datu(i,j)*ubt_trans(i,j) + uhbt0(i,j)
1783 if (cs%BT_OBC%apply_u_OBCs)
then
1785 do j=jsv,jev ;
do i=isv-1,iev ;
if (obc%segnum_u(i,j) /= obc_none)
then
1786 ubt(i,j) = ubt_prev(i,j); uhbt(i,j) = uhbt_prev(i,j)
1787 endif ;
enddo ;
enddo
1792 do j=jsv-1,jev+1 ;
do i=isv-1,iev
1793 cor_u(i,j) = ((azon(i,j) * vbt(i+1,j) + czon(i,j) * vbt(i,j-1)) + &
1794 (bzon(i,j) * vbt(i,j) + dzon(i,j) * vbt(i+1,j-1))) - &
1796 pfu(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_e(i,j) - &
1797 (eta_pf_bt(i+1,j)-eta_pf(i+1,j))*gtot_w(i+1,j)) * &
1798 dgeo_de * cs%IdxCu(i,j)
1801 if (cs%dynamic_psurf)
then
1803 do j=jsv-1,jev+1 ;
do i=isv-1,iev
1804 pfu(i,j) = pfu(i,j) + (p_surf_dyn(i,j) - p_surf_dyn(i+1,j)) * cs%IdxCu(i,j)
1808 if (cs%BT_OBC%apply_u_OBCs)
then
1810 do j=jsv,jev ;
do i=isv-1,iev ;
if (obc%segnum_u(i,j) /= obc_none)
then
1812 endif ;
enddo ;
enddo
1816 do j=jsv-1,jev+1 ;
do i=isv-1,iev
1818 ubt(i,j) = bt_rem_u(i,j) * (ubt(i,j) + &
1819 dtbt * ((bt_force_u(i,j) + cor_u(i,j)) + pfu(i,j)))
1820 if (abs(ubt(i,j)) < cs%vel_underflow) ubt(i,j) = 0.0
1821 ubt_trans(i,j) = trans_wt1*ubt(i,j) + trans_wt2*vel_prev
1823 if (cs%linear_wave_drag)
then
1824 u_accel_bt(i,j) = u_accel_bt(i,j) + wt_accel(n) * &
1825 ((cor_u(i,j) + pfu(i,j)) - ubt(i,j)*rayleigh_u(i,j))
1827 u_accel_bt(i,j) = u_accel_bt(i,j) + wt_accel(n) * (cor_u(i,j) + pfu(i,j))
1831 if (use_bt_cont)
then
1833 do j=jsv-1,jev+1 ;
do i=isv-1,iev
1834 uhbt(i,j) =
find_uhbt(ubt_trans(i,j), btcl_u(i,j), us) + uhbt0(i,j)
1838 do j=jsv-1,jev+1 ;
do i=isv-1,iev
1839 uhbt(i,j) = datu(i,j)*ubt_trans(i,j) + uhbt0(i,j)
1842 if (cs%BT_OBC%apply_u_OBCs)
then
1844 do j=jsv-1,jev+1 ;
do i=isv-1,iev ;
if (obc%segnum_u(i,j) /= obc_none)
then
1845 ubt(i,j) = ubt_prev(i,j); uhbt(i,j) = uhbt_prev(i,j)
1846 endif ;
enddo ;
enddo
1850 if (cs%use_old_coriolis_bracket_bug)
then
1852 do j=jsv-1,jev ;
do i=isv,iev
1853 cor_v(i,j) = -1.0*((amer(i-1,j) * ubt(i-1,j) + bmer(i,j) * ubt(i,j)) + &
1854 (cmer(i,j+1) * ubt(i,j+1) + dmer(i-1,j+1) * ubt(i-1,j+1))) - cor_ref_v(i,j)
1855 pfv(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_n(i,j) - &
1856 (eta_pf_bt(i,j+1)-eta_pf(i,j+1))*gtot_s(i,j+1)) * &
1857 dgeo_de * cs%IdyCv(i,j)
1861 do j=jsv-1,jev ;
do i=isv,iev
1862 cor_v(i,j) = -1.0*((amer(i-1,j) * ubt(i-1,j) + cmer(i,j+1) * ubt(i,j+1)) + &
1863 (bmer(i,j) * ubt(i,j) + dmer(i-1,j+1) * ubt(i-1,j+1))) - cor_ref_v(i,j)
1864 pfv(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_n(i,j) - &
1865 (eta_pf_bt(i,j+1)-eta_pf(i,j+1))*gtot_s(i,j+1)) * &
1866 dgeo_de * cs%IdyCv(i,j)
1870 if (cs%dynamic_psurf)
then
1872 do j=jsv-1,jev ;
do i=isv,iev
1873 pfv(i,j) = pfv(i,j) + (p_surf_dyn(i,j) - p_surf_dyn(i,j+1)) * cs%IdyCv(i,j)
1877 if (cs%BT_OBC%apply_v_OBCs)
then
1879 do j=jsv-1,jev ;
do i=isv-1,iev+1 ;
if (obc%segnum_v(i,j) /= obc_none)
then
1881 endif ;
enddo ;
enddo
1885 do j=jsv-1,jev ;
do i=isv,iev
1887 vbt(i,j) = bt_rem_v(i,j) * (vbt(i,j) + &
1888 dtbt * ((bt_force_v(i,j) + cor_v(i,j)) + pfv(i,j)))
1889 if (abs(vbt(i,j)) < cs%vel_underflow) vbt(i,j) = 0.0
1890 vbt_trans(i,j) = trans_wt1*vbt(i,j) + trans_wt2*vel_prev
1892 if (cs%linear_wave_drag)
then
1893 v_accel_bt(i,j) = v_accel_bt(i,j) + wt_accel(n) * &
1894 ((cor_v(i,j) + pfv(i,j)) - vbt(i,j)*rayleigh_v(i,j))
1896 v_accel_bt(i,j) = v_accel_bt(i,j) + wt_accel(n) * (cor_v(i,j) + pfv(i,j))
1899 if (use_bt_cont)
then
1901 do j=jsv-1,jev ;
do i=isv,iev
1902 vhbt(i,j) =
find_vhbt(vbt_trans(i,j), btcl_v(i,j), us) + vhbt0(i,j)
1906 do j=jsv-1,jev ;
do i=isv,iev
1907 vhbt(i,j) = datv(i,j)*vbt_trans(i,j) + vhbt0(i,j)
1910 if (cs%BT_OBC%apply_v_OBCs)
then
1912 do j=jsv-1,jev ;
do i=isv,iev ;
if (obc%segnum_v(i,j) /= obc_none)
then
1913 vbt(i,j) = vbt_prev(i,j); vhbt(i,j) = vhbt_prev(i,j)
1914 endif ;
enddo ;
enddo
1922 do j=js,je ;
do i=is-1,ie
1923 pfu_bt_sum(i,j) = pfu_bt_sum(i,j) + wt_accel2(n) * pfu(i,j)
1926 do j=js-1,je ;
do i=is,ie
1927 pfv_bt_sum(i,j) = pfv_bt_sum(i,j) + wt_accel2(n) * pfv(i,j)
1932 do j=js,je ;
do i=is-1,ie
1933 coru_bt_sum(i,j) = coru_bt_sum(i,j) + wt_accel2(n) * cor_u(i,j)
1936 do j=js-1,je ;
do i=is,ie
1937 corv_bt_sum(i,j) = corv_bt_sum(i,j) + wt_accel2(n) * cor_v(i,j)
1942 do j=js,je ;
do i=is-1,ie
1943 ubt_sum(i,j) = ubt_sum(i,j) + wt_trans(n) * ubt_trans(i,j)
1944 uhbt_sum(i,j) = uhbt_sum(i,j) + wt_trans(n) * uhbt(i,j)
1945 ubt_wtd(i,j) = ubt_wtd(i,j) + wt_vel(n) * ubt(i,j)
1948 do j=js-1,je ;
do i=is,ie
1949 vbt_sum(i,j) = vbt_sum(i,j) + wt_trans(n) * vbt_trans(i,j)
1950 vhbt_sum(i,j) = vhbt_sum(i,j) + wt_trans(n) * vhbt(i,j)
1951 vbt_wtd(i,j) = vbt_wtd(i,j) + wt_vel(n) * vbt(i,j)
1955 if (apply_obcs)
then
1956 if (cs%BT_OBC%apply_u_OBCs)
then
1958 do j=js,je ;
do i=is-1,ie
1959 if (obc%segnum_u(i,j) /= obc_none)
then
1960 ubt_sum(i,j) = ubt_sum_prev(i,j) ; uhbt_sum(i,j) = uhbt_sum_prev(i,j)
1961 ubt_wtd(i,j) = ubt_wtd_prev(i,j)
1966 if (cs%BT_OBC%apply_v_OBCs)
then
1968 do j=js-1,je ;
do i=is,ie
1969 if (obc%segnum_v(i,j) /= obc_none)
then
1970 vbt_sum(i,j) = vbt_sum_prev(i,j) ; vhbt_sum(i,j) = vhbt_sum_prev(i,j)
1971 vbt_wtd(i,j) = vbt_wtd_prev(i,j)
1977 ubt_trans, vbt_trans, eta, ubt_old, vbt_old, cs%BT_OBC, &
1978 g, ms, us, iev-ie, dtbt, bebt, use_bt_cont, datu, datv, btcl_u, btcl_v, &
1980 if (cs%BT_OBC%apply_u_OBCs)
then ;
do j=js,je ;
do i=is-1,ie
1981 if (obc%segnum_u(i,j) /= obc_none)
then
1982 ubt_sum(i,j) = ubt_sum(i,j) + wt_trans(n) * ubt_trans(i,j)
1983 uhbt_sum(i,j) = uhbt_sum(i,j) + wt_trans(n) * uhbt(i,j)
1984 ubt_wtd(i,j) = ubt_wtd(i,j) + wt_vel(n) * ubt(i,j)
1986 enddo ;
enddo ;
endif
1987 if (cs%BT_OBC%apply_v_OBCs)
then ;
do j=js-1,je ;
do i=is,ie
1988 if (obc%segnum_v(i,j) /= obc_none)
then
1989 vbt_sum(i,j) = vbt_sum(i,j) + wt_trans(n) * vbt_trans(i,j)
1990 vhbt_sum(i,j) = vhbt_sum(i,j) + wt_trans(n) * vhbt(i,j)
1991 vbt_wtd(i,j) = vbt_wtd(i,j) + wt_vel(n) * vbt(i,j)
1993 enddo ;
enddo ;
endif
1996 if (cs%debug_bt)
then
1997 call uvchksum(
"BT [uv]hbt just after OBC", uhbt, vhbt, cs%debug_BT_HI, haloshift=iev-ie, &
1998 scale=us%s_to_T*us%L_to_m**2*gv%H_to_m)
2002 do j=jsv,jev ;
do i=isv,iev
2003 eta(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * cs%IareaT(i,j)) * &
2004 ((uhbt(i-1,j) - uhbt(i,j)) + (vhbt(i,j-1) - vhbt(i,j)))
2005 eta_wtd(i,j) = eta_wtd(i,j) + eta(i,j) * wt_eta(n)
2009 if (do_hifreq_output)
then
2010 time_step_end = time_bt_start +
real_to_time(n*us%T_to_s*dtbt)
2012 if (cs%id_ubt_hifreq > 0)
call post_data(cs%id_ubt_hifreq, ubt(isdb:iedb,jsd:jed), cs%diag)
2013 if (cs%id_vbt_hifreq > 0)
call post_data(cs%id_vbt_hifreq, vbt(isd:ied,jsdb:jedb), cs%diag)
2014 if (cs%id_eta_hifreq > 0)
call post_data(cs%id_eta_hifreq, eta(isd:ied,jsd:jed), cs%diag)
2015 if (cs%id_uhbt_hifreq > 0)
call post_data(cs%id_uhbt_hifreq, uhbt(isdb:iedb,jsd:jed), cs%diag)
2016 if (cs%id_vhbt_hifreq > 0)
call post_data(cs%id_vhbt_hifreq, vhbt(isd:ied,jsdb:jedb), cs%diag)
2017 if (cs%id_eta_pred_hifreq > 0)
call post_data(cs%id_eta_pred_hifreq, eta_pf_bt(isd:ied,jsd:jed), cs%diag)
2020 if (cs%debug_bt)
then
2021 write(mesg,
'("BT step ",I4)') n
2022 call uvchksum(trim(mesg)//
" [uv]bt", ubt, vbt, cs%debug_BT_HI, haloshift=iev-ie, &
2023 scale=us%L_T_to_m_s)
2024 call hchksum(eta, trim(mesg)//
" eta", cs%debug_BT_HI, haloshift=iev-ie, scale=gv%H_to_m)
2032 if (do_hifreq_output)
call enable_averaging(time_int_in, time_end_in, cs%diag)
2034 i_sum_wt_vel = 1.0 / sum_wt_vel ; i_sum_wt_eta = 1.0 / sum_wt_eta
2035 i_sum_wt_accel = 1.0 / sum_wt_accel ; i_sum_wt_trans = 1.0 / sum_wt_trans
2037 if (find_etaav)
then ;
do j=js,je ;
do i=is,ie
2038 etaav(i,j) = eta_sum(i,j) * i_sum_wt_accel
2039 enddo ;
enddo ;
endif
2040 do j=js-1,je+1 ;
do i=is-1,ie+1 ; e_anom(i,j) = 0.0 ;
enddo ;
enddo
2041 if (interp_eta_pf)
then
2042 do j=js,je ;
do i=is,ie
2043 e_anom(i,j) = dgeo_de * (0.5 * (eta(i,j) + eta_in(i,j)) - &
2044 (eta_pf_1(i,j) + 0.5*d_eta_pf(i,j)))
2047 do j=js,je ;
do i=is,ie
2048 e_anom(i,j) = dgeo_de * (0.5 * (eta(i,j) + eta_in(i,j)) - eta_pf(i,j))
2051 if (apply_obcs)
then
2053 if (cs%BT_OBC%apply_u_OBCs)
then
2055 do j=js,je ;
do i=is-1,ie
2056 if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e)
then
2057 e_anom(i+1,j) = e_anom(i,j)
2058 elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w)
then
2059 e_anom(i,j) = e_anom(i+1,j)
2064 if (cs%BT_OBC%apply_v_OBCs)
then
2066 do j=js-1,je ;
do i=is,ie
2068 e_anom(i,j+1) = e_anom(i,j)
2069 elseif (obc%segment(obc%segnum_v(i,j))%direction ==
obc_direction_s)
then
2070 e_anom(i,j) = e_anom(i,j+1)
2077 do j=js,je ;
do i=is,ie
2078 eta_out(i,j) = eta_wtd(i,j) * i_sum_wt_eta
2083 if (g%nonblocking_updates)
then
2086 if (find_etaav)
call do_group_pass(cs%pass_etaav, g%Domain)
2087 call do_group_pass(cs%pass_e_anom, g%Domain)
2092 do j=js,je ;
do i=is-1,ie
2093 cs%ubtav(i,j) = ubt_sum(i,j) * i_sum_wt_trans
2094 uhbtav(i,j) = uhbt_sum(i,j) * i_sum_wt_trans
2097 ubt_wtd(i,j) = ubt_wtd(i,j) * i_sum_wt_vel
2100 do j=js-1,je ;
do i=is,ie
2101 cs%vbtav(i,j) = vbt_sum(i,j) * i_sum_wt_trans
2102 vhbtav(i,j) = vhbt_sum(i,j) * i_sum_wt_trans
2105 vbt_wtd(i,j) = vbt_wtd(i,j) * i_sum_wt_vel
2110 if (g%nonblocking_updates)
then
2115 call do_group_pass(cs%pass_ubta_uhbta, g%Domain)
2123 do j=js,je ;
do i=is-1,ie
2124 accel_layer_u(i,j,k) = (u_accel_bt(i,j) - &
2125 ((pbce(i+1,j,k) - gtot_w(i+1,j)) * e_anom(i+1,j) - &
2126 (pbce(i,j,k) - gtot_e(i,j)) * e_anom(i,j)) * cs%IdxCu(i,j) )
2127 if (abs(accel_layer_u(i,j,k)) < accel_underflow) accel_layer_u(i,j,k) = 0.0
2129 do j=js-1,je ;
do i=is,ie
2130 accel_layer_v(i,j,k) = (v_accel_bt(i,j) - &
2131 ((pbce(i,j+1,k) - gtot_s(i,j+1)) * e_anom(i,j+1) - &
2132 (pbce(i,j,k) - gtot_n(i,j)) * e_anom(i,j)) * cs%IdyCv(i,j) )
2133 if (abs(accel_layer_v(i,j,k)) < accel_underflow) accel_layer_v(i,j,k) = 0.0
2137 if (apply_obcs)
then
2140 if (cs%BT_OBC%apply_u_OBCs)
then ;
do j=js,je ;
do i=is-1,ie
2141 if (obc%segnum_u(i,j) /= obc_none)
then
2142 u_accel_bt(i,j) = (ubt_wtd(i,j) - ubt_first(i,j)) / dt
2143 do k=1,nz ; accel_layer_u(i,j,k) = u_accel_bt(i,j) ;
enddo
2145 enddo ;
enddo ;
endif
2146 if (cs%BT_OBC%apply_v_OBCs)
then ;
do j=js-1,je ;
do i=is,ie
2147 if (obc%segnum_v(i,j) /= obc_none)
then
2148 v_accel_bt(i,j) = (vbt_wtd(i,j) - vbt_first(i,j)) / dt
2149 do k=1,nz ; accel_layer_v(i,j,k) = v_accel_bt(i,j) ;
enddo
2151 enddo ;
enddo ;
endif
2157 if (query_averaging_enabled(cs%diag))
then
2159 do j=js,je ;
do i=is-1,ie ; cs%ubt_IC(i,j) = ubt_wtd(i,j) ;
enddo ;
enddo
2160 do j=js-1,je ;
do i=is,ie ; cs%vbt_IC(i,j) = vbt_wtd(i,j) ;
enddo ;
enddo
2161 if (use_bt_cont)
then
2162 do j=js,je ;
do i=is-1,ie
2163 cs%uhbt_IC(i,j) =
find_uhbt(ubt_wtd(i,j), btcl_u(i,j), us) + uhbt0(i,j)
2165 do j=js-1,je ;
do i=is,ie
2166 cs%vhbt_IC(i,j) =
find_vhbt(vbt_wtd(i,j), btcl_v(i,j), us) + vhbt0(i,j)
2169 do j=js,je ;
do i=is-1,ie
2170 cs%uhbt_IC(i,j) = ubt_wtd(i,j) * datu(i,j) + uhbt0(i,j)
2172 do j=js-1,je ;
do i=is,ie
2173 cs%vhbt_IC(i,j) = vbt_wtd(i,j) * datv(i,j) + vhbt0(i,j)
2178 if (cs%id_PFu_bt > 0)
then
2179 do j=js,je ;
do i=is-1,ie
2180 pfu_bt_sum(i,j) = pfu_bt_sum(i,j) * i_sum_wt_accel
2182 call post_data(cs%id_PFu_bt, pfu_bt_sum(isdb:iedb,jsd:jed), cs%diag)
2184 if (cs%id_PFv_bt > 0)
then
2185 do j=js-1,je ;
do i=is,ie
2186 pfv_bt_sum(i,j) = pfv_bt_sum(i,j) * i_sum_wt_accel
2188 call post_data(cs%id_PFv_bt, pfv_bt_sum(isd:ied,jsdb:jedb), cs%diag)
2190 if (cs%id_Coru_bt > 0)
then
2191 do j=js,je ;
do i=is-1,ie
2192 coru_bt_sum(i,j) = coru_bt_sum(i,j) * i_sum_wt_accel
2194 call post_data(cs%id_Coru_bt, coru_bt_sum(isdb:iedb,jsd:jed), cs%diag)
2196 if (cs%id_Corv_bt > 0)
then
2197 do j=js-1,je ;
do i=is,ie
2198 corv_bt_sum(i,j) = corv_bt_sum(i,j) * i_sum_wt_accel
2200 call post_data(cs%id_Corv_bt, corv_bt_sum(isd:ied,jsdb:jedb), cs%diag)
2202 if (cs%id_ubtforce > 0)
call post_data(cs%id_ubtforce, bt_force_u(isdb:iedb,jsd:jed), cs%diag)
2203 if (cs%id_vbtforce > 0)
call post_data(cs%id_vbtforce, bt_force_v(isd:ied,jsdb:jedb), cs%diag)
2204 if (cs%id_uaccel > 0)
call post_data(cs%id_uaccel, u_accel_bt(isdb:iedb,jsd:jed), cs%diag)
2205 if (cs%id_vaccel > 0)
call post_data(cs%id_vaccel, v_accel_bt(isd:ied,jsdb:jedb), cs%diag)
2207 if (cs%id_eta_cor > 0)
call post_data(cs%id_eta_cor, cs%eta_cor, cs%diag)
2208 if (cs%id_eta_bt > 0)
call post_data(cs%id_eta_bt, eta_out, cs%diag)
2209 if (cs%id_gtotn > 0)
call post_data(cs%id_gtotn, gtot_n(isd:ied,jsd:jed), cs%diag)
2210 if (cs%id_gtots > 0)
call post_data(cs%id_gtots, gtot_s(isd:ied,jsd:jed), cs%diag)
2211 if (cs%id_gtote > 0)
call post_data(cs%id_gtote, gtot_e(isd:ied,jsd:jed), cs%diag)
2212 if (cs%id_gtotw > 0)
call post_data(cs%id_gtotw, gtot_w(isd:ied,jsd:jed), cs%diag)
2213 if (cs%id_ubt > 0)
call post_data(cs%id_ubt, ubt_wtd(isdb:iedb,jsd:jed), cs%diag)
2214 if (cs%id_vbt > 0)
call post_data(cs%id_vbt, vbt_wtd(isd:ied,jsdb:jedb), cs%diag)
2215 if (cs%id_ubtav > 0)
call post_data(cs%id_ubtav, cs%ubtav, cs%diag)
2216 if (cs%id_vbtav > 0)
call post_data(cs%id_vbtav, cs%vbtav, cs%diag)
2217 if (cs%id_visc_rem_u > 0)
call post_data(cs%id_visc_rem_u, visc_rem_u, cs%diag)
2218 if (cs%id_visc_rem_v > 0)
call post_data(cs%id_visc_rem_v, visc_rem_v, cs%diag)
2220 if (cs%id_frhatu > 0)
call post_data(cs%id_frhatu, cs%frhatu, cs%diag)
2221 if (cs%id_uhbt > 0)
call post_data(cs%id_uhbt, uhbtav, cs%diag)
2222 if (cs%id_frhatv > 0)
call post_data(cs%id_frhatv, cs%frhatv, cs%diag)
2223 if (cs%id_vhbt > 0)
call post_data(cs%id_vhbt, vhbtav, cs%diag)
2224 if (cs%id_uhbt0 > 0)
call post_data(cs%id_uhbt0, uhbt0(isdb:iedb,jsd:jed), cs%diag)
2225 if (cs%id_vhbt0 > 0)
call post_data(cs%id_vhbt0, vhbt0(isd:ied,jsdb:jedb), cs%diag)
2227 if (cs%id_frhatu1 > 0)
call post_data(cs%id_frhatu1, cs%frhatu1, cs%diag)
2228 if (cs%id_frhatv1 > 0)
call post_data(cs%id_frhatv1, cs%frhatv1, cs%diag)
2230 if (use_bt_cont)
then
2231 if (cs%id_BTC_FA_u_EE > 0)
call post_data(cs%id_BTC_FA_u_EE, bt_cont%FA_u_EE, cs%diag)
2232 if (cs%id_BTC_FA_u_E0 > 0)
call post_data(cs%id_BTC_FA_u_E0, bt_cont%FA_u_E0, cs%diag)
2233 if (cs%id_BTC_FA_u_W0 > 0)
call post_data(cs%id_BTC_FA_u_W0, bt_cont%FA_u_W0, cs%diag)
2234 if (cs%id_BTC_FA_u_WW > 0)
call post_data(cs%id_BTC_FA_u_WW, bt_cont%FA_u_WW, cs%diag)
2235 if (cs%id_BTC_uBT_EE > 0)
call post_data(cs%id_BTC_uBT_EE, bt_cont%uBT_EE, cs%diag)
2236 if (cs%id_BTC_uBT_WW > 0)
call post_data(cs%id_BTC_uBT_WW, bt_cont%uBT_WW, cs%diag)
2237 if (cs%id_BTC_FA_v_NN > 0)
call post_data(cs%id_BTC_FA_v_NN, bt_cont%FA_v_NN, cs%diag)
2238 if (cs%id_BTC_FA_v_N0 > 0)
call post_data(cs%id_BTC_FA_v_N0, bt_cont%FA_v_N0, cs%diag)
2239 if (cs%id_BTC_FA_v_S0 > 0)
call post_data(cs%id_BTC_FA_v_S0, bt_cont%FA_v_S0, cs%diag)
2240 if (cs%id_BTC_FA_v_SS > 0)
call post_data(cs%id_BTC_FA_v_SS, bt_cont%FA_v_SS, cs%diag)
2241 if (cs%id_BTC_vBT_NN > 0)
call post_data(cs%id_BTC_vBT_NN, bt_cont%vBT_NN, cs%diag)
2242 if (cs%id_BTC_vBT_SS > 0)
call post_data(cs%id_BTC_vBT_SS, bt_cont%vBT_SS, cs%diag)
2245 if (cs%id_frhatu1 > 0) cs%frhatu1(:,:,:) = cs%frhatu(:,:,:)
2246 if (cs%id_frhatv1 > 0) cs%frhatv1(:,:,:) = cs%frhatv(:,:,:)
2249 if (g%nonblocking_updates)
then
2258 subroutine set_dtbt(G, GV, US, CS, eta, pbce, BT_cont, gtot_est, SSH_add)
2263 real,
dimension(SZI_(G),SZJ_(G)),
optional,
intent(in) :: eta
2265 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
optional,
intent(in) :: pbce
2271 real,
optional,
intent(in) :: gtot_est
2273 real,
optional,
intent(in) :: ssh_add
2278 real,
dimension(SZI_(G),SZJ_(G)) :: &
2285 real,
dimension(SZIBS_(G),SZJ_(G)) :: &
2288 real,
dimension(SZI_(G),SZJBS_(G)) :: &
2305 logical :: use_bt_cont
2308 character(len=200) :: mesg
2309 integer :: i, j, k, is, ie, js, je, nz
2311 if (.not.
associated(cs))
call mom_error(fatal, &
2312 "set_dtbt: Module MOM_barotropic must be initialized before it is used.")
2313 if (.not.cs%split)
return
2314 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
2315 ms%isdw = g%isd ; ms%iedw = g%ied ; ms%jsdw = g%jsd ; ms%jedw = g%jed
2317 if (.not.(
present(pbce) .or.
present(gtot_est)))
call mom_error(fatal, &
2318 "set_dtbt: Either pbce or gtot_est must be present.")
2320 add_ssh = 0.0 ;
if (
present(ssh_add)) add_ssh = ssh_add
2322 use_bt_cont = .false.
2323 if (
present(bt_cont)) use_bt_cont = (
associated(bt_cont))
2325 if (use_bt_cont)
then
2327 elseif (cs%Nonlinear_continuity .and.
present(eta))
then
2330 call find_face_areas(datu, datv, g, gv, us, cs, ms, halo=0, add_max=add_ssh)
2335 dgeo_de = 1.0 + max(0.0, det_de + cs%G_extra)
2336 if (
present(pbce))
then
2337 do j=js,je ;
do i=is,ie
2338 gtot_e(i,j) = 0.0 ; gtot_w(i,j) = 0.0
2339 gtot_n(i,j) = 0.0 ; gtot_s(i,j) = 0.0
2341 do k=1,nz ;
do j=js,je ;
do i=is,ie
2342 gtot_e(i,j) = gtot_e(i,j) + pbce(i,j,k) * cs%frhatu(i,j,k)
2343 gtot_w(i,j) = gtot_w(i,j) + pbce(i,j,k) * cs%frhatu(i-1,j,k)
2344 gtot_n(i,j) = gtot_n(i,j) + pbce(i,j,k) * cs%frhatv(i,j,k)
2345 gtot_s(i,j) = gtot_s(i,j) + pbce(i,j,k) * cs%frhatv(i,j-1,k)
2346 enddo ;
enddo ;
enddo
2348 do j=js,je ;
do i=is,ie
2349 gtot_e(i,j) = gtot_est * gv%H_to_Z ; gtot_w(i,j) = gtot_est * gv%H_to_Z
2350 gtot_n(i,j) = gtot_est * gv%H_to_Z ; gtot_s(i,j) = gtot_est * gv%H_to_Z
2354 min_max_dt2 = 1.0e38*us%s_to_T**2
2355 do j=js,je ;
do i=is,ie
2358 idt_max2 = 0.5 * (1.0 + 2.0*cs%bebt) * (g%IareaT(i,j) * &
2359 ((gtot_e(i,j)*datu(i,j)*g%IdxCu(i,j) + gtot_w(i,j)*datu(i-1,j)*g%IdxCu(i-1,j)) + &
2360 (gtot_n(i,j)*datv(i,j)*g%IdyCv(i,j) + gtot_s(i,j)*datv(i,j-1)*g%IdyCv(i,j-1))) + &
2361 ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
2362 (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2)))
2363 if (idt_max2 * min_max_dt2 > 1.0) min_max_dt2 = 1.0 / idt_max2
2365 dtbt_max = sqrt(min_max_dt2 / dgeo_de)
2367 call min_across_pes(dtbt_max)
2370 cs%dtbt = cs%dtbt_fraction * dtbt_max
2371 cs%dtbt_max = dtbt_max
2378 eta, ubt_old, vbt_old, BT_OBC, &
2379 G, MS, US, halo, dtbt, bebt, use_BT_cont, Datu, Datv, &
2380 BTCL_u, BTCL_v, uhbt0, vhbt0)
2385 real,
dimension(SZIBW_(MS),SZJW_(MS)),
intent(inout) :: ubt
2386 real,
dimension(SZIBW_(MS),SZJW_(MS)),
intent(inout) :: uhbt
2388 real,
dimension(SZIBW_(MS),SZJW_(MS)),
intent(inout) :: ubt_trans
2390 real,
dimension(SZIW_(MS),SZJBW_(MS)),
intent(inout) :: vbt
2392 real,
dimension(SZIW_(MS),SZJBW_(MS)),
intent(inout) :: vhbt
2394 real,
dimension(SZIW_(MS),SZJBW_(MS)),
intent(inout) :: vbt_trans
2396 real,
dimension(SZIW_(MS),SZJW_(MS)),
intent(in) :: eta
2398 real,
dimension(SZIBW_(MS),SZJW_(MS)),
intent(in) :: ubt_old
2400 real,
dimension(SZIW_(MS),SZJBW_(MS)),
intent(in) :: vbt_old
2406 integer,
intent(in) :: halo
2407 real,
intent(in) :: dtbt
2408 real,
intent(in) :: bebt
2410 logical,
intent(in) :: use_BT_cont
2412 real,
dimension(SZIBW_(MS),SZJW_(MS)),
intent(in) :: Datu
2414 real,
dimension(SZIW_(MS),SZJBW_(MS)),
intent(in) :: Datv
2422 real,
dimension(SZIBW_(MS),SZJW_(MS)),
intent(in) :: uhbt0
2426 real,
dimension(SZIW_(MS),SZJBW_(MS)),
intent(in) :: vhbt0
2441 real :: cff, Cx, Cy, tau
2442 real :: dhdt, dhdx, dhdy
2443 integer :: i, j, is, ie, js, je
2444 real,
dimension(SZIB_(G),SZJB_(G)) :: grad
2445 real,
parameter :: eps = 1.0e-20
2446 real :: rx_max, ry_max
2447 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
2448 rx_max = obc%rx_max ; ry_max = obc%rx_max
2450 if (bt_obc%apply_u_OBCs)
then
2451 do j=js,je ;
do i=is-1,ie ;
if (obc%segnum_u(i,j) /= obc_none)
then
2452 if (obc%segment(obc%segnum_u(i,j))%specified)
then
2453 uhbt(i,j) = bt_obc%uhbt(i,j)
2454 ubt(i,j) = bt_obc%ubt_outer(i,j)
2455 vel_trans = ubt(i,j)
2456 elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e)
then
2457 if (obc%segment(obc%segnum_u(i,j))%Flather)
then
2458 cfl = dtbt * bt_obc%Cg_u(i,j) * g%IdxCu(i,j)
2459 u_inlet = cfl*ubt_old(i-1,j) + (1.0-cfl)*ubt_old(i,j)
2460 h_in = eta(i,j) + (0.5-cfl)*(eta(i,j)-eta(i-1,j))
2461 h_u = bt_obc%H_u(i,j)
2463 ubt(i,j) = 0.5*((u_inlet + bt_obc%ubt_outer(i,j)) + &
2464 (bt_obc%Cg_u(i,j)/h_u) * (h_in-bt_obc%eta_outer_u(i,j)))
2465 vel_trans = (1.0-bebt)*vel_prev + bebt*ubt(i,j)
2466 elseif (obc%segment(obc%segnum_u(i,j))%gradient)
then
2467 ubt(i,j) = ubt(i-1,j)
2468 vel_trans = ubt(i,j)
2470 elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w)
then
2471 if (obc%segment(obc%segnum_u(i,j))%Flather)
then
2472 cfl = dtbt * bt_obc%Cg_u(i,j) * g%IdxCu(i,j)
2473 u_inlet = cfl*ubt_old(i+1,j) + (1.0-cfl)*ubt_old(i,j)
2474 h_in = eta(i+1,j) + (0.5-cfl)*(eta(i+1,j)-eta(i+2,j))
2476 h_u = bt_obc%H_u(i,j)
2478 ubt(i,j) = 0.5*((u_inlet + bt_obc%ubt_outer(i,j)) + &
2479 (bt_obc%Cg_u(i,j)/h_u) * (bt_obc%eta_outer_u(i,j)-h_in))
2481 vel_trans = (1.0-bebt)*vel_prev + bebt*ubt(i,j)
2482 elseif (obc%segment(obc%segnum_u(i,j))%gradient)
then
2483 ubt(i,j) = ubt(i+1,j)
2484 vel_trans = ubt(i,j)
2488 if (.not. obc%segment(obc%segnum_u(i,j))%specified)
then
2489 if (use_bt_cont)
then
2490 uhbt(i,j) =
find_uhbt(vel_trans, btcl_u(i,j), us) + uhbt0(i,j)
2492 uhbt(i,j) = datu(i,j)*vel_trans + uhbt0(i,j)
2496 ubt_trans(i,j) = vel_trans
2497 endif ;
enddo ;
enddo
2500 if (bt_obc%apply_v_OBCs)
then
2501 do j=js-1,je ;
do i=is,ie ;
if (obc%segnum_v(i,j) /= obc_none)
then
2502 if (obc%segment(obc%segnum_v(i,j))%specified)
then
2503 vhbt(i,j) = bt_obc%vhbt(i,j)
2504 vbt(i,j) = bt_obc%vbt_outer(i,j)
2505 vel_trans = vbt(i,j)
2506 elseif (obc%segment(obc%segnum_v(i,j))%direction ==
obc_direction_n)
then
2507 if (obc%segment(obc%segnum_v(i,j))%Flather)
then
2508 cfl = dtbt * bt_obc%Cg_v(i,j) * g%IdyCv(i,j)
2509 v_inlet = cfl*vbt_old(i,j-1) + (1.0-cfl)*vbt_old(i,j)
2510 h_in = eta(i,j) + (0.5-cfl)*(eta(i,j)-eta(i,j-1))
2512 h_v = bt_obc%H_v(i,j)
2514 vbt(i,j) = 0.5*((v_inlet + bt_obc%vbt_outer(i,j)) + &
2515 (bt_obc%Cg_v(i,j)/h_v) * (h_in-bt_obc%eta_outer_v(i,j)))
2517 vel_trans = (1.0-bebt)*vel_prev + bebt*vbt(i,j)
2518 elseif (obc%segment(obc%segnum_v(i,j))%gradient)
then
2519 vbt(i,j) = vbt(i,j-1)
2520 vel_trans = vbt(i,j)
2522 elseif (obc%segment(obc%segnum_v(i,j))%direction ==
obc_direction_s)
then
2523 if (obc%segment(obc%segnum_v(i,j))%Flather)
then
2524 cfl = dtbt * bt_obc%Cg_v(i,j) * g%IdyCv(i,j)
2525 v_inlet = cfl*vbt_old(i,j+1) + (1.0-cfl)*vbt_old(i,j)
2526 h_in = eta(i,j+1) + (0.5-cfl)*(eta(i,j+1)-eta(i,j+2))
2528 h_v = bt_obc%H_v(i,j)
2530 vbt(i,j) = 0.5*((v_inlet + bt_obc%vbt_outer(i,j)) + &
2531 (bt_obc%Cg_v(i,j)/h_v) * (bt_obc%eta_outer_v(i,j)-h_in))
2533 vel_trans = (1.0-bebt)*vel_prev + bebt*vbt(i,j)
2534 elseif (obc%segment(obc%segnum_v(i,j))%gradient)
then
2535 vbt(i,j) = vbt(i,j+1)
2536 vel_trans = vbt(i,j)
2540 if (.not. obc%segment(obc%segnum_v(i,j))%specified)
then
2541 if (use_bt_cont)
then
2542 vhbt(i,j) =
find_vhbt(vel_trans, btcl_v(i,j), us) + vhbt0(i,j)
2544 vhbt(i,j) = vel_trans*datv(i,j) + vhbt0(i,j)
2548 vbt_trans(i,j) = vel_trans
2549 endif ;
enddo ;
enddo
2556 subroutine set_up_bt_obc(OBC, eta, BT_OBC, BT_Domain, G, GV, US, MS, halo, use_BT_cont, Datu, Datv, BTCL_u, BTCL_v)
2560 real,
dimension(SZIW_(MS),SZJW_(MS)),
intent(in) :: eta
2569 integer,
intent(in) :: halo
2570 logical,
intent(in) :: use_BT_cont
2572 real,
dimension(SZIBW_(MS),SZJW_(MS)),
intent(in) :: Datu
2574 real,
dimension(SZIW_(MS),SZJBW_(MS)),
intent(in) :: Datv
2584 integer :: i, j, k, is, ie, js, je, n, nz, Isq, Ieq, Jsq, Jeq
2585 integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
2586 integer :: isdw, iedw, jsdw, jedw
2591 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
2592 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
2593 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
2594 isdw = ms%isdw ; iedw = ms%iedw ; jsdw = ms%jsdw ; jedw = ms%jedw
2597 if ((isdw < isd) .or. (jsdw < jsd))
then
2598 call mom_error(fatal,
"set_up_BT_OBC: Open boundary conditions are not "//&
2599 "yet fully implemented with wide barotropic halos.")
2602 if (.not. bt_obc%is_alloced)
then
2603 allocate(bt_obc%Cg_u(isdw-1:iedw,jsdw:jedw)) ; bt_obc%Cg_u(:,:) = 0.0
2604 allocate(bt_obc%H_u(isdw-1:iedw,jsdw:jedw)) ; bt_obc%H_u(:,:) = 0.0
2605 allocate(bt_obc%uhbt(isdw-1:iedw,jsdw:jedw)) ; bt_obc%uhbt(:,:) = 0.0
2606 allocate(bt_obc%ubt_outer(isdw-1:iedw,jsdw:jedw)) ; bt_obc%ubt_outer(:,:) = 0.0
2607 allocate(bt_obc%eta_outer_u(isdw-1:iedw,jsdw:jedw)) ; bt_obc%eta_outer_u(:,:) = 0.0
2609 allocate(bt_obc%Cg_v(isdw:iedw,jsdw-1:jedw)) ; bt_obc%Cg_v(:,:) = 0.0
2610 allocate(bt_obc%H_v(isdw:iedw,jsdw-1:jedw)) ; bt_obc%H_v(:,:) = 0.0
2611 allocate(bt_obc%vhbt(isdw:iedw,jsdw-1:jedw)) ; bt_obc%vhbt(:,:) = 0.0
2612 allocate(bt_obc%vbt_outer(isdw:iedw,jsdw-1:jedw)) ; bt_obc%vbt_outer(:,:) = 0.0
2613 allocate(bt_obc%eta_outer_v(isdw:iedw,jsdw-1:jedw)) ; bt_obc%eta_outer_v(:,:)=0.0
2614 bt_obc%is_alloced = .true.
2615 call create_group_pass(bt_obc%pass_uv, bt_obc%ubt_outer, bt_obc%vbt_outer, bt_domain)
2617 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)
2618 call create_group_pass(bt_obc%pass_h, bt_obc%H_u, bt_obc%H_v, bt_domain,to_all+scalar_pair)
2619 call create_group_pass(bt_obc%pass_cg, bt_obc%Cg_u, bt_obc%Cg_v, bt_domain,to_all+scalar_pair)
2622 if (bt_obc%apply_u_OBCs)
then
2623 if (obc%specified_u_BCs_exist_globally)
then
2624 do n = 1, obc%number_of_segments
2625 segment => obc%segment(n)
2626 if (segment%is_E_or_W .and. segment%specified)
then
2627 do j=segment%HI%jsd,segment%HI%jed ;
do i=segment%HI%IsdB,segment%HI%IedB
2628 bt_obc%uhbt(i,j) = 0.
2630 do k=1,nz ;
do j=segment%HI%jsd,segment%HI%jed ;
do i=segment%HI%IsdB,segment%HI%IedB
2631 bt_obc%uhbt(i,j) = bt_obc%uhbt(i,j) + segment%normal_trans(i,j,k)
2632 enddo ;
enddo ;
enddo
2636 do j=js,je ;
do i=is-1,ie ;
if (obc%segnum_u(i,j) /= obc_none)
then
2638 if (obc%segment(obc%segnum_u(i,j))%specified)
then
2639 if (use_bt_cont)
then
2640 bt_obc%ubt_outer(i,j) =
uhbt_to_ubt(bt_obc%uhbt(i,j), btcl_u(i,j), us)
2642 if (datu(i,j) > 0.0) bt_obc%ubt_outer(i,j) = bt_obc%uhbt(i,j) / datu(i,j)
2645 if (gv%Boussinesq)
then
2646 if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e)
then
2647 bt_obc%H_u(i,j) = g%bathyT(i,j)*gv%Z_to_H + eta(i,j)
2648 elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w)
then
2649 bt_obc%H_u(i,j) = g%bathyT(i+1,j)*gv%Z_to_H + eta(i+1,j)
2652 if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e)
then
2653 bt_obc%H_u(i,j) = eta(i,j)
2654 elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w)
then
2655 bt_obc%H_u(i,j) = eta(i+1,j)
2658 bt_obc%Cg_u(i,j) = sqrt(gv%g_prime(1) * gv%H_to_Z*bt_obc%H_u(i,j))
2660 endif ;
enddo ;
enddo
2661 if (obc%Flather_u_BCs_exist_globally)
then
2662 do n = 1, obc%number_of_segments
2663 segment => obc%segment(n)
2664 if (segment%is_E_or_W .and. segment%Flather)
then
2665 do j=segment%HI%jsd,segment%HI%jed ;
do i=segment%HI%IsdB,segment%HI%IedB
2666 bt_obc%ubt_outer(i,j) = segment%normal_vel_bt(i,j)
2667 bt_obc%eta_outer_u(i,j) = segment%eta(i,j)
2674 if (bt_obc%apply_v_OBCs)
then
2675 if (obc%specified_v_BCs_exist_globally)
then
2676 do n = 1, obc%number_of_segments
2677 segment => obc%segment(n)
2678 if (segment%is_N_or_S .and. segment%specified)
then
2679 do j=segment%HI%JsdB,segment%HI%JedB ;
do i=segment%HI%isd,segment%HI%ied
2680 bt_obc%vhbt(i,j) = 0.
2682 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB ;
do i=segment%HI%isd,segment%HI%ied
2683 bt_obc%vhbt(i,j) = bt_obc%vhbt(i,j) + segment%normal_trans(i,j,k)
2684 enddo ;
enddo ;
enddo
2688 do j=js-1,je ;
do i=is,ie ;
if (obc%segnum_v(i,j) /= obc_none)
then
2690 if (obc%segment(obc%segnum_v(i,j))%specified)
then
2691 if (use_bt_cont)
then
2692 bt_obc%vbt_outer(i,j) =
vhbt_to_vbt(bt_obc%vhbt(i,j), btcl_v(i,j), us)
2694 if (datv(i,j) > 0.0) bt_obc%vbt_outer(i,j) = bt_obc%vhbt(i,j) / datv(i,j)
2697 if (gv%Boussinesq)
then
2699 bt_obc%H_v(i,j) = g%bathyT(i,j)*gv%Z_to_H + eta(i,j)
2700 elseif (obc%segment(obc%segnum_v(i,j))%direction ==
obc_direction_s)
then
2701 bt_obc%H_v(i,j) = g%bathyT(i,j+1)*gv%Z_to_H + eta(i,j+1)
2705 bt_obc%H_v(i,j) = eta(i,j)
2706 elseif (obc%segment(obc%segnum_v(i,j))%direction ==
obc_direction_s)
then
2707 bt_obc%H_v(i,j) = eta(i,j+1)
2710 bt_obc%Cg_v(i,j) = sqrt(gv%g_prime(1) * gv%H_to_Z*bt_obc%H_v(i,j))
2712 endif ;
enddo ;
enddo
2713 if (obc%Flather_v_BCs_exist_globally)
then
2714 do n = 1, obc%number_of_segments
2715 segment => obc%segment(n)
2716 if (segment%is_N_or_S .and. segment%Flather)
then
2717 do j=segment%HI%JsdB,segment%HI%JedB ;
do i=segment%HI%isd,segment%HI%ied
2718 bt_obc%vbt_outer(i,j) = segment%normal_vel_bt(i,j)
2719 bt_obc%eta_outer_v(i,j) = segment%eta(i,j)
2726 call do_group_pass(bt_obc%pass_uv, bt_domain)
2727 call do_group_pass(bt_obc%pass_uhvh, bt_domain)
2728 call do_group_pass(bt_obc%pass_eta_outer, bt_domain)
2729 call do_group_pass(bt_obc%pass_h, bt_domain)
2730 call do_group_pass(bt_obc%pass_cg, bt_domain)
2740 if (bt_obc%is_alloced)
then
2741 deallocate(bt_obc%Cg_u)
2742 deallocate(bt_obc%H_u)
2743 deallocate(bt_obc%uhbt)
2744 deallocate(bt_obc%ubt_outer)
2745 deallocate(bt_obc%eta_outer_u)
2747 deallocate(bt_obc%Cg_v)
2748 deallocate(bt_obc%H_v)
2749 deallocate(bt_obc%vhbt)
2750 deallocate(bt_obc%vbt_outer)
2751 deallocate(bt_obc%eta_outer_v)
2752 bt_obc%is_alloced = .false.
2761 subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC)
2764 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
2768 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
2769 optional,
intent(in) :: h_u
2770 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
2771 optional,
intent(in) :: h_v
2772 logical,
optional,
intent(in) :: may_use_default
2781 real :: hatutot(szib_(g))
2782 real :: hatvtot(szi_(g))
2783 real :: ihatutot(szib_(g))
2784 real :: ihatvtot(szi_(g))
2792 real :: e_u(szib_(g),szk_(g)+1)
2793 real :: e_v(szi_(g),szk_(g)+1)
2794 real :: d_shallow_u(szi_(g))
2795 real :: d_shallow_v(szib_(g))
2799 logical :: use_default, test_dflt, apply_obcs
2800 integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz, i, j, k
2801 integer :: iss, ies, n
2805 if (.not.
associated(cs))
call mom_error(fatal, &
2806 "btcalc: Module MOM_barotropic must be initialized before it is used.")
2807 if (.not.cs%split)
return
2809 use_default = .false.
2810 test_dflt = .false. ;
if (
present(may_use_default)) test_dflt = may_use_default
2813 if (.not.((
present(h_u) .and.
present(h_v)) .or. &
2814 (cs%hvel_scheme ==
harmonic) .or. (cs%hvel_scheme ==
hybrid) .or.&
2815 (cs%hvel_scheme ==
arithmetic))) use_default = .true.
2817 if (.not.((
present(h_u) .and.
present(h_v)) .or. &
2818 (cs%hvel_scheme ==
harmonic) .or. (cs%hvel_scheme ==
hybrid) .or.&
2820 "btcalc: Inconsistent settings of optional arguments and hvel_scheme.")
2823 apply_obcs = .false.
2824 if (
present(obc))
then ;
if (
associated(obc))
then ;
if (obc%OBC_pe)
then
2827 apply_obcs = (obc%number_of_segments > 0)
2828 endif ;
endif ;
endif
2830 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
2831 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
2832 h_neglect = gv%H_subroundoff
2840 if (
present(h_u))
then
2841 do i=is-1,ie ; hatutot(i) = h_u(i,j,1) ;
enddo
2842 do k=2,nz ;
do i=is-1,ie
2843 hatutot(i) = hatutot(i) + h_u(i,j,k)
2845 do i=is-1,ie ; ihatutot(i) = g%mask2dCu(i,j) / (hatutot(i) + h_neglect) ;
enddo
2846 do k=1,nz ;
do i=is-1,ie
2847 cs%frhatu(i,j,k) = h_u(i,j,k) * ihatutot(i)
2852 cs%frhatu(i,j,1) = 0.5 * (h(i+1,j,1) + h(i,j,1))
2853 hatutot(i) = cs%frhatu(i,j,1)
2855 do k=2,nz ;
do i=is-1,ie
2856 cs%frhatu(i,j,k) = 0.5 * (h(i+1,j,k) + h(i,j,k))
2857 hatutot(i) = hatutot(i) + cs%frhatu(i,j,k)
2859 elseif (cs%hvel_scheme ==
hybrid .or. use_default)
then
2861 e_u(i,nz+1) = -0.5 * gv%Z_to_H * (g%bathyT(i+1,j) + g%bathyT(i,j))
2862 d_shallow_u(i) = -gv%Z_to_H * min(g%bathyT(i+1,j), g%bathyT(i,j))
2865 do k=nz,1,-1 ;
do i=is-1,ie
2866 e_u(i,k) = e_u(i,k+1) + 0.5 * (h(i+1,j,k) + h(i,j,k))
2867 h_arith = 0.5 * (h(i+1,j,k) + h(i,j,k))
2868 if (e_u(i,k+1) >= d_shallow_u(i))
then
2869 cs%frhatu(i,j,k) = h_arith
2871 h_harm = (h(i+1,j,k) * h(i,j,k)) / (h_arith + h_neglect)
2872 if (e_u(i,k) <= d_shallow_u(i))
then
2873 cs%frhatu(i,j,k) = h_harm
2875 wt_arith = (e_u(i,k) - d_shallow_u(i)) / (h_arith + h_neglect)
2876 cs%frhatu(i,j,k) = wt_arith*h_arith + (1.0-wt_arith)*h_harm
2879 hatutot(i) = hatutot(i) + cs%frhatu(i,j,k)
2881 elseif (cs%hvel_scheme ==
harmonic)
then
2883 cs%frhatu(i,j,1) = 2.0*(h(i+1,j,1) * h(i,j,1)) / &
2884 ((h(i+1,j,1) + h(i,j,1)) + h_neglect)
2885 hatutot(i) = cs%frhatu(i,j,1)
2887 do k=2,nz ;
do i=is-1,ie
2888 cs%frhatu(i,j,k) = 2.0*(h(i+1,j,k) * h(i,j,k)) / &
2889 ((h(i+1,j,k) + h(i,j,k)) + h_neglect)
2890 hatutot(i) = hatutot(i) + cs%frhatu(i,j,k)
2893 do i=is-1,ie ; ihatutot(i) = g%mask2dCu(i,j) / (hatutot(i) + h_neglect) ;
enddo
2894 do k=1,nz ;
do i=is-1,ie
2895 cs%frhatu(i,j,k) = cs%frhatu(i,j,k) * ihatutot(i)
2903 if (
present(h_v))
then
2904 do i=is,ie ; hatvtot(i) = h_v(i,j,1) ;
enddo
2905 do k=2,nz ;
do i=is,ie
2906 hatvtot(i) = hatvtot(i) + h_v(i,j,k)
2908 do i=is,ie ; ihatvtot(i) = g%mask2dCv(i,j) / (hatvtot(i) + h_neglect) ;
enddo
2909 do k=1,nz ;
do i=is,ie
2910 cs%frhatv(i,j,k) = h_v(i,j,k) * ihatvtot(i)
2915 cs%frhatv(i,j,1) = 0.5 * (h(i,j+1,1) + h(i,j,1))
2916 hatvtot(i) = cs%frhatv(i,j,1)
2918 do k=2,nz ;
do i=is,ie
2919 cs%frhatv(i,j,k) = 0.5 * (h(i,j+1,k) + h(i,j,k))
2920 hatvtot(i) = hatvtot(i) + cs%frhatv(i,j,k)
2922 elseif (cs%hvel_scheme ==
hybrid .or. use_default)
then
2924 e_v(i,nz+1) = -0.5 * gv%Z_to_H * (g%bathyT(i,j+1) + g%bathyT(i,j))
2925 d_shallow_v(i) = -gv%Z_to_H * min(g%bathyT(i,j+1), g%bathyT(i,j))
2928 do k=nz,1,-1 ;
do i=is,ie
2929 e_v(i,k) = e_v(i,k+1) + 0.5 * (h(i,j+1,k) + h(i,j,k))
2930 h_arith = 0.5 * (h(i,j+1,k) + h(i,j,k))
2931 if (e_v(i,k+1) >= d_shallow_v(i))
then
2932 cs%frhatv(i,j,k) = h_arith
2934 h_harm = (h(i,j+1,k) * h(i,j,k)) / (h_arith + h_neglect)
2935 if (e_v(i,k) <= d_shallow_v(i))
then
2936 cs%frhatv(i,j,k) = h_harm
2938 wt_arith = (e_v(i,k) - d_shallow_v(i)) / (h_arith + h_neglect)
2939 cs%frhatv(i,j,k) = wt_arith*h_arith + (1.0-wt_arith)*h_harm
2942 hatvtot(i) = hatvtot(i) + cs%frhatv(i,j,k)
2944 elseif (cs%hvel_scheme ==
harmonic)
then
2946 cs%frhatv(i,j,1) = 2.0*(h(i,j+1,1) * h(i,j,1)) / &
2947 ((h(i,j+1,1) + h(i,j,1)) + h_neglect)
2948 hatvtot(i) = cs%frhatv(i,j,1)
2950 do k=2,nz ;
do i=is,ie
2951 cs%frhatv(i,j,k) = 2.0*(h(i,j+1,k) * h(i,j,k)) / &
2952 ((h(i,j+1,k) + h(i,j,k)) + h_neglect)
2953 hatvtot(i) = hatvtot(i) + cs%frhatv(i,j,k)
2956 do i=is,ie ; ihatvtot(i) = g%mask2dCv(i,j) / (hatvtot(i) + h_neglect) ;
enddo
2957 do k=1,nz ;
do i=is,ie
2958 cs%frhatv(i,j,k) = cs%frhatv(i,j,k) * ihatvtot(i)
2963 if (apply_obcs)
then ;
do n=1,obc%number_of_segments
2964 if (.not. obc%segment(n)%on_pe) cycle
2966 j = obc%segment(n)%HI%JsdB
2967 if ((j >= js-1) .and. (j <= je))
then
2968 iss = max(is,obc%segment(n)%HI%isd) ; ies = min(ie,obc%segment(n)%HI%ied)
2969 do i=iss,ies ; hatvtot(i) = h(i,j,1) ;
enddo
2970 do k=2,nz ;
do i=iss,ies
2971 hatvtot(i) = hatvtot(i) + h(i,j,k)
2974 ihatvtot(i) = g%mask2dCv(i,j) / (hatvtot(i) + h_neglect)
2976 do k=1,nz ;
do i=iss,ies
2977 cs%frhatv(i,j,k) = h(i,j,k) * ihatvtot(i)
2981 j = obc%segment(n)%HI%JsdB
2982 if ((j >= js-1) .and. (j <= je))
then
2983 iss = max(is,obc%segment(n)%HI%isd) ; ies = min(ie,obc%segment(n)%HI%ied)
2984 do i=iss,ies ; hatvtot(i) = h(i,j+1,1) ;
enddo
2985 do k=2,nz ;
do i=iss,ies
2986 hatvtot(i) = hatvtot(i) + h(i,j+1,k)
2989 ihatvtot(i) = g%mask2dCv(i,j) / (hatvtot(i) + h_neglect)
2991 do k=1,nz ;
do i=iss,ies
2992 cs%frhatv(i,j,k) = h(i,j+1,k) * ihatvtot(i)
2995 elseif (obc%segment(n)%direction == obc_direction_e)
then
2996 i = obc%segment(n)%HI%IsdB
2997 if ((i >= is-1) .and. (i <= ie))
then
2998 do j = max(js,obc%segment(n)%HI%jsd), min(je,obc%segment(n)%HI%jed)
3000 do k=2,nz ; htot = htot + h(i,j,k) ;
enddo
3001 ihtot = g%mask2dCu(i,j) / (htot + h_neglect)
3002 do k=1,nz ; cs%frhatu(i,j,k) = h(i,j,k) * ihtot ;
enddo
3005 elseif (obc%segment(n)%direction == obc_direction_w)
then
3006 i = obc%segment(n)%HI%IsdB
3007 if ((i >= is-1) .and. (i <= ie))
then
3008 do j = max(js,obc%segment(n)%HI%jsd), min(je,obc%segment(n)%HI%jed)
3010 do k=2,nz ; htot = htot + h(i+1,j,k) ;
enddo
3011 ihtot = g%mask2dCu(i,j) / (htot + h_neglect)
3012 do k=1,nz ; cs%frhatu(i,j,k) = h(i+1,j,k) * ihtot ;
enddo
3016 call mom_error(fatal,
"btcalc encountered and OBC segment of indeterminate direction.")
3021 call uvchksum(
"btcalc frhat[uv]", cs%frhatu, cs%frhatv, g%HI, 0, .true., .true.)
3022 if (
present(h_u) .and.
present(h_v)) &
3023 call uvchksum(
"btcalc h_[uv]", h_u, h_v, g%HI, 0, .true., .true., scale=gv%H_to_m)
3024 call hchksum(h,
"btcalc h",g%HI, haloshift=1, scale=gv%H_to_m)
3030 function find_uhbt(u, BTC, US)
result(uhbt)
3031 real,
intent(in) :: u
3041 elseif (u < btc%uBT_EE)
then
3042 uhbt = (u - btc%uBT_EE) * btc%FA_u_EE + btc%uh_EE
3043 elseif (u < 0.0)
then
3044 uhbt = u * (btc%FA_u_E0 + btc%uh_crvE * u**2)
3045 elseif (u <= btc%uBT_WW)
then
3046 uhbt = u * (btc%FA_u_W0 + btc%uh_crvW * u**2)
3048 uhbt = (u - btc%uBT_WW) * btc%FA_u_WW + btc%uh_WW
3055 function uhbt_to_ubt(uhbt, BTC, US, guess)
result(ubt)
3056 real,
intent(in) :: uhbt
3062 real,
optional,
intent(in) :: guess
3067 real :: ubt_min, ubt_max, uhbt_err, derr_du
3068 real :: uherr_min, uherr_max
3069 real,
parameter :: tol = 1.0e-10
3072 real,
parameter :: vs1 = 1.25
3073 real,
parameter :: vs2 = 2.0
3075 integer :: itt, max_itt = 20
3078 if (uhbt == 0.0)
then
3080 elseif (uhbt < btc%uh_EE)
then
3081 ubt = btc%uBT_EE + (uhbt - btc%uh_EE) / btc%FA_u_EE
3082 elseif (uhbt < 0.0)
then
3085 ubt_min = btc%uBT_EE ; uherr_min = btc%uh_EE - uhbt
3086 ubt_max = 0.0 ; uherr_max = -uhbt
3088 ubt = btc%uBT_EE * (uhbt / btc%uh_EE)
3090 uhbt_err = ubt * (btc%FA_u_E0 + btc%uh_crvE * ubt**2) - uhbt
3092 if (abs(uhbt_err) < tol*abs(uhbt))
exit
3093 if (uhbt_err > 0.0)
then ; ubt_max = ubt ; uherr_max = uhbt_err ;
endif
3094 if (uhbt_err < 0.0)
then ; ubt_min = ubt ; uherr_min = uhbt_err ;
endif
3096 derr_du = btc%FA_u_E0 + 3.0 * btc%uh_crvE * ubt**2
3097 if ((uhbt_err >= derr_du*(ubt - ubt_min)) .or. &
3098 (-uhbt_err >= derr_du*(ubt_max - ubt)) .or. (derr_du <= 0.0))
then
3100 ubt = ubt_max + (ubt_min-ubt_max) * (uherr_max / (uherr_max-uherr_min))
3102 ubt = ubt - uhbt_err / derr_du
3103 if (abs(uhbt_err) < (0.01*tol)*abs(ubt_min*derr_du))
exit
3106 elseif (uhbt <= btc%uh_WW)
then
3108 ubt_min = 0.0 ; uherr_min = -uhbt
3109 ubt_max = btc%uBT_WW ; uherr_max = btc%uh_WW - uhbt
3111 ubt = btc%uBT_WW * (uhbt / btc%uh_WW)
3113 uhbt_err = ubt * (btc%FA_u_W0 + btc%uh_crvW * ubt**2) - uhbt
3115 if (abs(uhbt_err) < tol*abs(uhbt))
exit
3116 if (uhbt_err > 0.0)
then ; ubt_max = ubt ; uherr_max = uhbt_err ;
endif
3117 if (uhbt_err < 0.0)
then ; ubt_min = ubt ; uherr_min = uhbt_err ;
endif
3119 derr_du = btc%FA_u_W0 + 3.0 * btc%uh_crvW * ubt**2
3120 if ((uhbt_err >= derr_du*(ubt - ubt_min)) .or. &
3121 (-uhbt_err >= derr_du*(ubt_max - ubt)) .or. (derr_du <= 0.0))
then
3123 ubt = ubt_min + (ubt_max-ubt_min) * (-uherr_min / (uherr_max-uherr_min))
3125 ubt = ubt - uhbt_err / derr_du
3126 if (abs(uhbt_err) < (0.01*tol)*(ubt_max*derr_du))
exit
3130 ubt = btc%uBT_WW + (uhbt - btc%uh_WW) / btc%FA_u_WW
3133 if (
present(guess))
then
3134 dvel = abs(ubt) - vs1*abs(guess)
3135 if (dvel > 0.0)
then
3136 if (dvel < 40.0 * (abs(guess)*(vs2-vs1)) )
then
3137 vsr = vs2 - (vs2-vs1) * exp(-dvel / (abs(guess)*(vs2-vs1)))
3141 ubt = sign(vsr * guess, ubt)
3148 function find_vhbt(v, BTC, US)
result(vhbt)
3149 real,
intent(in) :: v
3158 elseif (v < btc%vBT_NN)
then
3159 vhbt = (v - btc%vBT_NN) * btc%FA_v_NN + btc%vh_NN
3160 elseif (v < 0.0)
then
3161 vhbt = v * (btc%FA_v_N0 + btc%vh_crvN * v**2)
3162 elseif (v <= btc%vBT_SS)
then
3163 vhbt = v * (btc%FA_v_S0 + btc%vh_crvS * v**2)
3165 vhbt = (v - btc%vBT_SS) * btc%FA_v_SS + btc%vh_SS
3172 function vhbt_to_vbt(vhbt, BTC, US, guess)
result(vbt)
3173 real,
intent(in) :: vhbt
3179 real,
optional,
intent(in) :: guess
3184 real :: vbt_min, vbt_max, vhbt_err, derr_dv
3185 real :: vherr_min, vherr_max
3186 real,
parameter :: tol = 1.0e-10
3189 real,
parameter :: vs1 = 1.25
3190 real,
parameter :: vs2 = 2.0
3192 integer :: itt, max_itt = 20
3195 if (vhbt == 0.0)
then
3197 elseif (vhbt < btc%vh_NN)
then
3198 vbt = btc%vBT_NN + (vhbt - btc%vh_NN) / btc%FA_v_NN
3199 elseif (vhbt < 0.0)
then
3202 vbt_min = btc%vBT_NN ; vherr_min = btc%vh_NN - vhbt
3203 vbt_max = 0.0 ; vherr_max = -vhbt
3205 vbt = btc%vBT_NN * (vhbt / btc%vh_NN)
3207 vhbt_err = vbt * (btc%FA_v_N0 + btc%vh_crvN * vbt**2) - vhbt
3209 if (abs(vhbt_err) < tol*abs(vhbt))
exit
3210 if (vhbt_err > 0.0)
then ; vbt_max = vbt ; vherr_max = vhbt_err ;
endif
3211 if (vhbt_err < 0.0)
then ; vbt_min = vbt ; vherr_min = vhbt_err ;
endif
3213 derr_dv = btc%FA_v_N0 + 3.0 * btc%vh_crvN * vbt**2
3214 if ((vhbt_err >= derr_dv*(vbt - vbt_min)) .or. &
3215 (-vhbt_err >= derr_dv*(vbt_max - vbt)) .or. (derr_dv <= 0.0))
then
3217 vbt = vbt_max + (vbt_min-vbt_max) * (vherr_max / (vherr_max-vherr_min))
3219 vbt = vbt - vhbt_err / derr_dv
3220 if (abs(vhbt_err) < (0.01*tol)*abs(derr_dv*vbt_min))
exit
3223 elseif (vhbt <= btc%vh_SS)
then
3225 vbt_min = 0.0 ; vherr_min = -vhbt
3226 vbt_max = btc%vBT_SS ; vherr_max = btc%vh_SS - vhbt
3228 vbt = btc%vBT_SS * (vhbt / btc%vh_SS)
3230 vhbt_err = vbt * (btc%FA_v_S0 + btc%vh_crvS * vbt**2) - vhbt
3232 if (abs(vhbt_err) < tol*abs(vhbt))
exit
3233 if (vhbt_err > 0.0)
then ; vbt_max = vbt ; vherr_max = vhbt_err ;
endif
3234 if (vhbt_err < 0.0)
then ; vbt_min = vbt ; vherr_min = vhbt_err ;
endif
3236 derr_dv = btc%FA_v_S0 + 3.0 * btc%vh_crvS * vbt**2
3237 if ((vhbt_err >= derr_dv*(vbt - vbt_min)) .or. &
3238 (-vhbt_err >= derr_dv*(vbt_max - vbt)) .or. (derr_dv <= 0.0))
then
3240 vbt = vbt_min + (vbt_max-vbt_min) * (-vherr_min / (vherr_max-vherr_min))
3242 vbt = vbt - vhbt_err / derr_dv
3243 if (abs(vhbt_err) < (0.01*tol)*(vbt_max*derr_dv))
exit
3247 vbt = btc%vBT_SS + (vhbt - btc%vh_SS) / btc%FA_v_SS
3250 if (
present(guess))
then
3251 dvel = abs(vbt) - vs1*abs(guess)
3252 if (dvel > 0.0)
then
3253 if (dvel < 40.0 * (abs(guess)*(vs2-vs1)) )
then
3254 vsr = vs2 - (vs2-vs1) * exp(-dvel / (abs(guess)*(vs2-vs1)))
3258 vbt = sign(guess * vsr, vbt)
3280 integer,
optional,
intent(in) :: halo
3283 real,
dimension(SZIBW_(MS),SZJW_(MS)) :: &
3284 u_polarity, uBT_EE, uBT_WW, FA_u_EE, FA_u_E0, FA_u_W0, FA_u_WW
3285 real,
dimension(SZIW_(MS),SZJBW_(MS)) :: &
3286 v_polarity, vBT_NN, vBT_SS, FA_v_NN, FA_v_N0, FA_v_S0, FA_v_SS
3287 real,
parameter :: C1_3 = 1.0/3.0
3288 integer :: i, j, is, ie, js, je, hs
3289 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
3290 hs = 1 ;
if (
present(halo)) hs = max(halo,0)
3297 do j=js-hs,je+hs ;
do i=is-hs-1,ie+hs
3298 u_polarity(i,j) = 1.0
3299 ubt_ee(i,j) = 0.0 ; ubt_ww(i,j) = 0.0
3300 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
3303 do j=js-hs-1,je+hs ;
do i=is-hs,ie+hs
3304 v_polarity(i,j) = 1.0
3305 vbt_nn(i,j) = 0.0 ; vbt_ss(i,j) = 0.0
3306 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
3309 do j=js,je;
do i=is-1,ie
3310 ubt_ee(i,j) = bt_cont%uBT_EE(i,j) ; ubt_ww(i,j) = bt_cont%uBT_WW(i,j)
3311 fa_u_ee(i,j) = bt_cont%FA_u_EE(i,j) ; fa_u_e0(i,j) = bt_cont%FA_u_E0(i,j)
3312 fa_u_w0(i,j) = bt_cont%FA_u_W0(i,j) ; fa_u_ww(i,j) = bt_cont%FA_u_WW(i,j)
3315 do j=js-1,je;
do i=is,ie
3316 vbt_nn(i,j) = bt_cont%vBT_NN(i,j) ; vbt_ss(i,j) = bt_cont%vBT_SS(i,j)
3317 fa_v_nn(i,j) = bt_cont%FA_v_NN(i,j) ; fa_v_n0(i,j) = bt_cont%FA_v_N0(i,j)
3318 fa_v_s0(i,j) = bt_cont%FA_v_S0(i,j) ; fa_v_ss(i,j) = bt_cont%FA_v_SS(i,j)
3325 call create_group_pass(bt_cont%pass_polarity_BT, u_polarity, v_polarity, bt_domain)
3329 call create_group_pass(bt_cont%pass_FA_uv, fa_u_ee, fa_v_nn, bt_domain, to_all+scalar_pair)
3330 call create_group_pass(bt_cont%pass_FA_uv, fa_u_e0, fa_v_n0, bt_domain, to_all+scalar_pair)
3331 call create_group_pass(bt_cont%pass_FA_uv, fa_u_w0, fa_v_s0, bt_domain, to_all+scalar_pair)
3332 call create_group_pass(bt_cont%pass_FA_uv, fa_u_ww, fa_v_ss, bt_domain, to_all+scalar_pair)
3335 call do_group_pass(bt_cont%pass_polarity_BT, bt_domain)
3336 call do_group_pass(bt_cont%pass_FA_uv, bt_domain)
3345 do j=js-hs,je+hs ;
do i=is-hs-1,ie+hs
3346 btcl_u(i,j)%FA_u_EE = fa_u_ee(i,j) ; btcl_u(i,j)%FA_u_E0 = fa_u_e0(i,j)
3347 btcl_u(i,j)%FA_u_W0 = fa_u_w0(i,j) ; btcl_u(i,j)%FA_u_WW = fa_u_ww(i,j)
3348 btcl_u(i,j)%uBT_EE = ubt_ee(i,j) ; btcl_u(i,j)%uBT_WW = ubt_ww(i,j)
3350 if (u_polarity(i,j) < 0.0)
then
3351 call swap(btcl_u(i,j)%FA_u_EE, btcl_u(i,j)%FA_u_WW)
3352 call swap(btcl_u(i,j)%FA_u_E0, btcl_u(i,j)%FA_u_W0)
3353 call swap(btcl_u(i,j)%uBT_EE, btcl_u(i,j)%uBT_WW)
3356 btcl_u(i,j)%uh_EE = btcl_u(i,j)%uBT_EE * &
3357 (c1_3 * (2.0*btcl_u(i,j)%FA_u_E0 + btcl_u(i,j)%FA_u_EE))
3358 btcl_u(i,j)%uh_WW = btcl_u(i,j)%uBT_WW * &
3359 (c1_3 * (2.0*btcl_u(i,j)%FA_u_W0 + btcl_u(i,j)%FA_u_WW))
3361 btcl_u(i,j)%uh_crvE = 0.0 ; btcl_u(i,j)%uh_crvW = 0.0
3362 if (abs(btcl_u(i,j)%uBT_WW) > 0.0) btcl_u(i,j)%uh_crvW = &
3363 (c1_3 * (btcl_u(i,j)%FA_u_WW - btcl_u(i,j)%FA_u_W0)) / btcl_u(i,j)%uBT_WW**2
3364 if (abs(btcl_u(i,j)%uBT_EE) > 0.0) btcl_u(i,j)%uh_crvE = &
3365 (c1_3 * (btcl_u(i,j)%FA_u_EE - btcl_u(i,j)%FA_u_E0)) / btcl_u(i,j)%uBT_EE**2
3368 do j=js-hs-1,je+hs ;
do i=is-hs,ie+hs
3369 btcl_v(i,j)%FA_v_NN = fa_v_nn(i,j) ; btcl_v(i,j)%FA_v_N0 = fa_v_n0(i,j)
3370 btcl_v(i,j)%FA_v_S0 = fa_v_s0(i,j) ; btcl_v(i,j)%FA_v_SS = fa_v_ss(i,j)
3371 btcl_v(i,j)%vBT_NN = vbt_nn(i,j) ; btcl_v(i,j)%vBT_SS = vbt_ss(i,j)
3373 if (v_polarity(i,j) < 0.0)
then
3374 call swap(btcl_v(i,j)%FA_v_NN, btcl_v(i,j)%FA_v_SS)
3375 call swap(btcl_v(i,j)%FA_v_N0, btcl_v(i,j)%FA_v_S0)
3376 call swap(btcl_v(i,j)%vBT_NN, btcl_v(i,j)%vBT_SS)
3379 btcl_v(i,j)%vh_NN = btcl_v(i,j)%vBT_NN * &
3380 (c1_3 * (2.0*btcl_v(i,j)%FA_v_N0 + btcl_v(i,j)%FA_v_NN))
3381 btcl_v(i,j)%vh_SS = btcl_v(i,j)%vBT_SS * &
3382 (c1_3 * (2.0*btcl_v(i,j)%FA_v_S0 + btcl_v(i,j)%FA_v_SS))
3384 btcl_v(i,j)%vh_crvN = 0.0 ; btcl_v(i,j)%vh_crvS = 0.0
3385 if (abs(btcl_v(i,j)%vBT_SS) > 0.0) btcl_v(i,j)%vh_crvS = &
3386 (c1_3 * (btcl_v(i,j)%FA_v_SS - btcl_v(i,j)%FA_v_S0)) / btcl_v(i,j)%vBT_SS**2
3387 if (abs(btcl_v(i,j)%vBT_NN) > 0.0) btcl_v(i,j)%vh_crvN = &
3388 (c1_3 * (btcl_v(i,j)%FA_v_NN - btcl_v(i,j)%FA_v_N0)) / btcl_v(i,j)%vBT_NN**2
3399 real,
dimension(SZIBW_(MS),SZJW_(MS)), &
3401 real,
dimension(SZIBW_(MS),SZJW_(MS)), &
3404 real,
dimension(SZIW_(MS),SZJBW_(MS)), &
3406 real,
dimension(SZIW_(MS),SZJBW_(MS)), &
3410 intent(out) :: BTCL_u
3412 intent(out) :: BTCL_v
3415 integer,
optional,
intent(in) :: halo
3418 real,
dimension(SZIBW_(MS),SZJW_(MS)) :: &
3419 u_polarity, uBT_EE, uBT_WW, FA_u_EE, FA_u_E0, FA_u_W0, FA_u_WW
3420 real,
dimension(SZIW_(MS),SZJBW_(MS)) :: &
3421 v_polarity, vBT_NN, vBT_SS, FA_v_NN, FA_v_N0, FA_v_S0, FA_v_SS
3422 real,
parameter :: C1_3 = 1.0/3.0
3423 integer :: i, j, is, ie, js, je, hs
3424 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
3425 hs = 1 ;
if (
present(halo)) hs = max(halo,0)
3428 do j=js-hs,je+hs ;
do i=is-hs-1,ie+hs
3429 if ((ubt(i,j) > btcl_u(i,j)%uBT_WW) .and. (uhbt(i,j) > btcl_u(i,j)%uh_WW))
then
3431 btcl_u(i,j)%ubt_WW = ubt(i,j)
3432 if (3.0*uhbt(i,j) < 2.0*ubt(i,j) * btcl_u(i,j)%FA_u_W0)
then
3434 btcl_u(i,j)%uh_crvW = (uhbt(i,j) - ubt(i,j) * btcl_u(i,j)%FA_u_W0) / ubt(i,j)**3
3436 btcl_u(i,j)%FA_u_W0 = 1.5*uhbt(i,j) / ubt(i,j)
3437 btcl_u(i,j)%uh_crvW = -0.5*uhbt(i,j) / ubt(i,j)**3
3439 btcl_u(i,j)%uh_WW = uhbt(i,j)
3442 elseif ((ubt(i,j) < btcl_u(i,j)%uBT_EE) .and. (uhbt(i,j) < btcl_u(i,j)%uh_EE))
then
3444 btcl_u(i,j)%ubt_EE = ubt(i,j)
3445 if (3.0*uhbt(i,j) < 2.0*ubt(i,j) * btcl_u(i,j)%FA_u_E0)
then
3447 btcl_u(i,j)%uh_crvE = (uhbt(i,j) - ubt(i,j) * btcl_u(i,j)%FA_u_E0) / ubt(i,j)**3
3449 btcl_u(i,j)%FA_u_E0 = 1.5*uhbt(i,j) / ubt(i,j)
3450 btcl_u(i,j)%uh_crvE = -0.5*uhbt(i,j) / ubt(i,j)**3
3452 btcl_u(i,j)%uh_EE = uhbt(i,j)
3458 do j=js-hs-1,je+hs ;
do i=is-hs,ie+hs
3459 if ((vbt(i,j) > btcl_v(i,j)%vBT_SS) .and. (vhbt(i,j) > btcl_v(i,j)%vh_SS))
then
3461 btcl_v(i,j)%vbt_SS = vbt(i,j)
3462 if (3.0*vhbt(i,j) < 2.0*vbt(i,j) * btcl_v(i,j)%FA_v_S0)
then
3464 btcl_v(i,j)%vh_crvS = (vhbt(i,j) - vbt(i,j) * btcl_v(i,j)%FA_v_S0) / vbt(i,j)**3
3466 btcl_v(i,j)%FA_v_S0 = 1.5*vhbt(i,j) / (vbt(i,j))
3467 btcl_v(i,j)%vh_crvS = -0.5*vhbt(i,j) / vbt(i,j)**3
3469 btcl_v(i,j)%vh_SS = vhbt(i,j)
3472 elseif ((vbt(i,j) < btcl_v(i,j)%vBT_NN) .and. (vhbt(i,j) < btcl_v(i,j)%vh_NN))
then
3474 btcl_v(i,j)%vbt_NN = vbt(i,j)
3475 if (3.0*vhbt(i,j) < 2.0*vbt(i,j) * btcl_v(i,j)%FA_v_N0)
then
3477 btcl_v(i,j)%vh_crvN = (vhbt(i,j) - vbt(i,j) * btcl_v(i,j)%FA_v_N0) / vbt(i,j)**3
3479 btcl_v(i,j)%FA_v_N0 = 1.5*vhbt(i,j) / (vbt(i,j))
3480 btcl_v(i,j)%vh_crvN = -0.5*vhbt(i,j) / vbt(i,j)**3
3482 btcl_v(i,j)%vh_NN = vhbt(i,j)
3497 real,
dimension(MS%isdw-1:MS%iedw,MS%jsdw:MS%jedw), &
3499 real,
dimension(MS%isdw:MS%iedw,MS%jsdw-1:MS%jedw), &
3503 integer,
optional,
intent(in) :: halo
3504 logical,
optional,
intent(in) :: maximize
3509 integer :: i, j, is, ie, js, je, hs
3510 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
3511 hs = 1 ;
if (
present(halo)) hs = max(halo,0)
3512 find_max = .false. ;
if (
present(maximize)) find_max = maximize
3515 do j=js-hs,je+hs ;
do i=is-1-hs,ie+hs
3516 datu(i,j) = max(bt_cont%FA_u_EE(i,j), bt_cont%FA_u_E0(i,j), &
3517 bt_cont%FA_u_W0(i,j), bt_cont%FA_u_WW(i,j))
3519 do j=js-1-hs,je+hs ;
do i=is-hs,ie+hs
3520 datv(i,j) = max(bt_cont%FA_v_NN(i,j), bt_cont%FA_v_N0(i,j), &
3521 bt_cont%FA_v_S0(i,j), bt_cont%FA_v_SS(i,j))
3524 do j=js-hs,je+hs ;
do i=is-1-hs,ie+hs
3525 datu(i,j) = 0.5 * (bt_cont%FA_u_E0(i,j) + bt_cont%FA_u_W0(i,j))
3527 do j=js-1-hs,je+hs ;
do i=is-hs,ie+hs
3528 datv(i,j) = 0.5 * (bt_cont%FA_v_N0(i,j) + bt_cont%FA_v_S0(i,j))
3535 subroutine swap(a,b)
3536 real,
intent(inout) :: a
3537 real,
intent(inout) :: b
3539 tmp = a ; a = b ; b = tmp
3544 subroutine find_face_areas(Datu, Datv, G, GV, US, CS, MS, eta, halo, add_max)
3546 real,
dimension(MS%isdw-1:MS%iedw,MS%jsdw:MS%jedw), &
3548 real,
dimension(MS%isdw:MS%iedw,MS%jsdw-1:MS%jedw), &
3555 real,
dimension(MS%isdw:MS%iedw,MS%jsdw:MS%jedw), &
3556 optional,
intent(in) :: eta
3558 integer,
optional,
intent(in) :: halo
3559 real,
optional,
intent(in) :: add_max
3564 integer :: i, j, is, ie, js, je, hs
3565 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
3566 hs = 1 ;
if (
present(halo)) hs = max(halo,0)
3570 if (
present(eta))
then
3572 if (gv%Boussinesq)
then
3574 do j=js-hs,je+hs ;
do i=is-1-hs,ie+hs
3575 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)
3576 datu(i,j) = 0.0 ;
if ((h1 > 0.0) .and. (h2 > 0.0)) &
3577 datu(i,j) = cs%dy_Cu(i,j) * (2.0 * h1 * h2) / (h1 + h2)
3581 do j=js-1-hs,je+hs ;
do i=is-hs,ie+hs
3582 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)
3583 datv(i,j) = 0.0 ;
if ((h1 > 0.0) .and. (h2 > 0.0)) &
3584 datv(i,j) = cs%dx_Cv(i,j) * (2.0 * h1 * h2) / (h1 + h2)
3589 do j=js-hs,je+hs ;
do i=is-1-hs,ie+hs
3590 datu(i,j) = 0.0 ;
if ((eta(i,j) > 0.0) .and. (eta(i+1,j) > 0.0)) &
3591 datu(i,j) = cs%dy_Cu(i,j) * (2.0 * eta(i,j) * eta(i+1,j)) / &
3592 (eta(i,j) + eta(i+1,j))
3596 do j=js-1-hs,je+hs ;
do i=is-hs,ie+hs
3597 datv(i,j) = 0.0 ;
if ((eta(i,j) > 0.0) .and. (eta(i,j+1) > 0.0)) &
3598 datv(i,j) = cs%dx_Cv(i,j) * (2.0 * eta(i,j) * eta(i,j+1)) / &
3599 (eta(i,j) + eta(i,j+1))
3603 elseif (
present(add_max))
then
3605 do j=js-hs,je+hs ;
do i=is-1-hs,ie+hs
3606 datu(i,j) = cs%dy_Cu(i,j) * gv%Z_to_H * &
3607 (max(cs%bathyT(i+1,j), cs%bathyT(i,j)) + add_max)
3610 do j=js-1-hs,je+hs ;
do i=is-hs,ie+hs
3611 datv(i,j) = cs%dx_Cv(i,j) * gv%Z_to_H * &
3612 (max(cs%bathyT(i,j+1), cs%bathyT(i,j)) + add_max)
3616 do j=js-hs,je+hs ;
do i=is-1-hs,ie+hs
3619 if (cs%bathyT(i+1,j)+cs%bathyT(i,j)>0.) &
3620 datu(i,j) = 2.0*cs%dy_Cu(i,j) * gv%Z_to_H * &
3621 (cs%bathyT(i+1,j) * cs%bathyT(i,j)) / &
3622 (cs%bathyT(i+1,j) + cs%bathyT(i,j))
3625 do j=js-1-hs,je+hs ;
do i=is-hs,ie+hs
3628 if (cs%bathyT(i,j+1)+cs%bathyT(i,j)>0.) &
3629 datv(i,j) = 2.0*cs%dx_Cv(i,j) * gv%Z_to_H * &
3630 (cs%bathyT(i,j+1) * cs%bathyT(i,j)) / &
3631 (cs%bathyT(i,j+1) + cs%bathyT(i,j))
3645 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
3646 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: eta
3648 logical,
intent(in) :: set_cor
3656 real :: h_tot(szi_(g))
3657 real :: eta_h(szi_(g))
3661 integer :: is, ie, js, je, nz, i, j, k
3662 real,
parameter :: frac_cor = 0.25
3663 real,
parameter :: slow_rate = 0.125
3665 if (.not.
associated(cs))
call mom_error(fatal,
"bt_mass_source: "// &
3666 "Module MOM_barotropic must be initialized before it is used.")
3667 if (.not.cs%split)
return
3669 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
3673 do i=is,ie ; h_tot(i) = h(i,j,1) ;
enddo
3674 if (gv%Boussinesq)
then
3675 do i=is,ie ; eta_h(i) = h(i,j,1) - g%bathyT(i,j)*gv%Z_to_H ;
enddo
3677 do i=is,ie ; eta_h(i) = h(i,j,1) ;
enddo
3679 do k=2,nz ;
do i=is,ie
3680 eta_h(i) = eta_h(i) + h(i,j,k)
3681 h_tot(i) = h_tot(i) + h(i,j,k)
3686 d_eta = eta_h(i) - eta(i,j)
3687 cs%eta_cor(i,j) = d_eta
3691 d_eta = eta_h(i) - eta(i,j)
3692 cs%eta_cor(i,j) = cs%eta_cor(i,j) + d_eta
3702 subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, &
3703 restart_CS, calc_dtbt, BT_cont, tides_CSp)
3707 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
3709 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
3711 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
3713 real,
dimension(SZI_(G),SZJ_(G)), &
3716 type(time_type),
target,
intent(in) :: time
3718 type(
diag_ctrl),
target,
intent(inout) :: diag
3723 logical,
intent(out) :: calc_dtbt
3730 pointer :: tides_csp
3734 #include "version_variable.h"
3736 character(len=40) :: mdl =
"MOM_barotropic"
3737 real :: datu(szibs_(g),szj_(g))
3738 real :: datv(szi_(g),szjbs_(g))
3739 real :: gtot_estimate
3744 real :: wave_drag_scale
3746 character(len=200) :: inputdir
3747 character(len=200) :: wave_drag_file
3749 character(len=80) :: wave_drag_var
3755 real,
allocatable,
dimension(:,:) :: lin_drag_h
3757 type(group_pass_type) :: pass_static_data, pass_q_d_cor
3758 type(group_pass_type) :: pass_bt_hbt_btav, pass_a_polarity
3759 logical :: apply_bt_drag, use_bt_cont_type
3760 character(len=48) :: thickness_units, flux_units
3761 character*(40) :: hvel_str
3762 integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz
3763 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
3764 integer :: isdw, iedw, jsdw, jedw
3766 integer :: wd_halos(2), bt_halo_sz
3767 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
3768 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
3769 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
3770 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
3771 ms%isdw = g%isd ; ms%iedw = g%ied ; ms%jsdw = g%jsd ; ms%jedw = g%jed
3773 if (cs%module_is_initialized)
then
3774 call mom_error(warning,
"barotropic_init called with a control structure "// &
3775 "that has already been initialized.")
3778 cs%module_is_initialized = .true.
3780 cs%diag => diag ; cs%Time => time
3781 if (
present(tides_csp))
then
3782 if (
associated(tides_csp)) cs%tides_CSp => tides_csp
3787 call get_param(param_file, mdl,
"SPLIT", cs%split, &
3788 "Use the split time stepping if true.", default=.true.)
3789 if (.not.cs%split)
return
3791 call get_param(param_file, mdl,
"BOUND_BT_CORRECTION", cs%bound_BT_corr, &
3792 "If true, the corrective pseudo mass-fluxes into the "//&
3793 "barotropic solver are limited to values that require "//&
3794 "less than maxCFL_BT_cont to be accommodated.",default=.false.)
3795 call get_param(param_file, mdl,
"BT_CONT_CORR_BOUNDS", cs%BT_cont_bounds, &
3796 "If true, and BOUND_BT_CORRECTION is true, use the "//&
3797 "BT_cont_type variables to set limits determined by "//&
3798 "MAXCFL_BT_CONT on the CFL number of the velocities "//&
3799 "that are likely to be driven by the corrective mass fluxes.", &
3801 call get_param(param_file, mdl,
"ADJUST_BT_CONT", cs%adjust_BT_cont, &
3802 "If true, adjust the curve fit to the BT_cont type "//&
3803 "that is used by the barotropic solver to match the "//&
3804 "transport about which the flow is being linearized.", default=.false.)
3805 call get_param(param_file, mdl,
"GRADUAL_BT_ICS", cs%gradual_BT_ICs, &
3806 "If true, adjust the initial conditions for the "//&
3807 "barotropic solver to the values from the layered "//&
3808 "solution over a whole timestep instead of instantly. "//&
3809 "This is a decent approximation to the inclusion of "//&
3810 "sum(u dh_dt) while also correcting for truncation errors.", &
3812 call get_param(param_file, mdl,
"BT_USE_VISC_REM_U_UH0", cs%visc_rem_u_uh0, &
3813 "If true, use the viscous remnants when estimating the "//&
3814 "barotropic velocities that were used to calculate uh0 "//&
3815 "and vh0. False is probably the better choice.", default=.false.)
3816 call get_param(param_file, mdl,
"BT_USE_WIDE_HALOS", cs%use_wide_halos, &
3817 "If true, use wide halos and march in during the "//&
3818 "barotropic time stepping for efficiency.", default=.true., &
3820 call get_param(param_file, mdl,
"BTHALO", bt_halo_sz, &
3821 "The minimum halo size for the barotropic solver.", default=0, &
3823 #ifdef STATIC_MEMORY_
3824 if ((bt_halo_sz > 0) .and. (bt_halo_sz /= bthalo_))
call mom_error(fatal, &
3825 "barotropic_init: Run-time values of BTHALO must agree with the "//&
3826 "macro BTHALO_ with STATIC_MEMORY_.")
3827 wd_halos(1) = whaloi_+nihalo_ ; wd_halos(2) = whaloj_+njhalo_
3829 wd_halos(1) = bt_halo_sz; wd_halos(2) = bt_halo_sz
3831 call log_param(param_file, mdl,
"!BT x-halo", wd_halos(1), &
3832 "The barotropic x-halo size that is actually used.", &
3834 call log_param(param_file, mdl,
"!BT y-halo", wd_halos(2), &
3835 "The barotropic y-halo size that is actually used.", &
3838 call get_param(param_file, mdl,
"USE_BT_CONT_TYPE", use_bt_cont_type, &
3839 "If true, use a structure with elements that describe "//&
3840 "effective face areas from the summed continuity solver "//&
3841 "as a function the barotropic flow in coupling between "//&
3842 "the barotropic and baroclinic flow. This is only used "//&
3843 "if SPLIT is true. \n", default=.true.)
3844 call get_param(param_file, mdl,
"NONLINEAR_BT_CONTINUITY", &
3845 cs%Nonlinear_continuity, &
3846 "If true, use nonlinear transports in the barotropic "//&
3847 "continuity equation. This does not apply if "//&
3848 "USE_BT_CONT_TYPE is true.", default=.false.)
3849 cs%Nonlin_cont_update_period = 1
3850 if (cs%Nonlinear_continuity) &
3851 call get_param(param_file, mdl,
"NONLIN_BT_CONT_UPDATE_PERIOD", &
3852 cs%Nonlin_cont_update_period, &
3853 "If NONLINEAR_BT_CONTINUITY is true, this is the number "//&
3854 "of barotropic time steps between updates to the face "//&
3855 "areas, or 0 to update only before the barotropic stepping.",&
3856 units=
"nondim", default=1)
3857 call get_param(param_file, mdl,
"BT_PROJECT_VELOCITY", cs%BT_project_velocity,&
3858 "If true, step the barotropic velocity first and project "//&
3859 "out the velocity tendency by 1+BEBT when calculating the "//&
3860 "transport. The default (false) is to use a predictor "//&
3861 "continuity step to find the pressure field, and then "//&
3862 "to do a corrector continuity step using a weighted "//&
3863 "average of the old and new velocities, with weights "//&
3864 "of (1-BEBT) and BEBT.", default=.false.)
3866 call get_param(param_file, mdl,
"DYNAMIC_SURFACE_PRESSURE", cs%dynamic_psurf, &
3867 "If true, add a dynamic pressure due to a viscous ice "//&
3868 "shelf, for instance.", default=.false.)
3869 if (cs%dynamic_psurf)
then
3870 call get_param(param_file, mdl,
"ICE_LENGTH_DYN_PSURF", cs%ice_strength_length, &
3871 "The length scale at which the Rayleigh damping rate due "//&
3872 "to the ice strength should be the same as if a Laplacian "//&
3873 "were applied, if DYNAMIC_SURFACE_PRESSURE is true.", &
3874 units=
"m", default=1.0e4, scale=us%m_to_L)
3875 call get_param(param_file, mdl,
"DEPTH_MIN_DYN_PSURF", cs%Dmin_dyn_psurf, &
3876 "The minimum depth to use in limiting the size of the "//&
3877 "dynamic surface pressure for stability, if "//&
3878 "DYNAMIC_SURFACE_PRESSURE is true..", &
3879 units=
"m", default=1.0e-6, scale=us%m_to_Z)
3880 call get_param(param_file, mdl,
"CONST_DYN_PSURF", cs%const_dyn_psurf, &
3881 "The constant that scales the dynamic surface pressure, "//&
3882 "if DYNAMIC_SURFACE_PRESSURE is true. Stable values "//&
3883 "are < ~1.0.", units=
"nondim", default=0.9)
3886 call get_param(param_file, mdl,
"TIDES", cs%tides, &
3887 "If true, apply tidal momentum forcing.", default=.false.)
3888 call get_param(param_file, mdl,
"SADOURNY", cs%Sadourny, &
3889 "If true, the Coriolis terms are discretized with the "//&
3890 "Sadourny (1975) energy conserving scheme, otherwise "//&
3891 "the Arakawa & Hsu scheme is used. If the internal "//&
3892 "deformation radius is not resolved, the Sadourny scheme "//&
3893 "should probably be used.", default=.true.)
3895 call get_param(param_file, mdl,
"BT_THICK_SCHEME", hvel_str, &
3896 "A string describing the scheme that is used to set the "//&
3897 "open face areas used for barotropic transport and the "//&
3898 "relative weights of the accelerations. Valid values are:\n"//&
3899 "\t ARITHMETIC - arithmetic mean layer thicknesses \n"//&
3900 "\t HARMONIC - harmonic mean layer thicknesses \n"//&
3901 "\t HYBRID (the default) - use arithmetic means for \n"//&
3902 "\t layers above the shallowest bottom, the harmonic \n"//&
3903 "\t mean for layers below, and a weighted average for \n"//&
3904 "\t layers that straddle that depth \n"//&
3905 "\t FROM_BT_CONT - use the average thicknesses kept \n"//&
3906 "\t in the h_u and h_v fields of the BT_cont_type", &
3908 select case (hvel_str)
3914 call mom_mesg(
'barotropic_init: BT_THICK_SCHEME ="'//trim(hvel_str)//
'"', 0)
3915 call mom_error(fatal,
"barotropic_init: Unrecognized setting "// &
3916 "#define BT_THICK_SCHEME "//trim(hvel_str)//
" found in input file.")
3918 if ((cs%hvel_scheme ==
from_bt_cont) .and. .not.use_bt_cont_type) &
3919 call mom_error(fatal,
"barotropic_init: BT_THICK_SCHEME FROM_BT_CONT "//&
3920 "can only be used if USE_BT_CONT_TYPE is defined.")
3922 call get_param(param_file, mdl,
"BT_STRONG_DRAG", cs%strong_drag, &
3923 "If true, use a stronger estimate of the retarding "//&
3924 "effects of strong bottom drag, by making it implicit "//&
3925 "with the barotropic time-step instead of implicit with "//&
3926 "the baroclinic time-step and dividing by the number of "//&
3927 "barotropic steps.", default=.false.)
3928 call get_param(param_file, mdl,
"BT_LINEAR_WAVE_DRAG", cs%linear_wave_drag, &
3929 "If true, apply a linear drag to the barotropic velocities, "//&
3930 "using rates set by lin_drag_u & _v divided by the depth of "//&
3931 "the ocean. This was introduced to facilitate tide modeling.", &
3933 call get_param(param_file, mdl,
"BT_WAVE_DRAG_FILE", wave_drag_file, &
3934 "The name of the file with the barotropic linear wave drag "//&
3935 "piston velocities.", default=
"", do_not_log=.not.cs%linear_wave_drag)
3936 call get_param(param_file, mdl,
"BT_WAVE_DRAG_VAR", wave_drag_var, &
3937 "The name of the variable in BT_WAVE_DRAG_FILE with the "//&
3938 "barotropic linear wave drag piston velocities at h points.", &
3939 default=
"rH", do_not_log=.not.cs%linear_wave_drag)
3940 call get_param(param_file, mdl,
"BT_WAVE_DRAG_SCALE", wave_drag_scale, &
3941 "A scaling factor for the barotropic linear wave drag "//&
3942 "piston velocities.", default=1.0, units=
"nondim", &
3943 do_not_log=.not.cs%linear_wave_drag)
3945 call get_param(param_file, mdl,
"CLIP_BT_VELOCITY", cs%clip_velocity, &
3946 "If true, limit any velocity components that exceed "//&
3947 "CFL_TRUNCATE. This should only be used as a desperate "//&
3948 "debugging measure.", default=.false.)
3949 call get_param(param_file, mdl,
"CFL_TRUNCATE", cs%CFL_trunc, &
3950 "The value of the CFL number that will cause velocity "//&
3951 "components to be truncated; instability can occur past 0.5.", &
3952 units=
"nondim", default=0.5, do_not_log=.not.cs%clip_velocity)
3953 call get_param(param_file, mdl,
"MAXVEL", cs%maxvel, &
3954 "The maximum velocity allowed before the velocity "//&
3955 "components are truncated.", units=
"m s-1", default=3.0e8, scale=us%m_s_to_L_T, &
3956 do_not_log=.not.cs%clip_velocity)
3957 call get_param(param_file, mdl,
"MAXCFL_BT_CONT", cs%maxCFL_BT_cont, &
3958 "The maximum permitted CFL number associated with the "//&
3959 "barotropic accelerations from the summed velocities "//&
3960 "times the time-derivatives of thicknesses.", units=
"nondim", &
3962 call get_param(param_file, mdl,
"VEL_UNDERFLOW", cs%vel_underflow, &
3963 "A negligibly small velocity magnitude below which velocity "//&
3964 "components are set to 0. A reasonable value might be "//&
3965 "1e-30 m/s, which is less than an Angstrom divided by "//&
3966 "the age of the universe.", units=
"m s-1", default=0.0, scale=us%m_s_to_L_T)
3968 call get_param(param_file, mdl,
"DT_BT_FILTER", cs%dt_bt_filter, &
3969 "A time-scale over which the barotropic mode solutions "//&
3970 "are filtered, in seconds if positive, or as a fraction "//&
3971 "of DT if negative. When used this can never be taken to "//&
3972 "be longer than 2*dt. Set this to 0 to apply no filtering.", &
3973 units=
"sec or nondim", default=-0.25)
3974 if (cs%dt_bt_filter > 0.0) cs%dt_bt_filter = us%s_to_T*cs%dt_bt_filter
3975 call get_param(param_file, mdl,
"G_BT_EXTRA", cs%G_extra, &
3976 "A nondimensional factor by which gtot is enhanced.", &
3977 units=
"nondim", default=0.0)
3978 call get_param(param_file, mdl,
"SSH_EXTRA", ssh_extra, &
3979 "An estimate of how much higher SSH might get, for use "//&
3980 "in calculating the safe external wave speed. The "//&
3981 "default is the minimum of 10 m or 5% of MAXIMUM_DEPTH.", &
3982 units=
"m", default=min(10.0,0.05*g%max_depth*us%Z_to_m), scale=us%m_to_Z)
3984 call get_param(param_file, mdl,
"DEBUG", cs%debug, &
3985 "If true, write out verbose debugging data.", &
3986 default=.false., debuggingparam=.true.)
3987 call get_param(param_file, mdl,
"DEBUG_BT", cs%debug_bt, &
3988 "If true, write out verbose debugging data within the "//&
3989 "barotropic time-stepping loop. The data volume can be "//&
3990 "quite large if this is true.", default=cs%debug, &
3991 debuggingparam=.true.)
3993 cs%linearized_BT_PV = .true.
3994 call get_param(param_file, mdl,
"BEBT", cs%bebt, &
3995 "BEBT determines whether the barotropic time stepping "//&
3996 "uses the forward-backward time-stepping scheme or a "//&
3997 "backward Euler scheme. BEBT is valid in the range from "//&
3998 "0 (for a forward-backward treatment of nonrotating "//&
3999 "gravity waves) to 1 (for a backward Euler treatment). "//&
4000 "In practice, BEBT must be greater than about 0.05.", &
4001 units=
"nondim", default=0.1)
4002 call get_param(param_file, mdl,
"DTBT", dtbt_input, &
4003 "The barotropic time step, in s. DTBT is only used with "//&
4004 "the split explicit time stepping. To set the time step "//&
4005 "automatically based the maximum stable value use 0, or "//&
4006 "a negative value gives the fraction of the stable value. "//&
4007 "Setting DTBT to 0 is the same as setting it to -0.98. "//&
4008 "The value of DTBT that will actually be used is an "//&
4009 "integer fraction of DT, rounding down.", units=
"s or nondim",&
4011 call get_param(param_file, mdl,
"BT_USE_OLD_CORIOLIS_BRACKET_BUG", &
4012 cs%use_old_coriolis_bracket_bug , &
4013 "If True, use an order of operations that is not bitwise "//&
4014 "rotationally symmetric in the meridional Coriolis term of "//&
4015 "the barotropic solver.", default=.false.)
4018 call clone_mom_domain(g%Domain, cs%BT_Domain, min_halo=wd_halos, symmetric=.true.)
4019 #ifdef STATIC_MEMORY_
4020 if (wd_halos(1) /= whaloi_+nihalo_)
call mom_error(fatal,
"barotropic_init: "//&
4021 "Barotropic x-halo sizes are incorrectly resized with STATIC_MEMORY_.")
4022 if (wd_halos(2) /= whaloj_+njhalo_)
call mom_error(fatal,
"barotropic_init: "//&
4023 "Barotropic y-halo sizes are incorrectly resized with STATIC_MEMORY_.")
4025 if (bt_halo_sz > 0)
then
4026 if (wd_halos(1) > bt_halo_sz) &
4027 call mom_mesg(
"barotropic_init: barotropic x-halo size increased.", 3)
4028 if (wd_halos(2) > bt_halo_sz) &
4029 call mom_mesg(
"barotropic_init: barotropic y-halo size increased.", 3)
4033 cs%isdw = g%isc-wd_halos(1) ; cs%iedw = g%iec+wd_halos(1)
4034 cs%jsdw = g%jsc-wd_halos(2) ; cs%jedw = g%jec+wd_halos(2)
4035 isdw = cs%isdw ; iedw = cs%iedw ; jsdw = cs%jsdw ; jedw = cs%jedw
4037 alloc_(cs%frhatu(isdb:iedb,jsd:jed,nz)) ; alloc_(cs%frhatv(isd:ied,jsdb:jedb,nz))
4038 alloc_(cs%eta_cor(isd:ied,jsd:jed))
4039 if (cs%bound_BT_corr)
then
4040 alloc_(cs%eta_cor_bound(isd:ied,jsd:jed)) ; cs%eta_cor_bound(:,:) = 0.0
4042 alloc_(cs%IDatu(isdb:iedb,jsd:jed)) ; alloc_(cs%IDatv(isd:ied,jsdb:jedb))
4044 alloc_(cs%ua_polarity(isdw:iedw,jsdw:jedw))
4045 alloc_(cs%va_polarity(isdw:iedw,jsdw:jedw))
4047 cs%frhatu(:,:,:) = 0.0 ; cs%frhatv(:,:,:) = 0.0
4048 cs%eta_cor(:,:) = 0.0
4049 cs%IDatu(:,:) = 0.0 ; cs%IDatv(:,:) = 0.0
4051 cs%ua_polarity(:,:) = 1.0 ; cs%va_polarity(:,:) = 1.0
4052 call create_group_pass(pass_a_polarity, cs%ua_polarity, cs%va_polarity, cs%BT_domain, to_all, agrid)
4053 call do_group_pass(pass_a_polarity, cs%BT_domain)
4055 if (use_bt_cont_type) &
4059 allocate(cs%debug_BT_HI)
4060 cs%debug_BT_HI%isc=g%isc
4061 cs%debug_BT_HI%iec=g%iec
4062 cs%debug_BT_HI%jsc=g%jsc
4063 cs%debug_BT_HI%jec=g%jec
4064 cs%debug_BT_HI%IscB=g%isc-1
4065 cs%debug_BT_HI%IecB=g%iec
4066 cs%debug_BT_HI%JscB=g%jsc-1
4067 cs%debug_BT_HI%JecB=g%jec
4068 cs%debug_BT_HI%isd=cs%isdw
4069 cs%debug_BT_HI%ied=cs%iedw
4070 cs%debug_BT_HI%jsd=cs%jsdw
4071 cs%debug_BT_HI%jed=cs%jedw
4072 cs%debug_BT_HI%IsdB=cs%isdw-1
4073 cs%debug_BT_HI%IedB=cs%iedw
4074 cs%debug_BT_HI%JsdB=cs%jsdw-1
4075 cs%debug_BT_HI%JedB=cs%jedw
4079 alloc_(cs%IareaT(cs%isdw:cs%iedw,cs%jsdw:cs%jedw)) ; cs%IareaT(:,:) = 0.0
4080 alloc_(cs%bathyT(cs%isdw:cs%iedw,cs%jsdw:cs%jedw)) ; cs%bathyT(:,:) = gv%Angstrom_m
4081 alloc_(cs%IdxCu(cs%isdw-1:cs%iedw,cs%jsdw:cs%jedw)) ; cs%IdxCu(:,:) = 0.0
4082 alloc_(cs%IdyCv(cs%isdw:cs%iedw,cs%jsdw-1:cs%jedw)) ; cs%IdyCv(:,:) = 0.0
4083 alloc_(cs%dy_Cu(cs%isdw-1:cs%iedw,cs%jsdw:cs%jedw)) ; cs%dy_Cu(:,:) = 0.0
4084 alloc_(cs%dx_Cv(cs%isdw:cs%iedw,cs%jsdw-1:cs%jedw)) ; cs%dx_Cv(:,:) = 0.0
4085 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
4086 cs%IareaT(i,j) = g%IareaT(i,j)
4087 cs%bathyT(i,j) = g%bathyT(i,j)
4092 do j=g%jsd,g%jed ;
do i=g%IsdB,g%IedB
4093 cs%IdxCu(i,j) = g%IdxCu(i,j) ; cs%dy_Cu(i,j) = g%dy_Cu(i,j)
4095 do j=g%JsdB,g%JedB ;
do i=g%isd,g%ied
4096 cs%IdyCv(i,j) = g%IdyCv(i,j) ; cs%dx_Cv(i,j) = g%dx_Cv(i,j)
4100 call create_group_pass(pass_static_data, cs%IdxCu, cs%IdyCv, cs%BT_domain, to_all+scalar_pair)
4101 call create_group_pass(pass_static_data, cs%dy_Cu, cs%dx_Cv, cs%BT_domain, to_all+scalar_pair)
4102 call do_group_pass(pass_static_data, cs%BT_domain)
4104 if (cs%linearized_BT_PV)
then
4105 alloc_(cs%q_D(cs%isdw-1:cs%iedw,cs%jsdw-1:cs%jedw))
4106 alloc_(cs%D_u_Cor(cs%isdw-1:cs%iedw,cs%jsdw:cs%jedw))
4107 alloc_(cs%D_v_Cor(cs%isdw:cs%iedw,cs%jsdw-1:cs%jedw))
4108 cs%q_D(:,:) = 0.0 ; cs%D_u_Cor(:,:) = 0.0 ; cs%D_v_Cor(:,:) = 0.0
4109 do j=js,je ;
do i=is-1,ie
4110 cs%D_u_Cor(i,j) = 0.5 * (g%bathyT(i+1,j) + g%bathyT(i,j))
4112 do j=js-1,je ;
do i=is,ie
4113 cs%D_v_Cor(i,j) = 0.5 * (g%bathyT(i,j+1) + g%bathyT(i,j))
4115 do j=js-1,je ;
do i=is-1,ie
4116 if (g%mask2dT(i,j)+g%mask2dT(i,j+1)+g%mask2dT(i+1,j)+g%mask2dT(i+1,j+1)>0.)
then
4117 cs%q_D(i,j) = 0.25 * g%CoriolisBu(i,j) * &
4118 ((g%areaT(i,j) + g%areaT(i+1,j+1)) + (g%areaT(i+1,j) + g%areaT(i,j+1))) / &
4119 ((g%areaT(i,j) * g%bathyT(i,j) + g%areaT(i+1,j+1) * g%bathyT(i+1,j+1)) + &
4120 (g%areaT(i+1,j) * g%bathyT(i+1,j) + g%areaT(i,j+1) * g%bathyT(i,j+1)) )
4127 call create_group_pass(pass_q_d_cor, cs%q_D, cs%BT_Domain, to_all, position=corner)
4130 call do_group_pass(pass_q_d_cor, cs%BT_Domain)
4133 if (cs%linear_wave_drag)
then
4134 alloc_(cs%lin_drag_u(isdb:iedb,jsd:jed)) ; cs%lin_drag_u(:,:) = 0.0
4135 alloc_(cs%lin_drag_v(isd:ied,jsdb:jedb)) ; cs%lin_drag_v(:,:) = 0.0
4137 if (len_trim(wave_drag_file) > 0)
then
4138 inputdir =
"." ;
call get_param(param_file, mdl,
"INPUTDIR", inputdir)
4139 wave_drag_file = trim(slasher(inputdir))//trim(wave_drag_file)
4140 call log_param(param_file, mdl,
"INPUTDIR/BT_WAVE_DRAG_FILE", wave_drag_file)
4142 allocate(lin_drag_h(isd:ied,jsd:jed)) ; lin_drag_h(:,:) = 0.0
4144 call mom_read_data(wave_drag_file, wave_drag_var, lin_drag_h, g%Domain, scale=us%m_to_Z*us%T_to_s)
4145 call pass_var(lin_drag_h, g%Domain)
4146 do j=js,je ;
do i=is-1,ie
4147 cs%lin_drag_u(i,j) = (gv%Z_to_H * wave_drag_scale) * &
4148 0.5 * (lin_drag_h(i,j) + lin_drag_h(i+1,j))
4150 do j=js-1,je ;
do i=is,ie
4151 cs%lin_drag_v(i,j) = (gv%Z_to_H * wave_drag_scale) * &
4152 0.5 * (lin_drag_h(i,j) + lin_drag_h(i,j+1))
4154 deallocate(lin_drag_h)
4158 cs%dtbt_fraction = 0.98 ;
if (dtbt_input < 0.0) cs%dtbt_fraction = -dtbt_input
4163 if ((us%s_to_T_restart /= 0.0) .and. (us%s_to_T_restart /= us%s_to_T)) &
4164 dtbt_tmp = (us%s_to_T / us%s_to_T_restart) * cs%dtbt
4169 do k=1,g%ke ; gtot_estimate = gtot_estimate + gv%g_prime(k) ;
enddo
4170 call set_dtbt(g, gv, us, cs, gtot_est=gtot_estimate, ssh_add=ssh_extra)
4172 if (dtbt_input > 0.0)
then
4173 cs%dtbt = us%s_to_T * dtbt_input
4174 elseif (dtbt_tmp > 0.0)
then
4177 if ((dtbt_tmp > 0.0) .and. (dtbt_input > 0.0)) calc_dtbt = .false.
4179 call log_param(param_file, mdl,
"DTBT as used", cs%dtbt*us%T_to_s)
4180 call log_param(param_file, mdl,
"estimated maximum DTBT", cs%dtbt_max*us%T_to_s)
4185 if (gv%Boussinesq)
then
4186 thickness_units =
"m" ; flux_units =
"m3 s-1"
4188 thickness_units =
"kg m-2" ; flux_units =
"kg s-1"
4191 cs%id_PFu_bt = register_diag_field(
'ocean_model',
'PFuBT', diag%axesCu1, time, &
4192 'Zonal Anomalous Barotropic Pressure Force Force Acceleration',
'm s-2', conversion=us%L_T2_to_m_s2)
4193 cs%id_PFv_bt = register_diag_field(
'ocean_model',
'PFvBT', diag%axesCv1, time, &
4194 'Meridional Anomalous Barotropic Pressure Force Acceleration',
'm s-2', conversion=us%L_T2_to_m_s2)
4195 cs%id_Coru_bt = register_diag_field(
'ocean_model',
'CoruBT', diag%axesCu1, time, &
4196 'Zonal Barotropic Coriolis Acceleration',
'm s-2', conversion=us%L_T2_to_m_s2)
4197 cs%id_Corv_bt = register_diag_field(
'ocean_model',
'CorvBT', diag%axesCv1, time, &
4198 'Meridional Barotropic Coriolis Acceleration',
'm s-2', conversion=us%L_T2_to_m_s2)
4199 cs%id_uaccel = register_diag_field(
'ocean_model',
'u_accel_bt', diag%axesCu1, time, &
4200 'Barotropic zonal acceleration',
'm s-2', conversion=us%L_T2_to_m_s2)
4201 cs%id_vaccel = register_diag_field(
'ocean_model',
'v_accel_bt', diag%axesCv1, time, &
4202 'Barotropic meridional acceleration',
'm s-2', conversion=us%L_T2_to_m_s2)
4203 cs%id_ubtforce = register_diag_field(
'ocean_model',
'ubtforce', diag%axesCu1, time, &
4204 'Barotropic zonal acceleration from baroclinic terms',
'm s-2', conversion=us%L_T2_to_m_s2)
4205 cs%id_vbtforce = register_diag_field(
'ocean_model',
'vbtforce', diag%axesCv1, time, &
4206 'Barotropic meridional acceleration from baroclinic terms',
'm s-2', conversion=us%L_T2_to_m_s2)
4208 cs%id_eta_bt = register_diag_field(
'ocean_model',
'eta_bt', diag%axesT1, time, &
4209 'Barotropic end SSH', thickness_units, conversion=gv%H_to_m)
4210 cs%id_ubt = register_diag_field(
'ocean_model',
'ubt', diag%axesCu1, time, &
4211 'Barotropic end zonal velocity',
'm s-1', conversion=us%L_T_to_m_s)
4212 cs%id_vbt = register_diag_field(
'ocean_model',
'vbt', diag%axesCv1, time, &
4213 'Barotropic end meridional velocity',
'm s-1', conversion=us%L_T_to_m_s)
4214 cs%id_eta_st = register_diag_field(
'ocean_model',
'eta_st', diag%axesT1, time, &
4215 'Barotropic start SSH', thickness_units, conversion=gv%H_to_m)
4216 cs%id_ubt_st = register_diag_field(
'ocean_model',
'ubt_st', diag%axesCu1, time, &
4217 'Barotropic start zonal velocity',
'm s-1', conversion=us%L_T_to_m_s)
4218 cs%id_vbt_st = register_diag_field(
'ocean_model',
'vbt_st', diag%axesCv1, time, &
4219 'Barotropic start meridional velocity',
'm s-1', conversion=us%L_T_to_m_s)
4220 cs%id_ubtav = register_diag_field(
'ocean_model',
'ubtav', diag%axesCu1, time, &
4221 'Barotropic time-average zonal velocity',
'm s-1', conversion=us%L_T_to_m_s)
4222 cs%id_vbtav = register_diag_field(
'ocean_model',
'vbtav', diag%axesCv1, time, &
4223 'Barotropic time-average meridional velocity',
'm s-1', conversion=us%L_T_to_m_s)
4224 cs%id_eta_cor = register_diag_field(
'ocean_model',
'eta_cor', diag%axesT1, time, &
4225 'Corrective mass flux',
'm s-1', conversion=gv%H_to_m)
4226 cs%id_visc_rem_u = register_diag_field(
'ocean_model',
'visc_rem_u', diag%axesCuL, time, &
4227 'Viscous remnant at u',
'nondim')
4228 cs%id_visc_rem_v = register_diag_field(
'ocean_model',
'visc_rem_v', diag%axesCvL, time, &
4229 'Viscous remnant at v',
'nondim')
4230 cs%id_gtotn = register_diag_field(
'ocean_model',
'gtot_n', diag%axesT1, time, &
4231 'gtot to North',
'm s-2', conversion=gv%m_to_H*(us%L_T_to_m_s**2))
4232 cs%id_gtots = register_diag_field(
'ocean_model',
'gtot_s', diag%axesT1, time, &
4233 'gtot to South',
'm s-2', conversion=gv%m_to_H*(us%L_T_to_m_s**2))
4234 cs%id_gtote = register_diag_field(
'ocean_model',
'gtot_e', diag%axesT1, time, &
4235 'gtot to East',
'm s-2', conversion=gv%m_to_H*(us%L_T_to_m_s**2))
4236 cs%id_gtotw = register_diag_field(
'ocean_model',
'gtot_w', diag%axesT1, time, &
4237 'gtot to West',
'm s-2', conversion=gv%m_to_H*(us%L_T_to_m_s**2))
4238 cs%id_eta_hifreq = register_diag_field(
'ocean_model',
'eta_hifreq', diag%axesT1, time, &
4239 'High Frequency Barotropic SSH', thickness_units, conversion=gv%H_to_m)
4240 cs%id_ubt_hifreq = register_diag_field(
'ocean_model',
'ubt_hifreq', diag%axesCu1, time, &
4241 'High Frequency Barotropic zonal velocity',
'm s-1', conversion=us%L_T_to_m_s)
4242 cs%id_vbt_hifreq = register_diag_field(
'ocean_model',
'vbt_hifreq', diag%axesCv1, time, &
4243 'High Frequency Barotropic meridional velocity',
'm s-1', conversion=us%L_T_to_m_s)
4244 cs%id_eta_pred_hifreq = register_diag_field(
'ocean_model',
'eta_pred_hifreq', diag%axesT1, time, &
4245 'High Frequency Predictor Barotropic SSH', thickness_units, &
4246 conversion=gv%H_to_m)
4247 cs%id_uhbt_hifreq = register_diag_field(
'ocean_model',
'uhbt_hifreq', diag%axesCu1, time, &
4248 'High Frequency Barotropic zonal transport',
'm3 s-1', &
4249 conversion=gv%H_to_m*us%L_to_m*us%L_T_to_m_s)
4250 cs%id_vhbt_hifreq = register_diag_field(
'ocean_model',
'vhbt_hifreq', diag%axesCv1, time, &
4251 'High Frequency Barotropic meridional transport',
'm3 s-1', &
4252 conversion=gv%H_to_m*us%L_to_m*us%L_T_to_m_s)
4253 cs%id_frhatu = register_diag_field(
'ocean_model',
'frhatu', diag%axesCuL, time, &
4254 'Fractional thickness of layers in u-columns',
'nondim')
4255 cs%id_frhatv = register_diag_field(
'ocean_model',
'frhatv', diag%axesCvL, time, &
4256 'Fractional thickness of layers in v-columns',
'nondim')
4257 cs%id_frhatu1 = register_diag_field(
'ocean_model',
'frhatu1', diag%axesCuL, time, &
4258 'Predictor Fractional thickness of layers in u-columns',
'nondim')
4259 cs%id_frhatv1 = register_diag_field(
'ocean_model',
'frhatv1', diag%axesCvL, time, &
4260 'Predictor Fractional thickness of layers in v-columns',
'nondim')
4261 cs%id_uhbt = register_diag_field(
'ocean_model',
'uhbt', diag%axesCu1, time, &
4262 'Barotropic zonal transport averaged over a baroclinic step',
'm3 s-1', &
4263 conversion=gv%H_to_m*us%L_to_m*us%L_T_to_m_s)
4264 cs%id_vhbt = register_diag_field(
'ocean_model',
'vhbt', diag%axesCv1, time, &
4265 'Barotropic meridional transport averaged over a baroclinic step',
'm3 s-1', &
4266 conversion=gv%H_to_m*us%L_to_m*us%L_T_to_m_s)
4268 if (use_bt_cont_type)
then
4269 cs%id_BTC_FA_u_EE = register_diag_field(
'ocean_model',
'BTC_FA_u_EE', diag%axesCu1, time, &
4270 'BTCont type far east face area',
'm2', conversion=us%L_to_m*gv%H_to_m)
4271 cs%id_BTC_FA_u_E0 = register_diag_field(
'ocean_model',
'BTC_FA_u_E0', diag%axesCu1, time, &
4272 'BTCont type near east face area',
'm2', conversion=us%L_to_m*gv%H_to_m)
4273 cs%id_BTC_FA_u_WW = register_diag_field(
'ocean_model',
'BTC_FA_u_WW', diag%axesCu1, time, &
4274 'BTCont type far west face area',
'm2', conversion=us%L_to_m*gv%H_to_m)
4275 cs%id_BTC_FA_u_W0 = register_diag_field(
'ocean_model',
'BTC_FA_u_W0', diag%axesCu1, time, &
4276 'BTCont type near west face area',
'm2', conversion=us%L_to_m*gv%H_to_m)
4277 cs%id_BTC_ubt_EE = register_diag_field(
'ocean_model',
'BTC_ubt_EE', diag%axesCu1, time, &
4278 'BTCont type far east velocity',
'm s-1', conversion=us%L_T_to_m_s)
4279 cs%id_BTC_ubt_WW = register_diag_field(
'ocean_model',
'BTC_ubt_WW', diag%axesCu1, time, &
4280 'BTCont type far west velocity',
'm s-1', conversion=us%L_T_to_m_s)
4281 cs%id_BTC_FA_v_NN = register_diag_field(
'ocean_model',
'BTC_FA_v_NN', diag%axesCv1, time, &
4282 'BTCont type far north face area',
'm2', conversion=us%L_to_m*gv%H_to_m)
4283 cs%id_BTC_FA_v_N0 = register_diag_field(
'ocean_model',
'BTC_FA_v_N0', diag%axesCv1, time, &
4284 'BTCont type near north face area',
'm2', conversion=us%L_to_m*gv%H_to_m)
4285 cs%id_BTC_FA_v_SS = register_diag_field(
'ocean_model',
'BTC_FA_v_SS', diag%axesCv1, time, &
4286 'BTCont type far south face area',
'm2', conversion=us%L_to_m*gv%H_to_m)
4287 cs%id_BTC_FA_v_S0 = register_diag_field(
'ocean_model',
'BTC_FA_v_S0', diag%axesCv1, time, &
4288 'BTCont type near south face area',
'm2', conversion=us%L_to_m*gv%H_to_m)
4289 cs%id_BTC_vbt_NN = register_diag_field(
'ocean_model',
'BTC_vbt_NN', diag%axesCv1, time, &
4290 'BTCont type far north velocity',
'm s-1', conversion=us%L_T_to_m_s)
4291 cs%id_BTC_vbt_SS = register_diag_field(
'ocean_model',
'BTC_vbt_SS', diag%axesCv1, time, &
4292 'BTCont type far south velocity',
'm s-1', conversion=us%L_T_to_m_s)
4294 cs%id_uhbt0 = register_diag_field(
'ocean_model',
'uhbt0', diag%axesCu1, time, &
4295 'Barotropic zonal transport difference',
'm3 s-1', conversion=gv%H_to_m*us%L_to_m**2*us%s_to_T)
4296 cs%id_vhbt0 = register_diag_field(
'ocean_model',
'vhbt0', diag%axesCv1, time, &
4297 'Barotropic meridional transport difference',
'm3 s-1', conversion=gv%H_to_m*us%L_to_m**2*us%s_to_T)
4299 if (cs%id_frhatu1 > 0)
call safe_alloc_ptr(cs%frhatu1, isdb,iedb,jsd,jed,nz)
4300 if (cs%id_frhatv1 > 0)
call safe_alloc_ptr(cs%frhatv1, isd,ied,jsdb,jedb,nz)
4304 call btcalc(h, g, gv, cs, may_use_default=.true.)
4305 cs%ubtav(:,:) = 0.0 ; cs%vbtav(:,:) = 0.0
4306 do k=1,nz ;
do j=js,je ;
do i=is-1,ie
4307 cs%ubtav(i,j) = cs%ubtav(i,j) + cs%frhatu(i,j,k) * u(i,j,k)
4308 enddo ;
enddo ;
enddo
4309 do k=1,nz ;
do j=js-1,je ;
do i=is,ie
4310 cs%vbtav(i,j) = cs%vbtav(i,j) + cs%frhatv(i,j,k) * v(i,j,k)
4311 enddo ;
enddo ;
enddo
4312 elseif ((us%s_to_T_restart*us%m_to_L_restart /= 0.0) .and. &
4313 (us%m_to_L*us%s_to_T_restart) /= (us%m_to_L_restart*us%s_to_T))
then
4314 vel_rescale = (us%m_to_L*us%s_to_T_restart) / (us%m_to_L_restart*us%s_to_T)
4315 do j=js,je ;
do i=is-1,ie ; cs%ubtav(i,j) = vel_rescale * cs%ubtav(i,j) ;
enddo ;
enddo
4316 do j=js-1,je ;
do i=is,ie ; cs%vbtav(i,j) = vel_rescale * cs%vbtav(i,j) ;
enddo ;
enddo
4321 do j=js,je ;
do i=is-1,ie ; cs%ubt_IC(i,j) = cs%ubtav(i,j) ;
enddo ;
enddo
4322 do j=js-1,je ;
do i=is,ie ; cs%vbt_IC(i,j) = cs%vbtav(i,j) ;
enddo ;
enddo
4323 elseif ((us%s_to_T_restart*us%m_to_L_restart /= 0.0) .and. &
4324 (us%m_to_L*us%s_to_T_restart) /= (us%m_to_L_restart*us%s_to_T))
then
4325 vel_rescale = (us%m_to_L*us%s_to_T_restart) / (us%m_to_L_restart*us%s_to_T)
4326 do j=js,je ;
do i=is-1,ie ; cs%ubt_IC(i,j) = vel_rescale * cs%ubt_IC(i,j) ;
enddo ;
enddo
4327 do j=js-1,je ;
do i=is,ie ; cs%vbt_IC(i,j) = vel_rescale * cs%vbt_IC(i,j) ;
enddo ;
enddo
4332 do j=js,je ;
do i=is-1,ie
4333 if (g%mask2dCu(i,j)>0.)
then
4334 cs%IDatu(i,j) = g%mask2dCu(i,j) * 2.0 / (g%bathyT(i+1,j) + g%bathyT(i,j))
4339 do j=js-1,je ;
do i=is,ie
4340 if (g%mask2dCv(i,j)>0.)
then
4341 cs%IDatv(i,j) = g%mask2dCv(i,j) * 2.0 / (g%bathyT(i,j+1) + g%bathyT(i,j))
4348 if (cs%bound_BT_corr)
then
4351 do j=js,je ;
do i=is,ie
4352 cs%eta_cor_bound(i,j) = gv%m_to_H * g%IareaT(i,j) * 0.1 * cs%maxvel * &
4353 ((datu(i-1,j) + datu(i,j)) + (datv(i,j) + datv(i,j-1)))
4359 do j=js,je ;
do i=is-1,ie ; cs%uhbt_IC(i,j) = cs%ubtav(i,j) * datu(i,j) ;
enddo ;
enddo
4360 do j=js-1,je ;
do i=is,ie ; cs%vhbt_IC(i,j) = cs%vbtav(i,j) * datv(i,j) ;
enddo ;
enddo
4361 elseif ((us%s_to_T_restart * us%m_to_L_restart * gv%m_to_H_restart /= 0.0) .and. &
4362 ((us%s_to_T_restart * us%m_to_L**2 * gv%m_to_H) /= &
4363 (us%s_to_T * us%m_to_L_restart**2 * gv%m_to_H_restart)))
then
4364 uh_rescale = (us%s_to_T_restart * us%m_to_L**2 * gv%m_to_H) / &
4365 (us%s_to_T * us%m_to_L_restart**2 * gv%m_to_H_restart)
4366 do j=js,je ;
do i=is-1,ie ; cs%uhbt_IC(i,j) = uh_rescale * cs%uhbt_IC(i,j) ;
enddo ;
enddo
4367 do j=js-1,je ;
do i=is,ie ; cs%vhbt_IC(i,j) = uh_rescale * cs%vhbt_IC(i,j) ;
enddo ;
enddo
4373 call do_group_pass(pass_bt_hbt_btav, g%Domain)
4376 id_clock_calc_pre = cpu_clock_id(
'(Ocean BT pre-calcs only)', grain=clock_routine)
4377 id_clock_pass_pre = cpu_clock_id(
'(Ocean BT pre-step halo updates)', grain=clock_routine)
4378 id_clock_calc = cpu_clock_id(
'(Ocean BT stepping calcs only)', grain=clock_routine)
4379 id_clock_pass_step = cpu_clock_id(
'(Ocean BT stepping halo updates)', grain=clock_routine)
4381 id_clock_pass_post = cpu_clock_id(
'(Ocean BT post-step halo updates)', grain=clock_routine)
4382 if (dtbt_input <= 0.0) &
4383 id_clock_sync = cpu_clock_id(
'(Ocean BT global synch)', grain=clock_routine)
4391 real,
dimension(SZIB_(G),SZJ_(G)),
intent(inout) :: ubtav
4393 real,
dimension(SZI_(G),SZJB_(G)),
intent(inout) :: vbtav
4399 do j=g%jsc,g%jec ;
do i=g%isc-1,g%iec
4400 ubtav(i,j) = cs%ubtav(i,j)
4403 do j=g%jsc-1,g%jec ;
do i=g%isc,g%iec
4404 vbtav(i,j) = cs%vbtav(i,j)
4413 dealloc_(cs%frhatu) ; dealloc_(cs%frhatv)
4414 dealloc_(cs%IDatu) ; dealloc_(cs%IDatv)
4415 dealloc_(cs%ubtav) ; dealloc_(cs%vbtav)
4416 dealloc_(cs%eta_cor)
4417 dealloc_(cs%ua_polarity) ; dealloc_(cs%va_polarity)
4418 if (cs%bound_BT_corr)
then
4419 dealloc_(cs%eta_cor_bound)
4440 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
4441 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed
4442 isdb = hi%IsdB ; iedb = hi%IedB ; jsdb = hi%JsdB ; jedb = hi%JedB
4444 if (
associated(cs))
then
4445 call mom_error(warning,
"register_barotropic_restarts called with an associated "// &
4446 "control structure.")
4451 alloc_(cs%ubtav(isdb:iedb,jsd:jed)) ; cs%ubtav(:,:) = 0.0
4452 alloc_(cs%vbtav(isd:ied,jsdb:jedb)) ; cs%vbtav(:,:) = 0.0
4453 alloc_(cs%ubt_IC(isdb:iedb,jsd:jed)) ; cs%ubt_IC(:,:) = 0.0
4454 alloc_(cs%vbt_IC(isd:ied,jsdb:jedb)) ; cs%vbt_IC(:,:) = 0.0
4455 alloc_(cs%uhbt_IC(isdb:iedb,jsd:jed)) ; cs%uhbt_IC(:,:) = 0.0
4456 alloc_(cs%vhbt_IC(isd:ied,jsdb:jedb)) ; cs%vhbt_IC(:,:) = 0.0
4458 vd(2) =
var_desc(
"ubtav",
"m s-1",
"Time mean barotropic zonal velocity", &
4459 hor_grid=
'u', z_grid=
'1')
4460 vd(3) =
var_desc(
"vbtav",
"m s-1",
"Time mean barotropic meridional velocity",&
4461 hor_grid=
'v', z_grid=
'1')
4465 vd(2) =
var_desc(
"ubt_IC",
"m s-1", &
4466 longname=
"Next initial condition for the barotropic zonal velocity", &
4467 hor_grid=
'u', z_grid=
'1')
4468 vd(3) =
var_desc(
"vbt_IC",
"m s-1", &
4469 longname=
"Next initial condition for the barotropic meridional velocity",&
4470 hor_grid=
'v', z_grid=
'1')
4474 if (gv%Boussinesq)
then
4475 vd(2) =
var_desc(
"uhbt_IC",
"m3 s-1", &
4476 longname=
"Next initial condition for the barotropic zonal transport", &
4477 hor_grid=
'u', z_grid=
'1')
4478 vd(3) =
var_desc(
"vhbt_IC",
"m3 s-1", &
4479 longname=
"Next initial condition for the barotropic meridional transport",&
4480 hor_grid=
'v', z_grid=
'1')
4482 vd(2) =
var_desc(
"uhbt_IC",
"kg s-1", &
4483 longname=
"Next initial condition for the barotropic zonal transport", &
4484 hor_grid=
'u', z_grid=
'1')
4485 vd(3) =
var_desc(
"vhbt_IC",
"kg s-1", &
4486 longname=
"Next initial condition for the barotropic meridional transport",&
4487 hor_grid=
'v', z_grid=
'1')
4493 longname=
"Barotropic timestep", units=
"seconds")