fv_treat_da_inc_nlm.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> @brief Treat DA increment
!> @author Xi.Chen <xi.chen@noaa.gov>
!> @date 02/12/2016
!
!  REVISION HISTORY:
!  02/12/2016 - Initial Version
!-------------------------------------------------------------------------------

#ifdef OVERLOAD_R4
#define _GET_VAR1 get_var1_real 
#else
#define _GET_VAR1 get_var1_double
#endif

module fv_treat_da_inc_nlm_mod

  use fms_mod,           only: file_exist, read_data, &
                               field_exist, write_version_number
  use mpp_mod,           only: mpp_error, fatal, note, mpp_pe
  use mpp_domains_mod,   only: mpp_get_tile_id, &
                               domain2d, &
                               mpp_update_domains, &
                               north, &
                               east
  use tracer_manager_mod,only: get_tracer_names, &
                               get_number_tracers, &
                               get_tracer_index
  use field_manager_mod, only: model_atmos

#ifndef MAPL_MODE
  use constants_mod,     only: pi=>pi_8, omega, grav, kappa, &
                               rdgas, rvgas, cp_air
#else
  use mapl_constantsmod
#endif
  use fv_arrays_nlm_mod,     only: fv_atmos_type, &
                               fv_grid_type, &
                               fv_grid_bounds_type, &
                               r_grid,fvprc
  use fv_grid_utils_nlm_mod, only: ptop_min, g_sum, &
                               mid_pt_sphere, get_unit_vect2, &
                               get_latlon_vector, inner_prod, &
                               cubed_to_latlon
  use fv_mp_nlm_mod,         only: ng, &
                               is_master, &
                               fill_corners, &
                               ydir, &
                               mp_reduce_min, &
                               mp_reduce_max
#ifndef MAPL_MODE
  use sim_nc_nlm_mod,        only: open_ncfile, &
                               close_ncfile, &
                               get_ncdim1, &
                               get_var1_double, &
                               get_var2_real,   &
                               get_var3_r4, &
                               get_var1_real
#endif
  implicit none
  private

#ifdef MAPL_MODE
  real, parameter :: RADIUS       = mapl_radius
  real, parameter :: PI           = mapl_pi_r8
  real, parameter :: RDGAS        = mapl_rgas
  real, parameter :: GRAV         = mapl_grav
  real, parameter :: HLV          = mapl_alhl
  real, parameter :: CP_AIR       = mapl_cp
  real, parameter :: RVGAS        = mapl_rvap
  real, parameter :: OMEGA        = mapl_omega
  real, parameter :: KAPPA        = mapl_kappa
#endif

#ifndef MAPL_MODE
  public :: read_da_inc
#else
  public :: geos_get_da_increments
#endif

contains

  subroutine geos_get_da_increments(Atm, fv_domain, lon,lat,im,jm,km, &
                  u_amb, v_amb, t_amb, dp_amb, q_amb, o3_amb,  &
                  u_inc, v_inc, t_inc, dp_inc, q_inc, o3_inc) 
    type(fv_atmos_type), intent(inout) :: Atm(:)
    type(domain2d),      intent(inout) :: fv_domain
    real,                      intent(inout) :: lon(im), lat(jm)
    real, dimension(im,jm,km), intent(inout) :: u_amb, v_amb, t_amb, dp_amb, q_amb, o3_amb
    real, dimension(Atm(1)%bd%is:Atm(1)%bd%ie,Atm(1)%bd%js:Atm(1)%bd%je,Atm(1)%npz), intent(inout):: &
                                                u_inc, v_inc, t_inc, dp_inc, q_inc, o3_inc

    real, dimension(:,:)  , allocatable:: ut_inc, vt_inc
    real, dimension(:,:,:), allocatable:: ud_inc, vd_inc, ua_inc, va_inc

    real, allocatable:: pt_c(:,:,:), pt_d(:,:,:)
    real:: s2c(Atm(1)%bd%is:Atm(1)%bd%ie,Atm(1)%bd%js:Atm(1)%bd%je,4)
    real:: s2c_c(Atm(1)%bd%is:Atm(1)%bd%ie+1,Atm(1)%bd%js:Atm(1)%bd%je,4)
    real:: s2c_d(Atm(1)%bd%is:Atm(1)%bd%ie,Atm(1)%bd%js:Atm(1)%bd%je+1,4)
    integer, dimension(Atm(1)%bd%is:Atm(1)%bd%ie,Atm(1)%bd%js:Atm(1)%bd%je):: &
        id1, id2, jdc
    integer, dimension(Atm(1)%bd%is:Atm(1)%bd%ie+1,Atm(1)%bd%js:Atm(1)%bd%je)::&
        id1_c, id2_c, jdc_c
    integer, dimension(Atm(1)%bd%is:Atm(1)%bd%ie,Atm(1)%bd%js:Atm(1)%bd%je+1)::&
        id1_d, id2_d, jdc_d
    integer, parameter :: update_uv_ghost_cells=1
    real :: tmp1d(im)

    integer:: i, j, k, im, jm, km, npz
    integer:: i1, i2, j1, IMsplit
    real(kind=R_GRID), dimension(2):: p1, p2, p3
    real(kind=R_GRID), dimension(3):: e1, e2, ex, ey

    integer :: is,  ie,  js,  je
    integer :: isd, ied, jsd, jed

    is  = atm(1)%bd%is
    ie  = atm(1)%bd%ie
    js  = atm(1)%bd%js
    je  = atm(1)%bd%je
    isd = atm(1)%bd%isd
    ied = atm(1)%bd%ied
    jsd = atm(1)%bd%jsd
    jed = atm(1)%bd%jed
    npz = atm(1)%npz

  ! FV3 code wants lon 0:360
    imsplit = im/2
  ! Lons
    tmp1d(        1:imsplit) =          lon(imsplit+1:im     )    
    tmp1d(imsplit+1:im     ) = 2.0*pi + lon(        1:imsplit)
    lon = tmp1d
  ! ANA-BKG
    do k=1,km
      do j=1,jm
      ! U
        tmp1d(        1:imsplit) = u_amb(imsplit+1:im     ,j,k)
        tmp1d(imsplit+1:im     ) = u_amb(        1:imsplit,j,k)
        u_amb(:,j,k) = tmp1d
      ! V
        tmp1d(        1:imsplit) = v_amb(imsplit+1:im     ,j,k)
        tmp1d(imsplit+1:im     ) = v_amb(        1:imsplit,j,k)
        v_amb(:,j,k) = tmp1d
      ! T
        tmp1d(        1:imsplit) = t_amb(imsplit+1:im     ,j,k)
        tmp1d(imsplit+1:im     ) = t_amb(        1:imsplit,j,k)
        t_amb(:,j,k) = tmp1d
      ! DP
        tmp1d(        1:imsplit) = dp_amb(imsplit+1:im     ,j,k)
        tmp1d(imsplit+1:im     ) = dp_amb(        1:imsplit,j,k)
        dp_amb(:,j,k) = tmp1d
      ! Q
        tmp1d(        1:imsplit) = q_amb(imsplit+1:im     ,j,k)
        tmp1d(imsplit+1:im     ) = q_amb(        1:imsplit,j,k)
        q_amb(:,j,k) = tmp1d
      ! O3
        tmp1d(        1:imsplit) = o3_amb(imsplit+1:im     ,j,k)
        tmp1d(imsplit+1:im     ) = o3_amb(        1:imsplit,j,k)
        o3_amb(:,j,k) = tmp1d
      enddo
    enddo

    ! Initialize lat-lon to Cubed bi-linear interpolation coeff:
    call remap_coef( is, ie, js, je, isd, ied, jsd, jed, &
        im, jm, lon, lat, id1, id2, jdc, s2c, &
        atm(1)%gridstruct%agrid)

    ! perform increments on scalars
    call get_inc_on_3d_scalar( t_amb, t_inc)
    call get_inc_on_3d_scalar(dp_amb,dp_inc)
    call get_inc_on_3d_scalar( q_amb, q_inc)
    call get_inc_on_3d_scalar(o3_amb,o3_inc)

    ! perform increments on winds
    allocate (ud_inc(isd:ied  , jsd:jed+1, km))
    allocate (vd_inc(isd:ied+1, jsd:jed  , km))
    allocate (  pt_c(isd:ied+1, jsd:jed  ,  2))
    allocate (  pt_d(isd:ied  , jsd:jed+1,  2))
    allocate (ut_inc(is :ie+1 , js :je       ))
    allocate (vt_inc(is :ie+1 , js :je       ))

    call get_staggered_grid( &
        is, ie, js, je, &
        isd, ied, jsd, jed, &
        atm(1)%gridstruct%grid, pt_c, pt_d)

    !------ pt_c part ------
    ! Initialize lat-lon to Cubed bi-linear interpolation coeff:
    call remap_coef( is, ie+1, js, je, isd, ied+1, jsd, jed, &
        im, jm, lon, lat, id1_c, id2_c, jdc_c, s2c_c, &
        pt_c)

    do k=1,km
      do j=js,je
        do i=is,ie+1
          i1 = id1_c(i,j)
          i2 = id2_c(i,j)
          j1 = jdc_c(i,j)
          ut_inc(i,j) = s2c_c(i,j,1)*u_amb(i1,j1  ,k) + &
                        s2c_c(i,j,2)*u_amb(i2,j1  ,k) + &
                        s2c_c(i,j,3)*u_amb(i2,j1+1,k) + &
                        s2c_c(i,j,4)*u_amb(i1,j1+1,k)
          vt_inc(i,j) = s2c_c(i,j,1)*v_amb(i1,j1  ,k) + &
                        s2c_c(i,j,2)*v_amb(i2,j1  ,k) + &
                        s2c_c(i,j,3)*v_amb(i2,j1+1,k) + &
                        s2c_c(i,j,4)*v_amb(i1,j1+1,k)
          p1(:) = atm(1)%gridstruct%grid(i,j  ,1:2)
          p2(:) = atm(1)%gridstruct%grid(i,j+1,1:2)
          call  mid_pt_sphere(p1, p2, p3)
          call get_unit_vect2(p1, p2, e2)
          call get_latlon_vector(p3, ex, ey)
          vd_inc(i,j,k) = ut_inc(i,j)*inner_prod(e2,ex) + &
                          vt_inc(i,j)*inner_prod(e2,ey)
        enddo
      enddo
    enddo

    deallocate ( ut_inc, vt_inc )

    !------ pt_d part ------
    ! Initialize lat-lon to Cubed bi-linear interpolation coeff:
    call remap_coef( is, ie, js, je+1, isd, ied, jsd, jed+1, &
        im, jm, lon, lat, id1_d, id2_d, jdc_d, s2c_d, &
        pt_d)

    allocate (  ut_inc(is:ie,js:je+1) )
    allocate (  vt_inc(is:ie,js:je+1) )

    do k=1,km
      do j=js,je+1
        do i=is,ie
          i1 = id1_d(i,j)
          i2 = id2_d(i,j)
          j1 = jdc_d(i,j)
          ut_inc(i,j) = s2c_d(i,j,1)*u_amb(i1,j1  ,k) + &
                        s2c_d(i,j,2)*u_amb(i2,j1  ,k) + &
                        s2c_d(i,j,3)*u_amb(i2,j1+1,k) + &
                        s2c_d(i,j,4)*u_amb(i1,j1+1,k)
          vt_inc(i,j) = s2c_d(i,j,1)*v_amb(i1,j1  ,k) + &
                        s2c_d(i,j,2)*v_amb(i2,j1  ,k) + &
                        s2c_d(i,j,3)*v_amb(i2,j1+1,k) + &
                        s2c_d(i,j,4)*v_amb(i1,j1+1,k)
          p1(:) = atm(1)%gridstruct%grid(i,  j,1:2)
          p2(:) = atm(1)%gridstruct%grid(i+1,j,1:2)
          call  mid_pt_sphere(p1, p2, p3)
          call get_unit_vect2(p1, p2, e1)
          call get_latlon_vector(p3, ex, ey)
          ud_inc(i,j,k) = ut_inc(i,j)*inner_prod(e1,ex) + &
                          vt_inc(i,j)*inner_prod(e1,ey)
        enddo
      enddo
    enddo

    deallocate ( ut_inc, vt_inc )

#if 0
    allocate (ua_inc(isd:ied  , jsd:jed  , km))
    allocate (va_inc(isd:ied  , jsd:jed  , km))
    call cubed_to_latlon(ud_inc, vd_inc, ua_inc, va_inc, atm(1)%gridstruct,   &
          atm(1)%npx, atm(1)%npy, atm(1)%npz, update_uv_ghost_cells, atm(1)%gridstruct%grid_type, &
          fv_domain, atm(1)%gridstruct%nested, atm(1)%flagstruct%c2l_ord, atm(1)%bd)
    u_inc(is:ie,js:je,1:npz) = ua_inc(is:ie,js:je,1:npz)
    v_inc(is:ie,js:je,1:npz) = va_inc(is:ie,js:je,1:npz)
    deallocate ( ua_inc, va_inc )
#else
    u_inc(is:ie,js:je,1:npz) = ud_inc(is:ie,js:je,1:npz)
    v_inc(is:ie,js:je,1:npz) = vd_inc(is:ie,js:je,1:npz)
#endif

    !------ winds clean up ------
    deallocate ( pt_c, pt_d )
    deallocate ( ud_inc, vd_inc )

  contains
    !---------------------------------------------------------------------------
    subroutine get_inc_on_3d_scalar(amb,inc)
      real, dimension(1:im, 1:jm, 1:km), intent(in   ) :: amb
      real, dimension(is:ie,js:je,1:km), intent(  out) :: inc

      do k=1,km
        do j=js,je
          do i=is,ie
            i1 = id1(i,j)
            i2 = id2(i,j)
            j1 = jdc(i,j)
            inc(i,j,k) = s2c(i,j,1)*amb(i1,j1  ,k) + &
                         s2c(i,j,2)*amb(i2,j1  ,k) + &
                         s2c(i,j,3)*amb(i2,j1+1,k) + &
                         s2c(i,j,4)*amb(i1,j1+1,k)
          enddo
        enddo
      enddo

    end subroutine get_inc_on_3d_scalar
    !---------------------------------------------------------------------------
  end subroutine geos_get_da_increments

#ifndef MAPL_MODE
  !=============================================================================
  !> @brief description
  !> @author Xi.Chen <xi.chen@noaa.gov>
  !> @date 02/12/2016
  subroutine read_da_inc(Atm, fv_domain)
    type(fv_atmos_type), intent(inout) :: atm(:)
    type(domain2d),      intent(inout) :: fv_domain
    ! local
    integer :: nq

    real :: deg2rad
    character(len=128) :: fname
    real(kind=4), allocatable:: wk1(:), wk2(:,:), wk3(:,:,:)
    real(kind=4), allocatable:: wk3_u(:,:,:), wk3_v(:,:,:)
    real, allocatable:: tp(:,:,:), qp(:,:,:)
    real, dimension(:,:,:), allocatable:: u_inc, v_inc, ud_inc, vd_inc
    real, allocatable:: lat(:), lon(:)
    real, allocatable:: pt_c(:,:,:), pt_d(:,:,:)
    real:: s2c(atm(1)%bd%is:atm(1)%bd%ie,atm(1)%bd%js:atm(1)%bd%je,4)
    real:: s2c_c(atm(1)%bd%is:atm(1)%bd%ie+1,atm(1)%bd%js:atm(1)%bd%je,4)
    real:: s2c_d(atm(1)%bd%is:atm(1)%bd%ie,atm(1)%bd%js:atm(1)%bd%je+1,4)
    integer, dimension(Atm(1)%bd%is:Atm(1)%bd%ie,Atm(1)%bd%js:Atm(1)%bd%je):: &
        id1, id2, jdc
    integer, dimension(Atm(1)%bd%is:Atm(1)%bd%ie+1,Atm(1)%bd%js:Atm(1)%bd%je)::&
        id1_c, id2_c, jdc_c
    integer, dimension(Atm(1)%bd%is:Atm(1)%bd%ie,Atm(1)%bd%js:Atm(1)%bd%je+1)::&
        id1_d, id2_d, jdc_d

    integer:: i, j, k, im, jm, km, npz, npt
    integer:: i1, i2, j1, ncid
    integer:: jbeg, jend
    integer tsize(3)
    real(kind=R_GRID), dimension(2):: p1, p2, p3
    real(kind=R_GRID), dimension(3):: e1, e2, ex, ey

    logical:: found
    integer :: is,  ie,  js,  je
    integer :: isd, ied, jsd, jed
    integer :: sphum, liq_wat, o3mr

    is  = atm(1)%bd%is
    ie  = atm(1)%bd%ie
    js  = atm(1)%bd%js
    je  = atm(1)%bd%je
    isd = atm(1)%bd%isd
    ied = atm(1)%bd%ied
    jsd = atm(1)%bd%jsd
    jed = atm(1)%bd%jed

    deg2rad = pi/180.

    npz = atm(1)%npz

    fname = 'INPUT/'//atm(1)%flagstruct%res_latlon_dynamics

    if( file_exist(fname) ) then
      call open_ncfile( fname, ncid )        ! open the file
      call get_ncdim1( ncid, 'lon',   tsize(1) )
      call get_ncdim1( ncid, 'lat',   tsize(2) )
      call get_ncdim1( ncid, 'lev', tsize(3) )

      im = tsize(1); jm = tsize(2); km = tsize(3)

      if (km.ne.npz) then
        if (is_master()) print *, 'km = ', km
        call mpp_error(fatal, &
            '==> Error in read_da_inc: km is not equal to npz')
      endif

      if(is_master())  write(*,*) fname, ' DA increment dimensions:', tsize

      allocate (  lon(im) )
      allocate (  lat(jm) )

      call _get_var1 (ncid, 'lon', im, lon )
      call _get_var1 (ncid, 'lat', jm, lat )

      ! Convert to radian
      do i=1,im
        lon(i) = lon(i) * deg2rad  ! lon(1) = 0.
      enddo
      do j=1,jm
        lat(j) = lat(j) * deg2rad
      enddo

    else
      call mpp_error(fatal,'==> Error in read_da_inc: Expected file '&
          //trim(fname)//' for DA increment does not exist')
    endif

    ! Initialize lat-lon to Cubed bi-linear interpolation coeff:
    call remap_coef( is, ie, js, je, isd, ied, jsd, jed, &
        im, jm, lon, lat, id1, id2, jdc, s2c, &
        atm(1)%gridstruct%agrid)

    ! Find bounding latitudes:
    jbeg = jm-1;         jend = 2
    do j=js,je
      do i=is,ie
          j1 = jdc(i,j)
        jbeg = min(jbeg, j1) 
        jend = max(jend, j1+1)
      enddo
    enddo

    sphum   = get_tracer_index(model_atmos, 'sphum')
    o3mr    = get_tracer_index(model_atmos, 'o3mr')
    liq_wat = get_tracer_index(model_atmos, 'liq_wat')

    ! perform increments on scalars
    allocate ( wk3(1:im,jbeg:jend, 1:km) )
    allocate (  tp(is:ie,js:je,km) )

    call apply_inc_on_3d_scalar('T_inc',atm(1)%pt)
    call apply_inc_on_3d_scalar('delp_inc',atm(1)%delp)
    call apply_inc_on_3d_scalar('sphum_inc',atm(1)%q(:,:,:,sphum))
    call apply_inc_on_3d_scalar('liq_wat_inc',atm(1)%q(:,:,:,liq_wat))
    call apply_inc_on_3d_scalar('o3mr_inc',atm(1)%q(:,:,:,o3mr))

    deallocate ( tp )
    deallocate ( wk3 )

    ! perform increments on winds
    allocate (pt_c(isd:ied+1,jsd:jed  ,2))
    allocate (pt_d(isd:ied  ,jsd:jed+1,2))
    allocate (ud_inc(is:ie  , js:je+1, km))
    allocate (vd_inc(is:ie+1, js:je  , km))

    call get_staggered_grid( &
        is, ie, js, je, &
        isd, ied, jsd, jed, &
        atm(1)%gridstruct%grid, pt_c, pt_d)

    !------ pt_c part ------
    ! Initialize lat-lon to Cubed bi-linear interpolation coeff:
    call remap_coef( is, ie+1, js, je, isd, ied+1, jsd, jed, &
        im, jm, lon, lat, id1_c, id2_c, jdc_c, s2c_c, &
        pt_c)

    ! Find bounding latitudes:
    jbeg = jm-1;         jend = 2
    do j=js,je
      do i=is,ie+1
          j1 = jdc_c(i,j)
        jbeg = min(jbeg, j1) 
        jend = max(jend, j1+1)
      enddo
    enddo

    allocate ( wk3_u(1:im,jbeg:jend, 1:km) )
    allocate ( wk3_v(1:im,jbeg:jend, 1:km) )
    allocate (  u_inc(is:ie+1,js:je,km) )
    allocate (  v_inc(is:ie+1,js:je,km) )

    call get_var3_r4( ncid, 'u_inc', 1,im, jbeg,jend, 1,km, wk3_u )
    call get_var3_r4( ncid, 'v_inc', 1,im, jbeg,jend, 1,km, wk3_v )

    do k=1,km
      do j=js,je
        do i=is,ie+1
          i1 = id1_c(i,j)
          i2 = id2_c(i,j)
          j1 = jdc_c(i,j)
          u_inc(i,j,k) = s2c_c(i,j,1)*wk3_u(i1,j1  ,k) + &
                         s2c_c(i,j,2)*wk3_u(i2,j1  ,k) + &
                         s2c_c(i,j,3)*wk3_u(i2,j1+1,k) + &
                         s2c_c(i,j,4)*wk3_u(i1,j1+1,k)
          v_inc(i,j,k) = s2c_c(i,j,1)*wk3_v(i1,j1  ,k) + &
                         s2c_c(i,j,2)*wk3_v(i2,j1  ,k) + &
                         s2c_c(i,j,3)*wk3_v(i2,j1+1,k) + &
                         s2c_c(i,j,4)*wk3_v(i1,j1+1,k)
          p1(:) = atm(1)%gridstruct%grid(i,j  ,1:2)
          p2(:) = atm(1)%gridstruct%grid(i,j+1,1:2)
          call  mid_pt_sphere(p1, p2, p3)
          call get_unit_vect2(p1, p2, e2)
          call get_latlon_vector(p3, ex, ey)
          vd_inc(i,j,k) = u_inc(i,j,k)*inner_prod(e2,ex) + &
                          v_inc(i,j,k)*inner_prod(e2,ey)
          atm(1)%v(i,j,k) = atm(1)%v(i,j,k) + vd_inc(i,j,k)
        enddo
      enddo
    enddo

    deallocate ( u_inc, v_inc )
    deallocate ( wk3_u, wk3_v )

    !------ pt_d part ------
    ! Initialize lat-lon to Cubed bi-linear interpolation coeff:
    call remap_coef( is, ie, js, je+1, isd, ied, jsd, jed+1, &
        im, jm, lon, lat, id1_d, id2_d, jdc_d, s2c_d, &
        pt_d)

    ! Find bounding latitudes:
    jbeg = jm-1;         jend = 2
    do j=js,je+1
      do i=is,ie
          j1 = jdc_d(i,j)
        jbeg = min(jbeg, j1) 
        jend = max(jend, j1+1)
      enddo
    enddo

    allocate ( wk3_u(1:im,jbeg:jend, 1:km) )
    allocate ( wk3_v(1:im,jbeg:jend, 1:km) )
    allocate (  u_inc(is:ie,js:je+1,km) )
    allocate (  v_inc(is:ie,js:je+1,km) )

    call get_var3_r4( ncid, 'u_inc', 1,im, jbeg,jend, 1,km, wk3_u )
    call get_var3_r4( ncid, 'v_inc', 1,im, jbeg,jend, 1,km, wk3_v )

    do k=1,km
      do j=js,je+1
        do i=is,ie
          i1 = id1_d(i,j)
          i2 = id2_d(i,j)
          j1 = jdc_d(i,j)
          u_inc(i,j,k) = s2c_d(i,j,1)*wk3_u(i1,j1  ,k) + &
                         s2c_d(i,j,2)*wk3_u(i2,j1  ,k) + &
                         s2c_d(i,j,3)*wk3_u(i2,j1+1,k) + &
                         s2c_d(i,j,4)*wk3_u(i1,j1+1,k)
          v_inc(i,j,k) = s2c_d(i,j,1)*wk3_v(i1,j1  ,k) + &
                         s2c_d(i,j,2)*wk3_v(i2,j1  ,k) + &
                         s2c_d(i,j,3)*wk3_v(i2,j1+1,k) + &
                         s2c_d(i,j,4)*wk3_v(i1,j1+1,k)
          p1(:) = atm(1)%gridstruct%grid(i,  j,1:2)
          p2(:) = atm(1)%gridstruct%grid(i+1,j,1:2)
          call  mid_pt_sphere(p1, p2, p3)
          call get_unit_vect2(p1, p2, e1)
          call get_latlon_vector(p3, ex, ey)
          ud_inc(i,j,k) = u_inc(i,j,k)*inner_prod(e1,ex) + &
                          v_inc(i,j,k)*inner_prod(e1,ey)
          atm(1)%u(i,j,k) = atm(1)%u(i,j,k) + ud_inc(i,j,k)
        enddo
      enddo
    enddo

    deallocate ( u_inc, v_inc )
    deallocate ( wk3_u, wk3_v )

!rab The following is not necessary as ua/va will be re-calculated during model startup
!rab    call cubed_to_latlon(Atm(1)%u, Atm(1)%v, Atm(1)%ua, Atm(1)%va, &
!rab        Atm(1)%gridstruct, Atm(1)%flagstruct%npx, Atm(1)%flagstruct%npy, &
!rab        Atm(1)%flagstruct%npz, 1, Atm(1)%gridstruct%grid_type, &
!rab        fv_domain, Atm(1)%gridstruct%nested, &
!rab        Atm(1)%flagstruct%c2l_ord, Atm(1)%bd)

    !------ winds clean up ------
    deallocate ( pt_c, pt_d, ud_inc, vd_inc )
    !------ all clean up ------
    deallocate ( lon, lat )

  contains
    !---------------------------------------------------------------------------
    subroutine apply_inc_on_3d_scalar(field_name,var)
      character(len=*), intent(in) :: field_name
      real, dimension(isd:ied,jsd:jed,1:km), intent(inout) :: var

      call get_var3_r4( ncid, field_name, 1,im, jbeg,jend, 1,km, wk3 )

      do k=1,km
        do j=js,je
          do i=is,ie
            i1 = id1(i,j)
            i2 = id2(i,j)
            j1 = jdc(i,j)
            tp(i,j,k) = s2c(i,j,1)*wk3(i1,j1  ,k) + s2c(i,j,2)*wk3(i2,j1  ,k)+&
                        s2c(i,j,3)*wk3(i2,j1+1,k) + s2c(i,j,4)*wk3(i1,j1+1,k)
            var(i,j,k) = var(i,j,k)+tp(i,j,k)
          enddo
        enddo
      enddo

    end subroutine apply_inc_on_3d_scalar
    !---------------------------------------------------------------------------
  end subroutine read_da_inc
#endif

  !=============================================================================
  subroutine remap_coef( is, ie, js, je, isd, ied, jsd, jed, &
      im, jm, lon, lat, id1, id2, jdc, s2c, agrid )

    integer, intent(in):: is, ie, js, je, isd, ied, jsd, jed
    integer, intent(in):: im, jm
    real,    intent(in):: lon(im), lat(jm)
    real,    intent(out):: s2c(is:ie,js:je,4)
    integer, intent(out), dimension(is:ie,js:je):: id1, id2, jdc
    real,    intent(in):: agrid(isd:ied,jsd:jed,2)
    ! local:
    real :: rdlon(im)
    real :: rdlat(jm)
    real:: a1, b1
    integer i,j, i1, i2, jc, i0, j0
  
    do i=1,im-1
      rdlon(i) = 1. / (lon(i+1) - lon(i))
    enddo
    rdlon(im) = 1. / (lon(1) + 2.*pi - lon(im))

    do j=1,jm-1
      rdlat(j) = 1. / (lat(j+1) - lat(j))
    enddo

    ! * Interpolate to cubed sphere cell center
    do 5000 j=js,je

      do i=is,ie

        if ( agrid(i,j,1)>lon(im) ) then
          i1 = im;     i2 = 1
          a1 = (agrid(i,j,1)-lon(im)) * rdlon(im)
        elseif ( agrid(i,j,1)<lon(1) ) then
          i1 = im;     i2 = 1
          a1 = (agrid(i,j,1)+2.*pi-lon(im)) * rdlon(im)
        else
          do i0=1,im-1
            if ( agrid(i,j,1)>=lon(i0) .and. agrid(i,j,1)<=lon(i0+1) ) then
              i1 = i0;  i2 = i0+1
              a1 = (agrid(i,j,1)-lon(i1)) * rdlon(i0)
              go to 111
            endif
          enddo
        endif
111     continue

        if ( agrid(i,j,2)<lat(1) ) then
          jc = 1
          b1 = 0.
        elseif ( agrid(i,j,2)>lat(jm) ) then
          jc = jm-1
          b1 = 1.
        else
          do j0=1,jm-1
            if ( agrid(i,j,2)>=lat(j0) .and. agrid(i,j,2)<=lat(j0+1) ) then
              jc = j0
              b1 = (agrid(i,j,2)-lat(jc)) * rdlat(jc)
              go to 222
            endif
          enddo
        endif
222     continue

        if ( a1<0.0 .or. a1>1.0 .or.  b1<0.0 .or. b1>1.0 ) then
             write(*,*) 'gid=', mpp_pe(), i,j,a1, b1
        endif

        s2c(i,j,1) = (1.-a1) * (1.-b1)
        s2c(i,j,2) =     a1  * (1.-b1)
        s2c(i,j,3) =     a1  *     b1
        s2c(i,j,4) = (1.-a1) *     b1
        id1(i,j) = i1
        id2(i,j) = i2
        jdc(i,j) = jc
      enddo   !i-loop
5000 continue   ! j-loop

  end subroutine remap_coef
  !=============================================================================
  subroutine get_staggered_grid( &
      is, ie, js, je, &
      isd, ied, jsd, jed, &
      pt_b, pt_c, pt_d)
    integer, intent(in) :: is, ie, js, je, isd, ied, jsd, jed
    real, dimension(isd:ied+1,jsd:jed+1,2), intent(in) :: pt_b
    real, dimension(isd:ied+1,jsd:jed  ,2), intent(out) :: pt_c
    real, dimension(isd:ied  ,jsd:jed+1,2), intent(out) :: pt_d
    ! local
    real(kind=R_GRID), dimension(2):: p1, p2, p3
    integer :: i, j

    do j = js,je+1
      do i = is,ie
        p1(:) = pt_b(i,  j,1:2)
        p2(:) = pt_b(i+1,j,1:2)
        call  mid_pt_sphere(p1, p2, p3)
        pt_d(i,j,1:2) = p3(:)
      enddo
    enddo
    do j = js,je
      do i = is,ie+1
        p1(:) = pt_b(i,j  ,1:2)
        p2(:) = pt_b(i,j+1,1:2)
        call  mid_pt_sphere(p1, p2, p3)
        pt_c(i,j,1:2) = p3(:)
      enddo
    enddo

  end subroutine get_staggered_grid
  !=============================================================================
end module fv_treat_da_inc_nlm_mod