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_lap_blp.F90 in NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/TRA – NEMO

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/TRA/traldf_lap_blp.F90 @ 14644

Last change on this file since 14644 was 14644, checked in by sparonuz, 3 years ago

Merge trunk -r14642:HEAD
  • Property svn:keywords set to Id
File size: 16.3 KB
Line 
1MODULE traldf_lap_blp
2   !!==============================================================================
3   !!                       ***  MODULE  traldf_lap_blp  ***
4   !! Ocean tracers:  lateral diffusivity trend  (laplacian and bilaplacian)
5   !!==============================================================================
6   !! History :  3.7  ! 2014-01  (G. Madec, S. Masson)  Original code, re-entrant laplacian
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   tra_ldf_lap   : tracer trend update with iso-level laplacian diffusive operator
11   !!   tra_ldf_blp   : tracer trend update with iso-level or iso-neutral bilaplacian operator
12   !!----------------------------------------------------------------------
13   USE oce            ! ocean dynamics and active tracers
14   USE dom_oce        ! ocean space and time domain
15   USE domutl, ONLY : is_tile
16   USE ldftra         ! lateral physics: eddy diffusivity
17   USE traldf_iso     ! iso-neutral lateral diffusion (standard operator)     (tra_ldf_iso   routine)
18   USE traldf_triad   ! iso-neutral lateral diffusion (triad    operator)     (tra_ldf_triad routine)
19   USE diaptr         ! poleward transport diagnostics
20   USE diaar5         ! AR5 diagnostics
21   USE trc_oce        ! share passive tracers/Ocean variables
22   USE zpshde         ! partial step: hor. derivative     (zps_hde routine)
23   !
24   USE in_out_manager ! I/O manager
25   USE iom            ! I/O library
26   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
27   USE lib_mpp        ! distribued memory computing library
28   USE timing         ! Timing
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   tra_ldf_lap   ! called by traldf.F90
34   PUBLIC   tra_ldf_blp   ! called by traldf.F90
35
36   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
37   LOGICAL  ::   l_hst   ! flag to compute heat transport
38
39   !! * Substitutions
40#  include "do_loop_substitute.h90"
41#  include "domzgr_substitute.h90"
42#  include "single_precision_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_lap( kt, Kmm, kit000, cdtype, pahu, pahv,             &
51      &                                             pgu , pgv , pgui, pgvi, &
52      &                                             pt, pt_rhs, kjpt, kpass )
53      !!
54      INTEGER                     , INTENT(in   ) ::   kt         ! ocean time-step index
55      INTEGER                     , INTENT(in   ) ::   kit000     ! first time step index
56      CHARACTER(len=3)            , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
57      INTEGER                     , INTENT(in   ) ::   kjpt       ! number of tracers
58      INTEGER                     , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
59      INTEGER                     , INTENT(in   ) ::   Kmm        ! ocean time level index
60      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
61      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
62      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
63      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pt         ! before tracer fields
64      REAL(dp), DIMENSION(:,:,:,:), INTENT(inout) ::   pt_rhs     ! tracer trend
65      !!
66      CALL tra_ldf_lap_t( kt, Kmm, kit000, cdtype, pahu, pahv, is_tile(pahu),                            &
67      &                                            pgu , pgv , is_tile(pgu) , pgui, pgvi, is_tile(pgui), &
68      &                                            pt, is_tile(pt), pt_rhs, is_tile(pt_rhs), kjpt, kpass )
69   END SUBROUTINE tra_ldf_lap
70
71
72   SUBROUTINE tra_ldf_lap_t( kt, Kmm, kit000, cdtype, pahu, pahv, ktah,                   &
73      &                                               pgu , pgv , ktg , pgui, pgvi, ktgi, &
74      &                                               pt, ktt, pt_rhs, ktt_rhs, kjpt, kpass )
75      !!----------------------------------------------------------------------
76      !!                  ***  ROUTINE tra_ldf_lap  ***
77      !!
78      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
79      !!      trend and add it to the general trend of tracer equation.
80      !!
81      !! ** Method  :   Second order diffusive operator evaluated using before
82      !!      fields (forward time scheme). The horizontal diffusive trends of
83      !!      the tracer is given by:
84      !!          difft = 1/(e1e2t*e3t) {  di-1[ pahu e2u*e3u/e1u di(tb) ]
85      !!                                 + dj-1[ pahv e1v*e3v/e2v dj(tb) ] }
86      !!      Add this trend to the general tracer trend pt_rhs :
87      !!          pt_rhs = pt_rhs + difft
88      !!
89      !! ** Action  : - Update pt_rhs arrays with the before iso-level
90      !!                harmonic mixing trend.
91      !!----------------------------------------------------------------------
92      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
93      INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index
94      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
95      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
96      INTEGER                              , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
97      INTEGER                              , INTENT(in   ) ::   Kmm        ! ocean time level index
98      INTEGER                              , INTENT(in   ) ::   ktah, ktg, ktgi, ktt, ktt_rhs
99      REAL(wp), DIMENSION(A2D_T(ktah),   JPK)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
100      REAL(wp), DIMENSION(A2D_T(ktg),        KJPT), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
101      REAL(wp), DIMENSION(A2D_T(ktgi),       KJPT), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
102      REAL(wp), DIMENSION(A2D_T(ktt),    JPK,KJPT), INTENT(in   ) ::   pt         ! before tracer fields
103      REAL(dp), DIMENSION(A2D_T(ktt_rhs),JPK,KJPT), INTENT(inout) ::   pt_rhs     ! tracer trend
104      !
105      INTEGER  ::   ji, jj, jk, jn      ! dummy loop indices
106      INTEGER  ::   isi, iei, isj, iej  ! local integers
107      REAL(wp) ::   zsign               ! local scalars
108      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   ztu, ztv, zaheeu, zaheev
109      !!----------------------------------------------------------------------
110      !
111      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
112         IF( kt == nit000 .AND. lwp )  THEN
113            WRITE(numout,*)
114            WRITE(numout,*) 'tra_ldf_lap : iso-level laplacian diffusion on ', cdtype, ', pass=', kpass
115            WRITE(numout,*) '~~~~~~~~~~~ '
116         ENDIF
117         !
118         l_hst = .FALSE.
119         l_ptr = .FALSE.
120         IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) )     l_ptr = .TRUE.
121         IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
122            &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE.
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      !                                !==  Initialization of metric arrays used for all tracers  ==!
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      IF( ntsi == Nis0 ) THEN ; isi = nn_hls - 1 ; ELSE ; isi = 0 ; ENDIF    ! Avoid double-counting when using tiling
137      IF( ntsj == Njs0 ) THEN ; isj = nn_hls - 1 ; ELSE ; isj = 0 ; ENDIF
138      IF( ntei == Nie0 ) THEN ; iei = nn_hls - 1 ; ELSE ; iei = 0 ; ENDIF
139      IF( ntej == Nje0 ) THEN ; iej = nn_hls - 1 ; ELSE ; iej = 0 ; ENDIF
140
141      DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 )            !== First derivative (gradient)  ==!
142         zaheeu(ji,jj,jk) = zsign * pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm)   !!gm   * umask(ji,jj,jk) pah masked!
143         zaheev(ji,jj,jk) = zsign * pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm)   !!gm   * vmask(ji,jj,jk)
144      END_3D
145      !
146      !                             ! =========== !
147      DO jn = 1, kjpt               ! tracer loop !
148         !                          ! =========== !
149         !
150         DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 )            !== First derivative (gradient)  ==!
151            ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) )
152            ztv(ji,jj,jk) = zaheev(ji,jj,jk) * ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) )
153         END_3D
154         IF( ln_zps ) THEN                             ! set gradient at bottom/top ocean level
155            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )                              ! bottom
156               ztu(ji,jj,mbku(ji,jj)) = zaheeu(ji,jj,mbku(ji,jj)) * pgu(ji,jj,jn)
157               ztv(ji,jj,mbkv(ji,jj)) = zaheev(ji,jj,mbkv(ji,jj)) * pgv(ji,jj,jn)
158            END_2D
159            IF( ln_isfcav ) THEN                             ! top in ocean cavities only
160               DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
161                  IF( miku(ji,jj) > 1 )   ztu(ji,jj,miku(ji,jj)) = zaheeu(ji,jj,miku(ji,jj)) * pgui(ji,jj,jn)
162                  IF( mikv(ji,jj) > 1 )   ztv(ji,jj,mikv(ji,jj)) = zaheev(ji,jj,mikv(ji,jj)) * pgvi(ji,jj,jn)
163               END_2D
164            ENDIF
165         ENDIF
166         !
167         DO_3D( isi, iei, isj, iej, 1, jpkm1 )            !== Second derivative (divergence) added to the general tracer trends  ==!
168            pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)     &
169               &                                      +    ztv(ji,jj,jk) - ztv(ji,jj-1,jk) )   &
170               &                                      / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) )
171         END_3D
172         !
173         !                             !== "Poleward" diffusive heat or salt transports  ==!
174         IF( ( kpass == 1 .AND. .NOT.ln_traldf_blp ) .OR.  &     !==  first pass only (  laplacian)  ==!
175             ( kpass == 2 .AND.      ln_traldf_blp ) ) THEN      !==  2nd   pass only (bilaplacian)  ==!
176
177            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', -ztv(:,:,:)  )
178            IF( l_hst )  CALL dia_ar5_hst( jn, 'ldf', -ztu(:,:,:), -ztv(:,:,:) )
179         ENDIF
180         !                          ! ==================
181      END DO                        ! end of tracer loop
182      !                             ! ==================
183      !
184   END SUBROUTINE tra_ldf_lap_t
185
186
187   SUBROUTINE tra_ldf_blp( kt, Kmm, kit000, cdtype, pahu, pahv  ,             &
188      &                                             pgu , pgv   , pgui, pgvi, &
189      &                                             pt  , pt_rhs, kjpt, kldf )
190      !!----------------------------------------------------------------------
191      !!                 ***  ROUTINE tra_ldf_blp  ***
192      !!
193      !! ** Purpose :   Compute the before lateral tracer diffusive
194      !!      trend and add it to the general trend of tracer equation.
195      !!
196      !! ** Method  :   The lateral diffusive trends is provided by a bilaplacian
197      !!      operator applied to before field (forward in time).
198      !!      It is computed by two successive calls to laplacian routine
199      !!
200      !! ** Action :   pta   updated with the before rotated bilaplacian diffusion
201      !!----------------------------------------------------------------------
202      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
203      INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index
204      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
205      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
206      INTEGER                              , INTENT(in   ) ::   kldf       ! type of operator used
207      INTEGER                              , INTENT(in   ) ::   Kmm        ! ocean time level indices
208      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
209      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
210      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top levels
211      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt         ! before and now tracer fields
212      REAL(dp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs     ! tracer trend
213      !
214      INTEGER ::   ji, jj, jk, jn   ! dummy loop indices
215      REAL(dp), DIMENSION(A2D(nn_hls),jpk,kjpt) :: zlap         ! laplacian at t-point
216      REAL(wp), DIMENSION(A2D(nn_hls),    kjpt) :: zglu, zglv   ! bottom GRADh of the laplacian (u- and v-points)
217      REAL(wp), DIMENSION(A2D(nn_hls),    kjpt) :: zgui, zgvi   ! top    GRADh of the laplacian (u- and v-points)
218      !!---------------------------------------------------------------------
219      !
220      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
221         IF( kt == kit000 .AND. lwp )  THEN
222            WRITE(numout,*)
223            SELECT CASE ( kldf )
224            CASE ( np_blp    )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-level   bilaplacian operator on ', cdtype
225            CASE ( np_blp_i  )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (Standard)'
226            CASE ( np_blp_it )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (triad)'
227            END SELECT
228            WRITE(numout,*) '~~~~~~~~~~~'
229         ENDIF
230      ENDIF
231
232      zlap(:,:,:,:) = 0._wp
233      !
234      SELECT CASE ( kldf )       !==  1st laplacian applied to pt (output in zlap)  ==!
235      !
236      CASE ( np_blp    )               ! iso-level bilaplacian
237         CALL tra_ldf_lap  ( kt, Kmm, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, pt,     zlap, kjpt, 1 )
238      CASE ( np_blp_i  )               ! rotated   bilaplacian : standard operator (Madec)
239         CALL tra_ldf_iso  ( kt, Kmm, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, pt, pt, zlap, kjpt, 1 )
240      CASE ( np_blp_it )               ! rotated  bilaplacian : triad operator (griffies)
241         CALL tra_ldf_triad( kt, Kmm, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, pt, pt, zlap, kjpt, 1 )
242      END SELECT
243      !
244      CALL lbc_lnk( 'traldf_lap_blp', zlap(:,:,:,:) , 'T', 1.0_wp )     ! Lateral boundary conditions (unchanged sign)
245      !                                               ! Partial top/bottom cell: GRADh( zlap )
246      IF( ln_isfcav .AND. ln_zps ) THEN   ;   CALL zps_hde_isf( kt, Kmm, kjpt, zlap, zglu, zglv, zgui, zgvi )  ! both top & bottom
247      ELSEIF(             ln_zps ) THEN   ;   CALL zps_hde    ( kt, Kmm, kjpt, zlap, zglu, zglv )              ! only bottom
248      ENDIF
249      !
250      SELECT CASE ( kldf )       !==  2nd laplacian applied to zlap (output in pt_rhs)  ==!
251      !
252      CASE ( np_blp    )               ! iso-level bilaplacian
253         CALL tra_ldf_lap  ( kt, Kmm, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, CASTWP(zlap), pt_rhs,         kjpt, 2 )
254      CASE ( np_blp_i  )               ! rotated   bilaplacian : standard operator (Madec)
255         CALL tra_ldf_iso  ( kt, Kmm, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, CASTWP(zlap), pt    , pt_rhs, kjpt, 2 )
256      CASE ( np_blp_it )               ! rotated  bilaplacian : triad operator (griffies)
257         CALL tra_ldf_triad( kt, Kmm, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, CASTWP(zlap), pt    , pt_rhs, kjpt, 2 )
258      END SELECT
259      !
260   END SUBROUTINE tra_ldf_blp
261
262   !!==============================================================================
263END MODULE traldf_lap_blp
Note: See TracBrowser for help on using the repository browser.