mpasjedi_linvarcha_c2a_mod.F90

Go to the documentation of this file.

! (C) Copyright 2018 UCAR
! 
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 
 
module mpasjedi_linvarcha_c2a_mod
 
use fckit_configuration_module, only: fckit_configuration
use fckit_log_module, only : fckit_log
use iso_c_binding
use kinds
 
!oops
use oops_variables_mod, only: oops_variables
 
!mpas-jedi
use mpas_constants_mod
use mpas_fields_mod, only: mpas_fields
use mpas_geom_mod, only: mpas_geom
 
!MPAS-Model
use mpas_derived_types
use mpas_field_routines
use mpas_pool_routines
use mpas_dmpar
 
implicit none
 
private
public :: mpasjedi_linvarcha_c2a
 
!> Fortran derived type to hold configuration data for the B mat variable change
type :: mpasjedi_linvarcha_c2a
  type (mpas_pool_type), pointer :: trajectories => null()
  contains
    procedure, public :: create
    procedure, public :: delete
    procedure, public :: multiply
    procedure, public :: multiplyadjoint
    procedure, public :: multiplyinverse
    procedure, public :: multiplyinverseadjoint
end type mpasjedi_linvarcha_c2a
 
character(len=1024) :: message
 
! ------------------------------------------------------------------------------
contains
! ------------------------------------------------------------------------------
 
subroutine create(self, geom, bg, fg, f_conf, vars)
 
   implicit none
   class(mpasjedi_linvarcha_c2a),intent(inout) :: self
   type(mpas_fields), target,    intent(in)    :: bg
   type(mpas_fields), target,    intent(in)    :: fg
   type(mpas_geom),              intent(in)    :: geom
   type(fckit_configuration),    intent(in)    :: f_conf
   type(oops_variables),         intent(in)    :: vars
 
   type (field2DReal), pointer :: fld2d_t, fld2d_p, fld2d_qs
 
   integer :: ngrid
 
   if ( vars % has ('relhum') ) then
      !-- set trajectories for linear variable change
      ngrid = geom % nCells ! local + halo
 
      call mpas_pool_create_pool(self % trajectories)
 
      if ( .not. fg % has ('temperature') .or. .not. fg % has ('pressure') ) &
         call abor1_ftn("LinVarChaC2A::LinVarChaC2A, mpasjedi_linvarcha_c2a_create: Trajectory failed")
 
      call fg%get('temperature', fld2d_t)
      call fg%get('pressure'   , fld2d_p)
 
      call mpas_duplicate_field(fld2d_t, fld2d_qs) ! for saturation specific humidity
 
      call da_tp_to_qs( fld2d_t % array(:,1:ngrid), fld2d_p % array(:,1:ngrid), fld2d_qs% array(:,1:ngrid))
 
      call mpas_pool_add_field(self % trajectories, 'spechum_sat', fld2d_qs)
   end if
 
end subroutine create
 
! ------------------------------------------------------------------------------
 
subroutine delete(self)
 
   implicit none
   class(mpasjedi_linvarcha_c2a),intent(inout) :: self
 
   if (associated(self % trajectories)) then
      call mpas_pool_destroy_pool(self % trajectories)
   end if
 
end subroutine delete
 
! ------------------------------------------------------------------------------
 
subroutine multiply(self,geom,xctl,xana)
 
   implicit none
   class(mpasjedi_linvarcha_c2a),intent(inout) :: self
   type(mpas_geom),              intent(inout) :: geom
   type(mpas_fields),            intent(inout) :: xctl
   type(mpas_fields),            intent(inout) :: xana
 
   type (mpas_pool_iterator_type) :: poolItr
   type (field2DReal), pointer :: fld2d_sf, fld2d_vp, fld2d_uRz, fld2d_uRm, fld2d_traj_qs
   type (field2DReal), pointer :: fld2d_v_src, fld2d_sf_v, fld2d_e_src, fld2d_u
 
   real(kind=kind_real), dimension(:), pointer :: ptrr1_ctl, ptrr1_ana
   real(kind=kind_real), dimension(:,:), pointer :: ptrr2_ctl, ptrr2_ana
   real(kind=kind_real), dimension(:,:), pointer :: ptrr2_sh, ptrr2_rh
 
   integer :: k, ngrid, nCells, nVertices, nEdges
 
   write(message,*) "DEBUG: mpasjedi_linvarcha_c2a_multiply: xana % fldnames(:) =",xana % fldnames(:)
   call fckit_log%debug(message)
   write(message,*) "DEBUG: mpasjedi_linvarcha_c2a_multiply: xctl % fldnames(:) =",xctl % fldnames(:)
   call fckit_log%debug(message)
 
   ngrid = geom % nCellsSolve ! local 
 
   ! no variable change
   call mpas_pool_begin_iteration(xctl%subFields)
   do while ( mpas_pool_get_next_member(xctl%subFields, poolitr) )
      ! NOTE: only 1d or 2d real fields are handled.
      if (poolitr % memberType == mpas_pool_field .and. poolitr % dataType == mpas_pool_real) then 
         if( xana%has(poolitr % memberName) ) then ! check if the variable names are the same.
            if (poolitr % nDims == 1) then
               call xctl%get(poolitr % memberName, ptrr1_ctl)
               call xana%get(poolitr % memberName, ptrr1_ana)
               ptrr1_ana(1:ngrid)=ptrr1_ctl(1:ngrid)
            else if (poolitr % nDims == 2) then
               call xctl%get(poolitr % memberName, ptrr2_ctl)
               call xana%get(poolitr % memberName, ptrr2_ana)
               ptrr2_ana(:,1:ngrid)=ptrr2_ctl(:,1:ngrid)
            else
               write(message,*) '--> mpasjedi_linvarcha_c2a::multiply: poolItr % nDims == ',poolitr % nDims,' not handled'
               call abor1_ftn(message)
            end if
         end if
      end if
   end do !- end of pool iteration
 
   if( xctl%has('relhum') .and. xana%has('spechum') ) then
      call xctl%get(     'relhum', ptrr2_rh)
      call xana%get(    'spechum', ptrr2_sh)
      call mpas_pool_get_field(self%trajectories, 'spechum_sat', fld2d_traj_qs)
      call pseudorh_to_spechum( ptrr2_rh(:,1:ngrid), ptrr2_sh(:,1:ngrid), fld2d_traj_qs%array(:,1:ngrid) )
   end if
 
   if( xctl%has('stream_function') .and. xctl%has('velocity_potential') .and. &
       xana%has('uReconstructZonal') .and. xana%has('uReconstructMeridional') ) then
      call xctl%get(       'stream_function', fld2d_sf)
      call xctl%get(    'velocity_potential', fld2d_vp)
      call xana%get(     'uReconstructZonal', fld2d_urz)
      call xana%get('uReconstructMeridional', fld2d_urm)
 
      call mpas_dmpar_exch_halo_field(fld2d_sf)
      call mpas_dmpar_exch_halo_field(fld2d_vp)
 
      ncells = geom % nCells    ! local + halo
      nvertices = geom % nVertices ! local + halo
      nedges = geom % nEdges    ! local + halo
 
      ! duplicate two temporary working spaces
      call mpas_pool_get_field( geom % domain % blocklist % allFields, 'vorticity', fld2d_v_src)
      call mpas_duplicate_field(fld2d_v_src, fld2d_sf_v)
      fld2d_sf_v % fieldName = 'stream_function at vertices'
      call mpas_pool_get_field( geom % domain % blocklist % allFields, 'u', fld2d_e_src)
      call mpas_duplicate_field(fld2d_e_src, fld2d_u)
      fld2d_u % fieldName = 'Horizontal normal velocity at edges'
 
      call psichi_to_uv_edge_step1( geom, fld2d_sf%array(:,1:ncells), fld2d_sf_v%array(:,1:nvertices) )
      call mpas_dmpar_exch_halo_field(fld2d_sf_v)
      call psichi_to_uv_edge_step2( geom, fld2d_sf_v%array(:,1:nvertices), fld2d_vp%array(:,1:ncells), &
                                   fld2d_u%array(:,1:nedges) )
      call mpas_dmpar_exch_halo_field(fld2d_u)
      call psichi_to_uv_edge_step3( geom, fld2d_u%array(:,1:nedges), &
                                   fld2d_urz%array(:,1:ncells), fld2d_urm%array(:,1:ncells) )
 
      call mpas_deallocate_field(fld2d_sf_v)
      call mpas_deallocate_field(fld2d_u)
   end if
 
end subroutine multiply
 
! ------------------------------------------------------------------------------
 
subroutine multiplyadjoint(self,geom,xana,xctl)
 
   implicit none
   class(mpasjedi_linvarcha_c2a),intent(inout) :: self
   type(mpas_geom),              intent(inout) :: geom
   type(mpas_fields),            intent(inout) :: xana
   type(mpas_fields),            intent(inout) :: xctl
 
   type (mpas_pool_iterator_type) :: poolItr
   type (field2DReal), pointer :: fld2d_sf, fld2d_vp, fld2d_uRz, fld2d_uRm, fld2d_traj_qs
   type (field2DReal), pointer :: fld2d_v_src, fld2d_sf_v, fld2d_e_src, fld2d_u
 
   real(kind=kind_real), dimension(:), pointer :: ptrr1_ctl, ptrr1_ana
   real(kind=kind_real), dimension(:,:), pointer :: ptrr2_ctl, ptrr2_ana
   real(kind=kind_real), dimension(:,:), pointer :: ptrr2_sh, ptrr2_rh
   integer :: k, ngrid, nCells, nVertices, nEdges
 
   write(message,*) "DEBUG: mpasjedi_linvarcha_c2a_multiplyadjoint: xana % fldnames(:) =",xana % fldnames(:)
   call fckit_log%debug(message)
   write(message,*) "DEBUG: mpasjedi_linvarcha_c2a_multiplyadjoint: xctl % fldnames(:) =",xctl % fldnames(:)
   call fckit_log%debug(message)
 
   ngrid = geom % nCellsSolve
 
   ! no variable change
   call mpas_pool_begin_iteration(xctl%subFields)
   do while ( mpas_pool_get_next_member(xctl%subFields, poolitr) )
      ! NOTE: only 1d or 2d real fields are handled.
      if (poolitr % memberType == mpas_pool_field .and. poolitr % dataType == mpas_pool_real) then
         if( xana%has(poolitr % memberName) ) then ! check if the variable names are the same.
            if (poolitr % nDims == 1) then
               call xctl%get(poolitr % memberName, ptrr1_ctl)
               call xana%get(poolitr % memberName, ptrr1_ana)
               ptrr1_ctl(1:ngrid)=ptrr1_ctl(1:ngrid)+ptrr1_ana(1:ngrid)
               ptrr1_ana(1:ngrid)=mpas_jedi_zero_kr
            else if (poolitr % nDims == 2) then
               call xctl%get(poolitr % memberName, ptrr2_ctl)
               call xana%get(poolitr % memberName, ptrr2_ana)
               ptrr2_ctl(:,1:ngrid)=ptrr2_ctl(:,1:ngrid)+ptrr2_ana(:,1:ngrid)
               ptrr2_ana(:,1:ngrid)=mpas_jedi_zero_kr
            else
               write(message,*) '--> mpasjedi_linvarcha_c2a::multiplyadjoint: poolItr % nDims == ',poolitr % nDims,' not handled'
               call abor1_ftn(message)
            end if
         end if
      end if
   end do !- end of pool iteration
 
   if( xctl%has('relhum') .and. xana%has('spechum') ) then
      call xctl%get(     'relhum', ptrr2_rh)
      call xana%get(    'spechum', ptrr2_sh)
      call mpas_pool_get_field(self%trajectories, 'spechum_sat', fld2d_traj_qs)
      call pseudorh_to_spechumad( ptrr2_rh(:,1:ngrid), ptrr2_sh(:,1:ngrid), fld2d_traj_qs%array(:,1:ngrid) )
   end if
 
   if( xctl%has('stream_function') .and. xctl%has('velocity_potential') .and. &
       xana%has('uReconstructZonal') .and. xana%has('uReconstructMeridional') ) then
      call xctl%get(       'stream_function', fld2d_sf)
      call xctl%get(    'velocity_potential', fld2d_vp)
      call xana%get(     'uReconstructZonal', fld2d_urz)
      call xana%get('uReconstructMeridional', fld2d_urm)
 
      ncells = geom % nCells    ! local + halo
      nvertices = geom % nVertices ! local + halo
      nedges = geom % nEdges    ! local + halo
 
      ! duplicate two temporary working spaces
      call mpas_pool_get_field( geom % domain % blocklist % allFields, 'vorticity', fld2d_v_src)
      call mpas_duplicate_field(fld2d_v_src, fld2d_sf_v)
      fld2d_sf_v % fieldName = 'stream_function at vertices'
      call mpas_pool_get_field( geom % domain % blocklist % allFields, 'u', fld2d_e_src)
      call mpas_duplicate_field(fld2d_e_src, fld2d_u)
      fld2d_u % fieldName = 'Horizontal normal velocity at edges'
 
      call psichi_to_uv_edge_step3ad( geom, fld2d_u%array(:,1:nedges), &
                                     fld2d_urz%array(:,1:ncells), fld2d_urm%array(:,1:ncells) )
      call mpas_dmpar_exch_halo_adj_field(fld2d_u)
      call psichi_to_uv_edge_step2ad( geom, fld2d_sf_v%array(:,1:nvertices), fld2d_vp%array(:,1:ncells), &
                                     fld2d_u%array(:,1:nedges) )
      call mpas_dmpar_exch_halo_adj_field(fld2d_sf_v)
      call psichi_to_uv_edge_step1ad( geom, fld2d_sf%array(:,1:ncells), fld2d_sf_v%array(:,1:nvertices) )
 
      call mpas_deallocate_field(fld2d_sf_v)
      call mpas_deallocate_field(fld2d_u)
 
      fld2d_urz%array(:,1:geom%nCellsSolve)=mpas_jedi_zero_kr
      fld2d_urm%array(:,1:geom%nCellsSolve)=mpas_jedi_zero_kr
 
      call mpas_dmpar_exch_halo_adj_field(fld2d_sf)
      call mpas_dmpar_exch_halo_adj_field(fld2d_vp)
   end if
 
end subroutine multiplyadjoint
 
! ------------------------------------------------------------------------------
 
subroutine multiplyinverse(self,geom,xana,xctl)
 
   implicit none
   class(mpasjedi_linvarcha_c2a),intent(inout) :: self
   type(mpas_geom),              intent(inout) :: geom
   type(mpas_fields),            intent(inout) :: xana
   type(mpas_fields),            intent(inout) :: xctl
 
   type (mpas_pool_iterator_type) :: poolItr
   type (field2DReal), pointer :: fld2d_traj_qs
   real(kind=kind_real), dimension(:), pointer :: ptrr1_ctl, ptrr1_ana
   real(kind=kind_real), dimension(:,:), pointer :: ptrr2_ctl, ptrr2_ana
   real(kind=kind_real), dimension(:,:), pointer :: ptrr2_sh, ptrr2_rh
   integer :: ngrid
 
   write(message,*) "DEBUG: mpasjedi_linvarcha_c2a_multiplyinverse: xana % fldnames(:) =",xana % fldnames(:)
   call fckit_log%debug(message)
   write(message,*) "DEBUG: mpasjedi_linvarcha_c2a_multiplyinverse: xctl % fldnames(:) =",xctl % fldnames(:)
   call fckit_log%debug(message)
 
   ngrid = geom % nCellsSolve ! local
 
   ! no variable change
   call mpas_pool_begin_iteration(xctl%subFields)
   do while ( mpas_pool_get_next_member(xctl%subFields, poolitr) )
      ! only 1d or 2d real fields are handled.
      if (poolitr % memberType == mpas_pool_field .and. poolitr % dataType == mpas_pool_real) then
         if( xana%has(poolitr % memberName) ) then ! check if the variable names are the same.
            if (poolitr % nDims == 1) then
               call xctl%get(poolitr % memberName, ptrr1_ctl)
               call xana%get(poolitr % memberName, ptrr1_ana)
               ptrr1_ctl(1:ngrid)=ptrr1_ana(1:ngrid)
            else if (poolitr % nDims == 2) then
               call xctl%get(poolitr % memberName, ptrr2_ctl)
               call xana%get(poolitr % memberName, ptrr2_ana)
               ptrr2_ctl(:,1:ngrid)=ptrr2_ana(:,1:ngrid)
            else
               write(message,*) '--> mpasjedi_linvarcha_c2a::multiplyinverse: poolItr % nDims == ',poolitr % nDims,' not handled'
               call abor1_ftn(message)
            end if
         end if
      end if
   end do !- end of pool iteration
 
   if( xctl%has('relhum') .and. xana%has('spechum') ) then
      call xana%get(    'spechum', ptrr2_sh)
      call xctl%get(     'relhum', ptrr2_rh)
      call mpas_pool_get_field(self%trajectories, 'spechum_sat', fld2d_traj_qs)
      call pseudorh_to_spechum_inverse( ptrr2_rh(:,1:ngrid), ptrr2_sh(:,1:ngrid), fld2d_traj_qs%array(:,1:ngrid) )
   end if
 
   !-- dummy inverse operator
   if( xctl%has('stream_function') .and. xctl%has('velocity_potential') .and. &
       xana%has('uReconstructZonal') .and. xana%has('uReconstructMeridional') ) then
      call xana%get('uReconstructZonal', ptrr2_ana)
      call xctl%get(  'stream_function', ptrr2_ctl)
      ptrr2_ctl(:,1:ngrid)=ptrr2_ana(:,1:ngrid)
 
      call xana%get('uReconstructMeridional', ptrr2_ana)
      call xctl%get(    'velocity_potential', ptrr2_ctl)
      ptrr2_ctl(:,1:ngrid)=ptrr2_ana(:,1:ngrid)
   end if
 
end subroutine multiplyinverse
 
! ------------------------------------------------------------------------------
 
subroutine multiplyinverseadjoint(self,geom,xctl,xana)
 
   implicit none
   class(mpasjedi_linvarcha_c2a),intent(inout) :: self
   type(mpas_geom),              intent(inout) :: geom
   type(mpas_fields),            intent(inout) :: xctl
   type(mpas_fields),            intent(inout) :: xana
 
   type (mpas_pool_iterator_type) :: poolItr
   type (field2DReal), pointer :: fld2d_traj_qs
   real(kind=kind_real), dimension(:), pointer :: ptrr1_ctl, ptrr1_ana
   real(kind=kind_real), dimension(:,:), pointer :: ptrr2_ctl, ptrr2_ana
   real(kind=kind_real), dimension(:,:), pointer :: ptrr2_sh, ptrr2_rh, ptrr2_traj_qs
   integer :: ngrid
 
   write(message,*) "DEBUG: mpasjedi_linvarcha_c2a_multiplyinverseadjoint: xana % fldnames(:) =",xana % fldnames(:)
   call fckit_log%debug(message)
   write(message,*) "DEBUG: mpasjedi_linvarcha_c2a_multiplyinverseadjoint: xctl % fldnames(:) =",xctl % fldnames(:)
   call fckit_log%debug(message)
 
   ngrid = geom % nCellsSolve ! local
 
   ! no variable change
   call mpas_pool_begin_iteration(xctl%subFields)
   do while ( mpas_pool_get_next_member(xctl%subFields, poolitr) )
      ! only 1d or 2d real fields are handled.
      if (poolitr % memberType == mpas_pool_field .and. poolitr % dataType == mpas_pool_real) then
         if( xana%has(poolitr % memberName) ) then ! check if the variable names are the same.
            if (poolitr % nDims == 1) then
               call xctl%get(poolitr % memberName, ptrr1_ctl)
               call xana%get(poolitr % memberName, ptrr1_ana)
               ptrr1_ana(1:ngrid)=ptrr1_ana(1:ngrid)+ptrr1_ctl(1:ngrid)
               ptrr1_ctl(1:ngrid)=mpas_jedi_zero_kr
            else if (poolitr % nDims == 2) then
               call xctl%get(poolitr % memberName, ptrr2_ctl)
               call xana%get(poolitr % memberName, ptrr2_ana)
               ptrr2_ana(:,1:ngrid)=ptrr2_ana(:,1:ngrid)+ptrr2_ctl(:,1:ngrid)
               ptrr2_ctl(:,1:ngrid)=mpas_jedi_zero_kr
            else
               write(message,*) '--> mpasjedi_linvarcha_c2a::multiplyinverse: poolItr % nDims == ',poolitr % nDims,' not handled'
               call abor1_ftn(message)
            end if
         end if
      end if
   end do !- end of pool iteration
 
   if( xctl%has('relhum') .and. xana%has('spechum') ) then
      call xana%get(    'spechum', ptrr2_sh)
      call xctl%get(     'relhum', ptrr2_rh)
      call mpas_pool_get_field(self%trajectories, 'spechum_sat', fld2d_traj_qs)
      call pseudorh_to_spechum_inversead( ptrr2_rh(:,1:ngrid), ptrr2_sh(:,1:ngrid), fld2d_traj_qs%array(:,1:ngrid) )
   end if
 
   !-- dummy inverseAD operator
   if( xctl%has('stream_function') .and. xctl%has('velocity_potential') .and. &
       xana%has('uReconstructZonal') .and. xana%has('uReconstructMeridional') ) then
      call xana%get('uReconstructZonal', ptrr2_ana)
      call xctl%get(  'stream_function', ptrr2_ctl)
      ptrr2_ana(:,1:ngrid)=ptrr2_ana(:,1:ngrid)+ptrr2_ctl(:,1:ngrid)
      ptrr2_ctl(:,1:ngrid)=mpas_jedi_zero_kr
 
      call xana%get('uReconstructMeridional', ptrr2_ana)
      call xctl%get(    'velocity_potential', ptrr2_ctl)
      ptrr2_ana(:,1:ngrid)=ptrr2_ana(:,1:ngrid)+ptrr2_ctl(:,1:ngrid)
      ptrr2_ctl(:,1:ngrid)=mpas_jedi_zero_kr
   end if
 
end subroutine multiplyinverseadjoint
 
!-------------------------------------------------------------------------------
! Adjoint code of subroutine mpas_reconstruct()
!                 in MPAS-Model/src/operators/mpas_vector_reconstruction.F
!
subroutine mpas_reconstruct_2dad(meshPool, u, uReconstructX, uReconstructY, uReconstructZ,&
                                 uReconstructZonal, uReconstructMeridional, includeHalos)!{{{
 
    implicit none
 
    type (mpas_pool_type), intent(in) :: meshPool !< Input: Mesh information
    real (kind=rkind), dimension(:,:), intent(inout) :: u !< Input: Velocity field on edges
    real (kind=rkind), dimension(:,:), intent(inout) :: ureconstructx !< Output: X Component of velocity reconstructed to cell centers
    real (kind=rkind), dimension(:,:), intent(inout) :: ureconstructy !< Output: Y Component of velocity reconstructed to cell centers
    real (kind=rkind), dimension(:,:), intent(inout) :: ureconstructz !< Output: Z Component of velocity reconstructed to cell centers
    real (kind=rkind), dimension(:,:), intent(inout) :: ureconstructzonal !< Output: Zonal Component of velocity reconstructed to cell centers
    real (kind=rkind), dimension(:,:), intent(inout) :: ureconstructmeridional !< Output: Meridional Component of velocity reconstructed to cell centers
    logical, optional, intent(in) :: includeHalos !< Input: Optional logical that allows reconstruction over halo regions
 
    !   temporary arrays needed in the compute procedure
    logical :: includeHalosLocal
    integer, pointer :: nCells
    integer, dimension(:,:), pointer :: edgesOnCell
    integer, dimension(:), pointer :: nEdgesOnCell
    integer :: iCell,iEdge, i
    real(kind=rkind), dimension(:), pointer :: latcell, loncell
 
    real (kind=rkind), dimension(:,:,:), pointer :: coeffs_reconstruct
 
    logical, pointer :: on_a_sphere
 
    real (kind=rkind) :: clat, slat, clon, slon
 
    if ( present(includehalos) ) then
       includehaloslocal = includehalos
    else
       includehaloslocal = .false.
    end if
 
    ! stored arrays used during compute procedure
    call mpas_pool_get_array(meshpool, 'coeffs_reconstruct', coeffs_reconstruct)
 
    ! temporary variables
    call mpas_pool_get_array(meshpool, 'nEdgesOnCell', nedgesoncell)
    call mpas_pool_get_array(meshpool, 'edgesOnCell', edgesoncell)
 
    if ( includehaloslocal ) then
       call mpas_pool_get_dimension(meshpool, 'nCells', ncells)
    else
       call mpas_pool_get_dimension(meshpool, 'nCellsSolve', ncells)
    end if
 
    call mpas_pool_get_array(meshpool, 'latCell', latcell)
    call mpas_pool_get_array(meshpool, 'lonCell', loncell)
 
    call mpas_pool_get_config(meshpool, 'on_a_sphere', on_a_sphere)
 
    ! initialize the reconstructed vectors
    ureconstructx = mpas_jedi_zero_kr
    ureconstructy = mpas_jedi_zero_kr
    ureconstructz = mpas_jedi_zero_kr
    if (on_a_sphere) then
      do icell = 1, ncells
        clat = cos(latcell(icell))
        slat = sin(latcell(icell))
        clon = cos(loncell(icell))
        slon = sin(loncell(icell))
 
        ureconstructx(:,icell) = ureconstructx(:,icell) - clon*slat * ureconstructmeridional(:,icell)
        ureconstructy(:,icell) = ureconstructy(:,icell) - slon*slat * ureconstructmeridional(:,icell)
        ureconstructz(:,icell) = ureconstructz(:,icell) + clat * ureconstructmeridional(:,icell)
 
        ureconstructx(:,icell) = ureconstructx(:,icell) - slon * ureconstructzonal(:,icell)
        ureconstructy(:,icell) = ureconstructy(:,icell) + clon * ureconstructzonal(:,icell)
      end do
    else
      do icell = 1, ncells
        ureconstructy(:,icell) = ureconstructy(:,icell) + ureconstructmeridional(:,icell)
        ureconstructx(:,icell) = ureconstructx(:,icell) + ureconstructzonal(:,icell)
      end do
    end if
 
    ! loop over cell centers
    do icell = 1, ncells
      ! a more efficient reconstruction where rbf_values*matrix_reconstruct
      ! has been precomputed in coeffs_reconstruct
      do i=1,nedgesoncell(icell)
        iedge = edgesoncell(i,icell)
        u(:,iedge) = u(:,iedge) + coeffs_reconstruct(1,i,icell) * ureconstructx(:,icell) &
                                + coeffs_reconstruct(2,i,icell) * ureconstructy(:,icell) &
                                + coeffs_reconstruct(3,i,icell) * ureconstructz(:,icell)
      enddo
    enddo   ! iCell
 
end subroutine mpas_reconstruct_2dad!}}}
 
!-------------------------------------------------------------------------------
 
subroutine mpas_reconstruct_1dad(meshPool, u, uReconstructX, uReconstructY, uReconstructZ,&
                                 uReconstructZonal, uReconstructMeridional, includeHalos)!{{{
 
    implicit none
 
    type (mpas_pool_type), intent(in) :: meshPool !< Input: Mesh information
    real (kind=rkind), dimension(:), intent(inout) :: u !< Input: Velocity field on edges
    real (kind=rkind), dimension(:), intent(inout) :: ureconstructx !< Output: X Component of velocity reconstructed to cell centers
    real (kind=rkind), dimension(:), intent(inout) :: ureconstructy !< Output: Y Component of velocity reconstructed to cell centers
    real (kind=rkind), dimension(:), intent(inout) :: ureconstructz !< Output: Z Component of velocity reconstructed to cell centers
    real (kind=rkind), dimension(:), intent(inout) :: ureconstructzonal !< Output: Zonal Component of velocity reconstructed to cell centers
    real (kind=rkind), dimension(:), intent(inout) :: ureconstructmeridional !< Output: Meridional Component of velocity reconstructed to cell centers
    logical, optional, intent(in) :: includeHalos !< Input: Optional logical that allows reconstruction over halo regions
 
    !   temporary arrays needed in the compute procedure
    logical :: includeHalosLocal
    integer, pointer :: nCells
    integer, dimension(:,:), pointer :: edgesOnCell
    integer, dimension(:), pointer :: nEdgesOnCell
    integer :: iCell,iEdge, i
    real(kind=rkind), dimension(:), pointer :: latcell, loncell
 
    real (kind=rkind), dimension(:,:,:), pointer :: coeffs_reconstruct
 
    logical, pointer :: on_a_sphere
 
    real (kind=rkind) :: clat, slat, clon, slon
 
    if ( present(includehalos) ) then
       includehaloslocal = includehalos
    else
       includehaloslocal = .false.
    end if
 
    ! stored arrays used during compute procedure
    call mpas_pool_get_array(meshpool, 'coeffs_reconstruct', coeffs_reconstruct)
 
    ! temporary variables
    call mpas_pool_get_array(meshpool, 'nEdgesOnCell', nedgesoncell)
    call mpas_pool_get_array(meshpool, 'edgesOnCell', edgesoncell)
 
    if ( includehaloslocal ) then
       call mpas_pool_get_dimension(meshpool, 'nCells', ncells)
    else
       call mpas_pool_get_dimension(meshpool, 'nCellsSolve', ncells)
    end if
 
    call mpas_pool_get_array(meshpool, 'latCell', latcell)
    call mpas_pool_get_array(meshpool, 'lonCell', loncell)
 
    call mpas_pool_get_config(meshpool, 'on_a_sphere', on_a_sphere)
 
    ! initialize the reconstructed vectors
    ureconstructx = mpas_jedi_zero_kr
    ureconstructy = mpas_jedi_zero_kr
    ureconstructz = mpas_jedi_zero_kr
    if (on_a_sphere) then
      do icell = 1, ncells
        clat = cos(latcell(icell))
        slat = sin(latcell(icell))
        clon = cos(loncell(icell))
        slon = sin(loncell(icell))
 
        ureconstructx(icell) = ureconstructx(icell) - clon*slat * ureconstructmeridional(icell)
        ureconstructy(icell) = ureconstructy(icell) - slon*slat * ureconstructmeridional(icell)
        ureconstructz(icell) = ureconstructz(icell) + clat * ureconstructmeridional(icell)
 
        ureconstructx(icell) = ureconstructx(icell) - slon * ureconstructzonal(icell)
        ureconstructy(icell) = ureconstructy(icell) + clon * ureconstructzonal(icell)
      end do
    else
      do icell = 1, ncells
        ureconstructy(icell) = ureconstructy(icell) + ureconstructmeridional(icell)
        ureconstructx(icell) = ureconstructx(icell) + ureconstructzonal(icell)
      end do
    end if
 
    ! loop over cell centers
    do icell = 1, ncells
      ! a more efficient reconstruction where rbf_values*matrix_reconstruct
      ! has been precomputed in coeffs_reconstruct
      do i=1,nedgesoncell(icell)
        iedge = edgesoncell(i,icell)
        u(iedge) = u(iedge) + coeffs_reconstruct(1,i,icell) * ureconstructx(icell) &
                                + coeffs_reconstruct(2,i,icell) * ureconstructy(icell) &
                                + coeffs_reconstruct(3,i,icell) * ureconstructz(icell)
      enddo
    enddo   ! iCell
 
end subroutine mpas_reconstruct_1dad!}}}
 
!-------------------------------------------------------------------------------
! input  : psi & chi @ cell center
! output : u & v @ cell center
subroutine psichi_to_uv_center(geom, psi, chi, u, v)
 
   implicit none
   type (mpas_geom),                                               intent(in)  :: geom         !< geometry
   real (kind=kind_real), dimension(geom%nVertLevels,geom%nCells), intent(in)  :: psi, chi
   real (kind=kind_real), dimension(geom%nVertLevels,geom%nCells), intent(out) :: u, v
 
   integer :: iC, iE, j
   real (kind=kind_real), dimension(geom%nVertLevels,geom%nCellsSolve) :: &
              psi_line_intg_dx, psi_line_intg_dy, chi_line_intg_dx, chi_line_intg_dy
 
   psi_line_intg_dx=mpas_jedi_zero_kr
   psi_line_intg_dy=mpas_jedi_zero_kr
   chi_line_intg_dx=mpas_jedi_zero_kr
   chi_line_intg_dy=mpas_jedi_zero_kr
 
   do ic = 1, geom%nCellsSolve
     do j = 1, geom%nEdgesOnCell(ic) ! or geom%maxEdges
       ie = geom%edgesOnCell(j,ic)
       psi_line_intg_dx(:,ic) = psi_line_intg_dx(:,ic) &
          + mpas_jedi_half_kr * ( psi(:,geom%cellsOnEdge(1,ie)) + psi(:,geom%cellsOnEdge(2,ie)) ) &
                              * geom%dvEdge(ie) * sin(mpas_jedi_zero_kr-geom%angleEdge(ie)) * geom%edgesOnCell_sign(j,ic)
       psi_line_intg_dy(:,ic) = psi_line_intg_dy(:,ic) &
          + mpas_jedi_half_kr * ( psi(:,geom%cellsOnEdge(1,ie)) + psi(:,geom%cellsOnEdge(2,ie)) ) &
                              * geom%dvEdge(ie) * cos(geom%angleEdge(ie)) * geom%edgesOnCell_sign(j,ic)
       chi_line_intg_dx(:,ic) = chi_line_intg_dx(:,ic) &
          + mpas_jedi_half_kr * ( chi(:,geom%cellsOnEdge(1,ie)) + chi(:,geom%cellsOnEdge(2,ie)) ) &
                              * geom%dvEdge(ie) * sin(mpas_jedi_zero_kr-geom%angleEdge(ie)) * geom%edgesOnCell_sign(j,ic)
       chi_line_intg_dy(:,ic) = chi_line_intg_dy(:,ic) &
          + mpas_jedi_half_kr * ( chi(:,geom%cellsOnEdge(1,ie)) + chi(:,geom%cellsOnEdge(2,ie)) ) &
                              * geom%dvEdge(ie) * cos(geom%angleEdge(ie)) * geom%edgesOnCell_sign(j,ic)
     enddo !- j
   enddo !- iC
 
   do ic=1, geom%nCellsSolve
     u(:,ic) = ( psi_line_intg_dx(:,ic) - chi_line_intg_dy(:,ic) ) / geom%areaCell(ic)
     v(:,ic) = ( psi_line_intg_dy(:,ic) + chi_line_intg_dx(:,ic) ) / geom%areaCell(ic)
   enddo
 
end subroutine psichi_to_uv_center
 
!-------------------------------------------------------------------------------
 
subroutine psichi_to_uv_centerad(geom, psi, chi, u, v)
 
   implicit none
   type (mpas_geom),                                               intent(in)    :: geom         !< geometry
   real (kind=kind_real), dimension(geom%nVertLevels,geom%nCells), intent(inout) :: psi, chi
   real (kind=kind_real), dimension(geom%nVertLevels,geom%nCells), intent(inout) :: u, v
 
   integer :: iC, iE, j
   real (kind=kind_real), dimension(geom%nVertLevels,geom%nCellsSolve) :: &
              psi_line_intg_dx, psi_line_intg_dy, chi_line_intg_dx, chi_line_intg_dy
 
   psi_line_intg_dx=mpas_jedi_zero_kr
   psi_line_intg_dy=mpas_jedi_zero_kr
   chi_line_intg_dx=mpas_jedi_zero_kr
   chi_line_intg_dy=mpas_jedi_zero_kr
 
   do ic=1, geom%nCellsSolve
     psi_line_intg_dx(:,ic) = psi_line_intg_dx(:,ic) + u(:,ic) / geom%areaCell(ic)
     chi_line_intg_dy(:,ic) = chi_line_intg_dy(:,ic) - u(:,ic) / geom%areaCell(ic)
     psi_line_intg_dy(:,ic) = psi_line_intg_dy(:,ic) + v(:,ic) / geom%areaCell(ic)
     chi_line_intg_dx(:,ic) = chi_line_intg_dx(:,ic) + v(:,ic) / geom%areaCell(ic)
   enddo
 
   do ic = 1, geom%nCellsSolve
     do j = 1, geom%nEdgesOnCell(ic) ! or geom%maxEdges
       ie = geom%edgesOnCell(j,ic)
 
       psi(:,geom%cellsOnEdge(1,ie)) = psi(:,geom%cellsOnEdge(1,ie)) + mpas_jedi_half_kr * geom%dvEdge(ie)&
 * sin(mpas_jedi_zero_kr-geom%angleEdge(ie)) * geom%edgesOnCell_sign(j,ic) * psi_line_intg_dx(:,ic)
       psi(:,geom%cellsOnEdge(2,ie)) = psi(:,geom%cellsOnEdge(2,ie)) + mpas_jedi_half_kr * geom%dvEdge(ie)&
 * sin(mpas_jedi_zero_kr-geom%angleEdge(ie)) * geom%edgesOnCell_sign(j,ic) * psi_line_intg_dx(:,ic)
       
       psi(:,geom%cellsOnEdge(1,ie)) = psi(:,geom%cellsOnEdge(1,ie)) + mpas_jedi_half_kr * geom%dvEdge(ie)&
 * cos(geom%angleEdge(ie)) * geom%edgesOnCell_sign(j,ic) * psi_line_intg_dy(:,ic)
       psi(:,geom%cellsOnEdge(2,ie)) = psi(:,geom%cellsOnEdge(2,ie)) + mpas_jedi_half_kr * geom%dvEdge(ie)&
 * cos(geom%angleEdge(ie)) * geom%edgesOnCell_sign(j,ic) * psi_line_intg_dy(:,ic)
 
       chi(:,geom%cellsOnEdge(1,ie)) = chi(:,geom%cellsOnEdge(1,ie)) + mpas_jedi_half_kr * geom%dvEdge(ie)&
 * sin(mpas_jedi_zero_kr-geom%angleEdge(ie)) * geom%edgesOnCell_sign(j,ic) * chi_line_intg_dx(:,ic)
       chi(:,geom%cellsOnEdge(2,ie)) = chi(:,geom%cellsOnEdge(2,ie)) + mpas_jedi_half_kr * geom%dvEdge(ie)&
 * sin(mpas_jedi_zero_kr-geom%angleEdge(ie)) * geom%edgesOnCell_sign(j,ic) * chi_line_intg_dx(:,ic)
 
       chi(:,geom%cellsOnEdge(1,ie)) = chi(:,geom%cellsOnEdge(1,ie)) + mpas_jedi_half_kr * geom%dvEdge(ie)&
 * cos(geom%angleEdge(ie)) * geom%edgesOnCell_sign(j,ic) * chi_line_intg_dy(:,ic)
       chi(:,geom%cellsOnEdge(2,ie)) = chi(:,geom%cellsOnEdge(2,ie)) + mpas_jedi_half_kr * geom%dvEdge(ie)&
 * cos(geom%angleEdge(ie)) * geom%edgesOnCell_sign(j,ic) * chi_line_intg_dy(:,ic)
     enddo !- j
   enddo !- iC
 
end subroutine psichi_to_uv_centerad
 
!-------------------------------------------------------------------------------
! input  : psi @ cell center
! output : psi @ vertices
subroutine psichi_to_uv_edge_step1(geom, psi, psi_v)
 
   implicit none
   type (mpas_geom),                                                  intent(in)    :: geom         !< geometry
   real (kind=kind_real), dimension(geom%nVertLevels,geom%nCells),    intent(in)    :: psi
   real (kind=kind_real), dimension(geom%nVertLevels,geom%nVertices), intent(inout) :: psi_v
 
   integer :: iC, iV, j
 
   psi_v=mpas_jedi_zero_kr
 
   ! Interpolate psi in cell center to vertice
   do iv = 1, geom%nVerticesSolve ! local
     do j = 1, geom%vertexDegree
       ic = geom%cellsOnVertex(j,iv)
       psi_v(:,iv) = psi_v(:,iv) + geom%kiteAreasOnVertex(j,iv) * psi(:,ic)
     enddo
     psi_v(:,iv) = psi_v(:,iv) / geom%areaTriangle(iv)
   enddo
 
end subroutine psichi_to_uv_edge_step1
 
!-------------------------------------------------------------------------------
! input  : psi @ vertices, chi @ cell center
! output : edge_normal_wind @ edges
subroutine psichi_to_uv_edge_step2(geom, psi_v, chi, edge_normal_wind)
 
   implicit none
   type (mpas_geom),                                                  intent(in)    :: geom         !< geometry
   real (kind=kind_real), dimension(geom%nVertLevels,geom%nVertices), intent(in)    :: psi_v
   real (kind=kind_real), dimension(geom%nVertLevels,geom%nCells),    intent(in)    :: chi
   real (kind=kind_real), dimension(geom%nVertLevels,geom%nEdges),    intent(inout) :: edge_normal_wind
 
   integer :: iE
 
   edge_normal_wind=mpas_jedi_zero_kr
 
   !get edge_normal_wind
   do ie = 1, geom%nEdgesSolve ! local
     edge_normal_wind(:,ie) = edge_normal_wind(:,ie) - &
                ( chi(:,geom%cellsOnEdge(2,ie)) - chi(:,geom%cellsOnEdge(1,ie)) ) / geom%dcEdge(ie) - &
                ( psi_v(:,geom%verticesOnEdge(2,ie)) - psi_v(:,geom%verticesOnEdge(1,ie)) ) / geom%dvEdge(ie)
   enddo
 
end subroutine psichi_to_uv_edge_step2
 
!-------------------------------------------------------------------------------
! input  : edge_normal_wind @ edges
! output : zonal- and meridional- wind @ cell center
subroutine psichi_to_uv_edge_step3(geom, edge_normal_wind, u, v)
 
   use mpas_vector_reconstruction
 
   implicit none
   type (mpas_geom),                                               intent(in)    :: geom         !< geometry
   real (kind=kind_real), dimension(geom%nVertLevels,geom%nEdges), intent(in)    :: edge_normal_wind
   real (kind=kind_real), dimension(geom%nVertLevels,geom%nCells), intent(inout) :: u, v
 
   real (kind=kind_real), dimension(:,:), allocatable :: &
                                        ureconstructx, ureconstructy, ureconstructz
   type (mpas_pool_type), pointer :: mesh
 
   allocate(ureconstructx(geom%nVertLevels,geom%nCells))
   allocate(ureconstructy(geom%nVertLevels,geom%nCells))
   allocate(ureconstructz(geom%nVertLevels,geom%nCells))
 
   !edge_wind -> u/v @ cell center
   call mpas_pool_get_subpool( geom % domain % blocklist % structs, 'mesh', mesh)
   call mpas_reconstruct(mesh, edge_normal_wind,    &
                         ureconstructx,             &
                         ureconstructy,             &
                         ureconstructz,             &
                         u,                         &
                         v, .false. ) ! local only, no halo calculation
 
   deallocate(ureconstructx, ureconstructy, ureconstructz)
 
end subroutine psichi_to_uv_edge_step3
 
!-------------------------------------------------------------------------------
subroutine psichi_to_uv_edge_step1ad(geom, psi, psi_v)
 
   implicit none
   type (mpas_geom),                                                  intent(in)    :: geom         !< geometry
   real (kind=kind_real), dimension(geom%nVertLevels,geom%nCells),    intent(inout) :: psi
   real (kind=kind_real), dimension(geom%nVertLevels,geom%nVertices), intent(inout) :: psi_v
 
   integer :: iC, iV, j
 
   psi=mpas_jedi_zero_kr
 
   ! Interpolate psi in cell center to vertice
   do iv = 1, geom%nVerticesSolve ! local
     psi_v(:,iv) = psi_v(:,iv) / geom%areaTriangle(iv)
     do j = 1, geom%vertexDegree
       ic = geom%cellsOnVertex(j,iv)
       psi(:,ic) = psi(:,ic) + geom%kiteAreasOnVertex(j,iv) * psi_v(:,iv)
     enddo
   enddo
 
end subroutine psichi_to_uv_edge_step1ad
 
!-------------------------------------------------------------------------------
subroutine psichi_to_uv_edge_step2ad(geom, psi_v, chi, edge_normal_wind)
 
   implicit none
   type (mpas_geom),                                                  intent(in)    :: geom         !< geometry
   real (kind=kind_real), dimension(geom%nVertLevels,geom%nVertices), intent(inout) :: psi_v
   real (kind=kind_real), dimension(geom%nVertLevels,geom%nCells),    intent(inout) :: chi
   real (kind=kind_real), dimension(geom%nVertLevels,geom%nEdges),    intent(inout) :: edge_normal_wind
 
   integer :: iE
 
   chi=mpas_jedi_zero_kr
   psi_v=mpas_jedi_zero_kr
 
   !get edge_normal_wind
   do ie = 1, geom%nEdgesSolve ! local
     chi(:,geom%cellsOnEdge(2,ie))      = chi(:,geom%cellsOnEdge(2,ie))      - edge_normal_wind(:,ie) / geom%dcEdge(ie)
     chi(:,geom%cellsOnEdge(1,ie))      = chi(:,geom%cellsOnEdge(1,ie))      + edge_normal_wind(:,ie) / geom%dcEdge(ie)
     psi_v(:,geom%verticesOnEdge(2,ie)) = psi_v(:,geom%verticesOnEdge(2,ie)) - edge_normal_wind(:,ie) / geom%dvEdge(ie)
     psi_v(:,geom%verticesOnEdge(1,ie)) = psi_v(:,geom%verticesOnEdge(1,ie)) + edge_normal_wind(:,ie) / geom%dvEdge(ie)
   enddo
 
end subroutine psichi_to_uv_edge_step2ad
 
!-------------------------------------------------------------------------------
subroutine psichi_to_uv_edge_step3ad(geom, edge_normal_wind, u, v)
 
   use mpas_vector_reconstruction
 
   implicit none
   type (mpas_geom),                                               intent(in)    :: geom         !< geometry
   real (kind=kind_real), dimension(geom%nVertLevels,geom%nEdges), intent(inout) :: edge_normal_wind
   real (kind=kind_real), dimension(geom%nVertLevels,geom%nCells), intent(inout) :: u, v
 
   real (kind=kind_real), dimension(:,:), allocatable :: &
                                        ureconstructx, ureconstructy, ureconstructz
   type (mpas_pool_type), pointer :: mesh
 
   edge_normal_wind=mpas_jedi_zero_kr
 
   allocate(ureconstructx(geom%nVertLevels,geom%nCells))
   allocate(ureconstructy(geom%nVertLevels,geom%nCells))
   allocate(ureconstructz(geom%nVertLevels,geom%nCells))
 
   call mpas_pool_get_subpool( geom % domain % blocklist % structs, 'mesh', mesh)
   call mpas_reconstruct_2dad(mesh, edge_normal_wind, &
                         ureconstructx,             &
                         ureconstructy,             &
                         ureconstructz,             &
                         u,                         &
                         v, .false. ) ! local only, no halo calculation
 
   deallocate(ureconstructx, ureconstructy, ureconstructz)
 
end subroutine psichi_to_uv_edge_step3ad
 
! ------------------------------------------------------------------------------
 
elemental subroutine pseudorh_to_spechum(pseudorh,spechum,saturation_spechum)
 
   implicit none
   real (kind=kind_real), intent(in)  :: pseudorh
   real (kind=kind_real), intent(out) :: spechum
   real (kind=kind_real), intent(in)  :: saturation_spechum
 
   spechum = pseudorh * saturation_spechum / mpas_jedi_hundred_kr
 
end subroutine pseudorh_to_spechum
 
elemental subroutine pseudorh_to_spechumad(pseudorh,spechum,saturation_spechum)
 
   implicit none
   real (kind=kind_real), intent(inout) :: pseudorh
   real (kind=kind_real), intent(inout) :: spechum
   real (kind=kind_real), intent(in)    :: saturation_spechum
 
   pseudorh = pseudorh + spechum * saturation_spechum / mpas_jedi_hundred_kr
   spechum = mpas_jedi_zero_kr
 
end subroutine pseudorh_to_spechumad
 
elemental subroutine pseudorh_to_spechum_inverse(pseudorh,spechum,saturation_spechum)
 
   implicit none
   real (kind=kind_real), intent(out) :: pseudorh
   real (kind=kind_real), intent(in)  :: spechum
   real (kind=kind_real), intent(in)  :: saturation_spechum
 
   pseudorh = spechum / saturation_spechum * mpas_jedi_hundred_kr
 
end subroutine pseudorh_to_spechum_inverse
 
elemental subroutine pseudorh_to_spechum_inversead(pseudorh,spechum,saturation_spechum)
 
   implicit none
   real (kind=kind_real), intent(inout) :: pseudorh
   real (kind=kind_real), intent(inout) :: spechum
   real (kind=kind_real), intent(in)    :: saturation_spechum
 
   spechum = spechum + pseudorh / saturation_spechum * mpas_jedi_hundred_kr
   pseudorh = mpas_jedi_zero_kr
 
end subroutine pseudorh_to_spechum_inversead
 
elemental subroutine da_tp_to_qs( t, p, qs)
 
   !---------------------------------------------------------------------------
   ! Purpose: Convert T/p to saturation specific humidity.
   !
   !  Method: qs = es_alpha * es / ( p - ( 1 - rd_over_rv ) * es ).
   !          use Rogers & Yau (1989) formula: es = a exp( bTc / (T_c + c) )
   !--------------------------------------------------------------------------
 
   implicit none
 
   real (kind=kind_real), intent(in)  :: t, p
   real (kind=kind_real), intent(out) :: qs
 
   real (kind=kind_real) :: es
   real (kind=kind_real) :: t_c              ! T in degreesC.
 
   !---------------------------------------------------------------------------
   ! [1.0] initialise:
   !---------------------------------------------------------------------------
   t_c = t - t_kelvin
 
   !---------------------------------------------------------------------------
   ! [2.0] Calculate saturation vapour pressure:
   !---------------------------------------------------------------------------
   es = es_alpha * exp( es_beta * t_c / ( t_c + es_gamma ) )
 
   !---------------------------------------------------------------------------
   ! [3.0] Calculate saturation specific humidity:
   !---------------------------------------------------------------------------
   qs = rd_over_rv * es / ( p - rd_over_rv1 * es )
 
end subroutine da_tp_to_qs
 
end module mpasjedi_linvarcha_c2a_mod