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

Last change on this file since 2636 was 2636, checked in by gm, 13 years ago

dynamic mem: #785 ; move ctl_stop & warn in lib_mpp to avoid a circular dependency + ctl_stop improvment
  • Property svn:keywords set to Id
File size: 15.4 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   USE lib_mpp         ! MPP library
36
37   IMPLICIT NONE
38   PRIVATE
39
40   PUBLIC   tra_zdf_imp   !  routine called by step.F90
41
42   !! * Substitutions
43#  include "domzgr_substitute.h90"
44#  include "ldftra_substitute.h90"
45#  include "zdfddm_substitute.h90"
46#  include "vectopt_loop_substitute.h90"
47   !!----------------------------------------------------------------------
48   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
49   !! $Id$
50   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
51   !!----------------------------------------------------------------------
52CONTAINS
53 
54   SUBROUTINE tra_zdf_imp( kt, cdtype, p2dt, ptb, pta, kjpt ) 
55      !!----------------------------------------------------------------------
56      !!                  ***  ROUTINE tra_zdf_imp  ***
57      !!
58      !! ** Purpose :   Compute the trend due to the vertical tracer diffusion
59      !!     including the vertical component of lateral mixing (only for 2nd
60      !!     order operator, for fourth order it is already computed and add
61      !!     to the general trend in traldf.F) and add it to the general trend
62      !!     of the tracer equations.
63      !!
64      !! ** Method  :   The vertical component of the lateral diffusive trends
65      !!      is provided by a 2nd order operator rotated along neutral or geo-
66      !!      potential surfaces to which an eddy induced advection can be
67      !!      added. It is computed using before fields (forward in time) and
68      !!      isopycnal or geopotential slopes computed in routine ldfslp.
69      !!
70      !!      Second part: vertical trend associated with the vertical physics
71      !!      ===========  (including the vertical flux proportional to dk[t]
72      !!                  associated with the lateral mixing, through the
73      !!                  update of avt)
74      !!      The vertical diffusion of the tracer t  is given by:
75      !!             difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
76      !!      It is computed using a backward time scheme (t=ta).
77      !!      Surface and bottom boundary conditions: no diffusive flux on
78      !!      both tracers (bottom, applied through the masked field avt).
79      !!      Add this trend to the general trend ta,sa :
80      !!         ta = ta + dz( avt dz(t) )
81      !!         if lk_zdfddm=T, use avs for salinity or for passive tracers
82      !!         (sa = sa + dz( avs dz(t) )
83      !!
84      !!      Third part: recover avt resulting from the vertical physics
85      !!      ==========  alone, for further diagnostics (for example to
86      !!                  compute the turbocline depth in zdfmxl.F90).
87      !!         if lk_zdfddm=T, use avt = zavt
88      !!         (avs = zavs if lk_zdfddm=T )
89      !!
90      !! ** Action  : - Update (ta) with before vertical diffusion trend
91      !!---------------------------------------------------------------------
92      USE oce    , ONLY :   zwd   => ua   ! ua used as workspace
93      USE oce    , ONLY :   zws   => va   ! va  -          -
94      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
95      USE wrk_nemo, ONLY: zwi => wrk_3d_1, zwt => wrk_3d_2  ! workspace arrays
96      !!
97      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
98      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
99      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
100      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt     ! vertical profile of tracer time-step
101      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb      ! before and now tracer fields
102      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend
103      !!
104      INTEGER  ::  ji, jj, jk, jn        ! dummy loop indices
105      REAL(wp) ::  zavi, zrhs, znvvl     ! local scalars
106      REAL(wp) ::  ze3tb, ze3tn, ze3ta   ! variable vertical scale factors
107      !!---------------------------------------------------------------------
108
109      IF(wrk_in_use(3, 1,2))THEN
110         CALL ctl_stop('tra_zdf_imp : requested workspace arrays unavailable.')
111         RETURN
112      END IF
113
114      IF( kt == nit000 )  THEN
115         IF(lwp)WRITE(numout,*)
116         IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing on ', cdtype
117         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ '
118         zavi = 0._wp      ! avoid warning at compilation phase when lk_ldfslp=F
119      ENDIF
120      !
121      ! I. Local initialization
122      ! -----------------------
123      zwd(1,:, : ) = 0._wp     ;     zwd(jpi,:,:) = 0._wp
124      zws(1,:, : ) = 0._wp     ;     zws(jpi,:,:) = 0._wp
125      zwi(1,:, : ) = 0._wp     ;     zwi(jpi,:,:) = 0._wp
126      zwt(1,:, : ) = 0._wp     ;     zwt(jpi,:,:) = 0._wp
127      zwt(:,:,jpk) = 0._wp     ;     zwt( : ,:,1) = 0._wp
128
129      ! I.1 Variable volume : to take into account vertical variable vertical scale factors
130      ! -------------------
131      IF( lk_vvl ) THEN   ;    znvvl = 1._wp
132      ELSE                ;    znvvl = 0._wp
133      ENDIF
134
135      ! II. Vertical trend associated with the vertical physics
136      ! =======================================================
137      !     (including the vertical flux proportional to dk[t] associated
138      !      with the lateral mixing, through the avt update)
139      !     dk[ avt dk[ (t,s) ] ] diffusive trends
140
141      !
142      ! II.0 Matrix construction
143      ! ------------------------
144      DO jn = 1, kjpt
145         !
146         !  Matrix construction
147         ! ------------------------
148         IF( cdtype == 'TRA' .AND. jn == jp_tem )  THEN 
149#if defined key_ldfslp
150            IF( ln_traldf_grif ) THEN
151               DO jk = 2, jpkm1
152                  DO jj = 2, jpjm1
153                     DO ji = fs_2, fs_jpim1   ! vector opt.
154                        ! zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk)         ! vertical mixing coef. due to lateral mixing
155                        zavi = ah_wslp2(ji,jj,jk)                ! vertical mixing coef. due to lateral mixing
156                        zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
157                     END DO
158                  END DO
159               END DO
160            ! update and save of avt (and avs if double diffusive mixing)
161            ELSE IF( l_traldf_rot ) THEN
162               DO jk = 2, jpkm1
163                  DO jj = 2, jpjm1
164                     DO ji = fs_2, fs_jpim1   ! vector opt.
165                        zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
166                           & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
167                           &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
168                        zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
169                     END DO
170                  END DO
171               END DO
172            ELSE                         ! no rotation but key_ldfslp defined
173               zwt(:,:,:) = avt(:,:,:)
174            ENDIF
175#else
176            ! No isopycnal diffusion
177            zwt(:,:,:) = avt(:,:,:)           
178#endif
179            ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
180            DO jk = 1, jpkm1
181               DO jj = 2, jpjm1
182                  DO ji = fs_2, fs_jpim1   ! vector opt.
183                     ze3ta =  ( 1. - znvvl ) +        znvvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point
184                     ze3tn =         znvvl   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point
185                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
186                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
187                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
188                 END DO
189               END DO
190            END DO
191            ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
192            DO jj = 2, jpjm1
193               DO ji = fs_2, fs_jpim1
194                  zwt(ji,jj,1) = zwd(ji,jj,1)
195               END DO
196            END DO
197            DO jk = 2, jpkm1
198               DO jj = 2, jpjm1
199                  DO ji = fs_2, fs_jpim1
200                    zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
201                  END DO
202               END DO
203            END DO
204            !
205         ELSE IF( ( cdtype == 'TRA' .AND. jn == jp_sal ) .OR. ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN
206#if defined key_ldfslp
207            IF( ln_traldf_grif ) THEN
208               DO jk = 2, jpkm1
209                  DO jj = 2, jpjm1
210                     DO ji = fs_2, fs_jpim1   ! vector opt.
211                        zavi = ah_wslp2(ji,jj,jk)                ! vertical mixing coef. due to lateral mixing
212                        ! zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk)         ! vertical mixing coef. due to lateral mixing
213                        zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
214                     END DO
215                  END DO
216               END DO
217            ELSE IF( l_traldf_rot ) THEN
218               DO jk = 2, jpkm1
219                  DO jj = 2, jpjm1
220                     DO ji = fs_2, fs_jpim1   ! vector opt.
221                        zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
222                           & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
223                           &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
224                        zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on salinity)
225                     END DO
226                  END DO
227               END DO
228            ELSE                         ! no rotation but key_ldfslp defined
229               zwt(:,:,:) = fsavs(:,:,:)
230            ENDIF
231#else
232            ! No isopycnal diffusion
233            zwt(:,:,:) = fsavs(:,:,:)           
234#endif
235            ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
236            DO jk = 1, jpkm1
237               DO jj = 2, jpjm1
238                  DO ji = fs_2, fs_jpim1   ! vector opt.
239                     ze3ta =  ( 1. - znvvl ) +        znvvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point
240                     ze3tn =         znvvl   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point
241                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
242                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
243                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
244                 END DO
245               END DO
246            END DO
247            ! Surface boudary conditions
248            DO jj = 2, jpjm1
249               DO ji = fs_2, fs_jpim1   ! vector opt.
250                 ze3ta = ( 1. - znvvl ) +  znvvl * fse3t_a(ji,jj,1)    ! after scale factor at T-point
251                 zwi(ji,jj,1) = 0._wp
252                 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
253               END DO
254            END DO
255            !
256            ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
257            DO jj = 2, jpjm1
258               DO ji = fs_2, fs_jpim1
259                  zwt(ji,jj,1) = zwd(ji,jj,1)
260               END DO
261            END DO
262            DO jk = 2, jpkm1
263               DO jj = 2, jpjm1
264                  DO ji = fs_2, fs_jpim1
265                    zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
266                  END DO
267               END DO
268            END DO
269            !
270         END IF 
271
272         ! II.1. Vertical diffusion on tracer
273         ! ---------------------------
274         
275         ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
276         DO jj = 2, jpjm1
277            DO ji = fs_2, fs_jpim1
278               ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,1)
279               ze3tn = ( 1. - znvvl ) + znvvl * fse3t(ji,jj,1)
280               pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) + p2dt(1) * ze3tn * pta(ji,jj,1,jn)
281            END DO
282         END DO
283         DO jk = 2, jpkm1
284            DO jj = 2, jpjm1
285               DO ji = fs_2, fs_jpim1
286                  ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,jk)
287                  ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,jk)
288                  zrhs = ze3tb * ptb(ji,jj,jk,jn) + p2dt(jk) * ze3tn * pta(ji,jj,jk,jn)   ! zrhs=right hand side
289                  pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn)
290               END DO
291            END DO
292         END DO
293
294         ! third recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk
295         DO jj = 2, jpjm1
296            DO ji = fs_2, fs_jpim1
297               pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
298            END DO
299         END DO
300         DO jk = jpk-2, 1, -1
301            DO jj = 2, jpjm1
302               DO ji = fs_2, fs_jpim1
303                  pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) )   &
304                     &             / zwt(ji,jj,jk) * tmask(ji,jj,jk)
305               END DO
306            END DO
307         END DO
308         !
309      END DO
310      !
311      IF(wrk_not_released(3, 1,2))THEN
312         CALL ctl_stop('tra_zdf_imp : failed to release workspace arrays.')
313      END IF
314      !
315   END SUBROUTINE tra_zdf_imp
316
317   !!==============================================================================
318END MODULE trazdf_imp
Note: See TracBrowser for help on using the repository browser.