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

Last change on this file since 12102 was 12102, checked in by dford, 14 months ago

Fix syntax error and relax diagnostic check. See Met Office utils ticket 195.

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