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.F90 in NEMO/branches/2020/tickets_icb_1900/src/OCE/TRA – NEMO

source: NEMO/branches/2020/tickets_icb_1900/src/OCE/TRA/trazdf.F90 @ 14012

Last change on this file since 14012 was 14012, checked in by mathiot, 3 years ago

ticket 1900: upgrade branch to 13991 before the next upgrade
  • Property svn:keywords set to Id
File size: 13.5 KB
Line 
1MODULE trazdf
2   !!==============================================================================
3   !!                 ***  MODULE  trazdf  ***
4   !! Ocean active tracers:  vertical component of the tracer mixing trend
5   !!==============================================================================
6   !! History :  1.0  !  2005-11  (G. Madec)  Original code
7   !!            3.0  !  2008-01  (C. Ethe, G. Madec)  merge TRC-TRA
8   !!            4.0  !  2017-06  (G. Madec)  remove explict time-stepping option
9   !!----------------------------------------------------------------------
10
11   !!----------------------------------------------------------------------
12   !!   tra_zdf       : Update the tracer trend with the vertical diffusion
13   !!----------------------------------------------------------------------
14   USE oce            ! ocean dynamics and tracers variables
15   USE dom_oce        ! ocean space and time domain variables
16   USE domvvl         ! variable volume
17   USE phycst         ! physical constant
18   USE zdf_oce        ! ocean vertical physics variables
19   USE sbc_oce        ! surface boundary condition: ocean
20   USE ldftra         ! lateral diffusion: eddy diffusivity
21   USE ldfslp         ! lateral diffusion: iso-neutral slope
22   USE trd_oce        ! trends: ocean variables
23   USE trdtra         ! trends: tracer trend manager
24   !
25   USE in_out_manager ! I/O manager
26   USE prtctl         ! Print control
27   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
28   USE lib_mpp        ! MPP library
29   USE timing         ! Timing
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   tra_zdf       ! called by step.F90
35   PUBLIC   tra_zdf_imp   ! called by trczdf.F90
36
37   !! * Substitutions
38#  include "do_loop_substitute.h90"
39#  include "domzgr_substitute.h90"
40   !!----------------------------------------------------------------------
41   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
42   !! $Id$
43   !! Software governed by the CeCILL license (see ./LICENSE)
44   !!----------------------------------------------------------------------
45CONTAINS
46
47   SUBROUTINE tra_zdf( kt, Kbb, Kmm, Krhs, pts, Kaa )
48      !!----------------------------------------------------------------------
49      !!                  ***  ROUTINE tra_zdf  ***
50      !!
51      !! ** Purpose :   compute the vertical ocean tracer physics.
52      !!---------------------------------------------------------------------
53      INTEGER                                  , INTENT(in)    :: kt                  ! ocean time-step index
54      INTEGER                                  , INTENT(in)    :: Kbb, Kmm, Krhs, Kaa ! time level indices
55      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts                 ! active tracers and RHS of tracer equation
56      !
57      INTEGER  ::   ji, jj, jk   ! Dummy loop indices
58      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdt, ztrds   ! 3D workspace
59      !!---------------------------------------------------------------------
60      !
61      IF( ln_timing )   CALL timing_start('tra_zdf')
62      !
63      IF( kt == nit000 )  THEN
64         IF( ntile == 0 .OR. ntile == 1 )  THEN                   ! Do only on the first tile
65            IF(lwp)WRITE(numout,*)
66            IF(lwp)WRITE(numout,*) 'tra_zdf : implicit vertical mixing on T & S'
67            IF(lwp)WRITE(numout,*) '~~~~~~~ '
68         ENDIF
69      ENDIF
70      !
71      IF( l_trdtra )   THEN                  !* Save ta and sa trends
72         ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) )
73         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kaa)
74         ztrds(:,:,:) = pts(:,:,:,jp_sal,Kaa)
75      ENDIF
76      !
77      !                                      !* compute lateral mixing trend and add it to the general trend
78      CALL tra_zdf_imp( kt, nit000, 'TRA', rDt, Kbb, Kmm, Krhs, pts, Kaa, jpts ) 
79
80!!gm WHY here !   and I don't like that !
81      ! DRAKKAR SSS control {
82      ! JMM avoid negative salinities near river outlet ! Ugly fix
83      ! JMM : restore negative salinities to small salinities:
84      WHERE( pts(A2D(0),:,jp_sal,Kaa) < 0._wp )   pts(A2D(0),:,jp_sal,Kaa) = 0.1_wp
85!!gm
86
87      IF( l_trdtra )   THEN                      ! save the vertical diffusive trends for further diagnostics
88         DO jk = 1, jpk
89            ztrdt(:,:,jk) = (   (  pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa)     &
90               &                 - pts(:,:,jk,jp_tem,Kbb)*e3t(:,:,jk,Kbb)  )  &
91               &              / (  e3t(:,:,jk,Kmm)*rDt  )   )                 &
92               &          - ztrdt(:,:,jk)
93            ztrds(:,:,jk) = (   (  pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa)     &
94               &                 - pts(:,:,jk,jp_sal,Kbb)*e3t(:,:,jk,Kbb)  )  &
95               &             / (   e3t(:,:,jk,Kmm)*rDt  )   )                 &
96               &          - ztrds(:,:,jk)
97         END DO
98         ! NOTE: [tiling-comms-merge] The diagnostic results change along the north fold if this is removed
99!!gm this should be moved in trdtra.F90 and done on all trends
100         CALL lbc_lnk_multi( 'trazdf', ztrdt, 'T', 1.0_wp , ztrds, 'T', 1.0_wp )
101!!gm
102         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_zdf, ztrdt )
103         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_zdf, ztrds )
104         DEALLOCATE( ztrdt , ztrds )
105      ENDIF
106      !                                          ! print mean trends (used for debugging)
107      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Kaa), clinfo1=' zdf  - Ta: ', mask1=tmask,               &
108         &                                  tab3d_2=pts(:,:,:,jp_sal,Kaa), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
109      !
110      IF( ln_timing )   CALL timing_stop('tra_zdf')
111      !
112   END SUBROUTINE tra_zdf
113
114 
115   SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, Kbb, Kmm, Krhs, pt, Kaa, kjpt ) 
116      !!----------------------------------------------------------------------
117      !!                  ***  ROUTINE tra_zdf_imp  ***
118      !!
119      !! ** Purpose :   Compute the after tracer through a implicit computation
120      !!     of the vertical tracer diffusion (including the vertical component
121      !!     of lateral mixing (only for 2nd order operator, for fourth order
122      !!     it is already computed and add to the general trend in traldf)
123      !!
124      !! ** Method  :  The vertical diffusion of a tracer ,t , is given by:
125      !!          difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
126      !!      It is computed using a backward time scheme (t=after field)
127      !!      which provide directly the after tracer field.
128      !!      If ln_zdfddm=T, use avs for salinity or for passive tracers
129      !!      Surface and bottom boundary conditions: no diffusive flux on
130      !!      both tracers (bottom, applied through the masked field avt).
131      !!      If iso-neutral mixing, add to avt the contribution due to lateral mixing.
132      !!
133      !! ** Action  : - pt(:,:,:,:,Kaa)  becomes the after tracer
134      !!---------------------------------------------------------------------
135      INTEGER                                  , INTENT(in   ) ::   kt       ! ocean time-step index
136      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs, Kaa  ! ocean time level indices
137      INTEGER                                  , INTENT(in   ) ::   kit000   ! first time step index
138      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
139      INTEGER                                  , INTENT(in   ) ::   kjpt     ! number of tracers
140      REAL(wp)                                 , INTENT(in   ) ::   p2dt     ! tracer time-step
141      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt       ! tracers and RHS of tracer equation
142      !
143      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
144      REAL(wp) ::  zrhs, zzwi, zzws ! local scalars
145      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::  zwi, zwt, zwd, zws
146      !!---------------------------------------------------------------------
147      !
148      !                                               ! ============= !
149      DO jn = 1, kjpt                                 !  tracer loop  !
150         !                                            ! ============= !
151         !  Matrix construction
152         ! --------------------
153         ! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer
154         !
155         IF(  ( cdtype == 'TRA' .AND. ( jn == jp_tem .OR. ( jn == jp_sal .AND. ln_zdfddm ) ) ) .OR.   &
156            & ( cdtype == 'TRC' .AND. jn == 1 )  )  THEN
157            !
158            ! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers
159            IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN
160               DO_3D( 1, 1, 1, 1, 2, jpk )
161                  zwt(ji,jj,jk) = avt(ji,jj,jk)
162               END_3D
163            ELSE
164               DO_3D( 1, 1, 1, 1, 2, jpk )
165                  zwt(ji,jj,jk) = avs(ji,jj,jk)
166               END_3D
167            ENDIF
168            zwt(:,:,1) = 0._wp
169            !
170            IF( l_ldfslp ) THEN            ! isoneutral diffusion: add the contribution
171               IF( ln_traldf_msc  ) THEN     ! MSC iso-neutral operator
172                  DO_3D( 0, 0, 0, 0, 2, jpkm1 )
173                     zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk) 
174                  END_3D
175               ELSE                          ! standard or triad iso-neutral operator
176                  DO_3D( 0, 0, 0, 0, 2, jpkm1 )
177                     zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk)
178                  END_3D
179               ENDIF
180            ENDIF
181            !
182            ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked)
183            IF( ln_zad_Aimp ) THEN         ! Adaptive implicit vertical advection
184               DO_3D( 0, 0, 0, 0, 1, jpkm1 )
185                  zzwi = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk  ,Kmm)
186                  zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
187                  zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws   &
188                     &                 + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) )
189                  zwi(ji,jj,jk) = zzwi + p2dt *   MIN( wi(ji,jj,jk  ) , 0._wp )
190                  zws(ji,jj,jk) = zzws - p2dt *   MAX( wi(ji,jj,jk+1) , 0._wp )
191               END_3D
192            ELSE
193               DO_3D( 0, 0, 0, 0, 1, jpkm1 )
194                  zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk,Kmm)
195                  zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
196                  zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jj,jk) - zws(ji,jj,jk)
197               END_3D
198            ENDIF
199            !
200            !! Matrix inversion from the first level
201            !!----------------------------------------------------------------------
202            !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
203            !
204            !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
205            !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
206            !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
207            !        (        ...               )( ...  ) ( ...  )
208            !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
209            !
210            !   m is decomposed in the product of an upper and lower triangular matrix.
211            !   The 3 diagonal terms are in 3d arrays: zwd, zws, zwi.
212            !   Suffices i,s and d indicate "inferior" (below diagonal), diagonal
213            !   and "superior" (above diagonal) components of the tridiagonal system.
214            !   The solution will be in the 4d array pta.
215            !   The 3d array zwt is used as a work space array.
216            !   En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then
217            !   used as a work space array: its value is modified.
218            !
219            DO_2D( 0, 0, 0, 0 )      !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k) ! done one for all passive tracers (so included in the IF instruction)
220               zwt(ji,jj,1) = zwd(ji,jj,1)
221            END_2D
222            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
223               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
224            END_3D
225            !
226         ENDIF 
227         !         
228         DO_2D( 0, 0, 0, 0 )         !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
229            pt(ji,jj,1,jn,Kaa) =        e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb)    &
230               &               + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs)
231         END_2D
232         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
233            zrhs =        e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb)    &
234               & + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs)   ! zrhs=right hand side
235            pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa)
236         END_3D
237         !
238         DO_2D( 0, 0, 0, 0 )         !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk   (result is the after tracer)
239            pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
240         END_2D
241         DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 )
242            pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jj,jk) * pt(ji,jj,jk+1,jn,Kaa) )   &
243               &             / zwt(ji,jj,jk) * tmask(ji,jj,jk)
244         END_3D
245         !                                            ! ================= !
246      END DO                                          !  end tracer loop  !
247      !                                               ! ================= !
248   END SUBROUTINE tra_zdf_imp
249
250   !!==============================================================================
251END MODULE trazdf
Note: See TracBrowser for help on using the repository browser.