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

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90 @ 2528

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago

Update NEMOGCM from branch nemo_v3_3_beta
  • Property svn:keywords set to Id
File size: 15.0 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 trc_oce         ! share passive tracers/ocean variables
27   USE domvvl          ! variable volume
28   USE ldftra_oce      ! ocean active tracers: lateral physics
29   USE ldftra          ! lateral mixing type
30   USE ldfslp          ! lateral physics: slope of diffusion
31   USE zdfddm          ! ocean vertical physics: double diffusion
32   USE traldf_iso_grif ! active tracers: Griffies operator
33   USE in_out_manager  ! I/O manager
34   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
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   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
48   !! $Id$
49   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
50   !!----------------------------------------------------------------------
51CONTAINS
52 
53   SUBROUTINE tra_zdf_imp( kt, cdtype, p2dt, ptb, pta, kjpt ) 
54      !!----------------------------------------------------------------------
55      !!                  ***  ROUTINE tra_zdf_imp  ***
56      !!
57      !! ** Purpose :   Compute the trend due to the vertical tracer diffusion
58      !!     including the vertical component of lateral mixing (only for 2nd
59      !!     order operator, for fourth order it is already computed and add
60      !!     to the general trend in traldf.F) and add it to the general trend
61      !!     of the tracer equations.
62      !!
63      !! ** Method  :   The vertical component of the lateral diffusive trends
64      !!      is provided by a 2nd order operator rotated along neutral or geo-
65      !!      potential surfaces to which an eddy induced advection can be
66      !!      added. It is computed using before fields (forward in time) and
67      !!      isopycnal or geopotential slopes computed in routine ldfslp.
68      !!
69      !!      Second part: vertical trend associated with the vertical physics
70      !!      ===========  (including the vertical flux proportional to dk[t]
71      !!                  associated with the lateral mixing, through the
72      !!                  update of avt)
73      !!      The vertical diffusion of the tracer t  is given by:
74      !!             difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
75      !!      It is computed using a backward time scheme (t=ta).
76      !!      Surface and bottom boundary conditions: no diffusive flux on
77      !!      both tracers (bottom, applied through the masked field avt).
78      !!      Add this trend to the general trend ta,sa :
79      !!         ta = ta + dz( avt dz(t) )
80      !!         if lk_zdfddm=T, use avs for salinity or for passive tracers
81      !!         (sa = sa + dz( avs dz(t) )
82      !!
83      !!      Third part: recover avt resulting from the vertical physics
84      !!      ==========  alone, for further diagnostics (for example to
85      !!                  compute the turbocline depth in zdfmxl.F90).
86      !!         if lk_zdfddm=T, use avt = zavt
87      !!         (avs = zavs if lk_zdfddm=T )
88      !!
89      !! ** Action  : - Update (ta) with before vertical diffusion trend
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._wp      ! avoid warning at compilation phase when lk_ldfslp=F
112      ENDIF
113      !
114      ! I. Local initialization
115      ! -----------------------
116      zwd(1,:, : ) = 0._wp     ;     zwd(jpi,:,:) = 0._wp
117      zws(1,:, : ) = 0._wp     ;     zws(jpi,:,:) = 0._wp
118      zwi(1,:, : ) = 0._wp     ;     zwi(jpi,:,:) = 0._wp
119      zwt(1,:, : ) = 0._wp     ;     zwt(jpi,:,:) = 0._wp
120      zwt(:,:,jpk) = 0._wp     ;     zwt( : ,:,1) = 0._wp
121
122      ! I.1 Variable volume : to take into account vertical variable vertical scale factors
123      ! -------------------
124      IF( lk_vvl ) THEN   ;    znvvl = 1._wp
125      ELSE                ;    znvvl = 0._wp
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                        zavi = ah_wslp2(ji,jj,jk)                ! vertical mixing coef. due to lateral mixing
149                        zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
150                     END DO
151                  END DO
152               END DO
153            ! update and save of avt (and avs if double diffusive mixing)
154            ELSE IF( l_traldf_rot ) THEN
155               DO jk = 2, jpkm1
156                  DO jj = 2, jpjm1
157                     DO ji = fs_2, fs_jpim1   ! vector opt.
158                        zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
159                           & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
160                           &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
161                        zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
162                     END DO
163                  END DO
164               END DO
165            ELSE                         ! no rotation but key_ldfslp defined
166               zwt(:,:,:) = avt(:,:,:)
167            ENDIF
168#else
169            ! No isopycnal diffusion
170            zwt(:,:,:) = avt(:,:,:)           
171#endif
172            ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
173            DO jk = 1, jpkm1
174               DO jj = 2, jpjm1
175                  DO ji = fs_2, fs_jpim1   ! vector opt.
176                     ze3ta =  ( 1. - znvvl ) +        znvvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point
177                     ze3tn =         znvvl   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point
178                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
179                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
180                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
181                 END DO
182               END DO
183            END DO
184            ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
185            DO jj = 2, jpjm1
186               DO ji = fs_2, fs_jpim1
187                  zwt(ji,jj,1) = zwd(ji,jj,1)
188               END DO
189            END DO
190            DO jk = 2, jpkm1
191               DO jj = 2, jpjm1
192                  DO ji = fs_2, fs_jpim1
193                    zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
194                  END DO
195               END DO
196            END DO
197            !
198         ELSE IF( ( cdtype == 'TRA' .AND. jn == jp_sal ) .OR. ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN
199#if defined key_ldfslp
200            IF( ln_traldf_grif ) THEN
201               DO jk = 2, jpkm1
202                  DO jj = 2, jpjm1
203                     DO ji = fs_2, fs_jpim1   ! vector opt.
204                        zavi = ah_wslp2(ji,jj,jk)                ! vertical mixing coef. due to lateral mixing
205                        ! zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk)         ! vertical mixing coef. due to lateral mixing
206                        zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
207                     END DO
208                  END DO
209               END DO
210            ELSE IF( l_traldf_rot ) THEN
211               DO jk = 2, jpkm1
212                  DO jj = 2, jpjm1
213                     DO ji = fs_2, fs_jpim1   ! vector opt.
214                        zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
215                           & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
216                           &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
217                        zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on salinity)
218                     END DO
219                  END DO
220               END DO
221            ELSE                         ! no rotation but key_ldfslp defined
222               zwt(:,:,:) = fsavs(:,:,:)
223            ENDIF
224#else
225            ! No isopycnal diffusion
226            zwt(:,:,:) = fsavs(:,:,:)           
227#endif
228            ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
229            DO jk = 1, jpkm1
230               DO jj = 2, jpjm1
231                  DO ji = fs_2, fs_jpim1   ! vector opt.
232                     ze3ta =  ( 1. - znvvl ) +        znvvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point
233                     ze3tn =         znvvl   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point
234                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
235                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
236                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
237                 END DO
238               END DO
239            END DO
240            ! Surface boudary conditions
241            DO jj = 2, jpjm1
242               DO ji = fs_2, fs_jpim1   ! vector opt.
243                 ze3ta = ( 1. - znvvl ) +  znvvl * fse3t_a(ji,jj,1)    ! after scale factor at T-point
244                 zwi(ji,jj,1) = 0._wp
245                 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
246               END DO
247            END DO
248            !
249            ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
250            DO jj = 2, jpjm1
251               DO ji = fs_2, fs_jpim1
252                  zwt(ji,jj,1) = zwd(ji,jj,1)
253               END DO
254            END DO
255            DO jk = 2, jpkm1
256               DO jj = 2, jpjm1
257                  DO ji = fs_2, fs_jpim1
258                    zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
259                  END DO
260               END DO
261            END DO
262            !
263         END IF 
264
265         ! II.1. Vertical diffusion on tracer
266         ! ---------------------------
267         
268         ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
269         DO jj = 2, jpjm1
270            DO ji = fs_2, fs_jpim1
271               ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,1)
272               ze3tn = ( 1. - znvvl ) + znvvl * fse3t(ji,jj,1)
273               pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) + p2dt(1) * ze3tn * pta(ji,jj,1,jn)
274            END DO
275         END DO
276         DO jk = 2, jpkm1
277            DO jj = 2, jpjm1
278               DO ji = fs_2, fs_jpim1
279                  ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,jk)
280                  ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,jk)
281                  zrhs = ze3tb * ptb(ji,jj,jk,jn) + p2dt(jk) * ze3tn * pta(ji,jj,jk,jn)   ! zrhs=right hand side
282                  pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn)
283               END DO
284            END DO
285         END DO
286
287         ! third recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk
288         DO jj = 2, jpjm1
289            DO ji = fs_2, fs_jpim1
290               pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
291            END DO
292         END DO
293         DO jk = jpk-2, 1, -1
294            DO jj = 2, jpjm1
295               DO ji = fs_2, fs_jpim1
296                  pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) )   &
297                     &             / zwt(ji,jj,jk) * tmask(ji,jj,jk)
298               END DO
299            END DO
300         END DO
301         !
302      END DO
303      !
304   END SUBROUTINE tra_zdf_imp
305
306   !!==============================================================================
307END MODULE trazdf_imp
Note: See TracBrowser for help on using the repository browser.