Parallel version can now save/restart simulations with IPOUT/IPIN

src/FLEXPART_MPI.f90

@@ -316,7 +316,8 @@ program flexpart
     if (verbosity.gt.0 .and. lroot) then
       print*,'call readpartpositions'
     endif
-    call readpartpositions
+    ! readwind process skips this step
+    if (lmp_use_reader.and..not.lmpreader) call readpartpositions
   else
     if (verbosity.gt.0 .and. lroot) then
       print*,'numpart=0, numparticlecount=0'

src/makefile

@@ -13,7 +13,7 @@ SHELL = /bin/bash
 #    also set environment variable LD_LIBRARY_PATH to point to compiler libraries
 # 
 #  USAGE
-#    Compile serial FLEXPART (ECMWF)a
+#    Compile serial FLEXPART (ECMWF)
 #      make [-j] ecmwf
 #
 #    Compile parallel FLEXPART (ECMWF)

src/mpi_mod.f90

@@ -164,7 +164,7 @@ contains
 !   mp_np       number of running processes, decided at run-time
 !***********************************************************************
     use par_mod, only: maxpart, numwfmem, dep_prec
-    use com_mod, only: mpi_mode
+    use com_mod, only: mpi_mode, verbosity
 
     implicit none
 
@@ -385,6 +385,9 @@ contains
 !*******************************************************************************    
 ! Dynamic allocation of arrays (in serial code these have size maxpart)
 !
+! Not in use anymore, moved to com_mod for interoperability with serial version
+! 
+! TODO: remove 
 !*******************************************************************************
     use com_mod
 
@@ -1042,7 +1045,7 @@ contains
 !
 !***********************************************************************
     use com_mod
-    use par_mod, only: numwfmem
+    use par_mod,only: numwfmem
 
     implicit none
 
@@ -1225,7 +1228,7 @@ contains
     integer :: d3s2 = nxmax*nymax*nuvzmax
     integer :: d2s1 = nxmax*nymax
     integer :: d2s2 = nxmax*nymax*maxspec
-    integer :: d1s1 = maxwf
+!    integer :: d1s1 = maxwf
 
 !*******************************************************************************
 
@@ -1429,6 +1432,10 @@ contains
     else if (memstat.eq.19) then
 ! last read was asynchronous to index 2, 3 is free
       mind=3
+    else
+! illegal state
+      write(*,*) 'mpi_mod> FLEXPART ERROR: Illegal memstat value. Exiting.'
+      stop
     end if
 
     if (mp_dev_mode) then
@@ -1601,13 +1608,12 @@ contains
 !
 !
 !*******************************************************************************
-    use com_mod, only: readclouds
+!    use com_mod,only: readclouds
 
     implicit none
 
 
-    integer :: n_req
-    integer :: i
+    integer :: n_req !,i
 
 !***********************************************************************
 
@@ -1957,7 +1963,7 @@ contains
 !***********************************************************************
     implicit none
 
-    integer :: ip,j,r
+    integer :: ip !,j,r
 
 !***********************************************************************
 

src/par_mod.f90

@@ -185,7 +185,7 @@ module par_mod
   ! Maximum number of particles, species, and similar
   !**************************************************
 
-  integer,parameter :: maxpart=10000
+  integer,parameter :: maxpart=4000000
   integer,parameter :: maxspec=1
   real,parameter :: minmass=0.0001
 
@@ -228,9 +228,8 @@ module par_mod
   ! Dimension of random number field
   !*********************************
 
-!  integer,parameter :: maxrand=120000000
   integer,parameter :: maxrand=200000000
-!
+
   ! maxrand                 number of random numbers used
   
 

src/readpartpositions_mpi.f90

@@ -44,91 +44,153 @@ subroutine readpartpositions
   use par_mod
   use com_mod
   use random_mod, only: ran1
-  use mpi_mod, only: mp_seed
+  use mpi_mod !, only: mp_seed
 
   implicit none
 
-  integer :: ibdatein,ibtimein,nspecin,itimein,numpointin,i,j,ix
+  integer :: ibdatein,ibtimein,nspecin,itimein,numpointin,i,j,ix,ip
   integer :: id1,id2,it1,it2
+  integer :: addone,numparticlecount_all,numpart_all,lbnd,ubnd
   real :: xlonin,ylatin,topo,hmixi,pvi,qvi,rhoi,tri,tti
   character :: specin*7
   real(kind=dp) :: julin,julpartin,juldate
 
   integer :: idummy = -8
 
+  ! These variables are allocated at the root process for all particles in file.
+  real,dimension(maxpart) :: xtra1_all,ytra1_all,ztra1_all
+  real,dimension(maxpart,maxspec) :: xmass1_all
+  integer,dimension(maxpart) :: npoint_all,itramem_all
+
   ! Different seed for each process
   idummy=idummy+mp_seed
 
   numparticlecount=0
+  numparticlecount_all=0
+  numpart_all=0
 
   ! Open header file of dumped particle data
-  !*****************************************
-
-  open(unitpartin,file=path(2)(1:length(2))//'header', &
-       form='unformatted',err=998)
-
-  read(unitpartin) ibdatein,ibtimein
-  read(unitpartin)
-  read(unitpartin)
+  ! Each MPI process sequentially access file (just in case)
+  !*********************************************************
 
-  read(unitpartin)
-  read(unitpartin)
-  read(unitpartin) nspecin
-  nspecin=nspecin/3
-  if ((ldirect.eq.1).and.(nspec.ne.nspecin)) goto 997
+  do ip=0, mp_partgroup_np-1
+    call mpif_mpi_barrier
+    open(unitpartin,file=path(2)(1:length(2))//'header', &
+         form='unformatted',err=998)
 
-  do i=1,nspecin
+    read(unitpartin) ibdatein,ibtimein
     read(unitpartin)
     read(unitpartin)
-    read(unitpartin) j,specin
-    if ((ldirect.eq.1).and.(species(i)(1:7).ne.specin)) goto 996
-  end do
 
-  read(unitpartin) numpointin
-  if (numpointin.ne.numpoint) goto 995
-999 continue 
-  do i=1,numpointin
     read(unitpartin)
     read(unitpartin)
-    read(unitpartin)
-    read(unitpartin)
-    do j=1,nspec
-      read(unitpartin)
+    read(unitpartin) nspecin
+    nspecin=nspecin/3
+    if ((ldirect.eq.1).and.(nspec.ne.nspecin)) goto 997
+
+    do i=1,nspecin
       read(unitpartin)
       read(unitpartin)
+      read(unitpartin) j,specin
+      if ((ldirect.eq.1).and.(species(i)(1:7).ne.specin)) goto 996
     end do
-  end do
-  read(unitpartin)
-  read(unitpartin)
+    
+    read(unitpartin) numpointin
+    if (numpointin.ne.numpoint) then ! goto 995
+
+! eso 2016: moved this warning here to avoid out-of-block goto
+!995   write(*,*) ' #### FLEXPART MODEL WARNING IN READPARTPOSITIONS#### '
+      write(*,*) ' #### FLEXPART MODEL WARNING IN READPARTPOSITIONS#### '
+      write(*,*) ' #### NUMBER OF RELEASE LOCATIONS DOES NOT     #### '
+      write(*,*) ' #### AGREE WITH CURRENT SETTINGS!             #### '
+!  stop
+      goto 999 
+    end if
 
-  do ix=0,numxgrid-1
+999 continue 
+    do i=1,numpointin
+      read(unitpartin)
+      read(unitpartin)
+      read(unitpartin)
+      read(unitpartin)
+      do j=1,nspec
+        read(unitpartin)
+        read(unitpartin)
+        read(unitpartin)
+      end do
+    end do
     read(unitpartin)
-  end do
+    read(unitpartin)
+
+    do ix=0,numxgrid-1
+      read(unitpartin)
+    end do
+    close(unitpartin)
 
 
   ! Open data file of dumped particle data
+  ! All processes read the whole file
   !***************************************
 
-  close(unitpartin)
   open(unitpartin,file=path(2)(1:length(2))//'partposit_end', &
-       form='unformatted',err=998)
-
-
-100   read(unitpartin,end=99) itimein
+         form='unformatted',err=998)
+    
+100 read(unitpartin,end=99) itimein
     i=0
-200   i=i+1
-    read(unitpartin) npoint(i),xlonin,ylatin,ztra1(i),itramem(i), &
-         topo,pvi,qvi,rhoi,hmixi,tri,tti,(xmass1(i,j),j=1,nspec)
-
+200 i=i+1
+    read(unitpartin) npoint_all(i),xlonin,ylatin,ztra1_all(i),itramem_all(i), &
+         topo,pvi,qvi,rhoi,hmixi,tri,tti,(xmass1_all(i,j),j=1,nspec)
+    
     if (xlonin.eq.-9999.9) goto 100
-    xtra1(i)=(xlonin-xlon0)/dx
-    ytra1(i)=(ylatin-ylat0)/dy
-    numparticlecount=max(numparticlecount,npoint(i))
+    xtra1_all(i)=(xlonin-xlon0)/dx
+    ytra1_all(i)=(ylatin-ylat0)/dy
+    numparticlecount_all=max(numparticlecount_all,npoint(i))
     goto 200
 
-99   numpart=i-1
+99  numpart_all=i-1
+    
+    close(unitpartin)
 
-  close(unitpartin)
+  end do
+
+
+  ! Distribute particles among processes
+  !**************************************
+  lbnd=1
+  ubnd=0
+
+  do ip=0, mp_partgroup_np-1
+    
+! Extra particle distributed in case remainder in the division
+    if (ip.lt.mod(numpart_all,mp_partgroup_np)) then
+      addone=1
+    else
+      addone=0
+    end if
+
+    if (ip==0) then 
+      ubnd=ubnd + numpart_all/mp_partgroup_np + addone
+    else 
+      ubnd=lbnd + numpart_all/mp_partgroup_np + addone - 1
+    end if
+    
+    if (ip==mp_pid) then
+      
+      numparticlecount=numparticlecount_all/mp_partgroup_np+addone
+      numpart=numpart_all/mp_partgroup_np+addone
+
+      xtra1(1:numpart) = xtra1_all(lbnd:ubnd)
+      ytra1(1:numpart) = ytra1_all(lbnd:ubnd)
+      ztra1(1:numpart) = ztra1_all(lbnd:ubnd)
+      xmass1(1:numpart,:) = xmass1_all(lbnd:ubnd,:)
+      npoint(1:numpart) = npoint_all(lbnd:ubnd)
+      itramem(1:numpart) = itramem_all(lbnd:ubnd)
+
+    end if
+
+    lbnd=ubnd+1
+
+  end do
 
   julin=juldate(ibdatein,ibtimein)+real(itimein,kind=dp)/86400._dp
   if (abs(julin-bdate).gt.1.e-5) goto 994
@@ -155,13 +217,6 @@ subroutine readpartpositions
   write(*,*) 'bdate: ',bdate,id2,it2
   stop
 
-!995   write(*,*) ' #### FLEXPART MODEL ERROR IN READPARTPOSITIONS#### '
-995   write(*,*) ' #### FLEXPART MODEL WARNING IN READPARTPOSITIONS#### '
-  write(*,*) ' #### NUMBER OF RELEASE LOCATIONS DOES NOT     #### '
-  write(*,*) ' #### AGREE WITH CURRENT SETTINGS!             #### '
-!  stop
-  goto 999 
-
 996   write(*,*) ' #### FLEXPART MODEL ERROR IN READPARTPOSITIONS#### '
   write(*,*) ' #### SPECIES NAMES TO BE READ IN DO NOT       #### '
   write(*,*) ' #### AGREE WITH CURRENT SETTINGS!             #### '

src/readspecies.f90

@@ -150,7 +150,7 @@ subroutine readspecies(id_spec,pos_spec)
     read(unitspecies,'(e18.1)',end=22) density(pos_spec)
 !  write(*,*) density(pos_spec)
     read(unitspecies,'(e18.1)',end=22) dquer(pos_spec)
-    write(*,*) 'dquer(pos_spec):', dquer(pos_spec)
+!    write(*,*) 'dquer(pos_spec):', dquer(pos_spec)
     read(unitspecies,'(e18.1)',end=22) dsigma(pos_spec)
 !  write(*,*) dsigma(pos_spec)
     read(unitspecies,'(f18.2)',end=22) dryvel(pos_spec)

src/timemanager_mpi.f90

@@ -103,7 +103,7 @@ subroutine timemanager
   implicit none
 
   logical :: reqv_state=.false. ! .true. if waiting for a MPI_Irecv to complete
-  integer :: j,ks,kp,l,n,itime=0,nstop,nstop1,memstat=0,mind
+  integer :: j,ks,kp,l,n,itime=0,nstop,nstop1,memstat=0 !,mind
 ! integer :: ksp
   integer :: ip
   integer :: loutnext,loutstart,loutend
@@ -114,7 +114,8 @@ subroutine timemanager
   real :: decfact
 
   real(sp) :: gridtotalunc
-  real(dep_prec) :: drygridtotalunc,wetgridtotalunc,drydeposit(maxspec)
+  real(dep_prec) :: drygridtotalunc=0_dep_prec,wetgridtotalunc=0_dep_prec,&
+       & drydeposit(maxspec)=0_dep_prec
   real :: xold,yold,zold,xmassfract
   real, parameter :: e_inv = 1.0/exp(1.0)
 
@@ -136,9 +137,9 @@ subroutine timemanager
 
 
 !  itime=0
-  if (lroot) then
-  !  write(*,45) itime,numpart*mp_partgroup_np,gridtotalunc,wetgridtotalunc,drygridtotalunc
-    write(*,46) float(itime)/3600,itime,numpart*mp_partgroup_np
+  if (lroot.or.mp_dev_mode) then
+    write(*,45) itime,numpart*mp_partgroup_np,gridtotalunc,wetgridtotalunc,drygridtotalunc
+  !  write(*,46) float(itime)/3600,itime,numpart*mp_partgroup_np
     
     if (verbosity.gt.0) then
       write (*,*) 'timemanager> starting simulation'
@@ -274,7 +275,7 @@ subroutine timemanager
 
 ! For validation and tests: call the function below to set all fields to simple
 ! homogeneous values
-    if (mp_dev_mode) call set_fields_synthetic
+!    if (mp_dev_mode) call set_fields_synthetic
 
 !*******************************************************************************
 
@@ -553,7 +554,7 @@ subroutine timemanager
         endif
         
         !CGZ-lifetime: output species lifetime
-        if (lroot) then
+        if (lroot.or.mp_dev_mode) then
         !   write(*,*) 'Overview species lifetime in days', &
         !        real((species_lifetime(:,1)/species_lifetime(:,2))/real(3600.0*24.0))
         !   write(*,*) 'all info:',species_lifetime