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.
dynldf_bilap.F90 in branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90 @ 14555

Last change on this file since 14555 was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.
File size: 9.4 KB
Line 
1MODULE dynldf_bilap
2   !!======================================================================
3   !!                     ***  MODULE  dynldf_bilap  ***
4   !! Ocean dynamics:  lateral viscosity trend
5   !!======================================================================
6   !! History :  OPA  ! 1990-09  (G. Madec)  Original code
7   !!            4.0  ! 1993-03  (M. Guyon)  symetrical conditions (M. Guyon)
8   !!            6.0  ! 1996-01  (G. Madec)  statement function for e3
9   !!            8.0  ! 1997-07  (G. Madec)  lbc calls
10   !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90: Free form and module
11   !!            2.0  ! 2004-08  (C. Talandier) New trends organization
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   dyn_ldf_bilap : update the momentum trend with the lateral diffusion
16   !!                   using an iso-level bilaplacian operator
17   !!----------------------------------------------------------------------
18   USE oce             ! ocean dynamics and tracers
19   USE dom_oce         ! ocean space and time domain
20   USE ldfdyn_oce      ! ocean dynamics: lateral physics
21   !
22   USE in_out_manager  ! I/O manager
23   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
24   USE wrk_nemo        ! Memory Allocation
25   USE timing          ! Timing
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   dyn_ldf_bilap   ! 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_bilap( kt )
44      !!----------------------------------------------------------------------
45      !!                     ***  ROUTINE dyn_ldf_bilap  ***
46      !!
47      !! ** Purpose :   Compute the before trend of the lateral momentum
48      !!      diffusion and add it to the general trend of momentum equation.
49      !!
50      !! ** Method  :   The before horizontal momentum diffusion trend is a
51      !!      bi-harmonic operator (bilaplacian type) which separates the
52      !!      divergent and rotational parts of the flow.
53      !!      Its horizontal components are computed as follow:
54      !!      laplacian:
55      !!          zlu = 1/e1u di[ hdivb ] - 1/(e2u*e3u) dj-1[ e3f rotb ]
56      !!          zlv = 1/e2v dj[ hdivb ] + 1/(e1v*e3v) di-1[ e3f rotb ]
57      !!      third derivative:
58      !!       * multiply by the eddy viscosity coef. at u-, v-point, resp.
59      !!          zlu = ahmu * zlu
60      !!          zlv = ahmv * zlv
61      !!       * curl and divergence of the laplacian
62      !!          zuf = 1/(e1f*e2f) ( di[e2v zlv] - dj[e1u zlu] )
63      !!          zut = 1/(e1t*e2t*e3t) ( di[e2u*e3u zlu] + dj[e1v*e3v zlv] )
64      !!      bilaplacian:
65      !!              diffu = 1/e1u di[ zut ] - 1/(e2u*e3u) dj-1[ e3f zuf ]
66      !!              diffv = 1/e2v dj[ zut ] + 1/(e1v*e3v) di-1[ e3f zuf ]
67      !!      If ln_sco=F and ln_zps=F, the vertical scale factors in the
68      !!      rotational part of the diffusion are simplified
69      !!      Add this before trend to the general trend (ua,va):
70      !!            (ua,va) = (ua,va) + (diffu,diffv)
71      !!
72      !! ** Action : - Update (ua,va) with the before iso-level biharmonic
73      !!               mixing trend.
74      !!----------------------------------------------------------------------
75      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
76      !
77      INTEGER  ::   ji, jj, jk                  ! dummy loop indices
78      REAL(wp) ::   zua, zva, zbt, ze2u, ze2v   ! temporary scalar
79      REAL(wp), POINTER, DIMENSION(:,:  ) :: zcu, zcv
80      REAL(wp), POINTER, DIMENSION(:,:,:) :: zuf, zut, zlu, zlv
81      !!----------------------------------------------------------------------
82      !
83      IF( nn_timing == 1 )  CALL timing_start('dyn_ldf_bilap')
84      !
85      CALL wrk_alloc( jpi, jpj,      zcu, zcv           )
86      CALL wrk_alloc( jpi, jpj, jpk, zuf, zut, zlu, zlv ) 
87      !
88      IF( kt == nit000 .AND. lwp ) THEN
89         WRITE(numout,*)
90         WRITE(numout,*) 'dyn_ldf_bilap : iso-level bilaplacian operator'
91         WRITE(numout,*) '~~~~~~~~~~~~~'
92         IF(lwp .AND. lflush) CALL flush(numout)
93      ENDIF
94
95!!bug gm this should be enough
96!!$      zuf(:,:,jpk) = 0.e0
97!!$      zut(:,:,jpk) = 0.e0
98!!$      zlu(:,:,jpk) = 0.e0
99!!$      zlv(:,:,jpk) = 0.e0
100      zuf(:,:,:) = 0._wp
101      zut(:,:,:) = 0._wp
102      zlu(:,:,:) = 0._wp
103      zlv(:,:,:) = 0._wp
104
105      !                                                ! ===============
106      DO jk = 1, jpkm1                                 ! Horizontal slab
107         !                                             ! ===============
108         ! Laplacian
109         ! ---------
110
111         IF( ln_sco .OR. ln_zps ) THEN   ! s-coordinate or z-coordinate with partial steps
112            zuf(:,:,jk) = rotb(:,:,jk) * fse3f(:,:,jk)
113            DO jj = 2, jpjm1
114               DO ji = fs_2, fs_jpim1   ! vector opt.
115                  zlu(ji,jj,jk) = - ( zuf(ji,jj,jk) - zuf(ji,jj-1,jk) ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )   &
116                     &         + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj,jk) ) / e1u(ji,jj)
117   
118                  zlv(ji,jj,jk) = + ( zuf(ji,jj,jk) - zuf(ji-1,jj,jk) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )   &
119                     &         + ( hdivb(ji,jj+1,jk) - hdivb(ji,jj,jk) ) / e2v(ji,jj)
120               END DO
121            END DO
122         ELSE                            ! z-coordinate - full step
123            DO jj = 2, jpjm1
124               DO ji = fs_2, fs_jpim1   ! vector opt.
125                  zlu(ji,jj,jk) = - ( rotb (ji  ,jj,jk) - rotb (ji,jj-1,jk) ) / e2u(ji,jj)   &
126                     &         + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj  ,jk) ) / e1u(ji,jj)
127   
128                  zlv(ji,jj,jk) = + ( rotb (ji,jj  ,jk) - rotb (ji-1,jj,jk) ) / e1v(ji,jj)   &
129                     &         + ( hdivb(ji,jj+1,jk) - hdivb(ji  ,jj,jk) ) / e2v(ji,jj)
130               END DO 
131            END DO 
132         ENDIF
133      END DO
134      CALL lbc_lnk( zlu, 'U', -1. )   ;   CALL lbc_lnk( zlv, 'V', -1. )   ! Boundary conditions
135
136         
137      DO jk = 1, jpkm1
138   
139         ! Third derivative
140         ! ----------------
141         
142         ! Multiply by the eddy viscosity coef. (at u- and v-points)
143         zlu(:,:,jk) = zlu(:,:,jk) * ( fsahmu(:,:,jk) * (1-nkahm_smag) + nkahm_smag)
144
145         zlv(:,:,jk) = zlv(:,:,jk) * ( fsahmv(:,:,jk) * (1-nkahm_smag) + nkahm_smag)
146         
147         ! Contravariant "laplacian"
148         zcu(:,:) = e1u(:,:) * zlu(:,:,jk)
149         zcv(:,:) = e2v(:,:) * zlv(:,:,jk)
150         
151         ! Laplacian curl ( * e3f if s-coordinates or z-coordinate with partial steps)
152         DO jj = 1, jpjm1
153            DO ji = 1, fs_jpim1   ! vector opt.
154               zuf(ji,jj,jk) = fmask(ji,jj,jk) * (  zcv(ji+1,jj  ) - zcv(ji,jj)      &
155                  &                            - zcu(ji  ,jj+1) + zcu(ji,jj)  )   &
156                  &       * fse3f(ji,jj,jk) / ( e1f(ji,jj)*e2f(ji,jj) )
157            END DO 
158         END DO 
159
160         ! Laplacian Horizontal fluxes
161         DO jj = 1, jpjm1
162            DO ji = 1, fs_jpim1   ! vector opt.
163               zlu(ji,jj,jk) = e2u(ji,jj) * fse3u(ji,jj,jk) * zlu(ji,jj,jk)
164               zlv(ji,jj,jk) = e1v(ji,jj) * fse3v(ji,jj,jk) * zlv(ji,jj,jk)
165            END DO
166         END DO
167
168         ! Laplacian divergence
169         DO jj = 2, jpj
170            DO ji = fs_2, jpi   ! vector opt.
171               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
172               zut(ji,jj,jk) = (  zlu(ji,jj,jk) - zlu(ji-1,jj  ,jk)   &
173                  &             + zlv(ji,jj,jk) - zlv(ji  ,jj-1,jk) ) / zbt
174            END DO
175         END DO
176      END DO
177
178
179      ! boundary conditions on the laplacian curl and div (zuf,zut)
180!!bug gm no need to do this 2 following lbc...
181      CALL lbc_lnk( zuf, 'F', 1. )
182      CALL lbc_lnk( zut, 'T', 1. )
183
184      DO jk = 1, jpkm1     
185   
186         ! Bilaplacian
187         ! -----------
188
189         DO jj = 2, jpjm1
190            DO ji = fs_2, fs_jpim1   ! vector opt.
191               ze2u = e2u(ji,jj) * fse3u(ji,jj,jk)
192               ze2v = e1v(ji,jj) * fse3v(ji,jj,jk)
193               ! horizontal biharmonic diffusive trends
194               zua = - ( zuf(ji  ,jj,jk) - zuf(ji,jj-1,jk) ) / ze2u   &
195                  &  + ( zut(ji+1,jj,jk) - zut(ji,jj  ,jk) ) / e1u(ji,jj)
196
197               zva = + ( zuf(ji,jj  ,jk) - zuf(ji-1,jj,jk) ) / ze2v   &
198                  &  + ( zut(ji,jj+1,jk) - zut(ji  ,jj,jk) ) / e2v(ji,jj)
199               ! add it to the general momentum trends
200               ua(ji,jj,jk) = ua(ji,jj,jk) + zua * ( fsahmu(ji,jj,jk)*nkahm_smag +(1 -nkahm_smag ))
201               va(ji,jj,jk) = va(ji,jj,jk) + zva * ( fsahmv(ji,jj,jk)*nkahm_smag +(1 -nkahm_smag ))
202            END DO
203         END DO
204
205         !                                             ! ===============
206      END DO                                           !   End of slab
207      !                                                ! ===============
208      CALL wrk_dealloc( jpi, jpj,      zcu, zcv           )
209      CALL wrk_dealloc( jpi, jpj, jpk, zuf, zut, zlu, zlv ) 
210      !
211      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_bilap')
212      !
213   END SUBROUTINE dyn_ldf_bilap
214
215   !!======================================================================
216END MODULE dynldf_bilap
Note: See TracBrowser for help on using the repository browser.