| 1 | MODULE nonoro_gwd_ran_mod |
|---|
| 2 | |
|---|
| 3 | IMPLICIT NONE |
|---|
| 4 | |
|---|
| 5 | REAL, allocatable, save :: du_nonoro_gwd(:, :) ! Zonal wind tendency due to GWD |
|---|
| 6 | REAL, allocatable, save :: dv_nonoro_gwd(:, :) ! Meridional wind tendency due to GWD |
|---|
| 7 | REAL, allocatable, save :: east_gwstress(:, :) ! Eastward stress profile |
|---|
| 8 | REAL, allocatable, save :: west_gwstress(:, :) ! Westward stress profile |
|---|
| 9 | |
|---|
| 10 | CONTAINS |
|---|
| 11 | |
|---|
| 12 | SUBROUTINE nonoro_gwd_ran(ngrid, nlayer, dtime, pp, & |
|---|
| 13 | zmax_therm, pt, pu, pv, pdt, pdu, pdv, & |
|---|
| 14 | zustr, zvstr, d_t, d_u, d_v) |
|---|
| 15 | |
|---|
| 16 | !-------------------------------------------------------------------------------- |
|---|
| 17 | ! Parametrization of the momentum flux deposition due to a discrete |
|---|
| 18 | ! number of gravity waves. |
|---|
| 19 | ! F. Lott |
|---|
| 20 | ! Version 14, Gaussian distribution of the source |
|---|
| 21 | ! LMDz model online version |
|---|
| 22 | ! ADAPTED FOR VENUS / F. LOTT + S. LEBONNOIS |
|---|
| 23 | ! Version adapted on 03/04/2013: |
|---|
| 24 | ! - input flux compensated in the deepest layers |
|---|
| 25 | ! |
|---|
| 26 | ! ADAPTED FOR MARS G.GILLI 02/2016 |
|---|
| 27 | ! Revision with F.Forget 06/2016 Variable EP-flux |
|---|
| 28 | ! according to |
|---|
| 29 | ! PBL variation (max |
|---|
| 30 | ! velocity thermals) |
|---|
| 31 | ! D.BARDET 01/2020 - reproductibility of |
|---|
| 32 | ! the launching altitude |
|---|
| 33 | ! calculation |
|---|
| 34 | ! - wave characteristic |
|---|
| 35 | ! calculation using MOD |
|---|
| 36 | ! - adding east_gwstress |
|---|
| 37 | ! and west_gwstress variables |
|---|
| 38 | ! |
|---|
| 39 | ! ADAPTED FOR GENERIC D.BARDET 01/2020 |
|---|
| 40 | !--------------------------------------------------------------------------------- |
|---|
| 41 | |
|---|
| 42 | |
|---|
| 43 | |
|---|
| 44 | use comcstfi_mod, only: g, pi, cpp, r |
|---|
| 45 | USE ioipsl_getin_p_mod, ONLY : getin_p |
|---|
| 46 | use assert_m, only : assert |
|---|
| 47 | use vertical_layers_mod, only : presnivs |
|---|
| 48 | use callkeys_mod, only : epflux_max, & ! Max EP flux value at launching altitude (previous name: RUWMAX) |
|---|
| 49 | gwd_convective_source |
|---|
| 50 | use inifis_mod |
|---|
| 51 | use xios_output_mod, only: send_xios_field |
|---|
| 52 | implicit none |
|---|
| 53 | |
|---|
| 54 | ! 0. DECLARATIONS: |
|---|
| 55 | |
|---|
| 56 | ! 0.1 INPUTS |
|---|
| 57 | |
|---|
| 58 | INTEGER, intent(in):: ngrid, nlayer |
|---|
| 59 | REAL, intent(in):: dtime ! Time step of the Physics |
|---|
| 60 | REAL, intent(in):: zmax_therm(ngrid) ! Altitude of max velocity thermals (m) |
|---|
| 61 | REAL, intent(in):: pp(ngrid, nlayer) ! Pressure at full levels |
|---|
| 62 | REAL, intent(in):: pt(ngrid, nlayer) ! Temperature at full levels |
|---|
| 63 | REAL, intent(in):: pu(ngrid, nlayer) ! Zonal wind at full levels |
|---|
| 64 | REAL, intent(in):: pv(ngrid, nlayer) ! Meridional wind at full levels |
|---|
| 65 | REAL, intent(in):: pdt(ngrid, nlayer) ! Tendency on temperature |
|---|
| 66 | REAL, intent(in):: pdu(ngrid, nlayer) ! Tendency on zonal wind |
|---|
| 67 | REAL, intent(in):: pdv(ngrid, nlayer) ! Tendency on meridional wind |
|---|
| 68 | |
|---|
| 69 | ! 0.2 OUTPUTS |
|---|
| 70 | |
|---|
| 71 | REAL, intent(out):: zustr(ngrid) ! Zonal surface stress |
|---|
| 72 | REAL, intent(out):: zvstr(ngrid) ! Meridional surface stress |
|---|
| 73 | REAL, intent(out):: d_t(ngrid, nlayer) ! Tendency on temperature |
|---|
| 74 | REAL, intent(out):: d_u(ngrid, nlayer) ! Tendency on zonal wind |
|---|
| 75 | REAL, intent(out):: d_v(ngrid, nlayer) ! Tendency on meridional wind |
|---|
| 76 | |
|---|
| 77 | ! 0.3 INTERNAL ARRAYS |
|---|
| 78 | |
|---|
| 79 | REAL :: tt(ngrid, nlayer) ! Temperature at full levels |
|---|
| 80 | REAL :: uu(ngrid, nlayer) ! Zonal wind at full levels |
|---|
| 81 | REAL :: vv(ngrid, nlayer) ! Meridional wind at full levels |
|---|
| 82 | REAL :: bvlow(ngrid) ! Qu est-ce ? |
|---|
| 83 | REAL :: dz |
|---|
| 84 | |
|---|
| 85 | INTEGER ii, jj, ll |
|---|
| 86 | |
|---|
| 87 | ! 0.3.0 Time scale of the like cycle of the waves parametrized |
|---|
| 88 | REAL, parameter:: deltat = 24. * 3600. |
|---|
| 89 | |
|---|
| 90 | ! 0.3.1 Gravity waves specifications |
|---|
| 91 | INTEGER, parameter:: nk = 2 ! number of horizontal wavenumbers |
|---|
| 92 | INTEGER, parameter:: np = 2 ! directions (eastward and westward) phase speed |
|---|
| 93 | INTEGER, parameter:: no = 2 ! absolute values of phase speed |
|---|
| 94 | INTEGER, parameter:: na = 5 ! Number of realizations to get the phase speed |
|---|
| 95 | INTEGER, parameter:: nw = nk * np *no ! Total numbers of gravity waves |
|---|
| 96 | INTEGER jk, jp, jo, jw |
|---|
| 97 | REAL, parameter:: kmax = 7.e-4 ! Max horizontal wavenumber |
|---|
| 98 | REAL, parameter:: kmin = 2.e-5 ! Min horizontal wavenumber |
|---|
| 99 | !REAL, parameter:: cmax = 30. ! Max horizontal absolute phase velocity |
|---|
| 100 | REAL, parameter:: cmax = 100. ! Test for Saturn: Max horizontal absolute phase velocity |
|---|
| 101 | REAL, parameter:: cmin = 1. ! Min horizontal absolute phase velocity |
|---|
| 102 | REAL cpha ! absolute phase velocity frequency |
|---|
| 103 | REAL zk(nw, ngrid) ! horizontal wavenumber amplitude |
|---|
| 104 | REAL zp(nw, ngrid) ! horizontal wavenumber angle |
|---|
| 105 | REAL zo(nw, ngrid) ! absolute frequency |
|---|
| 106 | |
|---|
| 107 | REAL intr_freq_m(nw, ngrid) ! Waves Intr. freq. at the 1/2 lev below the full level (previous name: ZOM) |
|---|
| 108 | REAL intr_freq_p(nw, ngrid) ! Waves Intr. freq. at the 1/2 lev above the full level (previous name: ZOP) |
|---|
| 109 | REAL wwm(nw, ngrid) ! Wave EP-fluxes at the 1/2 level below the full level |
|---|
| 110 | REAL wwp(nw, ngrid) ! Wave EP-fluxes at the 1/2 level above the full level |
|---|
| 111 | REAL u_epflux_p(nw, ngrid) ! Partial zonal flux (=for each wave) at the 1/2 level above the full level (previous name: RUWP) |
|---|
| 112 | REAL v_epflux_p(nw, ngrid) ! Partial meridional flux (=for each wave) at the 1/2 level above the full level (previous name: RVWP) |
|---|
| 113 | REAL u_epflux_tot(ngrid, nlayer + 1) ! Total zonal flux (=for all waves (nw)) at the 1/2 level above the full level (3D) (previous name: RUW) |
|---|
| 114 | REAL v_epflux_tot(ngrid, nlayer + 1) ! Total meridional flux (=for all waves (nw)) at the 1/2 level above the full level (3D) (previous name: RVW) |
|---|
| 115 | REAL epflux_0(nw, ngrid) ! Fluxes at launching level (previous name: RUW0) |
|---|
| 116 | INTEGER launch ! Launching altitude |
|---|
| 117 | REAL, parameter:: xlaunch = 0.4 ! Control the launching altitude |
|---|
| 118 | REAL, parameter:: zmaxth_top = 8000. ! Top of convective layer (approx.) |
|---|
| 119 | |
|---|
| 120 | REAL prec(ngrid) ! precipitations |
|---|
| 121 | REAL prmax ! Max value of precipitation |
|---|
| 122 | |
|---|
| 123 | ! 0.3.2 Parameters of waves dissipations |
|---|
| 124 | REAL, parameter:: sat = 1. ! saturation parameter |
|---|
| 125 | REAL, parameter:: rdiss = 1. ! coefficient of dissipation |
|---|
| 126 | REAL, parameter:: zoisec = 1.e-6 ! security for intrinsic freguency |
|---|
| 127 | |
|---|
| 128 | ! 0.3.3 Background flow at 1/2 levels and vertical coordinate |
|---|
| 129 | REAL H0bis(ngrid, nlayer) ! characteristic height of the atmosphere |
|---|
| 130 | REAL, save:: H0 ! characteristic height of the atmosphere |
|---|
| 131 | REAL, parameter:: pr = 250 ! Reference pressure [Pa] |
|---|
| 132 | REAL, parameter:: tr = 220. ! Reference temperature [K] |
|---|
| 133 | REAL zh(ngrid, nlayer + 1) ! Log-pressure altitude (constant H0) |
|---|
| 134 | REAL zhbis(ngrid, nlayer + 1) ! Log-pressure altitude (varying H) |
|---|
| 135 | REAL uh(ngrid, nlayer + 1) ! Zonal wind at 1/2 levels |
|---|
| 136 | REAL vh(ngrid, nlayer + 1) ! Meridional wind at 1/2 levels |
|---|
| 137 | REAL ph(ngrid, nlayer + 1) ! Pressure at 1/2 levels |
|---|
| 138 | REAL, parameter:: psec = 1.e-6 ! Security to avoid division by 0 pressure |
|---|
| 139 | REAL bv(ngrid, nlayer + 1) ! Brunt Vaisala freq. at 1/2 levels |
|---|
| 140 | REAL, parameter:: bvsec = 1.e-5 ! Security to avoid negative BV |
|---|
| 141 | REAL href(nlayer + 1) ! Reference altitude for launching altitude |
|---|
| 142 | |
|---|
| 143 | |
|---|
| 144 | |
|---|
| 145 | REAL ran_num_1, ran_num_2, ran_num_3 |
|---|
| 146 | |
|---|
| 147 | |
|---|
| 148 | LOGICAL, save :: firstcall = .true. |
|---|
| 149 | |
|---|
| 150 | !-------------------------------------------------------------------------------------------------- |
|---|
| 151 | ! 1. INITIALISATION |
|---|
| 152 | !------------------ |
|---|
| 153 | IF (firstcall) THEN |
|---|
| 154 | write(*,*) "nonoro_gwd_ran: FLott non-oro GW scheme is active!" |
|---|
| 155 | ! Characteristic vertical scale height |
|---|
| 156 | H0 = r * tr / g |
|---|
| 157 | ! Control |
|---|
| 158 | if (deltat .LT. dtime) THEN |
|---|
| 159 | call abort_physic("nonoro_gwd_ran","gwd random: deltat lower than dtime!",1) |
|---|
| 160 | endif |
|---|
| 161 | if (nlayer .LT. nw) THEN |
|---|
| 162 | call abort_physic("nonoro_gwd_ran","gwd random: nlayer lower than nw!",1) |
|---|
| 163 | endif |
|---|
| 164 | firstcall = .false. |
|---|
| 165 | ENDIF |
|---|
| 166 | gwd_convective_source = .false. |
|---|
| 167 | |
|---|
| 168 | ! Compute subroutine's current values of temperature and winds |
|---|
| 169 | tt(:,:) = pt(:,:) + dtime * pdt(:,:) |
|---|
| 170 | uu(:,:) = pu(:,:) + dtime * pdu(:,:) |
|---|
| 171 | vv(:,:) = pv(:,:) + dtime * pdv(:,:) |
|---|
| 172 | ! print*,'epflux_max just after firstcall:', epflux_max |
|---|
| 173 | |
|---|
| 174 | |
|---|
| 175 | ! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS |
|---|
| 176 | !---------------------------------------------------- |
|---|
| 177 | ! Pressure and Inv of pressure, temperature at 1/2 level |
|---|
| 178 | DO ll = 2, nlayer |
|---|
| 179 | ph(:, ll) = exp((log(pp(:, ll)) + log(pp(:, ll - 1))) / 2.) |
|---|
| 180 | end DO |
|---|
| 181 | |
|---|
| 182 | ph(:, nlayer + 1) = 0. |
|---|
| 183 | ph(:, 1) = 2. * pp(:, 1) - ph(:, 2) |
|---|
| 184 | |
|---|
| 185 | ! Launching altitude for reproductible case |
|---|
| 186 | DO ll = 2, nlayer |
|---|
| 187 | href(ll) = exp((log(presnivs(ll)) + log(presnivs(ll - 1))) / 2.) |
|---|
| 188 | end DO |
|---|
| 189 | href(nlayer + 1) = 0. |
|---|
| 190 | href(1) = 2. * presnivs(1) - href(2) |
|---|
| 191 | |
|---|
| 192 | launch = 0. |
|---|
| 193 | DO ll =1, nlayer |
|---|
| 194 | IF (href (ll) / href(1) > xlaunch) launch = ll |
|---|
| 195 | end DO |
|---|
| 196 | |
|---|
| 197 | ! Log-pressure vertical coordinate |
|---|
| 198 | DO ll = 1, nlayer + 1 |
|---|
| 199 | zh(:,ll) = H0 * log(pr / (ph(:,ll) + psec)) |
|---|
| 200 | end DO |
|---|
| 201 | |
|---|
| 202 | ! Winds and Brunt Vaisala frequency |
|---|
| 203 | DO ll = 2, nlayer |
|---|
| 204 | uh(:, ll) = 0.5 * (uu(:, ll) + uu(:, ll - 1)) ! Zonal wind |
|---|
| 205 | vh(:, ll) = 0.5 * (vv(:, ll) + vv(:, ll - 1)) ! Meridional wind |
|---|
| 206 | |
|---|
| 207 | bv(:, ll) = 0.5 * (tt(:, ll) + tt(:, ll - 1)) & |
|---|
| 208 | * r**2 / cpp / H0**2 & |
|---|
| 209 | + (tt(:, ll) - tt(:, ll - 1)) & |
|---|
| 210 | / (zh(:, ll) - zh(:, ll - 1)) & |
|---|
| 211 | * r / H0 |
|---|
| 212 | bv(:, ll) = sqrt(max(bvsec, bv(:,ll))) |
|---|
| 213 | end DO |
|---|
| 214 | |
|---|
| 215 | bv(:, 1) = bv(:, 2) |
|---|
| 216 | uh(:, 1) = 0. |
|---|
| 217 | vh(:, 1) = 0. |
|---|
| 218 | bv(:, nlayer + 1) = bv(:, nlayer) |
|---|
| 219 | uh(:, nlayer + 1) = uu(:, nlayer) |
|---|
| 220 | vh(:, nlayer + 1) = vv(:, nlayer) |
|---|
| 221 | |
|---|
| 222 | ! 3. WAVES CHARACTERISTICS CHOSEN RANDOMLY |
|---|
| 223 | !----------------------------------------- |
|---|
| 224 | ! The mod function of here a weird arguments |
|---|
| 225 | ! are used to produce the waves characteristics |
|---|
| 226 | ! in a stochastic way |
|---|
| 227 | DO jw = 1, nw |
|---|
| 228 | DO ii = 1, ngrid |
|---|
| 229 | |
|---|
| 230 | ran_num_1 = mod(tt(ii, jw) * 10., 1.) |
|---|
| 231 | ran_num_2 = mod(tt(ii, jw) * 100., 1.) |
|---|
| 232 | |
|---|
| 233 | ! angle (0 or pi so far) |
|---|
| 234 | zp(jw, ii) = (sign(1., 0.5 - ran_num_1) & |
|---|
| 235 | + 1.) * pi / 2. |
|---|
| 236 | ! horizontal wavenumber amplitude |
|---|
| 237 | zk(jw, ii) = kmin + (kmax - kmin) * ran_num_2 |
|---|
| 238 | ! horizontal phase speed |
|---|
| 239 | cpha = 0. |
|---|
| 240 | DO jj = 1, na |
|---|
| 241 | ran_num_3 = mod(tt(ii, jw + 3 * jj)**2, 1.) |
|---|
| 242 | cpha = cpha + 2. * cmax * & |
|---|
| 243 | (ran_num_3 - 0.5) * & |
|---|
| 244 | sqrt(3.) / sqrt(na * 1.) |
|---|
| 245 | end DO |
|---|
| 246 | IF (cpha < 0.) THEN |
|---|
| 247 | cpha = - 1. * cpha |
|---|
| 248 | zp (jw, ii) = zp(jw, ii) + pi |
|---|
| 249 | ENDIF |
|---|
| 250 | ! Intrinsic frequency |
|---|
| 251 | zo(jw, ii) = cpha * zk(jw, ii) |
|---|
| 252 | ! Intrinsic frequency is imposed |
|---|
| 253 | zo(jw, ii) = zo(jw, ii) & |
|---|
| 254 | + zk(jw, ii) * cos(zp(jw, ii)) * uh(ii, launch) & |
|---|
| 255 | + zk(jw, ii) * sin(zp(jw, ii)) * vh(ii, launch) |
|---|
| 256 | ! Momentum flux at launch level |
|---|
| 257 | epflux_0(jw, ii) = epflux_max & |
|---|
| 258 | * mod(100. * (uu(ii, jw)**2 + vv(ii, jw)**2), 1.) |
|---|
| 259 | |
|---|
| 260 | end DO |
|---|
| 261 | end DO |
|---|
| 262 | |
|---|
| 263 | ! print*,'epflux_max just after waves charac. ramdon:', epflux_max |
|---|
| 264 | ! print*,'epflux_0 just after waves charac. ramdon:', epflux_0 |
|---|
| 265 | ! 4. COMPUTE THE FLUXES |
|---|
| 266 | !---------------------- |
|---|
| 267 | ! 4.1 Vertical velocity at launching altitude to ensure |
|---|
| 268 | ! the correct value to the imposed fluxes. |
|---|
| 269 | !------------------------------------------------------ |
|---|
| 270 | DO jw = 1, nw |
|---|
| 271 | ! Evaluate intrinsic frequency at launching altutide: |
|---|
| 272 | intr_freq_p(jw, :) = zo(jw, :) & |
|---|
| 273 | - zk(jw, :) * cos(zp(jw, :)) * uh(:, launch) & |
|---|
| 274 | - zk(jw, :) * sin(zp(jw, :)) * vh(:, launch) |
|---|
| 275 | end DO |
|---|
| 276 | |
|---|
| 277 | IF (gwd_convective_source) THEN |
|---|
| 278 | DO jw = 1, nw |
|---|
| 279 | ! VERSION WITH CONVECTIVE SOURCE ! designed for Earth |
|---|
| 280 | |
|---|
| 281 | ! Vertical velocity at launch level, value to ensure the |
|---|
| 282 | ! imposed mmt flux factor related to the convective forcing: |
|---|
| 283 | ! precipitations. |
|---|
| 284 | |
|---|
| 285 | ! tanh limitation to values above prmax: |
|---|
| 286 | ! WWP(JW, :) = epflux_0(JW, :) & |
|---|
| 287 | ! * (r / cpp / H0 * RLVTT * PRMAX * TANH(PREC(:) / PRMAX))**2 |
|---|
| 288 | ! Here, we neglected the kinetic energy providing of the thermodynamic |
|---|
| 289 | ! phase change |
|---|
| 290 | |
|---|
| 291 | ! |
|---|
| 292 | |
|---|
| 293 | ! Factor related to the characteristics of the waves: |
|---|
| 294 | wwp(jw, :) = wwp(jw, :) * zk(jw, :)**3 / kmin / bvlow(:) & |
|---|
| 295 | / MAX(ABS(intr_freq_p(jw, :)), zoisec)**3 |
|---|
| 296 | |
|---|
| 297 | ! Moderation by the depth of the source (dz here): |
|---|
| 298 | wwp(jw, :) = wwp(jw, :) & |
|---|
| 299 | * exp(- bvlow(:)**2 / max(abs(intr_freq_p(jw, :)), zoisec)**2 & |
|---|
| 300 | * zk(jw, :)**2 * dz**2) |
|---|
| 301 | |
|---|
| 302 | ! Put the stress in the right direction: |
|---|
| 303 | u_epflux_p(jw, :) = intr_freq_p(jw, :) / max(abs(intr_freq_p(jw, :)), zoisec)**2 & |
|---|
| 304 | * bv(:, launch) * cos(zp(jw, :)) * wwp(jw, :)**2 |
|---|
| 305 | v_epflux_p(JW, :) = intr_freq_p(jw, :) / max(abs(intr_freq_p(jw, :)), zoisec)**2 & |
|---|
| 306 | * bv(:, launch) * sin(zp(jw, :)) * wwp(jw, :)**2 |
|---|
| 307 | end DO |
|---|
| 308 | ELSE ! VERSION WITHOUT CONVECTIVE SOURCE |
|---|
| 309 | ! Vertical velocity at launch level, value to ensure the imposed |
|---|
| 310 | ! mom flux: |
|---|
| 311 | DO jw = 1, nw |
|---|
| 312 | ! WW is directly a flux, here, not vertical velocity anymore |
|---|
| 313 | wwp(jw, :) = epflux_0(JW,:) |
|---|
| 314 | u_epflux_p(jw, :) = cos(zp(jw, :)) * sign(1., intr_freq_p(jw, :)) * epflux_0(jw, :) |
|---|
| 315 | v_epflux_p(jw, :) = sin(zp(jw, :)) * sign(1., intr_freq_p(jw, :)) * epflux_0(jw, :) |
|---|
| 316 | |
|---|
| 317 | end DO |
|---|
| 318 | ENDIF |
|---|
| 319 | |
|---|
| 320 | ! 4.2 Initial flux at launching altitude |
|---|
| 321 | !--------------------------------------- |
|---|
| 322 | u_epflux_tot(:, launch) = 0. |
|---|
| 323 | v_epflux_tot(:, launch) = 0. |
|---|
| 324 | DO jw = 1, nw |
|---|
| 325 | u_epflux_tot(:, launch) = u_epflux_tot(:, launch) + u_epflux_p(jw, :) |
|---|
| 326 | v_epflux_tot(:, launch) = v_epflux_tot(:, launch) + v_epflux_p(jw, :) |
|---|
| 327 | end DO |
|---|
| 328 | |
|---|
| 329 | ! 4.3 Loop over altitudes, with passage from one level to the next done by: |
|---|
| 330 | !-------------------------------------------------------------------------- |
|---|
| 331 | ! i) conserving the EP flux, |
|---|
| 332 | ! ii) dissipating a little, |
|---|
| 333 | ! iii) testing critical levels, |
|---|
| 334 | ! iv) testing the breaking. |
|---|
| 335 | !-------------------------------------------------------------------------- |
|---|
| 336 | DO ll = launch, nlayer - 1 |
|---|
| 337 | ! Warning! all the physics is here (passage from one level to the next |
|---|
| 338 | DO jw = 1, nw |
|---|
| 339 | intr_freq_m(jw, :) = intr_freq_p(jw, :) |
|---|
| 340 | wwm(jw, :) = wwp(jw, :) |
|---|
| 341 | ! Intrinsic frequency |
|---|
| 342 | intr_freq_p(jw, :) = zo(jw, :) - zk(jw, :) * cos(zp(jw, :)) * uh(:, ll + 1) & |
|---|
| 343 | - zk(jw, :) * sin(zp(jw, :)) * vh(:, ll + 1) |
|---|
| 344 | ! Vertical velocity in flux formulation |
|---|
| 345 | wwp(jw, :) = min( & |
|---|
| 346 | ! No breaking (Eq. 6): |
|---|
| 347 | wwm(jw, :) & |
|---|
| 348 | ! Dissipation (Eq. 8): |
|---|
| 349 | * exp(-rdiss * pr / (ph(:, ll + 1) + ph (:, ll)) & |
|---|
| 350 | * ((bv(:, ll + 1) + bv (:, ll)) / 2.)**3 & |
|---|
| 351 | / max(abs(intr_freq_p(jw, :) + intr_freq_m(jw, :)) / 2., zoisec)**4 & |
|---|
| 352 | * zk(jw, :)**3 * (zh(:, ll + 1) - zh(:, ll))) , & |
|---|
| 353 | ! Critical levels (forced to zero if intrinsic frequency |
|---|
| 354 | ! changes sign) |
|---|
| 355 | max(0., sign(1., intr_freq_p(jw, :) * intr_freq_m(jw, :))) & |
|---|
| 356 | ! Saturation (Eq. 12) |
|---|
| 357 | * abs(intr_freq_p(jw, :))**3 / bv(:, ll + 1) & |
|---|
| 358 | * exp(-zh(:, ll + 1) / H0) & |
|---|
| 359 | * sat**2 * kmin**2 / zk(jw, :)**4) |
|---|
| 360 | end DO |
|---|
| 361 | |
|---|
| 362 | ! Evaluate EP-flux from Eq. 7 and give the right orientation to the stress |
|---|
| 363 | DO jw = 1, nw |
|---|
| 364 | u_epflux_p(jw, :) = sign(1., intr_freq_p(jw, :)) * cos(zp(jw, :)) * wwp(jw, :) |
|---|
| 365 | v_epflux_p(jw, :) = sign(1., intr_freq_p(jw, :)) * sin(zp(jw, :)) * wwp(jw, :) |
|---|
| 366 | end DO |
|---|
| 367 | |
|---|
| 368 | u_epflux_tot(:, ll + 1) = 0. |
|---|
| 369 | v_epflux_tot(:, ll + 1) = 0. |
|---|
| 370 | DO jw = 1, nw |
|---|
| 371 | u_epflux_tot(:, ll + 1) = u_epflux_tot(:, ll + 1) + u_epflux_p(jw, :) |
|---|
| 372 | v_epflux_tot(:, ll + 1) = v_epflux_tot(:, ll + 1) + v_epflux_p(jw, :) |
|---|
| 373 | east_gwstress(:, ll) = east_gwstress(:, ll) & |
|---|
| 374 | + max(0., u_epflux_p(jw, :)) / float(nw) |
|---|
| 375 | west_gwstress(:, ll) = west_gwstress(:, ll) & |
|---|
| 376 | + min(0., u_epflux_p(jw, ll)) / float(nw) |
|---|
| 377 | end DO |
|---|
| 378 | end DO ! DO LL = LAUNCH, nlayer - 1 |
|---|
| 379 | |
|---|
| 380 | ! print*,'u_epflux_tot just after launching:', u_epflux_tot |
|---|
| 381 | ! print*,'v_epflux_tot just after launching:', v_epflux_tot |
|---|
| 382 | ! print*,'u_epflux_p just after launching:', maxval(u_epflux_p), minval(u_epflux_p) |
|---|
| 383 | ! print*,'v_epflux_p just after launching:', maxval(v_epflux_p), minval(v_epflux_p) |
|---|
| 384 | |
|---|
| 385 | ! 5. TENDENCY CALCULATIONS |
|---|
| 386 | !------------------------- |
|---|
| 387 | ! 5.1 Flow rectification at the top and in the low layers |
|---|
| 388 | ! -------------------------------------------------------- |
|---|
| 389 | ! Warning, this is the total on all GW... |
|---|
| 390 | |
|---|
| 391 | u_epflux_tot(:, nlayer + 1) = 0. |
|---|
| 392 | v_epflux_tot(:, nlayer + 1) = 0. |
|---|
| 393 | |
|---|
| 394 | ! Here, big change compared to FLott version: |
|---|
| 395 | ! We compensate (u_epflux_(:, LAUNCH), ie total emitted upward flux |
|---|
| 396 | ! over the layers max(1,LAUNCH-3) to LAUNCH-1 |
|---|
| 397 | DO LL = 1, max(1,LAUNCH-3) |
|---|
| 398 | u_epflux_tot(:, LL) = 0. |
|---|
| 399 | v_epflux_tot(:, LL) = 0. |
|---|
| 400 | end DO |
|---|
| 401 | DO LL = max(2,LAUNCH-2), LAUNCH-1 |
|---|
| 402 | u_epflux_tot(:, LL) = u_epflux_tot(:, LL - 1) & |
|---|
| 403 | + u_epflux_tot(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1)) & |
|---|
| 404 | / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3))) |
|---|
| 405 | v_epflux_tot(:, LL) = v_epflux_tot(:, LL - 1) & |
|---|
| 406 | + v_epflux_tot(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1)) & |
|---|
| 407 | / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3))) |
|---|
| 408 | east_gwstress(:,LL) = east_gwstress(:, LL - 1) & |
|---|
| 409 | + east_gwstress(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1)) & |
|---|
| 410 | / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3))) |
|---|
| 411 | west_gwstress(:,LL) = west_gwstress(:, LL - 1) & |
|---|
| 412 | + west_gwstress(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1)) & |
|---|
| 413 | / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3))) |
|---|
| 414 | end DO |
|---|
| 415 | ! This way, the total flux from GW is zero, but there is a net transport |
|---|
| 416 | ! (upward) that should be compensated by circulation |
|---|
| 417 | ! and induce additional friction at the surface |
|---|
| 418 | |
|---|
| 419 | ! 5.2 AR-1 RECURSIVE FORMULA (13) IN VERSION 4 |
|---|
| 420 | !--------------------------------------------- |
|---|
| 421 | DO LL = 1, nlayer |
|---|
| 422 | d_u(:, LL) = (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL))! & |
|---|
| 423 | !d_u(:, LL) = G * (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL))! & |
|---|
| 424 | ! / (PH(:, LL + 1) - PH(:, LL)) * DTIME |
|---|
| 425 | d_v(:, LL) = (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL))! & |
|---|
| 426 | !d_v(:, LL) = G * (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL))! & |
|---|
| 427 | ! / (PH(:, LL + 1) - PH(:, LL)) * DTIME |
|---|
| 428 | ENDDO |
|---|
| 429 | d_t(:,:) = 0. |
|---|
| 430 | |
|---|
| 431 | ! 5.3 Update tendency of wind with the previous (and saved) values |
|---|
| 432 | !----------------------------------------------------------------- |
|---|
| 433 | d_u(:,:) = DTIME/DELTAT/REAL(NW) * d_u(:,:) & |
|---|
| 434 | + (1.-DTIME/DELTAT) * du_nonoro_gwd(:,:) |
|---|
| 435 | d_v(:,:) = DTIME/DELTAT/REAL(NW) * d_v(:,:) & |
|---|
| 436 | + (1.-DTIME/DELTAT) * dv_nonoro_gwd(:,:) |
|---|
| 437 | du_nonoro_gwd(:,:) = d_u(:,:) |
|---|
| 438 | dv_nonoro_gwd(:,:) = d_v(:,:) |
|---|
| 439 | |
|---|
| 440 | ! print*,'u_epflux_tot just after tendency:', u_epflux_tot |
|---|
| 441 | ! print*,'v_epflux_tot just after tendency:', v_epflux_tot |
|---|
| 442 | ! print*,'d_u just after tendency:', maxval(d_u(:,:)), minval(d_u) |
|---|
| 443 | ! print*,'d_v just after tendency:', maxval(d_v(:,:)), minval(d_v) |
|---|
| 444 | ! Cosmetic: evaluation of the cumulated stress |
|---|
| 445 | |
|---|
| 446 | ZUSTR(:) = 0. |
|---|
| 447 | ZVSTR(:) = 0. |
|---|
| 448 | DO LL = 1, nlayer |
|---|
| 449 | ZUSTR(:) = ZUSTR(:) + D_U(:, LL) !/ g * (PH(:, LL + 1) - PH(:, LL)) |
|---|
| 450 | ZVSTR(:) = ZVSTR(:) + D_V(:, LL) !/ g * (PH(:, LL + 1) - PH(:, LL)) |
|---|
| 451 | ENDDO |
|---|
| 452 | |
|---|
| 453 | call send_xios_field("du_nonoro", d_u) |
|---|
| 454 | call send_xios_field("dv_nonoro", d_v) |
|---|
| 455 | call send_xios_field("east_gwstress", east_gwstress) |
|---|
| 456 | call send_xios_field("west_gwstress", west_gwstress) |
|---|
| 457 | |
|---|
| 458 | |
|---|
| 459 | END SUBROUTINE nonoro_gwd_ran |
|---|
| 460 | |
|---|
| 461 | ! =================================================================== |
|---|
| 462 | ! Subroutines used to write variables of memory in start files |
|---|
| 463 | ! =================================================================== |
|---|
| 464 | |
|---|
| 465 | SUBROUTINE ini_nonoro_gwd_ran(ngrid,nlayer) |
|---|
| 466 | |
|---|
| 467 | IMPLICIT NONE |
|---|
| 468 | |
|---|
| 469 | INTEGER, INTENT (in) :: ngrid ! number of atmospheric columns |
|---|
| 470 | INTEGER, INTENT (in) :: nlayer ! number of atmospheric layers |
|---|
| 471 | |
|---|
| 472 | allocate(du_nonoro_gwd(ngrid, nlayer)) |
|---|
| 473 | allocate(dv_nonoro_gwd(ngrid, nlayer)) |
|---|
| 474 | allocate(east_gwstress(ngrid, nlayer)) |
|---|
| 475 | allocate(west_gwstress(ngrid, nlayer)) |
|---|
| 476 | |
|---|
| 477 | END SUBROUTINE ini_nonoro_gwd_ran |
|---|
| 478 | ! ---------------------------------- |
|---|
| 479 | SUBROUTINE end_nonoro_gwd_ran |
|---|
| 480 | |
|---|
| 481 | IMPLICIT NONE |
|---|
| 482 | |
|---|
| 483 | if (allocated(du_nonoro_gwd)) deallocate(du_nonoro_gwd) |
|---|
| 484 | if (allocated(dv_nonoro_gwd)) deallocate(dv_nonoro_gwd) |
|---|
| 485 | if (allocated(east_gwstress)) deallocate(east_gwstress) |
|---|
| 486 | if (allocated(west_gwstress)) deallocate(west_gwstress) |
|---|
| 487 | |
|---|
| 488 | END SUBROUTINE end_nonoro_gwd_ran |
|---|
| 489 | |
|---|
| 490 | END MODULE nonoro_gwd_ran_mod |
|---|
| 491 | |
|---|
