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/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90

Last change on this file was 10774, checked in by andmirek, 5 years ago

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