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 1536 additions and 760 deletions
src/OCE/DIA/diawri.F90
View file @ 3efbcb17
......@@ -122,8 +122,9 @@ CONTAINS
REAL(wp):: zztmp , zztmpx ! local scalar
REAL(wp):: zztmp2, zztmpy ! - -
REAL(wp):: ze3
REAL(wp), DIMENSION(A2D( 0)) :: z2d ! 2D workspace
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: z3d ! 3D workspace
REAL(wp), DIMENSION(A2D(0)) :: z2d0 ! 2D workspace
REAL(wp), DIMENSION(A2D(0) ,jpk) :: z3d0 ! 3D workspace
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: z3d ! 3D workspace
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('dia_wri')
......@@ -135,8 +136,9 @@ CONTAINS
ENDIF
! initialize arrays
z2d(:,:) = 0._wp
z3d(:,:,:) = 0._wp
z2d0(:,:) = 0._wp
z3d0(:,:,:) = 0._wp
z3d (:,:,:) = 0._wp
! Output of initial vertical scale factor
CALL iom_put("e3t_0", e3t_0(:,:,:) )
......@@ -146,44 +148,44 @@ CONTAINS
!
IF ( iom_use("tpt_dep") ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = gdept(ji,jj,jk,Kmm)
z3d0(ji,jj,jk) = gdept(ji,jj,jk,Kmm)
END_3D
CALL iom_put( "tpt_dep", z3d )
CALL iom_put( "tpt_dep", z3d0 )
ENDIF
! --- vertical scale factors --- !
IF ( iom_use("e3t") .OR. iom_use("e3tdef") ) THEN ! time-varying e3t
DO_3D( 0, 0, 0, 0, 1, jpk )
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
z3d(ji,jj,jk) = e3t(ji,jj,jk,Kmm)
END_3D
CALL iom_put( "e3t", z3d )
IF ( iom_use("e3tdef") ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
z3d(ji,jj,jk) = ( ( z3d(ji,jj,jk) - e3t_0(ji,jj,jk) ) / e3t_0(ji,jj,jk) * 100._wp * tmask(ji,jj,jk) ) ** 2
END_3D
CALL iom_put( "e3tdef", z3d )
ENDIF
ENDIF
IF ( iom_use("e3u") ) THEN ! time-varying e3u
DO_3D( 0, 0, 0, 0, 1, jpk )
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
z3d(ji,jj,jk) = e3u(ji,jj,jk,Kmm)
END_3D
CALL iom_put( "e3u" , z3d )
ENDIF
IF ( iom_use("e3v") ) THEN ! time-varying e3v
DO_3D( 0, 0, 0, 0, 1, jpk )
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
z3d(ji,jj,jk) = e3v(ji,jj,jk,Kmm)
END_3D
CALL iom_put( "e3v" , z3d )
ENDIF
IF ( iom_use("e3w") ) THEN ! time-varying e3w
DO_3D( 0, 0, 0, 0, 1, jpk )
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
z3d(ji,jj,jk) = e3w(ji,jj,jk,Kmm)
END_3D
CALL iom_put( "e3w" , z3d )
ENDIF
IF ( iom_use("e3f") ) THEN ! time-varying e3f caution here at Kaa
DO_3D( 0, 0, 0, 0, 1, jpk )
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
z3d(ji,jj,jk) = e3f(ji,jj,jk)
END_3D
CALL iom_put( "e3f" , z3d )
......@@ -213,9 +215,9 @@ CONTAINS
IF ( iom_use("sbt") ) THEN
DO_2D( 0, 0, 0, 0 )
ikbot = mbkt(ji,jj)
z2d(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm)
z2d0(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm)
END_2D
CALL iom_put( "sbt", z2d ) ! bottom temperature
CALL iom_put( "sbt", z2d0 ) ! bottom temperature
ENDIF
CALL iom_put( "soce", ts(:,:,:,jp_sal,Kmm) ) ! 3D salinity
......@@ -223,9 +225,9 @@ CONTAINS
IF ( iom_use("sbs") ) THEN
DO_2D( 0, 0, 0, 0 )
ikbot = mbkt(ji,jj)
z2d(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm)
z2d0(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm)
END_2D
CALL iom_put( "sbs", z2d ) ! bottom salinity
CALL iom_put( "sbs", z2d0 ) ! bottom salinity
ENDIF
IF( .NOT.lk_SWE ) CALL iom_put( "rhop", rhop(:,:,:) ) ! 3D potential density (sigma0)
......@@ -233,16 +235,16 @@ CONTAINS
! --- momentum --- !
IF ( iom_use("taubot") ) THEN ! bottom stress
zztmp = rho0 * 0.25_wp
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
DO_2D( 0, 0, 0, 0 )
zztmp2 = ( ( rCdU_bot(ji+1,jj)+rCdU_bot(ji ,jj) ) * uu(ji ,jj,mbku(ji ,jj),Kmm) )**2 &
& + ( ( rCdU_bot(ji ,jj)+rCdU_bot(ji-1,jj) ) * uu(ji-1,jj,mbku(ji-1,jj),Kmm) )**2 &
& + ( ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj ) ) * vv(ji,jj ,mbkv(ji,jj ),Kmm) )**2 &
& + ( ( rCdU_bot(ji,jj )+rCdU_bot(ji,jj-1) ) * vv(ji,jj-1,mbkv(ji,jj-1),Kmm) )**2
z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1)
z2d0(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1)
!
END_2D
CALL iom_put( "taubot", z2d )
CALL iom_put( "taubot", z2d0 )
ENDIF
CALL iom_put( "uoce", uu(:,:,:,Kmm) ) ! 3D i-current
......@@ -250,9 +252,9 @@ CONTAINS
IF ( iom_use("sbu") ) THEN
DO_2D( 0, 0, 0, 0 )
ikbot = mbku(ji,jj)
z2d(ji,jj) = uu(ji,jj,ikbot,Kmm)
z2d0(ji,jj) = uu(ji,jj,ikbot,Kmm)
END_2D
CALL iom_put( "sbu", z2d ) ! bottom i-current
CALL iom_put( "sbu", z2d0 ) ! bottom i-current
ENDIF
CALL iom_put( "voce", vv(:,:,:,Kmm) ) ! 3D j-current
......@@ -260,18 +262,15 @@ CONTAINS
IF ( iom_use("sbv") ) THEN
DO_2D( 0, 0, 0, 0 )
ikbot = mbkv(ji,jj)
z2d(ji,jj) = vv(ji,jj,ikbot,Kmm)
z2d0(ji,jj) = vv(ji,jj,ikbot,Kmm)
END_2D
CALL iom_put( "sbv", z2d ) ! bottom j-current
CALL iom_put( "sbv", z2d0 ) ! bottom j-current
ENDIF
! ! vertical velocity
IF( ln_zad_Aimp ) THEN
IF( iom_use('woce') ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = ww(ji,jj,jk) + wi(ji,jj,jk)
END_3D
CALL iom_put( "woce", z3d ) ! explicit plus implicit parts
CALL iom_put( "woce", ww+wi ) ! explicit plus implicit parts
ENDIF
ELSE
CALL iom_put( "woce", ww )
......@@ -281,15 +280,15 @@ CONTAINS
! ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
IF( ln_zad_Aimp ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = rho0 * e1e2t(ji,jj) * ( ww(ji,jj,jk) + wi(ji,jj,jk) )
z3d0(ji,jj,jk) = rho0 * e1e2t(ji,jj) * ( ww(ji,jj,jk) + wi(ji,jj,jk) )
END_3D
ELSE
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = rho0 * e1e2t(ji,jj) * ww(ji,jj,jk)
z3d0(ji,jj,jk) = rho0 * e1e2t(ji,jj) * ww(ji,jj,jk)
END_3D
ENDIF
CALL iom_put( "w_masstr" , z3d )
IF( iom_use('w_masstr2') ) CALL iom_put( "w_masstr2", z3d * z3d )
CALL iom_put( "w_masstr" , z3d0 )
IF( iom_use('w_masstr2') ) CALL iom_put( "w_masstr2", z3d0 * z3d0 )
ENDIF
CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef.
......@@ -304,15 +303,15 @@ CONTAINS
zztmp = ts(ji,jj,1,jp_sal,Kmm)
zztmpx = (ts(ji+1,jj,1,jp_sal,Kmm) - zztmp) * r1_e1u(ji,jj) + (zztmp - ts(ji-1,jj ,1,jp_sal,Kmm)) * r1_e1u(ji-1,jj)
zztmpy = (ts(ji,jj+1,1,jp_sal,Kmm) - zztmp) * r1_e2v(ji,jj) + (zztmp - ts(ji ,jj-1,1,jp_sal,Kmm)) * r1_e2v(ji,jj-1)
z2d(ji,jj) = 0.25_wp * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
z2d0(ji,jj) = 0.25_wp * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
& * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * vmask(ji,jj-1,1)
END_2D
CALL iom_put( "sssgrad2", z2d ) ! square of module of sss gradient
CALL iom_put( "sssgrad2", z2d0 ) ! square of module of sss gradient
IF ( iom_use("sssgrad") ) THEN
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = SQRT( z2d(ji,jj) )
z2d0(ji,jj) = SQRT( z2d0(ji,jj) )
END_2D
CALL iom_put( "sssgrad", z2d ) ! module of sss gradient
CALL iom_put( "sssgrad", z2d0 ) ! module of sss gradient
ENDIF
ENDIF
......@@ -321,80 +320,80 @@ CONTAINS
zztmp = ts(ji,jj,1,jp_tem,Kmm)
zztmpx = ( ts(ji+1,jj,1,jp_tem,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj ,1,jp_tem,Kmm) ) * r1_e1u(ji-1,jj)
zztmpy = ( ts(ji,jj+1,1,jp_tem,Kmm) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - ts(ji ,jj-1,1,jp_tem,Kmm) ) * r1_e2v(ji,jj-1)
z2d(ji,jj) = 0.25_wp * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
z2d0(ji,jj) = 0.25_wp * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
& * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * vmask(ji,jj-1,1)
END_2D
CALL iom_put( "sstgrad2", z2d ) ! square of module of sst gradient
CALL iom_put( "sstgrad2", z2d0 ) ! square of module of sst gradient
IF ( iom_use("sstgrad") ) THEN
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = SQRT( z2d(ji,jj) )
z2d0(ji,jj) = SQRT( z2d0(ji,jj) )
END_2D
CALL iom_put( "sstgrad", z2d ) ! module of sst gradient
CALL iom_put( "sstgrad", z2d0 ) ! module of sst gradient
ENDIF
ENDIF
! heat and salt contents
IF( iom_use("heatc") ) THEN
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk)
z2d0(ji,jj) = z2d0(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk)
END_3D
CALL iom_put( "heatc", rho0_rcp * z2d ) ! vertically integrated heat content (J/m2)
CALL iom_put( "heatc", rho0_rcp * z2d0 ) ! vertically integrated heat content (J/m2)
ENDIF
IF( iom_use("saltc") ) THEN
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
z2d0(ji,jj) = z2d0(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
END_3D
CALL iom_put( "saltc", rho0 * z2d ) ! vertically integrated salt content (PSU*kg/m2)
CALL iom_put( "saltc", rho0 * z2d0 ) ! vertically integrated salt content (PSU*kg/m2)
ENDIF
!
IF( iom_use("salt2c") ) THEN
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
z2d0(ji,jj) = z2d0(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
END_3D
CALL iom_put( "salt2c", rho0 * z2d ) ! vertically integrated square of salt content (PSU2*kg/m2)
CALL iom_put( "salt2c", rho0 * z2d0 ) ! vertically integrated square of salt content (PSU2*kg/m2)
ENDIF
!
IF ( iom_use("ke") .OR. iom_use("ke_int") ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
zztmpx = uu(ji-1,jj ,jk,Kmm) + uu(ji,jj,jk,Kmm)
zztmpy = vv(ji ,jj-1,jk,Kmm) + vv(ji,jj,jk,Kmm)
z3d(ji,jj,jk) = 0.25_wp * ( zztmpx*zztmpx + zztmpy*zztmpy )
z3d0(ji,jj,jk) = 0.25_wp * ( zztmpx*zztmpx + zztmpy*zztmpy )
END_3D
CALL iom_put( "ke", z3d ) ! kinetic energy
CALL iom_put( "ke", z3d0 ) ! kinetic energy
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * e1e2t(ji,jj) * tmask(ji,jj,jk)
z2d0(ji,jj) = z2d0(ji,jj) + e3t(ji,jj,jk,Kmm) * z3d0(ji,jj,jk) * e1e2t(ji,jj) * tmask(ji,jj,jk)
END_3D
CALL iom_put( "ke_int", z2d ) ! vertically integrated kinetic energy
CALL iom_put( "ke_int", z2d0 ) ! vertically integrated kinetic energy
ENDIF
!
IF ( iom_use("sKE") ) THEN ! surface kinetic energy at T point
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = 0.25_wp * ( uu(ji ,jj,1,Kmm) * uu(ji ,jj,1,Kmm) * e1e2u(ji ,jj) * e3u(ji ,jj,1,Kmm) &
z2d0(ji,jj) = 0.25_wp * ( uu(ji ,jj,1,Kmm) * uu(ji ,jj,1,Kmm) * e1e2u(ji ,jj) * e3u(ji ,jj,1,Kmm) &
& + uu(ji-1,jj,1,Kmm) * uu(ji-1,jj,1,Kmm) * e1e2u(ji-1,jj) * e3u(ji-1,jj,1,Kmm) &
& + vv(ji,jj ,1,Kmm) * vv(ji,jj ,1,Kmm) * e1e2v(ji,jj ) * e3v(ji,jj ,1,Kmm) &
& + vv(ji,jj-1,1,Kmm) * vv(ji,jj-1,1,Kmm) * e1e2v(ji,jj-1) * e3v(ji,jj-1,1,Kmm) ) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,1,Kmm) * ssmask(ji,jj)
END_2D
IF ( iom_use("sKE" ) ) CALL iom_put( "sKE" , z2d )
IF ( iom_use("sKE" ) ) CALL iom_put( "sKE" , z2d0 )
ENDIF
!
IF ( iom_use("ssKEf") ) THEN ! surface kinetic energy at F point
z2d(:,:) = 0._wp ! CAUTION : only valid in SWE, not with bathymetry
z2d0(:,:) = 0._wp ! CAUTION : only valid in SWE, not with bathymetry
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = 0.25_wp * ( uu(ji,jj ,1,Kmm) * uu(ji,jj ,1,Kmm) * e1e2u(ji,jj ) * e3u(ji,jj ,1,Kmm) &
z2d0(ji,jj) = 0.25_wp * ( uu(ji,jj ,1,Kmm) * uu(ji,jj ,1,Kmm) * e1e2u(ji,jj ) * e3u(ji,jj ,1,Kmm) &
& + uu(ji,jj+1,1,Kmm) * uu(ji,jj+1,1,Kmm) * e1e2u(ji,jj+1) * e3u(ji,jj+1,1,Kmm) &
& + vv(ji ,jj,1,Kmm) * vv(ji,jj ,1,Kmm) * e1e2v(ji ,jj) * e3v(ji ,jj,1,Kmm) &
& + vv(ji+1,jj,1,Kmm) * vv(ji+1,jj,1,Kmm) * e1e2v(ji+1,jj) * e3v(ji+1,jj,1,Kmm) ) &
& * r1_e1e2f(ji,jj) / e3f(ji,jj,1) * ssfmask(ji,jj)
END_2D
CALL iom_put( "ssKEf", z2d )
CALL iom_put( "ssKEf", z2d0 )
ENDIF
!
CALL iom_put( "hdiv", hdiv ) ! Horizontal divergence
......@@ -402,31 +401,31 @@ CONTAINS
IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = rho0 * uu(ji,jj,jk,Kmm) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * umask(ji,jj,jk)
z3d0(ji,jj,jk) = rho0 * uu(ji,jj,jk,Kmm) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * umask(ji,jj,jk)
END_3D
CALL iom_put( "u_masstr" , z3d ) ! mass transport in i-direction
CALL iom_put( "u_masstr" , z3d0 ) ! mass transport in i-direction
IF( iom_use("u_masstr_vint") ) THEN
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk)
z2d0(ji,jj) = z2d0(ji,jj) + z3d0(ji,jj,jk)
END_3D
CALL iom_put( "u_masstr_vint", z2d ) ! mass transport in i-direction vertical sum
CALL iom_put( "u_masstr_vint", z2d0 ) ! mass transport in i-direction vertical sum
ENDIF
IF( iom_use("u_heattr") ) THEN
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
zztmp = 0.5_wp * rcp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z2d(ji,jj) = z2d(ji,jj) + zztmp * z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji+1,jj,jk,jp_tem,Kmm) )
z2d0(ji,jj) = z2d0(ji,jj) + zztmp * z3d0(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji+1,jj,jk,jp_tem,Kmm) )
END_3D
CALL iom_put( "u_heattr", z2d ) ! heat transport in i-direction
CALL iom_put( "u_heattr", z2d0 ) ! heat transport in i-direction
ENDIF
IF( iom_use("u_salttr") ) THEN
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z2d(ji,jj) = z2d(ji,jj) + 0.5 * z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji+1,jj,jk,jp_sal,Kmm) )
z2d0(ji,jj) = z2d0(ji,jj) + 0.5 * z3d0(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji+1,jj,jk,jp_sal,Kmm) )
END_3D
CALL iom_put( "u_salttr", z2d ) ! heat transport in i-direction
CALL iom_put( "u_salttr", z2d0 ) ! heat transport in i-direction
ENDIF
ENDIF
......@@ -434,41 +433,41 @@ CONTAINS
IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = rho0 * vv(ji,jj,jk,Kmm) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
z3d0(ji,jj,jk) = rho0 * vv(ji,jj,jk,Kmm) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
END_3D
CALL iom_put( "v_masstr", z3d ) ! mass transport in j-direction
CALL iom_put( "v_masstr", z3d0 ) ! mass transport in j-direction
IF( iom_use("v_heattr") ) THEN
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
zztmp = 0.5_wp * rcp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z2d(ji,jj) = z2d(ji,jj) + zztmp * z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji,jj+1,jk,jp_tem,Kmm) )
z2d0(ji,jj) = z2d0(ji,jj) + zztmp * z3d0(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji,jj+1,jk,jp_tem,Kmm) )
END_3D
CALL iom_put( "v_heattr", z2d ) ! heat transport in j-direction
CALL iom_put( "v_heattr", z2d0 ) ! heat transport in j-direction
ENDIF
IF( iom_use("v_salttr") ) THEN
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z2d(ji,jj) = z2d(ji,jj) + 0.5 * z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji,jj+1,jk,jp_sal,Kmm) )
z2d0(ji,jj) = z2d0(ji,jj) + 0.5 * z3d0(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji,jj+1,jk,jp_sal,Kmm) )
END_3D
CALL iom_put( "v_salttr", z2d ) ! heat transport in j-direction
CALL iom_put( "v_salttr", z2d0 ) ! heat transport in j-direction
ENDIF
ENDIF
IF( iom_use("tosmint") ) THEN
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z2d(ji,jj) = z2d(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm)
z2d0(ji,jj) = z2d0(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm)
END_3D
CALL iom_put( "tosmint", z2d ) ! Vertical integral of temperature
CALL iom_put( "tosmint", z2d0 ) ! Vertical integral of temperature
ENDIF
IF( iom_use("somint") ) THEN
z2d(:,:) = 0._wp
z2d0(:,:) = 0._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
z2d(ji,jj) = z2d(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm)
z2d0(ji,jj) = z2d0(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm)
END_3D
CALL iom_put( "somint", z2d ) ! Vertical integral of salinity
CALL iom_put( "somint", z2d0 ) ! Vertical integral of salinity
ENDIF
CALL iom_put( "bn2", rn2 ) ! Brunt-Vaisala buoyancy frequency (N^2)
......@@ -477,45 +476,47 @@ CONTAINS
! Output of surface vorticity terms
!
CALL iom_put( "ssplavor", ff_f ) ! planetary vorticity ( f )
!
IF ( iom_use("ssrelvor") .OR. iom_use("ssEns") .OR. &
IF ( iom_use("ssplavor") .OR. iom_use("ssrelvor") .OR. iom_use("ssEns") .OR. &
& iom_use("ssrelpotvor") .OR. iom_use("ssabspotvor") ) THEN
!
z2d(:,:) = 0._wp
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = ( e2v(ji+1,jj ) * vv(ji+1,jj ,1,Kmm) - e2v(ji,jj) * vv(ji,jj,1,Kmm) &
z2d0(ji,jj) = ff_f(ji,jj)
END_2D
CALL iom_put( "ssplavor", z2d0 ) ! planetary vorticity ( f )
DO_2D( 0, 0, 0, 0 )
z2d0(ji,jj) = ( e2v(ji+1,jj ) * vv(ji+1,jj ,1,Kmm) - e2v(ji,jj) * vv(ji,jj,1,Kmm) &
& - e1u(ji ,jj+1) * uu(ji ,jj+1,1,Kmm) + e1u(ji,jj) * uu(ji,jj,1,Kmm) ) * r1_e1e2f(ji,jj)
END_2D
CALL iom_put( "ssrelvor", z2d ) ! relative vorticity ( zeta )
CALL iom_put( "ssrelvor", z2d0 ) ! relative vorticity ( zeta )
!
IF ( iom_use("ssEns") .OR. iom_use("ssrelpotvor") .OR. iom_use("ssabspotvor") ) THEN
DO_2D( 0, 0, 0, 0 )
ze3 = ( e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1) &
& + e3t(ji,jj ,1,Kmm) * e1e2t(ji,jj ) + e3t(ji+1,jj ,1,Kmm) * e1e2t(ji+1,jj ) ) * r1_e1e2f(ji,jj)
& + e3t(ji,jj ,1,Kmm) * e1e2t(ji,jj ) + e3t(ji+1,jj ,1,Kmm) * e1e2t(ji+1,jj ) ) * r1_e1e2f(ji,jj)
IF( ze3 /= 0._wp ) THEN ; ze3 = 4._wp / ze3
ELSE ; ze3 = 0._wp
ENDIF
z2d(ji,jj) = ze3 * z2d(ji,jj)
z2d0(ji,jj) = ze3 * z2d0(ji,jj)
END_2D
CALL iom_put( "ssrelpotvor", z2d ) ! relative potential vorticity (zeta/h)
CALL iom_put( "ssrelpotvor", z2d0 ) ! relative potential vorticity (zeta/h)
!
IF ( iom_use("ssEns") .OR. iom_use("ssabspotvor") ) THEN
DO_2D( 0, 0, 0, 0 )
ze3 = ( e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1) &
& + e3t(ji,jj ,1,Kmm) * e1e2t(ji,jj ) + e3t(ji+1,jj ,1,Kmm) * e1e2t(ji+1,jj ) ) * r1_e1e2f(ji,jj)
& + e3t(ji,jj ,1,Kmm) * e1e2t(ji,jj ) + e3t(ji+1,jj ,1,Kmm) * e1e2t(ji+1,jj ) ) * r1_e1e2f(ji,jj)
IF( ze3 /= 0._wp ) THEN ; ze3 = 4._wp / ze3
ELSE ; ze3 = 0._wp
ENDIF
z2d(ji,jj) = ze3 * ff_f(ji,jj) + z2d(ji,jj)
z2d0(ji,jj) = ze3 * ff_f(ji,jj) + z2d0(ji,jj)
END_2D
CALL iom_put( "ssabspotvor", z2d ) ! absolute potential vorticity ( q )
CALL iom_put( "ssabspotvor", z2d0 ) ! absolute potential vorticity ( q )
!
IF ( iom_use("ssEns") ) THEN
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = 0.5_wp * z2d(ji,jj) * z2d(ji,jj)
z2d0(ji,jj) = 0.5_wp * z2d0(ji,jj) * z2d0(ji,jj)
END_2D
CALL iom_put( "ssEns", z2d ) ! potential enstrophy ( 1/2*q2 )
CALL iom_put( "ssEns", z2d0 ) ! potential enstrophy ( 1/2*q2 )
ENDIF
ENDIF
ENDIF
......@@ -582,8 +583,8 @@ CONTAINS
INTEGER :: jn, ierror ! local integers
REAL(wp) :: zsto, zout, zmax, zjulian ! local scalars
!
REAL(wp), DIMENSION(jpi,jpj ) :: z2d ! 2D workspace
REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d ! 3D workspace
REAL(wp), DIMENSION(jpi,jpj ) :: z2d0 ! 2D workspace
REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d0 ! 3D workspace
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl ! ABL 3D workspace
!!----------------------------------------------------------------------
!
......@@ -868,7 +869,11 @@ CONTAINS
CALL histdef( nid_T, "sohtc300", "Heat content 300 m" , "J/m2" , & ! htc3
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
#endif
CALL histdef( nid_T, "sozotaux", "Wind Stress along i-axis" , "N/m2" , & ! utau
& jpi, jpj, nh_T, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
CALL histdef( nid_T, "sometauy", "Wind Stress along j-axis" , "N/m2" , & ! vtau
& jpi, jpj, nh_T, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
!
CALL histend( nid_T, snc4chunks=snc4set )
! !!! nid_U : 3D
......@@ -878,10 +883,7 @@ CONTAINS
CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current" , "m/s" , & ! usd
& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
ENDIF
! !!! nid_U : 2D
CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis" , "N/m2" , & ! utau
& jpi, jpj, nh_U, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
!
CALL histend( nid_U, snc4chunks=snc4set )
! !!! nid_V : 3D
......@@ -891,10 +893,7 @@ CONTAINS
CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current" , "m/s" , & ! vsd
& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
ENDIF
! !!! nid_V : 2D
CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis" , "N/m2" , & ! vtau
& jpi, jpj, nh_V, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
!
CALL histend( nid_V, snc4chunks=snc4set )
! !!! nid_W : 3D
......@@ -937,21 +936,21 @@ CONTAINS
IF( .NOT.ln_linssh ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = ts(ji,jj,jk,jp_tem,Kmm) * e3t(ji,jj,jk,Kmm)
z3d0(ji,jj,jk) = ts(ji,jj,jk,jp_tem,Kmm) * e3t(ji,jj,jk,Kmm)
END_3D
CALL histwrite( nid_T, "votemper", it, z3d, ndim_T , ndex_T ) ! heat content
CALL histwrite( nid_T, "votemper", it, z3d0, ndim_T , ndex_T ) ! heat content
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = ts(ji,jj,jk,jp_sal,Kmm) * e3t(ji,jj,jk,Kmm)
z3d0(ji,jj,jk) = ts(ji,jj,jk,jp_sal,Kmm) * e3t(ji,jj,jk,Kmm)
END_3D
CALL histwrite( nid_T, "vosaline", it, z3d, ndim_T , ndex_T ) ! salt content
CALL histwrite( nid_T, "vosaline", it, z3d0, ndim_T , ndex_T ) ! salt content
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj ) = ts(ji,jj, 1,jp_tem,Kmm) * e3t(ji,jj, 1,Kmm)
z2d0(ji,jj ) = ts(ji,jj, 1,jp_tem,Kmm) * e3t(ji,jj, 1,Kmm)
END_2D
CALL histwrite( nid_T, "sosstsst", it, z2d, ndim_hT, ndex_hT ) ! sea surface heat content
CALL histwrite( nid_T, "sosstsst", it, z2d0, ndim_hT, ndex_hT ) ! sea surface heat content
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj ) = ts(ji,jj, 1,jp_sal,Kmm) * e3t(ji,jj, 1,Kmm)
z2d0(ji,jj ) = ts(ji,jj, 1,jp_sal,Kmm) * e3t(ji,jj, 1,Kmm)
END_2D
CALL histwrite( nid_T, "sosaline", it, z2d, ndim_hT, ndex_hT ) ! sea surface salinity content
CALL histwrite( nid_T, "sosaline", it, z2d0, ndim_hT, ndex_hT ) ! sea surface salinity content
ELSE
CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) , ndim_T , ndex_T ) ! temperature
CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) , ndim_T , ndex_T ) ! salinity
......@@ -960,41 +959,41 @@ CONTAINS
ENDIF
IF( .NOT.ln_linssh ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = e3t(ji,jj,jk,Kmm) ! 3D workspace for qco substitution
z3d0(ji,jj,jk) = e3t(ji,jj,jk,Kmm) ! 3D workspace for qco substitution
END_3D
CALL histwrite( nid_T, "vovvle3t", it, z3d , ndim_T , ndex_T ) ! level thickness
CALL histwrite( nid_T, "vovvle3t", it, z3d0 , ndim_T , ndex_T ) ! level thickness
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = gdept(ji,jj,jk,Kmm) ! 3D workspace for qco substitution
z3d0(ji,jj,jk) = gdept(ji,jj,jk,Kmm) ! 3D workspace for qco substitution
END_3D
CALL histwrite( nid_T, "vovvldep", it, z3d , ndim_T , ndex_T ) ! t-point depth
CALL histwrite( nid_T, "vovvldep", it, z3d0 , ndim_T , ndex_T ) ! t-point depth
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = ( ( e3t(ji,jj,jk,Kmm) - e3t_0(ji,jj,jk) ) / e3t_0(ji,jj,jk) * 100._wp * tmask(ji,jj,jk) ) ** 2
z3d0(ji,jj,jk) = ( ( e3t(ji,jj,jk,Kmm) - e3t_0(ji,jj,jk) ) / e3t_0(ji,jj,jk) * 100._wp * tmask(ji,jj,jk) ) ** 2
END_3D
CALL histwrite( nid_T, "vovvldef", it, z3d , ndim_T , ndex_T ) ! level thickness deformation
CALL histwrite( nid_T, "vovvldef", it, z3d0 , ndim_T , ndex_T ) ! level thickness deformation
ENDIF
CALL histwrite( nid_T, "sossheig", it, ssh(:,:,Kmm) , ndim_hT, ndex_hT ) ! sea surface height
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = emp(ji,jj) - rnf(ji,jj)
z2d0(ji,jj) = emp(ji,jj) - rnf(ji,jj)
END_2D
CALL histwrite( nid_T, "sowaflup", it, z2d , ndim_hT, ndex_hT ) ! upward water flux
CALL histwrite( nid_T, "sowaflup", it, z2d0 , ndim_hT, ndex_hT ) ! upward water flux
CALL histwrite( nid_T, "sorunoff", it, rnf , ndim_hT, ndex_hT ) ! river runoffs
CALL histwrite( nid_T, "sosfldow", it, sfx , ndim_hT, ndex_hT ) ! downward salt flux
! (includes virtual salt flux beneath ice
! in linear free surface case)
IF( ln_linssh ) THEN
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = emp (ji,jj) * ts(ji,jj,1,jp_tem,Kmm)
z2d0(ji,jj) = emp (ji,jj) * ts(ji,jj,1,jp_tem,Kmm)
END_2D
CALL histwrite( nid_T, "sosst_cd", it, z2d, ndim_hT, ndex_hT ) ! c/d term on sst
CALL histwrite( nid_T, "sosst_cd", it, z2d0, ndim_hT, ndex_hT ) ! c/d term on sst
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = emp (ji,jj) * ts(ji,jj,1,jp_sal,Kmm)
z2d0(ji,jj) = emp (ji,jj) * ts(ji,jj,1,jp_sal,Kmm)
END_2D
CALL histwrite( nid_T, "sosss_cd", it, z2d, ndim_hT, ndex_hT ) ! c/d term on sss
CALL histwrite( nid_T, "sosss_cd", it, z2d0, ndim_hT, ndex_hT ) ! c/d term on sss
ENDIF
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = qsr(ji,jj) + qns(ji,jj)
z2d0(ji,jj) = qsr(ji,jj) + qns(ji,jj)
END_2D
CALL histwrite( nid_T, "sohefldo", it, z2d , ndim_hT, ndex_hT ) ! total heat flux
CALL histwrite( nid_T, "sohefldo", it, z2d0 , ndim_hT, ndex_hT ) ! total heat flux
CALL histwrite( nid_T, "soshfldo", it, qsr , ndim_hT, ndex_hT ) ! solar heat flux
IF( ALLOCATED(hmld) ) THEN ! zdf_mxl not called by SWE
CALL histwrite( nid_T, "somixhgt", it, hmld , ndim_hT, ndex_hT ) ! turbocline depth
......@@ -1053,9 +1052,9 @@ CONTAINS
CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping
CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = erp(ji,jj) * ts(ji,jj,1,jp_sal,Kmm) * tmask(ji,jj,1)
z2d0(ji,jj) = erp(ji,jj) * ts(ji,jj,1,jp_sal,Kmm) * tmask(ji,jj,1)
END_2D
CALL histwrite( nid_T, "sosafldp", it, z2d , ndim_hT, ndex_hT ) ! salt flux damping
CALL histwrite( nid_T, "sosafldp", it, z2d0 , ndim_hT, ndex_hT ) ! salt flux damping
ENDIF
! zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
! CALL histwrite( nid_T, "sobowlin", it, zw2d , ndim_hT, ndex_hT ) ! ???
......@@ -1066,18 +1065,18 @@ CONTAINS
CALL histwrite( nid_T, "so28chgt", it, hd28 , ndim_hT, ndex_hT ) ! depth of the 28 isotherm
CALL histwrite( nid_T, "sohtc300", it, htc3 , ndim_hT, ndex_hT ) ! first 300m heaat content
#endif
CALL histwrite( nid_T, "sozotaux", it, utau , ndim_hT, ndex_hT ) ! i-wind stress
CALL histwrite( nid_T, "sometauy", it, vtau , ndim_hT, ndex_hT ) ! j-wind stress
CALL histwrite( nid_U, "vozocrtx", it, uu(:,:,:,Kmm) , ndim_U , ndex_U ) ! i-current
CALL histwrite( nid_U, "sozotaux", it, utau , ndim_hU, ndex_hU ) ! i-wind stress
CALL histwrite( nid_V, "vomecrty", it, vv(:,:,:,Kmm) , ndim_V , ndex_V ) ! j-current
CALL histwrite( nid_V, "sometauy", it, vtau , ndim_hV, ndex_hV ) ! j-wind stress
IF( ln_zad_Aimp ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = ww(ji,jj,jk) + wi(ji,jj,jk)
z3d0(ji,jj,jk) = ww(ji,jj,jk) + wi(ji,jj,jk)
END_3D
CALL histwrite( nid_W, "vovecrtz", it, z3d , ndim_T, ndex_T ) ! vert. current
CALL histwrite( nid_W, "vovecrtz", it, z3d0 , ndim_T, ndex_T ) ! vert. current
ELSE
CALL histwrite( nid_W, "vovecrtz", it, ww , ndim_T, ndex_T ) ! vert. current
ENDIF
......@@ -1126,8 +1125,8 @@ CONTAINS
!!
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: inum
REAL(wp), DIMENSION(jpi,jpj) :: z2d
REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d
REAL(wp), DIMENSION(A2D(0)) :: z2d0
REAL(wp), DIMENSION(A2D(0),jpk) :: z3d0
!!----------------------------------------------------------------------
!
IF(lwp) THEN
......@@ -1146,9 +1145,9 @@ CONTAINS
CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm) ) ! now j-velocity
IF( ln_zad_Aimp ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = ww(ji,jj,jk) + wi(ji,jj,jk)
z3d0(ji,jj,jk) = ww(ji,jj,jk) + wi(ji,jj,jk)
END_3D
CALL iom_rstput( 0, 0, inum, 'vovecrtz', z3d ) ! now k-velocity
CALL iom_rstput( 0, 0, inum, 'vovecrtz', z3d0 ) ! now k-velocity
ELSE
CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww ) ! now k-velocity
ENDIF
......@@ -1184,26 +1183,26 @@ CONTAINS
CALL iom_rstput( 0, 0, inum, 'ahmf', ahmf ) ! ahmf at v-point
ENDIF
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = emp(ji,jj) - rnf(ji,jj)
z2d0(ji,jj) = emp(ji,jj) - rnf(ji,jj)
END_2D
CALL iom_rstput( 0, 0, inum, 'sowaflup', z2d ) ! freshwater budget
CALL iom_rstput( 0, 0, inum, 'sowaflup', z2d0 ) ! freshwater budget
DO_2D( 0, 0, 0, 0 )
z2d(ji,jj) = qsr(ji,jj) + qns(ji,jj)
z2d0(ji,jj) = qsr(ji,jj) + qns(ji,jj)
END_2D
CALL iom_rstput( 0, 0, inum, 'sohefldo', z2d ) ! total heat flux
CALL iom_rstput( 0, 0, inum, 'sohefldo', z2d0 ) ! total heat flux
CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr ) ! solar heat flux
CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i ) ! ice fraction
CALL iom_rstput( 0, 0, inum, 'sozotaux', utau ) ! i-wind stress
CALL iom_rstput( 0, 0, inum, 'sometauy', vtau ) ! j-wind stress
IF( .NOT.ln_linssh ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = gdept(ji,jj,jk,Kmm) ! 3D workspace for qco substitution
z3d0(ji,jj,jk) = gdept(ji,jj,jk,Kmm) ! 3D workspace for qco substitution
END_3D
CALL iom_rstput( 0, 0, inum, 'vovvldep', z3d ) ! T-cell depth
CALL iom_rstput( 0, 0, inum, 'vovvldep', z3d0 ) ! T-cell depth
DO_3D( 0, 0, 0, 0, 1, jpk )
z3d(ji,jj,jk) = e3t(ji,jj,jk,Kmm) ! 3D workspace for qco substitution
z3d0(ji,jj,jk) = e3t(ji,jj,jk,Kmm) ! 3D workspace for qco substitution
END_3D
CALL iom_rstput( 0, 0, inum, 'vovvle3t', z3d ) ! T-cell thickness
CALL iom_rstput( 0, 0, inum, 'vovvle3t', z3d0 ) ! T-cell thickness
END IF
IF( ln_wave .AND. ln_sdw ) THEN
CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd ) ! now StokesDrift i-velocity
......
src/OCE/DIU/diu_bulk.F90
View file @ 3efbcb17
......@@ -64,7 +64,7 @@ CONTAINS
IF( ln_diurnal ) THEN
!
ALLOCATE( x_dsst(jpi,jpj), x_solfrac(jpi,jpj) )
ALLOCATE( x_dsst(A2D(0)), x_solfrac(A2D(0)) )
!
x_solfrac = 0._wp ! Initialise the solar fraction
x_dsst = 0._wp
......@@ -92,25 +92,25 @@ CONTAINS
!! ** Reference : Refinements to a prognostic scheme of skin sea surface
!! temperature, Takaya et al, JGR, 2010
!!----------------------------------------------------------------------
INTEGER , INTENT(in) :: kt ! time step
REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: psolflux ! solar flux (Watts)
REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: pqflux ! heat (non-solar) flux (Watts)
REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: ptauflux ! wind stress (kg/ m s^2)
REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: prho ! water density (kg/m^3)
REAL(wp) , INTENT(in) :: p_rdt ! time-step
REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: pLa ! Langmuir number
REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: pthick ! warm layer thickness (m)
REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: pcoolthick ! cool skin thickness (m)
REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: pmu ! mu parameter
REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: p_hflux_bkginc ! increment to the heat flux
REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) :: p_fvel_bkginc ! increment to the friction velocity
INTEGER , INTENT(in) :: kt ! time step
REAL(wp), DIMENSION(A2D(0)) , INTENT(in) :: psolflux ! solar flux (Watts)
REAL(wp), DIMENSION(A2D(0)) , INTENT(in) :: pqflux ! heat (non-solar) flux (Watts)
REAL(wp), DIMENSION(A2D(0)) , INTENT(in) :: ptauflux ! wind stress (kg/ m s^2)
REAL(wp), DIMENSION(jpi,jpj) , INTENT(in) :: prho ! water density (kg/m^3)
REAL(wp) , INTENT(in) :: p_rdt ! time-step
REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in) :: pLa ! Langmuir number
REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in) :: pthick ! warm layer thickness (m)
REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in) :: pcoolthick ! cool skin thickness (m)
REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in) :: pmu ! mu parameter
REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in) :: p_hflux_bkginc ! increment to the heat flux
REAL(wp), DIMENSION(A2D(0)), OPTIONAL, INTENT(in) :: p_fvel_bkginc ! increment to the friction velocity
!
INTEGER :: ji,jj
LOGICAL :: ll_calcfrac
REAL(wp), DIMENSION(jpi,jpj) :: z_fvel ! friction velocity
REAL(wp), DIMENSION(jpi,jpj) :: zthick, zcoolthick, zmu, zla
REAL(wp), DIMENSION(jpi,jpj) :: z_abflux ! absorbed flux
REAL(wp), DIMENSION(jpi,jpj) :: z_fla ! Langmuir function value
REAL(wp), DIMENSION(A2D(0)) :: z_fvel ! friction velocity
REAL(wp), DIMENSION(A2D(0)) :: zthick, zcoolthick, zmu, zla
REAL(wp), DIMENSION(A2D(0)) :: z_abflux ! absorbed flux
REAL(wp), DIMENSION(A2D(0)) :: z_fla ! Langmuir function value
!!----------------------------------------------------------------------
! Set optional arguments to their defaults
......@@ -129,14 +129,14 @@ CONTAINS
! If not done already, calculate the solar fraction
IF ( kt==nit000 ) THEN
DO_2D( 1, 1, 1, 1 )
IF( ( x_solfrac(ji,jj) == 0._wp ) .AND. ( tmask(ji,jj,1) == 1._wp ) ) &
DO_2D( 0, 0, 0, 0 )
IF( ( x_solfrac(ji,jj) == 0._wp ) .AND. ( smask0(ji,jj) == 1._wp ) ) &
& x_solfrac(ji,jj) = solfrac( zcoolthick(ji,jj),zthick(ji,jj) )
END_2D
ENDIF
! convert solar flux and heat flux to absorbed flux
WHERE ( tmask(:,:,1) == 1._wp)
WHERE ( smask0(:,:) == 1._wp)
z_abflux(:,:) = ( x_solfrac(:,:) * psolflux (:,:)) + pqflux(:,:)
ELSEWHERE
z_abflux(:,:) = 0._wp
......@@ -147,8 +147,8 @@ CONTAINS
ENDWHERE
! Calculate the friction velocity
WHERE ( (ptauflux /= 0) .AND. ( tmask(:,:,1) == 1.) )
z_fvel(:,:) = SQRT( ptauflux(:,:) / prho(:,:) )
WHERE ( (ptauflux(:,:) /= 0) .AND. ( smask0(:,:) == 1.) )
z_fvel(:,:) = SQRT( ptauflux(:,:) / prho(A2D(0)) )
ELSEWHERE
z_fvel(:,:) = 0._wp
ENDWHERE
......@@ -157,7 +157,7 @@ CONTAINS
! Calculate the Langmuir function value
WHERE ( tmask(:,:,1) == 1.)
WHERE ( smask0(:,:) == 1.)
z_fla(:,:) = MAX( 1._wp, zla(:,:)**( -2._wp / 3._wp ) )
ELSEWHERE
z_fla(:,:) = 0._wp
......@@ -176,16 +176,16 @@ CONTAINS
IMPLICIT NONE
! Function definition
REAL(wp), DIMENSION(jpi,jpj) :: t_imp
REAL(wp), DIMENSION(A2D(0)) :: t_imp
! Dummy variables
REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_dsst ! Delta SST
REAL(wp), INTENT(IN) :: p_rdt ! Time-step
REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_abflux ! Heat forcing
REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_fvel ! Friction velocity
REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_fla ! Langmuir number
REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pmu ! Structure parameter
REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: pthick ! Layer thickness
REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: prho ! Water density
REAL(wp), DIMENSION(A2D(0)), INTENT(IN) :: p_dsst ! Delta SST
REAL(wp), INTENT(IN) :: p_rdt ! Time-step
REAL(wp), DIMENSION(A2D(0)), INTENT(IN) :: p_abflux ! Heat forcing
REAL(wp), DIMENSION(A2D(0)), INTENT(IN) :: p_fvel ! Friction velocity
REAL(wp), DIMENSION(A2D(0)), INTENT(IN) :: p_fla ! Langmuir number
REAL(wp), DIMENSION(A2D(0)), INTENT(IN) :: pmu ! Structure parameter
REAL(wp), DIMENSION(A2D(0)), INTENT(IN) :: pthick ! Layer thickness
REAL(wp), DIMENSION(jpi,jpj),INTENT(IN) :: prho ! Water density
! Local variables
REAL(wp) :: z_olength ! Obukhov length
......@@ -198,10 +198,10 @@ CONTAINS
INTEGER :: ji,jj
DO_2D( 1, 1, 1, 1 )
DO_2D( 0, 0, 0, 0 )
! Only calculate outside tmask
IF ( tmask(ji,jj,1) /= 1._wp ) THEN
IF ( smask0(ji,jj) /= 1._wp ) THEN
t_imp(ji,jj) = 0._wp
CYCLE
END IF
......
src/OCE/DIU/diu_coolskin.F90
View file @ 3efbcb17
......@@ -59,7 +59,7 @@ MODULE diu_coolskin
!! ** Reference :
!!
!!----------------------------------------------------------------------
ALLOCATE( x_csdsst(jpi,jpj), x_csthick(jpi,jpj) )
ALLOCATE( x_csdsst(A2D(0)), x_csthick(A2D(0)) )
x_csdsst = 0.
x_csthick = 0.
!
......@@ -77,16 +77,16 @@ MODULE diu_coolskin
!! ** Reference :
!!----------------------------------------------------------------------
! Dummy variables
REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: psqflux ! Heat (non-solar)(Watts)
REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: pstauflux ! Wind stress (kg/ m s^2)
REAL(wp), INTENT(IN), DIMENSION(A2D(0)) :: psqflux ! Heat (non-solar)(Watts)
REAL(wp), INTENT(IN), DIMENSION(A2D(0)) :: pstauflux ! Wind stress (kg/ m s^2)
REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: psrho ! Water density (kg/m^3)
REAL(wp), INTENT(IN) :: pDt ! Time-step
! Local variables
REAL(wp), DIMENSION(jpi,jpj) :: z_fv ! Friction velocity
REAL(wp), DIMENSION(jpi,jpj) :: z_gamma ! Dimensionless function of wind speed
REAL(wp), DIMENSION(jpi,jpj) :: z_lamda ! Sauders (dimensionless) proportionality constant
REAL(wp), DIMENSION(jpi,jpj) :: z_wspd ! Wind speed (m/s)
REAL(wp), DIMENSION(A2D(0)) :: z_fv ! Friction velocity
REAL(wp), DIMENSION(A2D(0)) :: z_gamma ! Dimensionless function of wind speed
REAL(wp), DIMENSION(A2D(0)) :: z_lamda ! Sauders (dimensionless) proportionality constant
REAL(wp), DIMENSION(A2D(0)) :: z_wspd ! Wind speed (m/s)
REAL(wp) :: z_ztx ! Temporary u wind stress
REAL(wp) :: z_zty ! Temporary v wind stress
REAL(wp) :: z_zmod ! Temporary total wind stress
......@@ -96,10 +96,10 @@ MODULE diu_coolskin
!
IF( .NOT. (ln_blk .OR. ln_abl) ) CALL ctl_stop("diu_coolskin.f90: diurnal flux processing only implemented for bulk forcing")
!
DO_2D( 1, 1, 1, 1 )
DO_2D( 0, 0, 0, 0 )
!
! Calcualte wind speed from wind stress and friction velocity
IF( tmask(ji,jj,1) == 1. .AND. pstauflux(ji,jj) /= 0 .AND. psrho(ji,jj) /=0 ) THEN
IF( smask0(ji,jj) == 1. .AND. pstauflux(ji,jj) /= 0 .AND. psrho(ji,jj) /=0 ) THEN
z_fv(ji,jj) = SQRT( pstauflux(ji,jj) / psrho(ji,jj) )
z_wspd(ji,jj) = SQRT( pstauflux(ji,jj) / ( pp_cda * pp_rhoa ) )
ELSE
......@@ -108,28 +108,28 @@ MODULE diu_coolskin
ENDIF
!
! Calculate gamma function which is dependent upon wind speed
IF( tmask(ji,jj,1) == 1. ) THEN
IF( smask0(ji,jj) == 1. ) THEN
IF( ( z_wspd(ji,jj) <= 7.5 ) ) z_gamma(ji,jj) = ( 0.2 * z_wspd(ji,jj) ) + 0.5
IF( ( z_wspd(ji,jj) > 7.5 ) .AND. ( z_wspd(ji,jj) < 10. ) ) z_gamma(ji,jj) = ( 1.6 * z_wspd(ji,jj) ) - 10.
IF( ( z_wspd(ji,jj) >= 10. ) ) z_gamma(ji,jj) = 6.
ENDIF
!
! Calculate lamda function
IF( tmask(ji,jj,1) == 1. .AND. z_fv(ji,jj) /= 0 ) THEN
IF( smask0(ji,jj) == 1. .AND. z_fv(ji,jj) /= 0 ) THEN
z_lamda(ji,jj) = ( z_fv(ji,jj) * pp_k * pp_C ) / ( z_gamma(ji,jj) * psrho(ji,jj) * pp_cw * pp_h * pp_v )
ELSE
z_lamda(ji,jj) = 0.
ENDIF
!
! Calculate the cool skin thickness - only when heat flux is out of the ocean
IF( tmask(ji,jj,1) == 1. .AND. z_fv(ji,jj) /= 0 .AND. psqflux(ji,jj) < 0 ) THEN
IF( smask0(ji,jj) == 1. .AND. z_fv(ji,jj) /= 0 .AND. psqflux(ji,jj) < 0 ) THEN
x_csthick(ji,jj) = ( z_lamda(ji,jj) * pp_v ) / z_fv(ji,jj)
ELSE
x_csthick(ji,jj) = 0.
ENDIF
!
! Calculate the cool skin correction - only when the heat flux is out of the ocean
IF( tmask(ji,jj,1) == 1. .AND. x_csthick(ji,jj) /= 0. .AND. psqflux(ji,jj) < 0. ) THEN
IF( smask0(ji,jj) == 1. .AND. x_csthick(ji,jj) /= 0. .AND. psqflux(ji,jj) < 0. ) THEN
x_csdsst(ji,jj) = ( psqflux(ji,jj) * x_csthick(ji,jj) ) / pp_k
ELSE
x_csdsst(ji,jj) = 0.
......
src/OCE/DIU/step_diu.F90
View file @ 3efbcb17
......@@ -46,8 +46,7 @@ MODULE step_diu
!!----------------------------------------------------------------------
INTEGER :: jk ! dummy loop indices
INTEGER :: indic ! error indicator if < 0
REAL(wp), DIMENSION(jpi,jpj) :: z_fvel_bkginc, z_hflux_bkginc
INTEGER :: Nbb, Nnn, Naa, Nrhs ! local definitions as placeholders for now
INTEGER :: Nbb, Nnn, Naa, Nrhs ! local definitions as placeholders for now
!! ---------------------------------------------------------------------
IF(ln_diurnal_only) THEN
......
src/OCE/DOM/dom_oce.F90
View file @ 3efbcb17
......@@ -27,6 +27,8 @@ MODULE dom_oce
PUBLIC dom_oce_alloc ! Called from nemogcm.F90
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! time & space domain namelist
!! ----------------------------
......@@ -74,10 +76,10 @@ MODULE dom_oce
INTEGER :: nn_ltile_i, nn_ltile_j
! Domain tiling
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ntsi_a !: start of internal part of tile domain
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ntsj_a !
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ntei_a !: end of internal part of tile domain
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ntej_a !
INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntsi_a !: start of internal part of tile domain
INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntsj_a !
INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntei_a !: end of internal part of tile domain
INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntej_a !
LOGICAL, PUBLIC :: l_istiled ! whether tiling is currently active or not
! !: domain MPP decomposition parameters
......@@ -85,32 +87,30 @@ MODULE dom_oce
INTEGER , PUBLIC :: narea !: number for local area (starting at 1) = MPI rank + 1
INTEGER, PUBLIC :: nidom !: IOIPSL things...
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mig !: local ==> global domain, including halos (jpiglo), i-index
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mjg !: local ==> global domain, including halos (jpjglo), j-index
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mig0 !: local ==> global domain, excluding halos (Ni0glo), i-index
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mjg0 !: local ==> global domain, excluding halos (Nj0glo), j-index
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mi0, mi1 !: global, including halos (jpiglo) ==> local domain i-index
INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: mig !: local ==> global domain, i-index
INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: mjg !: local ==> global domain, j-index
INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: mi0, mi1 !: global ==> local domain, i-index
! !: (mi0=1 and mi1=0 if global index not in local domain)
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: mj0, mj1 !: global, including halos (jpjglo) ==> local domain j-index
INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: mj0, mj1 !: global ==> local domain, j-index
! !: (mj0=1 and mj1=0 if global index not in local domain)
INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: nfimpp, nfproc, nfjpi
INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: nfimpp, nfproc, nfjpi, nfni_0
!!----------------------------------------------------------------------
!! horizontal curvilinear coordinate and scale factors
!! ---------------------------------------------------------------------
REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: glamt , glamu, glamv , glamf !: longitude at t, u, v, f-points [degree]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: gphit , gphiu, gphiv , gphif !: latitude at t, u, v, f-points [degree]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1t , e2t , r1_e1t, r1_e2t !: t-point horizontal scale factors [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1u , e2u , r1_e1u, r1_e2u !: horizontal scale factors at u-point [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1v , e2v , r1_e1v, r1_e2v !: horizontal scale factors at v-point [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1f , e2f , r1_e1f, r1_e2f !: horizontal scale factors at f-point [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: glamt , glamu, glamv , glamf !: longitude at t, u, v, f-points [degree]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: gphit , gphiu, gphiv , gphif !: latitude at t, u, v, f-points [degree]
REAL(wp), PUBLIC, ALLOCATABLE, TARGET, DIMENSION(:,:) :: e1t , e2t , r1_e1t, r1_e2t !: t-point horizontal scale factors [m]
REAL(wp), PUBLIC, ALLOCATABLE, TARGET, DIMENSION(:,:) :: e1u , e2u , r1_e1u, r1_e2u !: horizontal scale factors at u-point [m]
REAL(wp), PUBLIC, ALLOCATABLE, TARGET, DIMENSION(:,:) :: e1v , e2v , r1_e1v, r1_e2v !: horizontal scale factors at v-point [m]
REAL(wp), PUBLIC, ALLOCATABLE, TARGET, DIMENSION(:,:) :: e1f , e2f , r1_e1f, r1_e2f !: horizontal scale factors at f-point [m]
!
REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: e1e2t , r1_e1e2t !: associated metrics at t-point
REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: e1e2u , r1_e1e2u , e2_e1u !: associated metrics at u-point
REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: e1e2v , r1_e1e2v , e1_e2v !: associated metrics at v-point
REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: e1e2f , r1_e1e2f !: associated metrics at f-point
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: e1e2t , r1_e1e2t !: associated metrics at t-point
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: e1e2u , r1_e1e2u , e2_e1u !: associated metrics at u-point
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: e1e2v , r1_e1e2v , e1_e2v !: associated metrics at v-point
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: e1e2f , r1_e1e2f !: associated metrics at f-point
!
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ff_f , ff_t !: Coriolis factor at f- & t-points [1/s]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: ff_f , ff_t !: Coriolis factor at f- & t-points [1/s]
!!----------------------------------------------------------------------
!! vertical coordinate and scale factors
......@@ -130,75 +130,77 @@ MODULE dom_oce
LOGICAL, PUBLIC :: ln_sco !: s-coordinate or hybrid z-s coordinate
LOGICAL, PUBLIC :: ln_isfcav !: presence of ISF
! ! reference scale factors
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_0 !: t- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3u_0 !: u- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3v_0 !: v- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3f_0 !: f- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3w_0 !: w- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3uw_0 !: uw-vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3vw_0 !: vw-vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: e3t_0 !: t- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: e3u_0 !: u- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: e3v_0 !: v- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: e3f_0 !: f- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: e3w_0 !: w- vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: e3uw_0 !: uw-vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: e3vw_0 !: vw-vert. scale factor [m]
! ! time-dependent scale factors (domvvl)
#if defined key_qco || defined key_linssh
#else
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e3t, e3u, e3v, e3w, e3uw, e3vw !: vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3f !: F-point vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:,:) :: e3t, e3u, e3v, e3w, e3uw, e3vw !: vert. scale factor [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: e3f !: F-point vert. scale factor [m]
#endif
! ! time-dependent ratio ssh / h_0 (domqco)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: r3t, r3u, r3v !: time-dependent ratio at t-, u- and v-point [-]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r3f !: mid-time-level ratio at f-point [-]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r3t_f, r3u_f, r3v_f !: now time-filtered ratio at t-, u- and v-point [-]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: r3t, r3u, r3v !: time-dependent ratio at t-, u- and v-point [-]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: r3f !: mid-time-level ratio at f-point [-]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: r3t_f, r3u_f, r3v_f !: now time-filtered ratio at t-, u- and v-point [-]
! ! reference depths of cells
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdept_0 !: t- depth [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdepw_0 !: w- depth [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: gdept_0 !: t- depth [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: gdepw_0 !: w- depth [m]
! ! time-dependent depths of cells (domvvl)
#if defined key_qco || defined key_linssh
#else
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gde3w_0, gde3w !: w- depth (sum of e3w) [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: gdept, gdepw
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: gde3w_0, gde3w !: w- depth (sum of e3w) [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:,:) :: gdept, gdepw
#endif
! ! reference heights of ocean water column and its inverse
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht_0, r1_ht_0 !: t-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_0, r1_hu_0 !: u-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hv_0, r1_hv_0 !: v-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hf_0, r1_hf_0 !: f-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: ht_0, r1_ht_0 !: t-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: hu_0, r1_hu_0 !: u-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: hv_0, r1_hv_0 !: v-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: hf_0, r1_hf_0 !: f-depth [m] and [1/m]
! ! time-dependent heights of ocean water column (domvvl)
#if defined key_qco || defined key_linssh
#else
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht !: t-points [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hu, r1_hu !: u-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hv, r1_hv !: v-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: ht !: t-points [m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: hu, r1_hu !: u-depth [m] and [1/m]
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: hv, r1_hv !: v-depth [m] and [1/m]
#endif
INTEGER, PUBLIC :: nla10 !: deepest W level Above ~10m (nlb10 - 1)
INTEGER, PUBLIC :: nlb10 !: shallowest W level Bellow ~10m (nla10 + 1)
!! 1D reference vertical coordinate
!! =-----------------====------
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: gdept_1d, gdepw_1d !: reference depth of t- and w-points (m)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: e3t_1d , e3w_1d !: reference vertical scale factors at T- and W-pts (m)
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: gdept_1d, gdepw_1d !: reference depth of t- and w-points (m)
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: e3t_1d , e3w_1d !: reference vertical scale factors at T- and W-pts (m)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: risfdep, bathy
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: risfdep, bathy
!!----------------------------------------------------------------------
!! masks, top and bottom ocean point position
!! ---------------------------------------------------------------------
!!gm Proposition of new name for top/bottom vertical indices
! INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mtk_t, mtk_u, mtk_v !: top first wet T-, U-, and V-level (ISF)
! INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbk_t, mbk_u, mbk_v !: bottom last wet T-, U-, and V-level
! INTEGER , PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: mtk_t, mtk_u, mtk_v !: top first wet T-, U-, and V-level (ISF)
! INTEGER , PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: mbk_t, mbk_u, mbk_v !: bottom last wet T-, U-, and V-level
!!gm
INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt, mbku, mbkv, mbkf !: bottom last wet T-, U-, V- and F-level
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_i !: interior (excluding halos+duplicated points) domain T-point mask
INTEGER , PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: mbkt, mbku, mbkv, mbkf !: bottom last wet T-, U-, V- and F-level
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tmask_i !: interior (excluding halos+duplicated points) domain T-point mask
INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mikt, miku, mikv, mikf !: top first wet T-, U-, V-, F-level (ISF)
INTEGER , PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: mikt, miku, mikv, mikf !: top first wet T-, U-, V-, F-level (ISF)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssmask, ssumask, ssvmask, ssfmask !: surface mask at T-,U-, V- and F-pts
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, wmask, fmask !: land/ocean mask at T-, U-, V-, W- and F-pts
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wumask, wvmask !: land/ocean mask at WU- and WV-pts
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: fe3mask !: land/ocean mask at F-pts
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_upd, umask_upd, vmask_upd !: land/ocean mask at F-pts
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: smask0 !: surface mask at T-pts on inner domain
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: smask0_i !: equivalent of tmask_i for inner domain
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: ssmask, ssumask, ssvmask, ssfmask !: surface mask at T-,U-, V- and F-pts
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, wmask, fmask !: land/ocean mask at T-, U-, V-, W- and F-pts
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:), TARGET :: wumask, wvmask !: land/ocean mask at WU- and WV-pts
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:), TARGET :: fe3mask !: land/ocean mask at F-pts
REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: tmask_upd, umask_upd, vmask_upd !: land/ocean mask at F-pts
!!----------------------------------------------------------------------
!! calendar variables
......@@ -326,7 +328,7 @@ CONTAINS
ALLOCATE( risfdep(jpi,jpj) , bathy(jpi,jpj) , STAT=ierr(ii) )
!
ii = ii+1
ALLOCATE( tmask_i(jpi,jpj) , &
ALLOCATE( tmask_i(jpi,jpj) , smask0(Nis0-(0):Nie0+(0),Njs0-(0):Nje0+(0)) , smask0_i(Nis0-(0):Nie0+(0),Njs0-(0):Nje0+(0)) , &
& ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) , &
& mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv (jpi,jpj) , mbkf(jpi,jpj) , STAT=ierr(ii) )
!
......
src/OCE/DOM/dommsk.F90
View file @ 3efbcb17
......@@ -144,7 +144,8 @@ CONTAINS
tmask(ji,jj,jk) = tmask(ji,jj,jk) * bdytmask(ji,jj)
END_3D
ENDIF
smask0(:,:) = tmask(A2D(0),1)
! Ocean/land mask at u-, v-, and f-points (computed from tmask)
! ----------------------------------------
! NB: at this point, fmask is designed for free slip lateral boundary condition
......@@ -194,7 +195,8 @@ CONTAINS
! --------------------
!
CALL dom_uniq( tmask_i, 'T' )
tmask_i(:,:) = ssmask(:,:) * tmask_i(:,:)
tmask_i (:,:) = ssmask(:,:) * tmask_i(:,:)
smask0_i(:,:) = tmask_i(A2D(0))
! Lateral boundary conditions on velocity (modify fmask)
! ---------------------------------------
......
src/OCE/DOM/domtile.F90
View file @ 3efbcb17
......@@ -52,16 +52,6 @@ CONTAINS
!!----------------------------------------------------------------------
IF( ln_tile .AND. nn_hls /= 2 ) CALL ctl_stop('dom_tile_init: Tiling is only supported for nn_hls = 2')
ntile = 0 ! Initialise to full domain
nijtile = 1
ntsi = Nis0
ntsj = Njs0
ntei = Nie0
ntej = Nje0
nthl = 0
nthr = 0
nthb = 0
ntht = 0
l_istiled = .FALSE.
IF( ln_tile ) THEN ! Calculate tile domain indices
......
src/OCE/DOM/domvvl.F90
View file @ 3efbcb17
......@@ -282,8 +282,8 @@ CONTAINS
IF( nn_cfg == 3 ) THEN ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2
ii0 = 103 + nn_hls - 1 ; ii1 = 111 + nn_hls - 1
ij0 = 128 + nn_hls ; ij1 = 135 + nn_hls
frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.0_wp
frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0_wp / rn_Dt
frq_rst_e3t( mi0(ii0,nn_hls):mi1(ii1,nn_hls) , mj0(ij0,nn_hls):mj1(ij1,nn_hls) ) = 0.0_wp
frq_rst_hdv( mi0(ii0,nn_hls):mi1(ii1,nn_hls) , mj0(ij0,nn_hls):mj1(ij1,nn_hls) ) = 1.e0_wp / rn_Dt
ENDIF
ENDIF
ENDIF
......
src/OCE/DOM/domzgr.F90
View file @ 3efbcb17
......@@ -130,14 +130,14 @@ CONTAINS
!
zmsk(:,:) = 1._wp ! default: no closed boundaries
IF( .NOT. l_Iperio ) THEN ! E-W closed:
zmsk( mi0( 1+nn_hls):mi1( 1+nn_hls),:) = 0._wp ! first column of inner global domain at 0
zmsk( mi0(jpiglo-nn_hls):mi1(jpiglo-nn_hls),:) = 0._wp ! last column of inner global domain at 0
zmsk( mi0( 1+nn_hls,nn_hls):mi1( 1+nn_hls,nn_hls),:) = 0._wp ! first column of inner global domain at 0
zmsk( mi0(jpiglo-nn_hls,nn_hls):mi1(jpiglo-nn_hls,nn_hls),:) = 0._wp ! last column of inner global domain at 0
ENDIF
IF( .NOT. l_Jperio ) THEN ! S closed:
zmsk(:,mj0( 1+nn_hls):mj1( 1+nn_hls) ) = 0._wp ! first line of inner global domain at 0
zmsk(:,mj0( 1+nn_hls,nn_hls):mj1( 1+nn_hls,nn_hls) ) = 0._wp ! first line of inner global domain at 0
ENDIF
IF( .NOT. ( l_Jperio .OR. l_NFold ) ) THEN ! N closed:
zmsk(:,mj0(jpjglo-nn_hls):mj1(jpjglo-nn_hls) ) = 0._wp ! last line of inner global domain at 0
zmsk(:,mj0(jpjglo-nn_hls,nn_hls):mj1(jpjglo-nn_hls,nn_hls) ) = 0._wp ! last line of inner global domain at 0
ENDIF
CALL lbc_lnk( 'usrdef_zgr', zmsk, 'T', 1. ) ! set halos
k_top(:,:) = k_top(:,:) * NINT( zmsk(:,:) )
......
src/OCE/DOM/dtatsd.F90
View file @ 3efbcb17
......@@ -161,8 +161,8 @@ CONTAINS
ij0 = 101 + nn_hls ; ij1 = 109 + nn_hls ! Reduced T & S in the Alboran Sea
ii0 = 141 + nn_hls - 1 ; ii1 = 155 + nn_hls - 1
IF( sf_tsd(jp_tem)%ln_tint .OR. irec_n(jp_tem) /= irec_b(jp_tem) ) THEN
DO jj = mj0(ij0), mj1(ij1)
DO ji = mi0(ii0), mi1(ii1)
DO jj = mj0(ij0,nn_hls), mj1(ij1,nn_hls)
DO ji = mi0(ii0,nn_hls), mi1(ii1,nn_hls)
sf_tsd(jp_tem)%fnow(ji,jj,13:13) = sf_tsd(jp_tem)%fnow(ji,jj,13:13) - 0.20_wp
sf_tsd(jp_tem)%fnow(ji,jj,14:15) = sf_tsd(jp_tem)%fnow(ji,jj,14:15) - 0.35_wp
sf_tsd(jp_tem)%fnow(ji,jj,16:25) = sf_tsd(jp_tem)%fnow(ji,jj,16:25) - 0.40_wp
......@@ -172,8 +172,8 @@ CONTAINS
ENDIF
!
IF( sf_tsd(jp_sal)%ln_tint .OR. irec_n(jp_sal) /= irec_b(jp_sal) ) THEN
DO jj = mj0(ij0), mj1(ij1)
DO ji = mi0(ii0), mi1(ii1)
DO jj = mj0(ij0,nn_hls), mj1(ij1,nn_hls)
DO ji = mi0(ii0,nn_hls), mi1(ii1,nn_hls)
sf_tsd(jp_sal)%fnow(ji,jj,13:13) = sf_tsd(jp_sal)%fnow(ji,jj,13:13) - 0.15_wp
sf_tsd(jp_sal)%fnow(ji,jj,14:15) = sf_tsd(jp_sal)%fnow(ji,jj,14:15) - 0.25_wp
sf_tsd(jp_sal)%fnow(ji,jj,16:17) = sf_tsd(jp_sal)%fnow(ji,jj,16:17) - 0.30_wp
......@@ -185,9 +185,9 @@ CONTAINS
!
ij0 = 87 + nn_hls ; ij1 = 96 + nn_hls ! Reduced temperature in Red Sea
ii0 = 148 + nn_hls - 1 ; ii1 = 160 + nn_hls - 1
sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 4:10 ) = 7.0_wp
sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 11:13 ) = 6.5_wp
sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 14:20 ) = 6.0_wp
sf_tsd(jp_tem)%fnow( mi0(ii0,nn_hls):mi1(ii1,nn_hls) , mj0(ij0,nn_hls):mj1(ij1,nn_hls) , 4:10 ) = 7.0_wp
sf_tsd(jp_tem)%fnow( mi0(ii0,nn_hls):mi1(ii1,nn_hls) , mj0(ij0,nn_hls):mj1(ij1,nn_hls) , 11:13 ) = 6.5_wp
sf_tsd(jp_tem)%fnow( mi0(ii0,nn_hls):mi1(ii1,nn_hls) , mj0(ij0,nn_hls):mj1(ij1,nn_hls) , 14:20 ) = 6.0_wp
ENDIF
ENDIF
!!gm end
......
src/OCE/DYN/dynadv.F90
View file @ 3efbcb17
......@@ -7,6 +7,7 @@ MODULE dynadv
!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase
!! 3.6 ! 2015-05 (N. Ducousso, G. Madec) add Hollingsworth scheme as an option
!! 4.0 ! 2017-07 (G. Madec) add a linear dynamics option
!! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -50,7 +51,7 @@ MODULE dynadv
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dyn_adv( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )
SUBROUTINE dyn_adv( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )
!!---------------------------------------------------------------------
!! *** ROUTINE dyn_adv ***
!!
......@@ -64,7 +65,6 @@ CONTAINS
!! (see dynvor module).
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt , Kbb, Kmm, Krhs ! ocean time step and level indices
INTEGER , OPTIONAL , INTENT(in ) :: no_zad ! no vertical advection compotation
REAL(wp), DIMENSION(:,:,:), OPTIONAL, TARGET, INTENT(in ) :: pau, pav, paw ! advective velocity
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), TARGET, INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum Eq.
!!----------------------------------------------------------------------
......@@ -73,12 +73,23 @@ CONTAINS
!
SELECT CASE( n_dynadv ) !== compute advection trend and add it to general trend ==!
CASE( np_VEC_c2 ) != vector form =!
CALL dyn_keg ( kt, nn_dynkeg , Kmm, puu, pvv, Krhs ) ! horizontal gradient of kinetic energy
CALL dyn_zad ( kt , Kmm, puu, pvv, Krhs ) ! vertical advection
! !* horizontal gradient of kinetic energy
IF (nn_hls==1) THEN ! halo 1 case
CALL dyn_keg_hls1( kt, nn_dynkeg , Kmm, puu, pvv, Krhs ) ! lbc needed with Hollingsworth scheme
ELSE ! halo 2 case
CALL dyn_keg ( kt, nn_dynkeg , Kmm, puu, pvv, Krhs )
ENDIF
CALL dyn_zad ( kt , Kmm, puu, pvv, Krhs ) !* vertical advection
!
CASE( np_FLX_c2 ) != flux form =!
CALL dyn_adv_cen2( kt , Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad ) ! 2nd order centered scheme
CASE( np_FLX_ubs )
CALL dyn_adv_ubs ( kt , Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad ) ! 3rd order UBS scheme (UP3)
CALL dyn_adv_cen2( kt , Kmm, puu, pvv, Krhs, pau, pav, paw ) !* 2nd order centered scheme
!
CASE( np_FLX_ubs ) !* 3rd order UBS scheme (UP3)
IF (nn_hls==1) THEN ! halo 1 case
CALL dyn_adv_ubs_hls1( kt , Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )
ELSE ! halo 2 case
CALL dyn_adv_ubs ( kt , Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )
ENDIF
END SELECT
!
IF( ln_timing ) CALL timing_stop( 'dyn_adv' )
......
src/OCE/DYN/dynadv_cen2.F90
View file @ 3efbcb17
......@@ -6,6 +6,7 @@ MODULE dynadv_cen2
!!======================================================================
!! History : 2.0 ! 2006-08 (G. Madec, S. Theetten) Original code
!! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option
!! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -35,7 +36,7 @@ MODULE dynadv_cen2
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dyn_adv_cen2( kt, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )
SUBROUTINE dyn_adv_cen2( kt, Kmm, puu, pvv, Krhs, pau, pav, paw )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_adv_cen2 ***
!!
......@@ -51,15 +52,17 @@ CONTAINS
!! ** Action : (puu,pvv)(:,:,:,Krhs) updated with the advective trend
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt , Kmm, Krhs ! ocean time-step and level indices
INTEGER , OPTIONAL , INTENT(in ) :: no_zad ! no vertical advection computation
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), TARGET, INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
REAL(wp), DIMENSION(:,:,:), OPTIONAL, TARGET, INTENT(in ) :: pau, pav, paw ! advective velocity
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zzu, zzv ! local scalars
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfu_t, zfu_f, zfu_uw, zfu
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfv_t, zfv_f, zfv_vw, zfv, zfw
REAL(wp) :: zzu, zzfu_kp1 ! local scalars
REAL(wp) :: zzv, zzfv_kp1 ! - -
REAL(wp), DIMENSION(A2D(1)) :: zfu_t, zfu_f, zfu
REAL(wp), DIMENSION(A2D(1)) :: zfv_t, zfv_f, zfv
REAL(wp), DIMENSION(A2D(1)) :: zfu_uw, zfv_vw, zfw
REAL(wp), DIMENSION(:,:,:) , POINTER :: zpt_u, zpt_v, zpt_w
REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: zu_trd, zv_trd
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
......@@ -71,8 +74,9 @@ CONTAINS
ENDIF
!
IF( l_trddyn ) THEN ! trends: store the input trends
zfu_uw(:,:,:) = puu(:,:,:,Krhs)
zfv_vw(:,:,:) = pvv(:,:,:,Krhs)
ALLOCATE( zu_trd(A2D(0),jpkm1), zv_trd(A2D(0),jpkm1) )
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs)
ENDIF
!
IF( PRESENT( pau ) ) THEN ! RK3: advective velocity (pau,pav,paw) /= advected velocity (puu,pvv,ww)
......@@ -89,84 +93,86 @@ CONTAINS
!
DO jk = 1, jpkm1 ! horizontal transport
DO_2D( 1, 1, 1, 1 )
zfu(ji,jj,jk) = 0.25_wp * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * zpt_u(ji,jj,jk)
zfv(ji,jj,jk) = 0.25_wp * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * zpt_v(ji,jj,jk)
zfu(ji,jj) = 0.25_wp * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * zpt_u(ji,jj,jk)
zfv(ji,jj) = 0.25_wp * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * zpt_v(ji,jj,jk)
END_2D
DO_2D( 1, 0, 1, 0 ) ! horizontal momentum fluxes (at T- and F-point)
zfu_t(ji+1,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj ,jk,Kmm) )
zfv_f(ji ,jj ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) )
zfu_f(ji ,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) )
zfv_t(ji ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji ,jj+1,jk,Kmm) )
zfu_t(ji+1,jj ) = ( zfu(ji,jj) + zfu(ji+1,jj) ) * ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj ,jk,Kmm) )
zfv_f(ji ,jj ) = ( zfv(ji,jj) + zfv(ji+1,jj) ) * ( puu(ji,jj,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) )
zfu_f(ji ,jj ) = ( zfu(ji,jj) + zfu(ji,jj+1) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) )
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,jk) - zfu_t(ji,jj ,jk) &
& + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) &
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) &
& / 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) &
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) &
& / e3v(ji,jj,jk,Kmm)
END_2D
END DO
!
IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic
zfu_uw(:,:,:) = puu(:,:,:,Krhs) - zfu_uw(:,:,:)
zfv_vw(:,:,:) = pvv(:,:,:,Krhs) - zfv_vw(:,:,:)
CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt, Kmm )
zfu_t(:,:,:) = puu(:,:,:,Krhs)
zfv_t(:,:,:) = pvv(:,:,:,Krhs)
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs) - zu_trd(A2D(0),:)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs) - zv_trd(A2D(0),:)
CALL trd_dyn( zu_trd, zv_trd, jpdyn_keg, kt, Kmm )
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs)
ENDIF
!
IF( PRESENT( no_zad ) ) THEN !== No vertical advection ==! (except if linear free surface)
! ==
IF( ln_linssh ) THEN ! linear free surface: advection through the surface z=0
DO_2D( 0, 0, 0, 0 )
zzu = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
zzv = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) - zzu * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,1,Kmm)
pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) - zzv * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,1,Kmm)
END_2D
ENDIF
!
ELSE !== Vertical advection ==!
! !== Vertical advection ==!
!
! ! surface vertical fluxes
!
IF( ln_linssh ) THEN ! linear free surface: advection through the surface z=0
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
zfv_vw(ji,jj) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
END_2D
ELSE ! non linear free: surface advective fluxes set to zero
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj) = 0._wp
zfv_vw(ji,jj) = 0._wp
END_2D
ENDIF
!
DO jk = 1, jpk-2 ! divergence of advective fluxes
!
DO_2D( 0, 0, 0, 0 ) ! surface/bottom advective fluxes set to zero
zfu_uw(ji,jj,jpk) = 0._wp ; zfv_vw(ji,jj,jpk) = 0._wp
zfu_uw(ji,jj, 1 ) = 0._wp ; zfv_vw(ji,jj, 1 ) = 0._wp
DO_2D( 0, 1, 0, 1 ) ! 1/4 * Vertical transport at level k+1
zfw(ji,jj) = 0.25_wp * e1e2t(ji,jj) * zpt_w(ji,jj,jk+1)
END_2D
IF( ln_linssh ) THEN ! linear free surface: advection through the surface z=0
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
END_2D
ENDIF
DO jk = 2, jpkm1 ! interior advective fluxes
DO_2D( 0, 1, 0, 1 ) ! 1/4 * Vertical transport
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) )
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) &
DO_2D( 0, 0, 0, 0 )
! ! vertical flux at level k+1
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,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,jk,Kmm)
END_3D
!
IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic
zfu_t(:,:,:) = puu(:,:,:,Krhs) - zfu_t(:,:,:)
zfv_t(:,:,:) = pvv(:,:,:,Krhs) - zfv_t(:,:,:)
CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt, Kmm )
ENDIF
! ! Control print
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' cen2 adv - Ua: ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
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
END_2D
END DO
!
jk = jpkm1
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)
END_2D
!
IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs) - zu_trd(A2D(0),:)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs) - zv_trd(A2D(0),:)
CALL trd_dyn( zu_trd, zv_trd, jpdyn_zad, kt, Kmm )
DEALLOCATE( zu_trd, zv_trd )
ENDIF
! ! Control print
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' cen2 adv - Ua: ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
END SUBROUTINE dyn_adv_cen2
......
src/OCE/DYN/dynadv_ubs.F90
View file @ 3efbcb17
......@@ -6,6 +6,7 @@ MODULE dynadv_ubs
!!======================================================================
!! History : 2.0 ! 2006-08 (R. Benshila, L. Debreu) Original code
!! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option
!! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -29,7 +30,8 @@ MODULE dynadv_ubs
REAL(wp), PARAMETER :: gamma1 = 1._wp/3._wp ! =1/4 quick ; =1/3 3rd order UBS
REAL(wp), PARAMETER :: gamma2 = 1._wp/32._wp ! =0 2nd order ; =1/32 4th order centred
PUBLIC dyn_adv_ubs ! routine called by step.F90
PUBLIC dyn_adv_ubs ! routine called by dynadv.F90
PUBLIC dyn_adv_ubs_hls1 ! routine called by dynadv.F90
!! * Substitutions
# include "do_loop_substitute.h90"
......@@ -41,7 +43,243 @@ MODULE dynadv_ubs
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dyn_adv_ubs( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )
SUBROUTINE dyn_adv_ubs( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_adv_ubs ***
!!
!! ** Purpose : Compute the now momentum advection trend in flux form
!! and the general trend of the momentum equation.
!!
!! ** Method : The scheme is the one implemeted in ROMS. It depends
!! on two parameter gamma1 and gamma2. The former control the
!! upstream baised part of the scheme and the later the centred
!! part: gamma1 = 0 pure centered (no diffusive part)
!! = 1/4 Quick scheme
!! = 1/3 3rd order Upstream biased scheme
!! gamma2 = 0 2nd order finite differencing
!! = 1/32 4th order finite differencing
!! For stability reasons, the first term of the fluxes which cor-
!! responds to a second order centered scheme is evaluated using
!! the now velocity (centered in time) while the second term which
!! is the diffusive part of the scheme, is evaluated using the
!! before velocity (forward in time).
!! Default value (hard coded in the begining of the module) are
!! gamma1=1/3 and gamma2=1/32.
!!
!! In RK3 time stepping case, the optional arguments
!! (pau,pav,paw) are present. They are used as advective velocity
!! while the advected velocity remains (puu,pvv).
!!
!! ** Action : (puu,pvv)(:,:,:,Krhs) updated with the advective trend
!!
!! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt , Kbb, Kmm, Krhs ! ocean time-step and level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), TARGET, INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
REAL(wp), DIMENSION(:,:,:), OPTIONAL, TARGET, INTENT(in ) :: pau, pav, paw ! advective velocity
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zzu, zui, zfuj, zl_u, zzfu_kp1 ! local scalars
REAL(wp) :: zzv, zvj, zfvi, zl_v, zzfv_kp1 ! - -
REAL(wp), DIMENSION(A2D(2)) :: zfu_t, zfu_f, zfu
REAL(wp), DIMENSION(A2D(2)) :: zfv_t, zfv_f, zfv
REAL(wp), DIMENSION(A2D(2),2) :: zlu_uu, zlu_uv
REAL(wp), DIMENSION(A2D(2),2) :: zlv_vv, zlv_vu
REAL(wp), DIMENSION(:,:,:) , POINTER :: zpt_u, zpt_v, zpt_w
REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: zu_trd, zv_trd
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
ENDIF
ENDIF
!
IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic
ALLOCATE( zu_trd(A2D(0),jpkm1), zv_trd(A2D(0),jpkm1) )
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs)
ENDIF
!
IF( PRESENT( pau ) ) THEN ! RK3: advective velocity (pau,pav,paw) /= advected velocity (puu,pvv,ww)
zpt_u => pau(:,:,:)
zpt_v => pav(:,:,:)
zpt_w => paw(:,:,:)
ELSE ! MLF: advective velocity = (puu,pvv,ww)
zpt_u => puu(:,:,:,Kmm)
zpt_v => pvv(:,:,:,Kmm)
zpt_w => ww (:,:,: )
ENDIF
!
! ! =========================== !
DO jk = 1, jpkm1 ! Laplacian of the velocity !
! ! =========================== !
! ! horizontal volume fluxes
DO_2D( 2, 2, 2, 2 )
zfu(ji,jj) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * zpt_u(ji,jj,jk)
zfv(ji,jj) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * zpt_v(ji,jj,jk)
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_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)
END_2D
!
! ! ====================== !
! ! Horizontal advection !
! ! ====================== !
! ! horizontal volume fluxes
DO_2D( 1, 1, 1, 1 )
zfu(ji,jj) = 0.25_wp * zfu(ji,jj)
zfv(ji,jj) = 0.25_wp * zfv(ji,jj)
END_2D
!
DO_2D( 1, 0, 1, 0 ) ! horizontal momentum fluxes at T- and F-point
zui = ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj ,jk,Kmm) )
zvj = ( pvv(ji,jj,jk,Kmm) + pvv(ji ,jj+1,jk,Kmm) )
!
IF( zui > 0 ) THEN ; zl_u = zlu_uu(ji ,jj,1)
ELSE ; zl_u = zlu_uu(ji+1,jj,1)
ENDIF
IF( zvj > 0 ) THEN ; zl_v = zlv_vv(ji,jj ,1)
ELSE ; zl_v = zlv_vv(ji,jj+1,1)
ENDIF
!
zfu_t(ji+1,jj ) = ( zfu(ji,jj) + zfu(ji+1,jj ) - gamma2 * ( zlu_uu(ji,jj,2) + zlu_uu(ji+1,jj ,2) ) ) &
& * ( zui - gamma1 * zl_u )
zfv_t(ji ,jj+1) = ( zfv(ji,jj) + zfv(ji ,jj+1) - gamma2 * ( zlv_vv(ji,jj,2) + zlv_vv(ji ,jj+1,2) ) ) &
& * ( zvj - gamma1 * zl_v )
!
zfuj = ( zfu(ji,jj) + zfu(ji ,jj+1) )
zfvi = ( zfv(ji,jj) + zfv(ji+1,jj ) )
IF( zfuj > 0 ) THEN ; zl_v = zlv_vu(ji ,jj,1)
ELSE ; zl_v = zlv_vu(ji+1,jj,1)
ENDIF
IF( zfvi > 0 ) THEN ; zl_u = zlu_uv(ji,jj ,1)
ELSE ; zl_u = zlu_uv(ji,jj+1,1)
ENDIF
!
zfv_f(ji ,jj ) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,2) + zlv_vu(ji+1,jj ,2) ) ) &
& * ( puu(ji,jj,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) - gamma1 * zl_u )
zfu_f(ji ,jj ) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,2) + zlu_uv(ji ,jj+1,2) ) ) &
& * ( 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) &
& / 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) &
& / e3v(ji,jj,jk,Kmm)
END_2D
END DO
!
IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs) - zu_trd(A2D(0),:)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs) - zv_trd(A2D(0),:)
CALL trd_dyn( zu_trd, zv_trd, jpdyn_keg, kt, Kmm )
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs)
ENDIF
! ! ==================== !
! ! Vertical advection !
! ! ==================== !
!
#define zfu_uw zfu_t
#define zfv_vw zfv_t
#define zfw zfu
!
! ! surface vertical fluxes
!
IF( ln_linssh ) THEN ! linear free surface: advection through the surface z=0
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
zfv_vw(ji,jj) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
END_2D
ELSE ! non linear free: surface advective fluxes set to zero
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj) = 0._wp
zfv_vw(ji,jj) = 0._wp
END_2D
ENDIF
!
DO jk = 1, jpk-2 ! divergence of advective fluxes
!
DO_2D( 0, 1, 0, 1 ) ! 1/4 * Vertical transport at level k+1
zfw(ji,jj) = 0.25_wp * e1e2t(ji,jj) * zpt_w(ji,jj,jk+1)
END_2D
DO_2D( 0, 0, 0, 0 )
! ! vertical flux at level k+1
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)
! ! store vertical flux for next level calculation
zfu_uw(ji,jj) = zzfu_kp1
zfv_vw(ji,jj) = zzfv_kp1
END_2D
END DO
!
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)
END_2D
!
IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs) - zu_trd(A2D(0),:)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs) - zv_trd(A2D(0),:)
CALL trd_dyn( zu_trd, zv_trd, jpdyn_zad, kt, Kmm )
DEALLOCATE( zu_trd, zv_trd )
ENDIF
!
#undef zfu_uw
#undef zfv_vw
#undef zfw
! ! Control print
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ubs2 adv - Ua: ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
END SUBROUTINE dyn_adv_ubs
SUBROUTINE dyn_adv_ubs_hls1( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_adv_ubs ***
!!
......@@ -73,7 +311,6 @@ CONTAINS
!! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt , Kbb, Kmm, Krhs ! ocean time-step and level indices
INTEGER , OPTIONAL , INTENT(in ) :: no_zad ! no vertical advection compotation
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), TARGET, INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
REAL(wp), DIMENSION(:,:,:), OPTIONAL, TARGET, INTENT(in ) :: pau, pav, paw ! advective velocity
!
......@@ -89,7 +326,7 @@ CONTAINS
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection'
IF(lwp) WRITE(numout,*) 'dyn_adv_ubs_hls1 : UBS flux form momentum advection'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
ENDIF
ENDIF
......@@ -230,63 +467,44 @@ CONTAINS
! ! Vertical advection !
! ! ==================== !
!
! ! ======================== !
IF( PRESENT( no_zad ) ) THEN ! No vertical advection ! (except if linear free surface)
! ! ======================== ! ------
!
IF( ln_linssh ) THEN ! linear free surface: advection through the surface z=0
DO_2D( 0, 0, 0, 0 )
zzu = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
zzv = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) - zzu * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,1,Kmm)
pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) - zzv * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,1,Kmm)
END_2D
ENDIF
! ! =================== !
ELSE ! Vertical advection !
! ! =================== !
DO_2D( 0, 0, 0, 0 ) ! surface/bottom advective fluxes set to zero
zfu_uw(ji,jj,jpk) = 0._wp
zfv_vw(ji,jj,jpk) = 0._wp
zfu_uw(ji,jj, 1 ) = 0._wp
zfv_vw(ji,jj, 1 ) = 0._wp
DO_2D( 0, 0, 0, 0 ) ! surface/bottom advective fluxes set to zero
zfu_uw(ji,jj,jpk) = 0._wp
zfv_vw(ji,jj,jpk) = 0._wp
zfu_uw(ji,jj, 1 ) = 0._wp
zfv_vw(ji,jj, 1 ) = 0._wp
END_2D
IF( ln_linssh ) THEN ! constant volume : advection through the surface
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
END_2D
IF( ln_linssh ) THEN ! constant volume : advection through the surface
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
END_2D
ENDIF
DO jk = 2, jpkm1 ! interior fluxes
DO_2D( 0, 1, 0, 1 )
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) )
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)
END_3D
!
IF( l_trddyn ) THEN ! save the vertical advection trend for diagnostic
zfu_t(:,:,:) = puu(:,:,:,Krhs) - zfu_t(:,:,:)
zfv_t(:,:,:) = pvv(:,:,:,Krhs) - zfv_t(:,:,:)
CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt, Kmm )
ENDIF
! ! Control print
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ubs2 adv - Ua: ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
ENDIF
DO jk = 2, jpkm1 ! interior fluxes
DO_2D( 0, 1, 0, 1 )
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) )
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)
END_3D
!
END SUBROUTINE dyn_adv_ubs
IF( l_trddyn ) THEN ! save the vertical advection trend for diagnostic
zfu_t(:,:,:) = puu(:,:,:,Krhs) - zfu_t(:,:,:)
zfv_t(:,:,:) = pvv(:,:,:,Krhs) - zfv_t(:,:,:)
CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt, Kmm )
ENDIF
! ! Control print
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ubs2 adv - Ua: ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
END SUBROUTINE dyn_adv_ubs_hls1
!!==============================================================================
END MODULE dynadv_ubs
src/OCE/DYN/dynatf.F90
View file @ 3efbcb17
......@@ -328,29 +328,29 @@ CONTAINS
!
IF ( iom_use("utau") ) THEN
IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
ALLOCATE(zutau(jpi,jpj))
ALLOCATE(zutau(A2D(0)))
DO_2D( 0, 0, 0, 0 )
jk = miku(ji,jj)
zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * puu(ji,jj,jk,Kaa)
zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * rCdU_top(ji,jj) * ( puu(ji-1,jj,jk,Kaa) + puu(ji,jj,jk,Kaa) )
END_2D
CALL iom_put( "utau", zutau(:,:) )
DEALLOCATE(zutau)
ELSE
CALL iom_put( "utau", utau(:,:) )
CALL iom_put( "utau", utau(A2D(0)) )
ENDIF
ENDIF
!
IF ( iom_use("vtau") ) THEN
IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
ALLOCATE(zvtau(jpi,jpj))
ALLOCATE(zvtau(A2D(0)))
DO_2D( 0, 0, 0, 0 )
jk = mikv(ji,jj)
zvtau(ji,jj) = vtau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * pvv(ji,jj,jk,Kaa)
zvtau(ji,jj) = vtau(ji,jj) + 0.5_wp * rho0 * rCdU_top(ji,jj) * ( pvv(ji,jj-1,jk,Kaa) + pvv(ji,jj,jk,Kaa) )
END_2D
CALL iom_put( "vtau", zvtau(:,:) )
DEALLOCATE(zvtau)
ELSE
CALL iom_put( "vtau", vtau(:,:) )
CALL iom_put( "vtau", vtau(A2D(0)) )
ENDIF
ENDIF
!
......
src/OCE/DYN/dynatf_qco.F90
View file @ 3efbcb17
......@@ -245,29 +245,29 @@ CONTAINS
!
IF ( iom_use("utau") ) THEN
IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
ALLOCATE(zutau(jpi,jpj))
ALLOCATE(zutau(A2D(0)))
DO_2D( 0, 0, 0, 0 )
jk = miku(ji,jj)
zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * puu(ji,jj,jk,Kaa)
zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * rCdU_top(ji,jj) * ( puu(ji-1,jj,jk,Kaa) + puu(ji,jj,jk,Kaa) )
END_2D
CALL iom_put( "utau", zutau(:,:) )
DEALLOCATE(zutau)
ELSE
CALL iom_put( "utau", utau(:,:) )
CALL iom_put( "utau", utau(A2D(0)) )
ENDIF
ENDIF
!
IF ( iom_use("vtau") ) THEN
IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
ALLOCATE(zvtau(jpi,jpj))
ALLOCATE(zvtau(A2D(0)))
DO_2D( 0, 0, 0, 0 )
jk = mikv(ji,jj)
zvtau(ji,jj) = vtau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * pvv(ji,jj,jk,Kaa)
zvtau(ji,jj) = vtau(ji,jj) + 0.5_wp * rho0 * rCdU_top(ji,jj) * ( pvv(ji,jj-1,jk,Kaa) + pvv(ji,jj,jk,Kaa) )
END_2D
CALL iom_put( "vtau", zvtau(:,:) )
DEALLOCATE(zvtau)
ELSE
CALL iom_put( "vtau", vtau(:,:) )
CALL iom_put( "vtau", vtau(A2D(0)) )
ENDIF
ENDIF
!
......
src/OCE/DYN/dynkeg.F90
View file @ 3efbcb17
......@@ -6,7 +6,8 @@ MODULE dynkeg
!! History : 1.0 ! 1987-09 (P. Andrich, M.-A. Foujols) Original code
!! 7.0 ! 1997-05 (G. Madec) Split dynber into dynkeg and dynhpg
!! NEMO 1.0 ! 2002-07 (G. Madec) F90: Free form and module
!! 3.6 ! 2015-05 (N. Ducousso, G. Madec) add Hollingsworth scheme as an option
!! 3.6 ! 2015-05 (N. Ducousso, G. Madec) add Hollingsworth scheme as an option
!! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -27,7 +28,8 @@ MODULE dynkeg
IMPLICIT NONE
PRIVATE
PUBLIC dyn_keg ! routine called by step module
PUBLIC dyn_keg ! routine called by step module
PUBLIC dyn_keg_hls1 ! routine called by step module
INTEGER, PARAMETER, PUBLIC :: nkeg_C2 = 0 !: 2nd order centered scheme (standard scheme)
INTEGER, PARAMETER, PUBLIC :: nkeg_HW = 1 !: Hollingsworth et al., QJRMS, 1983
......@@ -44,6 +46,123 @@ MODULE dynkeg
CONTAINS
SUBROUTINE dyn_keg( kt, kscheme, Kmm, puu, pvv, Krhs )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_keg ***
!!
!! ** Purpose : Compute the now momentum trend due to the horizontal
!! gradient of the horizontal kinetic energy and add it to the
!! general momentum trend.
!!
!! ** Method : * kscheme = nkeg_C2 : 2nd order centered scheme that
!! conserve kinetic energy. Compute the now horizontal kinetic energy
!! zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ]
!! * kscheme = nkeg_HW : Hollingsworth correction following
!! Arakawa (2001). The now horizontal kinetic energy is given by:
!! zhke = 1/6 [ mi-1( 2 * un^2 + ((u(j+1)+u(j-1))/2)^2 )
!! + mj-1( 2 * vn^2 + ((v(i+1)+v(i-1))/2)^2 ) ]
!!
!! Take its horizontal gradient and add it to the general momentum
!! trend.
!! u(rhs) = u(rhs) - 1/e1u di[ zhke ]
!! v(rhs) = v(rhs) - 1/e2v dj[ zhke ]
!!
!! ** Action : - Update the (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) with the hor. ke gradient trend
!! - send this trends to trd_dyn (l_trddyn=T) for post-processing
!!
!! ** References : Arakawa, A., International Geophysics 2001.
!! Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983.
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: kscheme ! =0/1 type of KEG scheme
INTEGER , INTENT(in ) :: Kmm, Krhs ! ocean time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zu, zv ! local scalars
REAL(wp), DIMENSION(:,: ) , ALLOCATABLE :: zhke
REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: zu_trd, zv_trd
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('dyn_keg')
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme
IF(lwp) WRITE(numout,*) '~~~~~~~'
ENDIF
ENDIF
!
IF( l_trddyn ) THEN ! Save the input trends
ALLOCATE( zu_trd(A2D(0),jpk), zv_trd(A2D(0),jpk) )
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs)
ENDIF
!
SELECT CASE ( kscheme )
!
CASE ( nkeg_C2 ) !== Standard scheme ==!
ALLOCATE( zhke(A2D(1)) )
DO jk = 1, jpkm1
DO_2D( 0, 1, 0, 1 ) !* Horizontal kinetic energy at T-point
zu = puu(ji-1,jj ,jk,Kmm) * puu(ji-1,jj ,jk,Kmm) &
& + puu(ji ,jj ,jk,Kmm) * puu(ji ,jj ,jk,Kmm)
zv = pvv(ji ,jj-1,jk,Kmm) * pvv(ji ,jj-1,jk,Kmm) &
& + pvv(ji ,jj ,jk,Kmm) * pvv(ji ,jj ,jk,Kmm)
zhke(ji,jj) = 0.25_wp * ( zv + zu )
END_2D
!
DO_2D( 0, 0, 0, 0 ) !* grad( KE ) added to the general momentum trends
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zhke(ji+1,jj ) - zhke(ji,jj) ) * r1_e1u(ji,jj)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zhke(ji ,jj+1) - zhke(ji,jj) ) * r1_e2v(ji,jj)
END_2D
END DO
DEALLOCATE( zhke )
!
CASE ( nkeg_HW ) !* Hollingsworth scheme
ALLOCATE( zhke(A2D(1)) )
DO jk = 1, jpkm1
DO_2D( 0, 1, 0, 1 )
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
zu = ( puu(ji-1,jj ,jk,Kmm) * puu(ji-1,jj ,jk,Kmm) &
& + puu(ji ,jj ,jk,Kmm) * puu(ji ,jj ,jk,Kmm) ) * 8._wp &
& + ( ( puu(ji-1,jj-1,jk,Kmm) + puu(ji-1,jj+1,jk,Kmm) ) * ( puu(ji-1,jj-1,jk,Kmm) + puu(ji-1,jj+1,jk,Kmm) ) &
& + ( puu(ji ,jj-1,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) ) * ( puu(ji ,jj-1,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) ) &
& ) ! bracket for halo 1 - halo 2 compatibility
zv = ( pvv(ji ,jj-1,jk,Kmm) * pvv(ji ,jj-1,jk,Kmm) &
& + pvv(ji ,jj ,jk,Kmm) * pvv(ji ,jj ,jk,Kmm) ) * 8._wp &
& + ( ( pvv(ji-1,jj-1,jk,Kmm) + pvv(ji+1,jj-1,jk,Kmm) ) * ( pvv(ji-1,jj-1,jk,Kmm) + pvv(ji+1,jj-1,jk,Kmm) ) &
& + ( pvv(ji-1,jj ,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) ) * ( pvv(ji-1,jj ,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) ) &
& ) ! bracket for halo 1 - halo 2 compatibility
zhke(ji,jj) = r1_48 * ( zv + zu )
END_2D
!
DO_2D( 0, 0, 0, 0 ) !* grad( KE ) added to the general momentum trends
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zhke(ji+1,jj ) - zhke(ji,jj) ) * r1_e1u(ji,jj)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zhke(ji ,jj+1) - zhke(ji,jj) ) * r1_e2v(ji,jj)
END_2D
END DO
DEALLOCATE( zhke )
!
END SELECT
!
IF( l_trddyn ) THEN ! save the Kinetic Energy trends for diagnostic
zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs) - zu_trd(A2D(0),:)
zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs) - zv_trd(A2D(0),:)
CALL trd_dyn( zu_trd, zv_trd, jpdyn_keg, kt, Kmm )
DEALLOCATE( zu_trd, zv_trd )
ENDIF
!
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' keg - Ua: ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
IF( ln_timing ) CALL timing_stop('dyn_keg')
!
END SUBROUTINE dyn_keg
SUBROUTINE dyn_keg_hls1( kt, kscheme, Kmm, puu, pvv, Krhs )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_keg ***
!!
......@@ -86,7 +205,7 @@ CONTAINS
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme
IF(lwp) WRITE(numout,*) 'dyn_keg_hls1 : kinetic energy gradient trend, scheme number=', kscheme
IF(lwp) WRITE(numout,*) '~~~~~~~'
ENDIF
ENDIF
......@@ -147,7 +266,7 @@ CONTAINS
!
IF( ln_timing ) CALL timing_stop('dyn_keg')
!
END SUBROUTINE dyn_keg
END SUBROUTINE dyn_keg_hls1
!!======================================================================
END MODULE dynkeg
src/OCE/DYN/dynspg_ts.F90
View file @ 3efbcb17
......@@ -334,14 +334,14 @@ CONTAINS
! ! ------------------ !
IF( ln_bt_fw ) THEN
DO_2D( 0, 0, 0, 0 )
zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kmm)
zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rho0 * vtau(ji,jj) * r1_hv(ji,jj,Kmm)
zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rho0 * utauU(ji,jj) * r1_hu(ji,jj,Kmm)
zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rho0 * vtauV(ji,jj) * r1_hv(ji,jj,Kmm)
END_2D
ELSE
zztmp = r1_rho0 * r1_2
DO_2D( 0, 0, 0, 0 )
zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu(ji,jj,Kmm)
zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kmm)
zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utauU(ji,jj) ) * r1_hu(ji,jj,Kmm)
zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtauV(ji,jj) ) * r1_hv(ji,jj,Kmm)
END_2D
ENDIF
!
......
src/OCE/DYN/dynvor.F90
View file @ 3efbcb17
......@@ -22,6 +22,7 @@ MODULE dynvor
!! - ! 2018-04 (G. Madec) add pre-computed gradient for metric term calculation
!! 4.x ! 2020-03 (G. Madec, A. Nasser) make ln_dynvor_msk truly efficient on relative vorticity
!! 4.2 ! 2020-12 (G. Madec, E. Clementi) add vortex force trends (ln_vortex_force=T)
!! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -172,11 +173,20 @@ CONTAINS
CALL vor_enT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
ENDIF
CASE( np_EET ) !* energy conserving scheme (een scheme using e3t)
IF( nn_hls == 1 ) THEN
CALL vor_eeT_hls1( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN
CALL vor_eeT_hls1( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN
CALL vor_eeT_hls1( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
ENDIF
ELSE
CALL vor_eeT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN
IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN
CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN
ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN
CALL vor_eeT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
ENDIF
ENDIF
CASE( np_ENE ) !* energy conserving scheme
CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
......@@ -199,11 +209,20 @@ CONTAINS
IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
IF( ln_vortex_force ) CALL vor_ens( kt, Kmm, nrvm, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add vortex force
CASE( np_EEN ) !* energy and enstrophy conserving scheme
IF( nn_hls == 1 ) THEN
CALL vor_een_hls1( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN
CALL vor_een_hls1( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN
CALL vor_een_hls1( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
ENDIF
ELSE
CALL vor_een( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN
IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN
CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN
ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN
CALL vor_een( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
ENDIF
ENDIF
END SELECT
!
......@@ -320,7 +339,122 @@ CONTAINS
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars
REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx, zwy, zwt ! 2D workspace
REAL(wp), DIMENSION(A2D(1)) :: zwx, zwy, zwt ! 2D workspace
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zwz ! 3D workspace, jpkm1 -> avoid lbc_lnk on jpk that is not defined
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'dyn:vor_enT : vorticity term: t-point energy conserving scheme'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
ENDIF
ENDIF
!
!
! ! ===============
DO jk = 1, jpkm1 ! Horizontal slab
! ! ===============
!
SELECT CASE( kvor ) !== relative vorticity considered ==!
!
CASE ( np_RVO , np_CRV ) !* relative vorticity at f-point is used
ALLOCATE( zwz(A2D(1)) )
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
& - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)
END_2D
IF( ln_dynvor_msk ) THEN ! mask relative vorticity
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
END_2D
ENDIF
!
END SELECT
!
SELECT CASE( kvor ) !== volume weighted vorticity considered ==!
!
CASE ( np_COR ) !* Coriolis (planetary vorticity)
DO_2D( 0, 1, 0, 1 )
zwt(ji,jj) = ff_t(ji,jj) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
END_2D
CASE ( np_RVO ) !* relative vorticity
DO_2D( 0, 1, 0, 1 )
zwt(ji,jj) = r1_4 * ( zwz(ji-1,jj ) + zwz(ji,jj ) &
& + zwz(ji-1,jj-1) + zwz(ji,jj-1) ) &
& * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
END_2D
CASE ( np_MET ) !* metric term
DO_2D( 0, 1, 0, 1 )
zwt(ji,jj) = ( ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) &
& - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) &
& * e3t(ji,jj,jk,Kmm)
END_2D
CASE ( np_CRV ) !* Coriolis + relative vorticity
DO_2D( 0, 1, 0, 1 )
zwt(ji,jj) = ( ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj ) + zwz(ji,jj ) &
& + zwz(ji-1,jj-1) + zwz(ji,jj-1) ) ) &
& * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
END_2D
CASE ( np_CME ) !* Coriolis + metric
DO_2D( 0, 1, 0, 1 )
zwt(ji,jj) = ( ff_t(ji,jj) * e1e2t(ji,jj) &
& + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) &
& - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) &
& * e3t(ji,jj,jk,Kmm)
END_2D
CASE DEFAULT ! error
CALL ctl_stop('STOP','dyn_vor: wrong value for kvor')
END SELECT
!
! !== compute and add the vorticity term trend =!
DO_2D( 0, 0, 0, 0 )
pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) &
& * ( zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) ) &
& + zwt(ji ,jj) * ( pv(ji ,jj,jk) + pv(ji ,jj-1,jk) ) )
!
pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) &
& * ( zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) ) &
& + zwt(ji,jj ) * ( pu(ji,jj ,jk) + pu(ji-1,jj ,jk) ) )
END_2D
! ! ===============
END DO ! End of slab
! ! ===============
!
SELECT CASE( kvor ) ! deallocate zwz if necessary
CASE ( np_RVO , np_CRV ) ; DEALLOCATE( zwz )
END SELECT
!
END SUBROUTINE vor_enT
SUBROUTINE vor_enT_hls1( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
!!----------------------------------------------------------------------
!! *** ROUTINE vor_enT ***
!!
!! ** Purpose : Compute the now total vorticity trend and add it to
!! the general trend of the momentum equation.
!!
!! ** Method : Trend evaluated using now fields (centered in time)
!! and t-point evaluation of vorticity (planetary and relative).
!! conserves the horizontal kinetic energy.
!! The general trend of momentum is increased due to the vorticity
!! term which is given by:
!! voru = 1/bu mj[ ( mi(mj(bf*rvor))+bt*f_t)/e3t mj[vn] ]
!! vorv = 1/bv mi[ ( mi(mj(bf*rvor))+bt*f_t)/e3f mj[un] ]
!! where rvor is the relative vorticity at f-point
!!
!! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kmm ! ocean time level index
INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars
REAL(wp), DIMENSION(A2D(1)) :: zwx, zwy, zwt ! 2D workspace
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zwz ! 3D workspace, jpkm1 -> avoid lbc_lnk on jpk that is not defined
!!----------------------------------------------------------------------
!
......@@ -336,7 +470,7 @@ CONTAINS
SELECT CASE( kvor ) !== relative vorticity considered ==!
!
CASE ( np_RVO , np_CRV ) !* relative vorticity at f-point is used
ALLOCATE( zwz(A2D(nn_hls),jpk) )
ALLOCATE( zwz(A2D(1),jpk) )
DO jk = 1, jpkm1 ! Horizontal slab
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
zwz(ji,jj,jk) = ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
......@@ -409,7 +543,7 @@ CONTAINS
CASE ( np_RVO , np_CRV ) ; DEALLOCATE( zwz )
END SELECT
!
END SUBROUTINE vor_enT
END SUBROUTINE vor_enT_hls1
SUBROUTINE vor_ene( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
......@@ -440,7 +574,7 @@ CONTAINS
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zx1, zy1, zx2, zy2, ze3f, zmsk ! local scalars
REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx, zwy, zwz ! 2D workspace
REAL(wp), DIMENSION(A2D(1)) :: zwx, zwy, zwz ! 2D workspace
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
......@@ -573,7 +707,7 @@ CONTAINS
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zuav, zvau, ze3f, zmsk ! local scalars
REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx, zwy, zwz ! 2D workspace
REAL(wp), DIMENSION(A2D(1)) :: zwx, zwy, zwz ! 2D workspace
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
......@@ -705,16 +839,176 @@ CONTAINS
INTEGER :: ierr ! local integer
REAL(wp) :: zua, zva ! local scalars
REAL(wp) :: zmsk, ze3f ! local scalars
REAL(wp), DIMENSION(A2D(nn_hls)) :: z1_e3f
REAL(wp), DIMENSION(A2D(1)) :: z1_e3f
#if defined key_loop_fusion
REAL(wp) :: ztne, ztnw, ztnw_ip1, ztse, ztse_jp1, ztsw_jp1, ztsw_ip1
REAL(wp) :: zwx, zwx_im1, zwx_jp1, zwx_im1_jp1
REAL(wp) :: zwy, zwy_ip1, zwy_jm1, zwy_ip1_jm1
#else
REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx , zwy
REAL(wp), DIMENSION(A2D(nn_hls)) :: ztnw, ztne, ztsw, ztse
REAL(wp), DIMENSION(A2D(1)) :: zwx , zwy
REAL(wp), DIMENSION(A2D(1)) :: ztnw, ztne, ztsw, ztse
#endif
REAL(wp), DIMENSION(A2D(nn_hls),jpkm1) :: zwz ! 3D workspace, jpkm1 -> jpkm1 -> avoid lbc_lnk on jpk that is not defined
REAL(wp), DIMENSION(A2D(1)) :: zwz ! 3D workspace, jpkm1 -> jpkm1 -> avoid lbc_lnk on jpk that is not defined
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
ENDIF
ENDIF
!
! ! ===============
DO jk = 1, jpkm1 ! Horizontal slab
! ! ===============
!
#if defined key_qco || defined key_linssh
DO_2D( 1, 1, 1, 1 ) ! == reciprocal of e3 at F-point (key_qco)
z1_e3f(ji,jj) = 1._wp / e3f_vor(ji,jj,jk)
END_2D
#else
SELECT CASE( nn_e3f_typ ) ! == reciprocal of e3 at F-point
CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
DO_2D( 1, 1, 1, 1 )
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
ze3f = ( (e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) &
& + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)) &
& + (e3t(ji ,jj ,jk,Kmm)*tmask(ji ,jj ,jk) &
& + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk)) )
IF( ze3f /= 0._wp ) THEN ; z1_e3f(ji,jj) = 4._wp / ze3f
ELSE ; z1_e3f(ji,jj) = 0._wp
ENDIF
END_2D
CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
DO_2D( 1, 1, 1, 1 )
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
ze3f = ( (e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) &
& + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)) &
& + (e3t(ji ,jj ,jk,Kmm)*tmask(ji ,jj ,jk) &
& + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk)) )
zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &
& + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) )
IF( ze3f /= 0._wp ) THEN ; z1_e3f(ji,jj) = zmsk / ze3f
ELSE ; z1_e3f(ji,jj) = 0._wp
ENDIF
END_2D
END SELECT
#endif
!
SELECT CASE( kvor ) !== vorticity considered ==!
!
CASE ( np_COR ) !* Coriolis (planetary vorticity)
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = ff_f(ji,jj) * z1_e3f(ji,jj)
END_2D
CASE ( np_RVO ) !* relative vorticity
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = ( e2v(ji+1,jj ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
& - e1u(ji ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)
END_2D
IF( ln_dynvor_msk ) THEN ! mask the relative vorticity
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
END_2D
ENDIF
CASE ( np_MET ) !* metric term
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = ( ( pv(ji+1,jj,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
& - ( pu(ji,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj)
END_2D
CASE ( np_CRV ) !* Coriolis + relative vorticity
DO_2D( 1, 1, 1, 1 )
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
zwz(ji,jj) = ( ff_f(ji,jj) + ( ( e2v(ji+1,jj ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
& ) & ! bracket for halo 1 - halo 2 compatibility
& - ( e1u(ji ,jj+1) * pu(ji,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk) &
& ) & ! bracket for halo 1 - halo 2 compatibility
& ) * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj)
END_2D
IF( ln_dynvor_msk ) THEN ! mask the relative vorticity
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
END_2D
ENDIF
CASE ( np_CME ) !* Coriolis + metric
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = ( ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj)
END_2D
CASE DEFAULT ! error
CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' )
END SELECT
!
! !== horizontal fluxes ==!
DO_2D( 1, 1, 1, 1 )
zwx(ji,jj) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * pu(ji,jj,jk)
zwy(ji,jj) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pv(ji,jj,jk)
END_2D
!
! !== compute and add the vorticity term trend =!
DO_2D( 0, 1, 0, 1 )
ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
END_2D
!
DO_2D( 0, 0, 0, 0 )
zua = + r1_12 * r1_e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) &
& + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
zva = - r1_12 * r1_e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) &
& + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) )
pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
END_2D
END DO
! ! ===============
! ! End of slab
! ! ===============
END SUBROUTINE vor_een
SUBROUTINE vor_een_hls1( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
!!----------------------------------------------------------------------
!! *** ROUTINE vor_een ***
!!
!! ** Purpose : Compute the now total vorticity trend and add it to
!! the general trend of the momentum equation.
!!
!! ** Method : Trend evaluated using now fields (centered in time)
!! and the Arakawa and Lamb (1980) flux form formulation : conserves
!! both the horizontal kinetic energy and the potential enstrophy
!! when horizontal divergence is zero (see the NEMO documentation)
!! Add this trend to the general momentum trend (pu_rhs,pv_rhs).
!!
!! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
!!
!! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kmm ! ocean time level index
INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
!
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ierr ! local integer
REAL(wp) :: zua, zva ! local scalars
REAL(wp) :: zmsk, ze3f ! local scalars
REAL(wp), DIMENSION(A2D(1)) :: z1_e3f
#if defined key_loop_fusion
REAL(wp) :: ztne, ztnw, ztnw_ip1, ztse, ztse_jp1, ztsw_jp1, ztsw_ip1
REAL(wp) :: zwx, zwx_im1, zwx_jp1, zwx_im1_jp1
REAL(wp) :: zwy, zwy_ip1, zwy_jm1, zwy_ip1_jm1
#else
REAL(wp), DIMENSION(A2D(1)) :: zwx , zwy
REAL(wp), DIMENSION(A2D(1)) :: ztnw, ztne, ztsw, ztse
#endif
REAL(wp), DIMENSION(A2D(1),jpkm1) :: zwz ! 3D workspace, jpkm1 -> jpkm1 -> avoid lbc_lnk on jpk that is not defined
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
......@@ -874,7 +1168,7 @@ CONTAINS
! ! ===============
! ! End of slab
! ! ===============
END SUBROUTINE vor_een
END SUBROUTINE vor_een_hls1
SUBROUTINE vor_eeT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
......@@ -904,9 +1198,129 @@ CONTAINS
INTEGER :: ierr ! local integer
REAL(wp) :: zua, zva ! local scalars
REAL(wp) :: zmsk, z1_e3t ! local scalars
REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx , zwy
REAL(wp), DIMENSION(A2D(nn_hls)) :: ztnw, ztne, ztsw, ztse
REAL(wp), DIMENSION(A2D(nn_hls),jpkm1) :: zwz ! 3D workspace, avoid lbc_lnk on jpk that is not defined
REAL(wp), DIMENSION(A2D(1)) :: zwx , zwy
REAL(wp), DIMENSION(A2D(1)) :: ztnw, ztne, ztsw, ztse
REAL(wp), DIMENSION(A2D(1)) :: zwz ! 3D workspace, avoid lbc_lnk on jpk that is not defined
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'dyn:vor_eeT : vorticity term: energy and enstrophy conserving scheme'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
ENDIF
ENDIF
!
! ! ===============
DO jk = 1, jpkm1 ! Horizontal slab
! ! ===============
!
!
SELECT CASE( kvor ) !== vorticity considered ==!
CASE ( np_COR ) !* Coriolis (planetary vorticity)
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = ff_f(ji,jj)
END_2D
CASE ( np_RVO ) !* relative vorticity
DO_2D( 1, 1, 1, 1 )
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
zwz(ji,jj) = ( (e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk)) &
& - (e1u(ji ,jj+1) * pu(ji ,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk)) ) &
& * r1_e1e2f(ji,jj)
END_2D
IF( ln_dynvor_msk ) THEN ! mask the relative vorticity
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
END_2D
ENDIF
CASE ( np_MET ) !* metric term
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
END_2D
CASE ( np_CRV ) !* Coriolis + relative vorticity
DO_2D( 1, 1, 1, 1 )
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
zwz(ji,jj) = ( ff_f(ji,jj) + ( (e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk)) &
& - (e1u(ji ,jj+1) * pu(ji ,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk)) ) &
& * r1_e1e2f(ji,jj) )
END_2D
IF( ln_dynvor_msk ) THEN ! mask the relative vorticity
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
END_2D
ENDIF
CASE ( np_CME ) !* Coriolis + metric
DO_2D( 1, 1, 1, 1 )
zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
& - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
END_2D
CASE DEFAULT ! error
CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' )
END SELECT
!
!
! !== horizontal fluxes ==!
DO_2D( 1, 1, 1, 1 )
zwx(ji,jj) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * pu(ji,jj,jk)
zwy(ji,jj) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pv(ji,jj,jk)
END_2D
!
! !== compute and add the vorticity term trend =!
DO_2D( 0, 1, 0, 1 )
z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
ztne(ji,jj) = ( zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) ) * z1_e3t
ztnw(ji,jj) = ( zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) ) * z1_e3t
ztse(ji,jj) = ( zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) ) * z1_e3t
ztsw(ji,jj) = ( zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) ) * z1_e3t
END_2D
!
DO_2D( 0, 0, 0, 0 )
zua = + r1_12 * r1_e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) &
& + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
zva = - r1_12 * r1_e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) &
& + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) )
pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
END_2D
! ! ===============
END DO ! End of slab
! ! ===============
END SUBROUTINE vor_eeT
SUBROUTINE vor_eeT_hls1( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
!!----------------------------------------------------------------------
!! *** ROUTINE vor_eeT ***
!!
!! ** Purpose : Compute the now total vorticity trend and add it to
!! the general trend of the momentum equation.
!!
!! ** Method : Trend evaluated using now fields (centered in time)
!! and the Arakawa and Lamb (1980) vector form formulation using
!! a modified version of Arakawa and Lamb (1980) scheme (see vor_een).
!! The change consists in
!! Add this trend to the general momentum trend (pu_rhs,pv_rhs).
!!
!! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
!!
!! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kmm ! ocean time level index
INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
!
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ierr ! local integer
REAL(wp) :: zua, zva ! local scalars
REAL(wp) :: zmsk, z1_e3t ! local scalars
REAL(wp), DIMENSION(A2D(1)) :: zwx , zwy
REAL(wp), DIMENSION(A2D(1)) :: ztnw, ztne, ztsw, ztse
REAL(wp), DIMENSION(A2D(1),jpkm1) :: zwz ! 3D workspace, avoid lbc_lnk on jpk that is not defined
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
......@@ -1003,7 +1417,7 @@ CONTAINS
! ! ===============
END DO ! End of slab
! ! ===============
END SUBROUTINE vor_eeT
END SUBROUTINE vor_eeT_hls1
SUBROUTINE dyn_vor_init
......
src/OCE/DYN/dynzdf.F90
View file @ 3efbcb17
......@@ -6,6 +6,7 @@ MODULE dynzdf
!! History : 1.0 ! 2005-11 (G. Madec) Original code
!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase
!! 4.0 ! 2017-06 (G. Madec) remove the explicit time-stepping option + avm at t-point
!! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -79,7 +80,7 @@ CONTAINS
REAL(wp) :: zWu , zWv ! - -
REAL(wp) :: zWui, zWvi ! - -
REAL(wp) :: zWus, zWvs ! - -
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwi, zwd, zws ! 3D workspace
REAL(wp), DIMENSION(A1Di(0),jpk) :: zwi, zwd, zws ! 2D workspace
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv ! - -
!!---------------------------------------------------------------------
!
......@@ -105,315 +106,329 @@ CONTAINS
ztrdv(:,:,:) = pvv(:,:,:,Krhs)
ENDIF
!
! !== RHS: Leap-Frog time stepping on all trends but the vertical mixing ==! (put in puu(:,:,:,Kaa),pvv(:,:,:,Kaa))
!
! ! time stepping except vertical diffusion
IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kbb) + rDt * puu(ji,jj,jk,Krhs) ) * umask(ji,jj,jk)
pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kbb) + rDt * pvv(ji,jj,jk,Krhs) ) * vmask(ji,jj,jk)
END_3D
ELSE ! applied on thickness weighted velocity
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
puu(ji,jj,jk,Kaa) = ( e3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb ) &
& + rDt * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Krhs) ) &
& / e3u(ji,jj,jk,Kaa) * umask(ji,jj,jk)
pvv(ji,jj,jk,Kaa) = ( e3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb ) &
& + rDt * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Krhs) ) &
& / e3v(ji,jj,jk,Kaa) * vmask(ji,jj,jk)
END_3D
ENDIF
! ! add top/bottom friction
! With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom.
! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does
! not lead to the effective stress seen over the whole barotropic loop.
! G. Madec : in linear free surface, e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) = e3u_0, so systematic use of e3u(:,:,:,Kaa)
IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! remove barotropic velocities
puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - uu_b(ji,jj,Kaa) ) * umask(ji,jj,jk)
pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - vv_b(ji,jj,Kaa) ) * vmask(ji,jj,jk)
END_3D
DO_2D( 0, 0, 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) &
& / 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) &
& / e3v(ji,jj,ikv,Kaa)
END_2D
IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! Ocean cavities (ISF)
DO_2D( 0, 0, 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) &
! ! ================= !
DO_1Dj( 0, 0 ) ! i-k slices loop !
! ! ================= !
!
! !== RHS: Leap-Frog time stepping on all trends but the vertical mixing ==! (put in puu(:,:,:,Kaa),pvv(:,:,:,Kaa))
!
! ! time stepping except vertical diffusion
IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity
DO_2Dik( 0, 0, 1, jpkm1, 1 )
puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kbb) + rDt * puu(ji,jj,jk,Krhs) ) * umask(ji,jj,jk)
pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kbb) + rDt * pvv(ji,jj,jk,Krhs) ) * vmask(ji,jj,jk)
END_2D
ELSE ! applied on thickness weighted velocity
DO_2Dik( 0, 0, 1, jpkm1, 1 )
puu(ji,jj,jk,Kaa) = ( e3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb ) &
& + rDt * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Krhs) ) &
& / e3u(ji,jj,jk,Kaa) * umask(ji,jj,jk)
pvv(ji,jj,jk,Kaa) = ( e3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb ) &
& + rDt * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Krhs) ) &
& / e3v(ji,jj,jk,Kaa) * vmask(ji,jj,jk)
END_2D
ENDIF
! ! add top/bottom friction
! With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom.
! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does
! not lead to the effective stress seen over the whole barotropic loop.
! G. Madec : in linear free surface, e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) = e3u_0, so systematic use of e3u(:,:,:,Kaa)
IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
DO_2Dik( 0, 0, 1, jpkm1, 1 ) ! remove barotropic velocities
puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - uu_b(ji,jj,Kaa) ) * umask(ji,jj,jk)
pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - vv_b(ji,jj,Kaa) ) * vmask(ji,jj,jk)
END_2D
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) &
& / 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_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) &
& / e3v(ji,jj,ikv,Kaa)
END_2D
END IF
ENDIF
!
! !== Vertical diffusion on u ==!
!
! !* Matrix construction
IF( ln_zad_Aimp ) THEN !!
SELECT CASE( nldf_dyn )
CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu)
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
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)
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,jj,jk) = zzwi + zDt_2 * MIN( zWui, 0._wp )
zws(ji,jj,jk) = zzws - zDt_2 * MAX( zWus, 0._wp )
zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zDt_2 * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
END_3D
CASE DEFAULT ! iso-level lateral mixing
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
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 )
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)
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,jj,jk) = zzwi + zDt_2 * MIN( zWui, 0._wp )
zws(ji,jj,jk) = zzws - zDt_2 * MAX( zWus, 0._wp )
zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zDt_2 * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
END_3D
END SELECT
DO_2D( 0, 0, 0, 0 ) !* Surface boundary conditions
zwi(ji,jj,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)
zWus = ( wi(ji ,jj,2) + wi(ji+1,jj,2) ) / e3u(ji,jj,1,Kaa)
zws(ji,jj,1 ) = zzws - zDt_2 * MAX( zWus, 0._wp )
zwd(ji,jj,1 ) = 1._wp - zzws - zDt_2 * ( MIN( zWus, 0._wp ) )
END_2D
ELSE
SELECT CASE( nldf_dyn )
CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu)
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
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)
zwi(ji,jj,jk) = zzwi
zws(ji,jj,jk) = zzws
zwd(ji,jj,jk) = 1._wp - zzwi - zzws
END_3D
CASE DEFAULT ! iso-level lateral mixing
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
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 )
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)
zwi(ji,jj,jk) = zzwi
zws(ji,jj,jk) = zzws
zwd(ji,jj,jk) = 1._wp - zzwi - zzws
END_3D
END SELECT
DO_2D( 0, 0, 0, 0 ) !* Surface boundary conditions
zwi(ji,jj,1) = 0._wp
zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
END_2D
ENDIF
!
!
! !== Apply semi-implicit bottom friction ==!
!
! Only needed for semi-implicit bottom friction setup. The explicit
! bottom friction has been included in "u(v)a" which act as the R.H.S
! column vector of the tri-diagonal matrix equation
!
IF ( ln_drgimp ) THEN ! implicit bottom friction
DO_2D( 0, 0, 0, 0 )
iku = mbku(ji,jj) ! ocean bottom level at u- and v-points
zwd(ji,jj,iku) = zwd(ji,jj,iku) - zDt_2 *( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) &
& / e3u(ji,jj,iku,Kaa)
END_2D
IF ( ln_isfcav.OR.ln_drgice_imp ) THEN ! top friction (always implicit)
DO_2D( 0, 0, 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,jj,iku) = zwd(ji,jj,iku) - zDt_2 *( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) &
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) &
& / 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) &
& / e3v(ji,jj,ikv,Kaa)
END_1D
END IF
ENDIF
!
! !== Vertical diffusion on u ==!
!
!
! !* Matrix construction
IF( ln_zad_Aimp ) THEN !- including terms associated with partly implicit vertical advection
SELECT CASE( nldf_dyn )
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)
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 )
zws(ji,jk) = zzws - zDt_2 * MAX( zWus, 0._wp )
zwd(ji,jk) = 1._wp - zzwi - zzws + zDt_2 * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
END_2D
CASE DEFAULT ! iso-level lateral mixing
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 )
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)
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 )
zws(ji,jk) = zzws - zDt_2 * MAX( zWus, 0._wp )
zwd(ji,jk) = 1._wp - zzwi - zzws + zDt_2 * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
END_2D
END SELECT
!
zwi(:,1) = 0._wp
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)
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 ) )
END_1D
ELSE !- only vertical diffusive terms
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)
zwi(ji,jk) = zzwi
zws(ji,jk) = zzws
zwd(ji,jk) = 1._wp - zzwi - zzws
END_2D
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 )
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)
zwi(ji,jk) = zzwi
zws(ji,jk) = zzws
zwd(ji,jk) = 1._wp - zzwi - zzws
END_2D
END SELECT
!
zwi(:,1) = 0._wp
DO_1Di( 0, 0 ) !* Surface boundary conditions
zwd(ji,1) = 1._wp - zws(ji,1)
END_1D
ENDIF
!
!
! !== Apply semi-implicit bottom friction ==!
!
! Only needed for semi-implicit bottom friction setup. The explicit
! bottom friction has been included in "u(v)a" which act as the R.H.S
! column vector of the tri-diagonal matrix equation
!
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)
END_2D
END IF
ENDIF
!
! Matrix inversion starting 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 a lower triangular matrix
! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
! The solution (the after velocity) is in puu(:,:,:,Kaa)
!-----------------------------------------------------------------------
!
DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
END_3D
!
DO_2D( 0, 0, 0, 0 ) !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
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)
END_1D
ENDIF
ENDIF
!
! Matrix inversion starting 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 a lower triangular matrix
! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
! The solution (the after velocity) is in puu(:,:,:,Kaa)
!-----------------------------------------------------------------------
!
DO_2Dik( 0, 0, 2, jpkm1, 1 ) !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
zwd(ji,jk) = zwd(ji,jk) - zwi(ji,jk) * zws(ji,jk-1) / zwd(ji,jk-1)
END_2D
!
DO_1Di( 0, 0 ) !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
#if defined key_RK3
! ! RK3: use only utau (not utau_b)
puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * utau(ji,jj) &
& / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1)
! ! 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)
#else
puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + zDt_2 * ( utau_b(ji,jj) + utau(ji,jj) ) &
& / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1)
! ! 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)
#endif
END_2D
DO_3D( 0, 0, 0, 0, 2, jpkm1 )
puu(ji,jj,jk,Kaa) = puu(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * puu(ji,jj,jk-1,Kaa)
END_3D
!
DO_2D( 0, 0, 0, 0 ) !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==!
puu(ji,jj,jpkm1,Kaa) = puu(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1)
END_2D
DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 )
puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - zws(ji,jj,jk) * puu(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk)
END_3D
!
! !== Vertical diffusion on v ==!
!
! !* Matrix construction
IF( ln_zad_Aimp ) THEN !!
SELECT CASE( nldf_dyn )
CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzv)
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
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)
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,jj,jk) = zzwi + zDt_2 * MIN( zWvi, 0._wp )
zws(ji,jj,jk) = zzws - zDt_2 * MAX( zWvs, 0._wp )
zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zDt_2 * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
END_3D
CASE DEFAULT ! iso-level lateral mixing
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
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 )
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)
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,jj,jk) = zzwi + zDt_2 * MIN( zWvi, 0._wp )
zws(ji,jj,jk) = zzws - zDt_2 * MAX( zWvs, 0._wp )
zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zDt_2 * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
END_3D
END SELECT
DO_2D( 0, 0, 0, 0 ) !* Surface boundary conditions
zwi(ji,jj,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)
zWvs = ( wi(ji,jj ,2) + wi(ji,jj+1,2) ) / e3v(ji,jj,1,Kaa)
zws(ji,jj,1 ) = zzws - zDt_2 * MAX( zWvs, 0._wp )
zwd(ji,jj,1 ) = 1._wp - zzws - zDt_2 * ( MIN( zWvs, 0._wp ) )
END_2D
ELSE
SELECT CASE( nldf_dyn )
CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu)
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
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)
zwi(ji,jj,jk) = zzwi
zws(ji,jj,jk) = zzws
zwd(ji,jj,jk) = 1._wp - zzwi - zzws
END_3D
CASE DEFAULT ! iso-level lateral mixing
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
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 )
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)
zwi(ji,jj,jk) = zzwi
zws(ji,jj,jk) = zzws
zwd(ji,jj,jk) = 1._wp - zzwi - zzws
END_3D
END SELECT
DO_2D( 0, 0, 0, 0 ) !* Surface boundary conditions
zwi(ji,jj,1) = 0._wp
zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
END_1D
DO_2Dik( 0, 0, 2, jpkm1, 1 )
puu(ji,jj,jk,Kaa) = puu(ji,jj,jk,Kaa) - zwi(ji,jk) / zwd(ji,jk-1) * puu(ji,jj,jk-1,Kaa)
END_2D
ENDIF
!
! !== Apply semi-implicit top/bottom friction ==!
!
! Only needed for semi-implicit bottom friction setup. The explicit
! bottom friction has been included in "u(v)a" which act as the R.H.S
! column vector of the tri-diagonal matrix equation
!
IF( ln_drgimp ) THEN
DO_2D( 0, 0, 0, 0 )
ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - zDt_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) &
& / e3v(ji,jj,ikv,Kaa)
!
DO_1Di( 0, 0 ) !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==!
puu(ji,jj,jpkm1,Kaa) = puu(ji,jj,jpkm1,Kaa) / zwd(ji,jpkm1)
END_1D
DO_2Dik( 0, 0, jpk-2, 1, -1 )
puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - zws(ji,jk) * puu(ji,jj,jk+1,Kaa) ) / zwd(ji,jk)
END_2D
IF ( ln_isfcav.OR.ln_drgice_imp ) THEN
DO_2D( 0, 0, 0, 0 )
ikv = mikv(ji,jj) ! (first wet ocean u- and v-points)
zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - zDt_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) &
!
!
! !== Vertical diffusion on v ==!
!
! !* Matrix construction
IF( ln_zad_Aimp ) THEN !!
SELECT CASE( nldf_dyn )
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)
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 )
zws(ji,jk) = zzws - zDt_2 * MAX( zWvs, 0._wp )
zwd(ji,jk) = 1._wp - zzwi - zzws - zDt_2 * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
END_2D
CASE DEFAULT ! iso-level lateral mixing
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 )
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)
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 )
zws(ji,jk) = zzws - zDt_2 * MAX( zWvs, 0._wp )
zwd(ji,jk) = 1._wp - zzwi - zzws - zDt_2 * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
END_2D
END SELECT
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)
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 ) )
END_1D
ELSE
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)
zwi(ji,jk) = zzwi
zws(ji,jk) = zzws
zwd(ji,jk) = 1._wp - zzwi - zzws
END_2D
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 )
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)
zwi(ji,jk) = zzwi
zws(ji,jk) = zzws
zwd(ji,jk) = 1._wp - zzwi - zzws
END_2D
END SELECT
DO_1Di( 0, 0 ) !* Surface boundary conditions
zwi(ji,1) = 0._wp
zwd(ji,1) = 1._wp - zws(ji,1)
END_1D
ENDIF
!
! !== Apply semi-implicit top/bottom friction ==!
!
! Only needed for semi-implicit bottom friction setup. The explicit
! bottom friction has been included in "u(v)a" which act as the R.H.S
! column vector of the tri-diagonal matrix equation
!
IF( ln_drgimp ) THEN
DO_1Di( 0, 0 )
ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
zwd(ji,ikv) = zwd(ji,ikv) - zDt_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) &
& / e3v(ji,jj,ikv,Kaa)
END_2D
END_1D
IF ( ln_isfcav.OR.ln_drgice_imp ) THEN
DO_1Di( 0, 0 )
ikv = mikv(ji,jj) ! (first wet ocean u- and v-points)
zwd(ji,ikv) = zwd(ji,ikv) - zDt_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) &
& / e3v(ji,jj,ikv,Kaa)
END_1D
ENDIF
ENDIF
ENDIF
! Matrix inversion
!-----------------------------------------------------------------------
! 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 2d arrays: zwd, zws, zwi
! The solution (after velocity) is in 2d array va
!-----------------------------------------------------------------------
!
DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
END_3D
!
DO_2D( 0, 0, 0, 0 ) !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
! Matrix inversion
!-----------------------------------------------------------------------
! 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 2d arrays: zwd, zws, zwi
! The solution (after velocity) is in 2d array va
!-----------------------------------------------------------------------
!
DO_2Dik( 0, 0, 2, jpkm1, 1 ) !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
zwd(ji,jk) = zwd(ji,jk) - zwi(ji,jk) * zws(ji,jk-1) / zwd(ji,jk-1)
END_2D
!
DO_1Di( 0, 0 ) !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
#if defined key_RK3
! ! RK3: use only vtau (not vtau_b)
pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * vtau(ji,jj) &
& / ( e3v(ji,jj,1,Kaa) * rho0 ) * vmask(ji,jj,1)
! ! RK3: use only vtau (not vtau_b)
pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * vtauV(ji,jj) &
& / ( e3v(ji,jj,1,Kaa) * rho0 ) * vmask(ji,jj,1)
#else
pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + zDt_2*( vtau_b(ji,jj) + vtau(ji,jj) ) &
& / ( e3v(ji,jj,1,Kaa) * rho0 ) * vmask(ji,jj,1)
! ! MLF: average of vtau and vtau_b
pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + zDt_2*( vtau_b(ji,jj) + vtauV(ji,jj) ) &
& / ( e3v(ji,jj,1,Kaa) * rho0 ) * vmask(ji,jj,1)
#endif
END_2D
DO_3D( 0, 0, 0, 0, 2, jpkm1 )
pvv(ji,jj,jk,Kaa) = pvv(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * pvv(ji,jj,jk-1,Kaa)
END_3D
!
DO_2D( 0, 0, 0, 0 ) !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==!
pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1)
END_2D
DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 )
pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - zws(ji,jj,jk) * pvv(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk)
END_3D
END_1D
DO_2Dik( 0, 0, 2, jpkm1, 1 )
pvv(ji,jj,jk,Kaa) = pvv(ji,jj,jk,Kaa) - zwi(ji,jk) / zwd(ji,jk-1) * pvv(ji,jj,jk-1,Kaa)
END_2D
!
DO_1Di( 0, 0 ) !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==!
pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jpkm1)
END_1D
DO_2Dik( 0, 0, jpk-2, 1, -1 )
pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - zws(ji,jk) * pvv(ji,jj,jk+1,Kaa) ) / zwd(ji,jk)
END_2D
! ! ================= !
END_1D ! i-k slices loop !
! ! ================= !
!
IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics
ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) )*r1_Dt - ztrdu(:,:,:)
......
src/OCE/DYN/sshwzv.F90
View file @ 3efbcb17
......@@ -175,7 +175,8 @@ CONTAINS
IF(lwp) WRITE(numout,*) 'wzv_MLF : now vertical velocity '
IF(lwp) WRITE(numout,*) '~~~~~~~'
!
pww(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)
pww(:,:,:) = 0._wp ! bottom boundary condition: w=0 (set once for all)
! ! needed over the halos for the output (ww+wi) in diawri.F90
ENDIF
! !------------------------------!
! ! Now Vertical Velocity !
......@@ -244,28 +245,28 @@ CONTAINS
! inside computational domain (cosmetic)
DO jk = 1, jpkm1
IF( lk_west ) THEN ! --- West --- !
DO ji = mi0(2+nn_hls), mi1(2+nn_hls)
DO ji = mi0(2+nn_hls,nn_hls), mi1(2+nn_hls,nn_hls)
DO jj = 1, jpj
pww(ji,jj,jk) = 0._wp
END DO
END DO
ENDIF
IF( lk_east ) THEN ! --- East --- !
DO ji = mi0(jpiglo-1-nn_hls), mi1(jpiglo-1-nn_hls)
DO ji = mi0(jpiglo-1-nn_hls,nn_hls), mi1(jpiglo-1-nn_hls,nn_hls)
DO jj = 1, jpj
pww(ji,jj,jk) = 0._wp
END DO
END DO
ENDIF
IF( lk_south ) THEN ! --- South --- !
DO jj = mj0(2+nn_hls), mj1(2+nn_hls)
DO jj = mj0(2+nn_hls,nn_hls), mj1(2+nn_hls,nn_hls)
DO ji = 1, jpi
pww(ji,jj,jk) = 0._wp
END DO
END DO
ENDIF
IF( lk_north ) THEN ! --- North --- !
DO jj = mj0(jpjglo-1-nn_hls), mj1(jpjglo-1-nn_hls)
DO jj = mj0(jpjglo-1-nn_hls,nn_hls), mj1(jpjglo-1-nn_hls,nn_hls)
DO ji = 1, jpi
pww(ji,jj,jk) = 0._wp
END DO
......@@ -375,28 +376,28 @@ CONTAINS
! inside computational domain (cosmetic)
DO jk = 1, jpkm1
IF( lk_west ) THEN ! --- West --- !
DO ji = mi0(2+nn_hls), mi1(2+nn_hls)
DO ji = mi0(2+nn_hls,nn_hls), mi1(2+nn_hls,nn_hls)
DO jj = 1, jpj
pww(ji,jj,jk) = 0._wp
END DO
END DO
ENDIF
IF( lk_east ) THEN ! --- East --- !
DO ji = mi0(jpiglo-1-nn_hls), mi1(jpiglo-1-nn_hls)
DO ji = mi0(jpiglo-1-nn_hls,nn_hls), mi1(jpiglo-1-nn_hls,nn_hls)
DO jj = 1, jpj
pww(ji,jj,jk) = 0._wp
END DO
END DO
ENDIF
IF( lk_south ) THEN ! --- South --- !
DO jj = mj0(2+nn_hls), mj1(2+nn_hls)
DO jj = mj0(2+nn_hls,nn_hls), mj1(2+nn_hls,nn_hls)
DO ji = 1, jpi
pww(ji,jj,jk) = 0._wp
END DO
END DO
ENDIF
IF( lk_north ) THEN ! --- North --- !
DO jj = mj0(jpjglo-1-nn_hls), mj1(jpjglo-1-nn_hls)
DO jj = mj0(jpjglo-1-nn_hls,nn_hls), mj1(jpjglo-1-nn_hls,nn_hls)
DO ji = 1, jpi
pww(ji,jj,jk) = 0._wp
END DO
......