source: NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/STO/stopar.F90 @ 13279

Last change on this file since 13279 was 13279, checked in by clem, 3 months ago

merge with r4.0-HEAD at r13278

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