Changeset 4598 for branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPATAM_SRC/mt19937ar.f90

Timestamp:

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

Author:

pabouttier

Message:

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

File:

• branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPATAM_SRC/mt19937ar.f90 (modified) (4 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