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/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90 @ 3152

Last change on this file since 3152 was 3116, checked in by cetlod, 13 years ago

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