Add ability to read increment files on native cubed sphere grid

CMakeLists.txt

@@ -102,7 +102,8 @@ list(APPEND tools_srcs
   tools/module_diag_hailcast.F90
   tools/sim_nc_mod.F90
   tools/sorted_index.F90
-  tools/test_cases.F90)
+  tools/test_cases.F90
+  tools/cubed_sphere_inc_mod.F90)
 
 list(APPEND tools_srcs_extra
   tools/fv_iau_mod.F90)

model/fv_arrays.F90

@@ -831,6 +831,7 @@ module fv_arrays_mod
                                        !< condition file (if nggps_ic or ecwmf_ic are .true.). This overrides the
                                        !< hard-coded levels in fv_eta. The default is .false.
    logical :: read_increment = .false.   !< read in analysis increment and add to restart
+   logical :: increment_file_on_native_grid = .false. !< increment is on native cubed sphere grid grid else on Gaussian grid
 ! following are namelist parameters for Stochastic Energy Baskscatter dissipation estimate
    logical :: do_skeb  = .false.         !< save dissipation estimate
    integer :: skeb_npass  = 11           !< Filter dissipation estimate "skeb_npass" times

model/fv_control.F90

@@ -373,6 +373,7 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split,
      logical , pointer :: external_ic
      logical , pointer :: external_eta
      logical , pointer :: read_increment
+     logical , pointer :: increment_file_on_native_grid
      logical , pointer :: hydrostatic
      logical , pointer :: phys_hydrostatic
      logical , pointer :: use_hydro_pressure
@@ -949,7 +950,7 @@ subroutine set_namelist_pointers(Atm)
        external_ic                   => Atm%flagstruct%external_ic
        external_eta                  => Atm%flagstruct%external_eta
        read_increment                => Atm%flagstruct%read_increment
-
+       increment_file_on_native_grid => Atm%flagstruct%increment_file_on_native_grid
        hydrostatic                   => Atm%flagstruct%hydrostatic
        phys_hydrostatic              => Atm%flagstruct%phys_hydrostatic
        use_hydro_pressure            => Atm%flagstruct%use_hydro_pressure
@@ -1091,7 +1092,7 @@ subroutine read_namelist_fv_core_nml(Atm)
             write_coarse_restart_files,&
             write_coarse_diagnostics,&
             write_only_coarse_intermediate_restarts, &
-            write_coarse_agrid_vel_rst, write_coarse_dgrid_vel_rst, ignore_rst_cksum
+            write_coarse_agrid_vel_rst, write_coarse_dgrid_vel_rst, ignore_rst_cksum, increment_file_on_native_grid
 
 
        ! Read FVCORE namelist

tools/cubed_sphere_inc_mod.F90

@@ -0,0 +1,69 @@
+module cubed_sphere_inc_mod
+
+  use tracer_manager_mod, only: get_tracer_index, get_number_tracers, get_tracer_names
+  use field_manager_mod,  only: MODEL_ATMOS
+  use fv_arrays_mod,      only: fv_atmos_type
+  use sim_nc_mod,         only: open_ncfile, get_var5_r4, check_var_exists, close_ncfile
+
+  implicit none
+  type increment_data_type
+    real, allocatable :: ua_inc(:,:,:)
+    real, allocatable :: va_inc(:,:,:)
+    real, allocatable :: temp_inc(:,:,:)
+    real, allocatable :: delp_inc(:,:,:)
+    real, allocatable :: delz_inc(:,:,:)
+    real, allocatable :: tracer_inc(:,:,:,:)
+  end type increment_data_type
+
+  public read_cubed_sphere_inc, increment_data_type
+
+  contains
+
+!----------------------------------------------------------------------------------------
+
+    subroutine read_cubed_sphere_inc(filename, increment_data, Atm)
+      character(*), intent(in)                 :: filename
+      type(increment_data_type), intent(inout) :: increment_data
+      type(fv_atmos_type), intent(in)          :: Atm
+
+      integer           :: isc, iec, jsc, jec, npz, itracer, ntracers, itile
+      character(len=64) :: tracer_name
+      integer           :: ncid, ncerr
+
+      ! Get various dimensions
+      call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers)
+      isc = Atm%bd%isc
+      iec = Atm%bd%iec
+      jsc = Atm%bd%jsc
+      jec = Atm%bd%jec
+      npz = Atm%npz
+      itile = Atm%tile_of_mosaic
+
+      ! Open increment file
+      call open_ncfile( trim(filename), ncid )
+
+      ! Read increments
+      call get_var5_r4( ncid, 'u_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%ua_inc )
+      call get_var5_r4( ncid, 'v_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%va_inc )
+      call get_var5_r4( ncid, 'T_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%temp_inc )
+      call get_var5_r4( ncid, 'delp_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%delp_inc )
+      if ( .not. Atm%flagstruct%hydrostatic ) then
+         call get_var5_r4( ncid, 'delz_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%delz_inc )
+      end if
+
+      ! Read tracer increments
+      do itracer = 1,ntracers
+         call get_tracer_names(MODEL_ATMOS, itracer, tracer_name)
+         call check_var_exists(ncid, trim(tracer_name)//'_inc', ncerr)
+         if ( ncerr == 0 ) then
+            call get_var5_r4( ncid, trim(tracer_name)//'_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%tracer_inc(:,:,:,itracer) )
+         end if
+      end do
+
+      ! Close increment file
+      call close_ncfile(ncid)
+
+  end subroutine read_cubed_sphere_inc
+
+!----------------------------------------------------------------------------------------
+end module cubed_sphere_inc_mod

tools/fv_iau_mod.F90

@@ -66,6 +66,7 @@ module fv_iau_mod
   use fv_treat_da_inc_mod, only: remap_coef
   use tracer_manager_mod,  only: get_tracer_names,get_tracer_index, get_number_tracers
   use field_manager_mod,   only: MODEL_ATMOS
+  use cubed_sphere_inc_mod, only : read_cubed_sphere_inc, increment_data_type
   implicit none
 
   private
@@ -84,14 +85,6 @@ module fv_iau_mod
   integer, allocatable :: tracer_indicies(:)
 
   real(kind=4), allocatable:: wk3(:,:,:)
-  type iau_internal_data_type
-    real,allocatable :: ua_inc(:,:,:)
-    real,allocatable :: va_inc(:,:,:)
-    real,allocatable :: temp_inc(:,:,:)
-    real,allocatable :: delp_inc(:,:,:)
-    real,allocatable :: delz_inc(:,:,:)
-    real,allocatable :: tracer_inc(:,:,:,:)
-  end type iau_internal_data_type
   type iau_external_data_type
     real,allocatable :: ua_inc(:,:,:)
     real,allocatable :: va_inc(:,:,:)
@@ -103,8 +96,8 @@ module fv_iau_mod
     logical          :: drymassfixer = .false.
   end type iau_external_data_type
   type iau_state_type
-      type(iau_internal_data_type):: inc1
-      type(iau_internal_data_type):: inc2
+      type(increment_data_type)   :: inc1
+      type(increment_data_type)   :: inc2
       real(kind=kind_phys)        :: hr1
       real(kind=kind_phys)        :: hr2
       real(kind=kind_phys)        :: wt
@@ -114,10 +107,11 @@ module fv_iau_mod
   public iau_external_data_type,IAU_initialize,getiauforcing
 
 contains
-subroutine IAU_initialize (IPD_Control, IAU_Data,Init_parm)
+subroutine IAU_initialize (IPD_Control, IAU_Data, Init_parm, Atm)
     type (IPD_control_type), intent(in) :: IPD_Control
     type (IAU_external_data_type), intent(inout) :: IAU_Data
     type (IPD_init_type),    intent(in) :: Init_parm
+    type (fv_atmos_type), intent(inout) :: Atm
     ! local
 
     character(len=128) :: fname
@@ -274,7 +268,11 @@ subroutine IAU_initialize (IPD_Control, IAU_Data,Init_parm)
        enddo
        iau_state%wt_normfact = (2*nstep+1)/normfact
     endif
-    call read_iau_forcing(IPD_Control,iau_state%inc1,'INPUT/'//trim(IPD_Control%iau_inc_files(1)))
+    if ( Atm%flagstruct%increment_file_on_native_grid ) then
+       call read_cubed_sphere_inc('INPUT/'//trim(IPD_Control%iau_inc_files(1)), iau_state%inc1, Atm)
+    else
+       call read_iau_forcing(IPD_Control,iau_state%inc1,'INPUT/'//trim(IPD_Control%iau_inc_files(1)))
+    endif
     if (nfiles.EQ.1) then  ! only need to get incrments once since constant forcing over window
        call setiauforcing(IPD_Control,IAU_Data,iau_state%wt)
     endif
@@ -286,18 +284,23 @@ subroutine IAU_initialize (IPD_Control, IAU_Data,Init_parm)
        allocate (iau_state%inc2%delz_inc (is:ie, js:je, km))
        allocate (iau_state%inc2%tracer_inc(is:ie, js:je, km,ntracers))
        iau_state%hr2=IPD_Control%iaufhrs(2)
-       call read_iau_forcing(IPD_Control,iau_state%inc2,'INPUT/'//trim(IPD_Control%iau_inc_files(2)))
+       if ( Atm%flagstruct%increment_file_on_native_grid ) then
+          call read_cubed_sphere_inc('INPUT/'//trim(IPD_Control%iau_inc_files(2)), iau_state%inc2, Atm)
+       else
+          call read_iau_forcing(IPD_Control,iau_state%inc2,'INPUT/'//trim(IPD_Control%iau_inc_files(2)))
+       endif
     endif
 !   print*,'in IAU init',dt,rdt
     IAU_data%drymassfixer = IPD_control%iau_drymassfixer
 
 end subroutine IAU_initialize
 
-subroutine getiauforcing(IPD_Control,IAU_Data)
+subroutine getiauforcing(IPD_Control,IAU_Data,Atm)
 
    implicit none
    type (IPD_control_type), intent(in) :: IPD_Control
    type(IAU_external_data_type),  intent(inout) :: IAU_Data
+   type (fv_atmos_type), intent(inout) :: Atm
    real(kind=kind_phys) t1,t2,sx,wx,wt,dtp
    integer n,i,j,k,sphum,kstep,nstep,itnext
 
@@ -371,7 +374,11 @@ subroutine getiauforcing(IPD_Control,IAU_Data)
             iau_state%hr2=IPD_Control%iaufhrs(itnext)
             iau_state%inc1=iau_state%inc2
             if (is_master()) print *,'reading next increment file',trim(IPD_Control%iau_inc_files(itnext))
-            call read_iau_forcing(IPD_Control,iau_state%inc2,'INPUT/'//trim(IPD_Control%iau_inc_files(itnext)))
+            if ( Atm%flagstruct%increment_file_on_native_grid ) then
+               call read_cubed_sphere_inc('INPUT/'//trim(IPD_Control%iau_inc_files(itnext)), iau_state%inc2, Atm)
+            else
+               call read_iau_forcing(IPD_Control,iau_state%inc2,'INPUT/'//trim(IPD_Control%iau_inc_files(itnext)))
+            endif
          endif
          call updateiauforcing(IPD_Control,IAU_Data,iau_state%wt)
       endif
@@ -434,7 +441,7 @@ end subroutine setiauforcing
 
 subroutine read_iau_forcing(IPD_Control,increments,fname)
     type (IPD_control_type), intent(in) :: IPD_Control
-    type(iau_internal_data_type), intent(inout):: increments
+    type(increment_data_type), intent(inout):: increments
     character(len=*),  intent(in) :: fname
 !locals
     real, dimension(:,:,:), allocatable:: u_inc, v_inc

tools/fv_restart.F90

@@ -88,7 +88,7 @@ module fv_restart_mod
 !   </tr>
 !   <tr>
 !     <td>fv_treat_da_inc_mod</td>
-!     <td>read_da_inc</td>
+!     <td>read_da_inc, read_da_inc_cubed_sphere</td>
 !   </tr>
 !   <tr>
 !     <td>fv_timing_mod</td>
@@ -166,7 +166,7 @@ module fv_restart_mod
   use mpp_mod,             only: mpp_send, mpp_recv, mpp_sync_self, mpp_set_current_pelist, mpp_get_current_pelist, mpp_npes, mpp_pe, mpp_sync
   use mpp_domains_mod,     only: CENTER, CORNER, NORTH, EAST,  mpp_get_C2F_index, WEST, SOUTH
   use mpp_domains_mod,     only: mpp_global_field
-  use fv_treat_da_inc_mod, only: read_da_inc
+  use fv_treat_da_inc_mod, only: read_da_inc, read_da_inc_cubed_sphere
   use fms2_io_mod,         only: file_exists, set_filename_appendix, FmsNetcdfFile_t, open_file, close_file
   use coarse_grained_restart_files_mod, only: fv_io_write_restart_coarse
   use fv_regional_mod,     only: write_full_fields
@@ -371,11 +371,19 @@ subroutine fv_restart(fv_domain, Atm, seconds, days, cold_start, grid_type, this
                    i = (Atm(n)%bd%isc + Atm(n)%bd%iec)/2
                    j = (Atm(n)%bd%jsc + Atm(n)%bd%jec)/2
                    k = Atm(n)%npz/2
-                   if( is_master() ) write(*,*) 'Calling read_da_inc',Atm(n)%pt(i,j,k)
-                   call read_da_inc(Atm(n), Atm(n)%domain, Atm(n)%bd, Atm(n)%npz, Atm(n)%ncnst, &
-                        Atm(n)%u, Atm(n)%v, Atm(n)%q, Atm(n)%delp, Atm(n)%pt, Atm(n)%delz, isd, jsd, ied, jed, &
-                        isc, jsc, iec, jec )
-                   if( is_master() ) write(*,*) 'Back from read_da_inc',Atm(n)%pt(i,j,k)
+                   if ( Atm(n)%flagstruct%increment_file_on_native_grid ) then
+                      if( is_master() ) write(*,*) 'Calling read_da_inc_cubed_sphere',Atm(n)%pt(i,j,k)
+                      call read_da_inc_cubed_sphere(Atm(n), Atm(n)%domain, Atm(n)%bd, Atm(n)%npz, Atm(n)%ncnst, &
+                           Atm(n)%u, Atm(n)%v, Atm(n)%q, Atm(n)%delp, Atm(n)%pt, Atm(n)%delz, isd, jsd, ied, jed, &
+                           isc, jsc, iec, jec )
+                      if( is_master() ) write(*,*) 'Back from read_da_inc_cubed_sphere',Atm(n)%pt(i,j,k)
+                   else
+                      if( is_master() ) write(*,*) 'Calling read_da_inc',Atm(n)%pt(i,j,k)
+                      call read_da_inc(Atm(n), Atm(n)%domain, Atm(n)%bd, Atm(n)%npz, Atm(n)%ncnst, &
+                           Atm(n)%u, Atm(n)%v, Atm(n)%q, Atm(n)%delp, Atm(n)%pt, Atm(n)%delz, isd, jsd, ied, jed, &
+                           isc, jsc, iec, jec )
+                      if( is_master() ) write(*,*) 'Back from read_da_inc',Atm(n)%pt(i,j,k)
+                   endif
                 endif
                 !====== end PJP added DA functionailty======
              endif