Changeset 268 for trunk/LMDZ.MARS/libf


Ignore:
Timestamp:
Aug 16, 2011, 11:06:45 AM (13 years ago)
Author:
acolaitis
Message:

--- AC 03/08/2011 ---
M 267 libf/phymars/physiq.F
<> Minor modification to pass Ch from vdifc to meso_inc_les

M 267 libf/phymars/surflayer_interpol.F
<> Major modification to the formulation of integrals

Now stable for most cases. Some cases with highly negative Monin Obukhov length
remain to be explored.

M 267 libf/phymars/vdif_cd.F
<> Added gustiness to the Richardson computation. Gustiness factor is for now of beta=1., after

several comparisons with LES aerodynamic conductances. May be subject to a minor change (+/- 0.1)
in the near future. (almost transparent for the user)

M 267 libf/phymars/vdifc.F
<> Minor modifications relative to variables.

M 267 libf/phymars/calltherm_mars.F90
<> Added a comment on a "sensitive" parameter that should not be changed without knowing the consequence !

M 267 libf/phymars/meso_inc/meso_inc_les.F
<> Changed the definition for HFX computation in the LES (to be discussed with Aymeric). New definition yields

very similar results to old one and follows a strict definition of what HFX should be.

Location:
trunk/LMDZ.MARS/libf/phymars
Files:
6 edited

Legend:

Unmodified
Added
Removed

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

    r190 r268  
    8787!  ===================
    8888 
    89          r_aspect_thermals=2.
     89         r_aspect_thermals=2.   ! ultimately conrols the amount of mass going through the thermals
     90                                ! decreasing it increases the thermals effect. Tests at gcm resolution
     91                                ! shows that too low values destabilize the model
     92                                ! when changing this value, one should check that the surface layer model
     93                                ! outputs the correct Cd*u and Ch*u through changing the gustiness coefficient beta
    9094         nsplit_thermals=40
    9195         call getin("nsplit_thermals",nsplit_thermals)

  • trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_les.F

    r265 r268  
    11         DO ig=1,ngrid
    22          !! sensible heat flux in W/m2
    3           hfx(ig) = zflubid(ig)-capcal(ig)*zdtsdif(ig)
     3!          hfx(ig) = zflubid(ig)-capcal(ig)*zdtsdif(ig)
     4
     5! New SL parametrization, correct formulation for hfx :
     6
     7          hfx(ig) = (pplay(ig,1)/(r*pt(ig,1)))*cpp
     8     &    *sqrt((pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1)))
     9     &    *zcdh(ig)*(tsurf(ig)-zh(ig,1))
     10
    411          !! u star in similarity theory in m/s
    512!          ust(ig) = 0.4

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

    r267 r268  
    327327      INTEGER lmax_th(ngridmx)
    328328      REAL dtke_th(ngridmx,nlayermx+1)
    329       REAL zcdv(ngridmx)
     329      REAL zcdv(ngridmx), zcdh(ngridmx)
    330330      REAL Teta_out(ngridmx),u_out(ngridmx)  ! Interpolated teta and u at z_out
    331331      REAL z_out                          ! height of interpolation between z0 and z1
    … …  
    346346         call zerophys(ngrid*nlayer,dtrad)
    347347         call zerophys(ngrid,fluxrad)
     348
     349         wmax_th(:)=0.
    348350
    349351c        read startfi
    … …  
    700702     $        zdum1,zdum2,zdh,pdq,zflubid,
    701703     $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
    702      &        zdqdif,zdqsdif,wmax_th,zcdv)
     704     &        zdqdif,zdqsdif,wmax_th,zcdv,zcdh)
    703705
    704706#ifdef MESOSCALE
    … …  
    16791681! THERMALS STUFF 1D
    16801682
     1683         z_out=0.
     1684         if (calltherm .and. (z_out .gt. 0.)) then
     1685
     1686         call surflayer_interpol(ngrid,nlayer,z0,g,zzlay,zu,zv,wmax_th
     1687     &              ,tsurf,zt(:,:)*(zplay(:,:)/zplev(:,:))**rcp
     1688     &              ,z_out,Teta_out,u_out,ustar,tstar)
     1689
     1690         zu2(:)=sqrt(zu(:,1)*zu(:,1)+zv(:,1)*zv(:,1))
     1691         call WRITEDIAGFI(ngridmx,'sqrt(zu2)',
     1692     &              'horizontal velocity norm','m/s',
     1693     &                         0,zu2)
     1694
     1695         call WRITEDIAGFI(ngridmx,'Teta_out',
     1696     &              'potential temperature at z_out','K',
     1697     &                         0,Teta_out)
     1698         call WRITEDIAGFI(ngridmx,'u_out',
     1699     &              'horizontal velocity norm at z_out','m/s',
     1700     &                         0,u_out)
     1701         call WRITEDIAGFI(ngridmx,'u*',
     1702     &              'friction velocity','m/s',
     1703     &                         0,ustar)
     1704         call WRITEDIAGFI(ngridmx,'teta*',
     1705     &              'friction potential temperature','K',
     1706     &                         0,tstar)
     1707         else
     1708           if((.not. calltherm).and.(z_out .gt. 0.)) then
     1709            print*, 'WARNING : no interpolation in surface-layer :'
     1710            print*, 'Outputing surface-layer quantities without thermals
     1711     & does not make sense'
     1712           endif
     1713         endif
     1714
     1715! ---
     1716
     1717
     1718
    16811719         if(calltherm) then
    16821720
    … …  
    16931731     &                         0,wmax_th)
    16941732
    1695 
    16961733         co2col(:)=0.
    16971734         if (tracer) then
    … …  
    17121749         call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
    17131750     &                  tsurf)
     1751!         call WRITEDIAGFI(ngrid,"u","u wind","m/s",1,zu)
     1752!         call WRITEDIAGFI(ngrid,"v","v wind","m/s",1,zv)
    17141753
    17151754         call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)

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

    r267 r268  
    55!
    66!   Subject: interpolation of temperature and velocity norm in the surface layer
    7 !   by recomputation of surface layer quantities (Richardson, Prandtl, u*, teta*)
     7!   by recomputation of surface layer quantities (Richardson, Prandtl, Reynolds, u*, teta*)
    88!   ------- 
    99!
    … …  
    4242      REAL, INTENT(IN) :: wmax(ngrid)
    4343      REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay)
    44       REAL, INTENT(INOUT) :: z_out                      ! output height (in m above surface)
    45       REAL, INTENT(OUT) :: Teta_out(ngrid),u_out(ngrid) ! interpolated fields at z_out : potential temperature and norm(uv)
     44      REAL, INTENT(INOUT) :: z_out                  ! output height (in m above surface)
     45      REAL, INTENT(OUT) :: Teta_out(ngrid),u_out(ngrid)! interpolated fields at z_out : potential temperature and norm(uv)
    4646      REAL, INTENT(OUT) :: ustar(ngrid), tstar(ngrid) ! u* and teta*
    4747
    … …  
    7979! For output :
    8080
    81       REAL zu2                 ! Large-scale wind at first layer
     81      REAL zu2(ngrid)                  ! Large-scale wind at first layer
    8282      REAL L_mo(ngrid)                ! Monin-Obukhov length
    8383!-----------------------------------------------------------------------
    … …  
    8888      ustar(:)=0.
    8989      reynolds(:)=0.
    90 
    91 ! Original formulation :
    92 
    93 !      DO ig=1,ngrid
    94 !         z1=1.E+0 + pz(ig,1)/pz0(ig)
    95 !         zcd0=karman/log(z1)
    96 !         zcd0=zcd0*zcd0
    97 !         pcdv(ig)=zcd0
    98 !         pcdh(ig)=zcd0
    99 !      ENDDO
    100 
    101 !      print*,'old : cd,ch; ',pcdv,pcdh
    10290
    10391! New formulation (AC) :
    … …  
    141129! Richardson number formulation proposed by D.E. England et al. (1995)
    142130
     131!        IF((pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) .lt. 1.)
     132!     &            .and. (wmax(ig) .gt. 0.)) THEN
     133          zu2(ig)=
     134!     &  (MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1),wmax(ig)**2))
     135     &  ( pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) + wmax(ig)**2)
     136!     &  pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1)
     137!        ELSE
     138!        zu2(ig)=pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1)
     139!        ENDIF
     140
    143141          rib(ig) = (pg/ph(ig,1))
    144      &      *pz(ig,1)*pz0(ig)/sqrt(pz(ig,1)*pz0t)
    145 !     &      *sqrt(pz(ig,1)*pz0(ig))
     142!     &      *pz(ig,1)*pz0(ig)/sqrt(pz(ig,1)*pz0t)
     143     &      *sqrt(pz(ig,1)*pz0(ig))
    146144     &      *(((log(pz(ig,1)/pz0(ig)))**2)/(log(pz(ig,1)/pz0t)))
    147      &      *(ph(ig,1)-pts(ig))
    148      &  /(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1))
     145     &      *(ph(ig,1)-pts(ig))/zu2(ig)
     146
     147!     &  /(MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1),wmax(ig)**2))
     148!     &  /( pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) + wmax(ig)**2)
    149149
    150150         else
    … …  
    208208      ENDDO
    209209
    210 ! Large-scale wind at first layer (with gustiness) and
     210! Large-scale wind at first layer (without gustiness) and
    211211! u* theta* computation
    212212      DO ig=1,ngrid
    … …  
    216216           tstar(ig)=0.
    217217         else
    218            zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)+(0.*wmax(ig))**2
    219            ustar(ig)=karman*sqrt(fm(ig)*zu2)/(log(pz(ig,1)/pz0(ig)))
    220            tstar(ig)=karman*fh(ig)*(ph(ig,1)-pts(ig))
    221      &                   /(log(pz(ig,1)/pz0tcomp(ig))*sqrt(fm(ig)))
     218
     219!           ustar(ig)=karman*sqrt(fm(ig)*zu2(ig))/(log(pz(ig,1)/pz0(ig)))
     220!           tstar(ig)=karman*fh(ig)*(ph(ig,1)-pts(ig))
     221!     &                   /(log(pz(ig,1)/pz0tcomp(ig))*sqrt(fm(ig)))
     222
     223!simpler definition of u* and teta*, equivalent to the formula above :
     224
     225            ustar(ig)=sqrt(pcdv(ig))*sqrt(zu2(ig))                   
     226            tstar(ig)=-pcdh(ig)*(pts(ig)-ph(ig,1))/sqrt(pcdv(ig))
     227
     228           if (tstar(ig) .lt. -50) then
     229             print*, fh(ig),rib(ig),(ph(ig,1)-pts(ig))
     230     &               ,log(pz(ig,1)/pz0tcomp(ig)),sqrt(fm(ig))
     231           endif
    222232         endif
    223233      ENDDO
    … …  
    233243      ENDDO
    234244
     245      DO ig=1,ngrid
    235246      IF(z_out .ge. pz(ig,1)) THEN
    236247           z_out=1.
    … …  
    252263     &                  and computation is resumed'
    253264      ENDIF
    254 
     265      ENDDO
    255266
    256267      print*, 'interpolation of u and teta at z_out=',z_out
    … …  
    265276     &                            *(z_out-pz0tcomp(ig))
    266277        ELSEIF (L_mo(ig) .lt. 0.) THEN
    267            u_out(ig)=(ustar(ig)/karman)*((
     278
     279          IF(L_mo(ig) .gt. -1000.) THEN
     280           
     281           u_out(ig)=(ustar(ig)/karman)*(
    268282     &  2.*atan((1.-16.*z_out/L_mo(ig))**0.25)
    269      &  +log((-1.+(1.-16.*z_out/L_mo(ig))**0.25)
    270      &  /(-1.+(1.-16.*z_out/L_mo(ig))**0.25)))-(
    271      &  2.*atan((1.-16.*pz0(ig)/L_mo(ig))**0.25)
    272      &  +log((-1.+(1.-16.*pz0(ig)/L_mo(ig))**0.25)
    273      &  /(-1.+(1.-16.*pz0(ig)/L_mo(ig))**0.25)))
     283     & -2.*atan((1.-16.*pz0(ig)/L_mo(ig))**0.25)
     284     & -2.*log(1.+(1.-16.*z_out/L_mo(ig))**0.25)
     285     & +2.*log(1.+(1.-16.*pz0(ig)/L_mo(ig))**0.25)
     286     & -   log(1.+sqrt(1.-16.*z_out/L_mo(ig)))
     287     & +   log(1.+sqrt(1.-16.*pz0(ig)/L_mo(ig)))
     288     & +   log(z_out/pz0(ig))
    274289     &                                  )
    275            Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*((
    276      &  log((-1.+sqrt(1.-16.*z_out/L_mo(ig)))
    277      &  /(1.+sqrt(1.-16.*z_out/L_mo(ig)))))-(
    278      &  log((-1.+sqrt(1.-16.*pz0tcomp(ig)/L_mo(ig)))
    279      &  /(1.+sqrt(1.-16.*pz0tcomp(ig)/L_mo(ig)))))
     290
     291           Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*(
     292     &  2.*log(1.+sqrt(1.-16.*pz0tcomp(ig)/L_mo(ig)))
     293     & -2.*log(1.+sqrt(1.-16.*z_out/L_mo(ig)))
     294     & +   log(z_out/pz0tcomp(ig))
    280295     &                                                        )
     296
     297          ELSE
     298
     299! We have to treat the case where L is very large and negative,
     300! corresponding to a very nearly stable atmosphere (but not quite !)
     301! i.e. teta* <0 and almost zero. We choose the cutoff
     302! corresponding to L_mo < -1000 and use a 3rd order taylor expansion. The difference
     303! between the two expression for z/L = -1/1000 is -1.54324*10^-9
     304! (we do that to avoid using r*4 precision, otherwise, we get -inf values)         
     305
     306           u_out(ig)=(ustar(ig)/karman)*(
     307     &     (4./L_mo(ig))*(z_out-pz0(ig))
     308     &   + (20./(L_mo(ig))**2)*(z_out**2-pz0(ig)**2)
     309     &   + (160./(L_mo(ig))**3)*(z_out**3-pz0(ig)**3)
     310     &   + log(z_out/pz0(ig))
     311     &                                   )
     312
     313           Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*(
     314     &     (8./L_mo(ig))*(z_out-pz0tcomp(ig))
     315     &   + (48./(L_mo(ig))**2)*(z_out**2-pz0tcomp(ig)**2)
     316     &   + (1280./(3.*(L_mo(ig))**3))*(z_out**3-pz0tcomp(ig)**3)
     317     &   + log(z_out/pz0tcomp(ig))
     318     &                                                           )
     319
     320          ENDIF
    281321        ELSE
    282322           u_out(ig)=0.
    283            Teta_out(ig)=0.
     323           Teta_out(ig)=pts(ig)
    284324        ENDIF
    285325      ENDDO
     326
     327! Usefull diagnostics for the interpolation routine :
    286328
    287329         call WRITEDIAGFI(ngrid,'L',

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

    r267 r268  
    11      SUBROUTINE vdif_cd(ngrid,nlay,pz0,
    2      &           pg,pz,pu,pv,pts,ph,pcdv,pcdh)
     2     &           pg,pz,pu,pv,wmax,pts,ph,pcdv,pcdh)
    33      IMPLICIT NONE
    44c=======================================================================
    … …  
    4545      REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
    4646      REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay)
     47      REAL, INTENT(IN) :: wmax(ngrid)
    4748      REAL, INTENT(OUT) :: pcdv(ngrid),pcdh(ngrid) ! momentum and heat drag coefficient
    4849
    … …  
    6263
    6364      REAL rib(ngrid)  ! Bulk Richardson number
     65      REAL rig(ngrid)  ! Gradient Richardson number
    6466      REAL fm(ngrid) ! stability function for momentum
    6567      REAL fh(ngrid) ! stability function for heat
    … …  
    7981      REAL ite
    8082      REAL residual
     83      REAL zu2
    8184c-----------------------------------------------------------------------
    8285c   couche de surface:
    … …  
    137140c Richardson number formulation proposed by D.E. England et al. (1995)
    138141
     142!         zu2=MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1),0.25*wmax(ig)**2)
     143!         zu2=pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1)
     144!         zu2=MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1),wmax(ig)**2)
     145         zu2= pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) + (1.*wmax(ig))**2
     146
     147               ! we add the wmax to simulate
     148               ! bulk Ri changes due to subgrid wind feeding the thermals
     149
     150!          rig(ig) = (pg/ph(ig,1))*((pz(ig,1)-pz0(ig))**2
     151!     &         /(pz(ig,1)-pz0t))*(ph(ig,1)-pts(ig))
     152!     &         /zu2
     153
    139154          rib(ig) = (pg/ph(ig,1))
    140      &      *pz(ig,1)*pz0(ig)/sqrt(pz(ig,1)*pz0t)
     155!     &      *pz(ig,1)*pz0(ig)/sqrt(pz(ig,1)*pz0t)
     156     &      *sqrt(pz(ig,1)*pz0(ig))
    141157     &      *(((log(pz(ig,1)/pz0(ig)))**2)/(log(pz(ig,1)/pz0t)))
    142158     &      *(ph(ig,1)-pts(ig))
    143      &  /(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1))
     159     &  /zu2
    144160
    145161         else
    … …  
    204220! Some useful diagnostics :
    205221
    206 !        call WRITEDIAGFI(ngrid,'Ri',
    207 !     &              'Richardson nb','no units',
     222!        call WRITEDIAGFI(ngrid,'RiB',
     223!     &              'Bulk Richardson nb','no units',
    208224!     &                         2,rib)
     225!                call WRITEDIAGFI(ngrid,'RiG',
     226!     &              'Grad Richardson nb','no units',
     227!     &                         2,rig)
    209228!        call WRITEDIAGFI(ngrid,'Pr',
    210229!     &              'Prandtl nb','no units',

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

    r267 r268  
    55     $                pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
    66     $                pdudif,pdvdif,pdhdif,pdtsrf,pq2,
    7      $                pdqdif,pdqsdif,wmax,zcdv_true)
     7     $                pdqdif,pdqsdif,wmax,zcdv_true,zcdh_true)
    88      IMPLICIT NONE
    99
    … …  
    7878      REAL zkv(ngridmx,nlayermx+1),zkh(ngridmx,nlayermx+1)
    7979      REAL zcdv(ngridmx),zcdh(ngridmx)
    80       REAL zcdv_true(ngridmx),zcdh_true(ngridmx)
     80      REAL zcdv_true(ngridmx),zcdh_true(ngridmx)    ! drag coeff are used by the LES to recompute u* and hfx
    8181      REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
    8282      REAL zh(ngridmx,nlayermx)
    … …  
    232232c       ---------------------
    233233
    234       CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pu,pv,ptsrf,ph
     234      CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pu,pv,wmax,ptsrf,ph
    235235     &             ,zcdv_true,zcdh_true)
    236236
    … …  
    238238
    239239        zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)
    240      &             +(1.2*wmax(ig))**2        !subgrid gustiness is used to enhance surface flux only, and not u*,t* computations
     240     &             +(1.*wmax(ig))**2        ! subgrid gustiness added by quadratic interpolation with a coeff beta, here beta=1. (tuned from
     241                                            ! LES comparisons. This parameter is linked to the thermals settings)
    241242       
    242         zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
    243         zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
    244 
     243        zcdv(ig)=zcdv_true(ig)*sqrt(zu2)     ! 1 / bulk aerodynamic momentum conductance
     244        zcdh(ig)=zcdh_true(ig)*sqrt(zu2)     ! 1 / bulk aerodynamic heat conductance
     245                                             ! these are the quantities to be looked at when comparing surface layers of different models
    245246      ENDDO
    246247
    247248! Some usefull diagnostics for the new surface layer parametrization :
    248249
    249 !        call WRITEDIAGFI(ngridmx,'Cd',
     250!         call WRITEDIAGFI(ngridmx,'pcdv',
    250251!     &              'momentum drag','no units',
    251252!     &                         2,zcdv_true)
    252 !        call WRITEDIAGFI(ngridmx,'Ch',
     253!        call WRITEDIAGFI(ngridmx,'pcdh',
    253254!     &              'heat drag','no units',
    254255!     &                         2,zcdh_true)
    255 !        call WRITEDIAGFI(ngridmx,'u*2/u',
    256 !     &              'friction velocity','m/s',
     256!        call WRITEDIAGFI(ngridmx,'u*',
     257!    &              'friction velocity','m/s',
     258!     &                         2,sqrt(zcdv_true(:))*sqrt(zu2))
     259!        call WRITEDIAGFI(ngridmx,'T*',
     260!     &              'friction temperature','K',
     261!     &          2,-zcdh_true(:)*(ptsrf(:)-ph(:,1))/sqrt(zcdv_true(:)))
     262!        call WRITEDIAGFI(ngridmx,'rm-1',
     263!     &              'aerodyn momentum conductance','m/s',
    257264!     &                         2,zcdv)
    258 !        call WRITEDIAGFI(ngridmx,'u*T*/(T0-T)',
    259 !     &              'friction temperature','m/s',
     265!        call WRITEDIAGFI(ngridmx,'rh-1',
     266!     &              'aerodyn heat conductance','m/s',
    260267!     &                         2,zcdh)
    261 
    262268
    263269c    ** schema de diffusion turbulente dans la couche limite
Note: See TracChangeset for help on using the changeset viewer.