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

source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 5 years ago

GMED 450 add flush after prints

File size: 17.8 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         IF(lwp .AND. lflush) CALL flush(numout)
91      ENDIF
92
93      ! 1. Laplacian of ptb * aht
94      ! -----------------------------
95      CALL ldfght( kt, cdtype, ptb, zwk1, kjpt, 1 )      ! rotated harmonic operator applied to ptb and multiply by aht
96      !                                                 ! output in wk1
97      !
98      DO jn = 1, kjpt
99         CALL lbc_lnk( zwk1(:,:,:,jn) , 'T', 1. )        ! Lateral boundary conditions on wk1   (unchanged sign)
100      END DO
101
102      ! 2. Bilaplacian of ptb
103      ! -------------------------
104      CALL ldfght( kt, cdtype, zwk1, zwk2, kjpt, 2 )      ! rotated harmonic operator applied to wk1 ; output in wk2
105
106
107      ! 3. Update the tracer trends                    (j-slab :   2, jpj-1)
108      ! ---------------------------
109      DO jn = 1, kjpt
110         DO jj = 2, jpjm1
111            DO jk = 1, jpkm1
112               DO ji = 2, jpim1
113                  ! add it to the general tracer trends
114                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zwk2(ji,jj,jk,jn)
115               END DO
116            END DO
117         END DO
118      END DO
119      !
120      CALL wrk_dealloc( jpi, jpj, jpk, kjpt, zwk1, zwk2 ) 
121      !
122      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_bilapg')
123      !
124   END SUBROUTINE tra_ldf_bilapg
125
126
127   SUBROUTINE ldfght ( kt, cdtype, pt, plt, kjpt, kaht )
128      !!----------------------------------------------------------------------
129      !!                  ***  ROUTINE ldfght  ***
130      !!         
131      !! ** Purpose :   Apply a geopotential harmonic operator to (pt) and
132      !!      multiply it by the eddy diffusivity coefficient (if kaht=1).
133      !!      Routine only used in s-coordinates (l_sco=T) with bilaplacian
134      !!      operator (ln_traldf_bilap=T) acting along geopotential surfaces
135      !!      (ln_traldf_hor).
136      !!
137      !! ** Method  :   The harmonic operator rotated along geopotential
138      !!      surfaces is applied to (pt) using the slopes of geopotential
139      !!      surfaces computed in inildf routine. The result is provided in
140      !!      (plt,pls) arrays. It is computed in 2 steps:
141      !!
142      !!      First step: horizontal part of the operator. It is computed on
143      !!      ==========  pt as follows (idem on ps)
144      !!      horizontal fluxes :
145      !!         zftu = e2u*e3u/e1u di[ pt ] - e2u*uslp dk[ mi(mk(pt)) ]
146      !!         zftv = e1v*e3v/e2v dj[ pt ] - e1v*vslp dk[ mj(mk(pt)) ]
147      !!      take the horizontal divergence of the fluxes (no divided by
148      !!      the volume element :
149      !!         plt  = di-1[ zftu ] +  dj-1[ zftv ]
150      !!
151      !!      Second step: vertical part of the operator. It is computed on
152      !!      ===========  pt as follows (idem on ps)
153      !!      vertical fluxes :
154      !!         zftw = e1t*e2t/e3w * (wslpi^2+wslpj^2)  dk-1[ pt ]
155      !!              -     e2t     *       wslpi        di[ mi(mk(pt)) ]
156      !!              -     e1t     *       wslpj        dj[ mj(mk(pt)) ]
157      !!      take the vertical divergence of the fluxes add it to the hori-
158      !!      zontal component, divide the result by the volume element and
159      !!      if kaht=1, multiply by the eddy diffusivity coefficient:
160      !!         plt  = aht / (e1t*e2t*e3t) { plt + dk[ zftw ] }
161      !!      else:
162      !!         plt  =  1  / (e1t*e2t*e3t) { plt + dk[ zftw ] }
163      !!
164      !!----------------------------------------------------------------------
165      USE oce     , ONLY:   zftv => ua       ! ua used as workspace
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      REAL(wp), POINTER, DIMENSION(:,:) ::  zftu, zdkt, zdk1t
183      REAL(wp), POINTER, DIMENSION(:,:) ::  zftw, zdit, zdjt, zdj1t
184      !!----------------------------------------------------------------------
185      !
186      IF( nn_timing == 1 )  CALL timing_start('ldfght')
187      !
188      CALL wrk_alloc( jpi, jpj, zftu, zdkt, zdk1t ) 
189      CALL wrk_alloc( jpi, jpk, zftw, zdit, zdjt, zdj1t ) 
190      !
191      DO jn = 1, kjpt
192         !                               ! ********** !   ! ===============
193         DO jk = 1, jpkm1                ! First step !   ! Horizontal slab
194            !                            ! ********** !   ! ===============
195           
196            ! I.1 Vertical gradient of pt and ps at level jk and jk+1
197            ! -------------------------------------------------------
198            !     surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
199           
200            zdk1t(:,:) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
201            IF( jk == 1 ) THEN
202               zdkt(:,:) = zdk1t(:,:)
203            ELSE
204               zdkt(:,:) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk)
205            ENDIF
206
207
208            ! I.2 Horizontal fluxes
209            ! ---------------------
210           
211            DO jj = 1, jpjm1
212               DO ji = 1, jpim1
213                  zabe1 = re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk)
214                  zabe2 = re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk)
215                 
216                  zmku = 1./MAX( tmask(ji+1,jj,jk  )+tmask(ji,jj,jk+1)   &
217                     &          +tmask(ji+1,jj,jk+1)+tmask(ji,jj,jk  ),1. )
218                  zmkv = 1./MAX( tmask(ji,jj+1,jk  )+tmask(ji,jj,jk+1)   &
219                     &          +tmask(ji,jj+1,jk+1)+tmask(ji,jj,jk  ),1. )
220                 
221                  zcof1 = -e2u(ji,jj) * uslp(ji,jj,jk) * zmku
222                  zcof2 = -e1v(ji,jj) * vslp(ji,jj,jk) * zmkv
223                 
224                  zftu(ji,jj)= umask(ji,jj,jk) *   &
225                     &         (  zabe1 *( pt   (ji+1,jj,jk,jn) - pt(ji,jj,jk,jn) )   &
226                     &          + zcof1 *( zdkt (ji+1,jj) + zdk1t(ji,jj)     &
227                     &                    +zdk1t(ji+1,jj) + zdkt (ji,jj) )  )
228                 
229                  zftv(ji,jj,jk)= vmask(ji,jj,jk) *   &
230                     &         (  zabe2 *( pt(ji,jj+1,jk,jn) - pt(ji,jj,jk,jn) )   &
231                     &          + zcof2 *( zdkt (ji,jj+1) + zdk1t(ji,jj)     &
232                     &                    +zdk1t(ji,jj+1) + zdkt (ji,jj) )  )                 
233               END DO
234            END DO
235
236
237            ! I.3 Second derivative (divergence) (not divided by the volume)
238            ! ---------------------
239           
240            DO jj = 2 , jpjm1
241               DO ji = 2 , jpim1
242                  ztah = zftu(ji,jj) - zftu(ji-1,jj) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)
243                  plt(ji,jj,jk,jn) = ztah
244               END DO
245            END DO
246            !                                             ! ===============
247         END DO                                           !   End of slab
248         !                                                ! ===============
249         ! "Poleward" diffusive heat or salt transport
250        ! note sign is reversed to give down-gradient diffusive transports (#1043)
251         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( kaht == 2 ) ) CALL dia_ptr_ohst_components( jn, 'ldf', -zftv(:,:,:) )
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 = e12t(ji,jj) / fse3w_n(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 / ( e12t(ji,jj) * fse3t_n(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 / ( e12t(ji,jj) * fse3t_n(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               WRITE(numout,*) ' ldfght: kaht= 1 or 2, here =', kaht
331               WRITE(numout,*) '         We stop'
332               CALL ctl_stop( 'STOP', 'ldfght : unexpected kaht value')
333            ENDIF
334            !                                             ! ===============
335         END DO                                           !   End of slab
336         !                                                ! ===============
337      END DO
338      !
339      CALL wrk_dealloc( jpi, jpj, zftu, zdkt, zdk1t ) 
340      CALL wrk_dealloc( jpi, jpk, zftw, zdit, zdjt, zdj1t ) 
341      !
342      IF( nn_timing == 1 )  CALL timing_stop('ldfght')
343      !
344   END SUBROUTINE ldfght
345
346#else 
347   !!----------------------------------------------------------------------
348   !!   Dummy module :             NO rotation of the lateral mixing tensor
349   !!----------------------------------------------------------------------
350CONTAINS
351   SUBROUTINE tra_ldf_bilapg( kt, kit000, cdtype, ptb, pta, kjpt )      ! Empty routine
352      INTEGER :: kt, kit000
353      CHARACTER(len=3) ::   cdtype
354      REAL, DIMENSION(:,:,:,:) ::   ptb, pta
355      WRITE(*,*) 'tra_ldf_iso: You should not have seen this print! error?', &
356        &         kt, kit000, cdtype, ptb(1,1,1,1), pta(1,1,1,1), kjpt
357   END SUBROUTINE tra_ldf_bilapg
358#endif
359
360   !!==============================================================================
361END MODULE traldf_bilapg
Note: See TracBrowser for help on using the repository browser.