Problems when restarting from a partposit file

[ahilboll]

Using the new FLEXPART 10.3, I have trouble restarting from a previous partposit file --- not partposit_end, but some intermediate timestep. (I needed this functionality while debugging #234 because I didn't want to wait for the model to reach the point where it broke).

While this seems to work, it would be more user-friendly if there were checks in the code like if restart_run .. release_not_required (suggestion by @pesei via e-mail). Optimally, the procedure to restart from a previous state would be documented in the current Pisso et al. GMDD paper.

I could prepare a patch implementing this, by just checking for IPIN/=1 around all the stop statements I mention below, if that would be helpful.

These are the steps I had to do when restarting:

1. I copied header and partposit file to a new output directory,

cp output_partdump_103/header output_debug/
cp output_partdump_103/partposit_20180821160000 output_debug/partposit_end

2. I set IPIN=1 in COMMAND and adapting the simulation start time to the time of my partposit file

3. I removed the release from the RELEASE file by commenting out everything inside the only &RELEASE section and setting NSPEC=0. If I wouldn't do this there were errors about the number of points / releases
4. I had to comment out several sanity checks in the code in order to get rid of several errors regarding number of releases / points / species

   diff --git a/src/advance.f90 b/src/advance.f90
   index 539a473..632e0f9 100644
   --- a/src/advance.f90
   +++ b/src/advance.f90
   @@ -540,8 +540,8 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
              end do
              if (nsp.gt.nspec) then
      ! This should never happen          
   -            write(*,*) 'advance.f90: ERROR: could not find releasepoint'
   -            stop
   +!           write(*,*) 'advance.f90: ERROR: could not find releasepoint'
   +!           stop
              end if
              if (density(nsp).gt.0.) then
                call get_settling(itime,real(xt),real(yt),zt,nsp,settling)  !bugfix
   @@ -710,8 +710,8 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
          end do
          if (nsp.gt.nspec) then
      ! This should never happen          
   -        write(*,*) 'advance.f90: ERROR: could not find releasepoint'
   -        stop
   +!       write(*,*) 'advance.f90: ERROR: could not find releasepoint'
   +!       stop
          end if
          if (density(nsp).gt.0.) then
            call get_settling(itime,real(xt),real(yt),zt,nsp,settling)  !bugfix
   @@ -920,8 +929,8 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
          end do
          if (nsp.gt.nspec) then
      ! This should never happen          
   -        write(*,*) 'advance.f90: ERROR: could not find releasepoint'
   -        stop
   +!       write(*,*) 'advance.f90: ERROR: could not find releasepoint'
   +!       stop
          end if
          if (density(nsp).gt.0.) then
            call get_settling(itime+ldt,real(xt),real(yt),zt,nsp,settling) !bugfix
   
   
   diff --git a/src/coordtrafo.f90 b/src/coordtrafo.f90
   index 2747c60..7d22fb4 100644
   --- a/src/coordtrafo.f90
   +++ b/src/coordtrafo.f90
   @@ -47,7 +47,7 @@ subroutine coordtrafo
      integer :: i,j,k
      real :: yrspc ! small real number relative to x
    
   -  if (numpoint.eq.0) goto 30
   +!  if (numpoint.eq.0) goto 30
    
      ! TRANSFORM X- AND Y- COORDINATES OF STARTING POINTS TO GRID COORDINATES
      !***********************************************************************
   @@ -107,11 +107,11 @@ subroutine coordtrafo
        endif
      end do
    
   -30   if(numpoint.eq.0) then
   -    write(*,*) ' FLEXPART MODEL SUBROUTINE COORDTRAFO: ERROR ! '
   -    write(*,*) ' NO PARTICLE RELEASES ARE DEFINED!'
   -    write(*,*) ' CHECK FILE RELEASES...'
   -    stop
   -  endif
   +!30   if(numpoint.eq.0) then
   +!    write(*,*) ' FLEXPART MODEL SUBROUTINE COORDTRAFO: ERROR ! '
   +!    write(*,*) ' NO PARTICLE RELEASES ARE DEFINED!'
   +!    write(*,*) ' CHECK FILE RELEASES...'
   +!    stop
   +!  endif
    
    end subroutine coordtrafo
   
   
   
   diff --git a/src/readpartpositions.f90 b/src/readpartpositions.f90
   index 3e98ded..180d6bd 100644
   --- a/src/readpartpositions.f90
   +++ b/src/readpartpositions.f90
   @@ -67,13 +67,13 @@ subroutine readpartpositions
      read(unitpartin)
      read(unitpartin) nspecin
      nspecin=nspecin/3
   -  if ((ldirect.eq.1).and.(nspec.ne.nspecin)) goto 997
   +!  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
   +!    if ((ldirect.eq.1).and.(species(i)(1:7).ne.specin)) goto 996
      end do
    
      read(unitpartin) numpointin
   

5. Now the model was running, but in the next timestep, all particles were terminated due to mass, so I commented out that part as well. This seems to be a bug on its own.

diff --git a/src/timemanager.f90 b/src/timemanager.f90
index 58d8a89..f515254 100644
--- a/src/timemanager.f90
+++ b/src/timemanager.f90
@@ -669,12 +669,12 @@ subroutine timemanager(metdata_format)
             end if
           end do
 
-          if (xmassfract.lt.minmass) then   ! terminate all particles carrying less mass
-            itra1(j)=-999999999
-            if (verbosity.gt.0) then
-              print*,'terminated particle ',j,' for small mass'
-            endif
-          endif
+!         if (xmassfract.lt.minmass) then   ! terminate all particles carrying less mass
+!           itra1(j)=-999999999
+!           if (verbosity.gt.0) then
+!             print*,'terminated particle ',j,' for small mass'
+!           endif
+!         endif
 
   !        Sabine Eckhardt, June 2008
   !        don't create depofield for backward runs

[ahilboll]

I have to add another problem: While the model runs and completes, I have trouble reading the output. The netCDF grid_conc_ file doesn't contain any spec_ variables. So while I can run the model (and use this run to verify that it doesn't crash), I cannot use the model output.

[pesei]

What's the status?

[ahilboll]

I didn't further look into this. Just needed this to get quicker to the breaking point in #234 in order to debug that one more efficiently. Once I fixed that one, I just re-ran my whole simulation and didn't need the restart functionality.

[pesei]

See also #257.