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