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

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

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