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/2019/dev_r11943_MERGE_2019/src/OCE/TRA – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/traldf_lap_blp.F90 @ 13700

Last change on this file since 13700 was 12353, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. Additions to the do loop macro implementation: converted a few loops previously missed because they used jpi-1 instead of jpim1 etc.; changed internal macro names in do_loop_substitute.h90 to strings that are much more unlikely to appear in any future code elsewhere and removed the key_vectopt_loop option (and all related code) since the do loop macros have suppressed this option. These changes have been fully SETTE-tested and this branch should now be ready to go back to the trunk.
  • Property svn:keywords set to Id
File size: 13.1 KB
RevLine 
[5883]1MODULE traldf_lap_blp
[3]2   !!==============================================================================
[5883]3   !!                       ***  MODULE  traldf_lap_blp  ***
[5836]4   !! Ocean tracers:  lateral diffusivity trend  (laplacian and bilaplacian)
[3]5   !!==============================================================================
[5883]6   !! History :  3.7  ! 2014-01  (G. Madec, S. Masson)  Original code, re-entrant laplacian
[2528]7   !!----------------------------------------------------------------------
[3]8
9   !!----------------------------------------------------------------------
[5883]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
[3]12   !!----------------------------------------------------------------------
[5883]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
[7646]19   USE diaar5         ! AR5 diagnostics
[5883]20   USE trc_oce        ! share passive tracers/Ocean variables
21   USE zpshde         ! partial step: hor. derivative     (zps_hde routine)
[5836]22   !
[5883]23   USE in_out_manager ! I/O manager
[9019]24   USE iom            ! I/O library
[5883]25   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
26   USE lib_mpp        ! distribued memory computing library
27   USE timing         ! Timing
[3]28
29   IMPLICIT NONE
30   PRIVATE
31
[5883]32   PUBLIC   tra_ldf_lap   ! called by traldf.F90
33   PUBLIC   tra_ldf_blp   ! called by traldf.F90
[3]34
[7646]35   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
36   LOGICAL  ::   l_hst   ! flag to compute heat transport
37
[3]38   !! * Substitutions
[12340]39#  include "do_loop_substitute.h90"
[3]40   !!----------------------------------------------------------------------
[9598]41   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[2528]42   !! $Id$
[10068]43   !! Software governed by the CeCILL license (see ./LICENSE)
[3]44   !!----------------------------------------------------------------------
45CONTAINS
46
[11949]47   SUBROUTINE tra_ldf_lap( kt, Kmm, kit000, cdtype, pahu, pahv  ,               &
48      &                                             pgu , pgv   , pgui, pgvi,   &
49      &                                             pt  , pt_rhs, kjpt, kpass ) 
[3]50      !!----------------------------------------------------------------------
51      !!                  ***  ROUTINE tra_ldf_lap  ***
52      !!                   
53      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
54      !!      trend and add it to the general trend of tracer equation.
55      !!
56      !! ** Method  :   Second order diffusive operator evaluated using before
57      !!      fields (forward time scheme). The horizontal diffusive trends of
[2528]58      !!      the tracer is given by:
[5836]59      !!          difft = 1/(e1e2t*e3t) {  di-1[ pahu e2u*e3u/e1u di(tb) ]
60      !!                                 + dj-1[ pahv e1v*e3v/e2v dj(tb) ] }
[11949]61      !!      Add this trend to the general tracer trend pt_rhs :
62      !!          pt_rhs = pt_rhs + difft
[3]63      !!
[11949]64      !! ** Action  : - Update pt_rhs arrays with the before iso-level
[3]65      !!                harmonic mixing trend.
[2528]66      !!----------------------------------------------------------------------
67      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
[5836]68      INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index
[2528]69      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
70      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
[5836]71      INTEGER                              , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
[11949]72      INTEGER                              , INTENT(in   ) ::   Kmm        ! ocean time level index
[5836]73      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
[2528]74      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
[5836]75      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
[11949]76      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt         ! before tracer fields
77      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs     ! tracer trend
[2715]78      !
[5836]79      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
80      REAL(wp) ::   zsign            ! local scalars
[9019]81      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztu, ztv, zaheeu, zaheev
[3]82      !!----------------------------------------------------------------------
[3294]83      !
[5836]84      IF( kt == nit000 .AND. lwp )  THEN
85         WRITE(numout,*)
86         WRITE(numout,*) 'tra_ldf_lap : iso-level laplacian diffusion on ', cdtype, ', pass=', kpass
87         WRITE(numout,*) '~~~~~~~~~~~ '
[3]88      ENDIF
[5836]89      !
[7646]90      l_hst = .FALSE.
91      l_ptr = .FALSE.
[12193]92      IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) )     l_ptr = .TRUE.
[7646]93      IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
94         &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE.
95      !
[5836]96      !                                !==  Initialization of metric arrays used for all tracers  ==!
97      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
98      ELSE                    ;   zsign = -1._wp
99      ENDIF
[12340]100      DO_3D_10_10( 1, jpkm1 )
101         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!
102         zaheev(ji,jj,jk) = zsign * pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm)   !!gm   * vmask(ji,jj,jk)
103      END_3D
[5836]104      !
105      !                             ! =========== !
106      DO jn = 1, kjpt               ! tracer loop !
107         !                          ! =========== !   
108         !                               
[12340]109         DO_3D_10_10( 1, jpkm1 )
110            ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) )
111            ztv(ji,jj,jk) = zaheev(ji,jj,jk) * ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) )
112         END_3D
[5836]113         IF( ln_zps ) THEN                ! set gradient at bottom/top ocean level
[12340]114            DO_2D_10_10
115               ztu(ji,jj,mbku(ji,jj)) = zaheeu(ji,jj,mbku(ji,jj)) * pgu(ji,jj,jn)
116               ztv(ji,jj,mbkv(ji,jj)) = zaheev(ji,jj,mbkv(ji,jj)) * pgv(ji,jj,jn)
117            END_2D
[5836]118            IF( ln_isfcav ) THEN                ! top in ocean cavities only
[12340]119               DO_2D_10_10
120                  IF( miku(ji,jj) > 1 )   ztu(ji,jj,miku(ji,jj)) = zaheeu(ji,jj,miku(ji,jj)) * pgui(ji,jj,jn) 
121                  IF( mikv(ji,jj) > 1 )   ztv(ji,jj,mikv(ji,jj)) = zaheev(ji,jj,mikv(ji,jj)) * pgvi(ji,jj,jn) 
122               END_2D
[2528]123            ENDIF
[5836]124         ENDIF
125         !
[12340]126         DO_3D_00_00( 1, jpkm1 )
127            pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)     &
128               &                                      +    ztv(ji,jj,jk) - ztv(ji,jj-1,jk) )   &
129               &                                      / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) )
130         END_3D
[2528]131         !
[5836]132         !                             !== "Poleward" diffusive heat or salt transports  ==!
133         IF( ( kpass == 1 .AND. .NOT.ln_traldf_blp ) .OR.  &     !==  first pass only (  laplacian)  ==!
134             ( kpass == 2 .AND.      ln_traldf_blp ) ) THEN      !==  2nd   pass only (bilaplacian)  ==!
[7646]135
136            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', -ztv(:,:,:)  )
137            IF( l_hst )  CALL dia_ar5_hst( jn, 'ldf', -ztu(:,:,:), -ztv(:,:,:) )
[457]138         ENDIF
[5836]139         !                          ! ==================
140      END DO                        ! end of tracer loop
141      !                             ! ==================
[3294]142      !
[3]143   END SUBROUTINE tra_ldf_lap
[5836]144   
[5883]145
[11949]146   SUBROUTINE tra_ldf_blp( kt, Kmm, kit000, cdtype, pahu, pahv  ,             &
147      &                                             pgu , pgv   , pgui, pgvi, &
148      &                                             pt  , pt_rhs, kjpt, kldf )
[5883]149      !!----------------------------------------------------------------------
150      !!                 ***  ROUTINE tra_ldf_blp  ***
151      !!                   
152      !! ** Purpose :   Compute the before lateral tracer diffusive
153      !!      trend and add it to the general trend of tracer equation.
154      !!
155      !! ** Method  :   The lateral diffusive trends is provided by a bilaplacian
156      !!      operator applied to before field (forward in time).
157      !!      It is computed by two successive calls to laplacian routine
158      !!
159      !! ** Action :   pta   updated with the before rotated bilaplacian diffusion
160      !!----------------------------------------------------------------------
161      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
162      INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index
163      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
164      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
165      INTEGER                              , INTENT(in   ) ::   kldf       ! type of operator used
[11949]166      INTEGER                              , INTENT(in   ) ::   Kmm        ! ocean time level indices
[5883]167      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
168      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
169      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top levels
[11949]170      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt         ! before and now tracer fields
171      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs     ! tracer trend
[5883]172      !
173      INTEGER ::   ji, jj, jk, jn   ! dummy loop indices
[9019]174      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt) :: zlap         ! laplacian at t-point
175      REAL(wp), DIMENSION(jpi,jpj,    kjpt) :: zglu, zglv   ! bottom GRADh of the laplacian (u- and v-points)
176      REAL(wp), DIMENSION(jpi,jpj,    kjpt) :: zgui, zgvi   ! top    GRADh of the laplacian (u- and v-points)
[5883]177      !!---------------------------------------------------------------------
178      !
179      IF( kt == kit000 .AND. lwp )  THEN
180         WRITE(numout,*)
181         SELECT CASE ( kldf )
182         CASE ( np_blp    )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-level   bilaplacian operator on ', cdtype
183         CASE ( np_blp_i  )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (Standard)'
184         CASE ( np_blp_it )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (triad)'
185         END SELECT
186         WRITE(numout,*) '~~~~~~~~~~~'
187      ENDIF
188
189      zlap(:,:,:,:) = 0._wp
190      !
[11949]191      SELECT CASE ( kldf )       !==  1st laplacian applied to pt (output in zlap)  ==!
[5883]192      !
193      CASE ( np_blp    )               ! iso-level bilaplacian
[11949]194         CALL tra_ldf_lap  ( kt, Kmm, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, pt,     zlap, kjpt, 1 )
[5883]195      CASE ( np_blp_i  )               ! rotated   bilaplacian : standard operator (Madec)
[11949]196         CALL tra_ldf_iso  ( kt, Kmm, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, pt, pt, zlap, kjpt, 1 )
[5883]197      CASE ( np_blp_it )               ! rotated  bilaplacian : triad operator (griffies)
[11949]198         CALL tra_ldf_triad( kt, Kmm, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, pt, pt, zlap, kjpt, 1 )
[5883]199      END SELECT
200      !
[10425]201      CALL lbc_lnk( 'traldf_lap_blp', zlap(:,:,:,:) , 'T', 1. )     ! Lateral boundary conditions (unchanged sign)
[5883]202      !                                               ! Partial top/bottom cell: GRADh( zlap ) 
[11949]203      IF( ln_isfcav .AND. ln_zps ) THEN   ;   CALL zps_hde_isf( kt, Kmm, kjpt, zlap, zglu, zglv, zgui, zgvi )  ! both top & bottom
204      ELSEIF(             ln_zps ) THEN   ;   CALL zps_hde    ( kt, Kmm, kjpt, zlap, zglu, zglv )              ! only bottom
[5883]205      ENDIF
206      !
[11949]207      SELECT CASE ( kldf )       !==  2nd laplacian applied to zlap (output in pt_rhs)  ==!
[5883]208      !
209      CASE ( np_blp    )               ! iso-level bilaplacian
[11949]210         CALL tra_ldf_lap  ( kt, Kmm, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pt_rhs,         kjpt, 2 )
[5883]211      CASE ( np_blp_i  )               ! rotated   bilaplacian : standard operator (Madec)
[11949]212         CALL tra_ldf_iso  ( kt, Kmm, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pt    , pt_rhs, kjpt, 2 )
[5883]213      CASE ( np_blp_it )               ! rotated  bilaplacian : triad operator (griffies)
[11949]214         CALL tra_ldf_triad( kt, Kmm, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pt    , pt_rhs, kjpt, 2 )
[5883]215      END SELECT
216      !
217   END SUBROUTINE tra_ldf_blp
218
[3]219   !!==============================================================================
[5883]220END MODULE traldf_lap_blp
Note: See TracBrowser for help on using the repository browser.