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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90 @ 4438

Last change on this file since 4438 was 4432, checked in by trackstand2, 10 years ago

mbkmax changes in dynldf_bilapg
  • Property svn:keywords set to Id
File size: 40.2 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!FTRANS zfuw zfvw zdiu zdiv :I :z
36   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zfuw, zfvw , zdiu, zdiv   ! 2D workspace (ldfguv)
37!FTRANS zdju zdj1u zdjv zdj1v :I :z
38   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zdju, zdj1u, zdjv, zdj1v  ! 2D workspace (ldfguv)
39
40   !! * Control permutation of array indices
41#  include "oce_ftrans.h90"
42#  include "dom_oce_ftrans.h90"
43#  include "ldfdyn_oce_ftrans.h90"
44#  include "zdf_oce_ftrans.h90"
45#  include "ldfslp_ftrans.h90"
46
47   !! * Substitutions
48#  include "domzgr_substitute.h90"
49#  include "ldfdyn_substitute.h90"
50   !!----------------------------------------------------------------------
51   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
52   !! $Id$
53   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
54   !!----------------------------------------------------------------------
55CONTAINS
56
57   INTEGER FUNCTION dyn_ldf_bilapg_alloc()
58      !!----------------------------------------------------------------------
59      !!               ***  ROUTINE dyn_ldf_bilapg_alloc  ***
60      !!----------------------------------------------------------------------
61      ALLOCATE( zfuw(jpi,jpk) , zfvw (jpi,jpk) , zdiu(jpi,jpk) , zdiv (jpi,jpk) ,     &
62         &      zdju(jpi,jpk) , zdj1u(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_bilapg_alloc )
63         !
64      IF( dyn_ldf_bilapg_alloc /= 0 )   CALL ctl_warn('dyn_ldf_bilapg_alloc: failed to allocate arrays')
65   END FUNCTION dyn_ldf_bilapg_alloc
66
67
68   SUBROUTINE dyn_ldf_bilapg( kt )
69      !!----------------------------------------------------------------------
70      !!                   ***  ROUTINE dyn_ldf_bilapg  ***
71      !!                     
72      !! ** Purpose :   Compute the before trend of the horizontal momentum
73      !!      diffusion and add it to the general trend of momentum equation.
74      !!
75      !! ** Method  :   The lateral momentum diffusive trends is provided by a
76      !!      a 4th order operator rotated along geopotential surfaces. It is
77      !!      computed using before fields (forward in time) and geopotential
78      !!      slopes computed in routine inildf.
79      !!         -1- compute the geopotential harmonic operator applied to
80      !!      (ub,vb) and multiply it by the eddy diffusivity coefficient
81      !!      (done by a call to ldfgpu and ldfgpv routines) The result is in
82      !!      (zwk1,zwk2) arrays. Applied the domain lateral boundary conditions
83      !!      by call to lbc_lnk.
84      !!         -2- applied to (zwk1,zwk2) the geopotential harmonic operator
85      !!      by a second call to ldfgpu and ldfgpv routines respectively. The
86      !!      result is in (zwk3,zwk4) arrays.
87      !!         -3- Add this trend to the general trend (ta,sa):
88      !!            (ua,va) = (ua,va) + (zwk3,zwk4)
89      !!      'key_trddyn' defined: the trend is saved for diagnostics.
90      !!
91      !! ** Action  : - Update (ua,va) arrays with the before geopotential
92      !!                biharmonic mixing trend.
93      !!              - save the trend in (zwk3,zwk4) ('key_trddyn')
94      !!----------------------------------------------------------------------
95      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
96      USE wrk_nemo, ONLY:   zwk1 => wrk_3d_3 , zwk2 => wrk_3d_4   ! 3D workspace
97      USE oce     , ONLY:   zwk3 => ta       , zwk4 => sa         ! ta, sa used as 3D workspace   
98      !! DCSE_NEMO: need additional directives for renamed module variables
99!FTRANS zwk1 :I :I :z
100!FTRANS zwk2 :I :I :z
101!FTRANS zwk3 :I :I :z
102!FTRANS zwk4 :I :I :z
103      !
104      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
105      !
106      INTEGER ::   ji, jj, jk                 ! dummy loop indices
107      !!----------------------------------------------------------------------
108
109      IF( wrk_in_use(3, 3,4) ) THEN
110         CALL ctl_stop('dyn_ldf_bilapg: requested workspace arrays unavailable')   ;   RETURN
111      ENDIF
112
113      IF( kt == nit000 ) THEN
114         IF(lwp) WRITE(numout,*)
115         IF(lwp) WRITE(numout,*) 'dyn_ldf_bilapg : horizontal biharmonic operator in s-coordinate'
116         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
117         zwk1(:,:,:) = 0.e0_wp   ;   zwk3(:,:,:) = 0.e0_wp
118         zwk2(:,:,:) = 0.e0_wp   ;   zwk4(:,:,:) = 0.e0_wp
119         !                                      ! allocate dyn_ldf_bilapg arrays
120         IF( dyn_ldf_bilapg_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_ldf_bilapg: failed to allocate arrays')
121      ENDIF
122
123      ! Laplacian of (ub,vb) multiplied by ahm
124      ! -------------------------------------- 
125      CALL ldfguv( ub, vb, zwk1, zwk2, 1 )      ! rotated harmonic operator applied to (ub,vb)
126      !                                         ! and multiply by ahmu, ahmv (output in (zwk1,zwk2) )
127      CALL lbc_lnk( zwk1, 'U', -1. )   ;   CALL lbc_lnk( zwk2, 'V', -1. )     ! Lateral boundary conditions
128
129      ! Bilaplacian of (ub,vb)
130      ! ----------------------
131      CALL ldfguv( zwk1, zwk2, zwk3, zwk4, 2 )  ! rotated harmonic operator applied to (zwk1,zwk2)
132      !                                         ! (output in (zwk3,zwk4) )
133
134      ! Update the momentum trends
135      ! --------------------------
136#if defined key_z_first
137      DO jj = 2, jpjm1               ! add the diffusive trend to the general momentum trends
138         DO ji = 2, jpim1
139            DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
140#else
141      DO jj = 2, jpjm1               ! add the diffusive trend to the general momentum trends
142         DO jk = 1, jpkm1
143            DO ji = 2, jpim1
144#endif
145               ua(ji,jj,jk) = ua(ji,jj,jk) + zwk3(ji,jj,jk)
146               va(ji,jj,jk) = va(ji,jj,jk) + zwk4(ji,jj,jk)
147            END DO
148         END DO
149      END DO
150      !
151      IF( wrk_not_released(3, 3,4) )   CALL ctl_stop('dyn_ldf_bilapg: failed to release workspace arrays')
152      !
153   END SUBROUTINE dyn_ldf_bilapg
154
155
156   SUBROUTINE ldfguv( pu, pv, plu, plv, kahm )
157      !!----------------------------------------------------------------------
158      !!                  ***  ROUTINE ldfguv  ***
159      !!     
160      !! ** Purpose :   Apply a geopotential harmonic operator to (pu,pv)
161      !!      (defined at u- and v-points) and multiply it by the eddy
162      !!      viscosity coefficient (if kahm=1).
163      !!
164      !! ** Method  :   The harmonic operator rotated along geopotential
165      !!      surfaces is applied to (pu,pv) using the slopes of geopotential
166      !!      surfaces computed in inildf routine. The result is provided in
167      !!      (plu,plv) arrays. It is computed in 2 stepv:
168      !!
169      !!      First step: horizontal part of the operator. It is computed on
170      !!      ==========  pu as follows (idem on pv)
171      !!      horizontal fluxes :
172      !!         zftu = e2u*e3u/e1u di[ pu ] - e2u*uslp dk[ mi(mk(pu)) ]
173      !!         zftv = e1v*e3v/e2v dj[ pu ] - e1v*vslp dk[ mj(mk(pu)) ]
174      !!      take the horizontal divergence of the fluxes (no divided by
175      !!      the volume element :
176      !!         plu  = di-1[ zftu ] +  dj-1[ zftv ]
177      !!
178      !!      Second step: vertical part of the operator. It is computed on
179      !!      ===========  pu as follows (idem on pv)
180      !!      vertical fluxes :
181      !!         zftw = e1t*e2t/e3w * (wslpi^2+wslpj^2)  dk-1[ pu ]
182      !!              -     e2t     *       wslpi        di[ mi(mk(pu)) ]
183      !!              -     e1t     *       wslpj        dj[ mj(mk(pu)) ]
184      !!      take the vertical divergence of the fluxes add it to the hori-
185      !!      zontal component, divide the result by the volume element and
186      !!      if kahm=1, multiply by the eddy diffusivity coefficient:
187      !!         plu  = aht / (e1t*e2t*e3t) { plu + dk[ zftw ] }
188      !!      else:
189      !!         plu  =  1  / (e1t*e2t*e3t) { plu + dk[ zftw ] }
190      !!
191      !! ** Action :
192      !!        plu, plv        : partial harmonic operator applied to
193      !!                          pu and pv (all the components except
194      !!                          second order vertical derivative term)
195      !!      'key_trddyn' defined: the trend is saved for diagnostics.
196      !!----------------------------------------------------------------------
197      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
198      USE wrk_nemo, ONLY:   ziut => wrk_2d_1 , zjuf  => wrk_2d_2 , zjvt  => wrk_2d_3
199      USE wrk_nemo, ONLY:   zivf => wrk_2d_4
200#if ! defined key_z_first
201      USE wrk_nemo, ONLY:   zdku => wrk_2d_5 , zdk1u => wrk_2d_6
202      USE wrk_nemo, ONLY:   zdkv => wrk_2d_7 , zdk1v => wrk_2d_8
203#endif
204      USE timing  , ONLY:   timing_start, timing_stop
205      !!
206!FTRANS pu :I :I :z
207!FTRANS pv :I :I :z
208!FTRANS plu :I :I :z
209!FTRANS plv :I :I :z
210!! DCSE_NEMO: work around deficiency in ftrans
211!     REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu , pv    ! 1st call: before horizontal velocity
212      !                                                               ! 2nd call: ahm x these fields
213      REAL(wp), INTENT(in   ) ::   pu(jpi,jpj,jpkorig) , pv(jpi,jpj,jpkorig)
214!     REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   plu, plv   ! partial harmonic operator applied to
215      !                                                               ! pu and pv (all the components except
216      !                                                               ! second order vertical derivative term)
217      REAL(wp), INTENT(  out) ::   plu(jpi,jpj,jpkorig), plv(jpi,jpj,jpkorig) ! partial harmonic operator applied to
218      INTEGER                         , INTENT(in   ) ::   kahm       ! =1 1st call ; =2 2nd call
219      !
220      INTEGER  ::   ji, jj, jk   ! dummy loop indices
221      INTEGER  ::   jif, jjf     ! dummy loop indices over full domain
222      REAL(wp) ::   zabe1 , zabe2 , zcof1 , zcof2        ! local scalar
223      REAL(wp) ::   zcoef0, zcoef3, zcoef4               !   -      -
224      REAL(wp) ::   zbur, zbvr, zmkt, zmkf, zuav, zvav   !   -      -
225      REAL(wp) ::   zuwslpi, zuwslpj, zvwslpi, zvwslpj   !   -      -
226#if defined key_z_first
227      ! Can use scalars instead of work arrays when built with z-first
228      REAL(wp) ::   zdku, zdkv, zdk1u, zdk1v
229#endif
230      !!----------------------------------------------------------------------
231
232      CALL timing_start('ldfguv')
233
234      IF( wrk_in_use(2, 1,2,3,4,5,6,7,8) ) THEN
235         CALL ctl_stop('dyn:ldfguv: requested workspace arrays unavailable')   ;   RETURN
236      END IF
237      CALL timing_start('ldfguv_1st')
238
239#if defined key_z_first
240      !                               ! ********** !
241                                      ! First step !
242         !                            ! ********** !
243         DO jj = 2, jpjm1
244            DO ji = 2, jpim1
245
246               ! Treat jk = 1 separately as is special case
247               jk = 1
248
249               plu(ji,jj,jk) = &
250! ------------- ziut (ji+1, jj) -
251                           tmask(ji+1,jj,jk) *   &
252                           (  (e2t(ji+1,jj) * fse3t(ji+1,jj,jk) / e1t(ji+1,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) )   &
253!                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
254                            + (-e2t(ji+1,jj) / MAX(  umask(ji,jj,jk  )+umask(ji+1,jj,jk+1)   &
255                                 + umask(ji,jj,jk+1)+umask(ji+1,jj,jk  ), 1._wp )   &
256                     * 0.5  * ( uslp(ji,jj,jk) + uslp(ji+1,jj,jk) ) ) * &
257( ( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) + &
258  ( pu(ji  ,jj,jk) - pu(ji  ,jj,jk+1) ) * umask(ji  ,jj,jk+1) + &
259!                                       +zdk1u(ji,jj) + zdku (ji-1,jj) )  )
260  ( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) + &
261  ( pu(ji  ,jj,jk) - pu(ji  ,jj,jk+1) ) * umask(ji  ,jj,jk+1) &
262)  ) - &
263! ------------- ziut (ji,jj  ) +
264                           tmask(ji,jj,jk) *   &
265                           (  (e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) )   &
266!                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
267                            + (-e2t(ji,jj) / MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
268                                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1._wp )   &
269                     * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) ) * &
270( ( pu(ji  ,jj,jk) - pu(ji  ,jj,jk+1) ) * umask(ji,jj,jk+1) + &
271  ( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) + &
272!                                       +zdk1u(ji,jj) + zdku (ji-1,jj) )  )
273  ( pu(ji  ,jj,jk) - pu(ji  ,jj,jk+1) ) * umask(ji,jj,jk+1) + &
274  ( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) &
275)  ) + &
276! ------------- zjuf (ji  ,jj) -
277                           fmask(ji,jj,jk) *   &
278                           (  (e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) ) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) )   &
279!                            + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
280                             + (-e1f(ji,jj) /MAX(umask(ji,jj+1,jk)+umask(ji,jj,jk+1)   &
281                                      + umask(ji,jj+1,jk+1)+umask(ji,jj,jk), 1. )   &
282                     * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) ) * &
283( &
284             ( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1)  &
285           + ( pu(ji,jj  ,jk) - pu(ji,jj  ,jk+1) ) * umask(ji,jj  ,jk+1)  &
286!                                       +zdk1u(ji,jj+1) + zdku (ji,jj) )  )
287           + ( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) &
288           + ( pu(ji,jj  ,jk) - pu(ji,jj  ,jk+1) ) * umask(ji,jj,jk+1)   &
289)  ) - &
290! ------------- zjuf (ji,jj-1)
291                           fmask(ji,jj-1,jk) *   &
292                           (  (e1f(ji,jj-1) * fse3f(ji,jj-1,jk) / e2f(ji,jj-1) ) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) )   &
293!                            + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
294                             + (-e1f(ji,jj-1) /MAX(umask(ji,jj,jk)+umask(ji,jj-1,jk+1)   &
295                                      + umask(ji,jj,jk+1)+umask(ji,jj-1,jk), 1. )   &
296                     * 0.5  * ( vslp(ji+1,jj-1,jk) + vslp(ji,jj-1,jk) ) ) * &
297( &
298             ( pu(ji,jj,jk) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1)  &
299           + ( pu(ji,jj-1  ,jk) - pu(ji,jj-1  ,jk+1) ) * umask(ji,jj-1  ,jk+1)  &
300!                                       +zdk1u(ji,jj+1) + zdku (ji,jj) )  )
301           + ( pu(ji,jj,jk) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1) &
302           + ( pu(ji,jj-1  ,jk) - pu(ji,jj-1  ,jk+1) ) * umask(ji,jj-1,jk+1)   &
303)  )
304
305
306                  plv(ji,jj,jk) = &
307! ------------- zivf (ji,jj  ) -
308fmask(ji,jj,jk) *   &
309                           (  (e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) )   &
310!                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj)     &
311                             + ((-e2f(ji,jj) / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
312                              + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. ))   &
313                     * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * ( &
314  ( pu(ji  ,jj,jk) - pu(ji  ,jj,jk+1) ) * umask(ji  ,jj,jk+1) + &
315  ( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) +    &
316!                                       +zdk1u(ji,jj) + zdku (ji+1,jj) )  )
317  ( pu(ji  ,jj,jk) - pu(ji  ,jj,jk+1) ) * umask(ji  ,jj,jk+1) +  &
318  ( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1)    &
319) ) - &
320! ------------- zivf (ji-1,jj) +
321fmask(ji-1,jj,jk) *   &
322                           (  (e2f(ji-1,jj) * fse3f(ji-1,jj,jk) / e1f(ji-1,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) )   &
323!                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj)     &
324                             + ((-e2f(ji-1,jj) / MAX(  vmask(ji,jj,jk  )+vmask(ji-1,jj,jk+1)   &
325                              + vmask(ji,jj,jk+1)+vmask(ji-1,jj,jk  ), 1. ))   &
326                     * 0.5 * ( uslp(ji-1,jj+1,jk) + uslp(ji-1,jj,jk) )) * ( &
327  ( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) + &
328  ( pu(ji  ,jj,jk) - pu(ji  ,jj,jk+1) ) * umask(ji  ,jj,jk+1) +    &
329!                                       +zdk1u(ji,jj) + zdku (ji+1,jj) )  )
330  ( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) +  &
331  ( pu(ji  ,jj,jk) - pu(ji  ,jj,jk+1) ) * umask(ji  ,jj,jk+1)    &
332) ) + &
333! ------------- zjvt (ji,jj+1) -
334tmask(ji,jj+1,jk) *   &
335                           (  (e1t(ji,jj+1) * fse3t(ji,jj+1,jk) / e2t(ji,jj+1)) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) )   &
336!                            + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj)     &
337                            + ((-e1t(ji,jj+1)/MAX(  vmask(ji,jj,jk  )+vmask(ji,jj+1,jk+1)   &
338                              + vmask(ji,jj,jk+1)+vmask(ji,jj+1,jk  ), 1. ) )   &
339                     * 0.5 * ( vslp(ji,jj,jk) + vslp(ji,jj+1,jk) )) * ( &
340  ( pu(ji,jj  ,jk) - pu(ji,jj  ,jk+1) ) * umask(ji,jj  ,jk+1) + &
341  ( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1)   +  &
342!                                       +zdk1u(ji,jj-1) + zdku (ji,jj) )  )
343  ( pu(ji,jj  ,jk) - pu(ji,jj  ,jk+1) ) * umask(ji,jj  ,jk+1) + &
344  ( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1)  & 
345)  ) - &
346! ------------- zjvt (ji,jj  )
347tmask(ji,jj,jk) *   &
348                           (  (e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) )   &
349!                            + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj)     &
350                            + ((-e1t(ji,jj)/MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
351                              + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. ) )   &
352                     * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * ( &
353  ( pu(ji,jj-1,jk) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + &
354  ( pu(ji,jj  ,jk) - pu(ji,jj  ,jk+1) ) * umask(ji,jj,jk+1)   +  &
355!                                       +zdk1u(ji,jj-1) + zdku (ji,jj) )  )
356  ( pu(ji,jj-1,jk) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + &
357  ( pu(ji,jj  ,jk) - pu(ji,jj  ,jk+1) ) * umask(ji,jj,jk+1)  &
358)  )
359
360
361               DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1
362
363!               plu(ji,jj,jk) = ziut (ji+1,jj) - &
364!                               ziut (ji,jj  ) + &
365!                               zjuf (ji  ,jj) - &
366!                               zjuf (ji,jj-1)
367                  plu(ji,jj,jk) = &
368! ------------- ziut (ji+1, jj  ) -
369                  tmask(ji+1,jj,jk) *   &
370                           (  (e2t(ji+1,jj) * fse3t(ji+1,jj,jk) / e1t(ji+1,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) )   &
371!                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
372                            + (-e2t(ji+1,jj) / MAX(umask(ji,jj,jk)+umask(ji+1,jj,jk+1)   &
373                                 + umask(ji,jj,jk+1)+umask(ji+1,jj,jk), 1._wp )   &
374                     * 0.5  * ( uslp(ji,jj,jk) + uslp(ji+1,jj,jk) ) ) * &
375( ( pu(ji+1,jj,jk-1) - pu(ji+1,jj,jk  ) ) * umask(ji+1,jj,jk) + &
376  ( pu(ji  ,jj,jk  ) - pu(ji  ,jj,jk+1) ) * umask(ji  ,jj,jk+1) +     &
377!                                       +zdk1u(ji,jj) + zdku (ji-1,jj) )  )
378  ( pu(ji+1,jj,jk  ) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) + &
379  ( pu(ji  ,jj,jk-1) - pu(ji  ,jj,jk  ) ) * umask(ji  ,jj,jk) &
380)  ) - &
381! ------------- ziut (ji  , jj  ) +
382                  tmask(ji,jj,jk) *   &
383                           (  (e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) )   &
384!                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
385                            + (-e2t(ji,jj) / MAX(umask(ji-1,jj,jk)+umask(ji,jj,jk+1)   &
386                                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk), 1._wp )   &
387                     * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) ) * &
388( ( pu(ji  ,jj,jk-1) - pu(ji  ,jj,jk  ) ) * umask(ji,jj,jk) + &
389  ( pu(ji-1,jj,jk  ) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) +     &
390!                                       +zdk1u(ji,jj) + zdku (ji-1,jj) )  )
391  ( pu(ji  ,jj,jk  ) - pu(ji  ,jj,jk+1) ) * umask(ji,jj,jk+1) + &
392  ( pu(ji-1,jj,jk-1) - pu(ji-1,jj,jk  ) ) * umask(ji-1,jj,jk) &
393)  ) + &
394! ------------- zjuf (ji  , jj  ) -
395fmask(ji,jj,jk) *   &
396                           (  (e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) ) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) )   &
397!                            + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
398                             + (-e1f(ji,jj) /MAX(umask(ji,jj+1,jk)+umask(ji,jj,jk+1)   &
399                                      + umask(ji,jj+1,jk+1)+umask(ji,jj,jk), 1. )   &
400                     * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) ) * &
401(  &
402  (pu(ji,jj+1,jk-1) - pu(ji,jj+1,jk  ) ) * umask(ji,jj+1,jk) +   &
403  (pu(ji,jj  ,jk  ) - pu(ji,jj  ,jk+1) ) * umask(ji,jj,jk+1) +   &
404!                                       +zdk1u(ji,jj+1) + zdku (ji,jj) )  )
405  (pu(ji,jj+1,jk  ) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) + &
406  (pu(ji,jj  ,jk-1) - pu(ji,jj  ,jk  ) ) * umask(ji,jj,jk) &
407 )  ) - &
408! ------------- zjuf (ji  , jj-1)
409fmask(ji,jj-1,jk) *   &
410                           (  (e1f(ji,jj-1) * fse3f(ji,jj-1,jk) / e2f(ji,jj-1) ) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) )   &
411!                            + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
412                             + (-e1f(ji,jj-1) /MAX(umask(ji,jj,jk)+umask(ji,jj-1,jk+1)   &
413                                      + umask(ji,jj,jk+1)+umask(ji,jj-1,jk), 1. )   &
414                     * 0.5  * ( vslp(ji+1,jj-1,jk) + vslp(ji,jj-1,jk) ) ) * &
415(  &
416  (pu(ji,jj,jk-1) - pu(ji,jj,jk  ) ) * umask(ji,jj,jk) +   &
417  (pu(ji,jj-1,jk  ) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) +   &
418!                                       +zdk1u(ji,jj+1) + zdku (ji,jj) )  )
419  (pu(ji,jj,jk  ) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1) + &
420  (pu(ji,jj-1,jk-1) - pu(ji,jj-1  ,jk  ) ) * umask(ji,jj-1,jk) &
421 )  )
422
423
424!               plv(ji,jj,jk) = zivf (ji,jj  ) - &
425!                               zivf (ji-1,jj) + &
426!                               zjvt (ji,jj+1) - &
427!                               zjvt (ji,jj  )
428                  plv(ji,jj,jk) = &
429! ------------- zivf (ji,jj  ) -
430fmask(ji,jj,jk) *   &
431                           (  (e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) )   &
432!                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj)     &
433+ ((-e2f(ji,jj) / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
434                              + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. ))   &
435                     * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * ( &
436  ( pu(ji  ,jj,jk-1) - pu(ji  ,jj,jk  ) ) * umask(ji  ,jj,jk  ) + &
437  ( pu(ji+1,jj,jk  ) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) +    &
438!                                       +zdk1u(ji,jj) + zdku (ji+1,jj) )  )
439  ( pu(ji  ,jj,jk  ) - pu(ji  ,jj,jk+1) ) * umask(ji  ,jj,jk+1) + &
440  ( pu(ji+1,jj,jk-1) - pu(ji+1,jj,jk  ) ) * umask(ji+1,jj,jk  )   &
441)  ) - &
442! ------------- zivf (ji-1,jj) +
443fmask(ji-1,jj,jk) *   &
444                           (  (e2f(ji-1,jj) * fse3f(ji-1,jj,jk) / e1f(ji-1,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) )   &
445!                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj)     &
446+ ((-e2f(ji-1,jj) / MAX(  vmask(ji,jj,jk  )+vmask(ji-1,jj,jk+1)   &
447                              + vmask(ji,jj,jk+1)+vmask(ji-1,jj,jk  ), 1. ))   &
448                     * 0.5 * ( uslp(ji-1,jj+1,jk) + uslp(ji-1,jj,jk) )) * ( &
449  ( pu(ji-1 ,jj,jk-1) - pu(ji-1  ,jj,jk  ) ) * umask(ji-1  ,jj,jk  ) + &
450  ( pu(ji,jj,jk  ) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1) +    &
451!                                       +zdk1u(ji,jj) + zdku (ji+1,jj) )  )
452  ( pu(ji-1  ,jj,jk  ) - pu(ji-1  ,jj,jk+1) ) * umask(ji-1  ,jj,jk+1) + &
453  ( pu(ji,jj,jk-1) - pu(ji,jj,jk  ) ) * umask(ji,jj,jk  )  &
454)  ) + &
455! ------------- zjvt (ji,jj+1) -
456tmask(ji,jj+1,jk) *   &
457                           (  (e1t(ji,jj+1) * fse3t(ji,jj+1,jk) / e2t(ji,jj+1)) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) )   &
458!                            + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj)     &
459                            + ((-e1t(ji,jj+1)/MAX(  vmask(ji,jj,jk  )+vmask(ji,jj+1,jk+1)   &
460                              + vmask(ji,jj,jk+1)+vmask(ji,jj+1,jk  ), 1. ) )   &
461                     * 0.5 * ( vslp(ji,jj,jk) + vslp(ji,jj+1,jk) )) * &
462( &
463  ( pu(ji,jj,jk-1) - pu(ji,jj,jk  ) ) * umask(ji,jj,jk) + &
464  ( pu(ji,jj+1  ,jk  ) - pu(ji,jj+1  ,jk+1) ) * umask(ji,jj+1,jk+1)  +    &
465!                                       +zdk1u(ji,jj-1) + zdku (ji,jj) )  )
466  ( pu(ji,jj,jk  ) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1) + &
467  ( pu(ji,jj+1  ,jk-1) - pu(ji,jj+1  ,jk  ) ) * umask(ji,jj+1,jk) &
468)  ) - &
469! ------------- zjvt (ji,jj  )
470tmask(ji,jj,jk) *   &
471                           (  (e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) )   &
472!                            + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj)     &
473                            + ((-e1t(ji,jj)/MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
474                              + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. ) )   &
475                     * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * &
476( &
477  ( pu(ji,jj-1,jk-1) - pu(ji,jj-1,jk  ) ) * umask(ji,jj-1,jk) + &
478  ( pu(ji,jj  ,jk  ) - pu(ji,jj  ,jk+1) ) * umask(ji,jj,jk+1)  +    &
479!                                       +zdk1u(ji,jj-1) + zdku (ji,jj) )  )
480  ( pu(ji,jj-1,jk  ) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + &
481  ( pu(ji,jj  ,jk-1) - pu(ji,jj  ,jk  ) ) * umask(ji,jj,jk) )  )
482
483            END DO
484         END DO
485
486         !                                       
487      END DO
488
489#else
490      !                               ! ********** !   ! ===============
491      DO jk = 1, jpkm1                ! First step !   ! Horizontal slab
492         !                            ! ********** !   ! ===============
493
494         ! I.1 Vertical gradient of pu and pv at level jk and jk+1
495         ! -------------------------------------------------------
496         ! surface boundary condition: zdku(jk=1)=zdku(jk=2)
497         !                             zdkv(jk=1)=zdkv(jk=2)
498
499         zdk1u(:,:) = ( pu(:,:,jk) - pu(:,:,jk+1) ) * umask(:,:,jk+1)
500         zdk1v(:,:) = ( pv(:,:,jk) - pv(:,:,jk+1) ) * vmask(:,:,jk+1)
501
502         IF( jk == 1 ) THEN
503            zdku(:,:) = zdk1u(:,:)
504            zdkv(:,:) = zdk1v(:,:)
505         ELSE
506            zdku(:,:) = ( pu(:,:,jk-1) - pu(:,:,jk) ) * umask(:,:,jk)
507            zdkv(:,:) = ( pv(:,:,jk-1) - pv(:,:,jk) ) * vmask(:,:,jk)
508         ENDIF
509
510         !                                -----f-----
511         ! I.2 Horizontal fluxes on U          |
512         ! ------------------------===     t   u   t
513         !                                     |
514         ! i-flux at t-point              -----f-----
515         DO jj = 1, jpjm1
516            DO ji = 2, jpi
517               zabe1 = e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)
518
519               zmkt  = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
520                              + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
521
522               zcof1 = -e2t(ji,jj) * zmkt   &
523                     * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
524
525               ziut(ji,jj) = tmask(ji,jj,jk) *   &
526                           (  zabe1 * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) )   &
527                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
528                                       +zdk1u(ji,jj) + zdku (ji-1,jj) )  )
529            END DO
530         END DO
531
532         ! j-flux at f-point
533         DO jj = 1, jpjm1
534            DO ji = 1, jpim1
535               zabe2 = e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)
536
537               zmkf  = 1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   &
538                              + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. )
539
540               zcof2 = -e1f(ji,jj) * zmkf   &
541                     * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
542
543               zjuf(ji,jj) = fmask(ji,jj,jk) *   &
544                           (  zabe2 * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) )   &
545                            + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
546                                       +zdk1u(ji,jj+1) + zdku (ji,jj) )  )
547            END DO
548         END DO
549
550         !                                 |   t   |
551         ! I.3 Horizontal fluxes on V      |       |
552         ! ------------------------===     f---v---f
553         !                                 |       |
554         ! i-flux at f-point               |   t   |
555         DO jj = 1, jpjm1
556            DO ji = 1, jpim1
557               zabe1 = e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
558
559               zmkf  = 1./MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
560                              + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. )
561
562               zcof1 = -e2f(ji,jj) * zmkf   &
563                     * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
564
565               zivf(ji,jj) = fmask(ji,jj,jk) *   &
566                           (  zabe1 * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) )   &
567                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj)     &
568                                       +zdk1u(ji,jj) + zdku (ji+1,jj) )  )
569            END DO
570         END DO
571
572         ! j-flux at t-point
573         DO jj = 2, jpj
574            DO ji = 1, jpim1
575               zabe2 = e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)
576
577               zmkt  = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
578                              + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
579
580               zcof2 = -e1t(ji,jj) * zmkt   &
581                     * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
582
583               zjvt(ji,jj) = tmask(ji,jj,jk) *   &
584                           (  zabe2 * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) )   &
585                            + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj)     &
586                                       +zdk1u(ji,jj-1) + zdku (ji,jj) )  )
587            END DO
588         END DO
589
590
591         ! I.4 Second derivative (divergence) (not divided by the volume)
592         ! ---------------------
593
594         DO jj = 2, jpjm1
595            DO ji = 2, jpim1
596               plu(ji,jj,jk) = ziut (ji+1,jj) - ziut (ji,jj  )   &
597                             + zjuf (ji  ,jj) - zjuf (ji,jj-1)
598               plv(ji,jj,jk) = zivf (ji,jj  ) - zivf (ji-1,jj)   &
599                             + zjvt (ji,jj+1) - zjvt (ji,jj  ) 
600            END DO
601         END DO
602
603         !                                             ! ===============
604      END DO                                           !   End of slab
605      !                                                ! ===============
606#endif
607      CALL timing_stop('ldfguv_1st','section')
608
609      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
610
611      CALL timing_start('ldfguv_2nd')
612      !                             ! ************ !   ! ===============
613      DO jj = 2, jpjm1              !  Second step !   ! Horizontal slab
614         !                          ! ************ !   ! ===============
615
616         ! II.1 horizontal (pu,pv) gradients
617         ! ---------------------------------
618
619#if defined key_z_first
620         DO ji = 2, jpi
621            DO jk = 1, mbkmax(ji,jj) ! jpk
622               ! i-gradient of u at jj
623               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( pu(ji,jj  ,jk) - pu(ji-1,jj  ,jk) )
624               ! j-gradient of u and v at jj
625               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( pu(ji,jj+1,jk) - pu(ji  ,jj  ,jk) )
626               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( pv(ji,jj  ,jk) - pv(ji  ,jj-1,jk) )
627               ! j-gradient of u and v at jj+1
628               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( pu(ji,jj  ,jk) - pu(ji  ,jj-1,jk) )
629               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pv(ji,jj+1,jk) - pv(ji  ,jj  ,jk) )
630            END DO
631         END DO
632         DO ji = 1, jpim1
633            DO jk = 1, mbkmax(ji,jj) ! jpk
634               ! i-gradient of v at jj
635               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( pv(ji+1,jj,jk) - pv(ji  ,jj  ,jk) )
636            END DO
637         END DO
638#else
639         DO jk = 1, jpk
640            DO ji = 2, jpi
641               ! i-gradient of u at jj
642               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( pu(ji,jj  ,jk) - pu(ji-1,jj  ,jk) )
643               ! j-gradient of u and v at jj
644               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( pu(ji,jj+1,jk) - pu(ji  ,jj  ,jk) )
645               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( pv(ji,jj  ,jk) - pv(ji  ,jj-1,jk) )
646               ! j-gradient of u and v at jj+1
647               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( pu(ji,jj  ,jk) - pu(ji  ,jj-1,jk) )
648               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pv(ji,jj+1,jk) - pv(ji  ,jj  ,jk) )
649            END DO
650         END DO
651         DO jk = 1, jpk
652            DO ji = 1, jpim1
653               ! i-gradient of v at jj
654               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( pv(ji+1,jj,jk) - pv(ji  ,jj  ,jk) )
655            END DO
656         END DO
657#endif
658
659         ! II.2 Vertical fluxes
660         ! --------------------
661
662         ! Surface and bottom vertical fluxes set to zero
663
664#if defined key_z_first
665         DO ji=1, jpi
666            zfuw(ji, 1               ) = 0.e0
667            zfvw(ji, 1               ) = 0.e0
668            zfuw(ji,mbkmax(ji,jj):jpk) = 0.e0
669            zfvw(ji,mbkmax(ji,jj):jpk) = 0.e0
670         END DO
671#else
672         zfuw(:, 1 ) = 0.e0
673         zfvw(:, 1 ) = 0.e0
674         zfuw(:,jpk) = 0.e0
675         zfvw(:,jpk) = 0.e0
676#endif
677         ! interior (2=<jk=<jpk-1) on pu field
678
679#if defined key_z_first
680         DO ji = 2, jpim1
681            DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1
682#else
683         DO jk = 2, jpkm1
684            DO ji = 2, jpim1
685#endif
686               ! i- and j-slopes at uw-point
687               zuwslpi = 0.5 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
688               zuwslpj = 0.5 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
689               ! coef. for the vertical dirative
690               zcoef0 = e1u(ji,jj) * e2u(ji,jj) / fse3u(ji,jj,jk)   &
691                      * ( zuwslpi * zuwslpi + zuwslpj * zuwslpj )
692               ! weights for the i-k, j-k averaging at t- and f-points, resp.
693               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)   &
694                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ), 1. )
695               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1)   &
696                             + fmask(ji,jj-1,jk  )+fmask(ji,jj,jk  ), 1. )
697               ! coef. for the horitontal derivative
698               zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi
699               zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj
700               ! vertical flux on u field
701               zfuw(ji,jk) = umask(ji,jj,jk) *   &
702                           (  zcoef0 * ( pu  (ji,jj,jk-1) - pu  (ji,jj,jk) )   &
703                            + zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1)     &
704                                        +zdiu (ji,jk  ) + zdiu (ji+1,jk  ) )   &
705                            + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji  ,jk-1)     &
706                                        +zdj1u(ji,jk  ) + zdju (ji  ,jk  ) ) )
707            END DO
708         END DO
709
710         ! interior (2=<jk=<jpk-1) on pv field
711
712#if defined key_z_first
713         DO ji = 2, jpim1
714            DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1
715#else
716         DO jk = 2, jpkm1
717            DO ji = 2, jpim1
718#endif
719               ! i- and j-slopes at vw-point
720               zvwslpi = 0.5 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
721               zvwslpj = 0.5 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
722               ! coef. for the vertical derivative
723               zcoef0 = e1v(ji,jj) * e2v(ji,jj) / fse3v(ji,jj,jk)   &
724                      * ( zvwslpi * zvwslpi + zvwslpj * zvwslpj )
725               ! weights for the i-k, j-k averaging at f- and t-points, resp.
726               zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)   &
727                             + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ), 1. )
728               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)   &
729                             + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ), 1. )
730               ! coef. for the horizontal derivatives
731               zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi
732               zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj
733               ! vertical flux on pv field
734               zfvw(ji,jk) = vmask(ji,jj,jk) *   &
735                           (  zcoef0 * ( pv  (ji,jj,jk-1) - pv  (ji,jj,jk) )   &
736                            + zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1)     &
737                                        +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   &
738                            + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     &
739                                        +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) )  )
740            END DO
741         END DO
742
743
744         ! II.3 Divergence of vertical fluxes added to the horizontal divergence
745         ! ---------------------------------------------------------------------
746
747         IF( kahm == 1 ) THEN
748            ! multiply the laplacian by the eddy viscosity coefficient
749#if defined key_z_first
750            DO ji = 2, jpim1
751               DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
752#else
753            DO jk = 1, jpkm1
754               DO ji = 2, jpim1
755#endif
756                  ! eddy coef. divided by the volume element
757                  zbur = fsahmu(ji,jj,jk) / ( e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) )
758                  zbvr = fsahmv(ji,jj,jk) / ( e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) )
759                  ! vertical divergence
760                  zuav = zfuw(ji,jk) - zfuw(ji,jk+1)
761                  zvav = zfvw(ji,jk) - zfvw(ji,jk+1)
762                  ! harmonic operator applied to (pu,pv) and multiply by ahm
763                  plu(ji,jj,jk) = ( plu(ji,jj,jk) + zuav ) * zbur
764                  plv(ji,jj,jk) = ( plv(ji,jj,jk) + zvav ) * zbvr
765               END DO
766            END DO
767         ELSEIF( kahm == 2 ) THEN
768            ! second call, no multiplication
769#if defined key_z_first
770            DO ji = 2, jpim1
771               DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
772#else
773            DO jk = 1, jpkm1
774               DO ji = 2, jpim1
775#endif
776                  ! inverse of the volume element
777                  zbur = 1. / ( e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) )
778                  zbvr = 1. / ( e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) )
779                  ! vertical divergence
780                  zuav = zfuw(ji,jk) - zfuw(ji,jk+1)
781                  zvav = zfvw(ji,jk) - zfvw(ji,jk+1)
782                  ! harmonic operator applied to (pu,pv)
783                  plu(ji,jj,jk) = ( plu(ji,jj,jk) + zuav ) * zbur
784                  plv(ji,jj,jk) = ( plv(ji,jj,jk) + zvav ) * zbvr
785               END DO
786            END DO
787         ELSE
788            IF(lwp)WRITE(numout,*) ' ldfguv: kahm= 1 or 2, here =', kahm
789            IF(lwp)WRITE(numout,*) '         We stop'
790            STOP 'ldfguv'
791         ENDIF
792         !                                             ! ===============
793      END DO                                           !   End of slab
794      !                                                ! ===============
795      CALL timing_stop('ldfguv_2nd','section')
796
797      IF( wrk_not_released(2, 1,2,3,4,5,6,7,8) )   CALL ctl_stop('dyn:ldfguv: failed to release workspace arrays')
798      !
799      CALL timing_stop('ldfguv','section')
800
801   END SUBROUTINE ldfguv
802
803#else
804   !!----------------------------------------------------------------------
805   !!   Dummy module :                         NO rotation of mixing tensor
806   !!----------------------------------------------------------------------
807CONTAINS
808   SUBROUTINE dyn_ldf_bilapg( kt )               ! Dummy routine
809      INTEGER, INTENT(in) :: kt
810      WRITE(*,*) 'dyn_ldf_bilapg: You should not have seen this print! error?', kt
811   END SUBROUTINE dyn_ldf_bilapg
812#endif
813
814   !!======================================================================
815END MODULE dynldf_bilapg
Note: See TracBrowser for help on using the repository browser.