!### 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" !----------------------------------------------------------------------------- ! TM5 ! !----------------------------------------------------------------------------- !BOP ! ! !MODULE: ADVECTM_CFL ! ! !DESCRIPTION: module to check for CFL violation. ! ! Exactly the same operations are performed on m as in the 'real' advection ! routines. However, rm, rxm, rym, rzm are not updated. The 'old' m is saved ! using store_masses and restored with 'restore_masses'. !\\ !\\ ! !INTERFACE: ! MODULE ADVECTM_CFL ! ! !USES: ! use GO, only : gol, goErr, goPr use dims, only : nregions, maxref use dims, only : revert use dims, only : okdebug use tm5_distgrid, only : dgrid, Get_DistGrid, update_halo, update_halo_jband IMPLICIT NONE PRIVATE ! ! !PUBLIC MEMBER FUNCTIONS: ! public :: Init_CFL, Done_CFL public :: Check_CFL public :: Setup_MassFlow ! ! !PRIVATE TYPES: ! type oldmass_data real, pointer :: m(:,:,:) real, dimension(:,:,:), pointer :: am real, dimension(:,:,:), pointer :: bm real, dimension(:,:,:), pointer :: cm end type oldmass_data ! ! !PRIVATE DATA MEMBERS: ! integer, parameter :: max_global_iteration = 32 integer, parameter :: n_operators = 3 ! xyz integer, parameter :: nsplitsteps = 6 ! xyzzyx character, parameter :: splitorder(nsplitsteps) = (/'x','y','z','z','y','x'/) character(len=*), parameter :: mname = 'advectm_cfl' integer :: global_iteration integer :: ndyn_save integer :: regionm_status(nregions) integer :: mloop_max(nregions,3) = 0 ! local CFL counters real :: xim (nregions,3) = 0.0 ! logical :: cfl_ok type(oldmass_data) :: oldmass_dat(nregions) ! alternate mass array (used on row_roots, to deal with reduced grid). With and w/o halo. real, dimension(:,:,:), allocatable, target :: m_alt1 real, dimension(:,:,:), allocatable :: m_alt2 character :: splitorderzoom(nregions,maxref*nsplitsteps) ! ! !REVISION HISTORY: ! 10 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ CONTAINS !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: INIT_CFL ! ! !DESCRIPTION: initialize data for CFL !\\ !\\ ! !INTERFACE: ! SUBROUTINE INIT_CFL ! ! !USES: ! use dims, only : lm, im use global_data, only : wind_dat use meteodata, only : m_dat use partools, only : isRowRoot, npe_lon use redgridZoom, only : nred ! ! !REVISION HISTORY: ! 10 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC integer :: region, imr, lmr, halo integer :: i1, i2, j1, j2 ! --- begin ------------------------------- DO region = 1, nregions CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) lmr = lm(region) halo=2 ALLOCATE (oldmass_dat(region)%m(i1-halo:i2+halo, j1-halo:j2+halo, lmr)) halo = wind_dat(region)%halo ALLOCATE (oldmass_dat(region)%am(i1-halo:i2+halo, j1-halo:j2+halo, 0:lmr+1)) ALLOCATE (oldmass_dat(region)%bm(i1-halo:i2+halo, j1-halo:j2+halo, 0:lmr+1)) ALLOCATE (oldmass_dat(region)%cm(i1-halo:i2+halo, j1-halo:j2+halo, 0:lmr+1)) END DO ! array to work on row_root (only if really needed) if ( nred(1)/=0 .and. (npe_lon/=1)) then if(isRowRoot)then imr = im(1) lmr = lm(1) halo=m_dat(1)%halo allocate(m_alt1( 1-halo:imr+halo, j1:j2, lmr)) allocate(m_alt2( imr, j1:j2, lmr)) else allocate(m_alt2(1,1,1)) endif endif ! set splitorder splitorderzoom = ' ' splitorderzoom(1,1:nsplitsteps) = splitorder END SUBROUTINE INIT_CFL !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: DONE_CFL ! ! !DESCRIPTION: deallocate data for CFL !\\ !\\ ! !INTERFACE: ! SUBROUTINE DONE_CFL ! ! !REVISION HISTORY: ! !EOP !------------------------------------------------------------------------ !BOC integer :: region ! info ... write (gol,'(" ")'); call goPr write (gol,'("CFL info from check_cfl:")'); call goPr write (gol,'(" x: mloop_max, xim",i4,f10.4)') mloop_max(1,1), xim(1,1); call goPr write (gol,'(" y: mloop_max, xim",i4,f10.4)') mloop_max(1,2), xim(1,2); call goPr write (gol,'(" z: mloop_max, xim",i4,f10.4)') mloop_max(1,3), xim(1,3); call goPr ! clear do region=1, nregions deallocate(oldmass_dat(region)%m ) deallocate(oldmass_dat(region)%am) deallocate(oldmass_dat(region)%bm) deallocate(oldmass_dat(region)%cm) end do if (allocated(m_alt1)) deallocate(m_alt1) if (allocated(m_alt2)) deallocate(m_alt2) END SUBROUTINE DONE_CFL !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: CHECK_CFL ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! SUBROUTINE CHECK_CFL( t1, t2, n, status ) ! ! !USES: ! use GO, only : TDate, IncrDate, operator(+) use global_data, only : wind_dat use dims, only : ndyn, nsrce, nconv, nchem use dims, only : nread ! <-- 'ntimestep' from rcfile, read in 'start' use datetime, only : new_valid_timestep ! ! !INPUT PARAMETERS: ! type(TDate), intent(in) :: t1, t2 ! ! !INPUT/OUTPUT PARAMETERS: ! integer, intent(inout) :: n ! number of time steps: n is increased until no cfl ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Check_CFL' type(TDate) :: tr(2) integer :: i, ndyn_old integer :: region real :: fraction ! --- begin ------------------------------- if ( okdebug ) then write (gol,'(" checking CFL : initial timestep of ",i5," sec with ", i3," timesteps")') ndyn, n; call goPr endif ! increase number of time steps until cfl ok or maximum exceeded global_iteration = 0 do ! increase number of time steps: global_iteration = global_iteration + 1 ! check ... if ( global_iteration > max_global_iteration ) then write (gol,'("exceeded maximum number of time steps")'); call goErr write (gol,'("in ",a)') rname; status=1; return; call goErr end if ! save current air masses in 'oldmass_dat' call store_masses( n ) ! loop over number of ndyn timesteps; ! apply advection over complete (large) time interval do i = 1, n ! small time interval tr(1) = t1 + IncrDate(sec=(i-1)*ndyn) tr(2) = t1 + IncrDate(sec= i *ndyn) ! fill am/bm/cm with balanced mass flows, eventually time interpolated call Setup_MassFlow( tr, status ) IF_ERROR_RETURN(status=1) ! first half step of advection: call determine_cfl_iter(1, status) ! xyz ... IF_NOTOK_RETURN(status=1) ! second half step if first is ok: if ( cfl_ok ) then call determine_cfl_iter(1, status) ! ... zyx IF_NOTOK_RETURN(status=1) end if if ( okdebug ) then write (gol,'(" checking CFL : time step ",i2," of ",i2," ; cfl_ok : ",l1)') i, n, cfl_ok; call goPr endif ! cfl not ok ? then leave loop and decrease time step: if (.not.cfl_ok) exit ! flag ... regionm_status(:) = 0 end do ! restore massa at begin of time interval call restore_masses if ( .not. cfl_ok) then ndyn_old = ndyn ! new ndyn should be a denominator of e.g. 3 hour call new_valid_timestep( ndyn, nread, nread ) if ( okdebug ) then write (gol,'(" checking CFL : reducing timestep to ",i5," sec")') ndyn; call goPr endif nsrce = ndyn nconv = ndyn nchem = ndyn n = nint( real(n*ndyn_old)/real(ndyn) ) ! should work! fraction = real(ndyn)/real(ndyn_old) do region = 1, nregions wind_dat(region)%am = wind_dat(region)%am*fraction wind_dat(region)%bm = wind_dat(region)%bm*fraction wind_dat(region)%cm = wind_dat(region)%cm*fraction enddo else exit end if end do ! global_iteration ! ok status = 0 end subroutine Check_CFL !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: SETUP_MASSFLOW ! ! !DESCRIPTION: Fill am/bm/cm : ! ! o interpolate balanced mass fluxes pu/pv in time (if necessary) ! o multiply with dt; store in am/bm/cm !\\ !\\ ! !INTERFACE: ! subroutine Setup_MassFlow( tr, status ) ! ! !USES: ! use GO, only : TDate use meteodata, only : pu_dat, pv_dat use meteo, only : TimeInterpolation use advect, only : dynam0 ! ! !INPUT PARAMETERS: ! type(TDate), intent(in) :: tr(2) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Setup_MassFlow' integer :: n ! --- begin ---------------------------------- ! ! time interpolation ! do n = 1, nregions ! call TimeInterpolation( pu_dat(n), tr, status ) IF_ERROR_RETURN(status=1) ! call TimeInterpolation( pv_dat(n), tr, status ) IF_ERROR_RETURN(status=1) ! end do ! ! fill am/bm/cm ! do n = 1, nregions ! pu/pv to am/bm/cm call dynam0( n, status ) IF_NOTOK_RETURN(status=1) end do status = 0 end subroutine Setup_MassFlow !EOC ! **************************************************************** ! *** ! *** mass data ! *** ! **************************************************************** !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: STORE_MASSES ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! SUBROUTINE STORE_MASSES(n) ! ! !USES: ! use dims, only : nregions, tref, lm use global_data, only : wind_dat use meteodata, only : m_dat ! ! !INPUT PARAMETERS: ! integer, intent(in) :: n ! number of ndym timesteps... ! ! !REVISION HISTORY: ! 10 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC real, dimension(:,:,:), pointer :: m, am, bm, cm integer :: region, lmr, halo integer :: i1, i2, j1, j2 integer :: ntimes ! --- begin -------------------------- do region = 1, nregions CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) lmr = lm(region) oldmass_dat(region)%m = m_dat(region)%data oldmass_dat(region)%am = wind_dat(region)%am oldmass_dat(region)%bm = wind_dat(region)%bm oldmass_dat(region)%cm = wind_dat(region)%cm regionm_status(region) = 0 if (revert == -1) then m => m_dat(region)%data am => wind_dat(region)%am bm => wind_dat(region)%bm cm => wind_dat(region)%cm ! backward advection (two times since xyz..zyx) ntimes = 2*n*tref(region) ! fluxes for ndyn/2 s m(i1:i2,j1:j2,1:lmr) = m( i1:i2 , j1:j2 , 1:lmr ) & - ntimes*am( i1-1:i2-1, j1:j2 , 1:lmr ) & !east + ntimes*am( i1:i2 , j1:j2 , 1:lmr ) & !west - ntimes*bm( i1:i2 , j1:j2 , 1:lmr ) & !south + ntimes*bm( i1:i2 , j1+1:j2+1, 1:lmr ) & !north - ntimes*cm( i1:i2 , j1:j2 , 0:lmr-1) & !lower + ntimes*cm( i1:i2 , j1:j2 , 1:lmr ) !upper nullify(am) nullify(bm) nullify(cm) nullify(m) endif end do END SUBROUTINE STORE_MASSES !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RESTORE_MASSES ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! SUBROUTINE RESTORE_MASSES ! ! !USES: ! use dims, only : nregions use global_data, only : wind_dat use meteodata, only : m_dat ! ! !REVISION HISTORY: ! !EOP !------------------------------------------------------------------------ !BOC integer :: region do region = 1, nregions m_dat(region)%data = oldmass_dat(region)%m wind_dat(region)%am = oldmass_dat(region)%am wind_dat(region)%bm = oldmass_dat(region)%bm wind_dat(region)%cm = oldmass_dat(region)%cm enddo END SUBROUTINE RESTORE_MASSES !EOC ! **************************************************************** ! *** ! *** mass advection ! *** ! **************************************************************** !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: DETERMINE_CFL_ITER ! ! !DESCRIPTION: ! ! Determine a 'global' iteration that is needed to prevent CFL violations. ! Normally the operator splitting sequence is: ! xyz vvzyx ! xyzvvzyx vzyxxyzv ! If a CFL violation will occur anywhere, the global timestep is reduced and ! the whole sequence is called two (or more) times. The advantage is that ! 'masses' are largely restored by advection from neighbouring cells ! !\\ !\\ ! !INTERFACE: ! SUBROUTINE DETERMINE_CFL_ITER( REGION, STATUS ) ! ! !USES: ! use dims, only: parent, tref, revert, ndyn ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/determine_cfl_iter' integer :: i, tref_ integer, dimension(3) :: ll,uu ! set flag to signal cfl violations cfl_ok = .true. ! determine refinement factor with respect to the parent tref_ = tref(region)/tref(parent(region)) do i=1,tref_ call do_steps(region, status) IF_NOTOK_RETURN(status=1) if (.not. cfl_ok) return !do the remaining steps if necessary...! if (mod(regionm_status(region),n_operators) /= 0 ) then call do_steps(region, status) IF_NOTOK_RETURN(status=1) if(.not. cfl_ok) return end if end do END SUBROUTINE DETERMINE_CFL_ITER !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: DO_STEPS ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! subroutine do_steps(region, status) ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/do_steps' integer :: i123, reg, rgi, j character :: tobedone character(len=1) :: next_step, prev_step ! --- begin --- status=0 if ( okdebug ) then write(gol,*) 'do_steps: region ',region call goPr end if do i123=1,n_operators next_step = splitorderzoom(region,regionm_status(region)+1) if ( regionm_status(region) /= 0 ) then prev_step = splitorderzoom(region,regionm_status(region)) else prev_step = ' ' end if if (okdebug) then ! it's only to make work of the code visible -- step to be done ! will be printed with capital letter (X,Y or Z) do reg=1,region tobedone = ' ' if ( reg == region ) tobedone = upper(next_step) write(gol,*) 'do_steps: ',reg,': ', & splitorderzoom(reg,1:regionm_status(reg)),tobedone call goPr enddo endif tobedone = upper(next_step) select case(next_step) case('x') call advectmx case('y') call dynamvm case('z') call dynamwm case default write(gol,*)'do_steps: strange value in splitorderzoom: ', & splitorderzoom(region,regionm_status(region)) ; call goErr write(gol,*)'do_steps: (must be x, y or z)'; call goErr status=1; TRACEBACK return end select if (.not. cfl_ok) return regionm_status(region) = regionm_status(region)+1 ! these statements are needed when ! more processes are involved that change rm ! for instance emissions. ! for m-advection only just leave when ready if ( mod(regionm_status(region),n_operators) == 0 ) then exit ! e.g.after zyx or xyz end if end do END SUBROUTINE DO_STEPS !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !FUNCTION: UPPER ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! character function upper(xyz) ! ! !INPUT PARAMETERS: ! character(1), intent(in) :: xyz ! ! !REVISION HISTORY: ! !EOP !------------------------------------------------------------------------ !BOC if (xyz=='x') then upper = 'X' else if (xyz=='y') then upper = 'Y' else if (xyz=='z') then upper = 'Z' else upper = '_' end if end function upper !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: ADVECTMX ! ! !DESCRIPTION: does reduced grid pre-/postprocessing and calls dynamu ! !\\ !\\ ! !INTERFACE: ! SUBROUTINE ADVECTMX ! ! !USES: ! use dims, only : im, jm, lm use redgridZoom, only : grid_reduced, nred, uni2red_mf, idle_xadv use global_data, only : wind_dat use meteodata, only : m_dat use partools, only : npe_lon, myid, par_broadcast ! ! !REVISION HISTORY: ! written by mike botchev, march-june 1999 ! 30 Jan 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! 20 Dec 2012 - P. Le Sager - fixed for reduced grid ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/advectmx' real, dimension(:,:,:), pointer :: m real, dimension(:,:,:), pointer :: am real, dimension(:,:,:), allocatable :: m_uni ! m backup real, dimension(:,:,:), allocatable :: am_uni ! am backup integer :: region, imr,jmr,lmr, i1, i2, j1, j2, halo, status ! Block bounds region=1 CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) lmr = lm(region) ! Transform to the reduced grid if ( grid_reduced ) then halo = m_dat(region)%halo allocate( m_uni( i1-halo:i2+halo, j1-halo:j2+halo, lmr )) halo = wind_dat(region)%halo allocate(am_uni( i1-halo:i2+halo, j1-halo:j2+halo, 0:lmr+1)) am => wind_dat(region)%am m => m_dat(region)%data ! take care of BC for uniform am, so that am_uni is correct when ! updating m_uni after the call to dynamun if (npe_lon==1) then am(0,j1:j2,:) = am(im(region),j1:j2,:) else call UPDATE_HALO_JBAND(dgrid(region), am, halo, status ) IF_NOTOK_RETURN(status=1) endif ! save non-reduced m and am in m_uni and am_uni: m_uni = m am_uni = am ! reduce m (simplified local routine) call uni2red(region,i1,i2,j1,j2,status) IF_NOTOK_RETURN(status=1) ! reduce am (and update halo of reduced zonal bands) call uni2red_mf(region,i1,i2,j1,j2, status) IF_NOTOK_RETURN(status=1) end if if (.not.idle_xadv) then CALL DYNAMUM (region,j1,j2,1,lmr) endif ! idle cores need to know about cfl_ok if (nred(region)/=0 .and. (npe_lon/=1)) then ! if on a row of many cores that has reduced grid call par_broadcast(cfl_ok, status, row=.true.) endif ! Transform from the reduced grid if ( grid_reduced ) then ! advection on uniform grid m_uni(i1:i2,j1:j2,1:lmr) = m_uni(i1:i2,j1:j2,1:lmr) + & am_uni(i1-1:i2-1,j1:j2,1:lmr) - am_uni(i1:i2,j1:j2,1:lmr) ! refill m with m_uni, but only boxes that have been reduced halo = m_dat(region)%halo call red2uni(region, i1, i2, j1, j2, halo, m_uni, status) IF_NOTOK_RETURN(status=1) ! recreate rm/rxm/rym/rzm by using equal masses (not m_uni) !call red2uni_em(region) ! recreate am: am = am_uni nullify(am) nullify(m) deallocate(am_uni) deallocate(m_uni) endif END SUBROUTINE ADVECTMX !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: UNI2RED ! ! !DESCRIPTION: transforms data (air mass only) from uniform grid to reduced grid ! !\\ !\\ ! !INTERFACE: ! subroutine uni2red(region,i1,i2,j1,j2,status) ! ! !USES: ! use dims, only : im,jm,lm use redgridZoom, only : nred, clustsize, jred, imredj use meteodata, only : m_dat use partools, only : isRowRoot, npe_lon use tm5_distgrid, only : dgrid, gather_j_band ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region,i1,i2,j1,j2 ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! written by mike botchev, march-june 1999 ; modified by Maarten Krol, dec 2002 ! ! 19 Dec 2012 - P. Le Sager - modified for TM5-MP (decomposition along latitudes only) ! ! !REMARKS: ! - simplified version of redgridZoom/uni2red : applied to air mass only, not to the tracers. ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/uni2red' real, dimension(:,:,:), pointer :: m, m_alt integer :: halo,i,ie,is,j,l,lrg,redfact,lmr,imr real :: summ ! start m => m_dat(region)%data lmr=lm(region) imr=im(region) ! Two scenarios if (npe_lon /= 1) then ! ! decomposition along longitudes => work on row_root ! if (nred(region)/=0) then ! some reduce grid on this tile ! gather m (need to remove i-halo) allocate(m_alt( i1:i2, j1:j2, lmr)) m_alt = m(i1:i2,j1:j2,:) CALL GATHER_J_BAND(dgrid(region), m_alt, m_alt2, status, jref=j1, rowcom=.true.) IF_NOTOK_RETURN(status=1) deallocate(m_alt) if (isRowRoot) then m_alt1( 1:imr,:,:) = m_alt2 m_alt1( 0,:,:) = m_alt1( imr,:,:) m_alt1( -1,:,:) = m_alt1(imr-1,:,:) do lrg=1,nred(region) j = jred(lrg,region) redfact=clustsize(j,region) do l=1,lmr do i = 1,imredj(j,region) ! the is:ie array section will be reduced to i is = (i-1)*redfact + 1 ie = i*redfact summ = sum(m_alt1(is:ie,j,l)) m_alt1(i,j,l) = summ end do !i end do !l enddo endif endif else ! ! Reduced grid without longitudinal decomposition ! do lrg=1,nred(region) j = jred(lrg,region) redfact=clustsize(j,region) do l=1,lmr do i = 1,imredj(j,region) ! the is:ie array section will be reduced to i is = (i-1)*redfact + 1 ie = i*redfact summ = sum(m(is:ie,j,l)) m(i,j,l) = summ end do !i end do !l enddo endif nullify(m) status=0 END SUBROUTINE UNI2RED !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RED2UNI ! ! !DESCRIPTION: transforms data from reduced grid back to uniform grid ! !\\ !\\ ! !INTERFACE: ! subroutine red2uni(region, i1, i2, j1, j2, halo, m_uni, status) ! ! !USES: ! use dims, only : im, lm use redgridZoom, only : nred, clustsize, jred, imredj use meteodata, only : m_dat use partools, only : isRowRoot, npe_lon use tm5_distgrid, only : dgrid, scatter_j_band ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region, i1, i2, j1, j2, halo real, intent(in) :: m_uni(i1-halo:,j1-halo:,1:) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! written by mike botchev, march-june 1999 ! 20 Dec 2012 - P. Le Sager - TM5-MP version ! ! !REMARKS: ! - simplified version of redgridZoom/red2uni : applied to air mass only, not to the tracers ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'_MOD_red2uni' real,dimension(:,:,:), pointer :: m, m_ integer :: i, j, l, ie, is, lrg, n, redfact, lmr, imr ! start m => m_dat(region)%data lmr=lm(region) imr=im(region) ! Two scenarios if ((npe_lon /= 1).and.(nred(region)>0)) then ! ! decomposition along longitudes => work on row_root ! ! (1) use m_alt1 to fill m in bands that are not reduced - only if necessary if (nred(region)< (j2-j1+1)) then allocate(m_( i1:i2, j1:j2, lmr)) if (isRowRoot) m_alt2 = m_alt1( 1:imr,:,:) CALL SCATTER_J_BAND(dgrid(region), m_, m_alt2, status, jref=j1, rowcom=.true.) IF_NOTOK_RETURN(status=1) m(i1:i2,j1:j2,:) = m_ deallocate(m_) endif ! (2) bands with reduced grid, just use m_uni if (nred(region)/=0) then do lrg=1,nred(region) j = jred(lrg,region) m(:,j,:)= m_uni(:,j,:) end do endif ! update halo call UPDATE_HALO_JBAND(dgrid(region), m, m_dat(region)%halo, status ) else ! ! Reduced grid without longitudinal decomposition ! do lrg=1,nred(region) j = jred(lrg,region) redfact=clustsize(j,region) m(i1:i2,j,:)= m_uni(i1:i2,j,:) m(0,j,:) = m(imr,j,:) end do endif nullify(m) status=0 end subroutine red2uni !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RED2UNI_EM ! ! !DESCRIPTION: transforms data from reduced grid back to uniform grid ! ! **************************************************************** ! ***** COMMENTED - NOT PORTED TO TM5MP WITH NPE_LON /=1 !! ****** ! **************************************************************** ! ! !\\ !\\ ! !INTERFACE: ! !PLS subroutine red2uni_em(region, i1, i2, j1, j2) !PLS ! !PLS ! !USES: !PLS ! !PLS use dims, only : im,jm,lm !PLS use redgridZoom, only : nred, clustsize, jred, imredj !PLS use meteodata, only : m_dat !PLS ! !PLS ! !INPUT PARAMETERS: !PLS ! !PLS integer, intent(in) :: region, i1, i2, j1, j2 !PLS ! !PLS ! !REVISION HISTORY: !PLS ! written by mike botchev, march-june 1999 !PLS ! 20 Dec 2012 - P. Le Sager - TM5-MP version !PLS ! !PLS ! !REMARKS: !PLS ! - simplified version of redgridZoom/red2uni_em : applied to air mass only, !PLS ! no the tracers !PLS ! !PLS !EOP !PLS !------------------------------------------------------------------------ !PLS !BOC !PLS !PLS ! local !PLS real,dimension(:,:,:), pointer :: m !PLS integer i,ie,is,j,l,lrg,redfact,lmr !PLS real mass !PLS ! start !PLS lmr=lm(region) !PLS m => m_dat(region)%data !PLS do lrg=1,nred(region) !PLS j = jred(lrg,region) !PLS redfact=clustsize(j,region) !PLS do l=1,lmr !PLS do i = imredj(j,region),1,-1 !PLS is = (i-1)*redfact + 1 !PLS ie = i*redfact !PLS mass=m(i,j,l); m(is:ie,j,l)= mass/(ie-is+1) !PLS enddo !PLS m(0,j,l) = m(im(region),j,l) !PLS end do !PLS end do !PLS nullify(m) !PLS end subroutine red2uni_em !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: DYNAMUM ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! subroutine dynamum(region,js,je,ls,le) ! ! !USES: ! use dims, only : parent, im, lm, zero, one use redgridZoom, only : grid_reduced, nred, imredj, mfl_alt1 use global_data, only : wind_dat use meteodata, only : m_dat use partools, only : Par_Reduce, npe_lon, par_broadcast, isRowRoot ! ! !INPUT PARAMETERS: ! integer,intent(in) :: region integer,intent(in) :: js integer,intent(in) :: je integer,intent(in) :: ls integer,intent(in) :: le ! ! !REMARKS: ! - no J-loop anymore to work on lat-lon domain decomposition [we could keep ! the L-loop if faster] !EOP !------------------------------------------------------------------------ !BOC real, dimension(:,:,:), pointer :: m, am real, dimension(:,:,:), allocatable :: mnew, mx integer :: i, is, j, l, n, i1, i2, j1, j2, status integer :: lmr, nloop, iloop, iemax integer, allocatable :: ie(:) real :: max_alpha real :: my_alpha, alpha, alpha2, x, rmold logical :: cfl_okl integer, parameter :: max_nloop = 10 ! ------------ ! Bounds ! ------------ CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2) lmr = lm(region) allocate(ie(j1:j2)) ! ------------ ! Work arrays ! ------------ if ( ( grid_reduced ) .and. (npe_lon /= 1) .and. (nred(region)/=0) ) then m => m_alt1 am => mfl_alt1 is=1 ie=im(region) iemax=maxval(imredj(j1:j2,region)) else m => m_dat(region)%data ! halo=2, 1:lmr am => wind_dat(region)%am ! halo=1; 0:lmr+1 is=i1 ie=i2 iemax=i2 endif allocate(mnew( is:iemax, j1:j2, ls:le)) ! halo=0; ls:le = 1:lmr allocate( mx(is-1:iemax+1, j1-1:j2+1, ls:le)) ! halo=1, temp arr to store initial mass ! ------------ ! Reduced grid & halo ! -------------- if ( ( grid_reduced .and. npe_lon==1 ).or.(grid_reduced .and. npe_lon/=1.and.nred(region)>0)) then ie = imredj(j1:j2,region) do j=j1,j2 m( 0, j,:) = m(ie(j),j,:) m(ie(j)+1, j,:) = m( 1,j,:) end do ! halo of am is filled in advectmx (thru uni2red_mf for reduced part) else CALL UPDATE_HALO_JBAND( dgrid(region), am, wind_dat(region)%halo, status) ! IF_NOTOK_RETURN(status=1) CALL UPDATE_HALO_JBAND( dgrid(region), m, m_dat(region)%halo, status) ! IF_NOTOK_RETURN(status=1) end if ! ------------ ! Determine number of loops required to avoid CFLs ! ----------------- nloop = 1 ! default run the loop one time! alpha = 2.0 max_alpha = 0.0 !CMK added oct2003: nloop limit.. do while ( abs(alpha) >= one .and. nloop < max_nloop ) ! copy intial mass to temp array do j=js,je mx(is-1:ie(j)+1,j,:) = m(is-1:ie(j)+1,j,ls:le) end do xloop: do iloop = 1, nloop cfl_okl = .true. lloop: do l=ls,le do j=js,je do i=is-1,ie(j) if (am(i,j,l)>=zero) then my_alpha=am(i,j,l)/mx(i,j,l) else my_alpha=am(i,j,l)/mx(i+1,j,l) end if max_alpha = max(max_alpha, abs(my_alpha)) if((abs(my_alpha)>=one)) then exit lloop end if end do end do end do lloop if (grid_reduced .and. npe_lon/=1) then if (nred(region)==0) then call Par_Reduce(max_alpha, 'MAX', alpha2, status, all=.true., row=.true.) max_alpha=alpha2 end if if(isRowRoot)then call Par_Reduce(max_alpha, 'MAX', alpha, status, all=.true., col=.true.) endif if (nred(region)==0) then call par_broadcast(alpha, status, row=.true.) endif else call Par_Reduce(max_alpha, 'MAX', alpha, status, all=.true.) endif if(status==1) then write(gol,*) "error par_all_reduce";call goPr end if if (alpha>=one) then cfl_okl = .false. exit xloop end if ! update mass for next iteration if (iloop/=nloop) then ! [is:ie] do j=js,je mx(is:ie(j),j,:) = mx(is:ie(j),j,:) + am(is-1:ie(j)-1,j,ls:le) - am(is:ie(j),j,ls:le) end do ! bc if ( ( grid_reduced .and. npe_lon==1 ).or.(grid_reduced .and. npe_lon/=1.and.nred(region)>0)) then do j=js,je mx( 0, j,:) = mx(ie(j),j,ls:le) mx(ie(j)+1, j,:) = mx( 1,j,ls:le) end do else call UPDATE_HALO_JBAND( dgrid(region), mx, 1, status) ! IF_NOTOK_RETURN(status=1) end if end if end do xloop if(.not.cfl_okl) then do j=js,je ! reduce mass flux am(is-1:ie(j),j,ls:le) = am(is-1:ie(j),j,ls:le)*nloop/(nloop+1) !FASTER: am(is-1:ie,:,:) = am(is-1:ie,:,:)*(nloop/(nloop+1)) end do if ( okdebug ) then write(gol,*) 'dynamu: nloop increased to', nloop + 1; call goPr write(gol,*) 'dynamu: alpha is', alpha ; call goPr end if nloop = nloop + 1 max_alpha = 0.0 if (nloop == max_nloop) then ! not good solution found: set flag and return ! note that am is restored in the calling routines ! and m will be restored also! cfl_ok = .false. deallocate(mx, ie) deallocate(mnew) nullify(am) nullify(m) return endif end if end do !while alpha>1 mloop_max(region,1) = max(mloop_max(region,1),nloop) xim(region,1) = max(xim(region,1),alpha) ! CFL loop ! do iloop = 1,nloop ! ! ! calculate new air mass distribution ! mnew(is:ie,js:je,ls:le) = m(is:ie,js:je,ls:le) + am(is-1:ie-1,js:je,ls:le)-am(is:ie,js:je,ls:le) ! m(is:ie,js:je,ls:le) = mnew(is:ie,js:je,ls:le) ! ! end do ! iloop : SHOULD BE THE SAME AS do j=js,je ! calculate new air mass distribution m(is:ie(j),j,ls:le) = m(is:ie(j),j,ls:le) + nloop*( am(is-1:ie(j)-1,j,ls:le)-am(is:ie(j),j,ls:le) ) ! restore 'old' am am(is-1:ie(j),j,ls:le) = am(is-1:ie(j),j,ls:le)*nloop end do !Done deallocate(mnew, mx, ie) nullify(am) nullify(m) end subroutine dynamum !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: DYNAMVM ! ! !DESCRIPTION: south-north tracer transport: calculate amount of tracer moved in a ! south-north advection substep ! ! method : slopes scheme ! reference : Russel and Lerner, 1979 !\\ !\\ ! !INTERFACE: ! SUBROUTINE DYNAMVM ! ! !USES: ! use dims, only : yref, zref, jm, lm use dims, only : parent, nregions use dims, only : zero, one use global_data, only : wind_dat use meteodata, only : m_dat use chem_param, only : ntracet use partools, only : Par_Reduce ! ! !REVISION HISTORY: ! programmed by mh mpi HH 23-feb-1994 ! 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 ! 30 Jan 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC real, dimension(:,:,:), pointer :: m,bm real, dimension(:,:), allocatable :: mx real, dimension(:,:,:), allocatable :: mnew integer :: i,j,je,js,l,n,iee,iss integer :: lmr integer :: yref_ real :: sfs,sfzs,sfn,sfzn real :: my_beta, beta, max_beta integer :: iloop, nloop logical :: cfl_okl, SouthPole, NorthPole, SouthBorder, NorthBorder integer :: lrg, redfact, ixe, ixs real :: summ integer :: i1, i2, j1, j2, jss, jee, status integer :: is,ie,ls,le integer, parameter :: region=1 integer, parameter :: max_nloop = 1 !--- start --- ! compute refinement factors with respect to the parent yref_ = yref(region)/yref(parent(region)) ! Indices CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2,& hasSouthPole=SouthPole, hasNorthPole=NorthPole, & hasSouthBorder=SouthBorder, hasNorthBorder=NorthBorder ) lmr = lm(region) js = j1 je = j2 is=i1 ie=i2 ls=1 le=lmr ! Work arrays allocate(mnew( i1:i2, j1:j2, ls:le)) ! halo=0; ls:le = 1:lmr m => m_dat(region)%data ! halo=2, 1:lmr bm => wind_dat(region)%bm ! halo=1; 0:lmr+1 call UPDATE_HALO( dgrid(region), m, m_dat(region)%halo, status) call UPDATE_HALO( dgrid(region), bm, wind_dat(region)%halo, status) ! Reduce work domain if border of zoom region, and take care of poles if (NorthBorder) then ! j2=jmr by design je = j2 - yref_ + 1 if (NorthPole) je=j2-1 endif if (SouthBorder) then ! j1=1 by design js = j1 + yref_ - 1 if (SouthPole) js=j1+1 endif ! indices for mass update jss=js jee=je if (SouthPole) jss=js-1 ! =j1=1 if (NorthPole) jee=je+1 ! =j2=jmr ! loop over vertical layers LEVELS: do l=ls,le ! compute all inner fluxes ! f(*,j,*) is flux entering j-th cell from below?cmk nloop = 1 ! determine nloop per layer! beta = 2.0 max_beta = 0.0 ! allocate(mx(is:ie,js-1:je+1)) allocate(mx( i1-1:i2+1, j1-1:j2+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) xloop: do iloop = 1, nloop cfl_okl = .true. uloop:do j=js, je+1 do i=is, ie if (bm(i,j,l)>=zero) then my_beta=bm(i,j,l)/mx(i,j-1) else my_beta=bm(i,j,l)/mx(i,j) end if max_beta = max(max_beta, abs(my_beta)) if (abs(my_beta)>=one) then exit uloop end if end do !i end do uloop !j call Par_Reduce(max_beta, 'MAX', beta, status, all=.true.) if(status==1) then write(gol,*) "error par_all_reduce";call goErr end if xim(region,2) = max(xim(region,2), beta) if (beta>=one) then cfl_okl = .false. exit xloop end if if (iloop/=nloop) then mx(is:ie,jss:jee) = mx(is:ie,jss:jee) + bm(is:ie,jss:jee,l) - bm(is:ie,jss+1:jee+1,l) CALL UPDATE_HALO( dgrid(region), mx, 1, status) end if end do xloop if ( .not. cfl_okl ) then ! allow NO iterations in y-dir cfl_ok = .false. if ( okdebug) print *,'dynamv: ', i,j,l, beta deallocate(mx) deallocate(mnew) nullify(m) nullify(bm) return end if end do !while deallocate(mx) mloop_max(region,2) = max(mloop_max(region,2),nloop) ! calculate new air mass distribution, and store it in m array do iloop = 1,nloop mnew(is:ie,jss:jee,l) = m(is:ie,jss:jee,l) + bm(is:ie,jss:jee,l) - bm(is:ie,jss+1:jee+1,l) m(is:ie,jss:jee,l) = mnew(is:ie,jss:jee,l) end do ! restore old bm's bm(:,:,l) = bm(:,:,l)*nloop end do LEVELS ! end of l-loop over vertical layers.... deallocate(mnew) nullify(m) nullify(bm) END SUBROUTINE DYNAMVM !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: DYNAMWM ! ! !DESCRIPTION: vertical tracer transport: calculate amount of tracer moved in a ! vertical advection substep ! ! method : slopes scheme ! reference : Russel and Lerner, 1979 !\\ !\\ ! !INTERFACE: ! SUBROUTINE DYNAMWM ! ! !USES: ! use dims, only : lm, zref use dims, only : zbeg, zend, epsz, zero, one use dims, only : parent use global_data, only : wind_dat use meteodata, only : m_dat use partools, only : Par_Reduce ! ! !REVISION HISTORY: ! ! - programmed by mh, mpi HH, 1994-1995 ! - zoom version written by mike botchev, march-june 1999 ! - 1 Feb 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC real, dimension(:,:,:), pointer :: m, cm real, dimension(:,:,:), allocatable :: mnew, cmold, mold real, dimension(:,:), allocatable :: gamma integer :: i,j,l,le,lee,ls,lss,n,lmr,i1,i2,j1,j2 integer :: zref_, halo real :: summ real :: mxval, gam, beta, max_beta integer :: n_advz, iter, lrg, ixs, ixe, redfact, status integer :: is,ie,js,je integer, parameter :: iter_limit=1 ! no iterations! integer, parameter :: region=1 !----------- START ! Indices & work arrays CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2,J_STRT=j1, J_STOP=j2) lmr=lm(region) allocate( mnew( i1-1:i2+1, j1-1:j2+1, lmr) ) ! halo=1 allocate(gamma( i1-1:i2+1, j1-1:j2+1 ) ) ! halo=1 halo = wind_dat(region)%halo allocate(cmold( i1-halo:i2+halo, j1-halo:j2+halo, 0:lmr+1) ) ! same halo as wind_dat) halo = m_dat(region)%halo allocate( mold( i1-halo:i2+halo, j1-halo:j2+halo, lmr ) ) ! same halo as m_dat) m => m_dat(region)%data cm => wind_dat(region)%cm is=i1 ; ie=i2 ; js=j1 ; je=j2 ! compute refinement factors with respect to the parent zref_ = zref(region)/zref(parent(region)) ! 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(i,j)=cm(i,j,l)/m(i,j,l) else gamma(i,j)=cm(i,j,l)/m(i,j,l+1) end if max_beta = max(max_beta, abs(gamma(i,j))) if ( abs(gamma(i,j)) >= one ) then if ( okdebug ) then write(gol,*) 'dynamwm: CFL violation'; call goErr write(gol,*) ' gamma, m, cm=', max_beta, m(i,j,l), m(i,j,l+1), cm(i,j,l), & ' in region,i,j,l: ',region,i,j,l; call goErr endif exit cloop end if end do end do end do cloop call Par_Reduce(max_beta, 'MAX', beta, status, all=.true.) xim(region,3)=max(xim(region,3), beta) if ( beta >= one ) exit advz ! do not allow iterations in Z ! no CFLs upto now.....possibly in next? m(is:ie,js:je,ls:le)=mnew(is:ie,js:je,ls:le) enddo advz ! if still CFL violation increase iteration counter and re-start cfl loop ! escape is called when the number of iterations becomes too large if ( beta >= one ) then cfl_ok = .false. deallocate(mnew) deallocate(gamma) deallocate(cmold) deallocate(mold) nullify(m) nullify(cm) return end if ! reset m and cm to original values: cm = cmold/float(n_advz) m = mold ! do the iteration (Could get rid of mnew) mloop_max(region,3) = max(mloop_max(region,3), n_advz) do iter=1,n_advz if (abs(zbeg(region)-zbeg(1))