This subroutine does 1-d flux-form advection using a monotonic piecewise linear scheme.
657 type(ocean_grid_type),
intent(inout) :: G
658 type(verticalGrid_type),
intent(in) :: GV
659 type(tracer_type),
dimension(ntr),
intent(inout) :: Tr
660 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: hprev
662 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: vhr
664 real,
dimension(SZI_(G),SZJB_(G)),
intent(inout) :: vh_neglect
666 type(ocean_OBC_type),
pointer :: OBC
667 logical,
dimension(SZJB_(G),SZK_(G)),
intent(inout) :: domore_v
669 real,
intent(in) :: Idt
670 integer,
intent(in) :: ntr
671 integer,
intent(in) :: is
672 integer,
intent(in) :: ie
673 integer,
intent(in) :: js
674 integer,
intent(in) :: je
675 integer,
intent(in) :: k
676 logical,
intent(in) :: usePPM
677 logical,
intent(in) :: useHuynh
680 real,
dimension(SZI_(G),ntr,SZJ_(G)) :: &
682 real,
dimension(SZI_(G),ntr,SZJB_(G)) :: &
686 real :: vhh(SZI_(G),SZJB_(G))
692 real,
dimension(SZIB_(G)) :: &
700 logical :: do_j_tr(SZJ_(G))
701 logical :: do_i(SZIB_(G))
703 integer :: i, j, j2, m, n, j_up, stencil
704 real :: aR, aL, dMx, dMn, Tp, Tc, Tm, dA, mA, a6
705 real :: fac1,v_L_in,v_L_out
706 integer :: jshift, jdir
708 type(OBC_segment_type),
pointer :: segment=>null()
709 logical :: usePLMslope
711 useplmslope = .not. (useppm .and. usehuynh)
714 if (useppm .and. .not. usehuynh) stencil = 2
716 min_h = 0.1*gv%Angstrom_H
717 h_neglect = gv%H_subroundoff
732 do j=js-1,je ;
if (domore_v(j,k))
then ;
do j2=1-stencil,stencil ; do_j_tr(j+j2) = .true. ;
enddo ;
endif ;
enddo
736 if (useplmslope)
then
737 do j=js-stencil,je+stencil ;
if (do_j_tr(j))
then ;
do m=1,ntr ;
do i=is,ie
752 tp = tr(m)%t(i,j+1,k) ; tc = tr(m)%t(i,j,k) ; tm = tr(m)%t(i,j-1,k)
753 dmx = max( tp, tc, tm ) - tc
754 dmn= tc - min( tp, tc, tm )
755 slope_y(i,m,j) = g%mask2dCv(i,j)*g%mask2dCv(i,j-1) * &
756 sign( min(0.5*abs(tp-tm), 2.0*dmx, 2.0*dmn), tp-tm )
757 enddo ;
enddo ;
endif ;
enddo
764 do j=js-1,je ;
if (domore_v(j,k))
then
765 domore_v(j,k) = .false.
768 if (vhr(i,j,k) == 0.0)
then
771 elseif (vhr(i,j,k) < 0.0)
then
772 hup = hprev(i,j+1,k) - g%areaT(i,j+1)*min_h
773 hlos = max(0.0,vhr(i,j+1,k))
774 if ((((hup - hlos) + vhr(i,j,k)) < 0.0) .and. &
775 ((0.5*hup + vhr(i,j,k)) < 0.0))
then
776 vhh(i,j) = min(-0.5*hup,-hup+hlos,0.0)
777 domore_v(j,k) = .true.
779 vhh(i,j) = vhr(i,j,k)
782 cfl(i) = - vhh(i,j) / (hprev(i,j+1,k)+h_neglect)
784 hup = hprev(i,j,k) - g%areaT(i,j)*min_h
785 hlos = max(0.0,-vhr(i,j-1,k))
786 if ((((hup - hlos) - vhr(i,j,k)) < 0.0) .and. &
787 ((0.5*hup - vhr(i,j,k)) < 0.0))
then
788 vhh(i,j) = max(0.5*hup,hup-hlos,0.0)
789 domore_v(j,k) = .true.
791 vhh(i,j) = vhr(i,j,k)
794 cfl(i) = vhh(i,j) / (hprev(i,j,k)+h_neglect)
799 do m=1,ntr ;
do i=is,ie
801 if (vhh(i,j) >= 0.0)
then
808 tp = tr(m)%t(i,j_up+1,k) ; tc = tr(m)%t(i,j_up,k) ; tm = tr(m)%t(i,j_up-1,k)
811 al = ( 5.*tc + ( 2.*tm - tp ) )/6.
812 al = max( min(tc,tm), al) ; al = min( max(tc,tm), al)
813 ar = ( 5.*tc + ( 2.*tp - tm ) )/6.
814 ar = max( min(tc,tp), ar) ; ar = min( max(tc,tp), ar)
816 al = 0.5 * ((tm + tc) + (slope_y(i,m,j_up-1) - slope_y(i,m,j_up)) / 3.)
817 ar = 0.5 * ((tc + tp) + (slope_y(i,m,j_up) - slope_y(i,m,j_up+1)) / 3.)
820 da = ar - al ; ma = 0.5*( ar + al )
821 if (g%mask2dCv(i,j_up)*g%mask2dCv(i,j_up-1)*(tp-tc)*(tc-tm) <= 0.)
then
823 elseif ( da*(tc-ma) > (da*da)/6. )
then
825 elseif ( da*(tc-ma) < - (da*da)/6. )
then
829 a6 = 6.*tc - 3. * (ar + al)
831 if (vhh(i,j) >= 0.0)
then
832 flux_y(i,m,j) = vhh(i,j)*( ar - 0.5 * cfl(i) * ( &
833 ( ar - al ) - a6 * ( 1. - 2./3. * cfl(i) ) ) )
835 flux_y(i,m,j) = vhh(i,j)*( al + 0.5 * cfl(i) * ( &
836 ( ar - al ) + a6 * ( 1. - 2./3. * cfl(i) ) ) )
840 do m=1,ntr ;
do i=is,ie
841 if (vhh(i,j) >= 0.0)
then
851 flux_y(i,m,j) = vhh(i,j)*( tc + 0.5 * slope_y(i,m,j) * ( 1. - cfl(i) ) )
863 tc = tr(m)%t(i,j+1,k)
864 flux_y(i,m,j) = vhh(i,j)*( tc - 0.5 * slope_y(i,m,j+1) * ( 1. - cfl(i) ) )
871 if (
associated(obc))
then ;
if (obc%OBC_pe)
then
872 if (obc%specified_v_BCs_exist_globally)
then
873 do n=1,obc%number_of_segments
874 segment=>obc%segment(n)
875 if (.not. segment%specified) cycle
876 if (.not.
associated(segment%tr_Reg)) cycle
877 if (obc%segment(n)%is_N_or_S)
then
878 if (j >= segment%HI%JsdB .and. j<= segment%HI%JedB)
then
879 do i=segment%HI%isd,segment%HI%ied
882 if ((vhr(i,j,k) > 0.0) .and. (segment%direction == obc_direction_s) .or. &
883 (vhr(i,j,k) < 0.0) .and. (segment%direction == obc_direction_n))
then
884 vhh(i,j) = vhr(i,j,k)
886 if (
associated(segment%tr_Reg%Tr(m)%t))
then
887 flux_y(i,m,j) = vhh(i,j)*obc%segment(n)%tr_Reg%Tr(m)%tres(i,j,k)
888 else ; flux_y(i,m,j) = vhh(i,j)*obc%segment(n)%tr_Reg%Tr(m)%OBC_inflow_conc ;
endif
898 if (obc%open_v_BCs_exist_globally)
then
899 do n=1,obc%number_of_segments
900 segment=>obc%segment(n)
901 if (segment%specified) cycle
902 if (.not.
associated(segment%tr_Reg)) cycle
903 if (segment%is_N_or_S .and. &
904 (j >= segment%HI%JsdB .and. j<= segment%HI%JedB))
then
907 if (segment%direction == obc_direction_s)
then
911 do i=segment%HI%isd,segment%HI%ied
915 if (
associated(segment%tr_Reg%Tr(m)%tres))
then
917 v_l_in=max(jdir*vhh(i,j)*segment%Tr_InvLscale3_in,0.)
918 v_l_out=min(jdir*vhh(i,j)*segment%Tr_InvLscale3_out,0.)
919 fac1=1.0+dt*(v_l_in-v_l_out)
920 segment%tr_Reg%Tr(m)%tres(i,j,k)= (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(i,j,k) + &
921 dt*v_l_in*tr(m)%t(i,j+jshift,k) - &
922 dt*v_l_out*segment%tr_Reg%Tr(m)%t(i,j,k))
926 if ((vhr(i,j,k) > 0.0) .and. (g%mask2dT(i,j) < 0.5) .or. &
927 (vhr(i,j,k) < 0.0) .and. (g%mask2dT(i,j+1) < 0.5))
then
928 vhh(i,j) = vhr(i,j,k)
930 if (
associated(segment%tr_Reg%Tr(m)%t))
then
931 flux_y(i,m,j) = vhh(i,j)*segment%tr_Reg%Tr(m)%tres(i,j,k)
932 else ; flux_y(i,m,j) = vhh(i,j)*segment%tr_Reg%Tr(m)%OBC_inflow_conc ;
endif
942 do i=is,ie ; vhh(i,j) = 0.0 ;
enddo
943 do m=1,ntr ;
do i=is,ie ; flux_y(i,m,j) = 0.0 ;
enddo ;
enddo
946 do j=js-1,je ;
do i=is,ie
947 vhr(i,j,k) = vhr(i,j,k) - vhh(i,j)
948 if (abs(vhr(i,j,k)) < vh_neglect(i,j)) vhr(i,j,k) = 0.0
953 do j=js,je ;
if (do_j_tr(j))
then
955 if ((vhh(i,j) /= 0.0) .or. (vhh(i,j-1) /= 0.0))
then
957 hlst(i) = hprev(i,j,k)
958 hprev(i,j,k) = max(hprev(i,j,k) - (vhh(i,j) - vhh(i,j-1)), 0.0)
959 if (hprev(i,j,k) <= 0.0)
then ; do_i(i) = .false.
960 elseif (hprev(i,j,k) < h_neglect*g%areaT(i,j))
then
961 hlst(i) = hlst(i) + (h_neglect*g%areaT(i,j) - hprev(i,j,k))
962 ihnew(i) = 1.0 / (h_neglect*g%areaT(i,j))
963 else ; ihnew(i) = 1.0 / hprev(i,j,k) ;
endif
964 else ; do_i(i) = .false. ;
endif
969 do i=is,ie ;
if (do_i(i))
then
970 tr(m)%t(i,j,k) = (tr(m)%t(i,j,k) * hlst(i) - &
971 (flux_y(i,m,j) - flux_y(i,m,j-1))) * ihnew(i)
975 if (
associated(tr(m)%ad_y))
then ;
do i=is,ie ;
if (do_i(i))
then
976 tr(m)%ad_y(i,j,k) = tr(m)%ad_y(i,j,k) + flux_y(i,m,j)*idt
977 endif ;
enddo ;
endif
978 if (
associated(tr(m)%ad2d_y))
then ;
do i=is,ie ;
if (do_i(i))
then
979 tr(m)%ad2d_y(i,j) = tr(m)%ad2d_y(i,j) + flux_y(i,m,j)*idt
980 endif ;
enddo ;
endif
984 if (
associated(tr(m)%advection_xy))
then
985 do i=is,ie ;
if (do_i(i))
then
986 tr(m)%advection_xy(i,j,k) = tr(m)%advection_xy(i,j,k) - (flux_y(i,m,j) - flux_y(i,m,j-1))* idt * g%IareaT(i,j)