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.
trczdf_imp.F90 in trunk/NEMO/TOP_SRC/TRP – NEMO

source: trunk/NEMO/TOP_SRC/TRP/trczdf_imp.F90 @ 1271

Last change on this file since 1271 was 1271, checked in by rblod, 15 years ago

Addapt AGRIF routines to the new TOP organization, clean some routines and add a sponge layer for passive tracers, see ticket #293
  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 12.2 KB
Line 
1MODULE trczdf_imp
2   !!==============================================================================
3   !!                    ***  MODULE  trczdf_imp  ***
4   !! Ocean passive tracers:  vertical component of the tracer mixing trend
5   !!==============================================================================
6   !! History :  6.0  !  90-10  (B. Blanke)  Original code
7   !!            7.0  !  91-11  (G. Madec)
8   !!                 !  92-06  (M. Imbard)  correction on tracer trend loops
9   !!                 !  96-01  (G. Madec)  statement function for e3
10   !!                 !  97-05  (G. Madec)  vertical component of isopycnal
11   !!                 !  97-07  (G. Madec)  geopotential diffusion in s-coord
12   !!                 !  98-03  (L. Bopp MA Foujols) passive tracer generalisation
13   !!                 !  00-05  (MA Foujols) add lbc for tracer trends
14   !!                 !  00-06  (O Aumont)  correct isopycnal scheme suppress
15   !!                 !                     avt multiple correction
16   !!                 !  00-08  (G. Madec)  double diffusive mixing
17   !!            8.5  !  02-08  (G. Madec)  F90: Free form and module
18   !!            9.0  !  04-03  (C. Ethe )  adapted for passive tracers
19   !!                 !  07-02  (C. Deltel)  Diagnose ML trends for passive tracers
20   !!----------------------------------------------------------------------
21#if defined key_top
22   !!----------------------------------------------------------------------
23   !!   'key_top'                                                TOP models
24   !!----------------------------------------------------------------------
25   !!   trc_zdf_imp  : update the tracer trend with the vertical diffusion
26   !!                  using an implicit time-stepping.
27   !!----------------------------------------------------------------------
28   !! * Modules used
29   USE oce_trc             ! ocean dynamics and active tracers variables
30   USE trp_trc                 ! ocean passive tracers variables
31   USE trctrp_lec          ! passive tracers transport
32   USE prtctl_trc
33   USE trdmld_trc
34   USE trdmld_trc_oce
35
36   IMPLICIT NONE
37   PRIVATE
38
39   !! * Routine accessibility
40   PUBLIC trc_zdf_imp          ! routine called by step.F90
41
42   !! * Module variable
43   REAL(wp), DIMENSION(jpk) ::   &
44      rdttrc                     ! vertical profile of 2 x tracer time-step
45
46   !! * Substitutions
47#  include "top_substitute.h90"
48   !!----------------------------------------------------------------------
49   !!   TOP 1.0 , LOCEAN-IPSL (2005)
50   !! $Id$
51   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
52   !!----------------------------------------------------------------------
53
54CONTAINS
55
56   SUBROUTINE trc_zdf_imp( kt )
57      !!----------------------------------------------------------------------
58      !!                  ***  ROUTINE trc_zdf_imp  ***
59      !!
60      !! ** Purpose :   Compute the trend due to the vertical tracer mixing
61      !!      using an implicit time stepping and add it to the general trend
62      !!      of the tracer equations.
63      !!
64      !! ** Method  :   The vertical diffusion of tracers tra is given by:
65      !!          difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(tra) )
66      !!      It is thus evaluated using a backward time scheme
67      !!      Surface and bottom boundary conditions: no diffusive flux on
68      !!      both tracers (bottom, applied through the masked field avt).
69      !!      Add this trend to the general trend tra :
70      !!          tra = tra + dz( avt dz(t) )
71      !!         (tra = tra + dz( avs dz(t) ) if lk_zdfddmtrc=T)
72      !!
73      !! ** Action  : - Update tra with the before vertical diffusion trend
74      !!              - save the trends
75      !!
76      !!---------------------------------------------------------------------
77      !!
78      !! * Arguments
79      USE oce, ONLY :   ztrtrd => ua    ! use ua as 3D workspace
80
81      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
82      INTEGER ::   ikst, ikenm2, ikstp1
83      !! * Local declarations
84      INTEGER ::   ji, jj, jk, jn             ! dummy loop indices
85      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
86         zwd, zws, zwi,          &  ! ???
87         zwx, zwy, zwt              ! ???
88#if defined key_trc_diatrd
89      REAL(wp) ::  ztra      ! temporary scalars
90#endif
91      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra) ::   &
92         ztrd
93      CHARACTER (len=22) :: charout
94      !!---------------------------------------------------------------------
95
96      IF( kt == nittrc000 ) THEN
97         WRITE(numout,*)
98         WRITE(numout,*) 'trc_zdf_implicit : vertical tracer mixing'
99         WRITE(numout,*) '~~~~~~~~~~~~~~~'
100      ENDIF
101
102      ! 0. Local constant initialization
103      ! --------------------------------
104      IF( ln_trcadv_cen2 .OR. ln_trcadv_tvd ) THEN
105         ! time step = 2 rdttra with Arakawa or TVD advection scheme
106         IF( neuler == 0 .AND. kt == nittrc000 ) THEN
107            rdttrc(:) =  rdttra(:) * FLOAT(ndttrc)             ! restarting with Euler time stepping
108         ELSEIF( kt <= nittrc000 + ndttrc ) THEN
109            rdttrc(:) = 2. * rdttra(:) * FLOAT(ndttrc)         ! leapfrog
110         ENDIF
111      ELSE
112         rdttrc(:) =  rdttra(:) * FLOAT(ndttrc)     
113      ENDIF
114     !                                                       ! ===========
115      DO jn = 1, jptra                                        ! tracer loop
116         !                                                    ! ===========
117         IF( l_trdtrc ) ztrtrd(:,:,:) = tra(:,:,:,jn)         ! ??? validation needed
118
119    ! Initialisation     
120    zwd( 1 ,:,:) = 0.e0     ;     zwd(jpi,:,:) = 0.e0
121    zws( 1 ,:,:) = 0.e0     ;     zws(jpi,:,:) = 0.e0
122    zwi( 1 ,:,:) = 0.e0     ;     zwi(jpi,:,:) = 0.e0
123    zwt( 1 ,:,:) = 0.e0     ;     zwt(jpi,:,:) = 0.e0     
124         zwt(  :,:,1) = 0.e0     ;     zwt(  :,:,jpk) = 0.e0
125         !                                         
126         ! 0. Matrix construction
127         ! ----------------------
128
129         ! Diagonal, inferior, superior
130         ! (including the bottom boundary condition via avs masked
131         DO jk = 1, jpkm1                                                     
132            DO jj = 2, jpjm1                                     
133               DO ji = fs_2, fs_jpim1   ! vector opt.
134                  zwi(ji,jj,jk) = - rdttrc(jk) * fstravs(ji,jj,jk  ) /( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
135                  zws(ji,jj,jk) = - rdttrc(jk) * fstravs(ji,jj,jk+1) /( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
136                  zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk)
137               END DO
138            END DO
139         END DO
140
141         ! Surface boudary conditions
142         DO jj = 2, jpjm1       
143            DO ji = fs_2, fs_jpim1
144               zwi(ji,jj,1) = 0.e0
145               zwd(ji,jj,1) = 1. - zws(ji,jj,1)
146            END DO
147         END DO
148         
149         ! Second member construction
150         DO jk = 1, jpkm1
151            DO jj = 2, jpjm1     
152               DO ji = fs_2, fs_jpim1
153                  zwy(ji,jj,jk) = trb(ji,jj,jk,jn) + rdttrc(jk) * tra(ji,jj,jk,jn)
154               END DO
155            END DO
156         END DO
157         
158         
159         ! Matrix inversion from the first level
160         !----------------------------------------------------------------------
161         !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
162         !
163         !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
164         !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
165         !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
166         !        (        ...               )( ...  ) ( ...  )
167         !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
168         !
169         !   m is decomposed in the product of an upper and lower triangular matrix
170         !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
171         !   The second member is in 2d array zwy
172         !   The solution is in 2d array zwx
173         !   The 3d arry zwt is a work space array
174         !   zwy is used and then used as a work space array : its value is modified!
175         !
176         !   N.B. the starting vertical index (ikst) is equal to 1 except for
177         !   the resolution of tke matrix where surface tke value is prescribed
178         !   so that ikstrt=2.
179        ikst = 1
180        ikstp1 = ikst + 1
181        ikenm2 = jpk - 2
182       
183        DO jj = 2, jpjm1
184           DO ji = fs_2, fs_jpim1
185              zwt(ji,jj,ikst) = zwd(ji,jj,ikst)
186           END DO
187        END DO
188       
189        DO jk = ikstp1, jpkm1
190           DO jj = 2, jpjm1
191              DO ji = fs_2, fs_jpim1
192                 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
193              END DO
194           END DO
195        END DO
196       
197        DO jk = ikstp1, jpkm1
198           DO jj = 2, jpjm1
199              DO ji = fs_2, fs_jpim1
200                 zwy(ji,jj,jk) = zwy(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * zwy(ji,jj,jk-1)
201              END DO
202           END DO
203        END DO
204       
205        DO jj = 2, jpjm1
206           DO ji = fs_2, fs_jpim1
207              zwx(ji,jj,jpkm1) = zwy(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
208           END DO
209        END DO
210       
211        DO jk = ikenm2, ikst, -1
212           DO jj = 2, jpjm1
213              DO ji = fs_2, fs_jpim1
214                 zwx(ji,jj,jk) = ( zwy(ji,jj,jk) - zws(ji,jj,jk) * zwx(ji,jj,jk+1) ) / zwt(ji,jj,jk)
215              END DO
216           END DO
217        END DO
218
219#if defined key_trc_diatrd
220         ! Compute and save the vertical diffusive of tracers trends
221#  if defined key_trcldf_iso
222         DO jk = 1, jpkm1
223            DO jj = 2, jpjm1
224               DO ji = fs_2, fs_jpim1   ! vector opt.
225                  ztra = ( zwx(ji,jj,jk) - trb(ji,jj,jk,jn) ) / rdttrc(jk)
226                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),6) = ztra - tra(ji,jj,jk,jn) + trtrd(ji,jj,jk,ikeep(jn),6)
227               END DO
228            END DO
229         END DO
230#  else
231         DO jk = 1, jpkm1
232            DO jj = 2, jpjm1
233               DO ji = fs_2, fs_jpim1   ! vector opt.
234                  ztra = ( zwx(ji,jj,jk) - trb(ji,jj,jk,jn) ) / rdttrc(jk)
235                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),6) = ztra - tra(ji,jj,jk,jn)
236               END DO
237            END DO
238         END DO
239#  endif
240#endif
241
242         ! Save the masked passive tracer after in tra
243         ! (c a u t i o n:  tracer not its trend, Leap-frog scheme done
244         !                  it will not be done in tranxt)
245         DO jk = 1, jpkm1
246            DO jj = 2, jpjm1
247               DO ji = fs_2, fs_jpim1
248                  tra(ji,jj,jk,jn) = zwx(ji,jj,jk) * tmask(ji,jj,jk)
249               END DO
250            END DO
251         END DO
252
253         IF( l_trdtrc ) THEN                            ! trends
254            DO jk = 1, jpkm1
255               ztrtrd(:,:,jk) = ( ( tra(:,:,jk,jn) - trb(:,:,jk,jn) ) / rdttrc(jk) ) - ztrtrd(:,:,jk)
256            END DO
257            IF (luttrd(jn)) CALL trd_mod_trc(ztrtrd, jn, jptrc_trd_zdf, kt)
258         END IF
259
260         IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
261            ztrd(:,:,:,:) = 0.
262            DO jk = 1, jpkm1
263               DO jj = 2, jpjm1
264                  DO ji = fs_2, fs_jpim1
265                     ztrd(ji,jj,jk,jn) = ( zwx(ji,jj,jk) - trb(ji,jj,jk,jn) ) / rdttrc(jk)
266                  END DO
267               END DO
268            END DO
269         ENDIF
270
271         !                                                    ! ===========
272      END DO                                                  ! tracer loop
273      !                                                       ! ===========
274
275      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
276         WRITE(charout, FMT="('zdf - imp')")
277         CALL prt_ctl_trc_info(charout)
278         CALL prt_ctl_trc(tab4d=ztrd, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
279      ENDIF
280
281   END SUBROUTINE trc_zdf_imp
282
283#else
284   !!----------------------------------------------------------------------
285   !!   Dummy module :                      NO passive tracer
286   !!----------------------------------------------------------------------
287CONTAINS
288   SUBROUTINE trc_zdf_imp (kt )              ! Empty routine
289      INTEGER, INTENT(in) :: kt
290      WRITE(*,*) 'trc_zdf_imp: You should not have seen this print! error?', kt
291   END SUBROUTINE trc_zdf_imp
292#endif
293   
294!!==============================================================================
295END MODULE trczdf_imp
Note: See TracBrowser for help on using the repository browser.