text:
program main

use m_precision
use m_constants_math
use m_trajectory
use m_globals
use m_config_io
use m_table
use m_qsort
use mkl

implicit none

type(itable_t) :: neighbor_table
type(itable_t) :: cluster_table
character(len=256) :: cla_buf
character(len=:), allocatable :: fn_sf, fn_probagg, fn_rg, fn_label
integer, dimension(:), allocatable :: label
real(rp), dimension(:,:), allocatable :: com_coordinates
real(rp), dimension(:,:), allocatable :: com_clusters
real(rp), dimension(:,:), allocatable :: molbuf
real(rp), dimension(:), allocatable :: r_bin
real(rp), dimension(:), allocatable :: gr_bin !Binned rdf
real(rp), dimension(:), allocatable :: normalization
! real(rp), dimension(:), allocatable :: prob_nn
real(rp), dimension(:), allocatable :: prob_aggnum
real(rp), dimension(:,:), allocatable :: dist_table
real(rp), dimension(:), allocatable :: kth_dist_table
real(rp), dimension(:), allocatable :: dist_buffer
integer, dimension(:), allocatable :: num_neighbor
real(rp), dimension(3) :: ri, rj, rij, com
real(rp), dimension(3) :: comi, comj, comij
real(rp) :: Rgsq
real(rp) :: rij_mag, rmax, rmin, rcrit
integer :: numframes
integer :: nchains
integer :: teq
integer :: istep, iatm, jatm, ibin, iframe, ierr, ios
integer :: num_solvobead, num_totalbead
integer :: fu_sf, fu_probagg, fu_rg, fu_label
integer :: ibr, ibeg, iend, ia_beg, ia_end
integer :: i, j, k, c
integer :: ichain, jchain
integer :: idx_k, minpts
integer :: num_nn
integer :: num_bins
integer :: cluster_count, ineighbor, neighbor_idx, nofn
integer :: agg_num, total_cluster_count
integer(ip_long) :: mode
real(rp) :: lbox
real(rp) :: bin_size
real(rp) :: rho
integer :: mflag
integer, dimension(:), pointer :: res

call get_command_argument(1, cla_buf)
fn_cfg = trim(adjustl(cla_buf))

call get_command_argument(2, cla_buf)
fn_traj = trim(adjustl(cla_buf))

call get_command_argument(3, cla_buf)
fn_sf = trim(adjustl(cla_buf))

call get_command_argument(4, cla_buf)
fn_probagg = trim(adjustl(cla_buf))

call get_command_argument(5, cla_buf)
fn_rg = trim(adjustl(cla_buf))

call get_command_argument(6, cla_buf)
teq = int(str_to_d(cla_buf))

call get_command_argument(7, cla_buf)
num_solvobead = int(str_to_d(cla_buf))

call get_command_argument(8, cla_buf)
num_totalbead = int(str_to_d(cla_buf))

call get_command_argument(9, cla_buf)
lbox = str_to_d(cla_buf)

call get_command_argument(10, cla_buf)
num_bins = str_to_d(cla_buf)

fn_label = 'label_all'

!Read final (can be initial as well) configuration file
call read_config(fn_cfg)
nchains = num_atoms/15
rmax = lbox/2
bin_size = rmax/num_bins

! allocate(r_bin(1:num_bins))
! allocate(gr_bin(1:num_bins))
! allocate(normalization(1:num_bins))

! allocate(prob_nn(1:100))

! allocate arrays for building distance table
allocate(dist_table(nchains,nchains))
! allocate arrays for building neighbor table
allocate(num_neighbor(1:nchains))
allocate(kth_dist_table(1:nchains))
allocate(dist_buffer(1:nchains))
! allocate arrays for recording cluster info
allocate(label(1:nchains))
allocate(prob_aggnum(1:nchains))
! allocate arrays for calculating RDF of clusters
allocate(r_bin(1:num_bins))
allocate(gr_bin(1:num_bins))
allocate(normalization(1:num_bins))

! initialize variables
numframes = 0
prob_aggnum = 0.0
total_cluster_count = 0
gr_bin = 0.0
r_bin = 0.0
normalization = 0.0
rmin = 0.0
rcrit = 4.2
idx_k = 4
minpts = 4

write(*,'(a,t15,g0.8)') 'BOXL', lbox
write(*,'(a,t15,g0.8)') 'RMAX', rmax
write(*,'(a,t15,g0.8)') 'NCHAINS', nchains
write(*,'(a,t15,g0.8)') 'BIN_SIZE', bin_size
write(*,'(a,t15,g0.8)') 'NUM_BINS', num_bins

allocate(com_coordinates(3,nchains))
allocate(com_clusters(3,nchains))
allocate( molbuf(3,num_solvobead) )

mflag = 1

!Open file for writing sf data & write header line
open(newunit=fu_sf, file=fn_sf, action='write', status='replace')

!Open file for Rg writing
open(newunit=fu_rg, file=fn_rg, action='write', status='replace')

!Open file for reading trajectory
call traj%open(fn_traj, 'r', ierr)
write(*,*) 'numframes ', traj%num_frames
do iframe = 1, traj%num_frames
    call traj%read(iframe, nts, ierr, mflag, coordinates)
    !write(*,*) 'nts ', nts
    if (nts < teq) cycle
    write(*,*) 'iframe ', iframe
    numframes = numframes + 1
    ! compress atoms within the same chain to the center-of-mass of
    ! the chain being checked
    do ichain = 1, nchains
        do iatm = 1, num_solvobead
            molbuf(:,iatm)= coordinates(:,num_totalbead*(ichain-1)+iatm)
        end do
        do i = 1, num_solvobead
            molbuf(:,i) = molbuf(:,i) - molbuf(:,1)
            call simbox%get_image(molbuf(:,i))
        end do
        com = sum(molbuf(:,1:num_solvobead), 2)/num_solvobead &
                            +coordinates(:,num_totalbead*(ichain-1)+1)
        !write(*,*) com
        com_coordinates(:,ichain) = com
    end do

    ! calculate g(r) for "com of chains", no longer used
    ! call simbox%wrap_all(com_coordinates)
    ! do ichain = 1, nchains-1
    !     ri = com_coordinates(:,ichain)
    !     do jchain = ichain+1, nchains
    !         rj = com_coordinates(:,jchain)
    !         rij = rj - ri
    !         rij(1) = rij(1) - nint(rij(1)/lbox)*lbox
    !         rij(2) = rij(2) - nint(rij(2)/lbox)*lbox
    !         rij(3) = rij(3) - nint(rij(3)/lbox)*lbox
    !         rij_mag = sqrt(rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3))
    !         if (rij_mag > rmax) cycle
    !         ibin = floor(rij_mag/bin_size) + 1
    !         !r_bin(ibin) = rij(ibin) + rij_mag
    !         gr_bin(ibin) = gr_bin(ibin) + 1
    !     end do
    ! end do

    kth_dist_table = 0.0
    dist_buffer = 0.0
    dist_table = rmax
    label = -1
    num_neighbor = 0
    ! build distance table for all atoms
    do ichain = 1, nchains - 1
        ri = com_coordinates(:,ichain)
        do jchain = ichain + 1, nchains
            rj = com_coordinates(:,jchain)
            rij = rj - ri
            rij(1) = rij(1) - nint(rij(1)/lbox)*lbox
            rij(2) = rij(2) - nint(rij(2)/lbox)*lbox
            rij(3) = rij(3) - nint(rij(3)/lbox)*lbox
            rij_mag = sqrt(rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3))
            dist_table(jchain,ichain) = rij_mag
            dist_table(ichain,jchain) = rij_mag
        end do
    end do

    ! ! Build the neighbor table. Get the k_th nearest neighbor for each
    ! ! "com of chain"(atom), and store the distance between this neighbor and
    ! ! the atom under consideration
    ! call itbl_init(neighbor_table, nchains)
    ! do ichain = 1, nchains
    !     ! Find neighbors
    !     do jchain = 1, nchains
    !         if (jchain == ichain) cycle
    !         if (dist_table(ichain,jchain) < rcrit) then
    !             num_neighbor(ichain) = num_neighbor(ichain) + 1
    !             call neighbor_table%append(ichain,jchain)
    !         end if
    !     end do
    !     ! sort the distance table to find the k_th neighbor of each atom
    !     dist_buffer = dist_table(:,ichain)
    !     !write(*,*) dist_buffer
    !     call dqsort(dist_buffer)
    !     kth_dist_table(ichain) = dist_buffer(idx_k)

    ! end do
    ! ! write(*,*) kth_dist_table

    ! create number of neighbor list
    call itbl_init(neighbor_table, nchains)
    call itbl_init(cluster_table, nchains)
    do ichain = 1, nchains
        ! Find neighbors
        do jchain = 1, nchains
            if (jchain == ichain) cycle
            if (dist_table(jchain,ichain) < rcrit) then
                num_neighbor(ichain) = num_neighbor(ichain) + 1
                call neighbor_table%append(ichain,jchain)
            end if
        end do
    end do
    ! write(*,*) num_neighbor
    ! call neighbor_table%print

    ! Cluster counter
    cluster_count = 0
    do ichain = 1, nchains
        ! Previously processed in inner loop
        if (label(ichain) /= -1) cycle
        ! Density check
        if (num_neighbor(ichain) < minPts) then
            ! Label as Noise
            label(ichain) = 0
            cycle
        end if
        ! next cluster label
        cluster_count = cluster_count + 1
        ! Label initial point
        label(ichain) = cluster_count
        ! Neighbors to expand
        do k = 1, num_neighbor(ichain)
            neighbor_idx = neighbor_table%get_val(ichain,k)
            call cluster_table%append(cluster_count,neighbor_idx)
        end do
        i = 1
        do while( .TRUE. )
            call cluster_table%get_row(cluster_count,res)
            if (size(res) < i) exit
            neighbor_idx = cluster_table%get_val(cluster_count,i)
            ! case if this neighbor is previously labeled as noise
            if (label(neighbor_idx) == 0) then
                ! change noise to border point
                label(neighbor_idx) = cluster_count
                ! Find neighbors of this neighbor
                if (num_neighbor(neighbor_idx) >= minpts) then
                    do j = 1, num_neighbor(neighbor_idx)
                        nofn = neighbor_table%get_val(neighbor_idx,j)
                        if (cluster_table%is_in(cluster_count,nofn) .OR. nofn == ichain) cycle
                        call cluster_table%append(cluster_count,nofn)
                    end do
                end if
            end if
            ! case if this neighbor is already in the cluster
            i = i + 1
            if (neighbor_idx == ichain) cycle
            ! case if this neighbor is unlabeled
            if (label(neighbor_idx) == -1) then
                ! Label neighbor
                label(neighbor_idx) = cluster_count
                ! Find neighbors of this neighbor
                if (num_neighbor(neighbor_idx) >= minpts) then
                    do j = 1, num_neighbor(neighbor_idx)
                        nofn = neighbor_table%get_val(neighbor_idx,j)
                        if (cluster_table%is_in(cluster_count,nofn) .OR. nofn == ichain) cycle
                        call cluster_table%append(cluster_count,nofn)
                    end do
                end if
            end if
        end do
    end do
    ! call cluster_table%print

    ! loop over all clusters, calculate statistics (aggregation number; Rg; etc.)
    do c = 1, cluster_count
        ! calculate aggregation number
        call cluster_table%get_row(c,res)
        agg_num = size(res)
        prob_aggnum(agg_num) = prob_aggnum(agg_num) + 1
        ! calculate COM of cluster based on COM's of chains
        com = 0.0
        do i = 1, agg_num
            ichain = cluster_table%get_val(c,i)
            com = com + com_coordinates(:,ichain)
        end do
        com = com/agg_num
        ! record COM of this cluster in an array
        com_clusters(:,c) = com
        ! calculate Rg
        Rgsq = 0.0
        do i = 1, agg_num
            ichain = cluster_table%get_val(c,i)
            do iatm = 1, num_solvobead
                molbuf(:,iatm)= coordinates(:,num_totalbead*(ichain-1)+iatm)
                ri = molbuf(:,iatm) - com
                Rgsq = Rgsq + ri(1)*ri(1) + ri(2)*ri(2) + ri(3)*ri(3)
            end do
        end do
        Rgsq = Rgsq/agg_num/num_solvobead
        write(fu_rg, '(es16.7)') Rgsq
    end do
    total_cluster_count = total_cluster_count + cluster_count

    ! loop over all clusters' COM, calculate RDF of clusters
    do c = 1, cluster_count - 1
        ri = com_clusters(:,c)
        do k = c + 1, cluster_count
            rj = com_clusters(:,k)
            rij = rj - ri
            rij(1) = rij(1) - nint(rij(1)/lbox)*lbox
            rij(2) = rij(2) - nint(rij(2)/lbox)*lbox
            rij(3) = rij(3) - nint(rij(3)/lbox)*lbox
            rij_mag = sqrt(rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3))
            if (rij_mag > rmax) cycle
            ibin = floor(rij_mag/bin_size) + 1
            gr_bin(ibin) = gr_bin(ibin) + 1
        end do
    end do

    open(newunit=fu_label, file=fn_label, action='write', status='replace')
    do i = 1, nchains
        do iatm = 1, num_solvobead
            k = (i-1)*num_solvobead+iatm
            rij = coordinates(:,num_totalbead*(i-1)+iatm)
            write(fu_label, '(i0,2x,i0,2x,*(g0.14,2x))') k, label(i)+1, rij
        end do
    end do
    close(fu_label)

    call neighbor_table%delete
    call cluster_table%delete

    ! probability function of neighbors, no longer used
    ! do ichain = 1, nchains
    !     ri = com_coordinates(:,ichain)
    !     num_nn = 0
    !     do jchain = 1, nchains
    !         if (jchain == ichain) cycle
    !         rj = com_coordinates(:,jchain)
    !         rij = rj - ri
    !         rij(1) = rij(1) - nint(rij(1)/lbox)*lbox
    !         rij(2) = rij(2) - nint(rij(2)/lbox)*lbox
    !         rij(3) = rij(3) - nint(rij(3)/lbox)*lbox
    !         rij_mag = sqrt(rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3))
    !         if (rij_mag < rcrit) then
    !             num_nn = num_nn + 1
    !         end if
    !     end do
    !     prob_nn(num_nn) = prob_nn(num_nn) + 1
    ! end do
    ! do i = 1, nchains
    !     !write(*,*) com_coordinates(:,i)
    !     write(fu_sf, '(i0,2x,i0,2x,*(g0.14,2x))') i, mflag, com_coordinates(:,i)
    ! end do
    ! write(*,*) prob_nn
end do

!gr_bin = gr_bin/normalization/numframes
!prob_nn = prob_nn/nchains/numframes
prob_aggnum = prob_aggnum/total_cluster_count

do i = 1, num_bins
    r_bin(i) = (i-1)*(rmax-rmin)/num_bins
end do

rho = total_cluster_count/numframes/lbox**3

do j = 1, num_bins-1
    normalization(j) = ((r_bin(j+1))**3-(r_bin(j))**3)*(4.0/3.0*math_pi*rho)*(total_cluster_count/numframes-1)/2
end do
normalization(num_bins) = (rmax**3-(rmax-bin_size)**3)*(4.0/3.0*math_pi*rho)*(total_cluster_count/numframes-1)/2

gr_bin = gr_bin/normalization/numframes
! write(*,*) 'gr_bin', gr_bin

write(fu_sf, '(es16.7)') rho
do i = 1, num_bins
    write(fu_sf, '(es16.7,2x,es16.7)') r_bin(i), gr_bin(i)
end do

! do i = 1, nchains
!     ! write(fu_sf, '(I4,2x,es16.7)') i, kth_dist_table(i)
!     write(fu_sf, '(i0,2x,i0,2x,*(g0.14,2x))') i, label(i)+1, com_coordinates(:,i)
! end do

!Close files
close(fu_sf)

close(fu_rg)

open(newunit=fu_probagg, file=fn_probagg, action='write', status='replace')
do i = 1, 100
    write(fu_probagg, '(I4,2x,es16.7)') i, prob_aggnum(i)
end do

close(fu_probagg)

!*******************************************************************************

end program

max_stars_repo_name: saridut/BROWNPAK

max_issues_repo_name: saridut/BROWNPAK

max_forks_repo_name: saridut/BROWNPAK


text:
!>  This module uses full-field binary wind files to determine the wind inflow.
!!  This module assumes that the origin, (0,0,0), is located at the tower centerline at ground level,
!!  and that all units are specified in the metric system (using meters and seconds).
!!  Data is assumed periodic in the X direction (and thus not shifted like FFWind files are).
MODULE IfW_HAWCWind
!!
!!  Created 25-June-2010 by B. Jonkman, National Renewable Energy Laboratory
!!     using subroutines and modules from AeroDyn v12.58
!!
!!----------------------------------------------------------------------------------------------------
!! Updated 8-Aug-2015 for InflowWind v3.0 in the FAST v8 Framework
!
!**********************************************************************************************************************************
! LICENSING
! Copyright (C) 2015-2016  National Renewable Energy Laboratory
!
!    This file is part of InflowWind.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.
!
!**********************************************************************************************************************************
   USE                                          NWTC_Library
   USE                                          IfW_HAWCWind_Types

   IMPLICIT                                     NONE
   PRIVATE

   TYPE(ProgDesc),   PARAMETER               :: IfW_HAWCWind_Ver = ProgDesc( 'IfW_HAWCWind', 'v1.01.00', '14-Dec-2015' )

   PUBLIC                                    :: IfW_HAWCWind_Init
   PUBLIC                                    :: IfW_HAWCWind_End
   PUBLIC                                    :: IfW_HAWCWind_CalcOutput

   INTEGER(IntKi), PARAMETER  :: nc = 3                           !< number of wind components
   INTEGER(IntKi), PARAMETER  :: WindProfileType_Constant = 0     !< constant wind
   INTEGER(IntKi), PARAMETER  :: WindProfileType_Log      = 1     !< logarithmic
   INTEGER(IntKi), PARAMETER  :: WindProfileType_PL       = 2     !< power law

   INTEGER(IntKi), PARAMETER  :: ScaleMethod_None         = 0     !< no scaling
   INTEGER(IntKi), PARAMETER  :: ScaleMethod_Direct       = 1     !< direct scaling factors
   INTEGER(IntKi), PARAMETER  :: ScaleMethod_StdDev       = 2     !< requested standard deviation

CONTAINS
!====================================================================================================
!>  This routine is used to initialize the parameters for using HAWC wind format files.
SUBROUTINE IfW_HAWCWind_Init(InitInp, p, MiscVars, Interval, InitOut, ErrStat, ErrMsg)

      ! Passed Variables
   TYPE(IfW_HAWCWind_InitInputType),         INTENT(IN   )  :: InitInp           !< Initialization data passed to the module
   TYPE(IfW_HAWCWind_ParameterType),         INTENT(  OUT)  :: p                 !< Parameters
   TYPE(IfW_HAWCWind_MiscVarType),           INTENT(  OUT)  :: MiscVars          !< Misc variables for optimization (not copied in glue code)
   TYPE(IfW_HAWCWind_InitOutputType),        INTENT(  OUT)  :: InitOut           !< Initialization output

   REAL(DbKi),                               INTENT(IN   )  :: Interval          !< Time Interval to use (passed through here)

      ! Error Handling
   INTEGER(IntKi),                           INTENT(  OUT)  :: ErrStat           !< determines if an error has been encountered
   CHARACTER(*),                             INTENT(  OUT)  :: ErrMsg            !< Message about errors

      ! Temporary variables for error handling
   INTEGER(IntKi)                                           :: TmpErrStat        ! temporary error status
   CHARACTER(ErrMsgLen)                                     :: TmpErrMsg         ! temporary error message
   CHARACTER(*),                             PARAMETER      :: RoutineName = 'IfW_HAWCWind_Init'

      ! Local Variables:

   ErrStat = ErrID_None
   ErrMsg  = ""

   ! validate init input data:
   call ValidateInput(InitInp, TmpErrStat, TmpErrMsg)
      CALL SetErrStat(TmpErrStat, TmpErrMsg, ErrStat, ErrMsg, RoutineName)
      IF (ErrStat >= AbortErrLev) RETURN

   !-------------------------------------------------------------------------------------------------
   ! Set some internal module parameters based on input file values
   !-------------------------------------------------------------------------------------------------
   p%nx           = InitInp%nx
   p%ny           = InitInp%ny
   p%nz           = InitInp%nz
   p%RefHt        = InitInp%RefHt
   p%URef         = InitInp%URef
   p%InitPosition = 0.0_ReKi  ! bjj: someday we may want to let the users give an offset time/position

   p%deltaXInv   = 1.0 / InitInp%dx
   p%deltaYInv   = 1.0 / InitInp%dy
   p%deltaZInv   = 1.0 / InitInp%dz

   p%LengthX     = InitInp%dx * p%nx !(nx-1)   !because the turbulence box is periodic in the X direction, we need to consider the length between point 1 and the next point 1 (instead of between points 1 and nx)
   p%LengthYHalf = 0.5*InitInp%dy * (p%ny-1)
   p%GridBase    = p%RefHt - 0.5*(p%nz-1)*InitInp%dz

   IF ( p%GridBase < 0.0_ReKi ) THEN
      call SetErrStat( ErrID_Fatal, 'The bottom of the grid is located at a height of '//&
                      TRIM( Num2LStr(p%GridBase) )//' meters, which is below the ground.', ErrStat,ErrMsg, RoutineName)
      RETURN
   END IF

   !-------------------------------------------------------------------------------------------------
   ! Read data files:
   !-------------------------------------------------------------------------------------------------
   call ReadTurbulenceData( p, InitInp, TmpErrStat, TmpErrMsg)
      CALL SetErrStat(TmpErrStat, TmpErrMsg, ErrStat, ErrMsg, RoutineName)
      IF (ErrStat >= AbortErrLev) RETURN

   !-------------------------------------------------------------------------------------------------
   ! scale to requested TI (or use requested scale factors)
   !-------------------------------------------------------------------------------------------------
   call ScaleTurbulence( p, InitInp, Interval, InitOut, MiscVars, TmpErrStat, TmpErrMsg)
      CALL SetErrStat(TmpErrStat, TmpErrMsg, ErrStat, ErrMsg, RoutineName)
      IF (ErrStat >= AbortErrLev) RETURN

   !-------------------------------------------------------------------------------------------------
   ! Add the mean wind speed to the u component.
   !-------------------------------------------------------------------------------------------------
   call AddMeanVelocity(p, InitInp, TmpErrStat, TmpErrMsg)
      CALL SetErrStat(TmpErrStat, TmpErrMsg, ErrStat, ErrMsg, RoutineName)
      IF (ErrStat >= AbortErrLev) RETURN

   !-------------------------------------------------------------------------------------------------
   ! write info to summary file, if necessary
   !-------------------------------------------------------------------------------------------------

   IF ( InitInp%SumFileUnit > 0 ) THEN
      WRITE(InitInp%SumFileUnit,'(A)', IOSTAT=TmpErrStat)
      WRITE(InitInp%SumFileUnit,'(A)', IOSTAT=TmpErrStat)    'HAWC wind type.  Read by InflowWind sub-module '//TRIM(GetNVD(IfW_HAWCWind_Ver))

      WRITE(InitInp%SumFileUnit,'(A34,G12.4)',IOSTAT=TmpErrStat)    '     Reference height (m):        ',p%RefHt
      WRITE(InitInp%SumFileUnit,'(A34,G12.4)',IOSTAT=TmpErrStat)    '     Timestep (s):                ',p%deltaXInv / p%URef
      WRITE(InitInp%SumFileUnit,'(A34,I12)',  IOSTAT=TmpErrStat)    '     Number of timesteps:         ',p%nx
      WRITE(InitInp%SumFileUnit,'(A34,G12.4)',IOSTAT=TmpErrStat)    '     Mean windspeed (m/s):        ',p%URef
      WRITE(InitInp%SumFileUnit,'(A)',        IOSTAT=TmpErrStat)    '     Time range (s):              [ '// &
                     TRIM(Num2LStr(0.0_ReKi))//' : '//TRIM(Num2LStr( p%LengthX / p%URef ))//' ]'
      WRITE(InitInp%SumFileUnit,'(A)',        IOSTAT=TmpErrStat)    '     X range (m):                 [ '// &
                     TRIM(Num2LStr(0.0_ReKi))//' : '//TRIM(Num2LStr( p%LengthX ))//' ]'
      WRITE(InitInp%SumFileUnit,'(A)',        IOSTAT=TmpErrStat)    '     Y range (m):                 [ '// &
                     TRIM(Num2LStr(-p%LengthYHalf))//' : '//TRIM(Num2LStr(p%LengthYHalf))//' ]'
      WRITE(InitInp%SumFileUnit,'(A)',        IOSTAT=TmpErrStat)    '     Z range (m):                 [ '// &
                     TRIM(Num2LStr(p%GridBase))//' : '//TRIM(Num2LStr(p%GridBase + p%nz / p%deltaZInv))//' ]'

      WRITE(InitInp%SumFileUnit,'(A)', IOSTAT=TmpErrStat)    'Scaling factors used:'
      WRITE(InitInp%SumFileUnit,'(A)', IOSTAT=TmpErrStat)    '  u           v           w       '
      WRITE(InitInp%SumFileUnit,'(A)', IOSTAT=TmpErrStat)    '----------  ----------  ----------'
      WRITE(InitInp%SumFileUnit,'(F10.3,2x,F10.3,2x,F10.3)',IOSTAT=TmpErrStat)   InitOut%sf
   ENDIF

   RETURN

END SUBROUTINE IfW_HAWCWind_Init
!====================================================================================================
!>  This routine is used to make sure the initInp data is valid.
SUBROUTINE ValidateInput(InitInp, ErrStat, ErrMsg)

      ! Passed Variables
   TYPE(IfW_HAWCWind_InitInputType),         INTENT(IN   )  :: InitInp           !< Initialization input data passed to the module

      ! Error Handling
   INTEGER(IntKi),                           INTENT(  OUT)  :: ErrStat           !< determines if an error has been encountered
   CHARACTER(*),                             INTENT(  OUT)  :: ErrMsg            !< Message about errors
   character(*), parameter                                  :: RoutineName = 'ValidateInput'

   integer(intki)                                           :: ic                ! loop counter

   ErrStat = ErrID_None
   ErrMsg  = ""

   IF ( InitInp%nx < 1 ) CALL SetErrStat( ErrID_Fatal, 'Number of grid points in the X direction must be at least 1.', ErrStat, ErrMsg, RoutineName )
   IF ( InitInp%ny < 1 ) CALL SetErrStat( ErrID_Fatal, 'Number of grid points in the Y direction must be at least 1.', ErrStat, ErrMsg, RoutineName )
   IF ( InitInp%nz < 1 ) CALL SetErrStat( ErrID_Fatal, 'Number of grid points in the Z direction must be at least 1.', ErrStat, ErrMsg, RoutineName )

   IF ( InitInp%dx < 0.0_ReKi .or. EqualRealNos( InitInp%dx, 0.0_ReKi ) ) CALL SetErrStat( ErrID_Fatal, 'The grid spacing in the X direction must be larger than 0.', ErrStat, ErrMsg, RoutineName )
   IF ( InitInp%dy < 0.0_ReKi .or. EqualRealNos( InitInp%dy, 0.0_ReKi ) ) CALL SetErrStat( ErrID_Fatal, 'The grid spacing in the Y direction must be larger than 0.', ErrStat, ErrMsg, RoutineName )
   IF ( InitInp%dz < 0.0_ReKi .or. EqualRealNos( InitInp%dz, 0.0_ReKi ) ) CALL SetErrStat( ErrID_Fatal, 'The grid spacing in the Z direction must be larger than 0.', ErrStat, ErrMsg, RoutineName )

   IF ( InitInp%RefHt < 0.0_ReKi .or. EqualRealNos( InitInp%RefHt, 0.0_ReKi ) ) call SetErrStat( ErrID_Fatal, 'The grid reference height must be larger than 0.', ErrStat, ErrMsg, RoutineName )

   if ( InitInp%ScaleMethod == ScaleMethod_Direct) then
      do ic=1,nc
         if ( InitInp%sf(ic) < 0.0_ReKi )  CALL SetErrStat( ErrID_Fatal, 'Turbulence scaling factors must not be negative.', ErrStat, ErrMsg, RoutineName )
      end do
   elseif ( InitInp%ScaleMethod == ScaleMethod_StdDev ) then
      do ic=1,nc
         if ( InitInp%sigmaf(ic) < 0.0_ReKi )  CALL SetErrStat( ErrID_Fatal, 'Turbulence standard deviations must not be negative.', ErrStat, ErrMsg, RoutineName )
      end do
#ifdef UNUSED_INPUTFILE_LINES
      if ( InitInp%TStart < 0.0_ReKi )  CALL SetErrStat( ErrID_Fatal, 'TStart for turbulence standard deviation calculations must not be negative.', ErrStat, ErrMsg, RoutineName )
      if ( InitInp%TEnd <= InitInp%TStart )  CALL SetErrStat( ErrID_Fatal, 'TEnd for turbulence standard deviation calculations must be after TStart.', ErrStat, ErrMsg, RoutineName )
#endif
   elseif ( InitInp%ScaleMethod /= ScaleMethod_None ) then
      CALL SetErrStat( ErrID_Fatal, 'HAWC scaling method must be 0 (none), 1 (direct scaling factors), or 2 (target standard deviation).', ErrStat, ErrMsg, RoutineName )
   end if

   if (InitInp%WindProfileType == WindProfileType_Log) then
      if ( InitInp%z0 < 0.0_ReKi .or. EqualRealNos( InitInp%z0, 0.0_ReKi ) ) &
         call SetErrStat( ErrID_Fatal, 'The surface roughness length, Z0, must be greater than zero', ErrStat, ErrMsg, RoutineName )
   elseif ( InitInp%WindProfileType < WindProfileType_Constant .or. InitInp%WindProfileType > WindProfileType_PL)  then
       call SetErrStat( ErrID_Fatal, 'The WindProfile type must be 0 (constant), 1 (logarithmic) or 2 (power law).', ErrStat, ErrMsg, RoutineName )
   end if

   IF ( InitInp%URef < 0.0_ReKi ) call SetErrStat( ErrID_Fatal, 'The reference wind speed must not be negative.', ErrStat, ErrMsg, RoutineName )

END SUBROUTINE ValidateInput
!====================================================================================================
!>  This routine is used read the full-field turbulence data stored in HAWC format.
SUBROUTINE ReadTurbulenceData(p, InitInp, ErrStat, ErrMsg)

      ! Passed Variables
   TYPE(IfW_HAWCWind_ParameterType),         INTENT(INOUT)  :: p                 !< Parameters
   TYPE(IfW_HAWCWind_InitInputType),         INTENT(IN   )  :: InitInp           !< Initialization input data passed to the module
   INTEGER(IntKi),                           INTENT(  OUT)  :: ErrStat           !< determines if an error has been encountered
   CHARACTER(*),                             INTENT(  OUT)  :: ErrMsg            !< Message about errors

      ! Local Variables:
   INTEGER                                                   :: IC               ! Loop counter for the number of wind components
   INTEGER                                                   :: IX               ! Loop counter for the number of grid points in the X direction
   INTEGER                                                   :: IY               ! Loop counter for the number of grid points in the Y direction
   INTEGER                                                   :: unWind           ! unit number for reading binary files

   INTEGER(IntKi)                                           :: TmpErrStat        ! temporary error status
   CHARACTER(ErrMsgLen)                                     :: TmpErrMsg         ! temporary error message
   CHARACTER(*),                             PARAMETER      :: RoutineName = 'ReadTurbulenceData'

   ErrStat = ErrID_None
   ErrMsg  = ""

   !-------------------------------------------------------------------------------------------------
   ! Allocate space for the wind array.
   !-------------------------------------------------------------------------------------------------

   CALL AllocAry( p%HAWCData, p%nz, p%ny, p%nx, nc, 'p%HAWCData', TmpErrStat, TmpErrMsg )
      CALL SetErrStat(TmpErrStat, TmpErrMsg, ErrStat, ErrMsg, RoutineName)
      IF (ErrStat >= AbortErrLev) RETURN

   !-------------------------------------------------------------------------------------------------
   ! Read the 3 files containg the turbulent wind speeds.
   !-------------------------------------------------------------------------------------------------
!bjj: check these indices... they do not seem to be very consistant between the WAsP IEC Turbulence
!     simulator and documentation of OC3 file formats... the current implementation is from the
!     OC3/Kenneth Thompson documentation.

      ! this could take a while, so we'll write a message indicating what's going on:

   CALL WrScr( NewLine//'   Reading HAWC wind files with grids of '//&
      TRIM( Num2LStr(p%nx) )//' x '//TRIM( Num2LStr(p%ny) )//' x '//TRIM( Num2LStr(p%nz) )//' points.' )

   CALL GetNewUnit( UnWind, TmpErrStat, TmpErrMsg )
      CALL SetErrStat(TmpErrStat, TmpErrMsg, ErrStat, ErrMsg, RoutineName)
      IF (ErrStat >= AbortErrLev) RETURN

      ! The array must be filled so that x(i) < x(i+1), y(i) < y(i+1), and z(i) < z(i+1)
      ! Also, note that the time axis is the negative x axis.

   DO IC = 1,NC

      CALL OpenBInpFile ( UnWind, InitInp%WindFileName(IC), TmpErrStat, TmpErrMsg )
         CALL SetErrStat(TmpErrStat, TmpErrMsg, ErrStat, ErrMsg, RoutineName)
         IF (ErrStat >= AbortErrLev) RETURN

      DO IX = p%nx,1,-1                  ! Time is the opposite of X ....
         DO IY = p%ny,1,-1
            !DO IZ = 1,p%nz

               READ( UnWind, IOSTAT=TmpErrStat ) p%HAWCData(:,iy,ix,ic)  ! note that HAWCData is SiKi (4-byte reals, not default kinds)

               IF (TmpErrStat /= 0) THEN
                  TmpErrMsg = ' Error reading binary data from "'//TRIM(InitInp%WindFileName(IC))//'". I/O error ' &
                                       //TRIM(Num2LStr(TmpErrStat))//' occurred at IY='//TRIM(Num2LStr(IY))//', IX='//TRIM(Num2LStr(IX))//'.'
                  CLOSE ( UnWind )
                  CALL SetErrStat(ErrID_Fatal, TmpErrMsg, ErrStat, ErrMsg, RoutineName)
                  RETURN
               END IF

            !END DO
         END DO
      END DO

      CLOSE ( UnWind )

   END DO

END SUBROUTINE ReadTurbulenceData
!====================================================================================================
!>  This routine is used read scale the full-field turbulence data stored in HAWC format.
SUBROUTINE ScaleTurbulence(p, InitInp, Interval, InitOut, MiscVars, ErrStat, ErrMsg)

      ! Passed Variables
   TYPE(IfW_HAWCWind_ParameterType),         INTENT(INOUT)  :: p                 !< Parameters
   TYPE(IfW_HAWCWind_InitInputType),         INTENT(IN   )  :: InitInp           !< Initialization input data passed to the module
   TYPE(IfW_HAWCWind_InitOutputType),        INTENT(INOUT)  :: InitOut           !< Initialization output data passed from the module
   REAL(DbKi),                               INTENT(IN   )  :: Interval          !< Time Interval to use (passed through here)
   TYPE(IfW_HAWCWind_MiscVarType),           INTENT(INOUT)  :: MiscVars          !< Misc variables for optimization (not copied in glue code)
   INTEGER(IntKi),                           INTENT(  OUT)  :: ErrStat           !< determines if an error has been encountered
   CHARACTER(*),                             INTENT(  OUT)  :: ErrMsg            !< Message about errors

      ! Local Variables:
   REAL(DbKi)                                               :: v(3)              ! instanteanous wind speed at target position
   REAL(DbKi)                                               :: vMean(3)          ! average wind speeds over time at target position
   REAL(DbKi)                                               :: vSum(3)           ! sum over time of wind speeds at target position
   REAL(DbKi)                                               :: vSum2(3)          ! sum of wind speeds squared
   REAL(ReKi)                                               :: ActualSigma(3)    ! computed standard deviation

   INTEGER                                                  :: ic                ! Loop counter for wind component
   INTEGER                                                  :: ix                ! Loop counter for x
   INTEGER                                                  :: iy                ! Loop counter for y
   INTEGER                                                  :: iz                ! Loop counter for z

   INTEGER(IntKi)                                           :: TmpErrStat        ! temporary error status
   CHARACTER(ErrMsgLen)                                     :: TmpErrMsg         ! temporary error message
   CHARACTER(*),                             PARAMETER      :: RoutineName = 'CalcScaleFactors'

   ErrStat = ErrID_None
   ErrMsg  = ""

   if ( InitInp%ScaleMethod == ScaleMethod_None ) then
      InitOut%sf = 1.0_ReKi
   else ! ScaleMethod_Direct or ScaleMethod_StdDev
      if ( InitInp%ScaleMethod == ScaleMethod_Direct ) then
         InitOut%sf = InitInp%sf
      else !if ( InitInp%ScaleMethod == ScaleMethod_StdDev ) then

#ifdef UNUSED_INPUTFILE_LINES
         ! calculate actual standard deviations, then compute scaling factor based on ratio of target to actual
         position(1) = 0.0_ReKi
         position(2) = 0.0_ReKi
         position(3) = p%RefHt

         vSum  = 0.0
         vSum2 = 0.0

         n = nint( ( InitInp%TEnd - InitInp%TStart ) / Interval )
         DO it=1,n
            t = InitInp%TStart + it*Interval
            v = FF_Interp(t,position,p,OtherStates,TmpErrStat,TmpErrMsg)
               CALL SetErrStat( TmpErrStat,TmpErrMsg, ErrStat, ErrMsg, RoutineName )

            vSum  = vSum  + v
            vSum2 = vSum2 + v**2
         ENDDO ! IT

         vMean = vSum/n
         ActualSigma = SQRT( ABS( (vSum2/n) - vMean**2 ) )

         !InitOut%sf = InitInp%SigmaF / ActualSigma  ! factor = Target / actual

         do ic=1,nc
            if ( EqualRealNos( ActualSigma(ic), 0.0_ReKi ) ) then
               InitOut%sf(ic) = 0.0_ReKi
               if ( .not. EqualRealNos( InitInp%SigmaF(ic), 0.0_ReKi ) ) then
                  call SetErrStat( ErrID_Fatal,"Computed standard deviation is zero; cannot scale to achieve target standard deviation.", ErrStat, ErrMsg, RoutineName )
               end if
            else
               InitOut%sf(ic) = InitInp%SigmaF(ic) / ActualSigma(ic)
            end if
         end do
#else

            ! roughly the point in the center of the grid
         iz = (p%nz + 1) / 2 ! integer division
         iy = (p%ny + 1) / 2 ! integer division

         DO ix=1,p%nx
            v = p%HAWCData(iz,iy,ix,:)

            vSum  = vSum  + v
            vSum2 = vSum2 + v**2
         ENDDO ! IT

         vMean = vSum/p%nx
         ActualSigma = SQRT( ABS( (vSum2/p%nx) - vMean**2 ) )

         !InitOut%sf = InitInp%SigmaF / ActualSigma  ! factor = Target / actual
         do ic=1,nc
            if ( EqualRealNos( ActualSigma(ic), 0.0_ReKi ) ) then
               InitOut%sf(ic) = 0.0_ReKi
               if ( .not. EqualRealNos( InitInp%SigmaF(ic), 0.0_ReKi ) ) then
                  call SetErrStat( ErrID_Fatal,"Computed standard deviation is zero; cannot scale to achieve target standard deviation.", ErrStat, ErrMsg, RoutineName )
               end if
            else
               InitOut%sf(ic) = InitInp%SigmaF(ic) / ActualSigma(ic)
            end if
         end do
#endif
      end if

         ! scale using our scaling factors:
      do ic = 1,nc
         p%HAWCData( :, :, :, ic ) = InitOut%sf(IC)*p%HAWCData( :, :, :, ic )
      end do !IC

   end if

END SUBROUTINE ScaleTurbulence
!====================================================================================================
!>  This routine is used to add a mean wind profile to the HAWC format turbulence data.
SUBROUTINE AddMeanVelocity(p, InitInp, ErrStat, ErrMsg)

      ! Passed Variables
   TYPE(IfW_HAWCWind_ParameterType),         INTENT(INOUT)  :: p                 !< Parameters
   TYPE(IfW_HAWCWind_InitInputType),         INTENT(IN   )  :: InitInp           !< Initialization input data passed to the module
   INTEGER(IntKi),                           INTENT(  OUT)  :: ErrStat           !< determines if an error has been encountered
   CHARACTER(*),                             INTENT(  OUT)  :: ErrMsg            !< Message about errors

      ! Local Variables:
   REAL(ReKi)                                               :: Z                 ! height
   REAL(ReKi)                                               :: U                 ! mean wind speed
   INTEGER(IntKi)                                           :: iz                ! loop counter
   INTEGER(IntKi)                                           :: TmpErrStat        ! temporary error status
   CHARACTER(ErrMsgLen)                                     :: TmpErrMsg         ! temporary error message
   CHARACTER(*),                             PARAMETER      :: RoutineName = 'AddMeanVelocity'

   ErrStat = ErrID_None
   ErrMsg  = ""

   DO IZ = 1,p%NZ

      Z = p%GridBase  + ( IZ - 1 )*InitInp%dz

      SELECT CASE ( InitInp%WindProfileType )

      CASE ( WindProfileType_PL )

            U = p%URef*( Z / p%RefHt )**InitInp%PLExp      ! [IEC 61400-1 6.3.1.2 (10)]

      CASE ( WindProfileType_Log )

            IF ( .not. EqualRealNos( p%RefHt, InitInp%Z0 ) .and. z > 0.0_ReKi ) THEN
               U = p%URef*( LOG( Z / InitInp%Z0 ) )/( LOG( p%RefHt / InitInp%Z0 ) )
            ELSE
               U = 0.0_ReKi
            ENDIF

      CASE ( WindProfileType_Constant )

           U = p%URef

      CASE DEFAULT

            U = 0.0_ReKi

      END SELECT

      p%HAWCData( IZ, :, :, 1 ) = p%HAWCData( IZ, :, :, 1 ) + U

   END DO ! IZ

END SUBROUTINE AddMeanVelocity
!====================================================================================================
!> This routine acts as a wrapper for the GetWindSpeed routine. It steps through the array of input
!! positions and calls the GetWindSpeed routine to calculate the velocities at each point.
!!
!! There are inefficiencies in how this set of routines is coded, but that is a problem for another
!! day. For now, it merely needs to be functional. It can be fixed up and made all pretty later.
!!
!!   16-Apr-2013 - A. Platt, NREL.  Converted to modular framework. Modified for NWTC_Library 2.0
SUBROUTINE IfW_HAWCWind_CalcOutput(Time, PositionXYZ, p, Velocity, DiskVel, MiscVars, ErrStat, ErrMsg)

   IMPLICIT NONE

      ! Passed Variables
   REAL(DbKi),                               INTENT(IN   )  :: Time              !< time from the start of the simulation
   REAL(ReKi),                               INTENT(IN   )  :: PositionXYZ(:,:)  !< Array of XYZ coordinates, 3xN
   TYPE(IfW_HAWCWind_ParameterType),         INTENT(IN   )  :: p                 !< Parameters
   REAL(ReKi),                               INTENT(INOUT)  :: Velocity(:,:)     !< Velocity output at Time    (Set to INOUT so that array does not get deallocated)
   REAL(ReKi),                               INTENT(  OUT)  :: DiskVel(3)        !< HACK for AD14: disk velocity output at Time
   TYPE(IfW_HAWCWind_MiscVarType),           INTENT(INOUT)  :: MiscVars          !< Misc variables for optimization (not copied in glue code)

      ! Error handling
   INTEGER(IntKi),                           INTENT(  OUT)  :: ErrStat           !< error status
   CHARACTER(*),                             INTENT(  OUT)  :: ErrMsg            !< The error message

      ! local variables
   INTEGER(IntKi)                                           :: NumPoints      ! Number of points specified by the PositionXYZ array

      ! local counters
   INTEGER(IntKi)                                           :: PointNum       ! a loop counter for the current point

      ! temporary variables
   INTEGER(IntKi)                                           :: TmpErrStat     ! temporary error status
   CHARACTER(ErrMsgLen)                                     :: TmpErrMsg      ! temporary error message
   CHARACTER(*), PARAMETER                                  :: RoutineName = 'IfW_HAWCWind_CalcOutput'

      !-------------------------------------------------------------------------------------------------
      ! Check that the module has been initialized.
      !-------------------------------------------------------------------------------------------------

   ErrStat     = ErrID_None
   ErrMsg      = ''

      !-------------------------------------------------------------------------------------------------
      ! Initialize some things
      !-------------------------------------------------------------------------------------------------

      ! The array is transposed so that the number of points is the second index, x/y/z is the first.
      ! This is just in case we only have a single point, the SIZE command returns the correct number of points.
   NumPoints   =  SIZE(PositionXYZ,2)

      ! Step through all the positions and get the velocities
   DO PointNum = 1, NumPoints

         ! Calculate the velocity for the position
      Velocity(:,PointNum) = FF_Interp(Time,PositionXYZ(:,PointNum),p,MiscVars,TmpErrStat,TmpErrMsg)

         ! Error handling
      IF (TmpErrStat /= ErrID_None) THEN  !  adding this so we don't have to convert numbers to strings every time
         CALL SetErrStat( TmpErrStat, TmpErrMsg, ErrStat, ErrMsg, RoutineName//" [position=("//   &
                                                      TRIM(Num2LStr(PositionXYZ(1,PointNum)))//", "// &
                                                      TRIM(Num2LStr(PositionXYZ(2,PointNum)))//", "// &
                                                      TRIM(Num2LStr(PositionXYZ(3,PointNum)))//") in wind-file coordinates]" )
         IF (ErrStat >= AbortErrLev) RETURN
      END IF

   ENDDO

      !REMOVE THIS for AeroDyn 15
      ! Return the average disk velocity values needed by AeroDyn 14.  This is the WindInf_ADhack_diskVel routine.
   DiskVel(1)   =  p%URef
   DiskVel(2:3) =  0.0_ReKi

   RETURN

END SUBROUTINE IfW_HAWCWind_CalcOutput
!====================================================================================================
   !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
   !>    This function is used to interpolate into the full-field wind array.  It receives X, Y, Z and
   !!    TIME from the calling routine.  It then computes a time shift due to a nonzero X based upon
   !!    the average windspeed.  The modified time is used to decide which pair of time slices to interpolate
   !!    within and between.  After finding the two time slices, it decides which four grid points bound the
   !!    (Y,Z) pair.  It does a 3-d linear interpolation.  This routine assumes that X is downwind, Y is to the left when
   !!    looking downwind and Z is up.  It also assumes that no extrapolation will be needed.
   FUNCTION FF_Interp(Time, Position, p, MiscVars, ErrStat, ErrMsg)

      REAL(DbKi),                         INTENT(IN   )  :: Time           !< Time at which to find wind speed
      REAL(ReKi),                         INTENT(IN   )  :: Position(3)    !< takes the place of XGrnd, YGrnd, ZGrnd
      TYPE(IfW_HAWCWind_ParameterType),   INTENT(IN   )  :: p              !< Parameters
      TYPE(IfW_HAWCWind_MiscVarType),     INTENT(INOUT)  :: MiscVars       !< Misc variables for optimization (not copied in glue code)
      REAL(ReKi)                                         :: FF_Interp(3)   !< The U, V, W velocities

      INTEGER(IntKi),                     INTENT(  OUT)  :: ErrStat        !< error status
      CHARACTER(*),                       INTENT(  OUT)  :: ErrMsg         !< error message

         ! Local Variables:

      CHARACTER(*),           PARAMETER                  :: RoutineName="FF_Interp"

      REAL(ReKi)                                         :: ShiftedXPosition
      REAL(ReKi),PARAMETER                               :: Tol = 1.0E-3   ! a tolerance for determining if two reals are the same (for extrapolation)
      REAL(ReKi)                                         :: X              ! value between -1 and 1 (where we are between 2 x grid points)
      REAL(ReKi)                                         :: XGRID
      REAL(ReKi)                                         :: Y
      REAL(ReKi)                                         :: YGRID
      REAL(ReKi)                                         :: Z
      REAL(ReKi)                                         :: ZGRID
      REAL(ReKi)                                         :: N(8)           ! array for holding scaling factors for the interpolation algorithm
      REAL(ReKi)                                         :: u(8)           ! array for holding the corner values for the interpolation algorithm across a cubic volume

      INTEGER(IntKi)                                     :: IDIM
      INTEGER(IntKi)                                     :: IXHI
      INTEGER(IntKi)                                     :: IXLO
      INTEGER(IntKi)                                     :: IYHI
      INTEGER(IntKi)                                     :: IYLO
      INTEGER(IntKi)                                     :: IZHI
      INTEGER(IntKi)                                     :: IZLO

      !-------------------------------------------------------------------------------------------------
      ! Initialize variables
      !-------------------------------------------------------------------------------------------------

      FF_Interp(:)         = 0.0_ReKi                         ! the output velocities (in case p%NFFComp /= 3)

      ErrStat              = ErrID_None
      ErrMsg               = ""

      !-------------------------------------------------------------------------------------------------
      ! Find the bounding time slices.
      !-------------------------------------------------------------------------------------------------

   ! bjj: should we shift by MIN(YHalfWid,FFZHWid)?

            ! Assume Taylor's Frozen Turbulence Hypothesis applies: u(X,Y,Z,t) = u( X-U*t, Y, Z, 0)

      ShiftedXPosition = Position(1) - TIME*p%URef - p%InitPosition(1)      !this puts the first X grid point at the undeflected tower centerline (or p%InitXPosition)

         ! The wind file is periodic so we'll translate this position to ( 0 <= ShiftedXPosition < p%LengthX )

      ShiftedXPosition = MODULO( ShiftedXPosition, p%LengthX )
      ! If ShiftedXPosition is a very small negative number, modulo returns the incorrect value due to internal rounding errors.
      ! See bug report #471
      IF (ShiftedXPosition == p%LengthX) ShiftedXPosition = 0.0_ReKi

      XGrid = ShiftedXPosition * p%deltaXInv

      IXLO = INT( XGrid )              ! convert REAL to INTEGER
      X = 2.0_ReKi * ( XGrid - REAL(IXLO, ReKi) ) - 1.0_ReKi     ! a value between -1 and 1 that indicates a relative position between IXLO and IXHI
      IXLO = IXLO + 1  ! Add 1 because our grids start at 1, not zero

      IF ( IXLO == p%NX ) THEN
         IXHI = 1
      ELSE
         IXHI = IXLO + 1
      ENDIF

      !-------------------------------------------------------------------------------------------------
      ! Find the bounding rows for the Z position. [The lower-left corner is (1,1) when looking upwind.]
      !-------------------------------------------------------------------------------------------------

      ZGRID = ( Position(3) - p%GridBase - p%InitPosition(3) )*p%deltaZInv

         ! Index for start and end slices
      IZLO = INT( ZGRID ) + 1             ! convert REAL to INTEGER, then add one since our grids start at 1, not 0
      IZHI = IZLO + 1

         ! Set Z as a value between -1 and 1 for the relative location between IZLO and IZHI.
         ! Subtract 1_IntKi from Z since the indices are starting at 1, not 0
      Z = 2.0_ReKi * (ZGRID - REAL(IZLO - 1_IntKi, ReKi)) - 1.0_ReKi

      IF ( IZLO < 1 ) THEN
         IF ( IZLO == 0 .AND. Z >= 1.0-TOL ) THEN
            Z    = -1.0_ReKi
            IZLO = 1
         ELSE
            ErrMsg   = ' FF wind array boundaries violated. Grid too small in Z direction (Z='//&
                        TRIM(Num2LStr(Position(3)))//' m is below the grid).'
            ErrStat  = ErrID_Fatal
            RETURN
         ENDIF
      ELSEIF ( IZLO >= p%nz ) THEN
         IF ( IZLO == p%nz .AND. Z <= TOL ) THEN
            Z    = -1.0_ReKi
            IZHI = IZLO                   ! We're right on the last point, which is still okay
         ELSE
            ErrMsg   = ' FF wind array boundaries violated. Grid too small in Z direction (Z='//&
                        TRIM(Num2LStr(Position(3)))//' m is above the grid).'
            ErrStat  = ErrID_Fatal
            RETURN
         ENDIF
      ENDIF

      !-------------------------------------------------------------------------------------------------
      ! Find the bounding columns for the Y position. [The lower-left corner is (1,1) when looking upwind.]
      !-------------------------------------------------------------------------------------------------
      YGRID = ( Position(2) + p%LengthYHalf - p%InitPosition(2) )*p%deltaYInv    ! really, it's (Position(2) - -1.0*p%LengthYHalf)

      IYLO = INT( YGRID ) + 1             ! convert REAL to INTEGER, then add one since our grids start at 1, not 0
      IYHI = IYLO + 1

         ! Set Y as a value between -1 and 1 for the relative location between IYLO and IYHI.  Used in the interpolation.
         ! Subtract 1_IntKi from IYLO since grids start at index 1, not 0
      Y = 2.0_ReKi * (YGRID - REAL(IYLO - 1_IntKi, ReKi)) - 1.0_ReKi

      IF ( IYLO >= p%ny .OR. IYLO < 1 ) THEN
         IF ( IYLO == 0 .AND. Y >= 1.0-TOL ) THEN
            Y    = -1.0_ReKi
            IYLO = 1
         ELSE IF ( IYLO == p%ny .AND. Y <= TOL ) THEN
            Y    = -1.0_ReKi
            IYHI = IYLO                   ! We're right on the last point, which is still okay
         ELSE
            ErrMsg   = ' FF wind array boundaries violated: Grid too small in Y direction. Y='// &
                        TRIM(Num2LStr(Position(2)))//'; Y boundaries = ['//TRIM(Num2LStr(-1.0*p%LengthYHalf))// &
                        ', '//TRIM(Num2LStr(p%LengthYHalf))//']'
            ErrStat = ErrID_Fatal         ! we don't return anything
            RETURN
         ENDIF
      ENDIF

      !-------------------------------------------------------------------------------------------------
      ! Interpolate on the grid
      !-------------------------------------------------------------------------------------------------

      DO IDIM=1,nc       ! all the components

         N(1)  = ( 1.0_ReKi + Z )*( 1.0_ReKi - Y )*( 1.0_ReKi - X )
         N(2)  = ( 1.0_ReKi + Z )*( 1.0_ReKi + Y )*( 1.0_ReKi - X )
         N(3)  = ( 1.0_ReKi - Z )*( 1.0_ReKi + Y )*( 1.0_ReKi - X )
         N(4)  = ( 1.0_ReKi - Z )*( 1.0_ReKi - Y )*( 1.0_ReKi - X )
         N(5)  = ( 1.0_ReKi + Z )*( 1.0_ReKi - Y )*( 1.0_ReKi + X )
         N(6)  = ( 1.0_ReKi + Z )*( 1.0_ReKi + Y )*( 1.0_ReKi + X )
         N(7)  = ( 1.0_ReKi - Z )*( 1.0_ReKi + Y )*( 1.0_ReKi + X )
         N(8)  = ( 1.0_ReKi - Z )*( 1.0_ReKi - Y )*( 1.0_ReKi + X )
         N     = N / REAL( SIZE(N), ReKi )  ! normalize

         u(1)  = p%HAWCData( IZHI, IYLO, IXLO, IDIM )
         u(2)  = p%HAWCData( IZHI, IYHI, IXLO, IDIM )
         u(3)  = p%HAWCData( IZLO, IYHI, IXLO, IDIM )
         u(4)  = p%HAWCData( IZLO, IYLO, IXLO, IDIM )
         u(5)  = p%HAWCData( IZHI, IYLO, IXHI, IDIM )
         u(6)  = p%HAWCData( IZHI, IYHI, IXHI, IDIM )
         u(7)  = p%HAWCData( IZLO, IYHI, IXHI, IDIM )
         u(8)  = p%HAWCData( IZLO, IYLO, IXHI, IDIM )

         FF_Interp(IDIM)  =  SUM ( N * u )

      END DO !IDIM

      RETURN

   END FUNCTION FF_Interp
!====================================================================================================
!>  This subroutine cleans up any data that is still allocated.  The (possibly) open files are
!!  closed in InflowWindMod.
SUBROUTINE IfW_HAWCWind_End( p, MiscVars, ErrStat, ErrMsg)

   IMPLICIT                                                 NONE

   CHARACTER(*),           PARAMETER                     :: RoutineName="IfW_HAWCWind_End"

      ! Passed Variables
   TYPE(IfW_HAWCWind_ParameterType),      INTENT(INOUT)  :: p                 !< Parameters
   TYPE(IfW_HAWCWind_MiscVarType),        INTENT(  OUT)  :: MiscVars          !< Misc variables for optimization (not copied in glue code)

      ! Error Handling
   INTEGER(IntKi),                        INTENT(  OUT)  :: ErrStat           !< determines if an error has been encountered
   CHARACTER(*),                          INTENT(  OUT)  :: ErrMsg            !< Message about errors

      ! Local Variables
   INTEGER(IntKi)                                        :: TmpErrStat     ! temporary error status
   CHARACTER(ErrMsgLen)                                  :: TmpErrMsg      ! temporary error message

      !-=- Initialize the routine -=-

   ErrMsg   = ''
   ErrStat  = ErrID_None

      ! Destroy parameter data

   CALL IfW_HAWCWind_DestroyParam(       p,     TmpErrStat, TmpErrMsg )
   CALL SetErrStat( TmpErrStat, TmpErrMsg, ErrStat, ErrMsg, RoutineName)

      ! Destroy the state data

   CALL IfW_HAWCWind_DestroyMisc(  MiscVars,   TmpErrStat, TmpErrMsg )
   CALL SetErrStat( TmpErrStat, TmpErrMsg, ErrStat, ErrMsg, RoutineName)

END SUBROUTINE IfW_HAWCWind_End

!====================================================================================================
END MODULE IfW_HAWCWind

max_stars_repo_name: sayerhs/OpenFAST

max_issues_repo_name: sayerhs/OpenFAST

max_forks_repo_name: sayerhs/OpenFAST


text:
module fitting
contains

  function polyfit(vx, vy, d)
    implicit none
    integer, intent(in)                   :: d
    integer, parameter                    :: dp = selected_real_kind(15, 307)
    real(dp), dimension(d+1)              :: polyfit
    real(dp), dimension(:), intent(in)    :: vx, vy

    real(dp), dimension(:,:), allocatable :: X
    real(dp), dimension(:,:), allocatable :: XT
    real(dp), dimension(:,:), allocatable :: XTX

    integer :: i, j

    integer     :: n, lda, lwork
    integer :: info
    integer, dimension(:), allocatable :: ipiv
    real(dp), dimension(:), allocatable :: work

    n = d+1
    lda = n
    lwork = n

    allocate(ipiv(n))
    allocate(work(lwork))
    allocate(XT(n, size(vx)))
    allocate(X(size(vx), n))
    allocate(XTX(n, n))

    ! prepare the matrix
    do i = 0, d
       do j = 1, size(vx)
          X(j, i+1) = vx(j)**i
       end do
    end do

    XT  = transpose(X)
    XTX = matmul(XT, X)

    ! calls to LAPACK subs DGETRF and DGETRI
    call DGETRF(n, n, XTX, lda, ipiv, info)
    if ( info /= 0 ) then
       print *, "problem"
       return
    end if
    call DGETRI(n, XTX, lda, ipiv, work, lwork, info)
    if ( info /= 0 ) then
       print *, "problem"
       return
    end if

    polyfit = matmul( matmul(XTX, XT), vy)

    deallocate(ipiv)
    deallocate(work)
    deallocate(X)
    deallocate(XT)
    deallocate(XTX)

  end function

end module

ext: f

max_stars_repo_name: mullikine/RosettaCodeData

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
module neutrino_lib
use neutrino_def
use superfluid_def, only: max_number_sf_types

contains

subroutine get_core_neutrino_emissivity(nn,np,T,Tcs,use_mode,epsilon_nu)
    use eval_core_neutrino
    real(dp), intent(in) :: nn,np,T,Tcs(max_number_sf_types)
    logical, intent(in), dimension(num_core_nu_channels) :: use_mode
    type(core_neutrino_emissivity_channels), intent(out) :: epsilon_nu

    call emissivity(nn,np,T,Tcs,use_mode,epsilon_nu)
end subroutine get_core_neutrino_emissivity

subroutine get_crust_neutrino_emissivity(rho,T,ionic,chi,Tcn,epsilon,use_mode)
    use nucchem_def, only: composition_info_type
    use eval_crust_neutrino
    real(dp), intent(in) :: rho,T,chi
    type(composition_info_type), intent(in) :: ionic
    real(dp), intent(in) :: Tcn ! neutron critical tempererature
    type(crust_neutrino_emissivity_channels), intent(out) :: epsilon
    logical, intent(in), dimension(num_crust_nu_channels), optional :: use_mode
    logical, dimension(num_crust_nu_channels) :: channels

    channels = crust_nu_full_cooling
    if (present(use_mode)) channels = use_mode

    call emissivity(rho,T,ionic,chi,Tcn,channels,epsilon)
end subroutine get_crust_neutrino_emissivity

end module neutrino_lib

max_stars_repo_name: LBJ-Wade/dStar

max_issues_repo_name: LBJ-Wade/dStar

max_forks_repo_name: LBJ-Wade/dStar


text:
module sgs

! module for original SAM subgrid-scale SGS closure (Smagorinsky or 1st-order TKE)
! Marat Khairoutdinov, 2012

use grid, only: nx,nxp1,ny,nyp1,YES3D,nzm,nz,dimx1_s,dimx2_s,dimy1_s,dimy2_s
use params, only: dosgs
implicit none

!----------------------------------------------------------------------
! Required definitions:

!!! prognostic scalar (need to be advected arround the grid):

integer, parameter :: nsgs_fields = 1   ! total number of prognostic sgs vars

real sgs_field(dimx1_s:dimx2_s, dimy1_s:dimy2_s, nzm, nsgs_fields)

!!! sgs diagnostic variables that need to exchange boundary information (via MPI):

integer, parameter :: nsgs_fields_diag = 2   ! total number of diagnostic sgs vars

! diagnostic fields' boundaries:
integer, parameter :: dimx1_d=0, dimx2_d=nxp1, dimy1_d=1-YES3D, dimy2_d=nyp1

real sgs_field_diag(dimx1_d:dimx2_d, dimy1_d:dimy2_d, nzm, nsgs_fields_diag)

logical:: advect_sgs = .false. ! advect prognostics or not, default - not (Smagorinsky)
logical, parameter:: do_sgsdiag_bound = .true.  ! exchange boundaries for diagnostics fields

! SGS fields that output by default (if =1).
integer, parameter :: flag_sgs3Dout(nsgs_fields) = (/0/)
integer, parameter :: flag_sgsdiag3Dout(nsgs_fields_diag) = (/0,0/)

real fluxbsgs (nx, ny, 1:nsgs_fields) ! surface fluxes
real fluxtsgs (nx, ny, 1:nsgs_fields) ! top boundary fluxes

!!! these arrays may be needed for output statistics:

real sgswle(nz,1:nsgs_fields)  ! resolved vertical flux
real sgswsb(nz,1:nsgs_fields)  ! SGS vertical flux
real sgsadv(nz,1:nsgs_fields)  ! tendency due to vertical advection
real sgslsadv(nz,1:nsgs_fields)  ! tendency due to large-scale vertical advection
real sgsdiff(nz,1:nsgs_fields)  ! tendency due to vertical diffusion

!------------------------------------------------------------------
! internal (optional) definitions:

! make aliases for prognostic variables:

real tke(dimx1_s:dimx2_s, dimy1_s:dimy2_s, nzm)   ! SGS TKE
equivalence (tke(dimx1_s,dimy1_s,1),sgs_field(dimx1_s,dimy1_s,1,1))

! make aliases for diagnostic variables:

real tk  (dimx1_d:dimx2_d, dimy1_d:dimy2_d, nzm) ! SGS eddy viscosity
real tkh (dimx1_d:dimx2_d, dimy1_d:dimy2_d, nzm) ! SGS eddy conductivity
equivalence (tk(dimx1_d,dimy1_d,1), sgs_field_diag(dimx1_d, dimy1_d,1,1))
equivalence (tkh(dimx1_d,dimy1_d,1), sgs_field_diag(dimx1_d, dimy1_d,1,2))

real grdf_x(nzm)! grid factor for eddy diffusion in x
real grdf_y(nzm)! grid factor for eddy diffusion in y
real grdf_z(nzm)! grid factor for eddy diffusion in z

logical:: dosmagor   ! if true, then use Smagorinsky closure

! Local diagnostics:

real tkesbbuoy(nz), tkesbshear(nz),tkesbdiss(nz), tkesbdiff(nz)

CONTAINS

! required microphysics subroutines and function:
!----------------------------------------------------------------------
!!! Read microphysics options from prm (namelist) file

subroutine sgs_setparm()

  use grid, only: case, caseid, masterproc
  implicit none

  integer ierr, ios, ios_missing_namelist, place_holder

  !======================================================================
  NAMELIST /SGS_TKE/ &
       dosmagor ! Diagnostic Smagorinsky closure

  NAMELIST /BNCUIODSBJCB/ place_holder

  dosmagor = .true. ! default

  !----------------------------------
  !  Read namelist for microphysics 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,SGS_TKE,IOSTAT=ios)

  advect_sgs = .not.dosmagor

  if (ios.ne.0) then
     !namelist error checking
     if(ios.ne.ios_missing_namelist) then
        write(*,*) '****** ERROR: bad specification in SGS_TKE namelist'
        call task_abort()
     end if
  end if
  close(55)
   ! write namelist values out to file for documentation
   if(masterproc) then
      open(unit=55,file='./'//trim(case)//'/'//trim(case)//'_'//trim(caseid)//'.nml', form='formatted', position='append')
      write (unit=55,nml=SGS_TKE,IOSTAT=ios)
      write(55,*) ' '
      close(unit=55)
   end if

  ! END UW ADDITION
  !======================================================================

end subroutine sgs_setparm

!----------------------------------------------------------------------
!!! Initialize sgs:

subroutine sgs_init()

  use grid, only: nrestart, dx, dy, dz, adz, masterproc
  use params, only: LES, min_anisotropy
  integer k

  if(nrestart.eq.0) then

     sgs_field = 0.
     sgs_field_diag = 0.

     fluxbsgs = 0.
     fluxtsgs = 0.

  end if

  if(masterproc) then
     if(dosmagor) then
        write(*,*) 'Smagorinsky SGS Closure'
     else
        write(*,*) 'Prognostic TKE 1.5-order SGS Closure'
     end if
  end if

  if (LES) then
    do k=1,nzm
       grdf_x(k) = dx**2/(adz(k)*dz)**2
       grdf_y(k) = dy**2/(adz(k)*dz)**2
       grdf_z(k) = 1.
    end do
  else
    do k=1,nzm
       grdf_x(k) = min(min_anisotropy, dx**2/(adz(k)*dz)**2)
       grdf_y(k) = min(min_anisotropy, dy**2/(adz(k)*dz)**2)
       grdf_z(k) = 1.
    end do
  end if

  sgswle = 0.
  sgswsb = 0.
  sgsadv = 0.
  sgsdiff = 0.
  sgslsadv = 0.

end subroutine sgs_init

!----------------------------------------------------------------------
!!! make some initial noise in sgs:
!
subroutine setperturb_sgs(ptype)

use vars, only: q0, z
integer, intent(in) :: ptype
integer i,j,k

select case (ptype)

  case(0)

     do k=1,nzm
      do j=1,ny
       do i=1,nx
         if(k.le.4.and..not.dosmagor) then
            tke(i,j,k)=0.04*(5-k)
         endif
       end do
      end do
     end do

  case(1)

     do k=1,nzm
      do j=1,ny
       do i=1,nx
         if(q0(k).gt.6.e-3.and..not.dosmagor) then
            tke(i,j,k)=1.
         endif
       end do
      end do
     end do

  case(2)

  case(3)   ! gcss wg1 smoke-cloud case

     do k=1,nzm
      do j=1,ny
       do i=1,nx
         if(q0(k).gt.0.5e-3.and..not.dosmagor) then
            tke(i,j,k)=1.
         endif
       end do
      end do
     end do

  case(4)  ! gcss wg1 arm case

     do k=1,nzm
      do j=1,ny
       do i=1,nx
         if(z(k).le.150..and..not.dosmagor) then
            tke(i,j,k)=0.15*(1.-z(k)/150.)
         endif
       end do
      end do
     end do

  case(5)  ! gcss wg1 BOMEX case

     do k=1,nzm
      do j=1,ny
       do i=1,nx
         if(z(k).le.3000..and..not.dosmagor) then
            tke(i,j,k)=1.-z(k)/3000.
         endif
       end do
      end do
     end do

  case(6)  ! GCSS Lagragngian ASTEX

     do k=1,nzm
      do j=1,ny
       do i=1,nx
         if(q0(k).gt.6.e-3.and..not.dosmagor) then
            tke(i,j,k)=1.
         endif
       end do
      end do
     end do

  case default

end select

end subroutine setperturb_sgs

!----------------------------------------------------------------------
!!! Estimate Courant number limit for SGS
!

subroutine kurant_sgs(cfl)

use grid, only: dt, dx, dy, dz, adz, adzw
implicit none

real, intent(out) :: cfl

integer k
real tkhmax(nz)

do k = 1,nzm
 tkhmax(k) = maxval(tkh(1:nx,1:ny,k))
end do

cfl = 0.
do k=1,nzm
  cfl = max(cfl,        &
     0.5*tkhmax(k)*grdf_z(k)*dt/(dz*adzw(k))**2, &
     0.5*tkhmax(k)*grdf_x(k)*dt/dx**2, &
     YES3D*0.5*tkhmax(k)*grdf_y(k)*dt/dy**2)
end do

end subroutine kurant_sgs

subroutine kurant_sgs_level_by_level(cfl)

use grid, only: dt, dx, dy, dz, adz, adzw
implicit none

real, intent(out) :: cfl(nzm)

integer k
real tkhmax(nz)

do k = 1,nzm
 tkhmax(k) = maxval(tkh(1:nx,1:ny,k))
 cfl(k) = tkhmax(k)*dt*( grdf_x(k)/dx**2 &
                        + YES3D*grdf_y(k)/dy**2 &
                        + grdf_z(k)/(dz*adzw(k))**2 )
end do

end subroutine kurant_sgs_level_by_level

!----------------------------------------------------------------------
!!! compute sgs diffusion of momentum:
!
subroutine sgs_mom()

   call diffuse_mom()

end subroutine sgs_mom

!----------------------------------------------------------------------
!!! compute sgs diffusion of scalars:
!
subroutine sgs_scalars()

  use vars
  use microphysics
  use tracers
  use params, only: dotracers, dosgsthermo
  implicit none

    real dummy(nz)
    real fluxbtmp(nx,ny), fluxttmp(nx,ny) !bloss
    integer k

      if(advect_sgs) then
         call diffuse_scalar(tke,fzero,fzero,dummy,sgswsb, &
                                    dummy,dummy,dummy,.false.)
      end if

!
!    diffusion of microphysics prognostics:
!
!bloss      call micro_flux()

      if (dosgsthermo) then

         call diffuse_scalar(t,fluxbt,fluxtt,tdiff,twsb, &
              t2lediff,t2lediss,twlediff,.true.)

         total_water_evap = total_water_evap - total_water()

         do k = 1,nmicro_fields
            if(   k.eq.index_water_vapor             &! transport water-vapor variable no metter what
                 .or. docloud.and.flag_precip(k).ne.1    & ! transport non-precipitation vars
                 .or. doprecip.and.flag_precip(k).eq.1 ) then
               fluxbtmp(1:nx,1:ny) = fluxbmk(1:nx,1:ny,k)
               fluxttmp(1:nx,1:ny) = fluxtmk(1:nx,1:ny,k)
               call diffuse_scalar(micro_field(:,:,:,k),fluxbtmp,fluxttmp, &
                    mkdiff(:,k),mkwsb(:,k), dummy,dummy,dummy,.false.)
            end if
         end do

         total_water_evap = total_water_evap + total_water()

      end if

 ! diffusion of tracers:

      if(dotracers) then

        call tracers_flux()

        do k = 1,ntracers

          fluxbtmp = fluxbtr(:,:,k)
          fluxttmp = fluxttr(:,:,k)
          call diffuse_scalar(tracer(:,:,:,k),fluxbtmp,fluxttmp, &
               trdiff(:,k),trwsb(:,k), &
               dummy,dummy,dummy,.false.)
!!$          call diffuse_scalar(tracer(:,:,:,k),fluxbtr(:,:,k),fluxttr(:,:,k),trdiff(:,k),trwsb(:,k), &
!!$                           dummy,dummy,dummy,.false.)

        end do

      end if

end subroutine sgs_scalars

!----------------------------------------------------------------------
!!! compute sgs processes (beyond advection):
!
subroutine sgs_proc()

   use grid, only: nstep,dt,icycle
   use params, only: dosmoke

!    SGS TKE equation:

     if(dosgs) call tke_full()

end subroutine sgs_proc

!----------------------------------------------------------------------
!!! Diagnose arrays nessesary for dynamical core and statistics:
!
subroutine sgs_diagnose()
! None

end subroutine sgs_diagnose

!----------------------------------------------------------------------
!!!! Collect microphysics history statistics (vertical profiles)
!
subroutine sgs_statistics()

  use vars
  use hbuffer, only: hbuf_put, hbuf_avg_put
  use params, only : lcond

  real tmp(2), factor_xy
  real tkz(nzm), tkhz(nzm)
  integer i,j,k,n
  character(LEN=6) :: statname  !bloss: for conditional averages

  call t_startf ('sgs_statistics')

  factor_xy = 1./float(nx*ny)

  do k=1,nzm
    tkz(k) = 0.
    tkhz(k) = 0.
    do j=1,ny
    do i=1,nx
      tkz(k)=tkz(k)+tk(i,j,k)
      tkhz(k)=tkhz(k)+tkh(i,j,k)
    end do
    end do
  end do

  call hbuf_avg_put('TKES',tke,dimx1_s,dimx2_s,dimy1_s,dimy2_s,nzm,1.)

  call hbuf_put('TK',tkz,factor_xy)
  call hbuf_put('TKH',tkhz,factor_xy)

!---------------------------------------------------------
! SGS TKE Budget:

         call hbuf_put('ADVTRS',sgswle(:,1),factor_xy)
         call hbuf_put('BUOYAS',tkesbbuoy,factor_xy)
         call hbuf_put('SHEARS',tkesbshear,factor_xy)
         call hbuf_put('DISSIPS',tkesbdiss,factor_xy)

  call t_stopf ('sgs_statistics')

end subroutine sgs_statistics

!----------------------------------------------------------------------
! called when stepout() called

subroutine sgs_print()

 call fminmax_print('tke:',tke,dimx1_s,dimx2_s,dimy1_s,dimy2_s,nzm)
 call fminmax_print('tk:',tk,0,nxp1,1-YES3D,nyp1,nzm)
 call fminmax_print('tkh:',tkh,0,nxp1,1-YES3D,nyp1,nzm)

end subroutine sgs_print

!----------------------------------------------------------------------
!!! Initialize the list of sgs statistics
!
subroutine sgs_hbuf_init(namelist,deflist,unitlist,status,average_type,count,sgscount)
character(*) namelist(*), deflist(*), unitlist(*)
integer status(*),average_type(*),count,sgscount

end subroutine sgs_hbuf_init

end module sgs

max_stars_repo_name: kbren/uwnet

max_issues_repo_name: kbren/uwnet

max_forks_repo_name: kbren/uwnet


text:
$ yes r | ./f                                                         # rock
rock, paper, scissors?  scoring computer choice (r) and your choice (r)
rock, paper, scissors?  scoring computer choice (s) and your choice (r)
rock, paper, scissors?  scoring computer choice (r) and your choice (r)
...
rock, paper, scissors?  scoring computer choice (p) and your choice (r)
 I'm bored out of my skull
 Who's keeping score anyway???
 25.5  4.5
yes: standard output: Broken pipe
yes: write error
$ yes p 2>/dev/null | ./f                                             # paper
rock, paper, scissors?  scoring computer choice (r) and your choice (p)
rock, paper, scissors?  scoring computer choice (s) and your choice (p)
rock, paper, scissors?  scoring computer choice (r) and your choice (p)
rock, paper, scissors?  scoring computer choice (s) and your choice (p)
...
rock, paper, scissors?  scoring computer choice (s) and your choice (p)
 I'm bored out of my skull
 Who's keeping score anyway???
 25.5  4.5
$ yes scissors 2>/dev/null | ./f                                      # scissors
rock, paper, scissors?  scoring computer choice (r) and your choice (s)
rock, paper, scissors?  scoring computer choice (r) and your choice (s)
...
rock, paper, scissors?  scoring computer choice (r) and your choice (s)
 I'm bored out of my skull
 Who's keeping score anyway???
 26.5  3.5
$


text:
module heat_solvers

contains

subroutine solve_heat_explicit(x, u0, t0, tfinal, nsteps, u)
    implicit none
    real(kind=8), intent(in) :: x(0:), u0(0:), t0, tfinal
    integer, intent(in) :: nsteps
    real(kind=8), intent(out) :: u(0:)
    integer :: n, i, nstep
    real(kind=8) :: dx, dt, dtdx2, um1, up1
    real(kind=8), allocatable :: uxx(:)

    n = size(x) - 2
    dx = x(2) - x(1)
    dt = (tfinal - t0) / float(nsteps)

    print 21, n, dx, dt
21  format("Solving heat equation with n = ",i6, " dx = ",e10.3," dt = ",e10.3)

    dtdx2 = dt / dx**2
    print 22, dtdx2
22  format("  dt/dx**2 = ",f7.3)
	if (dtdx2 > 0.5d0) then
		print *, "*** Warning: the explicit method is unstable"
		endif

    allocate(uxx(1:n))

    u = u0

    do nstep=1,nsteps

        do i=1,n
            uxx(i) = (u(i-1) - 2.d0*u(i) + u(i+1)) / dx**2
            enddo

        do i=1,n
            u(i) = u(i) + dt * uxx(i)
            enddo

        enddo

end subroutine solve_heat_explicit

end module heat_solvers

max_stars_repo_name: philipwangdk/HPC

max_issues_repo_name: philipwangdk/HPC

max_forks_repo_name: philipwangdk/HPC


text:
      subroutine compute_db(ctfm,ctfmp,n_A,n_B,n_beta,x,Nx,hx,myzero,phy
     &s_bdy,res)
      implicit none
      integer i
      integer Nx
      real*8 hx
      real*8 myzero
      real*8 ctfm(Nx)
      real*8 ctfmp(Nx)
      real*8 n_A(Nx)
      real*8 n_B(Nx)
      real*8 n_beta(Nx)
      real*8 x(Nx)
      integer phys_bdy(2)
      real*8 res(Nx)
      real*8 qb
      if (phys_bdy(1) .eq. 1) then
      do i=1, 1, 1
      qb = myzero * x(i)
      res(i)=qb
      end do
      endif
      do i=2, Nx-1, 1
      qb = 0.2D1 / ctfm(i) * n_beta(i) + (-0.5000000000000000D0 * n_beta
     #(i - 1) + 0.5000000000000000D0 * n_beta(i + 1) + 0.250000000000000
     #0D0 / n_A(i) * n_beta(i) * (-0.1D1 * n_A(i - 1) + n_A(i + 1)) + 0.
     #5000000000000000D0 / n_B(i) * (-0.1D1 * n_B(i - 1) + n_B(i + 1)) *
     # n_beta(i)) / ctfmp(i) / hx
      res(i)=qb
      end do
      if (phys_bdy(2) .eq. 1) then
      do i=Nx, Nx, 1
      qb = myzero * x(i)
      res(i)=qb
      end do
      endif
      END


text:
C
C $Id: kurv2dp.f,v 1.3 2008-07-27 03:10:11 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
C NOTE: If you make any changes to this software, please remember to
C make the same changes to the corresponding single precision routine.
C
      SUBROUTINE KURV2DP(T,XS,YS,N,X,Y,XP,YP,S,SIGMA)
      DOUBLE PRECISION TN
      DOUBLE PRECISION SIGMAP
      DOUBLE PRECISION DEL1
      DOUBLE PRECISION DEL2
      DOUBLE PRECISION DELS
      DOUBLE PRECISION SUMX
      DOUBLE PRECISION SUMY
      DOUBLE PRECISION D
      DOUBLE PRECISION C1
      DOUBLE PRECISION C2
      DOUBLE PRECISION SIGDEL
      DOUBLE PRECISION SS
      DOUBLE PRECISION DUMMY
      DOUBLE PRECISION S1
      DOUBLE PRECISION S2
c
      INTEGER N
      DOUBLE PRECISION T,XS,YS,X(N),Y(N),XP(N),YP(N),S(N),SIGMA
c
c                                 coded by alan kaylor cline
c                           from fitpack -- january 26, 1987
c                        a curve and surface fitting package
c                      a product of pleasant valley software
c                  8603 altus cove, austin, texas 78759, usa
c
c this subroutine performs the mapping of points in the
c interval (0.,1.) onto a curve in the plane. the subroutine
c kurv1 should be called earlier to determine certain
c necessary parameters. the resulting curve has a parametric
c representation both of whose components are splines under
c tension and functions of the polygonal arclength
c parameter.
c
c on input--
c
c   t contains a real value to be mapped to a point on the
c   curve. the interval (0.,1.) is mapped onto the entire
c   curve, with 0. mapping to (x(1),y(1)) and 1. mapping
c   to (x(n),y(n)). values outside this interval result in
c   extrapolation.
c
c   n contains the number of points which were specified
c   to determine the curve.
c
c   x and y are arrays containing the x- and y-coordinates
c   of the specified points.
c
c   xp and yp are the arrays output from kurv1 containing
c   curvature information.
c
c   s is an array containing the polygonal arclengths of
c   the curve.
c
c and
c
c   sigma contains the tension factor (its sign is ignored).
c
c the parameters n, x, y, xp, yp, s, and sigma should be
c input unaltered from the output of kurv1.
c
c on output--
c
c   xs and ys contain the x- and y-coordinates of the image
c   point on the curve.
c
c none of the input parameters are altered.
c
c this subroutine references package modules intrvldp and
c snhcshdp.
c
c-----------------------------------------------------------
c
c determine interval
c
      TN = S(N)*T
      IM1 = INTRVLDP(TN,S,N)
      I = IM1 + 1
c
c denormalize tension factor
c
      SIGMAP = ABS(SIGMA)*DBLE(N-1)/S(N)
c
c set up and perform interpolation
c
      DEL1 = TN - S(IM1)
      DEL2 = S(I) - TN
      DELS = S(I) - S(IM1)
      SUMX = (X(I)*DEL1+X(IM1)*DEL2)/DELS
      SUMY = (Y(I)*DEL1+Y(IM1)*DEL2)/DELS
      IF (SIGMAP.NE.0.D0) GO TO 1
      D = DEL1*DEL2/ (6.D0*DELS)
      C1 = (DEL1+DELS)*D
      C2 = (DEL2+DELS)*D
      XS = SUMX - XP(I)*C1 - XP(IM1)*C2
      YS = SUMY - YP(I)*C1 - YP(IM1)*C2
      RETURN
    1 SIGDEL = SIGMAP*DELS
      CALL SNHCSHDP(SS,DUMMY,SIGDEL,-1)
      CALL SNHCSHDP(S1,DUMMY,SIGMAP*DEL1,-1)
      CALL SNHCSHDP(S2,DUMMY,SIGMAP*DEL2,-1)
      D = SIGDEL*SIGMAP* (1.D0+SS)
      C1 = DEL1* (S1-SS)/D
      C2 = DEL2* (S2-SS)/D
      XS = SUMX + XP(I)*C1 + XP(IM1)*C2
      YS = SUMY + YP(I)*C1 + YP(IM1)*C2
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


!-------------------------------------------------------------------------------
!> module ATMOSPHERE driver
!!
!! @par Description
!!          Atmosphere module
!!
!! @author Team SCALE
!!
!! @par History
!! @li      2011-11-11 (H.Yashiro) [new]
!! @li      2012-02-15 (H.Yashiro) [add] Microphysics
!! @li      2012-03-23 (H.Yashiro) [mod] Cleaning
!! @li      2013-08-31 (T.Yamaura) [add] Coupler
!!
!<
!-------------------------------------------------------------------------------
module mod_atmos
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  use scale_prof
  use scale_grid_index
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: ATMOS_setup
  public :: ATMOS_step

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
contains

  !-----------------------------------------------------------------------------
  !> Setup atmosphere
  !-----------------------------------------------------------------------------
  subroutine ATMOS_setup
    use scale_time, only: &
       TIME_NOWDATE
    use mod_atmos_vars, only: &
       ATMOS_DYN_TYPE, &
       ATMOS_PHY_SF_TYPE, &
       ATMOS_PHY_TB_TYPE, &
       ATMOS_PHY_MP_TYPE, &
       ATMOS_PHY_RD_TYPE, &
       ATMOS_PHY_AE_TYPE, &
       ATMOS_PHY_CH_TYPE, &
       DENS, &
       MOMZ, &
       MOMX, &
       MOMY, &
       RHOT, &
       QTRC, &
       DENS_tp, &
       MOMZ_tp, &
       MOMX_tp, &
       MOMY_tp, &
       RHOT_tp, &
       QTRC_tp, &
       sw_dyn    => ATMOS_sw_dyn,    &
       sw_phy_sf => ATMOS_sw_phy_sf, &
       sw_phy_tb => ATMOS_sw_phy_tb, &
       sw_phy_mp => ATMOS_sw_phy_mp, &
       sw_phy_rd => ATMOS_sw_phy_rd, &
       sw_phy_ae => ATMOS_sw_phy_ae
    use scale_atmos_solarins, only: &
       ATMOS_SOLARINS_setup
    use scale_atmos_refstate, only: &
       ATMOS_REFSTATE_setup
    use scale_atmos_boundary, only: &
       ATMOS_BOUNDARY_setup
    use mod_atmos_dyn_driver, only: &
       ATMOS_DYN_setup => ATMOS_DYN_driver_setup
    use mod_atmos_phy_sf_driver, only: &
       ATMOS_PHY_SF_setup => ATMOS_PHY_SF_driver_setup
    use mod_atmos_phy_tb_driver, only: &
       ATMOS_PHY_TB_setup => ATMOS_PHY_TB_driver_setup
    use mod_atmos_phy_mp_driver, only: &
       ATMOS_PHY_MP_setup => ATMOS_PHY_MP_driver_setup
    use mod_atmos_phy_rd_driver, only: &
       ATMOS_PHY_RD_setup => ATMOS_PHY_RD_driver_setup
    use mod_atmos_phy_ae_driver, only: &
       ATMOS_PHY_AE_setup => ATMOS_PHY_AE_driver_setup
    implicit none
    !---------------------------------------------------------------------------

    !--- setup solar insolation
    call ATMOS_SOLARINS_setup( TIME_NOWDATE(1) )

    call ATMOS_REFSTATE_setup( DENS, RHOT, QTRC ) ! (in)

    call ATMOS_BOUNDARY_setup( &
         DENS, MOMZ, MOMX, MOMY, RHOT, QTRC ) ! (in)

    ! setup each components
    if ( sw_dyn    ) call ATMOS_DYN_setup   ( ATMOS_DYN_TYPE    )

    if ( sw_phy_sf ) call ATMOS_PHY_SF_setup( ATMOS_PHY_SF_TYPE )

    if ( sw_phy_tb ) call ATMOS_PHY_TB_setup( ATMOS_PHY_TB_TYPE )

    if ( sw_phy_mp ) call ATMOS_PHY_MP_setup( ATMOS_PHY_MP_TYPE )

    if ( sw_phy_rd .or. sw_phy_ae ) call ATMOS_PHY_AE_setup( ATMOS_PHY_AE_TYPE )

    if ( sw_phy_rd ) call ATMOS_PHY_RD_setup( ATMOS_PHY_RD_TYPE )

    !########## initialize tendencies ##########
    DENS_tp(:,:,:)   = 0.0_RP
    MOMZ_tp(:,:,:)   = 0.0_RP
    MOMX_tp(:,:,:)   = 0.0_RP
    MOMY_tp(:,:,:)   = 0.0_RP
    RHOT_tp(:,:,:)   = 0.0_RP
    QTRC_tp(:,:,:,:) = 0.0_RP

    return
  end subroutine ATMOS_setup

  !-----------------------------------------------------------------------------
  !> advance atmospheric state
  !-----------------------------------------------------------------------------
  subroutine ATMOS_step
    use scale_time, only: &
       do_dyn    => TIME_DOATMOS_DYN,    &
       do_phy_sf => TIME_DOATMOS_PHY_SF, &
       do_phy_tb => TIME_DOATMOS_PHY_TB, &
       do_phy_mp => TIME_DOATMOS_PHY_MP, &
       do_phy_rd => TIME_DOATMOS_PHY_RD, &
       do_phy_ae => TIME_DOATMOS_PHY_AE
    use mod_atmos_vars, only: &
       DENS, &
       RHOT, &
       QTRC, &
       DENS_tp, &
       MOMZ_tp, &
       MOMX_tp, &
       MOMY_tp, &
       RHOT_tp, &
       QTRC_tp, &
       sw_dyn    => ATMOS_sw_dyn,    &
       sw_phy_sf => ATMOS_sw_phy_sf, &
       sw_phy_tb => ATMOS_sw_phy_tb, &
       sw_phy_mp => ATMOS_sw_phy_mp, &
       sw_phy_rd => ATMOS_sw_phy_rd, &
       sw_phy_ae => ATMOS_sw_phy_ae, &
       ATMOS_vars_history
    use mod_atmos_dyn_driver, only: &
       ATMOS_DYN => ATMOS_DYN_driver
    use mod_atmos_phy_sf_driver, only: &
       ATMOS_PHY_SF => ATMOS_PHY_SF_driver, &
       ATMOS_PHY_SF_driver_first,  &
       ATMOS_PHY_SF_driver_final
    use mod_atmos_phy_tb_driver, only: &
       ATMOS_PHY_TB => ATMOS_PHY_TB_driver
    use mod_atmos_phy_mp_driver, only: &
       ATMOS_PHY_MP => ATMOS_PHY_MP_driver
    use mod_atmos_phy_rd_driver, only: &
       ATMOS_PHY_RD => ATMOS_PHY_RD_driver
    use mod_atmos_phy_ae_driver, only: &
       ATMOS_PHY_AE => ATMOS_PHY_AE_driver
    use scale_atmos_refstate, only: &
       ATMOS_REFSTATE_update, &
       ATMOS_REFSTATE_UPDATE_FLAG
    use mod_cpl_vars, only: &
       sw_AtmLnd => CPL_sw_AtmLnd, &
       sw_AtmOcn => CPL_sw_AtmOcn
    implicit none
    !---------------------------------------------------------------------------

    !########## Reference State ###########
    if ( ATMOS_REFSTATE_UPDATE_FLAG ) then
       call ATMOS_REFSTATE_update( DENS, RHOT, QTRC ) ! (in)
    endif

    !########## Surface First ##########
    if ( sw_phy_sf ) then
       call PROF_rapstart('ATM SurfaceFlux')
       if ( sw_AtmLnd .or. sw_AtmOcn ) then
          call ATMOS_PHY_SF_driver_first
       else
          call ATMOS_PHY_SF( do_phy_sf, .true. )
       endif
       call PROF_rapend  ('ATM SurfaceFlux')
    endif

    !########## Turbulence ##########
    if ( sw_phy_tb ) then
       call PROF_rapstart('ATM Turbulence')
       call ATMOS_PHY_TB( do_phy_tb, .true. )
       call PROF_rapend  ('ATM Turbulence')
    endif

    !########## Microphysics ##########
    if ( sw_phy_mp ) then
       call PROF_rapstart('ATM Microphysics')
       call ATMOS_PHY_MP( do_phy_mp, .true. )
       call PROF_rapend  ('ATM Microphysics')
    endif

    !########## Aerosol ##########
    if ( sw_phy_ae ) then
       call PROF_rapstart('ATM Aerosol')
       call ATMOS_PHY_AE( do_phy_ae, .true. )
       call PROF_rapend  ('ATM Aerosol')
    endif

    !########## Radiation ##########
    if ( sw_phy_rd ) then
       call PROF_rapstart('ATM Radiation')
       call ATMOS_PHY_RD( do_phy_rd, .true. )
       call PROF_rapend  ('ATM Radiation')
    endif

    !########## Dynamics ##########
    if ( sw_dyn ) then
       call PROF_rapstart('ATM Dynamics')
       call ATMOS_DYN( do_dyn )
       call PROF_rapend  ('ATM Dynamics')
    endif

    !########## Surface Final ##########
    if ( sw_phy_sf ) then
       call PROF_rapstart('ATM SurfaceFlux')
       if ( sw_AtmLnd .or. sw_AtmOcn ) then
          call ATMOS_PHY_SF_driver_final
       endif
       call PROF_rapend  ('ATM SurfaceFlux')
    endif

    !########## History & Monitor ##########
    call PROF_rapstart('ATM History Vars')
    call ATMOS_vars_history
    call PROF_rapend  ('ATM History Vars')

    !########## reset tendencies ##########
    DENS_tp(:,:,:)   = 0.0_RP
    MOMZ_tp(:,:,:)   = 0.0_RP
    MOMX_tp(:,:,:)   = 0.0_RP
    MOMY_tp(:,:,:)   = 0.0_RP
    RHOT_tp(:,:,:)   = 0.0_RP
    QTRC_tp(:,:,:,:) = 0.0_RP

    return
  end subroutine ATMOS_step

end module mod_atmos