Skip to content

Commit

Permalink
Add tuv photolysis data.
Browse files Browse the repository at this point in the history
  • Loading branch information
Mingzhao Liu committed Sep 16, 2023
1 parent 67c67ee commit e13c7e9
Show file tree
Hide file tree
Showing 4 changed files with 195 additions and 7 deletions.
Binary file added data/photolysis_rate.nc
View file
Binary file not shown.
125 changes: 125 additions & 0 deletions src/libtrac.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -359,6 +359,59 @@ double clim_zm(

/*****************************************************************************/

double phot_rate(
double rate[CP][50][30],
phot_t * phot,
double p,
double sza,
double o3c ){

/* Check pressure range... */
double p_help = p;
if (p < phot->p[phot->np - 1])
p_help = phot->p[phot->np - 1];
else if (p > phot->p[0])
p_help = phot->p[0];

/* Check sza range... */
double sza_help = sza;
if (sza < phot->sza[0])
sza_help = phot->sza[0];
else if (sza > phot->sza[phot->nsza - 1])
sza_help = phot->sza[phot->nsza - 1];

/* Check ozone column range... */
double o3c_help = o3c;
if (o3c < phot->o3c[0])
o3c_help = phot->o3c[0];
else if (o3c > phot->o3c[phot->no3c - 1])
o3c_help = phot->o3c[phot->no3c - 1];

/* Get indices... */
int ip = locate_irr(phot->p, phot->np, p_help);
int isza = locate_reg(phot->sza, phot->nsza, sza_help);
int io3c = locate_reg(phot->o3c, phot->no3c, o3c_help);

/* Interpolate climatology data... */
double aux00 = LIN(phot->p[ip], rate[ip][isza][io3c],
phot->p[ip + 1], rate[ip + 1][isza][io3c], p_help);
double aux01 = LIN(phot->p[ip], rate[ip][isza][io3c + 1],
phot->p[ip + 1], rate[ip + 1][isza][io3c + 1], p_help);
double aux10 = LIN(phot->p[ip], rate[ip][isza + 1][io3c],
phot->p[ip + 1], rate[ip + 1][isza + 1][io3c], p_help);
double aux11 = LIN(phot->p[ip], rate[ip][isza + 1][io3c + 1],
phot->p[ip + 1], rate[ip + 1][isza + 1][io3c + 1],
p_help);

aux00 = LIN(phot->o3c[io3c], aux00, phot->o3c[io3c + 1], aux01, o3c_help);
aux11 = LIN(phot->o3c[io3c], aux10, phot->o3c[io3c + 1], aux11, o3c_help);
aux00 = LIN(phot->sza[isza], aux00, phot->sza[isza + 1], aux11, sza_help);

return GSL_MAX(aux00, 0.0);
}

/*****************************************************************************/

void compress_pack(
char *varname,
float *array,
Expand Down Expand Up @@ -2175,6 +2228,78 @@ void read_clim_zm(

/*****************************************************************************/

void read_photol(
phot_t * phot) {

int ncid, varid, ip, is, io;

double *help1, *help2, *help3, *help4;

char *filename="../../data/photolysis_rate.nc";

/* Write info... */
LOG(1, "Read photolysis rate data: %s", filename);

/* Open netCDF file... */
if (nc_open(filename, NC_NOWRITE, &ncid) != NC_NOERR) {
WARN("Photolysis rate data are missing!");
return;
}

/* Read pressure data... */
NC_INQ_DIM("press", &phot->np, 2, CP);
NC_GET_DOUBLE("press", phot->p, 1);

/* Check ordering of pressure data... */
if (phot->p[0] < phot->p[1])
ERRMSG("Pressure data are not descending!");

/* Read latitudes... */
NC_INQ_DIM("total_o3col", &phot->no3c, 2, 30);
NC_GET_DOUBLE("total_o3col", phot->o3c, 1);

if (phot->o3c[0] > phot->o3c[1])
ERRMSG("Latitude data are not ascending!");

NC_INQ_DIM("sza", &phot->nsza, 2, 50);
NC_GET_DOUBLE("sza", phot->sza, 1);

/* Check ordering of latitude data... */
if (phot->sza[0] > phot->sza[1])
ERRMSG("Latitude data are not ascending!");

/* Read data... */
ALLOC(help1, double,
phot->np * phot->nsza * phot->no3c);
ALLOC(help2, double,
phot->np * phot->nsza * phot->no3c);
ALLOC(help3, double,
phot->np * phot->nsza * phot->no3c);
ALLOC(help4, double,
phot->np * phot->nsza * phot->no3c);
NC_GET_DOUBLE("n2o", help1, 1);
NC_GET_DOUBLE("ccl4", help2, 1);
NC_GET_DOUBLE("cfc11", help3, 1);
NC_GET_DOUBLE("cfc12", help4, 1);
for (ip = 0; ip < phot->np; ip++)
for (is = 0; is < phot->nsza; is++)
for (io = 0; io < phot->no3c; io++) {
phot->n2o[ip][is][io] = help1[ARRAY_3D(ip, is, phot->nsza, io, phot->no3c)];
phot->ccl4[ip][is][io] = help2[ARRAY_3D(ip, is, phot->nsza, io, phot->no3c)];
phot->ccl3f[ip][is][io] = help3[ARRAY_3D(ip, is, phot->nsza, io, phot->no3c)];
phot->ccl2f2[ip][is][io] = help4[ARRAY_3D(ip, is, phot->nsza, io, phot->no3c)];
}
free(help1);
free(help2);
free(help3);
free(help4);

/* Close netCDF file... */
NC(nc_close(ncid));
}

/*****************************************************************************/

void read_ctl(
const char *filename,
int argc,
Expand Down
50 changes: 50 additions & 0 deletions src/libtrac.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1779,6 +1779,42 @@ typedef struct {

} clim_t;

/*! Photolysis rate. */
typedef struct{

/*! Dimension number of pressure levels. */
int np;

/*! Dimension number of solar zenith angles. */
int nsza;

/*! Dimension number of total ozone column. */
int no3c;

/*! Pressure [hPa]. */
double p[CP];

/*! Solar zenith angle [deg]. */
double sza[50];

/*! Total Ozone column [DU]. */
double o3c[30];

/*! N2O photolysis rate. */
double n2o[CP][50][30];

/*! N2O photolysis rate. */
double ccl4[CP][50][30];

/*! N2O photolysis rate. */
double ccl3f[CP][50][30];

/*! N2O photolysis rate. */
double ccl2f2[CP][50][30];

} phot_t;


/*! Meteo data. */
typedef struct {

Expand Down Expand Up @@ -2009,6 +2045,16 @@ double clim_zm(
const double lat,
const double p);

#ifdef _OPENACC
#pragma acc routine (clim_zm)
#endif
double phot_rate(
double rate[CP][50][30],
phot_t * phot,
double p,
double sza,
double o3c );

/*! Pack or unpack array. */
void compress_pack(
char *varname,
Expand Down Expand Up @@ -2371,6 +2417,10 @@ void read_clim_zm(
char *varname,
clim_zm_t * zm);

/*! Read Photolysis rates. */
void read_photol(
phot_t * phot);

/*! Read control parameters. */
void read_ctl(
const char *filename,
Expand Down
27 changes: 20 additions & 7 deletions src/trac.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -208,6 +208,7 @@ void module_h2o2_chem(
void module_tracer_chem(
ctl_t * ctl,
clim_t * clim,
phot_t * phot,
met_t * met0,
met_t * met1,
atm_t * atm,
Expand Down Expand Up @@ -300,6 +301,8 @@ int main(

ctl_t ctl;

phot_t phot;

atm_t *atm;

cache_t *cache;
Expand Down Expand Up @@ -395,6 +398,9 @@ int main(
/* Read climatological data... */
read_clim(&ctl, clim);

/* Read photolysis rate data... */
read_photol(&phot);

/* Read atmospheric data... */
sprintf(filename, "%s/%s", dirname, argv[3]);
if (!read_atm(filename, &ctl, atm))
Expand Down Expand Up @@ -563,7 +569,7 @@ int main(

/* First-order tracer chemistry... */
if (ctl.tracer_chem)
module_tracer_chem(&ctl, clim, met0, met1, atm, dt);
module_tracer_chem(&ctl, clim, &phot, met0, met1, atm, dt);

/* KPP chemistry... */
if (ctl.kpp_chem) {
Expand Down Expand Up @@ -1780,6 +1786,7 @@ void module_h2o2_chem(
void module_tracer_chem(
ctl_t * ctl,
clim_t * clim,
phot_t * phot,
met_t * met0,
met_t * met1,
atm_t * atm,
Expand All @@ -1795,9 +1802,10 @@ void module_tracer_chem(
double sza = sza_calc(atm->time[ip], atm->lon[ip], atm->lat[ip]);

/* Get temperature... */
double t;
double t, o3c;
INTPOL_INIT;
INTPOL_3D(t, 1);
INTPOL_2D(o3c, 1)

/* Get molecular density... */
double M = MOLEC_DENS(atm->p[ip], t);
Expand All @@ -1808,28 +1816,32 @@ void module_tracer_chem(
/* Reactions for CFC-10... */
if (ctl->qnt_Cccl4 >= 0) {
double K_o1d = ARRHENIUS(3.30e-10, 0, t) * o1d * M;
double K_hv = ROETH_PHOTOL(7.79e-07, 6.32497, 0.75857, sza);
// double K_hv = ROETH_PHOTOL(7.79e-07, 6.32497, 0.75857, sza);
double K_hv = phot_rate(phot->ccl4, phot, atm->p[ip], sza, o3c);
atm->q[ctl->qnt_Cccl4][ip] *= exp(-dt[ip] * (K_hv + K_o1d));
}

/* Reactions for CFC-11... */
if (ctl->qnt_Cccl3f >= 0) {
double K_o1d = ARRHENIUS(2.30e-10, 0, t) * o1d * M;
double K_hv = ROETH_PHOTOL(6.79e-07, 6.25031, 0.75941, sza);
// double K_hv = ROETH_PHOTOL(6.79e-07, 6.25031, 0.75941, sza);
double K_hv = phot_rate(phot->ccl3f, phot, atm->p[ip], sza, o3c);
atm->q[ctl->qnt_Cccl3f][ip] *= exp(-dt[ip] * (K_hv + K_o1d));
}

/* Reactions for CFC-12... */
if (ctl->qnt_Cccl2f2 >= 0) {
double K_o1d = ARRHENIUS(1.40e-10, -25, t) * o1d * M;
double K_hv = ROETH_PHOTOL(2.81e-08, 6.47452, 0.75909, sza);
// double K_hv = ROETH_PHOTOL(2.81e-08, 6.47452, 0.75909, sza);
double K_hv = phot_rate(phot->ccl2f2, phot, atm->p[ip], sza, o3c);
atm->q[ctl->qnt_Cccl2f2][ip] *= exp(-dt[ip] * (K_hv + K_o1d));
}

/* Reactions for N2O... */
if (ctl->qnt_Cn2o >= 0) {
double K_o1d = ARRHENIUS(1.19e-10, -20, t) * o1d * M;
double K_hv = ROETH_PHOTOL(1.61e-08, 6.21077, 0.76015, sza);
// double K_hv = ROETH_PHOTOL(1.61e-08, 6.21077, 0.76015, sza);
double K_hv = phot_rate(phot->n2o, phot, atm->p[ip], sza, o3c);
atm->q[ctl->qnt_Cn2o][ip] *= exp(-dt[ip] * (K_hv + K_o1d));
}
}
Expand Down Expand Up @@ -2568,7 +2580,8 @@ void write_output(
#endif

/* Write atmospheric data... */
if (ctl->atm_basename[0] != '-' && fmod(t, ctl->atm_dt_out) == 0) {
if (ctl->atm_basename[0] != '-' && \
(fmod(t, ctl->atm_dt_out) == 0 || t == ctl->t_stop) ) {
if (ctl->atm_type_out == 0)
sprintf(ext, "tab");
else if (ctl->atm_type_out == 1)
Expand Down

0 comments on commit e13c7e9

Please sign in to comment.