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

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90 @ 4488

Last change on this file since 4488 was 4488, checked in by rfurner, 10 years ago

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