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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90 @ 3270

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

Stephen Pickles, 11 Dec 2011

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