Changeset 2595
- Timestamp:
- Dec 16, 2021, 3:09:04 PM (3 years ago)
- Location:
- trunk/LMDZ.GENERIC
- Files:
-
- 5 edited
-
README (modified) (1 diff)
-
libf/phystd/callkeys_mod.F90 (modified) (2 diffs)
-
libf/phystd/inifis_mod.F90 (modified) (1 diff)
-
libf/phystd/nonoro_gwd_ran_mod.F90 (modified) (17 diffs)
-
libf/phystd/physiq_mod.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/README
r2585 r2595 1700 1700 == 17/11/2021 == YJ 1701 1701 Correct missing initialisation for k-coefficients mixing on the fly 1702 1703 == 16/12/2021 == JL 1704 Fixes and improvements in the Non-orographic GW scheme, namely: 1705 - increments are not tendencies 1706 - missing rho factor in EP flux computation 1707 - missing rho at launch altitude 1708 - changed inputs, because R and Cp are needed to compute rho and BV -
trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90
r2542 r2595 64 64 !$OMP THREADPRIVATE(photochem,photoheat,jonline,depos) 65 65 logical,save :: calllott_nonoro 66 logical,save :: gwd_convective_source 67 !$OMP THREADPRIVATE(calllott_nonoro,gwd_convective_source) 66 !$OMP THREADPRIVATE(calllott_nonoro) 68 67 logical,save :: global1d 69 68 real,save :: szangle … … 145 144 real,save :: surfemis 146 145 !$OMP THREADPRIVATE(surfalbedo,surfemis) 147 real,save :: epflux_max148 !$OMP THREADPRIVATE(epflux_max)149 146 real,save :: noseason_day 150 147 !$OMP THREADPRIVATE(noseason_day) -
trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90
r2583 r2595 329 329 ": calllott_nonoro = ",calllott_nonoro 330 330 331 if (is_master) write(*,*)trim(rname)//&332 ": Max Value of EP flux at launching altitude"333 epflux_max=0.0 ! default value334 call getin_p("epflux_max",epflux_max)335 if (is_master) write(*,*)trim(rname)//": epflux_max = ",epflux_max336 337 331 if (is_master) write(*,*) trim(rname)//": call thermal plume model ?" 338 332 calltherm=.false. ! default value -
trunk/LMDZ.GENERIC/libf/phystd/nonoro_gwd_ran_mod.F90
r2550 r2595 10 10 CONTAINS 11 11 12 SUBROUTINE nonoro_gwd_ran(ngrid, nlayer, dtime, pp, &12 SUBROUTINE nonoro_gwd_ran(ngrid, nlayer, dtime, cpnew, rnew, pp, & 13 13 zmax_therm, pt, pu, pv, pdt, pdu, pdv, & 14 14 zustr, zvstr, d_t, d_u, d_v) … … 38 38 ! 39 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. 40 45 !--------------------------------------------------------------------------------- 41 46 42 47 43 48 44 use comcstfi_mod, only: g, pi, cpp,r49 use comcstfi_mod, only: g, pi, r 45 50 USE ioipsl_getin_p_mod, ONLY : getin_p 46 use assert_m, only : assert47 51 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_source50 use inifis_mod51 52 use geometry_mod, only: cell_area 52 53 #ifdef CPP_XIOS 53 54 use xios_output_mod, only: send_xios_field 54 55 #endif 55 implicit none56 57 ! 0. DECLARATIONS:56 implicit none 57 58 ! 0. DECLARATIONS: 58 59 59 60 ! 0.1 INPUTS 60 61 61 INTEGER, intent(in):: ngrid, nlayer 62 INTEGER, intent(in):: ngrid ! number of atmospheric columns 63 INTEGER, intent(in):: nlayer ! number of atmospheric columns 62 64 REAL, intent(in):: dtime ! Time step of the Physics 63 65 REAL, intent(in):: zmax_therm(ngrid) ! Altitude of max velocity thermals (m) 64 REAL, intent(in):: pp(ngrid, nlayer) ! Pressure at full levels 65 REAL, intent(in):: pt(ngrid, nlayer) ! Temperature at full levels 66 REAL, intent(in):: pu(ngrid, nlayer) ! Zonal wind at full levels 67 REAL, intent(in):: pv(ngrid, nlayer) ! Meridional wind at full levels 68 REAL, intent(in):: pdt(ngrid, nlayer) ! Tendency on temperature 69 REAL, intent(in):: pdu(ngrid, nlayer) ! Tendency on zonal wind 70 REAL, intent(in):: pdv(ngrid, nlayer) ! Tendency on meridional wind 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) 71 75 72 76 ! 0.2 OUTPUTS … … 74 78 REAL, intent(out):: zustr(ngrid) ! Zonal surface stress 75 79 REAL, intent(out):: zvstr(ngrid) ! Meridional surface stress 76 REAL, intent(out):: d_t(ngrid, nlayer) ! Tendency on temperature 77 REAL, intent(out):: d_u(ngrid, nlayer) ! Tendency on zonal wind 78 REAL, intent(out):: d_v(ngrid, nlayer) ! Tendency on meridional wind 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 79 83 80 84 ! 0.3 INTERNAL ARRAYS 81 85 82 86 REAL :: tt(ngrid, nlayer) ! Temperature at full levels 87 REAL :: rho(ngrid, nlayer) ! Mass density at full levels 83 88 REAL :: uu(ngrid, nlayer) ! Zonal wind at full levels 84 89 REAL :: vv(ngrid, nlayer) ! Meridional wind at full levels 85 REAL :: bvlow(ngrid) ! Qu est-ce ?86 REAL :: dz 90 REAL :: bvlow(ngrid) ! initial N at each grid (not used) 91 REAL :: dz ! depth of the GW source if gwd_convective_source=.true. 87 92 88 93 INTEGER ii, jj, ll … … 98 103 INTEGER, parameter:: nw = nk * np *no ! Total numbers of gravity waves 99 104 INTEGER jk, jp, jo, jw 100 REAL kstar ! Control value to constrain the min horizontal wavenumber by the horizontal resolution 101 REAL, parameter:: kmax = 1.e-4 ! Max horizontal wavenumber=N/u, u(=30) zonal wind velocity at launch 102 !REAL, parameter:: kmax = 4.e-5 ! Max horizontal wavenumber=N/u, u(=70) zonal wind velocity at launch 103 REAL, parameter:: kmin = 2.e-6 ! Min horizontal wavenumber=1/sqrt(DxDy) Dx and Dy horizontal grid 104 REAL, parameter:: cmax = 30. ! Max horizontal absolute phase velocity=zonal wind at launch 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) 105 113 !REAL, parameter:: cmax = 100. ! Test for Saturn: Max horizontal absolute phase velocity 106 114 !REAL, parameter:: cmax = 70. ! Test for Saturn: Max horizontal absolute phase velocity 107 REAL, parameter:: cmin = 1. ! Min horizontal absolute phase velocity108 REAL cpha ! absolute phase velocity frequency109 REAL zk(nw, ngrid) ! horizontal wavenumber amplitude110 REAL zp(nw, ngrid) ! horizontal wavenumber angle111 REAL zo(nw, ngrid) ! absolute frequency115 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 112 120 113 121 REAL intr_freq_m(nw, ngrid) ! Waves Intr. freq. at the 1/2 lev below the full level (previous name: ZOM) … … 120 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) 121 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) 122 132 INTEGER launch ! Launching altitude 123 REAL, parameter:: xlaunch = 0.2 ! Control the launching altitude 133 REAL, parameter:: xlaunch = 0.2 ! Control the launching altitude by pressure 124 134 REAL, parameter:: zmaxth_top = 8000. ! Top of convective layer in m (approx.) 125 135 … … 128 138 129 139 ! 0.3.2 Parameters of waves dissipations 130 REAL, parameter:: sat = 1. ! saturation parameter 131 REAL, parameter:: rdiss = 1. ! coefficient of dissipation 132 REAL, parameter:: zoisec = 1.e-6 ! security for intrinsic freguency 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 133 145 134 146 ! 0.3.3 Background flow at 1/2 levels and vertical coordinate 135 147 REAL H0bis(ngrid, nlayer) ! characteristic height of the atmosphere 136 148 REAL, save:: H0 ! characteristic height of the atmosphere 149 !$OMP THREADPRIVATE(H0) 137 150 !REAL, parameter:: pr = 250 ! Reference pressure [Pa] 138 151 !REAL, parameter:: tr = 220. ! Reference temperature [K] … … 145 158 REAL vh(ngrid, nlayer + 1) ! Meridional wind at 1/2 levels 146 159 REAL ph(ngrid, nlayer + 1) ! Pressure at 1/2 levels 147 REAL, parameter:: psec = 1.e- 6 ! Security to avoid division by 0 pressure160 REAL, parameter:: psec = 1.e-12 ! Security to avoid division by 0 pressure(!!IMPORTANT: should be lower than the topmost layer's pressure) 148 161 REAL bv(ngrid, nlayer + 1) ! Brunt Vaisala freq. at 1/2 levels 149 REAL, parameter:: bvsec = 1.e- 5 ! Security to avoid negative BV162 REAL, parameter:: bvsec = 1.e-6 ! Security to avoid negative BV (!!IMPORTANT: tunable, depends on which Planet) 150 163 REAL href(nlayer + 1) ! Reference altitude for launching altitude 151 164 … … 156 169 157 170 LOGICAL, save :: firstcall = .true. 171 !$OMP THREADPRIVATE(firstcall) 172 LOGICAL, save :: gwd_convective_source = .false. 173 !$OMP THREADPRIVATE(gwd_convective_source) 158 174 159 175 !-------------------------------------------------------------------------------------------------- 160 176 ! 1. INITIALISATION 161 !------------------ 177 !-------------------------------------------------------------------------------------------------- 162 178 IF (firstcall) THEN 163 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 164 196 ! Characteristic vertical scale height 165 197 H0 = r * tr / g … … 175 207 gwd_convective_source = .false. 176 208 177 209 ! Compute subroutine's current values of temperature and winds 178 210 tt(:,:) = pt(:,:) + dtime * pdt(:,:) 179 211 uu(:,:) = pu(:,:) + dtime * pdu(:,:) 180 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 181 219 ! print*,'epflux_max just after firstcall:', epflux_max 182 220 183 221 !-------------------------------------------------------------------------------------------------- 184 222 ! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS 185 !---------------------------------------------------- 223 !-------------------------------------------------------------------------------------------------- 186 224 ! Pressure and Inv of pressure, temperature at 1/2 level 187 225 DO ll = 2, nlayer … … 205 243 206 244 ! Log-pressure vertical coordinate 207 DO ll = 1, nlayer + 1 208 zh(:,ll) = H0 * log(pr / (ph(:,ll) + psec)) 209 end DO 210 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 211 257 ! Winds and Brunt Vaisala frequency 212 258 DO ll = 2, nlayer … … 214 260 vh(:, ll) = 0.5 * (vv(:, ll) + vv(:, ll - 1)) ! Meridional wind 215 261 216 bv(:, ll) = 0.5 * (tt(:, ll) + tt(:, ll - 1))&217 * r**2 / cpp / H0**2&218 + (tt(:, ll) - tt(:, ll - 1)) &219 / (zh(:, ll) - zh(:, ll - 1)) & 220 * r / H0221 bv(:, ll) = sqrt(max(bvsec, bv(:,ll)))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 222 268 end DO 223 269 … … 229 275 vh(:, nlayer + 1) = vv(:, nlayer) 230 276 277 !----------------------------------------------------------------------------------------------------------------------- 231 278 ! 3. WAVES CHARACTERISTICS CHOSEN RANDOMLY 232 !----------------------------------------- 233 ! The mod function of here a weird arguments 234 ! are used to produce the waves characteristics 235 ! in a stochastic way 279 !----------------------------------------------------------------------------------------------------------------------- 280 ! The mod function of here a weird arguments are used to produce the waves characteristics in a stochastic way 236 281 DO jw = 1, nw 237 DO ii = 1, ngrid 282 !Angle 283 DO ii = 1, ngrid 238 284 239 ran_num_1 = mod(tt(ii, jw) * 10., 1.) 240 ran_num_2 = mod(tt(ii, jw) * 100., 1.) 241 242 ! angle (0 or pi so far) 243 zp(jw, ii) = (sign(1., 0.5 - ran_num_1) & 244 + 1.) * pi / 2. 245 ! horizontal wavenumber amplitude 246 ! TN+GG April/June 2020 247 ! "Individual waves are not supposed 248 ! to occupy the entire domain, but 249 ! only a faction of it" Lott+2012 250 251 !zk(jw, ii) = kmin + (kmax - kmin) * ran_num_2 252 kstar = pi / sqrt(cell_area(ii)) 253 zk(jw, ii) = max(kmin,kstar) + (kmax - max(kmin,kstar)) * ran_num_2 254 ! horizontal phase speed 255 cpha = 0. 256 DO jj = 1, na 257 ran_num_3 = mod(tt(ii, jw + 3 * jj)**2, 1.) 258 cpha = cpha + 2. * cmax * & 259 (ran_num_3 - 0.5) * & 260 sqrt(3.) / sqrt(na * 1.) 261 end DO 262 IF (cpha < 0.) THEN 263 cpha = - 1. * cpha 264 zp (jw, ii) = zp(jw, ii) + pi 265 ENDIF 266 ! Intrinsic frequency 267 zo(jw, ii) = cpha * zk(jw, ii) 268 ! Intrinsic frequency is imposed 269 zo(jw, ii) = zo(jw, ii) & 270 + zk(jw, ii) * cos(zp(jw, ii)) * uh(ii, launch) & 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) & 271 320 + zk(jw, ii) * sin(zp(jw, ii)) * vh(ii, launch) 272 ! Momentum flux at launch level273 epflux_0(jw, ii) = epflux_max &321 ! Momentum flux at launch level 322 epflux_0(jw, ii) = epflux_max & 274 323 * mod(100. * (uu(ii, jw)**2 + vv(ii, jw)**2), 1.) 275 276 end DO 277 end DO 324 ENDDO 325 end DO !DO jw = 1, nw 278 326 279 327 ! print*,'epflux_max just after waves charac. ramdon:', epflux_max 280 328 ! print*,'epflux_0 just after waves charac. ramdon:', epflux_0 329 330 !----------------------------------------------------------------------------------------------------------------------- 281 331 ! 4. COMPUTE THE FLUXES 282 !---------------------- 283 ! 4.1 Vertical velocity at launching altitude to ensure 284 ! the correct value to the imposed fluxes. 332 !----------------------------------------------------------------------------------------------------------------------- 333 ! 4.1 Vertical velocity at launching altitude to ensure the correct value to the imposed fluxes. 285 334 !------------------------------------------------------ 286 335 DO jw = 1, nw … … 362 411 ! No breaking (Eq. 6): 363 412 wwm(jw, :) & 364 ! Dissipation (Eq. 8): 365 * exp(-rdiss * pr / (ph(:, ll + 1) + ph (:, ll)) & 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) & 366 417 * ((bv(:, ll + 1) + bv (:, ll)) / 2.)**3 & 367 418 / max(abs(intr_freq_p(jw, :) + intr_freq_m(jw, :)) / 2., zoisec)**4 & … … 370 421 ! changes sign) 371 422 max(0., sign(1., intr_freq_p(jw, :) * intr_freq_m(jw, :))) & 372 ! Saturation (Eq. 12) 373 * abs(intr_freq_p(jw, :))**3 / bv(:, ll + 1) & 374 * exp(-zh(:, ll + 1) / H0) & 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) & 375 427 * sat**2 * kmin**2 / zk(jw, :)**4) 376 428 end DO 429 ! END OF W(KB)ARNING 377 430 378 431 ! Evaluate EP-flux from Eq. 7 and give the right orientation to the stress … … 388 441 v_epflux_tot(:, ll + 1) = v_epflux_tot(:, ll + 1) + v_epflux_p(jw, :) 389 442 east_gwstress(:, ll) = east_gwstress(:, ll) & 390 + max(0., u_epflux_p(jw, :)) / float(nw)443 + max(0., u_epflux_p(jw, ll)) / float(nw) 391 444 west_gwstress(:, ll) = west_gwstress(:, ll) & 392 445 + min(0., u_epflux_p(jw, ll)) / float(nw) 393 446 end DO 394 447 end DO ! DO LL = LAUNCH, nlayer - 1 395 396 ! print*,'u_epflux_tot just after launching:', u_epflux_tot397 ! print*,'v_epflux_tot just after launching:', v_epflux_tot398 ! print*,'u_epflux_p just after launching:', maxval(u_epflux_p), minval(u_epflux_p)399 ! print*,'v_epflux_p just after launching:', maxval(v_epflux_p), minval(v_epflux_p)400 448 401 449 ! 5. TENDENCY CALCULATIONS … … 436 484 !--------------------------------------------- 437 485 DO LL = 1, nlayer 438 d_u(:, LL) = (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL))! & 439 !d_u(:, LL) = G * (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL))! & 440 ! / (PH(:, LL + 1) - PH(:, LL)) * DTIME 441 d_v(:, LL) = (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL))! & 442 !d_v(:, LL) = G * (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL))! & 443 ! / (PH(:, LL + 1) - PH(:, LL)) * DTIME 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))) 444 491 ENDDO 445 492 d_t(:,:) = 0. … … 454 501 dv_nonoro_gwd(:,:) = d_v(:,:) 455 502 456 ! print*,'u_epflux_tot just after tendency:', u_epflux_tot457 ! print*,'v_epflux_tot just after tendency:', v_epflux_tot458 ! print*,'d_u just after tendency:', maxval(d_u(:,:)), minval(d_u)459 ! print*,'d_v just after tendency:', maxval(d_v(:,:)), minval(d_v)460 503 ! Cosmetic: evaluation of the cumulated stress 461 462 504 ZUSTR(:) = 0. 463 505 ZVSTR(:) = 0. 464 506 DO LL = 1, nlayer 465 ZUSTR(:) = ZUSTR(:) + D_U(:, LL) !/ g * (PH(:, LL + 1) - PH(:, LL))466 ZVSTR(:) = ZVSTR(:) + D_V(:, LL) !/ g * (PH(:, LL + 1) - PH(:, LL))507 ZUSTR(:) = ZUSTR(:) + D_U(:, LL) / g * (PH(:, LL + 1) - PH(:, LL)) 508 ZVSTR(:) = ZVSTR(:) + D_V(:, LL) / g * (PH(:, LL + 1) - PH(:, LL)) 467 509 ENDDO 468 510 -
trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90
r2542 r2595 60 60 waterrain, global1d, szangle 61 61 use nonoro_gwd_ran_mod, only: nonoro_gwd_ran 62 use conc_mod 62 use conc_mod, only: rnew, cpnew, ini_conc_mod 63 63 use phys_state_var_mod 64 64 use callcorrk_mod, only: callcorrk … … 1269 1269 IF (calllott_nonoro) THEN 1270 1270 1271 CALL nonoro_gwd_ran(ngrid,nlayer,ptimestep,pplay, & 1271 CALL nonoro_gwd_ran(ngrid,nlayer,ptimestep, & 1272 cpnew, rnew, & 1273 pplay, & 1272 1274 zmax_th, &! max altitude reached by thermals (m) 1273 1275 pt, pu, pv, & … … 1279 1281 pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) & 1280 1282 + d_t_hin(1:ngrid,1:nlayer) 1281 ! d_t_hin(:,:)= d_t_hin(:,:)/ptimestep ! K/s (for outputs?)1282 1283 pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) & 1283 1284 + d_u_hin(1:ngrid,1:nlayer) 1284 ! d_u_hin(:,:)= d_u_hin(:,:)/ptimestep ! (m/s)/s (for outputs?)1285 1285 pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer) & 1286 1286 + d_v_hin(1:ngrid,1:nlayer) 1287 ! d_v_hin(:,:)= d_v_hin(:,:)/ptimestep ! (m/s)/s (for outputs?)1288 1287 print*,'du_nonoro: max ', maxval(d_u_hin), 'min ', minval(d_u_hin) 1289 1288 print*,'dv_nonoro: max ', maxval(d_v_hin), 'min ', minval(d_v_hin) … … 2176 2175 ! GW non-oro outputs 2177 2176 if (calllott_nonoro) then 2178 call WRITEDIAGFI(ngrid,"dugwno","GW non-oro dU","m/s2", 3,d_u_hin /ptimestep)2179 call WRITEDIAGFI(ngrid,"dvgwno","GW non-oro dV","m/s2", 3,d_v_hin /ptimestep)2177 call WRITEDIAGFI(ngrid,"dugwno","GW non-oro dU","m/s2", 3,d_u_hin) 2178 call WRITEDIAGFI(ngrid,"dvgwno","GW non-oro dV","m/s2", 3,d_v_hin) 2180 2179 endif 2181 2180
Note: See TracChangeset
for help on using the changeset viewer.
