Changeset 13513


Ignore:
Timestamp:
2020-09-24T14:17:18+02:00 (4 months ago)
Author:
techene
Message:

#2527 and #2385 add a symmetric operator for ocean viscosity

Location:
NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/cfgs/SHARED/namelist_ref

    r13286 r13513  
    10121012   !                       !  Type of the operator : 
    10131013   ln_dynldf_OFF = .false.     !  No operator (i.e. no explicit diffusion) 
     1014   nn_dynldf_typ = 0           !  =0 div-rot (default)   ;   =1 symmetric 
    10141015   ln_dynldf_lap = .false.     !    laplacian operator 
    10151016   ln_dynldf_blp = .false.     !  bilaplacian operator 
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynldf_lap_blp.F90

    r13295 r13513  
    55   !!====================================================================== 
    66   !! History : 3.7  ! 2014-01  (G. Madec, S. Masson)  Original code, re-entrant laplacian 
     7   !!           4.0  ! 2020-04  (A. Nasser, G. Madec)  Add symmetric mixing tensor  
    78   !!---------------------------------------------------------------------- 
    89 
     
    1920   USE in_out_manager ! I/O manager 
    2021   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    21  
     22   USE lib_mpp 
     23    
    2224   IMPLICIT NONE 
    2325   PRIVATE 
     
    4749      !! 
    4850      !! ** Action : - pu_rhs, pv_rhs increased by the harmonic operator applied on pu, pv. 
     51      !! 
     52      !! Reference : S.Griffies, R.Hallberg 2000 Mon.Wea.Rev., DOI:/  
    4953      !!---------------------------------------------------------------------- 
    5054      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index 
     
    5761      REAL(wp) ::   zsign        ! local scalars 
    5862      REAL(wp) ::   zua, zva     ! local scalars 
    59       REAL(wp), DIMENSION(jpi,jpj) ::   zcur, zdiv 
     63      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zcur, zdiv 
     64      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zten, zshe   ! tension (diagonal) and shearing (anti-diagonal) terms 
    6065      !!---------------------------------------------------------------------- 
    6166      ! 
     
    7075      ENDIF 
    7176      ! 
    72       !                                                ! =============== 
    73       DO jk = 1, jpkm1                                 ! Horizontal slab 
    74          !                                             ! =============== 
    75          DO_2D( 0, 1, 0, 1 ) 
    76             !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1) 
    77             zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &   ! ahmf already * by fmask 
    78                &     * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  & 
    79                &        - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  ) 
    80             !                                      ! ahm * div        (computed from 2 to jpi/jpj) 
    81             zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb)               &   ! ahmt already * by tmask 
    82                &     * (  e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk)  & 
    83                &        + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk)  ) 
    84          END_2D 
     77      SELECT CASE( nn_dynldf_typ )   
     78      !               
     79      CASE ( np_typ_rot )       !==  Vorticity-Divergence operator  ==! 
    8580         ! 
    86          DO_2D( 0, 0, 0, 0 ) 
    87             pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * umask(ji,jj,jk) * (    &    ! * by umask is mandatory for dyn_ldf_blp use 
    88                &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   & 
    89                &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                      ) 
     81         ALLOCATE( zcur(jpi,jpj) , zdiv(jpi,jpj) ) 
     82         ! 
     83         DO jk = 1, jpkm1                                 ! Horizontal slab 
     84            ! 
     85            DO_2D( 0, 1, 0, 1 ) 
     86               !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1) 
     87               zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &   ! ahmf already * by fmask 
     88                  &     * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  & 
     89                  &        - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  ) 
     90               !                                      ! ahm * div        (computed from 2 to jpi/jpj) 
     91               zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb)               &   ! ahmt already * by tmask 
     92                  &     * (  e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk)  & 
     93                  &        + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk)  ) 
     94            END_2D 
     95            ! 
     96            DO_2D( 0, 0, 0, 0 ) 
     97               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * umask(ji,jj,jk) * (    &    ! * by umask is mandatory for dyn_ldf_blp use 
     98                  &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   & 
     99                  &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                      ) 
    90100               ! 
    91             pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * vmask(ji,jj,jk) * (    &    ! * by vmask is mandatory for dyn_ldf_blp use 
    92                &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm)   & 
    93                &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                      ) 
    94          END_2D 
    95          !                                             ! =============== 
    96       END DO                                           !   End of slab 
    97       !                                                ! =============== 
     101               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * vmask(ji,jj,jk) * (    &    ! * by vmask is mandatory for dyn_ldf_blp use 
     102                  &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm)   & 
     103                  &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                      ) 
     104            END_2D 
     105            ! 
     106         END DO                                           !   End of slab 
     107         ! 
     108         DEALLOCATE( zcur , zdiv ) 
     109         ! 
     110      CASE ( np_typ_sym )       !==  Symmetric operator  ==! 
     111         ! 
     112         ALLOCATE( zten(jpi,jpj) , zshe(jpi,jpj) ) 
     113         ! 
     114         DO jk = 1, jpkm1                                 ! Horizontal slab 
     115            ! 
     116            DO_2D( 0, 1, 0, 1 ) 
     117               !                                      ! shearing stress component (F-point)   NB : ahmf has already been multiplied by fmask 
     118               zshe(ji-1,jj-1) = ahmf(ji-1,jj-1,jk)                                                              & 
     119                  &     * (    e1f(ji-1,jj-1)    * r1_e2f(ji-1,jj-1)                                             & 
     120                  &         * ( pu(ji-1,jj  ,jk) * r1_e1u(ji-1,jj  )  - pu(ji-1,jj-1,jk) * r1_e1u(ji-1,jj-1) )   & 
     121                  &         +  e2f(ji-1,jj-1)    * r1_e1f(ji-1,jj-1)                                             & 
     122                  &         * ( pv(ji  ,jj-1,jk) * r1_e2v(ji  ,jj-1)  - pv(ji-1,jj-1,jk) * r1_e2v(ji-1,jj-1) )   )  
     123               !                                      ! tension stress component (T-point)   NB : ahmt has already been multiplied by tmask 
     124               zten(ji,jj)    = ahmt(ji,jj,jk)                                                       & 
     125                  &     * (    e2t(ji,jj)    * r1_e1t(ji,jj)                                         & 
     126                  &         * ( pu(ji,jj,jk) * r1_e2u(ji,jj)  - pu(ji-1,jj,jk) * r1_e2u(ji-1,jj) )   & 
     127                  &         -  e1t(ji,jj)    * r1_e2t(ji,jj)                                         & 
     128                  &         * ( pv(ji,jj,jk) * r1_e1v(ji,jj)  - pv(ji,jj-1,jk) * r1_e1v(ji,jj-1) )   )    
     129            END_2D 
     130            ! 
     131            DO_2D( 0, 0, 0, 0 ) 
     132               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)                               & 
     133                  &    * (   (   zten(ji+1,jj  ) * e2t(ji+1,jj  )*e2t(ji+1,jj  ) * e3t(ji+1,jj  ,jk,Kmm)                       & 
     134                  &            - zten(ji  ,jj  ) * e2t(ji  ,jj  )*e2t(ji  ,jj  ) * e3t(ji  ,jj  ,jk,Kmm) ) * r1_e2u(ji,jj)     &                                                     
     135                  &        + (   zshe(ji  ,jj  ) * e1f(ji  ,jj  )*e1f(ji  ,jj  ) * e3f(ji  ,jj  ,jk)                           & 
     136                  &            - zshe(ji  ,jj-1) * e1f(ji  ,jj-1)*e1f(ji  ,jj-1) * e3f(ji  ,jj-1,jk)     ) * r1_e1u(ji,jj) )    
     137               ! 
     138               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)                               & 
     139                  &    * (   (   zshe(ji  ,jj  ) * e2f(ji  ,jj  )*e2f(ji  ,jj  ) * e3f(ji  ,jj  ,jk)                           & 
     140                  &            - zshe(ji-1,jj  ) * e2f(ji-1,jj  )*e2f(ji-1,jj  ) * e3f(ji-1,jj  ,jk)     ) * r1_e2v(ji,jj)     & 
     141                  &        - (   zten(ji  ,jj+1) * e1t(ji  ,jj+1)*e1t(ji  ,jj+1) * e3t(ji  ,jj+1,jk,Kmm)                       & 
     142                  &            - zten(ji  ,jj  ) * e1t(ji  ,jj  )*e1t(ji  ,jj  ) * e3t(ji  ,jj  ,jk,Kmm) ) * r1_e1v(ji,jj) ) 
     143               ! 
     144            END_2D 
     145            ! 
     146         END DO 
     147         ! 
     148         DEALLOCATE( zten , zshe ) 
     149         ! 
     150      END SELECT 
    98151      ! 
    99152   END SUBROUTINE dyn_ldf_lap 
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/LDF/ldfdyn.F90

    r13295 r13513  
    3434   !                                    !!* Namelist namdyn_ldf : lateral mixing on momentum * 
    3535   LOGICAL , PUBLIC ::   ln_dynldf_OFF   !: No operator (i.e. no explicit diffusion) 
     36   INTEGER , PUBLIC ::   nn_dynldf_typ   !: operator type (0: div-rot ; 1: symmetric) 
    3637   LOGICAL , PUBLIC ::   ln_dynldf_lap   !: laplacian operator 
    3738   LOGICAL , PUBLIC ::   ln_dynldf_blp   !: bilaplacian operator 
     
    5253 
    5354   !                                    !!* Parameter to control the type of lateral viscous operator 
    54    INTEGER, PARAMETER, PUBLIC ::   np_ERROR  =-10                       !: error in setting the operator 
    55    INTEGER, PARAMETER, PUBLIC ::   np_no_ldf = 00                       !: without operator (i.e. no lateral viscous trend) 
     55   INTEGER, PARAMETER, PUBLIC ::   np_ERROR   =-10                      !: error in setting the operator 
     56   INTEGER, PARAMETER, PUBLIC ::   np_no_ldf  = 00                      !: without operator (i.e. no lateral viscous trend) 
     57   ! 
     58   INTEGER, PARAMETER, PUBLIC ::   np_typ_rot = 0                       !: div-rot   operator 
     59   INTEGER, PARAMETER, PUBLIC ::   np_typ_sym = 1                       !: symmetric operator 
     60   ! 
    5661   !                          !!      laplacian     !    bilaplacian    ! 
    5762   INTEGER, PARAMETER, PUBLIC ::   np_lap    = 10   ,   np_blp    = 20  !: iso-level operator 
     
    109114      CHARACTER(len=5) ::   cl_Units               ! units (m2/s or m4/s) 
    110115      !! 
    111       NAMELIST/namdyn_ldf/ ln_dynldf_OFF, ln_dynldf_lap, ln_dynldf_blp,   &   ! type of operator 
    112          &                 ln_dynldf_lev, ln_dynldf_hor, ln_dynldf_iso,   &   ! acting direction of the operator 
    113          &                 nn_ahm_ijk_t , rn_Uv    , rn_Lv,   rn_ahm_b,   &   ! lateral eddy coefficient 
    114          &                 rn_csmc      , rn_minfac    , rn_maxfac            ! Smagorinsky settings 
     116      NAMELIST/namdyn_ldf/ ln_dynldf_OFF, nn_dynldf_typ, ln_dynldf_lap, ln_dynldf_blp,   &   ! type of operator 
     117         &                 ln_dynldf_lev, ln_dynldf_hor, ln_dynldf_iso,                  &   ! acting direction of the operator 
     118         &                 nn_ahm_ijk_t , rn_Uv        , rn_Lv        ,   rn_ahm_b,      &   ! lateral eddy coefficient 
     119         &                 rn_csmc      , rn_minfac    , rn_maxfac                           ! Smagorinsky settings 
    115120      !!---------------------------------------------------------------------- 
    116121      ! 
     
    130135         WRITE(numout,*) '      type :' 
    131136         WRITE(numout,*) '         no explicit diffusion                ln_dynldf_OFF = ', ln_dynldf_OFF 
     137         WRITE(numout,*) '         type of operator (div-rot or sym)    nn_dynldf_typ = ', nn_dynldf_typ 
    132138         WRITE(numout,*) '         laplacian operator                   ln_dynldf_lap = ', ln_dynldf_lap 
    133139         WRITE(numout,*) '         bilaplacian operator                 ln_dynldf_blp = ', ln_dynldf_blp 
     
    147153         WRITE(numout,*) '         Smagorinsky coefficient              rn_csmc       = ', rn_csmc 
    148154         WRITE(numout,*) '         factor multiplier for eddy visc.' 
    149          WRITE(numout,*) '            lower limit (default 1.0)         rn_minfac    = ', rn_minfac 
    150          WRITE(numout,*) '            upper limit (default 1.0)         rn_maxfac    = ', rn_maxfac 
     155         WRITE(numout,*) '            lower limit (default 1.0)         rn_minfac     = ', rn_minfac 
     156         WRITE(numout,*) '            upper limit (default 1.0)         rn_maxfac     = ', rn_maxfac 
    151157      ENDIF 
    152158 
     
    160166      IF( ln_dynldf_lap ) THEN   ;                              ioptio = ioptio + 1   ;   ENDIF 
    161167      IF( ln_dynldf_blp ) THEN   ;                              ioptio = ioptio + 1   ;   ENDIF 
    162       IF( ioptio /= 1   )   CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 operator options (NONE/lap/blp)' ) 
     168      IF( ioptio /= 1   )   CALL ctl_stop( 'ldf_dyn_init: use ONE of the 3 operator options (NONE/lap/blp)' ) 
    163169      ! 
    164170      IF(.NOT.ln_dynldf_OFF ) THEN     !==  direction ==>> type of operator  ==! 
     171         ! 
     172         SELECT CASE( nn_dynldf_typ )  ! div-rot or symmetric 
     173         CASE( np_typ_rot )   ;   WRITE(numout,*) '   ==>>>   use div-rot   operator ' 
     174         CASE( np_typ_sym )   ;   WRITE(numout,*) '   ==>>>   use symmetric operator ' 
     175         CASE DEFAULT                                     ! error 
     176            CALL ctl_stop('ldf_dyn_init: wrong value for nn_dynldf_typ (0 or 1)'  ) 
     177         END SELECT 
     178         ! 
    165179         ioptio = 0 
    166180         IF( ln_dynldf_lev )   ioptio = ioptio + 1 
    167181         IF( ln_dynldf_hor )   ioptio = ioptio + 1 
    168182         IF( ln_dynldf_iso )   ioptio = ioptio + 1 
    169          IF( ioptio /= 1   )   CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 direction options (level/hor/iso)' ) 
     183         IF( ioptio /= 1   )   CALL ctl_stop( 'ldf_dyn_init: use ONE of the 3 direction options (level/hor/iso)' ) 
    170184         ! 
    171185         !                             ! Set nldf_dyn, the type of lateral diffusion, from ln_dynldf_... logicals 
Note: See TracChangeset for help on using the changeset viewer.