Changeset 512

Timestamp:

Feb 3, 2012, 11:11:45 AM (13 years ago)

Author:

acolaitis

Message:

Updated thermals version, that produces quite good overshoots in the inversion region. Series of zmax matches very well the LES and underline one of the weaknesses of convadj. This version will change, as this inversion is only produced in 222 levels. Some resolution tricks must be found for it to work at 32 levels.

Location:

trunk/LMDZ.MARS/libf/phymars

Files:

• calltherm_interface.F90 (modified) (4 diffs)
• thermcell_main_mars.F90 (modified) (15 diffs)

trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90

@@ -59 +59 @@
       REAL detr_therm(ngridmx,nlayermx)
       REAL zw2(ngridmx,nlayermx+1)
-      REAL fraca(ngridmx,nlayermx+1)
+      REAL fraca(ngridmx,nlayermx+1),zfraca(ngridmx,nlayermx+1)
       REAL ztla(ngridmx,nlayermx)
       REAL q_therm(ngridmx,nlayermx), pq_therm(ngridmx,nlayermx,nqmx)
@@ -124 +124 @@
       zw2(:,:)=0.
       fraca(:,:)=0.
+      zfraca(:,:)=0.
       if (tracer) then
          pdq_th(:,:,:)=0.
@@ -264 +265 @@
             detr_therm(:,:)=detr_therm(:,:)  &
      &       +zdetr_therm(:,:)*fact
+            zfraca(:,:)=zfraca(:,:) + fraca(:,:)*fact
 
             heatFlux(:,:)=heatFlux(:,:) &
@@ -273 +275 @@
             buoyancyEst(:,:)=buoyancyEst(:,:) &
      &       +zbuoyancyEst(:,:)*fact
+  
 
             zw2(:,:)=zw2(:,:) + zzw2(:,:)*fact

trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90

@@ -113 +113 @@
       REAL zdz,zbuoy(ngridmx,nlayermx),zw2m
       LOGICAL activecell(ngridmx),activetmp(ngridmx)
-      REAL a1,b1,ae,be,ad,bd,fdfu
+      REAL a1,b1,ae,be,ad,bd,fdfu,b1inv,a1inv
       INTEGER tic
 
@@ -379 +379 @@
       fdfu = -1.3
 
+! Best configuration for 222 levels:
+      b1=0.
+      a1=0.7*a1
+      a1inv=0.25*a1
+      b1inv=0.0002
+
 ! --------------------------------------------------------------------------
 ! --------------------------------------------------------------------------
@@ -406 +412 @@
 !        do l=2,4
          do ig=1,ngridmx
-           if (ztv(ig,l)>(ztv(ig,l+1)+0.) .and. ztv(ig,1)>=ztv(ig,l) .and. (alim_star(ig,l-1) .ne. 0.)) then
+           if (ztv(ig,l)>(ztv(ig,l+1)+0.) .and. ztv(ig,1)>=ztv(ig,l) .and. (alim_star(ig,l-1) .ne. 0.) .and. (zlev(ig,l+1) .lt. 1000.)) then
                alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
      &                       *sqrt(zlev(ig,l+1))
@@ -499 +505 @@
      & -2.*zdz*w_est(ig,l)*ae*(a1*zbuoy(ig,l)/w_est(ig,l)-b1)**be)
                 else
-                w_est(ig,l+1)=Max(0.0001,w_est(ig,l)+2.*zdz*a1*zbuoy(ig,l)-2.*zdz*w_est(ig,l)*b1)
+                w_est(ig,l+1)=Max(0.0001,w_est(ig,l)+2.*zdz*a1inv*zbuoy(ig,l)-2.*zdz*w_est(ig,l)*b1inv)
                 endif
                 if (w_est(ig,l+1).lt.0.) then
@@ -518 +524 @@
           if((a1*(zbuoy(ig,l)/zw2m)-b1) .gt. 0.) then
           entr_star(ig,l)=f_star(ig,l)*zdz*  &
-!        &   MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be)
-         &  MAX(0.,log(1. + 0.027*sqrt(a1*(zbuoy(ig,l)/zw2m)-b1)))
+        &   MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be)
+!         &  MAX(0.,log(1. + 0.03*sqrt(a1*(zbuoy(ig,l)/zw2m)-b1)))
           else
           entr_star(ig,l)=0.
@@ -540 +546 @@
 ! new param from continuity eq with a fit on dfdz
                  detr_star(ig,l) = f_star(ig,l)*zdz*              &
-!            &  ad
-             &  Max(0., 0.001 - 0.45*zbuoy(ig,l)/zw2m)
+            &  ad
+!             &  Max(0., 0.0005 - 0.55*zbuoy(ig,l)/zw2m)
 
 !               &     MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005)   !svn baseline
@@ -553 +559 @@
           else
           detr_star(ig,l)=f_star(ig,l)*zdz*                        &
-!                &     bd*zbuoy(ig,l)/zw2m
-             &  Max(0., 0.001 - 0.45*zbuoy(ig,l)/zw2m)
+                &     bd*zbuoy(ig,l)/zw2m
+!             &  Max(0., 0.001 - 0.45*zbuoy(ig,l)/zw2m)
+!             &  Max(0., Min(0.001,0.0005 - 0.55*zbuoy(ig,l)/zw2m))
 
 
@@ -573 +580 @@
 ! alim_star et 0 sinon
 
-        if (l.lt.lalim(ig)) then
+       if (l.lt.lalim(ig)) then
           alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
-          entr_star(ig,l)=0.
-        endif
+          entr_star(ig,l)=0. 
+       endif
 
 ! Calcul du flux montant normalise
@@ -591 +598 @@
 !---------------------------------------------------------------------------
 
-      DO tic=0,1  ! internal convergence loop
+      DO tic=0,0  ! internal convergence loop
       activetmp(:)=activecell(:) .and. f_star(:,l+1)>1.e-10
       do ig=1,ngridmx
@@ -614 +621 @@
      & 2.*zdz*zw2(ig,l)*b1-2.*zdz*zw2(ig,l)*ae*(a1*zbuoy(ig,l)/zw2(ig,l)-b1)**be)
            else
-           zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.*zdz*a1*zbuoy(ig,l)-2.*zdz*zw2(ig,l)*b1)
+           zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.*zdz*a1inv*zbuoy(ig,l)-2.*zdz*zw2(ig,l)*b1inv)
            endif
       endif
@@ -630 +637 @@
           entr_star(ig,l)=f_star(ig,l)*zdz*  &
         &   MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be)
+!         &  MAX(0.,log(1. + 0.03*sqrt(a1*(zbuoy(ig,l)/zw2m)-b1)))
           else
           entr_star(ig,l)=0.
@@ -640 +648 @@
                  detr_star(ig,l) = f_star(ig,l)*zdz*              &
             &  ad
+!             &  Max(0., 0.0005 - 0.55*zbuoy(ig,l)/zw2m)
+
              endif
           else
           detr_star(ig,l)=f_star(ig,l)*zdz*                        &
                 &     bd*zbuoy(ig,l)/zw2m
+!             &  Max(0.,Min(0.001,0.0005 - 0.55*zbuoy(ig,l)/zw2m))
+
           endif
           else
@@ -736 +748 @@
          do l=nlayermx,lalim(ig)+1,-1
             if (zw2(ig,l).le.1.e-10) then
-               lmax(ig)=l-1
+               lmax(ig)=l-1 
             endif
          enddo
@@ -1165 +1177 @@
          if (lmax(ig) .gt. 1) then
          do l=1,lmax(ig)
-              if(zlay(ig,l) .le. 0.7*zmax(ig)) then
+!              if(zlay(ig,l) .le. 0.8*zmax(ig)) then
+              if(zlay(ig,l) .le. zmax(ig)) then
               fm_down(ig,l) =fm(ig,l)* &
      &      max(fdfu,-3*max(0.,(zlay(ig,l)/zmax(ig)))-0.9)
               endif
 
-             if(zlay(ig,l) .le. 0.06*zmax(ig)) then
-          ztvd(ig,l)=ztv(ig,l)*max(0.,(1.+(sqrt((zlay(ig,l)/zmax(ig))/0.122449) - 1.)*(ztva(ig,l)/ztv(ig,l) - 1.)))
-             elseif(zlay(ig,l) .le. 0.4*zmax(ig)) then
-          ztvd(ig,l)=ztv(ig,l)*max(0.,1.-0.25*(ztva(ig,l)/ztv(ig,l) - 1.))
-             elseif(zlay(ig,l) .le. 0.7*zmax(ig)) then
-          ztvd(ig,l)=ztv(ig,l)*max(0.,(1.+(((zlay(ig,l)/zmax(ig))-0.7)/1.)*(ztva(ig,l)/ztv(ig,l) - 1.)))
+!             if(zlay(ig,l) .le. 0.06*zmax(ig)) then
+!          ztvd(ig,l)=ztv(ig,l)*max(0.,(1.+(sqrt((zlay(ig,l)/zmax(ig))/0.122449) - 1.)*(ztva(ig,l)/ztv(ig,l) - 1.)))
+!             elseif(zlay(ig,l) .le. 0.4*zmax(ig)) then
+!          ztvd(ig,l)=ztv(ig,l)*max(0.,1.-0.25*(ztva(ig,l)/ztv(ig,l) - 1.))
+!             elseif(zlay(ig,l) .le. 0.7*zmax(ig)) then
+!          ztvd(ig,l)=ztv(ig,l)*max(0.,(1.+(((zlay(ig,l)/zmax(ig))-0.7)/1.)*(ztva(ig,l)/ztv(ig,l) - 1.)))
+!             else
+!          ztvd(ig,l)=ztv(ig,l)
+!            endif
+
+!            if(zlay(ig,l) .le. 0.6*zmax(ig)) then
+!          ztvd(ig,l)=ztv(ig,l)*((zlay(ig,l)/zmax(ig))/179.848 + 0.997832)
+!            elseif(zlay(ig,l) .le. 0.8*zmax(ig)) then
+!          ztvd(ig,l)=-ztv(ig,l)*(((zlay(ig,l)/zmax(ig))-171.74)/170.94)
+!             else
+!          ztvd(ig,l)=ztv(ig,l)
+!            endif
+
+
+!            if(zlay(ig,l) .le. 0.8*zmax(ig)) then
+!          ztvd(ig,l)=ztv(ig,l)*((zlay(ig,l)/zmax(ig))/224.81 + 0.997832)
+!            elseif(zlay(ig,l) .le. zmax(ig)) then
+!          ztvd(ig,l)=-ztv(ig,l)*(((zlay(ig,l)/zmax(ig))-144.685)/143.885)
+!             else
+!          ztvd(ig,l)=ztv(ig,l)
+!            endif
+
+
+            if(zlay(ig,l) .le. zmax(ig)) then
+          ztvd(ig,l)=ztv(ig,l)*((zlay(ig,l)/zmax(ig))/299.7 + 0.997832)
              else
           ztvd(ig,l)=ztv(ig,l)
             endif
-
 
          enddo
@@ -1432 +1468 @@
          do ig=1,ngridmx
          pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))*zpopsk(ig,l)*ratiom(ig,l)
+!         pdtadj(ig,l)=zdthladj(ig,l)*zpopsk(ig,l)*ratiom(ig,l)
          enddo
       enddo