MODULE dynvor !!====================================================================== !! *** MODULE dynvor *** !! Ocean dynamics: Update the momentum trend with the relative and !! planetary vorticity trends !!====================================================================== !! History : OPA ! 1989-12 (P. Andrich) vor_ens: Original code !! 5.0 ! 1991-11 (G. Madec) vor_ene, vor_mix: Original code !! 6.0 ! 1996-01 (G. Madec) s-coord, suppress work arrays !! NEMO 0.5 ! 2002-08 (G. Madec) F90: Free form and module !! 1.0 ! 2004-02 (G. Madec) vor_een: Original code !! - ! 2003-08 (G. Madec) add vor_ctl !! - ! 2005-11 (G. Madec) add dyn_vor (new step architecture) !! 2.0 ! 2006-11 (G. Madec) flux form advection: add metric term !! 3.2 ! 2009-04 (R. Benshila) vvl: correction of een scheme !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase !! 3.7 ! 2014-04 (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! dyn_vor : Update the momentum trend with the vorticity trend !! vor_ens : enstrophy conserving scheme (ln_dynvor_ens=T) !! vor_ene : energy conserving scheme (ln_dynvor_ene=T) !! vor_mix : mixed enstrophy/energy conserving (ln_dynvor_mix=T) !! vor_een : energy and enstrophy conserving (ln_dynvor_een=T) !! dyn_vor_init : set and control of the different vorticity option !!---------------------------------------------------------------------- USE oce ! ocean dynamics and tracers USE dom_oce ! ocean space and time domain USE dommsk ! ocean mask USE dynadv ! momentum advection (use ln_dynadv_vec value) USE trd_oce ! trends: ocean variables USE trddyn ! trend manager: dynamics USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE prtctl ! Print control USE in_out_manager ! I/O manager USE lib_mpp ! MPP library USE wrk_nemo ! Memory Allocation USE timing ! Timing IMPLICIT NONE PRIVATE PUBLIC dyn_vor ! routine called by step.F90 PUBLIC dyn_vor_init ! routine called by opa.F90 ! !!* Namelist namdyn_vor: vorticity term LOGICAL, PUBLIC :: ln_dynvor_ene !: energy conserving scheme LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme LOGICAL, PUBLIC :: ln_dynvor_een !: energy and enstrophy conserving scheme LOGICAL, PUBLIC :: ln_dynvor_een_old !: energy and enstrophy conserving scheme (original formulation) INTEGER :: nvor = 0 ! type of vorticity trend used INTEGER :: ncor = 1 ! coriolis INTEGER :: nrvm = 2 ! =2 relative vorticity ; =3 metric term INTEGER :: ntot = 4 ! =4 total vorticity (relative + planetary) ; =5 coriolis + metric term !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/OPA 3.3 , NEMO Consortium (2010) !! $Id: dynvor.F90 4990 2014-12-15 16:42:49Z timgraham $ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE dyn_vor( kt ) !!---------------------------------------------------------------------- !! !! ** Purpose : compute the lateral ocean tracer physics. !! !! ** Action : - Update (ua,va) with the now vorticity term trend !! - save the trends in (ztrdu,ztrdv) in 2 parts (relative !! and planetary vorticity trends) and send them to trd_dyn !! for futher diagnostics (l_trddyn=T) !!---------------------------------------------------------------------- INTEGER, INTENT( in ) :: kt ! ocean time-step index ! REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv !!---------------------------------------------------------------------- ! IF( nn_timing == 1 ) CALL timing_start('dyn_vor') ! IF( l_trddyn ) CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) ! ! ! vorticity term SELECT CASE ( nvor ) ! compute the vorticity trend and add it to the general trend ! CASE ( -1 ) ! esopa: test all possibility with control print CALL vor_ene( kt, ntot, ua, va ) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, & & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) CALL vor_ens( kt, ntot, ua, va ) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, & & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) CALL vor_mix( kt ) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, & & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) CALL vor_een( kt, ntot, ua, va ) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, & & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) ! CASE ( 0 ) ! energy conserving scheme IF( l_trddyn ) THEN ztrdu(:,:,:) = ua(:,:,:) ztrdv(:,:,:) = va(:,:,:) CALL vor_ene( kt, nrvm, ua, va ) ! relative vorticity or metric trend ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) ztrdu(:,:,:) = ua(:,:,:) ztrdv(:,:,:) = va(:,:,:) CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) ELSE CALL vor_ene( kt, ntot, ua, va ) ! total vorticity ENDIF ! CASE ( 1 ) ! enstrophy conserving scheme IF( l_trddyn ) THEN ztrdu(:,:,:) = ua(:,:,:) ztrdv(:,:,:) = va(:,:,:) CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) ztrdu(:,:,:) = ua(:,:,:) ztrdv(:,:,:) = va(:,:,:) CALL vor_ens( kt, ncor, ua, va ) ! planetary vorticity trend ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) ELSE CALL vor_ens( kt, ntot, ua, va ) ! total vorticity ENDIF ! CASE ( 2 ) ! mixed ene-ens scheme IF( l_trddyn ) THEN ztrdu(:,:,:) = ua(:,:,:) ztrdv(:,:,:) = va(:,:,:) CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend (ens) ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) ztrdu(:,:,:) = ua(:,:,:) ztrdv(:,:,:) = va(:,:,:) CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend (ene) ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) ELSE CALL vor_mix( kt ) ! total vorticity (mix=ens-ene) ENDIF ! CASE ( 3 ) ! energy and enstrophy conserving scheme IF( l_trddyn ) THEN ztrdu(:,:,:) = ua(:,:,:) ztrdv(:,:,:) = va(:,:,:) CALL vor_een( kt, nrvm, ua, va ) ! relative vorticity or metric trend ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) ztrdu(:,:,:) = ua(:,:,:) ztrdv(:,:,:) = va(:,:,:) CALL vor_een( kt, ncor, ua, va ) ! planetary vorticity trend ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) ELSE ! CALL vor_een( kt, ntot, ua, va ) ! total vorticity CALL vor_een( kt, nrvm, ua, va ) ! Do not use the above compact form CALL vor_een( kt, ncor, ua, va ) ! in order to preserve reproducibility ENDIF ! END SELECT ! ! ! print sum trends (used for debugging) IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor - Ua: ', mask1=umask, & & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) ! IF( l_trddyn ) CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) ! IF( nn_timing == 1 ) CALL timing_stop('dyn_vor') ! END SUBROUTINE dyn_vor SUBROUTINE vor_ene( kt, kvor, pua, pva ) !!---------------------------------------------------------------------- !! *** ROUTINE vor_ene *** !! !! ** 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 Sadourny (1975) flux form formulation : conserves the !! horizontal kinetic energy. !! The trend of the vorticity term is given by: !! * s-coordinate (ln_sco=T), the e3. are inside the derivatives: !! voru = 1/e1u mj-1[ (rotn+f)/e3f mi(e1v*e3v vn) ] !! vorv = 1/e2v mi-1[ (rotn+f)/e3f mj(e2u*e3u un) ] !! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: !! voru = 1/e1u mj-1[ (rotn+f) mi(e1v vn) ] !! vorv = 1/e2v mi-1[ (rotn+f) mj(e2u un) ] !! Add this trend to the general momentum trend (ua,va): !! (ua,va) = (ua,va) + ( voru , vorv ) !! !! ** Action : - Update (ua,va) with the now vorticity term trend !! !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. !!---------------------------------------------------------------------- ! INTEGER , INTENT(in ) :: kt ! ocean time-step index INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; ! ! =nrvm (relative vorticity or metric) REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend ! INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: zx1, zy1, zfact2, zx2, zy2 ! local scalars REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz !!---------------------------------------------------------------------- ! IF( nn_timing == 1 ) CALL timing_start('vor_ene') ! CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz ) ! IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' ENDIF zfact2 = 0.5 * 0.5 ! Local constant initialization !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz ) ! ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! ! =============== ! ! Potential vorticity and horizontal fluxes ! ----------------------------------------- SELECT CASE( kvor ) ! vorticity considered CASE ( 1 ) ; zwz(:,:) = ff(:,:) ! planetary vorticity (Coriolis) CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) ! relative vorticity CASE ( 3 ) ! metric term DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) END DO END DO CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ! total (relative + planetary vorticity) CASE ( 5 ) ! total (coriolis + metric) DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zwz(ji,jj) = ( ff (ji,jj) & & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) & & ) END DO END DO END SELECT IF( ln_sco ) THEN zwz(:,:) = zwz(:,:) / fse3f(:,:,jk) zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) ELSE zwx(:,:) = e2u(:,:) * un(:,:,jk) zwy(:,:) = e1v(:,:) * vn(:,:,jk) ENDIF ! Compute and add the vorticity term trend ! ---------------------------------------- DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) zy2 = zwy(ji,jj ) + zwy(ji+1,jj ) zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) END DO END DO ! ! =============== END DO ! End of slab ! ! =============== CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz ) ! IF( nn_timing == 1 ) CALL timing_stop('vor_ene') ! END SUBROUTINE vor_ene SUBROUTINE vor_mix( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE vor_mix *** !! !! ** 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) !! Mixte formulation : conserves the potential enstrophy of a hori- !! zontally non-divergent flow for (rotzu x uh), the relative vor- !! ticity term and the horizontal kinetic energy for (f x uh), the !! coriolis term. the now trend of the vorticity term is given by: !! * s-coordinate (ln_sco=T), the e3. are inside the derivatives: !! voru = 1/e1u mj-1(rotn/e3f) mj-1[ mi(e1v*e3v vn) ] !! +1/e1u mj-1[ f/e3f mi(e1v*e3v vn) ] !! vorv = 1/e2v mi-1(rotn/e3f) mi-1[ mj(e2u*e3u un) ] !! +1/e2v mi-1[ f/e3f mj(e2u*e3u un) ] !! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: !! voru = 1/e1u mj-1(rotn) mj-1[ mi(e1v vn) ] !! +1/e1u mj-1[ f mi(e1v vn) ] !! vorv = 1/e2v mi-1(rotn) mi-1[ mj(e2u un) ] !! +1/e2v mi-1[ f mj(e2u un) ] !! Add this now trend to the general momentum trend (ua,va): !! (ua,va) = (ua,va) + ( voru , vorv ) !! !! ** Action : - Update (ua,va) arrays with the now vorticity term trend !! !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. !!---------------------------------------------------------------------- ! INTEGER, INTENT(in) :: kt ! ocean timestep index ! INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: zfact1, zua, zcua, zx1, zy1 ! local scalars REAL(wp) :: zfact2, zva, zcva, zx2, zy2 ! - - REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww !!---------------------------------------------------------------------- ! IF( nn_timing == 1 ) CALL timing_start('vor_mix') ! CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz, zww ) ! IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn:vor_mix : vorticity term: mixed energy/enstrophy conserving scheme' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' ENDIF zfact1 = 0.5 * 0.25 ! Local constant initialization zfact2 = 0.5 * 0.5 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww ) ! ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! ! =============== ! ! Relative and planetary potential vorticity and horizontal fluxes ! ---------------------------------------------------------------- IF( ln_sco ) THEN IF( ln_dynadv_vec ) THEN zww(:,:) = rotn(:,:,jk) / fse3f(:,:,jk) ELSE DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zww(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & & * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) * fse3f(ji,jj,jk) ) END DO END DO ENDIF zwz(:,:) = ff (:,:) / fse3f(:,:,jk) zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) ELSE IF( ln_dynadv_vec ) THEN zww(:,:) = rotn(:,:,jk) ELSE DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zww(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & & * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) ) END DO END DO ENDIF zwz(:,:) = ff (:,:) zwx(:,:) = e2u(:,:) * un(:,:,jk) zwy(:,:) = e1v(:,:) * vn(:,:,jk) ENDIF ! Compute and add the vorticity term trend ! ---------------------------------------- DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj) zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj) zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj) zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj) ! enstrophy conserving formulation for relative vorticity term zua = zfact1 * ( zww(ji ,jj-1) + zww(ji,jj) ) * ( zy1 + zy2 ) zva =-zfact1 * ( zww(ji-1,jj ) + zww(ji,jj) ) * ( zx1 + zx2 ) ! energy conserving formulation for planetary vorticity term zcua = zfact2 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) zcva =-zfact2 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) ! mixed vorticity trend added to the momentum trends ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva END DO END DO ! ! =============== END DO ! End of slab ! ! =============== CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz, zww ) ! IF( nn_timing == 1 ) CALL timing_stop('vor_mix') ! END SUBROUTINE vor_mix SUBROUTINE vor_ens( kt, kvor, pua, pva ) !!---------------------------------------------------------------------- !! *** ROUTINE vor_ens *** !! !! ** 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 Sadourny (1975) flux FORM formulation : conserves the !! potential enstrophy of a horizontally non-divergent flow. the !! trend of the vorticity term is given by: !! * s-coordinate (ln_sco=T), the e3. are inside the derivative: !! voru = 1/e1u mj-1[ (rotn+f)/e3f ] mj-1[ mi(e1v*e3v vn) ] !! vorv = 1/e2v mi-1[ (rotn+f)/e3f ] mi-1[ mj(e2u*e3u un) ] !! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: !! voru = 1/e1u mj-1[ rotn+f ] mj-1[ mi(e1v vn) ] !! vorv = 1/e2v mi-1[ rotn+f ] mi-1[ mj(e2u un) ] !! Add this trend to the general momentum trend (ua,va): !! (ua,va) = (ua,va) + ( voru , vorv ) !! !! ** Action : - Update (ua,va) arrays with the now vorticity term trend !! !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. !!---------------------------------------------------------------------- ! INTEGER , INTENT(in ) :: kt ! ocean time-step index INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; ! ! =nrvm (relative vorticity or metric) REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend ! INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: zfact1, zuav, zvau ! temporary scalars REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww !!---------------------------------------------------------------------- ! IF( nn_timing == 1 ) CALL timing_start('vor_ens') ! CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz ) ! IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' ENDIF zfact1 = 0.5 * 0.25 ! Local constant initialization !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz ) ! ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! ! =============== ! ! Potential vorticity and horizontal fluxes ! ----------------------------------------- SELECT CASE( kvor ) ! vorticity considered CASE ( 1 ) ; zwz(:,:) = ff(:,:) ! planetary vorticity (Coriolis) CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) ! relative vorticity CASE ( 3 ) ! metric term DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) END DO END DO CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ! total (relative + planetary vorticity) CASE ( 5 ) ! total (coriolis + metric) DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zwz(ji,jj) = ( ff (ji,jj) & & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) & & ) END DO END DO END SELECT ! IF( ln_sco ) THEN DO jj = 1, jpj ! caution: don't use (:,:) for this loop DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk) zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk) zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk) END DO END DO ELSE DO jj = 1, jpj ! caution: don't use (:,:) for this loop DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk) zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk) END DO END DO ENDIF ! ! Compute and add the vorticity term trend ! ---------------------------------------- DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zuav = zfact1 / e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) END DO END DO ! ! =============== END DO ! End of slab ! ! =============== CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz ) ! IF( nn_timing == 1 ) CALL timing_stop('vor_ens') ! END SUBROUTINE vor_ens SUBROUTINE vor_een( kt, kvor, pua, pva ) !!---------------------------------------------------------------------- !! *** 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 (ua,va). !! !! ** Action : - Update (ua,va) 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 ) :: kvor ! =ncor (planetary) ; =ntot (total) ; ! ! =nrvm (relative vorticity or metric) REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend !! INTEGER :: ji, jj, jk ! dummy loop indices INTEGER :: ierr ! local integer REAL(wp) :: zfac12, zua, zva ! local scalars REAL(wp) :: zmsk, ze3 ! local scalars ! ! 3D workspace REAL(wp), POINTER , DIMENSION(:,: ) :: zwx, zwy, zwz REAL(wp), POINTER , DIMENSION(:,: ) :: ztnw, ztne, ztsw, ztse #if defined key_vvl REAL(wp), POINTER , DIMENSION(:,:,:) :: ze3f ! 3D workspace (lk_vvl=T) #else REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: ze3f ! lk_vvl=F, ze3f=1/e3f saved one for all #endif !!---------------------------------------------------------------------- ! IF( nn_timing == 1 ) CALL timing_start('vor_een') ! CALL wrk_alloc( jpi, jpj, zwx , zwy , zwz ) CALL wrk_alloc( jpi, jpj, ztnw, ztne, ztsw, ztse ) #if defined key_vvl CALL wrk_alloc( jpi, jpj, jpk, ze3f ) #endif ! 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,*) '~~~~~~~~~~~' #if ! defined key_vvl IF( .NOT.ALLOCATED(ze3f) ) THEN ALLOCATE( ze3f(jpi,jpj,jpk) , STAT=ierr ) IF( lk_mpp ) CALL mpp_sum ( ierr ) IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' ) ENDIF ze3f(:,:,:) = 0._wp #endif ENDIF IF( kt == nit000 .OR. lk_vvl ) THEN ! reciprocal of e3 at F-point (masked averaging of e3t over ocean points) IF( ln_dynvor_een_old ) THEN ! original formulation DO jk = 1, jpk DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ze3 = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) IF ( ze3 /= 0._wp ) THEN ; ze3f(ji,jj,jk) = 4.0_wp / ze3 ELSE ; ze3f(ji,jj,jk) = 0.0_wp ENDIF END DO END DO END DO ELSE ! new formulation from NEMO 3.6 DO jk = 1, jpk DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ze3 = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*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 ( ze3 /= 0._wp ) THEN ; ze3f(ji,jj,jk) = zmsk / ze3 ELSE ; ze3f(ji,jj,jk) = 0.0_wp ENDIF END DO END DO END DO ENDIF CALL lbc_lnk( ze3f, 'F', 1. ) ENDIF zfac12 = 1._wp / 12._wp ! Local constant initialization !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse ) ! ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! ! =============== ! Potential vorticity and horizontal fluxes ! ----------------------------------------- SELECT CASE( kvor ) ! vorticity considered CASE ( 1 ) ! planetary vorticity (Coriolis) zwz(:,:) = ff(:,:) * ze3f(:,:,jk) CASE ( 2 ) ! relative vorticity zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk) CASE ( 3 ) ! metric term DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk) END DO END DO CALL lbc_lnk( zwz, 'F', 1. ) CASE ( 4 ) ! total (relative + planetary vorticity) zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) CASE ( 5 ) ! total (coriolis + metric) DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zwz(ji,jj) = ( ff (ji,jj) & & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) & & ) * ze3f(ji,jj,jk) END DO END DO CALL lbc_lnk( zwz, 'F', 1. ) END SELECT zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) ! Compute and add the vorticity term trend ! ---------------------------------------- jj = 2 ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ; ztsw(1,:) = 0 DO ji = 2, jpi 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 DO DO jj = 3, jpj DO ji = fs_2, jpi ! vector opt. ok because we start at jj = 3 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 DO END DO DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zua = + zfac12 / 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 = - zfac12 / 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 ) ) pua(ji,jj,jk) = pua(ji,jj,jk) + zua pva(ji,jj,jk) = pva(ji,jj,jk) + zva END DO END DO ! ! =============== END DO ! End of slab ! ! =============== CALL wrk_dealloc( jpi, jpj, zwx , zwy , zwz ) CALL wrk_dealloc( jpi, jpj, ztnw, ztne, ztsw, ztse ) #if defined key_vvl CALL wrk_dealloc( jpi, jpj, jpk, ze3f ) #endif ! IF( nn_timing == 1 ) CALL timing_stop('vor_een') ! END SUBROUTINE vor_een SUBROUTINE dyn_vor_init !!--------------------------------------------------------------------- !! *** ROUTINE dyn_vor_init *** !! !! ** Purpose : Control the consistency between cpp options for !! tracer advection schemes !!---------------------------------------------------------------------- INTEGER :: ioptio ! local integer INTEGER :: ji, jj, jk ! dummy loop indices INTEGER :: ios ! Local integer output status for namelist read !! NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, ln_dynvor_een_old !!---------------------------------------------------------------------- REWIND( numnam_ref ) ! Namelist namdyn_vor in reference namelist : Vorticity scheme options READ ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901) 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in reference namelist', lwp ) REWIND( numnam_cfg ) ! Namelist namdyn_vor in configuration namelist : Vorticity scheme options READ ( numnam_cfg, namdyn_vor, IOSTAT = ios, ERR = 902 ) 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist', lwp ) IF(lwm) WRITE ( numond, namdyn_vor ) IF(lwp) THEN ! Namelist print WRITE(numout,*) WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency' WRITE(numout,*) '~~~~~~~~~~~~' WRITE(numout,*) ' Namelist namdyn_vor : choice of the vorticity term scheme' WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een WRITE(numout,*) ' enstrophy and energy conserving scheme (old) ln_dynvor_een_old= ', ln_dynvor_een_old ENDIF ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks ! at angles with three ocean points and one land point IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN DO jk = 1, jpk DO jj = 2, jpjm1 DO ji = 2, jpim1 IF( tmask(ji,jj,jk)+tmask(ji+1,jj,jk)+tmask(ji,jj+1,jk)+tmask(ji+1,jj+1,jk) == 3._wp ) & fmask(ji,jj,jk) = 1._wp END DO END DO END DO ! CALL lbc_lnk( fmask, 'F', 1._wp ) ! Lateral boundary conditions on fmask ! ENDIF ioptio = 0 ! Control of vorticity scheme options IF( ln_dynvor_ene ) ioptio = ioptio + 1 IF( ln_dynvor_ens ) ioptio = ioptio + 1 IF( ln_dynvor_mix ) ioptio = ioptio + 1 IF( ln_dynvor_een ) ioptio = ioptio + 1 IF( ln_dynvor_een_old ) ioptio = ioptio + 1 IF( lk_esopa ) ioptio = 1 IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' ) ! ! Set nvor (type of scheme for vorticity) IF( ln_dynvor_ene ) nvor = 0 IF( ln_dynvor_ens ) nvor = 1 IF( ln_dynvor_mix ) nvor = 2 IF( ln_dynvor_een .or. ln_dynvor_een_old ) nvor = 3 IF( lk_esopa ) nvor = -1 ! ! Set ncor, nrvm, ntot (type of vorticity) IF(lwp) WRITE(numout,*) ncor = 1 IF( ln_dynadv_vec ) THEN IF(lwp) WRITE(numout,*) ' Vector form advection : vorticity = Coriolis + relative vorticity' nrvm = 2 ntot = 4 ELSE IF(lwp) WRITE(numout,*) ' Flux form advection : vorticity = Coriolis + metric term' nrvm = 3 ntot = 5 ENDIF IF(lwp) THEN ! Print the choice WRITE(numout,*) IF( nvor == 0 ) WRITE(numout,*) ' vorticity scheme : energy conserving scheme' IF( nvor == 1 ) WRITE(numout,*) ' vorticity scheme : enstrophy conserving scheme' IF( nvor == 2 ) WRITE(numout,*) ' vorticity scheme : mixed enstrophy/energy conserving scheme' IF( nvor == 3 ) WRITE(numout,*) ' vorticity scheme : energy and enstrophy conserving scheme' IF( nvor == -1 ) WRITE(numout,*) ' esopa test: use all lateral physics options' ENDIF ! END SUBROUTINE dyn_vor_init !!============================================================================== END MODULE dynvor