diff --git a/src/oce_adv_tra_driver.F90 b/src/oce_adv_tra_driver.F90 index 4f2e9e1b8..522ed705a 100644 --- a/src/oce_adv_tra_driver.F90 +++ b/src/oce_adv_tra_driver.F90 @@ -125,59 +125,58 @@ subroutine do_oce_adv_tra(dt, vel, w, wi, we, tr_num, dynamics, tracers, partit, !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(e, enodes, el, nl1, nu1, nl2, nu2, nu12, nl12, nz) #if !defined(DISABLE_OPENACC_ATOMICS) - !$ACC PARALLEL LOOP GANG PRIVATE(enodes, el) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) + !$ACC PARALLEL LOOP GANG VECTOR COLLAPSE(2) PRIVATE(enodes, el) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) #else !$ACC UPDATE SELF(fct_lo, adv_flux_hor) #endif do e=1, myDim_edge2D - enodes=edges(:,e) - el=edge_tri(:,e) - nl1=nlevels(el(1))-1 - nu1=ulevels(el(1)) - nl2=0 - nu2=0 - if(el(2)>0) then - nl2=nlevels(el(2))-1 - nu2=ulevels(el(2)) - end if + do nz=1, nl + enodes=edges(:,e) + el=edge_tri(:,e) + nl1=nlevels(el(1))-1 + nu1=ulevels(el(1)) + nl2=0 + nu2=0 + if(el(2)>0) then + nl2=nlevels(el(2))-1 + nu2=ulevels(el(2)) + end if - nl12 = max(nl1,nl2) - nu12 = nu1 - if (nu2>0) nu12 = min(nu1,nu2) + nl12 = max(nl1,nl2) + nu12 = nu1 + if (nu2>0) nu12 = min(nu1,nu2) + + if(nu12 <= nz .and. nz <= nl12) then !!PS do nz=1, max(nl1, nl2) #if defined(_OPENMP) && !defined(__openmp_reproducible) - call omp_set_lock(partit%plock(enodes(1))) + call omp_set_lock(partit%plock(enodes(1))) #else !$OMP ORDERED #endif -#if !defined(DISABLE_OPENACC_ATOMICS) - !$ACC LOOP VECTOR -#endif - do nz=nu12, nl12 + ! do nz=nu12, nl12 #if !defined(DISABLE_OPENACC_ATOMICS) !$ACC ATOMIC UPDATE #endif - fct_LO(nz, enodes(1))=fct_LO(nz, enodes(1))+adv_flux_hor(nz, e) + fct_LO(nz, enodes(1))=fct_LO(nz, enodes(1))+adv_flux_hor(nz, e) #if defined(_OPENMP) && !defined(__openmp_reproducible) - end do - call omp_unset_lock(partit%plock(enodes(1))) - call omp_set_lock (partit%plock(enodes(2))) - do nz=nu12, nl12 + ! end do + call omp_unset_lock(partit%plock(enodes(1))) + call omp_set_lock (partit%plock(enodes(2))) + ! do nz=nu12, nl12 #endif #if !defined(DISABLE_OPENACC_ATOMICS) - !$ACC ATOMIC UPDATE -#endif - fct_LO(nz, enodes(2))=fct_LO(nz, enodes(2))-adv_flux_hor(nz, e) - end do -#if !defined(DISABLE_OPENACC_ATOMICS) - !$ACC END LOOP + !$ACC ATOMIC UPDATE #endif + fct_LO(nz, enodes(2))=fct_LO(nz, enodes(2))-adv_flux_hor(nz, e) + ! end do #if defined(_OPENMP) && !defined(__openmp_reproducible) - call omp_unset_lock(partit%plock(enodes(2))) + call omp_unset_lock(partit%plock(enodes(2))) #else !$OMP END ORDERED #endif + end if + end do end do #if !defined(DISABLE_OPENACC_ATOMICS) !$ACC END PARALLEL LOOP @@ -192,16 +191,20 @@ subroutine do_oce_adv_tra(dt, vel, w, wi, we, tr_num, dynamics, tracers, partit, ! update the LO solution for vertical contribution !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(n, nu1, nl1, nz) - !$ACC PARALLEL LOOP GANG PRESENT(fct_LO) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) + !$ACC PARALLEL LOOP GANG VECTOR COLLAPSE(2) PRESENT(fct_LO) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) do n=1, myDim_nod2D - nu1 = ulevels_nod2D(n) - nl1 = nlevels_nod2D(n) - !!PS do nz=1, nlevels_nod2D(n)-1 - !$ACC LOOP VECTOR - do nz= nu1, nl1-1 - fct_LO(nz,n)=(ttf(nz,n)*hnode(nz,n)+(fct_LO(nz,n)+(adv_flux_ver(nz, n)-adv_flux_ver(nz+1, n)))*dt/areasvol(nz,n))/hnode_new(nz,n) + do nz= 1, nl + nu1 = ulevels_nod2D(n) + nl1 = nlevels_nod2D(n) + !!PS do nz=1, nlevels_nod2D(n)-1 + if (nu1 <= nz .and. nz < nl1) then + ! do nz= nu1, nl1-1 + fct_LO(nz,n)=(ttf(nz,n)*hnode(nz,n)+(fct_LO(nz,n) & + +(adv_flux_ver(nz, n)-adv_flux_ver(nz+1, n))) & + *dt/areasvol(nz,n))/hnode_new(nz,n) + ! end do + end if end do - !$ACC END LOOP end do !$ACC END PARALLEL LOOP !$OMP END PARALLEL DO @@ -288,6 +291,7 @@ subroutine oce_tra_adv_flux2dtracer(dt, dttf_h, dttf_v, flux_h, flux_v, partit, real(kind=WP), intent(inout) :: flux_h(mesh%nl-1, partit%myDim_edge2D) real(kind=WP), intent(inout) :: flux_v(mesh%nl, partit%myDim_nod2D) logical, optional :: use_lo + logical :: use_lo_present real(kind=WP), optional :: lo (mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) real(kind=WP), optional :: ttf(mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) integer :: n, nz, k, elem, enodes(3), num, el(2), nu12, nl12, nu1, nu2, nl1, nl2, edge @@ -300,93 +304,100 @@ subroutine oce_tra_adv_flux2dtracer(dt, dttf_h, dttf_v, flux_h, flux_v, partit, ! Vertical !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(n, nz, k, elem, enodes, num, el, nu12, nl12, nu1, nu2, nl1, nl2, edge) - if (present(use_lo)) then - if (use_lo) then -!$OMP DO - !$ACC PARALLEL LOOP GANG DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) - do n=1, myDim_nod2d - nu1 = ulevels_nod2D(n) - nl1 = nlevels_nod2D(n) - !!PS do nz=1,nlevels_nod2D(n)-1 - !$ACC LOOP VECTOR - do nz=nu1, nl1-1 - dttf_v(nz,n)=dttf_v(nz,n)-ttf(nz,n)*hnode(nz,n)+LO(nz,n)*hnode_new(nz,n) - end do - !$ACC END LOOP - end do - !$ACC END PARALLEL LOOP -!$OMP END DO - end if - end if -!$OMP DO - !$ACC PARALLEL LOOP GANG DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) +! if (present(use_lo)) then +! if (use_lo) then +! !$OMP DO +! !$ACC PARALLEL LOOP GANG VECTOR COLLAPSE(2) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) +! do n=1, myDim_nod2d +! do nz=1, nl +! nu1 = ulevels_nod2D(n) +! nl1 = nlevels_nod2D(n) +! !!PS do nz=1,nlevels_nod2D(n)-1 +! if (nu1 <= nz .and. nz < nl1) then +! ! do nz=nu1, nl1-1 +! dttf_v(nz,n)=dttf_v(nz,n)-ttf(nz,n)*hnode(nz,n)+LO(nz,n)*hnode_new(nz,n) +! ! end do +! end if +! end do +! end do +! !$ACC END PARALLEL LOOP +! !$OMP END DO +! end if +! end if +! !$OMP DO + use_lo_present = present(use_lo) .and. use_lo + !$ACC PARALLEL LOOP GANG VECTOR COLLAPSE(2) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) do n=1, myDim_nod2d - nu1 = ulevels_nod2D(n) - nl1 = nlevels_nod2D(n) - !$ACC LOOP VECTOR - do nz=nu1,nl1-1 - dttf_v(nz,n)=dttf_v(nz,n) + (flux_v(nz,n)-flux_v(nz+1,n))*dt/areasvol(nz,n) + do nz=1, nl + nu1 = ulevels_nod2D(n) + nl1 = nlevels_nod2D(n) + if (nu1 <= nz .and. nz < nl1) then + ! do nz=nu1,nl1-1 + if (use_lo_present) then + dttf_v(nz,n)=dttf_v(nz,n)-ttf(nz,n)*hnode(nz,n)+LO(nz,n)*hnode_new(nz,n) + end if + dttf_v(nz,n)=dttf_v(nz,n) + (flux_v(nz,n)-flux_v(nz+1,n))*dt/areasvol(nz,n) + ! end do + end if end do - !$ACC END LOOP end do !$ACC END PARALLEL LOOP -!$OMP END DO +!!$OMP END DO ! Horizontal !$OMP DO #if !defined(DISABLE_OPENACC_ATOMICS) - !$ACC PARALLEL LOOP GANG PRIVATE(enodes, el) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) + !$ACC PARALLEL LOOP GANG VECTOR COLLAPSE(2) PRIVATE(enodes, el) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) #else !$ACC UPDATE SELF(dttf_h, flux_h) #endif do edge=1, myDim_edge2D - enodes(1:2)=edges(:,edge) - el=edge_tri(:,edge) - nl1=nlevels(el(1))-1 - nu1=ulevels(el(1)) + do nz=1, nl + enodes(1:2)=edges(:,edge) + el=edge_tri(:,edge) + nl1=nlevels(el(1))-1 + nu1=ulevels(el(1)) - nl2=0 - nu2=0 - if(el(2)>0) then - nl2=nlevels(el(2))-1 - nu2=ulevels(el(2)) - end if + nl2=0 + nu2=0 + if(el(2)>0) then + nl2=nlevels(el(2))-1 + nu2=ulevels(el(2)) + end if - nl12 = max(nl1,nl2) - nu12 = nu1 - if (nu2>0) nu12 = min(nu1,nu2) + nl12 = max(nl1,nl2) + nu12 = nu1 + if (nu2>0) nu12 = min(nu1,nu2) + + if(nu12 <= nz .and. nz <= nl12) then #if defined(_OPENMP) && !defined(__openmp_reproducible) - call omp_set_lock(partit%plock(enodes(1))) + call omp_set_lock(partit%plock(enodes(1))) #else !$OMP ORDERED #endif + ! do nz=nu12, nl12 #if !defined(DISABLE_OPENACC_ATOMICS) - !$ACC LOOP VECTOR -#endif - do nz=nu12, nl12 -#if !defined(DISABLE_OPENACC_ATOMICS) - !$ACC ATOMIC UPDATE + !$ACC ATOMIC UPDATE #endif - dttf_h(nz,enodes(1))=dttf_h(nz,enodes(1))+flux_h(nz,edge)*dt/areasvol(nz,enodes(1)) + dttf_h(nz,enodes(1))=dttf_h(nz,enodes(1))+flux_h(nz,edge)*dt/areasvol(nz,enodes(1)) #if defined(_OPENMP) && !defined(__openmp_reproducible) - end do - call omp_unset_lock(partit%plock(enodes(1))) - call omp_set_lock (partit%plock(enodes(2))) - do nz=nu12, nl12 -#endif -#if !defined(DISABLE_OPENACC_ATOMICS) - !$ACC ATOMIC UPDATE + ! end do + call omp_unset_lock(partit%plock(enodes(1))) + call omp_set_lock (partit%plock(enodes(2))) + ! do nz=nu12, nl12 #endif - dttf_h(nz,enodes(2))=dttf_h(nz,enodes(2))-flux_h(nz,edge)*dt/areasvol(nz,enodes(2)) - end do #if !defined(DISABLE_OPENACC_ATOMICS) - !$ACC END LOOP + !$ACC ATOMIC UPDATE #endif + dttf_h(nz,enodes(2))=dttf_h(nz,enodes(2))-flux_h(nz,edge)*dt/areasvol(nz,enodes(2)) + ! end do #if defined(_OPENMP) && !defined(__openmp_reproducible) - call omp_unset_lock(partit%plock(enodes(2))) + call omp_unset_lock(partit%plock(enodes(2))) #else !$OMP END ORDERED #endif + end if + end do end do #if !defined(DISABLE_OPENACC_ATOMICS) diff --git a/src/oce_adv_tra_fct.F90 b/src/oce_adv_tra_fct.F90 index 24a43bef2..07c2aaa92 100644 --- a/src/oce_adv_tra_fct.F90 +++ b/src/oce_adv_tra_fct.F90 @@ -116,16 +116,18 @@ subroutine oce_tra_adv_fct(dt, ttf, lo, adf_h, adf_v, fct_ttf_min, fct_ttf_max, !$ACC DATA CREATE(tvert_max, tvert_min) !$OMP DO - !$ACC PARALLEL LOOP GANG DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) + !$ACC PARALLEL LOOP GANG VECTOR COLLAPSE(2) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) do n=1,myDim_nod2D + edim_nod2d - nu1 = ulevels_nod2D(n) - nl1 = nlevels_nod2D(n) - !$ACC LOOP VECTOR - do nz=nu1, nl1-1 - fct_ttf_max(nz,n)=max(LO(nz,n), ttf(nz,n)) - fct_ttf_min(nz,n)=min(LO(nz,n), ttf(nz,n)) + do nz=1, nl + nu1 = ulevels_nod2D(n) + nl1 = nlevels_nod2D(n) + if(nu1 <= nz .and. nz < nl1) then + ! do nz=nu1, nl1-1 + fct_ttf_max(nz,n)=max(LO(nz,n), ttf(nz,n)) + fct_ttf_min(nz,n)=min(LO(nz,n), ttf(nz,n)) + ! end do + end if end do - !$ACC END LOOP end do !$ACC END PARALLEL LOOP !$OMP END DO @@ -134,25 +136,26 @@ subroutine oce_tra_adv_fct(dt, ttf, lo, adf_h, adf_v, fct_ttf_min, fct_ttf_max, ! (only layers below the first and above the last layer) ! look for max, min bounds for each element --> AUX here auxilary array !$OMP DO - !$ACC PARALLEL LOOP GANG PRIVATE(enodes) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) + !$ACC PARALLEL LOOP GANG VECTOR COLLAPSE(2) PRIVATE(enodes) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) do elem=1, myDim_elem2D - enodes=elem2D_nodes(:,elem) - nu1 = ulevels(elem) - nl1 = nlevels(elem) - !$ACC LOOP VECTOR - do nz=nu1, nl1-1 - AUX(1,nz,elem)=max(fct_ttf_max(nz,enodes(1)), fct_ttf_max(nz,enodes(2)), fct_ttf_max(nz,enodes(3))) - AUX(2,nz,elem)=min(fct_ttf_min(nz,enodes(1)), fct_ttf_min(nz,enodes(2)), fct_ttf_min(nz,enodes(3))) - end do - !$ACC END LOOP - if (nl1<=nl-1) then - !$ACC LOOP VECTOR - do nz=nl1,nl-1 - AUX(1,nz,elem)=-bignumber - AUX(2,nz,elem)= bignumber - end do - !$ACC END LOOP - endif + do nz=1, nl - 1 + enodes=elem2D_nodes(:,elem) + nl1 = nlevels(elem) + nu1 = ulevels(elem) + if(nu1 <= nz) then + if(nz < nl1) then + ! do nz=nu1, nl1-1 + AUX(1,nz,elem)=max(fct_ttf_max(nz,enodes(1)), fct_ttf_max(nz,enodes(2)), fct_ttf_max(nz,enodes(3))) + AUX(2,nz,elem)=min(fct_ttf_min(nz,enodes(1)), fct_ttf_min(nz,enodes(2)), fct_ttf_min(nz,enodes(3))) + ! end do + else + ! do nz=nl1,nl-1 + AUX(1,nz,elem)=-bignumber + AUX(2,nz,elem)= bignumber + ! end do + endif + endif + end do ! --> do elem=1, myDim_elem2D end do ! --> do elem=1, myDim_elem2D !$ACC END PARALLEL LOOP !$OMP END DO @@ -163,56 +166,62 @@ subroutine oce_tra_adv_fct(dt, ttf, lo, adf_h, adf_v, fct_ttf_min, fct_ttf_max, ! vertical gradients are larger. !Horizontal !$OMP DO - !$ACC PARALLEL LOOP GANG DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) + !$ACC PARALLEL LOOP GANG VECTOR COLLAPSE(2) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) do n=1, myDim_nod2D - nu1 = ulevels_nod2D(n) - nl1 = nlevels_nod2D(n) - !___________________________________________________________________ - !$ACC LOOP VECTOR - do nz=nu1,nl1-1 - ! max,min horizontal bound in cluster around node n in every - ! vertical layer - ! nod_in_elem2D --> elem indices of which node n is surrounded - ! nod_in_elem2D_num --> max number of surrounded elem - tvert_max(nz, n) = AUX(1,nz, nod_in_elem2D(1, n)) - tvert_min(nz, n) = AUX(2,nz, nod_in_elem2D(1, n)) - !$ACC LOOP SEQ - do elem=2,nod_in_elem2D_num(n) - tvert_max(nz, n) = dmax1(tvert_max(nz, n), AUX(1,nz, nod_in_elem2D(elem,n))) - tvert_min(nz, n) = dmin1(tvert_min(nz, n), AUX(2,nz, nod_in_elem2D(elem,n))) - end do - !$ACC END LOOP - end do - !$ACC END LOOP + do nz=1, nl + nu1 = ulevels_nod2D(n) + nl1 = nlevels_nod2D(n) + if (nu1 <= nz .and. nz < nl1) then + !___________________________________________________________________ + ! do nz=nu1,nl1-1 + ! max,min horizontal bound in cluster around node n in every + ! vertical layer + ! nod_in_elem2D --> elem indices of which node n is surrounded + ! nod_in_elem2D_num --> max number of surrounded elem + tvert_max(nz, n) = AUX(1,nz, nod_in_elem2D(1, n)) + tvert_min(nz, n) = AUX(2,nz, nod_in_elem2D(1, n)) + !$ACC LOOP SEQ + do elem=2,nod_in_elem2D_num(n) + tvert_max(nz, n) = dmax1(tvert_max(nz, n), AUX(1,nz, nod_in_elem2D(elem,n))) + tvert_min(nz, n) = dmin1(tvert_min(nz, n), AUX(2,nz, nod_in_elem2D(elem,n))) + end do + !$ACC END LOOP + ! end do + end if + end do end do !$ACC END PARALLEL LOOP !$OMP END DO !$OMP DO - !$ACC PARALLEL LOOP GANG DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) + !$ACC PARALLEL LOOP GANG VECTOR COLLAPSE(2) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) do n=1, myDim_nod2D - nu1 = ulevels_nod2D(n) - nl1 = nlevels_nod2D(n) - !___________________________________________________________________ - ! calc max,min increment of surface layer with respect to low order - ! solution - fct_ttf_max(nu1,n)=tvert_max(nu1, n)-LO(nu1,n) - fct_ttf_min(nu1,n)=tvert_min(nu1, n)-LO(nu1,n) - ! calc max,min increment from nz-1:nz+1 with respect to low order - ! solution at layer nz - !$ACC LOOP VECTOR - do nz=nu1+1,nl1-2 - fct_ttf_max(nz,n)=dmax1(tvert_max(nz-1, n), tvert_max(nz, n), tvert_max(nz+1, n))-LO(nz,n) - fct_ttf_min(nz,n)=dmin1(tvert_min(nz-1, n), tvert_min(nz, n), tvert_min(nz+1, n))-LO(nz,n) - end do - !$ACC END LOOP - ! calc max,min increment of bottom layer -1 with respect to low order - ! solution - nz=nl1-1 - fct_ttf_max(nz,n)=tvert_max(nz, n)-LO(nz,n) - fct_ttf_min(nz,n)=tvert_min(nz, n)-LO(nz,n) + do nz=1, nl + nu1 = ulevels_nod2D(n) + nl1 = nlevels_nod2D(n) + !___________________________________________________________________ + ! calc max,min increment of surface layer with respect to low order + ! solution + if (nu1 == nz) then + fct_ttf_max(nz,n)=tvert_max(nz, n)-LO(nz,n) + fct_ttf_min(nz,n)=tvert_min(nz, n)-LO(nz,n) + else if (nu1 < nz .and. nz < nl1 - 1) then + ! calc max,min increment from nz-1:nz+1 with respect to low order + ! solution at layer nz + ! do nz=nu1+1,nl1-2 + fct_ttf_max(nz,n)=dmax1(tvert_max(nz-1, n), tvert_max(nz, n), tvert_max(nz+1, n))-LO(nz,n) + fct_ttf_min(nz,n)=dmin1(tvert_min(nz-1, n), tvert_min(nz, n), tvert_min(nz+1, n))-LO(nz,n) + ! end do + ! calc max,min increment of bottom layer -1 with respect to low order + ! solution + else if (nz == nl1 - 1) then + ! nz=nl1-1 + fct_ttf_max(nz,n)=tvert_max(nz, n)-LO(nz,n) + fct_ttf_min(nz,n)=tvert_min(nz, n)-LO(nz,n) + end if + end do end do - !$ACC END PARALLEL LOOP + !$ACC END PARALLEL LOOP !$OMP END DO !___________________________________________________________________________ ! b1. Split positive and negative antidiffusive contributions @@ -220,33 +229,36 @@ subroutine oce_tra_adv_fct(dt, ttf, lo, adf_h, adf_v, fct_ttf_min, fct_ttf_max, ! horizontal element and vertical node contribution to node n and layer nz ! see. R. Löhner et al. "finite element flux corrected transport (FEM-FCT) ! for the euler and navier stoke equation -!$OMP DO - !$ACC PARALLEL LOOP GANG DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) - do n=1, myDim_nod2D - nu1 = ulevels_nod2D(n) - nl1 = nlevels_nod2D(n) - !$ACC LOOP VECTOR - do nz=nu1,nl1-1 - fct_plus(nz,n)=0._WP - fct_minus(nz,n)=0._WP - end do - !$ACC END LOOP - end do - !$ACC END PARALLEL LOOP -!$OMP END DO - +!!$OMP DO + ! !$ACC PARALLEL LOOP GANG VECTOR COLLAPSE(2) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) + ! do n=1, myDim_nod2D + ! do nz=1, nl + ! nu1 = ulevels_nod2D(n) + ! nl1 = nlevels_nod2D(n) + ! if(nu1 <= nz .and. nz < nl1) then + ! ! do nz=nu1,nl1-1 + ! fct_plus(nz,n)=0._WP + ! fct_minus(nz,n)=0._WP + ! ! end do + ! end if + ! end do + ! end do + ! !$ACC END PARALLEL LOOP +!!$OMP END DO !Vertical !$OMP DO - !$ACC PARALLEL LOOP GANG DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) + !$ACC PARALLEL LOOP GANG VECTOR COLLAPSE(2) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) do n=1, myDim_nod2D - nu1 = ulevels_nod2D(n) - nl1 = nlevels_nod2D(n) - !$ACC LOOP VECTOR - do nz=nu1,nl1-1 - fct_plus(nz,n) =fct_plus(nz,n) +(max(0.0_WP,adf_v(nz,n))+max(0.0_WP,-adf_v(nz+1,n))) - fct_minus(nz,n)=fct_minus(nz,n)+(min(0.0_WP,adf_v(nz,n))+min(0.0_WP,-adf_v(nz+1,n))) - end do - !$ACC END LOOP + do nz=1, nl + nu1 = ulevels_nod2D(n) + nl1 = nlevels_nod2D(n) + if(nu1 <= nz .and. nz < nl1) then + ! do nz=nu1,nl1-1 + fct_plus(nz,n) = 0._WP + (max(0.0_WP,adf_v(nz,n)) + max(0.0_WP,-adf_v(nz+1,n))) + fct_minus(nz,n) = 0._WP + (min(0.0_WP,adf_v(nz,n)) + min(0.0_WP,-adf_v(nz+1,n))) + ! end do + end if + end do end do !$ACC END PARALLEL LOOP !$OMP END DO @@ -254,65 +266,65 @@ subroutine oce_tra_adv_fct(dt, ttf, lo, adf_h, adf_v, fct_ttf_min, fct_ttf_max, !$OMP DO !Horizontal #if !defined(DISABLE_OPENACC_ATOMICS) - !$ACC PARALLEL LOOP GANG PRIVATE(enodes, el) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) + !$ACC PARALLEL LOOP GANG VECTOR COLLAPSE(2) PRIVATE(enodes, el) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) #else !$ACC UPDATE SELF(fct_plus, fct_minus, adf_h) #endif do edge=1, myDim_edge2D - enodes(1:2)=edges(:,edge) - el=edge_tri(:,edge) - nl1=nlevels(el(1))-1 - nu1=ulevels(el(1)) - nl2=0 - nu2=0 - if (el(2)>0) then - nl2=nlevels(el(2))-1 - nu2=ulevels(el(2)) - end if + do nz=1, nl + enodes(1:2)=edges(:,edge) + el=edge_tri(:,edge) + nl1=nlevels(el(1))-1 + nu1=ulevels(el(1)) + nl2=0 + nu2=0 + if (el(2)>0) then + nl2=nlevels(el(2))-1 + nu2=ulevels(el(2)) + end if + + nl12 = max(nl1,nl2) + nu12 = nu1 + if (nu2>0) nu12 = min(nu1,nu2) + + if(nu12 <= nz .and. nz <= nl12) then - nl12 = max(nl1,nl2) - nu12 = nu1 - if (nu2>0) nu12 = min(nu1,nu2) #if defined(_OPENMP) && !defined(__openmp_reproducible) - call omp_set_lock(partit%plock(enodes(1))) + call omp_set_lock(partit%plock(enodes(1))) #else !$OMP ORDERED #endif + ! do nz=nu12, nl12 #if !defined(DISABLE_OPENACC_ATOMICS) - !$ACC LOOP VECTOR -#endif - do nz=nu12, nl12 -#if !defined(DISABLE_OPENACC_ATOMICS) - !$ACC ATOMIC UPDATE + !$ACC ATOMIC UPDATE #endif - fct_plus (nz,enodes(1))=fct_plus (nz,enodes(1)) + max(0.0_WP, adf_h(nz,edge)) + fct_plus (nz,enodes(1))=fct_plus (nz,enodes(1)) + max(0.0_WP, adf_h(nz,edge)) #if !defined(DISABLE_OPENACC_ATOMICS) - !$ACC ATOMIC UPDATE + !$ACC ATOMIC UPDATE #endif - fct_minus(nz,enodes(1))=fct_minus(nz,enodes(1)) + min(0.0_WP, adf_h(nz,edge)) + fct_minus(nz,enodes(1))=fct_minus(nz,enodes(1)) + min(0.0_WP, adf_h(nz,edge)) #if defined(_OPENMP) && !defined(__openmp_reproducible) - end do - call omp_unset_lock(partit%plock(enodes(1))) - call omp_set_lock (partit%plock(enodes(2))) - do nz=nu12, nl12 -#endif -#if !defined(DISABLE_OPENACC_ATOMICS) - !$ACC ATOMIC UPDATE + ! end do + call omp_unset_lock(partit%plock(enodes(1))) + call omp_set_lock (partit%plock(enodes(2))) + ! do nz=nu12, nl12 #endif - fct_plus (nz,enodes(2))=fct_plus (nz,enodes(2)) + max(0.0_WP,-adf_h(nz,edge)) #if !defined(DISABLE_OPENACC_ATOMICS) - !$ACC ATOMIC UPDATE + !$ACC ATOMIC UPDATE #endif - fct_minus(nz,enodes(2))=fct_minus(nz,enodes(2)) + min(0.0_WP,-adf_h(nz,edge)) - end do + fct_plus (nz,enodes(2))=fct_plus (nz,enodes(2)) + max(0.0_WP,-adf_h(nz,edge)) #if !defined(DISABLE_OPENACC_ATOMICS) - !$ACC END LOOP + !$ACC ATOMIC UPDATE #endif + fct_minus(nz,enodes(2))=fct_minus(nz,enodes(2)) + min(0.0_WP,-adf_h(nz,edge)) + ! end do #if defined(_OPENMP) && !defined(__openmp_reproducible) - call omp_unset_lock(partit%plock(enodes(2))) + call omp_unset_lock(partit%plock(enodes(2))) #else !$OMP END ORDERED #endif + end if + end do end do #if !defined(DISABLE_OPENACC_ATOMICS) !$ACC END PARALLEL LOOP @@ -320,22 +332,23 @@ subroutine oce_tra_adv_fct(dt, ttf, lo, adf_h, adf_v, fct_ttf_min, fct_ttf_max, !$ACC UPDATE DEVICE(fct_plus, fct_minus) #endif !$OMP END DO - !___________________________________________________________________________ ! b2. Limiting factors !$OMP DO - !$ACC PARALLEL LOOP GANG DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) + !$ACC PARALLEL LOOP GANG VECTOR COLLAPSE(2) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) do n=1,myDim_nod2D - nu1=ulevels_nod2D(n) - nl1=nlevels_nod2D(n) - !$ACC LOOP VECTOR - do nz=nu1,nl1-1 - flux=fct_plus(nz,n)*dt/areasvol(nz,n)+flux_eps - fct_plus(nz,n)=min(1.0_WP,fct_ttf_max(nz,n)/flux) - flux=fct_minus(nz,n)*dt/areasvol(nz,n)-flux_eps - fct_minus(nz,n)=min(1.0_WP,fct_ttf_min(nz,n)/flux) + do nz=1, nl + nu1=ulevels_nod2D(n) + nl1=nlevels_nod2D(n) + if (nu1 <= nz .and. nz < nl1) then + ! do nz=nu1,nl1-1 + flux=fct_plus(nz,n)*dt/areasvol(nz,n)+flux_eps + fct_plus(nz,n)=min(1.0_WP,fct_ttf_max(nz,n)/flux) + flux=fct_minus(nz,n)*dt/areasvol(nz,n)-flux_eps + fct_minus(nz,n)=min(1.0_WP,fct_ttf_min(nz,n)/flux) + ! end do + end if end do - !$ACC END LOOP end do !$ACC END PARALLEL LOOP !$OMP END DO @@ -350,78 +363,80 @@ subroutine oce_tra_adv_fct(dt, ttf, lo, adf_h, adf_v, fct_ttf_min, fct_ttf_max, ! b3. Limiting !Vertical !$OMP DO - !$ACC PARALLEL LOOP GANG DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) + !$ACC PARALLEL LOOP GANG VECTOR COLLAPSE(2) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) do n=1, myDim_nod2D - nu1=ulevels_nod2D(n) - nl1=nlevels_nod2D(n) - - !_______________________________________________________________________ - nz=nu1 - ae=1.0_WP - flux=adf_v(nz,n) - if(flux>=0.0_WP) then - ae=min(ae,fct_plus(nz,n)) - else - ae=min(ae,fct_minus(nz,n)) - end if - adf_v(nz,n)=ae*adf_v(nz,n) - - !_______________________________________________________________________ - !$ACC LOOP VECTOR - do nz=nu1+1,nl1-1 + do nz=1, nl + nu1=ulevels_nod2D(n) + nl1=nlevels_nod2D(n) ae=1.0_WP flux=adf_v(nz,n) - if(flux>=0._WP) then - ae=min(ae,fct_minus(nz-1,n)) - ae=min(ae,fct_plus(nz,n)) - else - ae=min(ae,fct_plus(nz-1,n)) - ae=min(ae,fct_minus(nz,n)) + + !_______________________________________________________________________ + if(nu1 == nz) then + ! nz=nu1 + if(flux>=0.0_WP) then + ae=min(ae,fct_plus(nz,n)) + else + ae=min(ae,fct_minus(nz,n)) + end if + adf_v(nz,n)=ae*adf_v(nz,n) + else if (nu1 < nz .and. nz < nl1) then + !_______________________________________________________________________ + ! do nz=nu1+1,nl1-1 + if(flux>=0._WP) then + ae=min(ae,fct_minus(nz-1,n)) + ae=min(ae,fct_plus(nz,n)) + else + ae=min(ae,fct_plus(nz-1,n)) + ae=min(ae,fct_minus(nz,n)) + end if + adf_v(nz,n)=ae*adf_v(nz,n) + ! end do end if - adf_v(nz,n)=ae*adf_v(nz,n) + ! the bottom flux is always zero end do - !$ACC END LOOP - ! the bottom flux is always zero end do !$ACC END PARALLEL LOOP !$OMP END DO !Horizontal !$OMP DO - - !$ACC PARALLEL LOOP GANG PRIVATE(enodes, el) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) + !$ACC PARALLEL LOOP GANG VECTOR COLLAPSE(2) PRIVATE(enodes, el) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) do edge=1, myDim_edge2D - enodes(1:2)=edges(:,edge) - el=edge_tri(:,edge) - nu1=ulevels(el(1)) - nl1=nlevels(el(1))-1 - nl2=0 - nu2=0 - if(el(2)>0) then - nu2=ulevels(el(2)) - nl2=nlevels(el(2))-1 - end if + do nz=1, nl + el=edge_tri(:,edge) + nu1=ulevels(el(1)) + nl1=nlevels(el(1))-1 + nl2=0 + nu2=0 + if(el(2)>0) then + nu2=ulevels(el(2)) + nl2=nlevels(el(2))-1 + end if - nl12 = max(nl1,nl2) - nu12 = nu1 - if (nu2>0) nu12 = min(nu1,nu2) + nl12 = max(nl1,nl2) + nu12 = nu1 + if (nu2>0) nu12 = min(nu1,nu2) - !$ACC LOOP VECTOR - do nz=nu12, nl12 - ae=1.0_WP - flux=adf_h(nz,edge) + if (nu12 <= nz .and. nz <= nl12) then + ! do nz=nu12, nl12 + enodes(1:2)=edges(:,edge) - if(flux>=0._WP) then - ae=min(ae,fct_plus(nz,enodes(1))) - ae=min(ae,fct_minus(nz,enodes(2))) - else - ae=min(ae,fct_minus(nz,enodes(1))) - ae=min(ae,fct_plus(nz,enodes(2))) - endif + ae=1.0_WP + flux=adf_h(nz,edge) - adf_h(nz,edge)=ae*adf_h(nz,edge) + if(flux>=0._WP) then + ae=min(ae,fct_plus(nz,enodes(1))) + ae=min(ae,fct_minus(nz,enodes(2))) + else + ae=min(ae,fct_minus(nz,enodes(1))) + ae=min(ae,fct_plus(nz,enodes(2))) + endif + + adf_h(nz,edge)=ae*adf_h(nz,edge) + ! end do + end if end do - !$ACC END LOOP end do !$ACC END PARALLEL LOOP !$OMP END DO diff --git a/src/oce_adv_tra_hor.F90 b/src/oce_adv_tra_hor.F90 index 8795f5a15..602bfd0a2 100644 --- a/src/oce_adv_tra_hor.F90 +++ b/src/oce_adv_tra_hor.F90 @@ -88,17 +88,17 @@ subroutine adv_tra_hor_upw1(vel, ttf, partit, mesh, flux, o_init_zero) if (present(o_init_zero)) then l_init_zero=o_init_zero end if - if (l_init_zero) then -!$OMP PARALLEL DO - !$ACC PARALLEL LOOP GANG VECTOR COLLAPSE(2) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) - do edge=1, myDim_edge2D - do nz=1, mesh%nl-1 - flux(nz,edge)=0.0_WP - end do - end do - !$ACC END PARALLEL LOOP -!$OMP END PARALLEL DO - end if +! if (l_init_zero) then +! !$OMP PARALLEL DO +! !$ACC PARALLEL LOOP GANG VECTOR COLLAPSE(2) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) +! do edge=1, myDim_edge2D +! do nz=1, mesh%nl-1 +! flux(nz,edge)=0.0_WP +! end do +! end do +! !$ACC END PARALLEL LOOP +! !$OMP END PARALLEL DO +! end if ! The result is the low-order solution horizontal fluxes ! They are put into flux @@ -106,128 +106,133 @@ subroutine adv_tra_hor_upw1(vel, ttf, partit, mesh, flux, o_init_zero) !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(edge, deltaX1, deltaY1, deltaX2, deltaY2, & !$OMP a, vflux, el, enodes, nz, nu12, nl12, nl1, nl2, nu1, nu2) !$OMP DO - !$ACC PARALLEL LOOP GANG PRIVATE(enodes, el) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) + !$ACC PARALLEL LOOP GANG VECTOR COLLAPSE(2) PRIVATE(enodes, el) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) do edge=1, myDim_edge2D - ! local indice of nodes that span up edge ed - enodes=edges(:,edge) - - ! local index of element that contribute to edge - el=edge_tri(:,edge) - - ! number of layers -1 at elem el(1) - nl1=nlevels(el(1))-1 - - ! index off surface layer in case of cavity !=1 - nu1=ulevels(el(1)) - - ! edge_cross_dxdy(1:2,ed)... dx,dy distance from element centroid el(1) to - ! center of edge --> needed to calc flux perpedicular to edge from elem el(1) - deltaX1=edge_cross_dxdy(1,edge) - deltaY1=edge_cross_dxdy(2,edge) - a=r_earth*elem_cos(el(1)) - - !_______________________________________________________________________ - ! same parameter but for other element el(2) that contributes to edge ed - ! if el(2)==0 than edge is boundary edge - nl2=0 - nu2=0 - if(el(2)>0) then - deltaX2=edge_cross_dxdy(3,edge) - deltaY2=edge_cross_dxdy(4,edge) - ! number of layers -1 at elem el(2) - nl2=nlevels(el(2))-1 - nu2=ulevels(el(2)) - a=0.5_WP*(a+r_earth*elem_cos(el(2))) - end if - - !_______________________________________________________________________ - ! nl12 ... minimum number of layers -1 between element el(1) & el(2) that - ! contribute to edge ed - ! nu12 ... upper index of layers between element el(1) & el(2) that - ! contribute to edge ed - ! be carefull !!! --> if ed is a boundary edge than el(1)~=0 and el(2)==0 - ! that means nl1>0, nl2==0, n2=min(nl1,nl2)=0 !!! - nl12=min(nl1,nl2) - nu12=max(nu1,nu2) - - !_______________________________________________________________________ - ! (A) goes only into this loop when the edge has only facing element - ! el(1) --> so the edge is a boundary edge --> this is for ocean - ! surface in case of cavity - !$ACC LOOP VECTOR - do nz=nu1, nu12-1 - !____________________________________________________________________ - ! volume flux across the segments - vflux=(-VEL(2,nz,el(1))*deltaX1 + VEL(1,nz,el(1))*deltaY1)*helem(nz,el(1)) - - !____________________________________________________________________ - ! 1st. low order upwind solution - flux(nz, edge)=-0.5_WP*( & - ttf(nz, enodes(1))*(vflux+abs(vflux))+ & - ttf(nz, enodes(2))*(vflux-abs(vflux)) & - )-flux(nz, edge) - end do - !$ACC END LOOP - - !_______________________________________________________________________ - ! (B) goes only into this loop when the edge has only facing elemenmt - ! el(2) --> so the edge is a boundary edge --> this is for ocean - ! surface in case of cavity - if (nu2 > 0) then - !$ACC LOOP VECTOR - do nz=nu2, nu12-1 - !___________________________________________________________ + do nz=1, nl + ! do nz=1, mesh%nl-1 + if (l_init_zero .and. nz < nl) flux(nz,edge)=0.0_WP + ! end do + + ! local indice of nodes that span up edge ed + enodes=edges(:,edge) + + ! local index of element that contribute to edge + el=edge_tri(:,edge) + + ! number of layers -1 at elem el(1) + nl1=nlevels(el(1))-1 + + ! index off surface layer in case of cavity !=1 + nu1=ulevels(el(1)) + + ! edge_cross_dxdy(1:2,ed)... dx,dy distance from element centroid el(1) to + ! center of edge --> needed to calc flux perpedicular to edge from elem el(1) + deltaX1=edge_cross_dxdy(1,edge) + deltaY1=edge_cross_dxdy(2,edge) + a=r_earth*elem_cos(el(1)) + + !_______________________________________________________________________ + ! same parameter but for other element el(2) that contributes to edge ed + ! if el(2)==0 than edge is boundary edge + nl2=0 + nu2=0 + if(el(2)>0) then + deltaX2=edge_cross_dxdy(3,edge) + deltaY2=edge_cross_dxdy(4,edge) + ! number of layers -1 at elem el(2) + nl2=nlevels(el(2))-1 + nu2=ulevels(el(2)) + a=0.5_WP*(a+r_earth*elem_cos(el(2))) + end if + + !_______________________________________________________________________ + ! nl12 ... minimum number of layers -1 between element el(1) & el(2) that + ! contribute to edge ed + ! nu12 ... upper index of layers between element el(1) & el(2) that + ! contribute to edge ed + ! be carefull !!! --> if ed is a boundary edge than el(1)~=0 and el(2)==0 + ! that means nl1>0, nl2==0, n2=min(nl1,nl2)=0 !!! + nl12=min(nl1,nl2) + nu12=max(nu1,nu2) + + !_______________________________________________________________________ + ! (A) goes only into this loop when the edge has only facing element + ! el(1) --> so the edge is a boundary edge --> this is for ocean + ! surface in case of cavity + if (nu1 <= nz .and. nz < nu12) then + ! do nz=nu1, nu12-1 + !____________________________________________________________________ ! volume flux across the segments - vflux=(VEL(2,nz,el(2))*deltaX2 - VEL(1,nz,el(2))*deltaY2)*helem(nz,el(2)) + vflux=(-VEL(2,nz,el(1))*deltaX1 + VEL(1,nz,el(1))*deltaY1)*helem(nz,el(1)) - !___________________________________________________________ + !____________________________________________________________________ ! 1st. low order upwind solution - flux(nz, edge)=-0.5_WP*( & - ttf(nz, enodes(1))*(vflux+abs(vflux))+ & - ttf(nz, enodes(2))*(vflux-abs(vflux)))-flux(nz, edge) - end do - !$ACC END LOOP - end if - - !_______________________________________________________________________ - ! (C) Both segments - ! loop over depth layers from top (nu12) to nl12 - ! be carefull !!! --> if ed is a boundary edge, el(2)==0 than nl12=0 so - ! you wont enter in this loop - !$ACC LOOP VECTOR - do nz=nu12, nl12 - !___________________________________________________________________ - ! 1st. low order upwind solution - ! here already assumed that ed is NOT! a boundary edge so el(2) should exist - vflux=(-VEL(2,nz,el(1))*deltaX1 + VEL(1,nz,el(1))*deltaY1)*helem(nz,el(1)) & - +(VEL(2,nz,el(2))*deltaX2 - VEL(1,nz,el(2))*deltaY2)*helem(nz,el(2)) - - flux(nz, edge)=-0.5_WP*( & - ttf(nz, enodes(1))*(vflux+abs(vflux))+ & - ttf(nz, enodes(2))*(vflux-abs(vflux)))-flux(nz, edge) - end do - !$ACC END LOOP - - !_______________________________________________________________________ - ! (D) remaining segments on the left or on the right - !$ACC LOOP VECTOR - do nz=nl12+1, nl1 - !____________________________________________________________________ - ! volume flux across the segments - vflux=(-VEL(2,nz,el(1))*deltaX1 + VEL(1,nz,el(1))*deltaY1)*helem(nz,el(1)) - !____________________________________________________________________ - ! 1st. low order upwind solution - flux(nz, edge)=-0.5_WP*( & - ttf(nz, enodes(1))*(vflux+abs(vflux))+ & - ttf(nz, enodes(2))*(vflux-abs(vflux)) & - )-flux(nz, edge) - end do - !$ACC END LOOP + flux(nz, edge)=-0.5_WP*( & + ttf(nz, enodes(1))*(vflux+abs(vflux))+ & + ttf(nz, enodes(2))*(vflux-abs(vflux)) & + )-flux(nz, edge) + ! end do + end if + + !_______________________________________________________________________ + ! (B) goes only into this loop when the edge has only facing elemenmt + ! el(2) --> so the edge is a boundary edge --> this is for ocean + ! surface in case of cavity + if (nu2 > 0) then + if (nu2 <= nz .and. nz < nu12) then + ! do nz=nu2, nu12-1 + !___________________________________________________________ + ! volume flux across the segments + vflux=(VEL(2,nz,el(2))*deltaX2 - VEL(1,nz,el(2))*deltaY2)*helem(nz,el(2)) + + !___________________________________________________________ + ! 1st. low order upwind solution + flux(nz, edge)=-0.5_WP*( & + ttf(nz, enodes(1))*(vflux+abs(vflux))+ & + ttf(nz, enodes(2))*(vflux-abs(vflux)))-flux(nz, edge) + ! end do + end if + end if + + !_______________________________________________________________________ + ! (C) Both segments + ! loop over depth layers from top (nu12) to nl12 + ! be carefull !!! --> if ed is a boundary edge, el(2)==0 than nl12=0 so + ! you wont enter in this loop + if (nu12 <= nz .and. nz <= nl12) then + ! do nz=nu12, nl12 + !___________________________________________________________________ + ! 1st. low order upwind solution + ! here already assumed that ed is NOT! a boundary edge so el(2) should exist + vflux=(-VEL(2,nz,el(1))*deltaX1 + VEL(1,nz,el(1))*deltaY1)*helem(nz,el(1)) & + +(VEL(2,nz,el(2))*deltaX2 - VEL(1,nz,el(2))*deltaY2)*helem(nz,el(2)) - !_______________________________________________________________________ - ! (E) remaining segments on the left or on the right - !$ACC LOOP VECTOR - do nz=nl12+1, nl2 + flux(nz, edge)=-0.5_WP*( & + ttf(nz, enodes(1))*(vflux+abs(vflux))+ & + ttf(nz, enodes(2))*(vflux-abs(vflux)))-flux(nz, edge) + ! end do + end if + + !_______________________________________________________________________ + ! (D) remaining segments on the left or on the right + if(nl12 < nz .and. nz <= nl1) then + ! do nz=nl12+1, nl1 + !____________________________________________________________________ + ! volume flux across the segments + vflux=(-VEL(2,nz,el(1))*deltaX1 + VEL(1,nz,el(1))*deltaY1)*helem(nz,el(1)) + !____________________________________________________________________ + ! 1st. low order upwind solution + flux(nz, edge)=-0.5_WP*( & + ttf(nz, enodes(1))*(vflux+abs(vflux))+ & + ttf(nz, enodes(2))*(vflux-abs(vflux)) & + )-flux(nz, edge) + ! end do + end if + + !_______________________________________________________________________ + ! (E) remaining segments on the left or on the right + if (nl12 < nz .and. nz <= nl2) then + ! do nz=nl12+1, nl2 !_______________________________________________________________ ! volume flux across the segments vflux=(VEL(2,nz,el(2))*deltaX2 - VEL(1,nz,el(2))*deltaY2)*helem(nz,el(2)) @@ -236,8 +241,9 @@ subroutine adv_tra_hor_upw1(vel, ttf, partit, mesh, flux, o_init_zero) flux(nz, edge)=-0.5_WP*( & ttf(nz, enodes(1))*(vflux+abs(vflux))+ & ttf(nz, enodes(2))*(vflux-abs(vflux)))-flux(nz, edge) + ! end do + end if end do - !$ACC END LOOP end do !$ACC END PARALLEL LOOP !$OMP END DO @@ -578,231 +584,221 @@ subroutine adv_tra_hor_mfct(vel, ttf, partit, mesh, num_ord, flux, edge_up_dn_gr !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(edge, deltaX1, deltaY1, deltaX2, deltaY2, Tmean1, Tmean2, cHO, & !$OMP a, vflux, el, enodes, nz, nu12, nl12, nl1, nl2, nu1, nu2) !$OMP DO - !$ACC PARALLEL LOOP GANG PRIVATE(enodes, el) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) + !$ACC PARALLEL LOOP GANG VECTOR COLLAPSE(2) PRIVATE(enodes, el) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) do edge=1, myDim_edge2D - ! local indice of nodes that span up edge ed - enodes=edges(:,edge) - - ! local index of element that contribute to edge - el=edge_tri(:,edge) - - ! number of layers -1 at elem el(1) - nl1=nlevels(el(1))-1 - - ! index off surface layer in case of cavity !=1 - nu1=ulevels(el(1)) - - ! edge_cross_dxdy(1:2,ed)... dx,dy distance from element centroid el(1) to - ! center of edge --> needed to calc flux perpedicular to edge from elem el(1) - deltaX1=edge_cross_dxdy(1,edge) - deltaY1=edge_cross_dxdy(2,edge) - a=r_earth*elem_cos(el(1)) - - !_______________________________________________________________________ - ! same parameter but for other element el(2) that contributes to edge ed - ! if el(2)==0 than edge is boundary edge - nl2=0 - nu2=0 - if(el(2)>0) then - deltaX2=edge_cross_dxdy(3,edge) - deltaY2=edge_cross_dxdy(4,edge) - ! number of layers -1 at elem el(2) - nl2=nlevels(el(2))-1 - nu2=ulevels(el(2)) - a=0.5_WP*(a+r_earth*elem_cos(el(2))) - end if - - !_______________________________________________________________________ - ! n2 ... minimum number of layers -1 between element el(1) & el(2) that - ! contribute to edge ed - ! nu12 ... upper index of layers between element el(1) & el(2) that - ! contribute to edge ed - ! be carefull !!! --> if ed is a boundary edge than el(1)~=0 and el(2)==0 - ! that means nl1>0, nl2==0, n2=min(nl1,nl2)=0 !!! - nl12=min(nl1,nl2) - nu12=max(nu1,nu2) - - !_______________________________________________________________________ - ! (A) goes only into this loop when the edge has only facing element - ! el(1) --> so the edge is a boundary edge --> this is for ocean - ! surface in case of cavity - !$ACC LOOP VECTOR - do nz=nu1, nu12-1 - !____________________________________________________________________ - Tmean2=ttf(nz, enodes(2))- & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(2,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(4,nz,edge))/6.0_WP - - Tmean1=ttf(nz, enodes(1))+ & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(1,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(3,nz,edge))/6.0_WP - - !____________________________________________________________________ - ! volume flux across the segments - vflux=(-VEL(2,nz,el(1))*deltaX1 + VEL(1,nz,el(1))*deltaY1)*helem(nz,el(1)) - cHO=(vflux+abs(vflux))*Tmean1 + (vflux-abs(vflux))*Tmean2 - flux(nz,edge)=-0.5_WP*(1.0_WP-num_ord)*cHO - vflux*num_ord*0.5_WP*(Tmean1+Tmean2)-flux(nz,edge) - end do - !$ACC END LOOP - - !_______________________________________________________________________ - ! (B) goes only into this loop when the edge has only facing elemenmt - ! el(2) --> so the edge is a boundary edge --> this is for ocean - ! surface in case of cavity - if (nu2 > 0) then - !$ACC LOOP VECTOR - do nz=nu2,nu12-1 - !___________________________________________________________________ - Tmean2=ttf(nz, enodes(2))- & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(2,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(4,nz,edge))/6.0_WP - - Tmean1=ttf(nz, enodes(1))+ & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(1,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(3,nz,edge))/6.0_WP - !___________________________________________________________________ - ! volume flux across the segments - vflux=(VEL(2,nz,el(2))*deltaX2 - VEL(1,nz,el(2))*deltaY2)*helem(nz,el(2)) - cHO=(vflux+abs(vflux))*Tmean1 + (vflux-abs(vflux))*Tmean2 - flux(nz,edge)=-0.5_WP*(1.0_WP-num_ord)*cHO - vflux*num_ord*0.5_WP*(Tmean1+Tmean2)-flux(nz,edge) - end do - !$ACC END LOOP - end if - - !_______________________________________________________________________ - ! (C) Both segments - ! loop over depth layers from top to n2 - ! be carefull !!! --> if ed is a boundary edge, el(2)==0 than n2=0 so - ! you wont enter in this loop - !$ACC LOOP VECTOR - do nz=nu12, nl12 - !___________________________________________________________________ - ! MUSCL-type reconstruction - ! check if upwind or downwind triagle is necessary - ! - ! cross product between velocity vector and cross vector edge-elem-center - ! cross product > 0 --> angle vec_v and (dx,dy) --> [0 180] --> upwind triangle - ! cross product < 0 --> angle vec_v and (dx,dy) --> [180 360] --> downwind triangle - ! - ! o o ! o o - ! / \ / \ ! / \ / \ - ! / \ \ vec_v / \ ! / \ / / \ - ! / up \ \ / dn \ ! / up \ / / dn \ - ! o-------o----+---->o-------o ! o-------o----+---->o-------o - ! 1 / 2 ! 1 \vec_v - ! /vec_v ! \ - ! --> downwind triangle ! --> upwind triangle - ! - ! edge_up_dn_grad(1,nz,edge) ... gradTR_x upwind - ! edge_up_dn_grad(2,nz,edge) ... gradTR_x downwind - ! edge_up_dn_grad(3,nz,edge) ... gradTR_y upwind - ! edge_up_dn_grad(4,nz,edge) ... gradTR_y downwind - - !___________________________________________________________________ - ! use downwind triangle to interpolate Tracer to edge center with - ! fancy scheme --> Linear upwind reconstruction - ! T_n+0.5 = T_n+1 - 1/2*deltax*GRADIENT - ! --> GRADIENT = 2/3 GRAD_edgecenter + 1/3 GRAD_downwindtri - ! T_n+0.5 = T_n+1 - 2/6*(T_n+1-T_n) + 1/6*gradT_down - ! --> edge_up_dn_grad ... contains already elemental tracer gradient - ! of up and dn wind triangle - ! --> Tmean2 ... edge center interpolated Tracer using tracer - ! gradient info from upwind triangle - Tmean2=ttf(nz, enodes(2))- & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(2,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(4,nz,edge))/6.0_WP + do nz = 1, nl + ! local indice of nodes that span up edge ed + enodes=edges(:,edge) + + ! local index of element that contribute to edge + el=edge_tri(:,edge) + + ! number of layers -1 at elem el(1) + nl1=nlevels(el(1))-1 + + ! index off surface layer in case of cavity !=1 + nu1=ulevels(el(1)) + + ! edge_cross_dxdy(1:2,ed)... dx,dy distance from element centroid el(1) to + ! center of edge --> needed to calc flux perpedicular to edge from elem el(1) + deltaX1=edge_cross_dxdy(1,edge) + deltaY1=edge_cross_dxdy(2,edge) + a=r_earth*elem_cos(el(1)) + + !_______________________________________________________________________ + ! same parameter but for other element el(2) that contributes to edge ed + ! if el(2)==0 than edge is boundary edge + nl2=0 + nu2=0 + if(el(2)>0) then + deltaX2=edge_cross_dxdy(3,edge) + deltaY2=edge_cross_dxdy(4,edge) + ! number of layers -1 at elem el(2) + nl2=nlevels(el(2))-1 + nu2=ulevels(el(2)) + a=0.5_WP*(a+r_earth*elem_cos(el(2))) + end if + + !_______________________________________________________________________ + ! n2 ... minimum number of layers -1 between element el(1) & el(2) that + ! contribute to edge ed + ! nu12 ... upper index of layers between element el(1) & el(2) that + ! contribute to edge ed + ! be carefull !!! --> if ed is a boundary edge than el(1)~=0 and el(2)==0 + ! that means nl1>0, nl2==0, n2=min(nl1,nl2)=0 !!! + nl12=min(nl1,nl2) + nu12=max(nu1,nu2) + + !_______________________________________________________________________ + ! (A) goes only into this loop when the edge has only facing element + ! el(1) --> so the edge is a boundary edge --> this is for ocean + ! surface in case of cavity + if(nu1 <= nz .and. nz < nu12) then + !____________________________________________________________________ + Tmean2=ttf(nz, enodes(2))- & + (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & + edge_dxdy(1,edge)*a*edge_up_dn_grad(2,nz,edge)+ & + edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(4,nz,edge))/6.0_WP - ! use upwind triangle to interpolate Tracer to edge center with - ! fancy scheme --> Linear upwind reconstruction - ! T_n+0.5 = T_n + 1/2*deltax*GRADIENT - ! --> GRADIENT = 2/3 GRAD_edgecenter + 1/3 GRAD_downwindtri - ! T_n+0.5 = T_n + 2/6*(T_n+1-T_n) + 1/6*gradT_down - ! --> Tmean1 ... edge center interpolated Tracer using tracer - ! gradient info from downwind triangle - Tmean1=ttf(nz, enodes(1))+ & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(1,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(3,nz,edge))/6.0_WP - !___________________________________________________________________ - ! volume flux along the edge segment ed - ! netto volume flux along segment that comes from edge node 1 and 2 - ! - ! - ! C1 (centroid el(1)) --> (u1,v1) - ! x - ! ^ - ! (dx1,dy1) | - ! |---> vec_n1 (dy1,-dx1)--> project vec_u1 onto vec_n1 --> -v1*dx1+u1*dy1 --> - ! | | - ! enodes(1) o----------O---------o enodes(2) |-> calculate volume flux out of/in - ! vflux_________/| | the volume of enode1(enode2) through - ! |---> vec_n2 (dy2,-dx2)--> project vec_u2 onto vec_n2 --> -v2*dx2+u2*dy2 --> sections of dx1,dy1 and dx2,dy2 - ! (dx2,dy2) | --> vflux - ! v - ! x - ! C2 (centroid el(2)) --> (u2,v2) + Tmean1=ttf(nz, enodes(1))+ & + (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & + edge_dxdy(1,edge)*a*edge_up_dn_grad(1,nz,edge)+ & + edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(3,nz,edge))/6.0_WP - ! here already assumed that ed is NOT! a boundary edge so el(2) should exist - vflux=(-VEL(2,nz,el(1))*deltaX1 + VEL(1,nz,el(1))*deltaY1)*helem(nz,el(1)) & - +(VEL(2,nz,el(2))*deltaX2 - VEL(1,nz,el(2))*deltaY2)*helem(nz,el(2)) + !____________________________________________________________________ + ! volume flux across the segments + vflux=(-VEL(2,nz,el(1))*deltaX1 + VEL(1,nz,el(1))*deltaY1)*helem(nz,el(1)) + cHO=(vflux+abs(vflux))*Tmean1 + (vflux-abs(vflux))*Tmean2 + flux(nz,edge)=-0.5_WP*(1.0_WP-num_ord)*cHO - vflux*num_ord*0.5_WP*(Tmean1+Tmean2)-flux(nz,edge) + end if + + !_______________________________________________________________________ + ! (B) goes only into this loop when the edge has only facing elemenmt + ! el(2) --> so the edge is a boundary edge --> this is for ocean + ! surface in case of cavity + if(nu2 > 0 .and. nu2 <= nz .and. nz < nu12) then + !___________________________________________________________________ + Tmean2=ttf(nz, enodes(2))- & + (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & + edge_dxdy(1,edge)*a*edge_up_dn_grad(2,nz,edge)+ & + edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(4,nz,edge))/6.0_WP - !___________________________________________________________________ - ! (1-num_ord) is done with 3rd order upwind - cHO=(vflux+abs(vflux))*Tmean1 + (vflux-abs(vflux))*Tmean2 - flux(nz,edge)=-0.5_WP*(1.0_WP-num_ord)*cHO - vflux*num_ord*0.5_WP*(Tmean1+Tmean2)-flux(nz,edge) - end do - !$ACC END LOOP + Tmean1=ttf(nz, enodes(1))+ & + (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & + edge_dxdy(1,edge)*a*edge_up_dn_grad(1,nz,edge)+ & + edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(3,nz,edge))/6.0_WP + !___________________________________________________________________ + ! volume flux across the segments + vflux=(VEL(2,nz,el(2))*deltaX2 - VEL(1,nz,el(2))*deltaY2)*helem(nz,el(2)) + cHO=(vflux+abs(vflux))*Tmean1 + (vflux-abs(vflux))*Tmean2 + flux(nz,edge)=-0.5_WP*(1.0_WP-num_ord)*cHO - vflux*num_ord*0.5_WP*(Tmean1+Tmean2)-flux(nz,edge) + end if + + !_______________________________________________________________________ + ! (C) Both segments + ! loop over depth layers from top to n2 + ! be carefull !!! --> if ed is a boundary edge, el(2)==0 than n2=0 so + ! you wont enter in this loop + if (nu12 <= nz .and. nz <= nl12) then + !___________________________________________________________________ + ! MUSCL-type reconstruction + ! check if upwind or downwind triagle is necessary + ! + ! cross product between velocity vector and cross vector edge-elem-center + ! cross product > 0 --> angle vec_v and (dx,dy) --> [0 180] --> upwind triangle + ! cross product < 0 --> angle vec_v and (dx,dy) --> [180 360] --> downwind triangle + ! + ! o o ! o o + ! / \ / \ ! / \ / \ + ! / \ \ vec_v / \ ! / \ / / \ + ! / up \ \ / dn \ ! / up \ / / dn \ + ! o-------o----+---->o-------o ! o-------o----+---->o-------o + ! 1 / 2 ! 1 \vec_v + ! /vec_v ! \ + ! --> downwind triangle ! --> upwind triangle + ! + ! edge_up_dn_grad(1,nz,edge) ... gradTR_x upwind + ! edge_up_dn_grad(2,nz,edge) ... gradTR_x downwind + ! edge_up_dn_grad(3,nz,edge) ... gradTR_y upwind + ! edge_up_dn_grad(4,nz,edge) ... gradTR_y downwind + + !___________________________________________________________________ + ! use downwind triangle to interpolate Tracer to edge center with + ! fancy scheme --> Linear upwind reconstruction + ! T_n+0.5 = T_n+1 - 1/2*deltax*GRADIENT + ! --> GRADIENT = 2/3 GRAD_edgecenter + 1/3 GRAD_downwindtri + ! T_n+0.5 = T_n+1 - 2/6*(T_n+1-T_n) + 1/6*gradT_down + ! --> edge_up_dn_grad ... contains already elemental tracer gradient + ! of up and dn wind triangle + ! --> Tmean2 ... edge center interpolated Tracer using tracer + ! gradient info from upwind triangle + Tmean2=ttf(nz, enodes(2))- & + (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & + edge_dxdy(1,edge)*a*edge_up_dn_grad(2,nz,edge)+ & + edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(4,nz,edge))/6.0_WP + + ! use upwind triangle to interpolate Tracer to edge center with + ! fancy scheme --> Linear upwind reconstruction + ! T_n+0.5 = T_n + 1/2*deltax*GRADIENT + ! --> GRADIENT = 2/3 GRAD_edgecenter + 1/3 GRAD_downwindtri + ! T_n+0.5 = T_n + 2/6*(T_n+1-T_n) + 1/6*gradT_down + ! --> Tmean1 ... edge center interpolated Tracer using tracer + ! gradient info from downwind triangle + Tmean1=ttf(nz, enodes(1))+ & + (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & + edge_dxdy(1,edge)*a*edge_up_dn_grad(1,nz,edge)+ & + edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(3,nz,edge))/6.0_WP + !___________________________________________________________________ + ! volume flux along the edge segment ed + ! netto volume flux along segment that comes from edge node 1 and 2 + ! + ! + ! C1 (centroid el(1)) --> (u1,v1) + ! x + ! ^ + ! (dx1,dy1) | + ! |---> vec_n1 (dy1,-dx1)--> project vec_u1 onto vec_n1 --> -v1*dx1+u1*dy1 --> + ! | | + ! enodes(1) o----------O---------o enodes(2) |-> calculate volume flux out of/in + ! vflux_________/| | the volume of enode1(enode2) through + ! |---> vec_n2 (dy2,-dx2)--> project vec_u2 onto vec_n2 --> -v2*dx2+u2*dy2 --> sections of dx1,dy1 and dx2,dy2 + ! (dx2,dy2) | --> vflux + ! v + ! x + ! C2 (centroid el(2)) --> (u2,v2) + + ! here already assumed that ed is NOT! a boundary edge so el(2) should exist + vflux=(-VEL(2,nz,el(1))*deltaX1 + VEL(1,nz,el(1))*deltaY1)*helem(nz,el(1)) & + +(VEL(2,nz,el(2))*deltaX2 - VEL(1,nz,el(2))*deltaY2)*helem(nz,el(2)) + + !___________________________________________________________________ + ! (1-num_ord) is done with 3rd order upwind + cHO=(vflux+abs(vflux))*Tmean1 + (vflux-abs(vflux))*Tmean2 + flux(nz,edge)=-0.5_WP*(1.0_WP-num_ord)*cHO - vflux*num_ord*0.5_WP*(Tmean1+Tmean2)-flux(nz,edge) + end if - !_______________________________________________________________________ - ! (D) remaining segments on the left or on the right - !$ACC LOOP VECTOR - do nz=nl12+1, nl1 - !____________________________________________________________________ - Tmean2=ttf(nz, enodes(2))- & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(2,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(4,nz,edge))/6.0_WP + !_______________________________________________________________________ + ! (D) remaining segments on the left or on the right + if(nl12 < nz .and. nz <= nl1) then + !____________________________________________________________________ + Tmean2=ttf(nz, enodes(2))- & + (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & + edge_dxdy(1,edge)*a*edge_up_dn_grad(2,nz,edge)+ & + edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(4,nz,edge))/6.0_WP - Tmean1=ttf(nz, enodes(1))+ & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(1,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(3,nz,edge))/6.0_WP + Tmean1=ttf(nz, enodes(1))+ & + (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & + edge_dxdy(1,edge)*a*edge_up_dn_grad(1,nz,edge)+ & + edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(3,nz,edge))/6.0_WP - !____________________________________________________________________ - ! volume flux across the segments - vflux=(-VEL(2,nz,el(1))*deltaX1 + VEL(1,nz,el(1))*deltaY1)*helem(nz,el(1)) - cHO=(vflux+abs(vflux))*Tmean1 + (vflux-abs(vflux))*Tmean2 - flux(nz,edge)=-0.5_WP*(1.0_WP-num_ord)*cHO - vflux*num_ord*0.5_WP*(Tmean1+Tmean2)-flux(nz,edge) - end do - !$ACC END LOOP + !____________________________________________________________________ + ! volume flux across the segments + vflux=(-VEL(2,nz,el(1))*deltaX1 + VEL(1,nz,el(1))*deltaY1)*helem(nz,el(1)) + cHO=(vflux+abs(vflux))*Tmean1 + (vflux-abs(vflux))*Tmean2 + flux(nz,edge)=-0.5_WP*(1.0_WP-num_ord)*cHO - vflux*num_ord*0.5_WP*(Tmean1+Tmean2)-flux(nz,edge) + end if - !_______________________________________________________________________ - ! (E) remaining segments on the left or on the right - !$ACC LOOP VECTOR - do nz=nl12+1, nl2 - !____________________________________________________________________ - Tmean2=ttf(nz, enodes(2))- & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(2,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(4,nz,edge))/6.0_WP + !_______________________________________________________________________ + ! (E) remaining segments on the left or on the right + if (nl12 < nz .and. nz <= nl2) then + !____________________________________________________________________ + Tmean2=ttf(nz, enodes(2))- & + (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & + edge_dxdy(1,edge)*a*edge_up_dn_grad(2,nz,edge)+ & + edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(4,nz,edge))/6.0_WP - Tmean1=ttf(nz, enodes(1))+ & - (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & - edge_dxdy(1,edge)*a*edge_up_dn_grad(1,nz,edge)+ & - edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(3,nz,edge))/6.0_WP + Tmean1=ttf(nz, enodes(1))+ & + (2.0_WP*(ttf(nz, enodes(2))-ttf(nz,enodes(1)))+ & + edge_dxdy(1,edge)*a*edge_up_dn_grad(1,nz,edge)+ & + edge_dxdy(2,edge)*r_earth*edge_up_dn_grad(3,nz,edge))/6.0_WP - !____________________________________________________________________ - ! volume flux across the segments - vflux=(VEL(2,nz,el(2))*deltaX2 - VEL(1,nz,el(2))*deltaY2)*helem(nz,el(2)) - cHO=(vflux+abs(vflux))*Tmean1 + (vflux-abs(vflux))*Tmean2 - flux(nz,edge)=-0.5_WP*(1.0_WP-num_ord)*cHO - vflux*num_ord*0.5_WP*(Tmean1+Tmean2)-flux(nz,edge) + !____________________________________________________________________ + ! volume flux across the segments + vflux=(VEL(2,nz,el(2))*deltaX2 - VEL(1,nz,el(2))*deltaY2)*helem(nz,el(2)) + cHO=(vflux+abs(vflux))*Tmean1 + (vflux-abs(vflux))*Tmean2 + flux(nz,edge)=-0.5_WP*(1.0_WP-num_ord)*cHO - vflux*num_ord*0.5_WP*(Tmean1+Tmean2)-flux(nz,edge) + end if end do - !$ACC END LOOP end do !$ACC END PARALLEL LOOP !$OMP END DO diff --git a/src/oce_adv_tra_ver.F90 b/src/oce_adv_tra_ver.F90 index a4d2daf74..bcd55e451 100644 --- a/src/oce_adv_tra_ver.F90 +++ b/src/oce_adv_tra_ver.F90 @@ -266,49 +266,52 @@ subroutine adv_tra_ver_upw1(w, ttf, partit, mesh, flux, o_init_zero) if (present(o_init_zero)) then l_init_zero=o_init_zero end if - if (l_init_zero) then -!$OMP PARALLEL DO - !$ACC PARALLEL LOOP GANG VECTOR COLLAPSE(2) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) - do n=1, myDim_nod2D - do nz=1,mesh%nl - flux(nz, n)=0.0_WP - end do - end do - !$ACC END PARALLEL LOOP -!$OMP END PARALLEL DO - end if +! if (l_init_zero) then +! !$OMP PARALLEL DO +! !$ACC PARALLEL LOOP GANG VECTOR COLLAPSE(2) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) +! do n=1, myDim_nod2D +! do nz=1,mesh%nl +! flux(nz, n)=0.0_WP +! end do +! end do +! !$ACC END PARALLEL LOOP +! !$OMP END PARALLEL DO +! end if !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(tvert, n, nz, nzmax, nzmin) !$OMP DO - !$ACC PARALLEL LOOP GANG DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) + !$ACC PARALLEL LOOP GANG VECTOR COLLAPSE(2) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) do n=1, myDim_nod2D - !_______________________________________________________________________ - nzmax=nlevels_nod2D(n) - nzmin=ulevels_nod2D(n) - - !_______________________________________________________________________ - ! vert. flux at surface layer - nz=nzmin - flux(nz,n)=-W(nz,n)*ttf(nz,n)*area(nz,n)-flux(nz,n) - - !_______________________________________________________________________ - ! vert. flux at bottom layer --> zero bottom flux - nz=nzmax - flux(nz,n)= 0.0_WP-flux(nz,n) - - !_______________________________________________________________________ - ! Be carefull have to do vertical tracer advection here on old vertical grid - ! also horizontal advection is done on old mesh (see helem contains old - ! mesh information) - !_______________________________________________________________________ - ! vert. flux at remaining levels - !$ACC LOOP VECTOR - do nz=nzmin+1,nzmax-1 - flux(nz,n)=-0.5*( & - ttf(nz ,n)*(W(nz,n)+abs(W(nz,n)))+ & - ttf(nz-1,n)*(W(nz,n)-abs(W(nz,n))))*area(nz,n)-flux(nz,n) - end do - !$ACC END LOOP + do nz=1, nl + !_______________________________________________________________________ + nzmax=nlevels_nod2D(n) + nzmin=ulevels_nod2D(n) + if(l_init_zero) flux(nz, n) = 0.0_WP + + !_______________________________________________________________________ + ! Be carefull have to do vertical tracer advection here on old vertical grid + ! also horizontal advection is done on old mesh (see helem contains old + ! mesh information) + !_______________________________________________________________________ + ! vert. flux at remaining levels + if ( nzmin < nz .and. nz < nzmax ) then + ! do nz=nzmin+1,nzmax-1 + flux(nz,n)=-0.5*( & + ttf(nz ,n)*(W(nz,n)+abs(W(nz,n)))+ & + ttf(nz-1,n)*(W(nz,n)-abs(W(nz,n))))*area(nz,n)-flux(nz,n) + ! end do + else if ( nzmin == nz ) then + !_______________________________________________________________________ + ! vert. flux at surface layer + ! nz=nzmin + flux(nz,n)=-W(nz,n)*ttf(nz,n)*area(nz,n)-flux(nz,n) + else if ( nzmax == nz ) then + !_______________________________________________________________________ + ! vert. flux at bottom layer --> zero bottom flux + ! nz=nzmax + flux(nz,n)= 0.0_WP-flux(nz,n) + end if + end do end do !$ACC END PARALLEL LOOP !$OMP END DO @@ -359,50 +362,49 @@ subroutine adv_tra_ver_qr4c(w, ttf, partit, mesh, num_ord, flux, o_init_zero) end if !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(tvert,n, nz, nzmax, nzmin, Tmean, Tmean1, Tmean2, qc, qu,qd) !$OMP DO - !$ACC PARALLEL LOOP GANG DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) + !$ACC PARALLEL LOOP GANG VECTOR COLLAPSE(2) DEFAULT(PRESENT) VECTOR_LENGTH(acc_vl) do n=1, myDim_nod2D - !_______________________________________________________________________ - nzmax=nlevels_nod2D(n) - nzmin=ulevels_nod2D(n) - !_______________________________________________________________________ - ! vert. flux at surface layer - nz=nzmin - flux(nz,n)=-ttf(nz,n)*W(nz,n)*area(nz,n)-flux(nz,n) - - !_______________________________________________________________________ - ! vert. flux 2nd layer --> centered differences - nz=nzmin+1 - flux(nz,n)=-0.5_WP*(ttf(nz-1,n)+ttf(nz,n))*W(nz,n)*area(nz,n)-flux(nz,n) - - !_______________________________________________________________________ - ! vert. flux at bottom - 1 layer --> centered differences - nz=nzmax-1 - flux(nz,n)=-0.5_WP*(ttf(nz-1,n)+ttf(nz,n))*W(nz,n)*area(nz,n)-flux(nz,n) - - !_______________________________________________________________________ - ! vert. flux at bottom layer --> zero bottom flux - nz=nzmax - flux(nz,n)= 0.0_WP-flux(nz,n) - - !_______________________________________________________________________ - ! Be carefull have to do vertical tracer advection here on old vertical grid - ! also horizontal advection is done on old mesh (see helem contains old - ! mesh information) - !_______________________________________________________________________ - ! vert. flux at remaining levels - !$ACC LOOP VECTOR - do nz=nzmin+2,nzmax-2 - !centered (4th order) - qc=(ttf(nz-1,n)-ttf(nz ,n))/(Z_3d_n(nz-1,n)-Z_3d_n(nz ,n)) - qu=(ttf(nz ,n)-ttf(nz+1,n))/(Z_3d_n(nz ,n)-Z_3d_n(nz+1,n)) - qd=(ttf(nz-2,n)-ttf(nz-1,n))/(Z_3d_n(nz-2,n)-Z_3d_n(nz-1,n)) - - Tmean1=ttf(nz ,n)+(2*qc+qu)*(zbar_3d_n(nz,n)-Z_3d_n(nz ,n))/3.0_WP - Tmean2=ttf(nz-1,n)+(2*qc+qd)*(zbar_3d_n(nz,n)-Z_3d_n(nz-1,n))/3.0_WP - Tmean =(W(nz,n)+abs(W(nz,n)))*Tmean1+(W(nz,n)-abs(W(nz,n)))*Tmean2 - flux(nz,n)=(-0.5_WP*(1.0_WP-num_ord)*Tmean - num_ord*(0.5_WP*(Tmean1+Tmean2))*W(nz,n))*area(nz,n)-flux(nz,n) - end do - !$ACC END LOOP + do nz=1, nl + !_______________________________________________________________________ + nzmax=nlevels_nod2D(n) + nzmin=ulevels_nod2D(n) + + !_______________________________________________________________________ + ! Be carefull have to do vertical tracer advection here on old vertical grid + ! also horizontal advection is done on old mesh (see helem contains old + ! mesh information) + !_______________________________________________________________________ + ! vert. flux at remaining levels + if (nzmin + 2 <= nz .and. nz <= nzmax - 2) then + ! do nz=nzmin+2,nzmax-2 + !centered (4th order) + qc=(ttf(nz-1,n)-ttf(nz ,n))/(Z_3d_n(nz-1,n)-Z_3d_n(nz ,n)) + qu=(ttf(nz ,n)-ttf(nz+1,n))/(Z_3d_n(nz ,n)-Z_3d_n(nz+1,n)) + qd=(ttf(nz-2,n)-ttf(nz-1,n))/(Z_3d_n(nz-2,n)-Z_3d_n(nz-1,n)) + + Tmean1=ttf(nz ,n)+(2*qc+qu)*(zbar_3d_n(nz,n)-Z_3d_n(nz ,n))/3.0_WP + Tmean2=ttf(nz-1,n)+(2*qc+qd)*(zbar_3d_n(nz,n)-Z_3d_n(nz-1,n))/3.0_WP + Tmean =(W(nz,n)+abs(W(nz,n)))*Tmean1+(W(nz,n)-abs(W(nz,n)))*Tmean2 + flux(nz,n)=(-0.5_WP*(1.0_WP-num_ord)*Tmean - num_ord*(0.5_WP*(Tmean1+Tmean2))*W(nz,n))*area(nz,n)-flux(nz,n) + ! end do + else if (nzmin + 1 == nz .or. nz == nzmax - 1) then + !_______________________________________________________________________ + ! vert. flux 2nd layer --> centered differences + !_______________________________________________________________________ + ! vert. flux at bottom - 1 layer --> centered differences + flux(nz,n)=-0.5_WP*(ttf(nz-1,n)+ttf(nz,n))*W(nz,n)*area(nz,n)-flux(nz,n) + else if( nzmin == nz ) then + !_______________________________________________________________________ + ! vert. flux at surface layer + ! nz=nzmin + flux(nz,n)=-ttf(nz,n)*W(nz,n)*area(nz,n)-flux(nz,n) + else if( nzmax == nz ) then + !_______________________________________________________________________ + ! vert. flux at bottom layer --> zero bottom flux + ! nz=nzmax + flux(nz,n)= 0.0_WP-flux(nz,n) + end if + end do end do !$ACC END PARALLEL LOOP !$OMP END DO