!
!  Copyright 2019-2020 SALMON developers
!
!  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.
!
module structure_opt_sub
  implicit none
contains

  !==============================================================initilize
  subroutine  initialization_opt(Miopt,opt,system,flag_opt_conv,nopt_max)
    use structures, only: s_opt,s_dft_system
    use parallelization, only: nproc_id_global
    use communication, only: comm_is_root
    use salmon_global, only: natom,nopt,yn_restart
    use write_sub, only: write_xyz
    use checkpoint_restart_sub, only: restart_opt
    implicit none
    type(s_opt) :: opt
    type(s_dft_system) :: system
    integer :: nopt_max,Miopt,NA3
    logical :: flag_opt_conv
    character(100) :: comment_line

    flag_opt_conv = .false.
    nopt_max = nopt

    NA3 = 3*natom
    allocate( opt%a_dRion(NA3), opt%dFion(NA3) )
    allocate( opt%Hess_mat(NA3,NA3), opt%Hess_mat_last(NA3,NA3) )

    if(yn_restart == 'y') then
       call restart_opt(Miopt,opt)
    else
       opt%a_dRion(:) = 0d0
       opt%dFion(:)   = 0d0
       opt%Hess_mat(:,:)      = 0d0
       opt%Hess_mat_last(:,:) = 0d0
    endif

    if(comm_is_root(nproc_id_global))then
      write(*,*) "===== Grand State Optimization Start ====="
      write(*,*) "       (Quasi-Newton method using Force only)       "
    end if

    write(comment_line,10) 0
    call write_xyz(comment_line,"new","r  ",system)
10  format("#opt iteration step=",i5)

  end subroutine  initialization_opt

  !======================================================convergence check
  subroutine structure_opt_check(iopt,flag_opt_conv,Force)
    use structures
    use salmon_global, only: natom,convrg_opt_fmax,unit_system,flag_opt_atom
    use inputoutput, only: au_length_aa, au_energy_ev
    use parallelization, only: nproc_id_global,nproc_group_global
    use communication, only: comm_is_root,comm_bcast
    implicit none
    integer,intent(in) :: iopt
    real(8),intent(in) :: Force(3,natom)
    logical,intent(inout) :: flag_opt_conv
    integer :: iatom,iatom_count
    real(8) :: fabs,fmax,fave

    fmax =0d0
    fave =0d0
    iatom_count=0
    do iatom=1,natom
      if(flag_opt_atom(iatom)=='y') then
        iatom_count = iatom_count+1
        fabs = sum(Force(:,iatom)**2d0)
        fave = fave + fabs
        if(fabs>=fmax) fmax = fabs
      end if
    enddo

    select case(unit_system)
    case('au','a.u.')
      fmax = sqrt(fmax)
      fave = sqrt(fave/iatom_count)
    case('A_eV_fs')
      fmax = sqrt(fmax)*au_energy_ev/au_length_aa
      fave = sqrt(fave/iatom_count)*au_energy_ev/au_length_aa
    end select

    if(comm_is_root(nproc_id_global))then
      write(*,*) " Max-force=",fmax, "  Mean-force=",fave
      write(*,*) "==================================================="
      write(*,*) "Quasi-Newton Optimization Step = ", iopt
      if( fmax <= convrg_opt_fmax ) flag_opt_conv=.true.
    end if
    call comm_bcast(flag_opt_conv,nproc_group_global)

  end subroutine structure_opt_check

  !===========================================================optimization
  subroutine structure_opt(opt,iopt,system)
    use structures, only: s_dft_system, s_opt
    use salmon_global, only: natom, flag_opt_atom, max_step_len_adjust
    use communication, only: comm_bcast
    implicit none
    type(s_dft_system),intent(inout) :: system
    type(s_opt) :: opt
    integer,intent(in) :: iopt
    !theta_opt=0.0d0:DFP,theta_opt=1.0d0:BFGS in Quasi_Newton method
    real(8), parameter :: theta_opt=1.0d0  !alpha=1.0d0 -- this is now from input
    real(8) :: alpha
    integer :: ii,ij,jj,icount,iatom, NA3,ixyz
    real(8) :: const1,const2,rtmp
    real(8) :: dRion(3,natom), dRabs(natom), dRabs_max
    real(8) :: force_1d(3*natom),dRion_1d(3*natom),optmat_1d(3*natom)
    real(8) :: optmat1_2d(3*natom,3*natom),optmat2_2d(3*natom,3*natom),optmat3_2d(3*natom,3*natom)

    alpha = 1.0d0

    NA3 = 3*natom

    icount=1
    do iatom=1,natom
    do ixyz=1,3
       force_1d(icount) = system%Force(ixyz,iatom)
       icount = icount+1
    end do
    end do

    if(iopt==1)then
      !update opt%Hess_mat
      do ii=1,NA3
      do ij=1,NA3
         if(ii==ij)then
            opt%Hess_mat(ii,ij) = 1d0
         else
            opt%Hess_mat(ii,ij) = 0d0
         end if
         opt%Hess_mat_last(ii,ij) = opt%Hess_mat(ii,ij)
      end do
      end do
    else
      !update opt%dFion
      opt%dFion=-(force_1d-opt%dFion)
      !prepare const and matrix
      call dgemm('n','n',1,1,NA3,1.0d0,opt%a_dRion,1,opt%dFion,NA3,0d0,const1,1)
      call dgemm('n','n',NA3,1,NA3,1d0,opt%Hess_mat,NA3,opt%dFion,NA3,0d0,optmat_1d,NA3)
      call dgemm('n','n',1,1,NA3,1d0,opt%dFion,1,optmat_1d,NA3,0d0,const2,1)
      call dgemm('n','n',NA3,NA3,1,1d0,opt%a_dRion,NA3,opt%a_dRion,1,0d0,optmat1_2d,NA3)
      !update opt%Hess_mat
      rtmp = (const1+theta_opt*const2)/(const1**2d0)
      !$omp parallel do collapse(2) private(ii,jj)
      do ii=1,NA3
      do jj=1,NA3
         opt%Hess_mat(ii,jj) = opt%Hess_mat_last(ii,jj) + rtmp * optmat1_2d(ii,jj)
      enddo
      enddo
      !$omp end parallel do
      if(theta_opt==0.0d0)then
        !theta_opt=0.0d0:DFP
        call dgemm('n','n',NA3,NA3,1,1d0,optmat_1d,NA3,optmat_1d,1,0d0,optmat2_2d,NA3)
        !$omp parallel do collapse(2) private(ii,jj)
        do ii=1,NA3
        do jj=1,NA3
           opt%Hess_mat(ii,jj) = opt%Hess_mat(ii,jj)-(1d0/const2)*optmat2_2d(ii,jj)
        enddo
        enddo
        !$omp end parallel do
      elseif(theta_opt==1.0d0)then
        !theta_opt=1.0d0:BFGS
        call dgemm('n','n',NA3,NA3,1,1d0,optmat_1d,NA3,opt%a_dRion,1,0d0,optmat2_2d,NA3)
        call dgemm('n','n',NA3,NA3,1,1d0,opt%a_dRion,NA3,optmat_1d,1,0d0,optmat3_2d,NA3)
        rtmp = theta_opt/const1
        !$omp parallel do collapse(2) private(ii,jj)
        do ii=1,NA3
        do jj=1,NA3
           opt%Hess_mat(ii,jj) = opt%Hess_mat(ii,jj)- rtmp *(optmat2_2d(ii,jj)+optmat3_2d(ii,jj))
        enddo
        enddo
        !$omp end parallel do
      endif
      !update opt%Hess_mat_last
      !$omp parallel do collapse(2) private(ii,jj)
      do ii=1,NA3
      do jj=1,NA3
         opt%Hess_mat_last(ii,jj) = opt%Hess_mat(ii,jj)
      enddo
      enddo
      !$omp end parallel do
    end if

    !update dRion_1d and dRion
    dRion_1d(:) = 0d0
    call dgemm('n','n',NA3,1,NA3,1d0,opt%Hess_mat,NA3,force_1d,NA3,0d0,dRion_1d,NA3)
    !$omp parallel do collapse(2) private(iatom,ixyz)
    do iatom=1,natom
    do ixyz=1,3
       dRion(ixyz,iatom) = dRion_1d(ixyz+3*(iatom-1))
!      dRion(1:3,iatom) = dRion_1d((1+3*(iatom-1)):(3+3*(iatom-1)))
    end do
    end do
    !$omp end parallel do

    !adjust alpha if required from input
    if(max_step_len_adjust .gt. 0d0) then
       !$omp parallel do private(iatom)
       do iatom=1,natom
          if(flag_opt_atom(iatom)=='y') then
               dRabs(iatom) = sqrt( sum(dRion(:,iatom)**2d0) )
          else
               dRabs(iatom) = 0d0
          endif
       enddo
       !$omp end parallel do
       dRabs_max = maxval(dRabs(:))
       alpha = max_step_len_adjust / dRabs_max
    endif

    !update opt%a_dRion,opt%dFion
    !$omp parallel do private(ii)
    do ii=1,NA3
       opt%a_dRion(ii) = alpha * dRion_1d(ii)
       opt%dFion(ii)   = force_1d(ii)
    enddo
    !$omp end parallel do

    !update Rion
    !$omp parallel do private(iatom)
    do iatom=1,natom
       if(flag_opt_atom(iatom)=='y') &
          system%Rion(1:3,iatom) = system%Rion(1:3,iatom) +alpha*dRion(1:3,iatom)
    end do
    !$omp end parallel do

  end subroutine structure_opt

  !===============================================================finilize
  subroutine structure_opt_fin(opt)
    use structures, only: s_opt
    use parallelization, only: nproc_id_global
    use communication, only: comm_is_root
    implicit none
    type(s_opt) :: opt
    deallocate(opt%a_dRion,opt%dFion)
    deallocate(opt%Hess_mat,opt%Hess_mat_last)
    if(comm_is_root(nproc_id_global)) write(*,*) "Optimization Converged"
  end subroutine structure_opt_fin

end module structure_opt_sub


text:
! KGEN-generated Fortran source file
!
! Filename    : mo_lrtm_driver.f90
! Generated at: 2015-02-19 15:30:29
! KGEN version: 0.4.4

    MODULE mo_lrtm_driver
        USE mo_kind, ONLY: wp
        USE mo_physical_constants, ONLY: amw
        USE mo_physical_constants, ONLY: amd
        USE mo_physical_constants, ONLY: grav
        USE mo_rrtm_params, ONLY: nbndlw
        USE mo_rrtm_params, ONLY: ngptlw
        USE mo_radiation_parameters, ONLY: do_gpoint
        USE mo_radiation_parameters, ONLY: i_overlap
        USE mo_radiation_parameters, ONLY: l_do_sep_clear_sky
        USE mo_radiation_parameters, ONLY: rad_undef
        USE mo_lrtm_setup, ONLY: ngb
        USE mo_lrtm_setup, ONLY: ngs
        USE mo_lrtm_setup, ONLY: delwave
        USE mo_lrtm_setup, ONLY: nspa
        USE mo_lrtm_setup, ONLY: nspb
        USE rrlw_planck, ONLY: totplanck
        USE mo_rrtm_coeffs, ONLY: lrtm_coeffs
        USE mo_lrtm_gas_optics, ONLY: gas_optics_lw
        USE mo_lrtm_solver, ONLY: find_secdiff
        USE mo_lrtm_solver, ONLY: lrtm_solver
        USE mo_cld_sampling, ONLY: sample_cld_state
        USE mo_spec_sampling, ONLY: spec_sampling_strategy
        USE mo_spec_sampling, ONLY: get_gpoint_set
        USE mo_taumol03, ONLY: taumol03_lwr,taumol03_upr
        USE mo_taumol04, ONLY: taumol04_lwr,taumol04_upr
        USE rrlw_planck, ONLY: chi_mls
        USE rrlw_kg03, ONLY: selfref
        USE rrlw_kg03, ONLY: forref
        USE rrlw_kg03, ONLY: ka_mn2o
        USE rrlw_kg03, ONLY: absa
        USE rrlw_kg03, ONLY: fracrefa
        USE rrlw_kg03, ONLY: kb_mn2o
        USE rrlw_kg03, ONLY: absb
        USE rrlw_kg03, ONLY: fracrefb
        IMPLICIT NONE
        PRIVATE
        PUBLIC lrtm
        CONTAINS

        ! read subroutines
        !-----------------------------------------------------------------------------
        !>
        !! @brief Prepares information for radiation call
        !!
        !! @remarks: This program is the driver subroutine for the longwave radiative
        !! transfer routine.  This routine is adapted from the AER LW RRTMG_LW model
        !! that itself has been adapted from RRTM_LW for improved efficiency.  Our
        !! routine does the spectral integration externally (the solver is explicitly
        !! called for each g-point, so as to facilitate sampling of g-points
        !! This routine:
        !!    1) calls INATM to read in the atmospheric profile from GCM;
        !!       all layering in RRTMG is ordered from surface to toa.
        !!    2) calls COEFFS to calculate various quantities needed for
        !!       the radiative transfer algorithm.  This routine is called only once for
        !!       any given thermodynamic state, i.e., it does not change if clouds chanege
        !!    3) calls TAUMOL to calculate gaseous optical depths for each
        !!       of the 16 spectral bands, this is updated band by band.
        !!    4) calls SOLVER (for both clear and cloudy profiles) to perform the
        !!       radiative transfer calculation with a maximum-random cloud
        !!       overlap method, or calls RTRN to use random cloud overlap.
        !!    5) passes the necessary fluxes and cooling rates back to GCM
        !!
        !

        SUBROUTINE lrtm(kproma, kbdim, klev, play, psfc, tlay, tlev, tsfc, wkl, wx, coldry, emis, cldfr, taucld, tauaer, rnseeds, &
        strategy, n_gpts_ts, uflx, dflx, uflxc, dflxc)
            INTEGER, intent(in) :: klev
            INTEGER, intent(in) :: kbdim
            INTEGER, intent(in) :: kproma
            !< Maximum block length
            !< Number of horizontal columns
            !< Number of model layers
            REAL(KIND=wp), intent(in) :: wx(:,:,:)
            REAL(KIND=wp), intent(in) :: cldfr(kbdim,klev)
            REAL(KIND=wp), intent(in) :: taucld(kbdim,klev,nbndlw)
            REAL(KIND=wp), intent(in) :: wkl(:,:,:)
            REAL(KIND=wp), intent(in) :: coldry(kbdim,klev)
            REAL(KIND=wp), intent(in) :: play(kbdim,klev)
            REAL(KIND=wp), intent(in) :: tlay(kbdim,klev)
            REAL(KIND=wp), intent(in) :: tauaer(kbdim,klev,nbndlw)
            REAL(KIND=wp), intent(in) :: tlev(kbdim,klev+1)
            REAL(KIND=wp), intent(in) :: tsfc(kbdim)
            REAL(KIND=wp), intent(in) :: psfc(kbdim)
            REAL(KIND=wp), intent(in) :: emis(kbdim,nbndlw)
            !< Layer pressures [hPa, mb] (kbdim,klev)
            !< Surface pressure [hPa, mb] (kbdim)
            !< Layer temperatures [K] (kbdim,klev)
            !< Interface temperatures [K] (kbdim,klev+1)
            !< Surface temperature [K] (kbdim)
            !< Gas volume mixing ratios
            !< CFC type gas volume mixing ratios
            !< Column dry amount
            !< Surface emissivity  (kbdim,nbndlw)
            !< Cloud fraction  (kbdim,klev)
            !< Coud optical depth (kbdim,klev,nbndlw)
            !< Aerosol optical depth (kbdim,klev,nbndlw)
            ! Variables for sampling cloud state and spectral points
            INTEGER, intent(inout) :: rnseeds(:, :) !< Seeds for random number generator (kbdim,:)
            TYPE(spec_sampling_strategy), intent(in) :: strategy
            INTEGER, intent(in   ) :: n_gpts_ts
            REAL(KIND=wp), intent(out) :: uflx (kbdim,0:klev)
            REAL(KIND=wp), intent(out) :: dflx (kbdim,0:klev)
            REAL(KIND=wp), intent(out) :: uflxc(kbdim,0:klev)
            REAL(KIND=wp), intent(out) :: dflxc(kbdim,0:klev)
            !< Tot sky longwave upward   flux [W/m2], (kbdim,0:klev)
            !< Tot sky longwave downward flux [W/m2], (kbdim,0:klev)
            !< Clr sky longwave upward   flux [W/m2], (kbdim,0:klev)
            !< Clr sky longwave downward flux [W/m2], (kbdim,0:klev)
            REAL(KIND=wp) :: taug(klev) !< Properties for one column at a time >! gas optical depth
            REAL(KIND=wp) :: rrpk_taug03(kbdim,klev)
            REAL(KIND=wp) :: rrpk_taug04(kbdim,klev)
            REAL(KIND=wp) :: fracs(kbdim,klev,n_gpts_ts)
            REAL(KIND=wp) :: taut  (kbdim,klev,n_gpts_ts)
            REAL(KIND=wp) :: tautot(kbdim,klev,n_gpts_ts)
            REAL(KIND=wp) :: pwvcm(kbdim)
            REAL(KIND=wp) :: secdiff(kbdim)
            !< Planck fraction per g-point
            !< precipitable water vapor [cm]
            !< diffusivity angle for RT calculation
            !< gaseous + aerosol optical depths for all columns
            !< cloud + gaseous + aerosol optical depths for all columns
            REAL(KIND=wp) :: planklay(kbdim,  klev,nbndlw)
            REAL(KIND=wp) :: planklev(kbdim,0:klev,nbndlw)
            REAL(KIND=wp) :: plankbnd(kbdim,       nbndlw) ! Properties for all bands
            ! Planck function at mid-layer
            ! Planck function at level interfaces
            ! Planck function at surface
            REAL(KIND=wp) :: layplnk(kbdim,  klev)
            REAL(KIND=wp) :: levplnk(kbdim,0:klev)
            REAL(KIND=wp) :: bndplnk(kbdim)
            REAL(KIND=wp) :: srfemis(kbdim) ! Properties for a single set of columns/g-points
            ! Planck function at mid-layer
            ! Planck function at level interfaces
            ! Planck function at surface
            ! Surface emission
            REAL(KIND=wp) :: zgpfd(kbdim,0:klev)
            REAL(KIND=wp) :: zgpfu(kbdim,0:klev)
            REAL(KIND=wp) :: zgpcu(kbdim,0:klev)
            REAL(KIND=wp) :: zgpcd(kbdim,0:klev)
            ! < gpoint clearsky downward flux
            ! < gpoint clearsky downward flux
            ! < gpoint fullsky downward flux
            ! < gpoint fullsky downward flux
            ! -----------------
            ! Variables for gas optics calculations
            INTEGER :: jt1     (kbdim,klev)
            INTEGER :: indfor  (kbdim,klev)
            INTEGER :: indself (kbdim,klev)
            INTEGER :: indminor(kbdim,klev)
            INTEGER :: laytrop (kbdim     )
            INTEGER :: jp      (kbdim,klev)
            INTEGER :: rrpk_jp (klev,kbdim)
            INTEGER :: jt      (kbdim,klev)
            INTEGER :: rrpk_jt (kbdim,0:1,klev)
            !< tropopause layer index
            !< lookup table index
            !< lookup table index
            !< lookup table index
            REAL(KIND=wp) :: wbrodl      (kbdim,klev)
            REAL(KIND=wp) :: selffac     (kbdim,klev)
            REAL(KIND=wp) :: colh2o      (kbdim,klev)
            REAL(KIND=wp) :: colo3       (kbdim,klev)
            REAL(KIND=wp) :: coln2o      (kbdim,klev)
            REAL(KIND=wp) :: colco       (kbdim,klev)
            REAL(KIND=wp) :: selffrac    (kbdim,klev)
            REAL(KIND=wp) :: colch4      (kbdim,klev)
            REAL(KIND=wp) :: colo2       (kbdim,klev)
            REAL(KIND=wp) :: colbrd      (kbdim,klev)
            REAL(KIND=wp) :: minorfrac   (kbdim,klev)
            REAL(KIND=wp) :: scaleminorn2(kbdim,klev)
            REAL(KIND=wp) :: scaleminor  (kbdim,klev)
            REAL(KIND=wp) :: forfac      (kbdim,klev)
            REAL(KIND=wp) :: colco2      (kbdim,klev)
            REAL(KIND=wp) :: forfrac     (kbdim,klev)
            !< column amount (h2o)
            !< column amount (co2)
            !< column amount (o3)
            !< column amount (n2o)
            !< column amount (co)
            !< column amount (ch4)
            !< column amount (o2)
            !< column amount (broadening gases)
            REAL(KIND=wp) :: wx_loc(size(wx, 2), size(wx, 3))
            !< Normalized CFC amounts (molecules/cm^2)
            REAL(KIND=wp) :: fac00(kbdim,klev)
            REAL(KIND=wp) :: fac01(kbdim,klev)
            REAL(KIND=wp) :: fac10(kbdim,klev)
            REAL(KIND=wp) :: fac11(kbdim,klev)
            REAL(KIND=wp) :: rrpk_fac0(kbdim,0:1,klev)
            REAL(KIND=wp) :: rrpk_fac1(kbdim,0:1,klev)
            REAL(KIND=wp) :: rat_n2oco2  (kbdim,klev)
            REAL(KIND=wp) :: rat_o3co2   (kbdim,klev)
            REAL(KIND=wp) :: rat_h2on2o  (kbdim,klev)
            REAL(KIND=wp) :: rat_n2oco2_1(kbdim,klev)
            REAL(KIND=wp) :: rat_h2on2o_1(kbdim,klev)
            REAL(KIND=wp) :: rat_h2oco2_1(kbdim,klev)
            REAL(KIND=wp) :: rat_h2oo3   (kbdim,klev)
            REAL(KIND=wp) :: rat_h2och4  (kbdim,klev)
            REAL(KIND=wp) :: rat_h2oco2  (kbdim,klev)
            REAL(KIND=wp) :: rrpk_rat_h2oco2  (kbdim,0:1,klev)
            REAL(KIND=wp) :: rrpk_rat_o3co2  (kbdim,0:1,klev)
            REAL(KIND=wp) :: rat_h2oo3_1 (kbdim,klev)
            REAL(KIND=wp) :: rat_o3co2_1 (kbdim,klev)
            REAL(KIND=wp) :: rat_h2och4_1(kbdim,klev)
            ! -----------------
            INTEGER :: jl
            INTEGER :: ig
            INTEGER :: jk ! loop indicies
            INTEGER :: igs(kbdim, n_gpts_ts)
            INTEGER :: ibs(kbdim, n_gpts_ts)
            INTEGER :: ib
            INTEGER :: igpt
            ! minimum val for clouds
            ! Variables for sampling strategy
            REAL(KIND=wp) :: gpt_scaling
            REAL(KIND=wp) :: clrsky_scaling(1:kbdim)
            REAL(KIND=wp) :: smp_tau(kbdim, klev, n_gpts_ts)
            LOGICAL :: cldmask(kbdim, klev, n_gpts_ts)
            LOGICAL :: colcldmask(kbdim,       n_gpts_ts) !< cloud mask in each cell
            !< cloud mask for each column
            !
            ! --------------------------------
            !
            ! 1.0 Choose a set of g-points to do consistent with the spectral sampling strategy
            !
            ! --------------------------------
            gpt_scaling = real(ngptlw,kind=wp)/real(n_gpts_ts,kind=wp)
            ! Standalone logic
            IF (do_gpoint == 0) THEN
                igs(1:kproma,1:n_gpts_ts) = get_gpoint_set(kproma, kbdim, strategy, rnseeds)
                ELSE IF (n_gpts_ts == 1) THEN ! Standalone logic
                IF (do_gpoint > ngptlw) RETURN
                igs(:, 1:n_gpts_ts) = do_gpoint
                ELSE
                PRINT *, "Asking for gpoint fluxes for too many gpoints!"
                STOP
            END IF
            ! Save the band nunber associated with each gpoint
            DO jl = 1, kproma
                DO ig = 1, n_gpts_ts
                    ibs(jl, ig) = ngb(igs(jl, ig))
                END DO
            END DO
            !
            ! ---  2.0 Optical properties
            !
            ! ---  2.1 Cloud optical properties.
            ! --------------------------------
            ! Cloud optical depth is only saved for the band associated with this g-point
            !   We sample clouds first because we may want to adjust water vapor based
            !   on presence/absence of clouds
            !
            CALL sample_cld_state(kproma, kbdim, klev, n_gpts_ts, rnseeds(:,:), i_overlap, cldfr(:,:), cldmask(:,:,:))
            !IBM* ASSERT(NODEPS)
            DO ig = 1, n_gpts_ts
                DO jl = 1, kproma
                    smp_tau(jl,:,ig) = merge(taucld(jl,1:klev,ibs(jl,ig)), 0._wp, cldmask(jl,:,ig))
                END DO
            END DO  ! Loop over samples - done with cloud optical depth calculations
            !
            ! Cloud masks for sorting out clear skies - by cell and by column
            !
            IF (.not. l_do_sep_clear_sky) THEN
                !
                ! Are any layers cloudy?
                !
                colcldmask(1:kproma,       1:n_gpts_ts) = any(cldmask(1:kproma,1:klev,1:n_gpts_ts), dim=2)
                !
                ! Clear-sky scaling is gpt_scaling/frac_clr or 0 if all samples are cloudy
                !
                clrsky_scaling(1:kproma) = gpt_scaling *                                                                          &
                         merge(real(n_gpts_ts,kind=wp) /                                        (real(n_gpts_ts - count(&
                colcldmask(1:kproma,:),dim=2),kind=wp)),                                                                          &
                  0._wp,                                                          any(.not. colcldmask(1:kproma,:),dim=2))
            END IF
            !
            ! ---  2.2. Gas optical depth calculations
            !
            ! --------------------------------
            !
            ! 2.2.1  Calculate information needed by the radiative transfer routine
            ! that is specific to this atmosphere, especially some of the
            ! coefficients and indices needed to compute the optical depths
            ! by interpolating data from stored reference atmospheres.
            ! The coefficients are functions of temperature and pressure and remain the same
            ! for all g-point samples.
            ! If gas concentrations, temperatures, or pressures vary with sample (ig)
            !   the coefficients need to be calculated inside the loop over samples
            ! --------------------------------
            !
            ! Broadening gases -- the number of molecules per cm^2 of all gases not specified explicitly
            !   (water is excluded)
            wbrodl(1:kproma,1:klev) = coldry(1:kproma,1:klev) - sum(wkl(1:kproma,2:,1:klev), dim=2)
            CALL lrtm_coeffs(kproma, kbdim, klev, play, tlay, coldry, wkl, wbrodl, laytrop, jp, jt, jt1, colh2o, colco2, colo3, &
            coln2o, colco, colch4, colo2, colbrd, fac00, fac01, fac10, fac11, rat_h2oco2, rat_h2oco2_1, rat_h2oo3, rat_h2oo3_1, &
            rat_h2on2o, rat_h2on2o_1, rat_h2och4, rat_h2och4_1, rat_n2oco2, rat_n2oco2_1, rat_o3co2, rat_o3co2_1, selffac, &
            selffrac, indself, forfac, forfrac, indfor, minorfrac, scaleminor, scaleminorn2, indminor)
            !
            !  2.2.2 Loop over g-points calculating gas optical properties.
            !
            ! --------------------------------
            !IBM* ASSERT(NODEPS)
            rrpk_rat_h2oco2(:,0,:) = rat_h2oco2
            rrpk_rat_h2oco2(:,1,:) = (rat_h2oco2_1)
            rrpk_rat_o3co2(:,0,:) = rat_o3co2
            rrpk_rat_o3co2(:,1,:) = (rat_o3co2_1)
            rrpk_fac0(:,0,:) = fac00
            rrpk_fac0(:,1,:) = fac01
            rrpk_fac1(:,0,:) = fac10
            rrpk_fac1(:,1,:) = fac11
            rrpk_jt(:,0,:) = jt
            rrpk_jt(:,1,:) = jt1
            DO ig = 1, n_gpts_ts
                igpt=igs(1,ig)
                IF(ngb(igpt) == 3) Then
                        jl=kproma
                        call taumol03_lwr(jl,laytrop(1), klev,                          &
                            rrpk_rat_h2oco2(1:jl,:,:), colco2(1:jl,:), colh2o(1:jl,:), coln2o(1:jl,:), coldry(1:jl,:), &
                            rrpk_fac0(1:jl,:,:), rrpk_fac1(1:jl,:,:), minorfrac(1:jl,:), &
                            selffac(1:jl,:),selffrac(1:jl,:),forfac(1:jl,:),forfrac(1:jl,:), &
                            jp(1:jl,:), rrpk_jt(1:jl,:,:), (igpt-ngs(ngb(igpt)-1)), indself(1:jl,:), &
                            indfor(1:jl,:), indminor(1:jl,:), &
                            rrpk_taug03(1:jl,:),fracs(1:jl,:,ig))
                        call taumol03_upr(jl,laytrop(1), klev,                          &
                            rrpk_rat_h2oco2(1:jl,:,:), colco2(1:jl,:), colh2o(1:jl,:), coln2o(1:jl,:), coldry(1:jl,:), &
                            rrpk_fac0(1:jl,:,:), rrpk_fac1(1:jl,:,:), minorfrac(1:jl,:), &
                            forfac(1:jl,:),forfrac(1:jl,:),           &
                            jp(1:jl,:), rrpk_jt(1:jl,:,:), (igpt-ngs(ngb(igpt)-1)), &
                            indfor(1:jl,:), indminor(1:jl,:), &
                            rrpk_taug03(1:jl,:),fracs(1:jl,:,ig))
                ENDIF
                IF(ngb(igpt) == 4) Then
                        jl=kproma
                        call taumol04_lwr(jl,laytrop(1), klev,                          &
                            rrpk_rat_h2oco2(1:jl,:,:), colco2(1:jl,:), colh2o(1:jl,:),  coldry(1:jl,:), &
                            rrpk_fac0(1:jl,:,:), rrpk_fac1(1:jl,:,:), minorfrac(1:jl,:), &
                            selffac(1:jl,:),selffrac(1:jl,:),forfac(1:jl,:),forfrac(1:jl,:), &
                            jp(1:jl,:), rrpk_jt(1:jl,:,:), (igpt-ngs(ngb(igpt)-1)), indself(1:jl,:), &
                            indfor(1:jl,:), &
                            rrpk_taug04(1:jl,:),fracs(1:jl,:,ig))
                        call taumol04_upr(jl,laytrop(1), klev,                          &
                            rrpk_rat_o3co2(1:jl,:,:), colco2(1:jl,:), colo3(1:jl,:),  coldry(1:jl,:), &
                            rrpk_fac0(1:jl,:,:), rrpk_fac1(1:jl,:,:), minorfrac(1:jl,:), &
                            forfac(1:jl,:),forfrac(1:jl,:),           &
                            jp(1:jl,:), rrpk_jt(1:jl,:,:), (igpt-ngs(ngb(igpt)-1)), &
                            indfor(1:jl,:), &
                            rrpk_taug04(1:jl,:),fracs(1:jl,:,ig))
                ENDIF
                DO jl = 1, kproma
                    ib = ibs(jl, ig)
                    igpt = igs(jl, ig)
                    !
                    ! Gas concentrations in colxx variables are normalized by 1.e-20_wp in lrtm_coeffs
                    !   CFC gas concentrations (wx) need the same normalization
                    !   Per Eli Mlawer the k values used in gas optics tables have been multiplied by 1e20
                    wx_loc(:,:) = 1.e-20_wp * wx(jl,:,:)
                    IF (ngb(igpt) == 3) THEN
                        taug = rrpk_taug03(jl,:)
                    ELSEIF (ngb(igpt) == 4) THEN
                        taug = rrpk_taug04(jl,:)
                    ELSE
                        CALL gas_optics_lw(klev, igpt, play        (jl,:), wx_loc    (:,:), coldry      (jl,:), laytrop     (jl), jp  &
                            (jl,:), jt        (jl,:), jt1         (jl,:), colh2o      (jl,:), colco2      (jl,:), colo3     (jl,:)&
                            , coln2o      (jl,:), colco       (jl,:), colch4      (jl,:), colo2     (jl,:), colbrd      (jl,:), fac00     &
                            (jl,:), fac01       (jl,:), fac10     (jl,:), fac11       (jl,:), rat_h2oco2  (jl,:), rat_h2oco2_1(jl,:), &
                            rat_h2oo3 (jl,:), rat_h2oo3_1 (jl,:), rat_h2on2o  (jl,:), rat_h2on2o_1(jl,:), rat_h2och4(jl,:), rat_h2och4_1(&
                            jl,:), rat_n2oco2  (jl,:), rat_n2oco2_1(jl,:), rat_o3co2 (jl,:), rat_o3co2_1 (jl,:), selffac     (jl,:), &
                            selffrac    (jl,:), indself   (jl,:), forfac      (jl,:), forfrac     (jl,:), indfor      (jl,:), minorfrac (&
                            jl,:), scaleminor  (jl,:), scaleminorn2(jl,:), indminor    (jl,:), fracs     (jl,:,ig), taug)
                    END IF
                    DO jk = 1, klev
                        taut(jl,jk,ig) = taug(jk) + tauaer(jl,jk,ib)
                    END DO
                END DO  ! Loop over columns
            END DO  ! Loop over g point samples - done with gas optical depth calculations
            tautot(1:kproma,:,:) = taut(1:kproma,:,:) + smp_tau(1:kproma,:,:) ! All-sky optical depth. Mask for 0 cloud optical depth?
            !
            ! ---  3.0 Compute radiative transfer.
            ! --------------------------------
            !
            ! Initialize fluxes to zero
            !
            uflx(1:kproma,0:klev) = 0.0_wp
            dflx(1:kproma,0:klev) = 0.0_wp
            uflxc(1:kproma,0:klev) = 0.0_wp
            dflxc(1:kproma,0:klev) = 0.0_wp
            !
            ! Planck function in each band at layers and boundaries
            !
            !IBM* ASSERT(NODEPS)
            DO ig = 1, nbndlw
                planklay(1:kproma,1:klev,ig) = planckfunction(tlay(1:kproma,1:klev  ),ig)
                planklev(1:kproma,0:klev,ig) = planckfunction(tlev(1:kproma,1:klev+1),ig)
                plankbnd(1:kproma,       ig) = planckfunction(tsfc(1:kproma         ),ig)
            END DO
            !
            ! Precipitable water vapor in each column - this can affect the integration angle secdiff
            !
            pwvcm(1:kproma) = ((amw * sum(wkl(1:kproma,1,1:klev), dim=2)) /                        (amd * sum(coldry(1:kproma,&
            1:klev) + wkl(1:kproma,1,1:klev), dim=2))) *                       (1.e3_wp * psfc(1:kproma)) / (1.e2_wp * grav)
            !
            ! Compute radiative transfer for each set of samples
            !
            DO ig = 1, n_gpts_ts
                secdiff(1:kproma) = find_secdiff(ibs(1:kproma, ig), pwvcm(1:kproma))
                !IBM* ASSERT(NODEPS)
                DO jl = 1, kproma
                    ib = ibs(jl,ig)
                    layplnk(jl,1:klev) = planklay(jl,1:klev,ib)
                    levplnk(jl,0:klev) = planklev(jl,0:klev,ib)
                    bndplnk(jl) = plankbnd(jl,       ib)
                    srfemis(jl) = emis    (jl,       ib)
                END DO
                !
                ! All sky fluxes
                !
                CALL lrtm_solver(kproma, kbdim, klev, tautot(:,:,ig), layplnk, levplnk, fracs(:,:,ig), secdiff, bndplnk, srfemis, &
                zgpfu, zgpfd)
                uflx(1:kproma,0:klev) = uflx (1:kproma,0:klev)                              + zgpfu(1:kproma,0:klev) * gpt_scaling
                dflx(1:kproma,0:klev) = dflx (1:kproma,0:klev)                              + zgpfd(1:kproma,0:klev) * gpt_scaling
                !
                ! Clear-sky fluxes
                !
                IF (l_do_sep_clear_sky) THEN
                    !
                    ! Remove clouds and do second RT calculation
                    !
                    CALL lrtm_solver(kproma, kbdim, klev, taut (:,:,ig), layplnk, levplnk, fracs(:,:,ig), secdiff, bndplnk, &
                    srfemis, zgpcu, zgpcd)
                    uflxc(1:kproma,0:klev) = uflxc(1:kproma,0:klev) + zgpcu(1:kproma,0:klev) * gpt_scaling
                    dflxc(1:kproma,0:klev) = dflxc(1:kproma,0:klev) + zgpcd(1:kproma,0:klev) * gpt_scaling
                    ELSE
                    !
                    ! Accumulate fluxes by excluding cloudy subcolumns, weighting to account for smaller sample size
                    !
                    !IBM* ASSERT(NODEPS)
                    DO jk = 0, klev
                        uflxc(1:kproma,jk) = uflxc(1:kproma,jk)                                                                   &
                                           + merge(0._wp,                                                                         &
                                                           zgpfu(1:kproma,jk) * clrsky_scaling(1:kproma),                         &
                                                                         colcldmask(1:kproma,ig))
                        dflxc(1:kproma,jk) = dflxc(1:kproma,jk)                                                                   &
                                           + merge(0._wp,                                                                         &
                                                           zgpfd(1:kproma,jk) * clrsky_scaling(1:kproma),                         &
                                                                         colcldmask(1:kproma,ig))
                    END DO
                END IF
            END DO  ! Loop over samples
            !
            ! ---  3.1 If computing clear-sky fluxes from samples, flag any columns where all samples were cloudy
            !
            ! --------------------------------
            IF (.not. l_do_sep_clear_sky) THEN
                !IBM* ASSERT(NODEPS)
                DO jl = 1, kproma
                    IF (all(colcldmask(jl,:))) THEN
                        uflxc(jl,0:klev) = rad_undef
                        dflxc(jl,0:klev) = rad_undef
                    END IF
                END DO
            END IF
        END SUBROUTINE lrtm
        !----------------------------------------------------------------------------

        elemental FUNCTION planckfunction(temp, band)
            !
            ! Compute the blackbody emission in a given band as a function of temperature
            !
            REAL(KIND=wp), intent(in) :: temp
            INTEGER, intent(in) :: band
            REAL(KIND=wp) :: planckfunction
            INTEGER :: index
            REAL(KIND=wp) :: fraction
            index = min(max(1, int(temp - 159._wp)),180)
            fraction = temp - 159._wp - float(index)
            planckfunction = totplanck(index, band)                    + fraction * (totplanck(index+1, band) - totplanck(index, &
            band))
            planckfunction = planckfunction * delwave(band)
        END FUNCTION planckfunction
    END MODULE mo_lrtm_driver

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
c-----------------------------------------------------------------------
      program re2torea

c     this program will take a nek all-ascii rea file and
c     extract geometry, bcs, and curve data to a new binary re2 file, plus
c     an ascii rea file for just the parameters

      character*80 fin
      character*1  fin1(80)
      equivalence (fin1,fin)

      character*80 fout
      character*1  fout1(80)
      equivalence (fout1,fout)

      character*80 fbin
      character*1  fbin1(80)
      equivalence (fbin1,fbin)

      character*80 string
      character*1  string1(80)
      equivalence (string,string1)

      character*80 sstring

      parameter(nelm=MAXNEL)
      real*8 x(8),y(8),z(8)
      real*8 bc(5,6,nelm)
      integer id, jd, wdsizi
      logical ifflow,ifheat

      character*80 hdr
      real*4 test
      character*5 version

      integer buf(30)

      real*8 rbc,rcurve
      character*3 cbc(6,nelm)
      logical ifbswap,if_byte_swap_test,ifexist

c     for workstation:
      in = 5

c     Get file name
      write(6,*) 'Input .rea & .re2 name:'
      call blank(fin,80)
      read(in,80) fin
      len = ltrunc(fin,80)
   80 format(a80)

c     Get file name
      write(6,*) 'Input .rea output name:'
      call blank(fout,80)
      read(in,80) fout
      lou = ltrunc(fout,80)

      fbin = fin
      call chcopy(fin1(len+1),'.rea',4)
      call chcopy(fout1(lou+1),'.rea',4)
      call chcopy(fbin1(len+1),'.re2'//CHAR(0),5)

      ifexist = .false.
      inquire(file=fin, exist=ifexist)
      if(.not.ifexist)then
        write(6,*)'could not find file ',(fin1(i),i=1,len+4),', abort.'
        stop
      endif
      ifexist = .false.
      inquire(file=fbin, exist=ifexist)
      if(.not.ifexist)then
        write(6,*)'could not find file ',(fbin1(i),i=1,len+4),', abort.'
        stop
      endif

      open(unit=10, file=fin)
      open(unit=11, file=fout)
      call byte_open(fbin,ierr)

      npscal = 0
      call scanout(sstring,'PARAMETERS FOLL',15,10,11)
      read(sstring,*) nparam
      write(11,80) sstring
      do ii=1,nparam
         read(10,80) sstring
         if(ii.eq.23) read(sstring,*) rnpscal
         write(11,80) sstring
      enddo
      npscal=rnpscal

      call scanout(sstring,'LOGICAL',7,10,11)
      read(sstring,*) nlogical
      write(11,80) sstring

      ifflow=.false.
      ifheat=.false.
      do ii = 1,nlogical
         read(10,80) sstring
         if(indx1(sstring,'IFFLOW',6).ne.0) read(sstring,*) ifflow
         if(indx1(sstring,'IFHEAT',6).ne.0) read(sstring,*) ifheat
         write(11,80) sstring
      enddo

      call blank(hdr,80)
      call byte_read(hdr,20,ierr)
      read (hdr,1) version,nelt,ndim,nelv
      nelt=abs(nelt)

      wdsizi = 4
      if(version.eq.'#v002') wdsizi = 8
      if(version.eq.'#v003') wdsizi = 8

      write(11,*)
     & ' 30.0000  20.0000  -6.00000  -10.0000 XFAC,YFAC,XZERO,YZERO'
      write(11,*)
     & '**MESH DATA** 6 lines are X,Y,Z;X,Y,Z. Columns corners 1-4;5-8'

      write(11,*) nelt,ndim,nelv

 1    format(a5,i9,i3,i9)

      call byte_read(test,1,ierr)
      ifbswap=if_byte_swap_test(test)

      print *, nelt, ndim, nelv

c mesh

      do ie=1,nelt

         nwds = (1 + ndim*(2**ndim))*(wdsizi/4)
         call byte_read(buf,nwds,ierr)
         call buf_to_xyz(buf,igroup,x,y,z,ifbswap,ndim,wdsizi)

         if(nelt.lt.100000) then
           write (11,12) ie, ie, 'a', igroup
         else
           write (11,'(A,I12,A,I1)')
     &        '  ELEMENT ', ie, '  GROUP  ', igroup
         endif
         write (11,15)   (x(k),k=1,4)
         write (11,15)   (y(k),k=1,4)

         if(ndim.eq.3) then
           write (11,15)   (z(k),k=1,4)
           write (11,15)   (x(k),k=5,8)
           write (11,15)   (y(k),k=5,8)
           write (11,15)   (z(k),k=5,8)
         endif
      enddo

 12   format(
     $     '            ELEMENT',i12,' [',i5,a1,']    GROUP     ',i1)

 15   format(4g14.6)

c curved sides

      if(wdsizi.eq.8) then
         call byte_read(rcurve,2,ierr)
         if(ifbswap) call byte_reverse8(rcurve,2,ierr)
         ncurve=rcurve
      else
         call byte_read(ncurve,1,ierr)
         if(ifbswap) call byte_reverse(ncurve,nwds,ierr)
      endif
      write(6,*) 'number of curved sides', ncurve

      write(11,13)
      write(11,14) ncurve
 13   format(' ***** CURVED SIDE DATA *****')
 14   format(i12
     $     ,' Curved sides follow IEDGE,IEL,CURVE(I),I=1,5, CCURVE')

      nwds = (2 + 1 + 5)*(wdsizi/4)

      IF (NCURVE.GT.0) THEN
         DO 50 ICURVE=1,NCURVE
            call byte_read(buf,nwds,ierr)
            call buf_to_curve(buf,nelt,ifbswap,wdsizi)
 50      CONTINUE
      endif

c boundary conditions
      nface=4
      if(ndim.eq.3) nface=6

      init_bc=2
      if(ifflow) init_bc=1
      nfld=1
      if(ifheat) nfld=2+npscal

      write(11,*) '***** BOUNDARY CONDITIONS *****'

      do ifld=init_bc,nfld
         if(ifld.eq.1) write(11,*)
     &    '***** FLUID   BOUNDARY CONDITIONS *****'
         if(ifld.eq.2.and.init_bc.eq.2) write(11,*)
     &    '***** NO FLUID   BOUNDARY CONDITIONS *****'
         if(ifld.eq.2) write(11,*)
     &    '***** THERMAL  BOUNDARY CONDITIONS *****'
         if(ifld.gt.2) write(11,*)
     &    '***** PASSIVE SCALAR  BOUNDARY CONDITIONS *****'

         nel=nelv
         if(ifld.ge.2) nel=nelt

         do ie = 1,nel
            do k=1,nface
               cbc(k,ie) = 'E  '
               do j=1,5
                 bc(j,k,ie) =0.0
               enddo
            enddo
         enddo

         if(wdsizi.eq.8) then
           call byte_read(rbc,2,ierr)
           if (ifbswap) call byte_reverse8(rbc,2,ierr)
           nbc = rbc
         else
           call byte_read(nbc,1,ierr)
           if (ifbswap) call byte_reverse(nbc,1,ierr)
         endif

         write(6,*) ifld, nbc, 'number of bc'

         do ie = 1,nbc
            call byte_read(buf,nwds,ierr)
            call buf_to_bc(cbc,bc,buf,nelm,ifbswap,wdsizi)
         enddo

         do ie = 1,nel
            do k=1,nface

              if (nelt.lt.1 000) then
                  write (11,20) cbc(k,ie),ie,k,(bc(j,k,ie),j=1,5)
              elseif(nelt.lt.100 000) then
                  write (11,21) cbc(k,ie),ie,k,(bc(j,k,ie),j=1,5)
              elseif(nelt.lt.1 000 000) then
                  write (11,22) cbc(k,ie),ie,(bc(j,k,ie),j=1,5)
              else
                  write (11,23) cbc(k,ie),ie,(bc(j,k,ie),j=1,5)
              endif
            enddo
         enddo
      enddo

      if(nfld.eq.1.and.init_bc.eq.1) write(11,*)
     &    '***** NO THERMAL BOUNDARY CONDITIONS *****'

   20 FORMAT(1x,A3,2I3,5G14.6)
   21 FORMAT(1x,A3,i5,i1,5G14.6)
   22 FORMAT(1x,A3,i6,5G14.7)
   23 FORMAT(1x,A3,i12,5G18.11)

      rewind(10)
      call scanout(sstring,'PRESOLVE',8,10,99)
      write(11,80) sstring
      call scanout(string,'endendend',9,10,11)

      close (unit=10)
      close (unit=11)
      call byte_close (fbin,ierr)

  999 continue
      stop
 9999 continue
      write(6,*) 'could not find file "indat". abort.'
      stop
      end
c-----------------------------------------------------------------------
      subroutine exitt

      stop

      return
      end
c-----------------------------------------------------------------------
      logical function if_byte_swap_test(bytetest)

      real*4 bytetest,test2
      real*4 test_pattern
      save   test_pattern

      test_pattern = 6.54321
      eps          = 0.00020
      etest        = abs(test_pattern-bytetest)
      if_byte_swap_test = .true.
      if (etest.le.eps) if_byte_swap_test = .false.

      test2 = bytetest
      call byte_reverse(test2,1,ierr)
      write(6,*) 'byte swap:',if_byte_swap_test,bytetest,test2
      return
      end
c-----------------------------------------------------------------------
      subroutine buf_to_xyz(buf,igroup,x,y,z,ifbswap,ndim,wdsizi)

      logical ifbswap
      integer buf(0:49),wdsizi
      real*8 x(8),y(8),z(8)

      nwds = (1 + ndim*(2**ndim))*(wdsizi/4) ! group + 2x4 for 2d, 3x8 for 3d

      if     (ifbswap.and.wdsizi.eq.8) then
          call byte_reverse8(buf,nwds,ierr)
      elseif (ifbswap.and.wdsizi.eq.4) then
          call byte_reverse (buf,nwds,ierr)
      endif

      if(wdsizi.eq.8) then
         call copyi4(igroup,buf(0),1) !0-1
         if (ndim.eq.3) then
            call copy8(x,buf( 2),8) !2 --17
            call copy8(y,buf(18),8) !18--33
            call copy8(z,buf(34),8) !34--49
         else
            call copy8(x,buf( 2),4) !2 --9
            call copy8(y,buf(10),4) !10--17
         endif
      else
         igroup = buf(0)
         if (ndim.eq.3) then
            call copy48(x,buf( 1),8)
            call copy48(y,buf( 9),8)
            call copy48(z,buf(17),8)
         else
            call copy48(x,buf( 1),4)
            call copy48(y,buf( 5),4)
         endif
      endif

      return
      end

c-----------------------------------------------------------------------
      subroutine buf_to_curve(buf,nel,ifbswap,wdsizi)

      integer e,f,wdsizi,buf(30)
      logical ifbswap
      character*1 ans
      real*8 r1,r2,r3,r4,r5

      nwds = (2 + 1 + 5)*(wdsizi/4)

      if     (ifbswap.and.wdsizi.eq.8) then
          call byte_reverse8(buf,nwds-2,ierr)
      elseif (ifbswap.and.wdsizi.eq.4) then
          call byte_reverse (buf,nwds-1,ierr)
      endif

      if(wdsizi.eq.8) then
        call copyi4(ieg,buf(1),1) !1-2

        call copyi4(iedg,buf(3),1) !3-4

        call copy8(r1,buf(5) ,1)
        call copy8(r2,buf(7) ,1)
        call copy8(r3,buf(9) ,1)
        call copy8(r4,buf(11),1)
        call copy8(r5,buf(13),1)
        call chcopy(ans,buf(15),1)
      else
        ieg   = buf(1)
        iedg  = buf(2)

        call copy48(r1,buf(3))
        call copy48(r2,buf(4))
        call copy48(r3,buf(5))
        call copy48(r4,buf(6))
        call copy48(r5,buf(7))
        call chcopy(ans,buf(8),1)
      endif
      IF (nel.LT.1000) THEN
        write(11,60) IEDG,IEG,R1,R2,R3,R4,R5,ANS
      ELSEIF (nel.LT.1 000 000) then
        write(11,61) IEDG,IEG,R1,R2,R3,R4,R5,ANS
      ELSE
        write(11,62) IEDG,IEG,R1,R2,R3,R4,R5,ANS
      ENDIF
 60   FORMAT(I3,I3 ,5G14.6,1X,A1)
 61   FORMAT(I2,I6 ,5G14.6,1X,A1)
 62   format(i2,i12,5g14.6,1x,a1)

      return
      end
c-----------------------------------------------------------------------
      subroutine buf_to_bc(cbc,bc,buf,nelm,ifbswap,wdsizi)

      character*3 cbc(6,nelm)
      real*8       bc(5,6,nelm)

      integer e,f,buf(30),wdsizi
      logical ifbswap

      nwds = (2 + 1 + 5)*(wdsizi/4)
      if     (ifbswap.and.wdsizi.eq.8) then
          call byte_reverse8(buf,nwds-2,ierr)
      elseif (ifbswap.and.wdsizi.eq.4) then
          call byte_reverse (buf,nwds-1,ierr)
      endif

      if(wdsizi.eq.8) then
        call copyi4(e,buf(1),1) !1-2

        call copyi4(f,buf(3),1) !3-4

        call copy8(bc(1,f,e),buf(5),5) !5--14
        call chcopy(cbc( f,e),buf(15),3)!15-16

      else
        e  = buf(1)
        f  = buf(2)

        call copy48 ( bc(1,f,e),buf(3),5)
        call chcopy (cbc(  f,e),buf(8),3)

      endif
      return
      end
c-----------------------------------------------------------------------

ext: f

max_stars_repo_name: RonRahaman/nekRS

max_issues_repo_name: neams-th-coe/nekRS

max_forks_repo_name: neams-th-coe/nekRS