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

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

Phasing branch dev_r2586_dynamic_mem with revision 2675 off the trunk
  • Property svn:keywords set to Id
File size: 12.1 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   !!             -   !  2011-02  (A. Coward, C. Ethe, G. Madec) improvment of surface boundary condition
18   !!----------------------------------------------------------------------
19 
20   !!----------------------------------------------------------------------
21   !!   tra_zdf_imp : Update the tracer trend with the diagonal vertical 
22   !!                 part of the mixing tensor.
23   !!----------------------------------------------------------------------
24   USE oce             ! ocean dynamics and tracers variables
25   USE dom_oce         ! ocean space and time domain variables
26   USE zdf_oce         ! ocean vertical physics variables
27   USE trc_oce         ! share passive tracers/ocean variables
28   USE domvvl          ! variable volume
29   USE ldftra_oce      ! ocean active tracers: lateral physics
30   USE ldftra          ! lateral mixing type
31   USE ldfslp          ! lateral physics: slope of diffusion
32   USE zdfddm          ! ocean vertical physics: double diffusion
33   USE traldf_iso_grif ! active tracers: Griffies operator
34   USE in_out_manager  ! I/O manager
35   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
36   USE lib_mpp         ! MPP library
37
38   IMPLICIT NONE
39   PRIVATE
40
41   PUBLIC   tra_zdf_imp   !  routine called by step.F90
42
43   REAL(wp) ::  r_vvl     ! variable volume indicator, =1 if lk_vvl=T, =0 otherwise
44
45   !! * Substitutions
46#  include "domzgr_substitute.h90"
47#  include "ldftra_substitute.h90"
48#  include "zdfddm_substitute.h90"
49#  include "vectopt_loop_substitute.h90"
50   !!----------------------------------------------------------------------
51   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
52   !! $Id$
53   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
54   !!----------------------------------------------------------------------
55CONTAINS
56 
57   SUBROUTINE tra_zdf_imp( kt, cdtype, p2dt, ptb, pta, kjpt ) 
58      !!----------------------------------------------------------------------
59      !!                  ***  ROUTINE tra_zdf_imp  ***
60      !!
61      !! ** Purpose :   Compute the after tracer through a implicit computation
62      !!     of the vertical tracer diffusion (including the vertical component
63      !!     of lateral mixing (only for 2nd order operator, for fourth order
64      !!     it is already computed and add to the general trend in traldf)
65      !!
66      !! ** Method  :  The vertical diffusion of the tracer t  is given by:
67      !!                  difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
68      !!      It is computed using a backward time scheme (t=ta).
69      !!      If lk_zdfddm=T, use avs for salinity or for passive tracers
70      !!      Surface and bottom boundary conditions: no diffusive flux on
71      !!      both tracers (bottom, applied through the masked field avt).
72      !!      If iso-neutral mixing, add to avt the contribution due to lateral mixing.
73      !!
74      !! ** Action  : - pta  becomes the after tracer
75      !!---------------------------------------------------------------------
76      USE oce    , ONLY :   zwd   => ua   ! ua used as workspace
77      USE oce    , ONLY :   zws   => va   ! va  -          -
78      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
79      USE wrk_nemo, ONLY: zwi => wrk_3d_1, zwt => wrk_3d_2  ! workspace arrays
80      !!
81      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
82      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
83      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
84      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt     ! vertical profile of tracer time-step
85      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb      ! before and now tracer fields
86      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend
87      !!
88      INTEGER  ::  ji, jj, jk, jn        ! dummy loop indices
89      REAL(wp) ::  zrhs                  ! local scalars
90      REAL(wp) ::  ze3tb, ze3tn, ze3ta   ! variable vertical scale factors
91      !!---------------------------------------------------------------------
92
93      IF(wrk_in_use(3, 1,2))THEN
94         CALL ctl_stop('tra_zdf_imp : requested workspace arrays unavailable.')
95         RETURN
96      END IF
97
98      IF( kt == nit000 )  THEN
99         IF(lwp)WRITE(numout,*)
100         IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing on ', cdtype
101         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ '
102         !
103         IF( lk_vvl ) THEN   ;    r_vvl = 1._wp       ! Variable volume indicator
104         ELSE                ;    r_vvl = 0._wp       
105         ENDIF
106      ENDIF
107      !
108      !                                               ! ============= !
109      DO jn = 1, kjpt                                 !  tracer loop  !
110         !                                            ! ============= !
111         !
112         !  Matrix construction
113         ! --------------------
114         ! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer
115         !
116         IF(  ( cdtype == 'TRA' .AND. ( ( jn == jp_tem ) .OR. ( jn == jp_sal .AND. lk_zdfddm ) ) ) .OR. &
117            & ( cdtype == 'TRC' .AND. jn == 1 )  )  THEN
118            !
119            ! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers
120            IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN   ;   zwt(:,:,2:jpk) = avt  (:,:,2:jpk)
121            ELSE                                            ;   zwt(:,:,2:jpk) = fsavs(:,:,2:jpk)
122            ENDIF
123            zwt(:,:,1) = 0._wp
124            !
125#if defined key_ldfslp
126            ! isoneutral diffusion: add the contribution
127            IF( ln_traldf_grif    ) THEN     ! Griffies isoneutral diff
128               DO jk = 2, jpkm1
129                  DO jj = 2, jpjm1
130                     DO ji = fs_2, fs_jpim1   ! vector opt.
131                        zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk)       
132                     END DO
133                  END DO
134               END DO
135            ELSE IF( l_traldf_rot ) THEN     ! standard isoneutral diff
136               DO jk = 2, jpkm1
137                  DO jj = 2, jpjm1
138                     DO ji = fs_2, fs_jpim1   ! vector opt.
139                        zwt(ji,jj,jk) = zwt(ji,jj,jk) + fsahtw(ji,jj,jk)                       &
140                           &                          * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
141                           &                             + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
142                     END DO
143                  END DO
144               END DO
145            ENDIF
146#endif
147            ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked)
148            DO jk = 1, jpkm1
149               DO jj = 2, jpjm1
150                  DO ji = fs_2, fs_jpim1   ! vector opt.
151                     ze3ta =  ( 1. - r_vvl ) +        r_vvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point
152                     ze3tn =         r_vvl   + ( 1. - r_vvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point
153                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
154                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
155                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
156                 END DO
157               END DO
158            END DO
159            !
160            !! Matrix inversion from the first level
161            !!----------------------------------------------------------------------
162            !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
163            !
164            !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
165            !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
166            !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
167            !        (        ...               )( ...  ) ( ...  )
168            !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
169            !
170            !   m is decomposed in the product of an upper and lower triangular matrix.
171            !   The 3 diagonal terms are in 3d arrays: zwd, zws, zwi.
172            !   Suffices i,s and d indicate "inferior" (below diagonal), diagonal
173            !   and "superior" (above diagonal) components of the tridiagonal system.
174            !   The solution will be in the 4d array pta.
175            !   The 3d array zwt is used as a work space array.
176            !   En route to the solution pta is used a to evaluate the rhs and then
177            !   used as a work space array: its value is modified.
178            !
179            ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
180            ! done once for all passive tracers (so included in the IF instruction)
181            DO jj = 2, jpjm1
182               DO ji = fs_2, fs_jpim1
183                  zwt(ji,jj,1) = zwd(ji,jj,1)
184               END DO
185            END DO
186            DO jk = 2, jpkm1
187               DO jj = 2, jpjm1
188                  DO ji = fs_2, fs_jpim1
189                    zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
190                  END DO
191               END DO
192            END DO
193            !
194         END IF 
195         !         
196         ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
197         DO jj = 2, jpjm1
198            DO ji = fs_2, fs_jpim1
199               ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,1)
200               ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t(ji,jj,1)
201               pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) + p2dt(1) * ze3tn * pta(ji,jj,1,jn)
202            END DO
203         END DO
204         DO jk = 2, jpkm1
205            DO jj = 2, jpjm1
206               DO ji = fs_2, fs_jpim1
207                  ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,jk)
208                  ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t  (ji,jj,jk)
209                  zrhs = ze3tb * ptb(ji,jj,jk,jn) + p2dt(jk) * ze3tn * pta(ji,jj,jk,jn)   ! zrhs=right hand side
210                  pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn)
211               END DO
212            END DO
213         END DO
214
215         ! third recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk   (result is the after tracer)
216         DO jj = 2, jpjm1
217            DO ji = fs_2, fs_jpim1
218               pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
219            END DO
220         END DO
221         DO jk = jpk-2, 1, -1
222            DO jj = 2, jpjm1
223               DO ji = fs_2, fs_jpim1
224                  pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) )   &
225                     &             / zwt(ji,jj,jk) * tmask(ji,jj,jk)
226               END DO
227            END DO
228         END DO
229         !                                            ! ================= !
230      END DO                                          !  end tracer loop  !
231      !                                               ! ================= !
232      !
233      IF(wrk_not_released(3, 1,2))THEN
234         CALL ctl_stop('tra_zdf_imp : failed to release workspace arrays.')
235      END IF
236      !
237   END SUBROUTINE tra_zdf_imp
238
239   !!==============================================================================
240END MODULE trazdf_imp
Note: See TracBrowser for help on using the repository browser.