Changeset 4653 for LMDZ6


Ignore:
Timestamp:
Aug 30, 2023, 9:15:43 AM (15 months ago)
Author:
evignon
Message:

prise en compte de l'humidite pour le calcul du flux de flottabilite dans atke
+ petits ajustements

Location:
LMDZ6/trunk/libf
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90

    r4644 r4653  
    66
    77subroutine atke_compute_km_kh(ngrid,nlay,dtime, &
    8                         wind_u,wind_v,temp,play,pinterf,cdrag_uv, &
     8                        wind_u,wind_v,temp,qvap,play,pinterf,cdrag_uv, &
    99                        tke,Km_out,Kh_out)
    1010
    … …  
    2626
    2727
    28 USE atke_turbulence_ini_mod, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd, ok_vdiff_tke
    29 USE atke_turbulence_ini_mod, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, ctkes,rg, rd
     28USE atke_turbulence_ini_mod, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd, atke_ok_virtual
     29USE atke_turbulence_ini_mod, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, ctkes,rg, rd, rv, atke_ok_vdiff
    3030USE atke_turbulence_ini_mod, ONLY : viscom, viscoh, clmix, iflag_atke_lmix, lmin, smmin
    3131
    … …  
    4343REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: wind_v   ! meridional velocity (m/s)
    4444REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: temp   ! temperature (K)
     45REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: qvap   ! specific humidity (kg/kg)
    4546REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: play   ! pressure (Pa)
    4647REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: pinterf   ! pressure at interfaces(Pa)
    … …  
    8182REAL    :: netsource  ! net source term of tke
    8283REAL    :: ustar      ! friction velocity estimation
     84REAL    :: invtau     
     85REAL    :: rvap
    8386
    8487! Initializations:
    … …  
    9093END DO
    9194
    92 ! Calculation of potential temperature: (if vapor -> todo virtual potential temperature)
     95! Calculation of potential temperature: (if vapor -> virtual potential temperature)
    9396!=====================================
    9497
    9598preff=100000.
    96 ! The result should not depend on the choice of preff
     99! results should not depend on the choice of preff
    97100DO ilay=1,nlay
    98101     DO igrid = 1, ngrid
    … …  
    101104END DO
    102105
     106! account for water vapor mass for buoyancy calculation
     107IF (atke_ok_virtual) THEN
     108  DO ilay=1,nlay
     109     DO igrid = 1, ngrid
     110        rvap=max(0.,qvap(igrid,ilay)/(1.-qvap(igrid,ilay)))
     111        theta(igrid,ilay)=theta(igrid,ilay)*(1.+rvap/(RD/RV))/(1.+rvap)
     112     END DO
     113  END DO
     114ENDIF
    103115
    104116
    … …  
    287299
    288300
     301ELSE 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
    289326ELSE
    290327   call abort_physic("atke_compute_km_kh", &
    … …  
    312349! vertical diffusion of TKE
    313350!==========================
    314 IF (ok_vdiff_tke) THEN
     351IF (atke_ok_vdiff) THEN
    315352   CALL atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke)
    316353ENDIF
    … …  
    338375
    339376! routine that computes the vertical diffusion of TKE by the turbulence
    340 ! using an implicit resolution (See note by Dufresne and Ghattas
     377! using an implicit resolution (See note by Dufresne and Ghattas (2009))
    341378! E Vignon, July 2023
    342379

  • LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90

    r4644 r4653  
    1313  integer :: lunout,prt_level
    1414  !$OMP THREADPRIVATE(lunout,prt_level)
    15   real :: rg, rd, rpi, rcpd
    16   !$OMP THREADPRIVATE(rg, rd, rpi, rcpd)
     15  real :: rg, rd, rpi, rcpd, rv
     16  !$OMP THREADPRIVATE(rg, rd, rpi, rcpd, rv)
    1717
    1818  real :: viscom, viscoh
    … …  
    2222  !$OMP THREADPRIVATE(lmin)
    2323
    24   logical :: ok_vdiff_tke
    25   !$OMP THREADPRIVATE(ok_vdiff_tke)
     24  logical :: atke_ok_vdiff, atke_ok_virtual
     25  !$OMP THREADPRIVATE(atke_ok_vdiff,atke_ok_virtual)
    2626
    2727CONTAINS
    2828
    29 SUBROUTINE atke_ini(prt_level_in, lunout_in, rg_in, rd_in, rpi_in, rcpd_in)
     29SUBROUTINE atke_ini(prt_level_in, lunout_in, rg_in, rd_in, rpi_in, rcpd_in, rv_in)
    3030
    3131   USE ioipsl_getin_p_mod, ONLY : getin_p
    3232
    3333  integer, intent(in) :: lunout_in,prt_level_in
    34   real, intent(in) :: rg_in, rd_in, rpi_in, rcpd_in
     34  real, intent(in) :: rg_in, rd_in, rpi_in, rcpd_in, rv_in
    3535
    3636
    … …  
    4141  rpi=rpi_in
    4242  rcpd=rcpd_in
     43  rv=rv_in
    4344
    4445  viscom=1.46E-5
    … …  
    5960
    6061  ! activate vertical diffusion of TKE or not
    61   ok_vdiff_tke=.false.
    62   CALL getin_p('atke_ok_vdiff_tke',ok_vdiff_tke)
     62  atke_ok_vdiff=.false.
     63  CALL getin_p('atke_ok_vdiff',atke_ok_vdiff)
     64
     65
     66  ! account for vapor for flottability
     67  atke_ok_virtual=.true.
     68  CALL getin_p('atke_ok_virtual',atke_ok_virtual)
     69
    6370
    6471  ! flag that controls the numerical treatment of diffusion coeffiient calculation
    … …  
    8289
    8390  ! constant for tke dissipation calculation
    84   cepsilon=16.6/2./sqrt(2.) ! default value as in yamada4
     91  cepsilon=5.87 ! default value as in yamada4
    8592  CALL getin_p('atke_cepsilon',cepsilon)
     93
     94
     95 ! coefficient for surface TKE
     96 ! following Lenderink & Holtslag 2004, ctkes=sqrt(cepsilon)
     97 ! (provided by limit condition in neutral conditions)
     98  ctkes=sqrt(cepsilon)
     99
    86100
    87101  ! slope of Pr=f(Ri) for stable conditions
    … …  
    111125  CALL getin_p('atke_smmin',smmin)
    112126
    113  ! coefficient for surface TKE (default value from Arpege, see E. Bazile note)
    114   ctkes=3.75
    115   CALL getin_p('atke_ctkes',ctkes)
    116127
    117128  ! ratio between the eddy diffusivity coeff for tke wrt that for momentum

  • LMDZ6/trunk/libf/phylmd/call_atke_mod.F90

    r4644 r4653  
    99
    1010subroutine call_atke(dtime,ngrid,nlay,cdrag_uv,cdrag_t,u_surf,v_surf,temp_surf, &
    11                         wind_u,wind_v,temp,play,pinterf, &
     11                        wind_u,wind_v,temp,qvap,play,pinterf, &
    1212                        tke,Km_out,Kh_out)
    1313
    … …  
    3838REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: wind_v   ! meridional velocity (m/s)
    3939REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: temp   ! temperature (K)
     40REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: qvap   ! specific humidity (kg/kg)
    4041REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: play   ! pressure (Pa)
    4142REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: pinterf   ! pressure at interfaces(Pa)
    … …  
    5556
    5657call atke_compute_km_kh(ngrid,nlay,dtime,&
    57                         wind_u,wind_v,temp,play,pinterf, cdrag_uv,&
     58                        wind_u,wind_v,temp,qvap,play,pinterf,cdrag_uv,&
    5859                        tke,Km_out,Kh_out)
    5960
    … …  
    7475
    7576   call atke_compute_km_kh(ngrid,nlay,dtime,&
    76                         wind_u_predict,wind_v_predict,temp,play,pinterf,cdrag_uv, &
     77                        wind_u_predict,wind_v_predict,temp,qvap,play,pinterf,cdrag_uv, &
    7778                        tke,Km_out,Kh_out)
    7879

  • LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90

    r4619 r4653  
    17041704        IF (iflag_pbl>=50) THEN
    17051705
    1706         CALL call_atke(dtime,knon,klev,ycdragm, ycdragh,yus0,yvs0,yts,yu,yv,yt, &
     1706        CALL call_atke(dtime,knon,klev,ycdragm, ycdragh,yus0,yvs0,yts,yu,yv,yt,yq, &
    17071707             ypplay,ypaprs,ytke,ycoefm, ycoefh)
    17081708
    … …  
    17481748        IF (iflag_pbl>=50) THEN
    17491749     
    1750         CALL call_atke(dtime,knon,klev,ycdragm_x,ycdragh_x,yus0,yvs0,yts_x,yu_x,yv_x,yt_x, &
     1750        CALL call_atke(dtime,knon,klev,ycdragm_x,ycdragh_x,yus0,yvs0,yts_x,yu_x,yv_x,yt_x,yq_x, &
    17511751             ypplay,ypaprs,ytke_x,ycoefm_x, ycoefh_x)
    17521752
    … …  
    17871787        IF (iflag_pbl>=50) THEN
    17881788       
    1789         CALL call_atke(dtime,knon,klev,ycdragm_w,ycdragh_w,yus0,yvs0,yts_w,yu_w,yv_w,yt_w, &
     1789        CALL call_atke(dtime,knon,klev,ycdragm_w,ycdragh_w,yus0,yvs0,yts_w,yu_w,yv_w,yt_w,yq_w, &
    17901790             ypplay,ypaprs,ytke_w,ycoefm_w, ycoefh_w)
    17911791

  • LMDZ6/trunk/libf/phylmd/physiq_mod.F90

    r4651 r4653  
    18051805       CALL wake_ini(rg,rd,rv,prt_level)
    18061806       CALL yamada_ini(klon,lunout,prt_level)
    1807        CALL atke_ini(prt_level, lunout, RG, RD, RPI, RCPD)
     1807       CALL atke_ini(prt_level, lunout, RG, RD, RPI, RCPD, RV)
    18081808       CALL thermcell_ini(iflag_thermals,prt_level,tau_thermals,lunout, &
    18091809   &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)

  • LMDZ6/trunk/libf/phylmdiso/pbl_surface_mod.F90

    r4633 r4653  
    19811981        IF (iflag_pbl>=50) THEN
    19821982
    1983         CALL call_atke(dtime,knon,klev,ycdragm, ycdragh,yus0,yvs0,yts,yu,yv,yt, &
     1983        CALL call_atke(dtime,knon,klev,ycdragm, ycdragh,yus0,yvs0,yts,yu,yv,yt,yq, &
    19841984             ypplay,ypaprs,ytke,ycoefm, ycoefh)
    19851985
    … …  
    20242024        IF (iflag_pbl>=50) THEN
    20252025
    2026         CALL call_atke(dtime,knon,klev,ycdragm_x,ycdragh_x,yus0,yvs0,yts_x,yu_x,yv_x,yt_x, &
     2026        CALL call_atke(dtime,knon,klev,ycdragm_x,ycdragh_x,yus0,yvs0,yts_x,yu_x,yv_x,yt_x,yq_x, &
    20272027             ypplay,ypaprs,ytke_x,ycoefm_x, ycoefh_x)
    20282028
    … …  
    20642064        IF (iflag_pbl>=50) THEN
    20652065
    2066         CALL call_atke(dtime,knon,klev,ycdragm_w,ycdragh_w,yus0,yvs0,yts_w,yu_w,yv_w,yt_w, &
     2066        CALL call_atke(dtime,knon,klev,ycdragm_w,ycdragh_w,yus0,yvs0,yts_w,yu_w,yv_w,yt_w,yq_w, &
    20672067             ypplay,ypaprs,ytke_w,ycoefm_w, ycoefh_w)
    20682068

  • LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90

    r4651 r4653  
    19171917       CALL wake_ini(rg,rd,rv,prt_level)
    19181918       CALL yamada_ini(klon,lunout,prt_level)
    1919        CALL atke_ini(prt_level, lunout, RG, RD, RPI, RCPD)
     1919       CALL atke_ini(prt_level, lunout, RG, RD, RPI, RCPD, RV)
    19201920       CALL thermcell_ini(iflag_thermals,prt_level,tau_thermals,lunout, &
    19211921   &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)
Note: See TracChangeset for help on using the changeset viewer.