Passive tracer vertical diffusivity in trazdf.F90 (PISCES offline)

[mgali]

Context

When running PISCES offline simulations in 1D, I noticed that passive tracers were not being diffused vertically, leading to rapid depletion of surface ocean nutrients and plankton stocks. The problem was less evident in PISCES offline 3D simulations because it was masked by nutrient advection. Still, the lack of vertical diffusion caused a sizable effect (I can provide plots).

Analysis

Passive tracers should be vertically diffused with the same diffusivity coefficient as salinity. This is correctly specified in the namelist_cfg, where the vertical diffusivity profiles of salinity (variable 'voddmavs' in dyna_grid_W.nc input file) are assigned to the array 'sn_avt'. Then, the input field 'voddmavs' is read into the array 'sn_avt' in the routine dtadyn.F90. However, in the routine trazdf.F90, line 154, the the salinity/tracer diffusivity stored in the array sn_avt are not used in the offline case, because cdtype == 'TRC', not 'TRA':

! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers
IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN   ;   zwt(:,:,2:jpk) = avt(:,:,2:jpk)
ELSE                                            ;   zwt(:,:,2:jpk) = avs(:,:,2:jpk)
ENDIF

Fix

The bug can be fixed by replacing "avs" by "avt" in line 154 of trazdf.F90:

! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers
IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN   ;   zwt(:,:,2:jpk) = avt(:,:,2:jpk)
ELSE                                            ;   zwt(:,:,2:jpk) = avt(:,:,2:jpk)
ENDIF

• Additional notes:

1. the temperature and salinity vertical diffusivity fields have exactly the same values in the input files used for PISCES offline (dyna_grid_W.nc). The current configuration of the code is confusing because the salinity diffusivity ('voddmavs') is read into the variable avt, which a priori should carry the temperature diffusivity (with salinity/tracer diffusivity being carreid by avs).

2. It seems that another workaround could be possible, by enabling double diffusive mixing in the namelist (ln_zdfddm = .true.). I did not test this option. However, enabling double diffusive mixing seems to introduce unnecessary complexity, as it would require (i) using the zdfddm.F90 routine, (ii) reading an additional input variable (temperature -avt- and salinity -avs- diffusivities), and (iii) might interfere in ways I do not understand with the maximum diffusivity for salinity in the namelist (by default, rn_avts = 1.e-4 ; note this is much lower than the enhanced tracer diffusivity of 100 m2/s used in during convection events).

[mgali]

Temporal evolution over 10 years of nitrate concentration at surface and 30 m depth in the Labrador Sea. LEFT: with bug. RIGHT: bug fixed.

[mgali]

Mean annual surface nitrate concentration in year 3 of PISCES offline simulation. LEFT: with bug. RIGHT: bug fixed.

[mgali]

Mean annualconcentration of diatom Chl at surface in year 3 of PISCES offline simulation. LEFT: with bug. RIGHT: bug fixed.

[cetlod]

Hummmm not sure it is bug because in dtadyn.F90 the avs array is always defined and is equal to avt