New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 5404 – NEMO

Changeset 5404


Ignore:
Timestamp:
2015-06-11T09:20:46+02:00 (9 years ago)
Author:
pabouttier
Message:

make the stochastic param. code more compliant with NEMO code conventions compliant; Add documentation skeleton

Location:
trunk
Files:
2 added
3 copied

Legend:

Unmodified
Added
Removed

  • trunk/NEMOGCM/NEMO/OPA_SRC/STO/stopar.F90

    r5390 r5404  
    263263      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    264264 
     265      IF ( lk_agrif ) CALL ctl_stop('EOS stochastic parametrization is not compatible with AGRIF') 
    265266      ! Read namsto namelist : stochastic parameterization 
    266267      REWIND( numnam_ref )              ! Namelist namdyn_adv in reference namelist : Momentum advection scheme 

  • trunk/NEMOGCM/NEMO/OPA_SRC/STO/storng.F90

    r5390 r5404  
    3333   !!   kiss_save, kiss_load, kiss_check, kiss_gamma, kiss_sample 
    3434   !!---------------------------------------------------------------------- 
     35   USE par_kind 
     36   USE lib_mpp 
     37 
    3538   IMPLICIT NONE 
    3639   PRIVATE 
    3740 
    3841   ! Public functions/subroutines 
    39    PUBLIC :: kiss, kiss_seed, kiss_state, kiss_save, kiss_load, kiss_reset, kiss_check 
     42   PUBLIC :: kiss, kiss_seed, kiss_state, kiss_reset ! kiss_save, kiss_load, kiss_check 
    4043   PUBLIC :: kiss_uniform, kiss_gaussian, kiss_gamma, kiss_sample 
    4144 
    4245   ! Default/initial seeds 
    43    INTEGER(KIND=8) :: x=1234567890987654321_8 
    44    INTEGER(KIND=8) :: y=362436362436362436_8 
    45    INTEGER(KIND=8) :: z=1066149217761810_8 
    46    INTEGER(KIND=8) :: w=123456123456123456_8 
     46   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 
    4750 
    4851   ! Parameters to generate real random variates 
    49    REAL(KIND=8), PARAMETER :: huge64=9223372036854775808.0_8  ! +1 
    50    REAL(KIND=8), PARAMETER :: zero=0.0_8, half=0.5_8, one=1.0_8, two=2.0_8 
     52   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 
    5154 
    5255   ! Variables to store 2 Gaussian random numbers with current index (ig) 
    53    INTEGER(KIND=8), SAVE :: ig=1 
    54    REAL(KIND=8), SAVE :: gran1, gran2 
     56   INTEGER(KIND=i8), SAVE :: ig=1 
     57   REAL(KIND=wp), SAVE :: gran1, gran2 
    5558 
    5659   !!---------------------------------------------------------------------- 
    … …  
    7982      ! -------------------------------------------------------------------- 
    8083      IMPLICIT NONE 
    81       INTEGER(KIND=8) :: kiss, t 
     84      INTEGER(KIND=i8) :: kiss, t 
    8285 
    8386      t = ISHFT(x,58) + w 
    … …  
    9699 
    97100         FUNCTION s(k) 
    98             INTEGER(KIND=8) :: s, k 
     101            INTEGER(KIND=i8) :: s, k 
    99102            s = ISHFT(k,-63) 
    100103         END FUNCTION s 
    101104 
    102105         FUNCTION m(k, n) 
    103             INTEGER(KIND=8) :: m, k, n 
     106            INTEGER(KIND=i8) :: m, k, n 
    104107            m =  IEOR(k, ISHFT(k, n) ) 
    105108         END FUNCTION m 
    … …  
    116119      !! -------------------------------------------------------------------- 
    117120      IMPLICIT NONE 
    118       INTEGER(KIND=8) :: ix, iy, iz, iw 
     121      INTEGER(KIND=i8) :: ix, iy, iz, iw 
    119122 
    120123      x = ix 
    … …  
    134137      !! -------------------------------------------------------------------- 
    135138      IMPLICIT NONE 
    136       INTEGER(KIND=8) :: ix, iy, iz, iw 
     139      INTEGER(KIND=i8) :: ix, iy, iz, iw 
    137140 
    138141      ix = x 
    … …  
    161164 
    162165 
    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 
    256263 
    257264 
    … …  
    264271      !! -------------------------------------------------------------------- 
    265272      IMPLICIT NONE 
    266       REAL(KIND=8) :: uran 
    267  
    268       uran = half * ( one + REAL(kiss(),8) / huge64 ) 
     273      REAL(KIND=wp) :: uran 
     274 
     275      uran = half * ( one + REAL(kiss(),wp) / huge64 ) 
    269276 
    270277   END SUBROUTINE kiss_uniform 
    … …  
    284291      !! -------------------------------------------------------------------- 
    285292      IMPLICIT NONE 
    286       REAL(KIND=8) :: gran, u1, u2, rsq, fac 
     293      REAL(KIND=wp) :: gran, u1, u2, rsq, fac 
    287294 
    288295      IF (ig.EQ.1) THEN 
    289296         rsq = two 
    290297         DO WHILE ( (rsq.GE.one).OR. (rsq.EQ.zero) ) 
    291             u1 = REAL(kiss(),8) / huge64 
    292             u2 = REAL(kiss(),8) / huge64 
     298            u1 = REAL(kiss(),wp) / huge64 
     299            u2 = REAL(kiss(),wp) / huge64 
    293300            rsq = u1*u1 + u2*u2 
    294301         ENDDO 
    … …  
    316323      !! -------------------------------------------------------------------- 
    317324      IMPLICIT NONE 
    318       REAL(KIND=8), PARAMETER :: p1 = 4.5_8 
    319       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, ee 
     325      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 
    322329      LOGICAL :: accepted 
    323330 
    … …  
    382389      !! -------------------------------------------------------------------- 
    383390      IMPLICIT NONE 
    384       INTEGER(KIND=8), DIMENSION(:) :: a 
    385       INTEGER(KIND=8) :: n, k, i, j, atmp 
    386       REAL(KIND=8) :: uran 
     391      INTEGER(KIND=i8), DIMENSION(:) :: a 
     392      INTEGER(KIND=i8) :: n, k, i, j, atmp 
     393      REAL(KIND=wp) :: uran 
    387394 
    388395      ! Select the sample using the swapping method 
Note: See TracChangeset for help on using the changeset viewer.