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 tags/nemo_dev_x5/NEMO/OPA_SRC/TRA – NEMO

source: tags/nemo_dev_x5/NEMO/OPA_SRC/TRA/traldf_bilapg.F90 @ 1792

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

CT : UPDATE082 : Finalization of the poleward transport diagnostics
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 18.0 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 diaptr          ! 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      !! * Modules used
70      USE oce           , wk1 => ua,  &  ! use ua as workspace
71         &                wk2 => va      ! use va as workspace
72
73      !! * Arguments
74      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
75
76      !! * Local declarations
77      INTEGER ::   ji, jj, jk                 ! dummy loop indices
78      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
79         wk3, wk4                ! work array used for rotated biharmonic
80         !                       ! operator on tracers and/or momentum
81      !!----------------------------------------------------------------------
82
83      IF( kt == nit000 ) THEN
84         IF(lwp) WRITE(numout,*)
85         IF(lwp) WRITE(numout,*) 'tra_ldf_bilapg : horizontal biharmonic operator in s-coordinate'
86         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
87      ENDIF
88
89      ! 1. Laplacian of (tb,sb) * aht
90      ! -----------------------------
91      ! rotated harmonic operator applied to (tb,sb)
92      ! and multiply by aht (output in (wk1,wk2) )
93
94      CALL ldfght ( kt, tb, sb, wk1, wk2, 1 )
95
96
97      ! Lateral boundary conditions on (wk1,wk2)   (unchanged sign)
98      CALL lbc_lnk( wk1, 'T', 1. )   ;   CALL lbc_lnk( wk2, 'T', 1. )
99
100      ! 2. Bilaplacian of (tb,sb)
101      ! -------------------------
102      ! rotated harmonic operator applied to (wk1,wk2)
103      ! (output in (wk3,wk4) )
104
105      CALL ldfght ( kt, wk1, wk2, wk3, wk4, 2 )
106
107
108      ! 3. Update the tracer trends                    (j-slab :   2, jpj-1)
109      ! ---------------------------
110      !                                                ! ===============
111      DO jj = 2, jpjm1                                 !  Vertical slab
112         !                                             ! ===============
113         DO jk = 1, jpkm1
114            DO ji = 2, jpim1
115               ! add it to the general tracer trends
116               ta(ji,jj,jk) = ta(ji,jj,jk) + wk3(ji,jj,jk)
117               sa(ji,jj,jk) = sa(ji,jj,jk) + wk4(ji,jj,jk)
118#if defined key_trdtra || defined key_trdmld
119               ! save the horizontal diffusive trends
120               ttrd(ji,jj,jk,3) = wk3(ji,jj,jk)
121               strd(ji,jj,jk,3) = wk4(ji,jj,jk)
122#endif
123            END DO
124         END DO
125         !                                             ! ===============
126      END DO                                           !   End of slab
127      !                                                ! ===============
128   END SUBROUTINE tra_ldf_bilapg
129
130
131   SUBROUTINE ldfght ( kt, pt, ps, plt, pls, kaht )
132      !!----------------------------------------------------------------------
133      !!                  ***  ROUTINE ldfght  ***
134      !!         
135      !! ** Purpose :   Apply a geopotential harmonic operator to (pt,ps) 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,ps) using the slopes of geopotential
143      !!      surfaces computed in inildf routine. The result is provided in
144      !!      (plt,pls) 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      !!----------------------------------------------------------------------
175      !! * Arguments
176      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
177      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in  ) ::   &
178         pt, ps           ! tracer fields (before t and s for 1st call
179      !                   ! and laplacian of these fields for 2nd call.
180      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) ::   &
181         plt, pls         ! partial harmonic operator applied to
182      !                   ! pt & ps 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             ! dummy loop indices
190      REAL(wp) ::   &
191         zabe1, zabe2, zmku, zmkv,     &  ! temporary scalars
192         zbtr, ztah, zsah, ztav, zsav, &
193         zcof0, zcof1, zcof2,          &
194         zcof3, zcof4
195      REAL(wp), DIMENSION(jpi,jpj) ::  &
196         zftu, zfsu,                   &  ! workspace
197         zdkt, zdk1t,                  &
198         zdks, zdk1s
199      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
200         zftv, zfsv                       ! workspace (only v components for ptr)
201      REAL(wp), DIMENSION(jpi,jpk) ::  &
202         zftw, zfsw,                   &  ! workspace
203         zdit, zdjt, zdj1t,            &
204         zdis, zdjs, zdj1s
205      !!----------------------------------------------------------------------
206
207      !                               ! ********** !   ! ===============
208      DO jk = 1, jpkm1                ! First step !   ! Horizontal slab
209         !                            ! ********** !   ! ===============
210
211         ! I.1 Vertical gradient of pt and ps at level jk and jk+1
212         ! -------------------------------------------------------
213         !     surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
214
215         zdk1t(:,:) = ( pt(:,:,jk) - pt(:,:,jk+1) ) * tmask(:,:,jk+1)
216         zdk1s(:,:) = ( ps(:,:,jk) - ps(:,:,jk+1) ) * tmask(:,:,jk+1)
217
218         IF( jk == 1 ) THEN
219            zdkt(:,:) = zdk1t(:,:)
220            zdks(:,:) = zdk1s(:,:)
221         ELSE
222            zdkt(:,:) = ( pt(:,:,jk-1) - pt(:,:,jk) ) * tmask(:,:,jk)
223            zdks(:,:) = ( ps(:,:,jk-1) - ps(:,:,jk) ) * tmask(:,:,jk)
224         ENDIF
225
226
227         ! I.2 Horizontal fluxes
228         ! ---------------------
229
230         DO jj = 1, jpjm1
231            DO ji = 1, jpim1
232               zabe1 = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj)
233               zabe2 = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj)
234
235               zmku=1./MAX( tmask(ji+1,jj,jk  )+tmask(ji,jj,jk+1)   &
236                           +tmask(ji+1,jj,jk+1)+tmask(ji,jj,jk  ),1. )
237               zmkv=1./MAX( tmask(ji,jj+1,jk  )+tmask(ji,jj,jk+1)   &
238                           +tmask(ji,jj+1,jk+1)+tmask(ji,jj,jk  ),1. )
239
240               zcof1= -e2u(ji,jj) * uslp(ji,jj,jk) * zmku
241               zcof2= -e1v(ji,jj) * vslp(ji,jj,jk) * zmkv
242
243               zftu(ji,jj)= umask(ji,jj,jk) *   &
244                  (  zabe1 *( pt(ji+1,jj,jk) - pt(ji,jj,jk) )   &
245                   + zcof1 *( zdkt (ji+1,jj) + zdk1t(ji,jj)     &
246                             +zdk1t(ji+1,jj) + zdkt (ji,jj) )  )
247
248               zftv(ji,jj,jk)= vmask(ji,jj,jk) *   &
249                  (  zabe2 *( pt(ji,jj+1,jk) - pt(ji,jj,jk) )   &
250                   + zcof2 *( zdkt (ji,jj+1) + zdk1t(ji,jj)     &
251                             +zdk1t(ji,jj+1) + zdkt (ji,jj) )  )
252
253               zfsu(ji,jj)= umask(ji,jj,jk) *   &
254                  (  zabe1 *( ps(ji+1,jj,jk) - ps(ji,jj,jk) )   &
255                   + zcof1 *( zdks (ji+1,jj) + zdk1s(ji,jj)     &
256                             +zdk1s(ji+1,jj) + zdks (ji,jj) )  )
257
258               zfsv(ji,jj,jk)= vmask(ji,jj,jk) *   &
259                  (  zabe2 *( ps(ji,jj+1,jk) - ps(ji,jj,jk) )   &
260                   + zcof2 *( zdks (ji,jj+1) + zdk1s(ji,jj)     &
261                             +zdk1s(ji,jj+1) + zdks (ji,jj) )  )
262            END DO
263         END DO
264
265
266         ! I.3 Second derivative (divergence) (not divided by the volume)
267         ! ---------------------
268
269         DO jj = 2 , jpjm1
270            DO ji = 2 , jpim1
271               ztah = zftu(ji,jj) - zftu(ji-1,jj) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)
272               zsah = zfsu(ji,jj) - zfsu(ji-1,jj) + zfsv(ji,jj,jk) - zfsv(ji,jj-1,jk)
273               plt(ji,jj,jk) = ztah
274               pls(ji,jj,jk) = zsah
275            END DO
276         END DO
277         !                                             ! ===============
278      END DO                                           !   End of slab
279      !                                                ! ===============
280
281      ! "zonal" mean diffusive heat and salt transport
282      IF( ln_diaptr .AND. ( kaht == 2 ) .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
283         pht_ldf(:) = ptr_vj( zftv(:,:,:) )
284         pst_ldf(:) = ptr_vj( zfsv(:,:,:) )
285      ENDIF
286
287      !                             ! ************ !   ! ===============
288      DO jj = 2, jpjm1              !  Second step !   ! Horizontal slab
289         !                          ! ************ !   ! ===============
290
291         ! II.1 horizontal tracer gradient
292         ! -------------------------------
293
294         DO jk = 1, jpk
295            DO ji = 1, jpim1
296               zdit (ji,jk) = ( pt(ji+1,jj  ,jk) - pt(ji,jj  ,jk) ) * umask(ji,jj  ,jk)
297               zdis (ji,jk) = ( ps(ji+1,jj  ,jk) - ps(ji,jj  ,jk) ) * umask(ji,jj  ,jk)
298               zdjt (ji,jk) = ( pt(ji  ,jj+1,jk) - pt(ji,jj  ,jk) ) * vmask(ji,jj  ,jk)
299               zdjs (ji,jk) = ( ps(ji  ,jj+1,jk) - ps(ji,jj  ,jk) ) * vmask(ji,jj  ,jk)
300               zdj1t(ji,jk) = ( pt(ji  ,jj  ,jk) - pt(ji,jj-1,jk) ) * vmask(ji,jj-1,jk)
301               zdj1s(ji,jk) = ( ps(ji  ,jj  ,jk) - ps(ji,jj-1,jk) ) * vmask(ji,jj-1,jk)
302            END DO
303         END DO
304
305
306         ! II.2 Vertical fluxes
307         ! --------------------
308
309         ! Surface and bottom vertical fluxes set to zero
310         zftw(:, 1 ) = 0.e0
311         zfsw(:, 1 ) = 0.e0
312         zftw(:,jpk) = 0.e0
313         zfsw(:,jpk) = 0.e0
314
315         ! interior (2=<jk=<jpk-1)
316         DO jk = 2, jpkm1
317            DO ji = 2, jpim1
318               zcof0 = e1t(ji,jj) * e2t(ji,jj) / fse3w(ji,jj,jk)   &
319                     * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)        &
320                        + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
321
322               zmku =1./MAX(  umask(ji  ,jj,jk-1)+umask(ji-1,jj,jk)   &
323                             +umask(ji-1,jj,jk-1)+umask(ji  ,jj,jk), 1. )
324
325               zmkv =1./MAX(  vmask(ji,jj  ,jk-1)+vmask(ji,jj-1,jk)   &
326                             +vmask(ji,jj-1,jk-1)+vmask(ji,jj  ,jk), 1. )
327
328               zcof3 = - e2t(ji,jj) * wslpi (ji,jj,jk) * zmku
329               zcof4 = - e1t(ji,jj) * wslpj (ji,jj,jk) * zmkv
330
331               zftw(ji,jk) = tmask(ji,jj,jk) *   &
332                  (  zcof0 * ( pt  (ji,jj,jk-1) - pt  (ji,jj,jk) )   &
333                   + zcof3 * ( zdit (ji  ,jk-1) + zdit (ji-1,jk)     &
334                              +zdit (ji-1,jk-1) + zdit (ji  ,jk) )   &
335                   + zcof4 * ( zdjt (ji  ,jk-1) + zdj1t(ji  ,jk)     &
336                              +zdj1t(ji  ,jk-1) + zdjt (ji  ,jk) )  )
337
338               zfsw(ji,jk) = tmask(ji,jj,jk) *   &
339                  (  zcof0 * ( ps  (ji,jj,jk-1) - ps  (ji,jj,jk) )   &
340                   + zcof3 * ( zdis (ji  ,jk-1) + zdis (ji-1,jk)     &
341                              +zdis (ji-1,jk-1) + zdis (ji  ,jk) )   &
342                   + zcof4 * ( zdjs (ji  ,jk-1) + zdj1s(ji  ,jk)     &
343                              +zdj1s(ji  ,jk-1) + zdjs (ji  ,jk) )  )
344            END DO
345         END DO
346
347
348         ! II.3 Divergence of vertical fluxes added to the horizontal divergence
349         ! ---------------------------------------------------------------------
350
351         IF( kaht == 1 ) THEN
352            ! multiply the laplacian by the eddy diffusivity coefficient
353            DO jk = 1, jpkm1
354               DO ji = 2, jpim1
355                  ! eddy coef. divided by the volume element
356                  zbtr = fsahtt(ji,jj,jk) / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
357                  ! vertical divergence
358                  ztav = zftw(ji,jk) - zftw(ji,jk+1)
359                  zsav = zfsw(ji,jk) - zfsw(ji,jk+1)
360                  ! harmonic operator applied to (pt,ps) and multiply by aht
361                  plt(ji,jj,jk) = ( plt(ji,jj,jk) + ztav ) * zbtr
362                  pls(ji,jj,jk) = ( pls(ji,jj,jk) + zsav ) * zbtr
363               END DO
364            END DO
365         ELSEIF( kaht == 2 ) THEN
366            ! second call, no multiplication
367            DO jk = 1, jpkm1
368               DO ji = 2, jpim1
369                  ! inverse of the volume element
370                  zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
371                  ! vertical divergence
372                  ztav = zftw(ji,jk) - zftw(ji,jk+1)
373                  zsav = zfsw(ji,jk) - zfsw(ji,jk+1)
374                  ! harmonic operator applied to (pt,ps)
375                  plt(ji,jj,jk) = ( plt(ji,jj,jk) + ztav ) * zbtr
376                  pls(ji,jj,jk) = ( pls(ji,jj,jk) + zsav ) * zbtr
377               END DO
378            END DO
379         ELSE
380            IF(lwp) WRITE(numout,*) ' ldfght: kaht= 1 or 2, here =', kaht
381            IF(lwp) WRITE(numout,*) '         We stop'
382            STOP 'ldfght'
383         ENDIF
384         !                                             ! ===============
385      END DO                                           !   End of slab
386      !                                                ! ===============
387   END SUBROUTINE ldfght
388
389#else 
390   !!----------------------------------------------------------------------
391   !!   Dummy module :             NO rotation of the lateral mixing tensor
392   !!----------------------------------------------------------------------
393CONTAINS
394   SUBROUTINE tra_ldf_bilapg( kt )               ! Dummy routine
395      WRITE(*,*) 'tra_ldf_bilapg: You should not have seen this print! error?', kt
396   END SUBROUTINE tra_ldf_bilapg
397#endif
398
399   !!==============================================================================
400END MODULE traldf_bilapg
Note: See TracBrowser for help on using the repository browser.