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

source: branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90 @ 5989

Last change on this file since 5989 was 5989, checked in by deazer, 8 years ago

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