text:
! module for production of movie (visualization) files written in raw format
! and postprocessed into gif files.
! Original implementation by Peter Bogenschutz, U. Uta, April 2008
!

module movies

use grid

implicit none

  ! Variables required for production of movie files
  real cldttmp_xy(nx,ny), cldttmpl_xy(nx,ny), cwppath(nx,ny)
  real sfcflux(nx,ny), sfcpflog10dat(nx,ny)
  real dtmovie, amin, amax
  integer irecc
  logical :: files_opened=.false.

! limits defaults. Can be set in namelist MOVIES
  real :: &
  u_min=-30., u_max=30., & ! surface velocity in x
  v_min=-30. , v_max=30., & ! surface velocity in x
  cldtop_min=290., cldtop_max=310., & ! cloud-top temperature
  sst_min=285., sst_max=305., & ! surface temperature
  tasfc_min=280., tasfc_max=310., & ! surface air temperature
  qvsfc_min=0.008, qvsfc_max=0.025, & ! surface vapor mixing ratio
  prec_min=0.00000001, prec_max=0.1, & ! surface precipitation (will plot log scale)
  cwp_min=0.1, cwp_max=1000., & ! cloud water path
  iwp_min=0.1, iwp_max=1000.   ! ice water path

CONTAINS

subroutine init_movies()

  integer ierr, ios, ios_missing_namelist, place_holder

  NAMELIST /MOVIES/ &
       u_min, u_max, &
       v_min, v_max, &
       cldtop_min, cldtop_max, &
       sst_min, sst_max, &
       tasfc_min, tasfc_max, &
       qvsfc_min, qvsfc_max, &
       prec_min, prec_max, &
       cwp_min, cwp_max, &
       iwp_min, iwp_max

  NAMELIST /BNCUIODSBJCB/ place_holder
 !----------------------------------
  !  Read namelist for movie options from prm file:
  !------------
  open(55,file='./'//trim(case)//'/prm', status='old',form='formatted')

  read (UNIT=55,NML=BNCUIODSBJCB,IOSTAT=ios_missing_namelist)
  rewind(55) !note that one must rewind before searching for new namelists

  read (55,MOVIES,IOSTAT=ios)

  if (ios.ne.0) then
     !namelist error checking
     if(ios.ne.ios_missing_namelist) then
        write(*,*) '****** ERROR: bad specification in MOVIES namelist'
        call task_abort()
     end if
  end if
  close(55)

end subroutine init_movies

subroutine openmoviefiles()

use vars
use params
implicit none

character *4 rankchar
integer, external :: lenstr, bytes_in_rec

! Subroutine to open the files needed to create the movies.

if(files_opened) return

if(masterproc) open(79,file='./OUT_MOVIES/'//trim(case)//'_'//trim(caseid)//'.info.movie',form='formatted')

write(rankchar,'(i4)') rank

open(unit=80,file='./OUT_MOVIES/'//trim(case)//'_'//trim(caseid)//'_usfc.raw_'//rankchar(5-lenstr(rankchar):4), &
                                              form='unformatted',access='direct',recl=nx*ny/bytes_in_rec())

open(unit=81,file='./OUT_MOVIES/'//trim(case)//'_'//trim(caseid)//'_vsfc.raw_'//rankchar(5-lenstr(rankchar):4), &
                                              form='unformatted',access='direct',recl=nx*ny/bytes_in_rec())

open(unit=82,file='./OUT_MOVIES/'//trim(case)//'_'//trim(caseid)//'_cldtop.raw_'//rankchar(5-lenstr(rankchar):4),&
                                              form='unformatted',access='direct',recl=nx*ny/bytes_in_rec())

open(unit=83,file='./OUT_MOVIES/'//trim(case)//'_'//trim(caseid)//'_tasfc.raw_'//rankchar(5-lenstr(rankchar):4),&
                                              form='unformatted',access='direct',recl=nx*ny/bytes_in_rec())

open(unit=84,file='./OUT_MOVIES/'//trim(case)//'_'//trim(caseid)//'_qvsfc.raw_'//rankchar(5-lenstr(rankchar):4),&
                                              form='unformatted',access='direct',recl=nx*ny/bytes_in_rec())

open(unit=85,file='./OUT_MOVIES/'//trim(case)//'_'//trim(caseid)//'_sfcprec.raw_'//rankchar(5-lenstr(rankchar):4), &
                                              form='unformatted',access='direct',recl=nx*ny/bytes_in_rec())

open(unit=86,file='./OUT_MOVIES/'//trim(case)//'_'//trim(caseid)//'_cwp.raw_'//rankchar(5-lenstr(rankchar):4), &
                                              form='unformatted',access='direct',recl=nx*ny/bytes_in_rec())

open(unit=87,file='./OUT_MOVIES/'//trim(case)//'_'//trim(caseid)//'_iwp.raw_'//rankchar(5-lenstr(rankchar):4), &
                                              form='unformatted',access='direct',recl=nx*ny/bytes_in_rec())

open(unit=88,file='./OUT_MOVIES/'//trim(case)//'_'//trim(caseid)//'_sst.raw_'//rankchar(5-lenstr(rankchar):4), &
                                              form='unformatted',access='direct',recl=nx*ny/bytes_in_rec())

files_opened = .true.
return

end subroutine openmoviefiles

!----------------------------------------------------------------------

subroutine mvmovies()

use vars
use params
implicit none

call t_startf ('movies')
! Subroutine to call the subroutines required to make the movies.
! Modify amin and max as you see fit.

call openmoviefiles()

call cldtoptmp()

!  Surface U-Wind Component
amin=u_max
amax=u_min
call plotpix(u(1:nx,1:ny,1),nx,ny,nstep,irecc,amin,amax,80)

!  Surface V-Wind Component
amin=v_max
amax=v_min
call plotpix(v(1:nx,1:ny,1),nx,ny,nstep,irecc,amin,amax,81)

!  Cloud Top Temperature
amin=cldtop_max
amax=cldtop_min
call plotpix(cldttmp_xy(:,:),nx,ny,nstep,irecc,amin,amax,82)

!  Surface temperature
amin=sst_max
amax=sst_min
call plotpix(sstxy(1:nx,1:ny),nx,ny,nstep,irecc,amin,amax,88)

!  Surface air temperature
amin=tasfc_max
amax=tasfc_min
call plotpix(tabs(:,:,1),nx,ny,nstep,irecc,amin,amax,83)

!  Water vapor mixing ratio
amin=qvsfc_max
amax=qvsfc_min
call plotpix(qv(:,:,1),nx,ny,nstep,irecc,amin,amax,84)

!  Surface precip rate
call getvtrr(nx,ny,nz_gl,qpl(:,:,1),tabs(:,:,1),rho,sfcflux)
amin=alog10(prec_max)
amax=alog10(prec_min)
call convertlog10(sfcflux,nx,ny,amax,sfcpflog10dat)
call plotpix(sfcpflog10dat,nx,ny,nstep,irecc,amin,amax,85)

!  Cloud Water Path
amin=alog10(cwp_max)
amax=alog10(cwp_min)
sfcpflog10dat(:,:)=0.
call pathvar(nx,ny,nz_gl,z,qcl(:,:,:),cwppath)
call convertlog10(cwppath,nx,ny,amax,sfcpflog10dat)
call plotpix(sfcpflog10dat,nx,ny,nstep,irecc,amin,amax,86)

!  Ice Water Path
amin=alog10(iwp_max)
amax=alog10(iwp_min)
sfcpflog10dat(:,:)=0.
cwppath(:,:)=0.
call pathvar(nx,ny,nz_gl,z,qci(:,:,:)+qpi(:,:,:),cwppath)
call convertlog10(cwppath,nx,ny,amax,sfcpflog10dat)
call plotpix(sfcpflog10dat,nx,ny,nstep,irecc,amin,amax,87)

irecc = irecc+1
call close_files()

if(masterproc) print*, 'appending OUT_MOVIES/*.raw files. nstep=', nstep

call t_stopf ('movies')
return

end subroutine mvmovies

!---------------------------------------------------------------------

subroutine close_files()

! Close the movie files after last timestep of integration

!  close(unit=80)
!  close(unit=81)
!  close(unit=82)
!  close(unit=83)
!  close(unit=84)
!  close(unit=85)
!  close(unit=86)
!  close(unit=87)
!  close(unit=88)

  if(masterproc) then
!   open(79,file='./OUT_MOVIES/'//trim(case)//'_'//trim(caseid)//'.info.movie',form='formatted')
   rewind(79)
   write(79,*) nsubdomains_x,nsubdomains_y
   write(79,*) nx_gl/nsubdomains_x,ny_gl/nsubdomains_y
   write(79,*) irecc-1
   write(79,*) 'field: min, max values:'
   write(79,*) 'u:',u_min, u_max
   write(79,*) 'v:', v_min, v_max
   write(79,*) 'cldtop:', cldtop_min, cldtop_max
   write(79,*) 'sst:', sst_min, sst_max
   write(79,*) 'tasfc:', tasfc_min, tasfc_max
   write(79,*) 'qvsfc:', qvsfc_min, qvsfc_max
   write(79,*) 'log10(prec):', alog10(prec_min), alog10(prec_max)
   write(79,*) 'log10(cwp):', alog10(cwp_min), alog10(cwp_max)
   write(79,*) 'log10(iwp):', alog10(iwp_min), alog10(iwp_max)
!   close(79)
  end if

  return

end subroutine close_files

!---------------------------------------------------------------------

subroutine pixplt_raw(data,nx,ny,amin,amax,iun,irec)

  implicit none

  real amin, amax
  integer nx, ny, i, j, idata, iun, ncolor, irec, count
  parameter(ncolor = 254)
  character*1 iclrs(nx,ny)
  real data(nx,ny), rdata

  count=0
  do j=1,ny
    do i=1,nx
      rdata=data(i,j)
      idata=int((rdata-amin)/(amax-amin)*float(ncolor))+1
      idata=max(0,min(idata,ncolor+1))
      iclrs(i,j)=char(idata)
      if (idata .ge. 254) then
        count=count+1
      endif
!      write(*,*) idata, data(i,j), iclrs(i,j)
    enddo
  enddo

  write(iun,rec=irec) iclrs(:,:)

  return

end subroutine pixplt_raw

!----------------------------------------------------------------------

subroutine plotpix(data,nx,ny,n,irec,amin,amax,unitnum)

  implicit none

  integer nx, ny, n, unitnum, irec
  real data(nx,ny), amin, amax

    call pixplt_raw(data,nx,ny,amin,amax,unitnum,irec)

  return

end subroutine plotpix

!---------------------------------------------------------------------

subroutine convertlog10(data,nx,nz,minval,retdat)

  implicit none

  integer i, j, nx, nz
  real minval, data(nx,nz), retdat(nx,nz)

  do j=1,nz
    do i=1,nx
      if (data(i,j) .eq. 0 .or. alog10(data(i,j)) .lt. minval) then
        retdat(i,j)=minval
      else
        retdat(i,j)=alog10(data(i,j))
      endif
    enddo
  enddo

  return

end subroutine convertlog10

!----------------------------------------------------------------------

subroutine getvtrr(nx,ny,nz,qr,tabs,rhopro,flux)

  implicit none

  integer i, j, nx, ny, nz
  real qrho, act2, vconr, rho, vtrr, rhofac
  real qr(nx,ny), tabs(nx,ny), p(nx,ny), flux(nx,ny), rhopro(nz)

  real, parameter :: pie=3.141593
  real, parameter :: rnzr=8.e6
  real, parameter :: alin=841.99667
  real, parameter :: gam480=17.8379
  real, parameter :: rhor=1.e3

  act2=pie*rnzr*rhor
  vconr=alin*gam480/(6.*act2**(0.2))

  do j=1,ny
    do i=1,nx
      rho=rhopro(1)
      qrho=rho*qr(i,j)
      rhofac=sqrt(1.226/rho)
      vtrr=amin1(vconr*qrho**(0.2)*rhofac,10.0)
      flux(i,j)=rho*vtrr*qr(i,j)
    enddo
  enddo

  return

end subroutine getvtrr

!----------------------------------------------------------------------

subroutine pathvar(nx,ny,nz,z,var,path)

  implicit none

  integer i, j, k, nx, ny, nz
  real var(nx,ny,nz), path(nx,ny), z(nz+1), rho

  rho=1000.
  path(:,:)=0.
  do k=2,nz
    do j=1,ny
      do i=1,nx
        path(i,j)=path(i,j)+(rho*var(i,j,k))*(z(k)-z(k-1))
      enddo
    enddo
  enddo

  return

end subroutine pathvar

!----------------------------------------------------------------------

subroutine cldtoptmp()

use vars
use params
implicit none

integer chkind, chkind2, tcount, i, j, k
real coef1, downint

cldttmp_xy(:,:)=0.
cldttmpl_xy(:,:)=0.
tcount=0
do j=1,ny
  do i=1,nx

    downint=0.
    chkind=1.
    chkind2=1.
    do k=nzm,1,-1  ! Integrate downward
      coef1=rho(k)*dz*adz(k)*dtfactor
      downint=downint+(qcl(i,j,k)+qci(i,j,k))*coef1
      if (downint .ge. 0.1 .and. chkind .eq. 1) then
        cldttmp_xy(i,j)=tabs(i,j,k)
        chkind=2.
      endif
    end do

    do k=1,nzm-1
      if (qcl(i,j,k)+qci(i,j,k) .gt. 0 .and. qcl(i,j,k+1)+qci(i,j,k+1) .eq. 0 .and. chkind2 .eq. 1) then
        cldttmpl_xy(i,j)=tabs(i,j,k)
!        write(*,*) z
        chkind2=2.
      endif
    enddo

    if (chkind2 .eq. 1) then
      cldttmpl_xy(i,j)=sstxy(i,j)
    endif

    if (chkind .eq. 1) then
      cldttmp_xy(i,j)=sstxy(i,j)
      tcount=tcount+1
    endif

  end do
end do

return

end subroutine cldtoptmp

end module movies

max_stars_repo_name: kbren/uwnet

max_issues_repo_name: kbren/uwnet

max_forks_repo_name: kbren/uwnet


text:
subroutine navierstokes ( &
  fn_len, filename, &
  mesh_blk, mesh_n_blk, mesh_bnd_blk, mesh_dmn_blk, &

  mesh_n_dim, mesh_n_vrt, mesh_vrt_idx, geom_x_vrt, &

  mshf_n_elm, mshf_elm_idx, mshf_elm_ty, mshf_elm_dim, &
  mshf_elm_fce_ptr, mshf_fce_vrt_ptr, mshf_fce_vrt_lst, &

  mesh_n_bnd, mesh_bnd_idx, mesh_bnd_tag_ptr, &
  mesh_bnd_tag_lst, mesh_bnd_vrt_ptr, mesh_bnd_vrt_lst, &

  mesh_n_gbnd, mesh_gbnd_idx, mesh_gbnd_tag_ptr, &
  mesh_gbnd_tag_lst, mesh_gbnd_vrt_ptr, mesh_gbnd_vrt_lst, &

  mesh_n_dmn, mesh_dmn_idx, mesh_dmn_tag_ptr, &
  mesh_dmn_tag_lst, mesh_dmn_vrt_ptr, mesh_dmn_vrt_lst, &

  mesh_n_gdmn, mesh_gdmn_idx, mesh_gdmn_tag_ptr, &
  mesh_gdmn_tag_lst, mesh_gdmn_vrt_ptr, mesh_gdmn_vrt_lst, &

  mesh_n_xbnd, mesh_xbnd_vrt_ptr, mesh_xbnd_vrt_lst, &
  mesh_xbnd_elm_ptr, mesh_xbnd_elm_lst, &

  mesh_n_xdmn, mesh_xdmn_elm_ptr, mesh_xdmn_elm_lst, &
  mesh_xdmn_vrt_ptr, mesh_xdmn_vrt_lst, &

  mesh_bnd_neig_ptr, mesh_bnd_neig_lst, &
  mesh_dmn_neig_ptr, mesh_dmn_neig_lst &
)

!   use mpi
  use constants_m
  use alloc_m
  use solv_m
  use hash_m
  use geometry_m
  use transport_m
  use fson

  implicit none

  integer, intent(in) :: fn_len
  character(len=fn_len), intent(in) :: filename
  integer, intent(in) :: mesh_blk
  integer, intent(in) :: mesh_n_blk
  integer, dimension(mesh_n_blk + 1), intent(in), target :: mesh_bnd_blk
  integer, dimension(mesh_n_blk + 1), intent(in), target :: mesh_dmn_blk

  integer, intent(in) :: mesh_n_dim, mesh_n_vrt
  integer, dimension(mesh_n_vrt), intent(in), target :: mesh_vrt_idx
  real(dp), dimension(3, mesh_n_vrt), intent(in), target :: geom_x_vrt

  integer, intent(in) :: mshf_n_elm
  integer, dimension(20), intent(in), target :: mshf_elm_idx
  integer, dimension(mshf_n_elm), intent(in), target :: mshf_elm_ty
  integer, dimension(mshf_n_elm), intent(in), target :: mshf_elm_dim
  integer, dimension(mshf_n_elm+1), intent(in), target :: mshf_elm_fce_ptr
  integer, dimension(mshf_elm_fce_ptr(mshf_n_elm + 1) + 1), &
    intent(in), target :: mshf_fce_vrt_ptr
  integer, dimension(mshf_fce_vrt_ptr(mshf_elm_fce_ptr(mshf_n_elm + 1) + 1)), &
    intent(in), target :: mshf_fce_vrt_lst

  integer, intent(in) :: mesh_n_bnd
  integer, dimension(mesh_n_bnd), intent(in), target :: mesh_bnd_idx
  integer, dimension(mesh_n_bnd + 1), intent(in), target :: mesh_bnd_tag_ptr
  integer, dimension(mesh_bnd_tag_ptr(mesh_n_bnd + 1)), &
    intent(in), target :: mesh_bnd_tag_lst
  integer, dimension(mesh_n_bnd + 1), intent(in), target :: mesh_bnd_vrt_ptr
  integer, dimension(mesh_bnd_vrt_ptr(mesh_n_bnd + 1)), &
    intent(in), target :: mesh_bnd_vrt_lst

  integer, intent(in) :: mesh_n_gbnd
  integer, dimension(mesh_n_gbnd), intent(in), target :: mesh_gbnd_idx
  integer, dimension(mesh_n_gbnd + 1), intent(in), target :: mesh_gbnd_tag_ptr
  integer, dimension(mesh_gbnd_tag_ptr(mesh_n_gbnd + 1)), &
    intent(in), target :: mesh_gbnd_tag_lst
  integer, dimension(mesh_n_gbnd + 1), intent(in), target :: mesh_gbnd_vrt_ptr
  integer, dimension(mesh_gbnd_vrt_ptr(mesh_n_gbnd + 1)), &
    intent(in), target :: mesh_gbnd_vrt_lst

  integer, intent(in) :: mesh_n_dmn
  integer, dimension(mesh_n_dmn), intent(in), target :: mesh_dmn_idx
  integer, dimension(mesh_n_dmn + 1), intent(in), target :: mesh_dmn_tag_ptr
  integer, dimension(mesh_dmn_tag_ptr(mesh_n_dmn + 1)), &
    intent(in), target :: mesh_dmn_tag_lst
  integer, dimension(mesh_n_dmn + 1), intent(in), target :: mesh_dmn_vrt_ptr
  integer, dimension(mesh_dmn_vrt_ptr(mesh_n_dmn + 1)), &
    intent(in), target :: mesh_dmn_vrt_lst

  integer, intent(in) :: mesh_n_gdmn
  integer, dimension(mesh_n_gdmn), intent(in), target :: mesh_gdmn_idx
  integer, dimension(mesh_n_gdmn + 1), intent(in), target :: mesh_gdmn_tag_ptr
  integer, dimension(mesh_gdmn_tag_ptr(mesh_n_gdmn + 1)), &
    intent(in), target :: mesh_gdmn_tag_lst
  integer, dimension(mesh_n_gdmn + 1), intent(in), target :: mesh_gdmn_vrt_ptr
  integer, dimension(mesh_gdmn_vrt_ptr(mesh_n_gdmn + 1)), &
    intent(in), target :: mesh_gdmn_vrt_lst

  integer, intent(in) :: mesh_n_xbnd
  integer, dimension(mesh_n_xbnd + 1), intent(in), target :: mesh_xbnd_vrt_ptr
  integer, dimension(mesh_xbnd_vrt_ptr(mesh_n_xbnd + 1)), &
    intent(in), target :: mesh_xbnd_vrt_lst
  integer, dimension(mesh_n_xbnd + 1), intent(in), target :: mesh_xbnd_elm_ptr
  integer, dimension(mesh_xbnd_elm_ptr(mesh_n_xbnd + 1)), &
    intent(in), target :: mesh_xbnd_elm_lst

  integer, intent(in) :: mesh_n_xdmn
  integer, dimension(mesh_n_xdmn + 1), intent(in), target :: mesh_xdmn_vrt_ptr
  integer, dimension(mesh_xdmn_vrt_ptr(mesh_n_xdmn + 1)), &
    intent(in), target :: mesh_xdmn_vrt_lst
  integer, dimension(mesh_n_xdmn + 1), intent(in), target :: mesh_xdmn_elm_ptr
  integer, dimension(mesh_xdmn_elm_ptr(mesh_n_xdmn + 1)), &
    intent(in), target :: mesh_xdmn_elm_lst

  integer, dimension(mesh_n_bnd + 1), intent(in), target :: mesh_bnd_neig_ptr
  integer, dimension(mesh_bnd_neig_ptr(mesh_n_bnd + 1)), &
    intent(in), target :: mesh_bnd_neig_lst
  integer, dimension(mesh_n_dmn + 1), intent(in), target :: mesh_dmn_neig_ptr
  integer, dimension(mesh_dmn_neig_ptr(mesh_n_dmn + 1)), &
    intent(in), target :: mesh_dmn_neig_lst

  integer :: err = 0
  integer, dimension(:), allocatable :: mpi_rqst
  integer, dimension(:,:), allocatable :: mpi_stat
  integer :: blk, n_blk, oblk
  type (mshf_t) :: mshf
  type (mesh_t) :: bnd, dmn
  type (geom_t) :: xbnd, xdmn

!   integer :: i, ii, j, jj1, jjn, k, m, n
  integer, dimension(1) :: lst1
  integer, dimension(:), allocatable :: itmp
  real(dp) :: dnum, bl_c, bl_d, vel_urf, pp_urf
  real(rp) :: cc, cd, summ, res0(10), gl_res0, time(2), res_drop
  real(rp), dimension(:), allocatable :: stmp
  real(rp), dimension(:,:), allocatable :: vtmp
  type(pde_t) :: vel, pp, pres
  integer :: iter, n_iter, fp_exit, n_gg
  integer :: ii, m, n, nnz, ipp, n_pp, n_pp_mod
  integer :: i,j,k,kii,kij,kji,kjj,ic,ib,ph_ty
  integer :: j1, j11, j1l, jb, jb1, jbl, jn, jnl, jj1, jjn, jj
  integer :: bnd1,bndl,bndu,dmn1,dmnl,dmnu
  integer :: f_len
  integer :: write_freq, write_pres, write_vel, write_spd, write_wall_force
  character(sl) :: out_file
  type(sfld_t) :: area
  integer,dimension(:),allocatable::bvrt_idx

  character(len=256):: bnd_type
  integer :: bnd_id
  real(rp),dimension(4) :: bnd_value

  type(fson_value), pointer :: json

  json => fson_parse(trim(filename))
  call fson_get(json, "Material.Density", cc)
  call fson_get(json, "Material.Viscosity", cd)
  call fson_get(json, "Solver.Iteration-Limit", n_iter)
  call fson_get(json, "Solver.Residual-Drop", res_drop)
  call fson_get(json, "Solver.Green-Gauss-Iterations", n_gg)
  call fson_get(json, "Solver.Pressure-Correction-Iterations", n_pp)
  call fson_get(json, "Solver.Pressure-Velocity-Scheme", pv_scheme)
  call fson_get(json, "Equation.Velocity.Relaxation", vel%urf)
  call fson_get(json, "Equation.Velocity.Iteration-Limit", vel%n_iter)
  call fson_get(json, "Equation.Velocity.Residual-Drop", vel%reduc)
  call fson_get(json, "Equation.Pressure.Relaxation", pp%urf)
  call fson_get(json, "Equation.Pressure.Iteration-Limit", pp%n_iter)
  call fson_get(json, "Equation.Pressure.Residual-Drop", pp%reduc)
  call fson_get(json, "Model.Navier-Stokes.Convection-Blending", vel%bl_c)
  call fson_get(json, "Model.Navier-Stokes.Diffusion-Blending", vel%bl_d)
  call fson_get(json, "Boundary.Type", bnd_type)
  call fson_get(json, "Boundary.Index", bnd_id)
  call fson_get(json, "Boundary.Velocity.X-Direction", bnd_value(1))
  call fson_get(json, "Boundary.Velocity.Y-Direction", bnd_value(2))
  call fson_get(json, "Boundary.Velocity.Z-Direction", bnd_value(3))
  call fson_get(json, "Boundary.Velocity.Magnitude", bnd_value(4))
  call fson_get(json, "Output.Frequency", write_freq)
  call fson_get(json, "Output.Pressure", write_pres)
  call fson_get(json, "Output.Velocity", write_vel)
  call fson_get(json, "Output.Speed", write_spd)
  call fson_get(json, "Output.Wall-Force", write_wall_force)
  call fson_destroy(json)

  !! initialize communications
  blk = 0; n_blk = 1;
#ifdef USE_MPI
  call MPI_Comm_size (mpi_comm_world, n_blk, err)
  call MPI_Comm_rank (mpi_comm_world, blk, err)
#endif

  n_proc = n_blk; proc = blk
#ifdef USE_MPI
  allocate (mpi_rqst(n_blk)); mpi_rqst = 0
  allocate (mpi_stat(mpi_status_size, n_blk)); mpi_stat = 0
#endif

  !! mesh format
  mshf%n_elm = mshf_n_elm
  mshf%elm_idx => mshf_elm_idx
  mshf%elm_ty => mshf_elm_ty
  mshf%elm_dim => mshf_elm_dim
  mshf%fce_vrt%blk => mshf_elm_fce_ptr
  mshf%fce_vrt%row => mshf_fce_vrt_ptr
  mshf%fce_vrt%col => mshf_fce_vrt_lst

  !! adjust mesh format to Fortran indexing
  mshf%elm_idx = eoshift (mshf%elm_idx, shift = 1)
  mshf%elm_idx = mshf%elm_idx + 1
  mshf%fce_vrt%blk = mshf%fce_vrt%blk + 1
  mshf%fce_vrt%row = mshf%fce_vrt%row + 1
  mshf%fce_vrt%col = mshf%fce_vrt%col + 1

  !! basic geometry
  dmn%n_dim = mesh_n_dim
  dmn%n_vrt = mesh_n_vrt
  dmn%vrt_idx => mesh_vrt_idx

!   xdmn%x_vrt => geom_x_vrt
  allocate (xdmn%x_vrt(dmn%n_dim, dmn%n_vrt))
  xdmn%x_vrt = geom_x_vrt(:dmn%n_dim,:)

  !! mesh boundary
  bnd%blk = mesh_blk
  bnd%n_blk = mesh_n_blk
  bnd%elm_blk => mesh_bnd_blk
  bnd%n_elm = mesh_n_bnd
  bnd%elm_idx => mesh_bnd_idx
  bnd%elm_tag%row => mesh_bnd_tag_ptr
  bnd%elm_tag%col => mesh_bnd_tag_lst
  bnd%elm_vrt%row => mesh_bnd_vrt_ptr
  bnd%elm_vrt%col => mesh_bnd_vrt_lst
  bnd%elm_neig%row => mesh_bnd_neig_ptr
  bnd%elm_neig%col => mesh_bnd_neig_lst

  bnd%n_ghst = mesh_n_gbnd
  bnd%ghst_idx => mesh_gbnd_idx
  bnd%ghst_tag%row => mesh_gbnd_tag_ptr
  bnd%ghst_tag%col => mesh_gbnd_tag_lst
  bnd%ghst_vrt%row => mesh_gbnd_vrt_ptr
  bnd%ghst_vrt%col => mesh_gbnd_vrt_lst

  bnd%n_fce = mesh_n_xbnd
  bnd%fce_elm%row => mesh_xbnd_elm_ptr
  bnd%fce_elm%col => mesh_xbnd_elm_lst
  bnd%fce_vrt%row => mesh_xbnd_vrt_ptr
  bnd%fce_vrt%col => mesh_xbnd_vrt_lst

  !! adjust mesh boundary to Fortran indexing
  bnd%blk = bnd%blk + 1
  bnd%elm_blk = bnd%elm_blk + 1

  bnd%elm_tag%row = bnd%elm_tag%row + 1
  bnd%elm_vrt%row = bnd%elm_vrt%row + 1
  bnd%elm_neig%row = bnd%elm_neig%row + 1
  bnd%elm_neig%col = bnd%elm_neig%col + 1

  bnd%ghst_tag%row = bnd%ghst_tag%row + 1
  bnd%ghst_vrt%row = bnd%ghst_vrt%row + 1

  bnd%fce_elm%row = bnd%fce_elm%row + 1
  bnd%fce_vrt%row = bnd%fce_vrt%row + 1

  !! mesh domain
  dmn%blk = mesh_blk
  dmn%n_blk = mesh_n_blk
  dmn%elm_blk => mesh_dmn_blk
  dmn%n_elm = mesh_n_dmn
  dmn%elm_idx => mesh_dmn_idx
  dmn%elm_tag%row => mesh_dmn_tag_ptr
  dmn%elm_tag%col => mesh_dmn_tag_lst
  dmn%elm_vrt%row => mesh_dmn_vrt_ptr
  dmn%elm_vrt%col => mesh_dmn_vrt_lst
  dmn%elm_neig%row => mesh_dmn_neig_ptr
  dmn%elm_neig%col => mesh_dmn_neig_lst

  dmn%n_ghst = mesh_n_gdmn
  dmn%ghst_idx => mesh_gdmn_idx
  dmn%ghst_tag%row => mesh_gdmn_tag_ptr
  dmn%ghst_tag%col => mesh_gdmn_tag_lst
  dmn%ghst_vrt%row => mesh_gdmn_vrt_ptr
  dmn%ghst_vrt%col => mesh_gdmn_vrt_lst

  dmn%n_fce = mesh_n_xdmn
  dmn%fce_elm%row => mesh_xdmn_elm_ptr
  dmn%fce_elm%col => mesh_xdmn_elm_lst
  dmn%fce_vrt%row => mesh_xdmn_vrt_ptr
  dmn%fce_vrt%col => mesh_xdmn_vrt_lst

  !! adjust mesh domain to Fortran indexing
  dmn%blk = dmn%blk + 1
  dmn%elm_blk = dmn%elm_blk + bnd%elm_blk(bnd%n_blk + 1)

  dmn%elm_tag%row = dmn%elm_tag%row + 1
  dmn%elm_vrt%row = dmn%elm_vrt%row + 1
  dmn%elm_neig%row = dmn%elm_neig%row + 1
  dmn%elm_neig%col = dmn%elm_neig%col + 1

  dmn%ghst_tag%row = dmn%ghst_tag%row + 1
  dmn%ghst_vrt%row = dmn%ghst_vrt%row + 1

  dmn%fce_elm%row = dmn%fce_elm%row + 1
  dmn%fce_vrt%row = dmn%fce_vrt%row + 1

  !! setup vertex index mapping
  call hash_initialise (dmn%vi2i, dmn%n_vrt)
  do i = 1, dmn%n_vrt
    call hash_insert (dmn%vi2i, dmn%vrt_idx(i), i)
  end do

  !! setup ghost index and block mappings
  call hash_initialise (bnd%gi2i, bnd%n_ghst)
  call hash_initialise (bnd%gi2b, bnd%n_ghst)
  do i=1,bnd%n_ghst
    jjn = bnd%ghst_tag%row(i + 1) - 1
    oblk = (-1) * bnd%ghst_tag%col(jjn)
    call hash_insert (bnd%gi2b, bnd%ghst_idx(i), oblk)
    call hash_insert (bnd%gi2i, bnd%ghst_idx(i), i)
  end do

  call hash_initialise (dmn%gi2i, dmn%n_ghst)
  call hash_initialise (dmn%gi2b, dmn%n_ghst)
  do i=1,dmn%n_ghst
    jjn = dmn%ghst_tag%row(i + 1) - 1
    oblk = (-1) * dmn%ghst_tag%col(jjn)
    call hash_insert (dmn%gi2b, dmn%ghst_idx(i), oblk)
    call hash_insert (dmn%gi2i, dmn%ghst_idx(i), i)
  end do

#ifdef USE_MPI
  call init_ghost_update (bnd, dmn)
#endif

  !! initialize domain geometry
  allocate (xdmn%x_vc%dmn(dmn%n_dim, dmn%n_elm))
  allocate (xdmn%x_vc%ldmn(dmn%n_dim, dmn%n_locl))
  allocate (xdmn%x_vc%gdmn(dmn%n_dim, dmn%n_ghst))

  allocate (xdmn%vol%dmn(dmn%n_elm))
  allocate (xdmn%vol%ldmn(dmn%n_locl))
  allocate (xdmn%vol%gdmn(dmn%n_ghst))

  allocate (xdmn%gl_min%dmn(dmn%n_dim, dmn%n_elm))
  allocate (xdmn%gl_min%ldmn(dmn%n_dim, dmn%n_locl))
  allocate (xdmn%gl_min%gdmn(dmn%n_dim, dmn%n_ghst))
  allocate (xdmn%gl_max%dmn(dmn%n_dim, dmn%n_elm))
  allocate (xdmn%gl_max%ldmn(dmn%n_dim, dmn%n_locl))
  allocate (xdmn%gl_max%gdmn(dmn%n_dim, dmn%n_ghst))
  allocate (xdmn%gl%dmn(dmn%n_dim, dmn%n_elm))
  allocate (xdmn%gl%ldmn(dmn%n_dim, dmn%n_locl))
  allocate (xdmn%gl%gdmn(dmn%n_dim, dmn%n_ghst))

  allocate (xdmn%ls%dmn(dmn%n_dim, dmn%n_dim, dmn%n_elm))
  allocate (xdmn%ls%ldmn(dmn%n_dim, dmn%n_dim, dmn%n_locl))
  allocate (xdmn%ls%gdmn(dmn%n_dim, dmn%n_dim, dmn%n_ghst))
  allocate (xdmn%gg%dmn(dmn%n_dim, dmn%n_dim, dmn%n_elm))
  allocate (xdmn%gg%ldmn(dmn%n_dim, dmn%n_dim, dmn%n_locl))
  allocate (xdmn%gg%gdmn(dmn%n_dim, dmn%n_dim, dmn%n_ghst))

  allocate (xdmn%x_fc(dmn%n_dim, dmn%n_fce))
  allocate (xdmn%x_pc(dmn%n_dim, dmn%n_fce))
  allocate (xdmn%x_nc(dmn%n_dim, dmn%n_fce))
  allocate (xdmn%norm(dmn%n_dim, dmn%n_fce))
  allocate (xdmn%volr(dmn%n_elm))
  allocate (xdmn%area(dmn%n_fce))
  allocate (xdmn%l_pn(dmn%n_fce))
  allocate (xdmn%w_fp(dmn%n_fce))

  ! call print_range("bnd%elm_blk [f90]", bnd%n_blk, bnd%elm_blk);
  ! call print_range("dmn%elm_blk [f90]", dmn%n_blk, dmn%elm_blk);
  ! call print_graph("dmn%fce_elm [f90]", dmn%n_fce, dmn%fce_elm%row, dmn%fce_elm%col);

  !! calculate domain geometry
  call geometry (bnd, dmn, xbnd, xdmn)

  !! check total volume
#ifdef USE_MPI
  if (n_blk > 1) then
    summ = xdmn%vol_sum
    call MPI_Reduce (summ, xdmn%vol_sum, 1, mpi_rp, &
      mpi_sum, 0, mpi_comm_world, err)
  end if
#endif

  if (blk == 0) print "(a,1x,g12.4)","Volume", xdmn%vol_sum

  !! compute LS gradient matrix
  call ls_grad_matrix (bnd, dmn, xbnd, xdmn)

  !! initialize linear solver
  call parms_init (dmn%n_elm)

  !! initialize partial differential equations
  vel%n_cmp = dmn%n_dim

  allocate (vel%bnd_c(vel%n_cmp, 49))
  allocate (vel%bnd_d(vel%n_cmp, 49))

  allocate (vel%grad%dmn(dmn%n_dim, vel%n_cmp, dmn%n_elm))
  allocate (vel%grad%ldmn(dmn%n_dim, vel%n_cmp, dmn%n_locl))
  allocate (vel%grad%gdmn(dmn%n_dim, vel%n_cmp, dmn%n_ghst))

  allocate (vel%phi%dmn(vel%n_cmp, dmn%n_elm))
  allocate (vel%phi%ldmn(vel%n_cmp, dmn%n_locl))
  allocate (vel%phi%gdmn(vel%n_cmp, dmn%n_ghst))
  allocate (vel%phi%bnd(vel%n_cmp, bnd%n_elm))

  allocate (vel%rhs(vel%n_cmp, dmn%n_elm))
  allocate (vel%spd(dmn%n_fce))

  allocate (vel%cc%dmn(dmn%n_elm))
  allocate (vel%cc%ldmn(dmn%n_locl))
  allocate (vel%cc%gdmn(dmn%n_ghst))
  allocate (vel%cc%bnd(bnd%n_elm))

  allocate (vel%cd%dmn(dmn%n_elm))
  allocate (vel%cd%ldmn(dmn%n_locl))
  allocate (vel%cd%gdmn(dmn%n_ghst))
  allocate (vel%cd%bnd(bnd%n_elm))

  allocate (vel%diag%dmn(dmn%n_elm))
  allocate (vel%diag%ldmn(dmn%n_locl))
  allocate (vel%diag%gdmn(dmn%n_ghst))

  m = dmn%n_elm
  nnz = dmn%elm_neig%row(m + 1) - 1
  allocate (vel%mat%idx(m))
  allocate (vel%mat%val(nnz))
  vel%mat%idx = dmn%elm_idx - dmn%elm_blk(1) + 1
  vel%mat%row => dmn%elm_neig%row
  vel%mat%col => dmn%elm_neig%col

  vel%calc_spd = 1
  vel%spd = zero
  vel%cc%dmn = cc; vel%cc%gdmn = cc; vel%cc%bnd = cc
  vel%cd%dmn = cd; vel%cd%gdmn = cd; vel%cd%bnd = cd
  vel%phi%dmn = zero; vel%phi%gdmn = zero; vel%phi%bnd = zero

!   call random_number (vel%phi%dmn)
#ifdef USE_MPI
  if (n_blk > 1) call ghost_update_v (bnd, dmn, vel%phi, 0, 1)
#endif

  vel%rhs = zero
  vel%mat%val = zero

  !! initialize pressure correction equation
  pp%n_cmp = 1

  allocate (pp%grad%dmn(dmn%n_dim, pp%n_cmp, dmn%n_elm))
  allocate (pp%grad%ldmn(dmn%n_dim, pp%n_cmp, dmn%n_locl))
  allocate (pp%grad%gdmn(dmn%n_dim, pp%n_cmp, dmn%n_ghst))

  allocate (pp%phi%dmn(pp%n_cmp, dmn%n_elm))
  allocate (pp%phi%ldmn(pp%n_cmp, dmn%n_locl))
  allocate (pp%phi%gdmn(pp%n_cmp, dmn%n_ghst))
  allocate (pp%phi%bnd(pp%n_cmp, bnd%n_elm))

  allocate (pp%ref(pp%n_cmp))

  pp%spd => vel%spd
  pp%cc%dmn => vel%cc%dmn
  pp%cc%ldmn => vel%cc%ldmn
  pp%cc%gdmn => vel%cc%gdmn
  pp%cc%bnd => vel%cc%bnd
  pp%rhs => vel%rhs
  pp%diag%dmn => vel%diag%dmn
  pp%diag%ldmn => vel%diag%ldmn
  pp%diag%gdmn => vel%diag%gdmn

  pp%mat%idx => vel%mat%idx
  pp%mat%val => vel%mat%val
  pp%mat%row => vel%mat%row
  pp%mat%col => vel%mat%col

  pp%ref = zero
  pp%phi%dmn = zero; pp%phi%gdmn = zero; pp%phi%bnd = zero

  !! initialize pressure field
  pres%n_cmp = 1

  allocate (pres%grad%dmn(dmn%n_dim, pres%n_cmp, dmn%n_elm))
  allocate (pres%grad%ldmn(dmn%n_dim, pres%n_cmp, dmn%n_locl))
  allocate (pres%grad%gdmn(dmn%n_dim, pres%n_cmp, dmn%n_ghst))

  allocate (pres%phi%dmn(pres%n_cmp, dmn%n_elm))
  allocate (pres%phi%ldmn(pres%n_cmp, dmn%n_locl))
  allocate (pres%phi%gdmn(pres%n_cmp, dmn%n_ghst))
  allocate (pres%phi%bnd(pres%n_cmp, bnd%n_elm))

  pres%phi%dmn = zero; pres%phi%gdmn = zero; pres%phi%bnd = zero

  call set_boundary_conditions()

  ! workspace
  allocate (wrk(3))
  allocate (wrk(1)%d(size(vel%mat%val)))
  allocate (wrk(2)%d(dmn%n_elm))
  allocate (wrk(3)%d(dmn%n_elm))

  ! iteratively solve PDEs
  do iter = 1, n_iter

     ! initialize solution
    if (iter > write_freq) then
       vel%calc_d = 1
       vel%bl_c = bl_c;
       vel%bl_d = bl_d;
       vel%urf = vel_urf;
       n_pp_mod = n_pp
    else
       bl_c = vel%bl_c; vel%bl_c = 0;
       bl_d = vel%bl_d; vel%bl_d = 0;
       vel_urf = vel%urf; vel%urf = 0.9;
       vel%calc_d = 0
       n_pp_mod = 0
    end if

    ! initialize workspace
    vel%rhs = zero
    vel%mat%val = zero

    ! pressure gradient
    if (n_gg == 0) call ls_gradient (bnd, dmn, xbnd, xdmn, pres%phi, pres%grad, 1)
    if (n_gg > 0) call gg_gradient (bnd, dmn, xbnd, xdmn, pres%phi, pres%grad, 1, n_gg)
#ifdef USE_MPI
    if (n_blk > 1) call ghost_update_t (bnd, dmn, pres%grad, 0, 1)
#endif
    call gl_local_extrema (bnd, dmn, xbnd, xdmn, pres%phi)
    call grad_limiter (bnd, dmn, xbnd, xdmn, pres%phi, pres%grad, 2)
#ifdef USE_MPI
    if (n_blk > 1) call ghost_update_t (bnd, dmn, pres%grad, 0, 1)
#endif

    ! velocity gradient
    if (n_gg == 0) call ls_gradient (bnd, dmn, xbnd, xdmn, vel%phi, vel%grad, 0)
    if (n_gg > 0) call gg_gradient (bnd, dmn, xbnd, xdmn, vel%phi, vel%grad, 0, n_gg)
#ifdef USE_MPI
    if (n_blk > 1) call ghost_update_t (bnd, dmn, vel%grad, 0, 1)
#endif
    call gl_local_extrema (bnd, dmn, xbnd, xdmn, vel%phi)
    call grad_limiter (bnd, dmn, xbnd, xdmn, vel%phi, vel%grad, 2)
#ifdef USE_MPI
    if (n_blk > 1) call ghost_update_t (bnd, dmn, vel%grad, 0, 1)
#endif

    ! momentum equation
    call ns_source (bnd, dmn, xbnd, xdmn, pres, vel)
    call construct_pde (bnd, dmn, xbnd, xdmn, vel)

#ifdef USE_MPI
    if (n_blk > 1) call ghost_update_s (bnd, dmn, vel%diag, 0, 1)
#endif
    do i = 1, vel%n_cmp
      if (i == 1) wrk(1)%d = real (vel%mat%val, dp)
      wrk(2)%d = real (vel%rhs(i,:), dp)
      wrk(3)%d = real (vel%phi%dmn(i,:), dp)

      call parms_solv ( &
        dmn%n_elm, vel%mat%idx, vel%mat%row, vel%mat%col, &
        wrk(1)%d, wrk(2)%d, wrk(3)%d, dnum)

      vel%phi%dmn(i,:) = real (wrk(3)%d, rp)
      res0(i) = dnum

#ifndef USE_MPI
      call lin_eqns_solver ( &
        vel%mat%val,vel%phi%dmn(i,:),vel%rhs(i,:), &
        res0(i),vel%mat%col,vel%mat%row,vel%n_iter,vel%reduc,2,i)
#endif
    end do

#ifdef USE_MPI
    if (n_blk > 1) call ghost_update_v (bnd, dmn, vel%phi, 0, 1)
#endif

    ! initialize workspace
    pp%rhs = zero
    pp%mat%val = zero

    ! pressure correction equation
    call construct_pp (bnd, dmn, xbnd, xdmn, vel, pres, pp)

    do ipp = 1, n_pp_mod

      i = 4
      pp%phi%dmn(1,:) = zero

      wrk(1)%d = real (pp%mat%val, dp)
      wrk(2)%d = real (pp%rhs(1,:), dp)
      wrk(3)%d = real (pp%phi%dmn(1,:), dp)

      call parms_solv ( &
        dmn%n_elm, pp%mat%idx, pp%mat%row, pp%mat%col, &
        wrk(1)%d, wrk(2)%d, wrk(3)%d, dnum)

      pp%phi%dmn(1,:) = real (wrk(3)%d, rp)
      res0(i) = dnum

#ifndef USE_MPI
      call lin_eqns_solver ( &
        pp%mat%val,pp%phi%dmn(1,:),pp%rhs(1,:), &
        res0(i),vel%mat%col,vel%mat%row,pp%n_iter,pp%reduc,1,i)
#endif

#ifdef USE_MPI
      if (n_blk > 1) call ghost_update_v (bnd, dmn, pp%phi, 0, 1)
#endif

      ! pressure correction gradient
      if (n_gg == 0) call ls_gradient (bnd, dmn, xbnd, xdmn, pp%phi, pp%grad, 1)
      if (n_gg > 0) call gg_gradient (bnd, dmn, xbnd, xdmn, pp%phi, pp%grad, 1, n_gg)
#ifdef USE_MPI
      if (n_blk > 1) call ghost_update_t (bnd, dmn, pp%grad, 0, 1)
#endif
      call gl_local_extrema (bnd, dmn, xbnd, xdmn, pp%phi)
      call grad_limiter (bnd, dmn, xbnd, xdmn, pp%phi, pp%grad, 2)
#ifdef USE_MPI
      if (n_blk > 1) call ghost_update_t (bnd, dmn, pp%grad, 0, 1)
#endif

      ! correct velocity and pressure
      call ns_correction (bnd, dmn, xbnd, xdmn, pp, vel, pres)
#ifdef USE_MPI
      if (n_blk > 1) call ghost_update_v (bnd, dmn, vel%phi, 0, 1)
      if (n_blk > 1) call ghost_update_v (bnd, dmn, pres%phi, 0, 1)
#endif

      pp%rhs = zero
      call correct_skew_pp (bnd, dmn, xbnd, xdmn, vel, pres, pp)
    end do

    if (blk == 0) then
      call fp_iter_ctrl(iter, n_iter, fp_exit, res_drop, res0, gl_res0)
    end if

    if (((write_freq > 0) .and. mod(iter, write_freq) == 0)) then
       call write_fields ()
    end if

#ifdef USE_MPI
    call MPI_Bcast (fp_exit, 1, mpi_integer, 0, mpi_comm_world, err)
#endif

    if ((iter > write_freq) .and. (fp_exit > 0)) exit
  end do

contains
subroutine write_fields ()
  if (write_vel == 1) then
     out_file = dir (filename)
     out_file = trim(out_file)//"vel.msh"
     if (n_blk > 1) then
        out_file = trim(out_file)//"_p"//trim(itoa(blk))
     end if
     f_len = len_trim(out_file)

     call write_gmsh_vector_field (f_len, out_file, dmn%n_elm, vel%n_cmp, dmn%elm_idx, vel%phi%dmn)
  end if

  if (write_spd == 1) then
     out_file = dir (filename)
     out_file = trim(out_file)//"spd.msh"
     if (n_blk > 1) then
        out_file = trim(out_file)//"_p"//trim(itoa(blk))
     end if
     f_len = len_trim(out_file)

     wrk(3)%d = real (vvec_mag(vel%phi%dmn), dp)
     call write_gmsh_field (f_len, out_file, dmn%n_elm, dmn%elm_idx, wrk(3)%d)
  end if

  if (write_pres == 1) then
     out_file = dir (filename)
     out_file = trim(out_file)//"pres.msh"
     if (n_blk > 1) then
        out_file = trim(out_file)//"_p"//trim(itoa(blk))
     end if
     f_len = len_trim(out_file)

     wrk(3)%d = real (pres%phi%dmn(1,:), dp)
     call write_gmsh_field (f_len, out_file, dmn%n_elm, dmn%elm_idx, wrk(3)%d)
  end if

  if (write_wall_force > 0) then
     call wall_force(write_wall_force);
  end if
end subroutine

subroutine set_boundary_conditions()
  integer :: bnd1,bndl,bndu,dmn1,dmnl,dmnu
  integer :: i,j,k,kii,kij,kji,kjj,ic,ib,ph_ty
  integer :: j1, j11, j1l, jb, jb1, jbl, jn, jnl, jj1, jjn, jj
  integer :: bnd_os

  print "(a,i2,4g12.4)", trim(bnd_type), bnd_id, bnd_value

  call boundary_offset(bnd_type, bnd_os)
  bnd_id = bnd_id + bnd_os

  bnd1 = bnd%elm_blk(1)
  bndl = bnd%elm_blk(bnd%blk)
  bndu = bnd%elm_blk(bnd%blk + 1)

  dmn1 = dmn%elm_blk(1)
  dmnl = dmn%elm_blk(dmn%blk)
  dmnu = dmn%elm_blk(dmn%blk + 1)

  k = 0
  do i = 1, dmn%n_fce
    jj1 = 2 * (i - 1) + 1
    jjn = 2 * (i)

    j1 = dmn%fce_elm%col(jj1)
    j11 = (j1 - dmn1) + 1
    j1l = (j1 - dmnl) + 1
    jn = dmn%fce_elm%col(jjn)

    if (jn < 0) then
      jb1 = ((-1)*jn - bnd1) + 1
      jbl = ((-1)*jn - bndl) + 1
      ph_ty = bnd%elm_tag%row(jbl)
      ph_ty = bnd%elm_tag%col(ph_ty+1)

      if (ph_ty == bnd_id) then
         if (vec_mag(bnd_value(1:3)) < 1.0e-20) then
            vel%phi%bnd(:,jbl) = -xdmn%norm(:,i) * bnd_value(4)
         else
            vel%phi%bnd(1,jbl) = bnd_value(1) * bnd_value(4)
            vel%phi%bnd(2,jbl) = bnd_value(2) * bnd_value(4)
            if (bnd%n_dim == 3) vel%phi%bnd(3,jbl) = bnd_value(3) * bnd_value(4)
         end if
      end if
    end if
  end do
end subroutine

subroutine wall_force(wall_id)
  integer,intent(in)::wall_id

  integer :: bnd1,bndl,bndu,dmn1,dmnl,dmnu
  integer :: i,j,k,kii,kij,kji,kjj,ic,ib,ph_ty
  integer :: j1, j11, j1l, jb, jb1, jbl, jn, jnl, jj1, jjn, jj
  integer :: bnd_os
  character(len=256):: wall_type

  real(rp)::d,fpres,area,l_pn,cd
  real(rp),dimension(mesh_n_dim)::n,t,f_shear,f_pres,vvel,fvel,f_wall

  f_wall=0; f_pres=0; f_shear=0;

  wall_type = "Wall";
  call boundary_offset(wall_type, bnd_os)

  bnd1 = bnd%elm_blk(1)
  bndl = bnd%elm_blk(bnd%blk)
  bndu = bnd%elm_blk(bnd%blk + 1)

  dmn1 = dmn%elm_blk(1)
  dmnl = dmn%elm_blk(dmn%blk)
  dmnu = dmn%elm_blk(dmn%blk + 1)

  k = 0
  do i = 1, dmn%n_fce
    jj1 = 2 * (i - 1) + 1
    jjn = 2 * (i)

    j1 = dmn%fce_elm%col(jj1)
    j11 = (j1 - dmn1) + 1
    j1l = (j1 - dmnl) + 1
    jn = dmn%fce_elm%col(jjn)

    if (jn < 0) then
      jb1 = ((-1)*jn - bnd1) + 1
      jbl = ((-1)*jn - bndl) + 1
      ph_ty = bnd%elm_tag%row(jbl)
      ph_ty = bnd%elm_tag%col(ph_ty+1)

      if (ph_ty == (wall_id + bnd_os)) then
         area = xdmn%area(i)
         n = xdmn%norm(:,i)
         l_pn = xdmn%l_pn(i)
         cd = vel%cd%bnd(jbl)

         d = cd * area/l_pn
         vvel = vel%phi%dmn(:,j1l)
         fvel = vel%phi%bnd(:,jbl)
         fpres = pres%phi%bnd(1,jbl)

         t = fvel+vvel
         t = t - n*dot_product(t,n)
         t = t / (vec_mag(t)+small)

         f_shear = -d*(fvel-t*dot_product(vvel,t))
         f_pres = fpres*n * area
         f_wall = f_wall+(f_shear+f_pres)
      end if
    end if
  end do

  if (mesh_n_dim == 2) then
     print "(a,i1,a,2es11.3,a)", "Wall (",wall_id,") force = [", f_wall,"]"
  else
     print "(a,i1,a,3es11.3,a)", "Wall (",wall_id,") force = [", f_wall,"]"
  end if
end subroutine

end subroutine

max_stars_repo_name: DominicJones/pcns

max_issues_repo_name: DominicJones/pcns

max_forks_repo_name: DominicJones/pcns


text:
**************************************************************************
*                                                                        *
*                  NASA/Goddard Space Flight Center                      *
*                     Laboratory for Atmospheres                         *
*                      Data Assimilation Office                          *
*                            Code 910.3                                  *
*                                                                        *
*            PHYSICAL-SPACE STATISTICAL ANALYSIS SYSTEM (PSAS)           *
*                                                                        *
**************************************************************************

*........................... ROUTINE PROLOGUE .............................
* !BOP
*
* !FILE: gridxx.f
*
* !ROUTINE: gridxx0
*
* !DESCRIPTION: Initialize analysis grid parameters.
*
* !CALLING SEQUENCE:
*
*         call GRIXX0
*
* !INPUT PARAMETERS: none
*
* !OUTPUT PARAMETERS: none
*
* !SYSTEM ROUTINES: none.
*
* !FILES USED: none.
*
* !WRITTEN BY: Jim Pf, 02jun93
*
* !REVISION HISTORY:
!
!	21Oct99	- Jing Guo
!		. Replaced qea.h dependency with m_qea(.F90).
!		. Removed stdio.h dependency by removing write()
!		  statements.
!	01Oct99	- Jing Guo
!		. Removed gridxx.h, since it is not used in here.
!
*  18Jul95 - Jing G.	- Splited gridxx.h and qea.h for variable grid
*			  settings
*  16nov94   A. da Silva  Implemented quasi-equal area grid option.
*                         Initialization of mandatory levels moved
*                         to gridxx.h for portability.
*
* !EOP
*.........................................................................

      subroutine GRIDXX0

	use m_inpak90, only : lablin,fltget
	use m_qea,     only : eaytresh
	implicit none

*     Parameters and table for analysis grid.
*     --------------------------------------
      call LABLIN ( 'latitude_treshold_for_equal_area_grid:' )
      eaytresh = abs(FLTGET ( 90. ))

      return
      end

*........................... ROUTINE PROLOGUE .............................
* !BOP
*
* !FILE: gridxx.f
*
* !ROUTINE: eagrid_set
*
* !DESCRIPTION:
*
*     Returns quasi-equal area grid points on the sphere poleward
*  of ytresh, regular lat/lon grid elsewhere. If ytresh = 90, then
*  a regular lat/lon grid is returned.
*
* !CALLING SEQUENCE:
*
*      call EAGRID_set ( ier, glon, glat, ndecl, ngrid,
*    .               idim, jdim )
*
* !INPUT PARAMETERS:
*
*      integer  ndecl     ! declared dimension in calling program
*      integer  idim      ! no. of longitudes at equator
*      integer  jdim      ! no. of equally space latitudes
*      integer  npole     ! no. zonal grid points at pole
*      real     ytresh    ! only change grid poleward of ytresh
*
* !OUTPUT PARAMETERS:
*
*      integer  ier           ! 0 if OK, .ne. 0 otherwise
*      integer  ngrid         ! total number of grid points in
*                             ! quasi equal area grid.
*      real     glon(ndecl)   ! longitude in degree
*      real     glat(ndecl)   ! latitude in degree
*
* !SEE ALSO: n/a
*
* !SYSTEM ROUTINES: none.
*
* !FILES USED: none.
*
* !WRITTEN BY: A. da Silva, 16nov94.
*
* !REVISION HISTORY:
!
!	08Nov99	- Jing Guo
!		. renamed EAGRID() to EAGRID_set, since the original
!		  procedure name conflicts with the module datatype
!		  name of m_EAGrid(.F90).
!	01Oct99	- Jing Guo
!		. Removed gridxx.h, since it is not used here.
*
* !EOP
*.........................................................................

      subroutine EAGRID_set ( ier, glon, glat, ndecl, ngrid,
     .                    idim, jdim )

	use m_qea,only : qea_init
	use m_qea,only : lea_beg,lea_len,ea_lon
	use m_qea,only : npole,eaytresh
	use m_qea,only : j_south,j_north

	implicit none

*     On output
*     ---------
      integer  ier           ! 0 if OK, .ne. 0 otherwise
      integer  ngrid         ! total number of grid points in
                             ! quasi equal area grid.

      integer  ndecl     ! declared dimension in calling program

      real     glon(ndecl)   ! longitude in degree
      real     glat(ndecl)   ! latitude in degree

*     On input:
*     --------
      integer  idim      ! no. of longitudes at equator
      integer  jdim      ! no. of equally space latitudes

! Local variables

	integer :: i,j,l
	integer :: nlon,neq
	real :: dlon,dlat,ylat,alat
	real :: pi,ytresh,d2r
*
*     Returns quasi-equal area grid points on the sphere.
*     This version defines a peculiar quasi-equal area grid at
*     the polar caps (northward of ytresh below). To simplify the
*     interpolation, the poles have a full set of longitudes
*     (usually the longitudinal mesh size decreases as the poles
*     are approached.)

	call qea_init(idim,jdim)

      ytresh = eaytresh
      pi = 4. * atan ( 1. )
      d2r = pi / 180.

*     Generate quasi-equal area grid
*     ------------------------------
      neq  = idim - npole
      dlat = 180. / ( jdim - 1 )
      l = 0
      do 10 j = 1, jdim

         ylat = -90. + (j-1) * dlat
         alat = d2r * ylat

*        Decide how many long. gridpoints for this latitude
*        --------------------------------------------------
         if ( abs(ylat) .gt. ytresh ) then
              nlon = nint ( npole + neq * cos(alat) )
         else
              nlon = idim
         end if
         if (j.eq.1 .or. j.eq.jdim) nlon = idim ! no interp at poles
         dlon = 360. / nlon

         lea_beg(j) = l + 1    ! pointer for interpolation
         lea_len(j) = nlon

*        Create gridpoints for this latitude
*        -----------------------------------
         do 20 i = 1, nlon
            l = l + 1
            if ( l .gt. ndecl ) then
               ier = 1
               return
            end if
            glon(l) = ( i - 1 ) * dlon - 180.  ! this will be sorted
            glat(l) = ylat                     ! this will be sorted
            ea_lon(l) = glon(l)                ! this will be kept as is

 20      continue

 10   continue

*     No. of q.e.a. grid points
*     -------------------------
      ngrid = l

*     Compute beginning and end latitudinal index of
*      regular lat/lon grid (based on ytresh) needed
*      for interpolation
*     ----------------------------------------------
      do 50 j = 1, jdim
         ylat = -90. + (j-1) * dlat
         j_south = j
         if ( abs(ylat) .le. ytresh ) go to 51
 50   continue
 51   continue
      do 60 j = jdim, 1, -1
         ylat = -90. + (j-1) * dlat
         j_north = j
         if ( abs(ylat) .le. ytresh ) go to 61
 60   continue
 61   continue

      ier = 0

*     All done.
*     ---------
      ier = 0
      return
      end

*........................... ROUTINE PROLOGUE .............................
* !BOP
*
* !FILE: gridxx.f
*
* !ROUTINE: ea2ll
*
* !DESCRIPTION: Converts from quasi-equal area to a lat/lon grid.
*
* !CALLING SEQUENCE:
*
* !INPUT PARAMETERS:
*
* !OUTPUT PARAMETERS:
*
* !SEE ALSO:
*
* !SYSTEM ROUTINES:
*
* !FILES USED:
*
* !WRITTEN BY: A. da Silva, 16nov94.
*
* !REVISION HISTORY:
!
!	01Oct99	- Jing Guo
!		. Removed gridxx.h, since it is not used here.
*
* !EOP
*.........................................................................

      subroutine EA2LL ( All, idim, jdim, Aea, ngrid )
	use m_qea,only : lea_beg,lea_len,ea_lon
	use m_qea,only : eaytresh
	use m_qea,only : j_south,j_north
	use m_stdio,only : stderr
	use m_die,only : die
	implicit none

	integer,intent(in) :: idim
	integer,intent(in) :: jdim

*     Output:
*     ------
      real,intent(out) :: All(idim,jdim)      ! array on lat/lon grid

*     Input:
*     -----
      integer,intent(in) :: ngrid ! No. of quasi-equal area grid points
      real,   intent(in) :: Aea(ngrid) ! array on quasi-equal area grid

      character(len=*), parameter :: myname='ea2ll'

*     Local work space
*     ----------------
	integer nw,ier
	real, allocatable :: b(:),c(:),d(:),x(:),y(:)	! spline coef.
	real :: dlon,dlat,ytresh,u
	integer :: i,j,l,n

	real,external :: SEVAL

	nw=idim+1

      ytresh = eaytresh

*     For latitudes equatorward of ytresh no interpolation is needed
*     --------------------------------------------------------------
      l = lea_beg(j_south)
      do j = j_south, j_north
         do i = 1, idim
            All(i,j) = Aea(l)
            l = l + 1
         end do
      end do

*     No interpolation needed if entire grid is lat/lon
*     -------------------------------------------------

      if ( ytresh .ge. 90.-90./(jdim-1) ) return

	allocate( b(nw), c(nw), d(nw), x(nw), y(nw), stat=ier)
	if(ier.ne.0) then
	  write(stderr,'(2a,i5)') myname,
     &	    ': allocate() error, stat =',ier
	  call die(myname)
	endif

*     Southern polar cap
*     ------------------
      dlon = 360. / idim
      do j = 2, j_south - 1                 ! no south pole

c-       ylat = -90. + (j-1)*dlat

         l = lea_beg(j)
         n = lea_len(j)

*        Copy q.e.a. data to local vector repeating first
*         longitudinal point to ensure periodicity
*        ------------------------------------------------
         do i = 1, n
            x(i) = ea_lon(l+i-1)
            y(i) = Aea(l+i-1)
         end do
         n = n + 1
         x(n) = +180.      ! -180 = + 180 - date line.
         y(n) = Aea(l)

*        Compute spline coefficients
*        ---------------------------
         call SPLINE ( b, c, d, n, x, y )

*        Compute interpolated values
*        ---------------------------
         do i = 1, idim
            u = -180. + (i-1) * dlon
            All(i,j) =  SEVAL ( u, x, y, n, b, c, d )
         end do

      end do

*     Northern polar cap
*     ------------------
      do j = j_north+1, jdim-1                  ! no north pole

c-       ylat = -90. + (j-1)*dlat

         l = lea_beg(j)
         n = lea_len(j)

*        Copy q.e.a. data to local vector repeating first
*         longitudinal point to ensure periodicity
*        ------------------------------------------------
         do i = 1, n
            x(i) = ea_lon(l+i-1)
            y(i) = Aea(l+i-1)
         end do
         n = n + 1
         x(n) = +180.      ! -180 = + 180 - date line.
         y(n) = Aea(l)

*        Compute spline coefficients
*        ---------------------------
         call SPLINE ( b, c, d, n, x, y )

*        Compute interpolated values
*        ---------------------------
         do i = 1, idim
            u = -180. + (i-1) * dlon
            All(i,j) =  SEVAL ( u, x, y, n, b, c, d )
         end do

      end do

*     No interpolation needed at poles
*     --------------------------------
      do i = 1, idim
         All(i,1)    = Aea(i)                 ! south pole
         All(i,jdim) = Aea(ngrid-idim+i)      ! north pole
      end do

      deallocate(b,c,d,x,y)

      end

ext: f

max_stars_repo_name: GEOS-ESM/AeroApps

max_forks_repo_name: GEOS-ESM/AeroApps


text:
Address(El Capitan Street) is a residential street in West Davis.  It is named after the wiki:wikipedia: El Capitan granite rock formation in Yosemite National Park, which rises about 3,000 feet on its tallest face,

 Intersecting Streets

  Shasta Drive
  Tahoe Place
  Point Reyes Place
  Carlsbad Place
  Rushmore Lane
  White Sands Lane

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
My name is Brian.

Im a student, bartender, waiter, and lover of curry. I used to cohost The Civil Left (a political talkshow) on KDRT with Users/KenjiYamada.

I play Classical guitar, and now live in Sacramento but work in Davis.

Oi!

Its true.  Users/KenjiYamada

Brian, you should do an article for Thai Recipes!  Users/KenjiYamada

Brian licked our dusty fan once, it was quite funny.  Users/RohiniJasavala

20060731 11:07:18 nbsp Please refrain from http://daviswiki.org/urinating_while_sitting?actioninfo editing while shitfaced. At least I hope you were shitfaced. Users/ArlenAbraham


text:
       integer function ifeffit(string)
c
c//////////////////////////////////////////////////////////////////////
c Copyright (c) 1997--2000 Matthew Newville, The University of Chicago
c Copyright (c) 1992--1996 Matthew Newville, University of Washington
c
c Permission to use and redistribute the source code or binary forms of
c this software and its documentation, with or without modification is
c hereby granted provided that the above notice of copyright, these
c terms of use, and the disclaimer of warranty below appear in the
c source code and documentation, and that none of the names of The
c University of Chicago, The University of Washington, or the authors
c appear in advertising or endorsement of works derived from this
c software without specific prior written permission from all parties.
c
c THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
c EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
c MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
c IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
c CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
c TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
c SOFTWARE OR THE USE OR OTHER DEALINGS IN THIS SOFTWARE.
c//////////////////////////////////////////////////////////////////////
c
c purpose:    main ifeffit function.
c             parse and send string as a command to one of
c             the ifeffit subroutines
c arguments:  string    command line for ifeffit            [in]
c returns:    0   success
c            -2   in the middle of defining a macro
c            -1   get another line
c             1   normal exit
c            >1   failure of some sort
c
c notes:      see documentation for further discussion
c
c copyright (c) 1998  matt newville
c
c   Leaving there and proceeding for three days towards the east, you
c   reach Diomira, a city with sixty silver domes, bronze statues of
c   all the gods, streets paved with lead, a crystal theater, a golden
c   cock that crows each morning on a tower.  All these beauties will
c   already be familiar to the visitor, who has seen them also in
c   other cities.  But the special quality of this city for the man
c   who arrives there on a September evening, when the days are
c   growing shorter and the multicolored lamps are lighted all at once
c   at the doors of the food stalls and from a terrace a woman's voice
c   cries ooh!, is that he feels envy toward those who now believe
c   they have once before lived an evening identical to this and who
c   think they were happy, that time.
c
       implicit none
       include 'consts.h'
       include 'keywrd.h'
       include 'echo.h'
       save
       double precision getsca
       character*2048 line, argu, string*(*), val, key*128
       integer        ilen, ik1, ik2, init, istrln, jmac, iff_load
       integer        iffcmd_exec
       logical        wrthis, isamac
       external istrln, getsca, isamac, iff_load, iffcmd_exec
       data init /1/
c
c init>0 means to initialize
c  initial string may be used to alter behavior:
c   &quiet :   suppress all screen echo's, even startup message
c
       if (init .gt. 0) then
          init  = 0
          ilen  = 0
          jmac  = 0
          argu  = string
          call iff_init
          call echo_init
          call iff_plot_init(-1)
          call iff_config
          call triml(argu)
          if (argu(1:6) .eq. '&quiet') then
             call setsca( '&screen_echo',  zero)
             i_echo = 0
          endif
c
       end if
       line   = string
       argu   = ' '
       wrthis = .true.
cc       print*, 'str = ', line(1:40)
c
c  get next command line
c    getline converts blank lines and comment lines to "def",
c    so we check here that there really is an argument to
c    avoid repeated trivial processing
       ifeffit = -1
       call iff_getline(line, key, val, argu, ilen)
       ik1 = istrln(key)
       ik2 = istrln(argu)
cc       print*, ' =>', key(1:ik1), ' <', argu(1:ik2),'> ', ilen
       if ((key.eq.'def').and.(argu.eq.' ').and.
     $      (ilen.ge.0)) ilen = 0
c
c  if we're in the middle of a macro definition ...
       if (mac_define(1).ge.1) then
          call iff_macro_def(line)
cc          print*, ' setting ifeffit from ', ifeffit, ' to -2'
          ifeffit = -2
          call triml(line)
          call lower(line)
          if (line(1:9).eq.'end macro') then
             ifeffit = 0
             mac_define(1) = -1
          endif
          return
       endif
c
c  if ilen < 0, the line will be continued to the next line:
cc       print*, ' ilen = ', ilen
       if (ilen.eq.0) ifeffit = 0
       if (ilen.le.0) return
       ifeffit = 3
       iprint= int(getsca('&print_level',0))
c
c=== execute ifeffit commands based on keyword ===c
c

c
c first check for meta-commands that cannot occur within macros
c

cc       print*, 'IFF key = ', key(1:30), jmac, isamac(key,jmac)
c  execute a macro
       if (isamac(key,jmac)) then
cc          print*, ' jmac= ', isamac(key,jmac), jmac
          call iff_macro_do(jmac, argu, iprint,wrthis)
          wrthis = .false.
          key  = ' '
          argu = ' '

c  define a macro  (more precisly, __begin__ a macro definition)
       else
          ifeffit = iffcmd_exec(key,argu,wrthis)
       endif
c===  write history ===
c  and  information about argument line
       if (histry.and.wrthis)  then
          write(val, '(1x,4a)') key(1:ik1),' (',argu(1:ik2),')'
          if (iprint.ge.19)        call echo(val)
          ilen = istrln(val)
          write(iohist, '(a)') val(1:ilen)
       endif
c  clean up used line -- we're done
       line = '                                                   '//
     $      '                                                     '//
     $      '                                                     '//
     $      '                                                     '//
     $      '                                                     '

       call sclean(line)
cc       print*, 'ifeffit ret ', ifeffit
       if ((ifeffit.ne.-2).and.(ifeffit.ne.1)) ifeffit = 0

c end function ifeffit
       return
       end

       integer function iffcmd_exec(key,argu,wrthis)
       implicit none
       include 'consts.h'
       include 'keywrd.h'
       include 'echo.h'
       save

       character*2048 argu, key*128
       logical  wrthis

       integer  ilen, istrln, iff_load
       external istrln, iff_load

ccc       print*, ' iffcmd ', key(1:20)
       iffcmd_exec = 0
       if (key .eq. 'macro') then
          call iff_macro_def(argu)
          iffcmd_exec = -2
          wrthis = .false.
c  quit is a special case
       elseif ((key .eq. 'quit').or.(key .eq. 'exit')) then
          call iff_done
          iffcmd_exec = 1
          return
c  load: read a set of command lines from a file
       elseif (key .eq. 'load')  then
          iffcmd_exec = iff_load(argu)
          if (iffcmd_exec.eq.1) return
          wrthis = .false.
       elseif (key .eq. 'history')  then
          call iff_history(argu)
          wrthis = .false.

c  data i/o commands
       elseif (key .eq. '__test') then
          call iff_test(argu)
       elseif (key .eq. 'read_data') then
          call iff_rddata(argu)
       elseif (key .eq. 'write_data') then
          call iff_wrdata(argu)
       elseif (key .eq. 'save') then
          call iff_save(argu)
       elseif (key .eq. 'restore') then
          call iff_restore(argu)
c  set/show arrays and scalars
       elseif (key .eq. 'log') then
          call iff_log(argu)
       elseif (key .eq. 'show') then
          call iff_show(argu,0)
       elseif (key .eq. 'print') then
          call iff_print(argu)
       elseif (key .eq. 'echo') then
          call iff_echo(argu,0)
       elseif (key .eq. 'pause') then
          call iff_echo(argu,1)
       elseif ((key .eq. 'def').or.(key .eq. 'define').or.
     $         (key .eq. 'set').or.(key .eq. 'guess')) then
          call iff_set(key,argu,wrthis)
       elseif (key .eq. 'unguess') then
          call iff_unguess(argu)
       elseif (key .eq. 'erase') then
          call iff_erase(argu)
       elseif ((key .eq. 'rename')) then
          call iff_rename(argu)
c  data manipulation commands
       elseif ((key .eq. 'pre_edge')) then
          call iff_pre_edge(argu)
       elseif ((key .eq. 'bkg_cl')) then
          call iff_bkg_cl(argu)
       elseif ((key .eq. 'spline').or.(key .eq. 'autobk')) then
          call iff_spline(argu)
       elseif (key .eq. 'random') then
          call iff_random(argu)
       elseif ((key .eq. 'fft').or.(key .eq. 'fftf').or.
     $         (key .eq. 'bft').or.(key .eq. 'fftr')) then
          call iff_fft(key,argu)
       elseif (key .eq. 'window')  then
          call iff_window(argu)
       elseif (key .eq. 'path')  then
          call iff_path(argu)
       elseif (key .eq. 'get_path')  then
          call iff_getpath(argu)
       elseif (key .eq. 'ff2chi')  then
          call iff_ff2chi(argu)
       elseif (key .eq. 'feffit')  then
          call iff_feffit(argu)
       elseif (key .eq. 'chi_noise')  then
          call iff_chieps(argu)
       elseif (key .eq. 'f1f2')  then
          call iff_f1f2(argu)
       elseif (key .eq. 'diffkk')  then
          call iff_diffkk(argu)
       elseif (key .eq. 'minimize')  then
          call iff_minimize(argu)
       elseif (key .eq. 'correl')  then
          call iff_correl(argu)
       elseif (key .eq. 'error')  then
          call iff_uncert(argu)
       elseif (key .eq. 'sync')  then
          call iff_sync
       elseif (key .eq. 'reset') then
          call iff_init
          call echo_init
          call iff_plot_init(1)
          call iff_config
c
c  plotting commands
       elseif ((key .eq. 'plot').or.
     $         (key .eq. 'overplot').or.(key.eq.'replot'))  then
          call iff_plot(argu)
       elseif (key .eq. 'newplot') then
          call iff_plot('new=1, '//argu)
       elseif (key .eq. 'plot_marker') then
          call iff_plotmarker(argu)
       elseif (key .eq. 'plot_text') then
          call iff_plottext(argu)
       elseif (key .eq. 'plot_arrow') then
          call iff_plotarrow(argu)
       elseif (key .eq. 'zoom') then
          call iff_zoom(argu)
       elseif (key .eq. 'cursor') then
          call iff_cursor(argu)
       elseif (key .eq. 'color') then
          call iff_color(argu)
       elseif (key .eq. 'linestyle') then
          call iff_pstyle(argu)
c  comment is used to write argu to history file only as a comment
       elseif (key .eq. 'comment') then
          if (histry) then
             call rmquot(argu)
             ilen = istrln(argu)
             write(iohist, '(1x,2a)') '# '//argu(1:ilen)
          end if
          wrthis = .false.
       else
          ilen = istrln(key)
          call warn(1,'unknown command: '//key(1:ilen))
       endif
       return
       end

max_stars_repo_name: keechul/ifeffit

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_name: bruceravel/ifeffit

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_name: bruceravel/ifeffit

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


program test_spectral_interpolation_ipol_2

use kind_values, only:  wp => kind_swd_internal
use spectral_interpolation_def, only: spectral_interpolation

implicit none

! This program tests the ipol=2 temporal scheme

! Dummy polynomials: f(t) = a0 + a1*t + a2*t**2 + ... a7*t**7
complex(wp), parameter :: a0 = cmplx(0.156_wp, 0.17_wp, wp)
complex(wp), parameter :: a1 = cmplx(1.234_wp, -0.54_wp, wp)
complex(wp), parameter :: a2 = cmplx(-2.34_wp, 0.154_wp, wp)
complex(wp), parameter :: a3 = cmplx(-0.18_wp, -0.25_wp, wp)
complex(wp), parameter :: a4 = cmplx(0.58_wp, -0.74_wp, wp)
complex(wp), parameter :: a5 = cmplx(0.23_wp, 0.09_wp, wp)
complex(wp), parameter :: a6 = cmplx(-0.84_wp, 0.12_wp, wp)
complex(wp), parameter :: a7 = cmplx( 0.13_wp, -0.70_wp, wp)

real(wp), parameter :: delta_t = 0.34_wp
real(wp), parameter :: t_im1 = 0.15_wp
real(wp), parameter :: t_i = t_im1 + delta_t
real(wp), parameter :: t_ip1 = t_i + delta_t
real(wp), parameter :: t_ip2 = t_ip1 + delta_t

type(spectral_interpolation) :: si
real(wp) :: t, delta
complex(wp) :: c1, dc1, ddc1, c2, dc2, ddc2
complex(wp) :: f_im1, f_i, f_ip1, f_ip2
complex(wp) :: df_im1, df_i, df_ip1, df_ip2
complex(wp) :: ddf_i, ddf_ip1

integer :: ipol, i, ierr, n_error = 0
real(wp), parameter :: deltas(5) = [0.0_wp, 0.29_wp, 0.64_wp, 0.84_wp, 1.0_wp]

print*, " Test of ipol=2..."
ipol = 2
call si % construct(ipol, delta_t, ierr)
if (ierr /= 0) then
    print *, "Unknown interpolation scheme: ", ipol
    stop
end if

f_im1 = f7(t_im1)
f_i = f7(t_i)
f_ip1 = f7(t_ip1)
f_ip2 = f7(t_ip2)
df_im1 = df7(t_im1)
df_i = df7(t_i)
df_ip1 = df7(t_ip1)
df_ip2 = df7(t_ip2)
! Adjust f_im1 make ddf_i based on finite-difference match polynomial at t_i
! Adjust f_ip2 to make ddf_ip1 based on finite-difference match polynomial at t_ip1
stop 'Not implemented yet !!!'

ddf_ip1 = ddf7(t_ip1)
f_ip2 = 0.5_wp * ddf_ip1 * delta_t**2 + 2 * f_ip1 - f_i + (df_ip2 - df_i) * delta_t / 4

do i = 1, size(deltas)
    delta = deltas(i)
    t = t_i + delta * delta_t
    c1 = f7(t)
    dc1 = df7(t)
    ddc1 = ddf7(t)
    call si % scheme(delta, &
                     f_im1, f_i, f_ip1, f_ip2, &
                     df_im1, df_i, df_ip1, df_ip2, &
                     c2, dc2)
    print*, "Delta = ", delta
    print*, "Compare function values: ", c1, c2, compare(c1, c2)
    print*, "Compare function derivatives: ", dc1, dc2, compare(dc1, dc2)
    if (i==1) then
        print*, "f_i, df_i: ", f_i, df_i
        print*, "compare(c2, f_i)", compare(c2, f_i)
        print*, "compare(dc2, df_i)", compare(dc2, df_i)
        print*, "compare(ddc1, ddf_i)", compare(ddc1, ddf_i)
    else if (i==size(deltas)) then
        print*, "f_ip1, df_ip1: ", f_ip1, df_ip1
        print*, "compare(c2, f_ip1)", compare(c2, f_ip1)
        print*, "compare(dc2, df_ip1)", compare(dc2, df_ip1)
        print*, "compare(ddc1, ddf_ip1)", compare(ddc1, ddf_ip1)
    end if
end do

print*, "Number of detected ERRORs = ", n_error

contains

function f7(t) result(res)
real(wp), intent(in) :: t
complex(wp) res
res = a0 + a1*t + a2*t**2 + a3*t**3 + a4*t**4 + a5*t**5 + a6*t**6 + a7*t**7
end function f7

function df7(t) result(res)
real(wp), intent(in) :: t
complex(wp) res
res = a1 + 2.0_wp * a2*t + 3.0_wp * a3*t**2 + 4.0_wp * a4*t**3 + &
      5.0_wp * a5*t**4 + 6.0_wp * a6*t**5  + 7.0_wp * a7*t**6
end function df7

function ddf7(t) result(res)
real(wp), intent(in) :: t
complex(wp) res
res = 2.0_wp * a2 + 6.0_wp * a3*t + 12.0_wp * a4*t**2 + &
     20.0_wp * a5*t**3 + 30.0_wp * a6*t**4  + 42.0_wp * a7*t**5
end function ddf7

function compare(c1, c2) result(res)
complex(wp), intent(in) :: c1, c2
character(len=10) res
real(wp), parameter :: eps = 100*epsilon(eps)
if (abs(c1 % re - c2 % re) > eps * max(abs(c1 % re), abs(c2 % re)) .or. &
    abs(c1 % im - c2 % im) > eps * max(abs(c1 % im), abs(c2 % im))) then
    res = "ERROR!"
    n_error = n_error + 1
else
    res = "OK"
end if
end function compare

end program test_spectral_interpolation_ipol_2


text:
	SUBROUTINE SNHDTA  ( isnfln, vlevel, nlev, lvert, levtyp,
     +			     u, v,   npts,   utag, vtag,  tag,
     +			     ntag,   iret )
C************************************************************************
C* SNHDTA								*
C*									*
C* This subroutine gets data from the sounding file for SNHODO.		*
C*									*
C* SNHDTA  ( ISNFLN, VLEVEL, NLEV, LVERT, LEVTYP, U, V, NPTS,		*
C*            UTAG, VTAG,   TAG,    NTAG, IRET )				*
C*									*
C* Input parameters:							*
C*	ISNFLN		INTEGER		Sounding file number		*
C*	VLEVEL (NLEV)	REAL		Levels				*
C*	NLEV		INTEGER		Number of levels		*
C*	LVERT		INTEGER		Output vertical coordinate	*
C*	LEVTYP		INTEGER		Level type			*
C*					  1 = list			*
C*					  2 = range without increment	*
C*									*
C* Output parameters:							*
C*	U    (NPTS)	REAL		U for line			*
C*	V    (NPTS)	REAL		V for line			*
C*	NPTS		INTEGER		Number of line levels		*
C*	UTAG (NTAG)	REAL		U for tags			*
C*	VTAG (NTAG)	REAL		V for tags			*
C*	TAG  (NTAG)	REAL		Vertical coordinate tags	*
C*	NTAG		INTEGER		Number of tags			*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -9 = no data listed		*
C**									*
C* Log:									*
C* G. Huffman/USRA	 8/89	Adapted from SNLDTA, SNLLVD		*
C* S. Schotz/GSC	 7/90	Plot only mandatory and sig winds	*
C* S. Schotz/GSC	 8/90	Corrected error getting levels		*
C* K. Brill/NMC		 8/93	Change for 8-char ID			*
C* A. Hardy/GSC		 3/99	Added priority parameter to PC_SSTN     *
C* T. Piper/SAIC	 4/02	Fixed UMR; initialized selv, slat, slon	*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	REAL		vlevel (*), u (*), v (*), utag (*), vtag (*),
     +			tag (*)
	INTEGER		idtype (LLMXLV)
C*
	REAL		data ( LLMXDT ), rdata (3)
	CHARACTER	stntim*20, stid*8, cdata ( MMPARM )*8
C*
	INCLUDE		'ERMISS.FNC'
C------------------------------------------------------------------------
	iret = 0
	selv = 0.
	slat = 0.
	slon = 0.
C
C*	Get the data.
C
	CALL SN_RDAT  ( isnfln, ndlev, data, stntim, ier )
	IF  ( ier .ne. 0 )  THEN
	    iret = -9
	    CALL ER_WMSG  ( 'SN', ier, ' ', ier1 )
	    RETURN
	END IF
	CALL SN_RTYP ( isnfln, ndlev, idtype, iret )
C
	ispri = 0
	CALL PC_SSTN  ( stid,  istnm, slat, slon, selv, ispri, stntim,
     +			ndlev, ier )
C
C*	Loop through the levels and retrieve the winds.  They are stored
C*	twice, once for the line (U, V), and once for the vertical
C*	coordinate "tags".
C
C*	A list of levels to tag; set the bounds on the line ('ALL' is
C*	0 to -1) and do a separate retrieval for the tags.
C
	IF  ( levtyp .eq. 1 )  THEN
	    rstrt =  0.
	    rstop = -1.
	    ntag  = 0
	    DO  lev = 1, nlev
		CALL PC_CMVR  ( vlevel (lev), lvert, data, rdata,
     +				cdata, ier )
		IF  ( ( .not. ERMISS ( rdata (1) ) ) .and.
     +		      ( .not. ERMISS ( rdata (2) ) ) .and.
     +		      ( .not. ERMISS ( rdata (3) ) ) )  THEN
		    ntag       = ntag + 1
		    tag  (ntag) = rdata (1)
		    utag (ntag) = rdata (2)
		    vtag (ntag) = rdata (3)
		END IF
	    END DO
C
C*	  A range without increment of levels to tag; make the bounds
C*	  on the line search the same as the bounds on the tag search.
C
	  ELSE
	    rstrt = vlevel (1)
	    rstop = vlevel (2)
	END IF
C
C*	Find bottom and top levels for line (and tag range).
C
	CALL PC_FLVL  ( rstrt, lvert, data, rrr, lvnum1, lvnum2, itype,
     +								   ier )
	IF  ( itype .eq. 3 )  THEN
	    ilev = 1
	  ELSE IF  ( itype .eq. 4 )  THEN
	    ilev = nlev + 1
	  ELSE IF  ( itype .eq. 2 )  THEN
	    ilev = lvnum2
	  ELSE
	    ilev = lvnum1
	END IF
C
	CALL PC_FLVL  ( rstop, lvert, data, end, lvnum1, lvnum2, itype,
     +								   ier )
	IF  ( itype .eq. 3 )  THEN
	    jlev = 0
	  ELSE IF  ( itype .eq. 4 )  THEN
	    jlev = nlev
	  ELSE
	    jlev = lvnum1
	END IF
C
C*	Compute data at the intervening levels just for the line.
C
	IF  ( levtyp .eq. 1 )  THEN
	    npts = 0
	    DO  lev = ilev, jlev
		lll = lev - ilev + 1
		CALL PC_CMLV  ( lev, data, rdata, cdata, ier )
		IF  ( ( .not. ERMISS ( rdata (1) ) ) .and.
     +		      ( .not. ERMISS ( rdata (2) ) ) .and.
     +		      ( .not. ERMISS ( rdata (3) ) ) .and.
     +                ( idtype (lev) .ne. 2 ) )  THEN
		    npts     = npts + 1
		    u (npts) = rdata (2)
		    v (npts) = rdata (3)
		END IF
	    END DO
C
C*	  Compute data at the intervening levels for the line and
C*	  the tags.
C
	  ELSE
	    npts = 0
	    ntag = 0
	    DO  lev = ilev, jlev
		lll = lev - ilev + 1
		CALL PC_CMLV  ( lev, data, rdata, cdata, ier )
		IF  ( ( .not. ERMISS ( rdata (1) ) ) .and.
     +		      ( .not. ERMISS ( rdata (2) ) ) .and.
     +		      ( .not. ERMISS ( rdata (3) ) ) .and.
     +                ( idtype (lev) .ne. 2 ) )  THEN
		    npts        = npts + 1
		    u    (npts) = rdata (2)
		    v    (npts) = rdata (3)
		    ntag        = ntag + 1
		    tag  (ntag) = rdata (1)
		    utag (ntag) = rdata (2)
		    vtag (ntag) = rdata (3)
		END IF
	    END DO
	END IF
C
C*	Check that there are points to plot.
C
	IF  ( ( npts + ntag ) .eq. 0 )  THEN
	    iret = -9
	    CALL ER_WMSG  ( 'SNHODO', iret, ' ', ier )
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak