source: branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90 @ 3876

Last change on this file since 3876 was 3876, checked in by gm, 8 years ago

dev_r3858_CNRS3_Ediag: #927 phasing with 2011/dev_r3309_LOCEAN12_Ediag branche + mxl diag update

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