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

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

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