Changeset 13005


Ignore:
Timestamp:
2020-06-02T17:46:26+02:00 (4 months ago)
Author:
gm
Message:

ADE and more options to AM98 config

Location:
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater
Files:
12 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/EXPREF/namelist_cfg

    r12667 r13005  
    2323   nn_it000    =       1   !  first time step 
    2424   nn_itend    =  1        !  10 ans -  2 min - last  time step 
    25 !   nn_itend    =  2592000  !  10 ans -  2 min - last  time step 
    2625!   nn_itend    =  1036800  !  10 ans -  5 min 
    27 !   nn_itend    =    34560  !  1  an  - 15 min 
    28 !   nn_itend    =   172800  !  10 ans  - 30 min 
     26!   nn_itend    =   345600  !  10 ans  - 15 min - AM98 1/4deg + P.Marchand 
     27!   nn_itend    =   172800  !  10 ans  - 30 min - AM98 1/4deg 
    2928!   nn_itend    =    34560  !  10 ans  - 2h 
    3029!   nn_itend    =    86400  !  10 ans  - 1h 
     
    4342&namusr_def    !   AM98 user defined namelist   
    4443!----------------------------------------------------------------------- 
    45    nn_AM98     =     4     !  AM98 resolution [1/degrees] 
     44   nn_AM98     =     4     !  AM98 resolution [1/nn_AM98] 
    4645   jpkglo      =     2     !  number of model levels 
    4746   rn_theta    =    0.    !  rotation angle fo the grid [deg] 
     47   ! 
     48   nn_gc       =        2        ! number of ghostcells  
     49   rn_domsiz   =  2000000        ! size of the domain (default 2000km)  [m] 
     50   rn_dx       =   100000        ! gridspacing (default 100km)          [m] 
    4851   ! 
    4952   rn_tau           = 0.20    ! wind stress on the surface    [N/m2] 
     
    5356   rn_beta = 2.e-11   ! beta-plan coriolis parameter  [1/m.s]  
    5457   rn_f0   = 0.7e-4   ! f-plan coriolis parameter     [1/s]  
     58   ! 
     59   nn_dynldf_lap_typ = 1        ! choose type of laplacian (ideally from namelist) 
     60   !                            !       = 1   divrot    laplacian 
     61   !                            !       = 2   symmetric laplacian (Griffies&Hallberg 2000) 
     62   !                            !       = 3   symmetric laplacian (cartesian) 
     63   ln_dynldf_lap_PM  =.false.   ! if True - apply the P.Marchand boundary condition on the laplacian viscosity 
     64   ! 
    5565/ 
    5666!----------------------------------------------------------------------- 
     
    6070   ! 
    6171!   rn_Dt      = 7200.     !  2h    - 1 deg cfg - time step for the dynamics [s] 
    62 !   rn_Dt      = 3600.     !  1h    - 1 deg cfg  
    63    rn_Dt      = 1800.     !  30min    - 1/4 deg cfg 
    64 !   rn_Dt      = 900.      !  15min - 1/4 deg cfg 
     72!   rn_Dt      = 3600.     !  1h     - AM98 1 deg 
     73   rn_Dt      = 1800.     !  30min    - AM98 1/4deg 
     74!   rn_Dt      = 900.      !  15min - AM98 1/4deg + P.Marchand 
    6575!   rn_Dt      = 600.      !  10min - 1/4 deg cfg 
    6676!   rn_Dt      = 300.      !  5min - 1/4 deg cfg  
     
    109119&namlbc        !   lateral momentum boundary condition                  (default: NO selection) 
    110120!----------------------------------------------------------------------- 
    111    rn_shlat    =    0.     !  free slip 
     121   !                       !  free slip  !   partial slip  !   no slip   ! strong slip 
     122   rn_shlat    =    0.     !  shlat = 0  !  0 < shlat < 2  !  shlat = 2  !  2 < shlat 
    112123/ 
    113124!!====================================================================== 
     
    141152&namdyn_adv    !   formulation of the momentum advection                (default: NO selection) 
    142153!----------------------------------------------------------------------- 
    143    ln_dynadv_vec = .true.  !  vector form - 2nd centered scheme 
     154   ln_dynadv_OFF = .false. !  linear dynamics (no momentum advection) 
     155   ln_dynadv_vec = .false. !  vector form - 2nd centered scheme 
    144156     nn_dynkeg     = 0        ! grad(KE) scheme: =0   C2  ;  =1   Hollingsworth correction 
     157   ln_dynadv_cen2 = .false. !  flux form - 2nd order centered scheme 
     158   ln_dynadv_ubs = .false. !  flux form - 3rd order UBS      scheme 
    145159/ 
    146160!----------------------------------------------------------------------- 
    147161&namdyn_vor    !   Vorticity / Coriolis scheme                          (default: NO selection) 
    148162!----------------------------------------------------------------------- 
    149    ln_dynvor_ene = .true.   !  energy conserving scheme 
    150    ln_dynvor_ens = .false.  !  enstrophy conserving scheme 
    151    ln_dynvor_een = .false.  !  energy & enstrophy scheme 
     163   ln_dynvor_ene = .false. !  energy    conserving scheme 
     164   ln_dynvor_ens = .false. !  enstrophy conserving scheme 
     165   ln_dynvor_enT = .false. !  energy conserving scheme (T-point) 
     166   ln_dynvor_eeT = .false. !  energy conserving scheme (een using e3t) 
     167   ln_dynvor_een = .false. !  energy & enstrophy scheme 
    152168      nn_een_e3f = 0          ! =0  e3f = mi(mj(e3t))/4 
    153169      !                       ! =1  e3f = mi(mj(e3t))/mi(mj( tmask)) 
    154    ln_dynvor_msk = .false.  !  vorticity multiplied by fmask (=T) 
    155       !                     !  (f-point vorticity schemes only) 
     170   ln_dynvor_msk = .false. !  vorticity multiplied by fmask (=T)         
     171      !                    !  (f-point vorticity schemes only) 
     172      ! 
     173   ln_dynvor_ens_adVO     = .false.   ! Alternative Direction 
     174   ln_dynvor_ens_adKE     = .false. 
     175   ln_dynvor_ens_adKEVO   = .false. 
     176   ln_dynvor_ene_adVO     = .false. 
     177   ln_dynvor_ene_adKE     = .false. 
     178   ln_dynvor_ene_adKEVO   = .false. 
    156179/ 
    157180!----------------------------------------------------------------------- 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/usrdef_hgr.F90

    r12614 r13005  
    6868      INTEGER  ::   ji, jj               ! dummy loop indices 
    6969      REAL(wp) ::   zlam1, zlam0, zcos_theta, zim1 , zjm1 , ze1  , ze1deg ! local scalars 
    70       REAL(wp) ::   zphi1, zphi0, zsin_theta, zim05, zjm05, znorme      !   -      - 
    71       REAL(wp) ::   zl, zgl, zbl      !   -      - 
     70      REAL(wp) ::   zphi1, zphi0, zsin_theta, zim05, zjm05, znorme        !   -      - 
     71      REAL(wp) ::   zgl, zbl      !   -      - 
     72 
    7273      !!------------------------------------------------------------------------------- 
    7374      ! 
     
    8182      !                       !==  grid point position  ==! 
    8283      ! 
    83       ze1 =  100000._wp / REAL(nn_AM98, wp)       ! (100km) gridspacing              [m] 
    84        zl = 2000000._wp                           ! lentgh side square domain        [m]   
    85       zgl = 2000000._wp + 4._wp*ze1               ! length side square + ghost cells [m] 
    86                                                   !    2000km x 2000km + 2 cells around sides 
    87       !  
    88       ! biggest L = 3500 km (worst case) 
    89       !zbl = 3500000._wp                           ! length side bigger domain [m] 
    90       ! non rotated (0deg) - worst case (centered in the 3500kmx3500km) 
    91       !zlam0 = 875000._wp    ! m 
    92       !zphi0 = 875000._wp    ! m 
    93       ! rotated case (45eg) 
    94       ! origin (in meters)  of the 2000km x 2000km domain  
    95       !zlam0 = - 1000000._wp            ! - 1000km 
    96       !zphi0 =   1474873.7341529163_wp  ! 3500*sin(45deg) - 1000km 
    97  
    98  
     84      ze1 =  rn_dx / REAL(nn_AM98, wp)                   ! [m] gridspacing used 
     85      zgl =  rn_domsiz + 2._wp * REAL(nn_gc, wp) * ze1   ! [m] length of the square with ghostcells 
    9986      ! fit the best square around the square + ghost cells  
    10087      zbl = zgl * ( COS( rn_theta * rad ) + SIN( rn_theta * rad ) )   ! length side bigger domain [m]                                  
     
    11198 
    11299      ! exact origin in meters 
    113       zlam1 =  zbl * COS((rn_theta + 45 )* rad ) / SQRT( 2._wp )  - zl/2._wp  
    114       zphi1 =  zbl * SIN((rn_theta + 45 )* rad ) / SQRT( 2._wp )  - zl/2._wp  
     100      zlam1 =  zbl * COS((rn_theta + 45 )* rad ) / SQRT( 2._wp )  - rn_domsiz/2._wp  
     101      zphi1 =  zbl * SIN((rn_theta + 45 )* rad ) / SQRT( 2._wp )  - rn_domsiz/2._wp  
    115102      ! origin put in the true corner of a cell so there will be no cropping  
    116103      ! of the edge cells 
    117       zlam0 = REAL( anint( zlam1 / ze1 ) + 1, wp ) * ze1 
    118       zphi0 = REAl( anint( zphi1 / ze1 ) + 1, wp ) * ze1 
     104      zlam0 = REAL( anint( zlam1 / ze1 ), wp ) * ze1 
     105      zphi0 = REAl( anint( zphi1 / ze1 ), wp ) * ze1 
    119106 
    120107      IF(lwp) WRITE(numout,*) '                  origin position    zlam0   = ', zlam0/1000,   ' km' 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/usrdef_nam.F90

    r12614 r13005  
    2828   INTEGER , PUBLIC ::   nn_AM98            ! 1/nn_AM98 = the resolution chosen in degrees and thus defining the horizontal domain size 
    2929   REAL(wp), PUBLIC ::   rn_theta           ! rotation angle (in degree) of the grid  
    30    REAL(wp), PUBLIC ::   rn_tau             ! wind stress on the surface    [N/m2] 
    31    REAL(wp), PUBLIC ::   rn_f0              !    f-plan coriolis parameter  [1/s]  
    32    REAL(wp), PUBLIC ::   rn_beta            ! beta-plan coriolis parameter  [1/m.s]  
    33    REAL(wp), PUBLIC ::   rn_modified_grav   ! modified gravity              [m/s2]  
    34    REAL(wp), PUBLIC ::   rn_rfr             ! layer friction                [1/s] 
     30   INTEGER , PUBLIC ::   nn_gc              ! number of ghostcells 
     31   REAL(wp), PUBLIC ::   rn_domsiz          ! size of the domain (default 2000km)  [m] 
     32   REAL(wp), PUBLIC ::   rn_dx              ! gridspacing (default 100km)          [m] 
     33   REAL(wp), PUBLIC ::   rn_tau             ! wind stress on the surface           [N/m2] 
     34   REAL(wp), PUBLIC ::   rn_f0              !    f-plan coriolis parameter         [1/s]  
     35   REAL(wp), PUBLIC ::   rn_beta            ! beta-plan coriolis parameter         [1/m.s]  
     36   REAL(wp), PUBLIC ::   rn_modified_grav   ! modified gravity                     [m/s2]  
     37   REAL(wp), PUBLIC ::   rn_rfr             ! layer friction                       [1/s] 
     38   !                              !!* temporary *!! 
     39   INTEGER , PUBLIC ::   nn_dynldf_lap_typ            ! choose type of laplacian (ideally from namelist) 
     40   !                                                  ! = 1   divrot    laplacian 
     41   !                                                  ! = 2   symmetric laplacian (Griffies&Hallberg 2000) 
     42   !                                                  ! = 3   symmetric laplacian (cartesian) 
     43   LOGICAL , PUBLIC ::   ln_dynldf_lap_PM             ! if True - apply the P.Marchand boundary condition on the laplacian 
    3544   ! 
    3645   !!---------------------------------------------------------------------- 
     
    6170      REAL(wp) ::   ze1, zgl, zbl   ! gridspacing, length of the biggest square 
    6271      !! 
    63       NAMELIST/namusr_def/ nn_AM98, rn_theta, jpkglo, rn_tau,   &   !  
     72      NAMELIST/namusr_def/ nn_AM98, rn_theta, jpkglo,           &   !  
     73         &                 nn_gc ,rn_domsiz, rn_dx,             &   ! domain parameters 
    6474         &                 rn_f0 ,rn_beta,                      &   ! coriolis parameter 
    65          &                 rn_modified_grav, rn_rfr                 ! 
     75         &                 rn_modified_grav, rn_rfr , rn_tau,   &   ! reduced gravity, friction, wind 
     76         &                 nn_dynldf_lap_typ, ln_dynldf_lap_PM      ! temporary parameter 
    6677      !!---------------------------------------------------------------------- 
    6778      ! 
     
    7889 
    7990      kk_cfg = nn_AM98 
    80        
    81       ! square domain of 2000km x 2000km with dx = dy = 100km  
    82       ! kpi = 20 * nn_AM98 + 2        ! Global Domain size 
    83       ! kpj = 20 * nn_AM98 + 2 
    84       
    85       ! Number of cells 
    86       ! lenght of the biggest square + 2 ghost cells on each side  
    87       ze1 =  100000._wp / REAL(nn_AM98, wp)       ! (100km) gridspacing [m] 
    88       zgl = 2000000._wp + 4._wp*ze1               ! length side square + ghostcells [m] 
    89       !zbl = zgl / COS( rn_theta * rad )           ! length side bigger domain [m]   
     91       ! 
     92      ze1 =  rn_dx / REAL(nn_AM98, wp)                   ! [m] gridspacing used 
     93      zgl =  rn_domsiz + 2._wp * REAL(nn_gc, wp) * ze1   ! [m] length of the square with ghostcells 
     94      ! rotation 
    9095      zbl = zgl * ( COS( rn_theta * rad ) + SIN( rn_theta * rad ) )   ! length side bigger domain [m]   
    91          
     96      ! 
    9297      kpi = ceiling(zbl / ze1 )    ! Global Domain size + ghost cells                 
    9398      kpj = ceiling(zbl / ze1 )    ! Global Domain size + ghost cells               
     
    104109      IF( rn_modified_grav /= 0._wp) grav = rn_modified_grav   ! update the gravity  
    105110      ! 
     111      !        !== Temporary namelist parameter ==! 
     112      ! 
     113      ! ll_dynldf_lap_PM = ln_dynldf_lap_PM 
     114      ! 
    106115      kpk = jpkglo 
    107116      !                             ! Set the lateral boundary condition of the global domain 
     
    114123         WRITE(numout,*) '~~~~~~~~~~~ ' 
    115124         WRITE(numout,*) '   Namelist namusr_def : AM98 case' 
     125         WRITE(numout,*) '                                   domain size       rn_domsiz   = ', rn_domsiz, 'm' 
     126         WRITE(numout,*) '                                   gridspacing           rn_dx   = ', rn_dx, 'm' 
    116127         WRITE(numout,*) '      inverse resolution & implied domain size         nn_AM98   = ', nn_AM98 
     128         WRITE(numout,*) '                           implied gridspacing           rn_dx   = ', rn_dx, 'm' 
     129         WRITE(numout,*) '                          number of ghostcells           nn_gc   = ', nn_gc 
     130         WRITE(numout,*) '   ' 
     131         WRITE(numout,*) '                    rotation angle chosen             rn_theta   = ', rn_theta, 'deg' 
     132         WRITE(numout,*) '                    modified gravity          rn_modified_grav   = ', rn_modified_grav, 'm2/s' 
    117133#if defined key_agrif 
    118134         IF( Agrif_Root() ) THEN 
    119135#endif 
    120          WRITE(numout,*) '         jpiglo = Nx*nn_AM98+4                            jpiglo = ', kpi 
    121          WRITE(numout,*) '         jpjglo = Ny*nn_AM98+4                            jpjglo = ', kpj 
     136         WRITE(numout,*) '         jpiglo = Nx*nn_AM98 + 2*n_ghost                    jpiglo = ', kpi 
     137         WRITE(numout,*) '         jpjglo = Nx*nn_AM98 + 2*n_ghost                    jpjglo = ', kpj 
    122138#if defined key_agrif 
    123139         ENDIF 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/usrdef_sbc.F90

    r12614 r13005  
    9292      END DO 
    9393      
    94       
    9594      ! module of wind stress and wind speed at T-point 
    9695      zcoef = 1. / ( zrhoa * zcdrag )  
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/usrdef_zgr.F90

    r12614 r13005  
    33   !!                       ***  MODULE  usrdef_zgr  *** 
    44   !! 
    5    !!                       ===  GYRE configuration  === 
     5   !!                       ===  AM98 configuration  === 
    66   !! 
    77   !! User defined : vertical coordinate system of a user configuration 
     
    6464      ! 
    6565      IF(lwp) WRITE(numout,*) 
    66       IF(lwp) WRITE(numout,*) 'usr_def_zgr : GYRE configuration (z-coordinate closed flat box ocean without cavities)' 
     66      IF(lwp) WRITE(numout,*) 'usr_def_zgr : AM98 configuration (z-coordinate closed flat box ocean without cavities)' 
    6767      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    6868      ! 
     
    7070      ! type of vertical coordinate 
    7171      ! --------------------------- 
    72       ld_zco    = .FALSE.         ! GYRE case:  z-coordinate without ocean cavities 
     72      ld_zco    = .FALSE.         ! AM98 case:  z-coordinate without ocean cavities 
    7373      ld_zps    = .FALSE. 
    7474      ld_sco    = .TRUE. 
     
    170170      !! ** Purpose :   set the masked top and bottom ocean t-levels 
    171171      !! 
    172       !! ** Method  :   GYRE case = closed flat box ocean without ocean cavities 
     172      !! ** Method  :   AM98 case = closed flat box ocean without ocean cavities 
    173173      !!                   k_top = 1     except along north, south, east and west boundaries 
    174174      !!                   k_bot = jpk-1 except along north, south, east and west boundaries 
     
    182182      INTEGER  ::   ji, jj                    ! dummy loop indices 
    183183      REAL(wp) ::   zylim0, zylim1, zxlim0, zxlim1 ! limit of the domain [m] 
     184      REAL(WP) ::   zcoeff    ! local scalar 
    184185      !!---------------------------------------------------------------------- 
    185186      ! 
     
    187188      IF(lwp) WRITE(numout,*) '    zgr_top_bot : defines the top and bottom wet ocean levels.' 
    188189      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~' 
    189       IF(lwp) WRITE(numout,*) '       GYRE case : closed flat box ocean without ocean cavities' 
     190      IF(lwp) WRITE(numout,*) '       AM98 case : closed flat box ocean without ocean cavities' 
    190191      ! 
    191192      z2d(:,:) = REAL( jpkm1 , wp )          ! flat bottom 
     
    193194      CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. )           ! set surrounding land to zero (here jperio=0 ==>> closed) 
    194195      ! 
    195  
    196196      ! rotated domain (45deg) 
    197197      !zylim0 =  100000._wp * SQRT( 2._wp ) / 3    !  dy*sqrt(2)/2 * 2/3  
     
    214214      DO jj = 1, jpj 
    215215         DO ji = 1, jpi 
    216           
    217216         ! if T point in the 2000 km x 2000 km domain 
    218217         ! IF ( gphit(ji,jj) > zylim0 .AND. gphit(ji,jj) < zylim1 .AND. &  
     
    227226         k_bot(ji,jj) = 0 
    228227         END IF 
    229  
    230228         END DO 
    231229      END DO 
    232  
     230      ! mask the lonely corners 
     231      DO jj = 1, jpj 
     232         DO ji = 1, jpi 
     233         zcoeff = k_top(ji+1,jj) + k_top(ji,jj+1)   & 
     234            +     k_top(ji-1,jj) + k_top(ji,jj-1) 
     235         IF ( zcoeff <= 1._wp )   THEN 
     236            k_top(ji,jj) = 0    ! = land 
     237            k_bot(ji,jj) = 0 
     238         END IF 
     239         END DO 
     240      END DO 
    233241      ! k_bot(:,:) = NINT( z2d(:,:) )           ! =jpkm1 over the ocean point, =0 elsewhere 
    234242      ! 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domvvl.F90

    r12614 r13005  
    723723               &                    + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj) 
    724724         END_2D 
     725!!an   dans le cas tourné, hf augmente et trend VOR diminue 
     726!         DO_2D_10_10 
     727!            zc3(ji,jj) =           (  e1e2t(ji  ,jj  ) * pssh(ji  ,jj  )  & 
     728!               &                    + e1e2t(ji+1,jj  ) * pssh(ji+1,jj  )  & 
     729!               &                    + e1e2t(ji  ,jj+1) * pssh(ji  ,jj+1)  & 
     730!               &                    + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj) & 
     731!               &           / MAX( tmask(ji,jj) + tmask(ji+1,jj) + tmask(ji,jj+1) + tmask(ji+1,jj+1), 1._wp )  
     732!         END_2D 
     733          
    725734         CALL lbc_lnk( 'domvvl', zc3(:,:), 'F', 1._wp ) 
    726735         ! 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynadv.F90

    r12822 r13005  
    7474      CASE( np_VEC_c2  )      
    7575         CALL dyn_keg     ( kt, nn_dynkeg,      Kmm, puu, pvv, Krhs )    ! vector form : horizontal gradient of kinetic energy 
     76!!an SWE : w = 0 
    7677         CALL dyn_zad     ( kt,                 Kmm, puu, pvv, Krhs )    ! vector form : vertical advection 
    7778      CASE( np_FLX_c2  )  
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynkeg.F90

    r12822 r13005  
    2929 
    3030   PUBLIC   dyn_keg    ! routine called by step module 
     31   PUBLIC   dyn_kegAD   ! routine called by step module 
    3132    
    3233   INTEGER, PARAMETER, PUBLIC  ::   nkeg_C2     = 0   !: 2nd order centered scheme (standard scheme) 
     
    155156      ! 
    156157   END SUBROUTINE dyn_keg 
    157  
     158    
     159    
     160   SUBROUTINE dyn_kegAD( kt, kscheme, puu, pvv, pu_rhs, pv_rhs ) 
     161      !!---------------------------------------------------------------------- 
     162      !!                  ***  ROUTINE dyn_kegAD  *** 
     163      !! 
     164      !! ** Purpose :   Compute the now momentum trend due to the horizontal 
     165      !!      gradient of the horizontal kinetic energy and add it to the  
     166      !!      general momentum trend. 
     167      !! 
     168      !! ** Method  : * kscheme = nkeg_C2 : 2nd order centered scheme that  
     169      !!      conserve kinetic energy. Compute the now horizontal kinetic energy  
     170      !!         zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ] 
     171      !!              * kscheme = nkeg_HW : Hollingsworth correction following 
     172      !!      Arakawa (2001). The now horizontal kinetic energy is given by: 
     173      !!         zhke = 1/6 [ mi-1(  2 * un^2 + ((u(j+1)+u(j-1))/2)^2  ) 
     174      !!                    + mj-1(  2 * vn^2 + ((v(i+1)+v(i-1))/2)^2  ) ] 
     175      !!       
     176      !!      Take its horizontal gradient and add it to the general momentum 
     177      !!      trend. 
     178      !!         u(rhs) = u(rhs) - 1/e1u di[ zhke ] 
     179      !!         v(rhs) = v(rhs) - 1/e2v dj[ zhke ] 
     180      !! 
     181      !! ** Action : - Update the (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) with the hor. ke gradient trend 
     182      !!             - send this trends to trd_dyn (l_trddyn=T) for post-processing 
     183      !! 
     184      !! ** References : Arakawa, A., International Geophysics 2001. 
     185      !!                 Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983. 
     186      !!---------------------------------------------------------------------- 
     187      ! 
     188      INTEGER                                  , INTENT( in )  ::  kt               ! ocean time-step index 
     189      INTEGER                                  , INTENT( in )  ::  kscheme          ! =0/1/2   type of KEG scheme  
     190      REAL(wp), DIMENSION(jpi,jpj,jpk)         , INTENT(inout) ::  puu, pvv         ! ocean velocities at Kmm 
     191      REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) ::  pu_rhs, pv_rhs   ! RHS  
     192      ! 
     193      INTEGER  ::   ji, jj, jk             ! dummy loop indices 
     194      REAL(wp) ::   zu, zv                   ! local scalars 
     195      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zhke 
     196      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv  
     197      !!---------------------------------------------------------------------- 
     198      ! 
     199      IF( ln_timing )   CALL timing_start('dyn_keg') 
     200      ! 
     201      IF( kt == nit000 ) THEN 
     202         IF(lwp) WRITE(numout,*) 
     203         IF(lwp) WRITE(numout,*) 'dyn_kegAD : kinetic energy gradient trend, scheme number=', kscheme 
     204         IF(lwp) WRITE(numout,*) '~~~~~~~~~' 
     205      ENDIF 
     206       
     207      zhke(:,:,jpk) = 0._wp 
     208 
     209      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==! 
     210!!an45 to be ADDED : que cas C2 - "wet points only" il suffit de x2 le terme quadratic a la coast (nn_dynkeg_adv = 2) 
     211      CASE ( nkeg_C2_wpo )                          !--  Standard scheme  --! 
     212         DO_3D_01_01( 1, jpkm1 ) 
     213            zu =  (   puu(ji-1,jj  ,jk) * puu(ji-1,jj  ,jk)   & 
     214               &    + puu(ji  ,jj  ,jk) * puu(ji  ,jj  ,jk)   ) * ( 2._wp - umask(ji-1,jj,jk) * umask(ji,jj,jk) ) 
     215            zv =  (   pvv(ji  ,jj-1,jk) * pvv(ji  ,jj-1,jk)   & 
     216               &    + pvv(ji  ,jj  ,jk) * pvv(ji  ,jj  ,jk)   ) * ( 2._wp - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) ) 
     217            zhke(ji,jj,jk) = 0.25_wp * ( zv + zu ) 
     218         END_3D 
     219!!an45          
     220      ! 
     221      CASE ( nkeg_C2 )                          !--  Standard scheme  --! 
     222         DO_3D_01_01( 1, jpkm1 ) 
     223            zu =    puu(ji-1,jj  ,jk) * puu(ji-1,jj  ,jk)   & 
     224               &  + puu(ji  ,jj  ,jk) * puu(ji  ,jj  ,jk) 
     225            zv =    pvv(ji  ,jj-1,jk) * pvv(ji  ,jj-1,jk)   & 
     226               &  + pvv(ji  ,jj  ,jk) * pvv(ji  ,jj  ,jk) 
     227            zhke(ji,jj,jk) = 0.25_wp * ( zv + zu ) 
     228         END_3D 
     229!!an 00_00 ? 
     230      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --! 
     231         DO_3D_00_00( 1, jpkm1 ) 
     232            zu = 8._wp * ( puu(ji-1,jj  ,jk) * puu(ji-1,jj  ,jk)    & 
     233               &         + puu(ji  ,jj  ,jk) * puu(ji  ,jj  ,jk) )  & 
     234               &   +     ( puu(ji-1,jj-1,jk) + puu(ji-1,jj+1,jk) ) * ( puu(ji-1,jj-1,jk) + puu(ji-1,jj+1,jk) )   & 
     235               &   +     ( puu(ji  ,jj-1,jk) + puu(ji  ,jj+1,jk) ) * ( puu(ji  ,jj-1,jk) + puu(ji  ,jj+1,jk) ) 
     236               ! 
     237            zv = 8._wp * ( pvv(ji  ,jj-1,jk) * pvv(ji  ,jj-1,jk)    & 
     238               &         + pvv(ji  ,jj  ,jk) * pvv(ji  ,jj  ,jk) )  & 
     239               &  +      ( pvv(ji-1,jj-1,jk) + pvv(ji+1,jj-1,jk) ) * ( pvv(ji-1,jj-1,jk) + pvv(ji+1,jj-1,jk) )   & 
     240               &  +      ( pvv(ji-1,jj  ,jk) + pvv(ji+1,jj  ,jk) ) * ( pvv(ji-1,jj  ,jk) + pvv(ji+1,jj  ,jk) ) 
     241            zhke(ji,jj,jk) = r1_48 * ( zv + zu ) 
     242         END_3D 
     243         CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. ) 
     244         ! 
     245      END SELECT  
     246      ! 
     247         IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN     !***  NO alternating direction  ***! 
     248            ! 
     249            DO_3D_00_00( 1, jpkm1 ) 
     250               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) 
     251               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj) 
     252            END_3D 
     253            ! 
     254         ELSEIF(       PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN            !***  Alternating direction : i-component  ***! 
     255            ! 
     256            DO_3D_00_00( 1, jpkm1 ) 
     257               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) 
     258            END_3D 
     259            ! 
     260         ELSEIF( .NOT. PRESENT( pu_rhs ) .AND.       PRESENT( pv_rhs ) ) THEN            !***  Alternating direction : j-component  ***! 
     261            ! 
     262            DO_3D_00_00( 1, jpkm1 ) 
     263               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj) 
     264            END_3D 
     265            ! 
     266         ENDIF 
     267      IF( ln_timing )   CALL timing_stop('dyn_kegAD') 
     268      ! 
     269   END SUBROUTINE dyn_kegAD 
    158270   !!====================================================================== 
    159271END MODULE dynkeg 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynldf_lap_blp.F90

    r12822 r13005  
    2020   USE in_out_manager ! I/O manager 
    2121   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    22  
     22   USE lib_mpp        ! MPP library 
     23   ! 
     24   USE usrdef_nam , ONLY : nn_dynldf_lap_typ    ! use laplacian parameter 
     25   ! 
    2326   IMPLICIT NONE 
    2427   PRIVATE 
     
    3134   INTEGER, PUBLIC, PARAMETER ::   np_dynldf_lap_symN = 3         ! symmetric laplacian (cartesian) 
    3235    
    33    INTEGER, PUBLIC, PARAMETER ::   ln_dynldf_lap_typ = 1         ! choose type of laplacian (ideally from namelist) 
     36   !INTEGER, PUBLIC, PARAMETER ::   nn_dynldf_lap_typ = 1         ! choose type of laplacian (ideally from namelist) 
    3437!!anSYM 
    3538   !! * Substitutions 
     
    7780         WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator, pass=', kpass 
    7881         WRITE(numout,*) '~~~~~~~ ' 
    79          WRITE(numout,*) '                  ln_dynldf_lap_typ = ', ln_dynldf_lap_typ 
    80          SELECT CASE( ln_dynldf_lap_typ )             ! print the choice of operator 
     82         WRITE(numout,*) '                  nn_dynldf_lap_typ = ', nn_dynldf_lap_typ 
     83         SELECT CASE( nn_dynldf_lap_typ )             ! print the choice of operator 
    8184         CASE( np_dynldf_lap_rot )   ;   WRITE(numout,*) '   ==>>>   div-rot   laplacian' 
    8285         CASE( np_dynldf_lap_sym )   ;   WRITE(numout,*) '   ==>>>   symmetric laplacian (covariant form)' 
    83          CASE( np_dynldf_lap_symN)   ;   WRITE(numout,*) '   ==>>>   symmetric laplacian (simple form)' 
     86         CASE( np_dynldf_lap_symN)   ;   WRITE(numout,*) '   ==>>>   symmetric laplacian (cartesian form)' 
    8487         END SELECT 
    8588      ENDIF 
     
    8992      ENDIF 
    9093      ! 
    91       SELECT CASE( ln_dynldf_lap_typ )   
     94      SELECT CASE( nn_dynldf_lap_typ )   
    9295         !               
    9396         CASE ( np_dynldf_lap_rot )       !==  Vorticity-Divergence form  ==! 
     
    157160            END DO                                           !   End of slab 
    158161            ! 
    159          CASE ( np_dynldf_lap_symN )       !==  Symmetric form  ==!   (naive way) 
     162         CASE ( np_dynldf_lap_symN )       !==  Symmetric form  ==!   (cartesian way) 
    160163            ! 
    161164            DO jk = 1, jpkm1                                 ! Horizontal slab 
     
    190193            !   
    191194         CASE DEFAULT                                     ! error 
    192             CALL ctl_stop('STOP','dyn_ldf_lap: wrong value for ln_dynldf_lap_typ'  ) 
     195            CALL ctl_stop('STOP','dyn_ldf_lap: wrong value for nn_dynldf_lap_typ'  ) 
    193196         END SELECT 
    194197         ! 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynvor.F90

    r12667 r13005  
    2222   !!             -   ! 2018-04  (G. Madec)  add pre-computed gradient for metric term calculation 
    2323   !!            4.x  ! 2020-03  (G. Madec, A. Nasser)  make ln_dynvor_msk truly efficient on relative vorticity 
     24   !!            4.x  ! 2020-03  (G. Madec, A. Nasser)  alternate direction computation of vorticity tendancy 
     25   !!                 !                                 for ENS, ENE 
    2426   !!---------------------------------------------------------------------- 
    2527 
     
    4547   USE lib_mpp        ! MPP library 
    4648   USE timing         ! Timing 
     49!!anAD only 
     50   USE dynkeg, ONLY : dyn_kegAD 
     51   USE dynadv, ONLY : nn_dynkeg 
    4752 
    4853   IMPLICIT NONE 
     
    5358 
    5459   !                                   !!* Namelist namdyn_vor: vorticity term 
    55    LOGICAL, PUBLIC ::   ln_dynvor_ens   !: enstrophy conserving scheme          (ENS) 
    56    LOGICAL, PUBLIC ::   ln_dynvor_ene   !: f-point energy conserving scheme     (ENE) 
    57    LOGICAL, PUBLIC ::   ln_dynvor_enT   !: t-point energy conserving scheme     (ENT) 
    58    LOGICAL, PUBLIC ::   ln_dynvor_eeT   !: t-point energy conserving scheme     (EET) 
    59    LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy & enstrophy conserving scheme (EEN) 
    60    INTEGER, PUBLIC ::      nn_een_e3f      !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 
    61    LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme                         (MIX) 
    62    LOGICAL, PUBLIC ::   ln_dynvor_msk   !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) 
     60   LOGICAL, PUBLIC ::   ln_dynvor_ens      !: enstrophy conserving scheme          (ENS) 
     61   LOGICAL, PUBLIC ::   ln_dynvor_ens_adVO   = .FALSE.   !: AD enstrophy conserving scheme       (ENS_ad) 
     62   LOGICAL, PUBLIC ::   ln_dynvor_ens_adKE   = .FALSE.   !: AD enstrophy conserving scheme       (ENS_ad) 
     63   LOGICAL, PUBLIC ::   ln_dynvor_ens_adKEVO = .FALSE.   !: AD enstrophy conserving scheme       (ENS_ad) 
     64   LOGICAL, PUBLIC ::   ln_dynvor_ene      !: f-point energy conserving scheme     (ENE) 
     65   LOGICAL, PUBLIC ::   ln_dynvor_ene_adVO   = .FALSE.   !: f-point AD energy conserving scheme  (ENE_ad) 
     66   LOGICAL, PUBLIC ::   ln_dynvor_ene_adKE   = .FALSE.   !: f-point AD energy conserving scheme  (ENE_ad) 
     67   LOGICAL, PUBLIC ::   ln_dynvor_ene_adKEVO = .FALSE.   !: f-point AD energy conserving scheme  (ENE_ad) 
     68   LOGICAL, PUBLIC ::   ln_dynvor_enT      !: t-point energy conserving scheme     (ENT) 
     69   LOGICAL, PUBLIC ::   ln_dynvor_eeT      !: t-point energy conserving scheme     (EET) 
     70   LOGICAL, PUBLIC ::   ln_dynvor_een      !: energy & enstrophy conserving scheme (EEN) 
     71   INTEGER, PUBLIC ::      nn_een_e3f         !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 
     72   LOGICAL, PUBLIC ::   ln_dynvor_mix      !: mixed scheme                         (MIX) 
     73   LOGICAL, PUBLIC ::   ln_dynvor_msk      !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) 
    6374 
    6475   INTEGER, PUBLIC ::   nvor_scheme     !: choice of the type of advection scheme 
     
    7081   INTEGER, PUBLIC, PARAMETER ::   np_EEN = 4   ! EEN scheme 
    7182   INTEGER, PUBLIC, PARAMETER ::   np_MIX = 5   ! MIX scheme 
    72  
    73    INTEGER ::   ncor, nrvm, ntot   ! choice of calculated vorticity  
     83!!an 
     84   INTEGER, PUBLIC, PARAMETER ::   np_ENS_adKE   = 11   ! ENS scheme - AD scheme (KE only) 
     85   INTEGER, PUBLIC, PARAMETER ::   np_ENS_adVO   = 12   ! ENS scheme - AD scheme (VOR only) 
     86   INTEGER, PUBLIC, PARAMETER ::   np_ENS_adKEVO = 13   ! ENS scheme - AD scheme (KE+VOR) 
     87   INTEGER, PUBLIC, PARAMETER ::   np_ENE_adKE   = 21   ! ENE scheme - AD scheme (KE only) 
     88   INTEGER, PUBLIC, PARAMETER ::   np_ENE_adVO   = 22   ! ENE scheme - AD scheme (VOR only) 
     89   INTEGER, PUBLIC, PARAMETER ::   np_ENE_adKEVO = 23   ! ENE scheme - AD scheme (KE+VOR) 
     90!!an       
     91    
     92!!an ds step on pourra spécifier la valeur de ntot = np_COR ou np_COR + np_RVO 
     93   INTEGER, PUBLIC ::   ncor, nrvm, ntot   ! choice of calculated vorticity  
    7494   !                               ! associated indices: 
    7595   INTEGER, PUBLIC, PARAMETER ::   np_COR = 1         ! Coriolis (planetary) 
     
    97117CONTAINS 
    98118 
    99    SUBROUTINE dyn_vor( kt, Kmm, puu, pvv, Krhs ) 
     119   SUBROUTINE dyn_vor( kt, Kbb, Kmm, puu, pvv, Krhs) 
    100120      !!---------------------------------------------------------------------- 
    101121      !! 
     
    107127      !!               for futher diagnostics (l_trddyn=T) 
    108128      !!---------------------------------------------------------------------- 
    109       INTEGER                             , INTENT( in  ) ::   kt          ! ocean time-step index 
    110       INTEGER                             , INTENT( in  ) ::   Kmm, Krhs   ! ocean time level indices 
    111       REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv    ! ocean velocity field and RHS of momentum equation 
     129      INTEGER ::   ji, jj, jk   ! dummy loop indice 
     130      INTEGER                             , INTENT( in  ) ::   kt               ! ocean time-step index 
     131      INTEGER                             , INTENT( in  ) ::   Kmm, Krhs, Kbb   ! ocean time level indices 
     132      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv         ! ocean velocity field and RHS of momentum equation 
    112133      ! 
    113134      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdu, ztrdv 
     135      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  zuu, zvv 
    114136      !!---------------------------------------------------------------------- 
    115137      ! 
     
    165187            IF( ln_stcor )   CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend 
    166188         CASE( np_ENE )                        !* energy conserving scheme 
     189                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )   
     190                             CALL vor_ene(   kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend 
     191            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend 
     192             
     193         CASE( np_ENE_adVO )                     !* energy conserving scheme with alternating direction 
     194            IF( MOD( kt , 2 ) ==  1 ) THEN           ! even time step:  u-vor then v-vor components 
     195                               
     196                              !==  Alternative direction - VOR only  ==! 
     197 
     198                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     199 
     200                             ALLOCATE( zuu(jpi,jpj,jpk) ) 
     201                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     202                             zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 
     203                             CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 
     204                             CALL vor_ene( kt, Kmm, ntot, zuu(:,:,:)     , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend 
     205                             DEALLOCATE( zuu ) 
     206            ELSE 
     207                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     208 
     209                             ALLOCATE( zvv(jpi,jpj,jpk) )  
     210                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend 
     211                             zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 
     212                             CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 
     213                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     214                             DEALLOCATE( zvv ) 
     215            ENDIF 
     216            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend 
     217         CASE( np_ENE_adKE )                     !* energy conserving scheme with alternating direction 
     218            IF( MOD( kt , 2 ) ==  1 ) THEN           ! even time step:  u-vor then v-vor components 
     219                                  
     220                                 !==  Alternative direction - KEG only  ==! 
     221 
    167222                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend 
     223 
     224                             ALLOCATE( zuu(jpi,jpj,jpk) ) 
     225                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     226                             zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 
     227                             CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 
     228                             CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:)     , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend 
     229                             DEALLOCATE( zuu ) 
     230            ELSE 
     231 
     232                             CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! total vorticity trend                       
     233 
     234                             ALLOCATE( zvv(jpi,jpj,jpk) ) 
     235                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend 
     236                             zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 
     237                             CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 
     238                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm)  , zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     239                             DEALLOCATE( zvv ) 
     240            ENDIF 
    168241            IF( ln_stcor )   CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )   ! add the Stokes-Coriolis trend 
     242             
     243         CASE( np_ENE_adKEVO )                     !* energy conserving scheme with alternating direction 
     244            IF( MOD( kt , 2 ) ==  1 ) THEN           ! even time step:  u-vor then v-vor components 
     245                               
     246                              !==  Alternative direction - KE + VOR  ==! 
     247 
     248                             ALLOCATE( zuu(jpi,jpj,jpk) ) 
     249                             CALL vor_ene(   kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     250                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   !  
     251                             zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 
     252                             CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 
     253                             CALL vor_ene(   kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend 
     254                             CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )    
     255                             DEALLOCATE( zuu ) 
     256            ELSE 
     257 
     258                             ALLOCATE( zvv(jpi,jpj,jpk) )  
     259                             CALL vor_ene(   kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend 
     260                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )    
     261                             zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 
     262                             CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 
     263                             CALL vor_ene(   kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     264                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )    
     265                             DEALLOCATE( zvv ) 
     266            ENDIF    
    169267         CASE( np_ENS )                        !* enstrophy conserving scheme 
     268                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )   
     269                             CALL vor_ens(   kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! total vorticity trend 
     270            IF( ln_stcor )   CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend 
     271         CASE( np_ENS_adVO )                     !* enstrophy conserving scheme with alternating direction 
     272            IF( MOD( kt , 2 ) ==  1 ) THEN           ! even time step:  u-vor then v-vor components 
     273    
     274                                 !==  Alternative direction - VOR only  ==! 
     275 
     276                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )   
     277 
     278                             ALLOCATE( zuu(jpi,jpj,jpk) ) 
     279                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     280                             zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 
     281                             CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 
     282                             CALL vor_ens( kt, Kmm, ntot, zuu(:,:,:)     , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend 
     283                             DEALLOCATE( zuu ) 
     284            ELSE 
     285                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) 
     286 
     287                             ALLOCATE( zvv(jpi,jpj,jpk) ) 
     288                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend 
     289                             zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 
     290                             CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 
     291                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , zvv(:,:,:), pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     292                             DEALLOCATE( zvv ) 
     293            ENDIF 
     294            IF( ln_stcor )   CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend 
     295         CASE( np_ENS_adKE )                     !* enstrophy conserving scheme with alternating direction 
     296            IF( MOD( kt , 2 ) ==  1 ) THEN           ! even time step:  u-vor then v-vor components 
     297             
     298                                !==  Alternative direction - KEG only  ==! 
     299 
    170300                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! total vorticity trend 
     301 
     302                             ALLOCATE( zuu(jpi,jpj,jpk) ) 
     303                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     304                             zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 
     305                             CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 
     306                             CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:)     , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend 
     307                             DEALLOCATE( zuu ) 
     308            ELSE 
     309 
     310                             CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! total vorticity trend                                                           
     311 
     312                             ALLOCATE( zvv(jpi,jpj,jpk) ) 
     313                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend 
     314                             zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 
     315                             CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 
     316                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , zvv(:,:,:)  , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     317                             DEALLOCATE( zvv ) 
     318            ENDIF 
     319         CASE( np_ENS_adKEVO )                     !* enstrophy conserving scheme with alternating direction 
     320            IF( MOD( kt , 2 ) ==  1 ) THEN           ! even time step:  u-vor then v-vor components 
     321                               
     322                              !==  Alternative direction - KE + VOR  ==! 
     323 
     324                             ALLOCATE( zuu(jpi,jpj,jpk) ) 
     325                             CALL vor_ens(   kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     326                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) )   !  
     327                             zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 
     328                             CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 
     329                             CALL vor_ens(   kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend 
     330                             CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )    
     331                             DEALLOCATE( zuu ) 
     332            ELSE 
     333 
     334                             ALLOCATE( zvv(jpi,jpj,jpk) )  
     335                             CALL vor_ens(   kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )   ! compute and add vv-vorticity trend 
     336                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )    
     337                             zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 
     338                             CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 
     339                             CALL vor_ens(   kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )   ! compute and add uu-vorticity trend 
     340                             CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )    
     341                             DEALLOCATE( zvv ) 
     342            ENDIF    
    171343            IF( ln_stcor )   CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add the Stokes-Coriolis trend 
    172344         CASE( np_MIX )                        !* mixed ene-ens scheme 
     
    313485      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 
    314486      !!---------------------------------------------------------------------- 
    315       INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index 
    316       INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index 
    317       INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric 
    318       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities 
    319       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend 
     487      INTEGER                                  , INTENT(in   ) ::   kt          ! ocean time-step index 
     488      INTEGER                                  , INTENT(in   ) ::   Kmm              ! ocean time level index 
     489      INTEGER                                  , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric 
     490      REAL(wp), DIMENSION(jpi,jpj,jpk)         , INTENT(inout) ::   pu, pv    ! now velocities 
     491      REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend 
    320492      ! 
    321493      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
     
    371543         END SELECT 
    372544         ! 
    373          !                                   !==  horizontal fluxes and potential vorticity ==! 
    374          zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
    375          zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
    376          zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
    377          ! 
    378          !                                   !==  compute and add the vorticity term trend  =! 
    379          DO_2D_00_00 
    380             zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 
    381             zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  ) 
    382             zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 
    383             zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1) 
    384             pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
    385             pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )  
    386          END_2D 
     545         IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN     !***  NO alternating direction  ***! 
     546            ! 
     547            !                                   !==  horizontal fluxes and potential vorticity ==! 
     548            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
     549            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
     550            zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
     551            ! 
     552            !                                   !==  compute and add the vorticity term trend  =! 
     553            DO_2D_00_00 
     554               zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 
     555               zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  ) 
     556               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 
     557               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1) 
     558               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
     559               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )  
     560            END_2D 
     561            !             
     562            ! 
     563         ELSEIF(       PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN            !***  Alternating direction : i-component  ***! 
     564            ! 
     565            ! 
     566            !                                   !==  horizontal fluxes and potential vorticity ==! 
     567            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
     568            zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
     569            ! 
     570            !                                   !==  compute and add the vorticity term trend  =! 
     571            DO_2D_00_00 
     572               zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 
     573               zy2 = zwy(ji,jj  ) + zwy(ji+1,jj  ) 
     574               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
     575            END_2D 
     576            ! 
     577         ELSEIF( .NOT. PRESENT( pu_rhs ) .AND.       PRESENT( pv_rhs ) ) THEN            !***  Alternating direction : j-component  ***! 
     578            ! 
     579            ! 
     580            !                                   !==  horizontal fluxes and potential vorticity ==! 
     581            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
     582            zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
     583            ! 
     584            !                                   !==  compute and add the vorticity term trend  =! 
     585            DO_2D_00_00 
     586               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 
     587               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1) 
     588               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )  
     589            END_2D 
     590            ! 
     591         ENDIF 
    387592         !                                             ! =============== 
    388593      END DO                                           !   End of slab 
     
    411616      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 
    412617      !!---------------------------------------------------------------------- 
    413       INTEGER                         , INTENT(in   ) ::   kt          ! ocean time-step index 
    414       INTEGER                         , INTENT(in   ) ::   Kmm              ! ocean time level index 
    415       INTEGER                         , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric 
    416       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu, pv    ! now velocities 
    417       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend 
     618      INTEGER                                  , INTENT(in   ) ::   kt          ! ocean time-step index 
     619      INTEGER                                  , INTENT(in   ) ::   Kmm              ! ocean time level index 
     620      INTEGER                                  , INTENT(in   ) ::   kvor        ! total, planetary, relative, or metric 
     621      REAL(wp), DIMENSION(jpi,jpj,jpk)         , INTENT(inout) ::   pu, pv    ! now velocities 
     622      REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) ::   pu_rhs, pv_rhs    ! total v-trend 
    418623      ! 
    419624      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     
    469674         ! 
    470675         ! 
    471          !                                   !==  horizontal fluxes and potential vorticity ==! 
    472          zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
    473          zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
    474          zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
    475          ! 
    476          !                                   !==  compute and add the vorticity term trend  =! 
    477          DO_2D_00_00 
    478             zuav = r1_8 * r1_e1u(ji,jj) * (  zwy(ji  ,jj-1) + zwy(ji+1,jj-1)  & 
    479                &                           + zwy(ji  ,jj  ) + zwy(ji+1,jj  )  ) 
    480             zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  & 
    481                &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1)  ) 
    482             pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) ) 
    483             pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) ) 
    484          END_2D 
     676!!an wut ? v et u  
     677         IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN     !***  NO alternating direction  ***! 
     678            ! 
     679            !                                   !==  horizontal fluxes and potential vorticity ==! 
     680            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
     681            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
     682            zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
     683            ! 
     684            !                                   !==  compute and add the vorticity term trend  =! 
     685            DO_2D_00_00 
     686               zuav = r1_8 * r1_e1u(ji,jj) * (  zwy(ji  ,jj-1) + zwy(ji+1,jj-1)  & 
     687                  &                           + zwy(ji  ,jj  ) + zwy(ji+1,jj  )  ) 
     688               zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  & 
     689                  &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1)  ) 
     690               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) ) 
     691               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) ) 
     692            END_2D 
     693            ! 
     694         ELSEIF(       PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN            !***  Alternating direction : i-component  ***! 
     695            ! 
     696            ! 
     697            !                                   !==  horizontal fluxes and potential vorticity ==! 
     698            zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 
     699            zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
     700            ! 
     701            !                                   !==  compute and add the vorticity term trend  =! 
     702            DO_2D_00_00 
     703               zuav = r1_8 * r1_e1u(ji,jj) * (  zwy(ji  ,jj-1) + zwy(ji+1,jj-1)  & 
     704                  &                           + zwy(ji  ,jj  ) + zwy(ji+1,jj  )  ) 
     705               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) ) 
     706            END_2D 
     707            ! 
     708         ELSEIF( .NOT. PRESENT( pu_rhs ) .AND.       PRESENT( pv_rhs ) ) THEN            !***  Alternating direction : j-component  ***! 
     709            ! 
     710            ! 
     711            !                                   !==  horizontal fluxes and potential vorticity ==! 
     712            zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 
     713            zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 
     714            ! 
     715            !                                   !==  compute and add the vorticity term trend  =! 
     716            DO_2D_00_00 
     717               zvau =-r1_8 * r1_e2v(ji,jj) * (  zwx(ji-1,jj  ) + zwx(ji-1,jj+1)  & 
     718                  &                           + zwx(ji  ,jj  ) + zwx(ji  ,jj+1)  ) 
     719               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) ) 
     720            END_2D 
     721            ! 
     722         ENDIF 
    485723         !                                             ! =============== 
    486724      END DO                                           !   End of slab 
     
    7621000      !! 
    7631001      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT,   & 
    764          &                 ln_dynvor_een, nn_een_e3f   , ln_dynvor_mix, ln_dynvor_msk 
     1002         &                 ln_dynvor_een, nn_een_e3f   , ln_dynvor_mix, ln_dynvor_msk,   & 
     1003         &                 ln_dynvor_ens_adVO, ln_dynvor_ens_adKE, ln_dynvor_ens_adKEVO, &   ! Alternative direction parameters 
     1004         &                 ln_dynvor_ene_adVO, ln_dynvor_ene_adKE, ln_dynvor_ene_adKEVO 
    7651005      !!---------------------------------------------------------------------- 
    7661006      ! 
     
    7791019      IF(lwp) THEN                    ! Namelist print 
    7801020         WRITE(numout,*) '   Namelist namdyn_vor : choice of the vorticity term scheme' 
    781          WRITE(numout,*) '      enstrophy conserving scheme                    ln_dynvor_ens = ', ln_dynvor_ens 
    782          WRITE(numout,*) '      f-point energy conserving scheme               ln_dynvor_ene = ', ln_dynvor_ene 
    783          WRITE(numout,*) '      t-point energy conserving scheme               ln_dynvor_enT = ', ln_dynvor_enT 
    784          WRITE(numout,*) '      energy conserving scheme  (een using e3t)      ln_dynvor_eeT = ', ln_dynvor_eeT 
    785          WRITE(numout,*) '      enstrophy and energy conserving scheme         ln_dynvor_een = ', ln_dynvor_een 
    786          WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f = ', nn_een_e3f 
    787          WRITE(numout,*) '      mixed enstrophy/energy conserving scheme       ln_dynvor_mix = ', ln_dynvor_mix 
    788          WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk = ', ln_dynvor_msk 
     1021         WRITE(numout,*) '      enstrophy conserving scheme                    ln_dynvor_ens    = ', ln_dynvor_ens 
     1022         WRITE(numout,*) '      f-point energy conserving scheme               ln_dynvor_ene    = ', ln_dynvor_ene 
     1023         WRITE(numout,*) '      t-point energy conserving scheme               ln_dynvor_enT    = ', ln_dynvor_enT 
     1024         WRITE(numout,*) '      energy conserving scheme  (een using e3t)      ln_dynvor_eeT    = ', ln_dynvor_eeT 
     1025         WRITE(numout,*) '      enstrophy and energy conserving scheme         ln_dynvor_een    = ', ln_dynvor_een 
     1026         WRITE(numout,*) '         e3f = averaging /4 (=0) or /sum(tmask) (=1)    nn_een_e3f    = ', nn_een_e3f 
     1027         WRITE(numout,*) '      mixed enstrophy/energy conserving scheme       ln_dynvor_mix    = ', ln_dynvor_mix 
     1028         WRITE(numout,*) '      masked (=T) or unmasked(=F) vorticity          ln_dynvor_msk    = ', ln_dynvor_msk 
    7891029      ENDIF 
    7901030 
     
    7991039         DO_3D_10_10( 1, jpk ) 
    8001040            IF(    tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)              & 
    801                & + tmask(ji,jj  ,jk) + tmask(ji+1,jj+1,jk) == 3._wp )   fmask(ji,jj,jk) = 1._wp 
     1041               & + tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk) == 3._wp )   fmask(ji,jj,jk) = 1._wp 
    8021042         END_3D 
    8031043         ! 
     
    8091049      ioptio = 0                     ! type of scheme for vorticity (set nvor_scheme) 
    8101050      IF( ln_dynvor_ens ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENS   ;   ENDIF 
     1051      IF( ln_dynvor_ens_adVO ) THEN   ;   ioptio = ioptio + 1     ;   nvor_scheme = np_ENS_adVO   ;   ENDIF 
     1052      IF( ln_dynvor_ens_adKE ) THEN   ;   ioptio = ioptio + 1     ;   nvor_scheme = np_ENS_adKE   ;   ENDIF 
     1053      IF( ln_dynvor_ens_adKEVO ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENS_adKEVO   ;   ENDIF 
    8111054      IF( ln_dynvor_ene ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENE   ;   ENDIF 
     1055      IF( ln_dynvor_ene_adVO ) THEN   ;   ioptio = ioptio + 1     ;   nvor_scheme = np_ENE_adVO   ;   ENDIF 
     1056      IF( ln_dynvor_ene_adKE ) THEN   ;   ioptio = ioptio + 1     ;   nvor_scheme = np_ENE_adKE   ;   ENDIF 
     1057      IF( ln_dynvor_ene_adKEVO ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENE_adKEVO   ;   ENDIF 
    8121058      IF( ln_dynvor_enT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_ENT   ;   ENDIF 
    8131059      IF( ln_dynvor_eeT ) THEN   ;   ioptio = ioptio + 1   ;   nvor_scheme = np_EET   ;   ENDIF 
     
    8261072      CASE( np_VEC_c2  ) 
    8271073         IF(lwp) WRITE(numout,*) '   ==>>>   vector form dynamics : total vorticity = Coriolis + relative vorticity'  
    828          nrvm = np_RVO        ! relative vorticity 
     1074         nrvm = np_RVO        ! relative vorticity       
    8291075         ntot = np_CRV        ! relative + planetary vorticity          
    8301076      CASE( np_FLX_c2 , np_FLX_ubs  ) 
     
    8561102         WRITE(numout,*) 
    8571103         SELECT CASE( nvor_scheme ) 
    858          CASE( np_ENS )   ;   WRITE(numout,*) '   ==>>>   enstrophy conserving scheme (ENS)' 
    859          CASE( np_ENE )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at F-points) (ENE)' 
    860          CASE( np_ENT )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at T-points) (ENT)' 
    861          CASE( np_EET )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3t) (EET)' 
    862          CASE( np_EEN )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)' 
    863          CASE( np_MIX )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)' 
     1104          
     1105         CASE( np_ENS    )   ;   WRITE(numout,*) '   ==>>>   enstrophy conserving scheme (ENS)' 
     1106         CASE( np_ENS_adVO ) ;   WRITE(numout,*) '   ==>>>   AD enstrophy conserving scheme (ENS_adVO) on vorticity only' 
     1107         CASE( np_ENS_adKE ) ;   WRITE(numout,*) '   ==>>>   AD enstrophy conserving scheme (ENS_adKE) on kinetic energy only' 
     1108         CASE( np_ENS_adKEVO ) ;   WRITE(numout,*) '   ==>>>   AD enstrophy conserving scheme (ENS_adKEVO) on kinetic energy and vorticity' 
     1109          
     1110         CASE( np_ENE    )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at F-points) (ENE)' 
     1111         CASE( np_ENE_adVO ) ;   WRITE(numout,*) '   ==>>>   AD energy conserving scheme (Coriolis at F-points) (ENE_adVO) on vorticity only' 
     1112         CASE( np_ENE_adKE ) ;   WRITE(numout,*) '   ==>>>   AD energy conserving scheme (Coriolis at F-points) (ENE_adKE) on kinetic energy only' 
     1113         CASE( np_ENE_adKEVO ) ;   WRITE(numout,*) '   ==>>>   AD energy conserving scheme (Coriolis at F-points) (ENE_adKEVO) on kinetic energy and vorticity' 
     1114          
     1115         CASE( np_ENT    )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (Coriolis at T-points) (ENT)' 
     1116         CASE( np_EET    )   ;   WRITE(numout,*) '   ==>>>   energy conserving scheme (EEN scheme using e3t) (EET)' 
     1117         CASE( np_EEN    )   ;   WRITE(numout,*) '   ==>>>   energy and enstrophy conserving scheme (EEN)' 
     1118         CASE( np_MIX    )   ;   WRITE(numout,*) '   ==>>>   mixed enstrophy/energy conserving scheme (MIX)' 
    8641119         END SELECT          
    8651120      ENDIF 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/ldfdyn.F90

    r12822 r13005  
    2525   USE lib_mpp         ! distribued memory computing library 
    2626   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    27  
     27   ! 
     28   USE usrdef_nam , ONLY : ln_dynldf_lap_PM 
     29   ! 
    2830   IMPLICIT NONE 
    2931   PRIVATE 
     
    6062   INTEGER           , PUBLIC ::   nldf_dyn         !: type of lateral diffusion used defined from ln_dynldf_... (namlist logicals) 
    6163   LOGICAL           , PUBLIC ::   l_ldfdyn_time    !: flag for time variation of the lateral eddy viscosity coef. 
    62  
     64!!an 
     65   !LOGICAL           , PUBLIC ::   ll_dynldf_lap_PM !: flag for P.Marchand modification on viscosity 
     66!!an 
    6367   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ahmt, ahmf   !: eddy viscosity coef. at T- and F-points [m2/s or m4/s] 
    6468   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dtensq       !: horizontal tension squared         (Smagorinsky only) 
     
    323327         IF( .NOT.l_ldfdyn_time ) THEN             !* No time variation  
    324328            IF(     ln_dynldf_lap ) THEN                 !   laplacian operator (mask only) 
    325                ahmt(:,:,1:jpkm1) =       ahmt(:,:,1:jpkm1)   * tmask(:,:,1:jpkm1) 
    326                WRITE(numout,*) ' ahmt tmask ' 
     329!!an          ! 
     330            WRITE(numout,*) '   ln_dynldf_lap_PM = ',ln_dynldf_lap_PM  
     331               IF(     ln_dynldf_lap_PM ) THEN                 !   laplacian operator (mask only) 
    327332!! mask tension at the coast (equivalent of ghostpoints at T) 
    328 !              DO jk = 1, jpk 
    329 !                 DO jj = 1, jpjm1 
    330 !                    DO ji = 1, jpim1      ! NO vector opt. 
    331 !                       ! si sum(fmask)==3 = mouillé (on touche pas) 
    332 !                       ! si sum = 2 alors on met a 0 zsum = fmask + fmask + fmask,.. et si zsum < 2 -> 0 sinon = 1 
    333 !                       zsum =   fmask(ji,jj  ,jk) + fmask(ji+1,jj  ,jk)   & 
    334 !                          &   + fmask(ji,jj+1,jk) + fmask(ji+1,jj+1,jk) 
    335 !                       IF ( zsum < 2._wp )   ahmt(ji,jj,jk) = 0 
    336 !                       ! 
    337 !                       !ahmt(ji,jj,jk) = ahmt(ji,jj,jk) * fmask(ji,jj  ,jk) * fmask(ji+1,jj  ,jk)   & 
    338 !                       !   &                            * fmask(ji,jj+1,jk) * fmask(ji+1,jj+1,jk) 
    339 !                    END DO 
    340 !                 END DO 
    341 !              END DO 
    342 !              ahmt(jpi,:,1:jpkm1) =  0._wp 
    343 !              ahmt(:,jpj,1:jpkm1) =  0._wp 
    344 !              WRITE(numout,*) '   an45 ahmt x0' 
    345  
    346                ahmf(:,:,1:jpkm1) =       ahmf(:,:,1:jpkm1)   * fmask(:,:,1:jpkm1) 
    347                WRITE(numout,*) ' ahmf fmask ' 
    348 !!an apply no slip at the coast (ssfmask = 1 within the domain and at the coast contrary to fmask in free slip) 
    349 !               DO jk = 1, jpkm1 
    350 !                  ahmf(:,:,jk) =    ahmf(:,:,jk) * ( 2._wp * ssfmask(:,:) - fmask(:,:,jk) ) 
    351 !               END DO 
    352 !               WRITE(numout,*) '   an45 ahmf x2' 
    353  
     333                  DO jk = 1, jpk 
     334                     DO jj = 1, jpjm1 
     335                        DO ji = 1, jpim1      ! NO vector opt. 
     336                           ! si sum(fmask)==3 = mouillé (on touche pas) 
     337                           ! si sum = 2 alors on met a 0 zsum = fmask + fmask + fmask,.. et si zsum < 2 -> 0 sinon = 1 
     338                           zsum =   fmask(ji,jj  ,jk) + fmask(ji+1,jj  ,jk)   & 
     339                              &   + fmask(ji,jj+1,jk) + fmask(ji+1,jj+1,jk) 
     340                           IF ( zsum < 2._wp )   ahmt(ji,jj,jk) = 0 
     341                           ! 
     342                           !ahmt(ji,jj,jk) = ahmt(ji,jj,jk) * fmask(ji,jj  ,jk) * fmask(ji+1,jj  ,jk)   & 
     343                           !   &                            * fmask(ji,jj+1,jk) * fmask(ji+1,jj+1,jk) 
     344                        END DO 
     345                     END DO 
     346                  END DO 
     347                  ahmt(jpi,:,1:jpkm1) =  0._wp 
     348                  ahmt(:,jpj,1:jpkm1) =  0._wp 
     349                  WRITE(numout,*) '  ahmt x0' 
     350!! apply no slip at the coast (ssfmask = 1 within the domain and at the coast contrary to fmask in free slip) 
     351                   DO jk = 1, jpkm1 
     352                      ahmf(:,:,jk) =    ahmf(:,:,jk) * ( 2._wp * ssfmask(:,:) - fmask(:,:,jk) ) 
     353                   END DO 
     354                   WRITE(numout,*) '  ahmf x2' 
     355               ELSE 
     356               ! classic boundary condition on the viscosity coefficient 
     357                  ahmt(:,:,1:jpkm1) =       ahmt(:,:,1:jpkm1)   * tmask(:,:,1:jpkm1) 
     358                  WRITE(numout,*) ' ahmt tmasked ' 
     359                  ahmf(:,:,1:jpkm1) =       ahmf(:,:,1:jpkm1)   * fmask(:,:,1:jpkm1) 
     360                  WRITE(numout,*) ' ahmf fmasked ' 
     361               ENDIF 
     362!!an         !                  
    354363            ELSEIF( ln_dynldf_blp ) THEN                 ! bilaplacian operator (square root + mask) 
    355364               ahmt(:,:,1:jpkm1) = SQRT( ahmt(:,:,1:jpkm1) ) * tmask(:,:,1:jpkm1) 
  • NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/step.F90

    r12614 r13005  
    1313   USE phycst           ! physical constants 
    1414   USE usrdef_nam 
     15   USE lib_mpp        ! MPP library 
     16   USE dynvor , ONLY : ln_dynvor_ens_adVO, ln_dynvor_ens_adKE, ln_dynvor_ens_adKEVO, & 
     17    &                  ln_dynvor_ene_adVO, ln_dynvor_ene_adKE, ln_dynvor_ene_adKEVO     
    1518   ! 
    1619   USE iom              ! xIOs server  
     
    138141               &         CALL Agrif_Sponge_dyn        ! momentum sponge 
    139142#endif 
    140                          CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF)   ==> RHS 
    141   
    142                          CALL dyn_vor( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity              ==> RHS 
    143   
    144                          CALL dyn_ldf( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! lateral mixing 
    145143 
    146144!!an - calcul du gradient de pression horizontal (explicit) 
     
    150148      END_3D 
    151149      ! 
     150       
     151!      IF( kstp == nit000 .AND. lwp ) THEN 
     152!         WRITE(numout,*) 
     153!         WRITE(numout,*) 'step.F90 : classic script used' 
     154!         WRITE(numout,*) '~~~~~~~' 
     155!         IF(       ln_dynvor_ens_adVO .OR. ln_dynvor_ens_adKE .OR. ln_dynvor_ens_adKEVO   & 
     156!         &    .OR. ln_dynvor_ene_adVO .OR. ln_dynvor_ene_adKE .OR. ln_dynvor_ene_adKEVO   ) THEN 
     157!            CALL ctl_stop('STOP','step : alternative direction asked but classis step'  ) 
     158!         ENDIF 
     159!      ENDIF 
     160!!an      
     161!                         CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF)  ==> RHS 
     162!  
     163!                         CALL dyn_vor( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! vorticity             ==> RHS 
     164!  
     165!!an     In dynvor, dynkegAD is called even if not AD, so we keep the same step.F90 
     166   
     167                         CALL dyn_vor( kstp, Nbb, Nnn      , uu, vv, Nrhs)  ! vorticity            ==> RHS 
     168   
     169                         CALL dyn_ldf( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! lateral mixing 
     170 
    152171      ! add wind stress forcing and layer linear friction to the RHS  
    153172      z1_2rho0 = 0.5_wp * r1_rho0 
Note: See TracChangeset for help on using the changeset viewer.