Changeset 0a94e13 in flexpart.git


Ignore:
Timestamp:
May 6, 2019, 11:43:21 AM (5 years ago)
Author:
Espen Sollum ATMOS <eso@…>
Branches:
master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug
Children:
0c8c7f2
Parents:
328fdf9
Message:

Added ipout=3 option for time averaged particle output

Files:
15 edited

Legend:

Unmodified
Added
Removed
  • options/COMMAND

    r2753a5c r0a94e13  
    1919 IFINE=                 4, ! Reduction for time step in vertical transport, used only if CTL>1
    2020 IOUT=                  1, ! Output type: [1]mass 2]pptv 3]1&2 4]plume 5]1&4, +8 for NetCDF output     
    21  IPOUT=                 0, ! Particle position output: 0]no 1]every output 2]only at end  
     21 IPOUT=                 0, ! Particle position output: 0]no 1]every output 2]only at end 3]time averaged
    2222 LSUBGRID=              0, ! Increase of ABL heights due to sub-grid scale orographic variations;[0]off 1]on
    2323 LCONVECTION=           1, ! Switch for convection parameterization;0]off [1]on   
  • src/FLEXPART.f90

    r2753a5c r0a94e13  
    6868  integer :: detectformat
    6969
    70 
    71 
    72   ! Initialize arrays in com_mod
    73   !*****************************
    74   call com_mod_allocate_part(maxpart)
    7570 
    76 
    7771  ! Generate a large number of random numbers
    7872  !******************************************
     
    172166  endif
    173167
     168  ! Initialize arrays in com_mod
     169  !*****************************
     170  call com_mod_allocate_part(maxpart)
     171
     172
    174173  ! Read the age classes to be used
    175174  !********************************
  • src/FLEXPART_MPI.f90

    r20963b1 r0a94e13  
    7777  if (mp_measure_time) call mpif_mtime('flexpart',0)
    7878
    79   ! Initialize arrays in com_mod
    80   !*****************************
    81 
    82   if(.not.(lmpreader.and.lmp_use_reader)) call com_mod_allocate_part(maxpart_mpi)
    83 
    84 
     79 
    8580  ! Generate a large number of random numbers
    8681  !******************************************
     
    180175  endif
    181176
     177  ! Initialize arrays in com_mod
     178  !*****************************
     179
     180  if(.not.(lmpreader.and.lmp_use_reader)) call com_mod_allocate_part(maxpart_mpi)
     181
    182182
    183183! Read the age classes to be used
  • src/com_mod.f90

    r328fdf9 r0a94e13  
    1818
    1919  implicit none
     20
     21
    2022
    2123  !****************************************************************
     
    6971
    7072  real :: ctl,fine
    71   integer :: ifine,iout,ipout,ipin,iflux,mdomainfill
     73  integer :: ifine,iout,ipout,ipin,iflux,mdomainfill,ipoutfac
    7274  integer :: mquasilag,nested_output,ind_source,ind_receptor
    7375  integer :: ind_rel,ind_samp,ioutputforeachrelease,linit_cond,surf_only
     
    8284  ! iout     output options: 1 conc. output (ng/m3), 2 mixing ratio (pptv), 3 both
    8385  ! ipout    particle dump options: 0 no, 1 every output interval, 2 only at end
     86  ! ipoutfac increase particle dump interval by factor (default 1)
    8487  ! ipin     read in particle positions from dumped file from a previous run
    8588  ! fine     real(ifine)
     
    128131
    129132  logical :: gdomainfill
    130 
    131133  ! gdomainfill             .T., if domain-filling is global, .F. if not
    132134
     
    675677  real, allocatable, dimension(:,:) :: xscav_frac1
    676678
     679! Variables used for writing out interval averages for partoutput
     680!****************************************************************
     681
     682  integer, allocatable, dimension(:) :: npart_av
     683  real, allocatable, dimension(:) :: part_av_cartx,part_av_carty,part_av_cartz,part_av_z,part_av_topo
     684  real, allocatable, dimension(:) :: part_av_pv,part_av_qv,part_av_tt,part_av_rho,part_av_tro,part_av_hmix
     685  real, allocatable, dimension(:) :: part_av_uu,part_av_vv,part_av_energy
     686
    677687  ! eso: Moved from timemanager
    678688  real, allocatable, dimension(:) :: uap,ucp,uzp,us,vs,ws
     
    781791         & idt(nmpart),itramem(nmpart),itrasplit(nmpart),&
    782792         & xtra1(nmpart),ytra1(nmpart),ztra1(nmpart),&
    783          & xmass1(nmpart, maxspec),&
    784          & checklifetime(nmpart,maxspec), species_lifetime(maxspec,2))!CGZ-lifetime
     793         & xmass1(nmpart, maxspec))  ! ,&
     794!         & checklifetime(nmpart,maxspec), species_lifetime(maxspec,2))!CGZ-lifetime
     795
     796    if (ipout.eq.3) then
     797      allocate(npart_av(nmpart),part_av_cartx(nmpart),part_av_carty(nmpart),&
     798           & part_av_cartz(nmpart),part_av_z(nmpart),part_av_topo(nmpart))
     799      allocate(part_av_pv(nmpart),part_av_qv(nmpart),part_av_tt(nmpart),&
     800           & part_av_rho(nmpart),part_av_tro(nmpart),part_av_hmix(nmpart))
     801      allocate(part_av_uu(nmpart),part_av_vv(nmpart),part_av_energy(nmpart))
     802    end if
    785803
    786804
    787805    allocate(uap(nmpart),ucp(nmpart),uzp(nmpart),us(nmpart),&
    788806         & vs(nmpart),ws(nmpart),cbt(nmpart))
    789    
     807
    790808  end subroutine com_mod_allocate_part
    791809
  • src/init_domainfill.f90

    r328fdf9 r0a94e13  
    8686    endif
    8787  endif
     88
     89! Exit here if resuming a run from particle dump
     90!***********************************************
     91  if (gdomainfill.and.ipin.ne.0) return
    8892
    8993! Do not release particles twice (i.e., not at both in the leftmost and rightmost
  • src/makefile

    r7123c70 r0a94e13  
    118118OBJECTS_SERIAL = \
    119119        releaseparticles.o      partoutput.o \
     120        partoutput_average.o \
    120121        conccalc.o \
    121122        init_domainfill.o       concoutput.o  \
     
    132133## For MPI version
    133134OBJECTS_MPI = releaseparticles_mpi.o partoutput_mpi.o \
    134         conccalc_mpi.o \
     135        partoutput_average_mpi.o conccalc_mpi.o \
    135136        init_domainfill_mpi.o concoutput_mpi.o  \
    136137        timemanager_mpi.o FLEXPART_MPI.o        \
     
    149150advance.o               initialize.o            \
    150151writeheader.o           writeheader_txt.o       \
    151 writeprecip.o \
     152partpos_average.o       writeprecip.o \
    152153writeheader_surf.o      assignland.o\
    153154part0.o                 gethourlyOH.o\
     
    348349part0.o: par_mod.o
    349350partdep.o: par_mod.o
     351partpos_average.o: com_mod.o par_mod.o
    350352partoutput.o: com_mod.o par_mod.o
     353partoutput_average.o: com_mod.o par_mod.o
     354partoutput_average_mpi.o: com_mod.o par_mod.o mpi_mod.o
    351355partoutput_mpi.o: com_mod.o mpi_mod.o par_mod.o
    352356partoutput_short.o: com_mod.o par_mod.o
  • src/mpi_mod.f90

    • Property mode changed from 100644 to 100755
    r328fdf9 r0a94e13  
    8888! Variables for MPI processes in the 'particle' group
    8989  integer, allocatable, dimension(:) :: mp_partgroup_rank
     90  integer, allocatable, dimension(:) :: npart_per_process
    9091  integer :: mp_partgroup_comm, mp_partgroup_pid, mp_partgroup_np
    9192
     
    125126! mp_time_barrier   Measure MPI barrier time
    126127! mp_exact_numpart  Use an extra MPI communication to give the exact number of particles
    127 !                   to standard output (this does *not* otherwise affect the simulation)
    128 ! mp_rebalance      Attempt to rebalance particle between processes
     128!                   to standard output (this does not otherwise affect the simulation)
    129129  logical, parameter :: mp_dbg_mode = .false.
    130130  logical, parameter :: mp_dev_mode = .false.
     
    133133  logical, parameter :: mp_measure_time=.false.
    134134  logical, parameter :: mp_exact_numpart=.true.
    135   logical, parameter :: mp_rebalance=.true.
    136135
    137136! for measuring CPU/Wall time
     
    146145  real(dp),private :: mp_getfields_wtime_beg, mp_getfields_wtime_end, mp_getfields_wtime_total=0.
    147146  real(sp),private :: mp_getfields_time_beg, mp_getfields_time_end, mp_getfields_time_total=0.
    148 !  real(dp),private :: mp_readwind_wtime_beg, mp_readwind_wtime_end, mp_readwind_wtime_total=0.
    149 !  real(sp),private :: mp_readwind_time_beg, mp_readwind_time_end, mp_readwind_time_total=0.
     147  real(dp),private :: mp_readwind_wtime_beg, mp_readwind_wtime_end, mp_readwind_wtime_total=0.
     148  real(sp),private :: mp_readwind_time_beg, mp_readwind_time_end, mp_readwind_time_total=0.
    150149  real(dp),private :: mp_io_wtime_beg, mp_io_wtime_end, mp_io_wtime_total=0.
    151150  real(sp),private :: mp_io_time_beg, mp_io_time_end, mp_io_time_total=0.
     
    367366    end if
    368367
     368! Allocate array for number of particles per process   
     369    allocate(npart_per_process(0:mp_partgroup_np-1))
     370
    369371! Write whether MPI_IN_PLACE is used or not
    370372#ifdef USE_MPIINPLACE
     
    607609    integer :: i,jj,nn, num_part=1,m,imin, num_trans
    608610    logical :: first_iter
    609     integer,allocatable,dimension(:) :: numparticles_mpi, idx_arr
     611    integer,allocatable,dimension(:) :: idx_arr
    610612    real,allocatable,dimension(:) :: sorted ! TODO: we don't really need this
    611613
     
    616618! All processes exchange information on number of particles
    617619!****************************************************************************
    618     allocate(numparticles_mpi(0:mp_partgroup_np-1), &
    619          &idx_arr(0:mp_partgroup_np-1), sorted(0:mp_partgroup_np-1))
    620 
    621     call MPI_Allgather(numpart, 1, MPI_INTEGER, numparticles_mpi, &
     620    allocate( idx_arr(0:mp_partgroup_np-1), sorted(0:mp_partgroup_np-1))
     621
     622    call MPI_Allgather(numpart, 1, MPI_INTEGER, npart_per_process, &
    622623         & 1, MPI_INTEGER, mp_comm_used, mp_ierr)
    623624
     
    625626! Sort from lowest to highest
    626627! Initial guess: correct order
    627     sorted(:) = numparticles_mpi(:)
     628    sorted(:) = npart_per_process(:)
    628629    do i=0, mp_partgroup_np-1
    629630      idx_arr(i) = i
    630631    end do
     632
     633! Do not rebalance particles for ipout=3   
     634    if (ipout.eq.3) return
    631635
    632636! For each successive element in index array, see if a lower value exists
     
    654658    m=mp_partgroup_np-1 ! index for last sorted process (most particles)
    655659    do i=0,mp_partgroup_np/2-1
    656       num_trans = numparticles_mpi(idx_arr(m)) - numparticles_mpi(idx_arr(i))
     660      num_trans = npart_per_process(idx_arr(m)) - npart_per_process(idx_arr(i))
    657661      if (mp_partid.eq.idx_arr(m).or.mp_partid.eq.idx_arr(i)) then
    658         if ( numparticles_mpi(idx_arr(m)).gt.mp_min_redist.and.&
    659              & real(num_trans)/real(numparticles_mpi(idx_arr(m))).gt.mp_redist_fract) then
     662        if ( npart_per_process(idx_arr(m)).gt.mp_min_redist.and.&
     663             & real(num_trans)/real(npart_per_process(idx_arr(m))).gt.mp_redist_fract) then
    660664! DBG
    661           ! write(*,*) 'mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, numparticles_mpi', &
    662           !      &mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, numparticles_mpi
     665          ! write(*,*) 'mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, npart_per_process', &
     666          !      &mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, npart_per_process
    663667! DBG
    664668          call mpif_redist_part(itime, idx_arr(m), idx_arr(i), num_trans/2)
     
    668672    end do
    669673
    670     deallocate(numparticles_mpi, idx_arr, sorted)
     674    deallocate(idx_arr, sorted)
    671675
    672676  end subroutine mpif_calculate_part_redist
     
    27162720      end if
    27172721
    2718 !    case ('readwind')
    2719 !      if (imode.eq.0) then
    2720 !        call cpu_time(mp_readwind_time_beg)
    2721 !        mp_readwind_wtime_beg = mpi_wtime()
    2722 !      else
    2723 !        call cpu_time(mp_readwind_time_end)
    2724 !        mp_readwind_wtime_end = mpi_wtime()
    2725 !
    2726 !        mp_readwind_time_total = mp_readwind_time_total + &
    2727 !             &(mp_readwind_time_end - mp_readwind_time_beg)
    2728 !        mp_readwind_wtime_total = mp_readwind_wtime_total + &
    2729 !             &(mp_readwind_wtime_end - mp_readwind_wtime_beg)
    2730 !      end if
     2722   case ('readwind')
     2723     if (imode.eq.0) then
     2724       call cpu_time(mp_readwind_time_beg)
     2725       mp_readwind_wtime_beg = mpi_wtime()
     2726     else
     2727       call cpu_time(mp_readwind_time_end)
     2728       mp_readwind_wtime_end = mpi_wtime()
     2729
     2730       mp_readwind_time_total = mp_readwind_time_total + &
     2731            &(mp_readwind_time_end - mp_readwind_time_beg)
     2732       mp_readwind_wtime_total = mp_readwind_wtime_total + &
     2733            &(mp_readwind_wtime_end - mp_readwind_wtime_beg)
     2734     end if
    27312735
    27322736    case ('commtime')
  • src/netcdf_output_mod.f90

    r4ad96c5 r0a94e13  
    9393  character(len=255), parameter :: institution = 'NILU'
    9494
    95   integer            :: tpointer
     95  integer            :: tpointer=0
    9696  character(len=255) :: ncfname, ncfnamen
    9797
  • src/par_mod.f90

    rdf96ea65 r0a94e13  
    280280
    281281  integer,parameter :: unitpath=1, unitcommand=1, unitageclasses=1, unitgrid=1
    282   integer,parameter :: unitavailab=1, unitreleases=88, unitpartout=93
     282  integer,parameter :: unitavailab=1, unitreleases=88, unitpartout=93, unitpartout_average=105
    283283  integer,parameter :: unitpartin=93, unitflux=98, unitouttraj=96
    284284  integer,parameter :: unitvert=1, unitoro=1, unitpoin=1, unitreceptor=1
  • src/partoutput.f90

    rd2a5a83 r0a94e13  
    7171  !**************************************
    7272
    73   if (ipout.eq.1) then
     73  if (ipout.eq.1.or.ipout.eq.3) then
    7474    open(unitpartout,file=path(2)(1:length(2))//'partposit_'//adate// &
    7575         atime,form='unformatted')
  • src/partoutput_mpi.f90

    rd2a5a83 r0a94e13  
    7878  !**************************************
    7979
    80   if (ipout.eq.1) then
     80  if (ipout.eq.1.or.ipout.eq.3) then
    8181    open(unitpartout,file=path(2)(1:length(2))//'partposit_'//adate// &
    8282         atime,form='unformatted',status=file_stat,position='append')
  • src/readcommand.f90

    r20963b1 r0a94e13  
    5050  ! ipin                 1 continue simulation with dumped particle data, 0 no *
    5151  ! ipout                0 no particle dump, 1 every output time, 3 only at end*
     52  ! ipoutfac             increase particle dump interval by factor (default 1) *
    5253  ! itsplit [s]          time constant for particle splitting                  *
    5354  ! loutaver [s]         concentration output is an average over loutaver      *
     
    9798  iout, &
    9899  ipout, &
     100  ipoutfac, &
    99101  lsubgrid, &
    100102  lconvection, &
     
    129131  iout=3
    130132  ipout=0
     133  ipoutfac=1
    131134  lsubgrid=1
    132135  lconvection=1
     
    507510  !****************************************************************
    508511
    509   if ((ipout.ne.0).and.(ipout.ne.1).and.(ipout.ne.2)) then
     512  if ((ipout.ne.0).and.(ipout.ne.1).and.(ipout.ne.2).and.(ipout.ne.3)) then
    510513    write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND:     #### '
    511     write(*,*) ' #### IPOUT MUST BE 1, 2 OR 3!                #### '
     514    write(*,*) ' #### IPOUT MUST BE 0, 1, 2 OR 3!             #### '
    512515    stop
    513516  endif
  • src/timemanager.f90

    rc7d1052 r0a94e13  
    45145145      format(i13,' Seconds simulated: ',i13, ' Particles:    Uncertainty: ',3f7.3)
    45245246      format(' Simulated ',f7.1,' hours (',i13,' s), ',i13, ' particles')
    453         if (ipout.ge.1) call partoutput(itime)    ! dump particle positions
     453        if (ipout.ge.1) then
     454          if (mod(itime,ipoutfac*loutstep).eq.0) call partoutput(itime) ! dump particle positions
     455          if (ipout.eq.3) call partoutput_average(itime) ! dump particle positions
     456        endif
    454457        loutnext=loutnext+loutstep
    455458        loutstart=loutnext-loutaver/2
     
    609612!        write (*,*) 'advance: ',prob(1),xmass1(j,1),ztra1(j)
    610613
     614  ! Calculate average position for particle dump output
     615  !****************************************************
     616
     617        if (ipout.eq.3) call partpos_average(itime,j)
     618
     619
    611620  ! Calculate the gross fluxes across layer interfaces
    612621  !***************************************************
  • src/timemanager_mpi.f90

    r328fdf9 r0a94e13  
    113113  integer :: j,ks,kp,l,n,itime=0,nstop,nstop1,memstat=0
    114114! integer :: ksp
    115   integer :: ip
     115  integer :: ip,irec
    116116  integer :: loutnext,loutstart,loutend
    117117  integer :: ix,jy,ldeltat,itage,nage,idummy
     
    129129! Measure time spent in timemanager
    130130  if (mp_measure_time) call mpif_mtime('timemanager',0)
     131
    131132
    132133! First output for time 0
     
    344345! Check if particles should be redistributed among processes
    345346!***********************************************************
    346     if (mp_rebalance) call mpif_calculate_part_redist(itime)
     347    call mpif_calculate_part_redist(itime)
    347348
    348349
     
    59359446      format(' Simulated ',f7.1,' hours (',i13,' s), ',i13, ' particles')
    594595        if (ipout.ge.1) then
     596          irec=0
    595597          do ip=0, mp_partgroup_np-1
    596             if (ip.eq.mp_partid) call partoutput(itime) ! dump particle positions
     598            if (ip.eq.mp_partid) then
     599              if (mod(itime,ipoutfac*loutstep).eq.0) call partoutput(itime) ! dump particle positions
     600              if (ipout.eq.3) call partoutput_average(itime,irec) ! dump particle positions
     601            endif
     602            if (ipout.eq.3) irec=irec+npart_per_process(ip)
    597603            call mpif_mpi_barrier
    598604          end do
     
    757763        if (mp_measure_time) call mpif_mtime('advance',1)
    758764
     765  ! Calculate average position for particle dump output
     766  !****************************************************
     767
     768        if (ipout.eq.3) call partpos_average(itime,j)
     769
    759770
    760771! Calculate the gross fluxes across layer interfaces
     
    895906    do ip=0, mp_partgroup_np-1
    896907      if (ip.eq.mp_partid) then
    897         !if (mp_dbg_mode) write(*,*) 'call partoutput(itime), proc, mp_partid',ip,mp_partid
     908        if (mp_dbg_mode) write(*,*) 'call partoutput(itime), proc, mp_partid',ip,mp_partid
    898909        call partoutput(itime)    ! dump particle positions
    899910      end if
  • src/verttransform_ecmwf.f90

    r79e0349 r0a94e13  
    7373  use com_mod
    7474  use cmapf_mod, only: cc2gll
    75 !  use mpi_mod
    7675
    7776  implicit none
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG