Ignore:
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.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r508 r512  
    113113      REAL zdz,zbuoy(ngridmx,nlayermx),zw2m
    114114      LOGICAL activecell(ngridmx),activetmp(ngridmx)
    115       REAL a1,b1,ae,be,ad,bd,fdfu
     115      REAL a1,b1,ae,be,ad,bd,fdfu,b1inv,a1inv
    116116      INTEGER tic
    117117
    … …  
    379379      fdfu = -1.3
    380380
     381! Best configuration for 222 levels:
     382      b1=0.
     383      a1=0.7*a1
     384      a1inv=0.25*a1
     385      b1inv=0.0002
     386
    381387! --------------------------------------------------------------------------
    382388! --------------------------------------------------------------------------
    … …  
    406412!        do l=2,4
    407413         do ig=1,ngridmx
    408            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
     414           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
    409415               alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
    410416     &                       *sqrt(zlev(ig,l+1))
    … …  
    499505     & -2.*zdz*w_est(ig,l)*ae*(a1*zbuoy(ig,l)/w_est(ig,l)-b1)**be)
    500506                else
    501                 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)
     507                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)
    502508                endif
    503509                if (w_est(ig,l+1).lt.0.) then
    … …  
    518524          if((a1*(zbuoy(ig,l)/zw2m)-b1) .gt. 0.) then
    519525          entr_star(ig,l)=f_star(ig,l)*zdz*  &
    520 !        &   MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be)
    521          &  MAX(0.,log(1. + 0.027*sqrt(a1*(zbuoy(ig,l)/zw2m)-b1)))
     526        &   MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be)
     527!         &  MAX(0.,log(1. + 0.03*sqrt(a1*(zbuoy(ig,l)/zw2m)-b1)))
    522528          else
    523529          entr_star(ig,l)=0.
    … …  
    540546! new param from continuity eq with a fit on dfdz
    541547                 detr_star(ig,l) = f_star(ig,l)*zdz*              &
    542 !            &  ad
    543              &  Max(0., 0.001 - 0.45*zbuoy(ig,l)/zw2m)
     548            &  ad
     549!             &  Max(0., 0.0005 - 0.55*zbuoy(ig,l)/zw2m)
    544550
    545551!               &     MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005)   !svn baseline
    … …  
    553559          else
    554560          detr_star(ig,l)=f_star(ig,l)*zdz*                        &
    555 !                &     bd*zbuoy(ig,l)/zw2m
    556              &  Max(0., 0.001 - 0.45*zbuoy(ig,l)/zw2m)
     561                &     bd*zbuoy(ig,l)/zw2m
     562!             &  Max(0., 0.001 - 0.45*zbuoy(ig,l)/zw2m)
     563!             &  Max(0., Min(0.001,0.0005 - 0.55*zbuoy(ig,l)/zw2m))
    557564
    558565
    … …  
    573580! alim_star et 0 sinon
    574581
    575         if (l.lt.lalim(ig)) then
     582       if (l.lt.lalim(ig)) then
    576583          alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
    577           entr_star(ig,l)=0.
    578         endif
     584          entr_star(ig,l)=0. 
     585       endif
    579586
    580587! Calcul du flux montant normalise
    … …  
    591598!---------------------------------------------------------------------------
    592599
    593       DO tic=0,1  ! internal convergence loop
     600      DO tic=0,0  ! internal convergence loop
    594601      activetmp(:)=activecell(:) .and. f_star(:,l+1)>1.e-10
    595602      do ig=1,ngridmx
    … …  
    614621     & 2.*zdz*zw2(ig,l)*b1-2.*zdz*zw2(ig,l)*ae*(a1*zbuoy(ig,l)/zw2(ig,l)-b1)**be)
    615622           else
    616            zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.*zdz*a1*zbuoy(ig,l)-2.*zdz*zw2(ig,l)*b1)
     623           zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.*zdz*a1inv*zbuoy(ig,l)-2.*zdz*zw2(ig,l)*b1inv)
    617624           endif
    618625      endif
    … …  
    630637          entr_star(ig,l)=f_star(ig,l)*zdz*  &
    631638        &   MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be)
     639!         &  MAX(0.,log(1. + 0.03*sqrt(a1*(zbuoy(ig,l)/zw2m)-b1)))
    632640          else
    633641          entr_star(ig,l)=0.
    … …  
    640648                 detr_star(ig,l) = f_star(ig,l)*zdz*              &
    641649            &  ad
     650!             &  Max(0., 0.0005 - 0.55*zbuoy(ig,l)/zw2m)
     651
    642652             endif
    643653          else
    644654          detr_star(ig,l)=f_star(ig,l)*zdz*                        &
    645655                &     bd*zbuoy(ig,l)/zw2m
     656!             &  Max(0.,Min(0.001,0.0005 - 0.55*zbuoy(ig,l)/zw2m))
     657
    646658          endif
    647659          else
    … …  
    736748         do l=nlayermx,lalim(ig)+1,-1
    737749            if (zw2(ig,l).le.1.e-10) then
    738                lmax(ig)=l-1
     750               lmax(ig)=l-1 
    739751            endif
    740752         enddo
    … …  
    11651177         if (lmax(ig) .gt. 1) then
    11661178         do l=1,lmax(ig)
    1167               if(zlay(ig,l) .le. 0.7*zmax(ig)) then
     1179!              if(zlay(ig,l) .le. 0.8*zmax(ig)) then
     1180              if(zlay(ig,l) .le. zmax(ig)) then
    11681181              fm_down(ig,l) =fm(ig,l)* &
    11691182     &      max(fdfu,-3*max(0.,(zlay(ig,l)/zmax(ig)))-0.9)
    11701183              endif
    11711184
    1172              if(zlay(ig,l) .le. 0.06*zmax(ig)) then
    1173           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.)))
    1174              elseif(zlay(ig,l) .le. 0.4*zmax(ig)) then
    1175           ztvd(ig,l)=ztv(ig,l)*max(0.,1.-0.25*(ztva(ig,l)/ztv(ig,l) - 1.))
    1176              elseif(zlay(ig,l) .le. 0.7*zmax(ig)) then
    1177           ztvd(ig,l)=ztv(ig,l)*max(0.,(1.+(((zlay(ig,l)/zmax(ig))-0.7)/1.)*(ztva(ig,l)/ztv(ig,l) - 1.)))
     1185!             if(zlay(ig,l) .le. 0.06*zmax(ig)) then
     1186!          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.)))
     1187!             elseif(zlay(ig,l) .le. 0.4*zmax(ig)) then
     1188!          ztvd(ig,l)=ztv(ig,l)*max(0.,1.-0.25*(ztva(ig,l)/ztv(ig,l) - 1.))
     1189!             elseif(zlay(ig,l) .le. 0.7*zmax(ig)) then
     1190!          ztvd(ig,l)=ztv(ig,l)*max(0.,(1.+(((zlay(ig,l)/zmax(ig))-0.7)/1.)*(ztva(ig,l)/ztv(ig,l) - 1.)))
     1191!             else
     1192!          ztvd(ig,l)=ztv(ig,l)
     1193!            endif
     1194
     1195!            if(zlay(ig,l) .le. 0.6*zmax(ig)) then
     1196!          ztvd(ig,l)=ztv(ig,l)*((zlay(ig,l)/zmax(ig))/179.848 + 0.997832)
     1197!            elseif(zlay(ig,l) .le. 0.8*zmax(ig)) then
     1198!          ztvd(ig,l)=-ztv(ig,l)*(((zlay(ig,l)/zmax(ig))-171.74)/170.94)
     1199!             else
     1200!          ztvd(ig,l)=ztv(ig,l)
     1201!            endif
     1202
     1203
     1204!            if(zlay(ig,l) .le. 0.8*zmax(ig)) then
     1205!          ztvd(ig,l)=ztv(ig,l)*((zlay(ig,l)/zmax(ig))/224.81 + 0.997832)
     1206!            elseif(zlay(ig,l) .le. zmax(ig)) then
     1207!          ztvd(ig,l)=-ztv(ig,l)*(((zlay(ig,l)/zmax(ig))-144.685)/143.885)
     1208!             else
     1209!          ztvd(ig,l)=ztv(ig,l)
     1210!            endif
     1211
     1212
     1213            if(zlay(ig,l) .le. zmax(ig)) then
     1214          ztvd(ig,l)=ztv(ig,l)*((zlay(ig,l)/zmax(ig))/299.7 + 0.997832)
    11781215             else
    11791216          ztvd(ig,l)=ztv(ig,l)
    11801217            endif
    1181 
    11821218
    11831219         enddo
    … …  
    14321468         do ig=1,ngridmx
    14331469         pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))*zpopsk(ig,l)*ratiom(ig,l)
     1470!         pdtadj(ig,l)=zdthladj(ig,l)*zpopsk(ig,l)*ratiom(ig,l)
    14341471         enddo
    14351472      enddo
Note: See TracChangeset for help on using the changeset viewer.