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_iso.F90 in branches/UKMO/dev_r5107_restart_func_and_date/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/UKMO/dev_r5107_restart_func_and_date/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90 @ 5500

Last change on this file since 5500 was 5500, checked in by dancopsey, 9 years ago

Removed SVN keywords.

File size: 21.5 KB
Line 
1MODULE dynldf_iso
2   !!======================================================================
3   !!                     ***  MODULE  dynldf_iso  ***
4   !! Ocean dynamics:  lateral viscosity trend
5   !!======================================================================
6   !! History :  OPA  !  97-07  (G. Madec)  Original code
7   !!  NEMO      1.0  !  2002-08  (G. Madec)  F90: Free form and module
8   !!             -   !  2004-08  (C. Talandier) New trends organization
9   !!            2.0  !  2005-11  (G. Madec)  s-coordinate: horizontal diffusion
10   !!----------------------------------------------------------------------
11#if defined key_ldfslp   ||   defined key_esopa
12   !!----------------------------------------------------------------------
13   !!   'key_ldfslp'                      slopes of the direction of mixing
14   !!----------------------------------------------------------------------
15   !!   dyn_ldf_iso  : update the momentum trend with the horizontal part
16   !!                  of the lateral diffusion using isopycnal or horizon-
17   !!                  tal s-coordinate laplacian operator.
18   !!----------------------------------------------------------------------
19   USE oce             ! ocean dynamics and tracers
20   USE dom_oce         ! ocean space and time domain
21   USE ldfdyn_oce      ! ocean dynamics lateral physics
22   USE ldftra_oce      ! ocean tracer   lateral physics
23   USE zdf_oce         ! ocean vertical physics
24   USE ldfslp          ! iso-neutral slopes
25   !
26   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
27   USE in_out_manager  ! I/O manager
28   USE lib_mpp         ! MPP library
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_iso           ! called by step.F90
37   PUBLIC   dyn_ldf_iso_alloc     ! called by nemogcm.F90
38
39   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zdiu, zdju, zdj1u   ! 2D workspace (dyn_ldf_iso)
40   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfvw, zdiv, zdjv, zdj1v   !  -      -
41
42   !! * Substitutions
43#  include "domzgr_substitute.h90"
44#  include "ldfdyn_substitute.h90"
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !! NEMO/OPA 3.3 , NEMO Consortium (2011)
48   !! $Id$
49   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
50   !!----------------------------------------------------------------------
51CONTAINS
52
53   INTEGER FUNCTION dyn_ldf_iso_alloc()
54      !!----------------------------------------------------------------------
55      !!                  ***  ROUTINE dyn_ldf_iso_alloc  ***
56      !!----------------------------------------------------------------------
57      ALLOCATE( zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) ,     & 
58         &      zfvw(jpi,jpk) , zdiv(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_iso_alloc )
59         !
60      IF( dyn_ldf_iso_alloc /= 0 )   CALL ctl_warn('dyn_ldf_iso_alloc: array allocate failed.')
61   END FUNCTION dyn_ldf_iso_alloc
62
63
64   SUBROUTINE dyn_ldf_iso( kt )
65      !!----------------------------------------------------------------------
66      !!                     ***  ROUTINE dyn_ldf_iso  ***
67      !!                       
68      !! ** Purpose :   Compute the before trend of the rotated laplacian
69      !!      operator of lateral momentum diffusion except the diagonal
70      !!      vertical term that will be computed in dynzdf module. Add it
71      !!      to the general trend of momentum equation.
72      !!
73      !! ** Method :
74      !!        The momentum lateral diffusive trend is provided by a 2nd
75      !!      order operator rotated along neutral or geopotential surfaces
76      !!      (in s-coordinates).
77      !!      It is computed using before fields (forward in time) and isopyc-
78      !!      nal or geopotential slopes computed in routine ldfslp.
79      !!      Here, u and v components are considered as 2 independent scalar
80      !!      fields. Therefore, the property of splitting divergent and rota-
81      !!      tional part of the flow of the standard, z-coordinate laplacian
82      !!      momentum diffusion is lost.
83      !!      horizontal fluxes associated with the rotated lateral mixing:
84      !!      u-component:
85      !!         ziut = ( ahmt + ahmb0 ) e2t * e3t / e1t  di[ ub ]
86      !!               -      ahmt       e2t * mi-1(uslp) dk[ mi(mk(ub)) ]
87      !!         zjuf = ( ahmf + ahmb0 ) e1f * e3f / e2f  dj[ ub ]
88      !!               -      ahmf       e1f * mi(vslp)   dk[ mj(mk(ub)) ]
89      !!      v-component:
90      !!         zivf = ( ahmf + ahmb0 ) e2t * e3t / e1t  di[ vb ]
91      !!               -      ahmf       e2t * mj(uslp)   dk[ mi(mk(vb)) ]
92      !!         zjvt = ( ahmt + ahmb0 ) e1f * e3f / e2f  dj[ ub ]
93      !!               -      ahmt       e1f * mj-1(vslp) dk[ mj(mk(vb)) ]
94      !!      take the horizontal divergence of the fluxes:
95      !!         diffu = 1/(e1u*e2u*e3u) {  di  [ ziut ] + dj-1[ zjuf ]  }
96      !!         diffv = 1/(e1v*e2v*e3v) {  di-1[ zivf ] + dj  [ zjvt ]  }
97      !!      Add this trend to the general trend (ua,va):
98      !!         ua = ua + diffu
99      !!      CAUTION: here the isopycnal part is with a coeff. of aht. This
100      !!      should be modified for applications others than orca_r2 (!!bug)
101      !!
102      !! ** Action :
103      !!        Update (ua,va) arrays with the before geopotential biharmonic
104      !!      mixing trend.
105      !!        Update (avmu,avmv) to accompt for the diagonal vertical component
106      !!      of the rotated operator in dynzdf module
107      !!----------------------------------------------------------------------
108      !
109      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
110      !
111      INTEGER  ::   ji, jj, jk   ! dummy loop indices
112      REAL(wp) ::   zabe1, zabe2, zcof1, zcof2                       ! local scalars
113      REAL(wp) ::   zmskt, zmskf, zbu, zbv, zuah, zvah               !   -      -
114      REAL(wp) ::   zcoef0, zcoef3, zcoef4, zmkt, zmkf               !   -      -
115      REAL(wp) ::   zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj   !   -      -
116      !
117      REAL(wp), POINTER, DIMENSION(:,:) :: ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v
118      !!----------------------------------------------------------------------
119      !
120      IF( nn_timing == 1 )  CALL timing_start('dyn_ldf_iso')
121      !
122      CALL wrk_alloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v ) 
123      !
124      IF( kt == nit000 ) THEN
125         IF(lwp) WRITE(numout,*)
126         IF(lwp) WRITE(numout,*) 'dyn_ldf_iso : iso-neutral laplacian diffusive operator or '
127         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate horizontal diffusive operator'
128         !                                      ! allocate dyn_ldf_bilap arrays
129         IF( dyn_ldf_iso_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_ldf_iso: failed to allocate arrays')
130      ENDIF
131
132      ! s-coordinate: Iso-level diffusion on tracer, but geopotential level diffusion on momentum
133      IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
134         !
135         DO jk = 1, jpk         ! set the slopes of iso-level
136            DO jj = 2, jpjm1
137               DO ji = 2, jpim1
138                  uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk)
139                  vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk)
140                  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
141                  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
142               END DO
143            END DO
144         END DO
145         ! Lateral boundary conditions on the slopes
146         CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
147         CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
148 
149!!bug
150         IF( kt == nit000 ) then
151            IF(lwp) WRITE(numout,*) ' max slop: u', SQRT( MAXVAL(uslp*uslp)), ' v ', SQRT(MAXVAL(vslp)),  &
152               &                             ' wi', sqrt(MAXVAL(wslpi))     , ' wj', sqrt(MAXVAL(wslpj))
153         endif
154!!end
155      ENDIF
156
157      !                                                ! ===============
158      DO jk = 1, jpkm1                                 ! Horizontal slab
159         !                                             ! ===============
160
161         ! Vertical u- and v-shears at level jk and jk+1
162         ! ---------------------------------------------
163         ! surface boundary condition: zdku(jk=1)=zdku(jk=2)
164         !                             zdkv(jk=1)=zdkv(jk=2)
165
166         zdk1u(:,:) = ( ub(:,:,jk) -ub(:,:,jk+1) ) * umask(:,:,jk+1)
167         zdk1v(:,:) = ( vb(:,:,jk) -vb(:,:,jk+1) ) * vmask(:,:,jk+1)
168
169         IF( jk == 1 ) THEN
170            zdku(:,:) = zdk1u(:,:)
171            zdkv(:,:) = zdk1v(:,:)
172         ELSE
173            zdku(:,:) = ( ub(:,:,jk-1) - ub(:,:,jk) ) * umask(:,:,jk)
174            zdkv(:,:) = ( vb(:,:,jk-1) - vb(:,:,jk) ) * vmask(:,:,jk)
175         ENDIF
176
177         !                               -----f-----
178         ! Horizontal fluxes on U             | 
179         ! --------------------===        t   u   t
180         !                                    | 
181         ! i-flux at t-point             -----f-----
182
183         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
184            DO jj = 2, jpjm1
185               DO ji = fs_2, jpi   ! vector opt.
186                  zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji,jj)
187
188                  zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
189                     &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
190
191                  zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
192   
193                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   &
194                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
195                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
196               END DO
197            END DO
198         ELSE                   ! other coordinate system (zco or sco) : e3t
199            DO jj = 2, jpjm1
200               DO ji = fs_2, jpi   ! vector opt.
201                  zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)
202
203                  zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
204                     &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
205
206                  zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
207
208                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   &
209                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
210                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
211               END DO
212            END DO
213         ENDIF
214
215         ! j-flux at f-point
216         DO jj = 1, jpjm1
217            DO ji = 1, fs_jpim1   ! vector opt.
218               zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)
219
220               zmskf = 1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   &
221                  &           + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. )
222
223               zcof2 = - aht0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
224
225               zjuf(ji,jj) = (  zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) )   &
226                  &           + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
227                  &                      +zdk1u(ji,jj+1) + zdku (ji,jj) )  ) * fmask(ji,jj,jk)
228            END DO
229         END DO
230
231         !                                |   t   |
232         ! Horizontal fluxes on V         |       |
233         ! --------------------===        f---v---f
234         !                                |       |
235         ! i-flux at f-point              |   t   |
236
237         DO jj = 2, jpjm1
238            DO ji = 1, fs_jpim1   ! vector opt.
239               zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
240
241               zmskf = 1./MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
242                  &           + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. )
243
244               zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
245
246               zivf(ji,jj) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )   &
247                  &           + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj)     &
248                  &                      +zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk)
249            END DO
250         END DO
251
252         ! j-flux at t-point
253         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
254            DO jj = 2, jpj
255               DO ji = 1, fs_jpim1   ! vector opt.
256                  zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / e2t(ji,jj)
257
258                  zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
259                     &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
260
261                  zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
262
263                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
264                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     &
265                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
266               END DO
267            END DO
268         ELSE                   ! other coordinate system (zco or sco) : e3t
269            DO jj = 2, jpj
270               DO ji = 1, fs_jpim1   ! vector opt.
271                  zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)
272
273                  zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
274                     &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
275
276                  zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
277
278                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
279                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     &
280                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
281               END DO
282            END DO
283         ENDIF
284
285
286         ! Second derivative (divergence) and add to the general trend
287         ! -----------------------------------------------------------
288
289         DO jj = 2, jpjm1
290            DO ji = 2, jpim1          !! Question vectop possible??? !!bug
291               ! volume elements
292               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
293               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
294               ! horizontal component of isopycnal momentum diffusive trends
295               zuah =( ziut (ji+1,jj) - ziut (ji,jj  ) +   &
296                  &    zjuf (ji  ,jj) - zjuf (ji,jj-1)  ) / zbu
297               zvah =( zivf (ji,jj  ) - zivf (ji-1,jj) +   &
298                  &    zjvt (ji,jj+1) - zjvt (ji,jj  )  ) / zbv
299               ! add the trends to the general trends
300               ua (ji,jj,jk) = ua (ji,jj,jk) + zuah
301               va (ji,jj,jk) = va (ji,jj,jk) + zvah
302            END DO
303         END DO
304         !                                             ! ===============
305      END DO                                           !   End of slab
306      !                                                ! ===============
307
308      ! print sum trends (used for debugging)
309      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ldfh - Ua: ', mask1=umask, &
310         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
311
312
313      !                                                ! ===============
314      DO jj = 2, jpjm1                                 !  Vertical slab
315         !                                             ! ===============
316
317 
318         ! I. vertical trends associated with the lateral mixing
319         ! =====================================================
320         !  (excluding the vertical flux proportional to dk[t]
321
322
323         ! I.1 horizontal momentum gradient
324         ! --------------------------------
325
326         DO jk = 1, jpk
327            DO ji = 2, jpi
328               ! i-gradient of u at jj
329               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( ub(ji,jj  ,jk) - ub(ji-1,jj  ,jk) )
330               ! j-gradient of u and v at jj
331               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( ub(ji,jj+1,jk) - ub(ji  ,jj  ,jk) )
332               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( vb(ji,jj  ,jk) - vb(ji  ,jj-1,jk) )
333               ! j-gradient of u and v at jj+1
334               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( ub(ji,jj  ,jk) - ub(ji  ,jj-1,jk) )
335               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( vb(ji,jj+1,jk) - vb(ji  ,jj  ,jk) )
336            END DO
337         END DO
338         DO jk = 1, jpk
339            DO ji = 1, jpim1
340               ! i-gradient of v at jj
341               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( vb(ji+1,jj,jk) - vb(ji  ,jj  ,jk) )
342            END DO
343         END DO
344
345
346         ! I.2 Vertical fluxes
347         ! -------------------
348
349         ! Surface and bottom vertical fluxes set to zero
350         DO ji = 1, jpi
351            zfuw(ji, 1 ) = 0.e0
352            zfvw(ji, 1 ) = 0.e0
353            zfuw(ji,jpk) = 0.e0
354            zfvw(ji,jpk) = 0.e0
355         END DO
356
357         ! interior (2=<jk=<jpk-1) on U field
358         DO jk = 2, jpkm1
359            DO ji = 2, jpim1
360               zcoef0= 0.5 * aht0 * umask(ji,jj,jk)
361
362               zuwslpi = zcoef0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
363               zuwslpj = zcoef0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
364
365               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)   &
366                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ), 1. )
367               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1)   &
368                             + fmask(ji,jj-1,jk  )+fmask(ji,jj,jk  ), 1. )
369
370               zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi
371               zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj
372               ! vertical flux on u field
373               zfuw(ji,jk) = zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1)     &
374                                       +zdiu (ji,jk  ) + zdiu (ji+1,jk  ) )   &
375                           + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji  ,jk-1)     &
376                                       +zdj1u(ji,jk  ) + zdju (ji  ,jk  ) )
377               ! update avmu (add isopycnal vertical coefficient to avmu)
378               ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0
379               avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / aht0
380            END DO
381         END DO
382
383         ! interior (2=<jk=<jpk-1) on V field
384         DO jk = 2, jpkm1
385            DO ji = 2, jpim1
386               zcoef0= 0.5 * aht0 * vmask(ji,jj,jk)
387
388               zvwslpi = zcoef0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
389               zvwslpj = zcoef0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
390
391               zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)   &
392                             + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ), 1. )
393               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)   &
394                             + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ), 1. )
395
396               zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi
397               zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj
398               ! vertical flux on v field
399               zfvw(ji,jk) = zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1)     &
400                                       +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   &
401                           + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     &
402                                       +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) )
403               ! update avmv (add isopycnal vertical coefficient to avmv)
404               ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0
405               avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / aht0
406            END DO
407         END DO
408
409
410         ! I.3 Divergence of vertical fluxes added to the general tracer trend
411         ! -------------------------------------------------------------------
412         DO jk = 1, jpkm1
413            DO ji = 2, jpim1
414               ! volume elements
415               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
416               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
417               ! part of the k-component of isopycnal momentum diffusive trends
418               zuav = ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / zbu
419               zvav = ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / zbv
420               ! add the trends to the general trends
421               ua(ji,jj,jk) = ua(ji,jj,jk) + zuav
422               va(ji,jj,jk) = va(ji,jj,jk) + zvav
423            END DO
424         END DO
425         !                                             ! ===============
426      END DO                                           !   End of slab
427      !                                                ! ===============
428      CALL wrk_dealloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v ) 
429      !
430      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_iso')
431      !
432   END SUBROUTINE dyn_ldf_iso
433
434# else
435   !!----------------------------------------------------------------------
436   !!   Dummy module                           NO rotation of mixing tensor
437   !!----------------------------------------------------------------------
438CONTAINS
439   SUBROUTINE dyn_ldf_iso( kt )               ! Empty routine
440      INTEGER, INTENT(in) :: kt
441      WRITE(*,*) 'dyn_ldf_iso: You should not have seen this print! error?', kt
442   END SUBROUTINE dyn_ldf_iso
443#endif
444
445   !!======================================================================
446END MODULE dynldf_iso
Note: See TracBrowser for help on using the repository browser.