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

source: branches/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 2789

Last change on this file since 2789 was 2789, checked in by cetlod, 13 years ago

Implementation of the merge of TRA/TRP : first guess, see ticket #842
  • 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     0.5  !  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 diffusion using a implicit time-stepping
14   !!----------------------------------------------------------------------
15   USE oce             ! ocean dynamics and tracers
16   USE dom_oce         ! ocean space and time domain
17   USE sbc_oce         ! surface boundary condition: ocean
18   USE zdf_oce         ! ocean vertical physics
19   USE phycst          ! physical constants
20   USE in_out_manager  ! I/O manager
21   USE lib_mpp         ! MPP library
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 wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
57      USE oce     , ONLY:   tsa             ! tsa used as 2 3D workspace
58      USE wrk_nemo, ONLY:   zwi => wrk_3d_3                 ! 3D 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) ::   z1_p2dt, zcoef, zzwi, zzws, zrhs   ! local scalars
65      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwd, zws
66      !!----------------------------------------------------------------------
67
68      IF( wrk_in_use(3, 3) ) THEN
69         CALL ctl_stop('dyn_zdf_imp: requested workspace array unavailable')   ;   RETURN
70      END IF
71      !
72      zwd => tsa(:,:,:,1) 
73      zws => tsa(:,:,:,2) 
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) )   CALL ctl_stop('dyn_zdf_imp: failed to release workspace array')
263      !
264   END SUBROUTINE dyn_zdf_imp
265
266   !!==============================================================================
267END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.