initibm Subroutine

   subroutine initibm
     use modglobal, only : libm, xh, xf, yh, yf, zh, zf, xhat, yhat, zhat, vec0, &
                           ib, ie, ih, ihc, jb, je, jh, jhc, kb, ke, kh, khc, nsv, &
                           iwallmom, lmoist, ltempeq, cexpnr, nfcts, lwritefac
     use decomp_2d, only : exchange_halo_z
     use modmpi,    only : myid
     use modstat_nc,only: open_nc, define_nc, ncinfo, writestat_dims_nc

     real, allocatable :: rhs(:,:,:)

     if (.not. libm) return

     solid_info_u%nsolpts = nsolpts_u
     solid_info_v%nsolpts = nsolpts_v
     solid_info_w%nsolpts = nsolpts_w
     call initibmnorm('solid_u.txt', solid_info_u)
     call initibmnorm('solid_v.txt', solid_info_v)
     call initibmnorm('solid_w.txt', solid_info_w)

     ! Define (real) masks
     ! Hopefully this can be removed eventually if (integer) IIx halos can be communicated
     ! These are only used in modibm, to cancel subgrid term across solid boundaries
     allocate(mask_u(ib-ih:ie+ih,jb-jh:je+jh,kb-kh:ke+kh)); mask_u = 1.
     allocate(mask_v(ib-ih:ie+ih,jb-jh:je+jh,kb-kh:ke+kh)); mask_v = 1.
     allocate(mask_w(ib-ih:ie+ih,jb-jh:je+jh,kb-kh:ke+kh)); mask_w = 1.
     mask_w(:,:,kb) = 0.     ! In future this shouldn't be needed?
     mask_u(:,:,kb-kh) = 0.
     mask_v(:,:,kb-kh) = 0.
     mask_w(:,:,kb-kh) = 0.

     allocate(rhs(ib-ih:ie+ih,jb-jh:je+jh,kb:ke+kh))
     call solid(solid_info_u, mask_u, rhs, 0., ih, jh, kh)
     call solid(solid_info_v, mask_v, rhs, 0., ih, jh, kh)
     call solid(solid_info_w, mask_w, rhs, 0., ih, jh, kh)
     call exchange_halo_z(mask_u)!, opt_zlevel=(/ih,jh,0/))
     call exchange_halo_z(mask_v)!, opt_zlevel=(/ih,jh,0/))
     call exchange_halo_z(mask_w)!, opt_zlevel=(/ih,jh,0/))

     if (iwallmom > 1) then
       bound_info_u%nbndpts = nbndpts_u
       bound_info_v%nbndpts = nbndpts_v
       bound_info_w%nbndpts = nbndpts_w
       bound_info_u%nfctsecs = nfctsecs_u
       bound_info_v%nfctsecs = nfctsecs_v
       bound_info_w%nfctsecs = nfctsecs_w
       call initibmwallfun('fluid_boundary_u.txt', 'facet_sections_u.txt', xhat, bound_info_u)
       call initibmwallfun('fluid_boundary_v.txt', 'facet_sections_v.txt', yhat, bound_info_v)
       call initibmwallfun('fluid_boundary_w.txt', 'facet_sections_w.txt', zhat, bound_info_w)
     end if

     if (ltempeq .or. lmoist .or. nsv>0 .or. lwritefac) then
       solid_info_c%nsolpts = nsolpts_c
       call initibmnorm('solid_c.txt', solid_info_c)

       bound_info_c%nbndpts = nbndpts_c
       bound_info_c%nfctsecs = nfctsecs_c
       call initibmwallfun('fluid_boundary_c.txt', 'facet_sections_c.txt', vec0, bound_info_c)

       allocate(mask_c(ib-ih:ie+ih,jb-jh:je+jh,kb-kh:ke+kh)); mask_c = 1.
       mask_c(:,:,kb-kh) = 0.
       call solid(solid_info_c, mask_c, rhs, 0., ih, jh, kh)
       call exchange_halo_z(mask_c)!, opt_zlevel=(/ih,jh,0/))
     end if

     deallocate(rhs)

     ! write facet stresses and pressure to fac.xxx.nc
     if (lwritefac) then
       allocate(fac_tau_x(1:nfcts))
       allocate(fac_tau_y(1:nfcts))
       allocate(fac_tau_z(1:nfcts))
       allocate(fac_pres(1:nfcts))
       allocate(fac_htc(1:nfcts))
       allocate(fac_cth(1:nfcts))
       fac_tau_x = 0.
       fac_tau_y = 0.
       fac_tau_z = 0.
       fac_pres = 0.
       fac_htc = 0.
       fac_cth = 0.
       allocate(fac_tau_x_av(1:nfcts))
       allocate(fac_tau_y_av(1:nfcts))
       allocate(fac_tau_z_av(1:nfcts))
       allocate(fac_pres_av(1:nfcts))
       allocate(fac_htc_av(1:nfcts))
       allocate(fac_cth_av(1:nfcts))
       fac_tau_x_av = 0.
       fac_tau_y_av = 0.
       fac_tau_z_av = 0.
       fac_pres_av = 0.
       fac_htc_av = 0.
       fac_cth_av = 0.

       facname(5:7) = cexpnr
       allocate(ncstatfac(nstatfac,4))
       call ncinfo(tncstatfac(1,:),'t', 'Time', 's', 'time')
       call ncinfo(ncstatfac( 1,:),'tau_x', 'tau_x', 'Pa','ft')
       call ncinfo(ncstatfac( 2,:),'tau_y', 'tau_y', 'Pa','ft')
       call ncinfo(ncstatfac( 3,:),'tau_z', 'tau_z', 'Pa','ft')
       call ncinfo(ncstatfac( 4,:),'pres', 'pressure', 'Pa','ft')
       call ncinfo(ncstatfac( 5,:),'htc', 'heat transfer coefficient', '','ft')
       call ncinfo(ncstatfac( 6,:),'cth', 'heat transfer coefficient (Ivo)', '','ft')

       if (myid==0) then
         call open_nc(facname, ncidfac, nrecfac, nfcts=nfcts)
         if (nrecfac==0) then
           call define_nc(ncidfac, 1, tncstatfac)
           call writestat_dims_nc(ncidfac)
         end if
         call define_nc(ncidfac, nstatfac, ncstatfac)
       end if
     end if

   end subroutine initibm