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/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/STO – NEMO

source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/STO/stopar.F90 @ 11671

Last change on this file since 11671 was 11671, checked in by acc, 5 years ago

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