Changeset 4663 for LMDZ6/trunk


Ignore:
Timestamp:
Sep 1, 2023, 4:01:08 PM (14 months ago)
Author:
evignon
Message:

ajout d'une longueur de melange sensible au cisaillement de vent dans atke

Location:
LMDZ6/trunk/libf/phylmd
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90

    r4653 r4663  
    2828USE atke_turbulence_ini_mod, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd, atke_ok_virtual
    2929USE atke_turbulence_ini_mod, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, ctkes,rg, rd, rv, atke_ok_vdiff
    30 USE atke_turbulence_ini_mod, ONLY : viscom, viscoh, clmix, iflag_atke_lmix, lmin, smmin
     30USE atke_turbulence_ini_mod, ONLY : viscom, viscoh, clmix, clmixshear, iflag_atke_lmix, lmin, smmin
    3131
    3232implicit none
     
    6868REAL, DIMENSION(ngrid,nlay+1) :: Sm          ! Stability function for momentum (at interface)
    6969REAL, DIMENSION(ngrid,nlay+1) :: Sh          ! Stability function for heat (at interface)
    70 LOGICAL, DIMENSION(ngrid,nlay+1) :: switch_num ! switch of numerical integration in very stable cases
    7170
    7271INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid)
     
    174173!==============================================================
    175174
    176 switch_num(:,:)=.false.
    177175
    178176IF (iflag_atke_lmix .EQ. 1 ) THEN
     
    183181          IF (N2(igrid,ilay) .GT. 0.) THEN
    184182             lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay))
    185              IF (lstrat .LT. l_exchange(igrid,ilay)) THEN
    186                 l_exchange(igrid,ilay)=max(lstrat,lmin)
    187                 switch_num(igrid,ilay)=.true.
    188              ENDIF
     183             lstrat=max(lstrat,lmin)
     184             !Monin-Obukhov consistent interpolation, Van de Wiel et al. 2010
     185             l_exchange(igrid,ilay)=(1./(2.*l_exchange(igrid,ilay))+sqrt(1./(4.*l_exchange(igrid,ilay) &
     186                                   *l_exchange(igrid,ilay))+1./(2.*lstrat*lstrat)))**(-1.0)
    189187          ENDIF
    190188      ENDDO
     
    192190
    193191ELSE IF (iflag_atke_lmix .EQ. 2 ) THEN
    194 ! add effect of wind shear on lstrat following grisogono and belusic 2008, qjrms, eq 2
     192! add effect of wind shear on lstrat following grisogono and belusic 2008, qjrms
    195193DO ilay=2,nlay
    196194      DO igrid=1,ngrid
    197195          l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
    198           IF (N2(igrid,ilay) .GT. 0.) THEN
    199              lstrat=clmix*sqrt(tke(igrid,ilay))/(2.*max(sqrt(shear2(igrid,ilay)),1.E-10)*(1.+sqrt(Ri(igrid,ilay))/2.))
    200              IF (lstrat .LT. l_exchange(igrid,ilay)) THEN
    201                 l_exchange(igrid,ilay)=max(lstrat,lmin)
    202              ENDIF
     196          IF (N2(igrid,ilay) .GT. 0. .AND. shear2(igrid,ilay) .GT. 0.) THEN
     197             lstrat=min(clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay)), &
     198                    clmixshear*sqrt(tke(igrid,ilay))/sqrt(shear2(igrid,ilay)))
     199             lstrat=max(lstrat,lmin)
     200             !Monin-Obukhov consistent interpolation, Van de Wiel et al. 2010   
     201             l_exchange(igrid,ilay)=(1./(2.*l_exchange(igrid,ilay))+sqrt(1./(4.*l_exchange(igrid,ilay) &
     202                                   *l_exchange(igrid,ilay))+1./(2.*lstrat*lstrat)))**(-1.0)
     203
    203204          ENDIF
    204205      ENDDO
     
    299300
    300301
    301 ELSE IF (iflag_atke == 5) THEN
    302 ! semi implicit scheme from Arpege (V. Masson methodology with
    303 ! Taylor expansion of the dissipation term)
    304 ! and implicit resolution when switch num (when we use the mixing length as a function of tke)
    305     DO ilay=2,nlay
    306         DO igrid=1,ngrid
    307            qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
    308            if (switch_num(igrid,ilay) .and. N2(igrid,ilay)>0.) then
    309              invtau=clmix*Sm(igrid,ilay)/sqrt(N2(igrid,ilay))*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))  &
    310                    -sqrt(N2(igrid,ilay))/(cepsilon*clmix)
    311              !qq=qq/(1.-dtime*invtau)
    312              !qq=qq*exp(dtime*invtau)
    313              tke(igrid,ilay)=tke(igrid,ilay)/(1.-dtime*invtau)
    314              tke(igrid,ilay)=max(0.,tke(igrid,ilay))
    315            else
    316              qq=(l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) &
    317              +qq*(1.+dtime*qq/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))) &
    318              /(1.+2.*qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))
    319              qq=max(0.,qq)
    320              tke(igrid,ilay)=0.5*(qq**2)
    321            endif
    322         ENDDO
    323     ENDDO
    324 
    325 
    326302ELSE
    327303   call abort_physic("atke_compute_km_kh", &
  • LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90

    r4653 r4663  
    99  real :: kappa = 0.4 ! Von Karman constant
    1010  !$OMP THREADPRIVATE(kappa)
    11   real :: l0, ric, ri0, cinf, cepsilon, pr_slope, pr_asym, pr_neut, clmix, smmin, ctkes,cke
    12   !$OMP THREADPRIVATE(l0,ric,cinf,cepsilon,pr_slope,pr_asym,pr_neut,clmix,smmin,ctkes,cke)
     11  real :: l0,ric,ri0,cinf,cepsilon,pr_slope,pr_asym,pr_neut,clmix,clmixshear,smmin,ctkes,cke
     12  !$OMP THREADPRIVATE(l0,ric,cinf,cepsilon,pr_slope,pr_asym,pr_neut,clmix,clmixshear,smmin,ctkes,cke)
    1313  integer :: lunout,prt_level
    1414  !$OMP THREADPRIVATE(lunout,prt_level)
     
    9494
    9595 ! coefficient for surface TKE
    96  ! following Lenderink & Holtslag 2004, ctkes=sqrt(cepsilon)
     96 ! following Lenderink & Holtslag 2004, ctkes=(cepsilon)**(2/3)
    9797 ! (provided by limit condition in neutral conditions)
    98   ctkes=sqrt(cepsilon)
    99 
     98  ctkes=(cepsilon)**(2./3.)
    10099
    101100  ! slope of Pr=f(Ri) for stable conditions
     
    118117  clmix=0.5
    119118  CALL getin_p('atke_clmix',clmix)
    120  
     119 
     120  ! coefficient for mixing length depending on local wind shear
     121  clmixshear=0.5
     122  CALL getin_p('atke_clmixshear',clmixshear)
     123
     124
    121125  ! minimum anisotropy coefficient (defined here as minsqrt(Ez/Ek)) at large Ri.
    122126  ! From Zilitinkevich et al. 2013, it equals sqrt(0.03)~0.17
Note: See TracChangeset for help on using the changeset viewer.