Index: /trunk/LMDZ.GENERIC/README =================================================================== --- /trunk/LMDZ.GENERIC/README (revision 2176) +++ /trunk/LMDZ.GENERIC/README (revision 2177) @@ -1495,3 +1495,5 @@ == 12/11/2019 == AB - fm0, entr0 and detr0 are now allocatable variables in physiq_mod. That is necessary if tau_thermals > 0. - +- In thermcell_flux, "bidouilles" are modified: now the plumes stop when the updraft fraction is greater than alpha_max ; + e > e_max is no longer permitted ; b <= incoming mass flux is checked last +- Cleanup thermal plums model subroutines (thermcell_main, thermcell_env, thercell_dq, thermcell_dv2, thermcell_closure, thermcell_height) Index: /trunk/LMDZ.GENERIC/libf/phystd/thermcell_closure.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/thermcell_closure.F90 (revision 2176) +++ /trunk/LMDZ.GENERIC/libf/phystd/thermcell_closure.F90 (revision 2177) @@ -29,19 +29,20 @@ ! ------- - INTEGER ngrid, nlay - INTEGER lmax(ngrid) + INTEGER, INTENT(in) :: ngrid + INTEGER, INTENT(in) :: nlay + INTEGER, INTENT(in) :: lmax(ngrid) - REAL ptimestep - REAL rho(ngrid,nlay) - REAL zlev(ngrid,nlay) - REAL alim_star(ngrid,nlay) - REAL zmin(ngrid) - REAL zmax(ngrid) - REAL wmax(ngrid) + REAL, INTENT(in) :: ptimestep + REAL, INTENT(in) :: rho(ngrid,nlay) + REAL, INTENT(in) :: zlev(ngrid,nlay) + REAL, INTENT(in) :: alim_star(ngrid,nlay) + REAL, INTENT(in) :: zmin(ngrid) + REAL, INTENT(in) :: zmax(ngrid) + REAL, INTENT(in) :: wmax(ngrid) ! Outputs: ! -------- - REAL f(ngrid) + REAL, INTENT(out) :: f(ngrid) ! Local: @@ -90,9 +91,7 @@ DO l=1,llmax-1 DO ig=1,ngrid - IF (l < lmax(ig)) THEN alim_star2(ig) = alim_star2(ig) + alim_star(ig,l)**2 & - & / (rho(ig,l) * (zlev(ig,l+1) - zlev(ig,l))) ! => intergration is ok because alim_star = a* dz + & / (rho(ig,l) * (zlev(ig,l+1) - zlev(ig,l))) ! => integration is ok because alim_star = a* dz alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig,l) - ENDIF ENDDO ENDDO Index: /trunk/LMDZ.GENERIC/libf/phystd/thermcell_dq.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/thermcell_dq.F90 (revision 2176) +++ /trunk/LMDZ.GENERIC/libf/phystd/thermcell_dq.F90 (revision 2177) @@ -3,5 +3,5 @@ ! SUBROUTINE thermcell_dq(ngrid,nlay,ptimestep,fm,entr,detr,masse, & - q,dq,qa,lmin) + q,dq,qa,lmin,lmax) @@ -16,6 +16,6 @@ ! ! Modif 2019/04 (AB alexandre.boissinot@lmd.jussieu.fr) -! dqimpl = 1 : implicit scheme -! dqimpl = 0 : explicit scheme +! dqimpl = true : implicit scheme +! dqimpl = false : explicit scheme ! !=============================================================================== @@ -34,19 +34,21 @@ ! ------- - INTEGER ngrid, nlay - INTEGER lmin(ngrid) + INTEGER, INTENT(in) :: ngrid + INTEGER, INTENT(in) :: nlay + INTEGER, INTENT(in) :: lmin(ngrid) + INTEGER, INTENT(in) :: lmax(ngrid) - REAL ptimestep - REAL masse(ngrid,nlay) - REAL fm(ngrid,nlay+1) - REAL entr(ngrid,nlay) - REAL detr(ngrid,nlay) - REAL q(ngrid,nlay) + REAL, INTENT(in) :: ptimestep + REAL, INTENT(in) :: masse(ngrid,nlay) + REAL, INTENT(in) :: fm(ngrid,nlay+1) + REAL, INTENT(in) :: entr(ngrid,nlay) + REAL, INTENT(in) :: detr(ngrid,nlay) ! Outputs: ! -------- - REAL dq(ngrid,nlay) - REAL qa(ngrid,nlay) + REAL, INTENT(inout) :: q(ngrid,nlay) + REAL, INTENT(out) :: dq(ngrid,nlay) + REAL, INTENT(out) :: qa(ngrid,nlay) ! Local: @@ -82,8 +84,8 @@ cfl = max(cfl, fm(ig,l) / zzm) - IF (entr(ig,l).gt.zzm) THEN + IF (entr(ig,l) > zzm) THEN print *, 'ERROR: entrainment is greater than the layer mass!' print *, 'ig,l,entr', ig, l, entr(ig,l) - print *, 'lmin', lmin(ig) + print *, 'lmin,lmax', lmin(ig), lmax(ig) print *, '-------------------------------' print *, 'entr*dt,mass', entr(ig,l)*ptimestep, masse(ig,l) Index: /trunk/LMDZ.GENERIC/libf/phystd/thermcell_dv2.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/thermcell_dv2.F90 (revision 2176) +++ /trunk/LMDZ.GENERIC/libf/phystd/thermcell_dv2.F90 (revision 2177) @@ -16,5 +16,5 @@ !=============================================================================== - USE print_control_mod, ONLY: prt_level,lunout + USE print_control_mod, ONLY: prt_level, lunout IMPLICIT none @@ -28,26 +28,25 @@ ! ------- - INTEGER ngrid, nlay + INTEGER, INTENT(in) :: ngrid + INTEGER, INTENT(in) :: nlay - REAL ptimestep - REAL masse(ngrid,nlay) - REAL fm(ngrid,nlay+1) - REAL entr(ngrid,nlay) - REAL detr(ngrid,nlay) - REAL fraca(ngrid,nlay+1) - REAL zmax(ngrid) - REAL zmin(ngrid) - REAL u(ngrid,nlay) - REAL v(ngrid,nlay) + REAL, INTENT(in) :: ptimestep + REAL, INTENT(in) :: masse(ngrid,nlay) + REAL, INTENT(in) :: fm(ngrid,nlay+1) + REAL, INTENT(in) :: entr(ngrid,nlay) + REAL, INTENT(in) :: detr(ngrid,nlay) + REAL, INTENT(in) :: fraca(ngrid,nlay+1) + REAL, INTENT(in) :: zmax(ngrid) + REAL, INTENT(in) :: zmin(ngrid) + REAL, INTENT(in) :: u(ngrid,nlay) + REAL, INTENT(in) :: v(ngrid,nlay) ! Outputs: ! -------- - REAL dua(ngrid,nlay) - REAL dva(ngrid,nlay) - REAL ua(ngrid,nlay) - REAL va(ngrid,nlay) - REAL du(ngrid,nlay) - REAL dv(ngrid,nlay) + REAL, INTENT(out) :: ua(ngrid,nlay) + REAL, INTENT(out) :: va(ngrid,nlay) + REAL, INTENT(out) :: du(ngrid,nlay) + REAL, INTENT(out) :: dv(ngrid,nlay) ! Local: @@ -66,4 +65,6 @@ REAL ue(ngrid,nlay) REAL ve(ngrid,nlay) + REAL dua(ngrid,nlay) + REAL dva(ngrid,nlay) REAL plume_height(ngrid) @@ -121,8 +122,8 @@ DO iter=1,5 DO ig=1,ngrid -! Calcul prenant en compte la fraction reelle +! FH: Calcul prenant en compte la fraction reelle zf = 0.5 * (fraca(ig,l) + fraca(ig,l+1)) zf2 = 1./(1.-zf) -! Calcul avec fraction infiniement petite +! FH: Calcul avec fraction infiniement petite ! zf = 0. ! zf2 = 1. @@ -133,7 +134,7 @@ !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF (ltherm(ig,l)) THEN -! On choisit une relaxation lineaire. +! FH: On choisit une relaxation lineaire. ! gamma(ig,l) = gamma0(ig,l) -! On choisit une relaxation quadratique. +! FH: On choisit une relaxation quadratique. gamma(ig,l) = gamma0(ig,l) * sqrt(dua(ig,l)**2 + dva(ig,l)**2) Index: /trunk/LMDZ.GENERIC/libf/phystd/thermcell_env.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/thermcell_env.F90 (revision 2176) +++ /trunk/LMDZ.GENERIC/libf/phystd/thermcell_env.F90 (revision 2177) @@ -38,32 +38,34 @@ ! ------- - INTEGER ngrid, nlay, nq + INTEGER, INTENT(in) :: ngrid + INTEGER, INTENT(in) :: nlay + INTEGER, INTENT(in) :: nq - REAL pq(ngrid,nlay,nq) ! Large scale water - REAL pt(ngrid,nlay) ! Large scale temperature - REAL pu(ngrid,nlay) ! Large scale zonal wind - REAL pv(ngrid,nlay) ! Large scale meridional wind - REAL pplay(ngrid,nlay) ! Layers pressure - REAL pplev(ngrid,nlay+1) ! Levels pressure - REAL zpopsk(ngrid,nlay) ! Exner function + REAL, INTENT(in) :: pq(ngrid,nlay,nq) ! Large scale water + REAL, INTENT(in) :: pt(ngrid,nlay) ! Large scale temperature + REAL, INTENT(in) :: pu(ngrid,nlay) ! Large scale zonal wind + REAL, INTENT(in) :: pv(ngrid,nlay) ! Large scale meridional wind + REAL, INTENT(in) :: pplay(ngrid,nlay) ! Layers pressure + REAL, INTENT(in) :: pplev(ngrid,nlay+1) ! Levels pressure + REAL, INTENT(in) :: zpopsk(ngrid,nlay) ! Exner function ! Outputs: ! -------- - REAL zqt(ngrid,nlay) ! qt environment - REAL zql(ngrid,nlay) ! ql environment - REAL zt(ngrid,nlay) ! T environment - REAL ztv(ngrid,nlay) ! TRPV environment - REAL zhl(ngrid,nlay) ! TP environment - REAL zu(ngrid,nlay) ! u environment - REAL zv(ngrid,nlay) ! v environment - REAL zqs(ngrid,nlay) ! qsat environment + REAL, INTENT(out) :: zt(ngrid,nlay) ! T environment + REAL, INTENT(out) :: ztv(ngrid,nlay) ! TRPV environment + REAL, INTENT(out) :: zhl(ngrid,nlay) ! TP environment + REAL, INTENT(out) :: zu(ngrid,nlay) ! u environment + REAL, INTENT(out) :: zv(ngrid,nlay) ! v environment + REAL, INTENT(out) :: zqt(ngrid,nlay) ! qt environment + REAL, INTENT(out) :: zql(ngrid,nlay) ! ql environment + REAL, INTENT(out) :: zqs(ngrid,nlay) ! qsat environment ! Local: ! ------ - INTEGER ig, l + INTEGER ig, k - REAL psat ! Dummy argument for Psat_water() + REAL psat ! Dummy argument for Psat_water() !=============================================================================== @@ -85,16 +87,16 @@ IF (water) THEN - DO l=1,nlay + DO k=1,nlay DO ig=1,ngrid - CALL Psat_water(pt(ig,l), pplev(ig,l), psat, zqs(ig,l)) + CALL Psat_water(pt(ig,k), pplev(ig,k), psat, zqs(ig,k)) ENDDO ENDDO - DO l=1,nlay + DO k=1,nlay DO ig=1,ngrid - zql(ig,l) = max(0.,pq(ig,l,igcm_h2o_vap) - zqs(ig,l)) - zt(ig,l) = pt(ig,l) + RLvCp * zql(ig,l) - ztv(ig,l) = zt(ig,l) / zpopsk(ig,l) & - & * (1. + RETV * (zqt(ig,l)-zql(ig,l)) - zql(ig,l)) + zql(ig,k) = max(0.,pq(ig,k,igcm_h2o_vap) - zqs(ig,k)) + zt(ig,k) = pt(ig,k) + RLvCp * zql(ig,k) + ztv(ig,k) = zt(ig,k) / zpopsk(ig,k) & + & * (1. + RETV * (zqt(ig,k)-zql(ig,k)) - zql(ig,k)) ENDDO ENDDO Index: /trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90 (revision 2176) +++ /trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90 (revision 2177) @@ -27,23 +27,25 @@ ! ------- - INTEGER ngrid, nlay - INTEGER lmin(ngrid) - INTEGER lmax(ngrid) - - REAL entr_star(ngrid,nlay) - REAL detr_star(ngrid,nlay) - REAL zw2(ngrid,nlay+1) - REAL zlev(ngrid,nlay+1) - REAL masse(ngrid,nlay) - REAL ptimestep - REAL rhobarz(ngrid,nlay) - REAL f(ngrid) + INTEGER, INTENT(in) :: ngrid + INTEGER, INTENT(in) :: nlay + INTEGER, INTENT(in) :: lmin(ngrid) + + REAL, INTENT(in) :: entr_star(ngrid,nlay) + REAL, INTENT(in) :: detr_star(ngrid,nlay) + REAL, INTENT(in) :: zw2(ngrid,nlay+1) + REAL, INTENT(in) :: zlev(ngrid,nlay+1) + REAL, INTENT(in) :: masse(ngrid,nlay) + REAL, INTENT(in) :: ptimestep + REAL, INTENT(in) :: rhobarz(ngrid,nlay) + REAL, INTENT(in) :: f(ngrid) ! Outputs: -! -------- - - REAL entr(ngrid,nlay) - REAL detr(ngrid,nlay) - REAL fm(ngrid,nlay+1) +! -------- + + INTEGER, INTENT(inout) :: lmax(ngrid) + + REAL, INTENT(out) :: entr(ngrid,nlay) + REAL, INTENT(out) :: detr(ngrid,nlay) + REAL, INTENT(out) :: fm(ngrid,nlay+1) ! Local: @@ -53,11 +55,11 @@ INTEGER igout, lout ! Error grid point and level - REAL zfm ! Mass flux such as updraft fraction is equal to its maximal authorized value alpha_max - REAL f_old ! Save initial value of mass flux + REAL fmax ! Maximal authorized mass flux (alpha < alpha_max) + REAL fff0 ! Save initial value of mass flux + REAL emax ! Maximal authorized entrainment (entr*dt < mass_max) REAL eee0 ! Save initial value of entrainment REAL ddd0 ! Save initial value of detrainment REAL eee ! eee0 - layer mass * maximal authorized mass fraction / dt REAL ddd ! ddd0 - eee - REAL zzz ! Temporary variable set to mass flux REAL fact REAL test @@ -103,5 +105,5 @@ !------------------------------------------------------------------------------- -! Calcul du flux de masse +! Mass flux !------------------------------------------------------------------------------- @@ -112,4 +114,5 @@ ELSEIF (l == lmax(ig)) THEN fm(ig,l+1) = 0. + entr(ig,l) = 0. detr(ig,l) = fm(ig,l) + entr(ig,l) ELSE @@ -170,36 +173,15 @@ print *, '---------------------------------------------------------' print *, 'ERROR: mass flux has negative value(s)!' - print *, 'ig,l,entr', igout, lout, entr(igout,lout) + print *, 'ig,l,norm', igout, lout, f(igout) print *, 'lmin,lmax', lmin(igout), lmax(igout) print *, '- - - - - - - - - - - - - - - - - - - - - - - - - - - - -' DO k=nlay,1,-1 - print *, 'fm ', fm(igout,k+1) + print *, 'k, fm ', k+1, fm(igout,k+1) print *, 'entr,detr', entr(igout,k), detr(igout,k) ENDDO - print *, 'fm ', fm(igout,1) + print *, 'k, fm ', 1, fm(igout,1) print *, '---------------------------------------------------------' CALL abort ENDIF - -!------------------------------------------------------------------------------- -! Is detrainment lessser than incoming mass flux ? -!------------------------------------------------------------------------------- - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -! AB : Even if fm has no negative value, it can be lesser than detr. -! That's not suitable because when we will mix the plume with the -! environment, it will detrain more mass than it is physically able to do. -! When it occures, that imply that entr + fm is greater than detr, -! otherwise we get a negative outgoing mass flux (cf. above). -! That is why we correct the detrainment as follows. -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - DO ig=1,ngrid - IF (detr(ig,l).gt.fm(ig,l)) THEN - detr(ig,l) = fm(ig,l) - entr(ig,l) = fm(ig,l+1) - ncorecdetr = ncorecdetr + 1 - ENDIF - ENDDO !------------------------------------------------------------------------------- @@ -213,5 +195,5 @@ ! If it's not enough, we can increase entr in the layer above and decrease ! the outgoing mass flux in the current layer. -! If it's still unsufficient, code will abort (now commented). +! If it's still insufficient, code will abort (now commented). ! Else we reset entr to its intial value and we renormalize entrainment, ! detrainment and mass flux profiles such as we do not exceed the maximal @@ -222,8 +204,13 @@ eee0 = entr(ig,l) ddd0 = detr(ig,l) - eee = entr(ig,l) - masse(ig,l) * fomass_max / ptimestep - ddd = detr(ig,l) - eee - IF (eee > 0.) THEN - entr(ig,l) = entr(ig,l) - eee + emax = masse(ig,l) * fomass_max / ptimestep + IF (emax < 0.) THEN + print *, 'ERROR: layer mass is negative!' + print *, 'mass,emax', masse(ig,l), emax + print *, 'ig,l', ig, l + ENDIF + IF (eee0 > emax) THEN + entr(ig,l) = emax + ddd = ddd0 + emax - eee0 ncorecentr = ncorecentr + 1 IF (ddd > 0.) THEN @@ -235,13 +222,27 @@ fm(ig,l+1) = fm(ig,l) + entr(ig,l) entr(ig,l+1) = entr(ig,l+1) + ddd +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +! AB: Simulation abort (try to reduce the physical time step). +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ELSE + entr(ig,l) = entr(ig,l) + eee + igout = ig + lout = l + labort_physic = .true. +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +! AB: We can renormalize the plume mass fluxes. I think it does not work. +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ELSE -! entr(ig,l) = entr(ig,l) + eee -! igout = ig -! lout = l -! labort_physic = .true. - ELSE - fact = max(fact, eee0 / (eee0 - eee)) - entr(ig,l) = eee0 - ncorecfact = ncorecfact + 1 +! fact = max(fact, eee0 / emax) +! entr(ig,l) = eee0 +! ncorecfact = ncorecfact + 1 +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +! AB: The renormalization can be just applied in the local plume. +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +! DO k=lmin(ig),lmax(ig) +! entr(ig,k) = entr(ig,k) * emax / eee0 +! detr(ig,k) = detr(ig,k) * emax / eee0 +! fm(ig,k) = fm(ig,k) * emax / eee0 +! ENDDO ENDIF ENDIF @@ -272,19 +273,53 @@ !------------------------------------------------------------------------------- -! Is Updraft fraction lesser than alpha_max ? -!------------------------------------------------------------------------------- - - DO ig=1,ngrid - zfm = rhobarz(ig,l+1) * zw2(ig,l+1) * alpha_max - - IF (fm(ig,l+1) > zfm) THEN - f_old = fm(ig,l+1) - fm(ig,l+1) = zfm - detr(ig,l) = detr(ig,l) + f_old - fm(ig,l+1) +! Is updraft fraction lesser than alpha_max ? +!------------------------------------------------------------------------------- + + DO ig=1,ngrid + fff0 = fm(ig,l+1) + fmax = rhobarz(ig,l+1) * zw2(ig,l+1) * alpha_max +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +! AB: The plume mass flux can be reduced. +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +! IF (fff0 > fmax) THEN +! fm(ig,l+1) = fmax +! detr(ig,l) = detr(ig,l) + fff0 - fmax +! ncorecalpha = ncorecalpha + 1 +! entr(ig,l+1) = entr(ig,l+1) + fff0 - fmax +! ENDIF +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +! AB: The plume can be stopped here. +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + IF (fff0 > fmax) THEN ncorecalpha = ncorecalpha + 1 -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -! AB : The previous change doesn't observe the equation df/dz = e - d. -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - entr(ig,l+1) = entr(ig,l+1) + f_old - fm(ig,l+1) + DO k=l+1,lmax(ig) + entr(ig,k) = 0. + detr(ig,k) = 0. + fm(ig,k) = 0. + ENDDO + lmax(ig) = l + entr(ig,l) = 0. + detr(ig,l) = fm(ig,l) + ENDIF + ENDDO + +!------------------------------------------------------------------------------- +! Is detrainment lesser than incoming mass flux ? +!------------------------------------------------------------------------------- + +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +! AB : Even if fm has no negative value, it can be lesser than detr. +! That is not suitable because when we will mix the plume with the +! environment, it will detrain more mass than it is physically able to do. +! When it occures, that imply that entr + fm is greater than detr, +! otherwise we get a negative outgoing mass flux (cf. above). +! That is why we decrease entrainment and detrainment as follows. +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + DO ig=1,ngrid + IF (detr(ig,l) > fm(ig,l)) THEN + detr(ig,l) = fm(ig,l) + entr(ig,l) = fm(ig,l+1) + ncorecdetr = ncorecdetr + 1 ENDIF ENDDO @@ -292,9 +327,12 @@ ENDDO - IF (fact > 0.) THEN - entr(:,:) = entr(:,:) / fact - detr(:,:) = detr(:,:) / fact - fm(:,:) = fm(:,:) / fact - ENDIF +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +! AB: The renormalization can be applied everywhere. +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +! IF (fact > 0.) THEN +! entr(:,:) = entr(:,:) / fact +! detr(:,:) = detr(:,:) / fact +! fm(:,:) = fm(:,:) / fact +! ENDIF !------------------------------------------------------------------------------- @@ -316,6 +354,6 @@ DO ig=1,ngrid - IF (lmax(ig).GT.0) THEN - detr(ig,lmax(ig)) = fm(ig,lmax(ig)) + entr(ig,lmax(ig)) + IF (lmax(ig) > 0) THEN + detr(ig,lmax(ig)) = fm(ig,lmax(ig)) fm(ig,lmax(ig)+1) = 0. entr(ig,lmax(ig)) = 0. @@ -329,27 +367,28 @@ IF (prt_level > 0) THEN - IF (ncorecdetr.ge.1) THEN + IF (ncorecdetr > 0) THEN print *, 'WARNING: Detrainment is greater than mass flux!' print *, 'In', ncorecdetr, 'grid point(s) over', nlay, 'x', ngrid ENDIF - IF (ncorecentr.ge.1) THEN + IF (ncorecentr > 0) THEN print *, 'WARNING: Entrained mass is greater than maximal authorized value!' - print *, 'in', ncorecentr, 'grid point(s) over', nlay, 'x', ngrid - ENDIF - - IF (ncorecfact.ge.1) THEN + print *, 'In', ncorecentr, 'grid point(s) over', nlay, 'x', ngrid + ENDIF + + IF (ncorecfact > 0) THEN print *, 'WARNING: Entrained mass needs renormalization to be ok!' - print *, 'in', ncorecfact, 'grid point(s) over', nlay, 'x', ngrid - ENDIF - -! IF (nerrorequa.ge.1) THEN + print *, 'In', ncorecfact, 'grid point(s) over', nlay, 'x', ngrid +! print *, 'WARNING: Entr fact:', fact + ENDIF + +! IF (nerrorequa > 0) THEN ! print *, 'WARNING: !' ! print *, 'in', nerrorequa, 'grid point(s) over', nlay, 'x', ngrid ! ENDIF - IF (ncorecalpha.ge.1) THEN + IF (ncorecalpha > 0) THEN print *, 'WARNING: Updraft fraction is greater than maximal authorized value!' - print *, 'in', ncorecalpha, 'grid point(s) over', nlay, 'x', ngrid + print *, 'In', ncorecalpha, 'grid point(s) over', nlay, 'x', ngrid ENDIF Index: /trunk/LMDZ.GENERIC/libf/phystd/thermcell_height.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/thermcell_height.F90 (revision 2176) +++ /trunk/LMDZ.GENERIC/libf/phystd/thermcell_height.F90 (revision 2177) @@ -2,6 +2,6 @@ ! ! -SUBROUTINE thermcell_height(ngrid,nlay,lmin,linter,lmix,lmax,zw2, & - zlev,zmin,zmix,zmax,wmax,f_star) +SUBROUTINE thermcell_height(ngrid,nlay,lmin,lmix,lmax,zlev, & + zmin,zmix,zmax,zw2,wmax,f_star) @@ -20,29 +20,30 @@ ! ------- - INTEGER ngrid, nlay - INTEGER lmin(ngrid) + INTEGER, INTENT(in) :: ngrid + INTEGER, INTENT(in) :: nlay + INTEGER, INTENT(in) :: lmin(ngrid) - REAL zlev(ngrid,nlay+1) - REAL f_star(ngrid,nlay+1) + REAL, INTENT(in) :: zlev(ngrid,nlay+1) + REAL, INTENT(in) :: f_star(ngrid,nlay+1) ! Outputs: ! -------- - INTEGER lmix(ngrid) - INTEGER lmax(ngrid) + INTEGER, INTENT(out) :: lmix(ngrid) + INTEGER, INTENT(out) :: lmax(ngrid) - REAL linter(ngrid) - REAL zmin(ngrid) ! Altitude of the plume bottom - REAL zmax(ngrid) ! Altitude of the plume top - REAL zmix(ngrid) ! Altitude of maximal vertical speed - REAL wmax(ngrid) ! Maximal vertical speed - REAL zw2(ngrid,nlay+1) ! Square vertical speed (becomes its square root) + REAL, INTENT(out) :: zmin(ngrid) ! Altitude of the plume bottom + REAL, INTENT(out) :: zmix(ngrid) ! Altitude of maximal vertical speed + REAL, INTENT(out) :: zmax(ngrid) ! Altitude of the plume top + REAL, INTENT(out) :: wmax(ngrid) ! Maximal vertical speed + + REAL, INTENT(inout) :: zw2(ngrid,nlay+1) ! Square vertical speed (becomes its square root) ! Local: ! ------ - INTEGER ig, l + INTEGER ig, k - REAL zlevinter(ngrid) + REAL linter(ngrid) ! Level (continuous) of maximal vertical speed !=============================================================================== @@ -50,17 +51,13 @@ !=============================================================================== - DO ig=1,ngrid - lmax(ig) = lmin(ig) - lmix(ig) = lmin(ig) - ENDDO + linter(:) = lmin(:) + lmix(:) = lmin(:) + lmax(:) = lmin(:) - DO ig=1,ngrid - wmax(ig) = 0. - ENDDO + wmax(:) = 0. - DO ig=1,ngrid - zmax(ig) = 0. - zlevinter(ig) = zlev(ig,1) - ENDDO + zmin(:) = 0. + zmix(:) = 0. + zmax(:) = 0. !=============================================================================== @@ -72,8 +69,6 @@ !------------------------------------------------------------------------------- - DO l=1,nlay - DO ig=1,ngrid - zmin(ig) = zlev(ig,lmin(ig)) - ENDDO + DO ig=1,ngrid + zmin(ig) = zlev(ig,lmin(ig)) ENDDO @@ -83,7 +78,7 @@ DO ig=1,ngrid - DO l=nlay,lmin(ig)+1,-1 - IF ((zw2(ig,l) <= 0.).or.(f_star(ig,l) <= 0.)) THEN - lmax(ig) = l - 1 + DO k=nlay,lmin(ig)+1,-1 + IF ((zw2(ig,k) <= 0.).or.(f_star(ig,k) <= 0.)) THEN + lmax(ig) = k - 1 ENDIF ENDDO @@ -96,5 +91,5 @@ IF (zw2(ig,nlay) > 0.) THEN print *, 'WARNING: a thermal plume reaches the highest layer!' - print *, 'ig,l', ig, nlay + print *, 'ig,k', ig, nlay print *, 'zw2', zw2(ig,nlay) lmax(ig) = nlay @@ -103,24 +98,20 @@ !------------------------------------------------------------------------------- -! Calcul de zmax continu via le linter +! Calcul de zmax continu via linter !------------------------------------------------------------------------------- !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -! AB : lmin=lmax means the plume is not active and then zw2=0 at each level. -! Otherwise we have lmax>lmin. +! AB: lmin=lmax means the plume is not active and then zw2=0 at each level. +! Otherwise we have lmax>lmin. !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ DO ig=1,ngrid - l = lmax(ig) - IF (l == lmin(ig)) THEN - linter(ig) = l + k = lmax(ig) + IF (k == lmin(ig)) THEN + linter(ig) = k ELSE - linter(ig) = l - zw2(ig,l) / (zw2(ig,l+1) - zw2(ig,l)) + linter(ig) = k - zw2(ig,k) / (zw2(ig,k+1) - zw2(ig,k)) ENDIF - ENDDO - - DO ig=1,ngrid - zlevinter(ig) = zlev(ig,lmax(ig)) + (linter(ig) - lmax(ig)) & - & * (zlev(ig,lmax(ig)+1) - zlev(ig,lmax(ig))) - zmax(ig) = max(zmax(ig), zlevinter(ig)) + zmax(ig) = zlev(ig,lmax(ig)) + (linter(ig) - lmax(ig)) & + & * (zlev(ig,lmax(ig)+1) - zlev(ig,lmax(ig))) ENDDO @@ -133,25 +124,25 @@ !------------------------------------------------------------------------------- - DO l=1,nlay + DO k=1,nlay DO ig=1,ngrid - IF((l <= lmax(ig)).and.(l > lmin(ig))) THEN - IF (zw2(ig,l).lt.0.) THEN + IF((k <= lmax(ig)).and.(k > lmin(ig))) THEN + IF (zw2(ig,k) < 0.) THEN print *, 'ERROR: zw2 has negative value(s)!' - print *, 'ig,l', ig, l - print *, 'zw2', zw2(ig,l) + print *, 'ig,k', ig, k + print *, 'zw2', zw2(ig,k) CALL abort ENDIF !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -! AB : WARNING zw2 becomes its square root! +! AB: WARNING zw2 becomes its square root! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - zw2(ig,l) = sqrt(zw2(ig,l)) + zw2(ig,k) = sqrt(zw2(ig,k)) - IF (zw2(ig,l) > wmax(ig)) THEN - wmax(ig) = zw2(ig,l) - lmix(ig) = l + IF (zw2(ig,k) > wmax(ig)) THEN + wmax(ig) = zw2(ig,k) + lmix(ig) = k ENDIF ELSE - zw2(ig,l) = 0. + zw2(ig,k) = 0. ENDIF ENDDO @@ -186,10 +177,6 @@ ENDDO -!=============================================================================== -! Check consistence between zmax and zmix -!=============================================================================== - !------------------------------------------------------------------------------- -! +! Check consistency between zmax and zmix !------------------------------------------------------------------------------- @@ -199,17 +186,10 @@ print *, 'zmax,zmix_old,zmix_new', zmax(ig), zmix(ig), 0.9 * zmax(ig) zmix(ig) = 0.9 * zmax(ig) + DO k=1,nlay + IF ((zmix(ig) >= zlev(ig,k)).and.(zmix(ig) < zlev(ig,k+1))) THEN + lmix(ig) = k + ENDIF + ENDDO ENDIF - ENDDO - -!------------------------------------------------------------------------------- -! Calcul du nouveau lmix correspondant -!------------------------------------------------------------------------------- - - DO ig=1,ngrid - DO l=1,nlay - IF ((zmix(ig) >= zlev(ig,l)).and.(zmix(ig) < zlev(ig,l+1))) THEN - lmix(ig) = l - ENDIF - ENDDO ENDDO Index: /trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90 (revision 2176) +++ /trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90 (revision 2177) @@ -42,7 +42,5 @@ ! New detr and entre formulae (no longer alimentation) ! lmin can be greater than 1 -! Mix every tracer (EN COURS) -! Old version of thermcell_dq is removed -! Alternative version thermcell_dv2 is removed +! Mix every tracer ! !=============================================================================== @@ -61,30 +59,36 @@ ! ------- - INTEGER ngrid, nlay, nq - - REAL ptimestep - REAL pplay(ngrid,nlay) ! Layer pressure - REAL pplev(ngrid,nlay+1) ! Level pressure - REAL pphi(ngrid,nlay) ! Geopotential - - REAL pu(ngrid,nlay) ! Zonal wind - REAL pv(ngrid,nlay) ! Meridional wind - REAL pt(ngrid,nlay) ! Temperature - REAL pq(ngrid,nlay,nq) ! Tracers mass mixing ratio - - LOGICAL firstcall + INTEGER, INTENT(in) :: ngrid + INTEGER, INTENT(in) :: nlay + INTEGER, INTENT(in) :: nq + + REAL, INTENT(in) :: ptimestep + REAL, INTENT(in) :: pplay(ngrid,nlay) ! Layer pressure + REAL, INTENT(in) :: pplev(ngrid,nlay+1) ! Level pressure + REAL, INTENT(in) :: pphi(ngrid,nlay) ! Geopotential + + REAL, INTENT(in) :: pu(ngrid,nlay) ! Zonal wind + REAL, INTENT(in) :: pv(ngrid,nlay) ! Meridional wind + REAL, INTENT(in) :: pt(ngrid,nlay) ! Temperature + REAL, INTENT(in) :: pq(ngrid,nlay,nq) ! Tracers mass mixing ratio + + LOGICAL, INTENT(in) :: firstcall ! Outputs: ! -------- - REAL pduadj(ngrid,nlay) ! u convective variations - REAL pdvadj(ngrid,nlay) ! v convective variations - REAL pdtadj(ngrid,nlay) ! t convective variations - REAL pdqadj(ngrid,nlay,nq) ! q convective variations - - REAL f0(ngrid) ! mass flux norm (after possible time relaxation) - REAL fm0(ngrid,nlay+1) ! mass flux (after possible time relaxation) - REAL entr0(ngrid,nlay) ! entrainment (after possible time relaxation) - REAL detr0(ngrid,nlay) ! detrainment (after possible time relaxation) + INTEGER, INTENT(out) :: lmax(ngrid) ! Highest layer reached by the plume + INTEGER, INTENT(out) :: lmix(ngrid) ! Layer in which plume vertical speed is maximal + INTEGER, INTENT(out) :: lmin(ngrid) ! First unstable layer + + REAL, INTENT(out) :: pduadj(ngrid,nlay) ! u convective variations + REAL, INTENT(out) :: pdvadj(ngrid,nlay) ! v convective variations + REAL, INTENT(out) :: pdtadj(ngrid,nlay) ! t convective variations + REAL, INTENT(out) :: pdqadj(ngrid,nlay,nq) ! q convective variations + + REAL, INTENT(inout) :: f0(ngrid) ! mass flux norm (after possible time relaxation) + REAL, INTENT(inout) :: fm0(ngrid,nlay+1) ! mass flux (after possible time relaxation) + REAL, INTENT(inout) :: entr0(ngrid,nlay) ! entrainment (after possible time relaxation) + REAL, INTENT(inout) :: detr0(ngrid,nlay) ! detrainment (after possible time relaxation) ! Local: @@ -92,55 +96,50 @@ INTEGER ig, k, l, iq - INTEGER lmax(ngrid) ! Highest layer reached by the plume - INTEGER lmix(ngrid) ! Layer in which plume vertical speed is maximal - INTEGER lmin(ngrid) ! First unstable layer - - REAL zmix(ngrid) ! Altitude of maximal vertical speed - REAL zmax(ngrid) ! Maximal altitudes where plumes are active - REAL zmin(ngrid) ! Minimal altitudes where plumes are active - - REAL zlay(ngrid,nlay) ! Layers altitudes - REAL zlev(ngrid,nlay+1) ! Levels altitudes - REAL rho(ngrid,nlay) ! Layers densities - REAL rhobarz(ngrid,nlay) ! Levels densities - REAL masse(ngrid,nlay) ! Layers masses - REAL zpopsk(ngrid,nlay) ! Exner function - - REAL zu(ngrid,nlay) ! u environment - REAL zv(ngrid,nlay) ! v environment - REAL zt(ngrid,nlay) ! TR environment - REAL zqt(ngrid,nlay) ! qt environment - REAL zql(ngrid,nlay) ! ql environment - REAL zhl(ngrid,nlay) ! TP environment - REAL ztv(ngrid,nlay) ! TRPV environment - REAL zqs(ngrid,nlay) ! qsat environment - - REAL zua(ngrid,nlay) ! u plume - REAL zva(ngrid,nlay) ! v plume - REAL zta(ngrid,nlay) ! TR plume - REAL zqla(ngrid,nlay) ! qv plume - REAL zqta(ngrid,nlay) ! qt plume - REAL zhla(ngrid,nlay) ! TP plume - REAL ztva(ngrid,nlay) ! TRPV plume - REAL zqsa(ngrid,nlay) ! qsat plume - - REAL zqa(ngrid,nlay,nq) ! q plume (ql=0, qv=qt) - - REAL f_star(ngrid,nlay+1) ! Normalized mass flux - REAL entr_star(ngrid,nlay) ! Normalized entrainment (E* = e* dz) - REAL detr_star(ngrid,nlay) ! Normalized detrainment (D* = d* dz) - - REAL fm(ngrid,nlay+1) ! Mass flux - REAL entr(ngrid,nlay) ! Entrainment (E = e dz) - REAL detr(ngrid,nlay) ! Detrainment (D = d dz) - - REAL f(ngrid) ! Mass flux norm - REAL lambda ! Time relaxation coefficent - REAL fraca(ngrid,nlay+1) ! Updraft fraction - REAL linter(ngrid) ! Level (continuous) of maximal vertical speed - REAL wmax(ngrid) ! Maximal vertical speed - REAL zw2(ngrid,nlay+1) ! Plume vertical speed - REAL zdthladj(ngrid,nlay) ! Potential temperature variations - REAL dummy(ngrid,nlay) ! Dummy argument for thermcell_dq() + + REAL zmix(ngrid) ! Altitude of maximal vertical speed + REAL zmax(ngrid) ! Maximal altitudes where plumes are active + REAL zmin(ngrid) ! Minimal altitudes where plumes are active + + REAL zlay(ngrid,nlay) ! Layers altitudes + REAL zlev(ngrid,nlay+1) ! Levels altitudes + REAL rho(ngrid,nlay) ! Layers densities + REAL rhobarz(ngrid,nlay) ! Levels densities + REAL masse(ngrid,nlay) ! Layers masses + REAL zpopsk(ngrid,nlay) ! Exner function + + REAL zu(ngrid,nlay) ! u environment + REAL zv(ngrid,nlay) ! v environment + REAL zt(ngrid,nlay) ! TR environment + REAL zqt(ngrid,nlay) ! qt environment + REAL zql(ngrid,nlay) ! ql environment + REAL zhl(ngrid,nlay) ! TP environment + REAL ztv(ngrid,nlay) ! TRPV environment + REAL zqs(ngrid,nlay) ! qsat environment + + REAL zua(ngrid,nlay) ! u plume + REAL zva(ngrid,nlay) ! v plume + REAL zqla(ngrid,nlay) ! qv plume + REAL zqta(ngrid,nlay) ! qt plume + REAL zhla(ngrid,nlay) ! TP plume + REAL ztva(ngrid,nlay) ! TRPV plume + REAL zqsa(ngrid,nlay) ! qsat plume + + REAL zqa(ngrid,nlay,nq) ! q plume (ql=0, qv=qt) + + REAL f_star(ngrid,nlay+1) ! Normalized mass flux + REAL entr_star(ngrid,nlay) ! Normalized entrainment (E* = e* dz) + REAL detr_star(ngrid,nlay) ! Normalized detrainment (D* = d* dz) + + REAL fm(ngrid,nlay+1) ! Mass flux + REAL entr(ngrid,nlay) ! Entrainment (E = e dz) + REAL detr(ngrid,nlay) ! Detrainment (D = d dz) + + REAL f(ngrid) ! Mass flux norm + REAL lambda ! Time relaxation coefficent + REAL fraca(ngrid,nlay+1) ! Updraft fraction + REAL wmax(ngrid) ! Maximal vertical speed + REAL zw2(ngrid,nlay+1) ! Plume vertical speed + REAL zdthladj(ngrid,nlay) ! Potential temperature variations + REAL dummy(ngrid,nlay) ! Dummy argument for thermcell_dq() !=============================================================================== @@ -154,17 +153,8 @@ ENDIF - f_star(:,:) = 0. - entr_star(:,:) = 0. - detr_star(:,:) = 0. - - f(:) = 0. - - fm(:,:) = 0. - entr(:,:) = 0. - detr(:,:) = 0. - - lmax(:) = 1 - lmix(:) = 1 - lmin(:) = 1 + DO ig=1,ngrid +! AB: Minimal f0 value is set to 0. (instead of 1.e-2 in Earth version) + f0(ig) = MAX(f0(ig), 0.) + ENDDO pduadj(:,:) = 0.0 @@ -173,10 +163,5 @@ pdqadj(:,:,:) = 0.0 - DO ig=1,ngrid -! AB: Careful: Hard-coded value from Earth version! -! f0(ig) = max(f0(ig), 1.e-2) -! AB: No pescribed minimal value for f0 - f0(ig) = max(f0(ig), 0.) - ENDDO + zdthladj(:,:) = 0.0 !=============================================================================== @@ -212,11 +197,8 @@ rho(:,:) = pplay(:,:) / (zpopsk(:,:) * RD * ztv(:,:)) + rhobarz(:,1) = rho(:,1) IF (prt_level.ge.10) THEN print *, 'WARNING: density in the first layer is equal to density at the first level!' - print *, 'rhobarz(:,1)', rhobarz(:,1) - print *, 'rho(:,1) ', rho(:,1) - ENDIF - - rhobarz(:,1) = rho(:,1) + ENDIF DO l=2,nlay @@ -244,32 +226,32 @@ ! ! --------------------------- -! _ -! ----- F_lmax+1=0 ------zmax \ -! lmax | -! ------F_lmax>0------------- | -! | -! --------------------------- | -! | -! --------------------------- | -! | -! ------------------wmax,zmix | -! lmix | -! --------------------------- | -! | -! --------------------------- | -! | E, D -! --------------------------- | -! | +! _ +! ----- F_lmax+1=0 ------zmax \ +! lmax | +! ------F_lmax>0------------- | +! | +! --------------------------- | +! | +! --------------------------- | +! | +! ------------------wmax,zmix | +! lmix | +! --------------------------- | +! | +! --------------------------- | +! | E, D +! --------------------------- | +! | ! --------------------------- rhobarz, f_star, fm, fm0, zw2, fraca -! zt, zu, zv, zo, rho | -! --------------------------- | -! | -! --------------------------- | -! | -! --------------------------- | -! | -! ------F_lmin+1>0----------- | -! lmin | -! ----- F_lmin=0 ------------ _/ +! zt, zu, zv, zo, rho | +! --------------------------- | +! | +! --------------------------- | +! | +! --------------------------- | +! | +! ------F_lmin+1>0----------- | +! lmin | +! ----- F_lmin=0 ------------ _/ ! ! --------------------------- @@ -309,9 +291,9 @@ !------------------------------------------------------------------------------- - CALL thermcell_plume(ngrid,nlay,nq,ptimestep,ztv, & - & zhl,zqt,zql,rhobarz,zlev,pplev,pphi,zpopsk, & + CALL thermcell_plume(ngrid,nlay,nq,ptimestep, & + & ztv,zhl,zqt,zql,zlev,pplev,zpopsk, & & detr_star,entr_star,f_star, & - & ztva,zhla,zqla,zqta,zta,zqsa, & - & zw2,lmix,lmin) + & ztva,zhla,zqta,zqla,zqsa, & + & zw2,lmin) !------------------------------------------------------------------------------- @@ -320,6 +302,6 @@ ! AB: Careful, zw2 became its square root in thermcell_height! - CALL thermcell_height(ngrid,nlay,lmin,linter,lmix,lmax,zw2, & - & zlev,zmin,zmix,zmax,wmax,f_star) + CALL thermcell_height(ngrid,nlay,lmin,lmix,lmax,zlev, & + & zmin,zmix,zmax,zw2,wmax,f_star) !=============================================================================== @@ -341,7 +323,5 @@ ENDIF -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -! Test valable seulement en 1D mais pas genant -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +! FH: Test valable seulement en 1D mais pas genant IF (.not. (f0(1).ge.0.) ) THEN print *, 'ERROR: mass flux norm is not positive!' @@ -403,5 +383,5 @@ CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse, & - & zhl,zdthladj,dummy,lmin) + & zhl,zdthladj,dummy,lmin,lmax) DO l=1,nlay @@ -417,5 +397,5 @@ DO iq=1,nq CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse, & - & pq(:,:,iq),pdqadj(:,:,iq),zqa(:,:,iq),lmin) + & pq(:,:,iq),pdqadj(:,:,iq),zqa(:,:,iq),lmin,lmax) ENDDO @@ -429,7 +409,7 @@ ELSE CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse, & - & zu,pduadj,zua,lmin) + & zu,pduadj,zua,lmin,lmax) CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse, & - & zv,pdvadj,zva,lmin) + & zv,pdvadj,zva,lmin,lmax) ENDIF