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 branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90 @ 2399

Last change on this file since 2399 was 2399, checked in by gm, 13 years ago

v3.3beta: diaptr (poleward heat & salt transports) #759 : rewriting including dynamical allocation + DOCTOR names
  • Property svn:keywords set to Id
File size: 16.8 KB
Line 
1MODULE traldf_bilapg
2   !!==============================================================================
3   !!                       ***  MODULE  traldf_bilapg  ***
4   !! Ocean  tracers:  horizontal component of the lateral tracer mixing trend
5   !!==============================================================================
6   !! History : 8.0   ! 1997-07  (G. Madec)  Original code
7   !!   NEMO    1.0   ! 2002-08  (G. Madec)  F90: Free form and module
8   !!           3.3   ! 2010-06  (C. Ethe, G. Madec) Merge TRA-TRC
9   !!==============================================================================
10#if defined key_ldfslp   ||   defined key_esopa
11   !!----------------------------------------------------------------------
12   !!   'key_ldfslp'                  rotation of the lateral mixing tensor
13   !!----------------------------------------------------------------------
14   !!   tra_ldf_bilapg : update the tracer trend with the horizontal diffusion
15   !!                    using an horizontal biharmonic operator in s-coordinate
16   !!   ldfght         :  ???
17   !!----------------------------------------------------------------------
18   USE oce             ! ocean dynamics and tracers variables
19   USE dom_oce         ! ocean space and time domain variables
20   USE ldftra_oce      ! ocean active tracers: lateral physics
21   USE in_out_manager  ! I/O manager
22   USE ldfslp          ! iso-neutral slopes available
23   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
24   USE diaptr          ! poleward transport diagnostics
25   USE trc_oce         ! share passive tracers/Ocean variables
26
27   IMPLICIT NONE
28   PRIVATE
29
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   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
38   !! $Id$
39   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
40   !!----------------------------------------------------------------------
41CONTAINS
42
43   SUBROUTINE tra_ldf_bilapg( kt, cdtype, ptb, pta, kjpt )
44      !!----------------------------------------------------------------------
45      !!                 ***  ROUTINE tra_ldf_bilapg  ***
46      !!                   
47      !! ** Purpose :   Compute the before horizontal tracer 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      !!        ptb and multiply it by the eddy diffusivity coefficient
56      !!       (done by a call to ldfght routine, result in wk1 arrays).
57      !!      Applied the domain lateral boundary conditions by call to lbc_lnk
58      !!         -2- compute the geopotential harmonic operator applied to
59      !!      wk1 by a second call to ldfght routine (result in wk2)
60      !!      arrays).
61      !!         -3- Add this trend to the general trend
62      !!            pta = pta + wk2
63      !!
64      !! ** Action : - Update pta arrays with the before geopotential
65      !!               biharmonic mixing trend.
66      !!----------------------------------------------------------------------
67      INTEGER         , INTENT(in   )                      ::   kt       ! ocean time-step index
68      CHARACTER(len=3), INTENT(in   )                      ::   cdtype   ! =TRA or TRC (tracer indicator)
69      INTEGER         , INTENT(in   )                      ::   kjpt     ! number of tracers
70      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb      ! before and now tracer fields
71      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend
72      !!
73      INTEGER ::   ji, jj, jk, jn                 ! dummy loop indices
74      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt) ::   wk1, wk2   ! 4D workspace
75      !!----------------------------------------------------------------------
76
77      IF( kt == nit000 )  THEN
78         IF(lwp) WRITE(numout,*)
79         IF(lwp) WRITE(numout,*) 'tra_ldf_bilapg : horizontal biharmonic operator in s-coordinate on ', cdtype
80         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
81      ENDIF
82
83      ! 1. Laplacian of ptb * aht
84      ! -----------------------------
85      CALL ldfght( kt, cdtype, ptb, wk1, kjpt, 1 )      ! rotated harmonic operator applied to ptb and multiply by aht
86      !                                                 ! output in wk1
87      !
88      DO jn = 1, kjpt
89         CALL lbc_lnk( wk1(:,:,:,jn) , 'T', 1. )        ! Lateral boundary conditions on wk1   (unchanged sign)
90      END DO
91
92      ! 2. Bilaplacian of ptb
93      ! -------------------------
94      CALL ldfght( kt, cdtype, wk1, wk2, kjpt, 2 )      ! rotated harmonic operator applied to wk1 ; output in wk2
95
96
97      ! 3. Update the tracer trends                    (j-slab :   2, jpj-1)
98      ! ---------------------------
99      DO jn = 1, kjpt
100         DO jj = 2, jpjm1
101            DO jk = 1, jpkm1
102               DO ji = 2, jpim1
103                  ! add it to the general tracer trends
104                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + wk2(ji,jj,jk,jn)
105               END DO
106            END DO
107         END DO
108      END DO
109      !
110   END SUBROUTINE tra_ldf_bilapg
111
112
113   SUBROUTINE ldfght ( kt, cdtype, pt, plt, kjpt, kaht )
114      !!----------------------------------------------------------------------
115      !!                  ***  ROUTINE ldfght  ***
116      !!         
117      !! ** Purpose :   Apply a geopotential harmonic operator to (pt) and
118      !!      multiply it by the eddy diffusivity coefficient (if kaht=1).
119      !!      Routine only used in s-coordinates (l_sco=T) with bilaplacian
120      !!      operator (ln_traldf_bilap=T) acting along geopotential surfaces
121      !!      (ln_traldf_hor).
122      !!
123      !! ** Method  :   The harmonic operator rotated along geopotential
124      !!      surfaces is applied to (pt) using the slopes of geopotential
125      !!      surfaces computed in inildf routine. The result is provided in
126      !!      (plt,pls) arrays. It is computed in 2 steps:
127      !!
128      !!      First step: horizontal part of the operator. It is computed on
129      !!      ==========  pt as follows (idem on ps)
130      !!      horizontal fluxes :
131      !!         zftu = e2u*e3u/e1u di[ pt ] - e2u*uslp dk[ mi(mk(pt)) ]
132      !!         zftv = e1v*e3v/e2v dj[ pt ] - e1v*vslp dk[ mj(mk(pt)) ]
133      !!      take the horizontal divergence of the fluxes (no divided by
134      !!      the volume element :
135      !!         plt  = di-1[ zftu ] +  dj-1[ zftv ]
136      !!
137      !!      Second step: vertical part of the operator. It is computed on
138      !!      ===========  pt as follows (idem on ps)
139      !!      vertical fluxes :
140      !!         zftw = e1t*e2t/e3w * (wslpi^2+wslpj^2)  dk-1[ pt ]
141      !!              -     e2t     *       wslpi        di[ mi(mk(pt)) ]
142      !!              -     e1t     *       wslpj        dj[ mj(mk(pt)) ]
143      !!      take the vertical divergence of the fluxes add it to the hori-
144      !!      zontal component, divide the result by the volume element and
145      !!      if kaht=1, multiply by the eddy diffusivity coefficient:
146      !!         plt  = aht / (e1t*e2t*e3t) { plt + dk[ zftw ] }
147      !!      else:
148      !!         plt  =  1  / (e1t*e2t*e3t) { plt + dk[ zftw ] }
149      !!
150      !!----------------------------------------------------------------------
151      USE oce         , zftv => ua     ! use ua as workspace
152      !!
153      INTEGER         , INTENT(in )                              ::  kt      ! ocean time-step index
154      CHARACTER(len=3), INTENT(in )                              ::  cdtype  ! =TRA or TRC (tracer indicator)
155      INTEGER         , INTENT(in )                              ::  kjpt    !: dimension of
156      REAL(wp)        , INTENT(in ), DIMENSION(jpi,jpj,jpk,kjpt) ::  pt      ! tracer fields ( before for 1st call
157      !                                                         ! and laplacian of these fields for 2nd call.
158      REAL(wp)        , INTENT(out), DIMENSION(jpi,jpj,jpk,kjpt) ::  plt     !: partial harmonic operator applied to  pt  components except
159      !                                                             !: second order vertical derivative term 
160      INTEGER         , INTENT(in )                              ::  kaht    !: =1 multiply the laplacian by the eddy diffusivity coeff.
161      !                                                             !: =2 no multiplication
162      !!
163      INTEGER ::   ji, jj, jk,jn          ! dummy loop indices
164      !                                   ! temporary scalars
165      REAL(wp) ::  zabe1, zabe2, zmku, zmkv 
166      REAL(wp) ::  zbtr, ztah, ztav
167      REAL(wp) ::  zcof0, zcof1, zcof2, zcof3, zcof4
168      REAL(wp), DIMENSION(jpi,jpj) ::  zftu,  zdkt, zdk1t       ! workspace
169      REAL(wp), DIMENSION(jpi,jpk) ::  zftw, zdit, zdjt, zdj1t  !
170      !!----------------------------------------------------------------------
171
172      !
173      DO jn = 1, kjpt
174         !                               ! ********** !   ! ===============
175         DO jk = 1, jpkm1                ! First step !   ! Horizontal slab
176            !                            ! ********** !   ! ===============
177           
178            ! I.1 Vertical gradient of pt and ps at level jk and jk+1
179            ! -------------------------------------------------------
180            !     surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
181           
182            zdk1t(:,:) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
183            IF( jk == 1 ) THEN
184               zdkt(:,:) = zdk1t(:,:)
185            ELSE
186               zdkt(:,:) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk)
187            ENDIF
188
189
190            ! I.2 Horizontal fluxes
191            ! ---------------------
192           
193            DO jj = 1, jpjm1
194               DO ji = 1, jpim1
195                  zabe1 = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj)
196                  zabe2 = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj)
197                 
198                  zmku = 1./MAX( tmask(ji+1,jj,jk  )+tmask(ji,jj,jk+1)   &
199                     &          +tmask(ji+1,jj,jk+1)+tmask(ji,jj,jk  ),1. )
200                  zmkv = 1./MAX( tmask(ji,jj+1,jk  )+tmask(ji,jj,jk+1)   &
201                     &          +tmask(ji,jj+1,jk+1)+tmask(ji,jj,jk  ),1. )
202                 
203                  zcof1 = -e2u(ji,jj) * uslp(ji,jj,jk) * zmku
204                  zcof2 = -e1v(ji,jj) * vslp(ji,jj,jk) * zmkv
205                 
206                  zftu(ji,jj)= umask(ji,jj,jk) *   &
207                     &         (  zabe1 *( pt   (ji+1,jj,jk,jn) - pt(ji,jj,jk,jn) )   &
208                     &          + zcof1 *( zdkt (ji+1,jj) + zdk1t(ji,jj)     &
209                     &                    +zdk1t(ji+1,jj) + zdkt (ji,jj) )  )
210                 
211                  zftv(ji,jj,jk)= vmask(ji,jj,jk) *   &
212                     &         (  zabe2 *( pt(ji,jj+1,jk,jn) - pt(ji,jj,jk,jn) )   &
213                     &          + zcof2 *( zdkt (ji,jj+1) + zdk1t(ji,jj)     &
214                     &                    +zdk1t(ji,jj+1) + zdkt (ji,jj) )  )                 
215               END DO
216            END DO
217
218
219            ! I.3 Second derivative (divergence) (not divided by the volume)
220            ! ---------------------
221           
222            DO jj = 2 , jpjm1
223               DO ji = 2 , jpim1
224                  ztah = zftu(ji,jj) - zftu(ji-1,jj) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)
225                  plt(ji,jj,jk,jn) = ztah
226               END DO
227            END DO
228            !                                             ! ===============
229         END DO                                           !   End of slab
230         !                                                ! ===============
231         ! "Poleward" diffusive heat or salt transport
232         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( kaht == 2 ) .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN
233            IF( jn == jp_tem)   htr_ldf(:) = ptr_vj( zftv(:,:,:) )
234            IF( jn == jp_sal)   str_ldf(:) = ptr_vj( zftv(:,:,:) )
235         ENDIF
236
237         !                             ! ************ !   ! ===============
238         DO jj = 2, jpjm1              !  Second step !   ! Horizontal slab
239            !                          ! ************ !   ! ===============
240           
241            ! II.1 horizontal tracer gradient
242            ! -------------------------------
243           
244            DO jk = 1, jpk
245               DO ji = 1, jpim1
246                  zdit (ji,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj  ,jk,jn) ) * umask(ji,jj  ,jk)
247                  zdjt (ji,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj  ,jk,jn) ) * vmask(ji,jj  ,jk)
248                  zdj1t(ji,jk) = ( pt(ji  ,jj  ,jk,jn) - pt(ji,jj-1,jk,jn) ) * vmask(ji,jj-1,jk)
249               END DO
250            END DO
251
252
253            ! II.2 Vertical fluxes
254            ! --------------------
255           
256            ! Surface and bottom vertical fluxes set to zero
257            zftw(:, 1 ) = 0.e0
258            zftw(:,jpk) = 0.e0
259           
260            ! interior (2=<jk=<jpk-1)
261            DO jk = 2, jpkm1
262               DO ji = 2, jpim1
263                  zcof0 = e1t(ji,jj) * e2t(ji,jj) / fse3w(ji,jj,jk)   &
264                     &     * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)        &
265                     &        + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
266                 
267                  zmku = 1./MAX(  umask(ji  ,jj,jk-1)+umask(ji-1,jj,jk)   &
268                     &           +umask(ji-1,jj,jk-1)+umask(ji  ,jj,jk), 1. )
269                 
270                  zmkv = 1./MAX(  vmask(ji,jj  ,jk-1)+vmask(ji,jj-1,jk)   &
271                     &           +vmask(ji,jj-1,jk-1)+vmask(ji,jj  ,jk), 1. )
272                 
273                  zcof3 = - e2t(ji,jj) * wslpi (ji,jj,jk) * zmku
274                  zcof4 = - e1t(ji,jj) * wslpj (ji,jj,jk) * zmkv
275                 
276                  zftw(ji,jk) = tmask(ji,jj,jk) *   &
277                     &    (  zcof0 * ( pt   (ji,jj,jk-1,jn) - pt   (ji  ,jj,jk,jn) )   &
278                     &     + zcof3 * ( zdit (ji   ,jk-1   ) + zdit (ji-1,jk      )     &
279                     &                +zdit (ji-1 ,jk-1   ) + zdit (ji  ,jk      ) )   &
280                     &     + zcof4 * ( zdjt (ji   ,jk-1   ) + zdj1t(ji  ,jk)     &
281                     &                +zdj1t(ji   ,jk-1   ) + zdjt (ji  ,jk      ) )  )                 
282               END DO
283            END DO
284
285
286            ! II.3 Divergence of vertical fluxes added to the horizontal divergence
287            ! ---------------------------------------------------------------------
288           
289            IF( kaht == 1 ) THEN
290               ! multiply the laplacian by the eddy diffusivity coefficient
291               DO jk = 1, jpkm1
292                  DO ji = 2, jpim1
293                     ! eddy coef. divided by the volume element
294                     zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
295                     ! vertical divergence
296                     ztav = fsahtt(ji,jj,jk) * ( zftw(ji,jk) - zftw(ji,jk+1) )
297                     ! harmonic operator applied to (pt,ps) and multiply by aht
298                     plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr
299                  END DO
300               END DO
301            ELSEIF( kaht == 2 ) THEN
302               ! second call, no multiplication
303               DO jk = 1, jpkm1
304                  DO ji = 2, jpim1
305                     ! inverse of the volume element
306                     zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
307                     ! vertical divergence
308                     ztav = zftw(ji,jk) - zftw(ji,jk+1)
309                     ! harmonic operator applied to (pt,ps)
310                     plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr
311                  END DO
312               END DO
313            ELSE
314               IF(lwp) WRITE(numout,*) ' ldfght: kaht= 1 or 2, here =', kaht
315               IF(lwp) WRITE(numout,*) '         We stop'
316               STOP 'ldfght'
317            ENDIF
318            !                                             ! ===============
319         END DO                                           !   End of slab
320         !                                                ! ===============
321      END DO
322      !
323   END SUBROUTINE ldfght
324
325#else 
326   !!----------------------------------------------------------------------
327   !!   Dummy module :             NO rotation of the lateral mixing tensor
328   !!----------------------------------------------------------------------
329CONTAINS
330   SUBROUTINE tra_ldf_bilapg( kt )               ! Dummy routine
331      WRITE(*,*) 'tra_ldf_bilapg: You should not have seen this print! error?', kt
332   END SUBROUTINE tra_ldf_bilapg
333#endif
334
335   !!==============================================================================
336END MODULE traldf_bilapg
Note: See TracBrowser for help on using the repository browser.