Context Navigation

stopack.F90 @ 15690

Last change on this file since 15690 was 15690, checked in by mattmartin, 2 years ago
Fixed minor typos from Jonah's review in ticket https://code.metoffice.gov.uk/trac/nwpscience/ticket/1125
File size: 114.0 KB

Line
1	!#define NEMO_V34
2
3	! Portability to V34 is achieved through (aggressive) precompiler directives:
4	! Uncomment the first line for NEMO v3.4
5
6	#if defined NEMO_V34
7
8	! Change name of trend indexes
9	#define jptra_xad jptra_trd_xad
10	#define jptra_yad jptra_trd_yad
11	#define jptra_zad jptra_trd_zad
12	#define jptra_ldf jptra_trd_ldf
13	#define jptra_zdf jptra_trd_zdf
14	#define jptra_bbc jptra_trd_bbc
15	#define jptra_bbl jptra_trd_bbl
16	#define jptra_npc jptra_trd_npc
17	#define jptra_dmp jptra_trd_dmp
18	#define jptra_qsr jptra_trd_qsr
19	#define jptra_nsr jptra_trd_nsr
20	#define jptra_atf jptra_trd_atf
21
22	#define jptra_sad -999
23	#define jptra_zdfp -999
24	#define jptra_evd -999
25
26	#define jpdyn_hpg jpdyn_trd_hpg
27	#define jpdyn_keg jpdyn_trd_keg
28	#define jpdyn_rvo jpdyn_trd_rvo
29	#define jpdyn_pvo jpdyn_trd_pvo
30	#define jpdyn_ldf jpdyn_trd_ldf
31	#define jpdyn_had jpdyn_trd_had
32	#define jpdyn_zad jpdyn_trd_zad
33	#define jpdyn_zdf jpdyn_trd_zdf
34	#define jpdyn_spg jpdyn_trd_spg
35	#define jpdyn_dat jpdyn_trd_dat
36	#define jpdyn_swf jpdyn_trd_swf
37	#define jpdyn_bfr jpdyn_trd_bfr
38
39	#define jpdyn_spgflt -999
40	#define jpdyn_spgexp -999
41	#define jpdyn_atf -999
42	#define jpdyn_tau -999
43	#define jpdyn_bfri -999
44
45	! Change name of namelist units
46	#define numnam_ref numnam
47	#define numnam_cfg numnam
48	#define lwm lwp
49	#define numond numout
50
51	#define wmask tmask
52
53	#endif
54
55	MODULE stopack
56	!!======================================================================
57	!! * MODULE stopack *
58	!! Calculate and Apply sotchastic physics perturbations
59	!!======================================================================
60	!! History : 1.0 ! 2018-02 (A. Storto), Original SPPT code
61	!! for NEMO 3.6
62	!! 2.0 ! 2019-05 (A. Storto), upgrades and updates:
63	!! (SPP, SKEB and sea-ice)
64	!!----------------------------------------------------------------------
65	!!
66	!! stopack : Generate stochastic physics perturbations
67	!!
68	!! Method
69	!! ======
70	!! The module allows users to activate:
71	!! - SPPT (Stochastically perturbed parameterization
72	!! tendencies )scheme for user-selected trends for
73	!! tracers, momentum and sea-ice
74	!! - SPP (Schastically perturbed parameters) scheme
75	!! for some (namelist) parameters
76	!! - SEB (Stochastic Energy backscatter)
77	!! backscatter energy dissipated numerically or
78	!! through deep convection.
79	!!
80	!!
81	!! Acknowledgements: C3S funded ERGO project
82	!!
83	!!----------------------------------------------------------------------
84	USE par_kind
85	USE timing ! Timing
86	USE oce ! ocean dynamics and tracers variables
87	USE dom_oce ! ocean space and time domain variables
88	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
89	USE in_out_manager ! I/O manager
90	#ifdef NEMO_V34
91	USE trdmod_oce
92	USE restart
93	#else
94	USE trd_oce
95	#endif
96	USE iom ! I/O manager library
97	USE lib_mpp ! MPP library
98	#if defined key_lim3
99	USE ice ! LIM-3 variables
100	#endif
101	#if defined key_lim2
102	USE ice_2 ! LIM-2 variables
103	#endif
104	USE wrk_nemo
105	USE diaptr
106	USE zdf_oce
107	USE phycst
108
109	IMPLICIT NONE
110	PRIVATE
111
112	! Accessible routines
113	PUBLIC tra_sppt_collect
114	PUBLIC dyn_sppt_collect
115	PUBLIC tra_sppt_apply
116	PUBLIC dyn_sppt_apply
117	PUBLIC stopack_rst
118	PUBLIC stopack_init
119	PUBLIC stopack_pert
120	PUBLIC sppt_ice
121	PUBLIC spp_aht
122	PUBLIC spp_ahm
123	PUBLIC spp_gen
124	PUBLIC skeb_comp
125	PUBLIC skeb_apply
126
127	! Main logical switch
128	LOGICAL, PUBLIC :: ln_stopack
129
130	! Internal switches for SPPT
131	LOGICAL, PUBLIC, SAVE :: ln_sppt_tra = .FALSE.
132	LOGICAL, PUBLIC, SAVE :: ln_sppt_dyn = .FALSE.
133	LOGICAL, PUBLIC, SAVE :: ln_sppt_ice = .FALSE.
134
135	! SPPT Options
136	INTEGER :: sppt_filter_pass, sppt_step, nn_vwei, &
137	& nn_sppt_step_bound, nn_rndm_freq, nn_deftau
138	REAL(wp), SAVE :: rn_sppt_tau, rn_sppt_bound, rn_distcoast, rn_sppt_stdev
139	REAL(wp), SAVE :: rn_skeb_tau, rn_skeb_stdev
140	REAL(wp), SAVE :: rn_spp_tau, rn_spp_stdev
141	INTEGER :: skeb_filter_pass, spp_filter_pass
142
143	! SPPT Logical switches for individual tendencies
144	LOGICAL :: ln_sppt_taumap, ln_stopack_restart, ln_distcoast, &
145	ln_sppt_traxad, ln_sppt_trayad, ln_sppt_trazad, ln_sppt_trasad, ln_sppt_traldf, &
146	ln_sppt_trazdf, ln_sppt_trazdfp,ln_sppt_traevd, ln_sppt_trabbc, ln_sppt_trabbl, &
147	ln_sppt_tranpc, ln_sppt_tradmp, ln_sppt_traqsr, ln_sppt_transr, ln_sppt_traatf
148	LOGICAL :: &
149	ln_sppt_dynhpg, ln_sppt_dynspg, ln_sppt_dynkeg, ln_sppt_dynrvo, ln_sppt_dynpvo, ln_sppt_dynzad,&
150	ln_sppt_dynldf, ln_sppt_dynzdf, ln_sppt_dynbfr, ln_sppt_dynatf, ln_sppt_dyntau, ln_sppt_glocon
151	LOGICAL, PUBLIC :: &
152	ln_sppt_icehdf, ln_sppt_icelat, ln_sppt_icezdf
153
154	! Arrays
155	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,: ) :: gauss_n_2d
156	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,: ) :: gauss_n_2d_p
157	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,: ) :: gauss_n_2d_k
158	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,: ) :: coeff_bfr
159
160	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gauss_n, zdc
161	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gauss_nu, gauss_nv
162	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: spptt, sppts, spptu, spptv
163	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,: ) :: gauss_b, sppt_tau, sppt_a, sppt_b
164	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,: ) :: skeb_a, skeb_b
165	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,: ) :: spp_a, spp_b
166	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,: ) :: gauss_bp, gauss_bk
167	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,: ) :: g2d_save
168	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(: ) :: gauss_w
169	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zicewrk
170	REAL(wp), SAVE :: flt_fac,flt_fac_k,flt_fac_p
171	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,: ) :: spp_tau
172	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,: ) :: skeb_tau
173
174	INTEGER, PARAMETER, PUBLIC :: jk_spp_alb = 1
175	INTEGER, PARAMETER, PUBLIC :: jk_spp_rhg = 2
176	INTEGER, PARAMETER, PUBLIC :: jk_spp_relw = 3
177	INTEGER, PARAMETER, PUBLIC :: jk_spp_dqdt = 4
178	INTEGER, PARAMETER, PUBLIC :: jk_spp_deds = 5
179	INTEGER, PARAMETER, PUBLIC :: jk_spp_arnf = 6
180	INTEGER, PARAMETER, PUBLIC :: jk_spp_geot = 7
181	INTEGER, PARAMETER, PUBLIC :: jk_spp_qsi0 = 8
182	INTEGER, PARAMETER, PUBLIC :: jk_spp_bfr = 9
183	INTEGER, PARAMETER, PUBLIC :: jk_spp_aevd = 10
184	INTEGER, PARAMETER, PUBLIC :: jk_spp_avt = 11
185	INTEGER, PARAMETER, PUBLIC :: jk_spp_avm = 12
186	INTEGER, PARAMETER, PUBLIC :: jk_spp_tkelc = 13
187	INTEGER, PARAMETER, PUBLIC :: jk_spp_tkedf = 14
188	INTEGER, PARAMETER, PUBLIC :: jk_spp_tkeds = 15
189	INTEGER, PARAMETER, PUBLIC :: jk_spp_tkebb = 16
190	INTEGER, PARAMETER, PUBLIC :: jk_spp_tkefr = 17
191
192	INTEGER, PARAMETER, PUBLIC :: jk_spp_ahtu = 18
193	INTEGER, PARAMETER, PUBLIC :: jk_spp_ahtv = 19
194	INTEGER, PARAMETER, PUBLIC :: jk_spp_ahtw = 20
195	INTEGER, PARAMETER, PUBLIC :: jk_spp_ahtt = 21
196
197	INTEGER, PARAMETER, PUBLIC :: jk_spp_ahubbl = 22
198	INTEGER, PARAMETER, PUBLIC :: jk_spp_ahvbbl = 23
199
200	INTEGER, PARAMETER, PUBLIC :: jk_spp_ahm1 = 24
201	INTEGER, PARAMETER, PUBLIC :: jk_spp_ahm2 = 25
202	INTEGER, PARAMETER, PUBLIC :: jk_spp_ahm3 = 26
203	INTEGER, PARAMETER, PUBLIC :: jk_spp_ahm4 = 27
204
205	INTEGER, PARAMETER, PUBLIC :: jk_skeb_dnum = 31
206	INTEGER, PARAMETER, PUBLIC :: jk_skeb_dcon = 32
207	INTEGER, PARAMETER, PUBLIC :: jk_skeb_tot = 33
208
209	INTEGER, PARAMETER, PUBLIC :: jk_sppt_tem = 41
210	INTEGER, PARAMETER, PUBLIC :: jk_sppt_sal = 42
211	INTEGER, PARAMETER, PUBLIC :: jk_sppt_uvl = 43
212	INTEGER, PARAMETER, PUBLIC :: jk_sppt_vvl = 44
213
214	INTEGER, PARAMETER, PUBLIC :: jk_spp = 27
215	INTEGER, PARAMETER, PUBLIC :: jk_stpk_tot = 44
216	REAL(wp), SAVE :: rn_mmax( jk_stpk_tot ) = 0._wp
217
218	INTEGER, SAVE :: numdiag = 601
219	INTEGER, SAVE :: numrep = 602
220	INTEGER, SAVE :: lkt
221
222	! Randome generator seed
223	INTEGER, SAVE :: nn_stopack_seed(4)
224	LOGICAL, SAVE, PUBLIC :: ln_stopack_repr = .TRUE.
225
226	! SPP switches
227	INTEGER, PUBLIC, SAVE :: nn_spp_bfr,nn_spp_dqdt,nn_spp_dedt,nn_spp_avt,nn_spp_avm,nn_spp_qsi0,nn_spp_relw, nn_spp_arnf,nn_spp_geot,nn_spp_aevd,nn_spp_ahubbl,nn_spp_ahvbbl
228	INTEGER, PUBLIC, SAVE :: nn_spp_tkelc,nn_spp_tkedf,nn_spp_tkeds,nn_spp_tkebb,nn_spp_tkefr
229	INTEGER, PUBLIC, SAVE :: nn_spp_ahtu, nn_spp_ahtv, nn_spp_ahtw, nn_spp_ahtt
230	INTEGER, PUBLIC, SAVE :: nn_spp_ahm1, nn_spp_ahm2
231	! SPP Sea-ice switches
232	INTEGER, PUBLIC, SAVE :: nn_spp_icealb,nn_spp_icestr
233
234	! SPP parameters
235	REAL(wp), PUBLIC, SAVE :: rn_bfr_sd, rn_dqdt_sd, rn_dedt_sd, rn_avt_sd, rn_avm_sd, rn_qsi0_sd, rn_relw_sd, rn_arnf_sd, rn_geot_sd, rn_aevd_sd, rn_ahubbl_sd, rn_ahvbbl_sd
236	REAL(wp), PUBLIC, SAVE :: rn_tkelc_sd,rn_tkedf_sd,rn_tkeds_sd,rn_tkebb_sd,rn_tkefr_sd
237	REAL(wp), PUBLIC, SAVE :: rn_ahtu_sd, rn_ahtv_sd, rn_ahtw_sd, rn_ahtt_sd
238	REAL(wp), PUBLIC, SAVE :: rn_ahm1_sd, rn_ahm2_sd
239	REAL(wp), PUBLIC, SAVE :: rn_icestr_sd, rn_icealb_sd
240
241	! Internal switches for SPP
242	LOGICAL, SAVE :: ln_spp_own_gauss
243	INTEGER :: nn_spp
244
245	LOGICAL, SAVE :: ln_stopack_diags
246
247	LOGICAL, SAVE :: ln_skeb_own_gauss
248	LOGICAL, SAVE, PUBLIC :: ln_skeb
249	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,: ) :: dpsiv, dpsiu
250	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dnum , dcon
251	INTEGER, SAVE :: jpri,jprj
252	INTEGER, SAVE :: nn_dcom_freq = 1
253	REAL(wp), SAVE :: rn_skeb, rn_kh, rn_kc
254	LOGICAL, SAVE, PUBLIC :: ln_dpsiv = .FALSE.
255	LOGICAL, SAVE, PUBLIC :: ln_dpsiu = .FALSE.
256	REAL(wp), SAVE :: rn_beta_num, rn_beta_con
257	INTEGER, SAVE :: nn_dconv = 1
258
259	! Default/initial seeds
260	INTEGER(KIND=i8) :: x=1234567890987654321_8
261	INTEGER(KIND=i8) :: y=362436362436362436_8
262	INTEGER(KIND=i8) :: z=1066149217761810_8
263	INTEGER(KIND=i8) :: w=123456123456123456_8
264	! Parameters to generate real random variates
265	REAL(KIND=wp), PARAMETER :: huge64=9223372036854775808.0 ! +1
266	REAL(KIND=wp), PARAMETER :: zero=0.0, half=0.5, one=1.0, two=2.0
267	! Variables to store 2 Gaussian random numbers with current index (ig)
268	INTEGER(KIND=i8), SAVE :: ig=1
269	REAL(KIND=wp), SAVE :: gran1, gran2
270
271	!! * Substitutions
272	# include "domzgr_substitute.h90"
273	# include "vectopt_loop_substitute.h90"
274
275	!!----------------------------------------------------------------------
276	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
277	!! $Id: bdytra.F90 4292 2013-11-20 16:28:04Z cetlod $
278	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
279	!!----------------------------------------------------------------------
280
281	CONTAINS
282
283	#if defined NEMO_V34
284	SUBROUTINE tra_sppt_collect( ptrdx, ptrdy, ktrd )
285	#else
286	SUBROUTINE tra_sppt_collect( ptrdx, ptrdy, ktrd, kt )
287	!!----------------------------------------------------------------------
288	!! * ROUTINE tra_sppt_collect *
289	!!
290	!! ** Purpose : Collect tracer tendencies (additive)
291	!! This function is called by the tendency diagnostics
292	!! module
293	!!----------------------------------------------------------------------
294	INTEGER , INTENT(in ) :: kt ! time step
295	#endif
296	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ptrdx ! Temperature
297	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ptrdy ! Salinity
298	INTEGER , INTENT(in ) :: ktrd ! tracer trend index
299
300	IF( ktrd .eq. jptra_xad .AND. ln_sppt_traxad ) THEN
301	spptt = spptt + ptrdx ; sppts = sppts + ptrdy
302	ENDIF
303	IF( ktrd .eq. jptra_yad .AND. ln_sppt_trayad ) THEN
304	spptt = spptt + ptrdx ; sppts = sppts + ptrdy
305	ENDIF
306	IF( ktrd .eq. jptra_zad .AND. ln_sppt_trazad ) THEN
307	spptt = spptt + ptrdx ; sppts = sppts + ptrdy
308	ENDIF
309	IF( ktrd .eq. jptra_sad .AND. ln_sppt_trasad ) THEN
310	spptt = spptt + ptrdx ; sppts = sppts + ptrdy
311	ENDIF
312	IF( ktrd .eq. jptra_ldf .AND. ln_sppt_traldf ) THEN
313	spptt = spptt + ptrdx ; sppts = sppts + ptrdy
314	ENDIF
315	IF( ktrd .eq. jptra_zdf .AND. ln_sppt_trazdf ) THEN
316	spptt = spptt + ptrdx ; sppts = sppts + ptrdy
317	ENDIF
318	IF( ktrd .eq. jptra_zdfp.AND. ln_sppt_trazdfp) THEN
319	spptt = spptt + ptrdx ; sppts = sppts + ptrdy
320	ENDIF
321	IF( ktrd .eq. jptra_evd .AND. ln_sppt_traevd ) THEN
322	spptt = spptt + ptrdx ; sppts = sppts + ptrdy
323	ENDIF
324	IF( ktrd .eq. jptra_bbc .AND. ln_sppt_trabbc ) THEN
325	spptt = spptt + ptrdx ; sppts = sppts + ptrdy
326	ENDIF
327	IF( ktrd .eq. jptra_bbl .AND. ln_sppt_trabbl ) THEN
328	spptt = spptt + ptrdx ; sppts = sppts + ptrdy
329	ENDIF
330	IF( ktrd .eq. jptra_npc .AND. ln_sppt_tranpc ) THEN
331	spptt = spptt + ptrdx ; sppts = sppts + ptrdy
332	ENDIF
333	IF( ktrd .eq. jptra_dmp .AND. ln_sppt_tradmp ) THEN
334	spptt = spptt + ptrdx ; sppts = sppts + ptrdy
335	ENDIF
336	IF( ktrd .eq. jptra_qsr .AND. ln_sppt_traqsr ) THEN
337	spptt = spptt + ptrdx ; sppts = sppts + ptrdy
338	ENDIF
339	IF( ktrd .eq. jptra_nsr .AND. ln_sppt_transr ) THEN
340	spptt = spptt + ptrdx ; sppts = sppts + ptrdy
341	ENDIF
342	IF( ktrd .eq. jptra_atf .AND. ln_sppt_traatf ) THEN
343	spptt = spptt + ptrdx ; sppts = sppts + ptrdy
344	ENDIF
345
346	END SUBROUTINE tra_sppt_collect
347
348	#if defined NEMO_V34
349	SUBROUTINE dyn_sppt_collect( ptrdx, ptrdy, ktrd )
350	#else
351	SUBROUTINE dyn_sppt_collect( ptrdx, ptrdy, ktrd, kt )
352	!!----------------------------------------------------------------------
353	!! * ROUTINE dyn_sppt_collect *
354	!!
355	!! ** Purpose : Collect momentum tendencies (additive)
356	!! This function is called by the tendency diagnostics
357	!! module
358	!!----------------------------------------------------------------------
359	INTEGER , INTENT(in ) :: kt ! time step
360	#endif
361	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ptrdx ! Temperature
362	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ptrdy ! Salinity
363	INTEGER , INTENT(in ) :: ktrd ! tracer trend index
364
365	IF( ktrd .eq. jpdyn_hpg .AND. ln_sppt_dynhpg ) THEN
366	spptu = spptu + ptrdx ; spptv = spptv + ptrdy
367	ENDIF
368	IF( ktrd .eq. jpdyn_spg .AND. ln_sppt_dynspg ) THEN
369	spptu = spptu + ptrdx ; spptv = spptv + ptrdy
370	ENDIF
371	IF( ktrd .eq. jpdyn_spgflt .AND. ln_sppt_dynspg ) THEN
372	spptu = spptu + ptrdx ; spptv = spptv + ptrdy
373	ENDIF
374	IF( ktrd .eq. jpdyn_spgexp .AND. ln_sppt_dynspg ) THEN
375	spptu = spptu + ptrdx ; spptv = spptv + ptrdy
376	ENDIF
377	IF( ktrd .eq. jpdyn_keg .AND. ln_sppt_dynkeg ) THEN
378	spptu = spptu + ptrdx ; spptv = spptv + ptrdy
379	ENDIF
380	IF( ktrd .eq. jpdyn_rvo .AND. ln_sppt_dynrvo ) THEN
381	spptu = spptu + ptrdx ; spptv = spptv + ptrdy
382	ENDIF
383	IF( ktrd .eq. jpdyn_pvo .AND. ln_sppt_dynpvo ) THEN
384	spptu = spptu + ptrdx ; spptv = spptv + ptrdy
385	ENDIF
386	IF( ktrd .eq. jpdyn_zad .AND. ln_sppt_dynzad ) THEN
387	spptu = spptu + ptrdx ; spptv = spptv + ptrdy
388	ENDIF
389	IF( ktrd .eq. jpdyn_ldf .AND. ln_sppt_dynldf ) THEN
390	spptu = spptu + ptrdx ; spptv = spptv + ptrdy
391	ENDIF
392	IF( ktrd .eq. jpdyn_zdf .AND. ln_sppt_dynzdf ) THEN
393	spptu = spptu + ptrdx ; spptv = spptv + ptrdy
394	ENDIF
395	IF( ktrd .eq. jpdyn_bfr .AND. ln_sppt_dynbfr ) THEN
396	spptu = spptu + ptrdx ; spptv = spptv + ptrdy
397	ENDIF
398	IF( ktrd .eq. jpdyn_atf .AND. ln_sppt_dynatf ) THEN
399	spptu = spptu + ptrdx ; spptv = spptv + ptrdy
400	ENDIF
401	IF( ktrd .eq. jpdyn_tau .AND. ln_sppt_dyntau ) THEN
402	spptu = spptu + ptrdx ; spptv = spptv + ptrdy
403	ENDIF
404	IF( ktrd .eq. jpdyn_bfri.AND. ln_sppt_dynbfr ) THEN
405	spptu = spptu + ptrdx ; spptv = spptv + ptrdy
406	ENDIF
407
408	END SUBROUTINE dyn_sppt_collect
409
410	SUBROUTINE stopack_pert( kt )
411	!!----------------------------------------------------------------------
412	!! * ROUTINE stopack_pert *
413	!!
414	!! ** Purpose : Update perturbation every nn_rndm_freq timestep
415	!! Calculate for SPPT, eventually for SPP and SKEB
416	!! If they are activeated and with different error scales
417	!!----------------------------------------------------------------------
418	INTEGER, INTENT( in ) :: kt
419	INTEGER :: ji,jj
420
421	IF( ln_stopack_diags ) lkt = kt
422	IF( MOD( kt - 1, nn_rndm_freq ) == 0 .OR. kt == nit000 ) THEN
423
424	IF(lwp) THEN
425	WRITE(numout,*)
426	WRITE(numout,*) ' Calculating perturbation at timestep ', kt
427	WRITE(numout,*)
428	ENDIF
429
430	! Generate red noise
431	CALL gaussian_ar1_field ( gauss_n , sppt_filter_pass, gauss_w , gauss_b,&
432	sppt_a, sppt_b, gauss_n_2d,flt_fac,1)
433
434	! Generate red noise for SPP, SKEB if required
435	IF( ln_spp_own_gauss ) CALL gaussian_ar1_field ( gauss_n , spp_filter_pass, gauss_w , gauss_bp,&
436	spp_a, spp_b, gauss_n_2d_p,flt_fac_p,2)
437	IF( ln_skeb_own_gauss ) CALL gaussian_ar1_field ( gauss_n ,skeb_filter_pass, gauss_w , gauss_bk,&
438	skeb_a, skeb_b, gauss_n_2d_k,flt_fac_k,3)
439
440	! Staggering on U-,V- grids
441	! (later should account also for possibly different bottom levels)
442	DO jj=1,jpj-1
443	DO ji=1,jpi-1
444	gauss_nu(ji,jj,:) = 0.5_wp * (gauss_n(ji,jj,:)+gauss_n(ji+1,jj,:))
445	gauss_nv(ji,jj,:) = 0.5_wp * (gauss_n(ji,jj,:)+gauss_n(ji,jj+1,:))
446	ENDDO
447	ENDDO
448	CALL lbc_lnk( gauss_nu, 'U', -1._wp)
449	CALL lbc_lnk( gauss_nv, 'V', -1._wp)
450
451	ENDIF
452
453	#ifdef key_iomput
454	! Output the perturbation field
455	CALL iom_put( "sppt_ar1" , gauss_n )
456	IF(ln_spp_own_gauss) CALL iom_put( "spp_ar1" , gauss_n_2d_p)
457	IF(ln_skeb_own_gauss) CALL iom_put( "skeb_ar1" , gauss_n_2d_k)
458	#endif
459
460	END SUBROUTINE stopack_pert
461
462	#if defined key_lim2
463	SUBROUTINE sppt_ice( kstep, kopt, z1,z2,z3,z4,z5,z6,z7 )
464	#else
465	SUBROUTINE sppt_ice( kstep, kopt)
466	#endif
467	!!----------------------------------------------------------------------
468	!! * ROUTINE sppt_ice *
469	!!
470	!! ** Purpose : Apply collinear perturbation to ice fields
471	!! For specific processes coded in LIM2/LIM3
472	!!----------------------------------------------------------------------
473	!
474	INTEGER, INTENT(IN) :: kstep ! Start (1) or End (2) of ice routine
475	INTEGER, INTENT(IN) :: kopt ! Option for sea-ice routine
476	#if defined key_lim2
477	REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(INOUT) :: z1,z2,z3,z4,z5,z6,z7
478	#endif
479	#if defined key_lim3 \|\| defined key_lim2
480	INTEGER :: jmt ! Number of sea-ice variables (depends on LIM version, process)
481	INTEGER :: jm, jl, jk
482
483	#if defined key_lim3
484	jmt = 3jpl+jplnlay_i
485	IF( kopt == 2 ) jmt=jmt+3*jpl+1
486	#else
487	IF( kopt == 1 ) jmt=6
488	IF( kopt == 2 ) jmt=7
489	IF( kopt == 3 ) jmt=8
490	#endif
491	IF( kopt == 4 ) jmt=2
492
493	IF( kstep == 1 ) THEN ! Store the before values
494	IF ( ALLOCATED ( zicewrk) ) DEALLOCATE ( zicewrk )
495	ALLOCATE ( zicewrk(jpi,jpj,jmt) )
496	jm=1
497	#if defined key_lim3
498	DO jl = 1, jpl
499	zicewrk(:,:,jm) = a_i(:,:,jl) ; jm=jm+1
500	zicewrk(:,:,jm) = v_i(:,:,jl) ; jm=jm+1
501	zicewrk(:,:,jm) =smv_i(:,:,jl); jm=jm+1
502	DO jk = 1, nlay_i
503	zicewrk(:,:,jm) = e_i(:,:,jk,jl) ; jm=jm+1
504	END DO
505	END DO
506	IF( kopt .EQ. 2 ) THEN
507	DO jl = 1, jpl
508	zicewrk(:,:,jm) =oa_i(:,:,jl) ; jm=jm+1
509	zicewrk(:,:,jm) = e_s(:,:,1,jl) ; jm=jm+1
510	zicewrk(:,:,jm) = v_s(:,:,jl) ; jm=jm+1
511	ENDDO
512	zicewrk(:,:,jm) =ato_i(:,:); jm=jm+1
513	ENDIF
514	#else
515	IF( kopt .EQ. 1 ) THEN
516	zicewrk(:,:,jm) = frld ; jm=jm+1
517	zicewrk(:,:,jm) = hsnif ; jm=jm+1
518	zicewrk(:,:,jm) = hicif ; jm=jm+1
519	zicewrk(:,:,jm) = tbif(:,:,1) ; jm=jm+1
520	zicewrk(:,:,jm) = tbif(:,:,2) ; jm=jm+1
521	zicewrk(:,:,jm) = tbif(:,:,3) ; jm=jm+1
522	ENDIF
523	IF( kopt .EQ. 2 ) THEN
524	IF(.NOT. PRESENT(z1) ) CALL ctl_stop( 'sppt icehdf problem, step 1')
525	zicewrk(:,:,jm) = z1 ; jm=jm+1
526	zicewrk(:,:,jm) = z2 ; jm=jm+1
527	zicewrk(:,:,jm) = z3 ; jm=jm+1
528	zicewrk(:,:,jm) = z4 ; jm=jm+1
529	zicewrk(:,:,jm) = z5 ; jm=jm+1
530	zicewrk(:,:,jm) = z6 ; jm=jm+1
531	zicewrk(:,:,jm) = z7
532	ENDIF
533	IF( kopt .EQ. 3 ) THEN
534	zicewrk(:,:,jm) = frld ; jm=jm+1
535	zicewrk(:,:,jm) = hsnif ; jm=jm+1
536	zicewrk(:,:,jm) = hicif ; jm=jm+1
537	zicewrk(:,:,jm) = sist ; jm=jm+1
538	zicewrk(:,:,jm) = tbif(:,:,1) ; jm=jm+1
539	zicewrk(:,:,jm) = tbif(:,:,2) ; jm=jm+1
540	zicewrk(:,:,jm) = tbif(:,:,3) ; jm=jm+1
541	ENDIF
542	#endif
543	IF ( kopt .EQ. 4 ) THEN
544	zicewrk(:,:,jm) = u_ice ; jm=jm+1
545	zicewrk(:,:,jm) = v_ice ; jm=jm+1
546	ENDIF
547	ELSEIF ( kstep == 2 ) THEN ! Add collinear perturbation
548	jm=1
549	#if defined key_lim3
550	DO jl = 1, jpl
551	a_i(:,:,jl) = a_i(:,:,jl) + (a_i(:,:,jl)-zicewrk(:,:,jm))*gauss_n_2d ; jm=jm+1 ; CALL lbc_lnk(a_i(:,:,jl), 'T', 1. )
552	v_i(:,:,jl) = v_i(:,:,jl) + (v_i(:,:,jl)-zicewrk(:,:,jm))*gauss_n_2d ; jm=jm+1 ; CALL lbc_lnk(v_i(:,:,jl), 'T', 1. )
553	smv_i(:,:,jl)=smv_i(:,:,jl)+ (smv_i(:,:,jl)-zicewrk(:,:,jm))*gauss_n_2d ; jm=jm+1 ; CALL lbc_lnk(smv_i(:,:,jl), 'T', 1. )
554	DO jk = 1, nlay_i
555	e_i(:,:,jk,jl) = e_i(:,:,jk,jl)+(e_i(:,:,jk,jl)-zicewrk(:,:,jm))*gauss_n_2d ; jm=jm+1 ; CALL lbc_lnk(e_i(:,:,jk,jl), 'T', 1. )
556	END DO
557	END DO
558	IF( kopt .EQ. 2 ) THEN
559	DO jl = 1, jpl
560	oa_i(:,:,jl)=oa_i(:,:,jl)+(oa_i(:,:,jl)-zicewrk(:,:,jm))*gauss_n_2d ; jm=jm+1 ; CALL lbc_lnk(oa_i(:,:,jl), 'T', 1. )
561	e_s(:,:,1,jl)= e_s(:,:,1,jl)+( e_s(:,:,1,jl)-zicewrk(:,:,jm))*gauss_n_2d ; jm=jm+1 ; CALL lbc_lnk(e_s(:,:,1,jl), 'T', 1. )
562	v_s(:,:,jl)= v_s(:,:,jl)+( v_s(:,:,jl)-zicewrk(:,:,jm))*gauss_n_2d ; jm=jm+1 ; CALL lbc_lnk(v_s(:,:,jl), 'T', 1. )
563	ENDDO
564	ato_i(:,:)=ato_i(:,:)+(ato_i(:,:)-zicewrk(:,:,jm))*gauss_n_2d ; jm=jm+1 ; CALL lbc_lnk(ato_i(:,:), 'T', 1. )
565	ENDIF
566	DEALLOCATE ( zicewrk )
567	#else
568	IF( kopt .EQ. 1 ) THEN
569	frld = frld + gauss_n_2d * ( frld - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk(frld, 'T', 1. )
570	hsnif = hsnif + gauss_n_2d * ( hsnif - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk(hsnif, 'T', 1. )
571	hicif = hicif + gauss_n_2d * ( hicif - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk(hicif, 'T', 1. )
572	tbif(:,:,1)=tbif(:,:,1) + gauss_n_2d * ( tbif(:,:,1) - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk(tbif(:,:,1), 'T', 1. )
573	tbif(:,:,2)=tbif(:,:,2) + gauss_n_2d * ( tbif(:,:,2) - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk(tbif(:,:,2), 'T', 1. )
574	tbif(:,:,3)=tbif(:,:,3) + gauss_n_2d * ( tbif(:,:,3) - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk(tbif(:,:,3), 'T', 1. )
575	ENDIF
576	IF( kopt .EQ. 2 ) THEN
577	IF(.NOT. PRESENT(z1) ) CALL ctl_stop( 'sppt icehdf problem, step 2')
578	z1 = z1 + gauss_n_2d * ( z1 - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk( z1, 'T', 1. )
579	z2 = z2 + gauss_n_2d * ( z2 - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk( z2, 'T', 1. )
580	z3 = z3 + gauss_n_2d * ( z3 - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk( z3, 'T', 1. )
581	z4 = z4 + gauss_n_2d * ( z4 - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk( z4, 'T', 1. )
582	z5 = z5 + gauss_n_2d * ( z5 - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk( z5, 'T', 1. )
583	z6 = z6 + gauss_n_2d * ( z6 - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk( z6, 'T', 1. )
584	z7 = z7 + gauss_n_2d * ( z7 - zicewrk(:,:,jm) ) ; CALL lbc_lnk( z7, 'T', 1. )
585	ENDIF
586	IF( kopt .EQ. 3 ) THEN
587	frld = frld + gauss_n_2d * ( frld - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk(frld, 'T', 1. )
588	hsnif = hsnif + gauss_n_2d * ( hsnif - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk(hsnif, 'T', 1. )
589	hicif = hicif + gauss_n_2d * ( hicif - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk(hicif, 'T', 1. )
590	sist = sist + gauss_n_2d * ( sist - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk(sist, 'T', 1. )
591	tbif(:,:,1)=tbif(:,:,1) + gauss_n_2d * ( tbif(:,:,1) - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk(tbif(:,:,1), 'T', 1. )
592	tbif(:,:,2)=tbif(:,:,2) + gauss_n_2d * ( tbif(:,:,2) - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk(tbif(:,:,2), 'T', 1. )
593	tbif(:,:,3)=tbif(:,:,3) + gauss_n_2d * ( tbif(:,:,3) - zicewrk(:,:,jm) ) ; jm=jm+1 ; CALL lbc_lnk(tbif(:,:,3), 'T', 1. )
594	ENDIF
595	#endif
596	! EVP or VP rheology common to LIM2-EVP and LIM3
597	IF ( kopt .EQ. 4 ) THEN
598	u_ice = u_ice + (u_ice - zicewrk(:,:,jm) ) * gauss_n_2d ; jm=jm+1
599	v_ice = v_ice + (v_ice - zicewrk(:,:,jm) ) * gauss_n_2d ; jm=jm+1
600	#if defined key_lim3 \|\| ( defined key_lim2 && ! defined key_lim2_vp )
601	CALL lbc_lnk( u_ice, 'U', -1. )
602	CALL lbc_lnk( v_ice, 'V', -1. )
603	#endif
604	#if defined key_lim2 && defined key_lim2_vp
605	CALL lbc_lnk( u_ice(:,1:jpj), 'I', -1. )
606	CALL lbc_lnk( v_ice(:,1:jpj), 'I', -1. )
607	#endif
608	ENDIF
609	DEALLOCATE ( zicewrk )
610	ENDIF
611	#endif
612	END SUBROUTINE sppt_ice
613
614	SUBROUTINE spp_gen ( kt, coeff, nn_type, rn_sd, kspp, klev )
615	!!----------------------------------------------------------------------
616	!! * ROUTINE spp_gen *
617	!!
618	!! ** Purpose : Perturbing parameters (generic function)
619	!! Given a value of standard deviation, the 2D parameter
620	!! coeff is perturbed following an additive Normal,
621	!! a lognormal with mean the original coeff,
622	!! a lognormal with median the original coeff,
623	!! or the previous functions but with value bounded [0.1]
624	!!----------------------------------------------------------------------
625	INTEGER, INTENT( in ) :: kt
626	REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj) :: coeff
627	INTEGER, INTENT( in ) :: nn_type
628	REAL(wp), INTENT( in ) :: rn_sd
629	INTEGER, INTENT( in ) :: kspp
630	INTEGER, INTENT( in ), OPTIONAL :: klev
631	REAL(wp), POINTER, DIMENSION(:,:) :: gauss
632	REAL(wp) :: zsd,xme,mm
633	CHARACTER (LEN=99) :: cstrng
634	INTEGER :: jklev
635
636	#if defined key_traldf_c3d \|\| key_traldf_c2d
637
638	CALL wrk_alloc(jpi,jpj,gauss)
639
640	IF( ln_spp_own_gauss ) THEN
641	gauss = gauss_n_2d_p
642	ELSE
643	gauss = rn_spp_stdev * gauss_n_2d / rn_sppt_stdev
644	ENDIF
645
646	IF( nn_type == 1 ) THEN
647	gauss = gauss * rn_sd
648	coeff = coeff * ( 1._wp + gauss )
649	ELSEIF ( nn_type == 2 ) THEN
650	zsd = rn_sd
651	xme = -0.5_wpzsdzsd
652	gauss = gauss * zsd + xme
653	coeff = exp(gauss) * coeff
654	ELSEIF ( nn_type == 3 ) THEN
655	zsd = rn_sd
656	xme = 0._wp
657	gauss = gauss * zsd + xme
658	coeff = exp(gauss) * coeff
659	ELSEIF( nn_type == 4 ) THEN
660	gauss = gauss * rn_sd
661	coeff = coeff * ( 1._wp + gauss )
662	WHERE( coeff > 1._wp ) coeff = 1._wp
663	WHERE( coeff < 0._wp ) coeff = 0._wp
664	ELSEIF ( nn_type == 5 ) THEN
665	zsd = rn_sd
666	xme = -0.5_wpzsdzsd
667	gauss = gauss * zsd + xme
668	coeff = exp(gauss) * coeff
669	WHERE( coeff > 1._wp ) coeff = 1._wp
670	WHERE( coeff < 0._wp ) coeff = 0._wp
671	ELSEIF ( nn_type == 6 ) THEN
672	zsd = rn_sd
673	xme = 0._wp
674	gauss = gauss * zsd + xme
675	coeff = exp(gauss) * coeff
676	WHERE( coeff > 1._wp ) coeff = 1._wp
677	WHERE( coeff < 0._wp ) coeff = 0._wp
678	ELSE
679	CALL ctl_stop( 'spp dqdt wrong option')
680	ENDIF
681
682	#ifdef key_iomput
683	write(cstrng,'(A,I2.2)') 'spp_par',kspp
684	IF (iom_use(TRIM(cstrng)) ) CALL iom_put( TRIM(cstrng) , coeff )
685	#endif
686
687	IF( ln_stopack_diags ) THEN
688	IF(PRESENT(klev)) THEN
689	jklev = klev
690	ELSE
691	jklev = 0
692	ENDIF
693	CALL spp_stats(kt,kspp,jklev,coeff)
694	ENDIF
695
696	CALL wrk_dealloc(jpi,jpj,gauss)
697
698	#else
699	CALL ctl_stop('key_traldf_c1d is not a valid key for STO')
700	#endif
701
702	END SUBROUTINE
703
704	SUBROUTINE spp_stats(mt,kp,kl,rcf)
705	!!----------------------------------------------------------------------
706	!! * ROUTINE spp_stats *
707	!!
708	!! ** Purpose : Compute and print basic SPP statistics
709	!!----------------------------------------------------------------------
710	IMPLICIT NONE
711	INTEGER, INTENT(IN) :: mt,kp,kl
712	REAL(wp), INTENT(IN) :: rcf(jpi,jpj)
713	REAL(wp) :: mi,ma
714	CHARACTER(LEN=16) :: cstr = ' '
715	SELECT CASE ( kp )
716	CASE( jk_spp_alb )
717	cstr = 'ALBEDO '
718	CASE( jk_spp_rhg )
719	cstr = 'ICE RHEOLOGY'
720	CASE( jk_spp_relw )
721	cstr = 'RELATIVE WND'
722	CASE( jk_spp_dqdt )
723	cstr = 'SST RELAXAT.'
724	CASE( jk_spp_deds )
725	cstr = 'SSS RELAXAT.'
726	CASE( jk_spp_arnf )
727	cstr = 'RIVER MIXING'
728	CASE( jk_spp_geot )
729	cstr = 'GEOTHERM.FLX'
730	CASE( jk_spp_qsi0 )
731	cstr = 'SOLAR EXTIN.'
732	CASE( jk_spp_bfr )
733	cstr = 'BOTTOM FRICT'
734	CASE( jk_spp_aevd )
735	cstr = 'EDDY VISCDIF'
736	CASE( jk_spp_avt )
737	cstr = 'VERT. DIFFUS'
738	CASE( jk_spp_avm )
739	cstr = 'VERT. VISCOS'
740	CASE( jk_spp_tkelc )
741	cstr = 'TKE LANGMUIR'
742	CASE( jk_spp_tkedf )
743	cstr = 'TKE RN_EDIFF'
744	CASE( jk_spp_tkeds )
745	cstr = 'TKE RN_EDISS'
746	CASE( jk_spp_tkebb )
747	cstr = 'TKE RN_EBB '
748	CASE( jk_spp_tkefr )
749	cstr = 'TKE RN_EFR '
750	CASE( jk_spp_ahtu )
751	cstr = 'TRALDF AHTU '
752	CASE( jk_spp_ahtv )
753	cstr = 'TRALDF AHTV '
754	CASE( jk_spp_ahtw )
755	cstr = 'TRALDF AHTW '
756	CASE( jk_spp_ahtt )
757	cstr = 'TRALDF AHTT '
758	CASE( jk_spp_ahubbl )
759	cstr = 'BBL DIFFUSU '
760	CASE( jk_spp_ahvbbl )
761	cstr = 'BBL DIFFUSV '
762	CASE( jk_spp_ahm1 )
763	cstr = 'DYNLDF AHM1 '
764	CASE( jk_spp_ahm2 )
765	cstr = 'DYNLDF AHM2 '
766	CASE( jk_spp_ahm3 )
767	cstr = 'DYNLDF AHM3 '
768	CASE( jk_spp_ahm4 )
769	cstr = 'DYNLDF AHM4 '
770	CASE DEFAULT
771	CALL ctl_stop('Unrecognized SPP parameter: add it or turn off diagnostics')
772	END SELECT
773	mi = MINVAL(rcf)
774	ma = MAXVAL(rcf)
775	IF(lk_mpp) CALL mpp_min(mi)
776	IF(lk_mpp) CALL mpp_max(ma)
777	IF(lwp) THEN
778	IF ( kl > 0 ) write(cstr(14:16),'(I3.3)') kl
779	WRITE(numdiag,9300) lkt,cstr,mi,ma
780	ENDIF
781	9300 FORMAT(' it :', i8, ' ', A16, ' Min: ',d10.3 , ' Max: ',d10.3)
782	rn_mmax ( kp ) = MAX ( MAX( abs(mi), ma), rn_mmax ( kp ) )
783	END SUBROUTINE spp_stats
784
785	SUBROUTINE stopack_report
786	!!----------------------------------------------------------------------
787	!! * ROUTINE stopack_report *
788	!!
789	!! ** Purpose : Report at the end of the run
790	!!----------------------------------------------------------------------
791	IMPLICIT NONE
792	INTEGER :: jp, numrep
793	CHARACTER(LEN=16) :: cstr = ' '
794	REAL(wp) :: zmul
795
796	CALL ctl_opn(numrep, 'stopack.report', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
797
798	WRITE(numrep,'(A)') ' REPORT from STOPACK stochastic physics package'
799	WRITE(numrep,'(A)')
800	WRITE(numrep,'(A)') ' SPP part'
801
802	DO jp=1,jk_spp
803	SELECT CASE ( jp )
804	CASE( jk_spp_alb )
805	cstr = 'ALBEDO '
806	CASE( jk_spp_rhg )
807	cstr = 'ICE RHEOLOGY'
808	CASE( jk_spp_relw )
809	cstr = 'RELATIVE WND'
810	CASE( jk_spp_dqdt )
811	cstr = 'SST RELAXAT.'
812	CASE( jk_spp_deds )
813	cstr = 'SSS RELAXAT.'
814	CASE( jk_spp_arnf )
815	cstr = 'RIVER MIXING'
816	CASE( jk_spp_geot )
817	cstr = 'GEOTHERM.FLX'
818	CASE( jk_spp_qsi0 )
819	cstr = 'SOLAR EXTIN.'
820	CASE( jk_spp_bfr )
821	cstr = 'BOTTOM FRICT'
822	CASE( jk_spp_aevd )
823	cstr = 'EDDY VISCDIF'
824	CASE( jk_spp_avt )
825	cstr = 'VERT. DIFFUS'
826	CASE( jk_spp_avm )
827	cstr = 'VERT. VISCOS'
828	CASE( jk_spp_tkelc )
829	cstr = 'TKE LANGMUIR'
830	CASE( jk_spp_tkedf )
831	cstr = 'TKE RN_EDIFF'
832	CASE( jk_spp_tkeds )
833	cstr = 'TKE RN_EDISS'
834	CASE( jk_spp_tkebb )
835	cstr = 'TKE RN_EBB '
836	CASE( jk_spp_tkefr )
837	cstr = 'TKE RN_EFR '
838	CASE( jk_spp_ahtu )
839	cstr = 'TRALDF AHTU '
840	CASE( jk_spp_ahtv )
841	cstr = 'TRALDF AHTV '
842	CASE( jk_spp_ahtw )
843	cstr = 'TRALDF AHTW '
844	CASE( jk_spp_ahtt )
845	cstr = 'TRALDF AHTT '
846	CASE( jk_spp_ahubbl )
847	cstr = 'BBL DIFFUSU '
848	CASE( jk_spp_ahvbbl )
849	cstr = 'BBL DIFFUSV '
850	CASE( jk_spp_ahm1 )
851	cstr = 'DYNLDF AHM1 '
852	CASE( jk_spp_ahm2 )
853	cstr = 'DYNLDF AHM2 '
854	CASE( jk_spp_ahm3 )
855	cstr = 'DYNLDF AHM3 '
856	CASE( jk_spp_ahm4 )
857	cstr = 'DYNLDF AHM4 '
858	CASE DEFAULT
859	CALL ctl_stop('Unrecognized SPP parameter: add it or turn off diagnostics')
860	END SELECT
861	IF ( rn_mmax(jp) .GT. 0._wp ) &
862	WRITE(numrep,'(A,A,A,D10.3)') ' Maximum of absolute values for parameter ', trim(cstr), ' = ', rn_mmax(jp)
863	ENDDO
864
865	IF(sppt_step .eq. 2 ) THEN
866	zmul = rdt
867	ELSEIF (sppt_step .eq. 0 .or. sppt_step .eq. 1 ) THEN
868	zmul = 1._wp
869	ENDIF
870	rn_mmax(jk_sppt_tem:jk_sppt_vvl) = rn_mmax(jk_sppt_tem:jk_sppt_vvl) * zmul
871
872	WRITE(numrep,'(A)')
873	WRITE(numrep,'(A)') ' SPPT part'
874	WRITE(numrep,'(A,D10.3)') ' Maximum of absolute values for TEM ', rn_mmax(jk_sppt_tem)
875	IF( rn_mmax(jk_sppt_tem) > 0.5 ) WRITE(numrep,'(A)' ) 'Larger than 0.5, might be too big '
876	WRITE(numrep,'(A,D10.3)') ' Maximum of absolute values for SAL ', rn_mmax(jk_sppt_sal)
877	IF( rn_mmax(jk_sppt_sal) > 0.2 ) WRITE(numrep,'(A)' ) 'Larger than 0.2, might be too big '
878	WRITE(numrep,'(A,D10.3)') ' Maximum of absolute values for U-VEL', rn_mmax(jk_sppt_uvl)
879	IF( rn_mmax(jk_sppt_uvl) > 0.1 ) WRITE(numrep,'(A)' ) 'Larger than 0.1, might be too big '
880	WRITE(numrep,'(A,D10.3)') ' Maximum of absolute values for V-VEL', rn_mmax(jk_sppt_vvl)
881	IF( rn_mmax(jk_sppt_vvl) > 0.1 ) WRITE(numrep,'(A)' ) 'Larger than 0.1, might be too big '
882
883	WRITE(numrep,'(A)')
884	WRITE(numrep,'(A)') ' SKEB part'
885	WRITE(numrep,'(A,A,A,D10.3)') ' Maximum of absolute values for NUM ', trim(cstr), ' = ', rn_mmax(jk_skeb_dnum)
886	WRITE(numrep,'(A)' ) ' (Perturbation from numerical dissipation)'
887	IF( rn_mmax(jk_skeb_dnum) < 0.5e-4 ) WRITE(numrep,'(A)' ) ' Smaller than 0.5e-4, might be too small'
888	IF( rn_mmax(jk_skeb_dnum) > 0.2e-2 ) WRITE(numrep,'(A)' ) ' Larger than 0.2e-2, might be too big '
889	WRITE(numrep,'(A)')
890	WRITE(numrep,'(A,A,A,D10.3)') ' Maximum of absolute values for CONV ', trim(cstr), ' = ', rn_mmax(jk_skeb_dcon)
891	WRITE(numrep,'(A)' ) ' (Perturbation from convection dissipation)'
892	IF( rn_mmax(jk_skeb_dcon) < 0.5e-4 ) WRITE(numrep,'(A)' ) ' Smaller than 0.5e-4, might be too small'
893	IF( rn_mmax(jk_skeb_dcon) > 0.2e-2 ) WRITE(numrep,'(A)' ) ' Larger than 0.2e-2, might be too big '
894	WRITE(numrep,'(A)')
895	WRITE(numrep,'(A,A,A,D10.3)') ' Maximum of absolute values for TOTAL', trim(cstr), ' = ', rn_mmax(jk_skeb_tot)
896	WRITE(numrep,'(A)' ) ' (Perturbation from total energy dissipation)'
897	IF( rn_mmax(jk_skeb_tot) < 0.5e-4 ) WRITE(numrep,'(A)' ) ' Smaller than 0.5e-4, might be too small'
898	IF( rn_mmax(jk_skeb_tot) > 0.2e-2 ) WRITE(numrep,'(A)' ) ' Larger than 0.2e-2, might be too big '
899
900	CLOSE(numrep)
901	CLOSE(numdiag)
902
903	END SUBROUTINE stopack_report
904
905	SUBROUTINE spp_aht ( kt, coeff,nn_type, rn_sd, kspp )
906	!!----------------------------------------------------------------------
907	!! * ROUTINE spp_aht *
908	!!
909	!! ** Purpose : Perturbing diffusivity parameters
910	!! As spp_gen, but specifically designed for diffusivity
911	!! where coeff can be 2D or 3D
912	!!----------------------------------------------------------------------
913	INTEGER, INTENT( in ) :: kt
914	#if defined key_traldf_c3d
915	REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj,jpk) :: coeff
916	#elif defined key_traldf_c2d
917	REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj) :: coeff
918	#elif defined key_traldf_c1d
919	REAL(wp), INTENT( inout ), DIMENSION(jpk) :: coeff
920	#else
921	REAL(wp), INTENT( inout ) :: coeff
922	#endif
923	INTEGER, INTENT( in ) :: kspp
924	INTEGER, INTENT( in ) :: nn_type
925	REAL(wp), INTENT( in ) :: rn_sd
926	REAL(wp), POINTER, DIMENSION(:,:) :: gauss
927	REAL(wp) :: zsd,xme
928	INTEGER :: jk
929
930	#if defined key_traldf_c3d \|\| key_traldf_c2d
931
932	CALL wrk_alloc(jpi,jpj,gauss)
933
934	IF( ln_spp_own_gauss ) THEN
935	gauss = gauss_n_2d_p
936	ELSE
937	gauss = rn_spp_stdev * gauss_n_2d / rn_sppt_stdev
938	ENDIF
939
940	IF( nn_type == 1 ) THEN
941	gauss = gauss * rn_sd
942	#if defined key_traldf_c3d
943	DO jk=1,jpk
944	coeff(:,:,jk) = coeff(:,:,jk) * ( 1._wp + gauss )
945	ENDDO
946	#elif defined key_traldf_c2d
947	coeff = coeff * ( 1._wp + gauss )
948	#endif
949	ELSEIF ( nn_type == 2 ) THEN
950	zsd = rn_sd
951	xme = -0.5_wpzsdzsd
952	gauss = gauss * zsd + xme
953	#if defined key_traldf_c3d
954	DO jk=1,jpk
955	coeff(:,:,jk) = exp(gauss) * coeff(:,:,jk)
956	ENDDOF
957	#elif defined key_traldf_c2d
958	coeff = exp(gauss) * coeff
959	#endif
960	ELSEIF ( nn_type == 3 ) THEN
961	zsd = rn_sd
962	xme = 0._wp
963	gauss = gauss * zsd + xme
964	#if defined key_traldf_c3d
965	DO jk=1,jpk
966	coeff(:,:,jk) = exp(gauss) * coeff(:,:,jk)
967	ENDDOF
968	#elif defined key_traldf_c2d
969	coeff = exp(gauss) * coeff
970	#endif
971	ELSE
972	CALL ctl_stop( 'spp aht wrong option')
973	ENDIF
974
975	IF( ln_stopack_diags ) THEN
976	CALL spp_stats(kt,kspp,0,coeff)
977	ENDIF
978
979	#else
980	CALL ctl_stop('key_traldf_c1d is not a valid key for STO')
981
982	#endif
983
984	CALL wrk_dealloc(jpi,jpj,gauss)
985
986	END SUBROUTINE
987
988	SUBROUTINE spp_ahm ( kt, coeff,nn_type, rn_sd, kspp )
989	!!----------------------------------------------------------------------
990	!! * ROUTINE spp_ahm *
991	!!
992	!! ** Purpose : Perturbing viscosity parameters
993	!! As spp_gen, but specifically designed for viscosity
994	!! where coeff can be 2D or 3D
995	!!----------------------------------------------------------------------
996	INTEGER, INTENT( in ) :: kt
997	#if defined key_dynldf_c3d
998	REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj,jpk) :: coeff
999	#elif defined key_dynldf_c2d
1000	REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj) :: coeff
1001	#elif defined key_dynldf_c1d
1002	REAL(wp), INTENT( inout ), DIMENSION(jpk) :: coeff
1003	#else
1004	REAL(wp), INTENT( inout ) :: coeff
1005	#endif
1006	INTEGER, INTENT( in ) :: kspp
1007	INTEGER, INTENT( in ) :: nn_type
1008	REAL(wp), INTENT( in ) :: rn_sd
1009	REAL(wp), POINTER, DIMENSION(:,:) :: gauss
1010	REAL(wp) :: zsd,xme
1011	INTEGER :: jk
1012
1013	#if defined key_dynldf_c3d \|\| key_dynldf_c2d
1014
1015	CALL wrk_alloc(jpi,jpj,gauss)
1016
1017	IF( ln_spp_own_gauss ) THEN
1018	gauss = gauss_n_2d_p
1019	ELSE
1020	gauss = rn_spp_stdev * gauss_n_2d / rn_sppt_stdev
1021	ENDIF
1022
1023	IF( nn_type == 1 ) THEN
1024	gauss = gauss * rn_sd
1025	#if defined key_dynldf_c3d
1026	DO jk=1,jpk
1027	coeff(:,:,jk) = coeff(:,:,jk) * ( 1._wp + gauss )
1028	ENDDO
1029	#elif defined key_dynldf_c2d
1030	coeff = coeff * ( 1._wp + gauss )
1031	#endif
1032	ELSEIF ( nn_type == 2 ) THEN
1033	zsd = rn_sd
1034	xme = -0.5_wpzsdzsd
1035	gauss = gauss * zsd + xme
1036	#if defined key_dynldf_c3d
1037	DO jk=1,jpk
1038	coeff(:,:,jk) = exp(gauss) * coeff(:,:,jk)
1039	ENDDO
1040	#elif defined key_dynldf_c2d
1041	coeff = exp(gauss) * coeff
1042	#endif
1043	ELSEIF ( nn_type == 3 ) THEN
1044	zsd = rn_sd
1045	xme = 0._wp
1046	gauss = gauss * zsd + xme
1047	#if defined key_dynldf_c3d
1048	DO jk=1,jpk
1049	coeff(:,:,jk) = exp(gauss) * coeff(:,:,jk)
1050	ENDDO
1051	#elif defined key_dynldf_c2d
1052	coeff = exp(gauss) * coeff
1053	#endif
1054	ELSE
1055	CALL ctl_stop( 'spp ahm wrong option')
1056	ENDIF
1057
1058	IF( ln_stopack_diags ) THEN
1059	CALL spp_stats(kt,kspp,0,coeff)
1060	ENDIF
1061
1062	CALL wrk_dealloc(jpi,jpj,gauss)
1063
1064	#else
1065	CALL ctl_stop('key_traldf_c1d is not a valid key for STO')
1066	#endif
1067
1068	END SUBROUTINE
1069
1070	SUBROUTINE tra_sppt_apply ( kt , ks)
1071	!!----------------------------------------------------------------------
1072	!! * ROUTINE tra_sppt_apply *
1073	!!
1074	!! ** Purpose : Apply perturbation to tracers
1075	!! Step 0/1 for tendencies, step 2 for variables
1076	!! after timestepping
1077	!!----------------------------------------------------------------------
1078
1079	INTEGER, INTENT( in ) :: kt, ks ! Main time step counter
1080	REAL(wp) :: mi,ma
1081
1082	IF( ks .NE. sppt_step ) RETURN
1083
1084	IF(lwp) THEN
1085	WRITE(numout,*)
1086	WRITE(numout,*) ' Applying sppt on tracers at timestep ', kt
1087	WRITE(numout,*)
1088	ENDIF
1089
1090	! Modify the tendencies
1091	IF(sppt_step .eq. 2 ) THEN
1092	! At the end of the timestepping, after array swapping
1093	tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) + gauss_n * spptt * rdt
1094	tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) + gauss_n * sppts * rdt
1095	tsb(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) + gauss_n * spptt * rdt
1096	tsb(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) + gauss_n * sppts * rdt
1097	IF( ln_sppt_glocon ) CALL sppt_glocon( tsn(:,:,:,jp_tem), 'T' )
1098	IF( ln_sppt_glocon ) CALL sppt_glocon( tsn(:,:,:,jp_sal), 'T' )
1099	IF( ln_sppt_glocon ) CALL sppt_glocon( tsb(:,:,:,jp_tem), 'T' )
1100	IF( ln_sppt_glocon ) CALL sppt_glocon( tsb(:,:,:,jp_sal), 'T' )
1101	CALL lbc_lnk( tsb(:,:,:,jp_tem) , 'T', 1._wp)
1102	CALL lbc_lnk( tsb(:,:,:,jp_sal) , 'T', 1._wp)
1103	CALL lbc_lnk( tsn(:,:,:,jp_tem) , 'T', 1._wp)
1104	CALL lbc_lnk( tsn(:,:,:,jp_sal) , 'T', 1._wp)
1105	ELSEIF (sppt_step .eq. 0 .or. sppt_step .eq. 1 ) THEN
1106	! At the beginning / before vertical diffusion
1107	tsa(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) + gauss_n * spptt
1108	tsa(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) + gauss_n * sppts
1109	IF( ln_sppt_glocon ) CALL sppt_glocon( tsa(:,:,:,jp_tem), 'T' )
1110	IF( ln_sppt_glocon ) CALL sppt_glocon( tsa(:,:,:,jp_sal), 'T' )
1111	CALL lbc_lnk( tsa(:,:,:,jp_tem) , 'T', 1._wp)
1112	CALL lbc_lnk( tsa(:,:,:,jp_sal) , 'T', 1._wp)
1113	ENDIF
1114
1115	IF( ln_stopack_diags ) THEN
1116	mi = MINVAL(spptt)
1117	ma = MAXVAL(spptt)
1118	IF(lk_mpp) CALL mpp_min(mi)
1119	IF(lk_mpp) CALL mpp_max(ma)
1120	IF(lwp) WRITE(numdiag,9301) lkt,'SPPT TEMPERATURE' ,mi,ma
1121	rn_mmax ( jk_sppt_tem ) = MAX ( MAX( abs(mi), ma), rn_mmax ( jk_sppt_tem ) )
1122
1123	mi = MINVAL(sppts)
1124	ma = MAXVAL(sppts)
1125	IF(lk_mpp) CALL mpp_min(mi)
1126	IF(lk_mpp) CALL mpp_max(ma)
1127	IF(lwp) WRITE(numdiag,9301) lkt,'SPPT SALINITY ' ,mi,ma
1128	9301 FORMAT(' it :', i8, ' ', A16, ' Min: ',d10.3 , ' Max: ',d10.3)
1129	rn_mmax ( jk_sppt_sal ) = MAX ( MAX( abs(mi), ma), rn_mmax ( jk_sppt_sal ) )
1130	ENDIF
1131	! Reset the tendencies
1132	spptt = 0._wp ; sppts = 0._wp
1133
1134	END SUBROUTINE tra_sppt_apply
1135
1136	SUBROUTINE dyn_sppt_apply ( kt , ks)
1137	!!----------------------------------------------------------------------
1138	!! * ROUTINE dyn_sppt_apply *
1139	!!
1140	!! ** Purpose : Apply perturbation to momentum
1141	!! Step 0/1 for tendencies, step 2 for variables
1142	!! after timestepping
1143	!!----------------------------------------------------------------------
1144
1145	INTEGER, INTENT( in ) :: kt, ks ! Main time step counter
1146	REAL(wp) :: mi,ma
1147
1148	IF( ks .NE. sppt_step ) RETURN
1149
1150	IF(lwp) THEN
1151	WRITE(numout,*)
1152	WRITE(numout,*) ' Applying sppt on currents at timestep ', kt
1153	WRITE(numout,*)
1154	ENDIF
1155
1156	! Modify the tendencies
1157	IF(sppt_step .eq. 2 ) THEN
1158	! At the end of the timestepping, after array swapping
1159	un(:,:,:) = un(:,:,:) + gauss_nu* spptu * rdt
1160	vn(:,:,:) = vn(:,:,:) + gauss_nv* spptv * rdt
1161	ub(:,:,:) = ub(:,:,:) + gauss_nu* spptu * rdt
1162	vb(:,:,:) = vb(:,:,:) + gauss_nv* spptv * rdt
1163	IF( ln_sppt_glocon ) CALL sppt_glocon( ub, 'U' )
1164	IF( ln_sppt_glocon ) CALL sppt_glocon( vb, 'V' )
1165	IF( ln_sppt_glocon ) CALL sppt_glocon( un, 'U' )
1166	IF( ln_sppt_glocon ) CALL sppt_glocon( vn, 'V' )
1167	CALL lbc_lnk( un , 'U', -1._wp)
1168	CALL lbc_lnk( vn , 'V', -1._wp)
1169	CALL lbc_lnk( ub , 'U', -1._wp)
1170	CALL lbc_lnk( vb , 'V', -1._wp)
1171	ELSEIF (sppt_step .eq. 0 .or. sppt_step .eq. 1 ) THEN
1172	! At the beginning / before vertical diffusion
1173	ua(:,:,:) = ua(:,:,:) + gauss_nu* spptu
1174	va(:,:,:) = va(:,:,:) + gauss_nv* spptv
1175	IF( ln_sppt_glocon ) CALL sppt_glocon( ua, 'U' )
1176	IF( ln_sppt_glocon ) CALL sppt_glocon( va, 'V' )
1177	CALL lbc_lnk( ua , 'U', -1._wp)
1178	CALL lbc_lnk( va , 'V', -1._wp)
1179	ENDIF
1180
1181	IF( ln_stopack_diags ) THEN
1182	mi = MINVAL(spptu)
1183	ma = MAXVAL(spptu)
1184	IF(lk_mpp) CALL mpp_min(mi)
1185	IF(lk_mpp) CALL mpp_max(ma)
1186	IF(lwp) WRITE(numdiag,9301) lkt,'SPPT ZONAL CURR.' ,mi,ma
1187	rn_mmax ( jk_sppt_uvl ) = MAX ( MAX( abs(mi), ma), rn_mmax ( jk_sppt_uvl ) )
1188
1189	mi = MINVAL(spptv)
1190	ma = MAXVAL(spptv)
1191	IF(lk_mpp) CALL mpp_min(mi)
1192	IF(lk_mpp) CALL mpp_max(ma)
1193	IF(lwp) WRITE(numdiag,9301) lkt,'SPPT MERID CURR.' ,mi,ma
1194	9301 FORMAT(' it :', i8, ' ', A16, ' Min: ',d10.3 , ' Max: ',d10.3)
1195	rn_mmax ( jk_sppt_vvl ) = MAX ( MAX( abs(mi), ma), rn_mmax ( jk_sppt_vvl ) )
1196	ENDIF
1197
1198	! Reset the tendencies
1199	spptu = 0._wp ; spptv = 0._wp
1200
1201	END SUBROUTINE dyn_sppt_apply
1202
1203	SUBROUTINE sppt_glocon ( zts, cd_type )
1204	!!----------------------------------------------------------------------
1205	!! * ROUTINE sppt_glocon *
1206	!!
1207	!! ** Purpose : Apply global conservation constraint
1208	!!----------------------------------------------------------------------
1209	REAL(wp), INTENT(INOUT) :: zts(jpi,jpj,jpk)
1210	CHARACTER(len=1), INTENT(IN) :: cd_type
1211	INTEGER :: ji,jj,jk
1212	REAL(wp) :: zv, vl
1213	! Calculate volume tendencies and renormalize
1214	! Note: sshn should be staggered before being used.
1215	SELECT CASE ( cd_type )
1216	CASE ( 'T' )
1217	jk=1
1218	zv = SUM( tmask_i(:,:)tmask(:,:,jk)e1t(:,:)e2t(:,:)sshn(:,:)*zts(:,:,1) )
1219	vl = SUM( tmask_i(:,:)tmask(:,:,jk)e1t(:,:)e2t(:,:)sshn(:,:) )
1220	DO jk = 1, jpkm1
1221	DO jj=1,jpj
1222	DO ji=1,jpi
1223	zv = zv + SUM( tmask_i(:,:)tmask(:,:,jk)e1t(:,:)e2t(:,:)fse3t(:,:,jk)*zts(:,:,jk) )
1224	vl = vl + SUM( tmask_i(:,:)tmask(:,:,jk)e1t(:,:)e2t(:,:)fse3t(:,:,jk) )
1225	END DO
1226	END DO
1227	END DO
1228	IF(lk_mpp) CALL mpp_sum(zv)
1229	IF(lk_mpp) CALL mpp_sum(vl)
1230	zv = zv / vl
1231	zts = zts - zv*tmask
1232	CASE ( 'U' )
1233	jk=1
1234	#if defined NEMO_V34
1235	zv = SUM( umask(:,:,jk)e1u(:,:)e2u(:,:)sshn(:,:)zts(:,:,1) )
1236	vl = SUM( umask(:,:,jk)e1u(:,:)e2u(:,:)*sshn(:,:) )
1237	#else
1238	zv = SUM( umask_i(:,:)umask(:,:,jk)e1u(:,:)e2u(:,:)sshn(:,:)*zts(:,:,1) )
1239	vl = SUM( umask_i(:,:)umask(:,:,jk)e1u(:,:)e2u(:,:)sshn(:,:) )
1240	#endif
1241	DO jk = 1, jpkm1
1242	DO jj=1,jpj
1243	DO ji=1,jpi
1244	#if defined NEMO_V34
1245	zv = zv + SUM( umask(:,:,jk)e1u(:,:)e2u(:,:)fse3u(:,:,jk)zts(:,:,jk) )
1246	vl = vl + SUM( umask(:,:,jk)e1u(:,:)e2u(:,:)*fse3u(:,:,jk) )
1247	#else
1248	zv = zv + SUM( umask_i(:,:)umask(:,:,jk)e1u(:,:)e2u(:,:)fse3u(:,:,jk)*zts(:,:,jk) )
1249	vl = vl + SUM( umask_i(:,:)umask(:,:,jk)e1u(:,:)e2u(:,:)fse3u(:,:,jk) )
1250	#endif
1251	END DO
1252	END DO
1253	END DO
1254	IF(lk_mpp) CALL mpp_sum(zv)
1255	IF(lk_mpp) CALL mpp_sum(vl)
1256	zv = zv / vl
1257	zts = zts - zv*umask
1258	CASE ( 'V' )
1259	jk=1
1260	#if defined NEMO_V34
1261	zv = SUM( vmask(:,:,jk)e1v(:,:)e2v(:,:)sshn(:,:)zts(:,:,1) )
1262	vl = SUM( vmask(:,:,jk)e1v(:,:)e2v(:,:)*sshn(:,:) )
1263	#else
1264	zv = SUM( vmask_i(:,:)vmask(:,:,jk)e1v(:,:)e2v(:,:)sshn(:,:)*zts(:,:,1) )
1265	vl = SUM( vmask_i(:,:)vmask(:,:,jk)e1v(:,:)e2v(:,:)sshn(:,:) )
1266	#endif
1267	DO jk = 1, jpkm1
1268	DO jj=1,jpj
1269	DO ji=1,jpi
1270	#if defined NEMO_V34
1271	zv = zv + SUM( vmask(:,:,jk)e1v(:,:)e2v(:,:)fse3v(:,:,jk)zts(:,:,jk) )
1272	vl = vl + SUM( vmask(:,:,jk)e1v(:,:)e2v(:,:)*fse3v(:,:,jk) )
1273	#else
1274	zv = zv + SUM( vmask_i(:,:)vmask(:,:,jk)e1v(:,:)e2v(:,:)fse3v(:,:,jk)*zts(:,:,jk) )
1275	vl = vl + SUM( vmask_i(:,:)vmask(:,:,jk)e1v(:,:)e2v(:,:)fse3v(:,:,jk) )
1276	#endif
1277	END DO
1278	END DO
1279	END DO
1280	IF(lk_mpp) CALL mpp_sum(zv)
1281	IF(lk_mpp) CALL mpp_sum(vl)
1282	zv = zv / vl
1283	zts = zts - zv*vmask
1284	ENDSELECT
1285
1286	END SUBROUTINE sppt_glocon
1287
1288	SUBROUTINE gaussian_ar1_field ( gn, nk, wei, gb, a, b,gn0, fltf, istep)
1289	!!----------------------------------------------------------------------
1290	!! * ROUTINE gaussian_ar1_field *
1291	!!
1292	!! ** Purpose : Generate correlated perturbation field
1293	!!----------------------------------------------------------------------
1294	REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: a, b
1295	REAL(wp), DIMENSION(jpk) , INTENT(INOUT) :: wei
1296	REAL(wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: gb
1297	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(INOUT) :: gn
1298	REAL(wp), DIMENSION(jpi,jpj ), INTENT(OUT) :: gn0
1299	INTEGER, INTENT(IN) :: nk,istep
1300	REAL(wp), INTENT(IN) :: fltf
1301	REAL(wp) :: g2d(jpi,jpj)
1302	INTEGER :: jf, jk, ji, jj, nbl
1303
1304	! Random noise on 2d field
1305	IF ( istep == 1 ) THEN
1306	CALL sppt_rand2d( g2d )
1307	CALL lbc_lnk( g2d , 'T', 1._wp)
1308	g2d_save = g2d
1309	IF( nn_sppt_step_bound .EQ. 1 ) THEN
1310	WHERE ( g2d < -ABS(rn_sppt_bound) ) g2d = -ABS(rn_sppt_bound)
1311	WHERE ( g2d > ABS(rn_sppt_bound) ) g2d = ABS(rn_sppt_bound)
1312	ENDIF
1313	ELSEIF ( istep == 2 ) THEN
1314	g2d = rn_spp_stdev * g2d_save / rn_sppt_stdev
1315	ELSEIF ( istep == 3 ) THEN
1316	g2d = rn_skeb_stdev * g2d_save / rn_sppt_stdev
1317	ENDIF
1318
1319	! Laplacian filter and re-normalization
1320	DO jf = 1, nk
1321	CALL sppt_flt( g2d )
1322	END DO
1323	g2d = g2d * fltf
1324
1325	#ifdef key_iommput
1326	! Output the random field
1327	IF(istep==1) THEN
1328	CALL iom_put( "sppt_ran" , g2d )
1329	ELSEIF (istep==2) THEN
1330	CALL iom_put( "spp_ran" , g2d )
1331	ELSEIF (istep==3) THEN
1332	CALL iom_put( "skeb_ran" , g2d )
1333	ENDIF
1334	#endif
1335
1336	! AR(1) process and array swap
1337	g2d = agb + bg2d
1338	gb = g2d
1339	gn0 = g2d * zdc(:,:,1)
1340
1341	IF (istep==2 .or. istep==3) RETURN
1342
1343	! From 2- to 3-d and vertical weigth
1344	IF(nn_vwei .eq. 2) THEN
1345	DO jj=1,jpj
1346	DO ji=1,jpi
1347	nbl = mbathy(ji,jj)
1348	wei(:) = 0._wp
1349	IF(nbl>0) THEN
1350	wei(1:nbl) = 1._wp
1351	wei(1) = 0._wp
1352	wei(2) = 0.5_wp
1353	wei(nbl-1) = 0.5_wp
1354	wei(nbl) = 0._wp
1355	DO jk=1,jpk
1356	gn(ji,jj,jk) = g2d(ji,jj) * wei(jk) * zdc(ji,jj,jk)
1357	ENDDO
1358	ENDIF
1359	ENDDO
1360	ENDDO
1361	ELSEIF(nn_vwei .eq. 3) THEN
1362	DO jj=1,jpj
1363	DO ji=1,jpi
1364	nbl = mbathy(ji,jj)
1365	wei(:) = 0._wp
1366	IF(nbl>0) THEN
1367	wei(1:nbl) = 1._wp
1368	wei(nbl-1) = 0.5_wp
1369	wei(nbl) = 0._wp
1370	DO jk=1,jpk
1371	gn(ji,jj,jk) = g2d(ji,jj) * wei(jk) * zdc(ji,jj,jk)
1372	ENDDO
1373	ENDIF
1374	ENDDO
1375	ENDDO
1376	ELSE
1377	DO jk=1,jpk
1378	gn(:,:,jk) = g2d * wei(jk) * zdc(:,:,jk)
1379	ENDDO
1380	ENDIF
1381
1382	! Bound
1383	IF( nn_sppt_step_bound .EQ. 2 ) THEN
1384	WHERE ( gn < -ABS(rn_sppt_bound) ) gn = -ABS(rn_sppt_bound)
1385	WHERE ( gn > ABS(rn_sppt_bound) ) gn = ABS(rn_sppt_bound)
1386	WHERE ( gn0< -ABS(rn_sppt_bound) ) gn0= -ABS(rn_sppt_bound)
1387	WHERE ( gn0> ABS(rn_sppt_bound) ) gn0= ABS(rn_sppt_bound)
1388	ENDIF
1389
1390	END SUBROUTINE gaussian_ar1_field
1391	!
1392	SUBROUTINE stopack_init
1393	!!----------------------------------------------------------------------
1394	!! * ROUTINE stopack_init *
1395	!!
1396	!! ** Purpose : Initialize stopack
1397	!!----------------------------------------------------------------------
1398	!!
1399	INTEGER :: ios, inum
1400	INTEGER :: nn_sppt_tra, nn_sppt_dyn, nn_spp_aht, nn_sppt_ice
1401	INTEGER :: ivals(8)
1402	INTEGER(KIND=8) :: ziseed(4) ! RNG seeds in integer type
1403	REAL(KIND=8) :: zrseed(4) ! RNG seeds in real type
1404	REAL(KIND=8) :: rdt_sppt
1405	!!
1406	NAMELIST/namstopack/ ln_stopack, ln_sppt_taumap, rn_sppt_tau, &
1407	ln_stopack_restart, sppt_filter_pass, rn_sppt_bound, sppt_step, nn_deftau,rn_sppt_stdev,&
1408	ln_distcoast, rn_distcoast, nn_vwei, nn_stopack_seed, nn_sppt_step_bound, nn_rndm_freq, &
1409	ln_sppt_glocon, ln_stopack_repr, &
1410	rn_icestr_sd, rn_icealb_sd, &
1411	ln_sppt_traxad, ln_sppt_trayad, ln_sppt_trazad, ln_sppt_trasad, ln_sppt_traldf, &
1412	ln_sppt_trazdf, ln_sppt_trazdfp,ln_sppt_traevd, ln_sppt_trabbc, ln_sppt_trabbl, &
1413	ln_sppt_tranpc, ln_sppt_tradmp, ln_sppt_traqsr, ln_sppt_transr, ln_sppt_traatf ,&
1414	ln_sppt_dynhpg, ln_sppt_dynspg, ln_sppt_dynkeg, ln_sppt_dynrvo, ln_sppt_dynpvo, &
1415	ln_sppt_dynzad, ln_sppt_dynldf, ln_sppt_dynzdf, ln_sppt_dynbfr, ln_sppt_dynatf, ln_sppt_dyntau,&
1416	ln_sppt_icehdf, ln_sppt_icelat, ln_sppt_icezdf, & !ln_sppt_icealb, ln_sppt_icestr
1417	nn_spp_icealb, nn_spp_icestr,&
1418	spp_filter_pass,rn_spp_stdev,rn_spp_tau,&
1419	nn_spp_bfr,nn_spp_dqdt,nn_spp_dedt,nn_spp_avt,nn_spp_avm,nn_spp_qsi0,nn_spp_ahtu,nn_spp_ahtv,nn_spp_ahm1,nn_spp_ahm2,&
1420	nn_spp_ahtw,nn_spp_ahtt,nn_spp_relw,nn_spp_arnf,nn_spp_geot,nn_spp_aevd,nn_spp_ahubbl,nn_spp_ahvbbl,&
1421	nn_spp_tkelc,nn_spp_tkedf,nn_spp_tkeds,nn_spp_tkebb,nn_spp_tkefr,&
1422	rn_bfr_sd,rn_dqdt_sd,rn_dedt_sd,rn_avt_sd,rn_avm_sd,rn_qsi0_sd,rn_ahtu_sd,rn_ahtv_sd,rn_ahm1_sd,rn_ahm2_sd,&
1423	rn_ahtw_sd,rn_ahtt_sd, rn_relw_sd, rn_arnf_sd,rn_geot_sd, rn_aevd_sd,rn_ahubbl_sd,rn_ahvbbl_sd,&
1424	rn_tkelc_sd,rn_tkedf_sd,rn_tkeds_sd,rn_tkebb_sd,rn_tkefr_sd,&
1425	ln_skeb,rn_skeb,nn_dcom_freq,rn_kh,rn_kc,ln_stopack_diags,skeb_filter_pass,rn_skeb_stdev,rn_skeb_tau,&
1426	nn_dconv,rn_beta_num, rn_beta_con
1427
1428	! Default values
1429	rn_sppt_bound = 3._wp
1430	ln_stopack = .false.
1431	ln_sppt_taumap = .false.
1432	rn_sppt_tau = 86400._wp * 5._wp
1433	sppt_filter_pass = 30
1434	nn_vwei = 0
1435	ln_distcoast = .false.
1436	ln_sppt_glocon = .false.
1437	rn_distcoast = 10._wp
1438	ln_stopack_restart = .true.
1439	nn_stopack_seed = (/1,2,3,narea/)
1440	nn_rndm_freq = 1
1441	nn_sppt_step_bound = 2
1442	nn_deftau = 1
1443	ln_sppt_traxad = .false. ! Switch for x advection
1444	ln_sppt_trayad = .false. ! Switch for y advection
1445	ln_sppt_trazad = .false. ! Switch for z advection
1446	ln_sppt_trasad = .false. ! Switch for z advection (s- case)
1447	ln_sppt_traldf = .false. ! Switch for lateral diffusion
1448	ln_sppt_trazdf = .false. ! Switch for vertical diffusion
1449	ln_sppt_trazdfp = .false. ! Switch for pure vertical diffusion
1450	ln_sppt_traevd = .false. ! Switch for enhanced vertical diffusion
1451	ln_sppt_trabbc = .false. ! Switch for bottom boundary condition
1452	ln_sppt_trabbl = .false. ! Switch for bottom boundary layer
1453	ln_sppt_tranpc = .false. ! Switch for non-penetrative convection
1454	ln_sppt_tradmp = .false. ! Switch for tracer damping
1455	ln_sppt_traqsr = .false. ! Switch for solar radiation
1456	ln_sppt_transr = .false. ! Switch for non-solar radiation / freshwater flux
1457	ln_sppt_traatf = .false. ! Switch for Asselin time-filter
1458
1459	ln_sppt_dynhpg = .false. ! Switch for hydrost. press. grad.
1460	ln_sppt_dynspg = .false. ! Switch for surface press. grad.
1461	ln_sppt_dynkeg = .false. ! Switch for horiz. advcetion
1462	ln_sppt_dynrvo = .false. ! Switch for Relative vorticity
1463	ln_sppt_dynpvo = .false. ! Switch for planetary vortic.
1464	ln_sppt_dynzad = .false. ! Switch for vertical advection
1465	ln_sppt_dynldf = .false. ! Switch for lateral viscosity
1466	ln_sppt_dynzdf = .false. ! Switch for vertical viscosity
1467	ln_sppt_dynbfr = .false. ! Switch for bottom friction
1468	ln_sppt_dynatf = .false. ! Switch for Asselin filter
1469	ln_sppt_dyntau = .false. ! Switch for wind stress
1470
1471	ln_sppt_icehdf = .false.
1472	ln_sppt_icelat = .false.
1473	ln_sppt_icezdf = .false.
1474
1475	nn_spp_icestr = 0
1476	nn_spp_icealb = 0
1477	rn_icestr_sd = 0.30_wp
1478	rn_icealb_sd = 0.30_wp
1479
1480	nn_spp_bfr =0
1481	nn_spp_dqdt =0
1482	nn_spp_arnf =0
1483	nn_spp_aevd =0
1484	nn_spp_geot =0
1485	nn_spp_dedt =0
1486	nn_spp_avt =0
1487	nn_spp_avm =0
1488	nn_spp_qsi0 =0
1489	nn_spp_relw =0
1490	nn_spp_tkelc =0
1491	nn_spp_tkedf =0
1492	nn_spp_tkeds =0
1493	nn_spp_tkebb =0
1494	nn_spp_tkefr =0
1495
1496	ln_skeb = .false.
1497	ln_stopack_diags = .false.
1498	rn_skeb_stdev = 1.0_wp
1499	skeb_filter_pass = 50
1500	rn_skeb_tau = 50
1501
1502	#ifdef NEMO_V34
1503	REWIND( numnam )
1504	READ ( numnam, namstopack )
1505	#else
1506	REWIND( numnam_ref )
1507	READ ( numnam_ref, namstopack, IOSTAT = ios, ERR = 901)
1508	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namstopack in reference namelist', lwp )
1509
1510	REWIND( numnam_cfg )
1511	READ ( numnam_cfg, namstopack, IOSTAT = ios, ERR = 902 )
1512	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namstopack in configuration namelist', lwp )
1513	IF(lwm .AND. nprint > 2) WRITE ( numond, namstopack )
1514	#endif
1515
1516	IF(lwp) THEN
1517	WRITE(numout,*)
1518	WRITE(numout,*) 'init_stopack : Stochastic physics package'
1519	WRITE(numout,*) '~~~~~~~~~~~~'
1520	WRITE(numout,*) ' Namelist namstopack: '
1521	WRITE(numout,*)
1522	WRITE(numout,*) ' Switch on STOPACK ln_stopack : ',ln_stopack
1523	WRITE(numout,*) ' Verbose diagnostics for STOPACK ln_stopack_diags:', ln_stopack_diags
1524	WRITE(numout,*)
1525	WRITE(numout,) ' ** Random generation'
1526	WRITE(numout,*) ' Read from restart file previous perturbation ln_stopack_restart :', ln_stopack_restart
1527	WRITE(numout,*) ' Frequency of calls of the SPPT perturbation field nn_rndm_freq :', nn_rndm_freq
1528	WRITE(numout,*) ' Reproducibility of random generation ln_stopack_repr :', ln_stopack_repr
1529	WRITE(numout,*) ' Seed for random number generator (no restart) nn_stopack_seed :', nn_stopack_seed(:)
1530	WRITE(numout,*)
1531	WRITE(numout,) ' ** SPPT scheme '
1532	WRITE(numout,*) ' SPPT step (0: beginning; 1: before ZVD; 2: after ZVD) sppt_step : ',sppt_step
1533	WRITE(numout,*) ' Use map of decorr. time scale ln_sppt_taumap :', ln_sppt_taumap
1534	WRITE(numout,*) ' If ln_sppt_taumap FALSE, use this constant (in days) rn_sppt_tau :', rn_sppt_tau
1535	WRITE(numout,*) ' Number of filter passes to correlate in space sppt_filter_pass :', sppt_filter_pass
1536	WRITE(numout,*) ' Standard deviation of the white noise rn_sppt_stdev :', rn_sppt_stdev
1537	WRITE(numout,*) ' Apply global conservation constraints ln_sppt_glocon :', ln_sppt_glocon
1538	WRITE(numout,*)
1539	WRITE(numout,*) ' Perturbation on tracers:'
1540	WRITE(numout,*) ' Switch for x advection ln_sppt_traxad :', ln_sppt_traxad
1541	WRITE(numout,*) ' Switch for y advection ln_sppt_trayad :', ln_sppt_trayad
1542	WRITE(numout,*) ' Switch for z advection ln_sppt_trazad :', ln_sppt_trazad
1543	WRITE(numout,*) ' Switch for z advection (s- case) ln_sppt_trasad :', ln_sppt_trasad
1544	WRITE(numout,*) ' Switch for lateral diffusion ln_sppt_traldf :', ln_sppt_traldf
1545	WRITE(numout,*) ' Switch for vertical diffusion ln_sppt_trazdf :', ln_sppt_trazdf
1546	WRITE(numout,*) ' Switch for pure vertical diffusion ln_sppt_trazdfp :', ln_sppt_trazdfp
1547	WRITE(numout,*) ' Switch for enhanced vertical diffusion ln_sppt_traevd :', ln_sppt_traevd
1548	WRITE(numout,*) ' Switch for bottom boundary condition ln_sppt_trabbc :', ln_sppt_trabbc
1549	WRITE(numout,*) ' Switch for bottom boundary layer ln_sppt_trabbl :', ln_sppt_trabbl
1550	WRITE(numout,*) ' Switch for non-penetrative convection ln_sppt_tranpc :', ln_sppt_tranpc
1551	WRITE(numout,*) ' Switch for tracer damping ln_sppt_tradmp :', ln_sppt_tradmp
1552	WRITE(numout,*) ' Switch for solar radiation ln_sppt_traqsr :', ln_sppt_traqsr
1553	WRITE(numout,*) ' Switch for non-solar rad. / freshwater flx flux ln_sppt_transr :', ln_sppt_transr
1554	WRITE(numout,*) ' Switch for Asselin time-filter ln_sppt_traatf :', ln_sppt_traatf
1555	WRITE(numout,*)
1556	WRITE(numout,*) ' Perturbation on dynamics:'
1557	WRITE(numout,*) ' Switch for horiz advection ln_sppt_dynkeg :', ln_sppt_dynkeg
1558	WRITE(numout,*) ' Switch for z advection ln_sppt_dynzad :', ln_sppt_dynzad
1559	WRITE(numout,*) ' Switch for wind stress ln_sppt_dyntau :', ln_sppt_dyntau
1560	WRITE(numout,*) ' Switch for lateral diffusion ln_sppt_dynldf :', ln_sppt_dynldf
1561	WRITE(numout,*) ' Switch for vertical diffusion ln_sppt_dynzdf :', ln_sppt_dynzdf
1562	WRITE(numout,*) ' Switch for planet. vorticity ln_sppt_dynpvo :', ln_sppt_dynpvo
1563	WRITE(numout,*) ' Switch for relative vorticity ln_sppt_dynrvo :', ln_sppt_dynrvo
1564	WRITE(numout,*) ' Switch for hydrost. press. grad. ln_sppt_dynhpg :', ln_sppt_dynhpg
1565	WRITE(numout,*) ' Switch for surface press. grad. ln_sppt_dynspg :', ln_sppt_dynspg
1566	WRITE(numout,*) ' Switch for bottom friction ln_sppt_dynbfr :', ln_sppt_dynbfr
1567	WRITE(numout,*) ' Switch for Asselin time-filter ln_sppt_dynatf :', ln_sppt_dynatf
1568	WRITE(numout,*)
1569	WRITE(numout,*) ' Perturbation on sea-ice:'
1570	WRITE(numout,*) ' Switch for sea-ice diffusivity ln_sppt_icehdf :', ln_sppt_icehdf
1571	WRITE(numout,*) ' Switch for sea-ice lateral accretion ln_sppt_icelat :', ln_sppt_icelat
1572	WRITE(numout,*) ' Switch for sea-ice vertical thermodyn. ln_sppt_icezdf :', ln_sppt_icezdf
1573	WRITE(numout,*)
1574	WRITE(numout,*) ' Bound for gaussian random numbers rn_sppt_bound :', rn_sppt_bound
1575	WRITE(numout,*) ' Bound Step (0: nobound; 1: Gaussian; 2: Pert) nn_sppt_step_bound :', nn_sppt_step_bound
1576	WRITE(numout,*) ' Definition of Tau (1: days; 2: timesteps) nn_deftau :', nn_deftau
1577	WRITE(numout,*) ' Smoothing of perturbation close to coast ln_distcoast :', ln_distcoast
1578	WRITE(numout,*) ' Spatial scale of the smoothing near coasts (m) rn_distcoast :', rn_distcoast
1579	WRITE(numout,*) ' Type of vertical weight: nn_vwei :', nn_vwei
1580	WRITE(numout,*) ' 0 : No weight '
1581	WRITE(numout,*) ' 1 : First/Last level smoothing '
1582	WRITE(numout,*) ' 2 : Top/Bottom smoothing'
1583	WRITE(numout,*) ' 3 : Bottom smoothing'
1584	WRITE(numout,*)
1585	WRITE(numout,) ' ** SPP schemes (0: off ; 1: normal; 2:lognormal, mean as unpert.; 3: lognormal, median as unpert.'
1586	WRITE(numout,*) ' (the meaning of standard dev. depends on the distribution and parametr.)'
1587	WRITE(numout,*)
1588	WRITE(numout,*) ' Number of passes for spatial filter (AR1 field) spp_filter_pass:', spp_filter_pass
1589	WRITE(numout,*) ' Standard deviation of random generator (AR1 field) rn_spp_stdev :', rn_spp_stdev
1590	WRITE(numout,*) ' Decorr. time scale (AR1 field) rn_spp_tau :', rn_spp_tau
1591	WRITE(numout,*)
1592	WRITE(numout,*) ' SPP for bottom friction coeff nn_spp_bfr :', nn_spp_bfr
1593	WRITE(numout,*) ' STDEV rn_bfr_sd :', rn_bfr_sd
1594	WRITE(numout,*) ' SPP for SST relaxation coeff nn_spp_dqdt :', nn_spp_dqdt
1595	WRITE(numout,*) ' STDEV rn_dqdt_sd :', rn_dqdt_sd
1596	WRITE(numout,*) ' SPP for SSS relaxation coeff nn_spp_dedt :', nn_spp_dedt
1597	WRITE(numout,*) ' STDEV rn_dedt_sd :', rn_dedt_sd
1598	WRITE(numout,*) ' SPP for vertical tra mixing coeff (only TKE, GLS) nn_spp_avt :', nn_spp_avt
1599	WRITE(numout,*) ' STDEV rn_avt_sd :', rn_avt_sd
1600	WRITE(numout,*) ' SPP for vertical dyn mixing coeff (only TKE, GLS) nn_spp_avm :', nn_spp_avm
1601	WRITE(numout,*) ' STDEV rn_avm_sd :', rn_avm_sd
1602	WRITE(numout,*) ' SPP for solar penetration scheme (only RGB) nn_spp_qsi0 :', nn_spp_qsi0
1603	WRITE(numout,*) ' STDEV rn_qsi0_sd :', rn_qsi0_sd
1604	WRITE(numout,*) ' SPP for horiz. diffusivity U nn_spp_ahtu :', nn_spp_ahtu
1605	WRITE(numout,*) ' STDEV rn_ahtu_sd :', rn_ahtu_sd
1606	WRITE(numout,*) ' SPP for horiz. diffusivity V nn_spp_ahtv :', nn_spp_ahtv
1607	WRITE(numout,*) ' STDEV rn_ahtv_sd :', rn_ahtv_sd
1608	WRITE(numout,*) ' SPP for horiz. diffusivity W nn_spp_ahtw :', nn_spp_ahtw
1609	WRITE(numout,*) ' STDEV rn_ahtw_sd :', rn_ahtw_sd
1610	WRITE(numout,*) ' SPP for horiz. diffusivity T nn_spp_ahtt :', nn_spp_ahtt
1611	WRITE(numout,*) ' STDEV rn_ahtt_sd :', rn_ahtt_sd
1612	WRITE(numout,*) ' SPP for horiz. viscosity (1/3) nn_spp_ahm1 :', nn_spp_ahm1
1613	WRITE(numout,*) ' STDEV rn_ahm1_sd :', rn_ahm1_sd
1614	WRITE(numout,*) ' SPP for horiz. viscosity (2/4) nn_spp_ahm2 :', nn_spp_ahm2
1615	WRITE(numout,*) ' STDEV rn_ahm2_sd :', rn_ahm2_sd
1616	WRITE(numout,*) ' SPP for relative wind factor nn_spp_relw :', nn_spp_relw
1617	WRITE(numout,*) ' (use 4, 5, 6 for nn_spp_relw to have options 1, 2, 3 with limits bounded to [0,1]'
1618	WRITE(numout,*) ' STDEV rn_relw_sd :', rn_relw_sd
1619	WRITE(numout,*) ' SPP for mixing close to river mouth nn_spp_arnf :', nn_spp_arnf
1620	WRITE(numout,*) ' STDEV rn_arnf_sd :', rn_arnf_sd
1621	WRITE(numout,*) ' SPP for geothermal heating nn_spp_geot :', nn_spp_geot
1622	WRITE(numout,*) ' STDEV rn_geot_sd :', rn_geot_sd
1623	WRITE(numout,*) ' SPP for enhanced vertical diffusion nn_spp_aevd :', nn_spp_aevd
1624	WRITE(numout,*) ' STDEV rn_aevd_sd :', rn_aevd_sd
1625	WRITE(numout,*) ' SPP for TKE rn_lc Langmuir cell coefficient nn_spp_tkelc :', nn_spp_tkelc
1626	WRITE(numout,*) ' STDEV rn_tkelc_sd :', rn_tkelc_sd
1627	WRITE(numout,*) ' SPP for TKE rn_ediff Eddy diff. coefficient nn_spp_tkedf :', nn_spp_tkedf
1628	WRITE(numout,*) ' STDEV rn_tkedf_sd :', rn_tkedf_sd
1629	WRITE(numout,*) ' SPP for TKE rn_ediss Kolmogoroff dissipation coeff. nn_spp_tkeds :', nn_spp_tkeds
1630	WRITE(numout,*) ' STDEV rn_tkeds_sd :', rn_tkeds_sd
1631	WRITE(numout,*) ' SPP for TKE rn_ebb Surface input of tke nn_spp_tkebb :', nn_spp_tkebb
1632	WRITE(numout,*) ' STDEV rn_tkebb_sd :', rn_tkebb_sd
1633	WRITE(numout,*) ' SPP for TKE rn_efr Fraction of srf TKE below ML nn_spp_tkefr :', nn_spp_tkefr
1634	WRITE(numout,*) ' STDEV rn_tkefr_sd :', rn_tkefr_sd
1635	WRITE(numout,*) ' SPP for BBL U diffusivity nn_spp_ahubbl:', nn_spp_ahubbl
1636	WRITE(numout,*) ' STDEV rn_ahubbl_sd :', rn_ahubbl_sd
1637	WRITE(numout,*) ' SPP for BBL V diffusivity nn_spp_ahvbbl:', nn_spp_ahvbbl
1638	WRITE(numout,*) ' STDEV rn_ahvbbl_sd :', rn_ahvbbl_sd
1639	WRITE(numout,*)
1640	WRITE(numout,) ' ** SPP schemes for sea-ice '
1641	WRITE(numout,*) ' Albedo nn_spp_icealb:', nn_spp_icealb
1642	WRITE(numout,*) ' St. dev. for ice albedo rn_icealb_sd :', rn_icealb_sd
1643	WRITE(numout,*) ' Ice Strength nn_spp_icestr:', nn_spp_icestr
1644	WRITE(numout,*) ' St. dev. for ice strength rn_icestr_sd :', rn_icestr_sd
1645	WRITE(numout,*)
1646	WRITE(numout,*) ' SKEB Perturbation scheme '
1647	WRITE(numout,*) ' SKEB switch ln_skeb :', ln_skeb
1648	WRITE(numout,*) ' SKEB ratio of backscattered energy rn_skeb :', rn_skeb
1649	WRITE(numout,*) ' Frequency update for dissipation mask nn_dcom_freq :', nn_dcom_freq
1650	WRITE(numout,*) ' Numerical dissipation factor (resolut. dependent) rn_kh :', rn_kh
1651	WRITE(numout,*) ' Number of passes for spatial filter (AR1 field) skeb_filter_pass:', skeb_filter_pass
1652	WRITE(numout,*) ' Standard deviation of random generator (AR1 field) rn_skeb_stdev:', rn_skeb_stdev
1653	WRITE(numout,*) ' Decorr. time scale (AR1 field) rn_skeb_tau :', rn_skeb_tau
1654	WRITE(numout,*) ' Option of convection energy dissipation nn_dconv :', nn_dconv
1655	WRITE(numout,*) ' Convection dissipation factor (resolut. dependent) rn_kc :', rn_kc
1656	WRITE(numout,*) ' Multiplier for numerical dissipation rn_beta_num :', rn_beta_num
1657	WRITE(numout,*) ' Multiplier for convection dissipation rn_beta_con :', rn_beta_con
1658
1659	WRITE(numout,*)
1660	IF(lflush) CALL flush(numout)
1661	ENDIF
1662
1663	IF( .NOT. ln_stopack ) THEN
1664	IF(lwp) WRITE(numout,*) ' STOPACK is switched off'
1665	IF(lflush) CALL flush(numout)
1666	RETURN
1667	ENDIF
1668
1669	IF( MOD( nitend - nit000 + 1, nn_rndm_freq) /= 0 .OR. &
1670	( MOD( nstock, nn_rndm_freq) /= 0 .AND. nstock .GT. 0 ) ) THEN
1671	WRITE(ctmp1,*) 'experiment length (', nitend - nit000 + 1, ') or nstock (', nstock, &
1672	& ' is NOT a multiple of nn_rndm_freq (', nn_rndm_freq, ')'
1673	CALL ctl_stop( ctmp1, 'Impossible to properly setup STOPACK restart' )
1674	ENDIF
1675
1676	nn_sppt_tra = COUNT( (/ln_sppt_traxad, ln_sppt_trayad, ln_sppt_trazad, &
1677	ln_sppt_trasad, ln_sppt_traldf, ln_sppt_trazdf, &
1678	ln_sppt_trazdfp, ln_sppt_traevd, ln_sppt_trabbc, &
1679	ln_sppt_trabbl, ln_sppt_tranpc, ln_sppt_tradmp, &
1680	ln_sppt_traqsr, ln_sppt_transr, ln_sppt_traatf/) )
1681	nn_sppt_dyn = COUNT( (/ln_sppt_dynhpg, ln_sppt_dynspg, ln_sppt_dynkeg, &
1682	ln_sppt_dynrvo, ln_sppt_dynpvo, ln_sppt_dynzad, &
1683	ln_sppt_dynldf, ln_sppt_dynzdf, ln_sppt_dynbfr, &
1684	ln_sppt_dynatf, ln_sppt_dyntau/) )
1685	nn_sppt_ice = COUNT( (/ln_sppt_icehdf, ln_sppt_icelat, ln_sppt_icezdf/) )
1686
1687	ln_sppt_tra = ( nn_sppt_tra > 0 )
1688	ln_sppt_dyn = ( nn_sppt_dyn > 0 )
1689	ln_sppt_ice = ( nn_sppt_ice > 0 )
1690
1691	nn_spp = nn_spp_bfr+nn_spp_dqdt+nn_spp_dedt+nn_spp_avt+nn_spp_avm+nn_spp_qsi0+&
1692	& nn_spp_ahtu+nn_spp_ahtv+nn_spp_ahtw+nn_spp_ahtt+nn_spp_relw+nn_spp_arnf+nn_spp_geot+nn_spp_aevd+&
1693	& nn_spp_tkelc+nn_spp_tkedf+nn_spp_tkeds+nn_spp_tkebb+nn_spp_tkefr+nn_spp_ahubbl+nn_spp_ahvbbl+&
1694	& nn_spp_ahm1+nn_spp_ahm2
1695
1696	#ifdef NEMO_V34
1697
1698	#ifndef key_trdtra
1699	IF( ln_sppt_tra ) &
1700	& CALL ctl_stop( ' SPPT on tracers on .AND. .NOT. key_trdtra ', &
1701	& ' SPPT cannot work without tracer tendency computation')
1702	#endif
1703	#ifndef key_trddyn
1704	IF( ln_sppt_dyn ) &
1705	& CALL ctl_stop( ' SPPT on momentum on .AND. .NOT. key_trddyn ', &
1706	& ' SPPT cannot work without dynamic tendency computation')
1707	#endif
1708
1709	#else
1710	IF( ln_sppt_tra .AND. .NOT. ln_tra_trd ) &
1711	& CALL ctl_stop( ' SPPT on tracers on .AND. .NOT. ln_tra_trd ', &
1712	& ' SPPT cannot work without tracer tendency computation')
1713
1714	IF( ln_sppt_dyn .AND. .NOT. ln_dyn_trd ) &
1715	& CALL ctl_stop( ' SPPT on momentum on .AND. .NOT. ln_dyn_trd ', &
1716	& ' SPPT cannot work without dynamic tendency computation')
1717	#endif
1718
1719	IF( ln_sppt_glocon .AND. lk_vvl ) &
1720	& CALL ctl_stop( ' ln_sppt_glocal .AND. lk_vvl ', &
1721	& ' SPPT conservation not coded yet for VVL')
1722
1723	IF( nn_deftau .NE. 1 .AND. nn_deftau .NE. 2 ) &
1724	& CALL ctl_stop( ' nn_deftau must be 1 or 2 ')
1725
1726	nn_spp_aht = nn_spp_ahtu + nn_spp_ahtv + nn_spp_ahtw + nn_spp_ahtt
1727	IF(nn_spp_aht .GT. 0) THEN
1728	#if defined key_traldf_c3d
1729	IF(lwp) WRITE(numout,*) 'SPP : diffusivity perturbation with 3D coefficients'
1730	#elif defined key_traldf_c2d
1731	IF(lwp) WRITE(numout,*) 'SPP : diffusivity perturbation with 2D coefficients'
1732	#else
1733	CALL ctl_stop( 'SPP : diffusivity perturbation requires key_traldf_c3d or key_traldf_c2d')
1734	#endif
1735	ENDIF
1736
1737	IF(nn_spp_ahm1+nn_spp_ahm2 .GT. 0) THEN
1738	#if defined key_dynldf_c3d
1739	IF(lwp) WRITE(numout,*) 'SPP : viscosity perturbation with 3D coefficients'
1740	#elif defined key_dynldf_c2d
1741	IF(lwp) WRITE(numout,*) 'SPP : viscosity perturbation with 2D coefficients'
1742	#else
1743	CALL ctl_stop( 'SPP : viscosity perturbation requires key_dynldf_c3d or key_dynldf_c2d')
1744	#endif
1745	ENDIF
1746
1747	ln_skeb_own_gauss = .FALSE.
1748	IF( ln_skeb ) THEN
1749	IF( skeb_filter_pass > 0 .AND. skeb_filter_pass .NE. sppt_filter_pass ) ln_skeb_own_gauss = .TRUE.
1750	IF( rn_skeb_tau > 0._wp .AND. rn_skeb_tau .NE. rn_sppt_tau ) ln_skeb_own_gauss = .TRUE.
1751	ENDIF
1752
1753	ln_spp_own_gauss = .FALSE.
1754	IF( nn_spp > 0 ) THEN
1755	IF( spp_filter_pass > 0 .AND. spp_filter_pass .NE. sppt_filter_pass ) ln_spp_own_gauss = .TRUE.
1756	IF( rn_spp_tau > 0._wp .AND. rn_spp_tau .NE. rn_sppt_tau ) ln_spp_own_gauss = .TRUE.
1757	ENDIF
1758
1759	IF(lwp) THEN
1760	WRITE(numout,*) ' SPPT on tracers : ',ln_sppt_tra
1761	WRITE(numout,*) ' SPPT on dynamics : ',ln_sppt_dyn
1762	WRITE(numout,*) ' SPPT on sea-ice : ',ln_sppt_ice
1763	WRITE(numout,*) ' SPP perturbed params: ',nn_spp
1764	WRITE(numout,*) ' SKEB scheme : ',ln_skeb
1765	WRITE(numout,*) ' SPP with own scales : ',ln_spp_own_gauss
1766	WRITE(numout,*) ' SKEB with own scales: ',ln_skeb_own_gauss
1767	IF(lflush) CALL flush(numout)
1768	ENDIF
1769
1770	IF( sppt_tra_alloc() /= 0) CALL ctl_stop( 'STOP', 'sppt_alloc: unable to allocate arrays' )
1771
1772	spptt = 0._wp ; sppts = 0._wp
1773	spptu = 0._wp ; spptv = 0._wp
1774
1775	! Find filter attenuation factor
1776
1777	flt_fac = sppt_fltfac( sppt_filter_pass )
1778	rdt_sppt = nn_rndm_freq * rn_rdt
1779
1780	IF( ln_sppt_taumap ) THEN
1781	CALL iom_open ( 'sppt_tau_map', inum )
1782	CALL iom_get ( inum, jpdom_data, 'tau', sppt_tau )
1783	CALL iom_close( inum )
1784	ELSE
1785	sppt_tau(:,:) = rn_sppt_tau
1786	ENDIF
1787
1788	IF( nn_deftau .EQ. 1 ) THEN
1789	! Decorrelation time-scale defined in days
1790	sppt_tau(:,:) = exp( -rdt_sppt / (sppt_tau(:,:)*86400._wp) )
1791	ELSE
1792	! Decorrelation time-scale defined in timesteps
1793	sppt_tau(:,:) = exp( -REAL( nn_rndm_freq, wp) / sppt_tau(:,:) )
1794	ENDIF
1795
1796	IF( ln_distcoast ) THEN
1797	CALL iom_open('dist.coast', inum )
1798	CALL iom_get(inum,jpdom_data,'Tcoast',zdc)
1799	CALL iom_close( inum )
1800	zdc = 1._wp - exp(-zdc/rn_distcoast)
1801	CALL lbc_lnk( zdc , 'T', 1._wp)
1802	ELSE
1803	zdc = 1._wp
1804	ENDIF
1805
1806	! Initialize
1807	sppt_a = sppt_tau
1808	sppt_b = sqrt ( 1._wp - sppt_tau*sppt_tau )
1809
1810	IF(lwp) THEN
1811	WRITE(numout,*)
1812	WRITE(numout,) ' *** SPPT SCHEME'
1813	WRITE(numout,*) ' Definit. time-scale : ',nn_deftau
1814	WRITE(numout,*) ' Decorr. time-scale : ',rn_sppt_tau
1815	WRITE(numout,) ' Mean AR1 a term : ',SUM(sppt_a)/REAL(jpijpj)
1816	WRITE(numout,) ' Mean AR1 b term : ',SUM(sppt_b)/REAL(jpijpj)
1817	WRITE(numout,*)
1818	IF(lflush) CALL flush(numout)
1819	ENDIF
1820
1821	gauss_b = 0._wp
1822	! Weigths
1823	gauss_w(:) = 1.0_wp
1824	IF( nn_vwei .eq. 1 ) THEN
1825	gauss_w(1) = 0.0_wp
1826	gauss_w(2) = 0.5_wp
1827	gauss_w(jpk) = 0.0_wp
1828	gauss_w(jpkm1)= 0.5_wp
1829	ENDIF
1830
1831	IF( ln_spp_own_gauss ) THEN
1832	IF( spp_alloc() /= 0) CALL ctl_stop( 'STOP', 'spp_alloc: unable to allocate arrays' )
1833	flt_fac_p = sppt_fltfac( spp_filter_pass )
1834	spp_tau (:, :) = rn_spp_tau
1835	IF( nn_deftau .EQ. 1 ) THEN
1836	spp_tau(:,:) = spp_tau(:,:) * 86400._wp
1837	spp_tau(:,:) = exp( -rdt_sppt / (spp_tau) )
1838	ELSE
1839	spp_tau(:,:) = exp( -1._wp / spp_tau )
1840	ENDIF
1841	spp_a = spp_tau
1842	spp_b = sqrt ( 1._wp - spp_tau*spp_tau )
1843	gauss_bp = 0._wp
1844	IF(lwp) THEN
1845	WRITE(numout,*)
1846	WRITE(numout,) ' *** SPP SCHEME'
1847	WRITE(numout,*) ' Definit. time-scale : ',nn_deftau
1848	WRITE(numout,*) ' Decorr. time-scale : ',rn_spp_tau
1849	WRITE(numout,) ' Mean AR1 a term : ',SUM(spp_a)/REAL(jpijpj)
1850	WRITE(numout,) ' Mean AR1 b term : ',SUM(spp_b)/REAL(jpijpj)
1851	WRITE(numout,*)
1852	IF(lflush) CALL flush(numout)
1853	ENDIF
1854	ENDIF
1855
1856	IF( ln_skeb_own_gauss ) THEN
1857	IF( skeb_alloc() /= 0) CALL ctl_stop( 'STOP', 'skeb_alloc: unable to allocate arrays' )
1858	flt_fac_k = sppt_fltfac( skeb_filter_pass )
1859	skeb_tau (:, :) = rn_skeb_tau
1860	IF( nn_deftau .EQ. 1 ) THEN
1861	skeb_tau(:,:) = skeb_tau(:,:) * 86400._wp
1862	skeb_tau(:,:) = exp( -rdt_sppt / (skeb_tau) )
1863	ELSE
1864	skeb_tau(:,:) = exp( -1._wp / skeb_tau )
1865	ENDIF
1866	skeb_a = skeb_tau
1867	skeb_b = sqrt ( 1._wp - skeb_tau*skeb_tau )
1868	gauss_bk = 0._wp
1869	IF(lwp) THEN
1870	WRITE(numout,*)
1871	WRITE(numout,) ' *** SEB SCHEME'
1872	WRITE(numout,*) ' Definit. time-scale : ',nn_deftau
1873	WRITE(numout,*) ' Decorr. time-scale : ',rn_skeb_tau
1874	WRITE(numout,) ' Mean AR1 a term : ',SUM(skeb_a)/REAL(jpijpj)
1875	WRITE(numout,) ' Mean AR1 b term : ',SUM(skeb_b)/REAL(jpijpj)
1876	WRITE(numout,*)
1877	IF(lflush) CALL flush(numout)
1878	ENDIF
1879	ENDIF
1880
1881	IF( ln_skeb_own_gauss .OR. ln_spp_own_gauss ) &
1882	& ALLOCATE ( g2d_save(jpi,jpj) )
1883
1884	CALL stopack_rst( nit000, 'READ' )
1885
1886	IF(lwp .and. ln_stopack_diags) &
1887	CALL ctl_opn(numdiag, 'stopack.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
1888
1889	END SUBROUTINE stopack_init
1890	!
1891	SUBROUTINE stopack_rst( kt, cdrw )
1892	!!----------------------------------------------------------------------
1893	!! * ROUTINE stopack_rst *
1894	!!
1895	!! ** Purpose : Read/write from/to restarsts seed and perturbation field
1896	!!----------------------------------------------------------------------
1897	INTEGER , INTENT(in) :: kt ! ocean time-step
1898	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1899	INTEGER :: id1, jseed
1900	CHARACTER(LEN=10) :: clseed='spsd0_0000'
1901	INTEGER(KIND=8) :: ziseed(4) ! RNG seeds in integer type
1902	INTEGER(KIND=8) :: ivals(8)
1903	REAL(wp) :: zrseed4(4) ! RNG seeds in integer type
1904	REAL(wp) :: zrseed2d(jpi,jpj)
1905	!
1906	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
1907	! ! ---------------
1908	IF(lwp) WRITE(numout,) ' stopack-rst *'
1909	IF( ln_rstart ) THEN !* Read the restart file
1910	id1 = iom_varid( numror, 'sppt_b', ldstop = .FALSE. )
1911	!
1912	IF( id1 > 0 ) THEN
1913	CALL iom_get( numror, jpdom_autoglo, 'sppt_b', gauss_b )
1914	DO jseed = 1 , 4
1915	WRITE(clseed(5:5) ,'(i1.1)') jseed
1916	CALL iom_get( numror, jpdom_autoglo, clseed(1:5) , zrseed2d )
1917	zrseed4(jseed) = zrseed2d(2,2)
1918	ziseed(jseed) = TRANSFER( zrseed4(jseed) , ziseed(jseed) )
1919	END DO
1920	IF (lwp) THEN
1921	WRITE(numout,*) 'Read ziseed4 = ',zrseed4(:)
1922	WRITE(numout,*) 'Read ziseed = ',ziseed(:)
1923	ENDIF
1924	CALL kiss_seed( ziseed(1) , ziseed(2) , ziseed(3) , ziseed(4) )
1925	ELSE
1926	IF( lwp ) WRITE(numout,*) ' ===>>>> : previous run without sppt_b'
1927	IF( ln_stopack_restart ) CALL ctl_stop( ' ln_stopack_restart TRUE :',&
1928	& ' variable not found in restart file ')
1929	IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of sppt_b to 0.'
1930	gauss_b = 0._wp
1931	IF ( .NOT. ln_stopack_repr ) THEN
1932	CALL date_and_time(VALUES=ivals)
1933	DO jseed=1,4
1934	nn_stopack_seed(jseed) = nn_stopack_seed(jseed) + INT(ivals(4+jseed))
1935	ENDDO
1936	ENDIF
1937	DO jseed=1,4
1938	ziseed(jseed) = TRANSFER(nn_stopack_seed(jseed)+narea, ziseed(jseed))
1939	ENDDO
1940	IF(lwp) WRITE(numout,*) ' ===>>>> Seed, nn_stopack_seed+narea',nn_stopack_seed+narea
1941	IF(lwp) WRITE(numout,*) ' ===>>>> Seed, ziseed',ziseed
1942	CALL kiss_seed( ziseed(1) , ziseed(2) , ziseed(3) , ziseed(4) )
1943	ENDIF
1944	IF( ln_spp_own_gauss ) THEN
1945	id1 = iom_varid( numror, 'spp_b', ldstop = .FALSE. )
1946	IF( id1 > 0 ) THEN
1947	CALL iom_get( numror, jpdom_autoglo, 'spp_b', gauss_bp )
1948	ELSE
1949	IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of spp_b to 0.'
1950	gauss_bp = 0._wp
1951	ENDIF
1952	ENDIF
1953	IF( ln_skeb_own_gauss ) THEN
1954	id1 = iom_varid( numror, 'skeb_b', ldstop = .FALSE. )
1955	IF( id1 > 0 ) THEN
1956	CALL iom_get( numror, jpdom_autoglo, 'skeb_b', gauss_bk )
1957	ELSE
1958	IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of skeb_b to 0.'
1959	gauss_bk = 0._wp
1960	ENDIF
1961	ENDIF
1962	ELSE
1963	gauss_b = 0._wp
1964	gauss_bp = 0._wp
1965	gauss_bk = 0._wp
1966	IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of STOPACK to 0.'
1967	IF ( .NOT. ln_stopack_repr ) THEN
1968	CALL date_and_time(VALUES=ivals)
1969	DO jseed=1,4
1970	nn_stopack_seed(jseed) = nn_stopack_seed(jseed) + INT(ivals(4+jseed))
1971	ENDDO
1972	ENDIF
1973	DO jseed=1,4
1974	ziseed(jseed) = TRANSFER(nn_stopack_seed(jseed)+narea, ziseed(jseed))
1975	ENDDO
1976	IF(lwp) WRITE(numout,*) ' ===>>>> Seed, nn_stopack_seed+narea',nn_stopack_seed+narea
1977	IF(lwp) WRITE(numout,*) ' ===>>>> Seed, ziseed',ziseed
1978	CALL kiss_seed( ziseed(1) , ziseed(2) , ziseed(3) , ziseed(4) )
1979	ENDIF
1980	!
1981	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
1982	! ! -------------------
1983	IF(lwp) WRITE(numout,*) '---- stopack-rst ----'
1984	CALL iom_rstput( kt, nitrst, numrow, 'sppt_b', gauss_b )
1985	IF(ln_spp_own_gauss) CALL iom_rstput( kt, nitrst, numrow, 'spp_b', gauss_bp )
1986	IF(ln_skeb_own_gauss) CALL iom_rstput( kt, nitrst, numrow, 'skeb_b', gauss_bk )
1987	CALL kiss_state( ziseed(1) , ziseed(2) , ziseed(3) , ziseed(4) )
1988	DO jseed=1,4
1989	zrseed4(jseed) = TRANSFER(ziseed(jseed), zrseed4(jseed))
1990	ENDDO
1991	IF (lwp) THEN
1992	WRITE(numout,*) 'Write ziseed4 = ',zrseed4(:)
1993	WRITE(numout,*) 'Write ziseed = ',ziseed(:)
1994	ENDIF
1995	DO jseed = 1 , 4
1996	WRITE(clseed(5:5) ,'(i1.1)') jseed
1997	zrseed2d(:,:) = zrseed4(jseed)
1998	CALL iom_rstput( kt, nitrst, numrow, clseed(1:5) , zrseed2d )
1999	END DO
2000	!
2001	ENDIF
2002	!
2003	END SUBROUTINE stopack_rst
2004	!
2005	INTEGER FUNCTION sppt_tra_alloc()
2006	!!---------------------------------------------------------------------
2007	!! * FUNCTION sppt_tra_alloc *
2008	!!---------------------------------------------------------------------
2009	!
2010	ALLOCATE( spptt(jpi,jpj,jpk) , sppts(jpi,jpj,jpk) , gauss_n(jpi,jpj,jpk) ,&
2011	gauss_nu(jpi,jpj,jpk) , gauss_nv(jpi,jpj,jpk) , &
2012	spptu(jpi,jpj,jpk) , spptv(jpi,jpj,jpk) , gauss_n_2d(jpi,jpj) ,&
2013	gauss_b (jpi,jpj), sppt_tau(jpi,jpj), sppt_a(jpi,jpj), sppt_b(jpi,jpj), gauss_w(jpk),&
2014	zdc(jpi,jpj,jpk), STAT= sppt_tra_alloc )
2015	!
2016	IF( lk_mpp ) CALL mpp_sum ( sppt_tra_alloc )
2017	IF( sppt_tra_alloc /= 0) CALL ctl_warn('sppt_tra_alloc: failed to allocate arrays')
2018
2019	IF(nn_spp_bfr>0) THEN
2020	ALLOCATE(coeff_bfr(jpi,jpj), STAT= sppt_tra_alloc )
2021	IF( lk_mpp ) CALL mpp_sum ( sppt_tra_alloc )
2022	IF( sppt_tra_alloc /= 0) CALL ctl_warn('sppt_tra_alloc: failed to allocate arrays')
2023	coeff_bfr = 0._wp
2024	ENDIF
2025	gauss_b = 0._wp
2026	gauss_n = 0._wp
2027	gauss_nu= 0._wp
2028	gauss_nv= 0._wp
2029	gauss_n_2d = 0._wp
2030
2031	END FUNCTION sppt_tra_alloc
2032
2033	INTEGER FUNCTION skeb_alloc()
2034	!!---------------------------------------------------------------------
2035	!! * FUNCTION skeb_alloc *
2036	!!---------------------------------------------------------------------
2037	!
2038	ALLOCATE( gauss_n_2d_k(jpi,jpj) , gauss_bk (jpi,jpj),skeb_tau(jpi,jpj),&
2039	skeb_a(jpi,jpj), skeb_b(jpi,jpj), STAT= skeb_alloc )
2040	!
2041	IF( lk_mpp ) CALL mpp_sum ( skeb_alloc )
2042	IF( skeb_alloc /= 0) CALL ctl_warn('sppt_tra_alloc: failed to allocate arrays')
2043	gauss_n_2d_k = 0._wp
2044	gauss_bk = 0._wp
2045
2046	END FUNCTION skeb_alloc
2047
2048	INTEGER FUNCTION spp_alloc()
2049	!!---------------------------------------------------------------------
2050	!! * FUNCTION spp_alloc *
2051	!!---------------------------------------------------------------------
2052	!
2053	ALLOCATE( gauss_n_2d_p(jpi,jpj), gauss_bp (jpi,jpj),spp_tau(jpi,jpj), &
2054	spp_a(jpi,jpj), spp_b(jpi,jpj), STAT= spp_alloc )
2055	!
2056	IF( lk_mpp ) CALL mpp_sum ( spp_alloc )
2057	IF( spp_alloc /= 0) CALL ctl_warn('sppt_tra_alloc: failed to allocate arrays')
2058	gauss_n_2d_p = 0._wp
2059	gauss_bp = 0._wp
2060
2061	END FUNCTION spp_alloc
2062
2063	SUBROUTINE sppt_flt( psto )
2064	!!----------------------------------------------------------------------
2065	!! * ROUTINE sppt_flt *
2066	!!
2067	!! ** Purpose : apply horizontal Laplacian filter to input array
2068	!! Adapted from STO package
2069	!!----------------------------------------------------------------------
2070	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: psto
2071	REAL(wp), DIMENSION(jpi,jpj) :: pstm
2072	!!
2073	INTEGER :: ji, jj
2074
2075	pstm = psto
2076	DO jj = 2, jpj-1
2077	DO ji = 2, jpi-1
2078	psto(ji,jj) = 0.5_wp * pstm(ji,jj) + 0.125_wp * &
2079	& ( pstm(ji-1,jj) + pstm(ji+1,jj) + &
2080	& pstm(ji,jj-1) + pstm(ji,jj+1) )
2081	END DO
2082	END DO
2083	CALL lbc_lnk( psto , 'T', 1._wp )
2084
2085	END SUBROUTINE sppt_flt
2086
2087	SUBROUTINE sppt_rand2d( psto )
2088	!!----------------------------------------------------------------------
2089	!! * ROUTINE sppt_rand2d *
2090	!!
2091	!! ** Purpose : fill input array with white Gaussian noise
2092	!! Adapted from STO package
2093	!!----------------------------------------------------------------------
2094	REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: psto
2095	!!
2096	INTEGER :: ji, jj
2097	REAL(KIND=8) :: gran ! Gaussian random number (forced KIND=8 as in kiss_gaussian)
2098
2099	DO jj = 1, jpj
2100	DO ji = 1, jpi
2101	CALL kiss_gaussian( gran )
2102	psto(ji,jj) = gran*rn_sppt_stdev
2103	END DO
2104	END DO
2105
2106	END SUBROUTINE sppt_rand2d
2107
2108	FUNCTION sppt_fltfac( kpasses )
2109	!!----------------------------------------------------------------------
2110	!! * FUNCTION sppt_fltfac *
2111	!!
2112	!! ** Purpose : compute factor to restore standard deviation
2113	!! as a function of the number of passes
2114	!! of the Laplacian filter
2115	!! From STO package
2116	!!----------------------------------------------------------------------
2117	INTEGER, INTENT(in) :: kpasses
2118	REAL(wp) :: sppt_fltfac
2119	!!
2120	INTEGER :: jpasses, ji, jj, jflti, jfltj
2121	INTEGER, DIMENSION(-1:1,-1:1) :: pflt0
2122	! 16-b reals to avoid overflows
2123	INTEGER, PARAMETER :: qp = SELECTED_REAL_KIND(33,4931)
2124	REAL(qp), DIMENSION(:,:), ALLOCATABLE :: pfltb
2125	REAL(qp), DIMENSION(:,:), ALLOCATABLE :: pflta
2126	REAL(qp) :: ratio
2127	REAL(qp) :: aux0, aux1
2128	pflt0(-1,-1) = 0 ; pflt0(-1,0) = 1 ; pflt0(-1,1) = 0
2129	pflt0( 0,-1) = 1 ; pflt0( 0,0) = 4 ; pflt0( 0,1) = 1
2130	pflt0( 1,-1) = 0 ; pflt0( 1,0) = 1 ; pflt0( 1,1) = 0
2131	ALLOCATE(pfltb(-kpasses-1:kpasses+1,-kpasses-1:kpasses+1))
2132	ALLOCATE(pflta(-kpasses-1:kpasses+1,-kpasses-1:kpasses+1))
2133	pfltb(:,:) = 0
2134	pfltb(0,0) = 1
2135	DO jpasses = 1, kpasses
2136	pflta(:,:) = 0._qp
2137	DO jflti= -1, 1
2138	DO jfltj= -1, 1
2139	DO ji= -kpasses, kpasses
2140	DO jj= -kpasses, kpasses
2141	pflta(ji,jj) = pflta(ji,jj) + pfltb(ji+jflti,jj+jfltj) * pflt0(jflti,jfltj)
2142	ENDDO
2143	ENDDO
2144	ENDDO
2145	ENDDO
2146	pfltb(:,:) = pflta(:,:)
2147	ENDDO
2148	ratio = SUM(pfltb(:,:))
2149	aux0 = SUM(pfltb(:,:)*pfltb(:,:))
2150	aux1 = ratio*ratio
2151	ratio = aux1 / aux0
2152	ratio = SQRT(ratio)
2153	DEALLOCATE(pfltb,pflta)
2154	sppt_fltfac = REAL(ratio, wp)
2155	END FUNCTION sppt_fltfac
2156
2157	SUBROUTINE skeb_comp( lt )
2158	!!----------------------------------------------------------------------
2159	!! * ROUTINE skeb_comp *
2160	!!
2161	!! ** Purpose : Computation of energy dissipation terms
2162	!! This is a wrapper to the enrgy terms computation
2163	!! routines
2164	!!----------------------------------------------------------------------
2165	INTEGER, INTENT(IN) :: lt
2166	CALL skeb_dnum ( lt )
2167	CALL skeb_dcon ( lt )
2168	END SUBROUTINE skeb_comp
2169
2170	SUBROUTINE bsfcomp( kt )
2171	!!----------------------------------------------------------------------
2172	!! * ROUTINE bsfcomp *
2173	!!
2174	!! ** Purpose : Online computation of barotropic streamfunction
2175	!! For diagnostics purposes
2176	!! This routine is not explicitly called anywhere
2177	!! and utilizes low-level MPI routines
2178	!!----------------------------------------------------------------------
2179	USE MPI
2180	!
2181	INTEGER, INTENT(IN) :: kt
2182	!
2183	REAL(wp) :: dtrpv(jpi,jpj)
2184	INTEGER :: jk,ji,jj
2185	!
2186	#if defined key_mpp_mpi
2187	INTEGER, DIMENSION(1) :: ish
2188	INTEGER, DIMENSION(2) :: ish2
2189	INTEGER :: ijp,tag,ierr,jp,isend,irecv
2190	REAL(wp), DIMENSION(jpj) :: zwrk ! mask flux array at V-point
2191	REAL(wp), DIMENSION(jpi) :: zwr2 ! mask flux array at V-point
2192	INTEGER :: istatus(MPI_STATUS_SIZE)
2193	#endif
2194	IF(kt .EQ. nit000 ) THEN
2195	#ifdef NEMO_V34
2196	ln_dpsiu = .FALSE.
2197	ln_dpsiv = .FALSE.
2198	#else
2199	IF (iom_use('bsft') ) ln_dpsiu = .TRUE.
2200	IF (iom_use('bsftv') ) ln_dpsiv = .TRUE.
2201	#endif
2202	IF( ln_dpsiv ) ALLOCATE ( dpsiv(jpi,jpj) )
2203	IF( ln_dpsiu ) ALLOCATE ( dpsiu(jpi,jpj) )
2204	IF( ln_dpsiv .OR. ln_dpsiu ) THEN
2205	jpri = nint ( REAL(nimpp) / REAL(jpi) ) + 1
2206	jprj = nint ( REAL(njmpp) / REAL(jpj) ) + 1
2207	ENDIF
2208	ENDIF
2209
2210	IF( ln_dpsiv ) THEN
2211	dtrpv = 0._wp
2212	DO jk = 1,jpkm1
2213	dtrpv = dtrpv + fse3v(:,:,jk) * e1v(:,:) * vn(:,:,jk)
2214	ENDDO
2215	dpsiv(1,:)=0._wp
2216	DO ji = 2,jpi
2217	dpsiv(ji,:) = dpsiv(ji-1,:) + dtrpv(ji,:)
2218	END DO
2219	ENDIF
2220
2221	IF( ln_dpsiu ) THEN
2222	dtrpv = 0._wp
2223	DO jk = 1,jpkm1
2224	dtrpv = dtrpv + fse3u(:,:,jk) * e2u(:,:) * un(:,:,jk)
2225	ENDDO
2226	dpsiu(:,1)=0._wp
2227	DO jj = 2,jpj
2228	dpsiu(:,jj) = dpsiu(:,jj-1) + dtrpv(:,jj)
2229	END DO
2230	ENDIF
2231
2232	#if defined key_mpp_mpi
2233	IF ( ln_dpsiv ) THEN
2234	DO jp=1,jpni-1
2235	IF( jpri == jp ) THEN ! SEND TO EAST
2236	zwrk(1:jpj) = dpsiv(jpi-1,:)
2237	tag=2000+narea
2238	CALL mpi_isend(zwrk, jpj, mpi_double_precision, noea, tag, mpi_comm_opa, isend,ierr)
2239	ELSEIF ( jpri == jp+1 ) THEN ! RECEIVE FROM WEST
2240	CALL mpi_irecv(zwrk, jpj, mpi_double_precision, nowe, mpi_any_tag, mpi_comm_opa, irecv,ierr)
2241	ENDIF
2242	IF(jpri == jp) CALL mpi_wait(isend, istatus, ierr)
2243	IF(jpri == jp+1 ) THEN
2244	CALL mpi_wait(irecv, istatus, ierr)
2245	DO ji=1,jpi
2246	dpsiv(ji,:) = dpsiv(ji,:) + zwrk(:)
2247	ENDDO
2248	ENDIF
2249	ENDDO
2250	ENDIF
2251
2252	IF ( ln_dpsiv ) THEN
2253	DO jp=1,jpnj-1
2254	IF( jprj == jp ) THEN ! SEND TO NORTH
2255	zwr2(1:jpi) = dpsiu(:,jpj-1)
2256	tag=3000+narea
2257	CALL mpi_isend(zwr2, jpi, mpi_double_precision, nono, tag, mpi_comm_opa, isend,ierr)
2258	ELSEIF ( jprj == jp+1 ) THEN ! RECEIVE FROM SOUTH
2259	CALL mpi_irecv(zwr2, jpi, mpi_double_precision, noso, mpi_any_tag, mpi_comm_opa, irecv,ierr)
2260	ENDIF
2261	IF(jprj == jp) CALL mpi_wait(isend, istatus, ierr)
2262	IF(jprj == jp+1 ) THEN
2263	CALL mpi_wait(irecv, istatus, ierr)
2264	DO ji=1,jpj
2265	dpsiu(:,ji) = dpsiu(:,ji) + zwr2(:)
2266	ENDDO
2267	ENDIF
2268	ENDDO
2269	ENDIF
2270	#endif
2271	IF (ln_dpsiu ) THEN
2272	CALL lbc_lnk(dpsiu,'T',1._wp)
2273	ENDIF
2274	IF (ln_dpsiv ) THEN
2275	CALL lbc_lnk(dpsiv,'T',1._wp)
2276	ENDIF
2277
2278	IF(kt == nitend ) THEN
2279	IF (ln_dpsiv ) DEALLOCATE ( dpsiv )
2280	IF (ln_dpsiu ) DEALLOCATE ( dpsiu )
2281	ENDIF
2282
2283	END SUBROUTINE
2284
2285	SUBROUTINE skeb_dnum ( mt )
2286	!!----------------------------------------------------------------------
2287	!! * ROUTINE skeb_dnum *
2288	!!
2289	!! ** Purpose : Computation of numerical energy dissipation
2290	!! For later use in the SKEB scheme
2291	!!----------------------------------------------------------------------
2292	INTEGER, INTENT(IN) :: mt
2293	REAL :: ds,dt,dtot,kh2
2294	INTEGER :: ji,jj,jk
2295
2296	IF ( mt .eq. nit000 ) THEN
2297	ALLOCATE ( dnum(jpi,jpj,jpk) )
2298	dnum (:,:,: ) = 0._wp
2299	ENDIF
2300
2301	kh2 = rn_kh * rn_kh
2302
2303	IF( mt .eq. nit000 .OR. MOD( mt - 1, nn_dcom_freq ) == 0 ) THEN
2304
2305	DO jk=1,jpkm1
2306	DO jj=2,jpj
2307	DO ji=2,jpi
2308	! Shear
2309	ds = (vn(ji,jj,jk)-vn(ji-1,jj,jk))vmask(ji,jj,jk)vmask(ji-1,jj,jk)/ e1v(ji,jj) + &
2310	(un(ji,jj,jk)-un(ji,jj-1,jk))umask(ji,jj,jk)umask(ji,jj-1,jk)/ e2u(ji,jj)
2311	! Tension
2312	dt = (vn(ji,jj,jk)-vn(ji-1,jj,jk))vmask(ji,jj,jk)vmask(ji-1,jj,jk)/ e2v(ji,jj) + &
2313	(un(ji,jj,jk)-un(ji,jj-1,jk))umask(ji,jj,jk)umask(ji,jj-1,jk)/ e1u(ji,jj)
2314
2315	dtot = sqrt( dsds + dtdt ) * tmask(ji,jj,jk)
2316	dnum(ji,jj,jk) = dtotdtotdtotkh2e1t(ji,jj)*e2t(ji,jj)
2317	ENDDO
2318	ENDDO
2319	ENDDO
2320
2321	CALL lbc_lnk(dnum,'T',1._wp)
2322
2323	ENDIF
2324
2325	END SUBROUTINE
2326
2327	SUBROUTINE skeb_dcon ( mt )
2328	!!----------------------------------------------------------------------
2329	!! * ROUTINE skeb_dcon *
2330	!!
2331	!! ** Purpose : Computation of convective energy dissipation
2332	!! For later use in the SKEB scheme
2333	!! Two formulation are implemented, based on buoyancy or
2334	!! convective mass flux formulations, respectively
2335	!!----------------------------------------------------------------------
2336	INTEGER, INTENT(IN) :: mt
2337	REAL(wp) :: zz, mf1,mf2,kc2
2338	INTEGER :: ji,jj,jk
2339
2340	IF ( mt .eq. nit000 ) THEN
2341	ALLOCATE ( dcon(jpi,jpj,jpk) )
2342	dcon (:,:,: ) = 0._wp
2343	ENDIF
2344
2345	kc2 = rn_kc * rn_kc
2346
2347	IF( mt .eq. nit000 .OR. MOD( mt - 1, nn_dcom_freq ) == 0 ) THEN
2348
2349	IF(nn_dconv .eq. 1 ) THEN
2350
2351	DO jk=2,jpkm1
2352	DO jj=1,jpj
2353	DO ji=1,jpi
2354
2355	zz = - gravavt(ji,jj,jk) ( rhd(ji,jj,jk)-rhd(ji,jj,jk-1) ) * wmask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) &
2356	& / ( rau0 * fse3w(ji,jj,jk) )
2357
2358	dcon(ji,jj,jk) = kc2zze1t(ji,jj)e2t(ji,jj)rau0 / fse3w(ji,jj,jk)
2359
2360	ENDDO
2361	ENDDO
2362	ENDDO
2363
2364	ELSEIF (nn_dconv .eq. 2 ) THEN
2365
2366	DO jk=2,jpkm1
2367	DO jj=1,jpj
2368	DO ji=1,jpi
2369
2370	mf1 = wn(ji,jj,jk+1)e1t(ji,jj)e2t(ji,jj)
2371	mf2 = wn(ji,jj,jk-1)e1t(ji,jj)e2t(ji,jj)
2372	dcon(ji,jj,jk) = rn_kc * (mf1-mf2)(mf1-mf2) tmask(ji,jj,jk+1) / (fse3w(ji,jj,jk)rau0rau0)
2373
2374	ENDDO
2375	ENDDO
2376	ENDDO
2377
2378	ENDIF
2379
2380	CALL lbc_lnk(dcon,'T',1._wp)
2381
2382	ENDIF
2383
2384	END SUBROUTINE
2385
2386	SUBROUTINE skeb_apply ( mt)
2387	!!----------------------------------------------------------------------
2388	!! * ROUTINE skeb_apply *
2389	!!
2390	!! ** Purpose : Application of SKEB perturbation
2391	!! Convective and Numerical energy dissipation are
2392	!! multiplied by the beta terms
2393	!!----------------------------------------------------------------------
2394
2395	INTEGER, INTENT(IN) :: mt
2396	INTEGER :: ji,jj,jk
2397	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ustar,vstar,psi
2398	REAL(wp) :: mi,ma
2399
2400	IF( ln_stopack_diags ) THEN
2401
2402	psi = 0._wp
2403	IF(ln_skeb_own_gauss) THEN
2404	DO jk=1,jpkm1
2405	psi(:,:,jk) = rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk) ) * gauss_n_2d_k(:,:)
2406	ENDDO
2407	ELSE
2408	DO jk=1,jpkm1
2409	psi(:,:,jk) = rn_skeb_stdev * rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk) ) &
2410	& * gauss_n_2d(:,:) / rn_sppt_stdev
2411	ENDDO
2412	ENDIF
2413	ustar = 0._wp
2414	vstar = 0._wp
2415	DO jk=1,jpkm1
2416	DO jj=2,jpj
2417	DO ji=2,jpi
2418	ustar(ji,jj,jk) = ( psi(ji,jj,jk)-psi(ji,jj-1,jk) )umask(ji,jj,jk)umask(ji,jj-1,jk)/ e2u(ji,jj)
2419	vstar(ji,jj,jk) = - ( psi(ji,jj,jk)-psi(ji-1,jj,jk) )vmask(ji,jj,jk)vmask(ji-1,jj,jk)/ e1v(ji,jj)
2420	ENDDO
2421	ENDDO
2422	ENDDO
2423	mi = MAXVAL(ABS(ustar))
2424	ma = MAXVAL(ABS(vstar))
2425	IF(lk_mpp) CALL mpp_max(mi)
2426	IF(lk_mpp) CALL mpp_max(ma)
2427	IF(lwp) WRITE(numdiag,9304) lkt,'DNUM ABS CURRENT' ,mi,ma
2428
2429	rn_mmax ( jk_skeb_dnum ) = MAX ( MAX( mi, ma), rn_mmax ( jk_skeb_dnum ) )
2430
2431	psi = 0._wp
2432	IF(ln_skeb_own_gauss) THEN
2433	DO jk=1,jpkm1
2434	psi(:,:,jk) = rn_skeb * sqrt( rn_beta_con * dcon(:,:,jk) ) * gauss_n_2d_k(:,:)
2435	ENDDO
2436	ELSE
2437	DO jk=1,jpkm1
2438	psi(:,:,jk) = rn_skeb_stdev * rn_skeb * sqrt( rn_beta_con * dcon(:,:,jk) ) &
2439	& * gauss_n_2d(:,:) / rn_sppt_stdev
2440	ENDDO
2441	ENDIF
2442	ustar = 0._wp
2443	vstar = 0._wp
2444	DO jk=1,jpkm1
2445	DO jj=2,jpj
2446	DO ji=2,jpi
2447	ustar(ji,jj,jk) = ( psi(ji,jj,jk)-psi(ji,jj-1,jk) )umask(ji,jj,jk)umask(ji,jj-1,jk)/ e2u(ji,jj)
2448	vstar(ji,jj,jk) = - ( psi(ji,jj,jk)-psi(ji-1,jj,jk) )vmask(ji,jj,jk)vmask(ji-1,jj,jk)/ e1v(ji,jj)
2449	ENDDO
2450	ENDDO
2451	ENDDO
2452	mi = MAXVAL(ABS(ustar))
2453	ma = MAXVAL(ABS(vstar))
2454	IF(lk_mpp) CALL mpp_max(mi)
2455	IF(lk_mpp) CALL mpp_max(ma)
2456	IF(lwp) WRITE(numdiag,9304) lkt,'DCON ABS CURRENT' ,mi,ma
2457
2458	rn_mmax ( jk_skeb_dcon ) = MAX ( MAX( mi, ma), rn_mmax ( jk_skeb_dcon ) )
2459
2460	9304 FORMAT(' it :', i8, ' ', A16, ' Min: ',d10.3 , ' Max: ',d10.3)
2461
2462	ENDIF
2463
2464	psi = 0._wp
2465	IF(ln_skeb_own_gauss) THEN
2466	DO jk=1,jpkm1
2467	psi(:,:,jk) = rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk)+ rn_beta_con * dcon(:,:,jk) ) * gauss_n_2d_k(:,:)
2468	ENDDO
2469	ELSE
2470	DO jk=1,jpkm1
2471	psi(:,:,jk) = rn_skeb_stdev * rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk)+ rn_beta_con * dcon(:,:,jk) ) &
2472	& * gauss_n_2d(:,:) / rn_sppt_stdev
2473	ENDDO
2474	ENDIF
2475
2476	ustar = 0._wp
2477	vstar = 0._wp
2478
2479	DO jk=1,jpkm1
2480	DO jj=2,jpj
2481	DO ji=2,jpi
2482	ustar(ji,jj,jk) = ( psi(ji,jj,jk)-psi(ji,jj-1,jk) )umask(ji,jj,jk)umask(ji,jj-1,jk)/ e2u(ji,jj)
2483	vstar(ji,jj,jk) = - ( psi(ji,jj,jk)-psi(ji-1,jj,jk) )vmask(ji,jj,jk)vmask(ji-1,jj,jk)/ e1v(ji,jj)
2484	ENDDO
2485	ENDDO
2486	ENDDO
2487
2488	IF( ln_stopack_diags ) THEN
2489
2490	mi = MAXVAL(ABS(dnum))
2491	ma = MAXVAL(ABS(dcon))
2492	IF(lk_mpp) CALL mpp_max(mi)
2493	IF(lk_mpp) CALL mpp_max(ma)
2494	IF(lwp) WRITE(numdiag,9302) lkt,'DNUM AND DCON MX' ,mi,ma
2495
2496	mi = MINVAL(ustar)
2497	ma = MAXVAL(ustar)
2498	IF(lk_mpp) CALL mpp_min(mi)
2499	IF(lk_mpp) CALL mpp_max(ma)
2500	IF(lwp) WRITE(numdiag,9302) lkt,'SKEB U-CURRENT ' ,mi,ma
2501
2502	rn_mmax ( jk_skeb_tot ) = MAX ( MAX( abs(mi), ma), rn_mmax ( jk_skeb_tot ) )
2503
2504	mi = MINVAL(vstar)
2505	ma = MAXVAL(vstar)
2506	IF(lk_mpp) CALL mpp_min(mi)
2507	IF(lk_mpp) CALL mpp_max(ma)
2508	IF(lwp) WRITE(numdiag,9302) lkt,'SKEB V-CURRENT ' ,mi,ma
2509	9302 FORMAT(' it :', i8, ' ', A16, ' Min: ',d10.3 , ' Max: ',d10.3)
2510
2511	rn_mmax ( jk_skeb_tot ) = MAX ( MAX( abs(mi), ma), rn_mmax ( jk_skeb_tot ) )
2512
2513	ENDIF
2514
2515	CALL lbc_lnk(ustar,'U',-1._wp)
2516	CALL lbc_lnk(vstar,'V',-1._wp)
2517
2518	DO jk=1,jpkm1
2519	DO jj=1,jpj
2520	DO ji=1,jpi
2521	un(ji,jj,jk) = un(ji,jj,jk) + ustar(ji,jj,jk)
2522	vn(ji,jj,jk) = vn(ji,jj,jk) + vstar(ji,jj,jk)
2523	ENDDO
2524	ENDDO
2525	ENDDO
2526
2527	#ifdef key_iomput
2528	IF (iom_use('ustar_skeb') ) CALL iom_put( "ustar_skeb" , ustar)
2529	IF (iom_use('vstar_skeb') ) CALL iom_put( "vstar_skeb" , vstar)
2530	#endif
2531
2532	IF ( mt .eq. nitend ) THEN
2533	IF(ALLOCATED(dnum)) DEALLOCATE ( dnum )
2534	IF(ALLOCATED(dcon)) DEALLOCATE ( dcon )
2535	IF (ln_stopack_diags .AND. lwp) CALL stopack_report
2536	ENDIF
2537
2538	END SUBROUTINE
2539
2540	!!======================================================================
2541	!! Random number generator, used in NEMO stochastic parameterization
2542	!!
2543	!!=====================================================================
2544	!! History : 3.3 ! 2011-10 (J.-M. Brankart) Original code
2545	!!----------------------------------------------------------------------
2546
2547	!!----------------------------------------------------------------------
2548	!! The module is based on (and includes) the
2549	!! 64-bit KISS (Keep It Simple Stupid) random number generator
2550	!! distributed by George Marsaglia :
2551	!! http://groups.google.com/group/comp.lang.fortran/
2552	!! browse_thread/thread/a85bf5f2a97f5a55
2553	!!----------------------------------------------------------------------
2554
2555	!!----------------------------------------------------------------------
2556	!! kiss : 64-bit KISS random number generator (period ~ 2^250)
2557	!! kiss_seed : Define seeds for KISS random number generator
2558	!! kiss_state : Get current state of KISS random number generator
2559	!! kiss_save : Save current state of KISS (for future restart)
2560	!! kiss_load : Load the saved state of KISS
2561	!! kiss_reset : Reset the default seeds
2562	!! kiss_check : Check the KISS pseudo-random sequence
2563	!! kiss_uniform : Real random numbers with uniform distribution in [0,1]
2564	!! kiss_gaussian : Real random numbers with Gaussian distribution N(0,1)
2565	!! kiss_gamma : Real random numbers with Gamma distribution Gamma(k,1)
2566	!! kiss_sample : Select a random sample from a set of integers
2567	!!
2568	!! ---CURRENTLY NOT USED IN NEMO :
2569	!! kiss_save, kiss_load, kiss_check, kiss_gamma, kiss_sample
2570	!!----------------------------------------------------------------------
2571	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
2572	!! $Id: dynhpg.F90 2528 2010-12-27 17:33:53Z rblod $
2573	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
2574	!!----------------------------------------------------------------------
2575	FUNCTION kiss()
2576	!! --------------------------------------------------------------------
2577	!! * FUNCTION kiss *
2578	!!
2579	!! ** Purpose : 64-bit KISS random number generator
2580	!!
2581	!! ** Method : combine several random number generators:
2582	!! (1) Xorshift (XSH), period 2^64-1,
2583	!! (2) Multiply-with-carry (MWC), period (2^121+2^63-1)
2584	!! (3) Congruential generator (CNG), period 2^64.
2585	!!
2586	!! overall period:
2587	!! (2^250+2^192+2^64-2^186-2^129)/6
2588	!! ~= 2^(247.42) or 10^(74.48)
2589	!!
2590	!! set your own seeds with 'kiss_seed'
2591	! --------------------------------------------------------------------
2592	IMPLICIT NONE
2593	INTEGER(KIND=i8) :: kiss, t
2594
2595	t = ISHFT(x,58) + w
2596	IF (s(x).eq.s(t)) THEN
2597	w = ISHFT(x,-6) + s(x)
2598	ELSE
2599	w = ISHFT(x,-6) + 1 - s(x+t)
2600	ENDIF
2601	x = t + x
2602	y = m( m( m(y,13_8), -17_8 ), 43_8 )
2603	z = 6906969069_8 * z + 1234567_8
2604
2605	kiss = x + y + z
2606
2607	CONTAINS
2608
2609	FUNCTION s(k)
2610	INTEGER(KIND=i8) :: s, k
2611	s = ISHFT(k,-63)
2612	END FUNCTION s
2613
2614	FUNCTION m(k, n)
2615	INTEGER(KIND=i8) :: m, k, n
2616	m = IEOR(k, ISHFT(k, n) )
2617	END FUNCTION m
2618
2619	END FUNCTION kiss
2620
2621
2622	SUBROUTINE kiss_seed(ix, iy, iz, iw)
2623	!! --------------------------------------------------------------------
2624	!! * ROUTINE kiss_seed *
2625	!!
2626	!! ** Purpose : Define seeds for KISS random number generator
2627	!!
2628	!! --------------------------------------------------------------------
2629	IMPLICIT NONE
2630	INTEGER(KIND=i8) :: ix, iy, iz, iw
2631
2632	IF(lwp) WRITE(numout,) ' ** kiss_seed called with args:'
2633	IF(lwp) WRITE(numout,) ' ** ix = ',ix
2634	IF(lwp) WRITE(numout,) ' ** iy = ',iy
2635	IF(lwp) WRITE(numout,) ' ** iz = ',iz
2636	IF(lwp) WRITE(numout,) ' ** iw = ',iw
2637
2638	x = ix
2639	y = iy
2640	z = iz
2641	w = iw
2642
2643	END SUBROUTINE kiss_seed
2644
2645
2646	SUBROUTINE kiss_state(ix, iy, iz, iw)
2647	!! --------------------------------------------------------------------
2648	!! * ROUTINE kiss_state *
2649	!!
2650	!! ** Purpose : Get current state of KISS random number generator
2651	!!
2652	!! --------------------------------------------------------------------
2653	IMPLICIT NONE
2654	INTEGER(KIND=i8) :: ix, iy, iz, iw
2655
2656	ix = x
2657	iy = y
2658	iz = z
2659	iw = w
2660
2661	END SUBROUTINE kiss_state
2662
2663
2664	SUBROUTINE kiss_reset()
2665	!! --------------------------------------------------------------------
2666	!! * ROUTINE kiss_reset *
2667	!!
2668	!! ** Purpose : Reset the default seeds for KISS random number generator
2669	!!
2670	!! --------------------------------------------------------------------
2671	IMPLICIT NONE
2672
2673	x=1234567890987654321_8
2674	y=362436362436362436_8
2675	z=1066149217761810_8
2676	w=123456123456123456_8
2677
2678	END SUBROUTINE kiss_reset
2679
2680	SUBROUTINE kiss_uniform(uran)
2681	!! --------------------------------------------------------------------
2682	!! * ROUTINE kiss_uniform *
2683	!!
2684	!! ** Purpose : Real random numbers with uniform distribution in [0,1]
2685	!!
2686	!! --------------------------------------------------------------------
2687	IMPLICIT NONE
2688	REAL(KIND=wp) :: uran
2689
2690	uran = half * ( one + REAL(kiss(),wp) / huge64 )
2691
2692	END SUBROUTINE kiss_uniform
2693
2694
2695	SUBROUTINE kiss_gaussian(gran)
2696	!! --------------------------------------------------------------------
2697	!! * ROUTINE kiss_gaussian *
2698	!!
2699	!! ** Purpose : Real random numbers with Gaussian distribution N(0,1)
2700	!!
2701	!! ** Method : Generate 2 new Gaussian draws (gran1 and gran2)
2702	!! from 2 uniform draws on [-1,1] (u1 and u2),
2703	!! using the Marsaglia polar method
2704	!! (see Devroye, Non-Uniform Random Variate Generation, p. 235-236)
2705	!!
2706	!! --------------------------------------------------------------------
2707	IMPLICIT NONE
2708	REAL(KIND=wp) :: gran, u1, u2, rsq, fac
2709
2710	IF (ig.EQ.1) THEN
2711	rsq = two
2712	DO WHILE ( (rsq.GE.one).OR. (rsq.EQ.zero) )
2713	u1 = REAL(kiss(),wp) / huge64
2714	u2 = REAL(kiss(),wp) / huge64
2715	rsq = u1u1 + u2u2
2716	ENDDO
2717	fac = SQRT(-two*LOG(rsq)/rsq)
2718	gran1 = u1 * fac
2719	gran2 = u2 * fac
2720	ENDIF
2721
2722	! Output one of the 2 draws
2723	IF (ig.EQ.1) THEN
2724	gran = gran1 ; ig = 2
2725	ELSE
2726	gran = gran2 ; ig = 1
2727	ENDIF
2728
2729	END SUBROUTINE kiss_gaussian
2730
2731
2732	SUBROUTINE kiss_gamma(gamr,k)
2733	!! --------------------------------------------------------------------
2734	!! * ROUTINE kiss_gamma *
2735	!!
2736	!! ** Purpose : Real random numbers with Gamma distribution Gamma(k,1)
2737	!!
2738	!! --------------------------------------------------------------------
2739	IMPLICIT NONE
2740	REAL(KIND=wp), PARAMETER :: p1 = 4.5_8
2741	REAL(KIND=wp), PARAMETER :: p2 = 2.50407739677627_8 ! 1+LOG(9/2)
2742	REAL(KIND=wp), PARAMETER :: p3 = 1.38629436111989_8 ! LOG(4)
2743	REAL(KIND=wp) :: gamr, k, u1, u2, b, c, d, xx, yy, zz, rr, ee
2744	LOGICAL :: accepted
2745
2746	IF (k.GT.one) THEN
2747	! Cheng's rejection algorithm
2748	! (see Devroye, Non-Uniform Random Variate Generation, p. 413)
2749	b = k - p3 ; d = SQRT(two*k-one) ; c = k + d
2750
2751	accepted=.FALSE.
2752	DO WHILE (.NOT.accepted)
2753	CALL kiss_uniform(u1)
2754	yy = LOG(u1/(one-u1)) / d ! Mistake in Devroye: "* k" instead of "/ d"
2755	xx = k * EXP(yy)
2756	rr = b + c * yy - xx
2757	CALL kiss_uniform(u2)
2758	zz = u1 * u1 * u2
2759
2760	accepted = rr .GE. (zz*p1-p2)
2761	IF (.NOT.accepted) accepted = rr .GE. LOG(zz)
2762	ENDDO
2763
2764	gamr = xx
2765
2766	ELSEIF (k.LT.one) THEN
2767	! Rejection from the Weibull density
2768	! (see Devroye, Non-Uniform Random Variate Generation, p. 415)
2769	c = one/k ; d = (one-k) * EXP( (k/(one-k)) * LOG(k) )
2770
2771	accepted=.FALSE.
2772	DO WHILE (.NOT.accepted)
2773	CALL kiss_uniform(u1)
2774	zz = -LOG(u1)
2775	xx = EXP( c * LOG(zz) )
2776	CALL kiss_uniform(u2)
2777	ee = -LOG(u2)
2778
2779	accepted = (zz+ee) .GE. (d+xx) ! Mistake in Devroye: "LE" instead of "GE"
2780	ENDDO
2781
2782	gamr = xx
2783
2784	ELSE
2785	! Exponential distribution
2786	CALL kiss_uniform(u1)
2787	gamr = -LOG(u1)
2788
2789	ENDIF
2790
2791	END SUBROUTINE kiss_gamma
2792
2793
2794	SUBROUTINE kiss_sample(a,n,k)
2795	!! --------------------------------------------------------------------
2796	!! * ROUTINE kiss_sample *
2797	!!
2798	!! ** Purpose : Select a random sample of size k from a set of n integers
2799	!!
2800	!! ** Method : The sample is output in the first k elements of a
2801	!! Set k equal to n to obtain a random permutation
2802	!! of the whole set of integers
2803	!!
2804	!! --------------------------------------------------------------------
2805	IMPLICIT NONE
2806	INTEGER(KIND=i8), DIMENSION(:) :: a
2807	INTEGER(KIND=i8) :: n, k, i, j, atmp
2808	REAL(KIND=wp) :: uran
2809
2810	! Select the sample using the swapping method
2811	! (see Devroye, Non-Uniform Random Variate Generation, p. 612)
2812	DO i=1,k
2813	! Randomly select the swapping element between i and n (inclusive)
2814	CALL kiss_uniform(uran)
2815	j = i - 1 + CEILING( REAL(n-i+1,8) * uran )
2816	! Swap elements i and j
2817	atmp = a(i) ; a(i) = a(j) ; a(j) = atmp
2818	ENDDO
2819
2820	END SUBROUTINE kiss_sample
2821
2822	#ifdef NEMO_V34
2823	SUBROUTINE ctl_nam(ier,cstr,lout)
2824	INTEGER, INTENT(IN) :: ier
2825	LOGICAL, INTENT(IN) :: lout
2826	CHARACTER(LEN=*),INTENT(IN) :: cstr
2827	IF(lout) WRITE(numout,'(A,I4,X,A)') 'Error ',ier,TRIM(cstr)
2828	END SUBROUTINE
2829	#endif
2830
2831	END MODULE stopack

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_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/STO/stopack.F90 @ 15690

Download in other formats: