!### 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 redgridZoom use GO , only : gol, goErr, goPr use binas, only : pi use dims, only : nregions, im,jm implicit none private !public :: uni2red_mfs public :: uni2red_mf, initredgrid, calc_pdiff public :: uni2red, red2uni, red2uni_em #ifdef secmom public :: uni2red_2nd, red2uni_em_2nd #endif public :: nred, nredmax, clustsize, grid_reduced public :: jred, imred, imredj ! --- const -------------------------------------- character(len=*), parameter :: mname = 'redgridZoom' ! public parameters integer,parameter :: nredmax=60 ! variable filled in 'ModelIntegration'; ! set to a value here for safety ... logical :: grid_reduced = .true. ! the number of reduced latitude circles per region integer,dimension(nregions) :: nred ! number of joined cells per latitude circle and region integer,dimension(nredmax,nregions) :: clustsize ! reduced dimension...... integer,dimension(nredmax,nregions) :: imredj ! reduced dimension over the full jm integer,dimension(180 ,nregions) :: imred ! latitude numbers where the reduction applies... integer,dimension(nredmax,nregions) :: jred ! private parameters real,parameter :: dl = 2.0*pi/im(1) real,parameter :: dp = pi/jm(1) real,parameter :: epsx = epsilon(0.0) logical,parameter :: fromfile = .true. contains ! Fortran code for Split-rg (reduced grid) scheme for advection on the ! globe. ! ! VERSION: 1.1 ! DATE: November 12, 1998 ! ! Version 1.0 Dated October 19, 1998. ! Version 1.1 Corrected red2uni routine and added new uni2red_m ! routine. ! (Orography effects were erroneously not taken into ! account in the reduced to uniform grid conversion) ! ! This code solves the discrete advection equation for the tracer ! mass, with the wind field specified in air mass fluxes: ! ! dtracm/dt = ( ! mu chi (i-half) - mu chi (i+half) ! mv chi (j-half) - mv chi (j+half) ! mw chi (k-half) - mw chi (k+half) ) ! ! where ! ! tracm = tracer mass ! mu, mv, mw = air mass fluxes ! chi = tracm / mass = tracer mixing ratio ! mass = air mass ! ! The integer values of indices i, j, and k correspond to cell centers. ! ! See: ! ! Petersen, A.C., E.J. Spee, H. van Dop, and W. Hundsdorfer, ! An evaluation and intercomparison of four new advection schemes ! for use in global chemistry models, ! Journal of Geophysical Research, 103, 19253-19269, 1998. ! ! Written by Edwin Spee, CWI, Amsterdam, The Netherlands ! ! ! Parts of the code are written by ! Joke Blom and Willem Hundsdorfer (CWI) and ! ! The advection routine is called in the following way: ! ! call advect_split(tracm, mu, mv, mw, mass, dt) ! ! tracm = tracer mass real(im,jm,lm,ntracet) [kg] ! ! mu = longitudinal mass flux real(im,jm,lm) [kg/s] ! (defined on eastward side of grid ! cells; for all grid cells) ! ! mv = latitudinal mass flux real(im,jm,lm) [kg/s] ! (defined on southward side of grid ! cells; only j=2,..,jm are used) ! ! mw = vertical massflux real(im,jm,lm) [kg/s] ! (defined on upper side of grid cell; ! only l=1,..,lm-1 are used since it is ! assumed that there is no flux across ! the top of the model) ! ! The number of latitudinal elements actually used for mu, mv, and mw ! depends on the advection grid; for mv this number is the determined ! by the latitude row with the largest number of cells bordering a ! certain cell-bounding latitude circle. ! ! mass = air mass real(im,jm,lm) [kg] ! ! dt = timestep real [s] ! ! The physical units were chosen for compatibility with the TM3 model ! but can be changed for other global tracer models, provided that ! the calculation of the fluxes is done with an equivalent of tracer ! mixing ratios and not of tracer concentrations (make for other tracer ! models, if necessary, the appropriate changes in the routines ! split_lab, split_phi, and split_eta, where now tracer mixing ratios ! are calculated as tracm / mass; also the grid transformations as done ! in advect_split are different for tracer mass and concentration). ! ! The longitude is counted from west to east, the latitude from south ! to north, and the height from surface to top. ! ! The variables tracm, mu, mv, mw and mass should be declared in ! the calling program under any name; dt is a parameter. The ! parameters im, jm, and lm denote the number of grid cells in ! longitudinal, latitudinal, and vertical direction, respectively. ! ntracet is the number of transported tracers. The letters 'u', 'v', ! and 'w' for the three dimensions are in the code often replaced by ! 'l', 'p', and 'e', denoting the coordinates lambda, phi, and eta, ! respectively (see Petersen et al. 1998). The advection scheme works ! in principle with any grid and coordinate system of the main model. ! However, in the current implementation (made for the TM3 model) the ! variables tracm and mass are assumed to be defined in the main ! program on a uniform grid with polar cap, and the variables mu, mv, ! and mw on the advection grid (i.e., reduced grid with polar cap). ! Locally in the advection routines a transformation of tracm and mass, ! if necessary, is made to the advection grid. Please contact us for ! information on the changes that have to be made to work with other ! main model grids. Depending on the interests of users, we intend to ! generalize the code with respect to the main model grid in the next ! version. ! ! The following parameters and variables, related to the specification ! of the advection grid are used by the routines of the advection ! scheme: ! ! integer, parameter :: nredmax = 10 ! real, parameter :: dl = 2.0 * pi / im ! real, parameter :: dp = pi / jm ! integer nred ! integer clustsize(0:nredmax), ! jred(-nredmax:nredmax+1),imredj(jm) ! real maxtau ! integer advsub ! ! nredmax = maximum number of grid reductions per hemisphere ! (in this implementation the minimum number of grid ! reductions per hemisphere equals one, since use is ! made of polar caps) ! dl = longitudinal uniform grid cell size in lambda ! coordinates ! dp = latitudinal grid cell size in phi coordinates ! nred = number of grid reductions ! clustsize = array with the ratio im / imred, as a function of the ! absolute value of the latitude zone index, where imred ! is the actually used number of cells in the longitudinal ! direction for the given latitude zone ! jred = latitude number where grid reduction starts (seen from ! SP), as a function of the reduction zone index ! imredj(j) = number of cells for latitude #j ! maxtau = time step such that CFL <= 1 in all directions (given ! the directional splitting, see routine advect_split) ! advsub = largest integer such that dtadv / advsub <= maxtau, ! with dtadv the advection time step in operator splitting ! ! The advection grid variables for the reduced grid are automatically ! set by the following line: ! ! call initredgrid(clustsize, jred, imredj, nred, nredmax, im, jm, ! . pi, .false.) ! ! The variables maxtau and advsub should be set each time new mass ! fluxes (wind fields) are available in the model through ! ! call split_maxtau(mu, mv, mw, mass, dtadv) ! ! where dtadv is the advection time step in the main model. ! ! Furthermore, some routines in the advection scheme need information ! on the machine precision. The following variable is used for this ! purpose: ! ! real, parameter :: epsx = epsilon(0.0) ! ! The value of epsx is dependent on the computing platform that is used. ! epsilon is a Fortran 90 intrinsic function. For example, for our Cray ! C90 eps is about 1e-14. ! ! The above-mentioned parameters and variables are assumed to be ! declared in the module global_data (Fortran 90) or in some common ! block (Fortran 77). The current implementation (Fortran 90) uses the ! statements ! ! use global_data ! implicit none ! ! If Fortran 77 is preferred one should define a common block or ! use an existing common block to include the above-given parameters ! and variables, e.g. ! ! implicit none ! (declarations as given in the above) ! common/advgrid/nredmax,dl,dp,nred,clustsize,jred,maxtau,advsub ! ! Also the parameters im, jm, lm, ntracet, pi, twopi (= 2 * pi), and ! eps are assumed to be declared in the module global_data. ! ! -------------------------- ! Example of a reduced grid ! -------------------------- ! ! NP=North Pole ! SP=South Pole ! EQ=Equator ! | grid boundary in longitudal direction ! o cell center ! j imredj clustsize latitude zone ! NP| o | 13 1 12 3 ! | o | o | o | 12 3 4 2 ! | o | o | o | o | o | o | 11 6 2 1 ! | o | o | o | o | o | o | 10 6 2 1 ! |o|o|o|o|o|o|o|o|o|o|o|o| 9 12 1 0 ! |o|o|o|o|o|o|o|o|o|o|o|o| 8 12 1 0 ! EQ|o|o|o|o|o|o|o|o|o|o|o|o| 7 12 1 0 ! |o|o|o|o|o|o|o|o|o|o|o|o| 6 12 1 0 ! |o|o|o|o|o|o|o|o|o|o|o|o| 5 12 1 0 ! | o | o | o | o | o | o | 4 6 2 -1 ! | o | o | o | o | o | o | 3 6 2 -1 ! | o | o | o | 2 3 4 -2 ! SP| o | 1 1 12 -3 ! ! This grid would have ! im = 12, jm = 13, nred = 3 ! ! jred(-3) = 1 by definition ! jred(-2) = 2 ! jred(-1) = 3 ! jred( 0) = 5 ! jred( 1) =10 ! jred( 2) =12 ! jred( 3) =13 ! jred( 4) =14=jm+1 by definition ! ! clustsize(3) = 12 ! clustsize(2) = 4 ! clustsize(1) = 2 ! clustsize(0) = 1 !CMK total new definition to allow zoom regios.. ! nred(nregions) :: number of reduced latitudes circles < jm(region) ! nredmax :: maximum number of reduced latitude circles... ! clustsize(nredmax.nregions) :: number of joined cells for n'th circle ! jred(nredmax.nregions) :: corresponding j value in 1...jm(region) array ! imredj(nredmax.nregions) :: 'new' im value for the reduced grid... !CMK ! subroutine uni2red_mfs(region,afl,bfl,cfl) ! ! !converts am,bm, and cm to reduced grid. ! !am is compressed ! !bm is averaged (all mass-fluxes will be the same) ! !cm is added (since we just loop over the reduced grid) ! use dims, only : im, jm, lm ! use dims, only : okdebug ! use dims, only : one !#ifdef MPI ! use mpi_const,only : lmloc !#endif ! implicit none ! ! ! in/out ! integer,intent(in) :: region ! real,dimension(0:im(region)+1,0:jm(region)+1,0:lm(region)) :: afl ! real,dimension(0:im(region)+1,0:jm(region)+1,0:lm(region)) :: bfl ! real,dimension(0:im(region)+1,0:jm(region)+1,0:lm(region)) :: cfl ! ! ! local ! integer :: clust_size, i, j, l, i2, imax, jreg, is, ie ! integer :: lmr ! ! if (okdebug) then ! print *,'uni2red_mfs: region =',region ! end if ! !#ifdef MPI ! lmr=lmloc !#else ! lmr = lm(region) !#endif ! ! do l=1,lmr ! do jreg = 1,nred(region) !loop over the reduced latitudes.... ! j = jred(jreg,region) !latitide ! clust_size = clustsize(jreg,region) !clustersize.... ! do i=1,imredj(jreg,region) ! afl(i,j,l) = afl(i*clust_size,j,l) !compress afl ! is = (i-1)*clust_size + 1 ! ie = i*clust_size ! cfl(i,j,l) = sum(cfl(is:ie,j,l)) !put the total flux ! ! bfl(i,j+1,l) = sum(bfl(is:ie,j+1,l))/(ie-is+1) ! ! calculate the average of norther flux! ! end do ! do i=imredj(jreg,region)+1,im(region) ! afl(i,j,l) = -one ! cfl(i,j,l) = -one ! end do ! do i=imredj(jreg,region),1,-1 ! is = (i-1)*clust_size + 1 ! ie = i*clust_size ! ! bfl(is:ie,j+1,l)=bfl(i,j+1,l) ! ! put the same value everywhere on northern edge.. ! end do ! ! the next statement is only meaningful for region cyclic regions; ! ! this subroutine should be applied only in this region ! ! set eastern boundary condition... ! afl(0,j,l)=afl(imredj(jreg,region),j,l) ! end do ! end do ! ! end subroutine uni2red_mfs subroutine uni2red_mf(region) ! Converts mass fluxes mfl from ! uniform grid to reduced grid ! ! Written by Edwin Spee, CWI, Amsterdam, The Netherlands ! ! On input: ! mfl = mass fluxes on uniform grid ! ! On output: ! mfl = mass fluxes on reduced grid !cmk changed..... use dims, only : lm use dims, only : okdebug #ifdef MPI use mpi_const, only: lmloc #endif use global_data, only: wind_dat implicit none ! in/out integer,intent(in) :: region ! local real,dimension(:,:,:),pointer :: mfl integer :: clust_size, i, j, l, i2, imax, jreg integer :: lmr mfl => wind_dat(region)%am_t lmr=lm(region) if ( okdebug ) then print *,'uni2red_mf: region =',region end if ! selection for mfl do l=1,lmr do jreg = 1,nred(region) !loop over the reduced latitudes.... j = jred(jreg,region) !latitide clust_size = clustsize(jreg,region) !clustersize.... do i=1,imredj(jreg,region) mfl(i,j,l) = mfl(i*clust_size,j,l) !compress mfl end do do i=imredj(jreg,region)+1,im(region) mfl(i,j,l) = -1.0 end do ! the next statement is only meaningful for region cyclic regions; ! this subroutine should be applied only in this region ! set eastern boundary condition... mfl(0,j,l)=mfl(imredj(jreg,region),j,l) end do end do nullify(mfl) end subroutine uni2red_mf !----------------------------------------------------------------------- ! Below are some routines related to the initialization of the advection ! grid data structure (i.e. the reduced grid). !----------------------------------------------------------------------- subroutine initredgrid( region, imr, jmr, status ) ! The subroutine initredgrid fills the arrays clustsize, jred, and ! imredj with appropriate values for the use of a reduced grid. ! ! Written by Edwin Spee, CWI, Amsterdam, The Netherlands ! ! This implementation is based on work of Joke Blom. ! ! On input: ! clustsize,jred,imredj,nred undefined ! nredmax ! im,jm number of grids in lon, lat direction ! pi 3.141.... ! if fromfile=true: read grid definition from file, else generate one. ! ! On output: ! array clustsize,jred,imredj contain grid definition ! integer nred is number of grid reductions use dims , only : okdebug use GO , only : TrcFile, Init, Done, ReadRc use global_data, only : rcfile ! in/out integer,intent(in) :: region integer,intent(in) :: imr integer,intent(in) :: jmr integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/initredgrid' ! --- local ------------------------------ integer :: j, lrg character(len=256) :: fname type(TrcFile) :: rcF character(len=256) :: line ! --- begin ----------------------------- ! open rcfile: call Init( rcF, rcfile, status ) IF_NOTOK_RETURN(status=1) ! info ... write (gol,'(" scan rcfile for `regions` key ...")'); call goPr ! try to read the 'regions' variable; ! does only exist for codes based on the pycasso branch: call ReadRc( rcF, 'regions', line, status, default='' ) IF_ERROR_RETURN(status=1) ! check if the list is empty (or key was not found in rcfile): if ( len_trim(line) == 0 ) then ! info ... write (gol,'(" not found; read reduced grid settings from data file ...")'); call goPr ! read name of reduced grid definition file: call ReadRc( rcF, 'proces.advection.RedGrid', fname, status ) IF_NOTOK_RETURN(status=1) ! read reduced grid definition from data file: call read_redgrid( region, fname, status ) IF_NOTOK_RETURN(status=1) else ! info ... write (gol,'(" found; read reduced grid settings from rcfile ...")'); call goPr ! variable 'regions' found, so read reduced grid from rcfile: call Read_RedGrid_rc( region, status ) IF_NOTOK_RETURN(status=1) end if ! close: call Done( rcF, status ) IF_NOTOK_RETURN(status=1) ! fill: imred(1:jm(region),region) = im(region) ! the full jm array... do lrg=1,nred(region) imredj(lrg,region)=imr/clustsize(lrg,region) imred(jred(lrg,region),region)=imredj(lrg,region) end do if ( okdebug ) then do j=1,nred(region) print *,'initredgrid: imredj(',jred(j,region),')=',imredj(j,region) end do end if ! ok status = 0 end subroutine initredgrid subroutine read_redgrid( region, fname, status ) use GO , only : goGetFU use dims , only : okdebug use dims , only : region_name, children use dims , only : ybeg, yend, dy, yref ! --- in/out ---------------------------------- integer, intent(in) :: region character(len=*), intent(in) :: fname integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/read_redgrid' ! --- local ---------------------------------------- integer :: fu logical :: exist integer :: jmr character(len=32) :: rnam integer :: i,j integer :: nc, child,latps,latpn integer :: jband, ncomb integer, allocatable :: ncombs(:) ! --- begin ---------------------------------------- ! number of lats jmr = jm(region) ! select free file unit: call goGetFu( fu, status ) IF_NOTOK_RETURN(status=1) if ( okdebug ) then print *,'read_redgrid: Reading parameters from file RedGrid.' print *,'read_redgrid: Region '//trim(region_name(region)) end if ! check ... inquire( file=trim(fname), exist=exist ) if ( .not. exist ) then write (gol,'("reduced grid file not found :")'); call goErr write (gol,'(" ",a)') trim(fname); call goErr TRACEBACK; status=1; return end if ! open file: open( fu, file=trim(fname), status='old', iostat=status ) if ( status /= 0 ) then write (gol,'("from open reduced grid file:")'); call goErr write (gol,'(" file : ",a)') trim(fname); call goErr write (gol,'(" iostat : ",i6)') status; call goErr TRACEBACK; status=1; return end if ! search record with region name: do read(fu,'(a)',iostat=status) rnam if ( status < 0 ) then ! end of file ... close(fu) write (gol,'("WARNING - No reduced grid for region :",a)') trim(region_name(region)); call goPr nred(region) = 0 !default no reduced grid... status=0; return else if ( status > 0 ) then ! error ... write (gol,'("reading rnam from reduced grid file:")'); call goErr write (gol,'(" ",a)') trim(fname); call goErr TRACEBACK; status=1; return end if if ( trim(rnam) == trim(region_name(region)) ) exit end do ! read number of reduced lat bands: read (fu,*) nred(region) if ( nred(region) > nredmax ) then write (gol,'("problem with reduced grid. nred > nredmax")'); call goErr TRACEBACK; status=1; return end if ! storage for combined cells for each lat; 1 by default allocate( ncombs(jmr) ) ncombs = 1 ! loop over lat bands in redgrid file: do i = 1, nred(region) ! read latitude band and number of cells to combine: read (fu,*,iostat=status) jband, ncomb if ( status /= 0 ) then write (gol,'("from reading reduced grid file:")'); call goErr write (gol,'(" file : ",a)') trim(fname); call goErr write (gol,'(" region : ",a)') trim(rnam); call goErr write (gol,'(" lat band : ",i4)') i; call goErr TRACEBACK; status=1; return end if ! check ... if ( (jband < -1*jmr) .or. (jband > jmr) ) then write (gol,'("strange lat band index in redgrid file:")'); call goErr write (gol,'(" redgrid file : ",a)') trim(fname); call goErr write (gol,'(" region : ",i2," ",a)') region, trim(region_name(region)); call goErr write (gol,'(" lat band : ",i4)') i; call goErr write (gol,'(" combined cells : ",i4)') ncombs(j); call goErr write (gol,'(" im : ",i4)') im(region); call goErr TRACEBACK; status=1; return end if ! store band: if ( jband >= 0 ) then jred(i,region) = jband else jred(i,region) = jband + jm(region) + 1 end if ! clustsize is number of original cells combined into one reduced cell: clustsize(i,region) = ncomb ! idem for testing if ( jband >= 0 ) then ncombs(jmr+1-jband) = ncomb else ncombs(0-jband) = ncomb end if ! now check if the latitude circle overlaps with ! southermost latitude of gridcel parent ... latps = ybeg(region) + (jred(i,region)-1)*nint(dy)/yref(region) ! northermost latitude ... latpn = ybeg(region) + (jred(i,region))*nint(dy)/yref(region) do nc = 1, children(region,0) child = children(region,nc) if (latps < yend(child) .and. latpn > ybeg(child) ) then write (gol,'("overlap detected between child region and reduced grid parent")'); call goErr write (gol,'(" parent : ",i4)') region; call goErr write (gol,'(" child : ",i4)') child; call goErr TRACEBACK; status=1; return end if end do end do ! lat bands in file ! testing ... do j = 1, jmr ! check if number of combined cells matches with grid: if ( modulo(im(region),ncombs(j)) /= 0 ) then write (gol,'("number of combined cells not ok:")'); call goErr write (gol,'(" redgrid file : ",a)') trim(fname); call goErr write (gol,'(" region : ",i2," ",a)') region, trim(region_name(region)); call goErr write (gol,'(" lat band : ",i4)') j; call goErr write (gol,'(" combined cells : ",i4)') ncombs(j); call goErr write (gol,'(" im : ",i4)') im(region); call goErr TRACEBACK; status=1; return end if ! check with previous ... if ( j > 1 ) then if ( modulo( max(ncombs(j-1),ncombs(j)), min(ncombs(j-1),ncombs(j)) ) /= 0 ) then write (gol,'("number of combined cells does match with previous:")'); call goErr write (gol,'(" redgrid file : ",a)') trim(fname); call goErr write (gol,'(" region : ",i2," ",a)') region, trim(region_name(region)); call goErr write (gol,'(" lat band : ",i4)') j; call goErr write (gol,'(" combined cells : ",i4)') ncombs(j); call goErr write (gol,'(" previous : ",i4)') ncombs(j-1); call goErr TRACEBACK; status=1; return end if end if end do ! lats ! clear deallocate( ncombs ) ! close file: close( fu, iostat=status ) if ( status /= 0 ) then write (gol,'("from close reduced grid file:")'); call goErr write (gol,'(" ",a)') trim(fname); call goErr TRACEBACK; status=1; return end if ! ok status = 0 end subroutine read_redgrid ! * subroutine Read_RedGrid_rc( region, status ) use GO , only : TrcFile, Init, Done, ReadRc use global_data, only : rcfile use dims , only : okdebug use dims , only : region_name, children use dims , only : ybeg, yend, dy, yref ! --- in/out ---------------------------------- integer, intent(in) :: region integer, intent(out) :: status ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/Read_RedGrid_rc' ! --- local ---------------------------------------- type(TrcFile) :: rcF integer :: jmr integer :: i,j integer :: nc, child,latps,latpn integer :: jband, ncomb integer, allocatable :: ncombs(:) integer :: nred_nh integer, allocatable :: ncombs_nh(:) integer :: nred_sh integer, allocatable :: ncombs_sh(:) ! --- begin ---------------------------------------- ! number of lats jmr = jm(region) ! open settings: call Init( rcF, rcfile, status ) IF_NOTOK_RETURN(status=1) ! number of reduced lat bands: call ReadRc( rcF, 'region.'//trim(region_name(region))//'.redgrid.nh.n', nred_nh, status ) IF_NOTOK_RETURN(status=1) ! number of reduced lat bands: call ReadRc( rcF, 'region.'//trim(region_name(region))//'.redgrid.sh.n', nred_sh, status ) IF_NOTOK_RETURN(status=1) ! total number of reduced bands: nred(region) = nred_nh + nred_sh ! reduced grid ? if ( nred(region) > 0 ) then if ( okdebug ) then write (gol,'("read_redgrid: Reading parameters from file rcfile.")'); call goPr write (gol,'("read_redgrid: Region ",a)') trim(region_name(region)); call goPr end if if ( nred(region) > nredmax ) then write (gol,'("problem with reduced grid. nred > nredmax")'); call goErr TRACEBACK; status=1; return end if ! storage for combined cells for each lat; 1 by default allocate( ncombs(jmr) ) ncombs = 1 ! counter of reduced lat bands from south-pole to north-pole: i = 0 ! southern hemisphere ? if ( nred_sh > 0 ) then ! temporary storage: allocate( ncombs_sh(jmr) ) ! number of combined cells per lat band: call ReadRc( rcF, 'region.'//trim(region_name(region))//'.redgrid.sh.comb', ncombs_sh(1:nred_sh), status ) IF_NOTOK_RETURN(status=1) ! loop over lat bands: do jband = 1, nred_sh ! increase counter: i = i + 1 ! select: ncomb = ncombs_sh(jband) ! store band: jred(i,region) = jband ! clustsize is number of original cells combined into one reduced cell: clustsize(i,region) = ncomb ! idem for testing ncombs(jmr+1-jband) = ncomb ! now check if the latitude circle overlaps with ! southernmost latitude of gridcel parent ... latps = ybeg(region) + (jred(i,region)-1)*nint(dy)/yref(region) ! northermost latitude ... latpn = ybeg(region) + (jred(i,region))*nint(dy)/yref(region) do nc = 1, children(region,0) child = children(region,nc) if (latps < yend(child) .and. latpn > ybeg(child) ) then write (gol,'("overlap detected between child region and reduced grid parent")'); call goErr write (gol,'(" parent : ",i4)') region; call goErr write (gol,'(" child : ",i4)') child; call goErr TRACEBACK; status=1; return end if end do end do ! lat bands sh ! clear: deallocate( ncombs_sh ) end if ! sh ! northern hemisphere ? if ( nred_nh > 0 ) then ! temporary storage: allocate( ncombs_nh(jmr) ) ! number of combined cells per lat band: call ReadRc( rcF, 'region.'//trim(region_name(region))//'.redgrid.nh.comb', ncombs_nh(1:nred_nh), status ) IF_NOTOK_RETURN(status=1) ! loop over lat bands; read in order from pole to equator, thus reverse for original storage: do jband = nred_nh, 1, -1 ! increase counter: i = i + 1 ! select: ncomb = ncombs_nh(jband) ! store band: jred(i,region) = jm(region) + 1 - jband ! clustsize is number of original cells combined into one reduced cell: clustsize(i,region) = ncomb ! idem for testing ncombs(jband) = ncomb ! now check if the latitude circle overlaps with ! southermost latitude of gridcel parent ... latps = ybeg(region) + (jred(i,region)-1)*nint(dy)/yref(region) ! northermost latitude ... latpn = ybeg(region) + (jred(i,region))*nint(dy)/yref(region) do nc = 1, children(region,0) child = children(region,nc) if (latps < yend(child) .and. latpn > ybeg(child) ) then write (gol,'("overlap detected between child region and reduced grid parent")'); call goErr write (gol,'(" parent : ",i4)') region; call goErr write (gol,'(" child : ",i4)') child; call goErr TRACEBACK; status=1; return end if end do end do ! lat bands nh ! clear: deallocate( ncombs_nh ) end if ! nh ! testing ... do j = 1, jmr ! check if number of combined cells matches with grid: if ( modulo(im(region),ncombs(j)) /= 0 ) then write (gol,'("number of combined cells not ok:")'); call goErr write (gol,'(" region : ",i2," ",a)') region, trim(region_name(region)); call goErr write (gol,'(" lat band : ",i4)') j; call goErr write (gol,'(" combined cells : ",i4)') ncombs(j); call goErr write (gol,'(" im : ",i4)') im(region); call goErr TRACEBACK; status=1; return end if ! check with previous ... if ( j > 1 ) then if ( modulo( max(ncombs(j-1),ncombs(j)), min(ncombs(j-1),ncombs(j)) ) /= 0 ) then write (gol,'("number of combined cells does match with previous:")'); call goErr write (gol,'(" region : ",i2," ",a)') region, trim(region_name(region)); call goErr write (gol,'(" lat band : ",i4)') j; call goErr write (gol,'(" combined cells : ",i4)') ncombs(j); call goErr write (gol,'(" previous : ",i4)') ncombs(j-1); call goErr TRACEBACK; status=1; return end if end if end do ! lats ! clear deallocate( ncombs ) else write (gol,'("WARNING - No reduced grid for region :",a)') trim(region_name(region)); call goPr nred(region) = 0 !default no reduced grid... end if ! close: call Done( rcF, status ) IF_NOTOK_RETURN(status=1) ! ok status = 0 end subroutine Read_RedGrid_rc ! * subroutine calc_pdiff(region,p,pold,pdiffmax) !---------------------------------- ! !---------------------------------- use dims, only : xref, yref, parent, xcyc, touch_sp, touch_np implicit none ! in/out integer,intent(in) :: region real,dimension(-1:im(region)+2,-1:jm(region)+2),intent(in) :: pold,p real,intent(out) :: pdiffmax ! local integer :: lrg,i,j,ratio,idx,iu,xref_,yref_ real :: work,work1 pdiffmax = 0.0 yref_ = yref(region)/yref(parent(region)) xref_ = xref(region)/xref(parent(region)) if ( xcyc(region) /= 1 .or. nred(region) == 0 ) then ! no reduced grid: jloop: do j=1,jm(region) if (yref_>1) then !there is a refinement in the y direction if( (j < yref_ + 1) .and. & (touch_sp(region) /= 1) ) cycle jloop if( (j > jm(region) - yref_) .and. & (touch_np(region) /= 1) ) cycle jloop end if iloop: do i=1,im(region) if ( xref_ > 1 .and. xcyc(region) /= 1 ) then !there is refinement in the x direction and region not cyclic. if( (i < xref_+1) .or. ( i>im(region)-xref_) ) cycle iloop end if pdiffmax = max(pdiffmax,abs(p(i,j)-pold(i,j))) end do iloop end do jloop else ! gid is reduced: do something special: do lrg = 1,nred(region) !first maximum over the reduced grid...! ratio = clustsize(lrg,region) j = jred(lrg,region) do i = 1,imredj(lrg,region) work = 0.0 work1 = 0.0 do iu = 1,ratio idx = (i-1)*ratio+iu work = work + p(idx,j) work1 = work1 + pold(idx,j) end do pdiffmax=max(pdiffmax,abs(work-work1)/ratio) end do end do ! now the pressure difference in the non-reduced part... yref_ = yref(region)/yref(parent(region)) j2loop: do j=1,jm(region) if(imred(j,region).ne.im(region)) cycle j2loop !if reduced...skip if ( yref_ > 1 ) then ! there is a refinement in the y direction: skip the edges! if( (j < yref_ + 1) .and. & (touch_sp(region) /= 1) ) cycle j2loop if( (j > jm(region) - yref_) .and. & (touch_np(region) /= 1) ) cycle j2loop end if i2loop: do i=1,im(region) pdiffmax = max(pdiffmax,abs(p(i,j)-pold(i,j))) end do i2loop end do j2loop end if end subroutine calc_pdiff subroutine uni2red(region) ! ! transforms data from uniform grid to reduced grid ! written by mike botchev, march-june 1999 ! modified by Maarten Krol, dec 2002 ! use dims, only: im,jm,lm use meteodata , only: m_dat use global_data, only: mass_dat use omp_partools,only: my_omp_domain #ifdef MPI use mpi_const,only : ntracetloc #else use chem_param,only : ntracet #endif implicit none ! input integer,intent(in) :: region ! local real,dimension(:,:,:,:),pointer :: rm #ifdef slopes real,dimension(:,:,:,:),pointer :: rxm,rym,rzm #endif real,dimension(:,:,:), pointer :: m integer i,ie,is,j,l,lrg,redfact,n,lmr,nt real summ,sumrm ! Openmp paramters integer nts, nte ! start #ifdef MPI nt=ntracetloc #else nt=ntracet #endif m => m_dat(region)%data rm => mass_dat(region)%rm_t #ifdef slopes rxm => mass_dat(region)%rxm_t rym => mass_dat(region)%rym_t rzm => mass_dat(region)%rzm_t #endif lmr=lm(region) do lrg=1,nred(region) redfact=clustsize(lrg,region) j = jred(lrg,region) do l=1,lmr do i = 1,imredj(lrg,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 !cmkm_uni(is:ie,j,l) = m_uni(is:ie,j,l)/summ ! use as distribution function ! when transferring back from reduced--->uniform grid... ! these summations mean that mixing ratio and the ! the slopes are averaged out within the is:ie section ! with m(is:ie,...) taken as the weights call my_omp_domain(1,nt,nts,nte) do n=nts,nte do i = 1,imredj(lrg,region) is = (i-1)*redfact + 1 ie = i*redfact sumrm = sum(rm(is:ie,j,l,n)) rm(i,j,l,n) = sumrm #ifdef slopes sumrm = sum(rxm(is:ie,j,l,n)) rxm(i,j,l,n) = sumrm sumrm = sum(rym(is:ie,j,l,n)) rym(i,j,l,n) = sumrm sumrm = sum(rzm(is:ie,j,l,n)) rzm(i,j,l,n) = sumrm #endif end do ! JFM: set remaining masses to zero ! for consistency with adjoint do i = imredj(lrg,region)+1, im(region) rm(i,j,l,n) = 0. #ifdef slopes rxm(i,j,l,n) = 0. rym(i,j,l,n) = 0. rzm(i,j,l,n) = 0. #endif end do end do !put periodic boundary... end do !l end do !redgrid... nullify(m) nullify(rm) #ifdef slopes nullify(rxm) nullify(rym) nullify(rzm) #endif end subroutine uni2red #ifdef secmom subroutine uni2red_2nd(region) ! ! transforms data from uniform grid to reduced grid ! written by mike botchev, march-june 1999 ! modified by Maarten Krol, dec 2002 ! use dims use meteodata , only: m_dat use global_data, only: mass_dat #ifdef MPI use mpi_const,only : ntracetloc #else use chem_param,only : ntracet #endif implicit none ! input integer,intent(in) :: region ! local real,dimension(:,:,:,:),pointer :: rm, rxm,rym,rzm, rxxm, rxym, rxzm, ryym, ryzm, rzzm real,dimension(:,:,:), pointer :: m integer i,ie,is,j,l,lrg,redfact,n,lmr real summ,sumrm ! 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 rxxm => mass_dat(region)%rxxm_t rxym => mass_dat(region)%rxym_t rxzm => mass_dat(region)%rxzm_t ryym => mass_dat(region)%ryym_t ryzm => mass_dat(region)%ryzm_t rzzm => mass_dat(region)%rzzm_t lmr=lm(region) do lrg=1,nred(region) redfact=clustsize(lrg,region) j = jred(lrg,region) do l=1,lmr do i = 1,imredj(lrg,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 !cmkm_uni(is:ie,j,l) = m_uni(is:ie,j,l)/summ ! use as distribution function ! when transferring back from reduced--->uniform grid... ! these summations mean that mixing ratio and the ! the slopes are averaged out within the is:ie section ! with m(is:ie,...) taken as the weights #ifdef MPI do n=1,ntracetloc #else do n=1,ntracet #endif sumrm = sum(rm(is:ie,j,l,n)) rm(i,j,l,n) = sumrm sumrm = sum(rxm(is:ie,j,l,n)) rxm(i,j,l,n) = sumrm sumrm = sum(rym(is:ie,j,l,n)) rym(i,j,l,n) = sumrm sumrm = sum(rzm(is:ie,j,l,n)) rzm(i,j,l,n) = sumrm sumrm = sum(rxxm(is:ie,j,l,n)) rxxm(i,j,l,n) = sumrm sumrm = sum(rxym(is:ie,j,l,n)) rxym(i,j,l,n) = sumrm sumrm = sum(rxzm(is:ie,j,l,n)) rxzm(i,j,l,n) = sumrm sumrm = sum(ryym(is:ie,j,l,n)) ryym(i,j,l,n) = sumrm sumrm = sum(ryzm(is:ie,j,l,n)) ryzm(i,j,l,n) = sumrm sumrm = sum(rzzm(is:ie,j,l,n)) rzzm(i,j,l,n) = sumrm end do !n end do !i !put periodic boundary... end do !l end do !redgrid... nullify(m) nullify(rm) nullify(rxm) nullify(rym) nullify(rzm) nullify(rxxm) nullify(rxym) nullify(rxzm) nullify(ryym) nullify(ryzm) nullify(rzzm) end subroutine uni2red_2nd #endif subroutine red2uni(region,m_uni) ! ! transforms data from reduced grid back to uniform grid ! written by mike botchev, march-june 1999 ! use dims use global_data, only : mass_dat use meteodata, only : m_dat use omp_partools,only : my_omp_domain #ifdef MPI use mpi_const, only : ntracetloc #else use chem_param, only : ntracet #endif implicit none ! input integer,intent(in) :: region real,dimension(-1:im(region)+2,-1:jm(region)+2,lm(region)), & intent(in) :: m_uni ! m_uni: same declaration as m! ! local real,dimension(:,:,:,:),pointer :: rm #ifdef slopes real,dimension(:,:,:,:),pointer :: rxm,rym,rzm #endif real,dimension(:,:,:) ,pointer :: m integer :: i, ie, ii, is, j, l, lrg, n, redfact, lmr, nt real :: hi, mass, mass_coord, rmm, slope, m_old character(len=5) :: distr_mode ! Openmp paramters integer :: nts, nte ! start #ifdef MPI nt=ntracetloc #else nt=ntracet #endif lmr=lm(region) m => m_dat(region)%data rm => mass_dat(region)%rm_t #ifdef slopes rxm => mass_dat(region)%rxm_t rym => mass_dat(region)%rym_t rzm => mass_dat(region)%rzm_t #endif !vh : remove if-statement in loop! ! distr_mode = 'unfrm' ! 'slope' or 'unfrm' ! if (distr_mode=='unfrm') then ! ! mixing ratio and x-slope will be UNiFoRMly distributed ! ! within is:ie do lrg=1,nred(region) redfact=clustsize(lrg,region) j = jred(lrg,region) do l=1,lmr do i = imredj(lrg,region),1,-1 ! the i cell will be distributed within the is:ie array section is = (i-1)*redfact + 1 ie = i*redfact !m_uni is the mass-distribution in the non-reduced grid/divided by !the reduced_grid mass. This is used as distribution function!.... mass=m(i,j,l) m(is:ie,j,l)= m_uni(is:ie,j,l) ! if (distr_mode=='unfrm') then ! mixing ratio and x-slope will be UNiFoRMly distributed ! within is:ie do n=1,nt rmm = rm(i,j,l,n) rm(is:ie,j,l,n)= m_uni(is:ie,j,l)/mass* rmm #ifdef slopes rmm = rxm(i,j,l,n) rxm(is:ie,j,l,n)= m_uni(is:ie,j,l)/mass * rmm ! rym and rzm are always distributed uniformly: rmm = rym(i,j,l,n) rym(is:ie,j,l,n)= m_uni(is:ie,j,l)/mass * rmm rmm = rzm(i,j,l,n) rzm(is:ie,j,l,n)= m_uni(is:ie,j,l)/mass * rmm #endif end do ! end if !cmkelseif(distr_mode=='slope') then end do ! update cell(0,...) according to the periodic bc's: ! update cell(0,...) according to the periodic bc's: rm(0,j,l,1:nt) = rm(im(region),j,l,1:nt) #ifdef slopes rxm(0,j,l,:) = rxm(im(region),j,l,:) rym(0,j,l,:) = rym(im(region),j,l,:) rzm(0,j,l,:) = rzm(im(region),j,l,:) #endif m(0,j,l) = m(im(region),j,l) end do end do nullify(m) nullify(rm) #ifdef slopes nullify(rxm) nullify(rym) nullify(rzm) #endif end subroutine red2uni subroutine red2uni_em(region) ! ! transforms data from reduced grid back to uniform grid ! written by mike botchev, march-june 1999 ! use dims use meteodata , only: m_dat use global_data, only: mass_dat #ifdef MPI use mpi_const, only : ntracetloc #else use chem_param,only : ntracet #endif implicit none ! input integer,intent(in) :: region ! local real,dimension(:,:,:,:),pointer :: rm #ifdef slopes real,dimension(:,:,:,:),pointer :: rxm,rym,rzm #endif real,dimension(:,:,:) ,pointer :: m integer i,ie,ii,is,j,l,lrg,n,redfact,lmr real hi,mass,mass_coord,rmm,slope,m_old character(len=5) distr_mode ! start lmr=lm(region) m => m_dat(region)%data rm => mass_dat(region)%rm_t #ifdef slopes rxm => mass_dat(region)%rxm_t rym => mass_dat(region)%rym_t rzm => mass_dat(region)%rzm_t #endif distr_mode = 'unfrm' ! 'slope' or 'unfrm' do lrg=1,nred(region) redfact=clustsize(lrg,region) j = jred(lrg,region) do l=1,lmr do i = imredj(lrg,region),1,-1 ! the i cell will be distributed within the is:ie array section is = (i-1)*redfact + 1 ie = i*redfact !m_uni is the mass-distribution in the non-reduced grid/divided by !the reduced_grid mass. This is used as distribution function!.... mass=m(i,j,l); m(is:ie,j,l)= mass/(ie-is+1) if (distr_mode=='unfrm') then ! mixing ratio and x-slope will be UNiFoRMly distributed ! within is:ie #ifdef MPI do n=1,ntracetloc #else do n=1,ntracet #endif rmm = rm(i,j,l,n) rm(is:ie,j,l,n)= m(is:ie,j,l)/mass* rmm #ifdef slopes rmm = rxm(i,j,l,n) rxm(is:ie,j,l,n)= m(is:ie,j,l)/mass * rmm ! rym and rzm are always distributed uniformly: rmm = rym(i,j,l,n) rym(is:ie,j,l,n)= m(is:ie,j,l)/mass * rmm rmm = rzm(i,j,l,n) rzm(is:ie,j,l,n)= m(is:ie,j,l)/mass * rmm #endif end do end if !cmkelseif(distr_mode=='slope') then end do ! update cell(0,...) according to the periodic bc's: #ifdef MPI rm(0,j,l,1:ntracetloc) = rm(im(region),j,l,1:ntracetloc) #else rm(0,j,l,1:ntracet) = rm(im(region),j,l,1:ntracet) #endif #ifdef slopes rxm(0,j,l,:) = rxm(im(region),j,l,:) rym(0,j,l,:) = rym(im(region),j,l,:) rzm(0,j,l,:) = rzm(im(region),j,l,:) #endif m(0,j,l) = m(im(region),j,l) end do end do nullify(m) nullify(rm) #ifdef slopes nullify(rxm) nullify(rym) nullify(rzm) #endif end subroutine red2uni_em #ifdef secmom subroutine red2uni_em_2nd(region) ! ! transforms data from reduced grid back to uniform grid ! written by mike botchev, march-june 1999 ! use dims use meteodata, only: m_dat use global_data, only: mass_dat #ifdef MPI use mpi_const, only : ntracetloc #else use chem_param,only : ntracet #endif implicit none ! input integer,intent(in) :: region ! local real,dimension(:,:,:,:),pointer :: rm,rxm,rym,rzm, rxxm, rxym, rxzm, ryym, ryzm, rzzm real,dimension(:,:,:) ,pointer :: m integer i,ie,ii,is,j,l,lrg,n,redfact,lmr real hi,mass,mass_coord,rmm,slope,m_old character(len=5) distr_mode ! start lmr=lm(region) 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 rxxm => mass_dat(region)%rxxm_t rxym => mass_dat(region)%rxym_t rxzm => mass_dat(region)%rxzm_t ryym => mass_dat(region)%ryym_t ryzm => mass_dat(region)%ryzm_t rzzm => mass_dat(region)%rzzm_t distr_mode = 'unfrm' ! 'slope' or 'unfrm' do lrg=1,nred(region) redfact=clustsize(lrg,region) j = jred(lrg,region) do l=1,lmr do i = imredj(lrg,region),1,-1 ! the i cell will be distributed within the is:ie array section is = (i-1)*redfact + 1 ie = i*redfact !m_uni is the mass-distribution in the non-reduced grid/divided by !the reduced_grid mass. This is used as distribution function!.... mass=m(i,j,l); m(is:ie,j,l)= mass/(ie-is+1) if (distr_mode=='unfrm') then ! mixing ratio and x-slope will be UNiFoRMly distributed ! within is:ie #ifdef MPI do n=1,ntracetloc #else do n=1,ntracet #endif rmm = rm(i,j,l,n) rm(is:ie,j,l,n)= m(is:ie,j,l)/mass* rmm rmm = rxm(i,j,l,n) rxm(is:ie,j,l,n)= m(is:ie,j,l)/mass * rmm ! rym and rzm are always distributed uniformly: rmm = rym(i,j,l,n) rym(is:ie,j,l,n)= m(is:ie,j,l)/mass * rmm rmm = rzm(i,j,l,n) rzm(is:ie,j,l,n)= m(is:ie,j,l)/mass * rmm rmm = rxxm(i,j,l,n) rxxm(is:ie,j,l,n)= m(is:ie,j,l)/mass * rmm rmm = rxym(i,j,l,n) rxym(is:ie,j,l,n)= m(is:ie,j,l)/mass * rmm rmm = rxzm(i,j,l,n) rxzm(is:ie,j,l,n)= m(is:ie,j,l)/mass * rmm rmm = ryym(i,j,l,n) ryym(is:ie,j,l,n)= m(is:ie,j,l)/mass * rmm rmm = ryzm(i,j,l,n) ryzm(is:ie,j,l,n)= m(is:ie,j,l)/mass * rmm rmm = rzzm(i,j,l,n) rzzm(is:ie,j,l,n)= m(is:ie,j,l)/mass * rmm end do end if !cmkelseif(distr_mode=='slope') then end do ! update cell(0,...) according to the periodic bc's: #ifdef MPI rm(0,j,l,1:ntracetloc) = rm(im(region),j,l,1:ntracetloc) #else rm(0,j,l,1:ntracet) = rm(im(region),j,l,1:ntracet) #endif rxm(0,j,l,:) = rxm(im(region),j,l,:) rym(0,j,l,:) = rym(im(region),j,l,:) rzm(0,j,l,:) = rzm(im(region),j,l,:) m(0,j,l) = m(im(region),j,l) end do end do nullify(m) nullify(rm) nullify(rxm) nullify(rym) nullify(rzm) nullify(rxxm) nullify(rxym) nullify(rxzm) nullify(ryym) nullify(ryzm) nullify(rzzm) end subroutine red2uni_em_2nd #endif end module redgridZoom