Changeset 2178
- Timestamp:
- Nov 12, 2019, 2:57:22 PM (5 years ago)
- Location:
- trunk/LMDZ.GENERIC
- Files:
-
- 2 edited
-
README (modified) (1 diff)
-
libf/phystd/thermcell_plume.F90 (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/README
r2177 r2178 1498 1498 e > e_max is no longer permitted ; b <= incoming mass flux is checked last 1499 1499 - Cleanup thermal plums model subroutines (thermcell_main, thermcell_env, thercell_dq, thermcell_dv2, thermcell_closure, thermcell_height) 1500 - In thermcell_plume, restore initial formula to compute the vertical speed and cleanup the file -
trunk/LMDZ.GENERIC/libf/phystd/thermcell_plume.F90
r2145 r2178 2 2 ! 3 3 ! 4 SUBROUTINE thermcell_plume(ngrid,nlay,nq,ptimestep, ztv,&5 z hl,zqt,zql,rhobarz,zlev,pplev,pphi,zpopsk,&4 SUBROUTINE thermcell_plume(ngrid,nlay,nq,ptimestep, & 5 ztv,zhl,zqt,zql,zlev,pplev,zpopsk, & 6 6 detr_star,entr_star,f_star, & 7 ztva,zhla,zq la,zqta,zta,zqsa,&8 zw2,lmi x,lmin)7 ztva,zhla,zqta,zqla,zqsa, & 8 zw2,lmin) 9 9 10 10 … … 37 37 ! ------- 38 38 39 INTEGER ngrid, nlay, nq40 41 REAL ptimestep42 REAL rhobarz(ngrid,nlay) ! Levels density43 REAL zlev(ngrid,nlay+1) ! Levels altitude44 REAL pplev(ngrid,nlay+1) ! Levels pressure45 REAL pphi(ngrid,nlay) ! Geopotential46 REAL zpopsk(ngrid,nlay)! Exner function47 48 REAL ztv(ngrid,nlay)! TRPV environment49 REAL zhl(ngrid,nlay)! TP environment50 REAL zqt(ngrid,nlay)! qt environment51 REAL zql(ngrid,nlay)! ql environment39 INTEGER, INTENT(in) :: ngrid 40 INTEGER, INTENT(in) :: nlay 41 INTEGER, INTENT(in) :: nq 42 43 REAL, INTENT(in) :: ptimestep 44 REAL, INTENT(in) :: zlev(ngrid,nlay+1) ! Levels altitude 45 REAL, INTENT(in) :: pplev(ngrid,nlay+1) ! Levels pressure 46 REAL, INTENT(in) :: zpopsk(ngrid,nlay) ! Exner function 47 48 REAL, INTENT(in) :: ztv(ngrid,nlay) ! TRPV environment 49 REAL, INTENT(in) :: zhl(ngrid,nlay) ! TP environment 50 REAL, INTENT(in) :: zqt(ngrid,nlay) ! qt environment 51 REAL, INTENT(in) :: zql(ngrid,nlay) ! ql environment 52 52 53 53 ! Outputs: 54 54 ! -------- 55 55 56 INTEGER lmin(ngrid) ! plume base level (first unstable level) 57 INTEGER lmix(ngrid) ! maximum vertical speed level 58 59 REAL detr_star(ngrid,nlay) ! normalized detrainment 60 REAL entr_star(ngrid,nlay) ! normalized entrainment 61 REAL f_star(ngrid,nlay+1) ! normalized mass flux 62 63 REAL ztva(ngrid,nlay) ! TRPV plume (after mixing) 64 REAL zhla(ngrid,nlay) ! TP plume ? 65 REAL zqla(ngrid,nlay) ! ql plume (after mixing) 66 REAL zqta(ngrid,nlay) ! qt plume ? 67 REAL zqsa(ngrid,nlay) ! qsat plume (after mixing) 68 REAL zw2(ngrid,nlay+1) ! w plume (after mixing) 56 INTEGER, INTENT(out) :: lmin(ngrid) ! Plume bottom level (first unstable level) 57 58 REAL, INTENT(out) :: detr_star(ngrid,nlay) ! Normalized detrainment 59 REAL, INTENT(out) :: entr_star(ngrid,nlay) ! Normalized entrainment 60 REAL, INTENT(out) :: f_star(ngrid,nlay+1) ! Normalized mass flux 61 62 REAL, INTENT(out) :: ztva(ngrid,nlay) ! TRPV plume (after mixing) 63 REAL, INTENT(out) :: zhla(ngrid,nlay) ! TP plume (after mixing) 64 REAL, INTENT(out) :: zqla(ngrid,nlay) ! ql plume (after mixing) 65 REAL, INTENT(out) :: zqta(ngrid,nlay) ! qt plume (after mixing) 66 REAL, INTENT(out) :: zqsa(ngrid,nlay) ! qsat plume (after mixing) 67 REAL, INTENT(out) :: zw2(ngrid,nlay+1) ! w plumr (after mixing) 69 68 70 69 ! Local: … … 73 72 INTEGER ig, l, k 74 73 75 REAL ztva_est(ngrid,nlay) ! TRPV plume (before mixing)76 REAL zqla_est(ngrid,nlay) ! ql plume (before mixing)77 REAL zta_est(ngrid,nlay) ! TR plume (before mixing)78 REAL zqsa_est(ngrid) ! qsat plume (before mixing)79 REAL zw2_est(ngrid,nlay+1) ! w plume (before mixing)80 81 REAL zta(ngrid,nlay) ! TR plume (after mixing)82 83 REAL zbuoy(ngrid,nlay) ! Plume buoyancy84 REAL ztemp(ngrid) ! Temperature for saturation vapor pressure computation in plume85 REAL zdz ! Layers heights86 REAL ztv2(ngrid,nlay) ! ztv + d_temp * Dirac(l=linf)87 88 REAL z betalpha !89 REAL z dw2!90 REAL z dw2bis !91 REAL zw2fact !92 REAL zw2m ! Average vertical velocity between two successive levels93 REAL gamma ! Plume acceleration term (to compute vertical velocity)94 REAL test ! Test to know how to compute entrainment and detrainment95 96 REAL psat ! Dummy argument for Psat_water() 97 98 LOGICAL active(ngrid) ! If the plume is active at ig (speed and incoming mass flux > 0 or l=lmin)99 LOGICAL activetmp(ngrid) ! If the plume is active at ig (active=true and outgoing mass flux > 0)74 REAL ztva_est(ngrid,nlay) ! TRPV plume (before mixing) 75 REAL zqla_est(ngrid,nlay) ! ql plume (before mixing) 76 REAL zta_est(ngrid,nlay) ! TR plume (before mixing) 77 REAL zqsa_est(ngrid) ! qsat plume (before mixing) 78 REAL zw2_est(ngrid,nlay+1) ! w plume (before mixing) 79 80 REAL zta(ngrid,nlay) ! TR plume (after mixing) 81 82 REAL zbuoy(ngrid,nlay) ! Plume buoyancy 83 REAL ztemp(ngrid) ! Temperature for saturation vapor pressure computation in plume 84 REAL zdz ! Layers heights 85 REAL ztv2(ngrid,nlay) ! ztv + d_temp * Dirac(l=linf) 86 87 REAL zdw2 ! 88 REAL zw2fact ! 89 REAL zw2m ! Average vertical velocity between two successive levels 90 REAL gamma ! Plume acceleration term (to compute vertical velocity) 91 REAL test ! 92 93 REAL psat ! Dummy argument for Psat_water() 94 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) 100 99 101 100 !=============================================================================== … … 103 102 !=============================================================================== 104 103 105 zbetalpha = betalpha / (1. + betalpha) 106 107 ztva(:,:) = ztv(:,:) ! ztva is set to TPV environment 108 zhla(:,:) = zhl(:,:) ! zhla is set to TP environment 109 zqta(:,:) = zqt(:,:) ! zqta is set to qt environment 110 zqla(:,:) = zql(:,:) ! zqla is set to ql environment 111 112 zqsa_est(:) = 0. 113 zqsa(:,:) = 0. 114 115 zw2_est(:,:) = 0. 116 zw2(:,:) = 0. 117 118 zbuoy(:,:) = 0. 119 120 f_star(:,:) = 0. 121 detr_star(:,:) = 0. 122 entr_star(:,:) = 0. 123 124 lmix(:) = 1 125 lmin(:) = 1 126 127 ztv2(:,:) = ztv(:,:) 128 ztv2(:,linf) = ztv(:,linf) + d_temp 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. 114 115 zbuoy(:,:) = 0. 116 117 f_star(:,:) = 0. 118 detr_star(:,:) = 0. 119 entr_star(:,:) = 0. 120 121 lmin(:) = 1 122 123 ztv2(:,:) = ztv(:,:) 124 ztv2(:,linf) = ztv(:,linf) + d_temp 129 125 130 126 active(:) = .false. … … 138 134 DO WHILE (.not.active(ig).and.(pplev(ig,l+1) > pres_limit).and.(l < nlay)) 139 135 zbuoy(ig,l) = RG * (ztv2(ig,l) - ztv2(ig,l+1)) / ztv2(ig,l+1) 140 zdz = zlev(ig,l+1) - zlev(ig,l) 141 zw2m = afact * zbuoy(ig,l) * zdz / (1. + betalpha) 142 ! gamma = afact * zbuoy(ig,l) - fact_epsilon * zw2m 143 ! test = gamma / zw2m - nu 144 test = zbuoy(ig,l) 145 IF (test > 0.) THEN 136 IF (zbuoy(ig,l) > 0.) THEN 146 137 lmin(ig) = l 147 ! entr_star(ig,l) = zdz * f_star(ig,l) * zbetalpha * gamma / zw2m - nu ! Problem because f*(ig,l) = 0 148 ! detr_star(ig,l) = f_star(ig,l) * nu ! Problem because f*(ig,l) = 0 149 ! f_star(ig,l+1) = entr_star(ig,l) - detr_star(ig,l) 138 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 139 ! AB: entrainement and mass flux initial values are set to 1. The physical value 140 ! will be computed thanks to the closure relation in thermcell_closure. 141 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 150 142 entr_star(ig,l) = 1. 151 143 f_star(ig,l+1) = 1. 152 zw2_est(ig,l+1) = zw2m * 2. 144 zdz = zlev(ig,l+1) - zlev(ig,l) 145 zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha) 146 zdw2 = 2. * afact * zbuoy(ig,l) * zdz / (1. + betalpha) 147 zw2_est(ig,l+1) = exp(-zw2fact) * zdw2 153 148 zw2(ig,l+1) = zw2_est(ig,l+1) 154 149 active(ig) = .true. … … 162 157 !=============================================================================== 163 158 164 DO l= 2,nlay-1159 DO l=linf+1,nlay-1 165 160 166 161 !------------------------------------------------------------------------------- … … 169 164 170 165 DO ig=1,ngrid 171 active(ig) = (active(ig).or.(l == lmin(ig)+1)) & 172 & .and.(zw2(ig,l) > 1.e-10) & 173 & .and.(f_star(ig,l) > 1.e-10) 166 active(ig) = (zw2(ig,l) > 1.e-10).and.(f_star(ig,l) > 1.e-10) 174 167 ENDDO 175 168 … … 181 174 182 175 DO ig=1,ngrid 183 CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsa_est(ig)) 176 IF (active(ig)) THEN 177 CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsa_est(ig)) 178 ENDIF 184 179 ENDDO 185 180 … … 199 194 zdz = zlev(ig,l+1) - zlev(ig,l) 200 195 201 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~202 ! AB: initial formulae203 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~204 ! zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)205 ! zdw2 = afact * zbuoy(ig,l) / fact_epsilon206 ! zdw2bis = afact * zbuoy(ig,l-1) / fact_epsilon207 ! zw2_est(ig,l+1) = Max(0.0001,exp(-zw2fact)*(zw2_est(ig,l)-zdw2)+zdw2)208 ! zw2_est(ig,l+1) = Max(0.0001,exp(-zw2fact)*(zw2_est(ig,l)-zdw2bis)+zdw2)209 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~210 ! AB: own derivation for zw2_est (Rio et al. 2010)211 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~212 ! zw2fact = 2. * fact_epsilon * zdz213 ! zdw2 = 2. * afact * zbuoy(ig,l) * zdz214 196 zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha) 215 zdw2 = 2. * afact * zbuoy(ig,l) * zdz / (1. + betalpha)216 zw2_est(ig,l+1) = Max(0., exp(-zw2fact) * zw2_est(ig,l) + zdw2)197 zdw2 = afact * zbuoy(ig,l) / fact_epsilon 198 zw2_est(ig,l+1) = Max(0., exp(-zw2fact) * (zw2_est(ig,l) - zdw2) + zdw2) 217 199 ENDIF 218 200 ENDDO … … 240 222 IF (test > 0.) THEN 241 223 detr_star(ig,l) = zdz * f_star(ig,l) * nu 242 entr_star(ig,l) = zdz * f_star(ig,l) * ( zbetalpha * gamma / zw2m + nu)224 entr_star(ig,l) = zdz * f_star(ig,l) * (betalpha * gamma / zw2m + nu) / (betalpha + 1) 243 225 ELSE 244 detr_star(ig,l) = zdz * f_star(ig,l) * ( nu - betalpha * gamma / zw2m)226 detr_star(ig,l) = zdz * f_star(ig,l) * ((betalpha + 1) * nu - betalpha * gamma / zw2m) 245 227 entr_star(ig,l) = zdz * f_star(ig,l) * nu 246 228 ENDIF … … 295 277 zdz = zlev(ig,l+1) - zlev(ig,l) 296 278 297 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~298 ! AB: initial formula299 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~300 ! zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)301 ! zdw2 = afact * zbuoy(ig,l) / fact_epsilon302 ! zw2(ig,l+1) = Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)303 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~304 ! AB: own derivation for zw2 (Rio et al. 2010)305 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~306 ! zw2fact = 2. * (fact_epsilon * zdz + entr_star(ig,l) / f_star(ig,l))307 ! zdw2 = 2. * afact * zbuoy(ig,l) * zdz308 279 zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha) 309 zdw2 = 2. * afact * zbuoy(ig,l) * zdz / (1. + betalpha)310 zw2(ig,l+1) = Max(0., exp(-zw2fact) * zw2(ig,l) + zdw2)280 zdw2 = afact * zbuoy(ig,l) / fact_epsilon 281 zw2(ig,l+1) = Max(0., exp(-zw2fact) * (zw2(ig,l) - zdw2) + zdw2) 311 282 ENDIF 312 283 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
