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 branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap_blp.F90 @ 9124

Last change on this file since 9124 was 9124, checked in by gm, 6 years ago

dev_merge_2017: ln_timing instead of nn_timing + restricted timing to nemo_init and routine called by step in OPA_SRC

  • Property svn:keywords set to Id
File size: 14.1 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 ldftra         ! lateral physics: eddy diffusivity
16   USE traldf_iso     ! iso-neutral lateral diffusion (standard operator)     (tra_ldf_iso   routine)
17   USE traldf_triad   ! iso-neutral lateral diffusion (triad    operator)     (tra_ldf_triad routine)
18   USE diaptr         ! poleward transport diagnostics
19   USE diaar5         ! AR5 diagnostics
20   USE trc_oce        ! share passive tracers/Ocean variables
21   USE zpshde         ! partial step: hor. derivative     (zps_hde routine)
22   !
23   USE in_out_manager ! I/O manager
24   USE iom            ! I/O library
25   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
26   USE lib_mpp        ! distribued memory computing library
27   USE timing         ! Timing
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   tra_ldf_lap   ! called by traldf.F90
33   PUBLIC   tra_ldf_blp   ! called by traldf.F90
34
35   !                      ! Flag to control the type of lateral diffusive operator
36   INTEGER, PARAMETER, PUBLIC ::   np_ERROR  =-10   ! error in specification of lateral diffusion
37   INTEGER, PARAMETER, PUBLIC ::   np_no_ldf = 00   ! without operator (i.e. no lateral diffusive trend)
38   !                          !!      laplacian     !    bilaplacian    !
39   INTEGER, PARAMETER, PUBLIC ::   np_lap    = 10   ,   np_blp    = 20  ! iso-level operator
40   INTEGER, PARAMETER, PUBLIC ::   np_lap_i  = 11   ,   np_blp_i  = 21  ! standard iso-neutral or geopotential operator
41   INTEGER, PARAMETER, PUBLIC ::   np_lap_it = 12   ,   np_blp_it = 22  ! triad    iso-neutral or geopotential operator
42
43   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
44   LOGICAL  ::   l_hst   ! flag to compute heat transport
45
46   !! * Substitutions
47#  include "vectopt_loop_substitute.h90"
48   !!----------------------------------------------------------------------
49   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
50   !! $Id$
51   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
52   !!----------------------------------------------------------------------
53CONTAINS
54
55   SUBROUTINE tra_ldf_lap( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   &
56      &                                                    pgui, pgvi,   &
57      &                                        ptb , pta , kjpt, kpass ) 
58      !!----------------------------------------------------------------------
59      !!                  ***  ROUTINE tra_ldf_lap  ***
60      !!                   
61      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
62      !!      trend and add it to the general trend of tracer equation.
63      !!
64      !! ** Method  :   Second order diffusive operator evaluated using before
65      !!      fields (forward time scheme). The horizontal diffusive trends of
66      !!      the tracer is given by:
67      !!          difft = 1/(e1e2t*e3t) {  di-1[ pahu e2u*e3u/e1u di(tb) ]
68      !!                                 + dj-1[ pahv e1v*e3v/e2v dj(tb) ] }
69      !!      Add this trend to the general tracer trend pta :
70      !!          pta = pta + difft
71      !!
72      !! ** Action  : - Update pta arrays with the before iso-level
73      !!                harmonic mixing trend.
74      !!----------------------------------------------------------------------
75      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
76      INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index
77      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
78      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
79      INTEGER                              , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
80      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
81      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
82      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
83      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields
84      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
85      !
86      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
87      REAL(wp) ::   zsign            ! local scalars
88      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztu, ztv, zaheeu, zaheev
89      !!----------------------------------------------------------------------
90      !
91      IF( kt == nit000 .AND. lwp )  THEN
92         WRITE(numout,*)
93         WRITE(numout,*) 'tra_ldf_lap : iso-level laplacian diffusion on ', cdtype, ', pass=', kpass
94         WRITE(numout,*) '~~~~~~~~~~~ '
95      ENDIF
96      !
97      l_hst = .FALSE.
98      l_ptr = .FALSE.
99      IF( cdtype == 'TRA' .AND. ln_diaptr )                                                l_ptr = .TRUE. 
100      IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
101         &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE.
102      !
103      !                                !==  Initialization of metric arrays used for all tracers  ==!
104      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
105      ELSE                    ;   zsign = -1._wp
106      ENDIF
107      DO jk = 1, jpkm1
108         DO jj = 1, jpjm1
109            DO ji = 1, fs_jpim1   ! vector opt.
110               zaheeu(ji,jj,jk) = zsign * pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk)   !!gm   * umask(ji,jj,jk) pah masked!
111               zaheev(ji,jj,jk) = zsign * pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk)   !!gm   * vmask(ji,jj,jk)
112            END DO
113         END DO
114      END DO
115      !
116      !                             ! =========== !
117      DO jn = 1, kjpt               ! tracer loop !
118         !                          ! =========== !   
119         !                               
120         DO jk = 1, jpkm1              !== First derivative (gradient)  ==!
121            DO jj = 1, jpjm1
122               DO ji = 1, fs_jpim1
123                  ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) )
124                  ztv(ji,jj,jk) = zaheev(ji,jj,jk) * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
125               END DO
126            END DO
127         END DO 
128         IF( ln_zps ) THEN                ! set gradient at bottom/top ocean level
129            DO jj = 1, jpjm1                    ! bottom
130               DO ji = 1, fs_jpim1
131                  ztu(ji,jj,mbku(ji,jj)) = zaheeu(ji,jj,mbku(ji,jj)) * pgu(ji,jj,jn)
132                  ztv(ji,jj,mbkv(ji,jj)) = zaheev(ji,jj,mbkv(ji,jj)) * pgv(ji,jj,jn)
133               END DO
134            END DO 
135            IF( ln_isfcav ) THEN                ! top in ocean cavities only
136               DO jj = 1, jpjm1
137                  DO ji = 1, fs_jpim1   ! vector opt.
138                     IF( miku(ji,jj) > 1 )   ztu(ji,jj,miku(ji,jj)) = zaheeu(ji,jj,miku(ji,jj)) * pgui(ji,jj,jn) 
139                     IF( mikv(ji,jj) > 1 )   ztv(ji,jj,mikv(ji,jj)) = zaheev(ji,jj,mikv(ji,jj)) * pgvi(ji,jj,jn) 
140                  END DO
141               END DO
142            ENDIF
143         ENDIF
144         !
145         DO jk = 1, jpkm1              !== Second derivative (divergence) added to the general tracer trends  ==!
146            DO jj = 2, jpjm1
147               DO ji = fs_2, fs_jpim1
148                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)     &
149                     &                                   + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) )   &
150                     &                                / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) )
151               END DO
152            END DO
153         END DO 
154         !
155         !                             !== "Poleward" diffusive heat or salt transports  ==!
156         IF( ( kpass == 1 .AND. .NOT.ln_traldf_blp ) .OR.  &     !==  first pass only (  laplacian)  ==!
157             ( kpass == 2 .AND.      ln_traldf_blp ) ) THEN      !==  2nd   pass only (bilaplacian)  ==!
158
159            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', -ztv(:,:,:)  )
160            IF( l_hst )  CALL dia_ar5_hst( jn, 'ldf', -ztu(:,:,:), -ztv(:,:,:) )
161         ENDIF
162         !                          ! ==================
163      END DO                        ! end of tracer loop
164      !                             ! ==================
165      !
166   END SUBROUTINE tra_ldf_lap
167   
168
169   SUBROUTINE tra_ldf_blp( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   &
170      &                                                    pgui, pgvi,   &
171      &                                                    ptb , pta , kjpt, kldf )
172      !!----------------------------------------------------------------------
173      !!                 ***  ROUTINE tra_ldf_blp  ***
174      !!                   
175      !! ** Purpose :   Compute the before lateral tracer diffusive
176      !!      trend and add it to the general trend of tracer equation.
177      !!
178      !! ** Method  :   The lateral diffusive trends is provided by a bilaplacian
179      !!      operator applied to before field (forward in time).
180      !!      It is computed by two successive calls to laplacian routine
181      !!
182      !! ** Action :   pta   updated with the before rotated bilaplacian diffusion
183      !!----------------------------------------------------------------------
184      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
185      INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index
186      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
187      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
188      INTEGER                              , INTENT(in   ) ::   kldf       ! type of operator used
189      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
190      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
191      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top levels
192      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields
193      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
194      !
195      INTEGER ::   ji, jj, jk, jn   ! dummy loop indices
196      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt) :: zlap         ! laplacian at t-point
197      REAL(wp), DIMENSION(jpi,jpj,    kjpt) :: zglu, zglv   ! bottom GRADh of the laplacian (u- and v-points)
198      REAL(wp), DIMENSION(jpi,jpj,    kjpt) :: zgui, zgvi   ! top    GRADh of the laplacian (u- and v-points)
199      !!---------------------------------------------------------------------
200      !
201      IF( kt == kit000 .AND. lwp )  THEN
202         WRITE(numout,*)
203         SELECT CASE ( kldf )
204         CASE ( np_blp    )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-level   bilaplacian operator on ', cdtype
205         CASE ( np_blp_i  )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (Standard)'
206         CASE ( np_blp_it )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (triad)'
207         END SELECT
208         WRITE(numout,*) '~~~~~~~~~~~'
209      ENDIF
210
211      zlap(:,:,:,:) = 0._wp
212      !
213      SELECT CASE ( kldf )       !==  1st laplacian applied to ptb (output in zlap)  ==!
214      !
215      CASE ( np_blp    )               ! iso-level bilaplacian
216         CALL tra_ldf_lap  ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb,      zlap, kjpt, 1 )
217      CASE ( np_blp_i  )               ! rotated   bilaplacian : standard operator (Madec)
218         CALL tra_ldf_iso  ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1 )
219      CASE ( np_blp_it )               ! rotated  bilaplacian : triad operator (griffies)
220         CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1 )
221      END SELECT
222      !
223      CALL lbc_lnk( zlap(:,:,:,:) , 'T', 1. )     ! Lateral boundary conditions (unchanged sign)
224      !                                               ! Partial top/bottom cell: GRADh( zlap ) 
225      IF( ln_isfcav .AND. ln_zps ) THEN   ;   CALL zps_hde_isf( kt, kjpt, zlap, zglu, zglv, zgui, zgvi )  ! both top & bottom
226      ELSEIF(             ln_zps ) THEN   ;   CALL zps_hde    ( kt, kjpt, zlap, zglu, zglv )              ! only bottom
227      ENDIF
228      !
229      SELECT CASE ( kldf )       !==  2nd laplacian applied to zlap (output in pta)  ==!
230      !
231      CASE ( np_blp    )               ! iso-level bilaplacian
232         CALL tra_ldf_lap  ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pta,      kjpt, 2 )
233      CASE ( np_blp_i  )               ! rotated   bilaplacian : standard operator (Madec)
234         CALL tra_ldf_iso  ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2 )
235      CASE ( np_blp_it )               ! rotated  bilaplacian : triad operator (griffies)
236         CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2 )
237      END SELECT
238      !
239   END SUBROUTINE tra_ldf_blp
240
241   !!==============================================================================
242END MODULE traldf_lap_blp
Note: See TracBrowser for help on using the repository browser.