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

Last change on this file since 2633 was 2633, checked in by trackstand2, 13 years ago

Renamed wrk_use => wrk_in_use and wrk_release => wrk_not_released
  • 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
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      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
68      USE wrk_nemo, ONLY: wk1 => wrk_4d_1, wk2 => wrk_4d_2
69      INTEGER         , INTENT(in   )                      ::   kt       ! ocean time-step index
70      CHARACTER(len=3), INTENT(in   )                      ::   cdtype   ! =TRA or TRC (tracer indicator)
71      INTEGER         , INTENT(in   )                      ::   kjpt     ! number of tracers
72      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb      ! before and now tracer fields
73      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend
74      !!
75      INTEGER ::   ji, jj, jk, jn                 ! dummy loop indices
76      !!----------------------------------------------------------------------
77
78      IF(wrk_in_use(4, 1,2))THEN
79         CALL ctl_stop('tra_ldf_bilapg : requested workspace arrays unavailable.')
80         RETURN
81      END IF
82
83      IF( kt == nit000 )  THEN
84         IF(lwp) WRITE(numout,*)
85         IF(lwp) WRITE(numout,*) 'tra_ldf_bilapg : horizontal biharmonic operator in s-coordinate on ', cdtype
86         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
87      ENDIF
88
89      ! 1. Laplacian of ptb * aht
90      ! -----------------------------
91      CALL ldfght( kt, cdtype, ptb, wk1, kjpt, 1 )      ! rotated harmonic operator applied to ptb and multiply by aht
92      !                                                 ! output in wk1
93      !
94      DO jn = 1, kjpt
95         CALL lbc_lnk( wk1(:,:,:,jn) , 'T', 1. )        ! Lateral boundary conditions on wk1   (unchanged sign)
96      END DO
97
98      ! 2. Bilaplacian of ptb
99      ! -------------------------
100      CALL ldfght( kt, cdtype, wk1, wk2, kjpt, 2 )      ! rotated harmonic operator applied to wk1 ; output in wk2
101
102
103      ! 3. Update the tracer trends                    (j-slab :   2, jpj-1)
104      ! ---------------------------
105      DO jn = 1, kjpt
106         DO jj = 2, jpjm1
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         END DO
114      END DO
115      !
116      IF(wrk_not_released(4, 1,2))THEN
117         CALL ctl_stop('tra_ldf_bilapg : failed to release workspace arrays.')
118      END IF
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 oce         , zftv => ua     ! use ua as workspace
162      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released, wrk_in_use_xz, wrk_not_released_xz
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                          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.')
186         RETURN
187      END IF
188      !
189      DO jn = 1, kjpt
190         !                               ! ********** !   ! ===============
191         DO jk = 1, jpkm1                ! First step !   ! Horizontal slab
192            !                            ! ********** !   ! ===============
193           
194            ! I.1 Vertical gradient of pt and ps at level jk and jk+1
195            ! -------------------------------------------------------
196            !     surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
197           
198            zdk1t(:,:) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
199            IF( jk == 1 ) THEN
200               zdkt(:,:) = zdk1t(:,:)
201            ELSE
202               zdkt(:,:) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk)
203            ENDIF
204
205
206            ! I.2 Horizontal fluxes
207            ! ---------------------
208           
209            DO jj = 1, jpjm1
210               DO ji = 1, jpim1
211                  zabe1 = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj)
212                  zabe2 = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj)
213                 
214                  zmku = 1./MAX( tmask(ji+1,jj,jk  )+tmask(ji,jj,jk+1)   &
215                     &          +tmask(ji+1,jj,jk+1)+tmask(ji,jj,jk  ),1. )
216                  zmkv = 1./MAX( tmask(ji,jj+1,jk  )+tmask(ji,jj,jk+1)   &
217                     &          +tmask(ji,jj+1,jk+1)+tmask(ji,jj,jk  ),1. )
218                 
219                  zcof1 = -e2u(ji,jj) * uslp(ji,jj,jk) * zmku
220                  zcof2 = -e1v(ji,jj) * vslp(ji,jj,jk) * zmkv
221                 
222                  zftu(ji,jj)= umask(ji,jj,jk) *   &
223                     &         (  zabe1 *( pt   (ji+1,jj,jk,jn) - pt(ji,jj,jk,jn) )   &
224                     &          + zcof1 *( zdkt (ji+1,jj) + zdk1t(ji,jj)     &
225                     &                    +zdk1t(ji+1,jj) + zdkt (ji,jj) )  )
226                 
227                  zftv(ji,jj,jk)= vmask(ji,jj,jk) *   &
228                     &         (  zabe2 *( pt(ji,jj+1,jk,jn) - pt(ji,jj,jk,jn) )   &
229                     &          + zcof2 *( zdkt (ji,jj+1) + zdk1t(ji,jj)     &
230                     &                    +zdk1t(ji,jj+1) + zdkt (ji,jj) )  )                 
231               END DO
232            END DO
233
234
235            ! I.3 Second derivative (divergence) (not divided by the volume)
236            ! ---------------------
237           
238            DO jj = 2 , jpjm1
239               DO ji = 2 , jpim1
240                  ztah = zftu(ji,jj) - zftu(ji-1,jj) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)
241                  plt(ji,jj,jk,jn) = ztah
242               END DO
243            END DO
244            !                                             ! ===============
245         END DO                                           !   End of slab
246         !                                                ! ===============
247         ! "Poleward" diffusive heat or salt transport
248         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( kaht == 2 ) .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN
249            IF( jn == jp_tem)   htr_ldf(:) = ptr_vj( zftv(:,:,:) )
250            IF( jn == jp_sal)   str_ldf(:) = ptr_vj( zftv(:,:,:) )
251         ENDIF
252
253         !                             ! ************ !   ! ===============
254         DO jj = 2, jpjm1              !  Second step !   ! Horizontal slab
255            !                          ! ************ !   ! ===============
256           
257            ! II.1 horizontal tracer gradient
258            ! -------------------------------
259           
260            DO jk = 1, jpk
261               DO ji = 1, jpim1
262                  zdit (ji,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj  ,jk,jn) ) * umask(ji,jj  ,jk)
263                  zdjt (ji,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj  ,jk,jn) ) * vmask(ji,jj  ,jk)
264                  zdj1t(ji,jk) = ( pt(ji  ,jj  ,jk,jn) - pt(ji,jj-1,jk,jn) ) * vmask(ji,jj-1,jk)
265               END DO
266            END DO
267
268
269            ! II.2 Vertical fluxes
270            ! --------------------
271           
272            ! Surface and bottom vertical fluxes set to zero
273            zftw(:, 1 ) = 0.e0
274            zftw(:,jpk) = 0.e0
275           
276            ! interior (2=<jk=<jpk-1)
277            DO jk = 2, jpkm1
278               DO ji = 2, jpim1
279                  zcof0 = e1t(ji,jj) * e2t(ji,jj) / fse3w(ji,jj,jk)   &
280                     &     * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)        &
281                     &        + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
282                 
283                  zmku = 1./MAX(  umask(ji  ,jj,jk-1)+umask(ji-1,jj,jk)   &
284                     &           +umask(ji-1,jj,jk-1)+umask(ji  ,jj,jk), 1. )
285                 
286                  zmkv = 1./MAX(  vmask(ji,jj  ,jk-1)+vmask(ji,jj-1,jk)   &
287                     &           +vmask(ji,jj-1,jk-1)+vmask(ji,jj  ,jk), 1. )
288                 
289                  zcof3 = - e2t(ji,jj) * wslpi (ji,jj,jk) * zmku
290                  zcof4 = - e1t(ji,jj) * wslpj (ji,jj,jk) * zmkv
291                 
292                  zftw(ji,jk) = tmask(ji,jj,jk) *   &
293                     &    (  zcof0 * ( pt   (ji,jj,jk-1,jn) - pt   (ji  ,jj,jk,jn) )   &
294                     &     + zcof3 * ( zdit (ji   ,jk-1   ) + zdit (ji-1,jk      )     &
295                     &                +zdit (ji-1 ,jk-1   ) + zdit (ji  ,jk      ) )   &
296                     &     + zcof4 * ( zdjt (ji   ,jk-1   ) + zdj1t(ji  ,jk)     &
297                     &                +zdj1t(ji   ,jk-1   ) + zdjt (ji  ,jk      ) )  )                 
298               END DO
299            END DO
300
301
302            ! II.3 Divergence of vertical fluxes added to the horizontal divergence
303            ! ---------------------------------------------------------------------
304           
305            IF( kaht == 1 ) THEN
306               ! multiply the laplacian by the eddy diffusivity coefficient
307               DO jk = 1, jpkm1
308                  DO ji = 2, jpim1
309                     ! eddy coef. divided by the volume element
310                     zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
311                     ! vertical divergence
312                     ztav = fsahtt(ji,jj,jk) * ( zftw(ji,jk) - zftw(ji,jk+1) )
313                     ! harmonic operator applied to (pt,ps) and multiply by aht
314                     plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr
315                  END DO
316               END DO
317            ELSEIF( kaht == 2 ) THEN
318               ! second call, no multiplication
319               DO jk = 1, jpkm1
320                  DO ji = 2, jpim1
321                     ! inverse of the volume element
322                     zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
323                     ! vertical divergence
324                     ztav = zftw(ji,jk) - zftw(ji,jk+1)
325                     ! harmonic operator applied to (pt,ps)
326                     plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr
327                  END DO
328               END DO
329            ELSE
330               IF(lwp) WRITE(numout,*) ' ldfght: kaht= 1 or 2, here =', kaht
331               IF(lwp) WRITE(numout,*) '         We stop'
332               STOP 'ldfght'
333            ENDIF
334            !                                             ! ===============
335         END DO                                           !   End of slab
336         !                                                ! ===============
337      END DO
338      !
339      IF( wrk_not_released(2, 1,2,3) .OR. wrk_not_released_xz(1,2,3,4) )THEN
340         CALL ctl_stop('ldfght : failed to release workspace arrays.')
341      END IF
342      !
343   END SUBROUTINE ldfght
344
345#else 
346   !!----------------------------------------------------------------------
347   !!   Dummy module :             NO rotation of the lateral mixing tensor
348   !!----------------------------------------------------------------------
349CONTAINS
350   SUBROUTINE tra_ldf_bilapg( kt, cdtype, ptb, pta, kjpt )      ! Empty routine
351      CHARACTER(len=3) ::   cdtype
352      REAL, DIMENSION(:,:,:,:) ::   ptb, pta
353      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
354   END SUBROUTINE tra_ldf_bilapg
355#endif
356
357   !!==============================================================================
358END MODULE traldf_bilapg
Note: See TracBrowser for help on using the repository browser.