The version for illustration is Held-Suarez model v3.4f. The source code is in /glade/u/home/lchengji/fms/idealized/exp/hs/v3.4f/src. Most of the modifications have comments starting with cliu30/cliu31/cliu33. The goal is to calculate the differential advection term on the RHS of (17) in the manscript.
The overriding in the dynamical core is optional. The switch that turns on/off the overriding is a namelist variable:
character(len=200) :: ovrd_option = 'none' ! The overriding option, can be 'full' or 'ubt' (barotropic)When ovrd_option is set to be 'none', the overriding is turned off and the model is standard. When it is set to be 'ubt', only the barotropic component of zonal-mean zonal wind will be overrode. When it is set to be 'full', the full zonal-mean zonal wind would be overrode.
For the 'full' overriding, the overriding wind consists of two parts. The first part is a perturbation field that doesn't change over time. This part is optional and is configured in the spectral_dynamics_init(...) subroutine. The second part of the overriding is a 6-hourly field read from another simulation. This part is configured in the spectral_dynamics(...) subroutine.
The type of perturbation is set by a namelist variable:
character(len=200) :: ovrd_pert_option = 'none' ! whether to add a perturbation to the overriding field or not. Can be 'const', 'gauss', or 'input'When ovrd_pert_option = 'gauss', the parameters for the Gaussian perturbation are set by these namelist variables:
! parameters for ovrd_pert_option = 'gauss' / 'const'
real :: lato_center = 0.0 ! the center latitude for Gaussian perturbation
real :: lato_width = 10.0 ! the meridional width of the Gaussian perturbation
real :: sigo_center = 0.5 ! the center vertical level
real :: sigo_width = 0.075 ! the vertical width
real :: uo_pert = 4.0 ! the amplitude of the Gaussian perturbation
! also used as the amplitude of the perturbation for ovrd_pert_option = 'const'When ovrd_pert_option = 'const', the constant perturbation is set by:
real :: uo_pert = 4.0 ! the amplitude of the Gaussian perturbation
! also used as the amplitude of the perturbation for ovrd_pert_option = 'const'When ovrd_pert_option = 'input', it means the perturbation is read from the first time step of a netCDF file. The file name is set by:
character(len=200) :: ovrd_pfile_name = '' ! the perturbation file name for ovrd_pert_option = 'input'For the part that reads in data 6-hourly, the file information is set by these variables:
character(len=200) :: ovrd_freq_option = '6hourly' ! the frequency of overriding, currently only support 6hourly
character(len=200) :: ovrd_folder = '' ! the folder path for (6-hourly) overriding files
character(len=200) :: ovrd_file_prefix = 'fullrad.' ! the prefix for overriding files
character(len=200) :: ovrd_file_suffix = '.atmos_4xday.nc' ! the suffix of overriding files
character(len=200) :: ovrd_file_name = '' ! this is generated by the model as prefix + yr_str + suffixThe following variables are declared to do the actual work of overriding. Their meaning is indicated in the comments. Note that some obsolete variables are still retained in the source code but not listed here.
real, allocatable, dimension(:, :, :) :: ug_diff ! the gridded difference between specified overriding u wind and the current u wind produced by the model
real, allocatable, dimension(:, :, :) :: ug_pert ! the gridded perturbation field
complex, allocatable, dimension(:, :, :) :: vors_diff ! the spectral difference between the vorticity of specified u wind and the current vorticity produced by the model
complex, allocatable, dimension(:, :, :) :: dt_vors_tmp ! a temperary variable to store vors tendency
real, allocatable, dimension(:, :, :) :: u_spcf ! gridded field of current specified u wind
real, allocatable, dimension(:, :, :) :: u_spcf_prev ! gridded field of previous 6-hourly u wind
real, allocatable, dimension(:, :, :) :: us_zm ! zonal mean of specified u wind
real, allocatable, dimension(:, :, :) :: dt_vorg_tmp ! a temperary variable to store vorg
integer :: id_us_zm ! id for outputting us_zm fieldThe code that setup the perturbation is in:
subroutine spectral_dynamics_init(Time, Time_step_in, tracer_attributes, dry_model_out, nhum_out, ocean_mask, n_days, n_years)The block of code reads:
!-- configure the perturbation (not a function of time) --
!- Gaussian -
if (uppercase(trim(ovrd_pert_option)) == 'GAUSS') then
allocate(temp_full(num_levels))
allocate(sig1d(num_levels))
allocate(temp_half(num_levels + 1))
call get_deg_lat(dg_lat)
call pressure_variables(temp_half, temp_half, sig1d, temp_full, reference_sea_level_press)
u_spcf = zonal_mean(ug(:, :, :, 1))
sig1d = sig1d / reference_sea_level_press
do k = 1, num_levels
do j = js, je
ug_pert(:, j, k) = uo_pert * exp(-((dg_lat(j) - lato_center) / lato_width)**2 &
-((sig1d(k) - sigo_center) / sigo_width)**2)
ug_pert(:, j, k) = ug_pert(:, j, k) + uo_pert * exp(-((dg_lat(j) + lato_center) / lato_width)**2 &
-((sig1d(k) - sigo_center) / sigo_width)**2)
end do
end do
deallocate(temp_full)
deallocate(temp_half)
deallocate(sig1d)
!- Specified by an input file -
else if (uppercase(trim(ovrd_pert_option)) == 'INPUT') then
call read_data(trim(ovrd_pfile_name), 'diff', ug_pert(:, :, :), grid_domain, timelevel = 1)
!note that the variable name is assumed to be 'diff'
ug_pert = zonal_mean(ug_pert)
!- Constant -
else if (uppercase(trim(ovrd_pert_option)) == 'CONST') then
ug_pert(:, :, :) = uo_pert
!- No perturbation -
else if (uppercase(trim(ovrd_pert_option)) == 'NONE') then
ug_pert = 0
else
call error_mesg('spectral_dynamics', trim(ovrd_pert_option)//'is not valid!', FATAL)
end ifNote that for ovrd_pert_option = 'input', the variable name in the input file must be diff.
The input file for overriding zonal-mean zonal wind contains 6-hourly wind fields. However, the model time step is much smaller than 6 hour. So at each model time step, the overriding wind is obtained by linear interpolation using two consecutive 6-hourly input u_spcf_prev and u_spcf.
In the:
subroutine spectral_dynamics_init(Time, Time_step_in, tracer_attributes, dry_model_out, nhum_out, ocean_mask, n_days, n_years)The u_spcf is initialized as:
!-- initialize u_spcf --
if (uppercase(trim(ovrd_option)) /= 'NONE') then
write(year_str, "(i5)") curr_year
!-- construct ovrd_file_name as prefix + year + suffix --
ovrd_file_name = trim(ovrd_folder) // '/' // trim(ovrd_file_prefix) // &
trim(adjustl(year_str)) // trim(ovrd_file_suffix)
if (file_exist(trim(ovrd_file_name))) then
call read_data(trim(ovrd_file_name), 'ucomp', u_spcf(:, :, :), grid_domain, timelevel = 1)
else
call error_mesg('spectral_dynamics', trim(ovrd_file_name)//' does not exist!', FATAL)
end if
end ifAfter the initialization, every 6 hour, the old u_spcf is assigned to u_spcf_prev and the new u_spcf is read in
subroutine spectral_dynamics(Time, psg_final, ug_final, vg_final, tg_final, tracer_attributes, grid_tracers_final, &
time_level_out, dt_psg, dt_ug, dt_vg, dt_tg, dt_tracers, wg_full, p_full, p_half, z_full)The code reads:
call get_time(Time, seconds, days)
if (uppercase(trim(ovrd_freq_option)) == '6HOURLY') then
if(mod(seconds, time_int) == 0 .and. file_exist(trim(ovrd_file_name))) then
! time_int = 6 hour (6 * 60 * 60 seconds)
u_spcf_prev = u_spcf
time_step = mod(max(days, 0), num_days) * 4 + seconds/time_int + 1
call read_data(trim(ovrd_file_name), 'ucomp', u_spcf(:, :, :), grid_domain, timelevel = time_step)
else if (.not.file_exist(trim(ovrd_file_name))) then
call error_mesg('spectral_dynamics', trim(ovrd_file_name)//' does not exist!', FATAL)
end if
else
...
end ifOther options of reading overriding field can be added in the else block.
ovrd_option determines the way to calculate ug_diff. For ovrd_option = 'ubt' (barotropic), the code reads:
else if(uppercase(trim(ovrd_option)) == 'UBT') then
us_zm = zonal_mean_vm(u_spcf)
ug_diff = zonal_mean_vm(time_interp(u_spcf_prev, u_spcf, seconds) - ug(:, :, :, current))
vors_diff = 0For ovrd_option = 'full', the code reads:
else if (uppercase(trim(ovrd_option)) == 'FULL') then
us_zm = zonal_mean(u_spcf) + ug_pert
ug_diff = zonal_mean(time_interp(u_spcf_prev, u_spcf, seconds) - ug(:, :, :, current)) + ug_pert
vors_diff = 0After ug_diff is calculated, the differential advection term is added to the tendency of temperature and vorticity respectively (see equation (17) in the manuscript):
!--- temperature advection ---
call horizontal_advection(ts(:,:,:,current), ug_diff, vg(:,:,:,current) * 0, dt_tg_tmp)
dt_tg_tmp = dt_tg_tmp + tg_diff
!==============
call trans_grid_to_spherical(dt_tg_tmp, dt_ts)dt_vorg_tmp = 0.
!--- vorticity advection ---
call horizontal_advection(vors(:, :, :, current), ug_diff, ug_diff * 0, dt_vorg_tmp)
call trans_grid_to_spherical(dt_vorg_tmp, dt_vors_tmp)
dt_vors = dt_vors + dt_vors_tmpThis concludes the modifications for overriding.