Add thermodynamic sea ice code

exp/test_cases/thermodynamic_sea_ice/grey_thermo_ice_test_case.py

@@ -0,0 +1,238 @@
+import os
+
+import numpy as np
+
+from isca import GreyCodeBase, DiagTable, Experiment, Namelist, GFDL_BASE
+
+NCORES =16
+
+# Point to code as defined by $GFDL_BASE
+cb = GreyCodeBase.from_directory(GFDL_BASE)
+
+base_dir = os.path.dirname(os.path.realpath(__file__))
+
+cb.compile()  # compile the source code to working directory $GFDL_WORK/codebase
+
+# create an Experiment object to handle the configuration of model parameters
+# and output diagnostics
+exp = Experiment('thermo_ice_test_experiment', codebase=cb)
+
+
+#Tell model how to write diagnostics
+diag = DiagTable()
+diag.add_file('atmos_monthly', 30, 'days', time_units='days')
+
+#Tell model which diagnostics to write
+diag.add_field('dynamics', 'ps', time_avg=True)
+diag.add_field('dynamics', 'bk')
+diag.add_field('dynamics', 'pk')
+diag.add_field('atmosphere', 'temp_2m', time_avg=True)
+diag.add_field('dynamics', 'ucomp', time_avg=True)
+diag.add_field('dynamics', 'temp', time_avg=True)
+diag.add_field('mixed_layer', 'albedo', time_avg=True)
+diag.add_field('mixed_layer', 'h_thermo_ice', time_avg=True)
+diag.add_field('mixed_layer', 't_surf', time_avg=True)
+diag.add_field('mixed_layer', 'flux_oceanq', time_avg=True)
+diag.add_field('two_stream', 'swdn_toa', time_avg=True)
+
+
+
+exp.diag_table = diag # register diag table 
+
+
+#Empty the run directory ready to run
+exp.clear_rundir()
+
+#Define values for the 'core' namelist
+exp.namelist = namelist = Namelist({
+    'main_nml': {
+        'days'   : 360, # each run lasts one year, and then multiple runs are strung together below (loop on e.g. Line 310)
+        'hours'  : 0,   # a different output file is produced for each run (in this case, each year). Data 
+        'minutes': 0,   # output at the frequency specified in the diag table 
+        'seconds': 0,
+        'dt_atmos':900, 
+        'current_date' : [1,1,1,0,0,0],
+        'calendar' : 'thirty_day'
+    },
+
+    'idealized_moist_phys_nml': {
+        'two_stream_gray': True, 
+        'do_rrtm_radiation': False, 
+        'convection_scheme': 'SIMPLE_BETTS_MILLER',                                          
+        'do_damping': True,       
+        'turb':True,          
+        'mixed_layer_bc':True,                
+        'do_virtual' :True, 
+        'roughness_mom':5.e-3, 
+        'roughness_heat':1.e-5, 
+        'roughness_moist':1.e-5, 
+        'land_roughness_prefactor':1.0,
+        'do_simple':False,
+    },
+
+    'vert_turb_driver_nml': {
+        'do_mellor_yamada': False, 
+        'do_diffusivity': True, 
+        'do_edt':False, 
+        'constant_gust': 1.0,   
+        'use_tau': False, 
+        'do_entrain':False, 
+        'do_stable_bl':False, 
+        'do_shallow_conv':False, 
+        'do_simple':False 
+    },
+
+    'diffusivity_nml': {
+        'do_entrain':False, 
+        'entr_ratio': 0.0, 
+        'free_atm_diff':False, 
+        'do_simple': False, 
+        'parcel_buoy': 0.0, 
+        'frac_inner': 0.1,  
+        'fixed_depth': False, 
+    },
+
+    'surface_flux_nml': {
+        'use_virtual_temp': True, 
+        'do_simple': False, 
+        'old_dtaudv': False,  
+        'gust_const':1.0, 
+        'land_humidity_prefactor' : 1.0, 
+        'land_evap_prefactor': 1.0,
+    },
+
+    'atmosphere_nml': {
+        'idealized_moist_model': True 
+    },
+
+
+    'mixed_layer_nml': { 
+        'depth': 30., 
+        'albedo_value': 0.1, 
+        'prescribe_initial_dist':True, 
+        'tconst' : 305., 
+        'delta_T': 60., 
+        'evaporation':True, 
+        'do_qflux': True,
+        'load_qflux':False, 
+        'time_varying_qflux' : False, 
+        'do_thermo_ice':True, # turn on thermodynamic ice 
+        'thermo_ice_albedo':0.55, # ice albedo 
+        't_thermo_ice_base':273.15-2., # freezing point 
+        't_surf_freeze':273.15-2., # freezing point 
+        'do_var_thermo_ice_albedo':False,
+        'read_const_correct':False, 
+        'read_nudge_out':False, 
+
+    },
+
+    'qflux_nml':{
+        'qflux_amp':30.,
+    }, 
+
+
+
+    'qe_moist_convection_nml': { 
+        'rhbm':0.7, 
+        'tau_bm':7200., 
+        'Tmin':120., 
+        'Tmax':360., 
+        'val_inc': 0.01, 
+        'precision':1.e-6
+    },
+
+    'lscale_cond_nml': {
+        'do_simple':False, 
+        'do_evap':False 
+    },
+
+    'sat_vapor_pres_nml': {
+        'do_simple':False, 
+
+    },    
+
+
+
+    'damping_driver_nml': {
+        'do_rayleigh': True,
+        'trayfric': -0.5,              # neg. value: time in *days*
+        'sponge_pbottom':  2600.,
+        'do_conserve_energy': True,         
+    },
+
+
+
+
+    'two_stream_gray_rad_nml': {
+        'rad_scheme': 'frierson',            #Select radiation scheme to use
+        'do_seasonal': True,                
+        'use_time_average_coszen':True,
+        'dt_rad_avg':86400.,
+        'atm_abs': 0.22,                  
+        'solar_exponent':2,
+        'ir_tau_eq':7.2, 
+        'ir_tau_pole':3.6,#1.8, 
+        'del_sol':0.98, 
+        'solar_constant':1360, 
+        'linear_tau':0.2, 
+        'odp':1.4, 
+        'do_toa_albedo':True,
+    },
+
+
+    'spectral_dynamics_nml': {
+        'damping_order': 4, # Yields lap^8 damping 
+        'water_correction_limit': 200.e2, 
+        'reference_sea_level_press':1.0e5, 
+        'num_levels':30, 
+        'valid_range_t':[100.,800.],
+        'initial_sphum':[2.e-6], 
+        'use_virtual_temperature':True,
+        'vert_coord_option':'uneven_sigma', 
+        'robert_coeff':0.03, 
+        # set to T42 resolution 
+        'lon_max': 128, 
+        'lat_max': 64, 
+        'num_fourier': 42, 
+        'num_spherical':43, 
+        'surf_res': 0.05, 
+        'exponent': 3., 
+        'scale_heights': 5 
+
+    }, 
+
+    'spectral_init_cond_nml':{
+        'initial_temperature':280.
+    }, 
+
+
+
+    'diag_manager_nml': {
+        'mix_snapshot_average_fields': False  # time avg fields are labelled with time in middle of window
+    },
+
+    'fms_nml': {
+        'domains_stack_size': 600000                        # default: 0
+    },
+
+    'fms_io_nml': {
+        'threading_write': 'single',                         # default: multi
+        'fileset_write': 'single',                           # default: multi
+    },
+
+
+
+
+})
+
+# Lets do a run!
+if __name__=="__main__":
+
+
+    exp.run(1, use_restart=False, num_cores=NCORES, overwrite_data=False)
+
+    for i in range(1,11): # run for 10 years 
+        exp.run(i, num_cores=NCORES, overwrite_data=False)
+
+
+

src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90

@@ -138,6 +138,7 @@ module two_stream_gray_rad_mod
 character(len=256)                  :: co2_file='co2'       !  file name of co2 file to read
 character(len=256)                  :: co2_variable_name='co2'       !  file name of co2 file to read
 
+logical :: do_toa_albedo = .false.
 
 namelist/two_stream_gray_rad_nml/ solar_constant, del_sol, &
            ir_tau_eq, ir_tau_pole, odp, atm_abs, sw_diff, &
@@ -148,14 +149,15 @@ module two_stream_gray_rad_mod
 		   window, carbon_conc, rad_scheme, &
            do_read_co2, co2_file, co2_variable_name, solday, equinox_day, bog_a, bog_b, bog_mu, &
            use_time_average_coszen, dt_rad_avg,&
-           diabatic_acce !Schneider Liu values
+           diabatic_acce, & !Schneider Liu values 
+           do_toa_albedo 
 
 !==================================================================================
 !-------------------- diagnostics fields -------------------------------
 
 integer :: id_olr, id_swdn_sfc, id_swdn_toa, id_net_lw_surf, id_lwdn_sfc, id_lwup_sfc, &
            id_tdt_rad, id_tdt_solar, id_flux_rad, id_flux_lw, id_flux_sw, id_coszen, id_fracsun, &
-           id_lw_dtrans, id_lw_dtrans_win, id_sw_dtrans, id_co2
+           id_lw_dtrans, id_lw_dtrans_win, id_sw_dtrans, id_co2, id_nudge_flux
 
 character(len=10), parameter :: mod_name = 'two_stream'
 
@@ -347,7 +349,6 @@ subroutine two_stream_gray_rad_init(is, ie, js, je, num_levels, axes, Time, lonb
         register_diag_field ( mod_name, 'flux_sw', axes(half), Time, &
                'Net shortwave radiative flux (positive up)', &
                'W/m^2', missing_value=missing_value               )
-
     id_coszen  = &
                register_diag_field ( mod_name, 'coszen', axes(1:2), Time, &
                  'cosine of zenith angle', &
@@ -378,13 +379,17 @@ subroutine two_stream_gray_rad_init(is, ie, js, je, num_levels, axes, Time, lonb
                register_diag_field ( mod_name, 'lw_dtrans', axes(1:3), Time, &
                  'LW transmission (non window)', &
                  'none', missing_value=missing_value      )
+   id_nudge_flux = &
+   register_diag_field ( mod_name, 'nudge_flux', axes(1:2), Time, &
+               'Sea ice nudging flux', &
+               'W/m2', missing_value=missing_value               )
 return
 end subroutine two_stream_gray_rad_init
 
 ! ==================================================================================
 
 subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t,         &
-                           net_surf_sw_down, surf_lw_down, albedo, q)
+                           net_surf_sw_down, surf_lw_down, albedo, q, const_correct, nudge_out) 
 
 ! Begin the radiation calculation by computing downward fluxes.
 ! This part of the calculation does not depend on the surface temperature.
@@ -395,6 +400,7 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t,
 real, intent(out), dimension(:,:)   :: net_surf_sw_down
 real, intent(out), dimension(:,:)   :: surf_lw_down
 real, intent(in), dimension(:,:,:)  :: t, q,  p_half
+real, intent(in), dimension(:,:) :: const_correct, nudge_out
 integer :: i, j, k, n, dyofyr
 
 integer :: seconds, year_in_s, days
@@ -404,6 +410,7 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t,
 
 real ,dimension(size(q,1),size(q,2),size(q,3)) :: co2f
 
+
 n = size(t,3)
 
 
@@ -446,13 +453,21 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t,
   end if
 
      insolation = solar_constant * coszen
+
+  if (do_toa_albedo) then 
+      p2          = (1. - 3.*sin(lat)**2)/4.
+      insolation  = insolation * (0.75 + 0.15 * 2. * p2)
+  endif 
 
 else if (sw_scheme==B_SCHNEIDER_LIU) then
   insolation = (solar_constant/pi)*cos(lat)
 else
   ! Default: Averaged Earth insolation at all longitudes
   p2          = (1. - 3.*sin(lat)**2)/4.
   insolation  = 0.25 * solar_constant * (1.0 + del_sol * p2 + del_sw * sin(lat))
+  if (do_toa_albedo) then 
+      insolation  = insolation * (0.75 + 0.15 * 2. * p2)
+  endif 
 end if
 
 select case(sw_scheme)
@@ -476,21 +491,25 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t,
      sw_down(:,:,k+1)   = sw_down(:,:,k) * sw_dtrans(:,:,k)
   end do
 
+
 case(B_FRIERSON, B_BYRNE)
   ! Default: Frierson handling of SW radiation
   ! SW optical thickness
   sw_tau_0    = (1.0 - sw_diff*sin(lat)**2)*atm_abs
 
   ! compute optical depths for each model level
+
   do k = 1, n+1
     sw_tau(:,:,k) = sw_tau_0 * (p_half(:,:,k)/pstd_mks)**solar_exponent
   end do
 
   ! compute downward shortwave flux
+
   do k = 1, n+1
      sw_down(:,:,k)   = insolation(:,:) * exp(-sw_tau(:,:,k))
   end do
 
+
 case(B_SCHNEIDER_LIU)
   ! Schneider & Liu 2009 Giant planet scheme
   ! SW optical thickness
@@ -618,9 +637,9 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t,
 end select
 
 ! =================================================================================
-surf_lw_down     = lw_down(:, :, n+1)
+surf_lw_down     = lw_down(:, :, n+1) 
 toa_sw_in        = sw_down(:, :, 1)
-net_surf_sw_down = sw_down(:, :, n+1) * (1. - albedo)
+net_surf_sw_down = sw_down(:, :, n+1) * (1. - albedo) 
 ! =================================================================================
 
 if(lw_scheme.eq.B_SCHNEIDER_LIU) then
@@ -632,10 +651,18 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t,
 if ( id_lwdn_sfc > 0 ) then
    used = send_data ( id_lwdn_sfc, surf_lw_down, Time_diag)
 endif
+
 !------- incoming sw flux toa -------
 if ( id_swdn_toa > 0 ) then
    used = send_data ( id_swdn_toa, toa_sw_in, Time_diag)
 endif
+
+! NTL NUDGING
+if( id_nudge_flux > 0) then 
+   used = send_data ( id_nudge_flux, const_correct+nudge_out, Time_diag)
+endif 
+
+
 !------- downward sw flux surface -------
 if ( id_swdn_sfc > 0 ) then
    used = send_data ( id_swdn_sfc, net_surf_sw_down, Time_diag)