Ignore:
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:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_ifs_rrtm.F90

    r4773 r4853  
    4141
    4242    use radiation_config
     43    use radiation_spectral_definition, only &
     44         &  : SolarReferenceTemperature, TerrestrialReferenceTemperature
    4345
    4446    type(config_type), intent(inout), target :: config
    … …  
    6567          &   44, 107, 94, 14, 108, 15, 16, 109, 17, 18, 110, 111, 112 &
    6668          & /)
     69
     70    logical :: do_sw, do_lw
    6771   
    6872    real(jphook) :: hook_handle
    … …  
    7781    if (lhook) call dr_hook('radiation_ifs_rrtm:setup_gas_optics',0,hook_handle)
    7882
     83    do_sw = (config%do_sw .and. config%i_gas_model_sw == IGasModelIFSRRTMG)
     84    do_lw = (config%do_lw .and. config%i_gas_model_lw == IGasModelIFSRRTMG)
     85   
    7986    ! The IFS implementation of RRTMG uses many global variables.  In
    8087    ! the IFS these will have been set up already; otherwise set them
    … …  
    8592      call SURRTPK
    8693      call SURRTRF
    87       call RRTM_INIT_140GP(directory)
    88       call SRTM_INIT(directory)
     94      if (do_lw) then
     95        call RRTM_INIT_140GP(directory)
     96      end if
     97      if (do_sw) then
     98        call SRTM_INIT(directory)
     99      end if
    89100    end if
    90101
    91     ! Cloud and aerosol properties can only be defined per band
    92     config%do_cloud_aerosol_per_sw_g_point = .false.
    93     config%do_cloud_aerosol_per_lw_g_point = .false.
    94 
    95     config%n_g_sw = jpgsw
    96     config%n_g_lw = jpglw
    97     config%n_bands_sw = 14
    98     config%n_bands_lw = 16
    99 
    100     ! Wavenumber ranges of each band may be needed so that the user
    101     ! can compute UV and photosynthetically active radiation for a
    102     ! particular wavelength range
    103     call config%gas_optics_sw%spectral_def%allocate_bands_only( &
    104          &  [2600.0_jprb, 3250.0_jprb, 4000.0_jprb, 4650.0_jprb, 5150.0_jprb, 6150.0_jprb, 7700.0_jprb, &
    105          &   8050.0_jprb, 12850.0_jprb, 16000.0_jprb, 22650.0_jprb, 29000.0_jprb, 38000.0_jprb, 820.0_jprb], &
    106          &  [3250.0_jprb, 4000.0_jprb, 4650.0_jprb, 5150.0_jprb, 6150.0_jprb, 7700.0_jprb, 8050.0_jprb, &
    107          &   12850.0_jprb, 16000.0_jprb, 22650.0_jprb, 29000.0_jprb, 38000.0_jprb, 50000.0_jprb, 2600.0_jprb])
    108     call config%gas_optics_lw%spectral_def%allocate_bands_only( &
    109          &  [10.0_jprb, 350.0_jprb, 500.0_jprb, 630.0_jprb, 700.0_jprb, 820.0_jprb, 980.0_jprb, 1080.0_jprb, &
    110          &   1180.0_jprb, 1390.0_jprb, 1480.0_jprb, 1800.0_jprb, 2080.0_jprb, 2250.0_jprb, 2380.0_jprb, 2600.0_jprb], &
    111          &  [350.0_jprb, 500.0_jprb, 630.0_jprb, 700.0_jprb, 820.0_jprb, 980.0_jprb, 1080.0_jprb, 1180.0_jprb, &
    112          &   1390.0_jprb, 1480.0_jprb, 1800.0_jprb, 2080.0_jprb, 2250.0_jprb, 2380.0_jprb, 2600.0_jprb, 3250.0_jprb])
    113 
    114     allocate(config%i_band_from_g_sw          (config%n_g_sw))
    115     allocate(config%i_band_from_g_lw          (config%n_g_lw))
    116     allocate(config%i_band_from_reordered_g_sw(config%n_g_sw))
    117     allocate(config%i_band_from_reordered_g_lw(config%n_g_lw))
    118     allocate(config%i_g_from_reordered_g_sw(config%n_g_sw))
    119     allocate(config%i_g_from_reordered_g_lw(config%n_g_lw))
    120 
    121     ! Shortwave starts at 16: need to start at 1
    122     config%i_band_from_g_sw = ngb_sw - ngb_sw(1)+1
    123     config%i_band_from_g_lw = ngb_lw
    124 
    125     if (config%i_solver_sw == ISolverSpartacus) then
    126       ! SPARTACUS requires g points ordered in approximately
    127       ! increasing order of optical depth
    128       config%i_g_from_reordered_g_sw = RRTM_GPOINT_REORDERING_SW
    129     else
    130       ! Implied-do for no reordering
    131 !      config%i_g_from_reordered_g_sw = RRTM_GPOINT_REORDERING_SW
    132       config%i_g_from_reordered_g_sw = (/ (irep, irep=1,config%n_g_sw) /)
     102    if (do_sw) then
     103     
     104      ! Cloud and aerosol properties can only be defined per band
     105      config%do_cloud_aerosol_per_sw_g_point = .false.
     106      config%n_g_sw = jpgsw
     107      config%n_bands_sw = 14
     108      ! Wavenumber ranges of each band may be needed so that the user
     109      ! can compute UV and photosynthetically active radiation for a
     110      ! particular wavelength range
     111      call config%gas_optics_sw%spectral_def%allocate_bands_only(SolarReferenceTemperature, &
     112           &  [2600.0_jprb, 3250.0_jprb, 4000.0_jprb, 4650.0_jprb, 5150.0_jprb, 6150.0_jprb, 7700.0_jprb, &
     113           &   8050.0_jprb, 12850.0_jprb, 16000.0_jprb, 22650.0_jprb, 29000.0_jprb, 38000.0_jprb, 820.0_jprb], &
     114           &  [3250.0_jprb, 4000.0_jprb, 4650.0_jprb, 5150.0_jprb, 6150.0_jprb, 7700.0_jprb, 8050.0_jprb, &
     115           &   12850.0_jprb, 16000.0_jprb, 22650.0_jprb, 29000.0_jprb, 38000.0_jprb, 50000.0_jprb, 2600.0_jprb])
     116      allocate(config%i_band_from_g_sw          (config%n_g_sw))
     117      allocate(config%i_band_from_reordered_g_sw(config%n_g_sw))
     118      allocate(config%i_g_from_reordered_g_sw   (config%n_g_sw))
     119      ! Shortwave starts at 16: need to start at 1
     120      config%i_band_from_g_sw = ngb_sw - ngb_sw(1)+1
     121
     122      if (config%i_solver_sw == ISolverSpartacus) then
     123        ! SPARTACUS requires g points ordered in approximately
     124        ! increasing order of optical depth
     125        config%i_g_from_reordered_g_sw = RRTM_GPOINT_REORDERING_SW
     126      else
     127        ! Implied-do for no reordering
     128        !      config%i_g_from_reordered_g_sw = RRTM_GPOINT_REORDERING_SW
     129        config%i_g_from_reordered_g_sw = (/ (irep, irep=1,config%n_g_sw) /)
     130      end if
     131
     132      config%i_band_from_reordered_g_sw &
     133           = config%i_band_from_g_sw(config%i_g_from_reordered_g_sw)
     134
     135      ! The i_spec_* variables are used solely for storing spectral
     136      ! data, and this can either be by band or by g-point
     137      if (config%do_save_spectral_flux .or. config%do_toa_spectral_flux) then
     138        if (config%do_save_gpoint_flux) then
     139          config%n_spec_sw = config%n_g_sw
     140          config%i_spec_from_reordered_g_sw => config%i_g_from_reordered_g_sw
     141        else
     142          config%n_spec_sw = config%n_bands_sw
     143          config%i_spec_from_reordered_g_sw => config%i_band_from_reordered_g_sw
     144        end if
     145      else
     146        config%n_spec_sw = 0
     147        nullify(config%i_spec_from_reordered_g_sw)
     148      end if
     149     
    133150    end if
    134151
    135     if (config%i_solver_lw == ISolverSpartacus) then
    136       ! SPARTACUS requires g points ordered in approximately
    137       ! increasing order of optical depth
    138       config%i_g_from_reordered_g_lw = RRTM_GPOINT_REORDERING_LW
    139     else
    140       ! Implied-do for no reordering
    141       config%i_g_from_reordered_g_lw = (/ (irep, irep=1,config%n_g_lw) /)
     152    if (do_lw) then
     153      ! Cloud and aerosol properties can only be defined per band
     154      config%do_cloud_aerosol_per_lw_g_point = .false.
     155      config%n_g_lw = jpglw
     156      config%n_bands_lw = 16
     157      call config%gas_optics_lw%spectral_def%allocate_bands_only(TerrestrialReferenceTemperature, &
     158           &  [10.0_jprb, 350.0_jprb, 500.0_jprb, 630.0_jprb, 700.0_jprb, 820.0_jprb, 980.0_jprb, 1080.0_jprb, &
     159           &   1180.0_jprb, 1390.0_jprb, 1480.0_jprb, 1800.0_jprb, 2080.0_jprb, 2250.0_jprb, 2380.0_jprb, 2600.0_jprb], &
     160           &  [350.0_jprb, 500.0_jprb, 630.0_jprb, 700.0_jprb, 820.0_jprb, 980.0_jprb, 1080.0_jprb, 1180.0_jprb, &
     161           &   1390.0_jprb, 1480.0_jprb, 1800.0_jprb, 2080.0_jprb, 2250.0_jprb, 2380.0_jprb, 2600.0_jprb, 3250.0_jprb])
     162      allocate(config%i_band_from_g_lw          (config%n_g_lw))
     163      allocate(config%i_band_from_reordered_g_lw(config%n_g_lw))
     164      allocate(config%i_g_from_reordered_g_lw   (config%n_g_lw))
     165      config%i_band_from_g_lw = ngb_lw
     166
     167      if (config%i_solver_lw == ISolverSpartacus) then
     168        ! SPARTACUS requires g points ordered in approximately
     169        ! increasing order of optical depth
     170        config%i_g_from_reordered_g_lw = RRTM_GPOINT_REORDERING_LW
     171      else
     172        ! Implied-do for no reordering
     173        config%i_g_from_reordered_g_lw = (/ (irep, irep=1,config%n_g_lw) /)
     174      end if
     175
     176      config%i_band_from_reordered_g_lw &
     177           = config%i_band_from_g_lw(config%i_g_from_reordered_g_lw)
     178
     179      ! The i_spec_* variables are used solely for storing spectral
     180      ! data, and this can either be by band or by g-point
     181      if (config%do_save_spectral_flux .or. config%do_toa_spectral_flux) then
     182        if (config%do_save_gpoint_flux) then
     183          config%n_spec_lw = config%n_g_lw
     184          config%i_spec_from_reordered_g_lw => config%i_g_from_reordered_g_lw
     185        else
     186          config%n_spec_lw = config%n_bands_lw
     187          config%i_spec_from_reordered_g_lw => config%i_band_from_reordered_g_lw
     188        end if
     189      else
     190        config%n_spec_lw = 0
     191        nullify(config%i_spec_from_reordered_g_lw)
     192      end if
     193
    142194    end if
    143 
    144     config%i_band_from_reordered_g_sw &
    145          = config%i_band_from_g_sw(config%i_g_from_reordered_g_sw)
    146 
    147     config%i_band_from_reordered_g_lw &
    148          = config%i_band_from_g_lw(config%i_g_from_reordered_g_lw)
    149 
    150     ! The i_spec_* variables are used solely for storing spectral
    151     ! data, and this can either be by band or by g-point
    152     if (config%do_save_spectral_flux .or. config%do_toa_spectral_flux) then
    153       if (config%do_save_gpoint_flux) then
    154         config%n_spec_sw = config%n_g_sw
    155         config%n_spec_lw = config%n_g_lw
    156         config%i_spec_from_reordered_g_sw => config%i_g_from_reordered_g_sw
    157         config%i_spec_from_reordered_g_lw => config%i_g_from_reordered_g_lw
    158       else
    159         config%n_spec_sw = config%n_bands_sw
    160         config%n_spec_lw = config%n_bands_lw
    161         config%i_spec_from_reordered_g_sw => config%i_band_from_reordered_g_sw
    162         config%i_spec_from_reordered_g_lw => config%i_band_from_reordered_g_lw
    163       end if
    164     else
    165       config%n_spec_sw = 0
    166       config%n_spec_lw = 0
    167       nullify(config%i_spec_from_reordered_g_sw)
    168       nullify(config%i_spec_from_reordered_g_lw)
    169     end if
    170 
     195   
    171196    if (lhook) call dr_hook('radiation_ifs_rrtm:setup_gas_optics',1,hook_handle)
    172197
    … …  
    201226    use yomhook  , only : lhook, dr_hook, jphook
    202227
    203     use radiation_config,         only : config_type, ISolverSpartacus
     228    use radiation_config,         only : config_type, ISolverSpartacus, IGasModelIFSRRTMG
    204229    use radiation_thermodynamics, only : thermodynamics_type
    205230    use radiation_single_level,   only : single_level_type
    … …  
    334359    integer :: istartlev, iendlev
    335360
     361    logical :: do_sw, do_lw
     362   
    336363    integer :: jlev, jgreorder, jg, ig, iband, jcol
    337364
    … …  
    345372
    346373    if (lhook) call dr_hook('radiation_ifs_rrtm:gas_optics',0,hook_handle)
     374
     375    do_sw = (config%do_sw .and. config%i_gas_model_sw == IGasModelIFSRRTMG)
     376    do_lw = (config%do_lw .and. config%i_gas_model_lw == IGasModelIFSRRTMG)
    347377
    348378    ! Compute start and end levels for indexing the gas mixing ratio
    … …  
    392422         &  ZPAVEL , ZTAVEL , ZPZ , ZTZ, IREFLECT) 
    393423
    394     CALL RRTM_SETCOEF_140GP &
    395          &( istartcol, iendcol, nlev , ZCOLDRY  , ZWBRODL , ZWKL , &
    396          &  ZFAC00 , ZFAC01   , ZFAC10 , ZFAC11 , ZFORFAC,ZFORFRAC,INDFOR, JP, JT, JT1 , &
    397          &  ZCOLH2O, ZCOLCO2  , ZCOLO3 , ZCOLN2O, ZCOLCH4, ZCOLO2,ZCO2MULT , ZCOLBRD, &
    398          &  ILAYTROP,ILAYSWTCH, ILAYLOW, ZPAVEL , ZTAVEL , ZSELFFAC, ZSELFFRAC, INDSELF, &
    399          &  INDMINOR,ZSCALEMINOR,ZSCALEMINORN2,ZMINORFRAC,&
    400          &  ZRAT_H2OCO2, ZRAT_H2OCO2_1, ZRAT_H2OO3, ZRAT_H2OO3_1, &
    401          &  ZRAT_H2ON2O, ZRAT_H2ON2O_1, ZRAT_H2OCH4, ZRAT_H2OCH4_1, &
    402          &  ZRAT_N2OCO2, ZRAT_N2OCO2_1, ZRAT_O3CO2, ZRAT_O3CO2_1)   
    403 
    404     ZTAUAERL(istartcol:iendcol,:,:) = 0.0_jprb
    405 
    406     CALL RRTM_GAS_OPTICAL_DEPTH &
    407          &( istartcol, iendcol, nlev, ZOD_LW, ZPAVEL, ZCOLDRY, ZCOLBRD, ZWX ,&
    408          &  ZTAUAERL, ZFAC00 , ZFAC01, ZFAC10 , ZFAC11 , ZFORFAC,ZFORFRAC,INDFOR, &
    409          &  JP, JT, JT1, ZONEMINUS ,&
    410          &  ZCOLH2O , ZCOLCO2, ZCOLO3, ZCOLN2O, ZCOLCH4, ZCOLO2,ZCO2MULT ,&
    411          &  ILAYTROP, ILAYSWTCH,ILAYLOW, ZSELFFAC, ZSELFFRAC, INDSELF, ZPFRAC, &
    412          &  INDMINOR,ZSCALEMINOR,ZSCALEMINORN2,ZMINORFRAC,&
    413          &  ZRAT_H2OCO2, ZRAT_H2OCO2_1, ZRAT_H2OO3, ZRAT_H2OO3_1, &
    414          &  ZRAT_H2ON2O, ZRAT_H2ON2O_1, ZRAT_H2OCH4, ZRAT_H2OCH4_1, &
    415          &  ZRAT_N2OCO2, ZRAT_N2OCO2_1, ZRAT_O3CO2, ZRAT_O3CO2_1)     
    416 
    417     if (present(lw_albedo)) then
    418    
    419       call planck_function_atmos(nlev, istartcol, iendcol, config, &
    420            &                     thermodynamics, ZPFRAC, planck_hl)
    421 
    422       if (single_level%is_simple_surface) then
    423         call planck_function_surf(istartcol, iendcol, config, &
    424              &                    single_level%skin_temperature, ZPFRAC(:,:,1), &
    425              &                    lw_emission)
     424    if (do_lw) then
     425   
     426      CALL RRTM_SETCOEF_140GP &
     427           &( istartcol, iendcol, nlev , ZCOLDRY  , ZWBRODL , ZWKL , &
     428           &  ZFAC00 , ZFAC01   , ZFAC10 , ZFAC11 , ZFORFAC,ZFORFRAC,INDFOR, JP, JT, JT1 , &
     429           &  ZCOLH2O, ZCOLCO2  , ZCOLO3 , ZCOLN2O, ZCOLCH4, ZCOLO2,ZCO2MULT , ZCOLBRD, &
     430           &  ILAYTROP,ILAYSWTCH, ILAYLOW, ZPAVEL , ZTAVEL , ZSELFFAC, ZSELFFRAC, INDSELF, &
     431           &  INDMINOR,ZSCALEMINOR,ZSCALEMINORN2,ZMINORFRAC,&
     432           &  ZRAT_H2OCO2, ZRAT_H2OCO2_1, ZRAT_H2OO3, ZRAT_H2OO3_1, &
     433           &  ZRAT_H2ON2O, ZRAT_H2ON2O_1, ZRAT_H2OCH4, ZRAT_H2OCH4_1, &
     434           &  ZRAT_N2OCO2, ZRAT_N2OCO2_1, ZRAT_O3CO2, ZRAT_O3CO2_1)   
     435
     436      ZTAUAERL(istartcol:iendcol,:,:) = 0.0_jprb
     437
     438      CALL RRTM_GAS_OPTICAL_DEPTH &
     439           &( istartcol, iendcol, nlev, ZOD_LW, ZPAVEL, ZCOLDRY, ZCOLBRD, ZWX ,&
     440           &  ZTAUAERL, ZFAC00 , ZFAC01, ZFAC10 , ZFAC11 , ZFORFAC,ZFORFRAC,INDFOR, &
     441           &  JP, JT, JT1, ZONEMINUS ,&
     442           &  ZCOLH2O , ZCOLCO2, ZCOLO3, ZCOLN2O, ZCOLCH4, ZCOLO2,ZCO2MULT ,&
     443           &  ILAYTROP, ILAYSWTCH,ILAYLOW, ZSELFFAC, ZSELFFRAC, INDSELF, ZPFRAC, &
     444           &  INDMINOR,ZSCALEMINOR,ZSCALEMINORN2,ZMINORFRAC,&
     445           &  ZRAT_H2OCO2, ZRAT_H2OCO2_1, ZRAT_H2OO3, ZRAT_H2OO3_1, &
     446           &  ZRAT_H2ON2O, ZRAT_H2ON2O_1, ZRAT_H2OCH4, ZRAT_H2OCH4_1, &
     447           &  ZRAT_N2OCO2, ZRAT_N2OCO2_1, ZRAT_O3CO2, ZRAT_O3CO2_1)     
     448   
     449      if (present(lw_albedo)) then
     450   
     451        call planck_function_atmos(nlev, istartcol, iendcol, config, &
     452             &                     thermodynamics, ZPFRAC, planck_hl)
     453
     454        if (single_level%is_simple_surface) then
     455          call planck_function_surf(istartcol, iendcol, config, &
     456               &                    single_level%skin_temperature, ZPFRAC(:,:,1), &
     457               &                    lw_emission)
     458         
     459          ! The following can be used to extract the parameters defined at
     460          ! the top of the planck_function routine below:
     461          !write(*,'(a,140(e12.5,","),a)') 'ZPFRAC_surf=[', &
     462          !&  sum(ZPFRAC(istartcol:iendcol,:,1),1) / (iendcol+1-istartcol), ']'
    426463       
    427         ! The following can be used to extract the parameters defined at
    428         ! the top of the planck_function routine below:
    429         !write(*,'(a,140(e12.5,","),a)') 'ZPFRAC_surf=[', &
    430         !&  sum(ZPFRAC(istartcol:iendcol,:,1),1) / (iendcol+1-istartcol), ']'
    431        
    432         ! lw_emission at this point is actually the planck function of
    433         ! the surface
    434         lw_emission = lw_emission * (1.0_jprb - lw_albedo)
     464          ! lw_emission at this point is actually the planck function of
     465          ! the surface
     466          lw_emission = lw_emission * (1.0_jprb - lw_albedo)
     467        else
     468          ! Longwave emission has already been computed
     469          if (config%use_canopy_full_spectrum_lw) then
     470            lw_emission = transpose(single_level%lw_emission(istartcol:iendcol,:))
     471          else
     472            lw_emission = transpose(single_level%lw_emission(istartcol:iendcol, &
     473                 & config%i_emiss_from_band_lw(config%i_band_from_reordered_g_lw)))
     474          end if
     475        end if
     476
     477      end if
     478
     479      if (config%i_solver_lw == ISolverSpartacus) then
     480        ! We need to rearrange the gas optics info in memory: reordering
     481        ! the g points in order of approximately increasing optical
     482        ! depth (for efficient 3D processing on only the regions of the
     483        ! spectrum that are optically thin for gases) and reorder in
     484        ! pressure since the the functions above treat pressure
     485        ! decreasing with increasing index.  Note that the output gas
     486        ! arrays have dimensions in a different order to the inputs,
     487        ! so there is some inefficiency here.
     488        do jgreorder = 1,config%n_g_lw
     489          iband = config%i_band_from_reordered_g_lw(jgreorder)
     490          ig = config%i_g_from_reordered_g_lw(jgreorder)
     491         
     492          ! Top-of-atmosphere half level
     493          do jlev = 1,nlev
     494            do jcol = istartcol,iendcol
     495              ! Some g points can return negative optical depths;
     496              ! specifically original g points 54-56 which causes
     497              ! unphysical single-scattering albedo when combined with
     498              ! aerosol
     499              od_lw(jgreorder,jlev,jcol) &
     500                   &   = max(config%min_gas_od_lw, ZOD_LW(ig,nlev+1-jlev,jcol))
     501            end do
     502          end do
     503        end do
    435504      else
    436         ! Longwave emission has already been computed
    437         if (config%use_canopy_full_spectrum_lw) then
    438           lw_emission = transpose(single_level%lw_emission(istartcol:iendcol,:))
    439         else
    440           lw_emission = transpose(single_level%lw_emission(istartcol:iendcol, &
    441                & config%i_emiss_from_band_lw(config%i_band_from_reordered_g_lw)))
    442         end if
     505        ! G points have not been reordered
     506        do jcol = istartcol,iendcol
     507          do jlev = 1,nlev
     508            ! Check for negative optical depth
     509            od_lw(:,jlev,jcol) = max(config%min_gas_od_lw, ZOD_LW(:,nlev+1-jlev,jcol))
     510          end do
     511        end do
    443512      end if
    444513
    445514    end if
    446515
    447     if (config%i_solver_lw == ISolverSpartacus) then
    448       ! We need to rearrange the gas optics info in memory: reordering
    449       ! the g points in order of approximately increasing optical
    450       ! depth (for efficient 3D processing on only the regions of the
    451       ! spectrum that are optically thin for gases) and reorder in
    452       ! pressure since the the functions above treat pressure
    453       ! decreasing with increasing index.  Note that the output gas
    454       ! arrays have dimensions in a different order to the inputs,
    455       ! so there is some inefficiency here.
    456       do jgreorder = 1,config%n_g_lw
    457         iband = config%i_band_from_reordered_g_lw(jgreorder)
    458         ig = config%i_g_from_reordered_g_lw(jgreorder)
    459        
    460         ! Top-of-atmosphere half level
    461         do jlev = 1,nlev
    462           do jcol = istartcol,iendcol
    463             ! Some g points can return negative optical depths;
    464             ! specifically original g points 54-56 which causes
    465             ! unphysical single-scattering albedo when combined with
    466             ! aerosol
    467             od_lw(jgreorder,jlev,jcol) &
    468                  &   = max(config%min_gas_od_lw, ZOD_LW(ig,nlev+1-jlev,jcol))
     516    if (do_sw) then
     517   
     518      CALL SRTM_SETCOEF &
     519           & ( istartcol, iendcol, nlev,&
     520           & ZPAVEL  , ZTAVEL,&
     521           & ZCOLDRY , ZWKL,&
     522           & ILAYTROP,&
     523           & ZCOLCH4  , ZCOLCO2 , ZCOLH2O , ZCOLMOL  , ZCOLO2 , ZCOLO3,&
     524           & ZFORFAC , ZFORFRAC , INDFOR  , ZSELFFAC, ZSELFFRAC, INDSELF, &
     525           & ZFAC00  , ZFAC01   , ZFAC10  , ZFAC11,&
     526           & JP      , JT       , JT1     , single_level%cos_sza(istartcol:iendcol)  &
     527           & ) 
     528   
     529      ! SRTM_GAS_OPTICAL_DEPTH will not initialize profiles when the sun
     530      ! is below the horizon, so we do it here
     531      ZOD_SW(istartcol:iendcol,:,:)  = 0.0_jprb
     532      ZSSA_SW(istartcol:iendcol,:,:) = 0.0_jprb
     533      ZINCSOL(istartcol:iendcol,:)   = 0.0_jprb
     534
     535      CALL SRTM_GAS_OPTICAL_DEPTH &
     536           &( istartcol, iendcol , nlev  , ZONEMINUS_ARRAY,&
     537           & single_level%cos_sza(istartcol:iendcol), ILAYTROP,&
     538           & ZCOLCH4 , ZCOLCO2  , ZCOLH2O, ZCOLMOL , ZCOLO2   , ZCOLO3,&
     539           & ZFORFAC , ZFORFRAC , INDFOR , ZSELFFAC, ZSELFFRAC, INDSELF,&
     540           & ZFAC00  , ZFAC01   , ZFAC10 , ZFAC11  ,&
     541           & JP      , JT       , JT1    ,&
     542           & ZOD_SW  , ZSSA_SW  , ZINCSOL )
     543   
     544      ! Scale the incoming solar per band, if requested
     545      if (config%use_spectral_solar_scaling) then
     546        do jg = 1,JPGPT_SW
     547          do jcol = istartcol,iendcol
     548            ZINCSOL(jcol,jg) = ZINCSOL(jcol,jg) * &
     549                 &   single_level%spectral_solar_scaling(config%i_band_from_reordered_g_sw(jg))
    469550          end do
    470551        end do
    471       end do
    472     else
    473       ! G points have not been reordered
    474       do jcol = istartcol,iendcol
    475         do jlev = 1,nlev
    476           ! Check for negative optical depth
    477           od_lw(:,jlev,jcol) = max(config%min_gas_od_lw, ZOD_LW(:,nlev+1-jlev,jcol))
    478         end do
    479       end do
    480     end if
    481    
    482     CALL SRTM_SETCOEF &
    483          & ( istartcol, iendcol, nlev,&
    484          & ZPAVEL  , ZTAVEL,&
    485          & ZCOLDRY , ZWKL,&
    486          & ILAYTROP,&
    487          & ZCOLCH4  , ZCOLCO2 , ZCOLH2O , ZCOLMOL  , ZCOLO2 , ZCOLO3,&
    488          & ZFORFAC , ZFORFRAC , INDFOR  , ZSELFFAC, ZSELFFRAC, INDSELF, &
    489          & ZFAC00  , ZFAC01   , ZFAC10  , ZFAC11,&
    490          & JP      , JT       , JT1     , single_level%cos_sza(istartcol:iendcol)  &
    491          & ) 
    492    
    493     ! SRTM_GAS_OPTICAL_DEPTH will not initialize profiles when the sun
    494     ! is below the horizon, so we do it here
    495     ZOD_SW(istartcol:iendcol,:,:)  = 0.0_jprb
    496     ZSSA_SW(istartcol:iendcol,:,:) = 0.0_jprb
    497     ZINCSOL(istartcol:iendcol,:)   = 0.0_jprb
    498 
    499     CALL SRTM_GAS_OPTICAL_DEPTH &
    500          &( istartcol, iendcol , nlev  , ZONEMINUS_ARRAY,&
    501          & single_level%cos_sza(istartcol:iendcol), ILAYTROP,&
    502          & ZCOLCH4 , ZCOLCO2  , ZCOLH2O, ZCOLMOL , ZCOLO2   , ZCOLO3,&
    503          & ZFORFAC , ZFORFRAC , INDFOR , ZSELFFAC, ZSELFFRAC, INDSELF,&
    504          & ZFAC00  , ZFAC01   , ZFAC10 , ZFAC11  ,&
    505          & JP      , JT       , JT1    ,&
    506          & ZOD_SW  , ZSSA_SW  , ZINCSOL )
    507    
    508     ! Scale the incoming solar per band, if requested
    509     if (config%use_spectral_solar_scaling) then
    510       do jg = 1,JPGPT_SW
    511         do jcol = istartcol,iendcol
    512           ZINCSOL(jcol,jg) = ZINCSOL(jcol,jg) * &
    513             &   single_level%spectral_solar_scaling(config%i_band_from_reordered_g_sw(jg))
    514         end do
    515       end do
    516     end if
    517 
    518     ! Scaling factor to ensure that the total solar irradiance is as
    519     ! requested.  Note that if the sun is below the horizon then
    520     ! ZINCSOL will be zero.
    521     if (present(incoming_sw)) then
    522       do jcol = istartcol,iendcol
    523         if (single_level%cos_sza(jcol) > 0.0_jprb) then
     552      end if
     553
     554      ! Scaling factor to ensure that the total solar irradiance is as
     555      ! requested.  Note that if the sun is below the horizon then
     556      ! ZINCSOL will be zero.
     557      if (present(incoming_sw)) then
     558        do jcol = istartcol,iendcol
     559          if (single_level%cos_sza(jcol) > 0.0_jprb) then
    524560! Added for DWD (2020)
    525561!NEC$ nounroll
    526           incoming_sw_scale(jcol) = single_level%solar_irradiance / sum(ZINCSOL(jcol,:))
    527         else
    528           incoming_sw_scale(jcol) = 1.0_jprb
    529         end if
    530       end do
    531     end if
    532 
    533     if (config%i_solver_sw == ISolverSpartacus) then
    534 !    if (.true.) then
    535       ! Account for reordered g points
    536       do jgreorder = 1,config%n_g_sw
    537         ig = config%i_g_from_reordered_g_sw(jgreorder)
    538         do jlev = 1,nlev
    539           do jcol = istartcol,iendcol
    540             ! Check for negative optical depth
    541             od_sw (jgreorder,nlev+1-jlev,jcol) &
    542                  &  = max(config%min_gas_od_sw, ZOD_SW (jcol,jlev,ig))
    543             ssa_sw(jgreorder,nlev+1-jlev,jcol) = ZSSA_SW(jcol,jlev,ig)
    544            end do
     562            incoming_sw_scale(jcol) = single_level%solar_irradiance / sum(ZINCSOL(jcol,:))
     563          else
     564            incoming_sw_scale(jcol) = 1.0_jprb
     565          end if
    545566        end do
    546         if (present(incoming_sw)) then
    547           incoming_sw(jgreorder,:) &
    548                &  = incoming_sw_scale(:) * ZINCSOL(:,ig)
    549         end if
    550       end do
    551     else
    552       ! G points have not been reordered
    553       do jcol = istartcol,iendcol
    554         do jlev = 1,nlev
    555           do jg = 1,config%n_g_sw
    556             ! Check for negative optical depth
    557             od_sw (jg,nlev+1-jlev,jcol) = max(config%min_gas_od_sw, ZOD_SW(jcol,jlev,jg))
    558             ssa_sw(jg,nlev+1-jlev,jcol) = ZSSA_SW(jcol,jlev,jg)
    559           end do
     567      end if
     568
     569      if (config%i_solver_sw == ISolverSpartacus) then
     570!      if (.true.) then
     571        ! Account for reordered g points
     572        do jgreorder = 1,config%n_g_sw
     573          ig = config%i_g_from_reordered_g_sw(jgreorder)
     574          do jlev = 1,nlev
     575            do jcol = istartcol,iendcol
     576              ! Check for negative optical depth
     577              od_sw (jgreorder,nlev+1-jlev,jcol) &
     578                   &  = max(config%min_gas_od_sw, ZOD_SW (jcol,jlev,ig))
     579              ssa_sw(jgreorder,nlev+1-jlev,jcol) = ZSSA_SW(jcol,jlev,ig)
     580            end do
     581          end do
     582          if (present(incoming_sw)) then
     583            incoming_sw(jgreorder,:) &
     584                 &  = incoming_sw_scale(:) * ZINCSOL(:,ig)
     585          end if
    560586        end do
    561         if (present(incoming_sw)) then
    562           do jg = 1,config%n_g_sw
    563             incoming_sw(jg,jcol) = incoming_sw_scale(jcol) * ZINCSOL(jcol,jg)
    564           end do
    565         end if
    566       end do
     587      else
     588        ! G points have not been reordered
     589        do jcol = istartcol,iendcol
     590          do jlev = 1,nlev
     591            do jg = 1,config%n_g_sw
     592              ! Check for negative optical depth
     593              od_sw (jg,nlev+1-jlev,jcol) = max(config%min_gas_od_sw, ZOD_SW(jcol,jlev,jg))
     594              ssa_sw(jg,nlev+1-jlev,jcol) = ZSSA_SW(jcol,jlev,jg)
     595            end do
     596          end do
     597          if (present(incoming_sw)) then
     598            do jg = 1,config%n_g_sw
     599              incoming_sw(jg,jcol) = incoming_sw_scale(jcol) * ZINCSOL(jcol,jg)
     600            end do
     601          end if
     602        end do
     603      end if
     604
    567605    end if
    568606   
Note: See TracChangeset for help on using the changeset viewer.