Index: initdyn9.ftn
===================================================================
--- initdyn9.ftn	2004-11-15 19:18:49.000000000 +0000
+++ initdyn9.ftn	2004-11-15 19:18:18.000000000 +0000
@@ -37,16 +37,9 @@ copyright (C) 2001  MSC-RPN COMM  %%%MC2
       pointer (pas1 ,s(*)), (pad1 ,d(*))
       integer id,i,m,dim,pndtime
       real prdts,prepsi
+      real, dimension(minx:maxx,miny:maxy,0:gnk) :: ppOrig
 *---------------------------------------------------------------------
 *
-*     * Is total initialization time greater than interval
-*     * between 2 consecutive sets of nesting data?
-*
-      if ((grninit*gndtini).gt.(Pil_nesdt)) then
-         if (myproc.eq.0) write (6,400) grninit*gndtini,real(Pil_nesdt)
-         call mc2stop(-1)
-      endif
-*
       dim      = ndynvar*dim3d+dim2d
       prepsi   = grepsi
       grepsi   = 1.0
@@ -55,10 +48,14 @@ copyright (C) 2001  MSC-RPN COMM  %%%MC2
       pndtime  = vmh_stime
       vmh_stime= 2*gndtini
 *
-      do 20 id=1,2
-         do 10 m=1,gndtini
+* Save original mass field for re-population after each step
+      ppOrig = pp0
+
+      do 10 m=1,gndtini
+         do 20 id=1,2
+
             gnstepno=m
-            if (myproc.eq.0) write (6,910) (id*m)+(id-1)*(gndtini-m),
+            if (myproc.eq.0) write (6,910) id+(m-1)*2-1,
      $                                     2*gndtini
 *
             call step8(0,trname,nt)
@@ -67,16 +64,20 @@ copyright (C) 2001  MSC-RPN COMM  %%%MC2
 *
             call tfilt5 (myproc.eq.0)
 *PDIR RELEASE
- 10      continue
 *
-         pas1 = pappp
-         pad1 = pappm
-         do i = 1, dim
-            d(i) = s(i)
-         end do
+	 pas1 = pappp
+	 pad1 = pappm
+	 do i=1,dim
+	   d(i) = s(i)
+	 enddo
          grdt = -grdt
 *
- 20   continue
+ 20     continue 
+*
+* Repopulate the mass field after step (see Nitta and Hovermale 1969 - MWR)
+	ppp = ppOrig
+*
+ 10   continue
 *
       grdt      = prdts
       grepsi    = prepsi