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/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trazdf.F90 @ 10825

Last change on this file since 10825 was 10825, checked in by davestorkey, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : trazdf.F90 and dynzdf.F90
  • Property svn:keywords set to Id
File size: 14.6 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 "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
41   !! $Id$
42   !! Software governed by the CeCILL license (see ./LICENSE)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE tra_zdf( kt, ktlev1, ktlev2, ktlev3, kt2lev, pts_rhs )
47      !!----------------------------------------------------------------------
48      !!                  ***  ROUTINE tra_zdf  ***
49      !!
50      !! ** Purpose :   compute the vertical ocean tracer physics.
51      !!---------------------------------------------------------------------
52      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
53      INTEGER, INTENT(in) ::   ktlev1, ktlev2, ktlev3   ! time level indices for 3-time-level source terms
54      INTEGER, INTENT(in) ::   kt2lev                   ! time level index for 2-time-level source terms
55      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends
56      !
57      INTEGER  ::   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(lwp)WRITE(numout,*)
65         IF(lwp)WRITE(numout,*) 'tra_zdf : implicit vertical mixing on T & S'
66         IF(lwp)WRITE(numout,*) '~~~~~~~ '
67      ENDIF
68      !
69      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =      rdt   ! at nit000, =   rdt (restarting with Euler time stepping)
70      ELSEIF( kt <= nit000 + 1           ) THEN   ;   r2dt = 2. * rdt   ! otherwise, = 2 rdt (leapfrog)
71      ENDIF
72      !
73      IF( l_trdtra )   THEN                  !* Save ta and sa trends
74         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
75         ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem)
76         ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal)
77      ENDIF
78      !
79      !                                      !* compute lateral mixing trend and add it to the general trend
80      CALL tra_zdf_imp( kt, nit000, ktlev1, ktlev2, ktlev3, kt2lev, 'TRA', r2dt, ts(:,:,:,:,ktlev1), pts_rhs, jpts ) 
81
82!!gm WHY here !   and I don't like that !
83      ! DRAKKAR SSS control {
84      ! JMM avoid negative salinities near river outlet ! Ugly fix
85      ! JMM : restore negative salinities to small salinities:
86      WHERE( pts_rhs(:,:,:,jp_sal) < 0._wp )   pts_rhs(:,:,:,jp_sal) = 0.1_wp
87!!gm
88
89      IF( l_trdtra )   THEN                      ! save the vertical diffusive trends for further diagnostics
90         DO jk = 1, jpkm1
91            ztrdt(:,:,jk) = ( ( pts_rhs(:,:,jk,jp_tem)*e3t(:,:,jk,ktlev3) - ts(:,:,jk,jp_tem,ktlev1)*e3t(:,:,jk,ktlev1) ) &
92               &          / (e3t(:,:,jk,ktlev2)*r2dt) ) - ztrdt(:,:,jk)
93            ztrds(:,:,jk) = ( ( pts_rhs(:,:,jk,jp_sal)*e3t(:,:,jk,ktlev3) - ts(:,:,jk,jp_sal,ktlev1)*e3t(:,:,jk,ktlev1) ) &
94              &           / (e3t(:,:,jk,ktlev2)*r2dt) ) - ztrds(:,:,jk)
95         END DO
96!!gm this should be moved in trdtra.F90 and done on all trends
97         CALL lbc_lnk_multi( 'trazdf', ztrdt, 'T', 1. , ztrds, 'T', 1. )
98!!gm
99         CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt )
100         CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds )
101         DEALLOCATE( ztrdt , ztrds )
102      ENDIF
103      !                                          ! print mean trends (used for debugging)
104      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' zdf  - Ta: ', mask1=tmask,               &
105         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
106      !
107      IF( ln_timing )   CALL timing_stop('tra_zdf')
108      !
109   END SUBROUTINE tra_zdf
110
111 
112   SUBROUTINE tra_zdf_imp( kt, kit000, ktlev1, ktlev2, ktlev3, kt2lev, cdtype, p2dt, pt, pt_rhs, kjpt ) 
113      !!----------------------------------------------------------------------
114      !!                  ***  ROUTINE tra_zdf_imp  ***
115      !!
116      !! ** Purpose :   Compute the after tracer through a implicit computation
117      !!     of the vertical tracer diffusion (including the vertical component
118      !!     of lateral mixing (only for 2nd order operator, for fourth order
119      !!     it is already computed and add to the general trend in traldf)
120      !!
121      !! ** Method  :  The vertical diffusion of a tracer ,t , is given by:
122      !!          difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
123      !!      It is computed using a backward time scheme (t=after field)
124      !!      which provide directly the after tracer field.
125      !!      If ln_zdfddm=T, use avs for salinity or for passive tracers
126      !!      Surface and bottom boundary conditions: no diffusive flux on
127      !!      both tracers (bottom, applied through the masked field avt).
128      !!      If iso-neutral mixing, add to avt the contribution due to lateral mixing.
129      !!
130      !! ** Action  : - pt_rhs  becomes the after tracer
131      !!---------------------------------------------------------------------
132      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
133      INTEGER                              , INTENT(in   ) ::   kit000   ! first time step index
134      INTEGER                              , INTENT(in   ) ::   ktlev1, ktlev2, ktlev3   ! time level indices for 3-time-level source terms
135      INTEGER                              , INTENT(in   ) ::   kt2lev                   ! time level index for 2-time-level source terms
136      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
137      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
138      REAL(wp)                             , INTENT(in   ) ::   p2dt     ! tracer time-step
139      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt      ! before and now tracer fields
140      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs      ! in: tracer trend ; out: after tracer field
141      !
142      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
143      REAL(wp) ::  zrhs, zzwi, zzws ! local scalars
144      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zwi, zwt, zwd, zws
145      !!---------------------------------------------------------------------
146      !
147      !                                               ! ============= !
148      DO jn = 1, kjpt                                 !  tracer loop  !
149         !                                            ! ============= !
150         !  Matrix construction
151         ! --------------------
152         ! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer
153         !
154         IF(  ( cdtype == 'TRA' .AND. ( jn == jp_tem .OR. ( jn == jp_sal .AND. ln_zdfddm ) ) ) .OR.   &
155            & ( cdtype == 'TRC' .AND. jn == 1 )  )  THEN
156            !
157            ! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers
158            IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN   ;   zwt(:,:,2:jpk) = avt(:,:,2:jpk)
159            ELSE                                            ;   zwt(:,:,2:jpk) = avs(:,:,2:jpk)
160            ENDIF
161            zwt(:,:,1) = 0._wp
162            !
163            IF( l_ldfslp ) THEN            ! isoneutral diffusion: add the contribution
164               IF( ln_traldf_msc  ) THEN     ! MSC iso-neutral operator
165                  DO jk = 2, jpkm1
166                     DO jj = 2, jpjm1
167                        DO ji = fs_2, fs_jpim1   ! vector opt.
168                           zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk) 
169                        END DO
170                     END DO
171                  END DO
172               ELSE                          ! standard or triad iso-neutral operator
173                  DO jk = 2, jpkm1
174                     DO jj = 2, jpjm1
175                        DO ji = fs_2, fs_jpim1   ! vector opt.
176                           zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk)
177                        END DO
178                     END DO
179                  END DO
180               ENDIF
181            ENDIF
182            !
183            ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked)
184            IF( ln_zad_Aimp ) THEN         ! Adaptive implicit vertical advection
185               DO jk = 1, jpkm1
186                  DO jj = 2, jpjm1
187                     DO ji = fs_2, fs_jpim1   ! vector opt. (ensure same order of calculation as below if wi=0.)
188                        zzwi = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk  ,kt2lev)
189                        zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,kt2lev)
190                        zwd(ji,jj,jk) = e3t(ji,jj,jk,ktlev3) - zzwi - zzws   &
191                           &                 + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) )
192                        zwi(ji,jj,jk) = zzwi + p2dt *   MIN( wi(ji,jj,jk  ) , 0._wp )
193                        zws(ji,jj,jk) = zzws - p2dt *   MAX( wi(ji,jj,jk+1) , 0._wp )
194                    END DO
195                  END DO
196               END DO
197            ELSE
198               DO jk = 1, jpkm1
199                  DO jj = 2, jpjm1
200                     DO ji = fs_2, fs_jpim1   ! vector opt.
201                        zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk,kt2lev)
202                        zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,kt2lev)
203                        zwd(ji,jj,jk) = e3t(ji,jj,jk,ktlev3) - zwi(ji,jj,jk) - zws(ji,jj,jk)
204                    END DO
205                  END DO
206               END DO
207            ENDIF
208            !
209            !! Matrix inversion from the first level
210            !!----------------------------------------------------------------------
211            !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
212            !
213            !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
214            !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
215            !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
216            !        (        ...               )( ...  ) ( ...  )
217            !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
218            !
219            !   m is decomposed in the product of an upper and lower triangular matrix.
220            !   The 3 diagonal terms are in 3d arrays: zwd, zws, zwi.
221            !   Suffices i,s and d indicate "inferior" (below diagonal), diagonal
222            !   and "superior" (above diagonal) components of the tridiagonal system.
223            !   The solution will be in the 4d array pt_rhs.
224            !   The 3d array zwt is used as a work space array.
225            !   En route to the solution pt_rhs is used a to evaluate the rhs and then
226            !   used as a work space array: its value is modified.
227            !
228            DO jj = 2, jpjm1        !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
229               DO ji = fs_2, fs_jpim1            ! done one for all passive tracers (so included in the IF instruction)
230                  zwt(ji,jj,1) = zwd(ji,jj,1)
231               END DO
232            END DO
233            DO jk = 2, jpkm1
234               DO jj = 2, jpjm1
235                  DO ji = fs_2, fs_jpim1
236                     zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
237                  END DO
238               END DO
239            END DO
240            !
241         ENDIF 
242         !         
243         DO jj = 2, jpjm1           !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
244            DO ji = fs_2, fs_jpim1
245               pt_rhs(ji,jj,1,jn) = e3t(ji,jj,1,ktlev1) * pt(ji,jj,1,jn) + p2dt * e3t(ji,jj,1,ktlev2) * pt_rhs(ji,jj,1,jn)
246            END DO
247         END DO
248         DO jk = 2, jpkm1
249            DO jj = 2, jpjm1
250               DO ji = fs_2, fs_jpim1
251                  zrhs = e3t(ji,jj,jk,ktlev1) * pt(ji,jj,jk,jn) + p2dt * e3t(ji,jj,jk,ktlev2) * pt_rhs(ji,jj,jk,jn)   ! zrhs=right hand side
252                  pt_rhs(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt_rhs(ji,jj,jk-1,jn)
253               END DO
254            END DO
255         END DO
256         !
257         DO jj = 2, jpjm1           !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk   (result is the after tracer)
258            DO ji = fs_2, fs_jpim1
259               pt_rhs(ji,jj,jpkm1,jn) = pt_rhs(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
260            END DO
261         END DO
262         DO jk = jpk-2, 1, -1
263            DO jj = 2, jpjm1
264               DO ji = fs_2, fs_jpim1
265                  pt_rhs(ji,jj,jk,jn) = ( pt_rhs(ji,jj,jk,jn) - zws(ji,jj,jk) * pt_rhs(ji,jj,jk+1,jn) )   &
266                     &             / zwt(ji,jj,jk) * tmask(ji,jj,jk)
267               END DO
268            END DO
269         END DO
270         !                                            ! ================= !
271      END DO                                          !  end tracer loop  !
272      !                                               ! ================= !
273   END SUBROUTINE tra_zdf_imp
274
275   !!==============================================================================
276END MODULE trazdf
Note: See TracBrowser for help on using the repository browser.