| 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 | !$OMP THREADPRIVATE(du_nonoro_gwd,dv_nonoro_gwd,east_gwstress,west_gwstress) |
|---|
| 10 | CONTAINS |
|---|
| 11 | |
|---|
| 12 | SUBROUTINE nonoro_gwd_ran(ngrid, nlayer, dtime, pp, & |
|---|
| 13 | pn2, presnivs, pt, pu, pv, & |
|---|
| 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 | ! READAPTED FOR VENUS S.LEBONNOIS 08/2021 |
|---|
| 41 | !--------------------------------------------------------------------------------- |
|---|
| 42 | |
|---|
| 43 | USE ioipsl_getin_p_mod, ONLY : getin_p |
|---|
| 44 | USE geometry_mod, ONLY: cell_area, latitude_deg |
|---|
| 45 | USE YOMCST_mod |
|---|
| 46 | !#ifdef CPP_XIOS |
|---|
| 47 | ! use xios_output_mod, only: send_xios_field |
|---|
| 48 | !#endif |
|---|
| 49 | implicit none |
|---|
| 50 | !#include "YOMCST.h" |
|---|
| 51 | |
|---|
| 52 | ! 0. DECLARATIONS: |
|---|
| 53 | |
|---|
| 54 | ! 0.1 INPUTS |
|---|
| 55 | |
|---|
| 56 | INTEGER, intent(in):: ngrid, nlayer |
|---|
| 57 | REAL, intent(in):: dtime ! Time step of the Physics |
|---|
| 58 | REAL, intent(in):: pp(ngrid, nlayer) ! Pressure at full levels |
|---|
| 59 | REAL, intent(in):: pn2(ngrid, nlayer) ! N2 (BV^2) at 1/2 levels |
|---|
| 60 | REAL, intent(in):: presnivs(nlayer) ! Approximate pressure for reproductible launching altitude |
|---|
| 61 | REAL, intent(in):: pt(ngrid, nlayer) ! Temperature at full levels |
|---|
| 62 | REAL, intent(in):: pu(ngrid, nlayer) ! Zonal wind at full levels |
|---|
| 63 | REAL, intent(in):: pv(ngrid, nlayer) ! Meridional wind at full levels |
|---|
| 64 | |
|---|
| 65 | ! 0.2 OUTPUTS |
|---|
| 66 | |
|---|
| 67 | REAL, intent(out):: zustr(ngrid) ! Zonal surface stress |
|---|
| 68 | REAL, intent(out):: zvstr(ngrid) ! Meridional surface stress |
|---|
| 69 | REAL, intent(out):: d_t(ngrid, nlayer) ! Tendency on temperature |
|---|
| 70 | REAL, intent(out):: d_u(ngrid, nlayer) ! Tendency on zonal wind |
|---|
| 71 | REAL, intent(out):: d_v(ngrid, nlayer) ! Tendency on meridional wind |
|---|
| 72 | |
|---|
| 73 | ! 0.3 INTERNAL ARRAYS |
|---|
| 74 | |
|---|
| 75 | REAL :: tt(ngrid, nlayer) ! Temperature at full levels |
|---|
| 76 | REAL :: uu(ngrid, nlayer) ! Zonal wind at full levels |
|---|
| 77 | REAL :: vv(ngrid, nlayer) ! Meridional wind at full levels |
|---|
| 78 | |
|---|
| 79 | INTEGER ii, jj, ll |
|---|
| 80 | |
|---|
| 81 | ! 0.3.0 Time scale of the like cycle of the waves parametrized |
|---|
| 82 | REAL, parameter:: deltat = 24. * 3600. |
|---|
| 83 | |
|---|
| 84 | ! 0.3.1 Gravity waves specifications |
|---|
| 85 | INTEGER, parameter:: nk = 2 ! number of horizontal wavenumbers |
|---|
| 86 | INTEGER, parameter:: np = 2 ! directions (eastward and westward) phase speed |
|---|
| 87 | INTEGER, parameter:: no = 2 ! absolute values of phase speed |
|---|
| 88 | INTEGER, parameter:: na = 5 ! Number of realizations to get the phase speed |
|---|
| 89 | INTEGER, parameter:: nw = nk * np *no ! Total numbers of gravity waves |
|---|
| 90 | INTEGER jk, jp, jo, jw |
|---|
| 91 | REAL kstar ! Control value to constrain the min horizontal wavenumber by the horizontal resolution |
|---|
| 92 | |
|---|
| 93 | !!! TEST Kobs and BestFit Diogo |
|---|
| 94 | REAL, parameter:: kmax = 1.e-4 ! Max horizontal wavenumber (lambda min 50 km) |
|---|
| 95 | ! generic : kmax=N/u, u(=30) zonal wind velocity at launch |
|---|
| 96 | REAL, parameter:: kmin = 1.e-5 ! Min horizontal wavenumber (lambda max 628 km) |
|---|
| 97 | ! generic : kmin=1/sqrt(DxDy) Dx and Dy horizontal grid |
|---|
| 98 | |
|---|
| 99 | !---------------------------------------- |
|---|
| 100 | ! VCD 1.1 tuning |
|---|
| 101 | ! REAL, parameter:: cmax = 61. ! Max horizontal absolute phase velocity |
|---|
| 102 | !---------------------------------------- |
|---|
| 103 | ! VCD 2.0 tuning |
|---|
| 104 | REAL, parameter:: cmax = 61. ! Max horizontal absolute phase velocity |
|---|
| 105 | !---------------------------------------- |
|---|
| 106 | |
|---|
| 107 | ! generic : cmax=zonal wind at launch |
|---|
| 108 | REAL, parameter:: cmin = 1. ! Min horizontal absolute phase velocity |
|---|
| 109 | |
|---|
| 110 | !---------------------------------------- |
|---|
| 111 | ! Nominal as Gilli2017 |
|---|
| 112 | ! REAL, parameter:: epflux_max = 0.005 ! Max EP flux value at launching altitude (previous name: RUWMAX) |
|---|
| 113 | !---------------------------------------- |
|---|
| 114 | ! VCD 2.1 tuning |
|---|
| 115 | REAL, parameter:: epflux_max = 0.001 ! Max EP flux value at launching altitude (previous name: RUWMAX) |
|---|
| 116 | !---------------------------------------- |
|---|
| 117 | |
|---|
| 118 | REAL cpha ! absolute phase velocity frequency |
|---|
| 119 | REAL zk(nw, ngrid) ! horizontal wavenumber amplitude |
|---|
| 120 | REAL zp(nw, ngrid) ! horizontal wavenumber angle |
|---|
| 121 | REAL zo(nw, ngrid) ! absolute frequency |
|---|
| 122 | |
|---|
| 123 | REAL intr_freq_m(nw, ngrid) ! Waves Intr. freq. at the 1/2 lev below the full level (previous name: ZOM) |
|---|
| 124 | REAL intr_freq_p(nw, ngrid) ! Waves Intr. freq. at the 1/2 lev above the full level (previous name: ZOP) |
|---|
| 125 | REAL wwm(nw, ngrid) ! Wave EP-fluxes at the 1/2 level below the full level |
|---|
| 126 | REAL wwp(nw, ngrid) ! Wave EP-fluxes at the 1/2 level above the full level |
|---|
| 127 | REAL u_epflux_p(nw, ngrid) ! Partial zonal flux (=for each wave) at the 1/2 level above the full level (previous name: RUWP) |
|---|
| 128 | REAL v_epflux_p(nw, ngrid) ! Partial meridional flux (=for each wave) at the 1/2 level above the full level (previous name: RVWP) |
|---|
| 129 | 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) |
|---|
| 130 | 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) |
|---|
| 131 | REAL epflux_0(nw, ngrid) ! Fluxes at launching level (previous name: RUW0) |
|---|
| 132 | INTEGER launch ! Launching altitude |
|---|
| 133 | REAL, parameter:: xlaunch = 5e-3 ! Value for top of cloud convective region |
|---|
| 134 | |
|---|
| 135 | ! 0.3.2 Parameters of waves dissipations |
|---|
| 136 | !---------------------------------------- |
|---|
| 137 | ! VCD 1.1 tuning |
|---|
| 138 | ! REAL, parameter:: sat = 0.85 ! saturation parameter |
|---|
| 139 | ! REAL, parameter:: rdiss = 0.1 ! coefficient of dissipation |
|---|
| 140 | !---------------------------------------- |
|---|
| 141 | ! VCD 2.0 tuning |
|---|
| 142 | REAL, parameter:: sat = 0.6 ! saturation parameter |
|---|
| 143 | REAL, parameter:: rdiss = 8.e-4 ! coefficient of dissipation |
|---|
| 144 | !---------------------------------------- |
|---|
| 145 | REAL, parameter:: zoisec = 1.e-8 ! security for intrinsic freguency |
|---|
| 146 | |
|---|
| 147 | ! 0.3.3 Background flow at 1/2 levels and vertical coordinate |
|---|
| 148 | REAL H0bis(ngrid, nlayer) ! characteristic height of the atmosphere |
|---|
| 149 | REAL min_k(ngrid) ! min(kstar,kmin) |
|---|
| 150 | REAL, save:: H0 ! characteristic height of the atmosphere |
|---|
| 151 | REAL, save:: MDR ! characteristic mass density |
|---|
| 152 | !---------------------------------------- |
|---|
| 153 | ! VCD 1.1 tuning |
|---|
| 154 | ! REAL, parameter:: pr = 5e5 ! Reference pressure [Pa] ! VENUS: cloud layer |
|---|
| 155 | !---------------------------------------- |
|---|
| 156 | ! VCD 2.0 tuning |
|---|
| 157 | REAL, parameter:: pr = 5e4 ! Reference pressure [Pa] ! VENUS: cloud layer |
|---|
| 158 | !---------------------------------------- |
|---|
| 159 | REAL, parameter:: tr = 300. ! Reference temperature [K] ! VENUS: cloud layer |
|---|
| 160 | REAL zh(ngrid, nlayer + 1) ! Log-pressure altitude (constant H0) |
|---|
| 161 | REAL zhbis(ngrid, nlayer + 1) ! Log-pressure altitude (varying H) |
|---|
| 162 | REAL uh(ngrid, nlayer + 1) ! Zonal wind at 1/2 levels |
|---|
| 163 | REAL vh(ngrid, nlayer + 1) ! Meridional wind at 1/2 levels |
|---|
| 164 | REAL ph(ngrid, nlayer + 1) ! Pressure at 1/2 levels |
|---|
| 165 | REAL, parameter:: psec = 1.e-9 ! Security to avoid division by 0 pressure |
|---|
| 166 | REAL bv(ngrid, nlayer + 1) ! Brunt Vaisala freq. at 1/2 levels |
|---|
| 167 | REAL, parameter:: bvsec = 1.e-5 ! Security to avoid negative BV |
|---|
| 168 | REAL href(nlayer + 1) ! Reference altitude for launching altitude |
|---|
| 169 | |
|---|
| 170 | REAL ran_num_1, ran_num_2, ran_num_3 |
|---|
| 171 | |
|---|
| 172 | LOGICAL, save :: firstcall = .true. |
|---|
| 173 | |
|---|
| 174 | !-------------------------------------------------------------------------------------------------- |
|---|
| 175 | ! 1. INITIALISATION |
|---|
| 176 | !------------------ |
|---|
| 177 | IF (firstcall) THEN |
|---|
| 178 | write(*,*) "nonoro_gwd_ran: FLott non-oro GW scheme is active!" |
|---|
| 179 | ! Characteristic vertical scale height |
|---|
| 180 | H0 = RD * tr / RG |
|---|
| 181 | ! Characteristic mass density at launch altitude |
|---|
| 182 | MDR = pr / ( RD * tr ) |
|---|
| 183 | ! Control |
|---|
| 184 | if (deltat .LT. dtime) THEN |
|---|
| 185 | call abort_physic("nonoro_gwd_ran","gwd random: deltat lower than dtime!",1) |
|---|
| 186 | endif |
|---|
| 187 | if (nlayer .LT. nw) THEN |
|---|
| 188 | call abort_physic("nonoro_gwd_ran","gwd random: nlayer lower than nw!",1) |
|---|
| 189 | endif |
|---|
| 190 | firstcall = .false. |
|---|
| 191 | ENDIF |
|---|
| 192 | |
|---|
| 193 | ! For VENUS, we use *_seri variables, that already integrate the different previous tendencies |
|---|
| 194 | tt(:,:) = pt(:,:) |
|---|
| 195 | uu(:,:) = pu(:,:) |
|---|
| 196 | vv(:,:) = pv(:,:) |
|---|
| 197 | |
|---|
| 198 | ! print*,'epflux_max just after firstcall:', epflux_max |
|---|
| 199 | |
|---|
| 200 | ! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS |
|---|
| 201 | !---------------------------------------------------- |
|---|
| 202 | ! Pressure and Inv of pressure, temperature at 1/2 level |
|---|
| 203 | DO ll = 2, nlayer |
|---|
| 204 | ph(:, ll) = exp((log(pp(:, ll)) + log(pp(:, ll - 1))) / 2.) |
|---|
| 205 | end DO |
|---|
| 206 | |
|---|
| 207 | ph(:, nlayer + 1) = 0. |
|---|
| 208 | ph(:, 1) = 2. * pp(:, 1) - ph(:, 2) |
|---|
| 209 | |
|---|
| 210 | ! Launching altitude for reproductible case |
|---|
| 211 | DO ll = 2, nlayer |
|---|
| 212 | href(ll) = exp((log(presnivs(ll)) + log(presnivs(ll - 1))) / 2.) |
|---|
| 213 | end DO |
|---|
| 214 | href(nlayer + 1) = 0. |
|---|
| 215 | href(1) = 2. * presnivs(1) - href(2) |
|---|
| 216 | |
|---|
| 217 | launch = 0. |
|---|
| 218 | DO ll =1, nlayer |
|---|
| 219 | IF (href (ll) / href(1) > xlaunch) launch = ll |
|---|
| 220 | end DO |
|---|
| 221 | |
|---|
| 222 | ! Log-pressure vertical coordinate |
|---|
| 223 | DO ll = 1, nlayer + 1 |
|---|
| 224 | zh(:,ll) = H0 * log(pr / (ph(:,ll) + psec)) |
|---|
| 225 | end DO |
|---|
| 226 | |
|---|
| 227 | ! Winds and Brunt Vaisala frequency |
|---|
| 228 | DO ll = 2, nlayer |
|---|
| 229 | uh(:, ll) = 0.5 * (uu(:, ll) + uu(:, ll - 1)) ! Zonal wind |
|---|
| 230 | vh(:, ll) = 0.5 * (vv(:, ll) + vv(:, ll - 1)) ! Meridional wind |
|---|
| 231 | ! VENUS ATTENTION: CP VARIABLE PSTAB CALCULE EN AMONT DES PARAMETRISATIONS |
|---|
| 232 | bv(:, ll) = max(bvsec,sqrt(pn2(:,ll))) |
|---|
| 233 | end DO |
|---|
| 234 | |
|---|
| 235 | bv(:, 1) = bv(:, 2) |
|---|
| 236 | uh(:, 1) = 0. |
|---|
| 237 | vh(:, 1) = 0. |
|---|
| 238 | bv(:, nlayer + 1) = bv(:, nlayer) |
|---|
| 239 | uh(:, nlayer + 1) = uu(:, nlayer) |
|---|
| 240 | vh(:, nlayer + 1) = vv(:, nlayer) |
|---|
| 241 | |
|---|
| 242 | ! TN+GG April/June 2020 |
|---|
| 243 | ! "Individual waves are not supposed |
|---|
| 244 | ! to occupy the entire domain, but |
|---|
| 245 | ! only a faction of it" Lott+2012 |
|---|
| 246 | ! minimum value of k between kmin and kstar |
|---|
| 247 | DO ii = 1, ngrid |
|---|
| 248 | kstar = RPI / sqrt(cell_area(ii)) |
|---|
| 249 | min_k(ii) = max(kmin,kstar) |
|---|
| 250 | end DO |
|---|
| 251 | |
|---|
| 252 | ! 3. WAVES CHARACTERISTICS CHOSEN RANDOMLY |
|---|
| 253 | !----------------------------------------- |
|---|
| 254 | ! The mod function of here a weird arguments |
|---|
| 255 | ! are used to produce the waves characteristics |
|---|
| 256 | ! in a stochastic way |
|---|
| 257 | DO jw = 1, nw |
|---|
| 258 | DO ii = 1, ngrid |
|---|
| 259 | |
|---|
| 260 | ran_num_1 = mod(tt(ii, jw) * 10., 1.) |
|---|
| 261 | ran_num_2 = mod(tt(ii, jw) * 100., 1.) |
|---|
| 262 | |
|---|
| 263 | !! OPTIONS GENERIC DIFF VENUS !! |
|---|
| 264 | ! angle (random) - reference |
|---|
| 265 | zp(jw, ii) = (sign(1., 0.5 - ran_num_1) & |
|---|
| 266 | + 1.) * RPI / 2. |
|---|
| 267 | ! zp(jw, ii) = atan2(4.*vh(ii, launch),uh(ii, launch)) & |
|---|
| 268 | ! + (sign(1., 0.5 - ran_num_1) + 1.) * RPI / 2. |
|---|
| 269 | ! --------- TRY 00 ----------------------- |
|---|
| 270 | !IF((abs(latitude_deg(ii)) .le. 55.)) THEN |
|---|
| 271 | ! zp(jw, ii) = 0. + 2. * ( ran_num_1 - 0.5 ) * RPI /180. & |
|---|
| 272 | ! * (10.) ! +/- 10 deg par rapport equateur |
|---|
| 273 | !ELSE IF ((abs(latitude_deg(ii)) .le. 75.)) THEN |
|---|
| 274 | ! zp(jw, ii) = 0. + 2. * ( ran_num_1 - 0.5 ) * RPI /180. & |
|---|
| 275 | ! * (10. +(90.-10.) & |
|---|
| 276 | ! * (latitude_deg(ii)-55.)/20.) |
|---|
| 277 | !ELSE |
|---|
| 278 | ! zp(jw, ii) = ran_num_1 * RPI |
|---|
| 279 | !ENDIF |
|---|
| 280 | !zp(jw, ii) = mod(zp(jw, ii),2.*RPI) |
|---|
| 281 | ! ------ TRY 01------------------------------- |
|---|
| 282 | !IF((abs(latitude_deg(ii)) .le. 55)) THEN |
|---|
| 283 | ! zp(jw, ii) = (sign(1., 0.5 - ran_num_1) & |
|---|
| 284 | ! + 1.) * RPI / 2. |
|---|
| 285 | !ELSE |
|---|
| 286 | ! zp(jw, ii) = ran_num_1 * RPI |
|---|
| 287 | !ENDIF |
|---|
| 288 | ! ---- angle (0 or RPI) ----- |
|---|
| 289 | !!zp(jw, ii) = RPI*mod(jw,2) |
|---|
| 290 | |
|---|
| 291 | ! horizontal wavenumber amplitude |
|---|
| 292 | ! zk(jw, ii) = kmin + (kmax - kmin) * ran_num_2 |
|---|
| 293 | zk(jw, ii) = min_k(ii) + (kmax - min_k(ii)) * ran_num_2 |
|---|
| 294 | |
|---|
| 295 | ! horizontal phase speed |
|---|
| 296 | ! this computation allows to get a gaussian distribution centered on 0 (right ?) |
|---|
| 297 | ! then cmin is not useful, and it favors small values. |
|---|
| 298 | ! VCD 2.0 tuning |
|---|
| 299 | cpha = 0. |
|---|
| 300 | DO jj = 1, na |
|---|
| 301 | ran_num_3 = mod(tt(ii, jw + 3 * jj)**2, 1.) |
|---|
| 302 | cpha = cpha + 2. * cmax * & |
|---|
| 303 | (ran_num_3 - 0.5) * & |
|---|
| 304 | sqrt(3.) / sqrt(na * 1.) |
|---|
| 305 | end DO |
|---|
| 306 | !cpha = cpha - uh(ii, launch) |
|---|
| 307 | IF (cpha < 0.) THEN |
|---|
| 308 | cpha = - 1. * cpha |
|---|
| 309 | zp (jw, ii) = zp(jw, ii) + RPI |
|---|
| 310 | ENDIF |
|---|
| 311 | !IF (cpha < 1.) THEN |
|---|
| 312 | ! cpha = 1. |
|---|
| 313 | !ENDIF |
|---|
| 314 | ! otherwise, with the computation below, we get a uniform distribution between cmin and cmax. |
|---|
| 315 | ! ran_num_3 = mod(tt(ii, jw)**2, 1.) |
|---|
| 316 | ! cpha = cmin + (cmax - cmin) * ran_num_3 |
|---|
| 317 | |
|---|
| 318 | ! Intrinsic frequency |
|---|
| 319 | zo(jw, ii) = cpha * zk(jw, ii) |
|---|
| 320 | ! Intrinsic frequency is imposed |
|---|
| 321 | zo(jw, ii) = zo(jw, ii) & |
|---|
| 322 | + zk(jw, ii) * cos(zp(jw, ii)) * uh(ii, launch) & |
|---|
| 323 | + zk(jw, ii) * sin(zp(jw, ii)) * vh(ii, launch) |
|---|
| 324 | |
|---|
| 325 | ! Momentum flux at launch level |
|---|
| 326 | epflux_0(jw, ii) = epflux_max & |
|---|
| 327 | * mod(100. * (uu(ii, jw)**2 + vv(ii, jw)**2), 1.) |
|---|
| 328 | |
|---|
| 329 | end DO |
|---|
| 330 | end DO |
|---|
| 331 | ! print*,'epflux_0 just after waves charac. ramdon:', epflux_0 |
|---|
| 332 | |
|---|
| 333 | ! 4. COMPUTE THE FLUXES |
|---|
| 334 | !---------------------- |
|---|
| 335 | ! 4.1 Vertical velocity at launching altitude to ensure |
|---|
| 336 | ! the correct value to the imposed fluxes. |
|---|
| 337 | !------------------------------------------------------ |
|---|
| 338 | DO jw = 1, nw |
|---|
| 339 | ! Evaluate intrinsic frequency at launching altutide: |
|---|
| 340 | intr_freq_p(jw, :) = zo(jw, :) & |
|---|
| 341 | - zk(jw, :) * cos(zp(jw, :)) * uh(:, launch) & |
|---|
| 342 | - zk(jw, :) * sin(zp(jw, :)) * vh(:, launch) |
|---|
| 343 | end DO |
|---|
| 344 | |
|---|
| 345 | ! Vertical velocity at launch level, value to ensure the imposed |
|---|
| 346 | ! mom flux: |
|---|
| 347 | DO jw = 1, nw |
|---|
| 348 | ! WW is directly a flux, here, not vertical velocity anymore |
|---|
| 349 | wwp(jw, :) = epflux_0(jw,:) |
|---|
| 350 | u_epflux_p(jw, :) = cos(zp(jw, :)) * sign(1., intr_freq_p(jw, :)) * epflux_0(jw, :) |
|---|
| 351 | v_epflux_p(jw, :) = sin(zp(jw, :)) * sign(1., intr_freq_p(jw, :)) * epflux_0(jw, :) |
|---|
| 352 | end DO |
|---|
| 353 | |
|---|
| 354 | ! 4.2 Initial flux at launching altitude |
|---|
| 355 | !--------------------------------------- |
|---|
| 356 | u_epflux_tot(:, launch) = 0. |
|---|
| 357 | v_epflux_tot(:, launch) = 0. |
|---|
| 358 | DO jw = 1, nw |
|---|
| 359 | u_epflux_tot(:, launch) = u_epflux_tot(:, launch) + u_epflux_p(jw, :) |
|---|
| 360 | v_epflux_tot(:, launch) = v_epflux_tot(:, launch) + v_epflux_p(jw, :) |
|---|
| 361 | end DO |
|---|
| 362 | |
|---|
| 363 | ! 4.3 Loop over altitudes, with passage from one level to the next done by: |
|---|
| 364 | !-------------------------------------------------------------------------- |
|---|
| 365 | ! i) conserving the EP flux, |
|---|
| 366 | ! ii) dissipating a little, |
|---|
| 367 | ! iii) testing critical levels, |
|---|
| 368 | ! iv) testing the breaking. |
|---|
| 369 | !-------------------------------------------------------------------------- |
|---|
| 370 | DO ll = launch, nlayer - 1 |
|---|
| 371 | ! Warning! all the physics is here (passage from one level to the next |
|---|
| 372 | DO jw = 1, nw |
|---|
| 373 | intr_freq_m(jw, :) = intr_freq_p(jw, :) |
|---|
| 374 | wwm(jw, :) = wwp(jw, :) |
|---|
| 375 | ! Intrinsic frequency |
|---|
| 376 | intr_freq_p(jw, :) = zo(jw, :) - zk(jw, :) * cos(zp(jw, :)) * uh(:, ll + 1) & |
|---|
| 377 | - zk(jw, :) * sin(zp(jw, :)) * vh(:, ll + 1) |
|---|
| 378 | ! Vertical velocity in flux formulation |
|---|
| 379 | wwp(jw, :) = min( & |
|---|
| 380 | ! No breaking (Eq. 6): |
|---|
| 381 | wwm(jw, :) & |
|---|
| 382 | ! Dissipation (Eq. 8): |
|---|
| 383 | * exp(-rdiss * pr / (ph(:, ll + 1) + ph (:, ll)) & |
|---|
| 384 | * ((bv(:, ll + 1) + bv (:, ll)) / 2.)**3 & |
|---|
| 385 | / max(abs(intr_freq_p(jw, :) + intr_freq_m(jw, :)) / 2., zoisec)**4 & |
|---|
| 386 | * zk(jw, :)**3 * (zh(:, ll + 1) - zh(:, ll))) , & |
|---|
| 387 | ! Critical levels (forced to zero if intrinsic frequency |
|---|
| 388 | ! changes sign) |
|---|
| 389 | max(0., sign(1., intr_freq_p(jw, :) * intr_freq_m(jw, :))) & |
|---|
| 390 | ! Saturation (Eq. 12) |
|---|
| 391 | * abs(intr_freq_p(jw, :))**3 / bv(:, ll + 1) & |
|---|
| 392 | * exp(-zh(:, ll + 1) / H0) & |
|---|
| 393 | !! Correction for VCD 2.0 |
|---|
| 394 | ! * sat**2 * kmin**2 / zk(jw, :)**4) |
|---|
| 395 | * sat**2 * min_k(:)**2 * MDR / zk(jw, :)**4) |
|---|
| 396 | end DO |
|---|
| 397 | |
|---|
| 398 | ! Evaluate EP-flux from Eq. 7 and give the right orientation to the stress |
|---|
| 399 | DO jw = 1, nw |
|---|
| 400 | u_epflux_p(jw, :) = sign(1., intr_freq_p(jw, :)) * cos(zp(jw, :)) * wwp(jw, :) |
|---|
| 401 | v_epflux_p(jw, :) = sign(1., intr_freq_p(jw, :)) * sin(zp(jw, :)) * wwp(jw, :) |
|---|
| 402 | end DO |
|---|
| 403 | |
|---|
| 404 | u_epflux_tot(:, ll + 1) = 0. |
|---|
| 405 | v_epflux_tot(:, ll + 1) = 0. |
|---|
| 406 | DO jw = 1, nw |
|---|
| 407 | u_epflux_tot(:, ll + 1) = u_epflux_tot(:, ll + 1) + u_epflux_p(jw, :) |
|---|
| 408 | v_epflux_tot(:, ll + 1) = v_epflux_tot(:, ll + 1) + v_epflux_p(jw, :) |
|---|
| 409 | east_gwstress(:, ll) = east_gwstress(:, ll) & |
|---|
| 410 | + max(0., u_epflux_p(jw, :)) / float(nw) |
|---|
| 411 | west_gwstress(:, ll) = west_gwstress(:, ll) & |
|---|
| 412 | + min(0., u_epflux_p(jw, ll)) / float(nw) |
|---|
| 413 | end DO |
|---|
| 414 | end DO ! DO LL = LAUNCH, nlayer - 1 |
|---|
| 415 | |
|---|
| 416 | ! print*,'u_epflux_tot just after launching:', u_epflux_tot |
|---|
| 417 | ! print*,'v_epflux_tot just after launching:', v_epflux_tot |
|---|
| 418 | ! print*,'u_epflux_p just after launching:', maxval(u_epflux_p), minval(u_epflux_p) |
|---|
| 419 | ! print*,'v_epflux_p just after launching:', maxval(v_epflux_p), minval(v_epflux_p) |
|---|
| 420 | |
|---|
| 421 | ! 5. TENDENCY CALCULATIONS |
|---|
| 422 | !------------------------- |
|---|
| 423 | ! 5.1 Flow rectification at the top and in the low layers |
|---|
| 424 | ! -------------------------------------------------------- |
|---|
| 425 | ! Warning, this is the total on all GW... |
|---|
| 426 | |
|---|
| 427 | u_epflux_tot(:, nlayer + 1) = 0. |
|---|
| 428 | v_epflux_tot(:, nlayer + 1) = 0. |
|---|
| 429 | |
|---|
| 430 | ! Here, big change compared to FLott version: |
|---|
| 431 | ! We compensate (u_epflux_(:, LAUNCH), ie total emitted upward flux |
|---|
| 432 | ! over the layers max(1,LAUNCH-3) to LAUNCH-1 |
|---|
| 433 | DO LL = 1, max(1,LAUNCH-3) |
|---|
| 434 | u_epflux_tot(:, LL) = 0. |
|---|
| 435 | v_epflux_tot(:, LL) = 0. |
|---|
| 436 | end DO |
|---|
| 437 | DO LL = max(2,LAUNCH-2), LAUNCH-1 |
|---|
| 438 | u_epflux_tot(:, LL) = u_epflux_tot(:, LL - 1) & |
|---|
| 439 | + u_epflux_tot(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1)) & |
|---|
| 440 | / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3))) |
|---|
| 441 | v_epflux_tot(:, LL) = v_epflux_tot(:, LL - 1) & |
|---|
| 442 | + v_epflux_tot(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1)) & |
|---|
| 443 | / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3))) |
|---|
| 444 | east_gwstress(:,LL) = east_gwstress(:, LL - 1) & |
|---|
| 445 | + east_gwstress(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1)) & |
|---|
| 446 | / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3))) |
|---|
| 447 | west_gwstress(:,LL) = west_gwstress(:, LL - 1) & |
|---|
| 448 | + west_gwstress(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1)) & |
|---|
| 449 | / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3))) |
|---|
| 450 | end DO |
|---|
| 451 | ! This way, the total flux from GW is zero, but there is a net transport |
|---|
| 452 | ! (upward) that should be compensated by circulation |
|---|
| 453 | ! and induce additional friction at the surface |
|---|
| 454 | |
|---|
| 455 | ! 5.2 AR-1 RECURSIVE FORMULA (13) IN VERSION 4 |
|---|
| 456 | !--------------------------------------------- |
|---|
| 457 | DO LL = 1, nlayer |
|---|
| 458 | ! d_u(:, LL) = (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL)) |
|---|
| 459 | d_u(:, LL) = RG * (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL)) & |
|---|
| 460 | / (PH(:, LL + 1) - PH(:, LL)) * DTIME |
|---|
| 461 | ! d_v(:, LL) = (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL)) |
|---|
| 462 | d_v(:, LL) = RG * (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL)) & |
|---|
| 463 | / (PH(:, LL + 1) - PH(:, LL)) * DTIME |
|---|
| 464 | ENDDO |
|---|
| 465 | d_t(:,:) = 0. |
|---|
| 466 | |
|---|
| 467 | ! 5.3 Update tendency of wind with the previous (and saved) values |
|---|
| 468 | !----------------------------------------------------------------- |
|---|
| 469 | d_u(:,:) = DTIME/DELTAT/REAL(NW) * d_u(:,:) & |
|---|
| 470 | + (1.-DTIME/DELTAT) * du_nonoro_gwd(:,:) |
|---|
| 471 | d_v(:,:) = DTIME/DELTAT/REAL(NW) * d_v(:,:) & |
|---|
| 472 | + (1.-DTIME/DELTAT) * dv_nonoro_gwd(:,:) |
|---|
| 473 | du_nonoro_gwd(:,:) = d_u(:,:) |
|---|
| 474 | dv_nonoro_gwd(:,:) = d_v(:,:) |
|---|
| 475 | |
|---|
| 476 | ! print*,'u_epflux_tot just after tendency:', u_epflux_tot |
|---|
| 477 | ! print*,'v_epflux_tot just after tendency:', v_epflux_tot |
|---|
| 478 | ! print*,'d_u just after tendency:', maxval(d_u(:,:)), minval(d_u) |
|---|
| 479 | ! print*,'d_v just after tendency:', maxval(d_v(:,:)), minval(d_v) |
|---|
| 480 | |
|---|
| 481 | ! Cosmetic: evaluation of the cumulated stress |
|---|
| 482 | ZUSTR(:) = 0. |
|---|
| 483 | ZVSTR(:) = 0. |
|---|
| 484 | DO LL = 1, nlayer |
|---|
| 485 | ZUSTR(:) = ZUSTR(:) + D_U(:, LL) / RG * (PH(:, LL + 1) - PH(:, LL)) |
|---|
| 486 | ZVSTR(:) = ZVSTR(:) + D_V(:, LL) / RG * (PH(:, LL + 1) - PH(:, LL)) |
|---|
| 487 | ENDDO |
|---|
| 488 | |
|---|
| 489 | !#ifdef CPP_XIOS |
|---|
| 490 | ! call send_xios_field("du_nonoro", d_u) |
|---|
| 491 | ! call send_xios_field("dv_nonoro", d_v) |
|---|
| 492 | ! call send_xios_field("east_gwstress", east_gwstress) |
|---|
| 493 | ! call send_xios_field("west_gwstress", west_gwstress) |
|---|
| 494 | !#endif |
|---|
| 495 | |
|---|
| 496 | END SUBROUTINE nonoro_gwd_ran |
|---|
| 497 | |
|---|
| 498 | ! =================================================================== |
|---|
| 499 | ! Subroutines used to write variables of memory in start files |
|---|
| 500 | ! =================================================================== |
|---|
| 501 | |
|---|
| 502 | SUBROUTINE ini_nonoro_gwd_ran(ngrid,nlayer) |
|---|
| 503 | |
|---|
| 504 | IMPLICIT NONE |
|---|
| 505 | |
|---|
| 506 | INTEGER, INTENT (in) :: ngrid ! number of atmospheric columns |
|---|
| 507 | INTEGER, INTENT (in) :: nlayer ! number of atmospheric layers |
|---|
| 508 | |
|---|
| 509 | allocate(du_nonoro_gwd(ngrid, nlayer)) |
|---|
| 510 | allocate(dv_nonoro_gwd(ngrid, nlayer)) |
|---|
| 511 | allocate(east_gwstress(ngrid, nlayer)) |
|---|
| 512 | allocate(west_gwstress(ngrid, nlayer)) |
|---|
| 513 | |
|---|
| 514 | END SUBROUTINE ini_nonoro_gwd_ran |
|---|
| 515 | ! ---------------------------------- |
|---|
| 516 | SUBROUTINE end_nonoro_gwd_ran |
|---|
| 517 | |
|---|
| 518 | IMPLICIT NONE |
|---|
| 519 | |
|---|
| 520 | if (allocated(du_nonoro_gwd)) deallocate(du_nonoro_gwd) |
|---|
| 521 | if (allocated(dv_nonoro_gwd)) deallocate(dv_nonoro_gwd) |
|---|
| 522 | if (allocated(east_gwstress)) deallocate(east_gwstress) |
|---|
| 523 | if (allocated(west_gwstress)) deallocate(west_gwstress) |
|---|
| 524 | |
|---|
| 525 | END SUBROUTINE end_nonoro_gwd_ran |
|---|
| 526 | |
|---|
| 527 | END MODULE nonoro_gwd_ran_mod |
|---|
| 528 | |
|---|
| 529 | |
|---|
