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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90 @ 2590

Last change on this file since 2590 was 2590, checked in by trackstand2, 13 years ago

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