From bd4dc5430f8f80d88bd4b5e9135d215cce58a4a4 Mon Sep 17 00:00:00 2001 From: hersfeldtn <85223968+hersfeldtn@users.noreply.github.com> Date: Sat, 22 Jul 2023 19:57:41 +0100 Subject: [PATCH 1/4] Add allowlibrate and soldaysperorb parameters to __init__.py --- exoplasim/__init__.py | 35 +++++++++++++++++++++++++++++++++-- 1 file changed, 33 insertions(+), 2 deletions(-) diff --git a/exoplasim/__init__.py b/exoplasim/__init__.py index 58350f5..ad339fb 100755 --- a/exoplasim/__init__.py +++ b/exoplasim/__init__.py @@ -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, @@ -2037,6 +2037,13 @@ 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, allows for libration of the substellar point due to obliquity and, + if keplerian==True, 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 @@ -2579,6 +2586,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)) @@ -3056,6 +3069,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, @@ -3085,7 +3107,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() ` can be passed as arguments. @@ -3293,6 +3316,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)) @@ -4032,6 +4061,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) From 7a6def972aa058f2d1b515b067b37db272d38687 Mon Sep 17 00:00:00 2001 From: hersfeldtn <85223968+hersfeldtn@users.noreply.github.com> Date: Sat, 22 Jul 2023 20:09:28 +0100 Subject: [PATCH 2/4] alter radmod.f90 to allow ssp libration and non-1:1 resonance --- exoplasim/plasim/src/radmod.f90 | 38 ++++++++++++++++++++++++++------- 1 file changed, 30 insertions(+), 8 deletions(-) diff --git a/exoplasim/plasim/src/radmod.f90 b/exoplasim/plasim/src/radmod.f90 index c52440e..a293c48 100644 --- a/exoplasim/plasim/src/radmod.f90 +++ b/exoplasim/plasim/src/radmod.f90 @@ -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] @@ -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 & @@ -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 @@ -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) @@ -786,6 +795,7 @@ subroutine radini call mpbci(necham6) call mpbci(nradice) call mpbci(npbroaden) + call mpbci(allowlibrate) call mpbcr(starbbtemp) call mpbci(nstartemp) @@ -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) @@ -1395,11 +1408,20 @@ 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 .AND. ngenkeplerian > 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 @@ -1407,7 +1429,7 @@ end function ndayofyear 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 @@ -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 @@ -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 From ddd3ed426f687a7d804b0b3a692a353ca8fbb9da Mon Sep 17 00:00:00 2001 From: hersfeldtn <85223968+hersfeldtn@users.noreply.github.com> Date: Sat, 22 Jul 2023 20:22:11 +0100 Subject: [PATCH 3/4] eccentricity libration doesn't actually require ngenkeplerian --- exoplasim/plasim/src/radmod.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exoplasim/plasim/src/radmod.f90 b/exoplasim/plasim/src/radmod.f90 index a293c48..e2a7963 100644 --- a/exoplasim/plasim/src/radmod.f90 +++ b/exoplasim/plasim/src/radmod.f90 @@ -1416,7 +1416,7 @@ end function ndayofyear endif call mpbcr(syncdays) lonoffset = mod(syncdays,1.0) * (-360.0) - if (allowlibrate > 0 .AND. ngenkeplerian > 0) lonoffset = lonoffset + (orbnu*180.0/PI - (zcday*360.0 + meananom0)) + 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 From 7a8237bfda35ae2489f3c1d08dd1d63919c5fbc3 Mon Sep 17 00:00:00 2001 From: hersfeldtn <85223968+hersfeldtn@users.noreply.github.com> Date: Sat, 22 Jul 2023 20:24:36 +0100 Subject: [PATCH 4/4] Update documentation to match last commit --- exoplasim/__init__.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/exoplasim/__init__.py b/exoplasim/__init__.py index ad339fb..0fefa29 100755 --- a/exoplasim/__init__.py +++ b/exoplasim/__init__.py @@ -2038,8 +2038,7 @@ def configure(self,noutput=True,flux=1367.0,startemp=None,starradius=1.0,starspe substellarlon : float, optional The longitude of the substellar point, if synchronous==True. Default 180° allowlibrate : bool, optional - True/False. If True, allows for libration of the substellar point due to obliquity and, - if keplerian==True, eccentricity. True by default. + 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