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_set_spectrum.F90

@@ -135,6 +135,7 @@ subroutine set_spectrum(n_instances, spectrum, spectrum_name, spectral_file, &
 else if (present(spectrum_name).or.present(spectrum)) then
   ! DEPENDS ON: read_spectrum
   call read_spectrum(spectral_file, spec)
+  call set_weight_uv_index(spec)
   if (spec%solar%weight_blue(1) == rmdi) then
     allocate(weight_blue(spec%basic%n_band))
     call set_weight_blue(spec, weight_blue, wavelength_blue)
@@ -405,6 +406,69 @@ end function retain_absorber
 end subroutine compress_spectrum
 
 
+subroutine set_weight_uv_index(spec)
+
+use realtype_rd, only: RealK
+
+implicit none
+
+type(StrSpecData), intent(inout) :: spec
+
+! Parameters defining UV-index erythema action spectrum (CIE 1998)
+real(RealK), parameter :: k_er = 40.0_RealK  ! m2 W-1
+real(RealK), parameter :: wl_1 = 250.0_RealK ! nm
+real(RealK), parameter :: wl_2 = 298.0_RealK ! nm
+real(RealK), parameter :: wl_3 = 328.0_RealK ! nm
+real(RealK), parameter :: wl_4 = 400.0_RealK ! nm
+real(RealK), parameter :: a_1 = 0.094_RealK
+real(RealK), parameter :: a_2 = 0.015_RealK
+real(RealK), parameter :: b_1 = 298.0_RealK
+real(RealK), parameter :: b_2 = 140.0_RealK
+real(RealK), parameter :: m2nm = 1.0e9_RealK
+
+integer :: i_sub
+real(RealK) :: sub_band_min, sub_band_max, total_range
+real(RealK) :: range_min, range_max, range
+
+do i_sub=1, spec%var%n_sub_band
+  spec%var%weight_uv_index(i_sub) = 0.0_RealK
+  sub_band_min = spec%var%wavelength_sub_band(1, i_sub) * m2nm
+  sub_band_max = spec%var%wavelength_sub_band(2, i_sub) * m2nm
+  total_range = sub_band_max - sub_band_min
+
+  ! Range 1: S_er = 1, for wavelengths wl_1 - wl_2
+  range_min = max(sub_band_min, wl_1)
+  range_max = min(sub_band_max, wl_2)
+  range = range_max - range_min
+  if (range > 0.0_RealK) then
+    spec%var%weight_uv_index(i_sub) = spec%var%weight_uv_index(i_sub) &
+      + k_er * range / total_range
+  end if
+
+  ! Range 2: S_er = 10**(a_1*(b_1-wavelength)), for wavelengths wl_2 - wl_3
+  range_min = max(sub_band_min, wl_2)
+  range_max = min(sub_band_max, wl_3)
+  range = range_max - range_min
+  if (range > 0.0_RealK) then
+    spec%var%weight_uv_index(i_sub) = spec%var%weight_uv_index(i_sub) &
+      + ( 10**(a_1 * (b_1 - range_min)) - 10**(a_1 * (b_1 - range_max)) ) &
+      * k_er / ( a_1 * log(10.0_RealK) * total_range )
+  end if
+
+  ! Range 3: S_er = 10**(a_2*(b_2-wavelength)), for wavelengths wl_3 - wl_4
+  range_min = max(sub_band_min, wl_3)
+  range_max = min(sub_band_max, wl_4)
+  range = range_max - range_min
+  if (range > 0.0_RealK) then
+    spec%var%weight_uv_index(i_sub) = spec%var%weight_uv_index(i_sub) &
+      + ( 10**(a_2 * (b_2 - range_min)) - 10**(a_2 * (b_2 - range_max)) ) &
+      * k_er / ( a_2 * log(10.0_RealK) * total_range )
+  end if
+end do
+
+end subroutine set_weight_uv_index
+
+
 subroutine set_weight_blue(spec, weight_blue, wavelength_blue)
 
 use realtype_rd, only: RealK, RealExt

src/modules_gen/input_head_pcf.f90

@@ -30,7 +30,7 @@ MODULE input_head_pcf
 !   Maximum number of input profiles
   INTEGER, Parameter :: NPD_data_column      = 35
 !   Maximum number of columns of data
-  INTEGER, Parameter :: NPD_phys_type        = 62
+  INTEGER, Parameter :: NPD_phys_type        = 63
 !   Maximum number of types of data
   INTEGER, Parameter :: NPD_unit             = 20
 !   Number of physical units
@@ -182,6 +182,8 @@ MODULE input_head_pcf
 !   Actinic flux
   INTEGER, Parameter :: IP_photolysis_rate                       = 62
 !   Photolysis rate
+  INTEGER, Parameter :: IP_uv_index                              = 63
+!   UV-index
 
 !
   CHARACTER  (LEN=len_col_header), Parameter, &
@@ -201,7 +203,7 @@ MODULE input_head_pcf
     'AZIM      ', 'RADN      ', 'SRFCHR    ', 'OPWT      ', &
     'OPICE     ', 'OPSS      ', 'ISOS      ', 'GEOM      ', &
     'PHOTOL    ', 'PLEV      ', 'CFI       ', 'CFF       ', &
-    'AFLX      ', 'PHRATE    ' /) 
+    'AFLX      ', 'PHRATE    ', 'UVINDEX   ' /) 
 !   Headers for physical data
 !
   CHARACTER  (LEN=len_file_suffix), Parameter, &
@@ -221,7 +223,7 @@ MODULE input_head_pcf
     'vwazim      ', 'radn        ', 'surf        ', 'op_water    ', &
     'op_ice      ', 'ss          ', 'isos        ', 'view        ', &
     'photol      ', 'pl          ', 'cfi         ', 'cff         ', &
-    'aflx        ', 'ph_rate     '  /) 
+    'aflx        ', 'ph_rate     ', 'uvindex     '  /) 
 !   File suffixes for physical data
 !
 !
@@ -257,7 +259,8 @@ MODULE input_head_pcf
     ' Isotropic source             ', ' Viewing Geometry             ', &
     ' Rate of photolysis           ', ' Pressure on Levels           ', &
     ' Contribution function (inty) ', ' Contribution function (flux) ', &
-    ' Actinic flux                 ', ' Photolysis rate              '/)
+    ' Actinic flux                 ', ' Photolysis rate              ', &
+    ' UV-index                     '/)
 !   Long titles for physical data
 !
 !