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 @ 247

Last change on this file since 247 was 247, checked in by opalod, 19 years ago

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