<commit_before>!> ----------------------------------------------------------------------------
!!  Model Initialization includes allocating memory for boundary and domain 
!!      data structures.  It reads all of the options from the namelist
!!      file (or files).  It also reads in Lat/Lon and Terrain data.  This module
!!      also sets up geographic (and vertical) look uptables for the forcing data
!!      Finally, there is a driver routine to initialize all model physics packages
!! 
!!   The module has been updated to allow arbitrary named variables
!!       this allows the use of e.g. ERAi, but still is not as flexible as it could be
!!
!!   The use of various python wrapper scripts in helpers/ makes it easy to add new
!!       datasets, and make them conform to the expectations of the current system. 
!!      For now there are no plans to near term plans to substantially modify this. 
!!
!!  Author: <NAME> (<EMAIL>)
!!
!! ----------------------------------------------------------------------------
module init
    use data_structures
    use io_routines,                only : io_read2d, io_read2di, io_read3d, &
                                           io_write3d,io_write3di
    use geo,                        only : geo_LUT, geo_interp, geo_interp2d
    use vertical_interpolation,     only : vLUT
    use microphysics,               only : mp_init
    use convection,                 only : init_convection
    use planetary_boundary_layer,   only : pbl_init
    use radiation,                  only : radiation_init
    use land_surface,               only : lsm_init
    use wind,                       only : init_winds
    use initialize_options,         only : init_options
    
    implicit none
    private
    public::init_model,init_physics
    
contains
    subroutine init_model(options,domain,boundary)
        implicit none
        type(options_type), intent(inout) :: options
        type(domain_type), intent(inout):: domain
        type(bc_type), intent(inout):: boundary
        
!       read in options file
        write(*,*) "Initializing Options"
        call init_options(options)

!       allocate and initialize the domain
        write(*,*) "Initializing Domain"
        call init_domain(options,domain)
!       allocate and initialize the boundary conditions structure (includes 3D grids too...)
!       this might be more apropriately though of as a forcing data structure (for low res model)
        write(*,*) "Initializing Boundaries"
        call init_bc(options,domain,boundary)
!         write(*,*) ""
        write(*,'(/ A)') "Finished basic initialization"
        write(*,'(A /)') "---------------------------------------"
        
    end subroutine init_model
    
    subroutine init_physics(options,domain)
        implicit none
        type(options_type), intent(in) :: options
        type(domain_type), intent(inout) :: domain

        ! initialize microphysics code (e.g. compute look up tables in Thompson et al)
        call mp_init(options) !this could easily be moved to init_model...
        
        call init_convection(domain,options)
        
        call pbl_init(domain,options)
        
        call radiation_init(domain,options)
        
        call lsm_init(domain,options)
        
    end subroutine init_physics
    
!------------------------------------------------------------------------------------------
!-==== Model Domain Section ====
!
! Begining of section focused on allocating and initializeing the model domain data structures
!
! ------------------------------------------------------------------------------------------

!   Allow running over a sub-domain, by removing the outer N grid cells from all sides of the domain (lat,lon,terrain)
    subroutine remove_edges(domain,edgesize)
        implicit none
        type(domain_type), intent(inout) :: domain
        integer, intent(in)::edgesize
        integer::nx1,ny1,nx2,ny2,nz
        real,allocatable,dimension(:,:)::temp_data
        
        nx1=size(domain%lat,1)
        ny1=size(domain%lat,2)
        nx2=nx1-(edgesize*2)
        ny2=ny1-(edgesize*2)
        
        allocate(temp_data(nx1,ny1))
        
        temp_data=domain%lat
        deallocate(domain%lat)
        allocate(domain%lat(nx2,ny2))
        domain%lat=temp_data(1+edgesize:nx1-edgesize,1+edgesize:ny1-edgesize)
        
        temp_data=domain%lon
        deallocate(domain%lon)
        allocate(domain%lon(nx2,ny2))

        domain%lon=temp_data(1+edgesize:nx1-edgesize,1+edgesize:ny1-edgesize)

        temp_data=domain%terrain
        deallocate(domain%terrain)
        allocate(domain%terrain(nx2,ny2))
        domain%terrain=temp_data(1+edgesize:nx1-edgesize,1+edgesize:ny1-edgesize)

        temp_data=domain%landmask
        deallocate(domain%landmask)
        allocate(domain%landmask(nx2,ny2))
        domain%landmask=temp_data(1+edgesize:nx1-edgesize,1+edgesize:ny1-edgesize)

        
        deallocate(temp_data)
        
        nx1=size(domain%u_geo%lat,1)
        ny1=size(domain%u_geo%lat,2)
        nx2=nx1-(edgesize*2)
        ny2=ny1-(edgesize*2)
        
        allocate(temp_data(nx1,ny1))
        temp_data=domain%u_geo%lat
        deallocate(domain%u_geo%lat)
        allocate(domain%u_geo%lat(nx2,ny2))
        domain%u_geo%lat=temp_data(1+edgesize:nx1-edgesize,1+edgesize:ny1-edgesize)
        
        temp_data=domain%u_geo%lon
        deallocate(domain%u_geo%lon)
        allocate(domain%u_geo%lon(nx2,ny2))
        domain%u_geo%lon=temp_data(1+edgesize:nx1-edgesize,1+edgesize:ny1-edgesize)

        deallocate(temp_data)
        
        nx1=size(domain%v_geo%lat,1)
        ny1=size(domain%v_geo%lat,2)
        nx2=nx1-(edgesize*2)
        ny2=ny1-(edgesize*2)
        
        allocate(temp_data(nx1,ny1))
        temp_data=domain%v_geo%lat
        deallocate(domain%v_geo%lat)
        allocate(domain%v_geo%lat(nx2,ny2))
        domain%v_geo%lat=temp_data(1+edgesize:nx1-edgesize,1+edgesize:ny1-edgesize)
        
        temp_data=domain%v_geo%lon
        deallocate(domain%v_geo%lon)
        allocate(domain%v_geo%lon(nx2,ny2))
        domain%v_geo%lon=temp_data(1+edgesize:nx1-edgesize,1+edgesize:ny1-edgesize)
        
        deallocate(temp_data)
    end subroutine remove_edges
    
!   allocate all arrays in domain
    subroutine domain_allocation(domain,nx,nz,ny,nsoil)
        implicit none
        type(domain_type), intent(inout) :: domain
        integer,intent(in)::nx,nz,ny
        integer,intent(in),optional :: nsoil
        integer :: ns
        
        ns=4
        if (present(nsoil)) then
            ns=nsoil
        endif
        
        ! atmosphere allocation
        allocate(domain%p(nx,nz,ny))        ! air pressure [Pa]
        domain%p=100000
        allocate(domain%u(nx+1,nz,ny))      ! eastward wind [m/s]
        domain%u=0
        allocate(domain%v(nx,nz,ny+1))      ! northward wind [m/s]
        domain%v=0
        allocate(domain%w(nx,nz,ny))        ! vertical wind [grid/s]
        domain%w=0
        allocate(domain%ur(nx+1,nz,ny))     ! eastward wind * density [m/s kg/m^3]
        domain%ur=0
        allocate(domain%vr(nx,nz,ny+1))     ! northward wind * density[m/s kg/m^3]
        domain%vr=0
        allocate(domain%wr(nx,nz,ny))       ! vertical wind * density [grid/s kg/m^3]
        domain%wr=0
        allocate(domain%th(nx,nz,ny))       ! potential temperature [K]
        domain%th=280
        allocate(domain%qv(nx,nz,ny))       ! water vapor [kg/kg]
        domain%qv=0.0002
        allocate(domain%cloud(nx,nz,ny))    ! liquid cloud water content mixing ratio [kg/kg]
        domain%cloud=0
        allocate(domain%ice(nx,nz,ny))      ! frozen cloud water content mixing ratio [kg/kg]
        domain%ice=0
        allocate(domain%nice(nx,nz,ny))     ! cloud ice number concentration [cm-3]
        domain%nice=0
        allocate(domain%qrain(nx,nz,ny))    ! rain mixing ratio [kg/kg]
        domain%qrain=0
        allocate(domain%nrain(nx,nz,ny))    ! rain drop number concentration [cm-3]
        domain%nrain=0
        allocate(domain%qsnow(nx,nz,ny))    ! snow  mixing ratio [kg/kg]
        domain%qsnow=0
        allocate(domain%qgrau(nx,nz,ny))    ! graupel mixing ratio [kg/kg]
        domain%qgrau=0
        allocate(domain%pii(nx,nz,ny))      ! exner function
        domain%pii=1
        allocate(domain%qv_adv_tendency(nx,nz,ny)) ! advective qv tendency [qv/s]
        domain%qv_adv_tendency=0
        allocate(domain%qv_pbl_tendency(nx,nz,ny)) ! qv tendency from PBL scheme [qv/s]
        domain%qv_pbl_tendency=0
        allocate(domain%rho(nx,nz,ny))      ! air density [kg/m^3]
        domain%rho=1
        allocate(domain%cloudfrac(nx,ny))   ! cloud fraction
        domain%cloudfrac=0
        
        ! land-atm flux allocation
        allocate(domain%rain(nx,ny))        ! accumulated total rainfall [kg/m^2]
        domain%rain=0
        allocate(domain%crain(nx,ny))       ! accumulated convective rainfall
        domain%crain=0
        allocate(domain%snow(nx,ny))        ! accumulated snow fall
        domain%snow=0
        allocate(domain%graupel(nx,ny))     ! accumulated graupel fall
        domain%graupel=0
        allocate(domain%current_rain(nx,ny))! rain fall in current time step
        domain%current_rain=0
        allocate(domain%current_snow(nx,ny))! snow fall in current time step
        domain%current_snow=0
        allocate(domain%swdown(nx,ny))      ! shortwave down at surface
        domain%swdown=0
        allocate(domain%lwdown(nx,ny))      ! longwave down at surface
        domain%lwdown=0
        allocate(domain%lwup(nx,ny))        ! longwave up from surface
        domain%lwup=0

        allocate(domain%sensible_heat(nx,ny)) ! sensible heat flux from surface
        domain%sensible_heat=0
        allocate(domain%latent_heat(nx,ny)) ! latent heat flux from surface
        domain%latent_heat=0
        allocate(domain%ground_heat(nx,ny)) ! ground heat flux into ground
        domain%ground_heat=0
        allocate(domain%pbl_height(nx,ny))  ! planetary boundary layer height (not always used)
        domain%pbl_height=0
        
        ! land surface allocation
        allocate(domain%soil_t(nx,ns,ny))   ! 3D soil temperature
        domain%soil_t=280
        allocate(domain%soil_vwc(nx,ns,ny)) ! 3D soil volumetric water content
        domain%soil_vwc=0.25
        
        allocate(domain%soil_tdeep(nx,ny))      ! deep soil temperature
        domain%soil_tdeep=280
        allocate(domain%soil_totalmoisture(nx,ny))  ! soil column total moisture content
        domain%soil_totalmoisture=500 ! =2000mm * 0.25 (vwc)
        allocate(domain%skin_t(nx,ny))          ! skin temperature
        domain%skin_t=280
        allocate(domain%snow_swe(nx,ny))        ! snow water equivalent
        domain%snow_swe=0
        allocate(domain%vegfrac(nx,ny))         ! vegetation cover fraction (%)
        domain%vegfrac=50 !% veg cover
        allocate(domain%canopy_water(nx,ny))    ! canopy water content
        domain%canopy_water=0
        allocate(domain%soil_type(nx,ny))       ! USGS soil type
        domain%soil_type=6 ! Loam
        allocate(domain%veg_type(nx,ny))        ! Vegetation type
        domain%veg_type=7  ! grassland

    end subroutine domain_allocation
    
    subroutine copy_z(input,output,interpolate_dim)
        implicit none
        class(interpolable_type), intent(in) :: input
        class(interpolable_type), intent(inout) :: output
        integer,intent(in)::interpolate_dim
        
        integer::nxi,nyi,nzi,nxo,nyo,nzo
        
        ! dimensions of the input data
        nxi=size(input%z,1)
        nzi=size(input%z,2)
        nyi=size(input%z,3)
        ! dimensions of the output data
        nxo=size(output%lat,1)
        nzo=nzi
        nyo=size(output%lat,2)

        if (allocated(output%z)) then
            deallocate(output%z)
        endif
        allocate(output%z(nxo,nzo,nyo))
        if (interpolate_dim==1) then
            output%z(2:nxo-1,:,:)=(input%z(1:nxi-1,:,:) + input%z(2:nxi,:,:))/2
            output%z(1,:,:)=input%z(1,:,:)
            output%z(nxo,:,:)=input%z(nxi,:,:)
        else if (interpolate_dim==3) then
            output%z(:,:,2:nyo-1)=(input%z(:,:,1:nyi-1) + input%z(:,:,2:nyi))/2
            output%z(:,:,1)=input%z(:,:,1)
            output%z(:,:,nyo)=input%z(:,:,nyi)
        else
            write(*,*) "Can not interpolate z data over z dimension"
        endif
        
    end subroutine copy_z

    subroutine init_domain_land(domain,options)
        implicit none
        type(options_type), intent(in) :: options
        type(domain_type), intent(inout):: domain<commit_msg>added code to copy_z in init to check for errors in domain dimensions also added a comment at the top of copy_z<commit_after>!> ----------------------------------------------------------------------------
!!  Model Initialization includes allocating memory for boundary and domain 
!!      data structures.  It reads all of the options from the namelist
!!      file (or files).  It also reads in Lat/Lon and Terrain data.  This module
!!      also sets up geographic (and vertical) look uptables for the forcing data
!!      Finally, there is a driver routine to initialize all model physics packages
!! 
!!   The module has been updated to allow arbitrary named variables
!!       this allows the use of e.g. ERAi, but still is not as flexible as it could be
!!
!!   The use of various python wrapper scripts in helpers/ makes it easy to add new
!!       datasets, and make them conform to the expectations of the current system. 
!!      For now there are no plans to near term plans to substantially modify this. 
!!
!!  Author: <NAME> (<EMAIL>)
!!
!! ----------------------------------------------------------------------------
module init
    use data_structures
    use io_routines,                only : io_read2d, io_read2di, io_read3d, &
                                           io_write3d,io_write3di
    use geo,                        only : geo_LUT, geo_interp, geo_interp2d
    use vertical_interpolation,     only : vLUT
    use microphysics,               only : mp_init
    use convection,                 only : init_convection
    use planetary_boundary_layer,   only : pbl_init
    use radiation,                  only : radiation_init
    use land_surface,               only : lsm_init
    use wind,                       only : init_winds
    use initialize_options,         only : init_options
    
    implicit none
    private
    public::init_model,init_physics
    
contains
    subroutine init_model(options,domain,boundary)
        implicit none
        type(options_type), intent(inout) :: options
        type(domain_type), intent(inout):: domain
        type(bc_type), intent(inout):: boundary
        
!       read in options file
        write(*,*) "Initializing Options"
        call init_options(options)

!       allocate and initialize the domain
        write(*,*) "Initializing Domain"
        call init_domain(options,domain)
!       allocate and initialize the boundary conditions structure (includes 3D grids too...)
!       this might be more apropriately though of as a forcing data structure (for low res model)
        write(*,*) "Initializing Boundaries"
        call init_bc(options,domain,boundary)
!         write(*,*) ""
        write(*,'(/ A)') "Finished basic initialization"
        write(*,'(A /)') "---------------------------------------"
        
    end subroutine init_model
    
    subroutine init_physics(options,domain)
        implicit none
        type(options_type), intent(in) :: options
        type(domain_type), intent(inout) :: domain

        ! initialize microphysics code (e.g. compute look up tables in Thompson et al)
        call mp_init(options) !this could easily be moved to init_model...
        
        call init_convection(domain,options)
        
        call pbl_init(domain,options)
        
        call radiation_init(domain,options)
        
        call lsm_init(domain,options)
        
    end subroutine init_physics
    
!------------------------------------------------------------------------------------------
!-==== Model Domain Section ====
!
! Begining of section focused on allocating and initializeing the model domain data structures
!
! ------------------------------------------------------------------------------------------

!   Allow running over a sub-domain, by removing the outer N grid cells from all sides of the domain (lat,lon,terrain)
    subroutine remove_edges(domain,edgesize)
        implicit none
        type(domain_type), intent(inout) :: domain
        integer, intent(in)::edgesize
        integer::nx1,ny1,nx2,ny2,nz
        real,allocatable,dimension(:,:)::temp_data
        
        nx1=size(domain%lat,1)
        ny1=size(domain%lat,2)
        nx2=nx1-(edgesize*2)
        ny2=ny1-(edgesize*2)
        
        allocate(temp_data(nx1,ny1))
        
        temp_data=domain%lat
        deallocate(domain%lat)
        allocate(domain%lat(nx2,ny2))
        domain%lat=temp_data(1+edgesize:nx1-edgesize,1+edgesize:ny1-edgesize)
        
        temp_data=domain%lon
        deallocate(domain%lon)
        allocate(domain%lon(nx2,ny2))

        domain%lon=temp_data(1+edgesize:nx1-edgesize,1+edgesize:ny1-edgesize)

        temp_data=domain%terrain
        deallocate(domain%terrain)
        allocate(domain%terrain(nx2,ny2))
        domain%terrain=temp_data(1+edgesize:nx1-edgesize,1+edgesize:ny1-edgesize)

        temp_data=domain%landmask
        deallocate(domain%landmask)
        allocate(domain%landmask(nx2,ny2))
        domain%landmask=temp_data(1+edgesize:nx1-edgesize,1+edgesize:ny1-edgesize)

        
        deallocate(temp_data)
        
        nx1=size(domain%u_geo%lat,1)
        ny1=size(domain%u_geo%lat,2)
        nx2=nx1-(edgesize*2)
        ny2=ny1-(edgesize*2)
        
        allocate(temp_data(nx1,ny1))
        temp_data=domain%u_geo%lat
        deallocate(domain%u_geo%lat)
        allocate(domain%u_geo%lat(nx2,ny2))
        domain%u_geo%lat=temp_data(1+edgesize:nx1-edgesize,1+edgesize:ny1-edgesize)
        
        temp_data=domain%u_geo%lon
        deallocate(domain%u_geo%lon)
        allocate(domain%u_geo%lon(nx2,ny2))
        domain%u_geo%lon=temp_data(1+edgesize:nx1-edgesize,1+edgesize:ny1-edgesize)

        deallocate(temp_data)
        
        nx1=size(domain%v_geo%lat,1)
        ny1=size(domain%v_geo%lat,2)
        nx2=nx1-(edgesize*2)
        ny2=ny1-(edgesize*2)
        
        allocate(temp_data(nx1,ny1))
        temp_data=domain%v_geo%lat
        deallocate(domain%v_geo%lat)
        allocate(domain%v_geo%lat(nx2,ny2))
        domain%v_geo%lat=temp_data(1+edgesize:nx1-edgesize,1+edgesize:ny1-edgesize)
        
        temp_data=domain%v_geo%lon
        deallocate(domain%v_geo%lon)
        allocate(domain%v_geo%lon(nx2,ny2))
        domain%v_geo%lon=temp_data(1+edgesize:nx1-edgesize,1+edgesize:ny1-edgesize)
        
        deallocate(temp_data)
    end subroutine remove_edges
    
!   allocate all arrays in domain
    subroutine domain_allocation(domain,nx,nz,ny,nsoil)
        implicit none
        type(domain_type), intent(inout) :: domain
        integer,intent(in)::nx,nz,ny
        integer,intent(in),optional :: nsoil
        integer :: ns
        
        ns=4
        if (present(nsoil)) then
            ns=nsoil
        endif
        
        ! atmosphere allocation
        allocate(domain%p(nx,nz,ny))        ! air pressure [Pa]
        domain%p=100000
        allocate(domain%u(nx+1,nz,ny))      ! eastward wind [m/s]
        domain%u=0
        allocate(domain%v(nx,nz,ny+1))      ! northward wind [m/s]
        domain%v=0
        allocate(domain%w(nx,nz,ny))        ! vertical wind [grid/s]
        domain%w=0
        allocate(domain%ur(nx+1,nz,ny))     ! eastward wind * density [m/s kg/m^3]
        domain%ur=0
        allocate(domain%vr(nx,nz,ny+1))     ! northward wind * density[m/s kg/m^3]
        domain%vr=0
        allocate(domain%wr(nx,nz,ny))       ! vertical wind * density [grid/s kg/m^3]
        domain%wr=0
        allocate(domain%th(nx,nz,ny))       ! potential temperature [K]
        domain%th=280
        allocate(domain%qv(nx,nz,ny))       ! water vapor [kg/kg]
        domain%qv=0.0002
        allocate(domain%cloud(nx,nz,ny))    ! liquid cloud water content mixing ratio [kg/kg]
        domain%cloud=0
        allocate(domain%ice(nx,nz,ny))      ! frozen cloud water content mixing ratio [kg/kg]
        domain%ice=0
        allocate(domain%nice(nx,nz,ny))     ! cloud ice number concentration [cm-3]
        domain%nice=0
        allocate(domain%qrain(nx,nz,ny))    ! rain mixing ratio [kg/kg]
        domain%qrain=0
        allocate(domain%nrain(nx,nz,ny))    ! rain drop number concentration [cm-3]
        domain%nrain=0
        allocate(domain%qsnow(nx,nz,ny))    ! snow  mixing ratio [kg/kg]
        domain%qsnow=0
        allocate(domain%qgrau(nx,nz,ny))    ! graupel mixing ratio [kg/kg]
        domain%qgrau=0
        allocate(domain%pii(nx,nz,ny))      ! exner function
        domain%pii=1
        allocate(domain%qv_adv_tendency(nx,nz,ny)) ! advective qv tendency [qv/s]
        domain%qv_adv_tendency=0
        allocate(domain%qv_pbl_tendency(nx,nz,ny)) ! qv tendency from PBL scheme [qv/s]
        domain%qv_pbl_tendency=0
        allocate(domain%rho(nx,nz,ny))      ! air density [kg/m^3]
        domain%rho=1
        allocate(domain%cloudfrac(nx,ny))   ! cloud fraction
        domain%cloudfrac=0
        
        ! land-atm flux allocation
        allocate(domain%rain(nx,ny))        ! accumulated total rainfall [kg/m^2]
        domain%rain=0
        allocate(domain%crain(nx,ny))       ! accumulated convective rainfall
        domain%crain=0
        allocate(domain%snow(nx,ny))        ! accumulated snow fall
        domain%snow=0
        allocate(domain%graupel(nx,ny))     ! accumulated graupel fall
        domain%graupel=0
        allocate(domain%current_rain(nx,ny))! rain fall in current time step
        domain%current_rain=0
        allocate(domain%current_snow(nx,ny))! snow fall in current time step
        domain%current_snow=0
        allocate(domain%swdown(nx,ny))      ! shortwave down at surface
        domain%swdown=0
        allocate(domain%lwdown(nx,ny))      ! longwave down at surface
        domain%lwdown=0
        allocate(domain%lwup(nx,ny))        ! longwave up from surface
        domain%lwup=0

        allocate(domain%sensible_heat(nx,ny)) ! sensible heat flux from surface
        domain%sensible_heat=0
        allocate(domain%latent_heat(nx,ny)) ! latent heat flux from surface
        domain%latent_heat=0
        allocate(domain%ground_heat(nx,ny)) ! ground heat flux into ground
        domain%ground_heat=0
        allocate(domain%pbl_height(nx,ny))  ! planetary boundary layer height (not always used)
        domain%pbl_height=0
        
        ! land surface allocation
        allocate(domain%soil_t(nx,ns,ny))   ! 3D soil temperature
        domain%soil_t=280
        allocate(domain%soil_vwc(nx,ns,ny)) ! 3D soil volumetric water content
        domain%soil_vwc=0.25
        
        allocate(domain%soil_tdeep(nx,ny))      ! deep soil temperature
        domain%soil_tdeep=280
        allocate(domain%soil_totalmoisture(nx,ny))  ! soil column total moisture content
        domain%soil_totalmoisture=500 ! =2000mm * 0.25 (vwc)
        allocate(domain%skin_t(nx,ny))          ! skin temperature
        domain%skin_t=280
        allocate(domain%snow_swe(nx,ny))        ! snow water equivalent
        domain%snow_swe=0
        allocate(domain%vegfrac(nx,ny))         ! vegetation cover fraction (%)
        domain%vegfrac=50 !% veg cover
        allocate(domain%canopy_water(nx,ny))    ! canopy water content
        domain%canopy_water=0
        allocate(domain%soil_type(nx,ny))       ! USGS soil type
        domain%soil_type=6 ! Loam
        allocate(domain%veg_type(nx,ny))        ! Vegetation type
        domain%veg_type=7  ! grassland

    end subroutine domain_allocation
    
! interpolate intput%z to output%z assuming that input has one less grid cell 
! in the interpolate_dim dimension
    subroutine copy_z(input,output,interpolate_dim)
        
        implicit none
        class(interpolable_type), intent(in) :: input
        class(interpolable_type), intent(inout) :: output
        integer,intent(in)::interpolate_dim
        
        integer::nxi,nyi,nzi,nxo,nyo,nzo
        
        ! dimensions of the input data
        nxi=size(input%z,1)
        nzi=size(input%z,2)
        nyi=size(input%z,3)
        ! dimensions of the output data
        nxo=size(output%lat,1)
        nzo=nzi
        nyo=size(output%lat,2)

        if (allocated(output%z)) then
            deallocate(output%z)
        endif
        allocate(output%z(nxo,nzo,nyo))
        if (interpolate_dim==1) then
            if ((nxo/=(nxi+1)).or.(nyo/=nyi).or.(nzo/=nzi)) then
                stop("Error copying z levels from mass grid to U grid, dimensions don't match")
            endif
            output%z(2:nxo-1,:,:)=(input%z(1:nxi-1,:,:) + input%z(2:nxi,:,:))/2
            output%z(1,:,:)=input%z(1,:,:)
            output%z(nxo,:,:)=input%z(nxi,:,:)
        else if (interpolate_dim==3) then
            if ((nyo/=(nyi+1)).or.(nxo/=nxi).or.(nzo/=nzi)) then
                stop("Error copying z levels from mass grid to V grid, dimensions don't match")
            endif
            output%z(:,:,2:nyo-1)=(input%z(:,:,1:nyi-1) + input%z(:,:,2:nyi))/2
            output%z(:,:,1)=input%z(:,:,1)
            output%z(:,:,nyo)=input%z(:,:,nyi)
        else
            write(*,*) "Can not interpolate z data over z dimension"
        endif
        
    end subroutine copy_z

    subroutine init_domain_land(domain,options)
        implicit none
        type(options_type), intent(in) :: options
        type(domain_type), intent(inout):: domain
{'fraction_non_alphanumeric': 0.13346713497240342, 'fraction_numerical': 0.015901810181790113, 'mean_word_length': 2.28473630831643, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4958909', 'n_tokens_mistral': 8384, 'n_tokens_neox': 8081, 'n_words': 1622}