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

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

Updates for 1d runnig

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