source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_lap.F90 @ 8627

Last change on this file since 8627 was 8627, checked in by gm, 3 years ago

#1963 : Correct implementation of CMIP6 KE dissipation + move EKE⇔PE diag in traadv_eiv (in both cases, no more need of ln_trd_KE=T)

  • Property svn:keywords set to Id
File size: 6.6 KB
Line 
1MODULE dynldf_lap
2   !!======================================================================
3   !!                       ***  MODULE  dynldf_lap  ***
4   !! Ocean dynamics:  lateral viscosity trend
5   !!======================================================================
6   !! History :  OPA  ! 1990-09 (G. Madec) Original code
7   !!            4.0  ! 1991-11 (G. Madec)
8   !!            6.0  ! 1996-01 (G. Madec) statement function for e3 and ahm
9   !!   NEMO     1.0  ! 2002-06 (G. Madec)  F90: Free form and module
10   !!             -   ! 2004-08 (C. Talandier) New trends organization
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   dyn_ldf_lap  : update the momentum trend with the lateral diffusion
15   !!                  using an iso-level harmonic operator
16   !!----------------------------------------------------------------------
17   USE oce             ! ocean dynamics and tracers
18   USE dom_oce         ! ocean space and time domain
19   USE phycst          ! physical constants
20   USE ldfdyn_oce      ! ocean dynamics: lateral physics
21   USE zdf_oce         ! ocean vertical physics
22   !
23   USE in_out_manager  ! I/O manager
24   USE iom             ! I/O library
25   USE timing          ! Timing
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC dyn_ldf_lap  ! called by step.F90
31
32   !! * Substitutions
33#  include "domzgr_substitute.h90"
34#  include "ldfdyn_substitute.h90"
35#  include "vectopt_loop_substitute.h90"
36   !!----------------------------------------------------------------------
37   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
38   !! $Id$
39   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
40   !!----------------------------------------------------------------------
41CONTAINS
42
43   SUBROUTINE dyn_ldf_lap( kt )
44      !!----------------------------------------------------------------------
45      !!                     ***  ROUTINE dyn_ldf_lap  ***
46      !!                       
47      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
48      !!      trend and add it to the general trend of tracer equation.
49      !!
50      !! ** Method  :   The before horizontal momentum diffusion trend is an
51      !!      harmonic operator (laplacian type) which separates the divergent
52      !!      and rotational parts of the flow.
53      !!      Its horizontal components are computed as follow:
54      !!         difu = 1/e1u di[ahmt hdivb] - 1/(e2u*e3u) dj-1[e3f ahmf rotb]
55      !!         difv = 1/e2v dj[ahmt hdivb] + 1/(e1v*e3v) di-1[e3f ahmf rotb]
56      !!      in the rotational part of the diffusion.
57      !!      Add this before trend to the general trend (ua,va):
58      !!            (ua,va) = (ua,va) + (diffu,diffv)
59      !!
60      !! ** Action : - Update (ua,va) with the iso-level harmonic mixing trend
61      !!----------------------------------------------------------------------
62      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
63      !
64      INTEGER  ::   ji, jj, jk             ! dummy loop indices
65      REAL(wp) ::   zua, zva, ze2u, ze1v   ! local scalars
66      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z2d   ! 2D workspace
67      !!----------------------------------------------------------------------
68      !
69      IF( nn_timing == 1 )  CALL timing_start('dyn_ldf_lap')
70      !
71      IF( kt == nit000 ) THEN
72         IF(lwp) WRITE(numout,*)
73         IF(lwp) WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator'
74         IF(lwp) WRITE(numout,*) '~~~~~~~ '
75      ENDIF
76      !                                                ! ===============
77      DO jk = 1, jpkm1                                 ! Horizontal slab
78         !                                             ! ===============
79         DO jj = 2, jpjm1
80            DO ji = fs_2, fs_jpim1   ! vector opt.
81               ze2u = rotb (ji,jj,jk) * fsahmf(ji,jj,jk) * fse3f(ji,jj,jk)
82               ze1v = hdivb(ji,jj,jk) * fsahmt(ji,jj,jk)
83               ! horizontal diffusive trends
84               zua = - ( ze2u - rotb (ji,jj-1,jk)*fsahmf(ji,jj-1,jk)*fse3f(ji,jj-1,jk) ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )   &
85                     + ( hdivb(ji+1,jj,jk)*fsahmt(ji+1,jj,jk) - ze1v                   ) / e1u(ji,jj)
86
87               zva = + ( ze2u - rotb (ji-1,jj,jk)*fsahmf(ji-1,jj,jk)*fse3f(ji-1,jj,jk) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )   &
88                     + ( hdivb(ji,jj+1,jk)*fsahmt(ji,jj+1,jk) - ze1v                   ) / e2v(ji,jj)
89
90               ! add it to the general momentum trends
91               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
92               va(ji,jj,jk) = va(ji,jj,jk) + zva
93            END DO
94         END DO
95         !                                             ! ===============
96      END DO                                           !   End of slab
97      !                                                ! ===============
98     
99      IF( iom_use('dispkexyfo') ) THEN   ! ocean Kinetic Energy dissipation per unit area
100         !                               ! due to lateral friction (xy=lateral)
101         !   see NEMO_book appendix C, §C.7.2 (N.B. here averaged at t-points)
102         !   local dissipation of KE at t-point due to laplacian operator is given by :
103         !      - ahmt hdivb**2 - mi( mj(ahmf rotb**2 e1f*e2f*e3t) ) / (e1e2t*e3t)
104         !
105         ALLOCATE( z2d(jpi,jpj) )
106         z2d(:,:) = 0._wp
107         DO jk = 1, jpkm1
108            DO jj = 2, jpjm1
109               DO ji = 2, jpim1
110                  z2d(ji,jj) = z2d(ji,jj)          -   (                                                         &
111                     &   hdivb(ji,jj,jk)**2 * fsahmt(ji,jj,jk) * fse3t_n(ji,jj,jk)                               &
112                     & + 0.25_wp * (                                                                             &
113                     &   rotb (ji  ,jj  ,jk)**2 * fsahmf(ji  ,jj  ,jk) * e12f(ji  ,jj  ) * fse3f(ji  ,jj  ,jk)   &
114                     & + rotb (ji-1,jj  ,jk)**2 * fsahmf(ji-1,jj  ,jk) * e12f(ji-1,jj  ) * fse3f(ji-1,jj  ,jk)   &
115                     & + rotb (ji  ,jj-1,jk)**2 * fsahmf(ji  ,jj-1,jk) * e12f(ji  ,jj-1) * fse3f(ji  ,jj-1,jk)   &
116                     & + rotb (ji-1,jj-1,jk)**2 * fsahmf(ji-1,jj-1,jk) * e12f(ji-1,jj-1) * fse3f(ji-1,jj-1,jk)   &
117                     &             ) * r1_e12t(ji,jj)  ) * tmask(ji,jj,jk)
118               END DO
119            END DO
120         END DO
121         z2d(:,:) = rau0 * z2d(:,:) 
122         CALL lbc_lnk( z2d,'T', 1. )
123         CALL iom_put( 'dispkexyfo', z2d )
124         DEALLOCATE( z2d )
125      ENDIF
126
127      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_lap')
128      !
129   END SUBROUTINE dyn_ldf_lap
130
131   !!======================================================================
132END MODULE dynldf_lap
Note: See TracBrowser for help on using the repository browser.