source: branches/UKMO/dev_r7750_GO6_package_oasis_timers/NEMOGCM/NEMO/OPA_SRC/STO/stopar.F90 @ 9261

Last change on this file since 9261 was 9261, checked in by andmirek, 3 years ago

#1978 only 1 and 2 for nn_timing

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