Skip to content

Commit

Permalink
Fixed bugs in how the polar and equatorial geographies calculate tran…
Browse files Browse the repository at this point in the history 
…sitional latitudes. Modern geography is still broken.
  • Loading branch information
@RoryBarnes
RoryBarnes committed Sep 29, 2024
1 parent c9aa846 commit 3bcdcbd
Show file tree
Hide file tree
Showing 4 changed files with 89 additions and 69 deletions.
2 changes: 1 addition & 1 deletion examples/ClimateLand/equatorial.in
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
sName equatorial #name of planet
saModules poise #what vplanet modules you want to use
sGeography equitorial
dLandWaterLatitude 73.6
dLandWaterLatitude 16.4

dMass 3.00316726e-06 #mass of planet
dRadius -1.00 #radius (not important right now)
Expand Down
126 changes: 73 additions & 53 deletions src/poise.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1512,6 +1512,12 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) {
options[OPT_LANDFRACMEAN].iType = 2;
options[OPT_LANDFRACMEAN].bMultiFile = 1;
fnRead[OPT_LANDFRACMEAN] = &ReadLandFracMean;
fvFormattedString(
&options[OPT_LANDWATERLATITUDE].cLongDescr,
"The average fraction of land per latitude. Note that the actual "
"distribution of land will be normalized so that the global fraction "
"of land on the surface will equal %s",
options[OPT_LANDFRACMEAN].cName);

fvFormattedString(&options[OPT_LANDFRACAMP].cName, "dLandFracAmp");
fvFormattedString(&options[OPT_LANDFRACAMP].cDescr,
Expand Down Expand Up @@ -2091,7 +2097,7 @@ void InitializeLatGrid(BODY *body, int iBody) {
}
}

void fvInitializeLandWaterUniform(BODY *body, int iBody) {
void fvInitializeGeographyUniform(BODY *body, int iBody) {
int iLat;

for (iLat = 0; iLat < body[iBody].iNumLats; iLat++) {
Expand All @@ -2100,7 +2106,7 @@ void fvInitializeLandWaterUniform(BODY *body, int iBody) {
}
}

void fvInitializeLandWaterModern(BODY *body, int iBody) {
void fvInitializeGeographyModern(BODY *body, int iBody) {
int iLat;

for (iLat = 0; iLat < body[iBody].iNumLats; iLat++) {
Expand All @@ -2122,9 +2128,9 @@ void fvInitializeLandWaterModern(BODY *body, int iBody) {
}
}

void fvInitializeLandWaterRandom(BODY *body, int iBody) {
void fvInitializeGeographyRandom(BODY *body, int iBody) {
int iLat;
double dOffset;
double dOffset, dLandFrac, dLandFracNormFactor;

for (iLat = 0; iLat < body[iBody].iNumLats; iLat++) {
dOffset = body[iBody].dLandFracAmp *
Expand All @@ -2136,88 +2142,79 @@ void fvInitializeLandWaterRandom(BODY *body, int iBody) {
if (body[iBody].daLandFrac[iLat] < LANDFRACMIN) {
body[iBody].daLandFrac[iLat] = LANDFRACMIN;
}
}

dLandFrac = fdLandFracGlobal(body, iBody);
dLandFracNormFactor = body[iBody].dLandFracMean / dLandFrac;
for (iLat = 0; iLat < body[iBody].iNumLats; iLat++) {
body[iBody].daLandFrac[iLat] *= dLandFracNormFactor;
body[iBody].daWaterFrac[iLat] = 1.0 - body[iBody].daLandFrac[iLat];
}
}

void fvInitializeLandWaterPolar(BODY *body, int iBody) {
void fvInitializeGeographyPolar(BODY *body, int iBody) {
int iLat;

for (iLat = 0; iLat < body[iBody].iLatLandWater; iLat++) {
body[iBody].daLandFrac[iLat] = LANDFRACMAX;
body[iBody].daWaterFrac[iLat] = LANDFRACMIN;
}
for (iLat = body[iBody].iLatLandWater;
iLat < (body[iBody].iNumLats - body[iBody].iLatLandWater);
iLat++) {
body[iBody].daLandFrac[iLat] = LANDFRACMIN;
body[iBody].daWaterFrac[iLat] = LANDFRACMAX;
}
for (iLat = (body[iBody].iNumLats - body[iBody].iLatLandWater);
iLat < body[iBody].iNumLats;
iLat++) {
body[iBody].daLandFrac[iLat] = LANDFRACMIN;
body[iBody].daWaterFrac[iLat] = LANDFRACMAX;
for (iLat = 0; iLat < body[iBody].iNumLats; iLat++) {
if (body[iBody].daLats[iLat] <= -body[iBody].dLatLandWater ||
body[iBody].daLats[iLat] >= body[iBody].dLatLandWater) {
body[iBody].daLandFrac[iLat] = LANDFRACMAX;
body[iBody].daWaterFrac[iLat] = LANDFRACMIN;
} else {
body[iBody].daLandFrac[iLat] = LANDFRACMIN;
body[iBody].daWaterFrac[iLat] = LANDFRACMAX;
}
}
}

void fvInitializeLandWaterEquatorial(BODY *body, int iBody) {
void fvInitializeGeographyEquatorial(BODY *body, int iBody) {
int iLat;

for (iLat = 0; iLat < body[iBody].iLatLandWater; iLat++) {
body[iBody].daLandFrac[iLat] = LANDFRACMIN;
body[iBody].daWaterFrac[iLat] = LANDFRACMAX;
}
for (iLat = body[iBody].iLatLandWater;
iLat < (body[iBody].iNumLats - body[iBody].iLatLandWater);
iLat++) {
body[iBody].daLandFrac[iLat] = LANDFRACMAX;
body[iBody].daWaterFrac[iLat] = LANDFRACMIN;
}
for (iLat = (body[iBody].iNumLats - body[iBody].iLatLandWater);
iLat < body[iBody].iNumLats;
iLat++) {
body[iBody].daLandFrac[iLat] = LANDFRACMIN;
body[iBody].daWaterFrac[iLat] = LANDFRACMAX;
for (iLat = 0; iLat < body[iBody].iNumLats; iLat++) {
if (body[iBody].daLats[iLat] <= -body[iBody].dLatLandWater ||
body[iBody].daLats[iLat] >= body[iBody].dLatLandWater) {
body[iBody].daLandFrac[iLat] = LANDFRACMIN;
body[iBody].daWaterFrac[iLat] = LANDFRACMAX;
} else {
body[iBody].daLandFrac[iLat] = LANDFRACMAX;
body[iBody].daWaterFrac[iLat] = LANDFRACMIN;
}
}
}

void fvWriteLandFracFile(BODY *body, int iBody) {
void fvWriteGeographyFile(BODY *body, int iBody) {
int iLat;
double dInterval, dLatNow;
FILE *fp;
char *cOut = NULL;

dInterval = 180. / body[iBody].iNumLats;
fvFormattedString(&cOut, "%s.geography", body[iBody].cName);
fp = fopen(cOut, "w");
for (iLat = 0; iLat < body[iBody].iNumLats; iLat++) {
dLatNow = -90 + 0.5 * dInterval + iLat * dInterval;
fprintf(fp, "%.5lf %.5lf %.5lf\n", dLatNow, body[iBody].daLandFrac[iLat],
body[iBody].daWaterFrac[iLat]);
fprintf(fp, "%.5lf %.5lf %.5lf\n", (body[iBody].daLats[iLat] * 180. / PI),
body[iBody].daLandFrac[iLat], body[iBody].daWaterFrac[iLat]);
}
fclose(fp);
}

void InitializeLandWater(BODY *body, int iBody) {
void InitializeGeography(BODY *body, int iBody) {
int iLat;

body[iBody].daLandFrac = malloc(body[iBody].iNumLats * sizeof(double));
body[iBody].daWaterFrac = malloc(body[iBody].iNumLats * sizeof(double));

if (body[iBody].iGeography == GEOGRAPHYUNIFORM) {
fvInitializeLandWaterUniform(body, iBody);
fvInitializeGeographyUniform(body, iBody);
} else if (body[iBody].iGeography == GEOGRAPHYMODERN) {
fvInitializeLandWaterModern(body, iBody);
fvInitializeGeographyModern(body, iBody);
} else if (body[iBody].iGeography == GEOGRAPHYRANDOM) {
fvInitializeLandWaterRandom(body, iBody);
fvInitializeGeographyRandom(body, iBody);
} else if (body[iBody].iGeography == GEOGRAPHYPOLAR) {
fvInitializeLandWaterPolar(body, iBody);
fvInitializeGeographyPolar(body, iBody);
} else if (body[iBody].iGeography == GEOGRAPHYEQUATORIAL) {
fvInitializeLandWaterEquatorial(body, iBody);
fvInitializeGeographyEquatorial(body, iBody);
}

fvWriteLandFracFile(body, iBody);
fvWriteGeographyFile(body, iBody);
}

void DampTemp(BODY *body, double dTGlobalTmp, int iBody) {
Expand Down Expand Up @@ -2516,7 +2513,7 @@ void InitializeClimateParams(BODY *body, int iBody, int iVerbose) {
body[iBody].daIceAccumTot = malloc(body[iBody].iNumLats * sizeof(double));
body[iBody].daIceAblateTot = malloc(body[iBody].iNumLats * sizeof(double));

InitializeLandWater(body, iBody);
InitializeGeography(body, iBody);
body[iBody].dLatFHeatCp = 83.5; // CC sez this is about right
body[iBody].dLatentHeatIce = body[iBody].dHeatCapWater *
body[iBody].dLatFHeatCp /
Expand Down Expand Up @@ -3021,9 +3018,6 @@ void fvVerifyGeographyPolarEquatorial(BODY *body, CONTROL *control,
options[OPT_LANDWATERLATITUDE].cName);
}
}

body[iBody].iLatLandWater =
(int)(body[iBody].dLatLandWater * body[iBody].iNumLats / PI);
}

void VerifyGeography(BODY *body, CONTROL *control, OPTIONS *options,
Expand Down Expand Up @@ -4338,6 +4332,13 @@ void WriteEnergyResW(BODY *body, CONTROL *control, OUTPUT *output,
}
}

void WriteLandFracGlobal(BODY *body, CONTROL *control, OUTPUT *output,
SYSTEM *system, UNITS *units, UPDATE *update,
int iBody, double *dTmp, char **cUnit) {
*dTmp = fdLandFracGlobal(body, iBody);
fvFormattedString(cUnit, "");
}

void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) {
fvFormattedString(&output[OUT_TGLOBAL].cName, "TGlobal");
fvFormattedString(&output[OUT_TGLOBAL].cDescr,
Expand Down Expand Up @@ -4908,6 +4909,14 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) {
"edge. "
"If not present, return 0. Note that some ice belts may in fact have a "
"southern edge at the equator.");

fvFormattedString(&output[OUT_LANDFRACGLOBAL].cName, "LandFracGlobal");
fvFormattedString(&output[OUT_LANDFRACGLOBAL].cDescr,
"Fraction of planetary surface covered by land");
output[OUT_LANDFRACGLOBAL].bNeg = 0;
output[OUT_LANDFRACGLOBAL].iNum = 1;
output[OUT_LANDFRACGLOBAL].iModuleBit = POISE;
fnWrite[OUT_LANDFRACGLOBAL] = &WriteLandFracGlobal;
}

/************ POISE Logging Functions **************/
Expand Down Expand Up @@ -8221,3 +8230,14 @@ void PoiseIceSheets(BODY *body, EVOLVE *evolve, int iBody) {
}
}
}

double fdLandFracGlobal(BODY *body, int iBody) {
int iLat;
double dLandFrac = 0;

for (iLat = 0; iLat < body[iBody].iNumLats; iLat++) {
dLandFrac += body[iBody].daLandFrac[iLat];
}
dLandFrac /= body[iBody].iNumLats;
return dLandFrac;
}
3 changes: 2 additions & 1 deletion src/poise.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -186,6 +186,7 @@
#define OUT_NORTHICEBELTLATSEA 1973
#define OUT_SOUTHICEBELTLATLAND 1974
#define OUT_SOUTHICEBELTLATSEA 1975
#define OUT_LANDFRACGLOBAL 1976

/* @cond DOXYGEN_OVERRIDE */

Expand Down Expand Up @@ -310,5 +311,5 @@ double fdPoiseDIceMassDtFlow(BODY *, SYSTEM *, int *);

double fdEccTrueAnomaly(double, double);
double fdAlbedoTOA350(double, double, double, double);

double fdLandFracGlobal(BODY *, int);
/* @endcond */
27 changes: 13 additions & 14 deletions src/vplanet.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -634,19 +634,17 @@ struct BODY {
double dAlbedoWater; /**< Sets base albedo of water (sea model) */
double dLatLandWater; /**< Lattitude boundary between land and water in polar
and equatorial options */
int iLatLandWater; /**< Lattitude boundary between land and water in polar and
equatorial options */
int bAlbedoZA; /**< Use albedo based on zenith angle (ann model) */
double dAreaIceCov; /**< Tracks area of surface covered in permanent ice*/
double dAstroDist; /**< Distance between primary and planet */
int bCalcAB; /**< Calc A and B from Williams & Kasting 1997 */
int iClimateModel; /**< Which EBM to be used (ann or sea) */
int bColdStart; /**< Start from global glaciation (snowball) conditions */
double dCw_dt; /**< Heat capacity of water / EBM time step */
double dDiffCoeff; /**< Diffusion coefficient set by user */
int bDiffRot; /**< Adjust heat diffusion for rotation rate */
int bElevFB; /**< Apply elevation feedback to ice ablation */
double dFixIceLat; /**< Fixes ice line latitude to user set value */
int bAlbedoZA; /**< Use albedo based on zenith angle (ann model) */
double dAreaIceCov; /**< Tracks area of surface covered in permanent ice*/
double dAstroDist; /**< Distance between primary and planet */
int bCalcAB; /**< Calc A and B from Williams & Kasting 1997 */
int iClimateModel; /**< Which EBM to be used (ann or sea) */
int bColdStart; /**< Start from global glaciation (snowball) conditions */
double dCw_dt; /**< Heat capacity of water / EBM time step */
double dDiffCoeff; /**< Diffusion coefficient set by user */
int bDiffRot; /**< Adjust heat diffusion for rotation rate */
int bElevFB; /**< Apply elevation feedback to ice ablation */
double dFixIceLat; /**< Fixes ice line latitude to user set value */
double dFluxInGlobal; /**< Global mean of incoming flux */
double dFluxInGlobalTmp; /**< Copy of global mean incoming flux */
double dFluxOutGlobal; /**< Global mean of outgoing flux */
Expand Down Expand Up @@ -724,7 +722,8 @@ struct BODY {
double *daFlux; /**< Meridional surface heat flux */
double *daFluxIn; /**< Incoming surface flux (insolation) */
double *daFluxOut; /**< Outgoing surface flux (longwave) */
double *daLats; /**< Latitude of each cell (centered); South Pole is 0 */
double *daLats; /**< Latitude of each cell (centered), assuming equal areas
per bin; South Pole is 0 */
double *daPeakInsol; /**< Annually averaged insolation at each latitude */
double *daTGrad; /**< Gradient of temperature (meridional) */

Expand Down

0 comments on commit 3bcdcbd

Please sign in to comment.