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/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90 @ 8197

Last change on this file since 8197 was 8197, checked in by glong, 7 years ago

Changed id's to be chars e.g. hpg to more easily identify output (and updated field_def.xml accordingly). Also rearranged scaling factors in dyn_vrt_dia subroutines in divcur.F90
File size: 9.9 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   USE divcur          ! for dyn_vrt_dia_3d
22   !
23   USE in_out_manager  ! I/O manager
24   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
25   USE wrk_nemo        ! Memory Allocation
26   USE timing          ! Timing
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   dyn_ldf_bilap   ! called by step.F90
32
33   !! * Substitutions
34#  include "domzgr_substitute.h90"
35#  include "ldfdyn_substitute.h90"
36#  include "vectopt_loop_substitute.h90"
37   !!----------------------------------------------------------------------
38   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
39   !! $Id$
40   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42CONTAINS
43
44   SUBROUTINE dyn_ldf_bilap( kt )
45      !!----------------------------------------------------------------------
46      !!                     ***  ROUTINE dyn_ldf_bilap  ***
47      !!
48      !! ** Purpose :   Compute the before trend of the lateral momentum
49      !!      diffusion and add it to the general trend of momentum equation.
50      !!
51      !! ** Method  :   The before horizontal momentum diffusion trend is a
52      !!      bi-harmonic operator (bilaplacian type) which separates the
53      !!      divergent and rotational parts of the flow.
54      !!      Its horizontal components are computed as follow:
55      !!      laplacian:
56      !!          zlu = 1/e1u di[ hdivb ] - 1/(e2u*e3u) dj-1[ e3f rotb ]
57      !!          zlv = 1/e2v dj[ hdivb ] + 1/(e1v*e3v) di-1[ e3f rotb ]
58      !!      third derivative:
59      !!       * multiply by the eddy viscosity coef. at u-, v-point, resp.
60      !!          zlu = ahmu * zlu
61      !!          zlv = ahmv * zlv
62      !!       * curl and divergence of the laplacian
63      !!          zuf = 1/(e1f*e2f) ( di[e2v zlv] - dj[e1u zlu] )
64      !!          zut = 1/(e1t*e2t*e3t) ( di[e2u*e3u zlu] + dj[e1v*e3v zlv] )
65      !!      bilaplacian:
66      !!              diffu = 1/e1u di[ zut ] - 1/(e2u*e3u) dj-1[ e3f zuf ]
67      !!              diffv = 1/e2v dj[ zut ] + 1/(e1v*e3v) di-1[ e3f zuf ]
68      !!      If ln_sco=F and ln_zps=F, the vertical scale factors in the
69      !!      rotational part of the diffusion are simplified
70      !!      Add this before trend to the general trend (ua,va):
71      !!            (ua,va) = (ua,va) + (diffu,diffv)
72      !!
73      !! ** Action : - Update (ua,va) with the before iso-level biharmonic
74      !!               mixing trend.
75      !!----------------------------------------------------------------------
76      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
77      !
78      CHARACTER(len=3) :: id_vrt_dia_ldf = "ldf"  ! TODO remove once flags set properly
79      INTEGER  ::   ji, jj, jk                  ! dummy loop indices
80      REAL(wp) ::   zbt, ze2u, ze2v             ! temporary scalar
81      REAL(wp), POINTER, DIMENSION(:,:  ) :: zcu, zcv
82      REAL(wp), POINTER, DIMENSION(:,:,:) :: zuf, zut, zlu, zlv, zua, zva
83      !!----------------------------------------------------------------------
84      !
85      IF( nn_timing == 1 )  CALL timing_start('dyn_ldf_bilap')
86      !
87      CALL wrk_alloc( jpi, jpj,      zcu, zcv           )
88      CALL wrk_alloc( jpi, jpj, jpk, zuf, zut, zlu, zlv, zua, zva ) 
89      !
90      IF( kt == nit000 .AND. lwp ) THEN
91         WRITE(numout,*)
92         WRITE(numout,*) 'dyn_ldf_bilap : iso-level bilaplacian operator'
93         WRITE(numout,*) '~~~~~~~~~~~~~'
94      ENDIF
95
96!!bug gm this should be enough
97!!$      zuf(:,:,jpk) = 0.e0
98!!$      zut(:,:,jpk) = 0.e0
99!!$      zlu(:,:,jpk) = 0.e0
100!!$      zlv(:,:,jpk) = 0.e0
101      zuf(:,:,:) = 0._wp
102      zut(:,:,:) = 0._wp
103      zlu(:,:,:) = 0._wp
104      zlv(:,:,:) = 0._wp
105
106      !                                                ! ===============
107      DO jk = 1, jpkm1                                 ! Horizontal slab
108         !                                             ! ===============
109         ! Laplacian
110         ! ---------
111
112         IF( ln_sco .OR. ln_zps ) THEN   ! s-coordinate or z-coordinate with partial steps
113            zuf(:,:,jk) = rotb(:,:,jk) * fse3f(:,:,jk)
114            DO jj = 2, jpjm1
115               DO ji = fs_2, fs_jpim1   ! vector opt.
116                  zlu(ji,jj,jk) = - ( zuf(ji,jj,jk) - zuf(ji,jj-1,jk) ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )   &
117                     &         + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj,jk) ) / e1u(ji,jj)
118   
119                  zlv(ji,jj,jk) = + ( zuf(ji,jj,jk) - zuf(ji-1,jj,jk) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )   &
120                     &         + ( hdivb(ji,jj+1,jk) - hdivb(ji,jj,jk) ) / e2v(ji,jj)
121               END DO
122            END DO
123         ELSE                            ! z-coordinate - full step
124            DO jj = 2, jpjm1
125               DO ji = fs_2, fs_jpim1   ! vector opt.
126                  zlu(ji,jj,jk) = - ( rotb (ji  ,jj,jk) - rotb (ji,jj-1,jk) ) / e2u(ji,jj)   &
127                     &         + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj  ,jk) ) / e1u(ji,jj)
128   
129                  zlv(ji,jj,jk) = + ( rotb (ji,jj  ,jk) - rotb (ji-1,jj,jk) ) / e1v(ji,jj)   &
130                     &         + ( hdivb(ji,jj+1,jk) - hdivb(ji  ,jj,jk) ) / e2v(ji,jj)
131               END DO 
132            END DO 
133         ENDIF
134      END DO
135      CALL lbc_lnk( zlu, 'U', -1. )   ;   CALL lbc_lnk( zlv, 'V', -1. )   ! Boundary conditions
136
137         
138      DO jk = 1, jpkm1
139   
140         ! Third derivative
141         ! ----------------
142         
143         ! Multiply by the eddy viscosity coef. (at u- and v-points)
144         zlu(:,:,jk) = zlu(:,:,jk) * ( fsahmu(:,:,jk) * (1-nkahm_smag) + nkahm_smag)
145
146         zlv(:,:,jk) = zlv(:,:,jk) * ( fsahmv(:,:,jk) * (1-nkahm_smag) + nkahm_smag)
147         
148         ! Contravariant "laplacian"
149         zcu(:,:) = e1u(:,:) * zlu(:,:,jk)
150         zcv(:,:) = e2v(:,:) * zlv(:,:,jk)
151         
152         ! Laplacian curl ( * e3f if s-coordinates or z-coordinate with partial steps)
153         DO jj = 1, jpjm1
154            DO ji = 1, fs_jpim1   ! vector opt.
155               zuf(ji,jj,jk) = fmask(ji,jj,jk) * (  zcv(ji+1,jj  ) - zcv(ji,jj)      &
156                  &                            - zcu(ji  ,jj+1) + zcu(ji,jj)  )   &
157                  &       * fse3f(ji,jj,jk) / ( e1f(ji,jj)*e2f(ji,jj) )
158            END DO 
159         END DO 
160
161         ! Laplacian Horizontal fluxes
162         DO jj = 1, jpjm1
163            DO ji = 1, fs_jpim1   ! vector opt.
164               zlu(ji,jj,jk) = e2u(ji,jj) * fse3u(ji,jj,jk) * zlu(ji,jj,jk)
165               zlv(ji,jj,jk) = e1v(ji,jj) * fse3v(ji,jj,jk) * zlv(ji,jj,jk)
166            END DO
167         END DO
168
169         ! Laplacian divergence
170         DO jj = 2, jpj
171            DO ji = fs_2, jpi   ! vector opt.
172               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
173               zut(ji,jj,jk) = (  zlu(ji,jj,jk) - zlu(ji-1,jj  ,jk)   &
174                  &             + zlv(ji,jj,jk) - zlv(ji  ,jj-1,jk) ) / zbt
175            END DO
176         END DO
177      END DO
178
179
180      ! boundary conditions on the laplacian curl and div (zuf,zut)
181!!bug gm no need to do this 2 following lbc...
182      CALL lbc_lnk( zuf, 'F', 1. )
183      CALL lbc_lnk( zut, 'T', 1. )
184
185      DO jk = 1, jpkm1     
186   
187         ! Bilaplacian
188         ! -----------
189
190         DO jj = 2, jpjm1
191            DO ji = fs_2, fs_jpim1   ! vector opt.
192               ze2u = e2u(ji,jj) * fse3u(ji,jj,jk)
193               ze2v = e1v(ji,jj) * fse3v(ji,jj,jk)
194               ! horizontal biharmonic diffusive trends multiplied by the
195               ! eddy viscosity coef. (at u- and v-points)
196               zua(ji,jj,jk) = ( - ( zuf(ji  ,jj,jk) - zuf(ji,jj-1,jk) ) / ze2u   &
197                  &  + ( zut(ji+1,jj,jk) - zut(ji,jj  ,jk) ) / e1u(ji,jj) )       &
198                  &  * ( fsahmu(ji,jj,jk)*nkahm_smag +(1 -nkahm_smag ))
199               zva(ji,jj,jk) = ( + ( zuf(ji,jj  ,jk) - zuf(ji-1,jj,jk) ) / ze2v   &
200                  &  + ( zut(ji,jj+1,jk) - zut(ji  ,jj,jk) ) / e2v(ji,jj) )       &
201                  &  * ( fsahmv(ji,jj,jk)*nkahm_smag +(1 -nkahm_smag ))
202               ! add it to the general momentum trends
203               ua(ji,jj,jk) = ua(ji,jj,jk) + zua(ji,jj,jk)
204               va(ji,jj,jk) = va(ji,jj,jk) + zva(ji,jj,jk)
205            END DO
206         END DO
207
208         !                                             ! ===============
209      END DO                                           !   End of slab
210      !                                                ! ===============
211      IF ( id_vrt_dia_ldf == "ldf" ) THEN
212          ! TODO - remove kt only used for validation
213          CALL dyn_vrt_dia_3d(zua, zva, id_vrt_dia_ldf, kt)
214      END IF
215      !
216      CALL wrk_dealloc( jpi, jpj,      zcu, zcv           )
217      CALL wrk_dealloc( jpi, jpj, jpk, zuf, zut, zlu, zlv, zua, zva ) 
218      !
219      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_bilap')
220      !
221   END SUBROUTINE dyn_ldf_bilap
222
223   !!======================================================================
224END MODULE dynldf_bilap
Note: See TracBrowser for help on using the repository browser.