Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adrianhill vn0.10.0 rc1 #14

Open
wants to merge 6 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 8 additions & 12 deletions components/buoyancy/src/buoyancy.F90
Original file line number Diff line number Diff line change
Expand Up @@ -29,8 +29,6 @@ module buoyancy_mod
real(kind=DEFAULT_PRECISION), dimension(:), allocatable :: tend_pr_tot_w
logical :: l_tend_pr_tot_w

integer :: diagnostic_generation_frequency

public buoyancy_get_descriptor

contains
Expand Down Expand Up @@ -170,9 +168,6 @@ subroutine initialisation_callback(current_state)
allocate( tend_pr_tot_w(current_state%local_grid%size(Z_INDEX)) )
endif

! Save the sampling_frequency to force diagnostic calculation on select time steps
diagnostic_generation_frequency=options_get_integer(current_state%options_database, "sampling_frequency")

end subroutine initialisation_callback


Expand All @@ -197,6 +192,10 @@ subroutine timestep_callback(current_state)

integer :: k, n
integer :: current_x_index, current_y_index, target_x_index, target_y_index
logical :: calculate_diagnostics

calculate_diagnostics = current_state%diagnostic_sample_timestep &
.and. .not. current_state%halo_column

current_x_index=current_state%column_local_x
current_y_index=current_state%column_local_y
Expand All @@ -210,10 +209,8 @@ subroutine timestep_callback(current_state)
endif
endif ! zero totals

if (mod(current_state%timestep, diagnostic_generation_frequency) == 0 .and. .not. current_state%halo_column) then
call save_precomponent_tendencies(current_state, current_x_index, current_y_index, target_x_index, target_y_index)
end if

if (calculate_diagnostics) &
call save_precomponent_tendencies(current_state, current_x_index, current_y_index, target_x_index, target_y_index)

#ifdef W_ACTIVE
if (.not. current_state%passive_th .and. current_state%th%active) then
Expand Down Expand Up @@ -252,9 +249,8 @@ subroutine timestep_callback(current_state)
end if
#endif

if (mod(current_state%timestep, diagnostic_generation_frequency) == 0 .and. .not. current_state%halo_column) then
call compute_component_tendencies(current_state, current_x_index, current_y_index, target_x_index, target_y_index)
end if
if (calculate_diagnostics) &
call compute_component_tendencies(current_state, current_x_index, current_y_index, target_x_index, target_y_index)

end subroutine timestep_callback

Expand Down
25 changes: 21 additions & 4 deletions components/casim/src/casim.F90
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,15 @@ module casim_mod
, l_pssub & ! sublimation of snow
, l_pgsub & ! sublimation of graupel
, l_pisub & ! sublimation of ice
, l_pimlt ! ice melting
, l_pimlt & ! ice melting
! New switches for sedimentation (these are sort-of temporary)
, l_gamma_online & ! when true use standard vn0.3.3 sed, when false use precalced gamma
, l_subseds_maxv & ! Use a CFL criteria based on max terminal velocity
! and sed_1M_2M
, l_sed_eulexp & ! switch for eulexp sed based on UM. Default is false
! so standard casim sed used
, cfl_vt_max & ! cfl limit for sedimentation (default = 1.0)
, l_kfsm

use micro_main, only: shipway_microphysics
use generic_diagnostic_variables, ONLY: casdiags, allocate_diagnostic_space, &
Expand All @@ -93,7 +101,7 @@ module casim_mod
, nc(:,:,:), qr(:,:,:), nr(:,:,:), m3r(:,:,:),rho(:,:,:) &
, exner(:,:,:), w(:,:,:), tke(:,:,:) &
, qi(:,:,:), ni(:,:,:), qs(:,:,:), ns(:,:,:), m3s(:,:,:) &
, qg(:,:,:), ng(:,:,:), m3g(:,:,:)
, qg(:,:,:), ng(:,:,:), m3g(:,:,:), cfliq(:,:,:), cfice(:,:,:)

REAL(wp), allocatable :: AccumSolMass(:,:,:), AccumSolNumber(:,:,:) ! Accumulation mode aerosol
REAL(wp), allocatable :: ActiveSolLiquid(:,:,:) ! Activated aerosol
Expand Down Expand Up @@ -277,6 +285,8 @@ subroutine initialisation_callback(current_state)
allocate(qg(kte,1,1))
allocate(ng(kte,1,1))
allocate(m3g(kte,1,1))
allocate(cfliq(kte,1,1))
allocate(cfice(kte,1,1))

allocate(AccumSolMass(kte,1,1))
allocate(AccumSolNumber(kte,1,1))
Expand Down Expand Up @@ -581,6 +591,7 @@ subroutine timestep_callback(current_state)
iqx = iql
qc(:,1,1) = current_state%zq(iqx)%data(:,jcol,icol)
dqc(:,1,1) = current_state%sq(iqx)%data(:,jcol,icol)
cfliq(:,1,1) = 1.0
end IF
IF (nq_r > 0)then
iqx = iqr
Expand Down Expand Up @@ -608,6 +619,7 @@ subroutine timestep_callback(current_state)
iqx = iqi
qi(:,1,1) = current_state%zq(iqx)%data(:,jcol,icol)
dqi(:,1,1) = current_state%sq(iqx)%data(:,jcol,icol)
cfice(:,1,1) = 1.0
end IF
IF (nq_s > 0)then
iqx = iqs
Expand Down Expand Up @@ -712,7 +724,7 @@ subroutine timestep_callback(current_state)
pressure, rho, &
w, tke, &
z_half, z_centre, &
dz, &
dz, cfliq, cfice, &
! in/out
dqv, dqc, dqr, dnc, dnr, dm3r, &
dqi, dqs, dqg, dni, dns, dng, dm3s, dm3g, &
Expand Down Expand Up @@ -838,7 +850,7 @@ subroutine timestep_callback(current_state)
! and surface
! snow rate (precip_s), which is the sum of ice, snow and graupel (See micromain.F90 in casim for
! calculation).
if (l_warm) then
if (l_warm .or. .not. casdiags % l_surface_snow ) then
surface_precip(target_y_index,target_x_index) = &
casdiags % SurfaceRainR(1,1)
else
Expand Down Expand Up @@ -947,6 +959,11 @@ subroutine read_configuration(current_state)
l_pgsub = options_get_logical(current_state%options_database, 'l_pgsub')
l_pisub = options_get_logical(current_state%options_database, 'l_pisub')
l_pimlt = options_get_logical(current_state%options_database, 'l_pimlt')
l_gamma_online = options_get_logical(current_state%options_database, 'l_gamma_online')
l_subseds_maxv = options_get_logical(current_state%options_database, 'l_subseds_maxv')
l_sed_eulexp = options_get_logical(current_state%options_database, 'l_sed_eulexp')
cfl_vt_max = options_get_real(current_state%options_database, 'cfl_vt_max')
l_kfsm = options_get_logical(current_state%options_database, 'l_kfsm')

end subroutine read_configuration

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -354,9 +354,12 @@ subroutine populate_casim_monc_dg(current_state, casdiags )
casdiags % dqr(1,1,:)

if (.not. l_warm) then
if ( casdiags % l_precip ) &
casim_monc_dgs % precip(target_y_index,target_x_index) = &
if ( casdiags % l_precip .and. casdiags % l_surface_snow ) &
casim_monc_dgs % precip(target_y_index,target_x_index) = &
casdiags % SurfaceRainR(1,1) + casdiags % SurfaceSnowR(1,1)
if ( casdiags % l_precip .and. .not. casdiags % l_surface_snow ) &
casim_monc_dgs % precip(target_y_index,target_x_index) = &
casdiags % SurfaceRainR(1,1)
if ( casdiags % l_surface_snow ) &
casim_monc_dgs % SurfaceSnowR(target_y_index,target_x_index) = &
casdiags % SurfaceSnowR(1,1)
Expand Down
194 changes: 149 additions & 45 deletions components/cfltest/src/cfltest.F90
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,8 @@ module cfltest_mod
use state_mod, only : model_state_type, parallel_state_type
use collections_mod, only : map_type
use logging_mod, only : LOG_WARN, LOG_DEBUG, LOG_ERROR, LOG_INFO, &
log_log, log_get_logging_level, log_newline
log_log, log_get_logging_level, log_master_newline, &
log_master_log, log_is_master
use conversions_mod, only : conv_to_string
use optionsdatabase_mod, only : options_get_integer, options_get_real, options_get_logical
use grids_mod, only : Z_INDEX, Y_INDEX, X_INDEX
Expand All @@ -21,7 +22,7 @@ module cfltest_mod

!! Configuration options - all are optional and have default values
real(kind=DEFAULT_PRECISION) :: tollerance, cvismax, cvelmax, dtmmax, dtmmin, rincmax
logical l_monitor_cfl
logical l_monitor_cfl, l_constant_dtm

public cfltest_get_descriptor
contains
Expand Down Expand Up @@ -49,6 +50,7 @@ subroutine initialisation_callback(current_state)
rincmax=options_get_real(current_state%options_database, "cfl_rincmax")

l_monitor_cfl = options_get_logical(current_state%options_database,"cfl_monitor")
l_constant_dtm = options_get_logical(current_state%options_database,"l_constant_dtm")

allocate(current_state%abswmax(current_state%local_grid%local_domain_end_index(Z_INDEX)))
end subroutine initialisation_callback
Expand All @@ -61,31 +63,48 @@ subroutine timestep_callback(current_state)

real(kind=DEFAULT_PRECISION) :: cfl_number

if (mod(current_state%timestep, current_state%cfl_frequency) == 1 .or. &
current_state%timestep-current_state%start_timestep .le. current_state%cfl_frequency) then
current_state%cvel=0.0_DEFAULT_PRECISION
current_state%cvel_x=0.0_DEFAULT_PRECISION
current_state%cvel_y=0.0_DEFAULT_PRECISION
current_state%cvel_z=0.0_DEFAULT_PRECISION
if (current_state%normal_step) then

call perform_cfl_and_galilean_transformation_calculation(current_state)
if ((mod(current_state%timestep, current_state%cfl_frequency) == 1 .or. &
current_state%timestep-current_state%start_timestep .le. current_state%cfl_frequency) &
.or. current_state%timestep .ge. (current_state%last_cfl_timestep + current_state%cfl_frequency)) then

current_state%cvel=(current_state%cvel_x*current_state%global_grid%configuration%horizontal%cx+current_state%cvel_y*&
current_state%global_grid%configuration%horizontal%cy+current_state%cvel_z)*current_state%dtm
current_state%cvis=current_state%cvis*(current_state%dtm * 4)
current_state%last_cfl_timestep = current_state%timestep
current_state%cvel=0.0_DEFAULT_PRECISION
current_state%cvel_x=0.0_DEFAULT_PRECISION
current_state%cvel_y=0.0_DEFAULT_PRECISION
current_state%cvel_z=0.0_DEFAULT_PRECISION

call perform_cfl_and_galilean_transformation_calculation(current_state)

current_state%cvel=(current_state%cvel_x*current_state%global_grid%configuration%horizontal%cx+current_state%cvel_y*&
current_state%global_grid%configuration%horizontal%cy+current_state%cvel_z)*current_state%dtm
current_state%cvis=current_state%cvis*(current_state%dtm * 4)

cfl_number=current_state%cvis/cvismax+current_state%cvel/cvelmax

current_state%absolute_new_dtm=current_state%dtm
current_state%update_dtm=.false.
if (cfl_number .gt. 0.0_DEFAULT_PRECISION) then
if (cfl_number .lt. (1.0_DEFAULT_PRECISION-tollerance) .or. cfl_number .gt. (1.0_DEFAULT_PRECISION+tollerance)) then
current_state%absolute_new_dtm=current_state%dtm/cfl_number
end if
end if
end if ! evaluate the cfl

call update_dtm_based_on_absolute(current_state, cfl_number)

cfl_number=current_state%cvis/cvismax+current_state%cvel/cvelmax
if (current_state%time_basis .or. current_state%force_output_on_interval) &
call evaluate_time_basis(current_state)

else ! do not evaluate if taking reduced NON-normal_step

current_state%absolute_new_dtm=current_state%dtm
current_state%update_dtm=.false.
if (cfl_number .gt. 0.0_DEFAULT_PRECISION) then
if (cfl_number .lt. (1.0_DEFAULT_PRECISION-tollerance) .or. cfl_number .gt. (1.0_DEFAULT_PRECISION+tollerance)) then
current_state%absolute_new_dtm=current_state%dtm/cfl_number
end if
end if
end if
call update_dtm_based_on_absolute(current_state, cfl_number)

end if ! end check normal_step

current_state%cvis=0.0_DEFAULT_PRECISION

end subroutine timestep_callback

!> Updates the (new) dtm value, which is actioned after time step completion, based upon the absolute value. This is incremented
Expand All @@ -103,32 +122,47 @@ subroutine update_dtm_based_on_absolute(current_state, cfl_number)

current_state%dtm_new=min(current_state%dtm*(1.0_DEFAULT_PRECISION+rincmax), current_state%absolute_new_dtm, dtmmax)

!! --- Diagnostic Writing -----------------
if (current_state%parallel%my_rank==0) then
if (log_get_logging_level() .eq. LOG_DEBUG) then
call log_log(LOG_DEBUG, "dtm changed from "//trim(conv_to_string(current_state%dtm, 5))//" to "//&
trim(conv_to_string(current_state%dtm_new, 5)))
if (l_monitor_cfl) then
call log_master_log(LOG_INFO, " --- CFL Monitoring Information --- ")
if (l_constant_dtm) &
call log_master_log(LOG_INFO, " *** l_constant_dtm=.true. - NOT CHANGING ***")
call log_master_log(LOG_INFO, "dtm changed from "//trim(conv_to_string(current_state%dtm, 5))&
//" to "//trim(conv_to_string(current_state%dtm_new, 5)))
if (cfl_number .gt. 0.0) then
call log_master_log(LOG_INFO, "cfl_number : "//trim(conv_to_string(cfl_number))//" (change divisor)")
call log_master_log(LOG_INFO, "cvis : "//trim(conv_to_string(current_state%cvis)) )
call log_master_log(LOG_INFO, "cvel : "//trim(conv_to_string(current_state%cvel)))
!//", "//trim(conv_to_string(current_state%cvel_x))//", "//trim(conv_to_string(current_state%cvel_y))//", "//trim(conv_to_string(current_state%cvel_z)) )
else
call log_master_log(LOG_INFO, "dtm change due to ratcheting only. Target dtm unchanged.")
end if
if (current_state%dtm_new .lt. dtmmin) then
call log_log(LOG_ERROR, "Timestep too small, dtmnew="//trim(conv_to_string(current_state%dtm_new, 5))//&
" dtmmin="//trim(conv_to_string(dtmmin, 5)))
end if
if (l_monitor_cfl) then
call log_log(LOG_INFO, " --- CFL Monitoring Information --- ")
call log_log(LOG_INFO, "dtm changed from "//trim(conv_to_string(current_state%dtm, 5))//" to "//&
trim(conv_to_string(current_state%dtm_new, 5)))
if (cfl_number .gt. 0.0) then
call log_log(LOG_INFO, "cfl_number : "//trim(conv_to_string(cfl_number))//" (change divisor)")
call log_log(LOG_INFO, "cvis : "//trim(conv_to_string(current_state%cvis)) )
call log_log(LOG_INFO, "cvel : "//trim(conv_to_string(current_state%cvel)) )
else
call log_log(LOG_INFO, "dtm change due to ratcheting only. Target dtm unchanged.")
end if
call log_log(LOG_INFO, "target dtm : "//trim(conv_to_string(current_state%absolute_new_dtm)) )
call log_newline()
call log_master_log(LOG_INFO, "target dtm : "//trim(conv_to_string(current_state%absolute_new_dtm)) )
call log_master_newline()
end if ! l_monitor_cfl


!! --- Diagnostic Writing -----------------
call log_master_log(LOG_DEBUG, "dtm changed from "//trim(conv_to_string(current_state%dtm, 5))//" to "//&
trim(conv_to_string(current_state%dtm_new, 5)))

if (current_state%dtm_new .lt. dtmmin .and. .not. l_constant_dtm) then
call log_master_log(LOG_ERROR, "Timestep too small, dtmnew="//&
trim(conv_to_string(current_state%dtm_new, 5))//&
", dtmmin="//trim(conv_to_string(dtmmin, 5)))
end if

if (current_state%dtm_new .lt. current_state%dtm .and. l_constant_dtm) then
call log_master_log(LOG_WARN, "The CFL check would like to reduce the model timestep "//&
"to a value below the specified constant dtm: dtmnew="//&
trim(conv_to_string(current_state%dtm_new, 5))//&
", constant dtm="//trim(conv_to_string(current_state%dtm, 5)))
end if
end if

end if ! l_monitor_cfl
end if ! Diagnostic Writing
if (l_constant_dtm) then
current_state%dtm = options_get_real(current_state%options_database, "dtm")
current_state%dtm_new = current_state%dtm
current_state%update_dtm=.false.
end if
end subroutine update_dtm_based_on_absolute

Expand Down Expand Up @@ -212,4 +246,74 @@ subroutine get_global_values(local_zumin, local_zumax, local_zvmin, local_zvmax,
call mpi_allreduce(local_zumin, global_zumin, 1, PRECISION_TYPE, MPI_MIN, parallel_state%monc_communicator, ierr)
call mpi_allreduce(local_zvmin, global_zvmin, 1, PRECISION_TYPE, MPI_MIN, parallel_state%monc_communicator, ierr)
end subroutine get_global_values


!> Handle timestep adjustment for time_basis=.true.
!! @param current_state The current model state
subroutine evaluate_time_basis(current_state)
type(model_state_type), intent(inout), target :: current_state

real(kind=DEFAULT_PRECISION) :: projected_time, dtm_trial
integer :: sample_nts, next_sample_time, ts_to_next_cfl, next_step, interval

next_sample_time = minval(current_state%sampling(:)%next_time)

! Handle timestep adjustment for time_basis=.true.
! Reduces timestep when approaching the next sample time.
if (current_state%time_basis) then
ts_to_next_cfl = current_state%cfl_frequency &
- mod(current_state%timestep, current_state%cfl_frequency)
projected_time = current_state%time + current_state%dtm &
+ (current_state%dtm_new * ts_to_next_cfl) + dtmmin
if ( next_sample_time .gt. 0 .and. projected_time .ge. next_sample_time ) then
sample_nts = max(1, ceiling( (next_sample_time &
- (current_state%time + current_state%dtm) ) &
/ current_state%dtm_new ) )
current_state%dtm_new = (next_sample_time - (current_state%time + current_state%dtm) ) &
/ sample_nts
if (l_constant_dtm) then
current_state%dtm = options_get_real(current_state%options_database, "dtm")
current_state%dtm_new = current_state%dtm
sample_nts = nint((next_sample_time - (current_state%time + current_state%dtm) ) &
/ current_state%dtm_new )
end if

! Record the next sampling step for intervals matching the next time
where(next_sample_time .eq. current_state%sampling(:)%next_time) &
current_state%sampling(:)%next_step = current_state%timestep + sample_nts

current_state%update_dtm = .true.
current_state%normal_step = .false.

end if ! next_sample_time to be exceeded

! Handle force_output_on_interval by finding if time step can be reduced to hit
! output interval exactly on a sampling interval.
else if (current_state%force_output_on_interval) then

next_step = maxval(current_state%sampling(:)%next_step, &
(next_sample_time .eq. current_state%sampling(:)%next_time))
interval = maxval(current_state%sampling(:)%interval, &
(next_step .eq. current_state%sampling(:)%next_step) .and. &
(next_sample_time .eq. current_state%sampling(:)%next_time))
sample_nts = next_step - current_state%timestep + interval

if (sample_nts .eq. 0) sample_nts = interval
projected_time = (current_state%time + current_state%dtm) + sample_nts*current_state%dtm_new
dtm_trial = (next_sample_time - (current_state%time + current_state%dtm) ) &
/ sample_nts

if (dtm_trial .gt. current_state%dtm) then
return
else ! dtm should be reduced to hit output time
current_state%dtm_new = dtm_trial
current_state%update_dtm = .true.
current_state%normal_step = .false.
end if

end if ! time_basis or force_output_on_interval

end subroutine evaluate_time_basis


end module cfltest_mod
Loading