diff --git a/src/ABL/ablmod.F90 b/src/ABL/ablmod.F90
index 84a2a99d389cdd865aa3c9de3dc9d91c53ba637b..d14691e4e56216314d10e8143d85aca8b7cace4c 100644
--- a/src/ABL/ablmod.F90
+++ b/src/ABL/ablmod.F90
@@ -1241,21 +1241,21 @@ CONTAINS
DO_2D( 0, 0, 1, 0 )
zFY ( ji, jj ) = zdY ( ji, jj ) &
- & + smth_a* ( (zdX ( ji, jj+1 ) - zdX( ji-1, jj+1 )) &
+ & + smth_a* ( (zdX ( ji, jj+1 ) - zdX( ji-1, jj+1 )) & ! add () for NP repro
& - (zdX ( ji, jj ) - zdX( ji-1, jj )) )
END_2D
DO_2D( 1, 0, 0, 0 )
zFX( ji, jj ) = zdX( ji, jj ) &
- & + smth_a*( (zdY( ji+1, jj ) - zdY( ji+1, jj-1)) &
+ & + smth_a*( (zdY( ji+1, jj ) - zdY( ji+1, jj-1)) & ! add () for NP repro
& - (zdY( ji , jj ) - zdY( ji , jj-1)) )
END_2D
DO_2D( 0, 0, 0, 0 )
pvar2d( ji ,jj ) = pvar2d( ji ,jj ) &
& + msk(ji,jj) * smth_b * ( &
- & zFX( ji, jj ) - zFX( ji-1, jj ) &
- & +zFY( ji, jj ) - zFY( ji, jj-1 ) )
+ & ( zFX( ji, jj ) - zFX( ji-1, jj ) ) & ! add () for NP repro
+ & + ( zFY( ji, jj ) - zFY( ji, jj-1 ) ) )
END_2D
!---------------------------------------------------------------------------------------------------
diff --git a/src/ICE/icedyn.F90 b/src/ICE/icedyn.F90
index 687cd916ed23110bcd9930e961649d6ab11dbe6a..4fdd22b66c271455802e33073ca4d366515639f9 100644
--- a/src/ICE/icedyn.F90
+++ b/src/ICE/icedyn.F90
@@ -164,8 +164,8 @@ CONTAINS
ALLOCATE( zdivu_i(A2D(0)) )
DO_2D( 0, 0, 0, 0 )
- zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
- & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj)
+ zdivu_i(ji,jj) = ( ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) ) & ! add () for NP repro
+ & + ( e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) ) * r1_e1e2t(ji,jj)
END_2D
CALL iom_put( 'icediv' , zdivu_i )
DEALLOCATE( zdivu_i )
diff --git a/src/ICE/icedyn_adv.F90 b/src/ICE/icedyn_adv.F90
index a6d2e957735e1c0663619f52f663e237f61c3215..6afabac00bcf4adc66c8cbd98cfb577df3151843 100644
--- a/src/ICE/icedyn_adv.F90
+++ b/src/ICE/icedyn_adv.F90
@@ -25,7 +25,6 @@ MODULE icedyn_adv
USE lib_mpp ! MPP library
USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
USE timing ! Timing
- USE prtctl ! Print control
IMPLICIT NONE
PRIVATE
@@ -108,6 +107,8 @@ CONTAINS
! controls
IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icedyn_adv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ') ! prints
+ IF( sn_cfctl%l_prtctl ) &
+ & CALL ice_prt3D('icedyn_adv') ! prints
IF( ln_timing ) CALL timing_stop ('icedyn_adv') ! timing
!
END SUBROUTINE ice_dyn_adv
diff --git a/src/ICE/icedyn_adv_pra.F90 b/src/ICE/icedyn_adv_pra.F90
index 54a90afdc79fcc2a8de05e609900434e5927a3b1..763f71a6611c1f296982b755a5aeafb9b3de820e 100644
--- a/src/ICE/icedyn_adv_pra.F90
+++ b/src/ICE/icedyn_adv_pra.F90
@@ -318,8 +318,9 @@ CONTAINS
! derive open water from ice concentration
zati2(:,:) = SUM( pa_i(A2D(0),:), dim=3 )
DO_2D( 0, 0, 0, 0 )
- pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & !--- open water
- & - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
+ pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & !--- open water
+ & - ( ( zudy(ji,jj) - zudy(ji-1,jj) ) & ! add () for NP repro
+ & + ( zvdx(ji,jj) - zvdx(ji,jj-1) ) ) * r1_e1e2t(ji,jj) * zdt
END_2D
CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T', 1.0_wp )
!
@@ -390,158 +391,152 @@ CONTAINS
DO jl = 1, jcat ! loop on categories
!
! Limitation of moments.
- DO jj = Njs0 - jj0, Nje0 + jj0
-
- DO ji = Nis0 - 1, Nie0 + 1
-
- zpsm = psm (ji,jj,jl) ! optimization
- zps0 = ps0 (ji,jj,jl)
- zpsx = psx (ji,jj,jl)
- zpsxx = psxx(ji,jj,jl)
- zpsy = psy (ji,jj,jl)
- zpsyy = psyy(ji,jj,jl)
- zpsxy = psxy(ji,jj,jl)
-
- ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
- zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 )
- !
- zslpmax = MAX( 0._wp, zps0 )
- zs1max = 1.5 * zslpmax
- zs1new = MIN( zs1max, MAX( -zs1max, zpsx ) )
- zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, zpsxx ) )
- rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask
-
- zps0 = zslpmax
- zpsx = zs1new * rswitch
- zpsxx = zs2new * rswitch
- zpsy = zpsy * rswitch
- zpsyy = zpsyy * rswitch
- zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch
-
- ! Calculate fluxes and moments between boxes i<-->i+1
- ! ! Flux from i to i+1 WHEN u GT 0
- zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
- zalf = MAX( 0._wp, put(ji,jj) ) * pdt / zpsm
- zalfq = zalf * zalf
- zalf1 = 1.0 - zalf
- zalf1q = zalf1 * zalf1
- !
- zfm (ji,jj) = zalf * zpsm
- zf0 (ji,jj) = zalf * ( zps0 + zalf1 * ( zpsx + (zalf1 - zalf) * zpsxx ) )
- zfx (ji,jj) = zalfq * ( zpsx + 3.0 * zalf1 * zpsxx )
- zfxx(ji,jj) = zalf * zpsxx * zalfq
- zfy (ji,jj) = zalf * ( zpsy + zalf1 * zpsxy )
- zfxy(ji,jj) = zalfq * zpsxy
- zfyy(ji,jj) = zalf * zpsyy
-
- ! ! Readjust moments remaining in the box.
- zpsm = zpsm - zfm(ji,jj)
- zps0 = zps0 - zf0(ji,jj)
- zpsx = zalf1q * ( zpsx - 3.0 * zalf * zpsxx )
- zpsxx = zalf1 * zalf1q * zpsxx
- zpsy = zpsy - zfy (ji,jj)
- zpsyy = zpsyy - zfyy(ji,jj)
- zpsxy = zalf1q * zpsxy
- !
- psm (ji,jj,jl) = zpsm ! optimization
- ps0 (ji,jj,jl) = zps0
- psx (ji,jj,jl) = zpsx
- psxx(ji,jj,jl) = zpsxx
- psy (ji,jj,jl) = zpsy
- psyy(ji,jj,jl) = zpsyy
- psxy(ji,jj,jl) = zpsxy
- !
- END DO
-
- DO ji = Nis0 - 1, Nie0
- ! ! Flux from i+1 to i when u LT 0.
- zalf = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl)
- zalg (ji,jj) = zalf
- zalfq = zalf * zalf
- zalf1 = 1.0 - zalf
- zalg1 (ji,jj) = zalf1
- zalf1q = zalf1 * zalf1
- zalg1q(ji,jj) = zalf1q
- !
- zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj,jl)
- zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj,jl) &
- & - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) )
- zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) )
- zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji+1,jj,jl) * zalfq
- zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) )
- zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj,jl)
- zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj,jl)
- END DO
-
- DO ji = Nis0, Nie0
- !
- zpsm = psm (ji,jj,jl) ! optimization
- zps0 = ps0 (ji,jj,jl)
- zpsx = psx (ji,jj,jl)
- zpsxx = psxx(ji,jj,jl)
- zpsy = psy (ji,jj,jl)
- zpsyy = psyy(ji,jj,jl)
- zpsxy = psxy(ji,jj,jl)
- ! ! Readjust moments remaining in the box.
- zbt = zbet(ji-1,jj)
- zbt1 = 1.0 - zbet(ji-1,jj)
- !
- zpsm = zbt * zpsm + zbt1 * ( zpsm - zfm(ji-1,jj) )
- zps0 = zbt * zps0 + zbt1 * ( zps0 - zf0(ji-1,jj) )
- zpsx = zalg1q(ji-1,jj) * ( zpsx + 3.0 * zalg(ji-1,jj) * zpsxx )
- zpsxx = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * zpsxx
- zpsy = zbt * zpsy + zbt1 * ( zpsy - zfy (ji-1,jj) )
- zpsyy = zbt * zpsyy + zbt1 * ( zpsyy - zfyy(ji-1,jj) )
- zpsxy = zalg1q(ji-1,jj) * zpsxy
+ DO_2D( 1, 1, jj0, jj0 )
+ !
+ zpsm = psm (ji,jj,jl) ! optimization
+ zps0 = ps0 (ji,jj,jl)
+ zpsx = psx (ji,jj,jl)
+ zpsxx = psxx(ji,jj,jl)
+ zpsy = psy (ji,jj,jl)
+ zpsyy = psyy(ji,jj,jl)
+ zpsxy = psxy(ji,jj,jl)
+
+ ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
+ zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 )
+ !
+ zslpmax = MAX( 0._wp, zps0 )
+ zs1max = 1.5 * zslpmax
+ zs1new = MIN( zs1max, MAX( -zs1max, zpsx ) )
+ zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, zpsxx ) )
+ rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask
+
+ zps0 = zslpmax
+ zpsx = zs1new * rswitch
+ zpsxx = zs2new * rswitch
+ zpsy = zpsy * rswitch
+ zpsyy = zpsyy * rswitch
+ zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch
+
+ ! Calculate fluxes and moments between boxes i<-->i+1
+ ! ! Flux from i to i+1 WHEN u GT 0
+ zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
+ zalf = MAX( 0._wp, put(ji,jj) ) * pdt / zpsm
+ zalfq = zalf * zalf
+ zalf1 = 1.0 - zalf
+ zalf1q = zalf1 * zalf1
+ !
+ zfm (ji,jj) = zalf * zpsm
+ zf0 (ji,jj) = zalf * ( zps0 + zalf1 * ( zpsx + (zalf1 - zalf) * zpsxx ) )
+ zfx (ji,jj) = zalfq * ( zpsx + 3.0 * zalf1 * zpsxx )
+ zfxx(ji,jj) = zalf * zpsxx * zalfq
+ zfy (ji,jj) = zalf * ( zpsy + zalf1 * zpsxy )
+ zfxy(ji,jj) = zalfq * zpsxy
+ zfyy(ji,jj) = zalf * zpsyy
+
+ ! ! Readjust moments remaining in the box.
+ zpsm = zpsm - zfm(ji,jj)
+ zps0 = zps0 - zf0(ji,jj)
+ zpsx = zalf1q * ( zpsx - 3.0 * zalf * zpsxx )
+ zpsxx = zalf1 * zalf1q * zpsxx
+ zpsy = zpsy - zfy (ji,jj)
+ zpsyy = zpsyy - zfyy(ji,jj)
+ zpsxy = zalf1q * zpsxy
+ !
+ psm (ji,jj,jl) = zpsm ! optimization
+ ps0 (ji,jj,jl) = zps0
+ psx (ji,jj,jl) = zpsx
+ psxx(ji,jj,jl) = zpsxx
+ psy (ji,jj,jl) = zpsy
+ psyy(ji,jj,jl) = zpsyy
+ psxy(ji,jj,jl) = zpsxy
+ !
+ END_2D
- ! Put the temporary moments into appropriate neighboring boxes.
- ! ! Flux from i to i+1 IF u GT 0.
- zbt = zbet(ji-1,jj)
- zbt1 = 1.0 - zbet(ji-1,jj)
- zpsm = zbt * ( zpsm + zfm(ji-1,jj) ) + zbt1 * zpsm
- zalf = zbt * zfm(ji-1,jj) / zpsm
- zalf1 = 1.0 - zalf
- ztemp = zalf * zps0 - zalf1 * zf0(ji-1,jj)
- !
- zps0 = zbt * ( zps0 + zf0(ji-1,jj) ) + zbt1 * zps0
- zpsx = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * zpsx + 3.0 * ztemp ) + zbt1 * zpsx
- zpsxx = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * zpsxx &
- & + 5.0 * ( zalf * zalf1 * ( zpsx - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) &
- & + zbt1 * zpsxx
- zpsxy = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * zpsxy &
- & + 3.0 * (- zalf1*zfy(ji-1,jj) + zalf * zpsy ) ) &
- & + zbt1 * zpsxy
- zpsy = zbt * ( zpsy + zfy (ji-1,jj) ) + zbt1 * zpsy
- zpsyy = zbt * ( zpsyy + zfyy(ji-1,jj) ) + zbt1 * zpsyy
+ DO_2D( 1, 0, jj0, jj0 )
+ ! ! Flux from i+1 to i when u LT 0.
+ zalf = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl)
+ zalg (ji,jj) = zalf
+ zalfq = zalf * zalf
+ zalf1 = 1.0 - zalf
+ zalg1 (ji,jj) = zalf1
+ zalf1q = zalf1 * zalf1
+ zalg1q(ji,jj) = zalf1q
+ !
+ zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj,jl)
+ zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj,jl) &
+ & - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) )
+ zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) )
+ zfxx (ji,jj) = zfxx(ji,jj) + zalf * ( psxx(ji+1,jj,jl) * zalfq )
+ zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) )
+ zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj,jl)
+ zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj,jl)
+ END_2D
- ! ! Flux from i+1 to i IF u LT 0.
- zbt = zbet(ji,jj)
- zbt1 = 1.0 - zbet(ji,jj)
- zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) )
- zalf = zbt1 * zfm(ji,jj) / zpsm
- zalf1 = 1.0 - zalf
- ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj)
- !
- zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji,jj) )
- zpsx = zbt * zpsx + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * zpsx + 3.0 * ztemp )
- zpsxx = zbt * zpsxx + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * zpsxx &
- & + 5.0 * ( zalf * zalf1 * ( - zpsx + zfx(ji,jj) ) &
- & + ( zalf1 - zalf ) * ztemp ) )
- zpsxy = zbt * zpsxy + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * zpsxy &
- & + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * zpsy ) )
- zpsy = zbt * zpsy + zbt1 * ( zpsy + zfy (ji,jj) )
- zpsyy = zbt * zpsyy + zbt1 * ( zpsyy + zfyy(ji,jj) )
- !
- psm (ji,jj,jl) = zpsm ! optimization
- ps0 (ji,jj,jl) = zps0
- psx (ji,jj,jl) = zpsx
- psxx(ji,jj,jl) = zpsxx
- psy (ji,jj,jl) = zpsy
- psyy(ji,jj,jl) = zpsyy
- psxy(ji,jj,jl) = zpsxy
- END DO
+ DO_2D( 0, 0, jj0, jj0 )
!
- END DO
+ zpsm = psm (ji,jj,jl) ! optimization
+ zps0 = ps0 (ji,jj,jl)
+ zpsx = psx (ji,jj,jl)
+ zpsxx = psxx(ji,jj,jl)
+ zpsy = psy (ji,jj,jl)
+ zpsyy = psyy(ji,jj,jl)
+ zpsxy = psxy(ji,jj,jl)
+ ! ! Readjust moments remaining in the box.
+ zbt = zbet(ji-1,jj)
+ zbt1 = 1.0 - zbet(ji-1,jj)
+ !
+ zpsm = zbt * zpsm + zbt1 * ( zpsm - zfm(ji-1,jj) )
+ zps0 = zbt * zps0 + zbt1 * ( zps0 - zf0(ji-1,jj) )
+ zpsx = zalg1q(ji-1,jj) * ( zpsx + 3.0 * zalg(ji-1,jj) * zpsxx )
+ zpsxx = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * zpsxx
+ zpsy = zbt * zpsy + zbt1 * ( zpsy - zfy (ji-1,jj) )
+ zpsyy = zbt * zpsyy + zbt1 * ( zpsyy - zfyy(ji-1,jj) )
+ zpsxy = zalg1q(ji-1,jj) * zpsxy
+
+ ! Put the temporary moments into appropriate neighboring boxes.
+ ! ! Flux from i to i+1 IF u GT 0.
+ zbt = zbet(ji-1,jj)
+ zbt1 = 1.0 - zbet(ji-1,jj)
+ zpsm = zbt1 * zpsm + zbt * ( zpsm + zfm(ji-1,jj) )
+ zalf = zbt * zfm(ji-1,jj) / zpsm
+ zalf1 = 1.0 - zalf
+ ztemp = zalf * zps0 - zalf1 * zf0(ji-1,jj)
+ !
+ zps0 = zbt1 * zps0 + zbt * ( zps0 + zf0(ji-1,jj) )
+ zpsx = zbt1 * zpsx + zbt * ( ( zalf * zfx(ji-1,jj) + zalf1 * zpsx ) + 3.0 * ztemp )
+ zpsxx = zbt1 * zpsxx + zbt * ( ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * zpsxx ) &
+ & + 5.0 * ( zalf * zalf1 * ( zpsx - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) )
+ zpsxy = zbt1 * zpsxy + zbt * ( ( zalf * zfxy(ji-1,jj) + zalf1 * zpsxy ) &
+ & + 3.0 * ( -zalf1 * zfy(ji-1,jj) + zalf * zpsy ) )
+ zpsy = zbt1 * zpsy + zbt * ( zpsy + zfy (ji-1,jj) )
+ zpsyy = zbt1 * zpsyy + zbt * ( zpsyy + zfyy(ji-1,jj) )
+
+ ! ! Flux from i+1 to i IF u LT 0.
+ zbt = zbet(ji,jj)
+ zbt1 = 1.0 - zbet(ji,jj)
+ zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji ,jj) )
+ zalf = zbt1 * zfm(ji ,jj) / zpsm
+ zalf1 = 1.0 - zalf
+ ztemp = -zalf * zps0 + zalf1 * zf0(ji ,jj)
+ !
+ zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji ,jj) )
+ zpsx = zbt * zpsx + zbt1 * ( ( zalf * zfx(ji ,jj) + zalf1 * zpsx ) + 3.0 * ztemp )
+ zpsxx = zbt * zpsxx + zbt1 * ( ( zalf * zalf * zfxx(ji ,jj) + zalf1 * zalf1 * zpsxx ) &
+ & + 5.0 * ( zalf * zalf1 * ( -zpsx + zfx(ji ,jj) ) + ( zalf1 - zalf ) * ztemp ) )
+ zpsxy = zbt * zpsxy + zbt1 * ( ( zalf * zfxy(ji ,jj) + zalf1 * zpsxy ) &
+ & + 3.0 * ( zalf1 * zfy(ji ,jj) - zalf * zpsy ) )
+ zpsy = zbt * zpsy + zbt1 * ( zpsy + zfy (ji ,jj) )
+ zpsyy = zbt * zpsyy + zbt1 * ( zpsyy + zfyy(ji ,jj) )
+ !
+ psm (ji,jj,jl) = zpsm ! optimization
+ ps0 (ji,jj,jl) = zps0
+ psx (ji,jj,jl) = zpsx
+ psxx(ji,jj,jl) = zpsxx
+ psy (ji,jj,jl) = zpsy
+ psyy(ji,jj,jl) = zpsyy
+ psxy(ji,jj,jl) = zpsxy
+ !
+ END_2D
!
END DO
!
@@ -601,14 +596,14 @@ CONTAINS
zslpmax = MAX( 0._wp, zps0 )
zs1max = 1.5 * zslpmax
zs1new = MIN( zs1max, MAX( -zs1max, zpsy ) )
- zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new )-zslpmax, zpsyy ) )
+ zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new ) - zslpmax, zpsyy ) )
rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask
!
zps0 = zslpmax
- zpsx = zpsx * rswitch
- zpsxx = zpsxx * rswitch
- zpsy = zs1new * rswitch
- zpsyy = zs2new * rswitch
+ zpsx = zpsx * rswitch
+ zpsxx = zpsxx * rswitch
+ zpsy = zs1new * rswitch
+ zpsyy = zs2new * rswitch
zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch
! Calculate fluxes and moments between boxes j<-->j+1
@@ -619,19 +614,19 @@ CONTAINS
zalf1 = 1.0 - zalf
zalf1q = zalf1 * zalf1
!
- zfm (ji,jj) = zalf * zpsm
- zf0 (ji,jj) = zalf * ( zps0 + zalf1 * ( zpsy + (zalf1-zalf) * zpsyy ) )
- zfy (ji,jj) = zalfq *( zpsy + 3.0*zalf1*zpsyy )
- zfyy(ji,jj) = zalf * zalfq * zpsyy
- zfx (ji,jj) = zalf * ( zpsx + zalf1 * zpsxy )
- zfxy(ji,jj) = zalfq * zpsxy
- zfxx(ji,jj) = zalf * zpsxx
+ zfm (ji,jj) = zalf * zpsm
+ zf0 (ji,jj) = zalf * ( zps0 + zalf1 * ( zpsy + (zalf1 - zalf) * zpsyy ) )
+ zfy (ji,jj) = zalfq * ( zpsy + 3.0 * zalf1 * zpsyy )
+ zfyy(ji,jj) = zalf * zpsyy * zalfq
+ zfx (ji,jj) = zalf * ( zpsx + zalf1 * zpsxy )
+ zfxy(ji,jj) = zalfq * zpsxy
+ zfxx(ji,jj) = zalf * zpsxx
!
! ! Readjust moments remaining in the box.
zpsm = zpsm - zfm(ji,jj)
zps0 = zps0 - zf0(ji,jj)
zpsy = zalf1q * ( zpsy -3.0 * zalf * zpsyy )
- zpsyy = zalf1 * zalf1q * zpsyy
+ zpsyy = zalf1 * zalf1q * zpsyy
zpsx = zpsx - zfx(ji,jj)
zpsxx = zpsxx - zfxx(ji,jj)
zpsxy = zalf1q * zpsxy
@@ -659,16 +654,13 @@ CONTAINS
zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji,jj+1,jl) &
& - zalf1 * (psy(ji,jj+1,jl) - (zalf1 - zalf ) * psyy(ji,jj+1,jl) ) )
zfy (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1,jl) - 3.0 * zalf1 * psyy(ji,jj+1,jl) )
- zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji,jj+1,jl) * zalfq
+ zfyy (ji,jj) = zfyy(ji,jj) + zalf * ( psyy(ji,jj+1,jl) * zalfq )
zfx (ji,jj) = zfx (ji,jj) + zalf * ( psx (ji,jj+1,jl) - zalf1 * psxy(ji,jj+1,jl) )
zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1,jl)
zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji,jj+1,jl)
END_2D
DO_2D( ji0, ji0, 0, 0 )
- ! ! Readjust moments remaining in the box.
- zbt = zbet(ji,jj-1)
- zbt1 = ( 1.0 - zbet(ji,jj-1) )
!
zpsm = psm (ji,jj,jl) ! optimization
zps0 = ps0 (ji,jj,jl)
@@ -677,53 +669,52 @@ CONTAINS
zpsy = psy (ji,jj,jl)
zpsyy = psyy(ji,jj,jl)
zpsxy = psxy(ji,jj,jl)
+ ! ! Readjust moments remaining in the box.
+ zbt = zbet(ji,jj-1)
+ zbt1 = ( 1.0 - zbet(ji,jj-1) )
!
- zpsm = zbt * zpsm + zbt1 * ( zpsm - zfm(ji,jj-1) )
- zps0 = zbt * zps0 + zbt1 * ( zps0 - zf0(ji,jj-1) )
+ zpsm = zbt * zpsm + zbt1 * ( zpsm - zfm(ji,jj-1) )
+ zps0 = zbt * zps0 + zbt1 * ( zps0 - zf0(ji,jj-1) )
zpsy = zalg1q(ji,jj-1) * ( zpsy + 3.0 * zalg(ji,jj-1) * zpsyy )
zpsyy = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * zpsyy
- zpsx = zbt * zpsx + zbt1 * ( zpsx - zfx (ji,jj-1) )
- zpsxx = zbt * zpsxx + zbt1 * ( zpsxx - zfxx(ji,jj-1) )
+ zpsx = zbt * zpsx + zbt1 * ( zpsx - zfx (ji,jj-1) )
+ zpsxx = zbt * zpsxx + zbt1 * ( zpsxx - zfxx(ji,jj-1) )
zpsxy = zalg1q(ji,jj-1) * zpsxy
! Put the temporary moments into appropriate neighboring boxes.
! ! Flux from j to j+1 IF v GT 0.
zbt = zbet(ji,jj-1)
zbt1 = 1.0 - zbet(ji,jj-1)
- zpsm = zbt * ( zpsm + zfm(ji,jj-1) ) + zbt1 * zpsm
- zalf = zbt * zfm(ji,jj-1) / zpsm
+ zpsm = zbt1 * zpsm + zbt * ( zpsm + zfm(ji,jj-1) )
+ zalf = zbt * zfm(ji,jj-1) / zpsm
zalf1 = 1.0 - zalf
- ztemp = zalf * zps0 - zalf1 * zf0(ji,jj-1)
+ ztemp = zalf * zps0 - zalf1 * zf0(ji,jj-1)
!
- zps0 = zbt * ( zps0 + zf0(ji,jj-1) ) + zbt1 * zps0
- zpsy = zbt * ( zalf * zfy(ji,jj-1) + zalf1 * zpsy + 3.0 * ztemp ) &
- & + zbt1 * zpsy
- zpsyy = zbt * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * zpsyy &
- & + 5.0 * ( zalf * zalf1 * ( zpsy - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &
- & + zbt1 * zpsyy
- zpsxy = zbt * ( zalf * zfxy(ji,jj-1) + zalf1 * zpsxy &
- & + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * zpsx ) ) &
- & + zbt1 * zpsxy
- zpsx = zbt * ( zpsx + zfx (ji,jj-1) ) + zbt1 * zpsx
- zpsxx = zbt * ( zpsxx + zfxx(ji,jj-1) ) + zbt1 * zpsxx
+ zps0 = zbt1 * zps0 + zbt * ( zps0 + zf0(ji,jj-1) )
+ zpsy = zbt1 * zpsy + zbt * ( ( zalf * zfy(ji,jj-1) + zalf1 * zpsy ) + 3.0 * ztemp )
+ zpsyy = zbt1 * zpsyy + zbt * ( ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * zpsyy ) &
+ & + 5.0 * ( zalf * zalf1 * ( zpsy - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )
+ zpsxy = zbt1 * zpsxy + zbt * ( ( zalf * zfxy(ji,jj-1) + zalf1 * zpsxy ) &
+ & + 3.0 * ( -zalf1 * zfx(ji,jj-1) + zalf * zpsx ) )
+ zpsx = zbt1 * zpsx + zbt * ( zpsx + zfx (ji,jj-1) )
+ zpsxx = zbt1 * zpsxx + zbt * ( zpsxx + zfxx(ji,jj-1) )
! ! Flux from j+1 to j IF v LT 0.
zbt = zbet(ji,jj)
zbt1 = 1.0 - zbet(ji,jj)
- zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) )
- zalf = zbt1 * zfm(ji,jj) / zpsm
+ zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) )
+ zalf = zbt1 * zfm(ji,jj) / zpsm
zalf1 = 1.0 - zalf
- ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj)
+ ztemp = -zalf * zps0 + zalf1 * zf0(ji,jj)
!
- zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji,jj) )
- zpsy = zbt * zpsy + zbt1 * ( zalf * zfy(ji,jj) + zalf1 * zpsy + 3.0 * ztemp )
- zpsyy = zbt * zpsyy + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * zpsyy &
- & + 5.0 * ( zalf * zalf1 * ( - zpsy + zfy(ji,jj) ) &
- & + ( zalf1 - zalf ) * ztemp ) )
- zpsxy = zbt * zpsxy + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * zpsxy &
- & + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * zpsx ) )
- zpsx = zbt * zpsx + zbt1 * ( zpsx + zfx (ji,jj) )
- zpsxx = zbt * zpsxx + zbt1 * ( zpsxx + zfxx(ji,jj) )
+ zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji,jj ) )
+ zpsy = zbt * zpsy + zbt1 * ( ( zalf * zfy(ji,jj ) + zalf1 * zpsy ) + 3.0 * ztemp )
+ zpsyy = zbt * zpsyy + zbt1 * ( ( zalf * zalf * zfyy(ji,jj ) + zalf1 * zalf1 * zpsyy ) &
+ & + 5.0 * ( zalf * zalf1 * ( - zpsy + zfy(ji,jj ) ) + ( zalf1 - zalf ) * ztemp ) )
+ zpsxy = zbt * zpsxy + zbt1 * ( ( zalf * zfxy(ji,jj ) + zalf1 * zpsxy ) &
+ & + 3.0 * ( zalf1 * zfx(ji,jj ) - zalf * zpsx ) )
+ zpsx = zbt * zpsx + zbt1 * ( zpsx + zfx (ji,jj ) )
+ zpsxx = zbt * zpsxx + zbt1 * ( zpsxx + zfxx(ji,jj ) )
!
psm (ji,jj,jl) = zpsm ! optimization
ps0 (ji,jj,jl) = zps0
diff --git a/src/ICE/icedyn_adv_umx.F90 b/src/ICE/icedyn_adv_umx.F90
index d46b056b8fb94c4e8dbaaa5530e81e4390d050d8..f15c2349f2a5f7fa8df0ab4d6493a3f8b2ffc86a 100644
--- a/src/ICE/icedyn_adv_umx.F90
+++ b/src/ICE/icedyn_adv_umx.F90
@@ -375,7 +375,8 @@ CONTAINS
zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
DO_2D( 0, 0, 0, 0 )
pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &
- & - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
+ & - ( ( zudy(ji,jj) - zudy(ji-1,jj) ) & ! ad () for NP repro
+ & + ( zvdx(ji,jj) - zvdx(ji,jj-1) ) ) * r1_e1e2t(ji,jj) * zdt
END_2D
CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1._wp )
!
@@ -516,7 +517,8 @@ CONTAINS
! thus we calculate the upstream solution and apply a limiter again
DO jl = 1, jpl
DO_2D( 0, 0, 0, 0 )
- ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) + zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) )
+ ztra = - ( ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) & ! add () for NP repro
+ & + ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) )
!
zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1)
END_2D
@@ -551,7 +553,8 @@ CONTAINS
! ---------------------------------
DO jl = 1, jpl
DO_2D( 0, 0, 0, 0 )
- ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) )
+ ztra = - ( ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) ) & ! add () for NP repro
+ & + ( zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) )
!
ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1)
END_2D
@@ -643,10 +646,10 @@ CONTAINS
!
DO jl = 1, jpl !-- after tracer with upstream scheme
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
- ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj ,jl) &
- & + pfv_ups(ji,jj,jl) - pfv_ups(ji ,jj-1,jl) ) &
- & + ( pu (ji,jj ) - pu (ji-1,jj ) &
- & + pv (ji,jj ) - pv (ji ,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk)
+ ztra = - ( ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj ,jl) ) & ! add () for NP repro
+ & + ( pfv_ups(ji,jj,jl) - pfv_ups(ji ,jj-1,jl) ) ) &
+ & + ( ( pu (ji,jj ) - pu (ji-1,jj ) ) &
+ & + ( pv (ji,jj ) - pv (ji ,jj-1 ) ) ) * pt(ji,jj,jl) * (1.-pamsk)
!
pt_ups(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
END_2D
@@ -910,7 +913,7 @@ CONTAINS
!
DO jl = 1, jpl
DO_2D( 1, 0, kloop, kloop )
- pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) &
+ pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt(ji+1,jj,jl) + pt(ji,jj,jl) ) &
& - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
END_2D
END DO
@@ -920,7 +923,7 @@ CONTAINS
DO jl = 1, jpl
DO_2D( 1, 0, kloop, kloop )
zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
- pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) &
+ pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt(ji+1,jj,jl) + pt(ji,jj,jl) ) &
& - zcu * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
END_2D
END DO
@@ -932,9 +935,9 @@ CONTAINS
zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
zdx2 = e1u(ji,jj) * e1u(ji,jj)
!!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)
- pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) &
+ pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) ) &
& - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) &
- & + r1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) &
+ & + r1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) ) &
& - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) )
END_2D
END DO
@@ -946,9 +949,9 @@ CONTAINS
zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
zdx2 = e1u(ji,jj) * e1u(ji,jj)
!!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)
- pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) &
+ pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) ) &
& - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) &
- & + r1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) &
+ & + r1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) ) &
& - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) )
END_2D
END DO
@@ -963,11 +966,11 @@ CONTAINS
zdx2 = e1u(ji,jj) * e1u(ji,jj)
!!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj)
zdx4 = zdx2 * zdx2
- pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) &
+ pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) ) &
& - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) &
- & + r1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) &
+ & + r1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) ) &
& - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) &
- & + r1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl) &
+ & + r1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ((ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl) ) &
& - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) )
END_2D
END DO
@@ -981,7 +984,7 @@ CONTAINS
DO jl = 1, jpl
DO_2D( 1, 0, kloop, kloop )
IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN
- pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) &
+ pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt(ji+1,jj,jl) + pt(ji,jj,jl) ) &
& - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
ENDIF
END_2D
@@ -1048,7 +1051,7 @@ CONTAINS
CASE( 1 ) !== 1st order central TIM ==! (Eq. 21)
DO jl = 1, jpl
DO_2D( kloop, kloop, 1, 0 )
- pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) &
+ pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) &
& - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
END_2D
END DO
@@ -1057,7 +1060,7 @@ CONTAINS
DO jl = 1, jpl
DO_2D( kloop, kloop, 1, 0 )
zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
- pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) &
+ pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) &
& - zcv * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
END_2D
END DO
@@ -1068,9 +1071,9 @@ CONTAINS
zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
zdy2 = e2v(ji,jj) * e2v(ji,jj)
!!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)
- pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) &
+ pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) ) &
& - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) &
- & + r1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) &
+ & + r1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) ) &
& - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
END_2D
END DO
@@ -1081,9 +1084,9 @@ CONTAINS
zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
zdy2 = e2v(ji,jj) * e2v(ji,jj)
!!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)
- pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) &
+ pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) ) &
& - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) &
- & + r1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) &
+ & + r1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) ) &
& - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
END_2D
END DO
@@ -1098,11 +1101,11 @@ CONTAINS
zdy2 = e2v(ji,jj) * e2v(ji,jj)
!!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj)
zdy4 = zdy2 * zdy2
- pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) &
+ pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) ) &
& - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) &
- & + r1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) &
+ & + r1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) ) &
& - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) &
- & + r1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl) &
+ & + r1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ((ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl) ) &
& - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) )
END_2D
END DO
@@ -1242,10 +1245,10 @@ CONTAINS
zneg = MAX( 0._wp, pfu_ho(ji ,jj ,jl) ) - MIN( 0._wp, pfu_ho(ji-1,jj ,jl) ) &
& + MAX( 0._wp, pfv_ho(ji ,jj ,jl) ) - MIN( 0._wp, pfv_ho(ji ,jj-1,jl) )
!
- zpos = zpos - (pt(ji,jj,jl) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MIN( 0., pv(ji,jj) - pv(ji,jj-1) ) &
- & ) * ( 1. - pamsk )
- zneg = zneg + (pt(ji,jj,jl) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MAX( 0., pv(ji,jj) - pv(ji,jj-1) ) &
- & ) * ( 1. - pamsk )
+ zpos = zpos - ( pt(ji,jj,jl) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) &
+ & + pt(ji,jj,jl) * MIN( 0., pv(ji,jj) - pv(ji,jj-1) ) ) * ( 1. - pamsk )
+ zneg = zneg + ( pt(ji,jj,jl) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) &
+ & + pt(ji,jj,jl) * MAX( 0., pv(ji,jj) - pv(ji,jj-1) ) ) * ( 1. - pamsk )
!
! ! up & down beta terms
! clem: zbetup and zbetdo must be 0 for zpos>1.e-10 & zneg>1.e-10 (do not put 0 instead of 1.e-10 !!!)
diff --git a/src/ICE/icedyn_rdgrft.F90 b/src/ICE/icedyn_rdgrft.F90
index e083da9c5dab0cb933b894ee34a245e80046c911..214d6234996fd2ad8391cbbd4f8b70e86a0dbe61 100644
--- a/src/ICE/icedyn_rdgrft.F90
+++ b/src/ICE/icedyn_rdgrft.F90
@@ -960,8 +960,8 @@ CONTAINS
DO_2D( 0, 0, 0, 0 )
IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN
zworka(ji,jj) = ( 4._wp * strength(ji,jj) &
- & + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) &
- & + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
+ & + ( ( strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) ) &
+ & + ( strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) ) ) &
& ) / ( 4._wp + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
ELSE
zworka(ji,jj) = 0._wp
diff --git a/src/ICE/icedyn_rhg_eap.F90 b/src/ICE/icedyn_rhg_eap.F90
index 55b1392ac6a5e8ffeda9e8c44c516e5b2c009630..e87e879a335e79f09dbc60fe8e5d4d7ab1a4f11a 100644
--- a/src/ICE/icedyn_rhg_eap.F90
+++ b/src/ICE/icedyn_rhg_eap.F90
@@ -306,9 +306,9 @@ CONTAINS
zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
- ! Ocean currents at U-V points
- v_oceU(ji,jj) = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1)
- u_oceV(ji,jj) = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1)
+ ! Ocean currents at U-V point, warning: add () for the North Pole reproducibility
+ v_oceU(ji,jj) = 0.25_wp * ( ( v_oce(ji,jj) + v_oce(ji,jj-1) ) + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) ) * umask(ji,jj,1)
+ u_oceV(ji,jj) = 0.25_wp * ( ( u_oce(ji,jj) + u_oce(ji-1,jj) ) + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) ) * vmask(ji,jj,1)
! m/dt
zmU_t(ji,jj) = zmassU * z1_dtevp
@@ -390,14 +390,14 @@ CONTAINS
DO_2D( 0, 0, 0, 0 )
- ! shear**2 at T points (doc eq. A16)
- zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) &
- & + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) &
+ ! shear**2 at T points (doc eq. A16), warning: add () for the North Pole reproducibility
+ zds2 = ( ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) ) &
+ & + ( zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) ) &
& ) * 0.25_wp * r1_e1e2t(ji,jj)
- ! divergence at T points
- zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
- & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
+ ! divergence at T points, warning: add () for the North Pole reproducibility
+ zdiv = ( ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) ) &
+ & + ( e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) &
& ) * r1_e1e2t(ji,jj)
zdiv2 = zdiv * zdiv
@@ -418,14 +418,14 @@ CONTAINS
DO_2D( 0, 0, 0, 0 )
- ! shear at T points
- zdsT = ( zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * e1e2f(ji-1,jj ) &
- & + zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) &
+ ! shear at T points, warning: add () for the North Pole reproducibility
+ zdsT = ( ( zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * e1e2f(ji-1,jj ) ) &
+ & + ( zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) ) &
& ) * 0.25_wp * r1_e1e2t(ji,jj)
- ! divergence at T points (duplication to avoid communications)
- zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
- & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
+ ! divergence at T points (duplication to avoid communications), warning: add () for the North Pole reproducibility
+ zdiv = ( ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) ) &
+ & + ( e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) &
& ) * r1_e1e2t(ji,jj)
! tension at T points (duplication to avoid communications)
@@ -760,9 +760,9 @@ CONTAINS
zten_i(ji,jj) = zdt
- ! shear**2 at T points (doc eq. A16)
- zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) &
- & + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) &
+ ! shear**2 at T points (doc eq. A16), warning: add () for the North Pole reproducibility
+ zds2 = ( ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) ) &
+ & + ( zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) ) &
& ) * 0.25_wp * r1_e1e2t(ji,jj)
! maximum shear rate at T points (includes tension, output only)
@@ -771,9 +771,9 @@ CONTAINS
! shear at T-points
zshear(ji,jj) = SQRT( zds2 )
- ! divergence at T points
- pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
- & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
+ ! divergence at T points, warning: add () for the North Pole reproducibility
+ pdivu_i(ji,jj) = ( ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) ) &
+ & + ( e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) &
& ) * r1_e1e2t(ji,jj)
! delta at T points
diff --git a/src/ICE/icedyn_rhg_evp.F90 b/src/ICE/icedyn_rhg_evp.F90
index 49e7753c2ae939a3c6a5a6a9a5ea74665ecdb7bf..136f13d122b5d03a3416906f729bae42f3befb41 100644
--- a/src/ICE/icedyn_rhg_evp.F90
+++ b/src/ICE/icedyn_rhg_evp.F90
@@ -290,7 +290,7 @@ CONTAINS
zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
! Ocean currents at U-V points
- ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility)
+ ! (brackets added to fix the order of floating point operations for the North Pole reproducibility)
v_oceU(ji,jj) = 0.25_wp * ( (v_oce(ji,jj) + v_oce(ji,jj-1)) + (v_oce(ji+1,jj) + v_oce(ji+1,jj-1)) ) * umask(ji,jj,1)
u_oceV(ji,jj) = 0.25_wp * ( (u_oce(ji,jj) + u_oce(ji-1,jj)) + (u_oce(ji,jj+1) + u_oce(ji-1,jj+1)) ) * vmask(ji,jj,1)
@@ -377,13 +377,13 @@ CONTAINS
DO_2D( 0, 0, 0, 0 )
! shear**2 at T points (doc eq. A16)
- zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) &
- & + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) &
+ zds2 = ( ( zds(ji,jj )*zds(ji,jj )*e1e2f(ji,jj ) + zds(ji-1,jj )*zds(ji-1,jj )*e1e2f(ji-1,jj ) ) & ! add () for
+ & + ( zds(ji,jj-1)*zds(ji,jj-1)*e1e2f(ji,jj-1) + zds(ji-1,jj-1)*zds(ji-1,jj-1)*e1e2f(ji-1,jj-1) ) & ! NP repro
& ) * 0.25_wp * r1_e1e2t(ji,jj)
! divergence at T points
- zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
- & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
+ zdiv = ( ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) ) & ! add () for
+ & + ( e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) & ! NP repro
& ) * r1_e1e2t(ji,jj)
zdiv2 = zdiv * zdiv
@@ -406,9 +406,8 @@ CONTAINS
DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls ) ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2
! divergence at T points (duplication to avoid communications)
- ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility)
- zdiv = ( (e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)) &
- & + (e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)) &
+ zdiv = ( ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) ) & ! add () for
+ & + ( e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) & ! NP repro
& ) * r1_e1e2t(ji,jj)
! tension at T points (duplication to avoid communications)
@@ -458,8 +457,8 @@ CONTAINS
ENDIF
! P/delta at F points
- ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility)
- zp_delf = 0.25_wp * ( (zp_delt(ji,jj) + zp_delt(ji+1,jj)) + (zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1)) )
+ ! (brackets added to fix the order of floating point operations for the North Pole reproducibility)
+ zp_delf = 0.25_wp * ( ( zp_delt(ji,jj) + zp_delt(ji+1,jj) ) + ( zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) ) )
! stress at F points (zkt/=0 if landfast)
zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) &
@@ -471,24 +470,24 @@ CONTAINS
! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility)
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
! !--- U points
- zfU(ji,jj) = 0.5_wp * ( (( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) &
- & + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) &
- & ) * r1_e2u(ji,jj)) &
+ zfU(ji,jj) = 0.5_wp * ( ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) &
+ & + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) &
+ & ) * r1_e2u(ji,jj) ) &
& + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) &
& ) * 2._wp * r1_e1u(ji,jj) &
& ) * r1_e1e2u(ji,jj)
!
! !--- V points
- zfV(ji,jj) = 0.5_wp * ( (( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) &
- & - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) &
- & ) * r1_e1v(ji,jj)) &
+ zfV(ji,jj) = 0.5_wp * ( ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) &
+ & - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) &
+ & ) * r1_e1v(ji,jj) ) &
& + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) &
& ) * 2._wp * r1_e2v(ji,jj) &
& ) * r1_e1e2v(ji,jj)
!
- ! !--- ice currents at U-V point
- v_iceU(ji,jj) = 0.25_wp * ( (v_ice(ji,jj) + v_ice(ji,jj-1)) + (v_ice(ji+1,jj) + v_ice(ji+1,jj-1)) ) * umask(ji,jj,1)
- u_iceV(ji,jj) = 0.25_wp * ( (u_ice(ji,jj) + u_ice(ji-1,jj)) + (u_ice(ji,jj+1) + u_ice(ji-1,jj+1)) ) * vmask(ji,jj,1)
+ ! !--- ice currents at U-V point, warning: add () for NP repro
+ v_iceU(ji,jj) = 0.25_wp * ( ( v_ice(ji,jj) + v_ice(ji,jj-1) ) + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) ) * umask(ji,jj,1)
+ u_iceV(ji,jj) = 0.25_wp * ( ( u_ice(ji,jj) + u_ice(ji-1,jj) ) + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) ) * vmask(ji,jj,1)
!
END_2D
!
@@ -755,8 +754,8 @@ CONTAINS
zten_i(ji,jj) = zdt
! shear**2 at T points (doc eq. A16)
- zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) &
- & + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) &
+ zds2 = ( ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) ) & ! add ()
+ & + ( zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) ) & ! NP rep
& ) * 0.25_wp * r1_e1e2t(ji,jj)
! maximum shear rate at T points (includes tension, output only)
@@ -766,8 +765,8 @@ CONTAINS
zshear(ji,jj) = SQRT( zds2 ) * zmsk(ji,jj)
! divergence at T points
- pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
- & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
+ pdivu_i(ji,jj) = ( ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) ) & ! add () for NP repro
+ & + ( e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) &
& ) * r1_e1e2t(ji,jj) * zmsk(ji,jj)
! delta at T points
@@ -1127,7 +1126,8 @@ CONTAINS
ENDIF
!
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
- ztra = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) + zfv_ups(ji,jj) - zfv_ups(ji,jj-1) )
+ ztra = - ( ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) & ! add () for NP repro
+ & + ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) )
pt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
END_2D
END DO
diff --git a/src/ICE/icedyn_rhg_vp.F90 b/src/ICE/icedyn_rhg_vp.F90
index b2c344993dbf9b471bc971d48c4bcf9301560228..1553424ea429d9093bdefb34aeca86ab0887dddc 100644
--- a/src/ICE/icedyn_rhg_vp.F90
+++ b/src/ICE/icedyn_rhg_vp.F90
@@ -328,7 +328,7 @@ CONTAINS
zmU_t(ji,jj) = zmassU_t(ji,jj) * u_ice(ji,jj)
! Ocean currents at U-V points
- ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility)
+ ! (brackets added to fix the order of floating point operations for the North Pole reproducibility
v_oceU(ji,jj) = 0.25_wp * ( (v_oce(ji,jj) + v_oce(ji,jj-1)) + (v_oce(ji+1,jj) + v_oce(ji+1,jj-1)) ) * umask(ji,jj,1)
! Wind stress
@@ -368,7 +368,7 @@ CONTAINS
zmV_t(ji,jj) = zmassV_t(ji,jj) * v_ice(ji,jj)
! Ocean currents at U-V points
- ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility)
+ ! (brackets added to fix the order of floating point operations for the North Pole reproducibility
u_oceV(ji,jj) = 0.25_wp * ( (u_oce(ji,jj) + u_oce(ji-1,jj)) + (u_oce(ji,jj+1) + u_oce(ji-1,jj+1)) ) * vmask(ji,jj,1)
! Wind stress
@@ -436,14 +436,15 @@ CONTAINS
! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
! shear**2 at T points (doc eq. A16)
- zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) &
- & + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) &
+ ! (brackets added to fix the order of floating point operations for the North Pole reproducibility
+ zds2 = ( ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) ) &
+ & + ( zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) ) &
& ) * 0.25_wp * r1_e1e2t(ji,jj)
! divergence at T points
- ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility)
- zdiv = ( (e2u(ji,jj) * zu_c(ji,jj) - e2u(ji-1,jj) * zu_c(ji-1,jj)) &
- & + (e1v(ji,jj) * zv_c(ji,jj) - e1v(ji,jj-1) * zv_c(ji,jj-1)) &
+ ! (brackets added to fix the order of floating point operations for the North Pole reproducibility
+ zdiv = ( ( e2u(ji,jj) * zu_c(ji,jj) - e2u(ji-1,jj) * zu_c(ji-1,jj) ) &
+ & + ( e1v(ji,jj) * zv_c(ji,jj) - e1v(ji,jj-1) * zv_c(ji,jj-1) ) &
& ) * r1_e1e2t(ji,jj)
zdiv2 = zdiv * zdiv
@@ -474,7 +475,7 @@ CONTAINS
DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )! 1-> jpj-1; 1->jpi-1
! P/delta* at F points
- ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility)
+ ! (brackets added to fix the order of floating point operations for the North Pole reproducibility
zvisc_f = 0.25_wp * ( (zvisc_t(ji,jj) + zvisc_t(ji+1,jj)) + (zvisc_t(ji,jj+1) + zvisc_t(ji+1,jj+1)) )
! Temporary zef factor at F-point
@@ -489,7 +490,7 @@ CONTAINS
DO_2D( nn_hls, nn_hls-1, nn_hls-1, nn_hls )
!--- ice u-velocity @V points, v-velocity @U points (for non-linear drag computation)
- ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility)
+ ! (brackets added to fix the order of floating point operations for the North Pole reproducibility
zv_cU = 0.25_wp * ( (zv_c(ji,jj) + zv_c(ji,jj-1)) + (zv_c(ji+1,jj) + zv_c(ji+1,jj-1)) ) * umask(ji,jj,1)
!--- non-linear drag coefficients (need to be updated at each outer loop, see Lemieux and Tremblay JGR09, p.3, beginning of Section 3)
@@ -509,7 +510,7 @@ CONTAINS
DO_2D( nn_hls-1, nn_hls, nn_hls, nn_hls-1 )
!--- ice u-velocity @V points, v-velocity @U points (for non-linear drag computation)
- ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility)
+ ! (brackets added to fix the order of floating point operations for the North Pole reproducibility
zu_cV = 0.25_wp * ( (zu_c(ji,jj) + zu_c(ji-1,jj)) + (zu_c(ji,jj+1) + zu_c(ji-1,jj+1)) ) * vmask(ji,jj,1)
!--- non-linear drag coefficients (need to be updated at each outer loop, see Lemieux and Tremblay JGR09, p.3, beginning of Section 3)
@@ -1066,8 +1067,9 @@ CONTAINS
DO_2D( 0, 0, 0, 0 ) ! 2->jpj-1; 2->jpi-1
! shear**2 at T points (doc eq. A16)
- zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) &
- & + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) &
+ ! (brackets added to fix the order of floating point operations for the North Pole reproducibility
+ zds2 = ( ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) ) &
+ & + ( zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) ) &
& ) * 0.25_wp * r1_e1e2t(ji,jj)
! tension**2 at T points
@@ -1085,8 +1087,9 @@ CONTAINS
zshear(ji,jj) = SQRT( zds2 ) * zmsk(ji,jj)
! divergence at T points
- pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
- & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
+ ! (brackets added to fix the order of floating point operations for the North Pole reproducibility
+ pdivu_i(ji,jj) = ( ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) ) &
+ & + ( e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) &
& ) * r1_e1e2t(ji,jj) * zmsk(ji,jj)
! delta at T points
@@ -1157,22 +1160,22 @@ CONTAINS
!--- Recalculate oceanic stress at last inner iteration
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1
- !--- ice u-velocity @V points, v-velocity @U points (for non-linear drag computation)
- zu_cV = 0.25_wp * ( u_ice(ji,jj) + u_ice(ji-1,jj) + u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * vmask(ji,jj,1)
- zv_cU = 0.25_wp * ( v_ice(ji,jj) + v_ice(ji,jj-1) + v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * umask(ji,jj,1)
-
- !--- non-linear drag coefficients (need to be updated at each outer loop, see Lemieux and Tremblay JGR09, p.3, beginning of Section 3)
- zCwU(ji,jj) = za_iU(ji,jj) * zrhoco * SQRT( ( u_ice(ji,jj) - u_oce (ji,jj) ) * ( u_ice(ji,jj) - u_oce (ji,jj) ) &
- & + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) )
- zCwV(ji,jj) = za_iV(ji,jj) * zrhoco * SQRT( ( v_ice(ji,jj) - v_oce (ji,jj) ) * ( v_ice(ji,jj) - v_oce (ji,jj) ) &
- & + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) )
-
- !--- Ocean-ice stress
- ztaux_oi(ji,jj) = zCwU(ji,jj) * ( u_oce(ji,jj) - u_ice(ji,jj) )
- ztauy_oi(ji,jj) = zCwV(ji,jj) * ( v_oce(ji,jj) - v_ice(ji,jj) )
-
+ !--- ice u-velocity @V points, v-velocity @U points (for non-linear drag computation)
+ ! (brackets added to fix the order of floating point operations for the North Pole reproducibility
+ zu_cV = 0.25_wp * ( ( u_ice(ji,jj) + u_ice(ji-1,jj) ) + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) ) * vmask(ji,jj,1)
+ zv_cU = 0.25_wp * ( ( v_ice(ji,jj) + v_ice(ji,jj-1) ) + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) ) * umask(ji,jj,1)
+
+ !--- non-linear drag coefficients (need to be updated at each outer loop, see Lemieux and Tremblay JGR09, p.3, beginning of Section 3)
+ zCwU(ji,jj) = za_iU(ji,jj) * zrhoco * SQRT( ( u_ice(ji,jj) - u_oce (ji,jj) ) * ( u_ice(ji,jj) - u_oce (ji,jj) ) &
+ & + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) )
+ zCwV(ji,jj) = za_iV(ji,jj) * zrhoco * SQRT( ( v_ice(ji,jj) - v_oce (ji,jj) ) * ( v_ice(ji,jj) - v_oce (ji,jj) ) &
+ & + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) )
+
+ !--- Ocean-ice stress
+ ztaux_oi(ji,jj) = zCwU(ji,jj) * ( u_oce(ji,jj) - u_ice(ji,jj) )
+ ztauy_oi(ji,jj) = zCwV(ji,jj) * ( v_oce(ji,jj) - v_ice(ji,jj) )
+
END_2D
-
!
CALL lbc_lnk( 'icedyn_rhg_vp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1. ) !, &
! & ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. )
@@ -1567,17 +1570,18 @@ CONTAINS
DO_2D( 0, 0, 0, 0 ) !clem check bounds
- zu_res(ji,jj) = ( prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1) &
- & - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) )
- zv_res(ji,jj) = ( prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj) &
- & - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) )
-
-! zu_res(ji,jj) = pFU(ji,jj) - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj)
-! zv_res(ji,jj) = pFV(ji,jj) - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1)
-
- zu_res(ji,jj) = SQRT( zu_res(ji,jj) * zu_res(ji,jj) ) * umask(ji,jj,1) * pat_iu(ji,jj) * e1e2u(ji,jj) * z1_pglob_area
- zv_res(ji,jj) = SQRT( zv_res(ji,jj) * zv_res(ji,jj) ) * vmask(ji,jj,1) * pat_iv(ji,jj) * e1e2v(ji,jj) * z1_pglob_area
-
+ ! (brackets added to fix the order of floating point operations for the North Pole reproducibility
+ zu_res(ji,jj) = prhsu(ji,jj) + ( pDU(ji,jj) * pu(ji ,jj-1) + pEU(ji,jj) * pu(ji ,jj+1) ) &
+ & - pBU(ji,jj) * pu(ji,jj) - ( pAU(ji,jj) * pu(ji-1,jj ) + pCU(ji,jj) * pu(ji+1,jj ) )
+ zv_res(ji,jj) = prhsv(ji,jj) + ( pDV(ji,jj) * pv(ji-1,jj ) + pEV(ji,jj) * pv(ji+1,jj ) ) &
+ & - pBV(ji,jj) * pv(ji,jj) - ( pAV(ji,jj) * pv(ji ,jj-1) + pCV(ji,jj) * pv(ji ,jj+1) )
+
+ ! zu_res(ji,jj) = pFU(ji,jj) - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj)
+ ! zv_res(ji,jj) = pFV(ji,jj) - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1)
+
+ zu_res(ji,jj) = SQRT( zu_res(ji,jj) * zu_res(ji,jj) ) * umask(ji,jj,1) * pat_iu(ji,jj) * e1e2u(ji,jj) * z1_pglob_area
+ zv_res(ji,jj) = SQRT( zv_res(ji,jj) * zv_res(ji,jj) ) * vmask(ji,jj,1) * pat_iv(ji,jj) * e1e2v(ji,jj) * z1_pglob_area
+
END_2D
! Global ice-concentration, grid-cell-area weighted mean
@@ -1607,21 +1611,22 @@ CONTAINS
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
- zu_res(ji,jj) = ( prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1) &
- & - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) )
- zv_res(ji,jj) = ( prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj) &
- & - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) )
-
- zu_res(ji,jj) = SQRT( zu_res(ji,jj) * zu_res(ji,jj) ) * umask(ji,jj,1)
- zv_res(ji,jj) = SQRT( zv_res(ji,jj) * zv_res(ji,jj) ) * vmask(ji,jj,1)
-
+ ! (brackets added to fix the order of floating point operations for the North Pole reproducibility
+ zu_res(ji,jj) = prhsu(ji,jj) + ( pDU(ji,jj) * pu(ji ,jj-1) + pEU(ji,jj) * pu(ji ,jj+1) ) &
+ & - pBU(ji,jj) * pu(ji,jj) - ( pAU(ji,jj) * pu(ji-1,jj ) + pCU(ji,jj) * pu(ji+1,jj ) )
+ zv_res(ji,jj) = prhsv(ji,jj) + ( pDV(ji,jj) * pv(ji-1,jj ) + pEV(ji,jj) * pv(ji+1,jj ) ) &
+ & - pBV(ji,jj) * pv(ji,jj) - ( pAV(ji,jj) * pv(ji ,jj-1) + pCV(ji,jj) * pv(ji ,jj+1) )
+
+ zu_res(ji,jj) = SQRT( zu_res(ji,jj) * zu_res(ji,jj) ) * umask(ji,jj,1)
+ zv_res(ji,jj) = SQRT( zv_res(ji,jj) * zv_res(ji,jj) ) * vmask(ji,jj,1)
+
END_2D
IF( nn_hls == 1 ) CALL lbc_lnk( 'icedyn_rhg_cvg_vp', zu_res, 'U', 1., zv_res , 'V', 1. )
DO_2D( 0, 0, 0, 0 ) !clem check bounds
- pvel_res(ji,jj) = 0.25_wp * ( zu_res(ji-1,jj) + zu_res(ji,jj) + zv_res(ji,jj-1) + zv_res(ji,jj) )
+ pvel_res(ji,jj) = 0.25_wp * ( ( zu_res(ji-1,jj) + zu_res(ji,jj) ) + ( zv_res(ji,jj-1) + zv_res(ji,jj) ) )
END_2D
CALL lbc_lnk( 'icedyn_rhg_cvg_vp', pvel_res, 'T', 1. )
diff --git a/src/ICE/icesbc.F90 b/src/ICE/icesbc.F90
index d9406ce66d2b685a391ee1eb45ab06889c0574f2..1efec20999ec2afcb1d7f6103d63025536862e56 100644
--- a/src/ICE/icesbc.F90
+++ b/src/ICE/icesbc.F90
@@ -308,7 +308,8 @@ CONTAINS
zv_io = v_ice(ji ,jj ) - ssv_m(ji ,jj )
zv_iom1 = v_ice(ji ,jj-1) - ssv_m(ji ,jj-1)
!
- zfric(ji,jj) = rn_cio * ( 0.5_wp * ( zu_io*zu_io + zu_iom1*zu_iom1 + zv_io*zv_io + zv_iom1*zv_iom1 ) ) * smask0(ji,jj)
+ zfric(ji,jj) = rn_cio * ( 0.5_wp * ( ( zu_io*zu_io + zu_iom1*zu_iom1 ) & ! add () for NP repro
+ & + ( zv_io*zv_io + zv_iom1*zv_iom1 ) ) ) * smask0(ji,jj)
zvel (ji,jj) = 0.5_wp * SQRT( ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) + &
& ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) )
END_2D
diff --git a/src/ICE/iceupdate.F90 b/src/ICE/iceupdate.F90
index 8efb5c2d0b2f2eb2c7c6d5e2f00afaa15264bcd5..8a0f233cb892a74bb56e877d883191ce5ab2f808 100644
--- a/src/ICE/iceupdate.F90
+++ b/src/ICE/iceupdate.F90
@@ -361,8 +361,8 @@ CONTAINS
IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN !== Ice time-step only ==! (i.e. surface module time-step)
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* rhoco * |U_ice-U_oce| at T-point
! ! 2*(U_ice-U_oce) at T-point
- zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj) ! u_oce = ssu_m
- zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1) ! v_oce = ssv_m
+ zu_t = ( u_ice(ji,jj) + u_ice(ji-1,jj) ) - ( u_oce(ji,jj) + u_oce(ji-1,jj) ) ! u_oce = ssu_m add () for
+ zv_t = ( v_ice(ji,jj) + v_ice(ji,jj-1) ) - ( v_oce(ji,jj) + v_oce(ji,jj-1) ) ! v_oce = ssv_m NP repro
! ! |U_ice-U_oce|^2
zmodt = 0.25_wp * ( zu_t * zu_t + zv_t * zv_t )
!
@@ -372,8 +372,8 @@ CONTAINS
!
DO_2D( 0, 0, 0, 0 ) !* save the air-ocean stresses at ice time-step
! ! 2*(U_ice-U_oce) at T-point
- zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj) ! u_oce = ssu_m
- zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1) ! v_oce = ssv_m
+ zu_t = ( u_ice(ji,jj) + u_ice(ji-1,jj) ) - ( u_oce(ji,jj) + u_oce(ji-1,jj) ) ! u_oce = ssu_m add () for
+ zv_t = ( v_ice(ji,jj) + v_ice(ji,jj-1) ) - ( v_oce(ji,jj) + v_oce(ji,jj-1) ) ! v_oce = ssv_m NP repro
! ! |U_ice-U_oce|^2
zmodt = 0.25_wp * ( zu_t * zu_t + zv_t * zv_t )
! ! update the ocean stress modulus
@@ -398,8 +398,10 @@ CONTAINS
!
DO_2D( 0, 0, 0, 0 ) !* update the stress WITHOUT an ice-ocean rotation angle
! ! linearized quadratic drag formulation
- zutau_ice = 0.5_wp * tmod_io(ji,jj) * ( u_ice(ji,jj) + u_ice(ji-1,jj) - pu_oce(ji,jj) - pu_oce(ji-1,jj) )
- zvtau_ice = 0.5_wp * tmod_io(ji,jj) * ( v_ice(ji,jj) + v_ice(ji,jj-1) - pv_oce(ji,jj) - pv_oce(ji,jj-1) )
+ zutau_ice = 0.5_wp * tmod_io(ji,jj) &
+ & * ( ( u_ice(ji,jj) + u_ice(ji-1,jj) ) - ( pu_oce(ji,jj) + pu_oce(ji-1,jj) ) ) ! add () for
+ zvtau_ice = 0.5_wp * tmod_io(ji,jj) &
+ & * ( ( v_ice(ji,jj) + v_ice(ji,jj-1) ) - ( pv_oce(ji,jj) + pv_oce(ji,jj-1) ) ) ! NP repro
! ! stresses at the ocean surface
utau(ji,jj) = ( 1._wp - at_i(ji,jj) ) * utau_oce(ji,jj) + at_i(ji,jj) * zutau_ice
vtau(ji,jj) = ( 1._wp - at_i(ji,jj) ) * vtau_oce(ji,jj) + at_i(ji,jj) * zvtau_ice
diff --git a/src/NST/agrif_oce_sponge.F90 b/src/NST/agrif_oce_sponge.F90
index 7d9e58852fc14f1635a9a49b88111f68d50b9d9c..2067c64a462c2c59a576842a56b22f26321b2265 100644
--- a/src/NST/agrif_oce_sponge.F90
+++ b/src/NST/agrif_oce_sponge.F90
@@ -239,8 +239,8 @@ CONTAINS
fspt(:,:) = 0._wp
fspf(:,:) = 0._wp
DO_2D( 0, 0, 0, 0 )
- fspt(ji,jj) = 0.25_wp * ( ztabramp(ji ,jj ) + ztabramp(ji-1,jj ) &
- & +ztabramp(ji ,jj-1) + ztabramp(ji-1,jj-1) ) * ssmask(ji,jj)
+ fspt(ji,jj) = 0.25_wp * ( ( ztabramp(ji ,jj ) + ztabramp(ji-1,jj ) ) & ! add () for NP repro
+ & + ( ztabramp(ji ,jj-1) + ztabramp(ji-1,jj-1) ) ) * ssmask(ji,jj)
fspf(ji,jj) = ztabramp(ji,jj) * ssvmask(ji,jj) * ssvmask(ji,jj+1)
END_2D
@@ -372,8 +372,8 @@ CONTAINS
fspt_2d(:,:) = 0._wp
fspf_2d(:,:) = 0._wp
DO_2D( 0, 0, 0, 0 )
- fspt_2d(ji,jj) = 0.25_wp * ( ztabramp(ji ,jj ) + ztabramp(ji-1,jj ) &
- & +ztabramp(ji ,jj-1) + ztabramp(ji-1,jj-1) ) * ssmask(ji,jj)
+ fspt_2d(ji,jj) = 0.25_wp * ( ( ztabramp(ji ,jj ) + ztabramp(ji-1,jj ) ) & ! add () for NP repro
+ & + ( ztabramp(ji ,jj-1) + ztabramp(ji-1,jj-1) ) ) * ssmask(ji,jj)
fspf_2d(ji,jj) = ztabramp(ji,jj) * ssvmask(ji,jj) * ssvmask(ji,jj+1)
END_2D
CALL lbc_lnk( 'agrif_Sponge_2d', fspu_2d, 'U', 1._wp, fspv_2d, 'V', 1._wp, fspt_2d, 'T', 1._wp, fspf_2d, 'F', 1._wp )
@@ -593,8 +593,8 @@ CONTAINS
IF (.NOT. tabspongedone_tsn(ji,jj)) THEN
zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a)
! horizontal diffusive trends
- ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) &
- & + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) &
+ ztsa = zbtr * ( ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) & ! add () for NP repro
+ & + ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) ) &
& - rn_trelax_tra * r1_Dt * fspt(ji,jj) * tsbdiff(ji,jj,jk,jn)
! add it to the general tracer trends
ts(ji,jj,jk,jn,Krhs_a) = ts(ji,jj,jk,jn,Krhs_a) + ztsa
diff --git a/src/NST/agrif_top_sponge.F90 b/src/NST/agrif_top_sponge.F90
index 641b9de74e96a4e352885b4c0cadddabd210ab70..00b159c615198b9580084372885d668d3e4c5c61 100644
--- a/src/NST/agrif_top_sponge.F90
+++ b/src/NST/agrif_top_sponge.F90
@@ -265,8 +265,8 @@ CONTAINS
IF (.NOT. tabspongedone_trn(ji,jj)) THEN
zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a)
! horizontal diffusive trends
- ztra = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) &
- & + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) &
+ ztra = zbtr * ( ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) & ! add () for NP repro
+ & + ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) ) &
& - rn_trelax_tra * r1_Dt * fspt(ji,jj) * trbdiff(ji,jj,jk,jn)
! add it to the general tracer trends
tr(ji,jj,jk,jn,Krhs_a) = tr(ji,jj,jk,jn,Krhs_a) + ztra
diff --git a/src/OCE/DOM/domqco.F90 b/src/OCE/DOM/domqco.F90
index 5944f59c1057780287cc63828673d090415ce239..678fd864db49e27b9bca9dc2a2a07b04b88a7a86 100644
--- a/src/OCE/DOM/domqco.F90
+++ b/src/OCE/DOM/domqco.F90
@@ -180,9 +180,9 @@ CONTAINS
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
pr3f(ji,jj) = 0.25_wp * ( ( e1e2t(ji ,jj ) * pssh(ji ,jj ) &
- & + e1e2t(ji+1,jj ) * pssh(ji+1,jj ) ) & ! bracket for halo 1 - halo 2 compatibility
+ & + e1e2t(ji+1,jj ) * pssh(ji+1,jj ) ) & ! add () for NP reproducibility
& + ( e1e2t(ji ,jj+1) * pssh(ji ,jj+1) &
- & + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1) ) & ! bracket for halo 1 - halo 2 compatibility
+ & + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1) ) & ! add () for NP reproducibility
& ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj)
END_2D
! ! lbc on ratio at u-,v-,f-points
@@ -247,9 +247,9 @@ CONTAINS
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
pr3f(ji,jj) = 0.25_wp * ( ( e1e2t(ji ,jj ) * pssh(ji ,jj ) &
- & + e1e2t(ji+1,jj ) * pssh(ji+1,jj ) ) & ! bracket for halo 1 - halo 2 compatibility
+ & + e1e2t(ji+1,jj ) * pssh(ji+1,jj ) ) & ! add () for NP reproducibility
& + ( e1e2t(ji ,jj+1) * pssh(ji ,jj+1) &
- & + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1) ) & ! bracket for halo 1 - halo 2 compatibility
+ & + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1) ) & ! add () for NP reproducibility
& ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj)
END_2D
!!st ELSE !- Flux Form (simple averaging)
@@ -258,7 +258,7 @@ CONTAINS
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
pr3f(ji,jj) = 0.25_wp * ( ( pssh(ji,jj ) + pssh(ji+1,jj ) ) &
- & + ( pssh(ji,jj+1) + pssh(ji+1,jj+1) ) & ! bracket for halo 1 - halo 2 compatibility
+ & + ( pssh(ji,jj+1) + pssh(ji+1,jj+1) ) & ! add () for NP reproducibility
& ) * r1_hf_0(ji,jj)
END_2D
!!st ENDIF
diff --git a/src/OCE/DOM/domvvl.F90 b/src/OCE/DOM/domvvl.F90
index 2de0762a7223f14face860df5c22b6729158807d..d0d585fc147558d6aca8f815b434d01a29042a29 100644
--- a/src/OCE/DOM/domvvl.F90
+++ b/src/OCE/DOM/domvvl.F90
@@ -422,9 +422,8 @@ CONTAINS
vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
END_2D
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! c - second derivative: divergence of diffusive fluxes
- tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) &
- & + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) &
- & ) * r1_e1e2t(ji,jj)
+ tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) ) & ! add () for NP repro
+ & + ( vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) ) ) * r1_e1e2t(ji,jj)
END_3D
! ! d - thickness diffusion transport: boundary conditions
! (stored for tracer advction and continuity equation)
diff --git a/src/OCE/DYN/divhor.F90 b/src/OCE/DYN/divhor.F90
index 795277715a0371243c155a465409880efb24177b..5e4188a2bd785fb6e81b295bd716fe9f7bcdce1e 100644
--- a/src/OCE/DYN/divhor.F90
+++ b/src/OCE/DYN/divhor.F90
@@ -83,11 +83,11 @@ CONTAINS
ENDIF
!
DO_3D_OVR( 0, 0, 0, 0, 1, jpkm1 )
- hdiv(ji,jj,jk) = ( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) * puu(ji ,jj,jk) &
- & - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * puu(ji-1,jj,jk) &
- & + e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) * pvv(ji,jj ,jk) &
- & - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * pvv(ji,jj-1,jk) ) &
- & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+ hdiv(ji,jj,jk) = ( ( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) * puu(ji ,jj,jk) & ! add () for NP repro
+ & - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * puu(ji-1,jj,jk) ) &
+ & + ( e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) * pvv(ji,jj ,jk) &
+ & - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * pvv(ji,jj-1,jk) ) &
+ & ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
END_3D
!
IF( ln_rnf ) CALL sbc_rnf_div( hdiv, Kmm ) !== + runoffs divergence ==!
@@ -149,15 +149,11 @@ CONTAINS
ENDIF
!
DO_3D_OVR( 0, 0, 0, 0, 1, jpkm1 ) !== Horizontal divergence ==!
- ! round brackets added to fix the order of floating point operations
- ! needed to ensure halo 1 - halo 2 compatibility
- hdiv(ji,jj,jk) = ( ( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) * uu(ji ,jj,jk,Kmm) &
- & - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & + ( e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) * vv(ji,jj ,jk,Kmm) &
- & - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+ hdiv(ji,jj,jk) = ( ( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) * uu(ji ,jj,jk,Kmm) & ! add () for NP repro
+ & - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm) ) &
+ & + ( e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) * vv(ji,jj ,jk,Kmm) &
+ & - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm) ) &
+ & ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
END_3D
!
IF( ln_rnf ) CALL sbc_rnf_div( hdiv, Kmm ) !== runoffs ==! (update hdiv field)
diff --git a/src/OCE/DYN/dynadv_cen2.F90 b/src/OCE/DYN/dynadv_cen2.F90
index 07dd931129ed0b4474092915e3b6c83abcaf2ff5..38306ccd6a2ebd83c8aa04249852327a36235387 100644
--- a/src/OCE/DYN/dynadv_cen2.F90
+++ b/src/OCE/DYN/dynadv_cen2.F90
@@ -103,11 +103,11 @@ CONTAINS
zfv_t(ji ,jj+1) = ( zfv(ji,jj) + zfv(ji,jj+1) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji ,jj+1,jk,Kmm) )
END_2D
DO_2D( 0, 0, 0, 0 ) ! divergence of horizontal momentum fluxes
- puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_t(ji+1,jj) - zfu_t(ji,jj ) &
- & + zfv_f(ji ,jj) - zfv_f(ji,jj-1) ) * r1_e1e2u(ji,jj) &
+ puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( ( zfu_t(ji+1,jj) - zfu_t(ji,jj ) ) & ! add () for NP repro
+ & + ( zfv_f(ji ,jj) - zfv_f(ji,jj-1) ) ) * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,jk,Kmm)
- pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfu_f(ji,jj ) - zfu_f(ji-1,jj) &
- & + zfv_t(ji,jj+1) - zfv_t(ji ,jj) ) * r1_e1e2v(ji,jj) &
+ pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( ( zfu_f(ji,jj ) - zfu_f(ji-1,jj) ) & ! add () for NP repro
+ & + ( zfv_t(ji,jj+1) - zfv_t(ji ,jj) ) ) * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,jk,Kmm)
END_2D
END DO
diff --git a/src/OCE/DYN/dynadv_ubs.F90 b/src/OCE/DYN/dynadv_ubs.F90
index 8ab908124d5b85beb120d5726d2d6ec12b8d28f7..f4f08f75b0bc3ce6a92b10e5df7093c51eefba0a 100644
--- a/src/OCE/DYN/dynadv_ubs.F90
+++ b/src/OCE/DYN/dynadv_ubs.F90
@@ -123,37 +123,23 @@ CONTAINS
END_2D
!
DO_2D( 1, 1, 1, 1 ) ! laplacian
- ! round brackets added to fix the order of floating point operations
- ! needed to ensure halo 1 - halo 2 compatibility (north fold)
-!! zlu_uu(ji,jj,1) = ( ( puu(ji+1,jj ,jk,Kbb) + puu(ji-1,jj ,jk,Kbb) ) - 2._wp * puu(ji,jj,jk,Kbb) ) * umask(ji,jj,jk)
-!! zlv_vv(ji,jj,1) = ( ( pvv(ji ,jj+1,jk,Kbb) + pvv(ji ,jj-1,jk,Kbb) ) - 2._wp * pvv(ji,jj,jk,Kbb) ) * vmask(ji,jj,jk)
- zlu_uu(ji,jj,1) = ( ( puu (ji+1,jj ,jk,Kbb) - puu (ji ,jj ,jk,Kbb) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & + ( puu (ji-1,jj ,jk,Kbb) - puu (ji ,jj ,jk,Kbb) &
- & ) ) * umask(ji ,jj ,jk)
- zlv_vv(ji,jj,1) = ( ( pvv (ji ,jj+1,jk,Kbb) - pvv (ji ,jj ,jk,Kbb) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & + ( pvv (ji ,jj-1,jk,Kbb) - pvv (ji ,jj ,jk,Kbb) &
- & ) ) * vmask(ji ,jj ,jk)
+ zlu_uu(ji,jj,1) = ( ( puu (ji+1,jj ,jk,Kbb) - puu (ji ,jj ,jk,Kbb) ) & ! add () for NP repro
+ & + ( puu (ji-1,jj ,jk,Kbb) - puu (ji ,jj ,jk,Kbb) ) ) * umask(ji ,jj ,jk)
+ zlv_vv(ji,jj,1) = ( ( pvv (ji ,jj+1,jk,Kbb) - pvv (ji ,jj ,jk,Kbb) ) & ! add () for NP repro
+ & + ( pvv (ji ,jj-1,jk,Kbb) - pvv (ji ,jj ,jk,Kbb) ) ) * vmask(ji ,jj ,jk)
zlu_uv(ji,jj,1) = ( puu(ji ,jj+1,jk,Kbb) - puu(ji ,jj ,jk,Kbb) ) * fmask(ji ,jj ,jk) &
& - ( puu(ji ,jj ,jk,Kbb) - puu(ji ,jj-1,jk,Kbb) ) * fmask(ji ,jj-1,jk)
zlv_vu(ji,jj,1) = ( pvv(ji+1,jj ,jk,Kbb) - pvv(ji ,jj ,jk,Kbb) ) * fmask(ji ,jj ,jk) &
& - ( pvv(ji ,jj ,jk,Kbb) - pvv(ji-1,jj ,jk,Kbb) ) * fmask(ji-1,jj ,jk)
!
-!! zlu_uu(ji,jj,2) = ( ( zfu(ji+1,jj ) + zfu(ji-1,jj ) ) - 2._wp * zfu(ji,jj) ) * umask(ji ,jj ,jk)
-!! zlv_vv(ji,jj,2) = ( ( zfv(ji ,jj+1) + zfv(ji ,jj-1) ) - 2._wp * zfv(ji,jj) ) * vmask(ji ,jj ,jk)
- zlu_uu(ji,jj,2) = ( ( zfu(ji+1,jj ) - zfu(ji ,jj ) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & + ( zfu(ji-1,jj ) - zfu(ji ,jj ) &
- & ) ) * umask(ji ,jj ,jk)
- zlv_vv(ji,jj,2) = ( ( zfv(ji ,jj+1) - zfv(ji ,jj ) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & + ( zfv(ji ,jj-1) - zfv(ji ,jj ) &
- & ) ) * vmask(ji ,jj ,jk)
- zlu_uv(ji,jj,2) = ( zfu(ji ,jj+1) - zfu(ji ,jj ) ) * fmask(ji ,jj ,jk) &
- & - ( zfu(ji ,jj ) - zfu(ji ,jj-1) ) * fmask(ji ,jj-1,jk)
- zlv_vu(ji,jj,2) = ( zfv(ji+1,jj ) - zfv(ji ,jj ) ) * fmask(ji ,jj ,jk) &
- & - ( zfv(ji ,jj ) - zfv(ji-1,jj ) ) * fmask(ji-1,jj ,jk)
+ zlu_uu(ji,jj,2) = ( ( zfu(ji+1,jj ) - zfu(ji ,jj ) ) & ! add () for NP repro
+ & + ( zfu(ji-1,jj ) - zfu(ji ,jj ) ) ) * umask(ji ,jj ,jk)
+ zlv_vv(ji,jj,2) = ( ( zfv(ji ,jj+1) - zfv(ji ,jj ) ) & ! add () for NP repro
+ & + ( zfv(ji ,jj-1) - zfv(ji ,jj ) ) ) * vmask(ji ,jj ,jk)
+ zlu_uv(ji,jj,2) = ( zfu(ji ,jj+1) - zfu(ji ,jj ) ) * fmask(ji ,jj ,jk) &
+ & - ( zfu(ji ,jj ) - zfu(ji ,jj-1) ) * fmask(ji ,jj-1,jk)
+ zlv_vu(ji,jj,2) = ( zfv(ji+1,jj ) - zfv(ji ,jj ) ) * fmask(ji ,jj ,jk) &
+ & - ( zfv(ji ,jj ) - zfv(ji-1,jj ) ) * fmask(ji-1,jj ,jk)
END_2D
!
! ! ====================== !
@@ -196,11 +182,11 @@ CONTAINS
& * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) - gamma1 * zl_v )
END_2D
DO_2D( 0, 0, 0, 0 ) ! divergence of horizontal momentum fluxes
- puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_t(ji+1,jj) - zfu_t(ji,jj ) &
- & + zfv_f(ji ,jj) - zfv_f(ji,jj-1) ) * r1_e1e2u(ji,jj) &
+ puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( ( zfu_t(ji+1,jj) - zfu_t(ji,jj ) ) & ! add () for NP repro
+ & + ( zfv_f(ji ,jj) - zfv_f(ji,jj-1) ) ) * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,jk,Kmm)
- pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfu_f(ji,jj ) - zfu_f(ji-1,jj) &
- & + zfv_t(ji,jj+1) - zfv_t(ji ,jj) ) * r1_e1e2v(ji,jj) &
+ pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( ( zfu_f(ji,jj ) - zfu_f(ji-1,jj) ) & ! add () for NP repro
+ & + ( zfv_t(ji,jj+1) - zfv_t(ji ,jj) ) ) * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,jk,Kmm)
END_2D
END DO
@@ -244,10 +230,8 @@ CONTAINS
zzfu_kp1 = ( zfw(ji,jj) + zfw(ji+1,jj ) ) * ( puu(ji,jj,jk+1,Kmm) + puu(ji,jj,jk,Kmm) )
zzfv_kp1 = ( zfw(ji,jj) + zfw(ji ,jj+1) ) * ( pvv(ji,jj,jk+1,Kmm) + pvv(ji,jj,jk,Kmm) )
! ! divergence of vertical momentum flux
- puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj) - zzfu_kp1 ) * r1_e1e2u(ji,jj) &
- & / e3u(ji,jj,jk,Kmm)
- pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj) - zzfv_kp1 ) * r1_e1e2v(ji,jj) &
- & / e3v(ji,jj,jk,Kmm)
+ puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj) - zzfu_kp1 ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
+ pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj) - zzfv_kp1 ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
! ! store vertical flux for next level calculation
zfu_uw(ji,jj) = zzfu_kp1
zfv_vw(ji,jj) = zzfv_kp1
@@ -256,10 +240,8 @@ CONTAINS
!
jk = jpkm1 ! compute last level (zzfu_kp1 = 0)
DO_2D( 0, 0, 0, 0 )
- puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - zfu_uw(ji,jj) * r1_e1e2u(ji,jj) &
- & / e3u(ji,jj,jk,Kmm)
- pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - zfv_vw(ji,jj) * r1_e1e2v(ji,jj) &
- & / e3v(ji,jj,jk,Kmm)
+ puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - zfu_uw(ji,jj) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
+ pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - zfv_vw(ji,jj) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
END_2D
!
IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic
@@ -366,45 +348,31 @@ CONTAINS
END_2D
!
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! laplacian
- ! round brackets added to fix the order of floating point operations
- ! needed to ensure halo 1 - halo 2 compatibility
- zlu_uu(ji,jj,jk,1) = ( ( puu (ji+1,jj ,jk,Kbb) - puu (ji ,jj ,jk,Kbb) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & + ( puu (ji-1,jj ,jk,Kbb) - puu (ji ,jj ,jk,Kbb) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & ) * umask(ji ,jj ,jk)
- zlv_vv(ji,jj,jk,1) = ( ( pvv (ji ,jj+1,jk,Kbb) - pvv (ji ,jj ,jk,Kbb) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & + ( pvv (ji ,jj-1,jk,Kbb) - pvv (ji ,jj ,jk,Kbb) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & ) * vmask(ji ,jj ,jk)
- zlu_uv(ji,jj,jk,1) = ( puu (ji ,jj+1,jk,Kbb) - puu (ji ,jj ,jk,Kbb) ) * fmask(ji ,jj ,jk) &
- & - ( puu (ji ,jj ,jk,Kbb) - puu (ji ,jj-1,jk,Kbb) ) * fmask(ji ,jj-1,jk)
- zlv_vu(ji,jj,jk,1) = ( pvv (ji+1,jj ,jk,Kbb) - pvv (ji ,jj ,jk,Kbb) ) * fmask(ji ,jj ,jk) &
- & - ( pvv (ji ,jj ,jk,Kbb) - pvv (ji-1,jj ,jk,Kbb) ) * fmask(ji-1,jj ,jk)
+ zlu_uu(ji,jj,jk,1) = ( ( puu (ji+1,jj ,jk,Kbb) - puu (ji ,jj ,jk,Kbb) ) & ! add () for NP repro
+ & + ( puu (ji-1,jj ,jk,Kbb) - puu (ji ,jj ,jk,Kbb) ) ) * umask(ji ,jj ,jk)
+ zlv_vv(ji,jj,jk,1) = ( ( pvv (ji ,jj+1,jk,Kbb) - pvv (ji ,jj ,jk,Kbb) ) & ! add () for NP repro
+ & + ( pvv (ji ,jj-1,jk,Kbb) - pvv (ji ,jj ,jk,Kbb) ) ) * vmask(ji ,jj ,jk)
+ zlu_uv(ji,jj,jk,1) = ( puu (ji ,jj+1,jk,Kbb) - puu (ji ,jj ,jk,Kbb) ) * fmask(ji ,jj ,jk) &
+ & - ( puu (ji ,jj ,jk,Kbb) - puu (ji ,jj-1,jk,Kbb) ) * fmask(ji ,jj-1,jk)
+ zlv_vu(ji,jj,jk,1) = ( pvv (ji+1,jj ,jk,Kbb) - pvv (ji ,jj ,jk,Kbb) ) * fmask(ji ,jj ,jk) &
+ & - ( pvv (ji ,jj ,jk,Kbb) - pvv (ji-1,jj ,jk,Kbb) ) * fmask(ji-1,jj ,jk)
!
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
- zlu_uu(ji,jj,jk,2) = ( ( zfu(ji+1,jj ,jk) - zfu(ji ,jj ,jk) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & + ( zfu(ji-1,jj ,jk) - zfu(ji ,jj ,jk) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & ) * umask(ji ,jj ,jk)
- zlv_vv(ji,jj,jk,2) = ( ( zfv(ji ,jj+1,jk) - zfv(ji ,jj ,jk) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & + ( zfv(ji ,jj-1,jk) - zfv(ji ,jj ,jk) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & ) * vmask(ji ,jj ,jk)
- zlu_uv(ji,jj,jk,2) = ( zfu(ji ,jj+1,jk) - zfu(ji ,jj ,jk) ) * fmask(ji ,jj ,jk) &
- & - ( zfu(ji ,jj ,jk) - zfu(ji ,jj-1,jk) ) * fmask(ji ,jj-1,jk)
- zlv_vu(ji,jj,jk,2) = ( zfv(ji+1,jj ,jk) - zfv(ji ,jj ,jk) ) * fmask(ji ,jj ,jk) &
- & - ( zfv(ji ,jj ,jk) - zfv(ji-1,jj ,jk) ) * fmask(ji-1,jj ,jk)
+ zlu_uu(ji,jj,jk,2) = ( ( zfu(ji+1,jj ,jk) - zfu(ji ,jj ,jk) ) & ! add () for NP repro
+ & + ( zfu(ji-1,jj ,jk) - zfu(ji ,jj ,jk) ) ) * umask(ji ,jj ,jk)
+ zlv_vv(ji,jj,jk,2) = ( ( zfv(ji ,jj+1,jk) - zfv(ji ,jj ,jk) ) & ! add () for NP repro
+ & + ( zfv(ji ,jj-1,jk) - zfv(ji ,jj ,jk) ) ) * vmask(ji ,jj ,jk)
+ zlu_uv(ji,jj,jk,2) = ( zfu(ji ,jj+1,jk) - zfu(ji ,jj ,jk) ) * fmask(ji ,jj ,jk) &
+ & - ( zfu(ji ,jj ,jk) - zfu(ji ,jj-1,jk) ) * fmask(ji ,jj-1,jk)
+ zlv_vu(ji,jj,jk,2) = ( zfv(ji+1,jj ,jk) - zfv(ji ,jj ,jk) ) * fmask(ji ,jj ,jk) &
+ & - ( zfv(ji ,jj ,jk) - zfv(ji-1,jj ,jk) ) * fmask(ji-1,jj ,jk)
END_2D
END DO
IF( nn_hls == 1 ) CALL lbc_lnk( 'dynadv_ubs', zlu_uu(:,:,:,1), 'U', -1.0_wp , zlu_uv(:,:,:,1), 'U', -1.0_wp, &
- & zlu_uu(:,:,:,2), 'U', -1.0_wp , zlu_uv(:,:,:,2), 'U', -1.0_wp, &
- & zlv_vv(:,:,:,1), 'V', -1.0_wp , zlv_vu(:,:,:,1), 'V', -1.0_wp, &
- & zlv_vv(:,:,:,2), 'V', -1.0_wp , zlv_vu(:,:,:,2), 'V', -1.0_wp )
+ & zlu_uu(:,:,:,2), 'U', -1.0_wp , zlu_uv(:,:,:,2), 'U', -1.0_wp, &
+ & zlv_vv(:,:,:,1), 'V', -1.0_wp , zlv_vu(:,:,:,1), 'V', -1.0_wp, &
+ & zlv_vv(:,:,:,2), 'V', -1.0_wp , zlv_vu(:,:,:,2), 'V', -1.0_wp )
!
! ! ====================== !
! ! Horizontal advection !
@@ -426,16 +394,14 @@ CONTAINS
ELSE ; zl_v = zlv_vv(ji,jj+1,jk,1)
ENDIF
!
- zfu_t(ji+1,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj ,jk) &
- & - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj ,jk,2) ) ) &
- & * ( zui - gamma1 * zl_u)
- zfv_t(ji ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji ,jj+1,jk) &
- & - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji ,jj+1,jk,2) ) ) &
- & * ( zvj - gamma1 * zl_v)
+ zfu_t(ji+1,jj ,jk) = ( ( zfu(ji,jj,jk ) + zfu(ji+1,jj ,jk ) ) & ! add () for NP repro
+ & - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj ,jk,2) ) ) * ( zui - gamma1 * zl_u)
+ zfv_t(ji ,jj+1,jk) = ( ( zfv(ji,jj,jk ) + zfv(ji ,jj+1,jk ) ) & ! add () for NP repro
+ & - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji ,jj+1,jk,2) ) ) * ( zvj - gamma1 * zl_v)
!
zfuj = ( zfu(ji,jj,jk) + zfu(ji ,jj+1,jk) )
zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj ,jk) )
- IF( zfuj > 0 ) THEN ; zl_v = zlv_vu( ji ,jj ,jk,1)
+ IF( zfuj > 0 ) THEN ; zl_v = zlv_vu( ji ,jj,jk,1)
ELSE ; zl_v = zlv_vu( ji+1,jj,jk,1)
ENDIF
IF( zfvi > 0 ) THEN ; zl_u = zlu_uv( ji,jj ,jk,1)
@@ -443,17 +409,17 @@ CONTAINS
ENDIF
!
zfv_f(ji ,jj ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj ,jk,2) ) ) &
- & * ( puu(ji,jj,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) - gamma1 * zl_u )
+ & * ( ( puu(ji,jj,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) ) - gamma1 * zl_u ) ! add () for NP repro
zfu_f(ji ,jj ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji ,jj+1,jk,2) ) ) &
- & * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) - gamma1 * zl_v )
+ & * ( ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) ) - gamma1 * zl_v ) ! add () for NP repro
END_2D
DO_2D( 0, 0, 0, 0 ) ! divergence of horizontal momentum fluxes
- puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_t(ji+1,jj,jk) - zfu_t(ji,jj ,jk) &
- & + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) &
- & / e3u(ji,jj,jk,Kmm)
- pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfu_f(ji,jj ,jk) - zfu_f(ji-1,jj,jk) &
- & + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) &
- & / e3v(ji,jj,jk,Kmm)
+ puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( ( zfu_t(ji+1,jj,jk) - zfu_t(ji,jj ,jk) ) & ! add () for NP repro
+ & + ( zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) ) * r1_e1e2u(ji,jj) &
+ & / e3u(ji,jj,jk,Kmm)
+ pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( ( zfu_f(ji,jj ,jk) - zfu_f(ji-1,jj,jk) ) & ! add () for NP repro
+ & + ( zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) ) * r1_e1e2v(ji,jj) &
+ & / e3v(ji,jj,jk,Kmm)
END_2D
END DO
IF( l_trddyn ) THEN ! trends: send trends to trddyn for diagnostic
@@ -484,15 +450,13 @@ CONTAINS
zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * zpt_w(ji,jj,jk)
END_2D
DO_2D( 0, 0, 0, 0 )
- zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji,jj,jk-1,Kmm) )
- zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk-1,Kmm) )
+ zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji,jj,jk-1,Kmm) )
+ zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk-1,Kmm) )
END_2D
END DO
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! divergence of vertical momentum flux divergence
- puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) &
- & / e3u(ji,jj,jk,Kmm)
- pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) &
- & / e3v(ji,jj,jk,Kmm)
+ puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
+ pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
END_3D
!
IF( l_trddyn ) THEN ! save the vertical advection trend for diagnostic
diff --git a/src/OCE/DYN/dynhpg.F90 b/src/OCE/DYN/dynhpg.F90
index 145c7f9e17407167fc54774c8fe282f407e31f8b..6cdd91ccd2919700cc6f4f06061130f19ab72d8c 100644
--- a/src/OCE/DYN/dynhpg.F90
+++ b/src/OCE/DYN/dynhpg.F90
@@ -369,13 +369,13 @@ CONTAINS
IF( iku > 1 ) THEN ! on i-direction (level 2 or more)
puu (ji,jj,iku,Krhs) = puu(ji,jj,iku,Krhs) - zhpi(ji,jj,iku) ! subtract old value
zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) & ! compute the new one
- & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + zgru(ji,jj) ) * r1_e1u(ji,jj)
+ & + zcoef2 * ( ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) ) + zgru(ji,jj) ) * r1_e1u(ji,jj)
puu (ji,jj,iku,Krhs) = puu(ji,jj,iku,Krhs) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend
ENDIF
IF( ikv > 1 ) THEN ! on j-direction (level 2 or more)
pvv (ji,jj,ikv,Krhs) = pvv(ji,jj,ikv,Krhs) - zhpj(ji,jj,ikv) ! subtract old value
zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) & ! compute the new one
- & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + zgrv(ji,jj) ) * r1_e2v(ji,jj)
+ & + zcoef3 * ( ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) ) + zgrv(ji,jj) ) * r1_e2v(ji,jj)
pvv (ji,jj,ikv,Krhs) = pvv(ji,jj,ikv,Krhs) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend
ENDIF
END_2D
@@ -456,8 +456,8 @@ CONTAINS
zcpy(ji,jj) = 1.0_wp
ELSE IF(ll_tmp2) THEN
! no worries about ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm) = 0, it won't happen ! here
- zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
- & / (ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm)) )
+ zcpy(ji,jj) = ABS( ( ( ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) - ( ssh(ji,jj,Kmm) + ht_0(ji,jj) ) ) & ! add () for NP repro
+ & / ( ssh(ji,jj+1,Kmm) - ssh(ji,jj,Kmm) ) )
ELSE
zcpy(ji,jj) = 0._wp
END IF
@@ -473,9 +473,9 @@ CONTAINS
& * ( e3w(ji ,jj+1,1,Kmm) * rhd(ji ,jj+1,1) &
& - e3w(ji ,jj ,1,Kmm) * rhd(ji ,jj ,1) )
! ! s-coordinate pressure gradient correction
- zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) ) &
+ zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) ) &
& * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj)
- zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) ) &
+ zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) ) &
& * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj)
!
IF( ln_wd_il ) THEN
@@ -663,8 +663,8 @@ CONTAINS
zcpx(ji,jj) = 1.0_wp
ELSE IF(ll_tmp2) THEN
! no worries about ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm) = 0, it won't happen ! here
- zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
- & / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) )
+ zcpx(ji,jj) = ABS( ( ( ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) ) - ( ssh(ji,jj,Kmm) + ht_0(ji,jj) ) ) & ! add () for NP repro
+ & / ( ssh(ji+1,jj,Kmm) - ssh(ji,jj,Kmm) ) )
ELSE
zcpx(ji,jj) = 0._wp
END IF
@@ -681,8 +681,8 @@ CONTAINS
zcpy(ji,jj) = 1.0_wp
ELSE IF(ll_tmp2) THEN
! no worries about ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm) = 0, it won't happen ! here
- zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
- & / (ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm)) )
+ zcpy(ji,jj) = ABS( ( ( ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) - ( ssh(ji,jj,Kmm) + ht_0(ji,jj) ) ) & ! add () for NP repro
+ & / ( ssh(ji,jj+1,Kmm) - ssh(ji,jj,Kmm) ) )
ELSE
zcpy(ji,jj) = 0._wp
END IF
@@ -709,7 +709,7 @@ CONTAINS
!!bug gm Not a true bug, but... zdzz=e3w for zdzx, zdzy verify what it is really
DO_3D( 1, 1, 1, 1, 2, jpkm1 )
- zdrhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd (ji,jj,jk-1)
+ zdrhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd (ji,jj,jk-1)
zdzz (ji,jj,jk) = - gde3w(ji ,jj ,jk) + gde3w(ji,jj,jk-1)
END_3D
@@ -727,7 +727,7 @@ CONTAINS
z1_cff = zdrhoz(ji,jj,jk) + zdrhoz(ji,jj,jk+1)
zdrho_k(ji,jj,jk) = cffw / SIGN( MAX( ABS(z1_cff), zep ), z1_cff )
zdz_k(ji,jj,jk) = 2._wp * zdzz(ji,jj,jk) * zdzz(ji,jj,jk+1) &
- & / ( zdzz(ji,jj,jk) + zdzz(ji,jj,jk+1) )
+ & / ( zdzz(ji,jj,jk) + zdzz(ji,jj,jk+1) )
END_3D
!----------------------------------------------------------------------------------
@@ -743,7 +743,7 @@ CONTAINS
DO_2D( 1, 1, 1, 1 )
IF ( mbkt(ji,jj)>1 ) THEN
iktb = mbkt(ji,jj)
- zdrho_k(ji,jj,iktb) = aco_bc_vrt * ( rhd(ji,jj,iktb) - rhd(ji,jj,iktb-1) ) - bco_bc_vrt * zdrho_k(ji,jj,iktb-1)
+ zdrho_k(ji,jj,iktb) = aco_bc_vrt * ( rhd(ji,jj,iktb) - rhd(ji,jj,iktb-1) ) - bco_bc_vrt * zdrho_k(ji,jj,iktb-1)
zdz_k (ji,jj,iktb) = aco_bc_vrt * (-gde3w(ji,jj,iktb) + gde3w(ji,jj,iktb-1) ) - bco_bc_vrt * zdz_k (ji,jj,iktb-1)
END IF
END_2D
@@ -772,14 +772,14 @@ CONTAINS
!-------------------------------------------------------------
DO_3D( 0, 1, 0, 1, 2, jpkm1 )
- z_rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd (ji,jj,jk-1) ) &
+ z_rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd (ji,jj,jk-1) ) &
& * ( - gde3w(ji,jj,jk) + gde3w(ji,jj,jk-1) ) &
& + z_grav_10 * ( &
- & ( zdrho_k (ji,jj,jk) - zdrho_k (ji,jj,jk-1) ) &
- & * ( - gde3w(ji,jj,jk) + gde3w(ji,jj,jk-1) - z1_12 * ( zdz_k (ji,jj,jk) + zdz_k (ji,jj,jk-1) ) ) &
- & - ( zdz_k (ji,jj,jk) - zdz_k (ji,jj,jk-1) ) &
- & * ( rhd (ji,jj,jk) - rhd (ji,jj,jk-1) - z1_12 * ( zdrho_k(ji,jj,jk) + zdrho_k(ji,jj,jk-1) ) ) &
- & )
+ & ( zdrho_k(ji,jj,jk) - zdrho_k(ji,jj,jk-1) ) &
+ & * ( - ( gde3w( ji,jj,jk) - gde3w( ji,jj,jk-1) ) - z1_12 * ( zdz_k (ji,jj,jk) + zdz_k (ji,jj,jk-1) ) ) &
+ & - ( zdz_k( ji,jj,jk) - zdz_k( ji,jj,jk-1) ) &
+ & * ( ( rhd( ji,jj,jk) - rhd( ji,jj,jk-1) ) - z1_12 * ( zdrho_k(ji,jj,jk) + zdrho_k(ji,jj,jk-1) ) ) &
+ & )
END_3D
!----------------------------------------------------------------------------------------
@@ -835,28 +835,28 @@ CONTAINS
! Walls coming from left: should check from 2 to jpi-1 (and jpj=2-jpj)
DO_2D( 0, 0, 0, 1 )
IF ( umask(ji,jj,jk) > 0.5_wp .AND. umask(ji-1,jj,jk) < 0.5_wp .AND. umask(ji+1,jj,jk) > 0.5_wp) THEN
- zz_drho_i(ji,jj) = aco_bc_hor * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) ) - bco_bc_hor * zdrho_i(ji+1,jj,jk)
+ zz_drho_i(ji,jj) = aco_bc_hor * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) ) - bco_bc_hor * zdrho_i(ji+1,jj,jk)
zz_dz_i (ji,jj) = aco_bc_hor * (-gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_i (ji+1,jj,jk)
END IF
END_2D
! Walls coming from right: should check from 3 to jpi (and jpj=2-jpj)
DO_2D( -1, 1, 0, 1 )
IF ( umask(ji,jj,jk) < 0.5_wp .AND. umask(ji-1,jj,jk) > 0.5_wp .AND. umask(ji-2,jj,jk) > 0.5_wp) THEN
- zz_drho_i(ji,jj) = aco_bc_hor * ( rhd (ji,jj,jk) - rhd (ji-1,jj,jk) ) - bco_bc_hor * zdrho_i(ji-1,jj,jk)
+ zz_drho_i(ji,jj) = aco_bc_hor * ( rhd (ji,jj,jk) - rhd (ji-1,jj,jk) ) - bco_bc_hor * zdrho_i(ji-1,jj,jk)
zz_dz_i (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji-1,jj,jk) ) - bco_bc_hor * zdz_i (ji-1,jj,jk)
END IF
END_2D
! Walls coming from left: should check from 2 to jpj-1 (and jpi=2-jpi)
DO_2D( 0, 1, 0, 0 )
IF ( vmask(ji,jj,jk) > 0.5_wp .AND. vmask(ji,jj-1,jk) < 0.5_wp .AND. vmask(ji,jj+1,jk) > 0.5_wp) THEN
- zz_drho_j(ji,jj) = aco_bc_hor * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) ) - bco_bc_hor * zdrho_j(ji,jj+1,jk)
+ zz_drho_j(ji,jj) = aco_bc_hor * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) ) - bco_bc_hor * zdrho_j(ji,jj+1,jk)
zz_dz_j (ji,jj) = aco_bc_hor * (-gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_j (ji,jj+1,jk)
END IF
END_2D
! Walls coming from right: should check from 3 to jpj (and jpi=2-jpi)
DO_2D( 0, 1, -1, 1 )
IF ( vmask(ji,jj,jk) < 0.5_wp .AND. vmask(ji,jj-1,jk) > 0.5_wp .AND. vmask(ji,jj-2,jk) > 0.5_wp) THEN
- zz_drho_j(ji,jj) = aco_bc_hor * ( rhd (ji,jj,jk) - rhd (ji,jj-1,jk) ) - bco_bc_hor * zdrho_j(ji,jj-1,jk)
+ zz_drho_j(ji,jj) = aco_bc_hor * ( rhd (ji,jj,jk) - rhd (ji,jj-1,jk) ) - bco_bc_hor * zdrho_j(ji,jj-1,jk)
zz_dz_j (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji,jj-1,jk) ) - bco_bc_hor * zdz_j (ji,jj-1,jk)
END IF
END_2D
@@ -873,25 +873,25 @@ CONTAINS
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
! two -ve signs cancel in next two lines (within zcoef0 and because gde3w is a depth not a height)
z_rho_i(ji,jj,jk) = zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) &
- & * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) )
+ & * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) )
IF ( umask(ji-1, jj, jk) > 0.5 .OR. umask(ji+1, jj, jk) > 0.5 ) THEN
z_rho_i(ji,jj,jk) = z_rho_i(ji,jj,jk) - z_grav_10 * ( &
- & ( zdrho_i (ji+1,jj,jk) - zdrho_i (ji,jj,jk) ) &
- & * ( - gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) - z1_12 * ( zdz_i (ji+1,jj,jk) + zdz_i (ji,jj,jk) ) ) &
- & - ( zdz_i (ji+1,jj,jk) - zdz_i (ji,jj,jk) ) &
- & * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) - z1_12 * ( zdrho_i(ji+1,jj,jk) + zdrho_i(ji,jj,jk) ) ) &
- & )
+ & ( zdrho_i(ji+1,jj,jk) - zdrho_i(ji,jj,jk) ) &
+ & * ( - ( gde3w( ji+1,jj,jk) - gde3w( ji,jj,jk) ) - z1_12 * ( zdz_i (ji+1,jj,jk) + zdz_i (ji,jj,jk) ) ) &
+ & - ( zdz_i( ji+1,jj,jk) - zdz_i( ji,jj,jk) ) &
+ & * ( ( rhd( ji+1,jj,jk) - rhd( ji,jj,jk) ) - z1_12 * ( zdrho_i(ji+1,jj,jk) + zdrho_i(ji,jj,jk) ) ) &
+ & )
END IF
z_rho_j(ji,jj,jk) = zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) &
& * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) )
IF ( vmask(ji, jj-1, jk) > 0.5 .OR. vmask(ji, jj+1, jk) > 0.5 ) THEN
z_rho_j(ji,jj,jk) = z_rho_j(ji,jj,jk) - z_grav_10 * ( &
- & ( zdrho_j (ji,jj+1,jk) - zdrho_j (ji,jj,jk) ) &
- & * ( - gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) - z1_12 * ( zdz_j (ji,jj+1,jk) + zdz_j (ji,jj,jk) ) ) &
- & - ( zdz_j (ji,jj+1,jk) - zdz_j (ji,jj,jk) ) &
- & * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) - z1_12 * ( zdrho_j(ji,jj+1,jk) + zdrho_j(ji,jj,jk) ) ) &
- & )
+ & ( zdrho_j(ji,jj+1,jk) - zdrho_j(ji,jj,jk) ) &
+ & * ( - ( gde3w( ji,jj+1,jk) - gde3w( ji,jj,jk) ) - z1_12 * ( zdz_j (ji,jj+1,jk) + zdz_j (ji,jj,jk) ) ) &
+ & - ( zdz_j( ji,jj+1,jk) - zdz_j( ji,jj,jk) ) &
+ & * ( ( rhd( ji,jj+1,jk) - rhd( ji,jj,jk) ) - z1_12 * ( zdrho_j(ji,jj+1,jk) + zdrho_j(ji,jj,jk) ) ) &
+ & )
END IF
END_3D
@@ -903,8 +903,8 @@ CONTAINS
! Surface value
! ---------------
DO_2D( 0, 0, 0, 0 )
- zhpi(ji,jj,1) = ( z_rho_k(ji,jj,1) - z_rho_k(ji+1,jj ,1) - z_rho_i(ji,jj,1) ) * r1_e1u(ji,jj)
- zhpj(ji,jj,1) = ( z_rho_k(ji,jj,1) - z_rho_k(ji ,jj+1,1) - z_rho_j(ji,jj,1) ) * r1_e2v(ji,jj)
+ zhpi(ji,jj,1) = ( ( z_rho_k(ji,jj,1) - z_rho_k(ji+1,jj ,1) ) - z_rho_i(ji,jj,1) ) * r1_e1u(ji,jj) ! add () for NP repro
+ zhpj(ji,jj,1) = ( ( z_rho_k(ji,jj,1) - z_rho_k(ji ,jj+1,1) ) - z_rho_j(ji,jj,1) ) * r1_e2v(ji,jj)
IF( ln_wd_il ) THEN
zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj)
zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)
@@ -1007,8 +1007,8 @@ CONTAINS
zcpx(ji,jj) = 1.0_wp
ELSE IF(ll_tmp2) THEN
! no worries about ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm) = 0, it won't happen ! here
- zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
- & / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) )
+ zcpx(ji,jj) = ABS( ( ( ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) ) - ( ssh(ji,jj,Kmm) + ht_0(ji,jj) ) ) &
+ & / ( ssh(ji+1,jj,Kmm) - ssh(ji,jj,Kmm) ) )
zcpx(ji,jj) = MAX(MIN( zcpx(ji,jj) , 1.0_wp),0.0_wp)
ELSE
zcpx(ji,jj) = 0._wp
@@ -1026,8 +1026,8 @@ CONTAINS
zcpy(ji,jj) = 1.0_wp
ELSE IF(ll_tmp2) THEN
! no worries about ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm) = 0, it won't happen ! here
- zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
- & / (ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm)) )
+ zcpy(ji,jj) = ABS( ( ( ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) - ( ssh(ji,jj,Kmm) + ht_0(ji,jj) ) ) &
+ & / ( ssh(ji,jj+1,Kmm) - ssh(ji,jj,Kmm) ) )
zcpy(ji,jj) = MAX(MIN( zcpy(ji,jj) , 1.0_wp),0.0_wp)
ELSE
zcpy(ji,jj) = 0._wp
diff --git a/src/OCE/DYN/dynldf_iso.F90 b/src/OCE/DYN/dynldf_iso.F90
index 891d8842cc5955dd9f3ea41c56aaf819a160d321..8ee3d4f7611f0417bf0f5eef8acd08ab057995b4 100644
--- a/src/OCE/DYN/dynldf_iso.F90
+++ b/src/OCE/DYN/dynldf_iso.F90
@@ -200,8 +200,8 @@ CONTAINS
zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
ziut(ji,jj) = ( zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) ) &
- & + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &
- & +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk)
+ & + zcof1 * ( ( zdku (ji,jj) + zdk1u(ji-1,jj) ) & ! add () for NP repro
+ & +( zdk1u(ji,jj) + zdku (ji-1,jj) ) ) ) * tmask(ji,jj,jk)
END_2D
ELSE ! other coordinate system (zco or sco) : e3t
DO_2D( 0, 1, 0, 0 )
@@ -214,8 +214,8 @@ CONTAINS
zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
ziut(ji,jj) = ( zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) ) &
- & + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &
- & +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk)
+ & + zcof1 * ( ( zdku (ji,jj) + zdk1u(ji-1,jj) ) & ! add () for NP repro
+ & +( zdk1u(ji,jj) + zdku (ji-1,jj) ) ) ) * tmask(ji,jj,jk)
END_2D
ENDIF
@@ -230,8 +230,8 @@ CONTAINS
zcof2 = - zaht_0 * e1f(ji,jj) * zmskf * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
zjuf(ji,jj) = ( zabe2 * ( puu(ji,jj+1,jk,Kbb) - puu(ji,jj,jk,Kbb) ) &
- & + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) &
- & +zdk1u(ji,jj+1) + zdku (ji,jj) ) ) * fmask(ji,jj,jk)
+ & + zcof2 * ( ( zdku (ji,jj+1) + zdk1u(ji,jj) ) & ! add () for NP repro
+ & +( zdk1u(ji,jj+1) + zdku (ji,jj) ) ) ) * fmask(ji,jj,jk)
END_2D
! | t |
@@ -250,8 +250,8 @@ CONTAINS
zcof1 = - zaht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
zivf(ji,jj) = ( zabe1 * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji,jj,jk,Kbb) ) &
- & + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj) &
- & + zdk1v(ji,jj) + zdkv (ji+1,jj) ) ) * fmask(ji,jj,jk)
+ & + zcof1 * ( ( zdkv (ji,jj) + zdk1v(ji+1,jj) ) & ! add () for NP repro
+ & + ( zdk1v(ji,jj) + zdkv (ji+1,jj) ) ) ) * fmask(ji,jj,jk)
END_2D
! j-flux at t-point
@@ -267,8 +267,8 @@ CONTAINS
zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
zjvt(ji,jj) = ( zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) ) &
- & + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) &
- & +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk)
+ & + zcof2 * ( ( zdkv (ji,jj-1) + zdk1v(ji,jj) ) & ! add () for NP repro
+ & +( zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) ) * tmask(ji,jj,jk)
END_2D
ELSE ! other coordinate system (zco or sco) : e3t
DO_2D( 1, 0, 0, 1 )
@@ -281,8 +281,8 @@ CONTAINS
zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
zjvt(ji,jj) = ( zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) ) &
- & + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) &
- & +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk)
+ & + zcof2 * ( ( zdkv (ji,jj-1) + zdk1v(ji,jj) ) & ! add () for NP repro
+ & +( zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) ) * tmask(ji,jj,jk)
END_2D
ENDIF
@@ -290,12 +290,10 @@ CONTAINS
! Second derivative (divergence) and add to the general trend
! -----------------------------------------------------------
DO_2D( 0, 0, 0, 0 ) !!gm Question vectop possible??? !!bug
- puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( ziut(ji+1,jj) - ziut(ji,jj ) &
- & + zjuf(ji ,jj) - zjuf(ji,jj-1) ) * r1_e1e2u(ji,jj) &
- & / e3u(ji,jj,jk,Kmm)
- pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zivf(ji,jj ) - zivf(ji-1,jj) &
- & + zjvt(ji,jj+1) - zjvt(ji,jj ) ) * r1_e1e2v(ji,jj) &
- & / e3v(ji,jj,jk,Kmm)
+ puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( ( ziut(ji+1,jj) - ziut(ji,jj ) ) & ! add () for NP repro
+ & + ( zjuf(ji ,jj) - zjuf(ji,jj-1) ) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
+ pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( ( zivf(ji,jj ) - zivf(ji-1,jj) ) & ! add () for NP repro
+ & + ( zjvt(ji,jj+1) - zjvt(ji,jj ) ) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
END_2D
! ! ===============
END DO ! End of slab
@@ -366,10 +364,10 @@ CONTAINS
zcof3 = - e2u(ji,jj) * zmkt * zuwslpi
zcof4 = - e1u(ji,jj) * zmkf * zuwslpj
! vertical flux on u field
- zfuw(ji,jk) = zcof3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1) &
- & + zdiu (ji,jk ) + zdiu (ji+1,jk ) ) &
- & + zcof4 * ( zdj1u(ji,jk-1) + zdju (ji ,jk-1) &
- & + zdj1u(ji,jk ) + zdju (ji ,jk ) )
+ zfuw(ji,jk) = zcof3 * ( ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1) ) & ! add () for NP repro
+ & + ( zdiu (ji,jk ) + zdiu (ji+1,jk ) ) ) &
+ & + zcof4 * ( ( zdj1u(ji,jk-1) + zdju (ji ,jk-1) ) & ! add () for NP repro
+ & + ( zdj1u(ji,jk ) + zdju (ji ,jk ) ) )
! vertical mixing coefficient (akzu)
! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0
akzu(ji,jj,jk) = ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / zaht_0
@@ -392,10 +390,10 @@ CONTAINS
zcof3 = - e2v(ji,jj) * zmkf * zvwslpi
zcof4 = - e1v(ji,jj) * zmkt * zvwslpj
! vertical flux on v field
- zfvw(ji,jk) = zcof3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1) &
- & + zdiv (ji,jk ) + zdiv (ji-1,jk ) ) &
- & + zcof4 * ( zdjv (ji,jk-1) + zdj1v(ji ,jk-1) &
- & + zdjv (ji,jk ) + zdj1v(ji ,jk ) )
+ zfvw(ji,jk) = zcof3 * ( ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1) ) & ! add () for NP repro
+ & + ( zdiv (ji,jk ) + zdiv (ji-1,jk ) ) ) &
+ & + zcof4 * ( ( zdjv (ji,jk-1) + zdj1v(ji ,jk-1) ) & ! add () for NP repro
+ & + ( zdjv (ji,jk ) + zdj1v(ji ,jk ) ) )
! vertical mixing coefficient (akzv)
! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0
akzv(ji,jj,jk) = ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / zaht_0
@@ -407,10 +405,8 @@ CONTAINS
! -------------------------------------------------------------------
DO jk = 1, jpkm1
DO ji = ntsi, ntei
- puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) &
- & / e3u(ji,jj,jk,Kmm)
- pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) &
- & / e3v(ji,jj,jk,Kmm)
+ puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
+ pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
END DO
END DO
! ! ===============
diff --git a/src/OCE/DYN/dynldf_lap_blp.F90 b/src/OCE/DYN/dynldf_lap_blp.F90
index 09ae6fe78b107e89def6af8016a01409dffc6dea..2c23827f3ad18e7529099e3858c28f0e19d4b2bc 100644
--- a/src/OCE/DYN/dynldf_lap_blp.F90
+++ b/src/OCE/DYN/dynldf_lap_blp.F90
@@ -114,13 +114,13 @@ CONTAINS
!
DO_2D( iij-1, iij, iij-1, iij )
! ! ahm * e3 * curl (warning: computed for ji-1,jj-1)
- zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1) & ! ahmf already * by fmask
- & * ( e2v(ji ,jj-1) * pv(ji ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk) &
- & - e1u(ji-1,jj ) * pu(ji-1,jj ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk) )
+ zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1) & ! ahmf already * by fmask
+ & * ( ( e2v(ji ,jj-1) * pv(ji ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk) ) & ! add () for NP repro
+ & - ( e1u(ji-1,jj ) * pu(ji-1,jj ,jk) - e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk) ) )
! ! ahm * div (warning: computed for ji,jj)
zdiv(ji,jj) = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb) & ! ahmt already * by tmask
- & * ( e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk) &
- & + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk) )
+ & * ( ( e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk) ) &
+ & + ( e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk) ) )
END_2D
!
DO_2D( iij-1, iij-1, iij-1, iij-1 ) ! - curl( curl) + grad( div )
diff --git a/src/OCE/DYN/dynspg.F90 b/src/OCE/DYN/dynspg.F90
index ffd2219a2eff482503a8e24da3f039df87c80e55..3d635744ad477d9b9bdcca0f2a609693692b64ed 100644
--- a/src/OCE/DYN/dynspg.F90
+++ b/src/OCE/DYN/dynspg.F90
@@ -107,10 +107,10 @@ CONTAINS
IF( ln_apr_dyn .AND. .NOT.ln_dynspg_ts ) THEN !== Atmospheric pressure gradient (added later in time-split case) ==!
zg_2 = grav * 0.5
DO_2D( 0, 0, 0, 0 ) ! gradient of Patm using inverse barometer ssh
- zpgu(ji,jj) = zpgu(ji,jj) + zg_2 * ( ssh_ib (ji+1,jj) - ssh_ib (ji,jj) &
- & + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj)
- zpgv(ji,jj) = zpgv(ji,jj) + zg_2 * ( ssh_ib (ji,jj+1) - ssh_ib (ji,jj) &
- & + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj)
+ zpgu(ji,jj) = zpgu(ji,jj) + zg_2 * ( ( ssh_ib (ji+1,jj) - ssh_ib (ji,jj) ) & ! add () for NP repro
+ & + ( ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj) ) ) * r1_e1u(ji,jj)
+ zpgv(ji,jj) = zpgv(ji,jj) + zg_2 * ( ( ssh_ib (ji,jj+1) - ssh_ib (ji,jj) ) & ! add () for NP repro
+ & + ( ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj) ) ) * r1_e2v(ji,jj)
END_2D
ENDIF
!
diff --git a/src/OCE/DYN/dynspg_ts.F90 b/src/OCE/DYN/dynspg_ts.F90
index bd40bf5409e9d9b034118816df9ad6f59fc5cf35..ddb60a5068b70a06066150a1e3ef41f2e2f7edc0 100644
--- a/src/OCE/DYN/dynspg_ts.F90
+++ b/src/OCE/DYN/dynspg_ts.F90
@@ -559,7 +559,7 @@ CONTAINS
!-- ssh = ssh - delta_t' * [ frc + div( flux ) ] --!
!-------------------------------------------------------------------------!
DO_2D( 0, 0, 0, 0 )
- zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj)
+ zhdiv = ( ( zhU(ji,jj) - zhU(ji-1,jj) ) + ( zhV(ji,jj) - zhV(ji,jj-1) ) ) * r1_e1e2t(ji,jj)
ssha_e(ji,jj) = ( sshn_e(ji,jj) - rDt_e * ( ssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj)
END_2D
!
diff --git a/src/OCE/DYN/dynzdf.F90 b/src/OCE/DYN/dynzdf.F90
index fed10f23182fea2131de19cd81d668d71ed3aec5..b460e35fe2b099a6e111826cf931a84d060d34c2 100644
--- a/src/OCE/DYN/dynzdf.F90
+++ b/src/OCE/DYN/dynzdf.F90
@@ -141,18 +141,18 @@ CONTAINS
DO_1Di( 0, 0 ) ! Add bottom/top stress due to barotropic component only
iku = mbku(ji,jj) ! ocean bottom level at u- and v-points
ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
- puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + zDt_2 *( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) &
+ puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + zDt_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) &
& / e3u(ji,jj,iku,Kaa)
- pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + zDt_2 *( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) &
+ pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + zDt_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) &
& / e3v(ji,jj,ikv,Kaa)
END_1D
IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! Ocean cavities (ISF)
DO_1Di( 0, 0 )
iku = miku(ji,jj) ! top ocean level at u- and v-points
ikv = mikv(ji,jj) ! (first wet ocean u- and v-points)
- puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + zDt_2 *( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) &
+ puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + zDt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) &
& / e3u(ji,jj,iku,Kaa)
- pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + zDt_2 *( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) &
+ pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + zDt_2 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) &
& / e3v(ji,jj,ikv,Kaa)
END_1D
END IF
@@ -167,10 +167,10 @@ CONTAINS
CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu)
DO_2Dik( 0, 0, 1, jpkm1, 1 )
z1_e3ua = 1._wp / e3u(ji,jj,jk,Kaa) ! after scale factor at U-point
- zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) &
- & / e3uw(ji,jj,jk ,Kmm) * z1_e3ua * wumask(ji,jj,jk )
- zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) &
- & / e3uw(ji,jj,jk+1,Kmm) * z1_e3ua * wumask(ji,jj,jk+1)
+ zzwi = - zDt_2 * ( ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) + akzu(ji,jj,jk ) ) & ! add () for NP repro
+ & / e3uw(ji,jj,jk ,Kmm) * z1_e3ua * wumask(ji,jj,jk )
+ zzws = - zDt_2 * ( ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) + akzu(ji,jj,jk+1) ) & ! add () for NP repro
+ & / e3uw(ji,jj,jk+1,Kmm) * z1_e3ua * wumask(ji,jj,jk+1)
zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) * z1_e3ua
zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) * z1_e3ua
zwi(ji,jk) = zzwi + zDt_2 * MIN( zWui, 0._wp )
@@ -181,9 +181,9 @@ CONTAINS
DO_2Dik( 0, 0, 1, jpkm1, 1 )
z1_e3ua = 1._wp / e3u(ji,jj,jk,Kaa) ! after scale factor at U-point
zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) &
- & / e3uw(ji,jj,jk ,Kmm) * z1_e3ua * wumask(ji,jj,jk )
+ & / e3uw(ji,jj,jk ,Kmm) * z1_e3ua * wumask(ji,jj,jk )
zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) &
- & / e3uw(ji,jj,jk+1,Kmm) * z1_e3ua * wumask(ji,jj,jk+1)
+ & / e3uw(ji,jj,jk+1,Kmm) * z1_e3ua * wumask(ji,jj,jk+1)
zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) * z1_e3ua
zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) * z1_e3ua
zwi(ji,jk) = zzwi + zDt_2 * MIN( zWui, 0._wp )
@@ -196,7 +196,7 @@ CONTAINS
DO_1Di( 0, 0 ) !* Surface boundary conditions
zwi(ji,1) = 0._wp
zzws = - zDt_2 * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) &
- & / ( e3u(ji,jj,1,Kaa) * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2)
+ & / ( e3u(ji,jj,1,Kaa) * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2)
zWus = ( wi(ji ,jj,2) + wi(ji+1,jj,2) ) / e3u(ji,jj,1,Kaa)
zws(ji,1) = zzws - zDt_2 * MAX( zWus, 0._wp )
zwd(ji,1) = 1._wp - zzws - zDt_2 * ( MIN( zWus, 0._wp ) )
@@ -205,10 +205,10 @@ CONTAINS
SELECT CASE( nldf_dyn )
CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu)
DO_2Dik( 0, 0, 1, jpkm1, 1 )
- zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) &
- & / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk )
- zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) &
- & / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
+ zzwi = - zDt_2 * ( ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) + akzu(ji,jj,jk ) ) & ! add () for NP repro
+ & / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk )
+ zzws = - zDt_2 * ( ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) + akzu(ji,jj,jk+1) ) & ! add () for NP repro
+ & / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
zwi(ji,jk) = zzwi
zws(ji,jk) = zzws
zwd(ji,jk) = 1._wp - zzwi - zzws
@@ -216,9 +216,9 @@ CONTAINS
CASE DEFAULT ! iso-level lateral mixing
DO_2Dik( 0, 0, 1, jpkm1, 1 )
zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) &
- & / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk )
+ & / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk )
zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) &
- & / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
+ & / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
zwi(ji,jk) = zzwi
zws(ji,jk) = zzws
zwd(ji,jk) = 1._wp - zzwi - zzws
@@ -241,15 +241,13 @@ CONTAINS
IF ( ln_drgimp ) THEN ! implicit bottom friction
DO_1Di( 0, 0 )
iku = mbku(ji,jj) ! ocean bottom level at u- and v-points
- zwd(ji,iku) = zwd(ji,iku) - zDt_2 *( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) &
- & / e3u(ji,jj,iku,Kaa)
+ zwd(ji,iku) = zwd(ji,iku) - zDt_2 *( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / e3u(ji,jj,iku,Kaa)
END_1D
IF ( ln_isfcav.OR.ln_drgice_imp ) THEN ! top friction (always implicit)
DO_1Di( 0, 0 )
!!gm top Cd is masked (=0 outside cavities) no need of test on mik>=2 ==>> it has been suppressed
iku = miku(ji,jj) ! ocean top level at u- and v-points
- zwd(ji,iku) = zwd(ji,iku) - zDt_2 *( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) &
- & / e3u(ji,jj,iku,Kaa)
+ zwd(ji,iku) = zwd(ji,iku) - zDt_2 *( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / e3u(ji,jj,iku,Kaa)
END_1D
ENDIF
ENDIF
@@ -277,11 +275,11 @@ CONTAINS
#if defined key_RK3
! ! RK3: use only utau (not utau_b)
puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * utauU(ji,jj) &
- & / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1)
+ & / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1)
#else
! ! MLF: average of utau and utau_b
puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + zDt_2 * ( utau_b(ji,jj) + utauU(ji,jj) ) &
- & / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1)
+ & / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1)
#endif
END_1D
DO_2Dik( 0, 0, 2, jpkm1, 1 )
@@ -304,10 +302,10 @@ CONTAINS
CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzv)
DO_2Dik( 0, 0, 1, jpkm1, 1 )
z1_e3va = 1._wp / e3v(ji,jj,jk,Kaa) ! after scale factor at V-point
- zzwi = - zDt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) &
- & / e3vw(ji,jj,jk ,Kmm) * z1_e3va * wvmask(ji,jj,jk )
- zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) &
- & / e3vw(ji,jj,jk+1,Kmm) * z1_e3va * wvmask(ji,jj,jk+1)
+ zzwi = - zDt_2 * ( ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) + akzv(ji,jj,jk ) ) & ! add () for NP repro
+ & / e3vw(ji,jj,jk ,Kmm) * z1_e3va * wvmask(ji,jj,jk )
+ zzws = - zDt_2 * ( ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) + akzv(ji,jj,jk+1) ) & ! add () for NP repro
+ & / e3vw(ji,jj,jk+1,Kmm) * z1_e3va * wvmask(ji,jj,jk+1)
zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) * z1_e3va
zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * z1_e3va
zwi(ji,jk) = zzwi + zDt_2 * MIN( zWvi, 0._wp )
@@ -318,9 +316,9 @@ CONTAINS
DO_2Dik( 0, 0, 1, jpkm1, 1 )
z1_e3va = 1._wp / e3v(ji,jj,jk,Kaa) ! after scale factor at V-point
zzwi = - zDt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) &
- & / e3vw(ji,jj,jk ,Kmm) * z1_e3va * wvmask(ji,jj,jk )
+ & / e3vw(ji,jj,jk ,Kmm) * z1_e3va * wvmask(ji,jj,jk )
zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) &
- & / e3vw(ji,jj,jk+1,Kmm) * z1_e3va * wvmask(ji,jj,jk+1)
+ & / e3vw(ji,jj,jk+1,Kmm) * z1_e3va * wvmask(ji,jj,jk+1)
zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) * z1_e3va
zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * z1_e3va
zwi(ji,jk) = zzwi + zDt_2 * MIN( zWvi, 0._wp )
@@ -331,7 +329,7 @@ CONTAINS
DO_1Di( 0, 0 ) !* Surface boundary conditions
zwi(ji,1) = 0._wp
zzws = - zDt_2 * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) &
- & / ( e3v(ji,jj,1,Kaa) * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2)
+ & / ( e3v(ji,jj,1,Kaa) * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2)
zWvs = ( wi(ji,jj ,2) + wi(ji,jj+1,2) ) / e3v(ji,jj,1,Kaa)
zws(ji,1 ) = zzws - zDt_2 * MAX( zWvs, 0._wp )
zwd(ji,1 ) = 1._wp - zzws - zDt_2 * ( MIN( zWvs, 0._wp ) )
@@ -340,10 +338,10 @@ CONTAINS
SELECT CASE( nldf_dyn )
CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu)
DO_2Dik( 0, 0, 1, jpkm1, 1 )
- zzwi = - zDt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) &
- & / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk )
- zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) &
- & / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
+ zzwi = - zDt_2 * ( ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) + akzv(ji,jj,jk ) ) & ! add () for NP repro
+ & / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk )
+ zzws = - zDt_2 * ( ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) + akzv(ji,jj,jk+1) ) & ! add () for NP repro
+ & / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
zwi(ji,jk) = zzwi
zws(ji,jk) = zzws
zwd(ji,jk) = 1._wp - zzwi - zzws
@@ -351,9 +349,9 @@ CONTAINS
CASE DEFAULT ! iso-level lateral mixing
DO_2Dik( 0, 0, 1, jpkm1, 1 )
zzwi = - zDt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) &
- & / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk )
+ & / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk )
zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) &
- & / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
+ & / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
zwi(ji,jk) = zzwi
zws(ji,jk) = zzws
zwd(ji,jk) = 1._wp - zzwi - zzws
diff --git a/src/OCE/DYN/sshwzv.F90 b/src/OCE/DYN/sshwzv.F90
index 52027ff1beef6f8587a7d307267c539c48526a2f..e6660cf219f5744ecec084b171ab57c47fc0b624 100644
--- a/src/OCE/DYN/sshwzv.F90
+++ b/src/OCE/DYN/sshwzv.F90
@@ -191,7 +191,8 @@ CONTAINS
! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
- zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
+ zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) ) & ! add () for NP repro
+ & + ( vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) ) )
END_2D
END DO
IF( nn_hls == 1) CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp) ! - ML - Perhaps not necessary: not used for horizontal "connexions"
@@ -222,11 +223,10 @@ CONTAINS
#if defined key_qco
!!gm slightly faster
pww(ji,jj,jk) = pww(ji,jj,jk+1) - ( e3t(ji,jj,jk,Kmm) * hdiv(ji,jj,jk) &
- & + r1_Dt * e3t_0(ji,jj,jk) * ( r3t(ji,jj,Kaa) - r3t(ji,jj,Kbb) ) ) * tmask(ji,jj,jk)
+ & + r1_Dt * e3t_0(ji,jj,jk) * ( r3t(ji,jj,Kaa) - r3t(ji,jj,Kbb) ) ) * tmask(ji,jj,jk)
#else
pww(ji,jj,jk) = pww(ji,jj,jk+1) - ( e3t(ji,jj,jk,Kmm) * hdiv(ji,jj,jk) &
- & + r1_Dt * ( e3t(ji,jj,jk,Kaa) &
- & - e3t(ji,jj,jk,Kbb) ) ) * tmask(ji,jj,jk)
+ & + r1_Dt * ( e3t(ji,jj,jk,Kaa) - e3t(ji,jj,jk,Kbb) ) ) * tmask(ji,jj,jk)
#endif
END_3D
ENDIF
@@ -334,7 +334,8 @@ CONTAINS
! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
- zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
+ zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) ) & ! had () for NP repro
+ & + ( vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) ) )
END_2D
END DO
IF( nn_hls == 1) CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp) ! - ML - Perhaps not necessary: not used for horizontal "connexions"
diff --git a/src/OCE/DYN/wet_dry.F90 b/src/OCE/DYN/wet_dry.F90
index 465134503a2995b5d1544e27f1ee7727011e42ea..f1d5e0fdfcbb0029746dbd8b0e7301b6a0911170 100644
--- a/src/OCE/DYN/wet_dry.F90
+++ b/src/OCE/DYN/wet_dry.F90
@@ -179,10 +179,10 @@ CONTAINS
IF( tmask(ji,jj,1) < 0.5_wp ) CYCLE ! we don't care about land cells
IF( ht_0(ji,jj) - ssh_ref > zdepwd ) CYCLE ! and cells which are unlikely to dry
!
- zflxp(ji,jj) = MAX( zflxu(ji,jj) , 0._wp ) - MIN( zflxu(ji-1,jj ) , 0._wp ) &
- & + MAX( zflxv(ji,jj) , 0._wp ) - MIN( zflxv(ji, jj-1) , 0._wp )
- zflxn(ji,jj) = MIN( zflxu(ji,jj) , 0._wp ) - MAX( zflxu(ji-1,jj ) , 0._wp ) &
- & + MIN( zflxv(ji,jj) , 0._wp ) - MAX( zflxv(ji, jj-1) , 0._wp )
+ zflxp(ji,jj) = ( MAX( zflxu(ji,jj) , 0._wp ) - MIN( zflxu(ji-1,jj ) , 0._wp ) ) & ! add () for NP repro
+ & + ( MAX( zflxv(ji,jj) , 0._wp ) - MIN( zflxv(ji, jj-1) , 0._wp ) )
+ zflxn(ji,jj) = ( MIN( zflxu(ji,jj) , 0._wp ) - MAX( zflxu(ji-1,jj ) , 0._wp ) ) & ! add () for NP repro
+ & + ( MIN( zflxv(ji,jj) , 0._wp ) - MAX( zflxv(ji, jj-1) , 0._wp ) )
!
zdep2 = ht_0(ji,jj) + psshb1(ji,jj) - rn_wdmin1
IF( zdep2 <= 0._wp ) THEN ! add more safty, but not necessary
@@ -217,10 +217,10 @@ CONTAINS
!
ztmp = e1e2t(ji,jj)
!
- zzflxp = MAX( zflxu1(ji,jj) , 0._wp ) - MIN( zflxu1(ji-1,jj ) , 0._wp) &
- & + MAX( zflxv1(ji,jj) , 0._wp ) - MIN( zflxv1(ji, jj-1) , 0._wp)
- zzflxn = MIN( zflxu1(ji,jj) , 0._wp ) - MAX( zflxu1(ji-1,jj ) , 0._wp) &
- & + MIN( zflxv1(ji,jj) , 0._wp ) - MAX( zflxv1(ji, jj-1) , 0._wp)
+ zzflxp = ( MAX( zflxu1(ji,jj) , 0._wp ) - MIN( zflxu1(ji-1,jj ) , 0._wp) ) & ! add () for NP repro
+ & + ( MAX( zflxv1(ji,jj) , 0._wp ) - MIN( zflxv1(ji, jj-1) , 0._wp) )
+ zzflxn = ( MIN( zflxu1(ji,jj) , 0._wp ) - MAX( zflxu1(ji-1,jj ) , 0._wp) ) & ! add () for NP repro
+ & + ( MIN( zflxv1(ji,jj) , 0._wp ) - MAX( zflxv1(ji, jj-1) , 0._wp) )
!
zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
zdep2 = ht_0(ji,jj) + psshb1(ji,jj) - rn_wdmin1 - z2dt * psshemp(ji,jj)
@@ -312,10 +312,10 @@ CONTAINS
IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE ! we don't care about land cells
IF( ht_0(ji,jj) > zdepwd ) CYCLE ! and cells which are unlikely to dry
!
- zflxp(ji,jj) = MAX( zflxu(ji,jj) , 0._wp ) - MIN( zflxu(ji-1,jj ) , 0._wp ) &
- & + MAX( zflxv(ji,jj) , 0._wp ) - MIN( zflxv(ji, jj-1) , 0._wp )
- zflxn(ji,jj) = MIN( zflxu(ji,jj) , 0._wp ) - MAX( zflxu(ji-1,jj ) , 0._wp ) &
- & + MIN( zflxv(ji,jj) , 0._wp ) - MAX( zflxv(ji, jj-1) , 0._wp )
+ zflxp(ji,jj) = ( MAX( zflxu(ji,jj) , 0._wp ) - MIN( zflxu(ji-1,jj ) , 0._wp ) ) & ! add () for NP repro
+ & + ( MAX( zflxv(ji,jj) , 0._wp ) - MIN( zflxv(ji, jj-1) , 0._wp ) )
+ zflxn(ji,jj) = ( MIN( zflxu(ji,jj) , 0._wp ) - MAX( zflxu(ji-1,jj ) , 0._wp ) ) & ! add () for NP repro
+ & + ( MIN( zflxv(ji,jj) , 0._wp ) - MAX( zflxv(ji, jj-1) , 0._wp ) )
!
zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1
IF( zdep2 <= 0._wp ) THEN !add more safety, but not necessary
@@ -340,10 +340,10 @@ CONTAINS
!
ztmp = e1e2t(ji,jj)
!
- zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj), 0._wp) &
- & + max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji, jj-1), 0._wp)
- zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj), 0._wp) &
- & + min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji, jj-1), 0._wp)
+ zzflxp = ( MAX(zflxu1(ji,jj), 0._wp) - MIN(zflxu1(ji-1,jj), 0._wp) ) & ! add () for NP repro
+ & + ( MAX(zflxv1(ji,jj), 0._wp) - MIN(zflxv1(ji, jj-1), 0._wp) )
+ zzflxn = ( MIN(zflxu1(ji,jj), 0._wp) - MAX(zflxu1(ji-1,jj), 0._wp) ) & ! add () for NP repro
+ & + ( MIN(zflxv1(ji,jj), 0._wp) - MAX(zflxv1(ji, jj-1), 0._wp) )
zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj)
diff --git a/src/OCE/IOM/prtctl.F90 b/src/OCE/IOM/prtctl.F90
index 75103afdd6ab8138c3df99bb38af84622f9f95da..bb00edc5b9521d1caa08311ffaafe90d4571cc5a 100644
--- a/src/OCE/IOM/prtctl.F90
+++ b/src/OCE/IOM/prtctl.F90
@@ -132,8 +132,9 @@ CONTAINS
CHARACTER(len=*) , INTENT(in), OPTIONAL :: clinfo3
INTEGER , INTENT(in), OPTIONAL :: kdim
!
- CHARACTER(len=30) :: cl1, cl2
+ CHARACTER(len=30) :: cl1, cl2, cl3
CHARACTER(len=6) :: clfmt
+ CHARACTER(len=1) :: cli1
INTEGER :: jn, jl, kdir
INTEGER :: ipi1, ipi2, ipim1, ipim2
INTEGER :: isht1, isht2, ishtm1, ishtm2
@@ -295,8 +296,33 @@ CONTAINS
WRITE(inum, "(3x,a,' : ',"//clfmt//",3x,a,' : ',"//clfmt//")") cl1, zsum1, cl2, zsum2
ELSE
WRITE(inum, "(3x,a,' : ',"//clfmt//" )") cl1, zsum1
+
+ ENDIF
+ ! replace .false. by .true. to switch on theses prints of the last inner line
+ IF( .FALSE. .AND. l_IdoNFold .AND. jj1e == Nje0 - isht1 ) THEN
+ IF( PRESENT(tab2d_1) ) THEN
+ WRITE(cli1, '(i1)') INT(LOG10(REAL(ii1e-ii1s+1,wp))) + 1 ! how many digits to we need to write ?
+ WRITE(cl3, "(i"//cli1//")") ii1e-ii1s+1
+ WRITE(inum, "(a,"//TRIM(cl3)//clfmt//")") 'Last line '//TRIM(cl1)//' ', tab2d_1(ii1s:ii1e,jj1e)
+ ENDIF
+ IF( PRESENT(tab3d_1) ) THEN
+ WRITE(cli1, '(i1)') INT(LOG10(REAL((ii1e-ii1s+1)*kdir,wp))) + 1 ! how many digits to we need to write ?
+ WRITE(cl3, "(i"//cli1//")") (ii1e-ii1s+1)*kdir
+ WRITE(inum, "(a,"//TRIM(cl3)//clfmt//")") 'Last line '//TRIM(cl1)//' ', tab3d_1(ii1s:ii1e,jj1e,1:kdir)
+ ENDIF
+ ENDIF
+ IF( .FALSE. .AND. l_IdoNFold .AND. jj2e == Nje0 - isht2 ) THEN
+ IF( PRESENT(tab2d_2) ) THEN
+ WRITE(cli1, '(i1)') INT(LOG10(REAL(ii2e-ii2s+1,wp))) + 1 ! how many digits to we need to write ?
+ WRITE(cl3, "(i"//cli1//")") ii2e-ii2s+1
+ WRITE(inum, "(a,"//TRIM(cl3)//clfmt//")") 'Last line '//TRIM(cl2)//' ', tab2d_2(ii2s:ii2e,jj2e)
+ ENDIF
+ IF( PRESENT(tab3d_2) ) THEN
+ WRITE(cli1, '(i1)') INT(LOG10(REAL((ii2e-ii2s+1)*kdir,wp))) + 1 ! how many digits to we need to write ?
+ WRITE(cl3, "(i"//cli1//")") (ii2e-ii2s+1)*kdir
+ WRITE(inum, "(a,"//TRIM(cl3)//clfmt//")") 'Last line '//TRIM(cl2)//' ', tab3d_2(ii2s:ii2e,jj2e,1:kdir)
+ ENDIF
ENDIF
-
END DO
ENDIF
IF( jpnij == 1 ) CALL FLUSH(inum)
diff --git a/src/OCE/ISF/isfcav.F90 b/src/OCE/ISF/isfcav.F90
index 790324c4e570193c3da8dd2c3a0645526ed723d2..c352d594cd8d9cd705299331ee0298bc17cd8743 100644
--- a/src/OCE/ISF/isfcav.F90
+++ b/src/OCE/ISF/isfcav.F90
@@ -143,10 +143,10 @@ CONTAINS
!!gm better to do it like in the new zdfric.F90 i.e. avm weighted Ri computation
zcoef = 0.5_wp / e3w(ji,jj,ikt+1,Kmm)
! ! shear of horizontal velocity
- zdku = zcoef * ( uu(ji-1,jj ,ikt ,Kmm) + uu(ji,jj,ikt ,Kmm) &
- & -uu(ji-1,jj ,ikt+1,Kmm) - uu(ji,jj,ikt+1,Kmm) )
- zdkv = zcoef * ( vv(ji ,jj-1,ikt ,Kmm) + vv(ji,jj,ikt ,Kmm) &
- & -vv(ji ,jj-1,ikt+1,Kmm) - vv(ji,jj,ikt+1,Kmm) )
+ zdku = zcoef * ( ( uu(ji-1,jj ,ikt ,Kmm) + uu(ji,jj,ikt ,Kmm) ) & ! add () for NP repro
+ & - ( uu(ji-1,jj ,ikt+1,Kmm) + uu(ji,jj,ikt+1,Kmm) ) )
+ zdkv = zcoef * ( ( vv(ji ,jj-1,ikt ,Kmm) + vv(ji,jj,ikt ,Kmm) ) & ! add () for NP repro
+ & - ( vv(ji ,jj-1,ikt+1,Kmm) + vv(ji,jj,ikt+1,Kmm) ) )
! ! richardson number (minimum value set to zero)
zRc(ji,jj) = MAX( rn2(ji,jj,ikt+1), 1.e-20_wp ) / MAX( zdku*zdku + zdkv*zdkv, 1.e-20_wp )
END_2D
diff --git a/src/OCE/ISF/isfcpl.F90 b/src/OCE/ISF/isfcpl.F90
index 574875e3d23dcf0915a57eb96eb3b061e8115a06..a032701634ac8453529857a6a1768aa7cc6681fe 100644
--- a/src/OCE/ISF/isfcpl.F90
+++ b/src/OCE/ISF/isfcpl.F90
@@ -214,11 +214,9 @@ CONTAINS
!
zsummsk = zssmask0(jip1,jj) + zssmask0(jim1,jj) + zssmask0(ji,jjp1) + zssmask0(ji,jjm1)
!
- IF( zdssmask == 1._wp .AND. zsummsk /= 0._wp ) THEN
- ssh(ji,jj,Kmm)=( zssh(jip1,jj)*zssmask0(jip1,jj) &
- & + zssh(jim1,jj)*zssmask0(jim1,jj) &
- & + zssh(ji,jjp1)*zssmask0(ji,jjp1) &
- & + zssh(ji,jjm1)*zssmask0(ji,jjm1)) / zsummsk
+ IF (zdssmask == 1._wp .AND. zsummsk /= 0._wp) THEN
+ ssh(ji,jj,Kmm)=( ( zssh(jip1,jj)*zssmask0(jip1,jj) + zssh(jim1,jj)*zssmask0(jim1,jj) ) & ! add () for NP repro
+ & + ( zssh(ji,jjp1)*zssmask0(ji,jjp1) + zssh(ji,jjm1)*zssmask0(ji,jjm1) ) ) / zsummsk
zssmask_b(ji,jj) = 1._wp
ENDIF
END_2D
diff --git a/src/OCE/LDF/ldfc1d_c2d.F90 b/src/OCE/LDF/ldfc1d_c2d.F90
index efe9c3ade67439ac2d1804c950c096941b6b458d..0c01406bf03982ef7ce38ab64758b38281f2eda6 100644
--- a/src/OCE/LDF/ldfc1d_c2d.F90
+++ b/src/OCE/LDF/ldfc1d_c2d.F90
@@ -80,8 +80,8 @@ CONTAINS
pah1(:,:,jk) = pahs1(:,:) * ( zratio + zc * ( 1._wp + TANH( - ( gdept_0(:,:,jk) - zh ) * zw) ) )
END DO
DO_3DS( 0, 0, 0, 0, jpkm1, 1, -1 ) ! pah2 at F-point (zdep2 is an approximation in zps-coord.)
- zdep2 = ( gdept_0(ji,jj+1,jk) + gdept_0(ji+1,jj+1,jk) &
- & + gdept_0(ji,jj ,jk) + gdept_0(ji+1,jj ,jk) ) * r1_4
+ zdep2 = ( ( gdept_0(ji,jj+1,jk) + gdept_0(ji+1,jj+1,jk) ) & ! add () for NP repro
+ & + ( gdept_0(ji,jj ,jk) + gdept_0(ji+1,jj ,jk) ) ) * r1_4
pah2(ji,jj,jk) = pahs2(ji,jj) * ( zratio + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) ) )
END_3D
CALL lbc_lnk( 'ldfc1d_c2d', pah2, 'F', 1.0_wp ) ! Lateral boundary conditions
diff --git a/src/OCE/LDF/ldfdyn.F90 b/src/OCE/LDF/ldfdyn.F90
index 8f1e38f26ed1c596f2cce40c3e9d298036c53602..be4afd18cd463cb0395ab7e9a74ee18172286094 100644
--- a/src/OCE/LDF/ldfdyn.F90
+++ b/src/OCE/LDF/ldfdyn.F90
@@ -420,18 +420,14 @@ CONTAINS
DO jk = 1, jpkm1
!
DO_2D( 0, 1, 0, 1 )
- zdb = ( uu(ji,jj,jk,Kbb) * r1_e2u(ji,jj) - uu(ji-1,jj,jk,Kbb) * r1_e2u(ji-1,jj) ) &
- & * r1_e1t(ji,jj) * e2t(ji,jj) &
- & - ( vv(ji,jj,jk,Kbb) * r1_e1v(ji,jj) - vv(ji,jj-1,jk,Kbb) * r1_e1v(ji,jj-1) ) &
- & * r1_e2t(ji,jj) * e1t(ji,jj)
+ zdb = ( uu(ji,jj,jk,Kbb) * r1_e2u(ji,jj) - uu(ji-1,jj,jk,Kbb) * r1_e2u(ji-1,jj) ) * r1_e1t(ji,jj) * e2t(ji,jj) &
+ & - ( vv(ji,jj,jk,Kbb) * r1_e1v(ji,jj) - vv(ji,jj-1,jk,Kbb) * r1_e1v(ji,jj-1) ) * r1_e2t(ji,jj) * e1t(ji,jj)
dtensq(ji,jj,jk) = zdb * zdb * tmask(ji,jj,jk)
END_2D
!
DO_2D( 1, 0, 1, 0 )
- zdb = ( uu(ji,jj+1,jk,Kbb) * r1_e1u(ji,jj+1) - uu(ji,jj,jk,Kbb) * r1_e1u(ji,jj) ) &
- & * r1_e2f(ji,jj) * e1f(ji,jj) &
- & + ( vv(ji+1,jj,jk,Kbb) * r1_e2v(ji+1,jj) - vv(ji,jj,jk,Kbb) * r1_e2v(ji,jj) ) &
- & * r1_e1f(ji,jj) * e2f(ji,jj)
+ zdb = ( uu(ji,jj+1,jk,Kbb) * r1_e1u(ji,jj+1) - uu(ji,jj,jk,Kbb) * r1_e1u(ji,jj) ) * r1_e2f(ji,jj) * e1f(ji,jj) &
+ & + ( vv(ji+1,jj,jk,Kbb) * r1_e2v(ji+1,jj) - vv(ji,jj,jk,Kbb) * r1_e2v(ji,jj) ) * r1_e1f(ji,jj) * e2f(ji,jj)
dshesq(ji,jj,jk) = zdb * zdb * fmask(ji,jj,jk)
END_2D
!
@@ -445,9 +441,9 @@ CONTAINS
zu2pv2_ij_m1 = uu(ji-1,jj ,jk,Kbb) * uu(ji-1,jj ,jk,Kbb) + vv(ji ,jj-1,jk,Kbb) * vv(ji ,jj-1,jk,Kbb)
!
zdelta = zcmsmag * esqt(ji,jj) ! L^2 * (C_smag/pi)^2
- ahmt(ji,jj,jk) = zdelta * SQRT( dtensq(ji ,jj,jk) + &
- & r1_4 * ( dshesq(ji ,jj,jk) + dshesq(ji ,jj-1,jk) + &
- & dshesq(ji-1,jj,jk) + dshesq(ji-1,jj-1,jk) ) )
+ ahmt(ji,jj,jk) = zdelta * SQRT( dtensq(ji ,jj,jk) + &
+ & r1_4 * ( ( dshesq(ji ,jj,jk) + dshesq(ji ,jj-1,jk) ) + & ! add () for NP repro
+ & ( dshesq(ji-1,jj,jk) + dshesq(ji-1,jj-1,jk) ) ) )
ahmt(ji,jj,jk) = MAX( ahmt(ji,jj,jk), SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac * |U|L/2
ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == maxfac * L^2/(4*2dt)
!
@@ -455,13 +451,13 @@ CONTAINS
!
DO_2D( 0, 0, 0, 0 ) ! F-point value
!
- zu2pv2_ij_p1 = uu(ji ,jj+1,jk, kbb) * uu(ji ,jj+1,jk, kbb) + vv(ji+1,jj ,jk, kbb) * vv(ji+1,jj ,jk, kbb)
- zu2pv2_ij = uu(ji ,jj ,jk, kbb) * uu(ji ,jj ,jk, kbb) + vv(ji ,jj ,jk, kbb) * vv(ji ,jj ,jk, kbb)
+ zu2pv2_ij_p1 = uu(ji ,jj+1,jk,kbb) * uu(ji ,jj+1,jk,kbb) + vv(ji+1,jj ,jk,kbb) * vv(ji+1,jj ,jk,kbb)
+ zu2pv2_ij = uu(ji ,jj ,jk,kbb) * uu(ji ,jj ,jk,kbb) + vv(ji ,jj ,jk,kbb) * vv(ji ,jj ,jk,kbb)
!
zdelta = zcmsmag * esqf(ji,jj) ! L^2 * (C_smag/pi)^2
- ahmf(ji,jj,jk) = zdelta * SQRT( dshesq(ji ,jj,jk) + &
- & r1_4 * ( dtensq(ji ,jj,jk) + dtensq(ji ,jj+1,jk) + &
- & dtensq(ji+1,jj,jk) + dtensq(ji+1,jj+1,jk) ) )
+ ahmf(ji,jj,jk) = zdelta * SQRT( dshesq(ji ,jj,jk) + &
+ & r1_4 * ( ( dtensq(ji ,jj,jk) + dtensq(ji ,jj+1,jk) ) + & ! add () for NP repro
+ & ( dtensq(ji+1,jj,jk) + dtensq(ji+1,jj+1,jk) ) ) )
ahmf(ji,jj,jk) = MAX( ahmf(ji,jj,jk), SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac * |U|L/2
ahmf(ji,jj,jk) = MIN( ahmf(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == maxfac * L^2/(4*2dt)
!
diff --git a/src/OCE/LDF/ldfslp.F90 b/src/OCE/LDF/ldfslp.F90
index e6247b49f6124844a48ede07ed6c755638df4ffe..1f80e0f296ef671d4b27370420d5d83cb646d735 100644
--- a/src/OCE/LDF/ldfslp.F90
+++ b/src/OCE/LDF/ldfslp.F90
@@ -296,10 +296,10 @@ CONTAINS
& + 4.* zwz(ji ,jj ,jk) ) * zcofw
zwslpj = ( ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) ) & ! need additional () for
- & + ( zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) ) & ! reproducibility around NP
- & + 2.*( ( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) ) &
- & + ( zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) ) &
- & + 4.* zww(ji ,jj ,jk) ) * zcofw
+ & + ( zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) ) & ! reproducibility around NP
+ & + 2.*( ( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) ) &
+ & + ( zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) ) &
+ & + 4.* zww(ji ,jj ,jk) ) * zcofw
! !* decrease in vicinity of topography
zck = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) &
& * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25
diff --git a/src/OCE/LDF/ldftra.F90 b/src/OCE/LDF/ldftra.F90
index 22393cb006902429743325c755471b3bd3a130ff..b24ce4bae1642cfef003d7149e1e967c1b01f0fc 100644
--- a/src/OCE/LDF/ldftra.F90
+++ b/src/OCE/LDF/ldftra.F90
@@ -765,8 +765,8 @@ CONTAINS
pv(ji,jj,jk) = pv(ji,jj,jk) - ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) )
END_3D
DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
- pw(ji,jj,jk) = pw(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj ,jk) &
- & + zpsi_vw(ji,jj,jk) - zpsi_vw(ji ,jj-1,jk) )
+ pw(ji,jj,jk) = pw(ji,jj,jk) + ( ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj ,jk) ) & ! add () for NP repro
+ & + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji ,jj-1,jk) ) )
END_3D
!
! ! diagnose the eddy induced velocity and associated heat transport
@@ -838,16 +838,16 @@ CONTAINS
!
IF( iom_use('woce_eiv') ) THEN
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! e1 e2 w_eiv = dk[psix] + dk[psix]
- zw3d(ji,jj,jk) = ( psi_vw(ji,jj,jk) - psi_vw(ji ,jj-1,jk) &
- & + psi_uw(ji,jj,jk) - psi_uw(ji-1,jj ,jk) ) / e1e2t(ji,jj)
+ zw3d(ji,jj,jk) = ( ( psi_vw(ji,jj,jk) - psi_vw(ji ,jj-1,jk) ) & ! add () for NP repro
+ & + ( psi_uw(ji,jj,jk) - psi_uw(ji-1,jj ,jk) ) ) / e1e2t(ji,jj)
END_3D
CALL iom_put( "woce_eiv", zw3d )
ENDIF
!
IF( iom_use('weiv_masstr') ) THEN
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
- zw3d(ji,jj,jk) = rho0 * ( psi_vw(ji,jj,jk) - psi_vw(ji ,jj-1,jk) &
- & + psi_uw(ji,jj,jk) - psi_uw(ji-1,jj ,jk) )
+ zw3d(ji,jj,jk) = rho0 * ( ( psi_vw(ji,jj,jk) - psi_vw(ji ,jj-1,jk) ) & ! add () for NP repro
+ & + ( psi_uw(ji,jj,jk) - psi_uw(ji-1,jj ,jk) ) )
END_3D
CALL iom_put( "weiv_masstr" , zw3d ) ! mass transport in z-direction
ENDIF
diff --git a/src/OCE/SBC/sbcmod.F90 b/src/OCE/SBC/sbcmod.F90
index 4561f6a45b01b91c74a600b72664bbd1b38ae082..58b695eba89a0a820ce49c059f907d02b316b9ab 100644
--- a/src/OCE/SBC/sbcmod.F90
+++ b/src/OCE/SBC/sbcmod.F90
@@ -451,8 +451,8 @@ CONTAINS
!
ELSEIF( ln_wave .AND. ln_taw ) THEN ! Wave stress reduction
DO_2D( 0, 0, 0, 0 )
- utau(ji,jj) = utau(ji,jj) - tawx(ji,jj) + twox(ji,jj)
- vtau(ji,jj) = vtau(ji,jj) - tawy(ji,jj) + twoy(ji,jj)
+ utau(ji,jj) = utau(ji,jj) - ( tawx(ji,jj) - twox(ji,jj) )
+ vtau(ji,jj) = vtau(ji,jj) - ( tawy(ji,jj) - twoy(ji,jj) )
taum(ji,jj) = SQRT( utau(ji,jj)*utau(ji,jj) + vtau(ji,jj)*vtau(ji,jj) )
END_2D
!
diff --git a/src/OCE/SBC/sbcwave.F90 b/src/OCE/SBC/sbcwave.F90
index cb0b219ef973ccb16666824974e2887977e6908f..d27f0061e1869a3d9012bdb23c0ad42e8c92d1dc 100644
--- a/src/OCE/SBC/sbcwave.F90
+++ b/src/OCE/SBC/sbcwave.F90
@@ -221,11 +221,10 @@ CONTAINS
! !== vertical Stokes Drift 3D velocity ==!
!
DO_3D( 0, 1, 0, 1, 1, jpkm1 ) ! Horizontal e3*divergence
- ze3divh(ji,jj,jk) = ( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) * usd(ji ,jj,jk) &
- & - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * usd(ji-1,jj,jk) &
- & + e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) * vsd(ji,jj ,jk) &
- & - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vsd(ji,jj-1,jk) ) &
- & * r1_e1e2t(ji,jj)
+ ze3divh(ji,jj,jk) = ( ( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) * usd(ji ,jj,jk) & ! add () for NP repro
+ & - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * usd(ji-1,jj,jk) ) &
+ & + ( e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) * vsd(ji,jj ,jk) &
+ & - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vsd(ji,jj-1,jk) ) ) * r1_e1e2t(ji,jj)
END_3D
!
CALL lbc_lnk( 'sbcwave', ze3divh, 'T', 1.0_wp )
diff --git a/src/OCE/TRA/traadv_cen.F90 b/src/OCE/TRA/traadv_cen.F90
index 72761ac41ca8ec84f7c2008678e61d167a36b4d6..cb88ea15af0d165d7c427c616d07ac3ece4142a9 100644
--- a/src/OCE/TRA/traadv_cen.F90
+++ b/src/OCE/TRA/traadv_cen.F90
@@ -119,8 +119,8 @@ CONTAINS
END_2D
!
DO_2D( 0, 0, 0, 0 ) ! Horizontal divergence of advective fluxes
- pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zft_u(ji,jj) - zft_u(ji-1,jj ) &
- & + zft_v(ji,jj) - zft_v(ji ,jj-1) ) * r1_e1e2t(ji,jj) &
+ pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ( zft_u(ji,jj) - zft_u(ji-1,jj ) ) & ! add () for NP repro
+ & + ( zft_v(ji,jj) - zft_v(ji ,jj-1) ) ) * r1_e1e2t(ji,jj) &
& / e3t(ji,jj,jk,Kmm)
END_2D
END DO
@@ -251,8 +251,8 @@ CONTAINS
END_2D
!
DO_2D( 0, 0, 0, 0 ) ! Horizontal divergence of advective fluxes
- pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zft_u(ji,jj) - zft_u(ji-1,jj ) &
- & + zft_v(ji,jj) - zft_v(ji ,jj-1) ) * r1_e1e2t(ji,jj) &
+ pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ( zft_u(ji,jj) - zft_u(ji-1,jj ) ) & ! add () for NP repro
+ & + ( zft_v(ji,jj) - zft_v(ji ,jj-1) ) ) * r1_e1e2t(ji,jj) &
& / e3t(ji,jj,jk,Kmm)
END_2D
END DO
@@ -276,8 +276,8 @@ CONTAINS
END_2D
!
DO_2D( 0, 0, 0, 0 ) ! Horizontal divergence of advective fluxes
- pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zft_u(ji,jj) - zft_u(ji-1,jj ) &
- & + zft_v(ji,jj) - zft_v(ji ,jj-1) ) * r1_e1e2t(ji,jj) &
+ pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ( zft_u(ji,jj) - zft_u(ji-1,jj ) ) & ! add () for NP repro
+ & + ( zft_v(ji,jj) - zft_v(ji ,jj-1) ) ) * r1_e1e2t(ji,jj) &
& / e3t(ji,jj,jk,Kmm)
END_2D
END DO
@@ -459,8 +459,8 @@ CONTAINS
END SELECT
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !-- Divergence of advective fluxes --!
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) &
- & - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
- & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) &
+ & - ( ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) ) & ! add () for NP repro
+ & + ( zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) ) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
END_3D
!
diff --git a/src/OCE/TRA/traadv_fct.F90 b/src/OCE/TRA/traadv_fct.F90
index 7ab4ae15ef13c62a98290149174b68eccbfaaac1..bf99ad877249049dca123c762d789ddd2fc5e290 100644
--- a/src/OCE/TRA/traadv_fct.F90
+++ b/src/OCE/TRA/traadv_fct.F90
@@ -182,9 +182,9 @@ CONTAINS
!
DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) !* trend and after field with monotonic scheme
! ! total intermediate advective trends
- ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
- & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
- & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj)
+ ztra = - ( ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) ) & ! add () for NP reproducibility
+ & + ( zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) &
+ & + ( zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) ) * r1_e1e2t(ji,jj)
! ! update and guess with monotonic sheme
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra &
& / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk)
@@ -246,13 +246,9 @@ CONTAINS
zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm)
! ! C4 minus upstream advective fluxes
! round brackets added to fix the order of floating point operations
- ! needed to ensure halo 1 - halo 2 compatibility
- zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( zC2t_u + ( zltu(ji,jj,jk) - zltu(ji+1,jj,jk) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & ) - zwx(ji,jj,jk)
- zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( zC2t_v + ( zltv(ji,jj,jk) - zltv(ji,jj+1,jk) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & ) - zwy(ji,jj,jk)
+ ! needed to ensure the North Pole reproducibility
+ zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( zC2t_u + ( zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) ) - zwx(ji,jj,jk)
+ zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( zC2t_v + ( zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) ) - zwy(ji,jj,jk)
END_3D
!
CASE( 41 ) !- 4th order centered ==>> !!gm coding attempt need to be tested
@@ -306,9 +302,9 @@ CONTAINS
IF ( ll_zAimp ) THEN
DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) !* trend and after field with monotonic scheme
! ! total intermediate advective trends
- ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
- & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
- & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj)
+ ztra = - ( ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) ) & ! add () NP halo
+ & + ( zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) &
+ & + ( zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) ) * r1_e1e2t(ji,jj)
ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
END_3D
!
@@ -328,9 +324,9 @@ CONTAINS
! !== final trend with corrected fluxes ==!
!
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
- ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
- & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
- & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj)
+ ztra = - ( ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) ) & ! add () for NP reproducibility
+ & + ( zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) &
+ & + ( zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) ) * r1_e1e2t(ji,jj)
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra / e3t(ji,jj,jk,Kmm)
zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
END_3D
@@ -808,7 +804,9 @@ CONTAINS
tracer_flux_i(zwx_im1,zfp_ui_m1,zfm_ui_m1,ji-1,jj,jk)
tracer_flux_j(zwy_3d(ji,jj,jk),zfp_vj,zfm_vj,ji,jj,jk)
tracer_flux_j(zwy_jm1,zfp_vj_m1,zfm_vj_m1,ji,jj-1,jk)
- ztra = - ( zwx_3d(ji,jj,jk) - zwx_im1 + zwy_3d(ji,jj,jk) - zwy_jm1 + zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) * r1_e1e2t(ji,jj)
+ ztra = - ( ( zwx_3d(ji,jj,jk) - zwx_im1 ) & ! add () NP pole
+ & + ( zwy_3d(ji,jj,jk) - zwy_jm1 ) &
+ & + ( zwz( ji,jj,jk) - zwz(ji,jj,jk+1) ) ) * r1_e1e2t(ji,jj)
! ! update and guess with monotonic sheme
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra &
& / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk)
@@ -856,10 +854,10 @@ CONTAINS
! ! Laplacian
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! 2nd derivative * 1/ 6
! ! 1st derivative (gradient)
- ztu = ( pt(ji+1,jj,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
- ztu_im1 = ( pt(ji,jj,jk,jn,Kmm) - pt(ji-1,jj,jk,jn,Kmm) ) * umask(ji-1,jj,jk)
- ztv = ( pt(ji,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
- ztv_jm1 = ( pt(ji,jj,jk,jn,Kmm) - pt(ji,jj-1,jk,jn,Kmm) ) * vmask(ji,jj-1,jk)
+ ztu = ( pt(ji+1,jj,jk,jn,Kmm) - pt(ji ,jj,jk,jn,Kmm) ) * umask(ji ,jj,jk)
+ ztu_im1 = ( pt(ji ,jj,jk,jn,Kmm) - pt(ji-1,jj,jk,jn,Kmm) ) * umask(ji-1,jj,jk)
+ ztv = ( pt(ji,jj+1,jk,jn,Kmm) - pt(ji,jj ,jk,jn,Kmm) ) * vmask(ji,jj ,jk)
+ ztv_jm1 = ( pt(ji,jj ,jk,jn,Kmm) - pt(ji,jj-1,jk,jn,Kmm) ) * vmask(ji,jj-1,jk)
! ! 2nd derivative * 1/ 6
zltu_3d(ji,jj,jk) = ( ztu + ztu_im1 ) * r1_6
zltv_3d(ji,jj,jk) = ( ztv + ztv_jm1 ) * r1_6
@@ -873,13 +871,11 @@ CONTAINS
zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm)
! ! C4 minus upstream advective fluxes
! round brackets added to fix the order of floating point operations
- ! needed to ensure halo 1 - halo 2 compatibility
- zwx_3d(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( zC2t_u + ( zltu_3d(ji,jj,jk) - zltu_3d(ji+1,jj,jk) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & ) - zwx_3d(ji,jj,jk)
- zwy_3d(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( zC2t_v + ( zltv_3d(ji,jj,jk) - zltv_3d(ji,jj+1,jk) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & ) - zwy_3d(ji,jj,jk)
+ ! needed to ensure the North Pole reproducibility
+ zwx_3d(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( zC2t_u + ( zltu_3d(ji,jj,jk) - zltu_3d(ji+1,jj,jk) ) ) &
+ & - zwx_3d(ji,jj,jk)
+ zwy_3d(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( zC2t_v + ( zltv_3d(ji,jj,jk) - zltv_3d(ji,jj+1,jk) ) ) &
+ & - zwy_3d(ji,jj,jk)
END_3D
!
CASE( 41 ) !- 4th order centered ==>> !!gm coding attempt need to be tested
@@ -926,9 +922,9 @@ CONTAINS
IF ( ll_zAimp ) THEN
DO_3D( 1, 1, 1, 1, 1, jpkm1 ) !* trend and after field with monotonic scheme
! ! total intermediate advective trends
- ztra = - ( zwx_3d(ji,jj,jk) - zwx_3d(ji-1,jj ,jk ) &
- & + zwy_3d(ji,jj,jk) - zwy_3d(ji ,jj-1,jk ) &
- & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj)
+ ztra = - ( ( zwx_3d(ji,jj,jk) - zwx_3d(ji-1,jj ,jk ) ) & ! add () for NP reproducibility
+ & + ( zwy_3d(ji,jj,jk) - zwy_3d(ji ,jj-1,jk ) ) &
+ & + ( zwz( ji,jj,jk) - zwz( ji ,jj ,jk+1) ) ) * r1_e1e2t(ji,jj)
ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
END_3D
!
@@ -948,9 +944,9 @@ CONTAINS
! !== final trend with corrected fluxes ==!
!
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
- ztra = - ( zwx_3d(ji,jj,jk) - zwx_3d(ji-1,jj ,jk ) &
- & + zwy_3d(ji,jj,jk) - zwy_3d(ji ,jj-1,jk ) &
- & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj)
+ ztra = - ( ( zwx_3d(ji,jj,jk) - zwx_3d(ji-1,jj ,jk ) ) & ! add () for NP reproducibility
+ & + ( zwy_3d(ji,jj,jk) - zwy_3d(ji ,jj-1,jk ) ) &
+ & + ( zwz( ji,jj,jk) - zwz( ji ,jj ,jk+1) ) ) * r1_e1e2t(ji,jj)
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra / e3t(ji,jj,jk,Kmm)
zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk)
END_3D
diff --git a/src/OCE/TRA/traadv_mus.F90 b/src/OCE/TRA/traadv_mus.F90
index 314f4c992cdbfc742af4555212159e4c3a058354..0c63be9a6e7eac69e040ae9d7baa145928d75680 100644
--- a/src/OCE/TRA/traadv_mus.F90
+++ b/src/OCE/TRA/traadv_mus.F90
@@ -172,8 +172,8 @@ CONTAINS
END_2D
!
DO_2D( 0, 0, 0, 0 ) !-- Tracer advective trend
- pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj) - zwx(ji-1,jj ) &
- & + zwy(ji,jj) - zwy(ji ,jj-1) ) &
+ pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ( zwx(ji,jj) - zwx(ji-1,jj ) ) & ! ad () for NP repro
+ & + ( zwy(ji,jj) - zwy(ji ,jj-1) ) ) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
END_2D
END DO
@@ -391,8 +391,8 @@ CONTAINS
END_3D
!
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !-- Tracer advective trend
- pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
- & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) &
+ pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) ) & ! add () for NP repro
+ & + ( zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) ) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
END_3D
! ! trend diagnostics
diff --git a/src/OCE/TRA/traadv_ubs.F90 b/src/OCE/TRA/traadv_ubs.F90
index 94d7833301bd11b37ea18017eac8f3a15aaf2197..455701d3a4aaf3a2b3c084155c480e58f40c21fe 100644
--- a/src/OCE/TRA/traadv_ubs.F90
+++ b/src/OCE/TRA/traadv_ubs.F90
@@ -156,8 +156,8 @@ CONTAINS
zcenut = pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj ,jk,jn,Kmm) )
zcenvt = pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) )
! ! UBS advective fluxes
- ztu(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) )
- ztv(ji,jj,jk) = 0.5 * ( zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) )
+ ztu(ji,jj,jk) = 0.5 * ( zcenut - ( zfp_ui * zltu(ji,jj,jk) + zfm_ui * zltu(ji+1,jj,jk) ) ) ! add () for NP repro
+ ztv(ji,jj,jk) = 0.5 * ( zcenvt - ( zfp_vj * zltv(ji,jj,jk) + zfm_vj * zltv(ji,jj+1,jk) ) )
END_3D
!
DO_3D( 0, 0, 0, 0, 1, jpk )
@@ -166,9 +166,9 @@ CONTAINS
!
DO jk = 1, jpkm1 !== add the horizontal advective trend ==!
DO_2D( 0, 0, 0, 0 )
- pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) &
- & - ( ztu(ji,jj,jk) - ztu(ji-1,jj ,jk) &
- & + ztv(ji,jj,jk) - ztv(ji ,jj-1,jk) ) &
+ pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) &
+ & - ( ( ztu(ji,jj,jk) - ztu(ji-1,jj ,jk) ) & ! add () for NP repro
+ & + ( ztv(ji,jj,jk) - ztv(ji ,jj-1,jk) ) ) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
END_2D
!
diff --git a/src/OCE/TRA/trabbl.F90 b/src/OCE/TRA/trabbl.F90
index df82587f2f2e87d32267a669f0d5b5b9cc356bed..4c0408b740b7c9b5e0c612e9d6710e7024e25a01 100644
--- a/src/OCE/TRA/trabbl.F90
+++ b/src/OCE/TRA/trabbl.F90
@@ -196,10 +196,10 @@ CONTAINS
DO_2D( 0, 0, 0, 0 ) ! Compute the trend
ik = mbkt(ji,jj) ! bottom T-level index
pt_rhs(ji,jj,ik,jn) = pt_rhs(ji,jj,ik,jn) &
- & + ( ahu_bbl(ji ,jj ) * ( zptb(ji+1,jj ) - zptb(ji ,jj ) ) &
- & - ahu_bbl(ji-1,jj ) * ( zptb(ji ,jj ) - zptb(ji-1,jj ) ) &
- & + ahv_bbl(ji ,jj ) * ( zptb(ji ,jj+1) - zptb(ji ,jj ) ) &
- & - ahv_bbl(ji ,jj-1) * ( zptb(ji ,jj ) - zptb(ji ,jj-1) ) ) &
+ & + ( ( ahu_bbl(ji ,jj ) * ( zptb(ji+1,jj ) - zptb(ji ,jj ) ) & ! add () for NP repro
+ & - ahu_bbl(ji-1,jj ) * ( zptb(ji ,jj ) - zptb(ji-1,jj ) ) ) &
+ & + ( ahv_bbl(ji ,jj ) * ( zptb(ji ,jj+1) - zptb(ji ,jj ) ) &
+ & - ahv_bbl(ji ,jj-1) * ( zptb(ji ,jj ) - zptb(ji ,jj-1) ) ) ) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,ik,Kmm)
END_2D
! ! ===========
diff --git a/src/OCE/TRA/traldf_iso.F90 b/src/OCE/TRA/traldf_iso.F90
index d814e0584aeb6aa94991302f7fb1d57bc75258df..e8776e0d6482e6d95cbb55242a6f8b49a6a18ec1 100644
--- a/src/OCE/TRA/traldf_iso.F90
+++ b/src/OCE/TRA/traldf_iso.F90
@@ -178,18 +178,10 @@ CONTAINS
zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) &
& + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp )
!
- ! round brackets added to fix the order of floating point operations
- ! needed to ensure halo 1 - halo 2 compatibility
- zahu_w = ( ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & + ( pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & ) * zmsku
- zahv_w = ( ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & + ( pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & ) * zmskv
+ zahu_w = ( ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) ) & ! () for North Pole reproducibility
+ & + ( pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) ) * zmsku
+ zahv_w = ( ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) ) & ! () for North Pole reproducibility
+ & + ( pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) ) * zmskv
!
ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk) &
& + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk)
@@ -197,16 +189,14 @@ CONTAINS
!
IF( ln_traldf_msc ) THEN ! stabilizing vertical diffusivity coefficient
DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
- ! round brackets added to fix the order of floating point operations
- ! needed to ensure halo 1 - halo 2 compatibility
akz(ji,jj,jk) = 0.25_wp * ( &
- & ( ( pahu(ji ,jj,jk) + pahu(ji ,jj,jk-1) ) / ( e1u(ji ,jj) * e1u(ji ,jj) ) &
+ & ( ( pahu(ji ,jj,jk) + pahu(ji ,jj,jk-1) ) / ( e1u(ji ,jj) * e1u(ji ,jj) ) & ! () for NP repro
& + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) ) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & + ( ( pahv(ji,jj ,jk) + pahv(ji,jj ,jk-1) ) / ( e2v(ji,jj ) * e2v(ji,jj ) ) &
- & + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & )
+ & ) &
+ & + ( ( pahv(ji,jj ,jk) + pahv(ji,jj ,jk-1) ) / ( e2v(ji,jj ) * e2v(ji,jj ) ) & ! () for NP repro
+ & + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ) &
+ & ) &
+ & )
END_3D
!
IF( ln_traldf_blp ) THEN ! bilaplacian operator
@@ -288,30 +278,18 @@ CONTAINS
zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku
zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv
!
- ! round brackets added to fix the order of floating point operations
- ! needed to ensure halo 1 - halo 2 compatibility
- zftu(ji,jj,jk ) = ( zabe1 * zdit(ji,jj,jk) &
- & + zcof1 * ( ( zdkt (ji+1,jj) + zdk1t(ji,jj) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & + ( zdk1t(ji+1,jj) + zdkt (ji,jj) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & ) ) * umask(ji,jj,jk)
+ zftu(ji,jj,jk) = ( zabe1 * zdit(ji,jj,jk) &
+ & + zcof1 * ( ( zdkt (ji+1,jj) + zdk1t(ji,jj) ) & ! () for North Pole reproducibility
+ & + ( zdk1t(ji+1,jj) + zdkt (ji,jj) ) ) ) * umask(ji,jj,jk)
zftv(ji,jj,jk) = ( zabe2 * zdjt(ji,jj,jk) &
- & + zcof2 * ( ( zdkt (ji,jj+1) + zdk1t(ji,jj) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & + ( zdk1t(ji,jj+1) + zdkt (ji,jj) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & ) ) * vmask(ji,jj,jk)
+ & + zcof2 * ( ( zdkt (ji,jj+1) + zdk1t(ji,jj) ) & ! () for North Pole reproducibility
+ & + ( zdk1t(ji,jj+1) + zdkt (ji,jj) ) ) ) * vmask(ji,jj,jk)
END_2D
!
DO_2D( iij-1, iij-1, iij-1, iij-1 ) !== horizontal divergence and add to pta
- ! round brackets added to fix the order of floating point operations
- ! needed to ensure halo 1 - halo 2 compatibility
pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) &
- & + zsign * ( ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & + ( zftv(ji,jj,jk) - zftv(ji,jj-1,jk) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
+ & + zsign * ( ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) ) & ! () for North Pole reproducibility
+ & + ( zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) &
& ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
END_2D
END DO ! End of slab
@@ -332,26 +310,18 @@ CONTAINS
zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) &
& + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp )
!
- zahu_w = ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) &
- & + pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) * zmsku
- zahv_w = ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) &
- & + pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) * zmskv
+ zahu_w = ( ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) ) & ! () for North Pole reproducibility
+ & + ( pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) ) * zmsku
+ zahv_w = ( ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) ) & ! () for North Pole reproducibility
+ & + ( pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) ) * zmskv
!
zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk) !wslpi & j are already w-masked
zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk)
!
- ! round brackets added to fix the order of floating point operations
- ! needed to ensure halo 1 - halo 2 compatibility
- ztfw(ji,jj,jk) = zcoef3 * ( ( zdit(ji ,jj ,jk-1) + zdit(ji-1,jj ,jk) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & + ( zdit(ji-1,jj ,jk-1) + zdit(ji ,jj ,jk) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & ) &
- & + zcoef4 * ( ( zdjt(ji ,jj ,jk-1) + zdjt(ji ,jj-1,jk) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & + ( zdjt(ji ,jj-1,jk-1) + zdjt(ji ,jj ,jk) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & )
+ ztfw(ji,jj,jk) = zcoef3 * ( ( zdit(ji ,jj ,jk-1) + zdit(ji-1,jj ,jk) ) & ! () for North Pole reproducibility
+ & + ( zdit(ji-1,jj ,jk-1) + zdit(ji ,jj ,jk) ) ) &
+ & + zcoef4 * ( ( zdjt(ji ,jj ,jk-1) + zdjt(ji ,jj-1,jk) ) & ! () for North Pole reproducibility
+ & + ( zdjt(ji ,jj-1,jk-1) + zdjt(ji ,jj ,jk) ) )
END_3D
! !== add the vertical 33 flux ==!
IF( ln_traldf_lap ) THEN ! laplacian case: eddy coef = ah_wslp2 - akz
diff --git a/src/OCE/TRA/traldf_triad.F90 b/src/OCE/TRA/traldf_triad.F90
index 0b341db589874d92f39a7bf38cdde84856d575aa..b9697b43733d2f84d0a62643ee4ec3f3016c01cc 100644
--- a/src/OCE/TRA/traldf_triad.F90
+++ b/src/OCE/TRA/traldf_triad.F90
@@ -156,18 +156,18 @@ CONTAINS
zah = 0.25_wp * pahu(ji,jj,jk)
zah1 = 0.25_wp * pahu(ji-1,jj,jk)
! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth)
- zslope2 = triadi_g(ji,jj,jk,1,kp) + ( gdept(ji+1,jj,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)
+ zslope2 = triadi_g(ji,jj,jk,1,kp) + ( gdept(ji+1,jj,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) &
+ & * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)
zslope2 = zslope2 *zslope2
- zslope21 = triadi_g(ji,jj,jk,0,kp) + ( gdept(ji,jj,jk,Kmm) - gdept(ji-1,jj,jk,Kmm) ) * r1_e1u(ji-1,jj) * umask(ji-1,jj,jk+kp)
+ zslope21 = triadi_g(ji,jj,jk,0,kp) + ( gdept(ji,jj,jk,Kmm) - gdept(ji-1,jj,jk,Kmm) ) &
+ & * r1_e1u(ji-1,jj) * umask(ji-1,jj,jk+kp)
zslope21 = zslope21 *zslope21
! round brackets added to fix the order of floating point operations
- ! needed to ensure halo 1 - halo 2 compatibility
- ah_wslp2(ji,jj,jk+kp) = ah_wslp2(ji,jj,jk+kp) + ( zah * zbu * ze3wr * r1_e1e2t(ji,jj) * zslope2 &
- & + zah1 * zbu1 * ze3wr * r1_e1e2t(ji,jj) * zslope21 &
- & ) ! bracket for halo 1 - halo 2 compatibility
- akz (ji,jj,jk+kp) = akz (ji,jj,jk+kp) + ( zah * r1_e1u(ji,jj) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp) &
- + zah1 * r1_e1u(ji-1,jj) * r1_e1u(ji-1,jj) * umask(ji-1,jj,jk+kp) &
- & ) ! bracket for halo 1 - halo 2 compatibility
+ ! needed to ensure the North Pole reproducibility
+ ah_wslp2(ji,jj,jk+kp) = ah_wslp2(ji,jj,jk+kp) + ( zah * zbu * ze3wr * r1_e1e2t(ji,jj) * zslope2 & ! () for NP
+ & + zah1 * zbu1 * ze3wr * r1_e1e2t(ji,jj) * zslope21 ) ! repro
+ akz (ji,jj,jk+kp) = akz (ji,jj,jk+kp) + ( zah * r1_e1u(ji ,jj) * r1_e1u(ji ,jj) * umask(ji ,jj,jk+kp) &
+ + zah1 * r1_e1u(ji-1,jj) * r1_e1u(ji-1,jj) * umask(ji-1,jj,jk+kp) )
END_3D
END DO
!
@@ -180,18 +180,18 @@ CONTAINS
zah1 = 0.25_wp * pahv(ji,jj-1,jk)
! Subtract s-coordinate slope at t-points to give slope rel to s surfaces
! (do this by *adding* gradient of depth)
- zslope2 = triadj_g(ji,jj,jk,1,kp) + ( gdept(ji,jj+1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)
+ zslope2 = triadj_g(ji,jj,jk,1,kp) + ( gdept(ji,jj+1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) &
+ & * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)
zslope2 = zslope2 * zslope2
- zslope21 = triadj_g(ji,jj,jk,0,kp) + ( gdept(ji,jj,jk,Kmm) - gdept(ji,jj-1,jk,Kmm) ) * r1_e2v(ji,jj-1) * vmask(ji,jj-1,jk+kp)
+ zslope21 = triadj_g(ji,jj,jk,0,kp) + ( gdept(ji,jj,jk,Kmm) - gdept(ji,jj-1,jk,Kmm) ) &
+ & * r1_e2v(ji,jj-1) * vmask(ji,jj-1,jk+kp)
zslope21 = zslope21 * zslope21
! round brackets added to fix the order of floating point operations
- ! needed to ensure halo 1 - halo 2 compatibility
- ah_wslp2(ji,jj,jk+kp) = ah_wslp2(ji,jj,jk+kp) + ( zah * zbv * ze3wr * r1_e1e2t(ji,jj) * zslope2 &
- & + zah1 * zbv1 * ze3wr * r1_e1e2t(ji,jj) * zslope21 &
- & ) ! bracket for halo 1 - halo 2 compatibility
- akz (ji,jj,jk+kp) = akz (ji,jj,jk+kp) + ( zah * r1_e2v(ji,jj) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp) &
- & + zah1 * r1_e2v(ji,jj-1) * r1_e2v(ji,jj-1) * vmask(ji,jj-1,jk+kp) &
- & ) ! bracket for halo 1 - halo 2 compatibility
+ ! needed to ensure the North Pole reproducibility
+ ah_wslp2(ji,jj,jk+kp) = ah_wslp2(ji,jj,jk+kp) + ( zah * zbv * ze3wr * r1_e1e2t(ji,jj) * zslope2 & ! () for NP
+ & + zah1 * zbv1 * ze3wr * r1_e1e2t(ji,jj) * zslope21 ) ! repro
+ akz (ji,jj,jk+kp) = akz (ji,jj,jk+kp) + ( zah * r1_e2v(ji,jj ) * r1_e2v(ji,jj ) * vmask(ji,jj ,jk+kp) &
+ & + zah1 * r1_e2v(ji,jj-1) * r1_e2v(ji,jj-1) * vmask(ji,jj-1,jk+kp) )
END_3D
END DO
!
@@ -225,16 +225,12 @@ CONTAINS
DO kp = 0, 1
DO_3D( 1, 0, 1, 0, 1, jpkm1 )
- ! round brackets added to fix the order of floating point operations
- ! needed to ensure halo 1 - halo 2 compatibility
- zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp) &
- & + ( 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji,jj,jk,1,kp) &
- & + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji+1,jj,jk,0,kp) &
- & ) ! bracket for halo 1 - halo 2 compatibility
- zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp) &
- & + ( 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj,jk,1,kp) &
- & + 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj+1,jk,0,kp) &
- & ) ! bracket for halo 1 - halo 2 compatibility
+ zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp) &
+ & + ( 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji ,jj,jk,1,kp) & ! () for NP reproducibility
+ & + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji+1,jj,jk,0,kp) )
+ zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp) & ! () for NP reproducibility
+ & + ( 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj ,jk,1,kp) &
+ & + 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj+1,jk,0,kp) )
END_3D
END DO
CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm )
@@ -440,14 +436,10 @@ CONTAINS
ENDIF
! !== horizontal divergence and add to the general trend ==!
DO_2D( iij-1, iij-1, iij-1, iij-1 )
- ! round brackets added to fix the order of floating point operations
- ! needed to ensure halo 1 - halo 2 compatibility
- pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) &
- & + zsign * ( ( zftu(ji-1,jj ,jk) - zftu(ji,jj,jk) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & + ( zftv(ji,jj-1,jk) - zftv(ji,jj,jk) &
- & ) & ! bracket for halo 1 - halo 2 compatibility
- & ) / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) )
+ pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) &
+ & + zsign * ( ( zftu(ji-1,jj ,jk) - zftu(ji,jj,jk) ) & ! () for NP reproducibility
+ & + ( zftv(ji ,jj-1,jk) - zftv(ji,jj,jk) ) ) &
+ & / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) )
END_2D
!
END DO
@@ -456,7 +448,7 @@ CONTAINS
IF( ln_traldf_lap ) THEN ! laplacian case: eddy coef = ah_wslp2 - akz
DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 )
ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) &
- & * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) &
+ & * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) &
& * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )
END_3D
ELSE ! bilaplacian
@@ -466,9 +458,9 @@ CONTAINS
ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) &
& * ah_wslp2(ji,jj,jk) * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )
END_3D
- CASE( 2 ) ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt and pt2 gradients, resp.
+ CASE( 2 ) ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt and pt2 gradients, resp.
DO_3D( 0, 0, 0, 0, 2, jpkm1 )
- ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) &
+ ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) &
& * ( ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) ) &
& + akz (ji,jj,jk) * ( pt2(ji,jj,jk-1,jn) - pt2(ji,jj,jk,jn) ) )
END_3D
diff --git a/src/OCE/TRA/tramle.F90 b/src/OCE/TRA/tramle.F90
index a6f00b0bf72866c66e504760e7b3a4c00a1eae2e..ef4ff1eccc2e2591dfd31acb1ae46cb8b6abc28b 100644
--- a/src/OCE/TRA/tramle.F90
+++ b/src/OCE/TRA/tramle.F90
@@ -211,10 +211,10 @@ CONTAINS
ELSEIF( nn_mle == 1 ) THEN ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat)
DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
zpsim_u(ji,jj) = rc_f * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) &
- & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )
+ & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )
!
zpsim_v(ji,jj) = rc_f * zhv(ji,jj) * zhv(ji,jj) * e1_e2v(ji,jj) &
- & * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )
+ & * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )
END_2D
ENDIF
!
@@ -251,12 +251,12 @@ CONTAINS
! !== transport increased by the MLE induced transport ==!
DO jk = 1, ikmax
DO_2D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
- pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) )
+ pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) ) ! add () for NO repro
pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) )
END_2D
DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
- pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk) &
- & + zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj-1,jk) ) * wmask(ji,jj,1)
+ pw(ji,jj,jk) = pw(ji,jj,jk) - ( ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk) ) & ! add () for NO repro
+ & + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj-1,jk) ) ) * wmask(ji,jj,1)
END_2D
END DO
diff --git a/src/OCE/TRA/trazdf.F90 b/src/OCE/TRA/trazdf.F90
index 7123aeba37c7fe1419e207dd89407a28d569e25e..12e547056c12251c569b90372d191a48d42a9aaa 100644
--- a/src/OCE/TRA/trazdf.F90
+++ b/src/OCE/TRA/trazdf.F90
@@ -205,7 +205,7 @@ CONTAINS
DO_2Dik( 0, 0, 1, jpkm1, 1 )
zzwi = - rDt * zwt(ji,jk ) / e3w(ji,jj,jk ,Kmm)
zzws = - rDt * zwt(ji,jk+1) / e3w(ji,jj,jk+1,Kmm)
- zwd(ji,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws &
+ zwd(ji,jk) = e3t(ji,jj,jk,Kaa) - ( zzwi + zzws ) &
& + rDt * ( MAX( wi(ji,jj,jk ) , 0._wp ) &
& - MIN( wi(ji,jj,jk+1) , 0._wp ) )
zwi(ji,jk) = zzwi + rDt * MIN( wi(ji,jj,jk ) , 0._wp )
@@ -215,7 +215,7 @@ CONTAINS
DO_2Dik( 0, 0, 1, jpkm1, 1 )
zwi(ji,jk) = - rDt * zwt(ji,jk ) / e3w(ji,jj,jk,Kmm)
zws(ji,jk) = - rDt * zwt(ji,jk+1) / e3w(ji,jj,jk+1,Kmm)
- zwd(ji,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jk) - zws(ji,jk)
+ zwd(ji,jk) = e3t(ji,jj,jk,Kaa) - ( zwi(ji,jk) + zws(ji,jk) )
END_2D
ENDIF
!
diff --git a/src/OCE/nemogcm.F90 b/src/OCE/nemogcm.F90
index 9e5ebf364cb019878d34e1d4fd7fc44a94021891..36dce715459e1c908ca07fad9d4ec85273917b41 100644
--- a/src/OCE/nemogcm.F90
+++ b/src/OCE/nemogcm.F90
@@ -164,6 +164,7 @@ CONTAINS
!
DO WHILE( istp <= nitend .AND. nstop == 0 )
!
+ ncom_stp = istp
# if defined key_qco || defined key_linssh
# if defined key_RK3
CALL stp_RK3
diff --git a/src/OCE/stp2d.F90 b/src/OCE/stp2d.F90
index 06c08b84b7ce3051c0bbb3a5a99f7aa44b99c548..f8fc31d10c25680f80db7d6b398a3498d04cf1aa 100644
--- a/src/OCE/stp2d.F90
+++ b/src/OCE/stp2d.F90
@@ -172,10 +172,10 @@ CONTAINS
ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW)
zztmp = grav * r1_2
DO_2D( 0, 0, 0, 0 )
- Ue_rhs(ji,jj) = Ue_rhs(ji,jj) + zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) &
- & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj)
- Ve_rhs(ji,jj) = Ve_rhs(ji,jj) + zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) &
- & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj)
+ Ue_rhs(ji,jj) = Ue_rhs(ji,jj) + zztmp * ( ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) & ! add () for NP repro
+ & + ( ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) ) * r1_e1u(ji,jj)
+ Ve_rhs(ji,jj) = Ve_rhs(ji,jj) + zztmp * ( ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) & ! add () for NP repro
+ & + ( ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) ) * r1_e2v(ji,jj)
END_2D
ENDIF
ENDIF
diff --git a/src/OFF/nemogcm.F90 b/src/OFF/nemogcm.F90
index 9e83a72731fbf1f2ad7821634dbb064606b0e729..a0d33001c10db68199dea0b3dbf7a77191127795 100644
--- a/src/OFF/nemogcm.F90
+++ b/src/OFF/nemogcm.F90
@@ -121,6 +121,7 @@ CONTAINS
CALL iom_init( cxios_context ) ! iom_put initialization (must be done after nemo_init for AGRIF+XIOS+OASIS)
!
DO WHILE ( istp <= nitend .AND. nstop == 0 ) !== OFF time-stepping ==!
+ ncom_stp = istp
IF( ln_timing ) THEN
zstptiming = MPI_Wtime()
IF ( istp == ( nit000 + 1 ) ) elapsed_time = zstptiming
diff --git a/src/SAS/nemogcm.F90 b/src/SAS/nemogcm.F90
index 8e07126365d8773cbedae282431f64f3d001e0b3..8c940266602440197820b189ec8f0133fbb9f3ee 100644
--- a/src/SAS/nemogcm.F90
+++ b/src/SAS/nemogcm.F90
@@ -124,6 +124,7 @@ CONTAINS
#endif
!
DO WHILE( istp <= nitend .AND. nstop == 0 )
+ ncom_stp = istp
CALL stp
istp = istp + 1
END DO
diff --git a/tests/ADIAB_WAVE/MY_SRC/sbcwave.F90 b/tests/ADIAB_WAVE/MY_SRC/sbcwave.F90
index 070d04706dc0d7909ea39283425645813b6eeafc..7d54953254123cb0638dec6a2d1354cc2cdfe2a7 100644
--- a/tests/ADIAB_WAVE/MY_SRC/sbcwave.F90
+++ b/tests/ADIAB_WAVE/MY_SRC/sbcwave.F90
@@ -245,11 +245,11 @@ CONTAINS
! !== vertical Stokes Drift 3D velocity ==!
!
DO_3D( 0, 1, 0, 1, 1, jpkm1 ) ! Horizontal e3*divergence
- ze3divh(ji,jj,jk) = ( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) * usd(ji ,jj,jk) &
- & - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * usd(ji-1,jj,jk) &
- & + e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) * vsd(ji,jj ,jk) &
- & - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vsd(ji,jj-1,jk) ) &
- & * r1_e1e2t(ji,jj)
+ ze3divh(ji,jj,jk) = ( ( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) * usd(ji ,jj,jk) & ! add () for NP repro
+ & - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * usd(ji-1,jj,jk) ) &
+ & + ( e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) * vsd(ji,jj ,jk) &
+ & - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vsd(ji,jj-1,jk) ) &
+ & ) * r1_e1e2t(ji,jj)
END_3D
!
CALL lbc_lnk( 'sbcwave', ze3divh, 'T', 1.0_wp )