! ! ! SUBROUTINE thermcell_dq(ngrid,nlay,impl,ptimestep,fm,entr,detr,masse, & q,dq,qa,lmin,lev_out) !============================================================================== ! Calcul du transport verticale dans la couche limite en presence ! de "thermiques" explicitement representes ! calcul du dq/dt une fois qu'on connait les ascendances ! ! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr) ! Introduction of an implicit computation of vertical advection in ! the environment of thermal plumes in thermcell_dq ! impl = 0 : explicit, 1 : implicit, -1 : old version ! ! Modif 2018/11/05 (AB alexandre.boissinot@lmd.jussieu.fr) ! ! !============================================================================== USE print_control_mod, ONLY: prt_level IMPLICIT NONE !============================================================================== ! Declaration !============================================================================== ! inputs: ! ------- INTEGER ngrid, nlay INTEGER impl INTEGER lmin(ngrid) INTEGER lev_out ! niveau pour les print REAL ptimestep REAL masse(ngrid,nlay) REAL fm(ngrid,nlay+1) REAL entr(ngrid,nlay) REAL detr(ngrid,nlay) REAL q(ngrid,nlay) REAL qa(ngrid,nlay) ! outputs: ! -------- REAL dq(ngrid,nlay) ! local: ! ------ INTEGER ig, l INTEGER niter, iter REAL cfl REAL qold(ngrid,nlay) REAL fqa(ngrid,nlay+1) REAL wqd(ngrid,nlay+1) REAL zzm CHARACTER (LEN=20) :: modname CHARACTER (LEN=80) :: abort_message !============================================================================== ! Initialization !============================================================================== qold = q modname = 'thermcell_dq' IF (impl.lt.0) THEN print *, 'OLD EXPLICIT SCHEME IS USED' CALL thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr, & & masse,q,dq,qa,lev_out) RETURN ENDIF !------------------------------------------------------------------------------ ! CFL criterion computation for advection in downdraft !------------------------------------------------------------------------------ cfl = 0. DO l=1,nlay DO ig=1,ngrid zzm = masse(ig,l) / ptimestep cfl = max(cfl, fm(ig,l) / zzm) IF (entr(ig,l).gt.zzm) THEN print *, 'ERROR: entrainment is greater than the layer mass!' print *, 'ig,l,entr', ig, l, entr(ig,l) print *, '-------------------------------' print *, 'entr*dt,mass', entr(ig,l)*ptimestep, masse(ig,l) print *, '-------------------------------' print *, 'fm+', fm(ig,l+1) print *, 'entr,detr', entr(ig,l), detr(ig,l) print *, 'fm ', fm(ig,l) print *, 'entr,detr', entr(ig,l-1), detr(ig,l-1) print *, 'fm-', fm(ig,l-1) abort_message = 'entr dt > m, 1st' CALL abort_physic(modname,abort_message,1) ENDIF ENDDO ENDDO !------------------------------------------------------------------------------ ! Computation of tracer concentrations in the ascending plume !------------------------------------------------------------------------------ DO ig=1,ngrid DO l=1,lmin(ig) qa(ig,l) = q(ig,l) ENDDO ENDDO DO ig=1,ngrid DO l=lmin(ig)+1,nlay if ((fm(ig,l+1)+detr(ig,l))*ptimestep.gt.1.e-5*masse(ig,l)) then qa(ig,l) = (fm(ig,l) * qa(ig,l-1) + entr(ig,l) * q(ig,l)) & & / (fm(ig,l+1) + detr(ig,l)) else qa(ig,l) = q(ig,l) endif ! IF (qa(ig,l).lt.0.) THEN ! print *, 'WARNING: qa is negative!' ! print *, 'qa', qa(ig,l) ! ENDIF ! IF (q(ig,l).lt.0.) THEN ! print *, 'WARNING: q is negative!' ! print *, 'q', q(ig,l) ! ENDIF ENDDO ENDDO !------------------------------------------------------------------------------ ! Plume vertical flux of q !------------------------------------------------------------------------------ DO l=2,nlay-1 fqa(:,l) = fm(:,l) * qa(:,l-1) ENDDO fqa(:,1) = 0. fqa(:,nlay) = 0. !------------------------------------------------------------------------------ ! Trace species evolution !------------------------------------------------------------------------------ IF (impl==0) THEN do l=1,nlay-1 q(:,l) = q(:,l) + (fqa(:,l) - fqa(:,l+1) - fm(:,l) * q(:,l) & & + fm(:,l+1) * q(:,l+1)) * ptimestep / masse(:,l) enddo ELSE do l=nlay-1,1,-1 q(:,l) = ( q(:,l) + ptimestep / masse(:,l) & & * ( fqa(:,l) - fqa(:,l+1) + fm(:,l+1) * q(:,l+1) ) ) & & / ( 1. + fm(:,l) * ptimestep / masse(:,l) ) ENDDO ENDIF !============================================================================== ! Tendencies !============================================================================== DO l=1,nlay DO ig=1,ngrid dq(ig,l) = (q(ig,l) - qold(ig,l)) / ptimestep q(ig,l) = qold(ig,l) ENDDO ENDDO RETURN END !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Subroutine thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr,masse, & q,dq,qa,lev_out) !============================================================================== ! ! Calcul du transport verticale dans la couche limite en presence ! de "thermiques" explicitement representes ! calcul du dq/dt une fois qu'on connait les ascendances ! !============================================================================== USE print_control_mod, ONLY: prt_level implicit none !============================================================================== ! Declaration !============================================================================== integer ngrid,nlay,impl real ptimestep real masse(ngrid,nlay),fm(ngrid,nlay+1) real entr(ngrid,nlay) real q(ngrid,nlay) real dq(ngrid,nlay) integer lev_out ! niveau pour les print real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1) real zzm integer ig,l real cfl real qold(ngrid,nlay) real ztimestep integer niter,iter CHARACTER (LEN=20) :: modname='thermcell_dq' CHARACTER (LEN=80) :: abort_message ! Calcul du critere CFL pour l'advection dans la subsidence cfl = 0. do l=1,nlay do ig=1,ngrid zzm=masse(ig,l)/ptimestep cfl=max(cfl,fm(ig,l)/zzm) if (entr(ig,l).gt.zzm) then print*,'entr*dt>m,2',l,entr(ig,l)*ptimestep,masse(ig,l) abort_message = 'entr dt > m, 2nd' CALL abort_physic (modname,abort_message,1) endif enddo enddo !IM 090508 print*,'CFL CFL CFL CFL ',cfl #undef CFL #ifdef CFL ! On subdivise le calcul en niter pas de temps. niter=int(cfl)+1 #else niter=1 #endif ztimestep=ptimestep/niter qold=q do iter=1,niter if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0' ! calcul du detrainement do l=1,nlay do ig=1,ngrid detr(ig,l)=fm(ig,l)-fm(ig,l+1)+entr(ig,l) ! print*,'Q2 DQ ',detr(ig,l),fm(ig,l),entr(ig,l) !test if (detr(ig,l).lt.0.) then entr(ig,l)=entr(ig,l)-detr(ig,l) detr(ig,l)=0. ! print*,'detr2<0!!!','ig=',ig,'l=',l,'f=',fm(ig,l), ! s 'f+1=',fm(ig,l+1),'e=',entr(ig,l),'d=',detr(ig,l) endif ! if (fm(ig,l+1).lt.0.) then ! print*,'fm2<0!!!' ! endif ! if (entr(ig,l).lt.0.) then ! print*,'entr2<0!!!' ! endif enddo enddo ! calcul de la valeur dans les ascendances do ig=1,ngrid qa(ig,1)=q(ig,1) enddo do l=2,nlay do ig=1,ngrid if ((fm(ig,l+1)+detr(ig,l))*ztimestep.gt.1.e-5*masse(ig,l)) then qa(ig,l) = (fm(ig,l)*qa(ig,l-1)+entr(ig,l)*q(ig,l)) & & / (fm(ig,l+1)+detr(ig,l)) else qa(ig,l)=q(ig,l) endif if (qa(ig,l).lt.0.) then ! print*,'qa<0!!!' endif if (q(ig,l).lt.0.) then ! print*,'q<0!!!' endif enddo enddo ! Calcul du flux subsident do l=2,nlay do ig=1,ngrid #undef centre #ifdef centre wqd(ig,l)=fm(ig,l)*0.5*(q(ig,l-1)+q(ig,l)) #else #define plusqueun #ifdef plusqueun ! Schema avec advection sur plus qu'une maille. zzm=masse(ig,l)/ztimestep if (fm(ig,l)>zzm) then wqd(ig,l)=zzm*q(ig,l)+(fm(ig,l)-zzm)*q(ig,l+1) else wqd(ig,l)=fm(ig,l)*q(ig,l) endif #else wqd(ig,l)=fm(ig,l)*q(ig,l) #endif #endif if (wqd(ig,l).lt.0.) then ! print*,'wqd<0!!!' endif enddo enddo do ig=1,ngrid wqd(ig,1)=0. wqd(ig,nlay+1)=0. enddo ! Calcul des tendances do l=1,nlay do ig=1,ngrid q(ig,l)=q(ig,l)+(detr(ig,l)*qa(ig,l)-entr(ig,l)*q(ig,l) & & -wqd(ig,l)+wqd(ig,l+1)) & & *ztimestep/masse(ig,l) ! if (dq(ig,l).lt.0.) then ! print*,'dq<0!!!' ! endif enddo enddo enddo ! Calcul des tendances do l=1,nlay do ig=1,ngrid dq(ig,l) = (q(ig,l) - qold(ig,l)) / ptimestep q(ig,l) = qold(ig,l) enddo enddo return end