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/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/STO – NEMO

source: branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/STO/stopar.F90 @ 9366

Last change on this file since 9366 was 9366, checked in by andmirek, 6 years ago

#2050 first version. Compiled OK in moci test suite

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