casim.F90

Go to the documentation of this file.

 
module casim_mod
  use datadefn_mod, only : default_precision
  use monc_component_mod, only : component_descriptor_type, &
       component_array_field_type, component_double_data_type, &
       component_field_value_type, component_field_information_type
  use state_mod, only : forward_stepping, model_state_type
  use grids_mod, only : z_index, y_index, x_index
  use science_constants_mod
  use q_indices_mod, only: get_q_index, standard_q_names
  use optionsdatabase_mod, only : options_get_real, options_get_logical, options_get_integer
  use logging_mod, only : log_info, log_error, log_master_log
  use registry_mod, only : is_component_enabled
 
  ! casim modules...
  use variable_precision, ONLY: wp
  use initialize, only: mphys_init
  use mphys_switches, only: set_mphys_switches,  &
     l_warm,                                     &
     nq_l, nq_r, nq_i, nq_s, nq_g,               &
     l_2mc, l_2mr, l_2mi, l_2ms, l_2mg,          &
     l_3mr, l_3ms, l_3mg,                        &
     soluble_modes, active_cloud, active_rain,   &
     insoluble_modes, active_ice, active_number, &
     isol, iinsol, option, aerosol_option,        &
     iopt_act, option, aerosol_option        &
     ! All process rates are read in using read_configuration
     ! routine in this component
     , l_aaut, l_aacc, l_aevp, l_ased, l_warm            &
     , l_inuc, iopt_rcrit, iopt_inuc, l_iaut, l_iacw                &
     , l_rain, l_boussinesq, diag_mu_option   &
     , l_sed_3mdiff, l_cons, l_abelshipway, l_sed_icecloud_as_1m         &
     , l_active_inarg2000, process_level, l_separate_rain, l_idep    &
     , max_step_length, max_sed_length, l_sg, l_g, l_passive        &
     , l_passive3m, l_limit_psd, l_override_checks &
     , max_mu, fix_mu, l_raci_g, l_onlycollect, l_inhom_revp &
     , l_tidy_conserve_e , l_tidy_conserve_q & 
     , l_pcond & ! Condensation  
     , l_praut & ! Autoconversion cloud -> rain
     , l_pracw & ! Accretion  cloud -> rain
     , l_pracr & ! aggregation of rain drops
     , l_prevp & ! evaporation of rain
     , l_psedl & ! sedimentation of cloud
     , l_psedr & ! sedimentation of rain
     , l_ptidy & ! tidying term 1
     , l_ptidy2 & ! tidying term 2
     , l_pinuc & ! ice nucleation
     , l_pidep & ! ice deposition
     , l_piacw & ! ice accreting water
     , l_psaut & ! ice autoconversion ice -> snow
     , l_psdep & ! vapour deposition onto snow
     , l_psacw & ! snow accreting water
     , l_pgdep & ! vapour deposition onto graupel
     , l_pseds & ! snow sedimentation
     , l_psedi & ! ice sedimentation
     , l_psedg & ! graupel sedimentation
     , l_psaci & ! snow accreting ice
     , l_praci & ! rain accreting ice
     , l_psacr & ! snow accreting rain
     , l_pgacr & ! graupel accreting rain
     , l_pgacw & ! graupel accreting cloud water
     , l_pgaci & ! graupel accreting ice
     , l_pgacs & ! graupel accreting snow
     , l_piagg & ! aggregation of ice particles
     , l_psagg & ! aggregation of snow particles
     , l_pgagg & ! aggregation of graupel particles
     , l_psbrk & ! break up of snow flakes
     , l_pgshd & ! shedding of liquid from graupel
     , l_pihal & ! hallet mossop
     , l_psmlt & ! snow melting
     , l_pgmlt & ! graupel melting
     , l_phomr & ! homogeneous freezing of rain
     , l_phomc & ! homogeneous freezing of cloud droplets
     , l_pssub & ! sublimation of snow
     , l_pgsub & ! sublimation of graupel
     , l_pisub & ! sublimation of ice
     , l_pimlt   ! ice melting
 
  use micro_main, only: shipway_microphysics
  use generic_diagnostic_variables, ONLY: casdiags, allocate_diagnostic_space, &
       deallocate_diagnostic_space
  use casim_monc_dgs_space, only: casim_monc_dgs, allocate_casim_monc_dgs_space, &
       populate_casim_monc_dg
 
  implicit none
 
#ifndef TEST_MODE
  private
#endif
 
  REAL(wp), allocatable :: theta(:,:,:), pressure(:,:,:),  &
     z_half(:,:,:), z_centre(:,:,:), dz(:,:,:), qv(:,:,:),qc(:,:,:) &
     , nc(:,:,:), qr(:,:,:), nr(:,:,:), m3r(:,:,:),rho(:,:,:) &
     , exner(:,:,:), w(:,:,:), tke(:,:,:)                               &
     , qi(:,:,:), ni(:,:,:), qs(:,:,:), ns(:,:,:), m3s(:,:,:) &
     , qg(:,:,:), ng(:,:,:), m3g(:,:,:) 
 
  REAL(wp), allocatable :: accumsolmass(:,:,:), accumsolnumber(:,:,:) ! Accumulation mode aerosol
  REAL(wp), allocatable :: activesolliquid(:,:,:)                      ! Activated aerosol
  REAL(wp), allocatable :: aitkensolmass(:,:,:), aitkensolnumber(:,:,:) ! Aitken mode aerosol
  REAL(wp), allocatable :: coarsesolmass(:,:,:), coarsesolnumber(:,:,:) ! Course mode aerosol
  REAL(wp), allocatable :: activesolrain(:,:,:)                      ! Activeated aerosol in rain
  REAL(wp), allocatable :: coarsedustmass(:,:,:), coarsedustnumber(:,:,:) ! Coarse Dust
  REAL(wp), allocatable :: activeinsolice(:,:,:)                      ! Activeated dust
  REAL(wp), allocatable :: activesolice(:,:,:)                      ! Activeated aerosol in ice
  REAL(wp), allocatable :: activeinsolliquid(:,:,:)                      ! Activeated dust in cloud
  REAL(wp), allocatable :: accuminsolmass(:,:,:)                      ! Accum mode dust mass
  REAL(wp), allocatable :: accuminsolnumber(:,:,:)                      ! Accum mode dust number
  REAL(wp), allocatable :: activesolnumber(:,:,:)                      ! Activated soluble number (if we need a tracer)
  REAL(wp), allocatable :: activeinsolnumber(:,:,:)                      ! Activated insoluble number (if we need a tracer)
  !!AH - these are dummy arrays to accommodate the hygroscopicity from
  !!     UKCA-MODE. At present these are only required so that the shipway_microphysics
  !!     argument list matches between the UM and MONC interface. These may become non-dummy 
  !!     arguments in the future
  REAL(wp), allocatable :: aitkensolbk(:,:,:)
  REAL(wp), allocatable :: accumsolbk(:,:,:) 
  REAL(wp), allocatable :: coarsesolbk(:,:,:)
 
  ! Tendency from other physics/advection/forcing
  ! NB This is tendency (/s) not an increment over the timestep
  REAL(wp), allocatable :: dqv(:,:,:), dth(:,:,:), dqc(:,:,:), dnc(:,:,:)              &
     , dqr(:,:,:), dnr(:,:,:), dm3r(:,:,:)                                  &
     , dqi(:,:,:), dni(:,:,:), dqs(:,:,:), dns(:,:,:), dm3s(:,:,:)      & 
     , dqg(:,:,:), dng(:,:,:), dm3g(:,:,:) 
 
  REAL(wp), allocatable :: daccumsolmass(:,:,:), daccumsolnumber(:,:,:) ! Accumulation mode aerosol
  REAL(wp), allocatable :: dactivesolliquid(:,:,:)                       ! Activated aerosol
  REAL(wp), allocatable :: daitkensolmass(:,:,:), daitkensolnumber(:,:,:) ! Aitken mode aerosol
  REAL(wp), allocatable :: dcoarsesolmass(:,:,:), dcoarsesolnumber(:,:,:) ! Course mode aerosol
  REAL(wp), allocatable :: dactivesolrain(:,:,:)                       ! Activeated aerosol in rain
  REAL(wp), allocatable :: dcoarsedustmass(:,:,:), dcoarsedustnumber(:,:,:) ! Dust
  REAL(wp), allocatable :: dactiveinsolice(:,:,:)                       ! Activeated dust
  REAL(wp), allocatable :: dactivesolice(:,:,:)                       ! Activeated aerosol in ice
  REAL(wp), allocatable :: dactiveinsolliquid(:,:,:)                       ! Activeated dust in cloud
  REAL(wp), allocatable :: daccuminsolmass(:,:,:)                      ! Accum mode dust mass
  REAL(wp), allocatable :: daccuminsolnumber(:,:,:)                      ! Accum mode dust number
  REAL(wp), allocatable :: dactivesolnumber(:,:,:)                      ! Activated soluble number (if we need a tracer)
  REAL(wp), allocatable :: dactiveinsolnumber(:,:,:)                      ! Activated insoluble number (if we need a tracer)
 
  REAL(wp), allocatable :: surface_precip(:,:)
 
  INTEGER ::   ils,ile, jls,jle, kls,kle, &
     its,ite, jts,jte, kts,kte
 
  INTEGER :: iqv=0, iql=0, iqr=0,  iqi=0, iqs=0,  iqg=0
  INTEGER ::        inl=0, inr=0,  ini=0, ins=0,  ing=0
  INTEGER ::               i3mr=0,        i3ms=0, i3mg=0 
  INTEGER ::                &
     i_accumsolmass=0,      &
     i_accumsolnumber=0,    &
     i_activesolliquid=0,   &
     i_aitkensolmass=0,     &
     i_aitkensolnumber=0,   &
     i_coarsesolmass=0,     &
     i_coarsesolnumber=0,   &
     i_activesolrain=0,     &
     i_coarsedustmass=0,    &
     i_coarsedustnumber=0,  &
     i_activeinsolice=0,    &
     i_activesolice=0,      &
     i_activeinsolliquid=0, &
     i_accuminsolmass=0,    &
     i_accuminsolnumber=0,  &
     i_activesolnumber=0,   &
     i_activeinsolnumber=0
 
  real(kind=default_precision), dimension(:), allocatable ::     &
       phomc_tot, pinuc_tot, pidep_tot, psdep_tot, piacw_tot, psacw_tot, psacr_tot, pisub_tot,   &
       pssub_tot, pimlt_tot, psmlt_tot, psaut_tot, psaci_tot, praut_tot, pracw_tot, prevp_tot,   &
       pgacw_tot, pgacs_tot, pgmlt_tot, pgsub_tot, psedi_tot, pseds_tot, psedr_tot, psedg_tot,   &
       psedl_tot, pcond_tot
  
  public casim_get_descriptor
contains
 
  type(component_descriptor_type) function casim_get_descriptor()
    casim_get_descriptor%name="casim"
    casim_get_descriptor%version=0.1
    casim_get_descriptor%initialisation=>initialisation_callback
    casim_get_descriptor%timestep=>timestep_callback
 
    ! Set up fields to be published for diagnostics. These are made available to
    ! IO server in the field_value_retrieval_callback at the end of this module
    !
    casim_get_descriptor%field_value_retrieval=>field_value_retrieval_callback
    casim_get_descriptor%field_information_retrieval=>field_information_retrieval_callback
    
    allocate(casim_get_descriptor%published_fields(27))
    
    casim_get_descriptor%published_fields(1)="surface_precip"    
    casim_get_descriptor%published_fields(2)="homogeneous_freezing_rate"
    casim_get_descriptor%published_fields(3)="ice_nucleations_rate"
    casim_get_descriptor%published_fields(4)="ice_deposition_rate"
    casim_get_descriptor%published_fields(5)="snow_deposition_rate"
    casim_get_descriptor%published_fields(6)="ice_acc_cloud_rate"
    casim_get_descriptor%published_fields(7)="snow_acc_cloud_rate"
    casim_get_descriptor%published_fields(8)="snow_acc_rain_rate"
    casim_get_descriptor%published_fields(9)="ice_sublime_rate"
    casim_get_descriptor%published_fields(10)="snow_sublime_rate"
    casim_get_descriptor%published_fields(11)="ice_melt_rate"
    casim_get_descriptor%published_fields(12)="snow_melt_rate"
    casim_get_descriptor%published_fields(13)="snow_autoconversion_rate"
    casim_get_descriptor%published_fields(14)="snow_acc_ice_rate"
    casim_get_descriptor%published_fields(15)="rain_autoconversion_rate"
    casim_get_descriptor%published_fields(16)="rain_acc_cloud_rate"
    casim_get_descriptor%published_fields(17)="rain_evap_rate"
    casim_get_descriptor%published_fields(18)="graup_acc_cloud_rate"
    casim_get_descriptor%published_fields(19)="graup_acc_snow_rate"
    casim_get_descriptor%published_fields(20)="graup_melt_rate"
    casim_get_descriptor%published_fields(21)="graup_sublime_rate"
    casim_get_descriptor%published_fields(22)="ice_sed_rate"
    casim_get_descriptor%published_fields(23)="snow_sed_rate"
    casim_get_descriptor%published_fields(24)="rain_sed_rate"
    casim_get_descriptor%published_fields(25)="graup_sed_rate"
    casim_get_descriptor%published_fields(26)="cloud_sed_rate"
    casim_get_descriptor%published_fields(27)="condensation_rate"
    
  end function casim_get_descriptor
 
  subroutine initialisation_callback(current_state)
    type(model_state_type), target, intent(inout) :: current_state
 
    integer :: y_size_local, x_size_local
 
    if (is_component_enabled(current_state%options_database, "simplecloud")) then
      call log_master_log(log_error, "Casim and Simplecloud are enabled, this does not work yet. Please disable one")
   end if
 
   !allocate(psedl_tot(current_state%local_grid%size(Z_INDEX)), &
   !      pcond_tot(current_state%local_grid%size(Z_INDEX)))
    
    y_size_local = current_state%local_grid%size(y_index)
    x_size_local = current_state%local_grid%size(x_index)
    
    call read_configuration(current_state)   
 
    ils=1
    ile=1
    jls=1
    jle=1
    kls=2
    kle=current_state%local_grid%size(z_index)
    its=1
    ite=1
    jts=1
    jte=1
    kts=1
    kte=current_state%local_grid%size(z_index)
 
 
    allocate(pressure(kte,1,1))
    allocate(z_half(0:kte,1,1))
    allocate(z_centre(kte,1,1))
    allocate(dz(kte,1,1))
    allocate(rho(kte,1,1))
    allocate(exner(kte,1,1))
    allocate(w(kte,1,1))
    allocate(tke(kte,1,1))
 
    allocate(theta(kte,1,1))
    allocate(qv(kte,1,1))
    allocate(qc(kte,1,1))
    allocate(nc(kte,1,1))
    allocate(qr(kte,1,1))
    allocate(nr(kte,1,1))
    allocate(m3r(kte,1,1))
    allocate(qi(kte,1,1))
    allocate(ni(kte,1,1))
    allocate(qs(kte,1,1))
    allocate(ns(kte,1,1))
    allocate(m3s(kte,1,1))
    allocate(qg(kte,1,1))
    allocate(ng(kte,1,1))
    allocate(m3g(kte,1,1))
 
    allocate(accumsolmass(kte,1,1))
    allocate(accumsolnumber(kte,1,1))
    allocate(activesolliquid(kte,1,1))
    allocate(aitkensolmass(kte,1,1))
    allocate(aitkensolnumber(kte,1,1))
    allocate(coarsesolmass(kte,1,1))
    allocate(coarsesolnumber(kte,1,1))
    allocate(activesolrain(kte,1,1))
    allocate(coarsedustmass(kte,1,1))
    allocate(coarsedustnumber(kte,1,1))
    allocate(activeinsolice(kte,1,1))
    allocate(activesolice(kte,1,1))
    allocate(activeinsolliquid(kte,1,1))
    allocate(accuminsolmass(kte,1,1))
    allocate(accuminsolnumber(kte,1,1))
    allocate(activesolnumber(kte,1,1))
    allocate(activeinsolnumber(kte,1,1))
    ! allocate the hygoscopicity arrays
    allocate(aitkensolbk(kte,1,1))
    allocate(accumsolbk(kte,1,1))
    allocate(coarsesolbk(kte,1,1))
 
    allocate(dth(kte,1,1))
    allocate(dqv(kte,1,1))
    allocate(dqc(kte,1,1))
    allocate(dnc(kte,1,1))
    allocate(dqr(kte,1,1))
    allocate(dnr(kte,1,1))
    allocate(dm3r(kte,1,1))
    allocate(dqi(kte,1,1))
    allocate(dni(kte,1,1))
    allocate(dqs(kte,1,1))
    allocate(dns(kte,1,1))
    allocate(dm3s(kte,1,1))
    allocate(dqg(kte,1,1))
    allocate(dng(kte,1,1))
    allocate(dm3g(kte,1,1))
 
    allocate(daccumsolmass(kte,1,1))
    allocate(daccumsolnumber(kte,1,1))
    allocate(dactivesolliquid(kte,1,1))
    allocate(daitkensolmass(kte,1,1))
    allocate(daitkensolnumber(kte,1,1))
    allocate(dcoarsesolmass(kte,1,1))
    allocate(dcoarsesolnumber(kte,1,1))
    allocate(dactivesolrain(kte,1,1))
    allocate(dcoarsedustmass(kte,1,1))
    allocate(dcoarsedustnumber(kte,1,1))
    allocate(dactiveinsolice(kte,1,1))
    allocate(dactivesolice(kte,1,1))
    allocate(dactiveinsolliquid(kte,1,1))
    allocate(daccuminsolmass(kte,1,1))
    allocate(daccuminsolnumber(kte,1,1))
    allocate(dactivesolnumber(kte,1,1))
    allocate(dactiveinsolnumber(kte,1,1))
 
    call set_mphys_switches(option,aerosol_option)
    call mphys_init(its, ite, jts, jte, kts, kte, ils, ile, jls, jle, kls, kle, l_tendency=.true.)
 
    ! Need to allocate the appropriate indices, e.g. iqv, iql...
    ! This needs to be compatible with the rest of the model
    ! This essentially reproduces the switching in the main microphysics
    ! code already done above (set_mphys_switches), so could be combined 
    ! once the MONC method has been finalized.
    ! Note the numbers assigned here may be different from those assigned 
    ! in the microphysics since we share the q array with other components. 
 
    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
 
    ! Mass
    iqv = get_q_index(standard_q_names%VAPOUR, 'casim')
    if (nq_l>0)then
       iql = get_q_index(standard_q_names%CLOUD_LIQUID_MASS, 'casim')
       current_state%cq(iql) = -1.0
    end if
    if (nq_r>0)then
       iqr = get_q_index(standard_q_names%RAIN_MASS, 'casim')
       current_state%rain_water_mixing_ratio_index=iqr
      current_state%cq(iqr) = -1.0
    end if
    if (.not. l_warm)then
      if (nq_i>0)then
         iqi = get_q_index(standard_q_names%ICE_MASS, 'casim')
         current_state%ice_water_mixing_ratio_index=iqi
        current_state%cq(iqi) = -1.0
      end if
      if (nq_s>0)then
         iqs = get_q_index(standard_q_names%SNOW_MASS, 'casim')
         current_state%snow_water_mixing_ratio_index=iqs
        current_state%cq(iqs) = -1.0
      end if
      if (nq_g>0)then
         iqg = get_q_index(standard_q_names%GRAUPEL_MASS, 'casim')
         current_state%graupel_water_mixing_ratio_index=iqg
        current_state%cq(iqg) = -1.0
      end if
    end if
 
    ! Number
    if (l_2mc)inl = get_q_index(standard_q_names%CLOUD_LIQUID_NUMBER, 'casim')
    if (l_2mr)inr = get_q_index(standard_q_names%RAIN_NUMBER, 'casim')
    if (.not. l_warm)then
      if (l_2mi)ini = get_q_index(standard_q_names%ICE_NUMBER, 'casim')
      if (l_2ms)ins = get_q_index(standard_q_names%SNOW_NUMBER, 'casim')
      if (l_2mg)ing = get_q_index(standard_q_names%GRAUPEL_NUMBER, 'casim')
    end if
 
    ! Third moments
    if (l_3mr)i3mr = get_q_index(standard_q_names%RAIN_THIRD_MOMENT, 'casim')
    if (.not. l_warm)then
      if (l_3ms)i3ms = get_q_index(standard_q_names%SNOW_THIRD_MOMENT, 'casim')
      if (l_3mg)i3mg = get_q_index(standard_q_names%GRAUPEL_THIRD_MOMENT, 'casim')
    end if
 
    ! Aerosol
    if (soluble_modes(1) > 1)   i_aitkensolmass     = &
       get_q_index(standard_q_names%AITKEN_SOL_MASS,     'casim') 
    if (soluble_modes(1) > 0)   i_aitkensolnumber   = &
       get_q_index(standard_q_names%AITKEN_SOL_NUMBER,   'casim') 
    if (soluble_modes(2) > 1)   i_accumsolmass      = &
       get_q_index(standard_q_names%ACCUM_SOL_MASS,      'casim') 
    if (soluble_modes(2) > 0)   i_accumsolnumber    = &
       get_q_index(standard_q_names%ACCUM_SOL_NUMBER,    'casim') 
    if (soluble_modes(3) > 1)   i_coarsesolmass     = &
       get_q_index(standard_q_names%COARSE_SOL_MASS,     'casim') 
    if (soluble_modes(3) > 0)   i_coarsesolnumber   = &
       get_q_index(standard_q_names%COARSE_SOL_NUMBER,   'casim') 
    if (active_cloud(isol))     i_activesolliquid   = &
       get_q_index(standard_q_names%ACTIVE_SOL_LIQUID,   'casim') 
    if (active_rain(isol))      i_activesolrain     = &
       get_q_index(standard_q_names%ACTIVE_SOL_RAIN,     'casim') 
    if (insoluble_modes(2) > 1) i_coarsedustmass    = &
       get_q_index(standard_q_names%COARSE_DUST_MASS,    'casim') 
    if (insoluble_modes(2) > 0) i_coarsedustnumber  = &
       get_q_index(standard_q_names%COARSE_DUST_NUMBER,  'casim') 
    if (active_ice(iinsol))     i_activeinsolice    = &
       get_q_index(standard_q_names%ACTIVE_INSOL_ICE,    'casim') 
    if (active_ice(isol))       i_activesolice      = &
       get_q_index(standard_q_names%ACTIVE_SOL_ICE,      'casim') 
    if (active_cloud(iinsol))   i_activeinsolliquid = &
       get_q_index(standard_q_names%ACTIVE_INSOL_LIQUID, 'casim') 
    if (insoluble_modes(1) > 1) i_accuminsolmass    = &
       get_q_index(standard_q_names%ACCUM_INSOL_MASS,    'casim') 
    if (insoluble_modes(1) > 0) i_accuminsolnumber  = &
       get_q_index(standard_q_names%ACCUM_INSOL_NUMBER,  'casim') 
    if (active_number(isol))    i_activesolnumber   = &
       get_q_index(standard_q_names%ACTIVE_SOL_NUMBER,   'casim') 
    if (active_number(iinsol))  i_activeinsolnumber = &
         get_q_index(standard_q_names%ACTIVE_INSOL_NUMBER, 'casim')
 
    ! set logicals for the microphysics diagnostics: process rates
    casdiags % l_dth = .true.
    casdiags % l_dqv = .true.
    if ( nq_l>0 ) casdiags % l_dqc = .true.
    if ( nq_r>0 ) casdiags % l_dqr = .true.
    if ( l_pcond ) casdiags % l_pcond = .true.
    if ( l_psedl ) then
       casdiags % l_psedl = .true.
       casdiags % l_surface_rain = .true.
       casdiags % l_precip = .true.
    endif
    if ( l_praut ) casdiags % l_praut = .true.
    if ( l_pracw ) casdiags % l_pracw = .true.
    if ( l_prevp ) casdiags % l_prevp = .true.
    if ( l_psedr ) then
       casdiags % l_psedr = .true.
       casdiags % l_surface_rain = .true.
       casdiags % l_precip = .true.
    endif
    if (.not. l_warm) then
       if ( nq_i>0 ) casdiags % l_dqi = .true.
       if ( nq_s>0 ) casdiags % l_dqs = .true.
       if ( nq_g>0 ) casdiags % l_dqg = .true.
       if ( l_phomc ) casdiags % l_phomc = .true.
       if ( l_pinuc ) casdiags % l_pinuc = .true.
       if ( l_pidep ) casdiags % l_pidep = .true.
       if ( l_piacw ) casdiags % l_piacw = .true.
       if ( l_pisub ) casdiags % l_pisub = .true.
       if ( l_pimlt ) casdiags % l_pimlt = .true.
       if ( l_psedi ) then
          casdiags % l_psedi = .true.
          casdiags % l_surface_snow = .true.
          casdiags % l_precip = .true.
       endif
       if ( l_psmlt ) casdiags % l_psmlt = .true.
       if ( l_psaut ) casdiags % l_psaut = .true.
       if ( l_psaci ) casdiags % l_psaci = .true.
       if ( l_psacw ) casdiags % l_psacw = .true.
       if ( l_psacr ) casdiags % l_psacr = .true.
       if ( l_pssub ) casdiags % l_pssub = .true.
       if ( l_psdep ) casdiags % l_psdep = .true.
       if ( l_pseds ) then
          casdiags % l_pseds = .true.
          casdiags % l_surface_snow = .true.
          casdiags % l_precip = .true.
       endif
       if ( l_pgacw ) casdiags % l_pgacw = .true.
       if ( l_pgacs ) casdiags % l_pgacs = .true.
       if ( l_pgmlt ) casdiags % l_pgmlt = .true.
       if ( l_pgsub ) casdiags % l_pgsub = .true.
       if ( l_psedg ) then
          casdiags % l_psedg = .true.
          casdiags % l_surface_graup = .true.
          casdiags % l_precip = .true.
       endif
    endif
 
    ! allocate diagnostic space in casdiags depending on the logicals defined above
    CALL allocate_diagnostic_space(its, ite, jts, jte, kts, kte)
    ! this is no longer needed since can use cas_monc_dgs structure but keep for now
    allocate(surface_precip(y_size_local, x_size_local))
    ! allocate diagnostic space for MONC fields to export to IO server
    call allocate_casim_monc_dgs_space(current_state, casdiags)
 
  end subroutine initialisation_callback
 
  subroutine timestep_callback(current_state)
    type(model_state_type), target, intent(inout) :: current_state
 
    REAL(wp) :: dtwp
    INTEGER :: icol, jcol, iqx, target_x_index, target_y_index
 
    icol=current_state%column_local_x
    jcol=current_state%column_local_y
    target_y_index=jcol-current_state%local_grid%halo_size(y_index)
    target_x_index=icol-current_state%local_grid%halo_size(x_index)
 
    !if (current_state%first_timestep_column) then
    !   psedl_tot(:)= 0.0_DEFAULT_PRECISION
    !   pcond_tot(:)= 0.0_DEFAULT_PRECISION
    !endif
    
    if (current_state%halo_column .or. current_state%timestep < 2) return
 
    if (current_state%field_stepping == forward_stepping)then
      call log_master_log(log_error, 'Currently, CASIM assumes CENTERED_STEPPING')
      dtwp = current_state%dtm
    else
      dtwp = 2.0*current_state%dtm
    end if
 
    ! Initialize aerosol fields to zero...
    aitkensolmass = 0.0
    daitkensolmass = 0.0
    aitkensolnumber = 0.0
    daitkensolnumber = 0.0
    accumsolmass = 0.0
    daccumsolmass = 0.0
    accumsolnumber = 0.0
    daccumsolnumber = 0.0
    coarsesolmass = 0.0
    dcoarsesolmass = 0.0
    coarsesolnumber = 0.0
    dcoarsesolnumber = 0.0
    activesolliquid = 0.0
    dactivesolliquid = 0.0
    coarsedustmass = 0.0
    dcoarsedustmass = 0.0
    coarsedustnumber = 0.0
    dcoarsedustnumber = 0.0
    activeinsolice = 0.0
    dactiveinsolice = 0.0
    activesolice = 0.0
    dactivesolice = 0.0
    activeinsolliquid = 0.0
    dactiveinsolliquid = 0.0
    accuminsolmass = 0.0
    daccuminsolmass = 0.0
    accuminsolnumber = 0.0
    daccuminsolnumber = 0.0
    activesolnumber = 0.0
    dactivesolnumber = 0.0
    activeinsolnumber = 0.0
    dactiveinsolnumber = 0.0
    aitkensolbk = 0.0
    accumsolbk = 0.0 
    coarsesolbk = 0.0
 
    theta(:,1,1) = current_state%zth%data(:, jcol, icol) + current_state%global_grid%configuration%vertical%thref(:)
    dth(:,1,1) = current_state%sth%data(:, jcol, icol)
    exner(:,1,1) = current_state%global_grid%configuration%vertical%rprefrcp(:)
    pressure(:,1,1) = current_state%global_grid%configuration%vertical%prefn(:)
    z_centre(:,1,1) = current_state%global_grid%configuration%vertical%zn(:)
    dz(:,1,1) = current_state%global_grid%configuration%vertical%dz(:)
    z_half(:kte-1,1,1) = current_state%global_grid%configuration%vertical%z(:)
    rho(:,1,1) = current_state%global_grid%configuration%vertical%rhon(:)
    w(:,1,1) = current_state%zw%data(:, jcol, icol)
    tke(:,1,1) = 0.1 ! Test value
 
    iqx = iqv
    qv(:,1,1) = current_state%zq(iqx)%data(:,jcol,icol)
    dqv(:,1,1) = current_state%sq(iqx)%data(:,jcol,icol)
 
    ! Warm microphysical fields
    IF (nq_l > 0)then
      iqx = iql
      qc(:,1,1) = current_state%zq(iqx)%data(:,jcol,icol)
      dqc(:,1,1) = current_state%sq(iqx)%data(:,jcol,icol)
    end IF
    IF (nq_r > 0)then
      iqx = iqr
      qr(:,1,1) = current_state%zq(iqx)%data(:,jcol,icol)
      dqr(:,1,1) = current_state%sq(iqx)%data(:,jcol,icol)
    end IF
    IF (nq_l > 1)then
      iqx = inl
      nc(:,1,1) = current_state%zq(iqx)%data(:,jcol,icol)
      dnc(:,1,1) = current_state%sq(iqx)%data(:,jcol,icol)
    end IF
    IF (nq_r > 1)then
      iqx = inr
      nr(:,1,1) = current_state%zq(iqx)%data(:,jcol,icol)
      dnr(:,1,1) = current_state%sq(iqx)%data(:,jcol,icol)
    end IF
    IF (nq_r > 2)then
      iqx = i3mr
      m3r(:,1,1) = current_state%zq(iqx)%data(:,jcol,icol)
      dm3r(:,1,1) = current_state%sq(iqx)%data(:,jcol,icol)
    end IF
 
    ! Ice microphysical fields
    IF (nq_i > 0)then
      iqx = iqi
      qi(:,1,1) = current_state%zq(iqx)%data(:,jcol,icol)
      dqi(:,1,1) = current_state%sq(iqx)%data(:,jcol,icol)
    end IF
    IF (nq_s > 0)then
      iqx = iqs
      qs(:,1,1) = current_state%zq(iqx)%data(:,jcol,icol)
      dqs(:,1,1) = current_state%sq(iqx)%data(:,jcol,icol)
    end IF
    IF (nq_g > 0)then
      iqx = iqg
      qg(:,1,1) = current_state%zq(iqx)%data(:,jcol,icol)
      dqg(:,1,1) = current_state%sq(iqx)%data(:,jcol,icol)
    end IF
    IF (nq_i > 1)then
      iqx = ini
      ni(:,1,1) = current_state%zq(iqx)%data(:,jcol,icol)
      dni(:,1,1) = current_state%sq(iqx)%data(:,jcol,icol)
    end IF
    IF (nq_s > 1)then
      iqx = ins
      ns(:,1,1) = current_state%zq(iqx)%data(:,jcol,icol)
      dns(:,1,1) = current_state%sq(iqx)%data(:,jcol,icol)
    end IF
    IF (nq_g > 1)then
      iqx = ing
      ng(:,1,1) = current_state%zq(iqx)%data(:,jcol,icol)
      dng(:,1,1) = current_state%sq(iqx)%data(:,jcol,icol)
    end IF
    IF (nq_s > 2)then
      iqx = i3ms
      m3s(:,1,1) = current_state%zq(iqx)%data(:,jcol,icol)
      dm3s(:,1,1) = current_state%sq(iqx)%data(:,jcol,icol)
    end IF
    IF (nq_g > 2)then
      iqx = i3mg
      m3g(:,1,1) = current_state%zq(iqx)%data(:,jcol,icol)
      dm3g(:,1,1) = current_state%sq(iqx)%data(:,jcol,icol)
    end IF
 
    ! Aerosol fields
 
    if (i_aitkensolmass>0)     aitkensolmass(:,1,1)      = current_state%zq(i_aitkensolmass)%data(:,jcol,icol)
    if (i_aitkensolmass>0)    daitkensolmass(:,1,1)     = current_state%sq(i_aitkensolmass)%data(:,jcol,icol)
    if (i_aitkensolnumber>0)   aitkensolnumber(:,1,1)    = current_state%zq(i_aitkensolnumber)%data(:,jcol,icol)
    if (i_aitkensolnumber>0)  daitkensolnumber(:,1,1)   = current_state%sq(i_aitkensolnumber)%data(:,jcol,icol)
    if (i_accumsolmass>0)      accumsolmass(:,1,1)       = current_state%zq(i_accumsolmass)%data(:,jcol,icol)
    if (i_accumsolmass>0)     daccumsolmass(:,1,1)      = current_state%sq(i_accumsolmass)%data(:,jcol,icol)
    if (i_accumsolnumber>0)    accumsolnumber(:,1,1)     = current_state%zq(i_accumsolnumber)%data(:,jcol,icol)
    if (i_accumsolnumber>0)   daccumsolnumber(:,1,1)    = current_state%sq(i_accumsolnumber)%data(:,jcol,icol)
    if (i_coarsesolmass>0)     coarsesolmass(:,1,1)      = current_state%zq(i_coarsesolmass)%data(:,jcol,icol)
    if (i_coarsesolmass>0)    dcoarsesolmass(:,1,1)     = current_state%sq(i_coarsesolmass)%data(:,jcol,icol)
    if (i_coarsesolnumber>0)   coarsesolnumber(:,1,1)    = current_state%zq(i_coarsesolnumber)%data(:,jcol,icol)
    if (i_coarsesolnumber>0)  dcoarsesolnumber(:,1,1)   = current_state%sq(i_coarsesolnumber)%data(:,jcol,icol)
    if (i_activesolliquid>0)   activesolliquid(:,1,1)    = current_state%zq(i_activesolliquid)%data(:,jcol,icol)
    if (i_activesolliquid>0)  dactivesolliquid(:,1,1)   = current_state%sq(i_activesolliquid)%data(:,jcol,icol)
    if (i_coarsedustmass>0)    coarsedustmass(:,1,1)     = current_state%zq(i_coarsedustmass)%data(:,jcol,icol)
    if (i_coarsedustmass>0)   dcoarsedustmass(:,1,1)    = current_state%sq(i_coarsedustmass)%data(:,jcol,icol)
    if (i_coarsedustnumber>0)  coarsedustnumber(:,1,1)   = current_state%zq(i_coarsedustnumber)%data(:,jcol,icol)
    if (i_coarsedustnumber>0) dcoarsedustnumber(:,1,1)  = current_state%sq(i_coarsedustnumber)%data(:,jcol,icol)
    if (i_activeinsolice>0)    activeinsolice(:,1,1)     = current_state%zq(i_activeinsolice)%data(:,jcol,icol)
    if (i_activeinsolice>0)   dactiveinsolice(:,1,1)    = current_state%sq(i_activeinsolice)%data(:,jcol,icol)
    if (i_activesolice>0)      activesolice(:,1,1)       = current_state%zq(i_activesolice)%data(:,jcol,icol)
    if (i_activesolice>0)     dactivesolice(:,1,1)      = current_state%sq(i_activesolice)%data(:,jcol,icol)
    if (i_activeinsolliquid>0) activeinsolliquid(:,1,1)  = current_state%zq(i_activeinsolliquid)%data(:,jcol,icol)
    if (i_activeinsolliquid>0)dactiveinsolliquid(:,1,1) = current_state%sq(i_activeinsolliquid)%data(:,jcol,icol)
    if (i_accuminsolmass>0)    accuminsolmass(:,1,1)     = current_state%zq(i_accuminsolmass)%data(:,jcol,icol)
    if (i_accuminsolmass>0)   daccuminsolmass(:,1,1)    = current_state%sq(i_accuminsolmass)%data(:,jcol,icol)
    if (i_accuminsolnumber>0)  accuminsolnumber(:,1,1)   = current_state%zq(i_accuminsolnumber)%data(:,jcol,icol)
    if (i_accuminsolnumber>0) daccuminsolnumber(:,1,1)  = current_state%sq(i_accuminsolnumber)%data(:,jcol,icol)
    if (i_activesolnumber>0)   activesolnumber(:,1,1)    = current_state%zq(i_activesolnumber)%data(:,jcol,icol)
    if (i_activesolnumber>0)  dactivesolnumber(:,1,1)   = current_state%sq(i_activesolnumber)%data(:,jcol,icol)
    if (i_activeinsolnumber>0) activeinsolnumber(:,1,1)  = current_state%zq(i_activeinsolnumber)%data(:,jcol,icol)
    if (i_activeinsolnumber>0)dactiveinsolnumber(:,1,1) = current_state%sq(i_activeinsolnumber)%data(:,jcol,icol)
 
    CALL shipway_microphysics(                     &
                                ! in
       its, ite,                                   &
       jts, jte,                                   &
       kts, kte,                                   &
       dtwp,                                       &
       qv, qc, qr,                                 &
       nc, nr, m3r,                                &
       qi, qs, qg,                                 &
       ni, ns, ng,                                 &
       m3s, m3g,                                   &
       theta,                                      &
       aitkensolmass, aitkensolnumber,             &
       accumsolmass, accumsolnumber,               &
       coarsesolmass, coarsesolnumber,             &
       activesolliquid,                            &
       activesolrain,                              &
       coarsedustmass, coarsedustnumber,           &
       activeinsolice,                             &
       activesolice,                               &
       activeinsolliquid,                          &
       accuminsolmass,                             &
       accuminsolnumber,                           &
       activesolnumber,                            &
       activeinsolnumber,                          &
       aitkensolbk,                                &
       accumsolbk,                                 & 
       coarsesolbk,                                &
       exner,                                      &
       pressure, rho,                              &
       w, tke,                                     &
       z_half, z_centre,                           &
       dz,                                         &
                                ! in/out
       dqv, dqc, dqr, dnc, dnr, dm3r,              &
       dqi, dqs, dqg, dni, dns, dng, dm3s, dm3g,   &
       dth,                                        &
       daitkensolmass, daitkensolnumber,           &
       daccumsolmass, daccumsolnumber,             &
       dcoarsesolmass, dcoarsesolnumber,           &
       dactivesolliquid,                           &
       dactivesolrain,                             &
       dcoarsedustmass, dcoarsedustnumber,         &
       dactiveinsolice,                            &
       dactivesolice,                              &
       dactiveinsolliquid,                         &
       daccuminsolmass,                            &
       daccuminsolnumber,                          &
       dactivesolnumber,                           &
       dactiveinsolnumber,                         &
       ils, ile,                                   &
       jls, jle,                                   &
       kls, kle,                                   &
       l_tendency=.true.                           &
       )
 
    ! write back the tendencies
    current_state%sth%data(:,jcol,icol) = current_state%sth%data(:,jcol,icol) + dth(:,1,1)
 
    iqx = iqv
    current_state%sq(iqx)%data(:,jcol,icol) = current_state%sq(iqx)%data(:,jcol,icol) + dqv(:,1,1)
 
    ! Warm microphysical fields
    IF (nq_l > 0)then
      iqx = iql
      current_state%sq(iqx)%data(:,jcol,icol) = current_state%sq(iqx)%data(:,jcol,icol) + dqc(:,1,1)
    end IF
 
    IF (nq_r > 0)then
      iqx = iqr
      current_state%sq(iqx)%data(:,jcol,icol) = current_state%sq(iqx)%data(:,jcol,icol) + dqr(:,1,1)
    end IF
    IF (nq_l > 1)then
      iqx = inl
      current_state%sq(iqx)%data(:,jcol,icol) = current_state%sq(iqx)%data(:,jcol,icol) + dnc(:,1,1)
    end IF
    IF (nq_r > 1)then
      iqx = inr
      current_state%sq(iqx)%data(:,jcol,icol) = current_state%sq(iqx)%data(:,jcol,icol) + dnr(:,1,1)
    end IF
    IF (nq_r > 2)then
      iqx = i3mr
      current_state%sq(iqx)%data(:,jcol,icol) = current_state%sq(iqx)%data(:,jcol,icol) + dm3r(:,1,1)
    end IF
 
    ! Ice microphysical fields
    IF (nq_i > 0)then
      iqx = iqi
      current_state%sq(iqx)%data(:,jcol,icol) = current_state%sq(iqx)%data(:,jcol,icol) + dqi(:,1,1)
    end IF
    IF (nq_s > 0)then
      iqx = iqs
      current_state%sq(iqx)%data(:,jcol,icol) = current_state%sq(iqx)%data(:,jcol,icol) + dqs(:,1,1)
    end IF
    IF (nq_g > 0)then
      iqx = iqg
      current_state%sq(iqx)%data(:,jcol,icol) = current_state%sq(iqx)%data(:,jcol,icol) + dqg(:,1,1)
    end IF
    IF (nq_i > 1)then
      iqx = ini
      current_state%sq(iqx)%data(:,jcol,icol) = current_state%sq(iqx)%data(:,jcol,icol) + dni(:,1,1)
    end IF
    IF (nq_s > 1)then
      iqx = ins
      current_state%sq(iqx)%data(:,jcol,icol) = current_state%sq(iqx)%data(:,jcol,icol) + dns(:,1,1)
    end IF
    IF (nq_g > 1)then
      iqx = ing
      current_state%sq(iqx)%data(:,jcol,icol) = current_state%sq(iqx)%data(:,jcol,icol) + dng(:,1,1)
    end IF
    IF (nq_s > 2)then
      iqx = i3ms
      current_state%sq(iqx)%data(:,jcol,icol) = current_state%sq(iqx)%data(:,jcol,icol) + dm3s(:,1,1)
    end IF
    IF (nq_g > 2)then
      iqx = i3mg
      current_state%sq(iqx)%data(:,jcol,icol) = current_state%sq(iqx)%data(:,jcol,icol) + dm3g(:,1,1)
    end IF
 
    ! Aerosol fields
 
    if (i_aitkensolmass>0)     current_state%sq(i_aitkensolmass)%data(:,jcol,icol) &
       = current_state%sq(i_aitkensolmass)%data(:,jcol,icol) + daitkensolmass(:,1,1)
    if (i_aitkensolnumber>0)   current_state%sq(i_aitkensolnumber)%data(:,jcol,icol) &
       =  current_state%sq(i_aitkensolnumber)%data(:,jcol,icol) + daitkensolnumber(:,1,1)
    if (i_accumsolmass>0)      current_state%sq(i_accumsolmass)%data(:,jcol,icol) &
       =  current_state%sq(i_accumsolmass)%data(:,jcol,icol) + daccumsolmass(:,1,1)
    if (i_accumsolnumber>0)    current_state%sq(i_accumsolnumber)%data(:,jcol,icol) &
       =  current_state%sq(i_accumsolnumber)%data(:,jcol,icol) + daccumsolnumber(:,1,1)
    if (i_coarsesolmass>0)     current_state%sq(i_coarsesolmass)%data(:,jcol,icol) &
       =  current_state%sq(i_coarsesolmass)%data(:,jcol,icol) + dcoarsesolmass(:,1,1)
    if (i_coarsesolnumber>0)   current_state%sq(i_coarsesolnumber)%data(:,jcol,icol) &
       =  current_state%sq(i_coarsesolnumber)%data(:,jcol,icol) + dcoarsesolnumber(:,1,1)
    if (i_activesolliquid>0)   current_state%sq(i_activesolliquid)%data(:,jcol,icol) &
       =  current_state%sq(i_activesolliquid)%data(:,jcol,icol) + dactivesolliquid(:,1,1)
    if (i_coarsedustmass>0)    current_state%sq(i_coarsedustmass)%data(:,jcol,icol) &
       =  current_state%sq(i_coarsedustmass)%data(:,jcol,icol) + dcoarsedustmass(:,1,1)
    if (i_coarsedustnumber>0)  current_state%sq(i_coarsedustnumber)%data(:,jcol,icol) &
       =  current_state%sq(i_coarsedustnumber)%data(:,jcol,icol) + dcoarsedustnumber(:,1,1)
    if (i_activeinsolice>0)    current_state%sq(i_activeinsolice)%data(:,jcol,icol) &
       =  current_state%sq(i_activeinsolice)%data(:,jcol,icol) + dactiveinsolice(:,1,1)
    if (i_activesolice>0)      current_state%sq(i_activesolice)%data(:,jcol,icol) &
       =  current_state%sq(i_activesolice)%data(:,jcol,icol) + dactivesolice(:,1,1)
    if (i_activeinsolliquid>0) current_state%sq(i_activeinsolliquid)%data(:,jcol,icol) &
       =  current_state%sq(i_activeinsolliquid)%data(:,jcol,icol) + dactiveinsolliquid(:,1,1)
    if (i_accuminsolmass>0)    current_state%sq(i_accuminsolmass)%data(:,jcol,icol) &
       =  current_state%sq(i_accuminsolmass)%data(:,jcol,icol) + daccuminsolmass(:,1,1)
    if (i_accuminsolnumber>0)  current_state%sq(i_accuminsolnumber)%data(:,jcol,icol) &
       =  current_state%sq(i_accuminsolnumber)%data(:,jcol,icol) + daccuminsolnumber(:,1,1)
    if (i_activesolnumber>0)   current_state%sq(i_activesolnumber)%data(:,jcol,icol) &
       =  current_state%sq(i_activesolnumber)%data(:,jcol,icol) + dactivesolnumber(:,1,1)
    if (i_activeinsolnumber>0) current_state%sq(i_activeinsolnumber)%data(:,jcol,icol) &
       =  current_state%sq(i_activeinsolnumber)%data(:,jcol,icol) + dactiveinsolnumber(:,1,1)
 
    ! for total surface precipitation, sum the surface rain rate (cloud + rain which is precip_r)
    ! 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
       surface_precip(target_y_index,target_x_index) = &
            casdiags % SurfaceRainR(1,1)
    else
       surface_precip(target_y_index,target_x_index) = &
            casdiags % SurfaceRainR(1,1) + casdiags % SurfaceSnowR(1,1)
    endif
    call populate_casim_monc_dg(current_state, casdiags)
 
    
  end subroutine timestep_callback
 
 
  !! Reads the casim configuration 
  !! @param current_state The current model state
  subroutine read_configuration(current_state)
 
 
    Use mphys_parameters, only: p1, p2, p3, sp1, sp2, sp3   
 
    type(model_state_type), target, intent(inout) :: current_state
 
    integer :: ierr
 
    option          = options_get_integer(current_state%options_database, 'option')
    diag_mu_option  = options_get_integer(current_state%options_database, 'diag_mu_option')
    iopt_act        = options_get_integer(current_state%options_database, 'iopt_act')
    iopt_inuc       = options_get_integer(current_state%options_database, 'iopt_inuc')
    process_level   = options_get_integer(current_state%options_database, 'process_level')
    aerosol_option  = options_get_integer(current_state%options_database, 'aerosol_option')
    max_step_length = options_get_real(current_state%options_database, 'max_step_length')
    max_sed_length  = options_get_real(current_state%options_database, 'max_sed_length')
    p1              = options_get_real(current_state%options_database, 'p1')
    p2              = options_get_real(current_state%options_database, 'p2')
    p3              = options_get_real(current_state%options_database, 'p3')
    sp1             = options_get_real(current_state%options_database, 'sp1')
    sp2             = options_get_real(current_state%options_database, 'sp2')
    sp3             = options_get_real(current_state%options_database, 'sp3')
    max_mu          = options_get_real(current_state%options_database, 'max_mu')
    fix_mu          = options_get_real(current_state%options_database, 'fix_mu')
 
    l_aaut          = options_get_logical(current_state%options_database, 'l_aaut')
    l_aacc          = options_get_logical(current_state%options_database, 'l_aacc')
    l_aevp          = options_get_logical(current_state%options_database, 'l_aevp')
    l_ased          = options_get_logical(current_state%options_database, 'l_ased')
    l_warm          = options_get_logical(current_state%options_database, 'l_warm')
    l_inuc          = options_get_logical(current_state%options_database, 'l_inuc')
    l_iaut          = options_get_logical(current_state%options_database, 'l_iaut')
    l_idep          = options_get_logical(current_state%options_database, 'l_idep')
    l_iacw          = options_get_logical(current_state%options_database, 'l_iacw')
    l_active_inarg2000 = options_get_logical(current_state%options_database, 'l_active_inarg2000')
    l_separate_rain = options_get_logical(current_state%options_database, 'l_separate_rain')
    l_sg            = options_get_logical(current_state%options_database, 'l_sg')
    l_g             = options_get_logical(current_state%options_database, 'l_g')
    l_passive       = options_get_logical(current_state%options_database, 'l_passive')
    l_passive3m     = options_get_logical(current_state%options_database, 'l_passive3m')
    l_limit_psd     = options_get_logical(current_state%options_database, 'l_limit_psd')
    l_override_checks = options_get_logical(current_state%options_database, 'l_override_checks')
    l_raci_g        = options_get_logical(current_state%options_database, 'l_raci_g')
    l_onlycollect   = options_get_logical(current_state%options_database, 'l_onlycollect')
    l_abelshipway   = options_get_logical(current_state%options_database, 'l_abelshipway')
    l_cons          = options_get_logical(current_state%options_database, 'l_cons')
    l_rain          = options_get_logical(current_state%options_database, 'l_rain')
    l_sed_3mdiff    = options_get_logical(current_state%options_database, 'l_sed_3mdiff')
    l_sed_icecloud_as_1m = options_get_logical(current_state%options_database, 'l_sed_icecloud_as_1m')
    l_tidy_conserve_e = options_get_logical(current_state%options_database, 'l_tidy_conserve_E')
    l_tidy_conserve_q = options_get_logical(current_state%options_database, 'l_tidy_conserve_q')
 
    l_inhom_revp    = options_get_logical(current_state%options_database, 'l_inhom_revp')
    l_pcond         = options_get_logical(current_state%options_database, 'l_pcond')
    l_praut         = options_get_logical(current_state%options_database, 'l_praut')
    l_pracw         = options_get_logical(current_state%options_database, 'l_pracw')
    l_pracr         = options_get_logical(current_state%options_database, 'l_pracr')
    l_prevp         = options_get_logical(current_state%options_database, 'l_prevp')
    l_psedl         = options_get_logical(current_state%options_database, 'l_psedl')
    l_psedr         = options_get_logical(current_state%options_database, 'l_psedr')
    l_ptidy         = options_get_logical(current_state%options_database, 'l_ptidy')
    l_ptidy2        = options_get_logical(current_state%options_database, 'l_ptidy2')
    l_pinuc         = options_get_logical(current_state%options_database, 'l_pinuc')
    l_pidep         = options_get_logical(current_state%options_database, 'l_pidep')
    l_piacw         = options_get_logical(current_state%options_database, 'l_piacw')
    l_psaut         = options_get_logical(current_state%options_database, 'l_psaut')
    l_psdep         = options_get_logical(current_state%options_database, 'l_psdep')
    l_psacw         = options_get_logical(current_state%options_database, 'l_psacw')
    l_pgdep         = options_get_logical(current_state%options_database, 'l_pgdep')
    l_pseds         = options_get_logical(current_state%options_database, 'l_pseds')
    l_psedi         = options_get_logical(current_state%options_database, 'l_psedi')
    l_psedg         = options_get_logical(current_state%options_database, 'l_psedg')
    l_psaci         = options_get_logical(current_state%options_database, 'l_psaci')
    l_praci         = options_get_logical(current_state%options_database, 'l_praci')
    l_psacr         = options_get_logical(current_state%options_database, 'l_psacr')
    l_pgacr         = options_get_logical(current_state%options_database, 'l_pgacr')
    l_pgacw         = options_get_logical(current_state%options_database, 'l_pgacw')
    l_pgaci         = options_get_logical(current_state%options_database, 'l_pgaci')
    l_pgacs         = options_get_logical(current_state%options_database, 'l_pgacs')
    l_piagg         = options_get_logical(current_state%options_database, 'l_piagg')
    l_psagg         = options_get_logical(current_state%options_database, 'l_psagg')
    l_pgagg         = options_get_logical(current_state%options_database, 'l_pgagg')
    l_psbrk         = options_get_logical(current_state%options_database, 'l_psbrk')
    l_pgshd         = options_get_logical(current_state%options_database, 'l_pgshd')
    l_pihal         = options_get_logical(current_state%options_database, 'l_pihal')
    l_psmlt         = options_get_logical(current_state%options_database, 'l_psmlt')
    l_pgmlt         = options_get_logical(current_state%options_database, 'l_pgmlt')
    l_phomr         = options_get_logical(current_state%options_database, 'l_phomr')
    l_phomc         = options_get_logical(current_state%options_database, 'l_phomc')
    l_pssub         = options_get_logical(current_state%options_database, 'l_pssub')
    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')
 
  end subroutine read_configuration
 
  subroutine field_information_retrieval_callback(current_state, name, field_information)
    type(model_state_type), target, intent(inout) :: current_state
    character(len=*), intent(in) :: name
    type(component_field_information_type), intent(out) :: field_information
 
    field_information%field_type=component_array_field_type
    field_information%data_type=component_double_data_type
    if (name .eq. "surface_precip") then
       field_information%number_dimensions=2
       field_information%dimension_sizes(1)=current_state%local_grid%size(y_index)
       field_information%dimension_sizes(2)=current_state%local_grid%size(x_index)
    !else if (name .eq. "pcond_total" .or. name .eq. "psedl_total") then
    !   field_information%number_dimensions=1
    !   field_information%dimension_sizes(1)=current_state%local_grid%size(Z_INDEX)
    else
       field_information%number_dimensions=3
       field_information%dimension_sizes(1)=current_state%local_grid%size(z_index)
       field_information%dimension_sizes(2)=current_state%local_grid%size(y_index)
       field_information%dimension_sizes(3)=current_state%local_grid%size(x_index)
    endif
       
    field_information%enabled=.true.    
 
  end subroutine field_information_retrieval_callback
 
  subroutine field_value_retrieval_callback(current_state, name, field_value)
    type(model_state_type), target, intent(inout) :: current_state
    character(len=*), intent(in) :: name
    type(component_field_value_type), intent(out) :: field_value
    
    integer :: i
 
    if (name .eq. "surface_precip") then
       allocate(field_value%real_2d_array(current_state%local_grid%size(y_index), &
            current_state%local_grid%size(x_index)))
       field_value%real_2d_array(:,:)= surface_precip(:,:)
    else if (name .eq. "condensation_rate") then
       allocate(field_value%real_3d_array(current_state%local_grid%size(z_index),  &
            current_state%local_grid%size(y_index),                                &
            current_state%local_grid%size(x_index)))
       field_value%real_3d_array(:,:,:) = casim_monc_dgs % pcond(:,:,:)
!!$    else if (name .eq. "pcond_total") then
!!$       allocate(field_value%real_1d_array(current_state%local_grid%size(Z_INDEX)))
!!$       field_value%real_1d_array(:)=pcond_tot(:)
    end if
    
  end subroutine field_value_retrieval_callback
 
end module casim_mod