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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90 @ 2633

Last change on this file since 2633 was 2633, checked in by trackstand2, 13 years ago

Renamed wrk_use => wrk_in_use and wrk_release => wrk_not_released
  • Property svn:keywords set to Id
File size: 22.5 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 in_out_manager  ! I/O manager
23   USE trdmod          ! ocean dynamics trends
24   USE trdmod_oce      ! ocean variables trends
25   USE ldfslp          ! iso-neutral slopes available
26   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
27   USE prtctl          ! Print control
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   dyn_ldf_bilapg       ! called by step.F90
33
34   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zfuw, zfvw , zdiu, zdiv   ! 2D workspace (ldfguv)
35   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zdju, zdj1u, zdjv, zdj1v  ! 2D workspace (ldfguv)
36
37   !! * Substitutions
38#  include "domzgr_substitute.h90"
39#  include "ldfdyn_substitute.h90"
40   !!----------------------------------------------------------------------
41   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
42   !! $Id$
43   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
44   !!----------------------------------------------------------------------
45
46CONTAINS
47
48   INTEGER FUNCTION dyn_ldf_bilapg_alloc()
49      !!----------------------------------------------------------------------
50      !!               ***  ROUTINE dyn_ldf_bilapg_alloc  ***
51      !!----------------------------------------------------------------------
52      !
53      ALLOCATE( zfuw(jpi,jpk) , zfvw (jpi,jpk) , zdiu(jpi,jpk) , zdiv (jpi,jpk) ,     &
54         &      zdju(jpi,jpk) , zdj1u(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_bilapg_alloc)
55         !
56      IF( dyn_ldf_bilapg_alloc /= 0 )   CALL ctl_warn('dyn_ldf_bilapg_alloc: failed to allocate arrays')
57      !
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      !!      'key_trddyn' defined: the trend is saved for diagnostics.
83      !!
84      !! ** Action  : - Update (ua,va) arrays with the before geopotential
85      !!                biharmonic mixing trend.
86      !!              - save the trend in (zwk3,zwk4) ('key_trddyn')
87      !!----------------------------------------------------------------------
88      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
89      USE wrk_nemo, ONLY:   zwk1 => wrk_3d_1 , zwk2 => wrk_3d_2   ! 3D workspace
90      USE oce     , ONLY:   zwk3 => ta       , zwk4 => sa         ! ta, sa used as 3D workspace   
91      !
92      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
93      !
94      INTEGER ::   ji, jj, jk                 ! dummy loop indices
95      !!----------------------------------------------------------------------
96
97      IF( wrk_in_use(3, 1,2) ) THEN
98         CALL ctl_stop('dyn_ldf_bilapg: requested workspace arrays unavailable.')   ;   RETURN
99      END IF
100
101      IF( kt == nit000 ) THEN
102         IF(lwp) WRITE(numout,*)
103         IF(lwp) WRITE(numout,*) 'dyn_ldf_bilapg : horizontal biharmonic operator in s-coordinate'
104         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
105         zwk1(:,:,:) = 0.e0   ;   zwk3(:,:,:) = 0.e0
106         zwk2(:,:,:) = 0.e0   ;   zwk4(:,:,:) = 0.e0
107         !                                      ! allocate dyn_ldf_bilapg arrays
108         IF( dyn_ldf_bilapg_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_ldf_bilapg: failed to allocate arrays')
109      ENDIF
110
111      ! Laplacian of (ub,vb) multiplied by ahm
112      ! -------------------------------------- 
113      CALL ldfguv( ub, vb, zwk1, zwk2, 1 )      ! rotated harmonic operator applied to (ub,vb)
114      !                                         ! and multiply by ahmu, ahmv (output in (zwk1,zwk2) )
115      CALL lbc_lnk( zwk1, 'U', -1. )   ;   CALL lbc_lnk( zwk2, 'V', -1. )     ! Lateral boundary conditions
116
117      ! Bilaplacian of (ub,vb)
118      ! ----------------------
119      CALL ldfguv( zwk1, zwk2, zwk3, zwk4, 2 )  ! rotated harmonic operator applied to (zwk1,zwk2)
120      !                                         ! (output in (zwk3,zwk4) )
121
122      ! Update the momentum trends
123      ! --------------------------
124      DO jj = 2, jpjm1               ! add the diffusive trend to the general momentum trends
125         DO jk = 1, jpkm1
126            DO ji = 2, jpim1
127               ua(ji,jj,jk) = ua(ji,jj,jk) + zwk3(ji,jj,jk)
128               va(ji,jj,jk) = va(ji,jj,jk) + zwk4(ji,jj,jk)
129            END DO
130         END DO
131      END DO
132      !
133      IF( wrk_not_released(3, 1,2) )   CALL ctl_stop('dyn_ldf_bilapg: failed to release workspace arrays')
134      !
135   END SUBROUTINE dyn_ldf_bilapg
136
137
138   SUBROUTINE ldfguv( pu, pv, plu, plv, kahm )
139      !!----------------------------------------------------------------------
140      !!                  ***  ROUTINE ldfguv  ***
141      !!     
142      !! ** Purpose :   Apply a geopotential harmonic operator to (pu,pv)
143      !!      (defined at u- and v-points) and multiply it by the eddy
144      !!      viscosity coefficient (if kahm=1).
145      !!
146      !! ** Method  :   The harmonic operator rotated along geopotential
147      !!      surfaces is applied to (pu,pv) using the slopes of geopotential
148      !!      surfaces computed in inildf routine. The result is provided in
149      !!      (plu,plv) arrays. It is computed in 2 stepv:
150      !!
151      !!      First step: horizontal part of the operator. It is computed on
152      !!      ==========  pu as follows (idem on pv)
153      !!      horizontal fluxes :
154      !!         zftu = e2u*e3u/e1u di[ pu ] - e2u*uslp dk[ mi(mk(pu)) ]
155      !!         zftv = e1v*e3v/e2v dj[ pu ] - e1v*vslp dk[ mj(mk(pu)) ]
156      !!      take the horizontal divergence of the fluxes (no divided by
157      !!      the volume element :
158      !!         plu  = di-1[ zftu ] +  dj-1[ zftv ]
159      !!
160      !!      Second step: vertical part of the operator. It is computed on
161      !!      ===========  pu as follows (idem on pv)
162      !!      vertical fluxes :
163      !!         zftw = e1t*e2t/e3w * (wslpi^2+wslpj^2)  dk-1[ pu ]
164      !!              -     e2t     *       wslpi        di[ mi(mk(pu)) ]
165      !!              -     e1t     *       wslpj        dj[ mj(mk(pu)) ]
166      !!      take the vertical divergence of the fluxes add it to the hori-
167      !!      zontal component, divide the result by the volume element and
168      !!      if kahm=1, multiply by the eddy diffusivity coefficient:
169      !!         plu  = aht / (e1t*e2t*e3t) { plu + dk[ zftw ] }
170      !!      else:
171      !!         plu  =  1  / (e1t*e2t*e3t) { plu + dk[ zftw ] }
172      !!
173      !! ** Action :
174      !!        plu, plv        : partial harmonic operator applied to
175      !!                          pu and pv (all the components except
176      !!                          second order vertical derivative term)
177      !!      'key_trddyn' defined: the trend is saved for diagnostics.
178      !!----------------------------------------------------------------------
179      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
180      USE wrk_nemo, ONLY: ziut => wrk_2d_1, zjuf => wrk_2d_2, zjvt => wrk_2d_3
181      USE wrk_nemo, ONLY: zivf => wrk_2d_4, zdku => wrk_2d_5, zdk1u => wrk_2d_6
182      USE wrk_nemo, ONLY: zdkv => wrk_2d_7, zdk1v => wrk_2d_8
183      !!
184      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu , pv    ! 1st call: before horizontal velocity
185      !                                                               ! 2nd call: ahm x these fields
186      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   plu, plv   ! partial harmonic operator applied to
187      !                                                               ! pu and pv (all the components except
188      !                                                               ! second order vertical derivative term)
189      INTEGER                         , INTENT(in   ) ::   kahm       ! =1 1st call ; =2 2nd call
190      !
191      INTEGER  ::   ji, jj, jk   ! dummy loop indices
192      REAL(wp) ::   zabe1 , zabe2 , zcof1 , zcof2        ! local scalar
193      REAL(wp) ::   zcoef0, zcoef3, zcoef4               !   -      -
194      REAL(wp) ::   zbur, zbvr, zmkt, zmkf, zuav, zvav   !   -      -
195      REAL(wp) ::   zuwslpi, zuwslpj, zvwslpi, zvwslpj   !   -      -
196      !!----------------------------------------------------------------------
197
198      IF( wrk_in_use(2, 1,2,3,4,5,6,7,8) ) THEN
199         CALL ctl_stop('dyn:ldfguv : requested workspace arrays unavailable.')   ;   RETURN
200      END IF
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
417         IF( kahm == 1 ) THEN
418            ! multiply the laplacian by the eddy viscosity coefficient
419            DO jk = 1, jpkm1
420               DO ji = 2, jpim1
421                  ! eddy coef. divided by the volume element
422                  zbur = fsahmu(ji,jj,jk) / ( e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) )
423                  zbvr = fsahmv(ji,jj,jk) / ( e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) )
424                  ! vertical divergence
425                  zuav = zfuw(ji,jk) - zfuw(ji,jk+1)
426                  zvav = zfvw(ji,jk) - zfvw(ji,jk+1)
427                  ! harmonic operator applied to (pu,pv) and multiply by ahm
428                  plu(ji,jj,jk) = ( plu(ji,jj,jk) + zuav ) * zbur
429                  plv(ji,jj,jk) = ( plv(ji,jj,jk) + zvav ) * zbvr
430               END DO
431            END DO
432         ELSEIF( kahm == 2 ) THEN
433            ! second call, no multiplication
434            DO jk = 1, jpkm1
435               DO ji = 2, jpim1
436                  ! inverse of the volume element
437                  zbur = 1. / ( e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) )
438                  zbvr = 1. / ( e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) )
439                  ! vertical divergence
440                  zuav = zfuw(ji,jk) - zfuw(ji,jk+1)
441                  zvav = zfvw(ji,jk) - zfvw(ji,jk+1)
442                  ! harmonic operator applied to (pu,pv)
443                  plu(ji,jj,jk) = ( plu(ji,jj,jk) + zuav ) * zbur
444                  plv(ji,jj,jk) = ( plv(ji,jj,jk) + zvav ) * zbvr
445               END DO
446            END DO
447         ELSE
448            IF(lwp)WRITE(numout,*) ' ldfguv: kahm= 1 or 2, here =', kahm
449            IF(lwp)WRITE(numout,*) '         We stop'
450            STOP 'ldfguv'
451         ENDIF
452         !                                             ! ===============
453      END DO                                           !   End of slab
454      !                                                ! ===============
455
456      IF(wrk_not_released(2, 1,2,3,4,5,6,7,8))THEN
457         CALL ctl_stop('dyn:ldfguv : failed to release workspace arrays.')
458      END IF
459      !
460   END SUBROUTINE ldfguv
461
462#else
463   !!----------------------------------------------------------------------
464   !!   Dummy module :                         NO rotation of mixing tensor
465   !!----------------------------------------------------------------------
466CONTAINS
467   SUBROUTINE dyn_ldf_bilapg( kt )               ! Dummy routine
468      INTEGER, INTENT(in) :: kt
469      WRITE(*,*) 'dyn_ldf_bilapg: You should not have seen this print! error?', kt
470   END SUBROUTINE dyn_ldf_bilapg
471#endif
472
473   !!======================================================================
474END MODULE dynldf_bilapg
Note: See TracBrowser for help on using the repository browser.