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

Last change on this file since 2287 was 2287, checked in by smasson, 13 years ago

update licence of all NEMO files...

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