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 branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90 @ 4401

Last change on this file since 4401 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • Property svn:keywords set to Id
File size: 9.8 KB
Line 
1MODULE traldf_bilap
2   !!==============================================================================
3   !!                   ***  MODULE  traldf_bilap  ***
4   !! Ocean  tracers:  horizontal component of the lateral tracer mixing trend
5   !!==============================================================================
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   !!==============================================================================
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
24   USE ldftra_oce      ! ocean tracer   lateral physics
25   USE in_out_manager  ! I/O manager
26   USE ldfslp          ! iso-neutral slopes
27   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
28   USE diaptr          ! poleward transport diagnostics
29   USE trc_oce         ! share passive tracers/Ocean variables
30   USE lib_mpp         ! MPP library
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   tra_ldf_bilap   ! routine called by step.F90
36
37   !! * Control permutation of array indices
38#  include "oce_ftrans.h90"
39#  include "dom_oce_ftrans.h90"
40#  include "ldftra_oce_ftrans.h90"
41#  include "ldfslp_ftrans.h90"
42#  include "trc_oce_ftrans.h90"
43
44   !! * Substitutions
45#  include "domzgr_substitute.h90"
46#  include "ldftra_substitute.h90"
47#  include "ldfeiv_substitute.h90"
48#  include "vectopt_loop_substitute.h90"
49   !!----------------------------------------------------------------------
50   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
51   !! $Id$
52   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
53   !!----------------------------------------------------------------------
54CONTAINS
55 
56   SUBROUTINE tra_ldf_bilap( kt, cdtype, pgu, pgv,      &
57      &                                  ptb, pta, kjpt ) 
58      !!----------------------------------------------------------------------
59      !!                  ***  ROUTINE tra_ldf_bilap  ***
60      !!
61      !! ** Purpose :   Compute the before horizontal tracer diffusive
62      !!      trend and add it to the general trend of tracer equation.
63      !!
64      !! ** Method  :   4th order diffusive operator along model level surfaces
65      !!      evaluated using before fields (forward time scheme). The hor.
66      !!      diffusive trends  is given by:
67      !!      Laplacian of tb:
68      !!         zlt   = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(tb) ]
69      !!                                  + dj-1[ e1v*e3v/e2v dj(tb) ]  }
70      !!      Multiply by the eddy diffusivity coef. and insure lateral bc:
71      !!        zlt   = ahtt * zlt
72      !!        call to lbc_lnk
73      !!      Bilaplacian (laplacian of zlt):
74      !!         difft = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(zlt) ]
75      !!                                  + dj-1[ e1v*e3v/e2v dj(zlt) ]  }
76      !!
77      !!      Add this trend to the general trend
78      !!         (pta) = (pta) + ( difft )
79      !!
80      !! ** Action : - Update pta arrays with the before iso-level
81      !!               biharmonic mixing trend.
82      !!----------------------------------------------------------------------
83      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
84      USE oce     , ONLY:   ztu  => ua       , ztv  => va                           ! (ua,va) used as workspace
85      USE wrk_nemo, ONLY:   zeeu => wrk_2d_1 , zeev => wrk_2d_2 , zlt => wrk_2d_3   ! 2D workspace
86
87      !! DCSE_NEMO: need additional directives for renamed module variables
88!FTRANS ztu ztv :I :I :z
89
90      !!
91      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
92      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
93      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
94      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
95
96      !! DCSE_NEMO: This style defeats ftrans
97!     REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields
98!     REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
99!FTRANS ptb pta :I :I :z :
100      REAL(wp), INTENT(in   ) ::   ptb(jpi,jpj,jpk,kjpt)        ! before tracer fields
101      REAL(wp), INTENT(inout) ::   pta(jpi,jpj,jpk,kjpt)        ! tracer trend
102
103      !!
104      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
105      REAL(wp) ::  zbtr, ztra       ! local scalars
106      !!----------------------------------------------------------------------
107
108      IF( wrk_in_use(2, 1,2,3) ) THEN
109         CALL ctl_stop('tra_ldf_bilap: requested workspace arrays unavailable')   ;   RETURN
110      ENDIF
111
112      IF( kt == nit000 )  THEN
113         IF(lwp) WRITE(numout,*)
114         IF(lwp) WRITE(numout,*) 'tra_ldf_bilap : iso-level biharmonic operator on ', cdtype
115         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
116      ENDIF
117      !                                                          ! ===========
118      DO jn = 1, kjpt                                            ! tracer loop
119         !                                                       ! ===========
120         !                                               
121         DO jk = 1, jpkm1                                        ! Horizontal slab
122            !                                             
123            !                          !==  Initialization of metric arrays (for z- or s-coordinates)  ==!
124            DO jj = 1, jpjm1
125               DO ji = 1, fs_jpim1   ! vector opt.
126                  zeeu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
127                  zeev(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
128               END DO
129            END DO
130
131            !                          !==  Laplacian  ==!
132            !
133            DO jj = 1, jpjm1                 ! First derivative (gradient)
134               DO ji = 1, fs_jpim1   ! vector opt.
135                  ztu(ji,jj,jk) = zeeu(ji,jj) * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) )
136                  ztv(ji,jj,jk) = zeev(ji,jj) * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
137               END DO
138            END DO
139            IF( ln_zps ) THEN                ! set gradient at partial step level (last ocean level)
140               DO jj = 1, jpjm1
141                  DO ji = 1, jpim1
142                     IF( mbku(ji,jj) == jk )  ztu(ji,jj,jk) = zeeu(ji,jj) * pgu(ji,jj,jn)
143                     IF( mbkv(ji,jj) == jk )  ztv(ji,jj,jk) = zeev(ji,jj) * pgv(ji,jj,jn)
144                  END DO
145               END DO
146            ENDIF
147            DO jj = 2, jpjm1                 ! Second derivative (divergence) time the eddy diffusivity coefficient
148               DO ji = fs_2, fs_jpim1   ! vector opt.
149                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
150                  zlt(ji,jj) = fsahtt(ji,jj,jk) * zbtr * (   ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   &
151                     &                                     + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)   )
152               END DO
153            END DO
154            CALL lbc_lnk( zlt, 'T', 1. )     ! Lateral boundary conditions (unchanged sgn)
155
156            !                          !==  Bilaplacian  ==!
157            !
158            DO jj = 1, jpjm1                 ! third derivative (gradient)
159               DO ji = 1, fs_jpim1   ! vector opt.
160                  ztu(ji,jj,jk) = zeeu(ji,jj) * ( zlt(ji+1,jj  ) - zlt(ji,jj) )
161                  ztv(ji,jj,jk) = zeev(ji,jj) * ( zlt(ji  ,jj+1) - zlt(ji,jj) )
162               END DO
163            END DO
164            DO jj = 2, jpjm1                 ! fourth derivative (divergence) and add to the general tracer trend
165               DO ji = fs_2, fs_jpim1   ! vector opt.
166                  ! horizontal diffusive trends
167                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
168                  ztra = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
169                  ! add it to the general tracer trends
170                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
171               END DO
172            END DO
173            !                                             
174         END DO                                           ! Horizontal slab
175         !                                               
176         ! "zonal" mean lateral diffusive heat and salt transport
177         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
178           IF( jn == jp_tem )  htr_ldf(:) = ptr_vj( ztv(:,:,:) )
179           IF( jn == jp_sal )  str_ldf(:) = ptr_vj( ztv(:,:,:) )
180         ENDIF
181         !                                                ! ===========
182      END DO                                              ! tracer loop
183      !                                                   ! ===========
184      IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('tra_ldf_bilap: failed to release workspace arrays')
185      !
186   END SUBROUTINE tra_ldf_bilap
187
188   !!==============================================================================
189END MODULE traldf_bilap
Note: See TracBrowser for help on using the repository browser.