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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90 @ 4401

Last change on this file since 4401 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

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