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

source: branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 2872

Last change on this file since 2872 was 2872, checked in by hliu, 13 years ago

1) added the semi-implicit bottom friction code (3 in DYN, 1 in ZDF), 2) added par_AMM7/12.h90, 3) added two arch files for NOCL mobius and ubuntu linux
  • 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     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
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 wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
58      USE oce     , ONLY:  zwd  => ta       , zws   => sa   ! (ta,sa) used as 3D workspace
59      USE wrk_nemo, ONLY:   zwi => wrk_3d_3                 ! 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) ::   bfr_imp = 1._wp
71      !!----------------------------------------------------------------------
72      !!----------------------------------------------------------------------
73
74      IF( wrk_in_use(3, 3) ) THEN
75         CALL ctl_stop('dyn_zdf_imp: requested workspace array unavailable')   ;   RETURN
76      END IF
77
78      IF( kt == nit000 ) THEN
79         IF(lwp) WRITE(numout,*)
80         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
81         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
82         IF(.NOT.ln_bfrimp) bfr_imp = 0._wp
83      ENDIF
84
85      ! 0. Local constant initialization
86      ! --------------------------------
87      z1_p2dt = 1._wp / p2dt      ! inverse of the timestep
88
89      ! 1. Vertical diffusion on u
90
91      ! Vertical diffusion on u&v
92      ! ---------------------------
93      ! Matrix and second member construction
94      !! bottom boundary condition: both zwi and zws must be masked as avmu can take
95      !! non zero value at the ocean bottom depending on the bottom friction
96      !! used but the bottom velocities have already been updated with the bottom
97      !! friction velocity in dyn_bfr using values from the previous timestep. There
98      !! is no need to include these in the implicit calculation.
99
100      ! The code has been modified here to implicitly implement bottom
101      ! friction: u(v)mask is not necessary here anymore.
102      ! H. Liu, April 2010.
103
104      ! 1. Vertical diffusion on u
105      DO jj = 2, jpjm1 
106         DO ji = fs_2, fs_jpim1   ! vector opt.
107            ikbum1 = mbku(ji,jj)
108            ! Apply stability criteria on absolute value  : Min abs(bfr) => Max (bfr)
109!               zbfru = MAX( bfrua(ji,jj), fse3u(ji,jj,ikbum1)*zinv )
110               zbfru = bfrua(ji,jj)
111
112            DO jk = 1, ikbum1
113               zcoef = - p2dt / fse3u(ji,jj,jk)
114               zwi(ji,jj,jk) = zcoef * avmu(ji,jj,jk  ) / fse3uw(ji,jj,jk  )
115               zws(ji,jj,jk) = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
116               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk)
117            END DO
118
119      ! Surface boudary conditions
120            zwi(ji,jj,1) = 0._wp
121            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
122
123      ! Bottom boudary conditions  ! H. Liu, May, 2010
124!            !commented out to be consisent with v3.2, h.liu
125!            z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 2.0 * bfr_imp
126            z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 1.0_wp * bfr_imp
127            zws(ji,jj,ikbum1) = 0._wp
128            zwd(ji,jj,ikbum1) = 1._wp - zwi(ji,jj,ikbum1) - z2dtf 
129
130      ! Matrix inversion starting from the first level
131      !-----------------------------------------------------------------------
132      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
133      !
134      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
135      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
136      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
137      !        (        ...               )( ...  ) ( ...  )
138      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
139      !
140      !   m is decomposed in the product of an upper and a lower triangular matrix
141      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
142      !   The solution (the after velocity) is in ua
143      !-----------------------------------------------------------------------
144
145      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
146            DO jk = 2, ikbum1
147               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
148            END DO
149
150      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
151            z2dtf = 0.5_wp * p2dt / ( fse3u(ji,jj,1) * rau0 )
152            ua(ji,jj,1) = ub(ji,jj,1) + p2dt * ua(ji,jj,1) + z2dtf * (utau_b(ji,jj) + utau(ji,jj))
153           DO jk = 2, ikbum1   
154               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side
155               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
156            END DO
157
158
159      ! thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk
160            ua(ji,jj,ikbum1) = ua(ji,jj,ikbum1) / zwd(ji,jj,ikbum1)
161            DO jk = ikbum1-1, 1, -1
162               ua(ji,jj,jk) =( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
163            END DO
164         END DO
165      END DO
166
167      ! Normalization to obtain the general momentum trend ua
168      DO jk = 1, jpkm1
169         DO jj = 2, jpjm1   
170            DO ji = fs_2, fs_jpim1   ! vector opt.
171               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt
172            END DO
173         END DO
174      END DO
175
176
177      ! 2. Vertical diffusion on v
178      ! ---------------------------
179
180      DO ji = fs_2, fs_jpim1   ! vector opt.
181         DO jj = 2, jpjm1   
182            ikbvm1 = mbkv(ji,jj)
183            ! Apply stability criteria on absolute value  : Min abs(bfr) => Max (bfr)
184!               zbfrv = MAX( bfrva(ji,jj), fse3v(ji,jj,ikbvm1)*zinv )
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 boudary conditions
195            zwi(ji,jj,1) = 0._wp
196            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
197
198      ! Bottom boudary conditions 
199      ! Bottom boudary conditions  ! H. Liu, May, 2010
200!            !commented out to be consisent with v3.2, h.liu
201!            z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 2.0 * bfr_imp
202            z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 1.0_wp * bfr_imp
203            zws(ji,jj,ikbvm1) = 0._wp
204            zwd(ji,jj,ikbvm1) = 1._wp - zwi(ji,jj,ikbvm1) - z2dtf 
205
206      ! Matrix inversion
207      !-----------------------------------------------------------------------
208      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
209      !
210      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
211      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
212      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
213      !        (        ...               )( ...  ) ( ...  )
214      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
215      !
216      !   m is decomposed in the product of an upper and lower triangular
217      !   matrix
218      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
219      !   The solution (after velocity) is in 2d array va
220      !-----------------------------------------------------------------------
221
222      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
223            DO jk = 2, ikbvm1
224               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
225            END DO
226
227      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
228            z2dtf = 0.5_wp * p2dt / ( fse3v(ji,jj,1)*rau0 )
229            va(ji,jj,1) = vb(ji,jj,1) + p2dt * va(ji,jj,1) + z2dtf * (vtau_b(ji,jj) + vtau(ji,jj))
230            DO jk = 2, ikbvm1
231               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)   ! zrhs=right hand side
232               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
233            END DO
234
235      ! thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk
236            va(ji,jj,ikbvm1) = va(ji,jj,ikbvm1) / zwd(ji,jj,ikbvm1)
237
238            DO jk = ikbvm1-1, 1, -1
239               va(ji,jj,jk) =( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
240            END DO
241
242         END DO
243      END DO
244
245      ! Normalization to obtain the general momentum trend va
246      DO jk = 1, jpkm1
247         DO jj = 2, jpjm1   
248            DO ji = fs_2, fs_jpim1   ! vector opt.
249               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
250            END DO
251         END DO
252      END DO
253      !
254      IF( wrk_not_released(3, 3) )   CALL ctl_stop('dyn_zdf_imp: failed to release workspace array')
255      !
256
257   END SUBROUTINE dyn_zdf_imp
258
259   !!==============================================================================
260END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.