In the current 4.2 branch my colleague and I found the following bug.IT seems to be recently introduced with Update Momentum Trends to Work with Split Explicit Timestepping (6a1e8701) · Commits · NEMO Workspace / Nemo · GitLab
On first glance it doesn’t seem to be included in the 5 branch.
In the bdydyn.F90 there are two times calls to ALLOCATE zbdytrdu and zbdytrdv.
If one runs a configuration with l_trddyn and ll_dyn3d nemo crashes due to trying to allocate the fields twice.
Furthermore the definitions for :
maxdzT ,utrd_bdy,utrd_bdy2d,utrd_bdy_e3u,utrd_bdy2d_hu,vtrd_bdy,vtrd_bdy2d,vtrd_bdy_e3v,vtrd_bdy2d_hu
are forgotten to be provided in the cfgs/SHARED/field_def_nemo-oce.xml which are required for the xios calls in trddyn.F90
For our setup we now use the following modified version of bdydyn.f90 and field_def_nemo-oce.xml
bdydyn.f90
MODULE bdydyn
!!======================================================================
!! *** MODULE bdydyn ***
!! Unstructured Open Boundary Cond. : Apply boundary conditions to velocities
!!======================================================================
!! History : 1.0 ! 2005-02 (J. Chanut, A. Sellar) Original code
!! - ! 2007-07 (D. Storkey) Move Flather implementation to separate routine.
!! 3.0 ! 2008-04 (NEMO team) add in the reference version
!! 3.2 ! 2008-04 (R. Benshila) consider velocity instead of transport
!! 3.3 ! 2010-09 (E.O'Dea) modifications for Shelf configurations
!! 3.3 ! 2010-09 (D.Storkey) add ice boundary conditions
!! 3.4 ! 2011 (D. Storkey) rewrite in preparation for OBC-BDY merge
!!----------------------------------------------------------------------
!! bdy_dyn : split velocities into barotropic and baroclinic parts
!! and call bdy_dyn2d and bdy_dyn3d to apply boundary
!! conditions
!!----------------------------------------------------------------------
USE oce ! ocean dynamics and tracers
USE dom_oce ! ocean space and time domain
USE bdy_oce ! ocean open boundary conditions
USE bdydyn2d ! open boundary conditions for barotropic solution
USE bdydyn3d ! open boundary conditions for baroclinic velocities
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
USE in_out_manager !
USE domvvl ! variable volume
USE trd_oce ! trends: ocean variables
USE trddyn ! trend manager: dynamics
IMPLICIT NONE
PRIVATE
PUBLIC bdy_dyn ! routine called in dyn_nxt
!! * Substitutions
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: bdydyn.F90 13237 2020-07-03 09:12:53Z smasson $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE bdy_dyn( kt, Kbb, puu, pvv, Kaa, dyn3d_only )
!!----------------------------------------------------------------------
!! *** SUBROUTINE bdy_dyn ***
!!
!! ** Purpose : - Wrapper routine for bdy_dyn2d and bdy_dyn3d.
!!
!!----------------------------------------------------------------------
INTEGER , INTENT(in) :: kt ! Main time step counter
INTEGER , INTENT(in) :: Kbb, Kaa ! Ocean time level indices
REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! Ocean velocities (to be updated at open boundaries)
LOGICAL, OPTIONAL , INTENT(in) :: dyn3d_only ! T => only update baroclinic velocities
!
INTEGER :: jk, ii, ij, ib_bdy, ib, igrd ! Loop counter
LOGICAL :: ll_dyn2d, ll_dyn3d, ll_orlanski
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ua_tmp, va_tmp ! RDP: temporary save of 2d momentum
REAL(wp), DIMENSION(jpi,jpj) :: zua2d, zva2d ! after barotropic velocities
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zbdytrdu, zbdytrdv ! bdy trends
!!----------------------------------------------------------------------
!
ll_dyn2d = .true.
ll_dyn3d = .true.
!
IF( PRESENT(dyn3d_only) ) THEN
IF( dyn3d_only ) ll_dyn2d = .false.
ENDIF
!
IF ( l_trddyn ) THEN
IF ( ll_dyn3d ) THEN
! ALLOCATE( zbdytrdu(jpi,jpj,jpk), zbdytrdv(jpi,jpj,jpk) )
IF (.NOT. ALLOCATED ( zbdytrdu) ) THEN
ALLOCATE( zbdytrdu(jpi,jpj,jpk) )
zbdytrdu(:,:,:) = 0._wp
ENDIF
IF (.NOT. ALLOCATED ( zbdytrdv) ) THEN
ALLOCATE( zbdytrdv(jpi,jpj,jpk) )
zbdytrdv(:,:,:) = 0._wp
ENDIF
ENDIF
ENDIF
ll_orlanski = .false.
DO ib_bdy = 1, nb_bdy
IF ( cn_dyn2d(ib_bdy) == 'orlanski' .OR. cn_dyn2d(ib_bdy) == 'orlanski_npo' &
& .OR. cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo') ll_orlanski = .true.
END DO
IF ( l_trddyn ) THEN
IF (.NOT. ALLOCATED ( zbdytrdu) ) THEN
ALLOCATE( zbdytrdu(jpi,jpj,jpk) )
zbdytrdu(:,:,:) = puu(:,:,:,Kaa)
ENDIF
IF (.NOT. ALLOCATED ( zbdytrdv) ) THEN
ALLOCATE( zbdytrdv(jpi,jpj,jpk) )
zbdytrdv(:,:,:) = pvv(:,:,:,Kaa)
ENDIF
ENDIF
!-------------------------------------------------------
! Split velocities into barotropic and baroclinic parts
!-------------------------------------------------------
! ! "After" velocities:
zua2d(:,:) = 0._wp
zva2d(:,:) = 0._wp
DO jk = 1, jpkm1
zua2d(:,:) = zua2d(:,:) + e3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
zva2d(:,:) = zva2d(:,:) + e3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
END DO
zua2d(:,:) = zua2d(:,:) * r1_hu(:,:,Kaa)
zva2d(:,:) = zva2d(:,:) * r1_hv(:,:,Kaa)
DO jk = 1 , jpkm1
puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zua2d(:,:) ) * umask(:,:,jk)
pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zva2d(:,:) ) * vmask(:,:,jk)
END DO
IF( ll_orlanski ) THEN ! "Before" velocities (Orlanski condition only)
DO jk = 1 , jpkm1
puu(:,:,jk,Kbb) = ( puu(:,:,jk,Kbb) - uu_b(:,:,Kbb) ) * umask(:,:,jk)
pvv(:,:,jk,Kbb) = ( pvv(:,:,jk,Kbb) - vv_b(:,:,Kbb) ) * vmask(:,:,jk)
END DO
ENDIF
!-------------------------------------------------------
! Apply boundary conditions to barotropic and baroclinic
! parts separately
!-------------------------------------------------------
IF( ll_dyn2d ) CALL bdy_dyn2d( kt, zua2d, zva2d, uu_b(:,:,Kbb), vv_b(:,:,Kbb), r1_hu(:,:,Kaa), r1_hv(:,:,Kaa), ssh(:,:,Kaa) )
IF( ll_dyn3d ) CALL bdy_dyn3d( kt, Kbb, puu, pvv, Kaa )
!-------------------------------------------------------
! Recombine velocities
!-------------------------------------------------------
!
DO jk = 1 , jpkm1
puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) + zua2d(:,:) ) * umask(:,:,jk)
pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) + zva2d(:,:) ) * vmask(:,:,jk)
END DO
!
IF ( ll_orlanski ) THEN
DO jk = 1 , jpkm1
puu(:,:,jk,Kbb) = ( puu(:,:,jk,Kbb) + uu_b(:,:,Kbb) ) * umask(:,:,jk)
pvv(:,:,jk,Kbb) = ( pvv(:,:,jk,Kbb) + vv_b(:,:,Kbb) ) * vmask(:,:,jk)
END DO
END IF
!
IF ( l_trddyn ) THEN
zbdytrdu(:,:,:) = ( puu(:,:,:,Kaa) - zbdytrdu(:,:,:) ) * r1_Dt
zbdytrdv(:,:,:) = ( pvv(:,:,:,Kaa) - zbdytrdv(:,:,:) ) * r1_Dt
CALL trd_dyn( zbdytrdu, zbdytrdv, jpdyn_bdy, kt )
DEALLOCATE( zbdytrdu, zbdytrdv )
ENDIF
END SUBROUTINE bdy_dyn
!!======================================================================
END MODULE bdydyn
additions to field_def_nemo-oce.xml
<!-- next variables available with key_diahth -->
<field id="mlddzt" long_name="Thermocline Depth (depth of max dT/dz)" standard_name="depth_at_maximum_upward_derivative_of_sea_water_potential_temperature" unit="m" />
<field id="maxdzT" long_name="Maximum of dT/dz)" standard_name="maximum_upward_derivative_of_sea_water_potential_temperature" unit="degC/m" />
<field id="mldr10_3" long_name="Mixed Layer Depth (dsigma = 0.03 wrt 10m)" standard_name="ocean_mixed_layer_thickness_defined_by_sigma_theta" unit="m" />
<field_group id="trendU" grid_ref="grid_U_3D">
<!-- variables available with ln_dyn_trd -->
<field id="hu" long_name="total time-varying depth at U-points" unit="m" grid_ref="grid_U_2D" />
<field id="utrd_hpg" long_name="i-trend: hydrostatic pressure gradient" unit="m/s^2" />
<field id="utrd_spg2d" long_name="i-trend: surface pressure gradient: true trend" unit="m/s^2" grid_ref="grid_U_2D" />
<field id="utrd_keg" long_name="i-trend: KE gradient or hor. adv." unit="m/s^2" />
<field id="utrd_rvo" long_name="i-trend: relative vorticity or metric term" unit="m/s^2" />
<field id="utrd_pvo" long_name="i-trend: planetary vorticity: 3D component" unit="m/s^2" />
<field id="utrd_pvo2d" long_name="i-trend: planetary vorticity: 2D component" unit="m/s^2" grid_ref="grid_U_2D" />
<field id="utrd_zad" long_name="i-trend: vertical advection" unit="m/s^2" />
<field id="utrd_udx" long_name="i-trend: U.dx[U]" unit="m/s^2" />
<field id="utrd_ldf" long_name="i-trend: lateral diffusion" unit="m/s^2" />
<field id="utrd_zdf" long_name="i-trend: vertical diffusion" unit="m/s^2" />
<field id="utrd_frc2d" long_name="i-trend: barotropic constant forcing " unit="m/s^2" grid_ref="grid_U_2D" />
<field id="utrd_tau" long_name="i-trend: wind stress in top layer " unit="m/s^2" grid_ref="grid_U_2D" />
<field id="utrd_tau2d" long_name="i-trend: wind stress: 2D component " unit="m/s^2" grid_ref="grid_U_2D" />
<field id="utrd_bfr" long_name="i-trend: bottom friction (3D) in bottom layer" unit="m/s^2" />
<field id="utrd_bfr2d" long_name="i-trend: bottom friction: 2D component" unit="m/s^2" grid_ref="grid_U_2D" />
<field id="utrd_tfr" long_name="i-trend: top friction (3D) in top layer" unit="m/s^2" />
<field id="utrd_tfr2d" long_name="i-trend: top friction: 2D component" unit="m/s^2" grid_ref="grid_U_2D" />
<field id="utrd_tot" long_name="i-trend: total momentum trend before atf" unit="m/s^2" />
<field id="utrd_tot2d" long_name="i-trend: total momentum trend before atf: 2D component" unit="m/s^2" grid_ref="grid_U_2D" />
<field id="utrd_atf" long_name="i-trend: asselin time filter trend" unit="m/s^2" />
<field id="utrd_atm2d" long_name="i-trend: atmospheric pressure: 2D component " unit="m/s^2" grid_ref="grid_U_2D" />
<field id="utrd_bdy" long_name="i-trend: lateral boundary forcing" unit="m/s^2" />
<field id="utrd_bdy2d" long_name="i-trend: lateral boundary forcing: 2D component " unit="m/s^2" grid_ref="grid_U_2D" />
<!-- thickness weighted versions -->
<field id="utrd_hpg_e3u" unit="m2/s^2" > utrd_hpg * e3u </field>
<field id="utrd_keg_e3u" unit="m2/s^2" > utrd_keg * e3u </field>
<field id="utrd_rvo_e3u" unit="m2/s^2" > utrd_rvo * e3u </field>
<field id="utrd_pvo_e3u" unit="m2/s^2" > utrd_pvo * e3u </field>
<field id="utrd_zad_e3u" unit="m2/s^2" > utrd_zad * e3u </field>
<field id="utrd_ldf_e3u" unit="m2/s^2" > utrd_ldf * e3u </field>
<field id="utrd_zdf_e3u" unit="m2/s^2" > utrd_zdf * e3u </field>
<field id="utrd_tot_e3u" unit="m2/s^2" > utrd_tot * e3u </field>
<field id="utrd_atf_e3u" unit="m2/s^2" > utrd_atf * e3u </field>
<field id="utrd_bfr_e3u" unit="m2/s^2" > utrd_bfr * e3u </field>
<field id="utrd_tfr_e3u" unit="m2/s^2" > utrd_tfr * e3u </field>
<field id="utrd_bdy_e3u" unit="m2/s^2" > utrd_bdy * e3u </field>
<field id="utrd_spg2d_hu" unit="m2/s^2" grid_ref="grid_U_2D" > utrd_spg2d * hu </field>
<field id="utrd_pvo2d_hu" unit="m2/s^2" grid_ref="grid_U_2D" > utrd_pvo2d * hu </field>
<field id="utrd_frc2d_hu" unit="m2/s^2" grid_ref="grid_U_2D" > utrd_frc2d * hu </field>
<field id="utrd_tau2d_hu" unit="m2/s^2" grid_ref="grid_U_2D" > utrd_tau2d * hu </field>
<field id="utrd_tfr2d_hu" unit="m2/s^2" grid_ref="grid_U_2D" > utrd_tfr2d * hu </field>
<field id="utrd_bfr2d_hu" unit="m2/s^2" grid_ref="grid_U_2D" > utrd_bfr2d * hu </field>
<field id="utrd_tot2d_hu" unit="m2/s^2" grid_ref="grid_U_2D" > utrd_tot2d * hu </field>
<field id="utrd_atm2d_hu" unit="m2/s^2" grid_ref="grid_U_2D" > utrd_atm2d * hu </field>
<field id="utrd_bdy2d_hu" unit="m2/s^2" grid_ref="grid_U_2D" > utrd_bdy2d * hu </field>
</field_group>
<field_group id="trendV" grid_ref="grid_V_3D">
<!-- variables available with ln_dyn_trd -->
<field id="hv" long_name="total time-varying depth at V-points" unit="m" grid_ref="grid_V_2D" />
<field id="vtrd_hpg" long_name="j-trend: hydrostatic pressure gradient" unit="m/s^2" />
<field id="vtrd_spg2d" long_name="j-trend: surface pressure gradient: true trend" unit="m/s^2" grid_ref="grid_V_2D" />
<field id="vtrd_keg" long_name="j-trend: KE gradient or hor. adv." unit="m/s^2" />
<field id="vtrd_rvo" long_name="j-trend: relative vorticity or metric term" unit="m/s^2" />
<field id="vtrd_pvo" long_name="j-trend: planetary vorticity: 3D component" unit="m/s^2" />
<field id="vtrd_pvo2d" long_name="j-trend: planetary vorticity: 2D component" unit="m/s^2" grid_ref="grid_V_2D" />
<field id="vtrd_zad" long_name="j-trend: vertical advection" unit="m/s^2" />
<field id="vtrd_vdy" long_name="j-trend: V.dx[V]" unit="m/s^2" />
<field id="vtrd_ldf" long_name="j-trend: lateral diffusion" unit="m/s^2" />
<field id="vtrd_zdf" long_name="j-trend: vertical diffusion" unit="m/s^2" />
<field id="vtrd_frc2d" long_name="j-trend: barotropic constant forcing " unit="m/s^2" grid_ref="grid_V_2D" />
<field id="vtrd_tau" long_name="j-trend: wind stress in top layer " unit="m/s^2" grid_ref="grid_V_2D" />
<field id="vtrd_tau2d" long_name="j-trend: wind stress: 2D component " unit="m/s^2" grid_ref="grid_V_2D" />
<field id="vtrd_bfr" long_name="j-trend: bottom friction (3D) in bottom layer" unit="m/s^2" />
<field id="vtrd_bfr2d" long_name="j-trend: bottom friction: 2D component" unit="m/s^2" grid_ref="grid_V_2D" />
<field id="vtrd_tfr" long_name="j-trend: top friction (3D) in top layer" unit="m/s^2" />
<field id="vtrd_tfr2d" long_name="j-trend: top friction: 2D component" unit="m/s^2" grid_ref="grid_V_2D" />
<field id="vtrd_tot" long_name="j-trend: total momentum trend before atf" unit="m/s^2" />
<field id="vtrd_tot2d" long_name="j-trend: total momentum trend before atf: 2D component" unit="m/s^2" grid_ref="grid_V_2D" />
<field id="vtrd_atf" long_name="j-trend: asselin time filter trend" unit="m/s^2" />
<field id="vtrd_atm2d" long_name="j-trend: atmospheric pressure: 2D component " unit="m/s^2" grid_ref="grid_V_2D" />
<field id="vtrd_bdy" long_name="j-trend: lateral boundary forcing" unit="m/s^2" />
<field id="vtrd_bdy2d" long_name="j-trend: lateral boundary forcing: 2D component " unit="m/s^2" grid_ref="grid_V_2D" />
<!-- thickness weighted versions -->
<field id="vtrd_hpg_e3v" unit="m2/s^2" > vtrd_hpg * e3v </field>
<field id="vtrd_keg_e3v" unit="m2/s^2" > vtrd_keg * e3v </field>
<field id="vtrd_rvo_e3v" unit="m2/s^2" > vtrd_rvo * e3v </field>
<field id="vtrd_pvo_e3v" unit="m2/s^2" > vtrd_pvo * e3v </field>
<field id="vtrd_zad_e3v" unit="m2/s^2" > vtrd_zad * e3v </field>
<field id="vtrd_ldf_e3v" unit="m2/s^2" > vtrd_ldf * e3v </field>
<field id="vtrd_zdf_e3v" unit="m2/s^2" > vtrd_zdf * e3v </field>
<field id="vtrd_tot_e3v" unit="m2/s^2" > vtrd_tot * e3v </field>
<field id="vtrd_atf_e3v" unit="m2/s^2" > vtrd_atf * e3v </field>
<field id="vtrd_bfr_e3v" unit="m2/s^2" > vtrd_bfr * e3v </field>
<field id="vtrd_tfr_e3v" unit="m2/s^2" > vtrd_tfr * e3v </field>
<field id="vtrd_bdy_e3v" unit="m2/s^2" > vtrd_bdy * e3v </field>
<field id="vtrd_spg2d_hv" unit="m2/s^2" grid_ref="grid_V_2D" > vtrd_spg2d * hv </field>
<field id="vtrd_pvo2d_hv" unit="m2/s^2" grid_ref="grid_V_2D" > vtrd_pvo2d * hv </field>
<field id="vtrd_frc2d_hv" unit="m2/s^2" grid_ref="grid_V_2D" > vtrd_frc2d * hv </field>
<field id="vtrd_tau2d_hv" unit="m2/s^2" grid_ref="grid_V_2D" > vtrd_tau2d * hv </field>
<field id="vtrd_tfr2d_hv" unit="m2/s^2" grid_ref="grid_V_2D" > vtrd_tfr2d * hv </field>
<field id="vtrd_bfr2d_hv" unit="m2/s^2" grid_ref="grid_V_2D" > vtrd_bfr2d * hv </field>
<field id="vtrd_tot2d_hv" unit="m2/s^2" grid_ref="grid_V_2D" > vtrd_tot2d * hv </field>
<field id="vtrd_atm2d_hv" unit="m2/s^2" grid_ref="grid_V_2D" > vtrd_atm2d * hv </field>
<field id="vtrd_bdy2d_hu" unit="m2/s^2" grid_ref="grid_V_2D" > vtrd_bdy2d * hv </field>
</field_group>