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_bilapg.F90 in trunk/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90 @ 3

Last change on this file since 3 was 3, checked in by opalod, 20 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 21.2 KB
Line 
1MODULE dynldf_bilapg
2   !!======================================================================
3   !!                       ***  MODULE  dynldf_bilapg  ***
4   !! Ocean dynamics:  lateral viscosity trend
5   !!======================================================================
6#if defined key_ldfslp   ||   defined key_esopa
7   !!----------------------------------------------------------------------
8   !!   'key_ldfslp'                              Rotation of mixing tensor
9   !!----------------------------------------------------------------------
10   !!   dyn_ldf_bilapg : update the momentum trend with the horizontal part
11   !!                    of the horizontal s-coord. bilaplacian diffusion
12   !!   ldfguv         :
13   !!----------------------------------------------------------------------
14   !! * Modules used
15   USE oce             ! ocean dynamics and tracers
16   USE dom_oce         ! ocean space and time domain
17   USE ldfdyn_oce      ! ocean dynamics lateral physics
18   USE zdf_oce         ! ocean vertical physics
19   USE in_out_manager  ! I/O manager
20   USE trddyn_oce      ! dynamics trends diagnostics variables
21   USE ldfslp          ! iso-neutral slopes available
22   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
23
24   IMPLICIT NONE
25   PRIVATE
26
27   !! * Routine accessibility
28   PUBLIC dyn_ldf_bilapg ! called by step.F90
29
30   !! * Substitutions
31#  include "domzgr_substitute.h90"
32#  include "ldfdyn_substitute.h90"
33   !!----------------------------------------------------------------------
34   !!   OPA 9.0 , LODYC-IPSL  (2003)
35   !!----------------------------------------------------------------------
36
37CONTAINS
38
39   SUBROUTINE dyn_ldf_bilapg( kt )
40      !!----------------------------------------------------------------------
41      !!                   ***  ROUTINE dyn_ldf_bilapg  ***
42      !!                     
43      !! ** Purpose :   Compute the before trend of the horizontal momentum
44      !!      diffusion and add it to the general trend of momentum equation.
45      !!
46      !! ** Method  :   The lateral momentum diffusive trends is provided by a
47      !!      a 4th order operator rotated along geopotential surfaces. It is
48      !!      computed using before fields (forward in time) and geopotential
49      !!      slopes computed in routine inildf.
50      !!         -1- compute the geopotential harmonic operator applied to
51      !!      (ub,vb) and multiply it by the eddy diffusivity coefficient
52      !!      (done by a call to ldfgpu and ldfgpv routines) The result is in
53      !!      (wk1,wk2) arrays. Applied the domain lateral boundary conditions
54      !!      by call to lbc_lnk.
55      !!         -2- applied to (wk1,wk2) the geopotential harmonic operator
56      !!      by a second call to ldfgpu and ldfgpv routines respectively. The
57      !!      result is in (wk3,wk4) arrays.
58      !!         -3- Add this trend to the general trend (ta,sa):
59      !!            (ua,va) = (ua,va) + (wk3,wk4)
60      !!      'key_trddyn' defined: the trend is saved for diagnostics.
61      !!
62      !! ** Action  : - Update (ua,va) arrays with the before geopotential
63      !!                biharmonic mixing trend.
64      !!              - save the trend in (utrd,vtrd) ('key_diatrends')
65      !!
66      !! History :
67      !!   8.0  !  97-07  (G. Madec)  Original code
68      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
69      !!----------------------------------------------------------------------
70      !! * Arguments
71      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
72
73      !! * Local declarations
74      INTEGER ::   ji, jj, jk                 ! dummy loop indices
75      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
76         wk1, wk2,            &  ! work array used for rotated biharmonic
77         wk3, wk4                ! operator on tracers and/or momentum
78      !!----------------------------------------------------------------------
79
80      IF( kt == nit000 ) THEN
81         IF(lwp) WRITE(numout,*)
82         IF(lwp) WRITE(numout,*) 'dyn_ldf_bilapg : horizontal biharmonic operator in s-coordinate'
83         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
84         wk1(:,:,:) = 0.e0   ;   wk3(:,:,:) = 0.e0
85         wk2(:,:,:) = 0.e0   ;   wk4(:,:,:) = 0.e0
86      ENDIF
87
88      ! Laplacian of (ub,vb) multiplied by ahm
89      ! -------------------------------------- 
90      ! rotated harmonic operator applied to (ub,vb)
91      !     and multiply by ahmu, ahmv (output in (wk1,wk2) )
92
93      CALL ldfguv ( ub, vb, wk1, wk2, 1 )
94
95
96      ! Lateral boundary conditions on (wk1,wk2)
97      CALL lbc_lnk( wk1, 'U', -1. )
98      CALL lbc_lnk( wk2, 'V', -1. )
99
100
101      ! Bilaplacian of (ub,vb)
102      ! ----------------------
103      ! rotated harmonic operator applied to (wk1,wk2) (output in (wk3,wk4) )
104
105      CALL ldfguv ( wk1, wk2, wk3, wk4, 2 )
106
107
108      ! Update the momentum trends           (j-slab :   2, jpj-1)
109      ! --------------------------
110      !                                                ! ===============
111      DO jj = 2, jpjm1                                 !  Vertical slab
112         !                                             ! ===============
113         DO jk = 1, jpkm1
114            DO ji = 2, jpim1
115               ! add the diffusive trend to the general momentum trends
116               ua(ji,jj,jk) = ua(ji,jj,jk) + wk3(ji,jj,jk)
117               va(ji,jj,jk) = va(ji,jj,jk) + wk4(ji,jj,jk)
118#if defined key_trddyn 
119               ! save the horizontal diffusive trends
120               utrd(ji,jj,jk,3) = wk3(ji,jj,jk)
121               vtrd(ji,jj,jk,3) = wk4(ji,jj,jk)
122#endif
123            END DO
124         END DO
125         !                                             ! ===============
126      END DO                                           !   End of slab
127      !                                                ! ===============
128
129   END SUBROUTINE dyn_ldf_bilapg
130
131
132   SUBROUTINE ldfguv( pu, pv, plu, plv, kahm )
133      !!----------------------------------------------------------------------
134      !!                  ***  ROUTINE ldfguv  ***
135      !!     
136      !! ** Purpose :   Apply a geopotential harmonic operator to (pu,pv)
137      !!      (defined at u- and v-points) and multiply it by the eddy
138      !!      viscosity coefficient (if kahm=1).
139      !!
140      !! ** Method  :   The harmonic operator rotated along geopotential
141      !!      surfaces is applied to (pu,pv) using the slopes of geopotential
142      !!      surfaces computed in inildf routine. The result is provided in
143      !!      (plu,plv) arrays. It is computed in 2 stepv:
144      !!
145      !!      First step: horizontal part of the operator. It is computed on
146      !!      ==========  pu as follows (idem on pv)
147      !!      horizontal fluxes :
148      !!         zftu = e2u*e3u/e1u di[ pu ] - e2u*uslp dk[ mi(mk(pu)) ]
149      !!         zftv = e1v*e3v/e2v dj[ pu ] - e1v*vslp dk[ mj(mk(pu)) ]
150      !!      take the horizontal divergence of the fluxes (no divided by
151      !!      the volume element :
152      !!         plu  = di-1[ zftu ] +  dj-1[ zftv ]
153      !!
154      !!      Second step: vertical part of the operator. It is computed on
155      !!      ===========  pu as follows (idem on pv)
156      !!      vertical fluxes :
157      !!         zftw = e1t*e2t/e3w * (wslpi^2+wslpj^2)  dk-1[ pu ]
158      !!              -     e2t     *       wslpi        di[ mi(mk(pu)) ]
159      !!              -     e1t     *       wslpj        dj[ mj(mk(pu)) ]
160      !!      take the vertical divergence of the fluxes add it to the hori-
161      !!      zontal component, divide the result by the volume element and
162      !!      if kahm=1, multiply by the eddy diffusivity coefficient:
163      !!         plu  = aht / (e1t*e2t*e3t) { plu + dk[ zftw ] }
164      !!      else:
165      !!         plu  =  1  / (e1t*e2t*e3t) { plu + dk[ zftw ] }
166      !!
167      !! ** Action :
168      !!        plu, plv        : partial harmonic operator applied to
169      !!                          pu and pv (all the components except
170      !!                          second order vertical derivative term)
171      !!      'key_trddyn' defined: the trend is saved for diagnostics.
172      !!
173      !! History :
174      !!   8.0  !  97-07  (G. Madec)  Original code
175      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
176      !!----------------------------------------------------------------------
177      !! * Arguments
178      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) ::   &
179         pu, pv     ! momentum fields (before u and v for the 1st call, and
180      !             ! laplacian of these fields multiplied by ahm for the 2nd
181      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) ::   &
182         plu, plv   ! partial harmonic operator applied to
183      !             ! pu and pv (all the components except
184      !             ! second order vertical derivative term)
185      INTEGER, INTENT( in ) ::   &
186         kahm       ! =1 the laplacian is multiplied by the eddy diffusivity coef.
187      !             ! =2 no multiplication
188
189      !! * Local declarations
190      INTEGER  ::   ji, jj, jk       ! dummy loop indices
191      REAL(wp) ::   &
192         zabe1, zabe2, zcof1, zcof2,    &  ! temporary scalars
193         zcoef0, zcoef3, zcoef4
194      REAL(wp) ::   &
195         zbur, zbvr, zmkt, zmkf, zuav, zvav,    &
196         zuwslpi, zuwslpj, zvwslpi, zvwslpj
197      REAL(wp), DIMENSION(jpi,jpj) ::   &
198         ziut, zjuf , zjvt, zivf,       &  ! workspace
199         zdku, zdk1u, zdkv, zdk1v
200      REAL(wp), DIMENSION(jpi,jpk) ::   &
201         zfuw, zfvw, zdiu, zdiv,        &  ! workspace
202         zdju, zdj1u, zdjv, zdj1v 
203      !!----------------------------------------------------------------------
204
205      !                               ! ********** !   ! ===============
206      DO jk = 1, jpkm1                ! First step !   ! Horizontal slab
207         !                            ! ********** !   ! ===============
208
209         ! I.1 Vertical gradient of pu and pv at level jk and jk+1
210         ! -------------------------------------------------------
211         ! surface boundary condition: zdku(jk=1)=zdku(jk=2)
212         !                             zdkv(jk=1)=zdkv(jk=2)
213
214         zdk1u(:,:) = ( pu(:,:,jk) - pu(:,:,jk+1) ) * umask(:,:,jk+1)
215         zdk1v(:,:) = ( pv(:,:,jk) - pv(:,:,jk+1) ) * vmask(:,:,jk+1)
216
217         IF( jk == 1 ) THEN
218            zdku(:,:) = zdk1u(:,:)
219            zdkv(:,:) = zdk1v(:,:)
220         ELSE
221            zdku(:,:) = ( pu(:,:,jk-1) - pu(:,:,jk) ) * umask(:,:,jk)
222            zdkv(:,:) = ( pv(:,:,jk-1) - pv(:,:,jk) ) * vmask(:,:,jk)
223         ENDIF
224
225         !                                -----f-----
226         ! I.2 Horizontal fluxes on U          |
227         ! ------------------------===     t   u   t
228         !                                     |
229         ! i-flux at t-point              -----f-----
230         DO jj = 1, jpjm1
231            DO ji = 2, jpi
232               zabe1 = e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)
233
234               zmkt  = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
235                              + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
236
237               zcof1 = -e2t(ji,jj) * zmkt   &
238                     * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
239
240               ziut(ji,jj) = tmask(ji,jj,jk) *   &
241                           (  zabe1 * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) )   &
242                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
243                                       +zdk1u(ji,jj) + zdku (ji-1,jj) )  )
244            END DO
245         END DO
246
247         ! j-flux at f-point
248         DO jj = 1, jpjm1
249            DO ji = 1, jpim1
250               zabe2 = e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)
251
252               zmkf  = 1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   &
253                              + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. )
254
255               zcof2 = -e1f(ji,jj) * zmkf   &
256                     * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
257
258               zjuf(ji,jj) = fmask(ji,jj,jk) *   &
259                           (  zabe2 * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) )   &
260                            + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
261                                       +zdk1u(ji,jj+1) + zdku (ji,jj) )  )
262            END DO
263         END DO
264
265         !                                 |   t   |
266         ! I.3 Horizontal fluxes on V      |       |
267         ! ------------------------===     f---v---f
268         !                                 |       |
269         ! i-flux at f-point               |   t   |
270         DO jj = 1, jpjm1
271            DO ji = 1, jpim1
272               zabe1 = e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
273
274               zmkf  = 1./MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
275                              + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. )
276
277               zcof1 = -e2f(ji,jj) * zmkf   &
278                     * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
279
280               zivf(ji,jj) = fmask(ji,jj,jk) *   &
281                           (  zabe1 * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) )   &
282                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj)     &
283                                       +zdk1u(ji,jj) + zdku (ji+1,jj) )  )
284            END DO
285         END DO
286
287         ! j-flux at t-point
288         DO jj = 2, jpj
289            DO ji = 1, jpim1
290               zabe2 = e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)
291
292               zmkt  = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
293                              + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
294
295               zcof2 = -e1t(ji,jj) * zmkt   &
296                     * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
297
298               zjvt(ji,jj) = tmask(ji,jj,jk) *   &
299                           (  zabe2 * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) )   &
300                            + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj)     &
301                                       +zdk1u(ji,jj-1) + zdku (ji,jj) )  )
302            END DO
303         END DO
304
305
306         ! I.4 Second derivative (divergence) (not divided by the volume)
307         ! ---------------------
308
309         DO jj = 2, jpjm1
310            DO ji = 2, jpim1
311               plu(ji,jj,jk) = ziut (ji+1,jj) - ziut (ji,jj  )   &
312                             + zjuf (ji  ,jj) - zjuf (ji,jj-1)
313               plv(ji,jj,jk) = zivf (ji,jj  ) - zivf (ji-1,jj)   &
314                             + zjvt (ji,jj+1) - zjvt (ji,jj  ) 
315            END DO
316         END DO
317
318         !                                             ! ===============
319      END DO                                           !   End of slab
320      !                                                ! ===============
321
322      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
323
324      !                             ! ************ !   ! ===============
325      DO jj = 2, jpjm1              !  Second step !   ! Horizontal slab
326         !                          ! ************ !   ! ===============
327
328         ! II.1 horizontal (pu,pv) gradients
329         ! ---------------------------------
330
331         DO jk = 1, jpk
332            DO ji = 2, jpi
333               ! i-gradient of u at jj
334               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( pu(ji,jj  ,jk) - pu(ji-1,jj  ,jk) )
335               ! j-gradient of u and v at jj
336               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( pu(ji,jj+1,jk) - pu(ji  ,jj  ,jk) )
337               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( pv(ji,jj  ,jk) - pv(ji  ,jj-1,jk) )
338               ! j-gradient of u and v at jj+1
339               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( pu(ji,jj  ,jk) - pu(ji  ,jj-1,jk) )
340               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pv(ji,jj+1,jk) - pv(ji  ,jj  ,jk) )
341            END DO
342         END DO
343         DO jk = 1, jpk
344            DO ji = 1, jpim1
345               ! i-gradient of v at jj
346               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( pv(ji+1,jj,jk) - pv(ji  ,jj  ,jk) )
347            END DO
348         END DO
349
350
351         ! II.2 Vertical fluxes
352         ! --------------------
353
354         ! Surface and bottom vertical fluxes set to zero
355
356         zfuw(:, 1 ) = 0.e0
357         zfvw(:, 1 ) = 0.e0
358         zfuw(:,jpk) = 0.e0
359         zfvw(:,jpk) = 0.e0
360
361         ! interior (2=<jk=<jpk-1) on pu field
362
363         DO jk = 2, jpkm1
364            DO ji = 2, jpim1
365               ! i- and j-slopes at uw-point
366               zuwslpi = 0.5 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
367               zuwslpj = 0.5 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
368               ! coef. for the vertical dirative
369               zcoef0 = e1u(ji,jj) * e2u(ji,jj) / fse3u(ji,jj,jk)   &
370                      * ( zuwslpi * zuwslpi + zuwslpj * zuwslpj )
371               ! weights for the i-k, j-k averaging at t- and f-points, resp.
372               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)   &
373                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ), 1. )
374               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1)   &
375                             + fmask(ji,jj-1,jk  )+fmask(ji,jj,jk  ), 1. )
376               ! coef. for the horitontal derivative
377               zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi
378               zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj
379               ! vertical flux on u field
380               zfuw(ji,jk) = umask(ji,jj,jk) *   &
381                           (  zcoef0 * ( pu  (ji,jj,jk-1) - pu  (ji,jj,jk) )   &
382                            + zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1)     &
383                                        +zdiu (ji,jk  ) + zdiu (ji+1,jk  ) )   &
384                            + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji  ,jk-1)     &
385                                        +zdj1u(ji,jk  ) + zdju (ji  ,jk  ) ) )
386            END DO
387         END DO
388
389         ! interior (2=<jk=<jpk-1) on pv field
390
391         DO jk = 2, jpkm1
392            DO ji = 2, jpim1
393               ! i- and j-slopes at vw-point
394               zvwslpi = 0.5 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
395               zvwslpj = 0.5 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
396               ! coef. for the vertical derivative
397               zcoef0 = e1v(ji,jj) * e2v(ji,jj) / fse3v(ji,jj,jk)   &
398                      * ( zvwslpi * zvwslpi + zvwslpj * zvwslpj )
399               ! weights for the i-k, j-k averaging at f- and t-points, resp.
400               zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)   &
401                             + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ), 1. )
402               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)   &
403                             + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ), 1. )
404               ! coef. for the horizontal derivatives
405               zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi
406               zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj
407               ! vertical flux on pv field
408               zfvw(ji,jk) = vmask(ji,jj,jk) *   &
409                           (  zcoef0 * ( pv  (ji,jj,jk-1) - pv  (ji,jj,jk) )   &
410                            + zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1)     &
411                                        +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   &
412                            + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     &
413                                        +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) )  )
414            END DO
415         END DO
416
417
418         ! II.3 Divergence of vertical fluxes added to the horizontal divergence
419         ! ---------------------------------------------------------------------
420
421         IF( kahm == 1 ) THEN
422            ! multiply the laplacian by the eddy viscosity coefficient
423            DO jk = 1, jpkm1
424               DO ji = 2, jpim1
425                  ! eddy coef. divided by the volume element
426                  zbur = fsahmu(ji,jj,jk) / ( e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) )
427                  zbvr = fsahmv(ji,jj,jk) / ( e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) )
428                  ! vertical divergence
429                  zuav = zfuw(ji,jk) - zfuw(ji,jk+1)
430                  zvav = zfvw(ji,jk) - zfvw(ji,jk+1)
431                  ! harmonic operator applied to (pu,pv) and multiply by ahm
432                  plu(ji,jj,jk) = ( plu(ji,jj,jk) + zuav ) * zbur
433                  plv(ji,jj,jk) = ( plv(ji,jj,jk) + zvav ) * zbvr
434               END DO
435            END DO
436         ELSEIF( kahm == 2 ) THEN
437            ! second call, no multiplication
438            DO jk = 1, jpkm1
439               DO ji = 2, jpim1
440                  ! inverse of the volume element
441                  zbur = 1. / ( e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) )
442                  zbvr = 1. / ( e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) )
443                  ! vertical divergence
444                  zuav = zfuw(ji,jk) - zfuw(ji,jk+1)
445                  zvav = zfvw(ji,jk) - zfvw(ji,jk+1)
446                  ! harmonic operator applied to (pu,pv)
447                  plu(ji,jj,jk) = ( plu(ji,jj,jk) + zuav ) * zbur
448                  plv(ji,jj,jk) = ( plv(ji,jj,jk) + zvav ) * zbvr
449               END DO
450            END DO
451         ELSE
452            IF(lwp)WRITE(numout,*) ' ldfguv: kahm= 1 or 2, here =', kahm
453            IF(lwp)WRITE(numout,*) '         We stop'
454            STOP 'ldfguv'
455         ENDIF
456         !                                             ! ===============
457      END DO                                           !   End of slab
458      !                                                ! ===============
459   END SUBROUTINE ldfguv
460
461#else
462   !!----------------------------------------------------------------------
463   !!   Dummy module :                         NO rotation of mixing tensor
464   !!----------------------------------------------------------------------
465CONTAINS
466   SUBROUTINE dyn_ldf_bilapg( kt )               ! Dummy routine
467      WRITE(*,*) kt
468   END SUBROUTINE dyn_ldf_bilapg
469#endif
470
471   !!======================================================================
472END MODULE dynldf_bilapg
Note: See TracBrowser for help on using the repository browser.