mpasjedi_lvc_model2geovars_mod.F90

Go to the documentation of this file.

! (C) Copyright 2020 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_lvc_model2geovars_mod
 
use iso_c_binding
 
use fckit_configuration_module, only: fckit_configuration
use fckit_log_module, only: fckit_log
 
!oops
use kinds, only : kind_real
 
!ufo
use gnssro_mod_transform, only: geometric2geop
use ufo_vars_mod
 
!MPAS-Model
use mpas_derived_types
use mpas_field_routines
use mpas_kind_types, only: rkind
use mpas_pool_routines
 
!mpas-jedi
use mpas_constants_mod
use mpas_fields_mod, only: mpas_fields
use mpas_geom_mod, only: mpas_geom
use mpas4da_mod, only: da_template_pool
 
!TODO: package the variable changes somewhere
use mpas2ufo_vars_mod
!use mpas_vc_utils
 
implicit none
 
private
public :: mpasjedi_lvc_model2geovars
 
type :: mpasjedi_lvc_model2geovars
  type(mpas_pool_type), pointer :: trajectory => null()
  contains
    procedure, public :: create
    procedure, public :: delete
    procedure, public :: multiply
    procedure, public :: multiplyadjoint
end type mpasjedi_lvc_model2geovars
 
! --------------------------------------------------------------------------------------------------
 
character(len=1024) :: message
 
contains
 
! --------------------------------------------------------------------------------------------------
 
subroutine create(self, geom, bg, fg, conf)
 
  class(mpasjedi_lvc_model2geovars), intent(inout) :: self
  type(mpas_geom),                   intent(in)    :: geom
  type(mpas_fields),                 intent(in)    :: bg
  type(mpas_fields),                 intent(in)    :: fg
  type(fckit_configuration),         intent(in)    :: conf
 
  integer :: iVar
 
  character(len=MAXVARLEN), parameter :: &
    trajFieldNames(4) = &
      [ character(len=maxvarlen) :: &
        'temperature', &
        'spechum', &
        'pressure', &
        'surface_pressure' &
      ]
 
  call da_template_pool(geom, self%trajectory, size(trajfieldnames), trajfieldnames)
 
  do ivar = 1, size(trajfieldnames)
    call bg%copy_to(trajfieldnames(ivar), self%trajectory)
  end do
 
end subroutine create
 
! --------------------------------------------------------------------------------------------------
 
subroutine delete(self)
  class(mpasjedi_lvc_model2geovars), intent(inout) :: self
  if (associated(self%trajectory)) then
    call mpas_pool_destroy_pool(self%trajectory)
  end if
end subroutine delete
 
! --------------------------------------------------------------------------------------------------
 
subroutine multiply(self, geom, dxm, dxg)
  class(mpasjedi_lvc_model2geovars), intent(inout) :: self
  type(mpas_geom),                   intent(inout) :: geom !< mpas mesh descriptors
  class(mpas_fields),                intent(in)    :: dxm  !< model increment fields
  class(mpas_fields),                intent(inout) :: dxg  !< increment containing linear geovar fields
 
  ! pool-related pointers
  type(mpas_pool_type), pointer :: mFields_tl, gFields_tl
  type(mpas_pool_data_type), pointer :: gdata
 
  ! reusable arrays
  real(kind=kind_real), dimension(:), pointer :: ptrr1_a
  real(kind=kind_real), dimension(:,:), pointer :: ptrr2_a, ptrr2_b, traj_ptrr2_a, traj_ptrr2_b
  real(kind=kind_real), dimension(:,:), allocatable :: r2, trajr2
 
  ! iteration-specific variables
  character(len=MAXVARLEN) :: geovar
  integer :: nCells, nVertLevels, nVertLevelsP1
  integer :: iVar
 
  ! air pressure on w levels
  real(kind=kind_real), allocatable :: plevels(:,:)
 
  ! convenient local variables
  mfields_tl => dxm % subFields
  gfields_tl => dxg % subFields
  ncells = geom%nCellsSolve
  nvertlevels = geom%nVertLevels
  nvertlevelsp1 = geom%nVertLevelsP1
 
  ! pre-calculate pressure on w levels
  allocate(plevels(1:nvertlevelsp1,1:ncells))
  call mpas_pool_get_array(self%trajectory, 'pressure', ptrr2_a)
  call mpas_pool_get_array(self%trajectory, 'surface_pressure', ptrr1_a)
  call pressure_half_to_full(ptrr2_a(:,1:ncells), geom%zgrid(:,1:ncells), ptrr1_a(1:ncells), &
                             ncells, nvertlevels, plevels)
 
  ! populate dxg
  do ivar = 1, dxg % nf
 
    geovar = trim(dxg % fldnames(ivar))
 
    if (geom%has_identity(geovar)) then
      ! "identity" variable changes are controlled in an external configuration file
      ! through the geom object. Refer to the mpas_geom type for more information.
      if (dxm%has(geom%identity(geovar))) then
        call dxm%copy_to(geom%identity(geovar), dxg, geovar)
      else
        !note: this warning is specifically for pressure => air_pressure in GNSSRO
        write(message,'(2A)') &
           'WARNING: mpasjedi_lvc_model2geovars::multiply: '&
          &'increment missing identity field for geovar => ', trim(geovar)
        call fckit_log%warning(message)
      end if
    else
 
      call dxg%get(geovar, gdata)
 
      select case (trim(geovar))
 
        case ( var_tv ) !-virtual_temperature
          ! get TL variables
          call dxm%get('temperature', ptrr2_a)
          call dxm%get('spechum', ptrr2_b)
 
          ! temporary work fields
          allocate(r2(1:nvertlevels,1:ncells))     ! TL mixing_ratio
          allocate(trajr2(1:nvertlevels,1:ncells)) ! NL mixing_ratio
 
          ! get linearization state
          call mpas_pool_get_array(self%trajectory, 'temperature', traj_ptrr2_a)
          call mpas_pool_get_array(self%trajectory, 'spechum', traj_ptrr2_b)
          call q_to_w(traj_ptrr2_b(:,1:ncells), trajr2(:,1:ncells)) !NL coeff.
 
          ! calculations
          call q_to_w_tl(ptrr2_b(:,1:ncells), traj_ptrr2_b(:,1:ncells), r2(:,1:ncells))
          call tw_to_tv_tl(ptrr2_a(:,1:ncells), r2(:,1:ncells), &
                           traj_ptrr2_a(:,1:ncells), trajr2(:,1:ncells), &
                           gdata%r2%array(:,1:ncells))
 
          ! cleanup
          deallocate(r2, trajr2)
 
        case ( var_mixr ) !-humidity_mixing_ratio
          ! get TL variables
          call dxm%get('spechum', ptrr2_a)
 
          ! get linearization state
          call mpas_pool_get_array(self%trajectory, 'spechum', traj_ptrr2_a)
 
          ! calculations
          call q_to_w_tl(ptrr2_a(:,1:ncells), traj_ptrr2_a(:,1:ncells), gdata%r2%array(:,1:ncells))
          gdata%r2%array(:,1:ncells) = gdata%r2%array(:,1:ncells) * mpas_jedi_thousand_kr
 
          ! Ensure positive-definite mixing ratios
          !  with respect to precision of crtm::CRTM_Parameters::ZERO.
          ! TODO: this should be moved to the mpas_fields%read step
          !       only other place it should be needed is add_incr (already there)
          where (traj_ptrr2_a(:,1:ncells) <= mpas_jedi_zero_kr)
            gdata%r2%array(:,1:ncells) = mpas_jedi_zero_kr
          end where
 
        case ( var_clw ) !-mass_content_of_cloud_liquid_water_in_atmosphere_layer
          call q_fields_tl('qc', mfields_tl, gdata%r2, plevels, ncells, nvertlevels)
 
        case ( var_cli ) !-mass_content_of_cloud_ice_in_atmosphere_layer
          call q_fields_tl('qi', mfields_tl, gdata%r2, plevels, ncells, nvertlevels)
 
        case ( var_clr ) !-mass_content_of_rain_in_atmosphere_layer
          call q_fields_tl('qr', mfields_tl, gdata%r2, plevels, ncells, nvertlevels)
 
        case ( var_cls ) !-mass_content_of_snow_in_atmosphere_layer
          call q_fields_tl('qs', mfields_tl, gdata%r2, plevels, ncells, nvertlevels)
 
        case ( var_clg ) !-mass_content_of_graupel_in_atmosphere_layer
          call q_fields_tl('qg', mfields_tl, gdata%r2, plevels, ncells, nvertlevels)
 
        case ( var_clh ) !-mass_content_of_hail_in_atmosphere_layer
          call q_fields_tl('qh', mfields_tl, gdata%r2, plevels, ncells, nvertlevels)
 
      end select
 
    end if
 
  end do !iVar
 
  deallocate(plevels)
 
 
end subroutine multiply
 
! --------------------------------------------------------------------------------------------------
 
subroutine multiplyadjoint(self, geom, dxg, dxm)
  class(mpasjedi_lvc_model2geovars), intent(inout) :: self
  type(mpas_geom),                   intent(inout) :: geom
  class(mpas_fields),                intent(in)    :: dxg
  class(mpas_fields),                intent(inout) :: dxm
 
  ! pool-related pointers
  type(mpas_pool_type), pointer :: mFields_ad, gFields_ad
  type(mpas_pool_data_type), pointer :: gdata
 
  ! reusable array pointers
  real(kind=kind_real), dimension(:), pointer :: ptrr1_a
  real(kind=kind_real), dimension(:,:), pointer :: ptrr2_a, ptrr2_b, traj_ptrr2_a, traj_ptrr2_b
  real(kind=kind_real), dimension(:,:), allocatable :: r2, trajr2
 
  ! iteration-specific variables
  character(len=MAXVARLEN) :: geovar
  integer :: nCells, nVertLevels, nVertLevelsP1
  integer :: iVar
 
  ! air pressure on w levels
  real(kind=kind_real), allocatable :: plevels(:,:)
 
  ! convenient local variables
  mfields_ad => dxm % subFields
  gfields_ad => dxg % subFields
  ncells = geom%nCellsSolve
  nvertlevels = geom%nVertLevels
  nvertlevelsp1 = geom%nVertLevelsP1
 
  ! pre-calculate pressure on w levels
  allocate(plevels(1:nvertlevelsp1,1:ncells))
  call mpas_pool_get_array(self%trajectory, 'pressure', ptrr2_a)
  call mpas_pool_get_array(self%trajectory, 'surface_pressure', ptrr1_a)
  call pressure_half_to_full(ptrr2_a(:,1:ncells), geom%zgrid(:,1:ncells), ptrr1_a(1:ncells), &
                             ncells, nvertlevels, plevels)
 
  ! populate dxg
  do ivar = 1, dxg % nf
 
    geovar = trim(dxg % fldnames(ivar))
 
    if (geom%has_identity(geovar)) then
      ! "identity" variable changes are controlled in an external configuration file
      ! through the geom object. Refer to the mpas_geom type for more information.
      if (dxm%has(geom%identity(geovar))) then
        call dxm%copy_to_ad(geom%identity(geovar), dxg, geovar)
      else
        !note: this warning is specifically for pressure => air_pressure in GNSSRO
        write(message,'(2A)') &
           'WARNING: mpasjedi_lvc_model2geovars::multiplyadjoint: '&
          &'increment missing identity field for geovar => ', trim(geovar)
        call fckit_log%warning(message)
      end if
    else
 
      call dxg%get(geovar, gdata)
 
      select case (trim(geovar))
 
        case ( var_tv ) !-virtual_temperature
          ! get AD variables
          call dxm%get('temperature', ptrr2_a)
          call dxm%get('spechum', ptrr2_b)
 
          ! temporary work fields
          allocate(r2(1:nvertlevels,1:ncells))     ! AD of mixing_ratio
          allocate(trajr2(1:nvertlevels,1:ncells)) ! NL mixing_ratio
 
          ! get linearization state
          call mpas_pool_get_array(self%trajectory, 'temperature', traj_ptrr2_a)
          call mpas_pool_get_array(self%trajectory, 'spechum', traj_ptrr2_b)
          call q_to_w(traj_ptrr2_b(:,1:ncells), trajr2(:,1:ncells)) !NL coeff.
 
          ! calculations
          r2(:,1:ncells) = mpas_jedi_zero_kr !initialize local var.
          call tw_to_tv_ad(ptrr2_a(:,1:ncells), r2(:,1:ncells), &
                           traj_ptrr2_a(:,1:ncells), trajr2(:,1:ncells), &
                           gdata%r2%array(:,1:ncells) )
          call q_to_w_ad(ptrr2_b(:,1:ncells), traj_ptrr2_b(:,1:ncells), r2(:,1:ncells))
 
          ! cleanup
          deallocate(r2, trajr2)
 
        case ( var_mixr ) !-humidity_mixing_ratio
          ! get AD variables
          call dxm%get('spechum', ptrr2_a)
 
          ! temporary work fields
          allocate(r2(1:nvertlevels,1:ncells)) ! AD of var_mixr
 
          ! get linearization state
          call mpas_pool_get_array(self%trajectory, 'spechum', traj_ptrr2_a)
 
          ! Ensure positive-definite mixing ratios
          !  with respect to precision of crtm::CRTM_Parameters::ZERO.
          ! TODO: this should be moved to the mpas_fields%read step
          !       only other place it should be needed is add_incr (already there)
          r2 = gdata%r2%array(1:nvertlevels,1:ncells)
          where (traj_ptrr2_a(:,1:ncells) <= mpas_jedi_zero_kr)
            r2(:,1:ncells) = mpas_jedi_zero_kr
          end where
 
          ! calculations
          r2(:,1:ncells) = r2(:,1:ncells) * mpas_jedi_thousand_kr
          call q_to_w_ad(ptrr2_a(:,1:ncells), traj_ptrr2_a(:,1:ncells), r2(:,1:ncells))
 
          ! cleanup
          deallocate(r2)
 
        case ( var_clw ) !-mass_content_of_cloud_liquid_water_in_atmosphere_layer
          call q_fields_ad('qc', mfields_ad, gdata%r2, plevels, ncells, nvertlevels)
 
        case ( var_cli ) !-mass_content_of_cloud_ice_in_atmosphere_layer
          call q_fields_ad('qi', mfields_ad, gdata%r2, plevels, ncells, nvertlevels)
 
        case ( var_clr ) !-mass_content_of_rain_in_atmosphere_layer
          call q_fields_ad('qr', mfields_ad, gdata%r2, plevels, ncells, nvertlevels)
 
        case ( var_cls ) !-mass_content_of_snow_in_atmosphere_layer
          call q_fields_ad('qs', mfields_ad, gdata%r2, plevels, ncells, nvertlevels)
 
        case ( var_clg ) !-mass_content_of_graupel_in_atmosphere_layer
          call q_fields_ad('qg', mfields_ad, gdata%r2, plevels, ncells, nvertlevels)
 
        case ( var_clh ) !-mass_content_of_hail_in_atmosphere_layer
          call q_fields_ad('qh', mfields_ad, gdata%r2, plevels, ncells, nvertlevels)
 
      end select
 
    end if
 
  end do !iVar
 
  deallocate(plevels)
 
end subroutine multiplyadjoint
 
! --------------------------------------------------------------------------------------------------
 
end module mpasjedi_lvc_model2geovars_mod
! --------------------------------------------------------------------------------------------------