!        Generated by TAPENADE     (INRIA, Ecuador team)
!  tapenade 3.x
!
!  Differentiation of forward_iter in reverse (adjoint) mode (with options noISIZE i8):
!   gradient     of useful results: *sdata
!   with respect to varying inputs: *d *e *f *concold *g *temp
!                *w *headold *x *sdata *head *dbc_data *bcperiod
!                *tempold *propunit *presold *conc *pres *a *b
!                *c
!   Plus diff mem management of: d:in e:in f:in concold:in g:in
!                temp:in w:in headold:in x:in sdata:in head:in
!                dbc_data:in bcperiod:in tempold:in propunit:in
!                tsal:in presold:in conc:in pres:in simtime:in
!                a:in b:in c:in
! MIT License
!
! Copyright (c) 2020 SHEMAT-Suite
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.
!> @brief time discretisation loop
!> @param[in] iter_out inverse iteration, SM realisation
!> @param[in] simtime_run start time of the simulation
!> @param[in] simtime_end finish time of the simulation
!> @param[in] iseed 0: FW simulation, 1 .. <mpara>: AD seeding index
!> @param[in] ismpl local sample index
!> @details
!> In-a-Nutshell description of this subroutine: \n
!> - Preprocessing before time step loop, initial variable values,
!>   monitoring output, extra steady-state initialisation, status_log\n
!> - Time loop: \n
!>   - before computations: time stepping, saving old variable arrays,
!>     output \n
!>   - computation: calling `forward_wrapper.f90`
!>   - after computations: save simulated data, update simtime,
!>     output, check divergence for variable step size \n
!> - Postprocessing: standard output
SUBROUTINE FORWARD_ITER_AD(simtime_run, simtime_end, iter_out, iseed, &
& ismpl)
  use arrays

  USE ARRAYS_AD

  USE MOD_GENRL
  USE MOD_GENRLC
  use mod_time

  USE MOD_TIME_AD

  USE MOD_LINFOS
  IMPLICIT NONE
! wrapper for output
!          call write_outt(deltt,ismpl)
! monitoring output
! initial monitoring output
! local sample index
  INTEGER :: ismpl
! iter_out: inverse iteration, SM realisation
  INTEGER :: iter_out, iseed
! Time step index
  INTEGER :: itimestep
! Size of a time period
  DOUBLE PRECISION :: deltt
! Start time of the simulation
  DOUBLE PRECISION, INTENT(IN) :: simtime_run
! Finish time of the simulation
  DOUBLE PRECISION, INTENT(IN) :: simtime_end
  DOUBLE PRECISION, EXTERNAL :: DELTAT
  INTEGER, EXTERNAL :: LBLANK
  INTEGER :: branch
  INTEGER :: res
  INTEGER :: res0
  INTEGER :: ad_count
  INTEGER :: i
! Preprocessing
! -------------
! initial values for some variables/arrays
  flag_1st_timestep(ismpl) = 0
  itimestep = itimestep_0
  simtime(ismpl) = simtime_run
  deltt = DELTAT(simtime(ismpl), ismpl)
  tr_switch(ismpl) = .true.
  iter_nlold = maxiter_nl/2
!          call write_monitor(1,ismpl)
!          call write_monitor_user(1,ismpl)
! runmode 2: extra steady state initialisation
  IF (transient .AND. runmode .EQ. 2) THEN
    tr_switch(ismpl) = .false.
    IF (iseed .EQ. 0 .AND. linfos(2) .GE. 0) WRITE(*, '(1A)') &
&                            '  [I] : extra steady state initialisation'
    CALL PUSHINTEGER8(iter_nlold)
    IF (ALLOCATED(c)) THEN
      CALL PUSHREAL8ARRAY(c, SIZE(c, 1)*SIZE(c, 2)*SIZE(c, 3)*SIZE(c, 4)&
&                  )
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(b)) THEN
      CALL PUSHREAL8ARRAY(b, SIZE(b, 1)*SIZE(b, 2)*SIZE(b, 3)*SIZE(b, 4)&
&                  )
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(a)) THEN
      CALL PUSHREAL8ARRAY(a, SIZE(a, 1)*SIZE(a, 2)*SIZE(a, 3)*SIZE(a, 4)&
&                  )
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(simtime)) THEN
      CALL PUSHREAL8ARRAY(simtime, SIZE(simtime, 1))
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(pres)) THEN
      CALL PUSHREAL8ARRAY(pres, SIZE(pres, 1)*SIZE(pres, 2)*SIZE(pres, 3&
&                   )*SIZE(pres, 4))
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(conc)) THEN
      CALL PUSHREAL8ARRAY(conc, SIZE(conc, 1)*SIZE(conc, 2)*SIZE(conc, 3&
&                   )*SIZE(conc, 4)*SIZE(conc, 5))
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(presold)) THEN
      CALL PUSHREAL8ARRAY(presold, SIZE(presold, 1)*SIZE(presold, 2)*&
&                   SIZE(presold, 3))
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(tsal)) THEN
      CALL PUSHREAL8ARRAY(tsal, SIZE(tsal, 1)*SIZE(tsal, 2)*SIZE(tsal, 3&
&                   )*SIZE(tsal, 4))
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(delt_count)) THEN
      CALL PUSHINTEGER8ARRAY(delt_count, SIZE(delt_count, 1))
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(flag_delt)) THEN
      CALL PUSHINTEGER8ARRAY(flag_delt, SIZE(flag_delt, 1))
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(flag_1st_timestep)) THEN
      CALL PUSHINTEGER8ARRAY(flag_1st_timestep, SIZE(flag_1st_timestep, &
&                      1))
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(tempold)) THEN
      CALL PUSHREAL8ARRAY(tempold, SIZE(tempold, 1)*SIZE(tempold, 2)*&
&                   SIZE(tempold, 3))
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(omp_iglobal)) THEN
      CALL PUSHINTEGER8ARRAY(omp_iglobal, SIZE(omp_iglobal, 1)*SIZE(&
&                      omp_iglobal, 2)*SIZE(omp_iglobal, 3))
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(dbc_data)) THEN
      CALL PUSHREAL8ARRAY(dbc_data, SIZE(dbc_data, 1)*SIZE(dbc_data, 2)*&
&                   SIZE(dbc_data, 3))
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(delt_old)) THEN
      CALL PUSHREAL8ARRAY(delt_old, SIZE(delt_old, 1))
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(omp_dglobal)) THEN
      CALL PUSHREAL8ARRAY(omp_dglobal, SIZE(omp_dglobal, 1)*SIZE(&
&                   omp_dglobal, 2)*SIZE(omp_dglobal, 3))
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(head)) THEN
      CALL PUSHREAL8ARRAY(head, SIZE(head, 1)*SIZE(head, 2)*SIZE(head, 3&
&                   )*SIZE(head, 4))
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(conv_ipos)) THEN
      CALL PUSHINTEGER8ARRAY(conv_ipos, SIZE(conv_ipos, 1)*SIZE(&
&                      conv_ipos, 2))
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(x)) THEN
      CALL PUSHREAL8ARRAY(x, SIZE(x, 1)*SIZE(x, 2)*SIZE(x, 3)*SIZE(x, 4)&
&                  )
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(w)) THEN
      CALL PUSHREAL8ARRAY(w, SIZE(w, 1)*SIZE(w, 2)*SIZE(w, 3)*SIZE(w, 4)&
&                  )
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(lcon)) THEN
      CALL PUSHBOOLEANARRAY(lcon, SIZE(lcon, 1)*SIZE(lcon, 2))
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(conc_conv)) THEN
      CALL PUSHREAL8ARRAY(conc_conv, SIZE(conc_conv, 1)*SIZE(conc_conv, &
&                   2))
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(conv_chlen)) THEN
      CALL PUSHINTEGER8ARRAY(conv_chlen, SIZE(conv_chlen, 1)*SIZE(&
&                      conv_chlen, 2))
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(temp)) THEN
      CALL PUSHREAL8ARRAY(temp, SIZE(temp, 1)*SIZE(temp, 2)*SIZE(temp, 3&
&                   )*SIZE(temp, 4))
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(g)) THEN
      CALL PUSHREAL8ARRAY(g, SIZE(g, 1)*SIZE(g, 2)*SIZE(g, 3)*SIZE(g, 4)&
&                  )
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(concold)) THEN
      CALL PUSHREAL8ARRAY(concold, SIZE(concold, 1)*SIZE(concold, 2)*&
&                   SIZE(concold, 3)*SIZE(concold, 4))
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(f)) THEN
      CALL PUSHREAL8ARRAY(f, SIZE(f, 1)*SIZE(f, 2)*SIZE(f, 3)*SIZE(f, 4)&
&                  )
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(conv_history)) THEN
      CALL PUSHREAL8ARRAY(conv_history, SIZE(conv_history, 1)*SIZE(&
&                   conv_history, 2)*SIZE(conv_history, 3))
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(e)) THEN
      CALL PUSHREAL8ARRAY(e, SIZE(e, 1)*SIZE(e, 2)*SIZE(e, 3)*SIZE(e, 4)&
&                  )
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    IF (ALLOCATED(d)) THEN
      CALL PUSHREAL8ARRAY(d, SIZE(d, 1)*SIZE(d, 2)*SIZE(d, 3)*SIZE(d, 4)&
&                  )
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHCONTROL1B(0)
    END IF
    CALL FORWARD_WRAPPER(itimestep, iseed, ismpl)
    IF (iseed .EQ. 0 .AND. linfos(2) .GE. 0) WRITE(*, '(1A)') &
&                                     '  [I] : normal transient process'
    CALL PUSHBOOLEAN(tr_switch(ismpl))
    tr_switch(ismpl) = .true.
    CALL PUSHCONTROL1B(0)
  ELSE
    CALL PUSHCONTROL1B(1)
  END IF
! Write to status_log
  IF (transient .AND. iseed .EQ. 0 .AND. (.NOT.write_iter_disable)) THEN
    OPEN(76, file=status_log, status='unknown', position='append')
    WRITE(76, fmt='(I8,1e14.6,1e14.6)') itimestep, deltt, simtime(ismpl)&
&   /tunit
    CLOSE(76)
  END IF
  ad_count = 1
! Time step loop for forward modeling
! -----------------------------------
 1000 IF (transient) THEN
! Advance time step
    CALL PUSHINTEGER8(itimestep)
    itimestep = itimestep + 1
! Initialize flag for variable time step size
    CALL PUSHINTEGER8(flag_delt(ismpl))
    flag_delt(ismpl) = 0
! Time stepping info to standard out
    IF (linfos(1) .GE. 1) THEN
      WRITE(*, *) ' '
      WRITE(*, '(1A,1I6)') '  >>>> new time step: ', itimestep
      WRITE(*, '(1A,1e16.8,1A,1e16.8)') '  >>>>     cum. time= ', (&
&     simtime(ismpl)+deltt)/tunit, '/', max_simtime/tunit
      WRITE(*, '(1A,1e16.8)') '  >>>>     time step= ', deltt/tunit
    END IF
! Save old time level
    CALL OLD_SAVE(cgen_time, ismpl)
    CALL PUSHCONTROL1B(0)
  ELSE
    CALL PUSHCONTROL1B(1)
  END IF
! ######### Forward Iteration ######
  CALL PUSHINTEGER8(iter_nlold)
  IF (ALLOCATED(c)) THEN
    CALL PUSHREAL8ARRAY(c, SIZE(c, 1)*SIZE(c, 2)*SIZE(c, 3)*SIZE(c, 4))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(b)) THEN
    CALL PUSHREAL8ARRAY(b, SIZE(b, 1)*SIZE(b, 2)*SIZE(b, 3)*SIZE(b, 4))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(a)) THEN
    CALL PUSHREAL8ARRAY(a, SIZE(a, 1)*SIZE(a, 2)*SIZE(a, 3)*SIZE(a, 4))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(simtime)) THEN
    CALL PUSHREAL8ARRAY(simtime, SIZE(simtime, 1))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(pres)) THEN
    CALL PUSHREAL8ARRAY(pres, SIZE(pres, 1)*SIZE(pres, 2)*SIZE(pres, 3)*&
&                 SIZE(pres, 4))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(conc)) THEN
    CALL PUSHREAL8ARRAY(conc, SIZE(conc, 1)*SIZE(conc, 2)*SIZE(conc, 3)*&
&                 SIZE(conc, 4)*SIZE(conc, 5))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(presold)) THEN
    CALL PUSHREAL8ARRAY(presold, SIZE(presold, 1)*SIZE(presold, 2)*SIZE(&
&                 presold, 3))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(tsal)) THEN
    CALL PUSHREAL8ARRAY(tsal, SIZE(tsal, 1)*SIZE(tsal, 2)*SIZE(tsal, 3)*&
&                 SIZE(tsal, 4))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(delt_count)) THEN
    CALL PUSHINTEGER8ARRAY(delt_count, SIZE(delt_count, 1))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(flag_delt)) THEN
    CALL PUSHINTEGER8ARRAY(flag_delt, SIZE(flag_delt, 1))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(flag_1st_timestep)) THEN
    CALL PUSHINTEGER8ARRAY(flag_1st_timestep, SIZE(flag_1st_timestep, 1)&
&                   )
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(tempold)) THEN
    CALL PUSHREAL8ARRAY(tempold, SIZE(tempold, 1)*SIZE(tempold, 2)*SIZE(&
&                 tempold, 3))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(omp_iglobal)) THEN
    CALL PUSHINTEGER8ARRAY(omp_iglobal, SIZE(omp_iglobal, 1)*SIZE(&
&                    omp_iglobal, 2)*SIZE(omp_iglobal, 3))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(dbc_data)) THEN
    CALL PUSHREAL8ARRAY(dbc_data, SIZE(dbc_data, 1)*SIZE(dbc_data, 2)*&
&                 SIZE(dbc_data, 3))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(delt_old)) THEN
    CALL PUSHREAL8ARRAY(delt_old, SIZE(delt_old, 1))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(omp_dglobal)) THEN
    CALL PUSHREAL8ARRAY(omp_dglobal, SIZE(omp_dglobal, 1)*SIZE(&
&                 omp_dglobal, 2)*SIZE(omp_dglobal, 3))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(head)) THEN
    CALL PUSHREAL8ARRAY(head, SIZE(head, 1)*SIZE(head, 2)*SIZE(head, 3)*&
&                 SIZE(head, 4))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(conv_ipos)) THEN
    CALL PUSHINTEGER8ARRAY(conv_ipos, SIZE(conv_ipos, 1)*SIZE(conv_ipos&
&                    , 2))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(x)) THEN
    CALL PUSHREAL8ARRAY(x, SIZE(x, 1)*SIZE(x, 2)*SIZE(x, 3)*SIZE(x, 4))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(w)) THEN
    CALL PUSHREAL8ARRAY(w, SIZE(w, 1)*SIZE(w, 2)*SIZE(w, 3)*SIZE(w, 4))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(lcon)) THEN
    CALL PUSHBOOLEANARRAY(lcon, SIZE(lcon, 1)*SIZE(lcon, 2))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(conc_conv)) THEN
    CALL PUSHREAL8ARRAY(conc_conv, SIZE(conc_conv, 1)*SIZE(conc_conv, 2)&
&                )
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(conv_chlen)) THEN
    CALL PUSHINTEGER8ARRAY(conv_chlen, SIZE(conv_chlen, 1)*SIZE(&
&                    conv_chlen, 2))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(temp)) THEN
    CALL PUSHREAL8ARRAY(temp, SIZE(temp, 1)*SIZE(temp, 2)*SIZE(temp, 3)*&
&                 SIZE(temp, 4))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(g)) THEN
    CALL PUSHREAL8ARRAY(g, SIZE(g, 1)*SIZE(g, 2)*SIZE(g, 3)*SIZE(g, 4))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(concold)) THEN
    CALL PUSHREAL8ARRAY(concold, SIZE(concold, 1)*SIZE(concold, 2)*SIZE(&
&                 concold, 3)*SIZE(concold, 4))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(f)) THEN
    CALL PUSHREAL8ARRAY(f, SIZE(f, 1)*SIZE(f, 2)*SIZE(f, 3)*SIZE(f, 4))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(conv_history)) THEN
    CALL PUSHREAL8ARRAY(conv_history, SIZE(conv_history, 1)*SIZE(&
&                 conv_history, 2)*SIZE(conv_history, 3))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(e)) THEN
    CALL PUSHREAL8ARRAY(e, SIZE(e, 1)*SIZE(e, 2)*SIZE(e, 3)*SIZE(e, 4))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(d)) THEN
    CALL PUSHREAL8ARRAY(d, SIZE(d, 1)*SIZE(d, 2)*SIZE(d, 3)*SIZE(d, 4))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  CALL FORWARD_WRAPPER(itimestep, iseed, ismpl)
! ##################################
! save and collect the computed values for:
! - comparison with 'ddata(*,cid_pv)' (observed data)
! - data-output (write_data.f)
  IF (ALLOCATED(simtime)) THEN
    CALL PUSHREAL8ARRAY(simtime, SIZE(simtime, 1))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(delt_count)) THEN
    CALL PUSHINTEGER8ARRAY(delt_count, SIZE(delt_count, 1))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(flag_1st_timestep)) THEN
    CALL PUSHINTEGER8ARRAY(flag_1st_timestep, SIZE(flag_1st_timestep, 1)&
&                   )
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(dbc_data)) THEN
    CALL PUSHREAL8ARRAY(dbc_data, SIZE(dbc_data, 1)*SIZE(dbc_data, 2)*&
&                 SIZE(dbc_data, 3))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  IF (ALLOCATED(delt_old)) THEN
    CALL PUSHREAL8ARRAY(delt_old, SIZE(delt_old, 1))
    CALL PUSHCONTROL1B(1)
  ELSE
    CALL PUSHCONTROL1B(0)
  END IF
  CALL SAVE_DATA(ismpl)
  IF (transient) THEN
! Update simulation time
    simtime(ismpl) = simtime(ismpl) + deltt
!            call write_monitor(2,ismpl)
!            call write_monitor_user(2,ismpl)
! Write to status_log
    IF (.NOT.write_iter_disable .AND. iseed .EQ. 0) THEN
! Status log info to standard out
      IF (linfos(1) .GE. 1) THEN
        CALL PUSHCONTROL1B(0)
        res = LBLANK(status_log)
        WRITE(*, '(3A)') '  [W] : "', status_log(1:res), '"'
      ELSE
        CALL PUSHCONTROL1B(0)
      END IF
      OPEN(76, file=status_log, status='unknown', position='append')
      WRITE(76, fmt='(I8,1e14.6,1e14.6)') itimestep, deltt, simtime(&
&     ismpl)/tunit
      CLOSE(76)
    ELSE
      CALL PUSHCONTROL1B(1)
    END IF
! Check for variable time stepping divergence flag
    IF (flag_delt(ismpl) .EQ. -2) THEN
! Restore old values of variable arrays if time step was
! halfed (restoring simtime is handled in "deltat")
      IF (ALLOCATED(pres)) THEN
        CALL PUSHREAL8ARRAY(pres, SIZE(pres, 1)*SIZE(pres, 2)*SIZE(pres&
&                     , 3)*SIZE(pres, 4))
        CALL PUSHCONTROL1B(1)
      ELSE
        CALL PUSHCONTROL1B(0)
      END IF
      IF (ALLOCATED(temp)) THEN
        CALL PUSHREAL8ARRAY(temp, SIZE(temp, 1)*SIZE(temp, 2)*SIZE(temp&
&                     , 3)*SIZE(temp, 4))
        CALL PUSHCONTROL1B(1)
      ELSE
        CALL PUSHCONTROL1B(0)
      END IF
      CALL OLD_RESTORE(cgen_time, ismpl)
      CALL PUSHCONTROL1B(0)
    ELSE
      CALL PUSHCONTROL1B(1)
    END IF
! Set variable time stepping divergence flag to zero to
! avoid double calling of deltat.
    CALL PUSHINTEGER8(flag_delt(ismpl))
    flag_delt(ismpl) = 0
    deltt = DELTAT(simtime(ismpl), ismpl)
    CALL PUSHCONTROL1B(0)
  ELSE
    CALL PUSHCONTROL1B(1)
  END IF
! Important: This if statement (with the goto) needs to be
! outside of the other "transient" scopes, to generate
! reverse-mode code!
  IF (transient .AND. simtime(ismpl) .LT. simtime_end) THEN
    ad_count = ad_count + 1
    GOTO 1000
  END IF
  CALL PUSHINTEGER8(ad_count)
! --------------- 1000: return to next time step ----------------------
! Postprocessing
! --------------
! Write to status_log and status_log_inv
  IF (transient .AND. runmode .GE. 2 .AND. (.NOT.write_iter_disable) &
&     .AND. iseed .EQ. 0) THEN
    res0 = LBLANK(status_log)
    WRITE(*, '(3A)') '  [W] : "', status_log(1:res0), '"'
    OPEN(76, file=status_log, status='unknown', position='append')
    WRITE(76, '(1A)') key_char//' transient end'
    CLOSE(76)
    OPEN(76, file=status_log_inv, status='unknown', position='append')
    WRITE(76, '(1A,I8,1e14.6,1e14.6)') key_char//' transient: ', &
&   itimestep, deltt, simtime(ismpl)/tunit
    CLOSE(76)
  END IF
! Standard output: steady state
  IF (iter_out .GT. 0 .AND. linfos(1) .GE. 1 .AND. (.NOT.transient) &
&     .AND. iseed .EQ. 0) WRITE(*, '(29X,1A)') &
&                         ' ===> leaving nonlinear iteration'
! Standard output: transient
  IF (linfos(1) .GE. 0 .AND. transient .AND. iseed .EQ. 0) WRITE(*, &
&                                                    '(1A,I8,1A,1e14.6)'&
&                                                         ) &
&                                           '  [I] : final time step = '&
&                                                          , itimestep, &
&                                                    ', simulation time'&
&                                                          , simtime(&
&                                                          ismpl)/tunit
  CALL POPINTEGER8(ad_count)
  DO i=1,ad_count
    IF (i .EQ. 1) THEN
      IF (ALLOCATED(d_ad)) d_ad = 0.D0
      IF (ALLOCATED(e_ad)) e_ad = 0.D0
      IF (ALLOCATED(f_ad)) f_ad = 0.D0
      IF (ALLOCATED(concold_ad)) concold_ad = 0.D0
      IF (ALLOCATED(g_ad)) g_ad = 0.D0
      IF (ALLOCATED(temp_ad2)) temp_ad2 = 0.D0
      IF (ALLOCATED(w_ad)) w_ad = 0.D0
      IF (ALLOCATED(headold_ad)) headold_ad = 0.D0
      IF (ALLOCATED(x_ad)) x_ad = 0.D0
      IF (ALLOCATED(head_ad)) head_ad = 0.D0
      IF (ALLOCATED(dbc_data_ad)) dbc_data_ad = 0.D0
      IF (ALLOCATED(bcperiod_ad)) bcperiod_ad = 0.D0
      IF (ALLOCATED(tempold_ad)) tempold_ad = 0.D0
      IF (ALLOCATED(propunit_ad)) propunit_ad = 0.D0
      IF (ALLOCATED(tsal_ad)) tsal_ad = 0.D0
      IF (ALLOCATED(presold_ad)) presold_ad = 0.D0
      IF (ALLOCATED(conc_ad)) conc_ad = 0.D0
      IF (ALLOCATED(pres_ad)) pres_ad = 0.D0
      IF (ALLOCATED(a_ad)) a_ad = 0.D0
      IF (ALLOCATED(b_ad)) b_ad = 0.D0
      IF (ALLOCATED(c_ad)) c_ad = 0.D0
    END IF
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 0) THEN
      CALL POPINTEGER8(flag_delt(ismpl))
      CALL POPCONTROL1B(branch)
      IF (branch .EQ. 0) THEN
        CALL POPCONTROL1B(branch)
        IF (branch .EQ. 1) CALL POPREAL8ARRAY(temp, SIZE(temp, 1)*SIZE(&
&                                       temp, 2)*SIZE(temp, 3)*SIZE(temp&
&                                       , 4))
        CALL POPCONTROL1B(branch)
        IF (branch .EQ. 1) CALL POPREAL8ARRAY(pres, SIZE(pres, 1)*SIZE(&
&                                       pres, 2)*SIZE(pres, 3)*SIZE(pres&
&                                       , 4))
        CALL OLD_RESTORE_AD(cgen_time, ismpl)
      END IF
      CALL POPCONTROL1B(branch)
    END IF
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(delt_old, SIZE(delt_old, 1))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(dbc_data, SIZE(dbc_data, 1)*&
&                                   SIZE(dbc_data, 2)*SIZE(dbc_data, 3))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPINTEGER8ARRAY(flag_1st_timestep, SIZE(&
&                                      flag_1st_timestep, 1))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPINTEGER8ARRAY(delt_count, SIZE(delt_count&
&                                      , 1))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(simtime, SIZE(simtime, 1))
    CALL SAVE_DATA_AD(ismpl)
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(d, SIZE(d, 1)*SIZE(d, 2)*SIZE(&
&                                   d, 3)*SIZE(d, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(e, SIZE(e, 1)*SIZE(e, 2)*SIZE(&
&                                   e, 3)*SIZE(e, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(conv_history, SIZE(&
&                                   conv_history, 1)*SIZE(conv_history, &
&                                   2)*SIZE(conv_history, 3))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(f, SIZE(f, 1)*SIZE(f, 2)*SIZE(&
&                                   f, 3)*SIZE(f, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(concold, SIZE(concold, 1)*SIZE&
&                                   (concold, 2)*SIZE(concold, 3)*SIZE(&
&                                   concold, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(g, SIZE(g, 1)*SIZE(g, 2)*SIZE(&
&                                   g, 3)*SIZE(g, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(temp, SIZE(temp, 1)*SIZE(temp&
&                                   , 2)*SIZE(temp, 3)*SIZE(temp, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPINTEGER8ARRAY(conv_chlen, SIZE(conv_chlen&
&                                      , 1)*SIZE(conv_chlen, 2))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(conc_conv, SIZE(conc_conv, 1)*&
&                                   SIZE(conc_conv, 2))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPBOOLEANARRAY(lcon, SIZE(lcon, 1)*SIZE(&
&                                     lcon, 2))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(w, SIZE(w, 1)*SIZE(w, 2)*SIZE(&
&                                   w, 3)*SIZE(w, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(x, SIZE(x, 1)*SIZE(x, 2)*SIZE(&
&                                   x, 3)*SIZE(x, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPINTEGER8ARRAY(conv_ipos, SIZE(conv_ipos, &
&                                      1)*SIZE(conv_ipos, 2))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(head, SIZE(head, 1)*SIZE(head&
&                                   , 2)*SIZE(head, 3)*SIZE(head, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(omp_dglobal, SIZE(omp_dglobal&
&                                   , 1)*SIZE(omp_dglobal, 2)*SIZE(&
&                                   omp_dglobal, 3))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(delt_old, SIZE(delt_old, 1))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(dbc_data, SIZE(dbc_data, 1)*&
&                                   SIZE(dbc_data, 2)*SIZE(dbc_data, 3))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPINTEGER8ARRAY(omp_iglobal, SIZE(&
&                                      omp_iglobal, 1)*SIZE(omp_iglobal&
&                                      , 2)*SIZE(omp_iglobal, 3))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(tempold, SIZE(tempold, 1)*SIZE&
&                                   (tempold, 2)*SIZE(tempold, 3))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPINTEGER8ARRAY(flag_1st_timestep, SIZE(&
&                                      flag_1st_timestep, 1))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPINTEGER8ARRAY(flag_delt, SIZE(flag_delt, &
&                                      1))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPINTEGER8ARRAY(delt_count, SIZE(delt_count&
&                                      , 1))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(tsal, SIZE(tsal, 1)*SIZE(tsal&
&                                   , 2)*SIZE(tsal, 3)*SIZE(tsal, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(presold, SIZE(presold, 1)*SIZE&
&                                   (presold, 2)*SIZE(presold, 3))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(conc, SIZE(conc, 1)*SIZE(conc&
&                                   , 2)*SIZE(conc, 3)*SIZE(conc, 4)*&
&                                   SIZE(conc, 5))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(pres, SIZE(pres, 1)*SIZE(pres&
&                                   , 2)*SIZE(pres, 3)*SIZE(pres, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(simtime, SIZE(simtime, 1))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(a, SIZE(a, 1)*SIZE(a, 2)*SIZE(&
&                                   a, 3)*SIZE(a, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(b, SIZE(b, 1)*SIZE(b, 2)*SIZE(&
&                                   b, 3)*SIZE(b, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(c, SIZE(c, 1)*SIZE(c, 2)*SIZE(&
&                                   c, 3)*SIZE(c, 4))
    CALL POPINTEGER8(iter_nlold)
    CALL FORWARD_WRAPPER_AD(itimestep, iseed, ismpl)
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 0) THEN
      CALL OLD_SAVE_AD(cgen_time, ismpl)
      CALL POPINTEGER8(flag_delt(ismpl))
      CALL POPINTEGER8(itimestep)
    END IF
  END DO
  CALL POPCONTROL1B(branch)
  IF (branch .EQ. 0) THEN
    CALL POPBOOLEAN(tr_switch(ismpl))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(d, SIZE(d, 1)*SIZE(d, 2)*SIZE(&
&                                   d, 3)*SIZE(d, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(e, SIZE(e, 1)*SIZE(e, 2)*SIZE(&
&                                   e, 3)*SIZE(e, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(conv_history, SIZE(&
&                                   conv_history, 1)*SIZE(conv_history, &
&                                   2)*SIZE(conv_history, 3))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(f, SIZE(f, 1)*SIZE(f, 2)*SIZE(&
&                                   f, 3)*SIZE(f, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(concold, SIZE(concold, 1)*SIZE&
&                                   (concold, 2)*SIZE(concold, 3)*SIZE(&
&                                   concold, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(g, SIZE(g, 1)*SIZE(g, 2)*SIZE(&
&                                   g, 3)*SIZE(g, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(temp, SIZE(temp, 1)*SIZE(temp&
&                                   , 2)*SIZE(temp, 3)*SIZE(temp, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPINTEGER8ARRAY(conv_chlen, SIZE(conv_chlen&
&                                      , 1)*SIZE(conv_chlen, 2))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(conc_conv, SIZE(conc_conv, 1)*&
&                                   SIZE(conc_conv, 2))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPBOOLEANARRAY(lcon, SIZE(lcon, 1)*SIZE(&
&                                     lcon, 2))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(w, SIZE(w, 1)*SIZE(w, 2)*SIZE(&
&                                   w, 3)*SIZE(w, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(x, SIZE(x, 1)*SIZE(x, 2)*SIZE(&
&                                   x, 3)*SIZE(x, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPINTEGER8ARRAY(conv_ipos, SIZE(conv_ipos, &
&                                      1)*SIZE(conv_ipos, 2))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(head, SIZE(head, 1)*SIZE(head&
&                                   , 2)*SIZE(head, 3)*SIZE(head, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(omp_dglobal, SIZE(omp_dglobal&
&                                   , 1)*SIZE(omp_dglobal, 2)*SIZE(&
&                                   omp_dglobal, 3))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(delt_old, SIZE(delt_old, 1))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(dbc_data, SIZE(dbc_data, 1)*&
&                                   SIZE(dbc_data, 2)*SIZE(dbc_data, 3))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPINTEGER8ARRAY(omp_iglobal, SIZE(&
&                                      omp_iglobal, 1)*SIZE(omp_iglobal&
&                                      , 2)*SIZE(omp_iglobal, 3))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(tempold, SIZE(tempold, 1)*SIZE&
&                                   (tempold, 2)*SIZE(tempold, 3))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPINTEGER8ARRAY(flag_1st_timestep, SIZE(&
&                                      flag_1st_timestep, 1))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPINTEGER8ARRAY(flag_delt, SIZE(flag_delt, &
&                                      1))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPINTEGER8ARRAY(delt_count, SIZE(delt_count&
&                                      , 1))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(tsal, SIZE(tsal, 1)*SIZE(tsal&
&                                   , 2)*SIZE(tsal, 3)*SIZE(tsal, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(presold, SIZE(presold, 1)*SIZE&
&                                   (presold, 2)*SIZE(presold, 3))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(conc, SIZE(conc, 1)*SIZE(conc&
&                                   , 2)*SIZE(conc, 3)*SIZE(conc, 4)*&
&                                   SIZE(conc, 5))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(pres, SIZE(pres, 1)*SIZE(pres&
&                                   , 2)*SIZE(pres, 3)*SIZE(pres, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(simtime, SIZE(simtime, 1))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(a, SIZE(a, 1)*SIZE(a, 2)*SIZE(&
&                                   a, 3)*SIZE(a, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(b, SIZE(b, 1)*SIZE(b, 2)*SIZE(&
&                                   b, 3)*SIZE(b, 4))
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 1) CALL POPREAL8ARRAY(c, SIZE(c, 1)*SIZE(c, 2)*SIZE(&
&                                   c, 3)*SIZE(c, 4))
    CALL POPINTEGER8(iter_nlold)
    CALL FORWARD_WRAPPER_AD(itimestep, iseed, ismpl)
  END IF
END SUBROUTINE FORWARD_ITER_AD


! { dg-do compile }
! { dg-additional-options "-fdump-tree-original" }
!
! PR fortran/66549
! The resolution of CVN in the middle CLWF's OpenMP construct was
! making the DO loop (wrongly) interpreted as an OpenMP-managed loop, leading
! to an ICE.
!
! Contributed by Andrew Benson <abensonca@gmail.com>.

module smfa
  type :: sgc
   contains
     procedure :: sla => sa
  end type sgc
  class(sgc), pointer :: sg_
  double precision, allocatable, dimension(:) :: vni
contains
  double precision function sa(self,i)
    class(sgc), intent(in   ) :: self
  end function sa
  subroutine cvn(sg_,vn)
    class(sgc), intent(inout) :: sg_
    double precision, intent(  out), dimension(:) :: vn
    integer :: i
    do i=1,2
       vn(i)= sg_%sla(i)
    end do
  end subroutine cvn
  subroutine clwf()
    !$omp parallel
    call cvn(sg_,vni)
    !$omp end parallel
  end subroutine clwf
end module smfa

! { dg-final { scan-tree-dump-times "#pragma\\s+omp\\s+parallel\\n" 1 "original" } }


program sort_galaxy
  use io
  use misc
  implicit none

  character(LEN=128)                         :: gal_list_filename, outfile, fileinfo='none', filename
  character(LEN=128)                         :: dm_halo_list_filename, associations_filename
  character(LEN=128)                         :: ramses_output_end, ramses_output_start, brick_file
  character(LEN=128)                         :: info_file_start, info_file_end, output_dir
  logical                                    :: verbose, ok
  integer                                    :: ngal, gal_cols, i, j, k, index
  integer                                    :: ndm_halo, dm_halo_cols
  integer                                    :: nassoc, assoc_cols
  integer                                    :: ellipticals = 0, spirals = 0, others = 0
  real(kind=4), dimension(:,:), allocatable  :: data_gal, data_halo
  ! contains association flag / dark matter id / mass / gal id / mass
  real(kind=4), dimension(:, :), allocatable :: data_associations
  real(kind=4), dimension(:), allocatable    :: kind ! 0 for elliptical, 1 for spirals, 2 for other
  real(kind=4)                               :: sigma
  real(kind=4)                               :: elliptical_threshold = 1.5d0 ! FIXME
  real(kind=4)                               :: spiral_threshold = 0.8d0 ! FIXME

  real(kind=8), dimension(:), allocatable    :: mDM, rvirDM
  real(kind=8), dimension(:), allocatable    :: mvirDM, TvirDM, hlevel
  real(kind=8), dimension(:, :), allocatable :: LDM, posDM
  integer, dimension(:), allocatable         :: idDM, cpu_list

  real(kind=4)                               :: aexp_tmp, age_univ
  real(kind=8), dimension(:), allocatable    :: X0, X1
  integer                                    :: nbodies, nb_of_halos, nb_of_subhalos, nDM

  integer                                    :: gal_id, halo_id, id_in_brick, cpu
  character(len=20)                          :: tmp_str
  integer, dimension(:), allocatable         :: gal_to_halo

  type(MEMBERS_T), dimension(:), allocatable :: members
  type(INFOS_T)                              :: infos_start, infos_end

  integer                                    :: ndim, nparts, nstars, particles_to_find
  real(kind=8), dimension(:,:), allocatable  :: pos, vel, tmp_pos, tmp_vel
  real(kind=8), dimension(:), allocatable    :: tmp_m, tmp_birth_date
  integer, dimension(:), allocatable         :: tmp_ids
  integer(kind=8), dimension(:), allocatable         :: halo_mask

  !-------------------------------------
  ! Read parameters
  !-------------------------------------
  call read_params()

  !-------------------------------------
  ! Read file
  !-------------------------------------
  if (gal_list_filename == 'none') then
     stop
  end if

  print*, 'Reading file "' // trim(gal_list_filename) // '"'
  call read_list_header(trim(gal_list_filename), ngal, gal_cols)
  allocate(data_gal(ngal, gal_cols+1))
  call read_list_data(ngal, gal_cols, data_gal(1:ngal, 1:gal_cols))
  print*, 'Got', ngal, 'galaxies in da pocket'

  print*, ''
  print*, 'Reading file "' // trim(associations_filename) // '"'
  call read_list_header(trim(associations_filename), nassoc, assoc_cols)
  allocate(data_associations(nassoc, assoc_cols))
  call read_list_data(nassoc, assoc_cols, data_associations)
  print*, 'Got', nassoc, 'associations in da pocket'

  print*, ''
  print*, 'Reading file "' // trim(dm_halo_list_filename) // '"'
  call read_list_header(trim(dm_halo_list_filename), ndm_halo, dm_halo_cols)
  allocate(data_halo(ndm_halo, dm_halo_cols))
  call read_list_data(ndm_halo, dm_halo_cols, data_halo)
  print*, 'Got', ndm_halo, 'halos in da pocket'

  allocate(kind(ngal))

  !-------------------------------------
  ! Processing
  !-------------------------------------
  print*, 'Sorting galaxies (ellipticals/spirals) using v/sigma criterion…'
  do i = 1, ngal
     sigma = 1d0/3d0 * sqrt(data_gal(i, 3)**2 + data_gal(i, 4)**2 + data_gal(i, 5)**2)
     if (sigma / data_gal(i, 2) > elliptical_threshold) then
        data_gal(i, gal_cols+1) = 0
        ellipticals = ellipticals + 1
     else if (sigma / data_gal(i, 2) < spiral_threshold) then
        data_gal(i, gal_cols+1) = 0
        spirals = spirals + 1
     else
        others = others + 1
        kind(i) = 2
     end if

  end do

  print*, 'Associating galaxies to halo…'
  allocate(gal_to_halo(int(maxval(data_gal(:, 1))))) ! this is the ids
  gal_to_halo = -1
  do i = 1, nassoc
     gal_id = int(data_associations(i, 4))
     halo_id = int(data_associations(i, 1))
     if (gal_id > 0) then
        gal_to_halo(gal_id) = halo_id
     end if
  end do
  i = -1

  print*, 'Some statistics…'
  ! Sort the galaxies
  print*, int(ellipticals*100./ngal), "% ellipticals"
  print*, int(spirals*100./ngal), "% spirals"
  print*, int(100 - (ellipticals + spirals)*100./ngal), "% others"

  !-------------------------------------
  ! Read brick at start output
  !-------------------------------------
  call read_info_headers(info_file_start, infos_start)
  call read_info_headers(info_file_start, infos_end)
  print*, "Reading brick file " // brick_file
  call read_brick_header(brick_file, infos_start, nbodies, aexp_tmp, age_univ,&
       nb_of_halos, nb_of_subhalos)
  nDM = nb_of_halos + nb_of_subhalos
  allocate(idDM(nDM))
  allocate(posDM(infos_start%ndim, nDM))
  allocate(rvirDM(nDM))
  allocate(mDM(nDM))
  allocate(mvirDM(nDM))
  allocate(TvirDM(nDM))
  allocate(hlevel(nDM))
  allocate(LDM(infos_start%ndim, nDM))
  allocate(members(nDM))
  print*, infos_start%boxlen, infos_start%unit_l/3.085677581e+24,  infos_start%boxlen*infos_start%unit_l/3.085677581e+24
  call read_brick_data(nDM, infos_start, .true., &
       & mDM, posDM, rvirDM, mvirDM, TvirDM,&
       & hlevel, LDM, idDM, members)

  print*, '', 'read!'
  ! print*, maxval(posDM(1,:)), maxval(posDM(2,:)), maxval(posDM(3,:))
  ! print*, minval(posDM(1,:)), minval(posDM(2,:)), minval(posDM(3,:))
  !-------------------------------------
  ! Iterate over each galaxy
  !-------------------------------------
  allocate(cpu_list(infos_start%ncpu))
  allocate(X0(infos_start%ndim))
  allocate(X1(infos_start%ndim))

  allocate(halo_mask(nDM))
  halo_mask = 0

  !$OMP DO SCHEDULE(guided) REDUCTION(+:halo_mask)
  do cpu = 1, infos_start%ncpu
     !-------------------------------------
     ! reading output of cpu
     !-------------------------------------
     write(tmp_str, '(i0.4)') cpu
     filename = trim(ramses_output_start) // '/part_00032.out0' // trim(tmp_str)

     print*, 'Reading', trim(filename)

     call read_particle_header(filename, ndim, nparts)
     allocate(tmp_pos(ndim, nparts))
     allocate(tmp_vel(ndim, nparts))
     allocate(tmp_m(nparts))
     allocate(tmp_ids(nparts))
     allocate(tmp_birth_date(nparts))
     call read_particle_data(ndim, nparts, nstars, tmp_pos, tmp_vel, tmp_m, &
          tmp_ids, tmp_birth_date)

     X0 = minval(tmp_pos, 2)
     X1 = maxval(tmp_pos, 2)

     do halo_id = 1, nDm
        if ( X0(1) < posDM(1, halo_id) - rvirDM(halo_id) .and. &
             X1(1) > posDM(1, halo_id) + rvirDM(halo_id) .and. &
             X0(2) < posDM(2, halo_id) - rvirDM(halo_id) .and. &
             X1(2) > posDM(2, halo_id) + rvirDM(halo_id) .and. &
             X0(3) < posDM(3, halo_id) - rvirDM(halo_id) .and. &
             X1(3) > posDM(3, halo_id) + rvirDM(halo_id) ) then
           ! Store as an integer in 'base ncpu', so that any value between
           ! 1 and ncpu gives the first cpu, value/ncpu%ncpu → 2nd CPU, value/ncpu**2%ncpu → 3rd value, …
           halo_mask(halo_id) = cpu + halo_mask(halo_id)*infos_start%ncpu
        end if

     end do
     deallocate(tmp_pos, tmp_vel, tmp_m, tmp_ids, tmp_birth_date)
  end do
  !$OMP END DO

  open(12, file='halo_to_cpu')
  write(12, '(2a64)') 'halo_id', 'cpu'
  do i = 1, nDM
     write(12, '(2i64)') i, halo_mask(i)
  end do
  close(12)
  open(12, file='halo_to_cpu.bin')
  write(12) halo_mask
  close(12)

  deallocate(halo_mask)

  !-------------------------------------
  ! Cleanup
  !-------------------------------------
  deallocate(cpu_list)
  deallocate(data_gal)
  deallocate(data_halo)
  deallocate(data_associations)
  deallocate(gal_to_halo)
  deallocate(idDM)
  deallocate(posDM)
  deallocate(rvirDM)
  deallocate(mDM)
  deallocate(mvirDM)
  deallocate(TvirDM)
  deallocate(hlevel)
  deallocate(LDM)
  do i = 1, nDM
     deallocate(members(i)%ids)
  end do
  deallocate(members)

contains
  subroutine read_params ()
    integer            :: i,n
    integer            :: iargc
    character(len=4)   :: opt
    character(len=128) :: arg
    n = iargc()

    ! if (n < 4) then
    !    print*, 'Sort the galaxies between elliptics and spirals'
    !    print*, ''
    !    print*, 'Usage (all lists should come in binary format):'
    !    print*, '\t -fga File to pick galaxies from'
    !    print*, '\t -fdm File to pick dark matter halos from'
    !    print*, '\t -fas File to get associations between dark matter halos and galaxies'
    !    !TODO: write usage
    ! end if

    !-------------------------------------
    ! default values
    !-------------------------------------
    gal_list_filename = "lists/list_kingal_00782.dat"
    dm_halo_list_filename = "lists/list_halo.dat.bin"
    associations_filename = "lists/associated_halogal_782.dat.bin"
    ramses_output_start = "/data52/Horizon-AGN/OUTPUT_DIR/output_00032"
    ramses_output_end = "/data52/Horizon-AGN/OUTPUT_DIR/output_00761"
    brick_file = "/data40b/Horizon-AGN/TREE_DM_celldx2kpc_Rmax_guess9/tree_bricks032"
    info_file_end = '/data52/Horizon-AGN/OUTPUT_DIR/output_00761/info_00761.txt'
    info_file_start = '/data52/Horizon-AGN/OUTPUT_DIR/output_00032/info_00032.txt'
    output_dir = '/data74/cadiou/halo_files/'

    do i = 1, n, 2
       call getarg(i, opt)
       if (i == n) then
          print*, '("option ",a2," has no argument")', opt
          stop
       end if

       call getarg(i+1, arg)
       select case(opt)
       case ('-fga')
          gal_list_filename = trim(arg)
       case ('-fdm')
          dm_halo_list_filename = trim(arg)
       case('-fas')
          associations_filename = trim(arg)
       case('-ods') ! Ramses output start
          ramses_output_start = trim(arg)
       case('-ode') ! Ramses output start
          ramses_output_end = trim(arg)
       case('-bri')
          brick_file = trim(arg)
       case('-ifs')
          info_file_start = trim(arg)
       case('-ife')
          info_file_end = trim(arg)
       case ('-out')
          outfile = trim(arg)
       case ('-fin')
          fileinfo = trim(arg)
       case ('-v')
          read (arg,*) verbose
       case default
          print '("unknown option ",a2," ignored")', opt
       end select
    end do

  end subroutine read_params

end program sort_galaxy


text:
      SUBROUTINE DBMALB ( LENREQ, INDEX )
C********************************************************************
C     DBMALB - ALLOCATES A MEMORY BLOCK OF LENGTH "LENREQ"
C              FROM THE FREE CHAIN AND RETURNS THE
C              POINTER IN MEMORY FOR THE BLOCK IN "INDEX" (RELATIVE TO.
C              /DBM/.
C
C     EACH FREE BLOCK IN MEMORY HAS THE FOLLOWING FORMAT:
C         WORD 1  POINTER TO PREVIOUS FREE BLOCK (=0, IF FIRST)
C         WORD 2  POINTER TO NEXT FREE BLOCK (=0, IF END OF CHAIN)
C         WORD 3  NUMBER OF WORDS AVAILABLE IN THIS FREE BLOCK
C
C     NOTE:  IDBFRE POINTS TO THE FIRST FREE BLOCK OF THE CHAIN
C********************************************************************
      INCLUDE 'DSIOF.COM'
      INCLUDE 'ZZZZZZ.COM'
      NEXT  = IDBFRE
      IF ( IDBFRE .EQ. 0 ) GO TO 701
C  OBTAIN THE LENGTH OF THE FREE BLOCK
10    LENAVL = MEM( NEXT + 2 )
      IF ( LENAVL .GE. LENREQ ) GO TO 100
C  MEMORY NOT AVAILABLE IN THIS BLOCK, CHECK FOR OTHER BLOCKS
      NEXT  = MEM( NEXT + 1 )
C  IF NO MORE FREE BLOCKS, RETURN WITH INDEX SET TO -1
      IF ( NEXT .EQ. 0 ) GO TO 701
      GO TO 10
C  RETURN POINTER FOR THIS BLOCK
100   INDEX = NEXT
      IF ( LENREQ .NE. LENAVL ) GO TO 200
C  COME HERE WHEN REQUESTED BLOCK SAME SIZE AS FREE BLOCK
      NEXT  = MEM( NEXT+1 )
      IPREV = MEM(INDEX   )
      IF ( IPREV .EQ. 0 ) GO TO 110
      IF ( NEXT  .EQ. 0 ) GO TO 120
C  CONNECT THE PREVIOUS FREE BLOCK WITH THE NEXT FREE BLOCK
      MEM( IPREV+1 ) = NEXT
      MEM( NEXT    ) = IPREV
      GO TO 700
110   IF ( NEXT .EQ. 0 ) GO TO 130
C  NO PREVIOUS BLOCK, SET IDBFRE TO POINT TO NEW FIRST FREE BLOCK
      IDBFRE = NEXT
      MEM( NEXT ) = 0
      GO TO 700
C  PREVIOUS BLOCK EXITS BUT BLOCK ALLOCATED WAS LAST IN CHAIN
120   MEM( IPREV+1 ) = 0
      GO TO 700
C  NO MORE FREE BLOCKS EXIST, SET IDBFRE TO ZERO
130   IDBFRE = 0
      GO TO 700
C  COME HERE WHEN FREE BLOCK HAS MORE SPACE THEN REQUESTED
200   NEWIND = INDEX + LENREQ + 4
      IPREV  = MEM( INDEX  )
      NEXT   = MEM( INDEX+1)
C  CHECK TO DETERMINE IF ANY SPACE REMAINS
      IF ( ( LENAVL-LENREQ-4 ) .LE. 0 ) GO TO 240
C  RECOMPUTE FREE SPACE AND SET UP CHAIN WORDS
      MEM( NEWIND+2 ) = LENAVL - LENREQ - 4
      IF ( IPREV .EQ. 0 ) GO TO 210
      IF ( NEXT  .EQ. 0 ) GO TO 220
C  CONNECT TO PREVIOUS AND NEXT FREE BLOCK
      MEM( NEWIND  ) = IPREV
      MEM( NEWIND+1) = NEXT
      MEM( IPREV+1 ) = NEWIND
      MEM( NEXT    ) = NEWIND
      GO TO 700
210   IF ( NEXT .EQ. 0 ) GO TO 230
C  NO PREVIOUS BLOCK, NEWLY CREATED BLOCK BECOMES THE FIRST FREE BLOCK
      IDBFRE         = NEWIND
      MEM( NEWIND  ) = 0
      MEM( NEWIND+1) = NEXT
      MEM( NEXT    ) = NEWIND
      GO TO 700
C  PREVIOUS BLOCK EXISTS BUT THE NEWLY CREATED BLOCK IS LAST
220   MEM( IPREV+1 ) = NEWIND
      MEM( NEWIND  ) = IPREV
      MEM( NEWIND+1) = 0
      GO TO 700
C  NEW BLOCK IS THE ONLY FREE BLOCK
230   IDBFRE  = NEWIND
      MEM( NEWIND  ) = 0
      MEM( NEWIND+1) = 0
      GO TO 700
C  FREE CHAIN IS EXHAUSTED
240   IDBFRE = 0
701   INDEX = -1
      RETURN
700   CONTINUE
      RETURN
      END

ext: f


text:
!
! Copyright (C) 2001-2009 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
subroutine pz_polarized (rs, ec, vc)
  !-----------------------------------------------------------------------
  !     J.P. Perdew and A. Zunger, PRB 23, 5048 (1981)
  !     spin-polarized energy and potential
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rs, ec, vc
  real(DP) :: a, b, c, d, gc, b1, b2
  parameter (a = 0.01555d0, b = - 0.0269d0, c = 0.0007d0, d = &
       - 0.0048d0, gc = - 0.0843d0, b1 = 1.3981d0, b2 = 0.2611d0)
  real(DP) :: lnrs, rs12, ox, dox
  REAL(DP), PARAMETER :: xcprefact = 0.022575584d0, pi34 = 0.6203504908994d0
  ! REAL(DP) :: betha, etha, csi, prefact
  !
  if (rs.lt.1.0d0) then
     ! high density formula
     lnrs = log (rs)
     ec = a * lnrs + b + c * rs * lnrs + d * rs
     vc = a * lnrs + (b - a / 3.d0) + 2.d0 / 3.d0 * c * rs * lnrs + &
          (2.d0 * d-c) / 3.d0 * rs
  else
     ! interpolation formula
     rs12 = sqrt (rs)
     ox = 1.d0 + b1 * rs12 + b2 * rs
     dox = 1.d0 + 7.d0 / 6.d0 * b1 * rs12 + 4.d0 / 3.d0 * b2 * rs
     ec = gc / ox
     vc = ec * dox / ox
  endif
  !
!  IF ( lxc_rel ) THEN
!     betha = prefact * pi34 / rs
!     etha = DSQRT( 1 + betha**2 )
!     csi = betha + etha
!     prefact = 1.0D0 - (3.0D0/2.0D0) * ( (betha*etha - log(csi))/betha**2 )**2
!     ec = ec * prefact
!     vc = vc * prefact
!  ENDIF
  return
end subroutine pz_polarized
!
!-----------------------------------------------------------------------
subroutine pz_spin (rs, zeta, ec, vcup, vcdw)
  !-----------------------------------------------------------------------
  !     J.P. Perdew and Y. Wang, PRB 45, 13244 (1992)
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rs, zeta, ec, vcup, vcdw
  !
  real(DP) :: ecu, vcu, ecp, vcp, fz, dfz
  real(DP) :: p43, third
  parameter (p43 = 4.0d0 / 3.d0, third = 1.d0 / 3.d0)
  !
  ! unpolarized part (Perdew-Zunger formula)
  call pz (rs, 1, ecu, vcu)
  ! polarization contribution
  call pz_polarized (rs, ecp, vcp)
  !
  fz = ( (1.0d0 + zeta) **p43 + (1.d0 - zeta) **p43 - 2.d0) / &
       (2.d0**p43 - 2.d0)
  dfz = p43 * ( (1.0d0 + zeta) **third- (1.d0 - zeta) **third) &
       / (2.d0**p43 - 2.d0)
  !
  ec = ecu + fz * (ecp - ecu)
  vcup = vcu + fz * (vcp - vcu) + (ecp - ecu) * dfz * (1.d0 - zeta)
  vcdw = vcu + fz * (vcp - vcu) + (ecp - ecu) * dfz * ( - 1.d0 - &
       zeta)
  !
  return
end subroutine pz_spin
!
!---------

!-----------------------------------------------------------------------
subroutine pw_spin (rs, zeta, ec, vcup, vcdw)
  !-----------------------------------------------------------------------
  !     J.P. Perdew and Y. Wang, PRB 45, 13244 (1992)
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rs, zeta, ec, vcup, vcdw
  ! xc parameters, unpolarised
  real(DP) :: a, a1, b1, b2, b3, b4, c0, c1, c2, c3, d0, d1
  parameter (a = 0.031091d0, a1 = 0.21370d0, b1 = 7.5957d0, b2 = &
       3.5876d0, b3 = 1.6382d0, b4 = 0.49294d0, c0 = a, c1 = 0.046644d0, &
       c2 = 0.00664d0, c3 = 0.01043d0, d0 = 0.4335d0, d1 = 1.4408d0)
  ! xc parameters, polarised
  real(DP) :: ap, a1p, b1p, b2p, b3p, b4p, c0p, c1p, c2p, c3p, d0p, &
       d1p
  parameter (ap = 0.015545d0, a1p = 0.20548d0, b1p = 14.1189d0, b2p &
       = 6.1977d0, b3p = 3.3662d0, b4p = 0.62517d0, c0p = ap, c1p = &
       0.025599d0, c2p = 0.00319d0, c3p = 0.00384d0, d0p = 0.3287d0, d1p &
       = 1.7697d0)
  ! xc parameters, antiferro
  real(DP) :: aa, a1a, b1a, b2a, b3a, b4a, c0a, c1a, c2a, c3a, d0a, &
       d1a
  parameter (aa = 0.016887d0, a1a = 0.11125d0, b1a = 10.357d0, b2a = &
       3.6231d0, b3a = 0.88026d0, b4a = 0.49671d0, c0a = aa, c1a = &
       0.035475d0, c2a = 0.00188d0, c3a = 0.00521d0, d0a = 0.2240d0, d1a &
       = 0.3969d0)
  real(DP) :: fz0
  parameter (fz0 = 1.709921d0)
  real(DP) :: rs12, rs32, rs2, zeta2, zeta3, zeta4, fz, dfz
  real(DP) :: om, dom, olog, epwc, vpwc
  real(DP) :: omp, domp, ologp, epwcp, vpwcp
  real(DP) :: oma, doma, ologa, alpha, vpwca
  !
  !     if(rs.lt.0.5d0) then
  ! high density formula (not implemented)
  !
  !     else if(rs.gt.100.d0) then
  ! low density formula  (not implemented)
  !
  !     else
  ! interpolation formula
  zeta2 = zeta * zeta
  zeta3 = zeta2 * zeta
  zeta4 = zeta3 * zeta
  rs12 = sqrt (rs)
  rs32 = rs * rs12
  rs2 = rs**2
  ! unpolarised
  om = 2.d0 * a * (b1 * rs12 + b2 * rs + b3 * rs32 + b4 * rs2)
  dom = 2.d0 * a * (0.5d0 * b1 * rs12 + b2 * rs + 1.5d0 * b3 * rs32 &
       + 2.d0 * b4 * rs2)
  olog = log (1.d0 + 1.0d0 / om)
  epwc = - 2.d0 * a * (1.d0 + a1 * rs) * olog
  vpwc = - 2.d0 * a * (1.d0 + 2.d0 / 3.d0 * a1 * rs) * olog - 2.d0 / &
       3.d0 * a * (1.d0 + a1 * rs) * dom / (om * (om + 1.d0) )
  ! polarized
  omp = 2.d0 * ap * (b1p * rs12 + b2p * rs + b3p * rs32 + b4p * rs2)
  domp = 2.d0 * ap * (0.5d0 * b1p * rs12 + b2p * rs + 1.5d0 * b3p * &
       rs32 + 2.d0 * b4p * rs2)
  ologp = log (1.d0 + 1.0d0 / omp)
  epwcp = - 2.d0 * ap * (1.d0 + a1p * rs) * ologp
  vpwcp = - 2.d0 * ap * (1.d0 + 2.d0 / 3.d0 * a1p * rs) * ologp - &
       2.d0 / 3.d0 * ap * (1.d0 + a1p * rs) * domp / (omp * (omp + 1.d0) &
       )
  ! antiferro
  oma = 2.d0 * aa * (b1a * rs12 + b2a * rs + b3a * rs32 + b4a * rs2)
  doma = 2.d0 * aa * (0.5d0 * b1a * rs12 + b2a * rs + 1.5d0 * b3a * &
       rs32 + 2.d0 * b4a * rs2)
  ologa = log (1.d0 + 1.0d0 / oma)
  alpha = 2.d0 * aa * (1.d0 + a1a * rs) * ologa
  vpwca = + 2.d0 * aa * (1.d0 + 2.d0 / 3.d0 * a1a * rs) * ologa + &
       2.d0 / 3.d0 * aa * (1.d0 + a1a * rs) * doma / (oma * (oma + 1.d0) &
       )
  !
  fz = ( (1.d0 + zeta) ** (4.d0 / 3.d0) + (1.d0 - zeta) ** (4.d0 / &
       3.d0) - 2.d0) / (2.d0** (4.d0 / 3.d0) - 2.d0)
  dfz = ( (1.d0 + zeta) ** (1.d0 / 3.d0) - (1.d0 - zeta) ** (1.d0 / &
       3.d0) ) * 4.d0 / (3.d0 * (2.d0** (4.d0 / 3.d0) - 2.d0) )
  !
  ec = epwc + alpha * fz * (1.d0 - zeta4) / fz0 + (epwcp - epwc) &
       * fz * zeta4
  !
  vcup = vpwc + vpwca * fz * (1.d0 - zeta4) / fz0 + (vpwcp - vpwc) &
       * fz * zeta4 + (alpha / fz0 * (dfz * (1.d0 - zeta4) - 4.d0 * fz * &
       zeta3) + (epwcp - epwc) * (dfz * zeta4 + 4.d0 * fz * zeta3) ) &
       * (1.d0 - zeta)

  vcdw = vpwc + vpwca * fz * (1.d0 - zeta4) / fz0 + (vpwcp - vpwc) &
       * fz * zeta4 - (alpha / fz0 * (dfz * (1.d0 - zeta4) - 4.d0 * fz * &
       zeta3) + (epwcp - epwc) * (dfz * zeta4 + 4.d0 * fz * zeta3) ) &
       * (1.d0 + zeta)
  !      endif
  !
  return
end subroutine pw_spin
!
!---------------------------------------------------------------
subroutine pbec_spin (rho, zeta, grho, iflag, sc, v1cup, v1cdw, v2c)
  !---------------------------------------------------------------
  !
  ! PBE correlation (without LDA part) - spin-polarized
  ! iflag = 1: J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996).
  ! iflag = 2: J.P.Perdew et al., PRL 100, 136406 (2008)
  !
  USE kinds, ONLY : DP
  implicit none
  integer, intent(in) :: iflag
  real(DP) :: rho, zeta, grho, sc, v1cup, v1cdw, v2c
  real(DP) :: ga, be(2)
  parameter (ga = 0.031091d0)
  data be / 0.066725d0 ,  0.046d0 /
  real(DP) :: third, pi34, xkf, xks
  parameter (third = 1.d0 / 3.d0, pi34 = 0.6203504908994d0)
  parameter (xkf = 1.919158292677513d0, xks = 1.128379167095513d0)
  ! pi34=(3/4pi)^(1/3), xkf=(9 pi/4)^(1/3), xks= sqrt(4/pi)
  real(DP) :: kf, ks, rs, ec, vcup, vcdw, t, expe, af, y, xy, qy, &
       s1, h0, ddh0
  real(DP) :: fz, fz2, fz3, fz4, dfz, bfup, bfdw, dh0up, dh0dw, &
       dh0zup, dh0zdw
  !
  rs = pi34 / rho**third
  call pw_spin (rs, zeta, ec, vcup, vcdw)
  kf = xkf / rs
  ks = xks * sqrt (kf)
  fz = 0.5d0 * ( (1.d0 + zeta) ** (2.d0 / 3.d0) + (1.d0 - zeta) ** ( &
       2.d0 / 3.d0) )
  fz2 = fz * fz
  fz3 = fz2 * fz
  fz4 = fz3 * fz
  dfz = ( (1.d0 + zeta) ** ( - 1.d0 / 3.d0) - (1.d0 - zeta) ** ( - &
       1.d0 / 3.d0) ) / 3.d0
  t = sqrt (grho) / (2.d0 * fz * ks * rho)
  expe = exp ( - ec / (fz3 * ga) )
  af = be(iflag) / ga * (1.d0 / (expe-1.d0) )
  bfup = expe * (vcup - ec) / fz3
  bfdw = expe * (vcdw - ec) / fz3
  y = af * t * t
  xy = (1.d0 + y) / (1.d0 + y + y * y)
  qy = y * y * (2.d0 + y) / (1.d0 + y + y * y) **2
  s1 = 1.d0 + be(iflag) / ga * t * t * xy
  h0 = fz3 * ga * log (s1)
  dh0up = be(iflag) * t * t * fz3 / s1 * ( - 7.d0 / 3.d0 * xy - qy * &
       (af * bfup / be(iflag)-7.d0 / 3.d0) )
  dh0dw = be(iflag) * t * t * fz3 / s1 * ( - 7.d0 / 3.d0 * xy - qy * &
       (af * bfdw / be(iflag)-7.d0 / 3.d0) )
  dh0zup = (3.d0 * h0 / fz - be(iflag) * t * t * fz2 / s1 * (2.d0 * xy - &
  qy * (3.d0 * af * expe * ec / fz3 / be(iflag)+2.d0) ) ) * dfz * (1.d0 - zeta)
  dh0zdw = - (3.d0 * h0 / fz - be(iflag) * t * t * fz2 / s1 * (2.d0 * xy - &
  qy * (3.d0 * af * expe * ec / fz3 / be(iflag)+2.d0) ) ) * dfz * (1.d0 + zeta)

  ddh0 = be(iflag) * fz / (2.d0 * ks * ks * rho) * (xy - qy) / s1
  sc = rho * h0
  v1cup = h0 + dh0up + dh0zup
  v1cdw = h0 + dh0dw + dh0zdw
  v2c = ddh0
  return
end subroutine pbec_spin
!
!-----------------------------------------------------------------------
subroutine slater_spin (rho, zeta, ex, vxup, vxdw)
  !-----------------------------------------------------------------------
  !     Slater exchange with alpha=2/3, spin-polarized case
  !
  USE kinds, ONLY : DP
  implicit none
  real(DP) :: rho, zeta, ex, vxup, vxdw
  real(DP) :: f, alpha, third, p43
  parameter (f = - 1.10783814957303361d0, alpha = 2.0d0 / 3.0d0)
  ! f = -9/8*(3/pi)^(1/3)
  parameter (third = 1.d0 / 3.d0, p43 = 4.d0 / 3.d0)
  real(DP) :: exup, exdw, rho13
  !
  rho13 = ( (1.d0 + zeta) * rho) **third
  exup = f * alpha * rho13
  vxup = p43 * f * alpha * rho13
  rho13 = ( (1.d0 - zeta) * rho) **third
  exdw = f * alpha * rho13
  vxdw = p43 * f * alpha * rho13
  ex = 0.5d0 * ( (1.d0 + zeta) * exup + (1.d0 - zeta) * exdw)
  !
  return
end subroutine slater_spin

!
!-----------------------------------------------------------------------
function dpz_polarized (rs, iflg)
  !-----------------------------------------------------------------------
  !  derivative of the correlation potential with respect to local density
  !  Perdew and Zunger parameterization of the Ceperley-Alder functional
  !  spin-polarized case
  !
  USE kinds, only : DP
  USE constants, ONLY : pi, fpi
  !
  implicit none
  !
  real(DP), intent (in) :: rs
  integer, intent(in) :: iflg
  real(DP) :: dpz_polarized
  !
  !  local variables
  !  a,b,c,d,gc,b1,b2 are the parameters defining the functional
  !
  real(DP), parameter :: a = 0.01555d0, b = -0.0269d0, c = 0.0007d0, &
       d = -0.0048d0, gc = -0.0843d0, b1 = 1.3981d0, b2 = 0.2611d0,&
       a1 = 7.0d0 * b1 / 6.d0, a2 = 4.d0 * b2 / 3.d0
  real(DP) :: x, den, dmx, dmrs
  !
  !
  if (iflg == 1) then
     dmrs = a / rs + 2.d0 / 3.d0 * c * (log (rs) + 1.d0) + &
          (2.d0 * d-c) / 3.d0
  else
     x = sqrt (rs)
     den = 1.d0 + x * (b1 + x * b2)
     dmx = gc * ( (a1 + 2.d0 * a2 * x) * den - 2.d0 * (b1 + 2.d0 * &
          b2 * x) * (1.d0 + x * (a1 + x * a2) ) ) / den**3
     dmrs = 0.5d0 * dmx / x
  endif
  !
  dpz_polarized = - fpi * rs**4.d0 / 9.d0 * dmrs
  return
  !
end function dpz_polarized

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
      SUBROUTINE SAOA (T,ANS)
C
C     SAOA COMPUTES THE INTEGRAND FOR EACH OF THE 6
C     SOMMERFELD INTEGRALS FOR SOURCE AND OBSERVER ABOVE GROUND
C
      COMPLEX ANS,XL,DXL,CGAM1,CGAM2,B0,B0P,COM,CK1,CK1SQ,CKSM,CT1,CT2,C
     1T3,DGAM,DEN1,DEN2
      COMMON /EVLCOM/ CKSM,CT1,CT2,CT3,CK1,CK1SQ,CK2,CK2SQ,TKMAG,TSMAG,C
     1K1R,ZPH,RHO,JH
      DIMENSION ANS(6)
      CALL LAMBDA (T,XL,DXL)
      IF (JH.GT.0) GO TO 1
C     BESSEL FUNCTION FORM
      CALL BESSEL (XL*RHO,B0,B0P)
      B0=2.*B0
      B0P=2.*B0P
      CGAM1=CSQRT(XL*XL-CK1SQ)
      CGAM2=CSQRT(XL*XL-CK2SQ)
      IF (REAL(CGAM1).EQ.0.) CGAM1=DCMPLX(0.D0,-ABS(AIMAG(CGAM1)))
      IF (REAL(CGAM2).EQ.0.) CGAM2=DCMPLX(0.D0,-ABS(AIMAG(CGAM2)))
      GO TO 2
C     HANKEL FUNCTION FORM
1     CALL HANKEL (XL*RHO,B0,B0P)
      COM=XL-CK1
      CGAM1=CSQRT(XL+CK1)*CSQRT(COM)
      IF (REAL(COM).LT.0..AND.AIMAG(COM).GE.0.) CGAM1=-CGAM1
      COM=XL-CK2
      CGAM2=CSQRT(XL+CK2)*CSQRT(COM)
      IF (REAL(COM).LT.0..AND.AIMAG(COM).GE.0.) CGAM2=-CGAM2
2     XLR=XL*CONJG(XL)
      IF (XLR.LT.TSMAG) GO TO 3
      IF (AIMAG(XL).LT.0.) GO TO 4
      XLR=REAL(XL)
      IF (XLR.LT.CK2) GO TO 5
      IF (XLR.GT.CK1R) GO TO 4
3     DGAM=CGAM2-CGAM1
      GO TO 7
4     SIGN=1.
      GO TO 6
5     SIGN=-1.
6     DGAM=1./(XL*XL)
      DGAM=SIGN*((CT3*DGAM+CT2)*DGAM+CT1)/XL
7     DEN2=CKSM*DGAM/(CGAM2*(CK1SQ*CGAM2+CK2SQ*CGAM1))
      DEN1=1./(CGAM1+CGAM2)-CKSM/CGAM2
      COM=DXL*XL*CEXP(-CGAM2*ZPH)
      ANS(6)=COM*B0*DEN1/CK1
      COM=COM*DEN2
      IF (RHO.EQ.0.) GO TO 8
      B0P=B0P/RHO
      ANS(1)=-COM*XL*(B0P+B0*XL)
      ANS(4)=COM*XL*B0P
      GO TO 9
8     ANS(1)=-COM*XL*XL*.5
      ANS(4)=ANS(1)
9     ANS(2)=COM*CGAM2*CGAM2*B0
      ANS(3)=-ANS(4)*CGAM2*RHO
      ANS(5)=COM*B0
      RETURN
      END

max_stars_repo_name: ViperScience/ViperNEC

max_issues_repo_name: ViperScience/ViperNEC

max_forks_repo_name: ViperScience/ViperNEC


module psi_d_serial_mod
  use psb_const_mod, only : psb_ipk_, psb_dpk_

  interface psb_gelp
    ! 2-D version
    subroutine psb_dgelp(trans,iperm,x,info)
      import :: psb_ipk_, psb_dpk_
      implicit none
      real(psb_dpk_), intent(inout)     ::  x(:,:)
      integer(psb_ipk_), intent(in)      ::  iperm(:)
      integer(psb_ipk_), intent(out)     ::  info
      character, intent(in)              :: trans
    end subroutine psb_dgelp
    subroutine psb_dgelpv(trans,iperm,x,info)
      import :: psb_ipk_, psb_dpk_
      implicit none
      real(psb_dpk_), intent(inout)     ::  x(:)
      integer(psb_ipk_), intent(in)      ::  iperm(:)
      integer(psb_ipk_), intent(out)     ::  info
      character, intent(in)              :: trans
    end subroutine psb_dgelpv
  end interface psb_gelp

  interface psb_geaxpby
    subroutine psi_daxpby(m,n,alpha, x, beta, y, info)
      import :: psb_ipk_, psb_dpk_
      implicit none
      integer(psb_ipk_), intent(in)      :: m, n
      real(psb_dpk_), intent (in)       ::  x(:,:)
      real(psb_dpk_), intent (inout)    ::  y(:,:)
      real(psb_dpk_), intent (in)       ::  alpha, beta
      integer(psb_ipk_), intent(out)     :: info
    end subroutine psi_daxpby
    subroutine psi_daxpbyv(m,alpha, x, beta, y, info)
      import :: psb_ipk_, psb_dpk_
      implicit none
      integer(psb_ipk_), intent(in)      :: m
      real(psb_dpk_), intent (in)       ::  x(:)
      real(psb_dpk_), intent (inout)    ::  y(:)
      real(psb_dpk_), intent (in)       :: alpha, beta
      integer(psb_ipk_), intent(out)     :: info
    end subroutine psi_daxpbyv
  end interface psb_geaxpby

  interface psi_gth
    subroutine psi_dgthmv(n,k,idx,alpha,x,beta,y)
      import :: psb_ipk_, psb_dpk_
      implicit none
      integer(psb_ipk_) :: n, k, idx(:)
      real(psb_dpk_) :: x(:,:), y(:),alpha,beta
    end subroutine psi_dgthmv
    subroutine psi_dgthv(n,idx,alpha,x,beta,y)
      import :: psb_ipk_, psb_dpk_
      implicit none
      integer(psb_ipk_) :: n, idx(:)
      real(psb_dpk_) :: x(:), y(:),alpha,beta
    end subroutine psi_dgthv
    subroutine psi_dgthzmv(n,k,idx,x,y)
      import :: psb_ipk_, psb_dpk_
      implicit none
      integer(psb_ipk_) :: n, k, idx(:)
      real(psb_dpk_) :: x(:,:), y(:)

    end subroutine psi_dgthzmv
    subroutine psi_dgthzmm(n,k,idx,x,y)
      import :: psb_ipk_, psb_dpk_
      implicit none
      integer(psb_ipk_) :: n, k, idx(:)
      real(psb_dpk_) :: x(:,:), y(:,:)

    end subroutine psi_dgthzmm
    subroutine psi_dgthzv(n,idx,x,y)
      import :: psb_ipk_, psb_dpk_
      implicit none
      integer(psb_ipk_) :: n, idx(:)
      real(psb_dpk_) :: x(:), y(:)
    end subroutine psi_dgthzv
  end interface psi_gth

  interface psi_sct
    subroutine psi_dsctmm(n,k,idx,x,beta,y)
      import :: psb_ipk_, psb_dpk_
      implicit none
      integer(psb_ipk_) :: n, k, idx(:)
      real(psb_dpk_) :: beta, x(:,:), y(:,:)
    end subroutine psi_dsctmm
    subroutine psi_dsctmv(n,k,idx,x,beta,y)
      import :: psb_ipk_, psb_dpk_
      implicit none
      integer(psb_ipk_) :: n, k, idx(:)
      real(psb_dpk_) :: beta, x(:), y(:,:)
    end subroutine psi_dsctmv
    subroutine psi_dsctv(n,idx,x,beta,y)
      import :: psb_ipk_, psb_dpk_
      implicit none

      integer(psb_ipk_) :: n, idx(:)
      real(psb_dpk_) :: beta, x(:), y(:)
    end subroutine psi_dsctv
  end interface psi_sct

end module psi_d_serial_mod


text:
C **********************************************************************
C *                                                                    *
C * LOWER - CONVERT C TO LOWER-CASE, IF IT IS UPPER-CASE               *
C *                                                                    *
C **********************************************************************
C
      INTEGER FUNCTION LOWER(C)
      INTEGER C
C
      LOWER = C
      IF (C .GE. 65 .AND. C .LE. 90) LOWER = LOWER - 65 + 97
      RETURN
      END

ext: for


module SphereBVEModule
!> @file SphereBVE.f90
!> Data structure for representing solutions of the Barotropic Vorticity Equation (BVE) on the surface of a rotating sphere
!> @author Peter Bosler, Sandia National Laboratories Center for Computing Research
!>
!>
!> @defgroup SphereBVE SphereBVE
!> @brief Data structure for representing solutions of the Barotropic Vorticity Equation (BVE) on the surface of a rotating sphere.
!>
!> Combines a @ref PolyMesh2d with the data @ref Field objects for the BVE.
!> Allows users to optionally add @ref Field objects for passive tracers.
!>
!> In addition to the variables carried by the base @ref Particles object (e.g., physical and Lagrangian coordinates),
!> the BVEMesh data type adds fields for relative vorticity @f$ \zeta(x,y,z,t) @f$ and the materially conserved absolute vorticity @f$ \omega(x_0,y_0,z_0) @f$.
!> These variables are related by the equation
!> @f[
!> 		\omega = \zeta + f,
!> @f]
!> where @f$ f = \frac{2\Omega}{a}z(t) @f$ is the Coriolis parameter, @f$ \Omega @f$ is the constant angular velocity of the sphere about the z-axis,
!> and @f$ a @f$ is the radius of the sphere.
!>
!> Two stream functions are defined on the particles, @f$ \psi_r(x,y,z,t) @f$ is defined with respect to the rotating frame
!> and @f$ \psi_a(x,y,z,t) @f$ is defined with respect to the inertial frame.
!> Both are represented (and computed) by convolution with the Green's function kernel,
!> @f{align*}{
!> 		\psi_r(\vec{x},t) &= \frac{1}{4\pi a}\int_{S} \log( 1 - \vec{x} \cdot \vec{\tilde{x}})\zeta(\vec{\tilde{x}})\,dA, \\
!> 	    \psi_a(\vec{x},t) &= \frac{1}{4\pi a}\int_{S} \log( 1 - \vec{x} \cdot \vec{\tilde{x}})\omega(\vec{\tilde{x}})\,dA.
!> @f}
!>
!> The particles also have a @ref Field to track their velocity, which is given by convolution with the Biot-Savart kernel,
!> @f[
!> 		\vec{u}(\vec{x},t) = -\frac{1}{4\pi a}\int_{S} \frac{ \vec{x} \times \vec{\tilde{x}} }{a^2 - \vec{x}\cdot\vec{\tilde{x}}} \zeta(\vec{\tilde{x}})\,dA.
!> @f]
!>
!>	For a more detailed discussion, see @n
!> * P. Bosler, L. Wang, C. Jablonowski, and R. Krasny, A Lagrangian particle/panel method for the barotropic vorticity
!> 	equations on a rotating sphere, _Fluid Dynamics Research_ 46 (2014).  DOI: 10.1088/0169-5983/46/3/031406.
!>
!>
!> @{
use NumberKindsModule
use OutputWriterModule
use LoggerModule
use ParticlesModule
use EdgesModule
use FacesModule
use PolyMesh2dModule
use FieldModule
use MPISetupModule
use SphereGeomModule

implicit none

include 'mpif.h'

private
!
!----------------
! Module types and public declarations
!----------------
!
public BVEMesh, New, Delete, Copy
public AddTracers
public LogStats
public OutputToVTK
public SetInitialVorticityOnMesh, SetVelocityOnMesh, SetStreamFunctionsOnMesh
public TotalKE, TotalEnstrophy
public MaxCirculationMagnitudePerFace

!> @brief Defines a spherical mesh and the appropriate variables for solving the Barotropic Vorticity Equation (BVE) on a rotating sphere.
type BVEMesh
	type(PolyMesh2d) :: mesh !< @ref PolyMesh2d for spatial discretization
	type(Field) :: absVort !< scalar @ref Field
	type(Field) :: relVort !< scalar @ref Field
	type(Field) :: relStream !< scalar @ref Field
	type(Field) :: absStream !< scalar @ref Field
	type(Field) :: velocity !< vector @ref Field
	type(Field), dimension(:), allocatable :: tracers !< allocatable array of scalar and vector @ref Field objects
	real(kreal) :: radius = 1.0_kreal !< radius of sphere
	real(kreal) :: rotationRate = 0.0_kreal !< background rotation rate (angular velocity)
	type(MPISetup) :: mpiParticles !< @ref MPISetup to distribute all particles over the available MPI ranks

	contains
		final :: deletePrivate
end type

!
!----------------
! Module interfaces
!----------------
!

interface New
	module procedure newPrivate
end interface

interface Delete
	module procedure deletePrivate
end interface

interface Copy
	module procedure copyPrivate
end interface

interface LogStats
	module procedure logStatsPrivate
end interface

interface SetVelocityOnMesh
	module procedure setVelocityFromFunction
	module procedure setVelocityFromVorticity
end interface

interface OutputToVTK
	module procedure outputVTKPrivate
end interface

!
!----------------
! Logging
!----------------
!
logical(klog), save :: logInit = .FALSE.
type(Logger) :: log
character(len=28), save :: logKey = 'BVE'
integer(kint), parameter :: logLevel = DEBUG_LOGGING_LEVEL

contains
!
!----------------
! public methods
!----------------
!

!> @brief Allocates memory and initializes a spherical mesh for BVE applications.
!>
!> @param[out] self new BVE mesh
!> @param[in] meshSeed mesh seed integer, as defined in @ref NumberKinds
!> @param[in] initNest initial level of uniform refinement
!> @param[in] maxNest maximum level of uniform refinment
!> @param[in] amrLimit maximum number of times an individual face can be divided beyond initNest
!> @param[in] sphereRadius radius of sphere
!> @param[in] rotationRate background angular velocity of sphere
subroutine newPrivate( self, meshSeed, initNest, maxNest, amrLimit, sphereRadius, rotationRate )
	type(BVEMesh), intent(out) :: self
	integer(kint), intent(in) :: meshSeed
	integer(kint), intent(in) :: initNest
	integer(kint), intent(in) :: maxNest
	integer(kint), intent(in) :: amrLimit
	real(kreal), intent(in) :: sphereRadius
	real(kreal), intent(in) :: rotationRate

	if ( .NOT. logInit) call InitLogger(log, procRank)

	if ( meshSeed /= ICOS_TRI_SPHERE_SEED .AND. meshSeed /= CUBED_SPHERE_SEED ) then
		call LogMessage(log, ERROR_LOGGING_LEVEL, trim(logKey), "new BVEMesh ERROR : invalid meshSeed.")
		return
	endif

	call New(self%mesh, meshSeed, initNest, maxNest, amrLimit, sphereRadius)

	call New(self%absVort, 1, self%mesh%particles%N_Max, "absVort", "1/s")
	call New(self%relVort, 1, self%mesh%particles%N_Max, "relVort", "1/s")
	call New(self%absStream, 1, self%mesh%particles%N_Max, "absStreamFn", "m^2/s")
	call New(self%relStream, 1, self%mesh%particles%N_Max, "relStreamFn", "m^2/s")
	call New(self%velocity, 3, self%mesh%particles%N_Max, "velocity", "m/s")
	call New(self%mpiParticles, self%mesh%particles%N, numProcs)

	self%radius = sphereRadius
	self%rotationRate = rotationRate
end subroutine

!> @brief Deletes and frees memory associated with a spherical BVE Mesh
!> @param[inout] self Target BVE Mesh
subroutine deletePrivate( self )
	type(BVEMesh), intent(inout) :: self
	integer(kint) :: i

	call Delete(self%mpiParticles)
	call Delete(self%velocity)
	call Delete(self%relStream)
	call Delete(self%absStream)
	call Delete(self%absVort)
	call Delete(self%relVort)
	call Delete(self%mesh)
	if ( allocated(self%tracers) ) then
		do i = 1, size(self%tracers)
			call Delete(self%tracers(i))
		enddo
		deallocate(self%tracers)
	endif
end subroutine

!> @brief Performs a deep copy of a BVE mesh and all of its variables.
!> The target mesh must have been allocated prior to calling this subroutine.
!>
!> @param[inout] self Target BVE mesh
!> @param[in] other Source BVE Mesh
subroutine copyPrivate(self, other)
	type(BVEMesh), intent(inout) :: self
	type(BVEMesh), intent(in) :: other
	integer(kint) :: i

	call Copy(self%mesh, other%mesh)
	call Copy(self%relVort, other%relVort)
	call Copy(self%absVort, other%absVort)
	call Copy(self%velocity, other%velocity)
	call Copy(self%absStream, other%absStream)
	call Copy(self%relStream, other%relStream)
	call Copy(self%mpiParticles, other%mpiParticles)
	self%rotationRate = other%rotationRate
	self%radius = other%radius

	if ( allocated(other%tracers)) then
		if ( allocated(self%tracers)) then
			if ( size(self%tracers) /= size(other%tracers) ) then
				call LogMessage(log, ERROR_LOGGING_LEVEL, logkey, " copy BVEMesh ERROR : tracers not allocated properly.")
			else
				do i = 1, size(other%tracers)
					call Copy(self%tracers(i), other%tracers(i))
				enddo
			endif
		else
			call LogMessage(log, WARNING_LOGGING_LEVEL, logKey, " copy BVEMesh WARNING : tracers not copied.")
		endif
	endif
end subroutine

!> @brief Adds memory for passive tracers to a spherical BVE mesh.
!> @param[inout] self Target BVE mesh
!> @param[in] nTracers number of tracers to add
!> @param[in] tracerDims array of values (each 1, 2, or 3) corresponding to the dimensions of each tracer.
!> Scalar tracers will have dimension 1, vector tracers will have dimension 2 or 3, depending on their coordinate system.
subroutine AddTracers(self, nTracers, tracerDims)
	type(BVEMesh), intent(inout) :: self
	integer(kint), intent(in) :: nTracers
	integer(kint), dimension(:), intent(in) :: tracerDims
	integer(kint) :: i

	if ( size(tracerDims) /= nTracers ) then
		call LogMessage(log, ERROR_LOGGING_LEVEL, trim(logkey)//" AddTracers ERROR : ", &
			" must specify dimension of each tracer field.")
		return
	endif

	allocate(self%tracers(nTracers))
	do i = 1, nTracers
		call New(self%tracers(i), tracerDims(i), self%mesh%particles%N_Max)
	enddo
end subroutine

!> @brief Output basic information about a BVE mesh to a @ref Logger
!> @param[in] self BVE mesh
!> @param[inout] aLog @ref Logger
subroutine logStatsPrivate( self, aLog)
	type(BVEMEsh), intent(in) :: self
	type(Logger), intent(inout) :: aLog
	!
	integer(kint) :: i

	call LogMessage(aLog, TRACE_LOGGING_LEVEL, "BVEMesh ", "stats : ")
	call LogStats(self%mesh, aLog)
	call LogStats(self%velocity, aLog)
	call LogStats(self%relVort, aLog)
	call LogStats(self%absVort, aLog)
	call LogStats(self%relStream, aLog)
	call LogStats(self%absStream, aLog)
	if ( allocated( self%tracers ) ) then
		do i = 1, size(self%tracers)
			call LogStats(self%tracers(i), aLog)
		enddo
	endif
end subroutine

!> @brief Writes a BVE mesh to a legacy formatted .vtk file, including all @ref Field data for variables and tracers.
!> @param[in] self BVE mesh
!> @param[in] filename Name of output file
subroutine outputVTKPrivate(self, filename)
	type(BVEMesh), intent(in) :: self
	character(len=*), intent(in) :: filename
	integer(kint) :: i, writeStat

	open(unit=WRITE_UNIT_1, file=filename, status='REPLACE', action='WRITE', iostat=writeStat)
		if ( writeStat /= 0 ) then
			call LogMessage(log, ERROR_LOGGING_LEVEL, trim(logkey)//" OutputToVTK ERROR writing to file = ", trim(filename))
			return
		endif
		!
		!	write points and topology
		!
		call WriteVTKPoints(self%mesh%particles, WRITE_UNIT_1)
		call WriteFacesToVTKPolygons(self%mesh%faces, WRITE_UNIT_1)
		!
		!	write vtk point data
		!
		call WriteVTKPointDataSectionHeader(WRITE_UNIT_1, self%mesh%particles%N)
		call WriteVTKLagCoords( self%mesh%particles, WRITE_UNIT_1)
		call WriteFieldToVTKPointData( self%relVort, WRITE_UNIT_1)
		call WriteFieldToVTKPointData( self%absVort, WRITE_UNIT_1)
		call WriteFieldToVTKPointData( self%velocity, WRITE_UNIT_1)
		call WriteFieldToVTKPointData( self%relStream, WRITE_UNIT_1)
		call WriteFieldToVTKPointData( self%absStream, WRITE_UNIT_1)
		if ( allocated( self%tracers ) ) then
			do i = 1, size(self%tracers)
				call WriteFieldToVTKPointData( self%tracers(i), WRITE_UNIT_1)
			enddo
		endif
		!
		!	write vtk cell data
		!
		call WriteFaceAreaToVTKCellData( self%mesh%faces, self%mesh%particles, WRITE_UNIT_1)
	close(WRITE_UNIT_1)
end subroutine

!> @brief Defines an initial vorticity distribution on a BVE mesh.
!>
!> @param[inout] self BVE mesh
!> @param[in] relVortFn Vorticity distribution function.  Must have same interface as numberkindsmodule::scalarFnOf3DSpace
subroutine SetInitialVorticityOnMesh(self, relVortFn )
	type(BVEMesh), intent(inout) :: self
	procedure(scalarFnOf3DSpace) :: relVortFn
	integer(kint) :: i
	real(kreal) :: zeta

	self%relVort%N = self%mesh%particles%N
	self%absVort%N = self%mesh%particles%N
	do i = 1, self%mesh%particles%N
		zeta = relVortFn( self%mesh%particles%x0(i), self%mesh%particles%y0(i), self%mesh%particles%z0(i) )
		!call InsertScalarToField(self%relVort, zeta )
		!call InsertScalarToField(self%absVort, zeta + 2.0_kreal*self%rotationRate*self%mesh%particles%z0(i) / self%radius)
		self%relVort%scalar(i) = zeta
		self%absVort%scalar(i) = zeta + 2.0_kreal * self%rotationRate * self%mesh%particles%z0(i) / self%radius
	enddo
end subroutine

!> @brief Defines an initial velocity distribution on a BVE mesh.
!>
!> @param[inout] self BVE mesh
!> @param[in] velFn Velocity distribution function.  Must have same interface as numberkindsmodule::vectorFnOf3DSpace
subroutine setVelocityFromFunction(self, velFn )
	type(BVEMesh), intent(inout) :: self
	procedure(vectorFnOf3DSpace) :: velFn
	integer(kint) :: i

	do i = 1, self%mesh%particles%N
		call InsertVectorToField(self%velocity, &
			velFn( self%mesh%particles%x0(i), self%mesh%particles%y0(i), self%mesh%particles%z0(i)) )
	enddo
end subroutine

!> @brief Defines an initial tracer distribution on a BVE mesh.
!>
!> @param[inout] self BVE mesh
!> @param[in] tracerID index of tracer in `bveMesh%%tracers(:)`
!> @param[in] tracerFn Vorticity distribution function.  Must have same interface as numberkindsmodule::scalarFnOf3DSpace
subroutine SetScalarTracerOnMesh(self, tracerId, tracerFn )
	type(BVEMesh), intent(inout) :: self
	integer(kint), intent(in) :: tracerId
	procedure(scalarFnOf3DSpace) :: tracerFn
	integer(kint) :: i

	call SetFieldToZero(self%tracers(tracerID))
	self%tracers(tracerId)%N = self%mesh%particles%N
	do i = 1, self%mesh%particles%N
		self%tracers(tracerID)%scalar(i) = tracerFn( self%mesh%particles%x0(i), self%mesh%particles%y0(i), &
												     self%mesh%particles%z0(i) )
	enddo
end subroutine

!> @brief Defines an initial tracer distribution on a BVE mesh.
!>
!> @param[inout] self BVE mesh
!> @param[in] tracerID index of tracer in `bveMesh%%tracers(:)`
!> @param[in] tracerFn Vorticity distribution function.  Must have same interface as numberkindsmodule::vectorFnOf3DSpace
subroutine SetVectorTracerOnMesh( self, tracerID, tracerFn )
	type(BVEMesh), intent(inout) :: self
	integer(kint), intent(in) :: tracerId
	procedure(vectorFnOf3DSpace) :: tracerFn
	integer(kint) :: i
	real(kreal), dimension(3) :: vec

	call SetFieldToZero(self%tracers(tracerID))
	self%tracers(tracerID)%N = self%mesh%particles%N
	do i = 1, self%mesh%particles%N
		vec = tracerFn(self%mesh%particles%x0(i), self%mesh%particles%y0(i), self%mesh%particles%z0(i))
		self%tracers(tracerID)%xComp(i) = vec(1)
		self%tracers(tracerID)%yComp(i) = vec(2)
		self%tracers(tracerID)%zComp(i) = vec(3)
	enddo
end subroutine

!> @brief Computes the total kinetic energy (a conserved integral) on a BVE mesh.
!> @param[in] self BVE mesh
!> @return totalKE total kinetic energy
function TotalKE( self )
	real(kreal) :: TotaLKE
	type(BVEMesh), intent(in) :: self
	integer(kint) :: i
	real(kreal) :: magSq

	TotalKE = 0.0_kreal
	do i = 1, self%mesh%particles%N
		if ( self%mesh%particles%isActive(i) ) then
			magSq = self%velocity%xComp(i)**2 + self%velocity%yComp(i)**2 + self%velocity%zComp(i)**2
			TotalKE = TotalKE + magSq * self%mesh%particles%area(i)
		endif
	enddo
	TotalKE = 0.5_kreal * TotalKE
end function

!> @brief Computes the total enstrophy (a conserved integral) on a BVE mesh.
!> @param[in] self BVE mesh
!> @return totalKE total enstrophy
function TotalEnstrophy( self )
	real(kreal) :: TotalEnstrophy
	type(BVEMesh), intent(in) :: self
	integer(kint) :: i

	TotalEnstrophy = 0.0_kreal
	do i = 1, self%mesh%particles%N
		if ( self%mesh%particles%isActive(i) ) then
			TotalEnstrophy = TotalEnstrophy + self%relVort%scalar(i)**2 * self%mesh%particles%area(i)
		endif
	enddo
	TotalEnstrophy = 0.5_kreal * TotalEnstrophy
end function

!> @brief Defines the absolute and relative stream functions on a BVE mesh using the Green's function integral.
!> @param[inout] self BVE mesh
subroutine SetStreamFunctionsOnMesh( self )
	type(BVEMesh), intent(inout) :: self
	integer(kint) :: i, j, mpiErrCode
	real(kreal) :: greensKernel
	real(kreal), dimension(3) :: xi, xj

	self%relStream%n = self%mesh%particles%n
	self%absStream%n = self%mesh%particles%n
	do i = self%mpiParticles%indexStart(procRank), self%mpiParticles%indexEnd(procRank)
		self%relStream%scalar(i) = 0.0_kreal
		self%absStream%scalar(i) = 0.0_kreal
		xi = PhysCoord(self%mesh%particles, i)
		do j = 1, i - 1
			if ( self%mesh%particles%isActive(j) ) then
				xj = PhysCoord(self%mesh%particles, j)
				greensKernel = - dlog( self%radius * self%radius - sum(xi * xj)) / (4.0_kreal * PI )
				self%relStream%scalar(i) = self%relStream%scalar(i) + greensKernel * &
					self%relVort%scalar(j) * self%mesh%particles%area(j)
				self%absStream%scalar(i) = self%absStream%scalar(i) + greensKernel * &
					self%absVort%scalar(j) * self%mesh%particles%area(j)
			endif
		enddo
		do j = i + 1, self%mesh%particles%N
			if ( self%mesh%particles%isActive(j) ) then
				xj = PhysCoord(self%mesh%particles, j)
				greensKernel = - dlog( self%radius * self%radius - sum(xi * xj)) / (4.0_kreal * PI )
				self%relStream%scalar(i) = self%relStream%scalar(i) + greensKernel * &
					self%relVort%scalar(j) * self%mesh%particles%area(j)
				self%absStream%scalar(i) = self%absStream%scalar(i) + greensKernel * &
					self%absVort%scalar(j) * self%mesh%particles%area(j)
			endif
		enddo
	enddo

	do i = 0, numProcs - 1
		call MPI_BCAST(self%relStream%scalar(self%mpiParticles%indexStart(i):self%mpiParticles%indexEnd(i)), &
			self%mpiParticles%messageLength(i), MPI_DOUBLE_PRECISION, i, MPI_COMM_WORLD, mpiErrCode)
		call MPI_BCAST(self%absStream%scalar(self%mpiParticles%indexStart(i):self%mpiParticles%indexEnd(i)), &
			self%mpiParticles%messageLength(i), MPI_DOUBLE_PRECISION, i, MPI_COMM_WORLD, mpiErrCode)
	enddo
end subroutine

!> @brief Defines the velocity distribution on a BVE mesh using the Biot-Savart integral.
!> @param[inout] self BVE mesh
subroutine setVelocityFromVorticity( self )
	type(BVEMesh), intent(inout) :: self
	integer(kint) :: i, j, mpiErrCode
	real(kreal) :: strength
	real(kreal), dimension(3) :: xi, xj

	self%velocity%n = self%mesh%particles%N
	do i = self%mpiParticles%indexStart(procRank), self%mpiParticles%indexEnd(procRank)
		self%velocity%xComp(i) = 0.0_kreal
		self%velocity%yComp(i) = 0.0_kreal
		self%velocity%zComp(i) = 0.0_kreal
		xi = PhysCoord(self%mesh%particles, i)
		do j = 1, i - 1
			if ( self%mesh%particles%isActive(j) ) then
				xj = PhysCoord(self%mesh%particles, j)
				strength = -self%relVort%scalar(j)*self%mesh%particles%area(j) / &
					(4.0_kreal * PI * self%radius * (self%radius * self%radius - sum( xi * xj )))
				self%velocity%xComp(i) = self%velocity%xComp(i) + (xi(2)*xj(3) - xi(3)*xj(2)) * strength
				self%velocity%yComp(i) = self%velocity%yComp(i) + (xi(3)*xj(1) - xi(1)*xj(3)) * strength
				self%velocity%zComp(i) = self%velocity%zComp(i) + (xi(1)*xj(2) - xi(2)*xj(1)) * strength
			endif
		enddo
		do j = i + 1, self%mesh%particles%N
			if ( self%mesh%particles%isActive(j) ) then
				xj = PhysCoord(self%mesh%particles, j)
				strength = -self%relVort%scalar(j)*self%mesh%particles%area(j) / &
					(4.0_kreal * PI * self%radius * (self%radius * self%radius - sum( xi * xj )))
				self%velocity%xComp(i) = self%velocity%xComp(i) + (xi(2)*xj(3) - xi(3)*xj(2)) * strength
				self%velocity%yComp(i) = self%velocity%yComp(i) + (xi(3)*xj(1) - xi(1)*xj(3)) * strength
				self%velocity%zComp(i) = self%velocity%zComp(i) + (xi(1)*xj(2) - xi(2)*xj(1)) * strength
			endif
		enddo
	enddo

	do i = 0, numProcs - 1
		call MPI_BCAST(self%velocity%xComp(self%mpiParticles%indexStart(i):self%mpiParticles%indexEnd(i)), &
					   self%mpiParticles%messageLength(i), MPI_DOUBLE_PRECISION, i, MPI_COMM_WORLD, mpiErrCode)
		call MPI_BCAST(self%velocity%yComp(self%mpiParticles%indexStart(i):self%mpiParticles%indexEnd(i)), &
					   self%mpiParticles%messageLength(i), MPI_DOUBLE_PRECISION, i, MPI_COMM_WORLD, mpiErrCode)
		call MPI_BCAST(self%velocity%zComp(self%mpiParticles%indexStart(i):self%mpiParticles%indexEnd(i)), &
					   self%mpiParticles%messageLength(i), MPI_DOUBLE_PRECISION, i, MPI_COMM_WORLD, mpiErrCode)
	enddo
end subroutine

!> @brief Computes the maximum circulation magnitude about each face and returns the maximum value for the whole mesh.
!> @param[in] self BVE mesh
!> @return MaximumCirculation magnitude
function MaxCirculationMagnitudePerFace( self )
	real(kreal) :: MaxCirculationMagnitudePerFace
	type(BVEMesh), intent(in) :: self
	!
	integer(kint) :: i
	real(kreal) :: testCirc
	integer(kint) :: particleIndex

	MaxCirculationMagnitudePerFace = 0.0_kreal
	do i = 1, self%mesh%faces%N
		if ( .NOT. self%mesh%faces%hasChildren(i) ) then
			particleIndex = self%mesh%faces%centerParticle(i)
			testCirc = abs( self%relVort%scalar(particleIndex) ) * self%mesh%particles%area(particleIndex)
			if ( testCirc > MaxCirculationMagnitudePerFace ) MaxCirculationMagnitudePerFace = testCirc
		endif
	enddo
end function

!
!----------------
! private methods
!----------------
!

!> @brief Initializes a logger for the BVE module
!>
!> Output is controlled both by message priority and by MPI Rank
!> @param aLog Target Logger object
!> @param rank Rank of this processor
subroutine InitLogger(aLog,rank)
	type(Logger), intent(out) :: aLog
	integer(kint), intent(in) :: rank
	write(logKey,'(A,A,I0.2,A)') trim(logKey),'_',rank,' : '
	if ( rank == 0 ) then
		call New(aLog,logLevel)
	else
		call New(aLog,ERROR_LOGGING_LEVEL)
	endif
	logInit = .TRUE.
end subroutine

!> @}
end module


text:
!==============================================================================!
  module Gen_Mod
!------------------------------------------------------------------------------!
!   Global variable definitions for the mesh generator.                        !
!------------------------------------------------------------------------------!
  implicit none
!==============================================================================!

  ! Twin nodes
  integer, allocatable :: twin_n(:,:)

  ! For each face, which neighbours are to each other cells which meet there
  ! The name derives from "face cell to cell"
  integer, allocatable :: face_c_to_c(:,:)

  ! Periodic and copy boundaries
  integer              :: n_periodic_cond,      n_copy_cond
  integer, allocatable ::   periodic_cond(:,:),   copy_cond(:,:)

end module

max_stars_repo_name: MassimoBorrelliPhysicist/T-Flows

max_issues_repo_name: EdinSmartLab/T-Flows

max_forks_repo_name: EdinSmartLab/T-Flows


text:
program dump_grids
  use job_manage, only: job_fork, time_message
  use mp, only: init_mp, proc0, nproc, broadcast
  use file_utils, only: init_file_utils, run_name
  use kt_grids, only: init_kt_grids
  use theta_grid, only: init_theta_grid
  use le_grids, only: init_le_grids

  implicit none
  logical :: list
  logical :: accel_x, accel_y
  logical, parameter :: quiet=.false.
  character (500), target :: cbuff
  character (len=500) :: filename
  real :: time_measure(2)

  ! <doc> Initialize message passing </doc>
  call init_mp

  ! <doc> Report # of processors being used </doc>
  if (proc0) then
     if(.not.quiet)then
        if (nproc == 1) then
           write(*,'("Running on ",I0," processor")') nproc
        else
           write(*,'("Running on ",I0," processors")') nproc
        end if
     endif

     ! <doc> Call init_file_utils, ie. initialize the inputs and outputs, checking
     !  whether we are doing a run or a list of runs. </doc>
     call init_file_utils (list, name="gs")
  end if

  call broadcast (list)

  ! <doc> If given a list of jobs, fork </doc>
  if (list) call job_fork

  if (proc0) cbuff = trim(run_name)
  call broadcast (cbuff)
  if (.not. proc0) run_name => cbuff

  !Initialise
  time_measure=0.
  !/Theta
  if(proc0.and.(.not.quiet)) then
     write(*,'("Init theta_grids")',advance='no')
     call time_message(.false.,time_measure,'init-theta')
  endif
  call init_theta_grid
  if(proc0.and.(.not.quiet)) then
     call time_message(.false.,time_measure,'init-theta')
     write(*,'("  --> Done in : ",F12.6," seconds")') time_measure(1)
  endif
  time_measure=0.
  !/KT
  if(proc0.and.(.not.quiet)) then
     write(*,'("Init kt_grids   ")',advance='no')
     call time_message(.false.,time_measure,'init-kt')
  endif
  call init_kt_grids
  if(proc0.and.(.not.quiet)) then
     call time_message(.false.,time_measure,'init-kt')
     write(*,'("  --> Done in : ",F12.6," seconds")') time_measure(1)
  endif
  time_measure=0.
  !/LE
  if(proc0.and.(.not.quiet)) then
     write(*,'("Init le_grids   ")',advance='no')
     call time_message(.false.,time_measure,'init-le')
  endif
  call init_le_grids(accel_x,accel_y)
  if(proc0.and.(.not.quiet)) then
     call time_message(.false.,time_measure,'init-le')
     write(*,'("  --> Done in : ",F12.6," seconds")') time_measure(1)
  endif
  time_measure=0.

  !Now write to file
  if(proc0.and.(.not.quiet)) then
     write(*,'("Write to file   ")',advance='no')
     call time_message(.false.,time_measure,'write-file')
  endif
  filename=trim(adjustl(run_name))
  if(proc0) call write_grids_to_file(filename)
  if(proc0.and.(.not.quiet)) then
     call time_message(.false.,time_measure,'write-file')
     write(*,'("  --> Done in : ",F12.6," seconds")') time_measure(1)
  endif
  time_measure=0.

contains

#ifdef NETCDF
!Prefer to write to netcdf file if possible
  subroutine write_grids_to_file(fname)
    use netcdf, only: NF90_CLOBBER
    use netcdf, only: nf90_create, nf90_close
    use netcdf, only: nf90_def_dim, nf90_def_var, nf90_enddef
    use netcdf, only: nf90_put_var

    use netcdf_utils, only: get_netcdf_code_precision
    use netcdf_utils, only: netcdf_real

    use kt_grids, only: aky, akx, theta0
    use theta_grid, only: theta
    use le_grids, only: al, energy

    implicit none
    character(len=*), intent(in) :: fname
    integer :: ncid, ierr
    integer :: ky_dimid, kx_dimid, theta_dimid, al_dimid, energy_dimid
    integer :: ky_varid, kx_varid, theta_varid, al_varid, energy_varid, theta0_varid

    !First get precisions
    if(netcdf_real==0) netcdf_real=get_netcdf_code_precision()

    !First create a file
    ierr=nf90_create(trim(adjustl(fname))//".grids.nc",NF90_CLOBBER,ncid)

    !Define dimensions
    ierr=nf90_def_dim(ncid,'aky',size(aky),ky_dimid)
    ierr=nf90_def_dim(ncid,'akx',size(akx),kx_dimid)
    ierr=nf90_def_dim(ncid,'theta',size(theta),theta_dimid)
    ierr=nf90_def_dim(ncid,'lambda',size(al),al_dimid)
    ierr=nf90_def_dim(ncid,'energy',size(energy),energy_dimid)

    !Define variables
    ierr=nf90_def_var(ncid,'ky',netcdf_real,(/ky_dimid/),ky_varid)
    ierr=nf90_def_var(ncid,'kx',netcdf_real,(/kx_dimid/),kx_varid)
    ierr=nf90_def_var(ncid,'theta0',netcdf_real,(/kx_dimid,ky_dimid/),theta0_varid)
    ierr=nf90_def_var(ncid,'theta',netcdf_real,(/theta_dimid/),theta_varid)
    ierr=nf90_def_var(ncid,'lambda',netcdf_real,(/al_dimid/),al_varid)
    ierr=nf90_def_var(ncid,'energy',netcdf_real,(/energy_dimid/),energy_varid)

    !End definitions
    ierr=nf90_enddef(ncid)

    !Now fill in variable data
    ierr=nf90_put_var(ncid,ky_varid,aky)
    ierr=nf90_put_var(ncid,kx_varid,akx)
    ierr=nf90_put_var(ncid,theta0_varid,theta0)
    ierr=nf90_put_var(ncid,theta_varid,theta)
    ierr=nf90_put_var(ncid,al_varid,al)
    ierr=nf90_put_var(ncid,energy_varid,energy)

    !Now close file
    ierr=nf90_close(ncid)
  end subroutine write_grids_to_file
#else
!Fall back to binary output when netcdf unavailable
  subroutine write_grids_to_file(fname)
    use kt_grids, only: aky, akx, theta0
    use theta_grid, only: theta
    use le_grids, only: al, energy
    use file_utils, only: get_unused_unit

    implicit none
    character(len=*), intent(in) :: fname
    integer :: unit

    call get_unused_unit(unit)

    open(unit=unit,file=trim(adjustl(fname))//".ky",form="unformatted")
    write(unit) aky
    close(unit)

    open(unit=unit,file=trim(adjustl(fname))//".kx",form="unformatted")
    write(unit) akx
    close(unit)

    open(unit=unit,file=trim(adjustl(fname))//".theta0",form="unformatted")
    write(unit) theta0
    close(unit)

    open(unit=unit,file=trim(adjustl(fname))//".theta",form="unformatted")
    write(unit) theta
    close(unit)

    open(unit=unit,file=trim(adjustl(fname))//".lambda",form="unformatted")
    write(unit) al
    close(unit)

    open(unit=unit,file=trim(adjustl(fname))//".energy",form="unformatted")
    write(unit) energy
    close(unit)
  end subroutine write_grids_to_file
#endif
end program dump_grids

ext: fpp


text:
c{ dg-do run { target fd_truncate } }
c{ dg-options "-std=legacy" }
c
c This program repeats many of the same tests as test_nml_1 but for integer
c instead of real. It also tests repeat nulls, comma delimited character read,
c a triplet qualifier, a range with an assumed start, a quote delimited string,
c a qualifier with an assumed end and a fully explicit range.  It also tests
c that integers and characters are successfully read back by namelist.
c Provided by Paul Thomas - pault@gcc.gnu.org

      program namelist_12

      integer x(10)
      integer(kind=8) xx
      integer ier
      character*10 ch , check
      namelist /mynml/ x, xx, ch

c set debug = 0 or 1 in the namelist! (line 33)

      do i = 1 , 10
        x(i) = -1
      end do
      x(6) = 6
      x(10) = 10
      xx = 0
      ch ="zzzzzzzzzz"
      check="abcdefghij"

      open (10,status="scratch", delim="apostrophe")
      write (10, '(a)') "!mynml"
      write (10, '(a)') " "
      write (10, '(a)') "&mynml  x(7) =+99 x=1, 2 ,"
      write (10, '(a)') " 2*3, ,, 2* !comment"
      write (10, '(a)') " 9 ch='qqqdefghqq' , x(8:7:-1) = 8 , 7"
      write (10, '(a)') " ch(:3) =""abc"","
      write (10, '(a)') " ch(9:)='ij' x(4:5)=4 ,5 xx = 42/"
      rewind (10)

      read (10, nml=mynml, IOSTAT=ier)
      if (ier.ne.0) STOP 1
      rewind (10)

      write (10, nml=mynml, iostat=ier)
      if (ier.ne.0) STOP 2
      rewind (10)

      read (10, NML=mynml, IOSTAT=ier)
      if (ier.ne.0) STOP 3
      close (10)

      do i = 1 , 10
        if ( abs( x(i) - i ) .ne. 0 ) STOP 1
        if ( ch(i:i).ne.check(I:I) ) STOP 4
      end do
      if (xx.ne.42) STOP 2
      end program