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

source: branches/UKMO/dev_r5518_optim_GO6_alloc/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90 @ 7601

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

#1821 more changes to replace pointers with ALLOCATE in traldf_iso
File size: 19.2 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   !!----------------------------------------------------------------------
12#if   defined key_ldfslp   ||   defined key_esopa
13   !!----------------------------------------------------------------------
14   !!   'key_ldfslp'               slope of the lateral diffusive direction
15   !!----------------------------------------------------------------------
16   !!   tra_ldf_iso  : update the tracer trend with the horizontal
17   !!                  component of a iso-neutral laplacian operator
18   !!                  and with the vertical part of
19   !!                  the isopycnal or geopotential s-coord. operator
20   !!----------------------------------------------------------------------
21   USE oce             ! ocean dynamics and active tracers
22   USE dom_oce         ! ocean space and time domain
23   USE trc_oce         ! share passive tracers/Ocean variables
24   USE zdf_oce         ! ocean vertical physics
25   USE ldftra_oce      ! ocean active tracers: lateral physics
26   USE ldfslp          ! iso-neutral slopes
27   USE diaptr          ! poleward transport diagnostics
28   USE trd_oce         ! trends: ocean variables
29   USE trdtra          ! trends manager: tracers
30   USE in_out_manager  ! I/O manager
31   USE iom             ! I/O library
32   USE phycst          ! physical constants
33   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
34   USE timing          ! Timing
35
36   IMPLICIT NONE
37   PRIVATE
38
39   PUBLIC   tra_ldf_iso   ! routine called by step.F90
40
41   !! * Substitutions
42#  include "domzgr_substitute.h90"
43#  include "ldftra_substitute.h90"
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
47   !! $Id$
48   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50CONTAINS
51
52   SUBROUTINE tra_ldf_iso( kt, kit000, cdtype, pgu, pgv,              &
53      &                                pgui, pgvi,                    &
54      &                                ptb, pta, kjpt, pahtb0 )
55      !!----------------------------------------------------------------------
56      !!                  ***  ROUTINE tra_ldf_iso  ***
57      !!
58      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
59      !!      trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and
60      !!      add it to the general trend of tracer equation.
61      !!
62      !! ** Method  :   The horizontal component of the lateral diffusive trends
63      !!      is provided by a 2nd order operator rotated along neural or geopo-
64      !!      tential surfaces to which an eddy induced advection can be added
65      !!      It is computed using before fields (forward in time) and isopyc-
66      !!      nal or geopotential slopes computed in routine ldfslp.
67      !!
68      !!      1st part :  masked horizontal derivative of T  ( di[ t ] )
69      !!      ========    with partial cell update if ln_zps=T.
70      !!
71      !!      2nd part :  horizontal fluxes of the lateral mixing operator
72      !!      ========   
73      !!         zftu = (aht+ahtb0) e2u*e3u/e1u di[ tb ]
74      !!               - aht       e2u*uslp    dk[ mi(mk(tb)) ]
75      !!         zftv = (aht+ahtb0) e1v*e3v/e2v dj[ tb ]
76      !!               - aht       e2u*vslp    dk[ mj(mk(tb)) ]
77      !!      take the horizontal divergence of the fluxes:
78      !!         difft = 1/(e1t*e2t*e3t) {  di-1[ zftu ] +  dj-1[ zftv ]  }
79      !!      Add this trend to the general trend (ta,sa):
80      !!         ta = ta + difft
81      !!
82      !!      3rd part: vertical trends of the lateral mixing operator
83      !!      ========  (excluding the vertical flux proportional to dk[t] )
84      !!      vertical fluxes associated with the rotated lateral mixing:
85      !!         zftw =-aht {  e2t*wslpi di[ mi(mk(tb)) ]
86      !!                     + e1t*wslpj dj[ mj(mk(tb)) ]  }
87      !!      take the horizontal divergence of the fluxes:
88      !!         difft = 1/(e1t*e2t*e3t) dk[ zftw ]
89      !!      Add this trend to the general trend (ta,sa):
90      !!         pta = pta + difft
91      !!
92      !! ** Action :   Update pta arrays with the before rotated diffusion
93      !!----------------------------------------------------------------------
94      USE oce     , ONLY:   zftu => ua       , zftv  => va         ! (ua,va) used as workspace
95      !
96      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
97      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
98      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
99      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
100      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu , pgv    ! tracer gradient at pstep levels
101      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgui, pgvi   ! tracer gradient at pstep levels
102      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields
103      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
104      REAL(wp)                             , INTENT(in   ) ::   pahtb0     ! background diffusion coef
105      !
106      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
107      INTEGER  ::  ikt
108      REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3       ! local scalars
109      REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4       !   -      -
110      REAL(wp) ::  zcoef0, zbtr                      !   -      -
111      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) ::  z2d
112      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  zdkt, zdk1t, zdit, zdjt, ztfw 
113      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  ztrax, ztray, ztraz 
114      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::  ztrax_T, ztray_T, ztraz_T
115      !!----------------------------------------------------------------------
116      !
117      IF( nn_timing == 1 )  CALL timing_start('tra_ldf_iso')
118      !
119      ALLOCATE( z2d(jpi, jpj)) 
120      ALLOCATE( zdit(jpi, jpj, jpk))
121      ALLOCATE( zdjt(jpi, jpj, jpk)) 
122      ALLOCATE( ztfw(jpi, jpj, jpk)) 
123      ALLOCTAE( zdkt(jpi, jpj, jpk)) 
124      ALLOCATE( zdk1t(jpi, jpj, jpk)) 
125      ALLOCATE( ztrax(jpi,jpj,jpk), ztray(jpi,jpj,jpk), ztraz(jpi,jpj,jpk) ) 
126      IF( l_trdtra .and. cdtype == 'TRA' ) ALLOCATE( ztrax_T(jpi,jpj,jpk), ztray_T(jpi,jpj,jpk), ztraz_T(jpi,jpj,jpk) ) 
127      !
128
129      IF( kt == kit000 )  THEN
130         IF(lwp) WRITE(numout,*)
131         IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator on ', cdtype
132         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
133      ENDIF
134      !
135      !                                                          ! ===========
136      DO jn = 1, kjpt                                            ! tracer loop
137         !                                                       ! ===========
138         ztrax(:,:,:) = 0._wp ; ztray(:,:,:) = 0._wp ; ztraz(:,:,:) = 0._wp ; 
139         !                                               
140         !!----------------------------------------------------------------------
141         !!   I - masked horizontal derivative
142         !!----------------------------------------------------------------------
143         !!bug ajout.... why?   ( 1,jpj,:) and (jpi,1,:) should be sufficient....
144         zdit (1,:,:) = 0.e0     ;     zdit (jpi,:,:) = 0.e0
145         zdjt (1,:,:) = 0.e0     ;     zdjt (jpi,:,:) = 0.e0
146         !!end
147
148         ! Horizontal tracer gradient
149         DO jk = 1, jpkm1
150            DO jj = 1, jpjm1
151               DO ji = 1, fs_jpim1   ! vector opt.
152                  zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk)
153                  zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
154               END DO
155            END DO
156         END DO
157
158         ! partial cell correction
159         IF( ln_zps ) THEN      ! partial steps correction at the last ocean level
160            DO jj = 1, jpjm1
161               DO ji = 1, fs_jpim1   ! vector opt.
162! IF useless if zpshde defines pgu everywhere
163                  zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)         
164                  zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
165               END DO
166            END DO
167         ENDIF
168         IF( ln_zps .AND. ln_isfcav ) THEN      ! partial steps correction at the first wet level beneath a cavity
169            DO jj = 1, jpjm1
170               DO ji = 1, fs_jpim1   ! vector opt.
171                  IF (miku(ji,jj) > 1) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)         
172                  IF (mikv(ji,jj) > 1) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)     
173               END DO
174            END DO
175         END IF
176
177         !!----------------------------------------------------------------------
178         !!   II - horizontal trend  (full)
179         !!----------------------------------------------------------------------
180!!!!!!!!!!CDIR PARALLEL DO PRIVATE( zdk1t )
181            ! 1. Vertical tracer gradient at level jk and jk+1
182            ! ------------------------------------------------
183         !
184         ! interior value
185         DO jk = 2, jpkm1               
186            DO jj = 1, jpj
187               DO ji = 1, jpi   ! vector opt.
188                  zdk1t(ji,jj,jk) = ( ptb(ji,jj,jk,jn  ) - ptb(ji,jj,jk+1,jn) ) * wmask(ji,jj,jk+1)
189                  !
190                  zdkt(ji,jj,jk)  = ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn  ) ) * wmask(ji,jj,jk)
191               END DO
192            END DO
193         END DO
194         ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
195         zdk1t(:,:,1) = ( ptb(:,:,1,jn  ) - ptb(:,:,2,jn) ) * wmask(:,:,2)
196         zdkt (:,:,1) = zdk1t(:,:,1)
197         IF ( ln_isfcav ) THEN
198            DO jj = 1, jpj
199               DO ji = 1, jpi   ! vector opt.
200                  ikt = mikt(ji,jj) ! surface level
201                  zdk1t(ji,jj,ikt) = ( ptb(ji,jj,ikt,jn  ) - ptb(ji,jj,ikt+1,jn) ) * wmask(ji,jj,ikt+1)
202                  zdkt (ji,jj,ikt) = zdk1t(ji,jj,ikt)
203               END DO
204            END DO
205         END IF
206
207         ! 2. Horizontal fluxes
208         ! --------------------   
209         DO jk = 1, jpkm1
210            DO jj = 1 , jpjm1
211               DO ji = 1, fs_jpim1   ! vector opt.
212                  zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk)
213                  zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk)
214                  !
215                  zmsku = 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                  !
218                  zmskv = 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 = - fsahtu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku
222                  zcof2 = - fsahtv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv
223                  !
224                  zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)   &
225                     &              + zcof1 * (  zdkt (ji+1,jj,jk) + zdk1t(ji,jj,jk)      &
226                     &                         + zdk1t(ji+1,jj,jk) + zdkt (ji,jj,jk)  )  ) * umask(ji,jj,jk)
227                  zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   &
228                     &              + zcof2 * (  zdkt (ji,jj+1,jk) + zdk1t(ji,jj,jk)      &
229                     &                         + zdk1t(ji,jj+1,jk) + zdkt (ji,jj,jk)  )  ) * vmask(ji,jj,jk)                 
230               END DO
231            END DO
232
233            ! II.4 Second derivative (divergence) and add to the general trend
234            ! ----------------------------------------------------------------
235            DO jj = 2 , jpjm1
236               DO ji = fs_2, fs_jpim1   ! vector opt.
237                  zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) )
238                  ztrax(ji,jj,jk) = zbtr * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) )
239                  ztray(ji,jj,jk) = zbtr * ( zftv(ji,jj,jk) - zftv(ji,jj-1,jk) )
240               END DO
241            END DO
242            !                                          ! ===============
243         END DO                                        !   End of slab 
244         !                                             ! ===============
245         !
246         pta(:,:,:,jn) = pta(:,:,:,jn) + ztrax(:,:,:) + ztray(:,:,:)
247         !
248         ! "Poleward" diffusive heat or salt transports (T-S case only)
249            ! note sign is reversed to give down-gradient diffusive transports (#1043)
250         IF( cdtype == 'TRA' .AND. ln_diaptr ) CALL dia_ptr_ohst_components( jn, 'ldf', -zftv(:,:,:)  )
251 
252         IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN
253           !
254           IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN
255               z2d(:,:) = 0._wp 
256               DO jk = 1, jpkm1
257                  DO jj = 2, jpjm1
258                     DO ji = fs_2, fs_jpim1   ! vector opt.
259                        z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 
260                     END DO
261                  END DO
262               END DO
263               z2d(:,:) = - rau0_rcp * z2d(:,:)     ! note sign is reversed to give down-gradient diffusive transports (#1043)
264               CALL lbc_lnk( z2d, 'U', -1. )
265               CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction
266               !
267               z2d(:,:) = 0._wp 
268               DO jk = 1, jpkm1
269                  DO jj = 2, jpjm1
270                     DO ji = fs_2, fs_jpim1   ! vector opt.
271                        z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 
272                     END DO
273                  END DO
274               END DO
275               z2d(:,:) = - rau0_rcp * z2d(:,:)     ! note sign is reversed to give down-gradient diffusive transports (#1043)
276               CALL lbc_lnk( z2d, 'V', -1. )
277               CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in i-direction
278            END IF
279            !
280         ENDIF
281
282         !!----------------------------------------------------------------------
283         !!   III - vertical trend of T & S (extra diagonal terms only)
284         !!----------------------------------------------------------------------
285         
286         ! Local constant initialization
287         ! -----------------------------
288         ztfw(1,:,:) = 0.e0     ;     ztfw(jpi,:,:) = 0.e0
289         
290         ! Vertical fluxes
291         ! ---------------
292         
293         ! Surface and bottom vertical fluxes set to zero
294         ztfw(:,:, 1 ) = 0.e0      ;      ztfw(:,:,jpk) = 0.e0
295         
296         ! interior (2=<jk=<jpk-1)
297         DO jk = 2, jpkm1
298            DO jj = 2, jpjm1
299               DO ji = fs_2, fs_jpim1   ! vector opt.
300                  zcoef0 = - fsahtw(ji,jj,jk) * wmask(ji,jj,jk)
301                  !
302                  zmsku = 1./MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)      &
303                     &            + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk), 1.  )
304                  zmskv = 1./MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)      &
305                     &            + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk), 1.  )
306                  !
307                  zcoef3 = zcoef0 * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk)
308                  zcoef4 = zcoef0 * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk)
309                  !
310                  ztfw(ji,jj,jk) = zcoef3 * (   zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)      &
311                     &                        + zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk)  )   &
312                     &           + zcoef4 * (   zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk)      &
313                     &                        + zdjt(ji  ,jj-1,jk-1) + zdjt(ji  ,jj  ,jk)  )
314               END DO
315            END DO
316         END DO
317         
318         
319         ! I.5 Divergence of vertical fluxes added to the general tracer trend
320         ! -------------------------------------------------------------------
321         DO jk = 1, jpkm1
322            DO jj = 2, jpjm1
323               DO ji = fs_2, fs_jpim1   ! vector opt.
324                  zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) )
325                  ztraz(ji,jj,jk) = (  ztfw(ji,jj,jk) - ztfw(ji,jj,jk+1)  ) * zbtr
326               END DO
327            END DO
328         END DO
329         pta(:,:,:,jn) = pta(:,:,:,jn) + ztraz(:,:,:)
330         !
331         IF( l_trdtra .AND. cdtype == "TRA" .AND. jn .eq. 1 )  THEN      ! save the temperature trends
332            ztrax_T(:,:,:) = ztrax(:,:,:)
333            ztray_T(:,:,:) = ztray(:,:,:)
334            ztraz_T(:,:,:) = ztraz(:,:,:)
335         ENDIF
336         IF( l_trdtrc .AND. cdtype == "TRC" )   THEN      ! save the horizontal component of diffusive trends for further diagnostics
337            CALL trd_tra( kt, cdtype, jn, jptra_iso_x, ztrax )
338            CALL trd_tra( kt, cdtype, jn, jptra_iso_y, ztray ) 
339            CALL trd_tra( kt, cdtype, jn, jptra_iso_z1, ztraz )  ! This is the first part of the vertical component.
340         ENDIF
341      END DO
342      !
343      IF( l_trdtra .AND. cdtype == "TRA" )   THEN      ! save the horizontal component of diffusive trends for further diagnostics
344         CALL trd_tra( kt, cdtype, jp_tem, jptra_iso_x, ztrax_T )
345         CALL trd_tra( kt, cdtype, jp_sal, jptra_iso_x, ztrax )
346         CALL trd_tra( kt, cdtype, jp_tem, jptra_iso_y, ztray_T )
347         CALL trd_tra( kt, cdtype, jp_sal, jptra_iso_y, ztray )
348         CALL trd_tra( kt, cdtype, jp_tem, jptra_iso_z1, ztraz_T )  ! This is the first part of the vertical component
349         CALL trd_tra( kt, cdtype, jp_sal, jptra_iso_z1, ztraz )    !
350      ENDIF
351      !
352      DEALLOCATE( z2d ) 
353      DEALLOCATE( zdit) 
354      DEALLOCATE( zdjt)
355      DEALLOCATE( ztfw) 
356      DEALLOCATE( zdkt )
357      DEALLOCATE( zdk1t ) 
358      DEALLOCATE( ztrax, ztray, ztraz ) 
359      IF( l_trdtra  .and. cdtype == 'TRA' ) DEALLOCATE( ztrax_T, ztray_T, ztraz_T ) 
360      !
361      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_iso')
362      !
363   END SUBROUTINE tra_ldf_iso
364
365#else
366   !!----------------------------------------------------------------------
367   !!   default option :   Dummy code   NO rotation of the diffusive tensor
368   !!----------------------------------------------------------------------
369CONTAINS
370   SUBROUTINE tra_ldf_iso( kt, kit000,cdtype, pgu, pgv, pgui, pgvi, ptb, pta, kjpt, pahtb0 )      ! Empty routine
371      INTEGER:: kt, kit000
372      CHARACTER(len=3) ::   cdtype
373      REAL, DIMENSION(:,:,:) ::   pgu, pgv, pgui, pgvi    ! tracer gradient at pstep levels
374      REAL, DIMENSION(:,:,:,:) ::   ptb, pta
375      WRITE(*,*) 'tra_ldf_iso: You should not have seen this print! error?', kt, kit000, cdtype,   &
376         &                       pgu(1,1,1), pgv(1,1,1), ptb(1,1,1,1), pta(1,1,1,1), kjpt, pahtb0
377   END SUBROUTINE tra_ldf_iso
378#endif
379
380   !!==============================================================================
381END MODULE traldf_iso
Note: See TracBrowser for help on using the repository browser.