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

Last change on this file since 2690 was 2690, checked in by gm, 13 years ago

dynamic mem: #785 ; homogeneization of the coding style associated with dyn allocation

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