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

source: branches/dev_r1821_mixed_ldfdyn/NEMO/OPA_SRC/TRA/traldf_iso.F90 @ 1839

Last change on this file since 1839 was 1756, checked in by smasson, 14 years ago

implement AR5 diagnostics, see ticket:610
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 16.1 KB
Line 
1MODULE traldf_iso
2   !!======================================================================
3   !!                   ***  MODULE  traldf_iso  ***
4   !! Ocean active tracers:  horizontal component of the lateral tracer mixing trend
5   !!======================================================================
6   !! History :        !  94-08  (G. Madec, M. Imbard)
7   !!                  !  97-05  (G. Madec)  split into traldf and trazdf
8   !!             8.5  !  02-08  (G. Madec)  Free form, F90
9   !!             9.0  !  05-11  (G. Madec)  merge traldf and trazdf :-)
10   !!----------------------------------------------------------------------
11#if   defined key_ldfslp   ||   defined key_esopa
12   !!----------------------------------------------------------------------
13   !!   'key_ldfslp'               slope of the lateral diffusive direction
14   !!----------------------------------------------------------------------
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 ldftra_oce      ! ocean active tracers: lateral physics
24   USE trdmod          ! ocean active tracers trends
25   USE trdmod_oce      ! ocean variables trends
26   USE zdf_oce         ! ocean vertical physics
27   USE in_out_manager  ! I/O manager
28   USE iom             !
29   USE ldfslp          ! iso-neutral slopes
30   USE diaptr          ! poleward transport diagnostics
31   USE prtctl          ! Print control
32#if defined key_diaar5
33   USE phycst          ! physical constants
34   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
35#endif
36
37   IMPLICIT NONE
38   PRIVATE
39
40   PUBLIC   tra_ldf_iso   ! routine called by step.F90
41
42   !! * Substitutions
43#  include "domzgr_substitute.h90"
44#  include "ldftra_substitute.h90"
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !!   OPA 9.0 , LOCEAN-IPSL (2005)
48   !! $Id$
49   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
50   !!----------------------------------------------------------------------
51
52CONTAINS
53
54   SUBROUTINE tra_ldf_iso( kt )
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 & S ( 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      !!         ta = ta + difft
91      !!
92      !! ** Action :   Update (ta,sa) arrays with the before rotated diffusion
93      !!            trend (except the dk[ dk[.] ] term)
94      !!----------------------------------------------------------------------
95      USE oce           , zftv => ua   ! use ua as workspace
96      USE oce           , zfsv => va   ! use va as workspace
97      !!
98      INTEGER, INTENT( in ) ::   kt    ! ocean time-step index
99      !!
100      INTEGER  ::   ji, jj, jk   ! dummy loop indices
101      INTEGER  ::   iku, ikv     ! temporary integer
102      REAL(wp) ::   zmsku, zabe1, zcof1, zcoef3, zta   ! temporary scalars
103      REAL(wp) ::   zmskv, zabe2, zcof2, zcoef4, zsa   !    "         "
104      REAL(wp) ::   zcoef0, zbtr                       !    "         "
105      REAL(wp), DIMENSION(jpi,jpj)     ::   zdkt , zdk1t         ! 2D workspace
106      REAL(wp), DIMENSION(jpi,jpj)     ::   zdks , zdk1s, zfsu   !  "         "
107#if defined key_diaar5
108      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d                  !  "         "
109      REAL(wp)                         ::   zztmp                !  "         "
110      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zftu                 ! 3D workspace
111#else
112      REAL(wp), DIMENSION(jpi,jpj)     ::   zftu                 ! 2D workspace
113#endif
114      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdit, zdjt, ztfw     ! 3D workspace
115      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdis, zdjs, zsfw     !  "      "
116# if defined key_diaar5
117# endif 
118      !!----------------------------------------------------------------------
119
120      IF( kt == nit000 ) THEN
121         IF(lwp) WRITE(numout,*)
122         IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator'
123         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
124      ENDIF
125
126      !!----------------------------------------------------------------------
127      !!   I - masked horizontal derivative of T & S
128      !!----------------------------------------------------------------------
129!!bug ajout.... why?   ( 1,jpj,:) and (jpi,1,:) should be sufficient....
130      zdit (1,:,:) = 0.e0     ;     zdit (jpi,:,:) = 0.e0
131      zdis (1,:,:) = 0.e0     ;     zdis (jpi,:,:) = 0.e0
132      zdjt (1,:,:) = 0.e0     ;     zdjt (jpi,:,:) = 0.e0
133      zdjs (1,:,:) = 0.e0     ;     zdjs (jpi,:,:) = 0.e0
134!!end
135
136      ! Horizontal temperature and salinity gradient
137      DO jk = 1, jpkm1
138         DO jj = 1, jpjm1
139            DO ji = 1, fs_jpim1   ! vector opt.
140               zdit(ji,jj,jk) = ( tb(ji+1,jj  ,jk) - tb(ji,jj,jk) ) * umask(ji,jj,jk)
141               zdis(ji,jj,jk) = ( sb(ji+1,jj  ,jk) - sb(ji,jj,jk) ) * umask(ji,jj,jk)
142               zdjt(ji,jj,jk) = ( tb(ji  ,jj+1,jk) - tb(ji,jj,jk) ) * vmask(ji,jj,jk)
143               zdjs(ji,jj,jk) = ( sb(ji  ,jj+1,jk) - sb(ji,jj,jk) ) * vmask(ji,jj,jk)
144            END DO
145         END DO
146      END DO
147      IF( ln_zps ) THEN      ! partial steps correction at the last level
148         DO jj = 1, jpjm1
149            DO ji = 1, fs_jpim1   ! vector opt.
150               ! last level
151               iku = MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1
152               ikv = MIN( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1
153               zdit(ji,jj,iku) = gtu(ji,jj) 
154               zdis(ji,jj,iku) = gsu(ji,jj)               
155               zdjt(ji,jj,ikv) = gtv(ji,jj) 
156               zdjs(ji,jj,ikv) = gsv(ji,jj)               
157            END DO
158         END DO
159      ENDIF
160
161      !!----------------------------------------------------------------------
162      !!   II - horizontal trend of T & S (full)
163      !!----------------------------------------------------------------------
164     
165#if defined key_diaar5
166!CDIR PARALLEL DO PRIVATE( zdk1t, zdk1s, zfsu )
167#else
168!CDIR PARALLEL DO PRIVATE( zdk1t, zdk1s, zftu, zfsu )
169#endif
170      !                                                ! ===============
171      DO jk = 1, jpkm1                                 ! Horizontal slab
172         !                                             ! ===============
173         ! 1. Vertical tracer gradient at level jk and jk+1
174         ! ------------------------------------------------
175         ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
176
177         zdk1t(:,:) = ( tb(:,:,jk) - tb(:,:,jk+1) ) * tmask(:,:,jk+1)
178         zdk1s(:,:) = ( sb(:,:,jk) - sb(:,:,jk+1) ) * tmask(:,:,jk+1)
179
180         IF( jk == 1 ) THEN
181            zdkt(:,:) = zdk1t(:,:)
182            zdks(:,:) = zdk1s(:,:)
183         ELSE
184            zdkt(:,:) = ( tb(:,:,jk-1) - tb(:,:,jk) ) * tmask(:,:,jk)
185            zdks(:,:) = ( sb(:,:,jk-1) - sb(:,:,jk) ) * tmask(:,:,jk)
186         ENDIF
187
188
189         ! 2. Horizontal fluxes
190         ! --------------------
191
192         DO jj = 1 , jpjm1
193            DO ji = 1, fs_jpim1   ! vector opt.
194               zabe1 = ( fsahtu(ji,jj,jk) + ahtb0 ) * e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj)
195               zabe2 = ( fsahtv(ji,jj,jk) + ahtb0 ) * e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj)
196
197               zmsku = 1. / MAX(  tmask(ji+1,jj,jk  ) + tmask(ji,jj,jk+1)   &
198                  &             + tmask(ji+1,jj,jk+1) + tmask(ji,jj,jk  ), 1. )
199
200               zmskv = 1. / MAX(  tmask(ji,jj+1,jk  ) + tmask(ji,jj,jk+1)   &
201                  &             + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk  ), 1. )
202
203               zcof1 = -fsahtu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku
204               zcof2 = -fsahtv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv
205
206#if defined key_diaar5
207               zftu(ji,jj,jk) = (  zabe1 * zdit(ji,jj,jk)   &
208#else
209               zftu(ji,jj   ) = (  zabe1 * zdit(ji,jj,jk)   &
210#endif
211                  &              + zcof1 * (  zdkt (ji+1,jj) + zdk1t(ji,jj)      &
212                  &                         + zdk1t(ji+1,jj) + zdkt (ji,jj)  )  ) * umask(ji,jj,jk)
213               zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   &
214                  &              + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj)      &
215                  &                         + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk)
216               zfsu(ji,jj   ) = (  zabe1 * zdis(ji,jj,jk)   &
217                  &              + zcof1 * (  zdks (ji+1,jj) + zdk1s(ji,jj)      &
218                  &                         + zdk1s(ji+1,jj) + zdks (ji,jj)  )  ) * umask(ji,jj,jk)
219               zfsv(ji,jj,jk) = (  zabe2 * zdjs(ji,jj,jk)   &
220                  &              + zcof2 * (  zdks (ji,jj+1) + zdk1s(ji,jj)      &
221                  &                         + zdk1s(ji,jj+1) + zdks (ji,jj)  )  ) * vmask(ji,jj,jk)
222            END DO
223         END DO
224
225
226         ! II.4 Second derivative (divergence) and add to the general trend
227         ! ----------------------------------------------------------------
228         DO jj = 2 , jpjm1
229            DO ji = fs_2, fs_jpim1   ! vector opt.
230               zbtr= 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
231#if defined key_diaar5
232               zta = zbtr * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )
233#else
234               zta = zbtr * ( zftu(ji,jj   ) - zftu(ji-1,jj   ) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )
235#endif
236               zsa = zbtr * ( zfsu(ji,jj   ) - zfsu(ji-1,jj   ) + zfsv(ji,jj,jk) - zfsv(ji,jj-1,jk)  )
237               ta (ji,jj,jk) = ta (ji,jj,jk) + zta
238               sa (ji,jj,jk) = sa (ji,jj,jk) + zsa
239            END DO
240         END DO
241         !                                          ! ===============
242      END DO                                        !   End of slab 
243      !                                             ! ===============
244
245      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN   ! Poleward diffusive heat and salt transports
246         pht_ldf(:) = ptr_vj( zftv(:,:,:) )
247         pst_ldf(:) = ptr_vj( zfsv(:,:,:) )
248      ENDIF
249#if defined key_diaar5
250      zztmp = 0.5 * rau0 * rcp 
251      z2d(:,:) = 0.e0 
252      DO jk = 1, jpkm1
253         DO jj = 2, jpjm1
254            DO ji = fs_2, fs_jpim1   ! vector opt.
255               z2d(ji,jj) = z2d(ji,jj) + zztmp * zftu(ji,jj,jk) * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) ) * e1u(ji,jj) * fse3u(ji,jj,jk) 
256            END DO
257         END DO
258      END DO
259      CALL lbc_lnk( z2d, 'U', -1. )
260      CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction
261      z2d(:,:) = 0.e0 
262      DO jk = 1, jpkm1
263         DO jj = 2, jpjm1
264            DO ji = fs_2, fs_jpim1   ! vector opt.
265               z2d(ji,jj) = z2d(ji,jj) + zztmp * zftv(ji,jj,jk) * ( tn(ji,jj,jk) + tn(ji,jj+1,jk) ) * e2v(ji,jj) * fse3v(ji,jj,jk) 
266            END DO
267         END DO
268      END DO
269      CALL lbc_lnk( z2d, 'V', -1. )
270      CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in i-direction
271#endif
272
273      !!----------------------------------------------------------------------
274      !!   III - vertical trend of T & S (extra diagonal terms only)
275      !!----------------------------------------------------------------------
276
277      ! Local constant initialization
278      ! -----------------------------
279      ztfw(1,:,:) = 0.e0     ;     ztfw(jpi,:,:) = 0.e0
280      zsfw(1,:,:) = 0.e0     ;     zsfw(jpi,:,:) = 0.e0
281
282
283      ! Vertical fluxes
284      ! ---------------
285
286      ! Surface and bottom vertical fluxes set to zero
287      ztfw(:,:, 1 ) = 0.e0      ;      ztfw(:,:,jpk) = 0.e0
288      zsfw(:,:, 1 ) = 0.e0      ;      zsfw(:,:,jpk) = 0.e0
289
290      ! interior (2=<jk=<jpk-1)
291      DO jk = 2, jpkm1
292         DO jj = 2, jpjm1
293            DO ji = fs_2, fs_jpim1   ! vector opt.
294               zcoef0 = - fsahtw(ji,jj,jk) * tmask(ji,jj,jk)
295
296               zmsku = 1./MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)      &
297                  &            + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk), 1.  )
298
299               zmskv = 1./MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)      &
300                  &            + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk), 1.  )
301
302               zcoef3 = zcoef0 * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk)
303               zcoef4 = zcoef0 * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk)
304
305               ztfw(ji,jj,jk) = zcoef3 * (   zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)      &
306                  &                        + zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk)  )   &
307                  &           + zcoef4 * (   zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk)      &
308                  &                        + zdjt(ji  ,jj-1,jk-1) + zdjt(ji  ,jj  ,jk)  )
309
310               zsfw(ji,jj,jk) = zcoef3 * (   zdis(ji  ,jj  ,jk-1) + zdis(ji-1,jj  ,jk)      &
311                  &                        + zdis(ji-1,jj  ,jk-1) + zdis(ji  ,jj  ,jk)  )   &
312                  &           + zcoef4 * (   zdjs(ji  ,jj  ,jk-1) + zdjs(ji  ,jj-1,jk)      &
313                  &                        + zdjs(ji  ,jj-1,jk-1) + zdjs(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
322      DO jk = 1, jpkm1
323         DO jj = 2, jpjm1
324            DO ji = fs_2, fs_jpim1   ! vector opt.
325               zbtr =  1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
326               zta  = (  ztfw(ji,jj,jk) - ztfw(ji,jj,jk+1)  ) * zbtr
327               zsa  = (  zsfw(ji,jj,jk) - zsfw(ji,jj,jk+1)  ) * zbtr
328               ta(ji,jj,jk) = ta(ji,jj,jk) + zta
329               sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
330            END DO
331         END DO
332      END DO
333      !
334   END SUBROUTINE tra_ldf_iso
335
336#else
337   !!----------------------------------------------------------------------
338   !!   default option :   Dummy code   NO rotation of the diffusive tensor
339   !!----------------------------------------------------------------------
340CONTAINS
341   SUBROUTINE tra_ldf_iso( kt )               ! Empty routine
342      WRITE(*,*) 'tra_ldf_iso: You should not have seen this print! error?', kt
343   END SUBROUTINE tra_ldf_iso
344#endif
345
346   !!==============================================================================
347END MODULE traldf_iso
Note: See TracBrowser for help on using the repository browser.