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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

stopack.F90 in branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/STO – NEMO

Context Navigation

source: branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/STO/stopack.F90 @ 13191

Last change on this file since 13191 was 13191, checked in by jwhile, 4 years ago
Updates for 1d runnig
File size: 114.7 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	#if defined key_traldf_c3d
627	REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj,jpk) :: coeff
628	REAL(wp), POINTER, DIMENSION(:,:,:) :: gauss
629	#elif defined key_traldf_c2d
630	REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj) :: coeff
631	REAL(wp), POINTER, DIMENSION(:,:) :: gauss
632	#elif defined key_traldf_c1d
633	REAL(wp), INTENT( inout ), DIMENSION(jpk) :: coeff
634	REAL(wp), POINTER, DIMENSION(:) :: gauss
635	#else
636	REAL(wp), INTENT( inout ) :: coeff
637	REAL(wp), POINTER :: gauss
638	#endif
639	INTEGER, INTENT( in ) :: nn_type
640	REAL(wp), INTENT( in ) :: rn_sd
641	INTEGER, INTENT( in ) :: kspp
642	INTEGER, INTENT( in ), OPTIONAL :: klev
643	REAL(wp), POINTER, DIMENSION(:,:) :: gauss
644	REAL(wp) :: zsd,xme,mm
645	CHARACTER (LEN=99) :: cstrng
646	INTEGER :: jklev
647
648	#if defined key_traldf_c2d \|\| key_traldf_c3d
649	CALL wrk_alloc(jpi,jpj,gauss)
650
651	IF( ln_spp_own_gauss ) THEN
652	gauss = gauss_n_2d_p
653	ELSE
654	gauss = rn_spp_stdev * gauss_n_2d / rn_sppt_stdev
655	ENDIF
656
657	IF( nn_type == 1 ) THEN
658	gauss = gauss * rn_sd
659	coeff = coeff * ( 1._wp + gauss )
660	ELSEIF ( nn_type == 2 ) THEN
661	zsd = rn_sd
662	xme = -0.5_wpzsdzsd
663	gauss = gauss * zsd + xme
664	coeff = exp(gauss) * coeff
665	ELSEIF ( nn_type == 3 ) THEN
666	zsd = rn_sd
667	xme = 0._wp
668	gauss = gauss * zsd + xme
669	coeff = exp(gauss) * coeff
670	ELSEIF( nn_type == 4 ) THEN
671	gauss = gauss * rn_sd
672	coeff = coeff * ( 1._wp + gauss )
673	#if defined key_traldf_c3d \|\| key_traldf_c2d \|\| key_traldf_c1d
674	WHERE( coeff > 1._wp ) coeff = 1._wp
675	WHERE( coeff < 0._wp ) coeff = 0._wp
676	#else
677	IF( coeff > 1._wp ) coeff = 1._wp
678	IF( coeff < 0._wp ) coeff = 0._wp
679	#endif
680	ELSEIF ( nn_type == 5 ) THEN
681	zsd = rn_sd
682	xme = -0.5_wpzsdzsd
683	gauss = gauss * zsd + xme
684	coeff = exp(gauss) * coeff
685	#if defined key_traldf_c3d \|\| key_traldf_c2d \|\| key_traldf_c1d
686	WHERE( coeff > 1._wp ) coeff = 1._wp
687	WHERE( coeff < 0._wp ) coeff = 0._wp
688	#else
689	IF( coeff > 1._wp ) coeff = 1._wp
690	IF( coeff < 0._wp ) coeff = 0._wp
691	#endif
692	ELSEIF ( nn_type == 6 ) THEN
693	zsd = rn_sd
694	xme = 0._wp
695	gauss = gauss * zsd + xme
696	coeff = exp(gauss) * coeff
697	#if defined key_traldf_c3d \|\| key_traldf_c2d \|\| key_traldf_c1d
698	WHERE( coeff > 1._wp ) coeff = 1._wp
699	WHERE( coeff < 0._wp ) coeff = 0._wp
700	#else
701	IF( coeff > 1._wp ) coeff = 1._wp
702	IF( coeff < 0._wp ) coeff = 0._wp
703	#endif
704	ELSE
705	CALL ctl_stop( 'spp dqdt wrong option')
706	ENDIF
707
708	#ifdef key_iomput
709	write(cstrng,'(A,I2.2)') 'spp_par',kspp
710	IF (iom_use(TRIM(cstrng)) ) CALL iom_put( TRIM(cstrng) , coeff )
711	#endif
712
713	IF( ln_stopack_diags ) THEN
714	IF(PRESENT(klev)) THEN
715	jklev = klev
716	ELSE
717	jklev = 0
718	ENDIF
719	CALL spp_stats(kt,kspp,jklev,coeff)
720	ENDIF
721
722	CALL wrk_dealloc(jpi,jpj,gauss)
723
724	#else
725	CALL ctl_stop( 'spp_gen: parameter perturbation will only work with '// &
726	'key_traldf_c2d or key_traldf_c3d')
727	#endif
728
729
730	END SUBROUTINE
731
732	SUBROUTINE spp_stats(mt,kp,kl,rcf)
733	!!----------------------------------------------------------------------
734	!! * ROUTINE spp_stats *
735	!!
736	!! ** Purpose : Compute and print basic SPP statistics
737	!!----------------------------------------------------------------------
738	IMPLICIT NONE
739	INTEGER, INTENT(IN) :: mt,kp,kl
740	#if defined key_traldf_c3d
741	REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj,jpk) :: rcf
742	#elif defined key_traldf_c2d
743	REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj) :: rcf
744	#elif defined key_traldf_c1d
745	REAL(wp), INTENT( inout ), DIMENSION(jpk) :: rcf
746	#else
747	REAL(wp), INTENT( inout ) :: rcf
748	#endif
749	REAL(wp) :: mi,ma
750	CHARACTER(LEN=16) :: cstr = ' '
751	SELECT CASE ( kp )
752	CASE( jk_spp_alb )
753	cstr = 'ALBEDO '
754	CASE( jk_spp_rhg )
755	cstr = 'ICE RHEOLOGY'
756	CASE( jk_spp_relw )
757	cstr = 'RELATIVE WND'
758	CASE( jk_spp_dqdt )
759	cstr = 'SST RELAXAT.'
760	CASE( jk_spp_deds )
761	cstr = 'SSS RELAXAT.'
762	CASE( jk_spp_arnf )
763	cstr = 'RIVER MIXING'
764	CASE( jk_spp_geot )
765	cstr = 'GEOTHERM.FLX'
766	CASE( jk_spp_qsi0 )
767	cstr = 'SOLAR EXTIN.'
768	CASE( jk_spp_bfr )
769	cstr = 'BOTTOM FRICT'
770	CASE( jk_spp_aevd )
771	cstr = 'EDDY VISCDIF'
772	CASE( jk_spp_avt )
773	cstr = 'VERT. DIFFUS'
774	CASE( jk_spp_avm )
775	cstr = 'VERT. VISCOS'
776	CASE( jk_spp_tkelc )
777	cstr = 'TKE LANGMUIR'
778	CASE( jk_spp_tkedf )
779	cstr = 'TKE RN_EDIFF'
780	CASE( jk_spp_tkeds )
781	cstr = 'TKE RN_EDISS'
782	CASE( jk_spp_tkebb )
783	cstr = 'TKE RN_EBB '
784	CASE( jk_spp_tkefr )
785	cstr = 'TKE RN_EFR '
786	CASE( jk_spp_ahtu )
787	cstr = 'TRALDF AHTU '
788	CASE( jk_spp_ahtv )
789	cstr = 'TRALDF AHTV '
790	CASE( jk_spp_ahtw )
791	cstr = 'TRALDF AHTW '
792	CASE( jk_spp_ahtt )
793	cstr = 'TRALDF AHTT '
794	CASE( jk_spp_ahubbl )
795	cstr = 'BBL DIFFUSU '
796	CASE( jk_spp_ahvbbl )
797	cstr = 'BBL DIFFUSV '
798	CASE( jk_spp_ahm1 )
799	cstr = 'DYNLDF AHM1 '
800	CASE( jk_spp_ahm2 )
801	cstr = 'DYNLDF AHM2 '
802	CASE( jk_spp_ahm3 )
803	cstr = 'DYNLDF AHM3 '
804	CASE( jk_spp_ahm4 )
805	cstr = 'DYNLDF AHM4 '
806	CASE DEFAULT
807	CALL ctl_stop('Unrecognized SPP parameter: add it or turn off diagnostics')
808	END SELECT
809	#if defined key_traldf_c3d \|\| key_traldf_c2d \|\| key_traldf_c1d
810	mi = MINVAL(rcf)
811	ma = MAXVAL(rcf)
812	#else
813	mi = rcf
814	ma = rcf
815	#endif
816	IF(lk_mpp) CALL mpp_min(mi)
817	IF(lk_mpp) CALL mpp_max(ma)
818	IF(lwp) THEN
819	IF ( kl > 0 ) write(cstr(14:16),'(I3.3)') kl
820	WRITE(numdiag,9300) lkt,cstr,mi,ma
821	ENDIF
822	9300 FORMAT(' it :', i8, ' ', A16, ' Min: ',d10.3 , ' Max: ',d10.3)
823	rn_mmax ( kp ) = MAX ( MAX( abs(mi), ma), rn_mmax ( kp ) )
824	END SUBROUTINE spp_stats
825
826	SUBROUTINE stopack_report
827	!!----------------------------------------------------------------------
828	!! * ROUTINE spp_stats *
829	!!
830	!! ** Purpose : Report at the end of the run
831	!!----------------------------------------------------------------------
832	IMPLICIT NONE
833	INTEGER :: jp, numrep
834	CHARACTER(LEN=16) :: cstr = ' '
835	REAL(wp) :: zmul
836
837	CALL ctl_opn(numrep, 'stopack.report', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
838
839	WRITE(numrep,'(A)') ' REPORT from STOPACK stochastic physics package'
840	WRITE(numrep,'(A)')
841	WRITE(numrep,'(A)') ' SPP part'
842
843	DO jp=1,jk_spp
844	SELECT CASE ( jp )
845	CASE( jk_spp_alb )
846	cstr = 'ALBEDO '
847	CASE( jk_spp_rhg )
848	cstr = 'ICE RHEOLOGY'
849	CASE( jk_spp_relw )
850	cstr = 'RELATIVE WND'
851	CASE( jk_spp_dqdt )
852	cstr = 'SST RELAXAT.'
853	CASE( jk_spp_deds )
854	cstr = 'SSS RELAXAT.'
855	CASE( jk_spp_arnf )
856	cstr = 'RIVER MIXING'
857	CASE( jk_spp_geot )
858	cstr = 'GEOTHERM.FLX'
859	CASE( jk_spp_qsi0 )
860	cstr = 'SOLAR EXTIN.'
861	CASE( jk_spp_bfr )
862	cstr = 'BOTTOM FRICT'
863	CASE( jk_spp_aevd )
864	cstr = 'EDDY VISCDIF'
865	CASE( jk_spp_avt )
866	cstr = 'VERT. DIFFUS'
867	CASE( jk_spp_avm )
868	cstr = 'VERT. VISCOS'
869	CASE( jk_spp_tkelc )
870	cstr = 'TKE LANGMUIR'
871	CASE( jk_spp_tkedf )
872	cstr = 'TKE RN_EDIFF'
873	CASE( jk_spp_tkeds )
874	cstr = 'TKE RN_EDISS'
875	CASE( jk_spp_tkebb )
876	cstr = 'TKE RN_EBB '
877	CASE( jk_spp_tkefr )
878	cstr = 'TKE RN_EFR '
879	CASE( jk_spp_ahtu )
880	cstr = 'TRALDF AHTU '
881	CASE( jk_spp_ahtv )
882	cstr = 'TRALDF AHTV '
883	CASE( jk_spp_ahtw )
884	cstr = 'TRALDF AHTW '
885	CASE( jk_spp_ahtt )
886	cstr = 'TRALDF AHTT '
887	CASE( jk_spp_ahubbl )
888	cstr = 'BBL DIFFUSU '
889	CASE( jk_spp_ahvbbl )
890	cstr = 'BBL DIFFUSV '
891	CASE( jk_spp_ahm1 )
892	cstr = 'DYNLDF AHM1 '
893	CASE( jk_spp_ahm2 )
894	cstr = 'DYNLDF AHM2 '
895	CASE( jk_spp_ahm3 )
896	cstr = 'DYNLDF AHM3 '
897	CASE( jk_spp_ahm4 )
898	cstr = 'DYNLDF AHM4 '
899	CASE DEFAULT
900	CALL ctl_stop('Unrecognized SPP parameter: add it or turn off diagnostics')
901	END SELECT
902	IF ( rn_mmax(jp) .GT. 0._wp ) &
903	WRITE(numrep,'(A,A,A,D10.3)') ' Maximum of absolute values for parameter ', trim(cstr), ' = ', rn_mmax(jp)
904	ENDDO
905
906	IF(sppt_step .eq. 2 ) THEN
907	zmul = rdt
908	ELSEIF (sppt_step .eq. 0 .or. sppt_step .eq. 1 ) THEN
909	zmul = 1._wp
910	ENDIF
911	rn_mmax(jk_sppt_tem:jk_sppt_vvl) = rn_mmax(jk_sppt_tem:jk_sppt_vvl) * zmul
912
913	WRITE(numrep,'(A)')
914	WRITE(numrep,'(A)') ' SPPT part'
915	WRITE(numrep,'(A,D10.3)') ' Maximum of absolute values for TEM ', rn_mmax(jk_sppt_tem)
916	IF( rn_mmax(jk_sppt_tem) > 0.5 ) WRITE(numrep,'(A)' ) 'Larger than 0.5, might be too big '
917	WRITE(numrep,'(A,D10.3)') ' Maximum of absolute values for SAL ', rn_mmax(jk_sppt_sal)
918	IF( rn_mmax(jk_sppt_sal) > 0.2 ) WRITE(numrep,'(A)' ) 'Larger than 0.2, might be too big '
919	WRITE(numrep,'(A,D10.3)') ' Maximum of absolute values for U-VEL', rn_mmax(jk_sppt_uvl)
920	IF( rn_mmax(jk_sppt_uvl) > 0.1 ) WRITE(numrep,'(A)' ) 'Larger than 0.1, might be too big '
921	WRITE(numrep,'(A,D10.3)') ' Maximum of absolute values for V-VEL', rn_mmax(jk_sppt_vvl)
922	IF( rn_mmax(jk_sppt_vvl) > 0.1 ) WRITE(numrep,'(A)' ) 'Larger than 0.1, might be too big '
923
924	WRITE(numrep,'(A)')
925	WRITE(numrep,'(A)') ' SKEB part'
926	WRITE(numrep,'(A,A,A,D10.3)') ' Maximum of absolute values for NUM ', trim(cstr), ' = ', rn_mmax(jk_skeb_dnum)
927	WRITE(numrep,'(A)' ) ' (Perturbation from numerical dissipation)'
928	IF( rn_mmax(jk_skeb_dnum) < 0.5e-4 ) WRITE(numrep,'(A)' ) ' Smaller than 0.5e-4, might be too small'
929	IF( rn_mmax(jk_skeb_dnum) > 0.2e-2 ) WRITE(numrep,'(A)' ) ' Larger than 0.2e-2, might be too big '
930	WRITE(numrep,'(A)')
931	WRITE(numrep,'(A,A,A,D10.3)') ' Maximum of absolute values for CONV ', trim(cstr), ' = ', rn_mmax(jk_skeb_dcon)
932	WRITE(numrep,'(A)' ) ' (Perturbation from convection dissipation)'
933	IF( rn_mmax(jk_skeb_dcon) < 0.5e-4 ) WRITE(numrep,'(A)' ) ' Smaller than 0.5e-4, might be too small'
934	IF( rn_mmax(jk_skeb_dcon) > 0.2e-2 ) WRITE(numrep,'(A)' ) ' Larger than 0.2e-2, might be too big '
935	WRITE(numrep,'(A)')
936	WRITE(numrep,'(A,A,A,D10.3)') ' Maximum of absolute values for TOTAL', trim(cstr), ' = ', rn_mmax(jk_skeb_tot)
937	WRITE(numrep,'(A)' ) ' (Perturbation from total energy dissipation)'
938	IF( rn_mmax(jk_skeb_tot) < 0.5e-4 ) WRITE(numrep,'(A)' ) ' Smaller than 0.5e-4, might be too small'
939	IF( rn_mmax(jk_skeb_tot) > 0.2e-2 ) WRITE(numrep,'(A)' ) ' Larger than 0.2e-2, might be too big '
940
941	CLOSE(numrep)
942	CLOSE(numdiag)
943
944	END SUBROUTINE stopack_report
945
946	SUBROUTINE spp_aht ( kt, coeff,nn_type, rn_sd, kspp )
947	!!----------------------------------------------------------------------
948	!! * ROUTINE spp_aht *
949	!!
950	!! ** Purpose : Perturbing diffusivity parameters
951	!! As spp_gen, but specifically designed for diffusivity
952	!! where coeff can be 2D or 3D
953	!!----------------------------------------------------------------------
954	INTEGER, INTENT( in ) :: kt
955	#if defined key_traldf_c3d
956	REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj,jpk) :: coeff
957	#elif defined key_traldf_c2d
958	REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj) :: coeff
959	#elif defined key_traldf_c1d
960	REAL(wp), INTENT( inout ), DIMENSION(jpk) :: coeff
961	#else
962	REAL(wp), INTENT( inout ) :: coeff
963	#endif
964	INTEGER, INTENT( in ) :: kspp
965	INTEGER, INTENT( in ) :: nn_type
966	REAL(wp), INTENT( in ) :: rn_sd
967	REAL(wp), POINTER, DIMENSION(:,:) :: gauss
968	REAL(wp) :: zsd,xme
969	INTEGER :: jk
970
971	CALL wrk_alloc(jpi,jpj,gauss)
972
973	IF( ln_spp_own_gauss ) THEN
974	gauss = gauss_n_2d_p
975	ELSE
976	gauss = rn_spp_stdev * gauss_n_2d / rn_sppt_stdev
977	ENDIF
978
979	IF( nn_type == 1 ) THEN
980	gauss = gauss * rn_sd
981	#if defined key_traldf_c3d
982	DO jk=1,jpk
983	coeff(:,:,jk) = coeff(:,:,jk) * ( 1._wp + gauss )
984	ENDDO
985	#elif defined key_traldf_c2d
986	coeff = coeff * ( 1._wp + gauss )
987	#endif
988	ELSEIF ( nn_type == 2 ) THEN
989	zsd = rn_sd
990	xme = -0.5_wpzsdzsd
991	gauss = gauss * zsd + xme
992	#if defined key_traldf_c3d
993	DO jk=1,jpk
994	coeff(:,:,jk) = exp(gauss) * coeff(:,:,jk)
995	ENDDOF
996	#elif defined key_traldf_c2d
997	coeff = exp(gauss) * coeff
998	#endif
999	ELSEIF ( nn_type == 3 ) THEN
1000	zsd = rn_sd
1001	xme = 0._wp
1002	gauss = gauss * zsd + xme
1003	#if defined key_traldf_c3d
1004	DO jk=1,jpk
1005	coeff(:,:,jk) = exp(gauss) * coeff(:,:,jk)
1006	ENDDOF
1007	#elif defined key_traldf_c2d
1008	coeff = exp(gauss) * coeff
1009	#endif
1010	ELSE
1011	CALL ctl_stop( 'spp aht wrong option')
1012	ENDIF
1013
1014	IF( ln_stopack_diags ) THEN
1015	CALL spp_stats(kt,kspp,0,coeff)
1016	ENDIF
1017
1018	CALL wrk_dealloc(jpi,jpj,gauss)
1019
1020	END SUBROUTINE
1021
1022	SUBROUTINE spp_ahm ( kt, coeff,nn_type, rn_sd, kspp )
1023	!!----------------------------------------------------------------------
1024	!! * ROUTINE spp_ahm *
1025	!!
1026	!! ** Purpose : Perturbing viscosity parameters
1027	!! As spp_gen, but specifically designed for viscosity
1028	!! where coeff can be 2D or 3D
1029	!!----------------------------------------------------------------------
1030	INTEGER, INTENT( in ) :: kt
1031	#if defined key_dynldf_c3d
1032	REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj,jpk) :: coeff
1033	#elif defined key_dynldf_c2d
1034	REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj) :: coeff
1035	#elif defined key_dynldf_c1d
1036	REAL(wp), INTENT( inout ), DIMENSION(jpk) :: coeff
1037	#else
1038	REAL(wp), INTENT( inout ) :: coeff
1039	#endif
1040	INTEGER, INTENT( in ) :: kspp
1041	INTEGER, INTENT( in ) :: nn_type
1042	REAL(wp), INTENT( in ) :: rn_sd
1043	REAL(wp), POINTER, DIMENSION(:,:) :: gauss
1044	REAL(wp) :: zsd,xme
1045	INTEGER :: jk
1046
1047	CALL wrk_alloc(jpi,jpj,gauss)
1048
1049	IF( ln_spp_own_gauss ) THEN
1050	gauss = gauss_n_2d_p
1051	ELSE
1052	gauss = rn_spp_stdev * gauss_n_2d / rn_sppt_stdev
1053	ENDIF
1054
1055	IF( nn_type == 1 ) THEN
1056	gauss = gauss * rn_sd
1057	#if defined key_dynldf_c3d
1058	DO jk=1,jpk
1059	coeff(:,:,jk) = coeff(:,:,jk) * ( 1._wp + gauss )
1060	ENDDO
1061	#elif defined key_dynldf_c2d
1062	coeff = coeff * ( 1._wp + gauss )
1063	#endif
1064	ELSEIF ( nn_type == 2 ) THEN
1065	zsd = rn_sd
1066	xme = -0.5_wpzsdzsd
1067	gauss = gauss * zsd + xme
1068	#if defined key_dynldf_c3d
1069	DO jk=1,jpk
1070	coeff(:,:,jk) = exp(gauss) * coeff(:,:,jk)
1071	ENDDO
1072	#elif defined key_dynldf_c2d
1073	coeff = exp(gauss) * coeff
1074	#endif
1075	ELSEIF ( nn_type == 3 ) THEN
1076	zsd = rn_sd
1077	xme = 0._wp
1078	gauss = gauss * zsd + xme
1079	#if defined key_dynldf_c3d
1080	DO jk=1,jpk
1081	coeff(:,:,jk) = exp(gauss) * coeff(:,:,jk)
1082	ENDDO
1083	#elif defined key_dynldf_c2d
1084	coeff = exp(gauss) * coeff
1085	#endif
1086	ELSE
1087	CALL ctl_stop( 'spp ahm wrong option')
1088	ENDIF
1089
1090	IF( ln_stopack_diags ) THEN
1091	CALL spp_stats(kt,kspp,0,coeff)
1092	ENDIF
1093
1094	CALL wrk_dealloc(jpi,jpj,gauss)
1095
1096	END SUBROUTINE
1097
1098	SUBROUTINE tra_sppt_apply ( kt , ks)
1099	!!----------------------------------------------------------------------
1100	!! * ROUTINE tra_sppt_apply *
1101	!!
1102	!! ** Purpose : Apply perturbation to tracers
1103	!! Step 0/1 for tendencies, step 2 for variables
1104	!! after timestepping
1105	!!----------------------------------------------------------------------
1106
1107	INTEGER, INTENT( in ) :: kt, ks ! Main time step counter
1108	REAL(wp) :: mi,ma
1109
1110	IF( ks .NE. sppt_step ) RETURN
1111
1112	IF(lwp) THEN
1113	WRITE(numout,*)
1114	WRITE(numout,*) ' Applying sppt on tracers at timestep ', kt
1115	WRITE(numout,*)
1116	ENDIF
1117
1118	! Modify the tendencies
1119	IF(sppt_step .eq. 2 ) THEN
1120	! At the end of the timestepping, after array swapping
1121	tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) + gauss_n * spptt * rdt
1122	tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) + gauss_n * sppts * rdt
1123	tsb(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) + gauss_n * spptt * rdt
1124	tsb(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) + gauss_n * sppts * rdt
1125	IF( ln_sppt_glocon ) CALL sppt_glocon( tsn(:,:,:,jp_tem), 'T' )
1126	IF( ln_sppt_glocon ) CALL sppt_glocon( tsn(:,:,:,jp_sal), 'T' )
1127	IF( ln_sppt_glocon ) CALL sppt_glocon( tsb(:,:,:,jp_tem), 'T' )
1128	IF( ln_sppt_glocon ) CALL sppt_glocon( tsb(:,:,:,jp_sal), 'T' )
1129	CALL lbc_lnk( tsb(:,:,:,jp_tem) , 'T', 1._wp)
1130	CALL lbc_lnk( tsb(:,:,:,jp_sal) , 'T', 1._wp)
1131	CALL lbc_lnk( tsn(:,:,:,jp_tem) , 'T', 1._wp)
1132	CALL lbc_lnk( tsn(:,:,:,jp_sal) , 'T', 1._wp)
1133	ELSEIF (sppt_step .eq. 0 .or. sppt_step .eq. 1 ) THEN
1134	! At the beginning / before vertical diffusion
1135	tsa(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) + gauss_n * spptt
1136	tsa(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) + gauss_n * sppts
1137	IF( ln_sppt_glocon ) CALL sppt_glocon( tsa(:,:,:,jp_tem), 'T' )
1138	IF( ln_sppt_glocon ) CALL sppt_glocon( tsa(:,:,:,jp_sal), 'T' )
1139	CALL lbc_lnk( tsa(:,:,:,jp_tem) , 'T', 1._wp)
1140	CALL lbc_lnk( tsa(:,:,:,jp_sal) , 'T', 1._wp)
1141	ENDIF
1142
1143	IF( ln_stopack_diags ) THEN
1144	mi = MINVAL(spptt)
1145	ma = MAXVAL(spptt)
1146	IF(lk_mpp) CALL mpp_min(mi)
1147	IF(lk_mpp) CALL mpp_max(ma)
1148	IF(lwp) WRITE(numdiag,9301) lkt,'SPPT TEMPERATURE' ,mi,ma
1149	rn_mmax ( jk_sppt_tem ) = MAX ( MAX( abs(mi), ma), rn_mmax ( jk_sppt_tem ) )
1150
1151	mi = MINVAL(sppts)
1152	ma = MAXVAL(sppts)
1153	IF(lk_mpp) CALL mpp_min(mi)
1154	IF(lk_mpp) CALL mpp_max(ma)
1155	IF(lwp) WRITE(numdiag,9301) lkt,'SPPT SALINITY ' ,mi,ma
1156	9301 FORMAT(' it :', i8, ' ', A16, ' Min: ',d10.3 , ' Max: ',d10.3)
1157	rn_mmax ( jk_sppt_sal ) = MAX ( MAX( abs(mi), ma), rn_mmax ( jk_sppt_sal ) )
1158	ENDIF
1159	! Reset the tendencies
1160	spptt = 0._wp ; sppts = 0._wp
1161
1162	END SUBROUTINE tra_sppt_apply
1163
1164	SUBROUTINE dyn_sppt_apply ( kt , ks)
1165	!!----------------------------------------------------------------------
1166	!! * ROUTINE dyn_sppt_apply *
1167	!!
1168	!! ** Purpose : Apply perturbation to momentum
1169	!! Step 0/1 for tendencies, step 2 for variables
1170	!! after timestepping
1171	!!----------------------------------------------------------------------
1172
1173	INTEGER, INTENT( in ) :: kt, ks ! Main time step counter
1174	REAL(wp) :: mi,ma
1175
1176	IF( ks .NE. sppt_step ) RETURN
1177
1178	IF(lwp) THEN
1179	WRITE(numout,*)
1180	WRITE(numout,*) ' Applying sppt on currents at timestep ', kt
1181	WRITE(numout,*)
1182	ENDIF
1183
1184	! Modify the tendencies
1185	IF(sppt_step .eq. 2 ) THEN
1186	! At the end of the timestepping, after array swapping
1187	un(:,:,:) = un(:,:,:) + gauss_nu* spptu * rdt
1188	vn(:,:,:) = vn(:,:,:) + gauss_nv* spptv * rdt
1189	ub(:,:,:) = ub(:,:,:) + gauss_nu* spptu * rdt
1190	vb(:,:,:) = vb(:,:,:) + gauss_nv* spptv * rdt
1191	IF( ln_sppt_glocon ) CALL sppt_glocon( ub, 'U' )
1192	IF( ln_sppt_glocon ) CALL sppt_glocon( vb, 'V' )
1193	IF( ln_sppt_glocon ) CALL sppt_glocon( un, 'U' )
1194	IF( ln_sppt_glocon ) CALL sppt_glocon( vn, 'V' )
1195	CALL lbc_lnk( un , 'U', -1._wp)
1196	CALL lbc_lnk( vn , 'V', -1._wp)
1197	CALL lbc_lnk( ub , 'U', -1._wp)
1198	CALL lbc_lnk( vb , 'V', -1._wp)
1199	ELSEIF (sppt_step .eq. 0 .or. sppt_step .eq. 1 ) THEN
1200	! At the beginning / before vertical diffusion
1201	ua(:,:,:) = ua(:,:,:) + gauss_nu* spptu
1202	va(:,:,:) = va(:,:,:) + gauss_nv* spptv
1203	IF( ln_sppt_glocon ) CALL sppt_glocon( ua, 'U' )
1204	IF( ln_sppt_glocon ) CALL sppt_glocon( va, 'V' )
1205	CALL lbc_lnk( ua , 'U', -1._wp)
1206	CALL lbc_lnk( va , 'V', -1._wp)
1207	ENDIF
1208
1209	IF( ln_stopack_diags ) THEN
1210	mi = MINVAL(spptu)
1211	ma = MAXVAL(spptu)
1212	IF(lk_mpp) CALL mpp_min(mi)
1213	IF(lk_mpp) CALL mpp_max(ma)
1214	IF(lwp) WRITE(numdiag,9301) lkt,'SPPT ZONAL CURR.' ,mi,ma
1215	rn_mmax ( jk_sppt_uvl ) = MAX ( MAX( abs(mi), ma), rn_mmax ( jk_sppt_uvl ) )
1216
1217	mi = MINVAL(spptv)
1218	ma = MAXVAL(spptv)
1219	IF(lk_mpp) CALL mpp_min(mi)
1220	IF(lk_mpp) CALL mpp_max(ma)
1221	IF(lwp) WRITE(numdiag,9301) lkt,'SPPT MERID CURR.' ,mi,ma
1222	9301 FORMAT(' it :', i8, ' ', A16, ' Min: ',d10.3 , ' Max: ',d10.3)
1223	rn_mmax ( jk_sppt_vvl ) = MAX ( MAX( abs(mi), ma), rn_mmax ( jk_sppt_vvl ) )
1224	ENDIF
1225
1226	! Reset the tendencies
1227	spptu = 0._wp ; spptv = 0._wp
1228
1229	END SUBROUTINE dyn_sppt_apply
1230
1231	SUBROUTINE sppt_glocon ( zts, cd_type )
1232	!!----------------------------------------------------------------------
1233	!! * ROUTINE sppt_glocon *
1234	!!
1235	!! ** Purpose : Apply global conservation constraint
1236	!!----------------------------------------------------------------------
1237	REAL(wp), INTENT(INOUT) :: zts(jpi,jpj,jpk)
1238	CHARACTER(len=1), INTENT(IN) :: cd_type
1239	INTEGER :: ji,jj,jk
1240	REAL(wp) :: zv, vl
1241	! Calculate volume tendencies and renormalize
1242	! Note: sshn should be staggered before being used.
1243	SELECT CASE ( cd_type )
1244	CASE ( 'T' )
1245	jk=1
1246	zv = SUM( tmask_i(:,:)tmask(:,:,jk)e1t(:,:)e2t(:,:)sshn(:,:)*zts(:,:,1) )
1247	vl = SUM( tmask_i(:,:)tmask(:,:,jk)e1t(:,:)e2t(:,:)sshn(:,:) )
1248	DO jk = 1, jpkm1
1249	DO jj=1,jpj
1250	DO ji=1,jpi
1251	zv = zv + SUM( tmask_i(:,:)tmask(:,:,jk)e1t(:,:)e2t(:,:)fse3t(:,:,jk)*zts(:,:,jk) )
1252	vl = vl + SUM( tmask_i(:,:)tmask(:,:,jk)e1t(:,:)e2t(:,:)fse3t(:,:,jk) )
1253	END DO
1254	END DO
1255	END DO
1256	IF(lk_mpp) CALL mpp_sum(zv)
1257	IF(lk_mpp) CALL mpp_sum(vl)
1258	zv = zv / vl
1259	zts = zts - zv*tmask
1260	CASE ( 'U' )
1261	jk=1
1262	#if defined NEMO_V34
1263	zv = SUM( umask(:,:,jk)e1u(:,:)e2u(:,:)sshn(:,:)zts(:,:,1) )
1264	vl = SUM( umask(:,:,jk)e1u(:,:)e2u(:,:)*sshn(:,:) )
1265	#else
1266	zv = SUM( umask_i(:,:)umask(:,:,jk)e1u(:,:)e2u(:,:)sshn(:,:)*zts(:,:,1) )
1267	vl = SUM( umask_i(:,:)umask(:,:,jk)e1u(:,:)e2u(:,:)sshn(:,:) )
1268	#endif
1269	DO jk = 1, jpkm1
1270	DO jj=1,jpj
1271	DO ji=1,jpi
1272	#if defined NEMO_V34
1273	zv = zv + SUM( umask(:,:,jk)e1u(:,:)e2u(:,:)fse3u(:,:,jk)zts(:,:,jk) )
1274	vl = vl + SUM( umask(:,:,jk)e1u(:,:)e2u(:,:)*fse3u(:,:,jk) )
1275	#else
1276	zv = zv + SUM( umask_i(:,:)umask(:,:,jk)e1u(:,:)e2u(:,:)fse3u(:,:,jk)*zts(:,:,jk) )
1277	vl = vl + SUM( umask_i(:,:)umask(:,:,jk)e1u(:,:)e2u(:,:)fse3u(:,:,jk) )
1278	#endif
1279	END DO
1280	END DO
1281	END DO
1282	IF(lk_mpp) CALL mpp_sum(zv)
1283	IF(lk_mpp) CALL mpp_sum(vl)
1284	zv = zv / vl
1285	zts = zts - zv*umask
1286	CASE ( 'V' )
1287	jk=1
1288	#if defined NEMO_V34
1289	zv = SUM( vmask(:,:,jk)e1v(:,:)e2v(:,:)sshn(:,:)zts(:,:,1) )
1290	vl = SUM( vmask(:,:,jk)e1v(:,:)e2v(:,:)*sshn(:,:) )
1291	#else
1292	zv = SUM( vmask_i(:,:)vmask(:,:,jk)e1v(:,:)e2v(:,:)sshn(:,:)*zts(:,:,1) )
1293	vl = SUM( vmask_i(:,:)vmask(:,:,jk)e1v(:,:)e2v(:,:)sshn(:,:) )
1294	#endif
1295	DO jk = 1, jpkm1
1296	DO jj=1,jpj
1297	DO ji=1,jpi
1298	#if defined NEMO_V34
1299	zv = zv + SUM( vmask(:,:,jk)e1v(:,:)e2v(:,:)fse3v(:,:,jk)zts(:,:,jk) )
1300	vl = vl + SUM( vmask(:,:,jk)e1v(:,:)e2v(:,:)*fse3v(:,:,jk) )
1301	#else
1302	zv = zv + SUM( vmask_i(:,:)vmask(:,:,jk)e1v(:,:)e2v(:,:)fse3v(:,:,jk)*zts(:,:,jk) )
1303	vl = vl + SUM( vmask_i(:,:)vmask(:,:,jk)e1v(:,:)e2v(:,:)fse3v(:,:,jk) )
1304	#endif
1305	END DO
1306	END DO
1307	END DO
1308	IF(lk_mpp) CALL mpp_sum(zv)
1309	IF(lk_mpp) CALL mpp_sum(vl)
1310	zv = zv / vl
1311	zts = zts - zv*vmask
1312	ENDSELECT
1313
1314	END SUBROUTINE sppt_glocon
1315
1316	SUBROUTINE gaussian_ar1_field ( gn, nk, wei, gb, a, b,gn0, fltf, istep)
1317	!!----------------------------------------------------------------------
1318	!! * ROUTINE gaussian_ar1_field *
1319	!!
1320	!! ** Purpose : Generate correlated perturbation field
1321	!!----------------------------------------------------------------------
1322	REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: a, b
1323	REAL(wp), DIMENSION(jpk) , INTENT(INOUT) :: wei
1324	REAL(wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: gb
1325	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(INOUT) :: gn
1326	REAL(wp), DIMENSION(jpi,jpj ), INTENT(OUT) :: gn0
1327	INTEGER, INTENT(IN) :: nk,istep
1328	REAL(wp), INTENT(IN) :: fltf
1329	REAL(wp) :: g2d(jpi,jpj)
1330	INTEGER :: jf, jk, ji, jj, nbl
1331
1332	! Random noise on 2d field
1333	IF ( istep == 1 ) THEN
1334	CALL sppt_rand2d( g2d )
1335	CALL lbc_lnk( g2d , 'T', 1._wp)
1336	g2d_save = g2d
1337	IF( nn_sppt_step_bound .EQ. 1 ) THEN
1338	WHERE ( g2d < -ABS(rn_sppt_bound) ) g2d = -ABS(rn_sppt_bound)
1339	WHERE ( g2d > ABS(rn_sppt_bound) ) g2d = ABS(rn_sppt_bound)
1340	ENDIF
1341	ELSEIF ( istep == 2 ) THEN
1342	g2d = rn_spp_stdev * g2d_save / rn_sppt_stdev
1343	ELSEIF ( istep == 3 ) THEN
1344	g2d = rn_skeb_stdev * g2d_save / rn_sppt_stdev
1345	ENDIF
1346
1347	! Laplacian filter and re-normalization
1348	DO jf = 1, nk
1349	CALL sppt_flt( g2d )
1350	END DO
1351	g2d = g2d * fltf
1352
1353	#ifdef key_iommput
1354	! Output the random field
1355	IF(istep==1) THEN
1356	CALL iom_put( "sppt_ran" , g2d )
1357	ELSEIF (istep==2) THEN
1358	CALL iom_put( "spp_ran" , g2d )
1359	ELSEIF (istep==3) THEN
1360	CALL iom_put( "skeb_ran" , g2d )
1361	ENDIF
1362	#endif
1363
1364	! AR(1) process and array swap
1365	g2d = agb + bg2d
1366	gb = g2d
1367	gn0 = g2d * zdc(:,:,1)
1368
1369	IF (istep==2 .or. istep==3) RETURN
1370
1371	! From 2- to 3-d and vertical weigth
1372	IF(nn_vwei .eq. 2) THEN
1373	DO jj=1,jpj
1374	DO ji=1,jpi
1375	nbl = mbathy(ji,jj)
1376	wei(:) = 0._wp
1377	IF(nbl>0) THEN
1378	wei(1:nbl) = 1._wp
1379	wei(1) = 0._wp
1380	wei(2) = 0.5_wp
1381	wei(nbl-1) = 0.5_wp
1382	wei(nbl) = 0._wp
1383	DO jk=1,jpk
1384	gn(ji,jj,jk) = g2d(ji,jj) * wei(jk) * zdc(ji,jj,jk)
1385	ENDDO
1386	ENDIF
1387	ENDDO
1388	ENDDO
1389	ELSEIF(nn_vwei .eq. 3) THEN
1390	DO jj=1,jpj
1391	DO ji=1,jpi
1392	nbl = mbathy(ji,jj)
1393	wei(:) = 0._wp
1394	IF(nbl>0) THEN
1395	wei(1:nbl) = 1._wp
1396	wei(nbl-1) = 0.5_wp
1397	wei(nbl) = 0._wp
1398	DO jk=1,jpk
1399	gn(ji,jj,jk) = g2d(ji,jj) * wei(jk) * zdc(ji,jj,jk)
1400	ENDDO
1401	ENDIF
1402	ENDDO
1403	ENDDO
1404	ELSE
1405	DO jk=1,jpk
1406	gn(:,:,jk) = g2d * wei(jk) * zdc(:,:,jk)
1407	ENDDO
1408	ENDIF
1409
1410	! Bound
1411	IF( nn_sppt_step_bound .EQ. 2 ) THEN
1412	WHERE ( gn < -ABS(rn_sppt_bound) ) gn = -ABS(rn_sppt_bound)
1413	WHERE ( gn > ABS(rn_sppt_bound) ) gn = ABS(rn_sppt_bound)
1414	WHERE ( gn0< -ABS(rn_sppt_bound) ) gn0= -ABS(rn_sppt_bound)
1415	WHERE ( gn0> ABS(rn_sppt_bound) ) gn0= ABS(rn_sppt_bound)
1416	ENDIF
1417
1418	END SUBROUTINE gaussian_ar1_field
1419	!
1420	SUBROUTINE stopack_init
1421	!!----------------------------------------------------------------------
1422	!! * ROUTINE stopack_init *
1423	!!
1424	!! ** Purpose : Initialize stopack
1425	!!----------------------------------------------------------------------
1426	!!
1427	INTEGER :: ios, inum
1428	INTEGER :: nn_sppt_tra, nn_sppt_dyn, nn_spp_aht, nn_sppt_ice
1429	INTEGER :: ivals(8)
1430	INTEGER(KIND=8) :: ziseed(4) ! RNG seeds in integer type
1431	REAL(KIND=8) :: zrseed(4) ! RNG seeds in real type
1432	REAL(KIND=8) :: rdt_sppt
1433	!!
1434	NAMELIST/namstopack/ ln_stopack, ln_sppt_taumap, rn_sppt_tau, &
1435	ln_stopack_restart, sppt_filter_pass, rn_sppt_bound, sppt_step, nn_deftau,rn_sppt_stdev,&
1436	ln_distcoast, rn_distcoast, nn_vwei, nn_stopack_seed, nn_sppt_step_bound, nn_rndm_freq, &
1437	ln_sppt_glocon, ln_stopack_repr, &
1438	rn_icestr_sd, rn_icealb_sd, &
1439	ln_sppt_traxad, ln_sppt_trayad, ln_sppt_trazad, ln_sppt_trasad, ln_sppt_traldf, &
1440	ln_sppt_trazdf, ln_sppt_trazdfp,ln_sppt_traevd, ln_sppt_trabbc, ln_sppt_trabbl, &
1441	ln_sppt_tranpc, ln_sppt_tradmp, ln_sppt_traqsr, ln_sppt_transr, ln_sppt_traatf ,&
1442	ln_sppt_dynhpg, ln_sppt_dynspg, ln_sppt_dynkeg, ln_sppt_dynrvo, ln_sppt_dynpvo, &
1443	ln_sppt_dynzad, ln_sppt_dynldf, ln_sppt_dynzdf, ln_sppt_dynbfr, ln_sppt_dynatf, ln_sppt_dyntau,&
1444	ln_sppt_icehdf, ln_sppt_icelat, ln_sppt_icezdf, & !ln_sppt_icealb, ln_sppt_icestr
1445	nn_spp_icealb, nn_spp_icestr,&
1446	spp_filter_pass,rn_spp_stdev,rn_spp_tau,&
1447	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,&
1448	nn_spp_ahtw,nn_spp_ahtt,nn_spp_relw,nn_spp_arnf,nn_spp_geot,nn_spp_aevd,nn_spp_ahubbl,nn_spp_ahvbbl,&
1449	nn_spp_tkelc,nn_spp_tkedf,nn_spp_tkeds,nn_spp_tkebb,nn_spp_tkefr,&
1450	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,&
1451	rn_ahtw_sd,rn_ahtt_sd, rn_relw_sd, rn_arnf_sd,rn_geot_sd, rn_aevd_sd,rn_ahubbl_sd,rn_ahvbbl_sd,&
1452	rn_tkelc_sd,rn_tkedf_sd,rn_tkeds_sd,rn_tkebb_sd,rn_tkefr_sd,&
1453	ln_skeb,rn_skeb,nn_dcom_freq,rn_kh,rn_kc,ln_stopack_diags,skeb_filter_pass,rn_skeb_stdev,rn_skeb_tau,&
1454	nn_dconv,rn_beta_num, rn_beta_con
1455
1456	! Default values
1457	rn_sppt_bound = 3._wp
1458	ln_stopack = .false.
1459	ln_sppt_taumap = .false.
1460	rn_sppt_tau = 86400._wp * 5._wp
1461	sppt_filter_pass = 30
1462	nn_vwei = 0
1463	ln_distcoast = .false.
1464	ln_sppt_glocon = .false.
1465	rn_distcoast = 10._wp
1466	ln_stopack_restart = .true.
1467	nn_stopack_seed = (/1,2,3,narea/)
1468	nn_rndm_freq = 1
1469	nn_sppt_step_bound = 2
1470	nn_deftau = 1
1471	ln_sppt_traxad = .false. ! Switch for x advection
1472	ln_sppt_trayad = .false. ! Switch for y advection
1473	ln_sppt_trazad = .false. ! Switch for z advection
1474	ln_sppt_trasad = .false. ! Switch for z advection (s- case)
1475	ln_sppt_traldf = .false. ! Switch for lateral diffusion
1476	ln_sppt_trazdf = .false. ! Switch for vertical diffusion
1477	ln_sppt_trazdfp = .false. ! Switch for pure vertical diffusion
1478	ln_sppt_traevd = .false. ! Switch for enhanced vertical diffusion
1479	ln_sppt_trabbc = .false. ! Switch for bottom boundary condition
1480	ln_sppt_trabbl = .false. ! Switch for bottom boundary layer
1481	ln_sppt_tranpc = .false. ! Switch for non-penetrative convection
1482	ln_sppt_tradmp = .false. ! Switch for tracer damping
1483	ln_sppt_traqsr = .false. ! Switch for solar radiation
1484	ln_sppt_transr = .false. ! Switch for non-solar radiation / freshwater flux
1485	ln_sppt_traatf = .false. ! Switch for Asselin time-filter
1486
1487	ln_sppt_dynhpg = .false. ! Switch for hydrost. press. grad.
1488	ln_sppt_dynspg = .false. ! Switch for surface press. grad.
1489	ln_sppt_dynkeg = .false. ! Switch for horiz. advcetion
1490	ln_sppt_dynrvo = .false. ! Switch for Relative vorticity
1491	ln_sppt_dynpvo = .false. ! Switch for planetary vortic.
1492	ln_sppt_dynzad = .false. ! Switch for vertical advection
1493	ln_sppt_dynldf = .false. ! Switch for lateral viscosity
1494	ln_sppt_dynzdf = .false. ! Switch for vertical viscosity
1495	ln_sppt_dynbfr = .false. ! Switch for bottom friction
1496	ln_sppt_dynatf = .false. ! Switch for Asselin filter
1497	ln_sppt_dyntau = .false. ! Switch for wind stress
1498
1499	ln_sppt_icehdf = .false.
1500	ln_sppt_icelat = .false.
1501	ln_sppt_icezdf = .false.
1502
1503	nn_spp_icestr = 0
1504	nn_spp_icealb = 0
1505	rn_icestr_sd = 0.30_wp
1506	rn_icealb_sd = 0.30_wp
1507
1508	nn_spp_bfr =0
1509	nn_spp_dqdt =0
1510	nn_spp_arnf =0
1511	nn_spp_aevd =0
1512	nn_spp_geot =0
1513	nn_spp_dedt =0
1514	nn_spp_avt =0
1515	nn_spp_avm =0
1516	nn_spp_qsi0 =0
1517	nn_spp_relw =0
1518	nn_spp_tkelc =0
1519	nn_spp_tkedf =0
1520	nn_spp_tkeds =0
1521	nn_spp_tkebb =0
1522	nn_spp_tkefr =0
1523
1524	ln_skeb = .false.
1525	ln_stopack_diags = .false.
1526	rn_skeb_stdev = 1.0_wp
1527	skeb_filter_pass = 50
1528	rn_skeb_tau = 50
1529
1530	#ifdef NEMO_V34
1531	REWIND( numnam )
1532	READ ( numnam, namstopack )
1533	#else
1534	REWIND( numnam_ref )
1535	READ ( numnam_ref, namstopack, IOSTAT = ios, ERR = 901)
1536	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namstopack in reference namelist', lwp )
1537
1538	REWIND( numnam_cfg )
1539	READ ( numnam_cfg, namstopack, IOSTAT = ios, ERR = 902 )
1540	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namstopack in configuration namelist', lwp )
1541	IF(lwm) WRITE ( numond, namstopack )
1542	#endif
1543
1544	IF(lwp) THEN
1545	WRITE(numout,*)
1546	WRITE(numout,*) 'init_stopack : Stochastic physics package'
1547	WRITE(numout,*) '~~~~~~~~~~~~'
1548	WRITE(numout,*) ' Namelist namstopack: '
1549	WRITE(numout,*)
1550	WRITE(numout,*) ' Switch on STOPACK ln_stopack : ',ln_stopack
1551	WRITE(numout,*) ' Verbose diagnostics for STOPACK ln_stopack_diags:', ln_stopack_diags
1552	WRITE(numout,*)
1553	WRITE(numout,) ' ** Random generation'
1554	WRITE(numout,*) ' Read from restart file previous perturbation ln_stopack_restart :', ln_stopack_restart
1555	WRITE(numout,*) ' Frequency of calls of the SPPT perturbation field nn_rndm_freq :', nn_rndm_freq
1556	WRITE(numout,*) ' Reproducibility of random generation ln_stopack_repr :', ln_stopack_repr
1557	WRITE(numout,*) ' Seed for random number generator (no restart) nn_stopack_seed :', nn_stopack_seed(:)
1558	WRITE(numout,*)
1559	WRITE(numout,) ' ** SPPT scheme '
1560	WRITE(numout,*) ' SPPT step (0: beginning; 1: before ZVD; 2: after ZVD) sppt_step : ',sppt_step
1561	WRITE(numout,*) ' Use map of decorr. time scale ln_sppt_taumap :', ln_sppt_taumap
1562	WRITE(numout,*) ' If ln_sppt_taumap FALSE, use this constant (in days) rn_sppt_tau :', rn_sppt_tau
1563	WRITE(numout,*) ' Number of filter passes to correlate in space sppt_filter_pass :', sppt_filter_pass
1564	WRITE(numout,*) ' Standard deviation of the white noise rn_sppt_stdev :', rn_sppt_stdev
1565	WRITE(numout,*) ' Apply global conservation constraints ln_sppt_glocon :', ln_sppt_glocon
1566	WRITE(numout,*)
1567	WRITE(numout,*) ' Perturbation on tracers:'
1568	WRITE(numout,*) ' Switch for x advection ln_sppt_traxad :', ln_sppt_traxad
1569	WRITE(numout,*) ' Switch for y advection ln_sppt_trayad :', ln_sppt_trayad
1570	WRITE(numout,*) ' Switch for z advection ln_sppt_trazad :', ln_sppt_trazad
1571	WRITE(numout,*) ' Switch for z advection (s- case) ln_sppt_trasad :', ln_sppt_trasad
1572	WRITE(numout,*) ' Switch for lateral diffusion ln_sppt_traldf :', ln_sppt_traldf
1573	WRITE(numout,*) ' Switch for vertical diffusion ln_sppt_trazdf :', ln_sppt_trazdf
1574	WRITE(numout,*) ' Switch for pure vertical diffusion ln_sppt_trazdfp :', ln_sppt_trazdfp
1575	WRITE(numout,*) ' Switch for enhanced vertical diffusion ln_sppt_traevd :', ln_sppt_traevd
1576	WRITE(numout,*) ' Switch for bottom boundary condition ln_sppt_trabbc :', ln_sppt_trabbc
1577	WRITE(numout,*) ' Switch for bottom boundary layer ln_sppt_trabbl :', ln_sppt_trabbl
1578	WRITE(numout,*) ' Switch for non-penetrative convection ln_sppt_tranpc :', ln_sppt_tranpc
1579	WRITE(numout,*) ' Switch for tracer damping ln_sppt_tradmp :', ln_sppt_tradmp
1580	WRITE(numout,*) ' Switch for solar radiation ln_sppt_traqsr :', ln_sppt_traqsr
1581	WRITE(numout,*) ' Switch for non-solar rad. / freshwater flx flux ln_sppt_transr :', ln_sppt_transr
1582	WRITE(numout,*) ' Switch for Asselin time-filter ln_sppt_traatf :', ln_sppt_traatf
1583	WRITE(numout,*)
1584	WRITE(numout,*) ' Perturbation on dynamics:'
1585	WRITE(numout,*) ' Switch for horiz advection ln_sppt_dynkeg :', ln_sppt_dynkeg
1586	WRITE(numout,*) ' Switch for z advection ln_sppt_dynzad :', ln_sppt_dynzad
1587	WRITE(numout,*) ' Switch for wind stress ln_sppt_dyntau :', ln_sppt_dyntau
1588	WRITE(numout,*) ' Switch for lateral diffusion ln_sppt_dynldf :', ln_sppt_dynldf
1589	WRITE(numout,*) ' Switch for vertical diffusion ln_sppt_dynzdf :', ln_sppt_dynzdf
1590	WRITE(numout,*) ' Switch for planet. vorticity ln_sppt_dynpvo :', ln_sppt_dynpvo
1591	WRITE(numout,*) ' Switch for relative vorticity ln_sppt_dynrvo :', ln_sppt_dynrvo
1592	WRITE(numout,*) ' Switch for hydrost. press. grad. ln_sppt_dynhpg :', ln_sppt_dynhpg
1593	WRITE(numout,*) ' Switch for surface press. grad. ln_sppt_dynspg :', ln_sppt_dynspg
1594	WRITE(numout,*) ' Switch for bottom friction ln_sppt_dynbfr :', ln_sppt_dynbfr
1595	WRITE(numout,*) ' Switch for Asselin time-filter ln_sppt_dynatf :', ln_sppt_dynatf
1596	WRITE(numout,*)
1597	WRITE(numout,*) ' Perturbation on sea-ice:'
1598	WRITE(numout,*) ' Switch for sea-ice diffusivity ln_sppt_icehdf :', ln_sppt_icehdf
1599	WRITE(numout,*) ' Switch for sea-ice lateral accretion ln_sppt_icelat :', ln_sppt_icelat
1600	WRITE(numout,*) ' Switch for sea-ice vertical thermodyn. ln_sppt_icezdf :', ln_sppt_icezdf
1601	WRITE(numout,*)
1602	WRITE(numout,*) ' Bound for gaussian random numbers rn_sppt_bound :', rn_sppt_bound
1603	WRITE(numout,*) ' Bound Step (0: nobound; 1: Gaussian; 2: Pert) nn_sppt_step_bound :', nn_sppt_step_bound
1604	WRITE(numout,*) ' Definition of Tau (1: days; 2: timesteps) nn_deftau :', nn_deftau
1605	WRITE(numout,*) ' Smoothing of perturbation close to coast ln_distcoast :', ln_distcoast
1606	WRITE(numout,*) ' Spatial scale of the smoothing near coasts (m) rn_distcoast :', rn_distcoast
1607	WRITE(numout,*) ' Type of vertical weight: nn_vwei :', nn_vwei
1608	WRITE(numout,*) ' 0 : No weight '
1609	WRITE(numout,*) ' 1 : First/Last level smoothing '
1610	WRITE(numout,*) ' 2 : Top/Bottom smoothing'
1611	WRITE(numout,*) ' 3 : Bottom smoothing'
1612	WRITE(numout,*)
1613	WRITE(numout,) ' ** SPP schemes (0: off ; 1: normal; 2:lognormal, mean as unpert.; 3: lognormal, median as unpert.'
1614	WRITE(numout,*) ' (the meaning of standard dev. depends on the distribution and parametr.)'
1615	WRITE(numout,*)
1616	WRITE(numout,*) ' Number of passes for spatial filter (AR1 field) spp_filter_pass:', spp_filter_pass
1617	WRITE(numout,*) ' Standard deviation of random generator (AR1 field) rn_spp_stdev :', rn_spp_stdev
1618	WRITE(numout,*) ' Decorr. time scale (AR1 field) rn_spp_tau :', rn_spp_tau
1619	WRITE(numout,*)
1620	WRITE(numout,*) ' SPP for bottom friction coeff nn_spp_bfr :', nn_spp_bfr
1621	WRITE(numout,*) ' STDEV rn_bfr_sd :', rn_bfr_sd
1622	WRITE(numout,*) ' SPP for SST relaxation coeff nn_spp_dqdt :', nn_spp_dqdt
1623	WRITE(numout,*) ' STDEV rn_dqdt_sd :', rn_dqdt_sd
1624	WRITE(numout,*) ' SPP for SSS relaxation coeff nn_spp_dedt :', nn_spp_dedt
1625	WRITE(numout,*) ' STDEV rn_dedt_sd :', rn_dedt_sd
1626	WRITE(numout,*) ' SPP for vertical tra mixing coeff (only TKE, GLS) nn_spp_avt :', nn_spp_avt
1627	WRITE(numout,*) ' STDEV rn_avt_sd :', rn_avt_sd
1628	WRITE(numout,*) ' SPP for vertical dyn mixing coeff (only TKE, GLS) nn_spp_avm :', nn_spp_avm
1629	WRITE(numout,*) ' STDEV rn_avm_sd :', rn_avm_sd
1630	WRITE(numout,*) ' SPP for solar penetration scheme (only RGB) nn_spp_qsi0 :', nn_spp_qsi0
1631	WRITE(numout,*) ' STDEV rn_qsi0_sd :', rn_qsi0_sd
1632	WRITE(numout,*) ' SPP for horiz. diffusivity U nn_spp_ahtu :', nn_spp_ahtu
1633	WRITE(numout,*) ' STDEV rn_ahtu_sd :', rn_ahtu_sd
1634	WRITE(numout,*) ' SPP for horiz. diffusivity V nn_spp_ahtv :', nn_spp_ahtv
1635	WRITE(numout,*) ' STDEV rn_ahtv_sd :', rn_ahtv_sd
1636	WRITE(numout,*) ' SPP for horiz. diffusivity W nn_spp_ahtw :', nn_spp_ahtw
1637	WRITE(numout,*) ' STDEV rn_ahtw_sd :', rn_ahtw_sd
1638	WRITE(numout,*) ' SPP for horiz. diffusivity T nn_spp_ahtt :', nn_spp_ahtt
1639	WRITE(numout,*) ' STDEV rn_ahtt_sd :', rn_ahtt_sd
1640	WRITE(numout,*) ' SPP for horiz. viscosity (1/3) nn_spp_ahm1 :', nn_spp_ahm1
1641	WRITE(numout,*) ' STDEV rn_ahm1_sd :', rn_ahm1_sd
1642	WRITE(numout,*) ' SPP for horiz. viscosity (2/4) nn_spp_ahm2 :', nn_spp_ahm2
1643	WRITE(numout,*) ' STDEV rn_ahm2_sd :', rn_ahm2_sd
1644	WRITE(numout,*) ' SPP for relative wind factor nn_spp_relw :', nn_spp_relw
1645	WRITE(numout,*) ' (use 4, 5, 6 for nn_spp_relw to have options 1, 2, 3 with limits bounded to [0,1]'
1646	WRITE(numout,*) ' STDEV rn_relw_sd :', rn_relw_sd
1647	WRITE(numout,*) ' SPP for mixing close to river mouth nn_spp_arnf :', nn_spp_arnf
1648	WRITE(numout,*) ' STDEV rn_arnf_sd :', rn_arnf_sd
1649	WRITE(numout,*) ' SPP for geothermal heating nn_spp_geot :', nn_spp_geot
1650	WRITE(numout,*) ' STDEV rn_geot_sd :', rn_geot_sd
1651	WRITE(numout,*) ' SPP for enhanced vertical diffusion nn_spp_aevd :', nn_spp_aevd
1652	WRITE(numout,*) ' STDEV rn_aevd_sd :', rn_aevd_sd
1653	WRITE(numout,*) ' SPP for TKE rn_lc Langmuir cell coefficient nn_spp_tkelc :', nn_spp_tkelc
1654	WRITE(numout,*) ' STDEV rn_tkelc_sd :', rn_tkelc_sd
1655	WRITE(numout,*) ' SPP for TKE rn_ediff Eddy diff. coefficient nn_spp_tkedf :', nn_spp_tkedf
1656	WRITE(numout,*) ' STDEV rn_tkedf_sd :', rn_tkedf_sd
1657	WRITE(numout,*) ' SPP for TKE rn_ediss Kolmogoroff dissipation coeff. nn_spp_tkeds :', nn_spp_tkeds
1658	WRITE(numout,*) ' STDEV rn_tkeds_sd :', rn_tkeds_sd
1659	WRITE(numout,*) ' SPP for TKE rn_ebb Surface input of tke nn_spp_tkebb :', nn_spp_tkebb
1660	WRITE(numout,*) ' STDEV rn_tkebb_sd :', rn_tkebb_sd
1661	WRITE(numout,*) ' SPP for TKE rn_efr Fraction of srf TKE below ML nn_spp_tkefr :', nn_spp_tkefr
1662	WRITE(numout,*) ' STDEV rn_tkefr_sd :', rn_tkefr_sd
1663	WRITE(numout,*) ' SPP for BBL U diffusivity nn_spp_ahubbl:', nn_spp_ahubbl
1664	WRITE(numout,*) ' STDEV rn_ahubbl_sd :', rn_ahubbl_sd
1665	WRITE(numout,*) ' SPP for BBL V diffusivity nn_spp_ahvbbl:', nn_spp_ahvbbl
1666	WRITE(numout,*) ' STDEV rn_ahvbbl_sd :', rn_ahvbbl_sd
1667	WRITE(numout,*)
1668	WRITE(numout,) ' ** SPP schemes for sea-ice '
1669	WRITE(numout,*) ' Albedo nn_spp_icealb:', nn_spp_icealb
1670	WRITE(numout,*) ' St. dev. for ice albedo rn_icealb_sd :', rn_icealb_sd
1671	WRITE(numout,*) ' Ice Strength nn_spp_icestr:', nn_spp_icestr
1672	WRITE(numout,*) ' St. dev. for ice strength rn_icestr_sd :', rn_icestr_sd
1673	WRITE(numout,*)
1674	WRITE(numout,*) ' SKEB Perturbation scheme '
1675	WRITE(numout,*) ' SKEB switch ln_skeb :', ln_skeb
1676	WRITE(numout,*) ' SKEB ratio of backscattered energy rn_skeb :', rn_skeb
1677	WRITE(numout,*) ' Frequency update for dissipation mask nn_dcom_freq :', nn_dcom_freq
1678	WRITE(numout,*) ' Numerical dissipation factor (resolut. dependent) rn_kh :', rn_kh
1679	WRITE(numout,*) ' Number of passes for spatial filter (AR1 field) skeb_filter_pass:', skeb_filter_pass
1680	WRITE(numout,*) ' Standard deviation of random generator (AR1 field) rn_skeb_stdev:', rn_skeb_stdev
1681	WRITE(numout,*) ' Decorr. time scale (AR1 field) rn_skeb_tau :', rn_skeb_tau
1682	WRITE(numout,*) ' Option of convection energy dissipation nn_dconv :', nn_dconv
1683	WRITE(numout,*) ' Convection dissipation factor (resolut. dependent) rn_kc :', rn_kc
1684	WRITE(numout,*) ' Multiplier for numerical dissipation rn_beta_num :', rn_beta_num
1685	WRITE(numout,*) ' Multiplier for convection dissipation rn_beta_con :', rn_beta_con
1686
1687	WRITE(numout,*)
1688	ENDIF
1689
1690	IF( .NOT. ln_stopack ) THEN
1691	IF(lwp) WRITE(numout,*) ' STOPACK is switched off'
1692	RETURN
1693	ENDIF
1694
1695	IF( MOD( nitend - nit000 + 1, nn_rndm_freq) /= 0 .OR. &
1696	( MOD( nstock, nn_rndm_freq) /= 0 .AND. nstock .GT. 0 ) ) THEN
1697	WRITE(ctmp1,*) 'experiment length (', nitend - nit000 + 1, ') or nstock (', nstock, &
1698	& ' is NOT a multiple of nn_rndm_freq (', nn_rndm_freq, ')'
1699	CALL ctl_stop( ctmp1, 'Impossible to properly setup STOPACK restart' )
1700	ENDIF
1701
1702	nn_sppt_tra = COUNT( (/ln_sppt_traxad, ln_sppt_trayad, ln_sppt_trazad, &
1703	ln_sppt_trasad, ln_sppt_traldf, ln_sppt_trazdf, &
1704	ln_sppt_trazdfp, ln_sppt_traevd, ln_sppt_trabbc, &
1705	ln_sppt_trabbl, ln_sppt_tranpc, ln_sppt_tradmp, &
1706	ln_sppt_traqsr, ln_sppt_transr, ln_sppt_traatf/) )
1707	nn_sppt_dyn = COUNT( (/ln_sppt_dynhpg, ln_sppt_dynspg, ln_sppt_dynkeg, &
1708	ln_sppt_dynrvo, ln_sppt_dynpvo, ln_sppt_dynzad, &
1709	ln_sppt_dynldf, ln_sppt_dynzdf, ln_sppt_dynbfr, &
1710	ln_sppt_dynatf, ln_sppt_dyntau/) )
1711	nn_sppt_ice = COUNT( (/ln_sppt_icehdf, ln_sppt_icelat, ln_sppt_icezdf/) )
1712
1713	ln_sppt_tra = ( nn_sppt_tra > 0 )
1714	ln_sppt_dyn = ( nn_sppt_dyn > 0 )
1715	ln_sppt_ice = ( nn_sppt_ice > 0 )
1716
1717	nn_spp = nn_spp_bfr+nn_spp_dqdt+nn_spp_dedt+nn_spp_avt+nn_spp_avm+nn_spp_qsi0+&
1718	& nn_spp_ahtu+nn_spp_ahtv+nn_spp_ahtw+nn_spp_ahtt+nn_spp_relw+nn_spp_arnf+nn_spp_geot+nn_spp_aevd+&
1719	& nn_spp_tkelc+nn_spp_tkedf+nn_spp_tkeds+nn_spp_tkebb+nn_spp_tkefr+nn_spp_ahubbl+nn_spp_ahvbbl+&
1720	& nn_spp_ahm1+nn_spp_ahm2
1721
1722	#ifdef NEMO_V34
1723
1724	#ifndef key_trdtra
1725	IF( ln_sppt_tra ) &
1726	& CALL ctl_stop( ' SPPT on tracers on .AND. .NOT. key_trdtra ', &
1727	& ' SPPT cannot work without tracer tendency computation')
1728	#endif
1729	#ifndef key_trddyn
1730	IF( ln_sppt_dyn ) &
1731	& CALL ctl_stop( ' SPPT on momentum on .AND. .NOT. key_trddyn ', &
1732	& ' SPPT cannot work without dynamic tendency computation')
1733	#endif
1734
1735	#else
1736	IF( ln_sppt_tra .AND. .NOT. ln_tra_trd ) &
1737	& CALL ctl_stop( ' SPPT on tracers on .AND. .NOT. ln_tra_trd ', &
1738	& ' SPPT cannot work without tracer tendency computation')
1739
1740	IF( ln_sppt_dyn .AND. .NOT. ln_dyn_trd ) &
1741	& CALL ctl_stop( ' SPPT on momentum on .AND. .NOT. ln_dyn_trd ', &
1742	& ' SPPT cannot work without dynamic tendency computation')
1743	#endif
1744
1745	IF( ln_sppt_glocon .AND. lk_vvl ) &
1746	& CALL ctl_stop( ' ln_sppt_glocal .AND. lk_vvl ', &
1747	& ' SPPT conservation not coded yet for VVL')
1748
1749	IF( nn_deftau .NE. 1 .AND. nn_deftau .NE. 2 ) &
1750	& CALL ctl_stop( ' nn_deftau must be 1 or 2 ')
1751
1752	nn_spp_aht = nn_spp_ahtu + nn_spp_ahtv + nn_spp_ahtw + nn_spp_ahtt
1753	IF(nn_spp_aht .GT. 0) THEN
1754	#if defined key_traldf_c3d
1755	IF(lwp) WRITE(numout,*) 'SPP : diffusivity perturbation with 3D coefficients'
1756	#elif defined key_traldf_c2d
1757	IF(lwp) WRITE(numout,*) 'SPP : diffusivity perturbation with 2D coefficients'
1758	#else
1759	CALL ctl_stop( 'SPP : diffusivity perturbation requires key_traldf_c3d or key_traldf_c2d')
1760	#endif
1761	ENDIF
1762
1763	IF(nn_spp_ahm1+nn_spp_ahm2 .GT. 0) THEN
1764	#if defined key_dynldf_c3d
1765	IF(lwp) WRITE(numout,*) 'SPP : viscosity perturbation with 3D coefficients'
1766	#elif defined key_dynldf_c2d
1767	IF(lwp) WRITE(numout,*) 'SPP : viscosity perturbation with 2D coefficients'
1768	#else
1769	CALL ctl_stop( 'SPP : viscosity perturbation requires key_dynldf_c3d or key_dynldf_c2d')
1770	#endif
1771	ENDIF
1772
1773	ln_skeb_own_gauss = .FALSE.
1774	IF( ln_skeb ) THEN
1775	IF( skeb_filter_pass > 0 .AND. skeb_filter_pass .NE. sppt_filter_pass ) ln_skeb_own_gauss = .TRUE.
1776	IF( rn_skeb_tau > 0._wp .AND. rn_skeb_tau .NE. rn_sppt_tau ) ln_skeb_own_gauss = .TRUE.
1777	ENDIF
1778
1779	ln_spp_own_gauss = .FALSE.
1780	IF( nn_spp > 0 ) THEN
1781	IF( spp_filter_pass > 0 .AND. spp_filter_pass .NE. sppt_filter_pass ) ln_spp_own_gauss = .TRUE.
1782	IF( rn_spp_tau > 0._wp .AND. rn_spp_tau .NE. rn_sppt_tau ) ln_spp_own_gauss = .TRUE.
1783	ENDIF
1784
1785	IF(lwp) THEN
1786	WRITE(numout,*) ' SPPT on tracers : ',ln_sppt_tra
1787	WRITE(numout,*) ' SPPT on dynamics : ',ln_sppt_dyn
1788	WRITE(numout,*) ' SPPT on sea-ice : ',ln_sppt_ice
1789	WRITE(numout,*) ' SPP perturbed params: ',nn_spp
1790	WRITE(numout,*) ' SKEB scheme : ',ln_skeb
1791	WRITE(numout,*) ' SPP with own scales : ',ln_spp_own_gauss
1792	WRITE(numout,*) ' SKEB with own scales: ',ln_skeb_own_gauss
1793	ENDIF
1794
1795	IF( sppt_tra_alloc() /= 0) CALL ctl_stop( 'STOP', 'sppt_alloc: unable to allocate arrays' )
1796
1797	spptt = 0._wp ; sppts = 0._wp
1798	spptu = 0._wp ; spptv = 0._wp
1799
1800	! Find filter attenuation factor
1801
1802	flt_fac = sppt_fltfac( sppt_filter_pass )
1803	rdt_sppt = nn_rndm_freq * rn_rdt
1804
1805	IF( ln_sppt_taumap ) THEN
1806	CALL iom_open ( 'sppt_tau_map', inum )
1807	CALL iom_get ( inum, jpdom_data, 'tau', sppt_tau )
1808	CALL iom_close( inum )
1809	ELSE
1810	sppt_tau(:,:) = rn_sppt_tau
1811	ENDIF
1812
1813	IF( nn_deftau .EQ. 1 ) THEN
1814	! Decorrelation time-scale defined in days
1815	sppt_tau(:,:) = exp( -rdt_sppt / (sppt_tau(:,:)*86400._wp) )
1816	ELSE
1817	! Decorrelation time-scale defined in timesteps
1818	sppt_tau(:,:) = exp( -REAL( nn_rndm_freq, wp) / sppt_tau(:,:) )
1819	ENDIF
1820
1821	IF( ln_distcoast ) THEN
1822	CALL iom_open('dist.coast', inum )
1823	CALL iom_get(inum,jpdom_data,'Tcoast',zdc)
1824	CALL iom_close( inum )
1825	zdc = 1._wp - exp(-zdc/rn_distcoast)
1826	CALL lbc_lnk( zdc , 'T', 1._wp)
1827	ELSE
1828	zdc = 1._wp
1829	ENDIF
1830
1831	! Initialize
1832	sppt_a = sppt_tau
1833	sppt_b = sqrt ( 1._wp - sppt_tau*sppt_tau )
1834
1835	IF(lwp) THEN
1836	WRITE(numout,*)
1837	WRITE(numout,) ' *** SPPT SCHEME'
1838	WRITE(numout,*) ' Definit. time-scale : ',nn_deftau
1839	WRITE(numout,*) ' Decorr. time-scale : ',rn_sppt_tau
1840	WRITE(numout,) ' Mean AR1 a term : ',SUM(sppt_a)/REAL(jpijpj)
1841	WRITE(numout,) ' Mean AR1 b term : ',SUM(sppt_b)/REAL(jpijpj)
1842	WRITE(numout,*)
1843	ENDIF
1844
1845	gauss_b = 0._wp
1846	! Weigths
1847	gauss_w(:) = 1.0_wp
1848	IF( nn_vwei .eq. 1 ) THEN
1849	gauss_w(1) = 0.0_wp
1850	gauss_w(2) = 0.5_wp
1851	gauss_w(jpk) = 0.0_wp
1852	gauss_w(jpkm1)= 0.5_wp
1853	ENDIF
1854
1855	IF( ln_spp_own_gauss ) THEN
1856	IF( spp_alloc() /= 0) CALL ctl_stop( 'STOP', 'spp_alloc: unable to allocate arrays' )
1857	flt_fac_p = sppt_fltfac( spp_filter_pass )
1858	spp_tau (:, :) = rn_spp_tau
1859	IF( nn_deftau .EQ. 1 ) THEN
1860	spp_tau(:,:) = spp_tau(:,:) * 86400._wp
1861	spp_tau(:,:) = exp( -rdt_sppt / (spp_tau) )
1862	ELSE
1863	spp_tau(:,:) = exp( -1._wp / spp_tau )
1864	ENDIF
1865	spp_a = spp_tau
1866	spp_b = sqrt ( 1._wp - spp_tau*spp_tau )
1867	gauss_bp = 0._wp
1868	IF(lwp) THEN
1869	WRITE(numout,*)
1870	WRITE(numout,) ' *** SPP SCHEME'
1871	WRITE(numout,*) ' Definit. time-scale : ',nn_deftau
1872	WRITE(numout,*) ' Decorr. time-scale : ',rn_spp_tau
1873	WRITE(numout,) ' Mean AR1 a term : ',SUM(spp_a)/REAL(jpijpj)
1874	WRITE(numout,) ' Mean AR1 b term : ',SUM(spp_b)/REAL(jpijpj)
1875	WRITE(numout,*)
1876	ENDIF
1877	ENDIF
1878
1879	IF( ln_skeb_own_gauss ) THEN
1880	IF( skeb_alloc() /= 0) CALL ctl_stop( 'STOP', 'skeb_alloc: unable to allocate arrays' )
1881	flt_fac_k = sppt_fltfac( skeb_filter_pass )
1882	skeb_tau (:, :) = rn_skeb_tau
1883	IF( nn_deftau .EQ. 1 ) THEN
1884	skeb_tau(:,:) = skeb_tau(:,:) * 86400._wp
1885	skeb_tau(:,:) = exp( -rdt_sppt / (skeb_tau) )
1886	ELSE
1887	skeb_tau(:,:) = exp( -1._wp / skeb_tau )
1888	ENDIF
1889	skeb_a = skeb_tau
1890	skeb_b = sqrt ( 1._wp - skeb_tau*skeb_tau )
1891	gauss_bk = 0._wp
1892	IF(lwp) THEN
1893	WRITE(numout,*)
1894	WRITE(numout,) ' *** SEB SCHEME'
1895	WRITE(numout,*) ' Definit. time-scale : ',nn_deftau
1896	WRITE(numout,*) ' Decorr. time-scale : ',rn_skeb_tau
1897	WRITE(numout,) ' Mean AR1 a term : ',SUM(skeb_a)/REAL(jpijpj)
1898	WRITE(numout,) ' Mean AR1 b term : ',SUM(skeb_b)/REAL(jpijpj)
1899	WRITE(numout,*)
1900	ENDIF
1901	ENDIF
1902
1903	IF( ln_skeb_own_gauss .OR. ln_spp_own_gauss ) &
1904	& ALLOCATE ( g2d_save(jpi,jpj) )
1905
1906	CALL stopack_rst( nit000, 'READ' )
1907
1908	IF(lwp .and. ln_stopack_diags) &
1909	CALL ctl_opn(numdiag, 'stopack.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
1910
1911	END SUBROUTINE stopack_init
1912	!
1913	SUBROUTINE stopack_rst( kt, cdrw )
1914	!!----------------------------------------------------------------------
1915	!! * ROUTINE stopack_rst *
1916	!!
1917	!! ** Purpose : Read/write from/to restarsts seed and perturbation field
1918	!!----------------------------------------------------------------------
1919	INTEGER , INTENT(in) :: kt ! ocean time-step
1920	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
1921	INTEGER :: id1, jseed
1922	CHARACTER(LEN=10) :: clseed='spsd0_0000'
1923	INTEGER(KIND=8) :: ziseed(4) ! RNG seeds in integer type
1924	INTEGER(KIND=8) :: ivals(8)
1925	REAL(wp) :: zrseed4(4) ! RNG seeds in integer type
1926	REAL(wp) :: zrseed2d(jpi,jpj)
1927	!
1928	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
1929	! ! ---------------
1930	IF(lwp) WRITE(numout,) ' stopack-rst *'
1931	IF( ln_rstart ) THEN !* Read the restart file
1932	id1 = iom_varid( numror, 'sppt_b', ldstop = .FALSE. )
1933	!
1934	IF( id1 > 0 ) THEN
1935	CALL iom_get( numror, jpdom_autoglo, 'sppt_b', gauss_b )
1936	DO jseed = 1 , 4
1937	WRITE(clseed(5:5) ,'(i1.1)') jseed
1938	CALL iom_get( numror, jpdom_autoglo, clseed(1:5) , zrseed2d )
1939	zrseed4(jseed) = zrseed2d(2,2)
1940	ziseed(jseed) = TRANSFER( zrseed4(jseed) , ziseed(jseed) )
1941	END DO
1942	IF (lwp) THEN
1943	WRITE(numout,*) 'Read ziseed4 = ',zrseed4(:)
1944	WRITE(numout,*) 'Read ziseed = ',ziseed(:)
1945	ENDIF
1946	CALL kiss_seed( ziseed(1) , ziseed(2) , ziseed(3) , ziseed(4) )
1947	ELSE
1948	IF( lwp ) WRITE(numout,*) ' ===>>>> : previous run without sppt_b'
1949	IF( ln_stopack_restart ) CALL ctl_stop( ' ln_stopack_restart TRUE :',&
1950	& ' variable not found in restart file ')
1951	IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of sppt_b to 0.'
1952	gauss_b = 0._wp
1953	IF ( .NOT. ln_stopack_repr ) THEN
1954	CALL date_and_time(VALUES=ivals)
1955	DO jseed=1,4
1956	nn_stopack_seed(jseed) = nn_stopack_seed(jseed) + INT(ivals(4+jseed))
1957	ENDDO
1958	ENDIF
1959	DO jseed=1,4
1960	ziseed(jseed) = TRANSFER(nn_stopack_seed(jseed)+narea, ziseed(jseed))
1961	ENDDO
1962	IF(lwp) WRITE(numout,*) ' ===>>>> Seed, nn_stopack_seed+narea',nn_stopack_seed+narea
1963	IF(lwp) WRITE(numout,*) ' ===>>>> Seed, ziseed',ziseed
1964	CALL kiss_seed( ziseed(1) , ziseed(2) , ziseed(3) , ziseed(4) )
1965	ENDIF
1966	IF( ln_spp_own_gauss ) THEN
1967	id1 = iom_varid( numror, 'spp_b', ldstop = .FALSE. )
1968	IF( id1 > 0 ) THEN
1969	CALL iom_get( numror, jpdom_autoglo, 'spp_b', gauss_bp )
1970	ELSE
1971	IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of spp_b to 0.'
1972	gauss_bp = 0._wp
1973	ENDIF
1974	ENDIF
1975	IF( ln_skeb_own_gauss ) THEN
1976	id1 = iom_varid( numror, 'skeb_b', ldstop = .FALSE. )
1977	IF( id1 > 0 ) THEN
1978	CALL iom_get( numror, jpdom_autoglo, 'skeb_b', gauss_bk )
1979	ELSE
1980	IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of skeb_b to 0.'
1981	gauss_bk = 0._wp
1982	ENDIF
1983	ENDIF
1984	ELSE
1985	gauss_b = 0._wp
1986	gauss_bp = 0._wp
1987	gauss_bk = 0._wp
1988	IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of STOPACK to 0.'
1989	IF ( .NOT. ln_stopack_repr ) THEN
1990	CALL date_and_time(VALUES=ivals)
1991	DO jseed=1,4
1992	nn_stopack_seed(jseed) = nn_stopack_seed(jseed) + INT(ivals(4+jseed))
1993	ENDDO
1994	ENDIF
1995	DO jseed=1,4
1996	ziseed(jseed) = TRANSFER(nn_stopack_seed(jseed)+narea, ziseed(jseed))
1997	ENDDO
1998	IF(lwp) WRITE(numout,*) ' ===>>>> Seed, nn_stopack_seed+narea',nn_stopack_seed+narea
1999	IF(lwp) WRITE(numout,*) ' ===>>>> Seed, ziseed',ziseed
2000	CALL kiss_seed( ziseed(1) , ziseed(2) , ziseed(3) , ziseed(4) )
2001	ENDIF
2002	!
2003	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
2004	! ! -------------------
2005	IF(lwp) WRITE(numout,*) '---- stopack-rst ----'
2006	CALL iom_rstput( kt, nitrst, numrow, 'sppt_b', gauss_b )
2007	IF(ln_spp_own_gauss) CALL iom_rstput( kt, nitrst, numrow, 'spp_b', gauss_bp )
2008	IF(ln_skeb_own_gauss) CALL iom_rstput( kt, nitrst, numrow, 'skeb_b', gauss_bk )
2009	CALL kiss_state( ziseed(1) , ziseed(2) , ziseed(3) , ziseed(4) )
2010	DO jseed=1,4
2011	zrseed4(jseed) = TRANSFER(ziseed(jseed), zrseed4(jseed))
2012	ENDDO
2013	IF (lwp) THEN
2014	WRITE(numout,*) 'Write ziseed4 = ',zrseed4(:)
2015	WRITE(numout,*) 'Write ziseed = ',ziseed(:)
2016	ENDIF
2017	DO jseed = 1 , 4
2018	WRITE(clseed(5:5) ,'(i1.1)') jseed
2019	zrseed2d(:,:) = zrseed4(jseed)
2020	CALL iom_rstput( kt, nitrst, numrow, clseed(1:5) , zrseed2d )
2021	END DO
2022	!
2023	ENDIF
2024	!
2025	END SUBROUTINE stopack_rst
2026	!
2027	INTEGER FUNCTION sppt_tra_alloc()
2028	!!---------------------------------------------------------------------
2029	!! * FUNCTION sppt_tra_alloc *
2030	!!---------------------------------------------------------------------
2031	!
2032	ALLOCATE( spptt(jpi,jpj,jpk) , sppts(jpi,jpj,jpk) , gauss_n(jpi,jpj,jpk) ,&
2033	gauss_nu(jpi,jpj,jpk) , gauss_nv(jpi,jpj,jpk) , &
2034	spptu(jpi,jpj,jpk) , spptv(jpi,jpj,jpk) , gauss_n_2d(jpi,jpj) ,&
2035	gauss_b (jpi,jpj), sppt_tau(jpi,jpj), sppt_a(jpi,jpj), sppt_b(jpi,jpj), gauss_w(jpk),&
2036	zdc(jpi,jpj,jpk), STAT= sppt_tra_alloc )
2037	!
2038	IF( lk_mpp ) CALL mpp_sum ( sppt_tra_alloc )
2039	IF( sppt_tra_alloc /= 0) CALL ctl_warn('sppt_tra_alloc: failed to allocate arrays')
2040
2041	IF(nn_spp_bfr>0) THEN
2042	ALLOCATE(coeff_bfr(jpi,jpj), STAT= sppt_tra_alloc )
2043	IF( lk_mpp ) CALL mpp_sum ( sppt_tra_alloc )
2044	IF( sppt_tra_alloc /= 0) CALL ctl_warn('sppt_tra_alloc: failed to allocate arrays')
2045	coeff_bfr = 0._wp
2046	ENDIF
2047	gauss_b = 0._wp
2048	gauss_n = 0._wp
2049	gauss_nu= 0._wp
2050	gauss_nv= 0._wp
2051	gauss_n_2d = 0._wp
2052
2053	END FUNCTION sppt_tra_alloc
2054
2055	INTEGER FUNCTION skeb_alloc()
2056	!!---------------------------------------------------------------------
2057	!! * FUNCTION skeb_alloc *
2058	!!---------------------------------------------------------------------
2059	!
2060	ALLOCATE( gauss_n_2d_k(jpi,jpj) , gauss_bk (jpi,jpj),skeb_tau(jpi,jpj),&
2061	skeb_a(jpi,jpj), skeb_b(jpi,jpj), STAT= skeb_alloc )
2062	!
2063	IF( lk_mpp ) CALL mpp_sum ( skeb_alloc )
2064	IF( skeb_alloc /= 0) CALL ctl_warn('sppt_tra_alloc: failed to allocate arrays')
2065	gauss_n_2d_k = 0._wp
2066	gauss_bk = 0._wp
2067
2068	END FUNCTION skeb_alloc
2069
2070	INTEGER FUNCTION spp_alloc()
2071	!!---------------------------------------------------------------------
2072	!! * FUNCTION spp_alloc *
2073	!!---------------------------------------------------------------------
2074	!
2075	ALLOCATE( gauss_n_2d_p(jpi,jpj), gauss_bp (jpi,jpj),spp_tau(jpi,jpj), &
2076	spp_a(jpi,jpj), spp_b(jpi,jpj), STAT= spp_alloc )
2077	!
2078	IF( lk_mpp ) CALL mpp_sum ( spp_alloc )
2079	IF( spp_alloc /= 0) CALL ctl_warn('sppt_tra_alloc: failed to allocate arrays')
2080	gauss_n_2d_p = 0._wp
2081	gauss_bp = 0._wp
2082
2083	END FUNCTION spp_alloc
2084
2085	SUBROUTINE sppt_flt( psto )
2086	!!----------------------------------------------------------------------
2087	!! * ROUTINE sppt_flt *
2088	!!
2089	!! ** Purpose : apply horizontal Laplacian filter to input array
2090	!! Adapted from STO package
2091	!!----------------------------------------------------------------------
2092	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: psto
2093	REAL(wp), DIMENSION(jpi,jpj) :: pstm
2094	!!
2095	INTEGER :: ji, jj
2096
2097	pstm = psto
2098	DO jj = 2, jpj-1
2099	DO ji = 2, jpi-1
2100	psto(ji,jj) = 0.5_wp * pstm(ji,jj) + 0.125_wp * &
2101	& ( pstm(ji-1,jj) + pstm(ji+1,jj) + &
2102	& pstm(ji,jj-1) + pstm(ji,jj+1) )
2103	END DO
2104	END DO
2105	CALL lbc_lnk( psto , 'T', 1._wp )
2106
2107	END SUBROUTINE sppt_flt
2108
2109	SUBROUTINE sppt_rand2d( psto )
2110	!!----------------------------------------------------------------------
2111	!! * ROUTINE sppt_rand2d *
2112	!!
2113	!! ** Purpose : fill input array with white Gaussian noise
2114	!! Adapted from STO package
2115	!!----------------------------------------------------------------------
2116	REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: psto
2117	!!
2118	INTEGER :: ji, jj
2119	REAL(KIND=8) :: gran ! Gaussian random number (forced KIND=8 as in kiss_gaussian)
2120
2121	DO jj = 1, jpj
2122	DO ji = 1, jpi
2123	CALL kiss_gaussian( gran )
2124	psto(ji,jj) = gran*rn_sppt_stdev
2125	END DO
2126	END DO
2127
2128	END SUBROUTINE sppt_rand2d
2129
2130	FUNCTION sppt_fltfac( kpasses )
2131	!!----------------------------------------------------------------------
2132	!! * FUNCTION sppt_fltfac *
2133	!!
2134	!! ** Purpose : compute factor to restore standard deviation
2135	!! as a function of the number of passes
2136	!! of the Laplacian filter
2137	!! From STO package
2138	!!----------------------------------------------------------------------
2139	INTEGER, INTENT(in) :: kpasses
2140	REAL(wp) :: sppt_fltfac
2141	!!
2142	INTEGER :: jpasses, ji, jj, jflti, jfltj
2143	INTEGER, DIMENSION(-1:1,-1:1) :: pflt0
2144	! 16-b reals to avoid overflows
2145	INTEGER, PARAMETER :: qp = SELECTED_REAL_KIND(33,4931)
2146	REAL(qp), DIMENSION(:,:), ALLOCATABLE :: pfltb
2147	REAL(qp), DIMENSION(:,:), ALLOCATABLE :: pflta
2148	REAL(qp) :: ratio
2149	REAL(qp) :: aux0, aux1
2150	pflt0(-1,-1) = 0 ; pflt0(-1,0) = 1 ; pflt0(-1,1) = 0
2151	pflt0( 0,-1) = 1 ; pflt0( 0,0) = 4 ; pflt0( 0,1) = 1
2152	pflt0( 1,-1) = 0 ; pflt0( 1,0) = 1 ; pflt0( 1,1) = 0
2153	ALLOCATE(pfltb(-kpasses-1:kpasses+1,-kpasses-1:kpasses+1))
2154	ALLOCATE(pflta(-kpasses-1:kpasses+1,-kpasses-1:kpasses+1))
2155	pfltb(:,:) = 0
2156	pfltb(0,0) = 1
2157	DO jpasses = 1, kpasses
2158	pflta(:,:) = 0._qp
2159	DO jflti= -1, 1
2160	DO jfltj= -1, 1
2161	DO ji= -kpasses, kpasses
2162	DO jj= -kpasses, kpasses
2163	pflta(ji,jj) = pflta(ji,jj) + pfltb(ji+jflti,jj+jfltj) * pflt0(jflti,jfltj)
2164	ENDDO
2165	ENDDO
2166	ENDDO
2167	ENDDO
2168	pfltb(:,:) = pflta(:,:)
2169	ENDDO
2170	ratio = SUM(pfltb(:,:))
2171	aux0 = SUM(pfltb(:,:)*pfltb(:,:))
2172	aux1 = ratio*ratio
2173	ratio = aux1 / aux0
2174	ratio = SQRT(ratio)
2175	DEALLOCATE(pfltb,pflta)
2176	sppt_fltfac = REAL(ratio, wp)
2177	END FUNCTION sppt_fltfac
2178
2179	SUBROUTINE skeb_comp( lt )
2180	!!----------------------------------------------------------------------
2181	!! * ROUTINE skeb_comp *
2182	!!
2183	!! ** Purpose : Computation of energy dissipation terms
2184	!! This is a wrapper to the enrgy terms computation
2185	!! routines
2186	!!----------------------------------------------------------------------
2187	INTEGER, INTENT(IN) :: lt
2188	CALL skeb_dnum ( lt )
2189	CALL skeb_dcon ( lt )
2190	END SUBROUTINE skeb_comp
2191
2192	SUBROUTINE bsfcomp( kt )
2193	!!----------------------------------------------------------------------
2194	!! * ROUTINE bsfcomp *
2195	!!
2196	!! ** Purpose : Online computation of barotropic streamfunction
2197	!! For diagnostics purposes
2198	!! This routine is not explicitly called anywhere
2199	!! and utilizes low-level MPI routines
2200	!!----------------------------------------------------------------------
2201	USE MPI
2202	!
2203	INTEGER, INTENT(IN) :: kt
2204	!
2205	REAL(wp) :: dtrpv(jpi,jpj)
2206	INTEGER :: jk,ji,jj
2207	!
2208	#if defined key_mpp_mpi
2209	INTEGER, DIMENSION(1) :: ish
2210	INTEGER, DIMENSION(2) :: ish2
2211	INTEGER :: ijp,tag,ierr,jp,isend,irecv
2212	REAL(wp), DIMENSION(jpj) :: zwrk ! mask flux array at V-point
2213	REAL(wp), DIMENSION(jpi) :: zwr2 ! mask flux array at V-point
2214	INTEGER :: istatus(MPI_STATUS_SIZE)
2215	#endif
2216	IF(kt .EQ. nit000 ) THEN
2217	#ifdef NEMO_V34
2218	ln_dpsiu = .FALSE.
2219	ln_dpsiv = .FALSE.
2220	#else
2221	IF (iom_use('bsft') ) ln_dpsiu = .TRUE.
2222	IF (iom_use('bsftv') ) ln_dpsiv = .TRUE.
2223	#endif
2224	IF( ln_dpsiv ) ALLOCATE ( dpsiv(jpi,jpj) )
2225	IF( ln_dpsiu ) ALLOCATE ( dpsiu(jpi,jpj) )
2226	IF( ln_dpsiv .OR. ln_dpsiu ) THEN
2227	jpri = nint ( REAL(nimpp) / REAL(jpi) ) + 1
2228	jprj = nint ( REAL(njmpp) / REAL(jpj) ) + 1
2229	ENDIF
2230	ENDIF
2231
2232	IF( ln_dpsiv ) THEN
2233	dtrpv = 0._wp
2234	DO jk = 1,jpkm1
2235	dtrpv = dtrpv + fse3v(:,:,jk) * e1v(:,:) * vn(:,:,jk)
2236	ENDDO
2237	dpsiv(1,:)=0._wp
2238	DO ji = 2,jpi
2239	dpsiv(ji,:) = dpsiv(ji-1,:) + dtrpv(ji,:)
2240	END DO
2241	ENDIF
2242
2243	IF( ln_dpsiu ) THEN
2244	dtrpv = 0._wp
2245	DO jk = 1,jpkm1
2246	dtrpv = dtrpv + fse3u(:,:,jk) * e2u(:,:) * un(:,:,jk)
2247	ENDDO
2248	dpsiu(:,1)=0._wp
2249	DO jj = 2,jpj
2250	dpsiu(:,jj) = dpsiu(:,jj-1) + dtrpv(:,jj)
2251	END DO
2252	ENDIF
2253
2254	#if defined key_mpp_mpi
2255	IF ( ln_dpsiv ) THEN
2256	DO jp=1,jpni-1
2257	IF( jpri == jp ) THEN ! SEND TO EAST
2258	zwrk(1:jpj) = dpsiv(jpi-1,:)
2259	tag=2000+narea
2260	CALL mpi_isend(zwrk, jpj, mpi_double_precision, noea, tag, mpi_comm_opa, isend,ierr)
2261	ELSEIF ( jpri == jp+1 ) THEN ! RECEIVE FROM WEST
2262	CALL mpi_irecv(zwrk, jpj, mpi_double_precision, nowe, mpi_any_tag, mpi_comm_opa, irecv,ierr)
2263	ENDIF
2264	IF(jpri == jp) CALL mpi_wait(isend, istatus, ierr)
2265	IF(jpri == jp+1 ) THEN
2266	CALL mpi_wait(irecv, istatus, ierr)
2267	DO ji=1,jpi
2268	dpsiv(ji,:) = dpsiv(ji,:) + zwrk(:)
2269	ENDDO
2270	ENDIF
2271	ENDDO
2272	ENDIF
2273
2274	IF ( ln_dpsiv ) THEN
2275	DO jp=1,jpnj-1
2276	IF( jprj == jp ) THEN ! SEND TO NORTH
2277	zwr2(1:jpi) = dpsiu(:,jpj-1)
2278	tag=3000+narea
2279	CALL mpi_isend(zwr2, jpi, mpi_double_precision, nono, tag, mpi_comm_opa, isend,ierr)
2280	ELSEIF ( jprj == jp+1 ) THEN ! RECEIVE FROM SOUTH
2281	CALL mpi_irecv(zwr2, jpi, mpi_double_precision, noso, mpi_any_tag, mpi_comm_opa, irecv,ierr)
2282	ENDIF
2283	IF(jprj == jp) CALL mpi_wait(isend, istatus, ierr)
2284	IF(jprj == jp+1 ) THEN
2285	CALL mpi_wait(irecv, istatus, ierr)
2286	DO ji=1,jpj
2287	dpsiu(:,ji) = dpsiu(:,ji) + zwr2(:)
2288	ENDDO
2289	ENDIF
2290	ENDDO
2291	ENDIF
2292	#endif
2293	IF (ln_dpsiu ) THEN
2294	CALL lbc_lnk(dpsiu,'T',1._wp)
2295	ENDIF
2296	IF (ln_dpsiv ) THEN
2297	CALL lbc_lnk(dpsiv,'T',1._wp)
2298	ENDIF
2299
2300	IF(kt == nitend ) THEN
2301	IF (ln_dpsiv ) DEALLOCATE ( dpsiv )
2302	IF (ln_dpsiu ) DEALLOCATE ( dpsiu )
2303	ENDIF
2304
2305	END SUBROUTINE
2306
2307	SUBROUTINE skeb_dnum ( mt )
2308	!!----------------------------------------------------------------------
2309	!! * ROUTINE skeb_dnum *
2310	!!
2311	!! ** Purpose : Computation of numerical energy dissipation
2312	!! For later use in the SKEB scheme
2313	!!----------------------------------------------------------------------
2314	INTEGER, INTENT(IN) :: mt
2315	REAL :: ds,dt,dtot,kh2
2316	INTEGER :: ji,jj,jk
2317
2318	IF ( mt .eq. nit000 ) THEN
2319	ALLOCATE ( dnum(jpi,jpj,jpk) )
2320	dnum (:,:,: ) = 0._wp
2321	ENDIF
2322
2323	kh2 = rn_kh * rn_kh
2324
2325	IF( mt .eq. nit000 .OR. MOD( mt - 1, nn_dcom_freq ) == 0 ) THEN
2326
2327	DO jk=1,jpkm1
2328	DO jj=2,jpj
2329	DO ji=2,jpi
2330	! Shear
2331	ds = (vn(ji,jj,jk)-vn(ji-1,jj,jk))vmask(ji,jj,jk)vmask(ji-1,jj,jk)/ e1v(ji,jj) + &
2332	(un(ji,jj,jk)-un(ji,jj-1,jk))umask(ji,jj,jk)umask(ji,jj-1,jk)/ e2u(ji,jj)
2333	! Tension
2334	dt = (vn(ji,jj,jk)-vn(ji-1,jj,jk))vmask(ji,jj,jk)vmask(ji-1,jj,jk)/ e2v(ji,jj) + &
2335	(un(ji,jj,jk)-un(ji,jj-1,jk))umask(ji,jj,jk)umask(ji,jj-1,jk)/ e1u(ji,jj)
2336
2337	dtot = sqrt( dsds + dtdt ) * tmask(ji,jj,jk)
2338	dnum(ji,jj,jk) = dtotdtotdtotkh2e1t(ji,jj)*e2t(ji,jj)
2339	ENDDO
2340	ENDDO
2341	ENDDO
2342
2343	CALL lbc_lnk(dnum,'T',1._wp)
2344
2345	ENDIF
2346
2347	END SUBROUTINE
2348
2349	SUBROUTINE skeb_dcon ( mt )
2350	!!----------------------------------------------------------------------
2351	!! * ROUTINE skeb_dcon *
2352	!!
2353	!! ** Purpose : Computation of convective energy dissipation
2354	!! For later use in the SKEB scheme
2355	!! Two formulation are implemented, based on buoyancy or
2356	!! convective mass flux formulations, respectively
2357	!!----------------------------------------------------------------------
2358	INTEGER, INTENT(IN) :: mt
2359	REAL(wp) :: zz, mf1,mf2,kc2
2360	INTEGER :: ji,jj,jk
2361
2362	IF ( mt .eq. nit000 ) THEN
2363	ALLOCATE ( dcon(jpi,jpj,jpk) )
2364	dcon (:,:,: ) = 0._wp
2365	ENDIF
2366
2367	kc2 = rn_kc * rn_kc
2368
2369	IF( mt .eq. nit000 .OR. MOD( mt - 1, nn_dcom_freq ) == 0 ) THEN
2370
2371	IF(nn_dconv .eq. 1 ) THEN
2372
2373	DO jk=2,jpkm1
2374	DO jj=1,jpj
2375	DO ji=1,jpi
2376
2377	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) &
2378	& / ( rau0 * fse3w(ji,jj,jk) )
2379
2380	dcon(ji,jj,jk) = kc2zze1t(ji,jj)e2t(ji,jj)rau0 / fse3w(ji,jj,jk)
2381
2382	ENDDO
2383	ENDDO
2384	ENDDO
2385
2386	ELSEIF (nn_dconv .eq. 2 ) THEN
2387
2388	DO jk=2,jpkm1
2389	DO jj=1,jpj
2390	DO ji=1,jpi
2391
2392	mf1 = wn(ji,jj,jk+1)e1t(ji,jj)e2t(ji,jj)
2393	mf2 = wn(ji,jj,jk-1)e1t(ji,jj)e2t(ji,jj)
2394	dcon(ji,jj,jk) = rn_kc * (mf1-mf2)(mf1-mf2) tmask(ji,jj,jk+1) / (fse3w(ji,jj,jk)rau0rau0)
2395
2396	ENDDO
2397	ENDDO
2398	ENDDO
2399
2400	ENDIF
2401
2402	CALL lbc_lnk(dcon,'T',1._wp)
2403
2404	ENDIF
2405
2406	END SUBROUTINE
2407
2408	SUBROUTINE skeb_apply ( mt)
2409	!!----------------------------------------------------------------------
2410	!! * ROUTINE skeb_apply *
2411	!!
2412	!! ** Purpose : Application of SKEB perturbation
2413	!! Convective and Numerical energy dissipation are
2414	!! multiplied by the beta terms
2415	!!----------------------------------------------------------------------
2416
2417	INTEGER, INTENT(IN) :: mt
2418	INTEGER :: ji,jj,jk
2419	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ustar,vstar,psi
2420	REAL(wp) :: mi,ma
2421
2422	IF( ln_stopack_diags ) THEN
2423
2424	psi = 0._wp
2425	IF(ln_skeb_own_gauss) THEN
2426	DO jk=1,jpkm1
2427	psi(:,:,jk) = rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk) ) * gauss_n_2d_k(:,:)
2428	ENDDO
2429	ELSE
2430	DO jk=1,jpkm1
2431	psi(:,:,jk) = rn_skeb_stdev * rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk) ) &
2432	& * gauss_n_2d(:,:) / rn_sppt_stdev
2433	ENDDO
2434	ENDIF
2435	ustar = 0._wp
2436	vstar = 0._wp
2437	DO jk=1,jpkm1
2438	DO jj=2,jpj
2439	DO ji=2,jpi
2440	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)
2441	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)
2442	ENDDO
2443	ENDDO
2444	ENDDO
2445	mi = MAXVAL(ABS(ustar))
2446	ma = MAXVAL(ABS(vstar))
2447	IF(lk_mpp) CALL mpp_max(mi)
2448	IF(lk_mpp) CALL mpp_max(ma)
2449	IF(lwp) WRITE(numdiag,9304) lkt,'DNUM ABS CURRENT' ,mi,ma
2450
2451	rn_mmax ( jk_skeb_dnum ) = MAX ( MAX( mi, ma), rn_mmax ( jk_skeb_dnum ) )
2452
2453	psi = 0._wp
2454	IF(ln_skeb_own_gauss) THEN
2455	DO jk=1,jpkm1
2456	psi(:,:,jk) = rn_skeb * sqrt( rn_beta_con * dcon(:,:,jk) ) * gauss_n_2d_k(:,:)
2457	ENDDO
2458	ELSE
2459	DO jk=1,jpkm1
2460	psi(:,:,jk) = rn_skeb_stdev * rn_skeb * sqrt( rn_beta_con * dcon(:,:,jk) ) &
2461	& * gauss_n_2d(:,:) / rn_sppt_stdev
2462	ENDDO
2463	ENDIF
2464	ustar = 0._wp
2465	vstar = 0._wp
2466	DO jk=1,jpkm1
2467	DO jj=2,jpj
2468	DO ji=2,jpi
2469	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)
2470	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)
2471	ENDDO
2472	ENDDO
2473	ENDDO
2474	mi = MAXVAL(ABS(ustar))
2475	ma = MAXVAL(ABS(vstar))
2476	IF(lk_mpp) CALL mpp_max(mi)
2477	IF(lk_mpp) CALL mpp_max(ma)
2478	IF(lwp) WRITE(numdiag,9304) lkt,'DCON ABS CURRENT' ,mi,ma
2479
2480	rn_mmax ( jk_skeb_dcon ) = MAX ( MAX( mi, ma), rn_mmax ( jk_skeb_dcon ) )
2481
2482	9304 FORMAT(' it :', i8, ' ', A16, ' Min: ',d10.3 , ' Max: ',d10.3)
2483
2484	ENDIF
2485
2486	psi = 0._wp
2487	IF(ln_skeb_own_gauss) THEN
2488	DO jk=1,jpkm1
2489	psi(:,:,jk) = rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk)+ rn_beta_con * dcon(:,:,jk) ) * gauss_n_2d_k(:,:)
2490	ENDDO
2491	ELSE
2492	DO jk=1,jpkm1
2493	psi(:,:,jk) = rn_skeb_stdev * rn_skeb * sqrt( rn_beta_num * dnum(:,:,jk)+ rn_beta_con * dcon(:,:,jk) ) &
2494	& * gauss_n_2d(:,:) / rn_sppt_stdev
2495	ENDDO
2496	ENDIF
2497
2498	ustar = 0._wp
2499	vstar = 0._wp
2500
2501	DO jk=1,jpkm1
2502	DO jj=2,jpj
2503	DO ji=2,jpi
2504	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)
2505	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)
2506	ENDDO
2507	ENDDO
2508	ENDDO
2509
2510	IF( ln_stopack_diags ) THEN
2511
2512	mi = MAXVAL(ABS(dnum))
2513	ma = MAXVAL(ABS(dcon))
2514	IF(lk_mpp) CALL mpp_max(mi)
2515	IF(lk_mpp) CALL mpp_max(ma)
2516	IF(lwp) WRITE(numdiag,9302) lkt,'DNUM AND DCON MX' ,mi,ma
2517
2518	mi = MINVAL(ustar)
2519	ma = MAXVAL(ustar)
2520	IF(lk_mpp) CALL mpp_min(mi)
2521	IF(lk_mpp) CALL mpp_max(ma)
2522	IF(lwp) WRITE(numdiag,9302) lkt,'SKEB U-CURRENT ' ,mi,ma
2523
2524	rn_mmax ( jk_skeb_tot ) = MAX ( MAX( abs(mi), ma), rn_mmax ( jk_skeb_tot ) )
2525
2526	mi = MINVAL(vstar)
2527	ma = MAXVAL(vstar)
2528	IF(lk_mpp) CALL mpp_min(mi)
2529	IF(lk_mpp) CALL mpp_max(ma)
2530	IF(lwp) WRITE(numdiag,9302) lkt,'SKEB V-CURRENT ' ,mi,ma
2531	9302 FORMAT(' it :', i8, ' ', A16, ' Min: ',d10.3 , ' Max: ',d10.3)
2532
2533	rn_mmax ( jk_skeb_tot ) = MAX ( MAX( abs(mi), ma), rn_mmax ( jk_skeb_tot ) )
2534
2535	ENDIF
2536
2537	CALL lbc_lnk(ustar,'U',-1._wp)
2538	CALL lbc_lnk(vstar,'V',-1._wp)
2539
2540	DO jk=1,jpkm1
2541	DO jj=1,jpj
2542	DO ji=1,jpi
2543	un(ji,jj,jk) = un(ji,jj,jk) + ustar(ji,jj,jk)
2544	vn(ji,jj,jk) = vn(ji,jj,jk) + vstar(ji,jj,jk)
2545	ENDDO
2546	ENDDO
2547	ENDDO
2548
2549	#ifdef key_iomput
2550	IF (iom_use('ustar_skeb') ) CALL iom_put( "ustar_skeb" , ustar)
2551	IF (iom_use('vstar_skeb') ) CALL iom_put( "vstar_skeb" , vstar)
2552	#endif
2553
2554	IF ( mt .eq. nitend ) THEN
2555	IF(ALLOCATED(dnum)) DEALLOCATE ( dnum )
2556	IF(ALLOCATED(dcon)) DEALLOCATE ( dcon )
2557	IF (ln_stopack_diags .AND. lwp) CALL stopack_report
2558	ENDIF
2559
2560	END SUBROUTINE
2561
2562	!!======================================================================
2563	!! Random number generator, used in NEMO stochastic parameterization
2564	!!
2565	!!=====================================================================
2566	!! History : 3.3 ! 2011-10 (J.-M. Brankart) Original code
2567	!!----------------------------------------------------------------------
2568
2569	!!----------------------------------------------------------------------
2570	!! The module is based on (and includes) the
2571	!! 64-bit KISS (Keep It Simple Stupid) random number generator
2572	!! distributed by George Marsaglia :
2573	!! http://groups.google.com/group/comp.lang.fortran/
2574	!! browse_thread/thread/a85bf5f2a97f5a55
2575	!!----------------------------------------------------------------------
2576
2577	!!----------------------------------------------------------------------
2578	!! kiss : 64-bit KISS random number generator (period ~ 2^250)
2579	!! kiss_seed : Define seeds for KISS random number generator
2580	!! kiss_state : Get current state of KISS random number generator
2581	!! kiss_save : Save current state of KISS (for future restart)
2582	!! kiss_load : Load the saved state of KISS
2583	!! kiss_reset : Reset the default seeds
2584	!! kiss_check : Check the KISS pseudo-random sequence
2585	!! kiss_uniform : Real random numbers with uniform distribution in [0,1]
2586	!! kiss_gaussian : Real random numbers with Gaussian distribution N(0,1)
2587	!! kiss_gamma : Real random numbers with Gamma distribution Gamma(k,1)
2588	!! kiss_sample : Select a random sample from a set of integers
2589	!!
2590	!! ---CURRENTLY NOT USED IN NEMO :
2591	!! kiss_save, kiss_load, kiss_check, kiss_gamma, kiss_sample
2592	!!----------------------------------------------------------------------
2593	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
2594	!! $Id: dynhpg.F90 2528 2010-12-27 17:33:53Z rblod $
2595	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
2596	!!----------------------------------------------------------------------
2597	FUNCTION kiss()
2598	!! --------------------------------------------------------------------
2599	!! * FUNCTION kiss *
2600	!!
2601	!! ** Purpose : 64-bit KISS random number generator
2602	!!
2603	!! ** Method : combine several random number generators:
2604	!! (1) Xorshift (XSH), period 2^64-1,
2605	!! (2) Multiply-with-carry (MWC), period (2^121+2^63-1)
2606	!! (3) Congruential generator (CNG), period 2^64.
2607	!!
2608	!! overall period:
2609	!! (2^250+2^192+2^64-2^186-2^129)/6
2610	!! ~= 2^(247.42) or 10^(74.48)
2611	!!
2612	!! set your own seeds with 'kiss_seed'
2613	! --------------------------------------------------------------------
2614	IMPLICIT NONE
2615	INTEGER(KIND=i8) :: kiss, t
2616
2617	t = ISHFT(x,58) + w
2618	IF (s(x).eq.s(t)) THEN
2619	w = ISHFT(x,-6) + s(x)
2620	ELSE
2621	w = ISHFT(x,-6) + 1 - s(x+t)
2622	ENDIF
2623	x = t + x
2624	y = m( m( m(y,13_8), -17_8 ), 43_8 )
2625	z = 6906969069_8 * z + 1234567_8
2626
2627	kiss = x + y + z
2628
2629	CONTAINS
2630
2631	FUNCTION s(k)
2632	INTEGER(KIND=i8) :: s, k
2633	s = ISHFT(k,-63)
2634	END FUNCTION s
2635
2636	FUNCTION m(k, n)
2637	INTEGER(KIND=i8) :: m, k, n
2638	m = IEOR(k, ISHFT(k, n) )
2639	END FUNCTION m
2640
2641	END FUNCTION kiss
2642
2643
2644	SUBROUTINE kiss_seed(ix, iy, iz, iw)
2645	!! --------------------------------------------------------------------
2646	!! * ROUTINE kiss_seed *
2647	!!
2648	!! ** Purpose : Define seeds for KISS random number generator
2649	!!
2650	!! --------------------------------------------------------------------
2651	IMPLICIT NONE
2652	INTEGER(KIND=i8) :: ix, iy, iz, iw
2653
2654	IF(lwp) WRITE(numout,) ' ** kiss_seed called with args:'
2655	IF(lwp) WRITE(numout,) ' ** ix = ',ix
2656	IF(lwp) WRITE(numout,) ' ** iy = ',iy
2657	IF(lwp) WRITE(numout,) ' ** iz = ',iz
2658	IF(lwp) WRITE(numout,) ' ** iw = ',iw
2659
2660	x = ix
2661	y = iy
2662	z = iz
2663	w = iw
2664
2665	END SUBROUTINE kiss_seed
2666
2667
2668	SUBROUTINE kiss_state(ix, iy, iz, iw)
2669	!! --------------------------------------------------------------------
2670	!! * ROUTINE kiss_state *
2671	!!
2672	!! ** Purpose : Get current state of KISS random number generator
2673	!!
2674	!! --------------------------------------------------------------------
2675	IMPLICIT NONE
2676	INTEGER(KIND=i8) :: ix, iy, iz, iw
2677
2678	ix = x
2679	iy = y
2680	iz = z
2681	iw = w
2682
2683	END SUBROUTINE kiss_state
2684
2685
2686	SUBROUTINE kiss_reset()
2687	!! --------------------------------------------------------------------
2688	!! * ROUTINE kiss_reset *
2689	!!
2690	!! ** Purpose : Reset the default seeds for KISS random number generator
2691	!!
2692	!! --------------------------------------------------------------------
2693	IMPLICIT NONE
2694
2695	x=1234567890987654321_8
2696	y=362436362436362436_8
2697	z=1066149217761810_8
2698	w=123456123456123456_8
2699
2700	END SUBROUTINE kiss_reset
2701
2702	SUBROUTINE kiss_uniform(uran)
2703	!! --------------------------------------------------------------------
2704	!! * ROUTINE kiss_uniform *
2705	!!
2706	!! ** Purpose : Real random numbers with uniform distribution in [0,1]
2707	!!
2708	!! --------------------------------------------------------------------
2709	IMPLICIT NONE
2710	REAL(KIND=wp) :: uran
2711
2712	uran = half * ( one + REAL(kiss(),wp) / huge64 )
2713
2714	END SUBROUTINE kiss_uniform
2715
2716
2717	SUBROUTINE kiss_gaussian(gran)
2718	!! --------------------------------------------------------------------
2719	!! * ROUTINE kiss_gaussian *
2720	!!
2721	!! ** Purpose : Real random numbers with Gaussian distribution N(0,1)
2722	!!
2723	!! ** Method : Generate 2 new Gaussian draws (gran1 and gran2)
2724	!! from 2 uniform draws on [-1,1] (u1 and u2),
2725	!! using the Marsaglia polar method
2726	!! (see Devroye, Non-Uniform Random Variate Generation, p. 235-236)
2727	!!
2728	!! --------------------------------------------------------------------
2729	IMPLICIT NONE
2730	REAL(KIND=wp) :: gran, u1, u2, rsq, fac
2731
2732	IF (ig.EQ.1) THEN
2733	rsq = two
2734	DO WHILE ( (rsq.GE.one).OR. (rsq.EQ.zero) )
2735	u1 = REAL(kiss(),wp) / huge64
2736	u2 = REAL(kiss(),wp) / huge64
2737	rsq = u1u1 + u2u2
2738	ENDDO
2739	fac = SQRT(-two*LOG(rsq)/rsq)
2740	gran1 = u1 * fac
2741	gran2 = u2 * fac
2742	ENDIF
2743
2744	! Output one of the 2 draws
2745	IF (ig.EQ.1) THEN
2746	gran = gran1 ; ig = 2
2747	ELSE
2748	gran = gran2 ; ig = 1
2749	ENDIF
2750
2751	END SUBROUTINE kiss_gaussian
2752
2753
2754	SUBROUTINE kiss_gamma(gamr,k)
2755	!! --------------------------------------------------------------------
2756	!! * ROUTINE kiss_gamma *
2757	!!
2758	!! ** Purpose : Real random numbers with Gamma distribution Gamma(k,1)
2759	!!
2760	!! --------------------------------------------------------------------
2761	IMPLICIT NONE
2762	REAL(KIND=wp), PARAMETER :: p1 = 4.5_8
2763	REAL(KIND=wp), PARAMETER :: p2 = 2.50407739677627_8 ! 1+LOG(9/2)
2764	REAL(KIND=wp), PARAMETER :: p3 = 1.38629436111989_8 ! LOG(4)
2765	REAL(KIND=wp) :: gamr, k, u1, u2, b, c, d, xx, yy, zz, rr, ee
2766	LOGICAL :: accepted
2767
2768	IF (k.GT.one) THEN
2769	! Cheng's rejection algorithm
2770	! (see Devroye, Non-Uniform Random Variate Generation, p. 413)
2771	b = k - p3 ; d = SQRT(two*k-one) ; c = k + d
2772
2773	accepted=.FALSE.
2774	DO WHILE (.NOT.accepted)
2775	CALL kiss_uniform(u1)
2776	yy = LOG(u1/(one-u1)) / d ! Mistake in Devroye: "* k" instead of "/ d"
2777	xx = k * EXP(yy)
2778	rr = b + c * yy - xx
2779	CALL kiss_uniform(u2)
2780	zz = u1 * u1 * u2
2781
2782	accepted = rr .GE. (zz*p1-p2)
2783	IF (.NOT.accepted) accepted = rr .GE. LOG(zz)
2784	ENDDO
2785
2786	gamr = xx
2787
2788	ELSEIF (k.LT.one) THEN
2789	! Rejection from the Weibull density
2790	! (see Devroye, Non-Uniform Random Variate Generation, p. 415)
2791	c = one/k ; d = (one-k) * EXP( (k/(one-k)) * LOG(k) )
2792
2793	accepted=.FALSE.
2794	DO WHILE (.NOT.accepted)
2795	CALL kiss_uniform(u1)
2796	zz = -LOG(u1)
2797	xx = EXP( c * LOG(zz) )
2798	CALL kiss_uniform(u2)
2799	ee = -LOG(u2)
2800
2801	accepted = (zz+ee) .GE. (d+xx) ! Mistake in Devroye: "LE" instead of "GE"
2802	ENDDO
2803
2804	gamr = xx
2805
2806	ELSE
2807	! Exponential distribution
2808	CALL kiss_uniform(u1)
2809	gamr = -LOG(u1)
2810
2811	ENDIF
2812
2813	END SUBROUTINE kiss_gamma
2814
2815
2816	SUBROUTINE kiss_sample(a,n,k)
2817	!! --------------------------------------------------------------------
2818	!! * ROUTINE kiss_sample *
2819	!!
2820	!! ** Purpose : Select a random sample of size k from a set of n integers
2821	!!
2822	!! ** Method : The sample is output in the first k elements of a
2823	!! Set k equal to n to obtain a random permutation
2824	!! of the whole set of integers
2825	!!
2826	!! --------------------------------------------------------------------
2827	IMPLICIT NONE
2828	INTEGER(KIND=i8), DIMENSION(:) :: a
2829	INTEGER(KIND=i8) :: n, k, i, j, atmp
2830	REAL(KIND=wp) :: uran
2831
2832	! Select the sample using the swapping method
2833	! (see Devroye, Non-Uniform Random Variate Generation, p. 612)
2834	DO i=1,k
2835	! Randomly select the swapping element between i and n (inclusive)
2836	CALL kiss_uniform(uran)
2837	j = i - 1 + CEILING( REAL(n-i+1,8) * uran )
2838	! Swap elements i and j
2839	atmp = a(i) ; a(i) = a(j) ; a(j) = atmp
2840	ENDDO
2841
2842	END SUBROUTINE kiss_sample
2843
2844	#ifdef NEMO_V34
2845	SUBROUTINE ctl_nam(ier,cstr,lout)
2846	INTEGER, INTENT(IN) :: ier
2847	LOGICAL, INTENT(IN) :: lout
2848	CHARACTER(LEN=*),INTENT(IN) :: cstr
2849	IF(lout) WRITE(numout,'(A,I4,X,A)') 'Error ',ier,TRIM(cstr)
2850	END SUBROUTINE
2851	#endif
2852
2853	END MODULE stopack

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

Download in other formats: