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

Last change on this file since 941 was 941, checked in by cetlod, 13 years ago

phasing the passive tracer transport module to the new version of NEMO, see ticket 143

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