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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90 @ 2633

Last change on this file since 2633 was 2633, checked in by trackstand2, 13 years ago

Renamed wrk_use => wrk_in_use and wrk_release => wrk_not_released
  • Property svn:keywords set to Id
File size: 15.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 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      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
94      USE wrk_nemo, ONLY: zwi => wrk_3d_1, zwt => wrk_3d_2  ! workspace arrays
95      !!
96      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
97      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
98      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
99      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt     ! vertical profile of tracer time-step
100      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb      ! before and now tracer fields
101      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend
102      !!
103      INTEGER  ::  ji, jj, jk, jn        ! dummy loop indices
104      REAL(wp) ::  zavi, zrhs, znvvl     ! local scalars
105      REAL(wp) ::  ze3tb, ze3tn, ze3ta   ! variable vertical scale factors
106      !!---------------------------------------------------------------------
107
108      IF(wrk_in_use(3, 1,2))THEN
109         CALL ctl_stop('tra_zdf_imp : requested workspace arrays unavailable.')
110         RETURN
111      END IF
112
113      IF( kt == nit000 )  THEN
114         IF(lwp)WRITE(numout,*)
115         IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing on ', cdtype
116         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ '
117         zavi = 0._wp      ! avoid warning at compilation phase when lk_ldfslp=F
118      ENDIF
119      !
120      ! I. Local initialization
121      ! -----------------------
122      zwd(1,:, : ) = 0._wp     ;     zwd(jpi,:,:) = 0._wp
123      zws(1,:, : ) = 0._wp     ;     zws(jpi,:,:) = 0._wp
124      zwi(1,:, : ) = 0._wp     ;     zwi(jpi,:,:) = 0._wp
125      zwt(1,:, : ) = 0._wp     ;     zwt(jpi,:,:) = 0._wp
126      zwt(:,:,jpk) = 0._wp     ;     zwt( : ,:,1) = 0._wp
127
128      ! I.1 Variable volume : to take into account vertical variable vertical scale factors
129      ! -------------------
130      IF( lk_vvl ) THEN   ;    znvvl = 1._wp
131      ELSE                ;    znvvl = 0._wp
132      ENDIF
133
134      ! II. Vertical trend associated with the vertical physics
135      ! =======================================================
136      !     (including the vertical flux proportional to dk[t] associated
137      !      with the lateral mixing, through the avt update)
138      !     dk[ avt dk[ (t,s) ] ] diffusive trends
139
140      !
141      ! II.0 Matrix construction
142      ! ------------------------
143      DO jn = 1, kjpt
144         !
145         !  Matrix construction
146         ! ------------------------
147         IF( cdtype == 'TRA' .AND. jn == jp_tem )  THEN 
148#if defined key_ldfslp
149            IF( ln_traldf_grif ) THEN
150               DO jk = 2, jpkm1
151                  DO jj = 2, jpjm1
152                     DO ji = fs_2, fs_jpim1   ! vector opt.
153                        ! zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk)         ! vertical mixing coef. due to lateral mixing
154                        zavi = ah_wslp2(ji,jj,jk)                ! vertical mixing coef. due to lateral mixing
155                        zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
156                     END DO
157                  END DO
158               END DO
159            ! update and save of avt (and avs if double diffusive mixing)
160            ELSE IF( l_traldf_rot ) THEN
161               DO jk = 2, jpkm1
162                  DO jj = 2, jpjm1
163                     DO ji = fs_2, fs_jpim1   ! vector opt.
164                        zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
165                           & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
166                           &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
167                        zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
168                     END DO
169                  END DO
170               END DO
171            ELSE                         ! no rotation but key_ldfslp defined
172               zwt(:,:,:) = avt(:,:,:)
173            ENDIF
174#else
175            ! No isopycnal diffusion
176            zwt(:,:,:) = avt(:,:,:)           
177#endif
178            ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
179            DO jk = 1, jpkm1
180               DO jj = 2, jpjm1
181                  DO ji = fs_2, fs_jpim1   ! vector opt.
182                     ze3ta =  ( 1. - znvvl ) +        znvvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point
183                     ze3tn =         znvvl   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point
184                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
185                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
186                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
187                 END DO
188               END DO
189            END DO
190            ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
191            DO jj = 2, jpjm1
192               DO ji = fs_2, fs_jpim1
193                  zwt(ji,jj,1) = zwd(ji,jj,1)
194               END DO
195            END DO
196            DO jk = 2, jpkm1
197               DO jj = 2, jpjm1
198                  DO ji = fs_2, fs_jpim1
199                    zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
200                  END DO
201               END DO
202            END DO
203            !
204         ELSE IF( ( cdtype == 'TRA' .AND. jn == jp_sal ) .OR. ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN
205#if defined key_ldfslp
206            IF( ln_traldf_grif ) THEN
207               DO jk = 2, jpkm1
208                  DO jj = 2, jpjm1
209                     DO ji = fs_2, fs_jpim1   ! vector opt.
210                        zavi = ah_wslp2(ji,jj,jk)                ! vertical mixing coef. due to lateral mixing
211                        ! zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk)         ! vertical mixing coef. due to lateral mixing
212                        zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
213                     END DO
214                  END DO
215               END DO
216            ELSE IF( l_traldf_rot ) THEN
217               DO jk = 2, jpkm1
218                  DO jj = 2, jpjm1
219                     DO ji = fs_2, fs_jpim1   ! vector opt.
220                        zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
221                           & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
222                           &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
223                        zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on salinity)
224                     END DO
225                  END DO
226               END DO
227            ELSE                         ! no rotation but key_ldfslp defined
228               zwt(:,:,:) = fsavs(:,:,:)
229            ENDIF
230#else
231            ! No isopycnal diffusion
232            zwt(:,:,:) = fsavs(:,:,:)           
233#endif
234            ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
235            DO jk = 1, jpkm1
236               DO jj = 2, jpjm1
237                  DO ji = fs_2, fs_jpim1   ! vector opt.
238                     ze3ta =  ( 1. - znvvl ) +        znvvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point
239                     ze3tn =         znvvl   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point
240                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
241                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
242                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
243                 END DO
244               END DO
245            END DO
246            ! Surface boudary conditions
247            DO jj = 2, jpjm1
248               DO ji = fs_2, fs_jpim1   ! vector opt.
249                 ze3ta = ( 1. - znvvl ) +  znvvl * fse3t_a(ji,jj,1)    ! after scale factor at T-point
250                 zwi(ji,jj,1) = 0._wp
251                 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
252               END DO
253            END DO
254            !
255            ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
256            DO jj = 2, jpjm1
257               DO ji = fs_2, fs_jpim1
258                  zwt(ji,jj,1) = zwd(ji,jj,1)
259               END DO
260            END DO
261            DO jk = 2, jpkm1
262               DO jj = 2, jpjm1
263                  DO ji = fs_2, fs_jpim1
264                    zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
265                  END DO
266               END DO
267            END DO
268            !
269         END IF 
270
271         ! II.1. Vertical diffusion on tracer
272         ! ---------------------------
273         
274         ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
275         DO jj = 2, jpjm1
276            DO ji = fs_2, fs_jpim1
277               ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,1)
278               ze3tn = ( 1. - znvvl ) + znvvl * fse3t(ji,jj,1)
279               pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) + p2dt(1) * ze3tn * pta(ji,jj,1,jn)
280            END DO
281         END DO
282         DO jk = 2, jpkm1
283            DO jj = 2, jpjm1
284               DO ji = fs_2, fs_jpim1
285                  ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,jk)
286                  ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,jk)
287                  zrhs = ze3tb * ptb(ji,jj,jk,jn) + p2dt(jk) * ze3tn * pta(ji,jj,jk,jn)   ! zrhs=right hand side
288                  pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn)
289               END DO
290            END DO
291         END DO
292
293         ! third recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk
294         DO jj = 2, jpjm1
295            DO ji = fs_2, fs_jpim1
296               pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
297            END DO
298         END DO
299         DO jk = jpk-2, 1, -1
300            DO jj = 2, jpjm1
301               DO ji = fs_2, fs_jpim1
302                  pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) )   &
303                     &             / zwt(ji,jj,jk) * tmask(ji,jj,jk)
304               END DO
305            END DO
306         END DO
307         !
308      END DO
309      !
310      IF(wrk_not_released(3, 1,2))THEN
311         CALL ctl_stop('tra_zdf_imp : failed to release workspace arrays.')
312      END IF
313      !
314   END SUBROUTINE tra_zdf_imp
315
316   !!==============================================================================
317END MODULE trazdf_imp
Note: See TracBrowser for help on using the repository browser.