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.
411 type(ocean_grid_type),
intent(inout) :: G
412 type(verticalGrid_type),
intent(in) :: GV
413 type(unit_scale_type),
intent(in) :: US
414 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: U_in
416 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: V_in
418 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: eta_in
420 real,
intent(in) :: dt
421 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: bc_accel_u
423 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: bc_accel_v
425 type(mech_forcing),
intent(in) :: forces
426 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: pbce
429 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: eta_PF_in
435 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: U_Cor
437 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: V_Cor
439 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(out) :: accel_layer_u
441 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(out) :: accel_layer_v
443 real,
dimension(SZI_(G),SZJ_(G)),
intent(out) :: eta_out
445 real,
dimension(SZIB_(G),SZJ_(G)),
intent(out) :: uhbtav
448 real,
dimension(SZI_(G),SZJB_(G)),
intent(out) :: vhbtav
451 type(barotropic_CS),
pointer :: CS
453 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: visc_rem_u
459 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: visc_rem_v
460 real,
dimension(SZI_(G),SZJ_(G)),
optional,
intent(out) :: etaav
462 type(accel_diag_ptrs),
optional,
pointer :: ADp
463 type(ocean_OBC_type),
optional,
pointer :: OBC
464 type(BT_cont_type),
optional,
pointer :: BT_cont
467 real,
dimension(:,:),
optional,
pointer :: eta_PF_start
470 real,
dimension(:,:),
optional,
pointer :: taux_bot
472 real,
dimension(:,:),
optional,
pointer :: tauy_bot
474 real,
dimension(:,:,:),
optional,
pointer :: uh0
476 real,
dimension(:,:,:),
optional,
pointer :: u_uh0
478 real,
dimension(:,:,:),
optional,
pointer :: vh0
480 real,
dimension(:,:,:),
optional,
pointer :: v_vh0
484 real :: ubt_Cor(SZIB_(G),SZJ_(G))
485 real :: vbt_Cor(SZI_(G),SZJB_(G))
487 real :: wt_u(SZIB_(G),SZJ_(G),SZK_(G))
488 real :: wt_v(SZI_(G),SZJB_(G),SZK_(G))
491 real,
dimension(SZIB_(G),SZJ_(G)) :: &
496 real,
dimension(SZI_(G),SZJB_(G)) :: &
501 real,
dimension(SZI_(G),SZJ_(G)) :: &
508 real :: q(SZIBW_(CS),SZJBW_(CS))
510 real,
dimension(SZIBW_(CS),SZJW_(CS)) :: &
545 real,
dimension(SZIW_(CS),SZJBW_(CS)) :: &
578 real,
target,
dimension(SZIW_(CS),SZJW_(CS)) :: &
582 real,
dimension(:,:),
pointer :: &
585 real,
dimension(SZIW_(CS),SZJW_(CS)) :: &
607 type(local_BT_cont_u_type),
dimension(SZIBW_(CS),SZJW_(CS)) :: &
609 type(local_BT_cont_v_type),
dimension(SZIW_(CS),SZJBW_(CS)) :: &
613 real,
dimension(SZIBW_(CS),SZJW_(CS)) :: &
614 ubt_prev, ubt_sum_prev, ubt_wtd_prev, &
615 uhbt_prev, uhbt_sum_prev, &
618 real,
dimension(SZIW_(CS),SZJBW_(CS)) :: &
619 vbt_prev, vbt_sum_prev, vbt_wtd_prev, &
620 vhbt_prev, vhbt_sum_prev, &
624 real :: mass_accel_to_Z
650 logical :: do_hifreq_output
651 logical :: use_BT_cont, do_ave, find_etaav, find_PF, find_Cor
652 logical :: integral_BT_cont
654 logical :: ice_is_rigid, nonblock_setup, interp_eta_PF
655 logical :: project_velocity, add_uh0
659 real :: ice_strength = 0.0
666 real :: u_max_cor, v_max_cor
667 real :: uint_cor, vint_cor
670 real :: accel_underflow
675 real,
allocatable,
dimension(:) :: wt_vel, wt_eta, wt_accel, wt_trans, wt_accel2
676 real :: sum_wt_vel, sum_wt_eta, sum_wt_accel, sum_wt_trans
677 real :: I_sum_wt_vel, I_sum_wt_eta, I_sum_wt_accel, I_sum_wt_trans
679 real :: trans_wt1, trans_wt2
682 logical :: apply_OBCs, apply_OBC_flather, apply_OBC_open
683 type(memory_size_type) :: MS
684 character(len=200) :: mesg
685 integer :: isv, iev, jsv, jev
687 integer :: isvf, ievf, jsvf, jevf, num_cycles
688 integer :: i, j, k, n
689 integer :: is, ie, js, je, nz, Isq, Ieq, Jsq, Jeq
690 integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
691 integer :: ioff, joff
694 if (.not.
associated(cs))
call mom_error(fatal, &
695 "btstep: Module MOM_barotropic must be initialized before it is used.")
696 if (.not.cs%split)
return
697 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
698 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
699 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
700 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
701 ms%isdw = cs%isdw ; ms%iedw = cs%iedw ; ms%jsdw = cs%jsdw ; ms%jedw = cs%jedw
702 h_neglect = gv%H_subroundoff
705 accel_underflow = cs%vel_underflow * idt
707 use_bt_cont = .false.
708 if (
present(bt_cont)) use_bt_cont = (
associated(bt_cont))
709 integral_bt_cont = use_bt_cont .and. cs%integral_BT_cont
711 interp_eta_pf = .false.
712 if (
present(eta_pf_start)) interp_eta_pf = (
associated(eta_pf_start))
714 project_velocity = cs%BT_project_velocity
718 if ((.not.use_bt_cont) .and. cs%Nonlinear_continuity .and. &
719 (cs%Nonlin_cont_update_period > 0)) stencil = 2
721 do_ave = query_averaging_enabled(cs%diag)
722 find_etaav =
present(etaav)
723 find_pf = (do_ave .and. ((cs%id_PFu_bt > 0) .or. (cs%id_PFv_bt > 0)))
724 find_cor = (do_ave .and. ((cs%id_Coru_bt > 0) .or. (cs%id_Corv_bt > 0)))
727 if (
present(uh0)) add_uh0 =
associated(uh0)
728 if (add_uh0 .and. .not.(
present(vh0) .and.
present(u_uh0) .and. &
729 present(v_vh0)))
call mom_error(fatal, &
730 "btstep: vh0, u_uh0, and v_vh0 must be present if uh0 is used.")
731 if (add_uh0 .and. .not.(
associated(vh0) .and.
associated(u_uh0) .and. &
732 associated(v_vh0)))
call mom_error(fatal, &
733 "btstep: vh0, u_uh0, and v_vh0 must be associated if uh0 is used.")
736 nonblock_setup = g%nonblocking_updates
738 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
740 apply_obcs = .false. ; cs%BT_OBC%apply_u_OBCs = .false. ; cs%BT_OBC%apply_v_OBCs = .false.
741 apply_obc_open = .false.
742 apply_obc_flather = .false.
743 if (
present(obc))
then ;
if (
associated(obc))
then
744 cs%BT_OBC%apply_u_OBCs = obc%open_u_BCs_exist_globally .or. obc%specified_u_BCs_exist_globally
745 cs%BT_OBC%apply_v_OBCs = obc%open_v_BCs_exist_globally .or. obc%specified_v_BCs_exist_globally
746 apply_obc_flather = open_boundary_query(obc, apply_flather_obc=.true.)
747 apply_obc_open = open_boundary_query(obc, apply_open_obc=.true.)
748 apply_obcs = open_boundary_query(obc, apply_specified_obc=.true.) .or. &
749 apply_obc_flather .or. apply_obc_open
751 if (apply_obc_flather .and. .not.gv%Boussinesq)
call mom_error(fatal, &
752 "btstep: Flather open boundary conditions have not yet been "// &
753 "implemented for a non-Boussinesq model.")
757 if (cs%use_wide_halos) &
758 num_cycles = min((is-cs%isdw) / stencil, (js-cs%jsdw) / stencil)
759 isvf = is - (num_cycles-1)*stencil ; ievf = ie + (num_cycles-1)*stencil
760 jsvf = js - (num_cycles-1)*stencil ; jevf = je + (num_cycles-1)*stencil
762 nstep = ceiling(dt/cs%dtbt - 0.0001)
763 if (is_root_pe() .and. (nstep /= cs%nstep_last))
then
764 write(mesg,
'("btstep is using a dynamic barotropic timestep of ", ES12.6, &
765 & " seconds, max ", ES12.6, ".")') (us%T_to_s*dt/nstep), us%T_to_s*cs%dtbt_max
766 call mom_mesg(mesg, 3)
768 cs%nstep_last = nstep
771 instep = 1.0 / real(nstep)
776 mass_accel_to_z = 1.0 / gv%Rho0
777 mass_to_z = us%m_to_Z / gv%Rho0
780 if (project_velocity)
then
781 trans_wt1 = (1.0 + be_proj); trans_wt2 = -be_proj
783 trans_wt1 = bebt ; trans_wt2 = (1.0-bebt)
786 do_hifreq_output = .false.
787 if ((cs%id_ubt_hifreq > 0) .or. (cs%id_vbt_hifreq > 0) .or. &
788 (cs%id_eta_hifreq > 0) .or. (cs%id_eta_pred_hifreq > 0) .or. &
789 (cs%id_uhbt_hifreq > 0) .or. (cs%id_vhbt_hifreq > 0))
then
790 do_hifreq_output = query_averaging_enabled(cs%diag, time_int_in, time_end_in)
791 if (do_hifreq_output) &
792 time_bt_start = time_end_in - real_to_time(us%T_to_s*dt)
796 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
797 if (.not. cs%linearized_BT_PV)
then
798 call create_group_pass(cs%pass_q_DCor, q, cs%BT_Domain, to_all, position=corner)
799 call create_group_pass(cs%pass_q_DCor, dcor_u, dcor_v, cs%BT_Domain, &
802 if ((isq > is-1) .or. (jsq > js-1)) &
803 call create_group_pass(cs%pass_tmp_uv, tmp_u, tmp_v, g%Domain)
804 call create_group_pass(cs%pass_gtot, gtot_e, gtot_n, cs%BT_Domain, &
805 to_all+scalar_pair, agrid)
806 call create_group_pass(cs%pass_gtot, gtot_w, gtot_s, cs%BT_Domain, &
807 to_all+scalar_pair, agrid)
809 if (cs%dynamic_psurf) &
810 call create_group_pass(cs%pass_eta_bt_rem, dyn_coef_eta, cs%BT_Domain)
811 if (interp_eta_pf)
then
812 call create_group_pass(cs%pass_eta_bt_rem, eta_pf_1, cs%BT_Domain)
813 call create_group_pass(cs%pass_eta_bt_rem, d_eta_pf, cs%BT_Domain)
815 call create_group_pass(cs%pass_eta_bt_rem, eta_pf, cs%BT_Domain)
817 if (integral_bt_cont) &
818 call create_group_pass(cs%pass_eta_bt_rem, eta_ic, cs%BT_Domain)
819 call create_group_pass(cs%pass_eta_bt_rem, eta_src, cs%BT_Domain)
823 call create_group_pass(cs%pass_eta_bt_rem, bt_rem_u, bt_rem_v, &
824 cs%BT_Domain, to_all+scalar_pair)
825 if (cs%linear_wave_drag) &
826 call create_group_pass(cs%pass_eta_bt_rem, rayleigh_u, rayleigh_v, &
827 cs%BT_Domain, to_all+scalar_pair)
830 if (((g%isd > cs%isdw) .or. (g%jsd > cs%jsdw)) .or. (isq <= is-1) .or. (jsq <= js-1)) &
831 call create_group_pass(cs%pass_force_hbt0_Cor_ref, bt_force_u, bt_force_v, cs%BT_Domain)
832 if (add_uh0)
call create_group_pass(cs%pass_force_hbt0_Cor_ref, uhbt0, vhbt0, cs%BT_Domain)
833 call create_group_pass(cs%pass_force_hbt0_Cor_ref, cor_ref_u, cor_ref_v, cs%BT_Domain)
834 if (.not. use_bt_cont)
then
835 call create_group_pass(cs%pass_Dat_uv, datu, datv, cs%BT_Domain, to_all+scalar_pair)
837 call create_group_pass(cs%pass_eta_ubt, eta, cs%BT_Domain)
838 call create_group_pass(cs%pass_eta_ubt, ubt, vbt, cs%BT_Domain)
839 if (integral_bt_cont)
then
840 call create_group_pass(cs%pass_eta_ubt, ubt_int, vbt_int, cs%BT_Domain)
842 if (apply_obc_open) &
843 call create_group_pass(cs%pass_eta_ubt, uhbt_int, vhbt_int, cs%BT_Domain)
846 call create_group_pass(cs%pass_ubt_Cor, ubt_cor, vbt_cor, g%Domain)
850 call create_group_pass(cs%pass_etaav, etaav, g%Domain)
852 call create_group_pass(cs%pass_e_anom, e_anom, g%Domain)
853 call create_group_pass(cs%pass_ubta_uhbta, cs%ubtav, cs%vbtav, g%Domain)
854 call create_group_pass(cs%pass_ubta_uhbta, uhbtav, vhbtav, g%Domain)
856 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
862 if (cs%linearized_BT_PV)
then
864 do j=jsvf-2,jevf+1 ;
do i=isvf-2,ievf+1
868 do j=jsvf-1,jevf+1 ;
do i=isvf-2,ievf+1
869 dcor_u(i,j) = cs%D_u_Cor(i,j)
872 do j=jsvf-2,jevf+1 ;
do i=isvf-1,ievf+1
873 dcor_v(i,j) = cs%D_v_Cor(i,j)
876 q(:,:) = 0.0 ; dcor_u(:,:) = 0.0 ; dcor_v(:,:) = 0.0
877 if (gv%Boussinesq)
then
879 do j=js,je ;
do i=is-1,ie
880 dcor_u(i,j) = 0.5 * (max(gv%Z_to_H*g%bathyT(i+1,j) + eta_in(i+1,j), 0.0) + &
881 max(gv%Z_to_H*g%bathyT(i,j) + eta_in(i,j), 0.0) )
884 do j=js-1,je ;
do i=is,ie
885 dcor_v(i,j) = 0.5 * (max(gv%Z_to_H*g%bathyT(i,j+1) + eta_in(i+1,j), 0.0) + &
886 max(gv%Z_to_H*g%bathyT(i,j) + eta_in(i,j), 0.0) )
889 do j=js-1,je ;
do i=is-1,ie
890 q(i,j) = 0.25 * (cs%BT_Coriolis_scale * g%CoriolisBu(i,j)) * &
891 ((g%areaT(i,j) + g%areaT(i+1,j+1)) + (g%areaT(i+1,j) + g%areaT(i,j+1))) / &
892 (max((g%areaT(i,j) * max(gv%Z_to_H*g%bathyT(i,j) + eta_in(i,j), 0.0) + &
893 g%areaT(i+1,j+1) * max(gv%Z_to_H*g%bathyT(i+1,j+1) + eta_in(i+1,j+1), 0.0)) + &
894 (g%areaT(i+1,j) * max(gv%Z_to_H*g%bathyT(i+1,j) + eta_in(i+1,j), 0.0) + &
895 g%areaT(i,j+1) * max(gv%Z_to_H*g%bathyT(i,j+1) + eta_in(i,j+1), 0.0)), h_neglect) )
899 do j=js,je ;
do i=is-1,ie
900 dcor_u(i,j) = 0.5 * (eta_in(i+1,j) + eta_in(i,j))
903 do j=js-1,je ;
do i=is,ie
904 dcor_v(i,j) = 0.5 * (eta_in(i,j+1) + eta_in(i,j))
907 do j=js-1,je ;
do i=is-1,ie
908 q(i,j) = 0.25 * (cs%BT_Coriolis_scale * g%CoriolisBu(i,j)) * &
909 ((g%areaT(i,j) + g%areaT(i+1,j+1)) + (g%areaT(i+1,j) + g%areaT(i,j+1))) / &
910 (max((g%areaT(i,j) * eta_in(i,j) + g%areaT(i+1,j+1) * eta_in(i+1,j+1)) + &
911 (g%areaT(i+1,j) * eta_in(i+1,j) + g%areaT(i,j+1) * eta_in(i,j+1)), h_neglect) )
919 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
920 if (nonblock_setup)
then
921 call start_group_pass(cs%pass_q_DCor, cs%BT_Domain, clock=id_clock_pass_pre)
923 call do_group_pass(cs%pass_q_DCor, cs%BT_Domain, clock=id_clock_pass_pre)
925 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
930 do j=cs%jsdw,cs%jedw ;
do i=cs%isdw,cs%iedw
931 gtot_e(i,j) = 0.0 ; gtot_w(i,j) = 0.0
932 gtot_n(i,j) = 0.0 ; gtot_s(i,j) = 0.0
935 if (interp_eta_pf)
then
936 eta_pf_1(i,j) = 0.0 ; d_eta_pf(i,j) = 0.0
938 if (integral_bt_cont)
then
941 p_surf_dyn(i,j) = 0.0
942 if (cs%dynamic_psurf) dyn_coef_eta(i,j) = 0.0
948 do j=cs%jsdw,cs%jedw ;
do i=cs%isdw-1,cs%iedw
949 cor_ref_u(i,j) = 0.0 ; bt_force_u(i,j) = 0.0 ; ubt(i,j) = 0.0
950 datu(i,j) = 0.0 ; bt_rem_u(i,j) = 0.0 ; uhbt0(i,j) = 0.0
953 do j=cs%jsdw-1,cs%jedw ;
do i=cs%isdw,cs%iedw
954 cor_ref_v(i,j) = 0.0 ; bt_force_v(i,j) = 0.0 ; vbt(i,j) = 0.0
955 datv(i,j) = 0.0 ; bt_rem_v(i,j) = 0.0 ; vhbt0(i,j) = 0.0
959 if (interp_eta_pf)
then
961 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
962 eta(i,j) = eta_in(i,j)
963 eta_pf_1(i,j) = eta_pf_start(i,j)
964 d_eta_pf(i,j) = eta_pf_in(i,j) - eta_pf_start(i,j)
968 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
969 eta(i,j) = eta_in(i,j)
970 eta_pf(i,j) = eta_pf_in(i,j)
973 if (integral_bt_cont)
then
975 do j=g%jsd,g%jed ;
do i=g%isd,g%ied
976 eta_ic(i,j) = eta_in(i,j)
981 do k=1,nz ;
do j=js,je ;
do i=is-1,ie
984 if (visc_rem_u(i,j,k) <= 0.0)
then ; visc_rem = 0.0
985 elseif (visc_rem_u(i,j,k) >= 1.0)
then ; visc_rem = 1.0
986 elseif (visc_rem_u(i,j,k)**2 > visc_rem_u(i,j,k) - 0.5*instep)
then
987 visc_rem = visc_rem_u(i,j,k)
988 else ; visc_rem = 1.0 - 0.5*instep/visc_rem_u(i,j,k) ;
endif
989 wt_u(i,j,k) = cs%frhatu(i,j,k) * visc_rem
990 enddo ;
enddo ;
enddo
992 do k=1,nz ;
do j=js-1,je ;
do i=is,ie
994 if (visc_rem_v(i,j,k) <= 0.0)
then ; visc_rem = 0.0
995 elseif (visc_rem_v(i,j,k) >= 1.0)
then ; visc_rem = 1.0
996 elseif (visc_rem_v(i,j,k)**2 > visc_rem_v(i,j,k) - 0.5*instep)
then
997 visc_rem = visc_rem_v(i,j,k)
998 else ; visc_rem = 1.0 - 0.5*instep/visc_rem_v(i,j,k) ;
endif
999 wt_v(i,j,k) = cs%frhatv(i,j,k) * visc_rem
1000 enddo ;
enddo ;
enddo
1005 do j=js-1,je+1 ;
do i=is-1,ie ; ubt_cor(i,j) = 0.0 ;
enddo ;
enddo
1007 do j=js-1,je ;
do i=is-1,ie+1 ; vbt_cor(i,j) = 0.0 ;
enddo ;
enddo
1009 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
1010 ubt_cor(i,j) = ubt_cor(i,j) + wt_u(i,j,k) * u_cor(i,j,k)
1011 enddo ;
enddo ;
enddo
1013 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
1014 vbt_cor(i,j) = vbt_cor(i,j) + wt_v(i,j,k) * v_cor(i,j,k)
1015 enddo ;
enddo ;
enddo
1023 do k=1,nz ;
do i=is-1,ie
1024 gtot_e(i,j) = gtot_e(i,j) + pbce(i,j,k) * wt_u(i,j,k)
1025 gtot_w(i+1,j) = gtot_w(i+1,j) + pbce(i+1,j,k) * wt_u(i,j,k)
1030 do k=1,nz ;
do i=is,ie
1031 gtot_n(i,j) = gtot_n(i,j) + pbce(i,j,k) * wt_v(i,j,k)
1032 gtot_s(i,j+1) = gtot_s(i,j+1) + pbce(i,j+1,k) * wt_v(i,j,k)
1037 call tidal_forcing_sensitivity(g, cs%tides_CSp, det_de)
1038 dgeo_de = 1.0 + det_de + cs%G_extra
1040 dgeo_de = 1.0 + cs%G_extra
1043 if (nonblock_setup .and. .not.cs%linearized_BT_PV)
then
1044 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1045 call complete_group_pass(cs%pass_q_DCor, cs%BT_Domain, clock=id_clock_pass_pre)
1046 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1051 if (integral_bt_cont)
then
1052 call set_local_bt_cont_types(bt_cont, btcl_u, btcl_v, g, us, ms, cs%BT_Domain, 1+ievf-ie, dt_baroclinic=dt)
1053 elseif (use_bt_cont)
then
1054 call set_local_bt_cont_types(bt_cont, btcl_u, btcl_v, g, us, ms, cs%BT_Domain, 1+ievf-ie)
1056 if (cs%Nonlinear_continuity)
then
1057 call find_face_areas(datu, datv, g, gv, us, cs, ms, eta, 1)
1059 call find_face_areas(datu, datv, g, gv, us, cs, ms, halo=1)
1064 if (apply_obcs)
then
1065 call set_up_bt_obc(obc, eta, cs%BT_OBC, cs%BT_Domain, g, gv, us, ms, ievf-ie, use_bt_cont, &
1066 integral_bt_cont, dt, datu, datv, btcl_u, btcl_v)
1073 do j=js,je ;
do i=is-1,ie ; uhbt(i,j) = 0.0 ; ubt(i,j) = 0.0 ;
enddo ;
enddo
1075 do j=js-1,je ;
do i=is,ie ; vhbt(i,j) = 0.0 ; vbt(i,j) = 0.0 ;
enddo ;
enddo
1076 if (cs%visc_rem_u_uh0)
then
1078 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
1079 uhbt(i,j) = uhbt(i,j) + uh0(i,j,k)
1080 ubt(i,j) = ubt(i,j) + wt_u(i,j,k) * u_uh0(i,j,k)
1081 enddo ;
enddo ;
enddo
1083 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
1084 vhbt(i,j) = vhbt(i,j) + vh0(i,j,k)
1085 vbt(i,j) = vbt(i,j) + wt_v(i,j,k) * v_vh0(i,j,k)
1086 enddo ;
enddo ;
enddo
1089 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
1090 uhbt(i,j) = uhbt(i,j) + uh0(i,j,k)
1091 ubt(i,j) = ubt(i,j) + cs%frhatu(i,j,k) * u_uh0(i,j,k)
1092 enddo ;
enddo ;
enddo
1094 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
1095 vhbt(i,j) = vhbt(i,j) + vh0(i,j,k)
1096 vbt(i,j) = vbt(i,j) + cs%frhatv(i,j,k) * v_vh0(i,j,k)
1097 enddo ;
enddo ;
enddo
1099 if ((use_bt_cont .or. integral_bt_cont) .and. cs%adjust_BT_cont)
then
1104 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1105 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
1106 call pass_vector(ubt, vbt, cs%BT_Domain, complete=.false., halo=1+ievf-ie)
1107 call pass_vector(uhbt, vhbt, cs%BT_Domain, complete=.true., halo=1+ievf-ie)
1108 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
1109 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1111 if (integral_bt_cont)
then
1112 call adjust_local_bt_cont_types(ubt, uhbt, vbt, vhbt, btcl_u, btcl_v, &
1113 g, us, ms, halo=1+ievf-ie, dt_baroclinic=dt)
1115 call adjust_local_bt_cont_types(ubt, uhbt, vbt, vhbt, btcl_u, btcl_v, &
1116 g, us, ms, halo=1+ievf-ie)
1119 if (integral_bt_cont)
then
1121 do j=js,je ;
do i=is-1,ie
1122 uhbt0(i,j) = uhbt(i,j) - find_uhbt(dt*ubt(i,j), btcl_u(i,j)) * idt
1125 do j=js-1,je ;
do i=is,ie
1126 vhbt0(i,j) = vhbt(i,j) - find_vhbt(dt*vbt(i,j), btcl_v(i,j)) * idt
1128 elseif (use_bt_cont)
then
1130 do j=js,je ;
do i=is-1,ie
1131 uhbt0(i,j) = uhbt(i,j) - find_uhbt(ubt(i,j), btcl_u(i,j))
1134 do j=js-1,je ;
do i=is,ie
1135 vhbt0(i,j) = vhbt(i,j) - find_vhbt(vbt(i,j), btcl_v(i,j))
1139 do j=js,je ;
do i=is-1,ie
1140 uhbt0(i,j) = uhbt(i,j) - datu(i,j)*ubt(i,j)
1143 do j=js-1,je ;
do i=is,ie
1144 vhbt0(i,j) = vhbt(i,j) - datv(i,j)*vbt(i,j)
1147 if (cs%BT_OBC%apply_u_OBCs)
then
1149 do j=js,je ;
do i=is-1,ie ;
if (obc%segnum_u(i,j) /= obc_none)
then
1151 endif ;
enddo ;
enddo
1153 if (cs%BT_OBC%apply_v_OBCs)
then
1155 do j=js-1,je ;
do i=is,ie ;
if (obc%segnum_v(i,j) /= obc_none)
then
1157 endif ;
enddo ;
enddo
1162 if (integral_bt_cont)
then
1164 do j=jsvf-1,jevf+1 ;
do i=isvf-2,ievf+1
1165 ubt(i,j) = 0.0 ; uhbt(i,j) = 0.0 ; u_accel_bt(i,j) = 0.0
1166 ubt_int(i,j) = 0.0 ; uhbt_int(i,j) = 0.0
1169 do j=jsvf-2,jevf+1 ;
do i=isvf-1,ievf+1
1170 vbt(i,j) = 0.0 ; vhbt(i,j) = 0.0 ; v_accel_bt(i,j) = 0.0
1171 vbt_int(i,j) = 0.0 ; vhbt_int(i,j) = 0.0
1175 do j=jsvf-1,jevf+1 ;
do i=isvf-2,ievf+1
1176 ubt(i,j) = 0.0 ; uhbt(i,j) = 0.0 ; u_accel_bt(i,j) = 0.0
1179 do j=jsvf-2,jevf+1 ;
do i=isvf-1,ievf+1
1180 vbt(i,j) = 0.0 ; vhbt(i,j) = 0.0 ; v_accel_bt(i,j) = 0.0
1184 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
1185 ubt(i,j) = ubt(i,j) + wt_u(i,j,k) * u_in(i,j,k)
1186 enddo ;
enddo ;
enddo
1188 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
1189 vbt(i,j) = vbt(i,j) + wt_v(i,j,k) * v_in(i,j,k)
1190 enddo ;
enddo ;
enddo
1192 do j=js,je ;
do i=is-1,ie
1193 if (abs(ubt(i,j)) < cs%vel_underflow) ubt(i,j) = 0.0
1196 do j=js-1,je ;
do i=is,ie
1197 if (abs(vbt(i,j)) < cs%vel_underflow) vbt(i,j) = 0.0
1200 if (apply_obcs)
then
1201 ubt_first(:,:) = ubt(:,:) ; vbt_first(:,:) = vbt(:,:)
1211 do j=js,je ;
do i=is-1,ie ;
if (g%mask2dCu(i,j) > 0.0)
then
1212 if (cs%nonlin_stress)
then
1213 if (gv%Boussinesq)
then
1214 htot_avg = 0.5*(max(cs%bathyT(i,j)*gv%Z_to_H + eta(i,j), 0.0) + &
1215 max(cs%bathyT(i+1,j)*gv%Z_to_H + eta(i+1,j), 0.0))
1217 htot_avg = 0.5*(eta(i,j) + eta(i+1,j))
1219 if (htot_avg*cs%dy_Cu(i,j) <= 0.0)
then
1221 elseif (integral_bt_cont)
then
1222 cs%IDatu(i,j) = cs%dy_Cu(i,j) / (max(find_duhbt_du(ubt(i,j)*dt, btcl_u(i,j)), &
1223 cs%dy_Cu(i,j)*htot_avg) )
1224 elseif (use_bt_cont)
then
1225 cs%IDatu(i,j) = cs%dy_Cu(i,j) / (max(find_duhbt_du(ubt(i,j), btcl_u(i,j)), &
1226 cs%dy_Cu(i,j)*htot_avg) )
1228 cs%IDatu(i,j) = 1.0 / htot_avg
1232 bt_force_u(i,j) = forces%taux(i,j) * mass_accel_to_z * cs%IDatu(i,j)*visc_rem_u(i,j,1)
1234 bt_force_u(i,j) = 0.0
1235 endif ;
enddo ;
enddo
1237 do j=js-1,je ;
do i=is,ie ;
if (g%mask2dCv(i,j) > 0.0)
then
1238 if (cs%nonlin_stress)
then
1239 if (gv%Boussinesq)
then
1240 htot_avg = 0.5*(max(cs%bathyT(i,j)*gv%Z_to_H + eta(i,j), 0.0) + &
1241 max(cs%bathyT(i,j+1)*gv%Z_to_H + eta(i,j+1), 0.0))
1243 htot_avg = 0.5*(eta(i,j) + eta(i,j+1))
1245 if (htot_avg*cs%dx_Cv(i,j) <= 0.0)
then
1247 elseif (integral_bt_cont)
then
1248 cs%IDatv(i,j) = cs%dx_Cv(i,j) / (max(find_dvhbt_dv(vbt(i,j)*dt, btcl_v(i,j)), &
1249 cs%dx_Cv(i,j)*htot_avg) )
1250 elseif (use_bt_cont)
then
1251 cs%IDatv(i,j) = cs%dx_Cv(i,j) / (max(find_dvhbt_dv(vbt(i,j), btcl_v(i,j)), &
1252 cs%dx_Cv(i,j)*htot_avg) )
1254 cs%IDatv(i,j) = 1.0 / htot_avg
1258 bt_force_v(i,j) = forces%tauy(i,j) * mass_accel_to_z * cs%IDatv(i,j)*visc_rem_v(i,j,1)
1260 bt_force_v(i,j) = 0.0
1261 endif ;
enddo ;
enddo
1262 if (
present(taux_bot) .and.
present(tauy_bot))
then
1263 if (
associated(taux_bot) .and.
associated(tauy_bot))
then
1265 do j=js,je ;
do i=is-1,ie ;
if (g%mask2dCu(i,j) > 0.0)
then
1266 bt_force_u(i,j) = bt_force_u(i,j) - taux_bot(i,j) * mass_to_z * cs%IDatu(i,j)
1267 endif ;
enddo ;
enddo
1269 do j=js-1,je ;
do i=is,ie ;
if (g%mask2dCv(i,j) > 0.0)
then
1270 bt_force_v(i,j) = bt_force_v(i,j) - tauy_bot(i,j) * mass_to_z * cs%IDatv(i,j)
1271 endif ;
enddo ;
enddo
1278 do j=js,je ;
do k=1,nz ;
do i=isq,ieq
1279 bt_force_u(i,j) = bt_force_u(i,j) + wt_u(i,j,k) * bc_accel_u(i,j,k)
1280 enddo ;
enddo ;
enddo
1282 do j=jsq,jeq ;
do k=1,nz ;
do i=is,ie
1283 bt_force_v(i,j) = bt_force_v(i,j) + wt_v(i,j,k) * bc_accel_v(i,j,k)
1284 enddo ;
enddo ;
enddo
1286 if (cs%gradual_BT_ICs)
then
1288 do j=js,je ;
do i=is-1,ie
1289 bt_force_u(i,j) = bt_force_u(i,j) + (ubt(i,j) - cs%ubt_IC(i,j)) * idt
1290 ubt(i,j) = cs%ubt_IC(i,j)
1291 if (abs(ubt(i,j)) < cs%vel_underflow) ubt(i,j) = 0.0
1294 do j=js-1,je ;
do i=is,ie
1295 bt_force_v(i,j) = bt_force_v(i,j) + (vbt(i,j) - cs%vbt_IC(i,j)) * idt
1296 vbt(i,j) = cs%vbt_IC(i,j)
1297 if (abs(vbt(i,j)) < cs%vel_underflow) vbt(i,j) = 0.0
1301 if ((isq > is-1) .or. (jsq > js-1))
then
1304 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1305 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
1306 tmp_u(:,:) = 0.0 ; tmp_v(:,:) = 0.0
1307 do j=js,je ;
do i=isq,ieq ; tmp_u(i,j) = bt_force_u(i,j) ;
enddo ;
enddo
1308 do j=jsq,jeq ;
do i=is,ie ; tmp_v(i,j) = bt_force_v(i,j) ;
enddo ;
enddo
1309 if (nonblock_setup)
then
1310 call start_group_pass(cs%pass_tmp_uv, g%Domain)
1312 call do_group_pass(cs%pass_tmp_uv, g%Domain)
1313 do j=jsd,jed ;
do i=isdb,iedb ; bt_force_u(i,j) = tmp_u(i,j) ;
enddo ;
enddo
1314 do j=jsdb,jedb ;
do i=isd,ied ; bt_force_v(i,j) = tmp_v(i,j) ;
enddo ;
enddo
1316 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
1317 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1320 if (nonblock_setup)
then
1321 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1322 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
1323 call start_group_pass(cs%pass_gtot, cs%BT_Domain)
1324 call start_group_pass(cs%pass_ubt_Cor, g%Domain)
1325 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
1326 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1331 do j=jsvf-1,jevf ;
do i=isvf-1,ievf+1
1332 if (cs%Sadourny)
then
1333 amer(i-1,j) = dcor_u(i-1,j) * q(i-1,j)
1334 bmer(i,j) = dcor_u(i,j) * q(i,j)
1335 cmer(i,j+1) = dcor_u(i,j+1) * q(i,j)
1336 dmer(i-1,j+1) = dcor_u(i-1,j+1) * q(i-1,j)
1338 amer(i-1,j) = dcor_u(i-1,j) * &
1339 ((q(i,j) + q(i-1,j-1)) + q(i-1,j)) / 3.0
1340 bmer(i,j) = dcor_u(i,j) * &
1341 (q(i,j) + (q(i-1,j) + q(i,j-1))) / 3.0
1342 cmer(i,j+1) = dcor_u(i,j+1) * &
1343 (q(i,j) + (q(i-1,j) + q(i,j+1))) / 3.0
1344 dmer(i-1,j+1) = dcor_u(i-1,j+1) * &
1345 ((q(i,j) + q(i-1,j+1)) + q(i-1,j)) / 3.0
1350 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf
1351 if (cs%Sadourny)
then
1352 azon(i,j) = dcor_v(i+1,j) * q(i,j)
1353 bzon(i,j) = dcor_v(i,j) * q(i,j)
1354 czon(i,j) = dcor_v(i,j-1) * q(i,j-1)
1355 dzon(i,j) = dcor_v(i+1,j-1) * q(i,j-1)
1357 azon(i,j) = dcor_v(i+1,j) * &
1358 (q(i,j) + (q(i+1,j) + q(i,j-1))) / 3.0
1359 bzon(i,j) = dcor_v(i,j) * &
1360 (q(i,j) + (q(i-1,j) + q(i,j-1))) / 3.0
1361 czon(i,j) = dcor_v(i,j-1) * &
1362 ((q(i,j) + q(i-1,j-1)) + q(i,j-1)) / 3.0
1363 dzon(i,j) = dcor_v(i+1,j-1) * &
1364 ((q(i,j) + q(i+1,j-1)) + q(i,j-1)) / 3.0
1369 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1370 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
1371 if (nonblock_setup)
then
1372 if ((isq > is-1) .or. (jsq > js-1))
then
1373 call complete_group_pass(cs%pass_tmp_uv, g%Domain)
1374 do j=jsd,jed ;
do i=isdb,iedb ; bt_force_u(i,j) = tmp_u(i,j) ;
enddo ;
enddo
1375 do j=jsdb,jedb ;
do i=isd,ied ; bt_force_v(i,j) = tmp_v(i,j) ;
enddo ;
enddo
1377 call complete_group_pass(cs%pass_gtot, cs%BT_Domain)
1378 call complete_group_pass(cs%pass_ubt_Cor, g%Domain)
1380 call do_group_pass(cs%pass_gtot, cs%BT_Domain)
1381 call do_group_pass(cs%pass_ubt_Cor, g%Domain)
1385 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf+1
1386 if (cs%ua_polarity(i,j) < 0.0)
call swap(gtot_e(i,j), gtot_w(i,j))
1387 if (cs%va_polarity(i,j) < 0.0)
call swap(gtot_n(i,j), gtot_s(i,j))
1391 do j=js,je ;
do i=is-1,ie
1393 ((azon(i,j) * vbt_cor(i+1,j) + czon(i,j) * vbt_cor(i ,j-1)) + &
1394 (bzon(i,j) * vbt_cor(i ,j) + dzon(i,j) * vbt_cor(i+1,j-1)))
1397 do j=js-1,je ;
do i=is,ie
1398 cor_ref_v(i,j) = -1.0 * &
1399 ((amer(i-1,j) * ubt_cor(i-1,j) + cmer(i ,j+1) * ubt_cor(i ,j+1)) + &
1400 (bmer(i ,j) * ubt_cor(i ,j) + dmer(i-1,j+1) * ubt_cor(i-1,j+1)))
1404 if (nonblock_setup)
then
1405 if (.not.use_bt_cont) &
1406 call start_group_pass(cs%pass_Dat_uv, cs%BT_Domain)
1409 call start_group_pass(cs%pass_force_hbt0_Cor_ref, cs%BT_Domain)
1411 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
1412 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1415 do j=js-1,je+1 ;
do i=is-1,ie ; av_rem_u(i,j) = 0.0 ;
enddo ;
enddo
1417 do j=js-1,je ;
do i=is-1,ie+1 ; av_rem_v(i,j) = 0.0 ;
enddo ;
enddo
1419 do j=js,je ;
do k=1,nz ;
do i=is-1,ie
1420 av_rem_u(i,j) = av_rem_u(i,j) + cs%frhatu(i,j,k) * visc_rem_u(i,j,k)
1421 enddo ;
enddo ;
enddo
1423 do j=js-1,je ;
do k=1,nz ;
do i=is,ie
1424 av_rem_v(i,j) = av_rem_v(i,j) + cs%frhatv(i,j,k) * visc_rem_v(i,j,k)
1425 enddo ;
enddo ;
enddo
1426 if (cs%strong_drag)
then
1428 do j=js,je ;
do i=is-1,ie
1429 bt_rem_u(i,j) = g%mask2dCu(i,j) * &
1430 ((nstep * av_rem_u(i,j)) / (1.0 + (nstep-1)*av_rem_u(i,j)))
1433 do j=js-1,je ;
do i=is,ie
1434 bt_rem_v(i,j) = g%mask2dCv(i,j) * &
1435 ((nstep * av_rem_v(i,j)) / (1.0 + (nstep-1)*av_rem_v(i,j)))
1439 do j=js,je ;
do i=is-1,ie
1441 if (g%mask2dCu(i,j) * av_rem_u(i,j) > 0.0) &
1442 bt_rem_u(i,j) = g%mask2dCu(i,j) * (av_rem_u(i,j)**instep)
1445 do j=js-1,je ;
do i=is,ie
1447 if (g%mask2dCv(i,j) * av_rem_v(i,j) > 0.0) &
1448 bt_rem_v(i,j) = g%mask2dCv(i,j) * (av_rem_v(i,j)**instep)
1451 if (cs%linear_wave_drag)
then
1453 do j=js,je ;
do i=is-1,ie ;
if (cs%lin_drag_u(i,j) > 0.0)
then
1454 htot = 0.5 * (eta(i,j) + eta(i+1,j))
1455 if (gv%Boussinesq) &
1456 htot = htot + 0.5*gv%Z_to_H * (cs%bathyT(i,j) + cs%bathyT(i+1,j))
1457 bt_rem_u(i,j) = bt_rem_u(i,j) * (htot / (htot + cs%lin_drag_u(i,j) * dtbt))
1459 rayleigh_u(i,j) = cs%lin_drag_u(i,j) / htot
1460 endif ;
enddo ;
enddo
1462 do j=js-1,je ;
do i=is,ie ;
if (cs%lin_drag_v(i,j) > 0.0)
then
1463 htot = 0.5 * (eta(i,j) + eta(i,j+1))
1464 if (gv%Boussinesq) &
1465 htot = htot + 0.5*gv%Z_to_H * (cs%bathyT(i,j) + cs%bathyT(i+1,j+1))
1466 bt_rem_v(i,j) = bt_rem_v(i,j) * (htot / (htot + cs%lin_drag_v(i,j) * dtbt))
1468 rayleigh_v(i,j) = cs%lin_drag_v(i,j) / htot
1469 endif ;
enddo ;
enddo
1473 if (find_etaav)
then
1475 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf+1
1476 eta_sum(i,j) = 0.0 ; eta_wtd(i,j) = 0.0
1480 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf+1
1485 do j=jsvf-1,jevf+1 ;
do i=isvf-1,ievf
1486 ubt_sum(i,j) = 0.0 ; uhbt_sum(i,j) = 0.0
1487 pfu_bt_sum(i,j) = 0.0 ; coru_bt_sum(i,j) = 0.0
1488 ubt_wtd(i,j) = 0.0 ; ubt_trans(i,j) = 0.0
1491 do j=jsvf-1,jevf ;
do i=isvf-1,ievf+1
1492 vbt_sum(i,j) = 0.0 ; vhbt_sum(i,j) = 0.0
1493 pfv_bt_sum(i,j) = 0.0 ; corv_bt_sum(i,j) = 0.0
1494 vbt_wtd(i,j) = 0.0 ; vbt_trans(i,j) = 0.0
1499 do j=jsvf-1,jevf+1;
do i=isvf-1,ievf+1 ; eta_src(i,j) = 0.0 ;
enddo ;
enddo
1500 if (cs%bound_BT_corr)
then ;
if ((use_bt_cont.or.integral_bt_cont) .and. cs%BT_cont_bounds)
then
1501 do j=js,je ;
do i=is,ie ;
if (g%mask2dT(i,j) > 0.0)
then
1502 if (cs%eta_cor(i,j) > 0.0)
then
1505 if (integral_bt_cont)
then
1506 uint_cor = g%dxT(i,j) * cs%maxCFL_BT_cont
1507 vint_cor = g%dyT(i,j) * cs%maxCFL_BT_cont
1508 eta_cor_max = (cs%IareaT(i,j) * &
1509 (((find_uhbt(uint_cor, btcl_u(i,j)) + dt*uhbt0(i,j)) - &
1510 (find_uhbt(-uint_cor, btcl_u(i-1,j)) + dt*uhbt0(i-1,j))) + &
1511 ((find_vhbt(vint_cor, btcl_v(i,j)) + dt*vhbt0(i,j)) - &
1512 (find_vhbt(-vint_cor, btcl_v(i,j-1)) + dt*vhbt0(i,j-1))) ))
1514 u_max_cor = g%dxT(i,j) * (cs%maxCFL_BT_cont*idt)
1515 v_max_cor = g%dyT(i,j) * (cs%maxCFL_BT_cont*idt)
1516 eta_cor_max = dt * (cs%IareaT(i,j) * &
1517 (((find_uhbt(u_max_cor, btcl_u(i,j)) + uhbt0(i,j)) - &
1518 (find_uhbt(-u_max_cor, btcl_u(i-1,j)) + uhbt0(i-1,j))) + &
1519 ((find_vhbt(v_max_cor, btcl_v(i,j)) + vhbt0(i,j)) - &
1520 (find_vhbt(-v_max_cor, btcl_v(i,j-1)) + vhbt0(i,j-1))) ))
1522 cs%eta_cor(i,j) = min(cs%eta_cor(i,j), max(0.0, eta_cor_max))
1526 if (gv%Boussinesq) htot = cs%bathyT(i,j)*gv%Z_to_H + eta(i,j)
1528 cs%eta_cor(i,j) = max(cs%eta_cor(i,j), -max(0.0,htot))
1530 endif ;
enddo ;
enddo
1531 else ;
do j=js,je ;
do i=is,ie
1532 if (abs(cs%eta_cor(i,j)) > dt*cs%eta_cor_bound(i,j)) &
1533 cs%eta_cor(i,j) = sign(dt*cs%eta_cor_bound(i,j), cs%eta_cor(i,j))
1534 enddo ;
enddo ;
endif ;
endif
1536 do j=js,je ;
do i=is,ie
1537 eta_src(i,j) = g%mask2dT(i,j) * (instep * cs%eta_cor(i,j))
1541 if (cs%dynamic_psurf)
then
1542 ice_is_rigid = (
associated(forces%rigidity_ice_u) .and. &
1543 associated(forces%rigidity_ice_v))
1544 h_min_dyn = gv%Z_to_H * cs%Dmin_dyn_psurf
1545 if (ice_is_rigid .and. use_bt_cont) &
1546 call bt_cont_to_face_areas(bt_cont, datu, datv, g, us, ms, 0, .true.)
1547 if (ice_is_rigid)
then
1549 do j=js,je ;
do i=is,ie
1555 idt_max2 = 0.5 * (dgeo_de * (1.0 + 2.0*bebt)) * (g%IareaT(i,j) * &
1556 ((gtot_e(i,j) * (datu(i,j)*g%IdxCu(i,j)) + &
1557 gtot_w(i,j) * (datu(i-1,j)*g%IdxCu(i-1,j))) + &
1558 (gtot_n(i,j) * (datv(i,j)*g%IdyCv(i,j)) + &
1559 gtot_s(i,j) * (datv(i,j-1)*g%IdyCv(i,j-1)))) + &
1560 ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
1561 (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2)) * cs%BT_Coriolis_scale**2 )
1562 h_eff_dx2 = max(h_min_dyn * ((g%IdxT(i,j))**2 + (g%IdyT(i,j))**2), &
1564 ((datu(i,j)*g%IdxCu(i,j) + datu(i-1,j)*g%IdxCu(i-1,j)) + &
1565 (datv(i,j)*g%IdyCv(i,j) + datv(i,j-1)*g%IdyCv(i,j-1)) ) )
1566 dyn_coef_max = cs%const_dyn_psurf * max(0.0, 1.0 - dtbt**2 * idt_max2) / &
1567 (dtbt**2 * h_eff_dx2)
1570 ice_strength = ((forces%rigidity_ice_u(i,j) + forces%rigidity_ice_u(i-1,j)) + &
1571 (forces%rigidity_ice_v(i,j) + forces%rigidity_ice_v(i,j-1))) / &
1572 (cs%ice_strength_length**2 * dtbt)
1575 dyn_coef_eta(i,j) = min(dyn_coef_max, ice_strength * gv%H_to_Z)
1576 enddo ;
enddo ;
endif
1579 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1580 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
1581 if (nonblock_setup)
then
1582 call start_group_pass(cs%pass_eta_bt_rem, cs%BT_Domain)
1585 call do_group_pass(cs%pass_eta_bt_rem, cs%BT_Domain)
1586 if (.not.use_bt_cont) &
1587 call do_group_pass(cs%pass_Dat_uv, cs%BT_Domain)
1588 call do_group_pass(cs%pass_force_hbt0_Cor_ref, cs%BT_Domain)
1590 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
1591 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1594 if (nonblock_setup)
then
1595 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1596 if (id_clock_pass_pre > 0)
call cpu_clock_begin(id_clock_pass_pre)
1598 if (.not.use_bt_cont)
call complete_group_pass(cs%pass_Dat_uv, cs%BT_Domain)
1599 call complete_group_pass(cs%pass_force_hbt0_Cor_ref, cs%BT_Domain)
1600 call complete_group_pass(cs%pass_eta_bt_rem, cs%BT_Domain)
1602 if (id_clock_pass_pre > 0)
call cpu_clock_end(id_clock_pass_pre)
1603 if (id_clock_calc_pre > 0)
call cpu_clock_begin(id_clock_calc_pre)
1607 call uvchksum(
"BT [uv]hbt", uhbt, vhbt, cs%debug_BT_HI, haloshift=0, &
1608 scale=us%s_to_T*us%L_to_m**2*gv%H_to_m)
1609 call uvchksum(
"BT Initial [uv]bt", ubt, vbt, cs%debug_BT_HI, haloshift=0, scale=us%L_T_to_m_s)
1610 call hchksum(eta,
"BT Initial eta", cs%debug_BT_HI, haloshift=0, scale=gv%H_to_m)
1611 call uvchksum(
"BT BT_force_[uv]", bt_force_u, bt_force_v, &
1612 cs%debug_BT_HI, haloshift=0, scale=us%L_T2_to_m_s2)
1613 if (interp_eta_pf)
then
1614 call hchksum(eta_pf_1,
"BT eta_PF_1",cs%debug_BT_HI,haloshift=0, scale=gv%H_to_m)
1615 call hchksum(d_eta_pf,
"BT d_eta_PF",cs%debug_BT_HI,haloshift=0, scale=gv%H_to_m)
1617 call hchksum(eta_pf,
"BT eta_PF",cs%debug_BT_HI,haloshift=0, scale=gv%H_to_m)
1618 call hchksum(eta_pf_in,
"BT eta_PF_in",g%HI,haloshift=0, scale=gv%H_to_m)
1620 call uvchksum(
"BT Cor_ref_[uv]", cor_ref_u, cor_ref_v, cs%debug_BT_HI, haloshift=0, scale=us%L_T2_to_m_s2)
1621 call uvchksum(
"BT [uv]hbt0", uhbt0, vhbt0, cs%debug_BT_HI, haloshift=0, &
1622 scale=us%L_to_m**2*us%s_to_T*gv%H_to_m)
1623 if (.not. use_bt_cont)
then
1624 call uvchksum(
"BT Dat[uv]", datu, datv, cs%debug_BT_HI, haloshift=1, scale=us%L_to_m*gv%H_to_m)
1626 call uvchksum(
"BT wt_[uv]", wt_u, wt_v, g%HI, haloshift=0, &
1627 symmetric=.true., omit_corners=.true., scalar_pair=.true.)
1628 call uvchksum(
"BT frhat[uv]", cs%frhatu, cs%frhatv, g%HI, haloshift=0, &
1629 symmetric=.true., omit_corners=.true., scalar_pair=.true.)
1630 call uvchksum(
"BT bc_accel_[uv]", bc_accel_u, bc_accel_v, g%HI, haloshift=0, scale=us%L_T2_to_m_s2)
1631 call uvchksum(
"BT IDat[uv]", cs%IDatu, cs%IDatv, g%HI, haloshift=0, &
1632 scale=us%m_to_Z, scalar_pair=.true.)
1633 call uvchksum(
"BT visc_rem_[uv]", visc_rem_u, visc_rem_v, g%HI, &
1634 haloshift=1, scalar_pair=.true.)
1637 if (cs%id_ubtdt > 0)
then
1638 do j=js-1,je+1 ;
do i=is-1,ie
1639 ubt_st(i,j) = ubt(i,j)
1642 if (cs%id_vbtdt > 0)
then
1643 do j=js-1,je ;
do i=is-1,ie+1
1644 vbt_st(i,j) = vbt(i,j)
1648 if (query_averaging_enabled(cs%diag))
then
1649 if (cs%id_eta_st > 0)
call post_data(cs%id_eta_st, eta(isd:ied,jsd:jed), cs%diag)
1650 if (cs%id_ubt_st > 0)
call post_data(cs%id_ubt_st, ubt(isdb:iedb,jsd:jed), cs%diag)
1651 if (cs%id_vbt_st > 0)
call post_data(cs%id_vbt_st, vbt(isd:ied,jsdb:jedb), cs%diag)
1654 if (id_clock_calc_pre > 0)
call cpu_clock_end(id_clock_calc_pre)
1655 if (id_clock_calc > 0)
call cpu_clock_begin(id_clock_calc)
1657 if (project_velocity)
then ; eta_pf_bt => eta ;
else ; eta_pf_bt => eta_pred ;
endif
1659 if (cs%dt_bt_filter >= 0.0)
then
1660 dt_filt = 0.5 * max(0.0, min(cs%dt_bt_filter, 2.0*dt))
1662 dt_filt = 0.5 * max(0.0, dt * min(-cs%dt_bt_filter, 2.0))
1664 nfilter = ceiling(dt_filt / dtbt)
1666 if (nstep+nfilter==0 )
call mom_error(fatal, &
1667 "btstep: number of barotropic step (nstep+nfilter) is 0")
1670 sum_wt_vel = 0.0 ; sum_wt_eta = 0.0 ; sum_wt_accel = 0.0 ; sum_wt_trans = 0.0
1671 allocate(wt_vel(nstep+nfilter)) ;
allocate(wt_eta(nstep+nfilter))
1672 allocate(wt_trans(nstep+nfilter+1)) ;
allocate(wt_accel(nstep+nfilter+1))
1673 allocate(wt_accel2(nstep+nfilter+1))
1674 do n=1,nstep+nfilter
1678 if ( (n==nstep) .or. (dt_filt - abs(n-nstep)*dtbt >= 0.0))
then
1679 wt_vel(n) = 1.0 ; wt_eta(n) = 1.0
1680 elseif (dtbt + dt_filt - abs(n-nstep)*dtbt > 0.0)
then
1681 wt_vel(n) = 1.0 + (dt_filt / dtbt) - abs(n-nstep) ; wt_eta(n) = wt_vel(n)
1683 wt_vel(n) = 0.0 ; wt_eta(n) = 0.0
1690 sum_wt_vel = sum_wt_vel + wt_vel(n) ; sum_wt_eta = sum_wt_eta + wt_eta(n)
1692 wt_trans(nstep+nfilter+1) = 0.0 ; wt_accel(nstep+nfilter+1) = 0.0
1693 do n=nstep+nfilter,1,-1
1694 wt_trans(n) = wt_trans(n+1) + wt_eta(n)
1695 wt_accel(n) = wt_accel(n+1) + wt_vel(n)
1696 sum_wt_accel = sum_wt_accel + wt_accel(n) ; sum_wt_trans = sum_wt_trans + wt_trans(n)
1699 i_sum_wt_vel = 1.0 / sum_wt_vel ; i_sum_wt_accel = 1.0 / sum_wt_accel
1700 i_sum_wt_eta = 1.0 / sum_wt_eta ; i_sum_wt_trans = 1.0 / sum_wt_trans
1701 do n=1,nstep+nfilter
1702 wt_vel(n) = wt_vel(n) * i_sum_wt_vel
1703 if (cs%answers_2018)
then
1704 wt_accel2(n) = wt_accel(n)
1707 wt_accel2(n) = wt_accel(n) * i_sum_wt_accel
1708 wt_trans(n) = wt_trans(n) * i_sum_wt_trans
1710 wt_accel(n) = wt_accel(n) * i_sum_wt_accel
1711 wt_eta(n) = wt_eta(n) * i_sum_wt_eta
1714 sum_wt_vel = 0.0 ; sum_wt_eta = 0.0 ; sum_wt_accel = 0.0 ; sum_wt_trans = 0.0
1717 isv=is ; iev=ie ; jsv=js ; jev=je
1718 do n=1,nstep+nfilter
1720 sum_wt_vel = sum_wt_vel + wt_vel(n)
1721 sum_wt_eta = sum_wt_eta + wt_eta(n)
1722 sum_wt_accel = sum_wt_accel + wt_accel2(n)
1723 sum_wt_trans = sum_wt_trans + wt_trans(n)
1725 if (cs%clip_velocity)
then
1726 do j=jsv,jev ;
do i=isv-1,iev
1727 if ((ubt(i,j) * (dt * g%dy_Cu(i,j))) * g%IareaT(i+1,j) < -cs%CFL_trunc)
then
1729 ubt(i,j) = (-0.95*cs%CFL_trunc) * (g%areaT(i+1,j) / (dt * g%dy_Cu(i,j)))
1730 elseif ((ubt(i,j) * (dt * g%dy_Cu(i,j))) * g%IareaT(i,j) > cs%CFL_trunc)
then
1732 ubt(i,j) = (0.95*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dy_Cu(i,j)))
1735 do j=jsv-1,jev ;
do i=isv,iev
1736 if ((vbt(i,j) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j+1) < -cs%CFL_trunc)
then
1738 vbt(i,j) = (-0.9*cs%CFL_trunc) * (g%areaT(i,j+1) / (dt * g%dx_Cv(i,j)))
1739 elseif ((vbt(i,j) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j) > cs%CFL_trunc)
then
1741 vbt(i,j) = (0.9*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dx_Cv(i,j)))
1746 if ((iev - stencil < ie) .or. (jev - stencil < je))
then
1747 if (id_clock_calc > 0)
call cpu_clock_end(id_clock_calc)
1748 call do_group_pass(cs%pass_eta_ubt, cs%BT_Domain, clock=id_clock_pass_step)
1749 isv = isvf ; iev = ievf ; jsv = jsvf ; jev = jevf
1750 if (id_clock_calc > 0)
call cpu_clock_begin(id_clock_calc)
1752 isv = isv+stencil ; iev = iev-stencil
1753 jsv = jsv+stencil ; jev = jev-stencil
1756 if ((.not.use_bt_cont) .and. cs%Nonlinear_continuity .and. &
1757 (cs%Nonlin_cont_update_period > 0))
then
1758 if ((n>1) .and. (mod(n-1,cs%Nonlin_cont_update_period) == 0)) &
1759 call find_face_areas(datu, datv, g, gv, us, cs, ms, eta, 1+iev-ie)
1762 if (integral_bt_cont)
then
1764 do j=jsv-1,jev+1 ;
do i=isv-2,iev+1
1765 ubt_int_prev(i,j) = ubt_int(i,j) ; uhbt_int_prev(i,j) = uhbt_int(i,j)
1768 do j=jsv-2,jev+1 ;
do i=isv-1,iev+1
1769 vbt_int_prev(i,j) = vbt_int(i,j) ; vhbt_int_prev(i,j) = vhbt_int(i,j)
1774 if (cs%dynamic_psurf .or. .not.project_velocity)
then
1775 if (integral_bt_cont)
then
1777 do j=jsv-1,jev+1 ;
do i=isv-2,iev+1
1778 uhbt_int(i,j) = find_uhbt(ubt_int(i,j) + dtbt*ubt(i,j), btcl_u(i,j)) + n*dtbt*uhbt0(i,j)
1782 do j=jsv-2,jev+1 ;
do i=isv-1,iev+1
1783 vhbt_int(i,j) = find_vhbt(vbt_int(i,j) + dtbt*vbt(i,j), btcl_v(i,j)) + n*dtbt*vhbt0(i,j)
1786 do j=jsv-1,jev+1 ;
do i=isv-1,iev+1
1787 eta_pred(i,j) = (eta_ic(i,j) + n*eta_src(i,j)) + cs%IareaT(i,j) * &
1788 ((uhbt_int(i-1,j) - uhbt_int(i,j)) + (vhbt_int(i,j-1) - vhbt_int(i,j)))
1790 elseif (use_bt_cont)
then
1792 do j=jsv-1,jev+1 ;
do i=isv-2,iev+1
1793 uhbt(i,j) = find_uhbt(ubt(i,j), btcl_u(i,j)) + uhbt0(i,j)
1796 do j=jsv-2,jev+1 ;
do i=isv-1,iev+1
1797 vhbt(i,j) = find_vhbt(vbt(i,j), btcl_v(i,j)) + vhbt0(i,j)
1800 do j=jsv-1,jev+1 ;
do i=isv-1,iev+1
1801 eta_pred(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * cs%IareaT(i,j)) * &
1802 ((uhbt(i-1,j) - uhbt(i,j)) + (vhbt(i,j-1) - vhbt(i,j)))
1806 do j=jsv-1,jev+1 ;
do i=isv-1,iev+1
1807 eta_pred(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * cs%IareaT(i,j)) * &
1808 (((datu(i-1,j)*ubt(i-1,j) + uhbt0(i-1,j)) - &
1809 (datu(i,j)*ubt(i,j) + uhbt0(i,j))) + &
1810 ((datv(i,j-1)*vbt(i,j-1) + vhbt0(i,j-1)) - &
1811 (datv(i,j)*vbt(i,j) + vhbt0(i,j))))
1815 if (cs%dynamic_psurf)
then
1817 do j=jsv-1,jev+1 ;
do i=isv-1,iev+1
1818 p_surf_dyn(i,j) = dyn_coef_eta(i,j) * (eta_pred(i,j) - eta(i,j))
1826 if (find_etaav)
then
1828 do j=js,je ;
do i=is,ie
1829 eta_sum(i,j) = eta_sum(i,j) + wt_accel2(n) * eta_pf_bt(i,j)
1834 if (interp_eta_pf)
then
1837 do j=jsv-1,jev+1 ;
do i=isv-1,iev+1
1838 eta_pf(i,j) = eta_pf_1(i,j) + wt_end*d_eta_pf(i,j)
1842 if (apply_obc_flather .or. apply_obc_open)
then
1844 do j=jsv,jev ;
do i=isv-2,iev+1
1845 ubt_old(i,j) = ubt(i,j)
1848 do j=jsv-2,jev+1 ;
do i=isv,iev
1849 vbt_old(i,j) = vbt(i,j)
1853 if (apply_obcs)
then
1854 if (mod(n+g%first_direction,2)==1)
then
1860 if (cs%BT_OBC%apply_u_OBCs)
then
1862 do j=jsv-joff,jev+joff ;
do i=isv-1,iev
1863 ubt_prev(i,j) = ubt(i,j) ; uhbt_prev(i,j) = uhbt(i,j)
1864 ubt_sum_prev(i,j) = ubt_sum(i,j) ; uhbt_sum_prev(i,j) = uhbt_sum(i,j) ; ubt_wtd_prev(i,j) = ubt_wtd(i,j)
1868 if (cs%BT_OBC%apply_v_OBCs)
then
1870 do j=jsv-1,jev ;
do i=isv-ioff,iev+ioff
1871 vbt_prev(i,j) = vbt(i,j) ; vhbt_prev(i,j) = vhbt(i,j)
1872 vbt_sum_prev(i,j) = vbt_sum(i,j) ; vhbt_sum_prev(i,j) = vhbt_sum(i,j) ; vbt_wtd_prev(i,j) = vbt_wtd(i,j)
1877 if (mod(n+g%first_direction,2)==1)
then
1880 do j=jsv-1,jev ;
do i=isv-1,iev+1
1881 cor_v(i,j) = -1.0*((amer(i-1,j) * ubt(i-1,j) + cmer(i,j+1) * ubt(i,j+1)) + &
1882 (bmer(i,j) * ubt(i,j) + dmer(i-1,j+1) * ubt(i-1,j+1))) - cor_ref_v(i,j)
1883 pfv(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_n(i,j) - &
1884 (eta_pf_bt(i,j+1)-eta_pf(i,j+1))*gtot_s(i,j+1)) * &
1885 dgeo_de * cs%IdyCv(i,j)
1888 if (cs%dynamic_psurf)
then
1890 do j=jsv-1,jev ;
do i=isv-1,iev+1
1891 pfv(i,j) = pfv(i,j) + (p_surf_dyn(i,j) - p_surf_dyn(i,j+1)) * cs%IdyCv(i,j)
1896 if (cs%BT_OBC%apply_v_OBCs)
then
1898 do j=jsv-1,jev ;
do i=isv-1,iev+1 ;
if (obc%segnum_v(i,j) /= obc_none)
then
1900 endif ;
enddo ;
enddo
1904 do j=jsv-1,jev ;
do i=isv-1,iev+1
1906 vbt(i,j) = bt_rem_v(i,j) * (vbt(i,j) + &
1907 dtbt * ((bt_force_v(i,j) + cor_v(i,j)) + pfv(i,j)))
1908 vbt_trans(i,j) = trans_wt1*vbt(i,j) + trans_wt2*vel_prev
1910 if (cs%linear_wave_drag)
then
1911 v_accel_bt(i,j) = v_accel_bt(i,j) + wt_accel(n) * &
1912 ((cor_v(i,j) + pfv(i,j)) - vbt(i,j)*rayleigh_v(i,j))
1914 v_accel_bt(i,j) = v_accel_bt(i,j) + wt_accel(n) * (cor_v(i,j) + pfv(i,j))
1918 if (integral_bt_cont)
then
1920 do j=jsv-1,jev ;
do i=isv-1,iev+1
1921 vbt_int(i,j) = vbt_int(i,j) + dtbt * vbt_trans(i,j)
1922 vhbt_int(i,j) = find_vhbt(vbt_int(i,j), btcl_v(i,j)) + n*dtbt*vhbt0(i,j)
1924 vhbt(i,j) = (vhbt_int(i,j) - vhbt_int_prev(i,j)) * idtbt
1926 elseif (use_bt_cont)
then
1928 do j=jsv-1,jev ;
do i=isv-1,iev+1
1929 vhbt(i,j) = find_vhbt(vbt_trans(i,j), btcl_v(i,j)) + vhbt0(i,j)
1934 do j=jsv-1,jev ;
do i=isv-1,iev+1
1935 vhbt(i,j) = datv(i,j)*vbt_trans(i,j) + vhbt0(i,j)
1939 if (cs%BT_OBC%apply_v_OBCs)
then
1941 do j=jsv-1,jev ;
do i=isv-1,iev+1 ;
if (obc%segnum_v(i,j) /= obc_none)
then
1942 vbt(i,j) = vbt_prev(i,j) ; vhbt(i,j) = vhbt_prev(i,j)
1943 endif ;
enddo ;
enddo
1947 do j=jsv,jev ;
do i=isv-1,iev
1948 cor_u(i,j) = ((azon(i,j) * vbt(i+1,j) + czon(i,j) * vbt(i,j-1)) + &
1949 (bzon(i,j) * vbt(i,j) + dzon(i,j) * vbt(i+1,j-1))) - &
1951 pfu(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_e(i,j) - &
1952 (eta_pf_bt(i+1,j)-eta_pf(i+1,j))*gtot_w(i+1,j)) * &
1953 dgeo_de * cs%IdxCu(i,j)
1957 if (cs%dynamic_psurf)
then
1959 do j=jsv,jev ;
do i=isv-1,iev
1960 pfu(i,j) = pfu(i,j) + (p_surf_dyn(i,j) - p_surf_dyn(i+1,j)) * cs%IdxCu(i,j)
1965 if (cs%BT_OBC%apply_u_OBCs)
then
1967 do j=jsv,jev ;
do i=isv-1,iev ;
if (obc%segnum_u(i,j) /= obc_none)
then
1969 endif ;
enddo ;
enddo
1973 do j=jsv,jev ;
do i=isv-1,iev
1975 ubt(i,j) = bt_rem_u(i,j) * (ubt(i,j) + &
1976 dtbt * ((bt_force_u(i,j) + cor_u(i,j)) + pfu(i,j)))
1977 if (abs(ubt(i,j)) < cs%vel_underflow) ubt(i,j) = 0.0
1978 ubt_trans(i,j) = trans_wt1*ubt(i,j) + trans_wt2*vel_prev
1980 if (cs%linear_wave_drag)
then
1981 u_accel_bt(i,j) = u_accel_bt(i,j) + wt_accel(n) * &
1982 ((cor_u(i,j) + pfu(i,j)) - ubt(i,j)*rayleigh_u(i,j))
1984 u_accel_bt(i,j) = u_accel_bt(i,j) + wt_accel(n) * (cor_u(i,j) + pfu(i,j))
1989 if (integral_bt_cont)
then
1991 do j=jsv,jev ;
do i=isv-1,iev
1992 ubt_int(i,j) = ubt_int(i,j) + dtbt * ubt_trans(i,j)
1993 uhbt_int(i,j) = find_uhbt(ubt_int(i,j), btcl_u(i,j)) + n*dtbt*uhbt0(i,j)
1995 uhbt(i,j) = (uhbt_int(i,j) - uhbt_int_prev(i,j)) * idtbt
1997 elseif (use_bt_cont)
then
1999 do j=jsv,jev ;
do i=isv-1,iev
2000 uhbt(i,j) = find_uhbt(ubt_trans(i,j), btcl_u(i,j)) + uhbt0(i,j)
2004 do j=jsv,jev ;
do i=isv-1,iev
2005 uhbt(i,j) = datu(i,j)*ubt_trans(i,j) + uhbt0(i,j)
2008 if (cs%BT_OBC%apply_u_OBCs)
then
2010 do j=jsv,jev ;
do i=isv-1,iev ;
if (obc%segnum_u(i,j) /= obc_none)
then
2011 ubt(i,j) = ubt_prev(i,j) ; uhbt(i,j) = uhbt_prev(i,j)
2012 endif ;
enddo ;
enddo
2017 do j=jsv-1,jev+1 ;
do i=isv-1,iev
2018 cor_u(i,j) = ((azon(i,j) * vbt(i+1,j) + czon(i,j) * vbt(i,j-1)) + &
2019 (bzon(i,j) * vbt(i,j) + dzon(i,j) * vbt(i+1,j-1))) - &
2021 pfu(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_e(i,j) - &
2022 (eta_pf_bt(i+1,j)-eta_pf(i+1,j))*gtot_w(i+1,j)) * &
2023 dgeo_de * cs%IdxCu(i,j)
2027 if (cs%dynamic_psurf)
then
2029 do j=jsv-1,jev+1 ;
do i=isv-1,iev
2030 pfu(i,j) = pfu(i,j) + (p_surf_dyn(i,j) - p_surf_dyn(i+1,j)) * cs%IdxCu(i,j)
2035 if (cs%BT_OBC%apply_u_OBCs)
then
2037 do j=jsv,jev ;
do i=isv-1,iev ;
if (obc%segnum_u(i,j) /= obc_none)
then
2039 endif ;
enddo ;
enddo
2043 do j=jsv-1,jev+1 ;
do i=isv-1,iev
2045 ubt(i,j) = bt_rem_u(i,j) * (ubt(i,j) + &
2046 dtbt * ((bt_force_u(i,j) + cor_u(i,j)) + pfu(i,j)))
2047 if (abs(ubt(i,j)) < cs%vel_underflow) ubt(i,j) = 0.0
2048 ubt_trans(i,j) = trans_wt1*ubt(i,j) + trans_wt2*vel_prev
2050 if (cs%linear_wave_drag)
then
2051 u_accel_bt(i,j) = u_accel_bt(i,j) + wt_accel(n) * &
2052 ((cor_u(i,j) + pfu(i,j)) - ubt(i,j)*rayleigh_u(i,j))
2054 u_accel_bt(i,j) = u_accel_bt(i,j) + wt_accel(n) * (cor_u(i,j) + pfu(i,j))
2058 if (integral_bt_cont)
then
2060 do j=jsv-1,jev+1 ;
do i=isv-1,iev
2061 ubt_int(i,j) = ubt_int(i,j) + dtbt * ubt_trans(i,j)
2062 uhbt_int(i,j) = find_uhbt(ubt_int(i,j), btcl_u(i,j)) + n*dtbt*uhbt0(i,j)
2064 uhbt(i,j) = (uhbt_int(i,j) - uhbt_int_prev(i,j)) * idtbt
2066 elseif (use_bt_cont)
then
2068 do j=jsv-1,jev+1 ;
do i=isv-1,iev
2069 uhbt(i,j) = find_uhbt(ubt_trans(i,j), btcl_u(i,j)) + uhbt0(i,j)
2074 do j=jsv-1,jev+1 ;
do i=isv-1,iev
2075 uhbt(i,j) = datu(i,j)*ubt_trans(i,j) + uhbt0(i,j)
2079 if (cs%BT_OBC%apply_u_OBCs)
then
2081 do j=jsv-1,jev+1 ;
do i=isv-1,iev ;
if (obc%segnum_u(i,j) /= obc_none)
then
2082 ubt(i,j) = ubt_prev(i,j) ; uhbt(i,j) = uhbt_prev(i,j)
2083 endif ;
enddo ;
enddo
2087 if (cs%use_old_coriolis_bracket_bug)
then
2089 do j=jsv-1,jev ;
do i=isv,iev
2090 cor_v(i,j) = -1.0*((amer(i-1,j) * ubt(i-1,j) + bmer(i,j) * ubt(i,j)) + &
2091 (cmer(i,j+1) * ubt(i,j+1) + dmer(i-1,j+1) * ubt(i-1,j+1))) - cor_ref_v(i,j)
2092 pfv(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_n(i,j) - &
2093 (eta_pf_bt(i,j+1)-eta_pf(i,j+1))*gtot_s(i,j+1)) * &
2094 dgeo_de * cs%IdyCv(i,j)
2099 do j=jsv-1,jev ;
do i=isv,iev
2100 cor_v(i,j) = -1.0*((amer(i-1,j) * ubt(i-1,j) + cmer(i,j+1) * ubt(i,j+1)) + &
2101 (bmer(i,j) * ubt(i,j) + dmer(i-1,j+1) * ubt(i-1,j+1))) - cor_ref_v(i,j)
2102 pfv(i,j) = ((eta_pf_bt(i,j)-eta_pf(i,j))*gtot_n(i,j) - &
2103 (eta_pf_bt(i,j+1)-eta_pf(i,j+1))*gtot_s(i,j+1)) * &
2104 dgeo_de * cs%IdyCv(i,j)
2109 if (cs%dynamic_psurf)
then
2111 do j=jsv-1,jev ;
do i=isv,iev
2112 pfv(i,j) = pfv(i,j) + (p_surf_dyn(i,j) - p_surf_dyn(i,j+1)) * cs%IdyCv(i,j)
2117 if (cs%BT_OBC%apply_v_OBCs)
then
2119 do j=jsv-1,jev ;
do i=isv-1,iev+1 ;
if (obc%segnum_v(i,j) /= obc_none)
then
2121 endif ;
enddo ;
enddo
2125 do j=jsv-1,jev ;
do i=isv,iev
2127 vbt(i,j) = bt_rem_v(i,j) * (vbt(i,j) + &
2128 dtbt * ((bt_force_v(i,j) + cor_v(i,j)) + pfv(i,j)))
2129 if (abs(vbt(i,j)) < cs%vel_underflow) vbt(i,j) = 0.0
2130 vbt_trans(i,j) = trans_wt1*vbt(i,j) + trans_wt2*vel_prev
2132 if (cs%linear_wave_drag)
then
2133 v_accel_bt(i,j) = v_accel_bt(i,j) + wt_accel(n) * &
2134 ((cor_v(i,j) + pfv(i,j)) - vbt(i,j)*rayleigh_v(i,j))
2136 v_accel_bt(i,j) = v_accel_bt(i,j) + wt_accel(n) * (cor_v(i,j) + pfv(i,j))
2140 if (integral_bt_cont)
then
2142 do j=jsv-1,jev ;
do i=isv,iev
2143 vbt_int(i,j) = vbt_int(i,j) + dtbt * vbt_trans(i,j)
2144 vhbt_int(i,j) = find_vhbt(vbt_int(i,j), btcl_v(i,j)) + n*dtbt*vhbt0(i,j)
2146 vhbt(i,j) = (vhbt_int(i,j) - vhbt_int_prev(i,j)) * idtbt
2148 elseif (use_bt_cont)
then
2150 do j=jsv-1,jev ;
do i=isv,iev
2151 vhbt(i,j) = find_vhbt(vbt_trans(i,j), btcl_v(i,j)) + vhbt0(i,j)
2155 do j=jsv-1,jev ;
do i=isv,iev
2156 vhbt(i,j) = datv(i,j)*vbt_trans(i,j) + vhbt0(i,j)
2159 if (cs%BT_OBC%apply_v_OBCs)
then
2161 do j=jsv-1,jev ;
do i=isv,iev ;
if (obc%segnum_v(i,j) /= obc_none)
then
2162 vbt(i,j) = vbt_prev(i,j); vhbt(i,j) = vhbt_prev(i,j)
2163 endif ;
enddo ;
enddo
2168 if (cs%debug_bt)
then
2169 write(mesg,
'("BT vel update ",I4)') n
2170 call uvchksum(trim(mesg)//
" PF[uv]", pfu, pfv, cs%debug_BT_HI, haloshift=iev-ie, &
2171 scale=us%L_T_to_m_s*us%s_to_T)
2172 call uvchksum(trim(mesg)//
" Cor_[uv]", cor_u, cor_v, cs%debug_BT_HI, haloshift=iev-ie, &
2173 scale=us%L_T_to_m_s*us%s_to_T)
2174 call uvchksum(trim(mesg)//
" BT_force_[uv]", bt_force_u, bt_force_v, cs%debug_BT_HI, haloshift=iev-ie, &
2175 scale=us%L_T_to_m_s*us%s_to_T)
2176 call uvchksum(trim(mesg)//
" BT_rem_[uv]", bt_rem_u, bt_rem_v, cs%debug_BT_HI, haloshift=iev-ie)
2177 call uvchksum(trim(mesg)//
" [uv]bt", ubt, vbt, cs%debug_BT_HI, haloshift=iev-ie, &
2178 scale=us%L_T_to_m_s)
2179 call uvchksum(trim(mesg)//
" [uv]bt_trans", ubt_trans, vbt_trans, cs%debug_BT_HI, haloshift=iev-ie, &
2180 scale=us%L_T_to_m_s)
2181 call uvchksum(trim(mesg)//
" [uv]hbt", uhbt, vhbt, cs%debug_BT_HI, haloshift=iev-ie, &
2182 scale=us%s_to_T*us%L_to_m**2*gv%H_to_m)
2183 if (integral_bt_cont) &
2184 call uvchksum(trim(mesg)//
" [uv]hbt_int", uhbt_int, vhbt_int, cs%debug_BT_HI, haloshift=iev-ie, &
2185 scale=us%L_to_m**2*gv%H_to_m)
2190 do j=js,je ;
do i=is-1,ie
2191 pfu_bt_sum(i,j) = pfu_bt_sum(i,j) + wt_accel2(n) * pfu(i,j)
2195 do j=js-1,je ;
do i=is,ie
2196 pfv_bt_sum(i,j) = pfv_bt_sum(i,j) + wt_accel2(n) * pfv(i,j)
2202 do j=js,je ;
do i=is-1,ie
2203 coru_bt_sum(i,j) = coru_bt_sum(i,j) + wt_accel2(n) * cor_u(i,j)
2207 do j=js-1,je ;
do i=is,ie
2208 corv_bt_sum(i,j) = corv_bt_sum(i,j) + wt_accel2(n) * cor_v(i,j)
2214 do j=js,je ;
do i=is-1,ie
2215 ubt_sum(i,j) = ubt_sum(i,j) + wt_trans(n) * ubt_trans(i,j)
2216 uhbt_sum(i,j) = uhbt_sum(i,j) + wt_trans(n) * uhbt(i,j)
2217 ubt_wtd(i,j) = ubt_wtd(i,j) + wt_vel(n) * ubt(i,j)
2221 do j=js-1,je ;
do i=is,ie
2222 vbt_sum(i,j) = vbt_sum(i,j) + wt_trans(n) * vbt_trans(i,j)
2223 vhbt_sum(i,j) = vhbt_sum(i,j) + wt_trans(n) * vhbt(i,j)
2224 vbt_wtd(i,j) = vbt_wtd(i,j) + wt_vel(n) * vbt(i,j)
2228 if (apply_obcs)
then
2231 call apply_velocity_obcs(obc, ubt, vbt, uhbt, vhbt, &
2232 ubt_trans, vbt_trans, eta, ubt_old, vbt_old, cs%BT_OBC, &
2233 g, ms, us, iev-ie, dtbt, bebt, use_bt_cont, integral_bt_cont, &
2234 n*dtbt, datu, datv, btcl_u, btcl_v, uhbt0, vhbt0, &
2235 ubt_int_prev, vbt_int_prev, uhbt_int_prev, vhbt_int_prev)
2238 if (cs%BT_OBC%apply_u_OBCs)
then
2240 do j=js,je ;
do i=is-1,ie
2241 if (obc%segnum_u(i,j) /= obc_none)
then
2243 ubt_sum(i,j) = ubt_sum_prev(i,j) + wt_trans(n) * ubt_trans(i,j)
2244 uhbt_sum(i,j) = uhbt_sum_prev(i,j) + wt_trans(n) * uhbt(i,j)
2245 ubt_wtd(i,j) = ubt_wtd_prev(i,j) + wt_vel(n) * ubt(i,j)
2246 if (integral_bt_cont)
then
2247 uhbt_int(i,j) = uhbt_int_prev(i,j) + dtbt * uhbt(i,j)
2248 ubt_int(i,j) = ubt_int_prev(i,j) + dtbt * ubt_trans(i,j)
2253 if (cs%BT_OBC%apply_v_OBCs)
then
2255 do j=js-1,je ;
do i=is,ie
2256 if (obc%segnum_v(i,j) /= obc_none)
then
2258 vbt_sum(i,j) = vbt_sum_prev(i,j) + wt_trans(n) * vbt_trans(i,j)
2259 vhbt_sum(i,j) = vhbt_sum_prev(i,j) + wt_trans(n) * vhbt(i,j)
2260 vbt_wtd(i,j) = vbt_wtd_prev(i,j) + wt_vel(n) * vbt(i,j)
2261 if (integral_bt_cont)
then
2262 vbt_int(i,j) = vbt_int_prev(i,j) + dtbt * vbt_trans(i,j)
2263 vhbt_int(i,j) = vhbt_int_prev(i,j) + dtbt * vhbt(i,j)
2270 if (cs%debug_bt)
then
2271 call uvchksum(
"BT [uv]hbt just after OBC", uhbt, vhbt, cs%debug_BT_HI, haloshift=iev-ie, &
2272 scale=us%s_to_T*us%L_to_m**2*gv%H_to_m)
2273 if (integral_bt_cont) &
2274 call uvchksum(
"BT [uv]hbt_int just after OBC", uhbt_int, vhbt_int, cs%debug_BT_HI, &
2275 haloshift=iev-ie, scale=us%L_to_m**2*gv%H_to_m)
2278 if (integral_bt_cont)
then
2280 do j=jsv,jev ;
do i=isv,iev
2281 eta(i,j) = (eta_ic(i,j) + n*eta_src(i,j)) + cs%IareaT(i,j) * &
2282 ((uhbt_int(i-1,j) - uhbt_int(i,j)) + (vhbt_int(i,j-1) - vhbt_int(i,j)))
2283 eta_wtd(i,j) = eta_wtd(i,j) + eta(i,j) * wt_eta(n)
2287 do j=jsv,jev ;
do i=isv,iev
2288 eta(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * cs%IareaT(i,j)) * &
2289 ((uhbt(i-1,j) - uhbt(i,j)) + (vhbt(i,j-1) - vhbt(i,j)))
2290 eta_wtd(i,j) = eta_wtd(i,j) + eta(i,j) * wt_eta(n)
2295 if (do_hifreq_output)
then
2296 time_step_end = time_bt_start + real_to_time(n*us%T_to_s*dtbt)
2297 call enable_averaging(us%T_to_s*dtbt, time_step_end, cs%diag)
2298 if (cs%id_ubt_hifreq > 0)
call post_data(cs%id_ubt_hifreq, ubt(isdb:iedb,jsd:jed), cs%diag)
2299 if (cs%id_vbt_hifreq > 0)
call post_data(cs%id_vbt_hifreq, vbt(isd:ied,jsdb:jedb), cs%diag)
2300 if (cs%id_eta_hifreq > 0)
call post_data(cs%id_eta_hifreq, eta(isd:ied,jsd:jed), cs%diag)
2301 if (cs%id_uhbt_hifreq > 0)
call post_data(cs%id_uhbt_hifreq, uhbt(isdb:iedb,jsd:jed), cs%diag)
2302 if (cs%id_vhbt_hifreq > 0)
call post_data(cs%id_vhbt_hifreq, vhbt(isd:ied,jsdb:jedb), cs%diag)
2303 if (cs%id_eta_pred_hifreq > 0)
call post_data(cs%id_eta_pred_hifreq, eta_pf_bt(isd:ied,jsd:jed), cs%diag)
2306 if (cs%debug_bt)
then
2307 write(mesg,
'("BT step ",I4)') n
2308 call uvchksum(trim(mesg)//
" [uv]bt", ubt, vbt, cs%debug_BT_HI, haloshift=iev-ie, &
2309 scale=us%L_T_to_m_s)
2310 call hchksum(eta, trim(mesg)//
" eta", cs%debug_BT_HI, haloshift=iev-ie, scale=gv%H_to_m)
2313 if (gv%Boussinesq)
then
2314 do j=js,je ;
do i=is,ie
2315 if (eta(i,j) < -gv%Z_to_H*g%bathyT(i,j)) &
2316 call mom_error(warning,
"btstep: eta has dropped below bathyT.")
2319 do j=js,je ;
do i=is,ie
2320 if (eta(i,j) < 0.0) &
2321 call mom_error(warning,
"btstep: negative eta in a non-Boussinesq barotropic solver.")
2326 if (id_clock_calc > 0)
call cpu_clock_end(id_clock_calc)
2327 if (id_clock_calc_post > 0)
call cpu_clock_begin(id_clock_calc_post)
2330 if (do_hifreq_output)
call enable_averaging(time_int_in, time_end_in, cs%diag)
2332 if (cs%answers_2018)
then
2333 i_sum_wt_vel = 1.0 / sum_wt_vel ; i_sum_wt_eta = 1.0 / sum_wt_eta
2334 i_sum_wt_accel = 1.0 / sum_wt_accel ; i_sum_wt_trans = 1.0 / sum_wt_trans
2336 i_sum_wt_vel = 1.0 ; i_sum_wt_eta = 1.0 ; i_sum_wt_accel = 1.0 ; i_sum_wt_trans = 1.0
2339 if (find_etaav)
then ;
do j=js,je ;
do i=is,ie
2340 etaav(i,j) = eta_sum(i,j) * i_sum_wt_accel
2341 enddo ;
enddo ;
endif
2342 do j=js-1,je+1 ;
do i=is-1,ie+1 ; e_anom(i,j) = 0.0 ;
enddo ;
enddo
2343 if (interp_eta_pf)
then
2344 do j=js,je ;
do i=is,ie
2345 e_anom(i,j) = dgeo_de * (0.5 * (eta(i,j) + eta_in(i,j)) - &
2346 (eta_pf_1(i,j) + 0.5*d_eta_pf(i,j)))
2349 do j=js,je ;
do i=is,ie
2350 e_anom(i,j) = dgeo_de * (0.5 * (eta(i,j) + eta_in(i,j)) - eta_pf(i,j))
2353 if (apply_obcs)
then
2355 if (cs%BT_OBC%apply_u_OBCs)
then
2357 do j=js,je ;
do i=is-1,ie
2358 l_seg = obc%segnum_u(i,j)
2359 if (l_seg == obc_none) cycle
2361 if (obc%segment(l_seg)%direction == obc_direction_e)
then
2362 e_anom(i+1,j) = e_anom(i,j)
2363 elseif (obc%segment(l_seg)%direction == obc_direction_w)
then
2364 e_anom(i,j) = e_anom(i+1,j)
2369 if (cs%BT_OBC%apply_v_OBCs)
then
2371 do j=js-1,je ;
do i=is,ie
2372 l_seg = obc%segnum_v(i,j)
2373 if (l_seg == obc_none) cycle
2375 if (obc%segment(l_seg)%direction == obc_direction_n)
then
2376 e_anom(i,j+1) = e_anom(i,j)
2377 elseif (obc%segment(l_seg)%direction == obc_direction_s)
then
2378 e_anom(i,j) = e_anom(i,j+1)
2385 do j=js,je ;
do i=is,ie
2386 eta_out(i,j) = eta_wtd(i,j) * i_sum_wt_eta
2389 if (id_clock_calc_post > 0)
call cpu_clock_end(id_clock_calc_post)
2390 if (id_clock_pass_post > 0)
call cpu_clock_begin(id_clock_pass_post)
2391 if (g%nonblocking_updates)
then
2392 call start_group_pass(cs%pass_e_anom, g%Domain)
2394 if (find_etaav)
call do_group_pass(cs%pass_etaav, g%Domain)
2395 call do_group_pass(cs%pass_e_anom, g%Domain)
2397 if (id_clock_pass_post > 0)
call cpu_clock_end(id_clock_pass_post)
2398 if (id_clock_calc_post > 0)
call cpu_clock_begin(id_clock_calc_post)
2400 if (cs%answers_2018)
then
2401 do j=js,je ;
do i=is-1,ie
2402 cs%ubtav(i,j) = ubt_sum(i,j) * i_sum_wt_trans
2403 uhbtav(i,j) = uhbt_sum(i,j) * i_sum_wt_trans
2404 ubt_wtd(i,j) = ubt_wtd(i,j) * i_sum_wt_vel
2407 do j=js-1,je ;
do i=is,ie
2408 cs%vbtav(i,j) = vbt_sum(i,j) * i_sum_wt_trans
2409 vhbtav(i,j) = vhbt_sum(i,j) * i_sum_wt_trans
2410 vbt_wtd(i,j) = vbt_wtd(i,j) * i_sum_wt_vel
2413 do j=js,je ;
do i=is-1,ie
2414 cs%ubtav(i,j) = ubt_sum(i,j)
2415 uhbtav(i,j) = uhbt_sum(i,j)
2418 do j=js-1,je ;
do i=is,ie
2419 cs%vbtav(i,j) = vbt_sum(i,j)
2420 vhbtav(i,j) = vhbt_sum(i,j)
2425 if (id_clock_calc_post > 0)
call cpu_clock_end(id_clock_calc_post)
2426 if (id_clock_pass_post > 0)
call cpu_clock_begin(id_clock_pass_post)
2427 if (g%nonblocking_updates)
then
2428 call complete_group_pass(cs%pass_e_anom, g%Domain)
2429 if (find_etaav)
call start_group_pass(cs%pass_etaav, g%Domain)
2430 call start_group_pass(cs%pass_ubta_uhbta, g%DoMain)
2432 call do_group_pass(cs%pass_ubta_uhbta, g%Domain)
2434 if (id_clock_pass_post > 0)
call cpu_clock_end(id_clock_pass_post)
2435 if (id_clock_calc_post > 0)
call cpu_clock_begin(id_clock_calc_post)
2440 do j=js,je ;
do i=is-1,ie
2441 accel_layer_u(i,j,k) = (u_accel_bt(i,j) - &
2442 ((pbce(i+1,j,k) - gtot_w(i+1,j)) * e_anom(i+1,j) - &
2443 (pbce(i,j,k) - gtot_e(i,j)) * e_anom(i,j)) * cs%IdxCu(i,j) )
2444 if (abs(accel_layer_u(i,j,k)) < accel_underflow) accel_layer_u(i,j,k) = 0.0
2446 do j=js-1,je ;
do i=is,ie
2447 accel_layer_v(i,j,k) = (v_accel_bt(i,j) - &
2448 ((pbce(i,j+1,k) - gtot_s(i,j+1)) * e_anom(i,j+1) - &
2449 (pbce(i,j,k) - gtot_n(i,j)) * e_anom(i,j)) * cs%IdyCv(i,j) )
2450 if (abs(accel_layer_v(i,j,k)) < accel_underflow) accel_layer_v(i,j,k) = 0.0
2454 if (apply_obcs)
then
2457 if (cs%BT_OBC%apply_u_OBCs)
then ;
do j=js,je ;
do i=is-1,ie
2458 if (obc%segnum_u(i,j) /= obc_none)
then
2459 u_accel_bt(i,j) = (ubt_wtd(i,j) - ubt_first(i,j)) / dt
2460 do k=1,nz ; accel_layer_u(i,j,k) = u_accel_bt(i,j) ;
enddo
2462 enddo ;
enddo ;
endif
2463 if (cs%BT_OBC%apply_v_OBCs)
then ;
do j=js-1,je ;
do i=is,ie
2464 if (obc%segnum_v(i,j) /= obc_none)
then
2465 v_accel_bt(i,j) = (vbt_wtd(i,j) - vbt_first(i,j)) / dt
2466 do k=1,nz ; accel_layer_v(i,j,k) = v_accel_bt(i,j) ;
enddo
2468 enddo ;
enddo ;
endif
2471 if (id_clock_calc_post > 0)
call cpu_clock_end(id_clock_calc_post)
2474 if (query_averaging_enabled(cs%diag))
then
2476 if (cs%gradual_BT_ICs)
then
2477 do j=js,je ;
do i=is-1,ie ; cs%ubt_IC(i,j) = ubt_wtd(i,j) ;
enddo ;
enddo
2478 do j=js-1,je ;
do i=is,ie ; cs%vbt_IC(i,j) = vbt_wtd(i,j) ;
enddo ;
enddo
2482 if (cs%id_PFu_bt > 0)
then
2483 do j=js,je ;
do i=is-1,ie
2484 pfu_bt_sum(i,j) = pfu_bt_sum(i,j) * i_sum_wt_accel
2486 call post_data(cs%id_PFu_bt, pfu_bt_sum(isdb:iedb,jsd:jed), cs%diag)
2488 if (cs%id_PFv_bt > 0)
then
2489 do j=js-1,je ;
do i=is,ie
2490 pfv_bt_sum(i,j) = pfv_bt_sum(i,j) * i_sum_wt_accel
2492 call post_data(cs%id_PFv_bt, pfv_bt_sum(isd:ied,jsdb:jedb), cs%diag)
2494 if (cs%id_Coru_bt > 0)
then
2495 do j=js,je ;
do i=is-1,ie
2496 coru_bt_sum(i,j) = coru_bt_sum(i,j) * i_sum_wt_accel
2498 call post_data(cs%id_Coru_bt, coru_bt_sum(isdb:iedb,jsd:jed), cs%diag)
2500 if (cs%id_Corv_bt > 0)
then
2501 do j=js-1,je ;
do i=is,ie
2502 corv_bt_sum(i,j) = corv_bt_sum(i,j) * i_sum_wt_accel
2504 call post_data(cs%id_Corv_bt, corv_bt_sum(isd:ied,jsdb:jedb), cs%diag)
2506 if (cs%id_ubtdt > 0)
then
2507 do j=js,je ;
do i=is-1,ie
2508 ubt_dt(i,j) = (ubt_wtd(i,j) - ubt_st(i,j))*idt
2510 call post_data(cs%id_ubtdt, ubt_dt(isdb:iedb,jsd:jed), cs%diag)
2512 if (cs%id_vbtdt > 0)
then
2513 do j=js-1,je ;
do i=is,ie
2514 vbt_dt(i,j) = (vbt_wtd(i,j) - vbt_st(i,j))*idt
2516 call post_data(cs%id_vbtdt, vbt_dt(isd:ied,jsdb:jedb), cs%diag)
2519 if (cs%id_ubtforce > 0)
call post_data(cs%id_ubtforce, bt_force_u(isdb:iedb,jsd:jed), cs%diag)
2520 if (cs%id_vbtforce > 0)
call post_data(cs%id_vbtforce, bt_force_v(isd:ied,jsdb:jedb), cs%diag)
2521 if (cs%id_uaccel > 0)
call post_data(cs%id_uaccel, u_accel_bt(isdb:iedb,jsd:jed), cs%diag)
2522 if (cs%id_vaccel > 0)
call post_data(cs%id_vaccel, v_accel_bt(isd:ied,jsdb:jedb), cs%diag)
2524 if (cs%id_eta_cor > 0)
call post_data(cs%id_eta_cor, cs%eta_cor, cs%diag)
2525 if (cs%id_eta_bt > 0)
call post_data(cs%id_eta_bt, eta_out, cs%diag)
2526 if (cs%id_gtotn > 0)
call post_data(cs%id_gtotn, gtot_n(isd:ied,jsd:jed), cs%diag)
2527 if (cs%id_gtots > 0)
call post_data(cs%id_gtots, gtot_s(isd:ied,jsd:jed), cs%diag)
2528 if (cs%id_gtote > 0)
call post_data(cs%id_gtote, gtot_e(isd:ied,jsd:jed), cs%diag)
2529 if (cs%id_gtotw > 0)
call post_data(cs%id_gtotw, gtot_w(isd:ied,jsd:jed), cs%diag)
2530 if (cs%id_ubt > 0)
call post_data(cs%id_ubt, ubt_wtd(isdb:iedb,jsd:jed), cs%diag)
2531 if (cs%id_vbt > 0)
call post_data(cs%id_vbt, vbt_wtd(isd:ied,jsdb:jedb), cs%diag)
2532 if (cs%id_ubtav > 0)
call post_data(cs%id_ubtav, cs%ubtav, cs%diag)
2533 if (cs%id_vbtav > 0)
call post_data(cs%id_vbtav, cs%vbtav, cs%diag)
2534 if (cs%id_visc_rem_u > 0)
call post_data(cs%id_visc_rem_u, visc_rem_u, cs%diag)
2535 if (cs%id_visc_rem_v > 0)
call post_data(cs%id_visc_rem_v, visc_rem_v, cs%diag)
2537 if (cs%id_frhatu > 0)
call post_data(cs%id_frhatu, cs%frhatu, cs%diag)
2538 if (cs%id_uhbt > 0)
call post_data(cs%id_uhbt, uhbtav, cs%diag)
2539 if (cs%id_frhatv > 0)
call post_data(cs%id_frhatv, cs%frhatv, cs%diag)
2540 if (cs%id_vhbt > 0)
call post_data(cs%id_vhbt, vhbtav, cs%diag)
2541 if (cs%id_uhbt0 > 0)
call post_data(cs%id_uhbt0, uhbt0(isdb:iedb,jsd:jed), cs%diag)
2542 if (cs%id_vhbt0 > 0)
call post_data(cs%id_vhbt0, vhbt0(isd:ied,jsdb:jedb), cs%diag)
2544 if (cs%id_frhatu1 > 0)
call post_data(cs%id_frhatu1, cs%frhatu1, cs%diag)
2545 if (cs%id_frhatv1 > 0)
call post_data(cs%id_frhatv1, cs%frhatv1, cs%diag)
2547 if (use_bt_cont)
then
2548 if (cs%id_BTC_FA_u_EE > 0)
call post_data(cs%id_BTC_FA_u_EE, bt_cont%FA_u_EE, cs%diag)
2549 if (cs%id_BTC_FA_u_E0 > 0)
call post_data(cs%id_BTC_FA_u_E0, bt_cont%FA_u_E0, cs%diag)
2550 if (cs%id_BTC_FA_u_W0 > 0)
call post_data(cs%id_BTC_FA_u_W0, bt_cont%FA_u_W0, cs%diag)
2551 if (cs%id_BTC_FA_u_WW > 0)
call post_data(cs%id_BTC_FA_u_WW, bt_cont%FA_u_WW, cs%diag)
2552 if (cs%id_BTC_uBT_EE > 0)
call post_data(cs%id_BTC_uBT_EE, bt_cont%uBT_EE, cs%diag)
2553 if (cs%id_BTC_uBT_WW > 0)
call post_data(cs%id_BTC_uBT_WW, bt_cont%uBT_WW, cs%diag)
2554 if (cs%id_BTC_FA_u_rat0 > 0)
then
2556 do j=js,je ;
do i=is-1,ie
2557 if ((g%mask2dCu(i,j) > 0.0) .and. (bt_cont%FA_u_W0(i,j) > 0.0))
then
2558 tmp_u(i,j) = (bt_cont%FA_u_E0(i,j)/ bt_cont%FA_u_W0(i,j))
2563 call post_data(cs%id_BTC_FA_u_rat0, tmp_u, cs%diag)
2565 if (cs%id_BTC_FA_v_NN > 0)
call post_data(cs%id_BTC_FA_v_NN, bt_cont%FA_v_NN, cs%diag)
2566 if (cs%id_BTC_FA_v_N0 > 0)
call post_data(cs%id_BTC_FA_v_N0, bt_cont%FA_v_N0, cs%diag)
2567 if (cs%id_BTC_FA_v_S0 > 0)
call post_data(cs%id_BTC_FA_v_S0, bt_cont%FA_v_S0, cs%diag)
2568 if (cs%id_BTC_FA_v_SS > 0)
call post_data(cs%id_BTC_FA_v_SS, bt_cont%FA_v_SS, cs%diag)
2569 if (cs%id_BTC_vBT_NN > 0)
call post_data(cs%id_BTC_vBT_NN, bt_cont%vBT_NN, cs%diag)
2570 if (cs%id_BTC_vBT_SS > 0)
call post_data(cs%id_BTC_vBT_SS, bt_cont%vBT_SS, cs%diag)
2571 if (cs%id_BTC_FA_v_rat0 > 0)
then
2573 do j=js-1,je ;
do i=is,ie
2574 if ((g%mask2dCv(i,j) > 0.0) .and. (bt_cont%FA_v_S0(i,j) > 0.0))
then
2575 tmp_v(i,j) = (bt_cont%FA_v_N0(i,j)/ bt_cont%FA_v_S0(i,j))
2580 call post_data(cs%id_BTC_FA_v_rat0, tmp_v, cs%diag)
2582 if (cs%id_BTC_FA_h_rat0 > 0)
then
2584 do j=js,je ;
do i=is,ie
2586 if ((g%mask2dCu(i,j) > 0.0) .and. (bt_cont%FA_u_W0(i,j) > 0.0) .and. (bt_cont%FA_u_E0(i,j) > 0.0))
then
2587 if (bt_cont%FA_u_W0(i,j) > bt_cont%FA_u_E0(i,j))
then
2588 tmp_h(i,j) = max(tmp_h(i,j), (bt_cont%FA_u_W0(i,j)/ bt_cont%FA_u_E0(i,j)))
2590 tmp_h(i,j) = max(tmp_h(i,j), (bt_cont%FA_u_E0(i,j)/ bt_cont%FA_u_W0(i,j)))
2593 if ((g%mask2dCu(i-1,j) > 0.0) .and. (bt_cont%FA_u_W0(i-1,j) > 0.0) .and. (bt_cont%FA_u_E0(i-1,j) > 0.0))
then
2594 if (bt_cont%FA_u_W0(i-1,j) > bt_cont%FA_u_E0(i-1,j))
then
2595 tmp_h(i,j) = max(tmp_h(i,j), (bt_cont%FA_u_W0(i-1,j)/ bt_cont%FA_u_E0(i-1,j)))
2597 tmp_h(i,j) = max(tmp_h(i,j), (bt_cont%FA_u_E0(i-1,j)/ bt_cont%FA_u_W0(i-1,j)))
2600 if ((g%mask2dCv(i,j) > 0.0) .and. (bt_cont%FA_v_S0(i,j) > 0.0) .and. (bt_cont%FA_v_N0(i,j) > 0.0))
then
2601 if (bt_cont%FA_v_S0(i,j) > bt_cont%FA_v_N0(i,j))
then
2602 tmp_h(i,j) = max(tmp_h(i,j), (bt_cont%FA_v_S0(i,j)/ bt_cont%FA_v_N0(i,j)))
2604 tmp_h(i,j) = max(tmp_h(i,j), (bt_cont%FA_v_N0(i,j)/ bt_cont%FA_v_S0(i,j)))
2607 if ((g%mask2dCv(i,j-1) > 0.0) .and. (bt_cont%FA_v_S0(i,j-1) > 0.0) .and. (bt_cont%FA_v_N0(i,j-1) > 0.0))
then
2608 if (bt_cont%FA_v_S0(i,j-1) > bt_cont%FA_v_N0(i,j-1))
then
2609 tmp_h(i,j) = max(tmp_h(i,j), (bt_cont%FA_v_S0(i,j-1)/ bt_cont%FA_v_N0(i,j-1)))
2611 tmp_h(i,j) = max(tmp_h(i,j), (bt_cont%FA_v_N0(i,j-1)/ bt_cont%FA_v_S0(i,j-1)))
2615 call post_data(cs%id_BTC_FA_h_rat0, tmp_h, cs%diag)
2619 if (cs%id_frhatu1 > 0) cs%frhatu1(:,:,:) = cs%frhatu(:,:,:)
2620 if (cs%id_frhatv1 > 0) cs%frhatv1(:,:,:) = cs%frhatv(:,:,:)
2623 if ((
present(adp)) .and. (
associated(adp%diag_hfrac_u)))
then
2624 do k=1,nz ;
do j=js,je ;
do i=is-1,ie
2625 adp%diag_hfrac_u(i,j,k) = cs%frhatu(i,j,k)
2626 enddo ;
enddo ;
enddo
2628 if ((
present(adp)) .and. (
associated(adp%diag_hfrac_v)))
then
2629 do k=1,nz ;
do j=js-1,je ;
do i=is,ie
2630 adp%diag_hfrac_v(i,j,k) = cs%frhatv(i,j,k)
2631 enddo ;
enddo ;
enddo
2634 if (g%nonblocking_updates)
then
2635 if (find_etaav)
call complete_group_pass(cs%pass_etaav, g%Domain)
2636 call complete_group_pass(cs%pass_ubta_uhbta, g%Domain)