source: LMDZ6/branches/Amaury_dev/libf/dyn3dmem/leapfrog_loc.F90 @ 5134

! $Id: leapfrog_loc.F90 5134 2024-07-26 15:56:37Z abarral $

SUBROUTINE leapfrog_loc(ucov0, vcov0, teta0, ps0, &
        masse0, phis0, q0, time_0)

  USE misc_mod
  USE parallel_lmdz
  USE times
  USE mod_hallo
  USE Bands
  USE lmdz_strings, ONLY: int2str
  USE Write_Field_p
  USE lmdz_vampir
  USE lmdz_timer_filtre, ONLY: print_filtre_timer
  USE infotrac
  USE guide_loc_mod, ONLY: guide_main
  USE getparam
  USE control_mod
  USE lmdz_filtreg_p
  USE write_field_loc
  USE allocate_field_mod
  USE call_dissip_mod, ONLY: call_dissip
  USE lmdz_call_calfis, ONLY: call_calfis
  USE leapfrog_mod, ONLY: ucov, vcov, teta, ps, masse, phis, q, dq &
          , ucovm1, vcovm1, tetam1, massem1, psm1, p, pks, pk, pkf, flxw &
          , pbaru, pbarv, du, dv, dteta, phi, dp, w &
          , leapfrog_allocate, leapfrog_switch_caldyn, leapfrog_switch_dissip

  USE exner_hyb_loc_m, ONLY: exner_hyb_loc
  USE exner_milieu_loc_m, ONLY: exner_milieu_loc
  USE comconst_mod, ONLY: cpp, dtvr, ihf
  USE comvert_mod, ONLY: ap, bp, pressure_exner
  USE logic_mod, ONLY: iflag_phys, ok_guide, forward, leapf, apphys, &
          statcl, conser, apdiss, purmats, ok_strato
  USE temps_mod, ONLY: itaufin, jD_ref, jH_ref, day_ini, &
          day_ref, start_time, dt
  USE mod_xios_dyn3dmem, ONLY: dyn3d_ctx_handle
  USE lmdz_xios, ONLY: xios_update_calendar, &
          xios_set_current_context, &
          using_xios
  USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_DEBUGIO
  USE lmdz_description, ONLY: descript
  USE lmdz_iniprint, ONLY: lunout, prt_level
  USE lmdz_academic, ONLY: tetarappel, knewt_t, kfrict, knewt_g, clat4
  USE lmdz_comdissnew, ONLY: lstardis, nitergdiv, nitergrot, niterh, tetagdiv, &
          tetagrot, tetatemp, coefdis, vert_prof_dissip

  IMPLICIT NONE

  ! ......   Version  du 10/01/98    ..........

  !        avec  coordonnees  verticales hybrides
  !   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )

  !=======================================================================
  !
  !   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
  !   -------
  !
  !   Objet:
  !   ------
  !
  !   GCM LMD nouvelle grille
  !
  !=======================================================================
  !
  !  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
  !  et possibilite d'appeler une fonction f(y)  a derivee tangente
  !  hyperbolique a la  place de la fonction a derivee sinusoidale.

  !  ... Possibilite de choisir le shema pour l'advection de
  !    q  , en modifiant iadv dans traceur.def  (10/02) .
  !
  !  Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
  !  Pour Van-Leer iadv=10
  !
  !-----------------------------------------------------------------------
  !   Declarations:
  !   -------------

  INCLUDE "dimensions.h"
  INCLUDE "paramet.h"
  INCLUDE "comgeom.h"

  REAL, INTENT(IN) :: time_0 ! not used

  !   dynamical variables:
  REAL, INTENT(IN) :: ucov0(ijb_u:ije_u, llm)    ! zonal covariant wind
  REAL, INTENT(IN) :: vcov0(ijb_v:ije_v, llm)    ! meridional covariant wind
  REAL, INTENT(IN) :: teta0(ijb_u:ije_u, llm)    ! potential temperature
  REAL, INTENT(IN) :: q0(ijb_u:ije_u, llm, nqtot) ! advected tracers
  REAL, INTENT(IN) :: ps0(ijb_u:ije_u)          ! surface pressure (Pa)
  REAL, INTENT(IN) :: masse0(ijb_u:ije_u, llm)   ! air mass
  REAL, INTENT(IN) :: phis0(ijb_u:ije_u)        ! geopotentiat at the surface

  REAL :: zqmin, zqmax

  ! REAL,SAVE,ALLOCATABLE :: p (:,:  )               ! pression aux interfac.des couches
  ! REAL,SAVE,ALLOCATABLE :: pks(:)                      ! exner au  sol
  ! REAL,SAVE,ALLOCATABLE :: pk(:,:)                   ! exner au milieu des couches
  ! REAL,SAVE,ALLOCATABLE :: pkf(:,:)                  ! exner filt.au milieu des couches
  ! REAL,SAVE,ALLOCATABLE :: phi(:,:)                  ! geopotentiel
  ! REAL,SAVE,ALLOCATABLE :: w(:,:)                    ! vitesse verticale

  ! variables dynamiques intermediaire pour le transport
  ! REAL,SAVE,ALLOCATABLE :: pbaru(:,:),pbarv(:,:) !flux de masse

  !   variables dynamiques au pas -1
  ! REAL,SAVE,ALLOCATABLE :: vcovm1(:,:),ucovm1(:,:)
  !      REAL,SAVE,ALLOCATABLE :: tetam1(:,:),psm1(:)
  ! REAL,SAVE,ALLOCATABLE :: massem1(:,:)

  !   tendances dynamiques
  ! REAL,SAVE,ALLOCATABLE :: dv(:,:),du(:,:)
  ! REAL,SAVE,ALLOCATABLE :: dteta(:,:),dp(:)
  ! REAL,DIMENSION(:,:,:), ALLOCATABLE, SAVE :: dq

  !   tendances de la dissipation
  ! REAL,SAVE,ALLOCATABLE :: dvdis(:,:),dudis(:,:)
  ! REAL,SAVE,ALLOCATABLE :: dtetadis(:,:)

  !   tendances physiques
  REAL, SAVE, ALLOCATABLE :: dvfi(:, :), dufi(:, :)
  REAL, SAVE, ALLOCATABLE :: dtetafi(:, :)
  REAL, SAVE, ALLOCATABLE :: dpfi(:)
  REAL, DIMENSION(:, :, :), ALLOCATABLE, SAVE :: dqfi

  !   variables pour le fichier histoire
  REAL :: dtav      ! intervalle de temps elementaire

  REAL :: tppn(iim), tpps(iim), tpn, tps
  !
  INTEGER :: itau, itaufinp1, iav
  ! INTEGER  iday ! jour julien
  REAL :: time

  ! REAL,SAVE,ALLOCATABLE :: finvmaold(:,:)

  !ym      LOGICAL  lafin
  LOGICAL :: lafin
  INTEGER :: ij, iq, l
  INTEGER :: ik

  REAL :: time_step, t_wrt, t_ops

  ! jD_cur: jour julien courant
  ! jH_cur: heure julienne courante
  REAL :: jD_cur, jH_cur
  INTEGER :: an, mois, jour
  REAL :: secondes

  LOGICAL :: physic
  LOGICAL :: first, callinigrads

  data callinigrads/.TRUE./
  CHARACTER(LEN = 10) :: string10

  ! REAL,SAVE,ALLOCATABLE :: flxw(:,:) ! flux de masse verticale

  !+jld variables test conservation energie
  ! REAL,SAVE,ALLOCATABLE :: ecin(:,:),ecin0(:,:)
  ! Tendance de la temp. potentiel d (theta)/ d t due a la
  ! tansformation d'energie cinetique en energie thermique
  ! cree par la dissipation
  !  REAL,SAVE,ALLOCATABLE :: dtetaecdt(:,:)
  !  REAL,SAVE,ALLOCATABLE :: vcont(:,:),ucont(:,:)
  !  REAL,SAVE,ALLOCATABLE :: vnat(:,:),unat(:,:)
  REAL :: d_h_vcol, d_qt, d_qw, d_ql, d_ec
  CHARACTER(len = 15) :: ztit
  !!      INTEGER   ip_ebil_dyn  ! PRINT level for energy conserv. diag.
  ! SAVE      ip_ebil_dyn
  ! DATA      ip_ebil_dyn/0/
  !-jld

  CHARACTER(LEN = 80) :: dynhist_file, dynhistave_file
  CHARACTER(LEN = *), parameter :: modname = "leapfrog_loc"
  CHARACTER(LEN = 80) :: abort_message

  logical, PARAMETER :: dissip_conservative = .TRUE.

  INTEGER :: testita
  PARAMETER (testita = 9)

  logical, parameter :: flag_verif = .FALSE.

  ! declaration liees au parallelisme
  INTEGER :: ierr
  LOGICAL :: FirstCaldyn
  LOGICAL :: FirstPhysic
  INTEGER :: ijb, ije, j, i
  type(Request) :: TestRequest
  type(Request) :: Request_Dissip
  type(Request) :: Request_physic

  INTEGER :: true_itau
  INTEGER :: iapptrac
  INTEGER :: AdjustCount
  ! INTEGER :: var_time
  LOGICAL :: ok_start_timer = .FALSE.
  LOGICAL, SAVE :: firstcall = .TRUE.
  TYPE(distrib), SAVE :: new_dist

  CALL check_isotopes(q0, ijb_u, ije_u, 'leapfrog204: debut')

  !$OMP MASTER
  ItCount = 0
  !$OMP END MASTER
  true_itau = 0
  FirstCaldyn = .TRUE.
  FirstPhysic = .TRUE.
  iapptrac = 0
  AdjustCount = 0
  lafin = .FALSE.

  IF (nday>=0) THEN
    itaufin = nday * day_step
  else
    itaufin = -nday
  ENDIF

  itaufinp1 = itaufin + 1

  CALL check_isotopes(q0, ijb_u, ije_u, 'leapfrog 226')

  itau = 0
  physic = .TRUE.
  IF (iflag_phys==0.OR.iflag_phys==2) physic = .FALSE.
  CALL init_nan
  CALL leapfrog_allocate
  ucov = ucov0
  vcov = vcov0
  teta = teta0
  ps = ps0
  masse = masse0
  phis = phis0
  q = q0

  CALL check_isotopes(q, ijb_u, ije_u, 'leapfrog 239')

  ! iday = day_ini+itau/day_step
  ! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
  !    IF(time.GT.1.) THEN
  !     time = time-1.
  !     iday = iday+1
  !    ENDIF

  ! Allocate variables depending on dynamic variable nqtot
  !$OMP MASTER
  IF (firstcall) THEN
    ! ALLOCATE(p(ijb_u:ije_u,llmp1))
    !      ALLOCATE(pks(ijb_u:ije_u))
    ! ALLOCATE(pk(ijb_u:ije_u,llm))
    ! ALLOCATE(pkf(ijb_u:ije_u,llm))
    ! ALLOCATE(phi(ijb_u:ije_u,llm))
    ! ALLOCATE(w(ijb_u:ije_u,llm))
    ! ALLOCATE(pbaru(ip1jmp1,llm),pbarv(ip1jm,llm))
    ! ALLOCATE(vcovm1(ijb_v:ije_v,llm),ucovm1(ijb_u:ije_u,llm))
    ! ALLOCATE(tetam1(ijb_u:ije_u,llm),psm1(ijb_u:ije_u))
    ! ALLOCATE(massem1(ijb_u:ije_u,llm))
    ! ALLOCATE(dv(ijb_v:ije_v,llm),du(ijb_u:ije_u,llm))
    ! ALLOCATE(dteta(ijb_u:ije_u,llm),dp(ijb_u:ije_u))
    ! ALLOCATE(dvdis(ijb_v:ije_v,llm),dudis(ijb_u:ije_u,llm))
    ! ALLOCATE(dtetadis(ijb_u:ije_u,llm))
    ALLOCATE(dvfi(ijb_v:ije_v, llm), dufi(ijb_u:ije_u, llm))
    ALLOCATE(dtetafi(ijb_u:ije_u, llm))
    ALLOCATE(dpfi(ijb_u:ije_u))
    ! ALLOCATE(dq(ijb_u:ije_u,llm,nqtot))
    ALLOCATE(dqfi(ijb_u:ije_u, llm, nqtot))
    ! ALLOCATE(dqfi_tmp(iip1,llm,nqtot))
    ! ALLOCATE(finvmaold(ijb_u:ije_u,llm))
    ! ALLOCATE(flxw(ijb_u:ije_u,llm))
    ! ALLOCATE(ecin(ijb_u:ije_u,llm),ecin0(ijb_u:ije_u,llm))
    ! ALLOCATE(dtetaecdt(ijb_u:ije_u,llm))
    ! ALLOCATE(vcont(ijb_v:ije_v,llm),ucont(ijb_u:ije_u,llm))
    ! ALLOCATE(vnat(ijb_v:ije_v,llm),unat(ijb_u:ije_u,llm))
  ENDIF
  !$OMP END MASTER
  !$OMP BARRIER

  ! CALL dynredem1_loc("restart.nc",0.0,
  ! &                           vcov,ucov,teta,q,masse,ps)


  !-----------------------------------------------------------------------
  !   On initialise la pression et la fonction d'Exner :
  !   --------------------------------------------------

  !$OMP MASTER
  dq(:, :, :) = 0.
  CALL pression (ijnb_u, ap, bp, ps, p)
  !$OMP END MASTER
  IF (pressure_exner) THEN
    CALL exner_hyb_loc(ijnb_u, ps, p, pks, pk, pkf)
  else
    CALL exner_milieu_loc(ijnb_u, ps, p, pks, pk, pkf)
  ENDIF
  !-----------------------------------------------------------------------
  !   Debut de l'integration temporelle:
  !   ----------------------------------
  ! et du parallelisme !!

  1   CONTINUE ! Matsuno Forward step begins here

  !   date: (NB: date remains unchanged for Backward step)
  !   -----

  jD_cur = jD_ref + day_ini - day_ref + &
          (itau + 1) / day_step
  jH_cur = jH_ref + start_time + &
          mod(itau + 1, day_step) / float(day_step)
  IF (jH_cur > 1.0) THEN
    jD_cur = jD_cur + 1.
    jH_cur = jH_cur - 1.
  ENDIF

  CALL check_isotopes(q, ijb_u, ije_u, 'leapfrog 321')

  IF (ok_guide) THEN
    CALL guide_main(itau,ucov,vcov,teta,q,masse,ps)
!$OMP BARRIER
  ENDIF


  !
  ! IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
  !   CALL  test_period ( ucov,vcov,teta,q,p,phis )
  !   PRINT *,' ----   Test_period apres continue   OK ! -----', itau
  ! ENDIF
  !
  !ym      CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
  !ym      CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
  !ym      CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
  !ym      CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
  !ym      CALL SCOPY( ip1jmp1, ps  , 1,   psm1 , 1 )

  IF (FirstCaldyn) THEN
    !$OMP MASTER
    ucovm1 = ucov
    vcovm1 = vcov
    tetam1 = teta
    massem1 = masse
    psm1 = ps

    ! Ehouarn: finvmaold is actually not used
    ! finvmaold = masse
    !$OMP END MASTER
    !$OMP BARRIER
    ! CALL filtreg_p ( finvmaold ,jjb_u,jje_u,jjb_u,jje_u,jjp1, llm,
    ! &                    -2,2, .TRUE., 1 )
  else
    ! Save fields obtained at previous time step as '...m1'
    ijb = ij_begin
    ije = ij_end

    !$OMP MASTER
    psm1     (ijb:ije) = ps    (ijb:ije)
    !$OMP END MASTER

    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO l = 1, llm
      ije = ij_end
      ucovm1   (ijb:ije, l) = ucov  (ijb:ije, l)
      tetam1   (ijb:ije, l) = teta  (ijb:ije, l)
      massem1  (ijb:ije, l) = masse (ijb:ije, l)
      ! finvmaold(ijb:ije,l)=masse(ijb:ije,l)

      IF (pole_sud) ije = ij_end - iip1
      vcovm1(ijb:ije, l) = vcov  (ijb:ije, l)

    ENDDO
    !$OMP ENDDO


    ! Ehouarn: finvmaold not used
    ! CALL filtreg_p(finvmaold ,jjb_u,jje_u,jj_begin,jj_end,jjp1,
    ! .                    llm, -2,2, .TRUE., 1 )

  ENDIF ! of if (FirstCaldyn)

  forward = .TRUE.
  leapf = .FALSE.
  dt = dtvr

  !   ...    P.Le Van .26/04/94  ....

  !ym      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
  !ym      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )

  !ym  ne sert a rien
  !ym      CALL minmax(ijp1llm,q(:,:,3),zqmin,zqmax)

  CALL check_isotopes(q, ijb_u, ije_u, 'leapfrog 400')

  2   CONTINUE ! Matsuno backward or leapfrog step begins here

  CALL check_isotopes(q, ijb_u, ije_u, 'leapfrog 402')

  !$OMP MASTER
  ItCount = ItCount + 1
  IF (MOD(ItCount, 1)==1) THEN
    debug = .TRUE.
  else
    debug = .FALSE.
  ENDIF
  !$OMP END MASTER
  !-----------------------------------------------------------------------

  !   date: (NB: only leapfrog step requires recomputing date)
  !   -----

  IF (leapf) THEN
    jD_cur = jD_ref + day_ini - day_ref + &
            (itau + 1) / day_step
    jH_cur = jH_ref + start_time + &
            mod(itau + 1, day_step) / float(day_step)
    IF (jH_cur > 1.0) THEN
      jD_cur = jD_cur + 1.
      jH_cur = jH_cur - 1.
    endif
  ENDIF

  !   gestion des appels de la physique et des dissipations:
  !   ------------------------------------------------------
  !
  !   ...    P.Le Van  ( 6/02/95 )  ....

  apphys = .FALSE.
  statcl = .FALSE.
  conser = .FALSE.
  apdiss = .FALSE.

  IF(purmats) THEN
    ! Purely Matsuno time stepping
    IF(MOD(itau, iconser) ==0.AND.  forward) conser = .TRUE.
    IF(MOD(itau, dissip_period)==0.AND..NOT.forward) &
            apdiss = .TRUE.
    IF(MOD(itau, iphysiq)==0.AND..NOT.forward &
            .AND. physic) apphys = .TRUE.
  ELSE
    ! Leapfrog/Matsuno time stepping
    IF(MOD(itau, iconser) == 0) conser = .TRUE.
    IF(MOD(itau + 1, dissip_period)==0 .AND. .NOT. forward) &
            apdiss = .TRUE.
    IF(MOD(itau + 1, iphysiq)==0.AND.physic) apphys = .TRUE.
  END IF

  ! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
  ! supress dissipation step
  IF (llm==1) THEN
    apdiss = .FALSE.
  ENDIF

  !ym    ---> Pour le moment
  !ym      apphys = .FALSE.
  statcl = .FALSE.
  ! conser = .FALSE. ! ie: no output of control variables to stdout in //

  IF (firstCaldyn) THEN
    !$OMP MASTER
    CALL Set_Distrib(distrib_caldyn)
    !$OMP END MASTER
    !$OMP BARRIER
    firstCaldyn = .FALSE.
    !ym          CALL InitTime
    !$OMP MASTER
    CALL Init_timer
    !$OMP END MASTER
  ENDIF

  !$OMP MASTER
  IF (ok_start_timer) THEN
    CALL InitTime
    ok_start_timer = .FALSE.
  ENDIF
  !$OMP END MASTER

  CALL check_isotopes(q, ijb_u, ije_u, 'leapfrog 471')

  !ym  PAS D'AJUSTEMENT POUR LE MOMENT
  IF (Adjust) THEN
    AdjustCount = AdjustCount + 1
    ! if (iapptrac==iapp_tracvl .AND. (forward .OR.  leapf)
    ! &         .AND. itau/iphysiq>2 .AND. Adjustcount>30) THEN
    IF (Adjustcount>1) THEN
      AdjustCount = 0
      !$OMP MASTER
      CALL allgather_timer_average

      IF (prt_level > 9) THEN
        print *, '*********************************'
        print *, '******    TIMER CALDYN     ******'
        do i = 0, mpi_size - 1
          print *, 'proc', i, ' :   Nb Bandes  :', jj_nb_caldyn(i), &
                  '  : temps moyen :', &
                  timer_average(jj_nb_caldyn(i), timer_caldyn, i), &
                  '+-', timer_delta(jj_nb_caldyn(i), timer_caldyn, i)
        enddo

        print *, '*********************************'
        print *, '******    TIMER VANLEER    ******'
        do i = 0, mpi_size - 1
          print *, 'proc', i, ' :   Nb Bandes  :', jj_nb_vanleer(i), &
                  '  : temps moyen :', &
                  timer_average(jj_nb_vanleer(i), timer_vanleer, i), &
                  '+-', timer_delta(jj_nb_vanleer(i), timer_vanleer, i)
        enddo

        print *, '*********************************'
        print *, '******    TIMER DISSIP    ******'
        do i = 0, mpi_size - 1
          print *, 'proc', i, ' :   Nb Bandes  :', jj_nb_dissip(i), &
                  '  : temps moyen :', &
                  timer_average(jj_nb_dissip(i), timer_dissip, i), &
                  '+-', timer_delta(jj_nb_dissip(i), timer_dissip, i)
        enddo

        ! if (mpi_rank==0) CALL WriteBands

      endif

      CALL AdjustBands_caldyn(new_dist)
      !$OMP END MASTER
      !$OMP BARRIER
      CALL leapfrog_switch_caldyn(new_dist)
      !$OMP BARRIER


      !$OMP MASTER
      distrib_caldyn = new_dist
      CALL set_distrib(distrib_caldyn)
      !$OMP END MASTER
      !$OMP BARRIER
      ! CALL Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
      ! &                                jj_Nb_caldyn,0,0,TestRequest)
      !     CALL Register_SwapFieldHallo(ucovm1,ucovm1,ip1jmp1,llm,
      ! &                                jj_Nb_caldyn,0,0,TestRequest)
      !     CALL Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
      ! &                                jj_Nb_caldyn,0,0,TestRequest)
      !     CALL Register_SwapFieldHallo(vcovm1,vcovm1,ip1jm,llm,
      ! &                                jj_Nb_caldyn,0,0,TestRequest)
      !     CALL Register_SwapFieldHallo(teta,teta,ip1jmp1,llm,
      ! &                                jj_Nb_caldyn,0,0,TestRequest)
      !     CALL Register_SwapFieldHallo(tetam1,tetam1,ip1jmp1,llm,
      ! &                                jj_Nb_caldyn,0,0,TestRequest)
      !     CALL Register_SwapFieldHallo(masse,masse,ip1jmp1,llm,
      ! &                                jj_Nb_caldyn,0,0,TestRequest)
      !     CALL Register_SwapFieldHallo(massem1,massem1,ip1jmp1,llm,
      ! &                                jj_Nb_caldyn,0,0,TestRequest)
      !     CALL Register_SwapFieldHallo(ps,ps,ip1jmp1,1,
      ! &                                jj_Nb_caldyn,0,0,TestRequest)
      !     CALL Register_SwapFieldHallo(psm1,psm1,ip1jmp1,1,
      ! &                                jj_Nb_caldyn,0,0,TestRequest)
      !     CALL Register_SwapFieldHallo(pkf,pkf,ip1jmp1,llm,
      ! &                                jj_Nb_caldyn,0,0,TestRequest)
      !     CALL Register_SwapFieldHallo(pk,pk,ip1jmp1,llm,
      ! &                                jj_Nb_caldyn,0,0,TestRequest)
      !     CALL Register_SwapFieldHallo(pks,pks,ip1jmp1,1,
      ! &                                jj_Nb_caldyn,0,0,TestRequest)
      !     CALL Register_SwapFieldHallo(phis,phis,ip1jmp1,1,
      ! &                                jj_Nb_caldyn,0,0,TestRequest)
      !     CALL Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
      ! &                                jj_Nb_caldyn,0,0,TestRequest)
      !     CALL Register_SwapFieldHallo(finvmaold,finvmaold,ip1jmp1,llm,
      ! &                                jj_Nb_caldyn,0,0,TestRequest)

      !    do j=1,nqtot
      !     CALL Register_SwapFieldHallo(q(:,:,j),q(:,:,j),ip1jmp1,llm,
      ! &                                jj_nb_caldyn,0,0,TestRequest)
      !    enddo

      !     CALL Set_Distrib(distrib_caldyn)
      !     CALL SendRequest(TestRequest)
      !     CALL WaitRequest(TestRequest)

      !$OMP MASTER
      CALL AdjustBands_dissip(new_dist)
      !$OMP END MASTER
      !$OMP BARRIER
      CALL leapfrog_switch_dissip(new_dist)
      !$OMP BARRIER
      !$OMP MASTER
      distrib_dissip = new_dist
      !$OMP END MASTER
      !$OMP BARRIER
      ! CALL AdjustBands_physic

      !$OMP MASTER
      IF (mpi_rank==0) CALL WriteBands
      !$OMP END MASTER

    endif
  ENDIF

  CALL check_isotopes(q, ijb_u, ije_u, 'leapfrog 589')

  !-----------------------------------------------------------------------
  !   calcul des tendances dynamiques:
  !   --------------------------------
  !$OMP BARRIER
  !$OMP MASTER
  CALL VTb(VThallo)
  !$OMP END MASTER

  CALL Register_Hallo_u(ucov, llm, 1, 1, 1, 1, TestRequest)
  CALL Register_Hallo_v(vcov, llm, 1, 1, 1, 1, TestRequest)
  CALL Register_Hallo_u(teta, llm, 1, 1, 1, 1, TestRequest)
  CALL Register_Hallo_u(ps, 1, 1, 2, 2, 1, TestRequest)
  CALL Register_Hallo_u(pkf, llm, 1, 1, 1, 1, TestRequest)
  CALL Register_Hallo_u(pk, llm, 1, 1, 1, 1, TestRequest)
  CALL Register_Hallo_u(pks, 1, 1, 1, 1, 1, TestRequest)
  CALL Register_Hallo_u(p, llmp1, 1, 1, 1, 1, TestRequest)

  ! do j=1,nqtot
  !   CALL Register_Hallo(q(1,1,j),ip1jmp1,llm,1,1,1,1,
  ! *                       TestRequest)
  !    enddo

  CALL SendRequest(TestRequest)
  !$OMP BARRIER
  CALL WaitRequest(TestRequest)

  !$OMP MASTER
  CALL VTe(VThallo)
  !$OMP END MASTER
  !$OMP BARRIER

  IF (debug) THEN
    CALL WriteField_u('ucov', ucov)
    CALL WriteField_v('vcov', vcov)
    CALL WriteField_u('teta', teta)
    CALL WriteField_u('ps', ps)
    CALL WriteField_u('masse', masse)
    CALL WriteField_u('pk', pk)
    CALL WriteField_u('pks', pks)
    CALL WriteField_u('pkf', pkf)
    CALL WriteField_u('phis', phis)
    do iq = 1, nqtot
      CALL WriteField_u('q' // trim(int2str(iq)), &
              q(:, :, iq))
    enddo
  ENDIF

  True_itau = True_itau + 1

  !$OMP MASTER
  IF (prt_level>9) THEN
    WRITE(lunout, *)"leapfrog_p: Iteration No", True_itau
  ENDIF

  CALL start_timer(timer_caldyn)

  ! compute geopotential phi()
  CALL geopot_loc  (ip1jmp1, teta, pk, pks, phis, phi)

  CALL check_isotopes(q, ijb_u, ije_u, 'leapfrog 651')

  CALL VTb(VTcaldyn)
  !$OMP END MASTER
  ! var_time=time+iday-day_ini

  !$OMP BARRIER
  ! CALL FTRACE_REGION_BEGIN("caldyn")
  time = jD_cur + jH_cur

  CALL caldyn_loc &
          (itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, &
          phi, conser, du, dv, dteta, dp, w, pbaru, pbarv, time)

  ! CALL FTRACE_REGION_END("caldyn")

  !$OMP MASTER
  IF (mpi_rank==0.AND.conser) THEN
    WRITE(lunout, *) 'leapfrog_loc, Time step: ', itau, ' Day:', time
  ENDIF
  CALL VTe(VTcaldyn)
  !$OMP END MASTER

  IF (CPPKEY_DEBUGIO) THEN
    CALL WriteField_u('du', du)
    CALL WriteField_v('dv', dv)
    CALL WriteField_u('dteta', dteta)
    CALL WriteField_u('dp', dp)
    CALL WriteField_u('w', w)
    CALL WriteField_u('pbaru', pbaru)
    CALL WriteField_v('pbarv', pbarv)
    CALL WriteField_u('p', p)
    CALL WriteField_u('masse', masse)
    CALL WriteField_u('pk', pk)
  END IF
  !-----------------------------------------------------------------------
  !   calcul des tendances advection des traceurs (dont l'humidite)
  !   -------------------------------------------------------------

  CALL check_isotopes(q, ijb_u, ije_u, &
          'leapfrog 686: avant caladvtrac')

  IF(forward .OR.  leapf)  THEN
    ! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step
    !WRITE(*,*) 'leapfrog 679: avant CALL caladvtrac_loc'
    CALL caladvtrac_loc(q, pbaru, pbarv, &
            p, masse, dq, teta, &
            flxw, pk, iapptrac)

    ! CALL creation of mass flux
    IF (offline .AND. .NOT. adjust) THEN
      CALL fluxstokenc_p(pbaru, pbarv, masse, teta, phi)
    ENDIF

    !WRITE(*,*) 'leapfrog 719'
    CALL check_isotopes(q, ijb_u, ije_u, &
            'leapfrog 698: apres caladvtrac')

    ! do j=1,nqtot
    !   CALL WriteField_u('qadv'//trim(int2str(j)),q(:,:,j))
    ! enddo

    ! Ehouarn: Storage of mass flux for off-line tracers... not implemented...

  ENDIF ! of IF( forward .OR.  leapf )


  !-----------------------------------------------------------------------
  !   integrations dynamique et traceurs:
  !   ----------------------------------

  !$OMP MASTER
  CALL VTb(VTintegre)
  !$OMP END MASTER
  IF (CPPKEY_DEBUGIO) THEN
    IF (true_itau>20) THEN
      CALL WriteField_u('ucovm1', ucovm1)
      CALL WriteField_v('vcovm1', vcovm1)
      CALL WriteField_u('tetam1', tetam1)
      CALL WriteField_u('psm1', psm1)
      CALL WriteField_u('ucov_int', ucov)
      CALL WriteField_v('vcov_int', vcov)
      CALL WriteField_u('teta_int', teta)
      CALL WriteField_u('ps_int', ps)
    endif
  END IF
  !$OMP BARRIER
  ! CALL FTRACE_REGION_BEGIN("integrd")

  !WRITE(*,*) 'leapfrog 720'
  CALL check_isotopes(q, ijb_u, ije_u, 'leapfrog 756')

  ! CRisi: pourquoi aller jusqu'à 2 et non pas jusqu'à nqtot??
  CALL integrd_loc (nqtot, vcovm1, ucovm1, tetam1, psm1, massem1, &
          dv, du, dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis)
  ! $              finvmaold                                    )

  !  !WRITE(*,*) 'leapfrog 724'
  CALL check_isotopes(q, ijb_u, ije_u, 'leapfrog 762')

  ! CALL FTRACE_REGION_END("integrd")
  !$OMP BARRIER
  IF (CPPKEY_DEBUGIO) THEN
    CALL WriteField_u('ucovm1', ucovm1)
    CALL WriteField_v('vcovm1', vcovm1)
    CALL WriteField_u('tetam1', tetam1)
    CALL WriteField_u('psm1', psm1)
    CALL WriteField_u('ucov_int', ucov)
    CALL WriteField_v('vcov_int', vcov)
    CALL WriteField_u('teta_int', teta)
    CALL WriteField_u('ps_int', ps)
  END IF

  CALL check_isotopes(q, ijb_u, ije_u, 'leapfrog 775')

  ! do j=1,nqtot
  !   CALL WriteField_p('q'//trim(int2str(j)),
  ! .                reshape(q(:,:,j),(/iip1,jmp1,llm/)))
  !    CALL WriteField_p('dq'//trim(int2str(j)),
  ! .                reshape(dq(:,:,j),(/iip1,jmp1,llm/)))
  !  enddo


  !$OMP MASTER
  CALL VTe(VTintegre)
  !$OMP END MASTER
  ! .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
  !
  !-----------------------------------------------------------------------
  !   calcul des tendances physiques:
  !   -------------------------------
  !    ########   P.Le Van ( Modif le  6/02/95 )   ###########
  !
  IF(purmats)  THEN
    IF(itau==itaufin.AND..NOT.forward) lafin = .TRUE.
  ELSE
    IF(itau + 1 == itaufin)              lafin = .TRUE.
  ENDIF

  !c$OMP END PARALLEL

  !
  !
  IF(apphys)  THEN

    CALL call_calfis(itau, lafin, ucov, vcov, teta, masse, ps, phis, q, flxw)

    ! c-jld
    !$OMP MASTER
    IF (FirstPhysic) THEN
      ok_start_timer = .TRUE.
      FirstPhysic = .FALSE.
    endif
    !$OMP END MASTER
  ENDIF ! of IF( apphys )

  CALL check_isotopes(q, ijb_u, ije_u, 'leapfrog 1132')
  !WRITE(*,*) 'leapfrog 1134: iflag_phys=',iflag_phys

  IF(iflag_phys==2) THEN ! "Newtonian" case
    !$OMP MASTER
    IF (FirstPhysic) THEN
      ok_start_timer = .TRUE.
      FirstPhysic = .FALSE.
    endif
    !$OMP END MASTER


    !   Calcul academique de la physique = Rappel Newtonien + fritcion
    !   --------------------------------------------------------------
    !ym       teta(:,:)=teta(:,:)
    !ym     s  -iphysiq*dtvr*(teta(:,:)-tetarappel(:,:))/taurappel
    ijb = ij_begin
    ije = ij_end
    !LF       teta(ijb:ije,:)=teta(ijb:ije,:)
    !LF     s  -iphysiq*dtvr*(teta(ijb:ije,:)-tetarappel(ijb:ije,:))/taurappel
    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    do l = 1, llm
      teta(ijb:ije, l) = teta(ijb:ije, l) - dtvr * &
              (teta(ijb:ije, l) - tetarappel(ijb:ije, l)) * &
              (knewt_g + knewt_t(l) * clat4(ijb:ije))
    enddo
    !$OMP END DO

    !$OMP MASTER
    IF (planet_type=="giant") THEN
      ! add an intrinsic heat flux at the base of the atmosphere
      teta(ijb:ije, 1) = teta(ijb:ije, 1) &
              + dtvr * aire(ijb:ije) * ihf / cpp / masse(ijb:ije, 1)
    endif
    !$OMP END MASTER
    !$OMP BARRIER

    CALL Register_Hallo_u(ucov, llm, 0, 1, 1, 0, Request_Physic)
    CALL Register_Hallo_v(vcov, llm, 1, 1, 1, 1, Request_Physic)
    CALL SendRequest(Request_Physic)
    !$OMP BARRIER
    CALL WaitRequest(Request_Physic)
    !$OMP BARRIER
    CALL friction_loc(ucov, vcov, dtvr)
    !$OMP BARRIER

    ! Sponge layer (if any)
    IF (ok_strato) THEN
      CALL top_bound_loc(vcov, ucov, teta, masse, dtvr)
      !$OMP BARRIER
    ENDIF ! of IF (ok_strato)
  ENDIF ! of IF(iflag_phys.EQ.2)

  CALL pression_loc (ip1jmp1, ap, bp, ps, p)
  !$OMP BARRIER
  IF (pressure_exner) THEN
    CALL exner_hyb_loc(ijnb_u, ps, p, pks, pk, pkf)
  else
    CALL exner_milieu_loc(ijnb_u, ps, p, pks, pk, pkf)
  ENDIF
  !$OMP BARRIER
  CALL massdair_loc(p, masse)
  !$OMP BARRIER

  !c$OMP END PARALLEL
  CALL check_isotopes(q, ijb_u, ije_u, 'leapfrog 1196')

  !-----------------------------------------------------------------------
  !   dissipation horizontale et verticale  des petites echelles:
  !   ----------------------------------------------------------
  IF(apdiss) THEN
    CALL call_dissip(ucov, vcov, teta, p, pk, ps)
  END IF

  CALL check_isotopes(q, ijb_u, ije_u, 'leapfrog 1430')

  !   ********************************************************************
  !   .... fin de l'integration dynamique  et physique pour le pas itau ..
  !   ********************************************************************

  !   preparation du pas d'integration suivant  ......
  !$OMP MASTER
  CALL stop_timer(timer_caldyn)
  !$OMP END MASTER
  IF (itau==itaumax) THEN
    !$OMP MASTER
    CALL allgather_timer_average
    CALL barrier
    IF (mpi_rank==0) THEN
      print *, '*********************************'
      print *, '******    TIMER CALDYN     ******'
      do i = 0, mpi_size - 1
        print *, 'proc', i, ' :   Nb Bandes  :', jj_nb_caldyn(i), &
                '  : temps moyen :', &
                timer_average(jj_nb_caldyn(i), timer_caldyn, i)
      enddo

      print *, '*********************************'
      print *, '******    TIMER VANLEER    ******'
      do i = 0, mpi_size - 1
        print *, 'proc', i, ' :   Nb Bandes  :', jj_nb_vanleer(i), &
                '  : temps moyen :', &
                timer_average(jj_nb_vanleer(i), timer_vanleer, i)
      enddo

      print *, '*********************************'
      print *, '******    TIMER DISSIP    ******'
      do i = 0, mpi_size - 1
        print *, 'proc', i, ' :   Nb Bandes  :', jj_nb_dissip(i), &
                '  : temps moyen :', &
                timer_average(jj_nb_dissip(i), timer_dissip, i)
      enddo

      print *, '*********************************'
      print *, '******    TIMER PHYSIC    ******'
      do i = 0, mpi_size - 1
        print *, 'proc', i, ' :   Nb Bandes  :', jj_nb_physic(i), &
                '  : temps moyen :', &
                timer_average(jj_nb_physic(i), timer_physic, i)
      enddo

    endif
    CALL barrier
    print *, 'Taille du Buffer MPI (REAL*8)', MaxBufferSize
    print *, 'Taille du Buffer MPI utilise (REAL*8)', MaxBufferSize_Used
    print *, 'Temps total ecoule sur la parallelisation :', DiffTime()
    print *, 'Temps CPU ecoule sur la parallelisation :', DiffCpuTime()
    CALL print_filtre_timer
    !$OMP END MASTER
    CALL dynredem1_loc("restart.nc", 0.0, &
            vcov, ucov, teta, q, masse, ps)
    !$OMP MASTER
    CALL fin_getparam
    !$OMP END MASTER

    IF (ok_guide) THEN
      ! set ok_guide to false to avoid extra output
      ! in following forward step
      ok_guide = .FALSE.
    endif

      IF (ANY(type_trac == ['inca', 'inco'])) THEN
        CALL finalize_inca
        ! switching back to LMDZDYN context
        !$OMP MASTER
        IF (ok_dyn_xios) THEN
          CALL xios_set_current_context(dyn3d_ctx_handle)
        ENDIF
        !$OMP END MASTER
      ENDIF
#ifdef REPROBUS
     IF (type_trac == 'repr') CALL finalize_reprobus
#endif

    !$OMP MASTER
    CALL finalize_parallel
    !$OMP END MASTER
    !$OMP BARRIER
    RETURN
  ENDIF

  CALL check_isotopes(q, ijb_u, ije_u, 'leapfrog 1509')

  IF (.NOT.purmats) THEN
    ! ........................................................
    ! ..............  schema matsuno + leapfrog  ..............
    ! ........................................................

    IF(forward .OR. leapf) THEN
      itau = itau + 1
      ! iday= day_ini+itau/day_step
      ! time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
      !   IF(time.GT.1.) THEN
      !     time = time-1.
      !     iday = iday+1
      !   ENDIF
    ENDIF

    IF(itau == itaufinp1) THEN
      IF (flag_verif) THEN
        WRITE(79, *) 'ucov', ucov
        WRITE(80, *) 'vcov', vcov
        WRITE(81, *) 'teta', teta
        WRITE(82, *) 'ps', ps
        WRITE(83, *) 'q', q
        WRITE(85, *) 'q1 = ', q(:, :, 1)
        WRITE(86, *) 'q3 = ', q(:, :, 3)
      endif


      !$OMP MASTER
      CALL fin_getparam
      !$OMP END MASTER

        IF (ANY(type_trac == ['inca', 'inco'])) THEN
          CALL finalize_inca
          ! switching back to LMDZDYN context
          !$OMP MASTER
          IF (ok_dyn_xios) THEN
            CALL xios_set_current_context(dyn3d_ctx_handle)
          ENDIF
          !$OMP END MASTER
        ENDIF
#ifdef REPROBUS
          IF (type_trac == 'repr') CALL finalize_reprobus
#endif

      !$OMP MASTER
      CALL finalize_parallel
      !$OMP END MASTER
      abort_message = 'Simulation finished'
      CALL abort_gcm(modname, abort_message, 0)
      RETURN
    ENDIF
    !-----------------------------------------------------------------------
    !   ecriture du fichier histoire moyenne:
    !   -------------------------------------

    IF(MOD(itau, iperiod)==0 .OR. itau==itaufin) THEN
      !$OMP BARRIER
      IF(itau==itaufin) THEN
        iav = 1
      ELSE
        iav = 0
      ENDIF

      ! Ehouarn: re-compute geopotential for outputs
      !$OMP BARRIER
      !$OMP MASTER
      CALL geopot_loc(ip1jmp1, teta, pk, pks, phis, phi)
      !$OMP END MASTER
      !$OMP BARRIER

         IF (ok_dynzon) THEN

          CALL bilan_dyn_loc(2,dtvr*iperiod,dtvr*day_step*periodav, &
                ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)

          ENDIF !ok_dynzon

          IF (ok_dyn_ave) THEN
             CALL writedynav_loc(itau,vcov, &
                   ucov,teta,pk,phi,q,masse,ps,phis)
          ENDIF

    ENDIF

    CALL check_isotopes(q, ijb_u, ije_u, 'leapfrog 1584')

    !-----------------------------------------------------------------------
    !   ecriture de la bande histoire:
    !   ------------------------------

    IF(MOD(itau, iecri)==0) THEN
      ! Ehouarn: output only during LF or Backward Matsuno
      IF (leapf.OR.(.NOT.leapf.AND.(.NOT.forward))) THEN
        !$OMP BARRIER
        !$OMP MASTER
        CALL geopot_loc(ip1jmp1, teta, pk, pks, phis, phi)
        !$OMP END MASTER
        !$OMP BARRIER

         IF (ok_dyn_ins) THEN
             CALL writehist_loc(itau,vcov,ucov,teta,pk,phi,q, &
                   masse,ps,phis)
         endif

        IF (ok_dyn_xios) THEN
          !$OMP MASTER
          CALL xios_update_calendar(itau)
          !$OMP END MASTER
          !$OMP BARRIER
          CALL writedyn_xios(vcov, &
                  ucov, teta, pk, phi, q, masse, ps, phis)
        ENDIF

      endif                 ! of if (leapf.OR.(.NOT.leapf.AND.(.NOT.forward)))

    ENDIF ! of IF(MOD(itau,iecri).EQ.0)

    IF(itau==itaufin) THEN

      !$OMP BARRIER

      ! if (planet_type.EQ."earth") THEN
      ! Write an Earth-format restart file
      CALL dynredem1_loc("restart.nc", 0.0, &
              vcov, ucov, teta, q, masse, ps)
      ! END IF ! of if (planet_type.EQ."earth")
      IF (ok_guide) THEN
        ! set ok_guide to false to avoid extra output
        ! in following forward step
        ok_guide = .FALSE.
      endif

      ! CLOSE(99)
    ENDIF ! of IF (itau.EQ.itaufin)

    CALL check_isotopes(q, ijb_u, ije_u, 'leapfrog 1624')

    !-----------------------------------------------------------------------
    !   gestion de l'integration temporelle:
    !   ------------------------------------

    IF(MOD(itau, iperiod)==0)    THEN
      GO TO 1
    ELSE IF (MOD(itau - 1, iperiod) == 0) THEN

      IF(forward)  THEN
        ! fin du pas forward et debut du pas backward

        forward = .FALSE.
        leapf = .FALSE.
        GO TO 2

      ELSE
        ! fin du pas backward et debut du premier pas leapfrog

        leapf = .TRUE.
        dt = 2. * dtvr
        GO TO 2
      END IF
    ELSE

      ! ......   pas leapfrog  .....

      leapf = .TRUE.
      dt = 2. * dtvr
      GO TO 2
    END IF ! of IF (MOD(itau,iperiod).EQ.0)
    !    ELSEIF (MOD(itau-1,iperiod).EQ.0)

  ELSE ! of IF (.NOT.purmats)

    CALL check_isotopes(q, ijb_u, ije_u, 'leapfrog 1664')

    ! ........................................................
    ! ..............       schema  matsuno        ...............
    ! ........................................................
    IF(forward)  THEN

      itau = itau + 1
      ! iday = day_ini+itau/day_step
      ! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0

      !      IF(time.GT.1.) THEN
      !       time = time-1.
      !       iday = iday+1
      !      ENDIF

      forward = .FALSE.
      IF(itau == itaufinp1) THEN
        !$OMP MASTER
        CALL fin_getparam
        !$OMP END MASTER

          IF (ANY(type_trac == ['inca', 'inco'])) THEN
            CALL finalize_inca
            ! switching back to LMDZDYN context
            !$OMP MASTER
            IF (ok_dyn_xios) THEN
              CALL xios_set_current_context(dyn3d_ctx_handle)
            ENDIF
            !$OMP END MASTER
          ENDIF
#ifdef REPROBUS
             IF (type_trac == 'repr') CALL finalize_reprobus
#endif

        !$OMP MASTER
        CALL finalize_parallel
        !$OMP END MASTER
        abort_message = 'Simulation finished'
        CALL abort_gcm(modname, abort_message, 0)
        RETURN
      ENDIF
      GO TO 2

    ELSE ! of IF(forward) i.e. backward step

      CALL check_isotopes(q, ijb_u, ije_u, 'leapfrog 1698')

      IF(MOD(itau, iperiod)==0 .OR. itau==itaufin) THEN
        IF(itau==itaufin) THEN
          iav = 1
        ELSE
          iav = 0
        ENDIF

          ! Ehouarn: re-compute geopotential for outputs
!$OMP BARRIER
!$OMP MASTER
          CALL geopot_loc(ip1jmp1,teta,pk,pks,phis,phi)
!$OMP END MASTER
!$OMP BARRIER

           IF (ok_dynzon) THEN
           CALL bilan_dyn_loc(2,dtvr*iperiod,dtvr*day_step*periodav, &
                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
           ENDIF

           IF (ok_dyn_ave) THEN
             CALL writedynav_loc(itau,vcov, &
                   ucov,teta,pk,phi,q,masse,ps,phis)
           ENDIF

      ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)

      IF(MOD(itau, iecri)==0) THEN

        !$OMP BARRIER
        !$OMP MASTER
        CALL geopot_loc(ip1jmp1, teta, pk, pks, phis, phi)
        !$OMP END MASTER
        !$OMP BARRIER


          IF (ok_dyn_ins) THEN
             CALL writehist_loc(itau,vcov,ucov,teta,pk,phi,q, &
                   masse,ps,phis)
          endif ! of if (ok_dyn_ins)

        IF (ok_dyn_xios) THEN
          !$OMP MASTER
          CALL xios_update_calendar(itau)
          !$OMP END MASTER
          !$OMP BARRIER
          CALL writedyn_xios(vcov, &
                  ucov, teta, pk, phi, q, masse, ps, phis)
        ENDIF

      ENDIF                ! of IF(MOD(itau,iecri).EQ.0)

      IF(itau==itaufin) THEN
        ! if (planet_type.EQ."earth") THEN
        CALL dynredem1_loc("restart.nc", 0.0, &
                vcov, ucov, teta, q, masse, ps)
        ! END IF ! of if (planet_type.EQ."earth")
        IF (ok_guide) THEN
          ! set ok_guide to false to avoid extra output
          ! in following forward step
          ok_guide = .FALSE.
        endif

      ENDIF ! of IF(itau.EQ.itaufin)

      forward = .TRUE.
      GO TO  1

    ENDIF ! of IF (forward)

    CALL check_isotopes(q, ijb_u, ije_u, 'leapfrog 1750')

  END IF ! of IF(.NOT.purmats)
  !$OMP MASTER
  CALL fin_getparam
  !$OMP END MASTER

    IF (ANY(type_trac == ['inca', 'inco'])) THEN
      CALL finalize_inca
      ! switching back to LMDZDYN context
      !$OMP MASTER
      IF (ok_dyn_xios) THEN
        CALL xios_set_current_context(dyn3d_ctx_handle)
      ENDIF
      !$OMP END MASTER
    ENDIF
#ifdef REPROBUS
  IF (type_trac == 'repr') CALL finalize_reprobus
#endif

  !$OMP MASTER
  CALL finalize_parallel
  !$OMP END MASTER
  abort_message = 'Simulation finished'
  CALL abort_gcm(modname, abort_message, 0)

END SUBROUTINE leapfrog_loc