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.
dynzdf.F90 in NEMO/trunk/src/OCE/DYN – NEMO

source: NEMO/trunk/src/OCE/DYN/dynzdf.F90 @ 14834

Last change on this file since 14834 was 14834, checked in by hadcv, 3 years ago

#2600: Merge in dev_r14273_HPC-02_Daley_Tiling
  • Property svn:keywords set to Id
File size: 25.4 KB
RevLine 
[456]1MODULE dynzdf
2   !!==============================================================================
3   !!                 ***  MODULE  dynzdf  ***
4   !! Ocean dynamics :  vertical component of the momentum mixing trend
5   !!==============================================================================
[2528]6   !! History :  1.0  !  2005-11  (G. Madec)  Original code
7   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
[9019]8   !!            4.0  !  2017-06  (G. Madec) remove the explicit time-stepping option + avm at t-point
[456]9   !!----------------------------------------------------------------------
[503]10
11   !!----------------------------------------------------------------------
[9019]12   !!   dyn_zdf       : compute the after velocity through implicit calculation of vertical mixing
[456]13   !!----------------------------------------------------------------------
[5836]14   USE oce            ! ocean dynamics and tracers variables
[9019]15   USE phycst         ! physical constants
[5836]16   USE dom_oce        ! ocean space and time domain variables
[9019]17   USE sbc_oce        ! surface boundary condition: ocean
[5836]18   USE zdf_oce        ! ocean vertical physics variables
[9019]19   USE zdfdrg         ! vertical physics: top/bottom drag coef.
20   USE dynadv    ,ONLY: ln_dynadv_vec    ! dynamics: advection form
[14834]21#if defined key_loop_fusion
22   USE dynldf_iso_lf,ONLY: akzu, akzv       ! dynamics: vertical component of rotated lateral mixing
23#else
[9019]24   USE dynldf_iso,ONLY: akzu, akzv       ! dynamics: vertical component of rotated lateral mixing
[14834]25#endif
[9490]26   USE ldfdyn         ! lateral diffusion: eddy viscosity coef. and type of operator
[5836]27   USE trd_oce        ! trends: ocean variables
28   USE trddyn         ! trend manager: dynamics
29   !
30   USE in_out_manager ! I/O manager
31   USE lib_mpp        ! MPP library
32   USE prtctl         ! Print control
33   USE timing         ! Timing
[456]34
35   IMPLICIT NONE
36   PRIVATE
37
[9019]38   PUBLIC   dyn_zdf   !  routine called by step.F90
[456]39
[9019]40   REAL(wp) ::  r_vvl     ! non-linear free surface indicator: =0 if ln_linssh=T, =1 otherwise
[456]41
42   !! * Substitutions
[12377]43#  include "do_loop_substitute.h90"
[13237]44#  include "domzgr_substitute.h90"
[456]45   !!----------------------------------------------------------------------
[9598]46   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1152]47   !! $Id$
[10068]48   !! Software governed by the CeCILL license (see ./LICENSE)
[456]49   !!----------------------------------------------------------------------
50CONTAINS
51   
[12377]52   SUBROUTINE dyn_zdf( kt, Kbb, Kmm, Krhs, puu, pvv, Kaa )
[456]53      !!----------------------------------------------------------------------
54      !!                  ***  ROUTINE dyn_zdf  ***
55      !!
[9019]56      !! ** Purpose :   compute the trend due to the vert. momentum diffusion
57      !!              together with the Leap-Frog time stepping using an
58      !!              implicit scheme.
59      !!
60      !! ** Method  :  - Leap-Frog time stepping on all trends but the vertical mixing
[12377]61      !!         u(after) =         u(before) + 2*dt *       u(rhs)                vector form or linear free surf.
[13237]62      !!         u(after) = ( e3u_b*u(before) + 2*dt * e3u_n*u(rhs) ) / e3u_after   otherwise
[9019]63      !!               - update the after velocity with the implicit vertical mixing.
64      !!      This requires to solver the following system:
[13237]65      !!         u(after) = u(after) + 1/e3u_after  dk+1[ mi(avm) / e3uw_after dk[ua] ]
[9019]66      !!      with the following surface/top/bottom boundary condition:
67      !!      surface: wind stress input (averaged over kt-1/2 & kt+1/2)
68      !!      top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90)
69      !!
[12377]70      !! ** Action :   (puu(:,:,:,Kaa),pvv(:,:,:,Kaa))   after velocity
[456]71      !!---------------------------------------------------------------------
[12377]72      INTEGER                             , INTENT( in )  ::  kt                  ! ocean time-step index
73      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs, Kaa ! ocean time level indices
74      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv            ! ocean velocities and RHS of momentum equation
[3294]75      !
[9019]76      INTEGER  ::   ji, jj, jk         ! dummy loop indices
77      INTEGER  ::   iku, ikv           ! local integers
78      REAL(wp) ::   zzwi, ze3ua, zdt   ! local scalars
79      REAL(wp) ::   zzws, ze3va        !   -      -
[10364]80      REAL(wp) ::   z1_e3ua, z1_e3va   !   -      -
81      REAL(wp) ::   zWu , zWv          !   -      -
82      REAL(wp) ::   zWui, zWvi         !   -      -
83      REAL(wp) ::   zWus, zWvs         !   -      -
[14834]84      REAL(wp), DIMENSION(A2D(nn_hls),jpk)        ::  zwi, zwd, zws   ! 3D workspace
[9019]85      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv   !  -      -
[456]86      !!---------------------------------------------------------------------
[3294]87      !
[9019]88      IF( ln_timing )   CALL timing_start('dyn_zdf')
[3294]89      !
[14834]90      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
91         IF( kt == nit000 ) THEN       !* initialization
92            IF(lwp) WRITE(numout,*)
93            IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
94            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
95            !
96            If( ln_linssh ) THEN   ;    r_vvl = 0._wp    ! non-linear free surface indicator
97            ELSE                   ;    r_vvl = 1._wp
98            ENDIF
[9019]99         ENDIF
100      ENDIF
[9250]101      !                             !* explicit top/bottom drag case
[12377]102      IF( .NOT.ln_drgimp )   CALL zdf_drg_exp( kt, Kmm, puu(:,:,:,Kbb), pvv(:,:,:,Kbb), puu(:,:,:,Krhs), pvv(:,:,:,Krhs) )  ! add top/bottom friction trend to (puu(Kaa),pvv(Kaa))
[9250]103      !
104      !
[9019]105      IF( l_trddyn )   THEN         !* temporary save of ta and sa trends
106         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 
[12377]107         ztrdu(:,:,:) = puu(:,:,:,Krhs)
108         ztrdv(:,:,:) = pvv(:,:,:,Krhs)
[456]109      ENDIF
[9019]110      !
[12377]111      !              !==  RHS: Leap-Frog time stepping on all trends but the vertical mixing  ==!   (put in puu(:,:,:,Kaa),pvv(:,:,:,Kaa))
[9019]112      !
113      !                    ! time stepping except vertical diffusion
114      IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
[13295]115         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[13286]116            puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kbb) + rDt * puu(ji,jj,jk,Krhs) ) * umask(ji,jj,jk)
117            pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kbb) + rDt * pvv(ji,jj,jk,Krhs) ) * vmask(ji,jj,jk)
118         END_3D
[9019]119      ELSE                                      ! applied on thickness weighted velocity
[13295]120         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[13286]121            puu(ji,jj,jk,Kaa) = (         e3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb )  &
122               &                  + rDt * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Krhs)  ) &
123               &                        / e3u(ji,jj,jk,Kaa) * umask(ji,jj,jk)
124            pvv(ji,jj,jk,Kaa) = (         e3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb )  &
125               &                  + rDt * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Krhs)  ) &
126               &                        / e3v(ji,jj,jk,Kaa) * vmask(ji,jj,jk)
127         END_3D
[9019]128      ENDIF
129      !                    ! add top/bottom friction
130      !     With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom.
131      !     J. Chanut: The bottom stress is computed considering after barotropic velocities, which does
132      !                not lead to the effective stress seen over the whole barotropic loop.
[12377]133      !     G. Madec : in linear free surface, e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) = e3u_0, so systematic use of e3u(:,:,:,Kaa)
[9019]134      IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
[13295]135         DO_3D( 0, 0, 0, 0, 1, jpkm1 )      ! remove barotropic velocities
[13286]136            puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - uu_b(ji,jj,Kaa) ) * umask(ji,jj,jk)
137            pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - vv_b(ji,jj,Kaa) ) * vmask(ji,jj,jk)
138         END_3D
[13497]139         DO_2D( 0, 0, 0, 0 )      ! Add bottom/top stress due to barotropic component only
[12377]140            iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points
141            ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
[13237]142            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm)    &
143               &             + r_vvl   * e3u(ji,jj,iku,Kaa)
144            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    &
145               &             + r_vvl   * e3v(ji,jj,ikv,Kaa)
[12489]146            puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua
147            pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va
[12377]148         END_2D
[13472]149         IF( ln_isfcav.OR.ln_drgice_imp ) THEN    ! Ocean cavities (ISF)
[13295]150            DO_2D( 0, 0, 0, 0 )
[12377]151               iku = miku(ji,jj)         ! top ocean level at u- and v-points
152               ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
[13237]153               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm)    &
154                  &             + r_vvl   * e3u(ji,jj,iku,Kaa)
155               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    &
156                  &             + r_vvl   * e3v(ji,jj,ikv,Kaa)
[12489]157               puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + rDt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua
158               pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va
[12377]159            END_2D
[9019]160         END IF
161      ENDIF
162      !
163      !              !==  Vertical diffusion on u  ==!
164      !
165      !                    !* Matrix construction
[12489]166      zdt = rDt * 0.5
[10364]167      IF( ln_zad_Aimp ) THEN   !!
168         SELECT CASE( nldf_dyn )
169         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
[13295]170            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[13237]171               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    &
172                  &             + r_vvl   * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
[12377]173               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
174                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
175               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
176                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
177               zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua
178               zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua
179               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 
180               zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
181               zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
182            END_3D
[10364]183         CASE DEFAULT               ! iso-level lateral mixing
[13295]184            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[13237]185               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    &    ! after scale factor at U-point
186                  &             + r_vvl   * e3u(ji,jj,jk,Kaa)
187               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) )   &
188                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
189               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) )   &
190                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
[12377]191               zWui = ( wi(ji,jj,jk  ) + wi(ji+1,jj,jk  ) ) / ze3ua
192               zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua
193               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp )
194               zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
195               zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
196            END_3D
[10364]197         END SELECT
[13497]198         DO_2D( 0, 0, 0, 0 )     !* Surface boundary conditions
[12377]199            zwi(ji,jj,1) = 0._wp
[13237]200            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm)    &
201               &             + r_vvl   * e3u(ji,jj,1,Kaa)
202            zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) )   &
203               &         / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2)
[12377]204            zWus = ( wi(ji  ,jj,2) +  wi(ji+1,jj,2) ) / ze3ua
205            zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp )
206            zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ) )
207         END_2D
[10364]208      ELSE
209         SELECT CASE( nldf_dyn )
210         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
[13295]211            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[13237]212               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    &
213                  &             + r_vvl   * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
[12377]214               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
215                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
216               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
217                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
218               zwi(ji,jj,jk) = zzwi
219               zws(ji,jj,jk) = zzws
220               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
221            END_3D
[10364]222         CASE DEFAULT               ! iso-level lateral mixing
[13295]223            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[13237]224               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm)    &
225                  &             + r_vvl   * e3u(ji,jj,jk,Kaa)   ! after scale factor at U-point
226               zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) )    &
227                  &         / ( ze3ua * e3uw(ji,jj,jk  ,Kmm) ) * wumask(ji,jj,jk  )
228               zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) )    &
229                  &         / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
[12377]230               zwi(ji,jj,jk) = zzwi
231               zws(ji,jj,jk) = zzws
232               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
233            END_3D
[10364]234         END SELECT
[13497]235         DO_2D( 0, 0, 0, 0 )     !* Surface boundary conditions
[12377]236            zwi(ji,jj,1) = 0._wp
237            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
238         END_2D
[10364]239      ENDIF
[9019]240      !
241      !
242      !              !==  Apply semi-implicit bottom friction  ==!
243      !
244      !     Only needed for semi-implicit bottom friction setup. The explicit
245      !     bottom friction has been included in "u(v)a" which act as the R.H.S
246      !     column vector of the tri-diagonal matrix equation
247      !
248      IF ( ln_drgimp ) THEN      ! implicit bottom friction
[13295]249         DO_2D( 0, 0, 0, 0 )
[12377]250            iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points
[13237]251            ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm)    &
252               &             + r_vvl   * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point
[12489]253            zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
[12377]254         END_2D
[13472]255         IF ( ln_isfcav.OR.ln_drgice_imp ) THEN   ! top friction (always implicit)
[13295]256            DO_2D( 0, 0, 0, 0 )
[12377]257               !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed
258               iku = miku(ji,jj)       ! ocean top level at u- and v-points
[13237]259               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm)    &
260                  &             + r_vvl   * e3u(ji,jj,iku,Kaa)   ! after scale factor at T-point
[12489]261               zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
[12377]262            END_2D
[9019]263         END IF
264      ENDIF
265      !
266      ! Matrix inversion starting from the first level
267      !-----------------------------------------------------------------------
268      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
269      !
270      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
271      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
272      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
273      !        (        ...               )( ...  ) ( ...  )
274      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
275      !
276      !   m is decomposed in the product of an upper and a lower triangular matrix
277      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
[12377]278      !   The solution (the after velocity) is in puu(:,:,:,Kaa)
[9019]279      !-----------------------------------------------------------------------
280      !
[13472]281      DO_3D( 0, 0, 0, 0, 2, jpkm1 )   !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
[12377]282         zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
283      END_3D
[9019]284      !
[13472]285      DO_2D( 0, 0, 0, 0 )             !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
[13237]286         ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm)    &
287            &             + r_vvl   * e3u(ji,jj,1,Kaa) 
[12489]288         puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
289            &                                      / ( ze3ua * rho0 ) * umask(ji,jj,1) 
[12377]290      END_2D
[13295]291      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
[12377]292         puu(ji,jj,jk,Kaa) = puu(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * puu(ji,jj,jk-1,Kaa)
293      END_3D
[9019]294      !
[13472]295      DO_2D( 0, 0, 0, 0 )             !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==!
[12377]296         puu(ji,jj,jpkm1,Kaa) = puu(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1)
297      END_2D
[13295]298      DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 )
[12377]299         puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - zws(ji,jj,jk) * puu(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk)
300      END_3D
[9019]301      !
302      !              !==  Vertical diffusion on v  ==!
303      !
304      !                       !* Matrix construction
[12489]305      zdt = rDt * 0.5
[10364]306      IF( ln_zad_Aimp ) THEN   !!
307         SELECT CASE( nldf_dyn )
308         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzv)
[13295]309            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[13237]310               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    &
311                  &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
[12377]312               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
313                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
314               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
315                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
316               zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va
317               zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va
318               zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp )
319               zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp )
320               zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
321            END_3D
[10364]322         CASE DEFAULT               ! iso-level lateral mixing
[13295]323            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[13237]324               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    &
325                  &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
326               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) )    &
327                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
328               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) )    &
329                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
[12377]330               zWvi = ( wi(ji,jj,jk  ) + wi(ji,jj+1,jk  ) ) / ze3va
331               zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va
332               zwi(ji,jj,jk) = zzwi  + zdt * MIN( zWvi, 0._wp )
333               zws(ji,jj,jk) = zzws  - zdt * MAX( zWvs, 0._wp )
334               zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
335            END_3D
[10364]336         END SELECT
[13472]337         DO_2D( 0, 0, 0, 0 )   !* Surface boundary conditions
[12377]338            zwi(ji,jj,1) = 0._wp
[13237]339            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm)    &
340               &             + r_vvl   * e3v(ji,jj,1,Kaa)
341            zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) )    &
342               &         / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2)
[12377]343            zWvs = ( wi(ji,jj  ,2) +  wi(ji,jj+1,2) ) / ze3va
344            zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp )
345            zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) )
346         END_2D
[10364]347      ELSE
348         SELECT CASE( nldf_dyn )
349         CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
[13295]350            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[13237]351               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    &
352                  &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
[12377]353               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
354                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
355               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
356                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
357               zwi(ji,jj,jk) = zzwi
358               zws(ji,jj,jk) = zzws
359               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
360            END_3D
[10364]361         CASE DEFAULT               ! iso-level lateral mixing
[13295]362            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[13237]363               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm)    &
364                  &             + r_vvl   * e3v(ji,jj,jk,Kaa)   ! after scale factor at V-point
365               zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) )    &
366                  &         / ( ze3va * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
367               zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) )    &
368                  &         / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
[12377]369               zwi(ji,jj,jk) = zzwi
370               zws(ji,jj,jk) = zzws
371               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
372            END_3D
[10364]373         END SELECT
[13497]374         DO_2D( 0, 0, 0, 0 )        !* Surface boundary conditions
[12377]375            zwi(ji,jj,1) = 0._wp
376            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
377         END_2D
[10364]378      ENDIF
[9019]379      !
380      !              !==  Apply semi-implicit top/bottom friction  ==!
381      !
382      !     Only needed for semi-implicit bottom friction setup. The explicit
383      !     bottom friction has been included in "u(v)a" which act as the R.H.S
384      !     column vector of the tri-diagonal matrix equation
385      !
386      IF( ln_drgimp ) THEN
[13295]387         DO_2D( 0, 0, 0, 0 )
[12377]388            ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
[13237]389            ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    &
390               &             + r_vvl   * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point
[12489]391            zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va           
[12377]392         END_2D
[13472]393         IF ( ln_isfcav.OR.ln_drgice_imp ) THEN
[13295]394            DO_2D( 0, 0, 0, 0 )
[12377]395               ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
[13237]396               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm)    &
397                  &             + r_vvl   * e3v(ji,jj,ikv,Kaa)   ! after scale factor at T-point
[12489]398               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va
[12377]399            END_2D
[9019]400         ENDIF
401      ENDIF
[456]402
[9019]403      ! Matrix inversion
404      !-----------------------------------------------------------------------
405      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
[503]406      !
[9019]407      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
408      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
409      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
410      !        (        ...               )( ...  ) ( ...  )
411      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
[503]412      !
[9019]413      !   m is decomposed in the product of an upper and lower triangular matrix
414      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
415      !   The solution (after velocity) is in 2d array va
416      !-----------------------------------------------------------------------
417      !
[13472]418      DO_3D( 0, 0, 0, 0, 2, jpkm1 )   !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
[12377]419         zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
420      END_3D
[9019]421      !
[13472]422      DO_2D( 0, 0, 0, 0 )             !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
[13237]423         ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm)    &
424            &             + r_vvl   * e3v(ji,jj,1,Kaa) 
[12489]425         pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
426            &                                      / ( ze3va * rho0 ) * vmask(ji,jj,1) 
[12377]427      END_2D
[13295]428      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
[12377]429         pvv(ji,jj,jk,Kaa) = pvv(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * pvv(ji,jj,jk-1,Kaa)
430      END_3D
[9019]431      !
[13472]432      DO_2D( 0, 0, 0, 0 )             !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==!
[12377]433         pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1)
434      END_2D
[13295]435      DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 )
[12377]436         pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - zws(ji,jj,jk) * pvv(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk)
437      END_3D
[9019]438      !
[503]439      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
[12489]440         ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) / rDt - ztrdu(:,:,:)
441         ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) / rDt - ztrdv(:,:,:)
[12377]442         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt, Kmm )
[9019]443         DEALLOCATE( ztrdu, ztrdv ) 
[456]444      ENDIF
445      !                                          ! print mean trends (used for debugging)
[12377]446      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' zdf  - Ua: ', mask1=umask,               &
447         &                                  tab3d_2=pvv(:,:,:,Kaa), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
[5836]448         !
[9019]449      IF( ln_timing )   CALL timing_stop('dyn_zdf')
[503]450      !
[456]451   END SUBROUTINE dyn_zdf
452
453   !!==============================================================================
454END MODULE dynzdf
Note: See TracBrowser for help on using the repository browser.