Skip to content

Libration of SSP and additional spin-orbit resonances #5

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
hersfeldtn wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Libration of SSP and additional spin-orbit resonances #5

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
bd4dc54
Add allowlibrate and soldaysperorb parameters to __init__.py
hersfeldtn Jul 22, 2023
7a6def9
alter radmod.f90 to allow ssp libration and non-1:1 resonance
hersfeldtn Jul 22, 2023
ddd3ed4
eccentricity libration doesn't actually require ngenkeplerian
hersfeldtn Jul 22, 2023
7a8237b
Update documentation to match last commit
hersfeldtn Jul 22, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 32 additions & 2 deletions exoplasim/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1898,7 +1898,7 @@ def configure(self,noutput=True,flux=1367.0,startemp=None,starradius=1.0,starspe
pHe=None,pN2=None,pO2=None,pCO2=None,pCH4=None,pAr=None,pNe=None,
pKr=None,pH2O=None,gascon=None,pressure=None,pressurebroaden=True,
vtype=0,rotationperiod=1.0,synchronous=False,substellarlon=180.0,
keplerian=False,meananomaly0=None,
keplerian=False,meananomaly0=None,allowlibrate=True,soldaysperorb=0.0,
year=None,glaciers={"toggle":False,"mindepth":2.0,"initialh":-1.0},
restartfile=None,gravity=10.9,radius=1.12,eccentricity=0.0,
obliquity=0.0,lonvernaleq=None,fixedorbit=True,orography=None,
Expand Down Expand Up @@ -2037,6 +2037,12 @@ def configure(self,noutput=True,flux=1367.0,startemp=None,starradius=1.0,starspe
The rate of drift of the substellar point in degrees per minute. May be positive or negative.
substellarlon : float, optional
The longitude of the substellar point, if synchronous==True. Default 180°
allowlibrate : bool, optional
True/False. If True and synchronous==True, allows for libration of the substellar point due to obliquity and eccentricity. True by default.
soldaysperorb : float, optional
If synchronous==True, solar days experienced per orbit. 0.0 is typical tidal-locked behavior,
but other values can be used for other spin-orbit resonances, e.g. 0.5 for 3:2 resonance
or 1.0 for 2:1. In these cases, substellarlon marks the starting substellar longitude (before accounting for libration).
pressurebroaden : bool, optional
True/False. If False, pressure-broadening of absorbers no longer depends
on surface pressure. Default is True
Expand Down Expand Up @@ -2579,6 +2585,12 @@ def configure(self,noutput=True,flux=1367.0,startemp=None,starradius=1.0,starspe
if meananomaly0 is not None:
self._edit_namelist("planet_namelist","MEANANOM0",str(meananomaly0))
self.meananomaly0 = meananomaly0

self._edit_namelist("radmod_namelist","ALLOWLIBRATE",str(allowlibrate*1))
self.allowlibrate = allowlibrate

self._edit_namelist("radmod_namelist","DAYSPERORB",str(soldaysperorb))
self.soldaysperorb = soldaysperorb

if type(orography)!=type(None):
self._edit_namelist("landmod_namelist","OROSCALE",str(orography))
Expand Down Expand Up @@ -3056,6 +3068,15 @@ def loadconfig(self,configfile):
aerobulk = 1
aerorad = True
aerofile = None

try:
allowlibrate = bool(cfg[93])
except:
allowlibrate = True
try:
soldaysperorb = float(cfg[94])
except:
soldaysperorb = 0.0

self.configure(noutput=noutput,flux=flux,startemp=startemp,starspec=starspec,starradius=starradius,
gascon=gascon,pressure=pressure,pressurebroaden=pressurebroaden,
Expand Down Expand Up @@ -3085,7 +3106,8 @@ def loadconfig(self,configfile):
nstorms=nstorms,stormcapture=stormcapture,topomap=topomap,tlcontrast=tlcontrast,
otherargs=otherargs,glaciers=glaciers,threshold=threshold,keplerian=keplerian,
meananomaly0=meananomaly0,apart=apart,rhop=rhop,fcoeff=fcoeff,aerobulk=aerobulk,
aerorad=aerorad,aerosol=aerosol,asource=asource,aerofile=aerofile)
aerorad=aerorad,aerosol=aerosol,asource=asource,aerofile=aerofile,
allowlibrate=allowlibrate,soldaysperorb=soldaysperorb)

def modify(self,**kwargs):
"""Modify any already-configured parameters. All parameters accepted by :py:func:`configure() <exoplasim.Model.configure>` can be passed as arguments.
Expand Down Expand Up @@ -3293,6 +3315,12 @@ def modify(self,**kwargs):
self.meananomaly0=value
if self.meanomaly0 is not None:
self._edit_namelist("planet_namelist","MEANANOMALY0",str(self.meananomaly0))
if key=="allowlibrate":
self.allowlibrate=value
self._edit_namelist("radmod_namelist","ALLOWLIBRATE",str(allowlibrate*1))
if key=="soldaysperorb":
self.soldaysperorb=value
self._edit_namelist("radmod_namelist","DAYSPERORB",str(soldaysperorb))
if key=="tlcontrast":
self.tlcontrast=value
self._edit_namelist("plasim_namelist","DTTL",str(self.tlcontrast))
Expand Down Expand Up @@ -4032,6 +4060,8 @@ def exportcfg(self,filename=None):
cfg.append(str(self.aerobulk))
cfg.append(str(self.aerorad*1))
cfg.append(str(self.aerofile))
cfg.append(str(self.allowlibrate*1))
cfg.append(str(self.soldaysperorb))

print("Writing configuration....\n"+"\n".join(cfg))
print("Writing to %s...."%filename)
Expand Down
38 changes: 30 additions & 8 deletions exoplasim/plasim/src/radmod.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,8 +70,12 @@ module radmod

integer :: npbroaden = 1 ! Should pressure broadening depend on surface pressure (1/0)
integer :: nfixed = 0 ! Switch for fixed zenith angle (0/1=no/yes)
integer :: allowlibrate = 1 ! Switch for libration of substellar point
real :: slowdown = 1.0 ! Factor by which to change diurnal insolation cycle
real :: desync = 0.0 ! Degrees per minute by which substellar point drifts (+/-)
real :: daysperorb = 0.0 ! Solar days per orbit if synchronous
real :: syncdays = 0.0 ! Solar days completed if synchronous
real :: lonoffset = 0.0 ! Current offset of substellar point


real :: minwavel = 316.036116751 ! Minimum wavelength to use when computing spectra [nm]
Expand Down Expand Up @@ -604,7 +608,7 @@ subroutine radini
!** 0) define namelist
!
namelist/radmod_nl/ndcycle,ncstsol,solclat,solcdec,no3,co2 &
& ,iyrbp,nswr,nlwr,nfixed,slowdown,nradice,npbroaden,desync &
& ,iyrbp,nswr,nlwr,nfixed,slowdown,nradice,npbroaden,desync,daysperorb,allowlibrate &
& ,a0o3,a1o3,aco3,bo3,co3,toffo3,o3scale,newrsc,necham,necham6 &
& ,nsol,nclouds,nswrcl,nrscat,rcl1,rcl2,acl2,clgray,tpofmt &
& ,acllwr,tswr1,tswr2,tswr3,th2oc,dawn,starbbtemp,nstartemp &
Expand Down Expand Up @@ -738,7 +742,10 @@ subroutine radini
minwavel = minwavel*1.0e-9

oldfixedlon = fixedlon
if (nrestart > 0.) call get_restart_real('fixedlon',fixedlon)
if (nrestart > 0.) then
call get_restart_real('fixedlon',fixedlon)
if (daysperorb /= 0.0) call get_restart_real('syncdays',syncdays)
endif
if ((fixedlon .ne. oldfixedlon) .and. (desync .eq. 0.0)) fixedlon=oldfixedlon

if ((necham.eq.1).and.(necham6.eq.1)) necham=0 !necham6 overrides necham
Expand All @@ -752,6 +759,8 @@ subroutine radini
call mpbci(nfixed)
call mpbcr(fixedlon)
call mpbcr(desync)
call mpbcr(daysperorb)
call mpbcr(syncdays)
call mpbcr(slowdown)
call mpbcr(a0o3)
call mpbcr(a1o3)
Expand Down Expand Up @@ -786,6 +795,7 @@ subroutine radini
call mpbci(necham6)
call mpbci(nradice)
call mpbci(npbroaden)
call mpbci(allowlibrate)

call mpbcr(starbbtemp)
call mpbci(nstartemp)
Expand Down Expand Up @@ -1313,7 +1323,10 @@ subroutine radstop
!
! no PUMA variables are used
!
if (mypid==NROOT) call put_restart_real('fixedlon',fixedlon)
if (mypid==NROOT) then
call put_restart_real('fixedlon',fixedlon)
call put_restart_real('syncdays',syncdays)
endif
call mpputgp('zsolars',zsolars,2,1)


Expand Down Expand Up @@ -1395,19 +1408,28 @@ end function ndayofyear
zmins = ihou * 60 + imin

if (nfixed==1) then
if (daysperorb==0.0) then
syncdays = 0.0
else
if (abs(zcday * daysperorb) < abs(mod(syncdays,daysperorb))) syncdays = syncdays + daysperorb/n_steps_per_year !ensures partial rotations carry across years
syncdays = syncdays - mod(syncdays,daysperorb) + zcday * daysperorb
endif
call mpbcr(syncdays)
lonoffset = mod(syncdays,1.0) * (-360.0)
if (allowlibrate > 0) lonoffset = lonoffset + (orbnu*180.0/PI - (zcday*360.0 + meananom0))
if (mypid==NROOT) fixedlon = fixedlon + desync*mpstep
call mpbcr(fixedlon)
zrtim = TWOPI
zmins = 1.0 - (fixedlon/360.) !Think about how to fix this: there's a dep
zdecl = obliqr !on rotspd. Maybe zrtim = TWOPI/1440.0?
zmins = 1.0 - ((fixedlon + lonoffset)/360.) !Think about how to fix this: there's a dep
if (allowlibrate == 0) zdecl = obliqr !on rotspd. Maybe zrtim = TWOPI/1440.0?
endif
jhor = 0
if (ncstsol==0) then
do jlat = 1 , NLPP
do jlon = 0 , NLON-1
jhor = jhor + 1
zhangle = zmins * zrtim + jlon * zrlon - PI
if (ngenkeplerian==1) zhangle = zhangle - rasc
if (ngenkeplerian==1 .AND. nfixed==0) zhangle = zhangle - rasc
if (zhangle < -PI) zhangle = zhangle + TWOPI
if (zhangle > PI) zhangle = zhangle - TWOPI

Expand Down Expand Up @@ -2468,7 +2490,7 @@ end subroutine lwr
! will work with reasonable accuracy for any bound orbit (eccen<1.0).

subroutine gen_orb_decl(yearfraction, eccen, obliqr, mvelpp, trueanomaly, lamb, rasc, zdecl, eccf)
use radmod, only : TWOPI, meananom0r, nfixed
use radmod, only : TWOPI, meananom0r, nfixed, allowlibrate
!Inputs
real :: eccen ! Eccentricity
real :: yearfraction ! Elapsed fraction of the year
Expand All @@ -2488,7 +2510,7 @@ subroutine gen_orb_decl(yearfraction, eccen, obliqr, mvelpp, trueanomaly, lamb,
real :: lamb ! True anomaly - longitude of vernal equinox
real :: rasc ! Right ascension

if (nfixed > 0) then
if (nfixed > 0 .AND. allowlibrate==0) then
trueanomaly = 0.
eccf = 1.
else
Expand Down