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 NEMO/trunk/src/OCE/STO – NEMO

source: NEMO/trunk/src/OCE/STO/stopar.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

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