Changeset 5821

Timestamp:

Sep 18, 2025, 10:29:26 AM (2 months ago)

Author:

idelkadi

Message:

Added a formulation to prescribe effective cloud size as a hyperbolic tangent function of pressure for calculating radiative fluxes related to 3D cloud effects.
Activation is controlled in namelist_ecrad by the logical key ok_separation_tanh

Location:

LMDZ6/trunk/libf/phylmd/ecrad

Files:

• lmdz/radiation_scheme_mod.f90 (modified) (2 diffs)
• lmdz/setup_config_from_lmdz.f90 (modified) (5 diffs)
• radiation/radiation_cloud.F90 (modified) (2 diffs)

LMDZ6/trunk/libf/phylmd/ecrad/lmdz/radiation_scheme_mod.f90

@@ -428 +428 @@
                &  driver_config%high_inv_effective_size, 0.8_jprb, 0.45_jprb, &
                &  KIDIA, KFDIA)
-   else if (driver_config%ok_separation) then
+   else if (driver_config%ok_separation_eta) then
      call cloud%param_cloud_effective_separation_eta(klon, klev, &
+               &  thermodynamics%pressure_hl, &
+               &  driver_config%cloud_separation_scale_surface, &
+               &  driver_config%cloud_separation_scale_toa, &
+               &  driver_config%cloud_separation_scale_power, &
+               &  driver_config%cloud_inhom_separation_factor, &
+               &  KIDIA, KFDIA)
+   else if (driver_config%ok_separation_tanh) then
+     call cloud%param_cloud_effective_separation_tanh(klon, klev, &
                &  thermodynamics%pressure_hl, &
                &  driver_config%cloud_separation_scale_surface, &
@@ -1134 +1142 @@
                &  driver_config%high_inv_effective_size, 0.8_jprb, 0.45_jprb, &
                &  KIDIA, KFDIA)
-   else if (driver_config%ok_separation) then
+   else if (driver_config%ok_separation_eta) then
      call cloud%param_cloud_effective_separation_eta(klon, klev, &
+               &  thermodynamics%pressure_hl, &
+               &  driver_config%cloud_separation_scale_surface, &
+               &  driver_config%cloud_separation_scale_toa, &
+               &  driver_config%cloud_separation_scale_power, &
+               &  driver_config%cloud_inhom_separation_factor, &
+               &  KIDIA, KFDIA)
+   else if (driver_config%ok_separation_tanh) then
+     call cloud%param_cloud_effective_separation_tanh(klon, klev, &
                &  thermodynamics%pressure_hl, &
                &  driver_config%cloud_separation_scale_surface, &

LMDZ6/trunk/libf/phylmd/ecrad/lmdz/setup_config_from_lmdz.f90

@@ -9 +9 @@
   type driver_config_type
      logical    :: ok_effective_size  = .true.
-     logical    :: ok_separation = .false.   
+     logical    :: ok_separation_eta = .false.   
+     logical    :: ok_separation_tanh = .false.
      real(jprb) :: high_inv_effective_size   = -1.0_jprb ! m-1
      real(jprb) :: middle_inv_effective_size = -1.0_jprb ! m-1
@@ -44 +45 @@
     character(*), intent(in)          :: file_name
     logical, intent(out), optional    :: is_success
-    logical    :: ok_effective_size, ok_separation
+    logical    :: ok_effective_size, ok_separation_eta, ok_separation_tanh
     integer    :: iosopen ! Status after calling open
     real(jprb) :: high_inv_effective_size 
@@ -61 +62 @@
     logical :: do_save_inputs
 
-    namelist /radiation_driver/ ok_effective_size, ok_separation, &
+    namelist /radiation_driver/ ok_effective_size, ok_separation_eta, ok_separation_tanh, &
          &  frac_std, overlap_decorr_length, kdecolat, &
          &  low_decorrelation_length, mid_decorrelation_length, high_decorrelation_length, &
@@ -70 +71 @@
 
     ok_effective_size = .false.
-    ok_separation = .false.
+    ok_separation_eta = .false.
+    ok_separation_tanh = .false.
     high_inv_effective_size   = -1.0_jprb 
     middle_inv_effective_size = -1.0_jprb 
@@ -109 +111 @@
     ! Copy namelist data into configuration object
     this%ok_effective_size = ok_effective_size 
-    this%ok_separation = ok_separation
+    this%ok_separation_eta = ok_separation_eta
+    this%ok_separation_tanh = ok_separation_tanh
     this%frac_std = frac_std
     this%overlap_decorr_length = overlap_decorr_length

LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_cloud.F90

@@ -91 +91 @@
     procedure :: create_inv_cloud_effective_size_eta
     procedure :: param_cloud_effective_separation_eta
+    procedure :: param_cloud_effective_separation_tanh
     procedure :: crop_cloud_fraction
     procedure :: out_of_physical_bounds
@@ -774 +775 @@
   end subroutine param_cloud_effective_separation_eta
 
-
+  !AI 20 mai 2025
+ !Nouvelle param Lx=f(P) avec Lf(P) tgte hyperb
+ subroutine param_cloud_effective_separation_tanh(this, ncol, nlev, &
+     &  pressure_hl, &
+     &  lx_surf, lx_toa, gama, &
+     &  inhom_separation_factor, istartcol, iendcol)
+
+  use yomhook,                  only : lhook, dr_hook, jphook
+  use write_field_phy
+  implicit none
+
+  class(cloud_type), intent(inout) :: this
+  integer,           intent(in)    :: ncol, nlev
+  ! Pressure on half levels (Pa) : size should be (ncol, nlev+1)
+  real(jprb),        intent(in)    :: pressure_hl(:,:)
+  ! Separation distances at surface and top-of-atmosphere (m)
+  real(jprb),        intent(in)    :: lx_surf, lx_toa
+  ! Facteur de raideur
+  real(jprb),        intent(in)    :: gama
+  ! Optional factor
+  real(jprb), optional, intent(in) :: inhom_separation_factor
+  integer,    optional, intent(in) :: istartcol, iendcol
+
+  ! Parameters (AI attention a revoir)
+  real(kind=jprb), parameter :: eps = 1.0e-5_jprb
+
+  ! Locals
+  real(jprb) :: eff_separation(ncol,nlev)
+  real(jprb) :: Lmin, Lmax, Pmin, Pmax, p_mean, p0, inhom_sep_factor, alpha
+  integer :: jcol, k
+  integer :: i1, i2
+  real(jphook) :: hook_handle
+  real(jprb) :: frac, denom1, denom2
+  real(jprb) :: pdiff
+
+  ! Hook entry
+  if (lhook) call dr_hook('radiation_cloud:param_cloud_effective_separation_eta',0,hook_handle)
+
+  ! inhom separation factor default
+  if (present(inhom_separation_factor)) then
+    inhom_sep_factor = inhom_separation_factor
+  else
+    inhom_sep_factor = 1.0_jprb
+  end if
+
+  ! (re)allocation of output arrays
+  if (allocated(this%inv_cloud_effective_size)) then
+    deallocate(this%inv_cloud_effective_size)
+  end if
+  if (allocated(this%inv_inhom_effective_size)) then
+    deallocate(this%inv_inhom_effective_size)
+  end if
+
+  allocate(this%inv_cloud_effective_size(ncol, nlev))
+  allocate(this%inv_inhom_effective_size(ncol, nlev))
+
+  ! columns range
+  if (present(istartcol)) then
+    i1 = istartcol
+  else
+    i1 = 1
+  end if
+
+  if (present(iendcol)) then
+    i2 = iendcol
+  else
+    i2 = ncol
+  end if
+
+  ! map inputs to local names
+  Lmin = lx_surf
+  Lmax = lx_toa
+
+  ! Determine pressure range (expect pressure_hl(:,1) ~ surface, (:,nlev+1) ~ TOA)
+  ! Use safe guards if pressure_hl dimension is smaller/wrong => runtime will fail earlier
+  Pmin = minval(pressure_hl(:,nlev+1))
+  Pmax = maxval(pressure_hl(:,1))
+
+  ! compute centre and alpha (guard if Pmax==Pmin)
+  p0 = 0.5_jprb*(Pmax + Pmin)
+  pdiff = (Pmax - Pmin)
+  if (abs(pdiff) < eps) then
+    ! improbable but guard against division by zero
+    alpha = gama / eps
+    write(*,*) 'Warning: Pmax - Pmin trop petit -> alpha forcé'
+  else
+    alpha = gama / pdiff
+  end if
+
+  ! Main loops
+  do jcol = i1, i2
+    do k = 1, nlev
+      ! pression moyenne de la couche k (half-levels k and k+1 must exist)
+      p_mean = 0.5_jprb * (pressure_hl(jcol,k) + pressure_hl(jcol,k+1))
+
+      ! loi hyperbolique bornée entre Lmin et Lmax
+      eff_separation(jcol,k) = 0.5_jprb*(Lmax + Lmin) - 0.5_jprb*(Lmax - Lmin) * tanh(alpha*(p_mean - p0))
+
+      ! borne de sécurité (eviter 0 ou négatif)
+      eff_separation(jcol,k) = max(eps, eff_separation(jcol,k))
+
+      ! fraction safe clamp
+      frac = this%fraction(jcol,k)
+      if (frac < 0.0_jprb) then
+        frac = 0.0_jprb
+      else if (frac > 1.0_jprb) then
+        frac = 1.0_jprb
+      end if
+
+      ! denominators and protection
+      denom1 = eff_separation(jcol,k) * sqrt( max(eps, frac*(1.0_jprb - frac)) )
+      if (denom1 < eps) denom1 = eps
+      this%inv_cloud_effective_size(jcol,k) = 1.0_jprb / denom1
+
+      denom2 = eff_separation(jcol,k) * inhom_sep_factor * sqrt( max(eps, 0.5_jprb*frac*(1.0_jprb - 0.5_jprb*frac)) )
+      if (denom2 < eps) denom2 = eps
+      this%inv_inhom_effective_size(jcol,k) = 1.0_jprb / denom2
+
+    end do
+  end do
+
+  ! diagnostics
+  !CALL writefield_phy('pressure_hl',pressure_hl,nlev+1)
+  !CALL writefield_phy('eff_lx',eff_separation,nlev)
+
+  if (lhook) call dr_hook('radiation_cloud:param_cloud_effective_separation_eta',1,hook_handle)
+
+ end subroutine param_cloud_effective_separation_tanh
+  
   !---------------------------------------------------------------------
   ! Remove "ghost" clouds: those with a cloud fraction that is too