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

source: branches/UKMO/dev_r5518_GO6M_dev/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90 @ 8200

Last change on this file since 8200 was 8200, checked in by frrh, 7 years ago

Merge branches/UKMO/dev_r5518_GC3_couple_pkg@7985 using command:

svn merge -r 6574:7985 ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/UKMO/dev_r5518_GC3_couple_pkg
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   USE wrk_nemo        ! Memory Allocation
28   USE timing          ! Timing
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   tra_ldf_bilapg   ! routine called by step.F90
34
35   !! * Substitutions
36#  include "domzgr_substitute.h90"
37#  include "ldftra_substitute.h90"
38#  include "ldfeiv_substitute.h90"
39   !!----------------------------------------------------------------------
40   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
41   !! $Id$
42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE tra_ldf_bilapg( kt, kit000, cdtype, ptb, pta, kjpt )
47      !!----------------------------------------------------------------------
48      !!                 ***  ROUTINE tra_ldf_bilapg  ***
49      !!                   
50      !! ** Purpose :   Compute the before horizontal tracer diffusive
51      !!      trend and add it to the general trend of tracer equation.
52      !!
53      !! ** Method  :   The lateral diffusive trends is provided by a 4th order
54      !!      operator rotated along geopotential surfaces. It is computed
55      !!      using before fields (forward in time) and geopotential slopes
56      !!      computed in routine inildf.
57      !!         -1- compute the geopotential harmonic operator applied to
58      !!        ptb and multiply it by the eddy diffusivity coefficient
59      !!       (done by a call to ldfght routine, result in wk1 arrays).
60      !!      Applied the domain lateral boundary conditions by call to lbc_lnk
61      !!         -2- compute the geopotential harmonic operator applied to
62      !!      wk1 by a second call to ldfght routine (result in wk2)
63      !!      arrays).
64      !!         -3- Add this trend to the general trend
65      !!            pta = pta + wk2
66      !!
67      !! ** Action : - Update pta arrays with the before geopotential
68      !!               biharmonic mixing trend.
69      !!----------------------------------------------------------------------
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      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zwk1, zwk2 
80      !!----------------------------------------------------------------------
81      !
82      IF( nn_timing == 1 )  CALL timing_start('tra_ldf_bilapg')
83      !
84      CALL wrk_alloc( jpi, jpj, jpk, kjpt, zwk1, zwk2 ) 
85      !
86      IF( kt == kit000 )  THEN
87         IF(lwp) WRITE(numout,*)
88         IF(lwp) WRITE(numout,*) 'tra_ldf_bilapg : horizontal biharmonic operator in s-coordinate on ', cdtype
89         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
90      ENDIF
91
92      ! 1. Laplacian of ptb * aht
93      ! -----------------------------
94      CALL ldfght( kt, cdtype, ptb, zwk1, kjpt, 1 )      ! rotated harmonic operator applied to ptb and multiply by aht
95      !                                                 ! output in wk1
96      !
97      DO jn = 1, kjpt
98         CALL lbc_lnk( zwk1(:,:,:,jn) , 'T', 1. )        ! Lateral boundary conditions on wk1   (unchanged sign)
99      END DO
100
101      ! 2. Bilaplacian of ptb
102      ! -------------------------
103      CALL ldfght( kt, cdtype, zwk1, zwk2, kjpt, 2 )      ! rotated harmonic operator applied to wk1 ; output in wk2
104
105
106      ! 3. Update the tracer trends                    (j-slab :   2, jpj-1)
107      ! ---------------------------
108      DO jn = 1, kjpt
109         DO jj = 2, jpjm1
110            DO jk = 1, jpkm1
111               DO ji = 2, jpim1
112                  ! add it to the general tracer trends
113                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zwk2(ji,jj,jk,jn)
114               END DO
115            END DO
116         END DO
117      END DO
118      !
119      CALL wrk_dealloc( jpi, jpj, jpk, kjpt, zwk1, zwk2 ) 
120      !
121      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_bilapg')
122      !
123   END SUBROUTINE tra_ldf_bilapg
124
125
126   SUBROUTINE ldfght ( kt, cdtype, pt, plt, kjpt, kaht )
127      !!----------------------------------------------------------------------
128      !!                  ***  ROUTINE ldfght  ***
129      !!         
130      !! ** Purpose :   Apply a geopotential harmonic operator to (pt) and
131      !!      multiply it by the eddy diffusivity coefficient (if kaht=1).
132      !!      Routine only used in s-coordinates (l_sco=T) with bilaplacian
133      !!      operator (ln_traldf_bilap=T) acting along geopotential surfaces
134      !!      (ln_traldf_hor).
135      !!
136      !! ** Method  :   The harmonic operator rotated along geopotential
137      !!      surfaces is applied to (pt) using the slopes of geopotential
138      !!      surfaces computed in inildf routine. The result is provided in
139      !!      (plt,pls) arrays. It is computed in 2 steps:
140      !!
141      !!      First step: horizontal part of the operator. It is computed on
142      !!      ==========  pt as follows (idem on ps)
143      !!      horizontal fluxes :
144      !!         zftu = e2u*e3u/e1u di[ pt ] - e2u*uslp dk[ mi(mk(pt)) ]
145      !!         zftv = e1v*e3v/e2v dj[ pt ] - e1v*vslp dk[ mj(mk(pt)) ]
146      !!      take the horizontal divergence of the fluxes (no divided by
147      !!      the volume element :
148      !!         plt  = di-1[ zftu ] +  dj-1[ zftv ]
149      !!
150      !!      Second step: vertical part of the operator. It is computed on
151      !!      ===========  pt as follows (idem on ps)
152      !!      vertical fluxes :
153      !!         zftw = e1t*e2t/e3w * (wslpi^2+wslpj^2)  dk-1[ pt ]
154      !!              -     e2t     *       wslpi        di[ mi(mk(pt)) ]
155      !!              -     e1t     *       wslpj        dj[ mj(mk(pt)) ]
156      !!      take the vertical divergence of the fluxes add it to the hori-
157      !!      zontal component, divide the result by the volume element and
158      !!      if kaht=1, multiply by the eddy diffusivity coefficient:
159      !!         plt  = aht / (e1t*e2t*e3t) { plt + dk[ zftw ] }
160      !!      else:
161      !!         plt  =  1  / (e1t*e2t*e3t) { plt + dk[ zftw ] }
162      !!
163      !!----------------------------------------------------------------------
164      USE oce     , ONLY:   zftv => ua       ! ua used as workspace
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      REAL(wp), POINTER, DIMENSION(:,:) ::  zftu, zdkt, zdk1t
182      REAL(wp), POINTER, DIMENSION(:,:) ::  zftw, zdit, zdjt, zdj1t
183      !!----------------------------------------------------------------------
184      !
185      IF( nn_timing == 1 )  CALL timing_start('ldfght')
186      !
187      CALL wrk_alloc( jpi, jpj, zftu, zdkt, zdk1t ) 
188      CALL wrk_alloc( jpi, jpk, zftw, zdit, zdjt, zdj1t ) 
189      !
190      DO jn = 1, kjpt
191         !                               ! ********** !   ! ===============
192         DO jk = 1, jpkm1                ! First step !   ! Horizontal slab
193            !                            ! ********** !   ! ===============
194           
195            ! I.1 Vertical gradient of pt and ps at level jk and jk+1
196            ! -------------------------------------------------------
197            !     surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
198           
199            zdk1t(:,:) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
200            IF( jk == 1 ) THEN
201               zdkt(:,:) = zdk1t(:,:)
202            ELSE
203               zdkt(:,:) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk)
204            ENDIF
205
206
207            ! I.2 Horizontal fluxes
208            ! ---------------------
209           
210            DO jj = 1, jpjm1
211               DO ji = 1, jpim1
212                  zabe1 = re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk)
213                  zabe2 = re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk)
214                 
215                  zmku = 1./MAX( tmask(ji+1,jj,jk  )+tmask(ji,jj,jk+1)   &
216                     &          +tmask(ji+1,jj,jk+1)+tmask(ji,jj,jk  ),1. )
217                  zmkv = 1./MAX( tmask(ji,jj+1,jk  )+tmask(ji,jj,jk+1)   &
218                     &          +tmask(ji,jj+1,jk+1)+tmask(ji,jj,jk  ),1. )
219                 
220                  zcof1 = -e2u(ji,jj) * uslp(ji,jj,jk) * zmku
221                  zcof2 = -e1v(ji,jj) * vslp(ji,jj,jk) * zmkv
222                 
223                  zftu(ji,jj)= umask(ji,jj,jk) *   &
224                     &         (  zabe1 *( pt   (ji+1,jj,jk,jn) - pt(ji,jj,jk,jn) )   &
225                     &          + zcof1 *( zdkt (ji+1,jj) + zdk1t(ji,jj)     &
226                     &                    +zdk1t(ji+1,jj) + zdkt (ji,jj) )  )
227                 
228                  zftv(ji,jj,jk)= vmask(ji,jj,jk) *   &
229                     &         (  zabe2 *( pt(ji,jj+1,jk,jn) - pt(ji,jj,jk,jn) )   &
230                     &          + zcof2 *( zdkt (ji,jj+1) + zdk1t(ji,jj)     &
231                     &                    +zdk1t(ji,jj+1) + zdkt (ji,jj) )  )                 
232               END DO
233            END DO
234
235
236            ! I.3 Second derivative (divergence) (not divided by the volume)
237            ! ---------------------
238           
239            DO jj = 2 , jpjm1
240               DO ji = 2 , jpim1
241                  ztah = zftu(ji,jj) - zftu(ji-1,jj) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)
242                  plt(ji,jj,jk,jn) = ztah
243               END DO
244            END DO
245            !                                             ! ===============
246         END DO                                           !   End of slab
247         !                                                ! ===============
248         ! "Poleward" diffusive heat or salt transport
249        ! note sign is reversed to give down-gradient diffusive transports (#1043)
250         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( kaht == 2 ) ) CALL dia_ptr_ohst_components( jn, 'ldf', -zftv(:,:,:) )
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 = e12t(ji,jj) / fse3w_n(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 / ( e12t(ji,jj) * fse3t_n(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 / ( e12t(ji,jj) * fse3t_n(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               WRITE(numout,*) ' ldfght: kaht= 1 or 2, here =', kaht
330               WRITE(numout,*) '         We stop'
331               CALL ctl_stop( 'STOP', 'ldfght : unexpected kaht value')
332            ENDIF
333            !                                             ! ===============
334         END DO                                           !   End of slab
335         !                                                ! ===============
336      END DO
337      !
338      CALL wrk_dealloc( jpi, jpj, zftu, zdkt, zdk1t ) 
339      CALL wrk_dealloc( jpi, jpk, zftw, zdit, zdjt, zdj1t ) 
340      !
341      IF( nn_timing == 1 )  CALL timing_stop('ldfght')
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, kit000, cdtype, ptb, pta, kjpt )      ! Empty routine
351      INTEGER :: kt, kit000
352      CHARACTER(len=3) ::   cdtype
353      REAL, DIMENSION(:,:,:,:) ::   ptb, pta
354      WRITE(*,*) 'tra_ldf_iso: You should not have seen this print! error?', &
355        &         kt, kit000, cdtype, ptb(1,1,1,1), pta(1,1,1,1), kjpt
356   END SUBROUTINE tra_ldf_bilapg
357#endif
358
359   !!==============================================================================
360END MODULE traldf_bilapg
Note: See TracBrowser for help on using the repository browser.