Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • nemo/nemo
  • sparonuz/nemo
  • hatfield/nemo
  • extdevs/nemo
4 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 513 additions and 469 deletions
src/OCE/TRA/tradmp.F90
View file @ 4dbfc9bb
......@@ -94,7 +94,7 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation
!
INTEGER :: ji, jj, jk, jn ! dummy loop indices
REAL(wp), DIMENSION(A2D(nn_hls),jpk,jpts) :: zts_dta
REAL(wp), DIMENSION(A2D(nn_hls),jpk,jpts) :: zts_dta
REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: zwrk
REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ztrdts
!!----------------------------------------------------------------------
......@@ -146,7 +146,7 @@ CONTAINS
!
! outputs (clem trunk)
IF( iom_use('hflx_dmp_cea') .OR. iom_use('sflx_dmp_cea') ) THEN
ALLOCATE( zwrk(A2D(nn_hls),jpk) ) ! Needed to handle expressions containing e3t when using key_qco or key_linssh
ALLOCATE( zwrk(A2D(0),jpk) ) ! Needed to handle expressions containing e3t when using key_qco or key_linssh
zwrk(:,:,:) = 0._wp
IF( iom_use('hflx_dmp_cea') ) THEN
......
src/OCE/TRA/traldf_iso.F90
View file @ 4dbfc9bb
......@@ -388,7 +388,7 @@ CONTAINS
!
! ! "Poleward" diffusive heat or salt transports (T-S case only)
! note sign is reversed to give down-gradient diffusive transports )
IF( l_ptr ) CALL dia_ptr_hst( jn, 'ldf', -zftv(:,:,:) )
IF( l_ptr ) CALL dia_ptr_hst( jn, 'ldf', -zftv(A2D(0),:) )
! ! Diffusive heat transports
IF( l_hst ) CALL dia_ar5_hst( jn, 'ldf', -zftu(:,:,:), -zftv(:,:,:) )
!
......
src/OCE/TRA/traldf_lap_blp.F90
View file @ 4dbfc9bb
......@@ -171,7 +171,7 @@ CONTAINS
IF( ( kpass == 1 .AND. .NOT.ln_traldf_blp ) .OR. & !== first pass only ( laplacian) ==!
( kpass == 2 .AND. ln_traldf_blp ) ) THEN !== 2nd pass only (bilaplacian) ==!
IF( l_ptr ) CALL dia_ptr_hst( jn, 'ldf', -ztv(:,:,:) )
IF( l_ptr ) CALL dia_ptr_hst( jn, 'ldf', -ztv(A2D(0),:) )
IF( l_hst ) CALL dia_ar5_hst( jn, 'ldf', -ztu(:,:,:), -ztv(:,:,:) )
ENDIF
! ! ==================
......
src/OCE/TRA/traldf_triad.F90
View file @ 4dbfc9bb
......@@ -485,7 +485,7 @@ CONTAINS
( kpass == 2 .AND. ln_traldf_blp ) ) THEN !== 2nd pass (bilaplacian) ==!
!
! ! "Poleward" diffusive heat or salt transports (T-S case only)
IF( l_ptr ) CALL dia_ptr_hst( jn, 'ldf', zftv(:,:,:) )
IF( l_ptr ) CALL dia_ptr_hst( jn, 'ldf', zftv(A2D(0),:) )
! ! Diffusive heat transports
IF( l_hst ) CALL dia_ar5_hst( jn, 'ldf', zftu(:,:,:), zftv(:,:,:) )
!
......
src/OCE/TRA/trasbc.F90
View file @ 4dbfc9bb
......@@ -100,7 +100,7 @@ CONTAINS
!
!!gm This should be moved into sbcmod.F90 module ? (especially now that ln_traqsr is read in namsbc namelist)
IF( .NOT.ln_traqsr ) THEN ! no solar radiation penetration
DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D_OVR( 0, 0, 0, 0 )
qns(ji,jj) = qns(ji,jj) + qsr(ji,jj) ! total heat flux in qns
qsr(ji,jj) = 0._wp ! qsr set to zero
END_2D
......@@ -121,24 +121,24 @@ CONTAINS
ENDIF
ELSE ! No restart or restart not found: Euler forward time stepping
zfact = 1._wp
DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D_OVR( 0, 0, 0, 0 )
sbc_tsc(ji,jj,:) = 0._wp
sbc_tsc_b(ji,jj,:) = 0._wp
END_2D
ENDIF
ELSE !* other time-steps: swap of forcing fields
zfact = 0.5_wp
DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D_OVR( 0, 0, 0, 0 )
sbc_tsc_b(ji,jj,:) = sbc_tsc(ji,jj,:)
END_2D
ENDIF
! !== Now sbc tracer content fields ==!
DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D_OVR( 0, 0, 0, 0 )
sbc_tsc(ji,jj,jp_tem) = r1_rho0_rcp * qns(ji,jj) ! non solar heat flux
sbc_tsc(ji,jj,jp_sal) = r1_rho0 * sfx(ji,jj) ! salt flux due to freezing/melting
END_2D
IF( ln_linssh ) THEN !* linear free surface
DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) !==>> add concentration/dilution effect due to constant volume cell
DO_2D_OVR( 0, 0, 0, 0 ) !==>> add concentration/dilution effect due to constant volume cell
sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_tem,Kmm)
sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_sal,Kmm)
END_2D !==>> output c./d. term
......@@ -275,7 +275,7 @@ CONTAINS
!!gm This should be moved into sbcmod.F90 module ? (especially now that ln_traqsr is read in namsbc namelist)
IF( .NOT.ln_traqsr .AND. kstg == 1) THEN ! no solar radiation penetration
DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D_OVR( 0, 0, 0, 0 )
qns(ji,jj) = qns(ji,jj) + qsr(ji,jj) ! total heat flux in qns
qsr(ji,jj) = 0._wp ! qsr set to zero
END_2D
......
src/OCE/TRA/trazdf.F90
View file @ 4dbfc9bb
......@@ -6,6 +6,7 @@ MODULE trazdf
!! History : 1.0 ! 2005-11 (G. Madec) Original code
!! 3.0 ! 2008-01 (C. Ethe, G. Madec) merge TRC-TRA
!! 4.0 ! 2017-06 (G. Madec) remove explict time-stepping option
!! 4.5 ! 2022-06 (G. Madec) refactoring to reduce memory usage (j-k-i loops)
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -22,7 +23,7 @@ MODULE trazdf
USE ldfslp ! lateral diffusion: iso-neutral slope
USE trd_oce ! trends: ocean variables
USE trdtra ! trends: tracer trend manager
USE eosbn2, ONLY: ln_SEOS, rn_b0
USE eosbn2 , ONLY: ln_SEOS, rn_b0
!
USE in_out_manager ! I/O manager
USE prtctl ! Print control
......@@ -77,7 +78,7 @@ CONTAINS
ENDIF
!
! !* compute lateral mixing trend and add it to the general trend
CALL tra_zdf_imp( kt, nit000, 'TRA', rDt, Kbb, Kmm, Krhs, pts, Kaa, jpts )
CALL tra_zdf_imp( 'TRA', Kbb, Kmm, Krhs, pts, Kaa, jpts )
!!gm WHY here ! and I don't like that !
! DRAKKAR SSS control {
......@@ -116,7 +117,7 @@ CONTAINS
END SUBROUTINE tra_zdf
SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, Kbb, Kmm, Krhs, pt, Kaa, kjpt )
SUBROUTINE tra_zdf_imp( cdtype, Kbb, Kmm, Krhs, pt, Kaa, kjpt )
!!----------------------------------------------------------------------
!! *** ROUTINE tra_zdf_imp ***
!!
......@@ -136,128 +137,177 @@ CONTAINS
!!
!! ** Action : - pt(:,:,:,:,Kaa) becomes the after tracer
!!---------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs, Kaa ! ocean time level indices
INTEGER , INTENT(in ) :: kit000 ! first time step index
CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
INTEGER , INTENT(in ) :: kjpt ! number of tracers
REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation
!
INTEGER :: ji, jj, jk, jn ! dummy loop indices
REAL(wp) :: zrhs, zzwi, zzws ! local scalars
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwi, zwt, zwd, zws
REAL(wp), DIMENSION(A1Di(0),jpk) :: zwi, zwt, zwd, zws
!!---------------------------------------------------------------------
!
! ! ============= !
DO jn = 1, kjpt ! tracer loop !
! ! ============= !
! Matrix construction
! --------------------
! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer
!
IF( ( cdtype == 'TRA' .AND. ( jn == jp_tem .OR. ( jn == jp_sal .AND. ln_zdfddm ) ) ) .OR. &
& ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN
! ! ================= !
DO_1Dj( 0, 0 ) ! i-k slices loop !
! ! ================= !
DO jn = 1, kjpt ! tracer loop !
! ! ================= !
!
! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers
IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN
DO_3D( 1, 1, 1, 1, 2, jpk )
zwt(ji,jj,jk) = avt(ji,jj,jk)
END_3D
ELSE
DO_3D( 1, 1, 1, 1, 2, jpk )
zwt(ji,jj,jk) = avs(ji,jj,jk)
END_3D
ENDIF
zwt(:,:,1) = 0._wp
! Matrix construction
! --------------------
! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer
!
IF( l_ldfslp ) THEN ! isoneutral diffusion: add the contribution
IF( ln_traldf_msc ) THEN ! MSC iso-neutral operator
DO_3D( 0, 0, 0, 0, 2, jpkm1 )
zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk)
END_3D
ELSE ! standard or triad iso-neutral operator
DO_3D( 0, 0, 0, 0, 2, jpkm1 )
zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk)
END_3D
IF( ( cdtype == 'TRA' .AND. ( jn == jp_tem .OR. ( jn == jp_sal .AND. ln_zdfddm ) ) ) .OR. &
& ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN
!
! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers
!
IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN ! use avt for temperature
!
IF( l_ldfslp ) THEN ! use avt + isoneutral diffusion contribution
IF( ln_traldf_msc ) THEN ! MSC iso-neutral operator
DO_2Dik( 0, 0, 2, jpk, 1 )
zwt(ji,jk) = avt(ji,jj,jk) + akz(ji,jj,jk)
END_2D
ELSE ! standard or triad iso-neutral operator
DO_2Dik( 0, 0, 2, jpk, 1 )
zwt(ji,jk) = avt(ji,jj,jk) + ah_wslp2(ji,jj,jk)
END_2D
ENDIF
ELSE ! use avt only
DO_2Dik( 0, 0, 2, jpk, 1 )
zwt(ji,jk) = avt(ji,jj,jk)
END_2D
ENDIF
!
ELSE ! use avs for salinty or passive tracers
!
IF( l_ldfslp ) THEN ! use avs + isoneutral diffusion contribution
IF( ln_traldf_msc ) THEN ! MSC iso-neutral operator
DO_2Dik( 0, 0, 2, jpk, 1 )
zwt(ji,jk) = avs(ji,jj,jk) + akz(ji,jj,jk)
END_2D
ELSE ! standard or triad iso-neutral operator
DO_2Dik( 0, 0, 2, jpk, 1 )
zwt(ji,jk) = avs(ji,jj,jk) + ah_wslp2(ji,jj,jk)
END_2D
ENDIF
ELSE !
DO_2Dik( 0, 0, 2, jpk, 1 )
zwt(ji,jk) = avs(ji,jj,jk)
END_2D
ENDIF
ENDIF
zwt(:,1) = 0._wp
!
! Diagonal, lower (i), upper (s) (including the bottom boundary condition since avt is masked)
IF( ln_zad_Aimp ) THEN ! Adaptive implicit vertical advection
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 &
& + 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 )
zws(ji,jk) = zzws - rDt * MAX( wi(ji,jj,jk+1) , 0._wp )
END_2D
ELSE
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)
END_2D
ENDIF
!
!!gm BUG?? : if edmfm is equivalent to a w ==>>> just add +/- rDt * edmfm(ji,jj,jk+1/jk )
!! but edmfm is at t-point !!!! crazy??? why not keep it at w-point????
!
IF( ln_zdfmfc ) THEN ! add upward Mass Flux in the matrix
DO_2Dik( 0, 0, 1, jpkm1, 1 )
zws(ji,jk) = zws(ji,jk) + e3t(ji,jj,jk,Kaa) * rDt * edmfm(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
zwd(ji,jk) = zwd(ji,jk) - e3t(ji,jj,jk,Kaa) * rDt * edmfm(ji,jj,jk ) / e3w(ji,jj,jk+1,Kmm)
END_2D
ENDIF
! DO_3D( 0, 0, 0, 0, 1, jpkm1 )
! edmfa(ji,jj,jk) = 0._wp
! edmfb(ji,jj,jk) = -edmfm(ji,jj,jk ) / e3w(ji,jj,jk+1,Kmm)
! edmfc(ji,jj,jk) = edmfm(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
! END_3D
!!gm BUG : level jpk never used in the inversion
! DO_2D( 0, 0, 0, 0 )
! edmfa(ji,jj,jpk) = -edmfm(ji,jj,jpk-1) / e3w(ji,jj,jpk,Kmm)
! edmfb(ji,jj,jpk) = edmfm(ji,jj,jpk ) / e3w(ji,jj,jpk,Kmm)
! edmfc(ji,jj,jpk) = 0._wp
! END_2D
!!
!!gm BUG ??? below e3t_Kmm should be used ?
!! or even no multiplication by e3t unless there is a bug in wi calculation
!!
! DO_3D( 0, 0, 0, 0, 1, jpkm1 )
!!gm edmfa = 0._wp except at jpk which is not used ==>> zdiagi update is useless !
! zdiagi(ji,jj,jk) = zdiagi(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfa(ji,jj,jk)
! zdiags(ji,jj,jk) = zdiags(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfc(ji,jj,jk)
! zdiagd(ji,jj,jk) = zdiagd(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfb(ji,jj,jk)
! END_3D
!!gm CALL diag_mfc( zwi, zwd, zws, rDt, Kaa )
!!gm SUBROUTINE diag_mfc( zdiagi, zdiagd, zdiags, p2dt, Kaa )
!
!! Matrix inversion from the first level
!!----------------------------------------------------------------------
! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
!
! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
! ( ... )( ... ) ( ... )
! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
!
! m is decomposed in the product of an upper and lower triangular matrix.
! The 3 diagonal terms are in 3d arrays: zwd, zws, zwi.
! Suffices i,s and d indicate "inferior" (below diagonal), diagonal
! and "superior" (above diagonal) components of the tridiagonal system.
! The solution will be in the 4d array pta.
! The 3d array zwt is used as a work space array.
! En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then
! used as a work space array: its value is modified.
!
DO_1Di( 0, 0 ) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k) ! done one for all passive tracers (so included in the IF instruction)
zwt(ji,1) = zwd(ji,1)
END_1D
DO_2Dik( 0, 0, 2, jpkm1, 1 )
zwt(ji,jk) = zwd(ji,jk) - zwi(ji,jk) * zws(ji,jk-1) / zwt(ji,jk-1)
END_2D
!
ENDIF
!
! Diagonal, lower (i), upper (s) (including the bottom boundary condition since avt is masked)
IF( ln_zad_Aimp ) THEN ! Adaptive implicit vertical advection
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zzwi = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk ,Kmm)
zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws &
& + p2dt * ( MAX( wi(ji,jj,jk ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) )
zwi(ji,jj,jk) = zzwi + p2dt * MIN( wi(ji,jj,jk ) , 0._wp )
zws(ji,jj,jk) = zzws - p2dt * MAX( wi(ji,jj,jk+1) , 0._wp )
END_3D
ELSE
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk,Kmm)
zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jj,jk) - zws(ji,jj,jk)
END_3D
IF( ln_zdfmfc ) THEN ! add Mass Flux to the RHS
DO_2Dik( 0, 0, 1, jpkm1, 1 )
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + edmftra(ji,jj,jk,jn)
END_2D
!!gm CALL rhs_mfc( pt(:,:,:,jn,Krhs), jn )
ENDIF
!
! Modification of diagonal to add MF scheme
IF ( ln_zdfmfc ) THEN
CALL diag_mfc( zwi, zwd, zws, p2dt, Kaa )
END IF
!
!! Matrix inversion from the first level
!!----------------------------------------------------------------------
! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
!
! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
! ( ... )( ... ) ( ... )
! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
!
! m is decomposed in the product of an upper and lower triangular matrix.
! The 3 diagonal terms are in 3d arrays: zwd, zws, zwi.
! Suffices i,s and d indicate "inferior" (below diagonal), diagonal
! and "superior" (above diagonal) components of the tridiagonal system.
! The solution will be in the 4d array pta.
! The 3d array zwt is used as a work space array.
! En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then
! used as a work space array: its value is modified.
!
DO_2D( 0, 0, 0, 0 ) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k) ! done one for all passive tracers (so included in the IF instruction)
zwt(ji,jj,1) = zwd(ji,jj,1)
DO_1Di( 0, 0 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1
pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb ) &
& + rDt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs)
END_1D
DO_2Dik( 0, 0, 2, jpkm1, 1 )
zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb ) &
& + rDt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs) ! zrhs=right hand side
pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jk) / zwt(ji,jk-1) * pt(ji,jj,jk-1,jn,Kaa)
END_2D
DO_3D( 0, 0, 0, 0, 2, jpkm1 )
zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
END_3D
!
ENDIF
!
! Modification of rhs to add MF scheme
IF ( ln_zdfmfc ) THEN
CALL rhs_mfc( pt(:,:,:,jn,Krhs), jn )
END IF
!
DO_2D( 0, 0, 0, 0 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1
pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) &
& + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs)
END_2D
DO_3D( 0, 0, 0, 0, 2, jpkm1 )
zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) &
& + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs) ! zrhs=right hand side
pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa)
END_3D
!
DO_2D( 0, 0, 0, 0 ) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk (result is the after tracer)
pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
END_2D
DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 )
pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jj,jk) * pt(ji,jj,jk+1,jn,Kaa) ) &
& / zwt(ji,jj,jk) * tmask(ji,jj,jk)
END_3D
DO_1Di( 0, 0 ) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk (result is the after tracer)
pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jpkm1) * tmask(ji,jj,jpkm1)
END_1D
DO_2Dik( 0, 0, jpk-2, 1, -1 )
pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jk) * pt(ji,jj,jk+1,jn,Kaa) ) &
& / zwt(ji,jk) * tmask(ji,jj,jk)
END_2D
! ! ================= !
END DO ! tracer loop !
! ! ================= !
END DO ! end tracer loop !
END_1D ! i-k slices loop !
! ! ================= !
END SUBROUTINE tra_zdf_imp
......
src/OCE/TRD/trddyn.F90
View file @ 4dbfc9bb
......@@ -56,10 +56,13 @@ CONTAINS
INTEGER , INTENT(in ) :: ktrd ! trend index
INTEGER , INTENT(in ) :: kt ! time step
INTEGER , INTENT(in ) :: Kmm ! time level index
INTEGER :: ji, jj, jk ! lopp indices
!!----------------------------------------------------------------------
!
putrd(:,:,:) = putrd(:,:,:) * umask(:,:,:) ! mask the trends
pvtrd(:,:,:) = pvtrd(:,:,:) * vmask(:,:,:)
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
putrd(ji,jj,jk) = putrd(ji,jj,jk) * umask(ji,jj,jk) ! mask the trends
pvtrd(ji,jj,jk) = pvtrd(ji,jj,jk) * vmask(ji,jj,jk)
END_3D
!
!!gm NB : here a lbc_lnk should probably be added
......@@ -120,10 +123,10 @@ CONTAINS
CALL iom_put( "vtrd_rvo", pvtrd )
CASE( jpdyn_keg ) ; CALL iom_put( "utrd_keg", putrd ) ! Kinetic Energy gradient (or had)
CALL iom_put( "vtrd_keg", pvtrd )
ALLOCATE( z3dx(jpi,jpj,jpk) , z3dy(jpi,jpj,jpk) )
z3dx(:,:,:) = 0._wp ! U.dxU & V.dyV (approximation)
z3dy(:,:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! no mask as un,vn are masked
ALLOCATE( z3dx(A2D(0),jpk) , z3dy(A2D(0),jpk) ) ! U.dxU & V.dyV (approximation)
z3dx(A2D(0),jpk) = 0._wp
z3dy(A2D(0),jpk) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! no mask as un,vn are masked
z3dx(ji,jj,jk) = uu(ji,jj,jk,Kmm) * ( uu(ji+1,jj,jk,Kmm) - uu(ji-1,jj,jk,Kmm) ) / ( 2._wp * e1u(ji,jj) )
z3dy(ji,jj,jk) = vv(ji,jj,jk,Kmm) * ( vv(ji,jj+1,jk,Kmm) - vv(ji,jj-1,jk,Kmm) ) / ( 2._wp * e2v(ji,jj) )
END_3D
......@@ -139,9 +142,11 @@ CONTAINS
CALL iom_put( "vtrd_zdf", pvtrd )
!
! ! wind stress trends
ALLOCATE( z2dx(jpi,jpj) , z2dy(jpi,jpj) )
z2dx(:,:) = ( utau_b(:,:) + utau(:,:) ) / ( e3u(:,:,1,Kmm) * rho0 )
z2dy(:,:) = ( vtau_b(:,:) + vtau(:,:) ) / ( e3v(:,:,1,Kmm) * rho0 )
ALLOCATE( z2dx(A2D(0)) , z2dy(A2D(0)) )
DO_2D( 0, 0, 0, 0 )
z2dx(ji,jj) = ( utau_b(ji,jj) + utauU(ji,jj) ) / ( e3u(ji,jj,1,Kmm) * rho0 )
z2dy(ji,jj) = ( vtau_b(ji,jj) + vtauV(ji,jj) ) / ( e3v(ji,jj,1,Kmm) * rho0 )
END_2D
CALL iom_put( "utrd_tau", z2dx )
CALL iom_put( "vtrd_tau", z2dy )
DEALLOCATE( z2dx , z2dy )
......
src/OCE/TRD/trdglo.F90
View file @ 4dbfc9bb
......@@ -42,6 +42,7 @@ MODULE trdglo
REAL(wp) :: tvolv ! volume of the whole ocean computed at v-points
REAL(wp) :: rpktrd ! potential to kinetic energy conversion
REAL(wp) :: peke ! conversion potential energy - kinetic energy trend
REAL(wp) :: xcof
! !!! domain averaged trends
REAL(wp), DIMENSION(jptot_tra) :: tmo, smo ! temperature and salinity trends
......@@ -77,44 +78,47 @@ CONTAINS
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ikbu, ikbv ! local integers
REAL(wp):: zvm, zvt, zvs, z1_2rho0 ! local scalars
REAL(wp), DIMENSION(jpi,jpj) :: ztswu, ztswv, z2dx, z2dy ! 2D workspace
REAL(wp), DIMENSION(A2D(0)) :: ztswu, ztswv, z2dx, z2dy ! 2D workspace
!!----------------------------------------------------------------------
!
IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
!
SELECT CASE( ctype )
!
CASE( 'TRA' ) !== Tracers (T & S) ==!
DO_3D( 1, 1, 1, 1, 1, jpkm1 ) ! global sum of mask volume trend and trend*T (including interior mask)
zvm = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) * tmask_i(ji,jj)
zvt = ptrdx(ji,jj,jk) * zvm
zvs = ptrdy(ji,jj,jk) * zvm
tmo(ktrd) = tmo(ktrd) + zvt
smo(ktrd) = smo(ktrd) + zvs
t2 (ktrd) = t2(ktrd) + zvt * ts(ji,jj,jk,jp_tem,Kmm)
s2 (ktrd) = s2(ktrd) + zvs * ts(ji,jj,jk,jp_sal,Kmm)
END_3D
SELECT CASE( ctype )
!
CASE( 'TRA' ) !== Tracers (T & S) ==!
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! global sum of mask volume trend and trend*T (including interior mask)
zvm = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) * tmask_i(ji,jj)
zvt = ptrdx(ji,jj,jk) * zvm
zvs = ptrdy(ji,jj,jk) * zvm
tmo(ktrd) = tmo(ktrd) + zvt
smo(ktrd) = smo(ktrd) + zvs
t2 (ktrd) = t2(ktrd) + zvt * ts(ji,jj,jk,jp_tem,Kmm)
s2 (ktrd) = s2(ktrd) + zvs * ts(ji,jj,jk,jp_sal,Kmm)
END_3D
! ! linear free surface: diagnose advective flux trough the fixed k=1 w-surface
IF( ln_linssh .AND. ktrd == jptra_zad ) THEN
tmo(jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_tem,Kmm) * e1e2t(:,:) * tmask_i(:,:) )
smo(jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_sal,Kmm) * e1e2t(:,:) * tmask_i(:,:) )
t2 (jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_tem,Kmm) * ts(:,:,1,jp_tem,Kmm) * e1e2t(:,:) * tmask_i(:,:) )
s2 (jptra_sad) = SUM( ww(:,:,1) * ts(:,:,1,jp_sal,Kmm) * ts(:,:,1,jp_sal,Kmm) * e1e2t(:,:) * tmask_i(:,:) )
ENDIF
IF( ln_linssh .AND. ktrd == jptra_zad ) THEN
DO_2D( 0, 0, 0, 0 ) ! global sum of mask volume trend and trend*T (including interior mask)
zvm = ww(ji,jj,1) * e1e2t(ji,jj) * tmask_i(ji,jj)
tmo(jptra_sad) = tmo(jptra_sad) + ts(ji,jj,1,jp_tem,Kmm) * zvm
smo(jptra_sad) = smo(jptra_sad) + ts(ji,jj,1,jp_sal,Kmm) * zvm
t2 (jptra_sad) = t2 (jptra_sad) + ts(ji,jj,1,jp_tem,Kmm) * ts(ji,jj,1,jp_tem,Kmm) * zvm
s2 (jptra_sad) = s2 (jptra_sad) + ts(ji,jj,1,jp_sal,Kmm) * ts(ji,jj,1,jp_sal,Kmm) * zvm
END_2D
ENDIF
!
IF( ktrd == jptra_atf ) THEN ! last trend (asselin time filter)
!
CALL glo_tra_wri( kt ) ! print the results in ocean.output
!
tmo(:) = 0._wp ! prepare the next time step (domain averaged array reset to zero)
smo(:) = 0._wp
t2 (:) = 0._wp
s2 (:) = 0._wp
!
ENDIF
IF( ktrd == jptra_atf ) THEN ! last trend (asselin time filter)
!
CALL glo_tra_wri( kt ) ! print the results in ocean.output
!
tmo(:) = 0._wp ! prepare the next time step (domain averaged array reset to zero)
smo(:) = 0._wp
t2 (:) = 0._wp
s2 (:) = 0._wp
!
ENDIF
!
CASE( 'DYN' ) !== Momentum and KE ==!
DO_3D( 1, 0, 1, 0, 1, jpkm1 )
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zvt = ptrdx(ji,jj,jk) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk) &
& * e1e2u (ji,jj) * e3u(ji,jj,jk,Kmm)
zvs = ptrdy(ji,jj,jk) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk) &
......@@ -126,11 +130,11 @@ CONTAINS
!
IF( ktrd == jpdyn_zdf ) THEN ! zdf trend: compute separately the surface forcing trend
z1_2rho0 = 0.5_wp / rho0
DO_2D( 1, 0, 1, 0 )
zvt = ( utau_b(ji,jj) + utau(ji,jj) ) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk) &
& * z1_2rho0 * e1e2u(ji,jj)
zvs = ( vtau_b(ji,jj) + vtau(ji,jj) ) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk) &
& * z1_2rho0 * e1e2v(ji,jj)
DO_2D( 0, 0, 0, 0 )
zvt = ( utau_b(ji,jj) + utauU(ji,jj) ) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk) &
& * z1_2rho0 * e1e2u(ji,jj)
zvs = ( vtau_b(ji,jj) + vtauV(ji,jj) ) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk) &
& * z1_2rho0 * e1e2v(ji,jj)
umo(jpdyn_tau) = umo(jpdyn_tau) + zvt
vmo(jpdyn_tau) = vmo(jpdyn_tau) + zvs
hke(jpdyn_tau) = hke(jpdyn_tau) + uu(ji,jj,1,Kmm) * zvt + vv(ji,jj,1,Kmm) * zvs
......@@ -184,8 +188,7 @@ CONTAINS
INTEGER, INTENT(in) :: Kmm ! time level index
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zcof ! local scalar
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zkx, zky, zkz, zkepe
REAL(wp), DIMENSION(A2D(0),jpk) :: zkx, zky, zkz, zkepe
!!----------------------------------------------------------------------
! I. Momentum trends
......@@ -196,24 +199,24 @@ CONTAINS
! I.1 Conversion potential energy - kinetic energy
! --------------------------------------------------
! c a u t i o n here, trends are computed at kt+1 (now , but after the swap)
zkx (:,:,:) = 0._wp
zky (:,:,:) = 0._wp
zkz (:,:,:) = 0._wp
zkepe(:,:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! loop from top to bottom
zkx (ji,jj,jk) = 0._wp
zky (ji,jj,jk) = 0._wp
zkz (ji,jj,jk) = 0._wp
zkepe(ji,jj,jk) = 0._wp
END_3D
CALL eos( ts(:,:,:,:,Kmm), rhd, rhop ) ! now potential density
zcof = 0.5_wp / rho0 ! Density flux at w-point
zkz(:,:,1) = 0._wp
DO jk = 2, jpk
zkz(:,:,jk) = zcof * e1e2t(:,:) * ww(:,:,jk) * ( rhop(:,:,jk) + rhop(:,:,jk-1) ) * tmask_i(:,:)
END DO
zkz(A2D(0),1) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpk ) ! loop from top to bottom
zkz(ji,jj,jk) = xcof * e1e2t(ji,jj) * ww(ji,jj,jk) * ( rhop(ji,jj,jk) + rhop(ji,jj,jk-1) ) * tmask_i(ji,jj)
END_3D
zcof = 0.5_wp / rho0 ! Density flux at u and v-points
DO_3D( 1, 0, 1, 0, 1, jpkm1 )
zkx(ji,jj,jk) = zcof * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) &
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zkx(ji,jj,jk) = xcof * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) &
& * uu(ji,jj,jk,Kmm) * ( rhop(ji,jj,jk) + rhop(ji+1,jj,jk) )
zky(ji,jj,jk) = zcof * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) &
zky(ji,jj,jk) = xcof * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) &
& * vv(ji,jj,jk,Kmm) * ( rhop(ji,jj,jk) + rhop(ji,jj+1,jk) )
END_3D
......@@ -227,10 +230,9 @@ CONTAINS
! I.2 Basin averaged kinetic energy trend
! ----------------------------------------
peke = 0._wp
DO jk = 1, jpkm1
peke = peke + SUM( zkepe(:,:,jk) * gdept(:,:,jk,Kmm) * e1e2t(:,:) &
& * e3t(:,:,jk,Kmm) )
END DO
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! loop from top to bottom
peke = peke + zkepe(ji,jj,jk) * gdept(ji,jj,jk,Kmm) * e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm)
END_3D
peke = grav * peke
! I.3 Sums over the global domain
......@@ -509,11 +511,14 @@ CONTAINS
WRITE(numout,*) '~~~~~~~~~~~~~'
ENDIF
xcof = 0.5_wp / rho0
! Total volume at t-points:
tvolt = 0._wp
DO jk = 1, jpkm1
tvolt = tvolt + SUM( e1e2t(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk) * tmask_i(:,:) )
END DO
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! Density flux divergence at t-point
tvolt = tvolt + e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) * tmask_i(ji,jj)
END_3D
CALL mpp_sum( 'trdglo', tvolt ) ! sum over the global domain
IF(lwp) WRITE(numout,*) ' total ocean volume at T-point tvolt = ',tvolt
......
src/OCE/TRD/trdken.F90
View file @ 4dbfc9bb
......@@ -30,6 +30,10 @@ MODULE trdken
IMPLICIT NONE
PRIVATE
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
PUBLIC trd_ken ! called by trddyn module
PUBLIC trd_ken_init ! called by trdini module
......@@ -38,9 +42,6 @@ MODULE trdken
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: bu, bv ! volume of u- and v-boxes
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: r1_bt ! inverse of t-box volume
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: trdken.F90 15104 2021-07-07 14:36:00Z clem $
......@@ -52,7 +53,7 @@ CONTAINS
!!---------------------------------------------------------------------
!! *** FUNCTION trd_ken_alloc ***
!!---------------------------------------------------------------------
ALLOCATE( bu(jpi,jpj,jpk) , bv(jpi,jpj,jpk) , r1_bt(jpi,jpj,jpk) , STAT= trd_ken_alloc )
ALLOCATE( bu(A2D(0),jpk) , bv(A2D(0),jpk) , r1_bt(A2D(0),jpk) , STAT= trd_ken_alloc )
!
CALL mpp_sum ( 'trdken', trd_ken_alloc )
IF( trd_ken_alloc /= 0 ) CALL ctl_stop( 'STOP', 'trd_ken_alloc: failed to allocate arrays' )
......@@ -77,31 +78,32 @@ CONTAINS
!
!
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: putrd, pvtrd ! U and V masked trends
INTEGER , INTENT(in ) :: ktrd ! trend index
INTEGER , INTENT(in ) :: kt ! time step
INTEGER , INTENT(in ) :: Kmm ! time level index
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: putrd, pvtrd ! U and V masked trends
INTEGER , INTENT(in ) :: ktrd ! trend index
INTEGER , INTENT(in ) :: kt ! time step
INTEGER , INTENT(in ) :: Kmm ! time level index
!
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ikbu , ikbv ! local integers
INTEGER :: ikbum1, ikbvm1 ! - -
REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z2dx, z2dy, zke2d ! 2D workspace
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zke ! 3D workspace
REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zke2d ! 2D workspace
REAL(wp) :: z2dx, z2dy, z2dxm1, z2dym1 ! 2D workspace
REAL(wp), DIMENSION(A2D(0),jpk) :: zke ! 3D workspace
!!----------------------------------------------------------------------
!
CALL lbc_lnk( 'trdken', putrd, 'U', -1.0_wp , pvtrd, 'V', -1.0_wp ) ! lateral boundary conditions
!
nkstp = kt
DO jk = 1, jpkm1
bu (:,:,jk) = e1e2u(:,:) * e3u(:,:,jk,Kmm)
bv (:,:,jk) = e1e2v(:,:) * e3v(:,:,jk,Kmm)
r1_bt(:,:,jk) = r1_e1e2t(:,:) / e3t(:,:,jk,Kmm) * tmask(:,:,jk)
DO_2D( 0, 0, 0, 0 )
bu (ji,jj,jk) = e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
bv (ji,jj,jk) = e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
r1_bt(ji,jj,jk) = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
END_2D
END DO
!
zke(:,:,jpk) = 0._wp
zke(1:nn_hls,:, : ) = 0._wp
zke(:,1:nn_hls, : ) = 0._wp
DO_3D( 0, nn_hls, 0, nn_hls, 1, jpkm1 )
zke(A2D(0),:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zke(ji,jj,jk) = 0.5_wp * rho0 *( uu(ji ,jj,jk,Kmm) * putrd(ji ,jj,jk) * bu(ji ,jj,jk) &
& + uu(ji-1,jj,jk,Kmm) * putrd(ji-1,jj,jk) * bu(ji-1,jj,jk) &
& + vv(ji,jj ,jk,Kmm) * pvtrd(ji,jj ,jk) * bv(ji,jj ,jk) &
......@@ -118,35 +120,31 @@ CONTAINS
CASE( jpdyn_ldf ) ; CALL iom_put( "ketrd_ldf" , zke ) ! lateral diffusion
CASE( jpdyn_zdf ) ; CALL iom_put( "ketrd_zdf" , zke ) ! vertical diffusion
! ! ! wind stress trends
ALLOCATE( z2dx(jpi,jpj) , z2dy(jpi,jpj) , zke2d(jpi,jpj) )
z2dx(:,:) = uu(:,:,1,Kmm) * ( utau_b(:,:) + utau(:,:) ) * e1e2u(:,:) * umask(:,:,1)
z2dy(:,:) = vv(:,:,1,Kmm) * ( vtau_b(:,:) + vtau(:,:) ) * e1e2v(:,:) * vmask(:,:,1)
zke2d(1:nn_hls,:) = 0._wp ; zke2d(:,1:nn_hls) = 0._wp
DO_2D( 0, nn_hls, 0, nn_hls )
zke2d(ji,jj) = r1_rho0 * 0.5_wp * ( z2dx(ji,jj) + z2dx(ji-1,jj) &
& + z2dy(ji,jj) + z2dy(ji,jj-1) ) * r1_bt(ji,jj,1)
END_2D
ALLOCATE( zke2d(A2D(0)) ) ; zke2d(A2D(0)) = 0._wp
DO_2D( 0, 0, 0, 0 )
z2dx = uu(ji,jj,1,Kmm) * ( utau_b(ji,jj) + utauU(ji,jj) ) * e1e2u(ji,jj) * umask(ji,jj,1)
z2dy = vv(ji,jj,1,Kmm) * ( vtau_b(ji,jj) + vtauV(ji,jj) ) * e1e2v(ji,jj) * vmask(ji,jj,1)
z2dxm1 = uu(ji-1,jj,1,Kmm) * ( utau_b(ji-1,jj) + utauU(ji-1,jj) ) * e1e2u(ji-1,jj) * umask(ji-1,jj,1)
z2dym1 = vv(ji,jj-1,1,Kmm) * ( vtau_b(ji,jj-1) + vtauV(ji,jj-1) ) * e1e2v(ji,jj-1) * vmask(ji,jj-1,1)
zke2d(ji,jj) = r1_rho0 * 0.5_wp * ( z2dx + z2dxm1 + z2dy + z2dym1 ) * r1_bt(ji,jj,1)
END_2D
CALL iom_put( "ketrd_tau" , zke2d ) !
DEALLOCATE( z2dx , z2dy , zke2d )
DEALLOCATE( zke2d )
CASE( jpdyn_bfr ) ; CALL iom_put( "ketrd_bfr" , zke ) ! bottom friction (explicit case)
!!gm TO BE DONE properly
!!gm only valid if ln_drgimp=F otherwise the bottom stress as to be recomputed at the end of the computation....
! IF(.NOT. ln_drgimp) THEN
! DO jj = 1, jpj !
! DO ji = 1, jpi
! ikbu = mbku(ji,jj) ! deepest ocean u- & v-levels
! ikbv = mbkv(ji,jj)
! z2dx(ji,jj) = uu(ji,jj,ikbu,Kmm) * bfrua(ji,jj) * uu(ji,jj,ikbu,Kmm)
! z2dy(ji,jj) = vv(ji,jj,ikbu,Kmm) * bfrva(ji,jj) * vv(ji,jj,ikbv,Kmm)
! END DO
! END DO
! zke2d(A2D(0)) = 0._wp
! DO_2D( 0, 0, 0, 0 )
! ikbu = mbku(ji,jj) ! deepest ocean u- & v-levels
! ikbv = mbkv(ji,jj)
! z2dx = uu(ji,jj,ikbu,Kmm) * bfrua(ji,jj) * uu(ji,jj,ikbu,Kmm)
! z2dy = vv(ji,jj,ikbu,Kmm) * bfrva(ji,jj) * vv(ji,jj,ikbv,Kmm)
! z2dxm1 = uu(ji-1,jj,ikbu,Kmm) * bfrua(ji-1,jj) * uu(ji-1,jj,ikbu,Kmm)
! z2dym1 = vv(ji,jj-1,ikbu,Kmm) * bfrva(ji,jj-1) * vv(ji,jj-1,ikbv,Kmm)
! zke2d(ji,jj) = 0.5_wp * ( z2dx + z2dxm1 + z2dy + z2dym1 ) * r1_bt(ji,jj,1)
! END_2D
! zke2d(1,:) = 0._wp ; zke2d(:,1) = 0._wp
! DO jj = 2, jpj
! DO ji = 2, jpi
! zke2d(ji,jj) = 0.5_wp * ( z2dx(ji,jj) + z2dx(ji-1,jj) &
! & + z2dy(ji,jj) + z2dy(ji,jj-1) ) * r1_bt(ji,jj, BEURK!!!
! END DO
! END DO
! CALL iom_put( "ketrd_bfr" , zke2d ) ! bottom friction (explicit case)
! ENDIF
!!gm end
......@@ -175,11 +173,11 @@ CONTAINS
CASE( jpdyn_ken ) ; ! kinetic energy
! called in dynnxt.F90 before asselin time filter
! with putrd=uu(Krhs) and pvtrd=vv(Krhs)
zke(:,:,:) = 0.5_wp * zke(:,:,:)
zke(A2D(0),:) = 0.5_wp * zke(A2D(0),:)
CALL iom_put( "KE", zke )
!
CALL ken_p2k( kt , zke, Kmm )
CALL iom_put( "ketrd_convP2K", zke ) ! conversion -rau*g*w
CALL iom_put( "ketrd_convP2K", zke ) ! conversion -rau*g*w
!
END SELECT
!
......@@ -203,22 +201,23 @@ CONTAINS
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: iku, ikv ! local integers
REAL(wp) :: zcoef ! local scalars
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zconv ! 3D workspace
REAL(wp), DIMENSION(A2D(0),jpk) :: zconv ! 3D workspace
!!----------------------------------------------------------------------
!
! Local constant initialization
zcoef = - rho0 * grav * 0.5_wp
! Surface value (also valid in partial step case)
zconv(:,:,1) = zcoef * ( 2._wp * rhd(:,:,1) ) * ww(:,:,1) * e3w(:,:,1,Kmm)
DO_2D( 0, 0, 0, 0 )
zconv(ji,jj,1) = zcoef * ( 2._wp * rhd(ji,jj,1) ) * ww(ji,jj,1) * e3w(ji,jj,1,Kmm)
END_2D
! interior value (2=<jk=<jpkm1)
DO jk = 2, jpk
zconv(:,:,jk) = zcoef * ( rhd(:,:,jk) + rhd(:,:,jk-1) ) * ww(:,:,jk) * e3w(:,:,jk,Kmm)
END DO
DO_3D( 0, 0, 0, 0 , 2, jpk )
zconv(ji,jj,jk) = zcoef * ( rhd(ji,jj,jk) + rhd(ji,jj,jk-1) ) * ww(ji,jj,jk) * e3w(ji,jj,jk,Kmm)
END_3D
! conv value on T-point
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zcoef = 0.5_wp / e3t(ji,jj,jk,Kmm)
pconv(ji,jj,jk) = zcoef * ( zconv(ji,jj,jk) + zconv(ji,jj,jk+1) ) * tmask(ji,jj,jk)
END_3D
......
src/OCE/TRD/trdmxl.F90
View file @ 4dbfc9bb
......@@ -45,20 +45,7 @@ MODULE trdmxl
PUBLIC trd_mxl_init ! routine called by opa.F90
PUBLIC trd_mxl_zint ! routine called by tracers routines
INTEGER :: nkstp ! current time step
!!gm to be moved from trdmxl_oce
! REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: hml ! ML depth (sum of e3t over nmln-1 levels) [m]
! REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: tml , sml ! now ML averaged T & S
! REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: tmlb_nf, smlb_nf ! not filtered before ML averaged T & S
!
!
! REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: hmlb, hmln ! before, now, and after Mixed Layer depths [m]
!
! REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: tb_mlb, tb_mln ! before (not filtered) tracer averaged over before and now ML
!
! REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: tn_mln ! now tracer averaged over now ML
!!gm end
INTEGER :: nkstp ! current time step
CHARACTER (LEN=40) :: clhstnam ! name of the trends NetCDF file
INTEGER :: nh_t, nmoymltrd
......@@ -111,44 +98,41 @@ CONTAINS
IF ( kt /= nkstp ) THEN != 1st call at kt time step =!
! !==============================!
nkstp = kt
! !== reset trend arrays to zero ==!
tmltrd(:,:,:) = 0._wp ; smltrd(:,:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpktrd )
tmltrd(ji,jj,jk) = 0._wp
smltrd(ji,jj,jk) = 0._wp
END_3D
!
wkx(:,:,:) = 0._wp !== now ML weights for vertical averaging ==!
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpktrd ) ! initialize wkx with vertical scale factor in mixed-layer
! !== now ML weights for vertical averaging ==!
DO_3D( 0, 0, 0, 0, 1, jpktrd ) ! initialize wkx with vertical scale factor in mixed-layer
IF( jk - kmxln(ji,jj) < 0 ) THEN
wkx(ji,jj,jk) = e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
ELSE
wkx(ji,jj,jk) = 0._wp
ENDIF
END_3D
!
hmxl(:,:) = 0._wp ! NOW mixed-layer depth
DO jk = 1, jpktrd
hmxl(:,:) = hmxl(:,:) + wkx(:,:,jk)
END DO
DO jk = 1, jpktrd ! integration weights
wkx(:,:,jk) = wkx(:,:,jk) / MAX( 1.e-20_wp, hmxl(:,:) ) * tmask(:,:,1)
END DO
DO_3D( 0, 0, 0, 0, 1, jpktrd )
hmxl(ji,jj) = hmxl(ji,jj) + wkx(ji,jj,jk)
wkx (ji,jj,jk) = wkx (ji,jj,jk) / MAX( 1.e-20_wp, hmxl(ji,jj) ) * tmask(ji,jj,1)
END_3D
!
! !== Vertically averaged T and S ==!
tml(:,:) = 0._wp ; sml(:,:) = 0._wp
DO jk = 1, jpktrd
tml(:,:) = tml(:,:) + wkx(:,:,jk) * ts(:,:,jk,jp_tem,Kmm)
sml(:,:) = sml(:,:) + wkx(:,:,jk) * ts(:,:,jk,jp_sal,Kmm)
END DO
DO_3D( 0, 0, 0, 0, 1, jpktrd )
tml(ji,jj) = tml(ji,jj) + wkx(ji,jj,jk) * ts(ji,jj,jk,jp_tem,Kmm)
sml(ji,jj) = sml(ji,jj) + wkx(ji,jj,jk) * ts(ji,jj,jk,jp_sal,Kmm)
END_3D
!
ENDIF
! mean now trends over the now ML
tmltrd(:,:,ktrd) = tmltrd(:,:,ktrd) + ptrdx(:,:,jk) * wkx(:,:,jk) ! temperature
smltrd(:,:,ktrd) = smltrd(:,:,ktrd) + ptrdy(:,:,jk) * wkx(:,:,jk) ! salinity
DO_2D( 0, 0, 0, 0 )
tmltrd(ji,jj,ktrd) = tmltrd(ji,jj,ktrd) + ptrdx(ji,jj,jk) * wkx(ji,jj,jk) ! temperature
smltrd(ji,jj,ktrd) = smltrd(ji,jj,ktrd) + ptrdy(ji,jj,jk) * wkx(ji,jj,jk) ! salinity
END_2D
!!gm to be put juste before the output !
......@@ -249,11 +233,11 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER , INTENT( in ) :: ktrd ! ocean trend index
CHARACTER(len=2) , INTENT( in ) :: ctype ! 2D surface/bottom or 3D interior physics
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pttrdmxl ! temperature trend
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pstrdmxl ! salinity trend
REAL(wp), DIMENSION(A2D(0),jpk), INTENT( in ) :: pttrdmxl ! temperature trend
REAL(wp), DIMENSION(A2D(0),jpk), INTENT( in ) :: pstrdmxl ! salinity trend
!
INTEGER :: ji, jj, jk, isum
REAL(wp), DIMENSION(jpi,jpj) :: zvlmsk
REAL(wp), DIMENSION(A2D(0)) :: zvlmsk
!!----------------------------------------------------------------------
! I. Definition of control surface and associated fields
......@@ -268,11 +252,13 @@ CONTAINS
! ... Set nmxl(ji,jj) = index of first T point below control surf. or outside mixed-layer
IF( nn_ctls == 0 ) THEN ! * control surface = mixed-layer with density criterion
nmxl(:,:) = nmln(:,:) ! array nmln computed in zdfmxl.F90
ELSEIF( nn_ctls == 1 ) THEN ! * control surface = read index from file
IF( nn_ctls == 0 ) THEN ! * control surface = mixed-layer with density criterion
DO_2D( 0, 0, 0, 0 )
nmxl(ji,jj) = nmln(ji,jj) ! array nmln computed in zdfmxl.F90
END_2D
ELSEIF( nn_ctls == 1 ) THEN ! * control surface = read index from file
nmxl(:,:) = nbol(:,:)
ELSEIF( nn_ctls >= 2 ) THEN ! * control surface = model level
ELSEIF( nn_ctls >= 2 ) THEN ! * control surface = model level
nn_ctls = MIN( nn_ctls, jpktrd - 1 )
nmxl(:,:) = nn_ctls + 1
ENDIF
......@@ -335,9 +321,9 @@ CONTAINS
REAL(wp) :: zavt, zfn, zfn2
! ! z(ts)mltot : dT/dt over the anlysis window (including Asselin)
! ! z(ts)mlres : residual = dh/dt entrainment term
REAL(wp), DIMENSION(jpi,jpj ) :: ztmltot , zsmltot , ztmlres , zsmlres , ztmlatf , zsmlatf
REAL(wp), DIMENSION(jpi,jpj ) :: ztmltot2, zsmltot2, ztmlres2, zsmlres2, ztmlatf2, zsmlatf2, ztmltrdm2, zsmltrdm2
REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmltrd2, zsmltrd2 ! only needed for mean diagnostics
REAL(wp), DIMENSION(A2D(0) ) :: ztmltot , zsmltot , ztmlres , zsmlres , ztmlatf , zsmlatf
REAL(wp), DIMENSION(A2D(0) ) :: ztmltot2, zsmltot2, ztmlres2, zsmlres2, ztmlatf2, zsmlatf2, ztmltrdm2, zsmltrdm2
REAL(wp), DIMENSION(A2D(0),jpk) :: ztmltrd2, zsmltrd2 ! only needed for mean diagnostics
!!----------------------------------------------------------------------
! ======================================================================
......@@ -446,12 +432,7 @@ CONTAINS
itmod = kt - nit000 + 1
MODULO_NTRD : IF( MOD( itmod, nn_trd ) == 0 ) THEN ! nitend MUST be multiple of nn_trd
!
ztmltot (:,:) = 0.e0 ; zsmltot (:,:) = 0.e0 ! reset arrays to zero
ztmlres (:,:) = 0.e0 ; zsmlres (:,:) = 0.e0
ztmltot2(:,:) = 0.e0 ; zsmltot2(:,:) = 0.e0
ztmlres2(:,:) = 0.e0 ; zsmlres2(:,:) = 0.e0
!
zfn = REAL( nmoymltrd, wp ) ; zfn2 = zfn * zfn
! III.1 Prepare fields for output ("instantaneous" diagnostics)
......@@ -469,25 +450,18 @@ CONTAINS
ztmlatf(:,:) = tmlatfm(:,:) - tmlatfn(:,:) + tmlatfb(:,:)
zsmlatf(:,:) = smlatfm(:,:) - smlatfn(:,:) + smlatfb(:,:)
!-- Lateral boundary conditions
! ... temperature ... ... salinity ...
CALL lbc_lnk( 'trdmxl', ztmltot , 'T', 1.0_wp, zsmltot , 'T', 1.0_wp, &
& ztmlres , 'T', 1.0_wp, zsmlres , 'T', 1.0_wp, &
& ztmlatf , 'T', 1.0_wp, zsmlatf , 'T', 1.0_wp )
! III.2 Prepare fields for output ("mean" diagnostics)
! ----------------------------------------------------
!-- Update the ML depth time sum (to build the Leap-Frog time mean)
hmxl_sum(:,:) = hmxlbn(:,:) + 2 * ( hmxl_sum(:,:) - hmxl(:,:) ) + hmxl(:,:)
hmxl_sum(:,:) = hmxlbn(:,:) + 2._wp * ( hmxl_sum(:,:) - hmxl(:,:) ) + hmxl(:,:)
!-- Compute temperature total trends
tml_sum (:,:) = tmlbn(:,:) + 2 * ( tml_sum(:,:) - tml(:,:) ) + tml(:,:)
tml_sum (:,:) = tmlbn(:,:) + 2._wp * ( tml_sum(:,:) - tml(:,:) ) + tml(:,:)
ztmltot2(:,:) = ( tml_sum(:,:) - tml_sumb(:,:) ) / p2dt ! now in degC/s
!-- Compute salinity total trends
sml_sum (:,:) = smlbn(:,:) + 2 * ( sml_sum(:,:) - sml(:,:) ) + sml(:,:)
sml_sum (:,:) = smlbn(:,:) + 2._wp * ( sml_sum(:,:) - sml(:,:) ) + sml(:,:)
zsmltot2(:,:) = ( sml_sum(:,:) - sml_sumb(:,:) ) / p2dt ! now in psu/s
!-- Compute temperature residuals
......@@ -495,7 +469,7 @@ CONTAINS
ztmltrd2(:,:,jl) = tmltrd_csum_ub(:,:,jl) + tmltrd_csum_ln(:,:,jl)
END DO
ztmltrdm2(:,:) = 0.e0
ztmltrdm2(:,:) = 0._wp
DO jl = 1, jpltrd
ztmltrdm2(:,:) = ztmltrdm2(:,:) + ztmltrd2(:,:,jl)
END DO
......@@ -508,7 +482,7 @@ CONTAINS
zsmltrd2(:,:,jl) = smltrd_csum_ub(:,:,jl) + smltrd_csum_ln(:,:,jl)
END DO
zsmltrdm2(:,:) = 0.
zsmltrdm2(:,:) = 0._wp
DO jl = 1, jpltrd
zsmltrdm2(:,:) = zsmltrdm2(:,:) + zsmltrd2(:,:,jl)
END DO
......@@ -519,13 +493,6 @@ CONTAINS
!-- Diagnose Asselin trend over the analysis window
ztmlatf2(:,:) = ztmltrd2(:,:,jpmxl_atf) - tmltrd_sum(:,:,jpmxl_atf) + tmltrd_atf_sumb(:,:)
zsmlatf2(:,:) = zsmltrd2(:,:,jpmxl_atf) - smltrd_sum(:,:,jpmxl_atf) + smltrd_atf_sumb(:,:)
!-- Lateral boundary conditions
! ... temperature ... ... salinity ...
CALL lbc_lnk( 'trdmxl', ztmltot2, 'T', 1.0_wp, zsmltot2, 'T', 1.0_wp, &
& ztmlres2, 'T', 1.0_wp, zsmlres2, 'T', 1.0_wp )
!
CALL lbc_lnk( 'trdmxl', ztmltrd2(:,:,:), 'T', 1.0_wp, zsmltrd2(:,:,:), 'T', 1.0_wp ) ! / in the NetCDF trends file
! III.3 Time evolution array swap
! -------------------------------
......@@ -622,6 +589,8 @@ CONTAINS
! ======================================================================
!-- Write the trends for T/S instantaneous diagnostics
! clem => these fields do not exist in field_def
IF( ln_trdmxl_instant ) THEN
......
src/OCE/TRD/trdmxl_oce.F90
View file @ 4dbfc9bb
......@@ -80,6 +80,8 @@ MODULE trdmxl_oce
smltrd_csum_ln, & !: ( idem for salinity )
smltrd_csum_ub !:
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: trdmxl_oce.F90 10425 2018-12-19 21:54:16Z smasson $
......@@ -101,29 +103,29 @@ CONTAINS
ierr(:) = 0
ALLOCATE( nmxl (jpi,jpj) , nbol (jpi,jpj), &
& wkx (jpi,jpj,jpk), hmxl (jpi,jpj), &
& tml (jpi,jpj) , sml (jpi,jpj), &
& tmlb (jpi,jpj) , smlb (jpi,jpj), &
& tmlbb(jpi,jpj) , smlbb(jpi,jpj), STAT = ierr(1) )
ALLOCATE( tmlbn(jpi,jpj) , smlbn(jpi,jpj), &
& tmltrdm(jpi,jpj), smltrdm(jpi,jpj), &
& tml_sum(jpi,jpj), tml_sumb(jpi,jpj),&
& tmltrd_atf_sumb(jpi,jpj) , STAT=ierr(2) )
ALLOCATE( sml_sum(jpi,jpj), sml_sumb(jpi,jpj), &
& smltrd_atf_sumb(jpi,jpj), &
& hmxl_sum(jpi,jpj), hmxlbn(jpi,jpj), &
& tmlatfb(jpi,jpj), tmlatfn(jpi,jpj), STAT = ierr(3) )
ALLOCATE( smlatfb(jpi,jpj), smlatfn(jpi,jpj), &
& tmlatfm(jpi,jpj), smlatfm(jpi,jpj), &
& tmltrd(jpi,jpj,jpltrd), smltrd(jpi,jpj,jpltrd), STAT=ierr(4))
ALLOCATE( tmltrd_sum(jpi,jpj,jpltrd),tmltrd_csum_ln(jpi,jpj,jpltrd), &
& tmltrd_csum_ub(jpi,jpj,jpltrd), smltrd_sum(jpi,jpj,jpltrd), &
& smltrd_csum_ln(jpi,jpj,jpltrd), smltrd_csum_ub(jpi,jpj,jpltrd), STAT=ierr(5) )
ALLOCATE( nmxl (A2D(0)) , nbol (A2D(0)), &
& wkx (A2D(0),jpk), hmxl (A2D(0)), &
& tml (A2D(0)) , sml (A2D(0)), &
& tmlb (A2D(0)) , smlb (A2D(0)), &
& tmlbb(A2D(0)) , smlbb(A2D(0)), STAT = ierr(1) )
ALLOCATE( tmlbn(A2D(0)) , smlbn(A2D(0)), &
& tmltrdm(A2D(0)), smltrdm(A2D(0)), &
& tml_sum(A2D(0)), tml_sumb(A2D(0)),&
& tmltrd_atf_sumb(A2D(0)) , STAT=ierr(2) )
ALLOCATE( sml_sum(A2D(0)), sml_sumb(A2D(0)), &
& smltrd_atf_sumb(A2D(0)), &
& hmxl_sum(A2D(0)), hmxlbn(A2D(0)), &
& tmlatfb(A2D(0)), tmlatfn(A2D(0)), STAT = ierr(3) )
ALLOCATE( smlatfb(A2D(0)), smlatfn(A2D(0)), &
& tmlatfm(A2D(0)), smlatfm(A2D(0)), &
& tmltrd(A2D(0),jpltrd), smltrd(A2D(0),jpltrd), STAT=ierr(4))
ALLOCATE( tmltrd_sum(A2D(0),jpltrd),tmltrd_csum_ln(A2D(0),jpltrd), &
& tmltrd_csum_ub(A2D(0),jpltrd), smltrd_sum(A2D(0),jpltrd), &
& smltrd_csum_ln(A2D(0),jpltrd), smltrd_csum_ub(A2D(0),jpltrd), STAT=ierr(5) )
!
trdmxl_oce_alloc = MAXVAL( ierr )
CALL mpp_sum ( 'trdmxl_oce', trdmxl_oce_alloc )
......
src/OCE/TRD/trdpen.F90
View file @ 4dbfc9bb
......@@ -35,6 +35,7 @@ MODULE trdpen
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: rab_pe ! partial derivatives of PE anomaly with respect to T and S
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
......@@ -48,7 +49,7 @@ CONTAINS
!!---------------------------------------------------------------------
!! *** FUNCTION trd_tra_alloc ***
!!---------------------------------------------------------------------
ALLOCATE( rab_pe(jpi,jpj,jpk,jpts) , STAT= trd_pen_alloc )
ALLOCATE( rab_pe(A2D(0),jpk,jpts) , STAT= trd_pen_alloc )
!
CALL mpp_sum ( 'trdpen', trd_pen_alloc )
IF( trd_pen_alloc /= 0 ) CALL ctl_stop( 'STOP', 'trd_pen_alloc: failed to allocate arrays' )
......@@ -69,37 +70,40 @@ CONTAINS
INTEGER , INTENT(in) :: Kmm ! time level index
REAL(wp) , INTENT(in) :: pdt ! time step [s]
!
INTEGER :: jk ! dummy loop indices
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z2d ! 2D workspace
REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpe ! 3D workspace
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z2d ! 2D workspace
REAL(wp), DIMENSION(A2D(0),jpk) :: zpe ! 3D workspace
!!----------------------------------------------------------------------
!
zpe(:,:,:) = 0._wp
zpe(A2D(0),:) = 0._wp
!
IF( kt /= nkstp ) THEN ! full eos: set partial derivatives at the 1st call of kt time step
nkstp = kt
CALL eos_pen( ts(:,:,:,:,Kmm), rab_PE, zpe, Kmm )
CALL eos_pen( ts(A2D(0),:,:,Kmm), rab_pe, zpe, Kmm )
CALL iom_put( "alphaPE", rab_pe(:,:,:,jp_tem) )
CALL iom_put( "betaPE" , rab_pe(:,:,:,jp_sal) )
CALL iom_put( "PEanom" , zpe )
ENDIF
!
zpe(:,:,jpk) = 0._wp
DO jk = 1, jpkm1
zpe(:,:,jk) = ( - ( rab_n(:,:,jk,jp_tem) + rab_pe(:,:,jk,jp_tem) ) * ptrdx(:,:,jk) &
& + ( rab_n(:,:,jk,jp_sal) + rab_pe(:,:,jk,jp_sal) ) * ptrdy(:,:,jk) )
END DO
zpe(A2D(0),jpk) = 0._wp
!
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
zpe(ji,jj,jk) = ( - ( rab_n(ji,jj,jk,jp_tem) + rab_pe(ji,jj,jk,jp_tem) ) * ptrdx(ji,jj,jk) &
& + ( rab_n(ji,jj,jk,jp_sal) + rab_pe(ji,jj,jk,jp_sal) ) * ptrdy(ji,jj,jk) )
END_3D
SELECT CASE ( ktrd )
CASE ( jptra_xad ) ; CALL iom_put( "petrd_xad", zpe ) ! zonal advection
CASE ( jptra_yad ) ; CALL iom_put( "petrd_yad", zpe ) ! merid. advection
CASE ( jptra_zad ) ; CALL iom_put( "petrd_zad", zpe ) ! vertical advection
IF( ln_linssh ) THEN ! cst volume : adv flux through z=0 surface
ALLOCATE( z2d(jpi,jpj) )
z2d(:,:) = ww(:,:,1) * ( &
& - ( rab_n(:,:,1,jp_tem) + rab_pe(:,:,1,jp_tem) ) * ts(:,:,1,jp_tem,Kmm) &
& + ( rab_n(:,:,1,jp_sal) + rab_pe(:,:,1,jp_sal) ) * ts(:,:,1,jp_sal,Kmm) &
& ) / e3t(:,:,1,Kmm)
ALLOCATE( z2d(A2D(0)) )
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = ww(ji,jj,1) * ( &
& - ( rab_n(ji,jj,1,jp_tem) + rab_pe(ji,jj,1,jp_tem) ) * ts(ji,jj,1,jp_tem,Kmm) &
& + ( rab_n(ji,jj,1,jp_sal) + rab_pe(ji,jj,1,jp_sal) ) * ts(ji,jj,1,jp_sal,Kmm) &
& ) / e3t(ji,jj,1,Kmm)
END_2D
CALL iom_put( "petrd_sad" , z2d )
DEALLOCATE( z2d )
ENDIF
......
src/OCE/TRD/trdtra.F90
View file @ 4dbfc9bb
......@@ -53,7 +53,7 @@ CONTAINS
!!---------------------------------------------------------------------
!! *** FUNCTION trd_tra_alloc ***
!!---------------------------------------------------------------------
ALLOCATE( trdtx(jpi,jpj,jpk) , trdty(jpi,jpj,jpk) , trdt(jpi,jpj,jpk) , avt_evd(jpi,jpj,jpk), STAT= trd_tra_alloc )
ALLOCATE( trdtx(A2D(0),jpk), trdty(A2D(0),jpk), trdt(A2D(0),jpk), avt_evd(A2D(0),jpk), STAT= trd_tra_alloc )
!
CALL mpp_sum ( 'trdtra', trd_tra_alloc )
IF( trd_tra_alloc /= 0 ) CALL ctl_stop( 'STOP', 'trd_tra_alloc: failed to allocate arrays' )
......@@ -73,24 +73,25 @@ CONTAINS
!! - 'TRA' case : regroup T & S trends
!! - send the trends to trd_tra_mng (trdtrc) for further processing
!!----------------------------------------------------------------------
INTEGER , INTENT(in) :: kt ! time step
CHARACTER(len=3) , INTENT(in) :: ctype ! tracers trends type 'TRA'/'TRC'
INTEGER , INTENT(in) :: ktra ! tracer index
INTEGER , INTENT(in) :: ktrd ! tracer trend index
INTEGER , INTENT(in) :: Kmm, Krhs ! time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: ptrd ! tracer trend or flux
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL :: pu ! now velocity
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL :: ptra ! now tracer variable
INTEGER , INTENT(in) :: kt ! time step
CHARACTER(len=3) , INTENT(in) :: ctype ! tracers trends type 'TRA'/'TRC'
INTEGER , INTENT(in) :: ktra ! tracer index
INTEGER , INTENT(in) :: ktrd ! tracer trend index
INTEGER , INTENT(in) :: Kmm, Krhs ! time level indices
REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ptrd ! tracer trend or flux
REAL(wp), DIMENSION(:,:,:), INTENT(in), OPTIONAL :: pu ! now velocity
REAL(wp), DIMENSION(:,:,:), INTENT(in), OPTIONAL :: ptra ! now tracer variable
!
INTEGER :: jk ! loop indices
INTEGER :: i01 ! 0 or 1
REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztrds ! 3D workspace
INTEGER :: ji,jj,jk ! loop indices
INTEGER :: i01 ! 0 or 1
REAL(wp) :: z1e3w
REAL(wp), DIMENSION(A2D(0),jpk) :: ztrds ! 3D workspace
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zwt, zws, ztrdt ! 3D workspace
!!----------------------------------------------------------------------
!
IF( .NOT. ALLOCATED( trdtx ) ) THEN ! allocate trdtra arrays
IF( trd_tra_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'trd_tra : unable to allocate arrays' )
avt_evd(:,:,:) = 0._wp
avt_evd(A2D(0),:) = 0._wp
ENDIF
!
i01 = COUNT( (/ PRESENT(pu) .OR. ( ktrd /= jptra_xad .AND. ktrd /= jptra_yad .AND. ktrd /= jptra_zad ) /) )
......@@ -103,13 +104,13 @@ CONTAINS
CASE( jptra_yad ) ; CALL trd_tra_adv( ptrd, pu, ptra, 'Y', trdty, Kmm )
CASE( jptra_zad ) ; CALL trd_tra_adv( ptrd, pu, ptra, 'Z', trdt, Kmm )
CASE( jptra_bbc, & ! qsr, bbc: on temperature only, send to trd_tra_mng
& jptra_qsr ) ; trdt(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
ztrds(:,:,:) = 0._wp
& jptra_qsr ) ; trdt(A2D(0),:) = ptrd(A2D(0),:) * tmask(A2D(0),:)
ztrds(A2D(0),:) = 0._wp
CALL trd_tra_mng( trdt, ztrds, ktrd, kt, Kmm )
!!gm Gurvan, verify the jptra_evd trend please !
CASE( jptra_evd ) ; avt_evd(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
CASE( jptra_evd ) ; avt_evd(A2D(0),:) = ptrd(A2D(0),:) * tmask(A2D(0),:)
CASE DEFAULT ! other trends: masked trends
trdt(:,:,:) = ptrd(:,:,:) * tmask(:,:,:) ! mask & store
trdt(A2D(0),:) = ptrd(A2D(0),:) * tmask(A2D(0),:) ! mask & store
END SELECT
!
ENDIF
......@@ -127,44 +128,44 @@ CONTAINS
CALL trd_tra_mng( trdt , ztrds, ktrd, kt, Kmm )
CASE( jptra_zdfp ) ! diagnose the "PURE" Kz trend (here: just before the swap)
! ! iso-neutral diffusion case otherwise jptra_zdf is "PURE"
ALLOCATE( zwt(jpi,jpj,jpk), zws(jpi,jpj,jpk), ztrdt(jpi,jpj,jpk) )
ALLOCATE( zwt(A2D(0),jpk), zws(A2D(0),jpk), ztrdt(A2D(0),jpk) )
!
zwt(:,:, 1 ) = 0._wp ; zws(:,:, 1 ) = 0._wp ! vertical diffusive fluxes
zwt(:,:,jpk) = 0._wp ; zws(:,:,jpk) = 0._wp
DO jk = 2, jpk
zwt(:,:,jk) = avt(:,:,jk) * ( ts(:,:,jk-1,jp_tem,Krhs) - ts(:,:,jk,jp_tem,Krhs) ) &
& / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
zws(:,:,jk) = avs(:,:,jk) * ( ts(:,:,jk-1,jp_sal,Krhs) - ts(:,:,jk,jp_sal,Krhs) ) &
& / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
END DO
zwt(A2D(0), 1 ) = 0._wp ; zws(A2D(0), 1 ) = 0._wp ! vertical diffusive fluxes
zwt(A2D(0),jpk) = 0._wp ; zws(A2D(0),jpk) = 0._wp
DO_3D( 0, 0, 0, 0, 2, jpk )
z1e3w = 1._wp / ( e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) )
zwt(ji,jj,jk) = avt(ji,jj,jk) * ( ts(ji,jj,jk-1,jp_tem,Krhs) - ts(ji,jj,jk,jp_tem,Krhs) ) * z1e3w
zws(ji,jj,jk) = avs(ji,jj,jk) * ( ts(ji,jj,jk-1,jp_sal,Krhs) - ts(ji,jj,jk,jp_sal,Krhs) ) * z1e3w
END_3D
!
ztrdt(:,:,jpk) = 0._wp ; ztrds(:,:,jpk) = 0._wp
DO jk = 1, jpkm1
ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / e3t(:,:,jk,Kmm)
ztrds(:,:,jk) = ( zws(:,:,jk) - zws(:,:,jk+1) ) / e3t(:,:,jk,Kmm)
END DO
ztrdt(A2D(0),jpk) = 0._wp ; ztrds(A2D(0),jpk) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z1e3w = 1._wp / e3w(ji,jj,jk,Kmm)
ztrdt(ji,jj,jk) = ( zwt(ji,jj,jk) - zwt(ji,jj,jk+1) ) * z1e3w
ztrds(ji,jj,jk) = ( zws(ji,jj,jk) - zws(ji,jj,jk+1) ) * z1e3w
END_3D
CALL trd_tra_mng( ztrdt, ztrds, jptra_zdfp, kt, Kmm )
!
! ! Also calculate EVD trend at this point.
zwt(:,:,:) = 0._wp ; zws(:,:,:) = 0._wp ! vertical diffusive fluxes
DO jk = 2, jpk
zwt(:,:,jk) = avt_evd(:,:,jk) * ( ts(:,:,jk-1,jp_tem,Krhs) - ts(:,:,jk,jp_tem,Krhs) ) &
& / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
zws(:,:,jk) = avt_evd(:,:,jk) * ( ts(:,:,jk-1,jp_sal,Krhs) - ts(:,:,jk,jp_sal,Krhs) ) &
& / e3w(:,:,jk,Kmm) * tmask(:,:,jk)
END DO
zwt(A2D(0), :) = 0._wp ; zws(A2D(0), :) = 0._wp ! vertical diffusive fluxes
DO_3D( 0, 0, 0, 0, 2, jpk )
z1e3w = 1._wp / ( e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) )
zwt(ji,jj,jk) = avt_evd(ji,jj,jk) * ( ts(ji,jj,jk-1,jp_tem,Krhs) - ts(ji,jj,jk,jp_tem,Krhs) ) * z1e3w
zws(ji,jj,jk) = avt_evd(ji,jj,jk) * ( ts(ji,jj,jk-1,jp_sal,Krhs) - ts(ji,jj,jk,jp_sal,Krhs) ) * z1e3w
END_3D
!
ztrdt(:,:,jpk) = 0._wp ; ztrds(:,:,jpk) = 0._wp
DO jk = 1, jpkm1
ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / e3t(:,:,jk,Kmm)
ztrds(:,:,jk) = ( zws(:,:,jk) - zws(:,:,jk+1) ) / e3t(:,:,jk,Kmm)
END DO
ztrdt(A2D(0),jpk) = 0._wp ; ztrds(A2D(0),jpk) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z1e3w = 1._wp / e3w(ji,jj,jk,Kmm)
ztrdt(ji,jj,jk) = ( zwt(ji,jj,jk) - zwt(ji,jj,jk+1) ) * z1e3w
ztrds(ji,jj,jk) = ( zws(ji,jj,jk) - zws(ji,jj,jk+1) ) * z1e3w
END_3D
CALL trd_tra_mng( ztrdt, ztrds, jptra_evd, kt, Kmm )
!
DEALLOCATE( zwt, zws, ztrdt )
!
CASE DEFAULT ! other trends: mask and send T & S trends to trd_tra_mng
ztrds(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
ztrds(A2D(0),:) = ptrd(A2D(0),:) * tmask(A2D(0),:)
CALL trd_tra_mng( trdt, ztrds, ktrd, kt, Kmm )
END SELECT
ENDIF
......@@ -177,7 +178,7 @@ CONTAINS
CASE( jptra_yad ) ; CALL trd_tra_adv( ptrd , pu , ptra, 'Y', ztrds, Kmm )
CASE( jptra_zad ) ; CALL trd_tra_adv( ptrd , pu , ptra, 'Z', ztrds, Kmm )
CASE DEFAULT ! other trends: just masked
ztrds(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
ztrds(A2D(0),:) = ptrd(A2D(0),:) * tmask(A2D(0),:)
END SELECT
! ! send trend to trd_trc
CALL trd_trc( ztrds, ktra, ktrd, kt, Kmm )
......@@ -199,12 +200,12 @@ CONTAINS
!! k-advective trends = -un. di+1[T] = -( dk+1[fi] - tn dk+1[un] )
!! where fi is the incoming advective flux.
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pf ! advective flux in one direction
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pu ! now velocity in one direction
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pt ! now or before tracer
CHARACTER(len=1) , INTENT(in ) :: cdir ! X/Y/Z direction
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: ptrd ! advective trend in one direction
INTEGER, INTENT(in) :: Kmm ! time level index
REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pf ! advective flux in one direction
REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pu ! now velocity in one direction
REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pt ! now or before tracer
CHARACTER(len=1) , INTENT(in ) :: cdir ! X/Y/Z direction
REAL(wp), DIMENSION(:,:,:), INTENT( out) :: ptrd ! advective trend in one direction
INTEGER, INTENT(in) :: Kmm ! time level index
!
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ii, ij, ik ! index shift as function of the direction
......@@ -217,9 +218,7 @@ CONTAINS
END SELECT
!
! ! set to zero uncomputed values
ptrd(jpi,:,:) = 0._wp ; ptrd(1,:,:) = 0._wp
ptrd(:,jpj,:) = 0._wp ; ptrd(:,1,:) = 0._wp
ptrd(:,:,jpk) = 0._wp
ptrd(:,:,:) = 0._wp
!
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! advective trend
ptrd(ji,jj,jk) = - ( pf (ji,jj,jk) - pf (ji-ii,jj-ij,jk-ik) &
......@@ -248,7 +247,7 @@ CONTAINS
! ! 3D output of tracers trends using IOM interface
IF( ln_tra_trd ) CALL trd_tra_iom ( ptrdx, ptrdy, ktrd, kt, Kmm )
! ! Integral Constraints Properties for tracers trends !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! ! Integral Constraints Properties for tracers trends
IF( ln_glo_trd ) CALL trd_glo( ptrdx, ptrdy, ktrd, 'TRA', kt, Kmm )
! ! Potential ENergy trends
......@@ -305,8 +304,8 @@ CONTAINS
INTEGER , INTENT(in ) :: kt ! time step
INTEGER , INTENT(in ) :: Kmm ! time level index
!!
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ikbu, ikbv ! local integers
INTEGER :: ji, jj
REAL(wp) :: z1e3
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z2dx, z2dy ! 2D workspace
!!----------------------------------------------------------------------
!
......@@ -329,9 +328,12 @@ CONTAINS
CASE( jptra_zad ) ; CALL iom_put( "ttrd_zad" , ptrdx ) ! z- vertical advection
CALL iom_put( "strd_zad" , ptrdy )
IF( ln_linssh ) THEN ! cst volume : adv flux through z=0 surface
ALLOCATE( z2dx(jpi,jpj), z2dy(jpi,jpj) )
z2dx(:,:) = ww(:,:,1) * ts(:,:,1,jp_tem,Kmm) / e3t(:,:,1,Kmm)
z2dy(:,:) = ww(:,:,1) * ts(:,:,1,jp_sal,Kmm) / e3t(:,:,1,Kmm)
ALLOCATE( z2dx(A2D(0)), z2dy(A2D(0)) )
DO_2D( 0, 0, 0, 0 )
z1e3 = ww(ji,jj,1) / e3t(ji,jj,1,Kmm)
z2dx(ji,jj) = ts(ji,jj,1,jp_tem,Kmm) * z1e3
z2dy(ji,jj) = ts(ji,jj,1,jp_sal,Kmm) * z1e3
END_2D
CALL iom_put( "ttrd_sad", z2dx )
CALL iom_put( "strd_sad", z2dy )
DEALLOCATE( z2dx, z2dy )
......
src/OCE/TRD/trdvor.F90
View file @ 4dbfc9bb
......@@ -68,9 +68,9 @@ CONTAINS
!!----------------------------------------------------------------------------
!! *** ROUTINE trd_vor_alloc ***
!!----------------------------------------------------------------------------
ALLOCATE( vor_avr (jpi,jpj) , vor_avrb(jpi,jpj) , vor_avrbb (jpi,jpj) , &
& vor_avrbn (jpi,jpj) , rotot (jpi,jpj) , vor_avrtot(jpi,jpj) , &
& vor_avrres(jpi,jpj) , vortrd (jpi,jpj,jpltot_vor) , &
ALLOCATE( vor_avr (A2D(0)) , vor_avrb(A2D(0)) , vor_avrbb (A2D(0)) , &
& vor_avrbn (A2D(0)) , rotot (A2D(0)) , vor_avrtot(A2D(0)) , &
& vor_avrres(A2D(0)) , vortrd (A2D(0),jpltot_vor) , &
& ndexvor1 (jpi*jpj) , STAT= trd_vor_alloc )
!
CALL mpp_sum ( 'trdvor', trd_vor_alloc )
......@@ -105,9 +105,9 @@ CONTAINS
CASE( jpdyn_zad ) ; CALL trd_vor_zint( putrd, pvtrd, jpvor_zad, Kmm ) ! Vertical Advection
CASE( jpdyn_spg ) ; CALL trd_vor_zint( putrd, pvtrd, jpvor_spg, Kmm ) ! Surface Pressure Grad.
CASE( jpdyn_zdf ) ! Vertical Diffusion
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! wind stress trends
ztswu(ji,jj) = 0.5 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( e3u(ji,jj,1,Kmm) * rho0 )
ztswv(ji,jj) = 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( e3v(ji,jj,1,Kmm) * rho0 )
DO_2D( 0, 1, 0, 1 ) ! wind stress trends
ztswu(ji,jj) = 0.5 * ( utau_b(ji,jj) + utauU(ji,jj) ) / ( e3u(ji,jj,1,Kmm) * rho0 )
ztswv(ji,jj) = 0.5 * ( vtau_b(ji,jj) + vtauV(ji,jj) ) / ( e3v(ji,jj,1,Kmm) * rho0 )
END_2D
CALL trd_vor_zint( putrd, pvtrd, jpvor_zdf, Kmm ) ! zdf trend including surf./bot. stresses
CALL trd_vor_zint( ztswu, ztswv, jpvor_swf, Kmm ) ! surface wind stress
......@@ -165,7 +165,7 @@ CONTAINS
SELECT CASE( ktrd )
!
CASE( jpvor_bfr ) ! bottom friction
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
DO_2D( 0, 1, 0, 1 )
ikbu = mbkv(ji,jj)
ikbv = mbkv(ji,jj)
zudpvor(ji,jj) = putrdvor(ji,jj) * e3u(ji,jj,ikbu,Kmm) * e1u(ji,jj) * umask(ji,jj,ikbu)
......@@ -173,14 +173,18 @@ CONTAINS
END_2D
!
CASE( jpvor_swf ) ! wind stress
zudpvor(:,:) = putrdvor(:,:) * e3u(:,:,1,Kmm) * e1u(:,:) * umask(:,:,1)
zvdpvor(:,:) = pvtrdvor(:,:) * e3v(:,:,1,Kmm) * e2v(:,:) * vmask(:,:,1)
DO_2D( 0, 1, 0, 1 )
zudpvor(ji,jj) = putrdvor(ji,jj) * e3u(ji,jj,1,Kmm) * e1u(ji,jj) * umask(ji,jj,1)
zvdpvor(ji,jj) = pvtrdvor(ji,jj) * e3v(ji,jj,1,Kmm) * e2v(ji,jj) * vmask(ji,jj,1)
END_2D
!
END SELECT
! Average except for Beta.V
zudpvor(:,:) = zudpvor(:,:) * r1_hu(:,:,Kmm)
zvdpvor(:,:) = zvdpvor(:,:) * r1_hv(:,:,Kmm)
DO_2D( 0, 1, 0, 1 )
zudpvor(ji,jj) = zudpvor(ji,jj) * r1_hu(ji,jj,Kmm)
zvdpvor(ji,jj) = zvdpvor(ji,jj) * r1_hv(ji,jj,Kmm)
END_2D
! Curl
DO_2D( 0, 0, 0, 0 )
......@@ -238,10 +242,11 @@ CONTAINS
! I vertical integration of 3D trends
! =====================================
! putrdvor and pvtrdvor terms
DO jk = 1,jpk
zudpvor(:,:) = zudpvor(:,:) + putrdvor(:,:,jk) * e3u(:,:,jk,Kmm) * e1u(:,:) * umask(:,:,jk)
zvdpvor(:,:) = zvdpvor(:,:) + pvtrdvor(:,:,jk) * e3v(:,:,jk,Kmm) * e2v(:,:) * vmask(:,:,jk)
END DO
zudpvor(:,:) = 0._wp ; zvdpvor(:,:) = 0._wp
DO_3D( 0, 1, 0, 1, 1, jpk )
zudpvor(ji,jj) = zudpvor(ji,jj) + putrdvor(ji,jj,jk) * e3u(ji,jj,jk,Kmm) * e1u(ji,jj) * umask(ji,jj,jk)
zvdpvor(ji,jj) = zvdpvor(ji,jj) + pvtrdvor(ji,jj,jk) * e3v(ji,jj,jk,Kmm) * e2v(ji,jj) * vmask(ji,jj,jk)
END_3D
! Planetary vorticity: 2nd computation (Beta.V term) store the vertical sum
! as Beta.V term need intergration, not average
......@@ -254,8 +259,10 @@ CONTAINS
ENDIF
!
! Average
zudpvor(:,:) = zudpvor(:,:) * r1_hu(:,:,Kmm)
zvdpvor(:,:) = zvdpvor(:,:) * r1_hv(:,:,Kmm)
DO_2D( 0, 1, 0, 1 )
zudpvor(ji,jj) = zudpvor(ji,jj) * r1_hu(ji,jj,Kmm)
zvdpvor(ji,jj) = zvdpvor(ji,jj) * r1_hv(ji,jj,Kmm)
END_2D
!
! Curl
DO_2D( 0, 0, 0, 0 )
......@@ -368,10 +375,6 @@ CONTAINS
zmean = 1._wp / REAL( nmoydpvor, wp )
vor_avrres(:,:) = vor_avrtot(:,:) - rotot(:,:) / zmean
! Boundary conditions
CALL lbc_lnk( 'trdvor', vor_avrtot, 'F', 1.0_wp , vor_avrres, 'F', 1.0_wp )
! III.3 time evolution array swap
! ------------------------------
vor_avrbb(:,:) = vor_avrb(:,:)
......
src/OCE/USR/usrdef_fmask.F90
View file @ 4dbfc9bb
......@@ -69,25 +69,25 @@ CONTAINS
IF(lwp) WRITE(numout,*) ' Gibraltar '
ij0 = 101 + nn_hls ; ij1 = 101 + nn_hls ! Gibraltar strait : partial slip (pfmsk=0.5)
ii0 = 139 + nn_hls - 1 ; ii1 = 140 + nn_hls - 1
pfmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp
pfmsk( mi0(ii0,nn_hls):mi1(ii1,nn_hls) , mj0(ij0,nn_hls):mj1(ij1,nn_hls) , 1:jpk ) = 0.5_wp
ij0 = 102 + nn_hls ; ij1 = 102 + nn_hls
ii0 = 139 + nn_hls - 1 ; ii1 = 140 + nn_hls - 1
pfmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp
pfmsk( mi0(ii0,nn_hls):mi1(ii1,nn_hls) , mj0(ij0,nn_hls):mj1(ij1,nn_hls) , 1:jpk ) = 0.5_wp
!
IF(lwp) WRITE(numout,*) ' Bab el Mandeb '
ij0 = 87 + nn_hls ; ij1 = 88 + nn_hls ! Bab el Mandeb : partial slip (pfmsk=1)
ii0 = 160 + nn_hls - 1 ; ii1 = 160 + nn_hls - 1
pfmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp
pfmsk( mi0(ii0,nn_hls):mi1(ii1,nn_hls) , mj0(ij0,nn_hls):mj1(ij1,nn_hls) , 1:jpk ) = 1._wp
ij0 = 88 + nn_hls ; ij1 = 88 + nn_hls
ii0 = 159 + nn_hls - 1 ; ii1 = 159 + nn_hls - 1
pfmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp
pfmsk( mi0(ii0,nn_hls):mi1(ii1,nn_hls) , mj0(ij0,nn_hls):mj1(ij1,nn_hls) , 1:jpk ) = 1._wp
!
! We keep this as an example but it is instable in this case
!IF(lwp) WRITE(numout,*) ' Danish straits '
! ij0 = 115 ; ij1 = 115 ! Danish straits : strong slip (pfmsk > 2)
! ii0 = 145 ; ii1 = 146 ; pfmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
! ii0 = 145 ; ii1 = 146 ; pfmsk( mi0(ii0,nn_hls):mi1(ii1,nn_hls) , mj0(ij0,nn_hls):mj1(ij1,nn_hls) , 1:jpk ) = 4._wp
! ij0 = 116 ; ij1 = 116
! ii0 = 145 ; ii1 = 146 ; pfmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
! ii0 = 145 ; ii1 = 146 ; pfmsk( mi0(ii0,nn_hls):mi1(ii1,nn_hls) , mj0(ij0,nn_hls):mj1(ij1,nn_hls) , 1:jpk ) = 4._wp
!
CASE( 1 ) ! R1 case
IF(lwp) WRITE(numout,*)
......@@ -104,42 +104,42 @@ CONTAINS
IF(lwp) WRITE(numout,*) ' Gibraltar '
ii0 = 282 + nn_hls - 1 ; ii1 = 283 + nn_hls - 1 ! Gibraltar Strait
ij0 = 241 + nn_hls - isrow ; ij1 = 241 + nn_hls - isrow
pfmsk( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
pfmsk( mi0(ii0,nn_hls):mi1(ii1,nn_hls),mj0(ij0,nn_hls):mj1(ij1,nn_hls),1:jpk ) = 2._wp
!
IF(lwp) WRITE(numout,*) ' Bhosporus '
ii0 = 314 + nn_hls - 1 ; ii1 = 315 + nn_hls - 1 ! Bhosporus Strait
ij0 = 248 + nn_hls - isrow ; ij1 = 248 + nn_hls - isrow
pfmsk( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
pfmsk( mi0(ii0,nn_hls):mi1(ii1,nn_hls),mj0(ij0,nn_hls):mj1(ij1,nn_hls),1:jpk ) = 2._wp
!
IF(lwp) WRITE(numout,*) ' Makassar (Top) '
ii0 = 48 + nn_hls - 1 ; ii1 = 48 + nn_hls - 1 ! Makassar Strait (Top)
ij0 = 189 + nn_hls - isrow ; ij1 = 190 + nn_hls - isrow
pfmsk( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp
pfmsk( mi0(ii0,nn_hls):mi1(ii1,nn_hls),mj0(ij0,nn_hls):mj1(ij1,nn_hls),1:jpk ) = 3._wp
!
IF(lwp) WRITE(numout,*) ' Lombok '
ii0 = 44 + nn_hls - 1 ; ii1 = 44 + nn_hls - 1 ! Lombok Strait
ij0 = 164 + nn_hls - isrow ; ij1 = 165 + nn_hls - isrow
pfmsk( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
pfmsk( mi0(ii0,nn_hls):mi1(ii1,nn_hls),mj0(ij0,nn_hls):mj1(ij1,nn_hls),1:jpk ) = 2._wp
!
IF(lwp) WRITE(numout,*) ' Ombai '
ii0 = 53 + nn_hls - 1 ; ii1 = 53 + nn_hls - 1 ! Ombai Strait
ij0 = 164 + nn_hls - isrow ; ij1 = 165 + nn_hls - isrow
pfmsk( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
pfmsk( mi0(ii0,nn_hls):mi1(ii1,nn_hls),mj0(ij0,nn_hls):mj1(ij1,nn_hls),1:jpk ) = 2._wp
!
IF(lwp) WRITE(numout,*) ' Timor Passage '
ii0 = 56 + nn_hls - 1 ; ii1 = 56 + nn_hls - 1 ! Timor Passage
ij0 = 164 + nn_hls - isrow ; ij1 = 165 + nn_hls - isrow
pfmsk( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp
pfmsk( mi0(ii0,nn_hls):mi1(ii1,nn_hls),mj0(ij0,nn_hls):mj1(ij1,nn_hls),1:jpk ) = 2._wp
!
IF(lwp) WRITE(numout,*) ' West Halmahera '
ii0 = 58 + nn_hls - 1 ; ii1 = 58 + nn_hls - 1 ! West Halmahera Strait
ij0 = 181 + nn_hls - isrow ; ij1 = 182 + nn_hls - isrow
pfmsk( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp
pfmsk( mi0(ii0,nn_hls):mi1(ii1,nn_hls),mj0(ij0,nn_hls):mj1(ij1,nn_hls),1:jpk ) = 3._wp
!
IF(lwp) WRITE(numout,*) ' East Halmahera '
ii0 = 55 + nn_hls - 1 ; ii1 = 55 + nn_hls - 1 ! East Halmahera Strait
ij0 = 181 + nn_hls - isrow ; ij1 = 182 + nn_hls - isrow
pfmsk( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp
pfmsk( mi0(ii0,nn_hls):mi1(ii1,nn_hls),mj0(ij0,nn_hls):mj1(ij1,nn_hls),1:jpk ) = 3._wp
!
CASE DEFAULT
IF(lwp) WRITE(numout,*)
......
src/OCE/USR/usrdef_hgr.F90
View file @ 4dbfc9bb
......@@ -115,8 +115,8 @@ CONTAINS
ENDIF
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zim1 = REAL( mig0(ji), wp ) - 1. ; zim05 = REAL( mig0(ji), wp ) - 1.5
zjm1 = REAL( mjg0(jj), wp ) - 1. ; zjm05 = REAL( mjg0(jj), wp ) - 1.5
zim1 = REAL( mig(ji,0), wp ) - 1. ; zim05 = REAL( mig(ji,0), wp ) - 1.5
zjm1 = REAL( mjg(jj,0), wp ) - 1. ; zjm05 = REAL( mjg(jj,0), wp ) - 1.5
!
!glamt(i,j) longitude at T-point
!gphit(i,j) latitude at T-point
......
src/OCE/USR/usrdef_sbc.F90
View file @ 4dbfc9bb
......@@ -109,7 +109,7 @@ CONTAINS
ztrp= - 40.e0 ! retroaction term on heat fluxes (W/m2/K)
zconv = 3.16e-5 ! convertion factor: 1 m/yr => 3.16e-5 mm/s
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! emp and rnf used in sshwzv over the whole domain
DO_2D( 0, 0, 0, 0 )
! domain from 15 deg to 50 deg between 27 and 28 degC at 15N, -3
! and 13 degC at 50N 53.5 + or - 11 = 1/4 period :
! 64.5 in summer, 42.5 in winter
......@@ -119,6 +119,8 @@ CONTAINS
! 23.5 deg : tropics
qsr (ji,jj) = 230 * COS( 3.1415 * ( gphit(ji,jj) - 23.5 * zcos_sais1 ) / ( 0.9 * 180 ) )
qns (ji,jj) = ztrp * ( ts(ji,jj,1,jp_tem,Kbb) - t_star ) - qsr(ji,jj)
END_2D
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! emp and rnf used in sshwzv over the whole domain
IF( gphit(ji,jj) >= 14.845 .AND. 37.2 >= gphit(ji,jj) ) THEN ! zero at 37.8 deg, max at 24.6 deg
emp (ji,jj) = zemp_S * zconv &
& * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (24.6 - 37.2) ) &
......@@ -137,6 +139,8 @@ CONTAINS
! freshwater (mass flux) and update of qns with heat content of emp
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! emp used in sshwzv over the whole domain
emp (ji,jj) = emp(ji,jj) - zsumemp * tmask(ji,jj,1) ! freshwater flux (=0 in domain average)
END_2D
DO_2D( 0, 0, 0, 0 )
sfx (ji,jj) = 0.0_wp ! no salt flux
qns (ji,jj) = qns(ji,jj) - emp(ji,jj) * sst_m(ji,jj) * rcp ! evap and precip are at SST
END_2D
......@@ -166,19 +170,17 @@ CONTAINS
! seasonal oscillation intensity
ztau_sais = 0.015
ztaun = ztau - ztau_sais * COS( (ztime - ztimemax) / (ztimemin - ztimemax) * rpi )
DO_2D( 1, 1, 1, 1 )
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
! domain from 15deg to 50deg and 1/2 period along 14deg
! so 5/4 of half period with seasonal cycle
utau(ji,jj) = - ztaun * SIN( rpi * (gphiu(ji,jj) - 15.) / (29.-15.) )
vtau(ji,jj) = ztaun * SIN( rpi * (gphiv(ji,jj) - 15.) / (29.-15.) )
utau(ji,jj) = - ztaun * SIN( rpi * (gphit(ji,jj) - 15.) / (29.-15.) )
vtau(ji,jj) = ztaun * SIN( rpi * (gphit(ji,jj) - 15.) / (29.-15.) )
END_2D
! module of wind stress and wind speed at T-point
zcoef = 1. / ( zrhoa * zcdrag )
DO_2D( 0, 0, 0, 0 )
ztx = utau(ji-1,jj ) + utau(ji,jj)
zty = vtau(ji ,jj-1) + vtau(ji,jj)
zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
zmod = SQRT( utau(ji,jj) * utau(ji,jj) + vtau(ji,jj) * vtau(ji,jj) )
taum(ji,jj) = zmod
wndm(ji,jj) = SQRT( zmod * zcoef )
END_2D
......
src/OCE/ZDF/zdf_oce.F90
View file @ 4dbfc9bb
......@@ -53,6 +53,8 @@ MODULE zdf_oce
REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) :: avmb , avtb !: background profile of avm and avt [m2/s]
REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: avtb_2d !: horizontal shape of background Kz profile [-]
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: zdf_oce.F90 14072 2020-12-04 07:48:38Z laurent $
......@@ -65,9 +67,9 @@ CONTAINS
!! *** FUNCTION zdf_oce_alloc ***
!!----------------------------------------------------------------------
!
ALLOCATE( avm (jpi,jpj,jpk) , avm_k(jpi,jpj,jpk) , avs(jpi,jpj,jpk) , &
& avt (jpi,jpj,jpk) , avt_k(jpi,jpj,jpk) , en (jpi,jpj,jpk) , &
& avmb(jpk) , avtb(jpk) , avtb_2d(jpi,jpj) , STAT = zdf_oce_alloc )
ALLOCATE( avm (jpi,jpj,jpk) , avm_k(jpi,jpj,jpk) , avs(A2D(0),jpk) , &
& avt (A2D(0) ,jpk) , avt_k(A2D(0) ,jpk) , en (A2D(0),jpk) , &
& avmb(jpk) , avtb(jpk) , avtb_2d(A2D(0)) , STAT = zdf_oce_alloc )
!
IF( zdf_oce_alloc /= 0 ) CALL ctl_stop( 'STOP', 'zdf_oce_alloc: failed to allocate arrays' )
!
......
src/OCE/ZDF/zdfddm.F90
View file @ 4dbfc9bb
......@@ -70,9 +70,9 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER, INTENT(in ) :: kt ! ocean time-step index
INTEGER, INTENT(in ) :: Kmm ! ocean time level index
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm ! Kz on momentum (w-points)
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avt ! Kz on temperature (w-points)
REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_avs ! Kz on salinity (w-points)
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm ! Kz on momentum (w-points)
REAL(wp), DIMENSION(A2D(0),jpk), INTENT(inout) :: p_avt ! Kz on temperature (w-points)
REAL(wp), DIMENSION(A2D(0),jpk), INTENT( out) :: p_avs ! Kz on salinity (w-points)
!
INTEGER :: ji, jj , jk ! dummy loop indices
REAL(wp) :: zaw, zbw, zrw ! local scalars
......@@ -82,7 +82,7 @@ CONTAINS
REAL(wp) :: zavft ! - -
REAL(dp) :: zavfs ! - -
REAL(wp) :: zavdt, zavds ! - -
REAL(wp), DIMENSION(A2D(nn_hls)) :: zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3
REAL(wp), DIMENSION(A2D(0)) :: zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3
!!----------------------------------------------------------------------
!
! ! ===============
......@@ -94,7 +94,7 @@ CONTAINS
!!gm ==>>> test in the loop instead of use of mask arrays
!!gm and many acces in memory
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !== R=zrau = (alpha / beta) (dk[t] / dk[s]) ==!
DO_2D( 0, 0, 0, 0 ) !== R=zrau = (alpha / beta) (dk[t] / dk[s]) ==!
zrw = ( gdepw(ji,jj,jk ,Kmm) - gdept(ji,jj,jk,Kmm) ) &
!!gm please, use e3w at Kmm below
& / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) )
......@@ -110,7 +110,7 @@ CONTAINS
zrau(ji,jj) = MAX( 1.e-20, zdt / zds ) ! only retains positive value of zrau
END_2D
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !== indicators ==!
DO_2D( 0, 0, 0, 0 ) !== indicators ==!
! stability indicator: msks=1 if rn2>0; 0 elsewhere
IF( rn2(ji,jj,jk) + 1.e-12 <= 0. ) THEN ; zmsks(ji,jj) = 0._wp
ELSE ; zmsks(ji,jj) = 1._wp * wmask(ji,jj,jk) ! mask so avt and avs masked
......@@ -137,7 +137,7 @@ CONTAINS
! Update avt and avs
! ------------------
! Constant eddy coefficient: reset to the background value
DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
DO_2D_OVR( 0, 0, 0, 0 )
zinr = 1._wp / zrau(ji,jj)
! salt fingering
zrr = zrau(ji,jj) / rn_hsbfr
......
src/OCE/ZDF/zdfevd.F90
View file @ 4dbfc9bb
......@@ -56,9 +56,10 @@ CONTAINS
!!
!! ** Action : avt, avm enhanced where static instability occurs
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step indexocean time step
INTEGER , INTENT(in ) :: Kmm, Krhs ! time level indices
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points)
INTEGER , INTENT(in ) :: kt ! ocean time-step indexocean time step
INTEGER , INTENT(in ) :: Kmm, Krhs ! time level indices
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: p_avm ! momentum Kz (w-points)
REAL(wp), DIMENSION(A2D(0),jpk), INTENT(inout) :: p_avt ! tracer Kz (w-points)
!
INTEGER :: ji, jj, jk ! dummy loop indices
! NOTE: [tiling] use a SAVE array to store diagnostics, then send after all tiles are finished. This is necessary because p_avt/p_avm are modified on adjacent tiles when using nn_hls > 1. zavt_evd/zavm_evd are then zero on some points when subsequently calculated for these tiles.
......@@ -73,12 +74,12 @@ CONTAINS
IF(lwp) WRITE(numout,*)
ENDIF
ALLOCATE( zavt_evd(jpi,jpj,jpk) )
IF( nn_evdm == 1 ) ALLOCATE( zavm_evd(jpi,jpj,jpk) )
ALLOCATE( zavt_evd(A2D(0),jpk) )
IF( nn_evdm == 1 ) ALLOCATE( zavm_evd(A2D(0),jpk) )
ENDIF
!
!
DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
DO_3D_OVR( 0, 0, 0, 0, 1, jpk )
zavt_evd(ji,jj,jk) = p_avt(ji,jj,jk) ! set avt prior to evd application
END_3D
!
......@@ -86,7 +87,7 @@ CONTAINS
!
CASE ( 1 ) !== enhance tracer & momentum Kz ==! (if rn2<-1.e-12)
!
DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
DO_3D_OVR( 0, 0, 0, 0, 1, jpk )
zavm_evd(ji,jj,jk) = p_avm(ji,jj,jk) ! set avm prior to evd application
END_3D
!
......@@ -96,14 +97,14 @@ CONTAINS
! p_avm(2:jpi,2:jpj,2:jpkm1) = rn_evd * wmask(2:jpi,2:jpj,2:jpkm1)
! END WHERE
!
DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
DO_3D_OVR( 0, 0, 0, 0, 1, jpkm1 )
IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) THEN
p_avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk)
p_avm(ji,jj,jk) = rn_evd * wmask(ji,jj,jk)
ENDIF
END_3D
!
DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
DO_3D_OVR( 0, 0, 0, 0, 1, jpk )
zavm_evd(ji,jj,jk) = p_avm(ji,jj,jk) - zavm_evd(ji,jj,jk) ! change in avm due to evd
END_3D
IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile
......@@ -117,14 +118,14 @@ CONTAINS
! p_avt(2:jpi,2:jpj,2:jpkm1) = rn_evd * wmask(2:jpi,2:jpj,2:jpkm1)
! END WHERE
DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
DO_3D_OVR( 0, 0, 0, 0, 1, jpkm1 )
IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) &
p_avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk)
END_3D
!
END SELECT
!
DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
DO_3D_OVR( 0, 0, 0, 0, 1, jpk )
zavt_evd(ji,jj,jk) = p_avt(ji,jj,jk) - zavt_evd(ji,jj,jk) ! change in avt due to evd
END_3D
!
......