diff --git a/components/forcing/src/forcing.F90 b/components/forcing/src/forcing.F90 index 72fa7b0e..35cd99ee 100644 --- a/components/forcing/src/forcing.F90 +++ b/components/forcing/src/forcing.F90 @@ -13,7 +13,10 @@ module forcing_mod use science_constants_mod, only: seconds_in_a_day use naming_conventions_mod use registry_mod, only : is_component_enabled - use logging_mod, only : LOG_ERROR, log_master_log, LOG_DEBUG, log_get_logging_level, log_log + use logging_mod, only : LOG_ERROR, log_master_log, LOG_DEBUG, log_get_logging_level, log_log, LOG_INFO + use conversions_mod, only : conv_to_string + + ! In order to set forcing from a netcdf file, need the following netcdf modules use netcdf, only : nf90_noerr, nf90_global, nf90_nowrite, & nf90_inquire_attribute, nf90_open, nf90_strerror, & @@ -27,9 +30,12 @@ module forcing_mod #endif character(len=*), parameter :: & - TIME_KEY = "time", & !< NetCDF data time key + TIME_KEY = "time", & !< NetCDF data time key Z_KEY = "z", & !< NetCDF data height(z) key - WSUBS_KEY = "wsubs" !< NetCDF data subsidence velocity + LEV_KEY = "lev", & !< NetCDF data pressure level key + TH_KEY = "theta_tendency", & !< NetCDF data theta tendency key + Q_KEY = "q_tendency", & !< NetCDF data water vapour tendency key + WSUBS_KEY = "wsubs" !< NetCDF data subsidence velocity key integer, parameter :: MAX_FILE_LEN=200 !< Maximum length of surface condition input filename character(MAX_FILE_LEN) :: input_file @@ -83,6 +89,8 @@ module forcing_mod logical :: relax_to_initial_theta_profile ! For relaxation, use initial profile as the target logical :: use_time_varying_subsidence ! Use time dependent subsidence veocity (read from file) + logical :: use_time_varying_theta ! Use time dependent theta forcing (read from file) + logical :: use_time_varying_q ! Use time dependent water vapour forcing (read from file) logical :: l_subs_pl_theta ! if .true. then subsidence applied to theta field logical :: l_subs_pl_q ! if .true. then subsidence applied to q fields @@ -114,6 +122,17 @@ module forcing_mod integer :: diagnostic_generation_frequency + ! Contains time varying forcing profile information + type time_varying_forcing_profile + real(kind=DEFAULT_PRECISION), allocatable :: forcing_times(:) ! input forcing times + real(kind=DEFAULT_PRECISION), allocatable :: forcing_values(:,:) ! input forcing values, interpolated to MONC heights + end type time_varying_forcing_profile + + type(time_varying_forcing_profile), allocatable :: time_varying_subsidence, time_varying_theta, time_varying_q + + logical :: convert_input_theta_from_temperature=.false. ! If .true. input forcing data is for temperature and should + ! be converted to theta (potential temerature). + public forcing_get_descriptor contains @@ -198,47 +217,51 @@ subroutine field_information_retrieval_callback(current_state, name, field_infor else if (name .eq. "th_subsidence") then field_information%enabled=current_state%th%active .and. l_subs_pl_theta .and. & allocated(current_state%global_grid%configuration%vertical%olzthbar) - else if (name .eq. "vapour_mmr_subsidence" .or. name .eq. "vapour_mmr_subsidence" .or. & - name .eq. "cloud_mmr_subsidence" .or. name .eq. "cloud_mmr_subsidence" ) then - field_information%enabled=.not. current_state%passive_q .and. & - current_state%number_q_fields .gt. 0 .and. l_subs_pl_q .and. & - allocated(current_state%global_grid%configuration%vertical%olzqbar) - else if (name .eq. "rain_mmr_subsidence" ) then - field_information%enabled=current_state%rain_water_mixing_ratio_index .gt. 0 .and. & - allocated(current_state%global_grid%configuration%vertical%olzqbar) - else if (name .eq. "ice_mmr_subsidence" ) then - field_information%enabled= current_state%ice_water_mixing_ratio_index .gt. 0 .and. & - allocated(current_state%global_grid%configuration%vertical%olzqbar) - else if (name .eq. "snow_mmr_subsidence" ) then - field_information%enabled= current_state%snow_water_mixing_ratio_index .gt. 0 .and. & - allocated(current_state%global_grid%configuration%vertical%olzqbar) - else if (name .eq. "graupel_mmr_subsidence" ) then - field_information%enabled= current_state%graupel_water_mixing_ratio_index .gt. 0 .and. & - allocated(current_state%global_grid%configuration%vertical%olzqbar) + else if (l_subs_pl_q) then + if (name .eq. "vapour_mmr_subsidence" .or. name .eq. "cloud_mmr_subsidence" ) then + field_information%enabled=.not. current_state%passive_q .and. & + current_state%number_q_fields .gt. 0 .and. & + allocated(current_state%global_grid%configuration%vertical%olzqbar) + else if (name .eq. "rain_mmr_subsidence" ) then + field_information%enabled=current_state%rain_water_mixing_ratio_index .gt. 0 .and. & + allocated(current_state%global_grid%configuration%vertical%olzqbar) + else if (name .eq. "ice_mmr_subsidence" ) then + field_information%enabled= current_state%ice_water_mixing_ratio_index .gt. 0 .and. & + allocated(current_state%global_grid%configuration%vertical%olzqbar) + else if (name .eq. "snow_mmr_subsidence" ) then + field_information%enabled= current_state%snow_water_mixing_ratio_index .gt. 0 .and. & + allocated(current_state%global_grid%configuration%vertical%olzqbar) + else if (name .eq. "graupel_mmr_subsidence" ) then + field_information%enabled= current_state%graupel_water_mixing_ratio_index .gt. 0 .and. & + allocated(current_state%global_grid%configuration%vertical%olzqbar) + end if + else if (name .eq. "u_large_scale") then field_information%enabled=current_state%u%active .and. l_constant_forcing_u else if (name .eq. "v_large_scale") then field_information%enabled=current_state%v%active .and. l_constant_forcing_v else if (name .eq. "th_large_scale") then field_information%enabled=current_state%th%active .and. l_constant_forcing_theta - else if (name .eq. "vapour_mmr_large_scale" .or. name .eq. "vapour_mmr_large_scale" .or. & - name .eq. "cloud_mmr_large_scale" .or. name .eq. "cloud_mmr_large_scale" ) then - field_information%enabled=.not. current_state%passive_q .and. & - current_state%number_q_fields .gt. 0 .and. l_subs_pl_q .and. & - allocated(current_state%global_grid%configuration%vertical%olzqbar) - else if (name .eq. "rain_mmr_large_scale" ) then - field_information%enabled=current_state%rain_water_mixing_ratio_index .gt. 0 .and. & - allocated(current_state%global_grid%configuration%vertical%olzqbar) - else if (name .eq. "ice_mmr_large_scale" ) then - field_information%enabled= current_state%ice_water_mixing_ratio_index .gt. 0 .and. & - allocated(current_state%global_grid%configuration%vertical%olzqbar) - else if (name .eq. "snow_mmr_large_scale" ) then - field_information%enabled= current_state%snow_water_mixing_ratio_index .gt. 0 .and. & - allocated(current_state%global_grid%configuration%vertical%olzqbar) - else if (name .eq. "graupel_mmr_large_scale" ) then - field_information%enabled= current_state%graupel_water_mixing_ratio_index .gt. 0 .and. & - allocated(current_state%global_grid%configuration%vertical%olzqbar) - end if + + else if (l_constant_forcing_q) then + if (name .eq. "vapour_mmr_large_scale" .or. name .eq. "cloud_mmr_large_scale" ) then + field_information%enabled=.not. current_state%passive_q .and. & + current_state%number_q_fields .gt. 0 .and. & + allocated(current_state%global_grid%configuration%vertical%olzqbar) + else if (name .eq. "rain_mmr_large_scale" ) then + field_information%enabled=current_state%rain_water_mixing_ratio_index .gt. 0 .and. & + allocated(current_state%global_grid%configuration%vertical%olzqbar) + else if (name .eq. "ice_mmr_large_scale" ) then + field_information%enabled= current_state%ice_water_mixing_ratio_index .gt. 0 .and. & + allocated(current_state%global_grid%configuration%vertical%olzqbar) + else if (name .eq. "snow_mmr_large_scale" ) then + field_information%enabled= current_state%snow_water_mixing_ratio_index .gt. 0 .and. & + allocated(current_state%global_grid%configuration%vertical%olzqbar) + else if (name .eq. "graupel_mmr_large_scale" ) then + field_information%enabled= current_state%graupel_water_mixing_ratio_index .gt. 0 .and. & + allocated(current_state%global_grid%configuration%vertical%olzqbar) + end if + end if ! Field information for 3d strcomp=INDEX(name, "forcing_3d_local") @@ -455,7 +478,7 @@ subroutine init_callback(current_state) real(kind=DEFAULT_PRECISION), dimension(:), allocatable :: f_force_pl_u ! Forcing values for u variable real(kind=DEFAULT_PRECISION), dimension(:), allocatable :: z_force_pl_u ! Forcing height values for u variable real(kind=DEFAULT_PRECISION), dimension(:), allocatable :: f_force_pl_v ! Forcing values for v variable - real(kind=DEFAULT_PRECISION), dimension(:), allocatable :: z_force_pl_v ! Forcing height values for v variabl + real(kind=DEFAULT_PRECISION), dimension(:), allocatable :: z_force_pl_v ! Forcing height values for v variable ! Read from netcdf file if used - always 2D real(kind=DEFAULT_PRECISION), dimension(:,:), allocatable :: f_subs_2d ! subsidence node for q variables @@ -470,9 +493,6 @@ subroutine init_callback(current_state) character(len=STRING_LENGTH) :: units_u_force='unset' ! units of theta variable forcing character(len=STRING_LENGTH) :: units_v_force='unset' ! units of theta variable forcing - logical :: convert_input_theta_from_temperature=.false. ! If .true. input forcing data is for temperature and should - ! be converted to theta (potential temerature). - integer :: k logical :: l_qdiag @@ -512,18 +532,44 @@ subroutine init_callback(current_state) if (current_state%graupel_water_mixing_ratio_index > 0) & iqg = current_state%graupel_water_mixing_ratio_index - ! Subsidence forcing initialization - + ! time_varying forcing initialization use_time_varying_subsidence= & options_get_logical(current_state%options_database, "use_time_varying_subsidence") + if ((l_subs_pl_theta .or. l_subs_pl_q) .and. use_time_varying_subsidence) then + allocate(time_varying_subsidence) + call init_time_varying_forcing(current_state, time_varying_subsidence, WSUBS_KEY, & + options_get_string(current_state%options_database, "varying_subsidence_file"), & + options_get_string(current_state%options_database, "varying_subsidence_coordinate")) + end if + + use_time_varying_theta= & + options_get_logical(current_state%options_database, "use_time_varying_theta") + if (use_time_varying_theta) then + convert_input_theta_from_temperature=options_get_logical(current_state%options_database, & + "convert_input_theta_from_temperature") + allocate(time_varying_theta) + call init_time_varying_forcing(current_state, time_varying_theta, TH_KEY, & + options_get_string(current_state%options_database, "varying_theta_file"), & + options_get_string(current_state%options_database, "varying_theta_coordinate")) + end if + + use_time_varying_q= & + options_get_logical(current_state%options_database, "use_time_varying_q") + if (use_time_varying_q) then + allocate(time_varying_q) + call init_time_varying_forcing(current_state, time_varying_q, Q_KEY, & + options_get_string(current_state%options_database, "varying_q_file"), & + options_get_string(current_state%options_database, "varying_q_coordinate")) + end if + + ! Subsidence forcing initialization + l_subs_pl_theta=options_get_logical(current_state%options_database, "l_subs_pl_theta") l_subs_pl_q=options_get_logical(current_state%options_database, "l_subs_pl_q") subsidence_input_type=options_get_integer(current_state%options_database, "subsidence_input_type") l_subs_local_theta=options_get_logical(current_state%options_database, "subsidence_local_theta") l_subs_local_q=options_get_logical(current_state%options_database, "subsidence_local_q") - input_file=options_get_string(current_state%options_database, "forcing_file") - if ((l_subs_pl_theta .and. .not. l_subs_local_theta) .or. & (l_subs_pl_q .and. .not. l_subs_local_q))then if (.not. is_component_enabled(current_state%options_database, "mean_profiles")) then @@ -531,49 +577,20 @@ subroutine init_callback(current_state) end if end if - if (l_subs_pl_theta .or. l_subs_pl_q)then - ! Read in the input_file - if (trim(input_file)=='' .or. trim(input_file)=='None') then - if (.not. use_time_varying_subsidence) then - allocate(z_subs_pl(options_get_array_size(current_state%options_database, "z_subs_pl")), & - f_subs_pl(options_get_array_size(current_state%options_database, "f_subs_pl"))) - call options_get_real_array(current_state%options_database, "z_subs_pl", z_subs_pl) - call options_get_real_array(current_state%options_database, "f_subs_pl", f_subs_pl) - ! Get profiles - zgrid=current_state%global_grid%configuration%vertical%z(:) - call piecewise_linear_1d(z_subs_pl(1:size(z_subs_pl)), f_subs_pl(1:size(f_subs_pl)), zgrid, & - current_state%global_grid%configuration%vertical%w_subs) - if (subsidence_input_type==DIVERGENCE) then - current_state%global_grid%configuration%vertical%w_subs(:) = & - -1.0*current_state%global_grid%configuration%vertical%w_subs(:)*zgrid(:) - end if - else - call log_master_log(LOG_ERROR, "timevarying forcing selected but no forcing file provided - STOP") - endif - deallocate(z_subs_pl, f_subs_pl) - else - if (use_time_varying_subsidence) then - call check_forcing_status(nf90_open(path = trim(input_file), mode = nf90_nowrite, ncid = ncid)) - if (log_get_logging_level() .ge. LOG_DEBUG) then - call log_master_log(LOG_DEBUG, "Reading in subsidence velocity profile from:"//trim(input_file) ) - end if - - call read_2d_forcing_dimensions(ncid, time_dim,z_dim) - allocate(forcing_input_times(time_dim), z_subs_pl(z_dim), f_subs_2d(z_dim, time_dim)) - call read_2d_forcing_variables(trim(input_file), ncid, time_dim, forcing_input_times, & - z_dim, z_subs_pl, WSUBS_KEY, f_subs_2d) - call check_forcing_status(nf90_close(ncid)) - - zgrid=current_state%global_grid%configuration%vertical%z(:) - ! interpolate forcing levels onto the MONC vertical grid (zgrid), for all forcing times - allocate(current_state%global_grid%configuration%vertical%wsubs_time_vary(size(zgrid), time_dim)) - call piecewise_linear_2d(z_subs_pl(1:z_dim), forcing_input_times(1:time_dim), & - f_subs_2d(1:z_dim,1:time_dim), zgrid, & - current_state%global_grid%configuration%vertical%wsubs_time_vary) - else - call log_master_log(LOG_ERROR, "constant forcing from file selected but not coded - STOP") - endif - endif + if ((l_subs_pl_theta .or. l_subs_pl_q) .and. .not. use_time_varying_subsidence) then + allocate(z_subs_pl(options_get_array_size(current_state%options_database, "z_subs_pl")), & + f_subs_pl(options_get_array_size(current_state%options_database, "f_subs_pl"))) + call options_get_real_array(current_state%options_database, "z_subs_pl", z_subs_pl) + call options_get_real_array(current_state%options_database, "f_subs_pl", f_subs_pl) + ! Get profiles + zgrid=current_state%global_grid%configuration%vertical%z(:) + call piecewise_linear_1d(z_subs_pl(1:size(z_subs_pl)), f_subs_pl(1:size(f_subs_pl)), zgrid, & + current_state%global_grid%configuration%vertical%w_subs) + if (subsidence_input_type==DIVERGENCE) then + current_state%global_grid%configuration%vertical%w_subs(:) = & + -1.0*current_state%global_grid%configuration%vertical%w_subs(:)*zgrid(:) + end if + deallocate(z_subs_pl, f_subs_pl) end if ! Time independent large-scale forcing (proxy for e.g. advection/radiation) @@ -596,7 +613,7 @@ subroutine init_callback(current_state) if (l_constant_forcing_theta)then constant_forcing_type_theta=options_get_integer(current_state%options_database, "constant_forcing_type_theta") forcing_timescale_theta=options_get_real(current_state%options_database, "forcing_timescale_theta") - l_constant_forcing_theta_z2pressure=options_get_logical(current_state%options_database, "l_constant_forcing_theta_z2pressure") + l_constant_forcing_theta_z2pressure=options_get_logical(current_state%options_database,"l_constant_forcing_theta_z2pressure") allocate(z_force_pl_theta(options_get_array_size(current_state%options_database, "z_force_pl_theta")), & f_force_pl_theta(options_get_array_size(current_state%options_database, "f_force_pl_theta"))) @@ -871,7 +888,11 @@ end subroutine init_callBack !! @param current_state The current model state_mod subroutine timestep_callback(current_state) type(model_state_type), target, intent(inout) :: current_state - integer :: current_x_index, current_y_index, target_x_index, target_y_index + integer :: current_x_index, current_y_index, target_x_index, target_y_index, k + logical :: calculate_diagnostics + real(kind=DEFAULT_PRECISION), dimension(current_state%local_grid%size(Z_INDEX)) :: temp_prof + + calculate_diagnostics = (mod(current_state%timestep, diagnostic_generation_frequency) == 0) current_x_index=current_state%column_local_x current_y_index=current_state%column_local_y @@ -921,38 +942,77 @@ subroutine timestep_callback(current_state) if (current_state%halo_column .or. current_state%timestep<3) return - if (mod(current_state%timestep, diagnostic_generation_frequency) == 0) then - call save_precomponent_tendencies(current_state, current_x_index, current_y_index, target_x_index, target_y_index) - end if + if (calculate_diagnostics) & + call save_precomponent_tendencies(current_state, current_x_index, current_y_index, target_x_index, target_y_index) ! AH: perform subsidence calculation but first determine if time varying or constant ! If timevarying then work out the profile of subsidence for the given time and ! assign to w_subs, which is used in apply_subsidence_to... ! if (use_time_varying_subsidence) then - call interpolate_point_linear_2d(forcing_input_times, & - current_state%global_grid%configuration%vertical%wsubs_time_vary, & + call interpolate_point_linear_2d(time_varying_subsidence%forcing_times, & + time_varying_subsidence%forcing_values, & current_state%time, current_state%global_grid%configuration%vertical%w_subs) endif ! if not w_subs is constant and set in the init_callback + + ! Apply time-varying theta and q (vapour only). + ! This functionality permits the user to apply a constant forcing separately as long as the + ! theta forcing is consistently in theta or absolute temperature units in both cases because + ! they share the convert_input_theta_from_temperature logical. + if (use_time_varying_theta) then + ! Obtain the profile, interpolated to the current time + call interpolate_point_linear_2d(time_varying_theta%forcing_times, & + time_varying_theta%forcing_values, & + current_state%time, temp_prof) + + ! Unit conversions + if (convert_input_theta_from_temperature)then ! Input is temperature not theta + temp_prof = temp_prof * current_state%global_grid%configuration%vertical%prefrcp(:) + end if + + ! Record the diagnostic and apply the forcing + dtheta_profile_diag = dtheta_profile_diag + temp_prof + do k=2,current_state%local_grid%size(Z_INDEX)-1 + current_state%sth%data(k,current_state%column_local_y,current_state%column_local_x) = & + current_state%sth%data(k,current_state%column_local_y,current_state%column_local_x) & + + temp_prof(k) + end do + endif ! use_time_varying_theta + + if (use_time_varying_q) then + ! Obtain the profile, interpolated to the current time + call interpolate_point_linear_2d(time_varying_q%forcing_times, & + time_varying_q%forcing_values, & + current_state%time, temp_prof) + + ! Record the diagnostic and apply the forcing + dq_profile_diag(:,iqv) = dq_profile_diag(:,iqv) + temp_prof + do k=2,current_state%local_grid%size(Z_INDEX)-1 + current_state%sq(iqv)%data(k,current_state%column_local_y,current_state%column_local_x) = & + current_state%sq(iqv)%data(k,current_state%column_local_y,current_state%column_local_x) & + + temp_prof(k) + end do + endif ! use_time_varying_q + + if (l_subs_pl_theta) then call apply_subsidence_to_flow_fields(current_state) call apply_subsidence_to_theta(current_state) end if if (l_subs_pl_q) call apply_subsidence_to_q_fields(current_state) - if (l_constant_forcing_theta)call apply_time_independent_forcing_to_theta(current_state) + if (l_constant_forcing_theta) call apply_time_independent_forcing_to_theta(current_state) #ifdef U_ACTIVE - if (l_constant_forcing_u)call apply_time_independent_forcing_to_u(current_state) + if (l_constant_forcing_u) call apply_time_independent_forcing_to_u(current_state) #endif #ifdef V_ACTIVE - if (l_constant_forcing_v)call apply_time_independent_forcing_to_v(current_state) + if (l_constant_forcing_v) call apply_time_independent_forcing_to_v(current_state) #endif - if (l_constant_forcing_q)call apply_time_independent_forcing_to_q(current_state) + if (l_constant_forcing_q) call apply_time_independent_forcing_to_q(current_state) - if (mod(current_state%timestep, diagnostic_generation_frequency) == 0) then - call compute_component_tendencies(current_state, current_x_index, current_y_index, target_x_index, target_y_index) - end if + if (calculate_diagnostics) & + call compute_component_tendencies(current_state, current_x_index, current_y_index, target_x_index, target_y_index) end subroutine timestep_callback @@ -1402,9 +1462,12 @@ end subroutine check_forcing_status !> Reads the dimensions for forcing from the NetCDF file. This routine assumes the forcing uses only time and height. !! @param ncid The NetCDF file id + !! @param vert_key The vertical coordinate key of the input data !! @param time_dim Number of elements in the time dimension - subroutine read_2d_forcing_dimensions(ncid, time_dim, z_dim) + !! @param z_dim Number of elements in the vertical dimension + subroutine read_2d_forcing_dimensions(ncid, vert_key, time_dim, z_dim) integer, intent(in) :: ncid + character(len=*), intent(in) :: vert_key integer, intent(out) :: time_dim integer, intent(out) :: z_dim integer :: time_dimid, z_dimid @@ -1412,7 +1475,7 @@ subroutine read_2d_forcing_dimensions(ncid, time_dim, z_dim) call check_forcing_status(nf90_inq_dimid(ncid, TIME_KEY, time_dimid)) call check_forcing_status(nf90_inquire_dimension(ncid, time_dimid, len=time_dim)) - call check_forcing_status(nf90_inq_dimid(ncid, Z_KEY, z_dimid)) + call check_forcing_status(nf90_inq_dimid(ncid, vert_key, z_dimid)) call check_forcing_status(nf90_inquire_dimension(ncid, z_dimid, len=z_dim)) end subroutine read_2d_forcing_dimensions @@ -1421,14 +1484,15 @@ end subroutine read_2d_forcing_dimensions !! @param ncid The id of the NetCDF file !! @param time_dim The number of elements in the time dimension !! @param time The time data field that is to be read + !! @param vert_key The vertical coordinate key of the input data !! @param z_dim is the number of elements in the height dimension !! @param force_2d_key is the string that defines the forcing variable in teh NetCDF file !! @param force_2d_var is the forcing data field that is read with dimension (t_dim, z_dim) - subroutine read_2d_forcing_variables(filename, ncid, time_dim, time, z_dim, z_profile, & + subroutine read_2d_forcing_variables(filename, ncid, time_dim, time, vert_key, z_dim, z_profile, & force_2d_key, force_2d_var ) character(*), intent(in) :: filename - character(len=*), intent(in) :: force_2d_key + character(len=*), intent(in) :: force_2d_key, vert_key integer, intent(in) :: ncid, time_dim, z_dim real(kind=DEFAULT_PRECISION), intent(inout) :: time(:), z_profile(:) real(kind=DEFAULT_PRECISION), dimension(:,:), allocatable, intent(inout) :: force_2d_var @@ -1446,11 +1510,11 @@ subroutine read_2d_forcing_variables(filename, ncid, time_dim, time, z_dim, z_pr call log_log(LOG_ERROR, "No recognized time variable found in"//trim(filename)) end if - status=nf90_inq_varid(ncid, Z_KEY, variable_id) + status=nf90_inq_varid(ncid, vert_key, variable_id) if (status==nf90_noerr)then - call read_single_forcing_variable(ncid, Z_KEY, data1d=z_profile) + call read_single_forcing_variable(ncid, vert_key, data1d=z_profile) else - call log_log(LOG_ERROR, "No recognized height variable found in"//trim(filename)) + call log_log(LOG_ERROR, "No recognized '"//trim(vert_key)//"' vertical coordinate variable found in"//trim(filename)) end if status=nf90_inq_varid(ncid, force_2d_key, variable_id) @@ -1491,4 +1555,70 @@ subroutine read_single_forcing_variable(ncid, key, data1d, data2d) end if end subroutine read_single_forcing_variable + !> Sets up time-varying forcing profiles + !! @param current_state Current model state + !! @param tvdata The time-varying data structure + !! @param key The variable key (name) to access + !! @param filename The input NetCDF file name + !! @param coordinate The vertical coordinate of the input data [ height | pressure ] + subroutine init_time_varying_forcing(current_state, tvdata, key, filename, coordinate) + type(model_state_type), target, intent(inout) :: current_state + type(time_varying_forcing_profile), intent(inout) :: tvdata + character(len=*), intent(in) :: key, filename, coordinate + + character(STRING_LENGTH) :: vert_key + integer :: ncid, time_dim_len, vert_dim_len ! Input file parameters + real(kind=DEFAULT_PRECISION), dimension(:), allocatable :: vert_coords ! contains input vertical coordinates + real(kind=DEFAULT_PRECISION), dimension(:,:), allocatable :: input_forcing + real(kind=DEFAULT_PRECISION), dimension(current_state%local_grid%size(Z_INDEX)) :: vert_grid + + ! Check for valid vertical coordinate specification + if (trim(coordinate) .eq. 'pressure') then + vert_key = LEV_KEY + if (key .eq. WSUBS_KEY) then + call piecewise_linear_1d(current_state%global_grid%configuration%vertical%zn(:), & + current_state%global_grid%configuration%vertical%prefn(:), & + current_state%global_grid%configuration%vertical%z(:), & + vert_grid) ! get pressure values on w-levels (z) + else + vert_grid = current_state%global_grid%configuration%vertical%prefn(:) + end if + else if (trim(coordinate) .eq. 'height') then + vert_key = Z_KEY + if (key .eq. WSUBS_KEY) then + vert_grid = current_state%global_grid%configuration%vertical%z(:) + else + vert_grid = current_state%global_grid%configuration%vertical%zn(:) + end if + else + call log_log(LOG_ERROR, "Must specify vertical coordinate for forcing file as 'height' [m] or 'pressure' [Pa] "// & + "for '"//trim(key)//"' of file: "//trim(filename)) + end if + + ! open forcing file + call check_forcing_status(nf90_open(path = trim(filename), mode = nf90_nowrite, ncid = ncid)) + + ! read the dimension sizes and allocate space to receive the data + call read_2d_forcing_dimensions(ncid, trim(vert_key), time_dim_len, vert_dim_len) + allocate(tvdata%forcing_times(time_dim_len), vert_coords(vert_dim_len), input_forcing(vert_dim_len, time_dim_len), & + tvdata%forcing_values(size(vert_grid), time_dim_len) ) + + ! read the forcing coordinates and data, then close file + call read_2d_forcing_variables(trim(filename), ncid, time_dim_len, tvdata%forcing_times, & + trim(vert_key), vert_dim_len, vert_coords, & + key, input_forcing) + call check_forcing_status(nf90_close(ncid)) + + + ! interpolate forcing levels onto the MONC vertical grid (vert_grid), for all forcing times + ! Linear gradient extrapolation beyond the input height bounds is implicitly performed. + ! This behaviour is different from that in piecewise_linear_1d, which does not extrapolate. + call piecewise_linear_2d(vert_coords, tvdata%forcing_times, input_forcing, & ! input coordinates and data + vert_grid, tvdata%forcing_values) ! output vertical coords and data + ! clean up + deallocate(vert_coords, input_forcing) + + end subroutine init_time_varying_forcing + + end module forcing_mod diff --git a/testcases/GASS_diurnal/PECAN.mcf b/testcases/GASS_diurnal/PECAN.mcf new file mode 100644 index 00000000..2ad9fa16 --- /dev/null +++ b/testcases/GASS_diurnal/PECAN.mcf @@ -0,0 +1,330 @@ +# Global configuration +global_configuration=global_config + +# Override global component defaults +fftsolver_enabled=.true. +pw_advection_enabled=.true. +simplesetup_enabled=.true. +smagorinsky_enabled=.true. +lower_bc_enabled=.true. +setfluxlook_enabled=.true. #This must be set to true if running with lower_bc +viscosity_enabled=.true. +diffusion_enabled=.true. +coriolis_enabled=.true. +damping_enabled=.true. +forcing_enabled=.true. +randomnoise_enabled=.true. +mean_profiles_enabled=.true. #This must be set to true if running with damping +casim_enabled=.true. +casim_profile_dgs_enabled=.true. +socrates_couple_enabled=.true. +th_advection_enabled=.true. +iobridge_enabled=.true. +profile_diagnostics_enabled=.true. +scalar_diagnostics_enabled=.true. +conditional_diagnostics_column_enabled=.true. +conditional_diagnostics_whole_enabled=.true. +pdf_analysis_enabled=.true. +subgrid_profile_diagnostics_enabled=.true. +flux_budget_enabled=.true. +l_lem_dissipation_rate=.false. + +registered=.true. # Print registered components and their version numbers to stdout +showcallbacks=.true. # Print registered callbacks in calling order to stdout + +cfl_monitor=.false. # Print dtm changes and cfl info to stdout + +# Specific diagnostic switches +l_cloud_mask=.true. # Enables 3d binary SOCRATES-based total cloud mask diagnostic +cloud_mask_method=DEFAULT # Cloud mask/fraction calculation method [ DEFAULT, SOCRATES ] + # DEFAULT is based on exceeding qlcrit and qicrit +l_partial_liq_ice=.true. # Calculate partial cloud fracions of liquid and ice, otherwise homogeneous binary + +# -------------------------------------------------- +# Parameters used to compute conditional diagnostics +# -------------------------------------------------- +# ncond: (automatically calculated) number of conditions entered under cond_request +# ndiag: (automatically calculated) number of diagnostics entered under diag_request +# : MUST include 'area' +# Resulting diagnostic array dimensions are (time, ndiag, ncond*2, nz) +# from 1:ncond cond=.true.; ncond+1:ncond*2 cond=.false. +# Debugging: when running the model, built with the cray debugger, the model will fail at the point of the +# IO reduction if ncond*2*ndiag > ~830. Just run with fewer to get past this. +cond_request=ALL, BYu, BCu, NrBCu, AC, ACu, ACd, WG1, WL1, ALu, ALd, CLu, CLd, AH, AL, AI, PPd, VPd, PVd, MO, BM, AA, AV +diag_request=area, W, W2, TH, WTH, THP, WTHP, THVP, WTHVP, THP2, WTHSG, W3, RH, U, V, WU, WV, WUSG, WVSG, TEMP, THL, THLP, THLP2, QVLI, QVLIP, QVLIP2, QRSG, QRSGP, QRSGP2, WQVLIP, WQRSGP +# critical thv and up/downdraft thresholds +thvprcrit=0.0 +wSupcrit=1.0 +wSdwncrit=-1.0 +wupcrit=0.0 +wdwncrit=0.0 +# critical ql, qi and q_hydrometeor for conditional sampling +# Cloud liquid water mixing ratio critical minimum to define cloud [kg/kg] +qlcrit= 1.e-5 # also affects cloud diagnostics elsewhere in model +qicrit= 1.e-5 # also affects cloud diagnostics elsewhere in model +qpptcrit= 1.e-5 +vpptcrit= 1.e-4 +# Diff calculations of thv; .false. means thv=th(1+0.61qv) and .true. mean thv=th(1+0.61qv-(ql+qi)) +thv_from_th_with_liqice=.true. +# -------------------------------------------------- + +# ----------------------------------------------------------------- +# Parameters used to compute vertical velocity critical thresholds +# ----------------------------------------------------------------- +# used only when pdf_analysis_enabled=.true. +# fractional percentiles: +# uppercrit=0.05 --> set updraft threshold at top 5% of w +# dwnpercrit=0.05 --> set downdraft threshold at bottom 5% of w +show_critical_w=.false. # Diagnostic printing +uppercrit=0.05 +dwnpercrit=0.05 + +# vertical velocity histogram parameters +n_w_bins=750 # number of bins (max-min)/bin_size (max determined by these parameters) +w_bin_size=0.1 # size of bin [m/s] +w_bin_min=-25.0 # lower bound of historgram [m/s] +# ----------------------------------------------------------------- + + +# Control configuration +display_synopsis_frequency=20 # Status output write frequency [ts] +termination_time=86400. # Model run end time [s] +dtm=1.0 # Initial model time step [s] + +# IO server configuration +#ioserver_configuration_file="io/io_cfg_files/data_write_with_conditionals.xml" +ioserver_configuration_file="io/io_cfg_files/paracon_level0_control.xml" +moncs_per_io_server=8 # Set to 8 for MetOffice machines + +time_basis=.true. # logical for sampling and output intervals [ .false. (timestep, DEFAULT) | .true. (time) ] + # - either unit is an integer (no fractional-second time intervals) +sampling_frequency=300 # Sampling interval for time averaging [ s | ts ] +3d_sampling_frequency=3600 # 3d Sampling interval for time averaging [ s | ts ] +mm=900. # An output interval [ s | ts ] +mm1=3600. # An output interval [ s | ts ] +diag_write_freq=3600. # Reinitialization interval for diagnostic files [s] + +diagnostic_file="diagnostic_files/PECAN_diagnostic_data.nc" # diagnostic file location and prefix + +# CRMStyle configuration +#crms_file="diagnostic_files/crmstyle_ts.nc" # CRMStyle file location and prefix + +# Checkpoint configuration +checkpoint_frequency=0 # Checkpoint file creation frequency [ts] +checkpoint_file="checkpoint_files/PECAN_dump.nc" # Checkpoint file location and prefix + +# Internal walltime configuration +check_walltime_frequency=100 # Frequency to check wall clock against walltime_limit [ts] +walltime_limit=00:18:00 # Internal wall clock time limit on simulation [hh:mm:ss] + +# Advection choices (choose pw or tvd component schemes for flow, theta, and q advection) +advection_flow_fields=pw +advection_theta_field=tvd +advection_q_fields=tvd + +# CFL configuration +cfl_frequency=17 # Frequency for checking CFL conditions [ts] +cfl_cvismax=0.4 # 'worst case' viscous stability parameter (eq 153 of lemdoc2.pdf) +cfl_cvelmax=0.4 # Largest advective Courant number (eq 152 of lemdoc2.pdf) +cfl_dtmmax=10.0 # Maximum time step [s] +cfl_dtmmin=0.001 # Minimum time step [s] + +# not convinced this works as fixing the gal +# adds the gal to the wind, this is not correct. +# Set as false for now +fix_ugal=.false. +ugal=-5.0 +fix_vgal=.false. +vgal=0.0 + +# Simple setup configuration +thref0=298.7259 +surface_pressure=96000. +surface_reference_pressure=100000. +x_size=64 +y_size=64 +z_size=99 +dxx=1000 +dyy=1000 +zztop=40000.0 +kgd=9,17,75,99 +hgd=500.,1500.,16000.,40000. +nsmth=20 +galilean_transformation=.false. + +enable_theta=.true. +use_anelastic_equations=.true. +origional_vertical_grid_setup=.true. +passive_th=.false. +passive_q=.false. +backscatter=.false. +use_viscosity_and_diffusion=.true. + +# Initialization of fields +l_init_pl_theta=.true. +z_init_pl_theta=0., 800., 1200.,3500.,4100.,8200.,12500.,13500.,14200.,16000.,20000.,24000.,28000.,32000.,36000.,40000. +f_init_pl_theta=297.0,297.0,300.0,306.5,311.0,318.0,328.5, 333.0, 340.0, 371.0, 483.0, 610.0, 738.0, 928.0, 1227.0,1447.0 +l_init_pl_u=.true. +z_init_pl_u=0.0, 40000. +f_init_pl_u=-5.0, -5.0 +l_init_pl_v=.false. +l_init_pl_q=.true. +names_init_pl_q=vapour +z_init_pl_q=0., 680., 1300., 3500., 4150., 4850., 5200., 6100., 7000., 8150., 9500., 10500., 11500., 12250., 13000., 14000., 18000., 40000. +f_init_pl_q=13.0e-3,12.5e-3,8.5e-3,4.3e-3,2.44e-3,1.52e-3,1.31e-3,0.75e-3,0.48e-3,0.28e-3,0.080e-3,0.038e-3,0.012e-3,0.008e-3,0.003e-3,0.0005e-3,0.0, 0.0 + +l_matchthref=.true. +l_thref_zero_buoy=.false. + +# Smagorinsky configuration +# Default values for the smagorinsky subgrid constants +# smag-subb=40.0 +# smag-subc=16.0 +# The subgrid constant values for the 'conventional' subgrid model +# of Brown (1999) +smag-subb=1.43 +smag-subc=1.43 + +# Random noise +l_rand_pl_theta=.true. +z_rand_pl_theta=0.0, 7000.0, 7001.0, 40000. +f_rand_pl_theta=0.1, 0.1, 0.0000, 0.0000 +names_rand_pl_q=vapour +z_rand_pl_q=0.0, 7000.0, 7001.0, 40000. +f_rand_pl_q=0.025e-3, 0.025e-3, 0.0000, 0.0000 + +# Simple cloud +max_height_cloud=30000. + +# physical constants +z0=0.0002 +z0th=0.0002 + +# Coriolis +fcoriol=9.07875e-5 # f = 2 * 7.292e-5 * sin(latitude_deg) [see radiation latitude for consistency] +geostrophic_wind_rate_of_change_in_x=0.0 +geostrophic_wind_rate_of_change_in_y=0.0 +surface_geostrophic_wind_x=-5.0 +surface_geostrophic_wind_y=0.0 + +# Damping configuration +dmptim=0.0002 +zdmp=20000.0 +hdmp=5000.0 + +# Time-varying forcing -------------------------------------------------------------------------------------------- +# enter files as the base directory-relative path to file +# specify the vertical coordinate of the forcing as 'height' [m] or 'pressure' [Pa] +# - expecting NetCDF variable names to be "wsubs", "theta_tendency", or "q_tendency" for +# subsidence, theta, and water vapour forcing time-height profiles, respectively. +# - expecting temporal coordinate variable/dimension to be called "time". +# - expecting wsubs as m/s +# - This is a subsidence velocity only, not divergence rate. +# - When using use_time_varying_subsidence=.true., subsidence_input_type is ignored. +# - expecting theta_tendency as K/s +# - Allowed to be temperature (instead of theta), but requires convert_input_theta_from_temperature=.true. +# - convert_input_theta_from_temperature will also apply to any additional constant forcing. +# - expecting q_tendecy as kg/kg/s +use_time_varying_subsidence=.false. +varying_subsidence_file= +varying_subsidence_coordinate= # 'height' [m] or 'pressure' [Pa] + +use_time_varying_theta=.true. +varying_theta_file=testcases/GASS_diurnal/pecan60varanaPECANC1.c1.20150601.000000_reformed.nc +varying_theta_coordinate=pressure # 'height' [m] or 'pressure' [Pa] +convert_input_theta_from_temperature=.true. + +use_time_varying_q=.true. +varying_q_file=testcases/GASS_diurnal/pecan60varanaPECANC1.c1.20150601.000000_reformed.nc +varying_q_coordinate=pressure # 'height' [m] or 'pressure' [Pa] + + +# Subsidence profile +l_subs_pl_theta=.false. +l_subs_pl_q=.false. + +# Large-scale forcing +# Add om a component to force theta +l_constant_forcing_theta=.false. +l_constant_forcing_q=.false. +l_constant_forcing_u=.true. +l_constant_forcing_v=.true. + +# TENDENCY=0, RELAXATION=1, INCREMENTS=2 +constant_forcing_type_theta=0 +constant_forcing_type_q=0 +constant_forcing_type_u=1 +constant_forcing_type_v=1 + +relax_to_initial_u_profile=.true. +relax_to_initial_v_profile=.true. + +forcing_timescale_u=21600. +forcing_timescale_v=21600. + +# Forcing profiles + +# surface flux config +# type_of_surface_boundary_conditions=PRESCRIBED_FLUX=0, PRESCRIBED_VALUES=1 +use_surface_boundary_conditions=.true. +use_time_varying_surface_values= .true. +surface_conditions_file=testcases/GASS_diurnal/pecan60varanaPECANC1.c1.20150601.000000_reformed.nc +type_of_surface_boundary_conditions = 0 + +#CASIM options +option=22222 +l_warm=.false. +aerosol_option=0 +iopt_act=0 +iopt_inuc=0 +process_level=0 +l_override_checks = .true. +number_q_fields=11 + + +#---------------------------------------------------------------------------------------------- +# SOCRATES inputs +mcc_temperature_profile = components/socrates_couple/data/mcc_profiles/one_km/mls.t.nc +mcc_vapour_profile = components/socrates_couple/data/mcc_profiles/one_km/mls.q.nc +mcc_ozone_profile = components/socrates_couple/data/mcc_profiles/one_km/mls.o3.nc +# Add options for rad_cntrl +spectral_file_lw = /projects/monc/fra23/socrates_spectra/ga7/sp_lw_ga7 +spectral_file_sw = /projects/monc/fra23/socrates_spectra/ga7/sp_sw_ga7 + +# 5 is clear sky, 2 is cloud (ice and liquid no overlap), 1 (ice and liquid full overlap) +i_cloud_representation = 2 + +## Time and location variables for socrates +l_360 = .false. # 360 days in year as opposed to 365 (a UM thing + # in the LEM, is this still required??) +l_solar_fixed = .false. # true equals fixed insolation using value in sol_fixed +l_no_solar = .false. +solar_fixed = 1361.0 # prescribed insolation value +sec_fixed = 1.15470054 # prescribed 1/cos(solar_zenith_angle) +latitude = 38.5 # latitude for the location of radiation calc [match to coriolis!] +longitude = -98.75 # longitude for the location of radiation calc +rad_start_year = 2015.0 # simulation year for earth sun distance +rad_start_day = 152.0 # day number from January 1st +rad_int_time = 300.0 # Radiation integration timestep [s]; if <= 0, call every timestep + # [trj: if time_basis=.true. above, this should get some consideration based on those timings + # ...probably. + # might be best to set it to the sampling_frequency] +rad_start_time = 0.0 # Start time for the radiation + # [trj: I think this is seconds in model time, and it should probably always be zero] + +## Surface albedo variables for socrates +l_variable_srf_albedo = .false. # not coded yet but will allow variable + # surface albedo with solar zenith angle +surface_albedo = 0.04 # surface albedo (fixed in time) + +mphys_nq_l=1 # cloud liquid mass +mphys_nd_l=0 # cloud drop number +mphys_nq_r=1 # rain mass +mphys_nq_i=1 # ice mass +mphys_nq_s=1 # snow mass +mphys_nq_g=1 # graupel mass + +l_fix_re = .true. +fixed_cloud_re = 10.0 # effective radius for cloud droplets 10 microns +fixed_ice_re = 30.0 # effective radius for ice 30 microns +#---------------------------------------------------------------------------------------------