diff --git a/exp/test_cases/thermodynamic_sea_ice/grey_thermo_ice_test_case.py b/exp/test_cases/thermodynamic_sea_ice/grey_thermo_ice_test_case.py new file mode 100644 index 000000000..0e331ad48 --- /dev/null +++ b/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) + + + diff --git a/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 b/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 index 9849c30c0..efea06f25 100644 --- a/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 +++ b/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,6 +453,11 @@ 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) @@ -453,6 +465,9 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, ! 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) diff --git a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 index 81963dbc6..73764ea4d 100644 --- a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 +++ b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 @@ -174,6 +174,7 @@ module idealized_moist_phys_mod real, allocatable, dimension(:,:) :: & z_surf, & ! surface height t_surf, & ! surface temperature + t_ml, h_thermo_ice, const_correct, nudge_out, & ! NTL 01/23 thermodynamic sea ice q_surf, & ! surface moisture u_surf, & ! surface U wind v_surf, & ! surface V wind @@ -392,7 +393,7 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l r_conv_scheme = NO_CONV lwet_convection = .false. do_bm = .false. - do_ras = .false. + do_ras = .false. call error_mesg('idealized_moist_phys','No convective adjustment scheme used.', NOTE) else if(uppercase(trim(convection_scheme)) == 'SIMPLE_BETTS_MILLER') then @@ -421,7 +422,7 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l call error_mesg('idealized_moist_phys','Using dry convection scheme.', NOTE) lwet_convection = .false. do_bm = .false. - do_ras = .false. + do_ras = .false. else if(uppercase(trim(convection_scheme)) == 'UNSET') then call error_mesg('idealized_moist_phys','determining convection scheme from flags', NOTE) @@ -432,7 +433,7 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l if (do_bm) then r_conv_scheme = FULL_BETTS_MILLER_CONV call error_mesg('idealized_moist_phys','Using Betts-Miller convection scheme.', NOTE) - end if + end if if (do_ras) then r_conv_scheme = RAS_CONV call error_mesg('idealized_moist_phys','Using relaxed Arakawa Schubert convection scheme.', NOTE) @@ -475,6 +476,12 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l allocate(depth_change_conv(is:ie, js:je)) allocate(z_surf (is:ie, js:je)) allocate(t_surf (is:ie, js:je)) +! NTL 01/23 thermodynamic sea ice +allocate(h_thermo_ice(is:ie, js:je)) +allocate(t_ml (is:ie, js:je)) +allocate(const_correct(is:ie, js:je)) +allocate(nudge_out(is:ie, js:je)) +! allocate(q_surf (is:ie, js:je)); q_surf = 0.0 allocate(u_surf (is:ie, js:je)); u_surf = 0.0 allocate(v_surf (is:ie, js:je)); v_surf = 0.0 @@ -642,7 +649,10 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l ! to quickly enter the atmosphere avoiding problems with the convection scheme t_surf = t_surf_init + 1.0 - call mixed_layer_init(is, ie, js, je, num_levels, t_surf, bucket_depth, get_axis_id(), Time, albedo, rad_lonb_2d(:,:), rad_latb_2d(:,:), land, bucket) ! t_surf is intent(inout) ! albedo distribution set here. + call mixed_layer_init(is, ie, js, je, num_levels, t_surf, & + h_thermo_ice, t_ml, const_correct,nudge_out, & ! NTL 01/23 thermodynamic sea ice + bucket_depth, get_axis_id(), Time, albedo, rad_lonb_2d(:,:), rad_latb_2d(:,:), land, bucket) ! t_surf is intent(inout) !s albedo distribution set here. + elseif(gp_surface) then albedo=0.0 @@ -1058,7 +1068,7 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, ps tg(:,:,:,previous), & net_surf_sw_down(:,:), & surf_lw_down(:,:), albedo, & - grid_tracers(:,:,:,previous,nsphum)) + grid_tracers(:,:,:,previous,nsphum), const_correct(:,:), nudge_out(:,:)) end if if(.not.mixed_layer_bc) then @@ -1311,6 +1321,7 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, ps js, & je, & t_surf(:,:), & ! t_surf is intent(inout) + h_thermo_ice(:,:), t_ml(:,:), const_correct(:,:), nudge_out(:,:), & ! NTL 01/23 thermodynamic sea ice flux_t(:,:), & flux_q(:,:), & flux_r(:,:), & @@ -1336,9 +1347,12 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, ps endif ! if(turb) then -! Adding relative humidity calculation so as to allow comparison with Frierson's thesis. - call rh_calc (p_full(:,:,:,previous),tg_tmp,qg_tmp,RH) - if(id_rh >0) used = send_data(id_rh, RH*100., Time) +!s Adding relative humidity calculation so as to allow comparison with Frierson's thesis. + + if(id_rh >0) then + call rh_calc (p_full(:,:,:,previous),tg_tmp,qg_tmp,RH) + used = send_data(id_rh, RH*100., Time) + endif ! Add bucket ! Timestepping for bucket. @@ -1406,7 +1420,9 @@ subroutine idealized_moist_phys_end call vert_turb_driver_end endif call lscale_cond_end -if(mixed_layer_bc) call mixed_layer_end(t_surf, bucket_depth, bucket) +if(mixed_layer_bc) call mixed_layer_end(t_surf, & + h_thermo_ice, t_ml, const_correct, nudge_out, albedo, & ! NTL 01/23 thermodynamic sea ice + bucket_depth, bucket) if(do_damping) call damping_driver_end #ifdef SOC_NO_COMPILE diff --git a/src/atmos_spectral/driver/solo/mixed_layer.F90 b/src/atmos_spectral/driver/solo/mixed_layer.F90 index 051d985fe..e9ad6c90d 100644 --- a/src/atmos_spectral/driver/solo/mixed_layer.F90 +++ b/src/atmos_spectral/driver/solo/mixed_layer.F90 @@ -37,18 +37,19 @@ module mixed_layer_mod use tracer_manager_mod, only: get_tracer_names, get_number_tracers -use constants_mod, only: HLV, PI, RHO_CP, CP_AIR, KELVIN +use constants_mod, only: HLV, PI, RHO_CP, CP_AIR, KELVIN, & + HLF, DENS_ICE, TFREEZE ! NTL 01/23 for sea ice code use diag_manager_mod, only: register_diag_field, register_static_field, send_data use time_manager_mod, only: time_type #ifdef COLUMN_MODEL -use column_grid_mod, only: get_deg_lon, get_deg_lat +use column_grid_mod, only: get_deg_lon, get_deg_lat, area_weighted_global_mean use column_mod, only: get_surf_geopotential use spec_mpp_mod, only: grid_domain #else -use transforms_mod, only: get_deg_lat, get_deg_lon, grid_domain +use transforms_mod, only: get_deg_lat, get_deg_lon, grid_domain, area_weighted_global_mean ! mj know about surface topography use spectral_dynamics_mod,only: get_surf_geopotential #endif @@ -136,6 +137,28 @@ module mixed_layer_mod character(len=256) :: flux_lhe_anom_file_name = 'INPUT/flux_lhe_anom.nc' character(len=256) :: flux_lhe_anom_field_name = 'flux_lhe_anom' +! THERMO ICE OPTIONS - NTL 01/23 +logical :: do_explicit = .false. ! explicit or implicit timestepping for ocean covered surface +logical :: do_thermo_ice = .false. ! do thermodynamic sea ice? +real :: thermo_ice_albedo = 0.4 ! from Feldl & Merlis - original XZ is 0.8, M. England 0.6 +real :: L_thermo_ice = DENS_ICE * HLF ! latent heat of fusion (J/m^3) +real, parameter :: k_thermo_ice = 2 ! conductivity of ice (W/m/K) +real, parameter :: thermo_ice_basal_flux_const = 120 ! linear coefficient for heat flux from ocean to ice base (W/m^2/K) +real :: t_thermo_ice_base = TFREEZE ! temperature at base of ice, taken to be freshwater freezing point +real :: t_surf_freeze = TFREEZE +real :: thermo_ice_nudge_timescale = 7. * 86400. +logical :: do_nudge_ice = .false. +logical :: do_nudging_correction = .false. +character(len=256) :: nudge_ice_file_name = 'nudge_ice' +logical :: time_varying_nudge_ice = .false. +logical :: do_prescribe_albedo = .false. +character(len=256) :: albedo_file_name = 'albedo' +logical :: time_varying_albedo = .false. +logical :: prescribe_nh_albedo_only = .false. +logical :: read_const_correct = .true. ! set to false if doing nudging, and want to restart from simulation with no nudging +logical :: read_nudge_out = .true. +logical :: albedo_from_nudge_file = .false. + namelist/mixed_layer_nml/ evaporation, depth, qflux_amp, qflux_width, tconst,& delta_T, prescribe_initial_dist,albedo_value, & land_depth,trop_depth, & !mj @@ -153,7 +176,11 @@ module mixed_layer_mod ice_albedo_value, specify_sst_over_ocean_only, & ice_concentration_threshold, ice_albedo_method,& add_latent_heat_flux_anom,flux_lhe_anom_file_name,& - flux_lhe_anom_field_name, do_ape_sst, qflux_field_name + flux_lhe_anom_field_name, do_ape_sst, qflux_field_name,& + do_explicit, do_thermo_ice, thermo_ice_albedo, t_surf_freeze, t_thermo_ice_base, & ! NTL ice + do_nudge_ice, nudge_ice_file_name, thermo_ice_nudge_timescale, time_varying_nudge_ice, & + do_nudging_correction, do_prescribe_albedo, time_varying_albedo, albedo_file_name, & + prescribe_nh_albedo_only, read_const_correct, read_nudge_out, albedo_from_nudge_file !================================================================================================================================= @@ -171,7 +198,11 @@ module mixed_layer_mod id_heat_cap, & ! heat capacity id_albedo, & ! mj albedo id_ice_conc, & ! st ice concentration - id_delta_t_surf + id_delta_t_surf, & + ! NTL 01/23 thermodynamic sea ice + id_h_thermo_ice, & ! sea ice thickness + id_t_ml, & ! mixed layer temperature + id_flux_thermo_ice ! conductive heat flux through ice real, allocatable, dimension(:,:) :: & ocean_qflux, & ! Q-flux @@ -201,10 +232,23 @@ module mixed_layer_mod zsurf, & ! mj know about topography land_sea_heat_capacity,& sst_new, & ! mj input SST - albedo_initial + albedo_initial, & + ! NTL 01/23 thermodynamic sea ice + dFdt_surf, & ! d(corrected_flux)/d(t_surf), for calculation of ice sfc temperature + delta_t_ml, & ! Increment in mixed layer temperature + delta_h_thermo_ice, & ! Increment in sea ice thickness + delta_t_thermo_ice, & ! Increment in ice (steady-state) surface temperature + flux_thermo_ice, & ! Conductive heat flux through ice + corrected_flux_noq +real, allocatable, dimension(:,:) :: h_thermo_ice_in, nudge_flux, nudge_correction, albedo_in +real :: const_nudge_correction +type(interpolate_type),save :: nudge_ice_interp +type(interpolate_type),save :: albedo_interp logical, allocatable, dimension(:,:) :: land_mask + + !mj read sst from input file type(interpolate_type),save :: sst_interp type(interpolate_type),save :: qflux_interp @@ -217,10 +261,13 @@ module mixed_layer_mod contains !================================================================================================================================= -subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, bucket_depth, axes, Time, albedo, rad_lonb_2d,rad_latb_2d, land, restart_file_bucket_depth) +subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, & + h_thermo_ice, t_ml, const_correct, nudge_out, & ! NTL 01/23 thermodynamic sea ice + bucket_depth, axes, Time, albedo, rad_lonb_2d,rad_latb_2d, land, restart_file_bucket_depth) type(time_type), intent(in) :: Time -real, intent(out), dimension(:,:) :: t_surf, albedo +real, intent(out), dimension(:,:) :: t_surf, albedo +real, intent(out), dimension(:,:) :: h_thermo_ice, t_ml, const_correct, nudge_out ! NTL 01/23 thermodynamic sea ice real, intent(out), dimension(:,:,:) :: bucket_depth integer, intent(in), dimension(4) :: axes real, intent(in), dimension(:,:) :: rad_lonb_2d, rad_latb_2d @@ -284,12 +331,24 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, bucket_depth, ax allocate(delta_t_surf (is:ie, js:je)) allocate(eff_heat_capacity (is:ie, js:je)) allocate(corrected_flux (is:ie, js:je)) +allocate(corrected_flux_noq (is:ie, js:je)) ! NTL 01/23 thermodynamic sea ice allocate(t_surf_dependence (is:ie, js:je)) allocate (albedo_initial (is:ie, js:je)) allocate(land_sea_heat_capacity (is:ie, js:je)) allocate(zsurf (is:ie, js:je)) allocate(sst_new (is:ie, js:je)) allocate(land_mask (is:ie, js:je)); land_mask=land + +! NTL 01/23 thermodynamic sea ice +allocate(dFdt_surf (is:ie, js:je)) +allocate(delta_t_ml (is:ie, js:je)) +allocate(delta_h_thermo_ice (is:ie, js:je)) +allocate(delta_t_thermo_ice (is:ie, js:je)) +allocate(flux_thermo_ice (is:ie, js:je)) +allocate(h_thermo_ice_in (is:ie, js:je)) +allocate(albedo_in (is:ie, js:je)) +allocate(nudge_flux (is:ie, js:je)) +allocate(nudge_correction (is:ie, js:je)) ! !see if restart file exists for the surface temperature ! @@ -298,6 +357,12 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, bucket_depth, ax ! latitude will be needed for oceanic q flux !s Moved up slightly so that rad_lat_2d can be used in initial temperature distribution if necessary. +if (do_thermo_ice .and. (.not. do_explicit)) then + call error_mesg('mixed_layer','setting do_explicit .true. for consistency with sea ice', WARNING) + do_explicit = .true. +endif + + call get_deg_lat(deg_lat) do j=js,je rad_lat_2d(:,j) = deg_lat(j)*PI/180. @@ -325,6 +390,23 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, bucket_depth, ax call nullify_domain() call read_data(trim('INPUT/mixed_layer.res'), 't_surf', t_surf, grid_domain) + if (do_thermo_ice) then + call read_data(trim('INPUT/mixed_layer.res'), 'h_thermo_ice', h_thermo_ice, grid_domain) + call read_data(trim('INPUT/mixed_layer.res'), 't_ml', t_ml, grid_domain) + if (read_const_correct) then + call read_data(trim('INPUT/mixed_layer.res'), 'const_correct', const_correct, grid_domain) + else + const_correct(:,:) = 0.0 + endif + if (read_nudge_out) then + call read_data(trim('INPUT/mixed_layer.res'), 'nudge_out', nudge_out, grid_domain) + else + nudge_out(:,:) = 0.0 + endif + call read_data(trim('INPUT/mixed_layer.res'), 'albedo', albedo, grid_domain) + else if (do_prescribe_albedo) then + call read_data(trim('INPUT/mixed_layer.res'), 'albedo', albedo, grid_domain) + endif if (restart_file_bucket_depth) then call read_data(trim('INPUT/mixed_layer.res'), 'bucket_depth', bucket_depth, grid_domain) @@ -344,14 +426,28 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, bucket_depth, ax elseif (prescribe_initial_dist) then ! call error_mesg('mixed_layer','mixed_layer restart file not found - initializing from prescribed distribution', WARNING) + t_surf(:,:) = tconst - delta_T*((3.*sin(rad_lat_2d)**2.)-1.)/3. - + if (do_thermo_ice) then ! NTL 01/23 thermodynamic sea ice + h_thermo_ice(:,:) = 0.0 + t_ml(:,:) = t_surf(:,:) + const_correct(:,:) = 0.0 + nudge_out(:,:) = 0.0 + albedo(:,:) = albedo_value + where(land) albedo(:,:) = land_albedo_prefactor * albedo_value + albedo_initial(:,:) = albedo (:,:) + else if (do_prescribe_albedo) then + albedo(:,:) = albedo_value + where(land) albedo(:,:) = land_albedo_prefactor * albedo_value + albedo_initial(:,:) = albedo (:,:) + endif else call error_mesg('mixed_layer','mixed_layer restart file not found - initializing from lowest model level temp', WARNING) endif + if(trim(ice_albedo_method) == 'ramp_function') then call error_mesg('mixed_layer','Alternative method ramp_function used for ice albedo output.', NOTE) endif @@ -368,11 +464,21 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, bucket_depth, ax axes(1:2), 'mixed layer heat capacity','joules/m^2/deg C') id_delta_t_surf = register_diag_field(mod_name, 'delta_t_surf', & axes(1:2), Time, 'change in sst','K') -if (update_albedo_from_ice) then +if (update_albedo_from_ice .or. do_thermo_ice .or. do_prescribe_albedo) then id_albedo = register_diag_field(mod_name, 'albedo', & axes(1:2), Time, 'surface albedo', 'none') - id_ice_conc = register_diag_field(mod_name, 'ice_conc', & - axes(1:2), Time, 'ice_concentration', 'none') + if (do_thermo_ice) then + ! NTL 01/23 thermodynamic sea ice + id_h_thermo_ice = register_diag_field(mod_name, 'h_thermo_ice', & + axes(1:2), Time, 'sea ice thickness','m') + id_t_ml = register_diag_field(mod_name, 't_ml', & + axes(1:2), Time, 'mixed layer tempeature','K') + id_flux_thermo_ice = register_diag_field(mod_name, 'flux_thermo_ice', & + axes(1:2), Time, 'conductive heat flux through sea ice','watts/m2') + else + id_ice_conc = register_diag_field(mod_name, 'ice_conc', & + axes(1:2), Time, 'ice_concentration', 'none') + endif else id_albedo = register_static_field(mod_name, 'albedo', & axes(1:2), 'surface albedo', 'none') @@ -408,6 +514,64 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, bucket_depth, ax endif endif + +h_thermo_ice_in = 999. ! this is used to determine whether nudging should occur + ! setting a large number ensures this will do nothing if not initialised from external file +! load ice +if (do_nudge_ice .or. albedo_from_nudge_file) then + if (time_varying_nudge_ice .or. albedo_from_nudge_file) then + call interpolator_init( nudge_ice_interp, trim(nudge_ice_file_name)//'.nc', rad_lonb_2d, rad_latb_2d, data_out_of_bounds=(/CONSTANT/) ) + else + + if(file_exist(trim('INPUT/'//nudge_ice_file_name//'.nc'))) then + call mpp_get_global_domain(grid_domain, xsize=global_num_lon, ysize=global_num_lat) + call field_size(trim('INPUT/'//nudge_ice_file_name//'.nc'), trim(nudge_ice_file_name), siz) + if ( siz(1) == global_num_lon .or. siz(2) == global_num_lat ) then + call read_data(trim('INPUT/'//nudge_ice_file_name//'.nc'), trim(nudge_ice_file_name), h_thermo_ice_in, grid_domain) + else + write(ctmp1(1: 4),'(i4)') siz(1) + write(ctmp1(9:12),'(i4)') siz(2) + write(ctmp2(1: 4),'(i4)') global_num_lon + write(ctmp2(9:12),'(i4)') global_num_lat + call error_mesg ('get_ice','nudge_ice file contains data on a '// & + ctmp1//' grid, but atmos model grid is '//ctmp2, FATAL) + endif + else + call error_mesg('get_ice','do_nudge_ice="'//trim('True')//'"'// & + ' but '//trim(nudge_ice_file_name)//' does not exist', FATAL) + endif + + endif +endif + +albedo_in = 0.0 +if (do_prescribe_albedo) then + if (time_varying_albedo) then + call interpolator_init( albedo_interp, trim(albedo_file_name)//'.nc', rad_lonb_2d, rad_latb_2d, data_out_of_bounds=(/CONSTANT/) ) + else + + if(file_exist(trim('INPUT/'//albedo_file_name//'.nc'))) then + call mpp_get_global_domain(grid_domain, xsize=global_num_lon, ysize=global_num_lat) + call field_size(trim('INPUT/'//albedo_file_name//'.nc'), trim(albedo_file_name), siz) + if ( siz(1) == global_num_lon .or. siz(2) == global_num_lat ) then + call read_data(trim('INPUT/'//albedo_file_name//'.nc'), trim(albedo_file_name), albedo_in, grid_domain) + else + write(ctmp1(1: 4),'(i4)') siz(1) + write(ctmp1(9:12),'(i4)') siz(2) + write(ctmp2(1: 4),'(i4)') global_num_lon + write(ctmp2(9:12),'(i4)') global_num_lat + call error_mesg ('get_albedo','albedo file contains data on a '// & + ctmp1//' grid, but atmos model grid is '//ctmp2, FATAL) + endif + else + call error_mesg('get_albedo','do_prescribe_albedo="'//trim('True')//'"'// & + ' but '//trim(albedo_file_name)//' does not exist', FATAL) + endif + + endif +endif + + !s Adding MiMA options for qfluxes. if ( do_qflux .or. do_warmpool) then @@ -430,65 +594,70 @@ subroutine mixed_layer_init(is, ie, js, je, num_levels, t_surf, bucket_depth, ax !s Prescribe albedo distribution here so that it will be the same in both two_stream_gray later and rrtmg radiation. -albedo(:,:) = albedo_value +if ((.not. do_thermo_ice) .and. (.not. do_prescribe_albedo)) then ! NTL 01/23 thermodynamic sea ice -- this is set / loaded earlier -if(trim(land_option) .eq. 'input') then + albedo(:,:) = albedo_value -where(land) albedo = land_albedo_prefactor * albedo + if(trim(land_option) .eq. 'input') then -endif + where(land) albedo = land_albedo_prefactor * albedo -!mj MiMA albedo choices. -select case (albedo_choice) - case (2) ! higher_albedo northward (lat_glacier>0) or southward (lat_glacier <0 ) of lat_glacier - do j = 1, size(t_surf,2) - lat = deg_lat(js+j-1) - ! mj SH or NH only - if ( lat_glacier .ge. 0. ) then - if ( lat > lat_glacier ) then - albedo(:,j) = higher_albedo - else - albedo(:,j) = albedo_value - endif - else - if ( lat < lat_glacier ) then + endif + + !mj MiMA albedo choices. + select case (albedo_choice) + case (2) ! higher_albedo northward (lat_glacier>0) or southward (lat_glacier <0 ) of lat_glacier + do j = 1, size(t_surf,2) + lat = deg_lat(js+j-1) + ! mj SH or NH only + if ( lat_glacier .ge. 0. ) then + if ( lat > lat_glacier ) then + albedo(:,j) = higher_albedo + else + albedo(:,j) = albedo_value + endif + else + if ( lat < lat_glacier ) then + albedo(:,j) = higher_albedo + else + albedo(:,j) = albedo_value + endif + endif + enddo + case (3) ! higher_albedo poleward of lat_glacier + do j = 1, size(t_surf,2) + lat = deg_lat(js+j-1) + if ( abs(lat) > lat_glacier ) then albedo(:,j) = higher_albedo - else + else albedo(:,j) = albedo_value - endif - endif - enddo - case (3) ! higher_albedo poleward of lat_glacier - do j = 1, size(t_surf,2) - lat = deg_lat(js+j-1) - if ( abs(lat) > lat_glacier ) then - albedo(:,j) = higher_albedo - else - albedo(:,j) = albedo_value - endif - enddo - case (4) ! exponential increase with albedo_exp - do j = 1, size(t_surf,2) - lat = abs(deg_lat(js+j-1)) - albedo(:,j) = albedo_value + (higher_albedo-albedo_value)*(lat/90.)**albedo_exp - enddo - case (5) ! tanh increase around albedo_cntr with albedo_wdth - do j = 1, size(t_surf,2) - lat = abs(deg_lat(js+j-1)) - albedo(:,j) = albedo_value + (higher_albedo-albedo_value)*& - 0.5*(1+tanh((lat-albedo_cntr)/albedo_wdth)) - enddo -end select - -albedo_initial=albedo - -if (update_albedo_from_ice) then - call interpolator_init( ice_interp, trim(ice_file_name)//'.nc', rad_lonb_2d, rad_latb_2d, data_out_of_bounds=(/CONSTANT/) ) - call read_ice_conc(Time) - call albedo_calc(albedo,Time) -else - if ( id_albedo > 0 ) used = send_data ( id_albedo, albedo ) -endif + endif + enddo + case (4) ! exponential increase with albedo_exp + do j = 1, size(t_surf,2) + lat = abs(deg_lat(js+j-1)) + albedo(:,j) = albedo_value + (higher_albedo-albedo_value)*(lat/90.)**albedo_exp + enddo + case (5) ! tanh increase around albedo_cntr with albedo_wdth + do j = 1, size(t_surf,2) + lat = abs(deg_lat(js+j-1)) + albedo(:,j) = albedo_value + (higher_albedo-albedo_value)*& + 0.5*(1+tanh((lat-albedo_cntr)/albedo_wdth)) + enddo + end select + + albedo_initial=albedo + + if (update_albedo_from_ice) then + call interpolator_init( ice_interp, trim(ice_file_name)//'.nc', rad_lonb_2d, rad_latb_2d, data_out_of_bounds=(/CONSTANT/) ) + call read_ice_conc(Time) + call albedo_calc(albedo,Time) + else + if ( id_albedo > 0 ) used = send_data ( id_albedo, albedo ) + endif +else ! NTL 01/23 thermodynamic sea ice + if ( id_albedo > 0 ) used = send_data ( id_albedo, albedo, Time) +endif ! if(do_thermo_ice) ! Note: do_calc_eff_heat_cap is true by default and control when the surface ! heat capacity is calculated (land and ocean). @@ -570,6 +739,12 @@ subroutine mixed_layer ( & Time_next, & js, je, & t_surf, & + ! NTL 01/23 thermodynamic sea ice + h_thermo_ice, & + t_ml, & + const_correct, & + nudge_out, & + ! end ice addition flux_t, & flux_q, & flux_r, & @@ -583,7 +758,7 @@ subroutine mixed_layer ( & drdt_surf, & dhdt_atm, & dedq_atm, & - albedo_out) + albedo_out) ! ---- arguments ----------------------------------------------------------- type(time_type), intent(in) :: Time, Time_next @@ -591,6 +766,7 @@ subroutine mixed_layer ( & real, intent(in), dimension(:,:) :: net_surf_sw_down, surf_lw_down real, intent(in), dimension(:,:) :: flux_t, flux_q, flux_r real, intent(inout), dimension(:,:) :: t_surf +real, intent(inout), dimension(:,:) :: h_thermo_ice, t_ml, const_correct, nudge_out ! NTL 01/23 thermodynamic sea ice real, intent(in), dimension(:,:) :: dhdt_surf, dedt_surf, dedq_surf, & drdt_surf, dhdt_atm, dedq_atm real, intent(in) :: dt @@ -616,7 +792,9 @@ subroutine mixed_layer ( & land_ice_mask=land_mask endif -call albedo_calc(albedo_out,Time_next) +if (.not. do_thermo_ice) then ! NTL 01/23 thermodynamic sea ice + call albedo_calc(albedo_out,Time_next) +endif !s Add latent heat flux anomalies before any of the calculations take place @@ -665,17 +843,73 @@ subroutine mixed_layer ( & endif +if((do_nudge_ice.and.time_varying_nudge_ice).or. albedo_from_nudge_file) then + call interpolator( nudge_ice_interp, Time_next, h_thermo_ice_in, trim(nudge_ice_file_name) ) +endif + +if(do_prescribe_albedo.and.time_varying_albedo) then + call interpolator( albedo_interp, Time, albedo_in, trim(albedo_file_name) ) +endif + + + + + + ! ! Implement mixed layer surface boundary condition ! corrected_flux = - net_surf_sw_down - surf_lw_down + alpha_t * CP_AIR + alpha_lw - ocean_qflux t_surf_dependence = beta_t * CP_AIR + beta_lw + + + + + + +nudge_out = 0.0 +if (do_nudge_ice) then + where ((h_thermo_ice .gt. 0.0) .and. (h_thermo_ice_in .le. 0)) + nudge_out = h_thermo_ice / thermo_ice_nudge_timescale * L_thermo_ice + endwhere + corrected_flux = corrected_flux - nudge_out +endif + + + +! if (prescribe_nh_albedo_only) then +! where (rad_lat_2d .gt. 0) + +const_nudge_correction = 0.0 +const_correct = 0.0 +if (do_nudging_correction) then + const_nudge_correction = -area_weighted_global_mean(nudge_out) + if (.not. (const_nudge_correction .eq. 0.)) then + where (nudge_out .eq. 0.) + const_correct(:,:) = const_nudge_correction + endwhere + const_correct = const_correct * const_nudge_correction / area_weighted_global_mean(const_correct) + corrected_flux = corrected_flux - const_correct + endif +endif + if (evaporation) then corrected_flux = corrected_flux + alpha_q * HLV t_surf_dependence = t_surf_dependence + beta_q * HLV endif +if (do_thermo_ice) then + ! NTL 01/23 thermodynamic sea ice + corrected_flux_noq = corrected_flux + ocean_qflux + flux_thermo_ice = 0 ! Initialize conductive heat flux through sea ice + nudge_flux = 0. ! Initialize nudging flux + nudge_correction = 0. ! Initializing nudging correction + const_nudge_correction = 0. + ! for calculation of ice sfc temperature + dFdt_surf = dhdt_surf + drdt_surf + HLV * (dedt_surf) ! d(corrected_flux)/d(T_surf) +endif + !s Surface heat_capacity calculation based on that in MiMA by mj if(do_sc_sst) then !mj sst read from input file @@ -717,33 +951,169 @@ subroutine mixed_layer ( & endif -if (do_calc_eff_heat_cap) then - !s use the land_sea_heat_capacity calculated in mixed_layer_init +if (.not. do_thermo_ice) then ! NTL 01/23 thermodynamic sea ice switch - ! Now update the mixed layer surface temperature using an implicit step - ! - eff_heat_capacity = land_sea_heat_capacity + t_surf_dependence * dt !s need to investigate how this works + if (do_calc_eff_heat_cap) then + !s use the land_sea_heat_capacity calculated in mixed_layer_init - if (any(eff_heat_capacity .eq. 0.0)) then - write(*,*) 'mixed_layer: error', eff_heat_capacity - call error_mesg('mixed_layer', 'Avoiding division by zero',fatal) - end if + ! Now update the mixed layer surface temperature using an implicit step + ! + if (do_explicit) then ! NTL 01/23 thermodynamic sea ice -- add explicit option for consistency with thermo ice + eff_heat_capacity = land_sea_heat_capacity + else + eff_heat_capacity = land_sea_heat_capacity + t_surf_dependence * dt !s need to investigate how this works + endif - if(do_sc_sst.and.specify_sst_over_ocean_only) then - where (land_ice_mask) delta_t_surf = - corrected_flux * dt / eff_heat_capacity - where (land_ice_mask) t_surf = t_surf + delta_t_surf - else - delta_t_surf = - corrected_flux * dt / eff_heat_capacity - t_surf = t_surf + delta_t_surf - endif + if (any(eff_heat_capacity .eq. 0.0)) then + write(*,*) 'mixed_layer: error', eff_heat_capacity + call error_mesg('mixed_layer', 'Avoiding division by zero',fatal) + end if + + if(do_sc_sst.and.specify_sst_over_ocean_only) then + where (land_ice_mask) delta_t_surf = - corrected_flux * dt / eff_heat_capacity + where (land_ice_mask) t_surf = t_surf + delta_t_surf + else + delta_t_surf = - corrected_flux * dt / eff_heat_capacity + t_surf = t_surf + delta_t_surf + endif + + endif !s end of if(do_sc_sst). + +else if (do_thermo_ice) then + where (land_ice_mask) ! name misleading, for thermo sea ice, land_ice_mask is just land mask + delta_t_ml = - corrected_flux * dt / land_sea_heat_capacity + delta_h_thermo_ice = 0 + ! t_surf and t_ml are equal (they differ only when mixed layer is ice-covered) + t_surf = t_ml + delta_t_ml + ! = update state = + t_ml = t_ml + delta_t_ml + h_thermo_ice = h_thermo_ice + delta_h_thermo_ice + elsewhere + + ! NTL 01/23 thermodynamic sea ice + ! COPIED FROM Sea ice by Ian Eisenman, XZ 02/2018 + ! === (3) surface temperature increment is calculated with + ! explicit forward time step using flux corrected for implicit atm; + ! next, (4) surface state is updated === + + ! ====================================================================== + ! + ! Ocean mixed layer and sea ice model equations + ! + ! ====================================================================== + ! + ! Mixed layer temperature is evolved is t_ml, and sea ice thickness is + ! h_ice. Atmosphere cares about surface temparature (t_surf). Where + ! h_ice=0, t_surf=t_ml, but where h_ice>0, t_surf=t_ice which is the + ! ice surface temperature (assumed to be in steady-state: "zero layer + ! model"). So h_ice, t_ml, and t_surf are all saved at each time step + ! (t_ice does not need to be saved). + ! + ! This model is equivalent to the Semtner (1976) "zero-layer model", + ! with several differences (no snow, ice latent heat is same at base + ! as surface, ice surface temperature calculated by including sensible + ! and latent heat derivatives in dF/dT_surf, inclusion of frazil + ! growth). + + ! calculate values of delta_h_ice, delta_t_ml, and delta_t_ice + + + + + where ( h_thermo_ice .le. 0 ) ! = ice-free = [ should be equivalent to use .eq. instead of .le. ] + delta_t_ml = - ( corrected_flux_noq - ocean_qflux) * dt / land_sea_heat_capacity + delta_h_thermo_ice = 0 + elsewhere ! = ice-covered = ( h>0 ) + delta_t_ml = - ( thermo_ice_basal_flux_const * ( t_ml - t_thermo_ice_base ) - ocean_qflux ) * dt / land_sea_heat_capacity + delta_h_thermo_ice = ( corrected_flux_noq - thermo_ice_basal_flux_const * ( t_ml - t_thermo_ice_base ) ) * dt / L_thermo_ice + endwhere + + where ( t_ml + delta_t_ml .lt. t_surf_freeze ) ! = frazil growth = + delta_h_thermo_ice = delta_h_thermo_ice - ( t_ml + delta_t_ml - t_surf_freeze ) * (land_sea_heat_capacity)/L_thermo_ice + delta_t_ml = t_surf_freeze - t_ml + endwhere + + where ( ( h_thermo_ice .gt. 0 ) .and. ( h_thermo_ice + delta_h_thermo_ice .le. 0 ) ) ! = complete ablation = + delta_t_ml = delta_t_ml - ( h_thermo_ice + delta_h_thermo_ice ) * L_thermo_ice/land_sea_heat_capacity + delta_h_thermo_ice = - h_thermo_ice + endwhere + + ! = update surface temperature = + ! compute surface melt + where ( h_thermo_ice + delta_h_thermo_ice .gt. 0 ) ! surface is ice-covered + ! calculate increment in steady-state ice surface temperature + flux_thermo_ice = k_thermo_ice / (h_thermo_ice + delta_h_thermo_ice) * ( t_thermo_ice_base - t_surf ) + delta_t_thermo_ice = ( - corrected_flux_noq + flux_thermo_ice ) & + / ( k_thermo_ice / (h_thermo_ice + delta_h_thermo_ice) + dFdt_surf ) + where ( t_surf + delta_t_thermo_ice .gt. t_surf_freeze ) ! surface ablation + delta_t_thermo_ice = t_surf_freeze - t_surf + endwhere + ! surface is ice-covered, so update t_surf as ice surface temperature + delta_t_surf = delta_t_thermo_ice + t_surf = t_surf + delta_t_thermo_ice + elsewhere ! ice-free, so update t_surf as mixed layer temperature + delta_t_surf = delta_t_ml + t_surf = t_ml + delta_t_ml + endwhere + + + -endif !s end of if(do_sc_sst). + ! = update state = + t_ml = t_ml + delta_t_ml + h_thermo_ice = h_thermo_ice + delta_h_thermo_ice + + + endwhere + + + ! compute albedo + where (land_ice_mask) + albedo_out(:,:) = albedo_value * land_albedo_prefactor + elsewhere + where ( h_thermo_ice .gt. 0.0 ) + albedo_out(:,:) = thermo_ice_albedo + elsewhere + albedo_out(:,:) = albedo_value + endwhere + endwhere + + + if (do_prescribe_albedo) then + if (prescribe_nh_albedo_only) then + where (rad_lat_2d .gt. 0) + albedo_out(:,:) = albedo_in(:,:) + endwhere + else + albedo_out(:,:) = albedo_in(:,:) + endif + endif + + if ( id_albedo > 0 ) used = send_data ( id_albedo, albedo_out, Time_next ) +endif + +if ((do_prescribe_albedo).and.(.not.do_thermo_ice)) then + + albedo_out(:,:) = albedo_in(:,:) + if ( id_albedo > 0 ) used = send_data ( id_albedo, albedo_out, Time_next ) +endif ! if (do_thermo_ice) + +if (albedo_from_nudge_file) then + where ((h_thermo_ice .gt. 0) .and. (h_thermo_ice_in .le. 0)) + albedo_out(:,:) = albedo_value + endwhere +endif ! ! Finally calculate the increments for the lowest atmospheric layer ! -Tri_surf%delta_t = fn_t + en_t * delta_t_surf -if (evaporation) Tri_surf%delta_tr(:,:,nhum) = fn_q + en_q * delta_t_surf +if (do_explicit) then ! NTL 01/23 thermodynamic sea ice uses explicit step + Tri_surf%delta_t = fn_t + if (evaporation) Tri_surf%delta_tr(:,:,nhum) = fn_q +else + Tri_surf%delta_t = fn_t + en_t * delta_t_surf + if (evaporation) Tri_surf%delta_tr(:,:,nhum) = fn_q + en_q * delta_t_surf +endif ! ! Note: @@ -755,6 +1125,10 @@ subroutine mixed_layer ( & if(id_flux_oceanq > 0) used = send_data(id_flux_oceanq, ocean_qflux, Time_next) if(id_delta_t_surf > 0) used = send_data(id_delta_t_surf, delta_t_surf, Time_next) +! NTL 01/23 thermodynamic sea ice +if(id_h_thermo_ice > 0) used = send_data(id_h_thermo_ice, h_thermo_ice, Time_next) +if(id_t_ml > 0) used = send_data(id_t_ml, t_ml, Time_next) +if(id_flux_thermo_ice > 0) used = send_data(id_flux_thermo_ice, flux_thermo_ice, Time_next) end subroutine mixed_layer @@ -799,9 +1173,12 @@ subroutine read_ice_conc(Time) end subroutine read_ice_conc !================================================================================================================================= -subroutine mixed_layer_end(t_surf, bucket_depth, restart_file_bucket_depth) +subroutine mixed_layer_end(t_surf, & + h_thermo_ice, t_ml, const_correct, nudge_out, albedo, & ! NTL 01/23 thermodynamic sea ice + bucket_depth, restart_file_bucket_depth) real, intent(inout), dimension(:,:) :: t_surf +real, intent(inout), dimension(:,:) :: h_thermo_ice, t_ml, const_correct, nudge_out, albedo ! NTL 01/23 thermodynamic sea ice real, intent(inout), dimension(:,:,:) :: bucket_depth logical, intent(in) :: restart_file_bucket_depth integer:: unit @@ -811,6 +1188,15 @@ subroutine mixed_layer_end(t_surf, bucket_depth, restart_file_bucket_depth) ! write a restart file for the surface temperature call nullify_domain() call write_data(trim('RESTART/mixed_layer.res'), 't_surf', t_surf, grid_domain) +if (do_thermo_ice) then ! NTL 01/23 thermodynamic sea ice + call write_data(trim('RESTART/mixed_layer.res'), 'h_thermo_ice', h_thermo_ice, grid_domain) + call write_data(trim('RESTART/mixed_layer.res'), 't_ml', t_ml, grid_domain) + call write_data(trim('RESTART/mixed_layer.res'), 'const_correct', const_correct, grid_domain) + call write_data(trim('RESTART/mixed_layer.res'), 'nudge_out', nudge_out, grid_domain) + call write_data(trim('RESTART/mixed_layer.res'), 'albedo', albedo, grid_domain) +else if (do_prescribe_albedo) then + call write_data(trim('RESTART/mixed_layer.res'), 'albedo', albedo, grid_domain) +endif if (restart_file_bucket_depth) then call write_data(trim('RESTART/mixed_layer.res'), 'bucket_depth', bucket_depth, grid_domain) endif diff --git a/src/shared/constants/constants.F90 b/src/shared/constants/constants.F90 index 42ee900a5..f7fa83c2d 100644 --- a/src/shared/constants/constants.F90 +++ b/src/shared/constants/constants.F90 @@ -120,6 +120,7 @@ module constants_mod real, public, parameter :: RVGAS = 461.50 real, public, parameter :: CP_VAPOR = 4.0*RVGAS real, public, parameter :: DENS_H2O = 1000. +real, public, parameter :: DENS_ICE = 917. ! NTL 01/23 for sea ice code real, public, parameter :: HLV = 2.500e6 real, public, parameter :: HLF = 3.34e5 real, public, parameter :: HLS = HLV + HLF