Changeset 2046 for LMDZ5/trunk

Timestamp:

May 18, 2014, 3:18:19 PM (11 years ago)

Author:

fhourdin

Message:

Nouvelles options des thermiques en route vers les stratocumulus

+ sorties de la discretisation verticale dans le vieux disvert_old

New options in thermals for stratocumulus

Location:

LMDZ5/trunk/libf/phylmd

Files:

• 1DUTILS.h (modified) (2 diffs)
• thermcell_plume.F90 (modified) (15 diffs)

LMDZ5/trunk/libf/phylmd/1DUTILS.h

@@ -822 +822 @@
       INTEGER np,ierr
       REAL pi,x
+      REAL :: scaleheight=8.
  
 !-----------------------------------------------------------------------
@@ -924 +925 @@
       PRINT *,  ap
  
-      DO l = 1, llm
-       dpres(l) = bp(l) - bp(l+1)
-       presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
-      ENDDO
+   print*,'scaleheight=',scaleheight
+  DO l = 1, llm
+     dpres(l) = bp(l) - bp(l+1)
+     presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
+     write(*, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
+          log(preff/presnivs(l))*scaleheight &
+          , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ &
+          max(ap(l+1)+bp(l+1)*preff, 1.e-10))
+  ENDDO
+
  
       PRINT *,' PRESNIVS '

LMDZ5/trunk/libf/phylmd/thermcell_plume.F90

@@ -6 +6 @@
      &           lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
      &           ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
-     &           ,lev_out,lunout1,igout) 
+     &           ,lev_out,lunout1,igout)
 !--------------------------------------------------------------------------
 !thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
@@ -80 +80 @@
       real zbuoyjam(ngrid,klev)
       real zbuoybis,zdz2,zdz3,lmel,entrbis,zdzbis
+      real fact_shell
+      real ztv1,ztv2,factinv,zinv,zlmel
+      real ztv_est1,ztv_est2
       real zcor,zdelta,zcvm5,qlbef
       real betalpha,zbetalpha
@@ -92 +95 @@
 ! Initialisation
 
-!     print*,'THERMCELL PLUME OK'
       RLvCp = RLVTT/RCPD
       fact_epsilon=0.002
       betalpha=0.9 
-      afact=2./3.            
+      afact=2./3.    
+      fact_shell=0.85       
 
       zbetalpha=betalpha/(1.+betalpha)
@@ -162 +165 @@
                lalim(ig)=l+1
                alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
-!              print*,'alim2',l,ztv(ig,l),ztv(ig,l+1),alim_star(ig,l)
+!               print*,'alim2',l,ztv(ig,l),ztv(ig,l+1),alim_star(ig,l)
             endif
          enddo
@@ -232 +235 @@
    call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:))
     do ig=1,ngrid 
-!      print*,'active',active(ig),ig,l
+!       print*,'active',active(ig),ig,l
         if(active(ig)) then 
         zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig))
@@ -258 +261 @@
 !    &                     Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)
 !              w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))*zdw2+w_est(ig,l)*exp(-zw2fact))
+
+!--------------------------------------------------
+!AJ052014: J'ai remplac? w_est(ig,l) par zw2(ig,l)
+!--------------------------------------------------
+         if (iflag_thermals_ed==8) then
+! Ancienne version
              w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
     &                     (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
     &                     (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2))
-             if (w_est(ig,l+1).lt.0.) then
+ 
+! Nouvelle version Arnaud
+         else
+              w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
+    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
+    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
+         endif
+
+
+         if (iflag_thermals_ed<6) then
+             zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l)
+!              fact_epsilon=0.0005/(zalpha+0.025)**0.5
+!              fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
+              fact_epsilon=0.0002/(zalpha+0.1)
+              zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
+              zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
+              zdw2=afact*zbuoy(ig,l)/fact_epsilon
+              zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon
+              w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
+    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
+    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
+
+         endif
+!--------------------------------------------------
+!AJ052014: J'ai comment? ce if plus n?cessaire puisqu'
+!on fait max(0.0001,.....)
+!--------------------------------------------------      
+
+!             if (w_est(ig,l+1).lt.0.) then
 !               w_est(ig,l+1)=zw2(ig,l)
-                w_est(ig,l+1)=0.0001
-             endif
+!                w_est(ig,l+1)=0.0001
+!             endif
+
        endif
     enddo
@@ -296 +334 @@
           
         lmel=fact_thermals_ed_dz*zlev(ig,l)
+        zlmel=zlev(ig,l)+lmel
 !        lmel=2.5*(zlev(ig,l)-zlev(ig,l-1))
         lt=l+1
-         do it=1,klev-(l+1) 
-          zdz2=zlev(ig,lt)-zlev(ig,l)
-          if (zdz2.gt.lmel) then 
-          zdz3=zlev(ig,lt)-zlev(ig,lt-1) 
+        zdz2=zlev(ig,lt)-zlev(ig,l)
+!--------------------------------------------------
+!AJ052014: J'ai remplac? la boucle do par un do while
+! afin de faire moins de calcul dans la boucle
+!--------------------------------------------------
+         do while (zdz2.lt.lmel)
+         lt=lt+1
+         zdz2=zlev(ig,lt)-zlev(ig,l)
+         end do
+
+         zdz3=zlev(ig,lt)-zlev(ig,lt-1) 
+
+!--------------------------------------------------
+!AJ052014: Si iflag_thermals_ed<8 (par ex 6), alors
+! on cherche o? se trouve l'altitude d'inversion 
+! en calculant ztv1 (interpolation de la valeur de 
+! theta au niveau lt en utilisant les niveaux lt-1 et
+! lt-2) et ztv2 (interpolation avec les niveaux lt+1
+! et lt+2). Si theta r?ellement calcul?e au niveau lt
+! comprise entre ztv1 et ztv2, alors il y a inversion
+! et on calcule son altitude zinv en supposant que ztv(lt)
+! est une combinaison lin?aire de ztv1 et ztv2.
+! Ensuite, on calcule la flottabilit? en comparant 
+! la temp?rature de la couche l ? celle de l'air situ? 
+! l+lmel plus haut, ce qui n?cessite de savoir quel fraction 
+! de cet air est au-dessus ou en-dessous de l'inversion   
+!--------------------------------------------------
+
+
+         if (iflag_thermals_ed.lt.8) then
+
+          ztv1=(ztv(ig,lt-1)-ztv(ig,lt-2))*zlev(ig,lt)/(zlev(ig,lt-1)-zlev(ig,lt-2)) &
+    &          +(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) &
+    &          /(zlev(ig,lt-1)-zlev(ig,lt-2))
+
+          ztv2=(ztv(ig,lt+2)-ztv(ig,lt+1))*zlev(ig,lt)/(zlev(ig,lt+2)-zlev(ig,lt+1)) &
+    &          +(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) &
+    &          /(zlev(ig,lt+2)-zlev(ig,lt+1))
+
+          if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then  
+
+          factinv=(ztv2-ztv(ig,lt))/(ztv2-ztv1)
+          zinv=zlev(ig,lt-1)+factinv*(zlev(ig,lt)-zlev(ig,lt-1))
+
+          if (zlmel+0.5*zdz.ge.zinv) then
+           if (zlmel-0.5*zdz.ge.zinv) then
+
+          ztv_est(ig,l)=(ztv(ig,lt+2)-ztv(ig,lt+1))*(zlmel-0.*zdz)/(zlev(ig,lt+2)-zlev(ig,lt+1)) &
+    &          +(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) &
+    &          /(zlev(ig,lt+2)-zlev(ig,lt+1))
+
+          zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l)+(1.-fact_shell)*zbuoy(ig,l)
+
+          else
+
+          ztv_est1=(ztv(ig,lt+2)-ztv(ig,lt+1))*0.5*(zlmel+zinv+0.5*zdz)/(zlev(ig,lt+2)-zlev(ig,lt+1)) &
+    &          +(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) &
+    &          /(zlev(ig,lt+2)-zlev(ig,lt+1))
+          ztv_est2=(ztv(ig,lt-1)-ztv(ig,lt-2))*0.5*(zinv+zlmel-0.5*zdz)/(zlev(ig,lt-1)-zlev(ig,lt-2)) &
+    &          +(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) &
+    &          /(zlev(ig,lt-1)-zlev(ig,lt-2))
+          zbuoyjam(ig,l)=fact_shell*RG*(((zlmel+0.5*zdz-zinv)/zdz)*(ztva_est(ig,l)- & 
+    &          ztv_est1)/ztv_est1+((zinv-zlmel+0.5*zdz)/zdz)*(ztva_est(ig,l)- &
+    &          ztv_est2)/ztv_est2)+(1.-fact_shell)*zbuoy(ig,l)
+
+           endif
+
+          else 
+
+          ztv_est(ig,l)=(ztv(ig,lt-1)-ztv(ig,lt-2))*(zlmel-0.*zdz)/(zlev(ig,lt-1)-zlev(ig,lt-2)) &
+    &          +(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) &
+    &          /(zlev(ig,lt-1)-zlev(ig,lt-2))
+
+          zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l)+(1.-fact_shell)*zbuoy(ig,l)
+!         ztv_est1=(ztv(ig,lt+2)-ztv(ig,lt+1))*0.5*(zlmel+zinv+0.5*zdz)/(zlev(ig,lt+2)-zlev(ig,lt+1)) &
+!    &          +(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) &
+!    &          /(zlev(ig,lt+2)-zlev(ig,lt+1))
+!          ztv_est2=(ztv(ig,lt-1)-ztv(ig,lt-2))*0.5*(zinv+zlmel-0.5*zdz)/(zlev(ig,lt-1)-zlev(ig,lt-2)) &
+!    &          +(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) &
+!    &          /(zlev(ig,lt-1)-zlev(ig,lt-2))
+!          zbuoyjam(ig,l)=fact_shell*RG*(((zlmel+0.5*zdz-zinv)/zdz)*(ztva_est(ig,l)- & 
+!    &          ztv_est1)/ztv_est1+((zinv-zlmel+0.5*zdz)/zdz)*(ztva_est(ig,l)- &
+!    &          ztv_est2)/ztv_est2)+(1.-fact_shell)*zbuoy(ig,l)
+
+
+
+
+!          print*,'on est pass? par l?',l,lt,zbuoyjam(ig,l),zbuoy(ig,l)
+          endif
+
+
+          else
+          
 !          ztv_est(ig,l)=(lmel/zdz2)*(ztv(ig,lt)-ztv(ig,l))+ztv(ig,l)
 !          zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l)
+
+         zbuoyjam(ig,l)=fact_shell*RG*(((lmel+zdz3-zdz2)/zdz3)*(ztva_est(ig,l)- & 
+    &          ztv(ig,lt))/ztv(ig,lt)+((zdz2-lmel)/zdz3)*(ztva_est(ig,l)- &
+    &          ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l)
+
+!         zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- & 
+!    &          po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- &
+!     &          po(ig,lt-1))/po(ig,lt-1))
+
+          endif
+
+          else
 
          zbuoyjam(ig,l)=1.*RG*(((lmel+zdz3-zdz2)/zdz3)*(ztva_est(ig,l)- & 
@@ -309 +449 @@
     &          ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l)
 
-!         zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- & 
-!    &          po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- &
-!    &          po(ig,lt-1))/po(ig,lt-1))
+          endif
           
-          else
-          lt=lt+1
-          endif
-          enddo
 
 !          zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
@@ -325 +459 @@
 !          entrbis=entr_star(ig,l)
 
-
-          detr_star(ig,l)=f_star(ig,l)*zdz                        &
+          if (iflag_thermals_ed.lt.6) then
+          fact_epsilon=0.0002/(zalpha+0.1)
+          endif
+
+          detr_star(ig,l)=f_star(ig,l)*zdz             &
     &     *MAX(1.e-4, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m          &
-    &     + 0.012*(zdqt(ig,l)/zw2m)**0.5 )
+    &     + 0.012*(zdqt(ig,l)/zw2m)**0.5)
 
           zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
 
           entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0.,  &
-    &     afact*zbuoy(ig,l)/zw2m - fact_epsilon )
-!   &     afact*zbuoy(ig,l)/zw2m - fact_epsilon+ 0.012*(zdqt(ig,l)/zw2m)**0.5)
-
+    &     afact*zbuoy(ig,l)/zw2m - fact_epsilon)
 
 !          entr_star(ig,l)=Max(0.,f_star(ig,l)*zdz*zbetalpha*  &     
@@ -352 +487 @@
 !        endif
 
-!print*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l)
+!       print*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l)
 ! Calcul du flux montant normalise
       f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
@@ -395 +530 @@
            zdzbis=zlev(ig,l+1)-zlev(ig,l-1)
            zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
-
+           fact_epsilon=0.002
             zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
             zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
@@ -404 +539 @@
     &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
     &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
+
+           if (iflag_thermals_ed.lt.6) then
+           zalpha=f0(ig)*f_star(ig,l)/sqrt(zw2(ig,l+1))/rhobarz(ig,l)
+!           fact_epsilon=0.0005/(zalpha+0.025)**0.5
+!           fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5)
+           fact_epsilon=0.0002/(zalpha+0.1)**1
+            zw2fact=fact_epsilon*2.*zdz/(1.+betalpha)
+            zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha)
+            zdw2= afact*zbuoy(ig,l)/(fact_epsilon)
+            zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon)
+
+            zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* &
+    &                     (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* &
+    &                     (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2))
+
+           endif
+
+
       endif
    enddo
@@ -470 +623 @@
         if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
 
+#undef wrgrads_thermcell
+#ifdef wrgrads_thermcell
+         call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta      ','esta      ')
+         call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta      ','dsta      ')
+         call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy      ','buoy      ')
+         call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt      ','dqt      ')
+         call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est     ','w_est     ')
+         call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2     ','w_es2     ')
+         call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A      ','zw2A      ')
+#endif
+
+
      return 
      end
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 
 
@@ -831 +1027 @@
                 zw2(ig,l+1)=0.
                 linter(ig)=l+1
-!CR:04/05/12:calcul linter
-        elseif (f_star(ig,l+1).lt.0.) then
-           linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l))  &
-     &               -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l))
-           zw2(ig,l+1)=0.
-!fin CR:04/05/12
-        endif
-
+            endif
 
         if (zw2(ig,l+1).lt.0.) then 
@@ -881 +1070 @@
         if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
 
+#undef wrgrads_thermcell
+#ifdef wrgrads_thermcell
+         call wrgradsfi(1,klev,entr_star(igout,1:klev),'esta      ','esta      ')
+         call wrgradsfi(1,klev,detr_star(igout,1:klev),'dsta      ','dsta      ')
+         call wrgradsfi(1,klev,zbuoy(igout,1:klev),'buoy      ','buoy      ')
+         call wrgradsfi(1,klev,zdqt(igout,1:klev),'dqt      ','dqt      ')
+         call wrgradsfi(1,klev,w_est(igout,1:klev),'w_est     ','w_est     ')
+         call wrgradsfi(1,klev,w_est(igout,2:klev+1),'w_es2     ','w_es2     ')
+         call wrgradsfi(1,klev,zw2(igout,1:klev),'zw2A      ','zw2A      ')
+#endif