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/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90 @ 3634

Last change on this file since 3634 was 3634, checked in by acc, 11 years ago

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