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

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

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