SOCOLv3 bug fixes

Collection of bug fixes for various SOCOLv3 versions incl. SOCOL-AER after Nov 2020:

SOCOL/SOCOL-AER: bug fix in interactive wet deposition scheme (Mar 2023)

A wrong cloud cover variable is used in the convective scavenging of gas-phase species:

 
...
!!$                   mcov(jl)            = zoldcov(i)
                   mcov(jl)            = zwashfrac(i)
...

New code: mo_socol_wetdep_main.f90

SOCOL/SOCOL-AER: bugs in interactive wet deposition scheme (Nov 2022)

SOCOL-AER did not run on 48 CPUs anymore

By change I realized that SOCOL-AER does not run on 48 CPUs (not sure if this problem is new and related to new compiler revisions etc. or if nobody run SOCOL-AER on 48 CPUs in the past...)

Problem: shape of argument SVI in call to socol_wet_deposition (src/socol_mezon.f90), replace SVI(1:kproma,:,:,krow) [or SVI(1:kproma,1:klev,1:n_aer_bin,krow) ] by SVI(:,:,:,krow)

New code:

 
...
     if (linteractive_wetdep) CALL socol_wet_deposition(krow, kproma, pxtm1, pxtte, t, p, ph, paclc, &
                                   plwc, piwc, pfrain_no, pfsnow_no, prate_r,prate_s, pr_cover, & !large scale properties
                                   pmflxr, pmflxs, cv_precnew, cv_snownew, zmfu, kconbot, & !convective properties
                                   SVI(:,:,:,krow))
...

Problems in calculation of convective cloud cover

For runs with interactive wet deposition, the model crashed from time to time with the following error: Numerical Conversion Not Representable

Reason was a bug in the calculation of the convective cloud cover in mo_socol_wetdep_main.f90:

New code: mo_socol_wetdep_main.f90

SOCOL-AER: Reading volcanic SO2 emissions T${RES}_so2_emiss_volcYYYYMM (Nov 2022)

Error message src/socol_boundary_conditions.f90, array so2_erupt: Index 31 of dim 1 beyond upper bound of 29

Problem occurs at the beginning of February, as the model run is resumed on 31 Jan 23:30:00

Workaround: Add if-statement to src/socol_boundary_conditions.f90, subroutine set_socol_bcond_fluxes:

 
...
    fac = 1._dp/(atomweight*amso2) ! [kg/m^2/s] -> [molec/m^2/s]
    call get_date_components(current_date,yr,mo,dy)

    if (dy .gt. size(so2_erupt,1)) goto 100
...
       ENDDO

100    continue

    ! SO2 aircraft emissions:
...

Ozone decoupling: lo3_coupl = .FALSE., but lsrb_lya_heating = .TRUE. (Jun 2022)

The previous bug fix allowed to use an ozone climatology provided on pressure levels only. The current fix allows now to use all types of possible ozone climatologies for the calculation of the mesospheric heating rates.

New code [modules/mo_socol_sun.f90, subroutine srb_lya_heating]:

 
    USE mo_radiation,       ONLY: io3
    USE mo_o3_lwb,          ONLY: o3_lwb
    USE mo_o3clim,          ONLY: o3clim, o3clim_rawlev
...
    REAL(dp) :: zdp(kbdim,klev)
...
   ! g/g -> mol/mol
    oz: SELECT CASE (io3)
    CASE (0)
       vmr_o3(1:kproma,:) = EPSILON(1._dp)
    CASE (2)
     ! layer pressure thickness
       zdp(1:kproma,:)=paphm1(1:kproma,2:klev+1)-paphm1(1:kproma,1:klev)
       vmr_o3(1:kproma,:) = o3_lwb(krow,zdp,paphm1)*amd/amo3
    CASE (3)
       vmr_o3(1:kproma,:) = o3clim(krow,kproma,kbdim,klev,paphm1,papm1)*amd/amo3
    CASE (4)
       vmr_o3(1:kproma,:) = o3clim(krow,kproma,kbdim,klev,paphm1,papm1)*amd/amo3
    CASE (5)
       vmr_o3(1:kproma,:) = o3clim_rawlev(krow,kproma,kbdim,klev)*amd/amo3
    CASE (10)   ! coupling with CTM
       vmr_o3(:,:)=xtm1(:,:,idt_o3,krow)
    CASE default
       CALL finish('srb_lya_heating','o3: this "io3" is not supported')
    END SELECT oz

Ozone decoupling: lo3_coupl = .FALSE., but lsrb_lya_heating = .TRUE. (Jul 2021)

The extra heating rates due to absorption in the Schumann Runge bands (175-200 nm) and by the Lyman alpha line (121.6 nm) depend on the ozone field. In case the coupling between chemistry and radiation via ozone is switched off (lo3_coupl = .FALSE.), the ozone climatology used for the radiative calculations has to be used in the subroutine srb_lya_heating as well! For that, an additional if-statement has to added at the beginning of the subroutine. Details see below.

Old code [modules/mo_socol_sun.f90, subroutine srb_lya_heating]: There are different code versions available, some include already a first bug fix. Either:

    vmr_o3(:,:)=xtm1(:,:,idt_o3,krow)

or:

    USE mo_o3clim        , ONLY: o3clim
...
    if (lchem) then
      vmr_o3(:,:)=xtm1(:,:,idt_o3,krow)
    else
      ! g/g -> mol/mol
      vmr_o3(1:kproma,:) = o3clim(krow,kproma,kbdim,klev,paphm1,papm1)*amd/amo3
    end if

New code [modules/mo_socol_sun.f90, subroutine srb_lya_heating]:

    USE mo_o3clim        , ONLY: o3clim
...
    if (lchem) then
      if (lo3_coupl) then
         vmr_o3(:,:)=xtm1(:,:,idt_o3,krow)
      else
         ! g/g -> mol/mol
         vmr_o3(1:kproma,:) = o3clim(krow,kproma,kbdim,klev,paphm1,papm1)*amd/amo3
      end if
    else
      ! g/g -> mol/mol
      vmr_o3(1:kproma,:) = o3clim(krow,kproma,kbdim,klev,paphm1,papm1)*amd/amo3
    end if

Problems with LSOLARVAR = .FALSE. (Dec 2020)

The SOCOL namelist parameter LSOLARVAR determines, whether the solar input (irradiance, additional heating by Schumann-Runge-, Hartley-, Huggins- and Lyman-alpha-radiation, NO photolysis) varies with time or not. In the latter case monthly mean values averaged over two consecutive solar cycles (1977-1998) are used. These values are set in mo_socol_namelist.f90.

For CCMI-1, the parameterization for extraheating as well as the solar input files have been updated, but not the settings in mo_socol_namelist.f90 for lsolarvar = .false., at least not in all model versions. In some model versions, the settings had been updated, but the values for Schumann-Runge- and Lyman-alpha-bands have been swapped.

Furthermore, the values for the two Huggings-parameters have to be converted from W/m2 to W/m2/angstrom. This was missing for the case lsolarvar = .false.. Modified source code

And here the updated constant values for mo_socol_namelist.f90:

  REAL(dp):: sun_srb_const(12) = (/ &
       1.51E-01_dp, 1.50E-01_dp, 1.48E-01_dp, 1.46E-01_dp, &
       1.44E-01_dp, 1.42E-01_dp, 1.41E-01_dp, 1.42E-01_dp, &
       1.44E-01_dp, 1.46E-01_dp, 1.48E-01_dp, 1.50E-01_dp /)
                                           ! Monthly values of Schumann-Runge 
                                           ! parameter if lsolarvar=.FALSE.
                                           ! (Default: mean over 1977-1998)
  REAL(dp):: sun_lya_const(12) = (/ &
       7.03E-03_dp, 7.04E-03_dp, 6.92E-03_dp, 6.80E-03_dp, &
       6.66E-03_dp, 6.57E-03_dp, 6.54E-03_dp, 6.60E-03_dp, &
       6.68E-03_dp, 6.82E-03_dp, 6.91E-03_dp, 7.07E-03_dp /)
                                           ! Monthly values of Lyman alpha
                                           ! parameter if lsolarvar=.FALSE.
                                           ! (Default: mean over 1977-1998)
  REAL(dp):: sun_har_const(12) = (/ &
       3.28E-01_dp, 3.08E-01_dp, 2.52E-01_dp, 1.77E-01_dp, &
       1.02E-01_dp, 4.68E-02_dp, 2.71E-02_dp, 4.77E-02_dp, &
       1.02E-01_dp, 1.77E-01_dp, 2.52E-01_dp, 3.09E-01_dp /)
                                           ! Monthly values of Hartley
                                           ! parameter if lsolarvar=.FALSE.
                                           ! (Default: mean over 1977-1998)
  REAL(dp):: sun_hug1_const(12) = (/ &
       7.15E-01_dp, 6.67E-01_dp, 5.41E-01_dp, 3.69E-01_dp, &
       1.97E-01_dp, 7.16E-02_dp, 2.68E-02_dp, 7.30E-02_dp, &
       1.97E-01_dp, 3.68E-01_dp, 5.40E-01_dp, 6.69E-01_dp/)
                                           ! Monthly values of Huggins 1
                                           ! parameter if lsolarvar=.FALSE.
                                           ! (Default: mean over 1977-1998)
  REAL(dp):: sun_hug2_const(12) = (/ &
       3.11E+00_dp, 2.90E+00_dp, 2.34E+00_dp, 1.57E+00_dp, &
       8.06E-01_dp, 2.49E-01_dp, 4.92E-02_dp, 2.53E-01_dp, &
       8.06E-01_dp, 1.57E+00_dp, 2.33E+00_dp, 2.90E+00_dp/)
                                           ! Monthly values of Huggins 2
                                           ! parameter if lsolarvar=.FALSE.
                                           ! (Default: mean over 1977-1998)
  REAL(dp):: sun_nopho_const(12) = (/ &
       4.67E-06_dp, 4.66E-06_dp, 4.59E-06_dp, 4.52E-06_dp, &
       4.44E-06_dp, 4.39E-06_dp, 4.37E-06_dp, 4.39E-06_dp, &
       4.45E-06_dp, 4.52E-06_dp, 4.59E-06_dp, 4.66E-06_dp /)
                                           ! Monthly values of NO-photolysis 
                                           ! parameter if lsolarvar=.FALSE.
                                           ! (Default: mean over 1977-1998)

SOCOL-AER: aerosol-radiation coupling (Dec 2020)

Due to an erroneous do loop over vertical levels, the extinction coefficients for AER aerosols had not been set to zero for the lower most 6 model levels. Furthermore, the extinctions had not been converted to optical depth. That leads to a double-counting of AOD from sulfate aerosols in the lowermost 6 levels:

• Short bug report
• Fixed source code