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

source: branches/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90 @ 2789

Last change on this file since 2789 was 2789, checked in by cetlod, 13 years ago

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