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_iso.F90 in NEMO/branches/UKMO/BPC_miniapp/Master – NEMO

source: NEMO/branches/UKMO/BPC_miniapp/Master/traldf_iso.F90 @ 10831

Last change on this file since 10831 was 10831, checked in by wayne_gaudin, 5 years ago

Ticket #2197 - added the Master directory with code and makefile

File size: 20.0 KB
Line 
1MODULE traldf_iso
2   !!======================================================================
3   !!                   ***  MODULE  traldf_iso  ***
4   !! Ocean  tracers:  horizontal component of the lateral tracer mixing trend
5   !!======================================================================
6   !! History :  OPA  ! 1994-08  (G. Madec, M. Imbard)
7   !!            8.0  ! 1997-05  (G. Madec)  split into traldf and trazdf
8   !!            NEMO ! 2002-08  (G. Madec)  Free form, F90
9   !!            1.0  ! 2005-11  (G. Madec)  merge traldf and trazdf :-)
10   !!            3.3  ! 2010-09  (C. Ethe, G. Madec) Merge TRA-TRC
11   !!            3.7  ! 2014-01  (G. Madec, S. Masson)  restructuration/simplification of aht/aeiv specification
12   !!             -   ! 2014-02  (F. Lemarie, G. Madec)  triad operator (Griffies) + Method of Stabilizing Correction
13   !!----------------------------------------------------------------------
14
15   !!----------------------------------------------------------------------
16   !!   tra_ldf_iso   : update the tracer trend with the horizontal component of a iso-neutral laplacian operator
17   !!                   and with the vertical part of the isopycnal or geopotential s-coord. operator
18   !!----------------------------------------------------------------------
19   USE par_kind
20   USE len_oce            ! ocean dynamics and active tracers
21   !
22   USE in_out_manager ! I/O manager
23
24   IMPLICIT NONE
25
26   PUBLIC   tra_ldf_iso   ! routine called by step.F90
27
28   !! * Substitutions
29#  include "vectopt_loop_substitute.h90"
30   !!----------------------------------------------------------------------
31   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
32   !! $Id$
33   !! Software governed by the CeCILL licence     (./LICENSE)
34   !!----------------------------------------------------------------------
35CONTAINS
36
37  SUBROUTINE tra_ldf_iso( ln_traldf_msc, ln_traldf_blp, ln_traldf_lap, ln_zps, ln_isfcav, kjpt, kpass, r2dt, &
38      &                   e1t, e1u, e1v, e2t, e2u, e2v, e1e2t, e1_e2v, e2_e1u, r1_e1e2t, mbku, mbkv, miku, mikv, &
39      &                   umask, vmask, wmask, e3t_n, e3u_n, e3v_n, e3w_n, wslpi, wslpj, uslp, vslp, &
40      &                   pahu, pahv, pgu, pgv, pgui, pgvi, ptb , ptbb, &
41      &                   akz, ah_wslp2, pta )
42
43      !!----------------------------------------------------------------------
44      !!                  ***  ROUTINE tra_ldf_iso  ***
45      !!
46      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
47      !!      trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and
48      !!      add it to the general trend of tracer equation.
49      !!
50      !! ** Method  :   The horizontal component of the lateral diffusive trends
51      !!      is provided by a 2nd order operator rotated along neural or geopo-
52      !!      tential surfaces to which an eddy induced advection can be added
53      !!      It is computed using before fields (forward in time) and isopyc-
54      !!      nal or geopotential slopes computed in routine ldfslp.
55      !!
56      !!      1st part :  masked horizontal derivative of T  ( di[ t ] )
57      !!      ========    with partial cell update if ln_zps=T
58      !!                  with top     cell update if ln_isfcav
59      !!
60      !!      2nd part :  horizontal fluxes of the lateral mixing operator
61      !!      ========   
62      !!         zftu =  pahu e2u*e3u/e1u di[ tb ]
63      !!               - pahu e2u*uslp    dk[ mi(mk(tb)) ]
64      !!         zftv =  pahv e1v*e3v/e2v dj[ tb ]
65      !!               - pahv e2u*vslp    dk[ mj(mk(tb)) ]
66      !!      take the horizontal divergence of the fluxes:
67      !!         difft = 1/(e1e2t*e3t) {  di-1[ zftu ] +  dj-1[ zftv ]  }
68      !!      Add this trend to the general trend (ta,sa):
69      !!         ta = ta + difft
70      !!
71      !!      3rd part: vertical trends of the lateral mixing operator
72      !!      ========  (excluding the vertical flux proportional to dk[t] )
73      !!      vertical fluxes associated with the rotated lateral mixing:
74      !!         zftw = - {  mi(mk(pahu)) * e2t*wslpi di[ mi(mk(tb)) ]
75      !!                   + mj(mk(pahv)) * e1t*wslpj dj[ mj(mk(tb)) ]  }
76      !!      take the horizontal divergence of the fluxes:
77      !!         difft = 1/(e1e2t*e3t) dk[ zftw ]
78      !!      Add this trend to the general trend (ta,sa):
79      !!         pta = pta + difft
80      !!
81      !! ** Action :   Update pta arrays with the before rotated diffusion
82      !!----------------------------------------------------------------------
83      LOGICAL                              , INTENT(in   ) ::   ln_traldf_msc, ln_traldf_blp, ln_traldf_lap, ln_zps, ln_isfcav
84      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers passed in
85      INTEGER                              , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
86      REAL(wp)                             , INTENT(in   ) ::   r2dt       ! 2 . dt
87      REAL(wp), DIMENSION(jpi,jpj)         , INTENT(in   ) ::   e1t, e1u, e1v, e2t, e2u, e2v, e1e2t, e1_e2v, e2_e1u, r1_e1e2t     
88      INTEGER, DIMENSION(jpi,jpj)          , INTENT(in   ) ::   mbku, mbkv, miku, mikv
89      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   umask, vmask, wmask
90      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   e3t_n, e3u_n, e3v_n, e3w_n, wslpi, wslpj, uslp, vslp
91
92      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
93      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
94      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
95      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! tracer (kpass=1) or laplacian of tracer (kpass=2)
96      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptbb       ! tracer (only used in kpass=2)
97
98      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(inout) ::   akz, ah_wslp2 ! mainly OUT; IN only for unused points on edges 
99      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
100
101      !
102      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
103      INTEGER  ::  ikt
104      INTEGER  ::  ierr             ! local integer
105      REAL(wp) ::  zmsku, zahu_w, zabe1, zcof1, zcoef3   ! local scalars
106      REAL(wp) ::  zmskv, zahv_w, zabe2, zcof2, zcoef4   !   -      -
107      REAL(wp) ::  zcoef0, ze3w_2, zsign, z2dt, z1_2dt   !   -      -
108      REAL(wp) ,ALLOCATABLE, DIMENSION(:,:)     ::   zdkt, zdk1t, z2d      ! Make these use global scratch at some point
109      REAL(wp) ,ALLOCATABLE, DIMENSION(:,:,:)   ::   zdit, zdjt, zftu, zftv, ztfw 
110
111      ALLOCATE(zdkt(jpi,jpj),zdk1t(jpi,jpj),z2d(jpi,jpj))
112      ALLOCATE(zdit(jpi,jpj,jpk),zdjt(jpi,jpj,jpk),zftu(jpi,jpj,jpk),zftv(jpi,jpj,jpk),ztfw(jpi,jpj,jpk))
113      !!----------------------------------------------------------------------
114      !
115      !   
116      !
117      !                                            ! set time step size (Euler/Leapfrog)
118!$ACC KERNELS
119      z1_2dt = 1._wp / r2dt   !  note this was dependent on kt and neuler before
120      !
121      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
122      ELSE                    ;   zsign = -1._wp
123      ENDIF
124         
125      !!----------------------------------------------------------------------
126      !!   0 - calculate  ah_wslp2 and akz
127      !!----------------------------------------------------------------------
128      !
129      IF( kpass == 1 ) THEN                  !==  first pass only  ==!
130         !
131         DO jk = 2, jpkm1
132            DO jj = 2, jpjm1
133               DO ji = fs_2, fs_jpim1   ! vector opt.
134                  !
135                  zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          &
136                     &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  )
137                  zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          &
138                     &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  )
139                     !
140                  zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    &
141                     &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku
142                  zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    &
143                     &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv
144                     !
145                  ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
146                     &               + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk)
147               END DO
148            END DO
149         END DO
150         !
151         IF( ln_traldf_msc ) THEN                ! stabilizing vertical diffusivity coefficient
152            DO jk = 2, jpkm1
153               DO jj = 2, jpjm1
154                  DO ji = fs_2, fs_jpim1
155                     akz(ji,jj,jk) = 0.25_wp * (                                                                     &
156                        &              ( pahu(ji  ,jj,jk) + pahu(ji  ,jj,jk-1) ) / ( e1u(ji  ,jj) * e1u(ji  ,jj) )   &
157                        &            + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) )   &
158                        &            + ( pahv(ji,jj  ,jk) + pahv(ji,jj  ,jk-1) ) / ( e2v(ji,jj  ) * e2v(ji,jj  ) )   &
159                        &            + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) )   )
160                  END DO
161               END DO
162            END DO
163            !
164            IF( ln_traldf_blp ) THEN                ! bilaplacian operator
165               DO jk = 2, jpkm1
166                  DO jj = 1, jpjm1
167                     DO ji = 1, fs_jpim1
168                        akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   &
169                           &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) )  )
170                     END DO
171                  END DO
172               END DO
173            ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator
174               DO jk = 2, jpkm1
175                  DO jj = 1, jpjm1
176                     DO ji = 1, fs_jpim1
177                        ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk)
178                        zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
179                        akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt
180                     END DO
181                  END DO
182               END DO
183           ENDIF
184           !
185         ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit
186            akz(:,:,:) = ah_wslp2(:,:,:)     
187         ENDIF
188      ENDIF
189      !
190      !                                                          ! ===========
191      DO jn = 1, kjpt                                            ! tracer loop
192         !                                                       ! ===========
193         !                                               
194         !!----------------------------------------------------------------------
195         !!   I - masked horizontal derivative
196         !!----------------------------------------------------------------------
197!!gm : bug.... why (x,:,:)?   (1,jpj,:) and (jpi,1,:) should be sufficient....
198         zdit (1,:,:) = 0._wp     ;     zdit (jpi,:,:) = 0._wp
199         zdjt (1,:,:) = 0._wp     ;     zdjt (jpi,:,:) = 0._wp
200         !!end
201
202         ! Horizontal tracer gradient
203         DO jk = 1, jpkm1
204            DO jj = 1, jpjm1
205               DO ji = 1, fs_jpim1   ! vector opt.
206                  zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk)
207                  zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
208               END DO
209            END DO
210         END DO
211         IF( ln_zps ) THEN      ! botton and surface ocean correction of the horizontal gradient
212            DO jj = 1, jpjm1              ! bottom correction (partial bottom cell)
213               DO ji = 1, fs_jpim1   ! vector opt.
214                  zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)         
215                  zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
216               END DO
217            END DO
218            IF( ln_isfcav ) THEN      ! first wet level beneath a cavity
219               DO jj = 1, jpjm1
220                  DO ji = 1, fs_jpim1   ! vector opt.
221                     IF( miku(ji,jj) > 1 )   zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)         
222                     IF( mikv(ji,jj) > 1 )   zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)     
223                  END DO
224               END DO
225            ENDIF
226         ENDIF
227         !
228         !!----------------------------------------------------------------------
229         !!   II - horizontal trend  (full)
230         !!----------------------------------------------------------------------
231         !
232         DO jk = 1, jpkm1                                 ! Horizontal slab
233            !
234            !                             !== Vertical tracer gradient
235            zdk1t(:,:) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * wmask(:,:,jk+1)     ! level jk+1
236            !
237            IF( jk == 1 ) THEN   ;   zdkt(:,:) = zdk1t(:,:)                          ! surface: zdkt(jk=1)=zdkt(jk=2)
238            ELSE                 ;   zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk)
239            ENDIF
240            DO jj = 1 , jpjm1            !==  Horizontal fluxes
241               DO ji = 1, fs_jpim1   ! vector opt.
242                  zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk)
243                  zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk)
244                  !
245                  zmsku = 1._wp / MAX(  wmask(ji+1,jj,jk  ) + wmask(ji,jj,jk+1)   &
246                     &             + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk  ), 1._wp )
247                  !
248                  zmskv = 1._wp / MAX(  wmask(ji,jj+1,jk  ) + wmask(ji,jj,jk+1)   &
249                     &             + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk  ), 1._wp )
250                  !
251                  zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku
252                  zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv
253                  !
254                  zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)   &
255                     &               + zcof1 * (  zdkt (ji+1,jj) + zdk1t(ji,jj)      &
256                     &                          + zdk1t(ji+1,jj) + zdkt (ji,jj)  )  ) * umask(ji,jj,jk)
257                  zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   &
258                     &               + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj)      &
259                     &                          + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk)                 
260               END DO
261            END DO
262            !
263            DO jj = 2 , jpjm1          !== horizontal divergence and add to pta
264               DO ji = fs_2, fs_jpim1   ! vector opt.
265                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      &
266                     &                                           + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   &
267                     &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
268               END DO
269            END DO
270         END DO                                        !   End of slab 
271
272         !!----------------------------------------------------------------------
273         !!   III - vertical trend (full)
274         !!----------------------------------------------------------------------
275         !
276         ztfw(1,:,:) = 0._wp     ;     ztfw(jpi,:,:) = 0._wp
277         !
278         ! Vertical fluxes
279         ! ---------------
280         !                          ! Surface and bottom vertical fluxes set to zero
281         ztfw(:,:, 1 ) = 0._wp      ;      ztfw(:,:,jpk) = 0._wp
282         
283         DO jk = 2, jpkm1           ! interior (2=<jk=<jpk-1)
284            DO jj = 2, jpjm1
285               DO ji = fs_2, fs_jpim1   ! vector opt.
286                  !
287                  zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          &
288                     &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  )
289                  zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          &
290                     &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  )
291                     !
292                  zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    &
293                     &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku
294                  zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    &
295                     &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv
296                     !
297                  zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk)   !wslpi & j are already w-masked
298                  zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk)
299                  !
300                  ztfw(ji,jj,jk) = zcoef3 * (   zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)      &
301                     &                        + zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk)  )   &
302                     &           + zcoef4 * (   zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk)      &
303                     &                        + zdjt(ji  ,jj-1,jk-1) + zdjt(ji  ,jj  ,jk)  )
304               END DO
305            END DO
306         END DO
307         !                                !==  add the vertical 33 flux  ==!
308         IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz
309            DO jk = 2, jpkm1       
310               DO jj = 1, jpjm1
311                  DO ji = fs_2, fs_jpim1
312                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)   &
313                        &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             &
314                        &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
315                  END DO
316               END DO
317            END DO
318            !
319         ELSE                                   ! bilaplacian
320            SELECT CASE( kpass )
321            CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2
322               DO jk = 2, jpkm1 
323                  DO jj = 1, jpjm1
324                     DO ji = fs_2, fs_jpim1
325                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk)    &
326                           &           + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   &
327                           &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)
328                     END DO
329                  END DO
330               END DO
331            CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb  and ptbb gradients, resp.
332               DO jk = 2, jpkm1 
333                  DO jj = 1, jpjm1
334                     DO ji = fs_2, fs_jpim1
335                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)                      &
336                           &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   &
337                           &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   )
338                     END DO
339                  END DO
340               END DO
341            END SELECT
342         ENDIF
343         !         
344         DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pta  ==!
345            DO jj = 2, jpjm1
346               DO ji = fs_2, fs_jpim1   ! vector opt.
347                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   &
348                     &                                        * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
349               END DO
350            END DO
351         END DO
352         !
353         !                                                        ! ===============
354      END DO                                                      ! end tracer loop
355!$ACC END KERNELS
356      !
357   END SUBROUTINE tra_ldf_iso
358
359   !!==============================================================================
360END MODULE traldf_iso
Note: See TracBrowser for help on using the repository browser.