# Changeset 10094 for NEMO/branches/UKMO/dev_r10037_dynvor_EEUV

Ignore:
Timestamp:
2018-09-06T15:24:19+02:00 (3 years ago)
Message:

UKMO dev_r10037_dynvor_EEUV branch : add new dynvor scheme.

Location:
NEMO/branches/UKMO/dev_r10037_dynvor_EEUV
Files:
3 edited

Unmodified
Removed
• ## NEMO/branches/UKMO/dev_r10037_dynvor_EEUV/cfgs/SHARED/namelist_ref

 r9950 ln_dynvor_enT = .false. !  energy conserving scheme (T-point) ln_dynvor_eeT = .false. !  energy conserving scheme (een using e3t) ln_dynvor_eeUV = .false.!  energy conserving scheme (een using e3u and e3v) ln_dynvor_een = .false. !  energy & enstrophy scheme nn_een_e3f = 1          ! =0  e3f = mi(mj(e3t))/4

• ## NEMO/branches/UKMO/dev_r10037_dynvor_EEUV/src/OCE/DYN/dynvor.F90

 r9950 LOGICAL, PUBLIC ::   ln_dynvor_enT   !: t-point energy conserving scheme     (ENT) LOGICAL, PUBLIC ::   ln_dynvor_eeT   !: t-point energy conserving scheme     (EET) LOGICAL, PUBLIC ::   ln_dynvor_eeUV  !: uv-point energy conserving scheme    (EEUV) LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy & enstrophy conserving scheme (EEN) INTEGER, PUBLIC ::      nn_een_e3f      !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) INTEGER, PUBLIC, PARAMETER ::   np_ENT = 2   ! ENT scheme (t-point vorticity) INTEGER, PUBLIC, PARAMETER ::   np_EET = 3   ! EET scheme (EEN using e3t) INTEGER, PUBLIC, PARAMETER ::   np_EEN = 4   ! EEN scheme INTEGER, PUBLIC, PARAMETER ::   np_MIX = 5   ! MIX scheme INTEGER, PUBLIC, PARAMETER ::   np_EEUV = 4  ! EET scheme (EEN using e3u and e3v) INTEGER, PUBLIC, PARAMETER ::   np_EEN = 5   ! EEN scheme INTEGER, PUBLIC, PARAMETER ::   np_MIX = 6   ! MIX scheme INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity REAL(wp) ::   r1_4  = 0.250_wp         ! =1/4 REAL(wp) ::   r1_6  = 1._wp / 6._wp    ! =1/6 REAL(wp) ::   r1_8  = 0.125_wp         ! =1/8 REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! 1/12 CASE( np_EET )           ;   CALL vor_eeT( kt, ncor, un , vn , ua, va )   ! energy conserving scheme (een with e3t) IF( ln_stcor )            CALL vor_eeT( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend CASE( np_EEUV )          ;   CALL vor_eeUV( kt, ncor, un , vn , ua, va )  ! energy conserving scheme (een with e3u and e3v) IF( ln_stcor )            CALL vor_eeUV( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend CASE( np_EEN )           ;   CALL vor_een( kt, ncor, un , vn , ua, va )   ! energy & enstrophy scheme IF( ln_stcor )            CALL vor_een( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend CASE( np_ENT )           ;   CALL vor_enT( kt, nrvm, un , vn , ua, va )  ! energy conserving scheme (T-pts) CASE( np_EET )           ;   CALL vor_eeT( kt, nrvm, un , vn , ua, va )  ! energy conserving scheme (een with e3t) CASE( np_EEUV )          ;   CALL vor_eeUV( kt, nrvm, un , vn , ua, va ) ! energy conserving scheme (een with e3u and e3v) CASE( np_ENE )           ;   CALL vor_ene( kt, nrvm, un , vn , ua, va )  ! energy conserving scheme CASE( np_ENS, np_MIX )   ;   CALL vor_ens( kt, nrvm, un , vn , ua, va )  ! enstrophy conserving scheme CALL vor_eeT( kt, ntot, un , vn , ua, va )   ! total vorticity trend IF( ln_stcor )   CALL vor_eeT( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend CASE( np_EEUV )                       !* energy conserving scheme (een scheme using e3u and e3v) CALL vor_eeUV( kt, ntot, un , vn , ua, va )  ! total vorticity trend IF( ln_stcor )   CALL vor_eeUV( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend CASE( np_ENE )                        !* energy conserving scheme CALL vor_ene( kt, ntot, un , vn , ua, va )   ! total vorticity trend SUBROUTINE vor_eeUV( kt, kvor, pun, pvn, pua, pva ) !!---------------------------------------------------------------------- !!                ***  ROUTINE vor_eeUV  *** !! !! ** 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). !! !!      Modified version of EEN as suggested by Mike Bell. !! !! ** Action : - Update (ua,va) with the now vorticity term trend !! !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 !!              Bell, Johnson and Marshall, unpublished notes 2017 !!---------------------------------------------------------------------- INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! 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(jpi,jpj) ::   zwx , zwy , zwz REAL(wp), DIMENSION(jpi,jpj) ::   ztnw, ztne, ztsw, ztse !!---------------------------------------------------------------------- ! 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 ! !                                                ! =============== DO jk = 1, jpkm1                                 ! Horizontal slab !                                             ! =============== ! SELECT CASE( kvor )                 !==  vorticity considered  ==! CASE ( np_COR )                           !* Coriolis (planetary vorticity) DO jj = 1, jpjm1 DO ji = 1, fs_jpim1   ! vector opt. zwz(ji,jj) = ff_f(ji,jj) END DO END DO CASE ( np_RVO )                           !* relative vorticity DO jj = 1, jpjm1 DO ji = 1, fs_jpim1   ! vector opt. zwz(ji,jj) = ( e2v(ji+1,jj  ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)  & &         - e1u(ji  ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj) END DO END DO CASE ( np_MET )                           !* metric term DO jj = 1, jpjm1 DO ji = 1, fs_jpim1   ! vector opt. zwz(ji,jj) = (   ( pvn(ji+1,jj,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & &           - ( pun(ji,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) END DO END DO CASE ( np_CRV )                           !* Coriolis + relative vorticity DO jj = 1, jpjm1 DO ji = 1, fs_jpim1   ! vector opt. zwz(ji,jj) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)      & &                           - e1u(ji  ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  )   & &                        * r1_e1e2f(ji,jj)   ) END DO END DO CASE ( np_CME )                           !* Coriolis + metric DO jj = 1, jpjm1 DO ji = 1, fs_jpim1   ! vector opt. zwz(ji,jj) = (   ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & &                         - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   ) END DO END DO CASE DEFAULT                                             ! error CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  ) END SELECT ! IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==! DO jj = 1, jpjm1 DO ji = 1, fs_jpim1   ! vector opt. zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) END DO END DO ENDIF ! CALL lbc_lnk( zwz, 'F', 1. ) ! !                                   !==  horizontal fluxes  ==! zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,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          ! split in 2 parts due to vector opt. ztne(ji,jj) = ( zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1) ) / ( e3u_n(ji  , jj, jk) + e3v_n(ji, jj  , jk ) ) ztnw(ji,jj) = ( zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  ) ) / ( e3u_n(ji-1, jj, jk) + e3v_n(ji, jj  , jk ) ) ztse(ji,jj) = ( zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1) ) / ( e3u_n(ji  , jj, jk) + e3v_n(ji, jj-1, jk ) ) ztsw(ji,jj) = ( zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  ) ) / ( e3u_n(ji-1, jj, jk) + e3v_n(ji, jj-1, jk ) ) 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) ) / ( e3u_n(ji  , jj, jk) + e3v_n(ji, jj  , jk ) ) ztnw(ji,jj) = ( zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  ) ) / ( e3u_n(ji-1, jj, jk) + e3v_n(ji, jj  , jk ) ) ztse(ji,jj) = ( zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1) ) / ( e3u_n(ji  , jj, jk) + e3v_n(ji, jj-1, jk ) ) ztsw(ji,jj) = ( zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  ) ) / ( e3u_n(ji-1, jj, jk) + e3v_n(ji, jj-1, jk ) ) END DO END DO DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1   ! vector opt. zua = + r1_6 * 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_6 * 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  ) ) 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 !                                                ! =============== END SUBROUTINE vor_eeUV SUBROUTINE dyn_vor_init !!--------------------------------------------------------------------- !! NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT,   & &                 ln_dynvor_een, nn_een_e3f   , ln_dynvor_mix, ln_dynvor_msk &                 ln_dynvor_eeUV, ln_dynvor_een, nn_een_e3f   , ln_dynvor_mix, ln_dynvor_msk !!---------------------------------------------------------------------- ! WRITE(numout,*) '      t-point energy conserving scheme               ln_dynvor_enT = ', ln_dynvor_enT WRITE(numout,*) '      energy conserving scheme  (een using e3t)      ln_dynvor_eeT = ', ln_dynvor_eeT WRITE(numout,*) '      energy conserving scheme  (een using e3u and e3v)  ln_dynvor_eeUV = ', ln_dynvor_eeUV WRITE(numout,*) '      enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f IF( ln_dynvor_enT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENT   ;   ENDIF IF( ln_dynvor_eeT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EET   ;   ENDIF IF( ln_dynvor_eeUV ) THEN  ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EEUV  ;   ENDIF IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EEN   ;   ENDIF IF( ln_dynvor_mix ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_MIX   ;   ENDIF CASE( np_ENT )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at T-points) (ENT)' CASE( np_EET )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3t) (EET)' CASE( np_EEUV )  ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3u and e3v) (EEUV)' CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)' CASE( np_MIX )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)'
Note: See TracChangeset for help on using the changeset viewer.