Context Navigation

stopar.F90 @ 9383

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

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/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/STO/stopar.F90 @ 9383

Download in other formats: