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.
trcldf_bilap.F90 in tags/nemo_v3_2/nemo_v3_2/NEMO/TOP_SRC/TRP – NEMO

source: tags/nemo_v3_2/nemo_v3_2/NEMO/TOP_SRC/TRP/trcldf_bilap.F90 @ 1878

Last change on this file since 1878 was 1878, checked in by flavoni, 14 years ago

initial test for nemogcm
File size: 11.0 KB
Line 
1MODULE trcldf_bilap
2   !!======================================================================
3   !!                   ***  MODULE  trcldf_bilap  ***
4   !! Ocean passive tracers:  horiz. component of the lateral tracer mixing trend
5   !!======================================================================
6   !! History :       !  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   !!                 !  00-05  (MA Foujols) add lbc for tracer trends
13   !!                 !  00-10  (MA Foujols E. Kestenare) use passive tracer coefficient
14   !!            8.5  !  02-08  (G. Madec)  F90: Free form and module
15   !!            9.0  !  04-03  (C. Ethe )  F90: Free form and module
16   !!                 !  07-02  (C. Deltel)  Diagnose ML trends for passive tracers
17   !!----------------------------------------------------------------------
18#if defined key_top
19   !!----------------------------------------------------------------------
20   !!   trc_ldf_bilap : update the tracer trend with the horizontal diffusion
21   !!                   using a iso-level biharmonic operator
22   !!----------------------------------------------------------------------
23   USE oce_trc         ! ocean dynamics and active tracers variables
24   USE trp_trc         ! ocean passive tracers variables
25   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
26   USE prtctl_trc      ! Print control for debbuging
27   USE trdmld_trc
28   USE trdmld_trc_oce     
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC trc_ldf_bilap   ! routine called by step.F90
34
35   !! * Substitutions
36#  include "top_substitute.h90"
37   !!----------------------------------------------------------------------
38   !!   TOP 1.0 , LOCEAN-IPSL (2005)
39   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/TOP_SRC/TRP/trcldf_bilap.F90,v 1.10 2006/09/12 11:10:14 opalod Exp $
40   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
41   !!----------------------------------------------------------------------
42
43CONTAINS
44   
45   SUBROUTINE trc_ldf_bilap( kt )
46      !!----------------------------------------------------------------------
47      !!                  ***  ROUTINE trc_ldf_bilap  ***
48      !!
49      !! ** Purpose :   Compute the before horizontal tracer tra diffusive
50      !!      trend and add it to the general trend of tracer equation.
51      !!
52      !! ** Method  :   4th order diffusive operator along model level surfaces
53      !!      evaluated using before fields (forward time scheme). The hor.
54      !!      diffusive trends of passive tracer is given by:
55      !!       * s-coordinate, the vertical scale
56      !!      factors e3. are inside the derivatives:
57      !!      Laplacian of trb:
58      !!         zlt   = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(trb) ]
59      !!                                  + dj-1[ e1v*e3v/e2v dj(trb) ]  }
60      !!      Multiply by the eddy diffusivity coef. and insure lateral bc:
61      !!        zlt   = ahtt * zlt
62      !!        call to lbc_lnk
63      !!      Bilaplacian (laplacian of zlt):
64      !!         difft = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(zlt) ]
65      !!                                  + dj-1[ e1v*e3v/e2v dj(zlt) ]  }
66      !!       * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:
67      !!      Laplacian of trb:
68      !!         zlt   = 1/(e1t*e2t) {  di-1[ e2u/e1u di(trb) ]
69      !!                              + dj-1[ e1v/e2v dj(trb) ] }
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) {  di-1[ e2u/e1u di(zlt) ]
75      !!                              + dj-1[ e1v/e2v dj(zlt) ]  }
76      !!
77      !!      Add this trend to the general trend tra :
78      !!         tra = tra + difft
79      !!
80      !! ** Action : - Update tra arrays with the before iso-level
81      !!               biharmonic mixing trend.
82      !!             - Save the trends ('key_trdmld_trc')
83      !!----------------------------------------------------------------------
84      USE oce_trc,   ztrtrd => ua        ! use ua as workspace
85      !!
86      INTEGER, INTENT( in ) ::   kt                             ! ocean time-step index
87      INTEGER ::   ji, jj, jk, jn                               ! dummy loop indices
88      INTEGER ::   iku, ikv                                     ! temporary integers
89      REAL(wp) ::   ztra                                        ! temporary scalars
90      REAL(wp), DIMENSION(jpi,jpj) ::   zeeu, zeev, zbtr, zlt   ! workspace
91      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztu, ztv            ! workspace
92      CHARACTER (len=22) :: charout
93      !!----------------------------------------------------------------------
94
95      IF( kt == nittrc000 ) THEN
96         IF(lwp) WRITE(numout,*)
97         IF(lwp) WRITE(numout,*) 'trc_ldf_bilap : iso-level biharmonic operator'
98         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
99      ENDIF
100      !                                                          ! ===========
101      DO jn = 1, jptra                                           ! tracer loop
102         !                                                       ! ===========
103         IF( l_trdtrc ) ztrtrd(:,:,:) = tra(:,:,:,jn)   ! save trends
104
105                                                          ! ===============
106         DO jk = 1, jpkm1                                 ! Horizontal slab
107            !                                             ! ===============
108
109            ! 0. Initialization of metric arrays (for z- or s-coordinates)
110            ! ----------------------------------
111
112            DO jj = 1, jpjm1
113               DO ji = 1, fs_jpim1   ! vector opt.
114#if ! defined key_zco
115                  ! s-coordinates, vertical scale factor are used
116                  zbtr(ji,jj) = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
117                  zeeu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
118                  zeev(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
119#else
120                  ! z-coordinates, no vertical scale factors
121                  zbtr(ji,jj) = 1. / ( e1t(ji,jj)*e2t(ji,jj) )
122                  zeeu(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk)
123                  zeev(ji,jj) = e1v(ji,jj) / e2v(ji,jj) * vmask(ji,jj,jk)
124#endif
125               END DO
126            END DO
127
128
129            ! 1. Laplacian
130            ! ------------
131
132            ! First derivative (gradient)
133            DO jj = 1, jpjm1
134               DO ji = 1, fs_jpim1   ! vector opt.
135                  ztu(ji,jj,jk) = zeeu(ji,jj) * ( trb(ji+1,jj  ,jk,jn) - trb(ji,jj,jk,jn) )
136                  ztv(ji,jj,jk) = zeev(ji,jj) * ( trb(ji  ,jj+1,jk,jn) - trb(ji,jj,jk,jn) )
137               END DO
138            END DO
139
140            IF( ln_zps ) THEN
141               DO jj = 1, jpj-1
142                  DO ji = 1, jpi-1
143                     ! last level
144                     iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1
145                     ikv = MIN ( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1
146                     IF( iku == jk ) THEN
147                        ztu(ji,jj,jk) = zeeu(ji,jj) * gtru(ji,jj,jn)
148                     ENDIF
149                     IF( ikv == jk ) THEN
150                        ztv(ji,jj,jk) = zeev(ji,jj) * gtrv(ji,jj,jn)
151                     ENDIF
152                  END DO
153               END DO
154            ENDIF
155
156            ! Second derivative (divergence)
157            DO jj = 2, jpjm1
158               DO ji = fs_2, fs_jpim1   ! vector opt.
159                  zlt(ji,jj) = zbtr(ji,jj) * (  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            ! Multiply by the eddy diffusivity coefficient
164            DO jj = 2, jpjm1
165               DO ji = fs_2, fs_jpim1   ! vector opt.
166                  zlt(ji,jj) = fsahtrt(ji,jj,jk) * zlt(ji,jj)
167               END DO
168            END DO
169
170            ! Lateral boundary conditions on the laplacian zlt   (unchanged sgn)
171            CALL lbc_lnk( zlt, 'T', 1. ) 
172
173            ! 2. Bilaplacian
174            ! --------------
175
176            ! third derivative (gradient)
177            DO jj = 1, jpjm1
178               DO ji = 1, fs_jpim1   ! vector opt.
179                  ztu(ji,jj,jk) = zeeu(ji,jj) * ( zlt(ji+1,jj  ) - zlt(ji,jj) )
180                  ztv(ji,jj,jk) = zeev(ji,jj) * ( zlt(ji  ,jj+1) - zlt(ji,jj) )
181               END DO
182            END DO
183
184            ! fourth derivative (divergence) and add to the general tracer trend
185            DO jj = 2, jpjm1
186               DO ji = fs_2, fs_jpim1   ! vector opt.
187                  ! horizontal diffusive trends
188                  ztra = zbtr(ji,jj) * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
189                  ! add it to the general tracer trends
190                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra
191#if defined key_trc_diatrd
192                  ! save the horizontal diffusive trends
193                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),4) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zbtr(ji,jj)
194                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),5) = (  ztv(ji,jj,jk) - ztv(ji-1,jj,jk) ) * zbtr(ji,jj)
195#endif
196
197               END DO
198            END DO
199            !                                             ! ===============
200         END DO                                           ! Horizontal slab
201         !                                                ! ===============
202#if defined key_trc_diatrd
203         ! Lateral boundary conditions on the laplacian zlt   (unchanged sgn)
204         IF (luttrd(jn)) CALL lbc_lnk( trtrd(:,:,:,ikeep(jn),5), 'T', 1. )
205#endif
206
207         IF( l_trdtrc ) THEN
208            ztrtrd(:,:,:) = tra(:,:,:,jn) - ztrtrd(:,:,:)
209            IF (luttrd(jn)) CALL trd_mod_trc( ztrtrd, jn, jptrc_trd_ldf, kt )    ! trends diags
210         END IF
211         !                                                       ! ===========
212      END DO                                                     ! tracer loop               
213      !                                                          ! ===========
214
215      IF( ln_ctl ) THEN    ! print mean trends (used for debugging)
216         WRITE(charout, FMT="('ldf - bilap')")
217         CALL prt_ctl_trc_info( charout )
218         CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd' )
219      ENDIF
220
221   END SUBROUTINE trc_ldf_bilap
222
223#else
224   !!----------------------------------------------------------------------
225   !!   Default option                                         Empty module
226   !!----------------------------------------------------------------------
227CONTAINS
228   SUBROUTINE trc_ldf_bilap( kt ) 
229      INTEGER, INTENT(in) :: kt
230      WRITE(*,*) 'trc_ldf_bilap: You should not have seen this print! error?', kt
231   END SUBROUTINE trc_ldf_bilap
232#endif
233   !!==============================================================================
234END MODULE trcldf_bilap
Note: See TracBrowser for help on using the repository browser.