How the code for the Schwarz method implementation in NEMO is organized

Implementation of the Schwarz method

The time stepping loop

The essential part to understand is how the original time stepping loop of NEMO in nemogcm.F90 is modified to allow rewinding of the state of NEMO after a coupling window has passed.

Initially the code has one loop over the command CALL stp( istp ) for each time step.

With the Schwarz algorithm, the total amount of time steps is divided in coupling windows (called Schwarz loops in the code), each coupling window is repeated the number of schwarz iterations given as a parameter (mswr), and for each schwarz iteration a number of time steps (the length of a coupling window) is done.

For a computation of 5 days with 15 time steps per day, a coupling window of 1 day and 6 Schwarz iterations there will be 3 loops.

• One global loop over the coupling windows : 5 to iterate. Counter iswloop, values 1 to nsloops.
• One inner loop over the number of schwarz iteration for each coupling window : 6. Counter kswr, values 1 to mswr.
• One most inner loop over the time steps done in one Schwarz coupling window : 15 time steps. Counter istp, values nit000 + (iswloop-1)*ntsinswr to nit000+iswloop*ntsinswr-1

The counter given to subroutine stp(istp) is rewind to its value before entering the Schwarz loop when the Schwarz coupling window is repeated.

This is coded like this in nemogcm.F90:

         iswloop = 1
         IF (lwp) WRITE(numout,*) 'init iswloop =',iswloop

         DO WHILE ( iswloop <= nsloops .AND. nstop == 0 ) ! iterate on nsloops schwarz loops

         kswr = 1
         IF (lwp) THEN 
            WRITE(numout,*) '*** schwarz loops ***'
            WRITE(numout,*) 'iswloop =',iswloop
            WRITE(numout,*) 'kswr = ',kswr
         ENDIF 
         
         CALL swz_store  ! store Ocean state at first schwarz iteration before first time step
         CALL swz_store_lim3 ! store lim3 variables state
         CALL swz_store_limdiahsb ! store limdiahsb variables
         IF(lwp) WRITE(numout,*) 'store schwarz initial state '

         DO WHILE ( kswr <= mswr .AND. nstop == 0 )

         IF (lwp) WRITE(numout,*) '  kswr = ',kswr
         istp = nit000 + (iswloop-1)*ntsinswr ! set time step to first value for iswloop (current schwarz loop)

         IF ( kswr > 1 ) CALL swz_reinit(istp) ! reset Ocean state to first time step for new schwarz iteration 

         IF (lwp) WRITE(numout,*) 'maxistp  = ',nit000+iswloop*ntsinswr-1
         DO WHILE ( istp <= nit000+iswloop*ntsinswr-1 .AND. nstop == 0 )

            IF (lwp) WRITE(numout,*) '    istp = ',istp,';    istpswz = ',istpswz 
#if defined key_agrif
            CALL stp                         ! AGRIF: time stepping
#else
            CALL stp( istp )                 ! standard time stepping
#endif
            istp = istp + 1
            istpswz = istpswz + 1
            IF( lk_mpp )   CALL mpp_max( nstop )
         END DO ! istp
         
         kswr = kswr + 1

         END DO ! kswr

         iswloop = iswloop +1

         END DO ! iswloop

The variable iswloop is the current Schwarz loop and nsloops is the number of Schwarz loops for the given time of simulation and Schwarz window size. For 5 days with a window of 1 day nsloops=5.

Inside the Schwarz repeating loop DO WHILE ( kswr <= mswr .AND. nstop == 0 ), in the time stepping loop, istp starts at the value nit000 + (iswloop-1)*ntsinswr and finishes at nit000+iswloop*ntsinswr-1.

This allows istp to take the value it would have for the current Schwarz window if no sub-iteration was done: setting mswr=1 gives you the original time stepping.

The Schwarz parameters are initialized (or computed) before entering the time loops, before and after the call to nemo_init depending on what is needed:

      ntsinswr = 1 !! set to allow fldread in nemo_init before computation of schwarz indices
      kswr = 0     !! set here to allow day_init to set nitrst properly

      WRITE(numout,*) "ntsinswr = ",ntsinswr

      !                            !-----------------------!
      CALL nemo_init               !==  Initialisations  ==!
      !                            !-----------------------!

      !! set integer constants for computation
!      ncplfrq  = 86400 !cpl_freq( 'O_QnsMix' ) read in namrun (nemoinit)
      ntsinswr = ncplfrq / INT(rdt)
      nsloops  = (nitend-nit000+1)/ntsinswr

      !! outputs to check values
      IF(lwp) THEN                  ! control print
         WRITE(numout,*)
         WRITE(numout,*) ' computation of loop indices for schwarz computation '
         WRITE(numout,*) '   nit000   = ',nit000
         WRITE(numout,*) '   nitend   = ',nitend
         WRITE(numout,*) '   rdt      = ',rdt
         WRITE(numout,*) '   ncplfrq  = ',ncplfrq
         WRITE(numout,*) '   ntsinswr = ',ntsinswr
         WRITE(numout,*) '   nsloops  = ',nsloops 
         WRITE(numout,*) '   mswr     = ',mswr
         WRITE(numout,*) '   nitrst   = ',nitrst
      ENDIF

There is also an "absolute" time step counter : istpswz which is incremented at each time step including those repeated for the Schwarz algorithm. It is needed for the coupling with OASIS.

The Schwarz algorithm: storing and restoring

The Schwarz algorithm is implemented by storing and restoring the state of the ocean and ice at the proper position inside the three loops structure.

The routines used to store and restore a state are all in schwarz.F90 which is a new module.

The storing is done just before the schwarz looping is done: inside the iswloop loop, before the kswr loop.

         CALL swz_store  ! store Ocean state at first schwarz iteration before first time step
         CALL swz_store_lim3 ! store lim3 variables state
         CALL swz_store_limdiahsb ! store limdiahsb variables
         IF(lwp) WRITE(numout,*) 'store schwarz initial state '

         DO WHILE ( kswr <= mswr .AND. nstop == 0 )

Like this we are certain to store the initial condition for a Schwarz loop.

The restoring is done at the beginning of the kswr loop before the istp loop. If we are not in the first iteration the variables state is restored:

IF ( kswr > 1 ) CALL swz_reinit(istp) ! reset Ocean state to first time step for new schwarz iteration

The routine swz_reinit is reproducing the initialisation structure from nemo_init in nemogcm.F90 by selecting only the relevant subroutines and coding the corresponding restoring of variables.

(In case you wish to use PISCES you will have to add the relevant routines in swz_reinit and module schwarz.F90)

Unfortunately this is not the only place variables are stored and restored!

Because NEMO is modular some implementation details are repeated in the modules some modules have their own storing and restoring details.

The files concerned are:

sbcmod.F90
sbcrnf.F90
traqsr.F90
trasbc.F90

Coding details: modification of other files

The routines:

iom.F90
step.F90

have been modified to output only one time serie for a set of Schwarz loops. The parameter ksout selects which Schwarz loop is given for outputs.

Some routines have a particular behavior when istp=nit000 or/and istp=nitend. They where modified to behave correctly when doing Schwarz iterations:

closea.F90
diaar5.F90
diafwb.F90
diahth.F90
dynspg_ts.F90
dynvor.F90
fldread.F90
limdyn.F90
limitd_me.F90
limrst.F90
limthd.F90
limtrp.F90
limupdate1.F90
limupdate2.F90
restart.F90
sbcblk_core.F90
sbcice_lim.F90
sbcssm.F90
stpctl.F90
traadv_eiv.F90
traadv_tvd.F90
trazdf_imp.F90
zdftke.F90
zdftmx.F90

Some other routines which are not compiled have also been modified for the sake of consistency. These are the routines relative to PISCES present in ORCA1_LIM3_PISCES/MY_SRC which get copied in ORCA2_LIM3_PISCES.

They are:

p4zflx.F90
p4zsed.F90
p4zsms.F90
trcrst.F90
trcsms_age.F90
trcsms_cfc.F90
trcstp.F90
trcwri.F90

You will need these updates if you want to add PISCES support to the Schwarz algorithm.

Lastly, but much more important, some files have been modified to declare the variables needed to store the parameters of the Schwarz algorithm and the storing/restoring variables.

The parameters are in:

in_out_manager.F90

and read in

domain.F90

The "Schwarz fields" named with _s appended at the end of the original variable name are defined in these subroutines:

dom_oce.F90
dynspg_oce.F90
ice.F90
limdiahsb.F90
oce.F90
sbc_oce.F90
sbcrnf.F90
zdf_oce.F90
zdftke.F90

a few words on coupling with OASIS

The coupling with OASIS is straightforward.

The calls to sbc_cpl_snd in step.F90 and sbc_cpl_rcv in sbcmod.F90 are modified to use the absolute counter istpswz instead of istp. In step.F90:

      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      ! Coupled mode
      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      IF( lk_oasis         )   CALL sbc_cpl_snd( istpswz )     ! coupled mode : field exchanges

In sbcmod.F90:

      CASE( jp_purecpl )  ;  CALL sbc_cpl_rcv ( istpswz, nn_fsbc, nn_ice )   ! pure coupled formulation

The only case treated is the pure coupling mode.

a few last words on testing and validation

The NEMO part of the Schwarz algorithm was validated by running a forced configuration with given atmospheric boundary conditions. The code reproduces mswr times the same day over 5 days with a coupling and Schwarz window of one day. The result was compared to the original code with no Schwarz iterations in the same configuration. The outputs are identical for the same day and any Schwarz iteration this day. Hence the storing and restoring is reliable.