! ! ! SUBROUTINE thermcell_plume(ngrid,nlay,nq,ptimestep,ztv, & zhl,zqt,zql,rhobarz,zlev,pplev,pphi,zpopsk, & detr_star,entr_star,f_star, & ztva,zhla,zqla,zqta,zta, & zw2,zqsa,lmix,lmin) !=============================================================================== ! Purpose: calcule les valeurs de qt, thetal et w dans l ascendance ! ! Nota Bene ! ql means "non-gaseous water mass mixing ratio" (liquid and solid) ! qv means "vapor mass mixing ratio" ! qt means "total water mass mixing ratio" ! TP means "potential temperature" ! TRPV means "virtual potential temperature with latent heat release" ! TPV means "virtual potential temperature" ! TR means "temperature with latent heat release" !=============================================================================== USE print_control_mod, ONLY: prt_level USE watercommon_h, ONLY: RLvCp, RETV, Psat_water USE tracer_h, ONLY: igcm_h2o_vap USE thermcell_mod IMPLICIT NONE !=============================================================================== ! Declaration !=============================================================================== ! Inputs: ! ------- INTEGER ngrid, nlay, nq REAL ptimestep REAL rhobarz(ngrid,nlay) ! Levels density REAL zlev(ngrid,nlay+1) ! Levels altitude REAL pplev(ngrid,nlay+1) ! Levels pressure REAL pphi(ngrid,nlay) ! Geopotential REAL zpopsk(ngrid,nlay) ! Exner function REAL ztv(ngrid,nlay) ! TRPV environment REAL zhl(ngrid,nlay) ! TP environment REAL zqt(ngrid,nlay) ! qt environment REAL zql(ngrid,nlay) ! ql environment ! Outputs: ! -------- INTEGER lmin(ngrid) ! plume base level (first unstable level) INTEGER lmix(ngrid) ! maximum vertical speed level REAL detr_star(ngrid,nlay) ! normalized detrainment REAL entr_star(ngrid,nlay) ! normalized entrainment REAL f_star(ngrid,nlay+1) ! normalized mass flux REAL ztva(ngrid,nlay) ! TRPV plume (after mixing) REAL zhla(ngrid,nlay) ! TP plume (after mixing) REAL zqla(ngrid,nlay) ! ql plume (after mixing) REAL zqta(ngrid,nlay) ! qt plume (after mixing) REAL zqsa(ngrid,nlay) ! qsat plume (after mixing) REAL zw2(ngrid,nlay+1) ! w2 plume (after mixing) ! Local: ! ------ INTEGER ig, l, k REAL ztva_est(ngrid,nlay) ! TRPV plume (before mixing) REAL zqla_est(ngrid,nlay) ! ql plume (before mixing) REAL zta_est(ngrid,nlay) ! TR plume (before mixing) REAL zqsa_est(ngrid) ! qsat plume (before mixing) REAL zw2_est(ngrid,nlay+1) ! w2 plume (before mixing) REAL zta(ngrid,nlay) ! TR plume (after mixing) REAL zbuoy(ngrid,nlay) ! Plume buoyancy REAL ztemp(ngrid) ! Temperature for saturation vapor pressure computation in plume REAL zdz ! Layers heights REAL ztv2(ngrid,nlay) ! ztv + d_temp * Dirac(l=linf) REAL zbetalpha ! REAL zdw2 ! REAL zdw2bis ! REAL zw2fact ! REAL zw2m ! Average vertical velocity between two successive levels REAL gamma ! Plume acceleration term (to compute vertical velocity) REAL test ! Test to know how to compute entrainment and detrainment REAL psat ! Dummy argument for Psat_water() LOGICAL active(ngrid) ! If the plume is active at ig (speed and incoming mass flux > 0 or l=lmin) LOGICAL activetmp(ngrid) ! If the plume is active at ig (active=true and outgoing mass flux > 0) !=============================================================================== ! Initialization !=============================================================================== zbetalpha = betalpha / (1. + betalpha) ztva(:,:) = ztv(:,:) ! ztva is set to TPV environment zhla(:,:) = zhl(:,:) ! zhla is set to TP environment zqta(:,:) = zqt(:,:) ! zqta is set to qt environment zqla(:,:) = zql(:,:) ! zqla is set to ql environment zqsa_est(:) = 0. zqsa(:,:) = 0. zw2_est(:,:) = 0. zw2(:,:) = 0. zbuoy(:,:) = 0. f_star(:,:) = 0. detr_star(:,:) = 0. entr_star(:,:) = 0. lmix(:) = 1 lmin(:) = 1 ztv2(:,:) = ztv(:,:) ztv2(:,linf) = ztv(:,linf) + d_temp !=============================================================================== ! First layer computation !=============================================================================== DO ig=1,ngrid active(ig) = .false. l = linf DO WHILE (.not.active(ig).and.(pplev(ig,l+1).GT.pres_limit).and.(l.LT.nlay)) zdz = (zlev(ig,l+1) - zlev(ig,l)) zw2(ig,l+1) = 2. * afact * RG * zdz & ! Do we have to divide by 1+betalpha (consider entrainment) ? & * (ztv2(ig,l) - ztv2(ig,l+1)) & & / (ztv2(ig,l+1) * (1. + betalpha)) zw2m = zw2(ig,l+1) / 2. entr_star(ig,l) = zdz * zbetalpha * (afact * RG * (ztv2(ig,l) & & - ztv2(ig,l+1)) / ztv2(ig,l+1) / zw2m - fact_epsilon) IF ((ztv2(ig,l).gt.ztv2(ig,l+1)).and.(entr_star(ig,l).GT.0.)) THEN active(ig) = .true. lmin(ig) = l f_star(ig,l+1) = entr_star(ig,l) zw2_est(ig,l+1) = zw2(ig,l+1) ELSE zw2(ig,l+1) = 0. entr_star(ig,l) = 0. ENDIF l = l + 1 ENDDO ENDDO !=============================================================================== ! Thermal plumes computations !=============================================================================== DO l=2,nlay-1 !------------------------------------------------------------------------------- ! Check if thermal plume is (still) active !------------------------------------------------------------------------------- !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! AB: we decide here if the plume is still active or not. When the plume's ! first level is reached, we set active to "true". Otherwise, it is given ! by zw2, f_star and entr_star. !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ DO ig=1,ngrid IF (l==lmin(ig)+1) THEN active(ig) = .true. ENDIF active(ig) = active(ig) & & .and. (zw2(ig,l).GT.1.e-10) & & .and. (f_star(ig,l) + entr_star(ig,l)).GT.1.e-10 ENDDO ztemp(:) = zpopsk(:,l) * zhla(:,l-1) DO ig=1,ngrid CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsa_est(ig)) ENDDO !------------------------------------------------------------------------------- ! Vertical speed (before mixing between plume and environment) !------------------------------------------------------------------------------- DO ig=1,ngrid IF (active(ig)) THEN zqla_est(ig,l) = max(0.,zqta(ig,l-1)-zqsa_est(ig)) ! zqla_est set to ql plume zta_est(ig,l) = zhla(ig,l-1) * zpopsk(ig,l) & ! zta_est set to TR plume & + RLvCp * zqla_est(ig,l) ztva_est(ig,l) = zta_est(ig,l) / zpopsk(ig,l) & ! ztva_est set to TRPV plume & * (1. + RETV * (zqta(ig,l-1)-zqla_est(ig,l)) - zqla_est(ig,l)) zbuoy(ig,l) = RG * (ztva_est(ig,l) - ztv(ig,l)) / ztv(ig,l) zdz = zlev(ig,l+1) - zlev(ig,l) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! AB: initial formulae !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha) ! zdw2 = afact * zbuoy(ig,l) / fact_epsilon ! zdw2bis = afact * zbuoy(ig,l-1) / fact_epsilon ! zw2_est(ig,l+1) = Max(0.0001,exp(-zw2fact)*(zw2_est(ig,l)-zdw2)+zdw2) ! zw2_est(ig,l+1) = Max(0.0001,exp(-zw2fact)*(zw2_est(ig,l)-zdw2bis)+zdw2) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! AB: own derivation for zw2_est (Rio et al. 2010) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! zw2fact = 2. * fact_epsilon * zdz ! zdw2 = 2. * afact * zbuoy(ig,l) * zdz zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha) zdw2 = 2. * afact * zbuoy(ig,l) * zdz / (1. + betalpha) zw2_est(ig,l+1) = Max(0., exp(-zw2fact) * zw2_est(ig,l) + zdw2) ENDIF ENDDO !------------------------------------------------------------------------------- ! Mass flux, entrainment and detrainment !------------------------------------------------------------------------------- DO ig=1,ngrid IF (active(ig)) THEN zdz = zlev(ig,l+1) - zlev(ig,l) zw2m = (zw2_est(ig,l+1) + zw2_est(ig,l)) / 2. ! AB: est-ce la bonne vitesse a utiliser ? gamma = afact * zbuoy(ig,l) - fact_epsilon * zw2m IF (zw2_est(ig,l).GT.0.) THEN test = gamma / zw2m - nu ELSE print *, 'ERROR: zw2_est is negative while plume is active!' print *, 'ig,l', ig, l print *, 'zw2_est', zw2_est(ig,l) call abort ENDIF IF (test.gt.0.) THEN detr_star(ig,l) = zdz * nu entr_star(ig,l) = zdz * (zbetalpha * gamma / zw2m + nu) ELSE detr_star(ig,l) = zdz * (nu - betalpha * gamma / zw2m) entr_star(ig,l) = zdz * nu ENDIF ! IF (detr_star(ig,l).lt.0.) THEN ! print *, 'WARNING: detrainment is negative!' ! print *, 'l,detr', l, detr_star(ig,l) ! ENDIF ! IF (entr_star(ig,l).lt.0.) THEN ! print *, 'WARNING: entrainment is negative!' ! print *, 'l,entr', l, entr_star(ig,l) ! ENDIF f_star(ig,l+1) = f_star(ig,l) + entr_star(ig,l) - detr_star(ig,l) ! IF (f_star(ig,l+1).le.0.) THEN ! print *, 'WARNING: mass flux is negative!' ! print *, 'l,f_star', l+1, f_star(ig,l+1) ! ENDIF ENDIF ENDDO !------------------------------------------------------------------------------- ! Mixing between thermal plume and environment !------------------------------------------------------------------------------- activetmp(:) = active(:) .and. (f_star(:,l+1).GT.1.e-10) DO ig=1,ngrid IF (activetmp(ig)) THEN zhla(ig,l) = (f_star(ig,l) * zhla(ig,l-1) & ! zhla is set to TP in plume (mixed) & + entr_star(ig,l) * zhl(ig,l)) & & / (f_star(ig,l+1) + detr_star(ig,l)) zqta(ig,l) = (f_star(ig,l) * zqta(ig,l-1) + & ! zqta is set to qt in plume (mixed) & + entr_star(ig,l) * zqt(ig,l)) & & / (f_star(ig,l+1) + detr_star(ig,l)) ENDIF ENDDO ztemp(:) = zpopsk(:,l) * zhla(:,l) DO ig=1,ngrid IF (activetmp(ig)) THEN CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsa(ig,l)) ENDIF ENDDO !------------------------------------------------------------------------------- ! Vertical speed (after mixing between plume and environment) !------------------------------------------------------------------------------- DO ig=1,ngrid IF (activetmp(ig)) THEN zqla(ig,l) = max(0.,zqta(ig,l)-zqsa(ig,l)) ! zqla is set to ql plume (mixed) zta(ig,l) = zhla(ig,l) * zpopsk(ig,l) + RLvCp * zqla(ig,l) ! ztva is set to TR plume (mixed) ztva(ig,l) = zta(ig,l) / zpopsk(ig,l) & ! ztva is set to TRPV plume (mixed) & * (1. + RETV*(zqta(ig,l)-zqla(ig,l)) - zqla(ig,l)) zbuoy(ig,l) = RG * (ztva(ig,l) - ztv(ig,l)) / ztv(ig,l) zdz = zlev(ig,l+1) - zlev(ig,l) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! AB: initial formula !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha) ! zdw2 = afact * zbuoy(ig,l) / fact_epsilon ! zw2(ig,l+1) = Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! AB: own derivation for zw2 (Rio et al. 2010) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! zw2fact = 2. * (fact_epsilon * zdz + entr_star(ig,l) / f_star(ig,l)) ! zdw2 = 2. * afact * zbuoy(ig,l) * zdz zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha) zdw2 = 2. * afact * zbuoy(ig,l) * zdz / (1. + betalpha) zw2(ig,l+1) = Max(0., exp(-zw2fact) * zw2(ig,l) + zdw2) ENDIF ENDDO ENDDO RETURN END