source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/OPA_SRC/STO/stopar.F90 @ 5443

Last change on this file since 5443 was 5443, checked in by davestorkey, 5 years ago

Update 2015/dev_r5021_UKMO1_CICE_coupling branch to revision 5442 of the trunk.

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$
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      IF ( lk_agrif ) CALL ctl_stop('EOS stochastic parametrization is not compatible with AGRIF')
266      ! Read namsto namelist : stochastic parameterization
267      REWIND( numnam_ref )              ! Namelist namdyn_adv in reference namelist : Momentum advection scheme
268      READ  ( numnam_ref, namsto, IOSTAT = ios, ERR = 901)
269901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsto in reference namelist', lwp )
270
271      REWIND( numnam_cfg )              ! Namelist namdyn_adv in configuration namelist : Momentum advection scheme
272      READ  ( numnam_cfg, namsto, IOSTAT = ios, ERR = 902 )
273902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsto in configuration namelist', lwp )
274      IF(lwm) WRITE ( numond, namsto )
275
276      !IF(ln_ens_rst_in) cn_storst_in = cn_mem//cn_storst_in
277
278      ! Parameter print
279      IF(lwp) THEN
280         WRITE(numout,*)
281         WRITE(numout,*) 'sto_par_init : stochastic parameterization'
282         WRITE(numout,*) '~~~~~~~~~~~~'
283         WRITE(numout,*) '   Namelist namsto : stochastic parameterization'
284         WRITE(numout,*) '      restart stochastic parameters           ln_rststo     = ', ln_rststo
285         WRITE(numout,*) '      read seed of RNG from restart file      ln_rstseed    = ', ln_rstseed
286         WRITE(numout,*) '      suffix of sto restart name (input)      cn_storst_in  = ', cn_storst_in
287         WRITE(numout,*) '      suffix of sto restart name (output)     cn_storst_out = ', cn_storst_out
288
289         ! WRITE(numout,*) '      stochastic lateral diffusion            ln_sto_ldf    = ', ln_sto_ldf
290         ! WRITE(numout,*) '      lateral diffusion std (in percent)      rn_ldf_std    = ', rn_ldf_std
291         ! WRITE(numout,*) '      lateral diffusion tcor (in timesteps)   rn_ldf_tcor   = ', rn_ldf_tcor
292
293         ! WRITE(numout,*) '      stochastic horizontal pressure gradient ln_sto_hpg    = ', ln_sto_hpg
294         ! WRITE(numout,*) '      density gradient std (in percent)       rn_hpg_std    = ', rn_hpg_std
295         ! WRITE(numout,*) '      density gradient tcor (in timesteps)    rn_hpg_tcor   = ', rn_hpg_tcor
296
297         ! WRITE(numout,*) '      stochastic ice strength                 ln_sto_pstar  = ', ln_sto_pstar
298         ! WRITE(numout,*) '      ice strength std (in percent)           rn_pstar_std  = ', rn_pstar_std
299         ! WRITE(numout,*) '      ice strength tcor (in timesteps)        rn_pstar_tcor = ', rn_pstar_tcor
300         ! WRITE(numout,*) '      order of autoregressive  processes      nn_pstar_ord  = ', nn_pstar_ord
301         ! WRITE(numout,*) '      passes of Laplacian filter              nn_pstar_flt  = ', nn_pstar_flt
302
303         !WRITE(numout,*) '      stochastic trend                        ln_sto_trd    = ', ln_sto_trd
304         !WRITE(numout,*) '      trend std (in percent)                  rn_trd_std    = ', rn_trd_std
305         !WRITE(numout,*) '      trend tcor (in timesteps)               rn_trd_tcor   = ', rn_trd_tcor
306
307         WRITE(numout,*) '      stochastic equation of state            ln_sto_eos    = ', ln_sto_eos
308         WRITE(numout,*) '      number of degrees of freedom            nn_sto_eos    = ', nn_sto_eos
309         WRITE(numout,*) '      random walk horz. std (in grid points)  rn_eos_stdxy  = ', rn_eos_stdxy
310         WRITE(numout,*) '      random walk vert. std (in grid points)  rn_eos_stdz   = ', rn_eos_stdz
311         WRITE(numout,*) '      random walk tcor (in timesteps)         rn_eos_tcor   = ', rn_eos_tcor
312         WRITE(numout,*) '      order of autoregressive  processes      nn_eos_ord    = ', nn_eos_ord
313         WRITE(numout,*) '      passes of Laplacian filter              nn_eos_flt    = ', nn_eos_flt
314         WRITE(numout,*) '      limitation factor                       rn_eos_lim    = ', rn_eos_lim
315
316         ! WRITE(numout,*) '      stochastic tracers dynamics             ln_sto_trc    = ', ln_sto_trc
317         ! WRITE(numout,*) '      number of degrees of freedom            nn_sto_trc    = ', nn_sto_trc
318         ! WRITE(numout,*) '      random walk horz. std (in grid points)  rn_trc_stdxy  = ', rn_trc_stdxy
319         ! WRITE(numout,*) '      random walk vert. std (in grid points)  rn_trc_stdz   = ', rn_trc_stdz
320         ! WRITE(numout,*) '      random walk tcor (in timesteps)         rn_trc_tcor   = ', rn_trc_tcor
321         ! WRITE(numout,*) '      order of autoregressive  processes      nn_trc_ord    = ', nn_trc_ord
322         ! WRITE(numout,*) '      passes of Laplacian filter              nn_trc_flt    = ', nn_trc_flt
323         ! WRITE(numout,*) '      limitation factor                       rn_trc_lim    = ', rn_trc_lim
324
325      ENDIF
326
327      IF(lwp) WRITE(numout,*)
328      IF(lwp) WRITE(numout,*) '   stochastic parameterization :'
329
330      ! Set number of 2D stochastic arrays
331      jpsto2d = 0
332      IF( ln_sto_ldf ) THEN
333         IF(lwp) WRITE(numout,*) '       - stochastic lateral diffusion'
334         jpsto2d   = jpsto2d + 1
335         jsto_ldf  = jpsto2d
336      ENDIF
337      IF( ln_sto_pstar ) THEN
338         IF(lwp) WRITE(numout,*) '       - stochastic ice strength'
339         jpsto2d    = jpsto2d + 1 * nn_pstar_ord
340         jsto_pstar = jpsto2d
341      ENDIF
342      IF( ln_sto_eos ) THEN
343         IF(lwp) WRITE(numout,*) '       - stochastic equation of state'
344         ALLOCATE(jsto_eosi(nn_sto_eos))
345         ALLOCATE(jsto_eosj(nn_sto_eos))
346         ALLOCATE(jsto_eosk(nn_sto_eos))
347         DO jdof = 1, nn_sto_eos
348            jpsto2d   = jpsto2d + 3 * nn_eos_ord
349            jsto_eosi(jdof) = jpsto2d - 2 * nn_eos_ord
350            jsto_eosj(jdof) = jpsto2d - 1 * nn_eos_ord
351            jsto_eosk(jdof) = jpsto2d
352         END DO
353      ELSE
354         nn_sto_eos = 0
355      ENDIF
356      IF( ln_sto_trc ) THEN
357         IF(lwp) WRITE(numout,*) '       - stochastic tracers dynamics'
358         ALLOCATE(jsto_trci(nn_sto_trc))
359         ALLOCATE(jsto_trcj(nn_sto_trc))
360         ALLOCATE(jsto_trck(nn_sto_trc))
361         DO jdof = 1, nn_sto_trc
362            jpsto2d   = jpsto2d + 3 * nn_trc_ord
363            jsto_trci(jdof) = jpsto2d - 2 * nn_trc_ord
364            jsto_trcj(jdof) = jpsto2d - 1 * nn_trc_ord
365            jsto_trck(jdof) = jpsto2d
366         END DO
367      ELSE
368         nn_sto_trc = 0
369      ENDIF
370
371      ! Set number of 3D stochastic arrays
372      jpsto3d = 0
373      IF( ln_sto_hpg ) THEN
374         IF(lwp) WRITE(numout,*) '       - stochastic horizontal pressure gradient'
375         jpsto3d   = jpsto3d + 2
376         jsto_hpgi = jpsto3d - 1
377         jsto_hpgj = jpsto3d
378      ENDIF
379      IF( ln_sto_trd ) THEN
380         IF(lwp) WRITE(numout,*) '       - stochastic trend'
381         jpsto3d   = jpsto3d + 1
382         jsto_trd  = jpsto3d
383      ENDIF
384
385      ! Allocate 2D stochastic arrays
386      IF ( jpsto2d > 0 ) THEN
387         ALLOCATE ( sto2d(jpi,jpj,jpsto2d) )
388         ALLOCATE ( sto2d_abc(jpsto2d,3) )
389         ALLOCATE ( sto2d_ave(jpsto2d) )
390         ALLOCATE ( sto2d_std(jpsto2d) )
391         ALLOCATE ( sto2d_lim(jpsto2d) )
392         ALLOCATE ( sto2d_tcor(jpsto2d) )
393         ALLOCATE ( sto2d_ord(jpsto2d) )
394         ALLOCATE ( sto2d_typ(jpsto2d) )
395         ALLOCATE ( sto2d_sgn(jpsto2d) )
396         ALLOCATE ( sto2d_flt(jpsto2d) )
397         ALLOCATE ( sto2d_fac(jpsto2d) )
398      ENDIF
399
400      ! Allocate 3D stochastic arrays
401      IF ( jpsto3d > 0 ) THEN
402         ALLOCATE ( sto3d(jpi,jpj,jpk,jpsto3d) )
403         ALLOCATE ( sto3d_abc(jpsto3d,3) )
404         ALLOCATE ( sto3d_ave(jpsto3d) )
405         ALLOCATE ( sto3d_std(jpsto3d) )
406         ALLOCATE ( sto3d_lim(jpsto3d) )
407         ALLOCATE ( sto3d_tcor(jpsto3d) )
408         ALLOCATE ( sto3d_ord(jpsto3d) )
409         ALLOCATE ( sto3d_typ(jpsto3d) )
410         ALLOCATE ( sto3d_sgn(jpsto3d) )
411         ALLOCATE ( sto3d_flt(jpsto3d) )
412         ALLOCATE ( sto3d_fac(jpsto3d) )
413      ENDIF
414
415      ! Allocate temporary workspace
416      IF ( jpsto2d > 0 .OR. jpsto3d > 0 ) THEN
417         ALLOCATE ( sto_tmp(jpi,jpj) ) ; sto_tmp(:,:) = 0._wp
418      ENDIF
419
420      ! 1) For every stochastic parameter:
421      ! ----------------------------------
422      ! - set nature of grid point and control of the sign
423      !       across the north fold (sto2d_typ, sto2d_sgn)
424      ! - set number of passes of Laplacian filter (sto2d_flt)
425      ! - set order of every autoregressive process (sto2d_ord)
426      DO jsto = 1, jpsto2d
427         sto2d_typ(jsto) = 'T'
428         sto2d_sgn(jsto) = 1._wp
429         sto2d_flt(jsto) = 0
430         sto2d_ord(jsto) = 1
431         DO jord = 0, nn_pstar_ord-1
432            IF ( jsto+jord == jsto_pstar ) THEN ! Stochastic ice strength (ave=1)
433               sto2d_ord(jsto) = nn_pstar_ord - jord
434               sto2d_flt(jsto) = nn_pstar_flt
435            ENDIF
436         ENDDO
437         DO jdof = 1, nn_sto_eos
438         DO jord = 0, nn_eos_ord-1
439            IF ( jsto+jord == jsto_eosi(jdof) ) THEN ! Stochastic equation of state i (ave=0)
440               sto2d_ord(jsto) = nn_eos_ord - jord
441               sto2d_sgn(jsto) = -1._wp
442               sto2d_flt(jsto) = nn_eos_flt
443            ENDIF
444            IF ( jsto+jord == jsto_eosj(jdof) ) THEN ! Stochastic equation of state j (ave=0)
445               sto2d_ord(jsto) = nn_eos_ord - jord
446               sto2d_sgn(jsto) = -1._wp
447               sto2d_flt(jsto) = nn_eos_flt
448            ENDIF
449            IF ( jsto+jord == jsto_eosk(jdof) ) THEN ! Stochastic equation of state k (ave=0)
450               sto2d_ord(jsto) = nn_eos_ord - jord
451               sto2d_flt(jsto) = nn_eos_flt
452            ENDIF
453         END DO
454         END DO
455         DO jdof = 1, nn_sto_trc
456         DO jord = 0, nn_trc_ord-1
457            IF ( jsto+jord == jsto_trci(jdof) ) THEN ! Stochastic tracers dynamics i (ave=0)
458               sto2d_ord(jsto) = nn_trc_ord - jord
459               sto2d_sgn(jsto) = -1._wp
460               sto2d_flt(jsto) = nn_trc_flt
461            ENDIF
462            IF ( jsto+jord == jsto_trcj(jdof) ) THEN ! Stochastic tracers dynamics j (ave=0)
463               sto2d_ord(jsto) = nn_trc_ord - jord
464               sto2d_sgn(jsto) = -1._wp
465               sto2d_flt(jsto) = nn_trc_flt
466            ENDIF
467            IF ( jsto+jord == jsto_trck(jdof) ) THEN ! Stochastic tracers dynamics k (ave=0)
468               sto2d_ord(jsto) = nn_trc_ord - jord
469               sto2d_flt(jsto) = nn_trc_flt
470            ENDIF
471         END DO
472         END DO
473
474         sto2d_fac(jsto) = sto_par_flt_fac ( sto2d_flt(jsto) )
475      END DO
476      !
477      DO jsto = 1, jpsto3d
478         sto3d_typ(jsto) = 'T'
479         sto3d_sgn(jsto) = 1._wp
480         sto3d_flt(jsto) = 0
481         sto3d_ord(jsto) = 1
482         IF ( jsto == jsto_hpgi ) THEN ! Stochastic density gradient i (ave=1)
483            sto3d_typ(jsto) = 'U'
484         ENDIF
485         IF ( jsto == jsto_hpgj ) THEN ! Stochastic density gradient j (ave=1)
486            sto3d_typ(jsto) = 'V'
487         ENDIF
488         sto3d_fac(jsto) = sto_par_flt_fac ( sto3d_flt(jsto) )
489      END DO
490
491      ! 2) For every stochastic parameter:
492      ! ----------------------------------
493      ! set average, standard deviation and time correlation
494      DO jsto = 1, jpsto2d
495         sto2d_ave(jsto)  = 0._wp
496         sto2d_std(jsto)  = 1._wp
497         sto2d_tcor(jsto) = 1._wp
498         sto2d_lim(jsto)  = 3._wp
499         IF ( jsto == jsto_ldf  ) THEN ! Stochastic lateral diffusion (ave=1)
500            sto2d_ave(jsto)  = 1._wp
501            sto2d_std(jsto)  = rn_ldf_std
502            sto2d_tcor(jsto) = rn_ldf_tcor
503         ENDIF
504         DO jord = 0, nn_pstar_ord-1
505            IF ( jsto+jord == jsto_pstar ) THEN ! Stochastic ice strength (ave=1)
506               sto2d_std(jsto) = 1._wp
507               sto2d_tcor(jsto) = rn_pstar_tcor
508            ENDIF
509         ENDDO
510         DO jdof = 1, nn_sto_eos
511         DO jord = 0, nn_eos_ord-1
512            IF ( jsto+jord == jsto_eosi(jdof) ) THEN ! Stochastic equation of state i (ave=0)
513               sto2d_std(jsto)  = rn_eos_stdxy
514               sto2d_tcor(jsto) = rn_eos_tcor
515               sto2d_lim(jsto)  = rn_eos_lim
516            ENDIF
517            IF ( jsto+jord == jsto_eosj(jdof) ) THEN ! Stochastic equation of state j (ave=0)
518               sto2d_std(jsto)  = rn_eos_stdxy
519               sto2d_tcor(jsto) = rn_eos_tcor
520               sto2d_lim(jsto)  = rn_eos_lim
521            ENDIF
522            IF ( jsto+jord == jsto_eosk(jdof) ) THEN ! Stochastic equation of state k (ave=0)
523               sto2d_std(jsto)  = rn_eos_stdz
524               sto2d_tcor(jsto) = rn_eos_tcor
525               sto2d_lim(jsto)  = rn_eos_lim
526            ENDIF
527         END DO
528         END DO
529         DO jdof = 1, nn_sto_trc
530         DO jord = 0, nn_trc_ord-1
531            IF ( jsto+jord == jsto_trci(jdof) ) THEN ! Stochastic tracer dynamics i (ave=0)
532               sto2d_std(jsto)  = rn_trc_stdxy
533               sto2d_tcor(jsto) = rn_trc_tcor
534               sto2d_lim(jsto)  = rn_trc_lim
535            ENDIF
536            IF ( jsto+jord == jsto_trcj(jdof) ) THEN ! Stochastic tracer dynamics j (ave=0)
537               sto2d_std(jsto)  = rn_trc_stdxy
538               sto2d_tcor(jsto) = rn_trc_tcor
539               sto2d_lim(jsto)  = rn_trc_lim
540            ENDIF
541            IF ( jsto+jord == jsto_trck(jdof) ) THEN ! Stochastic tracer dynamics k (ave=0)
542               sto2d_std(jsto)  = rn_trc_stdz
543               sto2d_tcor(jsto) = rn_trc_tcor
544               sto2d_lim(jsto)  = rn_trc_lim
545            ENDIF
546         END DO
547         END DO
548
549      END DO
550      !
551      DO jsto = 1, jpsto3d
552         sto3d_ave(jsto)  = 0._wp
553         sto3d_std(jsto)  = 1._wp
554         sto3d_tcor(jsto) = 1._wp
555         sto3d_lim(jsto)  = 3._wp
556         IF ( jsto == jsto_hpgi ) THEN ! Stochastic density gradient i (ave=1)
557            sto3d_ave(jsto)  = 1._wp
558            sto3d_std(jsto)  = rn_hpg_std
559            sto3d_tcor(jsto) = rn_hpg_tcor
560         ENDIF
561         IF ( jsto == jsto_hpgj ) THEN ! Stochastic density gradient j (ave=1)
562            sto3d_ave(jsto)  = 1._wp
563            sto3d_std(jsto)  = rn_hpg_std
564            sto3d_tcor(jsto) = rn_hpg_tcor
565         ENDIF
566         IF ( jsto == jsto_trd ) THEN ! Stochastic trend (ave=1)
567            sto3d_ave(jsto)  = 1._wp
568            sto3d_std(jsto)  = rn_trd_std
569            sto3d_tcor(jsto) = rn_trd_tcor
570         ENDIF
571      END DO
572
573      ! 3) For every stochastic parameter:
574      ! ----------------------------------
575      ! - compute parameters (a, b, c) of the AR1 autoregressive process
576      !   from expected value (ave), standard deviation (std)
577      !   and time correlation (tcor):
578      !     a = EXP ( - 1 / tcor )           --> sto2d_abc(:,1)
579      !     b = std * SQRT( 1 - a * a )      --> sto2d_abc(:,2)
580      !     c = ave * ( 1 - a )              --> sto2d_abc(:,3)
581      ! - for higher order processes (ARn, n>1), use approximate formula
582      !   for the b parameter (valid for tcor>>1 time step)
583      DO jsto = 1, jpsto2d
584         IF ( sto2d_tcor(jsto) == 0._wp ) THEN
585            sto2d_abc(jsto,1) = 0._wp
586         ELSE
587            sto2d_abc(jsto,1) = EXP ( - 1._wp / sto2d_tcor(jsto) )
588         ENDIF
589         IF ( sto2d_ord(jsto) == 1 ) THEN      ! Exact formula for 1st order process
590            rinflate = sto2d_std(jsto)
591         ELSE
592            ! Approximate formula, valid for tcor >> 1
593            jordm1 = sto2d_ord(jsto) - 1
594            rinflate = SQRT ( REAL( jordm1 , wp ) / REAL( 2*(2*jordm1-1) , wp ) )
595         ENDIF
596         sto2d_abc(jsto,2) = rinflate * SQRT ( 1._wp - sto2d_abc(jsto,1) &
597                                                     * sto2d_abc(jsto,1) )
598         sto2d_abc(jsto,3) = sto2d_ave(jsto) * ( 1._wp - sto2d_abc(jsto,1) )
599      END DO
600      !
601      DO jsto = 1, jpsto3d
602         IF ( sto3d_tcor(jsto) == 0._wp ) THEN
603            sto3d_abc(jsto,1) = 0._wp
604         ELSE
605            sto3d_abc(jsto,1) = EXP ( - 1._wp / sto3d_tcor(jsto) )
606         ENDIF
607         IF ( sto3d_ord(jsto) == 1 ) THEN      ! Exact formula for 1st order process
608            rinflate = sto3d_std(jsto)
609         ELSE
610            ! Approximate formula, valid for tcor >> 1
611            jordm1 = sto3d_ord(jsto) - 1
612            rinflate = SQRT ( REAL( jordm1 , wp ) / REAL( 2*(2*jordm1-1) , wp ) )
613         ENDIF
614         sto3d_abc(jsto,2) = rinflate * SQRT ( 1._wp - sto3d_abc(jsto,1) &
615                                                     * sto3d_abc(jsto,1) )
616         sto3d_abc(jsto,3) = sto3d_ave(jsto) * ( 1._wp - sto3d_abc(jsto,1) )
617      END DO
618
619      ! 4) Initialize seeds for random number generator
620      ! -----------------------------------------------
621      ! using different seeds for different processors (jarea)
622      ! and different ensemble members (jmem)
623      CALL kiss_reset( )
624      DO jarea = 1, narea
625         !DO jmem = 0, nmember
626         zseed1 = kiss() ; zseed2 = kiss() ; zseed3 = kiss() ; zseed4 = kiss()
627         !END DO
628      END DO
629      CALL kiss_seed( zseed1, zseed2, zseed3, zseed4 )
630
631      ! 5) Initialize stochastic parameters to: ave + std * w
632      ! -----------------------------------------------------
633      DO jsto = 1, jpsto2d
634         ! Draw random numbers from N(0,1) --> w
635         CALL sto_par_white( sto2d(:,:,jsto) )
636         ! Apply horizontal Laplacian filter to w
637         DO jflt = 1, sto2d_flt(jsto)
638            CALL lbc_lnk( sto2d(:,:,jsto), sto2d_typ(jsto), sto2d_sgn(jsto) )
639            CALL sto_par_flt( sto2d(:,:,jsto) )
640         END DO
641         ! Factor to restore standard deviation after filtering
642         sto2d(:,:,jsto) = sto2d(:,:,jsto) * sto2d_fac(jsto)
643         ! Limit random parameter to the limitation factor
644         sto2d(:,:,jsto) = SIGN(MIN(sto2d_lim(jsto),ABS(sto2d(:,:,jsto))),sto2d(:,:,jsto))
645         ! Multiply by standard devation and add average value
646         sto2d(:,:,jsto) = sto2d(:,:,jsto) * sto2d_std(jsto) + sto2d_ave(jsto)
647      END DO
648      !
649      DO jsto = 1, jpsto3d
650         DO jk = 1, jpk
651            ! Draw random numbers from N(0,1) --> w
652            CALL sto_par_white( sto3d(:,:,jk,jsto) )
653            ! Apply horizontal Laplacian filter to w
654            DO jflt = 1, sto3d_flt(jsto)
655               CALL lbc_lnk( sto3d(:,:,jk,jsto), sto3d_typ(jsto), sto3d_sgn(jsto) )
656               CALL sto_par_flt( sto3d(:,:,jk,jsto) )
657            END DO
658            ! Factor to restore standard deviation after filtering
659            sto3d(:,:,jk,jsto) = sto3d(:,:,jk,jsto) * sto3d_fac(jsto)
660            ! Limit random parameter to the limitation factor
661            sto3d(:,:,jk,jsto) = SIGN(MIN(sto3d_lim(jsto),ABS(sto3d(:,:,jk,jsto))),sto3d(:,:,jk,jsto))
662            ! Multiply by standard devation and add average value
663            sto3d(:,:,jk,jsto) = sto3d(:,:,jk,jsto) * sto3d_std(jsto) + sto3d_ave(jsto)
664         END DO
665      END DO
666
667      ! 6) Restart stochastic parameters from file
668      ! ------------------------------------------
669      IF( ln_rststo ) CALL sto_rst_read
670
671      ! Allocate drho_ran
672      ALLOCATE(drho_ran(jpi,jpj,jpk))
673
674   END SUBROUTINE sto_par_init
675
676
677   SUBROUTINE sto_rst_read
678      !!----------------------------------------------------------------------
679      !!                  ***  ROUTINE sto_rst_read  ***
680      !!
681
682      !! ** Purpose :   read stochastic parameters from restart file
683      !!----------------------------------------------------------------------
684
685      INTEGER  :: jsto, jseed
686      INTEGER(KIND=8)     ::   ziseed(4)           ! RNG seeds in integer type
687      REAL(KIND=8)        ::   zrseed(4)           ! RNG seeds in real type (with same bits to save in restart)
688      CHARACTER(LEN=9)    ::   clsto2d='sto2d_000' ! stochastic parameter variable name
689      CHARACTER(LEN=9)    ::   clsto3d='sto3d_000' ! stochastic parameter variable name
690      CHARACTER(LEN=10)   ::   clseed='seed0_0000' ! seed variable name
691
692      IF ( jpsto2d > 0 .OR. jpsto3d > 0 ) THEN
693
694         IF(lwp) THEN
695            WRITE(numout,*)
696            WRITE(numout,*) 'sto_rst_read : read stochastic parameters from restart file'
697            WRITE(numout,*) '~~~~~~~~~~~~'
698         ENDIF
699
700         ! Open the restart file
701         CALL iom_open( cn_storst_in, numstor, kiolib = jprstlib )
702
703         ! Get stochastic parameters from restart file:
704         ! 2D stochastic parameters
705         DO jsto = 1 , jpsto2d
706            WRITE(clsto2d(7:9),'(i3.3)') jsto
707            CALL iom_get( numstor, jpdom_autoglo, clsto2d , sto2d(:,:,jsto) )
708         END DO
709         ! 3D stochastic parameters
710         DO jsto = 1 , jpsto3d
711            WRITE(clsto3d(7:9),'(i3.3)') jsto
712            CALL iom_get( numstor, jpdom_autoglo, clsto3d , sto3d(:,:,:,jsto) )
713         END DO
714
715         IF (ln_rstseed) THEN
716            ! Get saved state of the random number generator
717            DO jseed = 1 , 4
718               WRITE(clseed(5:5) ,'(i1.1)') jseed
719               WRITE(clseed(7:10),'(i4.4)') narea
720               CALL iom_get( numstor, clseed , zrseed(jseed) )
721            END DO
722            ziseed = TRANSFER( zrseed , ziseed)
723            CALL kiss_seed( ziseed(1) , ziseed(2) , ziseed(3) , ziseed(4) )
724         ENDIF
725
726         ! Close the restart file
727         CALL iom_close( numstor )
728
729      ENDIF
730
731   END SUBROUTINE sto_rst_read
732
733
734   SUBROUTINE sto_rst_write( kt )
735      !!----------------------------------------------------------------------
736      !!                  ***  ROUTINE sto_rst_write  ***
737      !!
738      !! ** Purpose :   write stochastic parameters in restart file
739      !!----------------------------------------------------------------------
740      INTEGER, INTENT(in) ::   kt     ! ocean time-step
741      !!
742      INTEGER  :: jsto, jseed
743      INTEGER(KIND=8)     ::   ziseed(4)           ! RNG seeds in integer type
744      REAL(KIND=8)        ::   zrseed(4)           ! RNG seeds in real type (with same bits to save in restart)
745      CHARACTER(LEN=20)   ::   clkt                ! ocean time-step defined as a character
746      CHARACTER(LEN=50)   ::   clname              ! restart file name
747      CHARACTER(LEN=9)    ::   clsto2d='sto2d_000' ! stochastic parameter variable name
748      CHARACTER(LEN=9)    ::   clsto3d='sto3d_000' ! stochastic parameter variable name
749      CHARACTER(LEN=10)   ::   clseed='seed0_0000' ! seed variable name
750
751      IF ( jpsto2d > 0 .OR. jpsto3d > 0 ) THEN
752
753         IF( kt == nitrst .OR. kt == nitend ) THEN
754            IF(lwp) THEN
755               WRITE(numout,*)
756               WRITE(numout,*) 'sto_rst_write : write stochastic parameters in restart file'
757               WRITE(numout,*) '~~~~~~~~~~~~~'
758            ENDIF
759         ENDIF
760
761         ! Put stochastic parameters in restart files
762         ! (as opened at previous timestep, see below)
763         IF( kt > nit000) THEN
764         IF( kt == nitrst .OR. kt == nitend ) THEN
765            ! get and save current state of the random number generator
766            CALL kiss_state( ziseed(1) , ziseed(2) , ziseed(3) , ziseed(4) )
767            zrseed = TRANSFER( ziseed , zrseed)
768            DO jseed = 1 , 4
769               WRITE(clseed(5:5) ,'(i1.1)') jseed
770               WRITE(clseed(7:10),'(i4.4)') narea
771               CALL iom_rstput( kt, nitrst, numstow, clseed , zrseed(jseed) )
772            END DO
773            ! 2D stochastic parameters
774            DO jsto = 1 , jpsto2d
775               WRITE(clsto2d(7:9),'(i3.3)') jsto
776               CALL iom_rstput( kt, nitrst, numstow, clsto2d , sto2d(:,:,jsto) )
777            END DO
778            ! 3D stochastic parameters
779            DO jsto = 1 , jpsto3d
780               WRITE(clsto3d(7:9),'(i3.3)') jsto
781               CALL iom_rstput( kt, nitrst, numstow, clsto3d , sto3d(:,:,:,jsto) )
782            END DO
783            ! Save drho_ran in restart file
784            CALL iom_rstput( kt, nitrst, numstow, 'drho' , drho_ran(:,:,:) )
785            ! close the restart file
786            CALL iom_close( numstow )
787         ENDIF
788         ENDIF
789
790         ! Open the restart file one timestep before writing restart
791         IF( kt < nitend) THEN
792         IF( kt == nitrst - 1 .OR. nstock == 1 .OR. kt == nitend-1 ) THEN
793            ! create the filename
794            IF( nitrst > 999999999 ) THEN   ;   WRITE(clkt, *       ) nitrst
795            ELSE                            ;   WRITE(clkt, '(i8.8)') nitrst
796            ENDIF
797            clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_"//TRIM(cn_storst_out)
798            ! print information
799            IF(lwp) THEN
800               WRITE(numout,*) '             open stochastic parameters restart file: '//clname
801               IF( kt == nitrst - 1 ) THEN
802                  WRITE(numout,*) '             kt = nitrst - 1 = ', kt
803               ELSE
804                  WRITE(numout,*) '             kt = '             , kt
805               ENDIF
806            ENDIF
807            ! open the restart file
808            CALL iom_open( clname, numstow, ldwrt = .TRUE., kiolib = jprstlib )
809         ENDIF
810         ENDIF
811
812      ENDIF
813
814   END SUBROUTINE sto_rst_write
815
816
817   SUBROUTINE sto_par_white( psto )
818      !!----------------------------------------------------------------------
819      !!                  ***  ROUTINE sto_par_white  ***
820      !!
821      !! ** Purpose :   fill input array with white Gaussian noise
822      !!----------------------------------------------------------------------
823      REAL(wp), DIMENSION(jpi,jpj), INTENT(out)           ::   psto
824      !!
825      INTEGER  :: ji, jj
826      REAL(KIND=8) :: gran   ! Gaussian random number (forced KIND=8 as in kiss_gaussian)
827
828      DO jj = 1, jpj
829         DO ji = 1, jpi
830            CALL kiss_gaussian( gran )
831            psto(ji,jj) = gran
832         END DO
833      END DO
834
835   END SUBROUTINE sto_par_white
836
837
838   SUBROUTINE sto_par_flt( psto )
839      !!----------------------------------------------------------------------
840      !!                  ***  ROUTINE sto_par_flt  ***
841      !!
842      !! ** Purpose :   apply horizontal Laplacian filter to input array
843      !!----------------------------------------------------------------------
844      REAL(wp), DIMENSION(jpi,jpj), INTENT(out)           ::   psto
845      !!
846      INTEGER  :: ji, jj
847
848      DO jj = 2, jpj-1
849         DO ji = 2, jpi-1
850            psto(ji,jj) = 0.5_wp * psto(ji,jj) + 0.125_wp * &
851                              &  ( psto(ji-1,jj) + psto(ji+1,jj) +  &
852                              &    psto(ji,jj-1) + psto(ji,jj+1) )
853         END DO
854      END DO
855
856   END SUBROUTINE sto_par_flt
857
858
859   REAL(wp) FUNCTION sto_par_flt_fac( kpasses )
860      !!----------------------------------------------------------------------
861      !!                  ***  FUNCTION sto_par_flt_fac  ***
862      !!
863      !! ** Purpose :   compute factor to restore standard deviation
864      !!                as a function of the number of passes
865      !!                of the Laplacian filter
866      !!----------------------------------------------------------------------
867      INTEGER, INTENT(in) :: kpasses
868      !!
869      INTEGER :: jpasses, ji, jj, jflti, jfltj
870      INTEGER, DIMENSION(-1:1,-1:1) :: pflt0
871      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pfltb
872      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pflta
873      REAL(wp) :: ratio
874
875      pflt0(-1,-1) = 0 ; pflt0(-1,0) = 1 ; pflt0(-1,1) = 0
876      pflt0( 0,-1) = 1 ; pflt0( 0,0) = 4 ; pflt0( 0,1) = 1
877      pflt0( 1,-1) = 0 ; pflt0( 1,0) = 1 ; pflt0( 1,1) = 0
878
879      ALLOCATE(pfltb(-kpasses-1:kpasses+1,-kpasses-1:kpasses+1))
880      ALLOCATE(pflta(-kpasses-1:kpasses+1,-kpasses-1:kpasses+1))
881
882      pfltb(:,:) = 0
883      pfltb(0,0) = 1
884      DO jpasses = 1, kpasses
885        pflta(:,:) = 0
886        DO jflti= -1, 1
887        DO jfltj= -1, 1
888          DO ji= -kpasses, kpasses
889          DO jj= -kpasses, kpasses
890            pflta(ji,jj) = pflta(ji,jj) + pfltb(ji+jflti,jj+jfltj) * pflt0(jflti,jfltj)
891          ENDDO
892          ENDDO
893        ENDDO
894        ENDDO
895        pfltb(:,:) = pflta(:,:)
896      ENDDO
897
898      ratio = SUM(pfltb(:,:))
899      ratio = ratio * ratio / SUM(pfltb(:,:)*pfltb(:,:))
900      ratio = SQRT(ratio)
901
902      DEALLOCATE(pfltb,pflta)
903
904      sto_par_flt_fac = ratio
905
906   END FUNCTION sto_par_flt_fac
907
908
909END MODULE stopar
910
Note: See TracBrowser for help on using the repository browser.