File:Aquaplanet temperature distribution 1 1 1 1.png

1. 3. climt gfs_dynamical_core aquaplanet
2. 1. python3 anaconda
   2. based on climt example
   4. 30.10.2023 v 0000.0000
   6. uses gfs_dynamical core
   8. gfs_dynamical_core
3. pip install gfs_dynamical core
4. 1. if from source download sit source tree
5. https://github.com/Ai33L/gfs_dynamical_core
   1. and mote to its directory "gfs_dynamical_core-develop"
6. and then in python 3.9 ## pip install .
   1. warning you must remoce climt and sympl first pip uninstall climt, pip uninstall sympl

import matplotlib.pyplot as plt import matplotlib as mpl import numpy as np import math

from datetime import timedelta from datetime import datetime

import climt import gfs_dynamical_core

from sympl import (

   TimeDifferencingWrapper, UpdateFrequencyWrapper,
   set_constant,
   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(1, 1, 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_function0(fig, state):

   ax = fig.add_subplot(2, 2, 1)
   state['specific_humidity'].mean(
       dim='lon').plot.contourf(
           ax=ax, levels=16, robust=True)
   ax.set_title('Specific Humidity')

   ax = fig.add_subplot(2, 2, 3)
   state['eastward_wind'].mean(dim='lon').plot.contourf(
       ax=ax, levels=16, robust=True)
   ax.set_title('Zonal Wind')

   ax = fig.add_subplot(2, 2, 2)
   state['air_temperature_tendency_from_convection'].transpose().mean(
       dim='lon').plot.contourf(
       ax=ax, levels=16, robust=True)
   ax.set_title('Conv. Heating Rate')

   ax = fig.add_subplot(2, 2, 4)
   state['air_temperature'].mean(dim='lon').plot.contourf(
       ax=ax, levels=16)
   ax.set_title('Temperature')

   fig.tight_layout()
   fig.canvas.draw()
   fig.canvas.flush_events()

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. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 3

1. 1. main code

S1=0.65 co2coeff=1.0 ## co2 x330

timestep_minutes=60*4

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'}})

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',
                  'convective_heating_rate', 'latitude', 'longitude']

1. Create plotting object

monitor = PlotFunctionMonitor(plot_function) netcdf_monitor = NetCDFMonitor('gcm_without_seasonal_cycle.nc',

                              write_on_store=True,
                              store_names=fields_to_store)

1. set_constant('stellar_irradiance', value=200, units='W m^-2')

model_time_step = timedelta(minutes=timestep_minutes)

1. Create components

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

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

radiation_step = timedelta(hours=1)

radiation_lw = UpdateFrequencyWrapper(

   climt.RRTMGLongwave(), radiation_step)

radiation_sw = UpdateFrequencyWrapper(

   climt.RRTMGShortwave(), radiation_step)

slab_surface = climt.SlabSurface()

convection2 = EmanuelConvection()

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. dycore=gfs_dynamical_core.GFSDynamicalCore(
2. [simple_physics, slab_surface, radiation_sw,
3. radiation_lw, convection], number_of_damped_levels=5
4. )

dycore=gfs_dynamical_core.GFSDynamicalCore(

   [simple_physics, slab_surface, radiation_sw,
    radiation_lw, convection], number_of_damped_levels=1

)

grid = climt.get_grid(nx=16, ny=16, nz=6)

1. Create model state

state = climt.get_default_state([instellation, berger, dycore], grid_state=grid)

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

state['time']=newtime

1. Set initial/boundary conditions

latitudes = state['latitude'].values longitudes = state['longitude'].values

zenith_angle = np.radians(latitudes) surface_shape = latitudes.shape

state['zenith_angle'].values = zenith_angle

state['eastward_wind'].values[:] = np.random.randn(

   *state['eastward_wind'].shape)

state['ocean_mixed_layer_thickness'].values[:] = 50

1. surf_temp_profile = 290 - (40*np.sin(zenith_angle)**2)

1. state['surface_temperature'].values = surf_temp_profile

state['surface_temperature'].values[:] = state['surface_temperature'].values[:]*0+273.15

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

print(" Running , it is really slow process ...")

for i in range(1500*24*6): diag, state = dycore(state, model_time_step) state.update(diag) diag=berger(state) state.update(diag) 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) state['time'] += model_time_step if ((i % 10) == 0): netcdf_monitor.store(state) monitor.store(state) print('max. zonal wind: ', np.amax(state['eastward_wind'].values)) print('max. humidity: ', np.amax(state['specific_humidity'].values)) print('max. surf temp: ', np.amax(state['surface_temperature'].values)) print(state['time'])