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/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA – NEMO

source: NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/trazdf.F90 @ 13517

Last change on this file since 13517 was 13517, checked in by hadcv, 4 years ago

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