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/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90 @ 4616

Last change on this file since 4616 was 4616, checked in by gm, 10 years ago

#1260 : see the associated wiki page for explanation
  • Property svn:keywords set to Id
File size: 22.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   !!            3.7  ! 2014-01  (F. Lemarie, G. Madec)  restructuration/simplification of ahm specification,
10   !!                 !                                  add velocity dependent coefficient and optional read in file
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   dyn_ldf_blpg : 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          ! lateral diffusion: eddy viscosity coef.
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   !
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_blpg       ! 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   !!----------------------------------------------------------------------
44   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
45   !! $Id$
46   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
47   !!----------------------------------------------------------------------
48CONTAINS
49
50   INTEGER FUNCTION dyn_ldf_blpg_alloc()
51      !!----------------------------------------------------------------------
52      !!               ***  ROUTINE dyn_ldf_blpg_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_blpg_alloc )
56         !
57      IF( dyn_ldf_blpg_alloc /= 0 )   CALL ctl_warn('dyn_ldf_blpg_alloc: failed to allocate arrays')
58   END FUNCTION dyn_ldf_blpg_alloc
59
60
61   SUBROUTINE dyn_ldf_blpg( kt )
62      !!----------------------------------------------------------------------
63      !!                   ***  ROUTINE dyn_ldf_blpg  ***
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      !!      'key_trddyn' defined: the trend is saved for diagnostics.
83      !!
84      !! ** Action  : - Update (ua,va) arrays with the before geopotential
85      !!                biharmonic mixing trend.
86      !!              - save the trend in (zwk3,zwk4) ('key_trddyn')
87      !!----------------------------------------------------------------------
88      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
89      !
90      INTEGER ::   ji, jj, jk                 ! dummy loop indices
91      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwk1, zwk2, zwk3, zwk4
92      !!----------------------------------------------------------------------
93      !
94      IF( nn_timing == 1 )  CALL timing_start('dyn_ldf_blpg')
95      !
96      CALL wrk_alloc( jpi, jpj, jpk, zwk1, zwk2, zwk3, zwk4 ) 
97      !
98      IF( kt == nit000 ) THEN
99         IF(lwp) WRITE(numout,*)
100         IF(lwp) WRITE(numout,*) 'dyn_ldf_blpg : horizontal biharmonic operator in s-coordinate'
101         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
102         !                                      ! allocate dyn_ldf_blpg arrays
103         IF( dyn_ldf_blpg_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_ldf_blpg: failed to allocate arrays')
104      ENDIF
105      !
106      zwk1(:,:,:) = 0.e0   ;   zwk3(:,:,:) = 0.e0
107      zwk2(:,:,:) = 0.e0   ;   zwk4(:,:,:) = 0.e0
108
109      ! Laplacian of (ub,vb) multiplied by ahm
110      ! -------------------------------------- 
111      CALL ldfguv( ub, vb, zwk1, zwk2, 1 )      ! rotated harmonic operator applied to (ub,vb)
112      !                                         ! and multiply by ahmu, ahmv (output in (zwk1,zwk2) )
113      CALL lbc_lnk( zwk1, 'U', -1. )   ;   CALL lbc_lnk( zwk2, 'V', -1. )     ! Lateral boundary conditions
114
115      ! Bilaplacian of (ub,vb)
116      ! ----------------------
117      CALL ldfguv( zwk1, zwk2, zwk3, zwk4, 2 )  ! rotated harmonic operator applied to (zwk1,zwk2)
118      !                                         ! (output in (zwk3,zwk4) )
119
120      ! Update the momentum trends
121      ! --------------------------
122      DO jj = 2, jpjm1               ! add the diffusive trend to the general momentum trends
123         DO jk = 1, jpkm1
124            DO ji = 2, jpim1
125               ua(ji,jj,jk) = ua(ji,jj,jk) + zwk3(ji,jj,jk)
126               va(ji,jj,jk) = va(ji,jj,jk) + zwk4(ji,jj,jk)
127            END DO
128         END DO
129      END DO
130      !
131      CALL wrk_dealloc( jpi, jpj, jpk, zwk1, zwk2, zwk3, zwk4 ) 
132      !
133      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_blpg')
134      !
135   END SUBROUTINE dyn_ldf_blpg
136
137
138   SUBROUTINE ldfguv( pu, pv, plu, plv, kahm )
139      !!----------------------------------------------------------------------
140      !!                  ***  ROUTINE ldfguv  ***
141      !!     
142      !! ** Purpose :   Apply a geopotential harmonic operator to (pu,pv)
143      !!      (defined at u- and v-points) and multiply it by the eddy
144      !!      viscosity coefficient (if kahm=1).
145      !!
146      !! ** Method  :   The harmonic operator rotated along geopotential
147      !!      surfaces is applied to (pu,pv) using the slopes of geopotential
148      !!      surfaces computed in inildf routine. The result is provided in
149      !!      (plu,plv) arrays. It is computed in 2 stepv:
150      !!
151      !!      First step: horizontal part of the operator. It is computed on
152      !!      ==========  pu as follows (idem on pv)
153      !!      horizontal fluxes :
154      !!         zftu = ahmt ( e2u*e3u/e1u di[ pu ] - e2u*uslp dk[ mi(mk(pu)) ] )
155      !!         zftv = ahmf ( e1v*e3v/e2v dj[ pu ] - e1v*vslp dk[ mj(mk(pu)) ] )
156      !!      take the horizontal divergence of the fluxes (no divided by
157      !!      the volume element :
158      !!         plu  = di-1[ zftu ] + dj-1[ zftv ]
159      !!
160      !!      Second step: vertical part of the operator. It is computed on
161      !!      ===========  pu as follows (idem on pv)
162      !!      vertical fluxes :
163      !!         zftw = e1e2t/e3w * (ahm*wslpi^2+ahm*wslpj^2)  dk-1[ pu ]
164      !!              -   e2t     *  ahm*wslpi                di[ mi(mk(pu)) ]
165      !!              -   e1t     *  ahm*wslpj                dj[ mj(mk(pu)) ]
166      !!      take the vertical divergence of the fluxes add it to the hori-
167      !!      zontal component, divide the result by the volume element :
168      !!         plu  =  zsign / (e1e2t*e3t) { plu + dk[ zftw ] }
169      !!      where  zsign=+1  if kahm =1 (laplacian or 1st pass of bilaplacian)
170      !!                  =-1  if kahm =2 (2nd pass in case of bilaplacian)
171      !!
172      !! ** Action :
173      !!        plu, plv        : partial harmonic operator applied to
174      !!                          pu and pv (all the components except
175      !!                          second order vertical derivative term)
176      !!      'key_trddyn' defined: the trend is saved for diagnostics.
177      !!----------------------------------------------------------------------
178      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu , pv    ! fields on which laplacian is applied
179      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   plu, plv   ! partial laplacian operator applied to
180      !                                                               ! pu and pv (all the components except
181      !                                                               ! second order vertical derivative term)
182      INTEGER                         , INTENT(in   ) ::   kahm       ! =1 1st call ; =2 2nd call
183      !
184      INTEGER  ::   ji, jj, jk   ! dummy loop indices
185      REAL(wp) ::   zabe1 , zabe2 , zcof1 , zcof2        ! local scalar
186      REAL(wp) ::   zcoef0, zcoef3, zcoef4, zsign        !   -      -
187      REAL(wp) ::   zbur, zbvr, zmkt, zmkf, zuav, zvav   !   -      -
188      REAL(wp) ::   zuwslpi, zuwslpj, zvwslpi, zvwslpj   !   -      -
189      REAL(wp) ::   zaht_uw, zahf_uw                     !   -      -
190      REAL(wp) ::   zaht_vw, zahf_vw                     !   -      -
191      !
192      REAL(wp), POINTER, DIMENSION(:,:) :: ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v
193      !!----------------------------------------------------------------------
194      !
195      IF( nn_timing == 1 )  CALL timing_start('ldfguv')
196      !
197      CALL wrk_alloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v ) 
198      !
199      IF    ( kahm == 1 ) THEN   ;   zsign = +1._wp
200      ELSEIF( kahm == 2 ) THEN   ;   zsign = -1._wp
201      ELSE
202         IF(lwp)WRITE(numout,*) ' ldfguv: kahm= 1 or 2, here =', kahm
203         IF(lwp)WRITE(numout,*) '         We stop'
204         STOP 'ldfguv'
205      ENDIF     
206      !
207      !                               ! ********** !   ! ===============
208      DO jk = 1, jpkm1                ! First step !   ! Horizontal slab
209         !                            ! ********** !   ! ===============
210
211         ! I.1 Vertical gradient of pu and pv at level jk and jk+1
212         ! -------------------------------------------------------
213         ! surface boundary condition: zdku(jk=1)=zdku(jk=2)
214         !                             zdkv(jk=1)=zdkv(jk=2)
215
216         zdk1u(:,:) = ( pu(:,:,jk) - pu(:,:,jk+1) ) * umask(:,:,jk+1)
217         zdk1v(:,:) = ( pv(:,:,jk) - pv(:,:,jk+1) ) * vmask(:,:,jk+1)
218
219         IF( jk == 1 ) THEN
220            zdku(:,:) = zdk1u(:,:)
221            zdkv(:,:) = zdk1v(:,:)
222         ELSE
223            zdku(:,:) = ( pu(:,:,jk-1) - pu(:,:,jk) ) * umask(:,:,jk)
224            zdkv(:,:) = ( pv(:,:,jk-1) - pv(:,:,jk) ) * vmask(:,:,jk)
225         ENDIF
226
227         !                                -----f-----
228         ! I.2 Horizontal fluxes on U          |
229         ! ------------------------===     t   u   t
230         !                                     |
231         ! i-flux at t-point              -----f-----
232         DO jj = 1, jpjm1
233            DO ji = 2, jpi
234               zabe1 = e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)
235               !
236               zmkt  = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
237                  &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
238                  !
239               zcof1 = -e2t(ji,jj) * zmkt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
240               !
241               ziut(ji,jj) = tmask(ji,jj,jk) * ahmt(ji,jj,jk) *           &
242                  &        (  zabe1 * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) )   &
243                  &         + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
244                  &                    +zdk1u(ji,jj) + zdku (ji-1,jj) )  )
245            END DO
246         END DO
247
248         ! j-flux at f-point
249         DO jj = 1, jpjm1
250            DO ji = 1, jpim1
251               zabe2 = e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)
252               !
253               zmkf  = 1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   &
254                  &           + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. )
255                  !
256               zcof2 = -e1f(ji,jj) * zmkf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
257
258!!gm caution here fmask multiplication already done in the def of ahmf...
259!!gm so in noslip.... with fmask value=2 at the coast  !!!!
260
261               !
262               zjuf(ji,jj) = fmask(ji,jj,jk) * ahmf(ji,jj,jk) *           &
263                  &        (  zabe2 * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) )   &
264                  &         + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
265                  &                    +zdk1u(ji,jj+1) + zdku (ji,jj) )  )
266            END DO
267         END DO
268
269         !                                 |   t   |
270         ! I.3 Horizontal fluxes on V      |       |
271         ! ------------------------===     f---v---f
272         !                                 |       |
273         ! i-flux at f-point               |   t   |
274         DO jj = 1, jpjm1
275            DO ji = 1, jpim1
276               zabe1 = e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
277               !
278               zmkf  = 1./MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
279                  &           + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. )
280                  !
281               zcof1 = -e2f(ji,jj) * zmkf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
282               !
283               zivf(ji,jj) = fmask(ji,jj,jk) * ahmf(ji,jj,jk) *           &
284                  &        (  zabe1 * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) )   &
285                  &         + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj)     &
286                  &                    +zdk1u(ji,jj) + zdku (ji+1,jj) )  )
287            END DO
288         END DO
289
290         ! j-flux at t-point
291         DO jj = 2, jpj
292            DO ji = 1, jpim1
293               zabe2 = e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)
294               !
295               zmkt  = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
296                  &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
297                  !
298               zcof2 = -e1t(ji,jj) * zmkt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
299               !
300               zjvt(ji,jj) = tmask(ji,jj,jk) * ahmt(ji,jj,jk) *            &
301                  &        (  zabe2 * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) )    &
302                  &         + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj)      &
303                  &                    +zdk1u(ji,jj-1) + zdku (ji,jj) )  )
304            END DO
305         END DO
306
307
308         ! I.4 Second derivative (divergence) (not divided by the volume)
309         ! ---------------------
310
311         DO jj = 2, jpjm1
312            DO ji = 2, jpim1
313               plu(ji,jj,jk) =   ziut (ji+1,jj) - ziut (ji,jj  )   &
314                  &            + zjuf (ji  ,jj) - zjuf (ji,jj-1)   
315               plv(ji,jj,jk) =   zivf (ji,jj  ) - zivf (ji-1,jj)   &
316                  &            + zjvt (ji,jj+1) - zjvt (ji,jj  )   
317            END DO
318         END DO
319         !                                             ! ===============
320      END DO                                           !   End of slab
321      !                                                ! ===============
322
323      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
324
325      !                             ! ************ !   ! ===============
326      DO jj = 2, jpjm1              !  Second step !   ! Horizontal slab
327         !                          ! ************ !   ! ===============
328
329         ! II.1 horizontal (pu,pv) gradients
330         ! ---------------------------------
331
332         DO jk = 1, jpk
333            DO ji = 2, jpi
334!!gm caution here fmask multiplication already done in the def of ahmf...
335!!gm so in noslip.... with fmask value=2 at the coast  !!!!
336               !                                ! i-gradient of u at jj
337               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( pu(ji,jj  ,jk) - pu(ji-1,jj  ,jk) )
338               !                                ! j-gradient of u and v at jj
339               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( pu(ji,jj+1,jk) - pu(ji  ,jj  ,jk) )
340               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( pv(ji,jj  ,jk) - pv(ji  ,jj-1,jk) )
341               !                                ! j-gradient of u and v at jj+1
342               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( pu(ji,jj  ,jk) - pu(ji  ,jj-1,jk) )
343               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pv(ji,jj+1,jk) - pv(ji  ,jj  ,jk) )
344            END DO
345         END DO
346         DO jk = 1, jpk
347            DO ji = 1, jpim1               ! i-gradient of v at jj
348               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( pv(ji+1,jj,jk) - pv(ji  ,jj  ,jk) )
349            END DO
350         END DO
351
352
353         ! II.2 Vertical fluxes
354         ! --------------------
355
356         ! Surface and bottom vertical fluxes set to zero
357
358         zfuw(:, 1 ) = 0._wp
359         zfvw(:, 1 ) = 0._wp
360         zfuw(:,jpk) = 0._wp
361         zfvw(:,jpk) = 0._wp
362
363         ! interior (2=<jk=<jpk-1) on pu field
364
365         DO jk = 2, jpkm1
366            DO ji = 2, jpim1
367               ! i- and j-slopes at uw-point
368               zuwslpi = 0.5 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
369               zuwslpj = 0.5 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
370               ! coef. for the vertical dirative
371               zcoef0 = e1e2u(ji,jj) / fse3u(ji,jj,jk)   &
372                  &   * ( zuwslpi * zuwslpi + zuwslpj * zuwslpj )
373               ! weights for the i-k, j-k averaging at t- and f-points, resp.
374               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)   &
375                  &          + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ), 1. )
376!!gm caution here fmask multiplication already done in the def of ahmf...
377!!gm so in noslip.... with fmask value=2 at the coast  !!!!
378               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1)   &
379                  &          + fmask(ji,jj-1,jk  )+fmask(ji,jj,jk  ), 1. )
380               zaht_uw = (   ahmt(ji,jj,jk-1) + ahmt(ji+1,jj,jk-1)    &
381                  &        + ahmt(ji,jj,jk  ) + ahmt(ji+1,jj,jk  )  ) * zmkt
382               zahf_uw = (   ahmf(ji,jj-1,jk-1) + ahmf(ji,jj,jk-1)    &
383                  &        + ahmf(ji,jj-1,jk  ) + ahmf(ji,jj,jk  )  ) * zmkf
384               ! coef. for the horitontal derivative
385               zcoef3 = - e2u(ji,jj) * zaht_uw * zuwslpi
386               zcoef4 = - e1u(ji,jj) * zahf_uw * zuwslpj
387               ! vertical flux on u field
388               zfuw(ji,jk) = umask(ji,jj,jk) *   &
389                  &        (  zcoef0 * ( pu  (ji,jj,jk-1) - pu  (ji,jj,jk) )   &
390                  &         + zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1)     &
391                  &                     +zdiu (ji,jk  ) + zdiu (ji+1,jk  ) )   &
392                  &         + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji  ,jk-1)     &
393                  &                     +zdj1u(ji,jk  ) + zdju (ji  ,jk  ) ) )
394            END DO
395         END DO
396
397         ! interior (2=<jk=<jpk-1) on pv field
398
399         DO jk = 2, jpkm1
400            DO ji = 2, jpim1
401               ! i- and j-slopes at vw-point
402               zvwslpi = 0.5 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
403               zvwslpj = 0.5 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
404               ! coef. for the vertical derivative
405               zcoef0 = e1e2v(ji,jj) / fse3v(ji,jj,jk)   &
406                  &   * ( zvwslpi * zvwslpi + zvwslpj * zvwslpj )
407!!gm caution here fmask multiplication already done in the def of ahmf...
408!!gm so in noslip.... with fmask value=2 at the coast  !!!!
409               ! weights for the i-k, j-k averaging at f- and t-points, resp.
410               zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)   &
411                  &          + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ), 1. )
412               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)   &
413                  &          + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ), 1. )
414               zahf_vw = (   ahmf(ji-1,jj,jk-1) + ahmf(ji,jj,jk-1)    &
415                  &        + ahmf(ji-1,jj,jk  ) + ahmf(ji,jj,jk  )  ) * zmkf
416               zaht_vw = (   ahmt(ji,jj,jk-1) + ahmt(ji,jj+1,jk-1)    &
417                  &        + ahmt(ji,jj,jk  ) + ahmt(ji,jj+1,jk  )  ) * zmkt
418               ! coef. for the horizontal derivatives
419               zcoef3 = - e2v(ji,jj) * zahf_vw * zvwslpi
420               zcoef4 = - e1v(ji,jj) * zaht_vw * zvwslpj
421               ! vertical flux on pv field
422               zfvw(ji,jk) = vmask(ji,jj,jk) *   &
423                  &        (  zcoef0 * ( pv  (ji,jj,jk-1) - pv  (ji,jj,jk) )   &
424                  &         + zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1)     &
425                  &                     +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   &
426                  &         + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     &
427                  &                     +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) )  )
428            END DO
429         END DO
430
431
432         ! II.3 Divergence of vertical fluxes added to the horizontal divergence
433         ! ---------------------------------------------------------------------
434         
435         DO jk = 1, jpkm1
436            DO ji = 2, jpim1
437               ! vertical divergence
438               zuav = zfuw(ji,jk) - zfuw(ji,jk+1)
439               zvav = zfvw(ji,jk) - zfvw(ji,jk+1)
440               ! harmonic operator applied to (pu,pv) and multiply by ahm
441               plu(ji,jj,jk) = zsign * ( plu(ji,jj,jk) + zuav ) / ( e1e2u(ji,jj)*fse3u(ji,jj,jk) )
442               plv(ji,jj,jk) = zsign * ( plv(ji,jj,jk) + zvav ) / ( e1e2v(ji,jj)*fse3v(ji,jj,jk) )
443            END DO
444         END DO
445         !                                             ! ===============
446      END DO                                           !   End of slab
447      !                                                ! ===============
448
449      CALL wrk_dealloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v ) 
450      !
451      IF( nn_timing == 1 )  CALL timing_stop('ldfguv')
452      !
453   END SUBROUTINE ldfguv
454
455   !!======================================================================
456END MODULE dynldf_bilapg
Note: See TracBrowser for help on using the repository browser.