src/ecckd/calc_cost_function_lw.cpp
adept::aReal
calc_cost_function_ckd_lw(const adept::Vector& pressure_hl,       ///< Pressure (Pa)
			  const adept::Matrix& planck_hl,         ///< Planck function (W m-2)
			  const adept::Vector& surf_emiss_orig,   ///< Surface emissivity
			  const adept::Vector& surf_planck,       ///< Surface spectral Planck function (W m-2)
			  const adept::aMatrix& optical_depth,    ///< Optical depth of gases
			  const adept::Matrix& flux_dn,           ///< True downwelling flux (W m-2)
			  const adept::Matrix& flux_up,           ///< True upwelling flux (W m-2)
			  const adept::Matrix& hr,                ///< True heating rate (K s-1)
			  const adept::Vector& spectral_flux_dn_surf, ///< g-point surface downward flux (W m-2)
			  const adept::Vector& spectral_flux_up_toa,  ///< g-point TOA upward flux (W m-2)
			  adept::Real flux_weight,                ///< Weight applied to TOA and surface fluxes
			  adept::Real flux_profile_weight,        ///< Weight applied to other fluxes
			  adept::Real broadband_weight,           ///< Weight of broadband vs spectral (0-1)
			  adept::Real spectral_boundary_weight,   ///< Weight of spectral boundary fluxes
			  const adept::Vector& layer_weight,      ///< Weight applied to heating rates in each layer
			  adept::Matrix* relative_ckd_flux_dn,    ///< Subtract relative-to flux dn, if not NULL
			  adept::Matrix* relative_ckd_flux_up,    ///< Subtract relative-to flux up, if not NULL
			  const adept::intVector& band_mapping)
{
  using namespace adept;

  // Convert K s-1 to K day-1
  static const Real hr_weight = 3600.0*24.0;

  int nlay = pressure_hl.size()-1;
  int ng   = optical_depth.dimension(1);

  aMatrix flux_dn_fwd_orig(nlay+1,ng);
  aMatrix flux_up_fwd_orig(nlay+1,ng);
  Vector surf_emissivity;

  if (band_mapping.empty()) {
    surf_emissivity = surf_emiss_orig;
  }
  else {
    surf_emissivity = surf_emiss_orig(band_mapping);
  }

  radiative_transfer_lw(planck_hl,
			optical_depth,
			surf_emissivity,
			surf_planck,
			flux_dn_fwd_orig,
			flux_up_fwd_orig);

  // Fluxes are to be computed relative to a reference scenario
  if (relative_ckd_flux_dn) {
    flux_dn_fwd_orig -= *relative_ckd_flux_dn;
    flux_up_fwd_orig -= *relative_ckd_flux_up;
  }


  aMatrix flux_dn_fwd, flux_up_fwd;

  int nband = ng;
  if (band_mapping.empty()) {
    flux_dn_fwd >>= flux_dn_fwd_orig;
    flux_up_fwd >>= flux_up_fwd_orig;
  }
  else {
    nband = maxval(band_mapping)+1;
    flux_dn_fwd.resize(nlay+1,nband);
    flux_up_fwd.resize(nlay+1,nband);
    for (int iband = 0; iband < nband; ++iband) {
      intVector index = find(band_mapping == iband);
      flux_dn_fwd(__,iband) = sum(flux_dn_fwd_orig(__,index),1);
      flux_up_fwd(__,iband) = sum(flux_up_fwd_orig(__,index),1);
    }
  }

  aMatrix heating_rate_fwd(nlay,nband);
  heating_rate(pressure_hl, flux_dn_fwd, flux_up_fwd, heating_rate_fwd);

  aReal cost_fn = 0.0;
  // Spectral contribution to cost function

  // Weighting for interior fluxes 
  Vector interface_weight = flux_profile_weight * 0.5*(layer_weight(range(0,end-1))+layer_weight(range(1,end)));

  for (int iband = 0; iband < nband; ++iband) {
    cost_fn += hr_weight*hr_weight*sum(layer_weight*((heating_rate_fwd(__,iband))-hr(__,iband))
				       *(heating_rate_fwd(__,iband)-hr(__,iband)))
      + flux_weight*((flux_dn_fwd(end,iband)-flux_dn(end,iband))*(flux_dn_fwd(end,iband)-flux_dn(end,iband))
		     +(flux_up_fwd(0,iband)-flux_up(0,iband))*(flux_up_fwd(0,iband)-flux_up(0,iband)));
    if (flux_profile_weight > 0.0) {
      cost_fn += sum(interface_weight*((flux_dn_fwd(range(1,end-1),iband)-flux_dn(range(1,end-1),iband))
				       *(flux_dn_fwd(range(1,end-1),iband)-flux_dn(range(1,end-1),iband))
				       +(flux_up_fwd(range(1,end-1),iband)-flux_up(range(1,end-1),iband))
				       *(flux_up_fwd(range(1,end-1),iband)-flux_up(range(1,end-1),iband))));
    }
  }

  // Broadband contribution to cost function
  cost_fn = (cost_fn*(1.0-broadband_weight))/nband
    + broadband_weight*hr_weight*hr_weight*sum(layer_weight*(sum(heating_rate_fwd-hr,1)
							     *sum(heating_rate_fwd-hr,1)))
    + broadband_weight*flux_weight*(sum(flux_dn_fwd(end,__)-flux_dn(end,__))*sum(flux_dn_fwd(end,__)-flux_dn(end,__))
				    +sum(flux_up_fwd(0,__)-flux_up(0,__))*sum(flux_up_fwd(0,__)-flux_up(0,__)));
  //LOG << sum(flux_dn_fwd(end,__)) << " " << sum(flux_dn_fwd_orig(end,__)) << " " << sum(flux_dn(end,__)) << "\n";

  if (flux_profile_weight > 0.0) {
    aVector flux_dn_error = sum(flux_dn_fwd(range(1,end-1),__)-flux_dn(range(1,end-1),__),1);
    aVector flux_up_error = sum(flux_up_fwd(range(1,end-1),__)-flux_up(range(1,end-1),__),1);
    cost_fn += broadband_weight*sum(interface_weight*(flux_dn_error*flux_dn_error
						      +flux_up_error*flux_up_error));
  }

  if (spectral_boundary_weight > 0.0 && !spectral_flux_dn_surf.empty() && !spectral_flux_up_toa.empty()) {
    cost_fn += spectral_boundary_weight
      * sum(  (flux_dn_fwd_orig(end,__)-spectral_flux_dn_surf)
	     *(flux_dn_fwd_orig(end,__)-spectral_flux_dn_surf)
	    + (flux_up_fwd_orig(0  ,__)-spectral_flux_up_toa )
	     *(flux_up_fwd_orig(0  ,__)-spectral_flux_up_toa ));
  }

  return cost_fn;
}
/// Compute the cost function, in the form of the mean-squared error
/// in heating rate, of a CKD scheme
Compute the cost function, in the form of the mean-squared error
in heating rate, of a CKD scheme
Compute the cost function in the form of the mean - squared error in heating rate of a CKD scheme
adept::aReal
calc_cost_function_ckd_lw(const adept::Vector& pressure_hl,       
			  const adept::Matrix& planck_hl,         
			  const adept::Vector& surf_emiss_orig,   
			  const adept::Vector& surf_planck,       
			  const adept::aMatrix& optical_depth,    
			  const adept::Matrix& flux_dn,           
			  const adept::Matrix& flux_up,           
			  const adept::Matrix& hr,                
			  const adept::Vector& spectral_flux_dn_surf, 
			  const adept::Vector& spectral_flux_up_toa,  
			  adept::Real flux_weight,                
			  adept::Real flux_profile_weight,        
			  adept::Real broadband_weight,           
			  adept::Real spectral_boundary_weight,   
			  const adept::Vector& layer_weight,      
			  adept::Matrix* relative_ckd_flux_dn,    
			  adept::Matrix* relative_ckd_flux_up,    
			  const adept::intVector& band_mapping)
{
  using namespace adept;
  static const Real hr_weight = 3600.0*24.0;
  int nlay = pressure_hl.size()-1;
  int ng   = optical_depth.dimension(1);
  aMatrix flux_dn_fwd_orig(nlay+1,ng);
  aMatrix flux_up_fwd_orig(nlay+1,ng);
  Vector surf_emissivity;
  if (band_mapping.empty()) {
    surf_emissivity = surf_emiss_orig;
  }
  else {
    surf_emissivity = surf_emiss_orig(band_mapping);
  }
  radiative_transfer_lw(planck_hl,
			optical_depth,
			surf_emissivity,
			surf_planck,
			flux_dn_fwd_orig,
			flux_up_fwd_orig);
  if (relative_ckd_flux_dn) {
    flux_dn_fwd_orig -= *relative_ckd_flux_dn;
    flux_up_fwd_orig -= *relative_ckd_flux_up;
  }
  aMatrix flux_dn_fwd, flux_up_fwd;
  int nband = ng;
  if (band_mapping.empty()) {
    flux_dn_fwd >>= flux_dn_fwd_orig;
    flux_up_fwd >>= flux_up_fwd_orig;
  }
  else {
    nband = maxval(band_mapping)+1;
    flux_dn_fwd.resize(nlay+1,nband);
    flux_up_fwd.resize(nlay+1,nband);
    for (int iband = 0; iband < nband; ++iband) {
      intVector index = find(band_mapping == iband);
      flux_dn_fwd(__,iband) = sum(flux_dn_fwd_orig(__,index),1);
      flux_up_fwd(__,iband) = sum(flux_up_fwd_orig(__,index),1);
    }
  }
  aMatrix heating_rate_fwd(nlay,nband);
  heating_rate(pressure_hl, flux_dn_fwd, flux_up_fwd, heating_rate_fwd);
  aReal cost_fn = 0.0;
  Vector interface_weight = flux_profile_weight * 0.5*(layer_weight(range(0,end-1))+layer_weight(range(1,end)));
  for (int iband = 0; iband < nband; ++iband) {
    cost_fn += hr_weight*hr_weight*sum(layer_weight*((heating_rate_fwd(__,iband))-hr(__,iband))
				       *(heating_rate_fwd(__,iband)-hr(__,iband)))
      + flux_weight*((flux_dn_fwd(end,iband)-flux_dn(end,iband))*(flux_dn_fwd(end,iband)-flux_dn(end,iband))
		     +(flux_up_fwd(0,iband)-flux_up(0,iband))*(flux_up_fwd(0,iband)-flux_up(0,iband)));
    if (flux_profile_weight > 0.0) {
      cost_fn += sum(interface_weight*((flux_dn_fwd(range(1,end-1),iband)-flux_dn(range(1,end-1),iband))
				       *(flux_dn_fwd(range(1,end-1),iband)-flux_dn(range(1,end-1),iband))
				       +(flux_up_fwd(range(1,end-1),iband)-flux_up(range(1,end-1),iband))
				       *(flux_up_fwd(range(1,end-1),iband)-flux_up(range(1,end-1),iband))));
    }
  }
  cost_fn = (cost_fn*(1.0-broadband_weight))/nband
    + broadband_weight*hr_weight*hr_weight*sum(layer_weight*(sum(heating_rate_fwd-hr,1)
							     *sum(heating_rate_fwd-hr,1)))
    + broadband_weight*flux_weight*(sum(flux_dn_fwd(end,__)-flux_dn(end,__))*sum(flux_dn_fwd(end,__)-flux_dn(end,__))
				    +sum(flux_up_fwd(0,__)-flux_up(0,__))*sum(flux_up_fwd(0,__)-flux_up(0,__)));
  if (flux_profile_weight > 0.0) {
    aVector flux_dn_error = sum(flux_dn_fwd(range(1,end-1),__)-flux_dn(range(1,end-1),__),1);
    aVector flux_up_error = sum(flux_up_fwd(range(1,end-1),__)-flux_up(range(1,end-1),__),1);
    cost_fn += broadband_weight*sum(interface_weight*(flux_dn_error*flux_dn_error
						      +flux_up_error*flux_up_error));
  }
  if (spectral_boundary_weight > 0.0 && !spectral_flux_dn_surf.empty() && !spectral_flux_up_toa.empty()) {
    cost_fn += spectral_boundary_weight
      * sum(  (flux_dn_fwd_orig(end,__)-spectral_flux_dn_surf)
	     *(flux_dn_fwd_orig(end,__)-spectral_flux_dn_surf)
	    + (flux_up_fwd_orig(0  ,__)-spectral_flux_up_toa )
	     *(flux_up_fwd_orig(0  ,__)-spectral_flux_up_toa ));
  }
  return cost_fn;
}
adept :: aReal calc_cost_function_ckd_lw ( const adept :: Vector & pressure_hl , const adept :: Matrix & planck_hl , const adept :: Vector & surf_emiss_orig , const adept :: Vector & surf_planck , const adept :: aMatrix & optical_depth , const adept :: Matrix & flux_dn , const adept :: Matrix & flux_up , const adept :: Matrix & hr , const adept :: Vector & spectral_flux_dn_surf , const adept :: Vector & spectral_flux_up_toa , adept :: Real flux_weight , adept :: Real flux_profile_weight , adept :: Real broadband_weight , adept :: Real spectral_boundary_weight , const adept :: Vector & layer_weight , adept :: Matrix * relative_ckd_flux_dn , adept :: Matrix * relative_ckd_flux_up , const adept :: intVector & band_mapping ) { using namespace adept ; static const Real hr_weight = 3600.0 * 24.0 ; int nlay = pressure_hl . size ( ) - 1 ; int ng = optical_depth . dimension ( 1 ) ; aMatrix flux_dn_fwd_orig ( nlay + 1 , ng ) ; aMatrix flux_up_fwd_orig ( nlay + 1 , ng ) ; Vector surf_emissivity ; if ( band_mapping . empty ( ) ) { surf_emissivity = surf_emiss_orig ; } else { surf_emissivity = surf_emiss_orig ( band_mapping ) ; } radiative_transfer_lw ( planck_hl , optical_depth , surf_emissivity , surf_planck , flux_dn_fwd_orig , flux_up_fwd_orig ) ; if ( relative_ckd_flux_dn ) { flux_dn_fwd_orig -= * relative_ckd_flux_dn ; flux_up_fwd_orig -= * relative_ckd_flux_up ; } aMatrix flux_dn_fwd , flux_up_fwd ; int nband = ng ; if ( band_mapping . empty ( ) ) { flux_dn_fwd >>= flux_dn_fwd_orig ; flux_up_fwd >>= flux_up_fwd_orig ; } else { nband = maxval ( band_mapping ) + 1 ; flux_dn_fwd . resize ( nlay + 1 , nband ) ; flux_up_fwd . resize ( nlay + 1 , nband ) ; for ( int iband = 0 ; iband < nband ; ++ iband ) { intVector index = find ( band_mapping == iband ) ; flux_dn_fwd ( __ , iband ) = sum ( flux_dn_fwd_orig ( __ , index ) , 1 ) ; flux_up_fwd ( __ , iband ) = sum ( flux_up_fwd_orig ( __ , index ) , 1 ) ; } } aMatrix heating_rate_fwd ( nlay , nband ) ; heating_rate ( pressure_hl , flux_dn_fwd , flux_up_fwd , heating_rate_fwd ) ; aReal cost_fn = 0.0 ; Vector interface_weight = flux_profile_weight * 0.5 * ( layer_weight ( range ( 0 , end - 1 ) ) + layer_weight ( range ( 1 , end ) ) ) ; for ( int iband = 0 ; iband < nband ; ++ iband ) { cost_fn += hr_weight * hr_weight * sum ( layer_weight * ( ( heating_rate_fwd ( __ , iband ) ) - hr ( __ , iband ) ) * ( heating_rate_fwd ( __ , iband ) - hr ( __ , iband ) ) ) + flux_weight * ( ( flux_dn_fwd ( end , iband ) - flux_dn ( end , iband ) ) * ( flux_dn_fwd ( end , iband ) - flux_dn ( end , iband ) ) + ( flux_up_fwd ( 0 , iband ) - flux_up ( 0 , iband ) ) * ( flux_up_fwd ( 0 , iband ) - flux_up ( 0 , iband ) ) ) ; if ( flux_profile_weight > 0.0 ) { cost_fn += sum ( interface_weight * ( ( flux_dn_fwd ( range ( 1 , end - 1 ) , iband ) - flux_dn ( range ( 1 , end - 1 ) , iband ) ) * ( flux_dn_fwd ( range ( 1 , end - 1 ) , iband ) - flux_dn ( range ( 1 , end - 1 ) , iband ) ) + ( flux_up_fwd ( range ( 1 , end - 1 ) , iband ) - flux_up ( range ( 1 , end - 1 ) , iband ) ) * ( flux_up_fwd ( range ( 1 , end - 1 ) , iband ) - flux_up ( range ( 1 , end - 1 ) , iband ) ) ) ) ; } } cost_fn = ( cost_fn * ( 1.0 - broadband_weight ) ) / nband + broadband_weight * hr_weight * hr_weight * sum ( layer_weight * ( sum ( heating_rate_fwd - hr , 1 ) * sum ( heating_rate_fwd - hr , 1 ) ) ) + broadband_weight * flux_weight * ( sum ( flux_dn_fwd ( end , __ ) - flux_dn ( end , __ ) ) * sum ( flux_dn_fwd ( end , __ ) - flux_dn ( end , __ ) ) + sum ( flux_up_fwd ( 0 , __ ) - flux_up ( 0 , __ ) ) * sum ( flux_up_fwd ( 0 , __ ) - flux_up ( 0 , __ ) ) ) ; if ( flux_profile_weight > 0.0 ) { aVector flux_dn_error = sum ( flux_dn_fwd ( range ( 1 , end - 1 ) , __ ) - flux_dn ( range ( 1 , end - 1 ) , __ ) , 1 ) ; aVector flux_up_error = sum ( flux_up_fwd ( range ( 1 , end - 1 ) , __ ) - flux_up ( range ( 1 , end - 1 ) , __ ) , 1 ) ; cost_fn += broadband_weight * sum ( interface_weight * ( flux_dn_error * flux_dn_error + flux_up_error * flux_up_error ) ) ; } if ( spectral_boundary_weight > 0.0 && ! spectral_flux_dn_surf . empty ( ) && ! spectral_flux_up_toa . empty ( ) ) { cost_fn += spectral_boundary_weight * sum ( ( flux_dn_fwd_orig ( end , __ ) - spectral_flux_dn_surf ) * ( flux_dn_fwd_orig ( end , __ ) - spectral_flux_dn_surf ) + ( flux_up_fwd_orig ( 0 , __ ) - spectral_flux_up_toa ) * ( flux_up_fwd_orig ( 0 , __ ) - spectral_flux_up_toa ) ) ; } return cost_fn ; }
Compute the cost function, in the form of the mean-squared error
in heating rate, of a CKD scheme
Compute the cost function in the form of the mean - squared error in heating rate of a CKD scheme
///< Pressure (Pa) ///< Planck function (W m-2) ///< Surface emissivity ///< Surface spectral Planck function (W m-2) ///< Optical depth of gases ///< True downwelling flux (W m-2) ///< True upwelling flux (W m-2) ///< True heating rate (K s-1) ///< g-point surface downward flux (W m-2) ///< g-point TOA upward flux (W m-2) ///< Weight applied to TOA and surface fluxes ///< Weight applied to other fluxes ///< Weight of broadband vs spectral (0-1) ///< Weight of spectral boundary fluxes ///< Weight applied to heating rates in each layer ///< Subtract relative-to flux dn, if not NULL ///< Subtract relative-to flux up, if not NULL // Convert K s-1 to K day-1 // Fluxes are to be computed relative to a reference scenario // Spectral contribution to cost function // Weighting for interior fluxes  // Broadband contribution to cost function //LOG << sum(flux_dn_fwd(end,__)) << " " << sum(flux_dn_fwd_orig(end,__)) << " " << sum(flux_dn(end,__)) << "\n";