Ignore:
Timestamp:
Feb 7, 2024, 4:07:16 PM (8 months ago)
Author:
evignon
Message:

travail de documentation et commentaire des routines ATKE par Clement Dehondt

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r4781 r4804  
    44
    55! declaration of constants and parameters
    6 save
    76
    8 integer :: iflag_atke, iflag_num_atke, iflag_atke_lmix
    9   !$OMP THREADPRIVATE(iflag_atke, iflag_num_atke, iflag_atke_lmix)
    10   real :: kappa = 0.4 ! Von Karman constant
    11   !$OMP THREADPRIVATE(kappa)
    12   real :: cinffac
    13   !$OMP THREADPRIVATE(cinffac)
    14   real :: l0,ric,ri0,ri1,cinf,cn,cepsilon,pr_slope,pr_asym,pr_neut,clmix,clmixshear,smmin,ctkes,cke
    15   !$OMP THREADPRIVATE(l0,ri0,ri1,ric,cinf,cn,cepsilon,pr_slope,pr_asym,pr_neut,clmix,clmixshear,smmin,ctkes,cke)
    16   integer :: lunout,prt_level
    17   !$OMP THREADPRIVATE(lunout,prt_level)
    18   real :: rg, rd, rpi, rcpd, rv
    19   !$OMP THREADPRIVATE(rg, rd, rpi, rcpd, rv)
    20   real :: viscom, viscoh
    21   !$OMP THREADPRIVATE(viscom,viscoh)
    22   real :: lmin=0.01              ! minimum mixing length corresponding to the Kolmogov dissipation scale
    23                                  ! in planetary atmospheres (Chen et al 2016, JGR Atmos)
    24   !$OMP THREADPRIVATE(lmin)
    25   logical :: atke_ok_vdiff, atke_ok_virtual
    26   !$OMP THREADPRIVATE(atke_ok_vdiff,atke_ok_virtual)
     7real, save, protected  :: rg, rd, rpi, rcpd, rv, viscom, viscoh       
     8!$OMP THREADPRIVATE(rg, rd, rpi, rcpd, rv, viscom, viscoh)     
     9
     10integer, save, protected :: iflag_atke        ! flag that controls options in atke_compute_km_kh
     11integer, save, protected :: iflag_atke_lmix   ! flag that controls the calculation of mixing length in atke
     12integer, save, protected :: iflag_num_atke    ! flag that controls the numerical treatment of diffusion coeffiient calculation
     13!$OMP THREADPRIVATE(iflag_atke, iflag_atke_lmix, iflag_num_atke)
     14
     15logical, save, protected :: atke_ok_vdiff    ! activate vertical diffusion of TKE or not
     16logical, save, protected :: atke_ok_virtual  ! account for vapor for flottability
     17!$OMP THREADPRIVATE(atke_ok_vdiff, atke_ok_virtual)
     18
     19real, save, protected :: kappa = 0.4         ! Von Karman constant
     20real, save, protected :: cn                  ! Sm value at Ri=0
     21real, save, protected :: cinf                ! Sm value at Ri=-Inf
     22real, save, protected :: ri0                 ! Richardson number near zero to guarantee continuity in slope of Sm (stability function) at Ri=0
     23real, save, protected :: ri1                 ! Richardson number near zero to guarantee continuity in slope of Pr (Prandlt's number) at Ri=0
     24real, save, protected :: lmin = 0.01         ! minimum mixing length corresponding to the Kolmogov dissipation scale  in planetary atmospheres (Chen et al 2016, JGR Atmos)
     25real, save, protected :: ctkes               ! coefficient for surface TKE
     26real, save, protected :: clmixshear          ! coefficient for mixing length depending on local wind shear
     27!$OMP THREADPRIVATE(kappa, cn, cinf, ri0, ri1, lmin, ctkes, clmixshear)
     28
     29
     30! Tunable parameters for the ATKE scheme and their range of values
     31!!-------------------------------------------------------------------------------------------------------------
     32real, save, protected :: cepsilon            ! controls the value of the dissipation length scale, range [1.2 - 10]
     33real, save, protected :: cke                 ! controls the value of the diffusion coefficient of TKE, range [1 - 5]
     34real, save, protected :: l0                  ! asymptotic mixing length far from the ground [m] (Sun et al 2011, JAMC), range [15 - 75]
     35real, save, protected :: clmix               ! controls the value of the mixing length in stratified conditions, range [0.1 - 2]
     36real, save, protected :: ric                 ! critical Richardson number controlling the slope of Sm in stable conditions, range [0.19 - 0.25]
     37real, save, protected :: smmin               ! minimum value of Sm in very stable conditions (defined here as minsqrt(Ez/Ek)) at large Ri, range [0.025 - 0.1]
     38real, save, protected :: pr_neut             ! neutral value of the Prandtl number (Ri=0), range [0.7 - 1]
     39real, save, protected :: pr_slope            ! linear slope of Pr with Ri in the very stable regime, range [3 - 5]
     40real, save, protected :: cinffac             ! ratio between cinf and cn controlling the convective limit of Sm, range [1.2 - 5.0]
     41real, save, protected :: pr_asym             ! value of Prandlt in the convective limit(Ri=-Inf), range [0.3 - 0.5]
     42!$OMP THREADPRIVATE(cepsilon, cke, l0, clmix, ric, smmin, pr_neut, pr_slope, cinffac, pr_asym)
     43!!-------------------------------------------------------------------------------------------------------------
    2744
    2845CONTAINS
    2946
    3047SUBROUTINE atke_ini(rg_in, rd_in, rpi_in, rcpd_in, rv_in, viscom_in, viscoh_in)
     48      !!----------------------------------------------------------------------
     49      !!                  ***  ROUTINE  atke_ini  ***
     50      !!
     51      !! ** Purpose :   Initialization of the atke module and choice of some constants
     52      !!               
     53      !!----------------------------------------------------------------------
    3154
    32    USE ioipsl_getin_p_mod, ONLY : getin_p
     55        USE ioipsl_getin_p_mod, ONLY : getin_p                                                                                                     
    3356
    34   real, intent(in) :: rg_in, rd_in, rpi_in, rcpd_in, rv_in, viscom_in, viscoh_in
     57      ! input arguments (universal constants for planet)
     58      !-------------------------------------------------
     59      real, intent(in) :: rg_in, rd_in, rpi_in, rcpd_in, rv_in, viscom_in, viscoh_in
     60      !!----------------------------------------------------------------------
     61     
     62      rg=rg_in          ! gravity acceleration
     63      rd=rd_in          ! dry gas constant (R/M, R=perfect gas constant and M is the molar mass of the fluid)
     64      rpi=rpi_in        ! Pi number
     65      rcpd=rcpd_in      ! cp per unit mass of the gas
     66      rv=rv_in          ! water vapor constant (for simulations in Earth's atmosphere)
     67      viscom=viscom_in  ! kinematic molecular viscosity for momentum
     68      viscoh=viscoh_in  ! kinematic molecular viscosity for heat
    3569
    3670
    37   ! input arguments (universal constants for planet)
    38   !-------------------------------------------------
    39  
    40   ! gravity acceleration
    41   rg=rg_in
    42   ! dry gas constant (R/M, R=perfect gas constant and M is the molar mass of the fluid)
    43   rd=rd_in
    44   ! Pi number
    45   rpi=rpi_in
    46   ! cp per unit mass of the gas
    47   rcpd=rcpd_in
    48   ! water vapor constant (for simulations in Earth's atmosphere)
    49   rv=rv_in
    50   ! kinematic molecular viscosity for momentum
    51   viscom=viscom_in
    52   ! kinematic molecular viscosity for heat
    53   viscoh=viscoh_in
     71      ! Read flag values in .def files
     72      !-------------------------------
     73
     74      ! flag that controls options in atke_compute_km_kh
     75      iflag_atke=0
     76      CALL getin_p('iflag_atke',iflag_atke)
     77
     78      ! flag that controls the calculation of mixing length in atke
     79      iflag_atke_lmix=0
     80      CALL getin_p('iflag_atke_lmix',iflag_atke_lmix)
     81
     82      if (iflag_atke .eq. 0 .and. iflag_atke_lmix>0) then
     83            call abort_physic("atke_turbulence_ini", &
     84            'stationary scheme must use mixing length formulation not depending on tke', 1)
     85      endif
     86
     87      ! activate vertical diffusion of TKE or not
     88      atke_ok_vdiff=.false.
     89      CALL getin_p('atke_ok_vdiff',atke_ok_vdiff)
     90
     91      ! account for vapor for flottability
     92      atke_ok_virtual=.true.
     93      CALL getin_p('atke_ok_virtual',atke_ok_virtual)
    5494
    5595
    56   !viscom=1.46E-5
    57   !viscoh=2.06E-5
     96      ! flag that controls the numerical treatment of diffusion coeffiient calculation
     97      iflag_num_atke=0
     98      CALL getin_p('iflag_num_atke',iflag_num_atke)
    5899
     100      ! asymptotic mixing length in neutral conditions [m]
     101      ! Sun et al 2011, JAMC
     102      ! between 10 and 40
     103      l0=15.0
     104      CALL getin_p('atke_l0',l0)
    59105
    60   ! Read flag values in .def files
    61   !-------------------------------
     106      ! critical Richardson number
     107      ric=0.25
     108      CALL getin_p('atke_ric',ric)
    62109
     110      ! constant for tke dissipation calculation
     111      cepsilon=5.87 ! default value as in yamada4
     112      CALL getin_p('atke_cepsilon',cepsilon)
    63113
    64   ! flag that controls options in atke_compute_km_kh
    65   iflag_atke=0
    66   CALL getin_p('iflag_atke',iflag_atke)
     114      ! calculation of cn = Sm value at Ri=0
     115      ! direct dependance on cepsilon to guarantee Fm=1 (first-order like stability function) at Ri=0
     116      cn=(1./sqrt(cepsilon))**(2/3)
    67117
    68   ! flag that controls the calculation of mixing length in atke
    69   iflag_atke_lmix=0
    70   CALL getin_p('iflag_atke_lmix',iflag_atke_lmix)
     118      ! asymptotic value of Sm for Ri=-Inf
     119      cinffac=2.0
     120      CALL getin_p('atke_cinffac',cinffac)
     121      cinf=cinffac*cn
     122      if (cinf .le. cn) then
     123            call abort_physic("atke_turbulence_ini", &
     124            'cinf cannot be lower than cn', 1)
     125      endif
    71126
    72   if (iflag_atke .eq. 0 .and. iflag_atke_lmix>0) then
    73         call abort_physic("atke_turbulence_ini", &
    74         'stationary scheme must use mixing length formulation not depending on tke', 1)
    75   endif
     127      ! coefficient for surface TKE
     128      ! following Lenderink & Holtslag 2004, ctkes=(cepsilon)**(2/3)
     129      ! (provided by limit condition in neutral conditions)
     130      ctkes=(cepsilon)**(2./3.)
    76131
    77   ! activate vertical diffusion of TKE or not
    78   atke_ok_vdiff=.false.
    79   CALL getin_p('atke_ok_vdiff',atke_ok_vdiff)
     132      ! slope of Pr=f(Ri) for stable conditions
     133      pr_slope=5.0 ! default value from Zilitinkevich et al. 2005
     134      CALL getin_p('atke_pr_slope',pr_slope)
     135      if (pr_slope .le. 1) then
     136            call abort_physic("atke_turbulence_ini", &
     137            'pr_slope has to be greater than 1 for consistency of the tke scheme', 1)
     138      endif
    80139
     140      ! value of turbulent prandtl number in neutral conditions (Ri=0)
     141      pr_neut=0.8
     142      CALL getin_p('atke_pr_neut',pr_neut)
    81143
    82   ! account for vapor for flottability
    83   atke_ok_virtual=.true.
    84   CALL getin_p('atke_ok_virtual',atke_ok_virtual)
     144      ! asymptotic turbulent prandt number value for Ri=-Inf
     145      pr_asym=0.4
     146      CALL getin_p('atke_pr_asym',pr_asym)
     147      if (pr_asym .ge. pr_neut) then
     148            call abort_physic("atke_turbulence_ini", &
     149            'pr_asym must be be lower than pr_neut', 1)
     150      endif
    85151
     152      ! coefficient for mixing length depending on local stratification
     153      clmix=0.5
     154      CALL getin_p('atke_clmix',clmix)
    86155
    87   ! flag that controls the numerical treatment of diffusion coeffiient calculation
    88   iflag_num_atke=0
    89   CALL getin_p('iflag_num_atke',iflag_num_atke)
     156      ! coefficient for mixing length depending on local wind shear
     157      clmixshear=0.5
     158      CALL getin_p('atke_clmixshear',clmixshear)
    90159
    91   ! asymptotic mixing length in neutral conditions [m]
    92   ! Sun et al 2011, JAMC
    93   ! between 10 and 40
    94   l0=15.0
    95   CALL getin_p('atke_l0',l0)
     160      ! minimum anisotropy coefficient (defined here as minsqrt(Ez/Ek)) at large Ri.
     161      ! From Zilitinkevich et al. 2013, it equals sqrt(0.03)~0.17 
     162      smmin=0.17
     163      CALL getin_p('atke_smmin',smmin)
    96164
    97   ! critical Richardson number
    98   ric=0.25
    99   CALL getin_p('atke_ric',ric)
     165      ! ratio between the eddy diffusivity coeff for tke wrt that for momentum
     166      ! default value from Lenderink et al. 2004
     167      cke=2.
     168      CALL getin_p('atke_cke',cke)
    100169
     170      ! calculation of Ri0 such that continuity in slope of Sm at Ri=0
     171      ri0=2./rpi*(cinf - cn)*ric/cn
    101172
    102   ! constant for tke dissipation calculation
    103   cepsilon=5.87 ! default value as in yamada4
    104   CALL getin_p('atke_cepsilon',cepsilon)
     173      ! calculation of Ri1 to guarantee continuity in slope of Prandlt number at Ri=0
     174      ri1 = -2./rpi * (pr_asym - pr_neut)
    105175
    106 
    107   ! calculation of cn = Sm value at Ri=0
    108   ! direct dependance on cepsilon to guarantee Fm=1 (first-order like stability function) at Ri=0
    109   cn=(1./sqrt(cepsilon))**(2/3)
    110 
    111   ! asymptotic value of Sm for Ri=-Inf
    112   cinffac=2.0
    113   CALL getin_p('atke_cinffac',cinffac)
    114   cinf=cinffac*cn
    115   if (cinf .le. cn) then
    116         call abort_physic("atke_turbulence_ini", &
    117         'cinf cannot be lower than cn', 1)
    118   endif
    119 
    120 
    121  ! coefficient for surface TKE
    122  ! following Lenderink & Holtslag 2004, ctkes=(cepsilon)**(2/3)
    123  ! (provided by limit condition in neutral conditions)
    124   ctkes=(cepsilon)**(2./3.)
    125 
    126   ! slope of Pr=f(Ri) for stable conditions
    127   pr_slope=5.0 ! default value from Zilitinkevich et al. 2005
    128   CALL getin_p('atke_pr_slope',pr_slope)
    129   if (pr_slope .le. 1) then
    130         call abort_physic("atke_turbulence_ini", &
    131         'pr_slope has to be greater than 1 for consistency of the tke scheme', 1)
    132   endif
    133 
    134   ! value of turbulent prandtl number in neutral conditions (Ri=0)
    135   pr_neut=0.8
    136   CALL getin_p('atke_pr_neut',pr_neut)
    137 
    138 
    139  ! asymptotic turbulent prandt number value for Ri=-Inf
    140   pr_asym=0.4
    141   CALL getin_p('atke_pr_asym',pr_asym)
    142   if (pr_asym .ge. pr_neut) then
    143         call abort_physic("atke_turbulence_ini", &
    144         'pr_asym must be be lower than pr_neut', 1)
    145   endif
    146 
    147 
    148 
    149   ! coefficient for mixing length depending on local stratification
    150   clmix=0.5
    151   CALL getin_p('atke_clmix',clmix)
    152  
    153   ! coefficient for mixing length depending on local wind shear
    154   clmixshear=0.5
    155   CALL getin_p('atke_clmixshear',clmixshear)
    156 
    157 
    158   ! minimum anisotropy coefficient (defined here as minsqrt(Ez/Ek)) at large Ri.
    159   ! From Zilitinkevich et al. 2013, it equals sqrt(0.03)~0.17 
    160   smmin=0.17
    161   CALL getin_p('atke_smmin',smmin)
    162 
    163 
    164   ! ratio between the eddy diffusivity coeff for tke wrt that for momentum
    165   ! default value from Lenderink et al. 2004
    166   cke=2.
    167   CALL getin_p('atke_cke',cke)
    168 
    169 
    170   ! calculation of Ri0 such that continuity in slope of Sm at Ri=0
    171   ri0=2./rpi*(cinf - cn)*ric/cn
    172 
    173   ! calculation of Ri1 to guarantee continuity in slope of Prandlt number at Ri=0
    174   ri1 = -2./rpi * (pr_asym - pr_neut)
    175 
    176 
    177  RETURN
     176RETURN
    178177
    179178END SUBROUTINE atke_ini
Note: See TracChangeset for help on using the changeset viewer.