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