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_bilapg.F90 in trunk/NEMO/TOP_SRC/TRP – NEMO

source: trunk/NEMO/TOP_SRC/TRP/trcldf_bilapg.F90 @ 197

Last change on this file since 197 was 186, checked in by opalod, 20 years ago

CL + CE : NEMO TRC_SRC start
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 15.7 KB
Line 
1MODULE trcldf_bilapg
2   !!==============================================================================
3   !!                       ***  MODULE  trcldf_bilapg  ***
4   !! Ocean passive tracers:  horizontal component of the lateral tracer mixing trend
5   !!==============================================================================
6#if key_passivetrc && defined key_ldfslp
7   !!----------------------------------------------------------------------
8   !!   'key_ldfslp'                  rotation of the lateral mixing tensor
9   !!----------------------------------------------------------------------
10   !!   trc_ldf_bilapg : update the tracer trend with the horizontal diffusion
11   !!                    using an horizontal biharmonic operator in s-coordinate
12   !!   ldfght         :  ???
13   !!----------------------------------------------------------------------
14   !! * Modules used
15   USE oce_trc             ! ocean dynamics and tracers variables
16   USE trc            ! ocean space and time domain variables
17   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
18
19   IMPLICIT NONE
20   PRIVATE
21
22   !! * Routine accessibility
23   PUBLIC trc_ldf_bilapg    ! routine called by step.F90
24
25   !! * Substitutions
26#  include "passivetrc_substitute.h90"
27   !!----------------------------------------------------------------------
28   !!   OPA 9.0 , LODYC-IPSL (2003)
29   !!----------------------------------------------------------------------
30   
31CONTAINS
32
33   SUBROUTINE trc_ldf_bilapg( kt )
34      !!----------------------------------------------------------------------
35      !!                 ***  ROUTINE trc_ldf_bilapg  ***
36      !!                   
37      !! ** Purpose :   Compute the before horizontal passive tracer diffusive
38      !!      trend and add it to the general trend of tracer equation.
39      !!
40      !! ** Method  :   The lateral diffusive trends is provided by a 4th order
41      !!      operator rotated along geopotential surfaces. It is computed
42      !!      using before fields (forward in time) and geopotential slopes
43      !!      computed in routine inildf.
44      !!         -1- compute the geopotential harmonic operator applied to
45      !!      trb and multiply it by the eddy diffusivity coefficient
46      !!      (done by a call to ldfght routine, result in wk1 array).
47      !!      Applied the domain lateral boundary conditions by call to lbc_lnk
48      !!         -2- compute the geopotential harmonic operator applied to
49      !!      wk1 by a second call to ldfght routine (result in wk3 array).
50      !!         -3- Add this trend to the general trend (ta,sa):
51      !!            tra = tra + wk3
52      !!
53      !! ** Action : - Update tra arrays with the before geopotential
54      !!               biharmonic mixing trend.
55      !!             - Save the trends  in trtrd ('key_trc_diatrd')
56      !!
57      !! History :
58      !!   8.0  !  97-07  (G. Madec)  Original code
59      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
60      !!   9.0  !  04-03  (C. Ethe)  adapted for passive tracers
61      !!----------------------------------------------------------------------
62      !! * Arguments
63      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
64
65      !! * Local declarations
66      INTEGER ::   ji, jj, jk,jn                 ! dummy loop indices
67      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra) ::   &
68         wk1, wk2               ! work array used for rotated biharmonic
69                                ! operator on tracers and/or momentum
70      !!----------------------------------------------------------------------
71
72      IF( kt == nittrc000 ) THEN
73         IF(lwp) WRITE(numout,*)
74         IF(lwp) WRITE(numout,*) 'trc_ldf_bilapg : horizontal biharmonic operator in s-coordinate'
75         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
76      ENDIF
77
78      ! 1. Laplacian of passive tracers trb * aht
79      ! -----------------------------
80      ! rotated harmonic operator applied to trb
81      ! and multiply by aht (output in wk1 )
82
83      CALL ldfght ( trb, wk1, 1 )
84
85      DO jn = 1, jptra
86      ! Lateral boundary conditions on wk1   (unchanged sign)
87         CALL lbc_lnk( wk1(:,:,:,jn) , 'T', 1. )
88      END DO
89
90      ! 2. Bilaplacian of trb
91      ! -------------------------
92      ! rotated harmonic operator applied to wk1 (output in wk2 )
93
94      CALL ldfght ( wk1, wk2, 2 )
95
96
97      DO jn = 1, jptra
98         ! 3. Update the tracer trends                    (j-slab :   2, jpj-1)
99         ! ---------------------------
100         !                                                ! ===============
101         DO jj = 2, jpjm1                                 !  Vertical slab
102            !                                             ! ===============
103            DO jk = 1, jpkm1
104               DO ji = 2, jpim1
105                  ! add it to the general tracer trends
106                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + wk2(ji,jj,jk,jn)
107#if defined key_trc_diatrd
108                  ! save the horizontal diffusive trends
109                  trtrd(ji,jj,jk,jn,3) = wk2(ji,jj,jk,jn)
110#endif
111               END DO
112            END DO
113            !                                             ! ===============
114         END DO                                           !   End of slab
115         !                                                ! ===============
116      END DO
117
118   END SUBROUTINE trc_ldf_bilapg
119
120
121   SUBROUTINE ldfght ( pt, plt, kaht )
122      !!----------------------------------------------------------------------
123      !!                  ***  ROUTINE ldfght  ***
124      !!         
125      !! ** Purpose :   Apply a geopotential harmonic operator to pt and
126      !!      multiply it by the eddy diffusivity coefficient (if kaht=1).
127      !!      Routine only used in s-coordinates (l_sco=T) with bilaplacian
128      !!      operator (ln_traldf_bilap=T) acting along geopotential surfaces
129      !!      (ln_traldf_hor).
130      !!
131      !! ** Method  :   The harmonic operator rotated along geopotential
132      !!      surfaces is applied to pt using the slopes of geopotential
133      !!      surfaces computed in inildf routine. The result is provided in
134      !!      plt arrays. It is computed in 2 steps:
135      !!
136      !!      First step: horizontal part of the operator. It is computed on
137      !!      ==========  pt as follows (idem on ps)
138      !!      horizontal fluxes :
139      !!         zftu = e2u*e3u/e1u di[ pt ] - e2u*uslp dk[ mi(mk(pt)) ]
140      !!         zftv = e1v*e3v/e2v dj[ pt ] - e1v*vslp dk[ mj(mk(pt)) ]
141      !!      take the horizontal divergence of the fluxes (no divided by
142      !!      the volume element :
143      !!         plt  = di-1[ zftu ] +  dj-1[ zftv ]
144      !!
145      !!      Second step: vertical part of the operator. It is computed on
146      !!      ===========  pt as follows (idem on ps)
147      !!      vertical fluxes :
148      !!         zftw = e1t*e2t/e3w * (wslpi^2+wslpj^2)  dk-1[ pt ]
149      !!              -     e2t     *       wslpi        di[ mi(mk(pt)) ]
150      !!              -     e1t     *       wslpj        dj[ mj(mk(pt)) ]
151      !!      take the vertical divergence of the fluxes add it to the hori-
152      !!      zontal component, divide the result by the volume element and
153      !!      if kaht=1, multiply by the eddy diffusivity coefficient:
154      !!         plt  = aht / (e1t*e2t*e3t) { plt + dk[ zftw ] }
155      !!      else:
156      !!         plt  =  1  / (e1t*e2t*e3t) { plt + dk[ zftw ] }
157      !!
158      !! * Action :
159      !!      'key_trdtra' defined: the trend is saved for diagnostics.
160      !!
161      !! History :
162      !!   8.0  !  97-07  (G. Madec)  Original code
163      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
164      !!   9.0  !  04-03  (C. Ethe)  adapted for passive tracers
165      !!----------------------------------------------------------------------
166      !! * Arguments
167      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra), INTENT( in  ) ::   &
168         pt               ! tracer fields before for 1st call
169      !                   ! and laplacian of these fields for 2nd call.
170      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra), INTENT( out ) ::   &
171         plt              ! partial harmonic operator applied to
172      !                   ! pt components except
173      !                   ! second order vertical derivative term)
174      INTEGER, INTENT( in ) ::   &
175         kaht             ! =1 multiply the laplacian by the eddy diffusivity coeff.
176      !                   ! =2 no multiplication
177
178      !! * Local declarations
179      INTEGER ::   ji, jj, jk,jn             ! dummy loop indices
180      REAL(wp) ::   &
181         zabe1, zabe2, zmku, zmkv,     &  ! temporary scalars
182         zbtr, ztah, ztav, &
183         zcof0, zcof1, zcof2,          &
184         zcof3, zcof4
185
186      REAL(wp), DIMENSION(jpi,jpj) ::   &
187         zftu, zftv ,                  &  ! workspace
188         zdkt, zdk1t
189
190      REAL(wp), DIMENSION(jpi,jpk) ::   &
191         zftw,                   &  ! workspace
192         zdit, zdjt, zdj1t
193
194      !!----------------------------------------------------------------------
195
196
197      DO jn = 1, jptra
198         !                               ! ********** !   ! ===============
199         DO jk = 1, jpkm1                ! First step !   ! Horizontal slab
200            !                            ! ********** !   ! ===============
201
202            ! I.1 Vertical gradient of pt and ps at level jk and jk+1
203            ! -------------------------------------------------------
204            !     surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
205
206            zdk1t(:,:) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
207
208            IF( jk == 1 ) THEN
209               zdkt(:,:) = zdk1t(:,:)
210            ELSE
211               zdkt(:,:) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk)
212            ENDIF
213
214
215            ! I.2 Horizontal fluxes
216            ! ---------------------
217
218            DO jj = 1, jpjm1
219               DO ji = 1, jpim1
220                  zabe1 = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj)
221                  zabe2 = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj)
222
223                  zmku=1./MAX( tmask(ji+1,jj,jk  )+tmask(ji,jj,jk+1)   &
224                     +tmask(ji+1,jj,jk+1)+tmask(ji,jj,jk  ),1. )
225                  zmkv=1./MAX( tmask(ji,jj+1,jk  )+tmask(ji,jj,jk+1)   &
226                     +tmask(ji,jj+1,jk+1)+tmask(ji,jj,jk  ),1. )
227
228                  zcof1= -e2u(ji,jj) * uslp(ji,jj,jk) * zmku
229                  zcof2= -e1v(ji,jj) * vslp(ji,jj,jk) * zmkv
230
231                  zftu(ji,jj)= umask(ji,jj,jk) *   &
232                     (  zabe1 *( pt(ji+1,jj,jk,jn) - pt(ji,jj,jk,jn) )   &
233                     + zcof1 *( zdkt (ji+1,jj) + zdk1t(ji,jj)     &
234                     +zdk1t(ji+1,jj) + zdkt (ji,jj) )  )
235
236                  zftv(ji,jj)= vmask(ji,jj,jk) *   &
237                     (  zabe2 *( pt(ji,jj+1,jk,jn) - pt(ji,jj,jk,jn) )   &
238                     + zcof2 *( zdkt (ji,jj+1) + zdk1t(ji,jj)     &
239                     +zdk1t(ji,jj+1) + zdkt (ji,jj) )  )
240               END DO
241            END DO
242
243
244            ! I.3 Second derivative (divergence) (not divided by the volume)
245            ! ---------------------
246
247            DO jj = 2 , jpjm1
248               DO ji = 2 , jpim1
249                  ztah = zftu(ji,jj) - zftu(ji-1,jj) + zftv(ji,jj) - zftv(ji,jj-1)
250                  plt(ji,jj,jk,jn) = ztah
251               END DO
252            END DO
253            !                                             ! ===============
254         END DO                                           !   End of slab
255         !                                                ! ===============
256
257
258         !                             ! ************ !   ! ===============
259         DO jj = 2, jpjm1              !  Second step !   ! Horizontal slab
260            !                          ! ************ !   ! ===============
261
262            ! II.1 horizontal tracer gradient
263            ! -------------------------------
264
265            DO jk = 1, jpk
266               DO ji = 1, jpim1
267                  zdit (ji,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj  ,jk,jn) ) * umask(ji,jj  ,jk)
268                  zdjt (ji,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj  ,jk,jn) ) * vmask(ji,jj  ,jk)
269                  zdj1t(ji,jk) = ( pt(ji  ,jj  ,jk,jn) - pt(ji,jj-1,jk,jn) ) * vmask(ji,jj-1,jk)
270               END DO
271            END DO
272
273
274            ! II.2 Vertical fluxes
275            ! --------------------
276
277            ! Surface and bottom vertical fluxes set to zero
278            zftw(:, 1 ) = 0.e0
279            zftw(:,jpk) = 0.e0
280
281            ! interior (2=<jk=<jpk-1)
282            DO jk = 2, jpkm1
283               DO ji = 2, jpim1
284                  zcof0 = e1t(ji,jj) * e2t(ji,jj) / fse3w(ji,jj,jk)   &
285                     * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)        &
286                     + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
287
288                  zmku =1./MAX(  umask(ji  ,jj,jk-1)+umask(ji-1,jj,jk)   &
289                     +umask(ji-1,jj,jk-1)+umask(ji  ,jj,jk), 1. )
290
291                  zmkv =1./MAX(  vmask(ji,jj  ,jk-1)+vmask(ji,jj-1,jk)   &
292                     +vmask(ji,jj-1,jk-1)+vmask(ji,jj  ,jk), 1. )
293
294                  zcof3 = - e2t(ji,jj) * wslpi (ji,jj,jk) * zmku
295                  zcof4 = - e1t(ji,jj) * wslpj (ji,jj,jk) * zmkv
296
297                  zftw(ji,jk) = tmask(ji,jj,jk) *   &
298                     (  zcof0 * ( pt  (ji,jj,jk-1,jn) - pt  (ji,jj,jk,jn) )   &
299                     + zcof3 * ( zdit (ji  ,jk-1) + zdit (ji-1,jk)     &
300                     +zdit (ji-1,jk-1) + zdit (ji  ,jk) )   &
301                     + zcof4 * ( zdjt (ji  ,jk-1) + zdj1t(ji  ,jk)     &
302                     +zdj1t(ji  ,jk-1) + zdjt (ji  ,jk) )  )
303               END DO
304            END DO
305
306
307            ! II.3 Divergence of vertical fluxes added to the horizontal divergence
308            ! ---------------------------------------------------------------------
309
310            IF( kaht == 1 ) THEN
311               ! multiply the laplacian by the eddy diffusivity coefficient
312               DO jk = 1, jpkm1
313                  DO ji = 2, jpim1
314                     ! eddy coef. divided by the volume element
315                     zbtr = fsahtrt(ji,jj,jk) / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
316                     ! vertical divergence
317                     ztav = zftw(ji,jk) - zftw(ji,jk+1)
318                     ! harmonic operator applied to (pt,ps) and multiply by aht
319                     plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr
320                  END DO
321               END DO
322            ELSEIF( kaht == 2 ) THEN
323               ! second call, no multiplication
324               DO jk = 1, jpkm1
325                  DO ji = 2, jpim1
326                     ! inverse of the volume element
327                     zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
328                     ! vertical divergence
329                     ztav = zftw(ji,jk) - zftw(ji,jk+1)
330                     ! harmonic operator applied to (pt,ps)
331                     plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr
332                  END DO
333               END DO
334            ELSE
335               IF(lwp) WRITE(numout,*) ' ldfght: kaht= 1 or 2, here =', kaht
336               IF(lwp) WRITE(numout,*) '         We stop'
337               STOP 'ldfght'
338            ENDIF
339            !                                             ! ===============
340         END DO                                           !   End of slab
341         !                                                ! ===============
342
343      END DO
344
345   END SUBROUTINE ldfght
346
347#else 
348   !!----------------------------------------------------------------------
349   !!   Dummy module :             NO rotation of the lateral mixing tensor
350   !!----------------------------------------------------------------------
351CONTAINS
352   SUBROUTINE trc_ldf_bilapg( kt )               ! Dummy routine
353      INTEGER, INTENT(in) :: kt
354      WRITE(*,*) 'trc_ldf_bilapg: You should not have seen this print! error?', kt
355   END SUBROUTINE trc_ldf_bilapg
356#endif
357
358   !!==============================================================================
359END MODULE trcldf_bilapg
Note: See TracBrowser for help on using the repository browser.