add conversion from geopotential height to regular height

schemes/musica/musica_ccpp.F90

@@ -43,7 +43,10 @@ end subroutine musica_ccpp_init
   !> \section arg_table_musica_ccpp_run Argument Table
   !! \htmlinclude musica_ccpp_run.html
   subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, constituent_props, &
-                      constituents, height_midpoints, height_interfaces, errmsg, errcode)
+                      constituents, geopotential_height_wrt_surface_at_midpoint, &
+                      geopotential_height_wrt_surface_at_interface, &
+                      surface_geopotential, reciprocal_of_gravitational_acceleration, &
+                      errmsg, errcode)
     use musica_ccpp_micm_util,     only: reshape_into_micm_arr, reshape_into_ccpp_arr
     use musica_ccpp_micm_util,     only: convert_to_mol_per_cubic_meter, convert_to_mass_mixing_ratio
     use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
@@ -55,8 +58,10 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co
     real(kind_phys), target,           intent(in)    :: dry_air_density(:,:)   ! kg m-3
     type(ccpp_constituent_prop_ptr_t), intent(in)    :: constituent_props(:)
     real(kind_phys), target,           intent(inout) :: constituents(:,:,:)    ! kg kg-1
-    real(kind_phys), target,           intent(in)    :: height_midpoints(:,:)  ! km
-    real(kind_phys), target,           intent(in)    :: height_interfaces(:,:) ! km
+    real(kind_phys), target,           intent(in)    :: geopotential_height_wrt_surface_at_midpoint(:,:)  ! m
+    real(kind_phys), target,           intent(in)    :: geopotential_height_wrt_surface_at_interface(:,:) ! m
+    real(kind_phys), target,           intent(in)    :: surface_geopotential(:)    ! m2 s-2
+    real(kind_phys), target,           intent(in)    :: reciprocal_of_gravitational_acceleration ! s2 m-1
     character(len=512),                intent(out)   :: errmsg
     integer,                           intent(out)   :: errcode
 
@@ -78,7 +83,9 @@ subroutine musica_ccpp_run(time_step, temperature, pressure, dry_air_density, co
                                     * 3) :: photolysis_rate_constants ! s-1
     integer :: i_elem
 
-    call tuvx_run(temperature, dry_air_density, height_midpoints, height_interfaces, &
+    call tuvx_run(temperature, dry_air_density, geopotential_height_wrt_surface_at_midpoint, &
+                  geopotential_height_wrt_surface_at_interface, &
+                  surface_geopotential, reciprocal_of_gravitational_acceleration, &
                   photolysis_rate_constants, errmsg, errcode)
 
     ! Get the molar_mass that is set in the call to instantiate()

schemes/musica/musica_ccpp.meta

@@ -71,18 +71,30 @@
   type = real | kind = kind_phys
   dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_ccpp_constituents)
   intent = inout
-[ height_midpoints ]
-  standard_name = ccpp_height_midpoints
+[ geopotential_height_wrt_surface_at_midpoint ]
+  standard_name = geopotential_height_wrt_surface
   units = km
   type = real | kind = kind_phys
   dimensions = (horizontal_loop_extent,vertical_layer_dimension)
   intent = in
-[ height_interfaces ]
-  standard_name = ccpp_height_interfaces
+[ geopotential_height_wrt_surface_at_interface ]
+  standard_name = geopotential_height_wrt_surface_at_interface
   units = km
   type = real | kind = kind_phys
   dimensions = (horizontal_loop_extent,vertical_layer_dimension)
   intent = in
+[ surface_geopotential ]
+  standard_name = surface_geopotential
+  type = real | kind = kind_phys
+  units = m2 s-2
+  dimensions = (horizontal_loop_extent)
+  intent = in
+[ reciprocal_of_gravitational_acceleration ]
+  standard_name = reciprocal_of_gravitational_acceleration
+  units = s2 m-1
+  type = real | kind = kind_phys
+  dimensions = ()
+  intent = in
 [ errmsg ]
   standard_name = ccpp_error_message
   units = none

schemes/musica/tuvx/musica_ccpp_tuvx.F90

@@ -96,15 +96,20 @@ subroutine tuvx_init(vertical_layer_dimension, &
   end subroutine tuvx_init
 
   !> Calculates photolysis rate constants for the current model conditions
-  subroutine tuvx_run( temperature, dry_air_density, height_midpoints, &
-      height_interfaces, photolysis_rate_constants, errmsg, errcode )
+  subroutine tuvx_run( temperature, dry_air_density, &
+      geopotential_height_wrt_surface_at_midpoint, &
+      geopotential_height_wrt_surface_at_interface, &
+      surface_geopotential, reciprocal_of_gravitational_acceleration, &
+      photolysis_rate_constants, errmsg, errcode )
     use musica_util, only: error_t
-    use musica_ccpp_tuvx_height_grid, only: set_height_grid_values
+    use musica_ccpp_tuvx_height_grid, only: set_height_grid_values, calculate_heights
 
     real(kind_phys),    intent(in)  :: temperature(:,:)       ! K (column, layer)
     real(kind_phys),    intent(in)  :: dry_air_density(:,:)   ! molecule cm-3 (column, layer)
-    real(kind_phys),    intent(in)  :: height_midpoints(:,:)  ! km (column, layer)
-    real(kind_phys),    intent(in)  :: height_interfaces(:,:) ! km (column, interface)
+    real(kind_phys),    intent(in)  :: geopotential_height_wrt_surface_at_midpoint(:,:)  ! m (column, layer)
+    real(kind_phys),    intent(in)  :: geopotential_height_wrt_surface_at_interface(:,:) ! m (column, interface)
+    real(kind_phys),    intent(in)  :: surface_geopotential(:)    ! m2 s-2
+    real(kind_phys),    intent(in)  :: reciprocal_of_gravitational_acceleration ! s2 m-1
     ! temporarily set to Chapman mechanism and 1 dimension
     ! until mapping between MICM and TUV-x is implemented
     real(kind_phys),    intent(out) :: photolysis_rate_constants(:) ! s-1 (column, reaction)
@@ -113,14 +118,20 @@ subroutine tuvx_run( temperature, dry_air_density, height_midpoints, &
 
     ! local variables
     type(error_t) :: error
+    real(kind_phys), dimension(size(geopotential_height_wrt_surface_at_midpoint, dim = 2)) :: height_midpoints
+    real(kind_phys), dimension(size(geopotential_height_wrt_surface_at_interface, dim = 2)) :: height_interfaces
     integer :: i_col
 
     errcode = 0
     errmsg = ''
 
     do i_col = 1, size(temperature, dim=1)
-      call set_height_grid_values( height_grid, height_midpoints(i_col,:), &
-          height_interfaces(i_col,:), errmsg, errcode )
+      call calculate_heights( geopotential_height_wrt_surface_at_midpoint(i_col,:), &
+          geopotential_height_wrt_surface_at_interface(i_col,:), &
+          surface_geopotential(i_col), reciprocal_of_gravitational_acceleration, &
+          height_midpoints, height_interfaces )
+      call set_height_grid_values( height_grid, height_midpoints, &
+          height_interfaces, errmsg, errcode )
     end do
     if (errcode /= 0) return
 

schemes/musica/tuvx/musica_ccpp_tuvx_height_grid.F90

@@ -3,7 +3,7 @@ module musica_ccpp_tuvx_height_grid
   implicit none
 
   private
-  public :: create_height_grid, set_height_grid_values
+  public :: create_height_grid, set_height_grid_values, calculate_heights
 
   ! Conversions between the CAM-SIMA height grid and the TUVX height grid
   !
@@ -135,4 +135,39 @@ subroutine set_height_grid_values( height_grid, host_midpoints, &
 
   end subroutine set_height_grid_values
 
+  !> Calculates the heights needed for the TUV-x grid based on available data
+  !!
+  !! Uses the reciprocal of gravitational acceleration, the surface geopotential,
+  !! and the geopotential height wrt the surface and midpoints/interfaces to calculate
+  !! the midpoint/interface height above sea level in km needed for the TUV-x grid.
+  !!
+  !! The equation used is taked from CAMChem
+  !! (see https://github.com/ESCOMP/CAM/blob/f0e489e9708ce7b91635f6d4997fbf1e390b0dbb/src/chemistry/mozart/mo_gas_phase_chemdr.F90#L514-L526)
+  subroutine calculate_heights( geopotential_height_wrt_surface_at_midpoint, &
+      geopotential_height_wrt_surface_at_interface, &
+      surface_geopotential, reciprocal_of_gravitational_acceleration, &
+      height_midpoints, height_interfaces )
+
+    use ccpp_kinds,       only: kind_phys
+    use musica_ccpp_util, only: has_error_occurred
+    use musica_util,      only: error_t
+
+    real(kind_phys),  intent(in)  :: geopotential_height_wrt_surface_at_midpoint(:) ! m
+    real(kind_phys),  intent(in)  :: geopotential_height_wrt_surface_at_interface(:) ! m
+    real(kind_phys),  intent(in)  :: surface_geopotential ! m2 s-2
+    real(kind_phys),  intent(in)  :: reciprocal_of_gravitational_acceleration ! s2 m-1
+    real(kind_phys),  intent(out) :: height_midpoints(:) ! Pa
+    real(kind_phys),  intent(out) :: height_interfaces(:) ! Pa
+    real(kind_phys) :: surface_height ! km
+
+    surface_height = surface_geopotential * reciprocal_of_gravitational_acceleration
+    height_midpoints(:) = 0.001_kind_phys * &
+                          ( geopotential_height_wrt_surface_at_midpoint(:) &
+                            + surface_height )
+    height_interfaces(:) = 0.001_kind_phys * &
+                           ( geopotential_height_wrt_surface_at_interface(:) &
+                             + surface_height )
+
+  end subroutine calculate_heights
+
 end module musica_ccpp_tuvx_height_grid

test/musica/test_musica_api.F90

@@ -19,8 +19,10 @@ subroutine test_musica_ccpp_api()
   integer                                                                :: errcode
   character(len=512)                                                     :: errmsg
   real(kind_phys)                                                        :: time_step         ! s
-  real(kind_phys),         dimension(NUM_COLUMNS,NUM_LAYERS)             :: height_midpoints  ! km
-  real(kind_phys),         dimension(NUM_COLUMNS,NUM_LAYERS+1)           :: height_interfaces ! km
+  real(kind_phys),         dimension(NUM_COLUMNS,NUM_LAYERS)             :: geopotential_height_wrt_surface_at_midpoint  ! m
+  real(kind_phys),         dimension(NUM_COLUMNS,NUM_LAYERS+1)           :: geopotential_height_wrt_surface_at_interface ! m
+  real(kind_phys),         dimension(NUM_COLUMNS)                        :: surface_geopotential ! m2 s-2
+  real(kind_phys)                                                        :: reciprocal_of_gravitational_acceleration ! s2 m-1
   real(kind_phys), target, dimension(NUM_COLUMNS,NUM_LAYERS)             :: temperature       ! K
   real(kind_phys), target, dimension(NUM_COLUMNS,NUM_LAYERS)             :: pressure          ! Pa
   real(kind_phys), target, dimension(NUM_COLUMNS,NUM_LAYERS)             :: dry_air_density   ! kg m-3
@@ -39,10 +41,12 @@ subroutine test_musica_ccpp_api()
   solver_type = Rosenbrock
   num_grid_cells = NUM_COLUMNS * NUM_LAYERS
   time_step = 60._kind_phys
-  height_midpoints(1,:) = (/ 3.0_kind_phys, 1.0_kind_phys /)
-  height_midpoints(2,:) = (/ 4.0_kind_phys, 1.5_kind_phys /)
-  height_interfaces(1,:) = (/ 4.0_kind_phys, 2.0_kind_phys, 1.0_kind_phys /)
-  height_interfaces(2,:) = (/ 5.0_kind_phys, 2.5_kind_phys, 0.5_kind_phys /)
+  geopotential_height_wrt_surface_at_midpoint(1,:) = (/ 2000.0_kind_phys, 500.0_kind_phys /)
+  geopotential_height_wrt_surface_at_midpoint(2,:) = (/ 2000.0_kind_phys, -500.0_kind_phys /)
+  geopotential_height_wrt_surface_at_interface(1,:) = (/ 3000.0_kind_phys, 1000.0_kind_phys, 0.0_kind_phys /)
+  geopotential_height_wrt_surface_at_interface(2,:) = (/ 3000.0_kind_phys, 500.0_kind_phys, -1500.0_kind_phys /)
+  surface_geopotential = (/ 100.0_kind_phys, 200.0_kind_phys /)
+  reciprocal_of_gravitational_acceleration = 10.0_kind_phys
   temperature(:,1) = (/ 100._kind_phys, 200._kind_phys /)
   temperature(:,2) = (/ 300._kind_phys, 400._kind_phys /)
   pressure(:,1) = (/ 6000.04_kind_phys, 7000.04_kind_phys /)
@@ -102,7 +106,10 @@ subroutine test_musica_ccpp_api()
   write(*,fmt="(4(3x,e13.6))") constituents
 
   call musica_ccpp_run(time_step, temperature, pressure, dry_air_density, constituent_props_ptr, &
-                       constituents, height_midpoints, height_interfaces, errmsg, errcode)
+                       constituents, geopotential_height_wrt_surface_at_midpoint, &
+                       geopotential_height_wrt_surface_at_interface, &
+                       surface_geopotential, reciprocal_of_gravitational_acceleration, &
+                       errmsg, errcode)
   if (errcode /= 0) then
     write(*,*) trim(errmsg)
     stop 3

test/musica/tuvx/test_tuvx_height_grid.F90

@@ -10,6 +10,7 @@ program test_tuvx_height_grid
 #define ASSERT_NEAR( a, b, abs_error ) if( (abs(a - b) >= abs_error) .and. (abs(a - b) /= 0.0) ) then; write(*,*) "Assertion failed[", __FILE__, ":", __LINE__, "]: a, b"; stop 1; endif
 
   call test_create_height_grid()
+  call test_calculate_height_grid_values()
 
 contains
 
@@ -73,4 +74,34 @@ subroutine test_create_height_grid()
 
   end subroutine test_create_height_grid
 
+  subroutine test_calculate_height_grid_values()
+
+    use ccpp_kinds, only: kind_phys
+
+    integer, parameter                       :: NUM_LAYERS = 2
+    real(kind_phys), dimension(NUM_LAYERS)   :: geopotential_height_wrt_surface_at_midpoint ! m
+    real(kind_phys), dimension(NUM_LAYERS+1) :: geopotential_height_wrt_surface_at_interface ! m
+    real(kind_phys)                          :: surface_geopotential ! m2 s-2
+    real(kind_phys)                          :: reciprocal_of_gravitational_acceleration ! s2 m-1
+    real(kind_phys), dimension(NUM_LAYERS)   :: height_midpoints
+    real(kind_phys), dimension(NUM_LAYERS+1) :: height_interfaces
+
+    geopotential_height_wrt_surface_at_midpoint(:) = (/ 2000.0_kind_phys, 500.0_kind_phys /)
+    geopotential_height_wrt_surface_at_interface(:) = (/ 3000.0_kind_phys, 1000.0_kind_phys, 0.0_kind_phys /)
+    surface_geopotential = 100.0_kind_phys
+    reciprocal_of_gravitational_acceleration = 10.0_kind_phys
+
+    call calculate_heights(geopotential_height_wrt_surface_at_midpoint, &
+                           geopotential_height_wrt_surface_at_interface, &
+                           surface_geopotential, reciprocal_of_gravitational_acceleration, &
+                           height_midpoints, height_interfaces)
+
+    ASSERT_NEAR(height_midpoints(1), 3.0, 1e-5)
+    ASSERT_NEAR(height_midpoints(2), 1.5, 1e-5)
+    ASSERT_NEAR(height_interfaces(1), 4.0, 1e-5)
+    ASSERT_NEAR(height_interfaces(2), 2.0, 1e-5)
+    ASSERT_NEAR(height_interfaces(3), 1.0, 1e-5)
+
+  end subroutine test_calculate_height_grid_values
+
 end program test_tuvx_height_grid