Changeset 496


Ignore:
Timestamp:
Jan 11, 2012, 2:21:08 PM (13 years ago)
Author:
acolaitis
Message:

M 489 libf/phymars/thermcell_main_mars.F90
--------------------> New parameters for thermals, following latest version of the relevant article

M 489 libf/phymars/physiq.F
--------------------> Minor modifications

D 489 libf/phymars/surflayer_interpol.F
A 0 libf/phymars/pbl_parameters.F
--------------------> Replaced surflayer_interpol.F by a cleaner and richer version called pbl_parameters.F

This routine is the base of what will later be implemented in the MCD and as a tool to
compute surface properties from output fields (so that it will ultimately disappear from libf
and might become an utility)

M 489 libf/phymars/vdif_cd.F
--------------------> Minor modification to the Richardson computation

M 489 libf/phymars/vdifc.F
--------------------> Minor modification to the aerodynamic conductances computation, replacing wmax by hfmax, which

is a better proxy for convective activity. Might be replaced by w_star later.

Location:
trunk/LMDZ.MARS
Files:
1 added
1 deleted
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/README

    r485 r496  
    13611361.true. will call functions geticecover to yield a nice evolving CO2 cap as measured by Titus (mostly useful for mesoscale)
    13621362.false. is the default option, so this change is transparent to the casual user
     1363
     1364== 11/01/2012 == AC
     1365>> Corrections to the thermals scheme following latest revisions of the related paper. Replaced the surface layer interpolator by a more
     1366complete routine, that will ultimately become a post-processing utility (and disappear from libf) and a subroutine in the MCD.

  • trunk/LMDZ.MARS/libf/phymars/physiq.F

    r485 r496  
    333333      REAL lmax_th_out(ngridmx),zmax_th(ngridmx)
    334334      REAL, SAVE :: wmax_th(ngridmx)
    335       REAL hfmax_th(ngridmx)
     335      REAL, SAVE :: hfmax_th(ngridmx)
    336336      REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx)
    337337      REAL pdt_th(ngridmx,nlayermx),pdq_th(ngridmx,nlayermx,nqmx)
    … …  
    342342      REAL z_out                          ! height of interpolation between z0 and z1 [meters]
    343343      REAL ustar(ngridmx),tstar(ngridmx)  ! friction velocity and friction potential temp
    344       REAL L_mo(ngridmx)
     344      REAL L_mo(ngridmx),wstar(ngridmx),vhf(ngridmx),vvv(ngridmx)
    345345      REAL zu2(ngridmx)
    346346c=======================================================================
    … …  
    569569                 enddo
    570570              endif
     571           else
     572              zdtnlte(:,:)=0.
    571573           endif
    572574
    … …  
    746748             IF (zh(ig,1) .lt. tsurf(ig)) THEN
    747749               wmax_th(ig)=1.
    748              ENDIF       
     750               hfmax_th(ig)=0.2
     751             ELSE
     752               wmax_th(ig)=0.
     753               hfmax_th(ig)=0.
     754             ENDIF     
    749755          ENDDO
    750756        ENDIF
    … …  
    762768     $        zdum1,zdum2,zdh,pdq,zflubid,
    763769     $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
    764      &        zdqdif,zdqsdif,wmax_th,zcdv,zcdh)
     770     &        zdqdif,zdqsdif,wmax_th,zcdv,zcdh,hfmax_th)
    765771
    766772#ifdef MESOSCALE
    … …  
    853859        lmax_th(:)=0
    854860        wmax_th(:)=0.
     861        hfmax_th(:)=0.
    855862        lmax_th_out(:)=0.
    856863        end if
    … …  
    18031810         z_out=0.
    18041811         if (calltherm .and. (z_out .gt. 0.)) then
    1805          call surflayer_interpol(ngrid,nlayer,z0,g,zzlay,zu,zv,wmax_th
    1806      &              ,tsurf,zh,z_out,Teta_out,u_out,ustar,tstar,L_mo)
    1807 
    1808          zu2(:)=sqrt(zu(:,1)*zu(:,1)+zv(:,1)*zv(:,1))
    1809          call WRITEDIAGFI(ngridmx,'sqrt(zu2)',
    1810      &              'horizontal velocity norm','m/s',
    1811      &                         2,zu2)
    1812 
    1813          call WRITEDIAGFI(ngridmx,'Teta_out',
    1814      &              'potential temperature at z_out','K',
    1815      &                         2,Teta_out)
    1816          call WRITEDIAGFI(ngridmx,'u_out',
    1817      &              'horizontal velocity norm at z_out','m/s',
    1818      &                         2,u_out)
    1819          call WRITEDIAGFI(ngridmx,'u*',
    1820      &              'friction velocity','m/s',
    1821      &                         2,ustar)
    1822          call WRITEDIAGFI(ngridmx,'teta*',
    1823      &              'friction potential temperature','K',
    1824      &                         2,tstar)
    1825          call WRITEDIAGFI(ngrid,'L',
    1826      &              'Monin Obukhov length','m',
    1827      &                         2,L_mo)
     1812
     1813         call pbl_parameters(ngrid,nlayer,z0,
     1814     & g,zzlay,zu,zv,wmax_th,hfmax_th,zmax_th,tsurf,zh,z_out,
     1815     & Teta_out,u_out,ustar,tstar,wstar,L_mo,vhf,vvv)
     1816
    18281817         else
    18291818           if((.not. calltherm).and.(z_out .gt. 0.)) then
    … …  
    19141903         z_out=0.
    19151904         if (calltherm .and. (z_out .gt. 0.)) then
    1916          call surflayer_interpol(ngrid,nlayer,z0,g,zzlay,zu,zv,wmax_th
    1917      &              ,tsurf,zh,z_out,Teta_out,u_out,ustar,tstar,L_mo)
    1918 
    1919          zu2(:)=sqrt(zu(:,1)*zu(:,1)+zv(:,1)*zv(:,1))
    1920          call WRITEDIAGFI(ngridmx,'sqrt(zu2)',
    1921      &              'horizontal velocity norm','m/s',
    1922      &                         0,zu2)
    1923 
    1924          call WRITEDIAGFI(ngridmx,'Teta_out',
    1925      &              'potential temperature at z_out','K',
    1926      &                         0,Teta_out)
    1927          call WRITEDIAGFI(ngridmx,'u_out',
    1928      &              'horizontal velocity norm at z_out','m/s',
    1929      &                         0,u_out)
    1930          call WRITEDIAGFI(ngridmx,'u*',
    1931      &              'friction velocity','m/s',
    1932      &                         0,ustar)
    1933          call WRITEDIAGFI(ngridmx,'teta*',
    1934      &              'friction potential temperature','K',
    1935      &                         0,tstar)
     1905
     1906         call pbl_parameters(ngrid,nlayer,z0,
     1907     & g,zzlay,zu,zv,wmax_th,hfmax_th,zmax_th,tsurf,zh,z_out,
     1908     & Teta_out,u_out,ustar,tstar,wstar,L_mo,vhf,vvv)
     1909
    19361910         else
    19371911           if((.not. calltherm).and.(z_out .gt. 0.)) then
    … …  
    19471921     &              'hauteur du thermique','point',
    19481922     &                         0,lmax_th_out)
     1923        call WRITEDIAGFI(ngridmx,'zmax_th',
     1924     &              'hauteur du thermique','m',
     1925     &                         0,zmax_th)
    19491926        call WRITEDIAGFI(ngridmx,'hfmax_th',
    19501927     &              'maximum TH heat flux','K.m/s',

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

    r356 r496  
    111111      REAL zdz,zbuoy(ngridmx,nlayermx),zw2m
    112112      LOGICAL active(ngridmx),activetmp(ngridmx)
    113       REAL a1,b1,ae,be,ad,bd
     113      REAL a1,b1,ae,be,ad,bd,fdfu
    114114      INTEGER tic
    115115
    … …  
    304304      linter(:)=1.
    305305
     306! --------------------------------------------------------------------------
     307! -------------- MAIN PARAMETERS FOR THERMALS MODEL ------------------------
     308! --------------  see thermiques.pro and getfit.py -------------------------
     309
    306310!      a1=2.5 ; b1=0.0015 ; ae=0.045 ; be = 0.6  ! svn baseline
    307311
    308       a1=1.60226 ; b1=0.0006 ; ae=0.0454 ; be = 0.57 !improved fits
    309       ad = 0.0005114  ; bd = -0.662
     312! Using broad downdraft selection
     313!      a1=1.60226 ; b1=0.0006 ; ae=0.0454 ; be = 0.57
     314!      ad = 0.0005114  ; bd = -0.662
     315!      fdfu = -1.9
     316
     317! Using conditional sampling downdraft selection
     318      a1=1.4716 ; b1=0.0005698 ; ae=0.03683 ; be = 0.57421
     319      ad = 0.00048088  ; bd = -0.6697
     320      fdfu = -1.3
     321
     322! --------------------------------------------------------------------------
     323! --------------------------------------------------------------------------
     324! --------------------------------------------------------------------------
    310325
    311326! Initialisation des variables entieres
    … …  
    513528!---------------------------------------------------------------------------
    514529
    515       DO tic=0,6  ! internal convergence loop
     530      DO tic=0,1  ! internal convergence loop
    516531      activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
    517532      do ig=1,ngridmx
    … …  
    640655! ================= FIN PLUME ===============================================
    641656! ===========================================================================
    642 
    643657
    644658! ===========================================================================
    … …  
    10891103         do l=1,lmax(ig)
    10901104              if(zlay(ig,l) .le. 0.7*zmax(ig)) then
    1091               fm_down(ig,l) =-1.9*fm(ig,l)
     1105              fm_down(ig,l) =fm(ig,l)* &
     1106     &      max(fdfu,-3*max(0.,(zlay(ig,l)/zmax(ig)))-0.9)
    10921107              endif
    10931108
    … …  
    11021117            endif
    11031118
    1104 !             if(zlay(ig,l) .le. 0.06*zmax(ig)) then
    1105 !          ztvd(ig,l)=ztv(ig,l)*max(0.,(zlay(ig,l)/zmax(ig))/19.231 + 0.9938)
    1106 !             elseif(zlay(ig,l) .le. 0.6*zmax(ig)) then
    1107 !          ztvd(ig,l)=ztv(ig,l)*max(0.,(zlay(ig,l)/zmax(ig)-0.075)/187.931 + 0.9982)
    1108 !             else
    1109 !!             elseif(zlay(ig,l) .le. 0.7*zmax(ig)) then
    1110 !          ztvd(ig,l)=ztv(ig,l)*max(0.,(zlay(ig,l)/zmax(ig) -0.60)/(-1333) + 1.00025)
    1111 !!             else
    1112 !!          ztvd(ig,l)=ztv(ig,l)
    1113 !            endif
    11141119
    11151120         enddo
    … …  
    11261131         if(lmax(ig) .gt. 1) then
    11271132! No downdraft in the very-near surface layer, we begin at k=3
    1128          do k=2,lmax(ig)
     1133 
     1134         do k=3,lmax(ig)
    11291135            zdthladj_down(ig,k)=(1./masse(ig,k))*(fm_down(ig,k+1)*ztv(ig,k+1)- &
    11301136     & fm_down(ig,k)*ztv(ig,k)+fm_down(ig,k)*ztvd(ig,k)-fm_down(ig,k+1)*ztvd(ig,k+1))
    … …  
    13321338         do ig=1,ngridmx
    13331339           pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))*zpopsk(ig,l)
     1340!           pdtadj(ig,l)=zdthladj(ig,l)*zpopsk(ig,l)
    13341341         enddo
    13351342      enddo

  • trunk/LMDZ.MARS/libf/phymars/vdif_cd.F

    r329 r496  
    148148!         zu2=MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1),0.25*wmax(ig)**2)
    149149!         zu2=pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1)
    150 !         zu2=MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1),wmax(ig)**2)
    151        zu2(ig)=pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) + (0.5*wmax(ig))**2
     150         zu2(ig)=MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1),wmax(ig)**2)
     151!       zu2(ig)=pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) + (0.5*wmax(ig))**2
    152152
    153153               ! we add the wmax to simulate

  • trunk/LMDZ.MARS/libf/phymars/vdifc.F

    r473 r496  
    55     $                pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
    66     $                pdudif,pdvdif,pdhdif,pdtsrf,pq2,
    7      $                pdqdif,pdqsdif,wmax,zcdv_true,zcdh_true)
     7     $                pdqdif,pdqsdif,wmax,zcdv_true,zcdh_true,hfmax)
    88      IMPLICIT NONE
    99
    … …  
    4444      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
    4545      REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
    46       REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
     46      REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay),pt(ngrid,nlay)
    4747      REAL ptsrf(ngrid),pemis(ngrid)
    4848      REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay)
    … …  
    6060c    Argument added to account for subgrid gustiness :
    6161
    62       REAL wmax(ngridmx)
     62      REAL wmax(ngridmx), hfmax(ngridmx)
    6363
    6464c    Traceurs :
    … …  
    7878      REAL z4st,zdplanck(ngridmx)
    7979      REAL zkv(ngridmx,nlayermx+1),zkh(ngridmx,nlayermx+1)
     80      REAL zkq(ngridmx,nlayermx+1)
    8081      REAL zcdv(ngridmx),zcdh(ngridmx)
    8182      REAL zcdv_true(ngridmx),zcdh_true(ngridmx)    ! drag coeff are used by the LES to recompute u* and hfx
    … …  
    114115      REAL kmixmin
    115116
    116 
    117117c    Mass-variation scheme :
    118118c    ~~~~~~~
    … …  
    124124      REAL pdtc(ngrid,nlayermx)
    125125      REAL zhs(ngridmx,nlayermx)
    126       REAL cpice,ccond
    127       SAVE ccond,cpice
    128       DATA cpice /1000./
     126      REAL ccond
     127      SAVE ccond
    129128
    130129c     Theta_m formulation for mass-variation scheme :
    … …  
    320319
    321320        IF (callrichsl) THEN
    322         zcdv(:)=zcdv_true(:)*sqrt(zu2(:)+(0.7*wmax(:))**2)     ! subgrid gustiness added by quadratic interpolation with a coeff beta
    323         zcdh(:)=zcdh_true(:)*sqrt(zu2(:)+(1.2*wmax(:))**2)     ! LES comparisons. This parameter is linked to the thermals settings)
     321       zcdv(:)=zcdv_true(:)*sqrt(zu2(:)+(6.*hfmax(:))**2)     ! subgrid gustiness added by quadratic interpolation with a coeff beta
     322       zcdh(:)=zcdh_true(:)*sqrt(zu2(:)+(10.*hfmax(:))**2)     ! LES comparisons. This parameter is linked to the thermals settings)
     323
     324!       zcdv(:)=zcdv_true(:)*sqrt(zu2(:)+(0.7*wmax(:))**2)     ! subgrid gustiness added by quadratic interpolation with a coeff beta
     325!       zcdh(:)=zcdh_true(:)*sqrt(zu2(:)+(1.2*wmax(:))**2)     ! LES comparisons. This parameter is linked to the thermals settings)
     326
     327        ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)+(6.*hfmax(:))**2)
     328        tst(:)=(ph(:,1)-ptsrf(:))*zcdh(:)/ust(:)
    324329!        ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)+(0.7*wmax(:))**2)
    325330!        tst(:)=(ph(:,1)-ptsrf(:))*zcdh(:)/ust(:)
    … …  
    342347!     &              'heat drag','no units',
    343348!     &                         2,zcdh_true)
    344 !        call WRITEDIAGFI(ngridmx,'u*',
     349!        call WRITEDIAGFI(ngridmx,'ust',
    345350!     &              'friction velocity','m/s',2,ust)
    346 !       call WRITEDIAGFI(ngridmx,'T*',
     351!       call WRITEDIAGFI(ngridmx,'tst',
    347352!     &              'friction temperature','K',2,tst)
    348353!        call WRITEDIAGFI(ngridmx,'rm-1',
    … …  
    355360c    ** schema de diffusion turbulente dans la couche limite
    356361c       ----------------------------------------------------
    357 
    358         CALL vdif_kc(ptimestep,g,pzlev,pzlay
     362!
     363       CALL vdif_kc(ptimestep,g,pzlev,pzlay
    359364     &              ,pu,pv,ph,zcdv_true
    360365     &              ,pq2,zkv,zkh,zq)
     366!
     367!      pt(:,:)=ph(:,:)*ppopsk(:,:)
     368!      CALL yamada4(ngrid,nlay,ptimestep,g,r,pplev,pt
     369!     s   ,pzlev,pzlay,pu,pv,ph,zcdv_true,pq2,zkv,zkh,zkq,ust
     370!     s   ,8)
    361371
    362372      if ((doubleq).and.(ngrid.eq.1)) then
Note: See TracChangeset for help on using the changeset viewer.