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_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA – NEMO

source: branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trazdf_imp.F90 @ 2024

Last change on this file since 2024 was 2024, checked in by cetlod, 14 years ago

Merge of active and passive tracer advection/diffusion modules, see ticket:693
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 14.8 KB
Line 
1MODULE trazdf_imp
2   !!======================================================================
3   !!                 ***  MODULE  trazdf_imp  ***
4   !! Ocean  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   !!                 !  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   !!            3.2  !  2009-03  (G. Madec)  heat and salt content trends
16   !!            3.3  !  2010-06  (C. Ethe, G. Madec) Merge TRA-TRC
17   !!----------------------------------------------------------------------
18 
19   !!----------------------------------------------------------------------
20   !!   tra_zdf_imp : Update the tracer trend with the diagonal vertical 
21   !!                 part of the mixing tensor.
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 zdfddm          ! ocean vertical physics: double diffusion
29   USE in_out_manager  ! I/O manager
30   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
31   USE domvvl          ! variable volume
32   USE ldftra          ! lateral mixing type
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   tra_zdf_imp   !  routine called by step.F90
38
39   !! * Substitutions
40#  include "domzgr_substitute.h90"
41#  include "ldftra_substitute.h90"
42#  include "zdfddm_substitute.h90"
43#  include "vectopt_loop_substitute.h90"
44   !!----------------------------------------------------------------------
45   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
46   !! $Id$
47   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
48   !!----------------------------------------------------------------------
49CONTAINS
50 
51   SUBROUTINE tra_zdf_imp( kt    , cdtype, p2dt,    &
52      &                    ptrab , ptraa , kjpt     )
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 the tracer t  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      !!         if lk_zdfddm=T, use avs for salinity or for passive tracers
80      !!         (sa = sa + dz( avs dz(t) )
81      !!
82      !!      Third part: recover avt resulting from the vertical physics
83      !!      ==========  alone, for further diagnostics (for example to
84      !!                  compute the turbocline depth in zdfmxl.F90).
85      !!         if lk_zdfddm=T, use avt = zavt
86      !!         (avs = zavs if lk_zdfddm=T )
87      !!
88      !! ** Action  : - Update (ta) with before vertical diffusion trend
89      !!
90      !!---------------------------------------------------------------------
91      !! * Modules used
92      USE oce    , ONLY :   zwd   => ua   ! ua used as workspace
93      USE oce    , ONLY :   zws   => va   ! va  -          -
94      !! * Arguments
95      INTEGER         , INTENT(in   )                                ::   kt             ! ocean time-step index
96      CHARACTER(len=3), INTENT(in   )                                ::   cdtype         ! =TRA or TRC (tracer indicator)
97      INTEGER         , INTENT(in   )                                ::   kjpt            ! number of tracers
98      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)                ::   p2dt        ! vertical profile of tracer time-step
99      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)   ::   ptrab          ! before and now tracer fields
100      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)   ::   ptraa          ! tracer trend
101      !!
102      INTEGER  ::  ji, jj, jk, jn        ! dummy loop indices
103      REAL(wp) ::  zavi, zrhs, znvvl     ! temporary scalars
104      REAL(wp) ::  ze3tb, ze3tn, ze3ta   ! variable vertical scale factors
105      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi, zwt   ! workspace arrays
106      !!---------------------------------------------------------------------
107
108      IF( kt == nit000 ) THEN
109         IF(lwp)WRITE(numout,*)
110         IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing'
111         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ '
112         zavi = 0.e0      ! avoid warning at compilation phase when lk_ldfslp=F
113      ENDIF
114      !
115      ! I. Local initialization
116      ! -----------------------
117      zwd(1,:, : ) = 0.e0     ;     zwd(jpi,:,:) = 0.e0
118      zws(1,:, : ) = 0.e0     ;     zws(jpi,:,:) = 0.e0
119      zwi(1,:, : ) = 0.e0     ;     zwi(jpi,:,:) = 0.e0
120      zwt(1,:, : ) = 0.e0     ;     zwt(jpi,:,:) = 0.e0
121      zwt(:,:,jpk) = 0.e0     ;     zwt( : ,:,1) = 0.e0
122
123      ! I.1 Variable volume : to take into account vertical variable vertical scale factors
124      ! -------------------
125      IF( lk_vvl ) THEN   ;    znvvl = 1.
126      ELSE                ;    znvvl = 0.e0
127      ENDIF
128
129      ! II. Vertical trend associated with the vertical physics
130      ! =======================================================
131      !     (including the vertical flux proportional to dk[t] associated
132      !      with the lateral mixing, through the avt update)
133      !     dk[ avt dk[ (t,s) ] ] diffusive trends
134
135      !
136      ! II.0 Matrix construction
137      ! ------------------------
138      DO jn = 1, kjpt
139         !
140         !  Matrix construction
141         ! ------------------------
142         IF( cdtype == 'TRA' .AND. jn == jp_tem )  THEN 
143#if defined key_ldfslp
144            ! update and save of avt (and avs if double diffusive mixing)
145            IF( l_traldf_rot ) THEN
146               DO jk = 2, jpkm1
147                  DO jj = 2, jpjm1
148                     DO ji = fs_2, fs_jpim1   ! vector opt.
149                        zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
150                           & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
151                           &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
152                        zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
153                     END DO
154                  END DO
155               END DO
156            ELSE                         ! no rotation but key_ldfslp defined
157               zwt  (:,:,:) = avt(:,:,:)
158            ENDIF
159#else
160            ! No isopycnal diffusion
161            zwt(:,:,:) = avt(:,:,:)           
162#endif
163            ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
164            DO jk = 1, jpkm1
165               DO jj = 2, jpjm1
166                  DO ji = fs_2, fs_jpim1   ! vector opt.
167                     ze3ta =  ( 1. - znvvl ) +        znvvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point
168                     ze3tn =         znvvl   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point
169                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
170                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
171                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
172                 END DO
173               END DO
174            END DO
175            ! Surface boudary conditions
176            DO jj = 2, jpjm1
177               DO ji = fs_2, fs_jpim1   ! vector opt.
178                 ze3ta = ( 1. - znvvl ) +  znvvl * fse3t_a(ji,jj,1)    ! after scale factor at T-point
179                 zwi(ji,jj,1) = 0.e0
180                 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
181               END DO
182            END DO
183            !
184         ELSE IF( ( cdtype == 'TRA' .AND. jn == jp_sal ) .OR. ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN
185#if defined key_ldfslp
186            ! update and save of avt (and avs if double diffusive mixing)
187            IF( l_traldf_rot ) THEN
188               DO jk = 2, jpkm1
189                  DO jj = 2, jpjm1
190                     DO ji = fs_2, fs_jpim1   ! vector opt.
191                        zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
192                           & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
193                           &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
194                        zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
195                     END DO
196                  END DO
197               END DO
198            ELSE                         ! no rotation but key_ldfslp defined
199               zwt  (:,:,:) = fsavs(:,:,:)
200            ENDIF
201#else
202            ! No isopycnal diffusion
203            zwt(:,:,:) = fsavs(:,:,:)           
204#endif
205            ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
206            DO jk = 1, jpkm1
207               DO jj = 2, jpjm1
208                  DO ji = fs_2, fs_jpim1   ! vector opt.
209                     ze3ta =  ( 1. - znvvl ) +        znvvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point
210                     ze3tn =         znvvl   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point
211                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
212                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
213                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
214                 END DO
215               END DO
216            END DO
217            ! Surface boudary conditions
218            DO jj = 2, jpjm1
219               DO ji = fs_2, fs_jpim1   ! vector opt.
220                 ze3ta = ( 1. - znvvl ) +  znvvl * fse3t_a(ji,jj,1)    ! after scale factor at T-point
221                 zwi(ji,jj,1) = 0.e0
222                 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
223               END DO
224            END DO
225            !
226         END IF
227
228         ! II.1. Vertical diffusion on tracer
229         ! ---------------------------
230
231         !! Matrix inversion from the first level
232         !!----------------------------------------------------------------------
233         !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
234         !
235         !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
236         !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
237         !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
238         !        (        ...               )( ...  ) ( ...  )
239         !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
240         !
241         !   m is decomposed in the product of an upper and lower triangular matrix
242         !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
243         !   The second member is in 2d array zwy
244         !   The solution is in 2d array zwx
245         !   The 3d arry zwt is a work space array
246         !   zwy is used and then used as a work space array : its value is modified!
247         
248         ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
249         DO jj = 2, jpjm1
250            DO ji = fs_2, fs_jpim1
251               zwt(ji,jj,1) = zwd(ji,jj,1)
252            END DO
253         END DO
254         DO jk = 2, jpkm1
255            DO jj = 2, jpjm1
256               DO ji = fs_2, fs_jpim1
257                  zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
258               END DO
259            END DO
260         END DO
261         
262         ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
263         DO jj = 2, jpjm1
264            DO ji = fs_2, fs_jpim1
265               ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,1)
266               ze3tn = ( 1. - znvvl ) + znvvl * fse3t(ji,jj,1)
267               ptraa(ji,jj,1,jn) = ze3tb * ptrab(ji,jj,1,jn) + p2dt(1) * ze3tn * ptraa(ji,jj,1,jn)
268            END DO
269         END DO
270         DO jk = 2, jpkm1
271            DO jj = 2, jpjm1
272               DO ji = fs_2, fs_jpim1
273                  ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,jk)
274                  ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,jk)
275                  zrhs = ze3tb * ptrab(ji,jj,jk,jn) + p2dt(jk) * ze3tn * ptraa(ji,jj,jk,jn)   ! zrhs=right hand side
276                  ptraa(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * ptraa(ji,jj,jk-1,jn)
277               END DO
278            END DO
279         END DO
280         
281         ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
282         ! Save the masked temperature after in ta
283         ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
284         DO jj = 2, jpjm1
285            DO ji = fs_2, fs_jpim1
286               ptraa(ji,jj,jpkm1,jn) = ptraa(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
287            END DO
288         END DO
289         DO jk = jpk-2, 1, -1
290            DO jj = 2, jpjm1
291               DO ji = fs_2, fs_jpim1
292                  ptraa(ji,jj,jk,jn) = ( ptraa(ji,jj,jk,jn) - zws(ji,jj,jk) * ptraa(ji,jj,jk+1,jn) ) &
293                  &                    / zwt(ji,jj,jk) * tmask(ji,jj,jk)
294               END DO
295            END DO
296         END DO
297         !
298      END DO
299      !
300   END SUBROUTINE tra_zdf_imp
301
302   !!==============================================================================
303END MODULE trazdf_imp
Note: See TracBrowser for help on using the repository browser.