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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90 @ 2696

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

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