Changeset 4780

Timestamp:

Dec 21, 2023, 10:17:25 PM (11 months ago)

Author:

evignon

Message:

petites modifs dans les routines atke pour etre en accord avec les notations du papier

Location:

LMDZ6/trunk/libf/phylmd

Files:

• lmdz_atke_exchange_coeff.F90 (modified) (6 diffs)
• lmdz_atke_turbulence_ini.F90 (modified) (6 diffs)

LMDZ6/trunk/libf/phylmd/lmdz_atke_exchange_coeff.F90

@@ -25 +25 @@
 
 
-USE lmdz_atke_turbulence_ini, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd, atke_ok_virtual
+USE lmdz_atke_turbulence_ini, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd, atke_ok_virtual, ri0, ri1
 USE lmdz_atke_turbulence_ini, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, ctkes, rg, rd, rv, atke_ok_vdiff
-USE lmdz_atke_turbulence_ini, ONLY : viscom, viscoh, clmix, clmixshear, iflag_atke_lmix, lmin, smmin
+USE lmdz_atke_turbulence_ini, ONLY : viscom, viscoh, clmix, clmixshear, iflag_atke_lmix, lmin, smmin, cn
 
 implicit none
@@ -69 +69 @@
 
 INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid)
-REAL    :: cn,Ri0,Ri1    ! parameter for Sm stability function and Prandlt
 REAL    :: preff      ! reference pressure for potential temperature calculations
 REAL    :: thetam     ! mean potential temperature at interface
@@ -132 +131 @@
 ! Computes the gradient Richardson's number and stability functions:
 !===================================================================
-
-! calculation of cn = Sm value at Ri=0
-! direct dependance on cepsilon to guarantee Fm=1 (first-order like stability function) at Ri=0
-cn=(1./sqrt(cepsilon))**(2/3)
-! calculation of Ri0 such that continuity in slope of Sm at Ri=0
-Ri0=2./rpi*(cinf - cn)*ric/cn
-! calculation of Ri1 to guarantee continuity in slope of Prandlt number at Ri=0
-Ri1 = -2./rpi * (pr_asym - pr_neut) 
-
 
 DO ilay=2,nlay
@@ -176 +166 @@
 
 IF (iflag_atke_lmix .EQ. 1 ) THEN
-
+! Blackadar formulation (~kappa l) + buoyancy length scale (Deardoff 1980) for very stable conditions
    DO ilay=2,nlay
       DO igrid=1,ngrid
@@ -223 +213 @@
 
 ELSE
-! default: neglect effect of local stratification and shear
+! default Blackadar formulation: neglect effect of local stratification and shear
 
    DO ilay=2,nlay+1
@@ -250 +240 @@
 
 ! full implicit scheme resolved with a second order polynomial equation
-
+! default solution which shows fair convergence properties
     DO ilay=2,nlay
         DO igrid=1,ngrid

LMDZ6/trunk/libf/phylmd/lmdz_atke_turbulence_ini.F90

@@ -10 +10 @@
   real :: kappa = 0.4 ! Von Karman constant
   !$OMP THREADPRIVATE(kappa)
-  real :: l0,ric,ri0,cinf,cepsilon,pr_slope,pr_asym,pr_neut,clmix,clmixshear,smmin,ctkes,cke
-  !$OMP THREADPRIVATE(l0,ric,cinf,cepsilon,pr_slope,pr_asym,pr_neut,clmix,clmixshear,smmin,ctkes,cke)
+  real :: cinffac
+  !$OMP THREADPRIVATE(cinffac)
+  real :: l0,ric,ri0,ri1,cinf,cn,cepsilon,pr_slope,pr_asym,pr_neut,clmix,clmixshear,smmin,ctkes,cke
+  !$OMP THREADPRIVATE(l0,ri0,ri1,ric,cinf,cn,cepsilon,pr_slope,pr_asym,pr_neut,clmix,clmixshear,smmin,ctkes,cke)
   integer :: lunout,prt_level
   !$OMP THREADPRIVATE(lunout,prt_level)
@@ -18 +20 @@
   real :: viscom, viscoh
   !$OMP THREADPRIVATE(viscom,viscoh)
-  real :: lmin=0.01              ! minimum mixing length
+  real :: lmin=0.01              ! minimum mixing length corresponding to the Kolmogov dissipation scale 
+                                 ! in planetary atmospheres (Chen et al 2016, JGR Atmos)
   !$OMP THREADPRIVATE(lmin)
   logical :: atke_ok_vdiff, atke_ok_virtual
@@ -89 +92 @@
   ! Sun et al 2011, JAMC 
   ! between 10 and 40
-
   l0=15.0 
   CALL getin_p('atke_l0',l0)
@@ -97 +99 @@
   CALL getin_p('atke_ric',ric)
 
-  ! asymptotic value of Sm for Ri=-Inf
-  cinf=1.5
-  CALL getin_p('atke_cinf',cinf)
 
   ! constant for tke dissipation calculation
   cepsilon=5.87 ! default value as in yamada4
   CALL getin_p('atke_cepsilon',cepsilon)
+
+
+  ! calculation of cn = Sm value at Ri=0
+  ! direct dependance on cepsilon to guarantee Fm=1 (first-order like stability function) at Ri=0
+  cn=(1./sqrt(cepsilon))**(2/3)
+
+  ! asymptotic value of Sm for Ri=-Inf
+  cinffac=2.0
+  CALL getin_p('atke_cinffac',cinffac)
+  cinf=cinffac*cn
+  if (cinf .le. cn) then
+        call abort_physic("atke_turbulence_ini", &
+        'cinf cannot be lower than cn', 1)
+  endif
 
 
@@ -119 +132 @@
   endif
 
-  ! asymptotic turbulent prandt number value for Ri=-Inf
-  pr_asym=0.4
-  CALL getin_p('atke_pr_asym',pr_asym)
-
   ! value of turbulent prandtl number in neutral conditions (Ri=0)
   pr_neut=0.8
   CALL getin_p('atke_pr_neut',pr_neut)
+
+
+ ! asymptotic turbulent prandt number value for Ri=-Inf
+  pr_asym=0.4
+  CALL getin_p('atke_pr_asym',pr_asym)
+  if (pr_asym .le. pr_neut) then
+        call abort_physic("atke_turbulence_ini", &
+        'pr_asym must be be greater than pr_neut', 1)
+  endif
+
+
 
   ! coefficient for mixing length depending on local stratification
@@ -147 +167 @@
   CALL getin_p('atke_cke',cke)
 
+
+  ! calculation of Ri0 such that continuity in slope of Sm at Ri=0
+  ri0=2./rpi*(cinf - cn)*ric/cn
+
+  ! calculation of Ri1 to guarantee continuity in slope of Prandlt number at Ri=0
+  ri1 = -2./rpi * (pr_asym - pr_neut)
+
+
  RETURN