init_vegetation Subroutine

public subroutine init_vegetation()

Uses

  • proc~~init_vegetation~~UsesGraph proc~init_vegetation init_vegetation decomp_2d decomp_2d proc~init_vegetation->decomp_2d module~modglobal modglobal proc~init_vegetation->module~modglobal module~modmpi modmpi proc~init_vegetation->module~modmpi module~readinput readinput proc~init_vegetation->module~readinput mpi mpi module~modmpi->mpi module~readinput->decomp_2d module~readinput->module~modglobal module~readinput->module~modmpi module~readinput->mpi

Arguments

None

Calls

proc~~init_vegetation~~CallsGraph proc~init_vegetation init_vegetation exchange_halo_z exchange_halo_z proc~init_vegetation->exchange_halo_z mpi_barrier mpi_barrier proc~init_vegetation->mpi_barrier mpi_bcast mpi_bcast proc~init_vegetation->mpi_bcast proc~init_vegetation_legacy init_vegetation_legacy proc~init_vegetation->proc~init_vegetation_legacy proc~read_sparse_ijk read_sparse_ijk proc~init_vegetation->proc~read_sparse_ijk proc~read_sparse_real read_sparse_real proc~init_vegetation->proc~read_sparse_real proc~read_sparse_ijk->mpi_bcast zend zend proc~read_sparse_ijk->zend zstart zstart proc~read_sparse_ijk->zstart proc~read_sparse_real->mpi_bcast

Called by

proc~~init_vegetation~~CalledByGraph proc~init_vegetation init_vegetation program~udales uDALES program~udales->proc~init_vegetation

Source Code

  subroutine init_vegetation
    use modglobal,  only : ltrees,ltreedump,itree_mode,TREE_MODE_DRAG_ONLY,TREE_MODE_SVEG,TREE_MODE_LEGACY_SEB,ib,ie,jb,je,kb,ke,ih,jh,kh,cexpnr,dzf,nsv
    use modmpi,     only : myid,comm3d,mpierr,MY_REAL
    use readinput,  only : read_sparse_ijk, read_sparse_real
    use decomp_2d,  only : exchange_halo_x, exchange_halo_y, exchange_halo_z
    implicit none
    integer :: i,j,k,m
    integer :: npts
    integer, allocatable :: ids_loc(:)
    integer, allocatable :: pts_in(:,:)
    character(200) :: filename
    integer :: ierr, ifinput
    character(256) :: line
    integer, allocatable :: id_all(:)
    real, allocatable :: lad_all(:), cd_all(:), ud_all(:), dec_all(:), lsize_all(:), rs_all(:), sveg_loc(:)
    integer :: nread
    integer :: idx
    real, allocatable :: dcoef_f(:,:,:)
    logical, allocatable :: mask_f(:,:,:)
    logical :: sveg_exists

    if (.not. ltrees) then
      ltreedump = .false.
      return
    end if
    vegetation_ready = .false.
    has_sveg = .false.

    ! count points in veg.inp.<expnr> to set npts
    ! TEMPORARY WHILE WE RUN TREES AND VEG TOGETHER
    ! npts is contained in ntrees namelist variable
    write(filename, '(A,A)') 'veg.inp.', trim(cexpnr)
    ifinput = 99
    npts = 0
    if (myid == 0) then
      open(ifinput, file=filename, status='old', iostat=ierr)
      if (ierr /= 0) then
        write(*, '(A,A)') 'ERROR: Cannot open file: ', trim(filename)
        stop 1
      end if
      read(ifinput, '(a256)', iostat=ierr) line  ! skip header
      do
        read(ifinput, *, iostat=ierr) i, j, k
        if (ierr /= 0) exit
        npts = npts + 1
      end do
      close(ifinput)
    end if
    call MPI_BCAST(npts, 1, MPI_INTEGER, 0, comm3d, mpierr)

    ! Read veg_params.<expnr> (aligned with veg.inp ordering)
    if (allocated(id_all)) deallocate(id_all, lad_all, cd_all, ud_all, dec_all, lsize_all, rs_all)
    allocate(id_all(npts))
    allocate(lad_all(npts))
    allocate(cd_all(npts))
    allocate(ud_all(npts))
    allocate(dec_all(npts))
    allocate(lsize_all(npts))
    allocate(rs_all(npts))
    if (myid == 0) then
      write(filename, '(A,A)') 'veg_params.inp.', trim(cexpnr)
      open(ifinput, file=filename, status='old', iostat=ierr)
      if (ierr /= 0) then
        write(*, '(A,A)') 'ERROR: Cannot open file: ', trim(filename)
        stop 1
      end if
      read(ifinput, '(a256)', iostat=ierr) line  ! skip header
      do nread = 1, npts
        read(ifinput, *, iostat=ierr) id_all(nread), lad_all(nread), cd_all(nread), ud_all(nread), dec_all(nread), lsize_all(nread), rs_all(nread)
        if (ierr /= 0) then
          write(*, '(A,I0)') 'ERROR reading veg_params line ', nread
          stop 1
        end if
      end do
      close(ifinput)

      write(filename, '(A,A)') 'sveg.inp.', trim(cexpnr)
      inquire(file=filename, exist=sveg_exists)
      if (itree_mode /= TREE_MODE_SVEG) then
        if (sveg_exists) then
          write(*,'(A,A)') 'NOTE: Found optional vegetation shortwave file: ', trim(filename)
        end if
      else
        if (.not. sveg_exists) then
          write(*,'(A,A)') 'ERROR: Missing vegetation shortwave file: ', trim(filename)
          stop 1
        end if
      end if

    end if

    call MPI_BCAST(id_all, npts, MPI_INTEGER, 0, comm3d, mpierr)
    call MPI_BCAST(lad_all, npts, MY_REAL, 0, comm3d, mpierr)
    call MPI_BCAST(cd_all, npts, MY_REAL, 0, comm3d, mpierr)
    call MPI_BCAST(ud_all, npts, MY_REAL, 0, comm3d, mpierr)
    call MPI_BCAST(dec_all, npts, MY_REAL, 0, comm3d, mpierr)
    call MPI_BCAST(lsize_all, npts, MY_REAL, 0, comm3d, mpierr)
    call MPI_BCAST(rs_all, npts, MY_REAL, 0, comm3d, mpierr)
    call MPI_BCAST(sveg_exists, 1, MPI_LOGICAL, 0, comm3d, mpierr)

    call read_sparse_ijk('veg.inp.'//trim(cexpnr), npts, veg%npts, ids_loc, pts_in, nskip=1)
    if (sveg_exists) then
      call read_sparse_real('sveg.inp.'//trim(cexpnr), npts, ids_loc, sveg_loc, nskip=1)
    end if

    if (allocated(veg%id)) then
      deallocate(veg%id, veg%gidx, veg%ijk, veg%lad, veg%cd, veg%ud, veg%dec, veg%lsize, veg%rs, veg%laiv)
    end if
    if (allocated(sveg)) deallocate(sveg)
    if (allocated(vegp%qt)) deallocate(vegp%qt, vegp%qtR, vegp%qtA, vegp%thl, vegp%omega, vegp%sv)

    allocate(veg%id(veg%npts))
    allocate(veg%gidx(veg%npts))
    allocate(veg%ijk(veg%npts,3))
    allocate(veg%lad(veg%npts))
    allocate(veg%cd(veg%npts))
    allocate(veg%ud(veg%npts))
    allocate(veg%dec(veg%npts))
    allocate(veg%lsize(veg%npts))
    allocate(veg%rs(veg%npts))
    allocate(veg%laiv(veg%npts))
    allocate(sveg(veg%npts))
    allocate(vegp%qt(veg%npts))
    allocate(vegp%qtR(veg%npts))
    allocate(vegp%qtA(veg%npts))
    allocate(vegp%thl(veg%npts))
    allocate(vegp%omega(veg%npts))
    allocate(vegp%sv(veg%npts,max(1,nsv)))

    veg%ijk = 0
    veg%laiv = 0.
    sveg = 0.
    vegp%qt = 0.
    vegp%qtR = 0.
    vegp%qtA = 0.
    vegp%thl = 0.
    vegp%omega = 0.
    vegp%sv = 0.
    has_sveg = sveg_exists

    do m = 1, veg%npts
      veg%id(m)    = id_all(ids_loc(m))
      veg%gidx(m)  = ids_loc(m)
      veg%lad(m)   = lad_all(ids_loc(m))
      veg%cd(m)    = cd_all(ids_loc(m))
      veg%ud(m)    = ud_all(ids_loc(m))
      veg%dec(m)   = dec_all(ids_loc(m))
      veg%lsize(m) = lsize_all(ids_loc(m))
      veg%rs(m)    = rs_all(ids_loc(m))
      veg%ijk(m,1:3) = pts_in(m,1:3)
      if (has_sveg) sveg(m) = sveg_loc(m)
    end do

    call MPI_BARRIER(comm3d, mpierr)

    deallocate(ids_loc)
    deallocate(pts_in)
    if (allocated(sveg_loc)) deallocate(sveg_loc)
    deallocate(id_all, lad_all, cd_all, ud_all, dec_all, lsize_all, rs_all)

    if (allocated(lad_3d)) deallocate(lad_3d)
    if (allocated(dcoef_3d)) deallocate(dcoef_3d)
    allocate(lad_3d(ib-ih:ie+ih, jb-jh:je+jh, kb-kh:ke+kh))
    allocate(dcoef_3d(ib-ih:ie+ih, jb-jh:je+jh, kb-kh:ke+kh))
    lad_3d = 0.
    dcoef_3d = 0.

    ! Build 3D LAD and drag parameter fields from sparse vegetation points.
    do m = 1, veg%npts
      i = veg%ijk(m,1)
      j = veg%ijk(m,2)
      k = veg%ijk(m,3)
      lad_3d(i, j, k) = veg%lad(m)
      dcoef_3d(i, j, k) = veg%lad(m) * veg%cd(m)
    end do

    ! Exchange halos so face averaging has neighbor values (2D decomposition)
    call exchange_halo_z(lad_3d)
    call exchange_halo_z(dcoef_3d)

    if (itree_mode == TREE_MODE_LEGACY_SEB) then
      call init_vegetation_legacy
    end if

    ! ========================================================================
    ! Precompute face lists for staggered u/v/w points 
    ! ========================================================================
    if (allocated(ijk_u)) deallocate(ijk_u, dcoef_u, veg_up)
    if (allocated(ijk_v)) deallocate(ijk_v, dcoef_v, veg_vp)
    if (allocated(ijk_w)) deallocate(ijk_w, dcoef_w, veg_wp)
    npts_u = 0
    npts_v = 0
    npts_w = 0

    allocate(dcoef_f(ib:ie, jb:je, kb:ke))
    allocate(mask_f(ib:ie, jb:je, kb:ke))

    ! u-faces (i-1/i)
    dcoef_f = 0.5*(dcoef_3d(ib-1:ie-1, jb:je, kb:ke) + dcoef_3d(ib:ie, jb:je, kb:ke))
    mask_f = (dcoef_f > 0.0)

    npts_u = count(mask_f)
    if (npts_u > 0) then
      allocate(ijk_u(npts_u,3), dcoef_u(npts_u), veg_up(npts_u))
      idx = 0
      do k = kb, ke
        do j = jb, je
          do i = ib, ie
            if (mask_f(i,j,k)) then
              idx = idx + 1
              ijk_u(idx,1:3) = (/ i, j, k /)
              dcoef_u(idx) = dcoef_f(i,j,k)
            end if
          end do
        end do
      end do
      veg_up = 0.
    end if

    ! v-faces (j-1/j)
    dcoef_f = 0.5*(dcoef_3d(ib:ie, jb-1:je-1, kb:ke) + dcoef_3d(ib:ie, jb:je, kb:ke))
    mask_f = (dcoef_f > 0.0)
    npts_v = count(mask_f)
    if (npts_v > 0) then
      allocate(ijk_v(npts_v,3), dcoef_v(npts_v), veg_vp(npts_v))
      idx = 0
      do k = kb, ke
        do j = jb, je
          do i = ib, ie
            if (mask_f(i,j,k)) then
              idx = idx + 1
              ijk_v(idx,1:3) = (/ i, j, k /)
              dcoef_v(idx) = dcoef_f(i,j,k)
            end if
          end do
        end do
      end do
      veg_vp = 0.
    end if

    ! w-faces (k-1/k); no z-halos, so start at kb+1
    dcoef_f = 0.0
    dcoef_f(:,:,kb+1:ke) = 0.5*(dcoef_3d(ib:ie, jb:je, kb:ke-1) + dcoef_3d(ib:ie, jb:je, kb+1:ke))
    mask_f = (dcoef_f > 0.0)
    npts_w = count(mask_f)
    if (npts_w > 0) then
      allocate(ijk_w(npts_w,3), dcoef_w(npts_w), veg_wp(npts_w))
      idx = 0
      do k = kb+1, ke
        do j = jb, je
          do i = ib, ie
            if (mask_f(i,j,k)) then
              idx = idx + 1
              ijk_w(idx,1:3) = (/ i, j, k /)
              dcoef_w(idx) = dcoef_f(i,j,k)
            end if
          end do
        end do
      end do
      veg_wp = 0.
    end if

    deallocate(dcoef_f, mask_f)

    vegetation_ready = .true.

  end subroutine init_vegetation