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_imp.F90 in branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 5104

Last change on this file since 5104 was 5104, checked in by mathiot, 9 years ago

correction of minor bug with isf + add some extra test on ln_isfcav + restore ssu/ssv

  • Property svn:keywords set to Id
File size: 17.3 KB
RevLine 
[3]1MODULE dynzdf_imp
[2715]2   !!======================================================================
[3]3   !!                    ***  MODULE  dynzdf_imp  ***
4   !! Ocean dynamics:  vertical component(s) of the momentum mixing trend
[2715]5   !!======================================================================
[2528]6   !! History :  OPA  !  1990-10  (B. Blanke)  Original code
7   !!            8.0  !  1997-05  (G. Madec)  vertical component of isopycnal
[2715]8   !!   NEMO     0.5  !  2002-08  (G. Madec)  F90: Free form and module
[2528]9   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
[3294]10   !!            3.4  !  2012-01  (H. Liu) Semi-implicit bottom friction
[503]11   !!----------------------------------------------------------------------
[3]12
13   !!----------------------------------------------------------------------
[2715]14   !!   dyn_zdf_imp  : update the momentum trend with the vertical diffusion using a implicit time-stepping
[3]15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and tracers
17   USE dom_oce         ! ocean space and time domain
[4292]18   USE domvvl          ! variable volume
[888]19   USE sbc_oce         ! surface boundary condition: ocean
20   USE zdf_oce         ! ocean vertical physics
[719]21   USE phycst          ! physical constants
[3]22   USE in_out_manager  ! I/O manager
[2715]23   USE lib_mpp         ! MPP library
[3294]24   USE zdfbfr          ! Bottom friction setup
25   USE wrk_nemo        ! Memory Allocation
26   USE timing          ! Timing
[4292]27   USE dynadv          ! dynamics: vector invariant versus flux form
[4354]28   USE dynspg_oce, ONLY: lk_dynspg_ts
[3]29
30   IMPLICIT NONE
31   PRIVATE
32
[2528]33   PUBLIC   dyn_zdf_imp   ! called by step.F90
[3]34
[4292]35   REAL(wp) ::  r_vvl     ! variable volume indicator, =1 if lk_vvl=T, =0 otherwise
36
[3]37   !! * Substitutions
38#  include "domzgr_substitute.h90"
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
[2528]41   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
[888]42   !! $Id$
[2528]43   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[3]44   !!----------------------------------------------------------------------
45CONTAINS
46
[503]47   SUBROUTINE dyn_zdf_imp( kt, p2dt )
[3]48      !!----------------------------------------------------------------------
49      !!                  ***  ROUTINE dyn_zdf_imp  ***
50      !!                   
51      !! ** Purpose :   Compute the trend due to the vert. momentum diffusion
52      !!      and the surface forcing, and add it to the general trend of
53      !!      the momentum equations.
54      !!
55      !! ** Method  :   The vertical momentum mixing trend is given by :
56      !!             dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
57      !!      backward time stepping
[2528]58      !!      Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2)
[3]59      !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F)
60      !!      Add this trend to the general trend ua :
61      !!         ua = ua + dz( avmu dz(u) )
62      !!
[2528]63      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend.
[3]64      !!---------------------------------------------------------------------
[3294]65      INTEGER , INTENT(in) ::  kt     ! ocean time-step index
[2715]66      REAL(wp), INTENT(in) ::  p2dt   ! vertical profile of tracer time-step
[2528]67      !!
[2715]68      INTEGER  ::   ji, jj, jk   ! dummy loop indices
[3294]69      INTEGER  ::   ikbu, ikbv   ! local integers
[2715]70      REAL(wp) ::   z1_p2dt, zcoef, zzwi, zzws, zrhs   ! local scalars
[4292]71      REAL(wp) ::   ze3ua, ze3va
[3294]72      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwd, zws
73      !!----------------------------------------------------------------------
74      !
75      IF( nn_timing == 1 )  CALL timing_start('dyn_zdf_imp')
76      !
77      CALL wrk_alloc( jpi,jpj,jpk, zwi, zwd, zws ) 
78      !
[3]79      IF( kt == nit000 ) THEN
80         IF(lwp) WRITE(numout,*)
81         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
82         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
[4292]83         !
84         IF( lk_vvl ) THEN   ;    r_vvl = 1._wp       ! Variable volume indicator
85         ELSE                ;    r_vvl = 0._wp       
86         ENDIF
[3]87      ENDIF
88
89      ! 0. Local constant initialization
90      ! --------------------------------
[2528]91      z1_p2dt = 1._wp / p2dt      ! inverse of the timestep
[455]92
[3294]93      ! 1. Apply semi-implicit bottom friction
94      ! --------------------------------------
95      ! Only needed for semi-implicit bottom friction setup. The explicit
96      ! bottom friction has been included in "u(v)a" which act as the R.H.S
97      ! column vector of the tri-diagonal matrix equation
98      !
99
100      IF( ln_bfrimp ) THEN
[4292]101         DO jj = 2, jpjm1
102            DO ji = 2, jpim1
103               ikbu = mbku(ji,jj)       ! ocean bottom level at u- and v-points
104               ikbv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
105               avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1)
106               avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1)
107            END DO
[3294]108         END DO
[5098]109         IF ( ln_isfcav ) THEN
110            DO jj = 2, jpjm1
111               DO ji = 2, jpim1
112                  ikbu = miku(ji,jj)       ! ocean top level at u- and v-points
113                  ikbv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
114                  IF (ikbu .GE. 2) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * fse3uw(ji,jj,ikbu)
115                  IF (ikbv .GE. 2) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * fse3vw(ji,jj,ikbv)
116               END DO
117            END DO
118         END IF
[3294]119      ENDIF
120
[4292]121#if defined key_dynspg_ts
122      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
123         DO jk = 1, jpkm1
124            ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk)
125            va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk)
126         END DO
127      ELSE                                            ! applied on thickness weighted velocity
128         DO jk = 1, jpkm1
129            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
130               &           + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
131               &                               / fse3u_a(:,:,jk) * umask(:,:,jk)
132            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
133               &           + p2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
134               &                               / fse3v_a(:,:,jk) * vmask(:,:,jk)
135         END DO
136      ENDIF
137
138      IF ( ln_bfrimp .AND.lk_dynspg_ts ) THEN
139         ! remove barotropic velocities:
140         DO jk = 1, jpkm1
141            ua(:,:,jk) = (ua(:,:,jk) - ua_b(:,:)) * umask(:,:,jk)
142            va(:,:,jk) = (va(:,:,jk) - va_b(:,:)) * vmask(:,:,jk)
[4990]143         END DO
144         ! Add bottom/top stress due to barotropic component only:
[4292]145         DO jj = 2, jpjm1       
146            DO ji = fs_2, fs_jpim1   ! vector opt.
147               ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
148               ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
149               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu)
150               ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl   * fse3v_a(ji,jj,ikbv)
151               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua
152               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va
153            END DO
154         END DO
[5104]155         IF ( ln_isfcav ) THEN
156            DO jj = 2, jpjm1       
157               DO ji = fs_2, fs_jpim1   ! vector opt.
158                  ikbu = miku(ji,jj)         ! top ocean level at u- and v-points
159                  ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
160                  ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu)
161                  ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl   * fse3v_a(ji,jj,ikbv)
162                  ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua
163                  va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va
164               END DO
165            END DO
166         END IF
[4292]167      ENDIF
168#endif
169
[3294]170      ! 2. Vertical diffusion on u
[3]171      ! ---------------------------
172      ! Matrix and second member construction
[1662]173      ! bottom boundary condition: both zwi and zws must be masked as avmu can take
[3294]174      ! non zero value at the ocean bottom depending on the bottom friction used.
[2528]175      !
176      DO jk = 1, jpkm1        ! Matrix
[3]177         DO jj = 2, jpjm1 
178            DO ji = fs_2, fs_jpim1   ! vector opt.
[4292]179               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,jk) + r_vvl   * fse3u_a(ji,jj,jk)   ! after scale factor at T-point
180               zcoef = - p2dt / ze3ua     
[5104]181               zzwi          = zcoef * avmu  (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
182               zwi(ji,jj,jk) = zzwi  * wumask(ji,jj,jk  )
183               zzws          = zcoef * avmu  (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 
[5098]184               zws(ji,jj,jk) = zzws  * wumask(ji,jj,jk+1)
185               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
[3]186            END DO
187         END DO
188      END DO
[4292]189      DO jj = 2, jpjm1        ! Surface boundary conditions
[3]190         DO ji = fs_2, fs_jpim1   ! vector opt.
[2528]191            zwi(ji,jj,1) = 0._wp
192            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
[3]193         END DO
194      END DO
195
196      ! Matrix inversion starting from the first level
197      !-----------------------------------------------------------------------
198      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
199      !
200      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
201      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
202      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
203      !        (        ...               )( ...  ) ( ...  )
204      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
205      !
206      !   m is decomposed in the product of an upper and a lower triangular matrix
207      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
208      !   The solution (the after velocity) is in ua
209      !-----------------------------------------------------------------------
[2528]210      !
[4990]211      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
[5098]212      DO jk = 2, jpkm1
213         DO jj = 2, jpjm1   
214            DO ji = fs_2, fs_jpim1   ! vector opt.
[3]215               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
216            END DO
217         END DO
218      END DO
[2528]219      !
220      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
[3]221         DO ji = fs_2, fs_jpim1   ! vector opt.
[4292]222#if defined key_dynspg_ts
[5098]223            ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl   * fse3u_a(ji,jj,1) 
224            ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
225               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
[4292]226#else
[5098]227            ua(ji,jj,1) = ub(ji,jj,1) &
228               &                   + p2dt *(ua(ji,jj,1) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
229               &                                      / ( fse3u(ji,jj,1) * rau0     ) * umask(ji,jj,1) ) 
[4292]230#endif
[5098]231         END DO
232      END DO
233      DO jk = 2, jpkm1
234         DO jj = 2, jpjm1
235            DO ji = fs_2, fs_jpim1
[4292]236#if defined key_dynspg_ts
237               zrhs = ua(ji,jj,jk)   ! zrhs=right hand side
238#else
239               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)
240#endif
[3]241               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
242            END DO
243         END DO
244      END DO
[2528]245      !
246      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==
[3]247         DO ji = fs_2, fs_jpim1   ! vector opt.
248            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
[5098]249         END DO
250      END DO
251      DO jk = jpk-2, 1, -1
252         DO jj = 2, jpjm1
253            DO ji = fs_2, fs_jpim1
[2528]254               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
[3]255            END DO
256         END DO
257      END DO
258
[4292]259#if ! defined key_dynspg_ts
[3]260      ! Normalization to obtain the general momentum trend ua
261      DO jk = 1, jpkm1
262         DO jj = 2, jpjm1   
263            DO ji = fs_2, fs_jpim1   ! vector opt.
[2528]264               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt
[3]265            END DO
266         END DO
267      END DO
[4292]268#endif
[3]269
[3294]270      ! 3. Vertical diffusion on v
[3]271      ! ---------------------------
272      ! Matrix and second member construction
[1662]273      ! bottom boundary condition: both zwi and zws must be masked as avmv can take
[3294]274      ! non zero value at the ocean bottom depending on the bottom friction used
[2528]275      !
276      DO jk = 1, jpkm1        ! Matrix
[3]277         DO jj = 2, jpjm1   
278            DO ji = fs_2, fs_jpim1   ! vector opt.
[4292]279               ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,jk)  + r_vvl * fse3v_a(ji,jj,jk)   ! after scale factor at T-point
280               zcoef = - p2dt / ze3va
[2528]281               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  )
[5098]282               zwi(ji,jj,jk) =  zzwi * wvmask(ji,jj,jk)
[2528]283               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
[5098]284               zws(ji,jj,jk) =  zzws * wvmask(ji,jj,jk+1)
285               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
[3]286            END DO
287         END DO
288      END DO
[4292]289      DO jj = 2, jpjm1        ! Surface boundary conditions
[3]290         DO ji = fs_2, fs_jpim1   ! vector opt.
[2528]291            zwi(ji,jj,1) = 0._wp
292            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
[3]293         END DO
294      END DO
295
296      ! Matrix inversion
297      !-----------------------------------------------------------------------
298      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
299      !
300      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
301      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
302      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
303      !        (        ...               )( ...  ) ( ...  )
304      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
305      !
[2528]306      !   m is decomposed in the product of an upper and lower triangular matrix
[3]307      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
308      !   The solution (after velocity) is in 2d array va
309      !-----------------------------------------------------------------------
[2528]310      !
[4990]311      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
[5098]312      DO jk = 2, jpkm1       
313         DO jj = 2, jpjm1   
314            DO ji = fs_2, fs_jpim1   ! vector opt.
[3]315               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
316            END DO
317         END DO
318      END DO
[2528]319      !
320      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
[3]321         DO ji = fs_2, fs_jpim1   ! vector opt.
[4292]322#if defined key_dynspg_ts           
[5098]323            ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl   * fse3v_a(ji,jj,1) 
324            va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
[4292]325               &                                      / ( ze3va * rau0 ) 
326#else
[5098]327            va(ji,jj,1) = vb(ji,jj,1) &
328               &                   + p2dt *(va(ji,jj,1) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
329               &                                                       / ( fse3v(ji,jj,1) * rau0     )  )
[4292]330#endif
[5098]331         END DO
332      END DO
333      DO jk = 2, jpkm1
334         DO jj = 2, jpjm1
335            DO ji = fs_2, fs_jpim1   ! vector opt.
[4292]336#if defined key_dynspg_ts
337               zrhs = va(ji,jj,jk)   ! zrhs=right hand side
338#else
339               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)
340#endif
[3]341               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
342            END DO
343         END DO
344      END DO
[2528]345      !
[4292]346      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
[3]347         DO ji = fs_2, fs_jpim1   ! vector opt.
348            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
[5098]349         END DO
350      END DO
351      DO jk = jpk-2, 1, -1
352         DO jj = 2, jpjm1
353            DO ji = fs_2, fs_jpim1
[2528]354               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
[3]355            END DO
356         END DO
357      END DO
358
359      ! Normalization to obtain the general momentum trend va
[4292]360#if ! defined key_dynspg_ts
[3]361      DO jk = 1, jpkm1
362         DO jj = 2, jpjm1   
363            DO ji = fs_2, fs_jpim1   ! vector opt.
[2528]364               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
[3]365            END DO
366         END DO
367      END DO
[4292]368#endif
[3294]369
[4292]370      ! J. Chanut: Lines below are useless ?
[3294]371      !! restore bottom layer avmu(v)
372      IF( ln_bfrimp ) THEN
[4990]373        DO jj = 2, jpjm1
374           DO ji = 2, jpim1
375              ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
376              ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
377              avmu(ji,jj,ikbu+1) = 0.e0
378              avmv(ji,jj,ikbv+1) = 0.e0
379           END DO
380        END DO
[5098]381        IF (ln_isfcav) THEN
382           DO jj = 2, jpjm1
383              DO ji = 2, jpim1
384                 ikbu = miku(ji,jj)         ! ocean top level at u- and v-points
385                 ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
386                 IF (ikbu > 1) avmu(ji,jj,ikbu) = 0.e0
387                 IF (ikbv > 1) avmv(ji,jj,ikbv) = 0.e0
388              END DO
389           END DO
390        END IF
[3294]391      ENDIF
[2528]392      !
[3294]393      CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws) 
[2715]394      !
[3294]395      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_imp')
396      !
[3]397   END SUBROUTINE dyn_zdf_imp
398
399   !!==============================================================================
400END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.