MODULE lmdz_thermcell_dq CONTAINS SUBROUTINE thermcell_dq(ngrid, nlay, impl, ptimestep, fm, entr, & masse, q, dq, qa, lev_out) USE lmdz_print_control, ONLY: prt_level USE lmdz_abort_physic, ONLY: abort_physic IMPLICIT NONE !======================================================================= ! 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 !======================================================================= ! arguments INTEGER, INTENT(IN) :: ngrid, nlay, impl REAL, INTENT(IN) :: ptimestep REAL, INTENT(IN), DIMENSION(ngrid, nlay) :: masse REAL, INTENT(INOUT), DIMENSION(ngrid, nlay) :: entr, q REAL, INTENT(IN), DIMENSION(ngrid, nlay + 1) :: fm REAL, INTENT(OUT), DIMENSION(ngrid, nlay) :: dq, qa INTEGER, INTENT(IN) :: lev_out ! niveau pour les print ! Local REAL, DIMENSION(ngrid, nlay) :: detr, qold REAL, DIMENSION(ngrid, nlay + 1) :: wqd, fqa REAL zzm INTEGER ig, k REAL cfl INTEGER niter, iter CHARACTER (LEN = 20) :: modname = 'thermcell_dq' CHARACTER (LEN = 80) :: abort_message ! Old explicite scheme IF (impl<=-1) THEN CALL thermcell_dq_o(ngrid, nlay, impl, ptimestep, fm, entr, & masse, q, dq, qa, lev_out) else ! Calcul du critere CFL pour l'advection dans la subsidence cfl = 0. DO k = 1, nlay DO ig = 1, ngrid zzm = masse(ig, k) / ptimestep cfl = max(cfl, fm(ig, k) / zzm) IF (entr(ig, k)>zzm) THEN PRINT*, 'entr*dt>m,1', k, entr(ig, k) * ptimestep, masse(ig, k) abort_message = 'entr dt > m, 1st' CALL abort_physic (modname, abort_message, 1) endif enddo enddo qold = q IF (prt_level>=1) PRINT*, 'Q2 THERMCEL_DQ 0' ! calcul du detrainement DO k = 1, nlay DO ig = 1, ngrid detr(ig, k) = fm(ig, k) - fm(ig, k + 1) + entr(ig, k) ! PRINT*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k) !test IF (detr(ig, k)<0.) THEN entr(ig, k) = entr(ig, k) - detr(ig, k) detr(ig, k) = 0. ! PRINT*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k), ! s 'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k) endif IF (fm(ig, k + 1)<0.) THEN ! PRINT*,'fm2<0!!!' endif IF (entr(ig, k)<0.) THEN ! PRINT*,'entr2<0!!!' endif enddo enddo ! Computation of tracer concentrations in the ascending plume DO ig = 1, ngrid qa(ig, 1) = q(ig, 1) enddo DO k = 2, nlay DO ig = 1, ngrid IF ((fm(ig, k + 1) + detr(ig, k)) * ptimestep> & 1.e-5 * masse(ig, k)) THEN qa(ig, k) = (fm(ig, k) * qa(ig, k - 1) + entr(ig, k) * q(ig, k)) & / (fm(ig, k + 1) + detr(ig, k)) else qa(ig, k) = q(ig, k) endif IF (qa(ig, k)<0.) THEN ! PRINT*,'qa<0!!!' endif IF (q(ig, k)<0.) THEN ! PRINT*,'q<0!!!' endif enddo enddo ! Plume vertical flux DO k = 2, nlay - 1 fqa(:, k) = fm(:, k) * qa(:, k - 1) enddo fqa(:, 1) = 0. ; fqa(:, nlay) = 0. ! Trace species evolution IF (impl==0) THEN DO k = 1, nlay - 1 q(:, k) = q(:, k) + (fqa(:, k) - fqa(:, k + 1) - fm(:, k) * q(:, k) + fm(:, k + 1) * q(:, k + 1)) & * ptimestep / masse(:, k) enddo else DO k = nlay - 1, 1, -1 ! FH debut de modif : le calcul ci dessous modifiait numériquement ! la concentration quand le flux de masse etait nul car on divisait ! puis multipliait par masse/ptimestep. ! q(:,k)=(masse(:,k)*q(:,k)/ptimestep+fqa(:,k)-fqa(:,k+1)+fm(:,k+1)*q(:,k+1)) & ! & /(fm(:,k)+masse(:,k)/ptimestep) q(:, k) = (q(:, k) + ptimestep / masse(:, k) * (fqa(:, k) - fqa(:, k + 1) + fm(:, k + 1) * q(:, k + 1))) & / (1. + fm(:, k) * ptimestep / masse(:, k)) ! FH fin de modif. enddo endif ! Tendencies DO k = 1, nlay DO ig = 1, ngrid dq(ig, k) = (q(ig, k) - qold(ig, k)) / ptimestep q(ig, k) = qold(ig, k) enddo enddo END IF ! impl=-1 RETURN end !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Obsolete version kept for convergence with Cmip5 NPv3.1 simulations !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE thermcell_dq_o(ngrid, nlay, impl, ptimestep, fm, entr, & masse, q, dq, qa, lev_out) USE lmdz_print_control, ONLY: prt_level USE lmdz_abort_physic, ONLY: abort_physic IMPLICIT NONE !======================================================================= ! 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 !======================================================================= 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, k 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 k = 1, nlay DO ig = 1, ngrid zzm = masse(ig, k) / ptimestep cfl = max(cfl, fm(ig, k) / zzm) IF (entr(ig, k)>zzm) THEN PRINT*, 'entr*dt>m,2', k, entr(ig, k) * ptimestep, masse(ig, k) 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>=1) PRINT*, 'Q2 THERMCEL_DQ 0' ! calcul du detrainement DO k = 1, nlay DO ig = 1, ngrid detr(ig, k) = fm(ig, k) - fm(ig, k + 1) + entr(ig, k) ! PRINT*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k) !test IF (detr(ig, k)<0.) THEN entr(ig, k) = entr(ig, k) - detr(ig, k) detr(ig, k) = 0. ! PRINT*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k), ! s 'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k) endif IF (fm(ig, k + 1)<0.) THEN ! PRINT*,'fm2<0!!!' endif IF (entr(ig, k)<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 k = 2, nlay DO ig = 1, ngrid IF ((fm(ig, k + 1) + detr(ig, k)) * ztimestep> & 1.e-5 * masse(ig, k)) THEN qa(ig, k) = (fm(ig, k) * qa(ig, k - 1) + entr(ig, k) * q(ig, k)) & / (fm(ig, k + 1) + detr(ig, k)) else qa(ig, k) = q(ig, k) endif IF (qa(ig, k)<0.) THEN ! PRINT*,'qa<0!!!' endif IF (q(ig, k)<0.) THEN ! PRINT*,'q<0!!!' endif enddo enddo ! Calcul du flux subsident DO k = 2, nlay DO ig = 1, ngrid #undef centre #ifdef centre wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k)) #else #define CFL_plus_grand_que_un #ifdef CFL_plus_grand_que_un ! Schema avec advection sur plus qu'une maille. zzm = masse(ig, k) / ztimestep IF (fm(ig, k)>zzm) THEN wqd(ig, k) = zzm * q(ig, k) + (fm(ig, k) - zzm) * q(ig, k + 1) else wqd(ig, k) = fm(ig, k) * q(ig, k) endif #else wqd(ig,k)=fm(ig,k)*q(ig,k) #endif #endif IF (wqd(ig, k)<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 k = 1, nlay DO ig = 1, ngrid q(ig, k) = q(ig, k) + (detr(ig, k) * qa(ig, k) - entr(ig, k) * q(ig, k) & - wqd(ig, k) + wqd(ig, k + 1)) & * ztimestep / masse(ig, k) ! if (dq(ig,k).lt.0.) THEN ! PRINT*,'dq<0!!!' ! endif enddo enddo END DO ! Calcul des tendances DO k = 1, nlay DO ig = 1, ngrid dq(ig, k) = (q(ig, k) - qold(ig, k)) / ptimestep q(ig, k) = qold(ig, k) enddo enddo RETURN END END MODULE lmdz_thermcell_dq