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

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

cosmetic changes
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 14.7 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.3 , LOCEAN-IPSL (2010)
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, ptb, pta, kjpt ) 
52      !!----------------------------------------------------------------------
53      !!                  ***  ROUTINE tra_zdf_imp  ***
54      !!
55      !! ** Purpose :   Compute the trend due to the vertical tracer diffusion
56      !!     including the vertical component of lateral mixing (only for 2nd
57      !!     order operator, for fourth order it is already computed and add
58      !!     to the general trend in traldf.F) and add it to the general trend
59      !!     of the tracer equations.
60      !!
61      !! ** Method  :   The vertical component of the lateral diffusive trends
62      !!      is provided by a 2nd order operator rotated along neutral or geo-
63      !!      potential surfaces to which an eddy induced advection can be
64      !!      added. It is computed using before fields (forward in time) and
65      !!      isopycnal or geopotential slopes computed in routine ldfslp.
66      !!
67      !!      Second part: vertical trend associated with the vertical physics
68      !!      ===========  (including the vertical flux proportional to dk[t]
69      !!                  associated with the lateral mixing, through the
70      !!                  update of avt)
71      !!      The vertical diffusion of the tracer t  is given by:
72      !!             difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
73      !!      It is computed using a backward time scheme (t=ta).
74      !!      Surface and bottom boundary conditions: no diffusive flux on
75      !!      both tracers (bottom, applied through the masked field avt).
76      !!      Add this trend to the general trend ta,sa :
77      !!         ta = ta + dz( avt dz(t) )
78      !!         if lk_zdfddm=T, use avs for salinity or for passive tracers
79      !!         (sa = sa + dz( avs dz(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      !!         if lk_zdfddm=T, use avt = zavt
85      !!         (avs = zavs if lk_zdfddm=T )
86      !!
87      !! ** Action  : - Update (ta) with before vertical diffusion trend
88      !!
89      !!---------------------------------------------------------------------
90      !!
91      USE oce    , ONLY :   zwd   => ua   ! ua used as workspace
92      USE oce    , ONLY :   zws   => va   ! va  -          -
93      !!
94      INTEGER         , INTENT(in   )                                ::   kt             ! ocean time-step index
95      CHARACTER(len=3), INTENT(in   )                                ::   cdtype         ! =TRA or TRC (tracer indicator)
96      INTEGER         , INTENT(in   )                                ::   kjpt            ! number of tracers
97      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)                ::   p2dt        ! vertical profile of tracer time-step
98      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)   ::   ptb          ! before and now tracer fields
99      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)   ::   pta          ! tracer trend
100      !!
101      INTEGER  ::  ji, jj, jk, jn        ! dummy loop indices
102      REAL(wp) ::  zavi, zrhs, znvvl     ! temporary scalars
103      REAL(wp) ::  ze3tb, ze3tn, ze3ta   ! variable vertical scale factors
104      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi, zwt   ! workspace arrays
105      !!---------------------------------------------------------------------
106
107      IF( kt == nit000 ) THEN
108         IF(lwp)WRITE(numout,*)
109         IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing'
110         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ '
111         zavi = 0.e0      ! avoid warning at compilation phase when lk_ldfslp=F
112      ENDIF
113      !
114      ! I. Local initialization
115      ! -----------------------
116      zwd(1,:, : ) = 0.e0     ;     zwd(jpi,:,:) = 0.e0
117      zws(1,:, : ) = 0.e0     ;     zws(jpi,:,:) = 0.e0
118      zwi(1,:, : ) = 0.e0     ;     zwi(jpi,:,:) = 0.e0
119      zwt(1,:, : ) = 0.e0     ;     zwt(jpi,:,:) = 0.e0
120      zwt(:,:,jpk) = 0.e0     ;     zwt( : ,:,1) = 0.e0
121
122      ! I.1 Variable volume : to take into account vertical variable vertical scale factors
123      ! -------------------
124      IF( lk_vvl ) THEN   ;    znvvl = 1.
125      ELSE                ;    znvvl = 0.e0
126      ENDIF
127
128      ! II. Vertical trend associated with the vertical physics
129      ! =======================================================
130      !     (including the vertical flux proportional to dk[t] associated
131      !      with the lateral mixing, through the avt update)
132      !     dk[ avt dk[ (t,s) ] ] diffusive trends
133
134      !
135      ! II.0 Matrix construction
136      ! ------------------------
137      DO jn = 1, kjpt
138         !
139         !  Matrix construction
140         ! ------------------------
141         IF( cdtype == 'TRA' .AND. jn == jp_tem )  THEN 
142#if defined key_ldfslp
143            ! update and save of avt (and avs if double diffusive mixing)
144            IF( l_traldf_rot ) THEN
145               DO jk = 2, jpkm1
146                  DO jj = 2, jpjm1
147                     DO ji = fs_2, fs_jpim1   ! vector opt.
148                        zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
149                           & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
150                           &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
151                        zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
152                     END DO
153                  END DO
154               END DO
155            ELSE                         ! no rotation but key_ldfslp defined
156               zwt  (:,:,:) = avt(:,:,:)
157            ENDIF
158#else
159            ! No isopycnal diffusion
160            zwt(:,:,:) = avt(:,:,:)           
161#endif
162            ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
163            DO jk = 1, jpkm1
164               DO jj = 2, jpjm1
165                  DO ji = fs_2, fs_jpim1   ! vector opt.
166                     ze3ta =  ( 1. - znvvl ) +        znvvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point
167                     ze3tn =         znvvl   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point
168                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
169                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
170                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
171                 END DO
172               END DO
173            END DO
174            ! Surface boudary conditions
175            DO jj = 2, jpjm1
176               DO ji = fs_2, fs_jpim1   ! vector opt.
177                 ze3ta = ( 1. - znvvl ) +  znvvl * fse3t_a(ji,jj,1)    ! after scale factor at T-point
178                 zwi(ji,jj,1) = 0.e0
179                 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
180               END DO
181            END DO
182            !
183         ELSE IF( ( cdtype == 'TRA' .AND. jn == jp_sal ) .OR. ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN
184#if defined key_ldfslp
185            ! update and save of avt (and avs if double diffusive mixing)
186            IF( l_traldf_rot ) THEN
187               DO jk = 2, jpkm1
188                  DO jj = 2, jpjm1
189                     DO ji = fs_2, fs_jpim1   ! vector opt.
190                        zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
191                           & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
192                           &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
193                        zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
194                     END DO
195                  END DO
196               END DO
197            ELSE                         ! no rotation but key_ldfslp defined
198               zwt  (:,:,:) = fsavs(:,:,:)
199            ENDIF
200#else
201            ! No isopycnal diffusion
202            zwt(:,:,:) = fsavs(:,:,:)           
203#endif
204            ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
205            DO jk = 1, jpkm1
206               DO jj = 2, jpjm1
207                  DO ji = fs_2, fs_jpim1   ! vector opt.
208                     ze3ta =  ( 1. - znvvl ) +        znvvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point
209                     ze3tn =         znvvl   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point
210                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
211                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
212                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
213                 END DO
214               END DO
215            END DO
216            ! Surface boudary conditions
217            DO jj = 2, jpjm1
218               DO ji = fs_2, fs_jpim1   ! vector opt.
219                 ze3ta = ( 1. - znvvl ) +  znvvl * fse3t_a(ji,jj,1)    ! 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         END IF
226
227         ! II.1. Vertical diffusion on tracer
228         ! ---------------------------
229
230         !! Matrix inversion from the first level
231         !!----------------------------------------------------------------------
232         !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
233         !
234         !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
235         !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
236         !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
237         !        (        ...               )( ...  ) ( ...  )
238         !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
239         !
240         !   m is decomposed in the product of an upper and lower triangular matrix
241         !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
242         !   The second member is in 2d array zwy
243         !   The solution is in 2d array zwx
244         !   The 3d arry zwt is a work space array
245         !   zwy is used and then used as a work space array : its value is modified!
246         
247         ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
248         DO jj = 2, jpjm1
249            DO ji = fs_2, fs_jpim1
250               zwt(ji,jj,1) = zwd(ji,jj,1)
251            END DO
252         END DO
253         DO jk = 2, jpkm1
254            DO jj = 2, jpjm1
255               DO ji = fs_2, fs_jpim1
256                  zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
257               END DO
258            END DO
259         END DO
260         
261         ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
262         DO jj = 2, jpjm1
263            DO ji = fs_2, fs_jpim1
264               ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,1)
265               ze3tn = ( 1. - znvvl ) + znvvl * fse3t(ji,jj,1)
266               pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) + p2dt(1) * ze3tn * pta(ji,jj,1,jn)
267            END DO
268         END DO
269         DO jk = 2, jpkm1
270            DO jj = 2, jpjm1
271               DO ji = fs_2, fs_jpim1
272                  ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,jk)
273                  ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,jk)
274                  zrhs = ze3tb * ptb(ji,jj,jk,jn) + p2dt(jk) * ze3tn * pta(ji,jj,jk,jn)   ! zrhs=right hand side
275                  pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn)
276               END DO
277            END DO
278         END DO
279         
280         ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
281         ! Save the masked temperature after in ta
282         ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
283         DO jj = 2, jpjm1
284            DO ji = fs_2, fs_jpim1
285               pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
286            END DO
287         END DO
288         DO jk = jpk-2, 1, -1
289            DO jj = 2, jpjm1
290               DO ji = fs_2, fs_jpim1
291                  pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) ) &
292                  &                    / zwt(ji,jj,jk) * tmask(ji,jj,jk)
293               END DO
294            END DO
295         END DO
296         !
297      END DO
298      !
299   END SUBROUTINE tra_zdf_imp
300
301   !!==============================================================================
302END MODULE trazdf_imp
Note: See TracBrowser for help on using the repository browser.