Changeset 4758 for LMDZ6/branches/LMDZ_ECRad/libf/phylmd/ecrad

Timestamp:

Dec 1, 2023, 10:09:29 PM (14 months ago)

Author:

idelkadi

Message:

• Update of Ecrad in LMDZ version 1.6.1 ​https://github.com/lguez/ecrad.git
• Continuation of work on the implementation of the double call of Ecrad in LMDZ

Location:

LMDZ6/branches/LMDZ_ECRad/libf/phylmd/ecrad

Files:

• CHANGELOG (modified) (1 diff)
• radiation/radiation_aerosol_optics.F90 (modified) (4 diffs)
• radiation/radiation_aerosol_optics_data.F90 (modified) (1 diff)
• radiation/radiation_cloud.F90 (modified) (2 diffs)
• radiation/radiation_cloud_generator.F90 (modified) (1 diff)
• radiation/radiation_ecckd_interface.F90 (modified) (2 diffs)
• radiation/radiation_general_cloud_optics.F90 (modified) (1 diff)
• radiation/radiation_interface.F90 (modified) (4 diffs)
• radiation/radiation_spartacus_sw.F90 (modified) (3 diffs)
• radiation/radiation_two_stream.F90 (modified) (7 diffs)

LMDZ6/branches/LMDZ_ECRad/libf/phylmd/ecrad/CHANGELOG

@@ -13 +13 @@
           saves look-up table averaged to the bands of the radiation scheme
           (general cloud optics only)
+        - Increased security value in single-precision SW
+          reflectance-transmittance calculation from 1e-12 to 1e-6
 
 version 1.6.0 (27 April 2023)

LMDZ6/branches/LMDZ_ECRad/libf/phylmd/ecrad/radiation/radiation_aerosol_optics.F90

@@ -37 +37 @@
     use yomhook,                       only : lhook, dr_hook, jphook
     use radiation_config,              only : config_type
-    use radiation_aerosol_optics_data, only : aerosol_optics_type, &
-         data_setup_aerosol_optics => setup_aerosol_optics
+    use radiation_aerosol_optics_data, only : aerosol_optics_type
     use radiation_io,                  only : nulerr, radiation_abort
 
@@ -58 +57 @@
         ! Read file containing optical properties already in the bands
         ! of the gas-optics scheme
-         if (.not. associated(config%aerosol_optics%setup)) &
-              config%aerosol_optics%setup => data_setup_aerosol_optics
         call config%aerosol_optics%setup(trim(config%aerosol_optics_file_name), &
              &                           iverbose=config%iverbosesetup)
@@ -346 +343 @@
     use easy_netcdf,                   only : netcdf_file
     use radiation_config,              only : config_type
-    use radiation_aerosol_optics_data, only : aerosol_optics_type, &
-         data_setup_aerosol_optics => setup_aerosol_optics
+    use radiation_aerosol_optics_data, only : aerosol_optics_type
     use radiation_spectral_definition, only : SolarReferenceTemperature, &
          &                                    TerrestrialReferenceTemperature
@@ -374 +370 @@
     if (lhook) call dr_hook('radiation_aerosol_optics:setup_general_aerosol_optics_legacy',0,hook_handle)
     ao => config%aerosol_optics
-    ao_legacy%setup => data_setup_aerosol_optics
 
     ! Load file into a local structure

LMDZ6/branches/LMDZ_ECRad/libf/phylmd/ecrad/radiation/radiation_aerosol_optics_data.F90

@@ -135 +135 @@
      logical :: use_monochromatic = .false.
 
-     procedure(setup_aerosol_optics), pointer:: setup => null()
-
    contains
+     procedure :: setup => setup_aerosol_optics
      procedure :: save  => save_aerosol_optics
      procedure :: allocate

LMDZ6/branches/LMDZ_ECRad/libf/phylmd/ecrad/radiation/radiation_cloud.F90

@@ -64 +64 @@
     ! gridbox area for use in representing 3D effects. This variable
     ! is dimensioned (ncol,nlev).
-    real(jprb), allocatable, dimension(:,:) :: inv_cloud_effective_size ! m-1
+    real(jprb), allocatable,  dimension(:,:) :: inv_cloud_effective_size ! m-1
 
     ! Similarly for the in-cloud heterogeneities, used to compute the
@@ -606 +606 @@
 
     use yomhook,                  only : lhook, dr_hook, jphook
+!    USE mod_phys_lmdz_para
 
     class(cloud_type), intent(inout) :: this

LMDZ6/branches/LMDZ_ECRad/libf/phylmd/ecrad/radiation/radiation_cloud_generator.F90

@@ -540 +540 @@
     use radiation_pdf_sampler, only : pdf_sampler_type
     implicit none
-#if defined(__GFORTRAN__) || defined(__PGI) || defined(__NEC__)
-#else
-    !$omp declare simd(sample_from_pdf_simd) uniform(this) &
-    !$omp linear(ref(fsd)) linear(ref(cdf))
-#endif
+!#if defined(__GFORTRAN__) || defined(__PGI) || defined(__NEC__)
+!#else
+!    !$omp declare simd(sample_from_pdf_simd) uniform(this) &
+!    !$omp linear(ref(fsd)) linear(ref(cdf))
+!#endif
     type(pdf_sampler_type), intent(in)  :: this
 

LMDZ6/branches/LMDZ_ECRad/libf/phylmd/ecrad/radiation/radiation_ecckd_interface.F90

@@ -57 +57 @@
       end if
 
-      if (allocated(config%i_band_from_g_sw)) deallocate(config%i_band_from_g_sw)
-       allocate(config%i_band_from_g_sw          (config%n_g_sw))
-      if (allocated(config%i_band_from_reordered_g_sw)) deallocate(config%i_band_from_reordered_g_sw)
-       allocate(config%i_band_from_reordered_g_sw(config%n_g_sw))
-      if (allocated(config%i_g_from_reordered_g_sw)) deallocate(config%i_g_from_reordered_g_sw) 
-       allocate(config%i_g_from_reordered_g_sw   (config%n_g_sw))
+!      if (allocated(config%i_band_from_g_sw)) deallocate(config%i_band_from_g_sw)
+!      allocate(config%i_band_from_g_sw(config%n_g_sw))
+!      if (allocated(config%i_band_from_reordered_g_sw)) deallocate(config%i_band_from_reordered_g_sw)
+!      allocate(config%i_band_from_reordered_g_sw(config%n_g_sw))
+!      if (allocated(config%i_g_from_reordered_g_sw)) deallocate(config%i_g_from_reordered_g_sw)
+!      allocate(config%i_g_from_reordered_g_sw(config%n_g_sw))
+      if (.not.allocated(config%i_band_from_g_sw)) &
+              allocate(config%i_band_from_g_sw(config%n_g_sw))
+      if (.not.allocated(config%i_band_from_reordered_g_sw)) &
+              allocate(config%i_band_from_reordered_g_sw(config%n_g_sw))
+      if (.not.allocated(config%i_g_from_reordered_g_sw)) &
+              allocate(config%i_g_from_reordered_g_sw(config%n_g_sw))
         
       if (config%do_cloud_aerosol_per_sw_g_point) then
@@ -105 +111 @@
       end if
 
-      if (allocated(config%i_band_from_g_lw)) deallocate(config%i_band_from_g_lw)
-       allocate(config%i_band_from_g_lw          (config%n_g_lw))
-      if (allocated(config%i_band_from_reordered_g_lw)) deallocate(config%i_band_from_reordered_g_lw) 
-       allocate(config%i_band_from_reordered_g_lw(config%n_g_lw))
-      if (allocated(config%i_g_from_reordered_g_lw)) deallocate(config%i_g_from_reordered_g_lw) 
-       allocate(config%i_g_from_reordered_g_lw   (config%n_g_lw))
+!      if (allocated(config%i_band_from_g_lw)) deallocate(config%i_band_from_g_lw)
+!      allocate(config%i_band_from_g_lw          (config%n_g_lw))
+!      if (allocated(config%i_band_from_reordered_g_lw)) deallocate(config%i_band_from_reordered_g_lw)
+!      allocate(config%i_band_from_reordered_g_lw(config%n_g_lw))
+!      if (allocated(config%i_g_from_reordered_g_lw)) deallocate(config%i_g_from_reordered_g_lw)
+!      allocate(config%i_g_from_reordered_g_lw   (config%n_g_lw))
+      if (.not.allocated(config%i_band_from_g_lw)) &
+         allocate(config%i_band_from_g_lw          (config%n_g_lw))
+      if (.not.allocated(config%i_band_from_reordered_g_lw)) &
+         allocate(config%i_band_from_reordered_g_lw(config%n_g_lw))
+      if (.not.allocated(config%i_g_from_reordered_g_lw)) &
+         allocate(config%i_g_from_reordered_g_lw   (config%n_g_lw))
+
 
       if (config%do_cloud_aerosol_per_lw_g_point) then

LMDZ6/branches/LMDZ_ECRad/libf/phylmd/ecrad/radiation/radiation_general_cloud_optics.F90

@@ -75 +75 @@
     ! Allocate structures
     if (config%do_sw) then
-     if (allocated(config%cloud_optics_sw)) deallocate(config%cloud_optics_sw)       
+      if (allocated(config%cloud_optics_sw)) deallocate(config%cloud_optics_sw)      
       allocate(config%cloud_optics_sw(config%n_cloud_types))
     end if
 
     if (config%do_lw) then
-     if (allocated(config%cloud_optics_lw)) deallocate(config%cloud_optics_lw)       
+      if (allocated(config%cloud_optics_lw)) deallocate(config%cloud_optics_lw)      
       allocate(config%cloud_optics_lw(config%n_cloud_types))
     end if

LMDZ6/branches/LMDZ_ECRad/libf/phylmd/ecrad/radiation/radiation_interface.F90

@@ -227 +227 @@
     use radiation_general_cloud_optics, only : general_cloud_optics
     use radiation_aerosol_optics, only : add_aerosol_optics
+    USE mod_phys_lmdz_para
+
 
     ! Inputs
@@ -236 +238 @@
     type(thermodynamics_type),intent(in)   :: thermodynamics
     type(gas_type),           intent(in)   :: gas
-    type(cloud_type),         intent(inout):: cloud
+    type(cloud_type),        intent(inout):: cloud
     type(aerosol_type),       intent(in)   :: aerosol
     ! Output
@@ -295 +297 @@
 
     real(jphook) :: hook_handle
+    integer :: jcol, jlev
 
     if (lhook) call dr_hook('radiation_interface:radiation',0,hook_handle)
+
+    if (config%i_solver_sw == ISolverSPARTACUS) then
+         print*,'Dans radiation, mpi_rank, omp_rank, size, chape inv_cloud = ',&
+              mpi_rank, omp_rank, &
+              shape(cloud%inv_cloud_effective_size), &
+              size(cloud%inv_cloud_effective_size)   
+!         do jcol=istartcol, iendcol
+!              do jlev=1,nlev
+!              print*,'Entree radiation_interf, mpi_rank, omp_rank, jcol, jlev &
+!           &  cloud%inv_cloud_effective_size =',mpi_rank, omp_rank, jcol, jlev, &
+!           &  cloud%inv_cloud_effective_size(jcol,jlev)
+!              enddo
+!          enddo
+    endif
+!    cloud%inv_cloud_effective_size=0.05_jprb
 
     if (thermodynamics%pressure_hl(istartcol,2) &
@@ -456 +474 @@
         else if (config%i_solver_sw == ISolverSPARTACUS) then
           ! Compute fluxes using the SPARTACUS shortwave solver
+!             cloud%inv_cloud_effective_size=0.05_jprb
+!             do jcol=istartcol, iendcol
+!              do jlev=1,nlev
+!              print*,'jcol, jlev, dans radiation_interf i &
+!               &       cloud%inv_cloud_effective_size =',jcol, jlev, &
+!               cloud%inv_cloud_effective_size(jcol,jlev)
+!              enddo
+!             enddo
           call solver_spartacus_sw(nlev,istartcol,iendcol, &
                &  config, single_level, thermodynamics, cloud, & 

LMDZ6/branches/LMDZ_ECRad/libf/phylmd/ecrad/radiation/radiation_spartacus_sw.F90

@@ -87 +87 @@
     use radiation_constants, only      : Pi, GasConstantDryAir, &
          &                               AccelDueToGravity
+    USE mod_phys_lmdz_para
 
     implicit none
@@ -326 +327 @@
       write(nulout,'(a)',advance='no') '  Processing columns'
     end if
+
+      print*,'Dans radiation_spartacus, mpi_rank, omp_rank, &
+              size, chape inv_cloud = ',&
+              mpi_rank, omp_rank, &
+              shape(cloud%inv_cloud_effective_size), &
+              size(cloud%inv_cloud_effective_size)
 
     ! Main loop over columns
@@ -497 +504 @@
                &  .not. (nreg == 2 .and. cloud%fraction(jcol,jlev) &
                &  > 1.0-config%cloud_fraction_threshold)) then
+!           print*,' Dans radiation_spartacus mpi_rank, omp_rank, jcol, jlev, &
+!            &     cloud%inv_cloud_effective_size =', mpi_rank, &
+!            &     omp_rank, jcol, jlev, &
+!            &     cloud%inv_cloud_effective_size(jcol,jlev)
             if (cloud%inv_cloud_effective_size(jcol,jlev) > 0.0_jprb) then
               ! 3D effects are only simulated if

LMDZ6/branches/LMDZ_ECRad/libf/phylmd/ecrad/radiation/radiation_two_stream.F90

@@ -20 +20 @@
 !   2021-02-19  R Hogan  Security for shortwave singularity
 !   2022-11-22  P Ukkonen/R Hogan  Single precision uses no double precision
+!   2023-09-28  R Hogan  Increased security for single-precision SW "k"
 
 module radiation_two_stream
@@ -212 +213 @@
       if (od(jg) > 1.0e-3_jprd) then
         k_exponent = sqrt(max((gamma1(jg) - gamma2(jg)) * (gamma1(jg) + gamma2(jg)), &
-             1.E-12_jprd)) ! Eq 18 of Meador & Weaver (1980)
+             1.0e-12_jprd)) ! Eq 18 of Meador & Weaver (1980)
         exponential = exp_fast(-k_exponent*od(jg))
         exponential2 = exponential*exponential
@@ -235 +236 @@
       else
         k_exponent = sqrt(max((gamma1(jg) - gamma2(jg)) * (gamma1(jg) + gamma2(jg)), &
-             1.E-12_jprd)) ! Eq 18 of Meador & Weaver (1980)
+             1.0e-12_jprd)) ! Eq 18 of Meador & Weaver (1980)
         reflectance(jg) = gamma2(jg) * od(jg)
         transmittance(jg) = (1.0_jprb - k_exponent*od(jg)) / (1.0_jprb + od(jg)*(gamma1(jg)-k_exponent))
@@ -312 +313 @@
       gamma2 = factor * (1.0_jprb - asymmetry(jg))
       k_exponent = sqrt(max((gamma1 - gamma2) * (gamma1 + gamma2), &
-           1.E-12_jprb)) ! Eq 18 of Meador & Weaver (1980)
+           1.0e-12_jprb)) ! Eq 18 of Meador & Weaver (1980)
       if (od(jg) > 1.0e-3_jprb) then
         exponential = exp_fast(-k_exponent*od(jg))
@@ -646 +647 @@
       alpha2(jg) = gamma1(jg)*gamma3(jg) + gamma2(jg)*gamma4(jg) ! Eq. 17
       ! The following line crashes inexplicably with gfortran 8.5.0 in
-      ! single precision - try a later version
+      ! single precision - try a later version. Note that the minimum
+      ! value is needed to produce correct results for single
+      ! scattering albedos very close to or equal to one.
+#ifdef PARKIND1_SINGLE
+      k_exponent(jg) = sqrt(max((gamma1(jg) - gamma2(jg)) * (gamma1(jg) + gamma2(jg)), &
+           &       1.0e-6_jprb)) ! Eq 18
+#else
       k_exponent(jg) = sqrt(max((gamma1(jg) - gamma2(jg)) * (gamma1(jg) + gamma2(jg)), &
            &       1.0e-12_jprb)) ! Eq 18
+#endif
     end do
 
@@ -665 +673 @@
       ! Meador & Weaver (1980) Eq. 25
       ref_diff(jg) = gamma2(jg) * (1.0_jprb - exponential2) * reftrans_factor
-        
-      ! Meador & Weaver (1980) Eq. 26
-      trans_diff(jg) = k_2_exponential * reftrans_factor
-        
+      !ref_diff(jg)       = max(0.0_jprb, min(ref_diff(jg)), 1.0_jprb)
+
+      ! Meador & Weaver (1980) Eq. 26, with security (which is
+      ! sometimes needed, but apparently not on ref_diff)
+      trans_diff(jg) = max(0.0_jprb, min(k_2_exponential * reftrans_factor, 1.0_jprb-ref_diff(jg)))
+
       ! Here we need mu0 even though it wasn't in Meador and Weaver
       ! because we are assuming the incoming direct flux is defined to
@@ -694 +704 @@
       ref_dir(jg)        = max(0.0_jprb, min(ref_dir(jg), mu0*(1.0_jprb-trans_dir_dir(jg))))
       trans_dir_diff(jg) = max(0.0_jprb, min(trans_dir_diff(jg), mu0*(1.0_jprb-trans_dir_dir(jg))-ref_dir(jg)))
-
     end do