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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 2528

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago

Update NEMOGCM from branch nemo_v3_3_beta
  • Property svn:keywords set to Id
File size: 11.4 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 , NEMO Consortium (2010)
33   !! $Id$
34   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
35   !!----------------------------------------------------------------------
36CONTAINS
37
38   SUBROUTINE dyn_zdf_imp( kt, p2dt )
39      !!----------------------------------------------------------------------
40      !!                  ***  ROUTINE dyn_zdf_imp  ***
41      !!                   
42      !! ** Purpose :   Compute the trend due to the vert. momentum diffusion
43      !!      and the surface forcing, and add it to the general trend of
44      !!      the momentum equations.
45      !!
46      !! ** Method  :   The vertical momentum mixing trend is given by :
47      !!             dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
48      !!      backward time stepping
49      !!      Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2)
50      !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F)
51      !!      Add this trend to the general trend ua :
52      !!         ua = ua + dz( avmu dz(u) )
53      !!
54      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend.
55      !!---------------------------------------------------------------------
56      USE oce, ONLY :  zwd   => ta      ! use ta as workspace
57      USE oce, ONLY :  zws   => sa      ! use sa as workspace
58      !!
59      INTEGER , INTENT( in ) ::   kt    ! ocean time-step index
60      REAL(wp), INTENT( in ) ::  p2dt   ! vertical profile of tracer time-step
61      !!
62      INTEGER  ::   ji, jj, jk             ! dummy loop indices
63      REAL(wp) ::   z1_p2dt, zcoef         ! temporary scalars
64      REAL(wp) ::   zzwi, zzws, zrhs       ! temporary scalars
65      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi        ! 3D workspace
66      !!----------------------------------------------------------------------
67
68      IF( kt == nit000 ) THEN
69         IF(lwp) WRITE(numout,*)
70         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
71         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
72      ENDIF
73
74      ! 0. Local constant initialization
75      ! --------------------------------
76      z1_p2dt = 1._wp / p2dt      ! inverse of the timestep
77
78      ! 1. Vertical diffusion on u
79      ! ---------------------------
80      ! Matrix and second member construction
81      ! bottom boundary condition: both zwi and zws must be masked as avmu can take
82      ! non zero value at the ocean bottom depending on the bottom friction
83      ! used but the bottom velocities have already been updated with the bottom
84      ! friction velocity in dyn_bfr using values from the previous timestep. There
85      ! is no need to include these in the implicit calculation.
86      !
87      DO jk = 1, jpkm1        ! Matrix
88         DO jj = 2, jpjm1 
89            DO ji = fs_2, fs_jpim1   ! vector opt.
90               zcoef = - p2dt / fse3u(ji,jj,jk)
91               zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
92               zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk)
93               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
94               zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1)
95               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
96            END DO
97         END DO
98      END DO
99      DO jj = 2, jpjm1        ! Surface boudary conditions
100         DO ji = fs_2, fs_jpim1   ! vector opt.
101            zwi(ji,jj,1) = 0._wp
102            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
103         END DO
104      END DO
105
106      ! Matrix inversion starting from the first level
107      !-----------------------------------------------------------------------
108      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
109      !
110      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
111      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
112      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
113      !        (        ...               )( ...  ) ( ...  )
114      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
115      !
116      !   m is decomposed in the product of an upper and a lower triangular matrix
117      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
118      !   The solution (the after velocity) is in ua
119      !-----------------------------------------------------------------------
120      !
121      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
122         DO jj = 2, jpjm1   
123            DO ji = fs_2, fs_jpim1   ! vector opt.
124               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
125            END DO
126         END DO
127      END DO
128      !
129      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
130         DO ji = fs_2, fs_jpim1   ! vector opt.
131            ua(ji,jj,1) = ub(ji,jj,1) + p2dt * (  ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
132               &                                                       / ( fse3u(ji,jj,1) * rau0       )  )
133         END DO
134      END DO
135      DO jk = 2, jpkm1
136         DO jj = 2, jpjm1   
137            DO ji = fs_2, fs_jpim1   ! vector opt.
138               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side
139               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
140            END DO
141         END DO
142      END DO
143      !
144      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==
145         DO ji = fs_2, fs_jpim1   ! vector opt.
146            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
147         END DO
148      END DO
149      DO jk = jpk-2, 1, -1
150         DO jj = 2, jpjm1   
151            DO ji = fs_2, fs_jpim1   ! vector opt.
152               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
153            END DO
154         END DO
155      END DO
156
157      ! Normalization to obtain the general momentum trend ua
158      DO jk = 1, jpkm1
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) ) * z1_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._wp - 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._wp
191            zwd(ji,jj,1) = 1._wp - 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) + p2dt * (  va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
221               &                                                       / ( 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      ! Normalization to obtain the general momentum trend va
247      DO jk = 1, jpkm1
248         DO jj = 2, jpjm1   
249            DO ji = fs_2, fs_jpim1   ! vector opt.
250               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
251            END DO
252         END DO
253      END DO
254      !
255   END SUBROUTINE dyn_zdf_imp
256
257   !!==============================================================================
258END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.