Use correct input file for general aerosol optics in IFS-style drivers

.github/workflows/build.yml

@@ -81,7 +81,7 @@ jobs:
           brew install ninja
         else
           sudo apt-get update
-          sudo apt-get install ninja-build libcurl4-openssl-dev
+          sudo apt-get install ninja-build libcurl4-openssl-dev cdo
         fi
 
         printenv
@@ -173,7 +173,7 @@ jobs:
         dependency_branch: develop
         dependency_cmake_options: |
           ecmwf-ifs/fiat: "-G Ninja -DCMAKE_BUILD_TYPE=${{ matrix.build_type }} -DENABLE_TESTS=OFF"
-        cmake_options:    "-G Ninja -DCMAKE_BUILD_TYPE=${{ matrix.build_type }} ${{ matrix.cmake_options }} -DENABLE_SINGLE_PRECISION=ON"
+        cmake_options:    "-G Ninja -DCMAKE_BUILD_TYPE=${{ matrix.build_type }} ${{ matrix.cmake_options }} -DENABLE_SINGLE_PRECISION=ON -DENABLE_BITIDENTITY_TESTING=ON"
         ctest_options: "${{ matrix.ctest_options }}"
 
 #    - name: Codecov Upload

CMakeLists.txt

@@ -93,6 +93,23 @@ if( HAVE_PRINT_ENTRAPMENT_DATA )
     list( APPEND ECRAD_COMPILE_DEFINITIONS PRINT_ENTRAPMENT_DATA )
 endif()
 
+ecbuild_add_option(
+    FEATURE BITIDENTITY_TESTING
+    DEFAULT OFF
+    DESCRIPTION "Prescribe effective radii, seed and cloud_fraction in blocked drivers for bit-identical results"
+)
+if( HAVE_BITIDENTITY_TESTING )
+    list( APPEND ECRAD_COMPILE_DEFINITIONS BITIDENTITY_TESTING )
+endif()
+
+find_program( CDO NAMES cdo )
+if( "${CDO}" STREQUAL "CDO-NOTFOUND" )
+    ecbuild_info( "Could NOT find cdo" )
+else()
+    ecbuild_info( "Found program cdo [${CDO}]" )
+    set( HAVE_CDO ON CACHE BOOL "cdo binary" )
+endif()
+
 include(ecrad_compile_flags)
 
 # Dummy drhook module (precision-independent)

Makefile

@@ -61,6 +61,12 @@ $(info *** Building with NetCDF4/HDF5 support)
 CPPFLAGS += -DNC_NETCDF4
 endif
 
+# Optionally overwrite effective radii, seed, and cloud fraction
+# in IFS variants to retain bit-identical results
+ifdef BITIDENTITY_TESTING
+CPPFLAGS += -DBITIDENTITY_TESTING
+endif
+
 # Consolidate flags
 export FC
 export FCFLAGS = $(WARNFLAGS) $(BASICFLAGS) $(CPPFLAGS) -I../include \

driver/ecrad_ifs_driver.F90

@@ -421,7 +421,15 @@ program ecrad_ifs_driver
              &  flux%sw_dn_direct(:,nlev+1), flux%sw_dn_direct_clear(:,nlev+1), flux_sw_direct_normal, &
              &  flux_uv, flux_par, &
              &  flux_par_clear, flux%sw_dn(:,1), emissivity_out, flux%lw_derivatives, flux_diffuse_band, &
-             &  flux_direct_band)
+             &  flux_direct_band &
+#ifdef BITIDENTITY_TESTING
+            ! To validate results against standalone ecrad, we overwrite effective
+            ! radii, cloud overlap and seed with input values
+             &  ,pre_liq=cloud%re_liq, pre_ice=cloud%re_ice, &
+             &  pcloud_overlap=cloud%overlap_param, &
+             &  iseed=single_level%iseed &
+#endif
+             & )
       end do
       !$OMP END PARALLEL DO
 

driver/ecrad_ifs_driver_blocked.F90

@@ -103,7 +103,9 @@ program ecrad_ifs_driver
   ! Bespoke data types to set-up the blocked memory layout
   type(ifs_config_type)        :: ifs_config
   real(kind=jprb), allocatable :: zrgp(:,:,:) ! monolithic IFS data structure
+#ifdef BITIDENTITY_TESTING
   integer, allocatable         :: iseed(:,:) ! Seed for random number generator
+#endif
 
   integer :: ncol, nlev         ! Number of columns and levels
   integer :: nproma             ! block size
@@ -365,7 +367,11 @@ program ecrad_ifs_driver
   call ifs_copy_inputs_to_blocked(driver_config, ifs_config, yradiation,&
         & ncol, nlev, single_level, thermodynamics, gas, cloud, aerosol,&
         & sin_latitude, longitude_rad, land_frac, pressure_fl, temperature_fl,&
-        & zrgp, iseed=iseed)
+        & zrgp &
+#ifdef BITIDENTITY_TESTING
+        &, iseed=iseed &
+#endif
+        & )
 
   ! --------------------------------------------------------
   ! Section 4b: Call radiation_scheme with blocked memory data
@@ -437,6 +443,14 @@ program ecrad_ifs_driver
              &  zrgp(1,ifs_config%iparcf,ib),zrgp(1,ifs_config%itincf,ib), &
              &  zrgp(1,ifs_config%iemit,ib) ,zrgp(1,ifs_config%ilwderivative,ib), &
              &  zrgp(1,ifs_config%iswdiffuseband,ib), zrgp(1,ifs_config%iswdirectband,ib)&
+#ifdef BITIDENTITY_TESTING
+            ! To validate results against standalone ecrad, we overwrite effective
+            ! radii, cloud overlap and seed with input values
+             &  ,pre_liq=zrgp(1,ifs_config%ire_liq,ib), &
+             &  pre_ice=zrgp(1,ifs_config%ire_ice,ib), &
+             &  pcloud_overlap=zrgp(1,ifs_config%ioverlap,ib), &
+             &  iseed=iseed(:,ib) &
+#endif
              & )
       end do
       !$OMP END PARALLEL DO

driver/ifs_blocking.F90

@@ -29,7 +29,9 @@ module ifs_blocking
         &     ifrted, ifrsodc, ifrtedc, iemit, isudu, iuvdf, iparf, iparcf, itincf, ifdir, ifdif,     &
         &     ilwderivative, iswdirectband, iswdiffuseband, ifrso, iswfc, ifrth, ilwfc, iaer,         &
         &     iich4, iin2o, ino2, ic11, ic12, igix, iico2, iccno, ic22, icl4
+#ifdef BITIDENTITY_TESTING
     integer :: ire_liq, ire_ice, ioverlap
+#endif
     integer :: ifldstot
   end type ifs_config_type
 
@@ -182,9 +184,13 @@ subroutine ifs_setup_indices (driver_config, ifs_config, yradiation, nlev)
                                 ! end of local variables
 
                                   ! start of standalone inputs workaround variables
+#ifdef BITIDENTITY_TESTING
+  ! To validate results against standalone ecrad, we overwrite effective
+  ! radii, cloud overlap and seed with input values
   ifs_config%ire_liq =indrad(inext,nlev,.true.)
   ifs_config%ire_ice =indrad(inext,nlev,.true.)
   ifs_config%ioverlap =indrad(inext,nlev-1,.true.)
+#endif
                                   ! end of standalone inputs workaround variables
 
   ifldsin = iinend - iinbeg +1
@@ -254,9 +260,11 @@ subroutine ifs_setup_indices (driver_config, ifs_config, yradiation, nlev)
     write(nulout,'("ic12   =",i0)')ifs_config%ic12
     write(nulout,'("ic22   =",i0)')ifs_config%ic22
     write(nulout,'("icl4   =",i0)')ifs_config%icl4
+#ifdef BITIDENTITY_TESTING
     write(nulout,'("ire_liq=",i0)')ifs_config%ire_liq
     write(nulout,'("ire_ice=",i0)')ifs_config%ire_ice
     write(nulout,'("ioverlap=",i0)')ifs_config%ioverlap
+#endif
     write(nulout,'("ifldsin =",i0)')ifldsin
     write(nulout,'("ifldsout=",i0)')ifldsout
     write(nulout,'("ifldstot=",i0)')ifs_config%ifldstot
@@ -390,7 +398,7 @@ subroutine ifs_copy_inputs_to_blocked ( &
           zrgp(1:il,ifs_config%ialp+jalb-1,ib)  =  single_level%sw_albedo(ibeg:iend,jalb)
         end do
       end if
-      
+
       do jlev=1,nlev
         zrgp(1:il,ifs_config%iti+jlev-1,ib)   = temperature_fl(ibeg:iend,jlev) ! full level temperature
         zrgp(1:il,ifs_config%ipr+jlev-1,ib)   = pressure_fl(ibeg:iend,jlev) ! full level pressure
@@ -451,6 +459,9 @@ subroutine ifs_copy_inputs_to_blocked ( &
       ! enddo
 
       ! local workaround variables for standalone input files
+#ifdef BITIDENTITY_TESTING
+      ! To validate results against standalone ecrad, we overwrite effective
+      ! radii, cloud overlap and seed with input values
       if (rad_config%do_clouds) then
         do jlev=1,nlev
           ! missing full-level temperature and pressure as well as land-sea-mask
@@ -474,6 +485,7 @@ subroutine ifs_copy_inputs_to_blocked ( &
         enddo
         if(present(iseed)) iseed(1:il,ib) = 0
       endif ! do_clouds
+#endif
     enddo
     !$OMP END PARALLEL DO
 

ifs/radiation_scheme.F90

@@ -15,7 +15,9 @@ SUBROUTINE RADIATION_SCHEME &
      &  PFLUX_DIR, PFLUX_DIR_CLEAR, PFLUX_DIR_INTO_SUN, &
      &  PFLUX_UV, PFLUX_PAR, PFLUX_PAR_CLEAR, &
      &  PFLUX_SW_DN_TOA, PEMIS_OUT, PLWDERIVATIVE, &
-     &  PSWDIFFUSEBAND, PSWDIRECTBAND)
+     &  PSWDIFFUSEBAND, PSWDIRECTBAND, &
+     ! OPTIONAL ARGUMENTS for bit-identical results in tests
+     &  PRE_LIQ, PRE_ICE, ISEED, PCLOUD_OVERLAP)
 
 ! RADIATION_SCHEME - Interface to modular radiation scheme
 !
@@ -189,6 +191,12 @@ SUBROUTINE RADIATION_SCHEME &
 REAL(KIND=JPRB),  INTENT(OUT) :: PSWDIFFUSEBAND(KLON,YRADIATION%YRERAD%NSW)
 REAL(KIND=JPRB),  INTENT(OUT) :: PSWDIRECTBAND (KLON,YRADIATION%YRERAD%NSW)
 
+! Optional input arguments (Added for validating against ecrad standalone!)
+REAL(KIND=JPRB), INTENT(IN), OPTIONAL :: PRE_LIQ(KLON, KLEV)
+REAL(KIND=JPRB), INTENT(IN), OPTIONAL :: PRE_ICE(KLON, KLEV)
+INTEGER,         INTENT(IN), OPTIONAL :: ISEED(KLON)
+REAL(KIND=JPRB), INTENT(IN), OPTIONAL :: PCLOUD_OVERLAP(KLON, KLEV-1)
+
 ! LOCAL VARIABLES
 TYPE(SINGLE_LEVEL_TYPE)   :: SINGLE_LEVEL
 TYPE(THERMODYNAMICS_TYPE) :: THERMODYNAMICS
@@ -352,6 +360,12 @@ SUBROUTINE RADIATION_SCHEME &
 ! Simple initialization of the seeds for the Monte Carlo scheme
 call single_level%init_seed_simple(kidia, kfdia)
 
+! Added for bit-identity validation against ecrad standalone:
+! Overwrite seed with user-specified values
+if (present(iseed)) then
+   single_level%iseed(kidia:kfdia) = iseed(kidia:kfdia)
+end if
+
 ! Set the solar spectrum scaling, if required
 IF (YRERAD%NSOLARSPECTRUM == 1) THEN
   ALLOCATE(SINGLE_LEVEL%SPECTRAL_SOLAR_SCALING(RAD_CONFIG%N_BANDS_SW))
@@ -367,17 +381,25 @@ SUBROUTINE RADIATION_SCHEME &
 YLCLOUD%FRACTION(KIDIA:KFDIA,:) = PCLOUD_FRAC(KIDIA:KFDIA,:)
 
 ! Compute effective radii and convert to metres
+IF (PRESENT(PRE_LIQ)) THEN
+  YLCLOUD%RE_LIQ(KIDIA:KFDIA,:) = PRE_LIQ(KIDIA:KFDIA,:)
+ELSE
 CALL LIQUID_EFFECTIVE_RADIUS(YRERAD, &
      &  KIDIA, KFDIA, KLON, KLEV, &
      &  PPRESSURE, PTEMPERATURE, PCLOUD_FRAC, PQ_LIQUID, PQ_RAIN, &
      &  PLAND_SEA_MASK, PCCN_LAND, PCCN_SEA, &
      &  ZRE_LIQUID_UM) !, PPERT=PPERT)
 YLCLOUD%RE_LIQ(KIDIA:KFDIA,:) = ZRE_LIQUID_UM(KIDIA:KFDIA,:) * 1.0E-6_JPRB
+ENDIF
 
+IF (PRESENT(PRE_ICE)) THEN
+  YLCLOUD%RE_ICE(KIDIA:KFDIA,:) = PRE_ICE(KIDIA:KFDIA,:)
+ELSE
 CALL ICE_EFFECTIVE_RADIUS(YRERAD, KIDIA, KFDIA, KLON, KLEV, &
      &  PPRESSURE, PTEMPERATURE, PCLOUD_FRAC, PQ_ICE, PQ_SNOW, PGEMU, &
      &  ZRE_ICE_UM) !, PPERT=PPERT)
 YLCLOUD%RE_ICE(KIDIA:KFDIA,:) = ZRE_ICE_UM(KIDIA:KFDIA,:) * 1.0E-6_JPRB
+ENDIF
 
 ! Get the cloud overlap decorrelation length (for cloud boundaries),
 ! in km, according to the parameterization specified by NDECOLAT,
@@ -396,6 +418,12 @@ SUBROUTINE RADIATION_SCHEME &
       &                       ISTARTCOL=JLON, IENDCOL=JLON)
 ENDDO
 
+! Added for bit-identity validation against ecrad standalone:
+! Overwrite overlap param with provided value
+if(present(PCLOUD_OVERLAP)) then
+  YLCLOUD%OVERLAP_PARAM(KIDIA:KFDIA,:) = PCLOUD_OVERLAP(KIDIA:KFDIA,:)
+endif
+
 ! Cloud water content fractional standard deviation is configurable
 ! from namelist NAERAD but must be globally constant. Before it was
 ! hard coded at 1.0.

ifs/radiation_setup.F90

@@ -292,7 +292,11 @@ SUBROUTINE SETUP_RADIATION_SCHEME(PRADIATION,LDOUTPUT,FILE_NAME)
 
       ! The default aerosol optics file is the following - please
       ! update here, not in radiation/module/radiation_config.F90
-      RAD_CONFIG%AEROSOL_OPTICS_OVERRIDE_FILE_NAME = 'aerosol_ifs_rrtm_46R1_with_NI_AM.nc'
+      IF (RAD_CONFIG%USE_GENERAL_AEROSOL_OPTICS) THEN
+        RAD_CONFIG%AEROSOL_OPTICS_OVERRIDE_FILE_NAME = 'aerosol_ifs_49R1_20230119.nc'
+      ELSE
+        RAD_CONFIG%AEROSOL_OPTICS_OVERRIDE_FILE_NAME = 'aerosol_ifs_rrtm_46R1_with_NI_AM.nc'
+      END IF
 
     ELSE
       ! Using Tegen climatology
@@ -439,7 +443,7 @@ SUBROUTINE SETUP_RADIATION_SCHEME(PRADIATION,LDOUTPUT,FILE_NAME)
     ! second emissivity to represent values within it.
     CALL YDERAD%YSPECTPLANCK%INIT(2, [ 8.0E-6_JPRB, 13.0E-6_JPRB ], &
          &  [ 1,2,1 ])
-      
+
     ! Populate the mapping between the 14 RRTM shortwave bands and the
     ! 6 albedo inputs.
     YDERAD%NSW = 6
@@ -518,25 +522,34 @@ SUBROUTINE SETUP_RADIATION_SCHEME(PRADIATION,LDOUTPUT,FILE_NAME)
          &  PRADIATION%NWEIGHT_PAR, PRADIATION%IBAND_PAR, PRADIATION%WEIGHT_PAR,&
          &  'photosynthetically active radiation, PAR')
 
-    IF (YDERAD%NAERMACC > 0) THEN
-      ! With the MACC aerosol climatology we need to add in the
-      ! background aerosol afterwards using the Tegen arrays.  In this
-      ! case we first configure the background aerosol mass-extinction
-      ! coefficient at 550 nm, which corresponds to the 10th RRTMG
-      ! shortwave band.
-      PRADIATION%TROP_BG_AER_MASS_EXT  = DRY_AEROSOL_MASS_EXTINCTION(RAD_CONFIG,&
-           &                                   ITYPE_TROP_BG_AER, 550.0E-9_JPRB)
-      PRADIATION%STRAT_BG_AER_MASS_EXT = DRY_AEROSOL_MASS_EXTINCTION(RAD_CONFIG,&
-           &                                   ITYPE_STRAT_BG_AER, 550.0E-9_JPRB)
-
-      WRITE(NULOUT,'(a,i0)') 'Tropospheric background uses aerosol type ',&
-           &                 ITYPE_TROP_BG_AER
-      WRITE(NULOUT,'(a,i0)') 'Stratospheric background uses aerosol type ',&
-           &                 ITYPE_STRAT_BG_AER
-    ELSE
-      PRADIATION%TROP_BG_AER_MASS_EXT  = 0.0_JPRB
-      PRADIATION%STRAT_BG_AER_MASS_EXT = 0.0_JPRB
-    ENDIF
+    ! PRADIATION%TROP_BG_AER_MASS_EXT  = 0.0_JPRB
+    ! PRADIATION%STRAT_BG_AER_MASS_EXT = 0.0_JPRB
+    ! IF (YDERAD%NAERMACC > 0) THEN
+    !   ! With the MACC aerosol climatology we need to add in the
+    !   ! background aerosol afterwards using the Tegen arrays.  In this
+    !   ! case we first configure the background aerosol mass-extinction
+    !   ! coefficient at 550 nm, which corresponds to the 10th RRTMG
+    !   ! shortwave band.
+    !   IF (ITYPE_TROP_BG_AER > 0) THEN
+    !     PRADIATION%TROP_BG_AER_MASS_EXT  = DRY_AEROSOL_MASS_EXTINCTION(RAD_CONFIG,&
+    !         &                                   ITYPE_TROP_BG_AER, 550.0e-9_JPRB)
+    !     WRITE(NULOUT,'(a,i2,a,e12.4,a)') 'Tropospheric background:  aerosol type ',&
+    !         &  ITYPE_TROP_BG_AER, ', 550-nm mass-extinction coefficient ', &
+    !         &  PRADIATION%TROP_BG_AER_MASS_EXT, ' m2 kg-1'
+    !   ELSE
+    !     WRITE(NULOUT,'(a)') 'No tropospheric background aerosol'
+    !   ENDIF
+
+    !   IF (ITYPE_STRAT_BG_AER > 0) THEN
+    !     PRADIATION%STRAT_BG_AER_MASS_EXT = DRY_AEROSOL_MASS_EXTINCTION(RAD_CONFIG,&
+    !         &                                   ITYPE_STRAT_BG_AER, 550.0e-9_JPRB)
+    !     WRITE(NULOUT,'(a,i2,a,e12.4,a)') 'Stratospheric background: aerosol type ',&
+    !         &  ITYPE_STRAT_BG_AER, ', 550-nm mass-extinction coefficient ', &
+    !         &  PRADIATION%STRAT_BG_AER_MASS_EXT, ' m2 kg-1'
+    !   ELSE
+    !     WRITE(NULOUT,'(a)') 'No stratospheric background aerosol'
+    !   ENDIF
+    ! ENDIF
 
     IF (IVERBOSESETUP > 1) THEN
       WRITE(NULOUT,'(a)') '-------------------------------------------------------------------------------'

include/radiation_scheme.intfb.h

@@ -16,7 +16,8 @@ SUBROUTINE RADIATION_SCHEME &
  & PFLUX_DIR, PFLUX_DIR_CLEAR, PFLUX_DIR_INTO_SUN, &
  & PFLUX_UV, PFLUX_PAR, PFLUX_PAR_CLEAR, &
  & PFLUX_SW_DN_TOA, PEMIS_OUT, PLWDERIVATIVE, &
- & PSWDIFFUSEBAND, PSWDIRECTBAND)
+ & PSWDIFFUSEBAND, PSWDIRECTBAND, &
+ & PRE_LIQ, PRE_ICE, ISEED, PCLOUD_OVERLAP)
 use parkind1 , only:&
  & jpim,&
  & jprb
@@ -79,5 +80,9 @@ REAL(KIND=JPRB), INTENT(OUT) :: PEMIS_OUT(KLON)
 REAL(KIND=JPRB), INTENT(OUT) :: PLWDERIVATIVE(KLON,KLEV+1)
 REAL(KIND=JPRB), INTENT(OUT) :: PSWDIFFUSEBAND(KLON,YRADIATION%YRERAD%NSW)
 REAL(KIND=JPRB), INTENT(OUT) :: PSWDIRECTBAND (KLON,YRADIATION%YRERAD%NSW)
+REAL(KIND=JPRB), INTENT(IN), OPTIONAL :: PRE_LIQ(KLON, KLEV)
+REAL(KIND=JPRB), INTENT(IN), OPTIONAL :: PRE_ICE(KLON, KLEV)
+INTEGER,         INTENT(IN), OPTIONAL :: ISEED(KLON)
+REAL(KIND=JPRB), INTENT(IN), OPTIONAL :: PCLOUD_OVERLAP(KLON, KLEV-1)
 END SUBROUTINE RADIATION_SCHEME
 end interface