Changeset 4631 for LMDZ6


Ignore:
Timestamp:
Jul 18, 2023, 11:36:49 AM (17 months ago)
Author:
evignon
Message:

commission pour l'atelier TKE avec introduction d'une resolution pronostique
de lequation de la tke

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

Legend:

Unmodified
Added
Removed

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

    r4481 r4631  
    55contains
    66
    7 subroutine atke_compute_km_kh(ngrid,nlay, &
     7subroutine atke_compute_km_kh(ngrid,nlay,dtime, &
    88                        wind_u,wind_v,temp,play,pinterf, &
    99                        tke,Km_out,Kh_out)
    … …  
    2828USE atke_turbulence_ini_mod, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd
    2929USE atke_turbulence_ini_mod, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, rg, rd
    30 USE atke_turbulence_ini_mod, ONLY : viscom, viscoh
     30USE atke_turbulence_ini_mod, ONLY : viscom, viscoh, clmix, iflag_atke_lmix, lmin
    3131
    3232implicit none
    … …  
    3737
    3838INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
    39 INTEGER, INTENT(IN) :: nlay ! number of vertical index 
    40 
     39INTEGER, INTENT(IN) :: nlay  ! number of vertical index 
     40
     41REAL, INTENT(IN)    :: dtime ! physics time step (s)
    4142REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: wind_u   ! zonal velocity (m/s)
    4243REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: wind_v   ! meridional velocity (m/s)
    … …  
    6061REAL, DIMENSION(ngrid,nlay)   :: dz_interf   ! distance between two consecutive interfaces
    6162REAL, DIMENSION(ngrid,nlay)   :: dz_lay      ! distance between two layer middles (NB: first and last are half layers)
     63REAL, DIMENSION(ngrid,nlay+1) :: N2          ! square of Brunt Vaisala pulsation (at interface)
     64REAL, DIMENSION(ngrid,nlay+1) :: shear2      ! square of wind shear (at interface)
    6265REAL, DIMENSION(ngrid,nlay+1) :: Ri          ! Richardson's number (at interface)
    6366REAL, DIMENSION(ngrid,nlay+1) :: Prandtl     ! Turbulent Prandtl's number (at interface)
    6467REAL, DIMENSION(ngrid,nlay+1) :: Sm          ! Stability function for momentum (at interface)
    6568REAL, DIMENSION(ngrid,nlay+1) :: Sh          ! Stability function for heat (at interface)
     69LOGICAL, DIMENSION(ngrid,nlay+1) :: switch_num ! switch of numerical integration in very stable cases
    6670
    6771INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid)
    … …  
    6973REAL    :: preff      ! reference pressure for potential temperature calculations
    7074REAL    :: thetam     ! mean potential temperature at interface
    71 
     75REAL    :: delta      ! discriminant of the second order polynomial
     76REAL    :: qq         ! tke=qq**2/2
     77REAL    :: shear      ! wind shear
     78REAL    :: lstrat     ! mixing length depending on local stratification
     79REAL    :: taustrat   ! caracteristic timescale for turbulence in very stable conditions
    7280
    7381! Initializations:
    … …  
    109117
    110118
    111 ! Computing the mixing length:
    112 ! so far, we have neglected the effect of local stratification
    113 !==============================================================
    114 
    115 
    116 DO ilay=2,nlay+1
    117     DO igrid=1,ngrid
    118         l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
    119     ENDDO
    120 ENDDO
    121 
    122119! Computes the gradient Richardson's number and stability functions:
    123120!===================================================================
    … …  
    136133        dz_lay(igrid,ilay)=z_lay(igrid,ilay)-z_lay(igrid,ilay-1)
    137134        thetam=0.5*(theta(igrid,ilay) + theta(igrid,ilay-1))
    138         Ri(igrid,ilay) = rg * (theta(igrid,ilay) - theta(igrid,ilay-1))/thetam / dz_lay(igrid,ilay) / &
    139         MAX(((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + &
    140         ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2,1E-10)
     135        N2(igrid,ilay) = rg * (theta(igrid,ilay) - theta(igrid,ilay-1))/thetam / dz_lay(igrid,ilay)
     136        shear2(igrid,ilay)= (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + &
     137            ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 )
     138        Ri(igrid,ilay) = N2(igrid,ilay) / MAX(shear2(igrid,ilay),1E-10)
    141139       
    142140        IF (Ri(igrid,ilay) < 0.) THEN ! unstable cases
    … …  
    157155ENDDO
    158156
     157
     158! Computing the mixing length:
     159!==============================================================
     160
     161switch_num(:,:)=.false.
     162
     163IF (iflag_atke_lmix .EQ. 1 ) THEN
     164
     165   DO ilay=2,nlay
     166      DO igrid=1,ngrid
     167          l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
     168          IF (N2(igrid,ilay) .GT. 0.) THEN
     169             lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay))
     170             IF (lstrat .LT. l_exchange(igrid,ilay)) THEN
     171                l_exchange(igrid,ilay)=max(lstrat,lmin)
     172                switch_num(igrid,ilay)=.true.
     173             ENDIF
     174          ENDIF
     175      ENDDO
     176   ENDDO
     177
     178ELSE
     179! default: neglect effect of local stratification
     180
     181   DO ilay=2,nlay+1
     182      DO igrid=1,ngrid
     183          l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
     184      ENDDO
     185
     186   ENDDO
     187ENDIF
     188
     189
    159190! Computing the TKE:
    160191!===================
    161192IF (iflag_atke == 0) THEN
    162193
    163 ! stationary solution neglecting the vertical transport of TKE by turbulence
     194   DO ilay=2,nlay
     195        DO igrid=1,ngrid
     196            tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * &
     197            shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))
     198        ENDDO
     199    ENDDO
     200
     201ELSE IF (iflag_atke == 1) THEN
     202
     203! full implicit scheme resolved with a second order polynomial equation
     204
    164205    DO ilay=2,nlay
    165206        DO igrid=1,ngrid
    166             tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * &
    167             (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + &
    168             ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 ) * &
    169             (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))
    170         ENDDO
    171     ENDDO
    172 
    173 ELSE ! TO DO
    174    
     207           qq=sqrt(2.*tke(igrid,ilay))
     208           delta=1.+4.*dtime/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2)) * &
     209                (qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay)*shear2(igrid,ilay) &
     210                *(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)))
     211           qq=(-1. + sqrt(delta))/dtime*cepsilon*sqrt(2.)*l_exchange(igrid,ilay)
     212           tke(igrid,ilay)=0.5*(qq**2)
     213        ENDDO
     214    ENDDO
     215
     216ELSE IF (iflag_atke == 2) THEN
     217
     218! full implicit scheme resolved with a second order polynomial equation
     219! + exact resolution for very stable cases
     220    DO ilay=2,nlay
     221        DO igrid=1,ngrid
     222           qq=sqrt(2.*tke(igrid,ilay))
     223           IF (switch_num(igrid,ilay)) THEN
     224           taustrat=clmix/2./sqrt(N2(igrid,ilay))*Sm(igrid,ilay)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) &
     225                    - sqrt(N2(igrid,ilay))/2./cepsilon/clmix
     226           taustrat=-1./taustrat
     227           qq=qq*exp(-dtime/taustrat)
     228           ELSE
     229           delta=1.+4.*dtime/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2)) * &
     230                (qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay) * shear2(igrid,ilay) &
     231                *(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)))
     232           qq=(-1. + sqrt(delta))/dtime*cepsilon*sqrt(2.)*l_exchange(igrid,ilay)
     233           ENDIF
     234           tke(igrid,ilay)=0.5*(qq**2)
     235        ENDDO
     236    ENDDO
     237
     238
     239
     240ELSE IF (iflag_atke == 3) THEN
     241
     242! semi implicit scheme when l does not depend on tke
     243! positive-guaranteed if pr slope in stable condition >1
     244    DO ilay=2,nlay
     245        DO igrid=1,ngrid
     246           qq=sqrt(2.*tke(igrid,ilay))
     247           qq=(qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) &
     248             /(1.+dtime*qq/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2)))     
     249           tke(igrid,ilay)=0.5*(qq**2)
     250        ENDDO
     251    ENDDO
     252
     253
     254ELSE
    175255   call abort_physic("atke_compute_km_kh", &
    176         'traitement non-stationnaire de la tke pas encore prevu', 1)
     256        'numerical treatment of TKE not possible yet', 1)
    177257
    178258END IF

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

    r4545 r4631  
    55save
    66
    7 integer :: iflag_atke, iflag_num_atke
    8 !$OMP THREADPRIVATE(iflag_atke, iflag_num_atke)
     7integer :: iflag_atke, iflag_num_atke, iflag_atke_lmix
     8  !$OMP THREADPRIVATE(iflag_atke, iflag_num_atke, iflag_atke_lmix)
    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
    12   !$OMP THREADPRIVATE(l0, ric, cinf, cepsilon, pr_slope, pr_asym, pr_neut)
     11  real :: l0, ric, ri0, cinf, cepsilon, pr_slope, pr_asym, pr_neut, clmix
     12  !$OMP THREADPRIVATE(l0, ric, cinf, cepsilon, pr_slope, pr_asym, pr_neut, clmix)
    1313  integer :: lunout,prt_level
    1414  !$OMP THREADPRIVATE(lunout,prt_level)
    … …  
    1919  !$OMP THREADPRIVATE(viscom,viscoh)
    2020
     21  real :: lmin=0.001              ! minimum mixing length
     22  !$OMP THREADPRIVATE(lmin)
    2123
    2224
    … …  
    4547  CALL getin_p('iflag_atke',iflag_atke)
    4648
     49  ! flag that controls the calculation of mixing length in atke
     50  iflag_atke_lmix=0
     51  CALL getin_p('iflag_atke_lmix',iflag_atke_lmix)
     52
     53  if (iflag_atke .eq. 0 .and. iflag_atke_lmix>0) then
     54        call abort_physic("atke_turbulence_ini", &
     55        'stationary scheme must use mixing length formulation not depending on tke', 1)
     56  endif
     57
     58
    4759  ! flag that controls the numerical treatment of diffusion coeffiient calculation
    4860  iflag_num_atke=0
    4961  CALL getin_p('iflag_num_atke',iflag_num_atke)
    5062
    51   ! asymptotic mixing length [m]
    52   l0=150.0
     63  ! asymptotic mixing length in neutral conditions [m]
     64  ! Sun et al 2011, JAMC
     65  ! between 10 and 40
     66
     67  l0=15.0
    5368  CALL getin_p('atke_l0',l0)
    5469
    … …  
    6277
    6378  ! constant for tke dissipation calculation
    64   cepsilon=16.6 ! default value as in yamada4
     79  cepsilon=16.6/2./sqrt(2.) ! default value as in yamada4
    6580  CALL getin_p('atke_cepsilon',cepsilon)
    6681
    … …  
    6883  pr_slope=5.0 ! default value from Zilitinkevich et al. 2005
    6984  CALL getin_p('atke_pr_slope',pr_slope)
     85  if (pr_slope .le. 1) then
     86        call abort_physic("atke_turbulence_ini", &
     87        'pr_slope has to be greater than 1 for consistency of the tke scheme', 1)
     88  endif
    7089
    7190  ! asymptotic turbulent prandt number value for Ri=-Inf
    … …  
    7796  CALL getin_p('atke_pr_neut',pr_neut)
    7897
     98  ! coefficient for mixing length depending on local stratification
     99  clmix=0.2
     100  CALL getin_p('atke_clmix',clmix)
    79101   
    80102 RETURN

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

    r4545 r4631  
    5454
    5555
    56 call atke_compute_km_kh(ngrid,nlay, &
     56call atke_compute_km_kh(ngrid,nlay,dtime,&
    5757                        wind_u,wind_v,temp,play,pinterf, &
    5858                        tke,Km_out,Kh_out)
    … …  
    7070
    7171
    72    call atke_compute_km_kh(ngrid,nlay, &
     72   call atke_compute_km_kh(ngrid,nlay,dtime,&
    7373                        wind_u_predict,wind_v_predict,temp,play,pinterf, &
    7474                        tke,Km_out,Kh_out)
Note: See TracChangeset for help on using the changeset viewer.