From efbb9f8dbf387f84cef8b8ebaa4214a2cb4394db Mon Sep 17 00:00:00 2001
From: Clement Rousset <clement.rousset@locean.ipsl.fr>
Date: Tue, 4 Apr 2023 16:19:41 +0000
Subject: [PATCH] Resolve "debugging the main"
---
cfgs/SHARED/field_def_nemo-oce.xml | 2 +-
src/ICE/icedyn_adv_pra.F90 | 10 +++++-----
src/ICE/icedyn_adv_umx.F90 | 10 +++++-----
src/OCE/DIA/diaar5.F90 | 15 ++++++++++-----
src/OCE/DIA/diacfl.F90 | 16 ++++++++--------
src/OCE/SBC/sbcrnf.F90 | 23 +++++++++++++----------
src/OCE/TRA/traldf_iso.F90 | 5 +++++
src/OCE/ZDF/zdfiwm.F90 | 4 ++--
src/OCE/ZDF/zdfphy.F90 | 2 +-
9 files changed, 50 insertions(+), 37 deletions(-)
diff --git a/cfgs/SHARED/field_def_nemo-oce.xml b/cfgs/SHARED/field_def_nemo-oce.xml
index d66e26cc..293469fd 100644
--- a/cfgs/SHARED/field_def_nemo-oce.xml
+++ b/cfgs/SHARED/field_def_nemo-oce.xml
@@ -772,7 +772,7 @@ that are available in the tidal-forcing implementation (see
<!-- dissipation diagnostics (note: ediss_k is only available with tke scheme) -->
<field id="avt_k" long_name="vertical eddy diffusivity from closure schemes" standard_name="ocean_vertical_eddy_diffusivity" unit="m2/s" grid_ref="grid_W_3D_inner" />
- <field id="avm_k" long_name="vertical eddy viscosity from closure schemes" standard_name="ocean_vertical_eddy_viscosity" unit="m2/s" />
+ <field id="avm_k" long_name="vertical eddy viscosity from closure schemes" standard_name="ocean_vertical_eddy_viscosity" unit="m2/s" grid_ref="grid_W_3D_inner" />
<field id="ediss_k" long_name="Kolmogorov energy dissipation (tke scheme)" standard_name="Kolmogorov_energy_dissipation" unit="W/kg" grid_ref="grid_W_3D_inner" />
<field id="eshear_k" long_name="energy source from vertical shear" standard_name="energy_source_from_shear" unit="W/kg" grid_ref="grid_W_3D_inner" />
<field id="estrat_k" long_name="energy sink from stratification" standard_name="energy_sink_from_stratification" unit="W/kg" grid_ref="grid_W_3D_inner" />
diff --git a/src/ICE/icedyn_adv_pra.F90 b/src/ICE/icedyn_adv_pra.F90
index 5a3ff037..4302e40a 100644
--- a/src/ICE/icedyn_adv_pra.F90
+++ b/src/ICE/icedyn_adv_pra.F90
@@ -1028,7 +1028,7 @@ CONTAINS
z1_dt = 1._wp / pdt
!
DO_2D( ihls, ihls, ihls, ihls )
- IF ( pv_i(ji,jj,jcat) > 0._wp ) THEN
+ IF ( pv_i(ji,jj,jcat) > 0._wp .AND. pa_i(ji,jj,jcat) > 0._wp ) THEN
!
! ! -- check h_ip -- !
! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
@@ -1065,7 +1065,7 @@ CONTAINS
! ! -- check s_i -- !
IF( nn_icesal == 4 ) THEN
DO_3D( ihls, ihls, ihls, ihls, 1, nlay_i )
- IF ( pv_i(ji,jj,jcat) > 0._wp ) THEN
+ IF ( pv_i(ji,jj,jcat) > 0._wp .AND. pa_i(ji,jj,jcat) > 0._wp ) THEN
! if szv_i/v_i is larger than the surrounding 9 pts => put the salt excess in the ocean
zsi = pszv_i(ji,jj,jk,jcat) / pv_i(ji,jj,jcat)
IF( zsi > pszi_max(ji,jj,jk) .AND. pa_i(ji,jj,jcat) < 0.15 ) THEN
@@ -1077,7 +1077,7 @@ CONTAINS
END_3D
ELSE
DO_2D( ihls, ihls, ihls, ihls )
- IF ( pv_i(ji,jj,jcat) > 0._wp ) THEN
+ IF ( pv_i(ji,jj,jcat) > 0._wp .AND. pa_i(ji,jj,jcat) > 0._wp ) THEN
! if s_i is larger than the surrounding 9 pts => put salt excess in the ocean
zsi = psv_i(ji,jj,jcat) / pv_i(ji,jj,jcat)
IF( zsi > psi_max(ji,jj) .AND. pa_i(ji,jj,jcat) < 0.15 ) THEN
@@ -1091,7 +1091,7 @@ CONTAINS
ENDIF
! ! -- check e_i/v_i -- !
DO_3D( ihls, ihls, ihls, ihls, 1, nlay_i )
- IF ( pv_i(ji,jj,jcat) > 0._wp ) THEN
+ IF ( pv_i(ji,jj,jcat) > 0._wp .AND. pa_i(ji,jj,jcat) > 0._wp ) THEN
! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean
zei = pe_i(ji,jj,jk,jcat) / pv_i(ji,jj,jcat)
IF( zei > pei_max(ji,jj,jk) .AND. pa_i(ji,jj,jcat) < 0.15 ) THEN
@@ -1103,7 +1103,7 @@ CONTAINS
END_3D
! ! -- check e_s/v_s -- !
DO_3D( ihls, ihls, ihls, ihls, 1, nlay_s )
- IF ( pv_s(ji,jj,jcat) > 0._wp ) THEN
+ IF ( pv_s(ji,jj,jcat) > 0._wp .AND. pa_i(ji,jj,jcat) > 0._wp ) THEN
! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean
zes = pe_s(ji,jj,jk,jcat) / pv_s(ji,jj,jcat)
IF( zes > pes_max(ji,jj,jk) .AND. pa_i(ji,jj,jcat) < 0.15 ) THEN
diff --git a/src/ICE/icedyn_adv_umx.F90 b/src/ICE/icedyn_adv_umx.F90
index 4fe7d9a0..24ef858b 100644
--- a/src/ICE/icedyn_adv_umx.F90
+++ b/src/ICE/icedyn_adv_umx.F90
@@ -1739,7 +1739,7 @@ CONTAINS
!
DO jl = 1, jpl
DO_2D( 0, 0, 0, 0 )
- IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
+ IF ( pv_i(ji,jj,jl) > 0._wp .AND. pa_i(ji,jj,jl) > 0._wp ) THEN
!
! ! -- check h_ip -- !
! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
@@ -1778,7 +1778,7 @@ CONTAINS
IF( nn_icesal == 4 ) THEN
DO jl = 1, jpl
DO_3D( 0, 0, 0, 0, 1, nlay_i )
- IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
+ IF ( pv_i(ji,jj,jl) > 0._wp .AND. pa_i(ji,jj,jl) > 0._wp ) THEN
! if szv_i/v_i is larger than the surrounding 9 pts => put the salt excess in the ocean
zsi = pszv_i(ji,jj,jk,jl) / pv_i(ji,jj,jl)
IF( zsi > pszi_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
@@ -1792,7 +1792,7 @@ CONTAINS
ELSE
DO jl = 1, jpl
DO_2D( 0, 0, 0, 0 )
- IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
+ IF ( pv_i(ji,jj,jl) > 0._wp .AND. pa_i(ji,jj,jl) > 0._wp ) THEN
! if s_i is larger than the surrounding 9 pts => put salt excess in the ocean
zsi = psv_i(ji,jj,jl) / pv_i(ji,jj,jl)
IF( zsi > psi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
@@ -1809,7 +1809,7 @@ CONTAINS
! ! -- check e_i/v_i -- !
DO jl = 1, jpl
DO_3D( 0, 0, 0, 0, 1, nlay_i )
- IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
+ IF ( pv_i(ji,jj,jl) > 0._wp .AND. pa_i(ji,jj,jl) > 0._wp ) THEN
! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean
zei = pe_i(ji,jj,jk,jl) / pv_i(ji,jj,jl)
IF( zei > pei_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
@@ -1823,7 +1823,7 @@ CONTAINS
! ! -- check e_s/v_s -- !
DO jl = 1, jpl
DO_3D( 0, 0, 0, 0, 1, nlay_s )
- IF ( pv_s(ji,jj,jl) > 0._wp ) THEN
+ IF ( pv_s(ji,jj,jl) > 0._wp .AND. pa_i(ji,jj,jl) > 0._wp ) THEN
! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean
zes = pe_s(ji,jj,jk,jl) / pv_s(ji,jj,jl)
IF( zes > pes_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
diff --git a/src/OCE/DIA/diaar5.F90 b/src/OCE/DIA/diaar5.F90
index d3002417..4e9a4f45 100644
--- a/src/OCE/DIA/diaar5.F90
+++ b/src/OCE/DIA/diaar5.F90
@@ -166,11 +166,11 @@ CONTAINS
!
ENDIF
- IF( iom_use( 'botpres' ) .OR. iom_use( 'sshthster' ) .OR. iom_use( 'sshsteric' ) ) THEN
+ IF( iom_use( 'sshthster' ) ) THEN
ALLOCATE( zrhd(T2D(0),jpk) )
!
DO_3D( 0, 0, 0, 0, 1, jpk )
- ztsn(ji,jj,jk,jp_tem,Kmm) = ts(ji,jj,jk,jp_tem,Kmm) ! thermosteric ssh
+ ztsn(ji,jj,jk,jp_tem,Kmm) = ts(ji,jj,jk,jp_tem,Kmm)
ztsn(ji,jj,jk,jp_sal,Kmm) = sn0(ji,jj,jk)
END_3D
CALL eos( ztsn, Kmm, zrhd, kbnd=0 ) ! now in situ density using initial salinity
@@ -202,6 +202,10 @@ CONTAINS
zztmp = - zarho / area_tot
CALL iom_put( 'sshthster', zztmp )
ENDIF
+
+ ENDIF
+
+ IF( iom_use( 'botpres' ) .OR. iom_use( 'sshsteric' ) .OR. iom_use( 'masstot' ) ) THEN
! ! steric sea surface height
z2d(:,:) = 0._wp ! no atmospheric surface pressure, levitating sea-ice
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
@@ -225,6 +229,9 @@ CONTAINS
zarho = glob_sum( 'diaar5', zsarho2 )
zztmp = - zarho / area_tot
CALL iom_put( 'sshsteric', zztmp )
+ !
+ zztmp = rho0 * ( zarho + zvol )
+ CALL iom_put( 'masstot', zztmp )
END IF
! ! ocean bottom pressure
zztmp = rho0 * grav * 1.e-4_wp ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa
@@ -235,7 +242,7 @@ CONTAINS
!
ENDIF
- IF( iom_use( 'masstot' ) .OR. iom_use( 'temptot' ) .OR. iom_use( 'saltot' ) ) THEN
+ IF( iom_use( 'temptot' ) .OR. iom_use( 'saltot' ) ) THEN
! ! Mean density anomalie, temperature and salinity
ztsn(:,:,:,:,Kmm) = 0._wp ! ztsn(:,:,1,jp_tem/sal) is used here as 2D Workspace for temperature & salinity
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
@@ -264,9 +271,7 @@ CONTAINS
IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only for the last tile
ztemp = glob_sum( 'diaar5', zstemp )
zsal = glob_sum( 'diaar5', zssal )
- zztmp = rho0 * ( zarho + zvol )
!
- CALL iom_put( 'masstot', zztmp )
CALL iom_put( 'temptot', ztemp / zvol )
CALL iom_put( 'saltot' , zsal / zvol )
ENDIF
diff --git a/src/OCE/DIA/diacfl.F90 b/src/OCE/DIA/diacfl.F90
index 6ce29566..202f0622 100644
--- a/src/OCE/DIA/diacfl.F90
+++ b/src/OCE/DIA/diacfl.F90
@@ -93,9 +93,9 @@ CONTAINS
CALL mpp_maxloc( 'diacfl', zCw_cfl, llmsk, zCw_max, iloc_w )
!
IF( lwp ) THEN ! write out to file
- WRITE(numcfl,FMT='(2x,i6,3x,a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zCu_max, iloc_u(1), iloc_u(2), iloc_u(3)
- WRITE(numcfl,FMT='(11x, a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') 'Max Cv', zCv_max, iloc_v(1), iloc_v(2), iloc_v(3)
- WRITE(numcfl,FMT='(11x, a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') 'Max Cw', zCw_max, iloc_w(1), iloc_w(2), iloc_w(3)
+ WRITE(numcfl,FMT='(2x,i6,3x,a6,4x,f8.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zCu_max, iloc_u(1), iloc_u(2), iloc_u(3)
+ WRITE(numcfl,FMT='(11x, a6,4x,f8.4,1x,i4,1x,i4,1x,i4)') 'Max Cv', zCv_max, iloc_v(1), iloc_v(2), iloc_v(3)
+ WRITE(numcfl,FMT='(11x, a6,4x,f8.4,1x,i4,1x,i4,1x,i4)') 'Max Cw', zCw_max, iloc_w(1), iloc_w(2), iloc_w(3)
ENDIF
!
! ! update maximum Courant numbers from whole run if applicable
@@ -107,13 +107,13 @@ CONTAINS
IF( kt == nitend .AND. lwp ) THEN
! to ascii file
WRITE(numcfl,*) '******************************************'
- WRITE(numcfl,FMT='(3x,a12,6x,f7.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cu', rCu_max, nCu_loc(1), nCu_loc(2), nCu_loc(3)
+ WRITE(numcfl,FMT='(3x,a12,6x,f8.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cu', rCu_max, nCu_loc(1), nCu_loc(2), nCu_loc(3)
WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', rDt/rCu_max
WRITE(numcfl,*) '******************************************'
- WRITE(numcfl,FMT='(3x,a12,6x,f7.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cv', rCv_max, nCv_loc(1), nCv_loc(2), nCv_loc(3)
+ WRITE(numcfl,FMT='(3x,a12,6x,f8.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cv', rCv_max, nCv_loc(1), nCv_loc(2), nCv_loc(3)
WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', rDt/rCv_max
WRITE(numcfl,*) '******************************************'
- WRITE(numcfl,FMT='(3x,a12,6x,f7.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cw', rCw_max, nCw_loc(1), nCw_loc(2), nCw_loc(3)
+ WRITE(numcfl,FMT='(3x,a12,6x,f8.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cw', rCw_max, nCw_loc(1), nCw_loc(2), nCw_loc(3)
WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', rDt/rCw_max
CLOSE( numcfl )
!
@@ -146,8 +146,8 @@ CONTAINS
!
! create output ascii file
CALL ctl_opn( numcfl, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp, 1 )
- WRITE(numcfl,*) 'Timestep Direction Max C i j k'
- WRITE(numcfl,*) '******************************************'
+ WRITE(numcfl,*) 'Timestep Direction Max C i j k'
+ WRITE(numcfl,*) '*******************************************'
ENDIF
!
rCu_max = 0._wp
diff --git a/src/OCE/SBC/sbcrnf.F90 b/src/OCE/SBC/sbcrnf.F90
index 4ceb8bca..1b268e79 100644
--- a/src/OCE/SBC/sbcrnf.F90
+++ b/src/OCE/SBC/sbcrnf.F90
@@ -273,16 +273,6 @@ CONTAINS
! !== allocate runoff arrays
IF( sbc_rnf_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_rnf_alloc : unable to allocate arrays' )
!
- IF( .NOT. ln_rnf ) THEN ! no specific treatment in vicinity of river mouths
- ln_rnf_mouth = .FALSE. ! default definition needed for example by sbc_ssr or by tra_adv_muscl
- nkrnf = 0
- rnf (:,:) = 0.0_wp
- rnf_b (:,:) = 0.0_wp
- rnfmsk (:,:) = 0.0_wp
- rnfmsk_z(:) = 0.0_wp
- RETURN
- ENDIF
- !
! ! ============
! ! Namelist
! ! ============
@@ -294,6 +284,19 @@ CONTAINS
902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namsbc_rnf in configuration namelist' )
IF(lwm) WRITE ( numond, namsbc_rnf )
!
+ IF( .NOT. ln_rnf ) THEN ! no specific treatment in vicinity of river mouths
+ ln_rnf_mouth = .FALSE. ! default definition needed for example by sbc_ssr or by tra_adv_muscl
+ ln_rnf_tem = .FALSE.
+ ln_rnf_sal = .FALSE.
+ ln_rnf_icb = .FALSE.
+ nkrnf = 0
+ rnf (:,:) = 0.0_wp
+ rnf_b (:,:) = 0.0_wp
+ rnfmsk (:,:) = 0.0_wp
+ rnfmsk_z(:) = 0.0_wp
+ RETURN
+ ENDIF
+ !
! ! Control print
IF(lwp) THEN
WRITE(numout,*)
diff --git a/src/OCE/TRA/traldf_iso.F90 b/src/OCE/TRA/traldf_iso.F90
index 6ed2242e..20e8f67d 100644
--- a/src/OCE/TRA/traldf_iso.F90
+++ b/src/OCE/TRA/traldf_iso.F90
@@ -269,6 +269,11 @@ CONTAINS
# undef INN
# undef pt_in
!
+ END DO != end ij slab =!
+ !
+ !
+ DO jk = 1, jpkm1 != ij slab =!
+ !
! !* 2nd pass : bilaplacian = laplacian of 1st pass (zlap)
! computed over the INNER domain
# define iso_blp_p2
diff --git a/src/OCE/ZDF/zdfiwm.F90 b/src/OCE/ZDF/zdfiwm.F90
index 240adb95..5061a401 100644
--- a/src/OCE/ZDF/zdfiwm.F90
+++ b/src/OCE/ZDF/zdfiwm.F90
@@ -233,9 +233,9 @@ CONTAINS
IF( ln_mevar ) THEN ! Variable mixing efficiency case : modify zav_wave in the
DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224) regimes
IF( zReb(ji,jj,jk) > 480.00_wp ) THEN
- zav_wave(ji,jj,jk) = 3.6515_wp * rnu * SQRT( zReb(ji,jj,jk) )
+ zav_wave(ji,jj,jk) = 3.6515_wp * rnu * SQRT( MAX( 0._wp, zReb(ji,jj,jk) ) )
ELSEIF( zReb(ji,jj,jk) < 10.224_wp ) THEN
- zav_wave(ji,jj,jk) = 0.052125_wp * rnu * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
+ zav_wave(ji,jj,jk) = 0.052125_wp * rnu * zReb(ji,jj,jk) * SQRT( MAX( 0._wp, zReb(ji,jj,jk) ) )
ENDIF
END_3D
ENDIF
diff --git a/src/OCE/ZDF/zdfphy.F90 b/src/OCE/ZDF/zdfphy.F90
index 39d2bfa6..c5e0bb9b 100644
--- a/src/OCE/ZDF/zdfphy.F90
+++ b/src/OCE/ZDF/zdfphy.F90
@@ -382,7 +382,7 @@ CONTAINS
ENDIF
IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile
IF( iom_use('avt_k' ) ) CALL iom_put( 'avt_k' , avt_k * wmask(A2D(0),:) )
- IF( iom_use('avm_k' ) ) CALL iom_put( 'avm_k' , avm_k * wmask(:,: ,:) )
+ IF( iom_use('avm_k' ) ) CALL iom_put( 'avm_k' , avm_k(A2D(0),:) * wmask(A2D(0),:) )
IF( iom_use('estrat_k') ) CALL iom_put( 'estrat_k', - avt_k * rn2(A2D(0),:) * wmask(A2D(0),:) )
IF( iom_use('eshear_k') ) CALL iom_put( 'eshear_k', sh2 )
DEALLOCATE( sh2 )
--
GitLab