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

Last change on this file since 1264 was 1264, checked in by cetlod, 15 years ago

clean TOP model routines to avoid warning when compiling, see ticket:303
  • 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), DIMENSION(jpi,jpj,jpk,jptra) ::   wk1, wk2  ! work array used for rotated biharmonic
68      !                                                     ! operator on tracers and/or momentum
69      CHARACTER (len=22) ::   charout
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      DO jn = 1, jptra                                           ! tracer loop
97         !                                                       ! ===========
98
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                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),3) = wk2(ji,jj,jk,jn)
111#endif
112               END DO
113            END DO
114            !                                             ! ===============
115         END DO                                           !   End of slab
116         !                                                ! ===============
117         IF( l_trdtrc .AND. luttrd(jn) ) CALL trd_mod_trc( wk2(:,:,:,jn), jn, jptrc_trd_ldf, kt )
118
119         !                                                       ! ===========
120      END DO                                                     ! tracer loop
121      !                                                          ! ===========
122
123
124      IF( ln_ctl ) THEN    ! print mean trends (used for debugging)
125         WRITE(charout, FMT="('ldf - bilapg')")
126         CALL prt_ctl_trc_info( charout )
127         CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd' )
128      ENDIF
129
130   END SUBROUTINE trc_ldf_bilapg
131
132
133   SUBROUTINE ldfght ( pt, plt, kaht )
134      !!----------------------------------------------------------------------
135      !!                  ***  ROUTINE ldfght  ***
136      !!         
137      !! ** Purpose :   Apply a geopotential harmonic operator to pt and
138      !!      multiply it by the eddy diffusivity coefficient (if kaht=1).
139      !!      Routine only used in s-coordinates (l_sco=T) with bilaplacian
140      !!      operator (ln_traldf_bilap=T) acting along geopotential surfaces
141      !!      (ln_traldf_hor).
142      !!
143      !! ** Method  :   The harmonic operator rotated along geopotential
144      !!      surfaces is applied to pt using the slopes of geopotential
145      !!      surfaces computed in inildf routine. The result is provided in
146      !!      plt arrays. It is computed in 2 steps:
147      !!
148      !!      First step: horizontal part of the operator. It is computed on
149      !!      ==========  pt as follows (idem on ps)
150      !!      horizontal fluxes :
151      !!         zftu = e2u*e3u/e1u di[ pt ] - e2u*uslp dk[ mi(mk(pt)) ]
152      !!         zftv = e1v*e3v/e2v dj[ pt ] - e1v*vslp dk[ mj(mk(pt)) ]
153      !!      take the horizontal divergence of the fluxes (no divided by
154      !!      the volume element :
155      !!         plt  = di-1[ zftu ] +  dj-1[ zftv ]
156      !!
157      !!      Second step: vertical part of the operator. It is computed on
158      !!      ===========  pt as follows (idem on ps)
159      !!      vertical fluxes :
160      !!         zftw = e1t*e2t/e3w * (wslpi^2+wslpj^2)  dk-1[ pt ]
161      !!              -     e2t     *       wslpi        di[ mi(mk(pt)) ]
162      !!              -     e1t     *       wslpj        dj[ mj(mk(pt)) ]
163      !!      take the vertical divergence of the fluxes add it to the hori-
164      !!      zontal component, divide the result by the volume element and
165      !!      if kaht=1, multiply by the eddy diffusivity coefficient:
166      !!         plt  = aht / (e1t*e2t*e3t) { plt + dk[ zftw ] }
167      !!      else:
168      !!         plt  =  1  / (e1t*e2t*e3t) { plt + dk[ zftw ] }
169      !!
170      !! * Action :
171      !!      'key_trdtra' defined: the trend is saved for diagnostics.
172      !!----------------------------------------------------------------------
173      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra), INTENT( in  ) ::   &
174         pt               ! tracer fields before for 1st call
175      !                   ! and laplacian of these fields for 2nd call.
176      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra), INTENT( out ) ::   &
177         plt              ! partial harmonic operator applied to
178      !                   ! pt components except
179      !                   ! second order vertical derivative term)
180      INTEGER, INTENT( in ) ::   &
181         kaht             ! =1 multiply the laplacian by the eddy diffusivity coeff.
182      !                   ! =2 no multiplication
183
184      INTEGER ::   ji, jj, jk,jn             ! dummy loop indices
185      REAL(wp) ::   &
186         zabe1, zabe2, zmku, zmkv,     &  ! temporary scalars
187         zbtr, ztah, ztav, &
188         zcof0, zcof1, zcof2,          &
189         zcof3, zcof4
190
191      REAL(wp), DIMENSION(jpi,jpj) ::   &
192         zftu, zftv ,                   &  ! workspace
193         zdkt, zdk1t
194
195      REAL(wp), DIMENSION(jpi,jpk) ::   &
196         zftw,                          &  ! workspace
197         zdit, zdjt, zdj1t
198      !!----------------------------------------------------------------------
199
200      DO jn = 1, jptra
201         !                               ! ********** !   ! ===============
202         DO jk = 1, jpkm1                ! First step !   ! Horizontal slab
203            !                            ! ********** !   ! ===============
204
205            ! I.1 Vertical gradient of pt and ps at level jk and jk+1
206            ! -------------------------------------------------------
207            !     surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
208
209            zdk1t(:,:) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
210
211            IF( jk == 1 ) THEN
212               zdkt(:,:) = zdk1t(:,:)
213            ELSE
214               zdkt(:,:) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk)
215            ENDIF
216
217
218            ! I.2 Horizontal fluxes
219            ! ---------------------
220
221            DO jj = 1, jpjm1
222               DO ji = 1, jpim1
223                  zabe1 = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj)
224                  zabe2 = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj)
225
226                  zmku=1./MAX( tmask(ji+1,jj,jk  )+tmask(ji,jj,jk+1)   &
227                              +tmask(ji+1,jj,jk+1)+tmask(ji,jj,jk  ),1. )
228                  zmkv=1./MAX( tmask(ji,jj+1,jk  )+tmask(ji,jj,jk+1)   &
229                              +tmask(ji,jj+1,jk+1)+tmask(ji,jj,jk  ),1. )
230
231                  zcof1= -e2u(ji,jj) * uslp(ji,jj,jk) * zmku
232                  zcof2= -e1v(ji,jj) * vslp(ji,jj,jk) * zmkv
233
234                  zftu(ji,jj)= umask(ji,jj,jk) *   &
235                     (  zabe1 *( pt(ji+1,jj,jk,jn) - pt(ji,jj,jk,jn) )   &
236                      + zcof1 *( zdkt (ji+1,jj) + zdk1t(ji,jj)           &
237                                +zdk1t(ji+1,jj) + zdkt (ji,jj) )  )
238
239                  zftv(ji,jj)= vmask(ji,jj,jk) *   &
240                     (  zabe2 *( pt(ji,jj+1,jk,jn) - pt(ji,jj,jk,jn) )   &
241                      + zcof2 *( zdkt (ji,jj+1) + zdk1t(ji,jj)           &
242                                +zdk1t(ji,jj+1) + zdkt (ji,jj) )  )
243               END DO
244            END DO
245
246
247            ! I.3 Second derivative (divergence) (not divided by the volume)
248            ! ---------------------
249
250            DO jj = 2 , jpjm1
251               DO ji = 2 , jpim1
252                  ztah = zftu(ji,jj) - zftu(ji-1,jj) + zftv(ji,jj) - zftv(ji,jj-1)
253                  plt(ji,jj,jk,jn) = ztah
254               END DO
255            END DO
256            !                                             ! ===============
257         END DO                                           !   End of slab
258         !                                                ! ===============
259
260
261         !                             ! ************ !   ! ===============
262         DO jj = 2, jpjm1              !  Second step !   ! Horizontal slab
263            !                          ! ************ !   ! ===============
264
265            ! II.1 horizontal tracer gradient
266            ! -------------------------------
267
268            DO jk = 1, jpk
269               DO ji = 1, jpim1
270                  zdit (ji,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj  ,jk,jn) ) * umask(ji,jj  ,jk)
271                  zdjt (ji,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj  ,jk,jn) ) * vmask(ji,jj  ,jk)
272                  zdj1t(ji,jk) = ( pt(ji  ,jj  ,jk,jn) - pt(ji,jj-1,jk,jn) ) * vmask(ji,jj-1,jk)
273               END DO
274            END DO
275
276
277            ! II.2 Vertical fluxes
278            ! --------------------
279
280            ! Surface and bottom vertical fluxes set to zero
281            zftw(:, 1 ) = 0.e0
282            zftw(:,jpk) = 0.e0
283
284            ! interior (2=<jk=<jpk-1)
285            DO jk = 2, jpkm1
286               DO ji = 2, jpim1
287                  zcof0 = e1t(ji,jj) * e2t(ji,jj) / fse3w(ji,jj,jk)   &
288                     * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)        &
289                        + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
290
291                  zmku =1./MAX(  umask(ji  ,jj,jk-1)+umask(ji-1,jj,jk)   &
292                                +umask(ji-1,jj,jk-1)+umask(ji  ,jj,jk), 1. )
293
294                  zmkv =1./MAX(  vmask(ji,jj  ,jk-1)+vmask(ji,jj-1,jk)   &
295                                +vmask(ji,jj-1,jk-1)+vmask(ji,jj  ,jk), 1. )
296
297                  zcof3 = - e2t(ji,jj) * wslpi (ji,jj,jk) * zmku
298                  zcof4 = - e1t(ji,jj) * wslpj (ji,jj,jk) * zmkv
299
300                  zftw(ji,jk) = tmask(ji,jj,jk) *   &
301                     (  zcof0 * ( pt  (ji,jj,jk-1,jn) - pt  (ji,jj,jk,jn) )   &
302                      + zcof3 * ( zdit (ji  ,jk-1) + zdit (ji-1,jk)           &
303                                 +zdit (ji-1,jk-1) + zdit (ji  ,jk) )         &
304                      + zcof4 * ( zdjt (ji  ,jk-1) + zdj1t(ji  ,jk)           &
305                                 +zdj1t(ji  ,jk-1) + zdjt (ji  ,jk) )  )
306               END DO
307            END DO
308
309
310            ! II.3 Divergence of vertical fluxes added to the horizontal divergence
311            ! ---------------------------------------------------------------------
312
313            IF( kaht == 1 ) THEN
314               ! multiply the laplacian by the eddy diffusivity coefficient
315               DO jk = 1, jpkm1
316                  DO ji = 2, jpim1
317                     ! eddy coef. divided by the volume element
318                     zbtr = fsahtrt(ji,jj,jk) / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
319                     ! vertical divergence
320                     ztav = zftw(ji,jk) - zftw(ji,jk+1)
321                     ! harmonic operator applied to (pt,ps) and multiply by aht
322                     plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr
323                  END DO
324               END DO
325            ELSEIF( kaht == 2 ) THEN
326               ! second call, no multiplication
327               DO jk = 1, jpkm1
328                  DO ji = 2, jpim1
329                     ! inverse of the volume element
330                     zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
331                     ! vertical divergence
332                     ztav = zftw(ji,jk) - zftw(ji,jk+1)
333                     ! harmonic operator applied to (pt,ps)
334                     plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr
335                  END DO
336               END DO
337            ELSE
338               IF(lwp) WRITE(numout,*) ' ldfght: kaht= 1 or 2, here =', kaht
339               IF(lwp) WRITE(numout,*) '         We stop'
340               STOP 'ldfght'
341            ENDIF
342            !                                             ! ===============
343         END DO                                           !   End of slab
344         !                                                ! ===============
345
346      END DO
347
348   END SUBROUTINE ldfght
349
350#else 
351   !!----------------------------------------------------------------------
352   !!   Dummy module :             NO rotation of the lateral mixing tensor
353   !!----------------------------------------------------------------------
354CONTAINS
355   SUBROUTINE trc_ldf_bilapg( kt )               ! Dummy routine
356      INTEGER, INTENT(in) :: kt
357      WRITE(*,*) 'trc_ldf_bilapg: You should not have seen this print! error?', kt
358   END SUBROUTINE trc_ldf_bilapg
359#endif
360
361   !!==============================================================================
362END MODULE trcldf_bilapg
Note: See TracBrowser for help on using the repository browser.