Changeset 4598 for branches

Timestamp:

2014-03-27T09:56:56+01:00 (10 years ago)

Author:

pabouttier

Message:

Add routines to allow restartability of PRNG, see Ticket #1284

Location:

branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPATAM_SRC

Files:

• mt19937ar.f90 (modified) (4 diffs)
• ran_num.F90 (modified) (3 diffs)

branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPATAM_SRC/mt19937ar.f90

@@ -59 +59 @@
    !-----------------------------------------------------------------------!
    USE par_kind, ONLY : wp
+   USE dom_oce , ONLY : nproc
+   USE lib_mpp , ONLY : get_unit
+
    PRIVATE
    INTEGER, PARAMETER :: jpn = 624           ! period parameters
@@ -85 +88 @@
       & mtrand_real53,&
       & init_by_array,&
-      & init_mtrand
+      & init_mtrand,  &
+      & mtrand_seeddump, &
+      & mtrand_seedread
 CONTAINS
    SUBROUTINE init_mtrand(ks)
@@ -93 +98 @@
       INTEGER :: ks
       !
-      mt(0) = IAND(INT(ks,kind=8),jpall)
-      DO mti = 1, jpn-1
-         mt(mti) = 1812433253 * IEOR(mt(mti-1),ISHFT(mt(mti-1),-30)) + mti
-         ! See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. 
-         ! In the previous versions, MSBs of the seed affect   
-         ! only MSBs of the array mt[].                        
-         ! modified by Makoto Matsumoto 
-         mt(mti) = IAND(mt(mti),jpall)
-         ! for >32 bit machines
-      END DO
-      l_mtinit = .TRUE.
-!
+      IF (.NOT.l_mtinit) THEN
+         mt(0) = IAND(INT(ks,kind=8),jpall)
+         DO mti = 1, jpn-1
+            mt(mti) = 1812433253 * IEOR(mt(mti-1),ISHFT(mt(mti-1),-30)) + mti
+            ! See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. 
+            ! In the previous versions, MSBs of the seed affect   
+            ! only MSBs of the array mt[].                        
+            ! modified by Makoto Matsumoto 
+            mt(mti) = IAND(mt(mti),jpall)
+            ! for >32 bit machines
+         END DO
+         l_mtinit = .TRUE.
+      ENDIF
    END SUBROUTINE init_mtrand
+
    SUBROUTINE init_by_array(kinit_key,key_length)
 !-----------------------------------------------------------------------
@@ -231 +238 @@
       mtrand_res53 = ( za * 67108864._wp + zb ) / 9007199254740992._wp
    END FUNCTION mtrand_res53
+
+
+   SUBROUTINE mtrand_seeddump(cdfile)
+!-----------------------------------------------------------------------
+!     Dumps the seed in a restart file.
+!-----------------------------------------------------------------------
+      CHARACTER(len=*) :: cdfile
+      INTEGER :: iunit
+
+      iunit = get_unit()
+      OPEN( unit=iunit, file = TRIM(cdfile), form = 'unformatted', &
+         &  access = 'sequential', status = 'new' )
+      WRITE( iunit ) mt
+      WRITE( iunit ) mti
+      CLOSE( iunit )
+
+   END SUBROUTINE mtrand_seeddump
+
+   SUBROUTINE mtrand_seedread(cdfile)
+!-----------------------------------------------------------------------
+!     Reads the seed from a restart file.
+!-----------------------------------------------------------------------
+      CHARACTER(len=*) :: cdfile
+      INTEGER :: iunit
+
+      iunit = get_unit()
+      OPEN( unit=iunit, file = TRIM(cdfile), form = 'unformatted', &
+         &  access = 'sequential', status = 'old' )
+      READ( iunit ) mt
+      READ( iunit ) mti
+      CLOSE( iunit )
+      l_mtinit = .TRUE.
+
+   END SUBROUTINE mtrand_seedre
+   
 END MODULE mt19937ar

branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPATAM_SRC/ran_num.F90

@@ -15 +15 @@
    USE mt19937ar, ONLY : &
     & init_mtrand,       &
-    & mtrand_real1
+    & mtrand_real1,      &
+    & mtrand_seeddump,   &
+    & mtrand_seedread
  
    IMPLICIT NONE
@@ -25 +27 @@
    PUBLIC &
       & gaustb, &
-      & gaustb_2d
+      & gaustb_2d, &
+      & random_seeddump, &
+      & random_seedread
 
 CONTAINS
@@ -260 +264 @@
    END FUNCTION psrandom
 
+   
+   SUBROUTINE random_seeddump(cdfile)
+      !!----------------------------------------------------------------------
+      !!               ***  ROUTINE psrandom ***
+      !!
+      !! ** Purpose : Pseudo-Random number generator.
+      !!
+      !! ** Method  : Dump the seed
+      !!
+      !! ** Action  :
+      !!
+      !! History :
+      !!        ! 2013-03 (K. Mogensen) Original code
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      CHARACTER(len=128) :: cdfile
+      CALL mtrand_seeddump(cdfile)
+
+   END SUBROUTINE random_seeddump
+
+   SUBROUTINE random_seedread(cdfile)
+      !!----------------------------------------------------------------------
+      !!               ***  ROUTINE psrandom ***
+      !!
+      !! ** Purpose : Pseudo-Random number generator.
+      !!
+      !! ** Method  : Read the seed
+      !!
+      !! ** Action  :
+      !!
+      !! History :
+      !!        ! 2013-03 (K. Mogensen) Original code
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      CHARACTER(len=128) :: cdfile
+      CALL mtrand_seedread(cdfile)
+
+   END SUBROUTINE random_seedread
 
 END MODULE ran_num