!### macro's ##################################################### ! #define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr #define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if #define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if ! #include "tm5.inc" ! !################################################################# module advecty use GO, only : gol, goPr, goErr implicit none ! --- interface private public :: advectyzoom, put_yedges contains subroutine advectyzoom(region) ! ! set parameters for advecty ! written by patrick berkvens and mike botchev, march-june 1999 ! updated and modified by MK, dec 2002 ! use dims, only : xref, yref, zref, tref, im, jm, lm use dims, only : status, nsplitsteps, n_operators, zoom2D use dims, only : touch_sp, touch_np, okdebug use dims, only : parent, splitorderzoom, zero, xcyc #ifdef with_budgets use budget_global, only : sum_advection #endif use global_data, only : wind_dat, mass_dat use toolbox, only : escape_tm #ifdef MPI use mpi_const, only : my_real,com_lev,ierr, & root_t,mpi_sum, pe_first_tracer use mpi_comm, only : barrier_k #endif use partools , only : myid, lmloc, ntracetloc implicit none ! input/output integer,intent(in) :: region ! local real,dimension(:,:,:),pointer :: bm real,dimension(:,:,:,:),pointer :: rm real,dimension(:,:),allocatable :: bm0,bm1 integer :: is,ie,ls,le,n,q integer :: imr,jmr,lmr,tref_,xref_,yref_,zref_ logical :: y_encountered character(len=1) :: dir real :: sum_old, sum_new,sum_old_all,sum_new_all ! start if ( ntracetloc==0) return !WP! only PE's with nonzero lmloc proceed call put_yedges(region) tref_ = tref(region)/tref(parent(region)) xref_ = xref(region)/xref(parent(region)) yref_ = yref(region)/yref(parent(region)) zref_ = zref(region)/zref(parent(region)) imr = im(region); jmr = jm(region); lmr = lm(region) allocate(bm0(0:im(region)+1,0:lmr+1)) allocate(bm1(0:im(region)+1,0:lmr+1)) bm => wind_dat(region)%bm_t rm => mass_dat(region)%rm_t sum_old = sum(rm(1:imr,1:jmr,1:lmr,1)) ! determine the scope for advecty: if ( region == 1 ) then xref_ = 0; zref_ = 0 ! to have is/ie and ls/le properly computed end if ! find q - the place in the splitorderzoom q=status(region)/((nsplitsteps/2)*tref_) ! corresponding to the begining of the last ! triple x-y-z of the parent y_encountered=.false. ! now track four following steps from q !wp! changed to four steps do n=1,n_operators dir=splitorderzoom(region,q*(nsplitsteps/2)*tref_+n) select case(dir) case('x') if ((.not.y_encountered).or.(xcyc(region) == 1)) then is=1 ! x-substep is before y => ie=im(region) ! full i-scope else is=xref_+1 ! y-substep is before x => ie=im(region)-xref_ ! restricted i-scope end if case('y') y_encountered=.true. case('z') if (.not.y_encountered) then ls=1 ! z-substep is before y => le=lmr ! full l-scope else ls=zref_+1 ! y-substep is before z => le=lmr-zref_ ! restricted l-scope end if if (zoom2D) then ls=1; le=lmr !WP! this is always the case end if case('c') case('v') case('s') case default print *,'advectyzoom: strange value in splitorderzoom(',region,',', & q*(nsplitsteps/2)*tref_+n,'): ',q call escape_tm('advectyzoom: Error in do_next3') end select end do #ifdef slopes #ifndef secmom if ( mod(status(region),(nsplitsteps/2)*tref_) >= n_operators ) then ! if (1) more than fourth substep and (2) slope scheme then ! at this substep no coarse fluxes will be applied ! for slopes only: zero out velocities via the y-edges of the region ! to assure that no fluxes will be applied if ( okdebug ) print *,'advectyzoom: zeroing out outer ', & 'mass fluxes (yref)',yref_ if ( touch_sp(region) /= 1) then bm0 = bm(:, yref_,:); bm(:, yref_,:) = zero else if ( okdebug) print *,'advectyzoom: skipping zeroing outer ', & 'mass flux at SP' end if if ( touch_np(region) /= 1) then bm1 = bm(:,jmr-yref_+2,:); bm(:,jmr-yref_+2,:) = zero else if ( okdebug) print *,'advectyzoom: skipping zeroing outer ', & 'mass flux at NP' end if end if #endif #endif call advecty_work(region,is,ie,ls,le) #ifdef slopes #ifndef secmom if (mod(status(region),(nsplitsteps/2)*tref_) >= n_operators ) then ! for slopes only: recreate velocities via the x-edges if ( touch_sp(region) /= 1 ) bm(:, yref_,:) = bm0 if ( touch_np(region) /= 1 ) bm(:,jmr-yref_+2,:) = bm1 end if #endif #endif #ifdef with_budgets sum_new = sum(rm(1:imr,1:jmr,1:lmr,1)) #ifdef MPI if ( myid==pe_first_tracer ) & sum_advection(region) = sum_advection(region) + sum_new - sum_old if ( okdebug .and. myid == root_t ) & print *, 'advectyzoom: advect_y, region, sum_advection ', & region,sum_advection(region),sum_new,sum_old #else sum_advection(region) = sum_advection(region) + sum_new - sum_old if ( okdebug) print *, 'advectyzoom: advect_y, region, sum_advection ', & region,sum_advection(region),sum_new,sum_old #endif #endif deallocate(bm0) deallocate(bm1) nullify(bm) end subroutine advectyzoom subroutine advecty_work(region,is,ie,ls,le) ! ! calls dynamv ! written by mike botchev, march-june 1999 ! use toolbox, only : escape_tm implicit none ! input/output integer,intent(in) :: region integer,intent(in) :: is integer,intent(in) :: ie integer,intent(in) :: ls integer,intent(in) :: le ! local #ifndef slopes call escape_tm('advecty: slopes advection not defined') #endif call dynamv(region,is,ie,ls,le) end subroutine advecty_work subroutine dynamv(region,is,ie,ls,le) !----------------------------------------------------------------------- ! !**** dynamv - south-north tracer transport v 9.1 ! ! programmed by mh mpi HH 23-feb-1994 ! ! purpose ! ------- ! calculate amount of tracer moved in a south-north advection ! substep ! ! interface ! --------- ! call dynamv ! ! method ! ------ ! slopes scheme ! ! externals ! --------- ! none ! ! reference ! --------- ! Russel and Lerner, 1979 !----------------------------------------------------------------------- ! fixed bug in calculating maximum courant number ! mh, Thu, Feb 24, 1994 12:49:27 ! included code for limits of slopes to prevent negative tracer ! masses mh, 20-jun-1994 ! ! zoom version written by mike botchev, march-june 1999 !----------------------------------------------------------------------- use dims, only : xref, yref, zref, tref, im, jm, lm use dims, only : okdebug, parent, nregions, xi, nxi use dims, only : touch_sp, touch_np, limits, nloop_max, zero, one use redgridZoom, only : nred, jred, imredj, clustsize use global_data, only : wind_dat use tracer_data, only : mass_dat use meteodata , only : m_dat use zoom_tools, only : mix_edges use toolbox, only : escape_tm #ifdef with_budgets use budget_global, only : budget_flux, jflux1, jflux2, apply_budget_global_flux #endif #ifdef MPI use mpi_comm, only : barrier_k #endif use partools , only : ntracetloc, myid, ntracet_ar use chem_param, only : ntracet, ra use omp_partools,only : my_omp_domain, my_omp_max, my_omp_sum implicit none ! input/output integer,intent(in) :: region integer,intent(in) :: is integer,intent(in) :: ie integer,intent(in) :: ls integer,intent(in) :: le ! local real,dimension(:,:,:,:),pointer :: rm,rxm,rym,rzm real,dimension(:,:,:), pointer :: m,bm real,dimension(:,:),allocatable :: mnew real,dimension(:,:),allocatable :: f,pf,fx,fz integer :: i,j,je,js,l,n,iee,iss integer :: imr,imr2,jmr,lmr integer :: tref_,xref_,yref_,zref_ real :: sfs,sfzs,sfn,sfzn real :: beta,mxval integer,dimension(3) :: mxloc integer :: iloop, nloop, nglob, offsetn integer,parameter :: max_nloop = 6 real,dimension(:,:), allocatable :: mx logical :: cfl_ok integer :: lrg, redfact, ixe, ixs real :: summ ! Omp parameters integer :: ns, ne, lss, lee real :: l_xi, g_xi, l_nxi, g_nxi, l_nl, g_nl ! start if ( okdebug ) print *,'dynamv: region=',region,' is,ie,ls,le=',is,ie,ls,le if ( ( region < 0 ) .or. ( region > nregions ) ) then call escape_tm('dynamv: STOP, illegal number of region !!!') end if ! compute refinement factors with respect to the parent xref_ = xref(region)/xref(parent(region)) yref_ = yref(region)/yref(parent(region)) zref_ = zref(region)/zref(parent(region)) tref_ = tref(region)/tref(parent(region)) imr=im(region);jmr=jm(region);lmr=lm(region) call compress_yedges(region,is,ie,ls,le) m => m_dat(region)%data rm => mass_dat(region)%rm_t rxm => mass_dat(region)%rxm_t rym => mass_dat(region)%rym_t rzm => mass_dat(region)%rzm_t bm => wind_dat(region)%bm_t ! check is,ie,ls,le: if ( ( is /= xref_+1 ) .and. ( is /= 1 ) ) & call escape_tm( 'dynamv: Wrong value for IS in dynamv') if ( ( ie /= imr-xref_ ) .and. ( ie /= imr ) ) & call escape_tm( 'dynamv: Wrong value for IE in dynamv') if ( ( ls /= zref_+1 ) .and. ( ls /= 1 ) ) & call escape_tm( 'dynamv: Wrong value for LS in dynamv') if ( ( le /= lmr-zref_ ) .and. ( le /= lmr ) ) & call escape_tm( 'dynamv: Wrong value for LE in dynamv') ! compute js/je -- cells is:ie,js:je:ls:le will be updated if (region==1) then js = 2 je = jmr-1 else js = yref_ je = jmr-yref_+1 if ( touch_sp(region) == 1) js = 2 !cmk....new option.. reduced grid if ( touch_np(region) == 1) je = jmr-1 !exclude poles end if ! if requested limit meridional slopes such that no negative if (limits) then !$OMP PARALLEL & !$OMP default (none) & !$OMP shared (is, ie, js, je, ls, le, ntracetloc, rm, rym) & !$OMP private (i, j, l, n, ns, ne) call my_omp_domain (1, ntracetloc, ns, ne) do n = ns, ne do l = ls, le do j = js - 1, je + 1 do i = is, ie ! Otherwise negative rym(i,j,l,n) = max (min (rym(i,j,l,n), rm(i,j,l,n)), & -rm(i,j,l,n)) end do end do end do end do !$OMP END PARALLEL end if ! loop over tracers and vertical layers !$OMP PARALLEL & !$OMP default (none) & #ifdef with_budgets !$OMP shared (apply_budget_global_flux, budget_flux, jflux1, jflux2,ra) & #endif !$OMP shared (region, is, ie, imr, js, je, jmr, ls, le, m, bm, & !$OMP okdebug, ntracetloc, ntracet_ar, myid, rm, rym, rxm, & !$OMP rzm, limits, g_nl, g_xi, g_nxi) & !$OMP private (i, j, l, lss, lee, l_nl, l_xi, l_nxi, f, pf, fx, fz, & !$OMP mnew, nloop, beta, mx, cfl_ok, offsetn, nglob) l_nl = 0.0 l_xi = 0.0 l_nxi = 0.0 allocate (f(im(region),0:jm(region))) allocate (pf(im(region),0:jm(region))) allocate (fx(im(region),0:jm(region))) allocate (fz(im(region),0:jm(region))) allocate (mnew(im(region),jm(region))) call my_omp_domain (ls, le, lss, lee) do l=lss,lee ! compute all inner fluxes ! f(*,j,*) is flux entering j-th cell from below?cmk nloop = 1 ! determine nloop per layer! beta = 2.0 allocate(mx(is:ie,js-1:je+1)) do while(abs(beta) >= one .and. nloop < max_nloop) mx(is:ie,js-1:je+1) = m(is:ie,js-1:je+1,l) do iloop = 1, nloop cfl_ok = .true. do j=js+1,je do i=is,ie if (bm(i,j,l)>=zero) then beta=bm(i,j,l)/mx(i,j-1) else beta=bm(i,j,l)/mx(i,j) end if if (abs(beta)>=one) then cfl_ok = .false. end if end do !i end do !j if ( region == 1.or.touch_sp(region) == 1) then do i=is,ie if (bm(i,2,l)>=zero) then beta=bm(i,2,l)/mx(i,1) else beta=bm(i,2,l)/mx(i,2) end if if (abs(beta)>=one) then cfl_ok = .false. end if end do else !zoom region not touching south pole j = js do i=is,ie !no reduced grid allowed if (bm(i,j,l)>=zero) then beta=bm(i,j,l)/mx(i,j-1) else beta=bm(i,j,l)/mx(i,j) end if if (abs(beta)>=one) then cfl_ok = .false. end if end do end if ! compute boundary fluxes south pole... if ( region == 1.or.touch_np(region) == 1) then ! north pole do i=is,ie if (bm(i,jmr,l)>=zero) then beta=bm(i,jmr,l)/mx(i,jmr-1) else beta=bm(i,jmr,l)/mx(i,jmr) end if if (abs(beta)>=one) then cfl_ok = .false. end if end do else !zoom region not touching north pole j = je+1 do i=is,ie !no reduced grid allowed if (bm(i,j,l)>=zero) then beta=bm(i,j,l)/mx(i,j-1) else beta=bm(i,j,l)/mx(i,j) end if if (abs(beta)>=one) then cfl_ok = .false. end if end do end if ! compute boundary fluxes north pole... if (cfl_ok) then if (region==1) then mx(1:imr,1:jmr) = mx(1:imr,1:jmr) + & bm(1:imr,1:jmr,l) - bm(1:imr,2:jmr+1,l) else if ( touch_sp(region) == 1) then mx(1:imr,1:je) = mx(1:imr,1:je) + & bm(1:imr,1:je,l) - bm(1:imr,2:je+1,l) else if ( touch_np(region) == 1) then mx(1:imr,js:jmr) = mx(1:imr,js:jmr) + & bm(1:imr,js:jmr,l) - bm(1:imr,js+1:jmr+1,l) else mx(is:ie,js:je) = mx(is:ie,js:je) + & bm(is:ie,js :je ,l) - bm(is:ie,js+1:je+1,l) end if end if end do if ( .not. cfl_ok ) then beta = 2.0 bm(:,:,l) = bm(:,:,l)*nloop/(nloop+1) ! reduce mass flux if ( okdebug) print *,'dynamv: ', i,j,l, & 'nloop increased to', nloop + 1, beta nloop = nloop + 1 if(nloop == max_nloop) call escape_tm('nloop too high in dynamv') end if end do !while deallocate(mx) l_nl = max(l_nl, REAL(nloop)) do iloop = 1,nloop ! calculate new air mass distribution if (region==1) then mnew(1:imr,1:jmr) = m(1:imr,1:jmr,l) + & bm(1:imr,1:jmr,l) - bm(1:imr,2:jmr+1,l) else if ( touch_sp(region) == 1) then mnew(1:imr,1:je) = m(1:imr,1:je,l) + & bm(1:imr,1:je,l) - bm(1:imr,2:je+1,l) else if ( touch_np(region) == 1) then mnew(1:imr,js:jmr) = m(1:imr,js:jmr,l) + & bm(1:imr,js:jmr,l) - bm(1:imr,js+1:jmr+1,l) else mnew(is:ie,js:je) = m(is:ie,js:je,l) + & bm(is:ie,js :je ,l) - bm(is:ie,js+1:je+1,l) end if do n=1,ntracetloc offsetn=sum(ntracet_ar(0:myid-1)) nglob = offsetn + n do j=js+1,je do i=is,ie if (bm(i,j,l)>=zero) then beta=bm(i,j,l)/m(i,j-1,l) f(i,j)=beta*(rm(i,j-1,l,n)+(one-beta)*rym(i,j-1,l,n)) pf(i,j)=bm(i,j,l)*(beta*beta*rym(i,j-1,l,n)-3.*f(i,j)) fx(i,j)=beta*rxm(i,j-1,l,n) fz(i,j)=beta*rzm(i,j-1,l,n) else beta=bm(i,j,l)/m(i,j,l) f(i,j)=beta*(rm(i,j,l,n)-(one+beta)*rym(i,j,l,n)) pf(i,j)=bm(i,j,l)*(beta*beta*rym(i,j,l,n)-3.*f(i,j)) fx(i,j)=beta*rxm(i,j,l,n) fz(i,j)=beta*rzm(i,j,l,n) end if l_xi=max(l_xi,abs(beta)) if (abs(beta)>=one) then l_nxi=l_nxi+1 end if end do end do !compute boundary fluxes if ( region == 1.or.touch_sp(region) == 1) then do i=is,ie fz(i,1)=zero fx(i,1)=zero pf(i,1)=zero f(i,1)=zero if (bm(i,2,l)>=zero) then beta=bm(i,2,l)/m(i,1,l) f(i,2)=beta*rm(i,1,l,n) pf(i,2)=-3.*bm(i,2,l)*f(i,2) fx(i,2)=zero fz(i,2)=beta*rzm(i,1,l,n) else beta=bm(i,2,l)/m(i,2,l) f(i,2)=beta*(rm(i,2,l,n)-(one+beta)*rym(i,2,l,n)) pf(i,2)=bm(i,2,l)*(beta*beta*rym(i,2,l,n)-3.*f(i,2)) fx(i,2)=beta*rxm(i,2,l,n) fz(i,2)=beta*rzm(i,2,l,n) end if l_xi=max(l_xi,abs(beta)) if (abs(beta)>=one) then l_nxi=l_nxi+1 end if end do else !zoom region not touching south pole j = js do i=is,ie !no reduced grid allowed if (bm(i,j,l)>=zero) then beta=bm(i,j,l)/m(i,j-1,l) f(i,j)=beta*(rm(i,j-1,l,n)+(one-beta)*rym(i,j-1,l,n)) pf(i,j)=bm(i,j,l)*(beta*beta*rym(i,j-1,l,n)-3.*f(i,j)) fx(i,j)=beta*rxm(i,j-1,l,n) fz(i,j)=beta*rzm(i,j-1,l,n) else beta=bm(i,j,l)/m(i,j,l) f(i,j)=beta*(rm(i,j,l,n)-(one+beta)*rym(i,j,l,n)) pf(i,j)=bm(i,j,l)*(beta*beta*rym(i,j,l,n)-3.*f(i,j)) fx(i,j)=beta*rxm(i,j,l,n) fz(i,j)=beta*rzm(i,j,l,n) end if l_xi=max(l_xi,abs(beta)) if (abs(beta)>=one) then l_nxi=l_nxi+1 end if end do end if ! compute boundary fluxes south pole... if ( region == 1.or.touch_np(region) == 1) then ! north pole do i=is,ie if (bm(i,jmr,l)>=zero) then beta=bm(i,jmr,l)/m(i,jmr-1,l) f(i,jmr)=beta*(rm(i,jmr-1,l,n)+(one-beta)*rym(i,jmr-1,l,n)) pf(i,jmr)=bm(i,jmr,l)*(beta*beta*rym(i,jmr-1,l,n) - & 3.*f(i,jmr)) fx(i,jmr)=beta*rxm(i,jmr-1,l,n) fz(i,jmr)=beta*rzm(i,jmr-1,l,n) else beta=bm(i,jmr,l)/m(i,jmr,l) ! for solid-body test and polar mixing f(i,jmr)=beta*rm(i,jmr,l,n) pf(i,jmr)=-3.*bm(i,jmr,l)*f(i,jmr) fx(i,jmr)=zero ! for solid-body test and polar mixing fz(i,jmr)=beta*rzm(i,jmr,l,n) end if l_xi=max(l_xi,abs(beta)) if (abs(beta)>=one) then l_nxi=l_nxi+1 end if end do else !zoom region not touching north pole j = je+1 do i=is,ie !no reduced grid allowed if (bm(i,j,l)>=zero) then beta=bm(i,j,l)/m(i,j-1,l) f(i,j)=beta*(rm(i,j-1,l,n)+(one-beta)*rym(i,j-1,l,n)) pf(i,j)=bm(i,j,l)*(beta*beta*rym(i,j-1,l,n)-3.*f(i,j)) fx(i,j)=beta*rxm(i,j-1,l,n) fz(i,j)=beta*rzm(i,j-1,l,n) else beta=bm(i,j,l)/m(i,j,l) f(i,j)=beta*(rm(i,j,l,n)-(one+beta)*rym(i,j,l,n)) pf(i,j)=bm(i,j,l)*(beta*beta*rym(i,j,l,n)-3.*f(i,j)) fx(i,j)=beta*rxm(i,j,l,n) fz(i,j)=beta*rzm(i,j,l,n) end if l_xi=max(l_xi,abs(beta)) if (abs(beta)>=one) then l_nxi=l_nxi+1 end if end do end if ! compute boundary fluxes north pole... #ifdef with_budgets !cmk finished computing fluxes...now apply... ! first compute INCOMING fluxes from the south... if (apply_budget_global_flux) then do i=is,ie budget_flux(region)%flux_y1(i,l,nglob) = budget_flux(region)%flux_y1(i,l,nglob) + & f(i,jflux1(region))*1e3/ra(nglob) budget_flux(region)%flux_y2(i,l,nglob) = budget_flux(region)%flux_y2(i,l,nglob) + & f(i,jflux2(region))*1e3/ra(nglob) enddo endif #endif do j=js,je do i=is,ie rm(i,j,l,n) =rm(i,j,l,n)+(f(i,j)-f(i,j+1)) rym(i,j,l,n)=rym(i,j,l,n)+(pf(i,j)-pf(i,j+1) & - (bm(i,j,l)-bm(i,j+1,l))*rym(i,j,l,n) & + 3.*((bm(i,j,l)+bm(i,j+1,l))*rm(i,j,l,n) & - (f(i,j)+f(i,j+1))*m(i,j,l)))/mnew(i,j) if ( limits) rym(i,j,l,n) = max( min(rym(i,j,l,n), & rm(i,j,l,n)), -rm(i,j,l,n) ) rxm(i,j,l,n)=rxm(i,j,l,n)+(fx(i,j)-fx(i,j+1)) rzm(i,j,l,n)=rzm(i,j,l,n)+(fz(i,j)-fz(i,j+1)) end do end do !forall ! poles were excluded..... if ( region==1 .or. touch_sp(region) == 1 ) then rm(1:imr, 1,l,n) = rm(1:imr, 1,l,n) - f(1:imr,2) rzm(1:imr,1 ,l,n) = rzm(1:imr, 1,l,n) - fz(1:imr,2) end if if ( region==1 .or. touch_np(region) == 1 ) then rm(1:imr,jmr,l,n) = rm(1:imr,jmr,l,n) + f(1:imr,jmr) rzm(1:imr,jmr,l,n) = rzm(1:imr,jmr,l,n) + fz(1:imr,jmr) end if end do ! & n-loop ! store new air mass in m array if ( region == 1) then m(1:imr,1:jmr,l) = mnew(1:imr,1:jmr) else if ( touch_sp(region) == 1) then m(1:imr,1:je,l) = mnew(1:imr,1:je) else if ( touch_np(region) == 1) then m(1:imr,js:jmr,l) = mnew(1:imr,js:jmr) else m(is:ie,js:je,l) = mnew(is:ie,js:je) end if ! end do ! iloop bm(:,:,l) = bm(:,:,l)*nloop ! restore old bm's end do! end of l-loop over vertical layers.... call my_omp_max (l_nl, g_nl) call my_omp_max (l_xi, g_xi) call my_omp_sum (l_nxi, g_nxi) deallocate (mnew) deallocate (f) deallocate (pf) deallocate (fx) deallocate (fz) !$OMP END PARALLEL nloop_max(region,2) = max (nloop_max(region,2), nint (g_nl)) xi(region,2) = max (xi(region,2), g_xi) nxi(region,2) = nxi(region,2) + nint (g_nxi) nullify(m) nullify(rm) nullify(rxm) nullify(rym) nullify(rzm) nullify(bm) call uncompress_yedges(region,is,ie,ls,le) call mix_edges(region) end subroutine dynamv subroutine compress_yedges(region,is,ie,ls,le) ! ! this is for slope only: condense data at the y-edges of the zoom region ! to allow for uniform work in dynamv ! written by mike botchev, march-june 1999 ! modified by MK, dec 2002 ! use dims, only : yref, parent, jm, okdebug, touch_sp, touch_np use tracer_data, only : mass_dat use meteodata , only : m_dat use partools , only : ntracetloc use chem_param, only : ntracet ! input/output integer,intent(in) :: region integer,intent(in) :: is integer,intent(in) :: ie integer,intent(in) :: ls integer,intent(in) :: le ! local real,dimension(:,:,:,:),pointer :: rm,rxm,rym,rzm real,dimension(:,:,:), pointer :: m integer :: jmr,i,l,n,yref_ ! start yref_ = yref(region)/yref(parent(region)) jmr = jm(region) if ((yref_==1).or.(region==1)) then if ( okdebug) print *,'compress_yedges: no refinement, nothing to do' return end if if ( okdebug) then print *,'compress_yedges: region=',region,' jmr=',jmr print *,'compress_yedges: cells to be updated: ', & yref_-1,yref_,jmr-yref_+1,jmr-yref_+2 if ( touch_sp(region) == 1) & print *,'compress_yedges: SP will be skipped here' if ( touch_np(region) == 1) & print *,'compress_yedges: NP will be skipped here' end if m => m_dat(region)%data rm => mass_dat(region)%rm_t rxm => mass_dat(region)%rxm_t rym => mass_dat(region)%rym_t rzm => mass_dat(region)%rzm_t do l=ls,le do i=is,ie if ( touch_sp(region) /= 1) then m (i,yref_, l) = sum(m(i,1:yref_,l)) m (i,yref_-1,l) = m (i,0,l) end if if ( touch_np(region) /= 1) then m (i,jmr-yref_+1,l) = sum(m(i,jmr-yref_+1:jmr,l)) m (i,jmr-yref_+2,l) = m (i,jmr+1,l) end if do n=1,ntracetloc if ( touch_sp(region) /= 1) then rm (i,yref_, l,n) = sum(rm(i,1:yref_,l,n)) rm (i,yref_-1,l,n) = rm (i,0,l,n) rxm(i,yref_-1,l,n) = rxm(i,0,l,n) rym(i,yref_-1,l,n) = rym(i,0,l,n) rzm(i,yref_-1,l,n) = rzm(i,0,l,n) end if if ( touch_np(region) /= 1) then rm (i,jmr-yref_+1,l,n) = sum(rm(i,jmr-yref_+1:jmr,l,n)) rm (i,jmr-yref_+2,l,n) = rm (i,jmr+1,l,n) rxm(i,jmr-yref_+2,l,n) = rxm(i,jmr+1,l,n) rym(i,jmr-yref_+2,l,n) = rym(i,jmr+1,l,n) rzm(i,jmr-yref_+2,l,n) = rzm(i,jmr+1,l,n) end if end do !forall end do end do !forall nullify(m) nullify(rm) nullify(rxm) nullify(rym) nullify(rzm) end subroutine compress_yedges subroutine uncompress_yedges(region,is,ie,ls,le) ! ! distribute data at the y-edges of the zoom region ! over the whole interface cell ! written by mike botchev, march-june 1999 ! use dims, only : yref, parent, jm, okdebug, touch_sp, touch_np use tracer_data, only : mass_dat use meteodata , only : m_dat use partools , only : ntracetloc use chem_param, only : ntracet ! input/output integer,intent(in) :: region integer,intent(in) :: is integer,intent(in) :: ie integer,intent(in) :: ls integer,intent(in) :: le ! local real,dimension(:,:,:,:),pointer :: rm,rxm,rym,rzm real,dimension(:,:,:), pointer :: m integer :: jmr,i,l,n,yref_ real :: m_edge,rm_edge ! start yref_ = yref(region)/yref(parent(region)) jmr = jm(region) if ((yref_==1).or.(region==1)) then if ( okdebug) print *,'uncompress_yedges: no refinement, nothnig to do' return end if if ( okdebug) then print *,'uncompress_yedges: region=',region,' jmr=',jmr print *,'uncompress_yedges: cells to be updated: ',1,':', & yref_,' ',jmr-yref_+1,':',jmr if ( touch_sp(region) == 1) & print *,'uncompress_yedges: SP will be skipped here' if ( touch_np(region) == 1) & print *,'uncompress_yedges: NP will be skipped here' end if m => m_dat(region)%data rm => mass_dat(region)%rm_t rxm => mass_dat(region)%rxm_t rym => mass_dat(region)%rym_t rzm => mass_dat(region)%rzm_t do l=ls,le do i=is,ie if ( touch_sp(region) /= 1) then m_edge = m (i,yref_,l) m (i,1:yref_,l) = m_edge/yref_ end if if ( touch_np(region) /= 1) then m_edge = m (i,jmr-yref_+1,l) m (i,jmr-yref_+1:jmr,l) = m_edge/yref_ end if !forall(n=1:ntracet) do n=1,ntracetloc if ( touch_sp(region) /= 1 ) then rm_edge = rm (i,yref_,l,n) rm (i,1:yref_,l,n) = rm_edge/yref_ !* rxm(i,1:yref_-1,l,n) = rxm(i,yref_,l,n) !* rym(i,1:yref_-1,l,n) = rym(i,yref_,l,n) !* rzm(i,1:yref_-1,l,n) = rzm(i,yref_,l,n) rxm(i,1:yref_,l,n) = rxm(i,yref_,l,n)/yref_ rym(i,1:yref_,l,n) = rym(i,yref_,l,n)/yref_ rzm(i,1:yref_,l,n) = rzm(i,yref_,l,n)/yref_ end if if ( touch_np(region) /= 1 ) then rm_edge = rm (i,jmr-yref_+1,l,n) rm (i,jmr-yref_+1:jmr,l,n) = rm_edge/yref_ !* rxm(i,jmr-yref_+2:jmr,l,n) = rxm(i,jmr-yref_+1,l,n) !* rym(i,jmr-yref_+2:jmr,l,n) = rym(i,jmr-yref_+1,l,n) !* rzm(i,jmr-yref_+2:jmr,l,n) = rzm(i,jmr-yref_+1,l,n) rxm(i,jmr-yref_+1:jmr,l,n) = rxm(i,jmr-yref_+1,l,n)/yref_ rym(i,jmr-yref_+1:jmr,l,n) = rym(i,jmr-yref_+1,l,n)/yref_ rzm(i,jmr-yref_+1:jmr,l,n) = rzm(i,jmr-yref_+1,l,n)/yref_ end if end do end do end do nullify(m) nullify(rm) nullify(rxm) nullify(rym) nullify(rzm) end subroutine uncompress_yedges subroutine put_yedges(region) ! ! passes values along the y-boundaries to all the children ! written by mike botchev, march-june 1999 ! adapted for cray-run Maarten Krol, march, 2000 ! touch_np, touch_sp implemented...MK, nov, 2000 ! use dims, only : xref, yref, zref, im, jm, lm use dims, only : touch_sp, touch_np, okdebug, children use dims, only : ibeg, jbeg, lbeg, jend, nregions use tracer_data , only : mass_dat use meteodata , only : m_dat #ifdef with_budgets use budget_global, only : sum_advection #endif use chem_param, only : ra, ntracet use partools , only : offsetl, lmloc use partools , only : ntracetloc use partools , only : which_par, previous_par #ifdef MPI use mpi_const #endif use tracer_data , only : Par_Check_Domain implicit none ! input/output integer,intent(in) :: region ! local real,dimension(:,:,:,:),pointer :: rm,rxm,rym,rzm,rmc,rxmc,rymc,rzmc real,dimension(:,:,:), pointer :: m, mc real,dimension(:,:),allocatable :: toc0,toc1,tocm,tocm1 integer :: child, ichild, j, ip, l, lp, imc, lmc, n, i integer :: xref_, yref_, zref_ real :: xzref, xyzref real :: sum_old, sum_new, sum_old_all, sum_new_all integer :: nzne, nzne_v, nznem integer :: communicator, root_id, lglob, nglob, lmr, nt !integer :: mlev0 ! start which_par=previous_par(region) !WP! make sure parents are on same domain call Par_Check_Domain( region, 'c', which_par ) if ( which_par == 'tracer' .and. ntracetloc == 0 ) return if ( which_par == 'levels' .and. lmloc == 0 ) return !WP! rm => mass_dat(region)%rm_t rxm => mass_dat(region)%rxm_t rym => mass_dat(region)%rym_t rzm => mass_dat(region)%rzm_t m => m_dat(region)%data ichild = 0 do while(ichild m_dat(region)%data if ( which_par=='tracer') then rm => mass_dat(region)%rm_t rxm => mass_dat(region)%rxm_t rym => mass_dat(region)%rym_t rzm => mass_dat(region)%rzm_t mc => m_dat(child)%data rmc => mass_dat(child)%rm_t rxmc => mass_dat(child)%rxm_t rymc => mass_dat(child)%rym_t rzmc => mass_dat(child)%rzm_t lmr = lm(region) nt=ntracetloc #ifdef MPI communicator=com_trac !WP! assign com_trac as communicator root_id=root_t #endif else if ( which_par=='levels') then rm => mass_dat(region)%rm_k rxm => mass_dat(region)%rxm_k rym => mass_dat(region)%rym_k rzm => mass_dat(region)%rzm_k mc => m_dat(child)%data rmc => mass_dat(child)%rm_k rxmc => mass_dat(child)%rxm_k rymc => mass_dat(child)%rym_k rzmc => mass_dat(child)%rzm_k lmr = lmloc nt=ntracet #ifdef MPI communicator=com_lev !WP! assign com_lev as communicator root_id=root_k #endif end if ! ~~ mass (all levels on each pe) imc = im(child) lmc = lm(child) allocate(toc0(imc,lmc)) allocate(toc1(imc,lmc)) allocate(tocm(imc,lmc)) allocate(tocm1(imc,lmc)) xzref = 1./(xref_*zref_) xyzref = 1./(xref_*yref_*zref_) ! loop through the cells of y-walls of the child do l=1,lmc lp = lbeg(child) + (l-1)/zref_ do i=1,imc ip = mod(ibeg(child)-1 + (i-1)/xref_,im(region)) + 1 toc0(i,l) = m(ip,jbeg(child)-1,lp)*xzref toc1(i,l) = m(ip,jbeg(child) ,lp)*xyzref tocm(i,l) = m(ip,jend(child) ,lp)*xyzref tocm1(i,l) = m(ip,jend(child)+1,lp)*xzref end do end do if ( touch_sp(child) == 0) mc(1:imc,0 ,1:lmc) = toc0 if ( touch_np(child) == 0) mc(1:imc,jm(child)+1,1:lmc) = tocm1 do j=1,yref_ if ( touch_sp(child) == 0) mc( 1:imc , j, 1:lmc ) = toc1 if ( touch_np(child) == 0) mc( 1:imc , jm(child)+1-j,1:lmc ) = tocm end do nullify(mc) deallocate(toc0) deallocate(toc1) deallocate(tocm) deallocate(tocm1) ! ~~ tracers imc = im(child) lmc = lmr allocate(toc0(imc,lmc)) allocate(toc1(imc,lmc)) allocate(tocm(imc,lmc)) allocate(tocm1(imc,lmc)) sum_old_all = 0.0 sum_new = 0.0 sum_new_all = 0.0 sum_old = 0.0 do n=1,nt do l=1,lmc lp = lbeg(child) + (l-1)/zref_ do i=1,imc ip = mod(ibeg(child)-1 + (i-1)/xref_,im(region)) + 1 toc0(i,l) = rm(ip,jbeg(child)-1,lp,n)*xzref toc1(i,l) = rm(ip,jbeg(child) ,lp,n)*xyzref tocm(i,l) = rm(ip,jend(child) ,lp,n)*xyzref tocm1(i,l) =rm(ip,jend(child)+1,lp,n)*xzref end do end do if ( touch_sp(child) == 0) rmc(1:imc,0 ,1:lmc,n) = toc0 if ( touch_np(child) == 0) rmc(1:imc,jm(child)+1,1:lmc,n) = tocm1 nglob=n #ifdef MPI if ( myid > 0 .and. which_par == 'tracer' ) & nglob=sum(ntracet_ar(0:myid-1))+n #endif do l=1,lmc lglob=l #ifdef MPI if ( myid > 0 .and. which_par == 'levels' ) & lglob=sum(lmar(0:myid-1))+l #endif do i=1,imc if ( touch_sp(child) == 0 .and. n == 1) & sum_new = sum_new + yref_*toc1(i,l) if ( touch_np(child) == 0 .and. n == 1) & sum_new = sum_new + yref_*tocm(i,l) do j=1,yref_ if ( touch_sp(child) == 0) then if ( n == 1) sum_old = sum_old + rmc( i, j, l ,1) rmc(i,j,l,n) = toc1(i,l) end if if ( touch_np(child) == 0 ) then if ( n == 1 ) sum_old = sum_old + & rmc(i,jm(child)+1-j,l,1) rmc(i,jm(child)+1-j,l,n) = tocm(i,l) end if end do !j end do !i end do !l end do !n #ifdef with_budgets #ifdef MPI !WP! only on one when in tracer if ( which_par == 'tracer' .and. myid /= pe_first_tracer ) sum_new=0.0 if ( which_par == 'tracer' .and. myid /= pe_first_tracer ) sum_old=0.0 #endif !WP! When in levels, each sum_new and sum_old carries the full !WP! sum of each level for tracer one. !WP! When in tracer,sum_new and sum_old is zero on all PE's !WP! except PE_first_tracer #ifdef MPI call mpi_allreduce(sum_new,sum_new_all,1,my_real, & mpi_sum,communicator,ierr) call mpi_allreduce(sum_old,sum_old_all,1,my_real, & mpi_sum,communicator,ierr) #else sum_new_all = sum_new sum_old_all = sum_old #endif sum_advection(child) = sum_advection(child) + sum_new_all - sum_old_all #ifdef MPI if ( okdebug .and. myid==root_id ) & print *, 'put_yedges: region, sum_advection ', & child,sum_advection(child) #else if ( okdebug) print *, 'put_yedges: region, sum_advection ', & child,sum_advection(child) #endif #endif nullify(rmc) do n=1,nt do l=1,lmc lp = lbeg(child) + (l-1)/zref_ do i=1,imc ip = mod(ibeg(child)-1 + (i-1)/xref_,im(region)) + 1 toc0(i,l) = rxm(ip,jbeg(child)-1,lp,n)*xzref toc1(i,l) = rxm(ip,jbeg(child) ,lp,n)*xzref !note the xzref... tocm(i,l) = rxm(ip,jend(child) ,lp,n)*xzref tocm1(i,l) =rxm(ip,jend(child)+1,lp,n)*xzref end do end do if ( touch_sp(child) == 0 ) rxmc(1:imc,0 ,1:lmc,n) = toc0 if ( touch_np(child) == 0 ) rxmc(1:imc,jm(child)+1,1:lmc,n) = tocm1 do j=1,yref_ if ( touch_sp(child) == 0 ) & rxmc( 1:imc , j ,1:lmc ,n) = toc1 if ( touch_np(child) == 0 ) & rxmc( 1:imc , jm(child)+1-j,1:lmc ,n) = tocm end do end do !n nullify(rxmc) do n=1,nt do l=1,lmc lp = lbeg(child) + (l-1)/zref_ do i=1,imc ip = mod(ibeg(child)-1 + (i-1)/xref_,im(region)) + 1 toc0(i,l) = rym(ip,jbeg(child)-1,lp,n)*xzref toc1(i,l) = rym(ip,jbeg(child) ,lp,n)*xzref !note the xzref... tocm(i,l) = rym(ip,jend(child) ,lp,n)*xzref tocm1(i,l) =rym(ip,jend(child)+1,lp,n)*xzref end do end do if ( touch_sp(child) == 0) rymc(1:imc,0 ,1:lmc,n) = toc0 if ( touch_np(child) == 0) rymc(1:imc,jm(child)+1,1:lmc,n) = tocm1 do j=1,yref_ if ( touch_sp(child) == 0) & rymc( 1:imc , j ,1:lmc ,n) = toc1 if ( touch_np(child) == 0) & rymc( 1:imc , jm(child)+1-j,1:lmc ,n) = tocm end do end do !n nullify(rymc) do n=1,nt do l=1,lmc lp = lbeg(child) + (l-1)/zref_ do i=1,imc ip = mod(ibeg(child)-1 + (i-1)/xref_,im(region)) + 1 toc0(i,l) = rzm(ip,jbeg(child)-1,lp,n)*xzref toc1(i,l) = rzm(ip,jbeg(child) ,lp,n)*xzref !note the xzref... tocm(i,l) = rzm(ip,jend(child) ,lp,n)*xzref tocm1(i,l) =rzm(ip,jend(child)+1,lp,n)*xzref end do end do if ( touch_sp(child) == 0) rzmc(1:imc,0 ,1:lmc,n) = toc0 if ( touch_np(child) == 0) rzmc(1:imc,jm(child)+1,1:lmc,n) = tocm1 do j=1,yref_ if ( touch_sp(child) == 0) & rzmc( 1:imc , j ,1:lmc ,n) = toc1 if ( touch_np(child) == 0) & rzmc( 1:imc , jm(child)+1-j,1:lmc ,n) = tocm end do end do !n nullify(rzmc) deallocate(toc0) deallocate(toc1) deallocate(tocm) deallocate(tocm1) nullify(rm) nullify(rxm) nullify(rym) nullify(rzm) nullify(m) end do end subroutine put_yedges end module advecty