MODULE dynldf_lap_blp !!====================================================================== !! *** MODULE dynldf_lap_blp *** !! Ocean dynamics: lateral viscosity trend (laplacian and bilaplacian) !!====================================================================== !! History : 3.7 ! 2014-01 (G. Madec, S. Masson) Original code, re-entrant laplacian !! 4.0 ! 2020-02 (C. Wilson, ...) add bhm coefficient for bi-Laplacian GM implementation via momentum !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! dyn_ldf_lap : update the momentum trend with the lateral viscosity using an iso-level laplacian operator !! dyn_ldf_blp : update the momentum trend with the lateral viscosity using an iso-level bilaplacian operator !!---------------------------------------------------------------------- USE oce ! ocean dynamics and tracers USE dom_oce ! ocean space and time domain USE ldfdyn ! lateral diffusion: eddy viscosity coef. USE ldfslp ! iso-neutral slopes USE zdf_oce ! ocean vertical physics USE phycst ! USE in_out_manager ! I/O manager USE lbclnk ! ocean lateral boundary conditions (or mpp link) IMPLICIT NONE PRIVATE PUBLIC dyn_ldf_lap ! called by dynldf.F90 PUBLIC dyn_ldf_blp ! called by dynldf.F90 PUBLIC dyn_ldf_bgm ! called by dynldf.F90 PRIVATE dyn_ldf_lap_no_ahm !needed by dyn_ldf_bgm !! * Substitutions # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/OCE 4.0 , NEMO Consortium (2018) !! $Id: dynldf_lap_blp.F90 10425 2018-12-19 21:54:16Z smasson $ !! Software governed by the CeCILL license (see ./LICENSE) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE dyn_ldf_lap( kt, pub, pvb, pua, pva, kpass ) !!---------------------------------------------------------------------- !! *** ROUTINE dyn_ldf_lap *** !! !! ** Purpose : Compute the before horizontal momentum diffusive !! trend and add it to the general trend of momentum equation. !! !! ** Method : The Laplacian operator apply on horizontal velocity is !! writen as : grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) ) !! !! ** Action : - pua, pva increased by the harmonic operator applied on pub, pvb. !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: kt ! ocean time-step index INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pub, pvb ! before velocity [m/s] REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pua, pva ! velocity trend [m/s2] ! INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: zsign ! local scalars REAL(wp) :: zua, zva ! local scalars REAL(wp), DIMENSION(jpi,jpj) :: zcur, zdiv !!---------------------------------------------------------------------- ! IF( kt == nit000 .AND. lwp ) THEN WRITE(numout,*) WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator, pass=', kpass WRITE(numout,*) '~~~~~~~ ' ENDIF ! IF( kpass == 1 ) THEN ; zsign = 1._wp ! bilaplacian operator require a minus sign ELSE ; zsign = -1._wp ! (eddy viscosity coef. >0) ENDIF ! ! ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! ! =============== DO jj = 2, jpj DO ji = fs_2, jpi ! vector opt. ! ! ahm * e3 * curl (computed from 1 to jpim1/jpjm1) !!gm open question here : e3f at before or now ? probably now... !!gm note that ahmf has already been multiplied by fmask zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f_n(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1) & & * ( e2v(ji ,jj-1) * pvb(ji ,jj-1,jk) - e2v(ji-1,jj-1) * pvb(ji-1,jj-1,jk) & & - e1u(ji-1,jj ) * pub(ji-1,jj ,jk) + e1u(ji-1,jj-1) * pub(ji-1,jj-1,jk) ) ! ! ahm * div (computed from 2 to jpi/jpj) !!gm note that ahmt has already been multiplied by tmask zdiv(ji,jj) = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk) & & * ( e2u(ji,jj)*e3u_b(ji,jj,jk) * pub(ji,jj,jk) - e2u(ji-1,jj)*e3u_b(ji-1,jj,jk) * pub(ji-1,jj,jk) & & + e1v(ji,jj)*e3v_b(ji,jj,jk) * pvb(ji,jj,jk) - e1v(ji,jj-1)*e3v_b(ji,jj-1,jk) * pvb(ji,jj-1,jk) ) END DO END DO ! DO jj = 2, jpjm1 ! - curl( curl) + grad( div ) DO ji = fs_2, fs_jpim1 ! vector opt. pua(ji,jj,jk) = pua(ji,jj,jk) + zsign * ( & & - ( zcur(ji ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk) & & + ( zdiv(ji+1,jj) - zdiv(ji,jj ) ) * r1_e1u(ji,jj) ) ! pva(ji,jj,jk) = pva(ji,jj,jk) + zsign * ( & & ( zcur(ji,jj ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk) & & + ( zdiv(ji,jj+1) - zdiv(ji ,jj) ) * r1_e2v(ji,jj) ) END DO END DO ! ! =============== END DO ! End of slab ! ! =============== ! END SUBROUTINE dyn_ldf_lap !CW: new subroutine for the Laplacian of velocity, copied from dyn_ldf_lap above, but without eddy viscosity ahmf, ahmt SUBROUTINE dyn_ldf_lap_no_ahm( kt, pub, pvb, pulap, pvlap, kpass ) !!---------------------------------------------------------------------- !! *** ROUTINE dyn_ldf_lap_no_ahm *** !! !! ** Purpose : Companion subroutine to dyn_ldf_bgm to calculate !! the before horizontal momentum Laplacian !! trend and add it to the general trend of momentum equation. !! Note: mimics dyn_ldf_lap but without eddy viscosity ahmf, ahmt, !! because biharmonic GM eddy diffusivity is applied in dyn_bgm. !! !! ** Method : The Laplacian operator apply on horizontal velocity is !! writen as : grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) ) !! !! ** Action : - pua, pva increased by the harmonic operator applied on pub, pvb. !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: kt ! ocean time-step index INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pub, pvb ! before velocity [m/s] REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(out) :: pulap, pvlap ! velocity trend [m/s2] ! INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: zsign ! local scalars REAL(wp) :: zua, zva ! local scalars REAL(wp), DIMENSION(jpi,jpj) :: zcur, zdiv !!---------------------------------------------------------------------- ! IF( kt == nit000 .AND. lwp ) THEN WRITE(numout,*) WRITE(numout,*) 'dyn_ldf_lap_no_ahm : companion iso-level harmonic (laplacian) operator to dyn_bgm, pass=', kpass WRITE(numout,*) '~~~~~~~ ' ENDIF ! IF( kpass == 1 ) THEN ; zsign = 1._wp ! bilaplacian operator require a minus sign ELSE ; zsign = -1._wp ! (eddy viscosity coef. >0) ENDIF ! ! ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! ! =============== DO jj = 2, jpj DO ji = fs_2, jpi ! vector opt. ! ! ahm * e3 * curl (computed from 1 to jpim1/jpjm1) !!gm open question here : e3f at before or now ? probably now... !!CW: note that the removed ahmf had already been multiplied by fmask, so we multiply by fmask explicitly here, !! with the i,j,k of fmask aligning with those of ahmf, as defined in ldfdyn.F90. zcur(ji-1,jj-1) = fmask(ji-1,jj-1,jk) * e3f_n(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1) & & * ( e2v(ji ,jj-1) * pvb(ji ,jj-1,jk) - e2v(ji-1,jj-1) * pvb(ji-1,jj-1,jk) & & - e1u(ji-1,jj ) * pub(ji-1,jj ,jk) + e1u(ji-1,jj-1) * pub(ji-1,jj-1,jk) ) ! ! ahm * div (computed from 2 to jpi/jpj) !!CW: note that the removed ahmt had already been multiplied by tmask, so we multiply by tmask explicitly here, !! with the i,j,k of tmask aligning with those of ahmt, as defined in ldfdyn.F90 zdiv(ji,jj) = tmask(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk) & & * ( e2u(ji,jj)*e3u_b(ji,jj,jk) * pub(ji,jj,jk) - e2u(ji-1,jj)*e3u_b(ji-1,jj,jk) * pub(ji-1,jj,jk) & & + e1v(ji,jj)*e3v_b(ji,jj,jk) * pvb(ji,jj,jk) - e1v(ji,jj-1)*e3v_b(ji,jj-1,jk) * pvb(ji,jj-1,jk) ) END DO END DO ! !CW: multiply pulap, pvlap by umask, vmask DO jj = 2, jpjm1 ! - curl( curl) + grad( div ) DO ji = fs_2, fs_jpim1 ! vector opt. pulap(ji,jj,jk) = umask(ji,jj,jk) * zsign * ( & & - ( zcur(ji ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk) & & + ( zdiv(ji+1,jj) - zdiv(ji,jj ) ) * r1_e1u(ji,jj) ) ! pvlap(ji,jj,jk) = vmask(ji,jj,jk) * zsign * ( & & ( zcur(ji,jj ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk) & & + ( zdiv(ji,jj+1) - zdiv(ji ,jj) ) * r1_e2v(ji,jj) ) END DO END DO ! ! =============== END DO ! End of slab ! ! =============== ! END SUBROUTINE dyn_ldf_lap_no_ahm !CW: new subroutine for biharmonic GM, copied from SUBROUTINE dyn_ldf_blp below SUBROUTINE dyn_ldf_bgm( kt, pub, pvb, pua, pva ) !!---------------------------------------------------------------------- !! *** ROUTINE dyn_bgm *** !! !! ** Purpose : Compute the before lateral momentum trend due to the bi-Laplacian GM parameterisation !! and add it to the general trend of momentum equation. !! !! ** Method : The bi-Laplacian implementation of GM is via a -d/dz(diffusivity x d/dz(Laplacian of velocity)) !! operator applied at the 'now' time level. The existing code for the Laplacian contains the 'before' time also in zdiv. !! It is computed by a call to dyn_ldf_lap routine and vertical differentiation applied twice. !! !! ** Action : pua, pva increased with the before bi-Laplacian GM momentum trend calculated from pub, pvb. !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: kt ! ocean time-step index REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pub, pvb ! before velocity fields REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pua, pva ! momentum trend ! INTEGER :: iku, ikv ! local integers !CW INTEGER :: ji, jj, jk ! dummy loop indices ! REAL(wp), DIMENSION(jpi,jpj,jpk) :: zulap, zvlap ! laplacian at u- and v-point REAL(wp), DIMENSION(jpi,jpj,jpk) :: zulapdz, zvlapdz ! -1*bhm * d/dz(del^2 u) at u- and v-point REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmu ! = bhm / avm !!---------------------------------------------------------------------- ! IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn_bgm : bi-Laplacian GM operator via momentum ' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' ENDIF ! ! Calculate ratio of bhm and avm for stabilising correction terms ! and add to avm to be included in the vertical diffusion calculation later. DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, jpim1 ! vector opt. zmu(ji,jj,jk) = ( bhm(ji,jj,jk) / (avm(ji,jj,jk) + rsmall) ) * wmask(ji,jj,jk) avm(ji,jj,jk) = avm(ji,jj,jk) + zmu(ji,jj,jk) ENDDO ENDDO ENDDO CALL lbc_lnk_multi( 'dyn_ldf_bgm', zmu, 'W', 1. , avm, 'W', 1. ) ! Calculate (del2 u) CALL dyn_ldf_lap_no_ahm( kt, pub, pvb, zulap, zvlap, 1 ) ! =============== !CW: calculate -bhm * d/dz(del^2 u) ! Use of wumask and wvmask to ensure diffusive fluxes at topography are zero. DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, jpim1 ! vector opt. zulapdz(ji,jj,jk) = ( -0.5_wp*(bhm(ji+1,jj ,jk)+bhm(ji,jj,jk))*(zulap(ji,jj,jk-1) - zulap(ji,jj,jk)) & & -0.5_wp*(zmu(ji+1,jj ,jk)+zmu(ji,jj,jk))*(pub (ji,jj,jk-1) - pub (ji,jj,jk)) ) * wumask(ji,jj,jk) / e3uw_n(ji,jj,jk) zvlapdz(ji,jj,jk) = ( -0.5_wp*(bhm(ji ,jj+1,jk)+bhm(ji,jj,jk))*(zvlap(ji,jj,jk-1) - zvlap(ji,jj,jk)) & & -0.5_wp*(zmu(ji ,jj+1,jk)+zmu(ji,jj,jk))*(pvb (ji,jj,jk-1) - pvb (ji,jj,jk)) ) * wumask(ji,jj,jk) / e3uw_n(ji,jj,jk) ENDDO ENDDO ENDDO !CW: set boundary conditions: d/dz(del^2 u) = 0 at top and bottom, so that eddy-induced velocity, w*=0 ! Surface zulapdz(:,:,1) = 0._wp ; zvlapdz(:,:,1) = 0._wp ! Flat bottom case zulapdz(:,:,jpk) = 0._wp ; zvlapdz(:,:,jpk) = 0._wp ! Variable bathymetry case, including z-partial steps, as in dynhpg.F90, subroutine hpg_zps DO jj = 2, jpjm1 DO ji = 2, jpim1 iku = mbku(ji,jj)+1 ikv = mbkv(ji,jj)+1 zulapdz(:,:,iku) = 0._wp zvlapdz(:,:,ikv) = 0._wp ENDDO ENDDO !! calculate d/dz(-bhm * d/dz(del^2 u)) ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! ! =============== DO jj = 2, jpjm1 DO ji = fs_2, jpim1 ! vector opt. pua(ji,jj,jk) = pua(ji,jj,jk) + (zulapdz(ji,jj,jk) - zulapdz(ji,jj,jk+1)) / e3u_n(ji,jj,jk) pva(ji,jj,jk) = pva(ji,jj,jk) + (zvlapdz(ji,jj,jk) - zvlapdz(ji,jj,jk+1)) / e3v_n(ji,jj,jk) ENDDO ENDDO ENDDO ! ----- END SUBROUTINE dyn_ldf_bgm SUBROUTINE dyn_ldf_blp( kt, pub, pvb, pua, pva ) !!---------------------------------------------------------------------- !! *** ROUTINE dyn_ldf_blp *** !! !! ** Purpose : Compute the before lateral momentum viscous trend !! and add it to the general trend of momentum equation. !! !! ** Method : The lateral viscous trends is provided by a bilaplacian !! operator applied to before field (forward in time). !! It is computed by two successive calls to dyn_ldf_lap routine !! !! ** Action : pta updated with the before rotated bilaplacian diffusion !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: kt ! ocean time-step index REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pub, pvb ! before velocity fields REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pua, pva ! momentum trend ! REAL(wp), DIMENSION(jpi,jpj,jpk) :: zulap, zvlap ! laplacian at u- and v-point !!---------------------------------------------------------------------- ! IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn_ldf_blp : bilaplacian operator momentum ' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' ENDIF ! zulap(:,:,:) = 0._wp zvlap(:,:,:) = 0._wp ! CALL dyn_ldf_lap( kt, pub, pvb, zulap, zvlap, 1 ) ! rotated laplacian applied to ptb (output in zlap) ! CALL lbc_lnk_multi( 'dynldf_lap_blp', zulap, 'U', -1., zvlap, 'V', -1. ) ! Lateral boundary conditions ! CALL dyn_ldf_lap( kt, zulap, zvlap, pua, pva, 2 ) ! rotated laplacian applied to zlap (output in pta) ! END SUBROUTINE dyn_ldf_blp !!====================================================================== END MODULE dynldf_lap_blp