- Timestamp:
- Dec 15, 2021, 3:08:53 PM (4 years ago)
- Location:
- trunk/LMDZ.MARS
- Files:
-
- 3 edited
-
README (modified) (1 diff)
-
libf/phymars/nonoro_gwd_ran_mod.F90 (modified) (16 diffs)
-
libf/phymars/physiq_mod.F (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/README
r2593 r2594 3542 3542 == 10/12/2021 == MV 3543 3543 hdo_surfex_mod.F, vdifc_mod.F: addings to account for the effect of kinetics in the fractionation by condensation of HDO at the surface 3544 3545 == 15/12/2021 == JL 3546 Fixes and improvements in the Non-orographic GW scheme, namely: 3547 - increments are not tendencies 3548 - missing rho factor in EP flux computation 3549 - missing rho at launch altitude 3550 - changed inputs, because R and Cp are needed to compute rho and BV -
trunk/LMDZ.MARS/libf/phymars/nonoro_gwd_ran_mod.F90
r2578 r2594 12 12 CONTAINS 13 13 14 SUBROUTINE NONORO_GWD_RAN(ngrid,nlayer,DTIME, pp, &14 SUBROUTINE NONORO_GWD_RAN(ngrid,nlayer,DTIME, cpnew, rnew, pp, & 15 15 zmax_therm, pt, pu, pv, pdt, pdu, pdv, & 16 16 zustr,zvstr,d_t, d_u, d_v) … … 34 34 ! calculation using MOD 35 35 ! - adding east_gwstress and 36 ! west_gwstress variables 36 ! west_gwstress variables 37 ! UPDATED J.LIU 12/2021 The rho (density) at the specific 38 ! locations is introduced. The equation 39 ! of EP-flux is corrected by adding the 40 ! term of density at launch (source) 41 ! altitude. 42 ! 37 43 !--------------------------------------------------------------------------------- 38 44 39 use comcstfi_h, only: g, pi, cpp,r45 use comcstfi_h, only: g, pi, r 40 46 USE ioipsl_getin_p_mod, ONLY : getin_p 41 47 use assert_m, only : assert 42 48 use vertical_layers_mod, only : presnivs 43 49 use geometry_mod, only: cell_area 50 #ifdef CPP_XIOS 51 use xios_output_mod, only: send_xios_field 52 #endif 44 53 implicit none 45 46 54 include "yoegwd.h" 47 55 include "callkeys.h" 48 56 49 CHARACTER (LEN=20) :: modname=' flott_gwd_rando'57 CHARACTER (LEN=20) :: modname='nonoro_gwd_ran' 50 58 CHARACTER (LEN=80) :: abort_message 51 59 … … 54 62 55 63 ! 0.1 INPUTS 56 INTEGER, intent(in):: ngrid, nlayer 57 REAL, intent(in):: DTIME ! Time step of the Physics 58 REAL, intent(in):: zmax_therm(ngrid) ! altitude of max velocity thermals (m) 59 REAL, intent(in):: pp(ngrid,nlayer) ! Pressure at full levels 60 REAL, intent(in):: pt(ngrid,nlayer) ! Temp at full levels 61 REAL, intent(in):: pu(ngrid,nlayer),pv(ngrid,nlayer) ! Hor winds at full levels 62 REAL,INTENT(in) :: pdt(ngrid,nlayer) ! tendency on temperature 63 REAL,INTENT(in) :: pdu(ngrid,nlayer) ! tendency on zonal wind 64 REAL,INTENT(in) :: pdv(ngrid,nlayer) ! tendency on meridional wind 64 INTEGER, intent(in):: ngrid ! number of atmospheric columns 65 INTEGER, intent(in):: nlayer ! number of atmospheric layers 66 REAL, intent(in):: DTIME ! Time step of the Physics(s) 67 REAL, intent(in):: zmax_therm(ngrid) ! Altitude of max velocity thermals (m) 68 REAL,INTENT(IN) :: cpnew(ngrid,nlayer)! Cp of the atmosphere 69 REAL,INTENT(IN) :: rnew(ngrid,nlayer) ! R of the atmosphere 70 REAL, intent(in):: pp(ngrid,nlayer) ! Pressure at full levels(Pa) 71 REAL, intent(in):: pt(ngrid,nlayer) ! Temp at full levels(K) 72 REAL, intent(in):: pu(ngrid,nlayer) ! Zonal wind at full levels(m/s) 73 REAL, intent(in):: pv(ngrid,nlayer) ! Meridional winds at full levels(m/s) 74 REAL,INTENT(in) :: pdt(ngrid,nlayer) ! Tendency on temperature (K/s) 75 REAL,INTENT(in) :: pdu(ngrid,nlayer) ! Tendency on zonal wind (m/s/s) 76 REAL,INTENT(in) :: pdv(ngrid,nlayer) ! Tendency on meridional wind (m/s/s) 65 77 66 78 ! 0.2 OUTPUTS 67 REAL, intent(out):: zustr(ngrid), zvstr(ngrid) ! Surface Stresses 68 REAL, intent(out):: d_t(ngrid, nlayer) ! Tendency on Temp. 69 REAL, intent(out):: d_u(ngrid, nlayer), d_v(ngrid, nlayer) ! tendency on winds 70 71 ! O.3 INTERNAL ARRAYS 72 REAL :: TT(ngrid, nlayer) ! Temp at full levels 73 REAL :: UU(ngrid, nlayer) , VV(ngrid, nlayer) ! Hor winds at full levels 74 REAL :: BVLOW(ngrid) 75 REAL :: DZ 79 REAL, intent(out):: zustr(ngrid) ! Zonal surface stress 80 REAL, intent(out):: zvstr(ngrid) ! Meridional surface stress 81 REAL, intent(out):: d_t(ngrid, nlayer) ! Tendency on temperature (K/s) due to gravity waves (not used set to zero) 82 REAL, intent(out):: d_u(ngrid, nlayer) ! Tendency on zonal wind (m/s/s) due to gravity waves 83 REAL, intent(out):: d_v(ngrid, nlayer) ! Tendency on meridional wind (m/s/s) due to gravity waves 84 85 ! 0.3 INTERNAL ARRAYS 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 winds at full levels 90 REAL :: BVLOW(ngrid) ! initial N at each grid (not used) 76 91 77 92 INTEGER II, JJ, LL 78 93 79 94 ! 0.3.0 TIME SCALE OF THE LIFE CYCLE OF THE WAVES PARAMETERIZED 80 81 95 REAL, parameter:: DELTAT = 24. * 3600. 82 96 83 97 ! 0.3.1 GRAVITY-WAVES SPECIFICATIONS 84 85 98 INTEGER, PARAMETER:: NK = 2 ! number of horizontal wavenumbers 86 99 INTEGER, PARAMETER:: NP = 2 ! directions (eastward and westward) phase speed 87 100 INTEGER, PARAMETER:: NO = 2 ! absolute values of phase speed 101 INTEGER, PARAMETER:: NA = 5 ! number of realizations to get the phase speed 88 102 INTEGER, PARAMETER:: NW = NK * NP * NO ! Total numbers of gravity waves 89 INTEGER JK, JP, JO, JW 90 INTEGER, PARAMETER:: NA = 5 ! number of realizations to get the phase speed 91 REAL, parameter:: kmax = 7.e-4 ! Max horizontal wavenumber 92 REAL, parameter:: kmin = 2.e-5 ! Min horizontal wavenumber 103 INTEGER JK, JP, JO, JW ! Loop index for NK,NP,NO, and NW 104 105 REAL, save :: kmax ! Max horizontal wavenumber (lambda min,lambda=2*pi/kmax),kmax=N/u, u(=30~50) zonal wind velocity at launch altitude 106 !$OMP THREADPRIVATE(kmax) 107 REAL, parameter:: kmin = 2.e-5 ! Min horizontal wavenumber (lambda max = 314 km,lambda=2*pi/kmin) 93 108 REAL kstar ! Min horizontal wavenumber constrained by horizontal resolution 94 REAL, parameter:: cmax = 30. ! Max horizontal absolute phase velocity 95 REAL, parameter:: cmin = 1. ! Min horizontal absolute phase velocity 96 REAL CPHA ! absolute PHASE VELOCITY frequency 97 REAL ZK(NW, ngrid) ! Horizontal wavenumber amplitude 98 REAL ZP(NW, ngrid) ! Horizontal wavenumber angle 99 REAL ZO(NW, ngrid) ! Absolute frequency 109 REAL, save :: cmax ! Max horizontal absolute phase velocity (at the maginitide of zonal wind u at the launch altitude) 110 !$OMP THREADPRIVATE(cmax) 111 REAL, parameter:: cmin = 1. ! Min horizontal absolute phase velocity (not used) 112 REAL CPHA ! absolute PHASE VELOCITY frequency 113 REAL ZK(NW, ngrid) ! Horizontal wavenumber amplitude 114 REAL ZP(NW, ngrid) ! Horizontal wavenumber angle 115 REAL ZO(NW, ngrid) ! Absolute frequency 100 116 101 117 REAL intr_freq_m(nw, ngrid) ! Waves Intr. freq. at the 1/2 lev below the full level (previous name: ZOM) … … 108 124 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) 109 125 REAL epflux_0(nw, ngrid) ! Fluxes at launching level (previous name: RUW0) 110 REAL, save :: epflux_max ! Max EP flux value at launching altitude (previous name: RUWMAX) 111 INTEGER LAUNCH ! Launching altitude 112 REAL, parameter:: xlaunch = 0.4 ! Control the launching altitude 113 REAL, parameter:: zmaxth_top = 8000. ! Top of convective layer (approx.) 114 115 116 REAL PREC(ngrid) 117 REAL PRMAX ! Maximum value of PREC, and for which our linear formula 126 REAL, save :: epflux_max ! Max EP flux value at launching altitude (previous name: RUWMAX, tunable) 127 !$OMP THREADPRIVATE(epflux_max) 128 INTEGER LAUNCH ! Launching altitude 129 REAL, parameter:: xlaunch = 0.4 ! Control the launching altitude by pressure 130 REAL, parameter:: zmaxth_top = 8000. ! Top of convective layer (approx. not used) 118 131 119 132 120 133 ! 0.3.2 PARAMETERS OF WAVES DISSIPATIONS 121 REAL, parameter:: sat = 1. ! saturation parameter 122 REAL, parameter:: rdiss = 1. ! coefficient of dissipation 123 REAL, parameter:: zoisec = 1.e-6 ! security for intrinsic freguency 134 REAL, save :: sat ! saturation parameter(tunable) 135 !$OMP THREADPRIVATE(sat) 136 REAL, save :: rdiss ! dissipation coefficient (tunable) 137 !$OMP THREADPRIVATE(rdiss) 138 REAL, parameter:: zoisec = 1.e-10 ! security for intrinsic frequency 124 139 125 140 ! 0.3.3 Background flow at 1/2 levels and vertical coordinate 126 REAL H0bis(ngrid, nlayer) ! Characteristic Height of the atmosphere 127 REAL, save:: H0 ! Characteristic Height of the atmosphere 128 REAL, parameter:: pr = 250 ! Reference pressure [Pa] 129 REAL, parameter:: tr = 220. ! Reference temperature [K] 141 REAL H0bis(ngrid, nlayer) ! Characteristic Height of the atmosphere (specific locations) 142 REAL, save:: H0 ! Characteristic Height of the atmosphere (constant) 143 !$OMP THREADPRIVATE(H0) 144 REAL, parameter:: pr = 250 ! Reference pressure [Pa] 145 REAL, parameter:: tr = 220. ! Reference temperature [K] 130 146 REAL ZH(ngrid, nlayer + 1) ! Log-pressure altitude (constant H0) 131 REAL ZHbis(ngrid, nlayer + 1) ! Log-pressure altitude (varying H )147 REAL ZHbis(ngrid, nlayer + 1) ! Log-pressure altitude (varying H0bis) 132 148 REAL UH(ngrid, nlayer + 1) ! zonal wind at 1/2 levels 133 149 REAL VH(ngrid, nlayer + 1) ! meridional wind at 1/2 levels 134 150 REAL PH(ngrid, nlayer + 1) ! Pressure at 1/2 levels 135 REAL, parameter:: psec = 1.e-6 ! Security to avoid division by 0 pressure 136 REAL BV(ngrid, nlayer + 1) ! Brunt Vaisala freq. (BVF) at 1/2 levels 137 REAL, parameter:: bvsec = 1.e-5 ! Security to avoid negative BV 138 REAL HREF(nlayer + 1) ! Reference altitude for launching alt. 139 140 141 ! COSMETICS TO DIAGNOSE EACH WAVES CONTRIBUTION. 142 logical,save :: output=.false. 143 ! CAUTION ! IF output is .true. THEN change NO to 10 at least ! 144 character*14 outform 145 character*2 str2 146 integer ieq 151 REAL, parameter:: psec = 1.e-20 ! Security to avoid division by 0 pressure(!!IMPORTANT: should be lower than the topmost layer's pressure) 152 REAL BV(ngrid, nlayer + 1) ! Brunt Vaisala freq. (N^2) at 1/2 levels 153 REAL, parameter:: bvsec = 1.e-5 ! Security to avoid negative BV 154 REAL HREF(nlayer + 1) ! Reference pressure for launching altitude 155 147 156 148 157 REAL RAN_NUM_1,RAN_NUM_2,RAN_NUM_3 … … 152 161 153 162 154 !-----------------------------------------------------------------155 156 163 !----------------------------------------------------------------------------------------------------------------------- 164 ! 1. INITIALISATIONS 165 !----------------------------------------------------------------------------------------------------------------------- 157 166 IF (firstcall) THEN 158 167 write(*,*) "nonoro_gwd_ran: FLott non-oro GW scheme is active!" … … 160 169 call getin_p("nonoro_gwd_epflux_max", epflux_max) 161 170 write(*,*) "nonoro_gwd_ran: epflux_max=", epflux_max 171 sat = 1. ! default gravity waves saturation value !! 172 call getin_p("nonoro_gwd_sat", sat) 173 write(*,*) "nonoro_gwd_ran: sat=", sat 174 cmax = 50. ! default gravity waves phase velocity value !! 175 call getin_p("nonoro_gwd_cmax", cmax) 176 write(*,*) "nonoro_gwd_ran: cmax=", cmax 177 rdiss=0.01 ! default coefficient of dissipation !! 178 call getin_p("nonoro_gwd_rdiss", rdiss) 179 write(*,*) "nonoro_gwd_ran: rdiss=", rdiss 180 kmax=7.E-4 ! default Max horizontal wavenumber !! 181 call getin_p("nonoro_gwd_kmax", kmax) 182 write(*,*) "nonoro_gwd_ran: kmax=", kmax 183 162 184 ! Characteristic vertical scale height 163 185 H0 = r * tr / g … … 172 194 ENDIF 173 195 174 gwd_convective_source=.false. 175 176 ! Compute current values of temperature and winds 196 ! Compute current values of temperature and winds 177 197 tt(:,:)=pt(:,:)+dtime*pdt(:,:) 178 198 uu(:,:)=pu(:,:)+dtime*pdu(:,:) 179 199 vv(:,:)=pv(:,:)+dtime*pdv(:,:) 180 181 182 183 ! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS 184 !------------------------------------------------------------- 185 186 !Online output 187 if (output) OPEN(11,file="impact-gwno.dat") 188 189 ! Pressure and Inv of pressure, Temperature / at 1/2 level 200 ! Compute the real mass density by rho=p/R(T)T 201 DO ll=1,nlayer 202 DO ii=1,ngrid 203 rho(ii,ll) = pp(ii,ll)/(rnew(ii,ll)*tt(ii,ll)) 204 ENDDO 205 ENDDO 206 ! print*,'epflux_max just after firstcall:', epflux_max 207 208 !----------------------------------------------------------------------------------------------------------------------- 209 ! 2. EVALUATION OF THE BACKGROUND FLOW AT SEMI-LEVELS 210 !----------------------------------------------------------------------------------------------------------------------- 211 ! Pressure and inverse of pressure at 1/2 level 190 212 DO LL = 2, nlayer 191 213 PH(:, LL) = EXP((LOG(PP(:, LL)) + LOG(PP(:, LL - 1))) / 2.) 192 214 end DO 193 194 215 PH(:, nlayer + 1) = 0. 195 216 PH(:, 1) = 2. * PP(:, 1) - PH(:, 2) 196 217 197 ! Launching altitude 198 218 call writediagfi(ngrid,'nonoro_pp','nonoro_pp', 'm',3,PP(:,1:nlayer)) 219 call writediagfi(ngrid,'nonoro_ph','nonoro_ph', 'm',3,PH(:,1:nlayer)) 220 221 ! Launching level for reproductible case 199 222 !Pour revenir a la version non reproductible en changeant le nombre de 200 223 !process … … 209 232 DO LL = 1, nlayer 210 233 IF (HREF(LL) / HREF(1) > XLAUNCH) LAUNCH = LL 211 ENDDO 212 213 if (output) print*, " WE ARE IN FLOTT GW SCHEME " 214 234 ENDDO 235 215 236 ! Log pressure vert. coordinate 216 DO LL = 1, nlayer + 1 217 ZH(:, LL) = H0 * LOG(PR / (PH(:, LL) + PSEC)) 218 end DO 219 220 if (output) then 221 ! altitude above surface 222 ZHbis(:,1) = 0. 223 DO LL = 2, nlayer + 1 224 H0bis(:, LL-1) = r * TT(:, LL-1) / g 225 ZHbis(:, LL) = ZHbis(:, LL-1) & 237 ZH(:,1) = 0. 238 DO LL = 2, nlayer + 1 239 !ZH(:, LL) = H0 * LOG(PR / (PH(:, LL) + PSEC)) 240 H0bis(:, LL-1) = r * TT(:, LL-1) / g 241 ZH(:, LL) = ZH(:, LL-1) & 226 242 + H0bis(:, LL-1)*(PH(:, LL-1)-PH(:,LL))/PP(:, LL-1) 227 end DO 228 endif 229 230 ! Winds and BV frequency 243 end DO 244 ZH(:, 1) = H0 * LOG(PR / (PH(:, 1) + PSEC)) 245 246 call writediagfi(ngrid,'nonoro_zh','nonoro_zh', 'm',3,ZH(:,2:nlayer+1)) 247 248 ! Winds and Brunt Vaisala frequency 231 249 DO LL = 2, nlayer 232 UH(:, LL) = 0.5 * (UU(:, LL) + UU(:, LL - 1)) ! Zonal wind 233 VH(:, LL) = 0.5 * (VV(:, LL) + VV(:, LL - 1)) ! Meridional wind 234 ! GG test 235 !print*, 'TT, UH, VH, ZH at launch', TT(ngrid/2,LAUNCH), UH(ngrid/2,LAUNCH),VH(ngrid/2, LAUNCH), ZH(ngrid/2,LAUNCH) 236 ! BVSEC: BV Frequency 237 BV(:, LL) = 0.5 * (TT(:, LL) + TT(:, LL - 1)) & 238 * r**2 / cpp / H0**2 + (TT(:, LL) & 239 - TT(:, LL - 1)) / (ZH(:, LL) - ZH(:, LL - 1)) * r / H0 240 BV(:,LL) =SQRT(MAX(BVSEC,BV(:,LL))) 241 end DO 242 !GG test 243 !print*, 'BV freq in flott_gwnoro:',LAUNCH, BV(ngrid/2, LAUNCH) 244 250 UH(:, LL) = 0.5 * (UU(:, LL) + UU(:, LL - 1)) ! Zonal wind 251 VH(:, LL) = 0.5 * (VV(:, LL) + VV(:, LL - 1)) ! Meridional wind 252 ! Brunt Vaisala frequency (=g/T*[dT/dz + g/cp] ) 253 BV(:, LL)= G/(0.5 * (TT(:, LL) + TT(:, LL - 1))) & 254 *((TT(:, LL) - TT(:, LL - 1)) / (ZH(:, LL) - ZH(:, LL - 1))+ & 255 G / cpnew(:,LL)) 256 257 BV(:,LL) =MAX(1.E-12,BV(:,LL)) ! to ensure that BV is positive 258 BV(:,LL) = MAX(BVSEC,SQRT(BV(:,LL))) ! make sure it is not too small 259 end DO 245 260 BV(:, 1) = BV(:, 2) 246 261 UH(:, 1) = 0. … … 250 265 VH(:, nlayer + 1) = VV(:, nlayer) 251 266 252 253 ! 3. WAVES CHARACTERISTICS CHOSEN RANDOMLY 254 !------------------------------------------- 255 256 ! The mod function of here a weird arguments 257 ! are used to produce the waves characteristics 258 ! in a stochastic way 259 267 call writediagfi(ngrid,'nonoro_bv','nonoro_bv', 'm',3,BV(:,2:nlayer+1)) 268 269 !----------------------------------------------------------------------------------------------------------------------- 270 ! 3. WAVES CHARACTERISTICS CHOSEN RANDOMLY 271 !----------------------------------------------------------------------------------------------------------------------- 272 ! The mod function of here a weird arguments are used to produce the waves characteristics in a stochastic way 260 273 DO JW = 1, NW 261 274 ! Angle … … 264 277 RAN_NUM_1=MOD(TT(II, JW) * 10., 1.) 265 278 RAN_NUM_2= MOD(TT(II, JW) * 100., 1.) 266 ZP(JW, II) = (SIGN(1., 0.5 - RAN_NUM_1) + 1.) &267 * PI / 2.279 ZP(JW, II) = (SIGN(1., 0.5 - RAN_NUM_1) + 1.)* PI / 2. 280 268 281 ! Horizontal wavenumber amplitude 269 282 ! From Venus model: TN+GG April/June 2020 (rev2381) 270 ! "Individual waves are not supposed to occupy the 271 ! entire domain, but only a fraction of it" (Lott+2012) 272 273 !ZK(JW, II) = KMIN + (KMAX - KMIN) *RAN_NUM_2 283 ! "Individual waves are not supposed to occupy the entire domain, but only a fraction of it" (Lott+2012) 284 ! ZK(JW, II) = KMIN + (KMAX - KMIN) *RAN_NUM_2 274 285 KSTAR = PI/SQRT(cell_area(II)) 275 286 ZK(JW, II) = MAX(KMIN,KSTAR) + (KMAX - MAX(KMIN,KSTAR)) *RAN_NUM_2 276 ! Horizontal phase speed 287 288 ! Horizontal phase speed 289 ! this computation allows to get a gaussian distribution centered on 0 (right ?) 290 ! then cmin is not useful, and it favors small values. 277 291 CPHA = 0. 278 292 DO JJ = 1, NA 279 293 RAN_NUM_3=MOD(TT(II, JW+3*JJ)**2, 1.) 280 CPHA = CPHA + & 281 CMAX*2.*(RAN_NUM_3 -0.5)*SQRT(3.)/SQRT(NA*1.) 294 CPHA = CPHA + 2.*CMAX* & 295 (RAN_NUM_3 -0.5)* & 296 SQRT(3.)/SQRT(NA*1.) 282 297 END DO 283 298 IF (CPHA.LT.0.) THEN … … 285 300 ZP(JW,II) = ZP(JW,II) + PI 286 301 ENDIF 287 !Online output: linear 288 if (output) CPHA = CMIN + (CMAX - CMIN) * (JO-1)/(NO-1) 302 ! otherwise, with the computation below, we get a uniform distribution between cmin and cmax. 303 ! ran_num_3 = mod(tt(ii, jw)**2, 1.) 304 ! cpha = cmin + (cmax - cmin) * ran_num_3 305 289 306 ! Intrinsic frequency 290 307 ZO(JW, II) = CPHA * ZK(JW, II) 291 308 ! Intrinsic frequency is imposed 292 ZO(JW, II) = ZO(JW, II) & 293 + ZK(JW, II) * COS(ZP(JW, II)) * UH(II, LAUNCH) & 294 + ZK(JW, II) * SIN(ZP(JW, II)) * VH(II, LAUNCH) 295 ! Momentum flux at launch lev 296 ! epflux_0(JW, II) = epflux_max / REAL(NW) & 297 epflux_0(JW, II) = epflux_max & 298 * MOD(100. * (UU(II, JW)**2 + VV(II, JW)**2), 1.) 299 !Online output: fixed to max 300 if (output) epflux_0(JW, II) = epflux_max 301 ENDDO 309 ZO(JW, II) = ZO(JW, II) & 310 + ZK(JW, II) * COS(ZP(JW, II)) * UH(II, LAUNCH) & 311 + ZK(JW, II) * SIN(ZP(JW, II)) * VH(II, LAUNCH) 312 313 ! Momentum flux at launch level 314 epflux_0(JW, II) = epflux_max & 315 * MOD(100. * (UU(II, JW)**2 + VV(II, JW)**2), 1.) 316 ENDDO 302 317 end DO 303 304 ! 4. COMPUTE THE FLUXES 305 !--------------------------306 307 ! 4.1 Vertical velocity at launching altitude to ensure 308 ! 309 ! 318 ! print*,'epflux_0 just after waves charac. ramdon:', epflux_0 319 320 !----------------------------------------------------------------------------------------------------------------------- 321 ! 4. COMPUTE THE FLUXES 322 !----------------------------------------------------------------------------------------------------------------------- 323 ! 4.1 Vertical velocity at launching altitude to ensure the correct value to the imposed fluxes. 324 !------------------------------------------------------ 310 325 DO JW = 1, NW 311 326 ! Evaluate intrinsic frequency at launching altitude: 312 intr_freq_p(JW, :) = ZO(JW, :) & 313 - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LAUNCH) & 314 - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LAUNCH) 315 end DO 316 317 IF (gwd_convective_source) THEN 318 DO JW = 1, NW 319 ! VERSION WITH CONVECTIVE SOURCE (designed for Earth) 320 321 ! Vertical velocity at launch level, value to ensure the 322 ! imposed mmt flux factor related to the convective forcing: 323 ! precipitations. 324 325 ! tanh limitation to values above prmax: 326 ! WWP(JW, :) = epflux_0(JW, :) & 327 ! * (r / cpp / H0 * RLVTT * PRMAX * TANH(PREC(:) / PRMAX))**2 328 ! Here, we neglected the kinetic energy providing of the thermodynamic 329 ! phase change 330 331 ! 332 333 ! Factor related to the characteristics of the waves: 334 WWP(JW, :) = WWP(JW, :) * ZK(JW, :)**3 / KMIN / BVLOW(:) & 335 / MAX(ABS(intr_freq_p(JW, :)), ZOISEC)**3 336 337 ! Moderation by the depth of the source (dz here): 338 WWP(JW, :) = WWP(JW, :) & 339 * EXP(- BVLOW(:)**2 / MAX(ABS(intr_freq_p(JW, :)), ZOISEC)**2 & 340 * ZK(JW, :)**2 * DZ**2) 341 342 ! Put the stress in the right direction: 343 u_epflux_p(JW, :) = intr_freq_p(JW, :) / MAX(ABS(intr_freq_p(JW, :)), ZOISEC)**2 & 344 * BV(:, LAUNCH) * COS(ZP(JW, :)) * WWP(JW, :)**2 345 v_epflux_p(JW, :) = intr_freq_p(JW, :) / MAX(ABS(intr_freq_p(JW, :)), ZOISEC)**2 & 346 * BV(:, LAUNCH) * SIN(ZP(JW, :)) * WWP(JW, :)**2 347 end DO 348 ELSE ! VERSION WITHOUT CONVECTIVE SOURCE 327 intr_freq_p(JW, :) = ZO(JW, :) & 328 - ZK(JW, :) * COS(ZP(JW, :)) * UH(:, LAUNCH) & 329 - ZK(JW, :) * SIN(ZP(JW, :)) * VH(:, LAUNCH) 330 end DO 331 332 ! VERSION WITHOUT CONVECTIVE SOURCE 349 333 ! Vertical velocity at launch level, value to ensure the imposed 350 334 ! mom flux: … … 354 338 u_epflux_p(JW, :) = COS(ZP(JW, :)) * SIGN(1., intr_freq_p(JW, :)) * epflux_0(JW, :) 355 339 v_epflux_p(JW, :) = SIN(ZP(JW, :)) * SIGN(1., intr_freq_p(JW, :)) * epflux_0(JW, :) 356 357 340 end DO 358 ENDIF341 359 342 ! 4.2 Initial flux at launching altitude 360 343 !------------------------------------------------------ 361 344 u_epflux_tot(:, LAUNCH) = 0 362 345 v_epflux_tot(:, LAUNCH) = 0 … … 366 349 end DO 367 350 368 ! 4.3 Loop over altitudes, with passage from one level to the 369 ! next done by i) conserving the EP flux, ii) dissipating 370 ! a little, iii) testing critical levels, and vi) testing 371 ! the breaking. 372 373 !Online output 374 if (output) then 375 ieq=ngrid/2+1 376 write(str2,'(i2)') NW+2 377 outform="("//str2//"(E12.4,1X))" 378 WRITE(11,outform) ZH(IEQ, 1) / 1000., ZHbis(IEQ, 1) / 1000., & 379 (ZO(JW, IEQ)/ZK(JW, IEQ)*COS(ZP(JW, IEQ)), JW = 1, NW) 380 endif 381 351 ! 4.3 Loop over altitudes, with passage from one level to the next done by: 352 !---------------------------------------------------------------------------- 353 ! i) conserving the EP flux, 354 ! ii) dissipating a little, 355 ! iii) testing critical levels, 356 ! iv) testing the breaking. 357 !---------------------------------------------------------------------------- 382 358 DO LL = LAUNCH, nlayer - 1 383 384 385 ! W(KB)ARNING: ALL THE PHYSICS IS HERE (PASSAGE FROM ONE LEVEL 386 ! TO THE NEXT) 359 ! W(KB)ARNING: ALL THE PHYSICS IS HERE (PASSAGE FROM ONE LEVEL TO THE NEXT) 387 360 DO JW = 1, NW 388 361 intr_freq_m(JW, :) = intr_freq_p(JW, :) 389 362 WWM(JW, :) = WWP(JW, :) 390 363 ! Intrinsic Frequency 391 intr_freq_p(JW, :) = ZO(JW, :) - ZK(JW, :) * COS(ZP(JW,:)) * UH(:, LL + 1) & 392 - ZK(JW, :) * SIN(ZP(JW,:)) * VH(:, LL + 1) 393 394 WWP(JW, :) = MIN( & 395 ! No breaking (Eq.6) 396 WWM(JW, :) & 397 ! Dissipation (Eq. 8): 398 * EXP(- RDISS * PR / (PH(:, LL + 1) + PH(:, LL)) & 399 * ((BV(:, LL + 1) + BV(:, LL)) / 2.)**3 & 400 / MAX(ABS(intr_freq_p(JW, :) + intr_freq_m(JW, :)) / 2., ZOISEC)**4 & 401 * ZK(JW, :)**3 * (ZH(:, LL + 1) - ZH(:, LL))), & 402 ! Critical levels (forced to zero if intrinsic frequency changes sign) 403 MAX(0., SIGN(1., intr_freq_p(JW, :) * intr_freq_m(JW, :))) & 404 ! Saturation (Eq. 12) 405 * ABS(intr_freq_p(JW, :))**3 /BV(:, LL+1) & 406 * EXP(-ZH(:, LL + 1)/H0) * SAT**2*KMIN**2/ZK(JW, :)**4) 364 intr_freq_p(JW, :) = ZO(JW, :) - ZK(JW, :) * COS(ZP(JW,:)) * UH(:, LL + 1) & 365 - ZK(JW, :) * SIN(ZP(JW,:)) * VH(:, LL + 1) 366 ! Vertical velocity in flux formulation 367 WWP(JW, :) = MIN( & 368 ! No breaking (Eq.6) 369 WWM(JW, :) & 370 ! Dissipation (Eq. 8)(real rho used here rather than pressure 371 ! parametration because the original code has a bug if the density of 372 ! the planet at the launch altitude not approximate 1): ! 373 * EXP(- RDISS*2./rho(:, LL + 1) & 374 * ((BV(:, LL + 1) + BV(:, LL)) / 2.)**3 & 375 / MAX(ABS(intr_freq_p(JW, :) + intr_freq_m(JW, :)) / 2., ZOISEC)**4 & 376 * ZK(JW, :)**3 * (ZH(:, LL + 1) - ZH(:, LL))), & 377 ! Critical levels (forced to zero if intrinsic frequency 378 ! changes sign) 379 MAX(0., SIGN(1., intr_freq_p(JW, :) * intr_freq_m(JW, :))) & 380 ! Saturation (Eq. 12) (rho at launch altitude is imposed by J.Liu. 381 ! Same reason with Eq. 8) 382 * ABS(intr_freq_p(JW, :))**3 / (2.*BV(:, LL+1)) & 383 * rho(:,launch)*exp(-zh(:, ll + 1) / H0) & 384 * SAT**2 *KMIN**2 / ZK(JW, :)**4) 407 385 end DO 408 409 386 ! END OF W(KB)ARNING 410 ! Evaluate EP-flux from Eq. 7 and 411 ! Give the right orientation to the stress387 388 ! Evaluate EP-flux from Eq. 7 and give the right orientation to the stress 412 389 DO JW = 1, NW 413 390 u_epflux_p(JW, :) = SIGN(1.,intr_freq_p(JW, :)) * COS(ZP(JW, :)) * WWP(JW, :) … … 417 394 u_epflux_tot(:, LL + 1) = 0. 418 395 v_epflux_tot(:, LL + 1) = 0. 419 420 396 DO JW = 1, NW 421 397 u_epflux_tot(:, LL + 1) = u_epflux_tot(:, LL + 1) + u_epflux_p(JW, :) … … 423 399 EAST_GWSTRESS(:, LL)=EAST_GWSTRESS(:, LL)+MAX(0.,u_epflux_p(JW,:))/FLOAT(NW) 424 400 WEST_GWSTRESS(:, LL)=WEST_GWSTRESS(:, LL)+MIN(0.,u_epflux_p(JW,:))/FLOAT(NW) 425 end DO 426 !Online output 427 if (output) then 428 do JW=1,NW 429 if(u_epflux_p(JW, IEQ).gt.0.) then 430 u_epflux_p(JW, IEQ) = max(u_epflux_p(JW, IEQ), 1.e-99) 431 else 432 u_epflux_p(JW, IEQ) = min(u_epflux_p(JW, IEQ), -1.e-99) 433 endif 434 enddo 435 WRITE(11,outform) ZH(IEQ, LL+1) / 1000., & 436 ZHbis(IEQ, LL+1) / 1000., & 437 (u_epflux_p(JW, IEQ), JW = 1, NW) 438 endif 439 440 end DO 441 442 ! 5 CALCUL DES TENDANCES: 443 !------------------------ 444 445 ! 5.1 Rectification des flux au sommet et dans les basses couches: 446 447 ! Attention, ici c'est le total sur toutes les ondes... 448 401 end DO 402 end DO ! DO LL = LAUNCH, nlayer - 1 403 404 !----------------------------------------------------------------------------------------------------------------------- 405 ! 5. TENDENCY CALCULATIONS 406 !----------------------------------------------------------------------------------------------------------------------- 407 408 ! 5.1 Flow rectification at the top and in the low layers 409 ! -------------------------------------------------------- 410 ! Warning, this is the total on all GW 449 411 u_epflux_tot(:, nlayer + 1) = 0. 450 412 v_epflux_tot(:, nlayer + 1) = 0. … … 458 420 end DO 459 421 DO LL = max(2,LAUNCH-2), LAUNCH-1 460 u_epflux_tot(:, LL) = u_epflux_tot(:, LL - 1) + u_epflux_tot(:, LAUNCH) * &461 (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))462 v_epflux_tot(:, LL) = v_epflux_tot(:, LL - 1) + v_epflux_tot(:, LAUNCH) * &463 (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))464 EAST_GWSTRESS(:,LL) = EAST_GWSTRESS(:, LL - 1) + &465 EAST_GWSTRESS(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1))/&466 (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))467 WEST_GWSTRESS(:,LL) = WEST_GWSTRESS(:, LL - 1) + &468 WEST_GWSTRESS(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1))/&469 (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))422 u_epflux_tot(:, LL) = u_epflux_tot(:, LL - 1) + u_epflux_tot(:, LAUNCH) * & 423 (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3))) 424 v_epflux_tot(:, LL) = v_epflux_tot(:, LL - 1) + v_epflux_tot(:, LAUNCH) * & 425 (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3))) 426 EAST_GWSTRESS(:,LL) = EAST_GWSTRESS(:, LL - 1) + & 427 EAST_GWSTRESS(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1))/ & 428 (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3))) 429 WEST_GWSTRESS(:,LL) = WEST_GWSTRESS(:, LL - 1) + & 430 WEST_GWSTRESS(:, LAUNCH) * (PH(:,LL)-PH(:,LL-1))/ & 431 (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3))) 470 432 end DO 471 433 ! This way, the total flux from GW is zero, but there is a net transport 472 434 ! (upward) that should be compensated by circulation 473 435 ! and induce additional friction at the surface 474 475 !Online output 476 if (output) then 477 DO LL = 1, nlayer - 1 478 WRITE(11,*) ZHbis(IEQ, LL)/1000.,u_epflux_tot(IEQ,LL) 479 end DO 480 CLOSE(11) 481 CALL abort_physic("nonoro_gwd_ran","stop here, the work is done",0) 482 endif 483 436 call writediagfi(ngrid,'nonoro_u_epflux_tot','nonoro_u_epflux_tot', '',3,u_epflux_tot(:,2:nlayer+1)) 437 call writediagfi(ngrid,'nonoro_v_epflux_tot','nonoro_v_epflux_tot', '',3,v_epflux_tot(:,2:nlayer+1)) 484 438 485 439 ! 5.2 AR-1 RECURSIVE FORMULA (13) IN VERSION 4 486 440 !--------------------------------------------- 487 441 DO LL = 1, nlayer 442 !Notice here the d_u and d_v are tendency (For Mars) but not increment(Venus). 488 443 d_u(:, LL) = G * (u_epflux_tot(:, LL + 1) - u_epflux_tot(:, LL)) & 489 / (PH(:, LL + 1) - PH(:, LL)) * DTIME444 / (PH(:, LL + 1) - PH(:, LL)) 490 445 d_v(:, LL) = G * (v_epflux_tot(:, LL + 1) - v_epflux_tot(:, LL)) & 491 / (PH(:, LL + 1) - PH(:, LL)) * DTIME446 / (PH(:, LL + 1) - PH(:, LL)) 492 447 ENDDO 493 448 d_t(:,:) = 0. 449 call writediagfi(ngrid,'nonoro_d_u','nonoro_d_u', '',3,d_u) 450 call writediagfi(ngrid,'nonoro_d_v','nonoro_d_v', '',3,d_v) 494 451 495 452 ! 5.3 Update tendency of wind with the previous (and saved) values 496 453 !----------------------------------------------------------------- 497 454 d_u(:,:) = DTIME/DELTAT/REAL(NW) * d_u(:,:) & 498 + (1.-DTIME/DELTAT) * du_nonoro_gwd(:,:)455 + (1.-DTIME/DELTAT) * du_nonoro_gwd(:,:) 499 456 d_v(:,:) = DTIME/DELTAT/REAL(NW) * d_v(:,:) & 500 + (1.-DTIME/DELTAT) * dv_nonoro_gwd(:,:)457 + (1.-DTIME/DELTAT) * dv_nonoro_gwd(:,:) 501 458 du_nonoro_gwd(:,:) = d_u(:,:) 502 459 dv_nonoro_gwd(:,:) = d_v(:,:) 503 460 461 call writediagfi(ngrid,'du_nonoro_gwd','du_nonoro_gwd', '',3,du_nonoro_gwd) 462 call writediagfi(ngrid,'dv_nonoro_gwd','dv_nonoro_gwd', '',3,dv_nonoro_gwd) 463 504 464 ! Cosmetic: evaluation of the cumulated stress 505 506 465 ZUSTR(:) = 0. 507 466 ZVSTR(:) = 0. 508 467 DO LL = 1, nlayer 509 ZUSTR(:) = ZUSTR(:) + D_U(:, LL) / g * (PH(:, LL + 1) - PH(:, LL))510 ZVSTR(:) = ZVSTR(:) + D_V(:, LL) / g * (PH(:, LL + 1) - PH(:, LL))468 ZUSTR(:) = ZUSTR(:) + D_U(:, LL) *rho(:, LL) * (ZH(:, LL + 1) - ZH(:, LL))*DTIME 469 ZVSTR(:) = ZVSTR(:) + D_V(:, LL) *rho(:, LL) * (ZH(:, LL + 1) - ZH(:, LL))*DTIME 511 470 ENDDO 512 471 472 !#ifdef CPP_XIOS 473 ! call send_xios_field("du_nonoro", d_u) 474 ! call send_xios_field("dv_nonoro", d_v) 475 ! call send_xios_field("east_gwstress", east_gwstress) 476 ! call send_xios_field("west_gwstress", west_gwstress) 477 !#endif 478 513 479 END SUBROUTINE NONORO_GWD_RAN 480 481 514 482 515 483 ! ======================================================== 516 484 ! Subroutines used to allocate/deallocate module variables 517 485 ! ======================================================== 518 519 486 SUBROUTINE ini_nonoro_gwd_ran(ngrid,nlayer) 520 487 -
trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
r2589 r2594 1524 1524 IF (calllott_nonoro) THEN 1525 1525 1526 CALL nonoro_gwd_ran(ngrid,nlayer,ptimestep,zplay, 1526 CALL nonoro_gwd_ran(ngrid,nlayer,ptimestep, 1527 & cpnew,rnew, 1528 & zplay, 1527 1529 & zmax_th, ! max altitude reached by thermals (m) 1528 1530 & pt, pu, pv, … … 1534 1536 pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer) 1535 1537 & +d_t_hin(1:ngrid,1:nlayer) 1536 ! d_t_hin(:,:)= d_t_hin(:,:)/ptimestep ! K/s (for outputs?)1537 1538 pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer) 1538 1539 & +d_u_hin(1:ngrid,1:nlayer) 1539 ! d_u_hin(:,:)= d_u_hin(:,:)/ptimestep ! (m/s)/s (for outputs?)1540 1540 pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer) 1541 1541 & +d_v_hin(1:ngrid,1:nlayer) 1542 ! d_v_hin(:,:)= d_v_hin(:,:)/ptimestep ! (m/s)/s (for outputs?)1543 1542 1544 1543 ENDIF ! of IF (calllott_nonoro)
Note: See TracChangeset
for help on using the changeset viewer.
