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.
traldf_iso.F90 in NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA – NEMO

source: NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TRA/traldf_iso.F90 @ 12397

Last change on this file since 12397 was 12397, checked in by davestorkey, 4 years ago

2020/KERNEL-03_Storkey_Coward_RK3_stage2 : Consolidation of code to
handle initial Euler timestep in the context of leapfrog
timestepping. This version passes all SETTE tests but fails to bit
compare with the control for several tests (ORCA2_ICE_PISCES, AMM12,
ISOMIP, AGRIF_DEMO, SPITZ12).
  • Property svn:keywords set to Id
File size: 18.9 KB
Line 
1MODULE traldf_iso
2   !!======================================================================
3   !!                   ***  MODULE  traldf_iso  ***
4   !! Ocean  tracers:  horizontal component of the lateral tracer mixing trend
5   !!======================================================================
6   !! History :  OPA  ! 1994-08  (G. Madec, M. Imbard)
7   !!            8.0  ! 1997-05  (G. Madec)  split into traldf and trazdf
8   !!            NEMO ! 2002-08  (G. Madec)  Free form, F90
9   !!            1.0  ! 2005-11  (G. Madec)  merge traldf and trazdf :-)
10   !!            3.3  ! 2010-09  (C. Ethe, G. Madec) Merge TRA-TRC
11   !!            3.7  ! 2014-01  (G. Madec, S. Masson)  restructuration/simplification of aht/aeiv specification
12   !!             -   ! 2014-02  (F. Lemarie, G. Madec)  triad operator (Griffies) + Method of Stabilizing Correction
13   !!----------------------------------------------------------------------
14
15   !!----------------------------------------------------------------------
16   !!   tra_ldf_iso   : update the tracer trend with the horizontal component of a iso-neutral laplacian operator
17   !!                   and with the vertical part of the isopycnal or geopotential s-coord. operator
18   !!----------------------------------------------------------------------
19   USE oce            ! ocean dynamics and active tracers
20   USE dom_oce        ! ocean space and time domain
21   USE trc_oce        ! share passive tracers/Ocean variables
22   USE zdf_oce        ! ocean vertical physics
23   USE ldftra         ! lateral diffusion: tracer eddy coefficients
24   USE ldfslp         ! iso-neutral slopes
25   USE diaptr         ! poleward transport diagnostics
26   USE diaar5         ! AR5 diagnostics
27   !
28   USE in_out_manager ! I/O manager
29   USE iom            ! I/O library
30   USE phycst         ! physical constants
31   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   tra_ldf_iso   ! routine called by step.F90
37
38   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
39   LOGICAL  ::   l_hst   ! flag to compute heat transport
40
41   !! * Substitutions
42#  include "do_loop_substitute.h90"
43   !!----------------------------------------------------------------------
44   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
45   !! $Id$
46   !! Software governed by the CeCILL license (see ./LICENSE)
47   !!----------------------------------------------------------------------
48CONTAINS
49
50  SUBROUTINE tra_ldf_iso( kt, Kmm, kit000, cdtype, pahu, pahv,                    &
51      &                                            pgu , pgv    ,   pgui, pgvi,   &
52      &                                       pt , pt2 , pt_rhs , kjpt  , kpass )
53      !!----------------------------------------------------------------------
54      !!                  ***  ROUTINE tra_ldf_iso  ***
55      !!
56      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
57      !!      trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and
58      !!      add it to the general trend of tracer equation.
59      !!
60      !! ** Method  :   The horizontal component of the lateral diffusive trends
61      !!      is provided by a 2nd order operator rotated along neural or geopo-
62      !!      tential surfaces to which an eddy induced advection can be added
63      !!      It is computed using before fields (forward in time) and isopyc-
64      !!      nal or geopotential slopes computed in routine ldfslp.
65      !!
66      !!      1st part :  masked horizontal derivative of T  ( di[ t ] )
67      !!      ========    with partial cell update if ln_zps=T
68      !!                  with top     cell update if ln_isfcav
69      !!
70      !!      2nd part :  horizontal fluxes of the lateral mixing operator
71      !!      ========   
72      !!         zftu =  pahu e2u*e3u/e1u di[ tb ]
73      !!               - pahu e2u*uslp    dk[ mi(mk(tb)) ]
74      !!         zftv =  pahv e1v*e3v/e2v dj[ tb ]
75      !!               - pahv e2u*vslp    dk[ mj(mk(tb)) ]
76      !!      take the horizontal divergence of the fluxes:
77      !!         difft = 1/(e1e2t*e3t) {  di-1[ zftu ] +  dj-1[ zftv ]  }
78      !!      Add this trend to the general trend (ta,sa):
79      !!         ta = ta + difft
80      !!
81      !!      3rd part: vertical trends of the lateral mixing operator
82      !!      ========  (excluding the vertical flux proportional to dk[t] )
83      !!      vertical fluxes associated with the rotated lateral mixing:
84      !!         zftw = - {  mi(mk(pahu)) * e2t*wslpi di[ mi(mk(tb)) ]
85      !!                   + mj(mk(pahv)) * e1t*wslpj dj[ mj(mk(tb)) ]  }
86      !!      take the horizontal divergence of the fluxes:
87      !!         difft = 1/(e1e2t*e3t) dk[ zftw ]
88      !!      Add this trend to the general trend (ta,sa):
89      !!         pt_rhs = pt_rhs + difft
90      !!
91      !! ** Action :   Update pt_rhs arrays with the before rotated diffusion
92      !!----------------------------------------------------------------------
93      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
94      INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index
95      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
96      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
97      INTEGER                              , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
98      INTEGER                              , INTENT(in   ) ::   Kmm        ! ocean time level index
99      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
100      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
101      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
102      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt         ! tracer (kpass=1) or laplacian of tracer (kpass=2)
103      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt2        ! tracer (only used in kpass=2)
104      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs     ! tracer trend
105      !
106      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
107      INTEGER  ::  ikt
108      INTEGER  ::  ierr             ! local integer
109      REAL(wp) ::  zmsku, zahu_w, zabe1, zcof1, zcoef3   ! local scalars
110      REAL(wp) ::  zmskv, zahv_w, zabe2, zcof2, zcoef4   !   -      -
111      REAL(wp) ::  zcoef0, ze3w_2, zsign                 !   -      -
112      REAL(wp), DIMENSION(jpi,jpj)     ::   zdkt, zdk1t, z2d
113      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdit, zdjt, zftu, zftv, ztfw 
114      !!----------------------------------------------------------------------
115      !
116      IF( kpass == 1 .AND. kt == kit000 )  THEN
117         IF(lwp) WRITE(numout,*)
118         IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator on ', cdtype
119         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
120         !
121         akz     (:,:,:) = 0._wp     
122         ah_wslp2(:,:,:) = 0._wp
123      ENDIF
124      !   
125      l_hst = .FALSE.
126      l_ptr = .FALSE.
127      IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) )     l_ptr = .TRUE. 
128      IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
129         &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE.
130      !
131      !
132      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
133      ELSE                    ;   zsign = -1._wp
134      ENDIF
135         
136      !!----------------------------------------------------------------------
137      !!   0 - calculate  ah_wslp2 and akz
138      !!----------------------------------------------------------------------
139      !
140      IF( kpass == 1 ) THEN                  !==  first pass only  ==!
141         !
142         DO_3D_00_00( 2, jpkm1 )
143            !
144            zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          &
145               &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  )
146            zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          &
147               &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  )
148               !
149            zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    &
150               &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku
151            zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    &
152               &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv
153               !
154            ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
155               &               + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk)
156         END_3D
157         !
158         IF( ln_traldf_msc ) THEN                ! stabilizing vertical diffusivity coefficient
159            DO_3D_00_00( 2, jpkm1 )
160               akz(ji,jj,jk) = 0.25_wp * (                                                                     &
161                  &              ( pahu(ji  ,jj,jk) + pahu(ji  ,jj,jk-1) ) / ( e1u(ji  ,jj) * e1u(ji  ,jj) )   &
162                  &            + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) )   &
163                  &            + ( pahv(ji,jj  ,jk) + pahv(ji,jj  ,jk-1) ) / ( e2v(ji,jj  ) * e2v(ji,jj  ) )   &
164                  &            + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) )   )
165            END_3D
166            !
167            IF( ln_traldf_blp ) THEN                ! bilaplacian operator
168               DO_3D_10_10( 2, jpkm1 )
169                  akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   &
170                     &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) )  )
171               END_3D
172            ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator
173               DO_3D_10_10( 2, jpkm1 )
174                  ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)
175                  zcoef0 = r2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
176                  akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * r1_2dt
177               END_3D
178           ENDIF
179           !
180         ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit
181            akz(:,:,:) = ah_wslp2(:,:,:)     
182         ENDIF
183      ENDIF
184      !
185      !                                                          ! ===========
186      DO jn = 1, kjpt                                            ! tracer loop
187         !                                                       ! ===========
188         !                                               
189         !!----------------------------------------------------------------------
190         !!   I - masked horizontal derivative
191         !!----------------------------------------------------------------------
192!!gm : bug.... why (x,:,:)?   (1,jpj,:) and (jpi,1,:) should be sufficient....
193         zdit (1,:,:) = 0._wp     ;     zdit (jpi,:,:) = 0._wp
194         zdjt (1,:,:) = 0._wp     ;     zdjt (jpi,:,:) = 0._wp
195         !!end
196
197         ! Horizontal tracer gradient
198         DO_3D_10_10( 1, jpkm1 )
199            zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)
200            zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
201         END_3D
202         IF( ln_zps ) THEN      ! botton and surface ocean correction of the horizontal gradient
203            DO_2D_10_10
204               zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)         
205               zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
206            END_2D
207            IF( ln_isfcav ) THEN      ! first wet level beneath a cavity
208               DO_2D_10_10
209                  IF( miku(ji,jj) > 1 )   zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)         
210                  IF( mikv(ji,jj) > 1 )   zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)     
211               END_2D
212            ENDIF
213         ENDIF
214         !
215         !!----------------------------------------------------------------------
216         !!   II - horizontal trend  (full)
217         !!----------------------------------------------------------------------
218         !
219         DO jk = 1, jpkm1                                 ! Horizontal slab
220            !
221            !                             !== Vertical tracer gradient
222            zdk1t(:,:) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * wmask(:,:,jk+1)     ! level jk+1
223            !
224            IF( jk == 1 ) THEN   ;   zdkt(:,:) = zdk1t(:,:)                          ! surface: zdkt(jk=1)=zdkt(jk=2)
225            ELSE                 ;   zdkt(:,:) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * wmask(:,:,jk)
226            ENDIF
227            DO_2D_10_10
228               zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm)
229               zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
230               !
231               zmsku = 1. / MAX(  wmask(ji+1,jj,jk  ) + wmask(ji,jj,jk+1)   &
232                  &             + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk  ), 1. )
233               !
234               zmskv = 1. / MAX(  wmask(ji,jj+1,jk  ) + wmask(ji,jj,jk+1)   &
235                  &             + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk  ), 1. )
236               !
237               zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku
238               zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv
239               !
240               zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)   &
241                  &               + zcof1 * (  zdkt (ji+1,jj) + zdk1t(ji,jj)      &
242                  &                          + zdk1t(ji+1,jj) + zdkt (ji,jj)  )  ) * umask(ji,jj,jk)
243               zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   &
244                  &               + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj)      &
245                  &                          + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk)                 
246            END_2D
247            !
248            DO_2D_00_00
249               pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      &
250                  &                                                 + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   &
251                  &                                              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
252            END_2D
253         END DO                                        !   End of slab 
254
255         !!----------------------------------------------------------------------
256         !!   III - vertical trend (full)
257         !!----------------------------------------------------------------------
258         !
259         ! Vertical fluxes
260         ! ---------------
261         !                          ! Surface and bottom vertical fluxes set to zero
262         ztfw(:,:, 1 ) = 0._wp      ;      ztfw(:,:,jpk) = 0._wp
263         
264         DO_3D_00_00( 2, jpkm1 )
265            !
266            zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          &
267               &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  )
268            zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          &
269               &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  )
270               !
271            zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    &
272               &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku
273            zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    &
274               &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv
275               !
276            zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk)   !wslpi & j are already w-masked
277            zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk)
278            !
279            ztfw(ji,jj,jk) = zcoef3 * (   zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)      &
280               &                        + zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk)  )   &
281               &           + zcoef4 * (   zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk)      &
282               &                        + zdjt(ji  ,jj-1,jk-1) + zdjt(ji  ,jj  ,jk)  )
283         END_3D
284         !                                !==  add the vertical 33 flux  ==!
285         IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz
286            DO_3D_00_00( 2, jpkm1 )
287               ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)   &
288                  &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )               &
289                  &                            * (  pt(ji,jj,jk-1,jn) -  pt(ji,jj,jk,jn) )
290            END_3D
291            !
292         ELSE                                   ! bilaplacian
293            SELECT CASE( kpass )
294            CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2
295               DO_3D_00_00( 2, jpkm1 )
296                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk)                       &
297                     &           + ah_wslp2(ji,jj,jk)  * e1e2t(ji,jj)   &
298                     &           * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)
299               END_3D
300            CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt2 gradients, resp.
301               DO_3D_00_00( 2, jpkm1 )
302                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)                  &
303                     &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   &
304                     &                            +         akz(ji,jj,jk) * ( pt2(ji,jj,jk-1,jn) - pt2(ji,jj,jk,jn) )   )
305               END_3D
306            END SELECT
307         ENDIF
308         !         
309         DO_3D_00_00( 1, jpkm1 )
310            pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   &
311               &                                              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
312         END_3D
313         !
314         IF( ( kpass == 1 .AND. ln_traldf_lap ) .OR.  &     !==  first pass only (  laplacian)  ==!
315             ( kpass == 2 .AND. ln_traldf_blp ) ) THEN      !==  2nd   pass      (bilaplacian)  ==!
316            !
317            !                             ! "Poleward" diffusive heat or salt transports (T-S case only)
318               ! note sign is reversed to give down-gradient diffusive transports )
319            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', -zftv(:,:,:)  )
320            !                          ! Diffusive heat transports
321            IF( l_hst )  CALL dia_ar5_hst( jn, 'ldf', -zftu(:,:,:), -zftv(:,:,:) )
322            !
323         ENDIF                                                    !== end pass selection  ==!
324         !
325         !                                                        ! ===============
326      END DO                                                      ! end tracer loop
327      !
328   END SUBROUTINE tra_ldf_iso
329
330   !!==============================================================================
331END MODULE traldf_iso
Note: See TracBrowser for help on using the repository browser.