Add the creation of a TUV-x photolysis calculator to the MUSICA wrapper

musica/micm/micm_util.F90

@@ -0,0 +1,159 @@
+module micm_util
+  implicit none
+
+  private
+  public :: reshape_into_micm_arr, reshape_into_ccpp_arr, convert_to_mol_per_cubic_meter, &
+            convert_to_mass_mixing_ratio
+
+contains
+
+  subroutine reshape_into_micm_arr(temperature, pressure, dry_air_density, constituents, &
+      rate_params, m_temperature, m_pressure, m_dry_air_density, m_constituents, m_rate_params)
+    use iso_c_binding, only: c_double
+    use ccpp_kinds,    only: kind_phys
+
+    real(kind_phys), target, intent(in)  :: temperature(:,:)     ! K
+    real(kind_phys), target, intent(in)  :: pressure(:,:)        ! Pa
+    real(kind_phys), target, intent(in)  :: dry_air_density(:,:) ! kg m-3
+    real(kind_phys), target, intent(in)  :: constituents(:,:,:)  ! kg kg-1
+    real(kind_phys), target, intent(in)  :: rate_params(:,:,:)
+    real(c_double),  target, intent(out) :: m_temperature(:)     ! K
+    real(c_double),  target, intent(out) :: m_pressure(:)        ! Pa
+    real(c_double),  target, intent(out) :: m_dry_air_density(:) ! kg m-3
+    real(c_double),  target, intent(out) :: m_constituents(:)    ! kg kg-1
+    real(c_double),  target, intent(out) :: m_rate_params(:)
+
+    ! local variables
+    integer :: num_columns, num_layers
+    integer :: num_constituents, num_rate_params
+    integer :: i_column, i_layer, i_elem, i_constituents, i_rate_params
+
+    num_columns = size(constituents, dim=1)
+    num_layers = size(constituents, dim=2)
+    num_constituents = size(constituents, dim=3)
+    num_rate_params = size(rate_params, dim=3)
+
+    ! Reshape into 1-D arry in species-column first order
+    ! refers to: state.variables_[i_cell][i_species] = concentrations[i_species_elem++]
+    i_elem = 1
+    i_constituents = 1
+    i_rate_params = 1
+    do i_layer = 1, num_layers
+      do i_column = 1, num_columns
+        m_temperature(i_elem) = real(temperature(i_column, i_layer), c_double)
+        m_pressure(i_elem) = real(pressure(i_column, i_layer), c_double)
+        m_dry_air_density(i_elem) = real(dry_air_density(i_column, i_layer), c_double)
+        m_constituents(i_constituents : i_constituents + num_constituents - 1) &
+                            = real(constituents(i_column, i_layer, :), c_double)
+        m_rate_params(i_rate_params : i_rate_params + num_rate_params - 1) &
+                            = real(rate_params(i_column, i_layer, :), c_double)
+        i_elem = i_elem + 1
+        i_constituents = i_constituents + num_constituents
+        i_rate_params = i_rate_params + num_rate_params
+      end do
+    end do
+
+  end subroutine reshape_into_micm_arr
+
+  subroutine reshape_into_ccpp_arr(temperature, pressure, dry_air_density, constituents, &
+      rate_params, m_temperature, m_pressure, m_dry_air_density, m_constituents, m_rate_params)
+    use iso_c_binding, only: c_double
+    use ccpp_kinds,    only: kind_phys
+    real(kind_phys), intent(out) :: temperature(:,:)     ! K
+    real(kind_phys), intent(out) :: pressure(:,:)        ! Pa
+    real(kind_phys), intent(out) :: dry_air_density(:,:) ! kg m-3
+    real(kind_phys), intent(out) :: constituents(:,:,:)  ! kg kg-1
+    real(kind_phys), intent(out) :: rate_params(:,:,:)
+    real(c_double),  intent(in)  :: m_temperature(:)     ! K
+    real(c_double),  intent(in)  :: m_pressure(:)        ! Pa
+    real(c_double),  intent(in)  :: m_dry_air_density(:) ! kg m-3
+    real(c_double),  intent(in)  :: m_constituents(:)    ! kg kg-1
+    real(c_double),  intent(in)  :: m_rate_params(:)
+
+    ! local variables
+    integer :: num_columns, num_layers
+    integer :: num_constituents, num_rate_params
+    integer :: i_column, i_layer, i_elem, i_constituents, i_rate_params
+
+    num_columns = size(constituents, dim=1)
+    num_layers = size(constituents, dim=2)
+    num_constituents = size(constituents, dim=3)
+    num_rate_params = size(rate_params, dim=3)
+
+    i_elem = 1
+    i_constituents = 1
+    i_rate_params = 1
+    do i_layer = 1, num_layers
+      do i_column = 1, num_columns
+        temperature(i_column, i_layer) = real(m_temperature(i_elem), kind_phys)
+        pressure(i_column, i_layer) = real(m_pressure(i_elem), kind_phys)
+        dry_air_density(i_column, i_layer) = real(m_dry_air_density(i_elem), kind_phys)
+        constituents(i_column, i_layer, :) &
+            = real(m_constituents(i_constituents : i_constituents + num_constituents - 1), kind_phys)
+        rate_params(i_column, i_layer, :) &
+            = real(m_rate_params(i_rate_params : i_rate_params + num_rate_params - 1), kind_phys)
+        i_elem = i_elem + 1
+        i_constituents = i_constituents + num_constituents
+        i_rate_params = i_rate_params + num_rate_params
+      end do
+    end do
+
+  end subroutine reshape_into_ccpp_arr
+
+  ! Convert CAM-SIMA unit to MICM unit (kg kg-1  ->  mol m-3)
+  subroutine convert_to_mol_per_cubic_meter(dry_air_density, molar_mass_arr, constituents)
+    use ccpp_kinds,    only: kind_phys
+
+    real(kind_phys), intent(in)    :: dry_air_density(:,:) ! kg m-3
+    real(kind_phys), intent(in)    :: molar_mass_arr(:)    ! kg mol-1
+    real(kind_phys), intent(inout) :: constituents(:,:,:)  ! in: kg kg-1 | out: mol m-3
+
+    integer         :: num_columns, num_layers, num_constituents
+    integer         :: i_column, i_layer, i_elem
+    real(kind_phys) :: val
+
+    num_columns = size(constituents, dim=1)
+    num_layers = size(constituents, dim=2)
+    num_constituents = size(constituents, dim=3)
+
+    do i_elem = 1, num_constituents
+      do i_layer = 1, num_layers
+        do i_column = 1, num_columns
+          val = constituents(i_column, i_layer, i_elem) * dry_air_density(i_column, i_layer) &
+                / molar_mass_arr(i_elem)
+          constituents(i_column, i_layer, i_elem) = val
+        end do
+      end do
+    end do
+
+  end subroutine convert_to_mol_per_cubic_meter
+
+  ! Convert MICM unit to CAM-SIMA unit (mol m-3  ->  kg kg-1)
+  subroutine convert_to_mass_mixing_ratio(dry_air_density, molar_mass_arr, constituents)
+    use ccpp_kinds,    only: kind_phys
+
+    real(kind_phys), intent(in)    :: dry_air_density(:,:) ! kg m-3
+    real(kind_phys), intent(in)    :: molar_mass_arr(:)    ! kg mol-1
+    real(kind_phys), intent(inout) :: constituents(:,:,:)  ! in: mol m-3 | out: kg kg-1
+
+    integer         :: num_columns, num_layers, num_constituents
+    integer         :: i_column, i_layer, i_elem
+    real(kind_phys) :: val
+
+    num_columns = size(constituents, dim=1)
+    num_layers = size(constituents, dim=2)
+    num_constituents = size(constituents, dim=3)
+
+    do i_elem = 1, num_constituents
+      do i_layer = 1, num_layers
+        do i_column = 1, num_columns
+          val = constituents(i_column, i_layer, i_elem) / dry_air_density(i_column, i_layer) &
+                * molar_mass_arr(i_elem)
+          constituents(i_column, i_layer, i_elem) = val
+        end do
+      end do
+    end do
+
+  end subroutine convert_to_mass_mixing_ratio
+
+end module micm_util

musica/micm/musica_ccpp_micm.F90

@@ -1,12 +1,12 @@
 module musica_ccpp_micm
   use iso_c_binding
 
-  ! Note: "micm_core" is included in an external pre-built MICM library that the host
+  ! Note: "micm_t" is included in an external pre-built MICM library that the host
   ! model is responsible for linking to during compilation
-  use musica_micm, only: micm_t
-  use musica_ccpp_util, only: has_error_occurred
-  use ccpp_kinds, only: kind_phys
+  use musica_micm,          only: micm_t
+  use musica_ccpp_util,     only: has_error_occurred
   use musica_ccpp_namelist, only: filename_of_micm_configuration
+  use ccpp_kinds,           only: kind_phys
 
   implicit none
   private
@@ -18,23 +18,27 @@ module musica_ccpp_micm
 contains
 
   !> Register MICM constituents with the CCPP
-  subroutine micm_register(constituents, errmsg, errcode)
+  subroutine micm_register(constituents, solver_type, num_grid_cells, errmsg, errcode)
     use ccpp_constituent_prop_mod, only: ccpp_constituent_properties_t
-    use musica_util, only: error_t, mapping_t
-    type(ccpp_constituent_properties_t), allocatable, intent(out)   :: constituents(:)
-    character(len=512), intent(out) :: errmsg
-    integer, intent(out)            :: errcode
+    use musica_micm, only: Rosenbrock, RosenbrockStandardOrder
+    use musica_util, only: error_t
+
+    type(ccpp_constituent_properties_t), allocatable, intent(out) :: constituents(:)
+    integer(c_int),                                   intent(in)  :: solver_type
+    integer(c_int),                                   intent(in)  :: num_grid_cells
+    character(len=512),                               intent(out) :: errmsg
+    integer,                                          intent(out) :: errcode
 
-    type(error_t) :: error
-    type(mapping_t) :: mapping
+    ! local variables
+    type(error_t)        :: error
     real(kind=kind_phys) :: molar_mass
-    logical :: is_advected
-    integer :: i
+    logical              :: is_advected
+    integer              :: i
 
     errcode = 0
     errmsg = ''
 
-    micm => micm_t(filename_of_micm_configuration, error)
+    micm => micm_t(filename_of_micm_configuration, solver_type, num_grid_cells, error)
     if (has_error_occurred(error, errmsg, errcode)) return
 
     allocate(constituents(size(micm%species_ordering)), stat=errcode)
@@ -82,87 +86,41 @@ subroutine micm_init(errmsg, errcode)
   end subroutine micm_init
 
   !> Solve chemistry at the current time step
-  subroutine micm_run(time_step, temperature, pressure, dry_air_density, constituent_props, &
-                      constituents, errmsg, errcode)
-    use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
-    use musica_util, only: error_t
-
-    real(kind_phys),                   intent(in)    :: time_step            ! s
-    real(kind_phys),                   intent(in)    :: temperature(:,:)     ! K
-    real(kind_phys),                   intent(in)    :: pressure(:,:)        ! Pa
-    real(kind_phys),                   intent(in)    :: dry_air_density(:,:) ! kg m-3
-    type(ccpp_constituent_prop_ptr_t), intent(in)    :: constituent_props(:)
-    real(kind_phys),                   intent(inout) :: constituents(:,:,:)  ! kg kg-1
-    character(len=512),                intent(out)   :: errmsg
-    integer,                           intent(out)   :: errcode
+  subroutine micm_run(time_step, temperature, pressure, dry_air_density, constituents, &
+                      rate_params, errmsg, errcode)
+    use musica_micm, only: solver_stats_t
+    use musica_util, only: string_t, error_t
+
+    real(kind_phys),        intent(in)    :: time_step          ! s
+    real(c_double), target, intent(in)    :: temperature(:)     ! K
+    real(c_double), target, intent(in)    :: pressure(:)        ! Pa
+    real(c_double), target, intent(in)    :: dry_air_density(:) ! kg m-3
+    real(c_double), target, intent(inout) :: constituents(:)    ! mol m-3
+    real(c_double), target, intent(inout) :: rate_params(:)
+    character(len=512),     intent(out)   :: errmsg
+    integer,                intent(out)   :: errcode
 
     ! local variables
-    real(c_double)                                         :: c_time_step
-    real(c_double), dimension(size(temperature, dim=1), &
-                              size(temperature, dim=2))    :: c_temperature
-    real(c_double), dimension(size(pressure, dim=1), &
-                              size(pressure, dim=2))       :: c_pressure
-    real(c_double), dimension(size(constituents, dim=1), &
-                              size(constituents, dim=2), &
-                              size(constituents, dim=3))   :: c_constituents
-    real(c_double), dimension(size(constituents, dim=1), &
-                              size(constituents, dim=2), &
-                              0)                           :: c_rate_params
-
-    real(kind_phys), dimension(size(constituents, dim=3))  :: molar_mass_arr ! kg mol-1
-    type(error_t) :: error
-
-    integer :: num_columns, num_layers, num_constituents
-    integer :: i_column, i_layer, i_elem
-
-    num_columns = size(constituents, dim=1)
-    num_layers = size(constituents, dim=2)
-    num_constituents = size(constituents, dim=3)
+    type(string_t)       :: solver_state
+    type(solver_stats_t) :: solver_stats
+    type(error_t)        :: error
+    real(c_double)       :: c_time_step
+    integer              :: i_elem
 
     errcode = 0
     errmsg = ''
-
-    ! Get the molar_mass that is set in the call to instantiate()
-    do i_elem = 1, num_constituents
-      call constituent_props(i_elem)%molar_mass(molar_mass_arr(i_elem), errcode, errmsg)
-
-      if (errcode /= 0) then
-        errmsg = "[error] [micm] Unable to get molar mass."
-        return
-      end if
-    end do
-
-    ! TODO(jiwon) Check molar mass is non zero as it becomes a denominator for unit converison
-    ! this code needs to go when ccpp framework does the check
-    do i_elem = 1, num_constituents
-      if (molar_mass_arr(i_elem) == 0) then
-        errcode = 1
-        errmsg = "[error] [micm] Molar mass must be a non zero value."
-        return
-      end if
-    end do
-
-    ! Convert CAM-SIMA unit to MICM unit (kg kg-1  ->  mol m-3)
-    call convert_to_mol_per_cubic_meter(dry_air_density, molar_mass_arr, constituents)
-
-    c_time_step = real(time_step, c_double)
-    c_temperature = real(temperature, c_double)
-    c_pressure = real(pressure, c_double)
-    c_constituents = real(constituents, c_double)
-
-    do i_column = 1, num_columns
-      do i_layer = 1, num_layers
-        call micm%solve(c_temperature(i_column, i_layer), c_pressure(i_column, i_layer), &
-                        c_time_step, num_constituents, c_constituents(i_column, i_layer, :), &
-                        0, c_rate_params(i_column, i_layer, :), error)
-        if (has_error_occurred(error, errmsg, errcode)) return
-      end do
-    end do
-
-    constituents = real(c_constituents, kind_phys)
-
-    ! Convert MICM unit back to CAM-SIMA unit (mol m-3  ->  kg kg-1)
-    call convert_to_mass_mixing_ratio(dry_air_density, molar_mass_arr, constituents)
+    c_time_step = real(time_step, c_double) 
+
+    call micm%solve(c_time_step,     &
+                    temperature,     &
+                    pressure,        &
+                    dry_air_density, &
+                    constituents,    &
+                    rate_params,     &
+                    solver_state,    &
+                    solver_stats,    &
+                    error)
+    if (has_error_occurred(error, errmsg, errcode)) return
 
   end subroutine micm_run
 
@@ -176,58 +134,4 @@ subroutine micm_final(errmsg, errcode)
 
   end subroutine micm_final
 
-  ! Convert CAM-SIMA unit to MICM unit (kg kg-1  ->  mol m-3)
-  subroutine convert_to_mol_per_cubic_meter(dry_air_density, molar_mass_arr, constituents)
-    real(kind_phys), intent(in)    :: dry_air_density(:,:) ! kg m-3
-    real(kind_phys), intent(in)    :: molar_mass_arr(:)    ! kg mol-1
-    real(kind_phys), intent(inout) :: constituents(:,:,:)  ! in: kg kg-1 | out: mol m-3
-
-    integer :: num_columns, num_layers, num_constituents
-    integer :: i_column, i_layer, i_elem
-
-    real(kind_phys) :: val
-
-    num_columns = size(constituents, dim=1)
-    num_layers = size(constituents, dim=2)
-    num_constituents = size(constituents, dim=3)
-
-    do i_column = 1, num_columns
-      do i_layer = 1, num_layers
-        do i_elem = 1, num_constituents
-          val = constituents(i_column, i_layer, i_elem) * dry_air_density(i_column, i_layer) &
-                / molar_mass_arr(i_elem)
-          constituents(i_column, i_layer, i_elem) = val
-        end do
-      end do
-    end do
-
-  end subroutine convert_to_mol_per_cubic_meter
-
-  ! Convert MICM unit to CAM-SIMA unit (mol m-3  ->  kg kg-1)
-  subroutine convert_to_mass_mixing_ratio(dry_air_density, molar_mass_arr, constituents)
-    real(kind_phys), intent(in)    :: dry_air_density(:,:) ! kg m-3
-    real(kind_phys), intent(in)    :: molar_mass_arr(:)    ! kg mol-1
-    real(kind_phys), intent(inout) :: constituents(:,:,:)  ! in: mol m-3 | out: kg kg-1
-
-    integer :: num_columns, num_layers, num_constituents
-    integer :: i_column, i_layer, i_elem
-
-    real(kind_phys) :: val
-
-    num_columns = size(constituents, dim=1)
-    num_layers = size(constituents, dim=2)
-    num_constituents = size(constituents, dim=3)
-
-    do i_column = 1, num_columns
-      do i_layer = 1, num_layers
-        do i_elem = 1, num_constituents
-          val = constituents(i_column, i_layer, i_elem) / dry_air_density(i_column, i_layer) &
-                * molar_mass_arr(i_elem)
-          constituents(i_column, i_layer, i_elem) = val
-        end do
-      end do
-    end do
-
-  end subroutine convert_to_mass_mixing_ratio
-
 end module musica_ccpp_micm