File:Cloudless desert planet surface temperature 1 1 1 1.png

1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 33
2. 1. climt planet test code 2, desert planet
3. 1. Importat note:
5. uses Ubuntu Anaconda Python 3.10.12 and gfs_dunamical_core
   1. must create 3.10 enviroment in conda, ant then: pip install gfs-dynamical-core
6. maybe nort function in other python versions!
7. 1. 09.10.2023 v 0000.0000c

import matplotlib as mpl import math import numpy as np

from datetime import timedelta from datetime import datetime

import climt

1. 1. basic params

steppi=1 # timestep hours dimx=16 ## dim of grid dimy=8

pix=8 ## temperature et al specimen point insige grid piy=4

S1=1 ## solar constant coeff

co2coeff=1.0 ## co2 x330 albedo=0.07 ## basalt surface like Moon

qq=0

qqlimit=int(24/steppi)

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

1. print (hourticks1)

1. quit(-1)

tempcache=np.zeros(qqlimit)+15

from sympl import (

   DataArray, PlotFunctionMonitor,NetCDFMonitor,
   AdamsBashforth

)

from climt import (

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

)

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()

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({
3. 'planetary_rotation_rate:': {'value': 7.292e-6, 'units': 's^-1'}})
4. climt.set_constants_from_dict({
5. 'seconds_per_day:': {'value': 864000, 'units': }})

climt.list_available_constants()

1. quit(-1)

monitor = PlotFunctionMonitor(plot_function)

instellation = climt.Instellation()

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()

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. )

state = get_default_state(

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

)

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[:] = 10

state['heat_capacity_of_soil'][:]=750 state['surface_thermal_capacity'][:]=750 state['soil_layer_thickness'][:]=0.1

1. state['surface_material_density'][:]=1600 ## sand ## clay loam 1280

state['surface_material_density'][:]=2700 ## stine, basalt 2900 state['area_type'].values[:] = 'land'

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.000

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)

newtime = datetime(2023, 6, 14, 11, 55, 59, 0) ## near mid summer

state['time']=newtime

print(state['time'])

1. quit(-1)

1. state['time']=3600*24*180

1. print(climt.get_constants_from_dict('obliquity'))
2. print(climt.constants['obliquity'])

1. print(get_constants_string())

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])

for i in range(320000):

   diagnostics, state = time_stepper(state, timestep)
   state.update(diagnostics)
   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)
   #state=new_state
   diag = instellation(state)
   state.update(diag)
   #diag, out = ice(state, timestep)
   #state.update(diag)
   #state.update(out)

   monitor.store(state)
   netcdf_monitor.store(state)
   #print('Zenith Angle at time: '+str(state['time']))
   state['time'] += timestep