Ignore:
Timestamp:
Aug 3, 2011, 4:49:20 PM (13 years ago)
Author:
acolaitis
Message:

--- AC 03/08/2011 ---

M 255 libf/phymars/physiq.F
<> Modified to interface new SL parametrization

M 255 libf/phymars/vdif_cd.F
<> New SL parametrization based on a bulk Richardson Monin-Obukhov theory formulation

Stability functions are taken from D.E. England et al. (95)
Similarity functions coefficients based on Dyer and Hicks (70).
Includes thermal roughness length computation, heat and momentum drag coefficient computation
Can be used to output hydrodynamic-related SL quantities (bulk Richardson, turbulent Prandtl number estimation, Reynolds number)

M 255 libf/phymars/thermcell_dqupdown.F90
<> Minor modification to suit picky compilers

M 255 libf/phymars/vdifc.F
<> Now takes into account sub-grid gustiness, evaluated from thermals activity (it's proxy being the maximum vertical velocity)

M 255 libf/phymars/meso_inc/meso_inc_les.F
<> Minor modification : u* is now taken from the new vdifc and not recomputed from a simple law

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

Legend:

Unmodified
Added
Removed

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

    r226 r256  
    33          hfx(ig) = zflubid(ig)-capcal(ig)*zdtsdif(ig)
    44          !! u star in similarity theory in m/s
    5           ust(ig) = 0.4
    6      .               * sqrt( pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) )
    7      .               / log( 1.E+0 + zzlay(ig,1)/z0_default )
     5!          ust(ig) = 0.4
     6!     .               * sqrt( pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) )
     7!     .               / log( 1.E+0 + zzlay(ig,1)/z0_default )
     8
     9! New SL parametrization, ust is more accurately computed in vdif_cd :
     10          ust(ig) = zcdv(ig)
     11
    812         ENDDO   
    9 
    1013!         write (*,*) 'PHYS HFX cp zdts', hfx(100), zflubid(100),
    1114!     .       capcal(100),

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

    r234 r256  
    321321
    322322      REAL lmax_th_out(ngridmx),zmax_th(ngridmx)
    323       REAL wmax_th(ngridmx)
     323      REAL, SAVE :: wmax_th(ngridmx)
    324324      REAL, SAVE :: hfmax_th(ngridmx)
    325325      REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx)
    … …  
    327327      INTEGER lmax_th(ngridmx)
    328328      REAL dtke_th(ngridmx,nlayermx+1)
    329       REAL dummycol(ngridmx)
    330 
     329      REAL zcdv(ngridmx)
    331330c=======================================================================
    332331
    … …  
    697696     $        zdum1,zdum2,zdh,pdq,zflubid,
    698697     $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
    699      &        zdqdif,zdqsdif)
     698     &        zdqdif,zdqsdif,wmax_th,zcdv)
    700699
    701700#ifdef MESOSCALE
    … …  
    16411640
    16421641         co2col(:)=0.
    1643          dummycol(:)=0.
    16441642         if (tracer) then
    16451643         do l=1,nlayermx
    … …  
    16471645             co2col(ig)=co2col(ig)+
    16481646     &                  zq(ig,l,1)*(pplev(ig,l)-pplev(ig,l+1))/g
    1649          if (nqmx .gt. 1) then
    1650              iq=2 ! to avoid out-of-bounds spotting by picky compilers
    1651                   ! when gcm is compiled with only one tracer
    1652              dummycol(ig)=dummycol(ig)+
    1653      &                  zq(ig,l,iq)*(pplev(ig,l)-pplev(ig,l+1))/g
    1654          endif
    16551647         enddo
    16561648         enddo
    … …  
    16591651         call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass'          &
    16601652     &                                      ,'kg/m-2',0,co2col)
    1661          call WRITEDIAGFI(ngrid,'dummycol','integrated dummy mass'      &
    1662      &                                      ,'kg/m-2',0,dummycol)
    16631653         endif
    16641654         call WRITEDIAGFI(ngrid,'w','vertical velocity'                 &

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

    r185 r256  
    6161! of fmd in the equations, so it has to be positive
    6262!
    63 !      fmd(:,:)=-fm_down(:,:)
     63      fmd(:,:)=-fm_down(:,:)
    6464!
    6565!! ========== Entrainment, Detrainement and Mass =================

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

    r224 r256  
    1       SUBROUTINE vdif_cd( ngrid,nlay,pz0,pg,pz,pu,pv,pts,ph,pcdv,pcdh)
     1      SUBROUTINE vdif_cd(ngrid,nlay,pz0,
     2     &           pg,pz,pu,pv,pts,ph,pcdv,pcdh)
    23      IMPLICIT NONE
    34c=======================================================================
    … …  
    78c
    89c   Author: Frederic Hourdin  15 /10 /93
     10c   Modified by : Arnaud Colaitis 03/08/11
    911c   -------
    1012c
    … …  
    1618c     ngrid            size of the horizontal grid
    1719c     pg               gravity (m s -2)
    18 c     pz(ngrid)        height of the first atmospheric layer
    19 c     pu(ngrid)        u component of the wind in that layer
    20 c     pv(ngrid)        v component of the wind in that layer
    21 c     pts(ngrid)       surfacte temperature
     20c     pz(ngrid,nlay)   height of layers
     21c     pu(ngrid,nlay)   u component of the wind
     22c     pv(ngrid,nlay)   v component of the wind
     23c     pts(ngrid)       surface temperature
    2224c     ph(ngrid)        potential temperature T*(p/ps)^kappa
    2325c
    … …  
    3335c   -------------
    3436
     37#include "comcstfi.h"
     38
    3539c   Arguments:
    3640c   ----------
    3741
    38       INTEGER ngrid,nlay
    39       REAL pz0(ngrid)
    40       REAL pg,pz(ngrid,nlay)
    41       REAL pu(ngrid,nlay),pv(ngrid,nlay)
    42       REAL pts(ngrid,nlay),ph(ngrid,nlay)
    43       REAL pcdv(ngrid),pcdh(ngrid)
     42      INTEGER, INTENT(IN) :: ngrid,nlay
     43      REAL, INTENT(IN) :: pz0(ngrid)
     44      REAL, INTENT(IN) :: pg,pz(ngrid,nlay)
     45      REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)
     46      REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay)
     47      REAL, INTENT(OUT) :: pcdv(ngrid),pcdh(ngrid) ! momentum and heat drag coefficient
    4448
    4549c   Local:
    … …  
    4852      INTEGER ig
    4953
    50       REAL zu2,z1,zri,zcd0,zz
    51 
    52       REAL karman,b,c,d,c2b,c3bc,c3b,umin2
     54      REAL karman
    5355      LOGICAL firstcal
    54       DATA karman,b,c,d,umin2/.4,5.,5.,5.,1.e-12/
     56      DATA karman/.41/
    5557      DATA firstcal/.true./
    56       SAVE b,c,d,karman,c2b,c3bc,c3b,firstcal,umin2
    57 
     58      SAVE karman
     59
     60c    Local(2):
     61c    ---------
     62
     63      REAL rib(ngrid)  ! Bulk Richardson number
     64      REAL fm(ngrid) ! stability function for momentum
     65      REAL fh(ngrid) ! stability function for heat
     66      REAL z1z0,z1z0t ! ratios z1/z0 and z1/z0T
     67
     68c phim = 1+betam*zeta   or   (1-bm*zeta)**am
     69c phih = alphah + betah*zeta    or   alphah(1.-bh*zeta)**ah
     70      REAL betam, betah, alphah, bm, bh, lambda
     71c ah and am are assumed to be -0.25 and -0.5 respectively
     72
     73      REAL cdn(ngrid),chn(ngrid)  ! neutral momentum and heat drag coefficient
     74      REAL pz0t        ! initial thermal roughness length. (local)
     75      REAL ric         ! critical richardson number
     76      REAL prandtl(ngrid)     ! prandtl number for UBL
     77      REAL reynolds(ngrid)    ! reynolds number for UBL
     78      REAL pz0tcomp(ngrid)     ! computed z0t
     79      REAL ite
     80      REAL residual
    5881c-----------------------------------------------------------------------
    5982c   couche de surface:
    6083c   ------------------
    6184
    62 c     DO ig=1,ngrid
    63 c        zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2
    64 c        pcdv(ig)=pz0(ig)*(1.+sqrt(zu2))
    65 c        pcdh(ig)=pcdv(ig)
    66 c     ENDDO
    67 c     RETURN
    68 c
    69 c      IF (firstcal) THEN
    70 c         c2b=2.*b
    71 c         c3bc=3.*b*c
    72 c         c3b=3.*b
    73 c         firstcal=.false.
    74 c      ENDIF
    75 c
    76 c!!!! WARNING, verifier la formule originale de Louis!
    77 c      DO ig=1,ngrid
    78 c         zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2
    79 c         zri=pg*pz(ig)*(ph(ig)-pts(ig))/(ph(ig)*zu2)
    80 c._.
    81 c            zri=0.E+0
    82 c._.
    83 c         z1=1.+pz(ig)/pz0(ig)
    84 c         zcd0=karman/log(z1)
    85 c._.         zcd0=zcd0*zcd0*sqrt(zu2)
    86 c         zcd0=zcd0*zcd0
    87 c         IF(zri.LT.0.) THEN
    88 c._.            z1=b*zri/(1.+c3bc*zcd0*sqrt(-z1*zri))
    89 c            z1=b*zri/(1.+c3bc*zcd0*sqrt(-z1*zri*zu2))
    90 c            pcdv(ig)=zcd0*(1.-2.*z1)
    91 c            pcdh(ig)=zcd0*(1.-3.*z1)
    92 c         ELSE
    93 c            zz=sqrt(1.+d*zri)
    94 c            pcdv(ig)=zcd0/(1.+c2b*zri/zz)
    95 c            pcdh(ig)=zcd0/(1.+c3b*zri*zz)
    96 c         ENDIF
    97 c      ENDDO
    98 
    99 
    100 c On calcule un VRAI cdrag tout bete
     85      reynolds(:)=0.
     86
     87c Original formulation :
     88
     89!      DO ig=1,ngrid
     90!         z1=1.E+0 + pz(ig,1)/pz0(ig)
     91!         zcd0=karman/log(z1)
     92!         zcd0=zcd0*zcd0
     93!         pcdv(ig)=zcd0
     94!         pcdh(ig)=zcd0
     95!      ENDDO
     96
     97!      print*,'old : cd,ch; ',pcdv,pcdh
     98
     99c New formulation (AC) :
     100
     101c phim = 1+betam*zeta   or   (1-bm*zeta)**am
     102c phih = alphah + betah*zeta    or   alphah(1.-bh*zeta)**ah
     103c am=-0.5, ah=-0.25
     104
     105      pz0t = 0.     ! for the sake of simplicity
     106      pz0tcomp(:) = 0.1*pz0(:)
     107      rib(:)=0.
     108      pcdv(:)=0.
     109      pcdh(:)=0.
     110
     111c this formulation assumes alphah=1., implying betah=betam
     112c We use Dyer et al. parameters, as they cover a broad range of Richardson numbers :
     113      bm=16.            !UBL
     114      bh=16.            !UBL
     115      alphah=1.
     116      betam=5.         !SBL
     117      betah=5.         !SBL
     118      lambda=(sqrt(bh/bm))/alphah
     119      ric=betah/(betam**2)
    101120
    102121      DO ig=1,ngrid
    103          z1=1.E+0 + pz(ig,1)/pz0(ig)
    104          zcd0=karman/log(z1)
    105          zcd0=zcd0*zcd0
    106          pcdv(ig)=zcd0
    107          pcdh(ig)=zcd0
     122
     123         ite=0.
     124         residual=abs(pz0tcomp(ig)-pz0t)
     125
     126         do while((residual .gt. 0.01*pz0(ig)) .and.  (ite .lt. 10.))
     127
     128         pz0t=pz0tcomp(ig)
     129
     130         if ((pu(ig,1) .ne. 0.) .or. (pv(ig,1) .ne. 0.)) then
     131
     132c Classical Richardson number formulation
     133
     134c         rib(ig) = (pg/ph(ig,1))*pz(ig,1)*(ph(ig,1)-pts(ig))
     135c     &           /(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1))
     136
     137c Richardson number formulation proposed by D.E. England et al. (1995)
     138
     139          rib(ig) = (pg/ph(ig,1))
     140     &      *pz(ig,1)*pz0(ig)/sqrt(pz(ig,1)*pz0t)
     141     &      *(((log(pz(ig,1)/pz0(ig)))**2)/(log(pz(ig,1)/pz0t)))
     142     &      *(ph(ig,1)-pts(ig))
     143     &  /(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1))
     144
     145         else
     146         print*,'warning, infinite Richardson at surface'
     147         print*,pu(ig,1),pv(ig,1)
     148         rib(ig) = ric          ! traiter ce cas !
     149         endif
     150
     151         z1z0=pz(ig,1)/pz0(ig)
     152         z1z0t=pz(ig,1)/pz0t
     153
     154         cdn(ig)=karman/log(z1z0)
     155         cdn(ig)=cdn(ig)*cdn(ig)
     156         chn(ig)=cdn(ig)*log(z1z0)/log(z1z0t)
     157
     158c Stable case :
     159      if (rib(ig) .gt. 0.) then
     160c From D.E. England et al. (95)
     161      prandtl(ig)=1.
     162         if(rib(ig) .lt. ric) then
     163c Assuming alphah=1. and bh=bm for stable conditions :
     164            fm(ig)=((ric-rib(ig))/ric)**2
     165            fh(ig)=((ric-rib(ig))/ric)**2
     166         else
     167            fm(ig)=0.
     168            fh(ig)=0.
     169         endif
     170c Unstable case :
     171      else
     172c From D.E. England et al. (95)
     173         fm(ig)=sqrt(1.-lambda*bm*rib(ig))
     174         fh(ig)=(1./alphah)*((1.-lambda*bh*rib(ig))**0.5)*
     175     &                     (1.-lambda*bm*rib(ig))**0.25
     176         prandtl(ig)=alphah*((1.-lambda*bm*rib(ig))**0.25)/
     177     &             ((1.-lambda*bh*rib(ig))**0.5)
     178      endif
     179
     180       reynolds(ig)=sqrt(fm(ig))*sqrt(pu(ig,1)**2 + pv(ig,1)**2)*pz0(ig)
     181     &                   /(log(z1z0))
     182       pz0tcomp(ig)=pz0(ig)*exp(-karman*7.3*
     183     &              (reynolds(ig)**0.25)*(prandtl(ig)**0.5))
     184
     185         
     186      residual = abs(pz0t-pz0tcomp(ig))
     187      ite = ite+1
     188!      print*, "iteration nnumber, residual",ite,residual
     189
     190      enddo  ! of while
     191
     192          pz0t=0.
     193
     194c Drag computation :
     195
     196         pcdv(ig)=cdn(ig)*fm(ig)
     197         pcdh(ig)=chn(ig)*fh(ig)
     198       
    108199      ENDDO
    109200
    110        
    111 
    112 
    113 c-----------------------------------------------------------------------
     201!      print*,'new : cd,ch; ',pcdv,pcdh
     202
     203! Some useful diagnostics :
     204
     205!        call WRITEDIAGFI(ngrid,'Ri',
     206!     &              'Richardson nb','no units',
     207!     &                         2,rib)
     208!        call WRITEDIAGFI(ngrid,'Pr',
     209!     &              'Prandtl nb','no units',
     210!     &                         0,prandtl)
     211!        call WRITEDIAGFI(ngrid,'Re',
     212!     &              'Reynolds nb','no units',
     213!     &                         0,reynolds)
     214!        call WRITEDIAGFI(ngrid,'z0tcomp',
     215!     &              'computed z0t','m',
     216!     &                         2,pz0tcomp)
    114217
    115218      RETURN

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

    r224 r256  
    55     $                pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf,
    66     $                pdudif,pdvdif,pdhdif,pdtsrf,pq2,
    7      $                pdqdif,pdqsdif)
     7     $                pdqdif,pdqsdif,wmax,zcdv)
    88      IMPLICIT NONE
    99
    … …  
    5959      REAL pz0(ngridmx) ! surface roughness length (m)
    6060
     61c    Argument added to account for subgrid gustiness :
     62
     63      REAL wmax(ngridmx)
     64
    6165c    Traceurs :
    6266      integer nq
    … …  
    228232c       ---------------------
    229233
    230       CALL vdif_cd( ngrid,nlay,pz0,g,pzlay,pu,pv,ptsrf,ph
     234      CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pu,pv,ptsrf,ph
    231235     &             ,zcdv_true,zcdh_true)
    232       DO ig=1,ngrid
     236
     237      DO ig=1,ngrid
     238
    233239        zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)
     240     &             +(1.2*wmax(ig))**2
     241
    234242        zcdv(ig)=zcdv_true(ig)*sqrt(zu2)
    235243        zcdh(ig)=zcdh_true(ig)*sqrt(zu2)
    236       ENDDO
     244
     245      ENDDO
     246
     247! Some usefull diagnostics for the new surface layer parametrization :
     248
     249!        call WRITEDIAGFI(ngridmx,'Cd',
     250!     &              'momentum drag','no units',
     251!     &                         2,zcdv_true)
     252!        call WRITEDIAGFI(ngridmx,'Ch',
     253!     &              'heat drag','no units',
     254!     &                         2,zcdh_true)
     255!        call WRITEDIAGFI(ngridmx,'u*',
     256!     &              'friction velocity','m/s',
     257!     &                         2,zcdv)
     258!        call WRITEDIAGFI(ngridmx,'T*',
     259!     &              'friction temperature','m/s',
     260!     &                         2,zcdh)
     261
    237262
    238263c    ** schema de diffusion turbulente dans la couche limite
Note: See TracChangeset for help on using the changeset viewer.