Changeset 2177 for trunk/LMDZ.GENERIC
- Timestamp:
- Nov 12, 2019, 11:55:58 AM (5 years ago)
- Location:
- trunk/LMDZ.GENERIC
- Files:
-
- 8 edited
-
README (modified) (1 diff)
-
libf/phystd/thermcell_closure.F90 (modified) (2 diffs)
-
libf/phystd/thermcell_dq.F90 (modified) (4 diffs)
-
libf/phystd/thermcell_dv2.F90 (modified) (5 diffs)
-
libf/phystd/thermcell_env.F90 (modified) (2 diffs)
-
libf/phystd/thermcell_flux.F90 (modified) (12 diffs)
-
libf/phystd/thermcell_height.F90 (modified) (10 diffs)
-
libf/phystd/thermcell_main.F90 (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/README
r2176 r2177 1495 1495 == 12/11/2019 == AB 1496 1496 - fm0, entr0 and detr0 are now allocatable variables in physiq_mod. That is necessary if tau_thermals > 0. 1497 1497 - In thermcell_flux, "bidouilles" are modified: now the plumes stop when the updraft fraction is greater than alpha_max ; 1498 e > e_max is no longer permitted ; b <= incoming mass flux is checked last 1499 - Cleanup thermal plums model subroutines (thermcell_main, thermcell_env, thercell_dq, thermcell_dv2, thermcell_closure, thermcell_height) -
trunk/LMDZ.GENERIC/libf/phystd/thermcell_closure.F90
r2145 r2177 29 29 ! ------- 30 30 31 INTEGER ngrid, nlay 32 INTEGER lmax(ngrid) 31 INTEGER, INTENT(in) :: ngrid 32 INTEGER, INTENT(in) :: nlay 33 INTEGER, INTENT(in) :: lmax(ngrid) 33 34 34 REAL ptimestep35 REAL rho(ngrid,nlay)36 REAL zlev(ngrid,nlay)37 REAL alim_star(ngrid,nlay)38 REAL zmin(ngrid)39 REAL zmax(ngrid)40 REAL wmax(ngrid)35 REAL, INTENT(in) :: ptimestep 36 REAL, INTENT(in) :: rho(ngrid,nlay) 37 REAL, INTENT(in) :: zlev(ngrid,nlay) 38 REAL, INTENT(in) :: alim_star(ngrid,nlay) 39 REAL, INTENT(in) :: zmin(ngrid) 40 REAL, INTENT(in) :: zmax(ngrid) 41 REAL, INTENT(in) :: wmax(ngrid) 41 42 42 43 ! Outputs: 43 44 ! -------- 44 45 45 REAL f(ngrid)46 REAL, INTENT(out) :: f(ngrid) 46 47 47 48 ! Local: … … 90 91 DO l=1,llmax-1 91 92 DO ig=1,ngrid 92 IF (l < lmax(ig)) THEN93 93 alim_star2(ig) = alim_star2(ig) + alim_star(ig,l)**2 & 94 & / (rho(ig,l) * (zlev(ig,l+1) - zlev(ig,l))) ! => inte rgration is ok because alim_star = a* dz94 & / (rho(ig,l) * (zlev(ig,l+1) - zlev(ig,l))) ! => integration is ok because alim_star = a* dz 95 95 alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig,l) 96 ENDIF97 96 ENDDO 98 97 ENDDO -
trunk/LMDZ.GENERIC/libf/phystd/thermcell_dq.F90
r2144 r2177 3 3 ! 4 4 SUBROUTINE thermcell_dq(ngrid,nlay,ptimestep,fm,entr,detr,masse, & 5 q,dq,qa,lmin )5 q,dq,qa,lmin,lmax) 6 6 7 7 … … 16 16 ! 17 17 ! Modif 2019/04 (AB alexandre.boissinot@lmd.jussieu.fr) 18 ! dqimpl = 1: implicit scheme19 ! dqimpl = 0: explicit scheme18 ! dqimpl = true : implicit scheme 19 ! dqimpl = false : explicit scheme 20 20 ! 21 21 !=============================================================================== … … 34 34 ! ------- 35 35 36 INTEGER ngrid, nlay 37 INTEGER lmin(ngrid) 36 INTEGER, INTENT(in) :: ngrid 37 INTEGER, INTENT(in) :: nlay 38 INTEGER, INTENT(in) :: lmin(ngrid) 39 INTEGER, INTENT(in) :: lmax(ngrid) 38 40 39 REAL ptimestep 40 REAL masse(ngrid,nlay) 41 REAL fm(ngrid,nlay+1) 42 REAL entr(ngrid,nlay) 43 REAL detr(ngrid,nlay) 44 REAL q(ngrid,nlay) 41 REAL, INTENT(in) :: ptimestep 42 REAL, INTENT(in) :: masse(ngrid,nlay) 43 REAL, INTENT(in) :: fm(ngrid,nlay+1) 44 REAL, INTENT(in) :: entr(ngrid,nlay) 45 REAL, INTENT(in) :: detr(ngrid,nlay) 45 46 46 47 ! Outputs: 47 48 ! -------- 48 49 49 REAL dq(ngrid,nlay) 50 REAL qa(ngrid,nlay) 50 REAL, INTENT(inout) :: q(ngrid,nlay) 51 REAL, INTENT(out) :: dq(ngrid,nlay) 52 REAL, INTENT(out) :: qa(ngrid,nlay) 51 53 52 54 ! Local: … … 82 84 cfl = max(cfl, fm(ig,l) / zzm) 83 85 84 IF (entr(ig,l) .gt.zzm) THEN86 IF (entr(ig,l) > zzm) THEN 85 87 print *, 'ERROR: entrainment is greater than the layer mass!' 86 88 print *, 'ig,l,entr', ig, l, entr(ig,l) 87 print *, 'lmin ', lmin(ig)89 print *, 'lmin,lmax', lmin(ig), lmax(ig) 88 90 print *, '-------------------------------' 89 91 print *, 'entr*dt,mass', entr(ig,l)*ptimestep, masse(ig,l) -
trunk/LMDZ.GENERIC/libf/phystd/thermcell_dv2.F90
r2143 r2177 16 16 !=============================================================================== 17 17 18 USE print_control_mod, ONLY: prt_level, lunout18 USE print_control_mod, ONLY: prt_level, lunout 19 19 20 20 IMPLICIT none … … 28 28 ! ------- 29 29 30 INTEGER ngrid, nlay 30 INTEGER, INTENT(in) :: ngrid 31 INTEGER, INTENT(in) :: nlay 31 32 32 REAL ptimestep33 REAL masse(ngrid,nlay)34 REAL fm(ngrid,nlay+1)35 REAL entr(ngrid,nlay)36 REAL detr(ngrid,nlay)37 REAL fraca(ngrid,nlay+1)38 REAL zmax(ngrid)39 REAL zmin(ngrid)40 REAL u(ngrid,nlay)41 REAL v(ngrid,nlay)33 REAL, INTENT(in) :: ptimestep 34 REAL, INTENT(in) :: masse(ngrid,nlay) 35 REAL, INTENT(in) :: fm(ngrid,nlay+1) 36 REAL, INTENT(in) :: entr(ngrid,nlay) 37 REAL, INTENT(in) :: detr(ngrid,nlay) 38 REAL, INTENT(in) :: fraca(ngrid,nlay+1) 39 REAL, INTENT(in) :: zmax(ngrid) 40 REAL, INTENT(in) :: zmin(ngrid) 41 REAL, INTENT(in) :: u(ngrid,nlay) 42 REAL, INTENT(in) :: v(ngrid,nlay) 42 43 43 44 ! Outputs: 44 45 ! -------- 45 46 46 REAL dua(ngrid,nlay) 47 REAL dva(ngrid,nlay) 48 REAL ua(ngrid,nlay) 49 REAL va(ngrid,nlay) 50 REAL du(ngrid,nlay) 51 REAL dv(ngrid,nlay) 47 REAL, INTENT(out) :: ua(ngrid,nlay) 48 REAL, INTENT(out) :: va(ngrid,nlay) 49 REAL, INTENT(out) :: du(ngrid,nlay) 50 REAL, INTENT(out) :: dv(ngrid,nlay) 52 51 53 52 ! Local: … … 66 65 REAL ue(ngrid,nlay) 67 66 REAL ve(ngrid,nlay) 67 REAL dua(ngrid,nlay) 68 REAL dva(ngrid,nlay) 68 69 REAL plume_height(ngrid) 69 70 … … 121 122 DO iter=1,5 122 123 DO ig=1,ngrid 123 ! Calcul prenant en compte la fraction reelle124 ! FH: Calcul prenant en compte la fraction reelle 124 125 zf = 0.5 * (fraca(ig,l) + fraca(ig,l+1)) 125 126 zf2 = 1./(1.-zf) 126 ! Calcul avec fraction infiniement petite127 ! FH: Calcul avec fraction infiniement petite 127 128 ! zf = 0. 128 129 ! zf2 = 1. … … 133 134 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 134 135 IF (ltherm(ig,l)) THEN 135 ! 136 ! FH: On choisit une relaxation lineaire. 136 137 ! gamma(ig,l) = gamma0(ig,l) 137 ! 138 ! FH: On choisit une relaxation quadratique. 138 139 gamma(ig,l) = gamma0(ig,l) * sqrt(dua(ig,l)**2 + dva(ig,l)**2) 139 140 -
trunk/LMDZ.GENERIC/libf/phystd/thermcell_env.F90
r2143 r2177 38 38 ! ------- 39 39 40 INTEGER ngrid, nlay, nq 40 INTEGER, INTENT(in) :: ngrid 41 INTEGER, INTENT(in) :: nlay 42 INTEGER, INTENT(in) :: nq 41 43 42 REAL pq(ngrid,nlay,nq)! Large scale water43 REAL pt(ngrid,nlay)! Large scale temperature44 REAL pu(ngrid,nlay)! Large scale zonal wind45 REAL pv(ngrid,nlay)! Large scale meridional wind46 REAL pplay(ngrid,nlay)! Layers pressure47 REAL pplev(ngrid,nlay+1)! Levels pressure48 REAL zpopsk(ngrid,nlay)! Exner function44 REAL, INTENT(in) :: pq(ngrid,nlay,nq) ! Large scale water 45 REAL, INTENT(in) :: pt(ngrid,nlay) ! Large scale temperature 46 REAL, INTENT(in) :: pu(ngrid,nlay) ! Large scale zonal wind 47 REAL, INTENT(in) :: pv(ngrid,nlay) ! Large scale meridional wind 48 REAL, INTENT(in) :: pplay(ngrid,nlay) ! Layers pressure 49 REAL, INTENT(in) :: pplev(ngrid,nlay+1) ! Levels pressure 50 REAL, INTENT(in) :: zpopsk(ngrid,nlay) ! Exner function 49 51 50 52 ! Outputs: 51 53 ! -------- 52 54 53 REAL zqt(ngrid,nlay) ! qtenvironment54 REAL zql(ngrid,nlay) ! qlenvironment55 REAL zt(ngrid,nlay) ! Tenvironment56 REAL ztv(ngrid,nlay) ! TRPVenvironment57 REAL zhl(ngrid,nlay) ! TPenvironment58 REAL zu(ngrid,nlay) ! uenvironment59 REAL zv(ngrid,nlay) ! venvironment60 REAL zqs(ngrid,nlay)! qsat environment55 REAL, INTENT(out) :: zt(ngrid,nlay) ! T environment 56 REAL, INTENT(out) :: ztv(ngrid,nlay) ! TRPV environment 57 REAL, INTENT(out) :: zhl(ngrid,nlay) ! TP environment 58 REAL, INTENT(out) :: zu(ngrid,nlay) ! u environment 59 REAL, INTENT(out) :: zv(ngrid,nlay) ! v environment 60 REAL, INTENT(out) :: zqt(ngrid,nlay) ! qt environment 61 REAL, INTENT(out) :: zql(ngrid,nlay) ! ql environment 62 REAL, INTENT(out) :: zqs(ngrid,nlay) ! qsat environment 61 63 62 64 ! Local: 63 65 ! ------ 64 66 65 INTEGER ig, l67 INTEGER ig, k 66 68 67 REAL psat ! Dummy argument for Psat_water()69 REAL psat ! Dummy argument for Psat_water() 68 70 69 71 !=============================================================================== … … 85 87 IF (water) THEN 86 88 87 DO l=1,nlay89 DO k=1,nlay 88 90 DO ig=1,ngrid 89 CALL Psat_water(pt(ig, l), pplev(ig,l), psat, zqs(ig,l))91 CALL Psat_water(pt(ig,k), pplev(ig,k), psat, zqs(ig,k)) 90 92 ENDDO 91 93 ENDDO 92 94 93 DO l=1,nlay95 DO k=1,nlay 94 96 DO ig=1,ngrid 95 zql(ig, l) = max(0.,pq(ig,l,igcm_h2o_vap) - zqs(ig,l))96 zt(ig, l) = pt(ig,l) + RLvCp * zql(ig,l)97 ztv(ig, l) = zt(ig,l) / zpopsk(ig,l) &98 & * (1. + RETV * (zqt(ig, l)-zql(ig,l)) - zql(ig,l))97 zql(ig,k) = max(0.,pq(ig,k,igcm_h2o_vap) - zqs(ig,k)) 98 zt(ig,k) = pt(ig,k) + RLvCp * zql(ig,k) 99 ztv(ig,k) = zt(ig,k) / zpopsk(ig,k) & 100 & * (1. + RETV * (zqt(ig,k)-zql(ig,k)) - zql(ig,k)) 99 101 ENDDO 100 102 ENDDO -
trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90
r2146 r2177 27 27 ! ------- 28 28 29 INTEGER ngrid, nlay30 INTEGER lmin(ngrid)31 INTEGER lmax(ngrid)32 33 REAL entr_star(ngrid,nlay)34 REAL detr_star(ngrid,nlay)35 REAL zw2(ngrid,nlay+1)36 REAL zlev(ngrid,nlay+1)37 REAL masse(ngrid,nlay)38 REAL ptimestep39 REAL rhobarz(ngrid,nlay)40 REAL f(ngrid)29 INTEGER, INTENT(in) :: ngrid 30 INTEGER, INTENT(in) :: nlay 31 INTEGER, INTENT(in) :: lmin(ngrid) 32 33 REAL, INTENT(in) :: entr_star(ngrid,nlay) 34 REAL, INTENT(in) :: detr_star(ngrid,nlay) 35 REAL, INTENT(in) :: zw2(ngrid,nlay+1) 36 REAL, INTENT(in) :: zlev(ngrid,nlay+1) 37 REAL, INTENT(in) :: masse(ngrid,nlay) 38 REAL, INTENT(in) :: ptimestep 39 REAL, INTENT(in) :: rhobarz(ngrid,nlay) 40 REAL, INTENT(in) :: f(ngrid) 41 41 42 42 ! Outputs: 43 ! -------- 44 45 REAL entr(ngrid,nlay) 46 REAL detr(ngrid,nlay) 47 REAL fm(ngrid,nlay+1) 43 ! -------- 44 45 INTEGER, INTENT(inout) :: lmax(ngrid) 46 47 REAL, INTENT(out) :: entr(ngrid,nlay) 48 REAL, INTENT(out) :: detr(ngrid,nlay) 49 REAL, INTENT(out) :: fm(ngrid,nlay+1) 48 50 49 51 ! Local: … … 53 55 INTEGER igout, lout ! Error grid point and level 54 56 55 REAL zfm ! Mass flux such as updraft fraction is equal to its maximal authorized value alpha_max 56 REAL f_old ! Save initial value of mass flux 57 REAL fmax ! Maximal authorized mass flux (alpha < alpha_max) 58 REAL fff0 ! Save initial value of mass flux 59 REAL emax ! Maximal authorized entrainment (entr*dt < mass_max) 57 60 REAL eee0 ! Save initial value of entrainment 58 61 REAL ddd0 ! Save initial value of detrainment 59 62 REAL eee ! eee0 - layer mass * maximal authorized mass fraction / dt 60 63 REAL ddd ! ddd0 - eee 61 REAL zzz ! Temporary variable set to mass flux62 64 REAL fact 63 65 REAL test … … 103 105 104 106 !------------------------------------------------------------------------------- 105 ! Calcul du flux de masse107 ! Mass flux 106 108 !------------------------------------------------------------------------------- 107 109 … … 112 114 ELSEIF (l == lmax(ig)) THEN 113 115 fm(ig,l+1) = 0. 116 entr(ig,l) = 0. 114 117 detr(ig,l) = fm(ig,l) + entr(ig,l) 115 118 ELSE … … 170 173 print *, '---------------------------------------------------------' 171 174 print *, 'ERROR: mass flux has negative value(s)!' 172 print *, 'ig,l, entr', igout, lout, entr(igout,lout)175 print *, 'ig,l,norm', igout, lout, f(igout) 173 176 print *, 'lmin,lmax', lmin(igout), lmax(igout) 174 177 print *, '- - - - - - - - - - - - - - - - - - - - - - - - - - - - -' 175 178 DO k=nlay,1,-1 176 print *, ' fm ', fm(igout,k+1)179 print *, 'k, fm ', k+1, fm(igout,k+1) 177 180 print *, 'entr,detr', entr(igout,k), detr(igout,k) 178 181 ENDDO 179 print *, ' fm ', fm(igout,1)182 print *, 'k, fm ', 1, fm(igout,1) 180 183 print *, '---------------------------------------------------------' 181 184 CALL abort 182 185 ENDIF 183 184 !-------------------------------------------------------------------------------185 ! Is detrainment lessser than incoming mass flux ?186 !-------------------------------------------------------------------------------187 188 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~189 ! AB : Even if fm has no negative value, it can be lesser than detr.190 ! That's not suitable because when we will mix the plume with the191 ! environment, it will detrain more mass than it is physically able to do.192 ! When it occures, that imply that entr + fm is greater than detr,193 ! otherwise we get a negative outgoing mass flux (cf. above).194 ! That is why we correct the detrainment as follows.195 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~196 197 DO ig=1,ngrid198 IF (detr(ig,l).gt.fm(ig,l)) THEN199 detr(ig,l) = fm(ig,l)200 entr(ig,l) = fm(ig,l+1)201 ncorecdetr = ncorecdetr + 1202 ENDIF203 ENDDO204 186 205 187 !------------------------------------------------------------------------------- … … 213 195 ! If it's not enough, we can increase entr in the layer above and decrease 214 196 ! the outgoing mass flux in the current layer. 215 ! If it's still unsufficient, code will abort (now commented).197 ! If it's still insufficient, code will abort (now commented). 216 198 ! Else we reset entr to its intial value and we renormalize entrainment, 217 199 ! detrainment and mass flux profiles such as we do not exceed the maximal … … 222 204 eee0 = entr(ig,l) 223 205 ddd0 = detr(ig,l) 224 eee = entr(ig,l) - masse(ig,l) * fomass_max / ptimestep 225 ddd = detr(ig,l) - eee 226 IF (eee > 0.) THEN 227 entr(ig,l) = entr(ig,l) - eee 206 emax = masse(ig,l) * fomass_max / ptimestep 207 IF (emax < 0.) THEN 208 print *, 'ERROR: layer mass is negative!' 209 print *, 'mass,emax', masse(ig,l), emax 210 print *, 'ig,l', ig, l 211 ENDIF 212 IF (eee0 > emax) THEN 213 entr(ig,l) = emax 214 ddd = ddd0 + emax - eee0 228 215 ncorecentr = ncorecentr + 1 229 216 IF (ddd > 0.) THEN … … 235 222 fm(ig,l+1) = fm(ig,l) + entr(ig,l) 236 223 entr(ig,l+1) = entr(ig,l+1) + ddd 224 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 225 ! AB: Simulation abort (try to reduce the physical time step). 226 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 227 ELSE 228 entr(ig,l) = entr(ig,l) + eee 229 igout = ig 230 lout = l 231 labort_physic = .true. 232 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 233 ! AB: We can renormalize the plume mass fluxes. I think it does not work. 234 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 237 235 ! ELSE 238 ! entr(ig,l) = entr(ig,l) + eee 239 ! igout = ig 240 ! lout = l 241 ! labort_physic = .true. 242 ELSE 243 fact = max(fact, eee0 / (eee0 - eee)) 244 entr(ig,l) = eee0 245 ncorecfact = ncorecfact + 1 236 ! fact = max(fact, eee0 / emax) 237 ! entr(ig,l) = eee0 238 ! ncorecfact = ncorecfact + 1 239 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 240 ! AB: The renormalization can be just applied in the local plume. 241 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 242 ! DO k=lmin(ig),lmax(ig) 243 ! entr(ig,k) = entr(ig,k) * emax / eee0 244 ! detr(ig,k) = detr(ig,k) * emax / eee0 245 ! fm(ig,k) = fm(ig,k) * emax / eee0 246 ! ENDDO 246 247 ENDIF 247 248 ENDIF … … 272 273 273 274 !------------------------------------------------------------------------------- 274 ! Is Updraft fraction lesser than alpha_max ? 275 !------------------------------------------------------------------------------- 276 277 DO ig=1,ngrid 278 zfm = rhobarz(ig,l+1) * zw2(ig,l+1) * alpha_max 279 280 IF (fm(ig,l+1) > zfm) THEN 281 f_old = fm(ig,l+1) 282 fm(ig,l+1) = zfm 283 detr(ig,l) = detr(ig,l) + f_old - fm(ig,l+1) 275 ! Is updraft fraction lesser than alpha_max ? 276 !------------------------------------------------------------------------------- 277 278 DO ig=1,ngrid 279 fff0 = fm(ig,l+1) 280 fmax = rhobarz(ig,l+1) * zw2(ig,l+1) * alpha_max 281 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 282 ! AB: The plume mass flux can be reduced. 283 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 284 ! IF (fff0 > fmax) THEN 285 ! fm(ig,l+1) = fmax 286 ! detr(ig,l) = detr(ig,l) + fff0 - fmax 287 ! ncorecalpha = ncorecalpha + 1 288 ! entr(ig,l+1) = entr(ig,l+1) + fff0 - fmax 289 ! ENDIF 290 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 291 ! AB: The plume can be stopped here. 292 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 293 IF (fff0 > fmax) THEN 284 294 ncorecalpha = ncorecalpha + 1 285 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 286 ! AB : The previous change doesn't observe the equation df/dz = e - d. 287 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 288 entr(ig,l+1) = entr(ig,l+1) + f_old - fm(ig,l+1) 295 DO k=l+1,lmax(ig) 296 entr(ig,k) = 0. 297 detr(ig,k) = 0. 298 fm(ig,k) = 0. 299 ENDDO 300 lmax(ig) = l 301 entr(ig,l) = 0. 302 detr(ig,l) = fm(ig,l) 303 ENDIF 304 ENDDO 305 306 !------------------------------------------------------------------------------- 307 ! Is detrainment lesser than incoming mass flux ? 308 !------------------------------------------------------------------------------- 309 310 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 311 ! AB : Even if fm has no negative value, it can be lesser than detr. 312 ! That is not suitable because when we will mix the plume with the 313 ! environment, it will detrain more mass than it is physically able to do. 314 ! When it occures, that imply that entr + fm is greater than detr, 315 ! otherwise we get a negative outgoing mass flux (cf. above). 316 ! That is why we decrease entrainment and detrainment as follows. 317 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 318 319 DO ig=1,ngrid 320 IF (detr(ig,l) > fm(ig,l)) THEN 321 detr(ig,l) = fm(ig,l) 322 entr(ig,l) = fm(ig,l+1) 323 ncorecdetr = ncorecdetr + 1 289 324 ENDIF 290 325 ENDDO … … 292 327 ENDDO 293 328 294 IF (fact > 0.) THEN 295 entr(:,:) = entr(:,:) / fact 296 detr(:,:) = detr(:,:) / fact 297 fm(:,:) = fm(:,:) / fact 298 ENDIF 329 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 330 ! AB: The renormalization can be applied everywhere. 331 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 332 ! IF (fact > 0.) THEN 333 ! entr(:,:) = entr(:,:) / fact 334 ! detr(:,:) = detr(:,:) / fact 335 ! fm(:,:) = fm(:,:) / fact 336 ! ENDIF 299 337 300 338 !------------------------------------------------------------------------------- … … 316 354 317 355 DO ig=1,ngrid 318 IF (lmax(ig) .GT.0) THEN319 detr(ig,lmax(ig)) = fm(ig,lmax(ig)) + entr(ig,lmax(ig))356 IF (lmax(ig) > 0) THEN 357 detr(ig,lmax(ig)) = fm(ig,lmax(ig)) 320 358 fm(ig,lmax(ig)+1) = 0. 321 359 entr(ig,lmax(ig)) = 0. … … 329 367 IF (prt_level > 0) THEN 330 368 331 IF (ncorecdetr .ge.1) THEN369 IF (ncorecdetr > 0) THEN 332 370 print *, 'WARNING: Detrainment is greater than mass flux!' 333 371 print *, 'In', ncorecdetr, 'grid point(s) over', nlay, 'x', ngrid 334 372 ENDIF 335 373 336 IF (ncorecentr .ge.1) THEN374 IF (ncorecentr > 0) THEN 337 375 print *, 'WARNING: Entrained mass is greater than maximal authorized value!' 338 print *, ' in', ncorecentr, 'grid point(s) over', nlay, 'x', ngrid339 ENDIF 340 341 IF (ncorecfact .ge.1) THEN376 print *, 'In', ncorecentr, 'grid point(s) over', nlay, 'x', ngrid 377 ENDIF 378 379 IF (ncorecfact > 0) THEN 342 380 print *, 'WARNING: Entrained mass needs renormalization to be ok!' 343 print *, 'in', ncorecfact, 'grid point(s) over', nlay, 'x', ngrid 344 ENDIF 345 346 ! IF (nerrorequa.ge.1) THEN 381 print *, 'In', ncorecfact, 'grid point(s) over', nlay, 'x', ngrid 382 ! print *, 'WARNING: Entr fact:', fact 383 ENDIF 384 385 ! IF (nerrorequa > 0) THEN 347 386 ! print *, 'WARNING: !' 348 387 ! print *, 'in', nerrorequa, 'grid point(s) over', nlay, 'x', ngrid 349 388 ! ENDIF 350 389 351 IF (ncorecalpha .ge.1) THEN390 IF (ncorecalpha > 0) THEN 352 391 print *, 'WARNING: Updraft fraction is greater than maximal authorized value!' 353 print *, ' in', ncorecalpha, 'grid point(s) over', nlay, 'x', ngrid392 print *, 'In', ncorecalpha, 'grid point(s) over', nlay, 'x', ngrid 354 393 ENDIF 355 394 -
trunk/LMDZ.GENERIC/libf/phystd/thermcell_height.F90
r2143 r2177 2 2 ! 3 3 ! 4 SUBROUTINE thermcell_height(ngrid,nlay,lmin,l inter,lmix,lmax,zw2,&5 z lev,zmin,zmix,zmax,wmax,f_star)4 SUBROUTINE thermcell_height(ngrid,nlay,lmin,lmix,lmax,zlev, & 5 zmin,zmix,zmax,zw2,wmax,f_star) 6 6 7 7 … … 20 20 ! ------- 21 21 22 INTEGER ngrid, nlay 23 INTEGER lmin(ngrid) 22 INTEGER, INTENT(in) :: ngrid 23 INTEGER, INTENT(in) :: nlay 24 INTEGER, INTENT(in) :: lmin(ngrid) 24 25 25 REAL zlev(ngrid,nlay+1)26 REAL f_star(ngrid,nlay+1)26 REAL, INTENT(in) :: zlev(ngrid,nlay+1) 27 REAL, INTENT(in) :: f_star(ngrid,nlay+1) 27 28 28 29 ! Outputs: 29 30 ! -------- 30 31 31 INTEGER lmix(ngrid)32 INTEGER lmax(ngrid)32 INTEGER, INTENT(out) :: lmix(ngrid) 33 INTEGER, INTENT(out) :: lmax(ngrid) 33 34 34 REAL linter(ngrid)35 REAL zmin(ngrid) ! Altitude of the plume bottom36 REAL zmax(ngrid)! Altitude of the plume top37 REAL zmix(ngrid) ! Altitude of maximal vertical speed38 REAL wmax(ngrid) ! Maximal vertical speed39 REAL zw2(ngrid,nlay+1)! Square vertical speed (becomes its square root)35 REAL, INTENT(out) :: zmin(ngrid) ! Altitude of the plume bottom 36 REAL, INTENT(out) :: zmix(ngrid) ! Altitude of maximal vertical speed 37 REAL, INTENT(out) :: zmax(ngrid) ! Altitude of the plume top 38 REAL, INTENT(out) :: wmax(ngrid) ! Maximal vertical speed 39 40 REAL, INTENT(inout) :: zw2(ngrid,nlay+1) ! Square vertical speed (becomes its square root) 40 41 41 42 ! Local: 42 43 ! ------ 43 44 44 INTEGER ig, l45 INTEGER ig, k 45 46 46 REAL zlevinter(ngrid)47 REAL linter(ngrid) ! Level (continuous) of maximal vertical speed 47 48 48 49 !=============================================================================== … … 50 51 !=============================================================================== 51 52 52 DO ig=1,ngrid 53 lmax(ig) = lmin(ig) 54 lmix(ig) = lmin(ig) 55 ENDDO 53 linter(:) = lmin(:) 54 lmix(:) = lmin(:) 55 lmax(:) = lmin(:) 56 56 57 DO ig=1,ngrid 58 wmax(ig) = 0. 59 ENDDO 57 wmax(:) = 0. 60 58 61 DO ig=1,ngrid 62 zmax(ig) = 0. 63 zlevinter(ig) = zlev(ig,1) 64 ENDDO 59 zmin(:) = 0. 60 zmix(:) = 0. 61 zmax(:) = 0. 65 62 66 63 !=============================================================================== … … 72 69 !------------------------------------------------------------------------------- 73 70 74 DO l=1,nlay 75 DO ig=1,ngrid 76 zmin(ig) = zlev(ig,lmin(ig)) 77 ENDDO 71 DO ig=1,ngrid 72 zmin(ig) = zlev(ig,lmin(ig)) 78 73 ENDDO 79 74 … … 83 78 84 79 DO ig=1,ngrid 85 DO l=nlay,lmin(ig)+1,-186 IF ((zw2(ig, l) <= 0.).or.(f_star(ig,l) <= 0.)) THEN87 lmax(ig) = l- 180 DO k=nlay,lmin(ig)+1,-1 81 IF ((zw2(ig,k) <= 0.).or.(f_star(ig,k) <= 0.)) THEN 82 lmax(ig) = k - 1 88 83 ENDIF 89 84 ENDDO … … 96 91 IF (zw2(ig,nlay) > 0.) THEN 97 92 print *, 'WARNING: a thermal plume reaches the highest layer!' 98 print *, 'ig, l', ig, nlay93 print *, 'ig,k', ig, nlay 99 94 print *, 'zw2', zw2(ig,nlay) 100 95 lmax(ig) = nlay … … 103 98 104 99 !------------------------------------------------------------------------------- 105 ! Calcul de zmax continu via l e linter100 ! Calcul de zmax continu via linter 106 101 !------------------------------------------------------------------------------- 107 102 108 103 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 109 ! AB 110 ! 104 ! AB: lmin=lmax means the plume is not active and then zw2=0 at each level. 105 ! Otherwise we have lmax>lmin. 111 106 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 112 107 DO ig=1,ngrid 113 l= lmax(ig)114 IF ( l== lmin(ig)) THEN115 linter(ig) = l108 k = lmax(ig) 109 IF (k == lmin(ig)) THEN 110 linter(ig) = k 116 111 ELSE 117 linter(ig) = l - zw2(ig,l) / (zw2(ig,l+1) - zw2(ig,l))112 linter(ig) = k - zw2(ig,k) / (zw2(ig,k+1) - zw2(ig,k)) 118 113 ENDIF 119 ENDDO 120 121 DO ig=1,ngrid 122 zlevinter(ig) = zlev(ig,lmax(ig)) + (linter(ig) - lmax(ig)) & 123 & * (zlev(ig,lmax(ig)+1) - zlev(ig,lmax(ig))) 124 zmax(ig) = max(zmax(ig), zlevinter(ig)) 114 zmax(ig) = zlev(ig,lmax(ig)) + (linter(ig) - lmax(ig)) & 115 & * (zlev(ig,lmax(ig)+1) - zlev(ig,lmax(ig))) 125 116 ENDDO 126 117 … … 133 124 !------------------------------------------------------------------------------- 134 125 135 DO l=1,nlay126 DO k=1,nlay 136 127 DO ig=1,ngrid 137 IF(( l <= lmax(ig)).and.(l> lmin(ig))) THEN138 IF (zw2(ig, l).lt.0.) THEN128 IF((k <= lmax(ig)).and.(k > lmin(ig))) THEN 129 IF (zw2(ig,k) < 0.) THEN 139 130 print *, 'ERROR: zw2 has negative value(s)!' 140 print *, 'ig, l', ig, l141 print *, 'zw2', zw2(ig, l)131 print *, 'ig,k', ig, k 132 print *, 'zw2', zw2(ig,k) 142 133 CALL abort 143 134 ENDIF 144 135 145 136 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 146 ! AB 137 ! AB: WARNING zw2 becomes its square root! 147 138 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 148 zw2(ig, l) = sqrt(zw2(ig,l))139 zw2(ig,k) = sqrt(zw2(ig,k)) 149 140 150 IF (zw2(ig, l) > wmax(ig)) THEN151 wmax(ig) = zw2(ig, l)152 lmix(ig) = l141 IF (zw2(ig,k) > wmax(ig)) THEN 142 wmax(ig) = zw2(ig,k) 143 lmix(ig) = k 153 144 ENDIF 154 145 ELSE 155 zw2(ig, l) = 0.146 zw2(ig,k) = 0. 156 147 ENDIF 157 148 ENDDO … … 186 177 ENDDO 187 178 188 !===============================================================================189 ! Check consistence between zmax and zmix190 !===============================================================================191 192 179 !------------------------------------------------------------------------------- 193 ! 180 ! Check consistency between zmax and zmix 194 181 !------------------------------------------------------------------------------- 195 182 … … 199 186 print *, 'zmax,zmix_old,zmix_new', zmax(ig), zmix(ig), 0.9 * zmax(ig) 200 187 zmix(ig) = 0.9 * zmax(ig) 188 DO k=1,nlay 189 IF ((zmix(ig) >= zlev(ig,k)).and.(zmix(ig) < zlev(ig,k+1))) THEN 190 lmix(ig) = k 191 ENDIF 192 ENDDO 201 193 ENDIF 202 ENDDO203 204 !-------------------------------------------------------------------------------205 ! Calcul du nouveau lmix correspondant206 !-------------------------------------------------------------------------------207 208 DO ig=1,ngrid209 DO l=1,nlay210 IF ((zmix(ig) >= zlev(ig,l)).and.(zmix(ig) < zlev(ig,l+1))) THEN211 lmix(ig) = l212 ENDIF213 ENDDO214 194 ENDDO 215 195 -
trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90
r2144 r2177 42 42 ! New detr and entre formulae (no longer alimentation) 43 43 ! lmin can be greater than 1 44 ! Mix every tracer (EN COURS) 45 ! Old version of thermcell_dq is removed 46 ! Alternative version thermcell_dv2 is removed 44 ! Mix every tracer 47 45 ! 48 46 !=============================================================================== … … 61 59 ! ------- 62 60 63 INTEGER ngrid, nlay, nq 64 65 REAL ptimestep 66 REAL pplay(ngrid,nlay) ! Layer pressure 67 REAL pplev(ngrid,nlay+1) ! Level pressure 68 REAL pphi(ngrid,nlay) ! Geopotential 69 70 REAL pu(ngrid,nlay) ! Zonal wind 71 REAL pv(ngrid,nlay) ! Meridional wind 72 REAL pt(ngrid,nlay) ! Temperature 73 REAL pq(ngrid,nlay,nq) ! Tracers mass mixing ratio 74 75 LOGICAL firstcall 61 INTEGER, INTENT(in) :: ngrid 62 INTEGER, INTENT(in) :: nlay 63 INTEGER, INTENT(in) :: nq 64 65 REAL, INTENT(in) :: ptimestep 66 REAL, INTENT(in) :: pplay(ngrid,nlay) ! Layer pressure 67 REAL, INTENT(in) :: pplev(ngrid,nlay+1) ! Level pressure 68 REAL, INTENT(in) :: pphi(ngrid,nlay) ! Geopotential 69 70 REAL, INTENT(in) :: pu(ngrid,nlay) ! Zonal wind 71 REAL, INTENT(in) :: pv(ngrid,nlay) ! Meridional wind 72 REAL, INTENT(in) :: pt(ngrid,nlay) ! Temperature 73 REAL, INTENT(in) :: pq(ngrid,nlay,nq) ! Tracers mass mixing ratio 74 75 LOGICAL, INTENT(in) :: firstcall 76 76 77 77 ! Outputs: 78 78 ! -------- 79 79 80 REAL pduadj(ngrid,nlay) ! u convective variations 81 REAL pdvadj(ngrid,nlay) ! v convective variations 82 REAL pdtadj(ngrid,nlay) ! t convective variations 83 REAL pdqadj(ngrid,nlay,nq) ! q convective variations 84 85 REAL f0(ngrid) ! mass flux norm (after possible time relaxation) 86 REAL fm0(ngrid,nlay+1) ! mass flux (after possible time relaxation) 87 REAL entr0(ngrid,nlay) ! entrainment (after possible time relaxation) 88 REAL detr0(ngrid,nlay) ! detrainment (after possible time relaxation) 80 INTEGER, INTENT(out) :: lmax(ngrid) ! Highest layer reached by the plume 81 INTEGER, INTENT(out) :: lmix(ngrid) ! Layer in which plume vertical speed is maximal 82 INTEGER, INTENT(out) :: lmin(ngrid) ! First unstable layer 83 84 REAL, INTENT(out) :: pduadj(ngrid,nlay) ! u convective variations 85 REAL, INTENT(out) :: pdvadj(ngrid,nlay) ! v convective variations 86 REAL, INTENT(out) :: pdtadj(ngrid,nlay) ! t convective variations 87 REAL, INTENT(out) :: pdqadj(ngrid,nlay,nq) ! q convective variations 88 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) 89 93 90 94 ! Local: … … 92 96 93 97 INTEGER ig, k, l, iq 94 INTEGER lmax(ngrid) ! Highest layer reached by the plume 95 INTEGER lmix(ngrid) ! Layer in which plume vertical speed is maximal 96 INTEGER lmin(ngrid) ! First unstable layer 97 98 REAL zmix(ngrid) ! Altitude of maximal vertical speed 99 REAL zmax(ngrid) ! Maximal altitudes where plumes are active 100 REAL zmin(ngrid) ! Minimal altitudes where plumes are active 101 102 REAL zlay(ngrid,nlay) ! Layers altitudes 103 REAL zlev(ngrid,nlay+1) ! Levels altitudes 104 REAL rho(ngrid,nlay) ! Layers densities 105 REAL rhobarz(ngrid,nlay) ! Levels densities 106 REAL masse(ngrid,nlay) ! Layers masses 107 REAL zpopsk(ngrid,nlay) ! Exner function 108 109 REAL zu(ngrid,nlay) ! u environment 110 REAL zv(ngrid,nlay) ! v environment 111 REAL zt(ngrid,nlay) ! TR environment 112 REAL zqt(ngrid,nlay) ! qt environment 113 REAL zql(ngrid,nlay) ! ql environment 114 REAL zhl(ngrid,nlay) ! TP environment 115 REAL ztv(ngrid,nlay) ! TRPV environment 116 REAL zqs(ngrid,nlay) ! qsat environment 117 118 REAL zua(ngrid,nlay) ! u plume 119 REAL zva(ngrid,nlay) ! v plume 120 REAL zta(ngrid,nlay) ! TR plume 121 REAL zqla(ngrid,nlay) ! qv plume 122 REAL zqta(ngrid,nlay) ! qt plume 123 REAL zhla(ngrid,nlay) ! TP plume 124 REAL ztva(ngrid,nlay) ! TRPV plume 125 REAL zqsa(ngrid,nlay) ! qsat plume 126 127 REAL zqa(ngrid,nlay,nq) ! q plume (ql=0, qv=qt) 128 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 133 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 140 REAL linter(ngrid) ! Level (continuous) of maximal vertical speed 141 REAL wmax(ngrid) ! Maximal vertical speed 142 REAL zw2(ngrid,nlay+1) ! Plume vertical speed 143 REAL zdthladj(ngrid,nlay) ! Potential temperature variations 144 REAL dummy(ngrid,nlay) ! Dummy argument for thermcell_dq() 98 99 REAL zmix(ngrid) ! Altitude of maximal vertical speed 100 REAL zmax(ngrid) ! Maximal altitudes where plumes are active 101 REAL zmin(ngrid) ! Minimal altitudes where plumes are active 102 103 REAL zlay(ngrid,nlay) ! Layers altitudes 104 REAL zlev(ngrid,nlay+1) ! Levels altitudes 105 REAL rho(ngrid,nlay) ! Layers densities 106 REAL rhobarz(ngrid,nlay) ! Levels densities 107 REAL masse(ngrid,nlay) ! Layers masses 108 REAL zpopsk(ngrid,nlay) ! Exner function 109 110 REAL zu(ngrid,nlay) ! u environment 111 REAL zv(ngrid,nlay) ! v environment 112 REAL zt(ngrid,nlay) ! TR environment 113 REAL zqt(ngrid,nlay) ! qt environment 114 REAL zql(ngrid,nlay) ! ql environment 115 REAL zhl(ngrid,nlay) ! TP environment 116 REAL ztv(ngrid,nlay) ! TRPV environment 117 REAL zqs(ngrid,nlay) ! qsat environment 118 119 REAL zua(ngrid,nlay) ! u plume 120 REAL zva(ngrid,nlay) ! v plume 121 REAL zqla(ngrid,nlay) ! qv plume 122 REAL zqta(ngrid,nlay) ! qt plume 123 REAL zhla(ngrid,nlay) ! TP plume 124 REAL ztva(ngrid,nlay) ! TRPV plume 125 REAL zqsa(ngrid,nlay) ! qsat plume 126 127 REAL zqa(ngrid,nlay,nq) ! q plume (ql=0, qv=qt) 128 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 133 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 140 REAL wmax(ngrid) ! Maximal vertical speed 141 REAL zw2(ngrid,nlay+1) ! Plume vertical speed 142 REAL zdthladj(ngrid,nlay) ! Potential temperature variations 143 REAL dummy(ngrid,nlay) ! Dummy argument for thermcell_dq() 145 144 146 145 !=============================================================================== … … 154 153 ENDIF 155 154 156 f_star(:,:) = 0. 157 entr_star(:,:) = 0. 158 detr_star(:,:) = 0. 159 160 f(:) = 0. 161 162 fm(:,:) = 0. 163 entr(:,:) = 0. 164 detr(:,:) = 0. 165 166 lmax(:) = 1 167 lmix(:) = 1 168 lmin(:) = 1 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 169 159 170 160 pduadj(:,:) = 0.0 … … 173 163 pdqadj(:,:,:) = 0.0 174 164 175 DO ig=1,ngrid 176 ! AB: Careful: Hard-coded value from Earth version! 177 ! f0(ig) = max(f0(ig), 1.e-2) 178 ! AB: No pescribed minimal value for f0 179 f0(ig) = max(f0(ig), 0.) 180 ENDDO 165 zdthladj(:,:) = 0.0 181 166 182 167 !=============================================================================== … … 212 197 rho(:,:) = pplay(:,:) / (zpopsk(:,:) * RD * ztv(:,:)) 213 198 199 rhobarz(:,1) = rho(:,1) 214 200 IF (prt_level.ge.10) THEN 215 201 print *, 'WARNING: density in the first layer is equal to density at the first level!' 216 print *, 'rhobarz(:,1)', rhobarz(:,1) 217 print *, 'rho(:,1) ', rho(:,1) 218 ENDIF 219 220 rhobarz(:,1) = rho(:,1) 202 ENDIF 221 203 222 204 DO l=2,nlay … … 244 226 ! 245 227 ! --------------------------- 246 ! 247 ! ----- F_lmax+1=0 ------zmax 248 ! lmax 249 ! ------F_lmax>0------------- 250 ! 251 ! --------------------------- 252 ! 253 ! --------------------------- 254 ! 255 ! ------------------wmax,zmix 256 ! lmix 257 ! --------------------------- 258 ! 259 ! --------------------------- 260 ! 261 ! --------------------------- 262 ! 228 ! _ 229 ! ----- F_lmax+1=0 ------zmax \ 230 ! lmax | 231 ! ------F_lmax>0------------- | 232 ! | 233 ! --------------------------- | 234 ! | 235 ! --------------------------- | 236 ! | 237 ! ------------------wmax,zmix | 238 ! lmix | 239 ! --------------------------- | 240 ! | 241 ! --------------------------- | 242 ! | E, D 243 ! --------------------------- | 244 ! | 263 245 ! --------------------------- rhobarz, f_star, fm, fm0, zw2, fraca 264 ! zt, zu, zv, zo, rho 265 ! --------------------------- 266 ! 267 ! --------------------------- 268 ! 269 ! --------------------------- 270 ! 271 ! ------F_lmin+1>0----------- 272 ! lmin 273 ! ----- F_lmin=0 ------------ 246 ! zt, zu, zv, zo, rho | 247 ! --------------------------- | 248 ! | 249 ! --------------------------- | 250 ! | 251 ! --------------------------- | 252 ! | 253 ! ------F_lmin+1>0----------- | 254 ! lmin | 255 ! ----- F_lmin=0 ------------ _/ 274 256 ! 275 257 ! --------------------------- … … 309 291 !------------------------------------------------------------------------------- 310 292 311 CALL thermcell_plume(ngrid,nlay,nq,ptimestep, ztv,&312 & z hl,zqt,zql,rhobarz,zlev,pplev,pphi,zpopsk,&293 CALL thermcell_plume(ngrid,nlay,nq,ptimestep, & 294 & ztv,zhl,zqt,zql,zlev,pplev,zpopsk, & 313 295 & detr_star,entr_star,f_star, & 314 & ztva,zhla,zq la,zqta,zta,zqsa,&315 & zw2,lmi x,lmin)296 & ztva,zhla,zqta,zqla,zqsa, & 297 & zw2,lmin) 316 298 317 299 !------------------------------------------------------------------------------- … … 320 302 321 303 ! AB: Careful, zw2 became its square root in thermcell_height! 322 CALL thermcell_height(ngrid,nlay,lmin,l inter,lmix,lmax,zw2,&323 & z lev,zmin,zmix,zmax,wmax,f_star)304 CALL thermcell_height(ngrid,nlay,lmin,lmix,lmax,zlev, & 305 & zmin,zmix,zmax,zw2,wmax,f_star) 324 306 325 307 !=============================================================================== … … 341 323 ENDIF 342 324 343 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 344 ! Test valable seulement en 1D mais pas genant 345 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 325 ! FH: Test valable seulement en 1D mais pas genant 346 326 IF (.not. (f0(1).ge.0.) ) THEN 347 327 print *, 'ERROR: mass flux norm is not positive!' … … 403 383 404 384 CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse, & 405 & zhl,zdthladj,dummy,lmin )385 & zhl,zdthladj,dummy,lmin,lmax) 406 386 407 387 DO l=1,nlay … … 417 397 DO iq=1,nq 418 398 CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse, & 419 & pq(:,:,iq),pdqadj(:,:,iq),zqa(:,:,iq),lmin )399 & pq(:,:,iq),pdqadj(:,:,iq),zqa(:,:,iq),lmin,lmax) 420 400 ENDDO 421 401 … … 429 409 ELSE 430 410 CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse, & 431 & zu,pduadj,zua,lmin )411 & zu,pduadj,zua,lmin,lmax) 432 412 CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse, & 433 & zv,pdvadj,zva,lmin )413 & zv,pdvadj,zva,lmin,lmax) 434 414 ENDIF 435 415
Note: See TracChangeset
for help on using the changeset viewer.
