From a78a4e38c7d40dcab0b2b8ebb9b2ef0c40063422 Mon Sep 17 00:00:00 2001 From: Rory Barnes Date: Fri, 16 Aug 2024 16:59:07 -0700 Subject: [PATCH 01/12] Created new options for sGeography: random, polar, and equatorial. Created new example ClimateLand to show differences between these choices. Still testing. Random may be OK, but polar and equatorial appear to be broken. --- examples/ClimateLand/equatorial.in | 71 ++++++++++ examples/ClimateLand/modern.in | 70 ++++++++++ examples/ClimateLand/polar.in | 71 ++++++++++ examples/ClimateLand/random.in | 71 ++++++++++ examples/ClimateLand/sun.in | 9 ++ examples/ClimateLand/uniform.in | 71 ++++++++++ examples/ClimateLand/vpl.in | 20 +++ src/galhabit.c | 28 ---- src/galhabit.h | 1 - src/options.c | 27 ++++ src/options.h | 3 + src/poise.c | 212 ++++++++++++++++++++++++----- src/poise.h | 8 +- src/vplanet.c | 2 + src/vplanet.h | 2 + 15 files changed, 604 insertions(+), 62 deletions(-) create mode 100755 examples/ClimateLand/equatorial.in create mode 100755 examples/ClimateLand/modern.in create mode 100755 examples/ClimateLand/polar.in create mode 100755 examples/ClimateLand/random.in create mode 100755 examples/ClimateLand/sun.in create mode 100755 examples/ClimateLand/uniform.in create mode 100755 examples/ClimateLand/vpl.in diff --git a/examples/ClimateLand/equatorial.in b/examples/ClimateLand/equatorial.in new file mode 100755 index 000000000..e3a267408 --- /dev/null +++ b/examples/ClimateLand/equatorial.in @@ -0,0 +1,71 @@ +sName equatorial #name of planet +saModules poise #what vplanet modules you want to use +sGeography equitorial +dLandWaterLatitude 45 + +dMass 3.00316726e-06 #mass of planet +dRadius -1.00 #radius (not important right now) +dRotPeriod -1.00000 #rotation period (minus = days) +dObliquity 23.5 +dSemi 1.02 +dEcc 0.0 #eccentricity of orbit +dLongP 0 #pericenter, wrt Earth's position at spring equinox +# note that this is the typical value +180, +# since that one is solar position +dDynEllip 0.0 #shape of planet (0 = a sphere) +dPrecA 0.0 #orientation of spin axis + +#_______addition disorb/distrot parameters (leave these alone for now)__________________ +#dInc 5e-5 #inclination of orbit +#dLongA 348.73936 #orientation of orbital plane +#bGRCorr 0 #use GR correction (not important) +#bInvPlane 1 #convert to invariable plane coords +#bOverrideMaxEcc 1 #override max ecc halt (not recommended) +#dHaltMaxEcc 0.4 #eccentricity at which to halt simulation + +#_______poise parameters (have fun with these!)_________________________________________ +iLatCellNum 151 #number of latitude cells +sClimateModel sea #use seasonal or annual model +dTGlobalInit 14.85 #initial guess at average surface temp +iNumYears 4 #number of years (orbits) to run clim model +iNStepInYear 80 #number of steps to take in a "year" +#dSurfAlbedo 0.35 #average surface albedo (annual model only) + +#__ice params_________ +bIceSheets 1 #enable ice sheets +dInitIceLat 90. #how low do initial ice sheet extend? +dInitIceHeight 0. #height of initial ice sheets +dIceDepRate 2.25e-5 #rate of snow build up (when T < 0) +dIceAlbedo 0.6 #albedo of ice +iIceDt 1 #time step of ice-sheet model (orbits) +iReRunSeas 500 #how often to re-run seasonal model +bSeaIceModel 0 #use sea ice model (slow!) +bSkipSeasEnabled 0 #can skip seasonal if snowball state present + +#__heat diffusion______ +#bMEPDiff 1 #calculate diffusion using max entropy production +#bHadley 1 #mimic hadley heat diffusion +dDiffusion 0.58 #diffusion coefficient (fixed) +dNuLandWater 0.8 #Heat diffusion coefficient between Land and Water + +#__outgoing flux_______ +dPlanckA 203.3 #offset for OLR calculation (greenhouse) +dPlanckB 2.09 #slope of OLR calc (water vapor feedback) +bCalcAB 0 #calculate A & B from Kasting model fits +#dpCO2 0.00028 #partial pressure of co2 + +#__surface properties__ +dAlbedoLand 0.363 #albedo of land +dAlbedoWater 0.263 #albedo of water +dHeatCapLand 1.55e7 #land heat capacity +dHeatCapWater 4.428e6 #water heat capacity +dMixingDepth 70 #mixing depth of ocean + + +#________output options!_____________________________________________ +saOutputOrder Time PrecA -TGlobal AlbedoGlobal -FluxOutGlobal $ + -TotIceMass -TotIceFlow -TotIceBalance DeltaTime AreaIceCov Snowball Obliq $ + IceBeltLand IceBeltSea Ecce + +saGridOutput Time -Latitude -TempLat AlbedoLat -AnnInsol -FluxIn -FluxOut IceMass -IceHeight DIceMassDt $ + -IceFlow -BedrockH -TempMaxLat -TempMinLat -FluxMerid -DivFlux diff --git a/examples/ClimateLand/modern.in b/examples/ClimateLand/modern.in new file mode 100755 index 000000000..a8c7304bd --- /dev/null +++ b/examples/ClimateLand/modern.in @@ -0,0 +1,70 @@ +sName modern #name of planet +saModules poise #what vplanet modules you want to use +sGeography modern + +dMass 3.00316726e-06 #mass of planet +dRadius -1.00 #radius (not important right now) +dRotPeriod -1.00000 #rotation period (minus = days) +dObliquity 23.5 +dSemi 1.02 +dEcc 0.0 #eccentricity of orbit +dLongP 0 #pericenter, wrt Earth's position at spring equinox +# note that this is the typical value +180, +# since that one is solar position +dDynEllip 0.0 #shape of planet (0 = a sphere) +dPrecA 0.0 #orientation of spin axis + +#_______addition disorb/distrot parameters (leave these alone for now)__________________ +#dInc 5e-5 #inclination of orbit +#dLongA 348.73936 #orientation of orbital plane +#bGRCorr 0 #use GR correction (not important) +#bInvPlane 1 #convert to invariable plane coords +#bOverrideMaxEcc 1 #override max ecc halt (not recommended) +#dHaltMaxEcc 0.4 #eccentricity at which to halt simulation + +#_______poise parameters (have fun with these!)_________________________________________ +iLatCellNum 151 #number of latitude cells +sClimateModel sea #use seasonal or annual model +dTGlobalInit 14.85 #initial guess at average surface temp +iNumYears 4 #number of years (orbits) to run clim model +iNStepInYear 80 #number of steps to take in a "year" +#dSurfAlbedo 0.35 #average surface albedo (annual model only) + +#__ice params_________ +bIceSheets 1 #enable ice sheets +dInitIceLat 90. #how low do initial ice sheet extend? +dInitIceHeight 0. #height of initial ice sheets +dIceDepRate 2.25e-5 #rate of snow build up (when T < 0) +dIceAlbedo 0.6 #albedo of ice +iIceDt 1 #time step of ice-sheet model (orbits) +iReRunSeas 500 #how often to re-run seasonal model +bSeaIceModel 0 #use sea ice model (slow!) +bSkipSeasEnabled 0 #can skip seasonal if snowball state present + +#__heat diffusion______ +#bMEPDiff 1 #calculate diffusion using max entropy production +#bHadley 1 #mimic hadley heat diffusion +dDiffusion 0.58 #diffusion coefficient (fixed) +dNuLandWater 0.8 #Heat diffusion coefficient between Land and Water + +#__outgoing flux_______ +dPlanckA 203.3 #offset for OLR calculation (greenhouse) +dPlanckB 2.09 #slope of OLR calc (water vapor feedback) +bCalcAB 0 #calculate A & B from Kasting model fits +#dpCO2 0.00028 #partial pressure of co2 + +#__surface properties__ +dAlbedoLand 0.363 #albedo of land +dAlbedoWater 0.263 #albedo of water +dHeatCapLand 1.55e7 #land heat capacity +dHeatCapWater 4.428e6 #water heat capacity +dMixingDepth 70 #mixing depth of ocean + + +#________output options!_____________________________________________ +saOutputOrder Time PrecA -TGlobal AlbedoGlobal -FluxOutGlobal $ + -TotIceMass -TotIceFlow -TotIceBalance DeltaTime AreaIceCov Snowball Obliq $ + IceBeltLand IceBeltSea Ecce + +saGridOutput Time -Latitude -TempLat AlbedoLat -AnnInsol -FluxIn -FluxOut IceMass -IceHeight DIceMassDt $ + -IceFlow -BedrockH -TempMaxLat -TempMinLat -FluxMerid -DivFlux diff --git a/examples/ClimateLand/polar.in b/examples/ClimateLand/polar.in new file mode 100755 index 000000000..1a0d84d41 --- /dev/null +++ b/examples/ClimateLand/polar.in @@ -0,0 +1,71 @@ +sName polar #name of planet +saModules poise #what vplanet modules you want to use +sGeography polar +dLandWaterLatitude 45 + +dMass 3.00316726e-06 #mass of planet +dRadius -1.00 #radius (not important right now) +dRotPeriod -1.00000 #rotation period (minus = days) +dObliquity 23.5 +dSemi 1.02 +dEcc 0.0 #eccentricity of orbit +dLongP 0 #pericenter, wrt Earth's position at spring equinox +# note that this is the typical value +180, +# since that one is solar position +dDynEllip 0.0 #shape of planet (0 = a sphere) +dPrecA 0.0 #orientation of spin axis + +#_______addition disorb/distrot parameters (leave these alone for now)__________________ +#dInc 5e-5 #inclination of orbit +#dLongA 348.73936 #orientation of orbital plane +#bGRCorr 0 #use GR correction (not important) +#bInvPlane 1 #convert to invariable plane coords +#bOverrideMaxEcc 1 #override max ecc halt (not recommended) +#dHaltMaxEcc 0.4 #eccentricity at which to halt simulation + +#_______poise parameters (have fun with these!)_________________________________________ +iLatCellNum 151 #number of latitude cells +sClimateModel sea #use seasonal or annual model +dTGlobalInit 14.85 #initial guess at average surface temp +iNumYears 4 #number of years (orbits) to run clim model +iNStepInYear 80 #number of steps to take in a "year" +#dSurfAlbedo 0.35 #average surface albedo (annual model only) + +#__ice params_________ +bIceSheets 1 #enable ice sheets +dInitIceLat 90. #how low do initial ice sheet extend? +dInitIceHeight 0. #height of initial ice sheets +dIceDepRate 2.25e-5 #rate of snow build up (when T < 0) +dIceAlbedo 0.6 #albedo of ice +iIceDt 1 #time step of ice-sheet model (orbits) +iReRunSeas 500 #how often to re-run seasonal model +bSeaIceModel 0 #use sea ice model (slow!) +bSkipSeasEnabled 0 #can skip seasonal if snowball state present + +#__heat diffusion______ +#bMEPDiff 1 #calculate diffusion using max entropy production +#bHadley 1 #mimic hadley heat diffusion +dDiffusion 0.58 #diffusion coefficient (fixed) +dNuLandWater 0.8 #Heat diffusion coefficient between Land and Water + +#__outgoing flux_______ +dPlanckA 203.3 #offset for OLR calculation (greenhouse) +dPlanckB 2.09 #slope of OLR calc (water vapor feedback) +bCalcAB 0 #calculate A & B from Kasting model fits +#dpCO2 0.00028 #partial pressure of co2 + +#__surface properties__ +dAlbedoLand 0.363 #albedo of land +dAlbedoWater 0.263 #albedo of water +dHeatCapLand 1.55e7 #land heat capacity +dHeatCapWater 4.428e6 #water heat capacity +dMixingDepth 70 #mixing depth of ocean + + +#________output options!_____________________________________________ +saOutputOrder Time PrecA -TGlobal AlbedoGlobal -FluxOutGlobal $ + -TotIceMass -TotIceFlow -TotIceBalance DeltaTime AreaIceCov Snowball Obliq $ + IceBeltLand IceBeltSea Ecce + +saGridOutput Time -Latitude -TempLat AlbedoLat -AnnInsol -FluxIn -FluxOut IceMass -IceHeight DIceMassDt $ + -IceFlow -BedrockH -TempMaxLat -TempMinLat -FluxMerid -DivFlux diff --git a/examples/ClimateLand/random.in b/examples/ClimateLand/random.in new file mode 100755 index 000000000..05139b93c --- /dev/null +++ b/examples/ClimateLand/random.in @@ -0,0 +1,71 @@ +sName random #name of planet +saModules poise #what vplanet modules you want to use +sGeography random +iRandSeed 4.567e8 + +dMass 3.00316726e-06 #mass of planet +dRadius -1.00 #radius (not important right now) +dRotPeriod -1.00000 #rotation period (minus = days) +dObliquity 23.5 +dSemi 1.02 +dEcc 0.0 #eccentricity of orbit +dLongP 0 #pericenter, wrt Earth's position at spring equinox +# note that this is the typical value +180, +# since that one is solar position +dDynEllip 0.0 #shape of planet (0 = a sphere) +dPrecA 0.0 #orientation of spin axis + +#_______addition disorb/distrot parameters (leave these alone for now)__________________ +#dInc 5e-5 #inclination of orbit +#dLongA 348.73936 #orientation of orbital plane +#bGRCorr 0 #use GR correction (not important) +#bInvPlane 1 #convert to invariable plane coords +#bOverrideMaxEcc 1 #override max ecc halt (not recommended) +#dHaltMaxEcc 0.4 #eccentricity at which to halt simulation + +#_______poise parameters (have fun with these!)_________________________________________ +iLatCellNum 151 #number of latitude cells +sClimateModel sea #use seasonal or annual model +dTGlobalInit 14.85 #initial guess at average surface temp +iNumYears 4 #number of years (orbits) to run clim model +iNStepInYear 80 #number of steps to take in a "year" +#dSurfAlbedo 0.35 #average surface albedo (annual model only) + +#__ice params_________ +bIceSheets 1 #enable ice sheets +dInitIceLat 90. #how low do initial ice sheet extend? +dInitIceHeight 0. #height of initial ice sheets +dIceDepRate 2.25e-5 #rate of snow build up (when T < 0) +dIceAlbedo 0.6 #albedo of ice +iIceDt 1 #time step of ice-sheet model (orbits) +iReRunSeas 500 #how often to re-run seasonal model +bSeaIceModel 0 #use sea ice model (slow!) +bSkipSeasEnabled 0 #can skip seasonal if snowball state present + +#__heat diffusion______ +#bMEPDiff 1 #calculate diffusion using max entropy production +#bHadley 1 #mimic hadley heat diffusion +dDiffusion 0.58 #diffusion coefficient (fixed) +dNuLandWater 0.8 #Heat diffusion coefficient between Land and Water + +#__outgoing flux_______ +dPlanckA 203.3 #offset for OLR calculation (greenhouse) +dPlanckB 2.09 #slope of OLR calc (water vapor feedback) +bCalcAB 0 #calculate A & B from Kasting model fits +#dpCO2 0.00028 #partial pressure of co2 + +#__surface properties__ +dAlbedoLand 0.363 #albedo of land +dAlbedoWater 0.263 #albedo of water +dHeatCapLand 1.55e7 #land heat capacity +dHeatCapWater 4.428e6 #water heat capacity +dMixingDepth 70 #mixing depth of ocean + + +#________output options!_____________________________________________ +saOutputOrder Time PrecA -TGlobal AlbedoGlobal -FluxOutGlobal $ + -TotIceMass -TotIceFlow -TotIceBalance DeltaTime AreaIceCov Snowball Obliq $ + IceBeltLand IceBeltSea Ecce + +saGridOutput Time -Latitude -TempLat AlbedoLat -AnnInsol -FluxIn -FluxOut IceMass -IceHeight DIceMassDt $ + -IceFlow -BedrockH -TempMaxLat -TempMinLat -FluxMerid -DivFlux diff --git a/examples/ClimateLand/sun.in b/examples/ClimateLand/sun.in new file mode 100755 index 000000000..a1f97c666 --- /dev/null +++ b/examples/ClimateLand/sun.in @@ -0,0 +1,9 @@ +# sun parameters +sName sun +dMass 1 +dSemi 0 +dEcc 0 +dRadius 0.00135 +dLuminosity -1 +sStellarModel none #sun does not change over time +saModules stellar #use stellar module (needed for luminosity) diff --git a/examples/ClimateLand/uniform.in b/examples/ClimateLand/uniform.in new file mode 100755 index 000000000..15c83aa2a --- /dev/null +++ b/examples/ClimateLand/uniform.in @@ -0,0 +1,71 @@ +sName uniform #name of planet +saModules poise #what vplanet modules you want to use +sGeography uniform +dLandFrac 0.5 + +dMass 3.00316726e-06 #mass of planet +dRadius -1.00 #radius (not important right now) +dRotPeriod -1.00000 #rotation period (minus = days) +dObliquity 23.5 +dSemi 1.02 +dEcc 0.0 #eccentricity of orbit +dLongP 0 #pericenter, wrt Earth's position at spring equinox +# note that this is the typical value +180, +# since that one is solar position +dDynEllip 0.0 #shape of planet (0 = a sphere) +dPrecA 0.0 #orientation of spin axis + +#_______addition disorb/distrot parameters (leave these alone for now)__________________ +#dInc 5e-5 #inclination of orbit +#dLongA 348.73936 #orientation of orbital plane +#bGRCorr 0 #use GR correction (not important) +#bInvPlane 1 #convert to invariable plane coords +#bOverrideMaxEcc 1 #override max ecc halt (not recommended) +#dHaltMaxEcc 0.4 #eccentricity at which to halt simulation + +#_______poise parameters (have fun with these!)_________________________________________ +iLatCellNum 151 #number of latitude cells +sClimateModel sea #use seasonal or annual model +dTGlobalInit 14.85 #initial guess at average surface temp +iNumYears 4 #number of years (orbits) to run clim model +iNStepInYear 80 #number of steps to take in a "year" +#dSurfAlbedo 0.35 #average surface albedo (annual model only) + +#__ice params_________ +bIceSheets 1 #enable ice sheets +dInitIceLat 90. #how low do initial ice sheet extend? +dInitIceHeight 0. #height of initial ice sheets +dIceDepRate 2.25e-5 #rate of snow build up (when T < 0) +dIceAlbedo 0.6 #albedo of ice +iIceDt 1 #time step of ice-sheet model (orbits) +iReRunSeas 500 #how often to re-run seasonal model +bSeaIceModel 0 #use sea ice model (slow!) +bSkipSeasEnabled 0 #can skip seasonal if snowball state present + +#__heat diffusion______ +#bMEPDiff 1 #calculate diffusion using max entropy production +#bHadley 1 #mimic hadley heat diffusion +dDiffusion 0.58 #diffusion coefficient (fixed) +dNuLandWater 0.8 #Heat diffusion coefficient between Land and Water + +#__outgoing flux_______ +dPlanckA 203.3 #offset for OLR calculation (greenhouse) +dPlanckB 2.09 #slope of OLR calc (water vapor feedback) +bCalcAB 0 #calculate A & B from Kasting model fits +#dpCO2 0.00028 #partial pressure of co2 + +#__surface properties__ +dAlbedoLand 0.363 #albedo of land +dAlbedoWater 0.263 #albedo of water +dHeatCapLand 1.55e7 #land heat capacity +dHeatCapWater 4.428e6 #water heat capacity +dMixingDepth 70 #mixing depth of ocean + + +#________output options!_____________________________________________ +saOutputOrder Time PrecA -TGlobal AlbedoGlobal -FluxOutGlobal $ + -TotIceMass -TotIceFlow -TotIceBalance DeltaTime AreaIceCov Snowball Obliq $ + IceBeltLand IceBeltSea Ecce + +saGridOutput Time -Latitude -TempLat AlbedoLat -AnnInsol -FluxIn -FluxOut IceMass -IceHeight DIceMassDt $ + -IceFlow -BedrockH -TempMaxLat -TempMinLat -FluxMerid -DivFlux diff --git a/examples/ClimateLand/vpl.in b/examples/ClimateLand/vpl.in new file mode 100755 index 000000000..50acb1c16 --- /dev/null +++ b/examples/ClimateLand/vpl.in @@ -0,0 +1,20 @@ +sSystemName climateland +iVerbose 5 #how much do you want vplanet to yell at you? +iDigits 6 #how many digits do you want in your numbers? +bOverwrite 1 #overwrite old files +sUnitMass solar #mass unit used for input +sUnitLength au #length unit used for input +sUnitTime y #time unit +sUnitAngle d #angle unit +bDoLog 1 #create log file +saBodyFiles sun.in $ #you must list all input files here (except vpl.in) + uniform.in $ # Uniform land distribution across all latitudes + modern.in $ # Modern Earth's land distribution + random.in $ # Random land distribution + polar.in $ # Land masses at the poles + equatorial.in # One land mass along the equator +bDoForward 1 #integrate forward in time +bVarDt 1 #use variable time stepping (not relevant to poise) +dEta 0.1 #how much to scale variable time step +dStopTime 1e6 #how long should the integration be +dOutputTime 1e3 #how much output you want diff --git a/src/galhabit.c b/src/galhabit.c index b85545188..bee152c8e 100644 --- a/src/galhabit.c +++ b/src/galhabit.c @@ -210,25 +210,6 @@ void ReadDMDensity(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, } } -void ReadRandSeed(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, - SYSTEM *system, int iFile) { - /* This parameter can exist in any file, but only once */ - int lTmp = -1; - int iTmp; - - AddOptionInt(files->Infile[iFile].cIn, options->cName, &iTmp, &lTmp, - control->Io.iVerbose); - if (lTmp >= 0) { - CheckDuplication(files, options, files->Infile[iFile].cIn, lTmp, - control->Io.iVerbose); - system->iSeed = iTmp; - UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); - } else { - AssignDefaultInt(options, &system->iSeed, files->iNumInputs); - } -} - - void ReadEncounterRad(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, SYSTEM *system, int iFile) { /* This parameter can exist in any file, but only once */ @@ -583,15 +564,6 @@ void InitializeOptionsGalHabit(OPTIONS *options, fnReadOption fnRead[]) { options[OPT_GALACDENSITY].bMultiFile = 0; fnRead[OPT_GALACDENSITY] = &ReadGalacDensity; - fvFormattedString(&options[OPT_RANDSEED].cName, "iRandSeed"); - fvFormattedString(&options[OPT_RANDSEED].cDescr, - "Seed for random number generator (stellar encounters)"); - fvFormattedString(&options[OPT_RANDSEED].cDefault, "42"); - options[OPT_RANDSEED].dDefault = 42; - options[OPT_RANDSEED].iType = 1; - options[OPT_RANDSEED].bMultiFile = 0; - fnRead[OPT_RANDSEED] = &ReadRandSeed; - fvFormattedString(&options[OPT_ENCOUNTERRAD].cName, "dEncounterRad"); fvFormattedString(&options[OPT_ENCOUNTERRAD].cDescr, "Radius at which stellar encounters occur"); diff --git a/src/galhabit.h b/src/galhabit.h index 8a49e326a..17658c721 100644 --- a/src/galhabit.h +++ b/src/galhabit.h @@ -15,7 +15,6 @@ #define OPTENDGALHABIT 2300 /* End of GALHABIT options */ #define OPT_GALACDENSITY 2201 -#define OPT_RANDSEED 2202 #define OPT_ENCOUNTERRAD 2203 #define OPT_RFORM 2204 #define OPT_TMIGRATION 2205 diff --git a/src/options.c b/src/options.c index 17b3409f8..c8d82c6d5 100644 --- a/src/options.c +++ b/src/options.c @@ -3151,6 +3151,24 @@ void ReadRadiusGyration(BODY *body, CONTROL *control, FILES *files, } } +void ReadRandSeed(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, + SYSTEM *system, int iFile) { + /* This parameter can exist in any file, but only once */ + int lTmp = -1; + int iTmp; + + AddOptionInt(files->Infile[iFile].cIn, options->cName, &iTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + CheckDuplication(files, options, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + system->iSeed = iTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else { + AssignDefaultInt(options, &system->iSeed, files->iNumInputs); + } +} + /* Rotation Period */ void ReadRotPeriod(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, @@ -4416,6 +4434,15 @@ void InitializeOptionsGeneral(OPTIONS *options, fnReadOption fnRead[]) { options[OPT_RG].iFileType = 1; fnRead[OPT_RG] = &ReadRadiusGyration; + fvFormattedString(&options[OPT_RANDSEED].cName, "iRandSeed"); + fvFormattedString(&options[OPT_RANDSEED].cDescr, + "Seed for random number generator (stellar encounters)"); + fvFormattedString(&options[OPT_RANDSEED].cDefault, "42"); + options[OPT_RANDSEED].dDefault = 42; + options[OPT_RANDSEED].iType = 1; + options[OPT_RANDSEED].bMultiFile = 0; + fnRead[OPT_RANDSEED] = &ReadRandSeed; + fvFormattedString(&options[OPT_ROTPER].cName, "dRotPeriod"); fvFormattedString(&options[OPT_ROTPER].cDescr, "Rotation Period"); fvFormattedString(&options[OPT_ROTPER].cDefault, "1 Day"); diff --git a/src/options.h b/src/options.h index 18f67ffd4..2459f3f9d 100644 --- a/src/options.h +++ b/src/options.h @@ -117,6 +117,9 @@ #define OPT_MINENVELOPEMASS \ 817 /**< Minimum envelope mass (evaporated below this) */ +#define OPT_RANDSEED 850 + + /* @cond DOXYGEN_OVERRIDE */ void GetWords(char cLine[], char[MAXARRAY][OPTLEN], int *, int *); diff --git a/src/poise.c b/src/poise.c index 612a4b998..b80629e97 100644 --- a/src/poise.c +++ b/src/poise.c @@ -606,15 +606,21 @@ void ReadGeography(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, control->Io.iVerbose); //CheckDuplication(files,options,files->Infile[iFile].cIn,lTmp,\ control->Io.iVerbose); - if (!memcmp(sLower(cTmp), "uni3", 4)) { - body[iFile - 1].iGeography = UNIFORM3; - } else if (!memcmp(sLower(cTmp), "modn", 4)) { - body[iFile - 1].iGeography = MODERN; + if (!memcmp(sLower(cTmp), "u", 1)) { + body[iFile - 1].iGeography = LANDWATERUNIFORM; + } else if (!memcmp(sLower(cTmp), "m", 1)) { + body[iFile - 1].iGeography = LANDWATERMODERN; + } else if (!memcmp(sLower(cTmp), "r", 1)) { + body[iFile - 1].iGeography = LANDWATERRANDOM; + } else if (!memcmp(sLower(cTmp), "p", 1)) { + body[iFile - 1].iGeography = LANDWATERPOLAR; + } else if (!memcmp(sLower(cTmp), "e", 1)) { + body[iFile - 1].iGeography = LANDWATEREQUATORIAL; } else { if (control->Io.iVerbose >= VERBERR) { fprintf(stderr, "ERROR: Unknown argument to %s: %s." - " Options are uni3 or modn.\n", + " Options are uniform, modern, random, polar, or equatorial.\n", options->cName, cTmp); } LineExit(files->Infile[iFile].cIn, lTmp); @@ -625,6 +631,25 @@ void ReadGeography(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, } } +void ReadLandWaterLatitude(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, + SYSTEM *system, int iFile) { + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + body[iFile - 1].dLatLandWater = dTmp * fdUnitsAngle(control->Units[iFile-1].iAngle); + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else { + if (iFile > 0) { + body[iFile - 1].dLatLandWater = options->dDefault; + } + } +} + void ReadInitIceLat(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, SYSTEM *system, int iFile) { @@ -1475,11 +1500,33 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fvFormattedString(&options[OPT_GEOGRAPHY].cName, "sGeography"); fvFormattedString(&options[OPT_GEOGRAPHY].cDescr, "Type of land distribution"); - fvFormattedString(&options[OPT_GEOGRAPHY].cDefault, "uni3"); - options[OPT_GEOGRAPHY].dDefault = UNIFORM3; + fvFormattedString(&options[OPT_GEOGRAPHY].cDefault, "uniform"); + options[OPT_GEOGRAPHY].dDefault = LANDWATERUNIFORM; options[OPT_GEOGRAPHY].iType = 3; options[OPT_GEOGRAPHY].bMultiFile = 1; fnRead[OPT_GEOGRAPHY] = &ReadGeography; + fvFormattedString(&options[OPT_LANDWATERLATITUDE].cLongDescr, + "The distribution of land and water across a planetary surface. \"Uniform\" sets " + "all latitudes to have the land fraction set by %s. \"Modern\" is modern Earth's " + "distribution. \"Random\" sets a random ratio for each latitude. \"Polar\" creates " + "land masses that extend %s degrees (radians) from the pole. \"Equatorial\" creates " + "land masses that extend %s degress (radians) from the equator.",options[OPT_LANDFRAC].cName, + options[OPT_LANDWATERLATITUDE].cName,options[OPT_LANDWATERLATITUDE].cName + ); + + fvFormattedString(&options[OPT_LANDWATERLATITUDE].cName, "dLandWaterLatitude"); + fvFormattedString(&options[OPT_LANDWATERLATITUDE].cDescr, "N+S Latitude that separates land mass and oceans"); + fvFormattedString(&options[OPT_LANDWATERLATITUDE].cDefault, "+/-45 degrees"); + fvFormattedString(&options[OPT_LANDWATERLATITUDE].cDimension, "angle"); + options[OPT_LANDWATERLATITUDE].dDefault = PI/4; + options[OPT_LANDWATERLATITUDE].iType = 2; + options[OPT_LANDWATERLATITUDE].bMultiFile = 1; + fnRead[OPT_LANDWATERLATITUDE] = &ReadLandWaterLatitude; + fvFormattedString(&options[OPT_LANDWATERLATITUDE].cLongDescr, + "If the polar or equatorial geographies is selected (%s), then this option " + "set the boundary between the land mass and oceans. The distributions are " + "symmetric about the equator.",options[OPT_GEOGRAPHY].cName + ); fvFormattedString(&options[OPT_SEASOUTPUTTIME].cName, "dSeasOutputTime"); fvFormattedString(&options[OPT_SEASOUTPUTTIME].cDescr, "Output interval for seasonal" @@ -1913,39 +1960,99 @@ void InitializeLatGrid(BODY *body, int iBody) { } } +void fvInitializeLandWaterUniform(BODY *body,int iBody) { + int iLat; + + for (iLat = 0; iLat < body[iBody].iNumLats; iLat++) { + body[iBody].daLandFrac[iLat] = body[iBody].dLandFrac; + body[iBody].daWaterFrac[iLat] = 1.0 - body[iBody].daLandFrac[iLat]; + } +} + +void fvInitializeLandWaterModern(BODY *body,int iBody) { + int iLat; + + for (iLat = 0; iLat < body[iBody].iNumLats; iLat++) { + if (body[iBody].daLats[iLat] * 180. / PI <= -60) { + body[iBody].daLandFrac[iLat] = 0.95 / 1.0094; + } else if (body[iBody].daLats[iLat] * 180. / PI > -60 && + body[iBody].daLats[iLat] * 180. / PI <= -40) { + body[iBody].daLandFrac[iLat] = 0.05 / 1.0094; + } else if (body[iBody].daLats[iLat] * 180. / PI > -40 && + body[iBody].daLats[iLat] * 180. / PI <= 20) { + body[iBody].daLandFrac[iLat] = 0.25 / 1.0094; + } else if (body[iBody].daLats[iLat] * 180. / PI > 20 && + body[iBody].daLats[iLat] * 180. / PI <= 70) { + body[iBody].daLandFrac[iLat] = 0.5 / 1.0094; + } else { + body[iBody].daLandFrac[iLat] = 0.38 / 1.0094; + } + body[iBody].daWaterFrac[iLat] = 1.0 - body[iBody].daLandFrac[iLat]; + } +} + +void fvInitializeLandWaterRandom(BODY *body,int iBody) { + int iLat; + + for (iLat = 0; iLat < body[iBody].iNumLats; iLat++) { + body[iBody].daLandFrac[iLat] = (double)rand()/(double)RAND_MAX; + body[iBody].daWaterFrac[iLat] = 1.0 - body[iBody].daLandFrac[iLat]; + } +} + +void fvInitializeLandWaterPolar(BODY *body,int iBody) { + int iLat; + + for (iLat=0;iLat -60 && - body[iBody].daLats[iLat] * 180. / PI <= -40) { - body[iBody].daLandFrac[iLat] = 0.05 / 1.0094; - } else if (body[iBody].daLats[iLat] * 180. / PI > -40 && - body[iBody].daLats[iLat] * 180. / PI <= 20) { - body[iBody].daLandFrac[iLat] = 0.25 / 1.0094; - } else if (body[iBody].daLats[iLat] * 180. / PI > 20 && - body[iBody].daLats[iLat] * 180. / PI <= 70) { - body[iBody].daLandFrac[iLat] = 0.5 / 1.0094; - } else { - body[iBody].daLandFrac[iLat] = 0.38 / 1.0094; - } - body[iBody].daWaterFrac[iLat] = 1.0 - body[iBody].daLandFrac[iLat]; - } + if (body[iBody].iGeography == LANDWATERUNIFORM) { + fvInitializeLandWaterUniform(body,iBody); + } else if (body[iBody].iGeography == LANDWATERMODERN) { + fvInitializeLandWaterModern(body,iBody); + } else if (body[iBody].iGeography == LANDWATERRANDOM) { + fvInitializeLandWaterRandom(body,iBody); + } else if (body[iBody].iGeography == LANDWATERPOLAR) { + fvInitializeLandWaterPolar(body,iBody); + } else if (body[iBody].iGeography == LANDWATEREQUATORIAL) { + fvInitializeLandWaterEquatorial(body,iBody); } } - void DampTemp(BODY *body, double dTGlobalTmp, int iBody) { int iLat; double deltaT; @@ -2624,6 +2731,47 @@ void NullPoiseDerivatives(BODY *body, EVOLVE *evolve, UPDATE *update, // Nothing here, because entire climate simulation is ran in ForceBehavior } +void VerifyGeography(BODY *body,CONTROL *control,OPTIONS *options,SYSTEM *system,int iBody) { + + if (body[iBody].iGeography == LANDWATERUNIFORM || body[iBody].iGeography == LANDWATERMODERN || + body[iBody].iGeography == LANDWATERRANDOM) { + if (options[OPT_LANDWATERLATITUDE].iLine[iBody+1] != -1) { + if (control->Io.iVerbose > VERBINPUT) { + fprintf(stderr,"ERROR: Cannot set %s with \"uniform\", \"modern\", or \"random\" geography. " + "Note that \"uniform\" is the default.\n", + options[OPT_LANDWATERLATITUDE].cName); + } + DoubleLineExit(options[OPT_LANDWATERLATITUDE].cFile[iBody + 1],options[OPT_GEOGRAPHY].cFile[iBody + 1], + options[OPT_LANDWATERLATITUDE].iLine[iBody+1],options[OPT_GEOGRAPHY].iLine[iBody+1]); + } + } + + if (body[iBody].iGeography == LANDWATERMODERN || body[iBody].iGeography == LANDWATERRANDOM || + body[iBody].iGeography == LANDWATERPOLAR || body[iBody].iGeography == LANDWATEREQUATORIAL) { + if (options[OPT_LANDFRAC].iLine[iBody+1] != -1) { + if (control->Io.iVerbose > VERBINPUT) { + fprintf(stderr,"ERROR: %s can only be set when %s is set to \"uniform\". Note that \"uniform\" is the default.\n", + options[OPT_LANDFRAC].cName,options[OPT_GEOGRAPHY].cName); + } + } + } + + if (body[iBody].iGeography == LANDWATERRANDOM && options[OPT_RANDSEED].iLine[iBody+1] == -1) { + if (control->Io.iVerbose > VERBINPUT) { + fprintf(stderr,"\nWARNING: %s is set to \"radnom\", but %s is not set.\n\n", + options[OPT_GEOGRAPHY].cName,options[OPT_RANDSEED].cName); + } + } + + if (body[iBody].iGeography == LANDWATERRANDOM) { + srand(system->iSeed); + } + + if (body[iBody].iGeography == LANDWATERPOLAR || body[iBody].iGeography == LANDWATEREQUATORIAL) { + body[iBody].iLatLandWater = (int)(body[iBody].dLatLandWater*body[iBody].iNumLats/PI); + } +} + void VerifyPoise(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, OUTPUT *output, SYSTEM *system, UPDATE *update, int iBody, int iModule) { @@ -2642,6 +2790,8 @@ void VerifyPoise(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, control->Io.iVerbose); } + VerifyGeography(body,control,options,system,iBody); + /* Initialize climate arrays */ InitializeLatGrid(body, iBody); InitializeClimateParams(body, iBody, control->Io.iVerbose); diff --git a/src/poise.h b/src/poise.h index 2e028c21e..ee76140d1 100644 --- a/src/poise.h +++ b/src/poise.h @@ -26,8 +26,11 @@ #define ALBTAYLOR 1 /* Land Geography */ -#define UNIFORM3 0 -#define MODERN 1 +#define LANDWATERUNIFORM 0 +#define LANDWATERMODERN 1 +#define LANDWATERRANDOM 2 +#define LANDWATERPOLAR 3 +#define LANDWATEREQUATORIAL 4 // Constants for the ice model #define LFICE 3.34e5 // ??? @@ -111,6 +114,7 @@ #define OPT_ECCAMP 1968 #define OPT_ECCPER 1969 #define OPT_MINICEHEIGHT 1970 +#define OPT_LANDWATERLATITUDE 1980 #define OPT_OLRMODEL 1998 #define OPT_CLIMATEMODEL 1999 diff --git a/src/vplanet.c b/src/vplanet.c index 3c6709c46..cdd27940c 100644 --- a/src/vplanet.c +++ b/src/vplanet.c @@ -29,9 +29,11 @@ We need this wrapper so we can call `main_impl` from Python. int main_impl(int argc, char *argv[]) { #ifdef DEBUG #ifdef __x86_64__ + /* _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_INVALID); _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_OVERFLOW); _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_DIV_ZERO); +*/ //_MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_UNDERFLOW); fprintf(stderr, "INFO: Floating point trapping enabled.\n"); #else diff --git a/src/vplanet.h b/src/vplanet.h index 309341c5b..0dc1badeb 100644 --- a/src/vplanet.h +++ b/src/vplanet.h @@ -630,6 +630,8 @@ struct BODY { needed) */ double dAlbedoLand; /**< Sets base albedo of land (sea model) */ 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 */ From 0404fbb08a716d491ba80a973e47bb294f972253 Mon Sep 17 00:00:00 2001 From: Rory Barnes Date: Fri, 16 Aug 2024 17:05:51 -0700 Subject: [PATCH 02/12] Changed land/water fractions to never be 1/0, but instead 0.99/0.01 and the code appears to be working. --- src/poise.c | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/poise.c b/src/poise.c index b80629e97..be4a48476 100644 --- a/src/poise.c +++ b/src/poise.c @@ -2004,16 +2004,16 @@ void fvInitializeLandWaterPolar(BODY *body,int iBody) { int iLat; for (iLat=0;iLat Date: Fri, 6 Sep 2024 15:13:19 -0700 Subject: [PATCH 03/12] Working on making the random mode of land distribution more flexible --- src/poise.c | 102 +++++++++++++++++++++++++++++++++++++++++++------- src/poise.h | 10 +++-- src/vplanet.h | 4 +- 3 files changed, 99 insertions(+), 17 deletions(-) diff --git a/src/poise.c b/src/poise.c index be4a48476..d78c393d9 100644 --- a/src/poise.c +++ b/src/poise.c @@ -958,6 +958,56 @@ void ReadLandFrac(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, } } +void ReadLandFracMean(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, + SYSTEM *system, int iFile) { + /* This parameter cannot exist in primary file */ + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + if (dTmp > 1 || dTmp < 0) { + fprintf(stderr, "ERROR: %s must be in the range [0,1].\n", + options->cName); + LineExit(files->Infile[iFile].cIn, lTmp); + } + body[iFile - 1].dLandFracMean = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else { + if (iFile > 0) { + body[iFile - 1].dLandFracMean = options->dDefault; + } + } +} + +void ReadLandFracAmp(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, + SYSTEM *system, int iFile) { + /* This parameter cannot exist in primary file */ + int lTmp = -1; + double dTmp; + + AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp, + control->Io.iVerbose); + if (lTmp >= 0) { + NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, + control->Io.iVerbose); + if (dTmp > 1 || dTmp < 0) { + fprintf(stderr, "ERROR: %s must be in the range [0,1].\n", + options->cName); + LineExit(files->Infile[iFile].cIn, lTmp); + } + body[iFile - 1].dLandFracAmp = dTmp; + UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile); + } else { + if (iFile > 0) { + body[iFile - 1].dLandFracAmp = options->dDefault; + } + } +} + void ReadNuLandWater(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, SYSTEM *system, int iFile) { /* This parameter cannot exist in primary file */ @@ -1425,6 +1475,24 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { options[OPT_LANDFRAC].bMultiFile = 1; fnRead[OPT_LANDFRAC] = &ReadLandFrac; + fvFormattedString(&options[OPT_LANDFRACMEAN].cName, "dLandFracMean"); + fvFormattedString(&options[OPT_LANDFRACMEAN].cDescr, "Average fraction of land per latitude"); + fvFormattedString(&options[OPT_LANDFRACMEAN].cDefault, "0.5"); + fvFormattedString(&options[OPT_LANDFRACMEAN].cDimension, "nd"); + options[OPT_LANDFRACMEAN].dDefault = 0.5; + options[OPT_LANDFRACMEAN].iType = 2; + options[OPT_LANDFRACMEAN].bMultiFile = 1; + fnRead[OPT_LANDFRACMEAN] = &ReadLandFracMean; + + fvFormattedString(&options[OPT_LANDFRACAMP].cName, "dLandFrac"); + fvFormattedString(&options[OPT_LANDFRACAMP].cDescr, "Max difference in land fraction from the mean"); + fvFormattedString(&options[OPT_LANDFRACAMP].cDefault, "0.34"); + fvFormattedString(&options[OPT_LANDFRACAMP].cDimension, "nd"); + options[OPT_LANDFRACAMP].dDefault = 0.34; + options[OPT_LANDFRACAMP].iType = 2; + options[OPT_LANDFRACAMP].bMultiFile = 1; + fnRead[OPT_LANDFRACAMP] = &ReadLandFracAmp; + fvFormattedString(&options[OPT_NSTEPINYEAR].cName, "iNStepInYear"); fvFormattedString(&options[OPT_NSTEPINYEAR].cDescr, "Number of time-steps/year in" " seasonal model"); @@ -1993,9 +2061,17 @@ void fvInitializeLandWaterModern(BODY *body,int iBody) { void fvInitializeLandWaterRandom(BODY *body,int iBody) { int iLat; + double dOffset; for (iLat = 0; iLat < body[iBody].iNumLats; iLat++) { - body[iBody].daLandFrac[iLat] = (double)rand()/(double)RAND_MAX; + dOffset = body[iBody].dLandFracAmp * (1 - 2*(double)rand()/(double)RAND_MAX); + body[iBody].daLandFrac[iLat] = body[iBody].dLandFracMean + dOffset; + if (body[iBody].daLandFrac[iLat] > LANDFRACMAX) { + body[iBody].daLandFrac[iLat] = LANDFRACMAX; + } + if (body[iBody].daLandFrac[iLat] < LANDFRACMIN) { + body[iBody].daLandFrac[iLat] = LANDFRACMIN; + } body[iBody].daWaterFrac[iLat] = 1.0 - body[iBody].daLandFrac[iLat]; } } @@ -2004,16 +2080,16 @@ void fvInitializeLandWaterPolar(BODY *body,int iBody) { int iLat; for (iLat=0;iLat=263K (Pa^-3 s^-1 - ref?) #define Q1ICE 6e4 // energy in ice deformation at T<263K (J/mol) #define Q2ICE 13.9e4 // energy in ice deformation at T>=263 (J/mol) +#define LANDFRACMAX 0.99 +#define LANDFRACMIN 0.01 // Constant for the lithospheric model #define nGLEN 3.0 // Glen's law coefficient @@ -85,9 +87,11 @@ #define OPT_FILEORBITOBLDATA 1925 #define OPT_LANDFRAC 1940 -#define OPT_HEATCAPLAND 1942 -#define OPT_HEATCAPWATER 1943 -#define OPT_FRZTSEAICE 1944 +#define OPT_LANDFRACMEAN 1941 +#define OPT_LANDFRACAMP 1942 +#define OPT_HEATCAPLAND 1943 +#define OPT_HEATCAPWATER 1944 +#define OPT_FRZTSEAICE 1945 //#define OPT_LATENTHEAT 1945 #define OPT_ICECONDUCT 1946 #define OPT_MIXINGDEPTH 1947 diff --git a/src/vplanet.h b/src/vplanet.h index 0dc1badeb..ea8e282af 100644 --- a/src/vplanet.h +++ b/src/vplanet.h @@ -800,7 +800,9 @@ struct BODY { double **daInvMSea; /**< Inverted matrix in seasonal EBM */ double *daLambdaSea; /**< Diffusion terms in seasonal EBM matrix */ double dLandFrac; /**< Land fraction input by user */ - double *daLandFrac; /**< Fraction of cell which is land */ + double dLandFracMean; /**< Average land fraction per latitude */ + double dLandFracAmp; /**< Land fraction input by user */ + double *daLandFrac; /**< Max land fraction difference from the mean */ double **daMDiffSea; /**< Diffusion only matrix in seasonal EBM */ double **daMEulerCopySea; /**< Temporary copy of Euler time step matrix (seasonal) */ From a52db31c139bbd28a3c4363f16b8d103f9a5c1ea Mon Sep 17 00:00:00 2001 From: Rory Barnes Date: Fri, 6 Sep 2024 15:20:20 -0700 Subject: [PATCH 04/12] Added new file .geography to write out land/water fractions as a function of latitude. --- src/poise.c | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/src/poise.c b/src/poise.c index d78c393d9..cd6094569 100644 --- a/src/poise.c +++ b/src/poise.c @@ -2110,6 +2110,23 @@ void fvInitializeLandWaterEquatorial(BODY *body,int iBody) { } } +void fvWriteLandFracFile(BODY *body,int iBody) { + int iLat; + double dInterval,dLatNow; + FILE *fp; + char *cOut; + + dInterval = 180/body[iBody].iNumLats; + fvFormattedString(&cOut, "%s.geography", body[iBody].cName); + fp = open(cOut,"w"); + for (iLat=0;iLat Date: Fri, 6 Sep 2024 16:56:25 -0700 Subject: [PATCH 05/12] Modified examples/ClimateLand so all planets at 1 AU, and all have ~29.2% of their surface area covered in land. --- .gitignore | 1 + examples/ClimateLand/equatorial.in | 4 ++-- examples/ClimateLand/modern.in | 2 +- examples/ClimateLand/polar.in | 4 ++-- examples/ClimateLand/random.in | 4 +++- examples/ClimateLand/uniform.in | 4 ++-- examples/ClimateLand/vpl.in | 6 +++--- src/poise.c | 18 ++++++++---------- 8 files changed, 22 insertions(+), 21 deletions(-) diff --git a/.gitignore b/.gitignore index 500e18555..c0d4e42f6 100644 --- a/.gitignore +++ b/.gitignore @@ -18,6 +18,7 @@ SeasonalClimateFiles *.hdf5 *.vscode +*.geography #files created from running tests *MP_Checkpoint* diff --git a/examples/ClimateLand/equatorial.in b/examples/ClimateLand/equatorial.in index e3a267408..bc57c4a39 100755 --- a/examples/ClimateLand/equatorial.in +++ b/examples/ClimateLand/equatorial.in @@ -1,13 +1,13 @@ sName equatorial #name of planet saModules poise #what vplanet modules you want to use sGeography equitorial -dLandWaterLatitude 45 +dLandWaterLatitude 73.6 dMass 3.00316726e-06 #mass of planet dRadius -1.00 #radius (not important right now) dRotPeriod -1.00000 #rotation period (minus = days) dObliquity 23.5 -dSemi 1.02 +dSemi 1 dEcc 0.0 #eccentricity of orbit dLongP 0 #pericenter, wrt Earth's position at spring equinox # note that this is the typical value +180, diff --git a/examples/ClimateLand/modern.in b/examples/ClimateLand/modern.in index a8c7304bd..b93eeb957 100755 --- a/examples/ClimateLand/modern.in +++ b/examples/ClimateLand/modern.in @@ -6,7 +6,7 @@ dMass 3.00316726e-06 #mass of planet dRadius -1.00 #radius (not important right now) dRotPeriod -1.00000 #rotation period (minus = days) dObliquity 23.5 -dSemi 1.02 +dSemi 1 dEcc 0.0 #eccentricity of orbit dLongP 0 #pericenter, wrt Earth's position at spring equinox # note that this is the typical value +180, diff --git a/examples/ClimateLand/polar.in b/examples/ClimateLand/polar.in index 1a0d84d41..02c19eead 100755 --- a/examples/ClimateLand/polar.in +++ b/examples/ClimateLand/polar.in @@ -1,13 +1,13 @@ sName polar #name of planet saModules poise #what vplanet modules you want to use sGeography polar -dLandWaterLatitude 45 +dLandWaterLatitude 44.9 dMass 3.00316726e-06 #mass of planet dRadius -1.00 #radius (not important right now) dRotPeriod -1.00000 #rotation period (minus = days) dObliquity 23.5 -dSemi 1.02 +dSemi 1 dEcc 0.0 #eccentricity of orbit dLongP 0 #pericenter, wrt Earth's position at spring equinox # note that this is the typical value +180, diff --git a/examples/ClimateLand/random.in b/examples/ClimateLand/random.in index 05139b93c..257f5b08d 100755 --- a/examples/ClimateLand/random.in +++ b/examples/ClimateLand/random.in @@ -7,7 +7,7 @@ dMass 3.00316726e-06 #mass of planet dRadius -1.00 #radius (not important right now) dRotPeriod -1.00000 #rotation period (minus = days) dObliquity 23.5 -dSemi 1.02 +dSemi 1 dEcc 0.0 #eccentricity of orbit dLongP 0 #pericenter, wrt Earth's position at spring equinox # note that this is the typical value +180, @@ -30,6 +30,8 @@ dTGlobalInit 14.85 #initial guess at average surface temp iNumYears 4 #number of years (orbits) to run clim model iNStepInYear 80 #number of steps to take in a "year" #dSurfAlbedo 0.35 #average surface albedo (annual model only) +dLandFracMean 0.34 +dLandFracAmp 0.2 #__ice params_________ bIceSheets 1 #enable ice sheets diff --git a/examples/ClimateLand/uniform.in b/examples/ClimateLand/uniform.in index 15c83aa2a..c7423f090 100755 --- a/examples/ClimateLand/uniform.in +++ b/examples/ClimateLand/uniform.in @@ -1,13 +1,13 @@ sName uniform #name of planet saModules poise #what vplanet modules you want to use sGeography uniform -dLandFrac 0.5 +dLandFrac 0.292 dMass 3.00316726e-06 #mass of planet dRadius -1.00 #radius (not important right now) dRotPeriod -1.00000 #rotation period (minus = days) dObliquity 23.5 -dSemi 1.02 +dSemi 1 dEcc 0.0 #eccentricity of orbit dLongP 0 #pericenter, wrt Earth's position at spring equinox # note that this is the typical value +180, diff --git a/examples/ClimateLand/vpl.in b/examples/ClimateLand/vpl.in index 50acb1c16..ec7a23a9b 100755 --- a/examples/ClimateLand/vpl.in +++ b/examples/ClimateLand/vpl.in @@ -15,6 +15,6 @@ saBodyFiles sun.in $ #you must list all input files here (except vpl equatorial.in # One land mass along the equator bDoForward 1 #integrate forward in time bVarDt 1 #use variable time stepping (not relevant to poise) -dEta 0.1 #how much to scale variable time step -dStopTime 1e6 #how long should the integration be -dOutputTime 1e3 #how much output you want +dEta 0.01 #how much to scale variable time step +dStopTime 1e5 #how long should the integration be +dOutputTime 1e2 #how much output you want diff --git a/src/poise.c b/src/poise.c index cd6094569..577ba63ec 100644 --- a/src/poise.c +++ b/src/poise.c @@ -1484,7 +1484,7 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { options[OPT_LANDFRACMEAN].bMultiFile = 1; fnRead[OPT_LANDFRACMEAN] = &ReadLandFracMean; - fvFormattedString(&options[OPT_LANDFRACAMP].cName, "dLandFrac"); + fvFormattedString(&options[OPT_LANDFRACAMP].cName, "dLandFracAmp"); fvFormattedString(&options[OPT_LANDFRACAMP].cDescr, "Max difference in land fraction from the mean"); fvFormattedString(&options[OPT_LANDFRACAMP].cDefault, "0.34"); fvFormattedString(&options[OPT_LANDFRACAMP].cDimension, "nd"); @@ -2114,17 +2114,17 @@ void fvWriteLandFracFile(BODY *body,int iBody) { int iLat; double dInterval,dLatNow; FILE *fp; - char *cOut; + char *cOut = NULL; - dInterval = 180/body[iBody].iNumLats; + dInterval = 180./body[iBody].iNumLats; fvFormattedString(&cOut, "%s.geography", body[iBody].cName); - fp = open(cOut,"w"); + fp = fopen(cOut,"w"); for (iLat=0;iLatbDoNeg[iBody]) { fvFormattedString(cUnit, output->cNeg); - } else { - //*dTmp /= fdUnitsMass(units->iMass)/pow(fdUnitsLength(units->iLength),2); - // fsUnitsEnergyFlux(units,cUnit); + *dTmp /= fdUnitsMass(units->iMass)/pow(fdUnitsLength(units->iLength),2); + fsUnitsMass(units->iMass,cUnit); } } From b676dde07746f800fcd70e0f92eb8d585edb2991 Mon Sep 17 00:00:00 2001 From: Rory Barnes Date: Fri, 6 Sep 2024 17:08:03 -0700 Subject: [PATCH 06/12] Set the mean land fraction in examples/ClimateLand/random.in to 0.292. --- examples/ClimateLand/random.in | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/ClimateLand/random.in b/examples/ClimateLand/random.in index 257f5b08d..3c28f276a 100755 --- a/examples/ClimateLand/random.in +++ b/examples/ClimateLand/random.in @@ -30,7 +30,7 @@ dTGlobalInit 14.85 #initial guess at average surface temp iNumYears 4 #number of years (orbits) to run clim model iNStepInYear 80 #number of steps to take in a "year" #dSurfAlbedo 0.35 #average surface albedo (annual model only) -dLandFracMean 0.34 +dLandFracMean 0.292 dLandFracAmp 0.2 #__ice params_________ From a1ede1a873307bac4283a8c1c332b7500255024e Mon Sep 17 00:00:00 2001 From: Rory Barnes Date: Fri, 13 Sep 2024 10:04:57 -0700 Subject: [PATCH 07/12] Added template makeplot.py file --- examples/ClimateLand/makeplot.py | 335 +++++++++++++++++++++++++++++++ 1 file changed, 335 insertions(+) create mode 100644 examples/ClimateLand/makeplot.py diff --git a/examples/ClimateLand/makeplot.py b/examples/ClimateLand/makeplot.py new file mode 100644 index 000000000..083bf35d9 --- /dev/null +++ b/examples/ClimateLand/makeplot.py @@ -0,0 +1,335 @@ +import glob +import pathlib +import re +import sys + +import astropy.units as u +import matplotlib as mpl +import matplotlib.pyplot as plt +import numpy as np +import vplot + +import vplanet + +# Path hacks +path = pathlib.Path(__file__).parents[0].absolute() +sys.path.insert(1, str(path.parents[0])) +from get_args import get_args + + +def clim_evol(plname, xrange=False, show=True): + """ + Creates plots of insolation, temperature, albedo, ice mass, + and bed rock height over the length of the simulation + + Parameters + ---------- + plname : string + The name of the planet with .Climate data + + Keyword Arguments + ----------------- + xrange : float tuple, list, or numpy array + Range of x-values (time) to restrict plot + (default = False (no restriction)) + show : bool + Show plot in Python (default = True) + + """ + fig = plt.figure(figsize=(6.5, 9)) + fig.subplots_adjust(wspace=0.5, hspace=0.5) + + # Run vplanet + out = vplanet.run(path / "vpl.in") + + ctmp = 0 + for p in range(len(out.bodies)): + if out.bodies[p].name == plname: + body = out.bodies[p] + ctmp = 1 + else: + if p == len(out.bodies) - 1 and ctmp == 0: + raise Exception("Planet %s not found." % plname) + + try: + ecc = body.Eccentricity + except: + ecc = np.zeros_like(body.Time) + getattr(out.log.initial, plname).Eccentricity + + try: + inc = body.Inc + except: + inc = np.zeros_like(body.Time) + + try: + obl = body.Obliquity + except: + obltmp = getattr(out.log.initial, plname).Obliquity + if obltmp.unit == "rad": + obltmp *= 180 / np.pi + obl = np.zeros_like(body.Time) + obltmp + + f = open(path / f"{plname}.in", "r") + lines = f.readlines() + f.close() + pco2 = 0 + + for i in range(len(lines)): + if lines[i].split() != []: + if lines[i].split()[0] == "dRotPeriod": + P = -1 * float(lines[i].split()[1]) + if lines[i].split()[0] == "dSemi": + semi = float(lines[i].split()[1]) + if semi < 0: + semi *= -1 + if lines[i].split()[0] == "dpCO2": + pco2 = float(lines[i].split()[1]) + + try: + longp = (body.ArgP + body.LongA + body.PrecA) * np.pi / 180.0 + except: + longp = body.PrecA * np.pi / 180.0 + + esinv = ecc * np.sin(longp) * np.sin(obl * np.pi / 180.0) + + lats = np.unique(body.Latitude) + nlats = len(lats) + ntimes = len(body.Time) + + # plot temperature + temp = np.reshape(body.TempLat, (ntimes, nlats)) + ax1 = plt.subplot(5, 1, 1) + + c = plt.contourf(body.Time, lats, temp.T, cmap="plasma") + plt.ylabel(r"Latitude [$^\circ$]", fontsize=10) + plt.title(r"Surface Temp [$^{\circ}$C]", fontsize=12) + plt.ylim(-90, 90) + plt.yticks([-60, -30, 0, 30, 60], fontsize=9) + plt.xticks(fontsize=9) + if xrange == False: + left = 0 + else: + left = xrange[0] + if xrange: + plt.xlim(xrange) + cbar = plt.colorbar(c, ax=ax1) + plt.setp(cbar.ax.yaxis.get_ticklabels(), fontsize=9) + + # plot albedo + alb = np.reshape(body.AlbedoLat, (ntimes, nlats)) + ax2 = plt.subplot(5, 1, 3) + c = plt.contourf(body.Time, lats, alb.T, cmap="Blues_r") + plt.ylabel(r"Latitude [$^\circ$]", fontsize=10) + plt.title("Albedo [TOA]", fontsize=12) + plt.ylim(-90, 90) + plt.yticks([-60, -30, 0, 30, 60], fontsize=9) + plt.xticks(fontsize=9) + if xrange: + plt.xlim(xrange) + cbar = plt.colorbar(c, ax=ax2) + plt.setp(cbar.ax.yaxis.get_ticklabels(), fontsize=9) + + # plot ice height + ice = np.reshape(body.IceHeight, (ntimes, nlats)) + ax3 = plt.subplot(5, 1, 4) + c = plt.contourf(body.Time, lats, ice.T, cmap="Blues_r") + plt.ylabel(r"Latitude [$^\circ$]", fontsize=10) + plt.title("Ice sheet height [m]", fontsize=12) + plt.ylim(-90, 90) + plt.yticks([-60, -30, 0, 30, 60], fontsize=9) + plt.xticks(fontsize=9) + if xrange: + plt.xlim(xrange) + cbar = plt.colorbar(c, ax=ax3) + plt.setp(cbar.ax.yaxis.get_ticklabels(), fontsize=9) + + # plot bedrock + brock = np.reshape(body.BedrockH, (ntimes, nlats)) + ax4 = plt.subplot(5, 1, 5) + c = plt.contourf(body.Time, lats, brock.T, cmap="Reds_r") + plt.ylabel(r"Latitude [$^\circ$]", fontsize=10) + plt.title("Bedrock height [m]", fontsize=12) + plt.ylim(-90, 90) + plt.yticks([-60, -30, 0, 30, 60], fontsize=9) + plt.xlabel("Time [years]", fontsize=10) + plt.xticks(fontsize=9) + if xrange: + plt.xlim(xrange) + cbar = plt.colorbar(c, ax=ax4) + plt.setp(cbar.ax.yaxis.get_ticklabels(), fontsize=9) + + # plot insolation + insol = np.reshape(body.AnnInsol, (ntimes, nlats)) + ax5 = plt.subplot(5, 1, 2) + c = plt.contourf(body.Time, lats, insol.T, cmap="plasma") + plt.ylabel(r"Latitude [$^\circ$]", fontsize=10) + plt.title(r"Annual average insolation [W/m$^2$]", fontsize=12) + plt.ylim(-90, 90) + plt.yticks([-60, -30, 0, 30, 60], fontsize=9) + plt.xticks(fontsize=9) + if xrange: + plt.xlim(xrange) + cbar = plt.colorbar(c, ax=ax5) + plt.setp(cbar.ax.yaxis.get_ticklabels(), fontsize=9) + + # Save + ext = get_args().ext + plt.savefig(path / f"IceBeltClimateEvol.{ext}", dpi=300) + + if show: + plt.show() + else: + plt.close() + + +def seasonal_maps(time, show=True): + """ + Creates plots of insolation, temperature, and ice balance + over the course of an orbit (4 orbits for temp) + + Parameters + ---------- + time : float + The time of the data you want to plot (see SeasonalClimateFiles directory) + + Keyword Arguments + ----------------- + show : bool + Show plot in Python (default = True) + + """ + + check = 0 + for f in glob.glob(str(path / "SeasonalClimateFiles" / "*.DailyInsol.*")): + + f1 = f.split(".") + + if len(f1) == 4: + timestamp = f1[3] + elif len(f1) == 5: + timestamp = f1[3] + "." + f1[4] + + time0 = float(timestamp) + + if time0 == time: + # get system and planet names + sysname = f1[0] + plname = f1[1] + insolf = ( + path + / "SeasonalClimateFiles" + / f"{sysname}.{plname}.DailyInsol.{timestamp}" + ) + tempf = ( + path + / "SeasonalClimateFiles" + / f"{sysname}.{plname}.SeasonalTemp.{timestamp}" + ) + icef = ( + path + / "SeasonalClimateFiles" + / f"{sysname}.{plname}.SeasonalIceBalance.{timestamp}" + ) + check = 1 + + if check == 0: + raise StandardError("Climate data not found for time %f" % time) + + insol = np.loadtxt(insolf, unpack=True) + temp = np.loadtxt(tempf, unpack=True) + ice = np.loadtxt(icef, unpack=True) + output = vplanet.run(path / "vpl.in") + + ctmp = 0 + for p in range(len(output.bodies)): + if output.bodies[p].name == plname: + body = output.bodies[p] + ctmp = 1 + else: + if p == len(output.bodies) - 1 and ctmp == 0: + raise Exception("Planet %s not found" % plname) + + lats = np.unique(body.Latitude) + try: + obl = body.Obliquity[np.where(body.Time == time)[0]] + except: + obl = getattr(output.log.initial, plname).Obliquity + if obl.unit == "rad": + obl *= 180 / np.pi + + try: + ecc = body.Eccentricity[np.where(body.Time == time)[0]] + except: + ecc = getattr(output.log.initial, plname).Eccentricity + + try: + longp = (body.LongP + body.PrecA)[np.where(body.Time == time)[0]] + except: + try: + longp = getattr(output.log.initial, plname).LongP + except: + try: + longp = ( + getattr(output.log.initial, plname).LongA + + getattr(out.log.initial, plname).ArgP + ) + if longp.unit == "rad": + longp *= 180 / np.pi + longp = longp % 360 + except: + longp = 0 + + fig = plt.figure(figsize=(9, 3.5)) + plt.subplot(1, 3, 1) + fig.subplots_adjust(wspace=0.5) + + plt.title(r"Insolation [W/m$^2$]", fontsize=12) + c1 = plt.contourf(np.arange(np.shape(insol)[1]), lats, insol, cmap="plasma") + cbar = plt.colorbar(c1) + plt.ylim(lats[0], lats[-1]) + plt.xlabel("Time [days]", fontsize=10) + plt.ylabel(r"Latitude [$^\circ$]", fontsize=10) + plt.xticks(fontsize=9) + plt.yticks(fontsize=9) + plt.setp(cbar.ax.yaxis.get_ticklabels(), fontsize=9) + + scale = 4 * np.shape(insol)[1] / np.shape(temp)[1] + plt.subplot(1, 3, 2) + c2 = plt.contourf(np.arange(np.shape(temp)[1]) * scale, lats, temp, cmap="plasma") + plt.title(r"Surface Temp [$^{\circ}$C]", fontsize=12) + cbar = plt.colorbar(c2) + plt.ylim(lats[0], lats[-1]) + plt.ylabel(r"Latitude [$^\circ$]", fontsize=10) + plt.xlabel("Time [days]", fontsize=10) + plt.xlim(0, np.shape(temp)[1] * scale / 4.0) + plt.xticks(fontsize=9) + plt.yticks(fontsize=9) + plt.setp(cbar.ax.yaxis.get_ticklabels(), fontsize=9) + + scale = np.shape(insol)[1] / np.shape(ice)[1] + plt.subplot(1, 3, 3) + c3 = plt.contourf(np.arange(np.shape(ice)[1]) * scale, lats, ice, cmap="Blues_r") + plt.title(r"Ice balance [kg/m$^2$/s]", fontsize=12) + cbar = plt.colorbar(c3) + plt.setp(cbar.ax.yaxis.get_ticklabels(), fontsize=9) + plt.ylim(lats[0], lats[-1]) + plt.xlabel("Time [days]", fontsize=10) + plt.ylabel(r"Latitude [$^\circ$]", fontsize=10) + plt.xticks(fontsize=9) + plt.yticks(fontsize=9) + + # Save + ext = get_args().ext + fig.savefig(path / f"IceBeltSeasonal.{ext}", dpi=300) + + if show: + plt.show() + else: + plt.close() + + +# makes the plots +print("Making evolution plot.") +clim_evol("earth", show=False) +print("Making seasonal plot.") +seasonal_maps(0, show=False) From 97954f766289edd5bac5b9617d853e00ca36728a Mon Sep 17 00:00:00 2001 From: RoryBarnes Date: Sat, 14 Sep 2024 17:29:38 -0700 Subject: [PATCH 08/12] Refactored VerifyGeography --- src/poise.c | 182 +++++++++++++++++++++++++++++++++++++--------------- src/poise.h | 52 +++++++-------- 2 files changed, 158 insertions(+), 76 deletions(-) diff --git a/src/poise.c b/src/poise.c index 7b580f851..112f6ad2e 100644 --- a/src/poise.c +++ b/src/poise.c @@ -607,15 +607,15 @@ void ReadGeography(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, //CheckDuplication(files,options,files->Infile[iFile].cIn,lTmp,\ control->Io.iVerbose); if (!memcmp(sLower(cTmp), "u", 1)) { - body[iFile - 1].iGeography = LANDWATERUNIFORM; + body[iFile - 1].iGeography = GEOGRAPHYUNIFORM; } else if (!memcmp(sLower(cTmp), "m", 1)) { - body[iFile - 1].iGeography = LANDWATERMODERN; + body[iFile - 1].iGeography = GEOGRAPHYMODERN; } else if (!memcmp(sLower(cTmp), "r", 1)) { - body[iFile - 1].iGeography = LANDWATERRANDOM; + body[iFile - 1].iGeography = GEOGRAPHYRANDOM; } else if (!memcmp(sLower(cTmp), "p", 1)) { - body[iFile - 1].iGeography = LANDWATERPOLAR; + body[iFile - 1].iGeography = GEOGRAPHYPOLAR; } else if (!memcmp(sLower(cTmp), "e", 1)) { - body[iFile - 1].iGeography = LANDWATEREQUATORIAL; + body[iFile - 1].iGeography = GEOGRAPHYEQUATORIAL; } else { if (control->Io.iVerbose >= VERBERR) { fprintf(stderr, @@ -1606,7 +1606,7 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fvFormattedString(&options[OPT_GEOGRAPHY].cDescr, "Type of land distribution"); fvFormattedString(&options[OPT_GEOGRAPHY].cDefault, "uniform"); - options[OPT_GEOGRAPHY].dDefault = LANDWATERUNIFORM; + options[OPT_GEOGRAPHY].dDefault = GEOGRAPHYUNIFORM; options[OPT_GEOGRAPHY].iType = 3; options[OPT_GEOGRAPHY].bMultiFile = 1; fnRead[OPT_GEOGRAPHY] = &ReadGeography; @@ -2205,15 +2205,15 @@ void InitializeLandWater(BODY *body, int iBody) { body[iBody].daLandFrac = malloc(body[iBody].iNumLats * sizeof(double)); body[iBody].daWaterFrac = malloc(body[iBody].iNumLats * sizeof(double)); - if (body[iBody].iGeography == LANDWATERUNIFORM) { + if (body[iBody].iGeography == GEOGRAPHYUNIFORM) { fvInitializeLandWaterUniform(body, iBody); - } else if (body[iBody].iGeography == LANDWATERMODERN) { + } else if (body[iBody].iGeography == GEOGRAPHYMODERN) { fvInitializeLandWaterModern(body, iBody); - } else if (body[iBody].iGeography == LANDWATERRANDOM) { + } else if (body[iBody].iGeography == GEOGRAPHYRANDOM) { fvInitializeLandWaterRandom(body, iBody); - } else if (body[iBody].iGeography == LANDWATERPOLAR) { + } else if (body[iBody].iGeography == GEOGRAPHYPOLAR) { fvInitializeLandWaterPolar(body, iBody); - } else if (body[iBody].iGeography == LANDWATEREQUATORIAL) { + } else if (body[iBody].iGeography == GEOGRAPHYEQUATORIAL) { fvInitializeLandWaterEquatorial(body, iBody); } @@ -2898,58 +2898,140 @@ void NullPoiseDerivatives(BODY *body, EVOLVE *evolve, UPDATE *update, // Nothing here, because entire climate simulation is ran in ForceBehavior } -void VerifyGeography(BODY *body, CONTROL *control, OPTIONS *options, - SYSTEM *system, int iBody) { +void fvGeographyExitLandWaterLatitude(CONTROL *control, OPTIONS *options, + int iBody) { + if (control->Io.iVerbose > VERBINPUT) { + fprintf(stderr, + "ERROR: Cannot set %s with \"uniform\", \"modern\", or " + "\"random\" geography. " + "Note that \"uniform\" is the default.\n", + options[OPT_LANDWATERLATITUDE].cName); + } + DoubleLineExit(options[OPT_LANDWATERLATITUDE].cFile[iBody + 1], + options[OPT_GEOGRAPHY].cFile[iBody + 1], + options[OPT_LANDWATERLATITUDE].iLine[iBody + 1], + options[OPT_GEOGRAPHY].iLine[iBody + 1]); +} - if (body[iBody].iGeography == LANDWATERUNIFORM || - body[iBody].iGeography == LANDWATERMODERN || - body[iBody].iGeography == LANDWATERRANDOM) { - if (options[OPT_LANDWATERLATITUDE].iLine[iBody + 1] != -1) { - if (control->Io.iVerbose > VERBINPUT) { - fprintf(stderr, - "ERROR: Cannot set %s with \"uniform\", \"modern\", or " - "\"random\" geography. " - "Note that \"uniform\" is the default.\n", - options[OPT_LANDWATERLATITUDE].cName); - } - DoubleLineExit(options[OPT_LANDWATERLATITUDE].cFile[iBody + 1], - options[OPT_GEOGRAPHY].cFile[iBody + 1], - options[OPT_LANDWATERLATITUDE].iLine[iBody + 1], - options[OPT_GEOGRAPHY].iLine[iBody + 1]); - } +void fvGeographyExitLandFracMean(CONTROL *control, OPTIONS *options, + int iBody) { + if (control->Io.iVerbose > VERBINPUT) { + fprintf(stderr, + "ERROR: Cannot set %s with \"uniform\", \"modern\", \"polar\" or " + "\"equatorial\" geography. " + "Note that \"uniform\" is the default.\n", + options[OPT_LANDWATERLATITUDE].cName); } + DoubleLineExit(options[OPT_LANDFRACMEAN].cFile[iBody + 1], + options[OPT_GEOGRAPHY].cFile[iBody + 1], + options[OPT_LANDFRACMEAN].iLine[iBody + 1], + options[OPT_GEOGRAPHY].iLine[iBody + 1]); +} - if (body[iBody].iGeography == LANDWATERMODERN || - body[iBody].iGeography == LANDWATERRANDOM || - body[iBody].iGeography == LANDWATERPOLAR || - body[iBody].iGeography == LANDWATEREQUATORIAL) { - if (options[OPT_LANDFRAC].iLine[iBody + 1] != -1) { - if (control->Io.iVerbose > VERBINPUT) { - fprintf(stderr, - "ERROR: %s can only be set when %s is set to \"uniform\". Note " - "that \"uniform\" is the default.\n", - options[OPT_LANDFRAC].cName, options[OPT_GEOGRAPHY].cName); - } - } +void fvGeographyExitLandFracAmp(CONTROL *control, OPTIONS *options, int iBody) { + if (control->Io.iVerbose > VERBINPUT) { + fprintf(stderr, + "ERROR: Cannot set %s with \"uniform\", \"modern\", \"polar\" or " + "\"equatorial\" geography. " + "Note that \"uniform\" is the default.\n", + options[OPT_LANDWATERLATITUDE].cName); } + DoubleLineExit(options[OPT_LANDFRACAMP].cFile[iBody + 1], + options[OPT_GEOGRAPHY].cFile[iBody + 1], + options[OPT_LANDFRACAMP].iLine[iBody + 1], + options[OPT_GEOGRAPHY].iLine[iBody + 1]); +} + +void fvGeographyExitLandFrac(CONTROL *control, OPTIONS *options, int iBody) { + if (control->Io.iVerbose > VERBINPUT) { + fprintf(stderr, + "ERROR: %s can only be set when %s is set to \"uniform\". Note " + "that \"uniform\" is the default.\n", + options[OPT_LANDFRAC].cName, options[OPT_GEOGRAPHY].cName); + } + DoubleLineExit(options[OPT_LANDFRAC].cFile[iBody + 1], + options[OPT_GEOGRAPHY].cFile[iBody + 1], + options[OPT_LANDFRAC].iLine[iBody + 1], + options[OPT_GEOGRAPHY].iLine[iBody + 1]); +} - if (body[iBody].iGeography == LANDWATERRANDOM && - options[OPT_RANDSEED].iLine[iBody + 1] == -1) { +void fvVerifyLatitudeMeanAndAmpNotSet(CONTROL *control, OPTIONS *options, + int iBody) { + if (options[OPT_LANDWATERLATITUDE].iLine[iBody + 1] != -1) { + fvGeographyExitLandWaterLatitude(control, options, iBody); + } + if (options[OPT_LANDFRACMEAN].iLine[iBody + 1] != -1) { + fvGeographyExitLandFracMean(control, options, iBody); + } + if (options[OPT_LANDFRACAMP].iLine[iBody + 1] != -1) { + fvGeographyExitLandFracAmp(control, options, iBody); + } +} + +void fvVerifyGeographyUniform(CONTROL *control, OPTIONS *options, int iBody) { + fvVerifyLatitudeMeanAndAmpNotSet(control, options, iBody); +} + +void fvVerifyGeographyModern(CONTROL *control, OPTIONS *options, int iBody) { + fvVerifyLatitudeMeanAndAmpNotSet(control, options, iBody); + if (options[OPT_LANDFRAC].iLine[iBody + 1] != -1) { + fvGeographyExitLandFrac(control, options, iBody); + } +} + +void fvVerifyGeographyRandom(CONTROL *control, OPTIONS *options, SYSTEM *system, + int iBody) { + fvVerifyLatitudeMeanAndAmpNotSet(control, options, iBody); + if (options[OPT_RANDSEED].iLine[iBody + 1] == -1) { if (control->Io.iVerbose > VERBINPUT) { fprintf(stderr, - "\nWARNING: %s is set to \"radnom\", but %s is not set.\n\n", + "\nWARNING: %s is set to \"random\", but %s is not set.\n\n", options[OPT_GEOGRAPHY].cName, options[OPT_RANDSEED].cName); } } + srand(system->iSeed); +} - if (body[iBody].iGeography == LANDWATERRANDOM) { - srand(system->iSeed); +void fvVerifyGeographyPolarEquatorial(BODY *body, CONTROL *control, + OPTIONS *options, int iBody) { + if (options[OPT_LANDFRACMEAN].iLine[iBody + 1] != -1) { + fvGeographyExitLandFracMean(control, options, iBody); + } + if (options[OPT_LANDFRACAMP].iLine[iBody + 1] != -1) { + fvGeographyExitLandFracAmp(control, options, iBody); + } + if (options[OPT_LANDFRAC].iLine[iBody + 1] != -1) { + fvGeographyExitLandFrac(control, options, iBody); } + if (options[OPT_LANDWATERLATITUDE].iLine[iBody + 1] != -1) { + if (control->Io.iVerbose > VERBINPUT) { + fprintf(stderr, "\nWARNING: %s set to \"", options[OPT_GEOGRAPHY].cName); + if (body[iBody].iGeography == GEOGRAPHYPOLAR) { + fprintf(stderr, "polar"); + } else if (body[iBody].iGeography == GEOGRAPHYEQUATORIAL) { + fprintf(stderr, "equatorial"); + } + fprintf(stderr, "\" but %s is not set.\n", + options[OPT_LANDWATERLATITUDE].cName); + } + } + + body[iBody].iLatLandWater = + (int)(body[iBody].dLatLandWater * body[iBody].iNumLats / PI); +} + +void VerifyGeography(BODY *body, CONTROL *control, OPTIONS *options, + SYSTEM *system, int iBody) { - if (body[iBody].iGeography == LANDWATERPOLAR || - body[iBody].iGeography == LANDWATEREQUATORIAL) { - body[iBody].iLatLandWater = - (int)(body[iBody].dLatLandWater * body[iBody].iNumLats / PI); + if (body[iBody].iGeography == GEOGRAPHYUNIFORM) { + fvVerifyGeographyUniform(control, options, iBody); + } else if (body[iBody].iGeography == GEOGRAPHYMODERN) { + fvVerifyGeographyModern(control, options, iBody); + } else if (body[iBody].iGeography == GEOGRAPHYRANDOM) { + fvVerifyGeographyRandom(control, options, system, iBody); + } else if (body[iBody].iGeography == GEOGRAPHYPOLAR || + body[iBody].iGeography == GEOGRAPHYEQUATORIAL) { + fvVerifyGeographyPolarEquatorial(body, control, options, iBody); } } diff --git a/src/poise.h b/src/poise.h index fffb544ba..3c0c13757 100644 --- a/src/poise.h +++ b/src/poise.h @@ -26,11 +26,11 @@ #define ALBTAYLOR 1 /* Land Geography */ -#define LANDWATERUNIFORM 0 -#define LANDWATERMODERN 1 -#define LANDWATERRANDOM 2 -#define LANDWATERPOLAR 3 -#define LANDWATEREQUATORIAL 4 +#define GEOGRAPHYUNIFORM 0 +#define GEOGRAPHYMODERN 1 +#define GEOGRAPHYRANDOM 2 +#define GEOGRAPHYPOLAR 3 +#define GEOGRAPHYEQUATORIAL 4 // Constants for the ice model #define LFICE 3.34e5 // ??? @@ -83,8 +83,8 @@ #define OPT_SKIPSEASENABLED 1921 #define OPT_DIFFROT 1922 #define OPT_SPINUPTOL 1923 -#define OPT_READORBITOBLDATA 1924 -#define OPT_FILEORBITOBLDATA 1925 +#define OPT_READORBITOBLDATA 1924 +#define OPT_FILEORBITOBLDATA 1925 #define OPT_LANDFRAC 1940 #define OPT_LANDFRACMEAN 1941 @@ -92,7 +92,7 @@ #define OPT_HEATCAPLAND 1943 #define OPT_HEATCAPWATER 1944 #define OPT_FRZTSEAICE 1945 -//#define OPT_LATENTHEAT 1945 +// #define OPT_LATENTHEAT 1945 #define OPT_ICECONDUCT 1946 #define OPT_MIXINGDEPTH 1947 #define OPT_NULANDWATER 1948 @@ -198,7 +198,7 @@ void HelpOptionsPoise(OPTIONS *); void InitializeOptionsPoise(OPTIONS *, fnReadOption[]); void ReadOptionsPoise(BODY *, CONTROL *, FILES *, OPTIONS *, SYSTEM *, fnReadOption[], int); -void ReadOrbitOblData(BODY*,CONTROL*,FILES*,OPTIONS*,SYSTEM*,int); +void ReadOrbitOblData(BODY *, CONTROL *, FILES *, OPTIONS *, SYSTEM *, int); /* Verify Functions */ void VerifyPoise(BODY *, CONTROL *, FILES *, OPTIONS *, OUTPUT *, SYSTEM *, @@ -217,39 +217,39 @@ void FinalizeUpdateIceMassPoise(BODY *, UPDATE *, int *, int, int, int); void HelpOutputPoise(OUTPUT *); void WriteTGlobal(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *, UPDATE *, int, - double *, char**); + double *, char **); void WriteAlbedoGlobal(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *, UPDATE *, - int, double *, char**); + int, double *, char **); void WriteTempLat(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *, UPDATE *, int, - double *, char**); + double *, char **); void WriteTempMinLW(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *, UPDATE *, - int, double *, char**); + int, double *, char **); void WriteTempMaxLW(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *, UPDATE *, - int, double *, char**); + int, double *, char **); void WriteAlbedoLat(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *, UPDATE *, - int, double *, char**); + int, double *, char **); void WriteAnnualInsol(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *, UPDATE *, - int, double *, char**); + int, double *, char **); void WriteDailyInsol(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *, UPDATE *, - int, double *, char**); + int, double *, char **); void WritePlanckB(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *, UPDATE *, int, - double *, char**); + double *, char **); void WritePlanckBAvg(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *, UPDATE *, - int, double *, char**); + int, double *, char **); void WriteSeasonalTemp(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *, UPDATE *, - int, double *, char**); + int, double *, char **); void WriteSeasonalFluxes(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *, - UPDATE *, int, double *, char**); + UPDATE *, int, double *, char **); void WriteSeasonalIceBalance(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *, - UPDATE *, int, double *, char**); + UPDATE *, int, double *, char **); void WriteFluxMerid(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *, UPDATE *, - int, double *, char**); + int, double *, char **); void WriteFluxIn(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *, UPDATE *, int, - double *, char**); + double *, char **); void WriteFluxOut(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *, UPDATE *, int, - double *, char**); + double *, char **); void WriteDivFlux(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *, UPDATE *, int, - double *, char**); + double *, char **); void InitializeOutputPoise(OUTPUT *, fnWriteOutput[]); /* Logging Functions */ From 7c1d8689e9c046e6db694d4e55879216995c8c00 Mon Sep 17 00:00:00 2001 From: RoryBarnes Date: Sat, 14 Sep 2024 17:48:46 -0700 Subject: [PATCH 09/12] Fixed a few bugs in VerifyGeography --- src/poise.c | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/poise.c b/src/poise.c index 112f6ad2e..6fbd41dec 100644 --- a/src/poise.c +++ b/src/poise.c @@ -2920,7 +2920,7 @@ void fvGeographyExitLandFracMean(CONTROL *control, OPTIONS *options, "ERROR: Cannot set %s with \"uniform\", \"modern\", \"polar\" or " "\"equatorial\" geography. " "Note that \"uniform\" is the default.\n", - options[OPT_LANDWATERLATITUDE].cName); + options[OPT_LANDFRACMEAN].cName); } DoubleLineExit(options[OPT_LANDFRACMEAN].cFile[iBody + 1], options[OPT_GEOGRAPHY].cFile[iBody + 1], @@ -2934,7 +2934,7 @@ void fvGeographyExitLandFracAmp(CONTROL *control, OPTIONS *options, int iBody) { "ERROR: Cannot set %s with \"uniform\", \"modern\", \"polar\" or " "\"equatorial\" geography. " "Note that \"uniform\" is the default.\n", - options[OPT_LANDWATERLATITUDE].cName); + options[OPT_LANDFRACAMP].cName); } DoubleLineExit(options[OPT_LANDFRACAMP].cFile[iBody + 1], options[OPT_GEOGRAPHY].cFile[iBody + 1], @@ -2981,7 +2981,13 @@ void fvVerifyGeographyModern(CONTROL *control, OPTIONS *options, int iBody) { void fvVerifyGeographyRandom(CONTROL *control, OPTIONS *options, SYSTEM *system, int iBody) { - fvVerifyLatitudeMeanAndAmpNotSet(control, options, iBody); + + if (options[OPT_LANDWATERLATITUDE].iLine[iBody + 1] != -1) { + fvGeographyExitLandWaterLatitude(control, options, iBody); + } + if (options[OPT_LANDFRAC].iLine[iBody + 1] != -1) { + fvGeographyExitLandFrac(control, options, iBody); + } if (options[OPT_RANDSEED].iLine[iBody + 1] == -1) { if (control->Io.iVerbose > VERBINPUT) { fprintf(stderr, @@ -3003,7 +3009,7 @@ void fvVerifyGeographyPolarEquatorial(BODY *body, CONTROL *control, if (options[OPT_LANDFRAC].iLine[iBody + 1] != -1) { fvGeographyExitLandFrac(control, options, iBody); } - if (options[OPT_LANDWATERLATITUDE].iLine[iBody + 1] != -1) { + if (options[OPT_LANDWATERLATITUDE].iLine[iBody + 1] == -1) { if (control->Io.iVerbose > VERBINPUT) { fprintf(stderr, "\nWARNING: %s set to \"", options[OPT_GEOGRAPHY].cName); if (body[iBody].iGeography == GEOGRAPHYPOLAR) { From 3bcdcbdfec60de16591b0a8aae6b9e6daaad0d32 Mon Sep 17 00:00:00 2001 From: RoryBarnes Date: Sun, 29 Sep 2024 10:07:00 -0700 Subject: [PATCH 10/12] Fixed bugs in how the polar and equatorial geographies calculate transitional latitudes. Modern geography is still broken. --- examples/ClimateLand/equatorial.in | 2 +- src/poise.c | 126 +++++++++++++++++------------ src/poise.h | 3 +- src/vplanet.h | 27 +++---- 4 files changed, 89 insertions(+), 69 deletions(-) diff --git a/examples/ClimateLand/equatorial.in b/examples/ClimateLand/equatorial.in index bc57c4a39..73e384558 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 6fbd41dec..05e21769a 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 3c0c13757..333a9487d 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 fad4f4be3..cb44d1ccf 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) */ From e70b0a8be455cd6b46dc1383291e23c9a33030a4 Mon Sep 17 00:00:00 2001 From: RoryBarnes Date: Sun, 29 Sep 2024 17:53:42 -0700 Subject: [PATCH 11/12] Added more columns to the .geography file. Updated calculation of LandFracGlobal, but the mods didn't change the result. --- src/poise.c | 46 +++++++++++++++++++++++++++++++++------------- src/vplanet.h | 7 +++++-- 2 files changed, 38 insertions(+), 15 deletions(-) diff --git a/src/poise.c b/src/poise.c index 05e21769a..b68b06c6c 100644 --- a/src/poise.c +++ b/src/poise.c @@ -2085,15 +2085,28 @@ void VerifyNStepSeasonal(BODY *body, int iBody) { } void InitializeLatGrid(BODY *body, int iBody) { - double dDelta_x, SinLat; - int iLats; - dDelta_x = 2.0 / body[iBody].iNumLats; + double dInterval, dSinMinLat, dSinMidLat, dSinMaxLat, dMaxLat; + int iLat; + + body[iBody].daLats = malloc(body[iBody].iNumLats * sizeof(double)); + body[iBody].daLatsMin = malloc((body[iBody].iNumLats) * sizeof(double)); + body[iBody].daLatsWidth = malloc((body[iBody].iNumLats) * sizeof(double)); + body[iBody].daLatsArea = malloc((body[iBody].iNumLats) * sizeof(double)); - body[iBody].daLats = malloc(body[iBody].iNumLats * sizeof(double)); + dInterval = 2.0 / body[iBody].iNumLats; + for (iLat = 0; iLat < body[iBody].iNumLats; iLat++) { + dSinMinLat = -1.0 + iLat * dInterval; + body[iBody].daLatsMin[iLat] = asin(dSinMinLat); + + dSinMidLat = (-1.0 + 0.5 * dInterval) + iLat * dInterval; + body[iBody].daLats[iLat] = asin(dSinMidLat); - for (iLats = 0; iLats < body[iBody].iNumLats; iLats++) { - SinLat = (-1.0 + dDelta_x / 2.) + iLats * dDelta_x; - body[iBody].daLats[iLats] = asin(SinLat); + dSinMaxLat = -1 + (iLat + 1) * dInterval; + dMaxLat = asin(dSinMaxLat); + body[iBody].daLatsWidth[iLat] = dMaxLat - body[iBody].daLatsMin[iLat]; + body[iBody].daLatsArea[iLat] = 2 * PI * body[iBody].dRadius * + body[iBody].dRadius * + (dSinMaxLat - dSinMinLat); } } @@ -2185,13 +2198,20 @@ void fvInitializeGeographyEquatorial(BODY *body, int iBody) { void fvWriteGeographyFile(BODY *body, int iBody) { int iLat; FILE *fp; - char *cOut = NULL; + char *cOut = NULL; + double dRad2Deg = 180 / PI; fvFormattedString(&cOut, "%s.geography", body[iBody].cName); fp = fopen(cOut, "w"); + fprintf(fp, + "MidLat[deg] LandFrac WaterFrac MinLat[deg] Width[deg] Area[km^2]\n"); for (iLat = 0; iLat < body[iBody].iNumLats; iLat++) { - fprintf(fp, "%.5lf %.5lf %.5lf\n", (body[iBody].daLats[iLat] * 180. / PI), - body[iBody].daLandFrac[iLat], body[iBody].daWaterFrac[iLat]); + fprintf(fp, "%.5lf %.5lf %.5lf %.5lf %lf %lf\n", + (body[iBody].daLats[iLat] * dRad2Deg), body[iBody].daLandFrac[iLat], + body[iBody].daWaterFrac[iLat], + (body[iBody].daLatsMin[iLat] * dRad2Deg), + (body[iBody].daLatsWidth[iLat] * dRad2Deg), + body[iBody].daLatsArea[iLat] / 1e6); } fclose(fp); } @@ -8233,11 +8253,11 @@ void PoiseIceSheets(BODY *body, EVOLVE *evolve, int iBody) { double fdLandFracGlobal(BODY *body, int iBody) { int iLat; - double dLandFrac = 0; + double dLandFrac, dLandArea = 0; for (iLat = 0; iLat < body[iBody].iNumLats; iLat++) { - dLandFrac += body[iBody].daLandFrac[iLat]; + dLandArea += (body[iBody].daLandFrac[iLat] * body[iBody].daLatsArea[iLat]); } - dLandFrac /= body[iBody].iNumLats; + dLandFrac = dLandArea / (4 * PI * body[iBody].dRadius * body[iBody].dRadius); return dLandFrac; } diff --git a/src/vplanet.h b/src/vplanet.h index cb44d1ccf..8e943312a 100644 --- a/src/vplanet.h +++ b/src/vplanet.h @@ -722,8 +722,11 @@ 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), assuming equal areas - per bin; South Pole is 0 */ + double *daLats; /**< Latitude of each cell (centered), assuming equal areas + per bin; South Pole is 0 */ + double *daLatsMin; /**< Lower bound of each latitudinal bin */ + double *daLatsWidth; /**< Angular width of each latitudinal bin */ + double *daLatsArea; /**< Total area of each latitudinal bin */ double *daPeakInsol; /**< Annually averaged insolation at each latitude */ double *daTGrad; /**< Gradient of temperature (meridional) */ From 9d823176295b3a4732769a74d6fba7f26862982d Mon Sep 17 00:00:00 2001 From: Rory Barnes Date: Wed, 9 Oct 2024 14:14:51 -0700 Subject: [PATCH 12/12] Added unit tests for Poise:SeasonalNC79*. Also temporarily set units to the following outputs to non-dimensional to enable testing: - AreaIceCov - IceMass - IceAblate - PlackBAvg - IceAccum --- src/poise.c | 13 +- .../SeasonalNC79Equatorial/equatorial.in | 71 ++++++ tests/Poise/SeasonalNC79Equatorial/sun.in | 9 + .../test_SeasonalNC79Equatorial.py | 228 ++++++++++++++++++ tests/Poise/SeasonalNC79Equatorial/vpl.in | 15 ++ tests/Poise/SeasonalNC79Polar/polar.in | 71 ++++++ tests/Poise/SeasonalNC79Polar/sun.in | 9 + .../test_SeasonalNC79Polar.py | 228 ++++++++++++++++++ tests/Poise/SeasonalNC79Polar/vpl.in | 15 ++ tests/Poise/SeasonalNC79Random/random.in | 73 ++++++ tests/Poise/SeasonalNC79Random/sun.in | 9 + .../test_SeasonalNC79Random.py | 228 ++++++++++++++++++ tests/Poise/SeasonalNC79Random/vpl.in | 15 ++ tests/Poise/SeasonalNC79Uniform/sun.in | 9 + .../test_SeasonalNC79Uniform.py | 228 ++++++++++++++++++ tests/Poise/SeasonalNC79Uniform/uniform.in | 71 ++++++ tests/Poise/SeasonalNC79Uniform/vpl.in | 15 ++ 17 files changed, 1300 insertions(+), 7 deletions(-) create mode 100755 tests/Poise/SeasonalNC79Equatorial/equatorial.in create mode 100755 tests/Poise/SeasonalNC79Equatorial/sun.in create mode 100644 tests/Poise/SeasonalNC79Equatorial/test_SeasonalNC79Equatorial.py create mode 100755 tests/Poise/SeasonalNC79Equatorial/vpl.in create mode 100755 tests/Poise/SeasonalNC79Polar/polar.in create mode 100755 tests/Poise/SeasonalNC79Polar/sun.in create mode 100644 tests/Poise/SeasonalNC79Polar/test_SeasonalNC79Polar.py create mode 100755 tests/Poise/SeasonalNC79Polar/vpl.in create mode 100755 tests/Poise/SeasonalNC79Random/random.in create mode 100755 tests/Poise/SeasonalNC79Random/sun.in create mode 100644 tests/Poise/SeasonalNC79Random/test_SeasonalNC79Random.py create mode 100755 tests/Poise/SeasonalNC79Random/vpl.in create mode 100755 tests/Poise/SeasonalNC79Uniform/sun.in create mode 100644 tests/Poise/SeasonalNC79Uniform/test_SeasonalNC79Uniform.py create mode 100755 tests/Poise/SeasonalNC79Uniform/uniform.in create mode 100755 tests/Poise/SeasonalNC79Uniform/vpl.in diff --git a/src/poise.c b/src/poise.c index b68b06c6c..3a695bac7 100644 --- a/src/poise.c +++ b/src/poise.c @@ -3700,13 +3700,7 @@ void WriteAreaIceCov(BODY *body, CONTROL *control, OUTPUT *output, fvAreaIceCovered(body, iBody); *dTmp = body[iBody].dAreaIceCov; - // if (output->bDoNeg[iBody]) { - // // Negative option is SI - // fvFormattedString(cUnit,output->cNeg); - // } else { - // *dTmp /= fdUnitsMass(units->iMass); - // fsUnitsMass(units->iMass,cUnit); - // } + fvFormattedString(cUnit,""); } void WriteIceBalanceTot(BODY *body, CONTROL *control, OUTPUT *output, @@ -4186,6 +4180,7 @@ void WriteIceMass(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, if (output->bDoNeg[iBody]) { fvFormattedString(cUnit, output->cNeg); } else { + fvFormattedString(cUnit,""); /* XXX Need to fix units! *dTmp /= fdUnitsMass(units->iMass)/pow(fdUnitsLength(units->iLength),2); fsUnitsMass(units->iMass,cUnit); @@ -4239,6 +4234,8 @@ void WritePlanckBAvg(BODY *body, CONTROL *control, OUTPUT *output, *dTmp = 0.0; } + fvFormattedString(cUnit,""); + // if (output->bDoNeg[iBody]) { // fvFormattedString(cUnit,output->cNeg); // } else { @@ -4277,6 +4274,7 @@ void WriteIceAccum(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, *dTmp = 0.0; } + fvFormattedString(cUnit,""); // if (output->bDoNeg[iBody]) { // fvFormattedString(cUnit,output->cNeg); // } else { @@ -4296,6 +4294,7 @@ void WriteIceAblate(BODY *body, CONTROL *control, OUTPUT *output, *dTmp = 0.0; } + fvFormattedString(cUnit,""); // if (output->bDoNeg[iBody]) { // fvFormattedString(cUnit,output->cNeg); // } else { diff --git a/tests/Poise/SeasonalNC79Equatorial/equatorial.in b/tests/Poise/SeasonalNC79Equatorial/equatorial.in new file mode 100755 index 000000000..73e384558 --- /dev/null +++ b/tests/Poise/SeasonalNC79Equatorial/equatorial.in @@ -0,0 +1,71 @@ +sName equatorial #name of planet +saModules poise #what vplanet modules you want to use +sGeography equitorial +dLandWaterLatitude 16.4 + +dMass 3.00316726e-06 #mass of planet +dRadius -1.00 #radius (not important right now) +dRotPeriod -1.00000 #rotation period (minus = days) +dObliquity 23.5 +dSemi 1 +dEcc 0.0 #eccentricity of orbit +dLongP 0 #pericenter, wrt Earth's position at spring equinox +# note that this is the typical value +180, +# since that one is solar position +dDynEllip 0.0 #shape of planet (0 = a sphere) +dPrecA 0.0 #orientation of spin axis + +#_______addition disorb/distrot parameters (leave these alone for now)__________________ +#dInc 5e-5 #inclination of orbit +#dLongA 348.73936 #orientation of orbital plane +#bGRCorr 0 #use GR correction (not important) +#bInvPlane 1 #convert to invariable plane coords +#bOverrideMaxEcc 1 #override max ecc halt (not recommended) +#dHaltMaxEcc 0.4 #eccentricity at which to halt simulation + +#_______poise parameters (have fun with these!)_________________________________________ +iLatCellNum 151 #number of latitude cells +sClimateModel sea #use seasonal or annual model +dTGlobalInit 14.85 #initial guess at average surface temp +iNumYears 4 #number of years (orbits) to run clim model +iNStepInYear 80 #number of steps to take in a "year" +#dSurfAlbedo 0.35 #average surface albedo (annual model only) + +#__ice params_________ +bIceSheets 1 #enable ice sheets +dInitIceLat 90. #how low do initial ice sheet extend? +dInitIceHeight 0. #height of initial ice sheets +dIceDepRate 2.25e-5 #rate of snow build up (when T < 0) +dIceAlbedo 0.6 #albedo of ice +iIceDt 1 #time step of ice-sheet model (orbits) +iReRunSeas 500 #how often to re-run seasonal model +bSeaIceModel 0 #use sea ice model (slow!) +bSkipSeasEnabled 0 #can skip seasonal if snowball state present + +#__heat diffusion______ +#bMEPDiff 1 #calculate diffusion using max entropy production +#bHadley 1 #mimic hadley heat diffusion +dDiffusion 0.58 #diffusion coefficient (fixed) +dNuLandWater 0.8 #Heat diffusion coefficient between Land and Water + +#__outgoing flux_______ +dPlanckA 203.3 #offset for OLR calculation (greenhouse) +dPlanckB 2.09 #slope of OLR calc (water vapor feedback) +bCalcAB 0 #calculate A & B from Kasting model fits +#dpCO2 0.00028 #partial pressure of co2 + +#__surface properties__ +dAlbedoLand 0.363 #albedo of land +dAlbedoWater 0.263 #albedo of water +dHeatCapLand 1.55e7 #land heat capacity +dHeatCapWater 4.428e6 #water heat capacity +dMixingDepth 70 #mixing depth of ocean + + +#________output options!_____________________________________________ +saOutputOrder Time PrecA -TGlobal AlbedoGlobal -FluxOutGlobal $ + -TotIceMass -TotIceFlow -TotIceBalance DeltaTime AreaIceCov Snowball Obliq $ + IceBeltLand IceBeltSea Ecce + +saGridOutput Time -Latitude -TempLat AlbedoLat -AnnInsol -FluxIn -FluxOut IceMass -IceHeight DIceMassDt $ + -IceFlow -BedrockH -TempMaxLat -TempMinLat -FluxMerid -DivFlux diff --git a/tests/Poise/SeasonalNC79Equatorial/sun.in b/tests/Poise/SeasonalNC79Equatorial/sun.in new file mode 100755 index 000000000..a1f97c666 --- /dev/null +++ b/tests/Poise/SeasonalNC79Equatorial/sun.in @@ -0,0 +1,9 @@ +# sun parameters +sName sun +dMass 1 +dSemi 0 +dEcc 0 +dRadius 0.00135 +dLuminosity -1 +sStellarModel none #sun does not change over time +saModules stellar #use stellar module (needed for luminosity) diff --git a/tests/Poise/SeasonalNC79Equatorial/test_SeasonalNC79Equatorial.py b/tests/Poise/SeasonalNC79Equatorial/test_SeasonalNC79Equatorial.py new file mode 100644 index 000000000..6c0512b41 --- /dev/null +++ b/tests/Poise/SeasonalNC79Equatorial/test_SeasonalNC79Equatorial.py @@ -0,0 +1,228 @@ +from benchmark import Benchmark, benchmark +import astropy.units as u + +@benchmark( + { + "log.initial.system.Age": {"value": 0.000000, "unit": u.sec}, + "log.initial.system.Time": {"value": 0.000000, "unit": u.sec}, + "log.initial.system.TotAngMom": {"value": 1.501063e+42, "unit": (u.kg * u.m ** 2) / u.sec}, + "log.initial.system.TotEnergy": {"value": -7.839372e+41, "unit": u.Joule}, + "log.initial.system.PotEnergy": {"value": -7.839908e+41, "unit": u.Joule}, + "log.initial.system.KinEnergy": {"value": 5.361272e+37, "unit": u.Joule}, + "log.initial.system.DeltaTime": {"value": 0.000000, "unit": u.sec}, + "log.initial.sun.Mass": {"value": 1.988416e+30, "unit": u.kg}, + "log.initial.sun.Radius": {"value": 2.019571e+08, "unit": u.m}, + "log.initial.sun.RadGyra": {"value": 0.500000}, + "log.initial.sun.RotAngMom": {"value": 1.474456e+42, "unit": (u.kg * u.m ** 2) / u.sec}, + "log.initial.sun.RotVel": {"value": 1.468674e+04, "unit": u.m / u.sec}, + "log.initial.sun.BodyType": {"value": 0.000000}, + "log.initial.sun.RotRate": {"value": 7.272205e-05, "unit": 1 / u.sec}, + "log.initial.sun.RotPer": {"value": 8.640000e+04, "unit": u.sec}, + "log.initial.sun.Density": {"value": 5.762900e+04, "unit": u.kg / u.m ** 3}, + "log.initial.sun.HZLimitDryRunaway": {"value": 1.357831e+11, "unit": u.m}, + "log.initial.sun.HZLimRecVenus": {"value": 1.118929e+11, "unit": u.m}, + "log.initial.sun.HZLimRunaway": {"value": 1.461108e+11, "unit": u.m}, + "log.initial.sun.HZLimMoistGreenhouse": {"value": 1.480517e+11, "unit": u.m}, + "log.initial.sun.HZLimMaxGreenhouse": {"value": 2.509538e+11, "unit": u.m}, + "log.initial.sun.HZLimEarlyMars": {"value": 2.738109e+11, "unit": u.m}, + "log.initial.sun.Instellation": {"value": -1.000000, "unit": u.kg / u.sec ** 3}, + "log.initial.sun.CriticalSemiMajorAxis": {"value": -1.000000, "unit": u.m}, + "log.initial.sun.LXUVTot": {"value": 3.846000e+23, "unit": u.kg / u.sec ** 3}, + "log.initial.sun.LostEnergy": {"value": 5.562685e-309, "unit": u.Joule}, + "log.initial.sun.LostAngMom": {"value": 5.562685e-309, "unit": (u.kg * u.m ** 2) / u.sec}, + "log.initial.sun.EscapeVelocity": {"value": 1.146413e+06, "unit": u.m / u.sec}, + "log.initial.sun.Luminosity": {"value": 3.846000e+26, "unit": u.W}, + "log.initial.sun.LXUVStellar": {"value": 3.846000e+23, "unit": u.W}, + "log.initial.sun.Temperature": {"value": 5778.000000, "unit": u.K}, + "log.initial.sun.LXUVFrac": {"value": 0.001000}, + "log.initial.sun.RossbyNumber": {"value": 0.078260}, + "log.initial.sun.DRotPerDtStellar": {"value": 6.558557e-13}, + "log.initial.sun.WindTorque": {"value": 1.119248e+25}, + "log.initial.equatorial.Mass": {"value": 5.971546e+24, "unit": u.kg}, + "log.initial.equatorial.Obliquity": {"value": 0.410152, "unit": u.rad}, + "log.initial.equatorial.PrecA": {"value": 0.000000, "unit": u.rad}, + "log.initial.equatorial.Radius": {"value": 6.378100e+06, "unit": u.m}, + "log.initial.equatorial.RadGyra": {"value": 0.500000}, + "log.initial.equatorial.BodyType": {"value": 0.000000}, + "log.initial.equatorial.Density": {"value": 5494.449526, "unit": u.kg / u.m ** 3}, + "log.initial.equatorial.HZLimitDryRunaway": {"value": 1.117992e+11, "unit": u.m}, + "log.initial.equatorial.HZLimRecVenus": {"value": 1.118929e+11, "unit": u.m}, + "log.initial.equatorial.HZLimRunaway": {"value": 1.461108e+11, "unit": u.m}, + "log.initial.equatorial.HZLimMoistGreenhouse": {"value": 1.480517e+11, "unit": u.m}, + "log.initial.equatorial.HZLimMaxGreenhouse": {"value": 2.509538e+11, "unit": u.m}, + "log.initial.equatorial.HZLimEarlyMars": {"value": 2.738109e+11, "unit": u.m}, + "log.initial.equatorial.Instellation": {"value": 1367.566935, "unit": u.kg / u.sec ** 3}, + "log.initial.equatorial.Eccentricity": {"value": 0.000000}, + "log.initial.equatorial.MeanMotion": {"value": 1.990987e-07, "unit": 1 / u.sec}, + "log.initial.equatorial.OrbPeriod": {"value": 3.155815e+07, "unit": u.sec}, + "log.initial.equatorial.SemiMajorAxis": {"value": 1.495979e+11, "unit": u.m}, + "log.initial.equatorial.COPP": {"value": 0.000000}, + "log.initial.equatorial.EscapeVelocity": {"value": 1.117931e+04, "unit": u.m / u.sec}, + "log.initial.equatorial.TGlobal": {"value": 16.637409, "unit": u.deg_C}, + "log.initial.equatorial.AlbedoGlobal": {"value": 0.322068}, + "log.initial.equatorial.FluxInGlobal": {"value": 237.998617, "unit": u.kg / u.sec ** 3}, + "log.initial.equatorial.FluxOutGlobal": {"value": 238.072185, "unit": u.W / u.m ** 2}, + "log.initial.equatorial.TotIceMass": {"value": 0.000000, "unit": u.kg}, + "log.initial.equatorial.TotIceFlow": {"value": 0.000000, "unit": u.kg}, + "log.initial.equatorial.TotIceBalance": {"value": 0.000000, "unit": u.kg}, + "log.initial.equatorial.SkipSeas": {"value": 0.000000}, + "log.initial.equatorial.AreaIceCov": {"value": 0.079338}, + "log.initial.equatorial.Latitude": {"value": -83.402352, "unit": u.deg}, + "log.initial.equatorial.TempLat": {"value": -10.684802, "unit": u.deg_C}, + "log.initial.equatorial.AlbedoLat": {"value": 0.600000}, + "log.initial.equatorial.AnnInsol": {"value": 176.064060, "unit": u.W / u.m ** 2}, + "log.initial.equatorial.FluxMerid": {"value": -1.632808, "unit": u.PW}, + "log.initial.equatorial.FluxIn": {"value": 70.432311, "unit": u.W / u.m ** 2}, + "log.initial.equatorial.FluxOut": {"value": 180.968764, "unit": u.W / u.m ** 2}, + "log.initial.equatorial.DivFlux": {"value": -110.461418, "unit": u.W / u.m ** 2}, + "log.initial.equatorial.IceMass": {"value": 0.000000}, + "log.initial.equatorial.IceHeight": {"value": 0.000000, "unit": u.m}, + "log.initial.equatorial.DIceMassDt": {"value": 0.000000, "unit": u.m}, + "log.initial.equatorial.IceFlow": {"value": 0.000000, "unit": u.m / u.sec}, + "log.initial.equatorial.EnergyResL": {"value": -15.700340, "unit": u.kg / u.sec ** 3}, + "log.initial.equatorial.EnergyResW": {"value": 0.154552, "unit": u.kg / u.sec ** 3}, + "log.initial.equatorial.BedrockH": {"value": 0.000000, "unit": u.m}, + "log.initial.equatorial.TempLandLat": {"value": 262.574234, "unit": u.sec}, + "log.initial.equatorial.TempWaterLat": {"value": 262.312582, "unit": u.sec}, + "log.initial.equatorial.AlbedoLandLat": {"value": 0.600000}, + "log.initial.equatorial.AlbedoWaterLat": {"value": 0.600000}, + "log.initial.equatorial.TempMinLat": {"value": -12.488589, "unit": u.deg_C}, + "log.initial.equatorial.TempMaxLat": {"value": -8.948297, "unit": u.deg_C}, + "log.initial.equatorial.Snowball": {"value": 0.000000}, + "log.initial.equatorial.PlanckBAvg": {"value": 2.090000}, + "log.initial.equatorial.IceAccum": {"value": 0.764899}, + "log.initial.equatorial.IceAblate": {"value": 0.000000}, + "log.initial.equatorial.TempMaxLand": {"value": 265.054711, "unit": u.sec}, + "log.initial.equatorial.TempMaxWater": {"value": 264.046713, "unit": u.sec}, + "log.initial.equatorial.PeakInsol": {"value": 541.704677, "unit": u.kg / u.sec ** 3}, + "log.initial.equatorial.IceCapNorthLand": {"value": 1.000000}, + "log.initial.equatorial.IceCapNorthSea": {"value": 1.000000}, + "log.initial.equatorial.IceCapSouthLand": {"value": 1.000000}, + "log.initial.equatorial.IceCapSouthSea": {"value": 1.000000}, + "log.initial.equatorial.IceBeltLand": {"value": 0.000000}, + "log.initial.equatorial.IceBeltSea": {"value": 0.000000}, + "log.initial.equatorial.SnowballLand": {"value": 0.000000}, + "log.initial.equatorial.SnowballSea": {"value": 0.000000}, + "log.initial.equatorial.IceFree": {"value": 0.000000}, + "log.initial.equatorial.IceCapNorthLatLand": {"value": 1.121291, "unit": u.rad}, + "log.initial.equatorial.IceCapNorthLatSea": {"value": 1.152808, "unit": u.rad}, + "log.initial.equatorial.IceCapSouthLatLand": {"value": -1.121291, "unit": u.rad}, + "log.initial.equatorial.IceCapSouthLatSea": {"value": -1.121291, "unit": u.rad}, + "log.initial.equatorial.IceBeltNorthLatLand": {"value": 100.000000, "unit": u.rad}, + "log.initial.equatorial.IceBeltNorthLatSea": {"value": 100.000000, "unit": u.rad}, + "log.initial.equatorial.IceBeltSouthLatLand": {"value": 100.000000, "unit": u.rad}, + "log.initial.equatorial.IceBeltSouthLatSea": {"value": 100.000000, "unit": u.rad}, + "log.initial.equatorial.LandFracGlobal": {"value": 0.289073}, + "log.final.system.Age": {"value": 3.155760e+07, "unit": u.sec, "rtol": 1e-4}, + "log.final.system.Time": {"value": 3.155760e+07, "unit": u.sec, "rtol": 1e-4}, + "log.final.system.TotAngMom": {"value": 1.501063e+42, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, + "log.final.system.TotEnergy": {"value": -7.839372e+41, "unit": u.Joule, "rtol": 1e-4}, + "log.final.system.PotEnergy": {"value": -7.839908e+41, "unit": u.Joule, "rtol": 1e-4}, + "log.final.system.KinEnergy": {"value": 5.361272e+37, "unit": u.Joule, "rtol": 1e-4}, + "log.final.sun.Mass": {"value": 1.988416e+30, "unit": u.kg, "rtol": 1e-4}, + "log.final.sun.Radius": {"value": 2.019571e+08, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.RadGyra": {"value": 0.500000, "rtol": 1e-4}, + "log.final.sun.RotAngMom": {"value": 1.474456e+42, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, + "log.final.sun.RotVel": {"value": 1.468674e+04, "unit": u.m / u.sec, "rtol": 1e-4}, + "log.final.sun.BodyType": {"value": 0.000000, "rtol": 1e-4}, + "log.final.sun.RotRate": {"value": 7.272205e-05, "unit": 1 / u.sec, "rtol": 1e-4}, + "log.final.sun.RotPer": {"value": 8.640000e+04, "unit": u.sec, "rtol": 1e-4}, + "log.final.sun.Density": {"value": 5.762900e+04, "unit": u.kg / u.m ** 3, "rtol": 1e-4}, + "log.final.sun.HZLimitDryRunaway": {"value": 1.357831e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.HZLimRecVenus": {"value": 1.118929e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.HZLimRunaway": {"value": 1.461108e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.HZLimMoistGreenhouse": {"value": 1.480517e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.HZLimMaxGreenhouse": {"value": 2.509538e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.HZLimEarlyMars": {"value": 2.738109e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.Instellation": {"value": -1.000000, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.sun.CriticalSemiMajorAxis": {"value": -1.000000, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.LXUVTot": {"value": 3.846000e+23, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.sun.LostEnergy": {"value": 2.568599e+28, "unit": u.Joule, "rtol": 1e-4}, + "log.final.sun.LostAngMom": {"value": 3.532078e+32, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, + "log.final.sun.EscapeVelocity": {"value": 1.146413e+06, "unit": u.m / u.sec, "rtol": 1e-4}, + "log.final.sun.Luminosity": {"value": 3.846000e+26, "unit": u.W, "rtol": 1e-4}, + "log.final.sun.LXUVStellar": {"value": 3.846000e+23, "unit": u.W, "rtol": 1e-4}, + "log.final.sun.Temperature": {"value": 5778.000000, "unit": u.K, "rtol": 1e-4}, + "log.final.sun.LXUVFrac": {"value": 0.001000, "rtol": 1e-4}, + "log.final.sun.RossbyNumber": {"value": 0.078260, "rtol": 1e-4}, + "log.final.sun.DRotPerDtStellar": {"value": 6.558557e-13, "rtol": 1e-4}, + "log.final.sun.WindTorque": {"value": 1.119248e+25, "rtol": 1e-4}, + "log.final.equatorial.Mass": {"value": 5.971546e+24, "unit": u.kg, "rtol": 1e-4}, + "log.final.equatorial.Obliquity": {"value": 0.410152, "unit": u.rad, "rtol": 1e-4}, + "log.final.equatorial.PrecA": {"value": 0.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.equatorial.Radius": {"value": 6.378100e+06, "unit": u.m, "rtol": 1e-4}, + "log.final.equatorial.RadGyra": {"value": 0.500000, "rtol": 1e-4}, + "log.final.equatorial.BodyType": {"value": 0.000000, "rtol": 1e-4}, + "log.final.equatorial.Density": {"value": 5494.449526, "unit": u.kg / u.m ** 3, "rtol": 1e-4}, + "log.final.equatorial.HZLimitDryRunaway": {"value": 1.117992e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.equatorial.HZLimRecVenus": {"value": 1.118929e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.equatorial.HZLimRunaway": {"value": 1.461108e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.equatorial.HZLimMoistGreenhouse": {"value": 1.480517e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.equatorial.HZLimMaxGreenhouse": {"value": 2.509538e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.equatorial.HZLimEarlyMars": {"value": 2.738109e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.equatorial.Instellation": {"value": 1367.566935, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.equatorial.Eccentricity": {"value": 0.000000, "rtol": 1e-4}, + "log.final.equatorial.MeanMotion": {"value": 1.990987e-07, "unit": 1 / u.sec, "rtol": 1e-4}, + "log.final.equatorial.OrbPeriod": {"value": 3.155815e+07, "unit": u.sec, "rtol": 1e-4}, + "log.final.equatorial.SemiMajorAxis": {"value": 1.495979e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.equatorial.COPP": {"value": 0.000000, "rtol": 1e-4}, + "log.final.equatorial.EscapeVelocity": {"value": 1.117931e+04, "unit": u.m / u.sec, "rtol": 1e-4}, + "log.final.equatorial.TGlobal": {"value": 16.638715, "unit": u.deg_C, "rtol": 1e-4}, + "log.final.equatorial.AlbedoGlobal": {"value": 0.322068, "rtol": 1e-4}, + "log.final.equatorial.FluxInGlobal": {"value": 237.998617, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.equatorial.FluxOutGlobal": {"value": 238.074915, "unit": u.W / u.m ** 2, "rtol": 1e-4}, + "log.final.equatorial.TotIceMass": {"value": 8.926243e+13, "unit": u.kg, "rtol": 1e-4}, + "log.final.equatorial.TotIceFlow": {"value": 0.000000, "unit": u.kg, "rtol": 1e-4}, + "log.final.equatorial.TotIceBalance": {"value": 9.114271e-08, "unit": u.kg, "rtol": 1e-4}, + "log.final.equatorial.SkipSeas": {"value": 0.000000, "rtol": 1e-4}, + "log.final.equatorial.AreaIceCov": {"value": 0.079735, "rtol": 1e-4}, + "log.final.equatorial.Latitude": {"value": 83.402352, "unit": u.deg, "rtol": 1e-4}, + "log.final.equatorial.TempLat": {"value": -10.684156, "unit": u.deg_C, "rtol": 1e-4}, + "log.final.equatorial.AlbedoLat": {"value": 0.600000, "rtol": 1e-4}, + "log.final.equatorial.AnnInsol": {"value": 176.064060, "unit": u.W / u.m ** 2, "rtol": 1e-4}, + "log.final.equatorial.FluxMerid": {"value": 1.632799, "unit": u.PW, "rtol": 1e-4}, + "log.final.equatorial.FluxIn": {"value": 70.432311, "unit": u.W / u.m ** 2, "rtol": 1e-4}, + "log.final.equatorial.FluxOut": {"value": 180.970113, "unit": u.W / u.m ** 2, "rtol": 1e-4}, + "log.final.equatorial.DivFlux": {"value": -110.460798, "unit": u.W / u.m ** 2, "rtol": 1e-4}, + "log.final.equatorial.IceMass": {"value": 176.402053, "rtol": 1e-4}, + "log.final.equatorial.IceHeight": {"value": 0.192432, "unit": u.m, "rtol": 1e-4}, + "log.final.equatorial.DIceMassDt": {"value": 5.589844e-06, "unit": u.m, "rtol": 1e-4}, + "log.final.equatorial.IceFlow": {"value": 0.000000, "unit": u.m / u.sec, "rtol": 1e-4}, + "log.final.equatorial.EnergyResL": {"value": -15.700340, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.equatorial.EnergyResW": {"value": 0.154552, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.equatorial.BedrockH": {"value": -1.046896e-05, "unit": u.m, "rtol": 1e-4}, + "log.final.equatorial.TempLandLat": {"value": 262.578337, "unit": u.sec, "rtol": 1e-4}, + "log.final.equatorial.TempWaterLat": {"value": 262.313192, "unit": u.sec, "rtol": 1e-4}, + "log.final.equatorial.AlbedoLandLat": {"value": 0.600000, "rtol": 1e-4}, + "log.final.equatorial.AlbedoWaterLat": {"value": 0.600000, "rtol": 1e-4}, + "log.final.equatorial.TempMinLat": {"value": -12.488965, "unit": u.deg_C, "rtol": 1e-4}, + "log.final.equatorial.TempMaxLat": {"value": -8.947302, "unit": u.deg_C, "rtol": 1e-4}, + "log.final.equatorial.Snowball": {"value": 0.000000, "rtol": 1e-4}, + "log.final.equatorial.PlanckBAvg": {"value": 2.090000, "rtol": 1e-4}, + "log.final.equatorial.IceAccum": {"value": 0.774581, "rtol": 1e-4}, + "log.final.equatorial.IceAblate": {"value": 0.000000, "rtol": 1e-4}, + "log.final.equatorial.TempMaxLand": {"value": 265.050604, "unit": u.sec, "rtol": 1e-4}, + "log.final.equatorial.TempMaxWater": {"value": 264.047615, "unit": u.sec, "rtol": 1e-4}, + "log.final.equatorial.PeakInsol": {"value": 541.684612, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.equatorial.IceCapNorthLand": {"value": 1.000000, "rtol": 1e-4}, + "log.final.equatorial.IceCapNorthSea": {"value": 1.000000, "rtol": 1e-4}, + "log.final.equatorial.IceCapSouthLand": {"value": 1.000000, "rtol": 1e-4}, + "log.final.equatorial.IceCapSouthSea": {"value": 1.000000, "rtol": 1e-4}, + "log.final.equatorial.IceBeltLand": {"value": 0.000000, "rtol": 1e-4}, + "log.final.equatorial.IceBeltSea": {"value": 0.000000, "rtol": 1e-4}, + "log.final.equatorial.SnowballLand": {"value": 0.000000, "rtol": 1e-4}, + "log.final.equatorial.SnowballSea": {"value": 0.000000, "rtol": 1e-4}, + "log.final.equatorial.IceFree": {"value": 0.000000, "rtol": 1e-4}, + "log.final.equatorial.IceCapNorthLatLand": {"value": 1.091712, "unit": u.rad, "rtol": 1e-4}, + "log.final.equatorial.IceCapNorthLatSea": {"value": 1.152808, "unit": u.rad, "rtol": 1e-4}, + "log.final.equatorial.IceCapSouthLatLand": {"value": -1.091712, "unit": u.rad, "rtol": 1e-4}, + "log.final.equatorial.IceCapSouthLatSea": {"value": -1.091712, "unit": u.rad, "rtol": 1e-4}, + "log.final.equatorial.IceBeltNorthLatLand": {"value": 100.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.equatorial.IceBeltNorthLatSea": {"value": 100.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.equatorial.IceBeltSouthLatLand": {"value": 100.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.equatorial.IceBeltSouthLatSea": {"value": 100.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.equatorial.LandFracGlobal": {"value": 0.289073, "rtol": 1e-4}, + } +) +class Test_SeasonalNC79Equatorial(Benchmark): + pass diff --git a/tests/Poise/SeasonalNC79Equatorial/vpl.in b/tests/Poise/SeasonalNC79Equatorial/vpl.in new file mode 100755 index 000000000..3afe6e0fa --- /dev/null +++ b/tests/Poise/SeasonalNC79Equatorial/vpl.in @@ -0,0 +1,15 @@ +sSystemName climateland +iVerbose 0 #how much do you want vplanet to yell at you? +iDigits 6 #how many digits do you want in your numbers? +bOverwrite 1 #overwrite old files +sUnitMass solar #mass unit used for input +sUnitLength au #length unit used for input +sUnitTime y #time unit +sUnitAngle d #angle unit +bDoLog 1 #create log file +saBodyFiles sun.in equatorial.in +bDoForward 1 #integrate forward in time +bVarDt 1 #use variable time stepping (not relevant to poise) +dEta 0.01 #how much to scale variable time step +dStopTime 1 #how long should the integration be +dOutputTime 1 #how much output you want diff --git a/tests/Poise/SeasonalNC79Polar/polar.in b/tests/Poise/SeasonalNC79Polar/polar.in new file mode 100755 index 000000000..02c19eead --- /dev/null +++ b/tests/Poise/SeasonalNC79Polar/polar.in @@ -0,0 +1,71 @@ +sName polar #name of planet +saModules poise #what vplanet modules you want to use +sGeography polar +dLandWaterLatitude 44.9 + +dMass 3.00316726e-06 #mass of planet +dRadius -1.00 #radius (not important right now) +dRotPeriod -1.00000 #rotation period (minus = days) +dObliquity 23.5 +dSemi 1 +dEcc 0.0 #eccentricity of orbit +dLongP 0 #pericenter, wrt Earth's position at spring equinox +# note that this is the typical value +180, +# since that one is solar position +dDynEllip 0.0 #shape of planet (0 = a sphere) +dPrecA 0.0 #orientation of spin axis + +#_______addition disorb/distrot parameters (leave these alone for now)__________________ +#dInc 5e-5 #inclination of orbit +#dLongA 348.73936 #orientation of orbital plane +#bGRCorr 0 #use GR correction (not important) +#bInvPlane 1 #convert to invariable plane coords +#bOverrideMaxEcc 1 #override max ecc halt (not recommended) +#dHaltMaxEcc 0.4 #eccentricity at which to halt simulation + +#_______poise parameters (have fun with these!)_________________________________________ +iLatCellNum 151 #number of latitude cells +sClimateModel sea #use seasonal or annual model +dTGlobalInit 14.85 #initial guess at average surface temp +iNumYears 4 #number of years (orbits) to run clim model +iNStepInYear 80 #number of steps to take in a "year" +#dSurfAlbedo 0.35 #average surface albedo (annual model only) + +#__ice params_________ +bIceSheets 1 #enable ice sheets +dInitIceLat 90. #how low do initial ice sheet extend? +dInitIceHeight 0. #height of initial ice sheets +dIceDepRate 2.25e-5 #rate of snow build up (when T < 0) +dIceAlbedo 0.6 #albedo of ice +iIceDt 1 #time step of ice-sheet model (orbits) +iReRunSeas 500 #how often to re-run seasonal model +bSeaIceModel 0 #use sea ice model (slow!) +bSkipSeasEnabled 0 #can skip seasonal if snowball state present + +#__heat diffusion______ +#bMEPDiff 1 #calculate diffusion using max entropy production +#bHadley 1 #mimic hadley heat diffusion +dDiffusion 0.58 #diffusion coefficient (fixed) +dNuLandWater 0.8 #Heat diffusion coefficient between Land and Water + +#__outgoing flux_______ +dPlanckA 203.3 #offset for OLR calculation (greenhouse) +dPlanckB 2.09 #slope of OLR calc (water vapor feedback) +bCalcAB 0 #calculate A & B from Kasting model fits +#dpCO2 0.00028 #partial pressure of co2 + +#__surface properties__ +dAlbedoLand 0.363 #albedo of land +dAlbedoWater 0.263 #albedo of water +dHeatCapLand 1.55e7 #land heat capacity +dHeatCapWater 4.428e6 #water heat capacity +dMixingDepth 70 #mixing depth of ocean + + +#________output options!_____________________________________________ +saOutputOrder Time PrecA -TGlobal AlbedoGlobal -FluxOutGlobal $ + -TotIceMass -TotIceFlow -TotIceBalance DeltaTime AreaIceCov Snowball Obliq $ + IceBeltLand IceBeltSea Ecce + +saGridOutput Time -Latitude -TempLat AlbedoLat -AnnInsol -FluxIn -FluxOut IceMass -IceHeight DIceMassDt $ + -IceFlow -BedrockH -TempMaxLat -TempMinLat -FluxMerid -DivFlux diff --git a/tests/Poise/SeasonalNC79Polar/sun.in b/tests/Poise/SeasonalNC79Polar/sun.in new file mode 100755 index 000000000..a1f97c666 --- /dev/null +++ b/tests/Poise/SeasonalNC79Polar/sun.in @@ -0,0 +1,9 @@ +# sun parameters +sName sun +dMass 1 +dSemi 0 +dEcc 0 +dRadius 0.00135 +dLuminosity -1 +sStellarModel none #sun does not change over time +saModules stellar #use stellar module (needed for luminosity) diff --git a/tests/Poise/SeasonalNC79Polar/test_SeasonalNC79Polar.py b/tests/Poise/SeasonalNC79Polar/test_SeasonalNC79Polar.py new file mode 100644 index 000000000..524eb1218 --- /dev/null +++ b/tests/Poise/SeasonalNC79Polar/test_SeasonalNC79Polar.py @@ -0,0 +1,228 @@ +from benchmark import Benchmark, benchmark +import astropy.units as u + +@benchmark( + { + "log.initial.system.Age": {"value": 0.000000, "unit": u.sec}, + "log.initial.system.Time": {"value": 0.000000, "unit": u.sec}, + "log.initial.system.TotAngMom": {"value": 1.501063e+42, "unit": (u.kg * u.m ** 2) / u.sec}, + "log.initial.system.TotEnergy": {"value": -7.839372e+41, "unit": u.Joule}, + "log.initial.system.PotEnergy": {"value": -7.839908e+41, "unit": u.Joule}, + "log.initial.system.KinEnergy": {"value": 5.361272e+37, "unit": u.Joule}, + "log.initial.system.DeltaTime": {"value": 0.000000, "unit": u.sec}, + "log.initial.sun.Mass": {"value": 1.988416e+30, "unit": u.kg}, + "log.initial.sun.Radius": {"value": 2.019571e+08, "unit": u.m}, + "log.initial.sun.RadGyra": {"value": 0.500000}, + "log.initial.sun.RotAngMom": {"value": 1.474456e+42, "unit": (u.kg * u.m ** 2) / u.sec}, + "log.initial.sun.RotVel": {"value": 1.468674e+04, "unit": u.m / u.sec}, + "log.initial.sun.BodyType": {"value": 0.000000}, + "log.initial.sun.RotRate": {"value": 7.272205e-05, "unit": 1 / u.sec}, + "log.initial.sun.RotPer": {"value": 8.640000e+04, "unit": u.sec}, + "log.initial.sun.Density": {"value": 5.762900e+04, "unit": u.kg / u.m ** 3}, + "log.initial.sun.HZLimitDryRunaway": {"value": 1.357831e+11, "unit": u.m}, + "log.initial.sun.HZLimRecVenus": {"value": 1.118929e+11, "unit": u.m}, + "log.initial.sun.HZLimRunaway": {"value": 1.461108e+11, "unit": u.m}, + "log.initial.sun.HZLimMoistGreenhouse": {"value": 1.480517e+11, "unit": u.m}, + "log.initial.sun.HZLimMaxGreenhouse": {"value": 2.509538e+11, "unit": u.m}, + "log.initial.sun.HZLimEarlyMars": {"value": 2.738109e+11, "unit": u.m}, + "log.initial.sun.Instellation": {"value": -1.000000, "unit": u.kg / u.sec ** 3}, + "log.initial.sun.CriticalSemiMajorAxis": {"value": -1.000000, "unit": u.m}, + "log.initial.sun.LXUVTot": {"value": 3.846000e+23, "unit": u.kg / u.sec ** 3}, + "log.initial.sun.LostEnergy": {"value": 5.562685e-309, "unit": u.Joule}, + "log.initial.sun.LostAngMom": {"value": 5.562685e-309, "unit": (u.kg * u.m ** 2) / u.sec}, + "log.initial.sun.EscapeVelocity": {"value": 1.146413e+06, "unit": u.m / u.sec}, + "log.initial.sun.Luminosity": {"value": 3.846000e+26, "unit": u.W}, + "log.initial.sun.LXUVStellar": {"value": 3.846000e+23, "unit": u.W}, + "log.initial.sun.Temperature": {"value": 5778.000000, "unit": u.K}, + "log.initial.sun.LXUVFrac": {"value": 0.001000}, + "log.initial.sun.RossbyNumber": {"value": 0.078260}, + "log.initial.sun.DRotPerDtStellar": {"value": 6.558557e-13}, + "log.initial.sun.WindTorque": {"value": 1.119248e+25}, + "log.initial.polar.Mass": {"value": 5.971546e+24, "unit": u.kg}, + "log.initial.polar.Obliquity": {"value": 0.410152, "unit": u.rad}, + "log.initial.polar.PrecA": {"value": 0.000000, "unit": u.rad}, + "log.initial.polar.Radius": {"value": 6.378100e+06, "unit": u.m}, + "log.initial.polar.RadGyra": {"value": 0.500000}, + "log.initial.polar.BodyType": {"value": 0.000000}, + "log.initial.polar.Density": {"value": 5494.449526, "unit": u.kg / u.m ** 3}, + "log.initial.polar.HZLimitDryRunaway": {"value": 1.115053e+11, "unit": u.m}, + "log.initial.polar.HZLimRecVenus": {"value": 1.118929e+11, "unit": u.m}, + "log.initial.polar.HZLimRunaway": {"value": 1.461108e+11, "unit": u.m}, + "log.initial.polar.HZLimMoistGreenhouse": {"value": 1.480517e+11, "unit": u.m}, + "log.initial.polar.HZLimMaxGreenhouse": {"value": 2.509538e+11, "unit": u.m}, + "log.initial.polar.HZLimEarlyMars": {"value": 2.738109e+11, "unit": u.m}, + "log.initial.polar.Instellation": {"value": 1367.566935, "unit": u.kg / u.sec ** 3}, + "log.initial.polar.Eccentricity": {"value": 0.000000}, + "log.initial.polar.MeanMotion": {"value": 1.990987e-07, "unit": 1 / u.sec}, + "log.initial.polar.OrbPeriod": {"value": 3.155815e+07, "unit": u.sec}, + "log.initial.polar.SemiMajorAxis": {"value": 1.495979e+11, "unit": u.m}, + "log.initial.polar.COPP": {"value": 0.000000}, + "log.initial.polar.EscapeVelocity": {"value": 1.117931e+04, "unit": u.m / u.sec}, + "log.initial.polar.TGlobal": {"value": 18.638682, "unit": u.deg_C}, + "log.initial.polar.AlbedoGlobal": {"value": 0.325628}, + "log.initial.polar.FluxInGlobal": {"value": 241.918401, "unit": u.kg / u.sec ** 3}, + "log.initial.polar.FluxOutGlobal": {"value": 242.254845, "unit": u.W / u.m ** 2}, + "log.initial.polar.TotIceMass": {"value": 0.000000, "unit": u.kg}, + "log.initial.polar.TotIceFlow": {"value": 0.000000, "unit": u.kg}, + "log.initial.polar.TotIceBalance": {"value": 0.000000, "unit": u.kg}, + "log.initial.polar.SkipSeas": {"value": 0.000000}, + "log.initial.polar.AreaIceCov": {"value": 0.000000}, + "log.initial.polar.Latitude": {"value": -83.402352, "unit": u.deg}, + "log.initial.polar.TempLat": {"value": -12.816826, "unit": u.deg_C}, + "log.initial.polar.AlbedoLat": {"value": 0.574553}, + "log.initial.polar.AnnInsol": {"value": 176.064060, "unit": u.W / u.m ** 2}, + "log.initial.polar.FluxMerid": {"value": -1.477182, "unit": u.PW}, + "log.initial.polar.FluxIn": {"value": 75.552308, "unit": u.W / u.m ** 2}, + "log.initial.polar.FluxOut": {"value": 176.512834, "unit": u.W / u.m ** 2}, + "log.initial.polar.DivFlux": {"value": -99.933152, "unit": u.W / u.m ** 2}, + "log.initial.polar.IceMass": {"value": 0.000000}, + "log.initial.polar.IceHeight": {"value": 0.000000, "unit": u.m}, + "log.initial.polar.DIceMassDt": {"value": 0.000000, "unit": u.m}, + "log.initial.polar.IceFlow": {"value": 0.000000, "unit": u.m / u.sec}, + "log.initial.polar.EnergyResL": {"value": -8.639695, "unit": u.kg / u.sec ** 3}, + "log.initial.polar.EnergyResW": {"value": -16.409422, "unit": u.kg / u.sec ** 3}, + "log.initial.polar.BedrockH": {"value": 0.000000, "unit": u.m}, + "log.initial.polar.TempLandLat": {"value": 260.182790, "unit": u.sec}, + "log.initial.polar.TempWaterLat": {"value": 260.221189, "unit": u.sec}, + "log.initial.polar.AlbedoLandLat": {"value": 0.574931}, + "log.initial.polar.AlbedoWaterLat": {"value": 0.537138}, + "log.initial.polar.TempMinLat": {"value": -29.929775, "unit": u.deg_C}, + "log.initial.polar.TempMaxLat": {"value": 10.778913, "unit": u.deg_C}, + "log.initial.polar.Snowball": {"value": 0.000000}, + "log.initial.polar.PlanckBAvg": {"value": 2.090000}, + "log.initial.polar.IceAccum": {"value": 0.561571}, + "log.initial.polar.IceAblate": {"value": -0.358381}, + "log.initial.polar.TempMaxLand": {"value": 283.860961, "unit": u.sec}, + "log.initial.polar.TempMaxWater": {"value": 277.769388, "unit": u.sec}, + "log.initial.polar.PeakInsol": {"value": 541.704677, "unit": u.kg / u.sec ** 3}, + "log.initial.polar.IceCapNorthLand": {"value": 0.000000}, + "log.initial.polar.IceCapNorthSea": {"value": 0.000000}, + "log.initial.polar.IceCapSouthLand": {"value": 0.000000}, + "log.initial.polar.IceCapSouthSea": {"value": 0.000000}, + "log.initial.polar.IceBeltLand": {"value": 0.000000}, + "log.initial.polar.IceBeltSea": {"value": 0.000000}, + "log.initial.polar.SnowballLand": {"value": 0.000000}, + "log.initial.polar.SnowballSea": {"value": 0.000000}, + "log.initial.polar.IceFree": {"value": 1.000000}, + "log.initial.polar.IceCapNorthLatLand": {"value": 100.000000, "unit": u.rad}, + "log.initial.polar.IceCapNorthLatSea": {"value": 100.000000, "unit": u.rad}, + "log.initial.polar.IceCapSouthLatLand": {"value": 100.000000, "unit": u.rad}, + "log.initial.polar.IceCapSouthLatSea": {"value": 100.000000, "unit": u.rad}, + "log.initial.polar.IceBeltNorthLatLand": {"value": 100.000000, "unit": u.rad}, + "log.initial.polar.IceBeltNorthLatSea": {"value": 100.000000, "unit": u.rad}, + "log.initial.polar.IceBeltSouthLatLand": {"value": 100.000000, "unit": u.rad}, + "log.initial.polar.IceBeltSouthLatSea": {"value": 100.000000, "unit": u.rad}, + "log.initial.polar.LandFracGlobal": {"value": 0.295563}, + "log.final.system.Age": {"value": 3.155760e+07, "unit": u.sec, "rtol": 1e-4}, + "log.final.system.Time": {"value": 3.155760e+07, "unit": u.sec, "rtol": 1e-4}, + "log.final.system.TotAngMom": {"value": 1.501063e+42, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, + "log.final.system.TotEnergy": {"value": -7.839372e+41, "unit": u.Joule, "rtol": 1e-4}, + "log.final.system.PotEnergy": {"value": -7.839908e+41, "unit": u.Joule, "rtol": 1e-4}, + "log.final.system.KinEnergy": {"value": 5.361272e+37, "unit": u.Joule, "rtol": 1e-4}, + "log.final.sun.Mass": {"value": 1.988416e+30, "unit": u.kg, "rtol": 1e-4}, + "log.final.sun.Radius": {"value": 2.019571e+08, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.RadGyra": {"value": 0.500000, "rtol": 1e-4}, + "log.final.sun.RotAngMom": {"value": 1.474456e+42, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, + "log.final.sun.RotVel": {"value": 1.468674e+04, "unit": u.m / u.sec, "rtol": 1e-4}, + "log.final.sun.BodyType": {"value": 0.000000, "rtol": 1e-4}, + "log.final.sun.RotRate": {"value": 7.272205e-05, "unit": 1 / u.sec, "rtol": 1e-4}, + "log.final.sun.RotPer": {"value": 8.640000e+04, "unit": u.sec, "rtol": 1e-4}, + "log.final.sun.Density": {"value": 5.762900e+04, "unit": u.kg / u.m ** 3, "rtol": 1e-4}, + "log.final.sun.HZLimitDryRunaway": {"value": 1.357831e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.HZLimRecVenus": {"value": 1.118929e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.HZLimRunaway": {"value": 1.461108e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.HZLimMoistGreenhouse": {"value": 1.480517e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.HZLimMaxGreenhouse": {"value": 2.509538e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.HZLimEarlyMars": {"value": 2.738109e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.Instellation": {"value": -1.000000, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.sun.CriticalSemiMajorAxis": {"value": -1.000000, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.LXUVTot": {"value": 3.846000e+23, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.sun.LostEnergy": {"value": 2.568599e+28, "unit": u.Joule, "rtol": 1e-4}, + "log.final.sun.LostAngMom": {"value": 3.532078e+32, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, + "log.final.sun.EscapeVelocity": {"value": 1.146413e+06, "unit": u.m / u.sec, "rtol": 1e-4}, + "log.final.sun.Luminosity": {"value": 3.846000e+26, "unit": u.W, "rtol": 1e-4}, + "log.final.sun.LXUVStellar": {"value": 3.846000e+23, "unit": u.W, "rtol": 1e-4}, + "log.final.sun.Temperature": {"value": 5778.000000, "unit": u.K, "rtol": 1e-4}, + "log.final.sun.LXUVFrac": {"value": 0.001000, "rtol": 1e-4}, + "log.final.sun.RossbyNumber": {"value": 0.078260, "rtol": 1e-4}, + "log.final.sun.DRotPerDtStellar": {"value": 6.558557e-13, "rtol": 1e-4}, + "log.final.sun.WindTorque": {"value": 1.119248e+25, "rtol": 1e-4}, + "log.final.polar.Mass": {"value": 5.971546e+24, "unit": u.kg, "rtol": 1e-4}, + "log.final.polar.Obliquity": {"value": 0.410152, "unit": u.rad, "rtol": 1e-4}, + "log.final.polar.PrecA": {"value": 0.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.polar.Radius": {"value": 6.378100e+06, "unit": u.m, "rtol": 1e-4}, + "log.final.polar.RadGyra": {"value": 0.500000, "rtol": 1e-4}, + "log.final.polar.BodyType": {"value": 0.000000, "rtol": 1e-4}, + "log.final.polar.Density": {"value": 5494.449526, "unit": u.kg / u.m ** 3, "rtol": 1e-4}, + "log.final.polar.HZLimitDryRunaway": {"value": 1.115053e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.polar.HZLimRecVenus": {"value": 1.118929e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.polar.HZLimRunaway": {"value": 1.461108e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.polar.HZLimMoistGreenhouse": {"value": 1.480517e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.polar.HZLimMaxGreenhouse": {"value": 2.509538e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.polar.HZLimEarlyMars": {"value": 2.738109e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.polar.Instellation": {"value": 1367.566935, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.polar.Eccentricity": {"value": 0.000000, "rtol": 1e-4}, + "log.final.polar.MeanMotion": {"value": 1.990987e-07, "unit": 1 / u.sec, "rtol": 1e-4}, + "log.final.polar.OrbPeriod": {"value": 3.155815e+07, "unit": u.sec, "rtol": 1e-4}, + "log.final.polar.SemiMajorAxis": {"value": 1.495979e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.polar.COPP": {"value": 0.000000, "rtol": 1e-4}, + "log.final.polar.EscapeVelocity": {"value": 1.117931e+04, "unit": u.m / u.sec, "rtol": 1e-4}, + "log.final.polar.TGlobal": {"value": 18.638845, "unit": u.deg_C, "rtol": 1e-4}, + "log.final.polar.AlbedoGlobal": {"value": 0.325628, "rtol": 1e-4}, + "log.final.polar.FluxInGlobal": {"value": 241.918401, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.polar.FluxOutGlobal": {"value": 242.255187, "unit": u.W / u.m ** 2, "rtol": 1e-4}, + "log.final.polar.TotIceMass": {"value": 0.000000, "unit": u.kg, "rtol": 1e-4}, + "log.final.polar.TotIceFlow": {"value": 0.000000, "unit": u.kg, "rtol": 1e-4}, + "log.final.polar.TotIceBalance": {"value": 0.000000, "unit": u.kg, "rtol": 1e-4}, + "log.final.polar.SkipSeas": {"value": 0.000000, "rtol": 1e-4}, + "log.final.polar.AreaIceCov": {"value": 0.000000, "rtol": 1e-4}, + "log.final.polar.Latitude": {"value": 83.402352, "unit": u.deg, "rtol": 1e-4}, + "log.final.polar.TempLat": {"value": -12.829124, "unit": u.deg_C, "rtol": 1e-4}, + "log.final.polar.AlbedoLat": {"value": 0.574549, "rtol": 1e-4}, + "log.final.polar.AnnInsol": {"value": 176.064060, "unit": u.W / u.m ** 2, "rtol": 1e-4}, + "log.final.polar.FluxMerid": {"value": 1.474909, "unit": u.PW, "rtol": 1e-4}, + "log.final.polar.FluxIn": {"value": 75.557721, "unit": u.W / u.m ** 2, "rtol": 1e-4}, + "log.final.polar.FluxOut": {"value": 176.487131, "unit": u.W / u.m ** 2, "rtol": 1e-4}, + "log.final.polar.DivFlux": {"value": -99.779362, "unit": u.W / u.m ** 2, "rtol": 1e-4}, + "log.final.polar.IceMass": {"value": 0.000000, "rtol": 1e-4}, + "log.final.polar.IceHeight": {"value": 0.000000, "unit": u.m, "rtol": 1e-4}, + "log.final.polar.DIceMassDt": {"value": 0.000000, "unit": u.m, "rtol": 1e-4}, + "log.final.polar.IceFlow": {"value": 0.000000, "unit": u.m / u.sec, "rtol": 1e-4}, + "log.final.polar.EnergyResL": {"value": -0.554280, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.polar.EnergyResW": {"value": 15.300611, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.polar.BedrockH": {"value": 0.000000, "unit": u.m, "rtol": 1e-4}, + "log.final.polar.TempLandLat": {"value": 260.170525, "unit": u.sec, "rtol": 1e-4}, + "log.final.polar.TempWaterLat": {"value": 260.205647, "unit": u.sec, "rtol": 1e-4}, + "log.final.polar.AlbedoLandLat": {"value": 0.574926, "rtol": 1e-4}, + "log.final.polar.AlbedoWaterLat": {"value": 0.537137, "rtol": 1e-4}, + "log.final.polar.TempMinLat": {"value": -29.935349, "unit": u.deg_C, "rtol": 1e-4}, + "log.final.polar.TempMaxLat": {"value": 10.711221, "unit": u.deg_C, "rtol": 1e-4}, + "log.final.polar.Snowball": {"value": 0.000000, "rtol": 1e-4}, + "log.final.polar.PlanckBAvg": {"value": 2.090000, "rtol": 1e-4}, + "log.final.polar.IceAccum": {"value": 0.561571, "rtol": 1e-4}, + "log.final.polar.IceAblate": {"value": -0.429123, "rtol": 1e-4}, + "log.final.polar.TempMaxLand": {"value": 283.792990, "unit": u.sec, "rtol": 1e-4}, + "log.final.polar.TempMaxWater": {"value": 277.742487, "unit": u.sec, "rtol": 1e-4}, + "log.final.polar.PeakInsol": {"value": 541.684612, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.polar.IceCapNorthLand": {"value": 0.000000, "rtol": 1e-4}, + "log.final.polar.IceCapNorthSea": {"value": 0.000000, "rtol": 1e-4}, + "log.final.polar.IceCapSouthLand": {"value": 0.000000, "rtol": 1e-4}, + "log.final.polar.IceCapSouthSea": {"value": 0.000000, "rtol": 1e-4}, + "log.final.polar.IceBeltLand": {"value": 0.000000, "rtol": 1e-4}, + "log.final.polar.IceBeltSea": {"value": 0.000000, "rtol": 1e-4}, + "log.final.polar.SnowballLand": {"value": 0.000000, "rtol": 1e-4}, + "log.final.polar.SnowballSea": {"value": 0.000000, "rtol": 1e-4}, + "log.final.polar.IceFree": {"value": 1.000000, "rtol": 1e-4}, + "log.final.polar.IceCapNorthLatLand": {"value": 100.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.polar.IceCapNorthLatSea": {"value": 100.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.polar.IceCapSouthLatLand": {"value": 100.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.polar.IceCapSouthLatSea": {"value": 100.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.polar.IceBeltNorthLatLand": {"value": 100.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.polar.IceBeltNorthLatSea": {"value": 100.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.polar.IceBeltSouthLatLand": {"value": 100.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.polar.IceBeltSouthLatSea": {"value": 100.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.polar.LandFracGlobal": {"value": 0.295563, "rtol": 1e-4}, + } +) +class Test_SeasonalNC79Polar(Benchmark): + pass diff --git a/tests/Poise/SeasonalNC79Polar/vpl.in b/tests/Poise/SeasonalNC79Polar/vpl.in new file mode 100755 index 000000000..fc34dd9c0 --- /dev/null +++ b/tests/Poise/SeasonalNC79Polar/vpl.in @@ -0,0 +1,15 @@ +sSystemName climateland +iVerbose 0 #how much do you want vplanet to yell at you? +iDigits 6 #how many digits do you want in your numbers? +bOverwrite 1 #overwrite old files +sUnitMass solar #mass unit used for input +sUnitLength au #length unit used for input +sUnitTime y #time unit +sUnitAngle d #angle unit +bDoLog 1 #create log file +saBodyFiles sun.in polar.in +bDoForward 1 #integrate forward in time +bVarDt 1 #use variable time stepping (not relevant to poise) +dEta 0.01 #how much to scale variable time step +dStopTime 1 #how long should the integration be +dOutputTime 1 #how much output you want diff --git a/tests/Poise/SeasonalNC79Random/random.in b/tests/Poise/SeasonalNC79Random/random.in new file mode 100755 index 000000000..3c28f276a --- /dev/null +++ b/tests/Poise/SeasonalNC79Random/random.in @@ -0,0 +1,73 @@ +sName random #name of planet +saModules poise #what vplanet modules you want to use +sGeography random +iRandSeed 4.567e8 + +dMass 3.00316726e-06 #mass of planet +dRadius -1.00 #radius (not important right now) +dRotPeriod -1.00000 #rotation period (minus = days) +dObliquity 23.5 +dSemi 1 +dEcc 0.0 #eccentricity of orbit +dLongP 0 #pericenter, wrt Earth's position at spring equinox +# note that this is the typical value +180, +# since that one is solar position +dDynEllip 0.0 #shape of planet (0 = a sphere) +dPrecA 0.0 #orientation of spin axis + +#_______addition disorb/distrot parameters (leave these alone for now)__________________ +#dInc 5e-5 #inclination of orbit +#dLongA 348.73936 #orientation of orbital plane +#bGRCorr 0 #use GR correction (not important) +#bInvPlane 1 #convert to invariable plane coords +#bOverrideMaxEcc 1 #override max ecc halt (not recommended) +#dHaltMaxEcc 0.4 #eccentricity at which to halt simulation + +#_______poise parameters (have fun with these!)_________________________________________ +iLatCellNum 151 #number of latitude cells +sClimateModel sea #use seasonal or annual model +dTGlobalInit 14.85 #initial guess at average surface temp +iNumYears 4 #number of years (orbits) to run clim model +iNStepInYear 80 #number of steps to take in a "year" +#dSurfAlbedo 0.35 #average surface albedo (annual model only) +dLandFracMean 0.292 +dLandFracAmp 0.2 + +#__ice params_________ +bIceSheets 1 #enable ice sheets +dInitIceLat 90. #how low do initial ice sheet extend? +dInitIceHeight 0. #height of initial ice sheets +dIceDepRate 2.25e-5 #rate of snow build up (when T < 0) +dIceAlbedo 0.6 #albedo of ice +iIceDt 1 #time step of ice-sheet model (orbits) +iReRunSeas 500 #how often to re-run seasonal model +bSeaIceModel 0 #use sea ice model (slow!) +bSkipSeasEnabled 0 #can skip seasonal if snowball state present + +#__heat diffusion______ +#bMEPDiff 1 #calculate diffusion using max entropy production +#bHadley 1 #mimic hadley heat diffusion +dDiffusion 0.58 #diffusion coefficient (fixed) +dNuLandWater 0.8 #Heat diffusion coefficient between Land and Water + +#__outgoing flux_______ +dPlanckA 203.3 #offset for OLR calculation (greenhouse) +dPlanckB 2.09 #slope of OLR calc (water vapor feedback) +bCalcAB 0 #calculate A & B from Kasting model fits +#dpCO2 0.00028 #partial pressure of co2 + +#__surface properties__ +dAlbedoLand 0.363 #albedo of land +dAlbedoWater 0.263 #albedo of water +dHeatCapLand 1.55e7 #land heat capacity +dHeatCapWater 4.428e6 #water heat capacity +dMixingDepth 70 #mixing depth of ocean + + +#________output options!_____________________________________________ +saOutputOrder Time PrecA -TGlobal AlbedoGlobal -FluxOutGlobal $ + -TotIceMass -TotIceFlow -TotIceBalance DeltaTime AreaIceCov Snowball Obliq $ + IceBeltLand IceBeltSea Ecce + +saGridOutput Time -Latitude -TempLat AlbedoLat -AnnInsol -FluxIn -FluxOut IceMass -IceHeight DIceMassDt $ + -IceFlow -BedrockH -TempMaxLat -TempMinLat -FluxMerid -DivFlux diff --git a/tests/Poise/SeasonalNC79Random/sun.in b/tests/Poise/SeasonalNC79Random/sun.in new file mode 100755 index 000000000..a1f97c666 --- /dev/null +++ b/tests/Poise/SeasonalNC79Random/sun.in @@ -0,0 +1,9 @@ +# sun parameters +sName sun +dMass 1 +dSemi 0 +dEcc 0 +dRadius 0.00135 +dLuminosity -1 +sStellarModel none #sun does not change over time +saModules stellar #use stellar module (needed for luminosity) diff --git a/tests/Poise/SeasonalNC79Random/test_SeasonalNC79Random.py b/tests/Poise/SeasonalNC79Random/test_SeasonalNC79Random.py new file mode 100644 index 000000000..c7993c336 --- /dev/null +++ b/tests/Poise/SeasonalNC79Random/test_SeasonalNC79Random.py @@ -0,0 +1,228 @@ +from benchmark import Benchmark, benchmark +import astropy.units as u + +@benchmark( + { + "log.initial.system.Age": {"value": 0.000000, "unit": u.sec}, + "log.initial.system.Time": {"value": 0.000000, "unit": u.sec}, + "log.initial.system.TotAngMom": {"value": 1.501063e+42, "unit": (u.kg * u.m ** 2) / u.sec}, + "log.initial.system.TotEnergy": {"value": -7.839372e+41, "unit": u.Joule}, + "log.initial.system.PotEnergy": {"value": -7.839908e+41, "unit": u.Joule}, + "log.initial.system.KinEnergy": {"value": 5.361272e+37, "unit": u.Joule}, + "log.initial.system.DeltaTime": {"value": 0.000000, "unit": u.sec}, + "log.initial.sun.Mass": {"value": 1.988416e+30, "unit": u.kg}, + "log.initial.sun.Radius": {"value": 2.019571e+08, "unit": u.m}, + "log.initial.sun.RadGyra": {"value": 0.500000}, + "log.initial.sun.RotAngMom": {"value": 1.474456e+42, "unit": (u.kg * u.m ** 2) / u.sec}, + "log.initial.sun.RotVel": {"value": 1.468674e+04, "unit": u.m / u.sec}, + "log.initial.sun.BodyType": {"value": 0.000000}, + "log.initial.sun.RotRate": {"value": 7.272205e-05, "unit": 1 / u.sec}, + "log.initial.sun.RotPer": {"value": 8.640000e+04, "unit": u.sec}, + "log.initial.sun.Density": {"value": 5.762900e+04, "unit": u.kg / u.m ** 3}, + "log.initial.sun.HZLimitDryRunaway": {"value": 1.357831e+11, "unit": u.m}, + "log.initial.sun.HZLimRecVenus": {"value": 1.118929e+11, "unit": u.m}, + "log.initial.sun.HZLimRunaway": {"value": 1.461108e+11, "unit": u.m}, + "log.initial.sun.HZLimMoistGreenhouse": {"value": 1.480517e+11, "unit": u.m}, + "log.initial.sun.HZLimMaxGreenhouse": {"value": 2.509538e+11, "unit": u.m}, + "log.initial.sun.HZLimEarlyMars": {"value": 2.738109e+11, "unit": u.m}, + "log.initial.sun.Instellation": {"value": -1.000000, "unit": u.kg / u.sec ** 3}, + "log.initial.sun.CriticalSemiMajorAxis": {"value": -1.000000, "unit": u.m}, + "log.initial.sun.LXUVTot": {"value": 3.846000e+23, "unit": u.kg / u.sec ** 3}, + "log.initial.sun.LostEnergy": {"value": 5.562685e-309, "unit": u.Joule}, + "log.initial.sun.LostAngMom": {"value": 5.562685e-309, "unit": (u.kg * u.m ** 2) / u.sec}, + "log.initial.sun.EscapeVelocity": {"value": 1.146413e+06, "unit": u.m / u.sec}, + "log.initial.sun.Luminosity": {"value": 3.846000e+26, "unit": u.W}, + "log.initial.sun.LXUVStellar": {"value": 3.846000e+23, "unit": u.W}, + "log.initial.sun.Temperature": {"value": 5778.000000, "unit": u.K}, + "log.initial.sun.LXUVFrac": {"value": 0.001000}, + "log.initial.sun.RossbyNumber": {"value": 0.078260}, + "log.initial.sun.DRotPerDtStellar": {"value": 6.558557e-13}, + "log.initial.sun.WindTorque": {"value": 1.119248e+25}, + "log.initial.random.Mass": {"value": 5.971546e+24, "unit": u.kg}, + "log.initial.random.Obliquity": {"value": 0.410152, "unit": u.rad}, + "log.initial.random.PrecA": {"value": 0.000000, "unit": u.rad}, + "log.initial.random.Radius": {"value": 6.378100e+06, "unit": u.m}, + "log.initial.random.RadGyra": {"value": 0.500000}, + "log.initial.random.BodyType": {"value": 0.000000}, + "log.initial.random.Density": {"value": 5494.449526, "unit": u.kg / u.m ** 3}, + "log.initial.random.HZLimitDryRunaway": {"value": 1.112428e+11, "unit": u.m}, + "log.initial.random.HZLimRecVenus": {"value": 1.118929e+11, "unit": u.m}, + "log.initial.random.HZLimRunaway": {"value": 1.461108e+11, "unit": u.m}, + "log.initial.random.HZLimMoistGreenhouse": {"value": 1.480517e+11, "unit": u.m}, + "log.initial.random.HZLimMaxGreenhouse": {"value": 2.509538e+11, "unit": u.m}, + "log.initial.random.HZLimEarlyMars": {"value": 2.738109e+11, "unit": u.m}, + "log.initial.random.Instellation": {"value": 1367.566935, "unit": u.kg / u.sec ** 3}, + "log.initial.random.Eccentricity": {"value": 0.000000}, + "log.initial.random.MeanMotion": {"value": 1.990987e-07, "unit": 1 / u.sec}, + "log.initial.random.OrbPeriod": {"value": 3.155815e+07, "unit": u.sec}, + "log.initial.random.SemiMajorAxis": {"value": 1.495979e+11, "unit": u.m}, + "log.initial.random.COPP": {"value": 0.000000}, + "log.initial.random.EscapeVelocity": {"value": 1.117931e+04, "unit": u.m / u.sec}, + "log.initial.random.TGlobal": {"value": 17.283545, "unit": u.deg_C}, + "log.initial.random.AlbedoGlobal": {"value": 0.328800}, + "log.initial.random.FluxInGlobal": {"value": 238.235013, "unit": u.kg / u.sec ** 3}, + "log.initial.random.FluxOutGlobal": {"value": 239.422609, "unit": u.W / u.m ** 2}, + "log.initial.random.TotIceMass": {"value": 0.000000, "unit": u.kg}, + "log.initial.random.TotIceFlow": {"value": 0.000000, "unit": u.kg}, + "log.initial.random.TotIceBalance": {"value": 0.000000, "unit": u.kg}, + "log.initial.random.SkipSeas": {"value": 0.000000}, + "log.initial.random.AreaIceCov": {"value": 0.055747}, + "log.initial.random.Latitude": {"value": -83.402352, "unit": u.deg}, + "log.initial.random.TempLat": {"value": -13.024810, "unit": u.deg_C}, + "log.initial.random.AlbedoLat": {"value": 0.596005}, + "log.initial.random.AnnInsol": {"value": 176.064060, "unit": u.W / u.m ** 2}, + "log.initial.random.FluxMerid": {"value": -2.255041, "unit": u.PW}, + "log.initial.random.FluxIn": {"value": 71.047086, "unit": u.W / u.m ** 2}, + "log.initial.random.FluxOut": {"value": 176.078148, "unit": u.W / u.m ** 2}, + "log.initial.random.DivFlux": {"value": -152.556245, "unit": u.W / u.m ** 2}, + "log.initial.random.IceMass": {"value": 0.000000}, + "log.initial.random.IceHeight": {"value": 0.000000, "unit": u.m}, + "log.initial.random.DIceMassDt": {"value": 0.000000, "unit": u.m}, + "log.initial.random.IceFlow": {"value": 0.000000, "unit": u.m / u.sec}, + "log.initial.random.EnergyResL": {"value": -0.715364, "unit": u.kg / u.sec ** 3}, + "log.initial.random.EnergyResW": {"value": 0.296959, "unit": u.kg / u.sec ** 3}, + "log.initial.random.BedrockH": {"value": 0.000000, "unit": u.m}, + "log.initial.random.TempLandLat": {"value": 258.698815, "unit": u.sec}, + "log.initial.random.TempWaterLat": {"value": 261.176041, "unit": u.sec}, + "log.initial.random.AlbedoLandLat": {"value": 0.591759}, + "log.initial.random.AlbedoWaterLat": {"value": 0.600000}, + "log.initial.random.TempMinLat": {"value": -21.298326, "unit": u.deg_C}, + "log.initial.random.TempMaxLat": {"value": -3.277995, "unit": u.deg_C}, + "log.initial.random.Snowball": {"value": 0.000000}, + "log.initial.random.PlanckBAvg": {"value": 2.090000}, + "log.initial.random.IceAccum": {"value": 0.597880}, + "log.initial.random.IceAblate": {"value": -0.429577}, + "log.initial.random.TempMaxLand": {"value": 276.707036, "unit": u.sec}, + "log.initial.random.TempMaxWater": {"value": 263.187326, "unit": u.sec}, + "log.initial.random.PeakInsol": {"value": 541.704677, "unit": u.kg / u.sec ** 3}, + "log.initial.random.IceCapNorthLand": {"value": 0.000000}, + "log.initial.random.IceCapNorthSea": {"value": 1.000000}, + "log.initial.random.IceCapSouthLand": {"value": 0.000000}, + "log.initial.random.IceCapSouthSea": {"value": 1.000000}, + "log.initial.random.IceBeltLand": {"value": 0.000000}, + "log.initial.random.IceBeltSea": {"value": 0.000000}, + "log.initial.random.SnowballLand": {"value": 0.000000}, + "log.initial.random.SnowballSea": {"value": 0.000000}, + "log.initial.random.IceFree": {"value": 0.000000}, + "log.initial.random.IceCapNorthLatLand": {"value": 100.000000, "unit": u.rad}, + "log.initial.random.IceCapNorthLatSea": {"value": 1.152808, "unit": u.rad}, + "log.initial.random.IceCapSouthLatLand": {"value": 100.000000, "unit": u.rad}, + "log.initial.random.IceCapSouthLatSea": {"value": 100.000000, "unit": u.rad}, + "log.initial.random.IceBeltNorthLatLand": {"value": 100.000000, "unit": u.rad}, + "log.initial.random.IceBeltNorthLatSea": {"value": 100.000000, "unit": u.rad}, + "log.initial.random.IceBeltSouthLatLand": {"value": 100.000000, "unit": u.rad}, + "log.initial.random.IceBeltSouthLatSea": {"value": 100.000000, "unit": u.rad}, + "log.initial.random.LandFracGlobal": {"value": 0.292000}, + "log.final.system.Age": {"value": 3.155760e+07, "unit": u.sec, "rtol": 1e-4}, + "log.final.system.Time": {"value": 3.155760e+07, "unit": u.sec, "rtol": 1e-4}, + "log.final.system.TotAngMom": {"value": 1.501063e+42, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, + "log.final.system.TotEnergy": {"value": -7.839372e+41, "unit": u.Joule, "rtol": 1e-4}, + "log.final.system.PotEnergy": {"value": -7.839908e+41, "unit": u.Joule, "rtol": 1e-4}, + "log.final.system.KinEnergy": {"value": 5.361272e+37, "unit": u.Joule, "rtol": 1e-4}, + "log.final.sun.Mass": {"value": 1.988416e+30, "unit": u.kg, "rtol": 1e-4}, + "log.final.sun.Radius": {"value": 2.019571e+08, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.RadGyra": {"value": 0.500000, "rtol": 1e-4}, + "log.final.sun.RotAngMom": {"value": 1.474456e+42, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, + "log.final.sun.RotVel": {"value": 1.468674e+04, "unit": u.m / u.sec, "rtol": 1e-4}, + "log.final.sun.BodyType": {"value": 0.000000, "rtol": 1e-4}, + "log.final.sun.RotRate": {"value": 7.272205e-05, "unit": 1 / u.sec, "rtol": 1e-4}, + "log.final.sun.RotPer": {"value": 8.640000e+04, "unit": u.sec, "rtol": 1e-4}, + "log.final.sun.Density": {"value": 5.762900e+04, "unit": u.kg / u.m ** 3, "rtol": 1e-4}, + "log.final.sun.HZLimitDryRunaway": {"value": 1.357831e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.HZLimRecVenus": {"value": 1.118929e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.HZLimRunaway": {"value": 1.461108e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.HZLimMoistGreenhouse": {"value": 1.480517e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.HZLimMaxGreenhouse": {"value": 2.509538e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.HZLimEarlyMars": {"value": 2.738109e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.Instellation": {"value": -1.000000, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.sun.CriticalSemiMajorAxis": {"value": -1.000000, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.LXUVTot": {"value": 3.846000e+23, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.sun.LostEnergy": {"value": 2.568599e+28, "unit": u.Joule, "rtol": 1e-4}, + "log.final.sun.LostAngMom": {"value": 3.532078e+32, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, + "log.final.sun.EscapeVelocity": {"value": 1.146413e+06, "unit": u.m / u.sec, "rtol": 1e-4}, + "log.final.sun.Luminosity": {"value": 3.846000e+26, "unit": u.W, "rtol": 1e-4}, + "log.final.sun.LXUVStellar": {"value": 3.846000e+23, "unit": u.W, "rtol": 1e-4}, + "log.final.sun.Temperature": {"value": 5778.000000, "unit": u.K, "rtol": 1e-4}, + "log.final.sun.LXUVFrac": {"value": 0.001000, "rtol": 1e-4}, + "log.final.sun.RossbyNumber": {"value": 0.078260, "rtol": 1e-4}, + "log.final.sun.DRotPerDtStellar": {"value": 6.558557e-13, "rtol": 1e-4}, + "log.final.sun.WindTorque": {"value": 1.119248e+25, "rtol": 1e-4}, + "log.final.random.Mass": {"value": 5.971546e+24, "unit": u.kg, "rtol": 1e-4}, + "log.final.random.Obliquity": {"value": 0.410152, "unit": u.rad, "rtol": 1e-4}, + "log.final.random.PrecA": {"value": 0.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.random.Radius": {"value": 6.378100e+06, "unit": u.m, "rtol": 1e-4}, + "log.final.random.RadGyra": {"value": 0.500000, "rtol": 1e-4}, + "log.final.random.BodyType": {"value": 0.000000, "rtol": 1e-4}, + "log.final.random.Density": {"value": 5494.449526, "unit": u.kg / u.m ** 3, "rtol": 1e-4}, + "log.final.random.HZLimitDryRunaway": {"value": 1.112478e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.random.HZLimRecVenus": {"value": 1.118929e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.random.HZLimRunaway": {"value": 1.461108e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.random.HZLimMoistGreenhouse": {"value": 1.480517e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.random.HZLimMaxGreenhouse": {"value": 2.509538e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.random.HZLimEarlyMars": {"value": 2.738109e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.random.Instellation": {"value": 1367.566935, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.random.Eccentricity": {"value": 0.000000, "rtol": 1e-4}, + "log.final.random.MeanMotion": {"value": 1.990987e-07, "unit": 1 / u.sec, "rtol": 1e-4}, + "log.final.random.OrbPeriod": {"value": 3.155815e+07, "unit": u.sec, "rtol": 1e-4}, + "log.final.random.SemiMajorAxis": {"value": 1.495979e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.random.COPP": {"value": 0.000000, "rtol": 1e-4}, + "log.final.random.EscapeVelocity": {"value": 1.117931e+04, "unit": u.m / u.sec, "rtol": 1e-4}, + "log.final.random.TGlobal": {"value": 17.294952, "unit": u.deg_C, "rtol": 1e-4}, + "log.final.random.AlbedoGlobal": {"value": 0.328739, "rtol": 1e-4}, + "log.final.random.FluxInGlobal": {"value": 238.245969, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.random.FluxOutGlobal": {"value": 239.446450, "unit": u.W / u.m ** 2, "rtol": 1e-4}, + "log.final.random.TotIceMass": {"value": 2.595204e+13, "unit": u.kg, "rtol": 1e-4}, + "log.final.random.TotIceFlow": {"value": 0.000000, "unit": u.kg, "rtol": 1e-4}, + "log.final.random.TotIceBalance": {"value": 1.798315e-09, "unit": u.kg, "rtol": 1e-4}, + "log.final.random.SkipSeas": {"value": 0.000000, "rtol": 1e-4}, + "log.final.random.AreaIceCov": {"value": 0.056722, "rtol": 1e-4}, + "log.final.random.Latitude": {"value": 83.402352, "unit": u.deg, "rtol": 1e-4}, + "log.final.random.TempLat": {"value": -11.330438, "unit": u.deg_C, "rtol": 1e-4}, + "log.final.random.AlbedoLat": {"value": 0.600000, "rtol": 1e-4}, + "log.final.random.AnnInsol": {"value": 176.064060, "unit": u.W / u.m ** 2, "rtol": 1e-4}, + "log.final.random.FluxMerid": {"value": 1.717496, "unit": u.PW, "rtol": 1e-4}, + "log.final.random.FluxIn": {"value": 70.432311, "unit": u.W / u.m ** 2, "rtol": 1e-4}, + "log.final.random.FluxOut": {"value": 179.619385, "unit": u.W / u.m ** 2, "rtol": 1e-4}, + "log.final.random.DivFlux": {"value": -116.190681, "unit": u.W / u.m ** 2, "rtol": 1e-4}, + "log.final.random.IceMass": {"value": 52.023175, "rtol": 1e-4}, + "log.final.random.IceHeight": {"value": 0.056750, "unit": u.m, "rtol": 1e-4}, + "log.final.random.DIceMassDt": {"value": 1.648515e-06, "unit": u.m, "rtol": 1e-4}, + "log.final.random.IceFlow": {"value": 0.000000, "unit": u.m / u.sec, "rtol": 1e-4}, + "log.final.random.EnergyResL": {"value": -1.438093, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.random.EnergyResW": {"value": 0.179448, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.random.BedrockH": {"value": -3.087429e-06, "unit": u.m, "rtol": 1e-4}, + "log.final.random.TempLandLat": {"value": 259.880387, "unit": u.sec, "rtol": 1e-4}, + "log.final.random.TempWaterLat": {"value": 261.978765, "unit": u.sec, "rtol": 1e-4}, + "log.final.random.AlbedoLandLat": {"value": 0.600000, "rtol": 1e-4}, + "log.final.random.AlbedoWaterLat": {"value": 0.600000, "rtol": 1e-4}, + "log.final.random.TempMinLat": {"value": -14.519534, "unit": u.deg_C, "rtol": 1e-4}, + "log.final.random.TempMaxLat": {"value": -7.644358, "unit": u.deg_C, "rtol": 1e-4}, + "log.final.random.Snowball": {"value": 0.000000, "rtol": 1e-4}, + "log.final.random.PlanckBAvg": {"value": 2.090000, "rtol": 1e-4}, + "log.final.random.IceAccum": {"value": 0.614824, "rtol": 1e-4}, + "log.final.random.IceAblate": {"value": -0.398709, "rtol": 1e-4}, + "log.final.random.TempMaxLand": {"value": 275.152752, "unit": u.sec, "rtol": 1e-4}, + "log.final.random.TempMaxWater": {"value": 263.843906, "unit": u.sec, "rtol": 1e-4}, + "log.final.random.PeakInsol": {"value": 541.684612, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.random.IceCapNorthLand": {"value": 1.000000, "rtol": 1e-4}, + "log.final.random.IceCapNorthSea": {"value": 1.000000, "rtol": 1e-4}, + "log.final.random.IceCapSouthLand": {"value": 0.000000, "rtol": 1e-4}, + "log.final.random.IceCapSouthSea": {"value": 1.000000, "rtol": 1e-4}, + "log.final.random.IceBeltLand": {"value": 0.000000, "rtol": 1e-4}, + "log.final.random.IceBeltSea": {"value": 0.000000, "rtol": 1e-4}, + "log.final.random.SnowballLand": {"value": 0.000000, "rtol": 1e-4}, + "log.final.random.SnowballSea": {"value": 0.000000, "rtol": 1e-4}, + "log.final.random.IceFree": {"value": 0.000000, "rtol": 1e-4}, + "log.final.random.IceCapNorthLatLand": {"value": 1.371128, "unit": u.rad, "rtol": 1e-4}, + "log.final.random.IceCapNorthLatSea": {"value": 1.152808, "unit": u.rad, "rtol": 1e-4}, + "log.final.random.IceCapSouthLatLand": {"value": 100.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.random.IceCapSouthLatSea": {"value": 100.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.random.IceBeltNorthLatLand": {"value": 100.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.random.IceBeltNorthLatSea": {"value": 100.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.random.IceBeltSouthLatLand": {"value": 100.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.random.IceBeltSouthLatSea": {"value": 100.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.random.LandFracGlobal": {"value": 0.292000, "rtol": 1e-4}, + } +) +class Test_SeasonalNC79Random(Benchmark): + pass diff --git a/tests/Poise/SeasonalNC79Random/vpl.in b/tests/Poise/SeasonalNC79Random/vpl.in new file mode 100755 index 000000000..e5a7ff93d --- /dev/null +++ b/tests/Poise/SeasonalNC79Random/vpl.in @@ -0,0 +1,15 @@ +sSystemName climateland +iVerbose 0 #how much do you want vplanet to yell at you? +iDigits 6 #how many digits do you want in your numbers? +bOverwrite 1 #overwrite old files +sUnitMass solar #mass unit used for input +sUnitLength au #length unit used for input +sUnitTime y #time unit +sUnitAngle d #angle unit +bDoLog 1 #create log file +saBodyFiles sun.in random.in +bDoForward 1 #integrate forward in time +bVarDt 1 #use variable time stepping (not relevant to poise) +dEta 0.01 #how much to scale variable time step +dStopTime 1 #how long should the integration be +dOutputTime 1 #how much output you want diff --git a/tests/Poise/SeasonalNC79Uniform/sun.in b/tests/Poise/SeasonalNC79Uniform/sun.in new file mode 100755 index 000000000..a1f97c666 --- /dev/null +++ b/tests/Poise/SeasonalNC79Uniform/sun.in @@ -0,0 +1,9 @@ +# sun parameters +sName sun +dMass 1 +dSemi 0 +dEcc 0 +dRadius 0.00135 +dLuminosity -1 +sStellarModel none #sun does not change over time +saModules stellar #use stellar module (needed for luminosity) diff --git a/tests/Poise/SeasonalNC79Uniform/test_SeasonalNC79Uniform.py b/tests/Poise/SeasonalNC79Uniform/test_SeasonalNC79Uniform.py new file mode 100644 index 000000000..bd7242a76 --- /dev/null +++ b/tests/Poise/SeasonalNC79Uniform/test_SeasonalNC79Uniform.py @@ -0,0 +1,228 @@ +from benchmark import Benchmark, benchmark +import astropy.units as u + +@benchmark( + { + "log.initial.system.Age": {"value": 0.000000, "unit": u.sec}, + "log.initial.system.Time": {"value": 0.000000, "unit": u.sec}, + "log.initial.system.TotAngMom": {"value": 1.501063e+42, "unit": (u.kg * u.m ** 2) / u.sec}, + "log.initial.system.TotEnergy": {"value": -7.839372e+41, "unit": u.Joule}, + "log.initial.system.PotEnergy": {"value": -7.839908e+41, "unit": u.Joule}, + "log.initial.system.KinEnergy": {"value": 5.361272e+37, "unit": u.Joule}, + "log.initial.system.DeltaTime": {"value": 0.000000, "unit": u.sec}, + "log.initial.sun.Mass": {"value": 1.988416e+30, "unit": u.kg}, + "log.initial.sun.Radius": {"value": 2.019571e+08, "unit": u.m}, + "log.initial.sun.RadGyra": {"value": 0.500000}, + "log.initial.sun.RotAngMom": {"value": 1.474456e+42, "unit": (u.kg * u.m ** 2) / u.sec}, + "log.initial.sun.RotVel": {"value": 1.468674e+04, "unit": u.m / u.sec}, + "log.initial.sun.BodyType": {"value": 0.000000}, + "log.initial.sun.RotRate": {"value": 7.272205e-05, "unit": 1 / u.sec}, + "log.initial.sun.RotPer": {"value": 8.640000e+04, "unit": u.sec}, + "log.initial.sun.Density": {"value": 5.762900e+04, "unit": u.kg / u.m ** 3}, + "log.initial.sun.HZLimitDryRunaway": {"value": 1.357831e+11, "unit": u.m}, + "log.initial.sun.HZLimRecVenus": {"value": 1.118929e+11, "unit": u.m}, + "log.initial.sun.HZLimRunaway": {"value": 1.461108e+11, "unit": u.m}, + "log.initial.sun.HZLimMoistGreenhouse": {"value": 1.480517e+11, "unit": u.m}, + "log.initial.sun.HZLimMaxGreenhouse": {"value": 2.509538e+11, "unit": u.m}, + "log.initial.sun.HZLimEarlyMars": {"value": 2.738109e+11, "unit": u.m}, + "log.initial.sun.Instellation": {"value": -1.000000, "unit": u.kg / u.sec ** 3}, + "log.initial.sun.CriticalSemiMajorAxis": {"value": -1.000000, "unit": u.m}, + "log.initial.sun.LXUVTot": {"value": 3.846000e+23, "unit": u.kg / u.sec ** 3}, + "log.initial.sun.LostEnergy": {"value": 5.562685e-309, "unit": u.Joule}, + "log.initial.sun.LostAngMom": {"value": 5.562685e-309, "unit": (u.kg * u.m ** 2) / u.sec}, + "log.initial.sun.EscapeVelocity": {"value": 1.146413e+06, "unit": u.m / u.sec}, + "log.initial.sun.Luminosity": {"value": 3.846000e+26, "unit": u.W}, + "log.initial.sun.LXUVStellar": {"value": 3.846000e+23, "unit": u.W}, + "log.initial.sun.Temperature": {"value": 5778.000000, "unit": u.K}, + "log.initial.sun.LXUVFrac": {"value": 0.001000}, + "log.initial.sun.RossbyNumber": {"value": 0.078260}, + "log.initial.sun.DRotPerDtStellar": {"value": 6.558557e-13}, + "log.initial.sun.WindTorque": {"value": 1.119248e+25}, + "log.initial.uniform.Mass": {"value": 5.971546e+24, "unit": u.kg}, + "log.initial.uniform.Obliquity": {"value": 0.410152, "unit": u.rad}, + "log.initial.uniform.PrecA": {"value": 0.000000, "unit": u.rad}, + "log.initial.uniform.Radius": {"value": 6.378100e+06, "unit": u.m}, + "log.initial.uniform.RadGyra": {"value": 0.500000}, + "log.initial.uniform.BodyType": {"value": 0.000000}, + "log.initial.uniform.Density": {"value": 5494.449526, "unit": u.kg / u.m ** 3}, + "log.initial.uniform.HZLimitDryRunaway": {"value": 1.108439e+11, "unit": u.m}, + "log.initial.uniform.HZLimRecVenus": {"value": 1.118929e+11, "unit": u.m}, + "log.initial.uniform.HZLimRunaway": {"value": 1.461108e+11, "unit": u.m}, + "log.initial.uniform.HZLimMoistGreenhouse": {"value": 1.480517e+11, "unit": u.m}, + "log.initial.uniform.HZLimMaxGreenhouse": {"value": 2.509538e+11, "unit": u.m}, + "log.initial.uniform.HZLimEarlyMars": {"value": 2.738109e+11, "unit": u.m}, + "log.initial.uniform.Instellation": {"value": 1367.566935, "unit": u.kg / u.sec ** 3}, + "log.initial.uniform.Eccentricity": {"value": 0.000000}, + "log.initial.uniform.MeanMotion": {"value": 1.990987e-07, "unit": 1 / u.sec}, + "log.initial.uniform.OrbPeriod": {"value": 3.155815e+07, "unit": u.sec}, + "log.initial.uniform.SemiMajorAxis": {"value": 1.495979e+11, "unit": u.m}, + "log.initial.uniform.COPP": {"value": 0.000000}, + "log.initial.uniform.EscapeVelocity": {"value": 1.117931e+04, "unit": u.m / u.sec}, + "log.initial.uniform.TGlobal": {"value": 16.177005, "unit": u.deg_C}, + "log.initial.uniform.AlbedoGlobal": {"value": 0.333604}, + "log.initial.uniform.FluxInGlobal": {"value": 236.924441, "unit": u.kg / u.sec ** 3}, + "log.initial.uniform.FluxOutGlobal": {"value": 237.109940, "unit": u.W / u.m ** 2}, + "log.initial.uniform.TotIceMass": {"value": 0.000000, "unit": u.kg}, + "log.initial.uniform.TotIceFlow": {"value": 0.000000, "unit": u.kg}, + "log.initial.uniform.TotIceBalance": {"value": 0.000000, "unit": u.kg}, + "log.initial.uniform.SkipSeas": {"value": 0.000000}, + "log.initial.uniform.AreaIceCov": {"value": 0.065642}, + "log.initial.uniform.Latitude": {"value": -83.402352, "unit": u.deg}, + "log.initial.uniform.TempLat": {"value": -13.449779, "unit": u.deg_C}, + "log.initial.uniform.AlbedoLat": {"value": 0.600000}, + "log.initial.uniform.AnnInsol": {"value": 176.064060, "unit": u.W / u.m ** 2}, + "log.initial.uniform.FluxMerid": {"value": -1.510862, "unit": u.PW}, + "log.initial.uniform.FluxIn": {"value": 70.432311, "unit": u.W / u.m ** 2}, + "log.initial.uniform.FluxOut": {"value": 175.189962, "unit": u.W / u.m ** 2}, + "log.initial.uniform.DivFlux": {"value": -102.211645, "unit": u.W / u.m ** 2}, + "log.initial.uniform.IceMass": {"value": 0.000000}, + "log.initial.uniform.IceHeight": {"value": 0.000000, "unit": u.m}, + "log.initial.uniform.DIceMassDt": {"value": 0.000000, "unit": u.m}, + "log.initial.uniform.IceFlow": {"value": 0.000000, "unit": u.m / u.sec}, + "log.initial.uniform.EnergyResL": {"value": -0.923722, "unit": u.kg / u.sec ** 3}, + "log.initial.uniform.EnergyResW": {"value": 0.216110, "unit": u.kg / u.sec ** 3}, + "log.initial.uniform.BedrockH": {"value": 0.000000, "unit": u.m}, + "log.initial.uniform.TempLandLat": {"value": 257.607045, "unit": u.sec}, + "log.initial.uniform.TempWaterLat": {"value": 260.351644, "unit": u.sec}, + "log.initial.uniform.AlbedoLandLat": {"value": 0.600000}, + "log.initial.uniform.AlbedoWaterLat": {"value": 0.600000}, + "log.initial.uniform.TempMinLat": {"value": -18.988426, "unit": u.deg_C}, + "log.initial.uniform.TempMaxLat": {"value": -7.191886, "unit": u.deg_C}, + "log.initial.uniform.Snowball": {"value": 0.000000}, + "log.initial.uniform.PlanckBAvg": {"value": 2.090000}, + "log.initial.uniform.IceAccum": {"value": 0.609982}, + "log.initial.uniform.IceAblate": {"value": -0.263629}, + "log.initial.uniform.TempMaxLand": {"value": 274.781733, "unit": u.sec}, + "log.initial.uniform.TempMaxWater": {"value": 262.300735, "unit": u.sec}, + "log.initial.uniform.PeakInsol": {"value": 541.704677, "unit": u.kg / u.sec ** 3}, + "log.initial.uniform.IceCapNorthLand": {"value": 0.000000}, + "log.initial.uniform.IceCapNorthSea": {"value": 1.000000}, + "log.initial.uniform.IceCapSouthLand": {"value": 0.000000}, + "log.initial.uniform.IceCapSouthSea": {"value": 1.000000}, + "log.initial.uniform.IceBeltLand": {"value": 0.000000}, + "log.initial.uniform.IceBeltSea": {"value": 0.000000}, + "log.initial.uniform.SnowballLand": {"value": 0.000000}, + "log.initial.uniform.SnowballSea": {"value": 0.000000}, + "log.initial.uniform.IceFree": {"value": 0.000000}, + "log.initial.uniform.IceCapNorthLatLand": {"value": 100.000000, "unit": u.rad}, + "log.initial.uniform.IceCapNorthLatSea": {"value": 1.121291, "unit": u.rad}, + "log.initial.uniform.IceCapSouthLatLand": {"value": 100.000000, "unit": u.rad}, + "log.initial.uniform.IceCapSouthLatSea": {"value": 100.000000, "unit": u.rad}, + "log.initial.uniform.IceBeltNorthLatLand": {"value": 100.000000, "unit": u.rad}, + "log.initial.uniform.IceBeltNorthLatSea": {"value": 100.000000, "unit": u.rad}, + "log.initial.uniform.IceBeltSouthLatLand": {"value": 100.000000, "unit": u.rad}, + "log.initial.uniform.IceBeltSouthLatSea": {"value": 100.000000, "unit": u.rad}, + "log.initial.uniform.LandFracGlobal": {"value": 0.292000}, + "log.final.system.Age": {"value": 3.155760e+07, "unit": u.sec, "rtol": 1e-4}, + "log.final.system.Time": {"value": 3.155760e+07, "unit": u.sec, "rtol": 1e-4}, + "log.final.system.TotAngMom": {"value": 1.501063e+42, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, + "log.final.system.TotEnergy": {"value": -7.839372e+41, "unit": u.Joule, "rtol": 1e-4}, + "log.final.system.PotEnergy": {"value": -7.839908e+41, "unit": u.Joule, "rtol": 1e-4}, + "log.final.system.KinEnergy": {"value": 5.361272e+37, "unit": u.Joule, "rtol": 1e-4}, + "log.final.sun.Mass": {"value": 1.988416e+30, "unit": u.kg, "rtol": 1e-4}, + "log.final.sun.Radius": {"value": 2.019571e+08, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.RadGyra": {"value": 0.500000, "rtol": 1e-4}, + "log.final.sun.RotAngMom": {"value": 1.474456e+42, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, + "log.final.sun.RotVel": {"value": 1.468674e+04, "unit": u.m / u.sec, "rtol": 1e-4}, + "log.final.sun.BodyType": {"value": 0.000000, "rtol": 1e-4}, + "log.final.sun.RotRate": {"value": 7.272205e-05, "unit": 1 / u.sec, "rtol": 1e-4}, + "log.final.sun.RotPer": {"value": 8.640000e+04, "unit": u.sec, "rtol": 1e-4}, + "log.final.sun.Density": {"value": 5.762900e+04, "unit": u.kg / u.m ** 3, "rtol": 1e-4}, + "log.final.sun.HZLimitDryRunaway": {"value": 1.357831e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.HZLimRecVenus": {"value": 1.118929e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.HZLimRunaway": {"value": 1.461108e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.HZLimMoistGreenhouse": {"value": 1.480517e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.HZLimMaxGreenhouse": {"value": 2.509538e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.HZLimEarlyMars": {"value": 2.738109e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.Instellation": {"value": -1.000000, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.sun.CriticalSemiMajorAxis": {"value": -1.000000, "unit": u.m, "rtol": 1e-4}, + "log.final.sun.LXUVTot": {"value": 3.846000e+23, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.sun.LostEnergy": {"value": 2.568599e+28, "unit": u.Joule, "rtol": 1e-4}, + "log.final.sun.LostAngMom": {"value": 3.532078e+32, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, + "log.final.sun.EscapeVelocity": {"value": 1.146413e+06, "unit": u.m / u.sec, "rtol": 1e-4}, + "log.final.sun.Luminosity": {"value": 3.846000e+26, "unit": u.W, "rtol": 1e-4}, + "log.final.sun.LXUVStellar": {"value": 3.846000e+23, "unit": u.W, "rtol": 1e-4}, + "log.final.sun.Temperature": {"value": 5778.000000, "unit": u.K, "rtol": 1e-4}, + "log.final.sun.LXUVFrac": {"value": 0.001000, "rtol": 1e-4}, + "log.final.sun.RossbyNumber": {"value": 0.078260, "rtol": 1e-4}, + "log.final.sun.DRotPerDtStellar": {"value": 6.558557e-13, "rtol": 1e-4}, + "log.final.sun.WindTorque": {"value": 1.119248e+25, "rtol": 1e-4}, + "log.final.uniform.Mass": {"value": 5.971546e+24, "unit": u.kg, "rtol": 1e-4}, + "log.final.uniform.Obliquity": {"value": 0.410152, "unit": u.rad, "rtol": 1e-4}, + "log.final.uniform.PrecA": {"value": 0.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.uniform.Radius": {"value": 6.378100e+06, "unit": u.m, "rtol": 1e-4}, + "log.final.uniform.RadGyra": {"value": 0.500000, "rtol": 1e-4}, + "log.final.uniform.BodyType": {"value": 0.000000, "rtol": 1e-4}, + "log.final.uniform.Density": {"value": 5494.449526, "unit": u.kg / u.m ** 3, "rtol": 1e-4}, + "log.final.uniform.HZLimitDryRunaway": {"value": 1.108439e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.uniform.HZLimRecVenus": {"value": 1.118929e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.uniform.HZLimRunaway": {"value": 1.461108e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.uniform.HZLimMoistGreenhouse": {"value": 1.480517e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.uniform.HZLimMaxGreenhouse": {"value": 2.509538e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.uniform.HZLimEarlyMars": {"value": 2.738109e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.uniform.Instellation": {"value": 1367.566935, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.uniform.Eccentricity": {"value": 0.000000, "rtol": 1e-4}, + "log.final.uniform.MeanMotion": {"value": 1.990987e-07, "unit": 1 / u.sec, "rtol": 1e-4}, + "log.final.uniform.OrbPeriod": {"value": 3.155815e+07, "unit": u.sec, "rtol": 1e-4}, + "log.final.uniform.SemiMajorAxis": {"value": 1.495979e+11, "unit": u.m, "rtol": 1e-4}, + "log.final.uniform.COPP": {"value": 0.000000, "rtol": 1e-4}, + "log.final.uniform.EscapeVelocity": {"value": 1.117931e+04, "unit": u.m / u.sec, "rtol": 1e-4}, + "log.final.uniform.TGlobal": {"value": 16.181001, "unit": u.deg_C, "rtol": 1e-4}, + "log.final.uniform.AlbedoGlobal": {"value": 0.333604, "rtol": 1e-4}, + "log.final.uniform.FluxInGlobal": {"value": 236.924441, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.uniform.FluxOutGlobal": {"value": 237.118291, "unit": u.W / u.m ** 2, "rtol": 1e-4}, + "log.final.uniform.TotIceMass": {"value": 1.293774e+14, "unit": u.kg, "rtol": 1e-4}, + "log.final.uniform.TotIceFlow": {"value": 0.000000, "unit": u.kg, "rtol": 1e-4}, + "log.final.uniform.TotIceBalance": {"value": 4.524064e-09, "unit": u.kg, "rtol": 1e-4}, + "log.final.uniform.SkipSeas": {"value": 0.000000, "rtol": 1e-4}, + "log.final.uniform.AreaIceCov": {"value": 0.071444, "rtol": 1e-4}, + "log.final.uniform.Latitude": {"value": 83.402352, "unit": u.deg, "rtol": 1e-4}, + "log.final.uniform.TempLat": {"value": -13.511537, "unit": u.deg_C, "rtol": 1e-4}, + "log.final.uniform.AlbedoLat": {"value": 0.600000, "rtol": 1e-4}, + "log.final.uniform.AnnInsol": {"value": 176.064060, "unit": u.W / u.m ** 2, "rtol": 1e-4}, + "log.final.uniform.FluxMerid": {"value": 1.506570, "unit": u.PW, "rtol": 1e-4}, + "log.final.uniform.FluxIn": {"value": 70.432311, "unit": u.W / u.m ** 2, "rtol": 1e-4}, + "log.final.uniform.FluxOut": {"value": 175.060887, "unit": u.W / u.m ** 2, "rtol": 1e-4}, + "log.final.uniform.DivFlux": {"value": -101.921232, "unit": u.W / u.m ** 2, "rtol": 1e-4}, + "log.final.uniform.IceMass": {"value": 65.838491, "rtol": 1e-4}, + "log.final.uniform.IceHeight": {"value": 0.071821, "unit": u.m, "rtol": 1e-4}, + "log.final.uniform.DIceMassDt": {"value": 2.086296e-06, "unit": u.m, "rtol": 1e-4}, + "log.final.uniform.IceFlow": {"value": 0.000000, "unit": u.m / u.sec, "rtol": 1e-4}, + "log.final.uniform.EnergyResL": {"value": -0.923722, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.uniform.EnergyResW": {"value": 0.216110, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.uniform.BedrockH": {"value": -3.907329e-06, "unit": u.m, "rtol": 1e-4}, + "log.final.uniform.TempLandLat": {"value": 257.571610, "unit": u.sec, "rtol": 1e-4}, + "log.final.uniform.TempWaterLat": {"value": 260.279029, "unit": u.sec, "rtol": 1e-4}, + "log.final.uniform.AlbedoLandLat": {"value": 0.600000, "rtol": 1e-4}, + "log.final.uniform.AlbedoWaterLat": {"value": 0.600000, "rtol": 1e-4}, + "log.final.uniform.TempMinLat": {"value": -19.044471, "unit": u.deg_C, "rtol": 1e-4}, + "log.final.uniform.TempMaxLat": {"value": -7.261100, "unit": u.deg_C, "rtol": 1e-4}, + "log.final.uniform.Snowball": {"value": 0.000000, "rtol": 1e-4}, + "log.final.uniform.PlanckBAvg": {"value": 2.090000, "rtol": 1e-4}, + "log.final.uniform.IceAccum": {"value": 0.622085, "rtol": 1e-4}, + "log.final.uniform.IceAblate": {"value": -0.369692, "rtol": 1e-4}, + "log.final.uniform.TempMaxLand": {"value": 274.740432, "unit": u.sec, "rtol": 1e-4}, + "log.final.uniform.TempMaxWater": {"value": 262.225513, "unit": u.sec, "rtol": 1e-4}, + "log.final.uniform.PeakInsol": {"value": 541.684612, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, + "log.final.uniform.IceCapNorthLand": {"value": 1.000000, "rtol": 1e-4}, + "log.final.uniform.IceCapNorthSea": {"value": 1.000000, "rtol": 1e-4}, + "log.final.uniform.IceCapSouthLand": {"value": 1.000000, "rtol": 1e-4}, + "log.final.uniform.IceCapSouthSea": {"value": 1.000000, "rtol": 1e-4}, + "log.final.uniform.IceBeltLand": {"value": 0.000000, "rtol": 1e-4}, + "log.final.uniform.IceBeltSea": {"value": 0.000000, "rtol": 1e-4}, + "log.final.uniform.SnowballLand": {"value": 0.000000, "rtol": 1e-4}, + "log.final.uniform.SnowballSea": {"value": 0.000000, "rtol": 1e-4}, + "log.final.uniform.IceFree": {"value": 0.000000, "rtol": 1e-4}, + "log.final.uniform.IceCapNorthLatLand": {"value": 1.312738, "unit": u.rad, "rtol": 1e-4}, + "log.final.uniform.IceCapNorthLatSea": {"value": 1.121291, "unit": u.rad, "rtol": 1e-4}, + "log.final.uniform.IceCapSouthLatLand": {"value": -1.371128, "unit": u.rad, "rtol": 1e-4}, + "log.final.uniform.IceCapSouthLatSea": {"value": -1.371128, "unit": u.rad, "rtol": 1e-4}, + "log.final.uniform.IceBeltNorthLatLand": {"value": 100.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.uniform.IceBeltNorthLatSea": {"value": 100.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.uniform.IceBeltSouthLatLand": {"value": 100.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.uniform.IceBeltSouthLatSea": {"value": 100.000000, "unit": u.rad, "rtol": 1e-4}, + "log.final.uniform.LandFracGlobal": {"value": 0.292000, "rtol": 1e-4}, + } +) +class Test_SeasonalNC79Uniform(Benchmark): + pass diff --git a/tests/Poise/SeasonalNC79Uniform/uniform.in b/tests/Poise/SeasonalNC79Uniform/uniform.in new file mode 100755 index 000000000..c7423f090 --- /dev/null +++ b/tests/Poise/SeasonalNC79Uniform/uniform.in @@ -0,0 +1,71 @@ +sName uniform #name of planet +saModules poise #what vplanet modules you want to use +sGeography uniform +dLandFrac 0.292 + +dMass 3.00316726e-06 #mass of planet +dRadius -1.00 #radius (not important right now) +dRotPeriod -1.00000 #rotation period (minus = days) +dObliquity 23.5 +dSemi 1 +dEcc 0.0 #eccentricity of orbit +dLongP 0 #pericenter, wrt Earth's position at spring equinox +# note that this is the typical value +180, +# since that one is solar position +dDynEllip 0.0 #shape of planet (0 = a sphere) +dPrecA 0.0 #orientation of spin axis + +#_______addition disorb/distrot parameters (leave these alone for now)__________________ +#dInc 5e-5 #inclination of orbit +#dLongA 348.73936 #orientation of orbital plane +#bGRCorr 0 #use GR correction (not important) +#bInvPlane 1 #convert to invariable plane coords +#bOverrideMaxEcc 1 #override max ecc halt (not recommended) +#dHaltMaxEcc 0.4 #eccentricity at which to halt simulation + +#_______poise parameters (have fun with these!)_________________________________________ +iLatCellNum 151 #number of latitude cells +sClimateModel sea #use seasonal or annual model +dTGlobalInit 14.85 #initial guess at average surface temp +iNumYears 4 #number of years (orbits) to run clim model +iNStepInYear 80 #number of steps to take in a "year" +#dSurfAlbedo 0.35 #average surface albedo (annual model only) + +#__ice params_________ +bIceSheets 1 #enable ice sheets +dInitIceLat 90. #how low do initial ice sheet extend? +dInitIceHeight 0. #height of initial ice sheets +dIceDepRate 2.25e-5 #rate of snow build up (when T < 0) +dIceAlbedo 0.6 #albedo of ice +iIceDt 1 #time step of ice-sheet model (orbits) +iReRunSeas 500 #how often to re-run seasonal model +bSeaIceModel 0 #use sea ice model (slow!) +bSkipSeasEnabled 0 #can skip seasonal if snowball state present + +#__heat diffusion______ +#bMEPDiff 1 #calculate diffusion using max entropy production +#bHadley 1 #mimic hadley heat diffusion +dDiffusion 0.58 #diffusion coefficient (fixed) +dNuLandWater 0.8 #Heat diffusion coefficient between Land and Water + +#__outgoing flux_______ +dPlanckA 203.3 #offset for OLR calculation (greenhouse) +dPlanckB 2.09 #slope of OLR calc (water vapor feedback) +bCalcAB 0 #calculate A & B from Kasting model fits +#dpCO2 0.00028 #partial pressure of co2 + +#__surface properties__ +dAlbedoLand 0.363 #albedo of land +dAlbedoWater 0.263 #albedo of water +dHeatCapLand 1.55e7 #land heat capacity +dHeatCapWater 4.428e6 #water heat capacity +dMixingDepth 70 #mixing depth of ocean + + +#________output options!_____________________________________________ +saOutputOrder Time PrecA -TGlobal AlbedoGlobal -FluxOutGlobal $ + -TotIceMass -TotIceFlow -TotIceBalance DeltaTime AreaIceCov Snowball Obliq $ + IceBeltLand IceBeltSea Ecce + +saGridOutput Time -Latitude -TempLat AlbedoLat -AnnInsol -FluxIn -FluxOut IceMass -IceHeight DIceMassDt $ + -IceFlow -BedrockH -TempMaxLat -TempMinLat -FluxMerid -DivFlux diff --git a/tests/Poise/SeasonalNC79Uniform/vpl.in b/tests/Poise/SeasonalNC79Uniform/vpl.in new file mode 100755 index 000000000..aee523207 --- /dev/null +++ b/tests/Poise/SeasonalNC79Uniform/vpl.in @@ -0,0 +1,15 @@ +sSystemName climateland +iVerbose 0 #how much do you want vplanet to yell at you? +iDigits 6 #how many digits do you want in your numbers? +bOverwrite 1 #overwrite old files +sUnitMass solar #mass unit used for input +sUnitLength au #length unit used for input +sUnitTime y #time unit +sUnitAngle d #angle unit +bDoLog 1 #create log file +saBodyFiles sun.in uniform.in +bDoForward 1 #integrate forward in time +bVarDt 1 #use variable time stepping (not relevant to poise) +dEta 0.01 #how much to scale variable time step +dStopTime 1 #how long should the integration be +dOutputTime 1 #how much output you want