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

source: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90 @ 2830

Last change on this file since 2830 was 2830, checked in by kpedwards, 13 years ago

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