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/trunk/src/OCE/DYN – NEMO

source: NEMO/trunk/src/OCE/DYN/dynldf_iso.F90 @ 15795

Last change on this file since 15795 was 15094, checked in by clem, 3 years ago

repair isoneutral diffusion on momentum (ln_dynldf_iso=T)
  • Property svn:keywords set to Id
File size: 21.4 KB
RevLine 
[3]1MODULE dynldf_iso
2   !!======================================================================
3   !!                     ***  MODULE  dynldf_iso  ***
[5836]4   !! Ocean dynamics:   lateral viscosity trend (rotated laplacian operator)
[3]5   !!======================================================================
[2715]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
[5836]10   !!            3.7  !  2014-01  (F. Lemarie, G. Madec)  restructuration/simplification of ahm specification,
11   !!                 !                                   add velocity dependent coefficient and optional read in file
[2715]12   !!----------------------------------------------------------------------
[5836]13
[3]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
[5836]21   USE ldfdyn          ! lateral diffusion: eddy viscosity coef.
22   USE ldftra          ! lateral physics: eddy diffusivity
[3]23   USE zdf_oce         ! ocean vertical physics
24   USE ldfslp          ! iso-neutral slopes
[4990]25   !
[3]26   USE in_out_manager  ! I/O manager
[2715]27   USE lib_mpp         ! MPP library
[5836]28   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
[258]29   USE prtctl          ! Print control
[14834]30#if defined key_loop_fusion
31   USE dynldf_iso_lf, ONLY: dyn_ldf_iso_lf   ! lateral mixing - loop fusion version (dyn_ldf_iso routine )
32#endif
[3]33
34   IMPLICIT NONE
35   PRIVATE
36
[2715]37   PUBLIC   dyn_ldf_iso           ! called by step.F90
38   PUBLIC   dyn_ldf_iso_alloc     ! called by nemogcm.F90
[3]39
[9019]40   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akzu, akzv   !: vertical component of rotated lateral viscosity
[2715]41
[3]42   !! * Substitutions
[12377]43#  include "do_loop_substitute.h90"
[13237]44#  include "domzgr_substitute.h90"
[3]45   !!----------------------------------------------------------------------
[9598]46   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1156]47   !! $Id$
[10068]48   !! Software governed by the CeCILL license (see ./LICENSE)
[3]49   !!----------------------------------------------------------------------
50CONTAINS
51
[2715]52   INTEGER FUNCTION dyn_ldf_iso_alloc()
53      !!----------------------------------------------------------------------
54      !!                  ***  ROUTINE dyn_ldf_iso_alloc  ***
55      !!----------------------------------------------------------------------
[14834]56      dyn_ldf_iso_alloc = 0
57      IF( .NOT. ALLOCATED( akzu ) ) THEN
58         ALLOCATE( akzu(jpi,jpj,jpk), akzv(jpi,jpj,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      ENDIF
[2715]62   END FUNCTION dyn_ldf_iso_alloc
63
64
[12377]65   SUBROUTINE dyn_ldf_iso( kt, Kbb, Kmm, puu, pvv, Krhs )
[3]66      !!----------------------------------------------------------------------
67      !!                     ***  ROUTINE dyn_ldf_iso  ***
68      !!                       
[455]69      !! ** Purpose :   Compute the before trend of the rotated laplacian
70      !!      operator of lateral momentum diffusion except the diagonal
71      !!      vertical term that will be computed in dynzdf module. Add it
72      !!      to the general trend of momentum equation.
[3]73      !!
74      !! ** Method :
[455]75      !!        The momentum lateral diffusive trend is provided by a 2nd
76      !!      order operator rotated along neutral or geopotential surfaces
77      !!      (in s-coordinates).
[3]78      !!      It is computed using before fields (forward in time) and isopyc-
[455]79      !!      nal or geopotential slopes computed in routine ldfslp.
[3]80      !!      Here, u and v components are considered as 2 independent scalar
81      !!      fields. Therefore, the property of splitting divergent and rota-
82      !!      tional part of the flow of the standard, z-coordinate laplacian
83      !!      momentum diffusion is lost.
84      !!      horizontal fluxes associated with the rotated lateral mixing:
85      !!      u-component:
[12377]86      !!         ziut = ( ahmt + rn_ahm_b ) e2t * e3t / e1t  di[ uu ]
87      !!               -  ahmt              e2t * mi-1(uslp) dk[ mi(mk(uu)) ]
88      !!         zjuf = ( ahmf + rn_ahm_b ) e1f * e3f / e2f  dj[ uu ]
89      !!               -  ahmf              e1f * mi(vslp)   dk[ mj(mk(uu)) ]
[3]90      !!      v-component:
[12377]91      !!         zivf = ( ahmf + rn_ahm_b ) e2t * e3t / e1t  di[ vv ]
92      !!               -  ahmf              e2t * mj(uslp)   dk[ mi(mk(vv)) ]
93      !!         zjvt = ( ahmt + rn_ahm_b ) e1f * e3f / e2f  dj[ vv ]
94      !!               -  ahmt              e1f * mj-1(vslp) dk[ mj(mk(vv)) ]
[3]95      !!      take the horizontal divergence of the fluxes:
96      !!         diffu = 1/(e1u*e2u*e3u) {  di  [ ziut ] + dj-1[ zjuf ]  }
97      !!         diffv = 1/(e1v*e2v*e3v) {  di-1[ zivf ] + dj  [ zjvt ]  }
[12377]98      !!      Add this trend to the general trend (uu(rhs),vv(rhs)):
99      !!         uu(rhs) = uu(rhs) + diffu
[455]100      !!      CAUTION: here the isopycnal part is with a coeff. of aht. This
101      !!      should be modified for applications others than orca_r2 (!!bug)
[3]102      !!
103      !! ** Action :
[12377]104      !!       -(puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) updated with the before geopotential harmonic mixing trend
[9019]105      !!       -(akzu,akzv) to accompt for the diagonal vertical component
106      !!                    of the rotated operator in dynzdf module
[3]107      !!----------------------------------------------------------------------
[12377]108      INTEGER                             , INTENT( in )  ::  kt               ! ocean time-step index
109      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs   ! ocean time level indices
110      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv         ! ocean velocities and RHS of momentum equation
[2715]111      !
112      INTEGER  ::   ji, jj, jk   ! dummy loop indices
[9019]113      REAL(wp) ::   zabe1, zmskt, zmkt, zuav, zuwslpi, zuwslpj   ! local scalars
114      REAL(wp) ::   zabe2, zmskf, zmkf, zvav, zvwslpi, zvwslpj   !   -      -
[9490]115      REAL(wp) ::   zcof0, zcof1, zcof2, zcof3, zcof4, zaht_0    !   -      -
[14834]116      REAL(wp), DIMENSION(A2D(nn_hls))      ::   ziut, zivf, zdku, zdk1u  ! 2D workspace
117      REAL(wp), DIMENSION(A2D(nn_hls))      ::   zjuf, zjvt, zdkv, zdk1v  !  -      -
118      REAL(wp), DIMENSION(A1Di(nn_hls),jpk) ::   zfuw, zdiu, zdju, zdj1u  !  -      -
119      REAL(wp), DIMENSION(A1Di(nn_hls),jpk) ::   zfvw, zdiv, zdjv, zdj1v  !  -      -
[2715]120      !!----------------------------------------------------------------------
[3294]121      !
[14834]122#if defined key_loop_fusion
123      CALL dyn_ldf_iso_lf( kt, Kbb, Kmm, puu, pvv, Krhs    )
124#else
125
126      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
127         IF( kt == nit000 ) THEN
128            IF(lwp) WRITE(numout,*)
129            IF(lwp) WRITE(numout,*) 'dyn_ldf_iso : iso-neutral laplacian diffusive operator or '
130            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate horizontal diffusive operator'
131            !                                      ! allocate dyn_ldf_iso arrays
132            IF( dyn_ldf_iso_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_ldf_iso: failed to allocate arrays')
[15094]133            !
134            DO_2D_OVR( 0, 0, 0, 0 )
135               akzu(ji,jj,1)   = 0._wp
136               akzu(ji,jj,jpk) = 0._wp
137               akzv(ji,jj,1)   = 0._wp
138               akzv(ji,jj,jpk) = 0._wp
139            END_2D
140            !
[14834]141         ENDIF
[3]142      ENDIF
[216]143
[5836]144!!gm bug is dyn_ldf_iso called before tra_ldf_iso ....   <<<<<===== TO BE CHECKED
145      ! s-coordinate: Iso-level diffusion on momentum but not on tracer
[455]146      IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
[2715]147         !
[14834]148         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )      ! set the slopes of iso-level
[12377]149            uslp (ji,jj,jk) = - ( gdept(ji+1,jj,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
150            vslp (ji,jj,jk) = - ( gdept(ji,jj+1,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
151            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
152            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
153         END_3D
[455]154         ! Lateral boundary conditions on the slopes
[14834]155         IF (nn_hls == 1) 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 )
[9019]156         !
[14834]157      ENDIF
[9490]158         
159      zaht_0 = 0.5_wp * rn_Ud * rn_Ld                  ! aht_0 from namtra_ldf = zaht_max
160     
[3]161      !                                                ! ===============
162      DO jk = 1, jpkm1                                 ! Horizontal slab
163         !                                             ! ===============
164
165         ! Vertical u- and v-shears at level jk and jk+1
166         ! ---------------------------------------------
167         ! surface boundary condition: zdku(jk=1)=zdku(jk=2)
168         !                             zdkv(jk=1)=zdkv(jk=2)
169
[14834]170         DO_2D( 1, 1, 1, 1 )
171            zdk1u(ji,jj) = ( puu(ji,jj,jk,Kbb) -puu(ji,jj,jk+1,Kbb) ) * umask(ji,jj,jk+1)
172            zdk1v(ji,jj) = ( pvv(ji,jj,jk,Kbb) -pvv(ji,jj,jk+1,Kbb) ) * vmask(ji,jj,jk+1)
173         END_2D
[3]174
175         IF( jk == 1 ) THEN
176            zdku(:,:) = zdk1u(:,:)
177            zdkv(:,:) = zdk1v(:,:)
178         ELSE
[14834]179            DO_2D( 1, 1, 1, 1 )
180               zdku(ji,jj) = ( puu(ji,jj,jk-1,Kbb) - puu(ji,jj,jk,Kbb) ) * umask(ji,jj,jk)
181               zdkv(ji,jj) = ( pvv(ji,jj,jk-1,Kbb) - pvv(ji,jj,jk,Kbb) ) * vmask(ji,jj,jk)
182            END_2D
[3]183         ENDIF
184
185         !                               -----f-----
186         ! Horizontal fluxes on U             | 
187         ! --------------------===        t   u   t
188         !                                    | 
189         ! i-flux at t-point             -----f-----
190
[455]191         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
[14215]192            DO_2D( 0, 1, 0, 0 )
[13237]193               zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj)   &
194                  &    * MIN( e3u(ji  ,jj,jk,Kmm),                &
195                  &           e3u(ji-1,jj,jk,Kmm) ) * r1_e1t(ji,jj)
[3]196
[12377]197               zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)     &
198                  &                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ) , 1._wp )
[3]199
[12377]200               zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
201
202               ziut(ji,jj) = (  zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) )    &
203                  &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)      &
204                  &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
205            END_2D
[455]206         ELSE                   ! other coordinate system (zco or sco) : e3t
[14215]207            DO_2D( 0, 1, 0, 0 )
[13237]208               zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b )   &
209                  &     * e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e1t(ji,jj)
[3]210
[12377]211               zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  ) + umask(ji,jj,jk+1)     &
212                  &                 + umask(ji-1,jj,jk+1) + umask(ji,jj,jk  ) , 1._wp )
[455]213
[12377]214               zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
[455]215
[12377]216               ziut(ji,jj) = (  zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) )   &
217                  &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
218                  &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
219            END_2D
[455]220         ENDIF
[3]221
222         ! j-flux at f-point
[13295]223         DO_2D( 1, 0, 1, 0 )
[13237]224            zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b )   &
225               &     * e1f(ji,jj) * e3f(ji,jj,jk) * r1_e2f(ji,jj)
[3]226
[12377]227            zmskf = 1._wp / MAX(   umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)     &
228               &                 + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ) , 1._wp )
[3]229
[12377]230            zcof2 = - zaht_0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
[3]231
[12377]232            zjuf(ji,jj) = (  zabe2 * ( puu(ji,jj+1,jk,Kbb) - puu(ji,jj,jk,Kbb) )   &
233               &           + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
234               &                      +zdk1u(ji,jj+1) + zdku (ji,jj) )  ) * fmask(ji,jj,jk)
235         END_2D
[3]236
237         !                                |   t   |
238         ! Horizontal fluxes on V         |       |
239         ! --------------------===        f---v---f
240         !                                |       |
241         ! i-flux at f-point              |   t   |
242
[14215]243         DO_2D( 1, 0, 0, 0 )
[13237]244            zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b )   &
245               &     * e2f(ji,jj) * e3f(ji,jj,jk) * r1_e1f(ji,jj)
[3]246
[12377]247            zmskf = 1._wp / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)     &
248               &                + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ) , 1._wp )
[3]249
[12377]250            zcof1 = - zaht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
[3]251
[12377]252            zivf(ji,jj) = (  zabe1 * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji,jj,jk,Kbb) )    &
253               &           + zcof1 * (  zdkv (ji,jj) + zdk1v(ji+1,jj)      &
254               &                      + zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk)
255         END_2D
[3]256
257         ! j-flux at t-point
[455]258         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
[14215]259            DO_2D( 1, 0, 0, 1 )
[13237]260               zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj)   &
261                  &     * MIN( e3v(ji,jj  ,jk,Kmm),                 &
262                  &            e3v(ji,jj-1,jk,Kmm) ) * r1_e2t(ji,jj)
[3]263
[12377]264               zmskt = 1._wp / MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)     &
265                  &                + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ) , 1._wp )
[3]266
[12377]267               zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
[3]268
[12377]269               zjvt(ji,jj) = (  zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) )    &
270                  &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)      &
271                  &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
272            END_2D
[455]273         ELSE                   ! other coordinate system (zco or sco) : e3t
[14215]274            DO_2D( 1, 0, 0, 1 )
[13237]275               zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b )   &
276                  &     * e1t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e2t(ji,jj)
[3]277
[12377]278               zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
279                  &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
[3]280
[12377]281               zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
[455]282
[12377]283               zjvt(ji,jj) = (  zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) )   &
284                  &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     &
285                  &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
286            END_2D
[455]287         ENDIF
288
289
[3]290         ! Second derivative (divergence) and add to the general trend
291         ! -----------------------------------------------------------
[13497]292         DO_2D( 0, 0, 0, 0 )      !!gm Question vectop possible??? !!bug
[12377]293            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (  ziut(ji+1,jj) - ziut(ji,jj  )    &
[13237]294               &                           + zjuf(ji  ,jj) - zjuf(ji,jj-1)  ) * r1_e1e2u(ji,jj)   &
295               &                           / e3u(ji,jj,jk,Kmm)
[12377]296            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (  zivf(ji,jj  ) - zivf(ji-1,jj)    &
[13237]297               &                           + zjvt(ji,jj+1) - zjvt(ji,jj  )  ) * r1_e1e2v(ji,jj)   &
298               &                           / e3v(ji,jj,jk,Kmm)
[12377]299         END_2D
[3]300         !                                             ! ===============
301      END DO                                           !   End of slab
302      !                                                ! ===============
[216]303
[455]304      ! print sum trends (used for debugging)
[12377]305      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ldfh - Ua: ', mask1=umask, &
306         &                                  tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
[216]307
308
[455]309      !                                                ! ===============
[14834]310      DO jj = ntsj, ntej                               !  Vertical slab
[455]311         !                                             ! ===============
312
313 
314         ! I. vertical trends associated with the lateral mixing
315         ! =====================================================
316         !  (excluding the vertical flux proportional to dk[t]
317
318
319         ! I.1 horizontal momentum gradient
320         ! --------------------------------
321
322         DO jk = 1, jpk
[14834]323            DO ji = ntsi, ntei + nn_hls
[455]324               ! i-gradient of u at jj
[12377]325               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( puu(ji,jj  ,jk,Kbb) - puu(ji-1,jj  ,jk,Kbb) )
[455]326               ! j-gradient of u and v at jj
[12377]327               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( puu(ji,jj+1,jk,Kbb) - puu(ji  ,jj  ,jk,Kbb) )
328               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( pvv(ji,jj  ,jk,Kbb) - pvv(ji  ,jj-1,jk,Kbb) )
[455]329               ! j-gradient of u and v at jj+1
[12377]330               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( puu(ji,jj  ,jk,Kbb) - puu(ji  ,jj-1,jk,Kbb) )
331               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pvv(ji,jj+1,jk,Kbb) - pvv(ji  ,jj  ,jk,Kbb) )
[455]332            END DO
333         END DO
334         DO jk = 1, jpk
[14834]335            DO ji = ntsi - nn_hls, ntei
[455]336               ! i-gradient of v at jj
[12377]337               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji  ,jj  ,jk,Kbb) )
[455]338            END DO
339         END DO
340
341
342         ! I.2 Vertical fluxes
343         ! -------------------
344
345         ! Surface and bottom vertical fluxes set to zero
[14834]346         DO ji = ntsi - nn_hls, ntei + nn_hls
[455]347            zfuw(ji, 1 ) = 0.e0
348            zfvw(ji, 1 ) = 0.e0
349            zfuw(ji,jpk) = 0.e0
350            zfvw(ji,jpk) = 0.e0
351         END DO
352
353         ! interior (2=<jk=<jpk-1) on U field
354         DO jk = 2, jpkm1
[14834]355            DO ji = ntsi, ntei
[9490]356               zcof0 = 0.5_wp * zaht_0 * umask(ji,jj,jk)
[5836]357               !
[9019]358               zuwslpi = zcof0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
359               zuwslpj = zcof0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
[5836]360               !
[9019]361               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)      &
362                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ) , 1. )
363               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1) + fmask(ji,jj,jk-1)      &
364                             + fmask(ji,jj-1,jk  ) + fmask(ji,jj,jk  ) , 1. )
[455]365
[9019]366               zcof3 = - e2u(ji,jj) * zmkt * zuwslpi
367               zcof4 = - e1u(ji,jj) * zmkf * zuwslpj
[455]368               ! vertical flux on u field
[9019]369               zfuw(ji,jk) = zcof3 * (  zdiu (ji,jk-1) + zdiu (ji+1,jk-1)      &
370                  &                   + zdiu (ji,jk  ) + zdiu (ji+1,jk  )  )   &
371                  &        + zcof4 * (  zdj1u(ji,jk-1) + zdju (ji  ,jk-1)      &
372                  &                   + zdj1u(ji,jk  ) + zdju (ji  ,jk  )  )
373               ! vertical mixing coefficient (akzu)
[9490]374               ! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0
375               akzu(ji,jj,jk) = ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / zaht_0
[455]376            END DO
377         END DO
378
379         ! interior (2=<jk=<jpk-1) on V field
380         DO jk = 2, jpkm1
[14834]381            DO ji = ntsi, ntei
[9490]382               zcof0 = 0.5_wp * zaht_0 * vmask(ji,jj,jk)
[9019]383               !
384               zvwslpi = zcof0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
385               zvwslpj = zcof0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
386               !
387               zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)      &
388                  &          + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ) , 1. )
389               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)      &
390                  &          + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ) , 1. )
[455]391
[9019]392               zcof3 = - e2v(ji,jj) * zmkf * zvwslpi
393               zcof4 = - e1v(ji,jj) * zmkt * zvwslpj
[455]394               ! vertical flux on v field
[9019]395               zfvw(ji,jk) = zcof3 * (  zdiv (ji,jk-1) + zdiv (ji-1,jk-1)      &
396                  &                   + zdiv (ji,jk  ) + zdiv (ji-1,jk  )  )   &
397                  &        + zcof4 * (  zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)      &
398                  &                   + zdjv (ji,jk  ) + zdj1v(ji  ,jk  )  )
399               ! vertical mixing coefficient (akzv)
[9490]400               ! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0
401               akzv(ji,jj,jk) = ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / zaht_0
[455]402            END DO
403         END DO
404
405
406         ! I.3 Divergence of vertical fluxes added to the general tracer trend
407         ! -------------------------------------------------------------------
408         DO jk = 1, jpkm1
[14834]409            DO ji = ntsi, ntei
[13237]410               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj)   &
411                  &               / e3u(ji,jj,jk,Kmm)
412               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj)   &
413                  &               / e3v(ji,jj,jk,Kmm)
[455]414            END DO
415         END DO
416         !                                             ! ===============
417      END DO                                           !   End of slab
418      !                                                ! ===============
[14834]419#endif
[3]420   END SUBROUTINE dyn_ldf_iso
421
422   !!======================================================================
423END MODULE dynldf_iso
Note: See TracBrowser for help on using the repository browser.