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/Makefile b/Makefile index a0fea5c7c..995a8da3a 100644 --- a/Makefile +++ b/Makefile @@ -4,10 +4,12 @@ GITVERSION := $(shell git describe --tags --abbrev=40 --always) UNAME_S := $(shell uname -s) ifeq ($(UNAME_S),Linux) +CC = gcc GCC_FLAGS1 = -fPIC -Wl,-Bsymbolic-functions -c GCC_FLAGS2 = -shared -Wl,-Bsymbolic-functions,-soname,vplanetlib.so endif ifeq ($(UNAME_S),Darwin) +CC = clang GCC_FLAGS1 = -fPIC -c GCC_FLAGS2 = -shared -Wl,-install_name,vplanetlib.so endif @@ -17,7 +19,7 @@ default: -python setup.py develop legacy: - -gcc -o bin/vplanet src/*.c -lm -DGITVERSION=\"$(GITVERSION)\" + -$(CC) -o bin/vplanet src/*.c -lm -DGITVERSION=\"$(GITVERSION)\" @echo "" @echo "==========================================================================================================" @echo 'To add vplanet to your $$PATH, please run the appropriate command for your shell type:' @@ -30,13 +32,13 @@ legacy: @echo "==========================================================================================================" debug: - -gcc -g -D DEBUG -o bin/vplanet src/*.c -lm -DGITVERSION=\"$(GITVERSION)\" + -$(CC) -g -D DEBUG -o bin/vplanet src/*.c -lm -DGITVERSION=\"$(GITVERSION)\" debug_no_AE: - -gcc -g -o bin/vplanet src/*.c -lm -DGITVERSION=\"$(GITVERSION)\" + -$(CC) -g -o bin/vplanet src/*.c -lm -DGITVERSION=\"$(GITVERSION)\" opt: - -gcc -o bin/vplanet src/*.c -lm -O3 -DGITVERSION=\"$(GITVERSION)\" + -$(CC) -o bin/vplanet src/*.c -lm -O3 -DGITVERSION=\"$(GITVERSION)\" @echo "" @echo "==========================================================================================================" @echo 'To add vplanet to your $$PATH, please run the appropriate command for your shell type:' @@ -52,28 +54,28 @@ cpp: g++ -o bin/vplanet src/*.c -lm -O3 -fopenmp -fpermissive -w -DGITVERSION=\"$(GITVERSION)\" warnings: - -gcc -g -D DEBUG -Wunused-but-set-variable -Wunused-variable -Wfloat-equal -o bin/vplanet src/*.c -lm -DGITVERSION=\"$(GITVERSION)\" + -$(CC)-g -D DEBUG -Wunused-but-set-variable -Wunused-variable -Wfloat-equal -o bin/vplanet src/*.c -lm -DGITVERSION=\"$(GITVERSION)\" parallel: - gcc -o bin/vplanet src/*.c -lm -O3 -fopenmp -DGITVERSION=\"$(GITVERSION)\" + $(CC) -o bin/vplanet src/*.c -lm -O3 -fopenmp -DGITVERSION=\"$(GITVERSION)\" profile: - -gcc -pg -o bin/vplanet src/*.c -lm -DGITVERSION=\"$(GITVERSION)\" + -$(CC) -pg -o bin/vplanet src/*.c -lm -DGITVERSION=\"$(GITVERSION)\" optprof: - -gcc -pg -o bin/vplanet src/*.c -lm -O3 -DGITVERSION=\"$(GITVERSION)\" + -$(CC) -pg -o bin/vplanet src/*.c -lm -O3 -DGITVERSION=\"$(GITVERSION)\" sanitize: - -gcc -g -fsanitize=address -o bin/vplanet src/*.c -lm -DGITVERSION=\"$(GITVERSION)\" + -$(CC) -g -fsanitize=address -o bin/vplanet src/*.c -lm -DGITVERSION=\"$(GITVERSION)\" test: - -gcc -o bin/vplanet src/*.c -lm -O3 -DGITVERSION=\"$(GITVERSION)\" + -$(CC) -o bin/vplanet src/*.c -lm -O3 -DGITVERSION=\"$(GITVERSION)\" -pytest --tb=short coverage: -rm -f ./gcov/*.gcda ./gcov/*.gcno ./.coverage -mkdir -p ./gcov - -cd gcov && gcc -coverage -o ./../bin/vplanet ./../src/*.c -lm + -cd gcov && $(CC) -coverage -o ./../bin/vplanet ./../src/*.c -lm -python -m pytest --tb=short tests --junitxml=./junit/test-results.xml -lcov --capture --directory ./gcov --output-file ./.coverage -genhtml ./.coverage --output-directory ./gcov/html @@ -82,8 +84,8 @@ docs: -make -C docs html && echo 'Documentation available at `docs/.build/html/index.html`.' shared: - -gcc ${GCC_FLAGS1} src/*.c - -gcc ${GCC_FLAGS2} -o bin/vplanetlib.so *.o -lc + -$(CC) ${GCC_FLAGS1} src/*.c + -$(CC) ${GCC_FLAGS2} -o bin/vplanetlib.so *.o -lc clean: rm -f bin/vplanet diff --git a/examples/ClimateLand/equatorial.in b/examples/ClimateLand/equatorial.in new file mode 100755 index 000000000..73e384558 --- /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 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/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) diff --git a/examples/ClimateLand/modern.in b/examples/ClimateLand/modern.in new file mode 100755 index 000000000..b93eeb957 --- /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 +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..02c19eead --- /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 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/examples/ClimateLand/random.in b/examples/ClimateLand/random.in new file mode 100755 index 000000000..3c28f276a --- /dev/null +++ b/examples/ClimateLand/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/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..c7423f090 --- /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.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/examples/ClimateLand/vpl.in b/examples/ClimateLand/vpl.in new file mode 100755 index 000000000..ec7a23a9b --- /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.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/galhabit.c b/src/galhabit.c index dae38ffd4..2a85a5943 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 */ @@ -584,15 +565,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 09ccee1ac..916257d76 100644 --- a/src/options.c +++ b/src/options.c @@ -1161,6 +1161,10 @@ void ReadInitialOptions(BODY **body, CONTROL *control, FILES *files, ReadVerbose(files, options, &control->Io.iVerbose, iFile); } + if (control->Io.iVerbose >= VERBERR) { + fprintf(stderr, "INFO: Running VPLanet %s.\n", control->sGitVersion); + } + /* Need units prior to any parameter read */ control->Units = malloc(files->iNumInputs * sizeof(UNITS)); @@ -2617,7 +2621,7 @@ void fvAllocateOutputArrays(char ****saMatch, char ***saOutput, int **baNeg, for (iIndex = 0; iIndex < iNumArgs; iIndex++) { (*saMatch)[iIndex] = - malloc(MAXARRAY * sizeof(char*)); // Could be this many matches + malloc(MAXARRAY * sizeof(char *)); // Could be this many matches for (iMatch = 0; iMatch < MAXARRAY; iMatch++) { (*saMatch)[iIndex][iMatch] = NULL; } @@ -3243,6 +3247,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, @@ -4509,6 +4531,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..3a695bac7 100644 --- a/src/poise.c +++ b/src/poise.c @@ -140,7 +140,7 @@ void ReadFileOrbitOblData(BODY *body, CONTROL *control, FILES *files, AddOptionString(files->Infile[iFile].cIn, options->cName, cTmp, &lTmp, control->Io.iVerbose); - body[iFile - 1].sFileOrbitOblData=NULL; + body[iFile - 1].sFileOrbitOblData = NULL; if (lTmp >= 0) { /* Cannot exist in primary input file -- Each body has an output file */ NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp, @@ -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 = GEOGRAPHYUNIFORM; + } else if (!memcmp(sLower(cTmp), "m", 1)) { + body[iFile - 1].iGeography = GEOGRAPHYMODERN; + } else if (!memcmp(sLower(cTmp), "r", 1)) { + body[iFile - 1].iGeography = GEOGRAPHYRANDOM; + } else if (!memcmp(sLower(cTmp), "p", 1)) { + body[iFile - 1].iGeography = GEOGRAPHYPOLAR; + } else if (!memcmp(sLower(cTmp), "e", 1)) { + body[iFile - 1].iGeography = GEOGRAPHYEQUATORIAL; } 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,26 @@ 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) { @@ -933,6 +959,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 */ @@ -1126,8 +1202,9 @@ void ReadMinIceSheetHeight(BODY *body, CONTROL *control, FILES *files, void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fvFormattedString(&options[OPT_LATCELLNUM].cName, "iLatCellNum"); - fvFormattedString(&options[OPT_LATCELLNUM].cDescr, "Number of latitude cells used in" - " climate model"); + fvFormattedString(&options[OPT_LATCELLNUM].cDescr, + "Number of latitude cells used in" + " climate model"); fvFormattedString(&options[OPT_LATCELLNUM].cDefault, "50"); options[OPT_LATCELLNUM].dDefault = 50; options[OPT_LATCELLNUM].iType = 1; @@ -1135,7 +1212,8 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_LATCELLNUM] = &ReadLatCellNum; fvFormattedString(&options[OPT_READORBITOBLDATA].cName, "bReadOrbitOblData"); - fvFormattedString(&options[OPT_READORBITOBLDATA].cDescr, "Read in orbital and obliquity \ + fvFormattedString(&options[OPT_READORBITOBLDATA].cDescr, + "Read in orbital and obliquity \ data and use with poise"); fvFormattedString(&options[OPT_READORBITOBLDATA].cDefault, "0"); options[OPT_READORBITOBLDATA].dDefault = 0; @@ -1144,14 +1222,16 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_READORBITOBLDATA] = &ReadOrbitOblData; fvFormattedString(&options[OPT_FILEORBITOBLDATA].cName, "sFileOrbitOblData"); - fvFormattedString(&options[OPT_FILEORBITOBLDATA].cDescr, "Name of file containing orbit\ + fvFormattedString(&options[OPT_FILEORBITOBLDATA].cDescr, + "Name of file containing orbit\ and obliquity time series"); fvFormattedString(&options[OPT_FILEORBITOBLDATA].cDefault, "Obl_data.txt"); options[OPT_FILEORBITOBLDATA].iType = 3; fnRead[OPT_FILEORBITOBLDATA] = &ReadFileOrbitOblData; fvFormattedString(&options[OPT_PLANCKA].cName, "dPlanckA"); - fvFormattedString(&options[OPT_PLANCKA].cDescr, "Constant 'A' used in OLR calculation"); + fvFormattedString(&options[OPT_PLANCKA].cDescr, + "Constant 'A' used in OLR calculation"); fvFormattedString(&options[OPT_PLANCKA].cDefault, "203.3"); fvFormattedString(&options[OPT_PLANCKA].cDimension, "nd"); options[OPT_PLANCKA].dDefault = 203.3; @@ -1161,7 +1241,7 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fvFormattedString(&options[OPT_PLANCKB].cName, "dPlanckB"); fvFormattedString(&options[OPT_PLANCKB].cDescr, "Sensitivity 'B' used in" - " OLR calculation"); + " OLR calculation"); fvFormattedString(&options[OPT_PLANCKB].cDefault, "2.09"); fvFormattedString(&options[OPT_PLANCKB].cDimension, "nd"); options[OPT_PLANCKB].dDefault = 2.09; @@ -1179,7 +1259,8 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_ICEALBEDO] = &ReadIceAlbedo; fvFormattedString(&options[OPT_SURFALBEDO].cName, "dSurfAlbedo"); - fvFormattedString(&options[OPT_SURFALBEDO].cDescr, "Albedo of (ice-free) surface"); + fvFormattedString(&options[OPT_SURFALBEDO].cDescr, + "Albedo of (ice-free) surface"); fvFormattedString(&options[OPT_SURFALBEDO].cDefault, "0.3"); fvFormattedString(&options[OPT_SURFALBEDO].cDimension, "nd"); options[OPT_SURFALBEDO].dDefault = 0.3; @@ -1188,8 +1269,9 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_SURFALBEDO] = &ReadSurfAlbedo; fvFormattedString(&options[OPT_TGLOBALEST].cName, "dTGlobalInit"); - fvFormattedString(&options[OPT_TGLOBALEST].cDescr, "Estimate of initial global" - " temperature"); + fvFormattedString(&options[OPT_TGLOBALEST].cDescr, + "Estimate of initial global" + " temperature"); fvFormattedString(&options[OPT_TGLOBALEST].cDefault, "14.85"); fvFormattedString(&options[OPT_TGLOBALEST].cDimension, "temperature"); options[OPT_TGLOBALEST].dDefault = 14.85; @@ -1198,7 +1280,8 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_TGLOBALEST] = &ReadTGlobalInit; fvFormattedString(&options[OPT_PCO2].cName, "dpCO2"); - fvFormattedString(&options[OPT_PCO2].cDescr, "Partial pressure of CO2 in atmosphere"); + fvFormattedString(&options[OPT_PCO2].cDescr, + "Partial pressure of CO2 in atmosphere"); fvFormattedString(&options[OPT_PCO2].cDefault, "3.3e-4"); fvFormattedString(&options[OPT_PCO2].cDimension, "nd"); options[OPT_PCO2].dDefault = 3.3e-4; @@ -1207,8 +1290,9 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_PCO2] = &ReadPCO2; fvFormattedString(&options[OPT_CALCAB].cName, "bCalcAB"); - fvFormattedString(&options[OPT_CALCAB].cDescr, "Calculate A and B in OLR function, from" - " (T & pCO2)"); + fvFormattedString(&options[OPT_CALCAB].cDescr, + "Calculate A and B in OLR function, from" + " (T & pCO2)"); fvFormattedString(&options[OPT_CALCAB].cDefault, "0"); fvFormattedString(&options[OPT_CALCAB].cDimension, "nd"); options[OPT_CALCAB].dDefault = 0; @@ -1217,7 +1301,8 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_CALCAB] = &ReadCalcAB; fvFormattedString(&options[OPT_DIFFUSION].cName, "dDiffusion"); - fvFormattedString(&options[OPT_DIFFUSION].cDescr, "Heat diffusion coefficient"); + fvFormattedString(&options[OPT_DIFFUSION].cDescr, + "Heat diffusion coefficient"); fvFormattedString(&options[OPT_DIFFUSION].cDefault, "0.44"); fvFormattedString(&options[OPT_DIFFUSION].cDimension, "nd"); options[OPT_DIFFUSION].dDefault = 0.44; @@ -1235,7 +1320,7 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fvFormattedString(&options[OPT_COLDSTART].cName, "bColdStart"); fvFormattedString(&options[OPT_COLDSTART].cDescr, "Start from snowball Earth" - " conditions"); + " conditions"); fvFormattedString(&options[OPT_COLDSTART].cDefault, "0"); options[OPT_COLDSTART].dDefault = 0; options[OPT_COLDSTART].iType = 0; @@ -1244,7 +1329,7 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fvFormattedString(&options[OPT_FIXICELAT].cName, "dFixIceLat"); fvFormattedString(&options[OPT_FIXICELAT].cDescr, - "Force ice cap latitude to this value"); + "Force ice cap latitude to this value"); fvFormattedString(&options[OPT_FIXICELAT].cDefault, "None"); fvFormattedString(&options[OPT_FIXICELAT].cDimension, "nd"); options[OPT_FIXICELAT].dDefault = 0; @@ -1253,7 +1338,8 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_FIXICELAT] = &ReadFixIceLat; fvFormattedString(&options[OPT_ALBEDOZA].cName, "bAlbedoZA"); - fvFormattedString(&options[OPT_ALBEDOZA].cDescr, "Use albedo based on zenith angle"); + fvFormattedString(&options[OPT_ALBEDOZA].cDescr, + "Use albedo based on zenith angle"); fvFormattedString(&options[OPT_ALBEDOZA].cDefault, "0"); options[OPT_ALBEDOZA].dDefault = 0; options[OPT_ALBEDOZA].iType = 0; @@ -1261,8 +1347,9 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_ALBEDOZA] = &ReadAlbedoZA; fvFormattedString(&options[OPT_MEPDIFF].cName, "bMEPDiff"); - fvFormattedString(&options[OPT_MEPDIFF].cDescr, "Calculate diffusion from max entropy" - " production (D=B/4)?"); + fvFormattedString(&options[OPT_MEPDIFF].cDescr, + "Calculate diffusion from max entropy" + " production (D=B/4)?"); fvFormattedString(&options[OPT_MEPDIFF].cDefault, "0"); options[OPT_MEPDIFF].dDefault = 0; options[OPT_MEPDIFF].iType = 0; @@ -1270,8 +1357,9 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_MEPDIFF] = &ReadMEPDiff; fvFormattedString(&options[OPT_HEATCAPANN].cName, "dHeatCapAnn"); - fvFormattedString(&options[OPT_HEATCAPANN].cDescr, "Surface heat capacity in annual" - " model"); + fvFormattedString(&options[OPT_HEATCAPANN].cDescr, + "Surface heat capacity in annual" + " model"); fvFormattedString(&options[OPT_HEATCAPANN].cDefault, "0.2"); fvFormattedString(&options[OPT_HEATCAPANN].cDimension, "energy/temperature"); options[OPT_HEATCAPANN].dDefault = 0.2; @@ -1280,8 +1368,9 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_HEATCAPANN] = &ReadHeatCapAnn; fvFormattedString(&options[OPT_ICEDEPRATE].cName, "dIceDepRate"); - fvFormattedString(&options[OPT_ICEDEPRATE].cDescr, "Deposition rate of ice/snow to form" - " ice sheets"); + fvFormattedString(&options[OPT_ICEDEPRATE].cDescr, + "Deposition rate of ice/snow to form" + " ice sheets"); fvFormattedString(&options[OPT_ICEDEPRATE].cDefault, "2.9e-5"); fvFormattedString(&options[OPT_ICEDEPRATE].cDimension, "length/time"); options[OPT_ICEDEPRATE].dDefault = 2.9e-5; @@ -1298,7 +1387,8 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_ICESHEETS] = &ReadIceSheets; fvFormattedString(&options[OPT_INITICELAT].cName, "dInitIceLat"); - fvFormattedString(&options[OPT_INITICELAT].cDescr, "Sets initial ice sheet latitude"); + fvFormattedString(&options[OPT_INITICELAT].cDescr, + "Sets initial ice sheet latitude"); fvFormattedString(&options[OPT_INITICELAT].cDefault, "90"); fvFormattedString(&options[OPT_INITICELAT].cDimension, "angle"); options[OPT_INITICELAT].dDefault = 90.0; @@ -1307,7 +1397,8 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_INITICELAT] = &ReadInitIceLat; fvFormattedString(&options[OPT_INITICEHEIGHT].cName, "dInitIceHeight"); - fvFormattedString(&options[OPT_INITICEHEIGHT].cDescr, "Sets initial ice sheet height"); + fvFormattedString(&options[OPT_INITICEHEIGHT].cDescr, + "Sets initial ice sheet height"); fvFormattedString(&options[OPT_INITICEHEIGHT].cDefault, "50"); // 50 meters fvFormattedString(&options[OPT_INITICEHEIGHT].cDimension, "length"); options[OPT_INITICEHEIGHT].dDefault = 50.0; @@ -1316,7 +1407,8 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_INITICEHEIGHT] = &ReadInitIceHeight; fvFormattedString(&options[OPT_CLIMATEMODEL].cName, "sClimateModel"); - fvFormattedString(&options[OPT_CLIMATEMODEL].cDescr, "Use annual or seasonal model"); + fvFormattedString(&options[OPT_CLIMATEMODEL].cDescr, + "Use annual or seasonal model"); fvFormattedString(&options[OPT_CLIMATEMODEL].cDefault, "ann"); options[OPT_CLIMATEMODEL].dDefault = ANN; options[OPT_CLIMATEMODEL].iType = 3; @@ -1324,7 +1416,8 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_CLIMATEMODEL] = &ReadClimateModel; fvFormattedString(&options[OPT_OLRMODEL].cName, "iOLRModel"); - fvFormattedString(&options[OPT_OLRMODEL].cDescr, "Outgoing longwave rad model"); + fvFormattedString(&options[OPT_OLRMODEL].cDescr, + "Outgoing longwave rad model"); fvFormattedString(&options[OPT_OLRMODEL].cDefault, "sms09"); options[OPT_OLRMODEL].dDefault = SMS09; options[OPT_OLRMODEL].iType = 1; @@ -1332,8 +1425,9 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_OLRMODEL] = &ReadOLRModel; fvFormattedString(&options[OPT_SKIPSEASENABLED].cName, "bSkipSeasEnabled"); - fvFormattedString(&options[OPT_SKIPSEASENABLED].cDescr, "Run annual before seasonal and" - " allow skip seas?"); + fvFormattedString(&options[OPT_SKIPSEASENABLED].cDescr, + "Run annual before seasonal and" + " allow skip seas?"); fvFormattedString(&options[OPT_SKIPSEASENABLED].cDefault, "0"); options[OPT_SKIPSEASENABLED].dDefault = 0; options[OPT_SKIPSEASENABLED].iType = 0; @@ -1341,9 +1435,11 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_SKIPSEASENABLED] = &ReadSkipSeasEnabled; fvFormattedString(&options[OPT_HEATCAPLAND].cName, "dHeatCapLand"); - fvFormattedString(&options[OPT_HEATCAPLAND].cDescr, "Land heat capacity in seasonal" - " model"); - fvFormattedString(&options[OPT_HEATCAPLAND].cDefault, "1.42e7"); // XXX What units? + fvFormattedString(&options[OPT_HEATCAPLAND].cDescr, + "Land heat capacity in seasonal" + " model"); + fvFormattedString(&options[OPT_HEATCAPLAND].cDefault, + "1.42e7"); // XXX What units? fvFormattedString(&options[OPT_HEATCAPLAND].cDimension, "energy/temperature"); options[OPT_HEATCAPLAND].dDefault = 1.42e7; options[OPT_HEATCAPLAND].iType = 2; @@ -1351,18 +1447,22 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_HEATCAPLAND] = &ReadHeatCapLand; fvFormattedString(&options[OPT_HEATCAPWATER].cName, "dHeatCapWater"); - fvFormattedString(&options[OPT_HEATCAPWATER].cDescr, "Water heat capacity per meter in" - " seasonal model"); - fvFormattedString(&options[OPT_HEATCAPWATER].cDefault, "4.2e6"); // XXX What units - fvFormattedString(&options[OPT_HEATCAPWATER].cDimension, "energy/temperature"); + fvFormattedString(&options[OPT_HEATCAPWATER].cDescr, + "Water heat capacity per meter in" + " seasonal model"); + fvFormattedString(&options[OPT_HEATCAPWATER].cDefault, + "4.2e6"); // XXX What units + fvFormattedString(&options[OPT_HEATCAPWATER].cDimension, + "energy/temperature"); options[OPT_HEATCAPWATER].dDefault = 4.2e6; options[OPT_HEATCAPWATER].iType = 2; options[OPT_HEATCAPWATER].bMultiFile = 1; fnRead[OPT_HEATCAPWATER] = &ReadHeatCapWater; fvFormattedString(&options[OPT_MIXINGDEPTH].cName, "dMixingDepth"); - fvFormattedString(&options[OPT_MIXINGDEPTH].cDescr, "Mixing depth of ocean in seasonal" - " model"); + fvFormattedString(&options[OPT_MIXINGDEPTH].cDescr, + "Mixing depth of ocean in seasonal" + " model"); fvFormattedString(&options[OPT_MIXINGDEPTH].cDefault, "70"); // meters fvFormattedString(&options[OPT_MIXINGDEPTH].cDimension, "length"); options[OPT_MIXINGDEPTH].dDefault = 70.; @@ -1371,8 +1471,9 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_MIXINGDEPTH] = &ReadMixingDepth; fvFormattedString(&options[OPT_FRZTSEAICE].cName, "dFrzTSeaIce"); - fvFormattedString(&options[OPT_FRZTSEAICE].cDescr, "Temp of sea ice formation in" - " seasonal model"); + fvFormattedString(&options[OPT_FRZTSEAICE].cDescr, + "Temp of sea ice formation in" + " seasonal model"); fvFormattedString(&options[OPT_FRZTSEAICE].cDefault, "-2 deg C"); fvFormattedString(&options[OPT_FRZTSEAICE].cDimension, "temperature"); options[OPT_FRZTSEAICE].dDefault = -2.; @@ -1381,8 +1482,9 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_FRZTSEAICE] = &ReadFrzTSeaIce; fvFormattedString(&options[OPT_NULANDWATER].cName, "dNuLandWater"); - fvFormattedString(&options[OPT_NULANDWATER].cDescr, "Coefficient of land-ocean heat" - " flux"); + fvFormattedString(&options[OPT_NULANDWATER].cDescr, + "Coefficient of land-ocean heat" + " flux"); fvFormattedString(&options[OPT_NULANDWATER].cDefault, "0.81"); fvFormattedString(&options[OPT_NULANDWATER].cDimension, "energy/length^2"); options[OPT_NULANDWATER].dDefault = 0.81; @@ -1391,8 +1493,9 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_NULANDWATER] = &ReadNuLandWater; fvFormattedString(&options[OPT_LANDFRAC].cName, "dLandFrac"); - fvFormattedString(&options[OPT_LANDFRAC].cDescr, "Fraction of land on the planetary" - " surface"); + fvFormattedString(&options[OPT_LANDFRAC].cDescr, + "Fraction of land on the planetary" + " surface"); fvFormattedString(&options[OPT_LANDFRAC].cDefault, "0.34"); fvFormattedString(&options[OPT_LANDFRAC].cDimension, "nd"); options[OPT_LANDFRAC].dDefault = 0.34; @@ -1400,9 +1503,36 @@ 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_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, + "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"); + fvFormattedString(&options[OPT_NSTEPINYEAR].cDescr, + "Number of time-steps/year in" + " seasonal model"); fvFormattedString(&options[OPT_NSTEPINYEAR].cDefault, "60"); options[OPT_NSTEPINYEAR].dDefault = 60; options[OPT_NSTEPINYEAR].iType = 1; @@ -1410,8 +1540,9 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_NSTEPINYEAR] = &ReadNStepInYear; fvFormattedString(&options[OPT_NUMYEARS].cName, "iNumYears"); - fvFormattedString(&options[OPT_NUMYEARS].cDescr, "Number of years to run seasonal" - " model"); + fvFormattedString(&options[OPT_NUMYEARS].cDescr, + "Number of years to run seasonal" + " model"); fvFormattedString(&options[OPT_NUMYEARS].cDefault, "10"); options[OPT_NUMYEARS].dDefault = 10; options[OPT_NUMYEARS].iType = 1; @@ -1419,8 +1550,9 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_NUMYEARS] = &ReadNumYears; fvFormattedString(&options[OPT_SEAICEMODEL].cName, "bSeaIceModel"); - fvFormattedString(&options[OPT_SEAICEMODEL].cDescr, "model sea ice dynamics and heat" - " flow?"); + fvFormattedString(&options[OPT_SEAICEMODEL].cDescr, + "model sea ice dynamics and heat" + " flow?"); fvFormattedString(&options[OPT_SEAICEMODEL].cDefault, "1"); options[OPT_SEAICEMODEL].dDefault = 1; options[OPT_SEAICEMODEL].iType = 0; @@ -1428,7 +1560,8 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_SEAICEMODEL] = &ReadSeaIceModel; fvFormattedString(&options[OPT_ICECONDUCT].cName, "dSeaIceConduct"); - fvFormattedString(&options[OPT_ICECONDUCT].cDescr, "Heat conductivity of sea ice"); + fvFormattedString(&options[OPT_ICECONDUCT].cDescr, + "Heat conductivity of sea ice"); fvFormattedString(&options[OPT_ICECONDUCT].cDefault, "2"); fvFormattedString(&options[OPT_ICECONDUCT].cDimension, "nd"); options[OPT_ICECONDUCT].dDefault = 2.; @@ -1455,8 +1588,9 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_ALBEDOWATER] = &ReadAlbedoWater; fvFormattedString(&options[OPT_ICEDT].cName, "iIceDt"); - fvFormattedString(&options[OPT_ICEDT].cDescr, "Minimum ice sheet timestep (unit orbital" - " period)"); + fvFormattedString(&options[OPT_ICEDT].cDescr, + "Minimum ice sheet timestep (unit orbital" + " period)"); fvFormattedString(&options[OPT_ICEDT].cDefault, "5"); options[OPT_ICEDT].dDefault = 5; options[OPT_ICEDT].iType = 1; @@ -1464,9 +1598,10 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_ICEDT] = &ReadIceDt; fvFormattedString(&options[OPT_RERUNSEAS].cName, "iReRunSeas"); - fvFormattedString(&options[OPT_RERUNSEAS].cDescr, "how often to rerun seasonal in ice" - " sheet model, in number of orbital" - " periods"); + fvFormattedString(&options[OPT_RERUNSEAS].cDescr, + "how often to rerun seasonal in ice" + " sheet model, in number of orbital" + " periods"); fvFormattedString(&options[OPT_RERUNSEAS].cDefault, "500"); options[OPT_RERUNSEAS].dDefault = 5; options[OPT_RERUNSEAS].iType = 1; @@ -1474,16 +1609,49 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_RERUNSEAS] = &ReadReRunSeas; 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].cDescr, + "Type of land distribution"); + fvFormattedString(&options[OPT_GEOGRAPHY].cDefault, "uniform"); + options[OPT_GEOGRAPHY].dDefault = GEOGRAPHYUNIFORM; 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" - " parameters"); + fvFormattedString(&options[OPT_SEASOUTPUTTIME].cDescr, + "Output interval for seasonal" + " parameters"); fvFormattedString(&options[OPT_SEASOUTPUTTIME].cDefault, "0"); fvFormattedString(&options[OPT_SEASOUTPUTTIME].cDimension, "nd"); options[OPT_SEASOUTPUTTIME].dDefault = 0; @@ -1493,7 +1661,7 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fvFormattedString(&options[OPT_FORCEOBLIQ].cName, "bForceObliq"); fvFormattedString(&options[OPT_FORCEOBLIQ].cDescr, "Force obliquity to evolve" - " sinusoidally?"); + " sinusoidally?"); fvFormattedString(&options[OPT_FORCEOBLIQ].cDefault, "0"); options[OPT_FORCEOBLIQ].dDefault = 0; options[OPT_FORCEOBLIQ].iType = 0; @@ -1501,8 +1669,9 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_FORCEOBLIQ] = &ReadForceObliq; fvFormattedString(&options[OPT_DIFFROT].cName, "bDiffRot"); - fvFormattedString(&options[OPT_DIFFROT].cDescr, "Adjust heat diffusion for rotation" - " rate?"); + fvFormattedString(&options[OPT_DIFFROT].cDescr, + "Adjust heat diffusion for rotation" + " rate?"); fvFormattedString(&options[OPT_DIFFROT].cDefault, "0"); options[OPT_DIFFROT].dDefault = 0; options[OPT_DIFFROT].iType = 0; @@ -1510,7 +1679,8 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_DIFFROT] = &ReadDiffRot; fvFormattedString(&options[OPT_OBLIQAMP].cName, "dObliqAmp"); - fvFormattedString(&options[OPT_OBLIQAMP].cDescr, "Amplitude of forced obliquity oscill"); + fvFormattedString(&options[OPT_OBLIQAMP].cDescr, + "Amplitude of forced obliquity oscill"); fvFormattedString(&options[OPT_OBLIQAMP].cDefault, "50 deg"); fvFormattedString(&options[OPT_OBLIQAMP].cDimension, "angle"); options[OPT_OBLIQAMP].dDefault = 50; @@ -1519,8 +1689,10 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_OBLIQAMP] = &ReadObliqAmp; fvFormattedString(&options[OPT_OBLIQPER].cName, "dObliqPer"); - fvFormattedString(&options[OPT_OBLIQPER].cDescr, "Period of forced obliquity oscill"); - fvFormattedString(&options[OPT_OBLIQPER].cDefault, "50000"); // XXX What units? + fvFormattedString(&options[OPT_OBLIQPER].cDescr, + "Period of forced obliquity oscill"); + fvFormattedString(&options[OPT_OBLIQPER].cDefault, + "50000"); // XXX What units? fvFormattedString(&options[OPT_OBLIQPER].cDimension, "time"); options[OPT_OBLIQPER].dDefault = 50000; options[OPT_OBLIQPER].iType = 2; @@ -1528,8 +1700,9 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_OBLIQPER] = &ReadObliqPer; fvFormattedString(&options[OPT_FORCEECC].cName, "bForceEcc"); - fvFormattedString(&options[OPT_FORCEECC].cDescr, "Force Eccentricity to evolve" - " sinusoidally?"); + fvFormattedString(&options[OPT_FORCEECC].cDescr, + "Force Eccentricity to evolve" + " sinusoidally?"); fvFormattedString(&options[OPT_FORCEECC].cDefault, "0"); options[OPT_FORCEECC].dDefault = 0; options[OPT_FORCEECC].iType = 0; @@ -1537,8 +1710,9 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_FORCEECC] = &ReadForceEcc; fvFormattedString(&options[OPT_ECCAMP].cName, "dEccAmp"); - fvFormattedString(&options[OPT_ECCAMP].cDescr, "Amplitude of forced eccentricity" - " oscill"); + fvFormattedString(&options[OPT_ECCAMP].cDescr, + "Amplitude of forced eccentricity" + " oscill"); fvFormattedString(&options[OPT_ECCAMP].cDefault, "0.1"); fvFormattedString(&options[OPT_ECCAMP].cDimension, "nd"); options[OPT_ECCAMP].dDefault = 0.1; @@ -1547,7 +1721,8 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_ECCAMP] = &ReadEccAmp; fvFormattedString(&options[OPT_ECCPER].cName, "dEccPer"); - fvFormattedString(&options[OPT_ECCPER].cDescr, "Period of forced eccentricity oscill"); + fvFormattedString(&options[OPT_ECCPER].cDescr, + "Period of forced eccentricity oscill"); fvFormattedString(&options[OPT_ECCPER].cDefault, "50000"); // !!! What units? fvFormattedString(&options[OPT_ECCPER].cDimension, "time"); options[OPT_ECCPER].dDefault = 50000; @@ -1557,7 +1732,7 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fvFormattedString(&options[OPT_ACCUMODE].cName, "bAccuracyMode"); fvFormattedString(&options[OPT_ACCUMODE].cDescr, - "Re-invert matrix every EBM time step?"); + "Re-invert matrix every EBM time step?"); fvFormattedString(&options[OPT_ACCUMODE].cDefault, "0"); options[OPT_ACCUMODE].dDefault = 0; options[OPT_ACCUMODE].iType = 0; @@ -1565,8 +1740,9 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_ACCUMODE] = &ReadAccuracyMode; fvFormattedString(&options[OPT_ELEVFB].cName, "bElevFB"); - fvFormattedString(&options[OPT_ELEVFB].cDescr, "Use elevation feedback for ice sheet" - " ablation?"); + fvFormattedString(&options[OPT_ELEVFB].cDescr, + "Use elevation feedback for ice sheet" + " ablation?"); fvFormattedString(&options[OPT_ELEVFB].cDefault, "0"); options[OPT_ELEVFB].dDefault = 0; options[OPT_ELEVFB].iType = 0; @@ -1574,8 +1750,9 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_ELEVFB] = &ReadElevFB; fvFormattedString(&options[OPT_LAPSER].cName, "dLapseR"); - fvFormattedString(&options[OPT_LAPSER].cDescr, "Dry adiabatic lapse rate (for elev" - " feedback)"); + fvFormattedString(&options[OPT_LAPSER].cDescr, + "Dry adiabatic lapse rate (for elev" + " feedback)"); fvFormattedString(&options[OPT_LAPSER].cDefault, "9.8e-3 C/m"); fvFormattedString(&options[OPT_LAPSER].cDimension, "temperature/length"); options[OPT_LAPSER].dDefault = 9.8e-3; @@ -1584,8 +1761,9 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_LAPSER] = &ReadLapseR; fvFormattedString(&options[OPT_REFHEIGHT].cName, "dRefHeight"); - fvFormattedString(&options[OPT_REFHEIGHT].cDescr, "Reference height of atmos temp (for" - " elev feedback)"); + fvFormattedString(&options[OPT_REFHEIGHT].cDescr, + "Reference height of atmos temp (for" + " elev feedback)"); fvFormattedString(&options[OPT_REFHEIGHT].cDefault, "1000 m"); fvFormattedString(&options[OPT_REFHEIGHT].cDimension, "length"); options[OPT_REFHEIGHT].dDefault = 1000.0; @@ -1603,7 +1781,8 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_ABLATEFF] = &ReadAblateFF; fvFormattedString(&options[OPT_SPINUPTOL].cName, "dSpinUpTol"); - fvFormattedString(&options[OPT_SPINUPTOL].cDescr, "Tolerance for spin up phase"); + fvFormattedString(&options[OPT_SPINUPTOL].cDescr, + "Tolerance for spin up phase"); fvFormattedString(&options[OPT_SPINUPTOL].cDefault, "0.1 deg C"); fvFormattedString(&options[OPT_SPINUPTOL].cDimension, "temperature"); options[OPT_SPINUPTOL].dDefault = 0.1; @@ -1612,10 +1791,12 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { fnRead[OPT_SPINUPTOL] = &ReadSpinUpTol; fvFormattedString(&options[OPT_MINICEHEIGHT].cName, "dMinIceSheetHeight"); - fvFormattedString(&options[OPT_MINICEHEIGHT].cDescr, "Minimum ice sheet height for a" - " latitude to be considered" - " ice-covered"); - fvFormattedString(&options[OPT_MINICEHEIGHT].cDefault, "0.001"); // What units? XXX + fvFormattedString(&options[OPT_MINICEHEIGHT].cDescr, + "Minimum ice sheet height for a" + " latitude to be considered" + " ice-covered"); + fvFormattedString(&options[OPT_MINICEHEIGHT].cDefault, + "0.001"); // What units? XXX fvFormattedString(&options[OPT_MINICEHEIGHT].cDimension, "length"); options[OPT_MINICEHEIGHT].dDefault = 0.001; options[OPT_MINICEHEIGHT].iType = 2; @@ -1623,15 +1804,16 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) { options[OPT_MINICEHEIGHT].dNeg = 1; // Convert to SI fvFormattedString(&options[OPT_MINICEHEIGHT].cNeg, "meters"); fnRead[OPT_MINICEHEIGHT] = &ReadMinIceSheetHeight; - fvFormattedString(&options[OPT_MINICEHEIGHT].cLongDescr, - "The minimum thickness of permanent ice in a latitude bin for it to " - "be\n" - "labeled ice-covered. In some cases, such as rapid thawing, a " - "latituden\n" - "can have a very low value of ice height, primarily for numerical\n" - "reasons. Parameter forces unphysically small values of the ice " - "height\n" - "to be ignored.\n"); + fvFormattedString( + &options[OPT_MINICEHEIGHT].cLongDescr, + "The minimum thickness of permanent ice in a latitude bin for it to " + "be\n" + "labeled ice-covered. In some cases, such as rapid thawing, a " + "latituden\n" + "can have a very low value of ice height, primarily for numerical\n" + "reasons. Parameter forces unphysically small values of the ice " + "height\n" + "to be ignored.\n"); } void ReadOptionsPoise(BODY *body, CONTROL *control, FILES *files, @@ -1826,10 +2008,12 @@ void VerifyOrbitOblData(BODY *body, CONTROL *control, OPTIONS *options, iLine = 0; while (feof(fileorb) == 0) { - if (fscanf(fileorb, "%lf %lf %lf %lf %lf %lf %lf", &dttmp, &datmp, &detmp, - &daptmp, &dlatmp, &dobltmp, &dprecatmp) != 7) { - fprintf(stderr,"ERROR: Incorrect number of columns in orbit-obliquity file."); - exit(EXIT_INPUT); + if (fscanf(fileorb, "%lf %lf %lf %lf %lf %lf %lf", &dttmp, &datmp, + &detmp, &daptmp, &dlatmp, &dobltmp, &dprecatmp) != 7) { + fprintf( + stderr, + "ERROR: Incorrect number of columns in orbit-obliquity file."); + exit(EXIT_INPUT); } body[iBody].daTimeSeries[iLine] = @@ -1901,50 +2085,157 @@ 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].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)); - for (iLats = 0; iLats < body[iBody].iNumLats; iLats++) { - SinLat = (-1.0 + dDelta_x / 2.) + iLats * dDelta_x; - body[iBody].daLats[iLats] = asin(SinLat); + 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); + + 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); } } -void InitializeLandWater(BODY *body, int iBody) { +void fvInitializeGeographyUniform(BODY *body, int iBody) { int iLat; - body[iBody].daLandFrac = malloc(body[iBody].iNumLats * sizeof(double)); - body[iBody].daWaterFrac = malloc(body[iBody].iNumLats * sizeof(double)); + 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]; + } +} - if (body[iBody].iGeography == UNIFORM3) { - 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 fvInitializeGeographyModern(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; } - } else if (body[iBody].iGeography == MODERN) { - 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]; + body[iBody].daWaterFrac[iLat] = 1.0 - body[iBody].daLandFrac[iLat]; + } +} + +void fvInitializeGeographyRandom(BODY *body, int iBody) { + int iLat; + double dOffset, dLandFrac, dLandFracNormFactor; + + for (iLat = 0; iLat < body[iBody].iNumLats; iLat++) { + 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; + } + } + + 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 fvInitializeGeographyPolar(BODY *body, int iBody) { + int iLat; + + 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 fvInitializeGeographyEquatorial(BODY *body, int iBody) { + int iLat; + + 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 fvWriteGeographyFile(BODY *body, int iBody) { + int iLat; + FILE *fp; + 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 %.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); +} + +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) { + fvInitializeGeographyUniform(body, iBody); + } else if (body[iBody].iGeography == GEOGRAPHYMODERN) { + fvInitializeGeographyModern(body, iBody); + } else if (body[iBody].iGeography == GEOGRAPHYRANDOM) { + fvInitializeGeographyRandom(body, iBody); + } else if (body[iBody].iGeography == GEOGRAPHYPOLAR) { + fvInitializeGeographyPolar(body, iBody); + } else if (body[iBody].iGeography == GEOGRAPHYEQUATORIAL) { + fvInitializeGeographyEquatorial(body, iBody); + } + + fvWriteGeographyFile(body, iBody); +} void DampTemp(BODY *body, double dTGlobalTmp, int iBody) { int iLat; @@ -2242,7 +2533,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 / @@ -2624,6 +2915,146 @@ void NullPoiseDerivatives(BODY *body, EVOLVE *evolve, UPDATE *update, // Nothing here, because entire climate simulation is ran in ForceBehavior } +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]); +} + +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_LANDFRACMEAN].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]); +} + +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_LANDFRACAMP].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]); +} + +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) { + + 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, + "\nWARNING: %s is set to \"random\", but %s is not set.\n\n", + options[OPT_GEOGRAPHY].cName, options[OPT_RANDSEED].cName); + } + } + 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); + } + } +} + +void VerifyGeography(BODY *body, CONTROL *control, OPTIONS *options, + SYSTEM *system, int iBody) { + + 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); + } +} + void VerifyPoise(BODY *body, CONTROL *control, FILES *files, OPTIONS *options, OUTPUT *output, SYSTEM *system, UPDATE *update, int iBody, int iModule) { @@ -2642,6 +3073,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); @@ -3267,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, @@ -3332,7 +3759,7 @@ void WriteDailyInsol(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, UNITS *units, UPDATE *update, int iBody, double *dTmp, char **cUnit) { - char *cOut=NULL; + char *cOut = NULL; FILE *fp; int iLat, iDay; double dTime; @@ -3351,18 +3778,18 @@ void WriteDailyInsol(BODY *body, CONTROL *control, OUTPUT *output, if (dTime == 0) { - fvFormattedString(&cOut, "SeasonalClimateFiles/%s.%s.DailyInsol.0", system->cName, - body[iBody].cName); + fvFormattedString(&cOut, "SeasonalClimateFiles/%s.%s.DailyInsol.0", + system->cName, body[iBody].cName); } else if (dTime < 10000) { - fvFormattedString(&cOut, "SeasonalClimateFiles/%s.%s.DailyInsol.%.0f", system->cName, - body[iBody].cName, dTime); + fvFormattedString(&cOut, "SeasonalClimateFiles/%s.%s.DailyInsol.%.0f", + system->cName, body[iBody].cName, dTime); } else { - fvFormattedString(&cOut, "SeasonalClimateFiles/%s.%s.DailyInsol.%.2e", system->cName, - body[iBody].cName, dTime); + fvFormattedString(&cOut, "SeasonalClimateFiles/%s.%s.DailyInsol.%.2e", + system->cName, body[iBody].cName, dTime); } fp = fopen(cOut, "w"); @@ -3386,7 +3813,7 @@ void WritePlanckB(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, UNITS *units, UPDATE *update, int iBody, double *dTmp, char **cUnit) { - char *cOut=NULL; + char *cOut = NULL; FILE *fp; int iLat, iDay; double dTime; @@ -3405,18 +3832,18 @@ void WritePlanckB(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, if (dTime == 0) { - fvFormattedString(&cOut, "SeasonalClimateFiles/%s.%s.PlanckB.0", system->cName, - body[iBody].cName); + fvFormattedString(&cOut, "SeasonalClimateFiles/%s.%s.PlanckB.0", + system->cName, body[iBody].cName); } else if (dTime < 10000) { - fvFormattedString(&cOut, "SeasonalClimateFiles/%s.%s.PlanckB.%.0f", system->cName, - body[iBody].cName, dTime); + fvFormattedString(&cOut, "SeasonalClimateFiles/%s.%s.PlanckB.%.0f", + system->cName, body[iBody].cName, dTime); } else { - fvFormattedString(&cOut, "SeasonalClimateFiles/%s.%s.PlanckB.%.2e", system->cName, - body[iBody].cName, dTime); + fvFormattedString(&cOut, "SeasonalClimateFiles/%s.%s.PlanckB.%.2e", + system->cName, body[iBody].cName, dTime); } fp = fopen(cOut, "w"); @@ -3441,7 +3868,7 @@ void WriteSeasonalTemp(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, UNITS *units, UPDATE *update, int iBody, double *dTmp, char **cUnit) { - char *cOut=NULL; + char *cOut = NULL; FILE *fp; int iLat, iDay; double dTime; @@ -3460,18 +3887,18 @@ void WriteSeasonalTemp(BODY *body, CONTROL *control, OUTPUT *output, if (dTime == 0) { - fvFormattedString(&cOut, "SeasonalClimateFiles/%s.%s.SeasonalTemp.0", system->cName, - body[iBody].cName); + fvFormattedString(&cOut, "SeasonalClimateFiles/%s.%s.SeasonalTemp.0", + system->cName, body[iBody].cName); } else if (dTime < 10000) { - fvFormattedString(&cOut, "SeasonalClimateFiles/%s.%s.SeasonalTemp.%.0f", system->cName, - body[iBody].cName, dTime); + fvFormattedString(&cOut, "SeasonalClimateFiles/%s.%s.SeasonalTemp.%.0f", + system->cName, body[iBody].cName, dTime); } else { - fvFormattedString(&cOut, "SeasonalClimateFiles/%s.%s.SeasonalTemp.%.2e", system->cName, - body[iBody].cName, dTime); + fvFormattedString(&cOut, "SeasonalClimateFiles/%s.%s.SeasonalTemp.%.2e", + system->cName, body[iBody].cName, dTime); } fp = fopen(cOut, "w"); @@ -3498,7 +3925,7 @@ void WriteSeasonalFluxes(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, UNITS *units, UPDATE *update, int iBody, double *dTmp, char **cUnit) { - char *cOutM=NULL, *cOutI=NULL, *cOutO=NULL, *cOutD=NULL; + char *cOutM = NULL, *cOutI = NULL, *cOutO = NULL, *cOutD = NULL; FILE *fpM, *fpI, *fpO, *fpD; int iLat, iDay; double dTime; @@ -3517,36 +3944,36 @@ void WriteSeasonalFluxes(BODY *body, CONTROL *control, OUTPUT *output, if (dTime == 0) { - fvFormattedString(&cOutM, "SeasonalClimateFiles/%s.%s.SeasonalFMerid.0", system->cName, - body[iBody].cName); - fvFormattedString(&cOutI, "SeasonalClimateFiles/%s.%s.SeasonalFIn.0", system->cName, - body[iBody].cName); - fvFormattedString(&cOutO, "SeasonalClimateFiles/%s.%s.SeasonalFOut.0", system->cName, - body[iBody].cName); - fvFormattedString(&cOutD, "SeasonalClimateFiles/%s.%s.SeasonalDivF.0", system->cName, - body[iBody].cName); + fvFormattedString(&cOutM, "SeasonalClimateFiles/%s.%s.SeasonalFMerid.0", + system->cName, body[iBody].cName); + fvFormattedString(&cOutI, "SeasonalClimateFiles/%s.%s.SeasonalFIn.0", + system->cName, body[iBody].cName); + fvFormattedString(&cOutO, "SeasonalClimateFiles/%s.%s.SeasonalFOut.0", + system->cName, body[iBody].cName); + fvFormattedString(&cOutD, "SeasonalClimateFiles/%s.%s.SeasonalDivF.0", + system->cName, body[iBody].cName); } else if (dTime < 10000) { fvFormattedString(&cOutM, "SeasonalClimateFiles/%s.%s.SeasonalFMerid.%.0f", - system->cName, body[iBody].cName, dTime); - fvFormattedString(&cOutI, "SeasonalClimateFiles/%s.%s.SeasonalFIn.%.0f", system->cName, - body[iBody].cName, dTime); + system->cName, body[iBody].cName, dTime); + fvFormattedString(&cOutI, "SeasonalClimateFiles/%s.%s.SeasonalFIn.%.0f", + system->cName, body[iBody].cName, dTime); fvFormattedString(&cOutO, "SeasonalClimateFiles/%s.%s.SeasonalFOut.%.0f", - system->cName, body[iBody].cName, dTime); + system->cName, body[iBody].cName, dTime); fvFormattedString(&cOutD, "SeasonalClimateFiles/%s.%s.SeasonalDivF.%.0f", - system->cName, body[iBody].cName, dTime); + system->cName, body[iBody].cName, dTime); } else { fvFormattedString(&cOutM, "SeasonalClimateFiles/%s.%s.SeasonalFMerid.%.2e", - system->cName, body[iBody].cName, dTime); - fvFormattedString(&cOutI, "SeasonalClimateFiles/%s.%s.SeasonalFIn.%.2e", system->cName, - body[iBody].cName, dTime); + system->cName, body[iBody].cName, dTime); + fvFormattedString(&cOutI, "SeasonalClimateFiles/%s.%s.SeasonalFIn.%.2e", + system->cName, body[iBody].cName, dTime); fvFormattedString(&cOutO, "SeasonalClimateFiles/%s.%s.SeasonalFOut.%.2e", - system->cName, body[iBody].cName, dTime); + system->cName, body[iBody].cName, dTime); fvFormattedString(&cOutD, "SeasonalClimateFiles/%s.%s.SeasonalDivF.%.2e", - system->cName, body[iBody].cName, dTime); + system->cName, body[iBody].cName, dTime); } fpM = fopen(cOutM, "w"); @@ -3598,7 +4025,7 @@ void WriteSeasonalIceBalance(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, UNITS *units, UPDATE *update, int iBody, double *dTmp, char **cUnit) { - char *cOut=NULL; + char *cOut = NULL; FILE *fp; int iLat, iDay; double dTime; @@ -3618,17 +4045,19 @@ void WriteSeasonalIceBalance(BODY *body, CONTROL *control, OUTPUT *output, if (dTime == 0) { fvFormattedString(&cOut, "SeasonalClimateFiles/%s.%s.SeasonalIceBalance.0", - system->cName, body[iBody].cName); + system->cName, body[iBody].cName); } else if (dTime < 10000) { - fvFormattedString(&cOut, "SeasonalClimateFiles/%s.%s.SeasonalIceBalance.%.0f", - system->cName, body[iBody].cName, dTime); + fvFormattedString(&cOut, + "SeasonalClimateFiles/%s.%s.SeasonalIceBalance.%.0f", + system->cName, body[iBody].cName, dTime); } else { - fvFormattedString(&cOut, "SeasonalClimateFiles/%s.%s.SeasonalIceBalance.%.2e", - system->cName, body[iBody].cName, dTime); + fvFormattedString(&cOut, + "SeasonalClimateFiles/%s.%s.SeasonalIceBalance.%.2e", + system->cName, body[iBody].cName, dTime); } fp = fopen(cOut, "w"); @@ -3744,17 +4173,18 @@ void WriteIceMass(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system, if (body[iBody].bIceSheets) { *dTmp = body[iBody].daIceMass[body[iBody].iWriteLat]; - } else { *dTmp = 0.0; } if (output->bDoNeg[iBody]) { fvFormattedString(cUnit, output->cNeg); - } else { - //*dTmp /= fdUnitsMass(units->iMass)/pow(fdUnitsLength(units->iLength),2); - // fsUnitsEnergyFlux(units,cUnit); + fvFormattedString(cUnit,""); + /* XXX Need to fix units! + *dTmp /= fdUnitsMass(units->iMass)/pow(fdUnitsLength(units->iLength),2); + fsUnitsMass(units->iMass,cUnit); + */ } } @@ -3804,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 { @@ -3842,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 { @@ -3861,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 { @@ -3917,9 +4351,17 @@ 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, "Global mean temperature from POISE"); + fvFormattedString(&output[OUT_TGLOBAL].cDescr, + "Global mean temperature from POISE"); fvFormattedString(&output[OUT_TGLOBAL].cNeg, "Celsius"); output[OUT_TGLOBAL].bNeg = 1; // conversion is hardcoded in write function @@ -3930,21 +4372,23 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fvFormattedString(&output[OUT_ALBEDOGLOBAL].cName, "AlbedoGlobal"); fvFormattedString(&output[OUT_ALBEDOGLOBAL].cDescr, - "Global mean bond albedo from POISE"); + "Global mean bond albedo from POISE"); output[OUT_ALBEDOGLOBAL].bNeg = 0; output[OUT_ALBEDOGLOBAL].iNum = 1; output[OUT_ALBEDOGLOBAL].iModuleBit = POISE; fnWrite[OUT_ALBEDOGLOBAL] = &WriteAlbedoGlobal; fvFormattedString(&output[OUT_SNOWBALL].cName, "Snowball"); - fvFormattedString(&output[OUT_SNOWBALL].cDescr, "Is the planet in a snowball state?"); + fvFormattedString(&output[OUT_SNOWBALL].cDescr, + "Is the planet in a snowball state?"); output[OUT_SNOWBALL].bNeg = 0; output[OUT_SNOWBALL].iNum = 1; output[OUT_SNOWBALL].iModuleBit = POISE; fnWrite[OUT_SNOWBALL] = &WriteSnowball; fvFormattedString(&output[OUT_TOTICEMASS].cName, "TotIceMass"); - fvFormattedString(&output[OUT_TOTICEMASS].cDescr, "Global total ice mass in ice sheets"); + fvFormattedString(&output[OUT_TOTICEMASS].cDescr, + "Global total ice mass in ice sheets"); fvFormattedString(&output[OUT_TOTICEMASS].cNeg, "kg"); output[OUT_TOTICEMASS].bNeg = 1; output[OUT_TOTICEMASS].iNum = 1; @@ -3953,7 +4397,7 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fvFormattedString(&output[OUT_TOTICEFLOW].cName, "TotIceFlow"); fvFormattedString(&output[OUT_TOTICEFLOW].cDescr, - "Global total ice flow in ice sheets (should = 0)"); + "Global total ice flow in ice sheets (should = 0)"); fvFormattedString(&output[OUT_TOTICEFLOW].cNeg, "kg"); output[OUT_TOTICEFLOW].bNeg = 1; output[OUT_TOTICEFLOW].iNum = 1; @@ -3962,7 +4406,7 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fvFormattedString(&output[OUT_TOTICEBALANCE].cName, "TotIceBalance"); fvFormattedString(&output[OUT_TOTICEBALANCE].cDescr, - "Global total ice balance in ice sheets (this time step)"); + "Global total ice balance in ice sheets (this time step)"); fvFormattedString(&output[OUT_TOTICEBALANCE].cNeg, "kg"); output[OUT_TOTICEBALANCE].bNeg = 1; output[OUT_TOTICEBALANCE].iNum = 1; @@ -3971,7 +4415,7 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fvFormattedString(&output[OUT_FLUXINGLOBAL].cName, "FluxInGlobal"); fvFormattedString(&output[OUT_FLUXINGLOBAL].cDescr, - "Global mean flux in (insol*(1-albedo)) from POISE"); + "Global mean flux in (insol*(1-albedo)) from POISE"); /* Sadly, Russell, we must set the negative option to W/m^2. fvFormattedString(output[OUT_FLUXINGLOBAL].cNeg,"pirate-ninjas/m^2"); output[OUT_FLUXINGLOBAL].dNeg = 1/40.55185; @@ -3984,7 +4428,8 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fnWrite[OUT_FLUXINGLOBAL] = &WriteFluxInGlobal; fvFormattedString(&output[OUT_FLUXOUTGLOBAL].cName, "FluxOutGlobal"); - fvFormattedString(&output[OUT_FLUXOUTGLOBAL].cDescr, "Global mean flux out from POISE"); + fvFormattedString(&output[OUT_FLUXOUTGLOBAL].cDescr, + "Global mean flux out from POISE"); /* Here, too fvFormattedString(output[OUT_FLUXOUTGLOBAL].cNeg,"pirate-ninjas/m^2"); output[OUT_FLUXOUTGLOBAL].dNeg = 1/40.55185; @@ -3997,7 +4442,8 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fnWrite[OUT_FLUXOUTGLOBAL] = &WriteFluxOutGlobal; fvFormattedString(&output[OUT_TEMPLAT].cName, "TempLat"); - fvFormattedString(&output[OUT_TEMPLAT].cDescr, "Surface temperature by latitude."); + fvFormattedString(&output[OUT_TEMPLAT].cDescr, + "Surface temperature by latitude."); fvFormattedString(&output[OUT_TEMPLAT].cNeg, "Celsius"); output[OUT_TEMPLAT].bNeg = 1; // conversion is hardcoded in write function @@ -4009,7 +4455,7 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fvFormattedString(&output[OUT_TEMPMINLAT].cName, "TempMinLat"); fvFormattedString(&output[OUT_TEMPMINLAT].cDescr, - "Minimum surface temperature over a year by latitude."); + "Minimum surface temperature over a year by latitude."); fvFormattedString(&output[OUT_TEMPMINLAT].cNeg, "Celsius"); output[OUT_TEMPMINLAT].bNeg = 1; // conversion is hardcoded in write function @@ -4021,7 +4467,7 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fvFormattedString(&output[OUT_TEMPMAXLAT].cName, "TempMaxLat"); fvFormattedString(&output[OUT_TEMPMAXLAT].cDescr, - "Maximum surface temperature over a year by latitude."); + "Maximum surface temperature over a year by latitude."); fvFormattedString(&output[OUT_TEMPMAXLAT].cNeg, "Celsius"); output[OUT_TEMPMAXLAT].bNeg = 1; // conversion is hardcoded in write function @@ -4033,7 +4479,7 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fvFormattedString(&output[OUT_TEMPMAXLAND].cName, "TempMaxLand"); fvFormattedString(&output[OUT_TEMPMAXLAND].cDescr, - "Maximum surface temperature on land"); + "Maximum surface temperature on land"); fvFormattedString(&output[OUT_TEMPMAXLAND].cNeg, "Celsius"); output[OUT_TEMPMAXLAND].bNeg = 1; // conversion is hardcoded in write function @@ -4045,7 +4491,7 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fvFormattedString(&output[OUT_TEMPMAXWATER].cName, "TempMaxWater"); fvFormattedString(&output[OUT_TEMPMAXWATER].cDescr, - "Maximum surface temperature on water"); + "Maximum surface temperature on water"); fvFormattedString(&output[OUT_TEMPMAXWATER].cNeg, "Celsius"); output[OUT_TEMPMAXWATER].bNeg = 1; // conversion is hardcoded in write function @@ -4057,7 +4503,7 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fvFormattedString(&output[OUT_TEMPLANDLAT].cName, "TempLandLat"); fvFormattedString(&output[OUT_TEMPLANDLAT].cDescr, - "Land surface temperature by latitude."); + "Land surface temperature by latitude."); fvFormattedString(&output[OUT_TEMPLANDLAT].cNeg, "Celsius"); output[OUT_TEMPLANDLAT].bNeg = 1; // conversion is hardcoded in write function @@ -4068,7 +4514,7 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fvFormattedString(&output[OUT_TEMPWATERLAT].cName, "TempWaterLat"); fvFormattedString(&output[OUT_TEMPWATERLAT].cDescr, - "Water surface temperature by latitude."); + "Water surface temperature by latitude."); fvFormattedString(&output[OUT_TEMPWATERLAT].cNeg, "Celsius"); output[OUT_TEMPWATERLAT].bNeg = 1; // conversion is hardcoded in write function @@ -4088,7 +4534,8 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fnWrite[OUT_LATITUDE] = &WriteLatitude; fvFormattedString(&output[OUT_ALBEDOLAT].cName, "AlbedoLat"); - fvFormattedString(&output[OUT_ALBEDOLAT].cDescr, "Surface albedo by latitude."); + fvFormattedString(&output[OUT_ALBEDOLAT].cDescr, + "Surface albedo by latitude."); output[OUT_ALBEDOLAT].bNeg = 0; output[OUT_ALBEDOLAT].iNum = 1; output[OUT_ALBEDOLAT].bGrid = 1; @@ -4096,7 +4543,8 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fnWrite[OUT_ALBEDOLAT] = &WriteAlbedoLat; fvFormattedString(&output[OUT_ALBEDOLANDLAT].cName, "AlbedoLandLat"); - fvFormattedString(&output[OUT_ALBEDOLANDLAT].cDescr, "Land surface albedo by latitude."); + fvFormattedString(&output[OUT_ALBEDOLANDLAT].cDescr, + "Land surface albedo by latitude."); output[OUT_ALBEDOLANDLAT].bNeg = 0; output[OUT_ALBEDOLANDLAT].iNum = 1; output[OUT_ALBEDOLANDLAT].bGrid = 1; @@ -4104,14 +4552,15 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fvFormattedString(&output[OUT_ALBEDOWATERLAT].cName, "AlbedoWaterLat"); fvFormattedString(&output[OUT_ALBEDOWATERLAT].cDescr, - "Water surface albedo by latitude."); + "Water surface albedo by latitude."); output[OUT_ALBEDOWATERLAT].bNeg = 0; output[OUT_ALBEDOWATERLAT].iNum = 1; output[OUT_ALBEDOWATERLAT].bGrid = 1; fnWrite[OUT_ALBEDOWATERLAT] = &WriteAlbedoWaterLat; fvFormattedString(&output[OUT_ANNUALINSOL].cName, "AnnInsol"); - fvFormattedString(&output[OUT_ANNUALINSOL].cDescr, "Annual insolation by latitude."); + fvFormattedString(&output[OUT_ANNUALINSOL].cDescr, + "Annual insolation by latitude."); fvFormattedString(&output[OUT_ANNUALINSOL].cNeg, "W/m^2"); output[OUT_ANNUALINSOL].bNeg = 1; output[OUT_ANNUALINSOL].dNeg = 1 / 40.55185; @@ -4121,7 +4570,8 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fnWrite[OUT_ANNUALINSOL] = &WriteAnnualInsol; fvFormattedString(&output[OUT_PEAKINSOL].cName, "PeakInsol"); - fvFormattedString(&output[OUT_PEAKINSOL].cDescr, "Peak insolation by latitude."); + fvFormattedString(&output[OUT_PEAKINSOL].cDescr, + "Peak insolation by latitude."); fvFormattedString(&output[OUT_PEAKINSOL].cNeg, "W/m^2"); output[OUT_PEAKINSOL].bNeg = 1; output[OUT_PEAKINSOL].dNeg = 1 / 40.55185; @@ -4132,7 +4582,7 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fvFormattedString(&output[OUT_FLUXMERID].cName, "FluxMerid"); fvFormattedString(&output[OUT_FLUXMERID].cDescr, - "Total meridional (northward) flux by latitude"); + "Total meridional (northward) flux by latitude"); fvFormattedString(&output[OUT_FLUXMERID].cNeg, "PW"); output[OUT_FLUXMERID].bNeg = 1; output[OUT_FLUXMERID].dNeg = 1e-15; @@ -4152,7 +4602,8 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fnWrite[OUT_FLUXIN] = &WriteFluxIn; fvFormattedString(&output[OUT_FLUXOUT].cName, "FluxOut"); - fvFormattedString(&output[OUT_FLUXOUT].cDescr, "Outgoing (spaceward) flux by latitude"); + fvFormattedString(&output[OUT_FLUXOUT].cDescr, + "Outgoing (spaceward) flux by latitude"); fvFormattedString(&output[OUT_FLUXOUT].cNeg, "W/m^2"); output[OUT_FLUXOUT].bNeg = 1; output[OUT_FLUXOUT].dNeg = 1; @@ -4162,8 +4613,9 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fnWrite[OUT_FLUXOUT] = &WriteFluxOut; fvFormattedString(&output[OUT_DIVFLUX].cName, "DivFlux"); - fvFormattedString(&output[OUT_DIVFLUX].cDescr, - "Divergence of flux (flow into adjacent cells) by latitude"); + fvFormattedString( + &output[OUT_DIVFLUX].cDescr, + "Divergence of flux (flow into adjacent cells) by latitude"); fvFormattedString(&output[OUT_DIVFLUX].cNeg, "W/m^2"); output[OUT_DIVFLUX].bNeg = 1; output[OUT_DIVFLUX].dNeg = 1; @@ -4173,7 +4625,8 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fnWrite[OUT_DIVFLUX] = &WriteDivFlux; fvFormattedString(&output[OUT_ICEMASS].cName, "IceMass"); - fvFormattedString(&output[OUT_ICEMASS].cDescr, "Mass of ice sheets/area by latitude"); + fvFormattedString(&output[OUT_ICEMASS].cDescr, + "Mass of ice sheets/area by latitude"); fvFormattedString(&output[OUT_ICEMASS].cNeg, "kg/m^2"); output[OUT_ICEMASS].bNeg = 1; output[OUT_ICEMASS].dNeg = 1; @@ -4194,7 +4647,7 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fvFormattedString(&output[OUT_DICEMASSDT].cName, "DIceMassDt"); fvFormattedString(&output[OUT_DICEMASSDT].cDescr, - "derivative of mass of ice sheets/area by latitude"); + "derivative of mass of ice sheets/area by latitude"); fvFormattedString(&output[OUT_DICEMASSDT].cNeg, "kg/m^2/s"); output[OUT_DICEMASSDT].bNeg = 1; output[OUT_DICEMASSDT].dNeg = 1; @@ -4205,7 +4658,7 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fvFormattedString(&output[OUT_ICEACCUM].cName, "IceAccum"); fvFormattedString(&output[OUT_ICEACCUM].cDescr, - "Ice growth per orbit (accumulation only)"); + "Ice growth per orbit (accumulation only)"); fvFormattedString(&output[OUT_ICEACCUM].cNeg, "m/orbit"); output[OUT_ICEACCUM].bNeg = 1; output[OUT_ICEACCUM].dNeg = 1; @@ -4215,7 +4668,8 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fnWrite[OUT_ICEACCUM] = &WriteIceAccum; fvFormattedString(&output[OUT_ICEABLATE].cName, "IceAblate"); - fvFormattedString(&output[OUT_ICEABLATE].cDescr, "ice decay per orbit (ablation only)"); + fvFormattedString(&output[OUT_ICEABLATE].cDescr, + "ice decay per orbit (ablation only)"); fvFormattedString(&output[OUT_ICEABLATE].cNeg, "m/orbit"); output[OUT_ICEABLATE].bNeg = 1; output[OUT_ICEABLATE].dNeg = 1; @@ -4225,7 +4679,8 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fnWrite[OUT_ICEABLATE] = &WriteIceAblate; fvFormattedString(&output[OUT_ICEFLOW].cName, "IceFlow"); - fvFormattedString(&output[OUT_ICEFLOW].cDescr, "flow of ice sheets/area by latitude"); + fvFormattedString(&output[OUT_ICEFLOW].cDescr, + "flow of ice sheets/area by latitude"); fvFormattedString(&output[OUT_ICEFLOW].cNeg, "m/s"); output[OUT_ICEFLOW].bNeg = 1; output[OUT_ICEFLOW].dNeg = 1; @@ -4255,7 +4710,8 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fnWrite[OUT_ENERGYRESL] = &WriteEnergyResL; fvFormattedString(&output[OUT_ENERGYRESW].cName, "EnergyResW"); - fvFormattedString(&output[OUT_ENERGYRESW].cDescr, "Energy residual over water"); + fvFormattedString(&output[OUT_ENERGYRESW].cDescr, + "Energy residual over water"); fvFormattedString(&output[OUT_ENERGYRESW].cNeg, "W/m^2"); output[OUT_ENERGYRESW].bNeg = 1; output[OUT_ENERGYRESW].dNeg = 1; @@ -4266,14 +4722,15 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fvFormattedString(&output[OUT_SKIPSEAS].cName, "SkipSeas"); fvFormattedString(&output[OUT_SKIPSEAS].cDescr, - "Is Seasonal model skipped due to RGH or snowball?"); + "Is Seasonal model skipped due to RGH or snowball?"); output[OUT_SKIPSEAS].bNeg = 0; output[OUT_SKIPSEAS].iNum = 1; output[OUT_SKIPSEAS].iModuleBit = POISE; fnWrite[OUT_SKIPSEAS] = &WriteSkipSeas; fvFormattedString(&output[OUT_PLANCKBAVG].cName, "PlanckBAvg"); - fvFormattedString(&output[OUT_PLANCKBAVG].cDescr, "Annually averaged Planck B coeff"); + fvFormattedString(&output[OUT_PLANCKBAVG].cDescr, + "Annually averaged Planck B coeff"); fvFormattedString(&output[OUT_PLANCKBAVG].cNeg, "W/m^2/C"); output[OUT_PLANCKBAVG].bNeg = 0; output[OUT_PLANCKBAVG].iNum = 1; @@ -4282,7 +4739,8 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fnWrite[OUT_PLANCKBAVG] = &WritePlanckBAvg; fvFormattedString(&output[OUT_AREAICECOV].cName, "AreaIceCov"); - fvFormattedString(&output[OUT_AREAICECOV].cDescr, "Fractional area ice covered"); + fvFormattedString(&output[OUT_AREAICECOV].cDescr, + "Fractional area ice covered"); fvFormattedString(&output[OUT_AREAICECOV].cNeg, " "); output[OUT_AREAICECOV].bNeg = 1; output[OUT_AREAICECOV].iNum = 1; @@ -4291,7 +4749,7 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fvFormattedString(&output[OUT_NORTHICECAPLAND].cName, "IceCapNorthLand"); fvFormattedString(&output[OUT_NORTHICECAPLAND].cDescr, - "Does the planet have a northern polar ice cap on land?"); + "Does the planet have a northern polar ice cap on land?"); output[OUT_NORTHICECAPLAND].bNeg = 0; output[OUT_NORTHICECAPLAND].iNum = 1; output[OUT_NORTHICECAPLAND].iModuleBit = POISE; @@ -4299,7 +4757,7 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fvFormattedString(&output[OUT_NORTHICECAPSEA].cName, "IceCapNorthSea"); fvFormattedString(&output[OUT_NORTHICECAPSEA].cDescr, - "Does the planet have a northern polar sea ice cap"); + "Does the planet have a northern polar sea ice cap"); output[OUT_NORTHICECAPSEA].bNeg = 0; output[OUT_NORTHICECAPSEA].iNum = 1; output[OUT_NORTHICECAPSEA].iModuleBit = POISE; @@ -4307,7 +4765,7 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fvFormattedString(&output[OUT_SOUTHICECAPLAND].cName, "IceCapSouthLand"); fvFormattedString(&output[OUT_SOUTHICECAPLAND].cDescr, - "Does the planet have a southern polar ice cap on land?"); + "Does the planet have a southern polar ice cap on land?"); output[OUT_SOUTHICECAPLAND].bNeg = 0; output[OUT_SOUTHICECAPLAND].iNum = 1; output[OUT_SOUTHICECAPLAND].iModuleBit = POISE; @@ -4315,7 +4773,7 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fvFormattedString(&output[OUT_SOUTHICECAPSEA].cName, "IceCapSouthSea"); fvFormattedString(&output[OUT_SOUTHICECAPSEA].cDescr, - "Does the planet have a southern polar sea ice cap?"); + "Does the planet have a southern polar sea ice cap?"); output[OUT_SOUTHICECAPSEA].bNeg = 0; output[OUT_SOUTHICECAPSEA].iNum = 1; output[OUT_SOUTHICECAPSEA].iModuleBit = POISE; @@ -4323,7 +4781,7 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fvFormattedString(&output[OUT_ICEBELTLAND].cName, "IceBeltLand"); fvFormattedString(&output[OUT_ICEBELTLAND].cDescr, - "Does the planet have a land ice belt?"); + "Does the planet have a land ice belt?"); output[OUT_ICEBELTLAND].bNeg = 0; output[OUT_ICEBELTLAND].iNum = 1; output[OUT_ICEBELTLAND].iModuleBit = POISE; @@ -4331,21 +4789,23 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fvFormattedString(&output[OUT_ICEBELTSEA].cName, "IceBeltSea"); fvFormattedString(&output[OUT_ICEBELTSEA].cDescr, - "Does the planet have a sea ice belt?"); + "Does the planet have a sea ice belt?"); output[OUT_ICEBELTSEA].bNeg = 0; output[OUT_ICEBELTSEA].iNum = 1; output[OUT_ICEBELTSEA].iModuleBit = POISE; fnWrite[OUT_ICEBELTSEA] = &WriteIceBeltSea; fvFormattedString(&output[OUT_SNOWBALLLAND].cName, "SnowballLand"); - fvFormattedString(&output[OUT_SNOWBALLLAND].cDescr, "Is all land covered in ice?"); + fvFormattedString(&output[OUT_SNOWBALLLAND].cDescr, + "Is all land covered in ice?"); output[OUT_SNOWBALLLAND].bNeg = 0; output[OUT_SNOWBALLLAND].iNum = 1; output[OUT_SNOWBALLLAND].iModuleBit = POISE; fnWrite[OUT_SNOWBALLLAND] = &WriteSnowballLand; fvFormattedString(&output[OUT_SNOWBALLSEA].cName, "SnowballSea"); - fvFormattedString(&output[OUT_SNOWBALLSEA].cDescr, "Is all sea covered in ice?"); + fvFormattedString(&output[OUT_SNOWBALLSEA].cDescr, + "Is all sea covered in ice?"); output[OUT_SNOWBALLSEA].bNeg = 0; output[OUT_SNOWBALLSEA].iNum = 1; output[OUT_SNOWBALLSEA].iModuleBit = POISE; @@ -4353,59 +4813,66 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { fvFormattedString(&output[OUT_ICEFREE].cName, "IceFree"); fvFormattedString(&output[OUT_ICEFREE].cDescr, - "Is the planet free of sea and land ice?"); + "Is the planet free of sea and land ice?"); output[OUT_ICEFREE].bNeg = 0; output[OUT_ICEFREE].iNum = 1; output[OUT_ICEFREE].iModuleBit = POISE; fnWrite[OUT_ICEFREE] = &WriteIceFree; - fvFormattedString(&output[OUT_NORTHICECAPLATLAND].cName, "IceCapNorthLatLand"); + fvFormattedString(&output[OUT_NORTHICECAPLATLAND].cName, + "IceCapNorthLatLand"); fvFormattedString(&output[OUT_NORTHICECAPLATLAND].cDescr, - "Southernmost extent of northern land ice cap."); + "Southernmost extent of northern land ice cap."); output[OUT_NORTHICECAPLATLAND].bNeg = 0; output[OUT_NORTHICECAPLATLAND].iNum = 1; output[OUT_NORTHICECAPLATLAND].iModuleBit = POISE; fnWrite[OUT_NORTHICECAPLATLAND] = &WriteIceCapNorthLatLand; - fvFormattedString(&output[OUT_NORTHICECAPLATLAND].cLongDescr, - "If a northern land ice cap is present, return the latitude of its " - "southern edge. If not present, return +100 degrees."); + fvFormattedString( + &output[OUT_NORTHICECAPLATLAND].cLongDescr, + "If a northern land ice cap is present, return the latitude of its " + "southern edge. If not present, return +100 degrees."); fvFormattedString(&output[OUT_NORTHICECAPLATSEA].cName, "IceCapNorthLatSea"); fvFormattedString(&output[OUT_NORTHICECAPLATSEA].cDescr, - "Southernmost extent of northern sea ice cap."); + "Southernmost extent of northern sea ice cap."); output[OUT_NORTHICECAPLATSEA].bNeg = 0; output[OUT_NORTHICECAPLATSEA].iNum = 1; output[OUT_NORTHICECAPLATSEA].iModuleBit = POISE; fnWrite[OUT_NORTHICECAPLATSEA] = &WriteIceCapNorthLatSea; - fvFormattedString(&output[OUT_NORTHICECAPLATSEA].cLongDescr, - "If a northern sea ice cap is present, return the latitude of its " - "southern edge. If not present, return +100 degrees."); + fvFormattedString( + &output[OUT_NORTHICECAPLATSEA].cLongDescr, + "If a northern sea ice cap is present, return the latitude of its " + "southern edge. If not present, return +100 degrees."); - fvFormattedString(&output[OUT_SOUTHICECAPLATLAND].cName, "IceCapSouthLatLand"); + fvFormattedString(&output[OUT_SOUTHICECAPLATLAND].cName, + "IceCapSouthLatLand"); fvFormattedString(&output[OUT_SOUTHICECAPLATLAND].cDescr, - "Northernmost extent of southern land ice cap."); + "Northernmost extent of southern land ice cap."); output[OUT_SOUTHICECAPLATLAND].bNeg = 0; output[OUT_SOUTHICECAPLATLAND].iNum = 1; output[OUT_SOUTHICECAPLATLAND].iModuleBit = POISE; fnWrite[OUT_SOUTHICECAPLATLAND] = &WriteIceCapSouthLatLand; - fvFormattedString(&output[OUT_SOUTHICECAPLATLAND].cLongDescr, - "If a southern land ice cap is present, return the latitude of its " - "northern edge. If not present, return -100 degrees."); + fvFormattedString( + &output[OUT_SOUTHICECAPLATLAND].cLongDescr, + "If a southern land ice cap is present, return the latitude of its " + "northern edge. If not present, return -100 degrees."); fvFormattedString(&output[OUT_SOUTHICECAPLATSEA].cName, "IceCapSouthLatSea"); fvFormattedString(&output[OUT_SOUTHICECAPLATSEA].cDescr, - "Northernmost extent of southern sea ice cap."); + "Northernmost extent of southern sea ice cap."); output[OUT_SOUTHICECAPLATSEA].bNeg = 0; output[OUT_SOUTHICECAPLATSEA].iNum = 1; output[OUT_SOUTHICECAPLATSEA].iModuleBit = POISE; fnWrite[OUT_SOUTHICECAPLATSEA] = &WriteIceCapSouthLatLand; - fvFormattedString(&output[OUT_SOUTHICECAPLATSEA].cLongDescr, - "If a southern sea ice cap is present, return the latitude of its " - "northern edge. If not present, return -100 degrees."); + fvFormattedString( + &output[OUT_SOUTHICECAPLATSEA].cLongDescr, + "If a southern sea ice cap is present, return the latitude of its " + "northern edge. If not present, return -100 degrees."); - fvFormattedString(&output[OUT_NORTHICEBELTLATLAND].cName, "IceBeltNorthLatLand"); + fvFormattedString(&output[OUT_NORTHICEBELTLATLAND].cName, + "IceBeltNorthLatLand"); fvFormattedString(&output[OUT_NORTHICEBELTLATLAND].cDescr, - "Northernmost extent of land ice belt."); + "Northernmost extent of land ice belt."); output[OUT_NORTHICEBELTLATLAND].bNeg = 0; output[OUT_NORTHICEBELTLATLAND].iNum = 1; output[OUT_NORTHICEBELTLATLAND].iModuleBit = POISE; @@ -4417,9 +4884,10 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { "If not present, return 0. Note that some ice belts may in fact have a " "northern edge at the equator."); - fvFormattedString(&output[OUT_NORTHICEBELTLATSEA].cName, "IceBeltNorthLatSea"); + fvFormattedString(&output[OUT_NORTHICEBELTLATSEA].cName, + "IceBeltNorthLatSea"); fvFormattedString(&output[OUT_NORTHICEBELTLATSEA].cDescr, - "Northernmost extent of sea ice belt."); + "Northernmost extent of sea ice belt."); output[OUT_NORTHICEBELTLATSEA].bNeg = 0; output[OUT_NORTHICEBELTLATSEA].iNum = 1; output[OUT_NORTHICEBELTLATSEA].iModuleBit = POISE; @@ -4431,9 +4899,10 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { "If not present, return 0. Note that some ice belts may in fact have a " "northern edge at the equator."); - fvFormattedString(&output[OUT_SOUTHICEBELTLATLAND].cName, "IceBeltSouthLatLand"); + fvFormattedString(&output[OUT_SOUTHICEBELTLATLAND].cName, + "IceBeltSouthLatLand"); fvFormattedString(&output[OUT_SOUTHICEBELTLATLAND].cDescr, - "Southernmost extent of land ice belt."); + "Southernmost extent of land ice belt."); output[OUT_SOUTHICEBELTLATLAND].bNeg = 0; output[OUT_SOUTHICEBELTLATLAND].iNum = 1; output[OUT_SOUTHICEBELTLATLAND].iModuleBit = POISE; @@ -4445,9 +4914,10 @@ void InitializeOutputPoise(OUTPUT *output, fnWriteOutput fnWrite[]) { "If not present, return 0. Note that some ice belts may in fact have a " "southern edge at the equator."); - fvFormattedString(&output[OUT_SOUTHICEBELTLATSEA].cName, "IceBeltSouthLatSea"); + fvFormattedString(&output[OUT_SOUTHICEBELTLATSEA].cName, + "IceBeltSouthLatSea"); fvFormattedString(&output[OUT_SOUTHICEBELTLATSEA].cDescr, - "Southernmost extent of sea ice belt."); + "Southernmost extent of sea ice belt."); output[OUT_SOUTHICEBELTLATSEA].bNeg = 0; output[OUT_SOUTHICEBELTLATSEA].iNum = 1; output[OUT_SOUTHICEBELTLATSEA].iModuleBit = POISE; @@ -4458,6 +4928,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 **************/ @@ -7771,3 +8249,14 @@ void PoiseIceSheets(BODY *body, EVOLVE *evolve, int iBody) { } } } + +double fdLandFracGlobal(BODY *body, int iBody) { + int iLat; + double dLandFrac, dLandArea = 0; + + for (iLat = 0; iLat < body[iBody].iNumLats; iLat++) { + dLandArea += (body[iBody].daLandFrac[iLat] * body[iBody].daLatsArea[iLat]); + } + dLandFrac = dLandArea / (4 * PI * body[iBody].dRadius * body[iBody].dRadius); + return dLandFrac; +} diff --git a/src/poise.h b/src/poise.h index 2e028c21e..333a9487d 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 GEOGRAPHYUNIFORM 0 +#define GEOGRAPHYMODERN 1 +#define GEOGRAPHYRANDOM 2 +#define GEOGRAPHYPOLAR 3 +#define GEOGRAPHYEQUATORIAL 4 // Constants for the ice model #define LFICE 3.34e5 // ??? @@ -39,6 +42,8 @@ 1.733e3 // coeff of ice deformability at T>=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 @@ -78,14 +83,16 @@ #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_HEATCAPLAND 1942 -#define OPT_HEATCAPWATER 1943 -#define OPT_FRZTSEAICE 1944 -//#define OPT_LATENTHEAT 1945 +#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 #define OPT_NULANDWATER 1948 @@ -111,6 +118,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 @@ -178,6 +186,7 @@ #define OUT_NORTHICEBELTLATSEA 1973 #define OUT_SOUTHICEBELTLATLAND 1974 #define OUT_SOUTHICEBELTLATSEA 1975 +#define OUT_LANDFRACGLOBAL 1976 /* @cond DOXYGEN_OVERRIDE */ @@ -190,7 +199,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 *, @@ -209,39 +218,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 */ @@ -302,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 ff52fb16a..8e943312a 100644 --- a/src/vplanet.h +++ b/src/vplanet.h @@ -632,17 +632,19 @@ struct BODY { needed) */ double dAlbedoLand; /**< Sets base albedo of land (sea model) */ double dAlbedoWater; /**< Sets base albedo of water (sea model) */ - 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 dLatLandWater; /**< 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 */ double dFluxInGlobal; /**< Global mean of incoming flux */ double dFluxInGlobalTmp; /**< Copy of global mean incoming flux */ double dFluxOutGlobal; /**< Global mean of outgoing flux */ @@ -720,7 +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); 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) */ @@ -800,7 +806,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) */ 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