text:
! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! H0 X
! H0 X   libAtoms+QUIP: atomistic simulation library
! H0 X
! H0 X   Portions of this code were written by
! H0 X     Albert Bartok-Partay, Silvia Cereda, Gabor Csanyi, James Kermode,
! H0 X     Ivan Solt, Wojciech Szlachta, Csilla Varnai, Steven Winfield.
! H0 X
! H0 X   Copyright 2006-2010.
! H0 X
! H0 X   These portions of the source code are released under the GNU General
! H0 X   Public License, version 2, http://www.gnu.org/copyleft/gpl.html
! H0 X
! H0 X   If you would like to license the source code under different terms,
! H0 X   please contact Gabor Csanyi, gabor@csanyi.net
! H0 X
! H0 X   Portions of this code were written by Noam Bernstein as part of
! H0 X   his employment for the U.S. Government, and are not subject
! H0 X   to copyright in the USA.
! H0 X
! H0 X
! H0 X   When using this software, please cite the following reference:
! H0 X
! H0 X   http://www.libatoms.org
! H0 X
! H0 X  Additional contributions by
! H0 X    Alessio Comisso, Chiara Gattinoni, and Gianpietro Moras
! H0 X
! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

! Calculates density in a grid of boxes (possibly smoothed) or a radial mesh
! with error bars
! outputs number density (#/vol)

program density_new
use libatoms_module
implicit none

  type(Dictionary) :: cli_params
  character(len=STRING_LENGTH) :: infilename
  logical :: infile_is_list
  character(len=STRING_LENGTH) :: outfilename
  type(Inoutput) :: outfile
  character(len=STRING_LENGTH) :: commandfilename
  type(Inoutput) :: commandfile
  character(len=1024) :: args_str

  character(len=STRING_LENGTH) :: mask_str
  integer :: decimation
  real(dp) :: min_time, max_time
  logical :: gaussian_smoothing
  real(dp) :: gaussian_sigma
  logical :: radial_histo, random_samples
  real(dp) :: radial_center(3)
  integer :: n_samples(2)
  logical :: sort_Time, no_Time_dups
  real(dp) :: min_p(3), bin_width(3)
  integer :: n_bins(3)
  logical :: quiet
  logical :: autocorrelation, mean
  integer :: autocorrelation_max_lag
  real(dp) :: mean_decorrelation_time
  logical :: have_params
  integer :: status, arg_line_no

  logical :: no_compute_index
  integer :: n_histos
  type(Atoms_ll) :: structure_ll
  real(dp), allocatable :: histograms(:,:,:,:), histo_grid(:,:,:,:), histo_mean(:,:,:), histo_var(:,:,:)
  integer :: i_lag, i1, i2, i3
  real(dp), allocatable :: autocorr(:,:,:,:)

  call system_initialise(PRINT_NORMAL)

  call initialise(cli_params)
  call register_cli_params(cli_params,.true., infilename, infile_is_list, outfilename, commandfilename, &
    mask_str, min_p, bin_width, n_bins, decimation, min_time, max_time, gaussian_smoothing, gaussian_sigma, &
    radial_histo, radial_center, random_samples, n_samples, &
    sort_Time, no_Time_dups, mean, mean_decorrelation_time, autocorrelation, autocorrelation_max_lag, quiet)
  if (.not. param_read_args(cli_params)) then
      call print_usage()
    call system_abort('could not parse argument line')
  end if

  call print("Reading configurations")
  call read_xyz(structure_ll, infilename, infile_is_list, decimation, min_time, max_time, sort_Time, no_Time_dups, quiet, no_compute_index)

  if (len_trim(commandfilename) == 0) then
    args_str = write_string(cli_params)
    call finalise(cli_params)
  else
    call finalise(cli_params)
    call print("reading commands from file '"//trim(commandfilename)//"'")
    call initialise(commandfile, trim(commandfilename), INPUT)
    arg_line_no = 1
    args_str = read_line(commandfile)
    call initialise(cli_params)
    call print("got arguments line '"//trim(args_str)//"'")
    call register_cli_params(cli_params,.false., infilename, infile_is_list, outfilename, commandfilename, &
      mask_str, min_p, bin_width, n_bins, decimation, min_time, max_time, gaussian_smoothing, gaussian_sigma, &
      radial_histo, radial_center, random_samples, n_samples, &
      sort_Time, no_Time_dups, mean, mean_decorrelation_time, autocorrelation, autocorrelation_max_lag, quiet)
    if (.not. param_read_line(cli_params, trim(args_str), ignore_unknown = .true.)) then
      call print_usage()
      call system_abort("could not parse argument line "//arg_line_no//" from file '"//trim(commandfilename//"'"))
    end if
    call finalise(cli_params)
  endif

  have_params = .true.
  do while (have_params)
    no_compute_index=.false.
    if ((.not. sort_Time) .and. (.not. no_Time_dups)) no_compute_index=.true.

    call print("infile " // trim(infilename) // " infile_is_list " // infile_is_list)
    call print("outfilename " // trim(outfilename))
    call print("decimation " // decimation // " min_time " // min_time // " max_time " // max_time)
    call print("AtomMask " // trim(mask_str))
    call print("radial_histo " // radial_histo)
    if (radial_histo) then
    call print("radial_center " // radial_center)
      call print("bin_width " // bin_width(1) // " n_bins " // n_bins(1))
      n_bins(2:3) = 1
      call print("random_samples " // random_samples)
      if (random_samples) then
        call print("n_samples " // n_samples(1))
      else
        call print("n_theta " // n_samples(1) // " n_phi " // n_samples(2))
      endif
    else
      call print("min_p " // min_p)
      call print("bin_width " // bin_width // " n_bins " // n_bins)
    endif
    call print("gaussian " // gaussian_smoothing // " sigma " // gaussian_sigma)
    call print("sort_Time " // sort_Time // " no_Time_dups " // no_Time_dups)
    call print("mean " // mean // " mean_decorrelation_time " // mean_decorrelation_time)
    call print("autocorrelation " // autocorrelation // " autocorrelation_max_lag " // autocorrelation_max_lag)

    call print("Calculating densities")
    call calc_histos(histograms, n_histos, histo_grid, min_p, bin_width, n_bins, structure_ll, &
      mean_decorrelation_time, gaussian_smoothing, gaussian_sigma, &
      radial_histo, radial_center, random_samples, n_samples, mask_str, quiet)

    call initialise(outfile, outfilename, OUTPUT)

    if (autocorrelation) then
      call print("Calculating autocorrelations")
      allocate(autocorr(n_bins(1),n_bins(2),n_bins(3),autocorrelation_max_lag+1))
      autocorr = autocorrelation_3array(histograms(:,:,:,1:n_histos), max_lag=autocorrelation_max_lag)

      call print("Printing autocorrelations")
      call print("# Autocorrelation", file=outfile)
      do i_lag=0, autocorrelation_max_lag
        call print(i_lag // " " // reshape(autocorr(1:n_bins(1),1:n_bins(2),1:n_bins(3),i_lag+1), &
                                              (/ n_bins(1)*n_bins(2)*n_bins(3) /)),  file=outfile)
      end do
    endif

    if (mean) then
      if (autocorrelation) then
        call print("", file=outfile)
        call print("", file=outfile)
      endif
      allocate(histo_mean(n_bins(1), n_bins(2), n_bins(3)))
      allocate(histo_var(n_bins(1), n_bins(2), n_bins(3)))
      call calc_mean_var_3array(histograms(:,:,:,1:n_histos), histo_mean, histo_var)

      if (radial_histo) then
        call print("# Density (#/vol) ", file=outfile)
        call print("# r mean var n_samples", file=outfile)
        do i1=1, n_bins(1)
          call print(histo_grid(1,i1,1,1) // " " // histo_mean(i1, 1, 1) // " " // &
          histo_var(i1, 1, 1) // " " // n_histos, file=outfile)
        end do
      else
        call print("# Density (#/vol)", file=outfile)
        call print("# x y z   mean   var  n_samples", file=outfile)
        do i1=1,n_bins(1)
        do i2=1,n_bins(2)
        do i3=1,n_bins(3)
          call print( histo_grid(:,i1,i2,i3) // &
                      " " // histo_mean(i1, i2, i3) // &
                      " " // histo_var(i1, i2, i3) // " " // n_histos, file=outfile)
        end do
        call print('', file=outfile)
        end do
        call print('', file=outfile)
        end do
      endif
    end if

    call finalise(outfile)

    if (allocated(histograms)) deallocate(histograms)
    if (allocated(histo_grid)) deallocate(histo_grid)
    if (allocated(autocorr)) deallocate(autocorr)
    if (allocated(histo_mean)) deallocate(histo_mean)
    if (allocated(histo_var)) deallocate(histo_var)

    if (len_trim(commandfilename) == 0) then
      have_params = .false.
    else
      arg_line_no = arg_line_no + 1
      args_str = read_line(commandfile, status=status)
      if (status == 0) then
        call print("got arguments line '"//trim(args_str)//"'")
        call initialise(cli_params)
        call register_cli_params(cli_params,.false., infilename, infile_is_list, outfilename, commandfilename, &
          mask_str, min_p, bin_width, n_bins, decimation, min_time, max_time, gaussian_smoothing, gaussian_sigma, &
          radial_histo, radial_center, random_samples, n_samples, &
          sort_Time, no_Time_dups, mean, mean_decorrelation_time, autocorrelation, autocorrelation_max_lag, quiet)
        if (.not. param_read_line(cli_params, trim(args_str), ignore_unknown = .true.)) then
          call print_usage()
          call system_abort("could not parse argument line "//arg_line_no//" from file '"//trim(commandfilename//"'"))
        end if
        call finalise(cli_params)
      else
        have_params = .false.
      endif
    endif

  end do

  call system_finalise()

contains

  subroutine print_usage()

    character(len=1024) my_exec_name

    if (EXEC_NAME == "<UNKNOWN>") then
      my_exec_name=""
    else
      my_exec_name=EXEC_NAME
    endif

    call print("Usage: " // trim(my_exec_name)//" infile=filename [infile_is_list=logical(F)]", PRINT_ALWAYS)
    call print("       outfile=filename [commandfile=filename()] [AtomMask=species()] [min_p={x y z}]", PRINT_ALWAYS)
    call print("       bin_width={x y z}(y,z ignored for radial_histo) n_bins={nx ny nz}(y,z ignored for radial_histo)", PRINT_ALWAYS)
    call print("       [mean=logical(T)] [mean_decorrelation_time=t(0.0)]", PRINT_ALWAYS)
    call print("       [autocorrelation=logical(F)] [autocorrelation_max_lag=N(10000)]", PRINT_ALWAYS)
    call print("       [decimation=n(1)] [min_time=t(-1.0)] [max_time=t(-1.0)]", PRINT_ALWAYS)
    call print("       [gaussian=logical(F)(always true for radial_histo)] [sigma=s(1.0)] [radial_histo=logical(F)]", PRINT_ALWAYS)
    call print("       [radial_center={x y z}] [random_samples=logical(F)(radial_histo only)] [n_samples={2 4})(component 2 ignored for random)]", PRINT_ALWAYS)
    call print("       [sort_Time(F)] [no_Time_dups(F)]", PRINT_ALWAYS)
    call print("       [quiet=logical(F)]", PRINT_ALWAYS)
  end subroutine print_usage

  subroutine register_cli_params(cli_params, initial, infilename, infile_is_list, outfilename, commandfilename, &
    mask_str, min_p, bin_width, n_bins, decimation, min_time, max_time, gaussian_smoothing, gaussian_sigma, &
    radial_histo, radial_center, random_samples, n_samples, &
    sort_Time, no_Time_dups, mean, mean_decorrelation_time, autocorrelation, autocorrelation_max_lag, quiet)
    type(Dictionary), intent(inout) :: cli_params
    logical, intent(in) :: initial
    character(len=*), intent(inout) :: infilename
    logical, intent(inout) :: infile_is_list
    character(len=*), intent(inout) :: outfilename
    character(len=*), intent(inout) :: commandfilename
    character(len=*), intent(inout) :: mask_str
    real(dp), intent(inout) :: min_p(3), bin_width(3)
    integer, intent(inout) :: n_bins(3), decimation
    real(dp), intent(inout) :: min_time, max_time
    logical, intent(inout) :: gaussian_smoothing
    real(dp), intent(inout) :: gaussian_sigma
    logical, intent(inout) :: radial_histo
    real(dp), intent(inout) :: radial_center(3)
    logical, intent(inout) :: random_samples
    integer, intent(inout) :: n_samples(2)
    logical, intent(inout) :: sort_Time, no_Time_dups
    logical, intent(inout) :: mean
    real(dp), intent(inout) :: mean_decorrelation_time
    logical, intent(inout) :: autocorrelation
    integer, intent(inout) :: autocorrelation_max_lag
    logical, intent(inout) :: quiet
    character(len=STRING_LENGTH) :: t_field

    if (initial) then
      call param_register(cli_params, 'infile', param_mandatory, infilename, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'infile_is_list', 'F', infile_is_list, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'outfile', 'stdout', outfilename, help_string="No help yet.  This source file was $LastChangedBy$")
      mask_str = ""
      commandfilename = ""
      call param_register(cli_params, 'commandfile', "", commandfilename, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'AtomMask', "", mask_str, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'min_p', "0.0 0.0 0.0", min_p, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'bin_width', PARAM_MANDATORY, bin_width, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'n_bins', PARAM_MANDATORY, n_bins, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'decimation', '1', decimation, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'min_time', '-1.0', min_time, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'max_time', '-1.0', max_time, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'gaussian', 'F', gaussian_smoothing, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'sigma', '1.0', gaussian_sigma, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'radial_histo', 'F', radial_histo, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'radial_center', '0.0 0.0 0.0', radial_center, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'random_samples', 'F', random_samples, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'n_samples', '2 4', n_samples, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'sort_Time', 'F', sort_Time, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'no_Time_dups', 'F', no_Time_dups, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'mean', 'T', mean, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'mean_decorrelation_time', '0.0', mean_decorrelation_time, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'autocorrelation', 'F', autocorrelation, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'autocorrelation_max_lag', '10000', autocorrelation_max_lag, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'quiet', 'F', quiet, help_string="No help yet.  This source file was $LastChangedBy$")
    else
      t_field=infilename
      call param_register(cli_params, 'infile', trim(t_field), infilename, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'infile_is_list', ""//infile_is_list, infile_is_list, help_string="No help yet.  This source file was $LastChangedBy$")
      t_field=outfilename
      call param_register(cli_params, 'outfile', trim(t_field), outfilename, help_string="No help yet.  This source file was $LastChangedBy$")
      t_field=commandfilename
      call param_register(cli_params, 'commandfile', trim(t_field), commandfilename, help_string="No help yet.  This source file was $LastChangedBy$")
      t_field=mask_str
      call param_register(cli_params, 'AtomMask', trim(t_field), mask_str, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'min_p', ""//min_p, min_p, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'bin_width', ""//bin_width, bin_width, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'n_bins', ""//n_bins, n_bins, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'min_time', ""//min_time, min_time, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'max_time', ""//max_time, max_time, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'gaussian', ""//gaussian_smoothing, gaussian_smoothing, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'sigma', ""//gaussian_sigma, gaussian_sigma, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'radial_histo', ""//radial_histo, radial_histo, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'radial_center', ""//radial_center, radial_center, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'random_samples', ""//random_samples, random_samples, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'n_samples', '2 4', n_samples, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'mean', ""//mean, mean, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'mean_decorrelation_time', ""//mean_decorrelation_time, mean_decorrelation_time, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'autocorrelation', ""//autocorrelation, autocorrelation, help_string="No help yet.  This source file was $LastChangedBy$")
      call param_register(cli_params, 'autocorrelation_max_lag', ""//autocorrelation_max_lag, autocorrelation_max_lag, help_string="No help yet.  This source file was $LastChangedBy$")
    endif
  end subroutine register_cli_params

  function autocorrelation_scalar(v, max_lag) result(autocorrelation)
    real(dp), intent(in) :: v(:)
    integer, intent(in) :: max_lag
    real(dp) :: autocorrelation(max_lag)

    integer :: lag, i
    real(dp) :: t_sum, mean
    integer :: nt

    nt = size(v,1)

    mean = sum(v)/real(nt,dp)
    do lag=1, max_lag
      t_sum = 0.0_dp
      do i=1, nt-lag
        t_sum = t_sum + (v(i)-mean)*(v(i+lag)-mean)
      end do
      autocorrelation(lag) = t_sum/real(nt-lag,dp)
    end do

  end function autocorrelation_scalar

  function autocorrelation_3array(v, max_lag) result(autocorrelation)
    real(dp), intent(in) :: v(:,:,:,:)
    integer, intent(in) :: max_lag
    real(dp) :: autocorrelation(size(v,1),size(v,2),size(v,3),max_lag+1)

    integer :: lag, i
    real(dp), allocatable :: t_sum(:,:,:), mean(:,:,:)
    integer :: n1, n2, n3, nt

    n1 = size(v,1)
    n2 = size(v,2)
    n3 = size(v,3)
    nt = size(v,4)
    allocate(t_sum(n1,n2,n3), mean(n1,n2,n3))

    mean = sum(v,4)/real(nt,dp)
    do lag=0, max_lag
      t_sum = 0.0_dp
      do i=1, nt-lag
        t_sum(1:n1,1:n2,1:n3) = t_sum(1:n1,1:n2,1:n3) + (v(1:n1,1:n2,1:n3,i)-mean)*(v(1:n1,1:n2,1:n3,i+lag)-mean)
      end do
      autocorrelation(1:n1,1:n2,1:n3,lag+1) = t_sum(1:n1,1:n2,1:n3)/real(nt-lag,dp)
    end do

    deallocate(t_sum, mean)

  end function autocorrelation_3array

  subroutine calc_histos(histograms, n_histos, histo_grid, min_p, bin_width, n_bins, &
                         structure_ll, interval, gaussian, gaussian_sigma, &
                         radial_histo, radial_center, random_samples, n_samples, mask_str, quiet)
    real(dp), intent(inout), allocatable :: histograms(:,:,:,:)
    integer, intent(out) :: n_histos
    real(dp), intent(inout), allocatable :: histo_grid(:,:,:,:)
    real(dp), intent(in) :: min_p(3), bin_width(3)
    integer, intent(in) :: n_bins(3)
    type(atoms_ll), intent(in) :: structure_ll
    real(dp), intent(in) :: interval
    logical, intent(in) :: gaussian
    real(dp), intent(in) :: gaussian_sigma
    logical, intent(in) :: radial_histo
    real(dp), intent(in) :: radial_center(3)
    logical :: random_samples
    integer :: n_samples(2)
    character(len=*), optional, intent(in) :: mask_str
    logical, optional :: quiet

    real(dp) :: last_time, cur_time
    type(atoms_ll_entry), pointer :: entry
    logical :: do_this_histo, first_histo
    real(dp), allocatable :: theta(:), phi(:), w(:)
    logical :: my_quiet

    my_quiet = optional_default(.false., quiet)
    n_histos = 0
    allocate(histograms(n_bins(1),n_bins(2),n_bins(3),10))
    if (radial_histo) then
      allocate(histo_grid(1,n_bins(1),n_bins(2),n_bins(3)))
    else
      allocate(histo_grid(3,n_bins(1),n_bins(2),n_bins(3)))
    endif
    histograms = 0.0_dp

    if (radial_histo) then
      if (random_samples) then
        call calc_angular_samples_random(n_samples(1), theta, phi, w)
      else
        call calc_angular_samples_grid(n_samples(1), n_samples(2), theta, phi, w)
      endif
    endif

    first_histo = .true.
    entry => structure_ll%first
    cur_time = -1.0_dp
    last_time = -1.0e38_dp
    do while (associated(entry))
      do_this_histo = .false.
      if (interval > 0.0_dp) then
        if (cur_time < last_time) &
          call system_abort("calc_histos called with interval="//interval//" > 0, but unsorted sequence: cur_time="// &
                            cur_time//" < last_time="//last_time)
        if (.not. get_value(entry%at%params, "Time", cur_time)) &
          call system_abort("calc_histos called with interval="//interval//" > 0, but no Time value in entry")
        if ((cur_time .feq. last_time) .or. (cur_time >= last_time+interval)) then
          do_this_histo = .true.
          if (cur_time >= last_time+interval) last_time = cur_time
        endif
      else
        do_this_histo = .true.
      endif
      if (do_this_histo) then
        ! if (get_value(entry%at%params, "Time", cur_time)) call print("doing histo for config with time " // cur_time)
        n_histos = n_histos + 1
        if (.not. my_quiet) then
          if (mod(n_histos,1000) == 1) write (mainlog%unit,'(I4,a)') int(n_histos/1000)," "
          if (mod(n_histos,10) == 0) write (mainlog%unit,'(I1,$)') mod(n_histos/10,10)
          if (mod(n_histos,1000) == 0) write (mainlog%unit,'(a)') " "
        endif
        call reallocate_histos(histograms, n_histos, n_bins)
        if (radial_histo) then
          if (first_histo) then
            call sample_radial_mesh_Gaussians(histograms(:,1,1,n_histos), entry%at, radial_center, bin_width(1), n_bins(1), gaussian_sigma, &
              theta, phi, w, mask_str, histo_grid(1,:,1,1))
          else
            call sample_radial_mesh_Gaussians(histograms(:,1,1,n_histos), entry%at, radial_center, bin_width(1), n_bins(1), gaussian_sigma, &
              theta, phi, w, mask_str)
          endif
        else ! rectilinear
          if (gaussian) then
            if (first_histo) then
              call sample_rectilinear_mesh_Gaussians(histograms(:,:,:,n_histos), entry%at, min_p, bin_width, n_bins, gaussian_sigma, mask_str, histo_grid(:,:,:,:))
            else
              call sample_rectilinear_mesh_Gaussians(histograms(:,:,:,n_histos), entry%at, min_p, bin_width, n_bins, gaussian_sigma, mask_str)
            endif
          else
            if (first_histo) then
              call calc_histogram(histograms(:,:,:,n_histos), entry%at, min_p, bin_width, n_bins, mask_str, histo_grid(:,:,:,:))
            else
              call calc_histogram(histograms(:,:,:,n_histos), entry%at, min_p, bin_width, n_bins, mask_str)
            endif
          endif
        endif
      endif
      entry => entry%next
    end do
  end subroutine calc_histos

  subroutine calc_angular_samples_grid(n_t, n_p, theta, phi, w)
    integer, intent(in) :: n_t, n_p
    real(dp), allocatable, intent(inout) :: theta(:), phi(:), w(:)

    integer :: t_i, p_i, ii

    if (allocated(theta)) deallocate(theta)
    if (allocated(phi)) deallocate(phi)
    if (allocated(w)) deallocate(w)
    allocate(theta(n_t*n_p))
    allocate(phi(n_t*n_p))
    allocate(w(n_t*n_p))

    ii = 1
    do t_i=1, n_t
    do p_i=1, n_p
      phi(ii) = (p_i-1)*2.0_dp*PI/real(n_p,dp)
      if (mod(n_t,2) == 0) then
        theta(ii) = (t_i-1.5_dp-floor((n_t-1)/2.0_dp))*PI/real(n_t,dp)
      else
        theta(ii) = (t_i-1-floor((n_t-1)/2.0_dp))*PI/real(n_t,dp)
      endif
      ! not oversample top and bottom of spehre
      w(ii) = cos(theta(ii))
      ii = ii + 1
    end do
    end do
    w = w / ( (gaussian_sigma*sqrt(2.0_dp*PI))**3 * sum(w) )
  end subroutine calc_angular_samples_grid

  subroutine calc_angular_samples_random(n, theta, phi, w)
    integer, intent(in) :: n
    real(dp), allocatable, intent(inout) :: theta(:), phi(:), w(:)

    integer :: ii
    real(dp) :: p(3), p_norm

    if (allocated(theta)) deallocate(theta)
    if (allocated(phi)) deallocate(phi)
    if (allocated(w)) deallocate(w)
    allocate(theta(n))
    allocate(phi(n))
    allocate(w(n))

    do ii = 1, n
      p_norm = 2.0_dp
      do while (p_norm > 1.0_dp)
        p = 0.0_dp
        call randomise(p, 2.0_dp)
        p_norm = norm(p)
      end do
      phi(ii) = atan2(p(2),p(1))
      theta(ii) = acos(p(3)/p_norm)
      w(ii) = 1.0_dp
    end do
    w = w / ( (gaussian_sigma*sqrt(2.0_dp*PI))**3 * sum(w) )

  end subroutine calc_angular_samples_random

  subroutine sample_radial_mesh_Gaussians(histogram, at, radial_center, bin_width, n_bins, gaussian_sigma, theta, phi, w, mask_str, histo_grid, accumulate)
    real(dp), intent(inout) :: histogram(:)
    type(Atoms), intent(in) :: at
    real(dp), intent(in) :: radial_center(3)
    real(dp), intent(in) :: bin_width
    integer, intent(in) :: n_bins
    real(dp), intent(in) :: gaussian_sigma, theta(:), phi(:), w(:)
    character(len=*), optional, intent(in) :: mask_str
    real(dp), intent(out), optional :: histo_grid(:)
    logical, optional, intent(in) :: accumulate

    logical :: my_accumulate
    real(dp) :: bin_r, p(3), dist, r
    logical, allocatable :: mask_a(:)
    integer at_i, sample_i, bin_i

    my_accumulate = optional_default(.false., accumulate)
    if (.not. my_accumulate) histogram = 0.0_dp

    allocate(mask_a(at%N))
    call is_in_mask(mask_a, at, mask_str)

! ratio of 20/4=5 is bad
! ratio of 20/3=6.66 is bad
! ratio of 20/2.5=8 is borderline (2e-4)
! ratio of 20/2.22=9 is fine  (error 2e-5)
! ratio of 20/2=10 is fine
    if ( (norm(at%lattice(:,1)) < 9.0_dp*gaussian_sigma) .or. &
         (norm(at%lattice(:,2)) < 9.0_dp*gaussian_sigma) .or. &
         (norm(at%lattice(:,3)) < 9.0_dp*gaussian_sigma) ) &
      call print("WARNING: at%lattice may be too small for sigma, errors (noticeably too low a density) may result", PRINT_ALWAYS)

    if (present(histo_grid)) then
      do bin_i=1, n_bins
        bin_r = (bin_i-1)*bin_width
        histo_grid(bin_i) = bin_r
      end do
    endif

    do at_i=1, at%N
      if (.not. mask_a(at_i)) cycle
      r = norm(at%pos(:,at_i))
      do bin_i=1, n_bins
        bin_r = (bin_i-1)*bin_width
        do sample_i=1,size(theta)
          p(1) = radial_center(1)+bin_r*cos(phi(sample_i))*cos(theta(sample_i))
          p(2) = radial_center(2)+bin_r*sin(phi(sample_i))*cos(theta(sample_i))
          p(3) = radial_center(3)+bin_r*sin(theta(sample_i))
          dist = distance_min_image(at,p,at%pos(:,at_i))
!Include all the atoms, slow but minimises error
!          if (dist > 4.0_dp*gaussian_sigma) cycle
          histogram(bin_i) = histogram(bin_i) + exp(-0.5_dp*(dist/(gaussian_sigma))**2)*w(sample_i)
        end do ! sample_i
      end do ! bin_i
    end do ! at_i

    deallocate(mask_a)
  end subroutine sample_radial_mesh_Gaussians

  subroutine sample_rectilinear_mesh_Gaussians(histogram, at, min_p, sample_dist, n_bins, gaussian_sigma, mask_str, histo_grid, accumulate)
    real(dp), intent(inout) :: histogram(:,:,:)
    type(Atoms), intent(in) :: at
    real(dp), intent(in) :: min_p(3), sample_dist(3)
    integer, intent(in) :: n_bins(3)
    real(dp), intent(in) :: gaussian_sigma
    character(len=*), optional, intent(in) :: mask_str
    real(dp), intent(out), optional :: histo_grid(:,:,:,:)
    logical, optional, intent(in) :: accumulate

    logical :: my_accumulate
    integer :: at_i, i1, i2, i3
    real(dp) :: p(3), w, dist
    logical, allocatable :: mask_a(:)

    my_accumulate = optional_default(.false., accumulate)
    if (.not. my_accumulate) histogram = 0.0_dp

    allocate(mask_a(at%N))
    call is_in_mask(mask_a, at, mask_str)

! ratio of 20/4=5 is bad
! ratio of 20/3=6.66 is bad
! ratio of 20/2.5=8 is borderline (2e-4)
! ratio of 20/2.22=9 is fine  (error 2e-5)
! ratio of 20/2=10 is fine
    if ( (norm(at%lattice(:,1)) < 9.0_dp*gaussian_sigma) .or. &
         (norm(at%lattice(:,2)) < 9.0_dp*gaussian_sigma) .or. &
         (norm(at%lattice(:,3)) < 9.0_dp*gaussian_sigma) ) &
      call print("WARNING: at%lattice may be too small for sigma, errors (noticeably too low a density) may result", PRINT_ALWAYS)

    w = 1.0_dp / (gaussian_sigma*sqrt(2.0_dp*PI))**3

    if (present(histo_grid)) then
      do i1=1, n_bins(1)
        p(1) = min_p(1) + (i1-1)*sample_dist(1)
        do i2=1, n_bins(2)
          p(2) = min_p(2) + (i2-1)*sample_dist(2)
          do i3=1, n_bins(3)
            p(3) = min_p(3) + (i3-1)*sample_dist(3)
            histo_grid(:,i1,i2,i3) = p
          end do
        end do
      end do
    endif

    do at_i=1, at%N
      if (.not. mask_a(at_i)) cycle
      do i1=1, n_bins(1)
        p(1) = min_p(1) + (i1-1)*sample_dist(1)
        do i2=1, n_bins(2)
          p(2) = min_p(2) + (i2-1)*sample_dist(2)
          do i3=1, n_bins(3)
            p(3) = min_p(3) + (i3-1)*sample_dist(3)
            dist = distance_min_image(at,p,at%pos(:,at_i))
            histogram(i1,i2,i3) = histogram(i1,i2,i3) + exp(-0.5_dp*(dist/(gaussian_sigma))**2)*w
          end do ! i3
        end do ! i2
      end do ! i1
    end do ! at_i

    deallocate(mask_a)
  end subroutine sample_rectilinear_mesh_Gaussians

  subroutine calc_histogram(histogram, at, min_p, bin_width, n_bins, mask_str, histo_grid, accumulate)
    real(dp), intent(inout) :: histogram(:,:,:)
    type(Atoms), intent(in) :: at
    real(dp), intent(in) :: min_p(3), bin_width(3)
    integer, intent(in) :: n_bins(3)
    character(len=*), optional, intent(in) :: mask_str
    real(dp), intent(out), optional :: histo_grid(:,:,:,:)
    logical, optional, intent(in) :: accumulate

    logical :: my_accumulate
    integer :: i, bin(3)
    logical, allocatable :: mask_a(:)
    real(dp) :: p(3)

    my_accumulate = optional_default(.false., accumulate)
    if (.not. my_accumulate) histogram = 0.0_dp

    if (present(histo_grid)) then
      do i1=1, n_bins(1)
        p(1) = min_p(1) + (real(i1,dp)-0.5_dp)*bin_width(1)
        do i2=1, n_bins(2)
          p(2) = min_p(2) + (real(i2,dp)-0.5_dp)*bin_width(2)
          do i3=1, n_bins(3)
            p(3) = min_p(3) + (real(i3,dp)-0.5_dp)*bin_width(3)
            histo_grid(:,i1,i2,i3) = p
          end do
        end do
      end do
    endif

    allocate(mask_a(at%N))
    call is_in_mask(mask_a, at, mask_str)

    do i=1, at%N
      if (.not. mask_a(i)) cycle
      bin = floor((at%pos(:,i)-min_p)/bin_width)+1
      if (all(bin >= 1) .and. all (bin <= n_bins)) histogram(bin(1),bin(2),bin(3)) = histogram(bin(1),bin(2),bin(3)) + 1.0_dp
    end do

    deallocate(mask_a)

  end subroutine calc_histogram

  subroutine is_in_mask(mask_a, at, mask_str)
    type(Atoms), intent(in) :: at
    logical, intent(out) :: mask_a(at%N)
    character(len=*), optional, intent(in) :: mask_str

    integer :: i
    integer :: Zmask
    type(Table) :: atom_indices

    if (.not. present(mask_str)) then
      mask_a = .true.
      return
    endif

    if (len_trim(mask_str) == 0) then
      mask_a = .true.
      return
    endif

    mask_a = .false.
    if (mask_str(1:1)=='@') then
      call parse_atom_mask(mask_str,atom_indices)
      do i=1, atom_indices%N
        mask_a(atom_indices%int(1,i)) = .true.
      end do
    else if (scan(mask_str,'=')/=0) then
      call system_abort("property type mask not supported yet")
    else
      Zmask = Atomic_Number(mask_str)
      do i=1, at%N
        if (at%Z(i) == Zmask) mask_a(i) = .true.
      end do
    end if
  end subroutine is_in_mask

  subroutine reallocate_histos(histos, n, n_bins)
    real(dp), allocatable, intent(inout) :: histos(:,:,:,:)
    integer, intent(in) :: n, n_bins(3)

    integer :: new_size
    real(dp), allocatable :: t_histos(:,:,:,:)

    if (allocated(histos)) then
      if (n <= size(histos,4)) return
      allocate(t_histos(size(histos,1),size(histos,2),size(histos,3),size(histos,4)))
      t_histos = histos
      deallocate(histos)
      if (size(t_histos,4) <= 0) then
        new_size = 10
      else if (size(t_histos,4) < 1000) then
        new_size = 2*size(t_histos,4)
      else
        new_size = floor(1.25*size(t_histos,4))
      endif
      allocate(histos(size(t_histos,1), size(t_histos,2), size(t_histos,3), new_size))
      histos(1:size(t_histos,1),1:size(t_histos,2),1:size(t_histos,3),1:size(t_histos,4)) = &
        t_histos(1:size(t_histos,1),1:size(t_histos,2),1:size(t_histos,3),1:size(t_histos,4))
      histos(1:size(t_histos,1),1:size(t_histos,2),1:size(t_histos,3),size(t_histos,4)+1:size(histos,4)) = 0.0_dp
      deallocate(t_histos)
    else
      allocate(histos(n_bins(1), n_bins(2), n_bins(3), n))
    endif
  end subroutine reallocate_histos

  subroutine calc_mean_var_3array(a, mean, var)
    real(dp), intent(in) :: a(:,:,:,:)
    real(dp), intent(out) :: mean(:,:,:)
    real(dp), optional, intent(out) :: var(:,:,:)

    integer :: i1, i2, i3, n

    n = size(a,4)
    mean = sum(a, 4)/real(n,dp)
    if (present(var)) then
      do i1=1, size(a,1)
      do i2=1, size(a,2)
      do i3=1, size(a,3)
        var(i1,i2,i3) = sum((a(i1,i2,i3,:)-mean(i1,i2,i3))**2)/real(n,dp)
      end do
      end do
      end do
    end if
  end subroutine calc_mean_var_3array

end program density_new

max_stars_repo_name: Sideboard/QUIP

max_issues_repo_name: Sideboard/QUIP

max_forks_repo_name: Sideboard/QUIP


text:
!
!     (c) 2019 Guide Star Engineering, LLC
!     This Software was developed for the US Nuclear Regulatory Commission (US NRC)
!     under contract "Multi-Dimensional Physics Implementation into Fuel Analysis under
!     Steady-state and Transients (FAST)", contract # NRC-HQ-60-17-C-0007
!
!
!    NEMO - Numerical Engine (for) Multiphysics Operators
! Copyright (c) 2007, Stefano Toninel
!                     Gian Marco Bianchi  University of Bologna
!              David P. Schmidt    University of Massachusetts - Amherst
!              Salvatore Filippone University of Rome Tor Vergata
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification,
! are permitted provided that the following conditions are met:
!
!     1. Redistributions of source code must retain the above copyright notice,
!        this list of conditions and the following disclaimer.
!     2. Redistributions in binary form must reproduce the above copyright notice,
!        this list of conditions and the following disclaimer in the documentation
!        and/or other materials provided with the distribution.
!     3. Neither the name of the NEMO project nor the names of its contributors
!        may be used to endorse or promote products derived from this software
!        without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
! ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
! WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
! ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
!---------------------------------------------------------------------------------
!
! $Id: tools_nemo.f90 3322 2008-08-28 13:23:56Z sfilippo $
!
! Description:
!    Provides TOOLS_* and CLASS_* modules for NEMO-based high-end applications
!
MODULE tools_nemo

    USE class_psblas
    !
    USE class_connectivity
    USE class_dimensions
    USE class_discretization
    USE class_scalar_field
    USE class_vector_field
    USE class_iterating
    USE class_material
    USE class_mesh
    USE class_output
    USE class_scalar_pde
    USE class_vector_pde
    USE class_scalar_source
    USE class_surface
    USE class_vector
    USE class_bc
    !
    USE op_ddt
    USE op_d2dt2
    USE op_div
    USE op_field
    USE op_grad
    USE op_laplacian
    USE op_source
    !
    USE tools_input
    USE tools_math
    USE tools_mesh_check
    USE tools_mesh_move
    USE tools_mesh_optimize
    USE tools_operators
    USE tools_output_basics
    USE tools_output

END MODULE tools_nemo

max_stars_repo_name: sourceryinstitute/MORFEUS-Source

max_issues_repo_name: sourceryinstitute/MORFEUS-Source

max_forks_repo_name: sourceryinstitute/MORFEUS-Source


program anaice2rst
!!$ /////////////////////////////////////////////////////////////////////////
!!$ *
!!$ * Title: anaice2rst.f90
!!$ *
!!$ * Description: Based on Bin Zhao's code that converts PIOMASS 12 itd into CICE 5 itd.
!!$ *              The core philosophy is to keep the background itd and conserve the Enthalpy per unit volume.
!!$ *
!!$ * Options:
!!$ *              -rst     : File name for saltwater restart.
!!$ *                         Default: Salt_fname  = 'saltwater_internal_rst.nc4'
!!$ *              -orst    : File name for ocean temp salt restart.
!!$ *                         Default: Ocean_fname  = 'ocean_temp_salt.res.nc'
!!$ *              -ana     : File name containing analysis, or sea-ice state on ocean grid.
!!$ *                         variables in the file should be AICE and/or HICE and/or HSNO on the ocean native grid.
!!$ *                         Default: Ana_fname  = 'ana.nc4'
!!$ *              -ana_fmt : Format of Ana_fname, 'oletkf' or 'regular'
!!$ *                         Default: ana_fmt='regular'
!!$ *              -do_aice : DO_AICE_ANA, Read AICE, .true or .false.  !!!!! Assumes analysis file contains logit(AICE)
!!$ *              -do_hice : DO_HICE_ANA,   "  HICE,   "
!!$ *              -do_hsno : DO_HSNO_ANA,   "  HSNO,   "
!!$ *              -do_all_1st : DO_ALL_1ST, .true or .false. !! apply increment starting from 1st category; Default(.false.)
!!$ *               NOTE: if true, this scheme creates some biases because ice may not be present in the 1st category when removed
!!$ *              -hi_new : real, new ice thickness when DO_ALL_1ST is .true.;
!!$ *                        Default(0.5m)
!!$ *
!!$ * Output:
!!$ *              saltwater_internal_rst_OUT.nc4
!!$ * Example:
!!$ *
!!$ *             saltwater_internal_rst_OUT is a copy of saltwater_internal_rst
!!$ *             ./anaice2rst.x saltwater_internal_rst_OUT saltwater_internal_rst ana-20030407.nc .T. .F. .F.
!!$ *                 will read AICE and rebalance CICE's restart, result written in saltwater_internal_rst_OUT
!!$ *

  use netcdf
  use ocnice_utils

  implicit none

  integer, parameter :: real_kind = selected_real_kind(6)
  integer, parameter :: dbl_kind = selected_real_kind(13)
  integer :: i, j, n, nt, k, argc, fid_in, arg_index, Nargs
  character*1028 :: Salt_fname  = 'saltwater_internal_rst.nc4'
  character*1028 :: Ocean_fname  = 'ocean_temp_salt.res.nc'
  character*1028 :: Ana_fname  = 'ana.nc4'
  character*1028 :: Outfile  = 'saltwater_internal_rst_OUT.nc4'
  character*1028 :: OceanOutfile  = 'ocean_temp_salt.res.OUT.nc'
  character*400  :: tilefile_fname  = 'tile.data'
  logical       :: DO_AICE_ANA = .true.
  logical       :: DO_HICE_ANA = .true.
  logical       :: DO_HSNO_ANA = .true.
  logical       :: DO_SKIN_FLUSH = .false.
  logical       :: DO_ALL_1ST  = .false.
  logical       :: DO_ADHOC_TW = .false.
  logical       :: DO_CONS_VOL = .false.
  logical       :: oob
  character*128 :: ana_fmt = 'regular' ! 'oletkf' or 'regular'

  character*128 :: fname
  integer :: varid, varidvi, varidlon, varidlat
  character*128 :: varname

  real*8, allocatable, dimension(:)    :: tmpx
  real*8, allocatable, dimension(:, :) :: aicen, vicen ,vsnon
  real*4, allocatable, dimension(:,:)    :: Tw  ! Skin temperature (under ice) [K]
  real*4, allocatable, dimension(:,:)    :: Sw  ! Skin salinity (under ice)    [psu]
  real*4, allocatable, dimension(:,:)    :: Slmask

  real*8, allocatable, dimension(:, :, :) ::   eicen    ! energy of melting for each ice layer  (J/m^2)
  real*8, allocatable, dimension(:, :, :) ::   esnon    ! energy of melting for each snow layer (J/m^2)

  real*4, allocatable, dimension(:,  :)    :: sic_ana, sit_ana, snd_ana
  real*8, allocatable, dimension(:,  :)    :: sst_incr       ! Re-analyzed sst increment
  real*4, allocatable, dimension(:,  :, :) :: Tocean      ! 3d ocean potential temperature [C]
  real*4, allocatable, dimension(:,  :, :) :: Socean      ! 3d ocean salinity [PSU]
  real*4, allocatable, dimension(:,  :, : , :)    :: DUMVAR

  real*4, allocatable, dimension(:) :: madjice
  real*4, allocatable, dimension(:, :) :: madjiceij, Toij

  real*4            :: beta                    ! Scaling of the background field
  real*4            :: mass_incr               ! Sea-ice mass increment
  real*4            :: sic_bkg                 ! Background ice Fraction
  real*4            :: sit_bkg                 ! Background ice thickness
  real*4            :: snd_bkg                 ! Background snow depth
  real*8            :: sic_incr                ! Ice Fraction increment
  real*4            :: sic_incr_max = 0.9      ! Max ice Fraction increment
  real*4            :: Tliquid                 ! Background Temperature of liquid ocean
  real*4            :: Tf                      ! Freeze temp
  real*4            :: Tincr                   ! Tw increment
  real*4, parameter :: Tincr_max = 1.0         ! Max Tw increment
  real*4, parameter :: rho_ice = 917.0         ! Sea-ice density [kg/m3]
  real*4, parameter :: rho_sno = 330.0         ! Snow density [kg/m3]
  real*4, parameter :: tiny = 1.0e-7           !
  real*4, parameter :: puny = 1.0e-11          !
  real*4, parameter :: beta_max = 2.0         !
  real*4, parameter :: rho_lwe = 1003.091    ! Surface water density. S=4 psu, T=-0.216 C (freeze point) [kg/m3]

  real (kind=dbl_kind), parameter :: &
         p5   = 0.5_dbl_kind, &
         c0   = 0.0_dbl_kind, &
         c1   = 1.0_dbl_kind, &
         c1p5 = 1.5_dbl_kind, &
         c2   = 2.0_dbl_kind, &
         c3   = 3.0_dbl_kind, &
         c15  = 15.0_dbl_kind, &
         pi = 3.14159265358979323846_dbl_kind ! pi

  real(kind=dbl_kind), parameter  :: &
                         rhoi      = 917.0_dbl_kind , &
                         rhow      = 1026.0_dbl_kind, &
                         rhos      = 330.0_dbl_kind , &
                         cp_ice    = 2106._dbl_kind   ,&! specific heat of fresh ice (J/kg/K)
                         cp_ocn    = 4218._dbl_kind  ,& ! specific heat of ocn    (J/kg/K)
                         depressT  = 0.054_dbl_kind   ,&! Tf:brine salinity ratio (C/ppt)
                         saltmax   = 3.2_dbl_kind,  & ! max salinity at ice base (ppt)
                         nsal      = 0.407_dbl_kind, &
                         msal      = 0.573_dbl_kind, &
                         Lsub      = 2.835e6_dbl_kind ,&! latent heat, sublimation freshwater (J/kg)
                         Lvap      = 2.501e6_dbl_kind ,&! latent heat, vaporization freshwater (J/kg)
                         Lfresh    = Lsub-Lvap        ! latent heat of melting of fresh ice (J/kg)

  character*128        :: Usage="Example: anaice2rst.x saltwater_internal_rst_OUT saltwater_internal_rst_IN ANALYSIS.nc .T. .T. .T."
  integer, parameter   :: ncat = 5
  integer, parameter   :: nilyr = 4
  integer, parameter   :: nslyr = 1
  integer, allocatable :: GIs(:), GJs(:)
  integer              :: targt = 0
  integer              :: tilemax

  real*4               :: hinew = 0.5             ! new ice thickness generated
  real*8               :: hi, hs, qi0, incr_remain
  real*4               :: mass1, mass2
  real*8               :: qin_save(nilyr), qsn_save(nslyr)
  real (kind=real_kind), dimension(nilyr+1) :: &
                       salin       , & ! salinity (ppt)
                       Tmlt            ! melting temp, -depressT * salinity
                         ! nilyr + 1 index is for bottom surface
  real (kind=dbl_kind)    :: zn       ! normalized ice thickness
  real (kind=dbl_kind)    :: cc1, cc2, cc3, x1, rncat, hin
  real (kind=dbl_kind) :: &
         hin_max(0:ncat) ! category limits (m)
  real (kind=dbl_kind) :: &
         hsn

  integer ::  ilyr1 (ncat), & ! array position of top ice layer in each cat
              ilyrn (ncat), & ! array position of bottom ice layer in each cat
              slyr1 (ncat), & ! array position of top snow layer in each cat
              slyrn (ncat)    ! array position of bottom snow layer in each cat

  character*300                          :: buffer, command, arg1, arg2
  integer :: ncid, dimids2d(2), cat_dimid, catp1_dimid, tile_dimid, dimids1d(1)

  type (ModelGrid)                      :: ogrid     ! Ocean grid, tripolar
  type (OceanRaster)                    :: oraster   ! Ocean tiles

  !PROCESS COMMAND LINE
  !====================
  Nargs = iargc()
  do arg_index = 1, Nargs, 2
     call get_command_argument(arg_index, arg1)
     select case(arg1)
     case('-rst')
        call get_command_argument(arg_index+1, arg2)
        Salt_fname = arg2
     case('-orst')
        call get_command_argument(arg_index+1, arg2)
        Ocean_fname = arg2
     case('-ana')
        call get_command_argument(arg_index+1, arg2)
        Ana_fname = arg2
     case('-tilefile')
        call get_command_argument(arg_index+1, arg2)
        tilefile_fname = arg2
     case('-do_aice')
        call get_command_argument(arg_index+1, arg2)
        DO_AICE_ANA = str2bool(arg2)
     case('-do_hice')
        call get_command_argument(arg_index+1, arg2)
        DO_HICE_ANA = str2bool(arg2)
     case('-do_hsno')
        call get_command_argument(arg_index+1, arg2)
        DO_HSNO_ANA = str2bool(arg2)
     case('-do_all_1st')
        call get_command_argument(arg_index+1, arg2)
        DO_ALL_1ST = str2bool(arg2)
     case('-do_skin_flush')
        call get_command_argument(arg_index+1, arg2)
        DO_SKIN_FLUSH = str2bool(arg2)
     case('-do_conserve_volume')
        call get_command_argument(arg_index+1, arg2)
        DO_CONS_VOL = str2bool(arg2)
     case('-hi_new')
        call get_command_argument(arg_index+1, arg2)
        read(arg2, *) hinew
     case('-target_tile')
        call get_command_argument(arg_index+1, arg2)
        read(arg2, *) targt
     case('-do_adhoc_tw_adjustment')
        call get_command_argument(arg_index+1, arg2)
        DO_ADHOC_TW = str2bool(arg2)
     case('-ana_fmt')
        call get_command_argument(arg_index+1, arg2)
        ana_fmt = arg2
     end select
  end do

  command = "cp "//trim(Salt_fname)//" "//trim(Outfile)
  call execute_command_line (command, exitstat=i)

  command = "cp "//trim(Ocean_fname)//" "//trim(OceanOutfile)
  call execute_command_line (command, exitstat=i)

  !RASTER INFO
  !===========
  call readraster(oraster)

  !OCEAN GRID
  !==========
  call read_grid(ogrid)

  print *,'ocean tiles area = ',sum(oraster%areas)
  print *,'ocean grid area = ',sum(ogrid%wet*ogrid%area)/(1000.0**2)

  !-----------------------------------------------------------------
  ! vectors identifying first and last layer in each category
  !-----------------------------------------------------------------
  ilyr1(1) = 1                       ! if nilyr  = 4
  ilyrn(1) = nilyr                   !   ilyr1 = { 1,5,9 }
  do n = 2, ncat                     !   ilyrn = { 4,8,12} etc
       ilyr1(n) = ilyrn(n-1) + 1
       ilyrn(n) = ilyrn(n-1) + nilyr
  enddo

  slyr1(1) = 1
  slyrn(1) = nslyr
  do n = 2, ncat
       slyr1(n) = slyrn(n-1) + 1
       slyrn(n) = slyrn(n-1) + nslyr
  enddo

  do k = 1, nilyr
      zn = (real(k,kind=dbl_kind)-p5) /  &
            real(nilyr,kind=dbl_kind)
      salin(k)=(saltmax/c2)*(c1-cos(pi*zn**(nsal/(msal+zn))))
!            salin(k)=saltmax ! for isosaline ice
      Tmlt(k) = -salin(k)*depressT
   enddo
   salin(nilyr+1) = saltmax
   Tmlt(nilyr+1) = -salin(nilyr+1)*depressT

    ! linear remapping itd category limits
    rncat = real(ncat, kind=dbl_kind)
    cc1 = c3/rncat
    cc2 = c15*cc1
    cc3 = c3
    hin_max(0) = c0     ! minimum ice thickness, m
    do n = 1, ncat
        x1 = real(n-1,kind=dbl_kind) / rncat
        hin_max(n) = hin_max(n-1) &
                    + cc1 + cc2*(c1 + tanh(cc3*(x1-c1)))
    enddo

  !READ ANALYSIS
  !=============
  if (DO_AICE_ANA) then
     print *,'Rebalancing state for sea-ice fraction'
     allocate(sic_ana(ogrid%Nx, ogrid%Ny))
     varname = 'AICE'
     call readseaicerst(Ana_fname, varname, ogrid, ana_fmt, sic_ana)
     where (sic_ana<0.0)
        sic_ana=0.0
     end where
     where (sic_ana>=1.0-tiny)
        sic_ana=1.0-tiny
     end where
  end if
  if (DO_HICE_ANA) then
     print *,'Rebalancing state for sea-ice thickness'
     allocate(sit_ana(ogrid%Nx, ogrid%Ny))
     varname = 'HICE'
     call readseaicerst(Ana_fname, varname, ogrid, ana_fmt, sit_ana)
     where (sit_ana<0.0)
        sit_ana=0.0
     end where
  end if
  if (DO_HSNO_ANA) then
     print *,'Rebalancing state for snow depth over sea-ice'
     allocate(snd_ana(ogrid%Nx, ogrid%Ny))
     varname = 'HSNO'
     call readseaicerst(Ana_fname, varname, ogrid, ana_fmt, snd_ana)
     where (snd_ana<0.0)
        snd_ana=0.0
     end where
  end if

  !READ SALTWATER RESTART
  !======================
  allocate( aicen(oraster%Ntiles,ncat+1) )
  varname = 'FR'
  call readrst2d(Salt_fname, varname, aicen)

  allocate(Tw(oraster%Ntiles, 1))
  varname = 'TSKINW'
  call readrst2d(Salt_fname, varname, Tw)

  allocate(Sw(oraster%Ntiles, 1))
  varname = 'SSKINW'
  call readrst2d(Salt_fname, varname, Sw)

  allocate(vicen(oraster%Ntiles,ncat))
  varname = 'VOLICE'
  call readrst2d(Salt_fname, varname, vicen)

  allocate(vsnon(oraster%Ntiles,ncat))
  varname = 'VOLSNO'
  call readrst2d(Salt_fname, varname, vsnon)

  allocate(eicen(oraster%Ntiles,nilyr,ncat))
  varname = 'ERGICE'
  call readrst3d(Salt_fname, varname, eicen)

  allocate(esnon(oraster%Ntiles,nslyr,ncat))
  varname = 'ERGSNO'
  call readrst3d(Salt_fname, varname, esnon)

  allocate(Slmask(oraster%Ntiles, 1))
  varname = 'SLMASK'
  call readrst2d(Salt_fname, varname, Slmask)

  allocate(madjice(oraster%Ntiles))
  allocate(madjiceij(ogrid%Nx, ogrid%Ny))
  allocate(Toij(ogrid%Nx, ogrid%Ny))
  madjice   = 0.0
  madjiceij = 0.0
  Toij      = 0.0

  !READ OCEAN TEMP SALT RESTART
  !============================
  varname =  'temp'
  allocate(Tocean(ogrid%Nx, ogrid%Ny, ogrid%Nz))
  call readoceanrst(Ocean_fname, varname, ogrid, Tocean)

  varname =  'salt'
  allocate(Socean(ogrid%Nx, ogrid%Ny, ogrid%Nz))
  call readoceanrst(Ocean_fname, varname, ogrid, Socean)

  !ANALYSIS UPDATE
  !===============
  ! The following conserves the background ice thickness distribution
  ! and rebalances the ice/snow internal energy so that Enthalpy/unit volume
  ! is conserved.
  ! "Craming" the analysis in the background restart is done incrementally:
  !    1 - Update aicen                  => balance ice and snow internal E
  !    2 - If ice thickness was analyzed => balance ice E
  !    3 - If snow depth was analyzed    => balance snow E
  !

  allocate(sst_incr(ogrid%Nx, ogrid%Ny))
  do i=1,oraster%Ntiles
     if(i==targt) then
         print*, 'lon, lat: ', oraster%lons(i), oraster%lats(i)
         print*, 'initial aice0: ', aicen(i,1)
         print*, 'initial aicen: ', (aicen(i,n), n=2,ncat+1)
         print*, 'initial vicen: ', (vicen(i,n), n=1,ncat)
     endif
  enddo

  if(DO_SKIN_FLUSH) then
     do i=1,oraster%Ntiles
     ! first flush skin layer with MOM top layer T & S to avoid
     ! initial shock
       if(i==targt) then
             print*, 'tripole i,j: ', oraster%GIs(i), oraster%GJs(i)
             print*, 'before flushing: ', Tw(i,1), Tocean(oraster%GIs(i), oraster%GJs(i),1)
             print*, 'before flushing: ', Sw(i,1)
       endif
       if(SLMASK(i,1) < 0.5) then ! only flush when it is a MOM tile
          SW(i,1) = Socean(oraster%GIs(i), oraster%GJs(i),1)
          TW(i,1) = Tocean(oraster%GIs(i), oraster%GJs(i),1)+273.15
       endif
       if(i==targt) then
             print*, 'after flushing: ', Tw(i,1), Tocean(oraster%GIs(i), oraster%GJs(i),1)
             print*, 'after flushing: ', Sw(i,1)
       endif
     enddo
  endif

  mass_incr = 0.0
  mass1     = 0.0
  do i=1,oraster%Ntiles
      if(SLMASK(i,1) > 0.5) cycle ! no blending for salt lakes
      if(SLMASK(i,1) < 0.5) then ! only count MOM tiles
        mass1 = mass1 + oraster%areas(i)*(rho_ice*sum(vicen(i,:))+rho_sno*sum(vsnon(i,:)))
      endif

     ! SEA-ICE FRACTION
     !-----------------
     if (DO_AICE_ANA) then
        sic_bkg = sum( aicen(i,2:) )
        sic_incr = sic_ana(oraster%GIs(i), oraster%GJs(i))-sic_bkg
        !if(sic_bkg > 0.9 .and. sic_ana(oraster%GIs(i), oraster%GJs(i)) > 0.9) cycle
        if(sic_bkg > puny .and. sic_ana(oraster%GIs(i), oraster%GJs(i)) > puny .and. &
           abs(sic_bkg-sic_ana(oraster%GIs(i), oraster%GJs(i))) < 0.01) cycle
        if (DO_ALL_1ST) then
           if(sic_ana(oraster%GIs(i), oraster%GJs(i)) < -puny .or. &
              sic_ana(oraster%GIs(i), oraster%GJs(i)) > 1.0+puny) cycle
           !safety check
           oob = .false.
           do n=1, ncat
               if(aicen(i,n+1) > puny) then
                 hi = vicen(i,n) / aicen(i,n+1)
                 if(hi > hin_max(n) .or. hi < hin_max(n-1)) then
                    oob = .true.
                    exit
                 endif
               endif
           enddo
           if(oob) cycle
           Tf = -depressT * Sw(i,1) ! Linear liquidus approximation (Bitz and Lipscomb, 1999)
           !if(sic_bkg>0.0) then
              if(sic_incr > 0.0) then  ! generate ice
                if(aicen(i,2) > puny) then
                   call  save_old_state(aicen(i,:),vicen(i,:),vsnon(i,:),&
                                        eicen(i,:,:),esnon(i,:,:), &
                                        nilyr, nslyr, 2, 1, &
                                        hi, hs, qin_save, qsn_save)
                   call  restore_ice_state(aicen(i,:),vicen(i,:),vsnon(i,:),&
                                          eicen(i,:,:),esnon(i,:,:), &
                                          nilyr, nslyr, 2, 1, sic_incr, &
                                          hi, hs, qin_save, qsn_save, himin=hinew)
                else
                   aicen(i,2) = aicen(i,2) + sic_incr
                   vicen(i,1) = aicen(i,2) * hinew
                   ! new ice enthalpy
                   do k=1,nilyr
                        qi0 = -rhoi * (cp_ice*(Tmlt(k)-Tf)   &
                           + Lfresh*(1.-Tmlt(k)/Tf) - cp_ocn*Tmlt(k))
                        eicen(i,k,1)                          &
                              = qi0 * vicen(i,1)                  &
                               / real(nilyr,kind=real_kind)
                   enddo
                endif
              else ! remove ice
                  incr_remain = sic_incr
                  do n=2, ncat+1
                     !if(incr_remain < 0.0) then
                        if(aicen(i,n) > puny) then
                          call  save_old_state(aicen(i,:),vicen(i,:),vsnon(i,:),&
                                               eicen(i,:,:),esnon(i,:,:), &
                                               nilyr, nslyr, n, n-1, &
                                               hi, hs, qin_save, qsn_save)
                           incr_remain = aicen(i,n) + sic_incr
                           if(i==targt) then
                             print*, 'during processing: '
                             print*, i, n,' incr = ', sic_incr
                             print*, i, n,' remain incr = ', incr_remain
                           endif
                           if(incr_remain > puny) then
                              call  restore_ice_state(aicen(i,:),vicen(i,:),vsnon(i,:),&
                                                      eicen(i,:,:),esnon(i,:,:), &
                                                      nilyr, nslyr, n, n-1, sic_incr, &
                                                      hi, hs, qin_save, qsn_save)
                              exit
                           else
                              ! all ice in category n is gone
                              aicen(i,n)   = 0.0
                              vicen(i,n-1) = 0.0
                              vsnon(i,n-1) = 0.0
                              eicen(i,:,n-1) = 0.0
                              esnon(i,:,n-1) = 0.0
                           endif
                        endif
                        sic_incr = incr_remain
                  enddo ! n
                endif ! generate ice or not
           !endif
           aicen(i,1) = 1.0 - sum( aicen(i,2:) )       ! Re-compute ocean fraction
           if(i==targt) then
              print*, 'after processing: '
              print*, i, oraster%lons(i), oraster%lats(i)
              print*, aicen(i,1), aicen(i,2), aicen(i,3)
              print*, aicen(i,4), aicen(i,5), aicen(i,6)
           endif
        else  ! .not. DO_ALL_1ST
        if (sic_bkg>0.0) then
           beta = real( (sic_ana(oraster%GIs(i), oraster%GJs(i))+tiny)/sic_bkg ,kind=real_kind)
           beta = min(beta, beta_max)
           !print*, beta
!           if(beta <= beta_max) then
              if(i==targt) then
                 print*, 'beta : ', beta
                 print*, 'before scaling vicen: ', (vicen(i,n), n=1,ncat)
              endif
              aicen(i,:) = beta * aicen(i,:)         ! Re-adjust background fraction to match analysis
              mass_incr = mass_incr + (beta-1.0)*oraster%areas(i)*( rho_ice*sum( vicen(i,:) )  + rho_sno*sum( vsnon(i,:) ) ) ! Update sea-ice mass increment
                                                                                ! due to sea-ice fraction analysis
                  vicen(i,:)   = beta * vicen(i,:)       ! Adjust volume of ice
                  vsnon(i,:)   = beta * vsnon(i,:)       ! and snow accordingly.
                  eicen(i,:,:) = beta * eicen(i,:,:)     ! This way of rebalancing the ice
                  esnon(i,:,:) = beta * esnon(i,:,:)     ! and snow internal energy conserves the enthalpy per unit volume
              if(i==targt) then
                 print*, 'after scaling vicen: ', (vicen(i,n), n=1,ncat)
              endif
        elseif (sic_ana(oraster%GIs(i), oraster%GJs(i)) > puny) then
           !beta=1.0
           if(i==targt) then
                 print*, 'sic_ana : ', sic_ana(oraster%GIs(i), oraster%GJs(i))
               print*, 'in step C before vicen: ', (vicen(i,n), n=1,ncat)
           endif
           Tf = -depressT * Sw(i,1) ! Linear liquidus approximation (Bitz and Lipscomb, 1999)
           aicen(i,2) = aicen(i,2) + sic_ana(oraster%GIs(i), oraster%GJs(i))
           vicen(i,1) = aicen(i,2) * hinew
           ! new ice enthalpy
           do k=1,nilyr
              qi0 = -rhoi * (cp_ice*(Tmlt(k)-Tf)   &
                   + Lfresh*(1.-Tmlt(k)/Tf) - cp_ocn*Tmlt(k))
              eicen(i,k,1)                          &
                   = qi0 * vicen(i,1)                  &
                   / real(nilyr,kind=real_kind)
           enddo
           if(i==targt) then
               print*, 'in step C after vicen: ', (vicen(i,n), n=1,ncat)
           endif
        end if
        aicen(i,1) = 1.0 - sum( aicen(i,2:) )       ! Re-compute ocean fraction

        end if ! DO_ALL_1ST

        Tf = -0.054 * Sw(i,1) + 273.15

        if( DO_ADHOC_TW ) then
           if ( sum( aicen(i,2:) ) > 0.01 .and. sic_bkg < 0.01 ) then
              Tw(i,1) = Tf
              madjice(i) = oraster%wghts(i)
              madjiceij(oraster%GIs(i), oraster%GJs(i)) = &
                        madjiceij(oraster%GIs(i), oraster%GJs(i)) + madjice(i)
              Toij(oraster%GIs(i), oraster%GJs(i)) = &
                        Toij(oraster%GIs(i), oraster%GJs(i)) + madjice(i)*Tw(i,1)
           endif
           if ( sum( aicen(i,2:) ) < 0.01 .and. sic_bkg > 0.01 ) then
              Tw(i,1) = Tf + 0.5
              madjice(i) = oraster%wghts(i)
              madjiceij(oraster%GIs(i), oraster%GJs(i)) = &
                        madjiceij(oraster%GIs(i), oraster%GJs(i)) + madjice(i)
              Toij(oraster%GIs(i), oraster%GJs(i)) = &
                        Toij(oraster%GIs(i), oraster%GJs(i)) + madjice(i)*Tw(i,1)
           endif
        else
        ! the following adjustment should be done regardless
        if ( (sum( aicen(i,2:) ) > 0.01) .and. (sic_bkg>0.01) .and. (sum( aicen(i,2:))<0.99) .and. (sic_bkg<0.99) ) then
           Tliquid= ( Tw(i,1) - sic_bkg*Tf )/(1.0-sic_bkg)            ! Keep liquid water temp unchanged
           if(i==targt) then
              print*, Tf, Tliquid, Tw(i,1), sic_bkg
           endif
           Tincr = sum( aicen(i,2:) ) * Tf + ( 1.0 - sum( aicen(i,2:) ) ) * Tliquid - Tw(i,1)
           Tincr = Tincr_max * tanh( Tincr/Tincr_max ) ! limit the increment to ]-Tincr_max, Tincr_max[
           Tw(i,1) = Tw(i,1) + Tincr
           !Tw(i,1) = sum( aicen(i,2:) ) * Tf + ( 1.0 - sum( aicen(i,2:) ) ) * Tliquid
           if(i==targt) then
              print*,  sum(aicen(i,2:)), (1.0 - sum( aicen(i,2:))), Tw(i,1), Tincr
           endif
           madjice(i) = oraster%wghts(i)
           madjiceij(oraster%GIs(i), oraster%GJs(i)) = &
                        madjiceij(oraster%GIs(i), oraster%GJs(i)) + madjice(i)
           Toij(oraster%GIs(i), oraster%GJs(i)) = &
                        Toij(oraster%GIs(i), oraster%GJs(i)) + madjice(i)*Tw(i,1)
           !sst_incr(oraster%GIs(i), oraster%GJs(i)) = Tincr    ! Save to diagnoze
           !Tocean(oraster%GIs(i), oraster%GJs(i), 1) = Tw(i,1) - 273.15 ! Zero out ice-ocean fluxes where sea-ice is present, investigate propagating the
                                                                        ! increment further down the water column

           !Tw(i,1) = sum( aicen(i,2:) ) * Tf + ( 1.0 - sum( aicen(i,2:) ) ) * Tliquid  ! Adjust skin temperature below sea-ice
           !write(55,*)Tw(i,1)-273.15, Tliquid - 273.15, Tf - 273.15, sic_bkg, Tincr, sum( aicen(i,2:) )
        end if
        end if
     end if

     ! SEA-ICE THICKNESS
     !------------------
     if (DO_HICE_ANA) then
     if(i==targt) then
         print*, 'lon, lat: ', oraster%lons(i), oraster%lats(i)
         print*, 'before ana_hice aice0: ', aicen(i,1)
         print*, 'before ana_hice aicen: ', (aicen(i,n), n=2,ncat+1)
         print*, 'before ana_hice vicen: ', (vicen(i,n), n=1,ncat)
     endif
        sit_bkg = sum( vicen(i,:) )
        if (sit_bkg>0.0) then
           beta=real( (sit_ana(oraster%GIs(i), oraster%GJs(i))+tiny)/sit_bkg ,kind=real_kind)
        else
           beta=1.0
        end if
        mass_incr = mass_incr + (beta-1.0)*oraster%areas(i)*rho_ice*sum( vicen(i,:) )   ! Update sea-ice mass increment
                                                                                ! due to sea-ice thickness analysis

        vicen(i,:)   = beta * vicen(i,:)            ! Readjust background volume
        eicen(i,:,:) = beta * eicen(i,:,:)          ! and ice energy
     if(i==targt) then
         print*, 'lon, lat: ', oraster%lons(i), oraster%lats(i)
         print*, 'after ana_hice aice0: ', aicen(i,1)
         print*, 'after ana_hice aicen: ', (aicen(i,n), n=2,ncat+1)
         print*, 'after ana_hice vicen: ', (vicen(i,n), n=1,ncat)
     endif
     end if

     ! SNOW DEPTH
     !-----------
     if (DO_HSNO_ANA) then
        snd_bkg = sum( vsnon(i,:) )
        if (snd_bkg>0.0) then
           beta=real( (snd_ana(oraster%GIs(i), oraster%GJs(i))+tiny)/snd_bkg ,kind=real_kind)
        else
           beta=1.0
        end if

        mass_incr = mass_incr + (beta-1.0)*oraster%areas(i)*rho_sno*sum( vsnon(i,:) )   ! Update sea-ice mass increment
                                                                                ! due to sea-ice thickness analysis

        vsnon(i,:)   = beta * vsnon(i,:)            ! Readjust background snow volume
        esnon(i,:,:) = beta * esnon(i,:,:)          ! and snow energy
     end if
  enddo

  if (DO_AICE_ANA .and. DO_ADHOC_TW) then

  do j=1,ogrid%Ny
  do i=1,ogrid%Nx
     if(madjiceij(i,j) /= 0.0) then
        Tocean(i,j,1) = Toij(i,j)/madjiceij(i,j) - 273.15 !
     endif
  enddo
  enddo

  endif

  mass2 = 0.0
  do i=1,oraster%Ntiles
      if(SLMASK(i,1) < 0.5) then ! only count MOM tiles
        mass2 = mass2 + oraster%areas(i)*(rho_ice*sum(vicen(i,:))+rho_sno*sum(vsnon(i,:)))
      endif
  enddo

  !allocate(tmpx(oraster%Ntiles))
  !tmpx = sum(aicen(:,2:), dim=2)
  !tilemax = maxloc(tmpx,dim=1)
  !print*, 'max at tile: ', tilemax
  !do n=1,ncat
  !  print*, aicen(tilemax,n+1), vicen(tilemax,n)
  !enddo
  !do i=1,oraster%Ntiles
  !  if(abs(oraster%lons(i)) < 0.5 .and.  abs(oraster%lats(i)-(-60.0)) < 0.5) then
  !    print*, i, oraster%lons(i), oraster%lats(i), tmpx(i)
  !  endif
  !enddo
  !deallocate(tmpx)

 ! if(DO_ALL_1ST) then
     print *, 'before mass : ', mass1
     print *, 'after  mass : ', mass2
     print *,'Sea-ice mass increment = ',mass2-mass1,' kg'
     mass_incr = mass2-mass1
 ! else
 !    print *,'Sea-ice mass increment = ',mass_incr,' kg'
     print *,'Sea-ice mass increment/m^2 = ',mass_incr/(sum(ogrid%wet*ogrid%area)),' kg/m2'
     print *,'Sea-ice mass increment in lwe = ',1000.0*mass_incr/(sum(ogrid%wet*ogrid%area)),' m'

   mass_incr = mass_incr/(rho_lwe*(sum(ogrid%wet*ogrid%area))) !Compute Eustatic sea-level increment
   open(10, file='seaice_mass_incr.txt',status='replace')
   write(10,*) mass_incr
   close(10)

  !endif

  !SAVE IN OCEAN TEMP SALT RESTART
  !===============================

  varname = 'temp'
  allocate( DUMVAR(ogrid%Nx, ogrid%Ny, ogrid%Nz, 1) )
  DUMVAR(:,:,:,1) = Tocean
  call check(nf90_open(OceanOutFile, NF90_WRITE, fid_in) )
  call check(nf90_inq_varid(fid_in, varname, varid))
  call check(nf90_put_var(fid_in, varid, DUMVAR))
  call check(nf90_close(fid_in))

  !SAVE IN SALTWATER RESTART
  !=========================

  varname = 'FR'
  call check(nf90_open(OutFile, NF90_WRITE, fid_in) )
  call check(nf90_inq_varid(fid_in, varname, varid))
  call check(nf90_put_var(fid_in, varid, aicen))

  varname = 'TSKINW'
  call check(nf90_inq_varid(fid_in, varname, varid))
  call check(nf90_put_var(fid_in, varid, Tw(:,1)))

  varname = 'SSKINW'
  call check(nf90_inq_varid(fid_in, varname, varid))
  call check(nf90_put_var(fid_in, varid, Sw(:,1)))

  varname = 'VOLICE'
  call check(nf90_inq_varid(fid_in, varname, varid))
  call check(nf90_put_var(fid_in, varid, vicen))

  varname = 'VOLSNO'
  call check(nf90_inq_varid(fid_in, varname, varid))
  call check(nf90_put_var(fid_in, varid, vsnon))

  varname = 'ERGICE'
  call check(nf90_inq_varid(fid_in, varname, varid))
  call check(nf90_put_var(fid_in, varid, eicen))

  varname = 'ERGSNO'
  call check(nf90_inq_varid(fid_in, varname, varid))
  call check(nf90_put_var(fid_in, varid, esnon))
  call check(nf90_close(fid_in))

  !SAVE TILE LOCATION
  !==================
  fname='tileloc.nc'
  call check( nf90_create(fname, NF90_CLOBBER, ncid) )
  call check( nf90_def_dim(ncid, "ncats", ncat, cat_dimid) )
  call check( nf90_def_dim(ncid, "ntiles", oraster%Ntiles, tile_dimid) )
  dimids1d =  (/ tile_dimid /)

  varname = 'lon'
  call check( nf90_def_var(ncid, varname, NF90_REAL, dimids1d, varidlon) )
  varname = 'lat'
  call check( nf90_def_var(ncid, varname, NF90_REAL, dimids1d, varidlat) )
  call check( nf90_enddef(ncid) )

  varname = 'lon'
  call check( nf90_inq_varid(ncid,varname,varidlon) )
  call check( nf90_put_var(ncid, varidlon, oraster%lons) )

  varname = 'lat'
  call check( nf90_inq_varid(ncid,varname,varidlat) )
  call check( nf90_put_var(ncid, varidlat, oraster%lats) )

  call check( nf90_close(ncid) )

  ! deallocate stuff here ...
  if(allocated(sic_ana)) deallocate(sic_ana)
  if(allocated(sit_ana)) deallocate(sit_ana)
  if(allocated(snd_ana)) deallocate(snd_ana)
  if(allocated(sst_incr))deallocate(sst_incr)

  if(allocated(Tw))    deallocate(Tw)
  if(allocated(Sw))    deallocate(Sw)
  if(allocated(Tocean))deallocate(Tocean)

  if(allocated(aicen)) deallocate(aicen)
  if(allocated(vicen)) deallocate(vicen)
  if(allocated(vsnon)) deallocate(vsnon)
  if(allocated(eicen)) deallocate(eicen)
  if(allocated(esnon)) deallocate(esnon)

  if(allocated(DUMVAR)) deallocate(DUMVAR)

contains
     subroutine save_old_state(aicen,vicen,vsnon,eicen,esnon, &
                               nilyr, nslyr, ia, ib, &
                                hi, hs, qins, qsns)

  real*8, dimension(:),   intent(in)  :: aicen, vicen, vsnon
  real*8, dimension(:,:), intent(in)  :: eicen, esnon
  integer*4,              intent(in)  :: nilyr, nslyr, ia, ib
  real*8,                 intent(out) :: hi, hs
  real*8, dimension(:),   intent(out) :: qins, qsns

  integer*4                           :: k

      hi = vicen(ib) / aicen(ia)
      hs = vsnon(ib) / aicen(ia)
      qin_save(:) = 0.
      qsn_save(:) = 0.
      if (vicen(ib) > 0.0) then
           do k=1,nilyr
              qins(k) = eicen(k,ib)         &
                        * real(nilyr,kind=real_kind)             &
                           / vicen(ib)
           enddo
      endif
      if (vsnon(ib) > 0.0) then
          do k=1,nslyr
             qsns(k) = esnon(k,ib)         &
                       * real(nslyr,kind=real_kind)             &
                          / vsnon(ib)
          enddo
      endif
      return
   end subroutine save_old_state

    subroutine restore_ice_state(aicen,vicen,vsnon,eicen,esnon, &
                                nilyr, nslyr, ia, ib, sic_incr, &
                                hi, hs, qins, qsns, himin)

  real*8, dimension(:),   intent(inout) :: aicen
  real*8, dimension(:),   intent(out)  :: vicen, vsnon
  real*8, dimension(:,:), intent(out)  :: eicen, esnon
  integer*4,              intent(in)  :: nilyr, nslyr, ia, ib
  real*8,                 intent(in) :: sic_incr
  real*8,                 intent(in) :: hi, hs
  real*8, dimension(:),   intent(in) :: qins, qsns
  real*4, intent(in), optional       :: himin

  integer*4                           :: k

     if(present(himin)) then
      if(hi >= himin) then
         aicen(ia) = aicen(ia) + sic_incr
         vicen(ib) = aicen(ia) * hi
      else
         vicen(ib) = aicen(ia) * hi + sic_incr * hinew
         aicen(ia) = aicen(ia) + sic_incr
      endif
     else
       aicen(ia) = aicen(ia) + sic_incr
       vicen(ib) = aicen(ia) * hi
     endif
       vsnon(ib) = aicen(ia) * hs
       do k=1,nilyr
            eicen(k,ib)                          &
            = qins(k) * vicen(ib)                  &
              / real(nilyr,kind=real_kind)
       enddo
       do k=1,nslyr
            esnon(k,ib)                          &
               = qsns(k) * vsnon(ib)                  &
               / real(nslyr,kind=real_kind)
       enddo

      return
   end subroutine restore_ice_state

end program anaice2rst


text:
      subroutine readfrag
      use derivheader

c  read the frags.out file
c which was created by the fragmentation program
      open(unit=3,file='frags.out',status='old')

      read(3,*)comlin
      read(3,*)
      read(3,*)comlin
      read(3,*)
      read(3,*)comlin
      do m=1,nfrag
c nfrag = the number of fragments (read in readdata)
      read(3,*)nat0(m)
c nat0 = the number of real atoms in the fragment (readdata)
      enddo
      read(3,*)comlin
      do m=1,nfrag
      read(3,*)(nat(m,i),i=1,nat0(m))
c nat = the whole-molecule atom numbers in each fragment (readdata)
      enddo
      read(3,*)comlin
      do m=1,nfrag
      ncap(m)=0
c ncap = the number of caps in each fragment
      enddo
c  read in the identity of the atoms connected to the caps
444   continue
      read(3,*,end=445)k1,m1,n1,fac
c k1 the fragment number
c m1 and n1 are the atoms which define the cap
c fac is the weight for m1 versus n1
c ncap, natcap and fact were allocated in readdata
      ncap(k1)=ncap(k1)+1
      natcap(k1,ncap(k1),1)=m1
      natcap(k1,ncap(k1),2)=n1
      fact(k1,ncap(k1))=fac
      go to 444
445   continue

c check the data in frags.out matches that read in from the
c ab initio fragment output, read in readdata
      do m=1,nfrag
      if(iabs(nat0(m)+ncap(m)-natomfrag(m)).gt.0)then
      write(6,*)' The number of atoms in fragment ',m
      write(6,*)' does not match the numbers in frags.out'
      stop
      endif
      enddo

c  find out how many real atoms are in the whole molecule

c     natom=0
c     do m=1,nfrag
c     do n=1,nat0(m)
c     if(nat(m,n).gt.natom)natom=nat(m,n)
c     enddo
c     enddo
c     close(unit=3)

      return
      end

ext: f

max_stars_repo_name: mickcollins/SMFA

max_issues_repo_name: mickcollins/SMFA

max_forks_repo_name: mickcollins/SMFAPAC


text:
! Copyright (C) 2015-2017, Simone Marsili
! All rights reserved.
! License: BSD 3 clause

module mpi_wrapper
  use kinds
!  use mpi
  implicit none

  include 'mpif.h'

  integer :: iproc
  integer :: nproc

contains

  subroutine mpi_wrapper_initialize(err)
    !   initialize the MPI environment
    implicit none
    ! local variables
    integer, intent(out) :: err

    CALL MPI_Init(err)

    CALL MPI_Comm_size(MPI_COMM_WORLD, nproc, err)

    CALL MPI_Comm_rank(MPI_COMM_WORLD, iproc, err)

  end subroutine mpi_wrapper_initialize

  subroutine mpi_wrapper_finalize(err)
    implicit none
    integer, intent(out) :: err

    call mpi_finalize(err)

  end subroutine mpi_wrapper_finalize

  subroutine mpi_wrapper_barrier(err)
    implicit none
    integer, intent(out) :: err

    call mpi_barrier(MPI_COMM_WORLD,err)

  end subroutine mpi_wrapper_barrier

  subroutine float_allreduce(dim,a)
    integer,    intent(in)    :: dim
    real(kflt), intent(inout) :: a(dim)
    real(kflt) :: b(dim)
    integer    :: err

    select case(kflt)
    case(4)
       call MPI_ALLREDUCE(a,b,dim,MPI_REAL4,MPI_SUM,MPI_COMM_WORLD,err)
    case(8)
       call MPI_ALLREDUCE(a,b,dim,MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,err)
    end select

    a = b

  end subroutine float_allreduce

  subroutine int_allreduce(dim,a)
    integer, intent(in)    :: dim
    integer, intent(inout) :: a(dim)
    integer :: b(dim)
    integer :: err

    call MPI_ALLREDUCE(a,b,dim,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err)

    a = b

  end subroutine int_allreduce

  subroutine float_reduce(dim,a)
    integer, intent(in)       :: dim
    real(kflt), intent(inout) :: a(dim)
    real(kflt) :: b(dim)
    integer    :: err

    b = a

    select case(kflt)
    case(4)
       call MPI_REDUCE(b,a,dim,MPI_REAL4,MPI_SUM,0,MPI_COMM_WORLD,err)
    case(8)
       call MPI_REDUCE(b,a,dim,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,err)
    end select

  end subroutine float_reduce

  subroutine int_reduce(dim,a)
    integer, intent(in)    :: dim
    integer, intent(inout) :: a(dim)
    integer :: b(dim)
    integer :: err

    b = a

    call MPI_REDUCE(b,a,dim,MPI_INTEGER,MPI_SUM,0,MPI_COMM_WORLD,err)

  end subroutine int_reduce

  subroutine float_bcast(dim,a)
    integer, intent(in)       :: dim
    real(kflt), intent(inout) :: a(dim)
    integer :: err

    select case(kflt)
    case(4)
       call MPI_BCAST(a,dim,MPI_REAL4,0,MPI_COMM_WORLD,err)
    case(8)
       call MPI_BCAST(a,dim,MPI_REAL8,0,MPI_COMM_WORLD,err)
    end select

  end subroutine float_bcast

  subroutine int_bcast(dim,a)
    integer, intent(in)    :: dim
    integer, intent(inout) :: a(dim)
    integer :: err

    call MPI_BCAST(a,dim,MPI_INTEGER,0,MPI_COMM_WORLD,err)

  end subroutine int_bcast

end module mpi_wrapper

max_stars_repo_name: simomarsili/elss

max_issues_repo_name: simomarsili/elss

max_forks_repo_name: simomarsili/elss


text:
#include "pixel.h"

ps.1.1

;----------------------------
#beginfragment tex_1
tex t1
#endfragment

#beginfragment tex_0
tex t0
#endfragment

#beginfragment tex_2
tex t2
#endfragment

#beginfragment tex_3
tex t3
#endfragment

//----------------------------

#beginfragment tex_1_bem
texbem t1, t0
#endfragment

#beginfragment tex_2_bem
texbem t2, t1
#endfragment

#beginfragment tex_3_bem
texbem t3, t2
#endfragment

//----------------------------

#beginfragment tex_1_beml
texbeml t1, t0
#endfragment

#beginfragment tex_2_beml
texbeml t2, t1
#endfragment

#beginfragment tex_3_beml
texbeml t3, t2
#endfragment

;----------------------------
#beginfragment texkill_0
texkill t0
#endfragment

;----------------------------
#beginfragment texkill_1
texkill t1
#endfragment

;----------------------------
#beginfragment texkill_2
texkill t2
#endfragment

;----------------------------
#beginfragment texkill_3
texkill t3
#endfragment

;----------------------------
; Simple copy of diffuse color.
#beginfragment v0_copy
mov r0, v0
#endfragment

;----------------------------
; Simple copy of diffuse color.
#beginfragment v0_mul2x
mov_x2 r0, v0
#endfragment

;----------------------------
; Simple copy of texture.
#beginfragment t0_copy
mov r0, t0
#endfragment

;----------------------------
; Copy of texture.
#beginfragment t0_copy_inv
mov r0, 1-t0
#endfragment

;----------------------------
; Simple copy of texture.
#beginfragment t1_copy
mov r0, t1
#endfragment

;----------------------------
; Modulation in stage 0.
#beginfragment mod_t0_v0
mul r0, v0, t0
#endfragment

;----------------------------
; Modulation with t1
#beginfragment mod_r0_t1
mul r0, r0, t1
#endfragment

;----------------------------
; Modulation in stage 0.
#beginfragment modx2_t0_v0
mul_x2 r0.rgb, v0, t0
+mul r0.a, v0.a, t0.a

#endfragment

;----------------------------
; Simple copy of constant color
#beginfragment color_copy
mov r0, CP_COLOR
#endfragment

;----------------------------
; Simple copy of constant color
#beginfragment mod_color_v0
mul r0, CP_COLOR, v0
#endfragment

;----------------------------
; Copy of black color
#beginfragment copy_black
mov r0.rgb, CP_ZERO
mov r0.a, CP_ONE.a
#endfragment

;----------------------------
; Copy of black color with alpha from t0
#beginfragment copy_black_t0a
mov r0.rgb, CP_ZERO
+mov r0.a, t0.a
#endfragment

;----------------------------
#beginfragment mod_t0_constcolor
mul r0, t0, CP_COLOR
#endfragment

;----------------------------
#beginfragment mod_r0_t1
mul r0, r0, t1
#endfragment

;----------------------------
#beginfragment grayscale
dp3 r_gray, r0, CP_FACTOR
lrp r0, CP_FACTOR.a, r0, r_gray
#endfragment

;----------------------------
#beginfragment modx2_r0_t0
mul_x2 r0.rgb, r0, t0
+mul r0.a, r0.a, t0.a
#endfragment

;----------------------------
#beginfragment modx2_r0_t1
mul_x2 r0.rgb, r0, t1
+mul r0.a, r0.a, t1.a
#endfragment

;----------------------------
#beginfragment modx2_r0_t2
mul_x2 r0.rgb, r0, t2
+mul r0.a, r0.a, t2.a
#endfragment

;----------------------------
#beginfragment modx2_r0_t3
mul_x2 r0.rgb, r0, t3
+mul r0.a, r0.a, t3.a
#endfragment

;----------------------------
#beginfragment modx2_t0_t1
mul_x2 r0.rgb, t0, t1
+mul r0.a, t0.a, t1.a
#endfragment

;----------------------------
#beginfragment modx2_r0_r1
mul_x2 r0.rgb, r0, r1
+mul r0.a, r0.a, r1.a
#endfragment

;----------------------------
#beginfragment modx2_t0_r1
mul_x2 r0.rgb, t0, r1
+mul r0.a, t0.a, r1.a
#endfragment

;----------------------------
#beginfragment add_t0_v0
add r0.rgb, t0, v0
mov r0.a, v0.a
#endfragment

;----------------------------
#beginfragment add_t1_v0
add r0, t1, v0
#endfragment

;----------------------------
; Shadow cast fragment.
#beginfragment shadow_cast

mov r0.rgb, CP_COLOR
+mov r0.a, t0.a

#endfragment

;----------------------------
; Shadow receive fragment.
#beginfragment shadow_receive

mov r0.rgb, t0
+mov r0.a, v0.a

#endfragment

;----------------------------
; Copy r0 blue to alpha.
#beginfragment r0_b_2_a

mov r0.a, r0.b

#endfragment

;----------------------------

#beginfragment blend_by_alpha
                              //modulate t1 and t2 using t0.a
mul r1, t1, 1-t0.a
mad r1, t2, t0.a, r1

#endfragment

//----------------------------
#beginfragment night_view

dp3 r0, t0, CP_COLOR1                      ; put to grayscale
;mad r1, r0, CP_FACTOR1, CP_FACTOR   ; select green or white, based on brightness
mul r0, 1-r0, 1-r0                  ; increase curve by multiplying 1-inversed values together few times
mul r0, r0, r0
mad r1, 1-r0, CP_FACTOR1, CP_FACTOR ; select green or white, based on brightness
;mul r0, r0, r0
mul r0, r0, r0
mul r0, r0, r0

; multiply by selected green-white color, and add some base color
mad r0, 1-r0, r1, CP_COLOR

#endfragment

;----------------------------
; Test fragment.
#beginfragment ps_test_bump

dp3_sat r1, t1_bx2, v1_bx2
;mul r1, t2, r1.a

;add r0.rgb, r0, r1
;mul r1, r1, t1.a
;mad r0, r1, t2, r0
;add_d2 r1, r1, CP_ONE

lrp r1, t1.a, r1, CP_ONE

mad r1, r1, t0, r0

                              //environment
;mul_x2 r1, r1, t2

mov r0, r1

;mov r0, t1.a

#endfragment

;----------------------------
; Test fragment.
#beginfragment ps_test

nop

#endfragment


text:
    !	Copyright 2011 Johns Hopkins University
    !
    !  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.

    program TurbTest

    use, intrinsic :: iso_c_binding, only: c_ptr, c_int, c_float, c_f_pointer
    implicit none

    integer, parameter :: RP=4
    character(*), parameter :: authtoken = 'edu.jhu.pha.turbulence.testing-201406' // CHAR(0)
    character(1) :: field

    integer, parameter :: time_step=1, x_start=1, y_start=10, z_start=1, x_end=2, y_end=11, z_end=2
    integer, parameter :: x_step=1, y_step=1, z_step=1, filter_width=1
    integer :: size, i
    real(RP), allocatable :: result(:)

    ! Declare the return type of the turblib functions as integer.
    ! This is required for custom error handling (see the README).
    integer :: getcutout

    ! return code
    integer :: rc

    ! Formatting rules
    character(*), parameter :: rawformat1='(i4,1(a,f10.6))'

    !
    ! Intialize the gSOAP runtime.
    ! This is required before any WebService routines are called.
    !
    CALL soapinit()

    ! Enable exit on error.  See README for details.
    CALL turblibSetExitOnError(1)

    print *, ".........getAnyCutoutWeb........."

    print *, "......isotropic1024coarse u......";
    field='u'
    if (field=='u') then
        size = (x_end-x_start+1) * (y_end-y_start+1) * (z_end-z_start+1)*3;
    else if (field=='p') then
        size = (x_end-x_start+1) * (y_end-y_start+1) * (z_end-z_start+1);
    end if
    allocate(result(size))
    print *, field=='u', size

    rc = getcutout (authtoken, "isotropic1024coarse", field, time_step, x_start, y_start, z_start, &
        x_end, y_end, z_end, x_step, y_step, z_step, filter_width, result)
    do i = 1, size
        write(*,rawformat1) i, ': ', result(i)
    end do
    deallocate(result)

    print *, ".........transition_bl p.........";
    field='p'
    if (field=='u') then
        size = (x_end-x_start+1) * (y_end-y_start+1) * (z_end-z_start+1)*3;
    else if (field=='p') then
        size = (x_end-x_start+1) * (y_end-y_start+1) * (z_end-z_start+1);
    end if
    allocate(result(size))
    print *, field=='p', size

    rc = getcutout (authtoken, "transition_bl", field, time_step, x_start, y_start, z_start, &
        x_end, y_end, z_end, x_step, y_step, z_step, filter_width, result)
    do i = 1, size
        write(*,rawformat1) i, ': ', result(i)
    end do
    deallocate(result)
    !
    ! Destroy the gSOAP runtime.
    ! No more WebService routines may be called.
    !
    CALL soapdestroy()

    end program TurbTest

max_stars_repo_name: cmhyett/turblib

max_issues_repo_name: cmhyett/turblib

max_forks_repo_name: cmhyett/turblib


text:
      subroutine mdm
     *     (vk,tail, v,l, last,next, mark)
c***********************************************************************
c  mdm -- form element from uneliminated neighbors of vk
c***********************************************************************
      integer  vk, tail,  v(*), l(*),   last(*), next(*),   mark(*),
     *   tag, s,ls,vs,es, b,lb,vb, blp,blpmax
      equivalence  (vs, es)
c
c----initialize tag and list of uneliminated neighbors
      tag = mark(vk)
      tail = vk
c
c----for each vertex/element vs/es in element list of vk
      ls = l(vk)
   1  s = ls
      if (s.eq.0)  go to 5
        ls = l(s)
        vs = v(s)
        if (next(vs).lt.0)  go to 2
c
c------if vs is uneliminated vertex, then tag and append to list of
c------uneliminated neighbors
          mark(vs) = tag
          l(tail) = s
          tail = s
          go to 4
c
c------if es is active element, then ...
c--------for each vertex vb in boundary list of element es
   2      lb = l(es)
          blpmax = last(es)
          do 3 blp=1,blpmax
            b = lb
            lb = l(b)
            vb = v(b)
c
c----------if vb is untagged vertex, then tag and append to list of
c----------uneliminated neighbors
            if (mark(vb).ge.tag)  go to 3
              mark(vb) = tag
              l(tail) = b
              tail = b
   3        continue
c
c--------mark es inactive
          mark(es) = tag
c
   4    go to 1
c
c----terminate list of uneliminated neighbors
   5  l(tail) = 0
c
      return
      end

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


MODULE traqsr
  USE oce
  USE phycst
  USE dom_oce
  USE sbc_oce
  USE trc_oce
  USE trd_oce
  USE trdtra
  USE in_out_manager
  USE prtctl
  USE iom
  USE fldread
  USE restart
  USE lib_mpp
  USE lbclnk
  USE timing
  IMPLICIT NONE
  PRIVATE
  PUBLIC :: tra_qsr
  PUBLIC :: tra_qsr_init
  LOGICAL, PUBLIC :: ln_traqsr
  LOGICAL, PUBLIC :: ln_qsr_rgb
  LOGICAL, PUBLIC :: ln_qsr_2bd
  LOGICAL, PUBLIC :: ln_qsr_bio
  INTEGER, PUBLIC :: nn_chldta
  REAL(KIND = wp), PUBLIC :: rn_abs
  REAL(KIND = wp), PUBLIC :: rn_si0
  REAL(KIND = wp), PUBLIC :: rn_si1
  INTEGER, PUBLIC :: nksr
  INTEGER, PARAMETER :: np_RGB = 1
  INTEGER, PARAMETER :: np_RGBc = 2
  INTEGER, PARAMETER :: np_2BD = 3
  INTEGER, PARAMETER :: np_BIO = 4
  INTEGER :: nqsr
  REAL(KIND = wp) :: xsi0r
  REAL(KIND = wp) :: xsi1r
  REAL(KIND = wp), DIMENSION(3, 61) :: rkrgb
  TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_chl
  CONTAINS
  SUBROUTINE tra_qsr(kt)
    INTEGER, INTENT(IN) :: kt
    INTEGER :: ji, jj, jk
    INTEGER :: irgb
    REAL(KIND = wp) :: zchl, zcoef, z1_2
    REAL(KIND = wp) :: zc0, zc1, zc2, zc3
    REAL(KIND = wp) :: zzc0, zzc1, zzc2, zzc3
    REAL(KIND = wp) :: zz0, zz1
    REAL(KIND = wp) :: zCb, zCmax, zze, zpsi, zpsimax, zdelpsi, zCtot, zCze
    REAL(KIND = wp) :: zlogc, zlogc2, zlogc3
    REAL(KIND = wp), ALLOCATABLE, DIMENSION(:, :) :: zekb, zekg, zekr
    REAL(KIND = wp), ALLOCATABLE, DIMENSION(:, :, :) :: ze0, ze1, ze2, ze3, zea, ztrdt
    REAL(KIND = wp), ALLOCATABLE, DIMENSION(:, :, :) :: zetot, zchl3d
    IF (ln_timing) CALL timing_start('tra_qsr')
    IF (kt == nit000) THEN
      IF (lwp) WRITE(numout, FMT = *)
      IF (lwp) WRITE(numout, FMT = *) 'tra_qsr : penetration of the surface solar radiation'
      IF (lwp) WRITE(numout, FMT = *) '~~~~~~~'
    END IF
    IF (l_trdtra) THEN
      ALLOCATE(ztrdt(jpi, jpj, jpk))
      !$ACC KERNELS
      ztrdt(:, :, :) = tsa(:, :, :, jp_tem)
      !$ACC END KERNELS
    END IF
    IF (kt == nit000) THEN
      IF (ln_rstart .AND. iom_varid(numror, 'qsr_hc_b', ldstop = .FALSE.) > 0) THEN
        IF (lwp) WRITE(numout, FMT = *) '          nit000-1 qsr tracer content forcing field read in the restart file'
        z1_2 = 0.5_wp
        CALL iom_get(numror, jpdom_autoglo, 'qsr_hc_b', qsr_hc_b, ldxios = lrxios)
      ELSE
        !$ACC KERNELS
        z1_2 = 1._wp
        qsr_hc_b(:, :, :) = 0._wp
        !$ACC END KERNELS
      END IF
    ELSE
      !$ACC KERNELS
      z1_2 = 0.5_wp
      qsr_hc_b(:, :, :) = qsr_hc(:, :, :)
      !$ACC END KERNELS
    END IF
    SELECT CASE (nqsr)
    CASE (np_BIO)
      !$ACC KERNELS
      DO jk = 1, nksr
        qsr_hc(:, :, jk) = r1_rau0_rcp * (etot3(:, :, jk) - etot3(:, :, jk + 1))
      END DO
      !$ACC END KERNELS
    CASE (np_RGB, np_RGBc)
      ALLOCATE(zekb(jpi, jpj), zekg(jpi, jpj), zekr(jpi, jpj), ze0(jpi, jpj, jpk), ze1(jpi, jpj, jpk), ze2(jpi, jpj, jpk), ze3(jpi, jpj, jpk), zea(jpi, jpj, jpk), zchl3d(jpi, jpj, jpk))
      IF (nqsr == np_RGBc) THEN
        CALL fld_read(kt, 1, sf_chl)
        DO jk = 1, nksr + 1
          DO jj = 2, jpjm1
            DO ji = 2, jpim1
              zchl = sf_chl(1) % fnow(ji, jj, 1)
              zCtot = 40.6 * zchl ** 0.459
              zze = 568.2 * zCtot ** (- 0.746)
              IF (zze > 102.) zze = 200.0 * zCtot ** (- 0.293)
              zpsi = gdepw_n(ji, jj, jk) / zze
              zlogc = LOG(zchl)
              zlogc2 = zlogc * zlogc
              zlogc3 = zlogc * zlogc * zlogc
              zCb = 0.768 + 0.087 * zlogc - 0.179 * zlogc2 - 0.025 * zlogc3
              zCmax = 0.299 - 0.289 * zlogc + 0.579 * zlogc2
              zpsimax = 0.6 - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3
              zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2
              zCze = 1.12 * (zchl) ** 0.803
              zchl3d(ji, jj, jk) = zCze * (zCb + zCmax * EXP(- ((zpsi - zpsimax) / zdelpsi) ** 2))
            END DO
          END DO
        END DO
      ELSE
        !$ACC KERNELS
        DO jk = 1, nksr + 1
          zchl3d(:, :, jk) = 0.05
        END DO
        !$ACC END KERNELS
      END IF
      !$ACC KERNELS
      zcoef = (1. - rn_abs) / 3._wp
      DO jj = 2, jpjm1
        DO ji = 2, jpim1
          ze0(ji, jj, 1) = rn_abs * qsr(ji, jj)
          ze1(ji, jj, 1) = zcoef * qsr(ji, jj)
          ze2(ji, jj, 1) = zcoef * qsr(ji, jj)
          ze3(ji, jj, 1) = zcoef * qsr(ji, jj)
          zea(ji, jj, 1) = qsr(ji, jj)
        END DO
      END DO
      !$ACC END KERNELS
      DO jk = 2, nksr + 1
        DO jj = 2, jpjm1
          DO ji = 2, jpim1
            zchl = MIN(10., MAX(0.03, zchl3d(ji, jj, jk)))
            irgb = NINT(41 + 20. * LOG10(zchl) + 1.E-15)
            zekb(ji, jj) = rkrgb(1, irgb)
            zekg(ji, jj) = rkrgb(2, irgb)
            zekr(ji, jj) = rkrgb(3, irgb)
          END DO
        END DO
        !$ACC KERNELS
        DO jj = 2, jpjm1
          DO ji = 2, jpim1
            zc0 = ze0(ji, jj, jk - 1) * EXP(- e3t_n(ji, jj, jk - 1) * xsi0r)
            zc1 = ze1(ji, jj, jk - 1) * EXP(- e3t_n(ji, jj, jk - 1) * zekb(ji, jj))
            zc2 = ze2(ji, jj, jk - 1) * EXP(- e3t_n(ji, jj, jk - 1) * zekg(ji, jj))
            zc3 = ze3(ji, jj, jk - 1) * EXP(- e3t_n(ji, jj, jk - 1) * zekr(ji, jj))
            ze0(ji, jj, jk) = zc0
            ze1(ji, jj, jk) = zc1
            ze2(ji, jj, jk) = zc2
            ze3(ji, jj, jk) = zc3
            zea(ji, jj, jk) = (zc0 + zc1 + zc2 + zc3) * wmask(ji, jj, jk)
          END DO
        END DO
        !$ACC END KERNELS
      END DO
      !$ACC KERNELS
      DO jk = 1, nksr
        DO jj = 2, jpjm1
          DO ji = 2, jpim1
            qsr_hc(ji, jj, jk) = r1_rau0_rcp * (zea(ji, jj, jk) - zea(ji, jj, jk + 1))
          END DO
        END DO
      END DO
      !$ACC END KERNELS
      DEALLOCATE(zekb, zekg, zekr, ze0, ze1, ze2, ze3, zea, zchl3d)
    CASE (np_2BD)
      !$ACC KERNELS
      zz0 = rn_abs * r1_rau0_rcp
      zz1 = (1. - rn_abs) * r1_rau0_rcp
      DO jk = 1, nksr
        DO jj = 2, jpjm1
          DO ji = 2, jpim1
            zc0 = zz0 * EXP(- gdepw_n(ji, jj, jk) * xsi0r) + zz1 * EXP(- gdepw_n(ji, jj, jk) * xsi1r)
            zc1 = zz0 * EXP(- gdepw_n(ji, jj, jk + 1) * xsi0r) + zz1 * EXP(- gdepw_n(ji, jj, jk + 1) * xsi1r)
            qsr_hc(ji, jj, jk) = qsr(ji, jj) * (zc0 * wmask(ji, jj, jk) - zc1 * wmask(ji, jj, jk + 1))
          END DO
        END DO
      END DO
      !$ACC END KERNELS
    END SELECT
    !$ACC KERNELS
    DO jk = 1, nksr
      DO jj = 2, jpjm1
        DO ji = 2, jpim1
          tsa(ji, jj, jk, jp_tem) = tsa(ji, jj, jk, jp_tem) + z1_2 * (qsr_hc_b(ji, jj, jk) + qsr_hc(ji, jj, jk)) / e3t_n(ji, jj, jk)
        END DO
      END DO
    END DO
    DO jj = 2, jpjm1
      DO ji = 2, jpim1
        IF (qsr(ji, jj) /= 0._wp) THEN
          fraqsr_1lev(ji, jj) = qsr_hc(ji, jj, 1) / (r1_rau0_rcp * qsr(ji, jj))
        ELSE
          fraqsr_1lev(ji, jj) = 1._wp
        END IF
      END DO
    END DO
    !$ACC END KERNELS
    CALL lbc_lnk(fraqsr_1lev(:, :), 'T', 1._wp)
    IF (iom_use('qsr3d')) THEN
      ALLOCATE(zetot(jpi, jpj, jpk))
      !$ACC KERNELS
      zetot(:, :, nksr + 1 : jpk) = 0._wp
      DO jk = nksr, 1, - 1
        zetot(:, :, jk) = zetot(:, :, jk + 1) + qsr_hc(:, :, jk) * rau0_rcp
      END DO
      !$ACC END KERNELS
      CALL iom_put('qsr3d', zetot)
      DEALLOCATE(zetot)
    END IF
    IF (lrst_oce) THEN
      IF (lwxios) CALL iom_swap(cwxios_context)
      CALL iom_rstput(kt, nitrst, numrow, 'qsr_hc_b', qsr_hc, ldxios = lwxios)
      CALL iom_rstput(kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev, ldxios = lwxios)
      IF (lwxios) CALL iom_swap(cxios_context)
    END IF
    IF (l_trdtra) THEN
      !$ACC KERNELS
      ztrdt(:, :, :) = tsa(:, :, :, jp_tem) - ztrdt(:, :, :)
      !$ACC END KERNELS
      CALL trd_tra(kt, 'TRA', jp_tem, jptra_qsr, ztrdt)
      DEALLOCATE(ztrdt)
    END IF
    IF (ln_ctl) CALL prt_ctl(tab3d_1 = tsa(:, :, :, jp_tem), clinfo1 = ' qsr  - Ta: ', mask1 = tmask, clinfo3 = 'tra-ta')
    IF (ln_timing) CALL timing_stop('tra_qsr')
  END SUBROUTINE tra_qsr
  SUBROUTINE tra_qsr_init
    INTEGER :: ji, jj, jk
    INTEGER :: ios, irgb, ierror, ioptio
    REAL(KIND = wp) :: zz0, zc0, zc1, zcoef
    REAL(KIND = wp) :: zz1, zc2, zc3, zchl
    CHARACTER(LEN = 100) :: cn_dir
    TYPE(FLD_N) :: sn_chl
    NAMELIST /namtra_qsr/ sn_chl, cn_dir, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, nn_chldta, rn_abs, rn_si0, rn_si1
    REWIND(UNIT = numnam_ref)
    READ(numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901)
901 IF (ios /= 0) CALL ctl_nam(ios, 'namtra_qsr in reference namelist', lwp)
    REWIND(UNIT = numnam_cfg)
    READ(numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902)
902 IF (ios > 0) CALL ctl_nam(ios, 'namtra_qsr in configuration namelist', lwp)
    IF (lwm) WRITE(numond, namtra_qsr)
    IF (lwp) THEN
      WRITE(numout, FMT = *)
      WRITE(numout, FMT = *) 'tra_qsr_init : penetration of the surface solar radiation'
      WRITE(numout, FMT = *) '~~~~~~~~~~~~'
      WRITE(numout, FMT = *) '   Namelist namtra_qsr : set the parameter of penetration'
      WRITE(numout, FMT = *) '      RGB (Red-Green-Blue) light penetration       ln_qsr_rgb = ', ln_qsr_rgb
      WRITE(numout, FMT = *) '      2 band               light penetration       ln_qsr_2bd = ', ln_qsr_2bd
      WRITE(numout, FMT = *) '      bio-model            light penetration       ln_qsr_bio = ', ln_qsr_bio
      WRITE(numout, FMT = *) '      RGB : Chl data (=1) or cst value (=0)        nn_chldta  = ', nn_chldta
      WRITE(numout, FMT = *) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs     = ', rn_abs
      WRITE(numout, FMT = *) '      RGB & 2 bands: shortess depth of extinction  rn_si0     = ', rn_si0
      WRITE(numout, FMT = *) '      2 bands: longest depth of extinction         rn_si1     = ', rn_si1
      WRITE(numout, FMT = *)
    END IF
    ioptio = 0
    IF (ln_qsr_rgb) ioptio = ioptio + 1
    IF (ln_qsr_2bd) ioptio = ioptio + 1
    IF (ln_qsr_bio) ioptio = ioptio + 1
    IF (ioptio /= 1) CALL ctl_stop('Choose ONE type of light penetration in namelist namtra_qsr', ' 2 bands, 3 RGB bands or bio-model light penetration')
    IF (ln_qsr_rgb .AND. nn_chldta == 0) nqsr = np_RGB
    IF (ln_qsr_rgb .AND. nn_chldta == 1) nqsr = np_RGBc
    IF (ln_qsr_2bd) nqsr = np_2BD
    IF (ln_qsr_bio) nqsr = np_BIO
    xsi0r = 1._wp / rn_si0
    xsi1r = 1._wp / rn_si1
    SELECT CASE (nqsr)
    CASE (np_RGB, np_RGBc)
      IF (lwp) WRITE(numout, FMT = *) '   ==>>>   R-G-B   light penetration '
      CALL trc_oce_rgb(rkrgb)
      nksr = trc_oce_ext_lev(r_si2, 33._wp)
      IF (lwp) WRITE(numout, FMT = *) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr + 1), ' m'
      IF (nqsr == np_RGBc) THEN
        IF (lwp) WRITE(numout, FMT = *) '   ==>>>   Chlorophyll read in a file'
        ALLOCATE(sf_chl(1), STAT = ierror)
        IF (ierror > 0) THEN
          CALL ctl_stop('tra_qsr_init: unable to allocate sf_chl structure')
          RETURN
        END IF
        ALLOCATE(sf_chl(1) % fnow(jpi, jpj, 1))
        IF (sn_chl % ln_tint) ALLOCATE(sf_chl(1) % fdta(jpi, jpj, 1, 2))
        CALL fld_fill(sf_chl, (/sn_chl/), cn_dir, 'tra_qsr_init', 'Solar penetration function of read chlorophyll', 'namtra_qsr', no_print)
      END IF
      IF (nqsr == np_RGB) THEN
        IF (lwp) WRITE(numout, FMT = *) '   ==>>>   Constant Chlorophyll concentration = 0.05'
      END IF
    CASE (np_2BD)
      IF (lwp) WRITE(numout, FMT = *) '   ==>>>   2 bands light penetration'
      nksr = trc_oce_ext_lev(rn_si1, 100._wp)
      IF (lwp) WRITE(numout, FMT = *) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr + 1), ' m'
    CASE (np_BIO)
      IF (lwp) WRITE(numout, FMT = *) '   ==>>>   bio-model light penetration'
      IF (.NOT. lk_top) CALL ctl_stop('No bio model : ln_qsr_bio = true impossible ')
    END SELECT
    !$ACC KERNELS
    qsr_hc(:, :, :) = 0._wp
    !$ACC END KERNELS
    IF (iom_varid(numror, 'fraqsr_1lev', ldstop = .FALSE.) > 0) THEN
      CALL iom_get(numror, jpdom_autoglo, 'fraqsr_1lev', fraqsr_1lev, ldxios = lrxios)
    ELSE
      !$ACC KERNELS
      fraqsr_1lev(:, :) = 1._wp
      !$ACC END KERNELS
    END IF
    IF (lwxios) THEN
      CALL iom_set_rstw_var_active('qsr_hc_b')
      CALL iom_set_rstw_var_active('fraqsr_1lev')
    END IF
  END SUBROUTINE tra_qsr_init
END MODULE traqsr