diff --git a/examples/ClimateLand/equatorial.in b/examples/ClimateLand/equatorial.in index bc57c4a3..73e38455 100755 --- a/examples/ClimateLand/equatorial.in +++ b/examples/ClimateLand/equatorial.in @@ -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) diff --git a/src/poise.c b/src/poise.c index 6fbd41de..05e21769 100644 --- a/src/poise.c +++ b/src/poise.c @@ -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, @@ -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++) { @@ -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++) { @@ -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 * @@ -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) { @@ -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 / @@ -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, @@ -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, @@ -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 **************/ @@ -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; +} diff --git a/src/poise.h b/src/poise.h index 3c0c1375..333a9487 100644 --- a/src/poise.h +++ b/src/poise.h @@ -186,6 +186,7 @@ #define OUT_NORTHICEBELTLATSEA 1973 #define OUT_SOUTHICEBELTLATLAND 1974 #define OUT_SOUTHICEBELTLATSEA 1975 +#define OUT_LANDFRACGLOBAL 1976 /* @cond DOXYGEN_OVERRIDE */ @@ -310,5 +311,5 @@ double fdPoiseDIceMassDtFlow(BODY *, SYSTEM *, int *); double fdEccTrueAnomaly(double, double); double fdAlbedoTOA350(double, double, double, double); - +double fdLandFracGlobal(BODY *, int); /* @endcond */ diff --git a/src/vplanet.h b/src/vplanet.h index fad4f4be..cb44d1cc 100644 --- a/src/vplanet.h +++ b/src/vplanet.h @@ -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 */ @@ -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) */