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

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

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