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/dev_001_GM/NEMO/OPA_SRC/TRA – NEMO

source: branches/dev_001_GM/NEMO/OPA_SRC/TRA/traldf_bilap.F90 @ 786

Last change on this file since 786 was 786, checked in by gm, 16 years ago

dev_001_GM - merge TRC-TRA on OPA only, trabbl & zpshde not done and trdmld not OK - compilation OK
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.6 KB
Line 
1MODULE traldf_bilap
2   !!==============================================================================
3   !!                   ***  MODULE  traldf_bilap  ***
4   !! Ocean active tracers:  horizontal component of the lateral tracer mixing trend
5   !!==============================================================================
6   !! History :  OPA  !  91-11  (G. Madec)  Original code
7   !!                 !  93-03  (M. Guyon)  symetrical conditions
8   !!                 !  95-11  (G. Madec)  suppress volumetric scale factors
9   !!                 !  96-01  (G. Madec)  statement function for e3
10   !!                 !  96-01  (M. Imbard)  mpp exchange
11   !!                 !  97-07  (G. Madec)  optimization, and ahtt
12   !!            NEMO !  02-08  (G. Madec)  F90: Free form and module
13   !!            1.0  !  04-08  (C. Talandier) New trends organization
14   !!                 !  05-11  (G. Madec)  zps or sco as default option
15   !!            2.4  !  08-01  (G. Madec) Merge TRA-TRC
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 dom_oce         ! ocean space and time domain
23   USE ldftra_oce      ! ocean tracer   lateral physics
24   USE in_out_manager  ! I/O manager
25   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
26   USE diaptr          ! poleward transport diagnostics
27   USE prtctl          ! Print control
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC tra_ldf_bilap   ! routine called by step.F90
33
34   !! * Substitutions
35#  include "domzgr_substitute.h90"
36#  include "ldftra_substitute.h90"
37#  include "ldfeiv_substitute.h90"
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40   !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008)
41   !! $Id:$
42   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44
45CONTAINS
46   
47   SUBROUTINE tra_ldf_bilap( kt, cdtype, ktra, pgtu, pgtv,   &
48      &                                        ptb , pta   )
49      !!----------------------------------------------------------------------
50      !!                  ***  ROUTINE tra_ldf_bilap  ***
51      !!
52      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
53      !!      trend and add it to the general trend of tracer equation.
54      !!
55      !! ** Method  :   4th order diffusive operator along model level surfaces
56      !!      evaluated using before fields (forward time scheme). The hor.
57      !!      diffusive trends of temperature (idem for salinity) is given by:
58      !!      Laplacian of tb:
59      !!         zlt   = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(tb) ]
60      !!                                  + dj-1[ e1v*e3v/e2v dj(tb) ]  }
61      !!      Multiply by the eddy diffusivity coef. and insure lateral bc:
62      !!        zlt   = ahtt * zlt
63      !!        call to lbc_lnk
64      !!      Bilaplacian (laplacian of zlt):
65      !!         difft = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(zlt) ]
66      !!                                  + dj-1[ e1v*e3v/e2v dj(zlt) ]  }
67      !!      Note: if key_zco defined, e3t=e3u=e3v, they are simplified.
68      !!
69      !!      Add this trend to the general trend (ta,sa):
70      !!         (ta,sa) = (ta,sa) + ( difft , diffs )
71      !!
72      !! ** Action : - Update (ta,sa) arrays with the before iso-level
73      !!               biharmonic mixing trend.
74      !!----------------------------------------------------------------------
75      INTEGER         , INTENT(in   )                         ::   kt              ! ocean time-step index
76      CHARACTER(len=3), INTENT(in   )                         ::   cdtype          ! =TRA or TRC (tracer indicator)
77      INTEGER         , INTENT(in   )                         ::   ktra            ! tracer index
78      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj)     ::   pgtu, pgtv      ! tracer gradient at pstep levels
79      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   ptb             ! before tracer field
80      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta             ! tracer trend
81      !!
82      INTEGER ::   ji, jj, jk             ! dummy loop indices
83      INTEGER ::   iku, ikv               ! temporary integers
84      REAL(wp), DIMENSION(jpi,jpj)     ::   zeeu, zeev, zbtr, zlt   ! 2D workspace
85      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztu, ztv          ! 3D workspace
86      !!----------------------------------------------------------------------
87
88      IF( kt == nit000 ) THEN
89         IF(lwp) WRITE(numout,*)
90         IF(lwp) WRITE(numout,*) 'tra_ldf_bilap : iso-level biharmonic operator'
91         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
92      ENDIF
93
94
95      !                                                ! ===============
96      DO jk = 1, jpkm1                                 ! Horizontal slab
97         !                                             ! ===============
98
99         ! 0. Initialization of metric arrays (for z- or s-coordinates)
100         ! ----------------------------------
101
102         DO jj = 1, jpjm1
103            DO ji = 1, fs_jpim1   ! vector opt.
104               zbtr(ji,jj) = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
105               zeeu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
106               zeev(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
107            END DO
108         END DO
109
110
111         ! 1. Laplacian
112         ! ------------
113
114         ! First derivative (gradient)
115         DO jj = 1, jpjm1
116            DO ji = 1, fs_jpim1   ! vector opt.
117               ztu(ji,jj,jk) = zeeu(ji,jj) * ( ptb(ji+1,jj  ,jk) - ptb(ji,jj,jk) )
118               ztv(ji,jj,jk) = zeev(ji,jj) * ( ptb(ji  ,jj+1,jk) - ptb(ji,jj,jk) )
119            END DO
120         END DO
121         IF( ln_zps ) THEN      ! set gradient at partial step level
122            DO jj = 1, jpjm1
123               DO ji = 1, jpim1
124                  ! last level
125                  iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1
126                  ikv = MIN ( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1
127                  IF( iku == jk )   ztu(ji,jj,jk) = zeeu(ji,jj) * pgtu(ji,jj)
128                  IF( ikv == jk )   ztv(ji,jj,jk) = zeev(ji,jj) * pgtv(ji,jj)
129               END DO
130            END DO
131         ENDIF
132
133         ! Second derivative (divergence) multiply by the eddy diffusivity coefficient
134         DO jj = 2, jpjm1
135            DO ji = fs_2, fs_jpim1   ! vector opt.
136               zlt(ji,jj) = fsahtt(ji,jj,jk) * zbtr(ji,jj)   &
137                  &                          * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
138            END DO
139         END DO
140
141!!gm   k-loop must be cut here and a 3D lbclnk used
142
143         ! Lateral boundary conditions on the laplacian (zlt)   (unchanged sgn)
144         CALL lbc_lnk( zlt, 'T', 1. ) 
145
146         ! 2. Bilaplacian
147         ! --------------
148
149         DO jj = 1, jpjm1                              ! third derivative (gradient)
150            DO ji = 1, fs_jpim1   ! vector opt.
151               ztu(ji,jj,jk) = zeeu(ji,jj) * ( zlt(ji+1,jj  ) - zlt(ji,jj) )
152               ztv(ji,jj,jk) = zeev(ji,jj) * ( zlt(ji  ,jj+1) - zlt(ji,jj) )
153            END DO
154         END DO
155
156         DO jj = 2, jpjm1                              ! 4th derivative (divergence) and add to the general tracer trend
157            DO ji = fs_2, fs_jpim1   ! vector opt.
158               pta(ji,jj,jk) = pta(ji,jj,jk) + zbtr(ji,jj)   &
159                  &                          * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
160            END DO
161         END DO
162         !                                             ! ===============
163      END DO                                           ! Horizontal slab
164      !                                                ! ===============
165
166
167      !                              ! "Poleward" lateral diffusive heat or salt transport
168      IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
169         IF( ktra == jp_tem)   pht_ldf(:) = ptr_vj( ztv(:,:,:) )
170         IF( ktra == jp_sal)   pst_ldf(:) = ptr_vj( ztv(:,:,:) )
171      ENDIF
172
173      !                              ! control print
174      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' ldf - bilap : ', mask1=tmask, clinfo3=cdtype )
175      !
176   END SUBROUTINE tra_ldf_bilap
177
178   !!==============================================================================
179END MODULE traldf_bilap
Note: See TracBrowser for help on using the repository browser.