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 NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/DYN – NEMO

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/DYN/dynldf_iso.F90 @ 14644

Last change on this file since 14644 was 14644, checked in by sparonuz, 3 years ago

Merge trunk -r14642:HEAD
  • Property svn:keywords set to Id
File size: 20.7 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
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   dyn_ldf_iso           ! called by step.F90
35   PUBLIC   dyn_ldf_iso_alloc     ! called by nemogcm.F90
36
37   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akzu, akzv   !: vertical component of rotated lateral viscosity
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 "do_loop_substitute.h90"
44#  include "domzgr_substitute.h90"
45#  include "single_precision_substitute.h90"
46   !!----------------------------------------------------------------------
47   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
48   !! $Id$
49   !! Software governed by the CeCILL license (see ./LICENSE)
50   !!----------------------------------------------------------------------
51CONTAINS
52
53   INTEGER FUNCTION dyn_ldf_iso_alloc()
54      !!----------------------------------------------------------------------
55      !!                  ***  ROUTINE dyn_ldf_iso_alloc  ***
56      !!----------------------------------------------------------------------
57      ALLOCATE( akzu(jpi,jpj,jpk) , zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) ,     & 
58         &      akzv(jpi,jpj,jpk) , 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, Kbb, Kmm, puu, pvv, Krhs )
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 + rn_ahm_b ) e2t * e3t / e1t  di[ uu ]
86      !!               -  ahmt              e2t * mi-1(uslp) dk[ mi(mk(uu)) ]
87      !!         zjuf = ( ahmf + rn_ahm_b ) e1f * e3f / e2f  dj[ uu ]
88      !!               -  ahmf              e1f * mi(vslp)   dk[ mj(mk(uu)) ]
89      !!      v-component:
90      !!         zivf = ( ahmf + rn_ahm_b ) e2t * e3t / e1t  di[ vv ]
91      !!               -  ahmf              e2t * mj(uslp)   dk[ mi(mk(vv)) ]
92      !!         zjvt = ( ahmt + rn_ahm_b ) e1f * e3f / e2f  dj[ vv ]
93      !!               -  ahmt              e1f * mj-1(vslp) dk[ mj(mk(vv)) ]
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 (uu(rhs),vv(rhs)):
98      !!         uu(rhs) = uu(rhs) + 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      !!       -(puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) updated with the before geopotential harmonic mixing trend
104      !!       -(akzu,akzv) 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      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs   ! ocean time level indices
109      REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv         ! ocean velocities and RHS of momentum equation
110      !
111      INTEGER  ::   ji, jj, jk   ! dummy loop indices
112      REAL(wp) ::   zabe1, zmskt, zmkt, zuav, zuwslpi, zuwslpj   ! local scalars
113      REAL(wp) ::   zabe2, zmskf, zmkf, zvav, zvwslpi, zvwslpj   !   -      -
114      REAL(wp) ::   zcof0, zcof1, zcof2, zcof3, zcof4, zaht_0    !   -      -
115      REAL(wp), DIMENSION(jpi,jpj) ::   ziut, zivf, zdku, zdk1u  ! 2D workspace
116      REAL(wp), DIMENSION(jpi,jpj) ::   zjuf, zjvt, zdkv, zdk1v  !  -      -
117      !!----------------------------------------------------------------------
118      !
119      IF( kt == nit000 ) THEN
120         IF(lwp) WRITE(numout,*)
121         IF(lwp) WRITE(numout,*) 'dyn_ldf_iso : iso-neutral laplacian diffusive operator or '
122         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate horizontal diffusive operator'
123         !                                      ! allocate dyn_ldf_bilap arrays
124         IF( dyn_ldf_iso_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_ldf_iso: failed to allocate arrays')
125      ENDIF
126
127!!gm bug is dyn_ldf_iso called before tra_ldf_iso ....   <<<<<===== TO BE CHECKED
128      ! s-coordinate: Iso-level diffusion on momentum but not on tracer
129      IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
130         !
131         DO_3D( 0, 0, 0, 0, 1, jpk )      ! set the slopes of iso-level
132            uslp (ji,jj,jk) = - ( gdept(ji+1,jj,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
133            vslp (ji,jj,jk) = - ( gdept(ji,jj+1,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
134            wslpi(ji,jj,jk) = - ( gdepw(ji+1,jj,jk,Kbb) - gdepw(ji-1,jj,jk,Kbb) ) * r1_e1t(ji,jj) * tmask(ji,jj,jk) * 0.5
135            wslpj(ji,jj,jk) = - ( gdepw(ji,jj+1,jk,Kbb) - gdepw(ji,jj-1,jk,Kbb) ) * r1_e2t(ji,jj) * tmask(ji,jj,jk) * 0.5
136         END_3D
137         ! Lateral boundary conditions on the slopes
138         CALL lbc_lnk( 'dynldf_iso', uslp , 'U', -1.0_wp, vslp , 'V', -1.0_wp, wslpi, 'W', -1.0_wp, wslpj, 'W', -1.0_wp )
139         !
140       ENDIF
141         
142      zaht_0 = 0.5_wp * rn_Ud * rn_Ld                  ! aht_0 from namtra_ldf = zaht_max
143     
144      !                                                ! ===============
145      DO jk = 1, jpkm1                                 ! Horizontal slab
146         !                                             ! ===============
147
148         ! Vertical u- and v-shears at level jk and jk+1
149         ! ---------------------------------------------
150         ! surface boundary condition: zdku(jk=1)=zdku(jk=2)
151         !                             zdkv(jk=1)=zdkv(jk=2)
152
153         zdk1u(:,:) = ( puu(:,:,jk,Kbb) -puu(:,:,jk+1,Kbb) ) * umask(:,:,jk+1)
154         zdk1v(:,:) = ( pvv(:,:,jk,Kbb) -pvv(:,:,jk+1,Kbb) ) * vmask(:,:,jk+1)
155
156         IF( jk == 1 ) THEN
157            zdku(:,:) = zdk1u(:,:)
158            zdkv(:,:) = zdk1v(:,:)
159         ELSE
160            zdku(:,:) = ( puu(:,:,jk-1,Kbb) - puu(:,:,jk,Kbb) ) * umask(:,:,jk)
161            zdkv(:,:) = ( pvv(:,:,jk-1,Kbb) - pvv(:,:,jk,Kbb) ) * vmask(:,:,jk)
162         ENDIF
163
164         !                               -----f-----
165         ! Horizontal fluxes on U             | 
166         ! --------------------===        t   u   t
167         !                                    | 
168         ! i-flux at t-point             -----f-----
169
170         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
171            DO_2D( 0, 1, 0, 0 )
172               zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj)   &
173                  &    * MIN( e3u(ji  ,jj,jk,Kmm),                &
174                  &           e3u(ji-1,jj,jk,Kmm) ) * r1_e1t(ji,jj)
175
176               zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)     &
177                  &                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ) , 1._wp )
178
179               zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
180
181               ziut(ji,jj) = (  zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) )    &
182                  &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)      &
183                  &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
184            END_2D
185         ELSE                   ! other coordinate system (zco or sco) : e3t
186            DO_2D( 0, 1, 0, 0 )
187               zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b )   &
188                  &     * e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e1t(ji,jj)
189
190               zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  ) + umask(ji,jj,jk+1)     &
191                  &                 + umask(ji-1,jj,jk+1) + umask(ji,jj,jk  ) , 1._wp )
192
193               zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
194
195               ziut(ji,jj) = (  zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) )   &
196                  &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
197                  &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
198            END_2D
199         ENDIF
200
201         ! j-flux at f-point
202         DO_2D( 1, 0, 1, 0 )
203            zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b )   &
204               &     * e1f(ji,jj) * e3f(ji,jj,jk) * r1_e2f(ji,jj)
205
206            zmskf = 1._wp / MAX(   umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)     &
207               &                 + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ) , 1._wp )
208
209            zcof2 = - zaht_0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
210
211            zjuf(ji,jj) = (  zabe2 * ( puu(ji,jj+1,jk,Kbb) - puu(ji,jj,jk,Kbb) )   &
212               &           + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
213               &                      +zdk1u(ji,jj+1) + zdku (ji,jj) )  ) * fmask(ji,jj,jk)
214         END_2D
215
216         !                                |   t   |
217         ! Horizontal fluxes on V         |       |
218         ! --------------------===        f---v---f
219         !                                |       |
220         ! i-flux at f-point              |   t   |
221
222         DO_2D( 1, 0, 0, 0 )
223            zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b )   &
224               &     * e2f(ji,jj) * e3f(ji,jj,jk) * r1_e1f(ji,jj)
225
226            zmskf = 1._wp / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)     &
227               &                + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ) , 1._wp )
228
229            zcof1 = - zaht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
230
231            zivf(ji,jj) = (  zabe1 * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji,jj,jk,Kbb) )    &
232               &           + zcof1 * (  zdkv (ji,jj) + zdk1v(ji+1,jj)      &
233               &                      + zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk)
234         END_2D
235
236         ! j-flux at t-point
237         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
238            DO_2D( 1, 0, 0, 1 )
239               zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj)   &
240                  &     * MIN( e3v(ji,jj  ,jk,Kmm),                 &
241                  &            e3v(ji,jj-1,jk,Kmm) ) * r1_e2t(ji,jj)
242
243               zmskt = 1._wp / MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)     &
244                  &                + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ) , 1._wp )
245
246               zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
247
248               zjvt(ji,jj) = (  zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) )    &
249                  &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)      &
250                  &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
251            END_2D
252         ELSE                   ! other coordinate system (zco or sco) : e3t
253            DO_2D( 1, 0, 0, 1 )
254               zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b )   &
255                  &     * e1t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e2t(ji,jj)
256
257               zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
258                  &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
259
260               zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
261
262               zjvt(ji,jj) = (  zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) )   &
263                  &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     &
264                  &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
265            END_2D
266         ENDIF
267
268
269         ! Second derivative (divergence) and add to the general trend
270         ! -----------------------------------------------------------
271         DO_2D( 0, 0, 0, 0 )      !!gm Question vectop possible??? !!bug
272            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (  ziut(ji+1,jj) - ziut(ji,jj  )    &
273               &                           + zjuf(ji  ,jj) - zjuf(ji,jj-1)  ) * r1_e1e2u(ji,jj)   &
274               &                           / e3u(ji,jj,jk,Kmm)
275            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (  zivf(ji,jj  ) - zivf(ji-1,jj)    &
276               &                           + zjvt(ji,jj+1) - zjvt(ji,jj  )  ) * r1_e1e2v(ji,jj)   &
277               &                           / e3v(ji,jj,jk,Kmm)
278         END_2D
279         !                                             ! ===============
280      END DO                                           !   End of slab
281      !                                                ! ===============
282
283      ! print sum trends (used for debugging)
284      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=CASTWP(puu(:,:,:,Krhs)), clinfo1=' ldfh - Ua: ', mask1=umask, &
285         &                                  tab3d_2=CASTWP(pvv(:,:,:,Krhs)), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
286
287
288      !                                                ! ===============
289      DO jj = 2, jpjm1                                 !  Vertical slab
290         !                                             ! ===============
291
292 
293         ! I. vertical trends associated with the lateral mixing
294         ! =====================================================
295         !  (excluding the vertical flux proportional to dk[t]
296
297
298         ! I.1 horizontal momentum gradient
299         ! --------------------------------
300
301         DO jk = 1, jpk
302            DO ji = 2, jpi
303               ! i-gradient of u at jj
304               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( puu(ji,jj  ,jk,Kbb) - puu(ji-1,jj  ,jk,Kbb) )
305               ! j-gradient of u and v at jj
306               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( puu(ji,jj+1,jk,Kbb) - puu(ji  ,jj  ,jk,Kbb) )
307               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( pvv(ji,jj  ,jk,Kbb) - pvv(ji  ,jj-1,jk,Kbb) )
308               ! j-gradient of u and v at jj+1
309               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( puu(ji,jj  ,jk,Kbb) - puu(ji  ,jj-1,jk,Kbb) )
310               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pvv(ji,jj+1,jk,Kbb) - pvv(ji  ,jj  ,jk,Kbb) )
311            END DO
312         END DO
313         DO jk = 1, jpk
314            DO ji = 1, jpim1
315               ! i-gradient of v at jj
316               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji  ,jj  ,jk,Kbb) )
317            END DO
318         END DO
319
320
321         ! I.2 Vertical fluxes
322         ! -------------------
323
324         ! Surface and bottom vertical fluxes set to zero
325         DO ji = 1, jpi
326            zfuw(ji, 1 ) = 0.e0
327            zfvw(ji, 1 ) = 0.e0
328            zfuw(ji,jpk) = 0.e0
329            zfvw(ji,jpk) = 0.e0
330         END DO
331
332         ! interior (2=<jk=<jpk-1) on U field
333         DO jk = 2, jpkm1
334            DO ji = 2, jpim1
335               zcof0 = 0.5_wp * zaht_0 * umask(ji,jj,jk)
336               !
337               zuwslpi = zcof0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
338               zuwslpj = zcof0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
339               !
340               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)      &
341                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ) , 1. )
342               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1) + fmask(ji,jj,jk-1)      &
343                             + fmask(ji,jj-1,jk  ) + fmask(ji,jj,jk  ) , 1. )
344
345               zcof3 = - e2u(ji,jj) * zmkt * zuwslpi
346               zcof4 = - e1u(ji,jj) * zmkf * zuwslpj
347               ! vertical flux on u field
348               zfuw(ji,jk) = zcof3 * (  zdiu (ji,jk-1) + zdiu (ji+1,jk-1)      &
349                  &                   + zdiu (ji,jk  ) + zdiu (ji+1,jk  )  )   &
350                  &        + zcof4 * (  zdj1u(ji,jk-1) + zdju (ji  ,jk-1)      &
351                  &                   + zdj1u(ji,jk  ) + zdju (ji  ,jk  )  )
352               ! vertical mixing coefficient (akzu)
353               ! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0
354               akzu(ji,jj,jk) = ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / zaht_0
355            END DO
356         END DO
357
358         ! interior (2=<jk=<jpk-1) on V field
359         DO jk = 2, jpkm1
360            DO ji = 2, jpim1
361               zcof0 = 0.5_wp * zaht_0 * vmask(ji,jj,jk)
362               !
363               zvwslpi = zcof0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
364               zvwslpj = zcof0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
365               !
366               zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)      &
367                  &          + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ) , 1. )
368               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)      &
369                  &          + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ) , 1. )
370
371               zcof3 = - e2v(ji,jj) * zmkf * zvwslpi
372               zcof4 = - e1v(ji,jj) * zmkt * zvwslpj
373               ! vertical flux on v field
374               zfvw(ji,jk) = zcof3 * (  zdiv (ji,jk-1) + zdiv (ji-1,jk-1)      &
375                  &                   + zdiv (ji,jk  ) + zdiv (ji-1,jk  )  )   &
376                  &        + zcof4 * (  zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)      &
377                  &                   + zdjv (ji,jk  ) + zdj1v(ji  ,jk  )  )
378               ! vertical mixing coefficient (akzv)
379               ! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0
380               akzv(ji,jj,jk) = ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / zaht_0
381            END DO
382         END DO
383
384
385         ! I.3 Divergence of vertical fluxes added to the general tracer trend
386         ! -------------------------------------------------------------------
387         DO jk = 1, jpkm1
388            DO ji = 2, jpim1
389               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj)   &
390                  &               / e3u(ji,jj,jk,Kmm)
391               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj)   &
392                  &               / e3v(ji,jj,jk,Kmm)
393            END DO
394         END DO
395         !                                             ! ===============
396      END DO                                           !   End of slab
397      !                                                ! ===============
398   END SUBROUTINE dyn_ldf_iso
399
400   !!======================================================================
401END MODULE dynldf_iso
402
Note: See TracBrowser for help on using the repository browser.