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

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

Changes to allow jpk to be modified to deepest level within a subdomain. jpkorig holds original value.
  • Property svn:keywords set to Id
File size: 51.1 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, 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!      DO jk = 1, jpkm1                ! First step !   ! Horizontal slab
242         !                            ! ********** !   ! ===============
243
244         ! I.1 Vertical gradient of pu and pv at level jk and jk+1
245         ! -------------------------------------------------------
246         ! surface boundary condition: zdku(jk=1)=zdku(jk=2)
247         !                             zdkv(jk=1)=zdkv(jk=2)
248#if 0
249         DO jjf = 1, jpj
250            DO jif = 1, jpi
251
252!!$               jj= jjf
253!!$               ji = jif
254!!$
255!!$               zdk1u = ( pu(ji,jj,jk) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1)
256!!$               zdk1v = ( pv(ji,jj,jk) - pv(ji,jj,jk+1) ) * vmask(ji,jj,jk+1)
257!!$
258!!$               IF( jk == 1 ) THEN
259!!$                  zdku = zdk1u
260!!$                  zdkv = zdk1v
261!!$               ELSE
262!!$                  zdku = ( pu(ji,jj,jk-1) - pu(ji,jj,jk) ) * umask(ji,jj,jk)
263!!$                  zdkv = ( pv(ji,jj,jk-1) - pv(ji,jj,jk) ) * vmask(ji,jj,jk)
264!!$               ENDIF
265!            END DO
266!         END DO
267         !                                -----f-----
268         ! I.2 Horizontal fluxes on U          |
269         ! ------------------------===     t   u   t
270         !                                     |
271         ! i-flux at t-point              -----f-----
272!         DO jj = 1, jpjm1
273!            DO ji = 2, jpi
274!         DO jjf = 1, jpj
275!            DO jif = 1, jpi
276
277               jj = MIN(jjf, jpjm1)
278               ji = MAX(2, jif)
279
280!               zabe1 = e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)
281
282!               zmkt  = 1._wp/MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
283!                                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1._wp )
284
285!               zcof1 = -e2t(ji,jj) / MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
286!                                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1._wp )   &
287!                     * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
288
289               IF( jk == 1 )THEN
290                  ziut(ji,jj) = tmask(ji,jj,jk) *   &
291                           (  (e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) )   &
292!                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
293                            + (-e2t(ji,jj) / MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
294                                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1._wp )   &
295                     * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) ) * &
296( ( pu(ji  ,jj,jk) - pu(ji  ,jj,jk+1) ) * umask(ji,jj,jk+1) + &
297  ( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) + &
298!                                       +zdk1u(ji,jj) + zdku (ji-1,jj) )  )
299  ( pu(ji  ,jj,jk) - pu(ji  ,jj,jk+1) ) * umask(ji,jj,jk+1) + &
300  ( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) &
301)  )
302               ELSE
303                  ziut(ji,jj) = tmask(ji,jj,jk) *   &
304                           (  (e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) )   &
305!                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
306                            + (-e2t(ji,jj) / MAX(umask(ji-1,jj,jk)+umask(ji,jj,jk+1)   &
307                                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk), 1._wp )   &
308                     * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) ) * &
309( ( pu(ji  ,jj,jk-1) - pu(ji  ,jj,jk  ) ) * umask(ji,jj,jk) + &
310  ( pu(ji-1,jj,jk  ) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) +     &
311!                                       +zdk1u(ji,jj) + zdku (ji-1,jj) )  )
312  ( pu(ji  ,jj,jk  ) - pu(ji  ,jj,jk+1) ) * umask(ji,jj,jk+1) + &
313  ( pu(ji-1,jj,jk-1) - pu(ji-1,jj,jk  ) ) * umask(ji-1,jj,jk) &
314)  )
315               END IF
316!            END DO
317!         END DO
318
319         ! j-flux at f-point
320!         DO jj = 1, jpjm1
321!            DO ji = 1, jpim1
322!         DO jjf = 1, jpj
323!            DO jif = 1, jpi
324
325               jj = MIN(jjf, jpjm1)
326               ji = MIN(jif, jpim1)
327
328               !zabe2 = e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)
329
330               !zmkf  = 1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   &
331               !               + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. )
332
333               !zcof2 = -e1f(ji,jj) /MAX(umask(ji,jj+1,jk)+umask(ji,jj,jk+1)   &
334               !                       + umask(ji,jj+1,jk+1)+umask(ji,jj,jk), 1. )   &
335               !      * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
336
337               IF(jk == 1)THEN
338                  ! zdku = zdk1u
339                  zjuf(ji,jj) = fmask(ji,jj,jk) *   &
340                           (  (e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) ) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) )   &
341!                            + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
342                             + (-e1f(ji,jj) /MAX(umask(ji,jj+1,jk)+umask(ji,jj,jk+1)   &
343                                      + umask(ji,jj+1,jk+1)+umask(ji,jj,jk), 1. )   &
344                     * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) ) * &
345( &
346             ( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1)  &
347           + ( pu(ji,jj  ,jk) - pu(ji,jj  ,jk+1) ) * umask(ji,jj  ,jk+1)  &
348!                                       +zdk1u(ji,jj+1) + zdku (ji,jj) )  )
349           + ( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) &
350           + ( pu(ji,jj  ,jk) - pu(ji,jj  ,jk+1) ) * umask(ji,jj,jk+1)   &
351)  )
352               
353            ELSE
354               zjuf(ji,jj) = fmask(ji,jj,jk) *   &
355                           (  (e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) ) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) )   &
356!                            + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
357                             + (-e1f(ji,jj) /MAX(umask(ji,jj+1,jk)+umask(ji,jj,jk+1)   &
358                                      + umask(ji,jj+1,jk+1)+umask(ji,jj,jk), 1. )   &
359                     * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) ) * &
360(  &
361  (pu(ji,jj+1,jk-1) - pu(ji,jj+1,jk  ) ) * umask(ji,jj+1,jk) +   &
362  (pu(ji,jj  ,jk  ) - pu(ji,jj  ,jk+1) ) * umask(ji,jj,jk+1) +   &
363!                                       +zdk1u(ji,jj+1) + zdku (ji,jj) )  )
364  (pu(ji,jj+1,jk  ) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) + &
365  (pu(ji,jj  ,jk-1) - pu(ji,jj  ,jk  ) ) * umask(ji,jj,jk) &
366 )  )
367            ENDIF
368!            END DO
369!         END DO
370
371         !                                 |   t   |
372         ! I.3 Horizontal fluxes on V      |       |
373         ! ------------------------===     f---v---f
374         !                                 |       |
375         ! i-flux at f-point               |   t   |
376!         DO jj = 1, jpjm1
377!            DO ji = 1, jpim1
378!         DO jjf = 1, jpj
379!            DO jif = 1, jpi
380
381               jj = MIN(jjf, jpjm1)
382               ji = MIN(jif, jpim1)
383
384!               zabe1 = e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
385
386!               zmkf  = 1./MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
387!                              + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. )
388
389!               zcof1 = (-e2f(ji,jj) / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
390!                              + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. ))   &
391!                     * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
392
393               IF (jk == 1)THEN
394                  ! zdku == zdk1u
395                  zivf(ji,jj) = fmask(ji,jj,jk) *   &
396                           (  (e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) )   &
397!                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj)     &
398                             + ((-e2f(ji,jj) / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
399                              + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. ))   &
400                     * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * ( &
401  ( pu(ji  ,jj,jk) - pu(ji  ,jj,jk+1) ) * umask(ji  ,jj,jk+1) + &
402  ( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) +    &
403!                                       +zdk1u(ji,jj) + zdku (ji+1,jj) )  )
404  ( pu(ji  ,jj,jk) - pu(ji  ,jj,jk+1) ) * umask(ji  ,jj,jk+1) +  &
405  ( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) ) )
406               ELSE
407                  zivf(ji,jj) = fmask(ji,jj,jk) *   &
408                           (  (e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) )   &
409!                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj)     &
410+ ((-e2f(ji,jj) / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
411                              + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. ))   &
412                     * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * ( &
413  ( pu(ji  ,jj,jk-1) - pu(ji  ,jj,jk  ) ) * umask(ji  ,jj,jk  ) + &
414  ( pu(ji+1,jj,jk  ) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) +    &
415!                                       +zdk1u(ji,jj) + zdku (ji+1,jj) )  )
416  ( pu(ji  ,jj,jk  ) - pu(ji  ,jj,jk+1) ) * umask(ji  ,jj,jk+1) + &
417  ( pu(ji+1,jj,jk-1) - pu(ji+1,jj,jk  ) ) * umask(ji+1,jj,jk  ) )  )
418               END IF
419!            END DO
420!         END DO
421
422         ! j-flux at t-point
423!         DO jj = 2, jpj
424!            DO ji = 1, jpim1
425!         DO jjf = 1, jpj
426!            DO jif = 1, jpi
427
428               jj = MAX(2,jjf)
429               ji = MIN(jif, jpim1)
430
431               !zabe2 = e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)
432
433               !zmkt  = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
434               !               + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
435
436               !zcof2 = (-e1t(ji,jj)/MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
437               !               + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. ) )   &
438               !      * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
439
440               IF( jk == 1 )THEN
441                  zjvt(ji,jj) = tmask(ji,jj,jk) *   &
442                           (  (e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) )   &
443!                            + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj)     &
444                            + ((-e1t(ji,jj)/MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
445                              + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. ) )   &
446                     * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * ( &
447  ( pu(ji,jj-1,jk) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + &
448  ( pu(ji,jj  ,jk) - pu(ji,jj  ,jk+1) ) * umask(ji,jj,jk+1)   +  &
449!                                       +zdk1u(ji,jj-1) + zdku (ji,jj) )  )
450  ( pu(ji,jj-1,jk) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + &
451  ( pu(ji,jj  ,jk) - pu(ji,jj  ,jk+1) ) * umask(ji,jj,jk+1)  )  )
452               ELSE
453                  zjvt(ji,jj) = tmask(ji,jj,jk) *   &
454                           (  (e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) )   &
455!                            + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj)     &
456                            + ((-e1t(ji,jj)/MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
457                              + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. ) )   &
458                     * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * &
459( &
460  ( pu(ji,jj-1,jk-1) - pu(ji,jj-1,jk  ) ) * umask(ji,jj-1,jk) + &
461  ( pu(ji,jj  ,jk  ) - pu(ji,jj  ,jk+1) ) * umask(ji,jj,jk+1)  +    &
462!                                       +zdk1u(ji,jj-1) + zdku (ji,jj) )  )
463  ( pu(ji,jj-1,jk  ) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + &
464  ( pu(ji,jj  ,jk-1) - pu(ji,jj  ,jk  ) ) * umask(ji,jj,jk) )  )
465               END IF
466            END DO
467         END DO
468#endif
469
470         ! I.4 Second derivative (divergence) (not divided by the volume)
471         ! ---------------------
472
473
474         DO jj = 2, jpjm1
475            DO ji = 2, jpim1
476
477               ! Treat jk = 1 separately as is special case
478               jk = 1
479
480               plu(ji,jj,jk) = &
481! ------------- ziut (ji+1, jj) -
482                           tmask(ji+1,jj,jk) *   &
483                           (  (e2t(ji+1,jj) * fse3t(ji+1,jj,jk) / e1t(ji+1,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) )   &
484!                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
485                            + (-e2t(ji+1,jj) / MAX(  umask(ji,jj,jk  )+umask(ji+1,jj,jk+1)   &
486                                 + umask(ji,jj,jk+1)+umask(ji+1,jj,jk  ), 1._wp )   &
487                     * 0.5  * ( uslp(ji,jj,jk) + uslp(ji+1,jj,jk) ) ) * &
488( ( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) + &
489  ( pu(ji  ,jj,jk) - pu(ji  ,jj,jk+1) ) * umask(ji  ,jj,jk+1) + &
490!                                       +zdk1u(ji,jj) + zdku (ji-1,jj) )  )
491  ( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) + &
492  ( pu(ji  ,jj,jk) - pu(ji  ,jj,jk+1) ) * umask(ji  ,jj,jk+1) &
493)  ) - &
494! ------------- ziut (ji,jj  ) +
495                           tmask(ji,jj,jk) *   &
496                           (  (e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) )   &
497!                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
498                            + (-e2t(ji,jj) / MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
499                                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1._wp )   &
500                     * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) ) * &
501( ( pu(ji  ,jj,jk) - pu(ji  ,jj,jk+1) ) * umask(ji,jj,jk+1) + &
502  ( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) + &
503!                                       +zdk1u(ji,jj) + zdku (ji-1,jj) )  )
504  ( pu(ji  ,jj,jk) - pu(ji  ,jj,jk+1) ) * umask(ji,jj,jk+1) + &
505  ( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) &
506)  ) + &
507! ------------- zjuf (ji  ,jj) -
508                           fmask(ji,jj,jk) *   &
509                           (  (e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) ) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) )   &
510!                            + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
511                             + (-e1f(ji,jj) /MAX(umask(ji,jj+1,jk)+umask(ji,jj,jk+1)   &
512                                      + umask(ji,jj+1,jk+1)+umask(ji,jj,jk), 1. )   &
513                     * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) ) * &
514( &
515             ( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1)  &
516           + ( pu(ji,jj  ,jk) - pu(ji,jj  ,jk+1) ) * umask(ji,jj  ,jk+1)  &
517!                                       +zdk1u(ji,jj+1) + zdku (ji,jj) )  )
518           + ( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) &
519           + ( pu(ji,jj  ,jk) - pu(ji,jj  ,jk+1) ) * umask(ji,jj,jk+1)   &
520)  ) - &
521! ------------- zjuf (ji,jj-1)
522                           fmask(ji,jj-1,jk) *   &
523                           (  (e1f(ji,jj-1) * fse3f(ji,jj-1,jk) / e2f(ji,jj-1) ) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) )   &
524!                            + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
525                             + (-e1f(ji,jj-1) /MAX(umask(ji,jj,jk)+umask(ji,jj-1,jk+1)   &
526                                      + umask(ji,jj,jk+1)+umask(ji,jj-1,jk), 1. )   &
527                     * 0.5  * ( vslp(ji+1,jj-1,jk) + vslp(ji,jj-1,jk) ) ) * &
528( &
529             ( pu(ji,jj,jk) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1)  &
530           + ( pu(ji,jj-1  ,jk) - pu(ji,jj-1  ,jk+1) ) * umask(ji,jj-1  ,jk+1)  &
531!                                       +zdk1u(ji,jj+1) + zdku (ji,jj) )  )
532           + ( pu(ji,jj,jk) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1) &
533           + ( pu(ji,jj-1  ,jk) - pu(ji,jj-1  ,jk+1) ) * umask(ji,jj-1,jk+1)   &
534)  )
535
536
537                  plv(ji,jj,jk) = &
538! ------------- zivf (ji,jj  ) -
539fmask(ji,jj,jk) *   &
540                           (  (e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) )   &
541!                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj)     &
542                             + ((-e2f(ji,jj) / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
543                              + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. ))   &
544                     * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * ( &
545  ( pu(ji  ,jj,jk) - pu(ji  ,jj,jk+1) ) * umask(ji  ,jj,jk+1) + &
546  ( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) +    &
547!                                       +zdk1u(ji,jj) + zdku (ji+1,jj) )  )
548  ( pu(ji  ,jj,jk) - pu(ji  ,jj,jk+1) ) * umask(ji  ,jj,jk+1) +  &
549  ( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1)    &
550) ) - &
551! ------------- zivf (ji-1,jj) +
552fmask(ji-1,jj,jk) *   &
553                           (  (e2f(ji-1,jj) * fse3f(ji-1,jj,jk) / e1f(ji-1,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) )   &
554!                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj)     &
555                             + ((-e2f(ji-1,jj) / MAX(  vmask(ji,jj,jk  )+vmask(ji-1,jj,jk+1)   &
556                              + vmask(ji,jj,jk+1)+vmask(ji-1,jj,jk  ), 1. ))   &
557                     * 0.5 * ( uslp(ji-1,jj+1,jk) + uslp(ji-1,jj,jk) )) * ( &
558  ( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) + &
559  ( pu(ji  ,jj,jk) - pu(ji  ,jj,jk+1) ) * umask(ji  ,jj,jk+1) +    &
560!                                       +zdk1u(ji,jj) + zdku (ji+1,jj) )  )
561  ( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) +  &
562  ( pu(ji  ,jj,jk) - pu(ji  ,jj,jk+1) ) * umask(ji  ,jj,jk+1)    &
563) ) + &
564! ------------- zjvt (ji,jj+1) -
565tmask(ji,jj+1,jk) *   &
566                           (  (e1t(ji,jj+1) * fse3t(ji,jj+1,jk) / e2t(ji,jj+1)) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) )   &
567!                            + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj)     &
568                            + ((-e1t(ji,jj+1)/MAX(  vmask(ji,jj,jk  )+vmask(ji,jj+1,jk+1)   &
569                              + vmask(ji,jj,jk+1)+vmask(ji,jj+1,jk  ), 1. ) )   &
570                     * 0.5 * ( vslp(ji,jj,jk) + vslp(ji,jj+1,jk) )) * ( &
571  ( pu(ji,jj  ,jk) - pu(ji,jj  ,jk+1) ) * umask(ji,jj  ,jk+1) + &
572  ( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1)   +  &
573!                                       +zdk1u(ji,jj-1) + zdku (ji,jj) )  )
574  ( pu(ji,jj  ,jk) - pu(ji,jj  ,jk+1) ) * umask(ji,jj  ,jk+1) + &
575  ( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1)  & 
576)  ) - &
577! ------------- zjvt (ji,jj  )
578tmask(ji,jj,jk) *   &
579                           (  (e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) )   &
580!                            + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj)     &
581                            + ((-e1t(ji,jj)/MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
582                              + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. ) )   &
583                     * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * ( &
584  ( pu(ji,jj-1,jk) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + &
585  ( pu(ji,jj  ,jk) - pu(ji,jj  ,jk+1) ) * umask(ji,jj,jk+1)   +  &
586!                                       +zdk1u(ji,jj-1) + zdku (ji,jj) )  )
587  ( pu(ji,jj-1,jk) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + &
588  ( pu(ji,jj  ,jk) - pu(ji,jj  ,jk+1) ) * umask(ji,jj,jk+1)  &
589)  )
590
591
592               DO jk = 2, jpkm1
593
594!               plu(ji,jj,jk) = ziut (ji+1,jj) - &
595!                               ziut (ji,jj  ) + &
596!                               zjuf (ji  ,jj) - &
597!                               zjuf (ji,jj-1)
598                  plu(ji,jj,jk) = &
599! ------------- ziut (ji+1, jj  ) -
600                  tmask(ji+1,jj,jk) *   &
601                           (  (e2t(ji+1,jj) * fse3t(ji+1,jj,jk) / e1t(ji+1,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) )   &
602!                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
603                            + (-e2t(ji+1,jj) / MAX(umask(ji,jj,jk)+umask(ji+1,jj,jk+1)   &
604                                 + umask(ji,jj,jk+1)+umask(ji+1,jj,jk), 1._wp )   &
605                     * 0.5  * ( uslp(ji,jj,jk) + uslp(ji+1,jj,jk) ) ) * &
606( ( pu(ji+1,jj,jk-1) - pu(ji+1,jj,jk  ) ) * umask(ji+1,jj,jk) + &
607  ( pu(ji  ,jj,jk  ) - pu(ji  ,jj,jk+1) ) * umask(ji  ,jj,jk+1) +     &
608!                                       +zdk1u(ji,jj) + zdku (ji-1,jj) )  )
609  ( pu(ji+1,jj,jk  ) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) + &
610  ( pu(ji  ,jj,jk-1) - pu(ji  ,jj,jk  ) ) * umask(ji  ,jj,jk) &
611)  ) - &
612! ------------- ziut (ji  , jj  ) +
613                  tmask(ji,jj,jk) *   &
614                           (  (e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) )   &
615!                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
616                            + (-e2t(ji,jj) / MAX(umask(ji-1,jj,jk)+umask(ji,jj,jk+1)   &
617                                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk), 1._wp )   &
618                     * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) ) * &
619( ( pu(ji  ,jj,jk-1) - pu(ji  ,jj,jk  ) ) * umask(ji,jj,jk) + &
620  ( pu(ji-1,jj,jk  ) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) +     &
621!                                       +zdk1u(ji,jj) + zdku (ji-1,jj) )  )
622  ( pu(ji  ,jj,jk  ) - pu(ji  ,jj,jk+1) ) * umask(ji,jj,jk+1) + &
623  ( pu(ji-1,jj,jk-1) - pu(ji-1,jj,jk  ) ) * umask(ji-1,jj,jk) &
624)  ) + &
625! ------------- zjuf (ji  , jj  ) -
626fmask(ji,jj,jk) *   &
627                           (  (e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) ) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) )   &
628!                            + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
629                             + (-e1f(ji,jj) /MAX(umask(ji,jj+1,jk)+umask(ji,jj,jk+1)   &
630                                      + umask(ji,jj+1,jk+1)+umask(ji,jj,jk), 1. )   &
631                     * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) ) * &
632(  &
633  (pu(ji,jj+1,jk-1) - pu(ji,jj+1,jk  ) ) * umask(ji,jj+1,jk) +   &
634  (pu(ji,jj  ,jk  ) - pu(ji,jj  ,jk+1) ) * umask(ji,jj,jk+1) +   &
635!                                       +zdk1u(ji,jj+1) + zdku (ji,jj) )  )
636  (pu(ji,jj+1,jk  ) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) + &
637  (pu(ji,jj  ,jk-1) - pu(ji,jj  ,jk  ) ) * umask(ji,jj,jk) &
638 )  ) - &
639! ------------- zjuf (ji  , jj-1)
640fmask(ji,jj-1,jk) *   &
641                           (  (e1f(ji,jj-1) * fse3f(ji,jj-1,jk) / e2f(ji,jj-1) ) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) )   &
642!                            + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
643                             + (-e1f(ji,jj-1) /MAX(umask(ji,jj,jk)+umask(ji,jj-1,jk+1)   &
644                                      + umask(ji,jj,jk+1)+umask(ji,jj-1,jk), 1. )   &
645                     * 0.5  * ( vslp(ji+1,jj-1,jk) + vslp(ji,jj-1,jk) ) ) * &
646(  &
647  (pu(ji,jj,jk-1) - pu(ji,jj,jk  ) ) * umask(ji,jj,jk) +   &
648  (pu(ji,jj-1,jk  ) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) +   &
649!                                       +zdk1u(ji,jj+1) + zdku (ji,jj) )  )
650  (pu(ji,jj,jk  ) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1) + &
651  (pu(ji,jj-1,jk-1) - pu(ji,jj-1  ,jk  ) ) * umask(ji,jj-1,jk) &
652 )  )
653
654
655!               plv(ji,jj,jk) = zivf (ji,jj  ) - &
656!                               zivf (ji-1,jj) + &
657!                               zjvt (ji,jj+1) - &
658!                               zjvt (ji,jj  )
659                  plv(ji,jj,jk) = &
660! ------------- zivf (ji,jj  ) -
661fmask(ji,jj,jk) *   &
662                           (  (e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) )   &
663!                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj)     &
664+ ((-e2f(ji,jj) / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
665                              + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. ))   &
666                     * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * ( &
667  ( pu(ji  ,jj,jk-1) - pu(ji  ,jj,jk  ) ) * umask(ji  ,jj,jk  ) + &
668  ( pu(ji+1,jj,jk  ) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) +    &
669!                                       +zdk1u(ji,jj) + zdku (ji+1,jj) )  )
670  ( pu(ji  ,jj,jk  ) - pu(ji  ,jj,jk+1) ) * umask(ji  ,jj,jk+1) + &
671  ( pu(ji+1,jj,jk-1) - pu(ji+1,jj,jk  ) ) * umask(ji+1,jj,jk  )   &
672)  ) - &
673! ------------- zivf (ji-1,jj) +
674fmask(ji-1,jj,jk) *   &
675                           (  (e2f(ji-1,jj) * fse3f(ji-1,jj,jk) / e1f(ji-1,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) )   &
676!                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj)     &
677+ ((-e2f(ji-1,jj) / MAX(  vmask(ji,jj,jk  )+vmask(ji-1,jj,jk+1)   &
678                              + vmask(ji,jj,jk+1)+vmask(ji-1,jj,jk  ), 1. ))   &
679                     * 0.5 * ( uslp(ji-1,jj+1,jk) + uslp(ji-1,jj,jk) )) * ( &
680  ( pu(ji-1 ,jj,jk-1) - pu(ji-1  ,jj,jk  ) ) * umask(ji-1  ,jj,jk  ) + &
681  ( pu(ji,jj,jk  ) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1) +    &
682!                                       +zdk1u(ji,jj) + zdku (ji+1,jj) )  )
683  ( pu(ji-1  ,jj,jk  ) - pu(ji-1  ,jj,jk+1) ) * umask(ji-1  ,jj,jk+1) + &
684  ( pu(ji,jj,jk-1) - pu(ji,jj,jk  ) ) * umask(ji,jj,jk  )  &
685)  ) + &
686! ------------- zjvt (ji,jj+1) -
687tmask(ji,jj+1,jk) *   &
688                           (  (e1t(ji,jj+1) * fse3t(ji,jj+1,jk) / e2t(ji,jj+1)) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) )   &
689!                            + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj)     &
690                            + ((-e1t(ji,jj+1)/MAX(  vmask(ji,jj,jk  )+vmask(ji,jj+1,jk+1)   &
691                              + vmask(ji,jj,jk+1)+vmask(ji,jj+1,jk  ), 1. ) )   &
692                     * 0.5 * ( vslp(ji,jj,jk) + vslp(ji,jj+1,jk) )) * &
693( &
694  ( pu(ji,jj,jk-1) - pu(ji,jj,jk  ) ) * umask(ji,jj,jk) + &
695  ( pu(ji,jj+1  ,jk  ) - pu(ji,jj+1  ,jk+1) ) * umask(ji,jj+1,jk+1)  +    &
696!                                       +zdk1u(ji,jj-1) + zdku (ji,jj) )  )
697  ( pu(ji,jj,jk  ) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1) + &
698  ( pu(ji,jj+1  ,jk-1) - pu(ji,jj+1  ,jk  ) ) * umask(ji,jj+1,jk) &
699)  ) - &
700! ------------- zjvt (ji,jj  )
701tmask(ji,jj,jk) *   &
702                           (  (e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) )   &
703!                            + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj)     &
704                            + ((-e1t(ji,jj)/MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
705                              + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. ) )   &
706                     * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * &
707( &
708  ( pu(ji,jj-1,jk-1) - pu(ji,jj-1,jk  ) ) * umask(ji,jj-1,jk) + &
709  ( pu(ji,jj  ,jk  ) - pu(ji,jj  ,jk+1) ) * umask(ji,jj,jk+1)  +    &
710!                                       +zdk1u(ji,jj-1) + zdku (ji,jj) )  )
711  ( pu(ji,jj-1,jk  ) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + &
712  ( pu(ji,jj  ,jk-1) - pu(ji,jj  ,jk  ) ) * umask(ji,jj,jk) )  )
713
714            END DO
715         END DO
716
717         !                                             ! ===============
718      END DO                                           !   End of slab
719      !                                                ! ===============
720#else
721      !                               ! ********** !   ! ===============
722      DO jk = 1, jpkm1                ! First step !   ! Horizontal slab
723         !                            ! ********** !   ! ===============
724
725         ! I.1 Vertical gradient of pu and pv at level jk and jk+1
726         ! -------------------------------------------------------
727         ! surface boundary condition: zdku(jk=1)=zdku(jk=2)
728         !                             zdkv(jk=1)=zdkv(jk=2)
729
730         zdk1u(:,:) = ( pu(:,:,jk) - pu(:,:,jk+1) ) * umask(:,:,jk+1)
731         zdk1v(:,:) = ( pv(:,:,jk) - pv(:,:,jk+1) ) * vmask(:,:,jk+1)
732
733         IF( jk == 1 ) THEN
734            zdku(:,:) = zdk1u(:,:)
735            zdkv(:,:) = zdk1v(:,:)
736         ELSE
737            zdku(:,:) = ( pu(:,:,jk-1) - pu(:,:,jk) ) * umask(:,:,jk)
738            zdkv(:,:) = ( pv(:,:,jk-1) - pv(:,:,jk) ) * vmask(:,:,jk)
739         ENDIF
740
741         !                                -----f-----
742         ! I.2 Horizontal fluxes on U          |
743         ! ------------------------===     t   u   t
744         !                                     |
745         ! i-flux at t-point              -----f-----
746         DO jj = 1, jpjm1
747            DO ji = 2, jpi
748               zabe1 = e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)
749
750               zmkt  = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
751                              + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
752
753               zcof1 = -e2t(ji,jj) * zmkt   &
754                     * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
755
756               ziut(ji,jj) = tmask(ji,jj,jk) *   &
757                           (  zabe1 * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) )   &
758                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
759                                       +zdk1u(ji,jj) + zdku (ji-1,jj) )  )
760            END DO
761         END DO
762
763         ! j-flux at f-point
764         DO jj = 1, jpjm1
765            DO ji = 1, jpim1
766               zabe2 = e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)
767
768               zmkf  = 1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   &
769                              + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. )
770
771               zcof2 = -e1f(ji,jj) * zmkf   &
772                     * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
773
774               zjuf(ji,jj) = fmask(ji,jj,jk) *   &
775                           (  zabe2 * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) )   &
776                            + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
777                                       +zdk1u(ji,jj+1) + zdku (ji,jj) )  )
778            END DO
779         END DO
780
781         !                                 |   t   |
782         ! I.3 Horizontal fluxes on V      |       |
783         ! ------------------------===     f---v---f
784         !                                 |       |
785         ! i-flux at f-point               |   t   |
786         DO jj = 1, jpjm1
787            DO ji = 1, jpim1
788               zabe1 = e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
789
790               zmkf  = 1./MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
791                              + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. )
792
793               zcof1 = -e2f(ji,jj) * zmkf   &
794                     * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
795
796               zivf(ji,jj) = fmask(ji,jj,jk) *   &
797                           (  zabe1 * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) )   &
798                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj)     &
799                                       +zdk1u(ji,jj) + zdku (ji+1,jj) )  )
800            END DO
801         END DO
802
803         ! j-flux at t-point
804         DO jj = 2, jpj
805            DO ji = 1, jpim1
806               zabe2 = e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)
807
808               zmkt  = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
809                              + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
810
811               zcof2 = -e1t(ji,jj) * zmkt   &
812                     * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
813
814               zjvt(ji,jj) = tmask(ji,jj,jk) *   &
815                           (  zabe2 * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) )   &
816                            + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj)     &
817                                       +zdk1u(ji,jj-1) + zdku (ji,jj) )  )
818            END DO
819         END DO
820
821
822         ! I.4 Second derivative (divergence) (not divided by the volume)
823         ! ---------------------
824
825         DO jj = 2, jpjm1
826            DO ji = 2, jpim1
827               plu(ji,jj,jk) = ziut (ji+1,jj) - ziut (ji,jj  )   &
828                             + zjuf (ji  ,jj) - zjuf (ji,jj-1)
829               plv(ji,jj,jk) = zivf (ji,jj  ) - zivf (ji-1,jj)   &
830                             + zjvt (ji,jj+1) - zjvt (ji,jj  ) 
831            END DO
832         END DO
833
834         !                                             ! ===============
835      END DO                                           !   End of slab
836      !                                                ! ===============
837#endif
838      CALL timing_stop('ldfguv_1st','section')
839
840      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
841
842      CALL timing_start('ldfguv_2nd')
843      !                             ! ************ !   ! ===============
844      DO jj = 2, jpjm1              !  Second step !   ! Horizontal slab
845         !                          ! ************ !   ! ===============
846
847         ! II.1 horizontal (pu,pv) gradients
848         ! ---------------------------------
849
850#if defined key_z_first
851         DO ji = 2, jpi
852            DO jk = 1, jpk
853               ! i-gradient of u at jj
854               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( pu(ji,jj  ,jk) - pu(ji-1,jj  ,jk) )
855               ! j-gradient of u and v at jj
856               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( pu(ji,jj+1,jk) - pu(ji  ,jj  ,jk) )
857               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( pv(ji,jj  ,jk) - pv(ji  ,jj-1,jk) )
858               ! j-gradient of u and v at jj+1
859               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( pu(ji,jj  ,jk) - pu(ji  ,jj-1,jk) )
860               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pv(ji,jj+1,jk) - pv(ji  ,jj  ,jk) )
861            END DO
862         END DO
863         DO ji = 1, jpim1
864            DO jk = 1, jpk
865               ! i-gradient of v at jj
866               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( pv(ji+1,jj,jk) - pv(ji  ,jj  ,jk) )
867            END DO
868         END DO
869#else
870         DO jk = 1, jpk
871            DO ji = 2, jpi
872               ! i-gradient of u at jj
873               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( pu(ji,jj  ,jk) - pu(ji-1,jj  ,jk) )
874               ! j-gradient of u and v at jj
875               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( pu(ji,jj+1,jk) - pu(ji  ,jj  ,jk) )
876               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( pv(ji,jj  ,jk) - pv(ji  ,jj-1,jk) )
877               ! j-gradient of u and v at jj+1
878               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( pu(ji,jj  ,jk) - pu(ji  ,jj-1,jk) )
879               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pv(ji,jj+1,jk) - pv(ji  ,jj  ,jk) )
880            END DO
881         END DO
882         DO jk = 1, jpk
883            DO ji = 1, jpim1
884               ! i-gradient of v at jj
885               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( pv(ji+1,jj,jk) - pv(ji  ,jj  ,jk) )
886            END DO
887         END DO
888#endif
889
890         ! II.2 Vertical fluxes
891         ! --------------------
892
893         ! Surface and bottom vertical fluxes set to zero
894
895         zfuw(:, 1 ) = 0.e0
896         zfvw(:, 1 ) = 0.e0
897         zfuw(:,jpk) = 0.e0
898         zfvw(:,jpk) = 0.e0
899
900         ! interior (2=<jk=<jpk-1) on pu field
901
902#if defined key_z_first
903         DO ji = 2, jpim1
904            DO jk = 2, jpkm1
905#else
906         DO jk = 2, jpkm1
907            DO ji = 2, jpim1
908#endif
909               ! i- and j-slopes at uw-point
910               zuwslpi = 0.5 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
911               zuwslpj = 0.5 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
912               ! coef. for the vertical dirative
913               zcoef0 = e1u(ji,jj) * e2u(ji,jj) / fse3u(ji,jj,jk)   &
914                      * ( zuwslpi * zuwslpi + zuwslpj * zuwslpj )
915               ! weights for the i-k, j-k averaging at t- and f-points, resp.
916               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)   &
917                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ), 1. )
918               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1)   &
919                             + fmask(ji,jj-1,jk  )+fmask(ji,jj,jk  ), 1. )
920               ! coef. for the horitontal derivative
921               zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi
922               zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj
923               ! vertical flux on u field
924               zfuw(ji,jk) = umask(ji,jj,jk) *   &
925                           (  zcoef0 * ( pu  (ji,jj,jk-1) - pu  (ji,jj,jk) )   &
926                            + zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1)     &
927                                        +zdiu (ji,jk  ) + zdiu (ji+1,jk  ) )   &
928                            + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji  ,jk-1)     &
929                                        +zdj1u(ji,jk  ) + zdju (ji  ,jk  ) ) )
930            END DO
931         END DO
932
933         ! interior (2=<jk=<jpk-1) on pv field
934
935#if defined key_z_first
936         DO ji = 2, jpim1
937            DO jk = 2, jpkm1
938#else
939         DO jk = 2, jpkm1
940            DO ji = 2, jpim1
941#endif
942               ! i- and j-slopes at vw-point
943               zvwslpi = 0.5 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
944               zvwslpj = 0.5 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
945               ! coef. for the vertical derivative
946               zcoef0 = e1v(ji,jj) * e2v(ji,jj) / fse3v(ji,jj,jk)   &
947                      * ( zvwslpi * zvwslpi + zvwslpj * zvwslpj )
948               ! weights for the i-k, j-k averaging at f- and t-points, resp.
949               zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)   &
950                             + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ), 1. )
951               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)   &
952                             + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ), 1. )
953               ! coef. for the horizontal derivatives
954               zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi
955               zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj
956               ! vertical flux on pv field
957               zfvw(ji,jk) = vmask(ji,jj,jk) *   &
958                           (  zcoef0 * ( pv  (ji,jj,jk-1) - pv  (ji,jj,jk) )   &
959                            + zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1)     &
960                                        +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   &
961                            + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     &
962                                        +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) )  )
963            END DO
964         END DO
965
966
967         ! II.3 Divergence of vertical fluxes added to the horizontal divergence
968         ! ---------------------------------------------------------------------
969
970         IF( kahm == 1 ) THEN
971            ! multiply the laplacian by the eddy viscosity coefficient
972#if defined key_z_first
973            DO ji = 2, jpim1
974               DO jk = 1, jpkm1
975#else
976            DO jk = 1, jpkm1
977               DO ji = 2, jpim1
978#endif
979                  ! eddy coef. divided by the volume element
980                  zbur = fsahmu(ji,jj,jk) / ( e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) )
981                  zbvr = fsahmv(ji,jj,jk) / ( e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) )
982                  ! vertical divergence
983                  zuav = zfuw(ji,jk) - zfuw(ji,jk+1)
984                  zvav = zfvw(ji,jk) - zfvw(ji,jk+1)
985                  ! harmonic operator applied to (pu,pv) and multiply by ahm
986                  plu(ji,jj,jk) = ( plu(ji,jj,jk) + zuav ) * zbur
987                  plv(ji,jj,jk) = ( plv(ji,jj,jk) + zvav ) * zbvr
988               END DO
989            END DO
990         ELSEIF( kahm == 2 ) THEN
991            ! second call, no multiplication
992#if defined key_z_first
993            DO ji = 2, jpim1
994               DO jk = 1, jpkm1
995#else
996            DO jk = 1, jpkm1
997               DO ji = 2, jpim1
998#endif
999                  ! inverse of the volume element
1000                  zbur = 1. / ( e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) )
1001                  zbvr = 1. / ( e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) )
1002                  ! vertical divergence
1003                  zuav = zfuw(ji,jk) - zfuw(ji,jk+1)
1004                  zvav = zfvw(ji,jk) - zfvw(ji,jk+1)
1005                  ! harmonic operator applied to (pu,pv)
1006                  plu(ji,jj,jk) = ( plu(ji,jj,jk) + zuav ) * zbur
1007                  plv(ji,jj,jk) = ( plv(ji,jj,jk) + zvav ) * zbvr
1008               END DO
1009            END DO
1010         ELSE
1011            IF(lwp)WRITE(numout,*) ' ldfguv: kahm= 1 or 2, here =', kahm
1012            IF(lwp)WRITE(numout,*) '         We stop'
1013            STOP 'ldfguv'
1014         ENDIF
1015         !                                             ! ===============
1016      END DO                                           !   End of slab
1017      !                                                ! ===============
1018      CALL timing_stop('ldfguv_2nd','section')
1019
1020      IF( wrk_not_released(2, 1,2,3,4,5,6,7,8) )   CALL ctl_stop('dyn:ldfguv: failed to release workspace arrays')
1021      !
1022      CALL timing_stop('ldfguv','section')
1023
1024   END SUBROUTINE ldfguv
1025
1026#else
1027   !!----------------------------------------------------------------------
1028   !!   Dummy module :                         NO rotation of mixing tensor
1029   !!----------------------------------------------------------------------
1030CONTAINS
1031   SUBROUTINE dyn_ldf_bilapg( kt )               ! Dummy routine
1032      INTEGER, INTENT(in) :: kt
1033      WRITE(*,*) 'dyn_ldf_bilapg: You should not have seen this print! error?', kt
1034   END SUBROUTINE dyn_ldf_bilapg
1035#endif
1036
1037   !!======================================================================
1038END MODULE dynldf_bilapg
Note: See TracBrowser for help on using the repository browser.