ROMS
Loading...
Searching...
No Matches
Functions/Subroutines
tl_t3dmix4_mod Module Reference

Functions/Subroutines

subroutine, public tl_t3dmix4 (ng, tile)
 
subroutine tl_t3dmix4_geo_tile (ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nrhs, nstp, nnew, umask, vmask, umask_wet, vmask_wet, om_v, on_u, pm, pn, hz, tl_hz, z_r, tl_z_r, diff3d_u, diff3d_v, diff3d_r, diff4, tclm, t, tl_t)
 
subroutine tl_t3dmix4_iso_tile (ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nrhs, nstp, nnew, umask, vmask, umask_wet, vmask_wet, om_v, on_u, pm, pn, hz, tl_hz, z_r, tl_z_r, diff3d_u, diff3d_v, diff3d_r, diff4, pden, tl_pden, tclm, t, tl_t)
 
subroutine tl_t3dmix4_s_tile (ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nrhs, nstp, nnew, umask, vmask, umask_wet, vmask_wet, hz, tl_hz, pmon_u, pnom_v, pm, pn, diff3d_u, diff3d_v, diff3d_r, diff4, tclm, t, tl_t)
 

Function/Subroutine Documentation

◆ tl_t3dmix4()

subroutine public tl_t3dmix4_mod::tl_t3dmix4 ( integer, intent(in) ng,
integer, intent(in) tile )

Definition at line 25 of file tl_t3dmix4_geo.h.

26!***********************************************************************
27!
28 USE mod_param
29#ifdef TS_MIX_CLIMA
30 USE mod_clima
31#endif
32#ifdef DIAGNOSTICS_TS
33!! USE mod_diags
34#endif
35 USE mod_grid
36 USE mod_mixing
37 USE mod_ocean
38 USE mod_stepping
39!
40! Imported variable declarations.
41!
42 integer, intent(in) :: ng, tile
43!
44! Local variable declarations.
45!
46 character (len=*), parameter :: MyFile = &
47 & __FILE__
48!
49#include "tile.h"
50!
51#ifdef PROFILE
52 CALL wclock_on (ng, itlm, 28, __line__, myfile)
53#endif
54 CALL tl_t3dmix4_geo_tile (ng, tile, &
55 & lbi, ubi, lbj, ubj, &
56 & imins, imaxs, jmins, jmaxs, &
57 & nrhs(ng), nstp(ng), nnew(ng), &
58#ifdef MASKING
59 & grid(ng) % umask, &
60 & grid(ng) % vmask, &
61#endif
62#ifdef WET_DRY_NOT_YET
63 & grid(ng) % umask_wet, &
64 & grid(ng) % vmask_wet, &
65#endif
66 & grid(ng) % om_v, &
67 & grid(ng) % on_u, &
68 & grid(ng) % pm, &
69 & grid(ng) % pn, &
70 & grid(ng) % Hz, &
71 & grid(ng) % tl_Hz, &
72 & grid(ng) % z_r, &
73 & grid(ng) % tl_z_r, &
74#ifdef DIFF_3DCOEF
75# ifdef TS_U3ADV_SPLIT
76 & mixing(ng) % diff3d_u, &
77 & mixing(ng) % diff3d_v, &
78# else
79 & mixing(ng) % diff3d_r, &
80# endif
81#else
82 & mixing(ng) % diff4, &
83#endif
84#ifdef TS_MIX_CLIMA
85 & clima(ng) % tclm, &
86#endif
87#ifdef DIAGNOSTICS_TS
88!! & DIAGS(ng) % DiaTwrk, &
89#endif
90 & ocean(ng) % t, &
91 & ocean(ng) % tl_t)
92#ifdef PROFILE
93 CALL wclock_off (ng, itlm, 28, __line__, myfile)
94#endif
95!
96 RETURN
mod_clima
Definition mod_clima.F:2
mod_clima::clima
type(t_clima), dimension(:), allocatable clima
Definition mod_clima.F:153
mod_grid
Definition mod_grid.F:2
mod_grid::grid
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
mod_mixing
Definition mod_mixing.F:2
mod_mixing::mixing
type(t_mixing), dimension(:), allocatable mixing
Definition mod_mixing.F:399
mod_ocean
Definition mod_ocean.F:2
mod_ocean::ocean
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
mod_param
Definition mod_param.F:2
mod_param::itlm
integer, parameter itlm
Definition mod_param.F:663
mod_stepping
Definition mod_stepping.F:2
mod_stepping::nrhs
integer, dimension(:), allocatable nrhs
Definition mod_stepping.F:70
mod_stepping::nnew
integer, dimension(:), allocatable nnew
Definition mod_stepping.F:69
mod_stepping::nstp
integer, dimension(:), allocatable nstp
Definition mod_stepping.F:71
wclock_off
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
wclock_on
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3

References mod_clima::clima, mod_grid::grid, mod_param::itlm, mod_mixing::mixing, mod_stepping::nnew, mod_stepping::nrhs, mod_stepping::nstp, mod_ocean::ocean, tl_t3dmix4_geo_tile(), wclock_off(), and wclock_on().

Referenced by tl_rhs3d_mod::tl_rhs3d().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ tl_t3dmix4_geo_tile()

subroutine tl_t3dmix4_mod::tl_t3dmix4_geo_tile ( integer, intent(in) ng,
integer, intent(in) tile,
integer, intent(in) lbi,
integer, intent(in) ubi,
integer, intent(in) lbj,
integer, intent(in) ubj,
integer, intent(in) imins,
integer, intent(in) imaxs,
integer, intent(in) jmins,
integer, intent(in) jmaxs,
integer, intent(in) nrhs,
integer, intent(in) nstp,
integer, intent(in) nnew,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) umask,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) vmask,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) umask_wet,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) vmask_wet,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) om_v,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) on_u,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pm,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pn,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) hz,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) tl_hz,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) z_r,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) tl_z_r,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) diff3d_u,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) diff3d_v,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) diff3d_r,
real(r8), dimension(lbi:ubi,lbj:ubj,nt(ng)), intent(in) diff4,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),nt(ng)), intent(in) tclm,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),3,nt(ng)), intent(in) t,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),3,nt(ng)), intent(inout) tl_t )
private

Definition at line 100 of file tl_t3dmix4_geo.h.

129!***********************************************************************
130!
131 USE mod_param
132 USE mod_ncparam
133 USE mod_scalars
134!
135! Imported variable declarations.
136!
137 integer, intent(in) :: ng, tile
138 integer, intent(in) :: LBi, UBi, LBj, UBj
139 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
140 integer, intent(in) :: nrhs, nstp, nnew
141
142#ifdef ASSUMED_SHAPE
143# ifdef MASKING
144 real(r8), intent(in) :: umask(LBi:,LBj:)
145 real(r8), intent(in) :: vmask(LBi:,LBj:)
146# endif
147# ifdef WET_DRY_NOT_YET
148 real(r8), intent(in) :: umask_wet(LBi:,LBj:)
149 real(r8), intent(in) :: vmask_wet(LBi:,LBj:)
150# endif
151# ifdef DIFF_3DCOEF
152# ifdef TS_U3ADV_SPLIT
153 real(r8), intent(in) :: diff3d_u(LBi:,LBj:,:)
154 real(r8), intent(in) :: diff3d_v(LBi:,LBj:,:)
155# else
156 real(r8), intent(in) :: diff3d_r(LBi:,LBj:,:)
157# endif
158# else
159 real(r8), intent(in) :: diff4(LBi:,LBj:,:)
160# endif
161 real(r8), intent(in) :: om_v(LBi:,LBj:)
162 real(r8), intent(in) :: on_u(LBi:,LBj:)
163 real(r8), intent(in) :: pm(LBi:,LBj:)
164 real(r8), intent(in) :: pn(LBi:,LBj:)
165 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
166 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
167 real(r8), intent(in) :: t(LBi:,LBj:,:,:,:)
168# ifdef TS_MIX_CLIMA
169 real(r8), intent(in) :: tclm(LBi:,LBj:,:,:)
170# endif
171 real(r8), intent(in) :: tl_Hz(LBi:,LBj:,:)
172 real(r8), intent(in) :: tl_z_r(LBi:,LBj:,:)
173# ifdef DIAGNOSTICS_TS
174 real(r8), intent(inout) :: DiaTwrk(LBi:,LBj:,:,:,:)
175# endif
176 real(r8), intent(inout) :: tl_t(LBi:,LBj:,:,:,:)
177#else
178# ifdef MASKING
179 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
180 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
181# endif
182# ifdef WET_DRY_NOT_YET
183 real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
184 real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
185# endif
186# ifdef DIFF_3DCOEF
187# ifdef TS_U3ADV_SPLIT
188 real(r8), intent(in) :: diff3d_u(LBi:UBi,LBj:UBj,N(ng))
189 real(r8), intent(in) :: diff3d_v(LBi:UBi,LBj:UBj,N(ng))
190# else
191 real(r8), intent(in) :: diff3d_r(LBi:UBi,LBj:UBj,N(ng))
192# endif
193# else
194 real(r8), intent(in) :: diff4(LBi:UBi,LBj:UBj,NT(ng))
195# endif
196 real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
197 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
198 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
199 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
200 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
201 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
202 real(r8), intent(in) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
203# ifdef TS_MIX_CLIMA
204 real(r8), intent(in) :: tclm(LBi:UBi,LBj:UBj,N(ng),NT(ng))
205# endif
206 real(r8), intent(in) :: tl_Hz(LBi:UBi,LBj:UBj,N(ng))
207 real(r8), intent(in) :: tl_z_r(LBi:UBi,LBj:UBj,N(ng))
208# ifdef DIAGNOSTICS_TS
209!! real(r8), intent(inout) :: DiaTwrk(LBi:UBi,LBj:UBj,N(ng),NT(ng), &
210!! & NDT)
211# endif
212 real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
213#endif
214!
215! Local variable declarations.
216!
217 integer :: Imin, Imax, Jmin, Jmax
218 integer :: i, itrc, j, k, k1, k2
219
220 real(r8) :: cff, cff1, cff2, cff3, cff4, dife, difx
221 real(r8) :: tl_cff, tl_cff1, tl_cff2, tl_cff3, tl_cff4
222
223 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: LapT
224
225 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: tl_LapT
226
227 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FE
228 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FX
229
230 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_FE
231 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_FX
232
233 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: FS
234 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dTde
235 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dTdx
236 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dTdz
237 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZde
238 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZdx
239
240 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_FS
241 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_dTde
242 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_dTdx
243 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_dTdz
244 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_dZde
245 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_dZdx
246
247#include "set_bounds.h"
248!
249!-----------------------------------------------------------------------
250! Compute tangent linear horizontal biharmonic diffusion along
251! geopotential surfaces. The biharmonic operator is computed by
252! applying the harmonic operator twice.
253!-----------------------------------------------------------------------
254!
255! Set local I- and J-ranges.
256!
257 IF (ewperiodic(ng)) THEN
258 imin=istr-1
259 imax=iend+1
260 ELSE
261 imin=max(istr-1,1)
262 imax=min(iend+1,lm(ng))
263 END IF
264 IF (nsperiodic(ng)) THEN
265 jmin=jstr-1
266 jmax=jend+1
267 ELSE
268 jmin=max(jstr-1,1)
269 jmax=min(jend+1,mm(ng))
270 END IF
271!
272! Compute horizontal and vertical gradients associated with the
273! first rotated harmonic operator. Notice the recursive blocking
274! sequence. The vertical placement of the gradients is:
275!
276! dTdx,dTde(:,:,k1) k rho-points
277! dTdx,dTde(:,:,k2) k+1 rho-points
278! FS,dTdz(:,:,k1) k-1/2 W-points
279! FS,dTdz(:,:,k2) k+1/2 W-points
280!
281#ifdef TS_MIX_STABILITY
282! In order to increase stability, the rotated biharmonic is applied
283! as: 3/4 t(:,:,:,nrhs,:) + 1/4 t(:,:,:,nstp,:).
284!
285#endif
286
287 t_loop : DO itrc=1,nt(ng)
288 k2=1
289 k_loop1 : DO k=0,n(ng)
290 k1=k2
291 k2=3-k1
292 IF (k.lt.n(ng)) THEN
293 DO j=jmin,jmax
294 DO i=imin,imax+1
295 cff=0.5_r8*(pm(i,j)+pm(i-1,j))
296#ifdef MASKING
297 cff=cff*umask(i,j)
298#endif
299#ifdef WET_DRY_NOT_YET
300 cff=cff*umask_wet(i,j)
301#endif
302 dzdx(i,j,k2)=cff*(z_r(i ,j,k+1)- &
303 & z_r(i-1,j,k+1))
304 tl_dzdx(i,j,k2)=cff*(tl_z_r(i ,j,k+1)- &
305 & tl_z_r(i-1,j,k+1))
306#if defined TS_MIX_STABILITY
307 dtdx(i,j,k2)=cff*(0.75_r8*(t(i ,j,k+1,nrhs,itrc)- &
308 & t(i-1,j,k+1,nrhs,itrc))+ &
309 & 0.25_r8*(t(i ,j,k+1,nstp,itrc)- &
310 & t(i-1,j,k+1,nstp,itrc)))
311 tl_dtdx(i,j,k2)=cff*
312 & (0.75_r8*(tl_t(i ,j,k+1,nrhs,itrc)- &
313 & tl_t(i-1,j,k+1,nrhs,itrc))+ &
314 & 0.25_r8*(tl_t(i ,j,k+1,nstp,itrc)- &
315 & tl_t(i-1,j,k+1,nstp,itrc)))
316#elif defined TS_MIX_CLIMA
317 IF (ltracerclm(itrc,ng)) THEN
318 dtdx(i,j,k2)=cff*((t(i ,j,k+1,nrhs,itrc)- &
319 & tclm(i ,j,k+1,itrc))- &
320 & (t(i-1,j,k+1,nrhs,itrc)- &
321 & tclm(i-1,j,k+1,itrc)))
322 ELSE
323 dtdx(i,j,k2)=cff*(t(i ,j,k+1,nrhs,itrc)- &
324 & t(i-1,j,k+1,nrhs,itrc))
325 END IF
326 tl_dtdx(i,j,k2)=cff*(tl_t(i ,j,k+1,nrhs,itrc)- &
327 & tl_t(i-1,j,k+1,nrhs,itrc))
328#else
329 dtdx(i,j,k2)=cff*(t(i ,j,k+1,nrhs,itrc)- &
330 & t(i-1,j,k+1,nrhs,itrc))
331 tl_dtdx(i,j,k2)=cff*(tl_t(i ,j,k+1,nrhs,itrc)- &
332 & tl_t(i-1,j,k+1,nrhs,itrc))
333#endif
334 END DO
335 END DO
336 DO j=jmin,jmax+1
337 DO i=imin,imax
338 cff=0.5_r8*(pn(i,j)+pn(i,j-1))
339#ifdef MASKING
340 cff=cff*vmask(i,j)
341#endif
342#ifdef WET_DRY_NOT_YET
343 cff=cff*vmask_wet(i,j)
344#endif
345 dzde(i,j,k2)=cff*(z_r(i,j ,k+1)- &
346 & z_r(i,j-1,k+1))
347 tl_dzde(i,j,k2)=cff*(tl_z_r(i,j ,k+1)- &
348 & tl_z_r(i,j-1,k+1))
349#if defined TS_MIX_STABILITY
350 dtde(i,j,k2)=cff*(0.75_r8*(t(i,j ,k+1,nrhs,itrc)- &
351 & t(i,j-1,k+1,nrhs,itrc))+ &
352 & 0.25_r8*(t(i,j ,k+1,nstp,itrc)- &
353 & t(i,j-1,k+1,nstp,itrc)))
354 tl_dtde(i,j,k2)=cff*
355 & (0.75_r8*(tl_t(i,j ,k+1,nrhs,itrc)- &
356 & tl_t(i,j-1,k+1,nrhs,itrc))+ &
357 & 0.25_r8*(tl_t(i,j ,k+1,nstp,itrc)- &
358 & tl_t(i,j-1,k+1,nstp,itrc)))
359#elif defined TS_MIX_CLIMA
360 IF (ltracerclm(itrc,ng)) THEN
361 dtde(i,j,k2)=cff*((t(i,j ,k+1,nrhs,itrc)- &
362 & tclm(i,j ,k+1,itrc))- &
363 & (t(i,j-1,k+1,nrhs,itrc)- &
364 & tclm(i,j-1,k+1,itrc)))
365 ELSE
366 dtde(i,j,k2)=cff*(t(i,j ,k+1,nrhs,itrc)- &
367 & t(i,j-1,k+1,nrhs,itrc))
368 END IF
369 tl_dtde(i,j,k2)=cff*(tl_t(i,j ,k+1,nrhs,itrc)- &
370 & tl_t(i,j-1,k+1,nrhs,itrc))
371#else
372 dtde(i,j,k2)=cff*(t(i,j ,k+1,nrhs,itrc)- &
373 & t(i,j-1,k+1,nrhs,itrc))
374 tl_dtde(i,j,k2)=cff*(tl_t(i,j ,k+1,nrhs,itrc)- &
375 & tl_t(i,j-1,k+1,nrhs,itrc))
376#endif
377 END DO
378 END DO
379 END IF
380 IF ((k.eq.0).or.(k.eq.n(ng))) THEN
381 DO j=jmin-1,jmax+1
382 DO i=imin-1,imax+1
383 dtdz(i,j,k2)=0.0_r8
384 tl_dtdz(i,j,k2)=0.0_r8
385 fs(i,j,k2)=0.0_r8
386 tl_fs(i,j,k2)=0.0_r8
387 END DO
388 END DO
389 ELSE
390 DO j=jmin-1,jmax+1
391 DO i=imin-1,imax+1
392 cff=1.0_r8/(z_r(i,j,k+1)- &
393 & z_r(i,j,k ))
394 tl_cff=-cff*cff*(tl_z_r(i,j,k+1)- &
395 & tl_z_r(i,j,k ))
396#if defined TS_MIX_STABILITY
397 dtdz(i,j,k2)=cff*(0.75_r8*(t(i,j,k+1,nrhs,itrc)- &
398 & t(i,j,k ,nrhs,itrc))+ &
399 & 0.25_r8*(t(i,j,k+1,nstp,itrc)- &
400 & t(i,j,k ,nstp,itrc)))
401 tl_dtdz(i,j,k2)=tl_cff*(0.75_r8*(t(i,j,k+1,nrhs,itrc)- &
402 & t(i,j,k ,nrhs,itrc))+ &
403 & 0.25_r8*(t(i,j,k+1,nstp,itrc)- &
404 & t(i,j,k ,nstp,itrc)))+&
405 & cff*(0.75_r8*(tl_t(i,j,k+1,nrhs,itrc)- &
406 & tl_t(i,j,k ,nrhs,itrc))+ &
407 & 0.25_r8*(tl_t(i,j,k+1,nstp,itrc)- &
408 & tl_t(i,j,k ,nstp,itrc)))
409#elif defined TS_MIX_CLIMA
410 IF (ltracerclm(itrc,ng)) THEN
411 dtdz(i,j,k2)=cff*((t(i,j,k+1,nrhs,itrc)- &
412 & tclm(i,j,k+1,itrc))- &
413 & (t(i,j,k ,nrhs,itrc)- &
414 & tclm(i,j,k ,itrc)))
415 tl_dtdz(i,j,k2)=tl_cff*((t(i,j,k+1,nrhs,itrc)- &
416 & tclm(i,j,k+1,itrc))- &
417 & (t(i,j,k ,nrhs,itrc)- &
418 & tclm(i,j,k ,itrc)))+ &
419 & cff*(tl_t(i,j,k+1,nrhs,itrc)- &
420 & tl_t(i,j,k ,nrhs,itrc))
421
422 ELSE
423 dtdz(i,j,k2)=cff*(t(i,j,k+1,nrhs,itrc)- &
424 & t(i,j,k ,nrhs,itrc))
425 tl_dtdz(i,j,k2)=tl_cff*(t(i,j,k+1,nrhs,itrc)- &
426 & t(i,j,k ,nrhs,itrc))+ &
427 & cff*(tl_t(i,j,k+1,nrhs,itrc)- &
428 & tl_t(i,j,k ,nrhs,itrc))
429 END IF
430#else
431 dtdz(i,j,k2)=cff*(t(i,j,k+1,nrhs,itrc)- &
432 & t(i,j,k ,nrhs,itrc))
433 tl_dtdz(i,j,k2)=tl_cff*(t(i,j,k+1,nrhs,itrc)- &
434 & t(i,j,k ,nrhs,itrc))+ &
435 & cff*(tl_t(i,j,k+1,nrhs,itrc)- &
436 & tl_t(i,j,k ,nrhs,itrc))
437#endif
438 END DO
439 END DO
440 END IF
441 IF (k.gt.0) THEN
442 DO j=jmin,jmax
443 DO i=imin,imax+1
444#ifdef DIFF_3DCOEF
445# ifdef TS_U3ADV_SPLIT
446 cff=0.5_r8*diff3d_u(i,j,k)*on_u(i,j)
447# else
448 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
449 & on_u(i,j)
450# endif
451#else
452 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
453 & on_u(i,j)
454#endif
455 fx(i,j)=cff* &
456 & (hz(i,j,k)+hz(i-1,j,k))* &
457 & (dtdx(i,j,k1)- &
458 & 0.5_r8*(min(dzdx(i,j,k1),0.0_r8)* &
459 & (dtdz(i-1,j,k1)+ &
460 & dtdz(i ,j,k2))+ &
461 & max(dzdx(i,j,k1),0.0_r8)* &
462 & (dtdz(i-1,j,k2)+ &
463 & dtdz(i ,j,k1))))
464 tl_fx(i,j)=cff* &
465 & ((tl_hz(i,j,k)+tl_hz(i-1,j,k))* &
466 & (dtdx(i,j,k1)- &
467 & 0.5_r8*(min(dzdx(i,j,k1),0.0_r8)* &
468 & (dtdz(i-1,j,k1)+ &
469 & dtdz(i ,j,k2))+ &
470 & max(dzdx(i,j,k1),0.0_r8)* &
471 & (dtdz(i-1,j,k2)+ &
472 & dtdz(i ,j,k1))))+ &
473 & (hz(i,j,k)+hz(i-1,j,k))* &
474 & (tl_dtdx(i,j,k1)- &
475 & 0.5_r8*(min(dzdx(i,j,k1),0.0_r8)* &
476 & (tl_dtdz(i-1,j,k1)+ &
477 & tl_dtdz(i ,j,k2))+ &
478 & max(dzdx(i,j,k1),0.0_r8)* &
479 & (tl_dtdz(i-1,j,k2)+ &
480 & tl_dtdz(i ,j,k1)))- &
481 & 0.5_r8*((0.5_r8+ &
482 & sign(0.5_r8,-dzdx(i,j,k1)))* &
483 & tl_dzdx(i,j,k1)* &
484 & (dtdz(i-1,j,k1)+dtdz(i,j,k2))+ &
485 & (0.5_r8+ &
486 & sign(0.5_r8, dzdx(i,j,k1)))* &
487 & tl_dzdx(i,j,k1)* &
488 & (dtdz(i-1,j,k2)+dtdz(i,j,k1)))))
489 END DO
490 END DO
491 DO j=jmin,jmax+1
492 DO i=imin,imax
493#ifdef DIFF_3DCOEF
494# ifdef TS_U3ADV_SPLIT
495 cff=0.5_r8*diff3d_v(i,j,k)*om_v(i,j)
496# else
497 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
498 & om_v(i,j)
499# endif
500#else
501 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
502 & om_v(i,j)
503#endif
504 fe(i,j)=cff* &
505 & (hz(i,j,k)+hz(i,j-1,k))* &
506 & (dtde(i,j,k1)- &
507 & 0.5_r8*(min(dzde(i,j,k1),0.0_r8)* &
508 & (dtdz(i,j-1,k1)+ &
509 & dtdz(i,j ,k2))+ &
510 & max(dzde(i,j,k1),0.0_r8)* &
511 & (dtdz(i,j-1,k2)+ &
512 & dtdz(i,j ,k1))))
513 tl_fe(i,j)=cff* &
514 & ((tl_hz(i,j,k)+tl_hz(i,j-1,k))* &
515 & (dtde(i,j,k1)- &
516 & 0.5_r8*(min(dzde(i,j,k1),0.0_r8)* &
517 & (dtdz(i,j-1,k1)+ &
518 & dtdz(i,j ,k2))+ &
519 & max(dzde(i,j,k1),0.0_r8)* &
520 & (dtdz(i,j-1,k2)+ &
521 & dtdz(i,j ,k1))))+ &
522 & (hz(i,j,k)+hz(i,j-1,k))* &
523 & (tl_dtde(i,j,k1)- &
524 & 0.5_r8*(min(dzde(i,j,k1),0.0_r8)* &
525 & (tl_dtdz(i,j-1,k1)+ &
526 & tl_dtdz(i,j ,k2))+ &
527 & max(dzde(i,j,k1),0.0_r8)* &
528 & (tl_dtdz(i,j-1,k2)+ &
529 & tl_dtdz(i,j ,k1)))- &
530 & 0.5_r8*((0.5_r8+ &
531 & sign(0.5_r8,-dzde(i,j,k1)))* &
532 & tl_dzde(i,j,k1)* &
533 & (dtdz(i,j-1,k1)+dtdz(i,j,k2))+ &
534 & (0.5_r8+ &
535 & sign(0.5_r8, dzde(i,j,k1)))* &
536 & tl_dzde(i,j,k1)* &
537 & (dtdz(i,j-1,k2)+dtdz(i,j,k1)))))
538 END DO
539 END DO
540 IF (k.lt.n(ng)) THEN
541 DO j=jmin,jmax
542 DO i=imin,imax
543#ifdef DIFF_3DCOEF
544# ifdef TS_U3ADV_SPLIT
545 difx=0.125_r8*(diff3d_u(i,j,k )+diff3d_u(i+1,j,k )+ &
546 & diff3d_u(i,j,k+1)+diff3d_u(i+1,j,k+1))
547 dife=0.125_r8*(diff3d_v(i,j,k )+diff3d_v(i,j+1,k )+ &
548 & diff3d_v(i,j,k+1)+diff3d_v(i,j+1,k+1))
549# else
550 difx=0.5_r8*diff3d_r(i,j,k)
551 dife=difx
552# endif
553#else
554 difx=0.5_r8*diff4(i,j,itrc)
555 dife=difx
556#endif
557 cff1=min(dzdx(i ,j,k1),0.0_r8)
558 cff2=min(dzdx(i+1,j,k2),0.0_r8)
559 cff3=max(dzdx(i ,j,k2),0.0_r8)
560 cff4=max(dzdx(i+1,j,k1),0.0_r8)
561 tl_cff1=(0.5_r8+sign(0.5_r8,-dzdx(i ,j,k1)))* &
562 & tl_dzdx(i ,j,k1)
563 tl_cff2=(0.5_r8+sign(0.5_r8,-dzdx(i+1,j,k2)))* &
564 & tl_dzdx(i+1,j,k2)
565 tl_cff3=(0.5_r8+sign(0.5_r8, dzdx(i ,j,k2)))* &
566 & tl_dzdx(i ,j,k2)
567 tl_cff4=(0.5_r8+sign(0.5_r8, dzdx(i+1,j,k1)))* &
568 & tl_dzdx(i+1,j,k1)
569 fs(i,j,k2)=difx* &
570 & (cff1*(cff1*dtdz(i,j,k2)- &
571 & dtdx(i ,j,k1))+ &
572 & cff2*(cff2*dtdz(i,j,k2)- &
573 & dtdx(i+1,j,k2))+ &
574 & cff3*(cff3*dtdz(i,j,k2)- &
575 & dtdx(i ,j,k2))+ &
576 & cff4*(cff4*dtdz(i,j,k2)- &
577 & dtdx(i+1,j,k1)))
578 tl_fs(i,j,k2)=difx* &
579 & (tl_cff1*(cff1*dtdz(i,j,k2)- &
580 & dtdx(i ,j,k1))+ &
581 & tl_cff2*(cff2*dtdz(i,j,k2)- &
582 & dtdx(i+1,j,k2))+ &
583 & tl_cff3*(cff3*dtdz(i,j,k2)- &
584 & dtdx(i ,j,k2))+ &
585 & tl_cff4*(cff4*dtdz(i,j,k2)- &
586 & dtdx(i+1,j,k1))+ &
587 & cff1*(tl_cff1*dtdz(i,j,k2)+ &
588 & cff1*tl_dtdz(i,j,k2)- &
589 & tl_dtdx(i ,j,k1))+ &
590 & cff2*(tl_cff2*dtdz(i,j,k2)+ &
591 & cff2*tl_dtdz(i,j,k2)- &
592 & tl_dtdx(i+1,j,k2))+ &
593 & cff3*(tl_cff3*dtdz(i,j,k2)+ &
594 & cff3*tl_dtdz(i,j,k2)- &
595 & tl_dtdx(i ,j,k2))+ &
596 & cff4*(tl_cff4*dtdz(i,j,k2)+ &
597 & cff4*tl_dtdz(i,j,k2)- &
598 & tl_dtdx(i+1,j,k1)))
599!
600 cff1=min(dzde(i,j ,k1),0.0_r8)
601 cff2=min(dzde(i,j+1,k2),0.0_r8)
602 cff3=max(dzde(i,j ,k2),0.0_r8)
603 cff4=max(dzde(i,j+1,k1),0.0_r8)
604 tl_cff1=(0.5_r8+sign(0.5_r8,-dzde(i,j ,k1)))* &
605 & tl_dzde(i,j ,k1)
606 tl_cff2=(0.5_r8+sign(0.5_r8,-dzde(i,j+1,k2)))* &
607 & tl_dzde(i,j+1,k2)
608 tl_cff3=(0.5_r8+sign(0.5_r8, dzde(i,j ,k2)))* &
609 & tl_dzde(i,j ,k2)
610 tl_cff4=(0.5_r8+sign(0.5_r8, dzde(i,j+1,k1)))* &
611 & tl_dzde(i,j+1,k1)
612 fs(i,j,k2)=fs(i,j,k2)+ &
613 & dife* &
614 & (cff1*(cff1*dtdz(i,j,k2)- &
615 & dtde(i,j ,k1))+ &
616 & cff2*(cff2*dtdz(i,j,k2)- &
617 & dtde(i,j+1,k2))+ &
618 & cff3*(cff3*dtdz(i,j,k2)- &
619 & dtde(i,j ,k2))+ &
620 & cff4*(cff4*dtdz(i,j,k2)- &
621 & dtde(i,j+1,k1)))
622 tl_fs(i,j,k2)=tl_fs(i,j,k2)+ &
623 & dife* &
624 & (tl_cff1*(cff1*dtdz(i,j,k2)- &
625 & dtde(i,j ,k1))+ &
626 & tl_cff2*(cff2*dtdz(i,j,k2)- &
627 & dtde(i,j+1,k2))+ &
628 & tl_cff3*(cff3*dtdz(i,j,k2)- &
629 & dtde(i,j ,k2))+ &
630 & tl_cff4*(cff4*dtdz(i,j,k2)- &
631 & dtde(i,j+1,k1))+ &
632 & cff1*(tl_cff1*dtdz(i,j,k2)+ &
633 & cff1*tl_dtdz(i,j,k2)- &
634 & tl_dtde(i,j ,k1))+ &
635 & cff2*(tl_cff2*dtdz(i,j,k2)+ &
636 & cff2*tl_dtdz(i,j,k2)- &
637 & tl_dtde(i,j+1,k2))+ &
638 & cff3*(tl_cff3*dtdz(i,j,k2)+ &
639 & cff3*tl_dtdz(i,j,k2)- &
640 & tl_dtde(i,j ,k2))+ &
641 & cff4*(tl_cff4*dtdz(i,j,k2)+ &
642 & cff4*tl_dtdz(i,j,k2)- &
643 & tl_dtde(i,j+1,k1)))
644 END DO
645 END DO
646 END IF
647!
648! Compute first harmonic operator, without mixing coefficient.
649! Multiply by the metrics of the second harmonic operator. Save
650! into work array "LapT".
651!
652 DO j=jmin,jmax
653 DO i=imin,imax
654 cff=pm(i,j)*pn(i,j)
655 cff1=1.0_r8/hz(i,j,k)
656 tl_cff1=-cff1*cff1*tl_hz(i,j,k)
657 lapt(i,j,k)=cff1*(cff* &
658 & (fx(i+1,j)-fx(i,j)+ &
659 & fe(i,j+1)-fe(i,j))+ &
660 & (fs(i,j,k2)-fs(i,j,k1)))
661 tl_lapt(i,j,k)=tl_cff1*(cff* &
662 & (fx(i+1,j)-fx(i,j)+ &
663 & fe(i,j+1)-fe(i,j))+ &
664 & (fs(i,j,k2)-fs(i,j,k1)))+ &
665 & cff1*(cff* &
666 & (tl_fx(i+1,j)-tl_fx(i,j)+ &
667 & tl_fe(i,j+1)-tl_fe(i,j))+ &
668 & (tl_fs(i,j,k2)-tl_fs(i,j,k1)))
669 END DO
670 END DO
671 END IF
672 END DO k_loop1
673!
674! Apply boundary conditions (except periodic; closed or gradient)
675! to the first harmonic operator.
676!
677 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
678 IF (domain(ng)%Western_Edge(tile)) THEN
679 IF (tl_lbc(iwest,istvar(itrc),ng)%closed) THEN
680 DO k=1,n(ng)
681 DO j=jmin,jmax
682 lapt(istr-1,j,k)=0.0_r8
683 tl_lapt(istr-1,j,k)=0.0_r8
684 END DO
685 END DO
686 ELSE
687 DO k=1,n(ng)
688 DO j=jmin,jmax
689 lapt(istr-1,j,k)=lapt(istr,j,k)
690 tl_lapt(istr-1,j,k)=tl_lapt(istr,j,k)
691 END DO
692 END DO
693 END IF
694 END IF
695 END IF
696!
697 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
698 IF (domain(ng)%Eastern_Edge(tile)) THEN
699 IF (tl_lbc(ieast,istvar(itrc),ng)%closed) THEN
700 DO k=1,n(ng)
701 DO j=jmin,jmax
702 lapt(iend+1,j,k)=0.0_r8
703 tl_lapt(iend+1,j,k)=0.0_r8
704 END DO
705 END DO
706 ELSE
707 DO k=1,n(ng)
708 DO j=jmin,jmax
709 lapt(iend+1,j,k)=lapt(iend,j,k)
710 tl_lapt(iend+1,j,k)=tl_lapt(iend,j,k)
711 END DO
712 END DO
713 END IF
714 END IF
715 END IF
716!
717 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
718 IF (domain(ng)%Southern_Edge(tile)) THEN
719 IF (tl_lbc(isouth,istvar(itrc),ng)%closed) THEN
720 DO k=1,n(ng)
721 DO i=imin,imax
722 lapt(i,jstr-1,k)=0.0_r8
723 tl_lapt(i,jstr-1,k)=0.0_r8
724 END DO
725 END DO
726 ELSE
727 DO k=1,n(ng)
728 DO i=imin,imax
729 lapt(i,jstr-1,k)=lapt(i,jstr,k)
730 tl_lapt(i,jstr-1,k)=tl_lapt(i,jstr,k)
731 END DO
732 END DO
733 END IF
734 END IF
735 END IF
736!
737 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
738 IF (domain(ng)%Northern_Edge(tile)) THEN
739 IF (tl_lbc(inorth,istvar(itrc),ng)%closed) THEN
740 DO k=1,n(ng)
741 DO i=imin,imax
742 lapt(i,jend+1,k)=0.0_r8
743 tl_lapt(i,jend+1,k)=0.0_r8
744 END DO
745 END DO
746 ELSE
747 DO k=1,n(ng)
748 DO i=imin,imax
749 lapt(i,jend+1,k)=lapt(i,jend,k)
750 tl_lapt(i,jend+1,k)=tl_lapt(i,jend,k)
751 END DO
752 END DO
753 END IF
754 END IF
755 END IF
756!
757 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
758 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
759 IF (domain(ng)%SouthWest_Corner(tile)) THEN
760 DO k=1,n(ng)
761 lapt(istr-1,jstr-1,k)=0.5_r8* &
762 & (lapt(istr ,jstr-1,k)+ &
763 & lapt(istr-1,jstr ,k))
764 tl_lapt(istr-1,jstr-1,k)=0.5_r8* &
765 & (tl_lapt(istr ,jstr-1,k)+ &
766 & tl_lapt(istr-1,jstr ,k))
767 END DO
768 END IF
769 END IF
770
771 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
772 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
773 IF (domain(ng)%SouthEast_Corner(tile)) THEN
774 DO k=1,n(ng)
775 lapt(iend+1,jstr-1,k)=0.5_r8* &
776 & (lapt(iend ,jstr-1,k)+ &
777 & lapt(iend+1,jstr ,k))
778 tl_lapt(iend+1,jstr-1,k)=0.5_r8* &
779 & (tl_lapt(iend ,jstr-1,k)+ &
780 & tl_lapt(iend+1,jstr ,k))
781 END DO
782 END IF
783 END IF
784
785 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
786 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
787 IF (domain(ng)%NorthWest_Corner(tile)) THEN
788 DO k=1,n(ng)
789 lapt(istr-1,jend+1,k)=0.5_r8* &
790 & (lapt(istr ,jend+1,k)+ &
791 & lapt(istr-1,jend ,k))
792 tl_lapt(istr-1,jend+1,k)=0.5_r8* &
793 & (tl_lapt(istr ,jend+1,k)+ &
794 & tl_lapt(istr-1,jend ,k))
795 END DO
796 END IF
797 END IF
798
799 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
800 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
801 IF (domain(ng)%NorthEast_Corner(tile)) THEN
802 DO k=1,n(ng)
803 lapt(iend+1,jend+1,k)=0.5_r8* &
804 & (lapt(iend ,jend+1,k)+ &
805 & lapt(iend+1,jend ,k))
806 tl_lapt(iend+1,jend+1,k)=0.5_r8* &
807 & (tl_lapt(iend ,jend+1,k)+ &
808 & tl_lapt(iend+1,jend ,k))
809 END DO
810 END IF
811 END IF
812!
813! Compute tangent linear horizontal and vertical gradients associated
814! with the second rotated harmonic operator.
815!
816 k2=1
817 k_loop2: DO k=0,n(ng)
818 k1=k2
819 k2=3-k1
820 IF (k.lt.n(ng)) THEN
821 DO j=jstr,jend
822 DO i=istr,iend+1
823 cff=0.5_r8*(pm(i,j)+pm(i-1,j))
824#ifdef MASKING
825 cff=cff*umask(i,j)
826#endif
827#ifdef WET_DRY_NOT_YET
828 cff=cff*umask_wet(i,j)
829#endif
830 dzdx(i,j,k2)=cff*(z_r(i ,j,k+1)- &
831 & z_r(i-1,j,k+1))
832 tl_dzdx(i,j,k2)=cff*(tl_z_r(i ,j,k+1)- &
833 & tl_z_r(i-1,j,k+1))
834 dtdx(i,j,k2)=cff*(lapt(i ,j,k+1)- &
835 & lapt(i-1,j,k+1))
836 tl_dtdx(i,j,k2)=cff*(tl_lapt(i ,j,k+1)- &
837 & tl_lapt(i-1,j,k+1))
838 END DO
839 END DO
840 DO j=jstr,jend+1
841 DO i=istr,iend
842 cff=0.5_r8*(pn(i,j)+pn(i,j-1))
843#ifdef MASKING
844 cff=cff*vmask(i,j)
845#endif
846#ifdef WET_DRY_NOT_YET
847 cff=cff*vmask_wet(i,j)
848#endif
849 dzde(i,j,k2)=cff*(z_r(i,j ,k+1)- &
850 & z_r(i,j-1,k+1))
851 tl_dzde(i,j,k2)=cff*(tl_z_r(i,j ,k+1)- &
852 & tl_z_r(i,j-1,k+1))
853 dtde(i,j,k2)=cff*(lapt(i,j ,k+1)- &
854 & lapt(i,j-1,k+1))
855 tl_dtde(i,j,k2)=cff*(tl_lapt(i,j ,k+1)- &
856 & tl_lapt(i,j-1,k+1))
857 END DO
858 END DO
859 END IF
860 IF ((k.eq.0).or.(k.eq.n(ng))) THEN
861 DO j=jstr-1,jend+1
862 DO i=istr-1,iend+1
863 dtdz(i,j,k2)=0.0_r8
864 tl_dtdz(i,j,k2)=0.0_r8
865 fs(i,j,k2)=0.0_r8
866 tl_fs(i,j,k2)=0.0_r8
867 END DO
868 END DO
869 ELSE
870 DO j=jstr-1,jend+1
871 DO i=istr-1,iend+1
872 cff=1.0_r8/(z_r(i,j,k+1)- &
873 & z_r(i,j,k ))
874 tl_cff=-cff*cff*(tl_z_r(i,j,k+1)- &
875 & tl_z_r(i,j,k ))
876 dtdz(i,j,k2)=cff*(lapt(i,j,k+1)- &
877 & lapt(i,j,k ))
878 tl_dtdz(i,j,k2)=tl_cff*(lapt(i,j,k+1)- &
879 & lapt(i,j,k ))+ &
880 & cff*(tl_lapt(i,j,k+1)- &
881 & tl_lapt(i,j,k ))
882 END DO
883 END DO
884 END IF
885!
886! Compute components of the rotated tracer flux (T m4/s) along
887! geopotential surfaces.
888!
889 IF (k.gt.0) THEN
890 DO j=jstr,jend
891 DO i=istr,iend+1
892#ifdef DIFF_3DCOEF
893# ifdef TS_U3ADV_SPLIT
894 cff=0.5_r8*diff3d_u(i,j,k)*on_u(i,j)
895# else
896 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
897 & on_u(i,j)
898# endif
899#else
900 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
901 & on_u(i,j)
902#endif
903!^ FX(i,j)=cff*
904!^ & (Hz(i,j,k)+Hz(i-1,j,k))* &
905!^ & (dTdx(i ,j,k1)- &
906!^ & 0.5_r8*(MIN(dZdx(i,j,k1),0.0_r8)* &
907!^ & (dTdz(i-1,j,k1)+ &
908!^ & dTdz(i ,j,k2))+ &
909!^ & MAX(dZdx(i,j,k1),0.0_r8)* &
910!^ & (dTdz(i-1,j,k2)+ &
911!^ & dTdz(i ,j,k1))))
912!^
913 tl_fx(i,j)=cff* &
914 & ((tl_hz(i,j,k)+tl_hz(i-1,j,k))* &
915 & (dtdx(i,j,k1)- &
916 & 0.5_r8*(min(dzdx(i,j,k1),0.0_r8)* &
917 & (dtdz(i-1,j,k1)+ &
918 & dtdz(i ,j,k2))+ &
919 & max(dzdx(i,j,k1),0.0_r8)* &
920 & (dtdz(i-1,j,k2)+ &
921 & dtdz(i ,j,k1))))+ &
922 & (hz(i,j,k)+hz(i-1,j,k))* &
923 & (tl_dtdx(i,j,k1)- &
924 & 0.5_r8*(min(dzdx(i,j,k1),0.0_r8)* &
925 & (tl_dtdz(i-1,j,k1)+ &
926 & tl_dtdz(i ,j,k2))+ &
927 & max(dzdx(i,j,k1),0.0_r8)* &
928 & (tl_dtdz(i-1,j,k2)+ &
929 & tl_dtdz(i ,j,k1)))- &
930 & 0.5_r8*((0.5_r8+ &
931 & sign(0.5_r8,-dzdx(i,j,k1)))* &
932 & tl_dzdx(i,j,k1)* &
933 & (dtdz(i-1,j,k1)+dtdz(i,j,k2))+ &
934 & (0.5_r8+ &
935 & sign(0.5_r8, dzdx(i,j,k1)))* &
936 & tl_dzdx(i,j,k1)* &
937 & (dtdz(i-1,j,k2)+dtdz(i,j,k1)))))
938 END DO
939 END DO
940 DO j=jstr,jend+1
941 DO i=istr,iend
942#ifdef DIFF_3DCOEF
943# ifdef TS_U3ADV_SPLIT
944 cff=0.5_r8*diff3d_v(i,j,k)*om_v(i,j)
945# else
946 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
947 & om_v(i,j)
948# endif
949#else
950 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
951 & om_v(i,j)
952#endif
953!^ FE(i,j)=cff* &
954!^ & (Hz(i,j,k)+Hz(i,j-1,k))* &
955!^ & (dTde(i,j,k1)- &
956!^ & 0.5_r8*(MIN(dZde(i,j,k1),0.0_r8)* &
957!^ & (dTdz(i,j-1,k1)+ &
958!^ & dTdz(i,j ,k2))+ &
959!^ & MAX(dZde(i,j,k1),0.0_r8)* &
960!^ & (dTdz(i,j-1,k2)+ &
961!^ & dTdz(i,j ,k1))))
962!^
963 tl_fe(i,j)=cff* &
964 & ((tl_hz(i,j,k)+tl_hz(i,j-1,k))* &
965 & (dtde(i,j,k1)- &
966 & 0.5_r8*(min(dzde(i,j,k1),0.0_r8)* &
967 & (dtdz(i,j-1,k1)+ &
968 & dtdz(i,j ,k2))+ &
969 & max(dzde(i,j,k1),0.0_r8)* &
970 & (dtdz(i,j-1,k2)+ &
971 & dtdz(i,j ,k1))))+ &
972 & (hz(i,j,k)+hz(i,j-1,k))* &
973 & (tl_dtde(i,j,k1)- &
974 & 0.5_r8*(min(dzde(i,j,k1),0.0_r8)* &
975 & (tl_dtdz(i,j-1,k1)+ &
976 & tl_dtdz(i,j ,k2))+ &
977 & max(dzde(i,j,k1),0.0_r8)* &
978 & (tl_dtdz(i,j-1,k2)+ &
979 & tl_dtdz(i,j ,k1)))- &
980 & 0.5_r8*((0.5_r8+ &
981 & sign(0.5_r8,-dzde(i,j,k1)))* &
982 & tl_dzde(i,j,k1)* &
983 & (dtdz(i,j-1,k1)+dtdz(i,j,k2))+ &
984 & (0.5_r8+ &
985 & sign(0.5_r8, dzde(i,j,k1)))* &
986 & tl_dzde(i,j,k1)* &
987 & (dtdz(i,j-1,k2)+dtdz(i,j,k1)))))
988 END DO
989 END DO
990 IF (k.lt.n(ng)) THEN
991 DO j=jstr,jend
992 DO i=istr,iend
993#ifdef DIFF_3DCOEF
994# ifdef TS_U3ADV_SPLIT
995 difx=0.125_r8*(diff3d_u(i,j,k )+diff3d_u(i+1,j,k )+ &
996 & diff3d_u(i,j,k+1)+diff3d_u(i+1,j,k+1))
997 dife=0.125_r8*(diff3d_v(i,j,k )+diff3d_v(i,j+1,k )+ &
998 & diff3d_v(i,j,k+1)+diff3d_v(i,j+1,k+1))
999# else
1000 difx=0.5_r8*diff3d_r(i,j,k)
1001 dife=difx
1002# endif
1003#else
1004 difx=0.5_r8*diff4(i,j,itrc)
1005 dife=difx
1006#endif
1007 cff1=min(dzdx(i ,j,k1),0.0_r8)
1008 cff2=min(dzdx(i+1,j,k2),0.0_r8)
1009 cff3=max(dzdx(i ,j,k2),0.0_r8)
1010 cff4=max(dzdx(i+1,j,k1),0.0_r8)
1011 tl_cff1=(0.5_r8+sign(0.5_r8,-dzdx(i ,j,k1)))* &
1012 & tl_dzdx(i ,j,k1)
1013 tl_cff2=(0.5_r8+sign(0.5_r8,-dzdx(i+1,j,k2)))* &
1014 & tl_dzdx(i+1,j,k2)
1015 tl_cff3=(0.5_r8+sign(0.5_r8, dzdx(i ,j,k2)))* &
1016 & tl_dzdx(i ,j,k2)
1017 tl_cff4=(0.5_r8+sign(0.5_r8, dzdx(i+1,j,k1)))* &
1018 & tl_dzdx(i+1,j,k1)
1019!^
1020!^ FS(i,j,k2)=difx* &
1021!^ & (cff1*(cff1*dTdz(i,j,k2)- &
1022!^ & dTdx(i ,j,k1))+ &
1023!^ & cff2*(cff2*dTdz(i,j,k2)- &
1024!^ & dTdx(i+1,j,k2))+ &
1025!^ & cff3*(cff3*dTdz(i,j,k2)- &
1026!^ & dTdx(i ,j,k2))+ &
1027!^ & cff4*(cff4*dTdz(i,j,k2)- &
1028!^ & dTdx(i+1,j,k1)))
1029!^
1030 tl_fs(i,j,k2)=difx* &
1031 & (tl_cff1*(cff1*dtdz(i,j,k2)- &
1032 & dtdx(i ,j,k1))+ &
1033 & tl_cff2*(cff2*dtdz(i,j,k2)- &
1034 & dtdx(i+1,j,k2))+ &
1035 & tl_cff3*(cff3*dtdz(i,j,k2)- &
1036 & dtdx(i ,j,k2))+ &
1037 & tl_cff4*(cff4*dtdz(i,j,k2)- &
1038 & dtdx(i+1,j,k1))+ &
1039 & cff1*(tl_cff1*dtdz(i,j,k2)+ &
1040 & cff1*tl_dtdz(i,j,k2)- &
1041 & tl_dtdx(i ,j,k1))+ &
1042 & cff2*(tl_cff2*dtdz(i,j,k2)+ &
1043 & cff2*tl_dtdz(i,j,k2)- &
1044 & tl_dtdx(i+1,j,k2))+ &
1045 & cff3*(tl_cff3*dtdz(i,j,k2)+ &
1046 & cff3*tl_dtdz(i,j,k2)- &
1047 & tl_dtdx(i ,j,k2))+ &
1048 & cff4*(tl_cff4*dtdz(i,j,k2)+ &
1049 & cff4*tl_dtdz(i,j,k2)- &
1050 & tl_dtdx(i+1,j,k1)))
1051!
1052 cff1=min(dzde(i,j ,k1),0.0_r8)
1053 cff2=min(dzde(i,j+1,k2),0.0_r8)
1054 cff3=max(dzde(i,j ,k2),0.0_r8)
1055 cff4=max(dzde(i,j+1,k1),0.0_r8)
1056 tl_cff1=(0.5_r8+sign(0.5_r8,-dzde(i,j ,k1)))* &
1057 & tl_dzde(i,j ,k1)
1058 tl_cff2=(0.5_r8+sign(0.5_r8,-dzde(i,j+1,k2)))* &
1059 & tl_dzde(i,j+1,k2)
1060 tl_cff3=(0.5_r8+sign(0.5_r8, dzde(i,j ,k2)))* &
1061 & tl_dzde(i,j ,k2)
1062 tl_cff4=(0.5_r8+sign(0.5_r8, dzde(i,j+1,k1)))* &
1063 & tl_dzde(i,j+1,k1)
1064!^ FS(i,j,k2)=FS(i,j,k2)+ &
1065!^ & dife* &
1066!^ & (cff1*(cff1*dTdz(i,j,k2)- &
1067!^ & dTde(i,j ,k1))+ &
1068!^ & cff2*(cff2*dTdz(i,j,k2)- &
1069!^ & dTde(i,j+1,k2))+ &
1070!^ & cff3*(cff3*dTdz(i,j,k2)- &
1071!^ & dTde(i,j ,k2))+ &
1072!^ & cff4*(cff4*dTdz(i,j,k2)- &
1073!^ & dTde(i,j+1,k1)))
1074!^
1075 tl_fs(i,j,k2)=tl_fs(i,j,k2)+ &
1076 & dife* &
1077 & (tl_cff1*(cff1*dtdz(i,j,k2)- &
1078 & dtde(i,j ,k1))+ &
1079 & tl_cff2*(cff2*dtdz(i,j,k2)- &
1080 & dtde(i,j+1,k2))+ &
1081 & tl_cff3*(cff3*dtdz(i,j,k2)- &
1082 & dtde(i,j ,k2))+ &
1083 & tl_cff4*(cff4*dtdz(i,j,k2)- &
1084 & dtde(i,j+1,k1))+ &
1085 & cff1*(tl_cff1*dtdz(i,j,k2)+ &
1086 & cff1*tl_dtdz(i,j,k2)- &
1087 & tl_dtde(i,j ,k1))+ &
1088 & cff2*(tl_cff2*dtdz(i,j,k2)+ &
1089 & cff2*tl_dtdz(i,j,k2)- &
1090 & tl_dtde(i,j+1,k2))+ &
1091 & cff3*(tl_cff3*dtdz(i,j,k2)+ &
1092 & cff3*tl_dtdz(i,j,k2)- &
1093 & tl_dtde(i,j ,k2))+ &
1094 & cff4*(tl_cff4*dtdz(i,j,k2)+ &
1095 & cff4*tl_dtdz(i,j,k2)- &
1096 & tl_dtde(i,j+1,k1)))
1097 END DO
1098 END DO
1099 END IF
1100!
1101! Time-step biharmonic, geopotential diffusion term (m Tunits).
1102!
1103 DO j=jstr,jend
1104 DO i=istr,iend
1105!^ cff=dt(ng)*pm(i,j)*pn(i,j)* &
1106!^ & (FX(i+1,j )-FX(i,j)+ &
1107!^ & FE(i ,j+1)-FE(i,j))+ &
1108!^ & dt(ng)*(FS(i,j,k2)-FS(i,j,k1))
1109!^
1110 tl_cff=dt(ng)*pm(i,j)*pn(i,j)* &
1111 & (tl_fx(i+1,j)-tl_fx(i,j)+ &
1112 & tl_fe(i,j+1)-tl_fe(i,j))+ &
1113 & dt(ng)*(tl_fs(i,j,k2)-tl_fs(i,j,k1))
1114!^ t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)-cff
1115!^
1116 tl_t(i,j,k,nnew,itrc)=tl_t(i,j,k,nnew,itrc)-tl_cff
1117#ifdef DIAGNOSTICS_TS
1118!! DiaTwrk(i,j,k,itrc,iThdif)=-cff
1119#endif
1120 END DO
1121 END DO
1122 END IF
1123 END DO k_loop2
1124 END DO t_loop
1125!
1126 RETURN
mod_ncparam
Definition mod_ncparam.F:2
mod_ncparam::istvar
integer, dimension(:), allocatable istvar
Definition mod_ncparam.F:520
mod_param::n
integer, dimension(:), allocatable n
Definition mod_param.F:479
mod_param::tl_lbc
type(t_lbc), dimension(:,:,:), allocatable tl_lbc
Definition mod_param.F:379
mod_param::domain
type(t_domain), dimension(:), allocatable domain
Definition mod_param.F:329
mod_param::lm
integer, dimension(:), allocatable lm
Definition mod_param.F:455
mod_param::nt
integer, dimension(:), allocatable nt
Definition mod_param.F:489
mod_param::mm
integer, dimension(:), allocatable mm
Definition mod_param.F:456
mod_scalars
Definition mod_scalars.F:2
mod_scalars::dt
real(dp), dimension(:), allocatable dt
Definition mod_scalars.F:269
mod_scalars::ewperiodic
logical, dimension(:), allocatable ewperiodic
Definition mod_scalars.F:510
mod_scalars::iwest
integer, parameter iwest
Definition mod_scalars.F:1432
mod_scalars::nsperiodic
logical, dimension(:), allocatable nsperiodic
Definition mod_scalars.F:511
mod_scalars::compositegrid
logical, dimension(:,:), allocatable compositegrid
Definition mod_scalars.F:493
mod_scalars::isouth
integer, parameter isouth
Definition mod_scalars.F:1433
mod_scalars::ieast
integer, parameter ieast
Definition mod_scalars.F:1434
mod_scalars::ltracerclm
logical, dimension(:,:), allocatable ltracerclm
Definition mod_scalars.F:528
mod_scalars::inorth
integer, parameter inorth
Definition mod_scalars.F:1435

References mod_scalars::compositegrid, mod_param::domain, mod_scalars::dt, mod_scalars::ewperiodic, mod_scalars::ieast, mod_scalars::inorth, mod_scalars::isouth, mod_ncparam::istvar, mod_scalars::iwest, mod_param::lm, mod_scalars::ltracerclm, mod_param::mm, mod_scalars::nsperiodic, and mod_param::tl_lbc.

Referenced by tl_t3dmix4().

Here is the caller graph for this function:

◆ tl_t3dmix4_iso_tile()

subroutine tl_t3dmix4_mod::tl_t3dmix4_iso_tile ( integer, intent(in) ng,
integer, intent(in) tile,
integer, intent(in) lbi,
integer, intent(in) ubi,
integer, intent(in) lbj,
integer, intent(in) ubj,
integer, intent(in) imins,
integer, intent(in) imaxs,
integer, intent(in) jmins,
integer, intent(in) jmaxs,
integer, intent(in) nrhs,
integer, intent(in) nstp,
integer, intent(in) nnew,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) umask,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) vmask,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) umask_wet,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) vmask_wet,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) om_v,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) on_u,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pm,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pn,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) hz,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) tl_hz,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) z_r,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) tl_z_r,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) diff3d_u,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) diff3d_v,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) diff3d_r,
real(r8), dimension(lbi:ubi,lbj:ubj,nt(ng)), intent(in) diff4,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) pden,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) tl_pden,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),nt(ng)), intent(in) tclm,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),3,nt(ng)), intent(in) t,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),3,nt(ng)), intent(inout) tl_t )
private

Definition at line 102 of file tl_t3dmix4_iso.h.

132!***********************************************************************
133!
134 USE mod_param
135 USE mod_ncparam
136 USE mod_scalars
137!
138! Imported variable declarations.
139!
140 integer, intent(in) :: ng, tile
141 integer, intent(in) :: LBi, UBi, LBj, UBj
142 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
143 integer, intent(in) :: nrhs, nstp, nnew
144
145#ifdef ASSUMED_SHAPE
146# ifdef MASKING
147 real(r8), intent(in) :: umask(LBi:,LBj:)
148 real(r8), intent(in) :: vmask(LBi:,LBj:)
149# endif
150# ifdef WET_DRY_NOT_YET
151 real(r8), intent(in) :: umask_wet(LBi:,LBj:)
152 real(r8), intent(in) :: vmask_wet(LBi:,LBj:)
153# endif
154# ifdef DIFF_3DCOEF
155# ifdef TS_U3ADV_SPLIT
156 real(r8), intent(in) :: diff3d_u(LBi:,LBj:,:)
157 real(r8), intent(in) :: diff3d_v(LBi:,LBj:,:)
158# else
159 real(r8), intent(in) :: diff3d_r(LBi:,LBj:,:)
160# endif
161# else
162 real(r8), intent(in) :: diff4(LBi:,LBj:,:)
163# endif
164 real(r8), intent(in) :: om_v(LBi:,LBj:)
165 real(r8), intent(in) :: on_u(LBi:,LBj:)
166 real(r8), intent(in) :: pm(LBi:,LBj:)
167 real(r8), intent(in) :: pn(LBi:,LBj:)
168 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
169 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
170 real(r8), intent(in) :: pden(LBi:,LBj:,:)
171 real(r8), intent(in) :: t(LBi:,LBj:,:,:,:)
172# ifdef TS_MIX_CLIMA
173 real(r8), intent(in) :: tclm(LBi:,LBj:,:,:)
174# endif
175 real(r8), intent(in) :: tl_Hz(LBi:,LBj:,:)
176 real(r8), intent(in) :: tl_z_r(LBi:,LBj:,:)
177 real(r8), intent(in) :: tl_pden(LBi:,LBj:,:)
178# ifdef DIAGNOSTICS_TS
179 real(r8), intent(inout) :: DiaTwrk(LBi:,LBj:,:,:,:)
180# endif
181
182 real(r8), intent(inout) :: tl_t(LBi:,LBj:,:,:,:)
183#else
184# ifdef MASKING
185 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
186 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
187# endif
188# ifdef WET_DRY_NOT_YET
189 real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
190 real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
191# endif
192# ifdef DIFF_3DCOEF
193# ifdef TS_U3ADV_SPLIT
194 real(r8), intent(in) :: diff3d_u(LBi:UBi,LBj:UBj,N(ng))
195 real(r8), intent(in) :: diff3d_v(LBi:UBi,LBj:UBj,N(ng))
196# else
197 real(r8), intent(in) :: diff3d_r(LBi:UBi,LBj:UBj,N(ng))
198# endif
199# else
200 real(r8), intent(in) :: diff4(LBi:UBi,LBj:UBj,NT(ng))
201# endif
202 real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
203 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
204 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
205 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
206 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
207 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
208 real(r8), intent(in) :: pden(LBi:UBi,LBj:UBj,N(ng))
209 real(r8), intent(in) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
210# ifdef TS_MIX_CLIMA
211 real(r8), intent(in) :: tclm(LBi:UBi,LBj:UBj,N(ng),NT(ng))
212# endif
213 real(r8), intent(in) :: tl_Hz(LBi:UBi,LBj:UBj,N(ng))
214 real(r8), intent(in) :: tl_z_r(LBi:UBi,LBj:UBj,N(ng))
215 real(r8), intent(in) :: tl_pden(LBi:UBi,LBj:UBj,N(ng))
216# ifdef DIAGNOSTICS_TS
217!! real(r8), intent(inout) :: DiaTwrk(LBi:UBi,LBj:UBj,N(ng),NT(ng), &
218!! & NDT)
219# endif
220 real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
221#endif
222!
223! Local variable declarations.
224!
225 integer :: Imin, Imax, Jmin, Jmax
226 integer :: i, itrc, j, k, k1, k2
227
228 real(r8), parameter :: eps = 0.5_r8
229 real(r8), parameter :: small = 1.0e-14_r8
230 real(r8), parameter :: slope_max = 0.0001_r8
231 real(r8), parameter :: strat_min = 0.1_r8
232
233 real(r8) :: cff, cff1, cff2, cff3, cff4, dife, difx
234 real(r8) :: tl_cff, tl_cff1, tl_cff2, tl_cff3, tl_cff4
235
236 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: LapT
237
238 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: tl_LapT
239
240 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FE
241 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FX
242
243 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_FE
244 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_FX
245
246 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: FS
247 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dRde
248 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dRdx
249 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dTde
250 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dTdr
251 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dTdx
252
253 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_FS
254 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_dRde
255 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_dRdx
256 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_dTde
257 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_dTdr
258 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: tl_dTdx
259
260#include "set_bounds.h"
261!
262!-----------------------------------------------------------------------
263! Compute tangent linear horizontal biharmonic diffusion along
264! isopycnic surfaces. The biharmonic operator is computed by
265! applying the harmonic operator twice.
266!-----------------------------------------------------------------------
267!
268! Set local I- and J-ranges.
269!
270 IF (ewperiodic(ng)) THEN
271 imin=istr-1
272 imax=iend+1
273 ELSE
274 imin=max(istr-1,1)
275 imax=min(iend+1,lm(ng))
276 END IF
277 IF (nsperiodic(ng)) THEN
278 jmin=jstr-1
279 jmax=jend+1
280 ELSE
281 jmin=max(jstr-1,1)
282 jmax=min(jend+1,mm(ng))
283 END IF
284!
285! Compute horizontal and density gradients. Notice the recursive
286! blocking sequence. The vertical placement of the gradients is:
287!
288! dTdx,dTde(:,:,k1) k rho-points
289! dTdx,dTde(:,:,k2) k+1 rho-points
290! FS,dTdr(:,:,k1) k-1/2 W-points
291! FS,dTdr(:,:,k2) k+1/2 W-points
292!
293#ifdef TS_MIX_STABILITY
294! In order to increase stability, the biharmonic operator is applied
295! as: 3/4 t(:,:,:,nrhs,:) + 1/4 t(:,:,:,nstp,:).
296!
297#endif
298
299 t_loop : DO itrc=1,nt(ng)
300 k2=1
301 k_loop1 : DO k=0,n(ng)
302 k1=k2
303 k2=3-k1
304 IF (k.lt.n(ng)) THEN
305 DO j=jmin,jmax
306 DO i=imin,imax+1
307 cff=0.5_r8*(pm(i,j)+pm(i-1,j))
308#ifdef MASKING
309 cff=cff*umask(i,j)
310#endif
311#ifdef WET_DRY_NOT_YET
312 cff=cff*umask_wet(i,j)
313#endif
314 drdx(i,j,k2)=cff*(pden(i ,j,k+1)- &
315 & pden(i-1,j,k+1))
316 tl_drdx(i,j,k2)=cff*(tl_pden(i ,j,k+1)- &
317 & tl_pden(i-1,j,k+1))
318#if defined TS_MIX_STABILITY
319 dtdx(i,j,k2)=cff*(0.75_r8*(t(i ,j,k+1,nrhs,itrc)- &
320 & t(i-1,j,k+1,nrhs,itrc))+ &
321 & 0.25_r8*(t(i ,j,k+1,nstp,itrc)- &
322 & t(i-1,j,k+1,nstp,itrc)))
323 tl_dtdx(i,j,k2)=cff* &
324 & (0.75_r8*(tl_t(i ,j,k+1,nrhs,itrc)- &
325 & tl_t(i-1,j,k+1,nrhs,itrc))+ &
326 & 0.25_r8*(tl_t(i ,j,k+1,nstp,itrc)- &
327 & tl_t(i-1,j,k+1,nstp,itrc)))
328#elif defined TS_MIX_CLIMA
329 IF (ltracerclm(itrc,ng)) THEN
330 dtdx(i,j,k2)=cff*((t(i ,j,k+1,nrhs,itrc)- &
331 & tclm(i ,j,k+1,itrc))- &
332 & (t(i-1,j,k+1,nrhs,itrc)- &
333 & tclm(i-1,j,k+1,itrc)))
334 ELSE
335 dtdx(i,j,k2)=cff*(t(i ,j,k+1,nrhs,itrc)- &
336 & t(i-1,j,k+1,nrhs,itrc))
337 END IF
338 tl_dtdx(i,j,k2)=cff*(tl_t(i ,j,k+1,nrhs,itrc)- &
339 & tl_t(i-1,j,k+1,nrhs,itrc))
340#else
341 dtdx(i,j,k2)=cff*(t(i ,j,k+1,nrhs,itrc)- &
342 & t(i-1,j,k+1,nrhs,itrc))
343 tl_dtdx(i,j,k2)=cff*(tl_t(i ,j,k+1,nrhs,itrc)- &
344 & tl_t(i-1,j,k+1,nrhs,itrc))
345#endif
346 END DO
347 END DO
348 DO j=jmin,jmax+1
349 DO i=imin,imax
350 cff=0.5_r8*(pn(i,j)+pn(i,j-1))
351#ifdef MASKING
352 cff=cff*vmask(i,j)
353#endif
354#ifdef WET_DRY_NOT_YET
355 cff=cff*vmask_wet(i,j)
356#endif
357 drde(i,j,k2)=cff*(pden(i,j ,k+1)- &
358 & pden(i,j-1,k+1))
359 tl_drde(i,j,k2)=cff*(tl_pden(i,j ,k+1)- &
360 & tl_pden(i,j-1,k+1))
361#if defined TS_MIX_STABILITY
362 dtde(i,j,k2)=cff*(0.75_r8*(t(i,j ,k+1,nrhs,itrc)- &
363 & t(i,j-1,k+1,nrhs,itrc))+ &
364 & 0.25_r8*(t(i,j ,k+1,nstp,itrc)- &
365 & t(i,j-1,k+1,nstp,itrc)))
366 tl_dtde(i,j,k2)=cff* &
367 & (0.75_r8*(tl_t(i,j ,k+1,nrhs,itrc)- &
368 & tl_t(i,j-1,k+1,nrhs,itrc))+ &
369 & 0.25_r8*(tl_t(i,j ,k+1,nstp,itrc)- &
370 & tl_t(i,j-1,k+1,nstp,itrc)))
371#elif defined TS_MIX_CLIMA
372 IF (ltracerclm(itrc,ng)) THEN
373 dtde(i,j,k2)=cff*((t(i,j ,k+1,nrhs,itrc)- &
374 & tclm(i,j ,k+1,itrc))- &
375 & (t(i,j-1,k+1,nrhs,itrc)- &
376 & tclm(i,j-1,k+1,itrc)))
377 ELSE
378 dtde(i,j,k2)=cff*(t(i,j ,k+1,nrhs,itrc)- &
379 & t(i,j-1,k+1,nrhs,itrc))
380 END IF
381 tl_dtde(i,j,k2)=cff*(tl_t(i,j ,k+1,nrhs,itrc)- &
382 & tl_t(i,j-1,k+1,nrhs,itrc))
383#else
384 dtde(i,j,k2)=cff*(t(i,j ,k+1,nrhs,itrc)- &
385 & t(i,j-1,k+1,nrhs,itrc))
386 tl_dtde(i,j,k2)=cff*(tl_t(i,j ,k+1,nrhs,itrc)- &
387 & tl_t(i,j-1,k+1,nrhs,itrc))
388#endif
389 END DO
390 END DO
391 END IF
392 IF ((k.eq.0).or.(k.eq.n(ng))) THEN
393 DO j=jmin-1,jmax+1
394 DO i=imin-1,imax+1
395 dtdr(i,j,k2)=0.0_r8
396 tl_dtdr(i,j,k2)=0.0_r8
397 fs(i,j,k2)=0.0_r8
398 tl_fs(i,j,k2)=0.0_r8
399 END DO
400 END DO
401 ELSE
402 DO j=jmin-1,jmax+1
403 DO i=imin-1,imax+1
404#if defined TS_MIX_MAX_SLOPE
405 cff1=sqrt(drdx(i,j,k2)**2+drdx(i+1,j,k2)**2+ &
406 & drdx(i,j,k1)**2+drdx(i+1,j,k1)**2+ &
407 & drde(i,j,k2)**2+drde(i,j+1,k2)**2+ &
408 & drde(i,j,k1)**2+drde(i,j+1,k1)**2)
409 IF (cff1.ne.0.0_r8) THEN
410 tl_cff1=(drdx(i ,j,k2)*tl_drdx(i ,j,k2)+ &
411 & drdx(i+1,j,k2)*tl_drdx(i+1,j,k2)+ &
412 & drdx(i ,j,k1)*tl_drdx(i ,j,k1)+ &
413 & drdx(i+1,j,k1)*tl_drdx(i+1,j,k1)+ &
414 & drde(i,j ,k2)*tl_drde(i,j ,k2)+ &
415 & drde(i,j+1,k2)*tl_drde(i,j+1,k2)+ &
416 & drde(i,j ,k1)*tl_drde(i,j ,k1)+ &
417 & drde(i,j+1,k1)*tl_drde(i,j+1,k1))/cff1
418 ELSE
419 tl_cff1=0.0_r8
420 END IF
421 cff2=0.25_r8*slope_max* &
422 & (z_r(i,j,k+1)-z_r(i,j,k))*cff1
423 tl_cff2=0.25_r8*slope_max* &
424 & ((tl_z_r(i,j,k+1)-tl_z_r(i,j,k))*cff1+ &
425 & (z_r(i,j,k+1)-z_r(i,j,k))*tl_cff1)
426 cff3=max(pden(i,j,k)-pden(i,j,k+1),small)
427 tl_cff3=(0.5_r8+sign(0.5_r8,pden(i,j,k)-pden(i,j,k+1)- &
428 & small))* &
429 & (tl_pden(i,j,k)-tl_pden(i,j,k+1))
430 cff4=max(cff2,cff3)
431 tl_cff4=(0.5_r8+sign(0.5_r8,cff2-cff3))*tl_cff2+ &
432 & (0.5_r8-sign(0.5_r8,cff2-cff3))*tl_cff3
433 cff=-1.0_r8/cff4
434 tl_cff=cff*cff*tl_cff4
435#elif defined TS_MIX_MIN_STRAT
436 cff1=max(pden(i,j,k)-pden(i,j,k+1), &
437 & strat_min*(z_r(i,j,k+1)-z_r(i,j,k)))
438 tl_cff1=(0.5_r8+sign(0.5_r8, &
439 & pden(i,j,k)-pden(i,j,k+1)- &
440 & strat_min*(z_r(i,j,k+1)- &
441 & z_r(i,j,k ))))* &
442 & (tl_pden(i,j,k)-tl_pden(i,j,k+1))+ &
443 & (0.5_r8-sign(0.5_r8, &
444 & pden(i,j,k)-pden(i,j,k+1)- &
445 & strat_min*(z_r(i,j,k+1)- &
446 & z_r(i,j,k ))))* &
447 & (strat_min*(tl_z_r(i,j,k+1)-tl_z_r(i,j,k )))
448 cff=-1.0_r8/cff1
449 tl_cff=cff*cff*tl_cff1
450#else
451 cff1=max(pden(i,j,k)-pden(i,j,k+1),eps)
452 tl_cff1=(0.5_r8+sign(0.5_r8, &
453 & pden(i,j,k)-pden(i,j,k+1)-eps))* &
454 & (tl_pden(i,j,k)-tl_pden(i,j,k+1))
455 cff=-1.0_r8/cff1
456 tl_cff=cff*cff*tl_cff1
457#endif
458#if defined TS_MIX_STABILITY
459 dtdr(i,j,k2)=cff*(0.75_r8*(t(i,j,k+1,nrhs,itrc)- &
460 & t(i,j,k ,nrhs,itrc))+ &
461 & 0.25_r8*(t(i,j,k+1,nstp,itrc)- &
462 & t(i,j,k ,nstp,itrc)))
463 tl_dtdr(i,j,k2)=tl_cff* &
464 & (0.75_r8*(t(i,j,k+1,nrhs,itrc)- &
465 & t(i,j,k ,nrhs,itrc))+ &
466 & 0.25_r8*(t(i,j,k+1,nstp,itrc)- &
467 & t(i,j,k ,nstp,itrc)))+ &
468 & cff*
469 & (0.75_r8*(tl_t(i,j,k+1,nrhs,itrc)- &
470 & tl_t(i,j,k ,nrhs,itrc))+ &
471 & 0.25_r8*(tl_t(i,j,k+1,nstp,itrc)- &
472 & tl_t(i,j,k ,nstp,itrc)))
473#elif defined TS_MIX_CLIMA
474 IF (ltracerclm(itrc,ng)) THEN
475 dtdr(i,j,k2)=cff*((t(i,j,k+1,nrhs,itrc)- &
476 & tclm(i,j,k+1,itrc))- &
477 & (t(i,j,k ,nrhs,itrc)- &
478 & tclm(i,j,k ,itrc)))
479 tl_dtdr(i,j,k2)=tl_cff*((t(i,j,k+1,nrhs,itrc)- &
480 & tclm(i,j,k+1,itrc))- &
481 & (t(i,j,k ,nrhs,itrc)- &
482 & tclm(i,j,k ,itrc)))+ &
483 & cff*(tl_t(i,j,k+1,nrhs,itrc)- &
484 & tl_t(i,j,k ,nrhs,itrc))
485 ELSE
486 dtdr(i,j,k2)=cff*(t(i,j,k+1,nrhs,itrc)- &
487 & t(i,j,k ,nrhs,itrc))
488 tl_dtdr(i,j,k2)=tl_cff*(t(i,j,k+1,nrhs,itrc)- &
489 & t(i,j,k ,nrhs,itrc))+ &
490 & cff*(tl_t(i,j,k+1,nrhs,itrc)- &
491 & tl_t(i,j,k ,nrhs,itrc))
492 END IF
493#else
494 dtdr(i,j,k2)=cff*(t(i,j,k+1,nrhs,itrc)- &
495 & t(i,j,k ,nrhs,itrc))
496 tl_dtdr(i,j,k2)=tl_cff*(t(i,j,k+1,nrhs,itrc)- &
497 & t(i,j,k ,nrhs,itrc))+ &
498 & cff*(tl_t(i,j,k+1,nrhs,itrc)- &
499 & tl_t(i,j,k ,nrhs,itrc))
500#endif
501 fs(i,j,k2)=cff*(z_r(i,j,k+1)- &
502 & z_r(i,j,k ))
503 tl_fs(i,j,k2)=tl_cff*(z_r(i,j,k+1)- &
504 & z_r(i,j,k ))+ &
505 & cff*(tl_z_r(i,j,k+1)- &
506 & tl_z_r(i,j,k ))
507 END DO
508 END DO
509 END IF
510 IF (k.gt.0) THEN
511 DO j=jmin,jmax
512 DO i=imin,imax+1
513#ifdef DIFF_3DCOEF
514# ifdef TS_U3ADV_SPLIT
515 cff=0.5_r8*diff3d_u(i,j,k)*on_u(i,j)
516# else
517 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
518 & on_u(i,j)
519# endif
520#else
521 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
522 & on_u(i,j)
523#endif
524 fx(i,j)=cff* &
525 & (hz(i,j,k)+hz(i-1,j,k))* &
526 & (dtdx(i,j,k1)- &
527 & 0.5_r8*(max(drdx(i,j,k1),0.0_r8)* &
528 & (dtdr(i-1,j,k1)+ &
529 & dtdr(i ,j,k2))+ &
530 & min(drdx(i,j,k1),0.0_r8)* &
531 & (dtdr(i-1,j,k2)+ &
532 & dtdr(i ,j,k1))))
533 tl_fx(i,j)=cff* &
534 & ((tl_hz(i,j,k)+tl_hz(i-1,j,k))* &
535 & (dtdx(i,j,k1)- &
536 & 0.5_r8*(max(drdx(i,j,k1),0.0_r8)* &
537 & (dtdr(i-1,j,k1)+ &
538 & dtdr(i ,j,k2))+ &
539 & min(drdx(i,j,k1),0.0_r8)* &
540 & (dtdr(i-1,j,k2)+ &
541 & dtdr(i ,j,k1))))+ &
542 & (hz(i,j,k)+hz(i-1,j,k))* &
543 & (tl_dtdx(i,j,k1)- &
544 & 0.5_r8*(max(drdx(i,j,k1),0.0_r8)* &
545 & (tl_dtdr(i-1,j,k1)+ &
546 & tl_dtdr(i ,j,k2))+ &
547 & min(drdx(i,j,k1),0.0_r8)* &
548 & (tl_dtdr(i-1,j,k2)+ &
549 & tl_dtdr(i ,j,k1)))- &
550 & 0.5_r8*((0.5_r8+ &
551 & sign(0.5_r8, drdx(i,j,k1)))* &
552 & tl_drdx(i,j,k1)* &
553 & (dtdr(i-1,j,k1)+dtdr(i,j,k2))+ &
554 & (0.5_r8+ &
555 & sign(0.5_r8,-drdx(i,j,k1)))* &
556 & tl_drdx(i,j,k1)* &
557 & (dtdr(i-1,j,k2)+dtdr(i,j,k1)))))
558 END DO
559 END DO
560 DO j=jmin,jmax+1
561 DO i=imin,imax
562#ifdef DIFF_3DCOEF
563# ifdef TS_U3ADV_SPLIT
564 cff=0.5_r8*diff3d_v(i,j,k)*om_v(i,j)
565# else
566 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
567 & om_v(i,j)
568# endif
569#else
570 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
571 & om_v(i,j)
572#endif
573 fe(i,j)=cff* &
574 & (hz(i,j,k)+hz(i,j-1,k))* &
575 & (dtde(i,j,k1)- &
576 & 0.5_r8*(max(drde(i,j,k1),0.0_r8)* &
577 & (dtdr(i,j-1,k1)+ &
578 & dtdr(i,j ,k2))+ &
579 & min(drde(i,j,k1),0.0_r8)* &
580 & (dtdr(i,j-1,k2)+ &
581 & dtdr(i,j ,k1))))
582 tl_fe(i,j)=cff* &
583 & ((tl_hz(i,j,k)+tl_hz(i,j-1,k))* &
584 & (dtde(i,j,k1)- &
585 & 0.5_r8*(max(drde(i,j,k1),0.0_r8)* &
586 & (dtdr(i,j-1,k1)+ &
587 & dtdr(i,j ,k2))+ &
588 & min(drde(i,j,k1),0.0_r8)* &
589 & (dtdr(i,j-1,k2)+ &
590 & dtdr(i,j ,k1))))+ &
591 & (hz(i,j,k)+hz(i,j-1,k))* &
592 & (tl_dtde(i,j,k1)- &
593 & 0.5_r8*(max(drde(i,j,k1),0.0_r8)* &
594 & (tl_dtdr(i,j-1,k1)+ &
595 & tl_dtdr(i,j ,k2))+ &
596 & min(drde(i,j,k1),0.0_r8)* &
597 & (tl_dtdr(i,j-1,k2)+ &
598 & tl_dtdr(i,j ,k1)))- &
599 & 0.5_r8*((0.5_r8+ &
600 & sign(0.5_r8, drde(i,j,k1)))* &
601 & tl_drde(i,j,k1)* &
602 & (dtdr(i,j-1,k1)+dtdr(i,j,k2))+ &
603 & (0.5_r8+ &
604 & sign(0.5_r8,-drde(i,j,k1)))* &
605 & tl_drde(i,j,k1)* &
606 & (dtdr(i,j-1,k2)+dtdr(i,j,k1)))))
607 END DO
608 END DO
609 IF (k.lt.n(ng)) THEN
610 DO j=jmin,jmax
611 DO i=imin,imax
612#ifdef DIFF_3DCOEF
613# ifdef TS_U3ADV_SPLIT
614 difx=0.125_r8*(diff3d_u(i,j,k )+diff3d_u(i+1,j,k )+ &
615 & diff3d_u(i,j,k+1)+diff3d_u(i+1,j,k+1))
616 dife=0.125_r8*(diff3d_v(i,j,k )+diff3d_v(i,j+1,k )+ &
617 & diff3d_v(i,j,k+1)+diff3d_v(i,j+1,k+1))
618# else
619 difx=0.5_r8*diff3d_r(i,j,k)
620 dife=difx
621# endif
622#else
623 difx=0.5_r8*diff4(i,j,itrc)
624 dife=difx
625#endif
626 cff1=max(drdx(i ,j,k1),0.0_r8)
627 cff2=max(drdx(i+1,j,k2),0.0_r8)
628 cff3=min(drdx(i ,j,k2),0.0_r8)
629 cff4=min(drdx(i+1,j,k1),0.0_r8)
630 tl_cff1=(0.5_r8+sign(0.5_r8, drdx(i ,j,k1)))* &
631 & tl_drdx(i ,j,k1)
632 tl_cff2=(0.5_r8+sign(0.5_r8, drdx(i+1,j,k2)))* &
633 & tl_drdx(i+1,j,k2)
634 tl_cff3=(0.5_r8+sign(0.5_r8,-drdx(i ,j,k2)))* &
635 & tl_drdx(i ,j,k2)
636 tl_cff4=(0.5_r8+sign(0.5_r8,-drdx(i+1,j,k1)))* &
637 & tl_drdx(i+1,j,k1)
638 cff=difx* &
639 & (cff1*(cff1*dtdr(i,j,k2)-dtdx(i ,j,k1))+ &
640 & cff2*(cff2*dtdr(i,j,k2)-dtdx(i+1,j,k2))+ &
641 & cff3*(cff3*dtdr(i,j,k2)-dtdx(i ,j,k2))+ &
642 & cff4*(cff4*dtdr(i,j,k2)-dtdx(i+1,j,k1)))
643 tl_cff=difx* &
644 & (tl_cff1*(cff1*dtdr(i ,j,k2)- &
645 & dtdx(i ,j,k1))+ &
646 & tl_cff2*(cff2*dtdr(i,j,k2)- &
647 & dtdx(i+1,j,k2))+ &
648 & tl_cff3*(cff3*dtdr(i,j,k2)- &
649 & dtdx(i ,j,k2))+ &
650 & tl_cff4*(cff4*dtdr(i,j,k2)- &
651 & dtdx(i+1,j,k1))+ &
652 & cff1*(tl_cff1*dtdr(i,j,k2)+ &
653 & cff1*tl_dtdr(i,j,k2)- &
654 & tl_dtdx(i ,j,k1))+ &
655 & cff2*(tl_cff2*dtdr(i,j,k2)+ &
656 & cff2*tl_dtdr(i,j,k2)- &
657 & tl_dtdx(i+1,j,k2))+ &
658 & cff3*(tl_cff3*dtdr(i,j,k2)+ &
659 & cff3*tl_dtdr(i,j,k2)- &
660 & tl_dtdx(i ,j,k2))+ &
661 & cff4*(tl_cff4*dtdr(i,j,k2)+ &
662 & cff4*tl_dtdr(i,j,k2)- &
663 & tl_dtdx(i+1,j,k1)))
664 cff1=max(drde(i,j ,k1),0.0_r8)
665 cff2=max(drde(i,j+1,k2),0.0_r8)
666 cff3=min(drde(i,j ,k2),0.0_r8)
667 cff4=min(drde(i,j+1,k1),0.0_r8)
668 tl_cff1=(0.5_r8+sign(0.5_r8, drde(i,j ,k1)))* &
669 & tl_drde(i,j ,k1)
670 tl_cff2=(0.5_r8+sign(0.5_r8, drde(i,j+1,k2)))* &
671 & tl_drde(i,j+1,k2)
672 tl_cff3=(0.5_r8+sign(0.5_r8,-drde(i,j ,k2)))* &
673 & tl_drde(i,j ,k2)
674 tl_cff4=(0.5_r8+sign(0.5_r8,-drde(i,j+1,k1)))* &
675 & tl_drde(i,j+1,k1)
676 cff=cff+ &
677 & dife* &
678 & (cff1*(cff1*dtdr(i,j,k2)-dtde(i,j ,k1))+ &
679 & cff2*(cff2*dtdr(i,j,k2)-dtde(i,j+1,k2))+ &
680 & cff3*(cff3*dtdr(i,j,k2)-dtde(i,j ,k2))+ &
681 & cff4*(cff4*dtdr(i,j,k2)-dtde(i,j+1,k1)))
682 tl_cff=tl_cff+ &
683 & dife* &
684 & (tl_cff1*(cff1*dtdr(i,j,k2)- &
685 & dtde(i,j ,k1))+ &
686 & tl_cff2*(cff2*dtdr(i,j,k2)- &
687 & dtde(i,j+1,k2))+ &
688 & tl_cff3*(cff3*dtdr(i,j,k2)- &
689 & dtde(i,j ,k2))+ &
690 & tl_cff4*(cff4*dtdr(i,j,k2)- &
691 & dtde(i,j+1,k1))+ &
692 & cff1*(tl_cff1*dtdr(i,j,k2)+ &
693 & cff1*tl_dtdr(i,j,k2)- &
694 & tl_dtde(i,j ,k1))+ &
695 & cff2*(tl_cff2*dtdr(i,j,k2)+ &
696 & cff2*tl_dtdr(i,j,k2)- &
697 & tl_dtde(i,j+1,k2))+ &
698 & cff3*(tl_cff3*dtdr(i,j,k2)+ &
699 & cff3*tl_dtdr(i,j,k2)- &
700 & tl_dtde(i,j ,k2))+ &
701 & cff4*(tl_cff4*dtdr(i,j,k2)+ &
702 & cff4*tl_dtdr(i,j,k2)- &
703 & tl_dtde(i,j+1,k1)))
704!^ FS(i,j,k2)=cff*FS(i,j,k2) ! recursive
705!^ ! compute after TLM
706!^
707 tl_fs(i,j,k2)=tl_cff*fs(i,j,k2)+ &
708 & cff*tl_fs(i,j,k2)
709 fs(i,j,k2)=cff*fs(i,j,k2) ! recursive
710 END DO
711 END DO
712 END IF
713!
714! Compute first harmonic operator, without mixing coefficient.
715! Multiply by the metrics of the second harmonic operator. Save
716! into work array "LapT".
717!
718 DO j=jmin,jmax
719 DO i=imin,imax
720 cff=pm(i,j)*pn(i,j)
721 cff1=1.0_r8/hz(i,j,k)
722 tl_cff1=-cff1*cff1*tl_hz(i,j,k)
723 lapt(i,j,k)=cff1*(cff* &
724 & (fx(i+1,j)-fx(i,j)+ &
725 & fe(i,j+1)-fe(i,j))+ &
726 & (fs(i,j,k2)-fs(i,j,k1)))
727 tl_lapt(i,j,k)=tl_cff1*(cff* &
728 & (fx(i+1,j)-fx(i,j)+ &
729 & fe(i,j+1)-fe(i,j))+ &
730 & (fs(i,j,k2)-fs(i,j,k1)))+ &
731 & cff1*(cff* &
732 & (tl_fx(i+1,j)-tl_fx(i,j)+ &
733 & tl_fe(i,j+1)-tl_fe(i,j))+ &
734 & (tl_fs(i,j,k2)-tl_fs(i,j,k1)))
735 END DO
736 END DO
737 END IF
738 END DO k_loop1
739!
740! Apply boundary conditions (except periodic; closed or gradient)
741! to the first harmonic operator.
742!
743 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
744 IF (domain(ng)%Western_Edge(tile)) THEN
745 IF (tl_lbc(iwest,istvar(itrc),ng)%closed) THEN
746 DO k=1,n(ng)
747 DO j=jmin,jmax
748 lapt(istr-1,j,k)=0.0_r8
749 tl_lapt(istr-1,j,k)=0.0_r8
750 END DO
751 END DO
752 ELSE
753 DO k=1,n(ng)
754 DO j=jmin,jmax
755 lapt(istr-1,j,k)=lapt(istr,j,k)
756 tl_lapt(istr-1,j,k)=tl_lapt(istr,j,k)
757 END DO
758 END DO
759 END IF
760 END IF
761 END IF
762!
763 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
764 IF (domain(ng)%Eastern_Edge(tile)) THEN
765 IF (tl_lbc(ieast,istvar(itrc),ng)%closed) THEN
766 DO k=1,n(ng)
767 DO j=jmin,jmax
768 lapt(iend+1,j,k)=0.0_r8
769 tl_lapt(iend+1,j,k)=0.0_r8
770 END DO
771 END DO
772 ELSE
773 DO k=1,n(ng)
774 DO j=jmin,jmax
775 lapt(iend+1,j,k)=lapt(iend,j,k)
776 tl_lapt(iend+1,j,k)=tl_lapt(iend,j,k)
777 END DO
778 END DO
779 END IF
780 END IF
781 END IF
782!
783 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
784 IF (domain(ng)%Southern_Edge(tile)) THEN
785 IF (tl_lbc(isouth,istvar(itrc),ng)%closed) THEN
786 DO k=1,n(ng)
787 DO i=imin,imax
788 lapt(i,jstr-1,k)=0.0_r8
789 tl_lapt(i,jstr-1,k)=0.0_r8
790 END DO
791 END DO
792 ELSE
793 DO k=1,n(ng)
794 DO i=imin,imax
795 lapt(i,jstr-1,k)=lapt(i,jstr,k)
796 tl_lapt(i,jstr-1,k)=tl_lapt(i,jstr,k)
797 END DO
798 END DO
799 END IF
800 END IF
801 END IF
802!
803 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
804 IF (domain(ng)%Northern_Edge(tile)) THEN
805 IF (tl_lbc(inorth,istvar(itrc),ng)%closed) THEN
806 DO k=1,n(ng)
807 DO i=imin,imax
808 lapt(i,jend+1,k)=0.0_r8
809 tl_lapt(i,jend+1,k)=0.0_r8
810 END DO
811 END DO
812 ELSE
813 DO k=1,n(ng)
814 DO i=imin,imax
815 lapt(i,jend+1,k)=lapt(i,jend,k)
816 tl_lapt(i,jend+1,k)=tl_lapt(i,jend,k)
817 END DO
818 END DO
819 END IF
820 END IF
821 END IF
822!
823 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
824 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
825 IF (domain(ng)%SouthWest_Corner(tile)) THEN
826 DO k=1,n(ng)
827 lapt(istr-1,jstr-1,k)=0.5_r8* &
828 & (lapt(istr ,jstr-1,k)+ &
829 & lapt(istr-1,jstr ,k))
830 tl_lapt(istr-1,jstr-1,k)=0.5_r8* &
831 & (tl_lapt(istr ,jstr-1,k)+ &
832 tl_lapt(istr-1,jstr ,k))
833 END DO
834 END IF
835 END IF
836
837 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
838 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
839 IF (domain(ng)%SouthEast_Corner(tile)) THEN
840 DO k=1,n(ng)
841 lapt(iend+1,jstr-1,k)=0.5_r8* &
842 & (lapt(iend ,jstr-1,k)+ &
843 & lapt(iend+1,jstr ,k))
844 tl_lapt(iend+1,jstr-1,k)=0.5_r8* &
845 & (tl_lapt(iend ,jstr-1,k)+ &
846 & tl_lapt(iend+1,jstr ,k))
847 END DO
848 END IF
849 END IF
850
851 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
852 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
853 IF (domain(ng)%NorthWest_Corner(tile)) THEN
854 DO k=1,n(ng)
855 lapt(istr-1,jend+1,k)=0.5_r8* &
856 & (lapt(istr ,jend+1,k)+ &
857 & lapt(istr-1,jend ,k))
858 tl_lapt(istr-1,jend+1,k)=0.5_r8* &
859 & (tl_lapt(istr ,jend+1,k)+ &
860 & tl_lapt(istr-1,jend ,k))
861 END DO
862 END IF
863 END IF
864
865 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
866 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
867 IF (domain(ng)%NorthEast_Corner(tile)) THEN
868 DO k=1,n(ng)
869 lapt(iend+1,jend+1,k)=0.5_r8* &
870 & (lapt(iend ,jend+1,k)+ &
871 & lapt(iend+1,jend ,k))
872 tl_lapt(iend+1,jend+1,k)=0.5_r8* &
873 & (tl_lapt(iend ,jend+1,k)+ &
874 & tl_lapt(iend+1,jend ,k))
875 END DO
876 END IF
877 END IF
878!
879! Compute horizontal and density gradients associated with the
880! second rotated harmonic operator.
881!
882 k2=1
883 k_loop2: DO k=0,n(ng)
884 k1=k2
885 k2=3-k1
886 IF (k.lt.n(ng)) THEN
887 DO j=jstr,jend
888 DO i=istr,iend+1
889 cff=0.5_r8*(pm(i,j)+pm(i-1,j))
890#ifdef MASKING
891 cff=cff*umask(i,j)
892#endif
893#ifdef WET_DRY_NOT_YET
894 cff=cff*umask_wet(i,j)
895#endif
896 drdx(i,j,k2)=cff*(pden(i ,j,k+1)- &
897 & pden(i-1,j,k+1))
898 tl_drdx(i,j,k2)=cff*(tl_pden(i ,j,k+1)- &
899 & tl_pden(i-1,j,k+1))
900 dtdx(i,j,k2)=cff*(lapt(i ,j,k+1)- &
901 & lapt(i-1,j,k+1))
902 tl_dtdx(i,j,k2)=cff*(tl_lapt(i ,j,k+1)- &
903 & tl_lapt(i-1,j,k+1))
904 END DO
905 END DO
906 DO j=jstr,jend+1
907 DO i=istr,iend
908 cff=0.5_r8*(pn(i,j)+pn(i,j-1))
909#ifdef MASKING
910 cff=cff*vmask(i,j)
911#endif
912#ifdef WET_DRY_NOT_YET
913 cff=cff*vmask_wet(i,j)
914#endif
915 drde(i,j,k2)=cff*(pden(i,j ,k+1)- &
916 & pden(i,j-1,k+1))
917 tl_drde(i,j,k2)=cff*(tl_pden(i,j ,k+1)- &
918 & tl_pden(i,j-1,k+1))
919 dtde(i,j,k2)=cff*(lapt(i,j ,k+1)- &
920 & lapt(i,j-1,k+1))
921 tl_dtde(i,j,k2)=cff*(tl_lapt(i,j ,k+1)- &
922 & tl_lapt(i,j-1,k+1))
923 END DO
924 END DO
925 END IF
926 IF ((k.eq.0).or.(k.eq.n(ng))) THEN
927 DO j=jstr-1,jend+1
928 DO i=istr-1,iend+1
929 dtdr(i,j,k2)=0.0_r8
930 tl_dtdr(i,j,k2)=0.0_r8
931 fs(i,j,k2)=0.0_r8
932 tl_fs(i,j,k2)=0.0_r8
933 END DO
934 END DO
935 ELSE
936 DO j=jstr-1,jend+1
937 DO i=istr-1,iend+1
938#if defined TS_MIX_MAX_SLOPE
939 cff1=sqrt(drdx(i,j,k2)**2+drdx(i+1,j,k2)**2+ &
940 & drdx(i,j,k1)**2+drdx(i+1,j,k1)**2+ &
941 & drde(i,j,k2)**2+drde(i,j+1,k2)**2+ &
942 & drde(i,j,k1)**2+drde(i,j+1,k1)**2)
943 IF (cff1.ne.0.0_r8) THEN
944 tl_cff1=(drdx(i ,j,k2)*tl_drdx(i ,j,k2)+ &
945 & drdx(i+1,j,k2)*tl_drdx(i+1,j,k2)+ &
946 & drdx(i ,j,k1)*tl_drdx(i ,j,k1)+ &
947 & drdx(i+1,j,k1)*tl_drdx(i+1,j,k1)+ &
948 & drde(i,j ,k2)*tl_drde(i,j ,k2)+ &
949 & drde(i,j+1,k2)*tl_drde(i,j+1,k2)+ &
950 & drde(i,j ,k1)*tl_drde(i,j ,k1)+ &
951 & drde(i,j+1,k1)*tl_drde(i,j+1,k1))/cff1
952 ELSE
953 tl_cff1=0.0_r8
954 END IF
955 cff2=0.25_r8*slope_max* &
956 & (z_r(i,j,k+1)-z_r(i,j,k))*cff1
957 tl_cff2=0.25_r8*slope_max* &
958 & ((tl_z_r(i,j,k+1)-tl_z_r(i,j,k))*cff1+ &
959 & (z_r(i,j,k+1)-z_r(i,j,k))*tl_cff1)
960 cff3=max(pden(i,j,k)-pden(i,j,k+1),small)
961 tl_cff3=(0.5_r8+sign(0.5_r8,pden(i,j,k)-pden(i,j,k+1)- &
962 & small))* &
963 & (tl_pden(i,j,k)-tl_pden(i,j,k+1))
964 cff4=max(cff2,cff3)
965 tl_cff4=(0.5_r8+sign(0.5_r8,cff2-cff3))*tl_cff2+ &
966 & (0.5_r8-sign(0.5_r8,cff2-cff3))*tl_cff3
967 cff=-1.0_r8/cff4
968 tl_cff=cff*cff*tl_cff4
969#elif defined TS_MIX_MIN_STRAT
970 cff1=max(pden(i,j,k)-pden(i,j,k+1), &
971 & strat_min*(z_r(i,j,k+1)-z_r(i,j,k)))
972 tl_cff1=(0.5_r8+sign(0.5_r8, &
973 & pden(i,j,k)-pden(i,j,k+1)- &
974 & strat_min*(z_r(i,j,k+1)- &
975 & z_r(i,j,k ))))* &
976 & (tl_pden(i,j,k)-tl_pden(i,j,k+1))+ &
977 & (0.5_r8-sign(0.5_r8, &
978 & pden(i,j,k)-pden(i,j,k+1)- &
979 & strat_min*(z_r(i,j,k+1)- &
980 & z_r(i,j,k ))))* &
981 & (strat_min*(tl_z_r(i,j,k+1)-tl_z_r(i,j,k )))
982 cff=-1.0_r8/cff1
983 tl_cff=cff*cff*tl_cff1
984#else
985 cff1=max(pden(i,j,k)-pden(i,j,k+1),eps)
986 tl_cff1=(0.5_r8+sign(0.5_r8, &
987 & pden(i,j,k)-pden(i,j,k+1)-eps))* &
988 & (tl_pden(i,j,k)-tl_pden(i,j,k+1))
989 cff=-1.0_r8/cff1
990 tl_cff=cff*cff*tl_cff1
991#endif
992 dtdr(i,j,k2)=cff*(lapt(i,j,k+1)- &
993 & lapt(i,j,k ))
994 tl_dtdr(i,j,k2)=tl_cff*(lapt(i,j,k+1)- &
995 & lapt(i,j,k ))+ &
996 & cff*(tl_lapt(i,j,k+1)- &
997 & tl_lapt(i,j,k ))
998 fs(i,j,k2)=cff*(z_r(i,j,k+1)- &
999 & z_r(i,j,k ))
1000 tl_fs(i,j,k2)=tl_cff*(z_r(i,j,k+1)- &
1001 & z_r(i,j,k ))+ &
1002 & cff*(tl_z_r(i,j,k+1)- &
1003 & tl_z_r(i,j,k ))
1004 END DO
1005 END DO
1006 END IF
1007!
1008! Compute components of the rotated tracer flux (T m4/s) along
1009! isopycnic surfaces.
1010!
1011 IF (k.gt.0) THEN
1012 DO j=jstr,jend
1013 DO i=istr,iend+1
1014#ifdef DIFF_3DCOEF
1015# ifdef TS_U3ADV_SPLIT
1016 cff=0.5_r8*diff3d_u(i,j,k)*on_u(i,j)
1017# else
1018 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
1019 & on_u(i,j)
1020# endif
1021#else
1022 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
1023 & on_u(i,j)
1024#endif
1025!^ FX(i,j)=cff* &
1026!^ & (Hz(i,j,k)+Hz(i-1,j,k))* &
1027!^ & (dTdx(i,j,k1)- &
1028!^ & 0.5_r8*(MAX(dRdx(i,j,k1),0.0_r8)* &
1029!^ & (dTdr(i-1,j,k1)+ &
1030!^ & dTdr(i ,j,k2))+ &
1031!^ & MIN(dRdx(i,j,k1),0.0_r8)* &
1032!^ & (dTdr(i-1,j,k2)+ &
1033!^ & dTdr(i ,j,k1))))
1034!^
1035 tl_fx(i,j)=cff* &
1036 & ((tl_hz(i,j,k)+tl_hz(i-1,j,k))* &
1037 & (dtdx(i,j,k1)- &
1038 & 0.5_r8*(max(drdx(i,j,k1),0.0_r8)* &
1039 & (dtdr(i-1,j,k1)+ &
1040 & dtdr(i ,j,k2))+ &
1041 & min(drdx(i,j,k1),0.0_r8)* &
1042 & (dtdr(i-1,j,k2)+ &
1043 & dtdr(i ,j,k1))))+ &
1044 & (hz(i,j,k)+hz(i-1,j,k))* &
1045 & (tl_dtdx(i,j,k1)- &
1046 & 0.5_r8*(max(drdx(i,j,k1),0.0_r8)* &
1047 & (tl_dtdr(i-1,j,k1)+ &
1048 & tl_dtdr(i ,j,k2))+ &
1049 & min(drdx(i,j,k1),0.0_r8)* &
1050 & (tl_dtdr(i-1,j,k2)+ &
1051 & tl_dtdr(i ,j,k1)))- &
1052 & 0.5_r8*((0.5_r8+ &
1053 & sign(0.5_r8, drdx(i,j,k1)))* &
1054 & tl_drdx(i,j,k1)* &
1055 & (dtdr(i-1,j,k1)+dtdr(i,j,k2))+ &
1056 & (0.5_r8+ &
1057 & sign(0.5_r8,-drdx(i,j,k1)))* &
1058 & tl_drdx(i,j,k1)* &
1059 & (dtdr(i-1,j,k2)+dtdr(i,j,k1)))))
1060 END DO
1061 END DO
1062 DO j=jstr,jend+1
1063 DO i=istr,iend
1064#ifdef DIFF_3DCOEF
1065# ifdef TS_U3ADV_SPLIT
1066 cff=0.5_r8*diff3d_v(i,j,k)*om_v(i,j)
1067# else
1068 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
1069 & om_v(i,j)
1070# endif
1071#else
1072 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
1073 & om_v(i,j)
1074#endif
1075!^ FE(i,j)=cff* &
1076!^ & (Hz(i,j,k)+Hz(i,j-1,k))* &
1077!^ & (dTde(i,j,k1)- &
1078!^ & 0.5_r8*(MAX(dRde(i,j,k1),0.0_r8)* &
1079!^ & (dTdr(i,j-1,k1)+ &
1080!^ & dTdr(i,j ,k2))+ &
1081!^ & MIN(dRde(i,j,k1),0.0_r8)* &
1082!^ & (dTdr(i,j-1,k2)+ &
1083!^ & dTdr(i,j ,k1))))
1084!^
1085 tl_fe(i,j)=cff* &
1086 & ((tl_hz(i,j,k)+tl_hz(i,j-1,k))* &
1087 & (dtde(i,j,k1)- &
1088 & 0.5_r8*(max(drde(i,j,k1),0.0_r8)* &
1089 & (dtdr(i,j-1,k1)+ &
1090 & dtdr(i,j ,k2))+ &
1091 & min(drde(i,j,k1),0.0_r8)* &
1092 & (dtdr(i,j-1,k2)+ &
1093 & dtdr(i,j ,k1))))+ &
1094 & (hz(i,j,k)+hz(i,j-1,k))* &
1095 & (tl_dtde(i,j,k1)- &
1096 & 0.5_r8*(max(drde(i,j,k1),0.0_r8)* &
1097 & (tl_dtdr(i,j-1,k1)+ &
1098 & tl_dtdr(i,j ,k2))+ &
1099 & min(drde(i,j,k1),0.0_r8)* &
1100 & (tl_dtdr(i,j-1,k2)+ &
1101 & tl_dtdr(i,j ,k1)))- &
1102 & 0.5_r8*((0.5_r8+ &
1103 & sign(0.5_r8, drde(i,j,k1)))* &
1104 & tl_drde(i,j,k1)* &
1105 & (dtdr(i,j-1,k1)+dtdr(i,j,k2))+ &
1106 & (0.5_r8+ &
1107 & sign(0.5_r8,-drde(i,j,k1)))* &
1108 & tl_drde(i,j,k1)* &
1109 & (dtdr(i,j-1,k2)+dtdr(i,j,k1)))))
1110 END DO
1111 END DO
1112 IF (k.lt.n(ng)) THEN
1113 DO j=jstr,jend
1114 DO i=istr,iend
1115#ifdef DIFF_3DCOEF
1116# ifdef TS_U3ADV_SPLIT
1117 difx=0.125_r8*(diff3d_u(i,j,k )+diff3d_u(i+1,j,k )+ &
1118 & diff3d_u(i,j,k+1)+diff3d_u(i+1,j,k+1))
1119 dife=0.125_r8*(diff3d_v(i,j,k )+diff3d_v(i,j+1,k )+ &
1120 & diff3d_v(i,j,k+1)+diff3d_v(i,j+1,k+1))
1121# else
1122 difx=0.5_r8*diff3d_r(i,j,k)
1123 dife=difx
1124# endif
1125#else
1126 difx=0.5_r8*diff4(i,j,itrc)
1127 dife=difx
1128#endif
1129 cff1=max(drdx(i ,j,k1),0.0_r8)
1130 cff2=max(drdx(i+1,j,k2),0.0_r8)
1131 cff3=min(drdx(i ,j,k2),0.0_r8)
1132 cff4=min(drdx(i+1,j,k1),0.0_r8)
1133 tl_cff1=(0.5_r8+sign(0.5_r8, drdx(i ,j,k1)))* &
1134 & tl_drdx(i ,j,k1)
1135 tl_cff2=(0.5_r8+sign(0.5_r8, drdx(i+1,j,k2)))* &
1136 & tl_drdx(i+1,j,k2)
1137 tl_cff3=(0.5_r8+sign(0.5_r8,-drdx(i ,j,k2)))* &
1138 & tl_drdx(i ,j,k2)
1139 tl_cff4=(0.5_r8+sign(0.5_r8,-drdx(i+1,j,k1)))* &
1140 & tl_drdx(i+1,j,k1)
1141 cff=difx* &
1142 & (cff1*(cff1*dtdr(i,j,k2)-dtdx(i ,j,k1))+ &
1143 & cff2*(cff2*dtdr(i,j,k2)-dtdx(i+1,j,k2))+ &
1144 & cff3*(cff3*dtdr(i,j,k2)-dtdx(i ,j,k2))+ &
1145 & cff4*(cff4*dtdr(i,j,k2)-dtdx(i+1,j,k1)))
1146 tl_cff=difx* &
1147 & (tl_cff1*(cff1*dtdr(i ,j,k2)- &
1148 & dtdx(i ,j,k1))+ &
1149 & tl_cff2*(cff2*dtdr(i,j,k2)- &
1150 & dtdx(i+1,j,k2))+ &
1151 & tl_cff3*(cff3*dtdr(i,j,k2)- &
1152 & dtdx(i ,j,k2))+ &
1153 & tl_cff4*(cff4*dtdr(i,j,k2)- &
1154 & dtdx(i+1,j,k1))+ &
1155 & cff1*(tl_cff1*dtdr(i,j,k2)+ &
1156 & cff1*tl_dtdr(i,j,k2)- &
1157 & tl_dtdx(i ,j,k1))+ &
1158 & cff2*(tl_cff2*dtdr(i,j,k2)+ &
1159 & cff2*tl_dtdr(i,j,k2)- &
1160 & tl_dtdx(i+1,j,k2))+ &
1161 & cff3*(tl_cff3*dtdr(i,j,k2)+ &
1162 & cff3*tl_dtdr(i,j,k2)- &
1163 & tl_dtdx(i ,j,k2))+ &
1164 & cff4*(tl_cff4*dtdr(i,j,k2)+ &
1165 & cff4*tl_dtdr(i,j,k2)- &
1166 & tl_dtdx(i+1,j,k1)))
1167 cff1=max(drde(i,j ,k1),0.0_r8)
1168 cff2=max(drde(i,j+1,k2),0.0_r8)
1169 cff3=min(drde(i,j ,k2),0.0_r8)
1170 cff4=min(drde(i,j+1,k1),0.0_r8)
1171 tl_cff1=(0.5_r8+sign(0.5_r8, drde(i,j ,k1)))* &
1172 & tl_drde(i,j ,k1)
1173 tl_cff2=(0.5_r8+sign(0.5_r8, drde(i,j+1,k2)))* &
1174 & tl_drde(i,j+1,k2)
1175 tl_cff3=(0.5_r8+sign(0.5_r8,-drde(i,j ,k2)))* &
1176 & tl_drde(i,j ,k2)
1177 tl_cff4=(0.5_r8+sign(0.5_r8,-drde(i,j+1,k1)))* &
1178 & tl_drde(i,j+1,k1)
1179 cff=cff+ &
1180 & dife* &
1181 & (cff1*(cff1*dtdr(i,j,k2)-dtde(i,j ,k1))+ &
1182 & cff2*(cff2*dtdr(i,j,k2)-dtde(i,j+1,k2))+ &
1183 & cff3*(cff3*dtdr(i,j,k2)-dtde(i,j ,k2))+ &
1184 & cff4*(cff4*dtdr(i,j,k2)-dtde(i,j+1,k1)))
1185 tl_cff=tl_cff+ &
1186 & dife* &
1187 & (tl_cff1*(cff1*dtdr(i,j,k2)- &
1188 & dtde(i,j ,k1))+ &
1189 & tl_cff2*(cff2*dtdr(i,j,k2)- &
1190 & dtde(i,j+1,k2))+ &
1191 & tl_cff3*(cff3*dtdr(i,j,k2)- &
1192 & dtde(i,j ,k2))+ &
1193 & tl_cff4*(cff4*dtdr(i,j,k2)- &
1194 & dtde(i,j+1,k1))+ &
1195 & cff1*(tl_cff1*dtdr(i,j,k2)+ &
1196 & cff1*tl_dtdr(i,j,k2)- &
1197 & tl_dtde(i,j ,k1))+ &
1198 & cff2*(tl_cff2*dtdr(i,j,k2)+ &
1199 & cff2*tl_dtdr(i,j,k2)- &
1200 & tl_dtde(i,j+1,k2))+ &
1201 & cff3*(tl_cff3*dtdr(i,j,k2)+ &
1202 & cff3*tl_dtdr(i,j,k2)- &
1203 & tl_dtde(i,j ,k2))+ &
1204 & cff4*(tl_cff4*dtdr(i,j,k2)+ &
1205 & cff4*tl_dtdr(i,j,k2)- &
1206 & tl_dtde(i,j+1,k1)))
1207!^ FS(i,j,k2)=cff*FS(i,j,k2) ! recursive
1208!^ ! compute after TLM
1209!^
1210 tl_fs(i,j,k2)=tl_cff*fs(i,j,k2)+ &
1211 & cff*tl_fs(i,j,k2)
1212 fs(i,j,k2)=cff*fs(i,j,k2) ! recursive
1213 END DO
1214 END DO
1215 END IF
1216!
1217! Time-step biharmonic, isopycnal diffusion term (m Tunits).
1218!
1219 DO j=jstr,jend
1220 DO i=istr,iend
1221!^ cff=dt(ng)*pm(i,j)*pn(i,j)* &
1222!^ & (FX(i+1,j )-FX(i,j)+ &
1223!^ & FE(i ,j+1)-FE(i,j))+ &
1224!^ & dt(ng)*(FS(i,j,k2)-FS(i,j,k1))
1225!^
1226 tl_cff=dt(ng)*pm(i,j)*pn(i,j)* &
1227 & (tl_fx(i+1,j)-tl_fx(i,j)+ &
1228 & tl_fe(i,j+1)-tl_fe(i,j))+ &
1229 & dt(ng)*(tl_fs(i,j,k2)-tl_fs(i,j,k1))
1230!^ t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)-cff
1231!^
1232 tl_t(i,j,k,nnew,itrc)=tl_t(i,j,k,nnew,itrc)-tl_cff
1233#ifdef DIAGNOSTICS_TS
1234!! DiaTwrk(i,j,k,itrc,iThdif)=-cff
1235#endif
1236 END DO
1237 END DO
1238 END IF
1239 END DO k_loop2
1240 END DO t_loop
1241!
1242 RETURN

References mod_scalars::compositegrid, mod_param::domain, mod_scalars::dt, mod_scalars::ewperiodic, mod_scalars::ieast, mod_scalars::inorth, mod_scalars::isouth, mod_ncparam::istvar, mod_scalars::iwest, mod_param::lm, mod_scalars::ltracerclm, mod_param::mm, mod_scalars::nsperiodic, and mod_param::tl_lbc.

◆ tl_t3dmix4_s_tile()

subroutine tl_t3dmix4_mod::tl_t3dmix4_s_tile ( integer, intent(in) ng,
integer, intent(in) tile,
integer, intent(in) lbi,
integer, intent(in) ubi,
integer, intent(in) lbj,
integer, intent(in) ubj,
integer, intent(in) imins,
integer, intent(in) imaxs,
integer, intent(in) jmins,
integer, intent(in) jmaxs,
integer, intent(in) nrhs,
integer, intent(in) nstp,
integer, intent(in) nnew,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) umask,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) vmask,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) umask_wet,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) vmask_wet,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) hz,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) tl_hz,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pmon_u,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pnom_v,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pm,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pn,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) diff3d_u,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) diff3d_v,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) diff3d_r,
real(r8), dimension(lbi:ubi,lbj:ubj,nt(ng)), intent(in) diff4,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),nt(ng)), intent(in) tclm,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),3,nt(ng)), intent(in) t,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),3,nt(ng)), intent(inout) tl_t )
private

Definition at line 96 of file tl_t3dmix4_s.h.

124!***********************************************************************
125!
126 USE mod_param
127 USE mod_ncparam
128 USE mod_scalars
129!
130! Imported variable declarations.
131!
132 integer, intent(in) :: ng, tile
133 integer, intent(in) :: LBi, UBi, LBj, UBj
134 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
135 integer, intent(in) :: nrhs, nstp, nnew
136
137#ifdef ASSUMED_SHAPE
138# ifdef MASKING
139 real(r8), intent(in) :: umask(LBi:,LBj:)
140 real(r8), intent(in) :: vmask(LBi:,LBj:)
141# endif
142# ifdef WET_DRY_NOT_YET
143 real(r8), intent(in) :: umask_wet(LBi:,LBj:)
144 real(r8), intent(in) :: vmask_wet(LBi:,LBj:)
145# endif
146# ifdef DIFF_3DCOEF
147# ifdef TS_U3ADV_SPLIT_NOT_YET
148 real(r8), intent(in) :: diff3d_u(LBi:,LBj:,:)
149 real(r8), intent(in) :: diff3d_v(LBi:,LBj:,:)
150# else
151 real(r8), intent(in) :: diff3d_r(LBi:,LBj:,:)
152# endif
153# else
154 real(r8), intent(in) :: diff4(LBi:,LBj:,:)
155# endif
156 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
157 real(r8), intent(in) :: tl_Hz(LBi:,LBj:,:)
158 real(r8), intent(in) :: pmon_u(LBi:,LBj:)
159 real(r8), intent(in) :: pnom_v(LBi:,LBj:)
160 real(r8), intent(in) :: pm(LBi:,LBj:)
161 real(r8), intent(in) :: pn(LBi:,LBj:)
162 real(r8), intent(in) :: t(LBi:,LBj:,:,:,:)
163# ifdef TS_MIX_CLIMA
164 real(r8), intent(in) :: tclm(LBi:,LBj:,:,:)
165# endif
166# ifdef DIAGNOSTICS_TS
167!! real(r8), intent(inout) :: DiaTwrk(LBi:,LBj:,:,:,:)
168# endif
169 real(r8), intent(inout) :: tl_t(LBi:,LBj:,:,:,:)
170#else
171# ifdef MASKING
172 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
173 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
174# endif
175# ifdef WET_DRY_NOT_YET
176 real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
177 real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
178# endif
179# ifdef DIFF_3DCOEF
180# ifdef TS_U3ADV_SPLIT_NOT_YET
181 real(r8), intent(in) :: diff3d_u(LBi:UBi,LBj:UBj,N(ng))
182 real(r8), intent(in) :: diff3d_v(LBi:UBi,LBj:UBj,N(ng))
183# else
184 real(r8), intent(in) :: diff3d_r(LBi:UBi,LBj:UBj,N(ng))
185# endif
186# else
187 real(r8), intent(in) :: diff4(LBi:UBi,LBj:UBj,NT(ng))
188# endif
189 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
190 real(r8), intent(in) :: tl_Hz(LBi:UBi,LBj:UBj,N(ng))
191 real(r8), intent(in) :: pmon_u(LBi:UBi,LBj:UBj)
192 real(r8), intent(in) :: pnom_v(LBi:UBi,LBj:UBj)
193 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
194 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
195 real(r8), intent(in) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
196# ifdef TS_MIX_CLIMA
197 real(r8), intent(in) :: tclm(LBi:UBi,LBj:UBj,N(ng),NT(ng))
198# endif
199# ifdef DIAGNOSTICS_TS
200!! real(r8), intent(inout) :: DiaTwrk(LBi:UBi,LBj:UBj,N(ng),NT(ng), &
201!! & NDT)
202# endif
203 real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
204#endif
205!
206! Local variable declarations.
207!
208 integer :: Imin, Imax, Jmin, Jmax
209 integer :: i, itrc, j, k
210
211 real(r8) :: cff, cff1, tl_cff, tl_cff1
212
213 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FE
214 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FX
215 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: LapT
216
217 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_FE
218 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_FX
219 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_LapT
220
221#include "set_bounds.h"
222!
223!-----------------------------------------------------------------------
224! Compute horizontal biharmonic diffusion along constant S-surfaces.
225! The biharmonic operator is computed by applying the harmonic
226! operator twice.
227#ifdef TS_MIX_STABILITY
228! In order to increase stability, the biharmonic operator is applied
229! as: 3/4 t(:,:,:,nrhs,:) + 1/4 t(:,:,:,nstp,:).
230#endif
231!-----------------------------------------------------------------------
232!
233! Set local I- and J-ranges.
234!
235 IF (ewperiodic(ng)) THEN
236 imin=istr-1
237 imax=iend+1
238 ELSE
239 imin=max(istr-1,1)
240 imax=min(iend+1,lm(ng))
241 END IF
242 IF (nsperiodic(ng)) THEN
243 jmin=jstr-1
244 jmax=jend+1
245 ELSE
246 jmin=max(jstr-1,1)
247 jmax=min(jend+1,mm(ng))
248 END IF
249!
250! Compute horizontal tracer flux in the XI- and ETA-directions.
251!
252 DO itrc=1,nt(ng)
253 DO k=1,n(ng)
254 DO j=jmin,jmax
255 DO i=imin,imax+1
256#ifdef DIFF_3DCOEF
257# ifdef TS_U3ADV_SPLIT_NOT_YET
258 cff=0.5_r8*diff3d_u(i,j,k)*pmon_u(i,j)
259# else
260 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
261 & pmon_u(i,j)
262# endif
263#else
264 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
265 & pmon_u(i,j)
266#endif
267#ifdef MASKING
268 cff=cff*umask(i,j)
269#endif
270#ifdef WET_DRY_NOT_YET
271 cff=cff*umask_wet(i,j)
272#endif
273#if defined TS_MIX_STABILITY
274 fx(i,j)=cff* &
275 & (hz(i,j,k)+hz(i-1,j,k))* &
276 & (0.75_r8*(t(i ,j,k,nrhs,itrc)- &
277 & t(i-1,j,k,nrhs,itrc))+ &
278 & 0.25_r8*(t(i ,j,k,nstp,itrc)- &
279 & t(i-1,j,k,nstp,itrc)))
280 tl_fx(i,j)=cff* &
281 & ((tl_hz(i,j,k)+tl_hz(i-1,j,k))* &
282 & (0.75_r8*(t(i ,j,k,nrhs,itrc)- &
283 & t(i-1,j,k,nrhs,itrc))+ &
284 & 0.25_r8*(t(i ,j,k,nstp,itrc)- &
285 & t(i-1,j,k,nstp,itrc)))+ &
286 & (hz(i,j,k)+hz(i-1,j,k))* &
287 & (0.75_r8*(tl_t(i ,j,k,nrhs,itrc)- &
288 & tl_t(i-1,j,k,nrhs,itrc))+ &
289 & 0.25_r8*(tl_t(i ,j,k,nstp,itrc)- &
290 & tl_t(i-1,j,k,nstp,itrc))))
291#elif defined TS_MIX_CLIMA
292 IF (ltracerclm(itrc,ng)) THEN
293 fx(i,j)=cff* &
294 & (hz(i,j,k)+hz(i-1,j,k))* &
295 & ((t(i ,j,k,nrhs,itrc)-tclm(i ,j,k,itrc))- &
296 & (t(i-1,j,k,nrhs,itrc)-tclm(i-1,j,k,itrc)))
297 tl_fx(i,j)=cff* &
298 & ((tl_hz(i,j,k)+tl_hz(i-1,j,k))* &
299 & ((t(i ,j,k,nrhs,itrc)- &
300 & tclm(i ,j,k,itrc))- &
301 & (t(i-1,j,k,nrhs,itrc)- &
302 & tclm(i-1,j,k,itrc)))+ &
303 & (hz(i,j,k)+hz(i-1,j,k))* &
304 & (tl_t(i ,j,k,nrhs,itrc)- &
305 & tl_t(i-1,j,k,nrhs,itrc)))
306 ELSE
307 fx(i,j)=cff* &
308 & (hz(i,j,k)+hz(i-1,j,k))* &
309 & (t(i ,j,k,nrhs,itrc)- &
310 & t(i-1,j,k,nrhs,itrc))
311 tl_fx(i,j)=cff* &
312 & ((tl_hz(i,j,k)+tl_hz(i-1,j,k))* &
313 & (t(i ,j,k,nrhs,itrc)- &
314 & t(i-1,j,k,nrhs,itrc))+ &
315 & (hz(i,j,k)+hz(i-1,j,k))* &
316 & (tl_t(i ,j,k,nrhs,itrc)- &
317 & tl_t(i-1,j,k,nrhs,itrc)))
318 END IF
319#else
320 fx(i,j)=cff* &
321 & (hz(i,j,k)+hz(i-1,j,k))* &
322 & (t(i ,j,k,nrhs,itrc)- &
323 & t(i-1,j,k,nrhs,itrc))
324 tl_fx(i,j)=cff* &
325 & ((tl_hz(i,j,k)+tl_hz(i-1,j,k))* &
326 & (t(i ,j,k,nrhs,itrc)- &
327 & t(i-1,j,k,nrhs,itrc))+ &
328 & (hz(i,j,k)+hz(i-1,j,k))* &
329 & (tl_t(i ,j,k,nrhs,itrc)- &
330 & tl_t(i-1,j,k,nrhs,itrc)))
331#endif
332 END DO
333 END DO
334 DO j=jmin,jmax+1
335 DO i=imin,imax
336#ifdef DIFF_3DCOEF
337# ifdef TS_U3ADV_SPLIT_NOT_YET
338 cff=0.5_r8*diff3d_v(i,j,k)*pnom_v(i,j)
339# else
340 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
341 & pnom_v(i,j)
342# endif
343#else
344 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
345 & pnom_v(i,j)
346#endif
347#ifdef MASKING
348 cff=cff*vmask(i,j)
349#endif
350#ifdef WET_DRY_NOT_YET
351 cff=cff*vmask_wet(i,j)
352#endif
353#if defined TS_MIX_STABILITY
354 fe(i,j)=cff* &
355 & (hz(i,j,k)+hz(i,j-1,k))* &
356 & (0.75_r8*(t(i,j ,k,nrhs,itrc)- &
357 & t(i,j-1,k,nrhs,itrc))+ &
358 & 0.25_r8*(t(i,j ,k,nstp,itrc)- &
359 & t(i,j-1,k,nstp,itrc)))
360 tl_fe(i,j)=cff* &
361 & ((tl_hz(i,j,k)+tl_hz(i,j-1,k))* &
362 & (0.75_r8*(t(i,j ,k,nrhs,itrc)- &
363 & t(i,j-1,k,nrhs,itrc))+ &
364 & 0.25_r8*(t(i,j ,k,nstp,itrc)- &
365 & t(i,j-1,k,nstp,itrc)))+ &
366 & (hz(i,j,k)+hz(i,j-1,k))* &
367 & (0.75_r8*(tl_t(i,j ,k,nrhs,itrc)- &
368 & tl_t(i,j-1,k,nrhs,itrc))+ &
369 & 0.25_r8*(tl_t(i,j ,k,nstp,itrc)- &
370 & tl_t(i,j-1,k,nstp,itrc))))
371#elif defined TS_MIX_CLIMA
372 IF (ltracerclm(itrc,ng)) THEN
373 fe(i,j)=cff* &
374 & (hz(i,j,k)+hz(i,j-1,k))* &
375 & ((t(i,j ,k,nrhs,itrc)-tclm(i,j ,k,itrc))- &
376 & (t(i,j-1,k,nrhs,itrc)-tclm(i,j-1,k,itrc)))
377 tl_fe(i,j)=cff* &
378 & ((tl_hz(i,j,k)+tl_hz(i,j-1,k))* &
379 & ((t(i,j ,k,nrhs,itrc)- &
380 & tclm(i,j ,k,itrc))- &
381 & (t(i,j-1,k,nrhs,itrc)- &
382 & tclm(i,j-1,k,itrc)))+ &
383 & (hz(i,j,k)+hz(i,j-1,k))* &
384 & (tl_t(i,j ,k,nrhs,itrc)- &
385 & tl_t(i,j-1,k,nrhs,itrc)))
386 ELSE
387 fe(i,j)=cff* &
388 & (hz(i,j,k)+hz(i,j-1,k))* &
389 & (t(i,j ,k,nrhs,itrc)- &
390 & t(i,j-1,k,nrhs,itrc))
391 tl_fe(i,j)=cff* &
392 & ((tl_hz(i,j,k)+tl_hz(i,j-1,k))* &
393 & (t(i,j ,k,nrhs,itrc)- &
394 & t(i,j-1,k,nrhs,itrc))+ &
395 & (hz(i,j,k)+hz(i,j-1,k))* &
396 & (tl_t(i,j ,k,nrhs,itrc)- &
397 & tl_t(i,j-1,k,nrhs,itrc)))
398 END IF
399#else
400 fe(i,j)=cff* &
401 & (hz(i,j,k)+hz(i,j-1,k))* &
402 & (t(i,j ,k,nrhs,itrc)- &
403 & t(i,j-1,k,nrhs,itrc))
404 tl_fe(i,j)=cff* &
405 & ((tl_hz(i,j,k)+tl_hz(i,j-1,k))* &
406 & (t(i,j ,k,nrhs,itrc)- &
407 & t(i,j-1,k,nrhs,itrc))+ &
408 & (hz(i,j,k)+hz(i,j-1,k))* &
409 & (tl_t(i,j ,k,nrhs,itrc)- &
410 & tl_t(i,j-1,k,nrhs,itrc)))
411#endif
412 END DO
413 END DO
414!
415! Compute first harmonic operator and multiply by the metrics of the
416! second harmonic operator.
417!
418 DO j=jmin,jmax
419 DO i=imin,imax
420 cff=1.0_r8/hz(i,j,k)
421 tl_cff=-cff*cff*tl_hz(i,j,k)
422 lapt(i,j)=pm(i,j)*pn(i,j)*cff* &
423 & (fx(i+1,j)-fx(i,j)+ &
424 & fe(i,j+1)-fe(i,j))
425 tl_lapt(i,j)=pm(i,j)*pn(i,j)* &
426 & (tl_cff* &
427 & (fx(i+1,j)-fx(i,j)+ &
428 & fe(i,j+1)-fe(i,j))+ &
429 & cff* &
430 & (tl_fx(i+1,j)-tl_fx(i,j)+ &
431 & tl_fe(i,j+1)-tl_fe(i,j)))
432 END DO
433 END DO
434!
435! Apply boundary conditions (except periodic; closed or gradient)
436! to the first harmonic operator.
437!
438 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
439 IF (domain(ng)%Western_Edge(tile)) THEN
440 IF (tl_lbc(iwest,istvar(itrc),ng)%closed) THEN
441 DO j=jmin,jmax
442 lapt(istr-1,j)=0.0_r8
443 tl_lapt(istr-1,j)=0.0_r8
444 END DO
445 ELSE
446 DO j=jmin,jmax
447 lapt(istr-1,j)=lapt(istr,j)
448 tl_lapt(istr-1,j)=tl_lapt(istr,j)
449 END DO
450 END IF
451 END IF
452 END IF
453!
454 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
455 IF (domain(ng)%Eastern_Edge(tile)) THEN
456 IF (tl_lbc(ieast,istvar(itrc),ng)%closed) THEN
457 DO j=jmin,jmax
458 lapt(iend+1,j)=0.0_r8
459 tl_lapt(iend+1,j)=0.0_r8
460 END DO
461 ELSE
462 DO j=jmin,jmax
463 lapt(iend+1,j)=lapt(iend,j)
464 tl_lapt(iend+1,j)=tl_lapt(iend,j)
465 END DO
466 END IF
467 END IF
468 END IF
469!
470 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
471 IF (domain(ng)%Southern_Edge(tile)) THEN
472 IF (tl_lbc(isouth,istvar(itrc),ng)%closed) THEN
473 DO i=imin,imax
474 lapt(i,jstr-1)=0.0_r8
475 tl_lapt(i,jstr-1)=0.0_r8
476 END DO
477 ELSE
478 DO i=imin,imax
479 lapt(i,jstr-1)=lapt(i,jstr)
480 tl_lapt(i,jstr-1)=tl_lapt(i,jstr)
481 END DO
482 END IF
483 END IF
484 END IF
485!
486 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
487 IF (domain(ng)%Northern_Edge(tile)) THEN
488 IF (tl_lbc(inorth,istvar(itrc),ng)%closed) THEN
489 DO i=imin,imax
490 lapt(i,jend+1)=0.0_r8
491 tl_lapt(i,jend+1)=0.0_r8
492 END DO
493 ELSE
494 DO i=imin,imax
495 lapt(i,jend+1)=lapt(i,jend)
496 tl_lapt(i,jend+1)=tl_lapt(i,jend)
497 END DO
498 END IF
499 END IF
500 END IF
501!
502! Compute FX=d(LapT)/d(xi) and FE=d(LapT)/d(eta) terms.
503!
504 DO j=jstr,jend
505 DO i=istr,iend+1
506#ifdef DIFF_3DCOEF
507# ifdef TS_U3ADV_SPLIT_NOT_YET
508 cff=0.5_r8*diff3d_u(i,j,k)*pmon_u(i,j)
509# else
510 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
511 & pmon_u(i,j)
512# endif
513#else
514 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
515 & pmon_u(i,j)
516#endif
517!^ FX(i,j)=cff* &
518!^ & (Hz(i,j,k)+Hz(i-1,j,k))* &
519!^ & (LapT(i,j)-LapT(i-1,j))
520!^
521 tl_fx(i,j)=cff* &
522 & ((tl_hz(i,j,k)+tl_hz(i-1,j,k))* &
523 & (lapt(i,j)-lapt(i-1,j))+ &
524 & (hz(i,j,k)+hz(i-1,j,k))* &
525 & (tl_lapt(i,j)-tl_lapt(i-1,j)))
526#ifdef MASKING
527!^ FX(i,j)=FX(i,j)*umask(i,j)
528!^
529 tl_fx(i,j)=tl_fx(i,j)*umask(i,j)
530#endif
531#ifdef WET_DRY_NOT_YET
532 fx(i,j)=fx(i,j)*umask_wet(i,j)
533#endif
534 END DO
535 END DO
536 DO j=jstr,jend+1
537 DO i=istr,iend
538#ifdef DIFF_3DCOEF
539# ifdef TS_U3ADV_SPLIT_NOT_YET
540 cff=0.5_r8*diff3d_v(i,j,k)*pnom_v(i,j)
541# else
542 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
543 & pnom_v(i,j)
544# endif
545#else
546 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
547 & pnom_v(i,j)
548#endif
549!^ FE(i,j)=cff* &
550!^ & (Hz(i,j,k)+Hz(i,j-1,k))* &
551!^ & (LapT(i,j)-LapT(i,j-1))
552!^
553 tl_fe(i,j)=cff* &
554 & ((tl_hz(i,j,k)+tl_hz(i,j-1,k))* &
555 & (lapt(i,j)-lapt(i,j-1))+ &
556 & (hz(i,j,k)+hz(i,j-1,k))* &
557 & (tl_lapt(i,j)-tl_lapt(i,j-1)))
558#ifdef MASKING
559!^ FE(i,j)=FE(i,j)*vmask(i,j)
560!^
561 tl_fe(i,j)=tl_fe(i,j)*vmask(i,j)
562#endif
563#ifdef WET_DRY_NOT_YET
564 fe(i,j)=fe(i,j)*vmask_wet(i,j)
565#endif
566 END DO
567 END DO
568!
569! Time-step biharmonic, S-surfaces diffusion term (m Tunits).
570!
571 DO j=jstr,jend
572 DO i=istr,iend
573!^ cff=dt(ng)*pm(i,j)*pn(i,j)* &
574!^ & (FX(i+1,j)-FX(i,j)+ &
575!^ & FE(i,j+1)-FE(i,j))
576!^
577 tl_cff=dt(ng)*pm(i,j)*pn(i,j)* &
578 & (tl_fx(i+1,j)-tl_fx(i,j)+ &
579 & tl_fe(i,j+1)-tl_fe(i,j))
580!^ t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)-cff
581!^
582 tl_t(i,j,k,nnew,itrc)=tl_t(i,j,k,nnew,itrc)-tl_cff
583#ifdef DIAGNOSTICS_TS
584!! DiaTwrk(i,j,k,itrc,iThdif)=-cff
585#endif
586 END DO
587 END DO
588 END DO
589 END DO
590!
591 RETURN

References mod_scalars::compositegrid, mod_param::domain, mod_scalars::dt, mod_scalars::ewperiodic, mod_scalars::ieast, mod_scalars::inorth, mod_scalars::isouth, mod_ncparam::istvar, mod_scalars::iwest, mod_param::lm, mod_scalars::ltracerclm, mod_param::mm, mod_scalars::nsperiodic, and mod_param::tl_lbc.