diff --git a/components/componentheaders.static b/components/componentheaders.static index e08b8e29..39cd7718 100644 --- a/components/componentheaders.static +++ b/components/componentheaders.static @@ -22,6 +22,7 @@ use meanprofiles_mod, only : meanprofiles_get_descriptor use modelsynopsis_mod, only : modelsynopsis_get_descriptor use petsc_solver_mod, only : petsc_solver_get_descriptor use pressuresource_mod, only : pressuresource_get_descriptor +use tvdadvection_mod, only : tvdadvection_get_descriptor use profile_diagnostics_mod, only : profile_diagnostics_get_descriptor use pstep_mod, only : pstep_get_descriptor use pwadvection_mod, only : pwadvection_get_descriptor @@ -32,19 +33,19 @@ use setfluxlook_mod, only : setfluxlook_get_descriptor use simplecloud_mod, only : simplecloud_get_descriptor use simplesetup_mod, only : simplesetup_get_descriptor use smagorinsky_mod, only : smagorinsky_get_descriptor +use socrates_couple_mod, only : socrates_couple_get_descriptor use stepfields_mod, only : stepfields_get_descriptor use steppingdirection_mod, only : steppingdirection_get_descriptor use subgrid_profile_diagnostics_mod, only : subgrid_profile_diagnostics_get_descriptor use swapsmooth_mod, only : swapsmooth_get_descriptor use terminationcheck_mod, only : terminationcheck_get_descriptor use thadvection_mod, only : thadvection_get_descriptor -use tvdadvection_mod, only : tvdadvection_get_descriptor use viscosity_mod, only : viscosity_get_descriptor use lwrad_exponential_mod, only : lwrad_exponential_get_descriptor use lateral_bcs_mod, only : lateral_bcs_get_descriptor use casim_mod, only : casim_get_descriptor use casim_profile_dgs_mod, only: casim_profile_dgs_get_descriptor -use socrates_couple_mod, only : socrates_couple_get_descriptor use conditional_diagnostics_column_mod, only : conditional_diagnostics_column_get_descriptor use conditional_diagnostics_whole_mod, only : conditional_diagnostics_whole_get_descriptor use pdf_analysis_mod, only : pdf_analysis_get_descriptor +use dephy_forcings_mod, only : dephy_forcings_get_descriptor diff --git a/components/componentregistrations.static b/components/componentregistrations.static index c6410161..7bad6af4 100644 --- a/components/componentregistrations.static +++ b/components/componentregistrations.static @@ -22,6 +22,7 @@ call add_component(component_descriptions, meanprofiles_get_descriptor()) call add_component(component_descriptions, modelsynopsis_get_descriptor()) call add_component(component_descriptions, petsc_solver_get_descriptor()) call add_component(component_descriptions, pressuresource_get_descriptor()) +call add_component(component_descriptions, tvdadvection_get_descriptor()) call add_component(component_descriptions, profile_diagnostics_get_descriptor()) call add_component(component_descriptions, pstep_get_descriptor()) call add_component(component_descriptions, pwadvection_get_descriptor()) @@ -32,19 +33,19 @@ call add_component(component_descriptions, setfluxlook_get_descriptor()) call add_component(component_descriptions, simplecloud_get_descriptor()) call add_component(component_descriptions, simplesetup_get_descriptor()) call add_component(component_descriptions, smagorinsky_get_descriptor()) +call add_component(component_descriptions, socrates_couple_get_descriptor()) call add_component(component_descriptions, stepfields_get_descriptor()) call add_component(component_descriptions, steppingdirection_get_descriptor()) call add_component(component_descriptions, subgrid_profile_diagnostics_get_descriptor()) call add_component(component_descriptions, swapsmooth_get_descriptor()) call add_component(component_descriptions, terminationcheck_get_descriptor()) call add_component(component_descriptions, thadvection_get_descriptor()) -call add_component(component_descriptions, tvdadvection_get_descriptor()) call add_component(component_descriptions, viscosity_get_descriptor()) call add_component(component_descriptions, lwrad_exponential_get_descriptor()) call add_component(component_descriptions, lateral_bcs_get_descriptor()) call add_component(component_descriptions, casim_get_descriptor()) call add_component(component_descriptions, casim_profile_dgs_get_descriptor()) -call add_component(component_descriptions, socrates_couple_get_descriptor()) call add_component(component_descriptions, conditional_diagnostics_column_get_descriptor()) call add_component(component_descriptions, conditional_diagnostics_whole_get_descriptor()) call add_component(component_descriptions, pdf_analysis_get_descriptor()) +call add_component(component_descriptions, dephy_forcings_get_descriptor()) diff --git a/components/conditional_diagnostics_whole/makefile b/components/conditional_diagnostics_whole/makefile index e73bfcdd..1c6fa06d 100644 --- a/components/conditional_diagnostics_whole/makefile +++ b/components/conditional_diagnostics_whole/makefile @@ -2,7 +2,8 @@ SRCSF = src/conditional_diagnostics_whole.F90 BUILDDIR=build COREDIR=../../model_core/build -FFLAGS=-I $(BUILDDIR) -I $(COREDIR) $(COMPILERFFLAGS) +COLUMNDIR= ../conditional_diagnostics_column/build/ +FFLAGS=-I $(BUILDDIR) -I $(COREDIR) -I $(COLUMNDIR) $(COMPILERFFLAGS) OBJS = $(patsubst %.F90,$(BUILDDIR)/%.o,$(SRCSF)) all: create-build-dirs $(OBJS) diff --git a/components/conditional_diagnostics_whole/src/conditional_diagnostics_whole.F90 b/components/conditional_diagnostics_whole/src/conditional_diagnostics_whole.F90 index 97f5935c..573a6bf0 100644 --- a/components/conditional_diagnostics_whole/src/conditional_diagnostics_whole.F90 +++ b/components/conditional_diagnostics_whole/src/conditional_diagnostics_whole.F90 @@ -11,7 +11,7 @@ module conditional_diagnostics_whole_mod use grids_mod, only : Z_INDEX use datadefn_mod, only : PRECISION_TYPE, DEFAULT_PRECISION use mpi, only : MPI_SUM, MPI_IN_PLACE, MPI_INT, MPI_REAL, MPI_DOUBLE - use missing_data_mod, only: rmdi + use optionsdatabase_mod, only : options_get_integer @@ -22,8 +22,8 @@ module conditional_diagnostics_whole_mod #endif integer :: diagnostic_generation_frequency - public conditional_diagnostics_whole_get_descriptor + Real, Parameter :: RMDI = -32768.0*32768.0 contains @@ -70,7 +70,7 @@ subroutine timestep_callback(current_state) PRECISION_TYPE, MPI_SUM, 0, current_state%parallel%monc_communicator, ierr) end if - !> Average summed diagnostics over the domain by dividing the total diagnostic for each condition + !> Average summed diagnostics over the domain by dividing the total diagnostic for each condition !! by the total number of points for the associated conditions. !! This is NOT done for the area diagnostic, identified by the requested_area position in the array. !! Note: CondDiags_tot(k, ncond, ndiag) @@ -106,7 +106,7 @@ subroutine timestep_callback(current_state) end subroutine timestep_callback - + !> Called on termination: currently doesn't need to do anything !! @param current_state The current model state subroutine finalisation_callback(current_state) diff --git a/components/dephy_forcings/makefile b/components/dephy_forcings/makefile new file mode 100644 index 00000000..78944064 --- /dev/null +++ b/components/dephy_forcings/makefile @@ -0,0 +1,17 @@ +SRCSF = src/dephy_forcings.F90 + +BUILDDIR=build +COREDIR=../../model_core/build +LOWERBCDIR=../lowerbc/build +SETFLUXLOOKDIR=../setfluxlook/build +GRIDMANAGERDIR=../gridmanager/build +FFLAGS=-I $(BUILDDIR) -I $(COREDIR) -I $(LOWERBCDIR) -I $(SETFLUXLOOKDIR) -I $(GRIDMANAGERDIR) $(COMPILERFFLAGS) +OBJS = $(patsubst %.F90,$(BUILDDIR)/%.o,$(SRCSF)) + +all: create-build-dirs $(OBJS) + +create-build-dirs: + mkdir -p $(BUILDDIR) + +$(OBJS) : $(BUILDDIR)/%.o : %.F90 + $(FTN) $(OPT) $(FFLAGS) $< -o $(BUILDDIR)/$(notdir $@) diff --git a/components/dephy_forcings/src/dephy_forcings.F90 b/components/dephy_forcings/src/dephy_forcings.F90 new file mode 100644 index 00000000..e119c506 --- /dev/null +++ b/components/dephy_forcings/src/dephy_forcings.F90 @@ -0,0 +1,1414 @@ +! Note that this subroutine replaces part of the profile initialisation +! as well as setfluxlook, coriolis, and "forcing from mcf" routines + +! NOTES +! - Note q in MONC is mixing ratio (rather than specific humidity, as is more usual) +! - Module needs to be initialised after gridmanager but before random noise +! - This may need to be checked + +! TODO +! - Check for possible problematic nature of modifying both current state and vertical grid simultaneously +! and check what "target" keyword does in this context. +! - Handle (evolving) surface pressure, and surface pressure initialisation from file +! - Deal with surface non-zero height above sea level? +! - Make code self-documenting with Doxygen +! - Add diagnostics +! - Implement lat/lon dependence for radiation (socrates_opt%latitude,socrates_opt%longitude,socrates_opt%l_variable_srf_albedo,socrates_opt%surface_albedo) +! - Handle damping in prescribed fashion (could be a question for DEPHY community)? +! - Improve interpolation routines? (replace by Steffen interpolation) +! - Code up finalisation callback (deallocation)? +! - Implement a better column mode check? +! - Discuss best way to do lowerbc and z0/z0th reinitialisation +! - Check grid vertical halo size +! - Discuss best way to refer to z-coordinates + +! TODO: RELATED ISSUES +! - Discuss possible problematic nature of modifying both current state and vertical grid simultaneously in general +! - Work on problems with "entire domain" setfluxlook for heterogeneous surface forcings. +! - More systematic implementation of utilities in a separate module. + +module dephy_forcings_mod + use datadefn_mod, only : STRING_LENGTH + use monc_component_mod, only : component_descriptor_type + ! Use PRESCRIBED_SURFACE_FLUXES, PRESCRIBED_SURFACE_VALUES to specify boundary conditions for use in other modules + use state_mod, only : model_state_type, PRESCRIBED_SURFACE_FLUXES, PRESCRIBED_SURFACE_VALUES + use optionsdatabase_mod, only : options_get_integer, options_get_logical, options_get_real, & + options_get_string + use grids_mod, only : vertical_grid_configuration_type, X_INDEX, Y_INDEX, Z_INDEX + use logging_mod, only : LOG_ERROR, log_master_log, log_log + use datadefn_mod, only : DEFAULT_PRECISION + ! note z0 and z0th are overwritten during the simulation + use science_constants_mod, only : cp, rlvap, z0, z0th, G, von_karman_constant, ratio_mol_wts,r_over_cp,& + alphah,betah,betam,gammah,gammam + use q_indices_mod, only: get_q_index, standard_q_names + use interpolation_mod, only: piecewise_linear_1d, interpolate_point_linear_1d, interpolate_point_linear_2d + use registry_mod, only : is_component_enabled + use netcdf, only : nf90_noerr, nf90_global, nf90_nowrite, & + nf90_inquire_attribute, nf90_open, nf90_strerror, & + nf90_inq_dimid, nf90_inquire_dimension, nf90_inq_varid, & + nf90_get_var, nf90_inquire, nf90_close, nf90_get_att, & + nf90_ebaddim, nf90_enotatt, nf90_enotvar + use configuration_checkpoint_netcdf_parser_mod, only : remove_null_terminator_from_string + ! use existing fluxlook functionality + use setfluxlook_mod, only : set_look, change_look + ! re-initialise lowerbc module as z0 and z0th change over time + use lowerbc_mod, only: tstrcona, rhmbc, ddbc, ddbc_x4, eecon, r2ddbc, rcmbc, tstrconb, & + x4con, xx0con, y2con, yy0con, viscous_courant_coefficient + use saturation_mod, only : qsaturation + ! supersede initialisation from mcf file in grid manager + use gridmanager_mod, only : set_up_vertical_reference_properties,set_anelastic_pressure, & + setup_reference_state_liquid_water_temperature_and_saturation, & + calculate_mixing_length_for_neutral_case, set_buoyancy_coefficient + + implicit none + +#ifndef TEST_MODE + private +#endif + + real(kind=DEFAULT_PRECISION), allocatable :: time_dephy(:) + real(kind=DEFAULT_PRECISION), allocatable :: height_dephy(:) + real(kind=DEFAULT_PRECISION), allocatable :: module_z(:) + real(kind=DEFAULT_PRECISION), allocatable :: module_zn(:) + real(kind=DEFAULT_PRECISION), allocatable :: full_theta(:,:,:) + + character(len=STRING_LENGTH) :: dephy_file + integer :: ncid_dephy + integer :: time_len_dephy + integer :: height_len_dephy + integer :: kkp !module wide parameter for vertical grid + logical :: l_verbose=.false. ! Temporary flag for dirty debugging + ! Three parameters below are meant to catch the model being in column mode. + integer :: column_check_x ! + integer :: column_check_y ! + integer :: n_dephy_passes=0 ! start checking after the initial two dephy passes + + ! Surface fields (time-dependent) + real(kind=DEFAULT_PRECISION), allocatable :: lat_traj_dephy(:), & + lon_traj_dephy(:), & + ps_forc_dephy(:), & + ts_dephy(:), & + sfc_sens_flx_dephy(:), & + sfc_lat_flx_dephy(:), & + z0_traj_dephy(:), & + z0th_traj_dephy(:), & + ustar_dephy(:), & + u_traj_dephy(:), & + v_traj_dephy(:), & + albedo_traj_dephy(:), & + q_skin_traj_dephy(:) + + ! Surface fields at time step + real(kind=DEFAULT_PRECISION) :: lat_traj, & + lon_traj, & + ps_forc, & + ts, & + sfc_sens_flx, & + sfc_lat_flx, & + z0_traj, & + z0th_traj, & + ustar, & + u_traj, & + v_traj, & + albedo_traj, & + q_skin_traj + + ! Initial fields on MONC vertical grid + real(kind=DEFAULT_PRECISION), allocatable :: u_dephy(:), & + v_dephy(:), & + theta_dephy(:), & + rv_dephy(:), & + tke_dephy(:) + + ! Forcing fields on MONC vertical grid (time-dependent) + real(kind=DEFAULT_PRECISION), allocatable :: height_forc_dephy(:,:), & + pressure_forc_dephy(:,:), & + ug_dephy(:,:), & + vg_dephy(:,:), & + u_adv_dephy(:,:), & + v_adv_dephy(:,:), & + theta_adv_dephy(:,:), & + theta_rad_dephy(:,:), & + rv_adv_dephy(:,:), & + w_dephy(:,:), & + theta_nudging_dephy(:,:), & + rv_nudging_dephy(:,:), & + u_nudging_dephy(:,:), & + v_nudging_dephy(:,:), & + nudging_inv_u_traj_dephy(:,:), & + nudging_inv_v_traj_dephy(:,:), & + nudging_inv_theta_traj_dephy(:,:), & + nudging_inv_rv_traj_dephy(:,:) + + ! Forcing fields during a time step + real(kind=DEFAULT_PRECISION), allocatable :: height_forc(:), & + pressure_forc(:), & + ug(:), & + vg(:), & + u_adv(:), & + v_adv(:), & + theta_adv(:), & + theta_rad(:), & + rv_adv(:), & + w(:), & + theta_nudging(:), & + rv_nudging(:), & + u_nudging(:), & + v_nudging(:), & + nudging_inv_u_traj(:), & + nudging_inv_v_traj(:), & + nudging_inv_theta_traj(:), & + nudging_inv_rv_traj(:) + + ! Dephy flags: applied only during simulation + integer :: int_adv_theta, & + int_adv_rv, & + int_rad_theta, & + int_forc_w, & + int_forc_geo, & + int_nudging_u, & + int_nudging_v, & + int_nudging_theta, & + int_nudging_rv + + ! Dephy strings + character(len=STRING_LENGTH):: str_surfaceType, & + str_surfaceForcing, & + str_surfaceForcingWind + + public dephy_forcings_get_descriptor + +contains + + !> Provides the descriptor back to the caller and is used in component registration + !! @returns The dephy_forcings component descriptor + type(component_descriptor_type) function dephy_forcings_get_descriptor() + dephy_forcings_get_descriptor%name="dephy_forcings" + dephy_forcings_get_descriptor%version=0.1 + dephy_forcings_get_descriptor%initialisation=>initialise_callback + dephy_forcings_get_descriptor%timestep=>timestep_callback + end function dephy_forcings_get_descriptor + + + !> Called during initialisation and will initialise the horizontal and vertical grid configurations + !! Note that the model state_mod (from a checkpoint or external file) must have been initialised already + !! @param current_state The current model state_mod + subroutine initialise_callback(current_state) + implicit none + type(model_state_type), target, intent(inout) :: current_state + type(vertical_grid_configuration_type) :: vertical_grid + integer :: alloc_z, alloc_y, alloc_x + + if(l_verbose) write(*,*) "initialising dephy" + + alloc_z=current_state%local_grid%size(Z_INDEX) + current_state%local_grid%halo_size(Z_INDEX) * 2 + alloc_y=current_state%local_grid%size(Y_INDEX) + current_state%local_grid%halo_size(Y_INDEX) * 2 + alloc_x=current_state%local_grid%size(X_INDEX) + current_state%local_grid%halo_size(X_INDEX) * 2 + + !!! Might need re-thinking: vertical grid passed into subroutines as well as current state + if (.not. current_state%initialised) then + call log_log(LOG_ERROR, "Must initialise the model state_mod before constructing the grid properties") + end if + + vertical_grid=current_state%global_grid%configuration%vertical + + ! get DEPHY filename, which needs to be trimmed + dephy_file=options_get_string(current_state%options_database, "dephy_file") + kkp=current_state%local_grid%size(Z_INDEX) + + ! allocate module arrays + allocate(module_z(kkp)) + allocate(module_zn(kkp)) + allocate(full_theta(alloc_z, alloc_y, alloc_x)) + allocate(height_forc(kkp)) + allocate(pressure_forc(kkp)) + allocate(ug(kkp)) + allocate(vg(kkp)) + allocate(u_adv(kkp)) + allocate(v_adv(kkp)) + allocate(theta_adv(kkp)) + allocate(theta_rad(kkp)) + allocate(rv_adv(kkp)) + allocate(w(kkp)) + allocate(theta_nudging(kkp)) + allocate(rv_nudging(kkp)) + allocate(u_nudging(kkp)) + allocate(v_nudging(kkp)) + allocate(nudging_inv_u_traj(kkp)) + allocate(nudging_inv_v_traj(kkp)) + allocate(nudging_inv_theta_traj(kkp)) + allocate(nudging_inv_rv_traj(kkp)) + + if(l_verbose) write(*,*) "initialised dephy 1" + + module_z=vertical_grid%z(:) + module_zn=vertical_grid%zn(:) + + if(l_verbose) write(*,*) "initialised dephy 2" + + call check_status(nf90_open(path = trim(dephy_file), mode = nf90_nowrite, ncid = ncid_dephy)) + + if(l_verbose) write(*,*) "initialised dephy 3" + + call dephy_read_dimension_variables() ! reads the forcing time and height variables + call dephy_read_profile_variables() ! does profile initialisation + call dephy_read_forcing_variables() ! reads and interpolates forcings + call dephy_read_surface_variables() ! reads and interpolates forcings + call dephy_read_integers() ! reads flags + call dephy_read_strings() ! reads strings + + if(l_verbose) write(*,*) "initialised dephy 4" + + call dephy_sanity_checks(current_state) ! checks for incompatible elements + + if(l_verbose) write(*,*) "initialised dephy 5" + + call dephy_time_interpolate(current_state) ! needed to get such things as z0 + + if(l_verbose) write(*,*) "initialised dephy 6" + + call dephy_profiles_etc(current_state, vertical_grid) ! initialises proiles/reference profiles + + if(l_verbose) write(*,*) "initialised dephy 7" + + call dephy_setfluxlook_init(current_state) ! initialises surface + + if(l_verbose) write(*,*) "initialised dephy 8" + + call check_status(nf90_close(ncid_dephy)) + + if(l_verbose) write(*,*) "initialised dephy 9" + + !call dephy_bughunting(current_state) + + end subroutine initialise_callback + + + subroutine timestep_callback(current_state) + implicit none + type(model_state_type), target, intent(inout) :: current_state + + ! update forcing and current fields + call dephy_time_interpolate(current_state) + call dephy_column_mode_check(current_state) + + ! Reset z0 and z0th, and related parameters + z0=z0_traj + z0th=z0th_traj + current_state%global_grid%configuration%vertical%zlogm=& + log(1.0_DEFAULT_PRECISION+current_state%global_grid%configuration%vertical%zn(2)/z0) + current_state%global_grid%configuration%vertical%zlogth=& + log((current_state%global_grid%configuration%vertical%zn(2)+z0)/z0th) + + if(l_verbose) write(*,*) "dephy timestep 1" + + call dephy_setfluxlook_timestep(current_state) + + if(l_verbose) write(*,*) "dephy timestep 2" + + call dephy_apply_forcings(current_state) + + if(l_verbose) write(*,*) "dephy timestep 3" + + ! re-initialisation needed as z0 and z0th means lowerbc needs reinitialisation + call lowerbc_reset_constants(current_state) + + if(l_verbose) write(*,*) "dephy timestep 4" + + call dephy_bughunting(current_state) + + end subroutine timestep_callback + + + subroutine dephy_sanity_checks(current_state) + implicit none + type(model_state_type), target, intent(inout) :: current_state + + !! Check for initialisation from mcf file (incompatible) + logical :: l_init_pl_u ! if .true. then initialize u field from mcf file + logical :: l_init_pl_v ! if .true. then initialize v field + logical :: l_init_pl_theta ! if .true. then initialize potential temperature field + logical :: l_init_pl_temp ! if .true. then initialize temperature field + logical :: l_init_pl_rh ! if .true. then initialize relative humidity field + logical :: l_init_pl_q ! if .true. then initialize q fields + real(kind=DEFAULT_PRECISION) :: termination_time, zztop, max_height_cloud + + !! Check if incompatible routines are enabled + !! Check if incompatible flags are set + + !! Check "standard" coriolis force not activated + !! Check "other" forcings routine not acticated + !! Check mean profiles are enabled if nudging is used + !! Check radiation scheme is compatible (ACTIVATED if 0, DEACTIVATED otherwise) + !! Check initialisation from profiles not activated + !! Check not current_state%passive_q .or. current_state%passive_th + !! Check current_state%number_q_fields > 0 + !! Check for inconsistent surface conditions. Make sure setfluxlook component is not active + !! Catch ustar-based setups (check earlier work on BOMEX for potential fix). + + if (is_component_enabled(current_state%options_database, "setfluxlook")) then + call log_master_log(LOG_ERROR, "DEPHY: setfluxlook component incompatible with dephy forcing") + endif + if (is_component_enabled(current_state%options_database, "forcing")) then + call log_master_log(LOG_ERROR, "DEPHY: forcing component incompatible with dephy forcing") + endif + if(is_component_enabled(current_state%options_database, "socrated_couple")) then + if(int_rad_theta==1) then + call log_master_log(LOG_ERROR, "DEPHY: socrates_couple component incompatible with dephy flag rad_theta==1") + endif + endif + if(.not. is_component_enabled(current_state%options_database, "socrated_couple")) then + if(int_rad_theta==0) then + call log_master_log(LOG_ERROR, "DEPHY: absence of socrates_couple component incompatible with dephy flag rad_theta==0") + endif + endif + if (is_component_enabled(current_state%options_database, "lwrad_exponential")) then + call log_master_log(LOG_ERROR, "DEPHY: lwrad_exponential component incompatible with dephy forcing") + endif + if (.not. is_component_enabled(current_state%options_database, "buoyancy")) then + call log_master_log(LOG_ERROR, "DEPHY: absence of buoyancy component incompatible with dephy forcing") + endif + if (.not. is_component_enabled(current_state%options_database, "lower_bc")) then + call log_master_log(LOG_ERROR, "DEPHY: absence of lower_bc component incompatible with dephy forcing") + endif + if (.not. is_component_enabled(current_state%options_database, "set_consistent_lowbc")) then + call log_master_log(LOG_ERROR, "DEPHY: absence of set_consistent_lowbc component incompatible with dephy forcing") + endif + if (.not. is_component_enabled(current_state%options_database, "mean_profiles")) then + call log_master_log(LOG_ERROR, "DEPHY: absence of mean_profiles component incompatible with dephy forcing") + endif + if(.not. (trim(str_surfaceForcing)=="surfaceFlux" .or. trim(str_surfaceForcing)=="ts")) then + call log_master_log(LOG_ERROR, "DEPHY: surfaceForcing (thermodynamics) not implemented") + endif + if(.not. trim(str_surfaceForcingWind)=="z0_traj") then + call log_master_log(LOG_ERROR, "DEPHY: surfaceForcingWind not implemented") + endif + if(.not. current_state%number_q_fields > 0) then + call log_master_log(LOG_ERROR, "DEPHY: dephy_forcings need current_state%number_q_fields > 0") + endif + if(current_state%passive_q) then + call log_master_log(LOG_ERROR, "DEPHY: dephy_forcings incompatible with passive_q") + endif + if(current_state%passive_th) then + call log_master_log(LOG_ERROR, "DEPHY: dephy_forcings incompatible with passive_th") + endif + l_init_pl_theta=options_get_logical(current_state%options_database, "l_init_pl_theta") + l_init_pl_temp=options_get_logical(current_state%options_database, "l_init_pl_temp") + l_init_pl_rh=options_get_logical(current_state%options_database, "l_init_pl_rh") + l_init_pl_q=options_get_logical(current_state%options_database, "l_init_pl_q") + l_init_pl_u=options_get_logical(current_state%options_database, "l_init_pl_u") + l_init_pl_v=options_get_logical(current_state%options_database, "l_init_pl_v") + if(l_init_pl_theta .or. l_init_pl_temp .or. l_init_pl_rh .or. l_init_pl_q .or. & + l_init_pl_u .or. l_init_pl_v) then + call log_master_log(LOG_ERROR, & + "DEPHY: dephy_forcings incompatible with initialisation of profiles using "//& + "l_init_pl_theta or l_init_pl_temp or l_init_pl_rh or l_init_pl_q or l_init_pl_u or l_init_pl_v") + endif + termination_time=options_get_real(current_state%options_database, "termination_time") + if(termination_time>time_dephy(time_len_dephy)) then + call log_master_log(LOG_ERROR, "DEPHY: termination time beyond last time in forcing file") + endif + zztop=options_get_real(current_state%options_database, "zztop") + if(zztop>height_dephy(height_len_dephy)) then + call log_master_log(LOG_ERROR, "DEPHY: zztop beyond highest level in forcing file") + endif + max_height_cloud=options_get_real(current_state%options_database, "max_height_cloud") + if(max_height_cloud0. .and. w_prof(kk+1)>0.) then + ! UPSIDENCE: GET TENDENCY USING LEVEL BELOW + tendency(kk,jj,ii)=tendency(kk,jj,ii)-w_prof(kk)*& + (field(kk,jj,ii)-field(kk-1,jj,ii))/(module_zn(kk)-module_zn(kk-1)) + else + ! NO CONSISTENT SIGN OF SUBSIDENCE PROFILE, USE MEAN GRADIENTS AND VELOCITIES? + tendency(kk,jj,ii)=tendency(kk,jj,ii)-0.5_DEFAULT_PRECISION*& + (w_prof(kk+1)+w_prof(kk))*(field(kk+1,jj,ii)-field(kk-1,jj,ii))/(module_zn(kk+1)-module_zn(kk-1)) + end if + end do + tendency(kkp,jj,ii)=tendency(kkp,jj,ii)-w_prof(kkp)*& + (field(kkp,jj,ii)-field(kkp-1,jj,ii))/(module_zn(kkp)-module_zn(kkp-1)) + end do + end do + + end subroutine + + ! Implements energy consistent (non-traditional) coriolis force using time-dependent geostrophic wind + subroutine dephy_coriolis(u,v,w,u_geo,v_geo,u_gal,v_gal,lat,su,sv,sw) + implicit none + real(kind=DEFAULT_PRECISION), intent(in) :: u(:,:,:),v(:,:,:),w(:,:,:) + real(kind=DEFAULT_PRECISION), intent(in) :: u_geo(:),v_geo(:) + real(kind=DEFAULT_PRECISION), intent(in) :: u_gal,v_gal + real(kind=DEFAULT_PRECISION), intent(in) :: lat + real(kind=DEFAULT_PRECISION), intent(inout) :: su(:,:,:),sv(:,:,:),sw(:,:,:) + real(kind=DEFAULT_PRECISION) :: fcoriol, fcoriol2 + real(kind=DEFAULT_PRECISION), parameter :: omega_earth=7.2921e-5 ! radial frecuency of earth's rotation + real(kind=DEFAULT_PRECISION), parameter :: proper_pi=atan(1.0_DEFAULT_PRECISION) * 4.0_DEFAULT_PRECISION + integer ii,jj,kk + + fcoriol=2.0_DEFAULT_PRECISION*omega_earth*sin(lat*proper_pi/180.0_DEFAULT_PRECISION) + ! Non-traditional coriolis terms, needed for energy-consistency + ! See e.g. Igel and Biello 2020 + ! Note geostrphic wind parametrises pressure gradients, and is therefore not in the non-traditional terms. + fcoriol2=2.0_DEFAULT_PRECISION*omega_earth*cos(lat*proper_pi/180.0_DEFAULT_PRECISION) + +#if defined(U_ACTIVE) && defined(V_ACTIVE) + do ii=2,size(su,3)-1 + do jj=2,size(su,2)-1 + do kk=2,kkp + su(kk, jj, ii)=su(kk, jj, ii)+fcoriol*& + (0.25_DEFAULT_PRECISION*(v(kk, jj, ii)+v(kk, jj, ii+1)+& + v(kk, jj-1, ii)+v(kk, jj-1, ii+1))+& + v_gal-v_geo(kk))-& + fcoriol2*& + (0.25_DEFAULT_PRECISION*(w(kk, jj, ii)+w(kk, jj, ii+1)+& + w(kk-1, jj, ii)+w(kk-1, jj, ii+1))) + + sv(kk, jj, ii)=sv(kk, jj, ii)-fcoriol*& + (0.25_DEFAULT_PRECISION*(u(kk, jj, ii)+u(kk, jj, ii-1)+& + u(kk, jj+1, ii)+u(kk, jj+1, ii-1))+& + u_gal-u_geo(kk)) + + end do + do kk=1,kkp-1 + sw(kk, jj, ii)=fcoriol2*& + (0.25_DEFAULT_PRECISION*(u(kk, jj, ii)+u(kk+1, jj, ii)+& + u(kk, jj, ii-1)+u(kk+1, jj, ii-1))+u_gal) + end do + end do + end do +#endif + + end subroutine dephy_coriolis + + + subroutine dephy_time_interpolate(current_state) + implicit none + type(model_state_type), target, intent(inout) :: current_state + + ! interpolate surface variables + call interpolate_point_linear_1d(time_dephy, lat_traj_dephy, current_state%time, lat_traj) + call interpolate_point_linear_1d(time_dephy, lon_traj_dephy, current_state%time, lon_traj) + call interpolate_point_linear_1d(time_dephy, ps_forc_dephy, current_state%time, ps_forc) + call interpolate_point_linear_1d(time_dephy, ts_dephy, current_state%time, ts) + call interpolate_point_linear_1d(time_dephy, sfc_sens_flx_dephy, current_state%time, sfc_sens_flx) + call interpolate_point_linear_1d(time_dephy, sfc_lat_flx_dephy, current_state%time, sfc_lat_flx) + call interpolate_point_linear_1d(time_dephy, z0_traj_dephy, current_state%time, z0_traj) + call interpolate_point_linear_1d(time_dephy, z0th_traj_dephy, current_state%time, z0th_traj) + call interpolate_point_linear_1d(time_dephy, ustar_dephy, current_state%time, ustar) + call interpolate_point_linear_1d(time_dephy, u_traj_dephy, current_state%time, u_traj) + call interpolate_point_linear_1d(time_dephy, v_traj_dephy, current_state%time, v_traj) + call interpolate_point_linear_1d(time_dephy, albedo_traj_dephy, current_state%time, albedo_traj) + call interpolate_point_linear_1d(time_dephy, q_skin_traj_dephy, current_state%time, q_skin_traj) + + ! interpolate height-dependent variables + call interpolate_point_linear_2d(time_dephy, height_forc_dephy, current_state%time, height_forc) + call interpolate_point_linear_2d(time_dephy, pressure_forc_dephy, current_state%time, pressure_forc) + call interpolate_point_linear_2d(time_dephy, ug_dephy, current_state%time, ug) + call interpolate_point_linear_2d(time_dephy, vg_dephy, current_state%time, vg) + call interpolate_point_linear_2d(time_dephy, u_adv_dephy, current_state%time, u_adv) + call interpolate_point_linear_2d(time_dephy, v_adv_dephy, current_state%time, v_adv) + call interpolate_point_linear_2d(time_dephy, theta_adv_dephy, current_state%time, theta_adv) + call interpolate_point_linear_2d(time_dephy, theta_rad_dephy, current_state%time, theta_rad) + call interpolate_point_linear_2d(time_dephy, rv_adv_dephy, current_state%time, rv_adv) + call interpolate_point_linear_2d(time_dephy, w_dephy, current_state%time, w) + call interpolate_point_linear_2d(time_dephy, theta_nudging_dephy, current_state%time, theta_nudging) + call interpolate_point_linear_2d(time_dephy, rv_nudging_dephy, current_state%time, rv_nudging) + call interpolate_point_linear_2d(time_dephy, u_nudging_dephy, current_state%time, u_nudging) + call interpolate_point_linear_2d(time_dephy, v_nudging_dephy, current_state%time, v_nudging) + call interpolate_point_linear_2d(time_dephy, nudging_inv_u_traj_dephy, current_state%time, nudging_inv_u_traj) + call interpolate_point_linear_2d(time_dephy, nudging_inv_v_traj_dephy, current_state%time, nudging_inv_v_traj) + call interpolate_point_linear_2d(time_dephy, nudging_inv_theta_traj_dephy, current_state%time, nudging_inv_theta_traj) + call interpolate_point_linear_2d(time_dephy, nudging_inv_rv_traj_dephy, current_state%time, nudging_inv_rv_traj) + end subroutine dephy_time_interpolate + + + subroutine dephy_apply_forcings(current_state) + implicit none + type(model_state_type), target, intent(inout) :: current_state + ! Surface fields, in time + integer:: iqv,nn + + iqv=get_q_index(standard_q_names%VAPOUR, 'dephy_forcings') + + ! CALCULATE COMPLETE THETA FOR CONVENIENCE + call dephy_add_profile(current_state%th%data,current_state%global_grid%configuration%vertical%thref,full_theta) + + ! APPLY FORCINGS + ! NUDGINGS NEED MEAN +#ifdef U_ACTIVE + if(int_nudging_u==1) then + call dephy_apply_nudging(current_state%global_grid%configuration%vertical%olubar,& + u_nudging,nudging_inv_u_traj,current_state%su%data) + end if +#endif +#ifdef V_ACTIVE + if(int_nudging_v==1) then + call dephy_apply_nudging(current_state%global_grid%configuration%vertical%olvbar,& + v_nudging,nudging_inv_v_traj,current_state%sv%data) + end if +#endif + if(int_nudging_theta==1) then + call dephy_apply_nudging(current_state%global_grid%configuration%vertical%olthbar+& + current_state%global_grid%configuration%vertical%thref,& + theta_nudging,nudging_inv_theta_traj,current_state%sth%data) + end if + if(int_nudging_rv==1) then + call dephy_apply_nudging(current_state%global_grid%configuration%vertical%olqbar(:,iqv),& + rv_nudging,nudging_inv_rv_traj,current_state%sq(iqv)%data) + end if + + ! PROFILES +#ifdef U_ACTIVE + call dephy_apply_tendency(u_adv,current_state%su%data) +#endif +#ifdef V_ACTIVE + call dephy_apply_tendency(v_adv,current_state%sv%data) +#endif + if(int_adv_theta==1) then + call dephy_apply_tendency(theta_adv,current_state%sth%data) + end if + if(int_adv_rv==1) then + call dephy_apply_tendency(rv_adv,current_state%sq(iqv)%data) + end if + + ! RADIATION TENDENCIES + if(int_rad_theta==1) then + call dephy_apply_tendency(theta_rad,current_state%sth%data) + end if + + ! LARGE-SCALE VERTICAL WIND + ! USE A DOWNWIND FORMULATION (AS IN DALES), + ! SO ADVECTION ONLY IN DIRECTION OF WIND + ! ALWAYS USE LOCAL GRADIENTS, AS NON-LOCAL ONES ARE UNPHYSICAL + ! APPLY TO ALL Q SPECIES + if(int_forc_w ==1) then +#ifdef U_ACTIVE + call dephy_apply_subsidence(w,current_state%u%data,current_state%su%data) +#endif +#ifdef V_ACTIVE + call dephy_apply_subsidence(w,current_state%v%data,current_state%sv%data) +#endif + end if + if((int_forc_w ==1) .OR. (int_forc_w==2)) then + call dephy_apply_subsidence(w,full_theta,current_state%sth%data) + DO nn=1,current_state%number_q_fields + call dephy_apply_subsidence(w,current_state%q(nn)%data,current_state%sq(nn)%data) + END DO + end if + + ! IMPLEMENTATION OF FULL CORIOLIS FORCE + if(int_forc_geo==1) then + call dephy_coriolis(current_state%u%data,current_state%v%data,current_state%w%data,& + ug,vg,current_state%ugal,current_state%vgal,lat_traj,& + current_state%su%data,current_state%sv%data,current_state%sw%data) + end if + + end subroutine dephy_apply_forcings + + + subroutine dephy_initial_profiles(current_state) + implicit none + type(model_state_type), target, intent(inout) :: current_state + integer :: iqv + logical :: l_matchthref + + iqv=get_q_index(standard_q_names%VAPOUR, 'dephy_forcings') + l_matchthref=options_get_logical(current_state%options_database, "l_matchthref") + if(l_matchthref) then + if(.not. current_state%use_anelastic_equations) then + call log_master_log(LOG_ERROR, "Non-anelastic equation set and l_maththref are incompatible") + end if + current_state%global_grid%configuration%vertical%thref(:)=theta_dephy + else + current_state%global_grid%configuration%vertical%thref(:)=current_state%thref0 + endif + + call dephy_set_profile(u_dephy,current_state%u%data) + call dephy_set_profile(u_dephy,current_state%zu%data) + call dephy_set_profile(v_dephy,current_state%v%data) + call dephy_set_profile(v_dephy,current_state%zv%data) + call dephy_set_profile(theta_dephy-current_state%global_grid%configuration%vertical%thref,current_state%th%data) + call dephy_set_profile(theta_dephy-current_state%global_grid%configuration%vertical%thref,current_state%zth%data) + ! Note q in MONC is mixing ratio (rather than specific humidity, as is more usual) + call dephy_set_profile(rv_dephy,current_state%q(iqv)%data) + call dephy_set_profile(rv_dephy,current_state%zq(iqv)%data) + + end subroutine dephy_initial_profiles + + !!! THIS SHOULD LIKELY COME FROM A UTILITIES MODULE + !> Will check a NetCDF status and write to log_log error any decoded statuses. Can be used to decode + !! whether a dimension or variable exists within the NetCDF data file + !! @param status The NetCDF status flag + !! @param foundFlag Whether the field has been found or not + subroutine check_status(status, found_flag) + integer, intent(in) :: status + logical, intent(out), optional :: found_flag + + if (present(found_flag)) then + found_flag = status /= nf90_ebaddim .and. status /= nf90_enotatt .and. status /= nf90_enotvar + if (.not. found_flag) return + end if + + if (status /= nf90_noerr) then + call log_log(LOG_ERROR, "NetCDF returned error code of "//trim(nf90_strerror(status))) + end if + end subroutine check_status + + !!! THIS SUBROUTINE REPLACES set_vertical_reference_profile IN gridmanager.F90 + !> Sets up the vertical grid reference profile at each point + !! @param current_state The current model state_mod + !! @param vertical_grid The vertical grid that we are working on + !! @param kkp Number of grid points in a vertical column + subroutine dephy_profiles_etc(current_state, vertical_grid) + implicit none + type(model_state_type), intent(inout) :: current_state + type(vertical_grid_configuration_type), intent(inout) :: vertical_grid + + integer :: k + + call dephy_initial_profiles(current_state) + call set_up_vertical_reference_properties(current_state, vertical_grid, current_state%global_grid%size(Z_INDEX)) + + if(l_verbose) write(*,*) "initialised dephy 6.3" + + call set_anelastic_pressure(current_state) + + ! THIS CRUCIAL STATEMENT IS HIDDEN IN set_qv_init_from_rh in gridmanager.F90 + vertical_grid=current_state%global_grid%configuration%vertical + + if(l_verbose) write(*,*) "initialised dephy 6.4" + + do k=2,kkp-1 + ! for diffusion onto p-level from below + vertical_grid%czb(k)=(vertical_grid%rho(k-1)/vertical_grid%rhon(k))/(vertical_grid%dz(k)*vertical_grid%dzn(k)) + ! for diffusion onto p-level from above + vertical_grid%cza(k)=(vertical_grid%rho(k)/vertical_grid%rhon(k))/(vertical_grid%dz(k)*vertical_grid%dzn(k+1)) + vertical_grid%czg(k)=-vertical_grid%czb(k)-vertical_grid%cza(k) + if (k .gt. 2) vertical_grid%czh(k)=vertical_grid%czb(k)*vertical_grid%cza(k-1) + end do + do k=2,kkp-1 + ! advection onto p-level from below + vertical_grid%tzc1(k)=0.25_DEFAULT_PRECISION*vertical_grid%rdz(k)*vertical_grid%rho(k-1)/vertical_grid%rhon(k) + ! advection onto p-level from above + vertical_grid%tzc2(k)=0.25_DEFAULT_PRECISION*vertical_grid%rdz(k)*vertical_grid%rho(k)/vertical_grid%rhon(k) + end do + do k=2,kkp-1 + ! advection onto w-level (K) from below + vertical_grid%tzd1(k)=0.25_DEFAULT_PRECISION*vertical_grid%rdzn(k+1)*vertical_grid%rhon(k)/vertical_grid%rho(k) + ! advection onto w-level (K) from above + vertical_grid%tzd2(k)=0.25_DEFAULT_PRECISION*vertical_grid%rdzn(k+1)*vertical_grid%rhon(k+1)/vertical_grid%rho(k) + end do + k=kkp + vertical_grid%czb(k)=(vertical_grid%rho(k-1)/vertical_grid%rhon(k))/(vertical_grid%dz(k)*vertical_grid%dzn(k)) + vertical_grid%cza(k)=0.0_DEFAULT_PRECISION + vertical_grid%czg(k)=-vertical_grid%czb(k) + vertical_grid%czh(k)=vertical_grid%czb(k)*vertical_grid%cza(k-1) + vertical_grid%tzc2(k)=0.25_DEFAULT_PRECISION*vertical_grid%rdz(k)*vertical_grid%rho(k)/vertical_grid%rhon(k) + vertical_grid%tzc1(k)=0.25_DEFAULT_PRECISION*vertical_grid%rdz(k)*vertical_grid%rho(k-1)/vertical_grid%rhon(k) + vertical_grid%czn=vertical_grid%dzn(2)*0.5_DEFAULT_PRECISION + vertical_grid%zlogm=log(1.0_DEFAULT_PRECISION+vertical_grid%zn(2)/z0) + vertical_grid%zlogth=log((vertical_grid%zn(2)+z0)/z0th) + vertical_grid%vk_on_zlogm=von_karman_constant/vertical_grid%zlogm + + if(l_verbose) write(*,*) "initialised dephy 6.5" + + call setup_reference_state_liquid_water_temperature_and_saturation(& + current_state, vertical_grid, current_state%global_grid%size(Z_INDEX)) + + if(l_verbose) write(*,*) "initialised dephy 6.6" + + call calculate_mixing_length_for_neutral_case(current_state, vertical_grid, current_state%global_grid%size(Z_INDEX)) + + if(l_verbose) write(*,*) "initialised dephy 6.7" + + call set_buoyancy_coefficient(current_state, vertical_grid, current_state%global_grid%size(Z_INDEX)) + + end subroutine dephy_profiles_etc + + !!! THIS SUBROUTINE REPLACES the initialisation_callback IN setfluxlook.F90 + subroutine dephy_setfluxlook_init(current_state) + type(model_state_type), intent(inout), target :: current_state + + integer, parameter :: LOOKUP_ENTRIES = 80 !< Number of entries for MO lookup tables + integer :: iqv + + current_state%lookup_table_entries=LOOKUP_ENTRIES + current_state%saturated_surface = .true. ! Copied from setfluxlook module + ! We will change this if we find some humidity data + allocate(current_state%lookup_table_velocity(current_state%lookup_table_entries), & + current_state%lookup_table_ustr(current_state%lookup_table_entries)) + + if(l_verbose) write(*,*) "initialised dephy 7.1" + + if (.not. allocated(current_state%cq))then + allocate(current_state%cq(current_state%number_q_fields)) + current_state%cq=0.0_DEFAULT_PRECISION + end if + + iqv = get_q_index(standard_q_names%VAPOUR, 'dephy_forcings') + current_state%cq(iqv) = ratio_mol_wts-1.0 + + if(l_verbose) write(*,*) "initialised dephy 7.2" + + if (trim(str_surfaceForcing) == "ts") then + current_state%type_of_surface_boundary_conditions = PRESCRIBED_SURFACE_VALUES + else if (trim(str_surfaceForcing) == "surfaceFlux") then + current_state%type_of_surface_boundary_conditions = PRESCRIBED_SURFACE_FLUXES + else + call log_master_log(LOG_ERROR, "Surface condition for dephy not implemented") + endif + + if(l_verbose) write(*,*) "initialised dephy 7.3" + + if (current_state%type_of_surface_boundary_conditions == PRESCRIBED_SURFACE_FLUXES) then + call dephy_set_flux(current_state) + + current_state%fbuoy=0. + current_state%fbuoy=& + current_state%global_grid%configuration%vertical%buoy_co(1)*current_state%surface_temperature_flux+& + current_state%cq(iqv)*current_state%surface_vapour_flux*G + call set_look(current_state) + current_state%theta_surf=0.0_DEFAULT_PRECISION + current_state%surface_vapour_mixing_ratio=0.0_DEFAULT_PRECISION + else if (current_state%type_of_surface_boundary_conditions == PRESCRIBED_SURFACE_VALUES) then ! Prescribed surface temperatures + call dephy_set_flux(current_state) + end if + + if(l_verbose) write(*,*) "initialised dephy 7.4" + + end subroutine dephy_setfluxlook_init + + !!! THIS SUBROUTINE REPLACES the timestep_callback IN setfluxlook.F90 + subroutine dephy_setfluxlook_timestep(current_state) + type(model_state_type), intent(inout), target :: current_state + + if(l_verbose) write(*,*) "dephy timestep 1.1" + call dephy_set_flux(current_state) + if(l_verbose) write(*,*) "dephy timestep 1.2" + call change_look(current_state) + if(l_verbose) write(*,*) "dephy timestep 1.3" + + end subroutine dephy_setfluxlook_timestep + + !!! THIS SUBROUTINE REPLACES the set_flux subroutine IN setfluxlook.F90 + subroutine dephy_set_flux(current_state) + type(model_state_type), intent(inout), target :: current_state + + real(kind=DEFAULT_PRECISION) :: surface_temp ! Surface temperature + integer :: iqv + iqv = get_q_index(standard_q_names%VAPOUR, 'dephy_forcings') + + if(l_verbose) write(*,*) "initialised dephy 7.3.1" + + if (current_state%type_of_surface_boundary_conditions == PRESCRIBED_SURFACE_FLUXES) then ! Prescribed surface fluxes + + if(l_verbose) write(*,*) "initialised dephy 7.3.1.0" + + current_state%surface_temperature_flux=sfc_sens_flx/(current_state%global_grid%configuration%vertical%rho(1)*cp) + current_state%surface_vapour_flux=sfc_lat_flx/(current_state%global_grid%configuration%vertical%rho(1)*rlvap) + + ! Update buoyancy flux... + current_state%fbuoynew=0.0_DEFAULT_PRECISION + current_state%fbuoynew=& + current_state%global_grid%configuration%vertical%buoy_co(1)*current_state%surface_temperature_flux + current_state%fbuoynew=current_state%fbuoynew+current_state%cq(iqv)*current_state%surface_vapour_flux*G + + else if (current_state%type_of_surface_boundary_conditions == PRESCRIBED_SURFACE_VALUES) then ! Prescribed surface temperatures + + if(l_verbose) write(*,*) "initialised dephy 7.3.1.1" + + if (current_state%saturated_surface)then + current_state%surface_vapour_mixing_ratio = qsaturation(ts,current_state%surface_pressure*0.01) + else + call log_master_log(LOG_ERROR, "DEPHY: prescribed surface vapout mixing ratio not implemented"//& + " (note q_skin is reservoir content!)") + end if + + if(l_verbose) write(*,*) "initialised dephy 7.3.1.2" + + ! Set theta_v + current_state%theta_surf = ts*& + (current_state%surface_reference_pressure/current_state%surface_pressure)**r_over_cp + current_state%theta_virtual_surf = current_state%theta_surf + current_state%theta_virtual_surf = current_state%theta_surf + & + current_state%global_grid%configuration%vertical%thref(2)* & + current_state%cq(iqv)*current_state%surface_vapour_mixing_ratio + + if(l_verbose) write(*,*) "initialised dephy 7.3.1.3" + + ! Finally set up new values of THVSURF dependent constants + current_state%cmbc=betam*current_state%global_grid%configuration%vertical%zn(2)*G*& + von_karman_constant/current_state%theta_virtual_surf + + if(l_verbose) write(*,*) "initialised dephy 7.3.1.4" + + current_state%rcmbc=1.0_DEFAULT_PRECISION/current_state%cmbc + current_state%ellmocon=current_state%theta_virtual_surf/(G*von_karman_constant) + + if(l_verbose) write(*,*) "initialised dephy 7.3.1.5" + + end if + + end subroutine dephy_set_flux + + ! lowerbc need to re-initialise due to evolving z0/z0th + ! It is possible to improve on this by making more changes to the lowerbc code + subroutine lowerbc_reset_constants(current_state) + type(model_state_type), target, intent(inout) :: current_state + + real(kind=DEFAULT_PRECISION) :: bhbc + + if ( current_state%use_surface_boundary_conditions .and. & + current_state%type_of_surface_boundary_conditions == PRESCRIBED_SURFACE_VALUES) then + ! variables below are only required when PRESCRIBED_SURFACE_VALUES are used. + tstrcona=von_karman_constant/alphah*current_state%global_grid%configuration%vertical%zlogth + bhbc=alphah*current_state%global_grid%configuration%vertical%zlogth + rhmbc=betah*(current_state%global_grid%configuration%vertical%zn(2)+z0-z0th)/& + (betam*current_state%global_grid%configuration%vertical%zn(2)) + ddbc=current_state%global_grid%configuration%vertical%zlogm*(bhbc-& + rhmbc*current_state%global_grid%configuration%vertical%zlogm) + ddbc_x4=4.*ddbc + r2ddbc=0.5_DEFAULT_PRECISION/ddbc + eecon=2.0_DEFAULT_PRECISION*rhmbc*current_state%global_grid%configuration%vertical%zlogm-bhbc + rcmbc=1.0_DEFAULT_PRECISION/current_state%cmbc + tstrconb=von_karman_constant/alphah + x4con=gammam*(current_state%global_grid%configuration%vertical%zn(2)+z0) + xx0con=gammam*z0 + y2con=gammah*(current_state%global_grid%configuration%vertical%zn(2)+z0) + yy0con=gammah*z0th + endif + + end subroutine lowerbc_reset_constants + + +!!! COPY OF PIECEWISE LINEAR INTERPOLATION ROUTINE, BUT INCLUDING THE K=1 LEVEL +!!! MAINLY TO BE SAFE + + !> Does a simple 1d linear interpolation to a point + !! @param zvals input z nodes + !! @param vals input nodal values + !! @param z location to interpolate onto + !! @param f output interpolated value + subroutine piecewise_linear_2d_k1(zvals, time_vals, vals, z, field) + + ! Assumes input variables (vals) are 2-D, with dims (z, time) + + real(kind=DEFAULT_PRECISION), intent(in) :: zvals(:), time_vals(:) + real(kind=DEFAULT_PRECISION), intent(in) :: vals(:,:) + real(kind=DEFAULT_PRECISION), intent(in) :: z(:) + real(kind=DEFAULT_PRECISION), intent(out) :: field(:,:) + + real(kind=DEFAULT_PRECISION) :: scale_tmp + + integer :: nn, k_monc, k_force ! loop counter + integer :: nz_force, nt_force, nz_monc, nt_monc ! time and height array sizes for forcing and monc grids + integer :: nnodes ! number of input values + + nz_force = size(zvals) + nt_force = size(time_vals) + nz_monc = size(z) + nt_monc = size(time_vals) ! time is intepolated in the timestep callback + + if ( zvals(1) .GT. zvals(nz_force) ) then ! pressure + call log_master_log(LOG_ERROR, "Input forcing uses pressure, this has not been coded"// & + " - please modify your forcing file to using height coordinates or modify the" // & + " interpolation routine in model_core to work with pressure coords - STOP") + else + do k_monc=1,nz_monc + do k_force=1,nz_force-1 + if( z(k_monc) >= zvals(k_force) .AND. z(k_monc) < zvals(k_force+1) ) then + scale_tmp = ( z(k_monc) - zvals(k_force) ) / & + ( zvals(k_force+1) - zvals(k_force) ) + do nn=1, nt_force + field(k_monc,nn) = vals(k_force,nn) + & + ( vals(k_force+1,nn) - vals(k_force,nn) ) & + * scale_tmp + enddo + endif + enddo + enddo + ! now examine the cases below and above forlevs(1) and forlevs(ktmfor + ! uses the local vertical gradient in the forcing to determine the + ! new values + do k_monc=1,nz_monc + if ( z(k_monc) >= zvals(nz_force) ) then + scale_tmp = ( z(k_monc) - zvals(nz_force) ) & + / ( zvals(nz_force) - zvals(nz_force-1) ) + do nn=1,nt_force + field(k_monc,nn) = vals(nz_force,nn) + & + ( vals(nz_force,nn) - vals(nz_force-1,nn) ) & + * scale_tmp + enddo + elseif ( z(k_monc) < zvals(1) )THEN + scale_tmp = ( z(k_monc) - zvals(1) ) & + / ( zvals(1) - zvals(2) ) + do nn=1,nt_force + field(k_monc,nn) = vals(1,nn) + & + ( vals(1,nn) - vals(2,nn) ) & + * scale_tmp + enddo + endif + enddo + ! + endif ! pressure or height + + end subroutine piecewise_linear_2d_k1 + + +integer function maxloc1(field) + real(kind=DEFAULT_PRECISION), intent(in), dimension(:,:,:) :: field + integer, dimension(3) :: maxloc_res + maxloc_res=maxloc(field) + maxloc1=maxloc_res(1) +end function maxloc1 + + +integer function minloc1(field) + real(kind=DEFAULT_PRECISION), intent(in), dimension(:,:,:) :: field + integer, dimension(3) :: minloc_res + minloc_res=minloc(field) + minloc1=minloc_res(1) +end function minloc1 + + +subroutine dephy_bughunting(current_state) + implicit none + type(model_state_type), intent(in) :: current_state + integer:: iqv + + iqv=get_q_index(standard_q_names%VAPOUR, 'dephy_forcings') + + ! tendency debugging + write (*,'(A)') 'DEPHY MANUAL DEBUGGING ROUTINE' + write (*,*) 'time ',current_state%time + + write (*,'(A)') ' su sv sw sth sqv' + write (*,'(A,5ES12.2)') 'max vals',maxval(current_state%su%data), & + maxval(current_state%sv%data), maxval(current_state%sw%data), & + maxval(current_state%sth%data), maxval(current_state%sq(iqv)%data) + write (*,'(A,5I12)') 'max loc',maxloc1(current_state%su%data), & + maxloc1(current_state%sv%data), maxloc1(current_state%sw%data), & + maxloc1(current_state%sth%data), maxloc1(current_state%sq(iqv)%data) + write (*,'(A)') ' ' + write (*,'(A,5ES12.2)') 'min zvals',minval(current_state%su%data), & + minval(current_state%sv%data), minval(current_state%sw%data), & + minval(current_state%sth%data), minval(current_state%sq(iqv)%data) + write (*,'(A,5I12)') 'min zloc',minloc1(current_state%su%data), & + minloc1(current_state%sv%data), minloc1(current_state%sw%data), & + minloc1(current_state%sth%data), minloc1(current_state%sq(iqv)%data) + write (*,'(A)') ' ' + write (*,'(A)') ' ' + + ! value debugging + write (*,'(A)') ' u v w th qv' + write (*,'(A,5ES12.2)') 'max vals',maxval(current_state%u%data), & + maxval(current_state%v%data), maxval(current_state%w%data), & + maxval(current_state%th%data), maxval(current_state%q(iqv)%data) + write (*,'(A,5I12)') 'max zloc',maxloc1(current_state%u%data), & + maxloc1(current_state%v%data), maxloc1(current_state%w%data), & + maxloc1(current_state%th%data), maxloc1(current_state%q(iqv)%data) + write (*,'(A)') ' ' + write (*,'(A,5ES12.2)') 'min vals',minval(current_state%u%data), & + minval(current_state%v%data), minval(current_state%w%data), & + minval(current_state%th%data), minval(current_state%q(iqv)%data) + write (*,'(A,5I12)') 'min zloc',minloc1(current_state%u%data), & + minloc1(current_state%v%data), minloc1(current_state%w%data), & + minloc1(current_state%th%data), minloc1(current_state%q(iqv)%data) + write (*,'(A)') ' ' + write (*,'(A)') ' ' + + !z value debugging + write (*,'(A)') ' zu zv zw zth zqv' + write (*,'(A,5ES12.2)') 'max vals',maxval(current_state%zu%data), & + maxval(current_state%zv%data), maxval(current_state%zw%data), & + maxval(current_state%zth%data), maxval(current_state%zq(iqv)%data) + write (*,'(A,5I12)') 'max zloc',maxloc1(current_state%zu%data), & + maxloc1(current_state%zv%data), maxloc1(current_state%zw%data), & + maxloc1(current_state%zth%data), maxloc1(current_state%zq(iqv)%data) + write (*,'(A)') ' ' + write (*,'(A,5ES12.2)') 'min vals',minval(current_state%zu%data), & + minval(current_state%zv%data), minval(current_state%zw%data), & + minval(current_state%zth%data), minval(current_state%zq(iqv)%data) + write (*,'(A,5I12)') 'min zloc',minloc1(current_state%zu%data), & + minloc1(current_state%zv%data), minloc1(current_state%zw%data), & + minloc1(current_state%zth%data), minloc1(current_state%zq(iqv)%data) + +end subroutine dephy_bughunting + +end module dephy_forcings_mod + +!~ !! SOME MORE DIRTY DIAGNOSTICS JUST ADDED AS COMMENTS + +!~ if(l_verbose) write(*,*) 'lat_traj' +!~ if(l_verbose) write(*,*) lat_traj +!~ if(l_verbose) write(*,*) 'lon_traj' +!~ if(l_verbose) write(*,*) lon_traj +!~ if(l_verbose) write(*,*) 'ps_forc' +!~ if(l_verbose) write(*,*) ps_forc +!~ if(l_verbose) write(*,*) 'ts' +!~ if(l_verbose) write(*,*) ts +!~ if(l_verbose) write(*,*) 'sfc_sens_flx' +!~ if(l_verbose) write(*,*) sfc_sens_flx +!~ if(l_verbose) write(*,*) 'sfc_lat_flx' +!~ if(l_verbose) write(*,*) sfc_lat_flx +!~ if(l_verbose) write(*,*) 'z0_traj' +!~ if(l_verbose) write(*,*) z0_traj +!~ if(l_verbose) write(*,*) 'z0th_traj' +!~ if(l_verbose) write(*,*) z0th_traj +!~ if(l_verbose) write(*,*) 'ustar' +!~ if(l_verbose) write(*,*) ustar +!~ if(l_verbose) write(*,*) 'u_traj' +!~ if(l_verbose) write(*,*) u_traj +!~ if(l_verbose) write(*,*) 'v_traj' +!~ if(l_verbose) write(*,*) v_traj +!~ if(l_verbose) write(*,*) 'albedo_traj' +!~ if(l_verbose) write(*,*) albedo_traj +!~ if(l_verbose) write(*,*) 'q_skin_traj' +!~ if(l_verbose) write(*,*) q_skin_traj + +!~ if(l_verbose) write(*,*) 'height_forc' +!~ if(l_verbose) write(*,*) height_forc +!~ if(l_verbose) write(*,*) 'pressure_forc' +!~ if(l_verbose) write(*,*) pressure_forc +!~ if(l_verbose) write(*,*) 'ug' +!~ if(l_verbose) write(*,*) ug +!~ if(l_verbose) write(*,*) 'vg' +!~ if(l_verbose) write(*,*) vg +!~ if(l_verbose) write(*,*) 'u_adv' +!~ if(l_verbose) write(*,*) u_adv +!~ if(l_verbose) write(*,*) 'v_adv' +!~ if(l_verbose) write(*,*) v_adv +!~ if(l_verbose) write(*,*) 'theta_adv' +!~ if(l_verbose) write(*,*) theta_adv +!~ if(l_verbose) write(*,*) 'theta_rad' +!~ if(l_verbose) write(*,*) theta_rad +!~ if(l_verbose) write(*,*) 'rv_adv' +!~ if(l_verbose) write(*,*) rv_adv +!~ if(l_verbose) write(*,*) 'w' +!~ if(l_verbose) write(*,*) w +!~ if(l_verbose) write(*,*) 'theta_nudging' +!~ if(l_verbose) write(*,*) theta_nudging +!~ if(l_verbose) write(*,*) 'rv_nudging' +!~ if(l_verbose) write(*,*) rv_nudging +!~ if(l_verbose) write(*,*) 'u_nudging' +!~ if(l_verbose) write(*,*) u_nudging +!~ if(l_verbose) write(*,*) 'v_nudging' +!~ if(l_verbose) write(*,*) v_nudging +!~ if(l_verbose) write(*,*) 'nudging_inv_u_traj' +!~ if(l_verbose) write(*,*) nudging_inv_u_traj +!~ if(l_verbose) write(*,*) 'nudging_inv_v_traj' +!~ if(l_verbose) write(*,*) nudging_inv_v_traj +!~ if(l_verbose) write(*,*) 'nudging_inv_theta_traj' +!~ if(l_verbose) write(*,*) nudging_inv_theta_traj +!~ if(l_verbose) write(*,*) 'nudging_inv_rv_traj' +!~ if(l_verbose) write(*,*) nudging_inv_rv_traj diff --git a/components/gridmanager/src/gridmanager.F90 b/components/gridmanager/src/gridmanager.F90 index 16d3d1ac..9fe9bf84 100644 --- a/components/gridmanager/src/gridmanager.F90 +++ b/components/gridmanager/src/gridmanager.F90 @@ -19,7 +19,7 @@ module gridmanager_mod #ifndef TEST_MODE private #endif - + ! 1 = No adjustment ! 2 = ensure P0=PSF by adjusting THREF profile by constant factor ! 3 = ensure P0=PSF by adjusting PSF (not advised) @@ -28,7 +28,9 @@ module gridmanager_mod real, parameter :: DEFAULT_SPACING = 1.E9 !< The default spacing used if no grid is active in a specific dimension real(kind=DEFAULT_PRECISION), dimension(:,:), allocatable :: qinit - public gridmanager_get_descriptor + public gridmanager_get_descriptor,set_up_vertical_reference_properties,set_anelastic_pressure, & + setup_reference_state_liquid_water_temperature_and_saturation, & + calculate_mixing_length_for_neutral_case, set_buoyancy_coefficient contains @@ -57,7 +59,7 @@ subroutine initialise_callback(current_state) call initialise_verticalgrid_configuration_type(current_state) dimensions=1 if (current_state%global_grid%active(X_INDEX)) dimensions = dimensions+1 - if (current_state%global_grid%active(Y_INDEX)) dimensions = dimensions+1 + if (current_state%global_grid%active(Y_INDEX)) dimensions = dimensions+1 call log_master_log(LOG_INFO, trim(conv_to_string(dimensions))//"D system; z="//& trim(conv_to_string(current_state%global_grid%size(Z_INDEX)))//", y="//& trim(conv_to_string(current_state%global_grid%size(Y_INDEX)))//", x="//& @@ -84,7 +86,7 @@ subroutine finalise_callback(current_state) vertical_grid%q_rand, vertical_grid%theta_rand, vertical_grid%w_subs, vertical_grid%w_rand, & vertical_grid%q_force, vertical_grid%theta_force, vertical_grid%u_force, vertical_grid%v_force & ) - end subroutine finalise_callback + end subroutine finalise_callback !> Will initialise the vertical grid configuration !! @param current_state The current model state_mod @@ -100,7 +102,7 @@ subroutine initialise_verticalgrid_configuration_type(current_state) current_state%origional_vertical_grid_setup, current_state%continuation_run) call set_vertical_reference_profile(current_state, current_state%global_grid%configuration%vertical, & current_state%global_grid%size(Z_INDEX)) - + end subroutine initialise_verticalgrid_configuration_type !> Sets up the vertical grid reference profile at each point @@ -117,7 +119,7 @@ subroutine set_vertical_reference_profile(current_state, vertical_grid, kkp) call calculate_initial_profiles(current_state, vertical_grid) call set_up_vertical_reference_properties(current_state, vertical_grid, current_state%global_grid%size(Z_INDEX)) call set_anelastic_pressure(current_state) - ! + ! call set_qv_init_from_rh(current_state) do k=2,kkp-1 @@ -130,9 +132,9 @@ subroutine set_vertical_reference_profile(current_state, vertical_grid, kkp) end do do k=2,kkp-1 ! advection onto p-level from below - vertical_grid%tzc1(k)=0.25_DEFAULT_PRECISION*vertical_grid%rdz(k)*vertical_grid%rho(k-1)/vertical_grid%rhon(k) + vertical_grid%tzc1(k)=0.25_DEFAULT_PRECISION*vertical_grid%rdz(k)*vertical_grid%rho(k-1)/vertical_grid%rhon(k) ! advection onto p-level from above - vertical_grid%tzc2(k)=0.25_DEFAULT_PRECISION*vertical_grid%rdz(k)*vertical_grid%rho(k)/vertical_grid%rhon(k) + vertical_grid%tzc2(k)=0.25_DEFAULT_PRECISION*vertical_grid%rdz(k)*vertical_grid%rho(k)/vertical_grid%rhon(k) end do do k=2,kkp-1 ! advection onto w-level (K) from below @@ -146,7 +148,7 @@ subroutine set_vertical_reference_profile(current_state, vertical_grid, kkp) vertical_grid%czg(k)=-vertical_grid%czb(k) vertical_grid%czh(k)=vertical_grid%czb(k)*vertical_grid%cza(k-1) vertical_grid%tzc2(k)=0.25_DEFAULT_PRECISION*vertical_grid%rdz(k)*vertical_grid%rho(k)/vertical_grid%rhon(k) - vertical_grid%tzc1(k)=0.25_DEFAULT_PRECISION*vertical_grid%rdz(k)*vertical_grid%rho(k-1)/vertical_grid%rhon(k) + vertical_grid%tzc1(k)=0.25_DEFAULT_PRECISION*vertical_grid%rdz(k)*vertical_grid%rho(k-1)/vertical_grid%rhon(k) vertical_grid%czn=vertical_grid%dzn(2)*0.5_DEFAULT_PRECISION vertical_grid%zlogm=log(1.0_DEFAULT_PRECISION+vertical_grid%zn(2)/z0) vertical_grid%zlogth=log((vertical_grid%zn(2)+z0)/z0th) @@ -191,7 +193,7 @@ subroutine calculate_initial_profiles(current_state, vertical_grid) logical :: l_matchthref ! if .true. then initialize thref to be the same as theta_init character(len=STRING_LENGTH), dimension(:), allocatable :: names_init_pl_q ! names of q variables to initialize - + real(kind=DEFAULT_PRECISION), allocatable :: f_init_pl_q_tmp(:) !temporary 1D storage of initial q field real(kind=DEFAULT_PRECISION), allocatable :: zgrid(:) ! z grid to use in interpolation @@ -199,7 +201,7 @@ subroutine calculate_initial_profiles(current_state, vertical_grid) real(kind=DEFAULT_PRECISION) :: qsat allocate(zgrid(current_state%local_grid%local_domain_end_index(Z_INDEX))) - + zztop = current_state%global_grid%top(Z_INDEX) ! Initialize everything to zero. This won't make sense for theta. @@ -209,13 +211,13 @@ subroutine calculate_initial_profiles(current_state, vertical_grid) vertical_grid%theta_init = 0.0_DEFAULT_PRECISION l_init_pl_theta=options_get_logical(current_state%options_database, "l_init_pl_theta") - l_init_pl_rh=options_get_logical(current_state%options_database, "l_init_pl_rh") + l_init_pl_rh=options_get_logical(current_state%options_database, "l_init_pl_rh") l_init_pl_q=options_get_logical(current_state%options_database, "l_init_pl_q") if (l_init_pl_q) then allocate(names_init_pl_q(options_get_array_size(current_state%options_database, "names_init_pl_q"))) call options_get_string_array(current_state%options_database, "names_init_pl_q", names_init_pl_q) do n = 1,size(names_init_pl_q) - if (trim(names_init_pl_q(n)) .eq. 'vapour' .and. l_init_pl_rh) then + if (trim(names_init_pl_q(n)) .eq. 'vapour' .and. l_init_pl_rh) then call log_master_log(LOG_ERROR, "Initialisation of vapour and RH - STOP") endif enddo @@ -262,7 +264,7 @@ subroutine calculate_initial_profiles(current_state, vertical_grid) do i=current_state%local_grid%local_domain_start_index(X_INDEX), current_state%local_grid%local_domain_end_index(X_INDEX) do j=current_state%local_grid%local_domain_start_index(Y_INDEX), current_state%local_grid%local_domain_end_index(Y_INDEX) current_state%th%data(:,j,i) = current_state%global_grid%configuration%vertical%theta_init(:) - & - current_state%global_grid%configuration%vertical%thref(:) + current_state%global_grid%configuration%vertical%thref(:) end do end do end if @@ -338,7 +340,7 @@ subroutine calculate_initial_profiles(current_state, vertical_grid) end do deallocate(f_init_pl_q_tmp, z_init_pl_q, f_init_pl_q, names_init_pl_q) end if - deallocate(zgrid) + deallocate(zgrid) end subroutine calculate_initial_profiles !> Calculates the mixing length for the neutral case @@ -356,8 +358,8 @@ subroutine calculate_mixing_length_for_neutral_case(current_state, vertical_grid vertical_grid%rneutml(k)=sqrt(1.0_DEFAULT_PRECISION/(1.0_DEFAULT_PRECISION/(von_karman_constant*& (vertical_grid%z(k)+z0))**2+1.0_DEFAULT_PRECISION/current_state%rmlmax**2) ) vertical_grid%rneutml_sq(k)=vertical_grid%rneutml(k)*vertical_grid%rneutml(k) - end do - end subroutine calculate_mixing_length_for_neutral_case + end do + end subroutine calculate_mixing_length_for_neutral_case !> Sets the buoyancy coefficient from the grid configuration and configuration !! @param current_state The current model state_mod @@ -368,15 +370,15 @@ subroutine set_buoyancy_coefficient(current_state, vertical_grid, kkp) type(vertical_grid_configuration_type), intent(inout) :: vertical_grid integer, intent(in) :: kkp - integer :: k + integer :: k if (.not. current_state%passive_th) then - if(current_state%use_anelastic_equations)then - do k=1, kkp-1 + if(current_state%use_anelastic_equations)then + do k=1, kkp-1 vertical_grid%buoy_co(k)=cp*(vertical_grid%prefn(k)**r_over_cp-vertical_grid%prefn(k+1)**r_over_cp)/& ((current_state%surface_reference_pressure**r_over_cp)*vertical_grid%dzn(k+1)) end do - else + else vertical_grid%buoy_co(1:kkp-1)=G/current_state%thref0 ! _Boussinesq end if ! Dummy value at top level @@ -397,14 +399,15 @@ subroutine setup_reference_state_liquid_water_temperature_and_saturation(current real(kind=DEFAULT_PRECISION) :: delta_t=1.0_DEFAULT_PRECISION, qlinit, tinit, qsatin, dqsatdtin, dsatfacin integer :: iter, k - + do k=1,kkp vertical_grid%tref(k)=vertical_grid%thref(k)*(vertical_grid%prefn(k)/current_state%surface_reference_pressure)**r_over_cp - vertical_grid%tstarpr(k)=0.0_DEFAULT_PRECISION - end do + vertical_grid%tstarpr(k)=0.0_DEFAULT_PRECISION + end do + if (current_state%th%active) then ! PREFRCP is used and hence calculated if theta is active - do k = 1,kkp + do k = 1,kkp vertical_grid%prefrcp(k)=(current_state%surface_reference_pressure/vertical_grid%prefn(k))**r_over_cp vertical_grid%rprefrcp(k)=1.0_DEFAULT_PRECISION/vertical_grid%prefrcp(k) ! Denotion between setup run and chain run in LEM - need to consider here too @@ -438,12 +441,12 @@ subroutine setup_reference_state_liquid_water_temperature_and_saturation(current ! Denotion between setup run and chain run in LEM - need to consider here too vertical_grid%tstarpr(k)= tinit-vertical_grid%tref(k) vertical_grid%qsat(k)=qsatin - vertical_grid%dqsatdt(k)=dqsatdtin + vertical_grid%dqsatdt(k)=dqsatdtin vertical_grid%qsatfac(k)= ( 1.0_DEFAULT_PRECISION/ ( 1.0_DEFAULT_PRECISION + rlvap_over_cp*dqsatdtin ) ) end do endif end if - end subroutine setup_reference_state_liquid_water_temperature_and_saturation + end subroutine setup_reference_state_liquid_water_temperature_and_saturation !> Sets up the reference properties for the vertical grid at each point !! @param current_state The current model state_mod @@ -453,12 +456,12 @@ subroutine set_up_vertical_reference_properties(current_state, vertical_grid, kk type(model_state_type), intent(inout) :: current_state type(vertical_grid_configuration_type), intent(inout) :: vertical_grid integer, intent(in) :: kkp - + integer :: k do k=1,kkp ! vertical_grid%thref(k)=current_state%thref0 -! vertical_grid%theta_init(k)=vertical_grid%thref(k) ! In LEM this can also be set from configuration (TODO) +! vertical_grid%theta_init(k)=vertical_grid%thref(k) ! In LEM this can also be set from configuration (TODO) vertical_grid%prefn(k)=0.0_DEFAULT_PRECISION vertical_grid%pdiff(k)=0.0_DEFAULT_PRECISION vertical_grid%rho(k)=current_state%rhobous @@ -468,7 +471,7 @@ subroutine set_up_vertical_reference_properties(current_state, vertical_grid, kk vertical_grid%dthref(k)=vertical_grid%thref(k+1)-vertical_grid%thref(k) end do vertical_grid%dthref(kkp)=0.0_DEFAULT_PRECISION - end subroutine set_up_vertical_reference_properties + end subroutine set_up_vertical_reference_properties !> Sets up and smooths the vertical grid. This is based upon the grid configuration already read in !! @param vertical_grid The vertical grid @@ -495,13 +498,13 @@ subroutine set_up_and_smooth_grid(vertical_grid, kgd, hgd, ninitp, kkp, zztop, n call new_vertical_grid_setup(vertical_grid, kgd, kkp, zztop) end if end if - + ! Regardless of the vertical grid computation method, set the level deltas do k=2,kkp vertical_grid%dz(k)=vertical_grid%z(k)-vertical_grid%z(k-1) - vertical_grid%dzn(k)= vertical_grid%zn(k)-vertical_grid%zn(k-1) - vertical_grid%rdz(k)=1./vertical_grid%dz(k) - vertical_grid%rdzn(k)=1./vertical_grid%dzn(k) + vertical_grid%dzn(k)= vertical_grid%zn(k)-vertical_grid%zn(k-1) + vertical_grid%rdz(k)=1./vertical_grid%dz(k) + vertical_grid%rdzn(k)=1./vertical_grid%dzn(k) end do vertical_grid%dzn(1)=0.d0 end subroutine set_up_and_smooth_grid @@ -521,7 +524,7 @@ subroutine original_vertical_grid_setup(vertical_grid, kgd, hgd, ninitp, kkp, zz real(kind=DEFAULT_PRECISION), intent(in) :: zztop integer :: n, k - + call create_linear_grid(vertical_grid, kgd, hgd, ninitp, kkp, zztop) ! Smooth grid vertical_grid%z(1)=0.0_DEFAULT_PRECISION @@ -585,8 +588,8 @@ subroutine new_vertical_grid_setup(vertical_grid, kgd, kkp, zztop) end if a(k)=a0+real(k-k0, kind=DEFAULT_PRECISION)*d0 end if - end do - + end do + do k=1, kkp vertical_grid%z(k)=a(k*2) vertical_grid%zn(k)=a(2*k-1) @@ -625,7 +628,7 @@ subroutine create_linear_grid(vertical_grid, kgd, hgd, ninitp, kkp, zztop) if(kgd(i) .gt. 0) then kmax=kgd(i) zmax=hgd(i) - + do k=kgd(i-1)+1,kgd(i) vertical_grid%z(k)=hgd(i-1)+(hgd(i)-hgd(i-1))*real(k-kgd(i-1), kind=DEFAULT_PRECISION)& /real(kgd(i)-kgd(i-1), kind=DEFAULT_PRECISION) @@ -667,7 +670,7 @@ subroutine allocate_vertical_grid_data(vertical_grid, n, nq) if (.not. allocated(vertical_grid%zn)) allocate(vertical_grid%zn(n)) allocate(vertical_grid%q_rand(n,nq), vertical_grid%q_init(n,nq), vertical_grid%q_force(n,nq)) - end subroutine allocate_vertical_grid_data + end subroutine allocate_vertical_grid_data !> Initialises the horizontal grid configurations !! @param current_state The current model state_mod @@ -689,7 +692,7 @@ subroutine initialise_horizontalgrid_configuration_types(current_state) 0.25_DEFAULT_PRECISION/current_state%global_grid%configuration%horizontal%dx current_state%global_grid%configuration%horizontal%tcy=& 0.25_DEFAULT_PRECISION/current_state%global_grid%configuration%horizontal%dy - end subroutine initialise_horizontalgrid_configuration_types + end subroutine initialise_horizontalgrid_configuration_types !> Set reference profile of potential temperature for the Boussinesq/Anelastic approximation !! Note that this is not in general the same as the profile defining the initial vertical distribution of potential @@ -698,7 +701,7 @@ end subroutine initialise_horizontalgrid_configuration_types !! @param current_state The current model state subroutine set_anelastic_pressure(current_state) type(model_state_type), intent(inout) :: current_state - + if (current_state%use_anelastic_equations) then call compute_anelastic_pressure_profile_and_density(current_state) else @@ -711,7 +714,8 @@ subroutine set_anelastic_pressure(current_state) current_state%global_grid%configuration%vertical%rhon=current_state%rhobous current_state%global_grid%configuration%vertical%pdiff=0.0_DEFAULT_PRECISION end if - end subroutine set_anelastic_pressure + + end subroutine set_anelastic_pressure !> Computes the anelastic pressure only - if we are using Boussinesq approximation !! @param current_state The current model state @@ -723,7 +727,7 @@ subroutine compute_anelastic_pressure_profile_only(current_state) , ptop &!pressure at z=ZN(KKP) , thprof(current_state%local_grid%size(Z_INDEX)) - + ! TODO: NOTE - we are mocking in thprof at the moment, this should be read from a configuration and used here instead thprof=0.0_DEFAULT_PRECISION ptop=0.0_DEFAULT_PRECISION @@ -749,7 +753,7 @@ subroutine compute_anelastic_pressure_profile_only(current_state) end do p0=0.5_DEFAULT_PRECISION*(current_state%global_grid%configuration%vertical%prefn(1)+& current_state%global_grid%configuration%vertical%prefn(2)) - if (ipass .eq. 1) then + if (ipass .eq. 1) then ptop=current_state%surface_pressure**r_over_cp+ptop**r_over_cp-p0**r_over_cp if (ptop .le. 0.0_DEFAULT_PRECISION .and. current_state%parallel%my_rank==0) then call log_log(LOG_ERROR, "Negative ptop in setup of anelastic. Need a warmer THREF or different setup options") @@ -757,7 +761,7 @@ subroutine compute_anelastic_pressure_profile_only(current_state) ptop=ptop**(1.0_DEFAULT_PRECISION/r_over_cp) end if end do - end subroutine compute_anelastic_pressure_profile_only + end subroutine compute_anelastic_pressure_profile_only !> Computes the anelastic pressure and density - if we are using Anelastic approximation !! @param current_state The current model state @@ -772,14 +776,14 @@ subroutine compute_anelastic_pressure_profile_and_density(current_state) ptop=0.0_DEFAULT_PRECISION current_state%global_grid%configuration%vertical%pdiff(current_state%local_grid%size(Z_INDEX))=0.0_DEFAULT_PRECISION do ipass=1,2 ! _after first pass, may adjust basic states - if (ipass .eq. 1 .or. ANELASTIC_PROFILE_MODE .gt. 1) then + if (ipass .eq. 1 .or. ANELASTIC_PROFILE_MODE .gt. 1) then current_state%global_grid%configuration%vertical%prefn(current_state%local_grid%size(Z_INDEX))=& (ptop/current_state%surface_reference_pressure)**r_over_cp do k=current_state%local_grid%size(Z_INDEX)-1,1,-1 current_state%global_grid%configuration%vertical%pdiff(k)=G*& current_state%global_grid%configuration%vertical%dzn(k+1)/(0.5_DEFAULT_PRECISION*cp*& (current_state%global_grid%configuration%vertical%thref(k)+& - current_state%global_grid%configuration%vertical%thref(k+1))) + current_state%global_grid%configuration%vertical%thref(k+1))) end do do k=current_state%local_grid%size(Z_INDEX)-1,1,-1 current_state%global_grid%configuration%vertical%prefn(k)=& @@ -789,8 +793,8 @@ subroutine compute_anelastic_pressure_profile_and_density(current_state) do k=current_state%local_grid%size(Z_INDEX),1,-1 current_state%global_grid%configuration%vertical%prefn(k)=current_state%surface_reference_pressure*& current_state%global_grid%configuration%vertical%prefn(k)**(1.0_DEFAULT_PRECISION/r_over_cp) - end do - end if + end do + end if p0=0.5_DEFAULT_PRECISION*(current_state%global_grid%configuration%vertical%prefn(1)+& current_state%global_grid%configuration%vertical%prefn(2)) ! !------------------------------------------------------------- @@ -802,7 +806,7 @@ subroutine compute_anelastic_pressure_profile_and_density(current_state) ! ! _Option 3 tends to give rather large changes in PSF, so ! ! I prefer 2 or 4 for most purposes ! !------------------------------------------------------------- - if (ipass .eq. 1 .and. ANELASTIC_PROFILE_MODE .eq. 2) then + if (ipass .eq. 1 .and. ANELASTIC_PROFILE_MODE .eq. 2) then ! ! _adjust THREF profile by constant factor to enforce ! ! P0 = (fixed) PSF thfactor=((p0/current_state%surface_reference_pressure)**r_over_cp-(ptop/current_state%surface_reference_pressure)**& @@ -813,14 +817,14 @@ subroutine compute_anelastic_pressure_profile_and_density(current_state) current_state%global_grid%configuration%vertical%thref(k)*thfactor end do end if - if (ipass .eq. 1 .and. ANELASTIC_PROFILE_MODE .eq. 4) then + if (ipass .eq. 1 .and. ANELASTIC_PROFILE_MODE .eq. 4) then ! ! _adjust PTOP so that P0 = (fixed) PSF ptop=current_state%surface_pressure**r_over_cp+ptop**r_over_cp-p0**r_over_cp if (ptop .le. 0.0_DEFAULT_PRECISION .and. current_state%parallel%my_rank==0) then call log_log(LOG_ERROR, "Negative ptop in setup of anelastic. Need a warmer THREF or different setup options") end if ptop=ptop**(1.0_DEFAULT_PRECISION/r_over_cp) - end if + end if end do ! !--------------------------------------- ! ! _Finally compute density from pressure @@ -832,12 +836,12 @@ subroutine compute_anelastic_pressure_profile_and_density(current_state) end do do k=1,current_state%local_grid%size(Z_INDEX)-1 current_state%global_grid%configuration%vertical%rho(k)=sqrt(current_state%global_grid%configuration%vertical%rhon(k)*& - current_state%global_grid%configuration%vertical%rhon(k+1)) + current_state%global_grid%configuration%vertical%rhon(k+1)) end do current_state%global_grid%configuration%vertical%rho(current_state%local_grid%size(Z_INDEX))=& current_state%global_grid%configuration%vertical%rhon(current_state%local_grid%size(Z_INDEX))**2/& - current_state%global_grid%configuration%vertical%rhon(current_state%local_grid%size(Z_INDEX)-1) - end subroutine compute_anelastic_pressure_profile_and_density + current_state%global_grid%configuration%vertical%rhon(current_state%local_grid%size(Z_INDEX)-1) + end subroutine compute_anelastic_pressure_profile_and_density subroutine check_top(zztop, z, info) @@ -862,7 +866,7 @@ subroutine check_input_levels(z_levels, field_levels, field) end if end subroutine check_input_levels - + subroutine set_qv_init_from_rh(current_state) type(model_state_type), intent(inout) :: current_state @@ -882,10 +886,10 @@ subroutine set_qv_init_from_rh(current_state) vertical_grid=current_state%global_grid%configuration%vertical allocate(zgrid(current_state%local_grid%local_domain_end_index(Z_INDEX))) - + zztop = current_state%global_grid%top(Z_INDEX) - l_init_pl_rh=options_get_logical(current_state%options_database, "l_init_pl_rh") + l_init_pl_rh=options_get_logical(current_state%options_database, "l_init_pl_rh") if (l_init_pl_rh)then allocate(z_init_pl_rh(options_get_array_size(current_state%options_database, "z_init_pl_rh")), & @@ -897,7 +901,7 @@ subroutine set_qv_init_from_rh(current_state) zgrid=current_state%global_grid%configuration%vertical%zn(:) call piecewise_linear_1d(z_init_pl_rh(1:size(z_init_pl_rh)), f_init_pl_rh(1:size(f_init_pl_rh)), zgrid, & current_state%global_grid%configuration%vertical%rh_init) - + if (.not. current_state%passive_q .and. current_state%th%active) then iq=get_q_index('vapour', 'piecewise_initialization') allocate(TdegK(current_state%local_grid%local_domain_end_index(Z_INDEX))) @@ -905,8 +909,8 @@ subroutine set_qv_init_from_rh(current_state) (vertical_grid%prefn(:)/current_state%surface_reference_pressure)**r_over_cp do k = current_state%local_grid%local_domain_start_index(Z_INDEX), & current_state%local_grid%local_domain_end_index(Z_INDEX) - qsat=qsaturation(TdegK(k), current_state%global_grid%configuration%vertical%prefn(k)/100.) - current_state%global_grid%configuration%vertical%q_init(k, iq) = & + qsat=qsaturation(TdegK(k), current_state%global_grid%configuration%vertical%prefn(k)/100.) + current_state%global_grid%configuration%vertical%q_init(k, iq) = & (current_state%global_grid%configuration%vertical%rh_init(k)/100.0)*qsat !print *, current_state%global_grid%configuration%vertical%rh_init(k), & ! current_state%global_grid%configuration%vertical%q_init(k, iq), & @@ -926,7 +930,7 @@ subroutine set_qv_init_from_rh(current_state) else call log_master_log(LOG_ERROR, "Initialising with RH but q and/or theta passive") end if - + deallocate(z_init_pl_rh, f_init_pl_rh) end if diff --git a/components/iobridge/makefile b/components/iobridge/makefile index d4a31254..90deb259 100644 --- a/components/iobridge/makefile +++ b/components/iobridge/makefile @@ -3,7 +3,9 @@ SRCSF = src/iobridge.F90 BUILDDIR=build COREDIR=../../model_core/build IODIR=../../io/build -FFLAGS=-I $(BUILDDIR) -I $(COREDIR) -I $(IODIR) $(COMPILERFFLAGS) +COLUMNDIR=../conditional_diagnostics_column/build + +FFLAGS=-I $(BUILDDIR) -I $(COREDIR) -I $(IODIR) -I $(COLUMNDIR) $(COMPILERFFLAGS) OBJS = $(patsubst %.F90,$(BUILDDIR)/%.o,$(SRCSF)) all: create-build-dirs $(OBJS) diff --git a/components/lowerbc/src/lowerbc.F90 b/components/lowerbc/src/lowerbc.F90 index 9d45e08d..1e08d8f6 100644 --- a/components/lowerbc/src/lowerbc.F90 +++ b/components/lowerbc/src/lowerbc.F90 @@ -20,7 +20,7 @@ module lowerbc_mod #endif integer, parameter :: CONVERGENCE_SUCCESS=1, CONVERGENCE_RICHARDSON_TOO_LARGE=2, CONVERGENCE_FAILURE=3 - + real(kind=DEFAULT_PRECISION), parameter :: smth = 0.05_DEFAULT_PRECISION,& ! Smoothing between iterations tolm=1.0E-4_DEFAULT_PRECISION, tolt=1.0E-4_DEFAULT_PRECISION ! Convergence tollerance for u and t star @@ -34,6 +34,9 @@ module lowerbc_mod integer :: wrapping_comm_requests(4), y_wrapping_target_id, x_wrapping_target_id public lowerbc_get_descriptor + public tstrcona, rhmbc, ddbc, ddbc_x4, eecon, r2ddbc, rcmbc, tstrconb, & + x4con, xx0con, y2con, yy0con, viscous_courant_coefficient + contains !> Descriptor of this component for registration @@ -45,7 +48,7 @@ type(component_descriptor_type) function lowerbc_get_descriptor() lowerbc_get_descriptor%timestep=>timestep_callback lowerbc_get_descriptor%finalisation=>finalisation_callback end function lowerbc_get_descriptor - + subroutine initialisation_callback(current_state) type(model_state_type), target, intent(inout) :: current_state @@ -56,7 +59,7 @@ subroutine initialisation_callback(current_state) ! are allocated in their respective components if (.not. is_component_enabled(current_state%options_database, "diffusion")) then call log_master_log(LOG_ERROR, "Lowerbc requires the diffusion component to be enabled") - end if + end if if (.not. is_component_enabled(current_state%options_database, "viscosity")) then call log_master_log(LOG_ERROR, "Lowerbc requires the viscosity component to be enabled") end if @@ -71,7 +74,7 @@ subroutine initialisation_callback(current_state) if (num_wrapped_fields .gt. 0) then if (current_state%parallel%my_coords(Y_INDEX) == 0 .or. & - current_state%parallel%my_coords(Y_INDEX) == current_state%parallel%dim_sizes(Y_INDEX)-1) then + current_state%parallel%my_coords(Y_INDEX) == current_state%parallel%dim_sizes(Y_INDEX)-1) then if (current_state%parallel%my_coords(Y_INDEX) == 0) then y_wrapping_target_id=current_state%local_grid%neighbours(Y_INDEX, 1) else @@ -84,7 +87,7 @@ subroutine initialisation_callback(current_state) end if if (current_state%parallel%my_coords(X_INDEX) == 0 .or. & - current_state%parallel%my_coords(X_INDEX) == current_state%parallel%dim_sizes(X_INDEX)-1) then + current_state%parallel%my_coords(X_INDEX) == current_state%parallel%dim_sizes(X_INDEX)-1) then if (current_state%parallel%my_coords(X_INDEX) == 0) then x_wrapping_target_id=current_state%local_grid%neighbours(X_INDEX, 1) else @@ -102,10 +105,10 @@ subroutine initialisation_callback(current_state) current_state%global_grid%configuration%horizontal%dx)+& 1.0_DEFAULT_PRECISION/(current_state%global_grid%configuration%horizontal%dy*& current_state%global_grid%configuration%horizontal%dy) - + if ( current_state%use_surface_boundary_conditions .and. & current_state%type_of_surface_boundary_conditions == PRESCRIBED_SURFACE_VALUES) then - ! variables below are only required when PRESCRIBED_SURFACE_VALUES are used. + ! variables below are only required when PRESCRIBED_SURFACE_VALUES are used. tstrcona=von_karman_constant/alphah*current_state%global_grid%configuration%vertical%zlogth bhbc=alphah*current_state%global_grid%configuration%vertical%zlogth rhmbc=betah*(current_state%global_grid%configuration%vertical%zn(2)+z0-z0th)/& @@ -135,10 +138,10 @@ subroutine initialisation_callback(current_state) xx0con=0.0 y2con=0.0 yy0con=0.0 - endif - + endif + ! Determine vapour index - if (.not. current_state%passive_q) then + if (.not. current_state%passive_q) then iqv = get_q_index(standard_q_names%VAPOUR, 'lowerbc') endif @@ -167,12 +170,12 @@ subroutine allocate_applicable_fields(current_state) do i=1,current_state%number_q_fields allocate(current_state%disq(i)%data(z_size, y_size, x_size)) - end do + end do end subroutine allocate_applicable_fields subroutine timestep_callback(current_state) type(model_state_type), target, intent(inout) :: current_state - + integer :: current_y_index, current_x_index current_y_index=current_state%column_local_y @@ -189,8 +192,8 @@ subroutine timestep_callback(current_state) call compute_lower_boundary_conditions(current_state, current_y_index, current_x_index, & current_state%zu, current_state%zv, current_state%zth, current_state%zth, current_state%zq, current_state%zq) end if - end if - end subroutine timestep_callback + end if + end subroutine timestep_callback subroutine compute_lower_boundary_conditions(current_state, current_y_index, current_x_index, zu, zv, zth, th, zq, q) type(model_state_type), target, intent(inout) :: current_state @@ -211,7 +214,7 @@ subroutine compute_lower_boundary_conditions(current_state, current_y_index, cur current_state%column_local_x .gt. current_state%local_grid%local_domain_end_index(X_INDEX)+1)) then !if (.not. current_state%halo_column) then - ! Include one halo to ensure that the halo is set in tvdadvection. This is done using the + ! Include one halo to ensure that the halo is set in tvdadvection. This is done using the ! logic from the timestep callback in tvdadvection in the timestep callback above horizontal_velocity_at_k2=0.0_DEFAULT_PRECISION #ifdef U_ACTIVE @@ -222,7 +225,7 @@ subroutine compute_lower_boundary_conditions(current_state, current_y_index, cur horizontal_velocity_at_k2=horizontal_velocity_at_k2+(0.5_DEFAULT_PRECISION*(zv%data(& 2,current_y_index,current_x_index)+zv%data(2,current_y_index-1,current_x_index))+current_state%vgal)**2 #endif - horizontal_velocity_at_k2=sqrt(horizontal_velocity_at_k2)+smallp + horizontal_velocity_at_k2=sqrt(horizontal_velocity_at_k2)+smallp if (current_state%type_of_surface_boundary_conditions == PRESCRIBED_SURFACE_FLUXES) then call compute_using_fixed_surface_fluxes(current_state, current_y_index, current_x_index, & @@ -240,7 +243,7 @@ subroutine compute_lower_boundary_conditions(current_state, current_y_index, cur current_state%disq(n)%data(1,current_y_index,current_x_index)=0.0_DEFAULT_PRECISION end do end if - + !----------------------- ! _return viscous number !----------------------- @@ -274,7 +277,7 @@ subroutine register_async_wrapping_recv_requests(current_state) x_wrapping_target_id, 0, current_state%parallel%neighbour_comm, wrapping_comm_requests(3), ierr) end if end subroutine register_async_wrapping_recv_requests - + !> Completes the asynchronous wrapping if required for periodic boundary conditions !! @param current_state The current model state !! @param zth Temperature field @@ -290,26 +293,26 @@ subroutine complete_async_wrapping(current_state, zth, zq) if (allocated(y_wrapping_send_buffer)) then if (current_state%parallel%my_coords(Y_INDEX) == 0) then call package_y_wrapping_send_buffer(current_state, zth, zq, current_state%local_grid%local_domain_start_index(Y_INDEX),& - current_state%local_grid%local_domain_start_index(Y_INDEX)+1) + current_state%local_grid%local_domain_start_index(Y_INDEX)+1) else call package_y_wrapping_send_buffer(current_state, zth, zq, current_state%local_grid%local_domain_end_index(Y_INDEX)-1,& - current_state%local_grid%local_domain_end_index(Y_INDEX)) - end if + current_state%local_grid%local_domain_end_index(Y_INDEX)) + end if call mpi_isend(y_wrapping_send_buffer, size(y_wrapping_send_buffer), PRECISION_TYPE, & y_wrapping_target_id, 0, current_state%parallel%neighbour_comm, & - wrapping_comm_requests(2), ierr) + wrapping_comm_requests(2), ierr) end if if (allocated(x_wrapping_send_buffer)) then if (current_state%parallel%my_coords(X_INDEX) == 0) then call package_x_wrapping_send_buffer(current_state, zth, zq, current_state%local_grid%local_domain_start_index(X_INDEX),& - current_state%local_grid%local_domain_start_index(X_INDEX)+1) + current_state%local_grid%local_domain_start_index(X_INDEX)+1) else call package_x_wrapping_send_buffer(current_state, zth, zq, current_state%local_grid%local_domain_end_index(X_INDEX)-1,& - current_state%local_grid%local_domain_end_index(X_INDEX)) - end if + current_state%local_grid%local_domain_end_index(X_INDEX)) + end if call mpi_isend(x_wrapping_send_buffer, size(x_wrapping_send_buffer), PRECISION_TYPE, & x_wrapping_target_id, 0, current_state%parallel%neighbour_comm, & - wrapping_comm_requests(4), ierr) + wrapping_comm_requests(4), ierr) end if ! If send buffer is allocated then recv buffer is allocated, therefore only test the send buffer here and assume recv @@ -321,7 +324,7 @@ subroutine complete_async_wrapping(current_state, zth, zq) else call unpackage_y_wrapping_recv_buffer(current_state, zth, zq, & current_state%local_grid%local_domain_end_index(Y_INDEX)+1, & - current_state%local_grid%local_domain_end_index(Y_INDEX)+2) + current_state%local_grid%local_domain_end_index(Y_INDEX)+2) end if end if if (allocated(x_wrapping_recv_buffer)) then @@ -330,10 +333,10 @@ subroutine complete_async_wrapping(current_state, zth, zq) else call unpackage_x_wrapping_recv_buffer(current_state, zth, zq, & current_state%local_grid%local_domain_end_index(X_INDEX)+1, & - current_state%local_grid%local_domain_end_index(X_INDEX)+2) + current_state%local_grid%local_domain_end_index(X_INDEX)+2) end if end if - end if + end if if (current_state%parallel%my_rank == y_wrapping_target_id) then if (current_state%th%active) then zth%data(1,1,:)=zth%data(1, current_state%local_grid%local_domain_end_index(Y_INDEX)-1, :) @@ -418,7 +421,7 @@ subroutine package_x_wrapping_send_buffer(current_state, zth, zq, first_x_index, integer, intent(in) :: first_x_index, second_x_index integer :: index_start, n - + index_start=0 if (current_state%th%active) then x_wrapping_send_buffer(:,1,1)=zth%data(1,:,first_x_index) @@ -510,8 +513,8 @@ subroutine handle_convective_fluxes(current_state, current_y_index, current_x_in /ustr**3))/ (1.+sqrt(1.+gammah*von_karman_constant*current_state%fbuoy*z0th/ustr**3)))) if (current_state%th%active) th%data(1, current_y_index, current_x_index)= & (current_state%surface_temperature_flux*current_state%global_grid%configuration%vertical%dzn(2)/& - current_state%diff_coefficient%data(1, current_y_index, current_x_index))+th%data(2, current_y_index, current_x_index)-& - current_state%global_grid%configuration%vertical%thref(1)+& + current_state%diff_coefficient%data(1, current_y_index, current_x_index))+th%data(2, current_y_index, current_x_index)-& + current_state%global_grid%configuration%vertical%thref(1)+& current_state%global_grid%configuration%vertical%thref(2) ! Surface Flux of vapour @@ -532,16 +535,16 @@ real(kind=DEFAULT_PRECISION) function look(current_state, vel) lookup_real_posn=1.0_DEFAULT_PRECISION+real(current_state%lookup_table_entries-1)*& log(vel/current_state%velmin)*current_state%aloginv - lookup_int_posn=int(lookup_real_posn) + lookup_int_posn=int(lookup_real_posn) - if (lookup_int_posn .ge. 1) then + if (lookup_int_posn .ge. 1) then if (lookup_int_posn .lt. current_state%lookup_table_entries) then ! Linear interpolation look=current_state%lookup_table_ustr(lookup_int_posn)+ (lookup_real_posn-real(lookup_int_posn))*& (current_state%lookup_table_ustr(lookup_int_posn+1)-current_state%lookup_table_ustr(lookup_int_posn)) else ! Near neutral look=vel*current_state%cneut end if - else ! Nearly free convection + else ! Nearly free convection look=vel**(-convective_limit)*current_state%cfc end if end function look @@ -557,7 +560,7 @@ subroutine handle_neutral_fluxes(current_state, current_y_index, current_x_index integer :: n real(kind=DEFAULT_PRECISION) :: ustr - + ustr=horizontal_velocity_at_k2*current_state%global_grid%configuration%vertical%vk_on_zlogm current_state%vis_coefficient%data(1, current_y_index, current_x_index)=current_state%global_grid%configuration%vertical%czn*& ustr**2/horizontal_velocity_at_k2 @@ -566,8 +569,8 @@ subroutine handle_neutral_fluxes(current_state, current_y_index, current_x_index current_state%global_grid%configuration%vertical%zlogm/(alphah*current_state%global_grid%configuration%vertical%zlogth) if (current_state%th%active) th%data(1, current_y_index, current_x_index)= & (current_state%surface_temperature_flux*current_state%global_grid%configuration%vertical%dzn(2)/& - current_state%diff_coefficient%data(1, current_y_index, current_x_index))+th%data(2, current_y_index, current_x_index)-& - current_state%global_grid%configuration%vertical%thref(1)+& + current_state%diff_coefficient%data(1, current_y_index, current_x_index))+th%data(2, current_y_index, current_x_index)-& + current_state%global_grid%configuration%vertical%thref(1)+& current_state%global_grid%configuration%vertical%thref(2) ! Flux of vapour @@ -617,8 +620,8 @@ subroutine handle_stable_fluxes(current_state, current_y_index, current_x_index, von_karman_constant*current_state%fbuoy* (current_state%global_grid%configuration%vertical%zn(2)+ z0-z0th)/ustr**3) if (current_state%th%active) th%data(1, current_y_index, current_x_index)= & (current_state%surface_temperature_flux*current_state%global_grid%configuration%vertical%dzn(2)/& - current_state%diff_coefficient%data(1, current_y_index, current_x_index))+th%data(2, current_y_index, current_x_index)-& - current_state%global_grid%configuration%vertical%thref(1)+& + current_state%diff_coefficient%data(1, current_y_index, current_x_index))+th%data(2, current_y_index, current_x_index)-& + current_state%global_grid%configuration%vertical%thref(1)+& current_state%global_grid%configuration%vertical%thref(2) @@ -648,7 +651,7 @@ subroutine compute_using_fixed_surface_fluxes(current_state, current_y_index, cu call handle_stable_fluxes(current_state, current_y_index, current_x_index, horizontal_velocity_at_k2, th, q) end if end subroutine compute_using_fixed_surface_fluxes - + subroutine compute_using_fixed_surface_temperature(current_state, current_y_index, current_x_index, horizontal_velocity_at_k2, & zth, th, zq, q) @@ -661,19 +664,19 @@ subroutine compute_using_fixed_surface_temperature(current_state, current_y_inde real(kind=DEFAULT_PRECISION) :: dthv_surf, ustr, thvstr integer :: convergence_status, n - if (current_state%passive_q) then ! i.e. q is not active + if (current_state%passive_q) then ! i.e. q is not active ! Assuming no liquid water at level 2 dthv_surf = zth%data(2, current_y_index, current_x_index) + & current_state%global_grid%configuration%vertical%thref(2) - current_state%theta_virtual_surf - else + else dthv_surf=zth%data(2, current_y_index, current_x_index) + current_state%global_grid%configuration%vertical%thref(2)& *(1.0_DEFAULT_PRECISION + current_state%cq(current_state%water_vapour_mixing_ratio_index)*& zq(current_state%water_vapour_mixing_ratio_index)%data(2,current_y_index,current_x_index)) - & - current_state%theta_virtual_surf + current_state%theta_virtual_surf end if convergence_status = mostbc(current_state, horizontal_velocity_at_k2, dthv_surf,& current_state%global_grid%configuration%vertical%zn(2), ustr, thvstr) - + current_state%vis_coefficient%data(1, current_y_index, current_x_index)=& current_state%global_grid%configuration%vertical%czn*ustr**2/horizontal_velocity_at_k2 current_state%diff_coefficient%data(1, current_y_index, current_x_index)=& @@ -682,7 +685,7 @@ subroutine compute_using_fixed_surface_temperature(current_state, current_y_inde (current_state%global_grid%configuration%vertical%thref(2)+current_state%global_grid%configuration%vertical%thref(1)) th%data(1, current_y_index, current_x_index)=2.0*current_state%theta_surf-th%data(2, current_y_index, current_x_index)-& (current_state%global_grid%configuration%vertical%thref(2)+current_state%global_grid%configuration%vertical%thref(1)) - + if (current_state%number_q_fields .gt. 0) then n=iqv zq(n)%data(1, current_y_index, current_x_index)=zq(n)%data(2, current_y_index, current_x_index) @@ -695,8 +698,8 @@ subroutine compute_using_fixed_surface_temperature(current_state, current_y_inde 2.0_DEFAULT_PRECISION*current_state%surface_vapour_mixing_ratio-& q(current_state%water_vapour_mixing_ratio_index)%data(2,current_y_index,current_x_index) endif - end if - end subroutine compute_using_fixed_surface_temperature + end if + end subroutine compute_using_fixed_surface_temperature subroutine simple_boundary_values(current_state, current_y_index, current_x_index, th, q) type(model_state_type), target, intent(inout) :: current_state @@ -724,8 +727,8 @@ end subroutine simple_boundary_values !> Solves the Monin-Obukhov equations in the case of specified surface values of temperature and mixing ratio,combined !! into a specified value of virtual temperature. It is a modified version of the subroutine described in Bull and - !! Derbyshire (TDN197) based on the assumption that the similarity functions and roughness lengths for temperature and - !! mixing ratio are the same. In that case, all the original theory can be used if we replace temperature by virtual + !! Derbyshire (TDN197) based on the assumption that the similarity functions and roughness lengths for temperature and + !! mixing ratio are the same. In that case, all the original theory can be used if we replace temperature by virtual !! temperature. !! The form of the non-dimensionalised wind shear used is 1.0 + BETAM z/L for the stable case and !! (1.0 - GAMMAM z/L)**1/4 for the unstable case. The form of the non-dimensionalised temperature gradient used is @@ -743,7 +746,7 @@ integer function mostbc(current_state, delu, delt, z1, ustrdg, tstrdg) real(kind=DEFAULT_PRECISION), intent(out) :: ustrdg, tstrdg if (delt .lt. 0.0_DEFAULT_PRECISION) then - if (delu .le. smallp) then + if (delu .le. smallp) then ustrdg=0.0_DEFAULT_PRECISION tstrdg=tstrcona*delt else @@ -758,7 +761,7 @@ integer function mostbc(current_state, delu, delt, z1, ustrdg, tstrdg) else ! Trivial neutral case ustrdg=current_state%global_grid%configuration%vertical%vk_on_zlogm*delu - tstrdg=0.0_DEFAULT_PRECISION + tstrdg=0.0_DEFAULT_PRECISION mostbc=CONVERGENCE_SUCCESS end if @@ -768,7 +771,7 @@ integer function mostbc(current_state, delu, delt, z1, ustrdg, tstrdg) else if(mostbc .eq. CONVERGENCE_FAILURE) then call log_log(LOG_ERROR, "Convergence failure after 200 iterations") end if - end if + end if end function mostbc integer function solve_monin_obukhov_unstable_case(delu, delt, ellmocon, ustrdg, tstrdg, vertical_grid) @@ -780,7 +783,7 @@ integer function solve_monin_obukhov_unstable_case(delu, delt, ellmocon, ustrdg, real(kind=DEFAULT_PRECISION) :: ellmo, psim, psih, x4, xx, xx0, y2, yy, yy0, err_ustr, err_tstr, & ustrl, tstrl, & ! U and T star at start of iteration ustrpl, tstrpl ! U and T star at end of iteration - + ! First set initial values ustrl=vertical_grid%vk_on_zlogm*delu tstrl=tstrcona*delt @@ -790,11 +793,11 @@ integer function solve_monin_obukhov_unstable_case(delu, delt, ellmocon, ustrdg, ellmo=ustrl*ustrl*ellmocon/tstrl ! Test for possible square root of negative quantity - x4=1.0_DEFAULT_PRECISION-x4con/ellmo - if (x4 .lt. 0.0_DEFAULT_PRECISION) call log_log(LOG_ERROR, "Negative square root in x4") + x4=1.0_DEFAULT_PRECISION-x4con/ellmo + if (x4 .lt. 0.0_DEFAULT_PRECISION) call log_log(LOG_ERROR, "Negative square root in x4") xx=sqrt(sqrt(x4)) - xx0=sqrt(sqrt(1.0_DEFAULT_PRECISION-xx0con / ellmo)) + xx0=sqrt(sqrt(1.0_DEFAULT_PRECISION-xx0con / ellmo)) psim=2.*( log((xx+1.0_DEFAULT_PRECISION)/(xx0+1.0_DEFAULT_PRECISION))-atan(xx)+atan(xx0) )+& log((xx*xx+1.0_DEFAULT_PRECISION)/(xx0*xx0+1.0_DEFAULT_PRECISION)) ustrpl=von_karman_constant*delu/(vertical_grid%zlogm-psim) @@ -807,13 +810,13 @@ integer function solve_monin_obukhov_unstable_case(delu, delt, ellmocon, ustrdg, psih=2.*log((1.0_DEFAULT_PRECISION+yy)/(1.0_DEFAULT_PRECISION+yy0)) tstrpl=tstrconb*delt/(vertical_grid%zlogth-psih) err_ustr=abs((ustrpl-ustrl)/ ustrl) - err_tstr=abs((tstrpl-tstrl)/ tstrl) - if ((err_tstr .lt. tolt) .and. (err_ustr .lt. tolm)) then + err_tstr=abs((tstrpl-tstrl)/ tstrl) + if ((err_tstr .lt. tolt) .and. (err_ustr .lt. tolm)) then ustrdg=ustrpl tstrdg=tstrpl solve_monin_obukhov_unstable_case=CONVERGENCE_SUCCESS return - else + else ustrl=(1.0_DEFAULT_PRECISION-smth)*ustrpl+smth*ustrl tstrl=(1.0_DEFAULT_PRECISION-smth)*tstrpl+smth*tstrl end if @@ -839,10 +842,10 @@ integer function solve_monin_obukhov_stable_case(delu, delt, zlogm, cmbc, ustrdg ustrdg=(-ee+sqrt(det))*r2ddbc tstrdg=ustrdg*(am-zlogm*ustrdg)*rcmbc else - solve_monin_obukhov_stable_case=CONVERGENCE_RICHARDSON_TOO_LARGE + solve_monin_obukhov_stable_case=CONVERGENCE_RICHARDSON_TOO_LARGE ustrdg=0.0_DEFAULT_PRECISION tstrdg=0.0_DEFAULT_PRECISION - end if + end if else if (ddbc .eq. 0.0_DEFAULT_PRECISION) then ! Degenerate case ustrdg=-ff/ee diff --git a/components/makefile b/components/makefile index 22f9d6ca..42cf4b51 100644 --- a/components/makefile +++ b/components/makefile @@ -1,5 +1,6 @@ # List of components to compile is the directories in this top level with makefiles in them -COMPONENTS=$(subst /,,$(sort $(dir $(wildcard */makefile)))) +COMPONENTS=buoyancy casim cfltest checkpointer clearsourceterms conditional_diagnostics_column conditional_diagnostics_whole coriolis damping debugger decomposition diffusion diverr fftsolver flux_budget forcing gridmanager haloswapper iobridge iterativesolver kidreader lowerbc lwrad_exponential meanprofiles modelsynopsis pdf_analysis petsc_solver pressuresource tvdadvection profile_diagnostics pstep pwadvection randomnoise scalar_diagnostics set_consistent_lowbc setfluxlook simplecloud simplesetup smagorinsky socrates_couple stepfields steppingdirection subgrid_profile_diagnostics swapsmooth terminationcheck thadvection viscosity dephy_forcings + COMPONENT_HEADER_FILE=componentheaders.autogen COMPONENT_REGISTRATIONS_FILE=componentregistrations.autogen diff --git a/components/pdf_analysis/src/pdf_analysis.F90 b/components/pdf_analysis/src/pdf_analysis.F90 index f8759879..5ad35cb7 100644 --- a/components/pdf_analysis/src/pdf_analysis.F90 +++ b/components/pdf_analysis/src/pdf_analysis.F90 @@ -1,4 +1,4 @@ -!> Calculates fields related to distributions of data on full-domain horizontal 2d slices +!> Calculates fields related to distributions of data on full-domain horizontal 2d slices module pdf_analysis_mod use monc_component_mod, only : COMPONENT_ARRAY_FIELD_TYPE, & COMPONENT_DOUBLE_DATA_TYPE, component_descriptor_type, & @@ -25,7 +25,7 @@ module pdf_analysis_mod ! these are available on all processes integer :: tpts ! total number of horizontal grid points on full domain - integer :: lpts ! local number of horizontal grid points on + integer :: lpts ! local number of horizontal grid points on logical :: show_critical_w ! stdout diagnostic logical @@ -78,7 +78,7 @@ subroutine init_callback(current_state) show_critical_w = options_get_logical(current_state%options_database, "show_critical_w") !> Allocate space for the global 2d field only on a single process -! if (current_state%parallel%my_rank == 0) +! if (current_state%parallel%my_rank == 0) allocate(tmp_all(tpts)) ! else ! allocate(tmp_all(1)) @@ -89,7 +89,7 @@ subroutine init_callback(current_state) call mpi_allgather(lpts, 1, MPI_INT, gpts_on_proc, 1, MPI_INT, current_state%parallel%monc_communicator, ierr) !> Allocate and initialize displacement values - allocate(displacements(current_state%parallel%processes)) + allocate(displacements(current_state%parallel%processes)) displacements(1) = 0 do inc = 2, current_state%parallel%processes displacements(inc) = displacements(inc-1) + gpts_on_proc(inc-1) @@ -113,7 +113,7 @@ subroutine timestep_callback(current_state) type(model_state_type), target, intent(inout) :: current_state !> Current forumulation only handles vertical velocity percentiles. - !! Future enhancements may employ this component to perform additional + !! Future enhancements may employ this component to perform additional !! operations that require access to full horizontal fields, such as !! pdf calculations. @@ -127,7 +127,7 @@ end subroutine timestep_callback subroutine finalisation_callback(current_state) type(model_state_type), target, intent(inout) :: current_state - if (allocated(tmp_all)) deallocate(tmp_all) + if (allocated(tmp_all)) deallocate(tmp_all) end subroutine finalisation_callback @@ -138,7 +138,7 @@ subroutine calculate_w_percentiles(current_state) real(kind=DEFAULT_PRECISION), dimension(lpts) :: tmp_var integer :: i, j, k, num_neg, num_pos, dd_thresh_pos, ud_thresh_pos - integer :: max_up_k, min_dwn_k + integer :: max_up_k, min_dwn_k real(kind=DEFAULT_PRECISION), dimension((lpts+1)/2) :: T real(kind=DEFAULT_PRECISION), dimension((tpts+1)/2) :: Tall real(kind=DEFAULT_PRECISION) :: max_up, min_dwn, & @@ -155,7 +155,7 @@ subroutine calculate_w_percentiles(current_state) min_dwn_k = 0 !> reset thresholds - current_state%global_grid%configuration%vertical%w_dwn(:) = 0.0_DEFAULT_PRECISION + current_state%global_grid%configuration%vertical%w_dwn(:) = 0.0_DEFAULT_PRECISION current_state%global_grid%configuration%vertical%w_up(:) = 0.0_DEFAULT_PRECISION !> Loop over levels @@ -185,7 +185,7 @@ subroutine calculate_w_percentiles(current_state) num_neg = count(tmp_all < 0.0_DEFAULT_PRECISION) num_pos = count(tmp_all > 0.0_DEFAULT_PRECISION) - dd_thresh_pos = int(num_neg * dwnpercrit) + dd_thresh_pos = int(num_neg * dwnpercrit) ud_thresh_pos = tpts - int(num_pos * uppercrit) + 1 if ( dd_thresh_pos == 0 ) dd_thresh_pos = 1 @@ -200,7 +200,7 @@ subroutine calculate_w_percentiles(current_state) min_dwn_th = tmp_all(dd_thresh_pos) min_dwn = tmp_all(1) ! sorted array min_dwn_k = k - end if + end if if ( tmp_all(ud_thresh_pos) > max_up_th ) then max_up_th = tmp_all(ud_thresh_pos) max_up = tmp_all(tpts) ! sorted array @@ -227,12 +227,12 @@ subroutine calculate_w_percentiles(current_state) call log_master_log(LOG_INFO, 'Maximum updraft: '& //trim(conv_to_string(max_up))//' at level '//trim(conv_to_string(max_up_k)) ) call log_master_log(LOG_INFO, 'Minimum downdraft threshold: '& - //trim(conv_to_string(min_dwn_th))//' found at level '//trim(conv_to_string(min_dwn_k)) ) + //trim(conv_to_string(min_dwn_th))//' found at level '//trim(conv_to_string(min_dwn_k)) ) call log_master_log(LOG_INFO, 'Minimum downdraft: '& //trim(conv_to_string(min_dwn))//' at level '//trim(conv_to_string(min_dwn_k)) ) end if ! show_critical_w - end subroutine calculate_w_percentiles + end subroutine calculate_w_percentiles !> Combines with MergeSort sorting algorithm taken from: @@ -240,14 +240,14 @@ end subroutine calculate_w_percentiles ! and modified to match local type and renamed to avoid confusion with intrinsic merge ! All parameters based on MergeSort. No need to modify anything. subroutine MergeSortMerge(A,NA,B,NB,C,NC) - + integer, intent(in) :: NA,NB,NC ! Normal usage: NA+NB = NC real(kind=DEFAULT_PRECISION), intent(in out) :: A(NA) ! B overlays C(NA+1:NC) real(kind=DEFAULT_PRECISION), intent(in) :: B(NB) real(kind=DEFAULT_PRECISION), intent(in out) :: C(NC) - + integer :: I,J,K - + I = 1; J = 1; K = 1; do while(I <= NA .and. J <= NB) if (A(I) <= B(J)) then @@ -265,9 +265,9 @@ subroutine MergeSortMerge(A,NA,B,NB,C,NC) K = K + 1 enddo return - + end subroutine mergesortmerge - + !> Combines with MergeSortMerge sorting algorithm taken from: ! https://rosettacode.org/wiki/Sorting_algorithms/Merge_sort#Fortran ! and modified to match local type @@ -275,14 +275,14 @@ end subroutine mergesortmerge !! @N size of A !! @T I don't really understand T recursive subroutine MergeSort(A,N,T) - + integer, intent(in) :: N real(kind=DEFAULT_PRECISION), dimension(N), intent(in out) :: A real(kind=DEFAULT_PRECISION), dimension((N+1)/2), intent (out) :: T - + integer :: NA,NB real(kind=DEFAULT_PRECISION) :: V - + if (N < 2) return if (N == 2) then if (A(1) > A(2)) then @@ -291,19 +291,19 @@ recursive subroutine MergeSort(A,N,T) A(2) = V endif return - endif + endif NA=(N+1)/2 NB=N-NA - + call MergeSort(A,NA,T) call MergeSort(A(NA+1),NB,T) - + if (A(NA) > A(NA+1)) then T(1:NA)=A(1:NA) call MergeSortMerge(T,NA,A(NA+1),NB,A,N) endif return - + end subroutine MergeSort @@ -341,7 +341,7 @@ subroutine field_value_retrieval_callback(current_state, name, field_value) else if (name .eq. "critical_downdraft_local") then allocate(field_value%real_1d_array(current_state%local_grid%size(Z_INDEX)), & source=current_state%global_grid%configuration%vertical%w_dwn(:)) - end if + end if end subroutine field_value_retrieval_callback end module pdf_analysis_mod diff --git a/components/profile_diagnostics/makefile b/components/profile_diagnostics/makefile index af28a2a4..6667c184 100644 --- a/components/profile_diagnostics/makefile +++ b/components/profile_diagnostics/makefile @@ -2,7 +2,9 @@ SRCSF = src/profile_diagnostics.F90 BUILDDIR=build COREDIR=../../model_core/build -FFLAGS=-I $(BUILDDIR) -I $(COREDIR) $(COMPILERFFLAGS) +TVDDIR=../tvdadvection/build + +FFLAGS=-I $(BUILDDIR) -I $(COREDIR) -I $(TVDDIR) $(COMPILERFFLAGS) OBJS = $(patsubst %.F90,$(BUILDDIR)/%.o,$(SRCSF)) all: create-build-dirs $(OBJS) diff --git a/components/setfluxlook/src/setfluxlook.F90 b/components/setfluxlook/src/setfluxlook.F90 index 0e9a118e..bad21d1c 100644 --- a/components/setfluxlook/src/setfluxlook.F90 +++ b/components/setfluxlook/src/setfluxlook.F90 @@ -27,12 +27,12 @@ module setfluxlook_mod SURFACE_HUMIDITIES_KEY = "surface_humidity", & !< NetCDF data surface_humidities SURFACE_SHF_KEY = "surface_sensible_heat_flux",& !< NetCDF data surface_sensible_heat_flux SURFACE_LHF_KEY = "surface_latent_heat_flux" !< NetCDF data surface_latent_heat_flux - + integer, parameter :: LOOKUP_ENTRIES = 80 !< Number of entries for MO lookup tables integer, parameter :: MAX_FILE_LEN=200 !< Maximum length of surface condition input filename integer, parameter :: MAX_SURFACE_INPUTS=750 !< Specifies the maximum number of surface inputs through configuration file !! Inputs through netcdf files are not limitted by this. - + character(MAX_FILE_LEN) :: input_file character(len=STRING_LENGTH) :: units_surface_temp='unset' ! units of theta variable forcing @@ -46,7 +46,7 @@ module setfluxlook_mod integer :: iqv ! index for vapour - public setfluxlook_get_descriptor + public setfluxlook_get_descriptor, set_look, change_look contains @@ -61,7 +61,7 @@ end function setfluxlook_get_descriptor subroutine initialisation_callback(current_state) type(model_state_type), intent(inout), target :: current_state - + current_state%lookup_table_entries=LOOKUP_ENTRIES call read_configuration(current_state) @@ -84,11 +84,11 @@ subroutine initialisation_callback(current_state) current_state%surface_vapour_flux = surface_latent_heat_flux(1) & /(current_state%global_grid%configuration%vertical%rho(1)*rlvap) end if - - current_state%fbuoy=0. + + current_state%fbuoy=0. if(.not. current_state%passive_th) current_state%fbuoy=& current_state%global_grid%configuration%vertical%buoy_co(1)*current_state%surface_temperature_flux - if(.not. current_state%passive_q .and. current_state%number_q_fields > 0)then + if(.not. current_state%passive_q .and. current_state%number_q_fields > 0)then current_state%fbuoy=current_state%fbuoy+current_state%cq(iqv)*current_state%surface_vapour_flux*G end if call set_look(current_state) ! _set M-O lookup table @@ -101,17 +101,17 @@ subroutine initialisation_callback(current_state) if (current_state%use_time_varying_surface_values) then call set_flux(current_state) else - ! If surface_values are constant then surface_temperatures prescribed in + ! If surface_values are constant then surface_temperatures prescribed in ! config and read in read_configuration but if humidity not set then ! surface vapour (surface_vapour_mixing_ratio) set to saturated value (see read_config) if (current_state%saturated_surface)then current_state%surface_vapour_mixing_ratio = qsaturation(surface_temperatures(1),current_state%surface_pressure*0.01) - else + else current_state%surface_vapour_mixing_ratio = & options_get_real(current_state%options_database, "surface_vapour_mixing_ratio") endif - - ! The code below copied from set_flux as these values need to be + + ! The code below copied from set_flux as these values need to be ! set for both time varying and constant surface values ! Set theta_v current_state%theta_surf = surface_temperatures(1)*& @@ -121,7 +121,7 @@ subroutine initialisation_callback(current_state) current_state%theta_virtual_surf = current_state%theta_surf + & current_state%global_grid%configuration%vertical%thref(2)* & current_state%cq(iqv)*current_state%surface_vapour_mixing_ratio - end if + end if ! Finally set up new values of THVSURF dependent constants current_state%cmbc=betam*current_state%global_grid%configuration%vertical%zn(2)*G*& @@ -160,10 +160,10 @@ subroutine set_look(current_state) smth=0.1_DEFAULT_PRECISION ! _relaxation parameter for unstable case do ik=1, current_state%lookup_table_entries current_state%lookup_table_velocity(ik)=current_state%velmin*(current_state%velmax/current_state%velmin)**& - (real(ik-1)/real(current_state%lookup_table_entries-1)) + (real(ik-1)/real(current_state%lookup_table_entries-1)) current_state%lookup_table_ustr(ik)=current_state%lookup_table_velocity(ik)*& current_state%global_grid%configuration%vertical%vk_on_zlogm - if (current_state%fbuoy .gt. 0.0_DEFAULT_PRECISION) then ! _unstable + if (current_state%fbuoy .gt. 0.0_DEFAULT_PRECISION) then ! _unstable iters(ik)=0 do i=1, 30 ! @check how many iterations needed!! iters(ik)=iters(ik)+1 @@ -193,6 +193,7 @@ subroutine change_look(current_state) dvelustr(current_state%lookup_table_entries), ob, x1, x0 integer n, ik, & ! Loop counters nit ! Number of iterations + if (current_state%fbuoynew .le. 0.0_DEFAULT_PRECISION) then current_state%fbuoy=current_state%fbuoynew return @@ -202,17 +203,20 @@ subroutine change_look(current_state) nit = 1+int(rnit) else nit=1 - end if + end if + if (nit .gt. 10 .or. (current_state%fbuoynew .gt. 0.0_DEFAULT_PRECISION .and. & current_state%fbuoy .le. 0.0_DEFAULT_PRECISION)) then current_state%fbuoy=current_state%fbuoynew call set_look(current_state) - return + return end if + crelax=1./sqrt(real(nit)) ! maybe better with crelax=1. dfb=(current_state%fbuoynew-current_state%fbuoy)/real(nit) + do n=1, nit - current_state%fbuoy=current_state%fbuoy+dfb + current_state%fbuoy=current_state%fbuoy+dfb do ik=2, current_state%lookup_table_entries-1 dvelustr(ik)=log(current_state%lookup_table_velocity(ik+1)/current_state%lookup_table_velocity(ik-1))& /log(current_state%lookup_table_ustr(ik+1)/current_state%lookup_table_ustr(ik-1)) @@ -223,40 +227,42 @@ subroutine change_look(current_state) current_state%lookup_table_velocity(current_state%lookup_table_entries-1))/& log(current_state%lookup_table_ustr(current_state%lookup_table_entries)/& current_state%lookup_table_ustr(current_state%lookup_table_entries-1)) - do ik=1, current_state%lookup_table_entries - ! compute new mean vel. based on old ustar and new fbuoy + do ik=1, current_state%lookup_table_entries + ! compute new mean vel. based on old ustar and new fbuoy ob=-current_state%lookup_table_ustr(ik)**3/(von_karman_constant*current_state%fbuoy) x1=sqrt(sqrt(1.-gammam*(current_state%global_grid%configuration%vertical%zn(2)+z0)/ob)) x0=sqrt(sqrt(1.-gammam*z0/ob)) velnew=(current_state%lookup_table_ustr(ik)/von_karman_constant)*(current_state%global_grid%configuration%vertical%zlogm-& (2.0_DEFAULT_PRECISION*log((x1+1.0_DEFAULT_PRECISION)/(x0+1.0_DEFAULT_PRECISION)) + & log((x1*x1+1.0_DEFAULT_PRECISION)/(x0*x0+1.0_DEFAULT_PRECISION)) + 2.0_DEFAULT_PRECISION*atan(x0) & - -2.0_DEFAULT_PRECISION*atan(x1))) - ! relax to new ustar + -2.0_DEFAULT_PRECISION*atan(x1))) + ! relax to new ustar current_state%lookup_table_ustr(ik)=current_state%lookup_table_ustr(ik)/& (velnew/current_state%lookup_table_velocity(ik))**(crelax/dvelustr(ik)) end do end do + current_state%cneut=current_state%lookup_table_ustr(current_state%lookup_table_entries)/& current_state%lookup_table_velocity(current_state%lookup_table_entries) current_state%cfc=current_state%lookup_table_ustr(1)*current_state%lookup_table_velocity(1)**convective_limit ! _Businger-Dyer current_state%fbuoy=current_state%fbuoynew + end subroutine change_look subroutine set_flux(current_state) type(model_state_type), intent(inout), target :: current_state real(kind=DEFAULT_PRECISION) :: surface_temp ! Surface temperature - + if (current_state%type_of_surface_boundary_conditions == PRESCRIBED_SURFACE_FLUXES) then ! Prescribed surface fluxes - + ! Linear interpolation of input data... - call interpolate_point_linear_1d(surface_boundary_input_times, & + call interpolate_point_linear_1d(surface_boundary_input_times, & surface_sensible_heat_flux/(current_state%global_grid%configuration%vertical%rho(1)*cp), & current_state%time, current_state%surface_temperature_flux, & - extrapolate='constant') - - call interpolate_point_linear_1d(surface_boundary_input_times, & + extrapolate='constant') + + call interpolate_point_linear_1d(surface_boundary_input_times, & surface_latent_heat_flux/(current_state%global_grid%configuration%vertical%rho(1)*rlvap), & current_state%time, current_state%surface_vapour_flux, & extrapolate='constant') @@ -265,18 +271,18 @@ subroutine set_flux(current_state) current_state%fbuoynew=0.0_DEFAULT_PRECISION if (.not. current_state%passive_th) current_state%fbuoynew=& current_state%global_grid%configuration%vertical%buoy_co(1)*current_state%surface_temperature_flux - if (.not. current_state%passive_q) then + if (.not. current_state%passive_q) then current_state%fbuoynew=current_state%fbuoynew+current_state%cq(iqv)*current_state%surface_vapour_flux*G end if else if (current_state%type_of_surface_boundary_conditions == PRESCRIBED_SURFACE_VALUES) then ! Prescribed surface temperatures ! Linear interpolation of input data... - call interpolate_point_linear_1d(surface_boundary_input_times, & + call interpolate_point_linear_1d(surface_boundary_input_times, & surface_temperatures, & current_state%time, surface_temp, & - extrapolate='constant') - + extrapolate='constant') + select case(trim(units_surface_temp)) case(degC) ! degrees C surface_temp = surface_temp + 273.15_DEFAULT_PRECISION @@ -286,10 +292,10 @@ subroutine set_flux(current_state) if (current_state%saturated_surface)then current_state%surface_vapour_mixing_ratio = qsaturation(surface_temp,current_state%surface_pressure*0.01) else - call interpolate_point_linear_1d(surface_boundary_input_times, & + call interpolate_point_linear_1d(surface_boundary_input_times, & surface_humidities, & current_state%time, current_state%surface_vapour_mixing_ratio, & - extrapolate='constant') + extrapolate='constant') end if ! Set theta_v @@ -310,7 +316,7 @@ subroutine set_flux(current_state) end if - end subroutine set_flux + end subroutine set_flux subroutine read_configuration(current_state) type(model_state_type), intent(inout), target :: current_state @@ -321,7 +327,7 @@ subroutine read_configuration(current_state) integer :: number_input_humidities - current_state%use_surface_boundary_conditions= & + current_state%use_surface_boundary_conditions= & options_get_logical(current_state%options_database, "use_surface_boundary_conditions") if (current_state%use_surface_boundary_conditions)then @@ -329,19 +335,19 @@ subroutine read_configuration(current_state) "type_of_surface_boundary_conditions") current_state%use_time_varying_surface_values=options_get_logical(current_state%options_database, & "use_time_varying_surface_values") - + current_state%saturated_surface = .true. ! We will change this if we find some humidity data input_file=options_get_string(current_state%options_database, "surface_conditions_file") ! Read in the input_file if (trim(input_file)=='' .or. trim(input_file)=='None')then - if (current_state%use_time_varying_surface_values)then + if (current_state%use_time_varying_surface_values)then allocate(surface_boundary_input_times(MAX_SURFACE_INPUTS)) surface_boundary_input_times=0.0 call options_get_real_array(current_state%options_database, "surface_boundary_input_times", surface_boundary_input_times) end if if (current_state%type_of_surface_boundary_conditions == PRESCRIBED_SURFACE_FLUXES)then - allocate(surface_sensible_heat_flux(MAX_SURFACE_INPUTS), & + allocate(surface_sensible_heat_flux(MAX_SURFACE_INPUTS), & surface_latent_heat_flux(MAX_SURFACE_INPUTS) & ) surface_sensible_heat_flux=0.0 @@ -349,7 +355,7 @@ subroutine read_configuration(current_state) call options_get_real_array(current_state%options_database, "surface_sensible_heat_flux", surface_sensible_heat_flux) call options_get_real_array(current_state%options_database, "surface_latent_heat_flux", surface_latent_heat_flux) number_input_humidities=0 - else if (current_state%type_of_surface_boundary_conditions == PRESCRIBED_SURFACE_VALUES) then + else if (current_state%type_of_surface_boundary_conditions == PRESCRIBED_SURFACE_VALUES) then allocate(surface_temperatures(MAX_SURFACE_INPUTS), & surface_humidities(MAX_SURFACE_INPUTS) & ) @@ -357,7 +363,7 @@ subroutine read_configuration(current_state) surface_humidities=0.0 call options_get_real_array(current_state%options_database, "surface_temperatures", surface_temperatures) units_surface_temp=options_get_string(current_state%options_database, "units_surface_temp") - + call options_get_real_array(current_state%options_database, "surface_humidities", surface_humidities) number_input_humidities=options_get_array_size(current_state%options_database, "surface_humidities") end if @@ -421,9 +427,9 @@ subroutine read_variables(filename, ncid, time_dim, time, surface_temperatures, integer :: status, variable_id - ! Do some checking on the variable contents so that we can deal with different + ! Do some checking on the variable contents so that we can deal with different ! variable names or missing variables - + ! time... status=nf90_inq_varid(ncid, TIME_KEY, variable_id) if (status==nf90_noerr)then @@ -432,25 +438,25 @@ subroutine read_variables(filename, ncid, time_dim, time, surface_temperatures, else call log_log(LOG_ERROR, "No recognized time variable found in"//trim(filename)) end if - + status=nf90_inq_varid(ncid, SURFACE_TEMPERATURES_KEY, variable_id) if (status==nf90_noerr)then allocate(surface_temperatures(time_dim)) call read_single_variable(ncid, SURFACE_TEMPERATURES_KEY, data1d=surface_temperatures) end if - + status=nf90_inq_varid(ncid, SURFACE_HUMIDITIES_KEY, variable_id) if (status==nf90_noerr)then allocate(surface_humidities(time_dim)) call read_single_variable(ncid, SURFACE_HUMIDITIES_KEY, data1d=surface_humidities) end if - + status=nf90_inq_varid(ncid, SURFACE_LHF_KEY, variable_id) if (status==nf90_noerr)then allocate(surface_latent_heat_flux(time_dim)) call read_single_variable(ncid, SURFACE_LHF_KEY, data1d=surface_latent_heat_flux) end if - + status=nf90_inq_varid(ncid, SURFACE_SHF_KEY, variable_id) if (status==nf90_noerr)then allocate(surface_sensible_heat_flux(time_dim)) diff --git a/components/socrates_couple/makefile b/components/socrates_couple/makefile index 47ce3ec6..a736e81a 100644 --- a/components/socrates_couple/makefile +++ b/components/socrates_couple/makefile @@ -1,35 +1,40 @@ -SRCSF = src/socrates_couple.F90 +SRCSF = src/socrates_couple_stub.F90 BUILDDIR=build -COREDIR=../../core/build +COREDIR=../../model_core/build FFLAGS=-I $(BUILDDIR) -I $(COREDIR) $(COMPILERFFLAGS) OBJS = $(patsubst %.F90,$(BUILDDIR)/%.o,$(SRCSF)) -ESRAD_RADIANCE_CORE=/data/local/fra23/MONC/es_radiation/src/radiance_core/ -ESRAD_MODULE_CORE=src/modules_core/ -ESRAD_CORE_DIR_LOCAL=src/um_edward_slingo_core_src -COPY_ESRAD_RADCORE=$(wildcard $(ESRAD_RADIANCE_CORE)/*.F90) -COPY_ESRAD_MODCORE=$(wildcard $(ESRAD_MODULE_CORE)/*.F90) -ESRAD_OBJS= $(patsubst %.F90, $(BUILDDIR)/%.o, $(notdir $(wildcard $(ESRAD_CORE_DIR_LOCAL)/*.F90))) +#ESRAD_RADIANCE_CORE=/data/local/fra23/MONC/es_radiation/src/radiance_core/ +#ESRAD_MODULE_CORE=src/modules_core/ +#ESRAD_CORE_DIR_LOCAL=src/um_edward_slingo_core_src +#COPY_ESRAD_RADCORE=$(wildcard $(ESRAD_RADIANCE_CORE)/*.F90) +#COPY_ESRAD_MODCORE=$(wildcard $(ESRAD_MODULE_CORE)/*.F90) +#ESRAD_OBJS= $(patsubst %.F90, $(BUILDDIR)/%.o, $(notdir $(wildcard $(ESRAD_CORE_DIR_LOCAL)/*.F90))) -all: create-esrad-dirs copy_radcore copy_modcore create-build-dirs $(ESRAD_OBJS) $(OBJS) +#all: create-esrad-dirs copy_radcore copy_modcore create-build-dirs $(ESRAD_OBJS) $(OBJS) -copy_radcore: - cp -f ${COPY_ESRAD_RADCORE} $(ESRAD_DIR_LOCAL) +#copy_radcore: +# cp -f ${COPY_ESRAD_RADCORE} $(ESRAD_DIR_LOCAL) -copy_modcore: - cp -f ${COPY_ESRAD_MODCORE} $(ESRAD_DIR_LOCAL) +#copy_modcore: +# cp -f ${COPY_ESRAD_MODCORE} $(ESRAD_DIR_LOCAL) -create-esrad-dirs: - mkdir -p $(ESRAD_DIR_LOCAL) +#create-esrad-dirs: +# mkdir -p $(ESRAD_DIR_LOCAL) -create-build-dirs: - mkdir -p $(BUILDDIR) +#create-build-dirs: +# mkdir -p $(BUILDDIR) -include $(ESRAD_DIR)/Mkdepend +#include $(ESRAD_DIR)/Mkdepend -$(ESRAD_OBJS) : $(BUILDDIR)/%.o : $(ESRAD_DIR_LOCAL)/%.F90 - $(FTN) $(OPT) $(FFLAGS) $< -o $(BUILDDIR)/$(notdir $@) +#$(ESRAD_OBJS) : $(BUILDDIR)/%.o : $(ESRAD_DIR_LOCAL)/%.F90 +# $(FTN) $(OPT) $(FFLAGS) $< -o $(BUILDDIR)/$(notdir $@) + +all: create-build-dirs $(OBJS) + +create-build-dirs: + mkdir -p $(BUILDDIR) $(OBJS) : $(BUILDDIR)/%.o : %.F90 $(FTN) $(OPT) $(FFLAGS) $< -o $(BUILDDIR)/$(notdir $@) diff --git a/components/subgrid_profile_diagnostics/makefile b/components/subgrid_profile_diagnostics/makefile index cdc9710c..d5a43e65 100644 --- a/components/subgrid_profile_diagnostics/makefile +++ b/components/subgrid_profile_diagnostics/makefile @@ -2,7 +2,8 @@ SRCSF = src/subgrid_profile_diagnostics.F90 BUILDDIR=build COREDIR=../../model_core/build -FFLAGS=-I $(BUILDDIR) -I $(COREDIR) $(COMPILERFFLAGS) +SMAGDIR=../smagorinsky/build +FFLAGS=-I $(BUILDDIR) -I $(COREDIR) -I $(SMAGDIR) $(COMPILERFFLAGS) OBJS = $(patsubst %.F90,$(BUILDDIR)/%.o,$(SRCSF)) all: create-build-dirs $(OBJS) diff --git a/components/tvdadvection/makefile b/components/tvdadvection/makefile index 60e74303..120287f0 100644 --- a/components/tvdadvection/makefile +++ b/components/tvdadvection/makefile @@ -1,4 +1,4 @@ -SRCSF = src/ultimateflux.F90 src/tvdadvection.F90 +SRCSF =src/def_tvd_diagnostic_terms.F90 src/ultimateflux.F90 src/tvdadvection.F90 BUILDDIR=build COREDIR=../../model_core/build diff --git a/dephy_20191209_12.mcf b/dephy_20191209_12.mcf new file mode 100644 index 00000000..b6bc3f06 --- /dev/null +++ b/dephy_20191209_12.mcf @@ -0,0 +1,155 @@ +# 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=.false. #Superseded by dephy +viscosity_enabled=.true. +diffusion_enabled=.true. +simplecloud_enabled=.true. +coriolis_enabled=.false. #Superseded by dephy +damping_enabled=.true. +forcing_enabled=.false. #Superseded by dephy +randomnoise_enabled=.true. +mean_profiles_enabled=.true. #This must be set to true if running with damping +th_advection_enabled=.true. +iobridge_enabled=.true. +scalar_diagnostics_enabled=.true. +profile_diagnostics_enabled=.true. +subgrid_profile_diagnostics_enabled=.true. +flux_budget_enabled=.true. +dephy_forcings_enabled=.true. +dephy_file="eurec4a_20191209_12_lag_local_prescribed.nc" +# Control configuration +display_synopsis_frequency=100 +termination_time=342000. +dtm=0.4 + +# IO server configuration +ioserver_configuration_file="io/io_cfg_files/data_write_1file.xml" +diagnostic_file="diagnostic_file_dephy.nc" +moncs_per_io_server=1 +sampling_frequency=60 +3d_sampling_frequency=600 +mm=600.0 +mm1=60.0 +diag_write_freq=7200.0 + +# Checkpoint configuration +checkpoint_frequency=0 +checkpoint_file="checkpoint_files_dephy_dump.nc" +check_walltime_frequency=200 +walltime_limit=30:00:00 + +# Advection choices +advection_flow_fields=pw +advection_theta_field=tvd +advection_q_fields=tvd + +# CFL configuration +cfl_frequency=10 +cfl_cvismax=0.4 +cfl_cvelmax=0.4 +cfl_dtmmax=2.0 +cfl_dtmmin=0.001 + +# Simple setup configuration +# We'll want to change this reference profile later +thref0=298.7 +surface_pressure=100000. +surface_reference_pressure=100000. +x_size=32 +y_size=32 +z_size=76 +dxx=200 +dyy=200 +zztop=5000.0 +kgd=1,76 +hgd=0.0,5000.0 +nsmth=80 +galilean_transformation=.false. + +enable_theta=.true. +number_q_fields=2 +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=.false. +z_init_pl_theta=0.0, 520.0, 1480., 2000., 3000. +f_init_pl_theta=298.7, 298.7, 302.4, 308.2, 311.85 +l_init_pl_u=.false. +z_init_pl_u=0.0, 700.0, 3000. +f_init_pl_u=-8.75, -8.75, -4.61 +l_init_pl_v=.false. +l_init_pl_q=.false. +names_init_pl_q=vapour +z_init_pl_q=0.0, 520.0, 1480., 2000., 3000. +f_init_pl_q=17.0e-3, 16.3e-3, 10.7e-3, 4.2e-3, 3.0e-3 + +l_matchthref=.false. + +# Random noise +l_rand_pl_theta=.true. +z_rand_pl_theta=0.0, 500.0, 501.0, 3000. +f_rand_pl_theta=0.5, 0.5, 0.0001, 0.0001 + +# Simple cloud +max_height_cloud=5000. + +# physical constants +z0=0.0002 +z0th=0.0002 + +# Coriolis +fcoriol=0.0000376 +baroclinicity_use_geostrophic_shear=.false. +geostrophic_wind_rate_of_change_in_x=0.0 +geostrophic_wind_rate_of_change_in_y=0.0 +surface_geostrophic_wind_x=-10. +surface_geostrophic_wind_y=0.0 + +# Damping configuration +dmptim=0.001 +zdmp=2300.0 +hdmp=2000.0 + +# Subsidence profile +l_subs_pl_theta=.false. +z_subs_pl=0.0, 1500.0, 2100.0, 3000. +f_subs_pl=0.0, -0.0065, 0.0, 0.0 +l_subs_pl_q=.false. + +#SUBSIDENCE=1, DIVERGENCE=0 +subsidence_input_type=1 +subsidence_local_theta=.true. +subsidence_local_q=.true. + +# Large-scale forcing +l_constant_forcing_theta=.true. +l_constant_forcing_q=.true. +l_constant_forcing_u=.false. +l_constant_forcing_v=.false. + +# Unit options are K/s or K/day +units_theta_force=K/day +l_constant_forcing_theta_z2pressure=.true. +z_force_pl_theta=0.0, 1500.0, 2500.0, 3000. +f_force_pl_theta=-2.0, -2.0, 0.0, 0.0 + +names_constant_forcing_q=vapour +z_force_pl_q=0.0, 300.0, 500.0, 3000. +f_force_pl_q=-1.2e-5, -1.2e-5, 0.0, 0.0 +# Unit options are kg/kg/s, kg/kg/day, g/kg/s or g/kg/day +units_q_force=g/kg/s + +# surface flux config +# type_of_surface_boundary_conditions=PRESCRIBED_FLUX=0 diff --git a/dephy_20191209_12_casim.mcf b/dephy_20191209_12_casim.mcf new file mode 100644 index 00000000..c9f29938 --- /dev/null +++ b/dephy_20191209_12_casim.mcf @@ -0,0 +1,166 @@ +# 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=.false. #Superseded by dephy +viscosity_enabled=.true. +diffusion_enabled=.true. +simplecloud_enabled=.false. +casim_enabled=.true. +coriolis_enabled=.false. #Superseded by dephy +damping_enabled=.true. +forcing_enabled=.false. #Superseded by dephy +randomnoise_enabled=.true. +mean_profiles_enabled=.true. #This must be set to true if running with damping +th_advection_enabled=.true. +iobridge_enabled=.true. +scalar_diagnostics_enabled=.true. +profile_diagnostics_enabled=.true. +subgrid_profile_diagnostics_enabled=.true. +flux_budget_enabled=.true. +dephy_forcings_enabled=.true. +dephy_file="eurec4a_20191209_12_lag_local_prescribed.nc" +# Control configuration +display_synopsis_frequency=100 +termination_time=342000. +dtm=0.4 + +# IO server configuration +ioserver_configuration_file="io/io_cfg_files/data_write_1file.xml" +diagnostic_file="diagnostic_file_dephy_casim.nc" +moncs_per_io_server=2 +sampling_frequency=60 +3d_sampling_frequency=600 +mm=600.0 +mm1=60.0 +diag_write_freq=7200.0 + +# Checkpoint configuration +checkpoint_frequency=0 +checkpoint_file="checkpoint_files_dephy_casim_dump.nc" +check_walltime_frequency=200 +walltime_limit=30:00:00 + +# Advection choices +advection_flow_fields=pw +advection_theta_field=tvd +advection_q_fields=tvd + +# CFL configuration +cfl_frequency=10 +cfl_cvismax=0.4 +cfl_cvelmax=0.4 +cfl_dtmmax=2.0 +cfl_dtmmin=0.001 + +# Simple setup configuration +# We'll want to change this reference profile later +thref0=298.7 +surface_pressure=100000. +surface_reference_pressure=100000. +x_size=32 +y_size=32 +z_size=76 +dxx=200 +dyy=200 +zztop=5000.0 +kgd=1,76 +hgd=0.0,5000.0 +nsmth=80 +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=.false. +z_init_pl_theta=0.0, 520.0, 1480., 2000., 3000. +f_init_pl_theta=298.7, 298.7, 302.4, 308.2, 311.85 +l_init_pl_u=.false. +z_init_pl_u=0.0, 700.0, 3000. +f_init_pl_u=-8.75, -8.75, -4.61 +l_init_pl_v=.false. +l_init_pl_q=.false. +names_init_pl_q=vapour +z_init_pl_q=0.0, 520.0, 1480., 2000., 3000. +f_init_pl_q=17.0e-3, 16.3e-3, 10.7e-3, 4.2e-3, 3.0e-3 + +l_matchthref=.false. + +# Random noise +l_rand_pl_theta=.true. +z_rand_pl_theta=0.0, 500.0, 501.0, 3000. +f_rand_pl_theta=0.5, 0.5, 0.0001, 0.0001 + +# Simple cloud +max_height_cloud=5000. + +# physical constants +z0=0.0002 +z0th=0.0002 + +# Coriolis +fcoriol=0.0000376 +baroclinicity_use_geostrophic_shear=.false. +geostrophic_wind_rate_of_change_in_x=0.0 +geostrophic_wind_rate_of_change_in_y=0.0 +surface_geostrophic_wind_x=-10. +surface_geostrophic_wind_y=0.0 + +# Damping configuration +dmptim=0.001 +zdmp=2300.0 +hdmp=2000.0 + +# Subsidence profile +l_subs_pl_theta=.false. +z_subs_pl=0.0, 1500.0, 2100.0, 3000. +f_subs_pl=0.0, -0.0065, 0.0, 0.0 +l_subs_pl_q=.false. + +#SUBSIDENCE=1, DIVERGENCE=0 +subsidence_input_type=1 +subsidence_local_theta=.true. +subsidence_local_q=.true. + +# Large-scale forcing +l_constant_forcing_theta=.true. +l_constant_forcing_q=.true. +l_constant_forcing_u=.false. +l_constant_forcing_v=.false. + +# Unit options are K/s or K/day +units_theta_force=K/day +l_constant_forcing_theta_z2pressure=.true. +z_force_pl_theta=0.0, 1500.0, 2500.0, 3000. +f_force_pl_theta=-2.0, -2.0, 0.0, 0.0 + +names_constant_forcing_q=vapour +z_force_pl_q=0.0, 300.0, 500.0, 3000. +f_force_pl_q=-1.2e-5, -1.2e-5, 0.0, 0.0 +# Unit options are kg/kg/s, kg/kg/day, g/kg/s or g/kg/day +units_q_force=g/kg/s + +# surface flux config +# type_of_surface_boundary_conditions=PRESCRIBED_FLUX=0 + +#CASIM options +number_q_fields=5 +option=22000 +l_warm=.true. + +aerosol_option=0 +iopt_act=0 +iopt_inuc=0 +process_level=0 +l_override_checks = .true. diff --git a/dephy_check.mcf b/dephy_check.mcf new file mode 100644 index 00000000..28122b7f --- /dev/null +++ b/dephy_check.mcf @@ -0,0 +1,155 @@ +# 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=.false. #Superseded by dephy +viscosity_enabled=.true. +diffusion_enabled=.true. +simplecloud_enabled=.true. +coriolis_enabled=.false. #Superseded by dephy +damping_enabled=.true. +forcing_enabled=.false. #Superseded by dephy +randomnoise_enabled=.true. +mean_profiles_enabled=.true. #This must be set to true if running with damping +th_advection_enabled=.true. +iobridge_enabled=.true. +scalar_diagnostics_enabled=.true. +profile_diagnostics_enabled=.true. +subgrid_profile_diagnostics_enabled=.true. +flux_budget_enabled=.true. +dephy_forcings_enabled=.true. +dephy_file="eurec4a_campaign_eulerian_eurec4a_dephy_prescribed_rad.nc" +# Control configuration +display_synopsis_frequency=100 +termination_time=7200. +dtm=0.4 + +# IO server configuration +ioserver_configuration_file="io/io_cfg_files/data_write_1file.xml" +diagnostic_file="diagnostic_file_dephy.nc" +moncs_per_io_server=1 +sampling_frequency=20 +3d_sampling_frequency=600 +mm=300.0 +mm1=30.0 +diag_write_freq=600.0 + +# Checkpoint configuration +checkpoint_frequency=0 +checkpoint_file="checkpoint_files_dephy_dump.nc" +check_walltime_frequency=100 +walltime_limit=01:00:00 + +# Advection choices +advection_flow_fields=pw +advection_theta_field=tvd +advection_q_fields=tvd + +# CFL configuration +cfl_frequency=10 +cfl_cvismax=0.4 +cfl_cvelmax=0.4 +cfl_dtmmax=0.4 +cfl_dtmmin=0.001 + +# Simple setup configuration +# We'll want to change this reference profile later +thref0=298.7 +surface_pressure=100000. +surface_reference_pressure=100000. +x_size=64 +y_size=64 +z_size=76 +dxx=100 +dyy=100 +zztop=5000.0 +kgd=1,76 +hgd=0.0,5000.0 +nsmth=80 +galilean_transformation=.false. + +enable_theta=.true. +number_q_fields=2 +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=.false. +z_init_pl_theta=0.0, 520.0, 1480., 2000., 3000. +f_init_pl_theta=298.7, 298.7, 302.4, 308.2, 311.85 +l_init_pl_u=.false. +z_init_pl_u=0.0, 700.0, 3000. +f_init_pl_u=-8.75, -8.75, -4.61 +l_init_pl_v=.false. +l_init_pl_q=.false. +names_init_pl_q=vapour +z_init_pl_q=0.0, 520.0, 1480., 2000., 3000. +f_init_pl_q=17.0e-3, 16.3e-3, 10.7e-3, 4.2e-3, 3.0e-3 + +l_matchthref=.false. + +# Random noise +l_rand_pl_theta=.true. +z_rand_pl_theta=0.0, 500.0, 501.0, 3000. +f_rand_pl_theta=0.5, 0.5, 0.0001, 0.0001 + +# Simple cloud +max_height_cloud=5000. + +# physical constants +z0=0.0002 +z0th=0.0002 + +# Coriolis +fcoriol=0.0000376 +baroclinicity_use_geostrophic_shear=.false. +geostrophic_wind_rate_of_change_in_x=0.0 +geostrophic_wind_rate_of_change_in_y=0.0 +surface_geostrophic_wind_x=-10. +surface_geostrophic_wind_y=0.0 + +# Damping configuration +dmptim=0.001 +zdmp=2300.0 +hdmp=2000.0 + +# Subsidence profile +l_subs_pl_theta=.false. +z_subs_pl=0.0, 1500.0, 2100.0, 3000. +f_subs_pl=0.0, -0.0065, 0.0, 0.0 +l_subs_pl_q=.false. + +#SUBSIDENCE=1, DIVERGENCE=0 +subsidence_input_type=1 +subsidence_local_theta=.true. +subsidence_local_q=.true. + +# Large-scale forcing +l_constant_forcing_theta=.true. +l_constant_forcing_q=.true. +l_constant_forcing_u=.false. +l_constant_forcing_v=.false. + +# Unit options are K/s or K/day +units_theta_force=K/day +l_constant_forcing_theta_z2pressure=.true. +z_force_pl_theta=0.0, 1500.0, 2500.0, 3000. +f_force_pl_theta=-2.0, -2.0, 0.0, 0.0 + +names_constant_forcing_q=vapour +z_force_pl_q=0.0, 300.0, 500.0, 3000. +f_force_pl_q=-1.2e-5, -1.2e-5, 0.0, 0.0 +# Unit options are kg/kg/s, kg/kg/day, g/kg/s or g/kg/day +units_q_force=g/kg/s + +# surface flux config +# type_of_surface_boundary_conditions=PRESCRIBED_FLUX=0 diff --git a/fcm-make/casim-dephy.cfg b/fcm-make/casim-dephy.cfg new file mode 100644 index 00000000..6913b984 --- /dev/null +++ b/fcm-make/casim-dephy.cfg @@ -0,0 +1,13 @@ +extract.ns = monc casim +#This will overide the default behaviour... + +extract.path-excl[monc] = / components/casim/src/casim_stub.F90 components/petsc_solver/src/petsc_solver.F90 components/casim_profile_dgs/src/casim_profile_dgs_stub.F90 components/socrates_couple/src/socrates_couple.F90 +extract.path-incl[monc] = components model_core io misc testcases monc_driver.F90 components/casim/src/casim_monc_diagnostics/casim_monc_dgs_space.F90 + +extract.location{primary}[casim] = fcm:casim.x_tr +$casim_revision{?} = 8166 +extract.location[casim] = @$casim_revision +extract.location{diff}[casim] = +extract.path-incl[casim] = src +extract.path-excl[casim] = / +preprocess.prop{fpp.defs}[casim] = DEF_MODEL=MODEL_MONC MODEL_MONC=4 diff --git a/global_config b/global_config index 57f840a1..64ede9a6 100644 --- a/global_config +++ b/global_config @@ -57,6 +57,7 @@ lateral_bcs_enabled=.false. conditional_diagnostics_column_enabled=.false. conditional_diagnostics_whole_enabled=.false. pdf_analysis_enabled=.false. +dephy_forcing_enabled=.false. # Default disable the test case components (individual user config will enable these) bubble_enabled=.false. @@ -73,7 +74,7 @@ solver_group_type=entire pressure-terms_group_type=column last_group_type=entire -start_group_contents=clearsourceterms, stepping_direction, halo_swapper, lateral_bcs, setfluxlook +start_group_contents=clearsourceterms, stepping_direction, halo_swapper, lateral_bcs, dephy_forcings, setfluxlook subgrid_group_contents=lower_bc, smagorinsky dynamics_group_contents=kidtestcase, pw_advection, tvd_advection, th_advection, diffusion, viscosity, coriolis, buoyancy, damping, forcing, set_consistent_lowbc, socrates_couple, lwrad_exponential, simplecloud, casim, flux_budget, subgrid_profile_diagnostics, diverr, psrce, diagnostics_3d, profile_diagnostics, casim_profile_dgs, scalar_diagnostics, stepfields @@ -83,7 +84,7 @@ last_group_contents=conditional_diagnostics_whole, checkpointer, model_synopsis, # Component ordering for other stages -initialisation_stage_ordering=decomposition, kidreader, kidtestcase, checkpointer, simplesetup, grid_manager, mean_profiles, swap_smooth, termination_check, simplecloud, casim, coriolis, buoyancy, cfltest, damping, diverr, fftsolver, halo_swapper, iterativesolver, setfluxlook, lower_bc, physicsa, psrce, pw_advection, diffusion, set_consistent_lowbc, viscosity, smagorinsky, stepfields, stepping_direction, tvd_advection, model_synopsis, socrates_couple, lwrad_exponential,th_advection, randomnoise, forcing, flux_budget, diagnostics_3d, profile_diagnostics, casim_profile_dgs, conditional_diagnostics_column, conditional_diagnostics_whole, pdf_analysis, subgrid_profile_diagnostics, scalar_diagnostics, lateral_bcs, petsc_solver, pstep, iobridge +initialisation_stage_ordering=decomposition, kidreader, kidtestcase, checkpointer, simplesetup, grid_manager, dephy_forcings, mean_profiles, swap_smooth, termination_check, simplecloud, casim, coriolis, buoyancy, cfltest, damping, diverr, fftsolver, halo_swapper, iterativesolver, setfluxlook, lower_bc, physicsa, psrce, pw_advection, diffusion, set_consistent_lowbc, viscosity, smagorinsky, stepfields, stepping_direction, tvd_advection, model_synopsis, socrates_couple, lwrad_exponential,th_advection, randomnoise, forcing, flux_budget, diagnostics_3d, profile_diagnostics, casim_profile_dgs, conditional_diagnostics_column, conditional_diagnostics_whole, pdf_analysis, subgrid_profile_diagnostics, scalar_diagnostics, lateral_bcs, petsc_solver, pstep, iobridge finalisation_stage_ordering=iobridge, checkpointer, diverr, fftsolver, grid_manager, halo_swapper, iterativesolver, physicsa, psrce, smagorinsky, tvd_advection,lwrad_exponential, model_synopsis, mean_profiles, pdf_analysis, forcing, stepfields, flux_budget, buoyancy, diffusion, lower_bc, viscosity, profile_diagnostics, pw_advection, th_advection, damping, simplecloud, pstep, conditional_diagnostics_whole, conditional_diagnostics_column # Control configuration diff --git a/makefile b/makefile index 733f4ec0..f491af56 100644 --- a/makefile +++ b/makefile @@ -46,7 +46,7 @@ endif COMPILERFFLAGS=-O3 COMPILERRECURSIVE= ACTIVE=-DU_ACTIVE -DV_ACTIVE -DW_ACTIVE -DUSE_MAKE -DEBUG_FLAGS=-g -fcheck=all -ffpe-trap=invalid,zero,overflow -fbacktrace -DDEBUG_MODE +DEBUG_FLAGS=-g -O0 -fcheck=all -ffpe-trap=invalid,zero,overflow -fbacktrace -DDEBUG_MODE FFLAGS=-I $(CORE_DIR)/$(BUILD_DIR) -I $(COMPONENTS_DIR)/$(BUILD_DIR) -I $(TESTCASE_DIR)/$(BUILD_DIR) -I $(IO_SERVER_DIR)/$(BUILD_DIR) $(COMPILERFFLAGS) LFLAGS=-L$(NETCDF_DIR)/lib -L./io -L misc/forthreads -L$(FFTW_DIR)/lib -L$(HDF5_DIR)/lib $(PETSC_LIB_PATH) -lnetcdff -lnetcdf -lhdf5 -lhdf5_hl -lz -lfftw3 -lpthread $(PETSC_LIBS) diff --git a/model_core/src/configuration/checkpointnetcdfparser.F90 b/model_core/src/configuration/checkpointnetcdfparser.F90 index 3a4d75fa..b7c6b0f6 100644 --- a/model_core/src/configuration/checkpointnetcdfparser.F90 +++ b/model_core/src/configuration/checkpointnetcdfparser.F90 @@ -19,7 +19,8 @@ module configuration_checkpoint_netcdf_parser_mod character(len=*), parameter :: OPTIONS_KEY="options_database", & !< The options key which references the configuration OPTIONS_DIM_KEY="number_options" !< Options dimension key - public parse_configuration_checkpoint_netcdf + public parse_configuration_checkpoint_netcdf,remove_null_terminator_from_string + contains !> Will parse the NetCDF checkpoint file and loads the configuration into the options database diff --git a/model_core/src/grid/interpolation.F90 b/model_core/src/grid/interpolation.F90 index e29924dc..b8c8c261 100644 --- a/model_core/src/grid/interpolation.F90 +++ b/model_core/src/grid/interpolation.F90 @@ -23,27 +23,27 @@ subroutine piecewise_linear_1d(zvals, vals, zgrid, field) real(kind=DEFAULT_PRECISION) :: field_lem(size(field)) real(kind=DEFAULT_PRECISION) :: verylow, veryhigh - + integer :: nn,k ! loop counters integer :: nnodes ! number of input values nnodes=size(zvals) - ! Code replicated from the LEM. This duplicates the interpolation of the + ! Code replicated from the LEM. This duplicates the interpolation of the ! LEM exactly - verylow = -1.0e5 - - initgd_lem(1) = vals(1) - zngd_lem(1) = verylow - DO nn=2,nnodes+1 + verylow = -1.0e5 + + initgd_lem(1) = vals(1) + zngd_lem(1) = verylow + DO nn=2,nnodes+1 initgd_lem(nn) = vals(nn-1) - zngd_lem(nn) = zvals(nn-1) + zngd_lem(nn) = zvals(nn-1) ENDDO do nn=1,nnodes DO k=1,size(field) - IF( zngd_lem(nn).LE.zgrid(k) .AND. zgrid(k).LT.zngd_lem(nn+1)) then + IF( zngd_lem(nn).LE.zgrid(k) .AND. zgrid(k).LT.zngd_lem(nn+1)) then field(k) = initgd_lem(nn) + ( initgd_lem(nn+1) - initgd_lem(nn))*(zgrid(k)- zngd_lem(nn))/ & (zngd_lem(nn+1) - zngd_lem(nn)) end IF @@ -54,7 +54,7 @@ subroutine piecewise_linear_1d(zvals, vals, zgrid, field) if (zgrid(size(field)) == zngd_lem(nnodes+1)) Then field(size(field)) = initgd_lem(nnodes+1) endif - + end subroutine piecewise_linear_1d !> Does a simple 1d linear interpolation to a point @@ -68,10 +68,10 @@ subroutine interpolate_point_linear_1d(zvals, vals, z, f, extrapolate) real(kind=DEFAULT_PRECISION), intent(in) :: z real(kind=DEFAULT_PRECISION), intent(out) :: f character(*), intent(in), optional :: extrapolate - + integer :: nn ! loop counter integer :: nnodes ! number of input values - + integer, parameter :: MAXCHARS=20 character(MAXCHARS) :: ext_type @@ -86,7 +86,7 @@ subroutine interpolate_point_linear_1d(zvals, vals, z, f, extrapolate) ext_type='linear' if (present(extrapolate))ext_type=trim(extrapolate) - if (z < zvals(1))then + if (z < zvals(1))then nn=1 select case (trim(ext_type)) case ('linear') @@ -97,10 +97,10 @@ subroutine interpolate_point_linear_1d(zvals, vals, z, f, extrapolate) end select return end if - + if (present(extrapolate))ext_type=trim(extrapolate) - if (z >= zvals(nnodes))then + if (z >= zvals(nnodes))then nn=nnodes select case (trim(ext_type)) case ('linear') @@ -129,7 +129,7 @@ end subroutine interpolate_point_linear_1d !! @param f output interpolated value subroutine piecewise_linear_2d(zvals, time_vals, vals, z, field) - ! Assumes input variables (vals) are 2-D, with dims (z, time) + ! Assumes input variables (vals) are 2-D, with dims (z, time) real(kind=DEFAULT_PRECISION), intent(in) :: zvals(:), time_vals(:) real(kind=DEFAULT_PRECISION), intent(in) :: vals(:,:) @@ -137,11 +137,11 @@ subroutine piecewise_linear_2d(zvals, time_vals, vals, z, field) real(kind=DEFAULT_PRECISION), intent(out) :: field(:,:) real(kind=DEFAULT_PRECISION) :: scale_tmp - + integer :: nn, k_monc, k_force ! loop counter integer :: nz_force, nt_force, nz_monc, nt_monc ! time and height array sizes for forcing and monc grids integer :: nnodes ! number of input values - + nz_force = size(zvals) nt_force = size(time_vals) nz_monc = size(z) @@ -150,45 +150,45 @@ subroutine piecewise_linear_2d(zvals, time_vals, vals, z, field) if ( zvals(1) .GT. zvals(nz_force) ) then ! pressure call log_master_log(LOG_ERROR, "Input forcing uses pressure, this has not been coded"// & " - please modify your forcing file to using height coordinates or modify the" // & - " interpolation routine in model_core to work with pressure coords - STOP") + " interpolation routine in model_core to work with pressure coords - STOP") else - do k_monc=2,nz_monc - do k_force=1,nz_force-1 - if( z(k_monc) >= zvals(k_force) .AND. z(k_monc) < zvals(k_force+1) ) then + do k_monc=2,nz_monc + do k_force=1,nz_force-1 + if( z(k_monc) >= zvals(k_force) .AND. z(k_monc) < zvals(k_force+1) ) then scale_tmp = ( z(k_monc) - zvals(k_force) ) / & - ( zvals(k_force+1) - zvals(k_force) ) - do nn=1, nt_force - field(k_monc,nn) = vals(k_force,nn) + & - ( vals(k_force+1,nn) - vals(k_force,nn) ) & - * scale_tmp + ( zvals(k_force+1) - zvals(k_force) ) + do nn=1, nt_force + field(k_monc,nn) = vals(k_force,nn) + & + ( vals(k_force+1,nn) - vals(k_force,nn) ) & + * scale_tmp enddo endif enddo enddo ! now examine the cases below and above forlevs(1) and forlevs(ktmfor - ! uses the local vertical gradient in the forcing to determine the - ! new values - do k_monc=2,nz_monc - if ( z(k_monc) >= zvals(nt_force) ) then - scale_tmp = ( z(k_monc) - zvals(nz_force) ) & - / ( zvals(nz_force) - zvals(nz_force-1) ) - do nn=1,nt_force - field(k_monc,nn) = vals(nz_force,nn) + & - ( vals(nz_force,nn) - vals(nz_force-1,nn) ) & - * scale_tmp + ! uses the local vertical gradient in the forcing to determine the + ! new values + do k_monc=2,nz_monc + if ( z(k_monc) >= zvals(nz_force) ) then + scale_tmp = ( z(k_monc) - zvals(nz_force) ) & + / ( zvals(nz_force) - zvals(nz_force-1) ) + do nn=1,nt_force + field(k_monc,nn) = vals(nz_force,nn) + & + ( vals(nz_force,nn) - vals(nz_force-1,nn) ) & + * scale_tmp enddo - elseif ( z(k_monc) < zvals(1) )THEN - scale_tmp = ( z(k_monc) - zvals(1) ) & - / ( zvals(1) - zvals(2) ) - do nn=1,nt_force - field(k_monc,nn) = vals(1,nn) + & - ( vals(1,nn) - vals(2,nn) ) & - * scale_tmp + elseif ( z(k_monc) < zvals(1) )THEN + scale_tmp = ( z(k_monc) - zvals(1) ) & + / ( zvals(1) - zvals(2) ) + do nn=1,nt_force + field(k_monc,nn) = vals(1,nn) + & + ( vals(1,nn) - vals(2,nn) ) & + * scale_tmp enddo endif enddo - ! - endif ! pressure or height + ! + endif ! pressure or height end subroutine piecewise_linear_2d @@ -200,10 +200,10 @@ subroutine interpolate_point_linear_2d(zvals, vals, z, f, extrapolate) real(kind=DEFAULT_PRECISION), intent(in) :: z real(kind=DEFAULT_PRECISION), intent(out) :: f(:) ! height character(*), intent(in), optional :: extrapolate - + integer :: nn ! loop counter integer :: nnodes ! number of input values - + integer, parameter :: MAXCHARS=20 character(MAXCHARS) :: ext_type @@ -217,7 +217,7 @@ subroutine interpolate_point_linear_2d(zvals, vals, z, f, extrapolate) ! ext_type='linear' if (present(extrapolate))ext_type=trim(extrapolate) - + ! test where model time is outside the bounds of the input forcing times ! 1) less than the lowest time level in forcing times @@ -226,11 +226,11 @@ subroutine interpolate_point_linear_2d(zvals, vals, z, f, extrapolate) f(:) = vals(:,nn) return end if - + !if (present(extrapolate))ext_type=trim(extrapolate) ! 2) greater than the last time level in the forcing times, make forcing constant - if (z >= zvals(nnodes))then + if (z >= zvals(nnodes))then nn=nnodes f(:) = vals(:,nn) return @@ -238,7 +238,7 @@ subroutine interpolate_point_linear_2d(zvals, vals, z, f, extrapolate) ! Time is within the bounds of the forcing input time - do nn = 1, nnodes-1 + do nn = 1, nnodes-1 if (zvals(nn) <= z .and. z < zvals(nn+1)) then select case (trim(ext_type)) case ('linear') @@ -251,7 +251,7 @@ subroutine interpolate_point_linear_2d(zvals, vals, z, f, extrapolate) endif enddo - + end subroutine interpolate_point_linear_2d end module interpolation_mod