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/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN – NEMO

source: branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 2148

Last change on this file since 2148 was 2148, checked in by cetlod, 14 years ago

merge LOCEAN 2010 developments branches
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 11.3 KB
Line 
1MODULE dynzdf_imp
2   !!==============================================================================
3   !!                    ***  MODULE  dynzdf_imp  ***
4   !! Ocean dynamics:  vertical component(s) of the momentum mixing trend
5   !!==============================================================================
6   !! History :  OPA  !  1990-10  (B. Blanke)  Original code
7   !!            8.0  !  1997-05  (G. Madec)  vertical component of isopycnal
8   !!   NEMO     1.0  !  2002-08  (G. Madec)  F90: Free form and module
9   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   dyn_zdf_imp  : update the momentum trend with the vertical diffu-
14   !!                  sion using a implicit time-stepping.
15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and tracers
17   USE dom_oce         ! ocean space and time domain
18   USE sbc_oce         ! surface boundary condition: ocean
19   USE zdf_oce         ! ocean vertical physics
20   USE phycst          ! physical constants
21   USE in_out_manager  ! I/O manager
22
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC   dyn_zdf_imp   ! called by step.F90
27
28   !! * Substitutions
29#  include "domzgr_substitute.h90"
30#  include "vectopt_loop_substitute.h90"
31   !!----------------------------------------------------------------------
32   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)
33   !! $Id$
34   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
35   !!----------------------------------------------------------------------
36
37CONTAINS
38
39   SUBROUTINE dyn_zdf_imp( kt, p2dt )
40      !!----------------------------------------------------------------------
41      !!                  ***  ROUTINE dyn_zdf_imp  ***
42      !!                   
43      !! ** Purpose :   Compute the trend due to the vert. momentum diffusion
44      !!      and the surface forcing, and add it to the general trend of
45      !!      the momentum equations.
46      !!
47      !! ** Method  :   The vertical momentum mixing trend is given by :
48      !!             dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
49      !!      backward time stepping
50      !!      Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2)
51      !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F)
52      !!      Add this trend to the general trend ua :
53      !!         ua = ua + dz( avmu dz(u) )
54      !!
55      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend.
56      !!---------------------------------------------------------------------
57      USE oce, ONLY :  zwd   => ta      ! use ta as workspace
58      USE oce, ONLY :  zws   => sa      ! use sa as workspace
59      !!
60      INTEGER , INTENT( in ) ::   kt    ! ocean time-step index
61      REAL(wp), INTENT( in ) ::  p2dt   ! vertical profile of tracer time-step
62      !!
63      INTEGER  ::   ji, jj, jk             ! dummy loop indices
64      REAL(wp) ::   zrau0r, zcoef         ! temporary scalars
65      REAL(wp) ::   zzwi, zzws, zrhs       ! temporary scalars
66      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi        ! 3D workspace
67      !!----------------------------------------------------------------------
68
69      IF( kt == nit000 ) THEN
70         IF(lwp) WRITE(numout,*)
71         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
72         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
73      ENDIF
74
75      ! 0. Local constant initialization
76      ! --------------------------------
77      zrau0r = 1. / rau0      ! inverse of the reference density
78
79      ! 1. Vertical diffusion on u
80      ! ---------------------------
81      ! Matrix and second member construction
82      ! bottom boundary condition: both zwi and zws must be masked as avmu can take
83      ! non zero value at the ocean bottom depending on the bottom friction
84      ! used but the bottom velocities have already been updated with the bottom
85      ! friction velocity in dyn_bfr using values from the previous timestep. There
86      ! is no need to include these in the implicit calculation.
87      !
88      DO jk = 1, jpkm1        ! Matrix
89         DO jj = 2, jpjm1 
90            DO ji = fs_2, fs_jpim1   ! vector opt.
91               zcoef = - p2dt / fse3u(ji,jj,jk)
92               zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
93               zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk)
94               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
95               zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1)
96               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws
97            END DO
98         END DO
99      END DO
100      DO jj = 2, jpjm1        ! Surface boudary conditions
101         DO ji = fs_2, fs_jpim1   ! vector opt.
102            zwi(ji,jj,1) = 0.
103            zwd(ji,jj,1) = 1. - zws(ji,jj,1)
104         END DO
105      END DO
106
107      ! Matrix inversion starting from the first level
108      !-----------------------------------------------------------------------
109      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
110      !
111      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
112      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
113      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
114      !        (        ...               )( ...  ) ( ...  )
115      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
116      !
117      !   m is decomposed in the product of an upper and a lower triangular matrix
118      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
119      !   The solution (the after velocity) is in ua
120      !-----------------------------------------------------------------------
121      !
122      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
123         DO jj = 2, jpjm1   
124            DO ji = fs_2, fs_jpim1   ! vector opt.
125               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
126            END DO
127         END DO
128      END DO
129      !
130      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
131         DO ji = fs_2, fs_jpim1   ! vector opt.
132            ua(ji,jj,1) = ub(ji,jj,1)   &
133               &        + p2dt * (  ua(ji,jj,1) + 0.5 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( fse3u(ji,jj,1) * rau0 )  )
134         END DO
135      END DO
136      DO jk = 2, jpkm1
137         DO jj = 2, jpjm1   
138            DO ji = fs_2, fs_jpim1   ! vector opt.
139               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side
140               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
141            END DO
142         END DO
143      END DO
144      !
145      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==
146         DO ji = fs_2, fs_jpim1   ! vector opt.
147            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
148         END DO
149      END DO
150      DO jk = jpk-2, 1, -1
151         DO jj = 2, jpjm1   
152            DO ji = fs_2, fs_jpim1   ! vector opt.
153               ua(ji,jj,jk) =( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
154            END DO
155         END DO
156      END DO
157      !
158      DO jk = 1, jpkm1        !==  Normalization to obtain the general momentum trend ua  ==
159         DO jj = 2, jpjm1   
160            DO ji = fs_2, fs_jpim1   ! vector opt.
161               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) / p2dt
162            END DO
163         END DO
164      END DO
165
166
167      ! 2. Vertical diffusion on v
168      ! ---------------------------
169      ! Matrix and second member construction
170      ! bottom boundary condition: both zwi and zws must be masked as avmv can take
171      ! non zero value at the ocean bottom depending on the bottom friction
172      ! used but the bottom velocities have already been updated with the bottom
173      ! friction velocity in dyn_bfr using values from the previous timestep. There
174      ! is no need to include these in the implicit calculation.
175      !
176      DO jk = 1, jpkm1        ! Matrix
177         DO jj = 2, jpjm1   
178            DO ji = fs_2, fs_jpim1   ! vector opt.
179               zcoef = -p2dt / fse3v(ji,jj,jk)
180               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  )
181               zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk)
182               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
183               zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1)
184               zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws
185            END DO
186         END DO
187      END DO
188      DO jj = 2, jpjm1        ! Surface boudary conditions
189         DO ji = fs_2, fs_jpim1   ! vector opt.
190            zwi(ji,jj,1) = 0.e0
191            zwd(ji,jj,1) = 1. - zws(ji,jj,1)
192         END DO
193      END DO
194
195      ! Matrix inversion
196      !-----------------------------------------------------------------------
197      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
198      !
199      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
200      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
201      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
202      !        (        ...               )( ...  ) ( ...  )
203      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
204      !
205      !   m is decomposed in the product of an upper and lower triangular matrix
206      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
207      !   The solution (after velocity) is in 2d array va
208      !-----------------------------------------------------------------------
209      !
210      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
211         DO jj = 2, jpjm1   
212            DO ji = fs_2, fs_jpim1   ! vector opt.
213               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
214            END DO
215         END DO
216      END DO
217      !
218      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
219         DO ji = fs_2, fs_jpim1   ! vector opt.
220            va(ji,jj,1) = vb(ji,jj,1)   &
221               &        + p2dt * (  va(ji,jj,1) + 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( fse3v(ji,jj,1) * rau0 )  )
222         END DO
223      END DO
224      DO jk = 2, jpkm1
225         DO jj = 2, jpjm1
226            DO ji = fs_2, fs_jpim1   ! vector opt.
227               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)   ! zrhs=right hand side
228               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
229            END DO
230         END DO
231      END DO
232      !
233      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
234         DO ji = fs_2, fs_jpim1   ! vector opt.
235            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
236         END DO
237      END DO
238      DO jk = jpk-2, 1, -1
239         DO jj = 2, jpjm1   
240            DO ji = fs_2, fs_jpim1   ! vector opt.
241               va(ji,jj,jk) =( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
242            END DO
243         END DO
244      END DO
245      !
246      DO jk = 1, jpkm1     !==  Normalization to obtain the general momentum trend va  ==
247         DO jj = 2, jpjm1   
248            DO ji = fs_2, fs_jpim1   ! vector opt.
249               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) / p2dt
250            END DO
251         END DO
252      END DO
253      !
254   END SUBROUTINE dyn_zdf_imp
255
256   !!==============================================================================
257END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.