Index: LMDZ5/trunk/libf/phylmd/thermcell_main.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/thermcell_main.F90	(revision 2351)
+++ LMDZ5/trunk/libf/phylmd/thermcell_main.F90	(revision 2384)
@@ -747,302 +747,24 @@
       enddo
 !
-
+! $Id$
+!
+      print*,'THERM_ALP sorti'
+      CALL thermcell_alp(ngrid,nlay,ptimestep  &
+     &                  ,pplay,pplev  &
+     &                  ,fm0,entr0,lmax  &
+     &                  ,Ale_bl,Alp_bl,lalim_conv,wght_th &
+     &                  ,zw2,fraca &
+!!! necessire en plus
+     &                  ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &
 !!! nrlmd le 10/04/2012
-
-!------------Test sur le LCL des thermiques
-    do ig=1,ngrid
-      ok_lcl(ig)=.false.
-      if ( (pcon(ig) .gt. pplay(ig,klev-1)) .and. (pcon(ig) .lt. pplay(ig,1)) ) ok_lcl(ig)=.true.
-    enddo
-
-!------------Localisation des niveaux entourant le LCL et du coef d'interpolation
-    do l=1,nlay-1
-      do ig=1,ngrid
-        if (ok_lcl(ig)) then 
-!ATTENTION,zw2 calcule en pplev
-!          if ((pplay(ig,l) .ge. pcon(ig)) .and. (pplay(ig,l+1) .le. pcon(ig))) then
-!          klcl(ig)=l
-!          interp(ig)=(pcon(ig)-pplay(ig,klcl(ig)))/(pplay(ig,klcl(ig)+1)-pplay(ig,klcl(ig)))
-!          endif
-          if ((pplev(ig,l) .ge. pcon(ig)) .and. (pplev(ig,l+1) .le. pcon(ig))) then
-          klcl(ig)=l
-          interp(ig)=(pcon(ig)-pplev(ig,klcl(ig)))/(pplev(ig,klcl(ig)+1)-pplev(ig,klcl(ig)))
-          endif
-        endif
-      enddo
-    enddo
-
-!------------Hauteur des thermiques
-!!jyg le 27/04/2012
-!!    do ig =1,ngrid
-!!    rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
-!! &               -rhobarz(ig,klcl(ig)))*interp(ig)
-!!    zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
-!!    zmax(ig)=pphi(ig,lmax(ig))/rg
-!!      if ( (.not.ok_lcl(ig)) .or. (zlcl(ig).gt.zmax(ig)) ) zlcl(ig)=zmax(ig) ! Si zclc > zmax alors on pose zlcl = zmax
-!!    enddo
-    do ig =1,ngrid
-!CR:REHABILITATION ZMAX CONTINU
-!     zmax(ig)=pphi(ig,lmax(ig))/rg
-     if (ok_lcl(ig)) then 
-      rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
- &               -rhobarz(ig,klcl(ig)))*interp(ig)
-      zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
-      zlcl(ig)=min(zlcl(ig),zmax(ig))   ! Si zlcl > zmax alors on pose zlcl = zmax
-     else
-      rhobarz0(ig)=0.
-      zlcl(ig)=zmax(ig)
-     endif
-    enddo
-!!jyg fin
-
-!------------Calcul des propriétés du thermique au LCL 
-  IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) THEN 
-
-  !-----Initialisation de la TKE moyenne 
-   do l=1,nlay
-    do ig=1,ngrid
-     pbl_tke_max(ig,l)=0.
-    enddo
-   enddo
-
-!-----Calcul de la TKE moyenne 
-   do nsrf=1,nbsrf
-    do l=1,nlay
-     do ig=1,ngrid
-     pbl_tke_max(ig,l)=pctsrf(ig,nsrf)*pbl_tke(ig,l,nsrf)+pbl_tke_max(ig,l)
-     enddo
-    enddo
-   enddo
-
-!-----Initialisations des TKE dans et hors des thermiques 
-   do l=1,nlay
-    do ig=1,ngrid
-    therm_tke_max(ig,l)=pbl_tke_max(ig,l)
-    env_tke_max(ig,l)=pbl_tke_max(ig,l)
-    enddo
-   enddo
-
-!-----Calcul de la TKE transportée par les thermiques : therm_tke_max
-   call thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,  &
-  &           rg,pplev,therm_tke_max)
-!   print *,' thermcell_tke_transport -> '   !!jyg
-
-!-----Calcul des profils verticaux de TKE hors thermiques : env_tke_max, et de la vitesse verticale grande échelle : W_ls
-   do l=1,nlay
-    do ig=1,ngrid
-     pbl_tke_max(ig,l)=fraca(ig,l)*therm_tke_max(ig,l)+(1.-fraca(ig,l))*env_tke_max(ig,l)         !  Recalcul de TKE moyenne aprés transport de TKE_TH
-     env_tke_max(ig,l)=(pbl_tke_max(ig,l)-fraca(ig,l)*therm_tke_max(ig,l))/(1.-fraca(ig,l))       !  Recalcul de TKE dans  l'environnement aprés transport de TKE_TH
-     w_ls(ig,l)=-1.*omega(ig,l)/(RG*rhobarz(ig,l))                                                !  Vitesse verticale de grande échelle
-    enddo
-   enddo
-!    print *,' apres w_ls = '   !!jyg
-
-  do ig=1,ngrid
-   if (ok_lcl(ig)) then
-     fraca0(ig)=fraca(ig,klcl(ig))+(fraca(ig,klcl(ig)+1) &
- &             -fraca(ig,klcl(ig)))*interp(ig)
-     w0(ig)=zw2(ig,klcl(ig))+(zw2(ig,klcl(ig)+1) &
- &         -zw2(ig,klcl(ig)))*interp(ig)
-     w_conv(ig)=w_ls(ig,klcl(ig))+(w_ls(ig,klcl(ig)+1) &
- &             -w_ls(ig,klcl(ig)))*interp(ig)
-     therm_tke_max0(ig)=therm_tke_max(ig,klcl(ig)) &
- &                     +(therm_tke_max(ig,klcl(ig)+1)-therm_tke_max(ig,klcl(ig)))*interp(ig)
-     env_tke_max0(ig)=env_tke_max(ig,klcl(ig))+(env_tke_max(ig,klcl(ig)+1) &
- &                   -env_tke_max(ig,klcl(ig)))*interp(ig)
-     pbl_tke_max0(ig)=pbl_tke_max(ig,klcl(ig))+(pbl_tke_max(ig,klcl(ig)+1) &
- &                   -pbl_tke_max(ig,klcl(ig)))*interp(ig)
-     if (therm_tke_max0(ig).ge.20.) therm_tke_max0(ig)=20.
-     if (env_tke_max0(ig).ge.20.) env_tke_max0(ig)=20.
-     if (pbl_tke_max0(ig).ge.20.) pbl_tke_max0(ig)=20.
-   else 
-     fraca0(ig)=0.
-     w0(ig)=0.
-!!jyg le 27/04/2012
-!!     zlcl(ig)=0.
-!!
-   endif
-  enddo
-
-  ENDIF ! IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) )
-!  print *,'ENDIF  ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) '    !!jyg
-
-!------------Triggering------------------
-  IF (iflag_trig_bl.ge.1) THEN 
-
-!-----Initialisations
-   depth(:)=0.
-   n2(:)=0.
-   s2(:)=100. ! some low value, arbitrary
-   s_max(:)=0.
-
-!-----Epaisseur du nuage (depth) et détermination de la queue du spectre de panaches (n2,s2) et du panache le plus gros (s_max)
-   do ig=1,ngrid
-     zmax_moy(ig)=zlcl(ig)+zmax_moy_coef*(zmax(ig)-zlcl(ig))
-     depth(ig)=zmax_moy(ig)-zlcl(ig)
-     hmin(ig)=hmincoef*zlcl(ig)
-     if (depth(ig).ge.10.) then 
-       s2(ig)=(hcoef*depth(ig)+hmin(ig))**2
-       n2(ig)=(1.-eps1)*fraca0(ig)*airephy(ig)/s2(ig)
-!!
-!!jyg le 27/04/2012
-!!       s_max(ig)=s2(ig)*log(n2(ig))
-!!       if (n2(ig) .lt. 1) s_max(ig)=0.
-       s_max(ig)=s2(ig)*log(max(n2(ig),1.))
-!!fin jyg
-     else
-       n2(ig)=0.
-       s_max(ig)=0.
-     endif
-   enddo
-!   print *,'avant Calcul de Wmax '    !!jyg
-
-!-----Calcul de Wmax et ALE_BL_STAT associée
-!!jyg le 30/04/2012
-!!   do ig=1,ngrid
-!!     if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.1.) ) then
-!!     w_max(ig)=w0(ig)*(1.+sqrt(2.*log(s_max(ig)/su)-log(2.*3.14)-log(2.*log(s_max(ig)/su)-log(2.*3.14))))
-!!     ale_bl_stat(ig)=0.5*w_max(ig)**2
-!!     else
-!!     w_max(ig)=0.
-!!     ale_bl_stat(ig)=0.
-!!     endif
-!!   enddo
-   susqr2pi=su*sqrt(2.*Rpi)
-   Reuler=exp(1.)
-   do ig=1,ngrid
-     if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.susqr2pi*Reuler) ) then
-      w_max(ig)=w0(ig)*(1.+sqrt(2.*log(s_max(ig)/susqr2pi)-log(2.*log(s_max(ig)/susqr2pi))))
-      ale_bl_stat(ig)=0.5*w_max(ig)**2
-     else
-      w_max(ig)=0.
-      ale_bl_stat(ig)=0.
-     endif
-   enddo
-
-  ENDIF ! iflag_trig_bl
-!  print *,'ENDIF  iflag_trig_bl'    !!jyg
-
-!------------Closure------------------
-
-  IF (iflag_clos_bl.ge.2) THEN 
-
-!-----Calcul de ALP_BL_STAT
-  do ig=1,ngrid
-  alp_bl_det(ig)=0.5*coef_m*rhobarz0(ig)*(w0(ig)**3)*fraca0(ig)*(1.-2.*fraca0(ig))/((1.-fraca0(ig))**2)
-  alp_bl_fluct_m(ig)=1.5*rhobarz0(ig)*fraca0(ig)*(w_conv(ig)+coef_m*w0(ig))* &
- &                   (w0(ig)**2)
-  alp_bl_fluct_tke(ig)=3.*coef_m*rhobarz0(ig)*w0(ig)*fraca0(ig)*(therm_tke_max0(ig)-env_tke_max0(ig)) &
- &                    +3.*rhobarz0(ig)*w_conv(ig)*pbl_tke_max0(ig)
-    if (iflag_clos_bl.ge.2) then 
-    alp_bl_conv(ig)=1.5*coef_m*rhobarz0(ig)*fraca0(ig)*(fraca0(ig)/(1.-fraca0(ig)))*w_conv(ig)* &
- &                   (w0(ig)**2)
-    else
-    alp_bl_conv(ig)=0.
-    endif
-  alp_bl_stat(ig)=alp_bl_det(ig)+alp_bl_fluct_m(ig)+alp_bl_fluct_tke(ig)+alp_bl_conv(ig)
-  enddo
-
-!-----Sécurité ALP infinie
-  do ig=1,ngrid
-   if (fraca0(ig).gt.0.98) alp_bl_stat(ig)=2.
-  enddo
-
-  ENDIF ! (iflag_clos_bl.ge.2)
-
+     &                  ,pbl_tke,pctsrf,omega,airephy &
+     &                  ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
+     &                  ,n2,s2,ale_bl_stat &
+     &                  ,therm_tke_max,env_tke_max &
+     &                  ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
+     &                  ,alp_bl_conv,alp_bl_stat &
 !!! fin nrlmd le 10/04/2012
-
-      if (prt_level.ge.10) then
-         ig=igout
-         do l=1,nlay
-            print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
-            print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
-         enddo
-      endif
-
-!      print*,'avant calcul ale et alp' 
-!calcul de ALE et ALP pour la convection
-      Alp_bl(:)=0.
-      Ale_bl(:)=0.
-!          print*,'ALE,ALP ,l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)'
-      do l=1,nlay
-      do ig=1,ngrid
-           Alp_bl(ig)=max(Alp_bl(ig),0.5*rhobarz(ig,l)*wth3(ig,l) )
-           Ale_bl(ig)=max(Ale_bl(ig),0.5*zw2(ig,l)**2)
-!          print*,'ALE,ALP',l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)
-      enddo
-      enddo
-
-! Ale sec (max de wmax/2 sous la zone d'inhibition) dans
-! le cas iflag_trig_bl=3
-      IF (iflag_trig_bl==3) Ale_bl(:)=0.5*wmax_sec(:)**2
-
-!test:calcul de la ponderation des couches pour KE
-!initialisations
-
-      fm_tot(:)=0.
-      wght_th(:,:)=1.
-      lalim_conv(:)=lalim(:)
-
-      do k=1,klev
-         do ig=1,ngrid
-            if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
-         enddo
-      enddo
-
-! assez bizarre car, si on est dans la couche d'alim et que alim_star et
-! plus petit que 1.e-10, on prend wght_th=1.
-      do k=1,klev
-         do ig=1,ngrid
-            if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then
-               wght_th(ig,k)=alim_star(ig,k)
-            endif
-         enddo
-      enddo
-
-!      print*,'apres wght_th'
-!test pour prolonger la convection
-      do ig=1,ngrid
-!v1d  if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
-      if ((alim_star(ig,1).lt.1.e-10)) then
-      lalim_conv(ig)=1
-      wght_th(ig,1)=1.
-!      print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
-      endif
-      enddo
-
-!------------------------------------------------------------------------
-! Modif CR/FH 20110310 : Alp integree sur la verticale.
-! Integrale verticale de ALP.
-! wth3 etant aux niveaux inter-couches, on utilise d play comme masse des
-! couches
-!------------------------------------------------------------------------
-
-      alp_int(:)=0.
-      dp_int(:)=0.
-      do l=2,nlay
-        do ig=1,ngrid
-           if(l.LE.lmax(ig)) THEN
-           zdp=pplay(ig,l-1)-pplay(ig,l)
-           alp_int(ig)=alp_int(ig)+0.5*rhobarz(ig,l)*wth3(ig,l)*zdp
-           dp_int(ig)=dp_int(ig)+zdp
-           endif
-        enddo
-      enddo
-
-      if (iflag_coupl>=3 .and. iflag_coupl<=5) then
-      do ig=1,ngrid
-!valeur integree de alp_bl * 0.5:
-        if (dp_int(ig)>0.) then
-        Alp_bl(ig)=alp_int(ig)/dp_int(ig)
-        endif
-      enddo!
-      endif
-
-
-! Facteur multiplicatif sur Alp_bl
-      Alp_bl(:)=alp_bl_k*Alp_bl(:)
-
-!------------------------------------------------------------------------
+     &                   )
+
 
 
