diff --git a/physics/SFC_Models/Land/Noahmp/lnd_iau_mod.F90 b/physics/SFC_Models/Land/Noahmp/lnd_iau_mod.F90
new file mode 100644
index 000000000..40f3eb8f7
--- /dev/null
+++ b/physics/SFC_Models/Land/Noahmp/lnd_iau_mod.F90
@@ -0,0 +1,788 @@
+!***********************************************************************
+!> TODO: replace with appropriate licence for CCPP
+!* GNU Lesser General Public License
+!* .
+!***********************************************************************
+
+!> @brief Land IAU (Incremental Analysis Update) module,
+!> for the NoahMP soil/snow temperature (can be extended to include soil moisture)
+
+!! \section land_iau_mod
+!> - reads settings from namelist file (which indicates if IAU increments are available or not)
+!> - reads in DA increments from GSI/JEDI DA at the start of (the DA) cycle
+!> - maps increments to FV3 grid points belonging to mpi process
+!> - interpolates temporally (with filter-weights if required by configuration)
+!> - updates states with the interpolated increments
+
+!> March, 2024: Tseganeh Z. Gichamo, (EMC) based on the FV3 IAU mod
+!> by Xi.Chen and Philip Pegion, PSL
+!-------------------------------------------------------------------------------
+
+!> \section arg_table_land_iau_mod Argument table
+!! \htmlinclude land_iau_mod.html
+!!
+module land_iau_mod
+
+ use machine, only: kind_phys, kind_dyn
+ use netcdf
+
+ implicit none
+
+ private
+
+!> \section arg_table_land_iau_external_data_type Argument Table
+!! \htmlinclude land_iau_external_data_type.html
+!!
+ 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 !< number of soil layers
+ integer :: lsnow_lsm !< maximum number of snow layers internal to land surface model
+ 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
+
+ integer :: me !< MPI rank designator
+ integer :: mpi_root !< MPI rank of master atmosphere processor
+ character(len=64) :: fn_nml !< namelist filename for surface data cycling
+ real(kind=kind_phys) :: dtp !< physics timestep in seconds
+ real(kind=kind_phys) :: fhour !< current forecast hour
+
+ 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_i, 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 !< namelist filename for surface data cycling
+ character(len=:), intent(in), dimension(:), pointer :: input_nml_file_i
+ integer, intent(in) :: me, mpi_root !< MPI rank of master atmosphere processor
+ 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 !< physics timestep in seconds
+ real(kind=kind_phys), intent(in) :: fhour !< current forecast hour
+ 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=:), pointer, dimension(:) :: input_nml_file => null()
+ character(len=4) :: iosstr
+
+ !> land iau setting read from namelist
+ logical :: do_land_iau = .false.
+ real(kind=kind_phys) :: land_iau_delthrs = 0 !< iau time interval (to scale increments)
+ character(len=240) :: land_iau_inc_files(7) = '' !< list of increment files
+ real(kind=kind_phys) :: land_iau_fhrs(7) = -1 !< forecast hours associated with increment files
+ logical :: land_iau_filter_increments = .false. !< filter IAU increments
+
+ 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
+
+ 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
+
+ !Errors messages handled through CCPP error handling variables
+ errmsg = ''
+ errflg = 0
+
+!3.11.24: copied from GFS_typedefs.F90
+#ifdef INTERNAL_FILE_NML
+ ! allocate required to work around GNU compiler bug 100886
+ ! https://gcc.gnu.org/bugzilla/show_bug.cgi?id=100886
+ allocate(input_nml_file, mold=input_nml_file_i)
+ input_nml_file => input_nml_file_i
+ 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
+ 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
+ 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
+
+ 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
+
+ status =nf90_close(ncid)
+ call netcdf_err(status, 'closing file '//trim(fname), errflg, errmsg)
+
+end subroutine read_iau_forcing_fv3
+
+ !> Calculate soil mask for land on model grid.
+ !! Output is 1 - soil, 2 - snow-covered, 0 - land ice, -1 not land.
+ !!
+ !! @param[in] lensfc Number of land points for this tile
+ !! @param[in] veg_type_landice Value of vegetion class that indicates land-ice
+ !! @param[in] stype Soil type
+ !! @param[in] swe Model snow water equivalent
+ !! @param[in] vtype Model vegetation type
+ !! @param[out] mask Land mask for increments
+ !! @author Clara Draper @date March 2021
+ !! @author Yuan Xue: introduce stype to make the mask calculation more generic
+ 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
+
+
diff --git a/physics/SFC_Models/Land/Noahmp/lnd_iau_mod.meta b/physics/SFC_Models/Land/Noahmp/lnd_iau_mod.meta
new file mode 100644
index 000000000..8541af659
--- /dev/null
+++ b/physics/SFC_Models/Land/Noahmp/lnd_iau_mod.meta
@@ -0,0 +1,58 @@
+[ccpp-table-properties]
+ name = land_iau_external_data_type
+ type = ddt
+ dependencies =
+
+[ccpp-arg-table]
+ name = land_iau_external_data_type
+ type = ddt
+
+########################################################################
+
+[ccpp-table-properties]
+ name = land_iau_state_type
+ type = ddt
+ dependencies =
+
+[ccpp-arg-table]
+ name = land_iau_state_type
+ type = ddt
+
+########################################################################
+
+[ccpp-table-properties]
+ name = land_iau_control_type
+ type = ddt
+ dependencies =
+
+[ccpp-arg-table]
+ name = land_iau_control_type
+ type = ddt
+
+########################################################################
+[ccpp-table-properties]
+ name = land_iau_mod
+ type = module
+ dependencies = machine.F
+
+[ccpp-arg-table]
+ name = land_iau_mod
+ type = module
+[land_iau_external_data_type]
+ standard_name = land_iau_external_data_type
+ long_name = definition of type land_iau_external_data_type
+ units = DDT
+ dimensions = ()
+ type = land_iau_external_data_type
+[land_iau_state_type]
+ standard_name = land_iau_state_type
+ long_name = definition of type land_iau_state_type
+ units = DDT
+ dimensions = ()
+ type = land_iau_state_type
+[land_iau_control_type]
+ standard_name = land_iau_control_type
+ long_name = definition of type land_iau_control_type
+ units = DDT
+ dimensions = ()
+ type = land_iau_control_type
diff --git a/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 b/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90
index 1313e9ff3..d4971efd9 100644
--- a/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90
+++ b/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90
@@ -13,13 +13,19 @@ module noahmpdrv
use module_sf_noahmplsm
+! These hold and apply Land IAU increments for soil temperature
+! (possibly will extend to soil moisture increments)
+ use land_iau_mod, only: land_iau_control_type, land_iau_external_data_type, land_iau_state_type, &
+ land_iau_mod_init, land_iau_mod_getiauforcing, land_iau_mod_finalize, calculate_landinc_mask
+
implicit none
integer, parameter :: psi_opt = 0 ! 0: MYNN or 1:GFS
private
- public :: noahmpdrv_init, noahmpdrv_run
+ public :: noahmpdrv_init, noahmpdrv_run, &
+ noahmpdrv_timestep_init, noahmpdrv_finalize
contains
@@ -32,7 +38,8 @@ module noahmpdrv
subroutine noahmpdrv_init(lsm, lsm_noahmp, me, isot, ivegsrc, &
nlunit, pores, resid, &
do_mynnsfclay,do_mynnedmf, &
- errmsg, errflg)
+ errmsg, errflg, &
+ Land_IAU_Control, Land_IAU_Data, Land_IAU_state)
use machine, only: kind_phys
use set_soilveg_mod, only: set_soilveg
@@ -53,6 +60,19 @@ subroutine noahmpdrv_init(lsm, lsm_noahmp, me, isot, ivegsrc, &
character(len=*), intent(out) :: errmsg
integer, intent(out) :: errflg
+ ! Land iau mod DDTs ! made optional to allow NoahMP Component model call this function without having to deal with IAU
+
+ ! Land IAU Control holds settings' information, maily read from namelist
+ ! (e.g., block of global domain that belongs to current process,
+ ! whether to do IAU increment at this time step, time step informatoin, etc)
+ type(land_iau_control_type), intent(inout), optional :: Land_IAU_Control
+
+ ! land iau state holds increment data read from file (before interpolation)
+ type(land_iau_state_type), intent(inout), optional :: Land_IAU_state
+
+ ! Land IAU Data holds spatially and temporally interpolated increments per time step
+ type(land_iau_external_data_type), intent(inout), optional :: Land_IAU_Data ! arry of (number of blocks):each proc holds nblks
+
! Initialize CCPP error handling variables
errmsg = ''
errflg = 0
@@ -100,9 +120,282 @@ subroutine noahmpdrv_init(lsm, lsm_noahmp, me, isot, ivegsrc, &
pores (:) = maxsmc (:)
resid (:) = drysmc (:)
+
+ if (present(Land_IAU_Control) .and. present(Land_IAU_Data) .and. present(Land_IAU_State)) then
+
+ ! Initialize IAU for land--land_iau_control was set by host model
+ if (.not. Land_IAU_Control%do_land_iau) return
+ call land_iau_mod_init (Land_IAU_Control, Land_IAU_Data, Land_IAU_State, errmsg, errflg)
+
+ endif
end subroutine noahmpdrv_init
+!> \ingroup NoahMP_LSM
+!! \brief This subroutine is called before noahmpdrv_run
+!! to update states with iau increments, if available
+!! \section arg_table_noahmpdrv_timestep_init Argument Table
+!! \htmlinclude noahmpdrv_timestep_init.html
+!!
+subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, &
+ isot, ivegsrc, soiltyp, vegtype, weasd, &
+ land_iau_control, land_iau_data, land_iau_state, &
+ stc, slc, smc, errmsg, errflg, &
+ con_g, con_t0c, con_hfus)
+
+ use machine, only: kind_phys
+ use namelist_soilveg
+ ! use set_soilveg_snippet_mod, only: set_soilveg_noahmp
+ use noahmp_tables
+
+ implicit none
+
+ integer , intent(in) :: itime !current forecast iteration
+ real(kind=kind_phys) , intent(in) :: fhour !current forecast time (hr)
+ real(kind=kind_phys) , intent(in) :: delt ! time interval [s]
+ integer , intent(in) :: km !vertical soil layer dimension
+ integer, intent(in) :: ncols
+ integer, intent(in) :: isot
+ integer, intent(in) :: ivegsrc
+
+ integer , dimension(:) , intent(in) :: soiltyp ! soil type (integer index)
+ integer , dimension(:) , intent(in) :: vegtype ! vegetation type (integer index)
+ real(kind=kind_phys), dimension(:) , intent(inout) :: weasd ! water equivalent accumulated snow depth [mm]
+
+ 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
+ real(kind=kind_phys), dimension(:,:) , intent(inout) :: stc ! soiltemp [K]
+ real(kind=kind_phys), dimension(:,:) , intent(inout) :: slc !liquid soil moisture [m3/m3]'
+ real(kind=kind_phys), dimension(:,:) , intent(inout) :: smc !
+ character(len=*), intent(out) :: errmsg
+ integer, intent(out) :: errflg
+ real(kind=kind_phys), intent(in) :: con_g ! grav
+ real(kind=kind_phys), intent(in) :: con_t0c ! tfreez
+ real(kind=kind_phys), intent(in) :: con_hfus ! hfus
+
+ ! IAU update
+ real(kind=kind_phys),allocatable, dimension(:,:) :: stc_inc_flat, slc_inc_flat
+ real(kind=kind_phys), dimension(km) :: dz ! layer thickness
+
+!TODO: This is hard-coded in noahmpdrv
+ real(kind=kind_phys) :: zsoil(4) = (/ -0.1, -0.4, -1.0, -2.0 /) !zsoil(km)
+
+ integer :: lsoil_incr
+ integer, allocatable :: mask_tile(:)
+ integer,allocatable :: stc_updated(:), slc_updated(:)
+ logical :: soil_freeze, soil_ice
+ integer :: soiltype, n_stc, n_slc
+ real(kind=kind_phys) :: slc_new
+
+ integer :: i, j, ij, l, k, ib
+ integer :: lensfc
+
+ real(kind=kind_phys) :: smp !< for computing supercooled water
+ real(kind=kind_phys) :: hc_incr
+
+ integer :: nother, nsnowupd
+ integer :: nstcupd, nslcupd, nfrozen, nfrozen_upd
+ logical :: print_update_stats = .False.
+
+ ! --- Initialize CCPP error handling variables
+ errmsg = ''
+ errflg = 0
+
+ if (.not. Land_IAU_Control%do_land_iau) return
+
+ !> update current forecast hour
+ Land_IAU_Control%fhour=fhour
+
+ !> read iau increments
+ call land_iau_mod_getiauforcing(Land_IAU_Control, Land_IAU_Data, Land_IAU_state, errmsg, errflg)
+ if (errflg .ne. 0) then
+ return
+ endif
+
+ !> If no increment at the current timestep simply proceed forward
+ if (.not. Land_IAU_Data%in_interval) then
+ return
+ endif
+
+ if(Land_IAU_Control%me == Land_IAU_Control%mpi_root) then
+ print*, "adding land iau increments "
+ endif
+
+ if (Land_IAU_Control%lsoil .ne. km) then
+ write(errmsg,*) 'noahmpdrv_timestep_init: Land_IAU_Data%lsoil ',Land_IAU_Control%lsoil,' not equal to km ',km
+ errflg = 1
+ return
+ endif
+
+ ! local variable to copy blocked data Land_IAU_Data%stc_inc
+ allocate(stc_inc_flat(Land_IAU_Control%nx * Land_IAU_Control%ny, km)) !GFS_Control%ncols
+ allocate(slc_inc_flat(Land_IAU_Control%nx * Land_IAU_Control%ny, km)) !GFS_Control%ncols
+ allocate(stc_updated(Land_IAU_Control%nx * Land_IAU_Control%ny))
+ allocate(slc_updated(Land_IAU_Control%nx * Land_IAU_Control%ny))
+
+ !copy background stc
+ stc_updated = 0
+ slc_updated = 0
+ ib = 1
+ do j = 1, Land_IAU_Control%ny
+ do k = 1, km
+ stc_inc_flat(ib:ib+Land_IAU_Control%nx-1, k) = Land_IAU_Data%stc_inc(:,j, k)
+ slc_inc_flat(ib:ib+Land_IAU_Control%nx-1, k) = Land_IAU_Data%slc_inc(:,j, k)
+ enddo
+ ib = ib + Land_IAU_Control%nx
+ enddo
+
+ if ((Land_IAU_Control%dtp - delt) > 0.0001) then
+ if(Land_IAU_Control%me == Land_IAU_Control%mpi_root) then
+ print*, "Warning! noahmpdrv_timestep_init delt ",delt," different from Land_IAU_Control%dtp ",Land_IAU_Control%dtp
+ endif
+ endif
+
+ lsoil_incr = Land_IAU_Control%lsoil_incr
+ lensfc = Land_IAU_Control%nx * Land_IAU_Control%ny
+
+ if(Land_IAU_Control%me == Land_IAU_Control%mpi_root) print*,' adjusting first ', lsoil_incr, ' surface layers only, delt ', delt
+ ! initialize variables for counts statitics to be zeros
+ nother = 0 ! grid cells not land
+ nsnowupd = 0 ! grid cells with snow (temperature not yet updated)
+ nstcupd = 0 ! grid cells that are updated stc
+ nslcupd = 0 ! grid cells that are updated slc
+ nfrozen = 0 ! not update as frozen soil
+ nfrozen_upd = 0 ! not update as frozen soil
+
+!TODO---if only fv3 increment files are used, this can be read from file
+ allocate(mask_tile(lensfc))
+ call calculate_landinc_mask(weasd, vegtype, soiltyp, lensfc, isice_table, mask_tile)
+
+ !IAU increments are in units of 1/sec !Land_IAU_Control%dtp
+ !* only updating soil temp for now
+ ij_loop : do ij = 1, lensfc
+ ! mask: 1 - soil, 2 - snow, 0 - land-ice, -1 - not land
+ if (mask_tile(ij) == 1) then
+
+ soil_freeze=.false.
+ soil_ice=.false.
+ do k = 1, lsoil_incr ! k = 1, km
+ if ( stc(ij,k) < con_t0c) soil_freeze=.true.
+ if ( smc(ij,k) - slc(ij,k) > 0.001 ) soil_ice=.true.
+
+ if (Land_IAU_Control%upd_stc) then
+ stc(ij,k) = stc(ij,k) + stc_inc_flat(ij,k)*delt !Land_IAU_Control%dtp
+ if (k==1) then
+ stc_updated(ij) = 1
+ nstcupd = nstcupd + 1
+ endif
+ endif
+
+ if ( (stc(ij,k) < con_t0c) .and. (.not. soil_freeze) .and. (k==1) ) nfrozen_upd = nfrozen_upd + 1
+
+ ! do not do updates if this layer or any above is frozen
+ if ( (.not. soil_freeze ) .and. (.not. soil_ice ) ) then
+ if (Land_IAU_Control%upd_slc) then
+ if (k==1) then
+ nslcupd = nslcupd + 1
+ slc_updated(ij) = 1
+ endif
+ ! apply zero limit here (higher, model-specific limits are later)
+ slc(ij,k) = max(slc(ij,k) + slc_inc_flat(ij,k)*delt, 0.0)
+ smc(ij,k) = max(smc(ij,k) + slc_inc_flat(ij,k)*delt, 0.0)
+ endif
+ else
+ if (k==1) nfrozen = nfrozen+1
+ endif
+ enddo
+ endif ! if soil/snow point
+ enddo ij_loop
+
+ deallocate(stc_inc_flat, slc_inc_flat)
+
+ !!do moisture/temperature adjustment for consistency after increment add
+ call read_mp_table_parameters(errmsg, errflg)
+ if (errflg .ne. 0) then
+ errmsg = 'FATAL ERROR in noahmpdrv_timestep_init: problem in set_soilveg_noahmp'
+ return
+ endif
+ n_stc = 0
+ n_slc = 0
+ if (Land_IAU_Control%do_stcsmc_adjustment) then
+ if (Land_IAU_Control%upd_stc) then
+ do i=1,lensfc
+ if (stc_updated(i) == 1 ) then ! soil-only location
+ n_stc = n_stc+1
+ soiltype = soiltyp(i)
+ do l = 1, lsoil_incr
+ !case 1: frz ==> frz, recalculate slc, smc remains
+ !case 2: unfrz ==> frz, recalculate slc, smc remains
+ !both cases are considered in the following if case
+ if (stc(i,l) .LT. con_t0c )then
+ !recompute supercool liquid water,smc_anl remain unchanged
+ smp = con_hfus*(con_t0c-stc(i,l))/(con_g*stc(i,l)) !(m)
+ slc_new=maxsmc(soiltype)*(smp/satpsi(soiltype))**(-1./bb(soiltype))
+ slc(i,l) = max( min( slc_new, smc(i,l)), 0.0 )
+ endif
+ !case 3: frz ==> unfrz, melt all soil ice (if any)
+ if (stc(i,l) .GT. con_t0c )then !do not rely on stc_bck
+ slc(i,l)=smc(i,l)
+ endif
+ enddo
+ endif
+ enddo
+ endif
+
+ if (Land_IAU_Control%upd_slc) then
+ dz(1) = -zsoil(1)
+ do l = 2, km
+ dz(l) = -zsoil(l) + zsoil(l-1)
+ enddo
+ do i=1,lensfc
+ if (slc_updated(i) == 1 ) then
+ n_slc = n_slc+1
+ ! apply SM bounds (later: add upper SMC limit)
+ do l = 1, lsoil_incr
+ ! noah-mp minimum is 1 mm per layer (in SMC)
+ ! no need to maintain frozen amount, would be v. small.
+ slc(i,l) = max( 0.001/dz(l), slc(i,l) )
+ smc(i,l) = max( 0.001/dz(l), smc(i,l) )
+ enddo
+ endif
+ enddo
+ endif
+ endif
+
+ deallocate(stc_updated, slc_updated)
+ deallocate(mask_tile)
+
+ write(*,'(a,i4,a,i8)') 'noahmpdrv_timestep_init rank ', Land_IAU_Control%me, ' # of cells with stc update ', nstcupd
+
+
+end subroutine noahmpdrv_timestep_init
+
+ !> \ingroup NoahMP_LSM
+!! \brief This subroutine mirrors noahmpdrv_init
+!! it calls land_iau_finalize which frees up allocated memory by IAU_init (in noahmdrv_init)
+!! \section arg_table_noahmpdrv_finalize Argument Table
+!! \htmlinclude noahmpdrv_finalize.html
+!!
+ subroutine noahmpdrv_finalize (Land_IAU_Control, Land_IAU_Data, Land_IAU_State, errmsg, errflg)
+
+ use machine, only: kind_phys
+ 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
+ integer :: j, k, ib
+ ! --- Initialize CCPP error handling variables
+ errmsg = ''
+ errflg = 0
+
+ if (.not. Land_IAU_Control%do_land_iau) return
+ call land_iau_mod_finalize(Land_IAU_Control, Land_IAU_Data, Land_IAU_State, errmsg, errflg)
+
+ end subroutine noahmpdrv_finalize
+
!> \ingroup NoahMP_LSM
!! \brief This subroutine is the main CCPP entry point for the NoahMP LSM.
!! \section arg_table_noahmpdrv_run Argument Table
diff --git a/physics/SFC_Models/Land/Noahmp/noahmpdrv.meta b/physics/SFC_Models/Land/Noahmp/noahmpdrv.meta
index 753550016..7d1150c80 100644
--- a/physics/SFC_Models/Land/Noahmp/noahmpdrv.meta
+++ b/physics/SFC_Models/Land/Noahmp/noahmpdrv.meta
@@ -96,6 +96,233 @@
dimensions = ()
type = integer
intent = out
+[land_iau_control]
+ standard_name = land_data_assimilation_control
+ long_name = land data assimilation control
+ units = mixed
+ dimensions = ()
+ type = land_iau_control_type
+ intent = inout
+ optional = True
+[land_iau_data]
+ standard_name = land_data_assimilation_data
+ long_name = land data assimilation data
+ units = mixed
+ dimensions = ()
+ type = land_iau_external_data_type
+ intent = inout
+ optional = True
+[land_iau_state]
+ standard_name = land_data_assimilation_interpolated_data
+ long_name = land data assimilation space- and time-interpolated
+ units = mixed
+ dimensions = ()
+ type = land_iau_state_type
+ intent = inout
+ optional = True
+
+########################################################################
+[ccpp-arg-table]
+ name = noahmpdrv_timestep_init
+ type = scheme
+[itime]
+ standard_name = index_of_timestep
+ long_name = current forecast iteration
+ units = index
+ dimensions = ()
+ type = integer
+ intent = in
+[fhour]
+ standard_name = forecast_time
+ long_name = current forecast time
+ units = h
+ dimensions = ()
+ type = real
+ kind = kind_phys
+ intent = in
+[delt]
+ standard_name = timestep_for_dynamics
+ long_name = dynamics timestep
+ units = s
+ dimensions = ()
+ type = real
+ kind = kind_phys
+ intent = in
+[km]
+ standard_name = vertical_dimension_of_soil
+ long_name = vertical dimension of soil layers
+ units = count
+ dimensions = ()
+ type = integer
+ intent = in
+[ncols]
+ standard_name = horizontal_dimension
+ long_name = horizontal dimension
+ units = count
+ dimensions = ()
+ type = integer
+ intent = in
+[isot]
+ standard_name = control_for_soil_type_dataset
+ long_name = soil type dataset choice
+ units = index
+ dimensions = ()
+ type = integer
+ intent = in
+[ivegsrc]
+ standard_name = control_for_vegetation_dataset
+ long_name = land use dataset choice
+ units = index
+ dimensions = ()
+ type = integer
+ intent = in
+[soiltyp]
+ standard_name = soil_type_classification
+ long_name = soil type at each grid cell
+ units = index
+ dimensions = (horizontal_dimension)
+ type = integer
+ intent= in
+[vegtype]
+ standard_name = vegetation_type_classification
+ long_name = vegetation type at each grid cell
+ units = index
+ dimensions = (horizontal_dimension)
+ type = integer
+ intent= in
+[weasd]
+ standard_name = water_equivalent_accumulated_snow_depth_over_land
+ long_name = water equivalent of accumulated snow depth over land
+ units = mm
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = inout
+[land_iau_control]
+ standard_name = land_data_assimilation_control
+ long_name = land data assimilation control
+ units = mixed
+ dimensions = ()
+ type = land_iau_control_type
+ intent = inout
+[land_iau_data]
+ standard_name = land_data_assimilation_data
+ long_name = land data assimilation data
+ units = mixed
+ dimensions = ()
+ type = land_iau_external_data_type
+ intent = inout
+[land_iau_state]
+ standard_name = land_data_assimilation_interpolated_data
+ long_name = land data assimilation space- and time-interpolated
+ units = mixed
+ dimensions = ()
+ type = land_iau_state_type
+ intent = inout
+[stc]
+ standard_name = soil_temperature
+ long_name = soil temperature
+ units = K
+ dimensions = (horizontal_dimension,vertical_dimension_of_soil)
+ type = real
+ kind = kind_phys
+ intent = inout
+[slc]
+ standard_name = volume_fraction_of_unfrozen_water_in_soil
+ long_name = liquid soil moisture
+ units = frac
+ dimensions = (horizontal_dimension,vertical_dimension_of_soil)
+ type = real
+ kind = kind_phys
+ intent = inout
+[smc]
+ standard_name = volume_fraction_of_condensed_water_in_soil
+ long_name = total soil moisture
+ units = frac
+ dimensions = (horizontal_dimension,vertical_dimension_of_soil)
+ type = real
+ kind = kind_phys
+ intent = inout
+[errmsg]
+ standard_name = ccpp_error_message
+ long_name = error message for error handling in CCPP
+ units = none
+ dimensions = ()
+ type = character
+ kind = len=*
+ intent = out
+[errflg]
+ standard_name = ccpp_error_code
+ long_name = error code for error handling in CCPP
+ units = 1
+ dimensions = ()
+ type = integer
+ intent = out
+[con_g]
+ standard_name = gravitational_acceleration
+ long_name = gravitational acceleration
+ units = m s-2
+ dimensions = ()
+ type = real
+ kind = kind_phys
+ intent = in
+[con_t0c]
+ standard_name = temperature_at_zero_celsius
+ long_name = temperature at 0 degree Celsius
+ units = K
+ dimensions = ()
+ type = real
+ kind = kind_phys
+ intent = in
+[con_hfus]
+ standard_name = latent_heat_of_fusion_of_water_at_0C
+ long_name = latent heat of fusion
+ units = J kg-1
+ dimensions = ()
+ type = real
+ kind = kind_phys
+ intent = in
+
+#######################################################################
+[ccpp-arg-table]
+ name = noahmpdrv_finalize
+ type = scheme
+[land_iau_control]
+ standard_name = land_data_assimilation_control
+ long_name = land data assimilation control
+ units = mixed
+ dimensions = ()
+ type = land_iau_control_type
+ intent = in
+[land_iau_data]
+ standard_name = land_data_assimilation_data
+ long_name = land data assimilation data
+ units = mixed
+ dimensions = ()
+ type = land_iau_external_data_type
+ intent = inout
+[land_iau_state]
+ standard_name = land_data_assimilation_interpolated_data
+ long_name = land data assimilation space- and time-interpolated
+ units = mixed
+ dimensions = ()
+ type = land_iau_state_type
+ intent = inout
+[errmsg]
+ standard_name = ccpp_error_message
+ long_name = error message for error handling in CCPP
+ units = none
+ dimensions = ()
+ type = character
+ kind = len=*
+ intent = out
+[errflg]
+ standard_name = ccpp_error_code
+ long_name = error code for error handling in CCPP
+ units = 1
+ dimensions = ()
+ type = integer
+ intent = out
########################################################################
[ccpp-arg-table]