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.
trazdf_imp.F90 in branches/dev_001_GM/NEMO/OPA_SRC/TRA – NEMO

source: branches/dev_001_GM/NEMO/OPA_SRC/TRA/trazdf_imp.F90 @ 786

Last change on this file since 786 was 786, checked in by gm, 16 years ago

dev_001_GM - merge TRC-TRA on OPA only, trabbl & zpshde not done and trdmld not OK - compilation OK

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.9 KB
Line 
1MODULE trazdf_imp
2   !!==============================================================================
3   !!                 ***  MODULE  trazdf_imp  ***
4   !! Ocean active tracers:  vertical component of the tracer mixing trend
5   !!==============================================================================
6   !! History :  OPA  !  1990-10  (B. Blanke)  Original code
7   !!            7.0  !  1991-11 (G. Madec)
8   !!             -   !  1992-06 (M. Imbard) correction on tracer trend loops
9   !!            8.0  !  1996-01 (G. Madec) statement function for e3
10   !!             -   !  1997-05 (G. Madec) vertical component of isopycnal
11   !!             -   !  1997-07 (G. Madec) geopotential diffusion in s-coord
12   !!             -   !  2000-08 (G. Madec) double diffusive mixing
13   !!   NEMO     1.0  !  2002-08 (G. Madec)  F90: Free form and module
14   !!            2.0  !  2006-11 (G. Madec)  New step reorganisation
15   !!            2.4  !  2008-01  (G. Madec) Merge TRA-TRC
16   !!----------------------------------------------------------------------
17
18   !!----------------------------------------------------------------------
19   !!   tra_zdf_imp : Update the tracer trend with the diagonal vertical 
20   !!                 part of the mixing tensor.
21   !!                 Vector optimization, use k-j-i loops.
22   !!----------------------------------------------------------------------
23   USE oce             ! ocean dynamics and tracers variables
24   USE dom_oce         ! ocean space and time domain variables
25   USE zdf_oce         ! ocean vertical physics variables
26   USE ldftra_oce      ! ocean active tracers: lateral physics
27   USE ldfslp          ! lateral physics: slope of diffusion
28   USE trdmod          ! ocean active tracers trends
29   USE trdmod_oce      ! ocean variables trends
30   USE zdfddm          ! ocean vertical physics: double diffusion
31   USE in_out_manager  ! I/O manager
32   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
33   USE prtctl          ! Print control
34   USE domvvl          ! variable volume
35
36   IMPLICIT NONE
37   PRIVATE
38
39   PUBLIC tra_zdf_imp   !  routine called by step.F90
40
41   !! * Substitutions
42#  include "domzgr_substitute.h90"
43#  include "ldftra_substitute.h90"
44#  include "zdfddm_substitute.h90"
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !!----------------------------------------------------------------------
48   !!  OPA 9.0 , LOCEAN-IPSL (2005)
49   !!----------------------------------------------------------------------
50CONTAINS
51   
52   SUBROUTINE tra_zdf_imp( kt, p2dt, cdtype )
53      !!----------------------------------------------------------------------
54      !!                  ***  ROUTINE tra_zdf_imp  ***
55      !!
56      !! ** Purpose :   Compute the trend due to the vertical tracer diffusion
57      !!     including the vertical component of lateral mixing (only for 2nd
58      !!     order operator, for fourth order it is already computed and add
59      !!     to the general trend in traldf.F) and add it to the general trend
60      !!     of the tracer equations.
61      !!
62      !! ** Method  :   The vertical component of the lateral diffusive trends
63      !!      is provided by a 2nd order operator rotated along neutral or geo-
64      !!      potential surfaces to which an eddy induced advection can be
65      !!      added. It is computed using before fields (forward in time) and
66      !!      isopycnal or geopotential slopes computed in routine ldfslp.
67      !!
68      !!      Second part: vertical trend associated with the vertical physics
69      !!      ===========  (including the vertical flux proportional to dk[t]
70      !!                  associated with the lateral mixing, through the
71      !!                  update of avt)
72      !!      The vertical diffusion of tracers (t & s) is given by:
73      !!             difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
74      !!      It is computed using a backward time scheme (t=ta).
75      !!      Surface and bottom boundary conditions: no diffusive flux on
76      !!      both tracers (bottom, applied through the masked field avt).
77      !!      Add this trend to the general trend ta,sa :
78      !!         ta = ta + dz( avt dz(t) )
79      !!         (sa = sa + dz( avs dz(t) ) if lk_zdfddm=T )
80      !!
81      !!      Third part: recover avt resulting from the vertical physics
82      !!      ==========  alone, for further diagnostics (for example to
83      !!                  compute the turbocline depth in zdfmxl.F90).
84      !!         avt = zavt
85      !!         (avs = zavs if lk_zdfddm=T )
86      !!
87      !! ** Action  : - Update (ta,sa) with before vertical diffusion trend
88      !!
89      !!---------------------------------------------------------------------
90      INTEGER         , INTENT(in   )                 ::   kt       ! ocean time-step index
91      CHARACTER(len=3), INTENT(in   )                 ::   cdtype   ! =TRA or TRC (tracer indicator)
92      REAL(wp)        , INTENT(in   ), DIMENSION(jpk) ::   p2dt     ! vertical profile of tracer time-step
93      !!
94      INTEGER  ::   ji, jj, jk                     ! dummy loop indices
95      REAL(wp) ::   znvvl                    ! temporary scalars
96      REAL(wp) ::   ze3tn, ze3ta, zvsfvvl   ! variable vertical scale factors
97      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwd, zws, zwi, zav   ! workspace arrays
98      !!---------------------------------------------------------------------
99
100      IF( kt == nit000 ) THEN
101         IF(lwp)WRITE(numout,*)
102         IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing (k-j-i loops)'
103         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ '
104         IF( cdtype /= 'TRA' ) THEN
105            IF(lwp)WRITE(numout,*) 'CAUTION TRC casei still not coded '
106         ENDIF
107      ENDIF
108
109      ! I. Local initialization
110      ! -----------------------
111      zwd  (1,:, : ) = 0.e0     ;     zwd  (jpi,:,:) = 0.e0
112      zws  (1,:, : ) = 0.e0     ;     zws  (jpi,:,:) = 0.e0
113      zwi  (1,:, : ) = 0.e0     ;     zwi  (jpi,:,:) = 0.e0
114      zav  (1,:, : ) = 0.e0     ;     zav  (jpi,:,:) = 0.e0
115      zav  (:,:,jpk) = 0.e0     ;     zav  ( : ,:,1) = 0.e0
116
117      ! I.1 Variable volume : to take into account vertical variable vertical scale factors
118      ! -------------------
119      IF( lk_vvl ) THEN   ;    znvvl = 1.
120      ELSE                ;    znvvl = 0.e0
121      ENDIF
122
123      ! II. Vertical trend associated with the vertical physics
124      ! =======================================================
125      !     (including the vertical flux proportional to dk[t] associated
126      !      with the lateral mixing, through the avt update)
127      !     dk[ avt dk[ (t,s) ] ] diffusive trends
128
129
130      ! II.0 Matrix construction
131      ! ------------------------
132
133      ! vertical mixing coef. put in zav
134      IF( ln_traldf_iso ) THEN      ! zav = avt + lateral mixing contribution
135#if defined key_ldfslp
136         DO jk = 2, jpkm1 
137            DO jj = 2, jpjm1
138               DO ji = fs_2, fs_jpim1   ! vector opt.
139                  zav(ji,jj,jk) = avt(ji,jj,jk) + fsahtw(ji,jj,jk) * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
140                     &                                                + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
141               END DO
142            END DO
143         END DO
144#endif
145      ELSE
146         zav(:,:,:) = avt(:,:,:)   ! zav = avt
147      ENDIF
148
149      ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
150      DO jk = 1, jpkm1
151         DO jj = 2, jpjm1
152            DO ji = fs_2, fs_jpim1   ! vector opt.
153               zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + ssha(ji,jj) * mut(ji,jj,jk) )
154               ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl                                ! after scale factor at T-point
155               ze3tn = ( 1. - znvvl )*fse3t(ji,jj,jk) + znvvl                        ! now   scale factor at T-point
156               zwi(ji,jj,jk) = - p2dt(jk) * zav(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
157               zws(ji,jj,jk) = - p2dt(jk) * zav(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
158               zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
159            END DO
160         END DO
161      END DO
162
163      ! Surface boudary conditions
164      DO jj = 2, jpjm1
165         DO ji = fs_2, fs_jpim1   ! vector opt.
166            zvsfvvl = fsve3t(ji,jj,1) * ( 1 + ssha(ji,jj) * mut(ji,jj,1) )
167            ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl                                   ! after scale factor at T-point
168            zwi(ji,jj,1) = 0.e0
169            zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
170         END DO
171      END DO
172
173
174      ! II.1. Vertical diffusion on t
175      ! ---------------------------
176      CALL zdf_imp( p2dt, zwd, zws, zwi, tb, ta )
177
178
179      ! II.2 Vertical diffusion on salinity
180      ! -----------------------------------
181
182!!gm    question: why surface value of zav is zero????  to be checked
183      IF( lk_zdfddm ) THEN
184         ! vertical mixing coef. on salintity including isopycnal slope (if necessary) put in zav
185         IF( ln_traldf_iso ) THEN      ! zav = avt + lateral mixing contribution
186#if defined key_ldfslp
187            DO jk = 2, jpkm1
188               DO jj = 2, jpjm1
189                  DO ji = fs_2, fs_jpim1   ! vector opt.
190                     zav(ji,jj,jk) = zav(ji,jj,jk) - avt(ji,jj,jk) + fsavs(ji,jj,jk) 
191                  END DO
192               END DO
193            END DO
194#endif
195         ELSE                             ! zav = avs
196            zav(:,:,:) = avs(:,:,:) 
197         ENDIF
198
199         ! Rebuild the Matrix as avt /= avs
200
201         ! Diagonal, inferior, superior  (including the bottom boundary condition via avs masked)
202         DO jk = 1, jpkm1
203            DO jj = 2, jpjm1
204               DO ji = fs_2, fs_jpim1   ! vector opt.
205                  zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + ssha(ji,jj) * mut(ji,jj,jk) )
206                  ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl                              ! after scale factor at T-point
207                  ze3tn = ( 1. - znvvl )*fse3t(ji,jj,jk) + znvvl                      ! now   scale factor at T-point
208                  zwi(ji,jj,jk) = - p2dt(jk) * zav(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
209                  zws(ji,jj,jk) = - p2dt(jk) * zav(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
210                  zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
211               END DO
212            END DO
213         END DO
214
215         ! Surface boudary conditions
216         DO jj = 2, jpjm1
217            DO ji = fs_2, fs_jpim1   ! vector opt.
218               zvsfvvl = fsve3t(ji,jj,1) * ( 1 + ssha(ji,jj) * mut(ji,jj,1) )
219               ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl                                  ! after scale factor at T-point
220               zwi(ji,jj,1) = 0.e0
221               zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
222            END DO
223         END DO
224         !
225      ENDIF
226
227      CALL zdf_imp( p2dt, zwd, zws, zwi, tb, ta )
228      !
229   END SUBROUTINE tra_zdf_imp
230
231
232   SUBROUTINE zdf_imp( p2dt, pwd, pws, pwi, ptb, pta )
233      !!----------------------------------------------------------------------
234      !!                  ***  ROUTINE tra_zdf_imp  ***
235      !!
236      !! ** Purpose :   Compute the trend due to the vertical tracer diffusion
237      !!     including the vertical component of lateral mixing (only for 2nd
238      !!     order operator, for fourth order it is already computed and add
239      !!     to the general trend in traldf.F) and add it to the general trend
240      !!     of the tracer equations.
241      !!
242      !! ** Method  :   The vertical component of the lateral diffusive trends
243      !!      is provided by a 2nd order operator rotated along neutral or geo-
244      !!      potential surfaces to which an eddy induced advection can be
245      !!      added. It is computed using before fields (forward in time) and
246      !!      isopycnal or geopotential slopes computed in routine ldfslp.
247      !!
248      !!      Second part: vertical trend associated with the vertical physics
249      !!      ===========  (including the vertical flux proportional to dk[t]
250      !!                  associated with the lateral mixing, through the
251      !!                  update of avt)
252      !!      The vertical diffusion of tracers (t & s) is given by:
253      !!             difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
254      !!      It is computed using a backward time scheme (t=ta).
255      !!      Surface and bottom boundary conditions: no diffusive flux on
256      !!      both tracers (bottom, applied through the masked field avt).
257      !!      Add this trend to the general trend ta,sa :
258      !!         ta = ta + dz( avt dz(t) )
259      !!         (sa = sa + dz( avs dz(t) ) if lk_zdfddm=T )
260      !!
261      !!      Third part: recover avt resulting from the vertical physics
262      !!      ==========  alone, for further diagnostics (for example to
263      !!                  compute the turbocline depth in zdfmxl.F90).
264      !!         avt = zavt
265      !!         (avs = zavs if lk_zdfddm=T )
266      !!
267      !! ** Action  : - Update (ta,sa) with before vertical diffusion trend
268      !!
269      !!---------------------------------------------------------------------
270      REAL(wp), INTENT(in   ), DIMENSION(jpk)         ::   p2dt     ! vertical profile of tracer time-step
271      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwd      ! matrix: diagnal terms
272      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pws      ! matrix: upper   terms
273      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwi      ! matrix: lower   terms
274      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   ptb      ! before tracer field
275      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta      ! in : tracer trend
276      !                                                                     ! out: after tracer
277      !!
278      INTEGER  ::   ji, jj, jk                     ! dummy loop indices
279      REAL(wp) ::   zrhs, znvvl              ! temporary scalars
280      REAL(wp) ::   ze3tb, ze3tn, zvsfvvl   ! variable vertical scale factors
281      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwt    ! workspace arrays
282      !!---------------------------------------------------------------------
283
284      ! Variable volume : to take into account vertical variable vertical scale factors
285      IF( lk_vvl ) THEN   ;    znvvl = 1.
286      ELSE                ;    znvvl = 0.e0
287      ENDIF
288
289      !! Matrix inversion from the first level
290      !!----------------------------------------------------------------------
291      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
292      !
293      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
294      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
295      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
296      !        (        ...               )( ...  ) ( ...  )
297      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
298      !
299      !   m is decomposed in the product of an upper and lower triangular matrix
300      !   The 3 diagonal terms are the input 3d arrays: zwd, zws, zwi
301      !   The second member and solution are put in pta array
302      !   The zwt array is a work space array
303
304      ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
305      DO jj = 2, jpjm1
306         DO ji = fs_2, fs_jpim1
307            zwt(ji,jj,1) = pwd(ji,jj,1)
308         END DO
309      END DO
310      DO jk = 2, jpkm1
311         DO jj = 2, jpjm1
312            DO ji = fs_2, fs_jpim1
313               zwt(ji,jj,jk) = pwd(ji,jj,jk) - pwi(ji,jj,jk) * pws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
314            END DO
315         END DO
316      END DO
317
318      ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
319      DO jj = 2, jpjm1
320         DO ji = fs_2, fs_jpim1
321            zvsfvvl = fsve3t(ji,jj,1) * ( 1 + sshb(ji,jj) * mut(ji,jj,1) )
322            ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl
323            ze3tn = ( 1. - znvvl ) + znvvl*fse3t (ji,jj,1)
324            pta(ji,jj,1) = ze3tb * ptb(ji,jj,1) + p2dt(1) * ze3tn * pta(ji,jj,1)
325         END DO
326      END DO
327      DO jk = 2, jpkm1
328         DO jj = 2, jpjm1
329            DO ji = fs_2, fs_jpim1
330               zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + sshb(ji,jj) * mut(ji,jj,jk) )
331               ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl
332               ze3tn = ( 1. - znvvl ) + znvvl*fse3t (ji,jj,jk)
333               zrhs = ze3tb * tb(ji,jj,jk) + p2dt(jk) * ze3tn * ta(ji,jj,jk)   ! zrhs=right hand side
334               pta(ji,jj,jk) = zrhs - pwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pta(ji,jj,jk-1)
335            END DO
336         END DO
337      END DO
338
339      ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
340      ! Save the masked temperature after in ta
341      ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
342      DO jj = 2, jpjm1
343         DO ji = fs_2, fs_jpim1
344            pta(ji,jj,jpkm1) = pta(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
345         END DO
346      END DO
347      DO jk = jpk-2, 1, -1
348         DO jj = 2, jpjm1
349            DO ji = fs_2, fs_jpim1
350               pta(ji,jj,jk) = ( pta(ji,jj,jk) - pws(ji,jj,jk) * ta(ji,jj,jk+1) )   &
351                  &          / zwt(ji,jj,jk) * tmask(ji,jj,jk)
352            END DO
353         END DO
354      END DO
355      !
356   END SUBROUTINE zdf_imp
357
358   !!==============================================================================
359END MODULE trazdf_imp
Note: See TracBrowser for help on using the repository browser.