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