Changeset 10094


Ignore:
Timestamp:
2018-09-06T15:24:19+02:00 (22 months ago)
Author:
davestorkey
Message:

UKMO dev_r10037_dynvor_EEUV branch : add new dynvor scheme.

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

Legend:

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

    r9950 r10094  
    868868   ln_dynvor_enT = .false. !  energy conserving scheme (T-point) 
    869869   ln_dynvor_eeT = .false. !  energy conserving scheme (een using e3t) 
     870   ln_dynvor_eeUV = .false.!  energy conserving scheme (een using e3u and e3v) 
    870871   ln_dynvor_een = .false. !  energy & enstrophy scheme 
    871872      nn_een_e3f = 1          ! =0  e3f = mi(mj(e3t))/4  
  • NEMO/branches/UKMO/dev_r10037_dynvor_EEUV/src/OCE/DYN/dynspg_ts.F90

    r9950 r10094  
    104104      ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) ) 
    105105      ! 
    106       IF( ln_dynvor_een .OR. ln_dynvor_eeT )   & 
     106      IF( ln_dynvor_een .OR. ln_dynvor_eeT .OR. ln_dynvor_eeUV)   & 
    107107         &     ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &  
    108108         &               ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) ) 
     
    156156      REAL(wp) ::   za0, za1, za2, za3              !   -      - 
    157157      REAL(wp) ::   zmdi, zztmp            , z1_ht  !   -      - 
     158      REAL(wp) ::   zmask_ne,zmask_nw,zmask_se,zmask_sw 
     159      REAL(wp) ::   z1_huv_ne,z1_huv_nw,z1_huv_se,z1_huv_sw 
    158160      REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e, zhf 
    159161      REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zssh_frc 
     
    283285                  ftse(ji,jj) = ( ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht 
    284286                  ftsw(ji,jj) = ( ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) ) * z1_ht 
     287               END DO 
     288            END DO 
     289         CASE( np_EEUV )                  != EEN scheme using e3u and e3v (energy conserving scheme) 
     290            ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 
     291            DO jj = 2, jpj 
     292               DO ji = 2, jpi 
     293                  zmask_ne = MAX(ssumask(ji  ,jj),ssvmask(ji,jj  )) 
     294                  zmask_nw = MAX(ssumask(ji-1,jj),ssvmask(ji,jj  )) 
     295                  zmask_se = MAX(ssumask(ji  ,jj),ssvmask(ji,jj-1)) 
     296                  zmask_sw = MAX(ssumask(ji-1,jj),ssvmask(ji,jj-1)) 
     297                  z1_huv_ne = 2._wp * zmask_ne / ( hu_n(ji  ,jj) + hv_n(ji,jj  ) + 1._wp - zmask_ne ) 
     298                  z1_huv_nw = 2._wp * zmask_nw / ( hu_n(ji-1,jj) + hv_n(ji,jj  ) + 1._wp - zmask_nw ) 
     299                  z1_huv_se = 2._wp * zmask_se / ( hu_n(ji  ,jj) + hv_n(ji,jj-1) + 1._wp - zmask_se ) 
     300                  z1_huv_sw = 2._wp * zmask_sw / ( hu_n(ji-1,jj) + hv_n(ji,jj-1) + 1._wp - zmask_sw ) 
     301                  ftne(ji,jj) = ( ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) ) * z1_huv_ne 
     302                  ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) ) * z1_huv_nw 
     303                  ftse(ji,jj) = ( ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) ) * z1_huv_se 
     304                  ftsw(ji,jj) = ( ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) ) * z1_huv_sw 
    285305               END DO 
    286306            END DO 
     
    426446         END DO 
    427447         ! 
    428       CASE( np_EET , np_EEN )      ! energy & enstrophy scheme (using e3t or e3f)          
     448      CASE( np_EET , np_EEN,  np_EEUV )      ! energy & enstrophy scheme (using e3t or e3f or e3u and e3v)          
    429449         DO jj = 2, jpjm1 
    430450            DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    10351055            END DO 
    10361056            ! 
    1037          CASE( np_EET , np_EEN )   ! energy & enstrophy scheme (using e3t or e3f) 
     1057         CASE( np_EET , np_EEN,  np_EEUV )   ! energy & enstrophy scheme (using e3t or e3f or e3u and e3v) 
    10381058            DO jj = 2, jpjm1 
    10391059               DO ji = fs_2, fs_jpim1   ! vector opt. 
  • NEMO/branches/UKMO/dev_r10037_dynvor_EEUV/src/OCE/DYN/dynvor.F90

    r9950 r10094  
    5656   LOGICAL, PUBLIC ::   ln_dynvor_enT   !: t-point energy conserving scheme     (ENT) 
    5757   LOGICAL, PUBLIC ::   ln_dynvor_eeT   !: t-point energy conserving scheme     (EET) 
     58   LOGICAL, PUBLIC ::   ln_dynvor_eeUV  !: uv-point energy conserving scheme    (EEUV) 
    5859   LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy & enstrophy conserving scheme (EEN) 
    5960   INTEGER, PUBLIC ::      nn_een_e3f      !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 
     
    6768   INTEGER, PUBLIC, PARAMETER ::   np_ENT = 2   ! ENT scheme (t-point vorticity) 
    6869   INTEGER, PUBLIC, PARAMETER ::   np_EET = 3   ! EET scheme (EEN using e3t) 
    69    INTEGER, PUBLIC, PARAMETER ::   np_EEN = 4   ! EEN scheme 
    70    INTEGER, PUBLIC, PARAMETER ::   np_MIX = 5   ! MIX scheme 
     70   INTEGER, PUBLIC, PARAMETER ::   np_EEUV = 4  ! EET scheme (EEN using e3u and e3v) 
     71   INTEGER, PUBLIC, PARAMETER ::   np_EEN = 5   ! EEN scheme 
     72   INTEGER, PUBLIC, PARAMETER ::   np_MIX = 6   ! MIX scheme 
    7173 
    7274   INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity  
     
    8486    
    8587   REAL(wp) ::   r1_4  = 0.250_wp         ! =1/4 
     88   REAL(wp) ::   r1_6  = 1._wp / 6._wp    ! =1/6 
    8689   REAL(wp) ::   r1_8  = 0.125_wp         ! =1/8 
    8790   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! 1/12 
     
    128131         CASE( np_EET )           ;   CALL vor_eeT( kt, ncor, un , vn , ua, va )   ! energy conserving scheme (een with e3t) 
    129132            IF( ln_stcor )            CALL vor_eeT( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
     133         CASE( np_EEUV )          ;   CALL vor_eeUV( kt, ncor, un , vn , ua, va )  ! energy conserving scheme (een with e3u and e3v) 
     134            IF( ln_stcor )            CALL vor_eeUV( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend 
    130135         CASE( np_EEN )           ;   CALL vor_een( kt, ncor, un , vn , ua, va )   ! energy & enstrophy scheme 
    131136            IF( ln_stcor )            CALL vor_een( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
     
    141146            CASE( np_ENT )           ;   CALL vor_enT( kt, nrvm, un , vn , ua, va )  ! energy conserving scheme (T-pts) 
    142147            CASE( np_EET )           ;   CALL vor_eeT( kt, nrvm, un , vn , ua, va )  ! energy conserving scheme (een with e3t) 
     148            CASE( np_EEUV )          ;   CALL vor_eeUV( kt, nrvm, un , vn , ua, va ) ! energy conserving scheme (een with e3u and e3v) 
    143149            CASE( np_ENE )           ;   CALL vor_ene( kt, nrvm, un , vn , ua, va )  ! energy conserving scheme 
    144150            CASE( np_ENS, np_MIX )   ;   CALL vor_ens( kt, nrvm, un , vn , ua, va )  ! enstrophy conserving scheme 
     
    161167                             CALL vor_eeT( kt, ntot, un , vn , ua, va )   ! total vorticity trend 
    162168            IF( ln_stcor )   CALL vor_eeT( kt, ncor, usd, vsd, ua, va )   ! add the Stokes-Coriolis trend 
     169         CASE( np_EEUV )                       !* energy conserving scheme (een scheme using e3u and e3v) 
     170                             CALL vor_eeUV( kt, ntot, un , vn , ua, va )  ! total vorticity trend 
     171            IF( ln_stcor )   CALL vor_eeUV( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend 
    163172         CASE( np_ENE )                        !* energy conserving scheme 
    164173                             CALL vor_ene( kt, ntot, un , vn , ua, va )   ! total vorticity trend 
     
    806815 
    807816 
     817   SUBROUTINE vor_eeUV( kt, kvor, pun, pvn, pua, pva ) 
     818      !!---------------------------------------------------------------------- 
     819      !!                ***  ROUTINE vor_eeUV  *** 
     820      !! 
     821      !! ** Purpose :   Compute the now total vorticity trend and add it to  
     822      !!      the general trend of the momentum equation. 
     823      !! 
     824      !! ** Method  :   Trend evaluated using now fields (centered in time)  
     825      !!      and the Arakawa and Lamb (1980) flux form formulation : conserves  
     826      !!      both the horizontal kinetic energy and the potential enstrophy 
     827      !!      when horizontal divergence is zero (see the NEMO documentation) 
     828      !!      Add this trend to the general momentum trend (ua,va). 
     829      !! 
     830      !!      Modified version of EEN as suggested by Mike Bell.  
     831      !! 
     832      !! ** Action : - Update (ua,va) with the now vorticity term trend 
     833      !! 
     834      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 
     835      !!              Bell, Johnson and Marshall, unpublished notes 2017 
     836      !!---------------------------------------------------------------------- 
     837      INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index 
     838      INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric 
     839      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun, pvn    ! now velocities 
     840      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva    ! total v-trend 
     841      ! 
     842      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     843      INTEGER  ::   ierr         ! local integer 
     844      REAL(wp) ::   zua, zva     ! local scalars 
     845      REAL(wp) ::   zmsk, ze3f   ! local scalars 
     846      REAL(wp), DIMENSION(jpi,jpj) ::   zwx , zwy , zwz  
     847      REAL(wp), DIMENSION(jpi,jpj) ::   ztnw, ztne, ztsw, ztse 
     848      !!---------------------------------------------------------------------- 
     849      ! 
     850      IF( kt == nit000 ) THEN 
     851         IF(lwp) WRITE(numout,*) 
     852         IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme' 
     853         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
     854      ENDIF 
     855      ! 
     856      !                                                ! =============== 
     857      DO jk = 1, jpkm1                                 ! Horizontal slab 
     858         !                                             ! =============== 
     859         ! 
     860         SELECT CASE( kvor )                 !==  vorticity considered  ==! 
     861         CASE ( np_COR )                           !* Coriolis (planetary vorticity) 
     862            DO jj = 1, jpjm1 
     863               DO ji = 1, fs_jpim1   ! vector opt. 
     864                  zwz(ji,jj) = ff_f(ji,jj)  
     865               END DO 
     866            END DO 
     867         CASE ( np_RVO )                           !* relative vorticity 
     868            DO jj = 1, jpjm1 
     869               DO ji = 1, fs_jpim1   ! vector opt. 
     870                  zwz(ji,jj) = ( e2v(ji+1,jj  ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)  & 
     871                     &         - e1u(ji  ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
     872               END DO 
     873            END DO 
     874         CASE ( np_MET )                           !* metric term 
     875            DO jj = 1, jpjm1 
     876               DO ji = 1, fs_jpim1   ! vector opt. 
     877                  zwz(ji,jj) = (   ( pvn(ji+1,jj,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     878                     &           - ( pun(ji,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   )  
     879               END DO 
     880            END DO 
     881         CASE ( np_CRV )                           !* Coriolis + relative vorticity 
     882            DO jj = 1, jpjm1 
     883               DO ji = 1, fs_jpim1   ! vector opt. 
     884                  zwz(ji,jj) = (  ff_f(ji,jj) + (  e2v(ji+1,jj  ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk)      & 
     885                     &                           - e1u(ji  ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk)  )   & 
     886                     &                        * r1_e1e2f(ji,jj)   )  
     887               END DO 
     888            END DO 
     889         CASE ( np_CME )                           !* Coriolis + metric 
     890            DO jj = 1, jpjm1 
     891               DO ji = 1, fs_jpim1   ! vector opt. 
     892                  zwz(ji,jj) = (   ff_f(ji,jj) + ( pvn(ji+1,jj  ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj)   & 
     893                     &                         - ( pun(ji  ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)   )  
     894               END DO 
     895            END DO 
     896         CASE DEFAULT                                             ! error 
     897            CALL ctl_stop('STOP','dyn_vor: wrong value for kvor'  ) 
     898         END SELECT 
     899         ! 
     900         IF( ln_dynvor_msk ) THEN          !==  mask/unmask vorticity ==! 
     901            DO jj = 1, jpjm1 
     902               DO ji = 1, fs_jpim1   ! vector opt. 
     903                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 
     904               END DO 
     905            END DO 
     906         ENDIF 
     907         ! 
     908         CALL lbc_lnk( zwz, 'F', 1. ) 
     909         ! 
     910         !                                   !==  horizontal fluxes  ==! 
     911         zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 
     912         zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 
     913 
     914         !                                   !==  compute and add the vorticity term trend  =! 
     915         jj = 2 
     916         ztne(1,:) = 0   ;   ztnw(1,:) = 0   ;   ztse(1,:) = 0   ;   ztsw(1,:) = 0 
     917         DO ji = 2, jpi          ! split in 2 parts due to vector opt. 
     918               ztne(ji,jj) = ( zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1) ) / ( e3u_n(ji  , jj, jk) + e3v_n(ji, jj  , jk ) ) 
     919               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 ) ) 
     920               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 ) ) 
     921               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 ) ) 
     922         END DO 
     923         DO jj = 3, jpj 
     924            DO ji = fs_2, jpi   ! vector opt. ok because we start at jj = 3 
     925               ztne(ji,jj) = ( zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1) ) / ( e3u_n(ji  , jj, jk) + e3v_n(ji, jj  , jk ) ) 
     926               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 ) ) 
     927               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 ) ) 
     928               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 ) ) 
     929            END DO 
     930         END DO 
     931         DO jj = 2, jpjm1 
     932            DO ji = fs_2, fs_jpim1   ! vector opt. 
     933               zua = + r1_6 * r1_e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   & 
     934                  &                             + ztse(ji,jj  ) * zwy(ji  ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) ) 
     935               zva = - r1_6 * r1_e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   & 
     936                  &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) ) 
     937               pua(ji,jj,jk) = pua(ji,jj,jk) + zua 
     938               pva(ji,jj,jk) = pva(ji,jj,jk) + zva 
     939            END DO   
     940         END DO   
     941         !                                             ! =============== 
     942      END DO                                           !   End of slab 
     943      !                                                ! =============== 
     944   END SUBROUTINE vor_eeUV 
     945 
     946 
    808947   SUBROUTINE dyn_vor_init 
    809948      !!--------------------------------------------------------------------- 
     
    817956      !! 
    818957      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT,   & 
    819          &                 ln_dynvor_een, nn_een_e3f   , ln_dynvor_mix, ln_dynvor_msk 
     958         &                 ln_dynvor_eeUV, ln_dynvor_een, nn_een_e3f   , ln_dynvor_mix, ln_dynvor_msk 
    820959      !!---------------------------------------------------------------------- 
    821960      ! 
     
    840979         WRITE(numout,*) '      t-point energy conserving scheme               ln_dynvor_enT = ', ln_dynvor_enT 
    841980         WRITE(numout,*) '      energy conserving scheme  (een using e3t)      ln_dynvor_eeT = ', ln_dynvor_eeT 
     981         WRITE(numout,*) '      energy conserving scheme  (een using e3u and e3v)  ln_dynvor_eeUV = ', ln_dynvor_eeUV 
    842982         WRITE(numout,*) '      enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een 
    843983         WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f 
     
    8731013      IF( ln_dynvor_enT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENT   ;   ENDIF 
    8741014      IF( ln_dynvor_eeT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EET   ;   ENDIF 
     1015      IF( ln_dynvor_eeUV ) THEN  ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EEUV  ;   ENDIF 
    8751016      IF( ln_dynvor_een ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EEN   ;   ENDIF 
    8761017      IF( ln_dynvor_mix ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_MIX   ;   ENDIF 
     
    9251066         CASE( np_ENT )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at T-points) (ENT)' 
    9261067         CASE( np_EET )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3t) (EET)' 
     1068         CASE( np_EEUV )  ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3u and e3v) (EEUV)' 
    9271069         CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)' 
    9281070         CASE( np_MIX )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)' 
Note: See TracChangeset for help on using the changeset viewer.