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

Last change on this file since 1239 was 1175, checked in by cetlod, 16 years ago

update transport modules to take into account new trends organization, see ticket:248
  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 16.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   !! History :  8.0  !  97-07  (G. Madec)  Original code
7   !!            8.5  !  02-08  (G. Madec)  F90: Free form and module
8   !!            9.0  !  04-03  (C. Ethe)  adapted for passive tracers
9   !!                 !  07-02  (C. Deltel)  Diagnose ML trends for passive tracers
10   !!----------------------------------------------------------------------
11#if key_top && defined key_ldfslp
12   !!----------------------------------------------------------------------
13   !!   'key_ldfslp'                  rotation of the lateral mixing tensor
14   !!----------------------------------------------------------------------
15   !!   trc_ldf_bilapg : update the tracer trend with the horizontal diffusion
16   !!                    using an horizontal biharmonic operator in s-coordinate
17   !!   ldfght         :  ???
18   !!----------------------------------------------------------------------
19   USE oce_trc             ! ocean dynamics and tracers variables
20   USE trc                 ! ocean passive tracers variables
21   USE lbclnk              ! ocean lateral boundary condition (or mpp link)
22   USE prtctl_trc          ! Print control for debbuging
23   USE trp_trc
24   USE trdmld_trc
25   USE trdmld_trc_oce     
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC trc_ldf_bilapg    ! routine called by step.F90
31
32   !! * Substitutions
33#  include "top_substitute.h90"
34   !!----------------------------------------------------------------------
35   !!   TOP 1.0 , LOCEAN-IPSL (2005)
36   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/TOP_SRC/TRP/trcldf_bilapg.F90,v 1.9 2006/04/10 15:38:54 opalod Exp $
37   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39   
40CONTAINS
41
42   SUBROUTINE trc_ldf_bilapg( kt )
43      !!----------------------------------------------------------------------
44      !!                 ***  ROUTINE trc_ldf_bilapg  ***
45      !!                   
46      !! ** Purpose :   Compute the before horizontal passive tracer diffusive
47      !!      trend and add it to the general trend of tracer equation.
48      !!
49      !! ** Method  :   The lateral diffusive trends is provided by a 4th order
50      !!      operator rotated along geopotential surfaces. It is computed
51      !!      using before fields (forward in time) and geopotential slopes
52      !!      computed in routine inildf.
53      !!         -1- compute the geopotential harmonic operator applied to
54      !!      trb and multiply it by the eddy diffusivity coefficient
55      !!      (done by a call to ldfght routine, result in wk1 array).
56      !!      Applied the domain lateral boundary conditions by call to lbc_lnk
57      !!         -2- compute the geopotential harmonic operator applied to
58      !!      wk1 by a second call to ldfght routine (result in wk3 array).
59      !!         -3- Add this trend to the general trend (ta,sa):
60      !!            tra = tra + wk3
61      !!
62      !! ** Action : - Update tra arrays with the before geopotential
63      !!               biharmonic mixing trend.
64      !!----------------------------------------------------------------------
65      INTEGER, INTENT( in ) ::   kt                         ! ocean time-step index
66      INTEGER ::   ji, jj, jk, jn                           ! dummy loop indices
67      REAL(wp) ::   ztra                                    ! workspace   
68      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra) ::   wk1, wk2  ! work array used for rotated biharmonic
69      !                                                     ! operator on tracers and/or momentum
70      CHARACTER (len=22) ::   charout
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      DO jn = 1, jptra                                           ! tracer loop
98         !                                                       ! ===========
99
100         ! 3. Update the tracer trends                    (j-slab :   2, jpj-1)
101         ! ---------------------------
102         !                                                ! ===============
103         DO jj = 2, jpjm1                                 !  Vertical slab
104            !                                             ! ===============
105            DO jk = 1, jpkm1
106               DO ji = 2, jpim1
107                  ! add it to the general tracer trends
108                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + wk2(ji,jj,jk,jn)
109#if defined key_trc_diatrd
110                  ! save the horizontal diffusive trends
111                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),3) = wk2(ji,jj,jk,jn)
112#endif
113               END DO
114            END DO
115            !                                             ! ===============
116         END DO                                           !   End of slab
117         !                                                ! ===============
118         IF( l_trdtrc .AND. luttrd(jn) ) CALL trd_mod_trc( wk2(:,:,:,jn), jn, jptrc_trd_ldf, kt )
119
120         !                                                       ! ===========
121      END DO                                                     ! tracer loop
122      !                                                          ! ===========
123
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      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra), INTENT( in  ) ::   &
175         pt               ! tracer fields before for 1st call
176      !                   ! and laplacian of these fields for 2nd call.
177      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra), INTENT( out ) ::   &
178         plt              ! partial harmonic operator applied to
179      !                   ! pt components except
180      !                   ! second order vertical derivative term)
181      INTEGER, INTENT( in ) ::   &
182         kaht             ! =1 multiply the laplacian by the eddy diffusivity coeff.
183      !                   ! =2 no multiplication
184
185      INTEGER ::   ji, jj, jk,jn             ! dummy loop indices
186      REAL(wp) ::   &
187         zabe1, zabe2, zmku, zmkv,     &  ! temporary scalars
188         zbtr, ztah, ztav, &
189         zcof0, zcof1, zcof2,          &
190         zcof3, zcof4
191
192      REAL(wp), DIMENSION(jpi,jpj) ::   &
193         zftu, zftv ,                   &  ! workspace
194         zdkt, zdk1t
195
196      REAL(wp), DIMENSION(jpi,jpk) ::   &
197         zftw,                          &  ! workspace
198         zdit, zdjt, zdj1t
199      !!----------------------------------------------------------------------
200
201      DO jn = 1, jptra
202         !                               ! ********** !   ! ===============
203         DO jk = 1, jpkm1                ! First step !   ! Horizontal slab
204            !                            ! ********** !   ! ===============
205
206            ! I.1 Vertical gradient of pt and ps at level jk and jk+1
207            ! -------------------------------------------------------
208            !     surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
209
210            zdk1t(:,:) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
211
212            IF( jk == 1 ) THEN
213               zdkt(:,:) = zdk1t(:,:)
214            ELSE
215               zdkt(:,:) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk)
216            ENDIF
217
218
219            ! I.2 Horizontal fluxes
220            ! ---------------------
221
222            DO jj = 1, jpjm1
223               DO ji = 1, jpim1
224                  zabe1 = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj)
225                  zabe2 = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj)
226
227                  zmku=1./MAX( tmask(ji+1,jj,jk  )+tmask(ji,jj,jk+1)   &
228                              +tmask(ji+1,jj,jk+1)+tmask(ji,jj,jk  ),1. )
229                  zmkv=1./MAX( tmask(ji,jj+1,jk  )+tmask(ji,jj,jk+1)   &
230                              +tmask(ji,jj+1,jk+1)+tmask(ji,jj,jk  ),1. )
231
232                  zcof1= -e2u(ji,jj) * uslp(ji,jj,jk) * zmku
233                  zcof2= -e1v(ji,jj) * vslp(ji,jj,jk) * zmkv
234
235                  zftu(ji,jj)= umask(ji,jj,jk) *   &
236                     (  zabe1 *( pt(ji+1,jj,jk,jn) - pt(ji,jj,jk,jn) )   &
237                      + zcof1 *( zdkt (ji+1,jj) + zdk1t(ji,jj)           &
238                                +zdk1t(ji+1,jj) + zdkt (ji,jj) )  )
239
240                  zftv(ji,jj)= vmask(ji,jj,jk) *   &
241                     (  zabe2 *( pt(ji,jj+1,jk,jn) - pt(ji,jj,jk,jn) )   &
242                      + zcof2 *( zdkt (ji,jj+1) + zdk1t(ji,jj)           &
243                                +zdk1t(ji,jj+1) + zdkt (ji,jj) )  )
244               END DO
245            END DO
246
247
248            ! I.3 Second derivative (divergence) (not divided by the volume)
249            ! ---------------------
250
251            DO jj = 2 , jpjm1
252               DO ji = 2 , jpim1
253                  ztah = zftu(ji,jj) - zftu(ji-1,jj) + zftv(ji,jj) - zftv(ji,jj-1)
254                  plt(ji,jj,jk,jn) = ztah
255               END DO
256            END DO
257            !                                             ! ===============
258         END DO                                           !   End of slab
259         !                                                ! ===============
260
261
262         !                             ! ************ !   ! ===============
263         DO jj = 2, jpjm1              !  Second step !   ! Horizontal slab
264            !                          ! ************ !   ! ===============
265
266            ! II.1 horizontal tracer gradient
267            ! -------------------------------
268
269            DO jk = 1, jpk
270               DO ji = 1, jpim1
271                  zdit (ji,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj  ,jk,jn) ) * umask(ji,jj  ,jk)
272                  zdjt (ji,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj  ,jk,jn) ) * vmask(ji,jj  ,jk)
273                  zdj1t(ji,jk) = ( pt(ji  ,jj  ,jk,jn) - pt(ji,jj-1,jk,jn) ) * vmask(ji,jj-1,jk)
274               END DO
275            END DO
276
277
278            ! II.2 Vertical fluxes
279            ! --------------------
280
281            ! Surface and bottom vertical fluxes set to zero
282            zftw(:, 1 ) = 0.e0
283            zftw(:,jpk) = 0.e0
284
285            ! interior (2=<jk=<jpk-1)
286            DO jk = 2, jpkm1
287               DO ji = 2, jpim1
288                  zcof0 = e1t(ji,jj) * e2t(ji,jj) / fse3w(ji,jj,jk)   &
289                     * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)        &
290                        + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
291
292                  zmku =1./MAX(  umask(ji  ,jj,jk-1)+umask(ji-1,jj,jk)   &
293                                +umask(ji-1,jj,jk-1)+umask(ji  ,jj,jk), 1. )
294
295                  zmkv =1./MAX(  vmask(ji,jj  ,jk-1)+vmask(ji,jj-1,jk)   &
296                                +vmask(ji,jj-1,jk-1)+vmask(ji,jj  ,jk), 1. )
297
298                  zcof3 = - e2t(ji,jj) * wslpi (ji,jj,jk) * zmku
299                  zcof4 = - e1t(ji,jj) * wslpj (ji,jj,jk) * zmkv
300
301                  zftw(ji,jk) = tmask(ji,jj,jk) *   &
302                     (  zcof0 * ( pt  (ji,jj,jk-1,jn) - pt  (ji,jj,jk,jn) )   &
303                      + zcof3 * ( zdit (ji  ,jk-1) + zdit (ji-1,jk)           &
304                                 +zdit (ji-1,jk-1) + zdit (ji  ,jk) )         &
305                      + zcof4 * ( zdjt (ji  ,jk-1) + zdj1t(ji  ,jk)           &
306                                 +zdj1t(ji  ,jk-1) + zdjt (ji  ,jk) )  )
307               END DO
308            END DO
309
310
311            ! II.3 Divergence of vertical fluxes added to the horizontal divergence
312            ! ---------------------------------------------------------------------
313
314            IF( kaht == 1 ) THEN
315               ! multiply the laplacian by the eddy diffusivity coefficient
316               DO jk = 1, jpkm1
317                  DO ji = 2, jpim1
318                     ! eddy coef. divided by the volume element
319                     zbtr = fsahtrt(ji,jj,jk) / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
320                     ! vertical divergence
321                     ztav = zftw(ji,jk) - zftw(ji,jk+1)
322                     ! harmonic operator applied to (pt,ps) and multiply by aht
323                     plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr
324                  END DO
325               END DO
326            ELSEIF( kaht == 2 ) THEN
327               ! second call, no multiplication
328               DO jk = 1, jpkm1
329                  DO ji = 2, jpim1
330                     ! inverse of the volume element
331                     zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
332                     ! vertical divergence
333                     ztav = zftw(ji,jk) - zftw(ji,jk+1)
334                     ! harmonic operator applied to (pt,ps)
335                     plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr
336                  END DO
337               END DO
338            ELSE
339               IF(lwp) WRITE(numout,*) ' ldfght: kaht= 1 or 2, here =', kaht
340               IF(lwp) WRITE(numout,*) '         We stop'
341               STOP 'ldfght'
342            ENDIF
343            !                                             ! ===============
344         END DO                                           !   End of slab
345         !                                                ! ===============
346
347      END DO
348
349   END SUBROUTINE ldfght
350
351#else 
352   !!----------------------------------------------------------------------
353   !!   Dummy module :             NO rotation of the lateral mixing tensor
354   !!----------------------------------------------------------------------
355CONTAINS
356   SUBROUTINE trc_ldf_bilapg( kt )               ! Dummy routine
357      INTEGER, INTENT(in) :: kt
358      WRITE(*,*) 'trc_ldf_bilapg: You should not have seen this print! error?', kt
359   END SUBROUTINE trc_ldf_bilapg
360#endif
361
362   !!==============================================================================
363END MODULE trcldf_bilapg
Note: See TracBrowser for help on using the repository browser.