New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
stopar.F90 in branches/UKMO/dev_r7750_GO6_package_oasis_timers/NEMOGCM/NEMO/OPA_SRC/STO – NEMO

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

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

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