Context Navigation

stopack.F90 @ 13355

Last change on this file since 13355 was 13355, checked in by jwhile, 4 years ago
Merged in changes to fix 1d running - documented in UKMO ocean utils ticket 367
File size: 113.7 KB

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

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

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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/STO/stopack.F90 @ 13355

Download in other formats: