Changeset 5404 for trunk/NEMOGCM
- Timestamp:
- 2015-06-11T09:20:46+02:00 (9 years ago)
- Location:
- trunk/NEMOGCM/NEMO/OPA_SRC/STO
- Files:
-
- 1 added
- 3 copied
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/STO/stopar.F90
r5390 r5404 263 263 INTEGER :: ios ! Local integer output status for namelist read 264 264 265 IF ( lk_agrif ) CALL ctl_stop('EOS stochastic parametrization is not compatible with AGRIF') 265 266 ! Read namsto namelist : stochastic parameterization 266 267 REWIND( numnam_ref ) ! Namelist namdyn_adv in reference namelist : Momentum advection scheme -
trunk/NEMOGCM/NEMO/OPA_SRC/STO/storng.F90
r5390 r5404 33 33 !! kiss_save, kiss_load, kiss_check, kiss_gamma, kiss_sample 34 34 !!---------------------------------------------------------------------- 35 USE par_kind 36 USE lib_mpp 37 35 38 IMPLICIT NONE 36 39 PRIVATE 37 40 38 41 ! Public functions/subroutines 39 PUBLIC :: kiss, kiss_seed, kiss_state, kiss_ save, kiss_load, kiss_reset, kiss_check42 PUBLIC :: kiss, kiss_seed, kiss_state, kiss_reset ! kiss_save, kiss_load, kiss_check 40 43 PUBLIC :: kiss_uniform, kiss_gaussian, kiss_gamma, kiss_sample 41 44 42 45 ! Default/initial seeds 43 INTEGER(KIND= 8) :: x=1234567890987654321_844 INTEGER(KIND= 8) :: y=362436362436362436_845 INTEGER(KIND= 8) :: z=1066149217761810_846 INTEGER(KIND= 8) :: w=123456123456123456_846 INTEGER(KIND=i8) :: x=1234567890987654321_8 47 INTEGER(KIND=i8) :: y=362436362436362436_8 48 INTEGER(KIND=i8) :: z=1066149217761810_8 49 INTEGER(KIND=i8) :: w=123456123456123456_8 47 50 48 51 ! Parameters to generate real random variates 49 REAL(KIND= 8), PARAMETER :: huge64=9223372036854775808.0_8! +150 REAL(KIND= 8), PARAMETER :: zero=0.0_8, half=0.5_8, one=1.0_8, two=2.0_852 REAL(KIND=wp), PARAMETER :: huge64=9223372036854775808.0 ! +1 53 REAL(KIND=wp), PARAMETER :: zero=0.0, half=0.5, one=1.0, two=2.0 51 54 52 55 ! Variables to store 2 Gaussian random numbers with current index (ig) 53 INTEGER(KIND= 8), SAVE :: ig=154 REAL(KIND= 8), SAVE :: gran1, gran256 INTEGER(KIND=i8), SAVE :: ig=1 57 REAL(KIND=wp), SAVE :: gran1, gran2 55 58 56 59 !!---------------------------------------------------------------------- … … 79 82 ! -------------------------------------------------------------------- 80 83 IMPLICIT NONE 81 INTEGER(KIND= 8) :: kiss, t84 INTEGER(KIND=i8) :: kiss, t 82 85 83 86 t = ISHFT(x,58) + w … … 96 99 97 100 FUNCTION s(k) 98 INTEGER(KIND= 8) :: s, k101 INTEGER(KIND=i8) :: s, k 99 102 s = ISHFT(k,-63) 100 103 END FUNCTION s 101 104 102 105 FUNCTION m(k, n) 103 INTEGER(KIND= 8) :: m, k, n106 INTEGER(KIND=i8) :: m, k, n 104 107 m = IEOR(k, ISHFT(k, n) ) 105 108 END FUNCTION m … … 116 119 !! -------------------------------------------------------------------- 117 120 IMPLICIT NONE 118 INTEGER(KIND= 8) :: ix, iy, iz, iw121 INTEGER(KIND=i8) :: ix, iy, iz, iw 119 122 120 123 x = ix … … 134 137 !! -------------------------------------------------------------------- 135 138 IMPLICIT NONE 136 INTEGER(KIND= 8) :: ix, iy, iz, iw139 INTEGER(KIND=i8) :: ix, iy, iz, iw 137 140 138 141 ix = x … … 161 164 162 165 163 SUBROUTINE kiss_check(check_type) 164 !! -------------------------------------------------------------------- 165 !! *** ROUTINE kiss_check *** 166 !! 167 !! ** Purpose : Check the KISS pseudo-random sequence 168 !! 169 !! ** Method : Check that it reproduces the correct sequence 170 !! from the default seed 171 !! 172 !! -------------------------------------------------------------------- 173 IMPLICIT NONE 174 INTEGER(KIND=8) :: iter, niter, correct, iran 175 CHARACTER(LEN=*) :: check_type 176 LOGICAL :: print_success 177 178 ! Save current state of KISS 179 CALL kiss_save() 180 ! Reset the default seed 181 CALL kiss_reset() 182 183 ! Select check type 184 SELECT CASE(check_type) 185 CASE('short') 186 niter = 5_8 187 correct = 542381058189297533_8 188 print_success = .FALSE. 189 CASE('long') 190 niter = 100000000_8 191 correct = 1666297717051644203_8 ! Check provided by G. Marsaglia 192 print_success = .TRUE. 193 CASE('default') 194 CASE DEFAULT 195 STOP 'Bad check type in kiss_check' 196 END SELECT 197 198 ! Run kiss for the required number of iterations (niter) 199 DO iter=1,niter 200 iran = kiss() 201 ENDDO 202 203 ! Check that last iterate is correct 204 IF (iran.NE.correct) THEN 205 STOP 'Check failed: KISS internal error !!' 206 ELSE 207 IF (print_success) PRINT *, 'Check successful: 100 million calls to KISS OK' 208 ENDIF 209 210 ! Reload the previous state of KISS 211 CALL kiss_load() 212 213 END SUBROUTINE kiss_check 214 215 216 SUBROUTINE kiss_save 217 !! -------------------------------------------------------------------- 218 !! *** ROUTINE kiss_save *** 219 !! 220 !! ** Purpose : Save current state of KISS random number generator 221 !! 222 !! -------------------------------------------------------------------- 223 IMPLICIT NONE 224 225 OPEN(UNIT=30,FILE='.kiss_restart') 226 WRITE(30,*) x 227 WRITE(30,*) y 228 WRITE(30,*) z 229 WRITE(30,*) w 230 CLOSE(30) 231 232 END SUBROUTINE kiss_save 233 234 235 SUBROUTINE kiss_load 236 !! -------------------------------------------------------------------- 237 !! *** ROUTINE kiss_load *** 238 !! 239 !! ** Purpose : Load the saved state of KISS random number generator 240 !! 241 !! -------------------------------------------------------------------- 242 IMPLICIT NONE 243 LOGICAL :: filexists 244 245 INQUIRE(FILE='.kiss_restart',EXIST=filexists) 246 IF (filexists) THEN 247 OPEN(UNIT=30,FILE='.kiss_restart') 248 READ(30,*) x 249 READ(30,*) y 250 READ(30,*) z 251 READ(30,*) w 252 CLOSE(30) 253 ENDIF 254 255 END SUBROUTINE kiss_load 166 ! SUBROUTINE kiss_check(check_type) 167 ! !! -------------------------------------------------------------------- 168 ! !! *** ROUTINE kiss_check *** 169 ! !! 170 ! !! ** Purpose : Check the KISS pseudo-random sequence 171 ! !! 172 ! !! ** Method : Check that it reproduces the correct sequence 173 ! !! from the default seed 174 ! !! 175 ! !! -------------------------------------------------------------------- 176 ! IMPLICIT NONE 177 ! INTEGER(KIND=i8) :: iter, niter, correct, iran 178 ! CHARACTER(LEN=*) :: check_type 179 ! LOGICAL :: print_success 180 181 ! ! Save current state of KISS 182 ! CALL kiss_save() 183 ! ! Reset the default seed 184 ! CALL kiss_reset() 185 186 ! ! Select check type 187 ! SELECT CASE(check_type) 188 ! CASE('short') 189 ! niter = 5_8 190 ! correct = 542381058189297533 191 ! print_success = .FALSE. 192 ! CASE('long') 193 ! niter = 100000000_8 194 ! correct = 1666297717051644203 ! Check provided by G. Marsaglia 195 ! print_success = .TRUE. 196 ! CASE('default') 197 ! CASE DEFAULT 198 ! STOP 'Bad check type in kiss_check' 199 ! END SELECT 200 201 ! ! Run kiss for the required number of iterations (niter) 202 ! DO iter=1,niter 203 ! iran = kiss() 204 ! ENDDO 205 206 ! ! Check that last iterate is correct 207 ! IF (iran.NE.correct) THEN 208 ! STOP 'Check failed: KISS internal error !!' 209 ! ELSE 210 ! IF (print_success) PRINT *, 'Check successful: 100 million calls to KISS OK' 211 ! ENDIF 212 213 ! ! Reload the previous state of KISS 214 ! CALL kiss_load() 215 216 ! END SUBROUTINE kiss_check 217 218 219 ! SUBROUTINE kiss_save 220 ! !! -------------------------------------------------------------------- 221 ! !! *** ROUTINE kiss_save *** 222 ! !! 223 ! !! ** Purpose : Save current state of KISS random number generator 224 ! !! 225 ! !! -------------------------------------------------------------------- 226 ! INTEGER :: inum !! Local integer 227 228 ! IMPLICIT NONE 229 230 ! CALL ctl_opn( inum, '.kiss_restart', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp ) 231 232 ! ! OPEN(UNIT=30,FILE='.kiss_restart') 233 ! WRITE(inum,*) x 234 ! WRITE(inum,*) y 235 ! WRITE(inum,*) z 236 ! WRITE(inum,*) w 237 ! CALL flush(inum) 238 239 ! END SUBROUTINE kiss_save 240 241 242 ! SUBROUTINE kiss_load 243 ! !! -------------------------------------------------------------------- 244 ! !! *** ROUTINE kiss_load *** 245 ! !! 246 ! !! ** Purpose : Load the saved state of KISS random number generator 247 ! !! 248 ! !! -------------------------------------------------------------------- 249 ! IMPLICIT NONE 250 ! LOGICAL :: filexists 251 ! Use ctl_opn routine rather than fortran intrinsic functions 252 ! INQUIRE(FILE='.kiss_restart',EXIST=filexists) 253 ! IF (filexists) THEN 254 ! OPEN(UNIT=30,FILE='.kiss_restart') 255 ! READ(30,*) x 256 ! READ(30,*) y 257 ! READ(30,*) z 258 ! READ(30,*) w 259 ! CLOSE(30) 260 ! ENDIF 261 262 ! END SUBROUTINE kiss_load 256 263 257 264 … … 264 271 !! -------------------------------------------------------------------- 265 272 IMPLICIT NONE 266 REAL(KIND= 8) :: uran267 268 uran = half * ( one + REAL(kiss(), 8) / huge64 )273 REAL(KIND=wp) :: uran 274 275 uran = half * ( one + REAL(kiss(),wp) / huge64 ) 269 276 270 277 END SUBROUTINE kiss_uniform … … 284 291 !! -------------------------------------------------------------------- 285 292 IMPLICIT NONE 286 REAL(KIND= 8) :: gran, u1, u2, rsq, fac293 REAL(KIND=wp) :: gran, u1, u2, rsq, fac 287 294 288 295 IF (ig.EQ.1) THEN 289 296 rsq = two 290 297 DO WHILE ( (rsq.GE.one).OR. (rsq.EQ.zero) ) 291 u1 = REAL(kiss(), 8) / huge64292 u2 = REAL(kiss(), 8) / huge64298 u1 = REAL(kiss(),wp) / huge64 299 u2 = REAL(kiss(),wp) / huge64 293 300 rsq = u1*u1 + u2*u2 294 301 ENDDO … … 316 323 !! -------------------------------------------------------------------- 317 324 IMPLICIT NONE 318 REAL(KIND= 8), PARAMETER :: p1 = 4.5_8319 REAL(KIND= 8), PARAMETER :: p2 = 2.50407739677627_8 ! 1+LOG(9/2)320 REAL(KIND= 8), PARAMETER :: p3 = 1.38629436111989_8 ! LOG(4)321 REAL(KIND= 8) :: gamr, k, u1, u2, b, c, d, xx, yy, zz, rr, ee325 REAL(KIND=wp), PARAMETER :: p1 = 4.5_8 326 REAL(KIND=wp), PARAMETER :: p2 = 2.50407739677627_8 ! 1+LOG(9/2) 327 REAL(KIND=wp), PARAMETER :: p3 = 1.38629436111989_8 ! LOG(4) 328 REAL(KIND=wp) :: gamr, k, u1, u2, b, c, d, xx, yy, zz, rr, ee 322 329 LOGICAL :: accepted 323 330 … … 382 389 !! -------------------------------------------------------------------- 383 390 IMPLICIT NONE 384 INTEGER(KIND= 8), DIMENSION(:) :: a385 INTEGER(KIND= 8) :: n, k, i, j, atmp386 REAL(KIND= 8) :: uran391 INTEGER(KIND=i8), DIMENSION(:) :: a 392 INTEGER(KIND=i8) :: n, k, i, j, atmp 393 REAL(KIND=wp) :: uran 387 394 388 395 ! Select the sample using the swapping method
Note: See TracChangeset
for help on using the changeset viewer.