Implement a sub-band weighted UV-index diagnostic

examples/sp_uv/ref_trop_411rd.uvindex

@@ -0,0 +1,54 @@
+netcdf cdl_trop_411rd.uvindex{                                                  
+
+dimensions:
+    lon            =   1;
+    lat            =   1;
+    plev           =  51;
+
+
+variables:
+    float lon(lon);                                                           
+             lon:units = "degree";                                            
+             lon:title = "longitude";                                        
+    float lat(lat);                                                           
+             lat:units = "degree";                                            
+             lat:title = "latitude";                                         
+    float plev(plev);                                                         
+             plev:units = "pa";                                               
+             plev:title = "pressure";                                        
+
+    float uvindex(plev,lat,lon);                                               
+             uvindex:units = "None";                                          
+             uvindex:title = "UV-index";                                                     
+
+data:
+              lon =  0.000000E+00;
+              lat =  0.000000E+00;
+             plev =  0.182812E-02, 0.292500E-02, 0.503500E-02, 0.973500E-02,
+                     0.209500E-01, 0.488500E-01, 0.120400E+00, 0.306000E+00,
+                     0.770000E+00, 0.185000E+01, 0.420000E+01, 0.895000E+01,
+                     0.180000E+02, 0.347500E+02, 0.655000E+02, 0.100700E+03,
+                     0.137500E+03, 0.189500E+03, 0.262500E+03, 0.365500E+03,
+                     0.513000E+03, 0.726000E+03, 0.103600E+04, 0.149150E+04,
+                     0.216650E+04, 0.278500E+04, 0.325000E+04, 0.379500E+04,
+                     0.444500E+04, 0.522500E+04, 0.615500E+04, 0.727500E+04,
+                     0.863000E+04, 0.102350E+05, 0.121500E+05, 0.144000E+05,
+                     0.169000E+05, 0.197500E+05, 0.230000E+05, 0.266500E+05,
+                     0.307500E+05, 0.353500E+05, 0.405000E+05, 0.462000E+05,
+                     0.525500E+05, 0.596000E+05, 0.674000E+05, 0.760000E+05,
+                     0.854500E+05, 0.958500E+05, 0.101300E+06;
+          uvindex =  0.485189E+03, 0.485189E+03, 0.485189E+03, 0.485189E+03,
+                     0.485188E+03, 0.485183E+03, 0.485167E+03, 0.485124E+03,
+                     0.485020E+03, 0.484859E+03, 0.484663E+03, 0.484017E+03,
+                     0.481408E+03, 0.473449E+03, 0.451196E+03, 0.416864E+03,
+                     0.380697E+03, 0.330980E+03, 0.272796E+03, 0.214722E+03,
+                     0.161315E+03, 0.114071E+03, 0.766958E+02, 0.502044E+02,
+                     0.334822E+02, 0.271063E+02, 0.241858E+02, 0.219101E+02,
+                     0.202427E+02, 0.189017E+02, 0.177712E+02, 0.169043E+02,
+                     0.163464E+02, 0.159793E+02, 0.156828E+02, 0.153671E+02,
+                     0.150544E+02, 0.147246E+02, 0.143829E+02, 0.140304E+02,
+                     0.136636E+02, 0.132721E+02, 0.128537E+02, 0.124030E+02,
+                     0.119144E+02, 0.113869E+02, 0.108193E+02, 0.102052E+02,
+                     0.955016E+01, 0.885217E+01, 0.849889E+01;
+
+}

examples/sp_uv/run_me

@@ -45,6 +45,11 @@ fmove tropical trop_320_sg_ch
 Cl_run_cdf -B tropical -s sp_uv_photol -S -R 1 320 -g 4 -r -C 5
 fmove tropical trop_320
 
+Cl_run_cdf -B tropical -s sp_411rd_photol -S -R 1 411 -g 4 -r -C 5 -sg -ch 411
+fmove tropical trop_411rd_sg_ch
+Cl_run_cdf -B tropical -s sp_411rd_photol -S -R 1 411 -g 4 -r -C 5
+fmove tropical trop_411rd
+
 cp tropical.surflw tropical.surf
 Cl_run_cdf -B tropical -s $RAD_DATA/spectra/ga7/sp_lw_ga7 -I -R 1 9 -g 4 -C 5 -c
 fmove tropical trop_lw_ga7
@@ -72,6 +77,8 @@ Ccdf2cdl -o cdl_trop_xsw_photol_nlte.hrts          trop_xsw_photol_nlte.hrts > /
 Ccdf2cdl -o cdl_trop_sw_ga7.hrts                   trop_sw_ga7.hrts > /dev/null
 Ccdf2cdl -o cdl_trop_320_sg_ch.ph_rate_2           trop_320_sg_ch.ph_rate_2 > /dev/null
 Ccdf2cdl -o cdl_trop_320.hrts                      trop_320.hrts > /dev/null
+Ccdf2cdl -o cdl_trop_411rd_sg_ch.uvindex           trop_411rd_sg_ch.uvindex > /dev/null
+Ccdf2cdl -o cdl_trop_411rd.uvindex                 trop_411rd.uvindex > /dev/null
 Ccdf2cdl -o cdl_trop_lw_ga7.hrts                   trop_lw_ga7.hrts > /dev/null
 Ccdf2cdl -o cdl_trop_lw_ga7_nlte.hrts              trop_lw_ga7_nlte.hrts > /dev/null
 
@@ -108,6 +115,10 @@ diff -q cdl_trop_sw_ga7.hrts              ref_trop_sw_ga7.hrts
   echo 'Matched ifort output')                                                   || \
     ierr=1
 diff -q cdl_trop_320.hrts                 ref_trop_320.hrts                      || ierr=1
+diff -q cdl_trop_411rd_sg_ch.uvindex      gfortran_trop_411rd_sg_ch.uvindex      || \
+diff -q cdl_trop_411rd_sg_ch.uvindex      ref_trop_411rd_sg_ch.uvindex           || \
+    ierr=1
+diff -q cdl_trop_411rd.uvindex            ref_trop_411rd.uvindex                 || ierr=1
 diff -q cdl_trop_lw_ga7.hrts              ref_trop_lw_ga7.hrts                   || ierr=1
 diff -q cdl_trop_lw_ga7_nlte.hrts         ref_trop_lw_ga7_nlte.hrts              || ierr=1
 

sbin/fmove

@@ -3,7 +3,7 @@
 # Script to change the base name of all output files from a run of the
 # radiation code.
 #
-for suffix in uflx dflx sflx nflx vflx hrts radn cfi cff aflx \
+for suffix in uflx dflx sflx nflx vflx hrts radn cfi cff aflx uvindex \
               `ls | grep $1.ph_rate | sed s/$1.//`
 do
    if [ -f $1.$suffix ] 

sbin/resrm

@@ -2,7 +2,7 @@
 #
 # Script to remove all files of results with a given base name.
 #
-for suffix in uflx dflx sflx nflx vflx hrts radn cfi cff aflx
+for suffix in uflx dflx sflx nflx vflx hrts radn cfi cff aflx uvindex
 do
    if [ -f $1.$suffix ]
       then rm -f $1.$suffix

src/aux/l_run_cdf.F90

@@ -170,6 +170,8 @@ PROGRAM l_run_cdf
 !       Contribution function (to the upwards intensity at TOA)
   REAL  (RealK), ALLOCATABLE :: contrib_func_f(:,:,:)
 !       Contribution function (to the outgoing flux at TOA)
+  REAL  (RealK), ALLOCATABLE :: uv_index(:,:,:)
+!       UV-index
 
   REAL  (RealK), ALLOCATABLE :: layer_heat_capacity(:,:)
 !       Heat capacity of layers
@@ -334,6 +336,11 @@ PROGRAM l_run_cdf
   control%l_flux_direct_clear_sph_band = .FALSE.
   control%l_flux_down_clear_band = .FALSE.
   control%l_flux_up_clear_band = .FALSE.
+  control%l_uv_index = spectrum%basic%l_present(17) &
+                 .AND. spectrum%basic%l_present(2)
+  control%l_uv_index_surf = control%l_uv_index
+  control%l_uv_index_clear = .FALSE.
+  control%l_uv_index_clear_surf = .FALSE.
   control%l_actinic_flux = spectrum%basic%l_present(2)
   control%l_actinic_flux_clear = .FALSE.
   control%l_actinic_flux_band = .FALSE.
@@ -411,6 +418,12 @@ PROGRAM l_run_cdf
               dimen%nd_channel)                                       )
   ALLOCATE( heating_rate(dimen%nd_profile, dimen%nd_layer,              &
               dimen%nd_channel)                                       )
+  IF (control%l_uv_index) THEN
+    ALLOCATE( uv_index(dimen%nd_profile, 0: dimen%nd_layer,             &
+              dimen%nd_channel)                                       )
+  ELSE
+    ALLOCATE( uv_index(1, 1, 1) )
+  END IF
   IF (control%l_actinic_flux) THEN
     ALLOCATE( actinic_flux(dimen%nd_profile, dimen%nd_layer,            &
                 dimen%nd_channel)                                     )
@@ -1493,6 +1506,12 @@ PROGRAM l_run_cdf
       END IF
     END IF
 
+    IF (control%l_uv_index_surf) THEN
+      WRITE(iu_stdout, "(a,(f5.2))") 'UV-index:', radout%uv_index_surf
+    END IF
+    IF (control%l_uv_index) THEN
+      uv_index = radout%uv_index
+    END IF
     IF (control%l_actinic_flux) THEN
       actinic_flux = radout%actinic_flux
     END IF
@@ -1514,8 +1533,8 @@ PROGRAM l_run_cdf
       atm%p, atm%p_level,                                               &
       n_channel,                                                        &
       flux_total, flux_diffuse_down,                                    &
-      radout%flux_up, flux_net,                                         &
-      flux_direct, actinic_flux, photolysis_rate, heating_rate,         &
+      radout%flux_up, flux_net, flux_direct,                            &
+      uv_index, actinic_flux, photolysis_rate, heating_rate,            &
       contrib_func_i, contrib_func_f,                                   &
       dimen%nd_profile, npd_latitude, npd_longitude, dimen%nd_layer,    &
       dimen%nd_channel,                                                 &

src/aux/l_run_cdl.F90

@@ -262,6 +262,8 @@ PROGRAM l_run_cdl
   control%l_flux_direct_clear_sph_band = .FALSE.
   control%l_flux_down_clear_band = .FALSE.
   control%l_flux_up_clear_band = .FALSE.
+  control%l_uv_index = .FALSE.
+  control%l_uv_index_clear = .FALSE.
   control%l_actinic_flux = .FALSE.
   control%l_actinic_flux_clear = .FALSE.
   control%l_actinic_flux_band = .FALSE.

src/aux/output_flux_cdf.f

@@ -10,9 +10,6 @@
 !   This subroutine receives fluxes and heating rates as input and
 !   calls routines to write them to netCDF-files.
 !
-! Method:
-!   Straightforward.
-!
 !- ---------------------------------------------------------------------
       SUBROUTINE output_flux_cdf(ierr
      &  , control, sp
@@ -21,7 +18,7 @@ SUBROUTINE output_flux_cdf(ierr
      &  , n_layer, name_vert_coord, len_vert_coord, p, p_level
      &  , n_channel
      &  , flux_down, flux_down_diffuse, flux_up, flux_net, flux_direct
-     &  , actinic_flux, photolysis_rate, heating_rate
+     &  , uv_index, actinic_flux, photolysis_rate, heating_rate
      &  , contrib_funci, contrib_funcf
      &  , nd_profile, nd_latitude, nd_longitude, nd_layer, nd_channel
      &  , nd_cdl_dimen, nd_cdl_dimen_size, nd_cdl_data, nd_cdl_var
@@ -114,6 +111,8 @@ SUBROUTINE output_flux_cdf(ierr
 !           Net fluxes
      &  , flux_direct(nd_profile, 0: nd_layer, nd_channel)
 !           Direct fluxes
+     &  , uv_index(nd_profile, 0: nd_layer, nd_channel)
+!           UV-index
      &  , actinic_flux(nd_profile, nd_layer, nd_channel)
 !           Actinic fluxes
      &  , photolysis_rate(nd_profile, nd_layer,
@@ -488,8 +487,40 @@ SUBROUTINE output_flux_cdf(ierr
         IF (ierr /= i_normal) RETURN
 !
       ENDIF
-!
-!
+
+      IF (control%l_uv_index) THEN
+!       UV-index:
+        file_name(1: length_name+1+len_file_suffix)
+     &    =base_name(1: length_name)//'.'
+     &    //phys_suffix(IP_uv_index)
+        var_name(n_var)='uvindex'
+        var_type(n_var)='float'
+        var_unit(n_var)='None'
+        var_long(n_var)='UV-index'
+        n_data(n_var)=n_profile*(n_layer+1)*n_channel
+        DO k=1, n_channel
+          DO i=1, n_layer+1
+            DO l=1, n_profile
+              point=l+(i-1+(k-1)*(n_layer+1))*n_profile
+              data_fl(point, n_var)=uv_index(l, i-1, k)
+            ENDDO
+          ENDDO
+        ENDDO
+        CALL write_cdf(ierr
+     &    , file_name(1: length_name+1+len_file_suffix)
+     &    , nd_cdl_dimen, nd_cdl_dimen_size, nd_cdl_data, nd_cdl_var
+     &    , n_dimension, dimension_name, dimension_type
+     &    , dimension_unit
+     &    , dimension_long, dimension_size
+     &    , dimension_array_int, dimension_array_fl
+     &    , n_var, var_name, var_type, var_unit, var_long
+     &    , n_dimension_var, list_dimension_var
+     &    , n_data, data_int, data_fl
+     &    )
+        IF (ierr /= i_normal) RETURN
+      END IF
+
+
 !     Heating Rates:
       file_name(1: length_name+1+len_file_suffix)
      &  =base_name(1: length_name)//'.'//phys_suffix(IP_heating_rate)

src/interface_core/socrates_def_diag.F90

@@ -139,6 +139,30 @@ module socrates_def_diag
 real(RealExt) :: wavelength_uv = 3.2e-07_RealExt
 ! Upper wavelength boundary for the UV flux diagnostics, defaults to 320nm
 
+real(RealExt), pointer :: uv_index(:, :) => null()
+! UV-index (n_profile, 0:n_layer)
+
+real(RealExt), pointer :: uv_index_clear(:, :) => null()
+! Clear-sky UV-index (n_profile, 0:n_layer)
+
+real(RealExt), pointer :: uv_index_surf(:) => null()
+! UV-index at the surface (n_profile)
+
+real(RealExt), pointer :: uv_index_clear_surf(:) => null()
+! Clear-sky UV-index at the surface (n_profile)
+
+real(RealExt), pointer :: uv_index_clean(:, :) => null()
+! Clean-air UV-index (n_profile, 0:n_layer)
+
+real(RealExt), pointer :: uv_index_clear_clean(:, :) => null()
+! Clear-clean UV-index (n_profile, 0:n_layer)
+
+real(RealExt), pointer :: uv_index_clean_surf(:) => null()
+! Clean-air UV-index at the surface (n_profile)
+
+real(RealExt), pointer :: uv_index_clear_clean_surf(:) => null()
+! Clear-clean UV-index at the surface (n_profile)
+
 real(RealExt), pointer :: flux_direct_toa(:) => null()
 ! Direct downwards flux at top-of-atmosphere, Wm-2 (n_profile)
 

src/interface_core/socrates_runes.F90

@@ -1024,6 +1024,13 @@ subroutine runes(n_profile, n_layer, diag, &
   control_clean%l_photolysis_rate = control%l_photolysis_rate_clean
   control_clean%l_photolysis_rate_clear = control%l_photolysis_rate_clear_clean
 
+  control_clean%l_uv_index = associated(diag%uv_index_clean)
+  control_clean%l_uv_index_clear = associated(diag%uv_index_clear_clean)
+  control_clean%l_uv_index_surf = associated(diag%uv_index_clean_surf)
+  control_clean%l_uv_index_clear_surf = associated( &
+             diag%uv_index_clear_clean_surf)
+
+  control_clean%l_blue_flux_surf          = .false.
   control_clean%l_aerosol_absorption_band = .false.
   control_clean%l_aerosol_scattering_band = .false.
   control_clean%l_aerosol_asymmetry_band  = .false.
@@ -1045,6 +1052,8 @@ subroutine runes(n_profile, n_layer, diag, &
     associated(diag%flux_up_toa) .or. &
     associated(diag%flux_up_tile) .or. &
     control%l_clear .or. &
+    control%l_uv_index .or. &
+    control%l_uv_index_surf .or. &
     control%l_blue_flux_surf .or. &
     control%l_aerosol_absorption_band .or. &
     control%l_aerosol_scattering_band .or. &

src/interface_core/socrates_set_control.F90

@@ -131,6 +131,11 @@ subroutine set_control(control, diag, spectrum, l_set_defaults, &
 control%l_flux_down_clear_band = associated(diag%flux_down_clear_band)
 control%l_flux_up_clear_band = associated(diag%flux_up_clear_band)
 
+control%l_uv_index = associated(diag%uv_index)
+control%l_uv_index_clear = associated(diag%uv_index_clear)
+control%l_uv_index_surf = associated(diag%uv_index_surf)
+control%l_uv_index_clear_surf = associated(diag%uv_index_clear_surf)
+
 control%l_photolysis_rate = associated(diag%photolysis_rate)
 do k=1, spectrum%photol%n_pathway
   if (associated(diag%photolysis(k)%rate)) &
@@ -155,6 +160,8 @@ subroutine set_control(control, diag, spectrum, l_set_defaults, &
     control%l_flux_direct_clear_band        .or. &
     control%l_flux_down_clear_band          .or. &
     control%l_flux_up_clear_band            .or. &
+    control%l_uv_index_clear                .or. &
+    control%l_uv_index_clear_surf           .or. &
     control%l_photolysis_rate_clear) then
   control%l_clear = .true.
 end if
@@ -171,6 +178,8 @@ subroutine set_control(control, diag, spectrum, l_set_defaults, &
     associated(diag%flux_direct_clean_band) .or. &
     associated(diag%flux_down_clean_band)   .or. &
     associated(diag%flux_up_clean_band)     .or. &
+    associated(diag%uv_index_clean)         .or. &
+    associated(diag%uv_index_clean_surf)    .or. &
     control%l_photolysis_rate_clean) then
   control%l_clean = .true.
 end if
@@ -187,6 +196,8 @@ subroutine set_control(control, diag, spectrum, l_set_defaults, &
     associated(diag%flux_direct_clear_clean_band) .or. &
     associated(diag%flux_down_clear_clean_band)   .or. &
     associated(diag%flux_up_clear_clean_band)     .or. &
+    associated(diag%uv_index_clear_clean)         .or. &
+    associated(diag%uv_index_clear_clean_surf)    .or. &
     control%l_photolysis_rate_clear_clean) then
   control%l_clear_clean = .true.
 end if

src/interface_core/socrates_set_diag.F90

@@ -315,6 +315,35 @@ subroutine set_diag(diag, control, dimen, spectrum, &
 call sum_surf_weight_bands(diag%flux_up_uv_clear_surf, &
                            radout%flux_up_clear_band)
 
+call sum_flux_channels(diag%uv_index, radout%uv_index)
+call sum_flux_channels(diag%uv_index_clear, radout%uv_index_clear)
+call sum_flux_channels(diag%uv_index_clean, radout_clean%uv_index)
+call sum_flux_channels(diag%uv_index_clear_clean, radout_clean%uv_index_clear)
+if (associated(diag%uv_index_surf)) then
+  do l=1, n_profile
+    diag%uv_index_surf(list(l)) &
+      = real(radout%uv_index_surf(l), RealExt)
+  end do
+end if
+if (associated(diag%uv_index_clear_surf)) then
+  do l=1, n_profile
+    diag%uv_index_clear_surf(list(l)) &
+      = real(radout%uv_index_clear_surf(l), RealExt)
+  end do
+end if
+if (associated(diag%uv_index_clean_surf)) then
+  do l=1, n_profile
+    diag%uv_index_clean_surf(list(l)) &
+      = real(radout_clean%uv_index_surf(l), RealExt)
+  end do
+end if
+if (associated(diag%uv_index_clear_clean_surf)) then
+  do l=1, n_profile
+    diag%uv_index_clear_clean_surf(list(l)) &
+      = real(radout_clean%uv_index_clear_surf(l), RealExt)
+  end do
+end if
+
 !------------------------------------------------------------------------------
 ! Photolysis rates
 !------------------------------------------------------------------------------