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_v3_2/nemo_v3_2/NEMO/OPA_SRC/TRA – NEMO

source: tags/nemo_v3_2/nemo_v3_2/NEMO/OPA_SRC/TRA/traldf_bilapg.F90 @ 1878

Last change on this file since 1878 was 1878, checked in by flavoni, 14 years ago

initial test for nemogcm
File size: 18.1 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   USE prtctl          ! Print control
25
26   IMPLICIT NONE
27   PRIVATE
28
29   !! * Routine accessibility
30   PUBLIC tra_ldf_bilapg    ! routine called by step.F90
31
32   !! * Substitutions
33#  include "domzgr_substitute.h90"
34#  include "ldftra_substitute.h90"
35#  include "ldfeiv_substitute.h90"
36   !!----------------------------------------------------------------------
37   !!   OPA 9.0 , LOCEAN-IPSL (2005)
38   !! $Id: traldf_bilapg.F90 1152 2008-06-26 14:11:13Z rblod $
39   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
40   !!----------------------------------------------------------------------
41   
42CONTAINS
43
44   SUBROUTINE tra_ldf_bilapg( kt )
45      !!----------------------------------------------------------------------
46      !!                 ***  ROUTINE tra_ldf_bilapg  ***
47      !!                   
48      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
49      !!      trend and add it to the general trend of tracer equation.
50      !!
51      !! ** Method  :   The lateral diffusive trends is provided by a 4th order
52      !!      operator rotated along geopotential surfaces. It is computed
53      !!      using before fields (forward in time) and geopotential slopes
54      !!      computed in routine inildf.
55      !!         -1- compute the geopotential harmonic operator applied to
56      !!      (tb,sb) and multiply it by the eddy diffusivity coefficient
57      !!      (done by a call to ldfght routine, result in (wk1,wk2) arrays).
58      !!      Applied the domain lateral boundary conditions by call to lbc_lnk
59      !!         -2- compute the geopotential harmonic operator applied to
60      !!      (wk1,wk2) by a second call to ldfght routine (result in (wk3,wk4)
61      !!      arrays).
62      !!         -3- Add this trend to the general trend (ta,sa):
63      !!            (ta,sa) = (ta,sa) + (wk3,wk4)
64      !!
65      !! ** Action : - Update (ta,sa) arrays with the before geopotential
66      !!               biharmonic mixing trend.
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), DIMENSION(jpi,jpj,jpk) ::   &
83         wk3, wk4                ! work array used for rotated biharmonic
84         !                       ! operator on tracers and/or momentum
85      !!----------------------------------------------------------------------
86
87      IF( kt == nit000 ) THEN
88         IF(lwp) WRITE(numout,*)
89         IF(lwp) WRITE(numout,*) 'tra_ldf_bilapg : horizontal biharmonic operator in s-coordinate'
90         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
91      ENDIF
92
93      ! 1. Laplacian of (tb,sb) * aht
94      ! -----------------------------
95      ! rotated harmonic operator applied to (tb,sb)
96      ! and multiply by aht (output in (wk1,wk2) )
97
98      CALL ldfght ( kt, tb, sb, wk1, wk2, 1 )
99
100
101      ! Lateral boundary conditions on (wk1,wk2)   (unchanged sign)
102      CALL lbc_lnk( wk1, 'T', 1. )   ;   CALL lbc_lnk( wk2, 'T', 1. )
103
104      ! 2. Bilaplacian of (tb,sb)
105      ! -------------------------
106      ! rotated harmonic operator applied to (wk1,wk2)
107      ! (output in (wk3,wk4) )
108
109      CALL ldfght ( kt, wk1, wk2, wk3, wk4, 2 )
110
111
112      ! 3. Update the tracer trends                    (j-slab :   2, jpj-1)
113      ! ---------------------------
114      !                                                ! ===============
115      DO jj = 2, jpjm1                                 !  Vertical slab
116         !                                             ! ===============
117         DO jk = 1, jpkm1
118            DO ji = 2, jpim1
119               ! add it to the general tracer trends
120               ta(ji,jj,jk) = ta(ji,jj,jk) + wk3(ji,jj,jk)
121               sa(ji,jj,jk) = sa(ji,jj,jk) + wk4(ji,jj,jk)
122            END DO
123         END DO
124         !                                             ! ===============
125      END DO                                           !   End of slab
126      !                                                ! ===============
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      !!but this should be done somewhere after
282      ! "zonal" mean diffusive heat and salt transport
283      IF( ln_diaptr .AND. ( kaht == 2 ) .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
284         pht_ldf(:) = ptr_vj( zftv(:,:,:) )
285         pst_ldf(:) = ptr_vj( zfsv(:,:,:) )
286      ENDIF
287
288      !                             ! ************ !   ! ===============
289      DO jj = 2, jpjm1              !  Second step !   ! Horizontal slab
290         !                          ! ************ !   ! ===============
291
292         ! II.1 horizontal tracer gradient
293         ! -------------------------------
294
295         DO jk = 1, jpk
296            DO ji = 1, jpim1
297               zdit (ji,jk) = ( pt(ji+1,jj  ,jk) - pt(ji,jj  ,jk) ) * umask(ji,jj  ,jk)
298               zdis (ji,jk) = ( ps(ji+1,jj  ,jk) - ps(ji,jj  ,jk) ) * umask(ji,jj  ,jk)
299               zdjt (ji,jk) = ( pt(ji  ,jj+1,jk) - pt(ji,jj  ,jk) ) * vmask(ji,jj  ,jk)
300               zdjs (ji,jk) = ( ps(ji  ,jj+1,jk) - ps(ji,jj  ,jk) ) * vmask(ji,jj  ,jk)
301               zdj1t(ji,jk) = ( pt(ji  ,jj  ,jk) - pt(ji,jj-1,jk) ) * vmask(ji,jj-1,jk)
302               zdj1s(ji,jk) = ( ps(ji  ,jj  ,jk) - ps(ji,jj-1,jk) ) * vmask(ji,jj-1,jk)
303            END DO
304         END DO
305
306
307         ! II.2 Vertical fluxes
308         ! --------------------
309
310         ! Surface and bottom vertical fluxes set to zero
311         zftw(:, 1 ) = 0.e0
312         zfsw(:, 1 ) = 0.e0
313         zftw(:,jpk) = 0.e0
314         zfsw(:,jpk) = 0.e0
315
316         ! interior (2=<jk=<jpk-1)
317         DO jk = 2, jpkm1
318            DO ji = 2, jpim1
319               zcof0 = e1t(ji,jj) * e2t(ji,jj) / fse3w(ji,jj,jk)   &
320                     * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)        &
321                        + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
322
323               zmku =1./MAX(  umask(ji  ,jj,jk-1)+umask(ji-1,jj,jk)   &
324                             +umask(ji-1,jj,jk-1)+umask(ji  ,jj,jk), 1. )
325
326               zmkv =1./MAX(  vmask(ji,jj  ,jk-1)+vmask(ji,jj-1,jk)   &
327                             +vmask(ji,jj-1,jk-1)+vmask(ji,jj  ,jk), 1. )
328
329               zcof3 = - e2t(ji,jj) * wslpi (ji,jj,jk) * zmku
330               zcof4 = - e1t(ji,jj) * wslpj (ji,jj,jk) * zmkv
331
332               zftw(ji,jk) = tmask(ji,jj,jk) *   &
333                  (  zcof0 * ( pt  (ji,jj,jk-1) - pt  (ji,jj,jk) )   &
334                   + zcof3 * ( zdit (ji  ,jk-1) + zdit (ji-1,jk)     &
335                              +zdit (ji-1,jk-1) + zdit (ji  ,jk) )   &
336                   + zcof4 * ( zdjt (ji  ,jk-1) + zdj1t(ji  ,jk)     &
337                              +zdj1t(ji  ,jk-1) + zdjt (ji  ,jk) )  )
338
339               zfsw(ji,jk) = tmask(ji,jj,jk) *   &
340                  (  zcof0 * ( ps  (ji,jj,jk-1) - ps  (ji,jj,jk) )   &
341                   + zcof3 * ( zdis (ji  ,jk-1) + zdis (ji-1,jk)     &
342                              +zdis (ji-1,jk-1) + zdis (ji  ,jk) )   &
343                   + zcof4 * ( zdjs (ji  ,jk-1) + zdj1s(ji  ,jk)     &
344                              +zdj1s(ji  ,jk-1) + zdjs (ji  ,jk) )  )
345            END DO
346         END DO
347
348
349         ! II.3 Divergence of vertical fluxes added to the horizontal divergence
350         ! ---------------------------------------------------------------------
351
352         IF( kaht == 1 ) THEN
353            ! multiply the laplacian by the eddy diffusivity coefficient
354            DO jk = 1, jpkm1
355               DO ji = 2, jpim1
356                  ! eddy coef. divided by the volume element
357                  zbtr = fsahtt(ji,jj,jk) / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
358                  ! vertical divergence
359                  ztav = zftw(ji,jk) - zftw(ji,jk+1)
360                  zsav = zfsw(ji,jk) - zfsw(ji,jk+1)
361                  ! harmonic operator applied to (pt,ps) and multiply by aht
362                  plt(ji,jj,jk) = ( plt(ji,jj,jk) + ztav ) * zbtr
363                  pls(ji,jj,jk) = ( pls(ji,jj,jk) + zsav ) * zbtr
364               END DO
365            END DO
366         ELSEIF( kaht == 2 ) THEN
367            ! second call, no multiplication
368            DO jk = 1, jpkm1
369               DO ji = 2, jpim1
370                  ! inverse of the volume element
371                  zbtr = 1. / ( 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)
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         ELSE
381            IF(lwp) WRITE(numout,*) ' ldfght: kaht= 1 or 2, here =', kaht
382            IF(lwp) WRITE(numout,*) '         We stop'
383            STOP 'ldfght'
384         ENDIF
385         !                                             ! ===============
386      END DO                                           !   End of slab
387      !                                                ! ===============
388   END SUBROUTINE ldfght
389
390#else 
391   !!----------------------------------------------------------------------
392   !!   Dummy module :             NO rotation of the lateral mixing tensor
393   !!----------------------------------------------------------------------
394CONTAINS
395   SUBROUTINE tra_ldf_bilapg( kt )               ! Dummy routine
396      WRITE(*,*) 'tra_ldf_bilapg: You should not have seen this print! error?', kt
397   END SUBROUTINE tra_ldf_bilapg
398#endif
399
400   !!==============================================================================
401END MODULE traldf_bilapg
Note: See TracBrowser for help on using the repository browser.