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

Last change on this file since 2371 was 2371, checked in by acc, 13 years ago

nemo_v3_3_beta. Improvement of the Griffies operator code (#680). Some aspects are still undergoing scientific assessment but the code structure is ready for release. Dynamical allocation of the large triad arrays has been introduced to reduce memory when the new operator is not in use.

  • 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 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 , NEMO Consortium (2010)
47   !! $Id$
48   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_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      USE traldf_iso_grif    , ONLY :   ah_wslp2
94      !!
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), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt     ! vertical profile of tracer time-step
99      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb      ! before and now tracer fields
100      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend
101      !!
102      INTEGER  ::  ji, jj, jk, jn        ! dummy loop indices
103      REAL(wp) ::  zavi, zrhs, znvvl     ! local 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 on ', cdtype
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            IF( ln_traldf_grif ) 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) * wslp2(ji,jj,jk)         ! vertical mixing coef. due to lateral mixing
149                        zavi = ah_wslp2(ji,jj,jk)                ! vertical mixing coef. due to lateral mixing
150                        zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
151                     END DO
152                  END DO
153               END DO
154            ! update and save of avt (and avs if double diffusive mixing)
155            ELSE IF( l_traldf_rot ) THEN
156               DO jk = 2, jpkm1
157                  DO jj = 2, jpjm1
158                     DO ji = fs_2, fs_jpim1   ! vector opt.
159                        zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
160                           & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
161                           &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
162                        zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
163                     END DO
164                  END DO
165               END DO
166            ELSE                         ! no rotation but key_ldfslp defined
167               zwt  (:,:,:) = avt(:,:,:)
168            ENDIF
169#else
170            ! No isopycnal diffusion
171            zwt(:,:,:) = avt(:,:,:)           
172#endif
173            ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
174            DO jk = 1, jpkm1
175               DO jj = 2, jpjm1
176                  DO ji = fs_2, fs_jpim1   ! vector opt.
177                     ze3ta =  ( 1. - znvvl ) +        znvvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point
178                     ze3tn =         znvvl   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point
179                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
180                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
181                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
182                 END DO
183               END DO
184            END DO
185            ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
186            DO jj = 2, jpjm1
187               DO ji = fs_2, fs_jpim1
188                  zwt(ji,jj,1) = zwd(ji,jj,1)
189               END DO
190            END DO
191            DO jk = 2, jpkm1
192               DO jj = 2, jpjm1
193                  DO ji = fs_2, fs_jpim1
194                    zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
195                  END DO
196               END DO
197            END DO
198            !
199         ELSE IF( ( cdtype == 'TRA' .AND. jn == jp_sal ) .OR. ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN
200#if defined key_ldfslp
201            IF( ln_traldf_grif ) THEN
202               DO jk = 2, jpkm1
203                  DO jj = 2, jpjm1
204                     DO ji = fs_2, fs_jpim1   ! vector opt.
205                        zavi = ah_wslp2(ji,jj,jk)                ! vertical mixing coef. due to lateral mixing
206                        ! zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk)         ! vertical mixing coef. due to lateral mixing
207                        zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
208                     END DO
209                  END DO
210               END DO
211            ELSE IF( l_traldf_rot ) THEN
212               DO jk = 2, jpkm1
213                  DO jj = 2, jpjm1
214                     DO ji = fs_2, fs_jpim1   ! vector opt.
215                        zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
216                           & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
217                           &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
218                        zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on salinity)
219                     END DO
220                  END DO
221               END DO
222            ELSE                         ! no rotation but key_ldfslp defined
223               zwt  (:,:,:) = fsavs(:,:,:)
224            ENDIF
225#else
226            ! No isopycnal diffusion
227            zwt(:,:,:) = fsavs(:,:,:)           
228#endif
229            ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
230            DO jk = 1, jpkm1
231               DO jj = 2, jpjm1
232                  DO ji = fs_2, fs_jpim1   ! vector opt.
233                     ze3ta =  ( 1. - znvvl ) +        znvvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point
234                     ze3tn =         znvvl   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point
235                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
236                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
237                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
238                 END DO
239               END DO
240            END DO
241            ! Surface boudary conditions
242            DO jj = 2, jpjm1
243               DO ji = fs_2, fs_jpim1   ! vector opt.
244                 ze3ta = ( 1. - znvvl ) +  znvvl * fse3t_a(ji,jj,1)    ! after scale factor at T-point
245                 zwi(ji,jj,1) = 0.e0
246                 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
247               END DO
248            END DO
249            !
250            ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
251            DO jj = 2, jpjm1
252               DO ji = fs_2, fs_jpim1
253                  zwt(ji,jj,1) = zwd(ji,jj,1)
254               END DO
255            END DO
256            DO jk = 2, jpkm1
257               DO jj = 2, jpjm1
258                  DO ji = fs_2, fs_jpim1
259                    zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
260                  END DO
261               END DO
262            END DO
263            !
264         END IF 
265
266         ! II.1. Vertical diffusion on tracer
267         ! ---------------------------
268         
269         ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
270         DO jj = 2, jpjm1
271            DO ji = fs_2, fs_jpim1
272               ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,1)
273               ze3tn = ( 1. - znvvl ) + znvvl * fse3t(ji,jj,1)
274               pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) + p2dt(1) * ze3tn * pta(ji,jj,1,jn)
275            END DO
276         END DO
277         DO jk = 2, jpkm1
278            DO jj = 2, jpjm1
279               DO ji = fs_2, fs_jpim1
280                  ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,jk)
281                  ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,jk)
282                  zrhs = ze3tb * ptb(ji,jj,jk,jn) + p2dt(jk) * ze3tn * pta(ji,jj,jk,jn)   ! zrhs=right hand side
283                  pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn)
284               END DO
285            END DO
286         END DO
287
288         ! third recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk
289         DO jj = 2, jpjm1
290            DO ji = fs_2, fs_jpim1
291               pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
292            END DO
293         END DO
294         DO jk = jpk-2, 1, -1
295            DO jj = 2, jpjm1
296               DO ji = fs_2, fs_jpim1
297                  pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) ) &
298                     &             / zwt(ji,jj,jk) * tmask(ji,jj,jk)
299               END DO
300            END DO
301         END DO
302         !
303      END DO
304      !
305   END SUBROUTINE tra_zdf_imp
306
307   !!==============================================================================
308END MODULE trazdf_imp
Note: See TracBrowser for help on using the repository browser.