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.
traldf_bilapg.F90 in trunk/NEMO/OPA_SRC/TRA – NEMO

source: trunk/NEMO/OPA_SRC/TRA/traldf_bilapg.F90 @ 74

Last change on this file since 74 was 74, checked in by opalod, 20 years ago

CT : BUGFIX048 : Bug correction of Poleward Transport diagnostics
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 17.8 KB
Line 
1MODULE traldf_bilapg
2   !!==============================================================================
3   !!                       ***  MODULE  traldf_bilapg  ***
4   !! Ocean active tracers:  horizontal component of the lateral tracer mixing trend
5   !!==============================================================================
6#if defined key_ldfslp   ||   defined key_esopa
7   !!----------------------------------------------------------------------
8   !!   'key_ldfslp'                  rotation of the lateral mixing tensor
9   !!----------------------------------------------------------------------
10   !!   tra_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             ! ocean dynamics and tracers variables
16   USE dom_oce         ! ocean space and time domain variables
17   USE ldftra_oce      ! ocean active tracers: lateral physics
18   USE trdtra_oce      ! ocean active tracers: trend variables
19   USE in_out_manager  ! I/O manager
20   USE ldfslp          ! iso-neutral slopes available
21   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
22   USE ptr             ! poleward transport diagnostics
23
24   IMPLICIT NONE
25   PRIVATE
26
27   !! * Routine accessibility
28   PUBLIC tra_ldf_bilapg    ! routine called by step.F90
29
30   !! * Substitutions
31#  include "domzgr_substitute.h90"
32#  include "ldftra_substitute.h90"
33#  include "ldfeiv_substitute.h90"
34   !!----------------------------------------------------------------------
35   !!   OPA 9.0 , LODYC-IPSL (2003)
36   !!----------------------------------------------------------------------
37   
38CONTAINS
39
40   SUBROUTINE tra_ldf_bilapg( kt )
41      !!----------------------------------------------------------------------
42      !!                 ***  ROUTINE tra_ldf_bilapg  ***
43      !!                   
44      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
45      !!      trend and add it to the general trend of tracer equation.
46      !!
47      !! ** Method  :   The lateral diffusive trends is provided by a 4th order
48      !!      operator rotated along geopotential surfaces. It is computed
49      !!      using before fields (forward in time) and geopotential slopes
50      !!      computed in routine inildf.
51      !!         -1- compute the geopotential harmonic operator applied to
52      !!      (tb,sb) and multiply it by the eddy diffusivity coefficient
53      !!      (done by a call to ldfght routine, result in (wk1,wk2) arrays).
54      !!      Applied the domain lateral boundary conditions by call to lbc_lnk
55      !!         -2- compute the geopotential harmonic operator applied to
56      !!      (wk1,wk2) by a second call to ldfght routine (result in (wk3,wk4)
57      !!      arrays).
58      !!         -3- Add this trend to the general trend (ta,sa):
59      !!            (ta,sa) = (ta,sa) + (wk3,wk4)
60      !!
61      !! ** Action : - Update (ta,sa) arrays with the before geopotential
62      !!               biharmonic mixing trend.
63      !!             - Save the trends  in (ttrd,strd) ('key_diatrends')
64      !!
65      !! History :
66      !!   8.0  !  97-07  (G. Madec)  Original code
67      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
68      !!----------------------------------------------------------------------
69      !! * Arguments
70      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
71
72      !! * Local declarations
73      INTEGER ::   ji, jj, jk                 ! dummy loop indices
74      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
75         wk1, wk2,            &  ! work array used for rotated biharmonic
76         wk3, wk4                ! operator on tracers and/or momentum
77      !!----------------------------------------------------------------------
78
79      IF( kt == nit000 ) THEN
80         IF(lwp) WRITE(numout,*)
81         IF(lwp) WRITE(numout,*) 'tra_ldf_bilapg : horizontal biharmonic operator in s-coordinate'
82         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
83      ENDIF
84
85      ! 1. Laplacian of (tb,sb) * aht
86      ! -----------------------------
87      ! rotated harmonic operator applied to (tb,sb)
88      ! and multiply by aht (output in (wk1,wk2) )
89
90      CALL ldfght ( kt, tb, sb, wk1, wk2, 1 )
91
92
93      ! Lateral boundary conditions on (wk1,wk2)   (unchanged sign)
94      CALL lbc_lnk( wk1, 'T', 1. )   ;   CALL lbc_lnk( wk2, 'T', 1. )
95
96      ! 2. Bilaplacian of (tb,sb)
97      ! -------------------------
98      ! rotated harmonic operator applied to (wk1,wk2)
99      ! (output in (wk3,wk4) )
100
101      CALL ldfght ( kt, wk1, wk2, wk3, wk4, 2 )
102
103
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               ta(ji,jj,jk) = ta(ji,jj,jk) + wk3(ji,jj,jk)
113               sa(ji,jj,jk) = sa(ji,jj,jk) + wk4(ji,jj,jk)
114#if defined key_trdtra || defined key_trdmld
115               ! save the horizontal diffusive trends
116               ttrd(ji,jj,jk,3) = wk3(ji,jj,jk)
117               strd(ji,jj,jk,3) = wk4(ji,jj,jk)
118#endif
119            END DO
120         END DO
121         !                                             ! ===============
122      END DO                                           !   End of slab
123      !                                                ! ===============
124   END SUBROUTINE tra_ldf_bilapg
125
126
127   SUBROUTINE ldfght ( kt, pt, ps, plt, pls, kaht )
128      !!----------------------------------------------------------------------
129      !!                  ***  ROUTINE ldfght  ***
130      !!         
131      !! ** Purpose :   Apply a geopotential harmonic operator to (pt,ps) and
132      !!      multiply it by the eddy diffusivity coefficient (if kaht=1).
133      !!      Routine only used in s-coordinates (l_sco=T) with bilaplacian
134      !!      operator (ln_traldf_bilap=T) acting along geopotential surfaces
135      !!      (ln_traldf_hor).
136      !!
137      !! ** Method  :   The harmonic operator rotated along geopotential
138      !!      surfaces is applied to (pt,ps) using the slopes of geopotential
139      !!      surfaces computed in inildf routine. The result is provided in
140      !!      (plt,pls) arrays. It is computed in 2 steps:
141      !!
142      !!      First step: horizontal part of the operator. It is computed on
143      !!      ==========  pt as follows (idem on ps)
144      !!      horizontal fluxes :
145      !!         zftu = e2u*e3u/e1u di[ pt ] - e2u*uslp dk[ mi(mk(pt)) ]
146      !!         zftv = e1v*e3v/e2v dj[ pt ] - e1v*vslp dk[ mj(mk(pt)) ]
147      !!      take the horizontal divergence of the fluxes (no divided by
148      !!      the volume element :
149      !!         plt  = di-1[ zftu ] +  dj-1[ zftv ]
150      !!
151      !!      Second step: vertical part of the operator. It is computed on
152      !!      ===========  pt as follows (idem on ps)
153      !!      vertical fluxes :
154      !!         zftw = e1t*e2t/e3w * (wslpi^2+wslpj^2)  dk-1[ pt ]
155      !!              -     e2t     *       wslpi        di[ mi(mk(pt)) ]
156      !!              -     e1t     *       wslpj        dj[ mj(mk(pt)) ]
157      !!      take the vertical divergence of the fluxes add it to the hori-
158      !!      zontal component, divide the result by the volume element and
159      !!      if kaht=1, multiply by the eddy diffusivity coefficient:
160      !!         plt  = aht / (e1t*e2t*e3t) { plt + dk[ zftw ] }
161      !!      else:
162      !!         plt  =  1  / (e1t*e2t*e3t) { plt + dk[ zftw ] }
163      !!
164      !! * Action :
165      !!      'key_trdtra' defined: the trend is saved for diagnostics.
166      !!
167      !! History :
168      !!   8.0  !  97-07  (G. Madec)  Original code
169      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
170      !!----------------------------------------------------------------------
171      !! * Arguments
172      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
173      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in  ) ::   &
174         pt, ps           ! tracer fields (before t and s for 1st call
175      !                   ! and laplacian of these fields for 2nd call.
176      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) ::   &
177         plt, pls         ! partial harmonic operator applied to
178      !                   ! pt & ps 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      !! * Local declarations
185      INTEGER ::   ji, jj, jk             ! dummy loop indices
186      REAL(wp) ::   &
187         zabe1, zabe2, zmku, zmkv,     &  ! temporary scalars
188         zbtr, ztah, zsah, ztav, zsav, &
189         zcof0, zcof1, zcof2,          &
190         zcof3, zcof4
191      REAL(wp), DIMENSION(jpi,jpj) ::   &
192         zftu, zftv, zfsu, zfsv,       &  ! workspace
193         zdkt, zdk1t,                  &
194         zdks, zdk1s
195      REAL(wp), DIMENSION(jpi,jpk) ::   &
196         zftw, zfsw,                   &  ! workspace
197         zdit, zdjt, zdj1t,            &
198         zdis, zdjs, zdj1s
199      !!----------------------------------------------------------------------
200
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) - pt(:,:,jk+1) ) * tmask(:,:,jk+1)
210         zdk1s(:,:) = ( ps(:,:,jk) - ps(:,:,jk+1) ) * tmask(:,:,jk+1)
211
212         IF( jk == 1 ) THEN
213            zdkt(:,:) = zdk1t(:,:)
214            zdks(:,:) = zdk1s(:,:)
215         ELSE
216            zdkt(:,:) = ( pt(:,:,jk-1) - pt(:,:,jk) ) * tmask(:,:,jk)
217            zdks(:,:) = ( ps(:,:,jk-1) - ps(:,:,jk) ) * tmask(:,:,jk)
218         ENDIF
219
220
221         ! I.2 Horizontal fluxes
222         ! ---------------------
223
224         DO jj = 1, jpjm1
225            DO ji = 1, jpim1
226               zabe1 = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj)
227               zabe2 = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj)
228
229               zmku=1./MAX( tmask(ji+1,jj,jk  )+tmask(ji,jj,jk+1)   &
230                           +tmask(ji+1,jj,jk+1)+tmask(ji,jj,jk  ),1. )
231               zmkv=1./MAX( tmask(ji,jj+1,jk  )+tmask(ji,jj,jk+1)   &
232                           +tmask(ji,jj+1,jk+1)+tmask(ji,jj,jk  ),1. )
233
234               zcof1= -e2u(ji,jj) * uslp(ji,jj,jk) * zmku
235               zcof2= -e1v(ji,jj) * vslp(ji,jj,jk) * zmkv
236
237               zftu(ji,jj)= umask(ji,jj,jk) *   &
238                  (  zabe1 *( pt(ji+1,jj,jk) - pt(ji,jj,jk) )   &
239                   + zcof1 *( zdkt (ji+1,jj) + zdk1t(ji,jj)     &
240                             +zdk1t(ji+1,jj) + zdkt (ji,jj) )  )
241
242               zftv(ji,jj)= vmask(ji,jj,jk) *   &
243                  (  zabe2 *( pt(ji,jj+1,jk) - pt(ji,jj,jk) )   &
244                   + zcof2 *( zdkt (ji,jj+1) + zdk1t(ji,jj)     &
245                             +zdk1t(ji,jj+1) + zdkt (ji,jj) )  )
246
247               zfsu(ji,jj)= umask(ji,jj,jk) *   &
248                  (  zabe1 *( ps(ji+1,jj,jk) - ps(ji,jj,jk) )   &
249                   + zcof1 *( zdks (ji+1,jj) + zdk1s(ji,jj)     &
250                             +zdk1s(ji+1,jj) + zdks (ji,jj) )  )
251
252               zfsv(ji,jj)= vmask(ji,jj,jk) *   &
253                  (  zabe2 *( ps(ji,jj+1,jk) - ps(ji,jj,jk) )   &
254                   + zcof2 *( zdks (ji,jj+1) + zdk1s(ji,jj)     &
255                             +zdk1s(ji,jj+1) + zdks (ji,jj) )  )
256            END DO
257         END DO
258
259
260         ! I.3 Second derivative (divergence) (not divided by the volume)
261         ! ---------------------
262
263         DO jj = 2 , jpjm1
264            DO ji = 2 , jpim1
265               ztah = zftu(ji,jj) - zftu(ji-1,jj) + zftv(ji,jj) - zftv(ji,jj-1)
266               zsah = zfsu(ji,jj) - zfsu(ji-1,jj) + zfsv(ji,jj) - zfsv(ji,jj-1)
267               plt(ji,jj,jk) = ztah
268               pls(ji,jj,jk) = zsah
269            END DO
270         END DO
271         !                                             ! ===============
272      END DO                                           !   End of slab
273      !                                                ! ===============
274
275      IF( kaht == 2 .AND. MOD( kt, nf_ptr ) == 0 ) THEN
276#if defined key_diaptr
277         ! "zonal" mean diffusive heat and salt transport
278!         pht_ldf(:) = prt_vj( zftv(:,:) )
279!         pst_ldf(:) = prt_vj( zfsv(:,:) )
280        write(numout,cform_err)
281        nstop = nstop + 1
282#endif
283      ENDIF
284
285
286      !                             ! ************ !   ! ===============
287      DO jj = 2, jpjm1              !  Second step !   ! Horizontal slab
288         !                          ! ************ !   ! ===============
289
290         ! II.1 horizontal tracer gradient
291         ! -------------------------------
292
293         DO jk = 1, jpk
294            DO ji = 1, jpim1
295               zdit (ji,jk) = ( pt(ji+1,jj  ,jk) - pt(ji,jj  ,jk) ) * umask(ji,jj  ,jk)
296               zdis (ji,jk) = ( ps(ji+1,jj  ,jk) - ps(ji,jj  ,jk) ) * umask(ji,jj  ,jk)
297               zdjt (ji,jk) = ( pt(ji  ,jj+1,jk) - pt(ji,jj  ,jk) ) * vmask(ji,jj  ,jk)
298               zdjs (ji,jk) = ( ps(ji  ,jj+1,jk) - ps(ji,jj  ,jk) ) * vmask(ji,jj  ,jk)
299               zdj1t(ji,jk) = ( pt(ji  ,jj  ,jk) - pt(ji,jj-1,jk) ) * vmask(ji,jj-1,jk)
300               zdj1s(ji,jk) = ( ps(ji  ,jj  ,jk) - ps(ji,jj-1,jk) ) * vmask(ji,jj-1,jk)
301            END DO
302         END DO
303
304
305         ! II.2 Vertical fluxes
306         ! --------------------
307
308         ! Surface and bottom vertical fluxes set to zero
309         zftw(:, 1 ) = 0.e0
310         zfsw(:, 1 ) = 0.e0
311         zftw(:,jpk) = 0.e0
312         zfsw(:,jpk) = 0.e0
313
314         ! interior (2=<jk=<jpk-1)
315         DO jk = 2, jpkm1
316            DO ji = 2, jpim1
317               zcof0 = e1t(ji,jj) * e2t(ji,jj) / fse3w(ji,jj,jk)   &
318                     * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)        &
319                        + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
320
321               zmku =1./MAX(  umask(ji  ,jj,jk-1)+umask(ji-1,jj,jk)   &
322                             +umask(ji-1,jj,jk-1)+umask(ji  ,jj,jk), 1. )
323
324               zmkv =1./MAX(  vmask(ji,jj  ,jk-1)+vmask(ji,jj-1,jk)   &
325                             +vmask(ji,jj-1,jk-1)+vmask(ji,jj  ,jk), 1. )
326
327               zcof3 = - e2t(ji,jj) * wslpi (ji,jj,jk) * zmku
328               zcof4 = - e1t(ji,jj) * wslpj (ji,jj,jk) * zmkv
329
330               zftw(ji,jk) = tmask(ji,jj,jk) *   &
331                  (  zcof0 * ( pt  (ji,jj,jk-1) - pt  (ji,jj,jk) )   &
332                   + zcof3 * ( zdit (ji  ,jk-1) + zdit (ji-1,jk)     &
333                              +zdit (ji-1,jk-1) + zdit (ji  ,jk) )   &
334                   + zcof4 * ( zdjt (ji  ,jk-1) + zdj1t(ji  ,jk)     &
335                              +zdj1t(ji  ,jk-1) + zdjt (ji  ,jk) )  )
336
337               zfsw(ji,jk) = tmask(ji,jj,jk) *   &
338                  (  zcof0 * ( ps  (ji,jj,jk-1) - ps  (ji,jj,jk) )   &
339                   + zcof3 * ( zdis (ji  ,jk-1) + zdis (ji-1,jk)     &
340                              +zdis (ji-1,jk-1) + zdis (ji  ,jk) )   &
341                   + zcof4 * ( zdjs (ji  ,jk-1) + zdj1s(ji  ,jk)     &
342                              +zdj1s(ji  ,jk-1) + zdjs (ji  ,jk) )  )
343            END DO
344         END DO
345
346
347         ! II.3 Divergence of vertical fluxes added to the horizontal divergence
348         ! ---------------------------------------------------------------------
349
350         IF( kaht == 1 ) THEN
351            ! multiply the laplacian by the eddy diffusivity coefficient
352            DO jk = 1, jpkm1
353               DO ji = 2, jpim1
354                  ! eddy coef. divided by the volume element
355                  zbtr = fsahtt(ji,jj,jk) / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
356                  ! vertical divergence
357                  ztav = zftw(ji,jk) - zftw(ji,jk+1)
358                  zsav = zfsw(ji,jk) - zfsw(ji,jk+1)
359                  ! harmonic operator applied to (pt,ps) and multiply by aht
360                  plt(ji,jj,jk) = ( plt(ji,jj,jk) + ztav ) * zbtr
361                  pls(ji,jj,jk) = ( pls(ji,jj,jk) + zsav ) * zbtr
362               END DO
363            END DO
364         ELSEIF( kaht == 2 ) THEN
365            ! second call, no multiplication
366            DO jk = 1, jpkm1
367               DO ji = 2, jpim1
368                  ! inverse of the volume element
369                  zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
370                  ! vertical divergence
371                  ztav = zftw(ji,jk) - zftw(ji,jk+1)
372                  zsav = zfsw(ji,jk) - zfsw(ji,jk+1)
373                  ! harmonic operator applied to (pt,ps)
374                  plt(ji,jj,jk) = ( plt(ji,jj,jk) + ztav ) * zbtr
375                  pls(ji,jj,jk) = ( pls(ji,jj,jk) + zsav ) * zbtr
376               END DO
377            END DO
378         ELSE
379            IF(lwp) WRITE(numout,*) ' ldfght: kaht= 1 or 2, here =', kaht
380            IF(lwp) WRITE(numout,*) '         We stop'
381            STOP 'ldfght'
382         ENDIF
383         !                                             ! ===============
384      END DO                                           !   End of slab
385      !                                                ! ===============
386   END SUBROUTINE ldfght
387
388#else 
389   !!----------------------------------------------------------------------
390   !!   Dummy module :             NO rotation of the lateral mixing tensor
391   !!----------------------------------------------------------------------
392CONTAINS
393   SUBROUTINE tra_ldf_bilapg( kt )               ! Dummy routine
394      WRITE(*,*) 'tra_ldf_bilapg: You should not have seen this print! error?', kt
395   END SUBROUTINE tra_ldf_bilapg
396#endif
397
398   !!==============================================================================
399END MODULE traldf_bilapg
Note: See TracBrowser for help on using the repository browser.