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

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

dev_NEMO_MERGE_2011: supress TARGET attribute for tsa and use work arrays
  • Property svn:keywords set to Id
File size: 10.9 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      INTEGER , INTENT(in) ::  kt     ! ocean time-step index
60      REAL(wp), INTENT(in) ::  p2dt   ! vertical profile of tracer time-step
61      !!
62      INTEGER  ::   ji, jj, jk     ! dummy loop indices
63      INTEGER  ::   ikbum1, ikbvm1 ! local variable
64      REAL(wp) ::   z1_p2dt, z2dtf, zcoef, zzwi, zzws, zrhs ! local scalars
65
66      !! * Local variables for implicit bottom friction.    H. Liu
67      REAL(wp) ::   zbfru, zbfrv 
68      REAL(wp) ::   zbfr_imp = 0._wp                        ! toggle (SAVE'd by assignment)
69      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwd, zws
70      !!----------------------------------------------------------------------
71      !
72      IF( nn_timing == 1 )  CALL timing_start('dyn_zdf_imp')
73      !
74      CALL wrk_alloc( jpi,jpj,jpk, zwi, zwd, zws ) 
75      !
76      IF( kt == nit000 ) THEN
77         IF(lwp) WRITE(numout,*)
78         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
79         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
80         IF(ln_bfrimp) zbfr_imp = 1._wp
81      ENDIF
82
83      ! 0. Local constant initialization
84      ! --------------------------------
85      z1_p2dt = 1._wp / p2dt      ! inverse of the timestep
86
87      ! 1. Vertical diffusion on u
88
89      ! Vertical diffusion on u&v
90      ! ---------------------------
91      ! Matrix and second member construction
92      !! bottom boundary condition: both zwi and zws must be masked as avmu can take
93      !! non zero value at the ocean bottom depending on the bottom friction
94      !! used but the bottom velocities have already been updated with the bottom
95      !! friction velocity in dyn_bfr using values from the previous timestep. There
96      !! is no need to include these in the implicit calculation.
97
98      ! The code has been modified here to implicitly implement bottom
99      ! friction: u(v)mask is not necessary here anymore.
100      ! H. Liu, April 2010.
101
102      ! 1. Vertical diffusion on u
103      DO jj = 2, jpjm1 
104         DO ji = fs_2, fs_jpim1   ! vector opt.
105            ikbum1 = mbku(ji,jj)
106               zbfru = bfrua(ji,jj)
107
108            DO jk = 1, ikbum1
109               zcoef = - p2dt / fse3u(ji,jj,jk)
110               zwi(ji,jj,jk) = zcoef * avmu(ji,jj,jk  ) / fse3uw(ji,jj,jk  )
111               zws(ji,jj,jk) = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
112               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk)
113            END DO
114
115      ! Surface boundary conditions
116            zwi(ji,jj,1) = 0._wp
117            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
118
119      ! Bottom boundary conditions  ! H. Liu, May, 2010
120!           !commented out to be consistent with v3.2, h.liu
121!           z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 2.0_wp * zbfr_imp
122            z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 1.0_wp * zbfr_imp
123            zws(ji,jj,ikbum1) = 0._wp
124            zwd(ji,jj,ikbum1) = 1._wp - zwi(ji,jj,ikbum1) - z2dtf 
125
126      ! Matrix inversion starting from the first level
127      !-----------------------------------------------------------------------
128      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
129      !
130      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
131      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
132      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
133      !        (        ...               )( ...  ) ( ...  )
134      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
135      !
136      !   m is decomposed in the product of an upper and a lower triangular matrix
137      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
138      !   The solution (the after velocity) is in ua
139      !-----------------------------------------------------------------------
140
141      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
142            DO jk = 2, ikbum1
143               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
144            END DO
145
146      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
147            z2dtf = 0.5_wp * p2dt / ( fse3u(ji,jj,1) * rau0 )
148            ua(ji,jj,1) = ub(ji,jj,1) + p2dt * ua(ji,jj,1) + z2dtf * (utau_b(ji,jj) + utau(ji,jj))
149           DO jk = 2, ikbum1   
150               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side
151               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
152            END DO
153
154
155      ! third recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk
156            ua(ji,jj,ikbum1) = ua(ji,jj,ikbum1) / zwd(ji,jj,ikbum1)
157            DO jk = ikbum1-1, 1, -1
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
176      DO ji = fs_2, fs_jpim1   ! vector opt.
177         DO jj = 2, jpjm1   
178            ikbvm1 = mbkv(ji,jj)
179               zbfrv = bfrva(ji,jj)
180
181            DO jk = 1, ikbvm1
182               zcoef = -p2dt / fse3v(ji,jj,jk)
183               zwi(ji,jj,jk) = zcoef * avmv(ji,jj,jk  ) / fse3vw(ji,jj,jk  )
184               zws(ji,jj,jk) = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
185               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk)
186            END DO
187
188      ! Surface boundary conditions
189            zwi(ji,jj,1) = 0._wp
190            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
191
192      ! Bottom boundary conditions  ! H. Liu, May, 2010
193!           !commented out to be consistent with v3.2, h.liu
194!           z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 2.0_wp * zbfr_imp
195            z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 1.0_wp * zbfr_imp
196            zws(ji,jj,ikbvm1) = 0._wp
197            zwd(ji,jj,ikbvm1) = 1._wp - zwi(ji,jj,ikbvm1) - z2dtf 
198
199      ! Matrix inversion
200      !-----------------------------------------------------------------------
201      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
202      !
203      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
204      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
205      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
206      !        (        ...               )( ...  ) ( ...  )
207      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
208      !
209      !   m is decomposed in the product of an upper and lower triangular
210      !   matrix
211      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
212      !   The solution (after velocity) is in 2d array va
213      !-----------------------------------------------------------------------
214
215      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)
216            DO jk = 2, ikbvm1
217               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
218            END DO
219
220      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1
221            z2dtf = 0.5_wp * p2dt / ( fse3v(ji,jj,1)*rau0 )
222            va(ji,jj,1) = vb(ji,jj,1) + p2dt * va(ji,jj,1) + z2dtf * (vtau_b(ji,jj) + vtau(ji,jj))
223            DO jk = 2, ikbvm1
224               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)   ! zrhs=right hand side
225               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
226            END DO
227
228      ! third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk
229            va(ji,jj,ikbvm1) = va(ji,jj,ikbvm1) / zwd(ji,jj,ikbvm1)
230
231            DO jk = ikbvm1-1, 1, -1
232               va(ji,jj,jk) =( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
233            END DO
234
235         END DO
236      END DO
237
238      ! Normalization to obtain the general momentum trend va
239      DO jk = 1, jpkm1
240         DO jj = 2, jpjm1   
241            DO ji = fs_2, fs_jpim1   ! vector opt.
242               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
243            END DO
244         END DO
245      END DO
246      !
247      CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws ) 
248      !
249      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_imp')
250      !
251   END SUBROUTINE dyn_zdf_imp
252
253   !!==============================================================================
254END MODULE dynzdf_imp
Note: See TracBrowser for help on using the repository browser.