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 branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/STO – NEMO

source: branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/STO/stopar.F90 @ 10246

Last change on this file since 10246 was 10246, checked in by kingr, 5 years ago

Cleared SVN keywords.

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