Changeset 2232
- Timestamp:
- Jan 30, 2020, 10:36:35 AM (5 years ago)
- Location:
- trunk/LMDZ.GENERIC
- Files:
-
- 5 edited
-
README (modified) (1 diff)
-
libf/phystd/convadj.F (modified) (5 diffs)
-
libf/phystd/physiq_mod.F90 (modified) (8 diffs)
-
libf/phystd/thermcell_main.F90 (modified) (16 diffs)
-
libf/phystd/thermcell_plume.F90 (modified) (18 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/README
r2231 r2232 1506 1506 - fix a bug in thermcell_env. Now zqt is correctly initialized when tracer h2o_vap is missing (consistency with flag water is assumed). 1507 1507 - clean up thermcell_dv2 and use plume height (zmin-zmax) instead of maximal altitude (zmax) in computations 1508 - the thermal plume model is able to manage several plumes in the same column and work without the convective adjustment. -
trunk/LMDZ.GENERIC/libf/phystd/convadj.F
r2127 r2232 3 3 & pu,pv,ph,pq, 4 4 & pdufi,pdvfi,pdhfi,pdqfi, 5 & pduadj,pdvadj,pdhadj, 6 & pdqadj,lmax) 5 & pduadj,pdvadj,pdhadj,pdqadj) 7 6 8 7 USE tracer_h … … 41 40 REAL pu(ngrid,nlay),pdufi(ngrid,nlay),pduadj(ngrid,nlay) 42 41 REAL pv(ngrid,nlay),pdvfi(ngrid,nlay),pdvadj(ngrid,nlay) 43 INTEGER lmax(ngrid)44 42 45 43 ! Tracers … … 157 155 DO l=2,nlay 158 156 DO ig=1,ngrid 159 IF (zhc(ig,l).LT.zhc(ig,l-1) .and.(l.GT.lmax(ig))) THEN157 IF (zhc(ig,l).LT.zhc(ig,l-1)) THEN 160 158 vtest(ig)=.true. 161 159 ENDIF … … 197 195 l2 = l2 + 1 198 196 IF (l2 .GT. nlay) EXIT 199 IF ( (zhc(i, l2).LT.zhc(i, l2-1)).and.(l2.GT.lmax(i))) THEN197 IF (zhc(i, l2).LT.zhc(i, l2-1)) THEN 200 198 201 199 ! l2 is the highest level of the unstable column … … 226 224 down = .false. 227 225 IF (l1 .ne. 1) then !-- and then 228 IF ( (zhmc.LT.zhc(i, l1-1)).and.(l1.GT.lmax(i))) then226 IF (zhmc.LT.zhc(i, l1-1)) then 229 227 down = .true. 230 228 END IF -
trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90
r2176 r2232 253 253 ! VARIABLES for the thermal plume model 254 254 255 INTEGER lmin(ngrid) ! Plume base level 256 INTEGER lmix(ngrid) ! Plume maximal velocity level 257 INTEGER lmax(ngrid) ! Maximal level reached by the plume 258 259 real,allocatable,save :: f0(:) ! Mass flux norm 260 real,allocatable,save :: fm0(:,:) ! Mass flux 261 real,allocatable,save :: entr0(:,:) ! Entrainment 262 real,allocatable,save :: detr0(:,:) ! Detrainment 255 real f(ngrid) ! Mass flux norm 256 real fm(ngrid,nlayer+1) ! Mass flux 257 real fm_bis(ngrid,nlayer) ! Recasted fm 258 real entr(ngrid,nlayer) ! Entrainment 259 real detr(ngrid,nlayer) ! Detrainment 263 260 real dqevap(ngrid,nlayer,nq) ! water tracer mass mixing ratio variations due to evaporation 264 261 real dtevap(ngrid,nlayer) ! temperature variation due to evaporation 265 262 real zqtherm(ngrid,nlayer,nq) ! vapor mass mixing ratio after evaporation 266 263 real zttherm(ngrid,nlayer) ! temperature after evaporation 267 real fraca(ngrid, nlayer+1) ! Fraction of the surface that plumes occupies 268 real zw2(ngrid, nlayer+1) ! Vertical speed 269 real zw2_bis(ngrid, nlayer) ! 270 real zqsatth(ngrid, nlayer) ! Water vapor mixing ratio at saturation 271 real zqta(ngrid, nlayer) ! Total water mass mixing ratio in the plume 272 real zqla(ngrid, nlayer) ! Liquid water mass mixing ratio in the plume 273 real ztv(ngrid, nlayer) ! Virtual potential temperature in the environment considering large scale condensation 274 real ztva(ngrid, nlayer) ! Virtual potential temperature in the plume 275 real zthl(ngrid, nlayer) ! Potential temperature in the environment without considering condensation 276 real ztla(ngrid, nlayer) ! Potential temperature in the plume 277 real ratqscth(ngrid, nlayer) ! 278 real ratqsdiff(ngrid, nlayer) ! 279 280 ! AB : variables used for outputs 281 REAL pmax(ngrid) ! Maximal pressure reached by the plume 282 REAL pmin(ngrid) ! Minimal pressure reached by the plume 283 264 real fraca(ngrid,nlayer+1) ! Fraction of the surface that plumes occupies 265 real zw2(ngrid,nlayer+1) ! Vertical speed 266 real zw2_bis(ngrid,nlayer) ! Recasted zw2 267 284 268 ! TENDENCIES due to various processes : 285 269 … … 688 672 if (calltherm) then 689 673 CALL init_thermcell_mod(g, rcp, r, pi, T_h2o_ice_liq, RV) 690 ALLOCATE(f0(ngrid))691 ALLOCATE(fm0(ngrid,nlayer+1))692 ALLOCATE(entr0(ngrid,nlayer))693 ALLOCATE(detr0(ngrid,nlayer))694 674 endif 695 675 … … 1188 1168 ENDIF 1189 1169 1190 ! AB: If a plume stops, the parametrization never look above if somewhere the atmosphere is still unstable!1191 ! Maybe it's a good idea to call convective adjustment too.1192 1170 CALL thermcell_main(ngrid, nlayer, nq, ptimestep, firstcall, & 1193 1171 pplay, pplev, pphi, zpopsk, & 1194 1172 pu, pv, zttherm, zqtherm, & 1195 1173 zdutherm, zdvtherm, zdttherm, zdqtherm, & 1196 f0, fm0, entr0, detr0, zw2, fraca, & 1197 zqta, zqla, ztv, ztva, ztla, zthl, zqsatth, & 1198 lmin, lmix, lmax) 1199 1200 DO ig=1,ngrid 1201 pmin(ig) = pplev(ig,lmin(ig)) 1202 pmax(ig) = pplev(ig,lmax(ig)) 1203 ENDDO 1174 fm, entr, detr, zw2, fraca) 1204 1175 1205 1176 pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zdutherm(1:ngrid,1:nlayer) … … 1210 1181 pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqtherm(1:ngrid,1:nlayer,1:nq) 1211 1182 ENDIF 1212 1213 ELSE1214 1215 lmax(1:ngrid) = 11216 1183 1217 1184 ENDIF ! end of 'calltherm' … … 1234 1201 pdu,pdv,zdh,pdq, & 1235 1202 zduadj,zdvadj,zdhadj, & 1236 zdqadj ,lmax)1203 zdqadj) 1237 1204 1238 1205 pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer) + zduadj(1:ngrid,1:nlayer) … … 1748 1715 DO l=1,nlayer 1749 1716 zw2_bis(ig,l) = zw2(ig,l) 1717 fm_bis(ig,l) = fm(ig,l) 1750 1718 ENDDO 1751 1719 ENDDO … … 2138 2106 ! Thermal plume model 2139 2107 if (calltherm) then 2140 call writediagfi(ngrid,'entr','Entrainment','kg m$^{-2}$ s$^{-1}$', 3, entr 0)2141 call writediagfi(ngrid,'detr','Detrainment','kg m$^{-2}$ s$^{-1}$', 3, detr 0)2142 call writediagfi(ngrid,'fm','Mass flux','kg m$^{-2}$ s$^{-1}$', 3, fm 0)2108 call writediagfi(ngrid,'entr','Entrainment','kg m$^{-2}$ s$^{-1}$', 3, entr) 2109 call writediagfi(ngrid,'detr','Detrainment','kg m$^{-2}$ s$^{-1}$', 3, detr) 2110 call writediagfi(ngrid,'fm','Mass flux','kg m$^{-2}$ s$^{-1}$', 3, fm_bis) 2143 2111 call writediagfi(ngrid,'w_plm','Squared vertical velocity','m s$^{-1}$', 3, zw2_bis) 2144 2112 call writediagfi(ngrid,'fraca','Updraft fraction','', 3, fraca) 2145 ! call writediagfi(ngrid,'pmin', 'pmin', 'Pa', 2, pmin)2146 ! call writediagfi(ngrid,'pmax', 'pmax', 'Pa', 2, pmax)2147 2113 endif 2148 2114 … … 2352 2318 IF (calltherm) THEN 2353 2319 CALL send_xios_field('w_plm',zw2_bis) 2354 CALL send_xios_field('entr',entr0) 2355 CALL send_xios_field('detr',detr0) 2320 CALL send_xios_field('entr',entr) 2321 CALL send_xios_field('detr',detr) 2322 ! CALL send_xios_field('fm',fm_bis) 2323 ! CALL send_xios_field('fraca',fraca) 2356 2324 ENDIF 2357 2325 -
trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90
r2229 r2232 6 6 pu,pv,pt,pq, & 7 7 pduadj,pdvadj,pdtadj,pdqadj, & 8 f0,fm0,entr0,detr0,zw2,fraca, & 9 zqta,zqla,ztv,ztva,zhla,zhl,zqsa, & 10 lmin,lmix,lmax) 8 fm_tot,entr_tot,detr_tot,zw2_tot,fraca) 11 9 12 10 … … 43 41 ! lmin can be greater than 1 44 42 ! Mix every tracer 43 ! Can stack verticaly multiple plumes (it makes thermcell_dv2 unusable for the moment) 45 44 ! 46 45 !=============================================================================== … … 67 66 REAL, INTENT(in) :: pplev(ngrid,nlay+1) ! Level pressure 68 67 REAL, INTENT(in) :: pphi(ngrid,nlay) ! Geopotential 68 REAL, INTENT(in) :: zpopsk(ngrid,nlay) ! Exner function 69 69 70 70 REAL, INTENT(in) :: pu(ngrid,nlay) ! Zonal wind … … 78 78 ! -------- 79 79 80 INTEGER, INTENT(out) :: lmax(ngrid) ! Highest layer reached by the plume81 INTEGER, INTENT(out) :: lmix(ngrid) ! Layer in which plume vertical speed is maximal82 INTEGER, INTENT(out) :: lmin(ngrid) ! First unstable layer83 84 80 REAL, INTENT(out) :: pduadj(ngrid,nlay) ! u convective variations 85 81 REAL, INTENT(out) :: pdvadj(ngrid,nlay) ! v convective variations … … 87 83 REAL, INTENT(out) :: pdqadj(ngrid,nlay,nq) ! q convective variations 88 84 89 REAL, INTENT(inout) :: f0(ngrid) ! mass flux norm (after possible time relaxation) 90 REAL, INTENT(inout) :: fm0(ngrid,nlay+1) ! mass flux (after possible time relaxation) 91 REAL, INTENT(inout) :: entr0(ngrid,nlay) ! entrainment (after possible time relaxation) 92 REAL, INTENT(inout) :: detr0(ngrid,nlay) ! detrainment (after possible time relaxation) 85 REAL, INTENT(inout) :: fm_tot(ngrid,nlay+1) ! Total mass flux 86 REAL, INTENT(inout) :: entr_tot(ngrid,nlay) ! Total entrainment 87 REAL, INTENT(inout) :: detr_tot(ngrid,nlay) ! Total detrainment 88 89 REAL, INTENT(out) :: fraca(ngrid,nlay+1) ! Updraft fraction 90 REAL, INTENT(out) :: zw2_tot(ngrid,nlay+1) ! Total plume vertical speed 93 91 94 92 ! Local: … … 96 94 97 95 INTEGER ig, k, l, iq 96 INTEGER lmax(ngrid) ! Highest layer reached by the plume 97 INTEGER lmix(ngrid) ! Layer in which plume vertical speed is maximal 98 INTEGER lmin(ngrid) ! First unstable layer 98 99 99 100 REAL zmix(ngrid) ! Altitude of maximal vertical speed … … 106 107 REAL rhobarz(ngrid,nlay) ! Levels densities 107 108 REAL masse(ngrid,nlay) ! Layers masses 108 REAL zpopsk(ngrid,nlay) ! Exner function109 109 110 110 REAL zu(ngrid,nlay) ! u environment … … 128 128 129 129 REAL f_star(ngrid,nlay+1) ! Normalized mass flux 130 REAL entr_star(ngrid,nlay) ! Normalized entrainment (E* = e* dz) 131 REAL detr_star(ngrid,nlay) ! Normalized detrainment (D* = d* dz) 132 130 REAL entr_star(ngrid,nlay) ! Normalized entrainment 131 REAL detr_star(ngrid,nlay) ! Normalized detrainment 132 133 REAL f(ngrid) ! Mass flux norm 133 134 REAL fm(ngrid,nlay+1) ! Mass flux 134 REAL entr(ngrid,nlay) ! Entrainment (E = e dz) 135 REAL detr(ngrid,nlay) ! Detrainment (D = d dz) 136 137 REAL f(ngrid) ! Mass flux norm 138 REAL lambda ! Time relaxation coefficent 139 REAL fraca(ngrid,nlay+1) ! Updraft fraction 135 REAL entr(ngrid,nlay) ! Entrainment 136 REAL detr(ngrid,nlay) ! Detrainment 137 138 REAL zw2(ngrid,nlay+1) ! Plume vertical speed 140 139 REAL wmax(ngrid) ! Maximal vertical speed 141 REAL zw2(ngrid,nlay+1) ! Plume vertical speed142 140 REAL zdthladj(ngrid,nlay) ! Potential temperature variations 143 141 REAL dummy(ngrid,nlay) ! Dummy argument for thermcell_dq() 144 142 143 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 144 INTEGER lbot(ngrid) 145 LOGICAL re_tpm 146 INTEGER while_loop_counter 147 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 148 145 149 !=============================================================================== 146 150 ! Initialization 147 151 !=============================================================================== 148 152 149 IF (firstcall) THEN 150 fm0(:,:) = 0. 151 entr0(:,:) = 0. 152 detr0(:,:) = 0. 153 ENDIF 154 155 DO ig=1,ngrid 156 ! AB: Minimal f0 value is set to 0. (instead of 1.e-2 in Earth version) 157 f0(ig) = MAX(f0(ig), 0.) 158 ENDDO 153 fm_tot(:,:) = 0. 154 entr_tot(:,:) = 0. 155 detr_tot(:,:) = 0. 156 zw2_tot(:,:) = 0. 159 157 160 158 pduadj(:,:) = 0.0 … … 164 162 165 163 zdthladj(:,:) = 0.0 164 165 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 166 re_tpm = .true. 167 lbot(:) = linf 168 while_loop_counter = 0. 169 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 166 170 167 171 !=============================================================================== … … 243 247 ! --------------------------- | 244 248 ! | 245 ! --------------------------- rhobarz, f_star, fm, fm0,zw2, fraca249 ! --------------------------- rhobarz, f_star, fm, zw2, fraca 246 250 ! zt, zu, zv, zo, rho | 247 251 ! --------------------------- | … … 287 291 !=============================================================================== 288 292 289 !------------------------------------------------------------------------------- 290 ! Thermal plumes speeds, fluxes, tracers and temperatures 291 !------------------------------------------------------------------------------- 292 293 CALL thermcell_plume(ngrid,nlay,nq,ptimestep, & 294 & ztv,zhl,zqt,zql,zlev,pplev,zpopsk, & 295 & detr_star,entr_star,f_star, & 296 & ztva,zhla,zqta,zqla,zqsa, & 297 & zw2,lmin) 298 293 DO WHILE (re_tpm.and.(while_loop_counter<nlay)) 294 while_loop_counter = while_loop_counter + 1 295 296 !------------------------------------------------------------------------------- 297 ! Thermal plumes speeds, normalized fluxes, tracers and temperatures 298 !------------------------------------------------------------------------------- 299 300 CALL thermcell_plume(ngrid,nlay,nq,ptimestep, & 301 & ztv,zhl,zqt,zql,zlev,pplev,zpopsk, & 302 & detr_star,entr_star,f_star, & 303 & ztva,zhla,zqta,zqla,zqsa, & 304 & zw2,lbot,lmin) 305 299 306 !------------------------------------------------------------------------------- 300 307 ! Thermal plumes characteristics: zmax, zmix, wmax 301 308 !------------------------------------------------------------------------------- 302 309 303 310 ! AB: Careful, zw2 became its square root in thermcell_height! 304 CALL thermcell_height(ngrid,nlay,lmin,lmix,lmax,zlev, & 305 & zmin,zmix,zmax,zw2,wmax,f_star) 306 307 !=============================================================================== 308 ! Closure and mass fluxes computation 309 !=============================================================================== 310 311 CALL thermcell_height(ngrid,nlay,lmin,lmix,lmax, & 312 & zlev,zmin,zmix,zmax,zw2,wmax,f_star) 313 311 314 !------------------------------------------------------------------------------- 312 315 ! Closure 313 316 !------------------------------------------------------------------------------- 314 315 CALL thermcell_closure(ngrid,nlay,ptimestep,rho,zlev, & 316 & lmax,entr_star,zmin,zmax,wmax,f) 317 318 IF (tau_thermals>1.) THEN 319 lambda = exp(-ptimestep/tau_thermals) 320 f0(:) = (1.-lambda) * f(:) + lambda * f0(:) 321 ELSE 322 f0(:) = f(:) 323 ENDIF 324 317 318 CALL thermcell_closure(ngrid,nlay,ptimestep,rho,zlev, & 319 & lmax,entr_star,zmin,zmax,wmax,f) 320 325 321 ! FH: Test valable seulement en 1D mais pas genant 326 IF (.not. (f0(1).ge.0.) ) THEN327 print *, 'ERROR: mass flux norm is not positive!'328 print *, 'f0 =', f0(1)329 CALL abort330 ENDIF331 322 IF (.not. (f(1).ge.0.) ) THEN 323 print *, 'ERROR: mass flux norm is not positive!' 324 print *, 'f =', f(1) 325 CALL abort 326 ENDIF 327 332 328 !------------------------------------------------------------------------------- 333 329 ! Mass fluxes 334 330 !------------------------------------------------------------------------------- 335 336 CALL thermcell_flux(ngrid,nlay,ptimestep,masse, & 337 & lmin,lmax,entr_star,detr_star, & 338 & f,rhobarz,zlev,zw2,fm,entr,detr,zqla) 339 340 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 341 ! On ne prend pas directement les profils issus des calculs precedents mais on 342 ! s'autorise genereusement une relaxation vers ceci avec une constante de temps 343 ! tau_thermals (typiquement 1800s sur Terre). 344 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 345 346 IF (tau_thermals>1.) THEN 347 lambda = exp(-ptimestep/tau_thermals) 348 fm0 = (1.-lambda) * fm + lambda * fm0 349 entr0 = (1.-lambda) * entr + lambda * entr0 350 detr0 = (1.-lambda) * detr + lambda * detr0 351 ELSE 352 fm0(:,:) = fm(:,:) 353 entr0(:,:) = entr(:,:) 354 detr0(:,:) = detr(:,:) 355 ENDIF 356 357 !------------------------------------------------------------------------------- 331 332 CALL thermcell_flux(ngrid,nlay,ptimestep,masse, & 333 & lmin,lmax,entr_star,detr_star, & 334 & f,rhobarz,zlev,zw2,fm,entr,detr) 335 336 !------------------------------------------------------------------------------- 337 ! 338 !------------------------------------------------------------------------------- 339 340 re_tpm = .false. 341 DO ig=1,ngrid 342 IF(lmax(ig) > lmin(ig)) THEN 343 lbot(ig) = lmax(ig) 344 re_tpm = .true. 345 ELSE 346 lbot(ig) = nlay 347 ENDIF 348 ENDDO 349 350 !------------------------------------------------------------------------------- 351 ! Thermal plumes stacking 352 !------------------------------------------------------------------------------- 353 354 zw2_tot(:,:) = zw2_tot(:,:) + zw2(:,:) 355 entr_tot(:,:) = entr_tot(:,:) + entr(:,:) 356 detr_tot(:,:) = detr_tot(:,:) + detr(:,:) 357 fm_tot(:,:) = fm_tot(:,:) + fm(:,:) 358 359 ENDDO 360 361 !=============================================================================== 358 362 ! Updraft fraction 359 ! -------------------------------------------------------------------------------363 !=============================================================================== 360 364 361 365 DO ig=1,ngrid … … 366 370 DO l=2,nlay 367 371 DO ig=1,ngrid 368 IF (zw2 (ig,l) > 0.) THEN369 fraca(ig,l) = fm (ig,l) / (rhobarz(ig,l) * zw2(ig,l))372 IF (zw2_tot(ig,l) > 0.) THEN 373 fraca(ig,l) = fm_tot(ig,l) / (rhobarz(ig,l) * zw2_tot(ig,l)) 370 374 ELSE 371 375 fraca(ig,l) = 0. … … 375 379 376 380 !=============================================================================== 377 ! Transport vertical381 ! Vertical transport 378 382 !=============================================================================== 379 383 … … 382 386 !------------------------------------------------------------------------------- 383 387 384 CALL thermcell_dq(ngrid,nlay,ptimestep,fm 0,entr0,detr0,masse,&385 & zhl,zdthladj,dummy)388 CALL thermcell_dq(ngrid,nlay,ptimestep,fm_tot,entr_tot,detr_tot, & 389 & masse,zhl,zdthladj,dummy) 386 390 387 391 DO l=1,nlay … … 396 400 397 401 DO iq=1,nq 398 CALL thermcell_dq(ngrid,nlay,ptimestep,fm 0,entr0,detr0,masse,&399 & pq(:,:,iq),pdqadj(:,:,iq),zqa(:,:,iq))402 CALL thermcell_dq(ngrid,nlay,ptimestep,fm_tot,entr_tot,detr_tot, & 403 & masse,pq(:,:,iq),pdqadj(:,:,iq),zqa(:,:,iq)) 400 404 ENDDO 401 405 … … 404 408 !------------------------------------------------------------------------------- 405 409 410 ! AB: Careful, thermcell_dv2 wasn't checked! It is not sure that it works 411 ! correctly with the plumes stacking (zmin, wmax doesn't make sense). 406 412 IF (dvimpl) THEN 407 CALL thermcell_dv2(ngrid,nlay,ptimestep,fm 0,entr0,detr0,masse,fraca, &408 & zmax,zmin,pu,pv,pduadj,pdvadj,zua,zva)413 CALL thermcell_dv2(ngrid,nlay,ptimestep,fm_tot,entr_tot,detr_tot, & 414 & masse,fraca,zmax,zmin,pu,pv,pduadj,pdvadj,zua,zva) 409 415 ELSE 410 CALL thermcell_dq(ngrid,nlay,ptimestep,fm 0,entr0,detr0,masse,&411 & zu,pduadj,zua)412 CALL thermcell_dq(ngrid,nlay,ptimestep,fm 0,entr0,detr0,masse,&413 & zv,pdvadj,zva)416 CALL thermcell_dq(ngrid,nlay,ptimestep,fm_tot,entr_tot,detr_tot, & 417 & masse,zu,pduadj,zua) 418 CALL thermcell_dq(ngrid,nlay,ptimestep,fm_tot,entr_tot,detr_tot, & 419 & masse,zv,pdvadj,zva) 414 420 ENDIF 415 421 -
trunk/LMDZ.GENERIC/libf/phystd/thermcell_plume.F90
r2178 r2232 6 6 detr_star,entr_star,f_star, & 7 7 ztva,zhla,zqta,zqla,zqsa, & 8 zw2,l min)8 zw2,lbot,lmin) 9 9 10 10 … … 41 41 INTEGER, INTENT(in) :: nq 42 42 43 INTEGER, INTENT(in) :: lbot(ngrid) ! First considered layer 44 43 45 REAL, INTENT(in) :: ptimestep 44 46 REAL, INTENT(in) :: zlev(ngrid,nlay+1) ! Levels altitude … … 65 67 REAL, INTENT(out) :: zqta(ngrid,nlay) ! qt plume (after mixing) 66 68 REAL, INTENT(out) :: zqsa(ngrid,nlay) ! qsat plume (after mixing) 67 REAL, INTENT(out) :: zw2(ngrid,nlay+1) ! w plum r(after mixing)69 REAL, INTENT(out) :: zw2(ngrid,nlay+1) ! w plume (after mixing) 68 70 69 71 ! Local: … … 71 73 72 74 INTEGER ig, l, k 75 INTEGER l_start 73 76 74 77 REAL ztva_est(ngrid,nlay) ! TRPV plume (before mixing) … … 81 84 82 85 REAL zbuoy(ngrid,nlay) ! Plume buoyancy 83 REAL ztemp(ngrid) ! Temperature for saturation vapor pressure computation in plume86 REAL ztemp(ngrid) ! Temperature to compute saturation vapor pressure 84 87 REAL zdz ! Layers heights 85 88 REAL ztv2(ngrid,nlay) ! ztv + d_temp * Dirac(l=linf) … … 93 96 REAL psat ! Dummy argument for Psat_water() 94 97 95 ! REAL alpha 96 97 LOGICAL active(ngrid) ! If the plume is active at ig (speed and incoming mass flux > 0) 98 LOGICAL activetmp(ngrid) ! If the plume is active at ig (active=true and outgoing mass flux > 0) 98 LOGICAL active(ngrid) ! If the plume is active (speed and incoming mass flux > 0) 99 LOGICAL activetmp(ngrid) ! If the plume is active (active=true and outgoing mass flux > 0) 99 100 100 101 !=============================================================================== … … 102 103 !=============================================================================== 103 104 104 ztva(:,:) = ztv(:,:) ! ztva is set to TPV environment 105 zhla(:,:) = zhl(:,:) ! zhla is set to TP environment 106 zqta(:,:) = zqt(:,:) ! zqta is set to qt environment 107 zqla(:,:) = zql(:,:) ! zqla is set to ql environment 108 109 zqsa_est(:) = 0. 110 zqsa(:,:) = 0. 111 112 zw2_est(:,:) = 0. 113 zw2(:,:) = 0. 105 ztva(:,:) = ztv(:,:) ! ztva is set to TPV environment 106 zhla(:,:) = zhl(:,:) ! zhla is set to TP environment 107 zqta(:,:) = zqt(:,:) ! zqta is set to qt environment 108 zqla(:,:) = zql(:,:) ! zqla is set to ql environment 109 zqsa(:,:) = 0. 110 zw2(:,:) = 0. 111 112 ztva_est(:,:) = ztv(:,:) ! ztva_est is set to TPV environment 113 zqla_est(:,:) = zql(:,:) ! zqla_est is set to ql environment 114 zqsa_est(:) = 0. 115 zw2_est(:,:) = 0. 114 116 115 117 zbuoy(:,:) = 0. … … 119 121 entr_star(:,:) = 0. 120 122 121 lmin(:) = 1123 lmin(:) = lbot(:) 122 124 123 125 ztv2(:,:) = ztv(:,:) … … 126 128 active(:) = .false. 127 129 130 l_start = nlay 131 128 132 !=============================================================================== 129 133 ! First layer computation … … 131 135 132 136 DO ig=1,ngrid 133 l = linf 137 l = lbot(ig) 138 l_start = MIN(l_start, lbot(ig)+1) 134 139 DO WHILE (.not.active(ig).and.(pplev(ig,l+1) > pres_limit).and.(l < nlay)) 135 140 zbuoy(ig,l) = RG * (ztv2(ig,l) - ztv2(ig,l+1)) / ztv2(ig,l+1) … … 157 162 !=============================================================================== 158 163 159 DO l=l inf+1,nlay-1164 DO l=l_start,nlay-1 160 165 161 166 !------------------------------------------------------------------------------- … … 164 169 165 170 DO ig=1,ngrid 166 active(ig) = (zw2(ig,l) > 1.e- 10).and.(f_star(ig,l) > 1.e-10)171 active(ig) = (zw2(ig,l) > 1.e-9).and.(f_star(ig,l) > 1.e-9) 167 172 ENDDO 168 173 … … 185 190 DO ig=1,ngrid 186 191 IF (active(ig)) THEN 187 zqla_est(ig,l) = MAX(0.,zqta(ig,l-1) - zqsa_est(ig)) ! zqla_est set to ql plume188 zta_est(ig,l) = zhla(ig,l-1) * zpopsk(ig,l) & ! zta_est set to TR plume192 zqla_est(ig,l) = MAX(0.,zqta(ig,l-1) - zqsa_est(ig)) ! zqla_est is set to ql plume 193 zta_est(ig,l) = zhla(ig,l-1) * zpopsk(ig,l) & ! zta_est is set to TR plume 189 194 & + RLvCp * zqla_est(ig,l) 190 ztva_est(ig,l) = zta_est(ig,l) / zpopsk(ig,l) & ! ztva_est set to TRPV plume195 ztva_est(ig,l) = zta_est(ig,l) / zpopsk(ig,l) & ! ztva_est is set to TRPV plume 191 196 & * (1. + RETV * (zqta(ig,l-1)-zqla_est(ig,l)) - zqla_est(ig,l)) 192 197 … … 196 201 zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha) 197 202 zdw2 = afact * zbuoy(ig,l) / fact_epsilon 198 zw2_est(ig,l+1) = M ax(0., exp(-zw2fact) * (zw2_est(ig,l) - zdw2) + zdw2)203 zw2_est(ig,l+1) = MAX(0., exp(-zw2fact) * (zw2_est(ig,l) - zdw2) + zdw2) 199 204 ENDIF 200 205 ENDDO … … 208 213 209 214 zdz = zlev(ig,l+1) - zlev(ig,l) 210 zw2m = (zw2_est(ig,l+1) + zw2 _est(ig,l)) / 2.215 zw2m = (zw2_est(ig,l+1) + zw2(ig,l)) / 2. 211 216 gamma = afact * zbuoy(ig,l) - fact_epsilon * zw2m 212 217 213 IF (zw2 _est(ig,l)> 0.) THEN218 IF (zw2m > 0.) THEN 214 219 test = gamma / zw2m - nu 215 220 ELSE 216 print *, 'ERROR: zw2_est is negative while plume is active!' 221 test = 0. 222 print *, 'WARNING: vertical speed is negative while plume is active!' 217 223 print *, 'ig,l', ig, l 218 print *, 'zw2_est', zw2_est(ig,l) 219 call abort 224 print *, 'zw2m', zw2m 220 225 ENDIF 221 226 … … 237 242 !------------------------------------------------------------------------------- 238 243 239 activetmp(:) = active(:).and.(f_star(:,l+1) > 1.e- 10)244 activetmp(:) = active(:).and.(f_star(:,l+1) > 1.e-9) 240 245 241 246 DO ig=1,ngrid … … 271 276 zta(ig,l) = zhla(ig,l) * zpopsk(ig,l) & ! ztva is set to TR plume (mixed) 272 277 & + RLvCp * zqla(ig,l) 273 ztva(ig,l) = zta(ig,l) / zpopsk(ig,l) & ! ztva is set to TRPV plume (mixed) 278 ztva(ig,l) = zta(ig,l) / zpopsk(ig,l) & ! ztva is set to TRPV plume (mixed) 274 279 & * (1. + RETV*(zqta(ig,l)-zqla(ig,l)) - zqla(ig,l)) 275 280 … … 279 284 zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha) 280 285 zdw2 = afact * zbuoy(ig,l) / fact_epsilon 281 zw2(ig,l+1) = M ax(0., exp(-zw2fact) * (zw2(ig,l) - zdw2) + zdw2)286 zw2(ig,l+1) = MAX(0., exp(-zw2fact) * (zw2(ig,l) - zdw2) + zdw2) 282 287 ENDIF 283 288 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
