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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90 @ 2281

Last change on this file since 2281 was 2281, checked in by smasson, 14 years ago

set proper svn properties to all files...
  • 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   USE trc_oce         ! share passive tracers/Ocean variables
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   tra_zdf_imp   !  routine called by step.F90
39
40   !! * Substitutions
41#  include "domzgr_substitute.h90"
42#  include "ldftra_substitute.h90"
43#  include "zdfddm_substitute.h90"
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)
47   !! $Id$
48   !! Software governed by the CeCILL licence  (NEMOGCM/License_CeCILL.txt)
49   !!----------------------------------------------------------------------
50CONTAINS
51 
52   SUBROUTINE tra_zdf_imp( kt, cdtype, p2dt, ptb, pta, 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      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), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt     ! vertical profile of tracer time-step
98      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb      ! before and now tracer fields
99      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend
100      !!
101      INTEGER  ::  ji, jj, jk, jn        ! dummy loop indices
102      REAL(wp) ::  zavi, zrhs, znvvl     ! local 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 on ', cdtype
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            IF( ln_traldf_grif ) THEN
144               DO jk = 2, jpkm1
145                  DO jj = 2, jpjm1
146                     DO ji = fs_2, fs_jpim1   ! vector opt.
147                        zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk)         ! vertical mixing coef. due to lateral mixing
148                        zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
149                     END DO
150                  END DO
151               END DO
152            ! update and save of avt (and avs if double diffusive mixing)
153            ELSE IF( l_traldf_rot ) THEN
154               DO jk = 2, jpkm1
155                  DO jj = 2, jpjm1
156                     DO ji = fs_2, fs_jpim1   ! vector opt.
157                        zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
158                           & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
159                           &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
160                        zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
161                     END DO
162                  END DO
163               END DO
164            ELSE                         ! no rotation but key_ldfslp defined
165               zwt  (:,:,:) = avt(:,:,:)
166            ENDIF
167#else
168            ! No isopycnal diffusion
169            zwt(:,:,:) = avt(:,:,:)           
170#endif
171            ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
172            DO jk = 1, jpkm1
173               DO jj = 2, jpjm1
174                  DO ji = fs_2, fs_jpim1   ! vector opt.
175                     ze3ta =  ( 1. - znvvl ) +        znvvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point
176                     ze3tn =         znvvl   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point
177                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
178                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
179                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
180                 END DO
181               END DO
182            END DO
183            ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
184            DO jj = 2, jpjm1
185               DO ji = fs_2, fs_jpim1
186                  zwt(ji,jj,1) = zwd(ji,jj,1)
187               END DO
188            END DO
189            DO jk = 2, jpkm1
190               DO jj = 2, jpjm1
191                  DO ji = fs_2, fs_jpim1
192                    zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
193                  END DO
194               END DO
195            END DO
196            !
197         ELSE IF( ( cdtype == 'TRA' .AND. jn == jp_sal ) .OR. ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN
198#if defined key_ldfslp
199            IF( ln_traldf_grif ) 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) * wslp2(ji,jj,jk)         ! vertical mixing coef. due to lateral mixing
204                        zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
205                     END DO
206                  END DO
207               END DO
208            ELSE IF( l_traldf_rot ) THEN
209               DO jk = 2, jpkm1
210                  DO jj = 2, jpjm1
211                     DO ji = fs_2, fs_jpim1   ! vector opt.
212                        zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
213                           & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
214                           &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
215                        zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on salinity)
216                     END DO
217                  END DO
218               END DO
219            ELSE                         ! no rotation but key_ldfslp defined
220               zwt  (:,:,:) = fsavs(:,:,:)
221            ENDIF
222#else
223            ! No isopycnal diffusion
224            zwt(:,:,:) = fsavs(:,:,:)           
225#endif
226            ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
227            DO jk = 1, jpkm1
228               DO jj = 2, jpjm1
229                  DO ji = fs_2, fs_jpim1   ! vector opt.
230                     ze3ta =  ( 1. - znvvl ) +        znvvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point
231                     ze3tn =         znvvl   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point
232                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
233                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
234                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
235                 END DO
236               END DO
237            END DO
238            ! Surface boudary conditions
239            DO jj = 2, jpjm1
240               DO ji = fs_2, fs_jpim1   ! vector opt.
241                 ze3ta = ( 1. - znvvl ) +  znvvl * fse3t_a(ji,jj,1)    ! after scale factor at T-point
242                 zwi(ji,jj,1) = 0.e0
243                 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
244               END DO
245            END DO
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         END IF 
262
263         ! II.1. Vertical diffusion on tracer
264         ! ---------------------------
265         
266         ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
267         DO jj = 2, jpjm1
268            DO ji = fs_2, fs_jpim1
269               ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,1)
270               ze3tn = ( 1. - znvvl ) + znvvl * fse3t(ji,jj,1)
271               pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) + p2dt(1) * ze3tn * pta(ji,jj,1,jn)
272            END DO
273         END DO
274         DO jk = 2, jpkm1
275            DO jj = 2, jpjm1
276               DO ji = fs_2, fs_jpim1
277                  ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,jk)
278                  ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,jk)
279                  zrhs = ze3tb * ptb(ji,jj,jk,jn) + p2dt(jk) * ze3tn * pta(ji,jj,jk,jn)   ! zrhs=right hand side
280                  pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn)
281               END DO
282            END DO
283         END DO
284
285         ! third recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk
286         DO jj = 2, jpjm1
287            DO ji = fs_2, fs_jpim1
288               pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
289            END DO
290         END DO
291         DO jk = jpk-2, 1, -1
292            DO jj = 2, jpjm1
293               DO ji = fs_2, fs_jpim1
294                  pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) ) &
295                     &             / zwt(ji,jj,jk) * tmask(ji,jj,jk)
296               END DO
297            END DO
298         END DO
299         !
300      END DO
301      !
302   END SUBROUTINE tra_zdf_imp
303
304   !!==============================================================================
305END MODULE trazdf_imp
Note: See TracBrowser for help on using the repository browser.