!### 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 advectx use GO, only : gol, goPr, goErr implicit none ! --- interface private public :: advectxzoom, put_xedges contains subroutine advectxzoom( region ) ! set parameters for advectx ! 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 #ifdef with_budgets use budget_global, only : sum_advection #endif use global_data , only : mass_dat, wind_dat use toolbox, only : escape_tm #ifdef MPI use mpi_const, only : my_real, com_lev, ierr use mpi_const, only : root_t, mpi_sum, pe_first_tracer use mpi_comm, only : barrier_t #endif use partools , only : myid, lmloc, ntracetloc implicit none ! input/output integer,intent(in) :: region ! local real,dimension(:,:,:), pointer :: am integer :: js,je,ls,le,n,q integer :: imr,jmr,lmr,tref_,xref_,yref_,zref_ real,dimension(:,:),allocatable :: am0,am1 ! for slopes only logical :: x_encountered character(len=1) :: dir real :: sum_old, sum_new real :: sum_old_all, sum_new_all ! start if ( ntracetloc == 0 ) return !WP! only PE's with nonzero lmloc proceed call put_xedges(region) ! write BC to cildren 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(am0(0:jm(region)+1, 0:lmr+1)) allocate(am1(0:jm(region)+1, 0:lmr+1)) am => wind_dat(region)%am_t sum_old = sum(mass_dat(region)%rm_t(1:im(region),1:jm(region),1:lmr,1)) ! determine the scope for advectx: if ( region == 1 ) then yref_ = 0 zref_ = 0 ! to have js/je and ls/le properly computed end if q=status(region)/((nsplitsteps/2)*tref_) ! find q - the place in the splitorderzoom ! corresponding to the begining of the last ! triple x-y-z of the parent x_encountered=.false. do n=1,n_operators ! now track four following steps from q !wp! changed to n_op.. steps dir=splitorderzoom(region,q*(nsplitsteps/2)*tref_+n) select case(dir) case('x') x_encountered=.true. case('y') if ( .not. x_encountered) then js=1 ! y-substep is before x => je=jm(region) ! full j-scope else if( touch_sp(region) == 1 ) then js=1 !CMK special case...touching SP else js=yref_+1 ! x-substep is before y => end if if( touch_np(region) == 1 ) then je=jm(region) ! CMK special case: touching NP else je=jm(region)-yref_ ! restricted j-scope end if end if case('z') if ( .not. x_encountered ) then ls=1 ! z-substep is before x => le=lmr ! full l-scope else ls=zref_+1 ! x-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 *,'advectxzoom: strange value in splitorderzoom(',region,',', & q*(nsplitsteps/2)*tref_+n,'): ',q call escape_tm( 'advectxzoom: error in x-advection ') end select end do ! n=1,n_operators #ifdef slopes #ifndef secmom if ( ( mod(status(region),(nsplitsteps/2)*tref_) >= n_operators ) .and. & ( im(region)/xref(region) < im(1) ) ) then ! IF (1) more than fourth substep, (2) slope scheme, and ! (3) the region is not [0;360] degrees wide THEN ! at this substep no coarse fluxes will be applied ! for slopes only: zero out velocities via the x-edges of the region ! to assure that no fluxes will be applied if ( okdebug ) print *,'advectxzoom: zeroing out outer mass fluxes' am0 = am( xref_-1,:,:) am( xref_-1,:,:) = zero am1 = am(imr-xref_+1,:,:) am(imr-xref_+1,:,:) = zero end if #endif #endif call advectx_work(region,js,je,ls,le) #ifdef slopes #ifndef secmom if ( ( mod(status(region),(nsplitsteps/2)*tref_) >= n_operators ).and. & ( im(region)/xref(region) < im(1) ) ) then ! for slopes only: recreate velocities via the x-edges am( xref_-1,:,:) = am0 am(imr-xref_+1,:,:) = am1 end if #endif #endif sum_new = sum(mass_dat(region)%rm_t(1:im(region),1:jm(region),1:lmr,1)) #ifdef with_budgets #ifdef MPI if( myid == pe_first_tracer ) & sum_advection(region) = sum_advection(region) + sum_new - sum_old if( okdebug .and. myid == root_t ) & print *, 'advectxzoom: 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 *, 'advectxzoom: region, sum_advection ', & region,sum_advection(region),sum_new,sum_old #endif #endif deallocate(am0) deallocate(am1) nullify(am) end subroutine advectxzoom subroutine advectx_work(region,js,je,ls,le) ! ! makes reduced grid pre-/postprocessing and switches between ! dynamu and dynamu1 ! written by mike botchev, march-june 1999 ! use dims, only : im, jm, lm use redgridZoom, only : grid_reduced, nred, uni2red, uni2red_mf, red2uni use global_data, only : wind_dat use global_data, only : mass_dat use meteodata , only : m_dat use toolbox, only : escape_tm use partools , only : lmloc use chem_param, only : ntracet ! input integer,intent(in) :: region integer,intent(in) :: js integer,intent(in) :: je integer,intent(in) :: ls integer,intent(in) :: le ! local real,dimension(:,:,:),pointer :: m real,dimension(:,:,:),pointer :: am real,dimension(:,:,:),allocatable :: m_uni ! for reduced grid... real,dimension(:,:,:),allocatable :: am_uni ! for reduced grid... real :: mxval integer,dimension(3) :: mxloc integer :: imr,jmr,lmr ! start #ifndef slopes call escape_tm( 'advectx: slopes advection not defined') #endif imr = im(region) ; jmr = jm(region) ; lmr = lm(region) allocate(m_uni(-1:im(region)+2,-1:jm(region)+2, lmr)) allocate(am_uni(0:im(region)+1, 0:jm(region)+1, 0:lmr+1)) am => wind_dat(region)%am_t m => m_dat(region)%data ! transform to the reduced grid: ! check for reduced grid in region if ( grid_reduced .and. (nred(region) /= 0) ) then ! save non-reduced m and am in m_uni and am_uni: m_uni = m am_uni = am ! reduce m,rm,rxm,rym,rzm: call uni2red(region) ! reduce am: call uni2red_mf(region) end if call dynamu (region,js,je,ls,le) ! transform from the reduced grid: if ( grid_reduced .and. nred(region) /= 0 ) then ! advection on uniform grid: m_uni(1:imr,1:jmr,1:lmr) = m_uni(1:imr,1:jmr,1:lmr) + & am_uni(0:imr-1,1:jmr,1:lmr) - am_uni(1:imr,1:jmr,1:lmr) ! redistribute rm,rxm,rym,rzm by using m_uni and m call red2uni(region,m_uni) ! recreate rm/rxm/rym/rzm by using equal masses (not m_uni) !call red2uni_em(region) ! recreate am: am = am_uni end if nullify(am) nullify(m) deallocate(am_uni) deallocate(m_uni) end subroutine advectx_work subroutine dynamu(region,js,je,ls,le) !----------------------------------------------------------------------- ! !**** dynamu - east-west tracer transport - v 9.1 ! ! programmed by mh mpi HH 23-feb-1995 ! ! purpose ! ------- ! calculate amount of tracer moved in an east-west advection ! substep ! ! interface ! --------- ! call dynamu ! ! 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 : okdebug, nregions, parent use dims, only : im, jm, lm, xref, yref, zref, tref use dims, only : zero, one, xi, nxi, nloop_max, limits, xcyc use redgridZoom, only : grid_reduced, nred, imred 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, iflux1, iflux2, apply_budget_global_flux #endif #ifdef MPI use mpi_const, only : ntracetloc,myid,root_t,com_trac use mpi_const, only : ierr,my_real,mpi_max,mpi_min, ntracet_ar use mpi_comm, only : barrier_t #endif 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) :: js integer,intent(in) :: je integer,intent(in) :: ls integer,intent(in) :: le ! local real,dimension(:,:,:,:),pointer :: rm,rxm,rym,rzm real,dimension(:,:,:), pointer :: m,am real,dimension(:),allocatable :: f,pf,fy,fz,mnew integer i,ie,is,j,l,n, nglob, offsetn integer imr,jmr,lmr,tref_,xref_,yref_,zref_,ic,nloop, iloop real :: min_one,max_one,min_all,max_all real :: alpha, max_alpha, x, rmold real,dimension(:), allocatable :: mx logical :: cfl_ok integer, parameter :: max_nloop = 50 ! Openmp parameters integer :: ns, ne, jss, jee, iie, iimr real :: l_xi, g_xi, l_nxi, g_nxi, l_nl, g_nl logical :: special_grid ! start if ( okdebug ) print *,'dynamu: region=',region, & ' js,je,ls,le=',js,je,ls,le if ( (region < 0) .or. (region > nregions) ) & call escape_tm( 'dynamu: STOP, illegal number of region !!!') ! prepare x-edges to allow for uniform work: call compress_xedges(region,js,je,ls,le) ! 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) am => wind_dat(region)%am_t 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 ! check js,je,ls,le: if((js/=yref_+1) .and.(js/=1)) stop 'dynamu: Wrong value for JS ' if((je/=jmr-yref_).and.(je/=jmr)) stop 'dynamu: Wrong value for JE ' if((ls/=zref_+1) .and.(ls/=1)) stop 'dynamu: Wrong value for LS ' if((le/=lmr-zref_).and.(le/=lmr)) stop 'dynamu: Wrong value for LE ' ! compute is/ie -- cells is:ie,js:je:ls:le will be updated is = xref_; ie = imr-xref_+1 if (im(region)/xref(region)==im(1)) then ! periodic boundary condition is = 1; ie = imr end if special_grid = grid_reduced .and. nred(region) /= 0 ! if requested limit zonal slopes such that no negative ! tracer masses should occur if ( limits ) then !$OMP PARALLEL & !$OMP default (none) & !$OMP shared (special_grid, region, imred, is, ie, js, je, ls, le, & !$OMP rm, rxm) & !$OMP private (i, j, l, n, ns, ne) call my_omp_domain (1, ubound (rxm, 4), ns, ne) if (special_grid) then do n = ns, ne do l = ls, le do j = js, je do i = is - 1, imred(j,region) + 1 rxm(i,j,l,n) = max (min (rxm(i,j,l,n), & rm(i,j,l,n)), -rm(i,j,l,n)) end do end do end do end do else do n = ns, ne do l = ls, le do j = js, je do i = is - 1, ie + 1 rxm(i,j,l,n) = max (min (rxm(i,j,l,n), & rm(i,j,l,n)), -rm(i,j,l,n)) end do end do end do end do end if !$OMP END PARALLEL end if if ( okdebug ) then max_one=maxval(m(1:imr,1:jmr,1:lmr)) min_one=minval(m(1:imr,1:jmr,1:lmr)) #ifdef MPI call barrier_t call mpi_reduce(max_one,max_all,1,my_real,mpi_max,root_t,com_trac,ierr) call mpi_reduce(min_one,min_all,1,my_real,mpi_min,root_t,com_trac,ierr) if ( myid == root_t ) print *,'dynamu: m_min,m_max: ',min_all,max_all max_one=maxval(rm(1:imr,1:jmr,1:lmr,:)) min_one=minval(rm(1:imr,1:jmr,1:lmr,:)) call mpi_reduce(max_one,max_all,1,my_real,mpi_max,root_t,com_trac,ierr) call mpi_reduce(min_one,min_all,1,my_real,mpi_min,root_t,com_trac,ierr) if ( myid == root_t ) print *,'dynamu: rm_min,rm_max: ',min_all,max_all #else print *,'dynamu: m_min,m_max: ',min_one, max_one max_one=maxval(rm(1:imr,1:jmr,1:lmr,:)) min_one=minval(rm(1:imr,1:jmr,1:lmr,:)) print *,'dynamu: rm_min,rm_max: ',min_one,max_one #endif end if ! loop over vertical layers and latitudes !$OMP PARALLEL & !$OMP default (none) & #if defined (MPI) !$OMP shared (ntracetloc, ntracet_ar, myid) & #endif #ifdef with_budgets !$OMP shared (apply_budget_global_flux, iflux1, iflux2, budget_flux,ra) & #endif !$OMP shared (region, is, ie, imr, js, je, ls, le, okdebug, limits, & !$OMP special_grid, nred, imred, m, am, rm, rxm, rym, rzm, & !$OMP g_nl, g_nxi, g_xi) & !$OMP private (i, j, l, jss, jee, iie, iimr, nloop, alpha, max_alpha, & !$OMP nglob, cfl_ok, mx, f, pf, fy, fz, mnew, l_nl, l_nxi, & !$OMP l_xi) iie = ie iimr = imr l_nl = 0.0 l_xi = 0.0 l_nxi = 0.0 allocate (f(0:im(region))) allocate (pf(0:im(region))) allocate (fy(0:im(region))) allocate (fz(0:im(region))) allocate (mnew(im(region))) call my_omp_domain (js, je, jss, jee) do j = jss, jee do l=ls,le ! reduced grid: less points will be handled if (special_grid) iie = imred(j,region) ! Determine number of loops required to avoid CFLs... nloop = 1 ! default run the loop one time! alpha = 2.0 allocate(mx(is-1:iie+1)) max_alpha = 0.0 !CMK added oct2003: nloop limit.. do while ( abs(alpha) >= one .and. nloop < max_nloop ) mx(is-1:iie+1) = m(is-1:iie+1,j,l) ! copy mass to temp array if(xcyc(region)==1) then mx(0) = mx(iie) mx(iie+1) = mx(1) end if xloop: do iloop = 1, nloop cfl_ok = .true. do i=is-1,iie if (am(i,j,l)>=zero) then alpha=am(i,j,l)/mx(i) else alpha=am(i,j,l)/mx(i+1) end if if((abs(alpha)>=one)) then !PB cfl_ok = .false. ! if ( okdebug ) & ! print *, 'dynamu: ',i,j,l, ' alpha > 1, since :', & ! am(i,j,l),mx(i) exit xloop end if max_alpha = max(max_alpha, abs(alpha)) end do ! update mass ! note: only is:iie updated, is-1, iie+1 are BC mx(is:iie) = mx(is:iie) + am(is-1:iie-1,j,l) - am(is:iie,j,l) if ( xcyc(region) == 1 ) then !except for periodic BC... mx(0) = mx(iie) mx(iie+1) = mx(1) end if end do xloop if(.not.cfl_ok) then ! reduce mass flux am(is-1:iie,j,l) = am(is-1:iie,j,l)*nloop/(nloop+1) if ( okdebug ) & print *, 'dynamu: ',j,l, 'nloop increased to', & nloop + 1, alpha nloop = nloop + 1 max_alpha = 0.0 if (nloop == max_nloop) call escape_tm('dynuamu: nloop too high') end if end do !while alpha>1 deallocate(mx) l_nl = MAX (l_nl, REAL(nloop)) l_xi = MAX (l_xi, max_alpha) do iloop = 1,nloop ! CFL loop ! calculate new air mass distribution mnew(is:iie)=m(is:iie,j,l) + am(is-1:iie-1,j,l)-am(is:iie,j,l) if (im(region)/xref(region)==im(1)) then m (0 ,j,l) = m (iie,j,l) !periodic BCs... m (iie+1,j,l) = m (1, j,l) end if ! loop over tracers #ifdef MPI do n=1,ntracetloc nglob = n + SUM(ntracet_ar(0:myid-1)) #else do n=1,ntracet nglob = n #endif if (im(region)/xref(region)==im(1)) then rm (0, j,l,n) = rm (iie,j,l,n) !periodic BCs rxm(0, j,l,n) = rxm(iie,j,l,n) rym(0, j,l,n) = rym(iie,j,l,n) rzm(0, j,l,n) = rzm(iie,j,l,n) rm (iie+1,j,l,n) = rm (1, j,l,n) rxm(iie+1,j,l,n) = rxm(1, j,l,n) rym(iie+1,j,l,n) = rym(1, j,l,n) rzm(iie+1,j,l,n) = rzm(1, j,l,n) end if ! calculate fluxes for rm,rxm,rym,rzm do i=is-1,iie if (am(i,j,l)>=zero) then alpha=am(i,j,l)/m(i,j,l) f(i)=alpha*(rm(i,j,l,n)+(one-alpha)*rxm(i,j,l,n)) pf(i)=am(i,j,l)*(alpha*alpha*rxm(i,j,l,n) - 3.*f(i)) fy(i)=alpha*rym(i,j,l,n) fz(i)=alpha*rzm(i,j,l,n) else alpha=am(i,j,l)/m(i+1,j,l) f(i)=alpha*(rm(i+1,j,l,n)-(one+alpha)*rxm(i+1,j,l,n)) pf(i)=am(i,j,l)*(alpha*alpha*rxm(i+1,j,l,n) - 3.*f(i)) fy(i)=alpha*rym(i+1,j,l,n) fz(i)=alpha*rzm(i+1,j,l,n) end if !xi(region,1)=max(xi(region,1),abs(alpha)) if(abs(alpha)>=one) then if (okdebug) then print *,'dynamu: CFL violation: alpha,m=',alpha, & m(i,j,l),' in region,i,j,l: ',region,i,j,l end if l_nxi=l_nxi+1.0 end if end do if (im(region)/xref(region)==im(1)) then ! redefine left fluxes using periodicity: ! reduced grid if (special_grid) iimr = imred(j,region) f(0) = f (iimr) pf(0) = pf(iimr) fy(0) = fy(iimr) fz(0) = fz(iimr) end if #ifdef with_budgets ! ! add up flux budget! ! note that this fails if grid is reduced (ie check) ! if ( apply_budget_global_flux ) then if ( (iflux1(region)-1 >= is) .and. (iflux1(region)-1 <= iie) ) then budget_flux(region)%flux_x1(j,l,nglob) = budget_flux(region)%flux_x1(j,l,nglob) & + f(iflux1(region)-1)*1e3/ra(nglob) ! moles endif if ( (iflux2(region)-1 >= is) .and. (iflux2(region)-1 <= iie) ) then budget_flux(region)%flux_x2(j,l,nglob) = budget_flux(region)%flux_x2(j,l,nglob) & + f(iflux2(region)-1)*1e3/ra(nglob) endif end if #endif ! ! calculate new tracer mass, and tracer mass slopes ! do i=is,iie rm(i,j,l,n)=rm(i,j,l,n)+(f(i-1)-f(i)) rxm(i,j,l,n)=rxm(i,j,l,n)+(pf(i-1)-pf(i) & -(am(i-1,j,l)-am(i,j,l))*rxm(i,j,l,n) & +3.*((am(i-1,j,l)+am(i,j,l))*rm(i,j,l,n) & -(f(i-1)+f(i))*m(i,j,l)))/mnew(i) !CMK: apply limits again: might be in nloop! if ( limits ) rxm(i,j,l,n) = & max(min(rxm(i,j,l,n),rm(i,j,l,n)),-rm(i,j,l,n)) rym(i,j,l,n)=rym(i,j,l,n)+(fy(i-1)-fy(i)) rzm(i,j,l,n)=rzm(i,j,l,n)+(fz(i-1)-fz(i)) end do !forall end do ! end of n-loop ! store new air mass in m array m(is:iie,j,l)=mnew(is:iie) end do ! iloop ! restore 'old' am am(is-1:iie,j,l) = am(is-1:iie,j,l)*nloop end do end do ! end of l-loop/j-loop call my_omp_max (l_nl, g_nl) call my_omp_max (l_xi, g_xi) call my_omp_sum (l_nxi, g_nxi) if (im(region)/xref(region)==im(1)) then ! periodic boundary condition do j=jss,jee ! reduced grid if (special_grid) iimr = imred(j,region) rm (0, j,:,:) = rm (iimr,j,:,:) rxm(0, j,:,:) = rxm(iimr,j,:,:) rym(0, j,:,:) = rym(iimr,j,:,:) rzm(0, j,:,:) = rzm(iimr,j,:,:) rm (iimr+1,j,:,:) = rm (1, j,:,:) rxm(iimr+1,j,:,:) = rxm(1, j,:,:) rym(iimr+1,j,:,:) = rym(1, j,:,:) rzm(iimr+1,j,:,:) = rzm(1, j,:,:) m (0 ,j,:) = m (iimr,j,:) m (iimr+1,j,:) = m (1, j,:) end do end if deallocate (f) deallocate (pf) deallocate (fy) deallocate (fz) deallocate (mnew) !$OMP END PARALLEL nloop_max(region,1) = max (nloop_max(region,1), nint (g_nl)) xi(region,1) = max (xi(region,1), g_xi) nxi(region,1) = nxi(region,1) + nint (g_nxi) nullify(am) nullify(m) nullify(rm) nullify(rxm) nullify(rym) nullify(rzm) call uncompress_xedges(region,js,je,ls,le) call mix_edges(region) end subroutine dynamu subroutine compress_xedges(region,js,je,ls,le) ! ! this is for slope only: condense data at the x-edges of the zoom region ! to allow for uniform work in dynamu/dynamu1 ! written by mike botchev, march-june 1999 ! use dims, only : xref, parent, im, okdebug use global_data, only : mass_dat use meteodata , only : m_dat use partools , only : ntracetloc use chem_param, only : ntracet ! input integer,intent(in) :: region integer,intent(in) :: js integer,intent(in) :: je integer,intent(in) :: ls integer,intent(in) :: le ! local real,dimension(:,:,:,:),pointer :: rm,rxm,rym,rzm real,dimension(:,:,:), pointer :: m integer :: imr,j,l,n,xref_ ! start 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 xref_ = xref(region)/xref(parent(region)) imr = im(region) if ( (xref_ == 1 ) .or. ( im(region)/xref(region) == im(1) ) ) then if ( okdebug ) print *,"compress_xedges:", & " no refinement or periodic bc's, nothing to do" return end if if ( okdebug ) then print *,'compress_xedges: region=',region,' imr=',imr print *,'compress_xedges: cells to be updated: ', & xref_-1,xref_,imr-xref_+1,imr-xref_+2 end if !forall(j=js:je,l=ls:le) do j=js,je do l=ls,le !sum_m = sum(m(1:xref_,j,l),1); m (xref_, j,l) = sum_m m (xref_, j,l) = sum(m(1:xref_,j,l)) m (xref_-1,j,l) = m (0,j,l) m (imr-xref_+1,j,l) = sum(m(imr-xref_+1:imr,j,l)) m (imr-xref_+2,j,l) = m (imr+1,j,l) !forall(n=1:ntracet) #ifdef MPI do n=1,ntracetloc #else do n=1,ntracet #endif rm (xref_, j,l,n) = sum(rm(1:xref_,j,l,n)) rm (xref_-1,j,l,n) = rm (0,j,l,n) rxm(xref_-1,j,l,n) = rxm(0,j,l,n) rym(xref_-1,j,l,n) = rym(0,j,l,n) rzm(xref_-1,j,l,n) = rzm(0,j,l,n) rm (imr-xref_+1,j,l,n) = sum(rm(imr-xref_+1:imr,j,l,n)) rm (imr-xref_+2,j,l,n) = rm (imr+1,j,l,n) rxm(imr-xref_+2,j,l,n) = rxm(imr+1,j,l,n) rym(imr-xref_+2,j,l,n) = rym(imr+1,j,l,n) rzm(imr-xref_+2,j,l,n) = rzm(imr+1,j,l,n) end do !forall end do end do !forall nullify(m) nullify(rm) nullify(rxm) nullify(rym) nullify(rzm) end subroutine compress_xedges subroutine uncompress_xedges(region,js,je,ls,le) ! ! distribute data at the x-edges of the zoom region ! over the whole interface cell ! written by mike botchev, march-june 1999 ! use dims, only : xref, parent, im, okdebug use meteodata , only : m_dat use global_data, only : mass_dat use partools , only : ntracetloc use chem_param, only : ntracet implicit none ! input/output integer,intent(in) :: region integer,intent(in) :: js integer,intent(in) :: je integer,intent(in) :: ls integer,intent(in) :: le ! local real,dimension(:,:,:,:),pointer :: rm,rxm,rym,rzm real,dimension(:,:,:), pointer :: m integer :: imr,j,l,n,xref_ real :: m_edge,rm_edge ! start 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 xref_ = xref(region)/xref(parent(region)) imr = im(region) if ( (xref_==1) .or. (im(region)/xref(region)==im(1)) ) then if ( okdebug ) print *,"uncompress_xedges:", & " no refinement or periodic bc's, nothing to do" return end if if ( okdebug ) then print *,'uncompress_xedges: region=',region,' imr=',imr print *,'uncompress_xedges: cells to be updated: ',1,':', & xref_,' ',imr-xref_+1,':',imr end if ! forall(j=js:je,l=ls:le) do j=js,je do l=ls,le m_edge = m (xref_,j,l) m (1:xref_,j,l) = m_edge/xref_ m_edge = m (imr-xref_+1,j,l) m (imr-xref_+1:imr,j,l) = m_edge/xref_ ! ! !forall(n=1:ntracet) #ifdef MPI do n=1,ntracetloc #else do n=1,ntracet #endif rm_edge = rm (xref_,j,l,n) rm (1:xref_,j,l,n) = rm_edge/xref_ ! !* rxm(1:xref_-1,j,l,n) = rxm(xref_,j,l,n) !* rym(1:xref_-1,j,l,n) = rym(xref_,j,l,n) !* rzm(1:xref_-1,j,l,n) = rzm(xref_,j,l,n) rxm(1:xref_,j,l,n) = rxm(xref_,j,l,n)/xref_ rym(1:xref_,j,l,n) = rym(xref_,j,l,n)/xref_ rzm(1:xref_,j,l,n) = rzm(xref_,j,l,n)/xref_ ! rm_edge = rm (imr-xref_+1,j,l,n) rm (imr-xref_+1:imr,j,l,n) = rm_edge/xref_ ! !* rxm(imr-xref_+2:imr,j,l,n) = rxm(imr-xref_+1,j,l,n) !* rym(imr-xref_+2:imr,j,l,n) = rym(imr-xref_+1,j,l,n) !* rzm(imr-xref_+2:imr,j,l,n) = rzm(imr-xref_+1,j,l,n) rxm(imr-xref_+1:imr,j,l,n) = rxm(imr-xref_+1,j,l,n)/xref_ rym(imr-xref_+1:imr,j,l,n) = rym(imr-xref_+1,j,l,n)/xref_ rzm(imr-xref_+1:imr,j,l,n) = rzm(imr-xref_+1,j,l,n)/xref_ end do end do !forall end do !forall nullify(m) nullify(rm) nullify(rxm) nullify(rym) nullify(rzm) end subroutine uncompress_xedges subroutine put_xedges(region) ! ! passes values along the x-boundaries to all the children ! written by mike botchev, march-june 1999 ! structure implemented by MK, dec 2002 ! use dims, only : im, jm, lm, xref, yref, zref use dims, only : ibeg, jbeg, lbeg, iend, nregions use dims, only : children, okdebug use meteodata , only : m_dat use global_data , only : mass_dat #ifdef with_budgets use budget_global, only : sum_advection #endif use chem_param, only : ra, ntracet #ifdef MPI use mpi_const #endif use partools , only : lmloc, offsetl use partools , only : ntracetloc use partools , only : which_par, previous_par use tracer_data , only : Par_Check_Domain implicit none ! in/out integer,intent(in) :: region ! local real,dimension(:,:,:,:),pointer :: rm,rxm,rym,rzm,rmc,rxmc,rymc,rzmc real,dimension(:,:,:), pointer :: m,mc integer :: child,ichild,j,jp,l,lp,jmc,lmc,n,i integer :: xref_,yref_,zref_ real :: yzref,xyzref real :: sum_old,sum_new,sum_old_all,sum_new_all integer :: nzne, nzne_v,nznem real,dimension(:,:),allocatable :: toc0,toc1,tocm,tocm1 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 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) jmc = jm(child) lmc = lm(child) allocate(toc0(jmc,lmc)) allocate(toc1(jmc,lmc)) allocate(tocm(jmc,lmc)) allocate(tocm1(jmc,lmc)) yzref = 1./(yref_*zref_) xyzref = 1./(xref_*yref_*zref_) ! loop through the cells of x-walls of the child do l=1,lmc lp = lbeg(child) + (l-1)/zref_ do j=1,jmc jp = jbeg(child) + (j-1)/yref_ toc0(j,l) = m(ibeg(child)-1,jp,lp)*yzref toc1(j,l) = m(ibeg(child) ,jp,lp)*xyzref tocm(j,l) = m(iend(child) ,jp,lp)*xyzref tocm1(j,l) = m(iend(child)+1,jp,lp)*yzref end do end do ! copy info to the child .... mc( 0,1:jmc,1:lmc) = toc0 mc(im(child)+1,1:jmc,1:lmc) = tocm1 do i=1,xref_ mc(i , 1:jmc , 1:lmc ) = toc1 mc(im(child)+1-i, 1:jmc , 1:lmc ) = tocm end do nullify(mc) deallocate(toc0) deallocate(toc1) deallocate(tocm) deallocate(tocm1) ! ~~ tracers jmc = jm(child) lmc = lmr allocate(toc0(jmc,lmc)) allocate(toc1(jmc,lmc)) allocate(tocm(jmc,lmc)) allocate(tocm1(jmc,lmc)) sum_old = 0.0 sum_old_all = 0.0 sum_new = 0.0 sum_new_all = 0.0 do n=1,nt do l=1,lmc lp = lbeg(child) + (l-1)/zref_ do j=1,jmc jp = jbeg(child) + (j-1)/yref_ toc0(j,l) = rm(ibeg(child)-1,jp,lp,n)*yzref toc1(j,l) = rm(ibeg(child) ,jp,lp,n)*xyzref tocm(j,l) = rm(iend(child) ,jp,lp,n)*xyzref tocm1(j,l) =rm(iend(child)+1,jp,lp,n)*yzref end do end do rmc( 0,1:jmc,1:lmc,n) = toc0 rmc(im(child)+1,1:jmc,1:lmc,n) = tocm1 ! calculate the mass that is associated with the 'implicit' !advection that is done here.... 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 j=1,jmc if (n == 1) sum_new = sum_new + xref_*(toc1(j,l) + tocm(j,l)) do i=1,xref_ if (n == 1) sum_old = sum_old + & rmc(i,j,l,1) + rmc(im(child)-i+1,j,l,1) rmc(i ,j,l,n) = toc1(j,l) rmc(im(child)+1-i ,j,l,n) = tocm(j,l) end do end do end do end do #ifdef MPI !WP! only on one when in tracer if ( which_par == 'tracer' .and. myid /= pe_first_tracer ) sum_new=0.0 !WP! only on one when in tracer 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 sum !WP! 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 with_budgets #ifdef MPI !add all 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_xedges: region, sum_advection ', & child,sum_advection(child) #else if ( okdebug ) print *, 'put_xedges: 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 j=1,jmc jp = jbeg(child) + (j-1)/yref_ toc0(j,l) = rxm(ibeg(child)-1,jp,lp,n)*yzref toc1(j,l) = rxm(ibeg(child) ,jp,lp,n)*yzref !note the yzref... tocm(j,l) = rxm(iend(child) ,jp,lp,n)*yzref tocm1(j,l) =rxm(iend(child)+1,jp,lp,n)*yzref end do end do rxmc( 0,1:jmc,1:lmc,n) = toc0 rxmc(im(child)+1,1:jmc,1:lmc,n) = tocm1 do i=1,xref_ rxmc(i , 1:jmc , 1:lmc ,n) = toc1 rxmc(im(child)+1-i, 1:jmc , 1:lmc ,n) = tocm end do end do nullify(rxmc) do n=1,nt do l=1,lmc lp = lbeg(child) + (l-1)/zref_ do j=1,jmc jp = jbeg(child) + (j-1)/yref_ toc0(j,l) = rym(ibeg(child)-1,jp,lp,n)*yzref toc1(j,l) = rym(ibeg(child) ,jp,lp,n)*yzref !note the yzref... tocm(j,l) = rym(iend(child) ,jp,lp,n)*yzref tocm1(j,l) =rym(iend(child)+1,jp,lp,n)*yzref end do end do rymc( 0,1:jmc,1:lmc,n) = toc0 rymc(im(child)+1,1:jmc,1:lmc,n) = tocm1 do i=1,xref_ rymc(i , 1:jmc , 1:lmc ,n) = toc1 rymc(im(child)+1-i, 1:jmc , 1:lmc ,n) = tocm end do end do nullify(rymc) do n=1,nt do l=1,lmc lp = lbeg(child) + (l-1)/zref_ do j=1,jmc jp = jbeg(child) + (j-1)/yref_ toc0(j,l) = rzm(ibeg(child)-1,jp,lp,n)*yzref toc1(j,l) = rzm(ibeg(child) ,jp,lp,n)*yzref !note the yzref... tocm(j,l) = rzm(iend(child) ,jp,lp,n)*yzref tocm1(j,l) =rzm(iend(child)+1,jp,lp,n)*yzref end do end do rzmc( 0,1:jmc,1:lmc,n) = toc0 rzmc(im(child)+1,1:jmc,1:lmc,n) = tocm1 do i=1,xref_ rzmc(i , 1:jmc , 1:lmc ,n) = toc1 rzmc(im(child)+1-i, 1:jmc , 1:lmc ,n) = tocm end do end do 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_xedges end module advectx