diff --git a/scm/etc/case_config/MOSAiC-AMPS.nml b/scm/etc/case_config/MOSAiC-AMPS.nml
new file mode 100644
index 00000000..ca3148a9
--- /dev/null
+++ b/scm/etc/case_config/MOSAiC-AMPS.nml
@@ -0,0 +1,17 @@
+$case_config
+case_name = 'MOSAiC-AMPS',
+runtime = 604800.0,
+thermo_forcing_type = 2,
+mom_forcing_type = 3,
+relax_time = 3600.0,
+sfc_flux_spec = .false.,
+sfc_type = 2,
+sfc_roughness_length_cm = 0.02,
+reference_profile_choice = 1,
+year = 2019,
+month = 10,
+day = 31,
+hour = 0,
+column_area = 2.0E9,
+lsm_ics = .true.,
+$end
diff --git a/scm/etc/case_config/MOSAiC-SS.nml b/scm/etc/case_config/MOSAiC-SS.nml
new file mode 100644
index 00000000..a2f5caa8
--- /dev/null
+++ b/scm/etc/case_config/MOSAiC-SS.nml
@@ -0,0 +1,17 @@
+$case_config
+case_name = 'MOSAiC-SS',
+runtime = 604800.0,
+thermo_forcing_type = 2,
+mom_forcing_type = 3,
+relax_time = 3600.0,
+sfc_flux_spec = .false.,
+sfc_type = 2,
+sfc_roughness_length_cm = 0.02,
+reference_profile_choice = 1,
+year = 2020,
+month = 3,
+day = 2,
+hour = 0,
+column_area = 2.0E9,
+lsm_ics = .true.,
+$end
diff --git a/scm/etc/scripts/MOSAiC_AMPS_forcing_file_generator.py b/scm/etc/scripts/MOSAiC_AMPS_forcing_file_generator.py
new file mode 100644
index 00000000..cfbb2d2b
--- /dev/null
+++ b/scm/etc/scripts/MOSAiC_AMPS_forcing_file_generator.py
@@ -0,0 +1,554 @@
+#!/usr/bin/env python
+
+from netCDF4 import Dataset
+import numpy as np
+import forcing_file_common as ffc
+import scipy.interpolate
+import scm_plotting_routines as spr
+
+#reload(ffc)
+
+#read in raw input file
+
+nc_fid = Dataset("../../data/raw_case_input/MOSAiC_31Oct20190Z_raw.nc", 'r')
+
+#ncdump to look at raw input file
+#nc_attrs, nc_dims, nc_vars = ncdump(nc_fid, False)
+
+#netCDF how-to
+#set attributes
+#file_id.setncattr(file_id.variables['variable'].ncattr(), nc_fid.variables['time'].getncattr(ncattr))
+#get attributes
+#file_id.variables['variable'].getncattr(index)
+
+#get raw input variables
+
+day = nc_fid.variables['day'][:]
+hour = nc_fid.variables['hour'][:]
+for t in range(day.size):
+    #find time index corresponding to October 29, 2019 at 0Z
+    if day[t] == 31 and hour[t] == 0:
+        start_t_index = t
+        break
+
+time = nc_fid.variables['time_offset'][:] #number of seconds since 00Z on 1/17/2006 (starts at 03Z)
+#subtract the initial time_offset from all values to get elapsed time since the start of the simulation
+time = time - time[start_t_index]
+levels = nc_fid.variables['levels'][:] #pressure levels in mb
+#convert levels to Pa
+levels = 100.0*levels
+ice_thickness = np.zeros((2),dtype=float)
+#height = nc_fid.variables['alt'][:]
+lat = nc_fid.variables['lat'][:] #degrees north
+#convert latitutde to degrees north
+lon = nc_fid.variables['lon'][:] #degrees east
+#use upstream T
+T_abs = nc_fid.variables['Tu'][:] #absolute temperature (time, lev)
+T_abs = np.swapaxes(T_abs, 0, 1)
+thetail = nc_fid.variables['thetail'][:] #theta_il (time, lev)
+thetail = np.swapaxes(thetail, 0, 1)
+#use unpstream theta instead of thetail
+thetailu = nc_fid.variables['thetailu'][:] #theta (time, lev)
+thetailu = np.swapaxes(thetailu, 0, 1)
+#calculate theta_il from absolute temperature (assuming no condensate)
+#thetal = np.zeros((levels.size,time.size),dtype=float)
+#for t in range(time.size):
+#    thetal[:,t] = (ffc.p0/levels)**(ffc.R_dry/ffc.c_p)*T_abs[:,t]
+qv = nc_fid.variables['qu'][:] #water vapor mixing ratio in kg/kg (time, lev)
+qv = np.swapaxes(qv, 0, 1)
+qt = np.zeros((levels.size,time.size),dtype=float)
+qt_mr = nc_fid.variables['qt'][:] #water vapor mixing ratio in kg/kg (time, lev)
+qt_mr = np.swapaxes(qt_mr, 0, 1) #swap the time and levels axis
+qt = qt_mr/(1.0 + qt_mr) #convert to specific humidity from mixing ratio
+qtu = np.zeros((levels.size,time.size),dtype=float)
+#use upstream qv instead of qt
+qtu_mr = nc_fid.variables['qtu'][:] #water vapor mixing ratio in kg/kg (time, lev)
+qtu_mr = np.swapaxes(qtu_mr, 0, 1) #swap the time and levels axis
+qtu = qtu_mr/(1.0 + qtu_mr) #convert to specific humidity from mixing ratio
+
+#ql and tke are not specified; set to zero
+#ql = np.zeros((levels.size,time.size),dtype=float)
+#qi = np.zeros((levels.size,time.size),dtype=float)
+ql = nc_fid.variables['ql'][:] #ql (time, lev)
+ql = np.swapaxes(ql, 0, 1)
+qi = nc_fid.variables['qi'][:] #ql (time, lev)
+qi = np.swapaxes(qi, 0, 1)
+tke = np.zeros((levels.size,time.size),dtype=float)
+# ozone_mmr = nc_fid.variables['o3mmr'][:]
+u_wind = nc_fid.variables['u'][:]
+u_wind = np.swapaxes(u_wind, 0, 1) #swap the time and levels axis
+v_wind = nc_fid.variables['v'][:]
+v_wind = np.swapaxes(v_wind, 0, 1) #swap the time and levels axis
+
+#w_sub = np.zeros((levels.size,time.size),dtype=float)
+omega = nc_fid.variables['omega'][:] #vertical pressure velocity in Pa/s
+omega = np.swapaxes(omega, 0, 1) #swap the time and levels axis
+w_sub = nc_fid.variables['w'][:] #vertical velocity in m/s
+w_sub = np.swapaxes(w_sub, 0, 1) #swap the time and levels axis
+#convert to w
+#for t in range(time.size):
+#    w_sub[:,t] = ffc.omega_to_w(omega[:,t],levels,T_abs[:,t])
+
+T_surf = nc_fid.variables['T_skin'][:]
+#T_surf = T_surf + 273.15 #convert to K
+#T_surf = (29 + 273.15)*np.ones((time.size),dtype=float) #forcing instructions specify time-invariant 29 deg C.
+p_surf = nc_fid.variables['p_srf'][:] #Pa
+#p_surf = p_surf*100.0 #convert to Pa
+
+#h_advec_thil = np.zeros((levels.size,time.size),dtype=float)
+h_advec_thil = nc_fid.variables['h_advec_thetail'][:] #K/s
+h_advec_thil = np.swapaxes(h_advec_thil, 0, 1) #swap the time and levels axis
+#for t in range(time.size):
+#    h_advec_thil[:,t] = (ffc.p0/levels)**(ffc.R_dry/ffc.c_p)*h_advec_T[:,t] #convert to potential temperature
+
+#v_advec_thil = np.zeros((levels.size,time.size),dtype=float)
+v_advec_thil = nc_fid.variables['v_advec_thetail'][:] #K/s
+v_advec_thil = np.swapaxes(v_advec_thil, 0, 1) #swap the time and levels axis
+#for t in range(time.size):
+#    v_advec_thil[:,t] = (ffc.p0/levels)**(ffc.R_dry/ffc.c_p)*v_advec_T[:,t] #convert to potential temperature
+
+dT_dt = np.zeros((levels.size,time.size),dtype=float)
+
+# h_advec_T = h_advec_T*86400.0
+# v_advec_T = v_advec_T*86400.0
+# dT_dt = dT_dt*86400.0
+# spr.contour_plot_firl(time, levels, h_advec_T, np.min(h_advec_T), np.max(h_advec_T), 'h advec T', 'time', 'pressure', 'h_advec_T.eps', y_inverted = True)
+# spr.contour_plot_firl(time, levels, v_advec_T, np.min(v_advec_T), np.max(v_advec_T), 'v advec T', 'time', 'pressure', 'v_advec_T.eps', y_inverted = True)
+# spr.contour_plot_firl(time, levels, dT_dt, np.min(dT_dt), np.max(dT_dt), 'total T tend', 'time', 'pressure', 'dT_dt.eps', y_inverted = True)
+
+#h_advec_qt = np.zeros((levels.size,time.size),dtype=float)
+h_advec_qt = nc_fid.variables['h_advec_qt'][:] #kg/kg/s
+h_advec_qt = np.swapaxes(h_advec_qt, 0, 1) #swap the time and levels axis
+h_advec_qt = h_advec_qt/(1.0 + qt_mr)**2 #convert to specific humidity
+
+h_advec_qi = np.zeros((levels.size,time.size),dtype=float)
+h_advec_ql = np.zeros((levels.size,time.size),dtype=float)
+
+#v_advec_qt = np.zeros((levels.size,time.size),dtype=float)
+v_advec_qt = nc_fid.variables['v_advec_qt'][:] #kg/kg/s
+v_advec_qt = np.swapaxes(v_advec_qt, 0, 1) #swap the time and levels axis
+v_advec_qt = v_advec_qt/(1.0 + qt_mr)**2 #convert to specific humidity
+
+dq_dt = np.zeros((levels.size,time.size),dtype=float)
+#dq_dt = nc_fid.variables['dqdt'][:]
+#dq_dt = np.swapaxes(dq_dt, 0, 1)*1.0E-3/3600.0 #swap the time and levels axis, convert to K/s
+
+#h_advec_qt = h_advec_qt*86400.0
+#v_advec_qt = v_advec_qt*86400.0
+#dq_dt = dq_dt*86400.0
+#spr.contour_plot_firl(time, levels, h_advec_qt, np.min(h_advec_qt), np.max(h_advec_qt), 'h advec q', 'time', 'pressure', 'h_advec_q.eps', y_inverted = True)
+#spr.contour_plot_firl(time, levels, v_advec_qt, np.min(v_advec_qt), np.max(v_advec_qt), 'v advec q', 'time', 'pressure', 'v_advec_q.eps', y_inverted = True)
+#spr.contour_plot_firl(time, levels, dq_dt, np.min(dq_dt), np.max(dq_dt), 'total q tend', 'time', 'pressure', 'dq_dt.eps', y_inverted = True)
+
+#phi_sfc = nc_fid.variables['phis'][:]
+#z_sfc = nc_fid.variables['alt'][:]
+height = nc_fid.variables['height'][:]
+#height = ffc.get_height_from_pres(T_abs[:,0],levels,z_sfc)
+
+#the following variables are not in this forcing file, but are included in other cases
+#rad_heating = nc_fid.variables['dT_dt_rad'][:] #K/s
+#rad_heating = np.swapaxes(rad_heating, 0, 1) #swap the time and levels axis
+rad_heating = np.zeros((levels.size,time.size),dtype=float)
+u_g = np.zeros((levels.size,time.size),dtype=float)
+v_g = np.zeros((levels.size,time.size),dtype=float)
+
+ozone = np.zeros((levels.size,time.size),dtype=float)
+ozone = ozone + 2.e-8
+
+tiice = np.zeros((ice_thickness.size),dtype=float)
+##tiice[0] = T_surf[start_t_index] 
+tiice[0] = 271.35 + .75*(T_surf[start_t_index] - 271.35)
+tiice[1] = 271.35 + .25*(T_surf[start_t_index] - 271.35)
+stc      = tiice[0]
+smc      = 0.33
+slc      = 0.33
+hice     = 0.3
+slmsk    = 2.0
+tsfco    = T_surf[start_t_index]
+weasd    = 200.0
+fice     = 1.0
+tisfc    = tsfco
+snwdph   = 2.e-4
+tg3      = 271.35
+zorl     = 15.0
+alvsf    = 0.23
+alnsf    = 0.23
+alvwf    = 0.23
+alnwf    = 0.23
+facsf    = 0.5055632
+facwf    = 0.4944368
+vegfrac  = 0.0
+canopy   = 0.0
+vegtyp   = 10
+soiltyp  = 12
+uustar   = 0.3828793
+shdmin   = 0.01
+shdmax   = 0.8
+slopetyp = 1
+snoalb   = 0.7287961
+
+# Open ozone file
+#f = open('../../data/raw_case_input/twpice_CRM_ozone.txt', 'r')
+
+# Read and ignore header lines
+#header1 = f.readline()
+
+#oz_pres = []
+#oz_data = []
+# Loop over lines and extract variables of interest
+#for line in f:
+#    line = line.strip()
+#    columns = line.split()
+#    oz_pres.append(float(columns[1]))
+#    oz_data.append(float(columns[2]))
+
+#f.close()
+
+#oz_pres = 100.0*np.array(oz_pres)
+#oz_data = np.array(oz_data)
+#oz_f = scipy.interpolate.interp1d(oz_pres, oz_data)
+#ozone_ppb = oz_f(levels[1:].tolist())
+#ozone_ppb = np.insert(ozone_ppb, 0, oz_data[0])
+#ozone_mmr = ozone_ppb*1.0E-9
+
+#
+#open processed input file for writing
+
+writefile_fid = Dataset('../../data/processed_case_input/MOSAiC.nc', 'w', format='NETCDF4')
+writefile_fid.description = "CCPP SCM forcing file for MOSAiC case"
+
+#create groups for scalars, intitialization, and forcing
+
+writefile_scalar_grp = writefile_fid.createGroup("scalars")
+writefile_initial_grp = writefile_fid.createGroup("initial")
+writefile_forcing_grp = writefile_fid.createGroup("forcing")
+
+#create dimensions and write them out
+
+writefile_time_dim = writefile_fid.createDimension('time', None)
+writefile_time_var = writefile_fid.createVariable('time', 'f4', ('time',))
+writefile_time_var[:] = time[start_t_index:]
+writefile_time_var.units = 's'
+writefile_time_var.description = 'elapsed time since the beginning of the simulation'
+
+writefile_levels_dim = writefile_fid.createDimension('levels', None)
+writefile_levels_var = writefile_fid.createVariable('levels', 'f4', ('levels',))
+writefile_levels_var[:] = levels
+writefile_levels_var.units = 'Pa'
+writefile_levels_var.description = 'pressure levels'
+
+writefile_ice_thickness_dim = writefile_fid.createDimension('ice_thickness', None)
+writefile_ice_thickness_var = writefile_fid.createVariable('ice_thickness', 'f4', ('ice_thickness',))
+writefile_ice_thickness_var[:] = ice_thickness
+writefile_ice_thickness_var.units = 'm'
+writefile_ice_thickness_var.description = 'depth of ice layers'
+
+#create variables and write them out
+
+#scalar group
+
+writefile_lat_var = writefile_scalar_grp.createVariable('lat', 'f4')
+writefile_lat_var[:] = lat
+writefile_lat_var.units = 'degrees N'
+writefile_lat_var.description = 'latitude of column'
+
+writefile_lon_var = writefile_scalar_grp.createVariable('lon', 'f4')
+writefile_lon_var[:] = lon
+writefile_lon_var.units = 'degrees E'
+writefile_lon_var.description = 'longitude of column'
+
+writefile_hice_var = writefile_scalar_grp.createVariable('hice', 'f4')
+writefile_hice_var[:] = hice
+writefile_hice_var.units = 'm'
+writefile_hice_var.description = 'sea ice thickness'
+
+writefile_slmsk_var = writefile_scalar_grp.createVariable('slmsk', 'f4')
+writefile_slmsk_var[:] = slmsk
+writefile_slmsk_var.units = ''
+writefile_slmsk_var.description = 'land-sea-ice mask'
+
+writefile_tsfco_var = writefile_scalar_grp.createVariable('tsfco', 'f4')
+writefile_tsfco_var[:] = tsfco
+writefile_tsfco_var.units = 'm'
+writefile_tsfco_var.description = 'sea ice surface skin temperature'
+
+writefile_weasd_var = writefile_scalar_grp.createVariable('weasd', 'f4')
+writefile_weasd_var[:] = weasd
+writefile_weasd_var.units = 'mm'
+writefile_weasd_var.description = 'water equivalent accumulated snow depth'
+
+writefile_fice_var = writefile_scalar_grp.createVariable('fice', 'f4')
+writefile_fice_var[:] = fice
+writefile_fice_var.units = '1'
+writefile_fice_var.description = 'ice fraction'
+
+writefile_tisfc_var = writefile_scalar_grp.createVariable('tisfc', 'f4')
+writefile_tisfc_var[:] = tisfc
+writefile_tisfc_var.units = 'K'
+writefile_tisfc_var.description = 'ice surface temperature'
+
+writefile_snwdph_var = writefile_scalar_grp.createVariable('snwdph', 'f4')
+writefile_snwdph_var[:] = snwdph
+writefile_snwdph_var.units = 'mm'
+writefile_snwdph_var.description = 'water equivalent snow depth'
+
+writefile_tg3_var = writefile_scalar_grp.createVariable('tg3', 'f4')
+writefile_tg3_var[:] = tg3
+writefile_tg3_var.units = 'K'
+writefile_tg3_var.description = 'deep soil temperature'
+
+writefile_zorl_var = writefile_scalar_grp.createVariable('zorl', 'f4')
+writefile_zorl_var[:] = zorl
+writefile_zorl_var.units = 'cm'
+writefile_zorl_var.description = 'composite surface roughness length'
+
+writefile_alvsf_var = writefile_scalar_grp.createVariable('alvsf', 'f4')
+writefile_alvsf_var[:] = alvsf
+writefile_alvsf_var.units = ''
+writefile_alvsf_var.description = '60 degree vis albedo with strong cosz dependency'
+
+writefile_alnsf_var = writefile_scalar_grp.createVariable('alnsf', 'f4')
+writefile_alnsf_var[:] = alnsf
+writefile_alnsf_var.units = ''
+writefile_alnsf_var.description = '60 degree nir albedo with strong cosz dependency'
+
+writefile_alvwf_var = writefile_scalar_grp.createVariable('alvwf', 'f4')
+writefile_alvwf_var[:] = alvwf
+writefile_alvwf_var.units = ''
+writefile_alvwf_var.description = '60 degree vis albedo with weak cosz dependency'
+
+writefile_alnwf_var = writefile_scalar_grp.createVariable('alnwf', 'f4')
+writefile_alnwf_var[:] = alnwf
+writefile_alnwf_var.units = ''
+writefile_alnwf_var.description = '60 degree nir albedo with weak cosz dependency'
+
+writefile_facsf_var = writefile_scalar_grp.createVariable('facsf', 'f4')
+writefile_facsf_var[:] = facsf
+writefile_facsf_var.units = ''
+writefile_facsf_var.description = 'fractional coverage with strong cosz dependency'
+
+writefile_facwf_var = writefile_scalar_grp.createVariable('facwf', 'f4')
+writefile_facwf_var[:] = facwf
+writefile_facwf_var.units = ''
+writefile_facwf_var.description = 'fractional coverage with weak cosz dependency'
+
+writefile_vegfrac_var = writefile_scalar_grp.createVariable('vegfrac', 'f4')
+writefile_vegfrac_var[:] = vegfrac
+writefile_vegfrac_var.units = '1'
+writefile_vegfrac_var.description = 'vegetation fraction'
+
+writefile_canopy_var = writefile_scalar_grp.createVariable('canopy', 'f4')
+writefile_canopy_var[:] = canopy
+writefile_canopy_var.units = 'kg m-2'
+writefile_canopy_var.description = 'amount of water stored in camopy'
+
+writefile_vegtyp_var = writefile_scalar_grp.createVariable('vegtyp', 'f4')
+writefile_vegtyp_var[:] = vegtyp
+writefile_vegtyp_var.units = ''
+writefile_vegtyp_var.description = 'vegetation type 1-12'
+
+writefile_soiltyp_var = writefile_scalar_grp.createVariable('soiltyp', 'f4')
+writefile_soiltyp_var[:] = soiltyp
+writefile_soiltyp_var.units = ''
+writefile_soiltyp_var.description = 'soil type 1-12'
+
+writefile_uustar_var = writefile_scalar_grp.createVariable('uustar', 'f4')
+writefile_uustar_var[:] = uustar
+writefile_uustar_var.units = 'm s-1'
+writefile_uustar_var.description = 'friction velocity'
+
+writefile_shdmin_var = writefile_scalar_grp.createVariable('shdmin', 'f4')
+writefile_shdmin_var[:] = shdmin
+writefile_shdmin_var.units = '1'
+writefile_shdmin_var.description = 'minimum vegetation fraction'
+
+writefile_shdmax_var = writefile_scalar_grp.createVariable('shdmax', 'f4')
+writefile_shdmax_var[:] = shdmax
+writefile_shdmax_var.units = '1'
+writefile_shdmax_var.description = 'maximum vegetation fraction'
+
+writefile_slopetyp_var = writefile_scalar_grp.createVariable('slopetyp', 'f4')
+writefile_slopetyp_var[:] = slopetyp
+writefile_slopetyp_var.units = ''
+writefile_slopetyp_var.description = 'slope type 1-9'
+
+writefile_snoalb_var = writefile_scalar_grp.createVariable('snoalb', 'f4')
+writefile_snoalb_var[:] = snoalb
+writefile_snoalb_var.units = '1'
+writefile_snoalb_var.description = 'maximum snow albedo'
+
+#initial group
+
+writefile_height_var = writefile_initial_grp.createVariable('height', 'f4', ('levels',))
+writefile_height_var[:] = height
+writefile_height_var.units = 'm'
+writefile_height_var.description = 'physical height at pressure levels'
+
+writefile_tiice_var = writefile_initial_grp.createVariable('tiice', 'f4', ('ice_thickness',))
+writefile_tiice_var[:] = tiice
+writefile_tiice_var.units = 'K'
+writefile_tiice_var.description = 'initial profile of sea ice internal temperature'
+
+writefile_stc_var = writefile_initial_grp.createVariable('stc', 'f4')
+writefile_stc_var[:] = stc
+writefile_stc_var.units = 'K'
+writefile_stc_var.description = 'initial profile of sea ice internal temperature'
+
+writefile_smc_var = writefile_initial_grp.createVariable('smc', 'f4')
+writefile_smc_var[:] = smc
+writefile_smc_var.units = 'm3 m-3'
+writefile_smc_var.description = 'initial profile of soil moisture'
+
+writefile_slc_var = writefile_initial_grp.createVariable('slc', 'f4')
+writefile_slc_var[:] = slc
+writefile_slc_var.units = 'm3 m-3'
+writefile_slc_var.description = 'initial profile of soil liquid water'
+
+writefile_thetail_var = writefile_initial_grp.createVariable('thetail', 'f4', ('levels',))
+writefile_thetail_var[:] = thetail[:,start_t_index]
+writefile_thetail_var.units = 'K'
+writefile_thetail_var.description = 'initial profile of ice-liquid water potential temperature'
+
+writefile_qt_var = writefile_initial_grp.createVariable('qt', 'f4', ('levels',))
+writefile_qt_var[:] = qt[:,start_t_index]
+writefile_qt_var.units = 'kg kg^-1'
+writefile_qt_var.description = 'initial profile of total water specific humidity'
+
+writefile_ql_var = writefile_initial_grp.createVariable('ql', 'f4', ('levels',))
+writefile_ql_var[:] = ql[:,start_t_index]
+writefile_ql_var.units = 'kg kg^-1'
+writefile_ql_var.description = 'initial profile of liquid water specific humidity'
+
+writefile_qi_var = writefile_initial_grp.createVariable('qi', 'f4', ('levels',))
+writefile_qi_var[:] = qi[:,start_t_index]
+writefile_qi_var.units = 'kg kg^-1'
+writefile_qi_var.description = 'initial profile of ice water specific humidity'
+
+writefile_u_var = writefile_initial_grp.createVariable('u', 'f4', ('levels',))
+writefile_u_var[:] = u_wind[:,start_t_index]
+writefile_u_var.units = 'm s^-1'
+writefile_u_var.description = 'initial profile of E-W horizontal wind'
+
+writefile_v_var = writefile_initial_grp.createVariable('v', 'f4', ('levels',))
+writefile_v_var[:] = v_wind[:,start_t_index]
+writefile_v_var.units = 'm s^-1'
+writefile_v_var.description = 'initial profile of N-S horizontal wind'
+
+writefile_tke_var = writefile_initial_grp.createVariable('tke', 'f4', ('levels',))
+writefile_tke_var[:] = tke[:,start_t_index]
+writefile_tke_var.units = 'm^2 s^-2'
+writefile_tke_var.description = 'initial profile of turbulence kinetic energy'
+
+writefile_ozone_var = writefile_initial_grp.createVariable('ozone', 'f4', ('levels',))
+writefile_ozone_var[:] = ozone[:,start_t_index]
+writefile_ozone_var.units = 'kg kg^-1'
+writefile_ozone_var.description = 'initial profile of ozone mass mixing ratio'
+
+#forcing group
+
+writefile_p_surf_var = writefile_forcing_grp.createVariable('p_surf', 'f4', ('time',))
+writefile_p_surf_var[:] = p_surf[start_t_index:]
+writefile_p_surf_var.units = 'Pa'
+writefile_p_surf_var.description = 'surface pressure'
+
+writefile_T_surf_var = writefile_forcing_grp.createVariable('T_surf', 'f4', ('time',))
+writefile_T_surf_var[:] = T_surf[start_t_index:]
+writefile_T_surf_var.units = 'K'
+writefile_T_surf_var.description = 'surface absolute temperature'
+
+writefile_w_ls_var = writefile_forcing_grp.createVariable('w_ls', 'f4', ('levels','time',))
+writefile_w_ls_var[:] = w_sub[:,start_t_index:]
+writefile_w_ls_var.units = 'm s^-1'
+writefile_w_ls_var.description = 'large scale vertical velocity'
+
+writefile_omega_var = writefile_forcing_grp.createVariable('omega', 'f4', ('levels','time',))
+writefile_omega_var[:] = omega[:,start_t_index:]
+writefile_omega_var.units = 'Pa s^-1'
+writefile_omega_var.description = 'large scale pressure vertical velocity'
+
+writefile_u_g_var = writefile_forcing_grp.createVariable('u_g', 'f4', ('levels','time',))
+writefile_u_g_var[:] = u_g[:,start_t_index:]
+writefile_u_g_var.units = 'm s^-1'
+writefile_u_g_var.description = 'large scale geostrophic E-W wind'
+
+writefile_v_g_var = writefile_forcing_grp.createVariable('v_g', 'f4', ('levels','time',))
+writefile_v_g_var[:] = v_g[:,start_t_index:]
+writefile_v_g_var.units = 'm s^-1'
+writefile_v_g_var.description = 'large scale geostrophic N-S wind'
+
+writefile_u_nudge_var = writefile_forcing_grp.createVariable('u_nudge', 'f4', ('levels','time',))
+writefile_u_nudge_var[:] = u_wind[:,start_t_index:]
+writefile_u_nudge_var.units = 'm s^-1'
+writefile_u_nudge_var.description = 'E-W wind to nudge toward'
+
+writefile_v_nudge_var = writefile_forcing_grp.createVariable('v_nudge', 'f4', ('levels','time',))
+writefile_v_nudge_var[:] = v_wind[:,start_t_index:]
+writefile_v_nudge_var.units = 'm s^-1'
+writefile_v_nudge_var.description = 'N-S wind to nudge toward'
+
+writefile_T_nudge_var = writefile_forcing_grp.createVariable('T_nudge', 'f4', ('levels','time',))
+writefile_T_nudge_var[:] = T_abs[:,start_t_index:]
+writefile_T_nudge_var.units = 'K'
+writefile_T_nudge_var.description = 'absolute temperature to nudge toward'
+
+writefile_thil_nudge_var = writefile_forcing_grp.createVariable('thil_nudge', 'f4', ('levels','time',))
+writefile_thil_nudge_var[:] = thetailu[:,start_t_index:]
+writefile_thil_nudge_var.units = 'K'
+writefile_thil_nudge_var.description = 'potential temperature to nudge toward'
+
+writefile_qt_nudge_var = writefile_forcing_grp.createVariable('qt_nudge', 'f4', ('levels','time',))
+writefile_qt_nudge_var[:] = qtu[:,start_t_index:]
+writefile_qt_nudge_var.units = 'kg kg^-1'
+writefile_qt_nudge_var.description = 'q_t to nudge toward'
+
+writefile_qi_nudge_var = writefile_forcing_grp.createVariable('qi_nudge', 'f4', ('levels','time',))
+writefile_qi_nudge_var[:] = qi[:,start_t_index:]
+writefile_qi_nudge_var.units = 'kg kg^-1'
+writefile_qi_nudge_var.description = 'q_i to nudge toward'
+
+writefile_ql_nudge_var = writefile_forcing_grp.createVariable('ql_nudge', 'f4', ('levels','time',))
+writefile_ql_nudge_var[:] = ql[:,start_t_index:]
+writefile_ql_nudge_var.units = 'kg kg^-1'
+writefile_ql_nudge_var.description = 'q_l to nudge toward'
+
+writefile_rad_heating_var = writefile_forcing_grp.createVariable('dT_dt_rad', 'f4', ('levels','time',))
+writefile_rad_heating_var[:] = rad_heating[:,start_t_index:]
+writefile_rad_heating_var.units = 'K s^-1'
+writefile_rad_heating_var.description = 'prescribed radiative heating rate'
+
+writefile_h_advec_thil_var = writefile_forcing_grp.createVariable('h_advec_thetail', 'f4', ('levels','time',))
+writefile_h_advec_thil_var[:] = h_advec_thil[:,start_t_index:]
+writefile_h_advec_thil_var.units = 'K s^-1'
+writefile_h_advec_thil_var.description = 'prescribed theta_il tendency due to horizontal advection'
+
+writefile_v_advec_thil_var = writefile_forcing_grp.createVariable('v_advec_thetail', 'f4', ('levels','time',))
+writefile_v_advec_thil_var[:] = v_advec_thil[:,start_t_index:]
+writefile_v_advec_thil_var.units = 'K s^-1'
+writefile_v_advec_thil_var.description = 'prescribed theta_il tendency due to vertical advection'
+
+writefile_h_advec_qt_var = writefile_forcing_grp.createVariable('h_advec_qt', 'f4', ('levels','time',))
+writefile_h_advec_qt_var[:] = h_advec_qt[:,start_t_index:]
+writefile_h_advec_qt_var.units = 'kg kg^-1 s^-1'
+writefile_h_advec_qt_var.description = 'prescribed q_t tendency due to horizontal advection'
+
+writefile_h_advec_qi_var = writefile_forcing_grp.createVariable('h_advec_qi', 'f4', ('levels','time',))
+writefile_h_advec_qi_var[:] = h_advec_qi[:,start_t_index:]
+writefile_h_advec_qi_var.units = 'kg kg^-1 s^-1'
+writefile_h_advec_qi_var.description = 'prescribed q_i tendency due to horizontal advection'
+
+writefile_h_advec_ql_var = writefile_forcing_grp.createVariable('h_advec_ql', 'f4', ('levels','time',))
+writefile_h_advec_ql_var[:] = h_advec_ql[:,start_t_index:]
+writefile_h_advec_ql_var.units = 'kg kg^-1 s^-1'
+writefile_h_advec_ql_var.description = 'prescribed q_l tendency due to horizontal advection'
+
+writefile_v_advec_qt_var = writefile_forcing_grp.createVariable('v_advec_qt', 'f4', ('levels','time',))
+writefile_v_advec_qt_var[:] = v_advec_qt[:,start_t_index:]
+writefile_v_advec_qt_var.units = 'kg kg^-1 s^-1'
+writefile_v_advec_qt_var.description = 'prescribed q_t tendency due to vertical advection'
+
+
+#close processed input file
+writefile_fid.close()
+
+#close raw input file
+nc_fid.close()
diff --git a/scm/etc/scripts/MOSAiC_SS_forcing_file_generator.py b/scm/etc/scripts/MOSAiC_SS_forcing_file_generator.py
new file mode 100644
index 00000000..918eea33
--- /dev/null
+++ b/scm/etc/scripts/MOSAiC_SS_forcing_file_generator.py
@@ -0,0 +1,554 @@
+#!/usr/bin/env python
+
+from netCDF4 import Dataset
+import numpy as np
+import forcing_file_common as ffc
+import scipy.interpolate
+import scm_plotting_routines as spr
+
+#reload(ffc)
+
+#read in raw input file
+
+nc_fid = Dataset("../../data/raw_case_input/MOSAiC_2Mar20200Z_raw.nc", 'r')
+
+#ncdump to look at raw input file
+#nc_attrs, nc_dims, nc_vars = ncdump(nc_fid, False)
+
+#netCDF how-to
+#set attributes
+#file_id.setncattr(file_id.variables['variable'].ncattr(), nc_fid.variables['time'].getncattr(ncattr))
+#get attributes
+#file_id.variables['variable'].getncattr(index)
+
+#get raw input variables
+
+day = nc_fid.variables['day'][:]
+hour = nc_fid.variables['hour'][:]
+for t in range(day.size):
+    #find time index corresponding to October 29, 2019 at 0Z
+    if day[t] == 2 and hour[t] == 0:
+        start_t_index = t
+        break
+
+time = nc_fid.variables['time_offset'][:] #number of seconds since 00Z on 1/17/2006 (starts at 03Z)
+#subtract the initial time_offset from all values to get elapsed time since the start of the simulation
+time = time - time[start_t_index]
+levels = nc_fid.variables['levels'][:] #pressure levels in mb
+#convert levels to Pa
+levels = 100.0*levels
+ice_thickness = np.zeros((2),dtype=float)
+#height = nc_fid.variables['alt'][:]
+lat = nc_fid.variables['lat'][:] #degrees north
+#convert latitutde to degrees north
+lon = nc_fid.variables['lon'][:] #degrees east
+#use upstream T
+T_abs = nc_fid.variables['Tu'][:] #absolute temperature (time, lev)
+T_abs = np.swapaxes(T_abs, 0, 1)
+thetailu = nc_fid.variables['thetailu'][:] #theta_il (time, lev)
+thetailu = np.swapaxes(thetailu, 0, 1)
+#use unpstream theta instead of thetail
+thetail = nc_fid.variables['thetail'][:] #theta (time, lev)
+thetail = np.swapaxes(thetail, 0, 1)
+#calculate theta_il from absolute temperature (assuming no condensate)
+#thetal = np.zeros((levels.size,time.size),dtype=float)
+#for t in range(time.size):
+#    thetal[:,t] = (ffc.p0/levels)**(ffc.R_dry/ffc.c_p)*T_abs[:,t]
+qv = nc_fid.variables['qu'][:] #water vapor mixing ratio in kg/kg (time, lev)
+qv = np.swapaxes(qv, 0, 1)
+qt = np.zeros((levels.size,time.size),dtype=float)
+qt_mr = nc_fid.variables['qt'][:] #water vapor mixing ratio in kg/kg (time, lev)
+qt_mr = np.swapaxes(qt_mr, 0, 1) #swap the time and levels axis
+qt = qt_mr/(1.0 + qt_mr) #convert to specific humidity from mixing ratio
+qtu = np.zeros((levels.size,time.size),dtype=float)
+#use upstream qv instead of qt
+qtu_mr = nc_fid.variables['qtu'][:] #water vapor mixing ratio in kg/kg (time, lev)
+qtu_mr = np.swapaxes(qtu_mr, 0, 1) #swap the time and levels axis
+qtu = qtu_mr/(1.0 + qtu_mr) #convert to specific humidity from mixing ratio
+
+#ql and tke are not specified; set to zero
+#ql = np.zeros((levels.size,time.size),dtype=float)
+#qi = np.zeros((levels.size,time.size),dtype=float)
+ql = nc_fid.variables['ql'][:] #ql (time, lev)
+ql = np.swapaxes(ql, 0, 1)
+qi = nc_fid.variables['qi'][:] #ql (time, lev)
+qi = np.swapaxes(qi, 0, 1)
+tke = np.zeros((levels.size,time.size),dtype=float)
+# ozone_mmr = nc_fid.variables['o3mmr'][:]
+u_wind = nc_fid.variables['u'][:]
+u_wind = np.swapaxes(u_wind, 0, 1) #swap the time and levels axis
+v_wind = nc_fid.variables['v'][:]
+v_wind = np.swapaxes(v_wind, 0, 1) #swap the time and levels axis
+
+#w_sub = np.zeros((levels.size,time.size),dtype=float)
+omega = nc_fid.variables['omega'][:] #vertical pressure velocity in Pa/s
+omega = np.swapaxes(omega, 0, 1) #swap the time and levels axis
+w_sub = nc_fid.variables['w'][:] #vertical velocity in m/s
+w_sub = np.swapaxes(w_sub, 0, 1) #swap the time and levels axis
+#convert to w
+#for t in range(time.size):
+#    w_sub[:,t] = ffc.omega_to_w(omega[:,t],levels,T_abs[:,t])
+
+T_surf = nc_fid.variables['T_skin'][:]
+#T_surf = T_surf + 273.15 #convert to K
+#T_surf = (29 + 273.15)*np.ones((time.size),dtype=float) #forcing instructions specify time-invariant 29 deg C.
+p_surf = nc_fid.variables['p_srf'][:] #Pa
+#p_surf = p_surf*100.0 #convert to Pa
+
+#h_advec_thil = np.zeros((levels.size,time.size),dtype=float)
+h_advec_thil = nc_fid.variables['h_advec_thetail'][:] #K/s
+h_advec_thil = np.swapaxes(h_advec_thil, 0, 1) #swap the time and levels axis
+#for t in range(time.size):
+#    h_advec_thil[:,t] = (ffc.p0/levels)**(ffc.R_dry/ffc.c_p)*h_advec_T[:,t] #convert to potential temperature
+
+#v_advec_thil = np.zeros((levels.size,time.size),dtype=float)
+v_advec_thil = nc_fid.variables['v_advec_thetail'][:] #K/s
+v_advec_thil = np.swapaxes(v_advec_thil, 0, 1) #swap the time and levels axis
+#for t in range(time.size):
+#    v_advec_thil[:,t] = (ffc.p0/levels)**(ffc.R_dry/ffc.c_p)*v_advec_T[:,t] #convert to potential temperature
+
+dT_dt = np.zeros((levels.size,time.size),dtype=float)
+
+# h_advec_T = h_advec_T*86400.0
+# v_advec_T = v_advec_T*86400.0
+# dT_dt = dT_dt*86400.0
+# spr.contour_plot_firl(time, levels, h_advec_T, np.min(h_advec_T), np.max(h_advec_T), 'h advec T', 'time', 'pressure', 'h_advec_T.eps', y_inverted = True)
+# spr.contour_plot_firl(time, levels, v_advec_T, np.min(v_advec_T), np.max(v_advec_T), 'v advec T', 'time', 'pressure', 'v_advec_T.eps', y_inverted = True)
+# spr.contour_plot_firl(time, levels, dT_dt, np.min(dT_dt), np.max(dT_dt), 'total T tend', 'time', 'pressure', 'dT_dt.eps', y_inverted = True)
+
+#h_advec_qt = np.zeros((levels.size,time.size),dtype=float)
+h_advec_qt = nc_fid.variables['h_advec_qt'][:] #kg/kg/s
+h_advec_qt = np.swapaxes(h_advec_qt, 0, 1) #swap the time and levels axis
+h_advec_qt = h_advec_qt/(1.0 + qt_mr)**2 #convert to specific humidity
+
+h_advec_qi = np.zeros((levels.size,time.size),dtype=float)
+h_advec_ql = np.zeros((levels.size,time.size),dtype=float)
+
+#v_advec_qt = np.zeros((levels.size,time.size),dtype=float)
+v_advec_qt = nc_fid.variables['v_advec_qt'][:] #kg/kg/s
+v_advec_qt = np.swapaxes(v_advec_qt, 0, 1) #swap the time and levels axis
+v_advec_qt = v_advec_qt/(1.0 + qt_mr)**2 #convert to specific humidity
+
+dq_dt = np.zeros((levels.size,time.size),dtype=float)
+#dq_dt = nc_fid.variables['dqdt'][:]
+#dq_dt = np.swapaxes(dq_dt, 0, 1)*1.0E-3/3600.0 #swap the time and levels axis, convert to K/s
+
+#h_advec_qt = h_advec_qt*86400.0
+#v_advec_qt = v_advec_qt*86400.0
+#dq_dt = dq_dt*86400.0
+#spr.contour_plot_firl(time, levels, h_advec_qt, np.min(h_advec_qt), np.max(h_advec_qt), 'h advec q', 'time', 'pressure', 'h_advec_q.eps', y_inverted = True)
+#spr.contour_plot_firl(time, levels, v_advec_qt, np.min(v_advec_qt), np.max(v_advec_qt), 'v advec q', 'time', 'pressure', 'v_advec_q.eps', y_inverted = True)
+#spr.contour_plot_firl(time, levels, dq_dt, np.min(dq_dt), np.max(dq_dt), 'total q tend', 'time', 'pressure', 'dq_dt.eps', y_inverted = True)
+
+#phi_sfc = nc_fid.variables['phis'][:]
+#z_sfc = nc_fid.variables['alt'][:]
+height = nc_fid.variables['height'][:]
+#height = ffc.get_height_from_pres(T_abs[:,0],levels,z_sfc)
+
+#the following variables are not in this forcing file, but are included in other cases
+#rad_heating = nc_fid.variables['dT_dt_rad'][:] #K/s
+#rad_heating = np.swapaxes(rad_heating, 0, 1) #swap the time and levels axis
+rad_heating = np.zeros((levels.size,time.size),dtype=float)
+u_g = np.zeros((levels.size,time.size),dtype=float)
+v_g = np.zeros((levels.size,time.size),dtype=float)
+
+ozone = np.zeros((levels.size,time.size),dtype=float)
+ozone = ozone + 2.e-8
+
+tiice = np.zeros((ice_thickness.size),dtype=float)
+##tiice[0] = T_surf[start_t_index] 
+tiice[0] = 271.35 + .75*(T_surf[start_t_index] - 271.35)
+tiice[1] = 271.35 + .25*(T_surf[start_t_index] - 271.35)
+stc      = tiice[0]
+smc      = 0.33
+slc      = 0.33
+hice     = 0.3
+slmsk    = 2.0
+tsfco    = T_surf[start_t_index]
+weasd    = 200.0
+fice     = 1.0
+tisfc    = tsfco
+snwdph   = 2.e-4
+tg3      = 271.35
+zorl     = 15.0
+alvsf    = 0.23
+alnsf    = 0.23
+alvwf    = 0.23
+alnwf    = 0.23
+facsf    = 0.5055632
+facwf    = 0.4944368
+vegfrac  = 0.0
+canopy   = 0.0
+vegtyp   = 10
+soiltyp  = 12
+uustar   = 0.3828793
+shdmin   = 0.01
+shdmax   = 0.8
+slopetyp = 1
+snoalb   = 0.7287961
+
+# Open ozone file
+#f = open('../../data/raw_case_input/twpice_CRM_ozone.txt', 'r')
+
+# Read and ignore header lines
+#header1 = f.readline()
+
+#oz_pres = []
+#oz_data = []
+# Loop over lines and extract variables of interest
+#for line in f:
+#    line = line.strip()
+#    columns = line.split()
+#    oz_pres.append(float(columns[1]))
+#    oz_data.append(float(columns[2]))
+
+#f.close()
+
+#oz_pres = 100.0*np.array(oz_pres)
+#oz_data = np.array(oz_data)
+#oz_f = scipy.interpolate.interp1d(oz_pres, oz_data)
+#ozone_ppb = oz_f(levels[1:].tolist())
+#ozone_ppb = np.insert(ozone_ppb, 0, oz_data[0])
+#ozone_mmr = ozone_ppb*1.0E-9
+
+#
+#open processed input file for writing
+
+writefile_fid = Dataset('../../data/processed_case_input/MOSAiC.nc', 'w', format='NETCDF4')
+writefile_fid.description = "CCPP SCM forcing file for MOSAiC case"
+
+#create groups for scalars, intitialization, and forcing
+
+writefile_scalar_grp = writefile_fid.createGroup("scalars")
+writefile_initial_grp = writefile_fid.createGroup("initial")
+writefile_forcing_grp = writefile_fid.createGroup("forcing")
+
+#create dimensions and write them out
+
+writefile_time_dim = writefile_fid.createDimension('time', None)
+writefile_time_var = writefile_fid.createVariable('time', 'f4', ('time',))
+writefile_time_var[:] = time[start_t_index:]
+writefile_time_var.units = 's'
+writefile_time_var.description = 'elapsed time since the beginning of the simulation'
+
+writefile_levels_dim = writefile_fid.createDimension('levels', None)
+writefile_levels_var = writefile_fid.createVariable('levels', 'f4', ('levels',))
+writefile_levels_var[:] = levels
+writefile_levels_var.units = 'Pa'
+writefile_levels_var.description = 'pressure levels'
+
+writefile_ice_thickness_dim = writefile_fid.createDimension('ice_thickness', None)
+writefile_ice_thickness_var = writefile_fid.createVariable('ice_thickness', 'f4', ('ice_thickness',))
+writefile_ice_thickness_var[:] = ice_thickness
+writefile_ice_thickness_var.units = 'm'
+writefile_ice_thickness_var.description = 'depth of ice layers'
+
+#create variables and write them out
+
+#scalar group
+
+writefile_lat_var = writefile_scalar_grp.createVariable('lat', 'f4')
+writefile_lat_var[:] = lat
+writefile_lat_var.units = 'degrees N'
+writefile_lat_var.description = 'latitude of column'
+
+writefile_lon_var = writefile_scalar_grp.createVariable('lon', 'f4')
+writefile_lon_var[:] = lon
+writefile_lon_var.units = 'degrees E'
+writefile_lon_var.description = 'longitude of column'
+
+writefile_hice_var = writefile_scalar_grp.createVariable('hice', 'f4')
+writefile_hice_var[:] = hice
+writefile_hice_var.units = 'm'
+writefile_hice_var.description = 'sea ice thickness'
+
+writefile_slmsk_var = writefile_scalar_grp.createVariable('slmsk', 'f4')
+writefile_slmsk_var[:] = slmsk
+writefile_slmsk_var.units = ''
+writefile_slmsk_var.description = 'land-sea-ice mask'
+
+writefile_tsfco_var = writefile_scalar_grp.createVariable('tsfco', 'f4')
+writefile_tsfco_var[:] = tsfco
+writefile_tsfco_var.units = 'm'
+writefile_tsfco_var.description = 'sea ice surface skin temperature'
+
+writefile_weasd_var = writefile_scalar_grp.createVariable('weasd', 'f4')
+writefile_weasd_var[:] = weasd
+writefile_weasd_var.units = 'mm'
+writefile_weasd_var.description = 'water equivalent accumulated snow depth'
+
+writefile_fice_var = writefile_scalar_grp.createVariable('fice', 'f4')
+writefile_fice_var[:] = fice
+writefile_fice_var.units = '1'
+writefile_fice_var.description = 'ice fraction'
+
+writefile_tisfc_var = writefile_scalar_grp.createVariable('tisfc', 'f4')
+writefile_tisfc_var[:] = tisfc
+writefile_tisfc_var.units = 'K'
+writefile_tisfc_var.description = 'ice surface temperature'
+
+writefile_snwdph_var = writefile_scalar_grp.createVariable('snwdph', 'f4')
+writefile_snwdph_var[:] = snwdph
+writefile_snwdph_var.units = 'mm'
+writefile_snwdph_var.description = 'water equivalent snow depth'
+
+writefile_tg3_var = writefile_scalar_grp.createVariable('tg3', 'f4')
+writefile_tg3_var[:] = tg3
+writefile_tg3_var.units = 'K'
+writefile_tg3_var.description = 'deep soil temperature'
+
+writefile_zorl_var = writefile_scalar_grp.createVariable('zorl', 'f4')
+writefile_zorl_var[:] = zorl
+writefile_zorl_var.units = 'cm'
+writefile_zorl_var.description = 'composite surface roughness length'
+
+writefile_alvsf_var = writefile_scalar_grp.createVariable('alvsf', 'f4')
+writefile_alvsf_var[:] = alvsf
+writefile_alvsf_var.units = ''
+writefile_alvsf_var.description = '60 degree vis albedo with strong cosz dependency'
+
+writefile_alnsf_var = writefile_scalar_grp.createVariable('alnsf', 'f4')
+writefile_alnsf_var[:] = alnsf
+writefile_alnsf_var.units = ''
+writefile_alnsf_var.description = '60 degree nir albedo with strong cosz dependency'
+
+writefile_alvwf_var = writefile_scalar_grp.createVariable('alvwf', 'f4')
+writefile_alvwf_var[:] = alvwf
+writefile_alvwf_var.units = ''
+writefile_alvwf_var.description = '60 degree vis albedo with weak cosz dependency'
+
+writefile_alnwf_var = writefile_scalar_grp.createVariable('alnwf', 'f4')
+writefile_alnwf_var[:] = alnwf
+writefile_alnwf_var.units = ''
+writefile_alnwf_var.description = '60 degree nir albedo with weak cosz dependency'
+
+writefile_facsf_var = writefile_scalar_grp.createVariable('facsf', 'f4')
+writefile_facsf_var[:] = facsf
+writefile_facsf_var.units = ''
+writefile_facsf_var.description = 'fractional coverage with strong cosz dependency'
+
+writefile_facwf_var = writefile_scalar_grp.createVariable('facwf', 'f4')
+writefile_facwf_var[:] = facwf
+writefile_facwf_var.units = ''
+writefile_facwf_var.description = 'fractional coverage with weak cosz dependency'
+
+writefile_vegfrac_var = writefile_scalar_grp.createVariable('vegfrac', 'f4')
+writefile_vegfrac_var[:] = vegfrac
+writefile_vegfrac_var.units = '1'
+writefile_vegfrac_var.description = 'vegetation fraction'
+
+writefile_canopy_var = writefile_scalar_grp.createVariable('canopy', 'f4')
+writefile_canopy_var[:] = canopy
+writefile_canopy_var.units = 'kg m-2'
+writefile_canopy_var.description = 'amount of water stored in camopy'
+
+writefile_vegtyp_var = writefile_scalar_grp.createVariable('vegtyp', 'f4')
+writefile_vegtyp_var[:] = vegtyp
+writefile_vegtyp_var.units = ''
+writefile_vegtyp_var.description = 'vegetation type 1-12'
+
+writefile_soiltyp_var = writefile_scalar_grp.createVariable('soiltyp', 'f4')
+writefile_soiltyp_var[:] = soiltyp
+writefile_soiltyp_var.units = ''
+writefile_soiltyp_var.description = 'soil type 1-12'
+
+writefile_uustar_var = writefile_scalar_grp.createVariable('uustar', 'f4')
+writefile_uustar_var[:] = uustar
+writefile_uustar_var.units = 'm s-1'
+writefile_uustar_var.description = 'friction velocity'
+
+writefile_shdmin_var = writefile_scalar_grp.createVariable('shdmin', 'f4')
+writefile_shdmin_var[:] = shdmin
+writefile_shdmin_var.units = '1'
+writefile_shdmin_var.description = 'minimum vegetation fraction'
+
+writefile_shdmax_var = writefile_scalar_grp.createVariable('shdmax', 'f4')
+writefile_shdmax_var[:] = shdmax
+writefile_shdmax_var.units = '1'
+writefile_shdmax_var.description = 'maximum vegetation fraction'
+
+writefile_slopetyp_var = writefile_scalar_grp.createVariable('slopetyp', 'f4')
+writefile_slopetyp_var[:] = slopetyp
+writefile_slopetyp_var.units = ''
+writefile_slopetyp_var.description = 'slope type 1-9'
+
+writefile_snoalb_var = writefile_scalar_grp.createVariable('snoalb', 'f4')
+writefile_snoalb_var[:] = snoalb
+writefile_snoalb_var.units = '1'
+writefile_snoalb_var.description = 'maximum snow albedo'
+
+#initial group
+
+writefile_height_var = writefile_initial_grp.createVariable('height', 'f4', ('levels',))
+writefile_height_var[:] = height
+writefile_height_var.units = 'm'
+writefile_height_var.description = 'physical height at pressure levels'
+
+writefile_tiice_var = writefile_initial_grp.createVariable('tiice', 'f4', ('ice_thickness',))
+writefile_tiice_var[:] = tiice
+writefile_tiice_var.units = 'K'
+writefile_tiice_var.description = 'initial profile of sea ice internal temperature'
+
+writefile_stc_var = writefile_initial_grp.createVariable('stc', 'f4')
+writefile_stc_var[:] = stc
+writefile_stc_var.units = 'K'
+writefile_stc_var.description = 'initial profile of sea ice internal temperature'
+
+writefile_smc_var = writefile_initial_grp.createVariable('smc', 'f4')
+writefile_smc_var[:] = smc
+writefile_smc_var.units = 'm3 m-3'
+writefile_smc_var.description = 'initial profile of soil moisture'
+
+writefile_slc_var = writefile_initial_grp.createVariable('slc', 'f4')
+writefile_slc_var[:] = slc
+writefile_slc_var.units = 'm3 m-3'
+writefile_slc_var.description = 'initial profile of soil liquid water'
+
+writefile_thetail_var = writefile_initial_grp.createVariable('thetail', 'f4', ('levels',))
+writefile_thetail_var[:] = thetail[:,start_t_index]
+writefile_thetail_var.units = 'K'
+writefile_thetail_var.description = 'initial profile of ice-liquid water potential temperature'
+
+writefile_qt_var = writefile_initial_grp.createVariable('qt', 'f4', ('levels',))
+writefile_qt_var[:] = qt[:,start_t_index]
+writefile_qt_var.units = 'kg kg^-1'
+writefile_qt_var.description = 'initial profile of total water specific humidity'
+
+writefile_ql_var = writefile_initial_grp.createVariable('ql', 'f4', ('levels',))
+writefile_ql_var[:] = ql[:,start_t_index]
+writefile_ql_var.units = 'kg kg^-1'
+writefile_ql_var.description = 'initial profile of liquid water specific humidity'
+
+writefile_qi_var = writefile_initial_grp.createVariable('qi', 'f4', ('levels',))
+writefile_qi_var[:] = qi[:,start_t_index]
+writefile_qi_var.units = 'kg kg^-1'
+writefile_qi_var.description = 'initial profile of ice water specific humidity'
+
+writefile_u_var = writefile_initial_grp.createVariable('u', 'f4', ('levels',))
+writefile_u_var[:] = u_wind[:,start_t_index]
+writefile_u_var.units = 'm s^-1'
+writefile_u_var.description = 'initial profile of E-W horizontal wind'
+
+writefile_v_var = writefile_initial_grp.createVariable('v', 'f4', ('levels',))
+writefile_v_var[:] = v_wind[:,start_t_index]
+writefile_v_var.units = 'm s^-1'
+writefile_v_var.description = 'initial profile of N-S horizontal wind'
+
+writefile_tke_var = writefile_initial_grp.createVariable('tke', 'f4', ('levels',))
+writefile_tke_var[:] = tke[:,start_t_index]
+writefile_tke_var.units = 'm^2 s^-2'
+writefile_tke_var.description = 'initial profile of turbulence kinetic energy'
+
+writefile_ozone_var = writefile_initial_grp.createVariable('ozone', 'f4', ('levels',))
+writefile_ozone_var[:] = ozone[:,start_t_index]
+writefile_ozone_var.units = 'kg kg^-1'
+writefile_ozone_var.description = 'initial profile of ozone mass mixing ratio'
+
+#forcing group
+
+writefile_p_surf_var = writefile_forcing_grp.createVariable('p_surf', 'f4', ('time',))
+writefile_p_surf_var[:] = p_surf[start_t_index:]
+writefile_p_surf_var.units = 'Pa'
+writefile_p_surf_var.description = 'surface pressure'
+
+writefile_T_surf_var = writefile_forcing_grp.createVariable('T_surf', 'f4', ('time',))
+writefile_T_surf_var[:] = T_surf[start_t_index:]
+writefile_T_surf_var.units = 'K'
+writefile_T_surf_var.description = 'surface absolute temperature'
+
+writefile_w_ls_var = writefile_forcing_grp.createVariable('w_ls', 'f4', ('levels','time',))
+writefile_w_ls_var[:] = w_sub[:,start_t_index:]
+writefile_w_ls_var.units = 'm s^-1'
+writefile_w_ls_var.description = 'large scale vertical velocity'
+
+writefile_omega_var = writefile_forcing_grp.createVariable('omega', 'f4', ('levels','time',))
+writefile_omega_var[:] = omega[:,start_t_index:]
+writefile_omega_var.units = 'Pa s^-1'
+writefile_omega_var.description = 'large scale pressure vertical velocity'
+
+writefile_u_g_var = writefile_forcing_grp.createVariable('u_g', 'f4', ('levels','time',))
+writefile_u_g_var[:] = u_g[:,start_t_index:]
+writefile_u_g_var.units = 'm s^-1'
+writefile_u_g_var.description = 'large scale geostrophic E-W wind'
+
+writefile_v_g_var = writefile_forcing_grp.createVariable('v_g', 'f4', ('levels','time',))
+writefile_v_g_var[:] = v_g[:,start_t_index:]
+writefile_v_g_var.units = 'm s^-1'
+writefile_v_g_var.description = 'large scale geostrophic N-S wind'
+
+writefile_u_nudge_var = writefile_forcing_grp.createVariable('u_nudge', 'f4', ('levels','time',))
+writefile_u_nudge_var[:] = u_wind[:,start_t_index:]
+writefile_u_nudge_var.units = 'm s^-1'
+writefile_u_nudge_var.description = 'E-W wind to nudge toward'
+
+writefile_v_nudge_var = writefile_forcing_grp.createVariable('v_nudge', 'f4', ('levels','time',))
+writefile_v_nudge_var[:] = v_wind[:,start_t_index:]
+writefile_v_nudge_var.units = 'm s^-1'
+writefile_v_nudge_var.description = 'N-S wind to nudge toward'
+
+writefile_T_nudge_var = writefile_forcing_grp.createVariable('T_nudge', 'f4', ('levels','time',))
+writefile_T_nudge_var[:] = T_abs[:,start_t_index:]
+writefile_T_nudge_var.units = 'K'
+writefile_T_nudge_var.description = 'absolute temperature to nudge toward'
+
+writefile_thil_nudge_var = writefile_forcing_grp.createVariable('thil_nudge', 'f4', ('levels','time',))
+writefile_thil_nudge_var[:] = thetailu[:,start_t_index:]
+writefile_thil_nudge_var.units = 'K'
+writefile_thil_nudge_var.description = 'potential temperature to nudge toward'
+
+writefile_qt_nudge_var = writefile_forcing_grp.createVariable('qt_nudge', 'f4', ('levels','time',))
+writefile_qt_nudge_var[:] = qtu[:,start_t_index:]
+writefile_qt_nudge_var.units = 'kg kg^-1'
+writefile_qt_nudge_var.description = 'q_t to nudge toward'
+
+writefile_qi_nudge_var = writefile_forcing_grp.createVariable('qi_nudge', 'f4', ('levels','time',))
+writefile_qi_nudge_var[:] = qi[:,start_t_index:]
+writefile_qi_nudge_var.units = 'kg kg^-1'
+writefile_qi_nudge_var.description = 'q_i to nudge toward'
+
+writefile_ql_nudge_var = writefile_forcing_grp.createVariable('ql_nudge', 'f4', ('levels','time',))
+writefile_ql_nudge_var[:] = ql[:,start_t_index:]
+writefile_ql_nudge_var.units = 'kg kg^-1'
+writefile_ql_nudge_var.description = 'q_l to nudge toward'
+
+writefile_rad_heating_var = writefile_forcing_grp.createVariable('dT_dt_rad', 'f4', ('levels','time',))
+writefile_rad_heating_var[:] = rad_heating[:,start_t_index:]
+writefile_rad_heating_var.units = 'K s^-1'
+writefile_rad_heating_var.description = 'prescribed radiative heating rate'
+
+writefile_h_advec_thil_var = writefile_forcing_grp.createVariable('h_advec_thetail', 'f4', ('levels','time',))
+writefile_h_advec_thil_var[:] = h_advec_thil[:,start_t_index:]
+writefile_h_advec_thil_var.units = 'K s^-1'
+writefile_h_advec_thil_var.description = 'prescribed theta_il tendency due to horizontal advection'
+
+writefile_v_advec_thil_var = writefile_forcing_grp.createVariable('v_advec_thetail', 'f4', ('levels','time',))
+writefile_v_advec_thil_var[:] = v_advec_thil[:,start_t_index:]
+writefile_v_advec_thil_var.units = 'K s^-1'
+writefile_v_advec_thil_var.description = 'prescribed theta_il tendency due to vertical advection'
+
+writefile_h_advec_qt_var = writefile_forcing_grp.createVariable('h_advec_qt', 'f4', ('levels','time',))
+writefile_h_advec_qt_var[:] = h_advec_qt[:,start_t_index:]
+writefile_h_advec_qt_var.units = 'kg kg^-1 s^-1'
+writefile_h_advec_qt_var.description = 'prescribed q_t tendency due to horizontal advection'
+
+writefile_h_advec_qi_var = writefile_forcing_grp.createVariable('h_advec_qi', 'f4', ('levels','time',))
+writefile_h_advec_qi_var[:] = h_advec_qi[:,start_t_index:]
+writefile_h_advec_qi_var.units = 'kg kg^-1 s^-1'
+writefile_h_advec_qi_var.description = 'prescribed q_i tendency due to horizontal advection'
+
+writefile_h_advec_ql_var = writefile_forcing_grp.createVariable('h_advec_ql', 'f4', ('levels','time',))
+writefile_h_advec_ql_var[:] = h_advec_ql[:,start_t_index:]
+writefile_h_advec_ql_var.units = 'kg kg^-1 s^-1'
+writefile_h_advec_ql_var.description = 'prescribed q_l tendency due to horizontal advection'
+
+writefile_v_advec_qt_var = writefile_forcing_grp.createVariable('v_advec_qt', 'f4', ('levels','time',))
+writefile_v_advec_qt_var[:] = v_advec_qt[:,start_t_index:]
+writefile_v_advec_qt_var.units = 'kg kg^-1 s^-1'
+writefile_v_advec_qt_var.description = 'prescribed q_t tendency due to vertical advection'
+
+
+#close processed input file
+writefile_fid.close()
+
+#close raw input file
+nc_fid.close()
diff --git a/scm/etc/scripts/plot_configs/MOSAiC-AMPS.ini b/scm/etc/scripts/plot_configs/MOSAiC-AMPS.ini
new file mode 100644
index 00000000..0ce84096
--- /dev/null
+++ b/scm/etc/scripts/plot_configs/MOSAiC-AMPS.ini
@@ -0,0 +1,63 @@
+scm_datasets = output_MOSAiC-AMPS_SCM_RRFS_v1beta/output.nc, output_MOSAiC-AMPS_SCM_GFS_v17_p8/output.nc, output_MOSAiC-AMPS_SCM_GFS_v16/output.nc
+scm_datasets_labels = RRFS_v1beta, GFSv17p8, GFSv16
+plot_dir = plots_MOSAiC-AMPS_all_suites/
+obs_file = ../data/raw_case_input/MOSAiC_31Oct20190Z_raw.nc
+obs_compare = True
+plot_ind_datasets = True
+time_series_resample = False
+
+[time_slices]
+  [[active]]
+    start = 2019, 11, 1, 0
+    end = 2019, 11, 2, 0
+
+[time_snapshots]
+
+[plots]
+  [[profiles_mean]]
+    vars = ql, qc, qi, qv, T, dT_dt_pbl, T_force_tend, dT_dt_conv, dT_dt_micro, dT_dt_lwrad, dT_dt_phys, qv_force_tend, dq_dt_pbl, dq_dt_shalconv, dq_dt_micro, dq_dt_phys, dq_dt_nonphys, rad_cloud_iwp, rad_cloud_swp, rad_cloud_rwp, rad_cloud_lwp
+    vars_labels = 'cloud water mixing ratio ($g$ $kg^{-1}$)', 'total cloud water mixing ratio ($g$ $kg^{-1}$)', 'ice+snow+graup mixing ratio ($g$ $kg^{-1}$)', 'specific humidity ($g$ $kg^{-1}$)', 'T (K)', 'PBL tendency (K/day)', 'force (K/day)', 'conv. tendency (K/day)', 'microphysics tendency (K/day)', 'LW tendency (K/day)', 'PHYS (K/day)', 'force (kg/kg/day)', 'PBL tendency (kg/kg/day)', 'ShalConv tendency (kg/kg/day)', 'microphysics tendency (kg/kg/day)', 'PHYS (kg/kg/day)', 'NONPHYS (kg/kg/day)', 'IWC (g/m3)', 'SWC (g/m3)', 'RWC (g/m3)', 'LWC (g/m3)'
+    vert_axis = pres_l
+    vert_axis_label = 'average pressure (Pa)'
+    y_inverted = True
+    y_log = False
+#    y_min_option = min             #min, max, val (if val, add y_min = float value)
+    y_min_option = val             #min, max, val (if val, add y_min = float value)
+    y_min = 70000.0             #min, max, val (if val, add y_min = float value)
+    y_max_option = max              #min, max, val (if val, add y_max = float value)
+    conversion_factor = 1000.0, 1000.0, 1000.0, 1000.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0
+
+  [[profiles_mean_multi]]
+    [[[T_forcing]]]
+      vars = T_force_tend, dT_dt_pbl, dT_dt_conv, dT_dt_micro, dT_dt_lwrad
+      vars_labels = 'force', 'PBL', 'Conv', 'MP', 'LW'
+      x_label = 'K/day'
+    [[[QV_forcing]]]
+      vars = qv_force_tend, dq_dt_pbl, dq_dt_shalconv, dq_dt_micro, dq_dt_phys, dq_dt_nonphys
+      vars_labels = 'force', 'PBL', 'ShalConv', 'MP', 'PHYS', 'NONPHYS'
+      x_label = 'kg/kg/day'
+   [[[conv_tendencies]]]
+      vars = dT_dt_deepconv, dT_dt_shalconv
+      vars_labels = 'deep', 'shallow'
+      x_label = 'K/day'
+      
+  [[profiles_instant]]
+
+  [[time_series]]
+    vars = 'pres_s','lhf','shf','tprcp_rate_inst','t2m','q2m','u10m','v10m','gflux','sfc_dwn_lw','tsfc'
+    vars_labels = 'surface pressure (Pa)','latent heat flux ($W$ $m^{-2}$)','sensible heat flux ($W$ $m^{-2}$)','surface rainfall rate ($mm$ $hr{-1}$)','2m temperature ($K$)','2m specific humidity ($g$ $kg{-1}$)','10m zonal wind ($m$ $s{-1}$)','10m meridional wind ($m$ $s{-1}$)','ground flux ($W$ $m{-2}$)','downward longwave flux ($W$ $m{-2}$)','surface temperature ($K$)'
+
+  [[contours]]
+    vars = qv, T, h_advec_qt, h_advec_thil, qi, ql
+    vars_labels = 'Water Vapor ($g$ $kg^{-1}$)','Temperature ($K$)','Horizontal total water advection ($g$ $kg^{-1}$ $s^{-1}$)','Horizontal thetail advection ($K$ $s^{-1}$)','QI ($g$ $kg^{-1}$)','QL ($g$ $kg^{-1}$)'
+    vert_axis = pres_l
+    vert_axis_label = 'p (Pa)'
+    y_inverted = True
+    y_log = False
+    y_min_option = val             #min, max, val (if val, add y_min = float value)
+    y_min = 10000.0
+    y_max_option = val              #min, max, val (if val, add y_max = float value)
+    y_max = 100000.0
+    x_ticks_num = 10
+    y_ticks_num = 10
+    conversion_factor = 1000.0, 1.0, 1.0, 1.0, 1000.0, 1000.0
diff --git a/scm/etc/scripts/plot_configs/MOSAiC-SS.ini b/scm/etc/scripts/plot_configs/MOSAiC-SS.ini
new file mode 100644
index 00000000..82672321
--- /dev/null
+++ b/scm/etc/scripts/plot_configs/MOSAiC-SS.ini
@@ -0,0 +1,65 @@
+scm_datasets = output_MOSAiC-SS_SCM_RRFS_v1beta/output.nc, output_MOSAiC-SS_SCM_GFS_v17_p8/output.nc, output_MOSAiC-SS_SCM_GFS_v16/output.nc
+##scm_datasets = output_MOSAiC-SS_SCM_GFS_v17_p8/output.nc, output_MOSAiC-SS_SCM_GFS_v16/output.nc
+##scm_datasets_labels = GFSv17p8, GFSv16
+scm_datasets_labels = RRFS_v1beta, GFSv17p8, GFSv16
+plot_dir = plots_MOSAiC-SS_all_suites/
+obs_file = ../data/raw_case_input/MOSAiC_2Mar20200Z_raw.nc
+obs_compare = True
+plot_ind_datasets = True
+time_series_resample = False
+
+[time_slices]
+  [[active]]
+    start = 2020, 3, 4, 0
+    end = 2020, 3, 5, 0
+
+[time_snapshots]
+
+[plots]
+  [[profiles_mean]]
+    vars = ql, qc, qi, qv, T, dT_dt_pbl, T_force_tend, dT_dt_conv, dT_dt_micro, dT_dt_lwrad, dT_dt_phys, qv_force_tend, dq_dt_pbl, dq_dt_shalconv, dq_dt_micro, dq_dt_phys, dq_dt_nonphys, rad_cloud_iwp, rad_cloud_swp, rad_cloud_rwp, rad_cloud_lwp
+    vars_labels = 'cloud water mixing ratio ($g$ $kg^{-1}$)', 'total cloud water mixing ratio ($g$ $kg^{-1}$)', 'ice+snow+graup mixing ratio ($g$ $kg^{-1}$)', 'specific humidity ($g$ $kg^{-1}$)', 'T (K)', 'PBL tendency (K/day)', 'force (K/day)', 'conv. tendency (K/day)', 'microphysics tendency (K/day)', 'LW tendency (K/day)', 'PHYS (K/day)', 'force (kg/kg/day)', 'PBL tendency (kg/kg/day)', 'ShalConv tendency (kg/kg/day)', 'microphysics tendency (kg/kg/day)', 'PHYS (kg/kg/day)', 'NONPHYS (kg/kg/day)', 'IWC (g/m3)', 'SWC (g/m3)', 'RWC (g/m3)', 'LWC (g/m3)'
+    vert_axis = pres_l
+    vert_axis_label = 'average pressure (Pa)'
+    y_inverted = True
+    y_log = False
+#    y_min_option = min             #min, max, val (if val, add y_min = float value)
+    y_min_option = val             #min, max, val (if val, add y_min = float value)
+    y_min = 70000.0             #min, max, val (if val, add y_min = float value)
+    y_max_option = max              #min, max, val (if val, add y_max = float value)
+    conversion_factor = 1000.0, 1000.0, 1000.0, 1000.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0
+
+  [[profiles_mean_multi]]
+    [[[T_forcing]]]
+      vars = T_force_tend, dT_dt_pbl, dT_dt_conv, dT_dt_micro, dT_dt_lwrad
+      vars_labels = 'force', 'PBL', 'Conv', 'MP', 'LW'
+      x_label = 'K/day'
+    [[[QV_forcing]]]
+      vars = qv_force_tend, dq_dt_pbl, dq_dt_shalconv, dq_dt_micro, dq_dt_phys, dq_dt_nonphys
+      vars_labels = 'force', 'PBL', 'ShalConv', 'MP', 'PHYS', 'NONPHYS'
+      x_label = 'kg/kg/day'
+   [[[conv_tendencies]]]
+      vars = dT_dt_deepconv, dT_dt_shalconv
+      vars_labels = 'deep', 'shallow'
+      x_label = 'K/day'
+      
+  [[profiles_instant]]
+
+  [[time_series]]
+    vars = 'pres_s','lhf','shf','tprcp_rate_inst','t2m','q2m','u10m','v10m','gflux','sfc_dwn_lw','tsfc'
+    vars_labels = 'surface pressure (Pa)','latent heat flux ($W$ $m^{-2}$)','sensible heat flux ($W$ $m^{-2}$)','surface rainfall rate ($mm$ $hr{-1}$)','2m temperature ($K$)','2m specific humidity ($g$ $kg{-1}$)','10m zonal wind ($m$ $s{-1}$)','10m meridional wind ($m$ $s{-1}$)','ground flux ($W$ $m{-2}$)','downward longwave flux ($W$ $m{-2}$)','surface temperature ($K$)'
+
+  [[contours]]
+    vars = qv, T, h_advec_qt, h_advec_thil, qi, ql
+    vars_labels = 'Water Vapor ($g$ $kg^{-1}$)','Temperature ($K$)','Horizontal total water advection ($g$ $kg^{-1}$ $s^{-1}$)','Horizontal thetail advection ($K$ $s^{-1}$)','QI ($g$ $kg^{-1}$)','QL ($g$ $kg^{-1}$)'
+    vert_axis = pres_l
+    vert_axis_label = 'p (Pa)'
+    y_inverted = True
+    y_log = False
+    y_min_option = val             #min, max, val (if val, add y_min = float value)
+    y_min = 10000.0
+    y_max_option = val              #min, max, val (if val, add y_max = float value)
+    y_max = 100000.0
+    x_ticks_num = 10
+    y_ticks_num = 10
+    conversion_factor = 1000.0, 1.0, 1.0, 1.0, 1000.0, 1000.0
diff --git a/scm/etc/scripts/scm_analysis.py b/scm/etc/scripts/scm_analysis.py
index 34c85c52..d6516dc7 100755
--- a/scm/etc/scripts/scm_analysis.py
+++ b/scm/etc/scripts/scm_analysis.py
@@ -856,6 +856,8 @@ def replace_fill_with_nan(nc_ds, var_name, var, group, time_diag, pres_l, datase
         obs_dict = sro.read_LASSO_obs(obs_file, time_slices, date_inst)
     elif('gabls3' in case_name.strip()):
         obs_dict = sro.read_gabls3_obs(obs_file, time_slices, date_inst)
+    elif('MOSAiC' in case_name.strip()):
+        obs_dict = sro.read_MOSAiC_obs(obs_file, time_slices, date_inst)
 
 try:
     os.makedirs(plot_dir)
@@ -1378,7 +1380,7 @@ def replace_fill_with_nan(nc_ds, var_name, var, group, time_diag, pres_l, datase
             else:
                 y_max_val = profiles_mean['y_max']
             y_lim_val = [y_min_val, y_max_val]
-
+            
             #plot mean profiles
             for k in range(len(profiles_mean['vars'])):
                 #get the python variable associated with the vars listed in the config file
diff --git a/scm/etc/scripts/scm_read_obs.py b/scm/etc/scripts/scm_read_obs.py
index 59ac43c3..ac6b8010 100644
--- a/scm/etc/scripts/scm_read_obs.py
+++ b/scm/etc/scripts/scm_read_obs.py
@@ -6,6 +6,71 @@
 import math
 import forcing_file_common as ffc
 
+def read_MOSAiC_obs(obs_file, time_slices, date):
+  obs_time_slice_indices = []
+
+  obs_fid = Dataset(obs_file, 'r')
+
+  obs_year = obs_fid.variables['year'][:]
+  obs_month = obs_fid.variables['month'][:]
+  obs_day = obs_fid.variables['day'][:]
+  obs_hour = obs_fid.variables['hour'][:]
+  obs_time = obs_fid.variables['time_offset'][:]
+
+  obs_date = []
+  for i in range(obs_hour.size):
+      obs_date.append(datetime.datetime(obs_year[i], obs_month[i], obs_day[i], obs_hour[i], 0, 0, 0))
+  obs_date = np.array(obs_date)
+
+  for time_slice in time_slices:
+      start_date = datetime.datetime(time_slices[time_slice]['start'][0], time_slices[time_slice]['start'][1],time_slices[time_slice]['start'][2], time_slices[time_slice]['start'][3], time_slices[time_slice]['start'][4])
+      end_date = datetime.datetime(time_slices[time_slice]['end'][0], time_slices[time_slice]['end'][1],time_slices[time_slice]['end'][2], time_slices[time_slice]['end'][3], time_slices[time_slice]['end'][4])
+      start_date_index = np.where(obs_date == start_date)[0][0]
+      end_date_index = np.where(obs_date == end_date)[0][0]
+      obs_time_slice_indices.append([start_date_index, end_date_index])
+
+#print(start_date, end_date, start_date_index, end_date_index, obs_date[start_date_index], obs_date[end_date_index])
+
+  #find the index corresponding to the start of the simulations
+  obs_start_index = np.where(obs_date == date[0][0])[0]
+  obs_time = obs_time - obs_time[obs_start_index]
+
+  obs_pres_l = obs_fid.variables['levels'][:]*100.0 #pressure levels in mb
+
+  obs_T = obs_fid.variables['T'][:]
+  obs_q = obs_fid.variables['q'][:]
+  obs_qi= obs_fid.variables['qi'][:]
+  obs_ql= obs_fid.variables['ql'][:]
+  obs_u = obs_fid.variables['u'][:]
+  obs_v = obs_fid.variables['v'][:]
+  obs_rad_net_srf = obs_fid.variables['rad_net_srf'][:]
+  obs_lw_dn_srf = obs_fid.variables['lw_dn_srf'][:]
+  obs_tsk = obs_fid.variables['T_skin'][:]
+  obs_shf = obs_fid.variables['SH'][:]  
+  obs_lhf = obs_fid.variables['LH'][:]
+  obs_t2m = obs_fid.variables['T_srf'][:]
+  obs_q2m = obs_fid.variables['q_srf'][:]
+
+  obs_time_h = obs_time/3600.0
+
+  Rd = 287.0
+  Rv = 461.0
+
+  e_s = 6.1078*np.exp(17.2693882*(obs_T - 273.16)/(obs_T - 35.86))*100.0 #Tetens formula produces e_s in mb (convert to Pa)
+  e = obs_q*obs_pres_l/(obs_q + (Rd/Rv)*(1.0 - obs_q)) #compute vapor pressure from specific humidity
+  obs_rh = np.clip(e/e_s, 0.0, 1.0)
+
+  return_dict = {'year': obs_year, 'month': obs_month, 'day': obs_day, 'hour': obs_hour,
+    'time': obs_time, 'date': obs_date, 'time_slice_indices': obs_time_slice_indices,
+    'pres_l': obs_pres_l, 'T': obs_T, 'q': obs_q, 'rh': obs_rh, 'u': obs_u, 'v': obs_v, 'shf': obs_shf,
+    'lhf': obs_lhf, 't2m': obs_t2m, 'q2m': obs_q2m, 'time_h': obs_time_h, 'tsfc': obs_tsk,
+    'qv': obs_q, 'qi': obs_qi, 'ql': obs_ql, 'rad_net_srf': obs_rad_net_srf, 'sfc_dwn_lw': obs_lw_dn_srf}
+#    'lwp': obs_lwp, 'T_force_tend': obs_T_forcing, 'qv_force_tend': obs_q_forcing}
+
+  obs_fid.close()
+
+  return return_dict
+
 def read_twpice_obs(obs_file, time_slices, date):
   obs_time_slice_indices = []
 
@@ -337,4 +402,4 @@ def read_gabls3_obs(obs_file, time_slices, date):
         'sfc_dwn_sw': obs_sw_dn, 'sfc_up_sw': obs_sw_up, 'sfc_rad_net_land': obs_sfc_rad_net, 'gflux': -1*obs_gflux, 't2m':obs_t2m, 'q2m':obs_q2m, 
         'ustar':obs_ustar,'u10m':obs_u10m, 'v10m':obs_v10m, 'hpbl':obs_hpbl, 'tsfc':obs_tsk}
     
-    return return_dict
\ No newline at end of file
+    return return_dict