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

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

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

Context Navigation

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

Last change on this file since 13320 was 13320, checked in by jwhile, 4 years ago
Changes for Matt's review
File size: 114.3 KB

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

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

Download in other formats: