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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 2633

Last change on this file since 2633 was 2633, checked in by trackstand2, 13 years ago

Renamed wrk_use => wrk_in_use and wrk_release => wrk_not_released
  • Property svn:keywords set to Id
File size: 11.7 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      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
59      USE wrk_nemo, ONLY: zwi => wrk_3d_3 ! workspace
60      !!
61      INTEGER , INTENT( in ) ::   kt    ! ocean time-step index
62      REAL(wp), INTENT( in ) ::  p2dt   ! vertical profile of tracer time-step
63      !!
64      INTEGER  ::   ji, jj, jk             ! dummy loop indices
65      REAL(wp) ::   z1_p2dt, zcoef         ! temporary scalars
66      REAL(wp) ::   zzwi, zzws, zrhs       ! temporary scalars
67      !!----------------------------------------------------------------------
68
69      IF(wrk_in_use(3, 3))THEN
70         CALL ctl_stop('dyn_zdf_imp : requested workspace array unavailable.')
71         RETURN
72      END IF
73
74      IF( kt == nit000 ) THEN
75         IF(lwp) WRITE(numout,*)
76         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
77         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
78      ENDIF
79
80      ! 0. Local constant initialization
81      ! --------------------------------
82      z1_p2dt = 1._wp / p2dt      ! inverse of the timestep
83
84      ! 1. Vertical diffusion on u
85      ! ---------------------------
86      ! Matrix and second member construction
87      ! bottom boundary condition: both zwi and zws must be masked as avmu can take
88      ! non zero value at the ocean bottom depending on the bottom friction
89      ! used but the bottom velocities have already been updated with the bottom
90      ! friction velocity in dyn_bfr using values from the previous timestep. There
91      ! is no need to include these in the implicit calculation.
92      !
93      DO jk = 1, jpkm1        ! Matrix
94         DO jj = 2, jpjm1 
95            DO ji = fs_2, fs_jpim1   ! vector opt.
96               zcoef = - p2dt / fse3u(ji,jj,jk)
97               zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
98               zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk)
99               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
100               zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1)
101               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
102            END DO
103         END DO
104      END DO
105      DO jj = 2, jpjm1        ! Surface boudary conditions
106         DO ji = fs_2, fs_jpim1   ! vector opt.
107            zwi(ji,jj,1) = 0._wp
108            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
109         END DO
110      END DO
111
112      ! Matrix inversion starting from the first level
113      !-----------------------------------------------------------------------
114      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
115      !
116      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
117      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
118      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
119      !        (        ...               )( ...  ) ( ...  )
120      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
121      !
122      !   m is decomposed in the product of an upper and a lower triangular matrix
123      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
124      !   The solution (the after velocity) is in ua
125      !-----------------------------------------------------------------------
126      !
127      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
128         DO jj = 2, jpjm1   
129            DO ji = fs_2, fs_jpim1   ! vector opt.
130               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
131            END DO
132         END DO
133      END DO
134      !
135      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
136         DO ji = fs_2, fs_jpim1   ! vector opt.
137            ua(ji,jj,1) = ub(ji,jj,1) + p2dt * (  ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
138               &                                                       / ( fse3u(ji,jj,1) * rau0       )  )
139         END DO
140      END DO
141      DO jk = 2, jpkm1
142         DO jj = 2, jpjm1   
143            DO ji = fs_2, fs_jpim1   ! vector opt.
144               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side
145               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
146            END DO
147         END DO
148      END DO
149      !
150      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==
151         DO ji = fs_2, fs_jpim1   ! vector opt.
152            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
153         END DO
154      END DO
155      DO jk = jpk-2, 1, -1
156         DO jj = 2, jpjm1   
157            DO ji = fs_2, fs_jpim1   ! vector opt.
158               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
159            END DO
160         END DO
161      END DO
162
163      ! Normalization to obtain the general momentum trend ua
164      DO jk = 1, jpkm1
165         DO jj = 2, jpjm1   
166            DO ji = fs_2, fs_jpim1   ! vector opt.
167               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt
168            END DO
169         END DO
170      END DO
171
172
173      ! 2. Vertical diffusion on v
174      ! ---------------------------
175      ! Matrix and second member construction
176      ! bottom boundary condition: both zwi and zws must be masked as avmv can take
177      ! non zero value at the ocean bottom depending on the bottom friction
178      ! used but the bottom velocities have already been updated with the bottom
179      ! friction velocity in dyn_bfr using values from the previous timestep. There
180      ! is no need to include these in the implicit calculation.
181      !
182      DO jk = 1, jpkm1        ! Matrix
183         DO jj = 2, jpjm1   
184            DO ji = fs_2, fs_jpim1   ! vector opt.
185               zcoef = -p2dt / fse3v(ji,jj,jk)
186               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  )
187               zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk)
188               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
189               zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1)
190               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
191            END DO
192         END DO
193      END DO
194      DO jj = 2, jpjm1        ! Surface boudary conditions
195         DO ji = fs_2, fs_jpim1   ! vector opt.
196            zwi(ji,jj,1) = 0._wp
197            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
198         END DO
199      END DO
200
201      ! Matrix inversion
202      !-----------------------------------------------------------------------
203      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
204      !
205      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
206      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
207      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
208      !        (        ...               )( ...  ) ( ...  )
209      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
210      !
211      !   m is decomposed in the product of an upper and lower triangular matrix
212      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
213      !   The solution (after velocity) is in 2d array va
214      !-----------------------------------------------------------------------
215      !
216      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
217         DO jj = 2, jpjm1   
218            DO ji = fs_2, fs_jpim1   ! vector opt.
219               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
220            END DO
221         END DO
222      END DO
223      !
224      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
225         DO ji = fs_2, fs_jpim1   ! vector opt.
226            va(ji,jj,1) = vb(ji,jj,1) + p2dt * (  va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
227               &                                                       / ( fse3v(ji,jj,1) * rau0       )  )
228         END DO
229      END DO
230      DO jk = 2, jpkm1
231         DO jj = 2, jpjm1
232            DO ji = fs_2, fs_jpim1   ! vector opt.
233               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)   ! zrhs=right hand side
234               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
235            END DO
236         END DO
237      END DO
238      !
239      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==
240         DO ji = fs_2, fs_jpim1   ! vector opt.
241            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
242         END DO
243      END DO
244      DO jk = jpk-2, 1, -1
245         DO jj = 2, jpjm1   
246            DO ji = fs_2, fs_jpim1   ! vector opt.
247               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
248            END DO
249         END DO
250      END DO
251
252      ! Normalization to obtain the general momentum trend va
253      DO jk = 1, jpkm1
254         DO jj = 2, jpjm1   
255            DO ji = fs_2, fs_jpim1   ! vector opt.
256               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
257            END DO
258         END DO
259      END DO
260      !
261      IF(wrk_not_released(3, 3))THEN
262         CALL ctl_stop('dyn_zdf_imp : failed to release workspace array.')
263      END IF
264      !
265   END SUBROUTINE dyn_zdf_imp
266
267   !!==============================================================================
268END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.