source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 2715

Last change on this file since 2715 was 2715, checked in by rblod, 10 years ago

First attempt to put dynamic allocation on the trunk

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