Changeset 4853 for LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_mcica_sw.F90

Timestamp:

Mar 19, 2024, 3:34:21 PM (2 months ago)

Author:

idelkadi

Message:

Ecrad update in LMDZ / Implementation of Ecrad double call in LMDZ

• version 1.6.1 (commit 210d7911380f53a788c3cad73b3cf9b4e022ef87)
• interface routines between lmdz and ecrad grouped in a new "lmdz" directory
• double call method redesigned so as to go through the Ecrad initialization and configuration part only once for the entire simulation
• clean-up in the read.F routine: delete unitules arguments
• modification of default gas model in namelist (default: ECCKD)

File:

• LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_mcica_sw.F90 (modified) (6 diffs)

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

@@ -18 +18 @@
 !   2017-10-23  R. Hogan  Renamed single-character variables
 
+#include "ecrad_config.h"
+
 module radiation_mcica_sw
 
@@ -119 +121 @@
     ! Total cloud cover output from the cloud generator
     real(jprb) :: total_cloud_cover
+
+    ! Temporary storage for more efficient summation
+#ifdef DWD_REDUCTION_OPTIMIZATIONS
+    real(jprb), dimension(nlev+1,3) :: sum_aux
+#else
+    real(jprb) :: sum_up, sum_dn_diff, sum_dn_dir
+#endif
 
     ! Number of g points
@@ -175 +184 @@
         
         ! Sum over g-points to compute and save clear-sky broadband
-        ! fluxes
-        flux%sw_up_clear(jcol,:) = sum(flux_up,1)
+        ! fluxes. Note that the built-in "sum" function is very slow,
+        ! and before being replaced by the alternatives below
+        ! accounted for around 40% of the total cost of this routine.
+#ifdef DWD_REDUCTION_OPTIMIZATIONS
+        ! Optimized summation for the NEC architecture
+        sum_aux(:,:) = 0.0_jprb
+        do jg = 1,ng
+          do jlev = 1,nlev+1
+            sum_aux(jlev,1) = sum_aux(jlev,1) + flux_up(jg,jlev)
+            sum_aux(jlev,2) = sum_aux(jlev,2) + flux_dn_direct(jg,jlev)
+            sum_aux(jlev,3) = sum_aux(jlev,3) + flux_dn_diffuse(jg,jlev)
+          end do
+        end do
+        flux%sw_up_clear(jcol,:) = sum_aux(:,1)
+        flux%sw_dn_clear(jcol,:) = sum_aux(:,2) + sum_aux(:,3)
         if (allocated(flux%sw_dn_direct_clear)) then
-          flux%sw_dn_direct_clear(jcol,:) &
-               &  = sum(flux_dn_direct,1)
-          flux%sw_dn_clear(jcol,:) = sum(flux_dn_diffuse,1) &
-               &  + flux%sw_dn_direct_clear(jcol,:)
-        else
-          flux%sw_dn_clear(jcol,:) = sum(flux_dn_diffuse,1) &
-               &  + sum(flux_dn_direct,1)
+          flux%sw_dn_direct_clear(jcol,:) = sum_aux(:,2)
         end if
+#else
+        ! Optimized summation for the x86-64 architecture
+        do jlev = 1,nlev+1
+          sum_up      = 0.0_jprb
+          sum_dn_diff = 0.0_jprb
+          sum_dn_dir  = 0.0_jprb
+          !$omp simd reduction(+:sum_up, sum_dn_diff, sum_dn_dir)
+          do jg = 1,ng
+            sum_up      = sum_up      + flux_up(jg,jlev)
+            sum_dn_diff = sum_dn_diff + flux_dn_diffuse(jg,jlev)
+            sum_dn_dir  = sum_dn_dir  + flux_dn_direct(jg,jlev)
+          end do
+          flux%sw_up_clear(jcol,jlev) = sum_up
+          flux%sw_dn_clear(jcol,jlev) = sum_dn_diff + sum_dn_dir
+          if (allocated(flux%sw_dn_direct_clear)) then
+            flux%sw_dn_direct_clear(jcol,jlev) = sum_dn_dir
+          end if
+        end do
+#endif
+        
         ! Store spectral downwelling fluxes at surface
-        flux%sw_dn_diffuse_surf_clear_g(:,jcol) = flux_dn_diffuse(:,nlev+1)
-        flux%sw_dn_direct_surf_clear_g(:,jcol)  = flux_dn_direct(:,nlev+1)
+        do jg = 1,ng
+          flux%sw_dn_diffuse_surf_clear_g(jg,jcol) = flux_dn_diffuse(jg,nlev+1)
+          flux%sw_dn_direct_surf_clear_g(jg,jcol)  = flux_dn_direct(jg,nlev+1)
+        end do
 
         ! Do cloudy-sky calculation
@@ -249 +287 @@
             else
               ! Clear-sky layer: copy over clear-sky values
-              reflectance(:,jlev) = ref_clear(:,jlev)
-              transmittance(:,jlev) = trans_clear(:,jlev)
-              ref_dir(:,jlev) = ref_dir_clear(:,jlev)
-              trans_dir_diff(:,jlev) = trans_dir_diff_clear(:,jlev)
-              trans_dir_dir(:,jlev) = trans_dir_dir_clear(:,jlev)
+              do jg = 1,ng
+                reflectance(jg,jlev) = ref_clear(jg,jlev)
+                transmittance(jg,jlev) = trans_clear(jg,jlev)
+                ref_dir(jg,jlev) = ref_dir_clear(jg,jlev)
+                trans_dir_diff(jg,jlev) = trans_dir_diff_clear(jg,jlev)
+                trans_dir_dir(jg,jlev) = trans_dir_dir_clear(jg,jlev)
+              end do
             end if
           end do
@@ -264 +304 @@
           
           ! Store overcast broadband fluxes
-          flux%sw_up(jcol,:) = sum(flux_up,1)
+#ifdef DWD_REDUCTION_OPTIMIZATIONS
+          sum_aux(:,:) = 0.0_jprb
+          do jg = 1,ng
+            do jlev = 1,nlev+1
+              sum_aux(jlev,1) = sum_aux(jlev,1) + flux_up(jg,jlev)
+              sum_aux(jlev,2) = sum_aux(jlev,2) + flux_dn_direct(jg,jlev)
+              sum_aux(jlev,3) = sum_aux(jlev,3) + flux_dn_diffuse(jg,jlev)
+            end do
+          end do
+          flux%sw_up(jcol,:) = sum_aux(:,1)
+          flux%sw_dn(jcol,:) = sum_aux(:,2) + sum_aux(:,3)
           if (allocated(flux%sw_dn_direct)) then
-            flux%sw_dn_direct(jcol,:) = sum(flux_dn_direct,1)
-            flux%sw_dn(jcol,:) = sum(flux_dn_diffuse,1) &
-                 &  + flux%sw_dn_direct(jcol,:)
-          else
-            flux%sw_dn(jcol,:) = sum(flux_dn_diffuse,1) &
-                 &  + sum(flux_dn_direct,1)
+            flux%sw_dn_direct(jcol,:) = sum_aux(:,2)
           end if
-
+#else
+          do jlev = 1,nlev+1
+            sum_up      = 0.0_jprb
+            sum_dn_diff = 0.0_jprb
+            sum_dn_dir  = 0.0_jprb
+            !$omp simd reduction(+:sum_up, sum_dn_diff, sum_dn_dir)
+            do jg = 1,ng
+              sum_up      = sum_up      + flux_up(jg,jlev)
+              sum_dn_diff = sum_dn_diff + flux_dn_diffuse(jg,jlev)
+              sum_dn_dir  = sum_dn_dir  + flux_dn_direct(jg,jlev)
+            end do
+            flux%sw_up(jcol,jlev) = sum_up
+            flux%sw_dn(jcol,jlev) = sum_dn_diff + sum_dn_dir
+            if (allocated(flux%sw_dn_direct)) then
+              flux%sw_dn_direct(jcol,jlev) = sum_dn_dir
+            end if
+          end do
+#endif
+          
           ! Cloudy flux profiles currently assume completely overcast
           ! skies; perform weighted average with clear-sky profile
-          flux%sw_up(jcol,:) =  total_cloud_cover *flux%sw_up(jcol,:) &
-               &  + (1.0_jprb - total_cloud_cover)*flux%sw_up_clear(jcol,:)
-          flux%sw_dn(jcol,:) =  total_cloud_cover *flux%sw_dn(jcol,:) &
-               &  + (1.0_jprb - total_cloud_cover)*flux%sw_dn_clear(jcol,:)
-          if (allocated(flux%sw_dn_direct)) then
-            flux%sw_dn_direct(jcol,:) = total_cloud_cover *flux%sw_dn_direct(jcol,:) &
-                 &  + (1.0_jprb - total_cloud_cover)*flux%sw_dn_direct_clear(jcol,:)
-          end if
+          do jlev = 1, nlev+1
+            flux%sw_up(jcol,jlev) =  total_cloud_cover *flux%sw_up(jcol,jlev) &
+                 &     + (1.0_jprb - total_cloud_cover)*flux%sw_up_clear(jcol,jlev)
+            flux%sw_dn(jcol,jlev) =  total_cloud_cover *flux%sw_dn(jcol,jlev) &
+                 &     + (1.0_jprb - total_cloud_cover)*flux%sw_dn_clear(jcol,jlev)
+            if (allocated(flux%sw_dn_direct)) then
+              flux%sw_dn_direct(jcol,jlev) = total_cloud_cover *flux%sw_dn_direct(jcol,jlev) &
+                   &  + (1.0_jprb - total_cloud_cover)*flux%sw_dn_direct_clear(jcol,jlev)
+            end if
+          end do
           ! Likewise for surface spectral fluxes
-          flux%sw_dn_diffuse_surf_g(:,jcol) = flux_dn_diffuse(:,nlev+1)
-          flux%sw_dn_direct_surf_g(:,jcol)  = flux_dn_direct(:,nlev+1)
-          flux%sw_dn_diffuse_surf_g(:,jcol) = total_cloud_cover *flux%sw_dn_diffuse_surf_g(:,jcol) &
-               &     + (1.0_jprb - total_cloud_cover)*flux%sw_dn_diffuse_surf_clear_g(:,jcol)
-          flux%sw_dn_direct_surf_g(:,jcol) = total_cloud_cover *flux%sw_dn_direct_surf_g(:,jcol) &
-               &     + (1.0_jprb - total_cloud_cover)*flux%sw_dn_direct_surf_clear_g(:,jcol)
-          
+          do jg = 1,ng
+            flux%sw_dn_diffuse_surf_g(jg,jcol) = flux_dn_diffuse(jg,nlev+1)
+            flux%sw_dn_direct_surf_g(jg,jcol)  = flux_dn_direct(jg,nlev+1)
+            flux%sw_dn_diffuse_surf_g(jg,jcol) = total_cloud_cover *flux%sw_dn_diffuse_surf_g(jg,jcol) &
+                 &                 + (1.0_jprb - total_cloud_cover)*flux%sw_dn_diffuse_surf_clear_g(jg,jcol)
+            flux%sw_dn_direct_surf_g(jg,jcol)  = total_cloud_cover *flux%sw_dn_direct_surf_g(jg,jcol) &
+                 &                 + (1.0_jprb - total_cloud_cover)*flux%sw_dn_direct_surf_clear_g(jg,jcol)
+          end do
+
         else
           ! No cloud in profile and clear-sky fluxes already
           ! calculated: copy them over
-          flux%sw_up(jcol,:) = flux%sw_up_clear(jcol,:)
-          flux%sw_dn(jcol,:) = flux%sw_dn_clear(jcol,:)
-          if (allocated(flux%sw_dn_direct)) then
-            flux%sw_dn_direct(jcol,:) = flux%sw_dn_direct_clear(jcol,:)
-          end if
-          flux%sw_dn_diffuse_surf_g(:,jcol) = flux%sw_dn_diffuse_surf_clear_g(:,jcol)
-          flux%sw_dn_direct_surf_g(:,jcol)  = flux%sw_dn_direct_surf_clear_g(:,jcol)
+          do jlev = 1, nlev+1
+            flux%sw_up(jcol,jlev) = flux%sw_up_clear(jcol,jlev)
+            flux%sw_dn(jcol,jlev) = flux%sw_dn_clear(jcol,jlev)
+            if (allocated(flux%sw_dn_direct)) then
+              flux%sw_dn_direct(jcol,jlev) = flux%sw_dn_direct_clear(jcol,jlev)
+            end if
+          end do
+          do jg = 1,ng
+            flux%sw_dn_diffuse_surf_g(jg,jcol) = flux%sw_dn_diffuse_surf_clear_g(jg,jcol)
+            flux%sw_dn_direct_surf_g(jg,jcol)  = flux%sw_dn_direct_surf_clear_g(jg,jcol)
+          end do
 
         end if ! Cloud is present in profile
@@ -307 +378 @@
       else
         ! Set fluxes to zero if sun is below the horizon
-        flux%sw_up(jcol,:) = 0.0_jprb
-        flux%sw_dn(jcol,:) = 0.0_jprb
-        if (allocated(flux%sw_dn_direct)) then
-          flux%sw_dn_direct(jcol,:) = 0.0_jprb
-        end if
-        flux%sw_up_clear(jcol,:) = 0.0_jprb
-        flux%sw_dn_clear(jcol,:) = 0.0_jprb
-        if (allocated(flux%sw_dn_direct_clear)) then
-          flux%sw_dn_direct_clear(jcol,:) = 0.0_jprb
-        end if
-        flux%sw_dn_diffuse_surf_g(:,jcol) = 0.0_jprb
-        flux%sw_dn_direct_surf_g(:,jcol)  = 0.0_jprb
-        flux%sw_dn_diffuse_surf_clear_g(:,jcol) = 0.0_jprb
-        flux%sw_dn_direct_surf_clear_g(:,jcol)  = 0.0_jprb
+        do jlev = 1, nlev+1
+          flux%sw_up(jcol,jlev) = 0.0_jprb
+          flux%sw_dn(jcol,jlev) = 0.0_jprb
+          if (allocated(flux%sw_dn_direct)) then
+            flux%sw_dn_direct(jcol,jlev) = 0.0_jprb
+          end if
+          flux%sw_up_clear(jcol,jlev) = 0.0_jprb
+          flux%sw_dn_clear(jcol,jlev) = 0.0_jprb
+          if (allocated(flux%sw_dn_direct_clear)) then
+            flux%sw_dn_direct_clear(jcol,jlev) = 0.0_jprb
+          end if
+        end do
+        do jg = 1,ng
+          flux%sw_dn_diffuse_surf_g(jg,jcol) = 0.0_jprb
+          flux%sw_dn_direct_surf_g(jg,jcol)  = 0.0_jprb
+          flux%sw_dn_diffuse_surf_clear_g(jg,jcol) = 0.0_jprb
+          flux%sw_dn_direct_surf_clear_g(jg,jcol)  = 0.0_jprb
+        end do
       end if ! Sun above horizon