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_bilap.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90 @ 2776

Last change on this file since 2776 was 2715, checked in by rblod, 13 years ago

First attempt to put dynamic allocation on the trunk
  • Property svn:keywords set to Id
File size: 9.3 KB
RevLine 
[3]1MODULE traldf_bilap
2   !!==============================================================================
3   !!                   ***  MODULE  traldf_bilap  ***
[2528]4   !! Ocean  tracers:  horizontal component of the lateral tracer mixing trend
[3]5   !!==============================================================================
[2528]6   !! History :  OPA  !  1991-11  (G. Madec)  Original code
7   !!                 !  1993-03  (M. Guyon)  symetrical conditions
8   !!                 !  1995-11  (G. Madec)  suppress volumetric scale factors
9   !!                 !  1996-01  (G. Madec)  statement function for e3
10   !!                 !  1996-01  (M. Imbard)  mpp exchange
11   !!                 !  1997-07  (G. Madec)  optimization, and ahtt
12   !!            8.5  !  2002-08  (G. Madec)  F90: Free form and module
13   !!   NEMO     1.0  !  2004-08  (C. Talandier) New trends organization
14   !!             -   !  2005-11  (G. Madec)  zps or sco as default option
15   !!            3.3  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
16   !!==============================================================================
[3]17
18   !!----------------------------------------------------------------------
19   !!   tra_ldf_bilap : update the tracer trend with the horizontal diffusion
20   !!                   using a iso-level biharmonic operator
21   !!----------------------------------------------------------------------
22   USE oce             ! ocean dynamics and active tracers
23   USE dom_oce         ! ocean space and time domain
[74]24   USE ldftra_oce      ! ocean tracer   lateral physics
[3]25   USE in_out_manager  ! I/O manager
26   USE ldfslp          ! iso-neutral slopes
27   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
[132]28   USE diaptr          ! poleward transport diagnostics
[2528]29   USE trc_oce         ! share passive tracers/Ocean variables
[2715]30   USE lib_mpp         ! MPP library
[3]31
32   IMPLICIT NONE
33   PRIVATE
34
[2528]35   PUBLIC   tra_ldf_bilap   ! routine called by step.F90
[3]36
37   !! * Substitutions
38#  include "domzgr_substitute.h90"
39#  include "ldftra_substitute.h90"
40#  include "ldfeiv_substitute.h90"
41#  include "vectopt_loop_substitute.h90"
42   !!----------------------------------------------------------------------
[2528]43   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
44   !! $Id$
45   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[3]46   !!----------------------------------------------------------------------
47CONTAINS
[2528]48 
49   SUBROUTINE tra_ldf_bilap( kt, cdtype, pgu, pgv,      &
50      &                                  ptb, pta, kjpt ) 
[3]51      !!----------------------------------------------------------------------
52      !!                  ***  ROUTINE tra_ldf_bilap  ***
53      !!
[2528]54      !! ** Purpose :   Compute the before horizontal tracer diffusive
[3]55      !!      trend and add it to the general trend of tracer equation.
56      !!
57      !! ** Method  :   4th order diffusive operator along model level surfaces
58      !!      evaluated using before fields (forward time scheme). The hor.
[2528]59      !!      diffusive trends  is given by:
[3]60      !!      Laplacian of tb:
61      !!         zlt   = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(tb) ]
62      !!                                  + dj-1[ e1v*e3v/e2v dj(tb) ]  }
63      !!      Multiply by the eddy diffusivity coef. and insure lateral bc:
64      !!        zlt   = ahtt * zlt
65      !!        call to lbc_lnk
66      !!      Bilaplacian (laplacian of zlt):
67      !!         difft = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(zlt) ]
68      !!                                  + dj-1[ e1v*e3v/e2v dj(zlt) ]  }
69      !!
[2528]70      !!      Add this trend to the general trend
71      !!         (pta) = (pta) + ( difft )
[3]72      !!
[2528]73      !! ** Action : - Update pta arrays with the before iso-level
[3]74      !!               biharmonic mixing trend.
[2528]75      !!----------------------------------------------------------------------
[2715]76      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
77      USE oce     , ONLY:   ztu  => ua       , ztv  => va                           ! (ua,va) used as workspace
78      USE wrk_nemo, ONLY:   zeeu => wrk_2d_1 , zeev => wrk_2d_2 , zlt => wrk_2d_3   ! 2D workspace
[3]79      !!
[2528]80      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
81      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
82      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
83      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
84      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields
85      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
86      !!
87      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
88      REAL(wp) ::  zbtr, ztra       ! local scalars
[3]89      !!----------------------------------------------------------------------
[132]90
[2715]91      IF( wrk_in_use(2, 1,2,3) ) THEN
92         CALL ctl_stop('tra_ldf_bilap: requested workspace arrays unavailable')   ;   RETURN
93      ENDIF
94
[2528]95      IF( kt == nit000 )  THEN
[3]96         IF(lwp) WRITE(numout,*)
[2528]97         IF(lwp) WRITE(numout,*) 'tra_ldf_bilap : iso-level biharmonic operator on ', cdtype
[3]98         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
99      ENDIF
[2528]100      !                                                          ! ===========
101      DO jn = 1, kjpt                                            ! tracer loop
102         !                                                       ! ===========
103         !                                               
104         DO jk = 1, jpkm1                                        ! Horizontal slab
105            !                                             
106            !                          !==  Initialization of metric arrays (for z- or s-coordinates)  ==!
[457]107            DO jj = 1, jpjm1
108               DO ji = 1, fs_jpim1   ! vector opt.
109                  zeeu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
110                  zeev(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
111               END DO
112            END DO
[3]113
[2528]114            !                          !==  Laplacian  ==!
115            !
116            DO jj = 1, jpjm1                 ! First derivative (gradient)
117               DO ji = 1, fs_jpim1   ! vector opt.
118                  ztu(ji,jj,jk) = zeeu(ji,jj) * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) )
119                  ztv(ji,jj,jk) = zeev(ji,jj) * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
120               END DO
[3]121            END DO
[2528]122            IF( ln_zps ) THEN                ! set gradient at partial step level (last ocean level)
123               DO jj = 1, jpjm1
124                  DO ji = 1, jpim1
125                     IF( mbku(ji,jj) == jk )  ztu(ji,jj,jk) = zeeu(ji,jj) * pgu(ji,jj,jn)
126                     IF( mbkv(ji,jj) == jk )  ztv(ji,jj,jk) = zeev(ji,jj) * pgv(ji,jj,jn)
127                  END DO
[457]128               END DO
[2528]129            ENDIF
130            DO jj = 2, jpjm1                 ! Second derivative (divergence) time the eddy diffusivity coefficient
131               DO ji = fs_2, fs_jpim1   ! vector opt.
132                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
133                  zlt(ji,jj) = fsahtt(ji,jj,jk) * zbtr * (   ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   &
134                     &                                     + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)   )
135               END DO
[3]136            END DO
[2528]137            CALL lbc_lnk( zlt, 'T', 1. )     ! Lateral boundary conditions (unchanged sgn)
[3]138
[2528]139            !                          !==  Bilaplacian  ==!
140            !
141            DO jj = 1, jpjm1                 ! third derivative (gradient)
142               DO ji = 1, fs_jpim1   ! vector opt.
143                  ztu(ji,jj,jk) = zeeu(ji,jj) * ( zlt(ji+1,jj  ) - zlt(ji,jj) )
144                  ztv(ji,jj,jk) = zeev(ji,jj) * ( zlt(ji  ,jj+1) - zlt(ji,jj) )
145               END DO
[3]146            END DO
[2528]147            DO jj = 2, jpjm1                 ! fourth derivative (divergence) and add to the general tracer trend
148               DO ji = fs_2, fs_jpim1   ! vector opt.
149                  ! horizontal diffusive trends
150                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
151                  ztra = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
152                  ! add it to the general tracer trends
153                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
[132]154               END DO
[3]155            END DO
[2528]156            !                                             
157         END DO                                           ! Horizontal slab
158         !                                               
159         ! "zonal" mean lateral diffusive heat and salt transport
160         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
161           IF( jn == jp_tem )  htr_ldf(:) = ptr_vj( ztv(:,:,:) )
162           IF( jn == jp_sal )  str_ldf(:) = ptr_vj( ztv(:,:,:) )
[457]163         ENDIF
[2528]164         !                                                ! ===========
165      END DO                                              ! tracer loop
166      !                                                   ! ===========
[2715]167      IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('tra_ldf_bilap: failed to release workspace arrays')
168      !
[3]169   END SUBROUTINE tra_ldf_bilap
170
171   !!==============================================================================
172END MODULE traldf_bilap
Note: See TracBrowser for help on using the repository browser.