Add ability to read increment files on native cubed sphere grid (#342)

* Initial commit

* Polishing up

* Debugging

* multiple updates

* Clean up

* Final updates (I hope)

* Spelling fix

* Deleted accidentally undeleted files

* Fix bug in GCC compilation related to float conversion

* Initial commit

* revert sim_nc_mod

* Implement FMS increment reads

* Fix filename issues

* Make sure to register axes

* Update x and y-axis names for increment to match restart axis names

* Last minute update of y-axis name

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
@@ -950,7 +951,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
@@ -1093,10 +1094,9 @@ 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, &
+            write_coarse_agrid_vel_rst, write_coarse_dgrid_vel_rst, increment_file_on_native_grid, &
             pass_full_omega_to_physics_in_non_hydrostatic_mode, ignore_rst_cksum
 
-
        ! Read FVCORE namelist
        read (input_nml_file,fv_core_nml,iostat=ios)
        ierr = check_nml_error(ios,'fv_core_nml')

tools/cubed_sphere_inc_mod.F90

@@ -0,0 +1,66 @@
+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 fms2_io_mod,        only: open_file, close_file, read_data, variable_exists, &
+                                FmsNetcdfDomainFile_t, register_axis
+
+  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(fname, increment_data, Atm)
+      character(*), intent(in)                 :: fname
+      type(increment_data_type), intent(inout) :: increment_data
+      type(fv_atmos_type), intent(in)          :: Atm
+
+      type(FmsNetcdfDomainFile_t) :: fileobj
+      integer                     :: itracer, ntracers
+      character(len=64)           :: tracer_name
+
+      ! Get various dimensions
+      call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers)
+
+      ! Open file
+      if ( open_file(fileobj, trim(fname), 'read', Atm%domain) ) then
+         ! Register axes
+         call register_axis(fileobj, 'xaxis_1', 'x')
+         call register_axis(fileobj, 'yaxis_1', 'y')
+
+         ! Read increments
+         call read_data(fileobj, 'u_inc', increment_data%ua_inc)
+         call read_data(fileobj, 'v_inc', increment_data%va_inc)
+         call read_data(fileobj, 'T_inc', increment_data%temp_inc)
+         call read_data(fileobj, 'delp_inc', increment_data%delp_inc)
+         if ( .not. Atm%flagstruct%hydrostatic ) then
+            call read_data(fileobj, 'delz_inc', increment_data%delz_inc)
+         end if
+
+         ! Read tracer increments
+         do itracer = 1,ntracers
+            call get_tracer_names(MODEL_ATMOS, itracer, tracer_name)
+            if ( variable_exists(fileobj, trim(tracer_name)//'_inc') ) then
+               call read_data(fileobj, trim(tracer_name)//'_inc', increment_data%tracer_inc(:,:,:,itracer))
+            end if
+         end do
+
+         call close_file(fileobj)
+      end if
+
+  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
@@ -471,11 +471,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 !n==this_grid