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

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 3186

Last change on this file since 3186 was 3186, checked in by smasson, 12 years ago

dev_NEMO_MERGE_2011: replace the old wrk_nemo with the new wrk_nemo
  • Property svn:keywords set to Id
File size: 11.1 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   USE zdfbfr          ! bottom friction setup
23   USE wrk_nemo        ! Memory Allocation
24   USE timing          ! Timing
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   dyn_zdf_imp   ! called by step.F90
30
31   !! * Substitutions
32#  include "domzgr_substitute.h90"
33#  include "vectopt_loop_substitute.h90"
34   !!----------------------------------------------------------------------
35   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
36   !! $Id$
37   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41   SUBROUTINE dyn_zdf_imp( kt, p2dt )
42      !!----------------------------------------------------------------------
43      !!                  ***  ROUTINE dyn_zdf_imp  ***
44      !!                   
45      !! ** Purpose :   Compute the trend due to the vert. momentum diffusion
46      !!      and the surface forcing, and add it to the general trend of
47      !!      the momentum equations.
48      !!
49      !! ** Method  :   The vertical momentum mixing trend is given by :
50      !!             dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
51      !!      backward time stepping
52      !!      Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2)
53      !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F)
54      !!      Add this trend to the general trend ua :
55      !!         ua = ua + dz( avmu dz(u) )
56      !!
57      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend.
58      !!---------------------------------------------------------------------
59      USE oce     , ONLY:   tsa       ! tsa used as 2 3D 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      INTEGER  ::   ikbum1, ikbvm1 ! local variable
66      REAL(wp) ::   z1_p2dt, z2dtf, zcoef, zzwi, zzws, zrhs ! local scalars
67
68      !! * Local variables for implicit bottom friction.    H. Liu
69      REAL(wp) ::   zbfru, zbfrv 
70      REAL(wp) ::   zbfr_imp = 0._wp                        ! toggle (SAVE'd by assignment)
71      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi
72      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwd, zws
73      !!----------------------------------------------------------------------
74      !
75      IF( nn_timing == 1 )  CALL timing_start('dyn_zdf_imp')
76      !
77      CALL wrk_alloc( jpi, jpj, jpk, zwi ) 
78      !
79      zwd => tsa(:,:,:,1) 
80      zws => tsa(:,:,:,2) 
81      !
82      IF( kt == nit000 ) THEN
83         IF(lwp) WRITE(numout,*)
84         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
85         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
86         IF(ln_bfrimp) zbfr_imp = 1._wp
87      ENDIF
88
89      ! 0. Local constant initialization
90      ! --------------------------------
91      z1_p2dt = 1._wp / p2dt      ! inverse of the timestep
92
93      ! 1. Vertical diffusion on u
94
95      ! Vertical diffusion on u&v
96      ! ---------------------------
97      ! Matrix and second member construction
98      !! bottom boundary condition: both zwi and zws must be masked as avmu can take
99      !! non zero value at the ocean bottom depending on the bottom friction
100      !! used but the bottom velocities have already been updated with the bottom
101      !! friction velocity in dyn_bfr using values from the previous timestep. There
102      !! is no need to include these in the implicit calculation.
103
104      ! The code has been modified here to implicitly implement bottom
105      ! friction: u(v)mask is not necessary here anymore.
106      ! H. Liu, April 2010.
107
108      ! 1. Vertical diffusion on u
109      DO jj = 2, jpjm1 
110         DO ji = fs_2, fs_jpim1   ! vector opt.
111            ikbum1 = mbku(ji,jj)
112               zbfru = bfrua(ji,jj)
113
114            DO jk = 1, ikbum1
115               zcoef = - p2dt / fse3u(ji,jj,jk)
116               zwi(ji,jj,jk) = zcoef * avmu(ji,jj,jk  ) / fse3uw(ji,jj,jk  )
117               zws(ji,jj,jk) = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
118               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk)
119            END DO
120
121      ! Surface boundary conditions
122            zwi(ji,jj,1) = 0._wp
123            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
124
125      ! Bottom boundary conditions  ! H. Liu, May, 2010
126!           !commented out to be consistent with v3.2, h.liu
127!           z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 2.0_wp * zbfr_imp
128            z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 1.0_wp * zbfr_imp
129            zws(ji,jj,ikbum1) = 0._wp
130            zwd(ji,jj,ikbum1) = 1._wp - zwi(ji,jj,ikbum1) - z2dtf 
131
132      ! Matrix inversion starting from the first level
133      !-----------------------------------------------------------------------
134      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
135      !
136      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
137      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
138      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
139      !        (        ...               )( ...  ) ( ...  )
140      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
141      !
142      !   m is decomposed in the product of an upper and a lower triangular matrix
143      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
144      !   The solution (the after velocity) is in ua
145      !-----------------------------------------------------------------------
146
147      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
148            DO jk = 2, ikbum1
149               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
150            END DO
151
152      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
153            z2dtf = 0.5_wp * p2dt / ( fse3u(ji,jj,1) * rau0 )
154            ua(ji,jj,1) = ub(ji,jj,1) + p2dt * ua(ji,jj,1) + z2dtf * (utau_b(ji,jj) + utau(ji,jj))
155           DO jk = 2, ikbum1   
156               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side
157               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
158            END DO
159
160
161      ! third recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk
162            ua(ji,jj,ikbum1) = ua(ji,jj,ikbum1) / zwd(ji,jj,ikbum1)
163            DO jk = ikbum1-1, 1, -1
164               ua(ji,jj,jk) =( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
165            END DO
166         END DO
167      END DO
168
169      ! Normalization to obtain the general momentum trend ua
170      DO jk = 1, jpkm1
171         DO jj = 2, jpjm1   
172            DO ji = fs_2, fs_jpim1   ! vector opt.
173               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt
174            END DO
175         END DO
176      END DO
177
178
179      ! 2. Vertical diffusion on v
180      ! ---------------------------
181
182      DO ji = fs_2, fs_jpim1   ! vector opt.
183         DO jj = 2, jpjm1   
184            ikbvm1 = mbkv(ji,jj)
185               zbfrv = bfrva(ji,jj)
186
187            DO jk = 1, ikbvm1
188               zcoef = -p2dt / fse3v(ji,jj,jk)
189               zwi(ji,jj,jk) = zcoef * avmv(ji,jj,jk  ) / fse3vw(ji,jj,jk  )
190               zws(ji,jj,jk) = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
191               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk)
192            END DO
193
194      ! Surface boundary conditions
195            zwi(ji,jj,1) = 0._wp
196            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
197
198      ! Bottom boundary conditions  ! H. Liu, May, 2010
199!           !commented out to be consistent with v3.2, h.liu
200!           z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 2.0_wp * zbfr_imp
201            z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 1.0_wp * zbfr_imp
202            zws(ji,jj,ikbvm1) = 0._wp
203            zwd(ji,jj,ikbvm1) = 1._wp - zwi(ji,jj,ikbvm1) - z2dtf 
204
205      ! Matrix inversion
206      !-----------------------------------------------------------------------
207      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
208      !
209      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
210      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
211      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
212      !        (        ...               )( ...  ) ( ...  )
213      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
214      !
215      !   m is decomposed in the product of an upper and lower triangular
216      !   matrix
217      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
218      !   The solution (after velocity) is in 2d array va
219      !-----------------------------------------------------------------------
220
221      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
222            DO jk = 2, ikbvm1
223               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
224            END DO
225
226      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
227            z2dtf = 0.5_wp * p2dt / ( fse3v(ji,jj,1)*rau0 )
228            va(ji,jj,1) = vb(ji,jj,1) + p2dt * va(ji,jj,1) + z2dtf * (vtau_b(ji,jj) + vtau(ji,jj))
229            DO jk = 2, ikbvm1
230               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)   ! zrhs=right hand side
231               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
232            END DO
233
234      ! third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk
235            va(ji,jj,ikbvm1) = va(ji,jj,ikbvm1) / zwd(ji,jj,ikbvm1)
236
237            DO jk = ikbvm1-1, 1, -1
238               va(ji,jj,jk) =( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
239            END DO
240
241         END DO
242      END DO
243
244      ! Normalization to obtain the general momentum trend va
245      DO jk = 1, jpkm1
246         DO jj = 2, jpjm1   
247            DO ji = fs_2, fs_jpim1   ! vector opt.
248               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
249            END DO
250         END DO
251      END DO
252      !
253      CALL wrk_dealloc( jpi, jpj, jpk, zwi ) 
254      !
255      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_imp')
256      !
257   END SUBROUTINE dyn_zdf_imp
258
259   !!==============================================================================
260END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.