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

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

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