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.
stopar.F90 in trunk/NEMOGCM/NEMO/OPA_SRC – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/stopar.F90 @ 5329

Last change on this file since 5329 was 5329, checked in by pabouttier, 9 years ago

Add stochastic parametrization of EOS
File size: 42.4 KB
Line 
1MODULE stopar
2   !!======================================================================
3   !!                       ***  MODULE  stopar  ***
4   !! Stochastic parameters : definition and time stepping
5   !!=====================================================================
6   !! History :  3.3  ! 2011-10 (J.-M. Brankart)  Original code
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   sto_par       : update the stochastic parameters
11   !!   sto_par_init  : define the stochastic parameterization
12   !!   sto_rst_read  : read restart file for stochastic parameters
13   !!   sto_rst_write : write restart file for stochastic parameters
14   !!   sto_par_white : fill input array with white Gaussian noise
15   !!   sto_par_flt   : apply horizontal Laplacian filter to input array
16   !!----------------------------------------------------------------------
17   USE storng          ! random number generator (external module)
18   USE par_oce         ! ocean parameters
19   USE dom_oce         ! ocean space and time domain variables
20   USE lbclnk          ! lateral boundary conditions (or mpp link)
21   USE in_out_manager  ! I/O manager
22   USE iom             ! I/O module
23   USE lib_mpp
24
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   sto_par_init    ! called by nemogcm.F90
30   PUBLIC   sto_par         ! called by step.F90
31   PUBLIC   sto_rst_write   ! called by step.F90
32
33   LOGICAL           :: ln_rststo = .FALSE.  ! restart stochastic parameters from restart file
34   LOGICAL           :: ln_rstseed = .FALSE. ! read seed of RNG from restart file
35   CHARACTER(len=32) :: cn_storst_in = "restart_sto"     ! suffix of sto restart name (input)
36   CHARACTER(len=32) :: cn_storst_out = "restart_sto"    ! suffix of sto restart name (output)
37   INTEGER           :: numstor, numstow     ! logical unit for restart (read and write)
38
39   INTEGER           :: jpsto2d = 0          ! number of 2D stochastic parameters
40   INTEGER           :: jpsto3d = 0          ! number of 3D stochastic parameters
41
42   REAL(wp), PUBLIC, DIMENSION(:,:,:),   ALLOCATABLE :: sto2d      ! 2D stochastic parameters
43   REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: sto3d      ! 3D stochastic parameters
44   REAL(wp),         DIMENSION(:,:),     ALLOCATABLE :: sto_tmp    ! temporary workspace
45   REAL(wp),         DIMENSION(:,:),     ALLOCATABLE :: sto2d_abc  ! a, b, c parameters (for 2D arrays)
46   REAL(wp),         DIMENSION(:,:),     ALLOCATABLE :: sto3d_abc  ! a, b, c parameters (for 3D arrays)
47   REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto2d_ave  ! mean value (for 2D arrays)
48   REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto3d_ave  ! mean value (for 3D arrays)
49   REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto2d_std  ! standard deviation (for 2D arrays)
50   REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto3d_std  ! standard deviation (for 3D arrays)
51   REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto2d_lim  ! limitation factor (for 2D arrays)
52   REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto3d_lim  ! limitation factor (for 3D arrays)
53   REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto2d_tcor ! time correlation (for 2D arrays)
54   REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto3d_tcor ! time correlation (for 3D arrays)
55   INTEGER,          DIMENSION(:),       ALLOCATABLE :: sto2d_ord  ! order of autoregressive process
56   INTEGER,          DIMENSION(:),       ALLOCATABLE :: sto3d_ord  ! order of autoregressive process
57
58   CHARACTER(len=1), DIMENSION(:),       ALLOCATABLE :: sto2d_typ  ! nature of grid point (T, U, V, W, F, I)
59   CHARACTER(len=1), DIMENSION(:),       ALLOCATABLE :: sto3d_typ  ! nature of grid point (T, U, V, W, F, I)
60   REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto2d_sgn  ! control of the sign accross the north fold
61   REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto3d_sgn  ! control of the sign accross the north fold
62   INTEGER,          DIMENSION(:),       ALLOCATABLE :: sto2d_flt  ! number of passes of Laplacian filter
63   INTEGER,          DIMENSION(:),       ALLOCATABLE :: sto3d_flt  ! number of passes of Laplacian filter
64   REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto2d_fac  ! factor to restore std after filtering
65   REAL(wp),         DIMENSION(:),       ALLOCATABLE :: sto3d_fac  ! factor to restore std after filtering
66
67   LOGICAL, PUBLIC :: ln_sto_ldf = .FALSE.    ! stochastic lateral diffusion
68   INTEGER, PUBLIC :: jsto_ldf                ! index of lateral diffusion stochastic parameter
69   REAL(wp)        :: rn_ldf_std              ! lateral diffusion standard deviation (in percent)
70   REAL(wp)        :: rn_ldf_tcor             ! lateral diffusion correlation timescale (in timesteps)
71
72   LOGICAL, PUBLIC :: ln_sto_hpg = .FALSE.    ! stochastic horizontal pressure gradient
73   INTEGER, PUBLIC :: jsto_hpgi               ! index of stochastic hpg parameter (i direction)
74   INTEGER, PUBLIC :: jsto_hpgj               ! index of stochastic hpg parameter (j direction)
75   REAL(wp)        :: rn_hpg_std              ! density gradient standard deviation (in percent)
76   REAL(wp)        :: rn_hpg_tcor             ! density gradient correlation timescale (in timesteps)
77
78   LOGICAL, PUBLIC :: ln_sto_pstar = .FALSE.  ! stochastic ice strength
79   INTEGER, PUBLIC :: jsto_pstar              ! index of stochastic ice strength
80   REAL(wp), PUBLIC:: rn_pstar_std            ! ice strength standard deviation (in percent)
81   REAL(wp)        :: rn_pstar_tcor           ! ice strength correlation timescale (in timesteps)
82   INTEGER         :: nn_pstar_flt = 0        ! number of passes of Laplacian filter
83   INTEGER         :: nn_pstar_ord = 1        ! order of autoregressive processes
84
85   LOGICAL, PUBLIC :: ln_sto_trd = .FALSE.    ! stochastic model trend
86   INTEGER, PUBLIC :: jsto_trd                ! index of stochastic trend parameter
87   REAL(wp)        :: rn_trd_std              ! trend standard deviation (in percent)
88   REAL(wp)        :: rn_trd_tcor             ! trend correlation timescale (in timesteps)
89
90   LOGICAL, PUBLIC :: ln_sto_eos = .FALSE.    ! stochastic equation of state
91   INTEGER, PUBLIC :: nn_sto_eos = 1          ! number of degrees of freedom in stochastic equation of state
92   INTEGER, PUBLIC, DIMENSION(:), ALLOCATABLE :: jsto_eosi ! index of stochastic eos parameter (i direction)
93   INTEGER, PUBLIC, DIMENSION(:), ALLOCATABLE :: jsto_eosj ! index of stochastic eos parameter (j direction)
94   INTEGER, PUBLIC, DIMENSION(:), ALLOCATABLE :: jsto_eosk ! index of stochastic eos parameter (k direction)
95   REAL(wp)        :: rn_eos_stdxy            ! random walk horz. standard deviation (in grid points)
96   REAL(wp)        :: rn_eos_stdz             ! random walk vert. standard deviation (in grid points)
97   REAL(wp)        :: rn_eos_tcor             ! random walk correlation timescale (in timesteps)
98   REAL(wp)        :: rn_eos_lim = 3.0_wp     ! limitation factor
99   INTEGER         :: nn_eos_flt = 0          ! number of passes of Laplacian filter
100   INTEGER         :: nn_eos_ord = 1          ! order of autoregressive processes
101
102   LOGICAL, PUBLIC :: ln_sto_trc = .FALSE.    ! stochastic tracer dynamics
103   INTEGER, PUBLIC :: nn_sto_trc = 1          ! number of degrees of freedom in stochastic tracer dynamics
104   INTEGER, PUBLIC, DIMENSION(:), ALLOCATABLE :: jsto_trci ! index of stochastic trc parameter (i direction)
105   INTEGER, PUBLIC, DIMENSION(:), ALLOCATABLE :: jsto_trcj ! index of stochastic trc parameter (j direction)
106   INTEGER, PUBLIC, DIMENSION(:), ALLOCATABLE :: jsto_trck ! index of stochastic trc parameter (k direction)
107   REAL(wp)        :: rn_trc_stdxy            ! random walk horz. standard deviation (in grid points)
108   REAL(wp)        :: rn_trc_stdz             ! random walk vert. standard deviation (in grid points)
109   REAL(wp)        :: rn_trc_tcor             ! random walk correlation timescale (in timesteps)
110   REAL(wp)        :: rn_trc_lim = 3.0_wp     ! limitation factor
111   INTEGER         :: nn_trc_flt = 0          ! number of passes of Laplacian filter
112   INTEGER         :: nn_trc_ord = 1          ! order of autoregressive processes
113
114   ! Public array with density correction
115   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: drho_ran
116
117   !!----------------------------------------------------------------------
118   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
119   !! $Id: dynhpg.F90 2528 2010-12-27 17:33:53Z rblod $
120   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
121   !!----------------------------------------------------------------------
122CONTAINS
123
124   SUBROUTINE sto_par( kt )
125      !!----------------------------------------------------------------------
126      !!                  ***  ROUTINE sto_par  ***
127      !!
128      !! ** Purpose :   update the stochastic parameters
129      !!
130      !! ** Method  :   model basic stochastic parameters
131      !!                as a first order autoregressive process AR(1),
132      !!                governed by the equation:
133      !!                   X(t) = a * X(t-1) + b * w + c
134      !!                where the parameters a, b and c are related
135      !!                to expected value, standard deviation
136      !!                and time correlation (all stationary in time) by:
137      !!                   E   [X(t)]        = c / ( 1 - a )
138      !!                   STD [X(t)]        = b / SQRT( 1 - a * a )
139      !!                   COR [X(t),X(t-k)] = a ** k
140      !!                and w is a Gaussian white noise.
141      !!
142      !!                Higher order autoregressive proces can be optionally generated
143      !!                by replacing the white noise by a lower order process.
144      !!
145      !!                1) The statistics of the stochastic parameters (X) are assumed
146      !!                constant in space (homogeneous) and time (stationary).
147      !!                This could be generalized by replacing the constant
148      !!                a, b, c parameters by functions of space and time.
149      !!
150      !!                2) The computation is performed independently for every model
151      !!                grid point, which corresponds to assume that the stochastic
152      !!                parameters are uncorrelated in space.
153      !!                This could be generalized by including a spatial filter: Y = Filt[ X ]
154      !!                (possibly non-homgeneous and non-stationary) in the computation,
155      !!                or by solving an elliptic equation: L[ Y ] = X.
156      !!
157      !!                3) The stochastic model for the parameters could also
158      !!                be generalized to depend on the current state of the ocean (not done here).
159      !!----------------------------------------------------------------------
160      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
161      !!
162      INTEGER  :: ji, jj, jk, jsto, jflt
163      REAL(wp) :: stomax
164
165      !
166      ! Update 2D stochastic arrays
167      !
168      DO jsto = 1, jpsto2d
169        ! Store array from previous time step
170        sto_tmp(:,:) = sto2d(:,:,jsto)
171
172        IF ( sto2d_ord(jsto) == 1 ) THEN
173          ! Draw new random numbers from N(0,1) --> w
174          CALL sto_par_white( sto2d(:,:,jsto) )
175          ! Apply horizontal Laplacian filter to w
176          DO jflt = 1, sto2d_flt(jsto)
177            CALL lbc_lnk( sto2d(:,:,jsto), sto2d_typ(jsto), sto2d_sgn(jsto) )
178            CALL sto_par_flt( sto2d(:,:,jsto) )
179          END DO
180          ! Factor to restore standard deviation after filtering
181          sto2d(:,:,jsto) = sto2d(:,:,jsto) * sto2d_fac(jsto)
182        ELSE
183          ! Use previous process (one order lower) instead of white noise
184          sto2d(:,:,jsto) = sto2d(:,:,jsto-1)
185        ENDIF
186
187        ! Multiply white noise (or lower order process) by b --> b * w
188        sto2d(:,:,jsto) = sto2d(:,:,jsto) * sto2d_abc(jsto,2)
189        ! Update autoregressive processes --> a * X(t-1) + b * w
190        sto2d(:,:,jsto) = sto2d(:,:,jsto) + sto_tmp(:,:) * sto2d_abc(jsto,1)
191        ! Add parameter c --> a * X(t-1) + b * w + c
192        sto2d(:,:,jsto) = sto2d(:,:,jsto) + sto2d_abc(jsto,3)
193        ! Limit random parameter anomalies to std times the limitation factor
194        stomax = sto2d_std(jsto) * sto2d_lim(jsto)
195        sto2d(:,:,jsto) = sto2d(:,:,jsto) - sto2d_ave(jsto)
196        sto2d(:,:,jsto) = SIGN(MIN(stomax,ABS(sto2d(:,:,jsto))),sto2d(:,:,jsto))
197        sto2d(:,:,jsto) = sto2d(:,:,jsto) + sto2d_ave(jsto)
198
199        ! Lateral boundary conditions on sto2d
200        CALL lbc_lnk( sto2d(:,:,jsto), sto2d_typ(jsto), sto2d_sgn(jsto) )
201      END DO
202      !
203      ! Update 3D stochastic arrays
204      !
205      DO jsto = 1, jpsto3d
206         DO jk = 1, jpk
207           ! Store array from previous time step
208           sto_tmp(:,:) = sto3d(:,:,jk,jsto)
209
210           IF ( sto3d_ord(jsto) == 1 ) THEN
211             ! Draw new random numbers from N(0,1) --> w
212             CALL sto_par_white( sto3d(:,:,jk,jsto) )
213             ! Apply horizontal Laplacian filter to w
214             DO jflt = 1, sto3d_flt(jsto)
215               CALL lbc_lnk( sto3d(:,:,jk,jsto), sto3d_typ(jsto), sto3d_sgn(jsto) )
216               CALL sto_par_flt( sto3d(:,:,jk,jsto) )
217             END DO
218             ! Factor to restore standard deviation after filtering
219             sto3d(:,:,jk,jsto) = sto3d(:,:,jk,jsto) * sto3d_fac(jsto)
220           ELSE
221             ! Use previous process (one order lower) instead of white noise
222             sto3d(:,:,jk,jsto) = sto3d(:,:,jk,jsto-1)
223           ENDIF
224
225           ! Multiply white noise by b --> b * w
226           sto3d(:,:,jk,jsto) = sto3d(:,:,jk,jsto) * sto3d_abc(jsto,2)
227           ! Update autoregressive processes --> a * X(t-1) + b * w
228           sto3d(:,:,jk,jsto) = sto3d(:,:,jk,jsto) + sto_tmp(:,:) * sto3d_abc(jsto,1)
229           ! Add parameter c --> a * X(t-1) + b * w + c
230           sto3d(:,:,jk,jsto) = sto3d(:,:,jk,jsto) + sto3d_abc(jsto,3)
231           ! Limit random parameters anomalies to std times the limitation factor
232           stomax = sto3d_std(jsto) * sto3d_lim(jsto)
233           sto3d(:,:,jk,jsto) = sto3d(:,:,jk,jsto) - sto3d_ave(jsto)
234           sto3d(:,:,jk,jsto) = SIGN(MIN(stomax,ABS(sto3d(:,:,jk,jsto))),sto3d(:,:,jk,jsto))
235           sto3d(:,:,jk,jsto) = sto3d(:,:,jk,jsto) + sto3d_ave(jsto)
236         END DO
237         ! Lateral boundary conditions on sto3d
238         CALL lbc_lnk( sto3d(:,:,:,jsto), sto3d_typ(jsto), sto3d_sgn(jsto) )
239      END DO
240
241   END SUBROUTINE sto_par
242
243
244   SUBROUTINE sto_par_init
245      !!----------------------------------------------------------------------
246      !!                  ***  ROUTINE sto_par_init  ***
247      !!
248      !! ** Purpose :   define the stochastic parameterization
249      !!----------------------------------------------------------------------
250      NAMELIST/namsto/ ln_sto_ldf, rn_ldf_std, rn_ldf_tcor, &
251        &              ln_sto_hpg, rn_hpg_std, rn_hpg_tcor, &
252        &              ln_sto_pstar, rn_pstar_std, rn_pstar_tcor, nn_pstar_flt, nn_pstar_ord, &
253        &              ln_sto_trd, rn_trd_std, rn_trd_tcor, &
254        &              ln_sto_eos, nn_sto_eos, rn_eos_stdxy, rn_eos_stdz, &
255        &              rn_eos_tcor, nn_eos_ord, nn_eos_flt, rn_eos_lim, &
256        &              ln_sto_trc, nn_sto_trc, rn_trc_stdxy, rn_trc_stdz, &
257        &              rn_trc_tcor, nn_trc_ord, nn_trc_flt, rn_trc_lim, &
258        &              ln_rststo, ln_rstseed, cn_storst_in, cn_storst_out
259      !!----------------------------------------------------------------------
260      INTEGER :: jsto, jmem, jarea, jdof, jord, jordm1, jk, jflt
261      INTEGER(KIND=8) :: zseed1, zseed2, zseed3, zseed4
262      REAL(wp) :: rinflate
263      INTEGER  ::   ios                 ! Local integer output status for namelist read
264
265      ! Read namsto namelist : stochastic parameterization
266      REWIND( numnam_ref )              ! Namelist namdyn_adv in reference namelist : Momentum advection scheme
267      READ  ( numnam_ref, namsto, IOSTAT = ios, ERR = 901)
268901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsto in reference namelist', lwp )
269
270      REWIND( numnam_cfg )              ! Namelist namdyn_adv in configuration namelist : Momentum advection scheme
271      READ  ( numnam_cfg, namsto, IOSTAT = ios, ERR = 902 )
272902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsto in configuration namelist', lwp )
273      IF(lwm) WRITE ( numond, namsto )
274
275      !IF(ln_ens_rst_in) cn_storst_in = cn_mem//cn_storst_in
276
277      ! Parameter print
278      IF(lwp) THEN
279         WRITE(numout,*)
280         WRITE(numout,*) 'sto_par_init : stochastic parameterization'
281         WRITE(numout,*) '~~~~~~~~~~~~'
282         WRITE(numout,*) '   Namelist namsto : stochastic parameterization'
283         WRITE(numout,*) '      restart stochastic parameters           ln_rststo     = ', ln_rststo
284         WRITE(numout,*) '      read seed of RNG from restart file      ln_rstseed    = ', ln_rstseed
285         WRITE(numout,*) '      suffix of sto restart name (input)      cn_storst_in  = ', cn_storst_in
286         WRITE(numout,*) '      suffix of sto restart name (output)     cn_storst_out = ', cn_storst_out
287
288         ! WRITE(numout,*) '      stochastic lateral diffusion            ln_sto_ldf    = ', ln_sto_ldf
289         ! WRITE(numout,*) '      lateral diffusion std (in percent)      rn_ldf_std    = ', rn_ldf_std
290         ! WRITE(numout,*) '      lateral diffusion tcor (in timesteps)   rn_ldf_tcor   = ', rn_ldf_tcor
291
292         ! WRITE(numout,*) '      stochastic horizontal pressure gradient ln_sto_hpg    = ', ln_sto_hpg
293         ! WRITE(numout,*) '      density gradient std (in percent)       rn_hpg_std    = ', rn_hpg_std
294         ! WRITE(numout,*) '      density gradient tcor (in timesteps)    rn_hpg_tcor   = ', rn_hpg_tcor
295
296         ! WRITE(numout,*) '      stochastic ice strength                 ln_sto_pstar  = ', ln_sto_pstar
297         ! WRITE(numout,*) '      ice strength std (in percent)           rn_pstar_std  = ', rn_pstar_std
298         ! WRITE(numout,*) '      ice strength tcor (in timesteps)        rn_pstar_tcor = ', rn_pstar_tcor
299         ! WRITE(numout,*) '      order of autoregressive  processes      nn_pstar_ord  = ', nn_pstar_ord
300         ! WRITE(numout,*) '      passes of Laplacian filter              nn_pstar_flt  = ', nn_pstar_flt
301
302         !WRITE(numout,*) '      stochastic trend                        ln_sto_trd    = ', ln_sto_trd
303         !WRITE(numout,*) '      trend std (in percent)                  rn_trd_std    = ', rn_trd_std
304         !WRITE(numout,*) '      trend tcor (in timesteps)               rn_trd_tcor   = ', rn_trd_tcor
305
306         WRITE(numout,*) '      stochastic equation of state            ln_sto_eos    = ', ln_sto_eos
307         WRITE(numout,*) '      number of degrees of freedom            nn_sto_eos    = ', nn_sto_eos
308         WRITE(numout,*) '      random walk horz. std (in grid points)  rn_eos_stdxy  = ', rn_eos_stdxy
309         WRITE(numout,*) '      random walk vert. std (in grid points)  rn_eos_stdz   = ', rn_eos_stdz
310         WRITE(numout,*) '      random walk tcor (in timesteps)         rn_eos_tcor   = ', rn_eos_tcor
311         WRITE(numout,*) '      order of autoregressive  processes      nn_eos_ord    = ', nn_eos_ord
312         WRITE(numout,*) '      passes of Laplacian filter              nn_eos_flt    = ', nn_eos_flt
313         WRITE(numout,*) '      limitation factor                       rn_eos_lim    = ', rn_eos_lim
314
315         ! WRITE(numout,*) '      stochastic tracers dynamics             ln_sto_trc    = ', ln_sto_trc
316         ! WRITE(numout,*) '      number of degrees of freedom            nn_sto_trc    = ', nn_sto_trc
317         ! WRITE(numout,*) '      random walk horz. std (in grid points)  rn_trc_stdxy  = ', rn_trc_stdxy
318         ! WRITE(numout,*) '      random walk vert. std (in grid points)  rn_trc_stdz   = ', rn_trc_stdz
319         ! WRITE(numout,*) '      random walk tcor (in timesteps)         rn_trc_tcor   = ', rn_trc_tcor
320         ! WRITE(numout,*) '      order of autoregressive  processes      nn_trc_ord    = ', nn_trc_ord
321         ! WRITE(numout,*) '      passes of Laplacian filter              nn_trc_flt    = ', nn_trc_flt
322         ! WRITE(numout,*) '      limitation factor                       rn_trc_lim    = ', rn_trc_lim
323
324      ENDIF
325
326      IF(lwp) WRITE(numout,*)
327      IF(lwp) WRITE(numout,*) '   stochastic parameterization :'
328
329      ! Set number of 2D stochastic arrays
330      jpsto2d = 0
331      IF( ln_sto_ldf ) THEN
332         IF(lwp) WRITE(numout,*) '       - stochastic lateral diffusion'
333         jpsto2d   = jpsto2d + 1
334         jsto_ldf  = jpsto2d
335      ENDIF
336      IF( ln_sto_pstar ) THEN
337         IF(lwp) WRITE(numout,*) '       - stochastic ice strength'
338         jpsto2d    = jpsto2d + 1 * nn_pstar_ord
339         jsto_pstar = jpsto2d
340      ENDIF
341      IF( ln_sto_eos ) THEN
342         IF(lwp) WRITE(numout,*) '       - stochastic equation of state'
343         ALLOCATE(jsto_eosi(nn_sto_eos))
344         ALLOCATE(jsto_eosj(nn_sto_eos))
345         ALLOCATE(jsto_eosk(nn_sto_eos))
346         DO jdof = 1, nn_sto_eos
347            jpsto2d   = jpsto2d + 3 * nn_eos_ord
348            jsto_eosi(jdof) = jpsto2d - 2 * nn_eos_ord
349            jsto_eosj(jdof) = jpsto2d - 1 * nn_eos_ord
350            jsto_eosk(jdof) = jpsto2d
351         END DO
352      ELSE
353         nn_sto_eos = 0
354      ENDIF
355      IF( ln_sto_trc ) THEN
356         IF(lwp) WRITE(numout,*) '       - stochastic tracers dynamics'
357         ALLOCATE(jsto_trci(nn_sto_trc))
358         ALLOCATE(jsto_trcj(nn_sto_trc))
359         ALLOCATE(jsto_trck(nn_sto_trc))
360         DO jdof = 1, nn_sto_trc
361            jpsto2d   = jpsto2d + 3 * nn_trc_ord
362            jsto_trci(jdof) = jpsto2d - 2 * nn_trc_ord
363            jsto_trcj(jdof) = jpsto2d - 1 * nn_trc_ord
364            jsto_trck(jdof) = jpsto2d
365         END DO
366      ELSE
367         nn_sto_trc = 0
368      ENDIF
369
370      ! Set number of 3D stochastic arrays
371      jpsto3d = 0
372      IF( ln_sto_hpg ) THEN
373         IF(lwp) WRITE(numout,*) '       - stochastic horizontal pressure gradient'
374         jpsto3d   = jpsto3d + 2
375         jsto_hpgi = jpsto3d - 1
376         jsto_hpgj = jpsto3d
377      ENDIF
378      IF( ln_sto_trd ) THEN
379         IF(lwp) WRITE(numout,*) '       - stochastic trend'
380         jpsto3d   = jpsto3d + 1
381         jsto_trd  = jpsto3d
382      ENDIF
383
384      ! Allocate 2D stochastic arrays
385      IF ( jpsto2d > 0 ) THEN
386         ALLOCATE ( sto2d(jpi,jpj,jpsto2d) )
387         ALLOCATE ( sto2d_abc(jpsto2d,3) )
388         ALLOCATE ( sto2d_ave(jpsto2d) )
389         ALLOCATE ( sto2d_std(jpsto2d) )
390         ALLOCATE ( sto2d_lim(jpsto2d) )
391         ALLOCATE ( sto2d_tcor(jpsto2d) )
392         ALLOCATE ( sto2d_ord(jpsto2d) )
393         ALLOCATE ( sto2d_typ(jpsto2d) )
394         ALLOCATE ( sto2d_sgn(jpsto2d) )
395         ALLOCATE ( sto2d_flt(jpsto2d) )
396         ALLOCATE ( sto2d_fac(jpsto2d) )
397      ENDIF
398
399      ! Allocate 3D stochastic arrays
400      IF ( jpsto3d > 0 ) THEN
401         ALLOCATE ( sto3d(jpi,jpj,jpk,jpsto3d) )
402         ALLOCATE ( sto3d_abc(jpsto3d,3) )
403         ALLOCATE ( sto3d_ave(jpsto3d) )
404         ALLOCATE ( sto3d_std(jpsto3d) )
405         ALLOCATE ( sto3d_lim(jpsto3d) )
406         ALLOCATE ( sto3d_tcor(jpsto3d) )
407         ALLOCATE ( sto3d_ord(jpsto3d) )
408         ALLOCATE ( sto3d_typ(jpsto3d) )
409         ALLOCATE ( sto3d_sgn(jpsto3d) )
410         ALLOCATE ( sto3d_flt(jpsto3d) )
411         ALLOCATE ( sto3d_fac(jpsto3d) )
412      ENDIF
413
414      ! Allocate temporary workspace
415      IF ( jpsto2d > 0 .OR. jpsto3d > 0 ) THEN
416         ALLOCATE ( sto_tmp(jpi,jpj) ) ; sto_tmp(:,:) = 0._wp
417      ENDIF
418
419      ! 1) For every stochastic parameter:
420      ! ----------------------------------
421      ! - set nature of grid point and control of the sign
422      !       across the north fold (sto2d_typ, sto2d_sgn)
423      ! - set number of passes of Laplacian filter (sto2d_flt)
424      ! - set order of every autoregressive process (sto2d_ord)
425      DO jsto = 1, jpsto2d
426         sto2d_typ(jsto) = 'T'
427         sto2d_sgn(jsto) = 1._wp
428         sto2d_flt(jsto) = 0
429         sto2d_ord(jsto) = 1
430         DO jord = 0, nn_pstar_ord-1
431            IF ( jsto+jord == jsto_pstar ) THEN ! Stochastic ice strength (ave=1)
432               sto2d_ord(jsto) = nn_pstar_ord - jord
433               sto2d_flt(jsto) = nn_pstar_flt
434            ENDIF
435         ENDDO
436         DO jdof = 1, nn_sto_eos
437         DO jord = 0, nn_eos_ord-1
438            IF ( jsto+jord == jsto_eosi(jdof) ) THEN ! Stochastic equation of state i (ave=0)
439               sto2d_ord(jsto) = nn_eos_ord - jord
440               sto2d_sgn(jsto) = -1._wp
441               sto2d_flt(jsto) = nn_eos_flt
442            ENDIF
443            IF ( jsto+jord == jsto_eosj(jdof) ) THEN ! Stochastic equation of state j (ave=0)
444               sto2d_ord(jsto) = nn_eos_ord - jord
445               sto2d_sgn(jsto) = -1._wp
446               sto2d_flt(jsto) = nn_eos_flt
447            ENDIF
448            IF ( jsto+jord == jsto_eosk(jdof) ) THEN ! Stochastic equation of state k (ave=0)
449               sto2d_ord(jsto) = nn_eos_ord - jord
450               sto2d_flt(jsto) = nn_eos_flt
451            ENDIF
452         END DO
453         END DO
454         DO jdof = 1, nn_sto_trc
455         DO jord = 0, nn_trc_ord-1
456            IF ( jsto+jord == jsto_trci(jdof) ) THEN ! Stochastic tracers dynamics i (ave=0)
457               sto2d_ord(jsto) = nn_trc_ord - jord
458               sto2d_sgn(jsto) = -1._wp
459               sto2d_flt(jsto) = nn_trc_flt
460            ENDIF
461            IF ( jsto+jord == jsto_trcj(jdof) ) THEN ! Stochastic tracers dynamics j (ave=0)
462               sto2d_ord(jsto) = nn_trc_ord - jord
463               sto2d_sgn(jsto) = -1._wp
464               sto2d_flt(jsto) = nn_trc_flt
465            ENDIF
466            IF ( jsto+jord == jsto_trck(jdof) ) THEN ! Stochastic tracers dynamics k (ave=0)
467               sto2d_ord(jsto) = nn_trc_ord - jord
468               sto2d_flt(jsto) = nn_trc_flt
469            ENDIF
470         END DO
471         END DO
472
473         sto2d_fac(jsto) = sto_par_flt_fac ( sto2d_flt(jsto) )
474      END DO
475      !
476      DO jsto = 1, jpsto3d
477         sto3d_typ(jsto) = 'T'
478         sto3d_sgn(jsto) = 1._wp
479         sto3d_flt(jsto) = 0
480         sto3d_ord(jsto) = 1
481         IF ( jsto == jsto_hpgi ) THEN ! Stochastic density gradient i (ave=1)
482            sto3d_typ(jsto) = 'U'
483         ENDIF
484         IF ( jsto == jsto_hpgj ) THEN ! Stochastic density gradient j (ave=1)
485            sto3d_typ(jsto) = 'V'
486         ENDIF
487         sto3d_fac(jsto) = sto_par_flt_fac ( sto3d_flt(jsto) )
488      END DO
489
490      ! 2) For every stochastic parameter:
491      ! ----------------------------------
492      ! set average, standard deviation and time correlation
493      DO jsto = 1, jpsto2d
494         sto2d_ave(jsto)  = 0._wp
495         sto2d_std(jsto)  = 1._wp
496         sto2d_tcor(jsto) = 1._wp
497         sto2d_lim(jsto)  = 3._wp
498         IF ( jsto == jsto_ldf  ) THEN ! Stochastic lateral diffusion (ave=1)
499            sto2d_ave(jsto)  = 1._wp
500            sto2d_std(jsto)  = rn_ldf_std
501            sto2d_tcor(jsto) = rn_ldf_tcor
502         ENDIF
503         DO jord = 0, nn_pstar_ord-1
504            IF ( jsto+jord == jsto_pstar ) THEN ! Stochastic ice strength (ave=1)
505               sto2d_std(jsto) = 1._wp
506               sto2d_tcor(jsto) = rn_pstar_tcor
507            ENDIF
508         ENDDO
509         DO jdof = 1, nn_sto_eos
510         DO jord = 0, nn_eos_ord-1
511            IF ( jsto+jord == jsto_eosi(jdof) ) THEN ! Stochastic equation of state i (ave=0)
512               sto2d_std(jsto)  = rn_eos_stdxy
513               sto2d_tcor(jsto) = rn_eos_tcor
514               sto2d_lim(jsto)  = rn_eos_lim
515            ENDIF
516            IF ( jsto+jord == jsto_eosj(jdof) ) THEN ! Stochastic equation of state j (ave=0)
517               sto2d_std(jsto)  = rn_eos_stdxy
518               sto2d_tcor(jsto) = rn_eos_tcor
519               sto2d_lim(jsto)  = rn_eos_lim
520            ENDIF
521            IF ( jsto+jord == jsto_eosk(jdof) ) THEN ! Stochastic equation of state k (ave=0)
522               sto2d_std(jsto)  = rn_eos_stdz
523               sto2d_tcor(jsto) = rn_eos_tcor
524               sto2d_lim(jsto)  = rn_eos_lim
525            ENDIF
526         END DO
527         END DO
528         DO jdof = 1, nn_sto_trc
529         DO jord = 0, nn_trc_ord-1
530            IF ( jsto+jord == jsto_trci(jdof) ) THEN ! Stochastic tracer dynamics i (ave=0)
531               sto2d_std(jsto)  = rn_trc_stdxy
532               sto2d_tcor(jsto) = rn_trc_tcor
533               sto2d_lim(jsto)  = rn_trc_lim
534            ENDIF
535            IF ( jsto+jord == jsto_trcj(jdof) ) THEN ! Stochastic tracer dynamics j (ave=0)
536               sto2d_std(jsto)  = rn_trc_stdxy
537               sto2d_tcor(jsto) = rn_trc_tcor
538               sto2d_lim(jsto)  = rn_trc_lim
539            ENDIF
540            IF ( jsto+jord == jsto_trck(jdof) ) THEN ! Stochastic tracer dynamics k (ave=0)
541               sto2d_std(jsto)  = rn_trc_stdz
542               sto2d_tcor(jsto) = rn_trc_tcor
543               sto2d_lim(jsto)  = rn_trc_lim
544            ENDIF
545         END DO
546         END DO
547
548      END DO
549      !
550      DO jsto = 1, jpsto3d
551         sto3d_ave(jsto)  = 0._wp
552         sto3d_std(jsto)  = 1._wp
553         sto3d_tcor(jsto) = 1._wp
554         sto3d_lim(jsto)  = 3._wp
555         IF ( jsto == jsto_hpgi ) THEN ! Stochastic density gradient i (ave=1)
556            sto3d_ave(jsto)  = 1._wp
557            sto3d_std(jsto)  = rn_hpg_std
558            sto3d_tcor(jsto) = rn_hpg_tcor
559         ENDIF
560         IF ( jsto == jsto_hpgj ) THEN ! Stochastic density gradient j (ave=1)
561            sto3d_ave(jsto)  = 1._wp
562            sto3d_std(jsto)  = rn_hpg_std
563            sto3d_tcor(jsto) = rn_hpg_tcor
564         ENDIF
565         IF ( jsto == jsto_trd ) THEN ! Stochastic trend (ave=1)
566            sto3d_ave(jsto)  = 1._wp
567            sto3d_std(jsto)  = rn_trd_std
568            sto3d_tcor(jsto) = rn_trd_tcor
569         ENDIF
570      END DO
571
572      ! 3) For every stochastic parameter:
573      ! ----------------------------------
574      ! - compute parameters (a, b, c) of the AR1 autoregressive process
575      !   from expected value (ave), standard deviation (std)
576      !   and time correlation (tcor):
577      !     a = EXP ( - 1 / tcor )           --> sto2d_abc(:,1)
578      !     b = std * SQRT( 1 - a * a )      --> sto2d_abc(:,2)
579      !     c = ave * ( 1 - a )              --> sto2d_abc(:,3)
580      ! - for higher order processes (ARn, n>1), use approximate formula
581      !   for the b parameter (valid for tcor>>1 time step)
582      DO jsto = 1, jpsto2d
583         IF ( sto2d_tcor(jsto) == 0._wp ) THEN
584            sto2d_abc(jsto,1) = 0._wp
585         ELSE
586            sto2d_abc(jsto,1) = EXP ( - 1._wp / sto2d_tcor(jsto) )
587         ENDIF
588         IF ( sto2d_ord(jsto) == 1 ) THEN      ! Exact formula for 1st order process
589            rinflate = sto2d_std(jsto)
590         ELSE
591            ! Approximate formula, valid for tcor >> 1
592            jordm1 = sto2d_ord(jsto) - 1
593            rinflate = SQRT ( REAL( jordm1 , wp ) / REAL( 2*(2*jordm1-1) , wp ) )
594         ENDIF
595         sto2d_abc(jsto,2) = rinflate * SQRT ( 1._wp - sto2d_abc(jsto,1) &
596                                                     * sto2d_abc(jsto,1) )
597         sto2d_abc(jsto,3) = sto2d_ave(jsto) * ( 1._wp - sto2d_abc(jsto,1) )
598      END DO
599      !
600      DO jsto = 1, jpsto3d
601         IF ( sto3d_tcor(jsto) == 0._wp ) THEN
602            sto3d_abc(jsto,1) = 0._wp
603         ELSE
604            sto3d_abc(jsto,1) = EXP ( - 1._wp / sto3d_tcor(jsto) )
605         ENDIF
606         IF ( sto3d_ord(jsto) == 1 ) THEN      ! Exact formula for 1st order process
607            rinflate = sto3d_std(jsto)
608         ELSE
609            ! Approximate formula, valid for tcor >> 1
610            jordm1 = sto3d_ord(jsto) - 1
611            rinflate = SQRT ( REAL( jordm1 , wp ) / REAL( 2*(2*jordm1-1) , wp ) )
612         ENDIF
613         sto3d_abc(jsto,2) = rinflate * SQRT ( 1._wp - sto3d_abc(jsto,1) &
614                                                     * sto3d_abc(jsto,1) )
615         sto3d_abc(jsto,3) = sto3d_ave(jsto) * ( 1._wp - sto3d_abc(jsto,1) )
616      END DO
617
618      ! 4) Initialize seeds for random number generator
619      ! -----------------------------------------------
620      ! using different seeds for different processors (jarea)
621      ! and different ensemble members (jmem)
622      CALL kiss_reset( )
623      DO jarea = 1, narea
624         !DO jmem = 0, nmember
625         zseed1 = kiss() ; zseed2 = kiss() ; zseed3 = kiss() ; zseed4 = kiss()
626         !END DO
627      END DO
628      CALL kiss_seed( zseed1, zseed2, zseed3, zseed4 )
629
630      ! 5) Initialize stochastic parameters to: ave + std * w
631      ! -----------------------------------------------------
632      DO jsto = 1, jpsto2d
633         ! Draw random numbers from N(0,1) --> w
634         CALL sto_par_white( sto2d(:,:,jsto) )
635         ! Apply horizontal Laplacian filter to w
636         DO jflt = 1, sto2d_flt(jsto)
637            CALL lbc_lnk( sto2d(:,:,jsto), sto2d_typ(jsto), sto2d_sgn(jsto) )
638            CALL sto_par_flt( sto2d(:,:,jsto) )
639         END DO
640         ! Factor to restore standard deviation after filtering
641         sto2d(:,:,jsto) = sto2d(:,:,jsto) * sto2d_fac(jsto)
642         ! Limit random parameter to the limitation factor
643         sto2d(:,:,jsto) = SIGN(MIN(sto2d_lim(jsto),ABS(sto2d(:,:,jsto))),sto2d(:,:,jsto))
644         ! Multiply by standard devation and add average value
645         sto2d(:,:,jsto) = sto2d(:,:,jsto) * sto2d_std(jsto) + sto2d_ave(jsto)
646      END DO
647      !
648      DO jsto = 1, jpsto3d
649         DO jk = 1, jpk
650            ! Draw random numbers from N(0,1) --> w
651            CALL sto_par_white( sto3d(:,:,jk,jsto) )
652            ! Apply horizontal Laplacian filter to w
653            DO jflt = 1, sto3d_flt(jsto)
654               CALL lbc_lnk( sto3d(:,:,jk,jsto), sto3d_typ(jsto), sto3d_sgn(jsto) )
655               CALL sto_par_flt( sto3d(:,:,jk,jsto) )
656            END DO
657            ! Factor to restore standard deviation after filtering
658            sto3d(:,:,jk,jsto) = sto3d(:,:,jk,jsto) * sto3d_fac(jsto)
659            ! Limit random parameter to the limitation factor
660            sto3d(:,:,jk,jsto) = SIGN(MIN(sto3d_lim(jsto),ABS(sto3d(:,:,jk,jsto))),sto3d(:,:,jk,jsto))
661            ! Multiply by standard devation and add average value
662            sto3d(:,:,jk,jsto) = sto3d(:,:,jk,jsto) * sto3d_std(jsto) + sto3d_ave(jsto)
663         END DO
664      END DO
665
666      ! 6) Restart stochastic parameters from file
667      ! ------------------------------------------
668      IF( ln_rststo ) CALL sto_rst_read
669
670      ! Allocate drho_ran
671      ALLOCATE(drho_ran(jpi,jpj,jpk))
672
673   END SUBROUTINE sto_par_init
674
675
676   SUBROUTINE sto_rst_read
677      !!----------------------------------------------------------------------
678      !!                  ***  ROUTINE sto_rst_read  ***
679      !!
680
681      !! ** Purpose :   read stochastic parameters from restart file
682      !!----------------------------------------------------------------------
683
684      INTEGER  :: jsto, jseed
685      INTEGER(KIND=8)     ::   ziseed(4)           ! RNG seeds in integer type
686      REAL(KIND=8)        ::   zrseed(4)           ! RNG seeds in real type (with same bits to save in restart)
687      CHARACTER(LEN=9)    ::   clsto2d='sto2d_000' ! stochastic parameter variable name
688      CHARACTER(LEN=9)    ::   clsto3d='sto3d_000' ! stochastic parameter variable name
689      CHARACTER(LEN=10)   ::   clseed='seed0_0000' ! seed variable name
690
691      IF ( jpsto2d > 0 .OR. jpsto3d > 0 ) THEN
692
693         IF(lwp) THEN
694            WRITE(numout,*)
695            WRITE(numout,*) 'sto_rst_read : read stochastic parameters from restart file'
696            WRITE(numout,*) '~~~~~~~~~~~~'
697         ENDIF
698
699         ! Open the restart file
700         CALL iom_open( cn_storst_in, numstor, kiolib = jprstlib )
701
702         ! Get stochastic parameters from restart file:
703         ! 2D stochastic parameters
704         DO jsto = 1 , jpsto2d
705            WRITE(clsto2d(7:9),'(i3.3)') jsto
706            CALL iom_get( numstor, jpdom_autoglo, clsto2d , sto2d(:,:,jsto) )
707         END DO
708         ! 3D stochastic parameters
709         DO jsto = 1 , jpsto3d
710            WRITE(clsto3d(7:9),'(i3.3)') jsto
711            CALL iom_get( numstor, jpdom_autoglo, clsto3d , sto3d(:,:,:,jsto) )
712         END DO
713
714         IF (ln_rstseed) THEN
715            ! Get saved state of the random number generator
716            DO jseed = 1 , 4
717               WRITE(clseed(5:5) ,'(i1.1)') jseed
718               WRITE(clseed(7:10),'(i4.4)') narea
719               CALL iom_get( numstor, clseed , zrseed(jseed) )
720            END DO
721            ziseed = TRANSFER( zrseed , ziseed)
722            CALL kiss_seed( ziseed(1) , ziseed(2) , ziseed(3) , ziseed(4) )
723         ENDIF
724
725         ! Close the restart file
726         CALL iom_close( numstor )
727
728      ENDIF
729
730   END SUBROUTINE sto_rst_read
731
732
733   SUBROUTINE sto_rst_write( kt )
734      !!----------------------------------------------------------------------
735      !!                  ***  ROUTINE sto_rst_write  ***
736      !!
737      !! ** Purpose :   write stochastic parameters in restart file
738      !!----------------------------------------------------------------------
739      INTEGER, INTENT(in) ::   kt     ! ocean time-step
740      !!
741      INTEGER  :: jsto, jseed
742      INTEGER(KIND=8)     ::   ziseed(4)           ! RNG seeds in integer type
743      REAL(KIND=8)        ::   zrseed(4)           ! RNG seeds in real type (with same bits to save in restart)
744      CHARACTER(LEN=20)   ::   clkt                ! ocean time-step defined as a character
745      CHARACTER(LEN=50)   ::   clname              ! restart file name
746      CHARACTER(LEN=9)    ::   clsto2d='sto2d_000' ! stochastic parameter variable name
747      CHARACTER(LEN=9)    ::   clsto3d='sto3d_000' ! stochastic parameter variable name
748      CHARACTER(LEN=10)   ::   clseed='seed0_0000' ! seed variable name
749
750      IF ( jpsto2d > 0 .OR. jpsto3d > 0 ) THEN
751
752         IF( kt == nitrst .OR. kt == nitend ) THEN
753            IF(lwp) THEN
754               WRITE(numout,*)
755               WRITE(numout,*) 'sto_rst_write : write stochastic parameters in restart file'
756               WRITE(numout,*) '~~~~~~~~~~~~~'
757            ENDIF
758         ENDIF
759
760         ! Put stochastic parameters in restart files
761         ! (as opened at previous timestep, see below)
762         IF( kt > nit000) THEN
763         IF( kt == nitrst .OR. kt == nitend ) THEN
764            ! get and save current state of the random number generator
765            CALL kiss_state( ziseed(1) , ziseed(2) , ziseed(3) , ziseed(4) )
766            zrseed = TRANSFER( ziseed , zrseed)
767            DO jseed = 1 , 4
768               WRITE(clseed(5:5) ,'(i1.1)') jseed
769               WRITE(clseed(7:10),'(i4.4)') narea
770               CALL iom_rstput( kt, nitrst, numstow, clseed , zrseed(jseed) )
771            END DO
772            ! 2D stochastic parameters
773            DO jsto = 1 , jpsto2d
774               WRITE(clsto2d(7:9),'(i3.3)') jsto
775               CALL iom_rstput( kt, nitrst, numstow, clsto2d , sto2d(:,:,jsto) )
776            END DO
777            ! 3D stochastic parameters
778            DO jsto = 1 , jpsto3d
779               WRITE(clsto3d(7:9),'(i3.3)') jsto
780               CALL iom_rstput( kt, nitrst, numstow, clsto3d , sto3d(:,:,:,jsto) )
781            END DO
782            ! Save drho_ran in restart file
783            CALL iom_rstput( kt, nitrst, numstow, 'drho' , drho_ran(:,:,:) )
784            ! close the restart file
785            CALL iom_close( numstow )
786         ENDIF
787         ENDIF
788
789         ! Open the restart file one timestep before writing restart
790         IF( kt < nitend) THEN
791         IF( kt == nitrst - 1 .OR. nstock == 1 .OR. kt == nitend-1 ) THEN
792            ! create the filename
793            IF( nitrst > 999999999 ) THEN   ;   WRITE(clkt, *       ) nitrst
794            ELSE                            ;   WRITE(clkt, '(i8.8)') nitrst
795            ENDIF
796            clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_storst_out)
797            ! print information
798            IF(lwp) THEN
799               WRITE(numout,*) '             open stochastic parameters restart file: '//clname
800               IF( kt == nitrst - 1 ) THEN
801                  WRITE(numout,*) '             kt = nitrst - 1 = ', kt
802               ELSE
803                  WRITE(numout,*) '             kt = '             , kt
804               ENDIF
805            ENDIF
806            ! open the restart file
807            CALL iom_open( clname, numstow, ldwrt = .TRUE., kiolib = jprstlib )
808         ENDIF
809         ENDIF
810
811      ENDIF
812
813   END SUBROUTINE sto_rst_write
814
815
816   SUBROUTINE sto_par_white( psto )
817      !!----------------------------------------------------------------------
818      !!                  ***  ROUTINE sto_par_white  ***
819      !!
820      !! ** Purpose :   fill input array with white Gaussian noise
821      !!----------------------------------------------------------------------
822      REAL(wp), DIMENSION(jpi,jpj), INTENT(out)           ::   psto
823      !!
824      INTEGER  :: ji, jj
825      REAL(KIND=8) :: gran   ! Gaussian random number (forced KIND=8 as in kiss_gaussian)
826
827      DO jj = 1, jpj
828         DO ji = 1, jpi
829            CALL kiss_gaussian( gran )
830            psto(ji,jj) = gran
831         END DO
832      END DO
833
834   END SUBROUTINE sto_par_white
835
836
837   SUBROUTINE sto_par_flt( psto )
838      !!----------------------------------------------------------------------
839      !!                  ***  ROUTINE sto_par_flt  ***
840      !!
841      !! ** Purpose :   apply horizontal Laplacian filter to input array
842      !!----------------------------------------------------------------------
843      REAL(wp), DIMENSION(jpi,jpj), INTENT(out)           ::   psto
844      !!
845      INTEGER  :: ji, jj
846
847      DO jj = 2, jpj-1
848         DO ji = 2, jpi-1
849            psto(ji,jj) = 0.5_wp * psto(ji,jj) + 0.125_wp * &
850                              &  ( psto(ji-1,jj) + psto(ji+1,jj) +  &
851                              &    psto(ji,jj-1) + psto(ji,jj+1) )
852         END DO
853      END DO
854
855   END SUBROUTINE sto_par_flt
856
857
858   REAL(wp) FUNCTION sto_par_flt_fac( kpasses )
859      !!----------------------------------------------------------------------
860      !!                  ***  FUNCTION sto_par_flt_fac  ***
861      !!
862      !! ** Purpose :   compute factor to restore standard deviation
863      !!                as a function of the number of passes
864      !!                of the Laplacian filter
865      !!----------------------------------------------------------------------
866      INTEGER, INTENT(in) :: kpasses
867      !!
868      INTEGER :: jpasses, ji, jj, jflti, jfltj
869      INTEGER, DIMENSION(-1:1,-1:1) :: pflt0
870      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pfltb
871      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pflta
872      REAL(wp) :: ratio
873
874      pflt0(-1,-1) = 0 ; pflt0(-1,0) = 1 ; pflt0(-1,1) = 0
875      pflt0( 0,-1) = 1 ; pflt0( 0,0) = 4 ; pflt0( 0,1) = 1
876      pflt0( 1,-1) = 0 ; pflt0( 1,0) = 1 ; pflt0( 1,1) = 0
877
878      ALLOCATE(pfltb(-kpasses-1:kpasses+1,-kpasses-1:kpasses+1))
879      ALLOCATE(pflta(-kpasses-1:kpasses+1,-kpasses-1:kpasses+1))
880
881      pfltb(:,:) = 0
882      pfltb(0,0) = 1
883      DO jpasses = 1, kpasses
884        pflta(:,:) = 0
885        DO jflti= -1, 1
886        DO jfltj= -1, 1
887          DO ji= -kpasses, kpasses
888          DO jj= -kpasses, kpasses
889            pflta(ji,jj) = pflta(ji,jj) + pfltb(ji+jflti,jj+jfltj) * pflt0(jflti,jfltj)
890          ENDDO
891          ENDDO
892        ENDDO
893        ENDDO
894        pfltb(:,:) = pflta(:,:)
895      ENDDO
896
897      ratio = SUM(pfltb(:,:))
898      ratio = ratio * ratio / SUM(pfltb(:,:)*pfltb(:,:))
899      ratio = SQRT(ratio)
900
901      DEALLOCATE(pfltb,pflta)
902
903      sto_par_flt_fac = ratio
904
905   END FUNCTION sto_par_flt_fac
906
907
908END MODULE stopar
909
Note: See TracBrowser for help on using the repository browser.