!### 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 advectz use GO, only : gol, goPr, goErr implicit none ! --- interface private public :: advectzzoom contains subroutine advectzzoom(region) ! ! set parameters for advectz ! written by patrick berkvens and mike botchev, march-june 1999 ! use dims, only : im, jm, lm, xref, yref, zref, tref use dims, only : parent, zoom2D, touch_sp, touch_np, xcyc use dims, only : splitorderzoom, n_operators, status, nsplitsteps use toolbox, only : escape_tm #ifdef MPI use mpi_const, only : my_real, com_trac, ierr, mpi_sum use mpi_comm, only : barrier_t #endif use partools , only : myid, root_t, ntracetloc implicit none !input/output integer,intent(in) :: region ! local integer :: is,ie,js,je,n,q integer :: imr,jmr,lmr,tref_,xref_,yref_,zref_ logical :: z_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 ! write cells along the z-walls of each child ! to rm(0,...),rm(imr+1,...) and to interface cells of the child ! ! this option should be there: zooming in z disabled.... if ( .not. zoom2D ) then call escape_tm('advectzzoom: ERROR zoom2D should be true:'// & ' zooming along z not allowed') !CMKCALL put_zedges(region,rm,rxm,rym,rzm,m, & !CMK rmgl,rxmgl,rymgl,rzmgl,mgl) end if 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) ! determine the scope for advectz: if ( region == 1 ) then xref_ = 0; yref_ = 0 ! to have is/ie and js/je 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 z_encountered=.false. do n=1,n_operators ! now track following steps from q dir=splitorderzoom(region,q*(nsplitsteps/2)*tref_+n) select case(dir) case('x') if ( (.not.z_encountered) .or. (xcyc(region) == 1) ) then is=1 ! x-substep is before z => ie=im(region) ! full i-scope else is=xref_+1 ! z-substep is before x => ie=im(region)-xref_ ! restricted i-scope end if case('y') if (.not.z_encountered) then js=1 ! y-substep is before z => je=jm(region) ! full j-scope else if(touch_sp(region) == 1 ) then js=1 ! sp is southern edge--->full scope. else js=yref_+1 ! x-substep is before y => end if if ( touch_np(region) == 1 ) then je=jm(region) ! np is northern edge--->full scope else je=jm(region)-yref_ ! restricted j-scope end if end if case('z') z_encountered=.true. case ('c') case ('v') case ('s') case default print *,'advectzzoom: strange value in splitorderzoom(',region,',', & q*(nsplitsteps/2)*tref_+n,'): ',q call escape_tm('advectzzoom: Error in advectz') end select end do call advectz_work(region,is,ie,js,je) end subroutine advectzzoom subroutine advectz_work(region,is,ie,js,je) ! ! calls dynamw ! written by mike botchev, march-june 1999 ! modified by MK, dec 2002 ! ! in/out integer,intent(in) :: region integer,intent(in) :: is integer,intent(in) :: ie integer,intent(in) :: js integer,intent(in) :: je ! start #ifndef slopes stop 'advectz: Wrong advection scheme: adv_scheme/=slope' #endif call dynamw(region,is,ie,js,je) end subroutine advectz_work subroutine dynamw(region,is,ie,js,je) !----------------------------------------------------------------------- ! !**** dynamw - vertical tracer transport v 9.1 ! ! programmed by mh mpi HH 23-feb-1995 ! ! purpose ! ------- ! calculate amount of tracer moved in a vertical advection ! substep ! ! interface ! --------- ! call dynamw ! ! 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 ! changed order of loops for increased performance ! mh, 11-jun-1994 ! 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 binas, only : xgamma=>gamma use dims, only : im, jm, lm, xref, yref, zref, tref use dims, only : zbeg, zend, epsz, nloop_max, zero, one use dims, only : okdebug, nregions, parent, limits, xi, nxi use global_data, only : wind_dat, mass_dat use meteodata , only : m_dat use chem_param, only : ntracet, ra use zoom_tools, only : mix_edges use toolbox, only : escape_tm #ifdef with_budgets use budget_global, only : budget_flux, lflux1, lflux2, apply_budget_global_flux #endif use redgridZoom, only : nred, clustsize, jred, imredj #ifdef MPI use mpi_const, only : my_real,com_trac,ierr, mpi_max, pe_first_tracer use mpi_comm, only : barrier_t #endif use partools , only : myid, root_t use partools , only : ntracetloc, ntracet_ar use omp_partools, only: my_omp_domain, my_omp_max implicit none ! input,output integer,intent(in) :: region integer,intent(in) :: is integer,intent(in) :: ie integer,intent(in) :: js integer,intent(in) :: je ! local real,dimension(:,:,:,:),pointer :: rm,rxm,rym,rzm real,dimension(:,:,:),pointer :: m,cm real,dimension(:,:,:),allocatable :: mnew, f, pf, fx, fy, cmold, mold integer :: i,j,l,le,lee,ls,lss,n,imr,jmr,lmr integer :: lrg integer :: tref_,xref_,yref_,zref_ real :: max_one,max_all, summ real :: mxval, gam, gamma integer,parameter :: iter_limit=15 integer :: n_advz, iter, nglob, offsetn, ixs, ixe, redfact ! Omp parameters integer :: ns, ne, jss, jee, jje, jjs ! start imr=im(region) ; jmr=jm(region) ; lmr=lm(region) !CMK next line: dec2003 reported by FD allocate(cmold(0:im(region)+1,0:jm(region)+1,0:lm(region)+1)) allocate(mold(-1:im(region)+2,-1:jm(region)+2,lm(region))) 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 cm => wind_dat(region)%cm_t if ( okdebug ) print *,'dynamw: region=',region,' is,ie,js,je=',is,ie,js,je if ( ( region < 0 ) .or. ( region > nregions ) ) & call escape_tm( 'dynamw: STOP, illegal number of region !!!') ! prepare z-edges to allow for uniform work: !CMK not any more allowed : ! call compress_zedges(region,is,ie,js,je,m,rm,rxm,rym,rzm) ! 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) ! check is,ie,js,je: if ( ( is /= xref_+1 ) .and. ( is /= 1 ) ) & call escape_tm('dynamw: Wrong value for IS in dynamw') if ( ( ie /= imr-xref_ ) .and. ( ie /= imr ) ) & call escape_tm('dynamw: Wrong value for IE in dynamw') if ( ( js /= yref_+1 ) .and. ( js /= 1 ) ) & call escape_tm('dynamw: Wrong value for JS in dynamw') if ( ( je /= jmr-yref_ ) .and. ( je /= jmr ) ) & call escape_tm('dynamw: Wrong value for JE in dynamw') ! compute ls/le -- cells is:ie,js:je:ls:le will be update ls = zref_; le = lmr-zref_+1 if(abs(zbeg(region)-zbeg(1))=zero) then gamma=cm(i,j,l)/m(i,j,l) else gamma=cm(i,j,l)/m(i,j,l+1) end if xi(region,3)=max(xi(region,3),abs(gamma)) if ( abs(gamma) >= one ) then nxi(region,3)=nxi(region,3)+1 gam=gamma exit advz end if end do end do end do ! no CFLs upto now.....possibly in next? !$OMP PARALLEL & !$OMP default (none) & !$OMP shared (region, is, ie, js, je, lmr, lss, lee, m, & !$OMP cm) & !$OMP private (jss, jee) call my_omp_domain (js, je, jss, jee) if (abs(zbeg(region)-zbeg(1))=one) then if(okdebug)then print *,'dynamw: CFL violation: gamma,m,cm=',& gam,m(i,j,l), m(i,j,l+1),cm(i,j,l), & ' in region,i,j,l,n: ',region,i,j,l,n print*,'dynamw: iterations:',n_advz end if n_advz = n_advz + 1 if ( n_advz > iter_limit ) & call escape_tm('dynamw: too many CFL violations!') cycle cfl else exit cfl end if end do cfl!______CFL ! n_advz is now determined: apply this number of iterations: ! calculate new air mass distribution nloop_max(region,3) = max(nloop_max(region,3), n_advz) ! reset m and cm to original values: !$OMP PARALLEL & !$OMP default (none) & #if defined (with_budgets) !$OMP shared (apply_budget_global_flux, budget_flux, lflux1, lflux2,ra) & #endif !$OMP shared (region, limits, n_advz, is, ie, imr, js, je, & !$OMP jmr, le, ls, lmr, lss, lee, m, mold, cm, cmold, rm, & !$OMP rzm, rxm, rym, xi, ntracetloc, ntracet_ar, myid) & !$OMP private (i, j, l, n, iter, jjs, jje, jss, jee, gamma, & !$OMP mnew, f, pf, fx, fy, nglob) call my_omp_domain (0, jmr + 1, jjs, jje) cm(:,jjs:jje,:) = cmold(:,jjs:jje,:) / float (n_advz) call my_omp_domain (-1, jmr + 2, jjs, jje) m(:,jjs:jje,:) = mold(:,jjs:jje,:) !$OMP BARRIER call my_omp_domain (js, je, jss, jee) allocate (mnew(is:ie,jss:jee,lm(region))) allocate (f(is:ie,jss:jee,0:lm(region))) allocate (pf(is:ie,jss:jee,0:lm(region))) allocate (fx(is:ie,jss:jee,0:lm(region))) allocate (fy(is:ie,jss:jee,0:lm(region))) ! do the iteration do iter=1,n_advz if (abs(zbeg(region)-zbeg(1))=zero) then gamma=cm(i,j,l)/m(i,j,l) f(i,j,l)=gamma*(rm(i,j,l,n)+(one-gamma)* & rzm(i,j,l,n) ) pf(i,j,l)=cm(i,j,l)*(gamma*gamma* & rzm(i,j,l,n)-3.*f(i,j,l)) fx(i,j,l)=gamma*rxm(i,j,l,n) fy(i,j,l)=gamma*rym(i,j,l,n) else gamma=cm(i,j,l)/m(i,j,l+1) f(i,j,l)=gamma*(rm(i,j,l+1,n)-(one+gamma)* & rzm(i,j,l+1,n)) pf(i,j,l)=cm(i,j,l)*(gamma*gamma* & rzm(i,j,l+1,n)-3.*f(i,j,l)) fx(i,j,l)=gamma*rxm(i,j,l+1,n) fy(i,j,l)=gamma*rym(i,j,l+1,n) end if end do end do end do !$OMP BARRIER ! calculate new tracer mass, and tracer mass slopes ! update rm, rzm, rxm and rym in interior layers of the column #ifdef with_budgets ! calculate flux flowing in from below: if (apply_budget_global_flux) then nglob = n + sum(ntracet_ar(0:myid-1)) do j=jss,jee do i=is,ie budget_flux(region)%flux_z1(i,j,nglob) = budget_flux(region)%flux_z1(i,j,nglob) + & f(i,j,lflux1(region)-1)*1e3/ra(nglob) ! moles budget_flux(region)%flux_z2(i,j,nglob) = budget_flux(region)%flux_z2(i,j,nglob) + & f(i,j,lflux2(region)-1)*1e3/ra(nglob) ! moles enddo enddo endif #endif do l=lss,lee !2,lmm1 rm(is:ie,jss:jee,l,n) = rm(is:ie,jss:jee,l,n) + & f(is:ie,jss:jee,l-1)-f(is:ie,jss:jee,l) rzm(is:ie,jss:jee,l,n) = rzm(is:ie,jss:jee,l,n) + & ( pf(is:ie,jss:jee,l-1) - pf(is:ie,jss:jee,l) & - (cm(is:ie,jss:jee,l-1)-cm(is:ie,jss:jee,l))* & rzm(is:ie,jss:jee,l,n) & + 3.0*( (cm(is:ie,jss:jee,l-1)+cm(is:ie,jss:jee,l))* & rm(is:ie,jss:jee,l,n) & - (f(is:ie,jss:jee,l-1)+f(is:ie,jss:jee,l))* & m(is:ie,jss:jee,l) ) ) / mnew(is:ie,jss:jee,l) if ( limits ) then !CMK added rzm(is:ie,jss:jee,l,n) = max(min(rzm(is:ie,jss:jee,l,n), & rm(is:ie,jss:jee,l,n)), & -rm(is:ie,jss:jee,l,n)) end if rxm(is:ie,jss:jee,l,n)=rxm(is:ie,jss:jee,l,n)+ & (fx(is:ie,jss:jee,l-1)-fx(is:ie,jss:jee,l)) rym(is:ie,jss:jee,l,n)=rym(is:ie,jss:jee,l,n)+ & (fy(is:ie,jss:jee,l-1)-fy(is:ie,jss:jee,l)) end do if ( abs(zbeg(region)-zbeg(1)) < epsz ) then ! bottom ! update rm, rzm, rxm and rym near the ground surface rm(is:ie,jss:jee,1,n)=(rm(is:ie,jss:jee,1,n)-f(is:ie,jss:jee,1)) rzm(is:ie,jss:jee,1,n)=rzm(is:ie,jss:jee,1,n)+(-pf(is:ie,jss:jee,1) & +cm(is:ie,jss:jee,1)*rzm(is:ie,jss:jee,1,n) & +3.*(cm(is:ie,jss:jee,1)*rm(is:ie,jss:jee,1,n) & -f(is:ie,jss:jee,1)*m(is:ie,jss:jee,1)))/mnew(is:ie,jss:jee,1) if ( limits ) then !CMK added rzm(is:ie,jss:jee,1,n) = max(min(rzm(is:ie,jss:jee,1,n), & rm(is:ie,jss:jee,1,n)), & -rm(is:ie,jss:jee,1,n)) end if rxm(is:ie,jss:jee,1,n)=rxm(is:ie,jss:jee,1,n)-fx(is:ie,jss:jee,1) rym(is:ie,jss:jee,1,n)=rym(is:ie,jss:jee,1,n)-fy(is:ie,jss:jee,1) end if if (abs(zend(region)-zend(1))