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

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

Improve the merge TRA-TRC, see ticket:701
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 14.3 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            ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
183            DO jj = 2, jpjm1
184               DO ji = fs_2, fs_jpim1
185                  zwt(ji,jj,1) = zwd(ji,jj,1)
186               END DO
187            END DO
188            DO jk = 2, jpkm1
189               DO jj = 2, jpjm1
190                  DO ji = fs_2, fs_jpim1
191                    zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
192                  END DO
193               END DO
194            END DO
195            !
196         ELSE IF( ( cdtype == 'TRA' .AND. jn == jp_sal ) .OR. ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN
197#if defined key_ldfslp
198            ! update and save of avt (and avs if double diffusive mixing)
199            IF( l_traldf_rot ) THEN
200               DO jk = 2, jpkm1
201                  DO jj = 2, jpjm1
202                     DO ji = fs_2, fs_jpim1   ! vector opt.
203                        zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
204                           & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
205                           &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
206                        zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on salinity)
207                     END DO
208                  END DO
209               END DO
210            ELSE                         ! no rotation but key_ldfslp defined
211               zwt  (:,:,:) = fsavs(:,:,:)
212            ENDIF
213#else
214            ! No isopycnal diffusion
215            zwt(:,:,:) = fsavs(:,:,:)           
216#endif
217            ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
218            DO jk = 1, jpkm1
219               DO jj = 2, jpjm1
220                  DO ji = fs_2, fs_jpim1   ! vector opt.
221                     ze3ta =  ( 1. - znvvl ) +        znvvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point
222                     ze3tn =         znvvl   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point
223                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
224                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
225                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
226                 END DO
227               END DO
228            END DO
229            ! Surface boudary conditions
230            DO jj = 2, jpjm1
231               DO ji = fs_2, fs_jpim1   ! vector opt.
232                 ze3ta = ( 1. - znvvl ) +  znvvl * fse3t_a(ji,jj,1)    ! after scale factor at T-point
233                 zwi(ji,jj,1) = 0.e0
234                 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
235               END DO
236            END DO
237            !
238            ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
239            DO jj = 2, jpjm1
240               DO ji = fs_2, fs_jpim1
241                  zwt(ji,jj,1) = zwd(ji,jj,1)
242               END DO
243            END DO
244            DO jk = 2, jpkm1
245               DO jj = 2, jpjm1
246                  DO ji = fs_2, fs_jpim1
247                    zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
248                  END DO
249               END DO
250            END DO
251            !
252         END IF 
253
254         ! II.1. Vertical diffusion on tracer
255         ! ---------------------------
256         
257         ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
258         DO jj = 2, jpjm1
259            DO ji = fs_2, fs_jpim1
260               ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,1)
261               ze3tn = ( 1. - znvvl ) + znvvl * fse3t(ji,jj,1)
262               pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) + p2dt(1) * ze3tn * pta(ji,jj,1,jn)
263            END DO
264         END DO
265         DO jk = 2, jpkm1
266            DO jj = 2, jpjm1
267               DO ji = fs_2, fs_jpim1
268                  ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,jk)
269                  ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,jk)
270                  zrhs = ze3tb * ptb(ji,jj,jk,jn) + p2dt(jk) * ze3tn * pta(ji,jj,jk,jn)   ! zrhs=right hand side
271                  pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn)
272               END DO
273            END DO
274         END DO
275         
276         ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
277         ! Save the masked temperature after in ta
278         ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
279         DO jj = 2, jpjm1
280            DO ji = fs_2, fs_jpim1
281               pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
282            END DO
283         END DO
284         DO jk = jpk-2, 1, -1
285            DO jj = 2, jpjm1
286               DO ji = fs_2, fs_jpim1
287                  pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) ) &
288                  &                    / zwt(ji,jj,jk) * tmask(ji,jj,jk)
289               END DO
290            END DO
291         END DO
292         !
293      END DO
294      !
295   END SUBROUTINE tra_zdf_imp
296
297   !!==============================================================================
298END MODULE trazdf_imp
Note: See TracBrowser for help on using the repository browser.