Documentation on d4 developments in ALADIN NH 
Radmila Brozkova (01/08/2005)

Summary
The aim of these developments was two-fold: i) to clean out all four sub-variants of d4 implementation (option ND4SYS): these subvariants handled artificial residuals of the equality d4 = d3 + X as regards the time-stepping scheme; ii) to optimise the implementation of d4. Developments were made in three consecutive steps:

1. To clean the ND4SYS option and using strictly d4 = d3 + X at all (x,y,eta,t).

2. To suppress the YAUX structure for the auxilary vertical divergence variable and introduce it as an additional thermodynamical variable NFTHER in a standard GMV structure.

3. To remove the biperiodicisation step on the auxilary variable at every time step and ensure the biperiodicity of this field by lateral coupling as it is done for other variables.

Definitions, naming convention and data flow
The variable d4 is defined as d3 + X, where d3 is vertical divergence g * dw/d_phi (caution, d_phi is "dry" here) and X is the horizontal wind-shear part (X=dV/d_phi dot grad_phi) of the total three dimensional divergence.

When using d4 variable we treat the evolution of its two parts d3 and X in a different way. While the material derivative D (d3)/Dt can be treated by a classical approach, D (X)/Dt is approximated simply by a first order scheme using values of X at the time levels t and/or t-dt pending the time scheme (see ALADIN NH documentation for details).
This however makes an important difference: X ougth to be freshly re-diagnosed within each time-step for time levels t and/or t-dt to restore a reasonable precision. The precision, as well as impliciteness, is naturally enhanced when an iterative time scheme is used (the predictor-corrector scheme in a short jargon).

Further, when advecting d4 variable (main option), there is a need to compute the horizontal derivatives of X in order to compute the horizontal derivatives of w (vertical velocity Dz/Dt) from relationships between d3 and w. These w-derivatives are needed to evaluate a so-called "Z" term present at the right-hand-side of d4 evolution equation. Since we cannot diagnose the derivatives of X technically as easy as X itself, we get them from the approximated X at t+dt computed at the previous time-step of the model. This step requires to transform the estimate of X at t+dt to the spectral space, where no other computations with X are made except the horizontal diffusion, and within the inverse transforms to obtain the derivatives of X.

The above described step is not needed when advecting w (LGWADV option) since in such a case we do not need to compute the "Z" term. This option has however also some drawbacks and is very likely unusable in the Eulerian advection scheme together with the current vertical staggering.

The particular treatment of X has a straight impact on the data flow:
- At a given time level X is systematically diagnosed, therefore it is not present in the history data files.
- Provided the "Z" term is not computed (LGWADV), there is only the diagnostics of X made within the grid-point computations of a given time-step (PNHXT0 and/or PNHXT9 computed at the CPG_GP level). When the time-scheme iterations are used (PC scheme), we need to transfer the quantity of X between the predictor and correctors calls. This is ensured by the GMV array of YT9, denoted as PGMVNHXT9 to distinguish from the true X at t-dt (PNHXT9).
- In a nominal configuration we compute the Z term, implying the spectral treatment of X. This means to allocate spectral array SPNHX, done now within the thermodynamical variables category NFTHER. At grid-point space the corresponding GMV of T0 and T1 type structures include X. It means as well to perform the lateral coupling of X in case of ALADIN.

To recall the X variable, we do have:
SPNHX spectral array within SPHV: allocated only when (NVDVAR==4 .AND. .NOT. LGWADV);
GMV of YT9%MNHX allocated when  NVDVAR==4 for both 2TL and 3TL schemes;
GMV of YT9%MNHXNL and YT9%MCNHXNL when NVDVAR==4 and Eulerian advection to account for the Asselin filter when PC scheme. It is used only at the level of CPEULDYN;
GMV of YT0%MNHX and derivatives YT0%MNHXL and YT0%MNHXM when (NVDVAR==4 .AND. .NOT. LGWADV);
GMV of YT1%MNHX when (NVDVAR==4 .AND. .NOT. LGWADV).
GMVCPL: this array serves for the large-scale coupling (ALADIN) and it has a same structure like GMV of YT0. Automatically the X variable is available at GMVCPL when allocated within GMV of YT0.

Locally at GP_MODEL level we do have:
PGMVNHXT9 corresponding to GMV of YT9%MNHX. It is used only as a storage to carry X from predictor to corrector within one time-step (2TL extrapolating schemes).
PGMVNHXT0 corresponding to GMV of YT0%MNHX. It is not used.
PGMVNHXT0L and PGMVNHXT0M derivatives corresponding to GMV of YT0%MNHXL and YT0%MNHXM. These are used to compute the Z term and advection term in the Eulerian scheme.
PGMVNHXT1 corresponding to GMV of YT1%MNHX. It contains tendency of X 
ZNHXT0 allocated at CPG level. It contains X diagnosed at the time level t and is used to compute the tendency of X.
ZNHXT9 allocated at CPG level. It contains X diagnosed at the time level t-dt (3TL schemes only) and is used to compute the tendency of X.
PNHXL9: part of SLB1 buffer (MSLB1NHX9). It contains the X quantity to be interpolated (case of 2TL extrapolating schemes).
  
Difference with respect to the old solution using the so-called auxilary variable „VDAUX“

There is no special structure of GMV YAUX. The special call to spectral trasnforms does not exist any more. In ALADIN, the biperiodic extension (EBIPAUX) is not made any more at every time step. It is ensured by lateral coupling of X as for other variables. The large scale values of X are present in GMVCPL array which is initialised when reading the large scale coupling files (both d4 – SPSVD and X – SPNHX are diagnosed at SPNH_CONV_NHVAR level).

In the old solution the content of VDAUX was changing within the time-step. At its beginning it contained d3 variable computed as d4 (result of the SI solver) – X (first order estimate) at the level of ESPNHSI and subjected to the horizontal diffusion (ESPCHOR). Therefore the derivatives available were those of d3 and not X. But the VDAUX at t+dt contained (-X), which was biperiodicised by EBIPAUX and ready for the subtraction at ESPNHSI to provide d3. This change of the content of VDAUX within a time-step was very likely convenient for building the ND4SYS options, where approximated values of d3 were re-used. These ND4SYS subvariants are now suppressed. The extra variable contains always X term. It provides a more readable code. For example, the tendency of X is directly comparable to the tendency of d4 (as concerns its X part) which was not the case before and was making the code illogical and less clear.

Validations made so far

New code was validated with a help of Patrick Saez running the mitraillette tool for both ARPEGE and ALADIN. There are no crashes except when LGWADV=.TRUE. but this option does not work in the reference code (cy29t4 main) either. This option should be debugged later.
Beside, 2D potential flow case was run and results checked for some basic configurations:
2TL PC with d3 and d4; 3TL SI with d3 and d4, Euler PC with d3 and d4. All runs with d3 variable give bit identical results with respect to the reference. The d4 runs do not give bit identical results, which is to be expected, but the solution is OK with respect to the reference.

Concerning the code efficieny, we may assume that the goal is reached. In the adiabatic 2TL PC run the CPU cost gets to 60% of the reference code and in the Euler PC run it gets even to 40% of the original cost. We may say that the overhead of using d4 with respect to d3 became negligible. 

List of modifications in the code

The modifications are in the pack: mrpe684_CY29T4_dctyri
It is made public.

Added decks
None.

Modified decks

Arp/adiab/call_sl.F90
Remove PB2AUX structure of VDAUX. Add SLB1 pointer MSLB1NHX9. Accordingly modified call to LAPINEB.

Arp/adiab/cpeuldyn.F90
Remove ND4SYS option and VDAUX structure. Change of the d4 and X tendency computation: use diagnosed values only. Proper name convention.

Arp/adiab/cpg.F90
Remove all VDAUX-type structures. Allocating ZNHXT0 and ZNHXT9 to diagnose X properly. These two arrays should be dimensioned by NPROMA since they are used also in the hydrostatic model. Modified calls to CPG_GP and CPG_DYN.

Arp/adiab/cpg_dyn.F90
Modified calling sequence. Remove VDAUX structures. Use of diagnostic PNHX0 and PNHX9. Modified call to LACDYN (also adding new SLB1 pointer MSLB1NHX9) and CPEULDYN.

Arp/adiab/cpg_end.F90
Remove VDAUX structure and modifying call to GPMPFC.

Arp/adiab/cpg_gp.F90
Modified call: suppressing PB2AUX and introducing PNHXT0 and PNHXT9, which are diagnosed at this level, replacing local ZX0 and ZX9. Modified call to GPSPNG, GPMPFC, GPINISLB. Since we have now the derivatives of X available, the derivatives of d3, needed for the Z term computation, are set as derivatives of d4 minus derivatives of X. Proper name convention is used. The X term is now properly distinguished from the Y residual (LPCOLD option).

Arp/adiab/cpglag.F90
Remove PB2AUX from the call. Modified call to GPENDTR due to PB2AUX.

Arp/adiab/gnh_conv_nhvar.F90
Replacing PSVDAUX by PSNHX. Proper diagnostics of PSNHX.

Arp/adiab/gpendtr.F90
Complete suppression of VDAUX structures.

Arp/adiab/gpinislb.F90
Remove VDAUX structures. Initialisation of PGMVNHXT1 by zero or by PNHXT9 diagnosed (Eulerian case).

Arp/adiab/gpmpfc.F90
Treating derivatives of X at the time t instead of those of VDAUX.

Arp/adiab/gptf1.F90
Remove VDAUX as well as the time filter of X: it is not needed.

Arp/adiab/gptf1pc.F90
Proper naming convention for X not to mislead it with the Y residual (used only if LPCOLD)

Arp/adiab/gptf2.F90
Remove VDAUX as well as the  time filter on X: it is not needed.

Arp/adiab/gptf2pc.F90
Proper naming convention for X not to mislead it with the Y residual (used only if LPCOLD)

Arp/adiab/lacdyn.F90
Proper naming convention, suppressing the VDAUX strucure. Passing correct X arguments to LATTEX and LANHSIB.

Arp/adiab/lanhsib.F90
Using proper quantity X not to mislead it with the Y residual (used only if LPCOLD)

Arp/adiab/lapineb.F90
Replacing the VDAUX arrays by proper NHX arrays.

Arp/adiab/larcinb.F90
Proper name convention: VDAUX => NHX

Arp/adiab/lattex.F90
Remove VDAUX structure and ND4SYS options. Tendency of X is ensured by the diagnosed quantities only and is in agreement with the „X“ tendency part of d4. The GMV of YT9%MNHX is used to transport the X quantity between the predictor and corrector step when extrapolating 2TL schemes. 

Arp/adiab/spchor.F90
Horizontal diffusion on X. Proper name convention VDAUX => NHX.

Arp/adiab/spnh_conv_nhvar.F90
Proper name convention: VDAUX => NHX

Arp/adiab/spnhsi.F90
Proper name convention: VDAUX => NHX. Now there is no real computation done with X: the swap to d3 was suppressed. The structure is simply maintained in case some SI computations on X are introduced later although it looks improbable by now.

Arp/control/gp_model.F90
Removal of all VDAUX features, especially PB2 buffer of VDAUX (ZSLB2VDAUX). Because of this removed PB2 type of buffer for VDAUX, calls to CPG, CALL_SL and CPGLAG are modified.

Arp/control/scan2mdm.F90
Remove all VDAUX features and suppress call to INI2SCAN2M

Arp/control/sim4d.F90
Replace SPVDAUX by SPNH in declarations and call to SUSPEC.

Arp/control/spcm.F90
Replace SPVDAUX by SPNHX etc. in declarations and calls to subroutines.

Arp/control/spcmad.F90
Idem (spcm.F90 modset).

Arp/control/testlievol.F90
Replace SPVDAUX by SPNHX in declarations and call to SUSPEC.

Arp/dfi/edfi2.F90
Replace SPVDAUX by SPNHX in declarations and call to SUSPEC.

Arp/dia/wrmlppadm.F90
Replace VDAUX by NHX.

Arp/dia/wrspeca.F90
Fix the number of fields to be written to the file. The 3D field X is a part of NF3D when (NVDVAR==4 .AND. .NOT.LGWADV) but it is not written to the file. NFLEVG fields ought to be subtracted to get the right number.

Arp/module/gmv_subs.F90
Remove all YOMAUX declarations, remove YAUX and all AUX pointers, as well as SETUP_AUX subroutine. Introduce new pointers: MNHX, MNHXL, MNHXM in T0,  MNHX in T9 and T1. Replace pointers MD3YNL, MCD3YNL by MNHXNL, MCNHXNL in T9.

Arp/module/ptrslb1.F90
Introduce proper MSLB1NHX9 pointer to SLB1 buffer.

Arp/module/type_gmvs.F90
Remove AUX pointers, introduce MNHX familly of pointers. Remove TYPE_AUX.

Arp/module/yemct0.F90
Introduce NBICNHX coupling switch for X.

Arp/module/yemspbc.F90
Remove SPA3AUXBC array.

Arp/module/yomdyna.F90
Remove ND4SYS option.

Arp/module/yomgmv.F90
Remove TYPE_AUX.

Arp/module/yomsp.F90
Introduce SPNHX (SP5NHX, SP7NHX) among SPHV spectral arrays.

Arp/namelist/namdyna.h
Remove ND4SYS option.

Arp/parallel/trmtos.F90
Replace VDAUX by NHX.

Arp/parallel/trstom.F90
Replace VDAUX by NHX.

Arp/setup/su1yom.F90
Replace SPVDAUX by SPNHX in declarations and call to SUSPEC.

Arp/setup/sudim1.F90
Remove YOMAUX and LVDAUX.

Arp/setup/sudim2.F90
Remove all YOMAUX setup. Increase NFTHER by one (X) when NVDVAR==4 .AND. .NOT.LGWADV.

Arp/setup/sudyn.F90
Remove YOMAUX and LVDAUX.

Arp/setup/sudyna.F90
Remove ND4SYS option.

Arp/setup/suinif.F90
Replace SPVDAUX by SPNHX in declarations and call to SUSPEC.

Arp/setup/suslb.F90
Remove YOMAUX and all related setup. Introduce MSLB1NHX9 pointer and set it up when NVDVAR==4.

Arp/setup/suspec.F90
Replace PSPVDAUX by PSPNHX.

Arp/setup/suspeca.F90
Replace PSPVDAUX by PSPNHX. Fix the number of fields to be read from the file. The 3D field X is a part of NF3D when (NVDVAR==4 .AND. .NOT.LGWADV) but it is not present in the file. NFLEVG fields ought to be subtracted to get the right number. 

Arp/sinvect/cun3.F90
Replace SPVDAUX by SPNHX in declarations and call to SUSPEC.

Arp/transform/transdir_nhconv.F90
Replace SPVDAUX by SPNHX in declarations and calls to transforms.

Arp/utility/dealspa.F90
Remove SPAUX, SPA3AUX, SPVDAUX. Introduce SPNHX, SP5NHX, SP7NHX.

Arp/utility/sualspa.F90
Remove SPAUX, SPA3AUX, SPVDAUX. Introduce SPNHX, SP5NHX, SP7NHX.

Arp/var/avarct.F90
Replace ZCVDAUX by ZCNHX in declarations and call to SUSPEC.

Arp/var/rd801.F90
Replace ZSPVDAUX by ZSPNHX in declarations and call to SUSPEC.

Arp/var/suecges.F90
Replace ZSPVDAUX by ZSPNHX in declarations and call to SUSPEC.

Arp/var/suvazx.F90
Replace SPVDAUX by SPNHX in declarations and call to SUSPEC.

Ald/adiab/espchor.F90
Idem SPCHOR.

Ald/adiab/espcm.F90
Idem SPCM.

Ald/adiab/espcmad.F90
Idem SPCMAD.

Ald/adiab/espnhsi.F90
Idem SPNHSI.

Ald/adiab/gpspng.F90
VDAUX replaced by PGMVNHX*.

Ald/coupling/ecoupl1.F90
X at T1 interface for ALADIN coupling.

Ald/coupling/esrlxt1.F90
Insert the lateral coupling for X.

Ald/inidata/elsirf.F90
Replace SPVDAUX by SPNHX in declarations and call to SUSPEC.

Ald/inidata/erlbc.F90
Remove all YOMAUX business.

Ald/parallel/etrmtos.F90
Idem TRMTOS

Ald/parallel/etrstom.F90
Idem TRSTOM

Ald/setup/suebicu.F90
Setup of the alpha relaxation function for X.

Ald/setup/suect0.F90
Initialise NBICNHX to drive the lateral coupling of X. It is not in the namelist; it depends on NVDVAR==4 .AND. .NOT.LGWADV.

Ald/setup/suedim.F90
Setup dimension of the relaxation function array pending NBICNHX.

Ald/transform/etransdir_mdl.F90
Remove YOMAUX business and external EBIPAUX, remove the extra AUX call to transforms.

Ald/transform/etransdir_nhconv.F90
Replace SPVDAUX by SPNHX in declarations and calls to transforms.

Ald/transform/etransinv_mdl.F90
Remove PSVDAUX, YOMAUX, YAUX business.

Ald/transform/etransinvh.F90
Remove YOMAUX bussines and SPVDAUX argument in the call to ETRANSINV_MDL.

Ald/var/evarjkini.F90
Replace ZCVDAUX by ZCNHX in declarations and call to SUSPEC.

Decks to be deleted

Arp/control/ini2scan2m.F90

Arp/module/yomaux.F90

Ald/coupling/ebipaux.F90