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

source: branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90 @ 3318

Last change on this file since 3318 was 3318, checked in by gm, 12 years ago

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