lnd_iau_mod.F90

!***********************************************************************
!*   GNU Lesser General Public License
!*   <http://www.gnu.org/licenses/>.
!***********************************************************************
 
 
!! \section land_iau_mod
 
!-------------------------------------------------------------------------------
 
module land_iau_mod
  
  use machine,                  only: kind_phys, kind_dyn
  use netcdf
 
  implicit none
 
  private
 
  type land_iau_external_data_type
      real(kind=kind_phys),allocatable :: stc_inc(:,:,:)   
      real(kind=kind_phys),allocatable :: slc_inc(:,:,:)   
      logical                          :: in_interval = .false.
      real(kind=kind_phys)              :: hr1        
      real(kind=kind_phys)              :: hr2      
      real(kind=kind_phys)              :: wt
      real(kind=kind_phys)              :: wt_normfact
      real(kind=kind_phys)              :: rdt      
      integer                           :: itnext ! track the increment steps here
  end type land_iau_external_data_type
 
!!> \section arg_table_land_iau_state_type Argument Table
!! \htmlinclude land_iau_state_type.html
!!
  ! land_iau_state_type holds 'raw' (not interpolated) inrements, 
  ! read during land_iau_mod_init
  type land_iau_state_type
      real(kind=kind_phys),allocatable :: stc_inc(:,:,:,:)
      real(kind=kind_phys),allocatable :: slc_inc(:,:,:,:) 
  end type land_iau_state_type
 
 
!!!> \section arg_table_land_iau_control_type Argument Table
!! \htmlinclude land_iau_control_type.html
!!
  type land_iau_control_type      
      integer :: isc
      integer :: jsc
      integer :: nx
      integer :: ny
      integer :: tile_num
      integer :: nblks
      integer, allocatable :: blksz(:)        ! this could vary for the last block
      integer, allocatable :: blk_strt_indx(:)
 
      integer              :: lsoil
      integer              :: lsnow_lsm
      logical              :: do_land_iau
      real(kind=kind_phys) :: iau_delthrs     ! iau time interval (to scale increments) in hours
      character(len=240)   :: iau_inc_files(7)     ! list of increment files
      real(kind=kind_phys) :: iaufhrs(7)      ! forecast hours associated with increment files
      logical              :: iau_filter_increments   
      integer              :: lsoil_incr      ! soil layers (from top) updated by DA   
      logical              :: upd_stc
      logical              :: upd_slc
      logical              :: do_stcsmc_adjustment  !do moisture/temperature adjustment for consistency after increment add
      real(kind=kind_phys) :: min_t_increment
      real(kind=kind_phys) :: min_slc_increment
 
      integer              :: me
      integer              :: mpi_root
      character(len=64)    :: fn_nml
      real(kind=kind_phys) :: dtp             
      real(kind=kind_phys) :: fhour           
 
      integer              :: ntimes
 
  end type land_iau_control_type
 
  public land_iau_control_type, land_iau_external_data_type, land_iau_state_type, land_iau_mod_set_control, &
         land_iau_mod_init, land_iau_mod_getiauforcing, land_iau_mod_finalize, calculate_landinc_mask
 
contains
 
subroutine land_iau_mod_set_control(Land_IAU_Control,fn_nml,input_nml_file, me, mpi_root, &
                                   isc, jsc, nx, ny, tile_num, nblks, blksz, &
                                   lsoil, lsnow_lsm, dtp, fhour, errmsg, errflg)          
 
   type (land_iau_control_type), intent(inout) :: land_iau_control
   character(*), intent(in)                    :: fn_nml
   character(len=:), intent(in), dimension(:), pointer :: input_nml_file
   integer, intent(in)                        :: me, mpi_root
   integer, intent(in)                        :: isc, jsc, nx, ny, tile_num, nblks, lsoil, lsnow_lsm
   integer, dimension(:),          intent(in) :: blksz                 !(one:) !GFS_Control%blksz
   real(kind=kind_phys), intent(in)           :: dtp                   
   real(kind=kind_phys), intent(in)           :: fhour                 
   character(len=*),              intent(out) :: errmsg
   integer,                       intent(out) :: errflg
   
   integer                                    :: nb, ix
   integer                                    :: nlunit = 360          ! unit for namelist  !, intent(in)
   integer                                    :: ios
   logical                                    :: exists
   character(len=512)                         :: ioerrmsg
 
   character(len=4)                           :: iosstr
 
   logical               :: do_land_iau               = .false.
   real(kind=kind_phys)  :: land_iau_delthrs              = 0           
   character(len=240)    :: land_iau_inc_files(7)         = ''
   real(kind=kind_phys)  :: land_iau_fhrs(7)               = -1         
   logical               :: land_iau_filter_increments    = .false.     
 
   integer               :: lsoil_incr = 4
   logical               :: land_iau_upd_stc = .false.
   logical               :: land_iau_upd_slc = .false.
   logical               :: land_iau_do_stcsmc_adjustment = .false.
   real(kind=kind_phys)  :: land_iau_min_t_increment = 0.0001
   real(kind=kind_phys)  :: land_iau_min_slc_increment = 0.000001
 
   namelist /land_iau_nml/ do_land_iau, land_iau_delthrs, land_iau_inc_files, land_iau_fhrs,   &  
                        land_iau_filter_increments, lsoil_incr, land_iau_upd_stc, land_iau_upd_slc, &
                        land_iau_do_stcsmc_adjustment, land_iau_min_t_increment, land_iau_min_slc_increment                                    
   !Errors messages handled through CCPP error handling variables
   errmsg = ''
   errflg = 0
 
!3.11.24: copied from GFS_typedefs.F90 
#ifdef INTERNAL_FILE_NML
    read(input_nml_file, nml=land_iau_nml, err=888, END=999, iostat=ios)
#else
   inquire (file=trim(fn_nml), exist=exists)    ! TODO: this maybe be replaced by nlunit passed from ccpp
   if (.not. exists) then    
      errmsg = 'lnd_iau_mod_set_control: namelist file '//trim(fn_nml)//' does not exist'
      errflg = 1
      return
   else
      land_iau_control%fn_nml = trim(fn_nml)  
      open (unit=nlunit, file=trim(fn_nml), action='READ', status='OLD', iostat=ios, iomsg=ioerrmsg)
      rewind(nlunit)
      read (nlunit, nml=land_iau_nml, err=888, END=999, iostat=ios)
      close (nlunit)
      if (ios /= 0) then      
         errmsg = 'lnd_iau_mod_set_control: error reading namelist file '//trim(fn_nml)  &
                  // 'the error message from file handler:' //trim(ioerrmsg) 
         errflg = 1
         return
       end if
   endif 
#endif
 
888 if (ios /= 0) then  
         write(iosstr, '(I0)') ios         
         errmsg = 'lnd_iau_mod_set_control: I/O error code '//trim(iosstr)//' at land_iau namelist read'  
         errflg = 1
#ifndef INTERNAL_FILE_NML
         close(nlunit)
#endif
         return
    end if 
       
999 if (ios /= 0) then  
        write(iosstr, '(I0)') ios
        if (me == mpi_root) then
          WRITE(6, * ) 'lnd_iau_mod_set_control: Warning! EoF ('//trim(iosstr)//') while reading land_iau namelist,' &
                  // ' likely because land_iau_nml was not found in input.nml. It will be set to default.' 
        endif
#ifndef INTERNAL_FILE_NML
        close(nlunit)
#endif
   endif
 
   if (me == mpi_root) then
      write(6,*) "land_iau_nml"
      write(6, land_iau_nml)
   endif
   
   land_iau_control%do_land_iau = do_land_iau
   land_iau_control%iau_delthrs = land_iau_delthrs
   land_iau_control%iau_inc_files = land_iau_inc_files
   land_iau_control%iaufhrs = land_iau_fhrs   
   land_iau_control%iau_filter_increments = land_iau_filter_increments
   land_iau_control%lsoil_incr = lsoil_incr
 
   land_iau_control%me = me
   land_iau_control%mpi_root = mpi_root
   land_iau_control%isc = isc
   land_iau_control%jsc = jsc
   land_iau_control%nx = nx
   land_iau_control%ny = ny
   land_iau_control%tile_num = tile_num
   land_iau_control%nblks = nblks
   land_iau_control%lsoil = lsoil
   land_iau_control%lsnow_lsm = lsnow_lsm
   land_iau_control%dtp = dtp
   land_iau_control%fhour = fhour
 
   land_iau_control%upd_stc = land_iau_upd_stc
   land_iau_control%upd_slc = land_iau_upd_slc
   land_iau_control%do_stcsmc_adjustment = land_iau_do_stcsmc_adjustment
   land_iau_control%min_T_increment = land_iau_min_t_increment
   land_iau_control%min_SLC_increment = land_iau_min_slc_increment
 
   allocate(land_iau_control%blksz(nblks))
   allocate(land_iau_control%blk_strt_indx(nblks))
 
   ! Land_IAU_Control%blk_strt_indx = start index of each block, for flattened (ncol=nx*ny) arrays 
   ! It's required in noahmpdriv_run to get subsection of the stc array for each proces/thread
   ix = 1
   do nb=1, nblks
      land_iau_control%blksz(nb) = blksz(nb)
      land_iau_control%blk_strt_indx(nb) = ix
      ix = ix + blksz(nb)
   enddo
 
end subroutine land_iau_mod_set_control
 
subroutine land_iau_mod_init (Land_IAU_Control, Land_IAU_Data, Land_IAU_State, errmsg, errflg)    
   type (land_iau_control_type),       intent(inout) :: land_iau_control
   type (land_iau_external_data_type), intent(inout) :: land_iau_data
   type(land_iau_state_type),          intent(inout) :: land_iau_state  
   character(len=*),                   intent(  out) :: errmsg
   integer,                            intent(  out) :: errflg
 
   ! local
   character(len=128)   :: fname
   real(kind=kind_phys) :: sx, wx, wt, normfact, dtp
   integer              :: k, nstep, kstep
   integer              :: nfilesall, ntimesall
   integer, allocatable :: idt(:) 
   integer              :: nlon, nlat
   logical              :: exists
   integer              :: ncid, dimid, varid, status, idim
   
   real(kind=kind_phys) :: dt   !, rdt
   integer              :: im, jm, km, nfiles, ntimes
 
   integer :: is,  ie,  js,  je
   integer :: npz
   integer :: i, j  
 
   !Errors messages handled through CCPP error handling variables
   errmsg = ''
   errflg = 0
 
   npz = land_iau_control%lsoil
   km = land_iau_control%lsoil
   
   is  = land_iau_control%isc
   ie  = is + land_iau_control%nx-1
   js  = land_iau_control%jsc
   je  = js + land_iau_control%ny-1
   nlon = land_iau_control%nx
   nlat = land_iau_control%ny
 
   ! allocate arrays that will hold iau state
   allocate(land_iau_data%stc_inc(nlon, nlat, km))
   allocate(land_iau_data%slc_inc(nlon, nlat, km))
 
   land_iau_data%hr1=land_iau_control%iaufhrs(1)
   land_iau_data%wt = 1.0 ! IAU increment filter weights (default 1.0)
   land_iau_data%wt_normfact = 1.0
   if (land_iau_control%iau_filter_increments) then
      ! compute increment filter weights, sum to obtain normalization factor
      dtp=land_iau_control%dtp
      nstep = 0.5*land_iau_control%iau_delthrs*3600/dtp
      ! compute normalization factor for filter weights
      normfact = 0.
      do k=1,2*nstep+1
         kstep = k-1-nstep
         sx     = acos(-1.)*kstep/nstep
         wx     = acos(-1.)*kstep/(nstep+1)
         if (kstep .ne. 0) then
            wt = sin(wx)/wx*sin(sx)/sx
         else
            wt = 1.0
         endif
         normfact = normfact + wt
         if (land_iau_control%me == land_iau_control%mpi_root) then 
            print *,'Land IAU init: IAU filter weights params k, kstep, wt ',k, kstep, wt
         endif
      enddo
      land_iau_data%wt_normfact = (2*nstep+1)/normfact
   endif
 
   ! increment files are in fv3 tiles 
   if (trim(land_iau_control%iau_inc_files(1)) .eq. '' .or. land_iau_control%iaufhrs(1) .lt. 0) then ! only 1 file expected
      errmsg = "Error! in Land IAU init: increment file name is empty or iaufhrs(1) is negative"
      errflg = 1
      return   
   endif    
   if (land_iau_control%me == land_iau_control%mpi_root) then
      print*,"Land_iau_init: Increment file name: ", trim(adjustl(land_iau_control%iau_inc_files(1)))
   endif
 
   ! determine number of valid forecast hours; read from the increment file ("Time" dim)
   if (land_iau_control%me == land_iau_control%mpi_root) then
      print *, "Land_iau_init: timesetps and forecast times (in hours) with valid increment values"
   endif
   ntimesall = size(land_iau_control%iaufhrs)
   ntimes = 0
   do k=1,ntimesall
      if (land_iau_control%iaufhrs(k) .lt. 0) exit
      if (land_iau_control%me == land_iau_control%mpi_root) then
         print *,k, " fhour ", land_iau_control%iaufhrs(k)
      endif
      ntimes = ntimes + 1
   enddo
 
   land_iau_control%ntimes = ntimes
   if (ntimes < 1) then
      errmsg = "Error! in Land IAU init: ntimes < 1 (no valid hour with increments); do_land_iau should not be .true."
      errflg = 1
      return
   endif
   if (ntimes > 1) then
      allocate(idt(ntimes-1))
      idt = land_iau_control%iaufhrs(2:ntimes)-land_iau_control%iaufhrs(1:ntimes-1)
      do k=1,ntimes-1
         if (idt(k) .ne. land_iau_control%iaufhrs(2)-land_iau_control%iaufhrs(1)) then
            errmsg = 'Fatal error in land_iau_init. forecast intervals in iaufhrs must be constant'
            errflg = 1
            return
         endif
      enddo
      deallocate(idt)
   endif
   dt = (land_iau_control%iau_delthrs*3600.)
   land_iau_data%rdt = 1.0/dt   !rdt
   
   ! Read all increment files at iau init time (at beginning of cycle) 
   ! increments are already in the fv3 grid--no need for interpolation     
   call read_iau_forcing_fv3(land_iau_control, land_iau_state%stc_inc, land_iau_state%slc_inc, errmsg, errflg)  
   if (errflg .ne. 0) return
 
   if (ntimes.EQ.1) then  ! only need to get incrments once since constant forcing over window
      call setiauforcing(land_iau_control, land_iau_data, land_iau_state)
      land_iau_data%itnext = 0
   endif
   if (ntimes.GT.1) then  !have increments at multiple forecast hours, 
      ! but only need 2 at a time and interpoalte for timesteps between them  
      ! interpolation is done in land_iau_mod_getiauforcing 
      land_iau_data%hr2=land_iau_control%iaufhrs(2)   
      land_iau_data%itnext = 2
   endif
 
end subroutine land_iau_mod_init
 
subroutine land_iau_mod_finalize(Land_IAU_Control, Land_IAU_Data, Land_IAU_state, errmsg, errflg)
 
   implicit none
 
   type(land_iau_control_type),           intent(in) :: land_iau_control
   type(land_iau_external_data_type),  intent(inout) :: land_iau_data
   type(land_iau_state_type),          intent(inout) :: land_iau_state
   character(len=*),                     intent(out) :: errmsg
   integer,                              intent(out) :: errflg
 
   ! Initialize CCPP error handling variables
   errmsg = ''
   errflg = 0
 
   if (allocated(land_iau_data%stc_inc)) deallocate (land_iau_data%stc_inc)
   if (allocated(land_iau_data%slc_inc)) deallocate (land_iau_data%slc_inc)
 
   if (allocated(land_iau_state%stc_inc)) deallocate(land_iau_state%stc_inc)
   if (allocated(land_iau_state%slc_inc)) deallocate(land_iau_state%slc_inc)
 
end subroutine land_iau_mod_finalize
 
 subroutine land_iau_mod_getiauforcing(Land_IAU_Control, Land_IAU_Data, Land_IAU_State, errmsg, errflg)
 
   implicit none
   type(land_iau_control_type),        intent(inout) :: land_iau_control
   type(land_iau_external_data_type),  intent(inout) :: land_iau_data
   type(land_iau_state_type),             intent(in) :: land_iau_state
   character(len=*),                     intent(out) :: errmsg
   integer,                              intent(out) :: errflg
   real(kind=kind_phys) t1,t2,sx,wx,wt,dtp
   integer n,i,j,k,kstep,nstep     
   integer :: ntimes
 
    ! Initialize CCPP error handling variables
   errmsg = ''
   errflg = 0
 
   ntimes = land_iau_control%ntimes
 
   land_iau_data%in_interval=.false.
   if (ntimes.LE.0) then
       errmsg = 'called land_iau_mod_getiauforcing, but ntimes <=0, probably there is no increment file. Exiting.'
       errflg = 1
       return
   endif
 
   if (ntimes .eq. 1) then 
       t1 = land_iau_control%iaufhrs(1)-0.5*land_iau_control%iau_delthrs
       t2 = land_iau_control%iaufhrs(1)+0.5*land_iau_control%iau_delthrs
   else
       t1 = land_iau_control%iaufhrs(1)
       t2 = land_iau_control%iaufhrs(ntimes)
   endif
   if (land_iau_control%iau_filter_increments) then
      ! compute increment filter weight
      ! t1 is beginning of window, t2 end of window, and Land_IAU_Control%fhour is current time
      ! in window kstep=-nstep,nstep (2*nstep+1 total) with time step of Land_IAU_Control%dtp
      dtp=land_iau_control%dtp
      nstep = 0.5*land_iau_control%iau_delthrs*3600/dtp
      ! compute normalized filter weight
      kstep = ((land_iau_control%fhour-t1) - 0.5*land_iau_control%iau_delthrs)*3600./dtp
      if (land_iau_control%fhour >= t1 .and. land_iau_control%fhour < t2) then
         sx     = acos(-1.)*kstep/nstep
         wx     = acos(-1.)*kstep/(nstep+1)
         if (kstep .ne. 0) then
            wt = (sin(wx)/wx*sin(sx)/sx)
         else
            wt = 1.
         endif
         land_iau_data%wt = land_iau_data%wt_normfact*wt
      else
         land_iau_data%wt = 0.
      endif
   endif
 
   if (ntimes.EQ.1) then
   !  check to see if we are in the IAU window, no need to update the states since they are fixed over the window
      if ( land_iau_control%fhour <= t1 .or. land_iau_control%fhour > t2 ) then
         land_iau_data%in_interval=.false.
      else
         land_iau_data%in_interval=.true.
         if (land_iau_control%me == land_iau_control%mpi_root) then 
            print *,'land_iau_mod_getiauforcing: applying forcing at t for t1,t,t2,filter wt rdt ', &
            t1,land_iau_control%fhour,t2,land_iau_data%wt/land_iau_data%wt_normfact,land_iau_data%rdt
         endif         
         if (land_iau_control%iau_filter_increments) call setiauforcing(land_iau_control, land_iau_data, land_iau_state)       
      endif
      return
   endif
 
   if (ntimes > 1) then
      if ( land_iau_control%fhour <= t1 .or. land_iau_control%fhour > t2 ) then
         land_iau_data%in_interval=.false.
      else
         land_iau_data%in_interval=.true.
         if (land_iau_control%me == land_iau_control%mpi_root) then 
            print *,'land_iau_mod_getiauforcing: applying forcing at t for t1,t,t2,filter wt rdt ', &
            t1,land_iau_control%fhour,t2,land_iau_data%wt/land_iau_data%wt_normfact,land_iau_data%rdt
         endif                 
         if (land_iau_control%fhour > land_iau_data%hr2) then ! need to read in next increment file
            land_iau_data%itnext = land_iau_data%itnext + 1
            land_iau_data%hr1=land_iau_data%hr2
            land_iau_data%hr2=land_iau_control%iaufhrs(land_iau_data%itnext)
         endif
         
         call updateiauforcing(land_iau_control, land_iau_data, land_iau_state)
      endif
   endif
 
 end subroutine land_iau_mod_getiauforcing
 
subroutine updateiauforcing(Land_IAU_Control, Land_IAU_Data, Land_IAU_State)
 
   implicit none
   
   type (land_iau_control_type),          intent(in) :: Land_IAU_Control
   type(land_iau_external_data_type),  intent(inout) :: Land_IAU_Data
   type(land_iau_state_type),             intent(in) :: Land_IAU_State
   real(kind=kind_phys) delt_t  
   integer i,j,k
   integer :: is,  ie,  js,  je, npz, t1, t2
  
   t2 = land_iau_data%itnext
   t1  = t2 - 1
   is  = 1  ! Land_IAU_Control%isc
   ie  = is + land_iau_control%nx-1
   js  = 1  ! Land_IAU_Control%jsc
   je  = js + land_iau_control%ny-1
   npz = land_iau_control%lsoil
 
   delt_t = (land_iau_data%hr2-(land_iau_control%fhour))/(land_iau_data%hr2-land_iau_data%hr1)
   
   do j = js,je
      do i = is,ie
         do k = 1,npz  ! do k = 1,n_soill    !         
            land_iau_data%stc_inc(i,j,k)  =(delt_t*land_iau_state%stc_inc(t1,i,j,k)  + (1.-delt_t)* land_iau_state%stc_inc(t2,i,j,k))*land_iau_data%rdt*land_iau_data%wt
            land_iau_data%slc_inc(i,j,k)  =(delt_t*land_iau_state%slc_inc(t1,i,j,k)  + (1.-delt_t)* land_iau_state%slc_inc(t2,i,j,k))*land_iau_data%rdt*land_iau_data%wt
         end do
       enddo
   enddo
 end subroutine updateiauforcing
 
 subroutine setiauforcing(Land_IAU_Control, Land_IAU_Data, Land_IAU_State)
 
   implicit none
   type(land_iau_control_type),       intent(in   ) :: Land_IAU_Control
   type(land_iau_external_data_type), intent(inout) :: Land_IAU_Data
   type(land_iau_state_type),         intent(in   ) :: Land_IAU_State
   real(kind=kind_phys) delt
   integer i, j, k
   integer :: is,  ie,  js,  je, npz
   
   is  = 1 
   ie  = is + land_iau_control%nx-1
   js  = 1 
   je  = js + land_iau_control%ny-1
   npz = land_iau_control%lsoil
 
   do j = js, je
      do i = is, ie
         do k = 1, npz          
            land_iau_data%stc_inc(i,j,k) = land_iau_data%wt*land_iau_state%stc_inc(1,i,j,k)*land_iau_data%rdt
            land_iau_data%slc_inc(i,j,k) = land_iau_data%wt*land_iau_state%slc_inc(1,i,j,k)*land_iau_data%rdt
         end do
      enddo
   enddo
 
 end subroutine setiauforcing
 
subroutine read_iau_forcing_fv3(Land_IAU_Control, wk3_stc, wk3_slc, errmsg, errflg)  
 
   type (land_iau_control_type),       intent(in) :: Land_IAU_Control
   real(kind=kind_phys), allocatable, intent(out) :: wk3_stc(:, :, :, :), wk3_slc(:, :, :, :)
   character(len=*),          intent(inout) :: errmsg
   integer,                   intent(inout) :: errflg
 
   integer  :: i, it, km 
   logical  :: exists
   integer  :: ncid, status, varid
   integer  :: ierr
   character(len=500)  :: fname
   character(len=2)    :: tile_str
   integer             :: n_t, n_y, n_x
 
   character(len=32), dimension(4) :: stc_vars = [character(len=32) :: 'soilt1_inc', 'soilt2_inc', 'soilt3_inc', 'soilt4_inc']
   character(len=32), dimension(4) :: slc_vars = [character(len=32) :: 'slc1_inc', 'slc2_inc', 'slc3_inc', 'slc4_inc']
   character(len=32) :: slsn_mask = "soilsnow_mask"
 
   !Errors messages handled through CCPP error handling variables
   errmsg = ''
   errflg = 0
 
   km = land_iau_control%lsoil
 
   write(tile_str, '(I0)') land_iau_control%tile_num
 
   fname = 'INPUT/'//trim(land_iau_control%iau_inc_files(1))//".tile"//trim(tile_str)//".nc"
   
   inquire (file=trim(fname), exist=exists)    
   if (exists) then
      status = nf90_open(trim(fname), nf90_nowrite, ncid)  ! open the file
      call netcdf_err(status, ' opening file '//trim(fname), errflg, errmsg) 
      if (errflg .ne. 0) return
   else
      errmsg = 'FATAL Error in land iau read_iau_forcing_fv3: Expected file '//trim(fname)//' for DA increment does not exist'
      errflg = 1
      return
   endif
   ! var stored as soilt1_inc(Time, yaxis_1, xaxis_1)
   call get_nc_dimlen(ncid, "Time", n_t, errflg, errmsg) 
   if (errflg .ne. 0) return
   call get_nc_dimlen(ncid, "yaxis_1", n_y, errflg, errmsg) 
   if (errflg .ne. 0) return
   call get_nc_dimlen(ncid, "xaxis_1", n_x, errflg, errmsg) 
   if (errflg .ne. 0) return
 
   if (n_x .lt. land_iau_control%nx) then 
      errmsg = 'Error in land iau read_iau_forcing_fv3: Land_IAU_Control%nx bigger than dim xaxis_1 in file '//trim(fname)
      errflg = 1
      return
   endif
   if (n_y .lt. land_iau_control%ny) then 
      errmsg = 'Error in land iau read_iau_forcing_fv3: Land_IAU_Control%ny bigger than dim yaxis_1 in file '//trim(fname)
      errflg = 1
      return
   endif
 
   allocate(wk3_stc(n_t, land_iau_control%nx, land_iau_control%ny, km))
   allocate(wk3_slc(n_t, land_iau_control%nx, land_iau_control%ny, km))
   
   do i = 1, size(stc_vars)
      if (land_iau_control%me == land_iau_control%mpi_root) print *, trim(stc_vars(i))
      status = nf90_inq_varid(ncid, trim(stc_vars(i)), varid)
      if (status == nf90_noerr) then   
         do it = 1, n_t
            ! var stored as soilt1_inc(Time, yaxis_1, xaxis_1)
            call get_var3d_values(ncid, varid, trim(stc_vars(i)), land_iau_control%isc, land_iau_control%nx, land_iau_control%jsc, land_iau_control%ny, &
                 it, 1, wk3_stc(it,:, :, i), status, errflg, errmsg)
            if (errflg .ne. 0) return 
         enddo
      else
         if (land_iau_control%me == land_iau_control%mpi_root) then 
            print *, 'warning! No increment for ',trim(stc_vars(i)),' found, assuming zero'
         endif
         wk3_stc(:, :, :, i) = 0.
      endif
   enddo
   do i = 1, size(slc_vars)
      if (land_iau_control%me == land_iau_control%mpi_root) print *, trim(slc_vars(i))
      status = nf90_inq_varid(ncid, trim(slc_vars(i)), varid)
      if (status == nf90_noerr) then   
         do it = 1, n_t
            call get_var3d_values(ncid, varid, trim(slc_vars(i)), land_iau_control%isc, land_iau_control%nx, land_iau_control%jsc, land_iau_control%ny, &
                 it, 1, wk3_slc(it, :, :, i), status, errflg, errmsg)
            if (errflg .ne. 0) return   
         end do       
      else
         if (land_iau_control%me == land_iau_control%mpi_root) then 
            print *, 'warning! No increment for ',trim(slc_vars(i)),' found, assuming zero'
         endif
         wk3_slc(:, :, :, i) = 0.
      endif
   enddo
   
   !set too small increments to zero
   where(abs(wk3_stc) < land_iau_control%min_T_increment) wk3_stc = 0.0
   where(abs(wk3_slc) < land_iau_control%min_SLC_increment) wk3_slc = 0.0
 
   status =nf90_close(ncid) 
   call netcdf_err(status, 'closing file '//trim(fname), errflg, errmsg) 
 
end subroutine read_iau_forcing_fv3
 
   subroutine calculate_landinc_mask(swe,vtype,stype,lensfc,veg_type_landice, mask)
   
      implicit none
 
      integer, intent(in)           :: lensfc, veg_type_landice
      real(kind=kind_phys), intent(in)              :: swe(lensfc)
      integer, intent(in)           :: vtype(lensfc),stype(lensfc)
      integer, intent(out)          :: mask(lensfc)
 
      integer :: i
 
      mask = -1 ! not land
 
      ! land (but not land-ice)
      do i=1,lensfc
         if (stype(i) .GT. 0) then
            if (swe(i) .GT. 0.001) then ! snow covered land
                  mask(i) = 2
            else                        ! non-snow covered land
                  mask(i) = 1
            endif
         end if ! else should work here too
         if ( vtype(i) ==  veg_type_landice  ) then ! land-ice
                  mask(i) = 0
         endif
      end do
 
   end subroutine calculate_landinc_mask
 
   subroutine netcdf_err(ERR, STRING, errflg, errmsg_out)
 
   !--------------------------------------------------------------
   ! Process the error flag from a NETCDF call and return it as (human readable) MESSAGE
   !--------------------------------------------------------------
      IMPLICIT NONE
 
      include 'mpif.h'
 
      INTEGER, INTENT(IN) :: ERR
      CHARACTER(LEN=*), INTENT(IN) :: STRING
      CHARACTER(LEN=80) :: ERRMSG
      integer :: errflg
      character(len=*) :: errmsg_out
 
      !Errors messages handled through CCPP error handling variables
      errmsg_out = ''
      errflg = 0
 
      IF (err == nf90_noerr) RETURN
      errmsg = nf90_strerror(err)
      errmsg_out = 'FATAL ERROR in Land IAU '//trim(string)//': '//trim(errmsg)
      errflg = 1
      return
 
   end subroutine netcdf_err
 
   subroutine get_nc_dimlen(ncid, dim_name, dim_len, errflg, errmsg_out )
      integer, intent(in):: ncid
      character(len=*), intent(in)::  dim_name
      integer, intent(out):: dim_len
      integer :: dimid
      integer :: errflg
      character(len=*) :: errmsg_out
      integer :: status
 
      !Errors messages handled through CCPP error handling variables
      errmsg_out = ''
      errflg = 0
 
      status = nf90_inq_dimid(ncid, dim_name, dimid)
      CALL netcdf_err(status, 'reading dim id '//trim(dim_name), errflg, errmsg_out)
      if (errflg .ne. 0) return
      status = nf90_inquire_dimension(ncid, dimid, len = dim_len)
      CALL netcdf_err(status, 'reading dim length '//trim(dim_name), errflg, errmsg_out)
 
   end subroutine get_nc_dimlen
 
   subroutine get_var1d(ncid, dim_len, var_name, var_arr, errflg, errmsg_out)
      integer, intent(in):: ncid, dim_len
      character(len=*), intent(in)::  var_name
      real(kind=kind_phys), intent(out):: var_arr(dim_len)
      integer :: errflg
      character(len=*) :: errmsg_out
      integer :: varid, status
 
      !Errors messages handled through CCPP error handling variables
      errmsg_out = ''
      errflg = 0
 
      status = nf90_inq_varid(ncid, trim(var_name), varid)
      call netcdf_err(status, 'getting varid: '//trim(var_name), errflg, errmsg_out)
      if (errflg .ne. 0) return
 
      status = nf90_get_var(ncid, varid, var_arr)
      call netcdf_err(status, 'reading var: '//trim(var_name), errflg, errmsg_out)
 
   end subroutine get_var1d
 
   subroutine get_var3d_values(ncid, varid, var_name, is,ix, js,jy, ks,kz, var3d, status, errflg, errmsg_out)
      integer, intent(in):: ncid, varid
      integer, intent(in):: is, ix, js, jy, ks,kz
      character(len=*), intent(in)::  var_name
      real(kind=kind_phys), intent(out):: var3d(ix, jy, kz)   
      integer, intent(out):: status 
      integer :: errflg
      character(len=*) :: errmsg_out
 
      !Errors messages handled through CCPP error handling variables
      errmsg_out = ''
      errflg = 0
 
      status = nf90_get_var(ncid, varid, var3d, &  
               start = (/is, js, ks/), count = (/ix, jy, kz/))
 
      call netcdf_err(status, 'get_var3d_values '//trim(var_name), errflg, errmsg_out)
      
 
   end subroutine get_var3d_values
 
   subroutine get_var3d_values_int(ncid, varid, var_name, is,ix, js,jy, ks,kz, var3d, status, errflg, errmsg_out)
      integer, intent(in):: ncid, varid
      integer, intent(in):: is, ix, js, jy, ks,kz
      character(len=*), intent(in)::  var_name
      integer, intent(out):: var3d(ix, jy, kz)   
      integer, intent(out):: status 
      integer :: errflg
      character(len=*) :: errmsg_out
 
      !Errors messages handled through CCPP error handling variables
      errmsg_out = ''
      errflg = 0
 
      status = nf90_get_var(ncid, varid, var3d, &  !start = start, count = nreco)
               start = (/is, js, ks/), count = (/ix, jy, kz/))
      
      call netcdf_err(status, 'get_var3d_values_int '//trim(var_name), errflg, errmsg_out)
 
   end subroutine get_var3d_values_int
  
end module land_iau_mod