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

source: branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90 @ 2887

Last change on this file since 2887 was 2887, checked in by hliu, 13 years ago

addition and modification of files for Smagorinsky method. for Maria Luneva
  • Property svn:keywords set to Id
File size: 10.1 KB
Line 
1MODULE dynldf_bilap
2   !!======================================================================
3   !!                     ***  MODULE  dynldf_bilap  ***
4   !! Ocean dynamics:  lateral viscosity trend
5   !!======================================================================
6
7   !!----------------------------------------------------------------------
8   !!   dyn_ldf_bilap : update the momentum trend with the lateral diffusion
9   !!                   using an iso-level bilaplacian operator
10   !!----------------------------------------------------------------------
11   !! * Modules used
12   USE oce             ! ocean dynamics and tracers
13   USE dom_oce         ! ocean space and time domain
14   USE ldfdyn_oce      ! ocean dynamics: lateral physics
15   USE in_out_manager  ! I/O manager
16   USE trdmod          ! ocean dynamics trends
17   USE trdmod_oce      ! ocean variables trends
18   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
19
20   IMPLICIT NONE
21   PRIVATE
22
23   !! * Routine accessibility
24   PUBLIC dyn_ldf_bilap  ! called by step.F90
25
26   !! * Substitutions
27#  include "domzgr_substitute.h90"
28#  include "ldfdyn_substitute.h90"
29#  include "vectopt_loop_substitute.h90"
30   !!----------------------------------------------------------------------
31   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
32   !! $Id$
33   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
34   !!----------------------------------------------------------------------
35
36CONTAINS
37
38   SUBROUTINE dyn_ldf_bilap( kt )
39      !!----------------------------------------------------------------------
40      !!                     ***  ROUTINE dyn_ldf_bilap  ***
41      !!
42      !! ** Purpose :   Compute the before trend of the lateral momentum
43      !!      diffusion and add it to the general trend of momentum equation.
44      !!
45      !! ** Method  :   The before horizontal momentum diffusion trend is a
46      !!      bi-harmonic operator (bilaplacian type) which separates the
47      !!      divergent and rotational parts of the flow.
48      !!      Its horizontal components are computed as follow:
49      !!      laplacian:
50      !!          zlu = 1/e1u di[ hdivb ] - 1/(e2u*e3u) dj-1[ e3f rotb ]
51      !!          zlv = 1/e2v dj[ hdivb ] + 1/(e1v*e3v) di-1[ e3f rotb ]
52      !!      third derivative:
53      !!       * multiply by the eddy viscosity coef. at u-, v-point, resp.
54      !!          zlu = ahmu * zlu
55      !!          zlv = ahmv * zlv
56      !!       * curl and divergence of the laplacian
57      !!          zuf = 1/(e1f*e2f) ( di[e2v zlv] - dj[e1u zlu] )
58      !!          zut = 1/(e1t*e2t*e3t) ( di[e2u*e3u zlu] + dj[e1v*e3v zlv] )
59      !!      bilaplacian:
60      !!              diffu = 1/e1u di[ zut ] - 1/(e2u*e3u) dj-1[ e3f zuf ]
61      !!              diffv = 1/e2v dj[ zut ] + 1/(e1v*e3v) di-1[ e3f zuf ]
62      !!      If ln_sco=F and ln_zps=F, the vertical scale factors in the
63      !!      rotational part of the diffusion are simplified
64      !!      Add this before trend to the general trend (ua,va):
65      !!            (ua,va) = (ua,va) + (diffu,diffv)
66      !!      'key_trddyn' defined: the two components of the horizontal
67      !!                               diffusion trend are saved.
68      !!
69      !! ** Action : - Update (ua,va) with the before iso-level biharmonic
70      !!               mixing trend.
71      !!
72      !! History :
73      !!        !  90-09  (G. Madec)  Original code
74      !!        !  91-11  (G. Madec)
75      !!        !  93-03  (M. Guyon)  symetrical conditions (M. Guyon)
76      !!        !  96-01  (G. Madec)  statement function for e3
77      !!        !  97-07  (G. Madec)  lbc calls
78      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
79      !!   9.0  !  04-08  (C. Talandier) New trends organization
80      !!----------------------------------------------------------------------
81      !! * Arguments
82      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
83
84      !! * Local declarations
85      INTEGER  ::   ji, jj, jk                ! dummy loop indices
86      REAL(wp) ::   zua, zva, zbt, ze2u, ze2v ! temporary scalar
87      REAL(wp), DIMENSION(jpi,jpj) ::   &
88         zcu, zcv                             ! temporary workspace
89      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
90         zuf, zut, zlu, zlv                   ! temporary workspace
91      !!----------------------------------------------------------------------
92      !!  OPA 8.5, LODYC-IPSL (2002)
93      !!----------------------------------------------------------------------
94
95      IF( kt == nit000 ) THEN
96         IF(lwp) WRITE(numout,*)
97         IF(lwp) WRITE(numout,*) 'dyn_ldf_bilap : iso-level bilaplacian operator'
98         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
99      ENDIF
100
101!!bug gm this should be enough
102!!$      zuf(:,:,jpk) = 0.e0
103!!$      zut(:,:,jpk) = 0.e0
104!!$      zlu(:,:,jpk) = 0.e0
105!!$      zlv(:,:,jpk) = 0.e0
106      zuf(:,:,:) = 0.e0
107      zut(:,:,:) = 0.e0
108      zlu(:,:,:) = 0.e0
109      zlv(:,:,:) = 0.e0
110
111      !                                                ! ===============
112      DO jk = 1, jpkm1                                 ! Horizontal slab
113         !                                             ! ===============
114         ! Laplacian
115         ! ---------
116
117         IF( ln_sco .OR. ln_zps ) THEN   ! s-coordinate or z-coordinate with partial steps
118            zuf(:,:,jk) = rotb(:,:,jk) * fse3f(:,:,jk)
119            DO jj = 2, jpjm1
120               DO ji = fs_2, fs_jpim1   ! vector opt.
121                  zlu(ji,jj,jk) = - ( zuf(ji,jj,jk) - zuf(ji,jj-1,jk) ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )   &
122                     &         + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj,jk) ) / e1u(ji,jj)
123   
124                  zlv(ji,jj,jk) = + ( zuf(ji,jj,jk) - zuf(ji-1,jj,jk) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )   &
125                     &         + ( hdivb(ji,jj+1,jk) - hdivb(ji,jj,jk) ) / e2v(ji,jj)
126               END DO
127            END DO
128         ELSE                            ! z-coordinate - full step
129            DO jj = 2, jpjm1
130               DO ji = fs_2, fs_jpim1   ! vector opt.
131                  zlu(ji,jj,jk) = - ( rotb (ji  ,jj,jk) - rotb (ji,jj-1,jk) ) / e2u(ji,jj)   &
132                     &         + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj  ,jk) ) / e1u(ji,jj)
133   
134                  zlv(ji,jj,jk) = + ( rotb (ji,jj  ,jk) - rotb (ji-1,jj,jk) ) / e1v(ji,jj)   &
135                     &         + ( hdivb(ji,jj+1,jk) - hdivb(ji  ,jj,jk) ) / e2v(ji,jj)
136               END DO 
137            END DO 
138         ENDIF
139      ENDDO
140
141      ! Boundary conditions on the laplacian  (zlu,zlv)
142      CALL lbc_lnk( zlu, 'U', -1. )
143      CALL lbc_lnk( zlv, 'V', -1. )
144         
145         
146      DO jk = 1, jpkm1
147   
148         ! Third derivative
149         ! ----------------
150         
151         ! Multiply by the eddy viscosity coef. (at u- and v-points)
152#if ! defined key_dynldf_smag
153         zlu(:,:,jk) = zlu(:,:,jk) * fsahmu(:,:,jk)
154         zlv(:,:,jk) = zlv(:,:,jk) * fsahmv(:,:,jk)
155#endif         
156         ! Contravariant "laplacian"
157         zcu(:,:) = e1u(:,:) * zlu(:,:,jk)
158         zcv(:,:) = e2v(:,:) * zlv(:,:,jk)
159         
160         ! Laplacian curl ( * e3f if s-coordinates or z-coordinate with partial steps)
161         DO jj = 1, jpjm1
162            DO ji = 1, fs_jpim1   ! vector opt.
163               zuf(ji,jj,jk) = fmask(ji,jj,jk) * (  zcv(ji+1,jj  ) - zcv(ji,jj)      &
164                  &                            - zcu(ji  ,jj+1) + zcu(ji,jj)  )   &
165                  &       * fse3f(ji,jj,jk) / ( e1f(ji,jj)*e2f(ji,jj) )
166            END DO 
167         END DO 
168
169         ! Laplacian Horizontal fluxes
170         DO jj = 1, jpjm1
171            DO ji = 1, fs_jpim1   ! vector opt.
172               zlu(ji,jj,jk) = e2u(ji,jj) * fse3u(ji,jj,jk) * zlu(ji,jj,jk)
173               zlv(ji,jj,jk) = e1v(ji,jj) * fse3v(ji,jj,jk) * zlv(ji,jj,jk)
174            END DO
175         END DO
176
177         ! Laplacian divergence
178         DO jj = 2, jpj
179            DO ji = fs_2, jpi   ! vector opt.
180               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
181               zut(ji,jj,jk) = (  zlu(ji,jj,jk) - zlu(ji-1,jj  ,jk)   &
182                  &             + zlv(ji,jj,jk) - zlv(ji  ,jj-1,jk) ) / zbt
183            END DO
184         END DO
185      END DO
186
187
188      ! boundary conditions on the laplacian curl and div (zuf,zut)
189!!bug gm no need to do this 2 following lbc...
190      CALL lbc_lnk( zuf, 'F', 1. )
191      CALL lbc_lnk( zut, 'T', 1. )
192
193      DO jk = 1, jpkm1     
194   
195         ! Bilaplacian
196         ! -----------
197#if !defined key_dynldf_smag
198         DO jj = 2, jpjm1
199            DO ji = fs_2, fs_jpim1   ! vector opt.
200               ze2u = e2u(ji,jj) * fse3u(ji,jj,jk)
201               ze2v = e1v(ji,jj) * fse3v(ji,jj,jk)
202               ! horizontal biharmonic diffusive trends
203               zua = - ( zuf(ji  ,jj,jk) - zuf(ji,jj-1,jk) ) / ze2u   &
204                  &  + ( zut(ji+1,jj,jk) - zut(ji,jj  ,jk) ) / e1u(ji,jj)
205
206               zva = + ( zuf(ji,jj  ,jk) - zuf(ji-1,jj,jk) ) / ze2v   &
207                  &  + ( zut(ji,jj+1,jk) - zut(ji  ,jj,jk) ) / e2v(ji,jj)
208               ! add it to the general momentum trends
209               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
210               va(ji,jj,jk) = va(ji,jj,jk) + zva
211            END DO
212         END DO
213#else
214
215         DO jj = 2, jpjm1
216            DO ji = fs_2, fs_jpim1   ! vector opt.
217               ze2u = e2u(ji,jj) * fse3u(ji,jj,jk)
218               ze2v = e1v(ji,jj) * fse3v(ji,jj,jk)
219               ! horizontal biharmonic diffusive trends
220               zua = - ( zuf(ji  ,jj,jk) - zuf(ji,jj-1,jk) ) / ze2u   &
221                  &  + ( zut(ji+1,jj,jk) - zut(ji,jj  ,jk) ) / e1u(ji,jj)
222
223               zva = + ( zuf(ji,jj  ,jk) - zuf(ji-1,jj,jk) ) / ze2v   &
224                  &  + ( zut(ji,jj+1,jk) - zut(ji  ,jj,jk) ) / e2v(ji,jj)
225               ! add it to the general momentum trends
226               ua(ji,jj,jk) = ua(ji,jj,jk) + zua*fsahmu(ji,jj,jk)
227               va(ji,jj,jk) = va(ji,jj,jk) + zva*fsahmv(ji,jj,jk)
228            END DO
229         END DO
230
231#endif
232         !                                             ! ===============
233      END DO                                           !   End of slab
234      !                                                ! ===============
235
236   END SUBROUTINE dyn_ldf_bilap
237
238   !!======================================================================
239END MODULE dynldf_bilap
Note: See TracBrowser for help on using the repository browser.