source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90 @ 4990

Last change on this file since 4990 was 4990, checked in by timgraham, 6 years ago

Merged branches/2014/dev_MERGE_2014 back onto the trunk as follows:

In the working copy of branch ran:
svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk@HEAD
1 conflict in LIM_SRC_3/limdiahsb.F90
Resolved by keeping the version from dev_MERGE_2014 branch
and commited at r4989

In working copy run:
svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
to switch working copy

Run:
svn merge —reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2014/dev_MERGE_2014
to merge the branch into the trunk - no conflicts at this stage.

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