MODULE timeloop_gcm_mod

   USE icosa

   IMPLICIT NONE !(TYPE, EXTERNAL)
   PRIVATE

   INTEGER, PARAMETER :: sync_it = 10
   TYPE(t_message), SAVE :: req_ps0, req_rhodz0, req_mass0, req_theta_rhodz0, req_u0, req_q0, req_W0, req_geopot0
   LOGICAL, SAVE :: positive_theta
   INTEGER :: itau_prof, id_timeloop, id_dyn, id_phys, id_dissip, id_adv, id_diags
   PUBLIC :: init_timeloop, timeloop

CONTAINS

   SUBROUTINE init_timeloop
      USE abort_mod, ONLY : abort_acc
      USE profiling_mod, ONLY: register_id
      USE dissip_gcm_mod, ONLY: init_dissip
      USE observable_mod, ONLY: init_observable
      USE etat0_mod, ONLY: etat0
      USE guided_mod, ONLY: init_guided
      USE advect_tracer_mod, ONLY: init_advect_tracer
      USE check_conserve_mod, ONLY: init_check_conserve
      USE output_field_mod, ONLY: enable_io
      USE theta2theta_rhodz_mod, ONLY: init_theta2theta_rhodz
      USE sponge_mod, ONLY: init_sponge
      USE time_scheme_base_mod, ONLY: scheme_family, scheme, nb_stage, nb_storage, matsuno_period
      USE time_scheme_base_mod, ONLY: explicit, hevi, euler, rk4, mlf, rk25, ark23, ark25, ark33, ark11
      USE time_scheme_base_mod, ONLY: f_phis, f_ps, f_mass, f_rhodz, f_theta_rhodz, f_u
      USE time_scheme_base_mod, ONLY: f_geopot, f_W, f_q, f_hflux, f_hfluxt, f_wflux, f_wfluxt
      USE time_scheme_base_mod, ONLY: f_dps, f_dmass, f_dtheta_rhodz, f_du
      USE time_scheme_base_mod, ONLY: f_psm1, f_massm1, f_theta_rhodzm1, f_um1
      USE time_scheme_base_mod, ONLY: f_dps_slow, f_dmass_slow, f_dtheta_rhodz_slow
      USE time_scheme_base_mod, ONLY: f_dtheta_rhodz_slow, f_du_slow, f_dW_slow, f_dPhi_slow
      USE time_scheme_base_mod, ONLY: f_du_fast, f_dW_fast, f_dPhi_fast
      USE time_scheme_base_mod, ONLY: f_psm2, f_massm2, f_theta_rhodzm2, f_um2
      USE hevi_scheme_mod, ONLY: set_coefs_ark33, set_coefs_ark23, set_coefs_ark11, set_coefs_ark25
      USE caldyn_vars_mod, ONLY: caldyn_conserv, conserv_gassmann

      CHARACTER(len=255) :: def

      CALL register_id('timeloop', id_timeloop)
      CALL register_id('dyn', id_dyn)
      CALL register_id('dissip', id_dissip)
      CALL register_id('phys', id_phys)
      CALL register_id('adv', id_adv)
      CALL register_id('diags', id_diags)

!    IF (xios_output) itau_out=1
      IF (.NOT. enable_io) itau_out = HUGE(itau_out)

      itau_prof = 100
      CALL getin('itau_profiling', itau_prof)

      positive_theta = .FALSE.
      CALL getin('positive_theta', positive_theta)
      IF (positive_theta .AND. nqtot < 1) THEN
         PRINT *, 'nqtot must be >0 if positive_theta is .TRUE.'
         STOP
      END IF

      nb_storage = 0 ! default value, set later to nb_stage unless set explicitly below
      def = 'ARK2.3'
      CALL getin('time_scheme', def)
      SELECT CASE (TRIM(def))
      CASE ('euler')
         scheme_family = explicit
         scheme = euler
         nb_stage = 1
      CASE ('RK4.4')
         scheme_family = explicit
         scheme = rk4
         nb_stage = 4
      CASE ('RK2.5')
         scheme_family = explicit
         scheme = rk25
         nb_stage = 5
      CASE ('leapfrog_matsuno')
         scheme_family = explicit
         scheme = mlf
         matsuno_period = 5
         CALL getin('matsuno_period', matsuno_period)
         nb_stage = matsuno_period + 1
      CASE ('ARK1.1')
         scheme_family = hevi
         scheme = ark11
         nb_stage = 1
         CALL set_coefs_ark11(dt)
      CASE ('ARK2.3')
         scheme_family = hevi
         scheme = ark23
         nb_stage = 3
         CALL set_coefs_ark23(dt)
      CASE ('ARK2.5')
         scheme_family = hevi
         scheme = ark25
         nb_stage = 5
         nb_storage = 2 ! this is a low-storage scheme
         CALL set_coefs_ark25(dt)
      CASE ('ARK3.3')
         scheme_family = hevi
         scheme = ark33
         nb_stage = 3
         CALL set_coefs_ark33(dt)
      CASE ('none')
         nb_stage = 0
      CASE default
         PRINT *, 'Bad selector for variable time_scheme : <', TRIM(def), &
            ' > options are <none>, <euler>, <leapfrog_matsuno>, <RK2.5>, <RK4.4>, <ARK1.1>, <ARK2.3>, <ARK2.5>, <ARK3.3>'
         STOP
      END SELECT

      IF (scheme_family /= hevi) THEN
         CALL abort_acc("scheme_family /= hevi")
         IF(grid_type/=grid_ico) THEN
            PRINT *, "Selected time scheme is time_scheme=", TRIM(def)
            PRINT *, "Unstructured meshes require time_scheme=ARK2.3 or ARK2.5"
            STOP
         END IF
      END IF

      ! Time-independant orography
      CALL allocate_field(f_phis, field_t, type_real, name='phis')
      ! Model state at current time step
      CALL allocate_field(f_ps, field_t, type_real, name='ps')
      CALL allocate_field(f_mass, field_t, type_real, llm, name='mass')
      CALL allocate_field(f_rhodz, field_t, type_real, llm, name='rhodz')
      CALL allocate_field(f_theta_rhodz, field_t, type_real, llm, nqdyn, name='theta_rhodz')
      CALL allocate_field(f_u, field_u, type_real, llm, name='u')
      CALL allocate_field(f_geopot, field_t, type_real, llm + 1, name='geopot')
      CALL allocate_field(f_W, field_t, type_real, llm + 1, name='W') ! used only if .not. hydrostatic
      CALL allocate_field(f_q, field_t, type_real, llm, nqtot, 'q')
      ! Mass fluxes
      CALL allocate_field(f_hflux, field_u, type_real, llm, ondevice=.TRUE.)    ! instantaneous mass fluxes
      CALL allocate_field(f_hfluxt, field_u, type_real, llm, ondevice=.TRUE.)   ! mass "fluxes" accumulated in time
      CALL allocate_field(f_wflux, field_t, type_real, llm + 1)  ! vertical mass fluxes
      CALL allocate_field(f_wfluxt, field_t, type_real, llm + 1, name='wfluxt', ondevice=.TRUE.)

      SELECT CASE (scheme_family)
      CASE (explicit)
         IF (caldyn_conserv == conserv_gassmann) THEN
            PRINT *, 'caldyn_conserv = energy_gassmann is not implemented for RK time schemes, use time_scheme=ARK2.3 .'
            STOP
         END IF
         ! Trends
         CALL allocate_field(f_dps, field_t, type_real, name='dps')
         CALL allocate_field(f_dmass, field_t, type_real, llm, name='dmass')
         CALL allocate_field(f_dtheta_rhodz, field_t, type_real, llm, nqdyn, name='dtheta_rhodz')
         CALL allocate_field(f_du, field_u, type_real, llm, name='du')
         ! Model state at previous time step (RK/MLF)
         CALL allocate_field(f_psm1, field_t, type_real, name='psm1')
         CALL allocate_field(f_massm1, field_t, type_real, llm, name='massm1')
         CALL allocate_field(f_theta_rhodzm1, field_t, type_real, llm, nqdyn, name='theta_rhodzm1')
         CALL allocate_field(f_um1, field_u, type_real, llm, name='um1')
      CASE (hevi)
         if (nb_storage == 0) nb_storage = nb_stage
         ! Trends
         CALL allocate_fields(nb_storage, f_dps_slow, field_t, type_real, name='dps_slow', ondevice=.TRUE.)
         CALL allocate_fields(nb_storage, f_dmass_slow, field_t, type_real, llm, name='dmass_slow', ondevice=.TRUE.)
         CALL allocate_fields(nb_storage, f_dtheta_rhodz_slow, field_t, type_real, llm, nqdyn, &
                              name='dtheta_rhodz_fast', ondevice=.TRUE.)
         CALL allocate_fields(nb_storage, f_du_slow, field_u, type_real, llm, name='du_slow', ondevice=.TRUE.)
         CALL allocate_fields(nb_storage, f_du_fast, field_u, type_real, llm, name='du_fast', ondevice=.TRUE.)
         CALL allocate_fields(nb_storage, f_dW_slow, field_t, type_real, llm + 1, name='dW_slow')
         CALL allocate_fields(nb_storage, f_dW_fast, field_t, type_real, llm + 1, name='dW_fast')
         CALL allocate_fields(nb_storage, f_dPhi_slow, field_t, type_real, llm + 1, name='dPhi_slow')
         CALL allocate_fields(nb_storage, f_dPhi_fast, field_t, type_real, llm + 1, name='dPhi_fast')
         f_dps => f_dps_slow(:, 1)
         f_du => f_du_slow(:, 1)
         f_dtheta_rhodz => f_dtheta_rhodz_slow(:, 1)
      END SELECT

      SELECT CASE (scheme)
      CASE (mlf)
         ! Model state 2 time steps ago (MLF)
         CALL allocate_field(f_psm2, field_t, type_real)
         CALL allocate_field(f_massm2, field_t, type_real, llm)
         CALL allocate_field(f_theta_rhodzm2, field_t, type_real, llm, nqdyn)
         CALL allocate_field(f_um2, field_u, type_real, llm)
      END SELECT

      CALL init_theta2theta_rhodz
      CALL init_dissip
      CALL init_sponge
      CALL init_observable
      CALL init_advect_tracer
      CALL init_check_conserve

      CALL etat0(f_ps, f_mass, f_phis, f_theta_rhodz, f_u, f_geopot, f_W, f_q)
      CALL init_guided(f_u, f_theta_rhodz, f_ps)

      CALL transfert_request(f_phis, req_i0)
      CALL transfert_request(f_phis, req_i1)

      CALL init_message(f_ps, req_i0, req_ps0)
      CALL init_message(f_mass, req_i0, req_mass0)
      CALL init_message(f_rhodz, req_i0, req_rhodz0)
      CALL init_message(f_theta_rhodz, req_i0, req_theta_rhodz0)
      CALL init_message(f_u, req_e0_vect, req_u0)
      CALL init_message(f_q, req_i0, req_q0)
      CALL init_message(f_geopot, req_i0, req_geopot0)
      CALL init_message(f_W, req_i0, req_W0)

   END SUBROUTINE init_timeloop

   SUBROUTINE timeloop
      USE abort_mod, ONLY : abort_acc
      USE trace, ONLY: trace_on
      USE omp_para, ONLY: is_master, is_omp_master, ll_begin, ll_end
      USE omp_para, ONLY: switch_omp_distrib_level
      USE profiling_mod, ONLY: enter_profile, exit_profile
      USE grid_param, ONLY: grid_type, grid_ico
      USE disvert_mod, ONLY: caldyn_eta, eta_mass
      USE dissip_gcm_mod, ONLY: dissip
      USE sponge_mod, ONLY: iflag_sponge, sponge
      USE observable_mod, ONLY: f_theta, f_buf_i, write_output_fields_basic
      USE compute_diagnostics_mod, ONLY: compute_rhodz
      USE guided_mod, ONLY: guided
      USE caldyn_mod, ONLY: caldyn_BC
      USE advect_tracer_mod, ONLY: advect_tracer
      USE diagflux_mod, ONLY: diagflux_on, diagflux_energy
      USE diagflux_mod, ONLY: f_masst, f_qmasst, f_massfluxt, f_qfluxt
      USE physics_mod, ONLY: physics
      USE mpipara, ONLY: mpi_rank
      USE check_conserve_mod, ONLY: check_conserve, check_conserve_detailed
      USE check_conserve_mod, ONLY: AAM_dyn, AAM_dissip, AAM_phys
      USE xios_mod, ONLY: xios_update_calendar
      USE output_field_mod, ONLY: xios_output
      USE write_etat0_mod, ONLY: write_etat0
      USE restart_mod, ONLY: is_write_start
      USE euler_scheme_mod, ONLY: euler_scheme
      USE time_scheme_base_mod, ONLY: scheme_family
      USE time_scheme_base_mod, ONLY: explicit, hevi
      USE time_scheme_base_mod, ONLY: f_phis, f_ps, f_mass, f_rhodz, f_theta_rhodz, f_u
      USE time_scheme_base_mod, ONLY: f_geopot, f_W, f_q, f_hflux, f_hfluxt, f_wflux, f_wfluxt
      USE time_scheme_base_mod, ONLY: f_dps, f_dtheta_rhodz, f_du
      USE explicit_scheme_mod, ONLY: explicit_scheme
      USE hevi_scheme_mod, ONLY: HEVI_scheme
      USE caldyn_vars_mod, ONLY: hydrostatic, caldyn_conserv, conserv_gassmann

      REAL(rstd), POINTER :: rhodz(:, :), mass(:, :), ps(:)

      REAL(rstd) :: adv_over_out ! ratio itau_adv/itau_out
      INTEGER :: ind, it, l
      LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating mass fluxes in time
      LOGICAL, PARAMETER :: check_rhodz = .FALSE.
      INTEGER(kind=8) :: start_clock, stop_clock, rate_clock
      INTEGER :: i, itau_sync   ! best iteration for synchronisation and ensure 1+1=2
      LOGICAL, SAVE :: first_physic = .TRUE.
      !$OMP THREADPRIVATE(first_physic)

      itau_sync = 1
      DO i = 2, 3*sync_it
         IF (MOD(86400, INT(i*dt)) == 0 .AND. &
             ABS((sync_it - itau_sync)*1./sync_it)/sync_it < (sync_it - itau_sync)*1./sync_it) itau_sync = i
      END DO
      CALL getin("itau_sync", itau_sync)
      IF (is_master) PRINT *, "Synchronize frontier every itau_sync =", itau_sync

      CALL switch_omp_distrib_level
      CALL caldyn_BC(f_phis, f_geopot, f_wflux) ! set constant values in first/last interfaces

      !$OMP BARRIER
      DO ind = 1, ndomain
         IF (.NOT. assigned_domain(ind)) CYCLE
         CALL swap_dimensions(ind)
         CALL swap_geometry(ind)
         rhodz = f_rhodz(ind); mass = f_mass(ind); ps = f_ps(ind)
         IF (caldyn_eta == eta_mass) THEN
            CALL compute_rhodz(.TRUE., ps, rhodz) ! save rhodz for transport scheme before dynamics update ps
         ELSE
            IF (grid_type == grid_ico) THEN
               DO l = ll_begin, ll_end
                  rhodz(:, l) = mass(:, l)
               END DO
            ELSE
               BLOCK
                  INTEGER :: ij
                  !$OMP DO SCHEDULE(STATIC)
                  DO ij = 1, primal_num
                     rhodz(:, ij) = mass(:, ij)
                  END DO
                  !$OMP END DO
               END BLOCK
            END IF
         END IF
      END DO
      !$OMP BARRIER

      fluxt_zero = .TRUE.

      IF (positive_theta) CALL copy_theta_to_q(f_theta_rhodz, f_rhodz, f_q)
      IF (diagflux_on) THEN
         adv_over_out = itau_adv*(1./itau_out)
      ELSE
         adv_over_out = 0.
      END IF

      CALL check_conserve(f_ps, f_dps, f_u, f_theta_rhodz, f_phis, f_q, itau0)

      Call trace_on

      IF (xios_output) THEN ! we must call update_calendar before any XIOS output
         IF (is_omp_master) CALL xios_update_calendar(1)
      END IF

      IF (is_write_start) THEN
         CALL write_etat0(itau0, f_ps, f_phis, f_theta_rhodz, f_u, f_q, f_geopot, f_W)
      ENDIF

      CALL write_output_fields_basic(.TRUE., f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q)

      !$OMP MASTER
      CALL SYSTEM_CLOCK(start_clock, rate_clock)
      !$OMP END MASTER
      call update_device_field(f_ps)
      call update_device_field(f_mass)
      CALL update_device_field(f_theta_rhodz)
      CALL update_device_field(f_u)
      CALL update_device_field(f_q)
      CALL update_device_field(f_geopot)
      CALL update_device_field(f_wflux)
      CALL update_device_field(f_rhodz)

      DO it = itau0 + 1, itau0 + itaumax
         CALL print_iteration(it, itau0, itaumax, start_clock, rate_clock)

         CALL enter_profile(id_timeloop)

         IF (xios_output) THEN

            IF (it > itau0 + 1 .AND. is_omp_master) CALL xios_update_calendar(it - itau0)
         ELSE
            CALL update_time_counter(dt*it)
         END IF

         IF (it == itau0 + 1 .OR. MOD(it - 1, itau_sync) == 0) THEN
            CALL send_message(f_ps, req_ps0)
            CALL wait_message(req_ps0)
            CALL send_message(f_mass, req_mass0)
            CALL wait_message(req_mass0)
            CALL send_message(f_rhodz, req_rhodz0)
            CALL wait_message(req_rhodz0)
            CALL send_message(f_theta_rhodz, req_theta_rhodz0)
            CALL wait_message(req_theta_rhodz0)
            CALL send_message(f_u, req_u0)
            CALL wait_message(req_u0)
            CALL send_message(f_q, req_q0)
            CALL wait_message(req_q0)
            IF (.NOT. hydrostatic) THEN
               CALL send_message(f_geopot, req_geopot0)
               CALL wait_message(req_geopot0)
               CALL send_message(f_W, req_W0)
               CALL wait_message(req_W0)
            END IF
         END IF

         CALL guided(it*dt, f_ps, f_theta_rhodz, f_u, f_q)

         CALL enter_profile(id_dyn)
         SELECT CASE (scheme_family)
         CASE (explicit)
            CALL abort_acc("explicit_scheme")
            CALL explicit_scheme(it, fluxt_zero)
         CASE (hevi)
            CALL HEVI_scheme(it, fluxt_zero)
         END SELECT
         CALL exit_profile(id_dyn)

         CALL enter_profile(id_dissip)
         IF (MOD(it, itau_dissip) == 0) THEN

            IF (caldyn_eta == eta_mass) THEN
               !ym flush ps
               !$OMP BARRIER
               DO ind = 1, ndomain
                  IF (.NOT. assigned_domain(ind)) CYCLE
                  CALL swap_dimensions(ind)
                  CALL swap_geometry(ind)
                  mass = f_mass(ind); ps = f_ps(ind);
                  CALL compute_rhodz(.TRUE., ps, mass)
               END DO
            END IF

            CALL enter_profile(id_diags)
            CALL check_conserve_detailed(it, AAM_dyn, &
                                         f_ps, f_dps, f_u, f_theta_rhodz, f_phis)
            CALL exit_profile(id_diags)

            CALL dissip(f_ps, f_mass, f_phis, f_geopot, f_theta_rhodz, f_u, f_dtheta_rhodz, f_du)

            CALL euler_scheme(.FALSE.)  ! update only u, theta
            IF (iflag_sponge > 0) THEN
               CALL sponge(f_u, f_du, f_theta_rhodz, f_dtheta_rhodz)
               CALL euler_scheme(.FALSE.)  ! update only u, theta
            END IF

            CALL enter_profile(id_diags)
            CALL check_conserve_detailed(it, AAM_dissip, &
                                         f_ps, f_dps, f_u, f_theta_rhodz, f_phis)
            CALL exit_profile(id_diags)
         END IF
         CALL exit_profile(id_dissip)

         CALL enter_profile(id_adv)
         IF (MOD(it, itau_adv) == 0) THEN
            CALL advect_tracer(f_hfluxt, f_wfluxt, f_u, f_q, f_rhodz, & ! update q and rhodz after RK step
                               adv_over_out, f_masst, f_qmasst, f_massfluxt, f_qfluxt,  f_theta_rhodz, f_ps)  ! accumulate mass and fluxes if diagflux_on
            fluxt_zero = .TRUE. ! restart accumulation of hfluxt and qfluxt at next call to explicit_scheme / HEVI_scheme

            ! At this point advect_tracer has obtained the halos of u and rhodz,
            ! needed for correct computation of kinetic energy
            IF (grid_type == grid_ico) THEN ! FIXME unst

               IF (diagflux_on) CALL abort_acc("diagflux_on")
               IF (diagflux_on) CALL diagflux_energy(adv_over_out, f_phis, f_rhodz, f_theta_rhodz, &
                                                     f_u, f_geopot, f_theta, f_buf_i, f_hfluxt)
            END IF

            IF (check_rhodz) THEN ! check that rhodz is consistent with ps
               DO ind = 1, ndomain
                  IF (.NOT. assigned_domain(ind)) CYCLE
                  CALL swap_dimensions(ind)
                  CALL swap_geometry(ind)
                  rhodz = f_rhodz(ind); ps = f_ps(ind);
                  CALL compute_rhodz(.FALSE., ps, rhodz)
               END DO
            END IF
            IF (positive_theta) CALL copy_q_to_theta(f_theta_rhodz, f_rhodz, f_q)
         END IF
         CALL exit_profile(id_adv)

         CALL enter_profile(id_diags)
         IF (MOD(it, itau_physics) == 0) THEN
            CALL check_conserve_detailed(it, AAM_dyn, &
                                         f_ps, f_dps, f_u, f_theta_rhodz, f_phis)
            CALL enter_profile(id_phys)
            CALL physics(it, f_phis, f_geopot, f_ps, f_theta_rhodz, f_u, f_wflux, f_q)
            CALL exit_profile(id_phys)
            CALL check_conserve_detailed(it, AAM_phys, &
                                         f_ps, f_dps, f_u, f_theta_rhodz, f_phis)
            !$OMP MASTER
            IF (first_physic) CALL SYSTEM_CLOCK(start_clock)
            !$OMP END MASTER
            first_physic = .FALSE.
         END IF

         IF (MOD(it, itau_check_conserv) == 0) THEN
            CALL check_conserve_detailed(it, AAM_dyn, &
                                         f_ps, f_dps, f_u, f_theta_rhodz, f_phis)
            CALL check_conserve(f_ps, f_dps, f_u, f_theta_rhodz, f_phis, f_q, it)
         END IF

         IF (mod(it, itau_out) == 0) THEN
            CALL transfert_request(f_u, req_e1_vect)
            CALL update_host_field(f_ps)
            CALL update_host_field(f_mass)
            CALL update_host_field(f_theta_rhodz)
            CALL update_host_field(f_geopot)
            CALL update_host_field(f_u)
            CALL update_host_field(f_q)
            CALL write_output_fields_basic(.FALSE., f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q)
         END IF
         CALL exit_profile(id_diags)

         CALL exit_profile(id_timeloop)
      END DO

      CALL update_host_field(f_ps)
      CALL update_host_field(f_theta_rhodz)
      CALL update_host_field(f_u)
      CALL update_host_field(f_q)
      CALL update_host_field(f_geopot)

      CALL write_etat0(itau0 + itaumax, f_ps, f_phis, f_theta_rhodz, f_u, f_q, f_geopot, f_W)

      CALL update_host_field(f_dps)
      CALL check_conserve_detailed(it, AAM_dyn, &
                                   f_ps, f_dps, f_u, f_theta_rhodz, f_phis)
      CALL check_conserve(f_ps, f_dps, f_u, f_theta_rhodz, f_phis, f_q, it)

      !$OMP MASTER
      CALL SYSTEM_CLOCK(stop_clock)

      IF (mpi_rank == 0) THEN
         PRINT *, "Time elapsed : ", (stop_clock - start_clock)*1./rate_clock
      END IF
      !$OMP END MASTER

   END SUBROUTINE timeloop

   SUBROUTINE print_iteration(it, itau0, itaumax, start_clock, rate_clock)
      USE omp_para, ONLY: is_master
      USE profiling_mod, ONLY: print_profile

      INTEGER :: it, itau0, itaumax, throughput
      INTEGER(kind=8) :: start_clock, stop_clock, rate_clock
      REAL :: per_step, total, elapsed

      IF (is_master) THEN
         WRITE (*, '(A,I7,A,F14.1)') "It No :", it, "   t :", dt*it
         IF (MOD(it, 10) == 0) THEN
            CALL SYSTEM_CLOCK(stop_clock)
            elapsed = (stop_clock - start_clock)*1./rate_clock
            per_step = elapsed/(it - itau0)
            throughput = INT(dt/per_step)
            total = per_step*itaumax
            WRITE (*, '(A,I5,A,F6.2,A,I6)') 'Time spent (s):', INT(elapsed), &
               '  -- ms/step : ', 1000*per_step, &
               '  -- Throughput :', throughput
            WRITE (*, '(A,I5,A,I5)') 'Whole job (min) :', INT(total/60.), &
               '  -- Completion in (min) : ', INT((total - elapsed)/60.)
         END IF
      END IF
      IF (MOD(it, itau_prof) == 0) CALL print_profile

   END SUBROUTINE print_iteration

   SUBROUTINE copy_theta_to_q(f_theta_rhodz, f_rhodz, f_q)
      TYPE(t_field), POINTER :: f_theta_rhodz(:), f_rhodz(:), f_q(:)
      REAL(rstd), POINTER :: theta_rhodz(:, :, :), rhodz(:, :), q(:, :, :)
      INTEGER :: ind, iq
      DO ind = 1, ndomain
         IF (.NOT. assigned_domain(ind)) CYCLE
         CALL swap_dimensions(ind)
         CALL swap_geometry(ind)
         theta_rhodz = f_theta_rhodz(ind)
         rhodz = f_rhodz(ind)
         q = f_q(ind)
         DO iq = 1, nqdyn
            q(:, :, iq) = theta_rhodz(:, :, iq)/rhodz(:, :)
         END DO
      END DO
   END SUBROUTINE copy_theta_to_q

   SUBROUTINE copy_q_to_theta(f_theta_rhodz, f_rhodz, f_q)
      TYPE(t_field), POINTER :: f_theta_rhodz(:), f_rhodz(:), f_q(:)
      REAL(rstd), POINTER :: theta_rhodz(:, :, :), rhodz(:, :), q(:, :, :)
      INTEGER :: ind, iq
      DO ind = 1, ndomain
         IF (.NOT. assigned_domain(ind)) CYCLE
         CALL swap_dimensions(ind)
         CALL swap_geometry(ind)
         theta_rhodz = f_theta_rhodz(ind)
         rhodz = f_rhodz(ind)
         q = f_q(ind)
         DO iq = 1, nqdyn
            theta_rhodz(:, :, iq) = rhodz(:, :)*q(:, :, iq)
         END DO
      END DO
   END SUBROUTINE copy_q_to_theta

END MODULE timeloop_gcm_mod
