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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90 @ 2281

Last change on this file since 2281 was 2281, checked in by smasson, 14 years ago

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