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

source: branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traldf_bilapg.F90 @ 2034

Last change on this file since 2034 was 2034, checked in by cetlod, 14 years ago

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