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 @ 202

Last change on this file since 202 was 202, checked in by opalod, 19 years ago

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