New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 5321 – NEMO

Changeset 5321


Ignore:
Timestamp:
2015-06-01T09:24:31+02:00 (9 years ago)
Author:
pabouttier
Message:

Add Hollingsworth correction in the computation of the gradient of Kinetic Energy

Location:
trunk/NEMOGCM
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMOGCM/CONFIG/SHARED/namelist_ref

    r5190 r5321  
    2828   nn_it000    =       1   !  first time step 
    2929   nn_itend    =    5475   !  last  time step (std 5475) 
    30    nn_date0    =  010101   !  date at nit_0000 (format yyyymmdd) used if ln_rstart=F or (ln_rstart=T and nn_rstctl=0 or 1) 
     30   nn_date0    =  20061101   !  date at nit_0000 (format yyyymmdd) used if ln_rstart=F or (ln_rstart=T and nn_rstctl=0 or 1) 
    3131   nn_leapy    =       0   !  Leap year calendar (1) or not (0) 
    3232   ln_rstart   = .false.   !  start from rest (F) or from a restart file (T) 
     
    118118   rn_bathy    =    0.     !  value of the bathymetry. if (=0) bottom flat at jpkm1 
    119119   nn_closea   =    0      !  remove (=0) or keep (=1) closed seas and lakes (ORCA) 
    120    nn_msh      =    0      !  create (=1) a mesh file or not (=0) 
     120   nn_msh      =    1      !  create (=1) a mesh file or not (=0) 
    121121   rn_hmin     =   -3.     !  min depth of the ocean (>0) or min number of ocean level (<0) 
    122122   rn_e3zps_min=   20.     !  partial step thickness is set larger than the minimum of 
     
    803803!----------------------------------------------------------------------- 
    804804   ln_dynadv_vec = .true.  !  vector form (T) or flux form (F) 
     805   nn_dynkeg     = 0       ! scheme for grad(KE): =0   C2  ;  =1   Hollingsworth correction 
    805806   ln_dynadv_cen2= .false. !  flux form - 2nd order centered scheme 
    806807   ln_dynadv_ubs = .false. !  flux form - 3rd order UBS      scheme 
  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90

    r4990 r5321  
    44   !! Ocean dynamics:  kinetic energy gradient trend 
    55   !!====================================================================== 
    6    !! History :  1.0  !  87-09  (P. Andrich, m.-a. Foujols)  Original code 
    7    !!            7.0  !  97-05  (G. Madec)  Split dynber into dynkeg and dynhpg 
    8    !!            9.0  !  02-07  (G. Madec)  F90: Free form and module 
     6   !! History :  1.0  !  1987-09  (P. Andrich, M.-A. Foujols)  Original code 
     7   !!            7.0  !  1997-05  (G. Madec)  Split dynber into dynkeg and dynhpg 
     8   !!  NEMO      1.0  !  2002-07  (G. Madec)  F90: Free form and module 
     9   !!            3.6  !  2015-05  (N. Ducousso, G. Madec)  add Hollingsworth scheme as an option  
    910   !!---------------------------------------------------------------------- 
    1011    
     
    1819   ! 
    1920   USE in_out_manager  ! I/O manager 
     21   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2022   USE lib_mpp         ! MPP library 
    2123   USE prtctl          ! Print control 
     
    2830   PUBLIC   dyn_keg    ! routine called by step module 
    2931    
     32   INTEGER, PARAMETER, PUBLIC  ::   nkeg_C2  = 0   !: 2nd order centered scheme (standard scheme) 
     33   INTEGER, PARAMETER, PUBLIC  ::   nkeg_HW  = 1   !: Hollingsworth et al., QJRMS, 1983 
     34   ! 
     35   REAL(wp) ::   r1_48 = 1._wp / 48._wp   !: =1/(4*2*6) 
     36    
    3037   !! * Substitutions 
    3138#  include "vectopt_loop_substitute.h90" 
    3239   !!---------------------------------------------------------------------- 
    33    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     40   !! NEMO/OPA 3.6 , NEMO Consortium (2015) 
    3441   !! $Id$  
    3542   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    3744CONTAINS 
    3845 
    39    SUBROUTINE dyn_keg( kt ) 
     46   SUBROUTINE dyn_keg( kt, kscheme ) 
    4047      !!---------------------------------------------------------------------- 
    4148      !!                  ***  ROUTINE dyn_keg  *** 
     
    4552      !!      general momentum trend. 
    4653      !! 
    47       !! ** Method  :   Compute the now horizontal kinetic energy  
     54      !! ** Method  : * kscheme = nkeg_C2 : 2nd order centered scheme that  
     55      !!      conserve kinetic energy. Compute the now horizontal kinetic energy  
    4856      !!         zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ] 
     57      !!              * kscheme = nkeg_HW : Hollingsworth correction following 
     58      !!      Arakawa (2001). The now horizontal kinetic energy is given by: 
     59      !!         zhke = 1/6 [ mi-1(  2 * un^2 + ((un(j+1)+un(j-1))/2)^2  ) 
     60      !!                    + mj-1(  2 * vn^2 + ((vn(i+1)+vn(i-1))/2)^2  ) ] 
     61      !!       
    4962      !!      Take its horizontal gradient and add it to the general momentum 
    5063      !!      trend (ua,va). 
     
    5467      !! ** Action : - Update the (ua, va) with the hor. ke gradient trend 
    5568      !!             - send this trends to trd_dyn (l_trddyn=T) for post-processing 
     69      !! 
     70      !! ** References : Arakawa, A., International Geophysics 2001. 
     71      !!                 Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983. 
    5672      !!---------------------------------------------------------------------- 
    57       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     73      INTEGER, INTENT( in ) ::   kt        ! ocean time-step index 
     74      INTEGER, INTENT( in ) ::   kscheme   ! =0/1   type of KEG scheme  
    5875      ! 
    5976      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     
    6380      !!---------------------------------------------------------------------- 
    6481      ! 
    65       IF( nn_timing == 1 )  CALL timing_start('dyn_keg') 
     82      IF( nn_timing == 1 )   CALL timing_start('dyn_keg') 
    6683      ! 
    67       CALL wrk_alloc( jpi, jpj, jpk, zhke ) 
     84      CALL wrk_alloc( jpi,jpj,jpk,  zhke ) 
    6885      ! 
    6986      IF( kt == nit000 ) THEN 
    7087         IF(lwp) WRITE(numout,*) 
    71          IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend' 
     88         IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme 
    7289         IF(lwp) WRITE(numout,*) '~~~~~~~' 
    7390      ENDIF 
    7491 
    7592      IF( l_trddyn ) THEN           ! Save ua and va trends 
    76          CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 
     93         CALL wrk_alloc( jpi,jpj,jpk,   ztrdu, ztrdv ) 
    7794         ztrdu(:,:,:) = ua(:,:,:)  
    7895         ztrdv(:,:,:) = va(:,:,:)  
    7996      ENDIF 
    8097       
    81       !                                                ! =============== 
    82       DO jk = 1, jpkm1                                 ! Horizontal slab 
    83          !                                             ! =============== 
    84          DO jj = 2, jpj         ! Horizontal kinetic energy at T-point 
    85             DO ji = fs_2, jpi   ! vector opt. 
    86                zu = 0.25 * (  un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   & 
    87                   &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)  ) 
    88                zv = 0.25 * (  vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   & 
    89                   &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)  ) 
    90                zhke(ji,jj,jk) = zv + zu 
    91 !!gm simplier coding  ==>>   ~ faster 
    92 !    don't forget to suppress local zu zv scalars 
    93 !               zhke(ji,jj,jk) = 0.25 * (   un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   & 
    94 !                  &                      + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)   & 
    95 !                  &                      + vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   & 
    96 !                  &                      + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)   ) 
    97 !!gm end <<== 
    98             END DO   
    99          END DO   
    100          DO jj = 2, jpjm1       ! add the gradient of kinetic energy to the general momentum trends 
     98      zhke(:,:,jpk) = 0._wp 
     99       
     100      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==! 
     101      ! 
     102      CASE ( nkeg_C2 )                          !--  Standard scheme  --! 
     103         DO jk = 1, jpkm1 
     104            DO jj = 2, jpj 
     105               DO ji = fs_2, jpi   ! vector opt. 
     106                  zu =    un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   & 
     107                     &  + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) 
     108                  zv =    vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   & 
     109                     &  + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) 
     110                  zhke(ji,jj,jk) = 0.25_wp * ( zv + zu ) 
     111               END DO   
     112            END DO 
     113         END DO 
     114         ! 
     115      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --! 
     116         DO jk = 1, jpkm1 
     117            DO jj = 2, jpjm1        
     118               DO ji = fs_2, jpim1   ! vector opt. 
     119                  zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    & 
     120                     &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  & 
     121                     &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   & 
     122                     &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) 
     123                     ! 
     124                  zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    & 
     125                     &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  & 
     126                     &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   & 
     127                     &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) 
     128                  zhke(ji,jj,jk) = r1_48 * ( zv + zu ) 
     129               END DO   
     130            END DO 
     131         END DO 
     132         CALL lbc_lnk( zhke, 'T', 1. ) 
     133         ! 
     134      END SELECT 
     135      ! 
     136      DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==! 
     137         DO jj = 2, jpjm1 
    101138            DO ji = fs_2, fs_jpim1   ! vector opt. 
    102139               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) 
     
    104141            END DO  
    105142         END DO 
    106 !!gm idea to be tested  ==>>   is it faster on scalar computers ? 
    107 !         DO jj = 2, jpjm1       ! add the gradient of kinetic energy to the general momentum trends 
    108 !            DO ji = fs_2, fs_jpim1   ! vector opt. 
    109 !               ua(ji,jj,jk) = ua(ji,jj,jk) - 0.25 * ( + un(ji+1,jj  ,jk) * un(ji+1,jj  ,jk)   & 
    110 !                  &                                   + vn(ji+1,jj-1,jk) * vn(ji+1,jj-1,jk)   & 
    111 !                  &                                   + vn(ji+1,jj  ,jk) * vn(ji+1,jj  ,jk)   & 
    112 !                  ! 
    113 !                  &                                   - un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   & 
    114 !                  &                                   - vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   & 
    115 !                  &                                   - vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)   ) / e1u(ji,jj) 
    116 !                  ! 
    117 !               va(ji,jj,jk) = va(ji,jj,jk) - 0.25 * (   un(ji-1,jj+1,jk) * un(ji-1,jj+1,jk)   & 
    118 !                  &                                   + un(ji  ,jj+1,jk) * un(ji  ,jj+1,jk)   & 
    119 !                  &                                   + vn(ji  ,jj+1,jk) * vn(ji  ,jj+1,jk)   & 
    120 !                  ! 
    121 !                  &                                   - un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   & 
    122 !                  &                                   - un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)   & 
    123 !                  &                                   - vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)   ) / e2v(ji,jj) 
    124 !            END DO  
    125 !         END DO 
    126 !!gm en idea            <<== 
    127          !                                             ! =============== 
    128       END DO                                           !   End of slab 
    129       !                                                ! =============== 
    130  
    131       IF( l_trddyn ) THEN      ! save the Kinetic Energy trends for diagnostic 
     143      END DO 
     144      ! 
     145      IF( l_trddyn ) THEN                 ! save the Kinetic Energy trends for diagnostic 
    132146         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    133147         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    134148         CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt ) 
    135          CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 
     149         CALL wrk_dealloc( jpi,jpj,jpk,   ztrdu, ztrdv ) 
    136150      ENDIF 
    137151      ! 
     
    139153         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    140154      ! 
    141       CALL wrk_dealloc( jpi, jpj, jpk, zhke ) 
     155      CALL wrk_dealloc( jpi,jpj,jpk,  zhke ) 
    142156      ! 
    143       IF( nn_timing == 1 )  CALL timing_stop('dyn_keg') 
     157      IF( nn_timing == 1 )   CALL timing_stop('dyn_keg') 
    144158      ! 
    145159   END SUBROUTINE dyn_keg 
Note: See TracChangeset for help on using the changeset viewer.