! $Id: thermcell.F90 2346 2015-08-21 15:13:46Z crio $ SUBROUTINE calcul_sec(ngrid, nlay, ptimestep, pplay, pplev, pphi, zlev, pu, & pv, pt, po, zmax, wmax, zw2, lmix & ! s ! ,pu_therm,pv_therm , r_aspect, l_mix, w2di, tho) USE dimphy IMPLICIT NONE ! ======================================================================= ! Calcul du transport verticale dans la couche limite en presence ! de "thermiques" explicitement representes ! Réécriture à partir d'un listing papier à Habas, le 14/02/00 ! le thermique est supposé homogène et dissipé par mélange avec ! son environnement. la longueur l_mix contrôle l'efficacité du ! mélange ! Le calcul du transport des différentes espèces se fait en prenant ! en compte: ! 1. un flux de masse montant ! 2. un flux de masse descendant ! 3. un entrainement ! 4. un detrainement ! ======================================================================= ! ----------------------------------------------------------------------- ! declarations: ! ------------- include "YOMCST.h" ! arguments: ! ---------- INTEGER ngrid, nlay, w2di REAL tho REAL ptimestep, l_mix, r_aspect REAL pt(ngrid, nlay), pdtadj(ngrid, nlay) REAL pu(ngrid, nlay), pduadj(ngrid, nlay) REAL pv(ngrid, nlay), pdvadj(ngrid, nlay) REAL po(ngrid, nlay), pdoadj(ngrid, nlay) REAL pplay(ngrid, nlay), pplev(ngrid, nlay+1) REAL pphi(ngrid, nlay) INTEGER idetr SAVE idetr DATA idetr/3/ !$OMP THREADPRIVATE(idetr) ! local: ! ------ INTEGER ig, k, l, lmaxa(klon), lmix(klon) REAL zsortie1d(klon) ! CR: on remplace lmax(klon,klev+1) INTEGER lmax(klon), lmin(klon), lentr(klon) REAL linter(klon) REAL zmix(klon), fracazmix(klon) ! RC REAL zmax(klon), zw, zw2(klon, klev+1), ztva(klon, klev) REAL zlev(klon, klev+1), zlay(klon, klev) REAL zh(klon, klev), zdhadj(klon, klev) REAL ztv(klon, klev) REAL zu(klon, klev), zv(klon, klev), zo(klon, klev) REAL wh(klon, klev+1) REAL wu(klon, klev+1), wv(klon, klev+1), wo(klon, klev+1) REAL zla(klon, klev+1) REAL zwa(klon, klev+1) REAL zld(klon, klev+1) ! real zwd(klon,klev+1) REAL zsortie(klon, klev) REAL zva(klon, klev) REAL zua(klon, klev) REAL zoa(klon, klev) REAL zha(klon, klev) REAL wa_moy(klon, klev+1) REAL fraca(klon, klev+1) REAL fracc(klon, klev+1) REAL zf, zf2 REAL thetath2(klon, klev), wth2(klon, klev) ! common/comtherm/thetath2,wth2 REAL count_time ! integer isplit,nsplit INTEGER isplit, nsplit, ialt PARAMETER (nsplit=10) DATA isplit/0/ SAVE isplit !$OMP THREADPRIVATE(isplit) LOGICAL sorties REAL rho(klon, klev), rhobarz(klon, klev+1), masse(klon, klev) REAL zpspsk(klon, klev) ! real wmax(klon,klev),wmaxa(klon) REAL wmax(klon), wmaxa(klon) REAL wa(klon, klev, klev+1) REAL wd(klon, klev+1) REAL larg_part(klon, klev, klev+1) REAL fracd(klon, klev+1) REAL xxx(klon, klev+1) REAL larg_cons(klon, klev+1) REAL larg_detr(klon, klev+1) REAL fm0(klon, klev+1), entr0(klon, klev), detr(klon, klev) REAL pu_therm(klon, klev), pv_therm(klon, klev) REAL fm(klon, klev+1), entr(klon, klev) REAL fmc(klon, klev+1) ! CR:nouvelles variables REAL f_star(klon, klev+1), entr_star(klon, klev) REAL entr_star_tot(klon), entr_star2(klon) REAL zalim(klon) INTEGER lalim(klon) REAL norme(klon) REAL f(klon), f0(klon) REAL zlevinter(klon) LOGICAL therm LOGICAL first DATA first/.FALSE./ SAVE first !$OMP THREADPRIVATE(first) ! RC CHARACTER *2 str2 CHARACTER *10 str10 CHARACTER (LEN=20) :: modname = 'calcul_sec' CHARACTER (LEN=80) :: abort_message ! LOGICAL vtest(klon),down EXTERNAL scopy INTEGER ncorrec SAVE ncorrec DATA ncorrec/0/ !$OMP THREADPRIVATE(ncorrec) ! ----------------------------------------------------------------------- ! initialisation: ! --------------- sorties = .TRUE. IF (ngrid/=klon) THEN PRINT * PRINT *, 'STOP dans convadj' PRINT *, 'ngrid =', ngrid PRINT *, 'klon =', klon END IF ! ----------------------------------------------------------------------- ! incrementation eventuelle de tendances precedentes: ! --------------------------------------------------- ! print*,'0 OK convect8' DO l = 1, nlay DO ig = 1, ngrid zpspsk(ig, l) = (pplay(ig,l)/pplev(ig,1))**rkappa zh(ig, l) = pt(ig, l)/zpspsk(ig, l) zu(ig, l) = pu(ig, l) zv(ig, l) = pv(ig, l) zo(ig, l) = po(ig, l) ztv(ig, l) = zh(ig, l)*(1.+0.61*zo(ig,l)) END DO END DO ! print*,'1 OK convect8' ! -------------------- ! + + + + + + + + + + + ! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz ! wh,wt,wo ... ! + + + + + + + + + + + zh,zu,zv,zo,rho ! -------------------- zlev(1) ! \\\\\\\\\\\\\\\\\\\\ ! ----------------------------------------------------------------------- ! Calcul des altitudes des couches ! ----------------------------------------------------------------------- DO l = 2, nlay DO ig = 1, ngrid zlev(ig, l) = 0.5*(pphi(ig,l)+pphi(ig,l-1))/rg END DO END DO DO ig = 1, ngrid zlev(ig, 1) = 0. zlev(ig, nlay+1) = (2.*pphi(ig,klev)-pphi(ig,klev-1))/rg END DO DO l = 1, nlay DO ig = 1, ngrid zlay(ig, l) = pphi(ig, l)/rg END DO END DO ! print*,'2 OK convect8' ! ----------------------------------------------------------------------- ! Calcul des densites ! ----------------------------------------------------------------------- DO l = 1, nlay DO ig = 1, ngrid rho(ig, l) = pplay(ig, l)/(zpspsk(ig,l)*rd*zh(ig,l)) END DO END DO DO l = 2, nlay DO ig = 1, ngrid rhobarz(ig, l) = 0.5*(rho(ig,l)+rho(ig,l-1)) END DO END DO DO k = 1, nlay DO l = 1, nlay + 1 DO ig = 1, ngrid wa(ig, k, l) = 0. END DO END DO END DO ! print*,'3 OK convect8' ! ------------------------------------------------------------------ ! Calcul de w2, quarre de w a partir de la cape ! a partir de w2, on calcule wa, vitesse de l'ascendance ! ATTENTION: Dans cette version, pour cause d'economie de memoire, ! w2 est stoke dans wa ! ATTENTION: dans convect8, on n'utilise le calcule des wa ! independants par couches que pour calculer l'entrainement ! a la base et la hauteur max de l'ascendance. ! Indicages: ! l'ascendance provenant du niveau k traverse l'interface l avec ! une vitesse wa(k,l). ! -------------------- ! + + + + + + + + + + ! wa(k,l) ---- -------------------- l ! /\ ! /||\ + + + + + + + + + + ! || ! || -------------------- ! || ! || + + + + + + + + + + ! || ! || -------------------- ! ||__ ! |___ + + + + + + + + + + k ! -------------------- ! ------------------------------------------------------------------ ! CR: ponderation entrainement des couches instables ! def des entr_star tels que entr=f*entr_star DO l = 1, klev DO ig = 1, ngrid entr_star(ig, l) = 0. END DO END DO ! determination de la longueur de la couche d entrainement DO ig = 1, ngrid lentr(ig) = 1 END DO ! on ne considere que les premieres couches instables therm = .FALSE. DO k = nlay - 2, 1, -1 DO ig = 1, ngrid IF (ztv(ig,k)>ztv(ig,k+1) .AND. ztv(ig,k+1)<=ztv(ig,k+2)) THEN lentr(ig) = k + 1 therm = .TRUE. END IF END DO END DO ! limitation de la valeur du lentr ! do ig=1,ngrid ! lentr(ig)=min(5,lentr(ig)) ! enddo ! determination du lmin: couche d ou provient le thermique DO ig = 1, ngrid lmin(ig) = 1 END DO DO ig = 1, ngrid DO l = nlay, 2, -1 IF (ztv(ig,l-1)>ztv(ig,l)) THEN lmin(ig) = l - 1 END IF END DO END DO ! initialisations DO ig = 1, ngrid zalim(ig) = 0. norme(ig) = 0. lalim(ig) = 1 END DO DO k = 1, klev - 1 DO ig = 1, ngrid zalim(ig) = zalim(ig) + zlev(ig, k)*max(0., (ztv(ig,k)-ztv(ig, & k+1))/(zlev(ig,k+1)-zlev(ig,k))) ! s *(zlev(ig,k+1)-zlev(ig,k)) norme(ig) = norme(ig) + max(0., (ztv(ig,k)-ztv(ig,k+1))/(zlev(ig, & k+1)-zlev(ig,k))) ! s *(zlev(ig,k+1)-zlev(ig,k)) END DO END DO DO ig = 1, ngrid IF (norme(ig)>1.E-10) THEN zalim(ig) = max(10.*zalim(ig)/norme(ig), zlev(ig,2)) ! zalim(ig)=min(zalim(ig),zlev(ig,lentr(ig))) END IF END DO ! détermination du lalim correspondant DO k = 1, klev - 1 DO ig = 1, ngrid IF ((zalim(ig)>zlev(ig,k)) .AND. (zalim(ig)<=zlev(ig,k+1))) THEN lalim(ig) = k END IF END DO END DO ! definition de l'entrainement des couches DO l = 1, klev - 1 DO ig = 1, ngrid IF (ztv(ig,l)>ztv(ig,l+1) .AND. l>=lmin(ig) .AND. lztv(ig,l+1) .AND. l>=lmin(ig) .AND. l<=lalim(ig) .AND. & zalim(ig)>1.E-10) THEN ! if (l.le.lentr(ig)) then ! entr_star(ig,l)=zlev(ig,l+1)*(1.-(zlev(ig,l+1) ! s /zalim(ig)))**(3./2.) ! write(10,*)zlev(ig,l),entr_star(ig,l) END IF END DO END DO ! endif ! pas de thermique si couche 1 stable DO ig = 1, ngrid IF (lmin(ig)>5) THEN DO l = 1, klev entr_star(ig, l) = 0. END DO END IF END DO ! calcul de l entrainement total DO ig = 1, ngrid entr_star_tot(ig) = 0. END DO DO ig = 1, ngrid DO k = 1, klev entr_star_tot(ig) = entr_star_tot(ig) + entr_star(ig, k) END DO END DO ! Calcul entrainement normalise DO ig = 1, ngrid IF (entr_star_tot(ig)>1.E-10) THEN ! do l=1,lentr(ig) DO l = 1, klev ! def possibles pour entr_star: zdthetadz, dthetadz, zdtheta entr_star(ig, l) = entr_star(ig, l)/entr_star_tot(ig) END DO END IF END DO ! print*,'fin calcul entr_star' DO k = 1, klev DO ig = 1, ngrid ztva(ig, k) = ztv(ig, k) END DO END DO ! RC ! print*,'7 OK convect8' DO k = 1, klev + 1 DO ig = 1, ngrid zw2(ig, k) = 0. fmc(ig, k) = 0. ! CR f_star(ig, k) = 0. ! RC larg_cons(ig, k) = 0. larg_detr(ig, k) = 0. wa_moy(ig, k) = 0. END DO END DO ! print*,'8 OK convect8' DO ig = 1, ngrid linter(ig) = 1. lmaxa(ig) = 1 lmix(ig) = 1 wmaxa(ig) = 0. END DO ! CR: DO l = 1, nlay - 2 DO ig = 1, ngrid IF (ztv(ig,l)>ztv(ig,l+1) .AND. entr_star(ig,l)>1.E-10 .AND. & zw2(ig,l)<1E-10) THEN f_star(ig, l+1) = entr_star(ig, l) ! test:calcul de dteta zw2(ig, l+1) = 2.*rg*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig, l+1)* & (zlev(ig,l+1)-zlev(ig,l))*0.4*pphi(ig, l)/(pphi(ig,l+1)-pphi(ig,l)) larg_detr(ig, l) = 0. ELSE IF ((zw2(ig,l)>=1E-10) .AND. (f_star(ig,l)+entr_star(ig, & l)>1.E-10)) THEN f_star(ig, l+1) = f_star(ig, l) + entr_star(ig, l) ztva(ig, l) = (f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l)*ztv(ig,l))/ & f_star(ig, l+1) zw2(ig, l+1) = zw2(ig, l)*(f_star(ig,l)/f_star(ig,l+1))**2 + & 2.*rg*(ztva(ig,l)-ztv(ig,l))/ztv(ig, l)*(zlev(ig,l+1)-zlev(ig,l)) END IF ! determination de zmax continu par interpolation lineaire IF (zw2(ig,l+1)<0.) THEN ! test IF (abs(zw2(ig,l+1)-zw2(ig,l))<1E-10) THEN ! print*,'pb linter' END IF linter(ig) = (l*(zw2(ig,l+1)-zw2(ig,l))-zw2(ig,l))/(zw2(ig,l+1)-zw2( & ig,l)) zw2(ig, l+1) = 0. lmaxa(ig) = l ELSE IF (zw2(ig,l+1)<0.) THEN ! print*,'pb1 zw2<0' END IF wa_moy(ig, l+1) = sqrt(zw2(ig,l+1)) END IF IF (wa_moy(ig,l+1)>wmaxa(ig)) THEN ! lmix est le niveau de la couche ou w (wa_moy) est maximum lmix(ig) = l + 1 wmaxa(ig) = wa_moy(ig, l+1) END IF END DO END DO ! print*,'fin calcul zw2' ! Calcul de la couche correspondant a la hauteur du thermique DO ig = 1, ngrid lmax(ig) = lentr(ig) ! lmax(ig)=lalim(ig) END DO DO ig = 1, ngrid DO l = nlay, lentr(ig) + 1, -1 ! do l=nlay,lalim(ig)+1,-1 IF (zw2(ig,l)<=1.E-10) THEN lmax(ig) = l - 1 END IF END DO END DO ! pas de thermique si couche 1 stable DO ig = 1, ngrid IF (lmin(ig)>5) THEN lmax(ig) = 1 lmin(ig) = 1 lentr(ig) = 1 lalim(ig) = 1 END IF END DO ! Determination de zw2 max DO ig = 1, ngrid wmax(ig) = 0. END DO DO l = 1, nlay DO ig = 1, ngrid IF (l<=lmax(ig)) THEN IF (zw2(ig,l)<0.) THEN ! print*,'pb2 zw2<0' END IF zw2(ig, l) = sqrt(zw2(ig,l)) wmax(ig) = max(wmax(ig), zw2(ig,l)) ELSE zw2(ig, l) = 0. END IF END DO END DO ! Longueur caracteristique correspondant a la hauteur des thermiques. DO ig = 1, ngrid zmax(ig) = 0. zlevinter(ig) = zlev(ig, 1) END DO DO ig = 1, ngrid ! calcul de zlevinter zlevinter(ig) = (zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*linter(ig) + & zlev(ig, lmax(ig)) - lmax(ig)*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig))) zmax(ig) = max(zmax(ig), zlevinter(ig)-zlev(ig,lmin(ig))) END DO DO ig = 1, ngrid ! write(8,*)zmax(ig),lmax(ig),lentr(ig),lmin(ig) END DO ! on stope après les calculs de zmax et wmax RETURN ! print*,'avant fermeture' ! Fermeture,determination de f ! Attention! entrainement normalisé ou pas? DO ig = 1, ngrid entr_star2(ig) = 0. END DO DO ig = 1, ngrid IF (entr_star_tot(ig)<1.E-10) THEN f(ig) = 0. ELSE DO k = lmin(ig), lentr(ig) ! do k=lmin(ig),lalim(ig) entr_star2(ig) = entr_star2(ig) + entr_star(ig, k)**2/(rho(ig,k)*( & zlev(ig,k+1)-zlev(ig,k))) END DO ! Nouvelle fermeture f(ig) = wmax(ig)/(max(500.,zmax(ig))*r_aspect*entr_star2(ig)) ! s *entr_star_tot(ig) ! test ! if (first) then f(ig) = f(ig) + (f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)*wmax(ig)) ! endif END IF f0(ig) = f(ig) ! first=.true. END DO ! print*,'apres fermeture' ! on stoppe après la fermeture RETURN ! Calcul de l'entrainement DO k = 1, klev DO ig = 1, ngrid entr(ig, k) = f(ig)*entr_star(ig, k) END DO END DO ! on stoppe après le calcul de entr ! RETURN ! CR:test pour entrainer moins que la masse ! do ig=1,ngrid ! do l=1,lentr(ig) ! if ((entr(ig,l)*ptimestep).gt.(0.9*masse(ig,l))) then ! entr(ig,l+1)=entr(ig,l+1)+entr(ig,l) ! s -0.9*masse(ig,l)/ptimestep ! entr(ig,l)=0.9*masse(ig,l)/ptimestep ! endif ! enddo ! enddo ! CR: fin test ! Calcul des flux DO ig = 1, ngrid DO l = 1, lmax(ig) - 1 fmc(ig, l+1) = fmc(ig, l) + entr(ig, l) END DO END DO ! RC ! print*,'9 OK convect8' ! print*,'WA1 ',wa_moy ! determination de l'indice du debut de la mixed layer ou w decroit ! calcul de la largeur de chaque ascendance dans le cas conservatif. ! dans ce cas simple, on suppose que la largeur de l'ascendance provenant ! d'une couche est égale à la hauteur de la couche alimentante. ! La vitesse maximale dans l'ascendance est aussi prise comme estimation ! de la vitesse d'entrainement horizontal dans la couche alimentante. DO l = 2, nlay DO ig = 1, ngrid IF (l<=lmaxa(ig)) THEN zw = max(wa_moy(ig,l), 1.E-10) larg_cons(ig, l) = zmax(ig)*r_aspect*fmc(ig, l)/(rhobarz(ig,l)*zw) END IF END DO END DO DO l = 2, nlay DO ig = 1, ngrid IF (l<=lmaxa(ig)) THEN ! if (idetr.eq.0) then ! cette option est finalement en dur. IF ((l_mix*zlev(ig,l))<0.) THEN ! print*,'pb l_mix*zlev<0' END IF ! CR: test: nouvelle def de lambda ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l)) IF (zw2(ig,l)>1.E-10) THEN larg_detr(ig, l) = sqrt((l_mix/zw2(ig,l))*zlev(ig,l)) ELSE larg_detr(ig, l) = sqrt(l_mix*zlev(ig,l)) END IF ! RC ! else if (idetr.eq.1) then ! larg_detr(ig,l)=larg_cons(ig,l) ! s *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig)) ! else if (idetr.eq.2) then ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l)) ! s *sqrt(wa_moy(ig,l)) ! else if (idetr.eq.4) then ! larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l)) ! s *wa_moy(ig,l) ! endif END IF END DO END DO ! print*,'10 OK convect8' ! print*,'WA2 ',wa_moy ! calcul de la fraction de la maille concernée par l'ascendance en tenant ! compte de l'epluchage du thermique. ! CR def de zmix continu (profil parabolique des vitesses) DO ig = 1, ngrid IF (lmix(ig)>1.) THEN ! test IF (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)))- & (zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))- & zw2(ig,lmix(ig)+1))*((zlev(ig,lmix(ig)-1))- & (zlev(ig,lmix(ig)))))>1E-10) THEN zmix(ig) = ((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)) & )**2-(zlev(ig,lmix(ig)+1))**2)-(zw2(ig,lmix(ig))-zw2(ig, & lmix(ig)+1))*((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))/ & (2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))*((zlev(ig,lmix(ig)))- & (zlev(ig,lmix(ig)+1)))-(zw2(ig,lmix(ig))- & zw2(ig,lmix(ig)+1))*((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig)))))) ELSE zmix(ig) = zlev(ig, lmix(ig)) ! print*,'pb zmix' END IF ELSE zmix(ig) = 0. END IF ! test IF ((zmax(ig)-zmix(ig))<0.) THEN zmix(ig) = 0.99*zmax(ig) ! print*,'pb zmix>zmax' END IF END DO ! calcul du nouveau lmix correspondant DO ig = 1, ngrid DO l = 1, klev IF (zmix(ig)>=zlev(ig,l) .AND. zmix(ig)1.) THEN ! print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK' fraca(ig, l) = (larg_cons(ig,l)-larg_detr(ig,l))/(r_aspect*zmax(ig)) ! test fraca(ig, l) = max(fraca(ig,l), 0.) fraca(ig, l) = min(fraca(ig,l), 0.5) fracd(ig, l) = 1. - fraca(ig, l) fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig)) ELSE ! wa_moy(ig,l)=0. fraca(ig, l) = 0. fracc(ig, l) = 0. fracd(ig, l) = 1. END IF END DO END DO ! CR: calcul de fracazmix DO ig = 1, ngrid fracazmix(ig) = (fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/ & (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig) + & fraca(ig, lmix(ig)) - zlev(ig, lmix(ig))*(fraca(ig,lmix(ig)+1)-fraca(ig & ,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig))) END DO DO l = 2, nlay DO ig = 1, ngrid IF (larg_cons(ig,l)>1.) THEN IF (l>lmix(ig)) THEN ! test IF (zmax(ig)-zmix(ig)<1.E-10) THEN ! print*,'pb xxx' xxx(ig, l) = (lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig)) ELSE xxx(ig, l) = (zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig)) END IF IF (idetr==0) THEN fraca(ig, l) = fracazmix(ig) ELSE IF (idetr==1) THEN fraca(ig, l) = fracazmix(ig)*xxx(ig, l) ELSE IF (idetr==2) THEN fraca(ig, l) = fracazmix(ig)*(1.-(1.-xxx(ig,l))**2) ELSE fraca(ig, l) = fracazmix(ig)*xxx(ig, l)**2 END IF ! print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL' fraca(ig, l) = max(fraca(ig,l), 0.) fraca(ig, l) = min(fraca(ig,l), 0.5) fracd(ig, l) = 1. - fraca(ig, l) fracc(ig, l) = larg_cons(ig, l)/(r_aspect*zmax(ig)) END IF END IF END DO END DO ! print*,'fin calcul fraca' ! print*,'11 OK convect8' ! print*,'Ea3 ',wa_moy ! ------------------------------------------------------------------ ! Calcul de fracd, wd ! somme wa - wd = 0 ! ------------------------------------------------------------------ DO ig = 1, ngrid fm(ig, 1) = 0. fm(ig, nlay+1) = 0. END DO DO l = 2, nlay DO ig = 1, ngrid fm(ig, l) = fraca(ig, l)*wa_moy(ig, l)*rhobarz(ig, l) ! CR:test IF (entr(ig,l-1)<1E-10 .AND. fm(ig,l)>fm(ig,l-1) .AND. l>lmix(ig)) THEN fm(ig, l) = fm(ig, l-1) ! write(1,*)'ajustement fm, l',l END IF ! write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l) ! RC END DO DO ig = 1, ngrid IF (fracd(ig,l)<0.1) THEN abort_message = 'fracd trop petit' CALL abort_physic(modname, abort_message, 1) ELSE ! vitesse descendante "diagnostique" wd(ig, l) = fm(ig, l)/(fracd(ig,l)*rhobarz(ig,l)) END IF END DO END DO DO l = 1, nlay DO ig = 1, ngrid ! masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) masse(ig, l) = (pplev(ig,l)-pplev(ig,l+1))/rg END DO END DO ! print*,'12 OK convect8' ! print*,'WA4 ',wa_moy ! c------------------------------------------------------------------ ! calcul du transport vertical ! ------------------------------------------------------------------ GO TO 4444 ! print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep DO l = 2, nlay - 1 DO ig = 1, ngrid IF (fm(ig,l+1)*ptimestep>masse(ig,l) .AND. fm(ig,l+1)*ptimestep>masse( & ig,l+1)) THEN ! print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM=' ! s ,fm(ig,l+1)*ptimestep ! s ,' M=',masse(ig,l),masse(ig,l+1) END IF END DO END DO DO l = 1, nlay DO ig = 1, ngrid IF (entr(ig,l)*ptimestep>masse(ig,l)) THEN ! print*,'WARN!!! E>M ig=',ig,' l=',l,' E==' ! s ,entr(ig,l)*ptimestep ! s ,' M=',masse(ig,l) END IF END DO END DO DO l = 1, nlay DO ig = 1, ngrid IF (.NOT. fm(ig,l)>=0. .OR. .NOT. fm(ig,l)<=10.) THEN ! print*,'WARN!!! fm exagere ig=',ig,' l=',l ! s ,' FM=',fm(ig,l) END IF IF (.NOT. masse(ig,l)>=1.E-10 .OR. .NOT. masse(ig,l)<=1.E4) THEN ! print*,'WARN!!! masse exagere ig=',ig,' l=',l ! s ,' M=',masse(ig,l) ! print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)', ! s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l) ! print*,'zlev(ig,l+1),zlev(ig,l)' ! s ,zlev(ig,l+1),zlev(ig,l) ! print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)' ! s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1) END IF IF (.NOT. entr(ig,l)>=0. .OR. .NOT. entr(ig,l)<=10.) THEN ! print*,'WARN!!! entr exagere ig=',ig,' l=',l ! s ,' E=',entr(ig,l) END IF END DO END DO 4444 CONTINUE ! CR:redefinition du entr DO l = 1, nlay DO ig = 1, ngrid detr(ig, l) = fm(ig, l) + entr(ig, l) - fm(ig, l+1) IF (detr(ig,l)<0.) THEN ! entr(ig,l)=entr(ig,l)-detr(ig,l) fm(ig, l+1) = fm(ig, l) + entr(ig, l) detr(ig, l) = 0. ! print*,'WARNING !!! detrainement negatif ',ig,l END IF END DO END DO ! RC IF (w2di==1) THEN fm0 = fm0 + ptimestep*(fm-fm0)/tho entr0 = entr0 + ptimestep*(entr-entr0)/tho ELSE fm0 = fm entr0 = entr END IF IF (1==1) THEN CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zh, zdhadj, & zha) CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zo, pdoadj, & zoa) ELSE CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zh, & zdhadj, zha) CALL dqthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zo, & pdoadj, zoa) END IF IF (1==0) THEN CALL dvthermcell2(ngrid, nlay, ptimestep, fm0, entr0, masse, fraca, zmax, & zu, zv, pduadj, pdvadj, zua, zva) ELSE CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zu, pduadj, & zua) CALL dqthermcell(ngrid, nlay, ptimestep, fm0, entr0, masse, zv, pdvadj, & zva) END IF DO l = 1, nlay DO ig = 1, ngrid zf = 0.5*(fracc(ig,l)+fracc(ig,l+1)) zf2 = zf/(1.-zf) thetath2(ig, l) = zf2*(zha(ig,l)-zh(ig,l))**2 wth2(ig, l) = zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2 END DO END DO ! print*,'13 OK convect8' ! print*,'WA5 ',wa_moy DO l = 1, nlay DO ig = 1, ngrid pdtadj(ig, l) = zdhadj(ig, l)*zpspsk(ig, l) END DO END DO ! do l=1,nlay ! do ig=1,ngrid ! if(abs(pdtadj(ig,l))*86400..gt.500.) then ! print*,'WARN!!! ig=',ig,' l=',l ! s ,' pdtadj=',pdtadj(ig,l) ! endif ! if(abs(pdoadj(ig,l))*86400..gt.1.) then ! print*,'WARN!!! ig=',ig,' l=',l ! s ,' pdoadj=',pdoadj(ig,l) ! endif ! enddo ! enddo ! print*,'14 OK convect8' ! ------------------------------------------------------------------ ! Calculs pour les sorties ! ------------------------------------------------------------------ IF (sorties) THEN DO l = 1, nlay DO ig = 1, ngrid zla(ig, l) = (1.-fracd(ig,l))*zmax(ig) zld(ig, l) = fracd(ig, l)*zmax(ig) IF (1.-fracd(ig,l)>1.E-10) zwa(ig, l) = wd(ig, l)*fracd(ig, l)/ & (1.-fracd(ig,l)) END DO END DO ! deja fait ! do l=1,nlay ! do ig=1,ngrid ! detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1) ! if (detr(ig,l).lt.0.) then ! entr(ig,l)=entr(ig,l)-detr(ig,l) ! detr(ig,l)=0. ! print*,'WARNING !!! detrainement negatif ',ig,l ! endif ! enddo ! enddo ! print*,'15 OK convect8' isplit = isplit + 1 ! #define und GO TO 123 #ifdef und CALL writeg1d(1, nlay, wd, 'wd ', 'wd ') CALL writeg1d(1, nlay, zwa, 'wa ', 'wa ') CALL writeg1d(1, nlay, fracd, 'fracd ', 'fracd ') CALL writeg1d(1, nlay, fraca, 'fraca ', 'fraca ') CALL writeg1d(1, nlay, wa_moy, 'wam ', 'wam ') CALL writeg1d(1, nlay, zla, 'la ', 'la ') CALL writeg1d(1, nlay, zld, 'ld ', 'ld ') CALL writeg1d(1, nlay, pt, 'pt ', 'pt ') CALL writeg1d(1, nlay, zh, 'zh ', 'zh ') CALL writeg1d(1, nlay, zha, 'zha ', 'zha ') CALL writeg1d(1, nlay, zu, 'zu ', 'zu ') CALL writeg1d(1, nlay, zv, 'zv ', 'zv ') CALL writeg1d(1, nlay, zo, 'zo ', 'zo ') CALL writeg1d(1, nlay, wh, 'wh ', 'wh ') CALL writeg1d(1, nlay, wu, 'wu ', 'wu ') CALL writeg1d(1, nlay, wv, 'wv ', 'wv ') CALL writeg1d(1, nlay, wo, 'w15uo ', 'wXo ') CALL writeg1d(1, nlay, zdhadj, 'zdhadj ', 'zdhadj ') CALL writeg1d(1, nlay, pduadj, 'pduadj ', 'pduadj ') CALL writeg1d(1, nlay, pdvadj, 'pdvadj ', 'pdvadj ') CALL writeg1d(1, nlay, pdoadj, 'pdoadj ', 'pdoadj ') CALL writeg1d(1, nlay, entr, 'entr ', 'entr ') CALL writeg1d(1, nlay, detr, 'detr ', 'detr ') CALL writeg1d(1, nlay, fm, 'fm ', 'fm ') CALL writeg1d(1, nlay, pdtadj, 'pdtadj ', 'pdtadj ') CALL writeg1d(1, nlay, pplay, 'pplay ', 'pplay ') CALL writeg1d(1, nlay, pplev, 'pplev ', 'pplev ') ! recalcul des flux en diagnostique... ! print*,'PAS DE TEMPS ',ptimestep CALL dt2f(pplev, pplay, pt, pdtadj, wh) CALL writeg1d(1, nlay, wh, 'wh2 ', 'wh2 ') #endif 123 CONTINUE END IF ! if(wa_moy(1,4).gt.1.e-10) stop ! print*,'19 OK convect8' RETURN END SUBROUTINE calcul_sec SUBROUTINE fermeture_seche(ngrid, nlay, pplay, pplev, pphi, zlev, rhobarz, & f0, zpspsk, alim_star, zh, zo, lentr, lmin, nu_min, nu_max, r_aspect, & zmax, wmax) USE dimphy IMPLICIT NONE include "YOMCST.h" INTEGER ngrid, nlay REAL pplay(ngrid, nlay), pplev(ngrid, nlay+1) REAL pphi(ngrid, nlay) REAL zlev(klon, klev+1) REAL alim_star(klon, klev) REAL f0(klon) INTEGER lentr(klon) INTEGER lmin(klon) REAL zmax(klon) REAL wmax(klon) REAL nu_min REAL nu_max REAL r_aspect REAL rhobarz(klon, klev+1) REAL zh(klon, klev) REAL zo(klon, klev) REAL zpspsk(klon, klev) INTEGER ig, l REAL f_star(klon, klev+1) REAL detr_star(klon, klev) REAL entr_star(klon, klev) REAL zw2(klon, klev+1) REAL linter(klon) INTEGER lmix(klon) INTEGER lmax(klon) REAL zlevinter(klon) REAL wa_moy(klon, klev+1) REAL wmaxa(klon) REAL ztv(klon, klev) REAL ztva(klon, klev) REAL nu(klon, klev) ! real zmax0_sec(klon) ! save zmax0_sec REAL, SAVE, ALLOCATABLE :: zmax0_sec(:) !$OMP THREADPRIVATE(zmax0_sec) LOGICAL, SAVE :: first = .TRUE. !$OMP THREADPRIVATE(first) IF (first) THEN ALLOCATE (zmax0_sec(klon)) first = .FALSE. END IF DO l = 1, nlay DO ig = 1, ngrid ztv(ig, l) = zh(ig, l)/zpspsk(ig, l) ztv(ig, l) = ztv(ig, l)*(1.+retv*zo(ig,l)) END DO END DO DO l = 1, nlay - 2 DO ig = 1, ngrid IF (ztv(ig,l)>ztv(ig,l+1) .AND. alim_star(ig,l)>1.E-10 .AND. & zw2(ig,l)<1E-10) THEN f_star(ig, l+1) = alim_star(ig, l) ! test:calcul de dteta zw2(ig, l+1) = 2.*rg*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig, l+1)* & (zlev(ig,l+1)-zlev(ig,l))*0.4*pphi(ig, l)/(pphi(ig,l+1)-pphi(ig,l)) ELSE IF ((zw2(ig,l)>=1E-10) .AND. (f_star(ig,l)+alim_star(ig, & l))>1.E-10) THEN ! estimation du detrainement a partir de la geometrie du pas ! precedent ! tests sur la definition du detr nu(ig, l) = (nu_min+nu_max)/2.*(1.-(nu_max-nu_min)/(nu_max+nu_min)* & tanh((((ztva(ig,l-1)-ztv(ig,l))/ztv(ig,l))/0.0005))) detr_star(ig, l) = rhobarz(ig, l)*sqrt(zw2(ig,l))/ & (r_aspect*zmax0_sec(ig))* & ! s ! /(r_aspect*zmax0(ig))* (sqrt(nu(ig,l)*zlev(ig,l+1)/sqrt(zw2(ig,l)))-sqrt(nu(ig,l)*zlev(ig, & l)/sqrt(zw2(ig,l)))) detr_star(ig, l) = detr_star(ig, l)/f0(ig) IF ((detr_star(ig,l))>f_star(ig,l)) THEN detr_star(ig, l) = f_star(ig, l) END IF entr_star(ig, l) = 0.9*detr_star(ig, l) IF ((l1.E-10) THEN ! AM on melange Tl et qt du thermique ztva(ig, l) = (f_star(ig,l)*ztva(ig,l-1)+(entr_star(ig, & l)+alim_star(ig,l))*ztv(ig,l))/(f_star(ig,l+1)+detr_star(ig,l)) zw2(ig, l+1) = zw2(ig, l)*(f_star(ig,l)/(f_star(ig, & l+1)+detr_star(ig,l)))**2 + 2.*rg*(ztva(ig,l)-ztv(ig,l))/ztv(ig, & l)*(zlev(ig,l+1)-zlev(ig,l)) END IF END IF IF (zw2(ig,l+1)<0.) THEN linter(ig) = (l*(zw2(ig,l+1)-zw2(ig,l))-zw2(ig,l))/(zw2(ig,l+1)-zw2( & ig,l)) zw2(ig, l+1) = 0. ! print*,'linter=',linter(ig) ELSE wa_moy(ig, l+1) = sqrt(zw2(ig,l+1)) END IF IF (wa_moy(ig,l+1)>wmaxa(ig)) THEN ! lmix est le niveau de la couche ou w (wa_moy) est maximum lmix(ig) = l + 1 wmaxa(ig) = wa_moy(ig, l+1) END IF END DO END DO ! print*,'fin calcul zw2' ! Calcul de la couche correspondant a la hauteur du thermique DO ig = 1, ngrid lmax(ig) = lentr(ig) END DO DO ig = 1, ngrid DO l = nlay, lentr(ig) + 1, -1 IF (zw2(ig,l)<=1.E-10) THEN lmax(ig) = l - 1 END IF END DO END DO ! pas de thermique si couche 1 stable DO ig = 1, ngrid IF (lmin(ig)>1) THEN lmax(ig) = 1 lmin(ig) = 1 lentr(ig) = 1 END IF END DO ! Determination de zw2 max DO ig = 1, ngrid wmax(ig) = 0. END DO DO l = 1, nlay DO ig = 1, ngrid IF (l<=lmax(ig)) THEN IF (zw2(ig,l)<0.) THEN ! print*,'pb2 zw2<0' END IF zw2(ig, l) = sqrt(zw2(ig,l)) wmax(ig) = max(wmax(ig), zw2(ig,l)) ELSE zw2(ig, l) = 0. END IF END DO END DO ! Longueur caracteristique correspondant a la hauteur des thermiques. DO ig = 1, ngrid zmax(ig) = 0. zlevinter(ig) = zlev(ig, 1) END DO DO ig = 1, ngrid ! calcul de zlevinter zlevinter(ig) = (zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*linter(ig) + & zlev(ig, lmax(ig)) - lmax(ig)*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig))) ! pour le cas ou on prend tjs lmin=1 ! zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig))) zmax(ig) = max(zmax(ig), zlevinter(ig)-zlev(ig,1)) zmax0_sec(ig) = zmax(ig) END DO RETURN END SUBROUTINE fermeture_seche