Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

stopar.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 5 years ago
The Dr Hook changes from my perl code.
File size: 44.9 KB

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

Note: See TracBrowser for help on using the repository browser.

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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/STO/stopar.F90 @ 11738

Download in other formats: