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, 5 months 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
55MODULE 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
281CONTAINS
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 = 3*jpl+jpl*nlay_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_wp*zsd*zsd
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_wp*zsd*zsd
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
8009300  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_wp*zsd*zsd
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_wp*zsd*zsd
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
11389301  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
12049301  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 = a*gb + b*g2d
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)
1518901   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 )
1522902   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(jpi*jpj)
1823         WRITE(numout,*) '        Mean AR1 b term : ',SUM(sppt_b)/REAL(jpi*jpj)
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(jpi*jpj)
1856          WRITE(numout,*) '        Mean AR1 b term : ',SUM(spp_b)/REAL(jpi*jpj)
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(jpi*jpj)
1880          WRITE(numout,*) '        Mean AR1 b term : ',SUM(skeb_b)/REAL(jpi*jpj)
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( ds*ds + dt*dt ) * tmask(ji,jj,jk)
2320          dnum(ji,jj,jk) = dtot*dtot*dtot*kh2*e1t(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 = - grav*avt(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) = kc2*zz*e1t(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)*rau0*rau0)
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
24649304 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
25139302  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 = u1*u1 + u2*u2
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
2827SUBROUTINE ctl_nam(ier,cstr,lout)
2828INTEGER, INTENT(IN)            :: ier
2829LOGICAL, INTENT(IN)            :: lout
2830CHARACTER(LEN=*),INTENT(IN)    :: cstr
2831IF(lout) WRITE(numout,'(A,I4,X,A)') 'Error ',ier,TRIM(cstr)
2832END SUBROUTINE
2833#endif
2834
2835END MODULE stopack
Note: See TracBrowser for help on using the repository browser.