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_r5518_GO6_package_FOAMv14_Kw660/NEMOGCM/NEMO/OPA_SRC/STO – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_FOAMv14_Kw660/NEMOGCM/NEMO/OPA_SRC/STO/stopack.F90 @ 11842

Last change on this file since 11842 was 11842, checked in by dford, 4 years ago
Fix to stochastic physics code to compile with key_dynldf_c3d, and relax MEDUSA's Kw660 check.
File size: 113.6 KB

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

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

Download in other formats: