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 @ 2636

Last change on this file since 2636 was 2636, checked in by gm, 13 years ago

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