simplecloud.F90

Go to the documentation of this file.

 
module simplecloud_mod
  use datadefn_mod, only : default_precision
  use monc_component_mod, only : component_array_field_type, &
      component_double_data_type, component_descriptor_type, &
      component_field_value_type, component_field_information_type
  use optionsdatabase_mod, only : options_get_real, options_get_integer
  use state_mod, only : forward_stepping, model_state_type
  use grids_mod, only : z_index, x_index, y_index
  use science_constants_mod, only : r_over_cp, rlvap_over_cp
  use saturation_mod, only: qsaturation, dqwsatdt
  use q_indices_mod, only: get_q_index, standard_q_names
  use registry_mod, only : is_component_enabled
  use logging_mod, only : log_error, log_master_log
 
implicit none
 
#ifndef TEST_MODE
  private
#endif
    
  ! Indices for vapour and cloud
  integer :: iqv, iql
 
  integer :: k_cloudmax ! max k index for height
  real(kind=default_precision) :: max_height_cloud
 
  ! Local tendency diagnostic variables for this component
  ! 3D tendency fields and logicals for their use
  real(kind=default_precision), dimension(:,:,:), allocatable ::     &
                   tend_3d_th,      tend_3d_qv,      tend_3d_ql,      tend_3d_tabs
  logical ::     l_tend_3d_th,    l_tend_3d_qv,    l_tend_3d_ql,    l_tend_3d_tabs
  ! Local mean tendency profile fields and logicals for their use
  real(kind=default_precision), dimension(:), allocatable ::         &
               tend_pr_tot_th,  tend_pr_tot_qv,  tend_pr_tot_ql,  tend_pr_tot_tabs
  logical :: l_tend_pr_tot_th,l_tend_pr_tot_qv,l_tend_pr_tot_ql,l_tend_pr_tot_tabs
  integer :: diagnostic_generation_frequency
 
 
  public simplecloud_get_descriptor
 
contains
 
  type(component_descriptor_type) function simplecloud_get_descriptor()
    simplecloud_get_descriptor%name="simplecloud"
    simplecloud_get_descriptor%version=0.1
    simplecloud_get_descriptor%initialisation=>initialisation_callback
    simplecloud_get_descriptor%timestep=>timestep_callback
    simplecloud_get_descriptor%finalisation=>finalisation_callback
 
    simplecloud_get_descriptor%field_value_retrieval=>field_value_retrieval_callback
    simplecloud_get_descriptor%field_information_retrieval=>field_information_retrieval_callback
    allocate(simplecloud_get_descriptor%published_fields(4+4))
 
    simplecloud_get_descriptor%published_fields(1)="tend_th_simplecloud_3d_local"
    simplecloud_get_descriptor%published_fields(2)="tend_qv_simplecloud_3d_local"
    simplecloud_get_descriptor%published_fields(3)="tend_ql_simplecloud_3d_local"
    simplecloud_get_descriptor%published_fields(4)="tend_tabs_simplecloud_3d_local"
 
    simplecloud_get_descriptor%published_fields(4+1)="tend_th_simplecloud_profile_total_local"
    simplecloud_get_descriptor%published_fields(4+2)="tend_qv_simplecloud_profile_total_local"
    simplecloud_get_descriptor%published_fields(4+3)="tend_ql_simplecloud_profile_total_local"
    simplecloud_get_descriptor%published_fields(4+4)="tend_tabs_simplecloud_profile_total_local"
 
  end function simplecloud_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
    integer :: strcomp
 
    ! Field information for 3d
    strcomp=index(name, "simplecloud_3d_local")
    if (strcomp .ne. 0) then
      field_information%field_type=component_array_field_type
      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)
      field_information%data_type=component_double_data_type
 
      if      (name .eq. "tend_th_simplecloud_3d_local") then
        field_information%enabled=l_tend_3d_th
      else if (name .eq. "tend_qv_simplecloud_3d_local") then
        field_information%enabled=l_tend_3d_qv
      else if (name .eq. "tend_ql_simplecloud_3d_local") then
        field_information%enabled=l_tend_3d_ql
      else if (name .eq. "tend_tabs_simplecloud_3d_local") then
        field_information%enabled=l_tend_3d_tabs
      else
        field_information%enabled=.true.
      end if
 
    end if !end 3d check
 
    ! Field information for profiles
    strcomp=index(name, "simplecloud_profile_total_local")
    if (strcomp .ne. 0) then
      field_information%field_type=component_array_field_type
      field_information%number_dimensions=1
      field_information%dimension_sizes(1)=current_state%local_grid%size(z_index)
      field_information%data_type=component_double_data_type
 
      if      (name .eq. "tend_th_simplecloud_profile_total_local") then
        field_information%enabled=l_tend_pr_tot_th
      else if (name .eq. "tend_qv_simplecloud_profile_total_local") then
        field_information%enabled=l_tend_pr_tot_qv
      else if (name .eq. "tend_ql_simplecloud_profile_total_local") then
        field_information%enabled=l_tend_pr_tot_ql
      else if (name .eq. "tend_tabs_simplecloud_profile_total_local") then
        field_information%enabled=l_tend_pr_tot_tabs
      else
        field_information%enabled=.true.
      end if
 
    end if !end profile check
 
  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
 
    ! 3d Tendency Fields
    if      (name .eq. "tend_th_simplecloud_3d_local" .and. allocated(tend_3d_th)) then
      call set_published_field_value(field_value, real_3d_field=tend_3d_th)
    else if (name .eq. "tend_qv_simplecloud_3d_local" .and. allocated(tend_3d_qv)) then
      call set_published_field_value(field_value, real_3d_field=tend_3d_qv)
    else if (name .eq. "tend_ql_simplecloud_3d_local" .and. allocated(tend_3d_ql)) then
      call set_published_field_value(field_value, real_3d_field=tend_3d_ql)
    else if (name .eq. "tend_tabs_simplecloud_3d_local" .and. allocated(tend_3d_tabs)) then
      call set_published_field_value(field_value, real_3d_field=tend_3d_tabs)
 
    ! Profile Tendency Fields
    else if (name .eq. "tend_th_simplecloud_profile_total_local" .and. allocated(tend_pr_tot_th)) then
      call set_published_field_value(field_value, real_1d_field=tend_pr_tot_th)
    else if (name .eq. "tend_qv_simplecloud_profile_total_local" .and. allocated(tend_pr_tot_qv)) then
      call set_published_field_value(field_value, real_1d_field=tend_pr_tot_qv)
    else if (name .eq. "tend_ql_simplecloud_profile_total_local" .and. allocated(tend_pr_tot_ql)) then
      call set_published_field_value(field_value, real_1d_field=tend_pr_tot_ql)
    else if (name .eq. "tend_tabs_simplecloud_profile_total_local" .and. allocated(tend_pr_tot_tabs)) then
      call set_published_field_value(field_value, real_1d_field=tend_pr_tot_tabs)
    end if
 
  end subroutine field_value_retrieval_callback
 
 
  subroutine initialisation_callback(current_state)
    type(model_state_type), target, intent(inout) :: current_state
 
    integer :: k ! look counter
    logical :: l_qdiag
 
    if (is_component_enabled(current_state%options_database, "casim")) then
      call log_master_log(log_error, "Casim and Simplecloud are enabled, this does not work yet. Please disable one")
    end if 
 
    iqv=get_q_index(standard_q_names%VAPOUR, 'simplecloud')
    iql=get_q_index(standard_q_names%CLOUD_LIQUID_MASS, 'simplecloud')
 
    ! set buoyancy coefficient (value for vapour should be set
    ! elsewhere for a moist model
    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
    current_state%cq(iql) = -1.0 
 
    max_height_cloud=options_get_real(current_state%options_database, "max_height_cloud")
    do k=2, current_state%local_grid%size(z_index)-1
      if (current_state%global_grid%configuration%vertical%zn(k) > max_height_cloud) exit
    end do
    k_cloudmax=k-1
 
    ! Set tendency diagnostic logicals based on availability
    ! Need to use 3d tendencies to compute the profiles, so they will be allocated
    !      in the case where profiles are available
    l_qdiag =  (.not. current_state%passive_q .and. current_state%number_q_fields .gt. 0)
 
    l_tend_pr_tot_th  = current_state%th%active
    l_tend_pr_tot_qv  = l_qdiag .and. current_state%number_q_fields .ge. 1
    l_tend_pr_tot_ql  = l_qdiag .and. current_state%number_q_fields .ge. 2
    l_tend_pr_tot_tabs  = l_tend_pr_tot_th
 
    l_tend_3d_th  = current_state%th%active .or. l_tend_pr_tot_th
    l_tend_3d_qv  = (l_qdiag .and. current_state%number_q_fields .ge. 1) .or. l_tend_pr_tot_qv
    l_tend_3d_ql  = (l_qdiag .and. current_state%number_q_fields .ge. 2) .or. l_tend_pr_tot_ql
    l_tend_3d_tabs = l_tend_3d_th
 
    ! Allocate 3d tendency fields upon availability
    if (l_tend_3d_th) then
      allocate( tend_3d_th(current_state%local_grid%size(z_index),  &
                           current_state%local_grid%size(y_index),  &
                           current_state%local_grid%size(x_index)   )    )
    endif
    if (l_tend_3d_qv) then
      allocate( tend_3d_qv(current_state%local_grid%size(z_index),  &
                           current_state%local_grid%size(y_index),  &
                           current_state%local_grid%size(x_index)   )    )
    endif
    if (l_tend_3d_ql) then
      allocate( tend_3d_ql(current_state%local_grid%size(z_index),  &
                           current_state%local_grid%size(y_index),  &
                           current_state%local_grid%size(x_index)   )    )
    endif
    if (l_tend_3d_tabs) then
      allocate( tend_3d_tabs(current_state%local_grid%size(z_index),  &
                             current_state%local_grid%size(y_index),  &
                             current_state%local_grid%size(x_index)   )    )
    endif
 
    ! Allocate profile tendency fields upon availability
    if (l_tend_pr_tot_th) then
      allocate( tend_pr_tot_th(current_state%local_grid%size(z_index)) )
    endif
    if (l_tend_pr_tot_qv) then
      allocate( tend_pr_tot_qv(current_state%local_grid%size(z_index)) )
    endif
    if (l_tend_pr_tot_ql) then
      allocate( tend_pr_tot_ql(current_state%local_grid%size(z_index)) )
    endif
    if (l_tend_pr_tot_tabs) then
      allocate( tend_pr_tot_tabs(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  
 
 
  subroutine finalisation_callback(current_state)
    type(model_state_type), target, intent(inout) :: current_state
 
    if (allocated(tend_3d_th)) deallocate(tend_3d_th)
    if (allocated(tend_3d_qv)) deallocate(tend_3d_qv)
    if (allocated(tend_3d_ql)) deallocate(tend_3d_ql)
    if (allocated(tend_3d_tabs)) deallocate(tend_3d_tabs)
 
    if (allocated(tend_pr_tot_th)) deallocate(tend_pr_tot_th)
    if (allocated(tend_pr_tot_qv)) deallocate(tend_pr_tot_qv)
    if (allocated(tend_pr_tot_ql)) deallocate(tend_pr_tot_ql)
    if (allocated(tend_pr_tot_tabs)) deallocate(tend_pr_tot_tabs)
 
  end subroutine finalisation_callback
 
 
  subroutine timestep_callback(current_state)
    type(model_state_type), target, intent(inout) :: current_state
 
    real(DEFAULT_PRECISION) :: TdegK   ! Temperature in Kelvin
    real(DEFAULT_PRECISION) :: Pmb     ! Pressure in mb
    real(DEFAULT_PRECISION) :: exner   ! Exner pressure
    real(DEFAULT_PRECISION) :: one_over_exner ! Reciprocal of Exner pressure
    real(DEFAULT_PRECISION) :: qv,qc   ! Shorthand for vapour and cloud mass mixing ratio
    real(DEFAULT_PRECISION) :: qs      ! Saturation mixing ratio
    real(DEFAULT_PRECISION) :: dqsdT   ! Rate of change of qs with temperature
    real(DEFAULT_PRECISION) :: qsatfac ! Multiplicative factor
    real(DEFAULT_PRECISION) :: dmass   ! Mass transfer mixing ratio
    
    integer :: k          ! Loop counter
    integer :: icol, jcol ! Shorthand column indices
 
    real(DEFAULT_PRECISION) :: dtm  ! Local timestep variable
    integer :: target_x_index, target_y_index
 
    ! Zero profile tendency totals on first instance in the sum
    if (current_state%first_timestep_column) then
      if (l_tend_pr_tot_th) then
        tend_pr_tot_th(:)=0.0_default_precision
      endif
      if (l_tend_pr_tot_qv) then
        tend_pr_tot_qv(:)=0.0_default_precision
      endif
      if (l_tend_pr_tot_ql) then
        tend_pr_tot_ql(:)=0.0_default_precision
      endif
      if (l_tend_pr_tot_tabs) then
        tend_pr_tot_tabs(:)=0.0_default_precision
      endif
    endif  ! zero totals
 
 
    if (current_state%halo_column) return
 
 
    dtm = current_state%dtm*2.0
    if (current_state%field_stepping == forward_stepping) dtm=current_state%dtm! Should this be revised to scalar_stepping
 
    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 (mod(current_state%timestep, diagnostic_generation_frequency) == 0) then
      call save_precomponent_tendencies(current_state, icol, jcol, target_x_index, target_y_index)
    end if
 
    do k=2,k_cloudmax
 
      exner = current_state%global_grid%configuration%vertical%rprefrcp(k)
      one_over_exner = current_state%global_grid%configuration%vertical%prefrcp(k)
      pmb   = (current_state%global_grid%configuration%vertical%prefn(k)/100.) 
 
      if (current_state%field_stepping == forward_stepping) then  ! Should this be revised to scalar_stepping
        qv    = current_state%q(iqv)%data(k, jcol, icol) + current_state%sq(iqv)%data(k, jcol, icol)*dtm
        qc    = current_state%q(iql)%data(k, jcol, icol) + current_state%sq(iql)%data(k, jcol, icol)*dtm
        tdegk = (current_state%th%data(k, jcol, icol) + current_state%sth%data(k, jcol, icol)*dtm &
           + current_state%global_grid%configuration%vertical%thref(k))*exner
      else
        qv    = current_state%zq(iqv)%data(k, jcol, icol) + current_state%sq(iqv)%data(k, jcol, icol)*dtm
        qc    = current_state%zq(iql)%data(k, jcol, icol) + current_state%sq(iql)%data(k, jcol, icol)*dtm
        tdegk = (current_state%zth%data(k, jcol, icol) + current_state%sth%data(k, jcol, icol)*dtm &
           + current_state%global_grid%configuration%vertical%thref(k))*exner
      end if
      ! Calculate the cloud/vapour increments
 
      qs = qsaturation(tdegk, pmb)
 
      if (qv > qs .or. qc >0.0)then
        dqsdt = dqwsatdt(qs, tdegk)
        
        qsatfac = 1.0/(1.0 + rlvap_over_cp*dqsdt)
      
        dmass = max(-qc,(qv-qs)*qsatfac)/dtm
      
        current_state%sq(iqv)%data(k, jcol, icol) = current_state%sq(iqv)%data(k, jcol, icol) - dmass
        current_state%sq(iql)%data(k, jcol, icol) = current_state%sq(iql)%data(k, jcol, icol) + dmass
 
        current_state%sth%data(k, jcol, icol) = current_state%sth%data(k, jcol, icol) &
           + rlvap_over_cp*dmass*one_over_exner
 
      end if
 
    end do
 
    ! If there's any cloud above then evaporate it
    do k=k_cloudmax+1, current_state%local_grid%size(z_index)
      if (current_state%scalar_stepping == forward_stepping) then
        qv    = current_state%q(iqv)%data(k, jcol, icol) + current_state%sq(iqv)%data(k, jcol, icol)*dtm
        qc    = current_state%q(iql)%data(k, jcol, icol) + current_state%sq(iql)%data(k, jcol, icol)*dtm
        tdegk = (current_state%th%data(k, jcol, icol) + current_state%sth%data(k, jcol, icol)*dtm &
           + current_state%global_grid%configuration%vertical%thref(k))*exner
      else
        qv    = current_state%zq(iqv)%data(k, jcol, icol) + current_state%sq(iqv)%data(k, jcol, icol)*dtm
        qc    = current_state%zq(iql)%data(k, jcol, icol) + current_state%sq(iql)%data(k, jcol, icol)*dtm
        tdegk = (current_state%zth%data(k, jcol, icol) + current_state%sth%data(k, jcol, icol)*dtm &
           + current_state%global_grid%configuration%vertical%thref(k))*exner
      end if
      if (qc >0.0)then
        dmass = -qc/dtm
      
        current_state%sq(iqv)%data(k, jcol, icol) = current_state%sq(iqv)%data(k, jcol, icol) - dmass
        current_state%sq(iql)%data(k, jcol, icol) = current_state%sq(iql)%data(k, jcol, icol) + dmass
 
        current_state%sth%data(k, jcol, icol) = current_state%sth%data(k, jcol, icol) &
           + rlvap_over_cp*dmass*one_over_exner
 
      end if
    end do
 
    if (mod(current_state%timestep, diagnostic_generation_frequency) == 0) then
      call compute_component_tendencies(current_state, icol, jcol, target_x_index, target_y_index)
    end if
 
  end subroutine timestep_callback
 
 
  subroutine save_precomponent_tendencies(current_state, cxn, cyn, txn, tyn)
    type(model_state_type), target, intent(in) :: current_state
    integer, intent(in) ::  cxn, cyn, txn, tyn
 
    ! Save 3d tendency fields upon request (of 3d or profiles) and availability
    if (l_tend_3d_th) then
      tend_3d_th(:,tyn,txn)=current_state%sth%data(:,cyn,cxn)
    endif
    if (l_tend_3d_qv) then
      tend_3d_qv(:,tyn,txn)=current_state%sq(iqv)%data(:,cyn,cxn)
    endif
    if (l_tend_3d_ql) then
      tend_3d_ql(:,tyn,txn)=current_state%sq(iql)%data(:,cyn,cxn)
    endif
    if (l_tend_3d_tabs) then
      tend_3d_tabs(:,tyn,txn)=current_state%sth%data(:,cyn,cxn) * current_state%global_grid%configuration%vertical%rprefrcp(:)
    endif
 
  end subroutine save_precomponent_tendencies
 
 
  subroutine compute_component_tendencies(current_state, cxn, cyn, txn, tyn)
    type(model_state_type), target, intent(inout) :: current_state
    integer, intent(in) ::  cxn, cyn, txn, tyn
 
    ! Calculate change in tendency due to component
    if (l_tend_3d_th) then
      tend_3d_th(:,tyn,txn)=current_state%sth%data(:,cyn,cxn)     - tend_3d_th(:,tyn,txn)
    endif
    if (l_tend_3d_qv) then
      tend_3d_qv(:,tyn,txn)=current_state%sq(iqv)%data(:,cyn,cxn) - tend_3d_qv(:,tyn,txn)
    endif
    if (l_tend_3d_ql) then
      tend_3d_ql(:,tyn,txn)=current_state%sq(iql)%data(:,cyn,cxn) - tend_3d_ql(:,tyn,txn)
    endif
    if (l_tend_3d_tabs) then
      tend_3d_tabs(:,tyn,txn)=   &
       current_state%sth%data(:,cyn,cxn) * current_state%global_grid%configuration%vertical%rprefrcp(:)   &
        - tend_3d_tabs(:,tyn,txn)
    endif
 
   ! Add local tendency fields to the profile total
    if (l_tend_pr_tot_th) then
      tend_pr_tot_th(:)=tend_pr_tot_th(:) + tend_3d_th(:,tyn,txn)
    endif
    if (l_tend_pr_tot_qv) then
      tend_pr_tot_qv(:)=tend_pr_tot_qv(:) + tend_3d_qv(:,tyn,txn)
    endif
    if (l_tend_pr_tot_ql) then
      tend_pr_tot_ql(:)=tend_pr_tot_ql(:) + tend_3d_ql(:,tyn,txn)
    endif
    if (l_tend_pr_tot_tabs) then
      tend_pr_tot_tabs(:)=tend_pr_tot_tabs(:) + tend_3d_tabs(:,tyn,txn)
    endif
 
  end subroutine compute_component_tendencies
 
 
  subroutine set_published_field_value(field_value, real_1d_field, real_2d_field, real_3d_field)
    type(component_field_value_type), intent(inout) :: field_value
    real(kind=default_precision), dimension(:), optional :: real_1d_field
    real(kind=default_precision), dimension(:,:), optional :: real_2d_field
    real(kind=default_precision), dimension(:,:,:), optional :: real_3d_field
 
    if (present(real_1d_field)) then
      allocate(field_value%real_1d_array(size(real_1d_field)), source=real_1d_field)
    else if (present(real_2d_field)) then
      allocate(field_value%real_2d_array(size(real_2d_field, 1), size(real_2d_field, 2)), source=real_2d_field)
    else if (present(real_3d_field)) then
      allocate(field_value%real_3d_array(size(real_3d_field, 1), size(real_3d_field, 2), size(real_3d_field, 3)), &
               source=real_3d_field)
    end if
  end subroutine set_published_field_value
 
end module simplecloud_mod