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

  use micro_main, only: shipway_microphysics
  use generic_diagnostic_variables, ONLY: casdiags, allocate_diagnostic_space, &
       deallocate_diagnostic_space

  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)


  ! 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

  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
    
    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(1))

    casim_get_descriptor%published_fields(1)="surface_precip" 

  end function casim_get_descriptor

  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
    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)

    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(:,:)
    end if
    
  end subroutine field_value_retrieval_callback

  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 
    
    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(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))

    casdiags % l_surface_rain = .true.
    casdiags % l_surface_snow = .true.
    casdiags % l_surface_graup = .true.

    CALL allocate_diagnostic_space(its, ite, jts, jte, kts, kte)

    allocate(surface_precip(y_size_local, x_size_local))

    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%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%cq(iqi) = -1.0
      end if
      if (nq_s>0)then
        iqs = get_q_index(standard_q_names%SNOW_MASS, 'casim')
        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%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') 

  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%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

    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,                          &
       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). 
    surface_precip(target_y_index,target_x_index) = &
         casdiags % SurfaceRainR(1,1) + casdiags % SurfaceSnowR(1,1) 
    
  end subroutine timestep_callback


  !! Reads the casim configuration 
  !! @param current_state The current model state
  subroutine read_configuration(current_state)


    Use mphys_switches, only: iopt_act, option, aerosol_option        &
       , 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 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

end module casim_mod