This subroutine time steps the barotropic equations explicitly. For gravity waves, anything between a forwards-backwards scheme and a simulated backwards Euler scheme is used, with bebt between 0.0 and 1.0 determining the scheme. In practice, bebt must be of order 0.2 or greater. A forwards-backwards treatment of the Coriolis terms is always used.
390 type(ocean_grid_type),
intent(inout) :: G
391 type(verticalGrid_type),
intent(in) :: GV
392 type(unit_scale_type),
intent(in) :: US
393 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: U_in
394 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: V_in
395 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: eta_in
397 real,
intent(in) :: dt
398 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: bc_accel_u
399 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: bc_accel_v
401 type(mech_forcing),
intent(in) :: forces
402 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: pbce
405 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: eta_PF_in
411 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: U_Cor
413 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: V_Cor
414 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(out) :: accel_layer_u
416 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(out) :: accel_layer_v
418 real,
dimension(SZI_(G),SZJ_(G)),
intent(out) :: eta_out
420 real,
dimension(SZIB_(G),SZJ_(G)),
intent(out) :: uhbtav
423 real,
dimension(SZI_(G),SZJB_(G)),
intent(out) :: vhbtav
426 type(barotropic_CS),
pointer :: CS
428 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: visc_rem_u
434 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: visc_rem_v
435 real,
dimension(SZI_(G),SZJ_(G)),
optional,
intent(out) :: etaav
437 type(ocean_OBC_type),
optional,
pointer :: OBC
438 type(BT_cont_type),
optional,
pointer :: BT_cont
441 real,
dimension(:,:),
optional,
pointer :: eta_PF_start
444 real,
dimension(:,:),
optional,
pointer :: taux_bot
446 real,
dimension(:,:),
optional,
pointer :: tauy_bot
448 real,
dimension(:,:,:),
optional,
pointer :: uh0
450 real,
dimension(:,:,:),
optional,
pointer :: u_uh0
451 real,
dimension(:,:,:),
optional,
pointer :: vh0
453 real,
dimension(:,:,:),
optional,
pointer :: v_vh0
456 real :: ubt_Cor(SZIB_(G),SZJ_(G))
457 real :: vbt_Cor(SZI_(G),SZJB_(G))
459 real :: wt_u(SZIB_(G),SZJ_(G),SZK_(G))
460 real :: wt_v(SZI_(G),SZJB_(G),SZK_(G))
463 real,
dimension(SZIB_(G),SZJ_(G)) :: &
466 real,
dimension(SZI_(G),SZJB_(G)) :: &
469 real,
dimension(SZI_(G),SZJ_(G)) :: &
475 real :: q(SZIBW_(CS),SZJBW_(CS))
476 real,
dimension(SZIBW_(CS),SZJW_(CS)) :: &
509 real,
dimension(SZIW_(CS),SZJBW_(CS)) :: &
540 real,
target,
dimension(SZIW_(CS),SZJW_(CS)) :: &
544 real,
dimension(:,:),
pointer :: &
547 real,
dimension(SZIW_(CS),SZJW_(CS)) :: &
568 type(local_BT_cont_u_type),
dimension(SZIBW_(CS),SZJW_(CS)) :: &
570 type(local_BT_cont_v_type),
dimension(SZIW_(CS),SZJBW_(CS)) :: &
574 real,
dimension(SZIBW_(CS),SZJW_(CS)) :: &
575 ubt_prev, uhbt_prev, ubt_sum_prev, uhbt_sum_prev, ubt_wtd_prev
576 real,
dimension(SZIW_(CS),SZJBW_(CS)) :: &
577 vbt_prev, vhbt_prev, vbt_sum_prev, vhbt_sum_prev, vbt_wtd_prev
605 logical :: do_hifreq_output
606 logical :: use_BT_cont, do_ave, find_etaav, find_PF, find_Cor
607 logical :: ice_is_rigid, nonblock_setup, interp_eta_PF
608 logical :: project_velocity, add_uh0
612 real :: ice_strength = 0.0
620 real :: u_max_cor, v_max_cor
623 real :: accel_underflow
625 real,
allocatable,
dimension(:) :: wt_vel, wt_eta, wt_accel, wt_trans, wt_accel2
626 real :: sum_wt_vel, sum_wt_eta, sum_wt_accel, sum_wt_trans
627 real :: I_sum_wt_vel, I_sum_wt_eta, I_sum_wt_accel, I_sum_wt_trans
629 real :: trans_wt1, trans_wt2
632 logical :: apply_OBCs, apply_OBC_flather, apply_OBC_open
633 type(memory_size_type) :: MS
634 character(len=200) :: mesg
635 integer :: isv, iev, jsv, jev
637 integer :: isvf, ievf, jsvf, jevf, num_cycles
638 integer :: i, j, k, n
639 integer :: is, ie, js, je, nz, Isq, Ieq, Jsq, Jeq
640 integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
641 integer :: ioff, joff
643 if (.not.
associated(cs))
call mom_error(fatal, &
644 "btstep: Module MOM_barotropic must be initialized before it is used.")
645 if (.not.cs%split)
return
646 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
647 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
648 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
649 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
650 ms%isdw = cs%isdw ; ms%iedw = cs%iedw ; ms%jsdw = cs%jsdw ; ms%jedw = cs%jedw
651 dt_in_t = us%s_to_T*dt
653 accel_underflow = us%L_T_to_m_s*cs%vel_underflow * us%s_to_T*idt
655 use_bt_cont = .false.
656 if (
present(bt_cont)) use_bt_cont = (
associated(bt_cont))
658 interp_eta_pf = .false.
659 if (
present(eta_pf_start)) interp_eta_pf = (
associated(eta_pf_start))
661 project_velocity = cs%BT_project_velocity
665 if ((.not.use_bt_cont) .and. cs%Nonlinear_continuity .and. &
666 (cs%Nonlin_cont_update_period > 0)) stencil = 2
668 do_ave = query_averaging_enabled(cs%diag)
669 find_etaav =
present(etaav)
670 find_pf = (do_ave .and. ((cs%id_PFu_bt > 0) .or. (cs%id_PFv_bt > 0)))
671 find_cor = (do_ave .and. ((cs%id_Coru_bt > 0) .or. (cs%id_Corv_bt > 0)))
674 if (
present(uh0)) add_uh0 =
associated(uh0)
675 if (add_uh0 .and. .not.(
present(vh0) .and.
present(u_uh0) .and. &
676 present(v_vh0)))
call mom_error(fatal, &
677 "btstep: vh0, u_uh0, and v_vh0 must be present if uh0 is used.")
678 if (add_uh0 .and. .not.(
associated(vh0) .and.
associated(u_uh0) .and. &
679 associated(v_vh0)))
call mom_error(fatal, &
680 "btstep: vh0, u_uh0, and v_vh0 must be associated if uh0 is used.")
683 nonblock_setup = g%nonblocking_updates
685 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
687 apply_obcs = .false. ; cs%BT_OBC%apply_u_OBCs = .false. ; cs%BT_OBC%apply_v_OBCs = .false.
688 apply_obc_open = .false.
689 apply_obc_flather = .false.
690 if (
present(obc))
then ;
if (
associated(obc))
then
691 cs%BT_OBC%apply_u_OBCs = obc%open_u_BCs_exist_globally .or. obc%specified_u_BCs_exist_globally
692 cs%BT_OBC%apply_v_OBCs = obc%open_v_BCs_exist_globally .or. obc%specified_v_BCs_exist_globally
693 apply_obc_flather = open_boundary_query(obc, apply_flather_obc=.true.)
694 apply_obc_open = open_boundary_query(obc, apply_open_obc=.true.)
695 apply_obcs = open_boundary_query(obc, apply_specified_obc=.true.) .or. &
696 apply_obc_flather .or. apply_obc_open
698 if (apply_obc_flather .and. .not.gv%Boussinesq)
call mom_error(fatal, &
699 "btstep: Flather open boundary conditions have not yet been "// &
700 "implemented for a non-Boussinesq model.")
704 if (cs%use_wide_halos) &
705 num_cycles = min((is-cs%isdw) / stencil, (js-cs%jsdw) / stencil)
706 isvf = is - (num_cycles-1)*stencil ; ievf = ie + (num_cycles-1)*stencil
707 jsvf = js - (num_cycles-1)*stencil ; jevf = je + (num_cycles-1)*stencil
709 nstep = ceiling(dt/cs%dtbt - 0.0001)
710 if (is_root_pe() .and. (nstep /= cs%nstep_last))
then
711 write(mesg,
'("btstep is using a dynamic barotropic timestep of ", ES12.6, &
712 & " seconds, max ", ES12.6, ".")') (dt/nstep), cs%dtbt_max
713 call mom_mesg(mesg, 3)
715 cs%nstep_last = nstep
718 instep = 1.0 / real(nstep)
719 dtbt = dt_in_t * instep
722 mass_to_z = us%m_to_L*us%T_to_s**2 * us%m_to_Z / gv%Rho0
725 if (project_velocity)
then
726 trans_wt1 = (1.0 + be_proj); trans_wt2 = -be_proj
728 trans_wt1 = bebt ; trans_wt2 = (1.0-bebt)
731 do_hifreq_output = .false.
732 if ((cs%id_ubt_hifreq > 0) .or. (cs%id_vbt_hifreq > 0) .or. &
733 (cs%id_eta_hifreq > 0) .or. (cs%id_eta_pred_hifreq > 0) .or. &
734 (cs%id_uhbt_hifreq > 0) .or. (cs%id_vhbt_hifreq > 0))
then
735 do_hifreq_output = query_averaging_enabled(cs%diag, time_int_in, time_end_in)
736 if (do_hifreq_output) &
737 time_bt_start = time_end_in - real_to_time(dt)
741 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
742 if (.not. cs%linearized_BT_PV)
then
743 call create_group_pass(cs%pass_q_DCor, q, cs%BT_Domain, to_all, position=corner)
744 call create_group_pass(cs%pass_q_DCor, dcor_u, dcor_v, cs%BT_Domain, &
747 if ((isq > is-1) .or. (jsq > js-1)) &
748 call create_group_pass(cs%pass_tmp_uv, tmp_u, tmp_v, g%Domain)
749 call create_group_pass(cs%pass_gtot, gtot_e, gtot_n, cs%BT_Domain, &
750 to_all+scalar_pair, agrid)
751 call create_group_pass(cs%pass_gtot, gtot_w, gtot_s, cs%BT_Domain, &
752 to_all+scalar_pair, agrid)
754 if (cs%dynamic_psurf) &
755 call create_group_pass(cs%pass_eta_bt_rem, dyn_coef_eta, cs%BT_Domain)
756 if (interp_eta_pf)
then
757 call create_group_pass(cs%pass_eta_bt_rem, eta_pf_1, cs%BT_Domain)
758 call create_group_pass(cs%pass_eta_bt_rem, d_eta_pf, cs%BT_Domain)
760 call create_group_pass(cs%pass_eta_bt_rem, eta_pf, cs%BT_Domain)
762 call create_group_pass(cs%pass_eta_bt_rem, eta_src, cs%BT_Domain)
766 call create_group_pass(cs%pass_eta_bt_rem, bt_rem_u, bt_rem_v, &
767 cs%BT_Domain, to_all+scalar_pair)
768 if (cs%linear_wave_drag) &
769 call create_group_pass(cs%pass_eta_bt_rem, rayleigh_u, rayleigh_v, &
770 cs%BT_Domain, to_all+scalar_pair)
773 if (((g%isd > cs%isdw) .or. (g%jsd > cs%jsdw)) .or. (isq <= is-1) .or. (jsq <= js-1)) &
774 call create_group_pass(cs%pass_force_hbt0_Cor_ref, bt_force_u, bt_force_v, cs%BT_Domain)
775 if (add_uh0)
call create_group_pass(cs%pass_force_hbt0_Cor_ref, uhbt0, vhbt0, cs%BT_Domain)
776 call create_group_pass(cs%pass_force_hbt0_Cor_ref, cor_ref_u, cor_ref_v, cs%BT_Domain)
777 if (.not. use_bt_cont)
then
778 call create_group_pass(cs%pass_Dat_uv, datu, datv, cs%BT_Domain, to_all+scalar_pair)
780 call create_group_pass(cs%pass_eta_ubt, eta, cs%BT_Domain)
781 call create_group_pass(cs%pass_eta_ubt, ubt, vbt, cs%BT_Domain)
783 call create_group_pass(cs%pass_ubt_Cor, ubt_cor, vbt_cor, g%Domain)
787 call create_group_pass(cs%pass_etaav, etaav, g%Domain)
789 call create_group_pass(cs%pass_e_anom, e_anom, g%Domain)
790 call create_group_pass(cs%pass_ubta_uhbta, cs%ubtav, cs%vbtav, g%Domain)
791 call create_group_pass(cs%pass_ubta_uhbta, uhbtav, vhbtav, g%Domain)
793 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
799 if (cs%linearized_BT_PV)
then
801 do j=jsvf-2,jevf+1 ;
do i=isvf-2,ievf+1
805 do j=jsvf-1,jevf+1 ;
do i=isvf-2,ievf+1
806 dcor_u(i,j) = cs%D_u_Cor(i,j)
809 do j=jsvf-2,jevf+1 ;
do i=isvf-1,ievf+1
810 dcor_v(i,j) = cs%D_v_Cor(i,j)
813 q(:,:) = 0.0 ; dcor_u(:,:) = 0.0 ; dcor_v(:,:) = 0.0
817 do j=js,je ;
do i=is-1,ie
818 dcor_u(i,j) = 0.5 * (g%bathyT(i+1,j) + g%bathyT(i,j))
821 do j=js-1,je ;
do i=is,ie
822 dcor_v(i,j) = 0.5 * (g%bathyT(i,j+1) + g%bathyT(i,j))
825 do j=js-1,je ;
do i=is-1,ie
826 q(i,j) = 0.25 * g%CoriolisBu(i,j) * &
827 ((g%areaT(i,j) + g%areaT(i+1,j+1)) + (g%areaT(i+1,j) + g%areaT(i,j+1))) / &
828 ((g%areaT(i,j) * g%bathyT(i,j) + g%areaT(i+1,j+1) * g%bathyT(i+1,j+1)) + &
829 (g%areaT(i+1,j) * g%bathyT(i+1,j) + g%areaT(i,j+1) * g%bathyT(i,j+1)) )
836 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
837 if (nonblock_setup)
then
838 call start_group_pass(cs%pass_q_DCor, cs%BT_Domain, clock=id_clock_pass_pre)
840 call do_group_pass(cs%pass_q_DCor, cs%BT_Domain, clock=id_clock_pass_pre)
842 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
847 do j=cs%jsdw,cs%jedw ;
do i=cs%isdw,cs%iedw
848 gtot_e(i,j) = 0.0 ; gtot_w(i,j) = 0.0
849 gtot_n(i,j) = 0.0 ; gtot_s(i,j) = 0.0
852 if (interp_eta_pf)
then
853 eta_pf_1(i,j) = 0.0 ; d_eta_pf(i,j) = 0.0
855 p_surf_dyn(i,j) = 0.0
856 if (cs%dynamic_psurf) dyn_coef_eta(i,j) = 0.0
862 do j=cs%jsdw,cs%jedw ;
do i=cs%isdw-1,cs%iedw
863 cor_ref_u(i,j) = 0.0 ; bt_force_u(i,j) = 0.0 ; ubt(i,j) = 0.0
864 datu(i,j) = 0.0 ; bt_rem_u(i,j) = 0.0 ; uhbt0(i,j) = 0.0
867 do j=cs%jsdw-1,cs%jedw ;
do i=cs%isdw,cs%iedw
868 cor_ref_v(i,j) = 0.0 ; bt_force_v(i,j) = 0.0 ; vbt(i,j) = 0.0
869 datv(i,j) = 0.0 ; bt_rem_v(i,j) = 0.0 ; vhbt0(i,j) = 0.0
873 if (interp_eta_pf)
then
875 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
876 eta(i,j) = eta_in(i,j)
877 eta_pf_1(i,j) = eta_pf_start(i,j)
878 d_eta_pf(i,j) = eta_pf_in(i,j) - eta_pf_start(i,j)
882 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
883 eta(i,j) = eta_in(i,j)
884 eta_pf(i,j) = eta_pf_in(i,j)
889 do k=1,nz ;
do j=js,je ;
do i=is-1,ie
892 if (visc_rem_u(i,j,k) <= 0.0)
then ; visc_rem = 0.0
893 elseif (visc_rem_u(i,j,k) >= 1.0)
then ; visc_rem = 1.0
894 elseif (visc_rem_u(i,j,k)**2 > visc_rem_u(i,j,k) - 0.5*instep)
then
895 visc_rem = visc_rem_u(i,j,k)
896 else ; visc_rem = 1.0 - 0.5*instep/visc_rem_u(i,j,k) ;
endif
897 wt_u(i,j,k) = cs%frhatu(i,j,k) * visc_rem
898 enddo ;
enddo ;
enddo
900 do k=1,nz ;
do j=js-1,je ;
do i=is,ie
902 if (visc_rem_v(i,j,k) <= 0.0)
then ; visc_rem = 0.0
903 elseif (visc_rem_v(i,j,k) >= 1.0)
then ; visc_rem = 1.0
904 elseif (visc_rem_v(i,j,k)**2 > visc_rem_v(i,j,k) - 0.5*instep)
then
905 visc_rem = visc_rem_v(i,j,k)
906 else ; visc_rem = 1.0 - 0.5*instep/visc_rem_v(i,j,k) ;
endif
907 wt_v(i,j,k) = cs%frhatv(i,j,k) * visc_rem
908 enddo ;
enddo ;
enddo
913 do j=js-1,je+1 ;
do i=is-1,ie ; ubt_cor(i,j) = 0.0 ;
enddo ;
enddo
915 do j=js-1,je ;
do i=is-1,ie+1 ; vbt_cor(i,j) = 0.0 ;
enddo ;
enddo
917 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
918 ubt_cor(i,j) = ubt_cor(i,j) + wt_u(i,j,k) * us%m_s_to_L_T*u_cor(i,j,k)
919 enddo ;
enddo ;
enddo
921 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
922 vbt_cor(i,j) = vbt_cor(i,j) + wt_v(i,j,k) * us%m_s_to_L_T*v_cor(i,j,k)
923 enddo ;
enddo ;
enddo
931 do k=1,nz ;
do i=is-1,ie
932 gtot_e(i,j) = gtot_e(i,j) + pbce(i,j,k) * wt_u(i,j,k)
933 gtot_w(i+1,j) = gtot_w(i+1,j) + pbce(i+1,j,k) * wt_u(i,j,k)
938 do k=1,nz ;
do i=is,ie
939 gtot_n(i,j) = gtot_n(i,j) + pbce(i,j,k) * wt_v(i,j,k)
940 gtot_s(i,j+1) = gtot_s(i,j+1) + pbce(i,j+1,k) * wt_v(i,j,k)
945 call tidal_forcing_sensitivity(g, cs%tides_CSp, det_de)
946 dgeo_de = 1.0 + det_de + cs%G_extra
948 dgeo_de = 1.0 + cs%G_extra
951 if (nonblock_setup .and. .not.cs%linearized_BT_PV)
then
952 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
953 call complete_group_pass(cs%pass_q_DCor, cs%BT_Domain, clock=id_clock_pass_pre)
954 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
959 if (use_bt_cont)
then
960 call set_local_bt_cont_types(bt_cont, btcl_u, btcl_v, g, us, ms, cs%BT_Domain, 1+ievf-ie)
962 if (cs%Nonlinear_continuity)
then
963 call find_face_areas(datu, datv, g, gv, us, cs, ms, eta, 1)
965 call find_face_areas(datu, datv, g, gv, us, cs, ms, halo=1)
971 call set_up_bt_obc(obc, eta, cs%BT_OBC, cs%BT_Domain, g, gv, us, ms, ievf-ie, use_bt_cont, &
972 datu, datv, btcl_u, btcl_v)
982 do j=js,je ;
do i=is-1,ie
986 bt_force_u(i,j) = forces%taux(i,j) * mass_to_z * cs%IDatu(i,j)*visc_rem_u(i,j,1)
989 do j=js-1,je ;
do i=is,ie
993 bt_force_v(i,j) = forces%tauy(i,j) * mass_to_z * cs%IDatv(i,j)*visc_rem_v(i,j,1)
995 if (
present(taux_bot) .and.
present(tauy_bot))
then
996 if (
associated(taux_bot) .and.
associated(tauy_bot))
then
998 do j=js,je ;
do i=is-1,ie
999 bt_force_u(i,j) = bt_force_u(i,j) - taux_bot(i,j) * mass_to_z * cs%IDatu(i,j)
1002 do j=js-1,je ;
do i=is,ie
1003 bt_force_v(i,j) = bt_force_v(i,j) - tauy_bot(i,j) * mass_to_z * cs%IDatv(i,j)
1011 do j=js,je ;
do k=1,nz ;
do i=isq,ieq
1012 bt_force_u(i,j) = bt_force_u(i,j) + wt_u(i,j,k) * us%m_to_L*us%T_to_s**2*bc_accel_u(i,j,k)
1013 enddo ;
enddo ;
enddo
1015 do j=jsq,jeq ;
do k=1,nz ;
do i=is,ie
1016 bt_force_v(i,j) = bt_force_v(i,j) + wt_v(i,j,k) * us%m_to_L*us%T_to_s**2*bc_accel_v(i,j,k)
1017 enddo ;
enddo ;
enddo
1023 do j=js,je ;
do i=is-1,ie ; uhbt(i,j) = 0.0 ; ubt(i,j) = 0.0 ;
enddo ;
enddo
1025 do j=js-1,je ;
do i=is,ie ; vhbt(i,j) = 0.0 ; vbt(i,j) = 0.0 ;
enddo ;
enddo
1026 if (cs%visc_rem_u_uh0)
then
1028 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
1029 uhbt(i,j) = uhbt(i,j) + us%T_to_s*us%m_to_L**2*uh0(i,j,k)
1030 ubt(i,j) = ubt(i,j) + wt_u(i,j,k) * us%m_s_to_L_T*u_uh0(i,j,k)
1031 enddo ;
enddo ;
enddo
1033 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
1034 vhbt(i,j) = vhbt(i,j) + us%T_to_s*us%m_to_L**2*vh0(i,j,k)
1035 vbt(i,j) = vbt(i,j) + wt_v(i,j,k) * us%m_s_to_L_T*v_vh0(i,j,k)
1036 enddo ;
enddo ;
enddo
1039 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
1040 uhbt(i,j) = uhbt(i,j) + us%T_to_s*us%m_to_L**2*uh0(i,j,k)
1041 ubt(i,j) = ubt(i,j) + cs%frhatu(i,j,k) * us%m_s_to_L_T*u_uh0(i,j,k)
1042 enddo ;
enddo ;
enddo
1044 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
1045 vhbt(i,j) = vhbt(i,j) + us%T_to_s*us%m_to_L**2*vh0(i,j,k)
1046 vbt(i,j) = vbt(i,j) + cs%frhatv(i,j,k) * us%m_s_to_L_T*v_vh0(i,j,k)
1047 enddo ;
enddo ;
enddo
1049 if (use_bt_cont)
then
1050 if (cs%adjust_BT_cont)
then
1055 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1056 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
1057 call pass_vector(ubt, vbt, cs%BT_Domain, complete=.false., halo=1+ievf-ie)
1058 call pass_vector(uhbt, vhbt, cs%BT_Domain, complete=.true., halo=1+ievf-ie)
1059 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
1060 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1062 call adjust_local_bt_cont_types(ubt, uhbt, vbt, vhbt, btcl_u, btcl_v, &
1063 g, us, ms, 1+ievf-ie)
1066 do j=js,je ;
do i=is-1,ie
1067 uhbt0(i,j) = uhbt(i,j) - find_uhbt(ubt(i,j), btcl_u(i,j), us)
1070 do j=js-1,je ;
do i=is,ie
1071 vhbt0(i,j) = vhbt(i,j) - find_vhbt(vbt(i,j), btcl_v(i,j), us)
1075 do j=js,je ;
do i=is-1,ie
1076 uhbt0(i,j) = uhbt(i,j) - datu(i,j)*ubt(i,j)
1079 do j=js-1,je ;
do i=is,ie
1080 vhbt0(i,j) = vhbt(i,j) - datv(i,j)*vbt(i,j)
1083 if (cs%BT_OBC%apply_u_OBCs)
then
1085 do j=js,je ;
do i=is-1,ie ;
if (obc%segnum_u(i,j) /= obc_none)
then
1087 endif ;
enddo ;
enddo
1089 if (cs%BT_OBC%apply_v_OBCs)
then
1091 do j=js-1,je ;
do i=is,ie ;
if (obc%segnum_v(i,j) /= obc_none)
then
1093 endif ;
enddo ;
enddo
1099 do j=jsvf-1,jevf+1 ;
do i=isvf-2,ievf+1
1100 ubt(i,j) = 0.0 ; uhbt(i,j) = 0.0 ; u_accel_bt(i,j) = 0.0
1103 do j=jsvf-2,jevf+1 ;
do i=isvf-1,ievf+1
1104 vbt(i,j) = 0.0 ; vhbt(i,j) = 0.0 ; v_accel_bt(i,j) = 0.0
1107 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
1108 ubt(i,j) = ubt(i,j) + wt_u(i,j,k) * us%m_s_to_L_T*u_in(i,j,k)
1109 enddo ;
enddo ;
enddo
1111 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
1112 vbt(i,j) = vbt(i,j) + wt_v(i,j,k) * us%m_s_to_L_T*v_in(i,j,k)
1113 enddo ;
enddo ;
enddo
1115 do j=js,je ;
do i=is-1,ie
1116 if (abs(ubt(i,j)) < cs%vel_underflow) ubt(i,j) = 0.0
1119 do j=js-1,je ;
do i=is,ie
1120 if (abs(vbt(i,j)) < cs%vel_underflow) vbt(i,j) = 0.0
1123 if (apply_obcs)
then
1124 ubt_first(:,:) = ubt(:,:) ; vbt_first(:,:) = vbt(:,:)
1127 if (cs%gradual_BT_ICs)
then
1129 do j=js,je ;
do i=is-1,ie
1130 bt_force_u(i,j) = bt_force_u(i,j) + (ubt(i,j) - cs%ubt_IC(i,j)) * idt
1131 ubt(i,j) = cs%ubt_IC(i,j)
1132 if (abs(ubt(i,j)) < cs%vel_underflow) ubt(i,j) = 0.0
1135 do j=js-1,je ;
do i=is,ie
1136 bt_force_v(i,j) = bt_force_v(i,j) + (vbt(i,j) - cs%vbt_IC(i,j)) * idt
1137 vbt(i,j) = cs%vbt_IC(i,j)
1138 if (abs(vbt(i,j)) < cs%vel_underflow) vbt(i,j) = 0.0
1142 if ((isq > is-1) .or. (jsq > js-1))
then
1145 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1146 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
1147 tmp_u(:,:) = 0.0 ; tmp_v(:,:) = 0.0
1148 do j=js,je ;
do i=isq,ieq ; tmp_u(i,j) = bt_force_u(i,j) ;
enddo ;
enddo
1149 do j=jsq,jeq ;
do i=is,ie ; tmp_v(i,j) = bt_force_v(i,j) ;
enddo ;
enddo
1150 if (nonblock_setup)
then
1151 call start_group_pass(cs%pass_tmp_uv, g%Domain)
1153 call do_group_pass(cs%pass_tmp_uv, g%Domain)
1154 do j=jsd,jed ;
do i=isdb,iedb ; bt_force_u(i,j) = tmp_u(i,j) ;
enddo ;
enddo
1155 do j=jsdb,jedb ;
do i=isd,ied ; bt_force_v(i,j) = tmp_v(i,j) ;
enddo ;
enddo
1157 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
1158 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1161 if (nonblock_setup)
then
1162 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1163 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
1164 call start_group_pass(cs%pass_gtot, cs%BT_Domain)
1165 call start_group_pass(cs%pass_ubt_Cor, g%Domain)
1166 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
1167 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1172 do j=jsvf-1,jevf ;
do i=isvf-1,ievf+1
1173 if (cs%Sadourny)
then
1174 amer(i-1,j) = dcor_u(i-1,j) * q(i-1,j)
1175 bmer(i,j) = dcor_u(i,j) * q(i,j)
1176 cmer(i,j+1) = dcor_u(i,j+1) * q(i,j)
1177 dmer(i-1,j+1) = dcor_u(i-1,j+1) * q(i-1,j)
1179 amer(i-1,j) = dcor_u(i-1,j) * &
1180 ((q(i,j) + q(i-1,j-1)) + q(i-1,j)) / 3.0
1181 bmer(i,j) = dcor_u(i,j) * &
1182 (q(i,j) + (q(i-1,j) + q(i,j-1))) / 3.0
1183 cmer(i,j+1) = dcor_u(i,j+1) * &
1184 (q(i,j) + (q(i-1,j) + q(i,j+1))) / 3.0
1185 dmer(i-1,j+1) = dcor_u(i-1,j+1) * &
1186 ((q(i,j) + q(i-1,j+1)) + q(i-1,j)) / 3.0
1191 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf
1192 if (cs%Sadourny)
then
1193 azon(i,j) = dcor_v(i+1,j) * q(i,j)
1194 bzon(i,j) = dcor_v(i,j) * q(i,j)
1195 czon(i,j) = dcor_v(i,j-1) * q(i,j-1)
1196 dzon(i,j) = dcor_v(i+1,j-1) * q(i,j-1)
1198 azon(i,j) = dcor_v(i+1,j) * &
1199 (q(i,j) + (q(i+1,j) + q(i,j-1))) / 3.0
1200 bzon(i,j) = dcor_v(i,j) * &
1201 (q(i,j) + (q(i-1,j) + q(i,j-1))) / 3.0
1202 czon(i,j) = dcor_v(i,j-1) * &
1203 ((q(i,j) + q(i-1,j-1)) + q(i,j-1)) / 3.0
1204 dzon(i,j) = dcor_v(i+1,j-1) * &
1205 ((q(i,j) + q(i+1,j-1)) + q(i,j-1)) / 3.0
1210 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1211 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
1212 if (nonblock_setup)
then
1213 if ((isq > is-1) .or. (jsq > js-1))
then
1214 call complete_group_pass(cs%pass_tmp_uv, g%Domain)
1215 do j=jsd,jed ;
do i=isdb,iedb ; bt_force_u(i,j) = tmp_u(i,j) ;
enddo ;
enddo
1216 do j=jsdb,jedb ;
do i=isd,ied ; bt_force_v(i,j) = tmp_v(i,j) ;
enddo ;
enddo
1218 call complete_group_pass(cs%pass_gtot, cs%BT_Domain)
1219 call complete_group_pass(cs%pass_ubt_Cor, g%Domain)
1221 call do_group_pass(cs%pass_gtot, cs%BT_Domain)
1222 call do_group_pass(cs%pass_ubt_Cor, g%Domain)
1226 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf+1
1227 if (cs%ua_polarity(i,j) < 0.0)
call swap(gtot_e(i,j), gtot_w(i,j))
1228 if (cs%va_polarity(i,j) < 0.0)
call swap(gtot_n(i,j), gtot_s(i,j))
1232 do j=js,je ;
do i=is-1,ie
1234 ((azon(i,j) * vbt_cor(i+1,j) + czon(i,j) * vbt_cor(i ,j-1)) + &
1235 (bzon(i,j) * vbt_cor(i ,j) + dzon(i,j) * vbt_cor(i+1,j-1)))
1238 do j=js-1,je ;
do i=is,ie
1239 cor_ref_v(i,j) = -1.0 * &
1240 ((amer(i-1,j) * ubt_cor(i-1,j) + cmer(i ,j+1) * ubt_cor(i ,j+1)) + &
1241 (bmer(i ,j) * ubt_cor(i ,j) + dmer(i-1,j+1) * ubt_cor(i-1,j+1)))
1245 if (nonblock_setup)
then
1246 if (.not.use_bt_cont) &
1247 call start_group_pass(cs%pass_Dat_uv, cs%BT_Domain)
1250 call start_group_pass(cs%pass_force_hbt0_Cor_ref, cs%BT_Domain)
1252 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
1253 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1264 do j=js-1,je+1 ;
do i=is-1,ie ; av_rem_u(i,j) = 0.0 ;
enddo ;
enddo
1266 do j=js-1,je ;
do i=is-1,ie+1 ; av_rem_v(i,j) = 0.0 ;
enddo ;
enddo
1268 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
1269 av_rem_u(i,j) = av_rem_u(i,j) + cs%frhatu(i,j,k) * visc_rem_u(i,j,k)
1270 enddo ;
enddo ;
enddo
1272 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
1273 av_rem_v(i,j) = av_rem_v(i,j) + cs%frhatv(i,j,k) * visc_rem_v(i,j,k)
1274 enddo ;
enddo ;
enddo
1275 if (cs%strong_drag)
then
1277 do j=js,je ;
do i=is-1,ie
1278 bt_rem_u(i,j) = g%mask2dCu(i,j) * &
1279 ((nstep * av_rem_u(i,j)) / (1.0 + (nstep-1)*av_rem_u(i,j)))
1282 do j=js-1,je ;
do i=is,ie
1283 bt_rem_v(i,j) = g%mask2dCv(i,j) * &
1284 ((nstep * av_rem_v(i,j)) / (1.0 + (nstep-1)*av_rem_v(i,j)))
1288 do j=js,je ;
do i=is-1,ie
1290 if (g%mask2dCu(i,j) * av_rem_u(i,j) > 0.0) &
1291 bt_rem_u(i,j) = g%mask2dCu(i,j) * (av_rem_u(i,j)**instep)
1294 do j=js-1,je ;
do i=is,ie
1296 if (g%mask2dCv(i,j) * av_rem_v(i,j) > 0.0) &
1297 bt_rem_v(i,j) = g%mask2dCv(i,j) * (av_rem_v(i,j)**instep)
1300 if (cs%linear_wave_drag)
then
1302 do j=js,je ;
do i=is-1,ie ;
if (cs%lin_drag_u(i,j) > 0.0)
then
1303 htot = 0.5 * (eta(i,j) + eta(i+1,j))
1304 if (gv%Boussinesq) &
1305 htot = htot + 0.5*gv%Z_to_H * (cs%bathyT(i,j) + cs%bathyT(i+1,j))
1306 bt_rem_u(i,j) = bt_rem_u(i,j) * (htot / (htot + cs%lin_drag_u(i,j) * dtbt))
1308 rayleigh_u(i,j) = cs%lin_drag_u(i,j) / htot
1309 endif ;
enddo ;
enddo
1311 do j=js-1,je ;
do i=is,ie ;
if (cs%lin_drag_v(i,j) > 0.0)
then
1312 htot = 0.5 * (eta(i,j) + eta(i,j+1))
1313 if (gv%Boussinesq) &
1314 htot = htot + 0.5*gv%Z_to_H * (cs%bathyT(i,j) + cs%bathyT(i+1,j+1))
1315 bt_rem_v(i,j) = bt_rem_v(i,j) * (htot / (htot + cs%lin_drag_v(i,j) * dtbt))
1317 rayleigh_v(i,j) = cs%lin_drag_v(i,j) / htot
1318 endif ;
enddo ;
enddo
1322 if (find_etaav)
then
1324 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf+1
1325 eta_sum(i,j) = 0.0 ; eta_wtd(i,j) = 0.0
1329 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf+1
1334 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf
1335 ubt_sum(i,j) = 0.0 ; uhbt_sum(i,j) = 0.0
1336 pfu_bt_sum(i,j) = 0.0 ; coru_bt_sum(i,j) = 0.0
1337 ubt_wtd(i,j) = 0.0 ; ubt_trans(i,j) = 0.0
1340 do j=jsvf-1,jevf ;
do i=isvf-1,ievf+1
1341 vbt_sum(i,j) = 0.0 ; vhbt_sum(i,j) = 0.0
1342 pfv_bt_sum(i,j) = 0.0 ; corv_bt_sum(i,j) = 0.0
1343 vbt_wtd(i,j) = 0.0 ; vbt_trans(i,j) = 0.0
1348 do j=jsvf-1,jevf+1;
do i=isvf-1,ievf+1 ; eta_src(i,j) = 0.0 ;
enddo ;
enddo
1349 if (cs%bound_BT_corr)
then ;
if (use_bt_cont .and. cs%BT_cont_bounds)
then
1350 do j=js,je ;
do i=is,ie ;
if (g%mask2dT(i,j) > 0.0)
then
1351 if (cs%eta_cor(i,j) > 0.0)
then
1355 u_max_cor = us%m_to_L*g%dxT(i,j) * (cs%maxCFL_BT_cont*idt)
1356 v_max_cor = us%m_to_L*g%dyT(i,j) * (cs%maxCFL_BT_cont*idt)
1357 eta_cor_max = dt_in_t * (cs%IareaT(i,j) * &
1358 (((find_uhbt(u_max_cor, btcl_u(i,j), us) + uhbt0(i,j)) - &
1359 (find_uhbt(-u_max_cor, btcl_u(i-1,j), us) + uhbt0(i-1,j))) + &
1360 ((find_vhbt(v_max_cor, btcl_v(i,j), us) + vhbt0(i,j)) - &
1361 (find_vhbt(-v_max_cor, btcl_v(i,j-1), us) + vhbt0(i,j-1))) ))
1362 cs%eta_cor(i,j) = min(cs%eta_cor(i,j), max(0.0, eta_cor_max))
1367 if (gv%Boussinesq) htot = cs%bathyT(i,j)*gv%Z_to_H + eta(i,j)
1369 cs%eta_cor(i,j) = max(cs%eta_cor(i,j), -max(0.0,htot))
1371 endif ;
enddo ;
enddo
1372 else ;
do j=js,je ;
do i=is,ie
1373 if (abs(cs%eta_cor(i,j)) > dt_in_t*cs%eta_cor_bound(i,j)) &
1374 cs%eta_cor(i,j) = sign(dt_in_t*cs%eta_cor_bound(i,j), cs%eta_cor(i,j))
1375 enddo ;
enddo ;
endif ;
endif
1377 do j=js,je ;
do i=is,ie
1378 eta_src(i,j) = g%mask2dT(i,j) * (instep * cs%eta_cor(i,j))
1382 if (cs%dynamic_psurf)
then
1383 ice_is_rigid = (
associated(forces%rigidity_ice_u) .and. &
1384 associated(forces%rigidity_ice_v))
1385 h_min_dyn = gv%Z_to_H * cs%Dmin_dyn_psurf
1386 if (ice_is_rigid .and. use_bt_cont) &
1387 call bt_cont_to_face_areas(bt_cont, datu, datv, g, us, ms, 0, .true.)
1388 if (ice_is_rigid)
then
1390 do j=js,je ;
do i=is,ie
1396 idt_max2 = 0.5 * (dgeo_de * (1.0 + 2.0*bebt)) * (us%L_to_m**2*g%IareaT(i,j) * &
1397 ((gtot_e(i,j) * (datu(i,j)*us%L_to_m*g%IdxCu(i,j)) + &
1398 gtot_w(i,j) * (datu(i-1,j)*us%L_to_m*g%IdxCu(i-1,j))) + &
1399 (gtot_n(i,j) * (datv(i,j)*us%L_to_m*g%IdyCv(i,j)) + &
1400 gtot_s(i,j) * (datv(i,j-1)*us%L_to_m*g%IdyCv(i,j-1)))) + &
1401 ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
1402 (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2)))
1403 h_eff_dx2 = max(h_min_dyn * ((us%L_to_m*g%IdxT(i,j))**2 + (us%L_to_m*g%IdyT(i,j))**2), &
1404 us%L_to_m**2*g%IareaT(i,j) * &
1405 ((datu(i,j)*us%L_to_m*g%IdxCu(i,j) + datu(i-1,j)*us%L_to_m*g%IdxCu(i-1,j)) + &
1406 (datv(i,j)*us%L_to_m*g%IdyCv(i,j) + datv(i,j-1)*us%L_to_m*g%IdyCv(i,j-1)) ) )
1407 dyn_coef_max = cs%const_dyn_psurf * max(0.0, 1.0 - dtbt**2 * idt_max2) / &
1408 (dtbt**2 * h_eff_dx2)
1411 ice_strength = us%m_to_L**4*us%Z_to_m*us%T_to_s* &
1412 ((forces%rigidity_ice_u(i,j) + forces%rigidity_ice_u(i-1,j)) + &
1413 (forces%rigidity_ice_v(i,j) + forces%rigidity_ice_v(i,j-1))) / &
1414 (cs%ice_strength_length**2 * dtbt)
1417 dyn_coef_eta(i,j) = min(dyn_coef_max, ice_strength * gv%H_to_Z)
1418 enddo ;
enddo ;
endif
1421 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1422 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
1423 if (nonblock_setup)
then
1424 call start_group_pass(cs%pass_eta_bt_rem, cs%BT_Domain)
1427 call do_group_pass(cs%pass_eta_bt_rem, cs%BT_Domain)
1428 if (.not.use_bt_cont) &
1429 call do_group_pass(cs%pass_Dat_uv, cs%BT_Domain)
1430 call do_group_pass(cs%pass_force_hbt0_Cor_ref, cs%BT_Domain)
1432 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
1433 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1436 if (nonblock_setup)
then
1437 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1438 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
1440 if (.not.use_bt_cont) &
1441 call complete_group_pass(cs%pass_Dat_uv, cs%BT_Domain)
1442 call complete_group_pass(cs%pass_force_hbt0_Cor_ref, cs%BT_Domain)
1443 call complete_group_pass(cs%pass_eta_bt_rem, cs%BT_Domain)
1445 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
1446 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1450 call uvchksum(
"BT [uv]hbt", uhbt, vhbt, cs%debug_BT_HI, haloshift=0, &
1451 scale=us%s_to_T*us%L_to_m**2*gv%H_to_m)
1452 call uvchksum(
"BT Initial [uv]bt", ubt, vbt, cs%debug_BT_HI, haloshift=0, scale=us%L_T_to_m_s)
1453 call hchksum(eta,
"BT Initial eta", cs%debug_BT_HI, haloshift=0, scale=gv%H_to_m)
1454 call uvchksum(
"BT BT_force_[uv]", bt_force_u, bt_force_v, &
1455 cs%debug_BT_HI, haloshift=0, scale=us%L_T2_to_m_s2)
1456 if (interp_eta_pf)
then
1457 call hchksum(eta_pf_1,
"BT eta_PF_1",cs%debug_BT_HI,haloshift=0, scale=gv%H_to_m)
1458 call hchksum(d_eta_pf,
"BT d_eta_PF",cs%debug_BT_HI,haloshift=0, scale=gv%H_to_m)
1460 call hchksum(eta_pf,
"BT eta_PF",cs%debug_BT_HI,haloshift=0, scale=gv%H_to_m)
1461 call hchksum(eta_pf_in,
"BT eta_PF_in",g%HI,haloshift=0, scale=gv%H_to_m)
1463 call uvchksum(
"BT Cor_ref_[uv]", cor_ref_u, cor_ref_v, cs%debug_BT_HI, haloshift=0, scale=us%L_T2_to_m_s2)
1464 call uvchksum(
"BT [uv]hbt0", uhbt0, vhbt0, cs%debug_BT_HI, haloshift=0, &
1465 scale=us%L_to_m**2*us%s_to_T*gv%H_to_m)
1466 if (.not. use_bt_cont)
then
1467 call uvchksum(
"BT Dat[uv]", datu, datv, cs%debug_BT_HI, haloshift=1, scale=us%L_to_m*gv%H_to_m)
1469 call uvchksum(
"BT wt_[uv]", wt_u, wt_v, g%HI, 0, .true., .true.)
1470 call uvchksum(
"BT frhat[uv]", cs%frhatu, cs%frhatv, g%HI, 0, .true., .true.)
1471 call uvchksum(
"BT bc_accel_[uv]", bc_accel_u, bc_accel_v, g%HI, haloshift=0)
1472 call uvchksum(
"BT IDat[uv]", cs%IDatu, cs%IDatv, g%HI, haloshift=0, scale=us%m_to_Z)
1473 call uvchksum(
"BT visc_rem_[uv]", visc_rem_u, visc_rem_v, g%HI, haloshift=1)
1476 if (query_averaging_enabled(cs%diag))
then
1477 if (cs%id_eta_st > 0)
call post_data(cs%id_eta_st, eta(isd:ied,jsd:jed), cs%diag)
1478 if (cs%id_ubt_st > 0)
call post_data(cs%id_ubt_st, ubt(isdb:iedb,jsd:jed), cs%diag)
1479 if (cs%id_vbt_st > 0)
call post_data(cs%id_vbt_st, vbt(isd:ied,jsdb:jedb), cs%diag)
1482 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1483 if (id_clock_calc > 0)
call cpu_clock_begin(id_clock_calc)
1485 if (project_velocity)
then ; eta_pf_bt => eta ;
else ; eta_pf_bt => eta_pred ;
endif
1487 if (cs%dt_bt_filter >= 0.0)
then
1488 dt_filt = 0.5 * max(0.0, min(cs%dt_bt_filter, 2.0*dt_in_t))
1490 dt_filt = 0.5 * max(0.0, dt_in_t * min(-cs%dt_bt_filter, 2.0))
1492 nfilter = ceiling(dt_filt / dtbt)
1494 if (nstep+nfilter==0 )
call mom_error(fatal, &
1495 "btstep: number of barotropic step (nstep+nfilter) is 0")
1498 sum_wt_vel = 0.0 ; sum_wt_eta = 0.0 ; sum_wt_accel = 0.0 ; sum_wt_trans = 0.0
1499 allocate(wt_vel(nstep+nfilter)) ;
allocate(wt_eta(nstep+nfilter))
1500 allocate(wt_trans(nstep+nfilter+1)) ;
allocate(wt_accel(nstep+nfilter+1))
1501 allocate(wt_accel2(nstep+nfilter+1))
1502 do n=1,nstep+nfilter
1504 if ( (n==nstep) .or. (dt_filt - abs(n-nstep)*dtbt >= 0.0))
then
1505 wt_vel(n) = 1.0 ; wt_eta(n) = 1.0
1506 elseif (dtbt + dt_filt - abs(n-nstep)*dtbt > 0.0)
then
1507 wt_vel(n) = 1.0 + (dt_filt / dtbt) - abs(n-nstep) ; wt_eta(n) = wt_vel(n)
1509 wt_vel(n) = 0.0 ; wt_eta(n) = 0.0
1515 sum_wt_vel = sum_wt_vel + wt_vel(n) ; sum_wt_eta = sum_wt_eta + wt_eta(n)
1517 wt_trans(nstep+nfilter+1) = 0.0 ; wt_accel(nstep+nfilter+1) = 0.0
1518 do n=nstep+nfilter,1,-1
1519 wt_trans(n) = wt_trans(n+1) + wt_eta(n)
1520 wt_accel(n) = wt_accel(n+1) + wt_vel(n)
1521 sum_wt_accel = sum_wt_accel + wt_accel(n) ; sum_wt_trans = sum_wt_trans + wt_trans(n)
1524 i_sum_wt_vel = 1.0 / sum_wt_vel ; i_sum_wt_accel = 1.0 / sum_wt_accel
1525 i_sum_wt_eta = 1.0 / sum_wt_eta ; i_sum_wt_trans = 1.0 / sum_wt_trans
1526 do n=1,nstep+nfilter
1527 wt_vel(n) = wt_vel(n) * i_sum_wt_vel
1528 wt_accel2(n) = wt_accel(n)
1529 wt_accel(n) = wt_accel(n) * i_sum_wt_accel
1530 wt_eta(n) = wt_eta(n) * i_sum_wt_eta
1535 sum_wt_vel = 0.0 ; sum_wt_eta = 0.0 ; sum_wt_accel = 0.0 ; sum_wt_trans = 0.0
1538 isv=is ; iev=ie ; jsv=js ; jev=je
1539 do n=1,nstep+nfilter
1541 sum_wt_vel = sum_wt_vel + wt_vel(n)
1542 sum_wt_eta = sum_wt_eta + wt_eta(n)
1543 sum_wt_accel = sum_wt_accel + wt_accel2(n)
1544 sum_wt_trans = sum_wt_trans + wt_trans(n)
1546 if (cs%clip_velocity)
then
1547 do j=jsv,jev ;
do i=isv-1,iev
1548 if ((ubt(i,j) * (dt_in_t * us%m_to_L*g%dy_Cu(i,j))) * us%L_to_m**2*g%IareaT(i+1,j) < -cs%CFL_trunc)
then
1550 ubt(i,j) = (-0.95*cs%CFL_trunc) * (us%m_to_L**2*g%areaT(i+1,j) / (dt_in_t * us%m_to_L*g%dy_Cu(i,j)))
1551 elseif ((ubt(i,j) * (dt_in_t * us%m_to_L*g%dy_Cu(i,j))) * us%L_to_m**2*g%IareaT(i,j) > cs%CFL_trunc)
then
1553 ubt(i,j) = (0.95*cs%CFL_trunc) * (us%m_to_L**2*g%areaT(i,j) / (dt_in_t * us%m_to_L*g%dy_Cu(i,j)))
1556 do j=jsv-1,jev ;
do i=isv,iev
1557 if ((vbt(i,j) * (dt_in_t * us%m_to_L*g%dx_Cv(i,j))) * us%L_to_m**2*g%IareaT(i,j+1) < -cs%CFL_trunc)
then
1559 vbt(i,j) = (-0.9*cs%CFL_trunc) * (us%m_to_L**2*g%areaT(i,j+1) / (dt_in_t * us%m_to_L*g%dx_Cv(i,j)))
1560 elseif ((vbt(i,j) * (dt_in_t * us%m_to_L*g%dx_Cv(i,j))) * us%L_to_m**2*g%IareaT(i,j) > cs%CFL_trunc)
then
1562 vbt(i,j) = (0.9*cs%CFL_trunc) * (us%m_to_L**2*g%areaT(i,j) / (dt_in_t * us%m_to_L*g%dx_Cv(i,j)))
1567 if ((iev - stencil < ie) .or. (jev - stencil < je))
then
1568 if (id_clock_calc > 0)
call cpu_clock_end(id_clock_calc)
1569 call do_group_pass(cs%pass_eta_ubt, cs%BT_Domain, clock=id_clock_pass_step)
1570 isv = isvf ; iev = ievf ; jsv = jsvf ; jev = jevf
1571 if (id_clock_calc > 0)
call cpu_clock_begin(id_clock_calc)
1573 isv = isv+stencil ; iev = iev-stencil
1574 jsv = jsv+stencil ; jev = jev-stencil
1577 if ((.not.use_bt_cont) .and. cs%Nonlinear_continuity .and. &
1578 (cs%Nonlin_cont_update_period > 0))
then
1579 if ((n>1) .and. (mod(n-1,cs%Nonlin_cont_update_period) == 0)) &
1580 call find_face_areas(datu, datv, g, gv, us, cs, ms, eta, 1+iev-ie)
1584 if (cs%dynamic_psurf .or. .not.project_velocity)
then
1585 if (use_bt_cont)
then
1587 do j=jsv-1,jev+1 ;
do i=isv-2,iev+1
1588 uhbt(i,j) = find_uhbt(ubt(i,j), btcl_u(i,j), us) + uhbt0(i,j)
1591 do j=jsv-2,jev+1 ;
do i=isv-1,iev+1
1592 vhbt(i,j) = find_vhbt(vbt(i,j), btcl_v(i,j), us) + vhbt0(i,j)
1595 do j=jsv-1,jev+1 ;
do i=isv-1,iev+1
1596 eta_pred(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * cs%IareaT(i,j)) * &
1597 ((uhbt(i-1,j) - uhbt(i,j)) + (vhbt(i,j-1) - vhbt(i,j)))
1601 do j=jsv-1,jev+1 ;
do i=isv-1,iev+1
1602 eta_pred(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * cs%IareaT(i,j)) * &
1603 (((datu(i-1,j)*ubt(i-1,j) + uhbt0(i-1,j)) - &
1604 (datu(i,j)*ubt(i,j) + uhbt0(i,j))) + &
1605 ((datv(i,j-1)*vbt(i,j-1) + vhbt0(i,j-1)) - &
1606 (datv(i,j)*vbt(i,j) + vhbt0(i,j))))
1610 if (cs%dynamic_psurf)
then
1612 do j=jsv-1,jev+1 ;
do i=isv-1,iev+1
1613 p_surf_dyn(i,j) = dyn_coef_eta(i,j) * (eta_pred(i,j) - eta(i,j))
1621 if (find_etaav)
then
1623 do j=js,je ;
do i=is,ie
1624 eta_sum(i,j) = eta_sum(i,j) + wt_accel2(n) * eta_pf_bt(i,j)
1628 if (interp_eta_pf)
then
1631 do j=jsv-1,jev+1 ;
do i=isv-1,iev+1
1632 eta_pf(i,j) = eta_pf_1(i,j) + wt_end*d_eta_pf(i,j)
1636 if (apply_obc_flather .or. apply_obc_open)
then
1638 do j=jsv,jev ;
do i=isv-2,iev+1
1639 ubt_old(i,j) = ubt(i,j)
1642 do j=jsv-2,jev+1 ;
do i=isv,iev
1643 vbt_old(i,j) = vbt(i,j)
1648 if (apply_obcs)
then
1649 if (mod(n+g%first_direction,2)==1)
then
1655 if (cs%BT_OBC%apply_u_OBCs)
then
1657 do j=jsv-joff,jev+joff ;
do i=isv-1,iev
1658 ubt_prev(i,j) = ubt(i,j) ; uhbt_prev(i,j) = uhbt(i,j)
1659 ubt_sum_prev(i,j) = ubt_sum(i,j) ; uhbt_sum_prev(i,j) = uhbt_sum(i,j) ; ubt_wtd_prev(i,j) = ubt_wtd(i,j)
1663 if (cs%BT_OBC%apply_v_OBCs)
then
1665 do j=jsv-1,jev ;
do i=isv-ioff,iev+ioff
1666 vbt_prev(i,j) = vbt(i,j) ; vhbt_prev(i,j) = vhbt(i,j)
1667 vbt_sum_prev(i,j) = vbt_sum(i,j) ; vhbt_sum_prev(i,j) = vhbt_sum(i,j) ; vbt_wtd_prev(i,j) = vbt_wtd(i,j)
1673 if (mod(n+g%first_direction,2)==1)
then
1676 do j=jsv-1,jev ;
do i=isv-1,iev+1
1677 cor_v(i,j) = -1.0*((amer(i-1,j) * ubt(i-1,j) + cmer(i,j+1) * ubt(i,j+1)) + &
1678 (bmer(i,j) * ubt(i,j) + dmer(i-1,j+1) * ubt(i-1,j+1))) - cor_ref_v(i,j)
1679 pfv(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_n(i,j) - &
1680 (eta_pf_bt(i,j+1)-eta_pf(i,j+1))*gtot_s(i,j+1)) * &
1681 dgeo_de * cs%IdyCv(i,j)
1683 if (cs%dynamic_psurf)
then
1685 do j=jsv-1,jev ;
do i=isv-1,iev+1
1686 pfv(i,j) = pfv(i,j) + (p_surf_dyn(i,j) - p_surf_dyn(i,j+1)) * cs%IdyCv(i,j)
1690 if (cs%BT_OBC%apply_v_OBCs)
then
1692 do j=jsv-1,jev ;
do i=isv-1,iev+1 ;
if (obc%segnum_v(i,j) /= obc_none)
then
1694 endif ;
enddo ;
enddo
1697 do j=jsv-1,jev ;
do i=isv-1,iev+1
1699 vbt(i,j) = bt_rem_v(i,j) * (vbt(i,j) + &
1700 dtbt * ((bt_force_v(i,j) + cor_v(i,j)) + pfv(i,j)))
1701 vbt_trans(i,j) = trans_wt1*vbt(i,j) + trans_wt2*vel_prev
1703 if (cs%linear_wave_drag)
then
1704 v_accel_bt(i,j) = v_accel_bt(i,j) + wt_accel(n) * &
1705 ((cor_v(i,j) + pfv(i,j)) - vbt(i,j)*rayleigh_v(i,j))
1707 v_accel_bt(i,j) = v_accel_bt(i,j) + wt_accel(n) * (cor_v(i,j) + pfv(i,j))
1711 if (use_bt_cont)
then
1713 do j=jsv-1,jev ;
do i=isv-1,iev+1
1714 vhbt(i,j) = find_vhbt(vbt_trans(i,j), btcl_v(i,j), us) + vhbt0(i,j)
1718 do j=jsv-1,jev ;
do i=isv-1,iev+1
1719 vhbt(i,j) = datv(i,j)*vbt_trans(i,j) + vhbt0(i,j)
1722 if (cs%BT_OBC%apply_v_OBCs)
then
1724 do j=jsv-1,jev ;
do i=isv-1,iev+1 ;
if (obc%segnum_v(i,j) /= obc_none)
then
1725 vbt(i,j) = vbt_prev(i,j) ; vhbt(i,j) = vhbt_prev(i,j)
1726 endif ;
enddo ;
enddo
1730 do j=jsv,jev ;
do i=isv-1,iev
1731 cor_u(i,j) = ((azon(i,j) * vbt(i+1,j) + czon(i,j) * vbt(i,j-1)) + &
1732 (bzon(i,j) * vbt(i,j) + dzon(i,j) * vbt(i+1,j-1))) - &
1734 pfu(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_e(i,j) - &
1735 (eta_pf_bt(i+1,j)-eta_pf(i+1,j))*gtot_w(i+1,j)) * &
1736 dgeo_de * cs%IdxCu(i,j)
1739 if (cs%dynamic_psurf)
then
1741 do j=jsv,jev ;
do i=isv-1,iev
1742 pfu(i,j) = pfu(i,j) + (p_surf_dyn(i,j) - p_surf_dyn(i+1,j)) * cs%IdxCu(i,j)
1746 if (cs%BT_OBC%apply_u_OBCs)
then
1748 do j=jsv,jev ;
do i=isv-1,iev ;
if (obc%segnum_u(i,j) /= obc_none)
then
1750 endif ;
enddo ;
enddo
1753 do j=jsv,jev ;
do i=isv-1,iev
1755 ubt(i,j) = bt_rem_u(i,j) * (ubt(i,j) + &
1756 dtbt * ((bt_force_u(i,j) + cor_u(i,j)) + pfu(i,j)))
1757 if (abs(ubt(i,j)) < cs%vel_underflow) ubt(i,j) = 0.0
1758 ubt_trans(i,j) = trans_wt1*ubt(i,j) + trans_wt2*vel_prev
1760 if (cs%linear_wave_drag)
then
1761 u_accel_bt(i,j) = u_accel_bt(i,j) + wt_accel(n) * &
1762 ((cor_u(i,j) + pfu(i,j)) - ubt(i,j)*rayleigh_u(i,j))
1764 u_accel_bt(i,j) = u_accel_bt(i,j) + wt_accel(n) * (cor_u(i,j) + pfu(i,j))
1768 if (use_bt_cont)
then
1770 do j=jsv,jev ;
do i=isv-1,iev
1771 uhbt(i,j) = find_uhbt(ubt_trans(i,j), btcl_u(i,j), us) + uhbt0(i,j)
1775 do j=jsv,jev ;
do i=isv-1,iev
1776 uhbt(i,j) = datu(i,j)*ubt_trans(i,j) + uhbt0(i,j)
1779 if (cs%BT_OBC%apply_u_OBCs)
then
1781 do j=jsv,jev ;
do i=isv-1,iev ;
if (obc%segnum_u(i,j) /= obc_none)
then
1782 ubt(i,j) = ubt_prev(i,j); uhbt(i,j) = uhbt_prev(i,j)
1783 endif ;
enddo ;
enddo
1788 do j=jsv-1,jev+1 ;
do i=isv-1,iev
1789 cor_u(i,j) = ((azon(i,j) * vbt(i+1,j) + czon(i,j) * vbt(i,j-1)) + &
1790 (bzon(i,j) * vbt(i,j) + dzon(i,j) * vbt(i+1,j-1))) - &
1792 pfu(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_e(i,j) - &
1793 (eta_pf_bt(i+1,j)-eta_pf(i+1,j))*gtot_w(i+1,j)) * &
1794 dgeo_de * cs%IdxCu(i,j)
1797 if (cs%dynamic_psurf)
then
1799 do j=jsv-1,jev+1 ;
do i=isv-1,iev
1800 pfu(i,j) = pfu(i,j) + (p_surf_dyn(i,j) - p_surf_dyn(i+1,j)) * cs%IdxCu(i,j)
1804 if (cs%BT_OBC%apply_u_OBCs)
then
1806 do j=jsv,jev ;
do i=isv-1,iev ;
if (obc%segnum_u(i,j) /= obc_none)
then
1808 endif ;
enddo ;
enddo
1812 do j=jsv-1,jev+1 ;
do i=isv-1,iev
1814 ubt(i,j) = bt_rem_u(i,j) * (ubt(i,j) + &
1815 dtbt * ((bt_force_u(i,j) + cor_u(i,j)) + pfu(i,j)))
1816 if (abs(ubt(i,j)) < cs%vel_underflow) ubt(i,j) = 0.0
1817 ubt_trans(i,j) = trans_wt1*ubt(i,j) + trans_wt2*vel_prev
1819 if (cs%linear_wave_drag)
then
1820 u_accel_bt(i,j) = u_accel_bt(i,j) + wt_accel(n) * &
1821 ((cor_u(i,j) + pfu(i,j)) - ubt(i,j)*rayleigh_u(i,j))
1823 u_accel_bt(i,j) = u_accel_bt(i,j) + wt_accel(n) * (cor_u(i,j) + pfu(i,j))
1827 if (use_bt_cont)
then
1829 do j=jsv-1,jev+1 ;
do i=isv-1,iev
1830 uhbt(i,j) = find_uhbt(ubt_trans(i,j), btcl_u(i,j), us) + uhbt0(i,j)
1834 do j=jsv-1,jev+1 ;
do i=isv-1,iev
1835 uhbt(i,j) = datu(i,j)*ubt_trans(i,j) + uhbt0(i,j)
1838 if (cs%BT_OBC%apply_u_OBCs)
then
1840 do j=jsv-1,jev+1 ;
do i=isv-1,iev ;
if (obc%segnum_u(i,j) /= obc_none)
then
1841 ubt(i,j) = ubt_prev(i,j); uhbt(i,j) = uhbt_prev(i,j)
1842 endif ;
enddo ;
enddo
1846 if (cs%use_old_coriolis_bracket_bug)
then
1848 do j=jsv-1,jev ;
do i=isv,iev
1849 cor_v(i,j) = -1.0*((amer(i-1,j) * ubt(i-1,j) + bmer(i,j) * ubt(i,j)) + &
1850 (cmer(i,j+1) * ubt(i,j+1) + dmer(i-1,j+1) * ubt(i-1,j+1))) - cor_ref_v(i,j)
1851 pfv(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_n(i,j) - &
1852 (eta_pf_bt(i,j+1)-eta_pf(i,j+1))*gtot_s(i,j+1)) * &
1853 dgeo_de * cs%IdyCv(i,j)
1857 do j=jsv-1,jev ;
do i=isv,iev
1858 cor_v(i,j) = -1.0*((amer(i-1,j) * ubt(i-1,j) + cmer(i,j+1) * ubt(i,j+1)) + &
1859 (bmer(i,j) * ubt(i,j) + dmer(i-1,j+1) * ubt(i-1,j+1))) - cor_ref_v(i,j)
1860 pfv(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_n(i,j) - &
1861 (eta_pf_bt(i,j+1)-eta_pf(i,j+1))*gtot_s(i,j+1)) * &
1862 dgeo_de * cs%IdyCv(i,j)
1866 if (cs%dynamic_psurf)
then
1868 do j=jsv-1,jev ;
do i=isv,iev
1869 pfv(i,j) = pfv(i,j) + (p_surf_dyn(i,j) - p_surf_dyn(i,j+1)) * cs%IdyCv(i,j)
1873 if (cs%BT_OBC%apply_v_OBCs)
then
1875 do j=jsv-1,jev ;
do i=isv-1,iev+1 ;
if (obc%segnum_v(i,j) /= obc_none)
then
1877 endif ;
enddo ;
enddo
1881 do j=jsv-1,jev ;
do i=isv,iev
1883 vbt(i,j) = bt_rem_v(i,j) * (vbt(i,j) + &
1884 dtbt * ((bt_force_v(i,j) + cor_v(i,j)) + pfv(i,j)))
1885 if (abs(vbt(i,j)) < cs%vel_underflow) vbt(i,j) = 0.0
1886 vbt_trans(i,j) = trans_wt1*vbt(i,j) + trans_wt2*vel_prev
1888 if (cs%linear_wave_drag)
then
1889 v_accel_bt(i,j) = v_accel_bt(i,j) + wt_accel(n) * &
1890 ((cor_v(i,j) + pfv(i,j)) - vbt(i,j)*rayleigh_v(i,j))
1892 v_accel_bt(i,j) = v_accel_bt(i,j) + wt_accel(n) * (cor_v(i,j) + pfv(i,j))
1895 if (use_bt_cont)
then
1897 do j=jsv-1,jev ;
do i=isv,iev
1898 vhbt(i,j) = find_vhbt(vbt_trans(i,j), btcl_v(i,j), us) + vhbt0(i,j)
1902 do j=jsv-1,jev ;
do i=isv,iev
1903 vhbt(i,j) = datv(i,j)*vbt_trans(i,j) + vhbt0(i,j)
1906 if (cs%BT_OBC%apply_v_OBCs)
then
1908 do j=jsv-1,jev ;
do i=isv,iev ;
if (obc%segnum_v(i,j) /= obc_none)
then
1909 vbt(i,j) = vbt_prev(i,j); vhbt(i,j) = vhbt_prev(i,j)
1910 endif ;
enddo ;
enddo
1918 do j=js,je ;
do i=is-1,ie
1919 pfu_bt_sum(i,j) = pfu_bt_sum(i,j) + wt_accel2(n) * pfu(i,j)
1922 do j=js-1,je ;
do i=is,ie
1923 pfv_bt_sum(i,j) = pfv_bt_sum(i,j) + wt_accel2(n) * pfv(i,j)
1928 do j=js,je ;
do i=is-1,ie
1929 coru_bt_sum(i,j) = coru_bt_sum(i,j) + wt_accel2(n) * cor_u(i,j)
1932 do j=js-1,je ;
do i=is,ie
1933 corv_bt_sum(i,j) = corv_bt_sum(i,j) + wt_accel2(n) * cor_v(i,j)
1938 do j=js,je ;
do i=is-1,ie
1939 ubt_sum(i,j) = ubt_sum(i,j) + wt_trans(n) * ubt_trans(i,j)
1940 uhbt_sum(i,j) = uhbt_sum(i,j) + wt_trans(n) * uhbt(i,j)
1941 ubt_wtd(i,j) = ubt_wtd(i,j) + wt_vel(n) * ubt(i,j)
1944 do j=js-1,je ;
do i=is,ie
1945 vbt_sum(i,j) = vbt_sum(i,j) + wt_trans(n) * vbt_trans(i,j)
1946 vhbt_sum(i,j) = vhbt_sum(i,j) + wt_trans(n) * vhbt(i,j)
1947 vbt_wtd(i,j) = vbt_wtd(i,j) + wt_vel(n) * vbt(i,j)
1951 if (apply_obcs)
then
1952 if (cs%BT_OBC%apply_u_OBCs)
then
1954 do j=js,je ;
do i=is-1,ie
1955 if (obc%segnum_u(i,j) /= obc_none)
then
1956 ubt_sum(i,j) = ubt_sum_prev(i,j) ; uhbt_sum(i,j) = uhbt_sum_prev(i,j)
1957 ubt_wtd(i,j) = ubt_wtd_prev(i,j)
1962 if (cs%BT_OBC%apply_v_OBCs)
then
1964 do j=js-1,je ;
do i=is,ie
1965 if (obc%segnum_v(i,j) /= obc_none)
then
1966 vbt_sum(i,j) = vbt_sum_prev(i,j) ; vhbt_sum(i,j) = vhbt_sum_prev(i,j)
1967 vbt_wtd(i,j) = vbt_wtd_prev(i,j)
1972 call apply_velocity_obcs(obc, ubt, vbt, uhbt, vhbt, &
1973 ubt_trans, vbt_trans, eta, ubt_old, vbt_old, cs%BT_OBC, &
1974 g, ms, us, iev-ie, dtbt, bebt, use_bt_cont, datu, datv, btcl_u, btcl_v, &
1976 if (cs%BT_OBC%apply_u_OBCs)
then ;
do j=js,je ;
do i=is-1,ie
1977 if (obc%segnum_u(i,j) /= obc_none)
then
1978 ubt_sum(i,j) = ubt_sum(i,j) + wt_trans(n) * ubt_trans(i,j)
1979 uhbt_sum(i,j) = uhbt_sum(i,j) + wt_trans(n) * uhbt(i,j)
1980 ubt_wtd(i,j) = ubt_wtd(i,j) + wt_vel(n) * ubt(i,j)
1982 enddo ;
enddo ;
endif
1983 if (cs%BT_OBC%apply_v_OBCs)
then ;
do j=js-1,je ;
do i=is,ie
1984 if (obc%segnum_v(i,j) /= obc_none)
then
1985 vbt_sum(i,j) = vbt_sum(i,j) + wt_trans(n) * vbt_trans(i,j)
1986 vhbt_sum(i,j) = vhbt_sum(i,j) + wt_trans(n) * vhbt(i,j)
1987 vbt_wtd(i,j) = vbt_wtd(i,j) + wt_vel(n) * vbt(i,j)
1989 enddo ;
enddo ;
endif
1992 if (cs%debug_bt)
then
1993 call uvchksum(
"BT [uv]hbt just after OBC", uhbt, vhbt, cs%debug_BT_HI, haloshift=iev-ie, &
1994 scale=us%s_to_T*us%L_to_m**2*gv%H_to_m)
1998 do j=jsv,jev ;
do i=isv,iev
1999 eta(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * cs%IareaT(i,j)) * &
2000 ((uhbt(i-1,j) - uhbt(i,j)) + (vhbt(i,j-1) - vhbt(i,j)))
2001 eta_wtd(i,j) = eta_wtd(i,j) + eta(i,j) * wt_eta(n)
2005 if (do_hifreq_output)
then
2006 time_step_end = time_bt_start + real_to_time(n*us%T_to_s*dtbt)
2007 call enable_averaging(us%T_to_s*dtbt, time_step_end, cs%diag)
2008 if (cs%id_ubt_hifreq > 0)
call post_data(cs%id_ubt_hifreq, ubt(isdb:iedb,jsd:jed), cs%diag)
2009 if (cs%id_vbt_hifreq > 0)
call post_data(cs%id_vbt_hifreq, vbt(isd:ied,jsdb:jedb), cs%diag)
2010 if (cs%id_eta_hifreq > 0)
call post_data(cs%id_eta_hifreq, eta(isd:ied,jsd:jed), cs%diag)
2011 if (cs%id_uhbt_hifreq > 0)
call post_data(cs%id_uhbt_hifreq, uhbt(isdb:iedb,jsd:jed), cs%diag)
2012 if (cs%id_vhbt_hifreq > 0)
call post_data(cs%id_vhbt_hifreq, vhbt(isd:ied,jsdb:jedb), cs%diag)
2013 if (cs%id_eta_pred_hifreq > 0)
call post_data(cs%id_eta_pred_hifreq, eta_pf_bt(isd:ied,jsd:jed), cs%diag)
2016 if (cs%debug_bt)
then
2017 write(mesg,
'("BT step ",I4)') n
2018 call uvchksum(trim(mesg)//
" [uv]bt", ubt, vbt, cs%debug_BT_HI, haloshift=iev-ie, &
2019 scale=us%L_T_to_m_s)
2020 call hchksum(eta, trim(mesg)//
" eta", cs%debug_BT_HI, haloshift=iev-ie, scale=gv%H_to_m)
2024 if (id_clock_calc > 0)
call cpu_clock_end(id_clock_calc)
2025 if (id_clock_calc_post > 0)
call cpu_clock_begin(id_clock_calc_post)
2028 if (do_hifreq_output)
call enable_averaging(time_int_in, time_end_in, cs%diag)
2030 i_sum_wt_vel = 1.0 / sum_wt_vel ; i_sum_wt_eta = 1.0 / sum_wt_eta
2031 i_sum_wt_accel = 1.0 / sum_wt_accel ; i_sum_wt_trans = 1.0 / sum_wt_trans
2033 if (find_etaav)
then ;
do j=js,je ;
do i=is,ie
2034 etaav(i,j) = eta_sum(i,j) * i_sum_wt_accel
2035 enddo ;
enddo ;
endif
2036 do j=js-1,je+1 ;
do i=is-1,ie+1 ; e_anom(i,j) = 0.0 ;
enddo ;
enddo
2037 if (interp_eta_pf)
then
2038 do j=js,je ;
do i=is,ie
2039 e_anom(i,j) = dgeo_de * (0.5 * (eta(i,j) + eta_in(i,j)) - &
2040 (eta_pf_1(i,j) + 0.5*d_eta_pf(i,j)))
2043 do j=js,je ;
do i=is,ie
2044 e_anom(i,j) = dgeo_de * (0.5 * (eta(i,j) + eta_in(i,j)) - eta_pf(i,j))
2047 if (apply_obcs)
then
2049 if (cs%BT_OBC%apply_u_OBCs)
then
2051 do j=js,je ;
do i=is-1,ie
2052 if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e)
then
2053 e_anom(i+1,j) = e_anom(i,j)
2054 elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w)
then
2055 e_anom(i,j) = e_anom(i+1,j)
2060 if (cs%BT_OBC%apply_v_OBCs)
then
2062 do j=js-1,je ;
do i=is,ie
2063 if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n)
then
2064 e_anom(i,j+1) = e_anom(i,j)
2065 elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s)
then
2066 e_anom(i,j) = e_anom(i,j+1)
2073 do j=js,je ;
do i=is,ie
2074 eta_out(i,j) = eta_wtd(i,j) * i_sum_wt_eta
2077 if (id_clock_calc_post > 0)
call cpu_clock_end(id_clock_calc_post)
2078 if (id_clock_pass_post > 0)
call cpu_clock_begin(id_clock_pass_post)
2079 if (g%nonblocking_updates)
then
2080 call start_group_pass(cs%pass_e_anom, g%Domain)
2082 if (find_etaav)
call do_group_pass(cs%pass_etaav, g%Domain)
2083 call do_group_pass(cs%pass_e_anom, g%Domain)
2085 if (id_clock_pass_post > 0)
call cpu_clock_end(id_clock_pass_post)
2086 if (id_clock_calc_post > 0)
call cpu_clock_begin(id_clock_calc_post)
2088 do j=js,je ;
do i=is-1,ie
2089 cs%ubtav(i,j) = ubt_sum(i,j) * i_sum_wt_trans
2090 uhbtav(i,j) = us%s_to_T*us%L_to_m**2*uhbt_sum(i,j) * i_sum_wt_trans
2093 ubt_wtd(i,j) = ubt_wtd(i,j) * i_sum_wt_vel
2096 do j=js-1,je ;
do i=is,ie
2097 cs%vbtav(i,j) = vbt_sum(i,j) * i_sum_wt_trans
2098 vhbtav(i,j) = us%s_to_T*us%L_to_m**2*vhbt_sum(i,j) * i_sum_wt_trans
2101 vbt_wtd(i,j) = vbt_wtd(i,j) * i_sum_wt_vel
2104 if (id_clock_calc_post > 0)
call cpu_clock_end(id_clock_calc_post)
2105 if (id_clock_pass_post > 0)
call cpu_clock_begin(id_clock_pass_post)
2106 if (g%nonblocking_updates)
then
2107 call complete_group_pass(cs%pass_e_anom, g%Domain)
2108 if (find_etaav)
call start_group_pass(cs%pass_etaav, g%Domain)
2109 call start_group_pass(cs%pass_ubta_uhbta, g%DoMain)
2111 call do_group_pass(cs%pass_ubta_uhbta, g%Domain)
2113 if (id_clock_pass_post > 0)
call cpu_clock_end(id_clock_pass_post)
2114 if (id_clock_calc_post > 0)
call cpu_clock_begin(id_clock_calc_post)
2119 do j=js,je ;
do i=is-1,ie
2120 accel_layer_u(i,j,k) = us%L_T2_to_m_s2 * (u_accel_bt(i,j) - &
2121 ((pbce(i+1,j,k) - gtot_w(i+1,j)) * e_anom(i+1,j) - &
2122 (pbce(i,j,k) - gtot_e(i,j)) * e_anom(i,j)) * cs%IdxCu(i,j) )
2123 if (abs(accel_layer_u(i,j,k)) < accel_underflow) accel_layer_u(i,j,k) = 0.0
2125 do j=js-1,je ;
do i=is,ie
2126 accel_layer_v(i,j,k) = us%L_T2_to_m_s2 * (v_accel_bt(i,j) - &
2127 ((pbce(i,j+1,k) - gtot_s(i,j+1)) * e_anom(i,j+1) - &
2128 (pbce(i,j,k) - gtot_n(i,j)) * e_anom(i,j)) * cs%IdyCv(i,j) )
2129 if (abs(accel_layer_v(i,j,k)) < accel_underflow) accel_layer_v(i,j,k) = 0.0
2133 if (apply_obcs)
then
2136 if (cs%BT_OBC%apply_u_OBCs)
then ;
do j=js,je ;
do i=is-1,ie
2137 if (obc%segnum_u(i,j) /= obc_none)
then
2138 u_accel_bt(i,j) = (ubt_wtd(i,j) - ubt_first(i,j)) / dt_in_t
2139 do k=1,nz ; accel_layer_u(i,j,k) = us%L_T2_to_m_s2*u_accel_bt(i,j) ;
enddo
2141 enddo ;
enddo ;
endif
2142 if (cs%BT_OBC%apply_v_OBCs)
then ;
do j=js-1,je ;
do i=is,ie
2143 if (obc%segnum_v(i,j) /= obc_none)
then
2144 v_accel_bt(i,j) = (vbt_wtd(i,j) - vbt_first(i,j)) / dt_in_t
2145 do k=1,nz ; accel_layer_v(i,j,k) = us%L_T2_to_m_s2*v_accel_bt(i,j) ;
enddo
2147 enddo ;
enddo ;
endif
2150 if (id_clock_calc_post > 0)
call cpu_clock_end(id_clock_calc_post)
2153 if (query_averaging_enabled(cs%diag))
then
2155 do j=js,je ;
do i=is-1,ie ; cs%ubt_IC(i,j) = ubt_wtd(i,j) ;
enddo ;
enddo
2156 do j=js-1,je ;
do i=is,ie ; cs%vbt_IC(i,j) = vbt_wtd(i,j) ;
enddo ;
enddo
2157 if (use_bt_cont)
then
2158 do j=js,je ;
do i=is-1,ie
2159 cs%uhbt_IC(i,j) = find_uhbt(ubt_wtd(i,j), btcl_u(i,j), us) + uhbt0(i,j)
2161 do j=js-1,je ;
do i=is,ie
2162 cs%vhbt_IC(i,j) = find_vhbt(vbt_wtd(i,j), btcl_v(i,j), us) + vhbt0(i,j)
2165 do j=js,je ;
do i=is-1,ie
2166 cs%uhbt_IC(i,j) = ubt_wtd(i,j) * datu(i,j) + uhbt0(i,j)
2168 do j=js-1,je ;
do i=is,ie
2169 cs%vhbt_IC(i,j) = vbt_wtd(i,j) * datv(i,j) + vhbt0(i,j)
2174 if (cs%id_PFu_bt > 0)
then
2175 do j=js,je ;
do i=is-1,ie
2176 pfu_bt_sum(i,j) = pfu_bt_sum(i,j) * i_sum_wt_accel
2178 call post_data(cs%id_PFu_bt, pfu_bt_sum(isdb:iedb,jsd:jed), cs%diag)
2180 if (cs%id_PFv_bt > 0)
then
2181 do j=js-1,je ;
do i=is,ie
2182 pfv_bt_sum(i,j) = pfv_bt_sum(i,j) * i_sum_wt_accel
2184 call post_data(cs%id_PFv_bt, pfv_bt_sum(isd:ied,jsdb:jedb), cs%diag)
2186 if (cs%id_Coru_bt > 0)
then
2187 do j=js,je ;
do i=is-1,ie
2188 coru_bt_sum(i,j) = coru_bt_sum(i,j) * i_sum_wt_accel
2190 call post_data(cs%id_Coru_bt, coru_bt_sum(isdb:iedb,jsd:jed), cs%diag)
2192 if (cs%id_Corv_bt > 0)
then
2193 do j=js-1,je ;
do i=is,ie
2194 corv_bt_sum(i,j) = corv_bt_sum(i,j) * i_sum_wt_accel
2196 call post_data(cs%id_Corv_bt, corv_bt_sum(isd:ied,jsdb:jedb), cs%diag)
2198 if (cs%id_ubtforce > 0)
call post_data(cs%id_ubtforce, bt_force_u(isdb:iedb,jsd:jed), cs%diag)
2199 if (cs%id_vbtforce > 0)
call post_data(cs%id_vbtforce, bt_force_v(isd:ied,jsdb:jedb), cs%diag)
2200 if (cs%id_uaccel > 0)
call post_data(cs%id_uaccel, u_accel_bt(isdb:iedb,jsd:jed), cs%diag)
2201 if (cs%id_vaccel > 0)
call post_data(cs%id_vaccel, v_accel_bt(isd:ied,jsdb:jedb), cs%diag)
2203 if (cs%id_eta_cor > 0)
call post_data(cs%id_eta_cor, cs%eta_cor, cs%diag)
2204 if (cs%id_eta_bt > 0)
call post_data(cs%id_eta_bt, eta_out, cs%diag)
2205 if (cs%id_gtotn > 0)
call post_data(cs%id_gtotn, gtot_n(isd:ied,jsd:jed), cs%diag)
2206 if (cs%id_gtots > 0)
call post_data(cs%id_gtots, gtot_s(isd:ied,jsd:jed), cs%diag)
2207 if (cs%id_gtote > 0)
call post_data(cs%id_gtote, gtot_e(isd:ied,jsd:jed), cs%diag)
2208 if (cs%id_gtotw > 0)
call post_data(cs%id_gtotw, gtot_w(isd:ied,jsd:jed), cs%diag)
2209 if (cs%id_ubt > 0)
call post_data(cs%id_ubt, ubt_wtd(isdb:iedb,jsd:jed), cs%diag)
2210 if (cs%id_vbt > 0)
call post_data(cs%id_vbt, vbt_wtd(isd:ied,jsdb:jedb), cs%diag)
2211 if (cs%id_ubtav > 0)
call post_data(cs%id_ubtav, cs%ubtav, cs%diag)
2212 if (cs%id_vbtav > 0)
call post_data(cs%id_vbtav, cs%vbtav, cs%diag)
2213 if (cs%id_visc_rem_u > 0)
call post_data(cs%id_visc_rem_u, visc_rem_u, cs%diag)
2214 if (cs%id_visc_rem_v > 0)
call post_data(cs%id_visc_rem_v, visc_rem_v, cs%diag)
2216 if (cs%id_frhatu > 0)
call post_data(cs%id_frhatu, cs%frhatu, cs%diag)
2217 if (cs%id_uhbt > 0)
call post_data(cs%id_uhbt, uhbtav, cs%diag)
2218 if (cs%id_frhatv > 0)
call post_data(cs%id_frhatv, cs%frhatv, cs%diag)
2219 if (cs%id_vhbt > 0)
call post_data(cs%id_vhbt, vhbtav, cs%diag)
2220 if (cs%id_uhbt0 > 0)
call post_data(cs%id_uhbt0, uhbt0(isdb:iedb,jsd:jed), cs%diag)
2221 if (cs%id_vhbt0 > 0)
call post_data(cs%id_vhbt0, vhbt0(isd:ied,jsdb:jedb), cs%diag)
2223 if (cs%id_frhatu1 > 0)
call post_data(cs%id_frhatu1, cs%frhatu1, cs%diag)
2224 if (cs%id_frhatv1 > 0)
call post_data(cs%id_frhatv1, cs%frhatv1, cs%diag)
2226 if (use_bt_cont)
then
2227 if (cs%id_BTC_FA_u_EE > 0)
call post_data(cs%id_BTC_FA_u_EE, bt_cont%FA_u_EE, cs%diag)
2228 if (cs%id_BTC_FA_u_E0 > 0)
call post_data(cs%id_BTC_FA_u_E0, bt_cont%FA_u_E0, cs%diag)
2229 if (cs%id_BTC_FA_u_W0 > 0)
call post_data(cs%id_BTC_FA_u_W0, bt_cont%FA_u_W0, cs%diag)
2230 if (cs%id_BTC_FA_u_WW > 0)
call post_data(cs%id_BTC_FA_u_WW, bt_cont%FA_u_WW, cs%diag)
2231 if (cs%id_BTC_uBT_EE > 0)
call post_data(cs%id_BTC_uBT_EE, bt_cont%uBT_EE, cs%diag)
2232 if (cs%id_BTC_uBT_WW > 0)
call post_data(cs%id_BTC_uBT_WW, bt_cont%uBT_WW, cs%diag)
2233 if (cs%id_BTC_FA_v_NN > 0)
call post_data(cs%id_BTC_FA_v_NN, bt_cont%FA_v_NN, cs%diag)
2234 if (cs%id_BTC_FA_v_N0 > 0)
call post_data(cs%id_BTC_FA_v_N0, bt_cont%FA_v_N0, cs%diag)
2235 if (cs%id_BTC_FA_v_S0 > 0)
call post_data(cs%id_BTC_FA_v_S0, bt_cont%FA_v_S0, cs%diag)
2236 if (cs%id_BTC_FA_v_SS > 0)
call post_data(cs%id_BTC_FA_v_SS, bt_cont%FA_v_SS, cs%diag)
2237 if (cs%id_BTC_vBT_NN > 0)
call post_data(cs%id_BTC_vBT_NN, bt_cont%vBT_NN, cs%diag)
2238 if (cs%id_BTC_vBT_SS > 0)
call post_data(cs%id_BTC_vBT_SS, bt_cont%vBT_SS, cs%diag)
2241 if (cs%id_frhatu1 > 0) cs%frhatu1(:,:,:) = cs%frhatu(:,:,:)
2242 if (cs%id_frhatv1 > 0) cs%frhatv1(:,:,:) = cs%frhatv(:,:,:)
2245 if (g%nonblocking_updates)
then
2246 if (find_etaav)
call complete_group_pass(cs%pass_etaav, g%Domain)
2247 call complete_group_pass(cs%pass_ubta_uhbta, g%Domain)