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

Last change on this file since 13322 was 13322, checked in by jwhile, 5 months ago

Updates to how STO deals with 1d case.

Refactor of stopack, as per review

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