= 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.