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