File:Temperature of gliese 12b as locked desert planet 1.png

3. tidally locked desert planet temperature simulation python3 climt code
6. 1. Important note:
8. uses Ubuntu Anaconda Python 3.10.11 and sometimes gfs_dynamical_core
9. 1. 27.05.2023 v 0000.0005
12. 1. if you use gfs_dynamical_core
13. 1. must create Python 3.10 enviroment in conda, ant then: pip install gfs-dynamical-core
14. maybe not function in other python versions!
    1. gfs_dynamical_core
15. pip install gfs_dynamical core
    1. of from source download sit source tree
16. https://github.com/Ai33L/gfs_dynamical_core
    1. and mote to its directory "gfs_dynamical_core-develop"
17. and then in python 3.9 ## pip install .
    1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 3

import matplotlib as mpl import math import numpy as np import scipy.ndimage

from datetime import timedelta from datetime import datetime

from netCDF4 import Dataset

1. import sympl
2. import gfs_dynamical_core

import climt

from sympl import (

   DataArray, PlotFunctionMonitor,NetCDFMonitor,
   AdamsBashforth

)

from climt import (

   EmanuelConvection, RRTMGShortwave, RRTMGLongwave, SlabSurface,
   DryConvectiveAdjustment, SimplePhysics, BucketHydrology, get_default_state

)

1. 1. basic params

steppi=4 # timestep hours

1. dimx=16 ## dim of grid
2. dimy=8

dimz=8

dimx=16 ## dim of grid dimy=8 dimz=0

1. pix=8 ## temperature et al specimen point insige grid
2. piy=4
3. piz=3

1. S1=1.0

S1=1.63 ## solar constant coeff

co2coeff=1.0 ## co2 x330

1. albedo=0.07 ## basalt surface like Moon
2. albedo=0.30 ## like Earth

albedo=0.3

qq=0

qqlimit=int(24/steppi)

hourticks1=np.arange(0,24, steppi)

1. print (hourticks1)

1. quit(-1)

caption1="Desert S0 = "+str(S1)

tempcache=np.zeros(qqlimit)+15

def save_netcdf_x(name1, blok1): sheip1=np.shape(blok1) dimx=sheip1[1] dimy=sheip1[0] ds = Dataset(name1, 'w', format='NETCDF4') lat = ds.createDimension('lat', dimy) lon = ds.createDimension('lon', dimx) lats = ds.createVariable('lat', 'f4', ('lat',)) lons = ds.createVariable('lon', 'f4', ('lon',)) Tempe = ds.createVariable('Temperature', 'f4', ('lat', 'lon',)) Tempe.units = 'deg C' lats[:] = np.arange(90, -90, -180/(dimy+0)) lons[:] = np.arange(-180.0, 180, 360/(dimx+0)) Tempe[:, :] = blok1 ds.close()

def plot_function(fig, state):

   global qq, qqlimit
   global tempcache
   global pix, piy
   global hourticks1
   vmin1=-100
   vmax1=100
   Tsurf=np.array(state['surface_temperature']-273.15)
   lons1=np.mean(state['longitude'], axis=0)
   lats1=np.mean(state['latitude'], axis=1)
   zonetemps=np.array(np.mean( Tsurf, axis=1))
   lats2=np.array(np.mean( state['latitude'], axis=1))
   print(zonetemps)
   print(lats2)
   lats3=lats2*math.pi/10
   coslats=np.cos(lats2)
   weightemps=zonetemps*np.abs(coslats)*np.abs(coslats)
   meantemp=np.mean(weightemps)
   print(" T mean ",meantemp)
   ax = fig.add_subplot(2, 2, 1)

   ax.contourf(state['longitude'], state['latitude'],
               Tsurf, cmap=mpl.colormaps['turbo'], origin='lower',levels=10,extent=[0,360,-90,90])
   #ax.imshow(Tsurf, origin='lower', vmin=vmin1, vmax=vmax1, cmap=mpl.colormaps['turbo'],extent=[0,360,-90,90])
   cs1=ax.contour(state['longitude'], state['latitude'],
               state['surface_temperature']-273.15, colors="k",levels=10, origin='lower', vmin=vmin1, vmax=vmax1, extent=[0,360,-90,90])
   
   ax.clabel(cs1, cs1.levels, inline=True, fmt="%.1f", fontsize=16)
   ax.set_xlabel('Longitude')
   ax.set_ylabel('Latitude')

   fig.suptitle('Surface temperature '+str(round(meantemp,2))+' at time: '+str(state['time']))
   fig.canvas.draw()
   fig.canvas.flush_events()
   ax2 = fig.add_subplot(2, 2, 2)
   #fig.suptitle('Tsurf lon profile: '+str(state['time']))
   ax2.plot(lons1, np.mean(Tsurf, axis=0))
   fig.canvas.draw()
   fig.canvas.flush_events()  
   ax3 = fig.add_subplot(2, 2, 3)
   #fig.suptitle('Tsurf lat profile: '+str(state['time']))
   ax3.plot(lats1,np.mean( Tsurf, axis=1))
   ax3.plot(lats1,zonetemps)
   fig.canvas.draw()
   fig.canvas.flush_events()
   qq=qq+1
   if(qq>=qqlimit): qq=0
   tempcache[qq]=Tsurf[piy, pix]
   ax4 = fig.add_subplot(2, 2, 4)
   #fig.suptitle('Tsurf time profile: '+str(qq) )
   ax4.plot(hourticks1, tempcache)
   fig.canvas.draw()
   fig.canvas.flush_events()

def plot_function_2(fig, state):

   global qq, qqlimit
   global tempcache
   global pix, piy
   global hourticks1
   global caption1
   vmin1=-50
   vmax1=50
   Tsurf=np.array(state['surface_temperature']-273.15)
   lons1=np.mean(state['longitude'], axis=0)
   lats1=np.mean(state['latitude'], axis=1)
   zonetemps=np.array(np.mean( Tsurf, axis=1))
   lats2=np.array(np.mean( state['latitude'], axis=1))
   #print(zonetemps)
   #print(lats2)
   lats3=lats2*math.pi/10
   coslats=np.cos(lats2)
   weightemps=zonetemps*np.abs(coslats)*np.abs(coslats)
   meantemp1=round(np.mean(weightemps),2)
   meantemp0=round(np.mean(Tsurf),2)
   maxtemp0=round(np.max(Tsurf),2)
   mintemp0=round(np.min(Tsurf),2)
   deltemp0=maxtemp0-mintemp0        
   #print(" T mean min max ",meantemp0, mintemp0, maxtemp0)
   
   ax = fig.add_subplot(1, 1, 1)

   levels1=np.array([-250,-220,-200,-150,-120,-100,-80,-60,-40,-30,-20,-10,-5,0,5,10,15,20,25,30,35,40,45,50,60,70,80,90,100,120,150,180,200,250,300,400,500])
   Tsurf2 = scipy.ndimage.zoom(Tsurf, 4)

   #ax.contourf(state['longitude'], state['latitude'],Tsurf, cmap=mpl.colormaps['turbo'], origin='lower',levels=5,extent=[0,360,-90,90])
   ax.imshow(Tsurf2,interpolation="bicubic", origin='lower', vmin=vmin1, vmax=vmax1, cmap='turbo',extent=[0,360,-90,90])
   cs1=ax.contour(Tsurf2, alpha=0.5,colors="k",levels=levels1, origin='lower', vmin=vmin1, vmax=vmax1, extent=[0,360,-90,90])
   
   ax.clabel(cs1, cs1.levels, inline=True, fmt="%.1f", fontsize=16)
   ax.set_xlabel('Longitude', fontsize=16)
   ax.set_ylabel('Latitude', fontsize=16)	
   ax.set_xticklabels(labels=ax.get_xticklabels(),fontsize=12)
   ax.set_yticklabels(labels=ax.get_yticklabels(),fontsize=12)
   #ax.set_xticklabels(labels=list(range(0,360, 45)),fontsize=16)
   #fig.suptitle('Surface temperature '+str(round(meantemp,2))+' at time: '+str(state['time']))
   fig.suptitle(caption1,fontsize=18)
   #fig.subtitle('maxtemp '+str(maxtemp),fontsize=11)
   fig.canvas.draw()
   fig.canvas.flush_events()

1. climt.list_available_constants()

1. quit(-1)

climt.set_constants_from_dict({

   'stellar_irradiance': {'value': S1*1365.2, 'units': 'W m^-2'}})

climt.set_constants_from_dict({

   'obliquity': {'value': 0, 'units': 'degree'}})

climt.set_constants_from_dict({

   'omega_tilde': {'value': 0, 'units': 'degree'}})

climt.set_constants_from_dict({

   'lambda_m0': {'value': 0, 'units': 'degree'}})

1. gravitational_acceleration: 9.80665 m s^-2
2. planetary_radius: 6371000.0 m
3. planetary_rotation_rate: 7.292e-05 s^-1
4. seconds_per_day: 86400.0

1. 1. NON setting this 10 days
2. climt.set_constants_from_dict({'planetary_rotation_rate:': {'value': 7.292e-6/365.25, 'units': 's^-1'}})
3. climt.set_constants_from_dict({'seconds_per_day:': {'value': 864000*365.25, 'units': }})

1. climt.list_available_constants()

1. quit(-1)

1. 1. 4 subwin
2. monitor = PlotFunctionMonitor(plot_function)
   1. one lat, lon subwin

monitor = PlotFunctionMonitor(plot_function_2)

instellation = climt.Instellation() berger=climt.BergerSolarInsolation()

1. print(instellation)

1. quit(-1)

timestep = timedelta(hours=steppi)

convection = EmanuelConvection() radiation_sw = RRTMGShortwave() radiation_lw = RRTMGLongwave() slab = SlabSurface() simple_physics = SimplePhysics() dry_convection = DryConvectiveAdjustment()

hydrology=BucketHydrology()

opticaldepth=climt.Frierson06LongwaveOpticalDepth() gray=climt.GrayLongwaveRadiation() condensation=climt.GridScaleCondensation()

1. dcmip=climt.DcmipInitialConditions()

heldsuarez=climt.HeldSuarez()

ice = climt.IceSheet(maximum_snow_ice_height=30.) ## 30

1. state = get_default_state(
2. [simple_physics, convection, dry_convection,
3. radiation_lw, radiation_sw, slab]
4. )

1. state = get_default_state(
2. [instellation,
3. radiation_lw, radiation_sw,simple_physics,convection, dry_convection,hydrology, slab],
4. grid_state=climt.get_grid(nx=dimx, ny=dimy,latitude_grid='regular')
5. )

state = get_default_state(

   [instellation,berger,
    radiation_lw, radiation_sw,simple_physics,convection, dry_convection,hydrology, slab, ice],
grid_state=climt.get_grid(nx=dimx, ny=dimy,latitude_grid='regular')

)

1. 2. dycore

1. 1. dycore code ....

1. convection = climt.EmanuelConvection()
2. simple_physics = TimeDifferencingWrapper(climt.SimplePhysics())

1. radiation_step = timedelta(hours=1)

1. radiation_lw = UpdateFrequencyWrapper(
2. climt.RRTMGLongwave(), radiation_step)

1. radiation_sw = UpdateFrequencyWrapper(
2. climt.RRTMGShortwave(), radiation_step)

1. slab_surface = climt.SlabSurface()

1. dycore = gfs_dynamical_core.GFSDynamicalCore([simple_physics, slab_surface, radiation_sw,radiation_lw, convection], number_of_damped_levels=5)

1. grid = climt.get_grid(nx=32, ny=16, nz=8)

1. Create model state
2. state = climt.get_default_state([dycore], grid_state=grid)

co2=state['mole_fraction_of_carbon_dioxide_in_air'].values[:]

print("Co2 ", round(np.mean(co2)*1e6,3))

1. quit(-1)

co2=co2*co2coeff

1. quit(-1)

state['mole_fraction_of_carbon_dioxide_in_air'].values[:]=co2

state['air_temperature'].values[:] = 288 state['surface_albedo_for_direct_shortwave'].values[:] = albedo state['surface_albedo_for_direct_near_infrared'].values[:] = albedo state['surface_albedo_for_diffuse_shortwave'].values[:] = albedo state['surface_albedo_for_diffuse_near_infrared'].values[:] = albedo

1. state['zenith_angle'].values[:] = np.pi/2.5

state['surface_temperature'].values[:] = 288. state['ocean_mixed_layer_thickness'].values[:] = 0

state['heat_capacity_of_soil'][:]=750 state['surface_thermal_capacity'][:]=750 state['soil_layer_thickness'][:]=10 state['surface_material_density'][:]=1600 ## sand ## clay loam 1280

1. state['surface_material_density'][:]=2700 ## stone, basalt 2900

state['area_type'].values[:] = 'land'

1. state['area_type'].values[:] = 'sea'

1. state['area_type'].values[:] = 'sea_ice'
2. state['sea_ice_thickness'].values[:] = 2.
3. state['surface_snow_thickness'].values[:] = 2.
4. state['surface_temperature'].values[:] = 260.
5. state['surface_upward_sensible_heat_flux'].values[:] = -0.5

1. 1. clouds, testing!

1. state['mass_content_of_cloud_liquid_water_in_atmosphere_layer'].loc[dict(mid_levels=slice(4, 8))] = 0.03
2. state['cloud_area_fraction_in_atmosphere_layer'].loc[dict(mid_levels=slice(4, 8))] = 0.15.
   1. clouds, testing!

state['mass_content_of_cloud_liquid_water_in_atmosphere_layer'][:]= 0.1

1. state['cloud_area_fraction_in_atmosphere_layer'][:] = 1.0

state['cloud_area_fraction_in_atmosphere_layer'][:] = 0.0

print(state)

1. quit(-1)

1. print(state['time'])

1. datetime(year, month, day, hour, minute, second, microsecond)
2. newtime = datetime(2023, 7, 1, 11, 55, 59, 0)
3. newtime = datetime(2023, 6, 14, 11, 55, 59, 0) ## near mid summer

newtime = datetime(1, 3, 20, 0, 0, 0, 0) ## near spring equinox

state['time']=newtime

print(state['time'])

1. quit(-1)

diag = instellation(state)

print("Instellation")

print(diag)

1. quit(-1)

fields_to_store = ['air_temperature', 'air_pressure', 'eastward_wind',

                  'northward_wind', 'air_pressure_on_interface_levels',
                  'surface_pressure', 'upwelling_longwave_flux_in_air',
                  'specific_humidity', 'surface_temperature',
                  'latitude', 'longitude',
                  'convective_heating_rate']

netcdf_monitor = NetCDFMonitor('climt_out.nc',

                              write_on_store=True,
                              store_names=fields_to_store)

time_stepper = AdamsBashforth([radiation_lw, radiation_sw, slab])

1. 1. dycore runner
2. for i in range(1500*24):
3. diag, state = dycore(state, model_time_step)
4. state.update(diag)
5. state['time'] += model_time_step
7. if i % (20) == 0:
8. #netcdf_monitor.store(state)
9. monitor.store(state)
10. print('max. zonal wind: ', np.amax(state['eastward_wind'].values))
11. print('max. humidity: ', np.amax(state['specific_humidity'].values))
12. print('max. surf temp: ', np.amax(state['surface_temperature'].values))
14. print(state['time'])

print("Runnin simulation, insolation")

for i in range(32000000): diagnostics, state = time_stepper(state, timestep) state.update(diagnostics) diag=berger(state) state.update(diag) diag, out = simple_physics(state, timestep) state.update(diag) state.update(out) #diag, out = hydrology(state, timestep) #state.update(diag) #state.update(out) diag, out = dry_convection(state, timestep) state.update(diag) state.update(out) diag = instellation(state) if(i==0): czenit=diag['zenith_angle'].values[:] #if(i>0): state['zenith_angle'].values[:]=czenit if(i>0): diag['zenith_angle'].values[:]=czenit state.update(diag) diag, out = ice(state, timestep) state.update(diag) state.update(out)

if((i%10)==0): tempslice1=state['air_temperature'][0] save_netcdf_x("temperature1.nc", tempslice1) monitor.store(state) netcdf_monitor.store(state)

state['time'] += timestep