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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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
Note: See TracChangeset for help on using the changeset viewer.