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

Functions/Subroutines

subroutine, public rp_t3dmix4 (ng, tile)
 
subroutine rp_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 rp_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 rp_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

◆ rp_t3dmix4()

subroutine public rp_t3dmix4_mod::rp_t3dmix4 ( integer, intent(in) ng,
integer, intent(in) tile )

Definition at line 25 of file rp_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, irpm, 28, __line__, myfile)
53#endif
54 CALL rp_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, irpm, 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::irpm
integer, parameter irpm
Definition mod_param.F:664
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::irpm, mod_mixing::mixing, mod_stepping::nnew, mod_stepping::nrhs, mod_stepping::nstp, mod_ocean::ocean, rp_t3dmix4_geo_tile(), wclock_off(), and wclock_on().

Referenced by rp_rhs3d_mod::rp_rhs3d().

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

◆ rp_t3dmix4_geo_tile()

subroutine rp_t3dmix4_mod::rp_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 rp_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#ifdef TL_IOMS
397 & 2.0_r8*cff
398#endif
399#if defined TS_MIX_STABILITY
400 dtdz(i,j,k2)=cff*(0.75_r8*(t(i,j,k+1,nrhs,itrc)- &
401 & t(i,j,k ,nrhs,itrc))+ &
402 & 0.25_r8*(t(i,j,k+1,nstp,itrc)- &
403 & t(i,j,k ,nstp,itrc)))
404 tl_dtdz(i,j,k2)=tl_cff*(0.75_r8*(t(i,j,k+1,nrhs,itrc)- &
405 & t(i,j,k ,nrhs,itrc))+ &
406 & 0.25_r8*(t(i,j,k+1,nstp,itrc)- &
407 & t(i,j,k ,nstp,itrc)))+&
408 & cff*(0.75_r8*(tl_t(i,j,k+1,nrhs,itrc)- &
409 & tl_t(i,j,k ,nrhs,itrc))+ &
410 & 0.25_r8*(tl_t(i,j,k+1,nstp,itrc)- &
411 & tl_t(i,j,k ,nstp,itrc)))-&
412# ifdef TL_IOMS
413 & dtdz(i,j,k2)
414# endif
415#elif defined TS_MIX_CLIMA
416 IF (ltracerclm(itrc,ng)) THEN
417 dtdz(i,j,k2)=cff*((t(i,j,k+1,nrhs,itrc)- &
418 & tclm(i,j,k+1,itrc))- &
419 & (t(i,j,k ,nrhs,itrc)- &
420 & tclm(i,j,k ,itrc)))
421 tl_dtdz(i,j,k2)=tl_cff*((t(i,j,k+1,nrhs,itrc)- &
422 & tclm(i,j,k+1,itrc))- &
423 & (t(i,j,k ,nrhs,itrc)- &
424 & tclm(i,j,k ,itrc)))+ &
425 & cff*(tl_t(i,j,k+1,nrhs,itrc)- &
426 & tl_t(i,j,k ,nrhs,itrc))- &
427# ifdef TL_IOMS
428 & dtdz(i,j,k2)
429# endif
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# ifdef TL_IOMS
438 & dtdz(i,j,k2)
439# endif
440 END IF
441#else
442 dtdz(i,j,k2)=cff*(t(i,j,k+1,nrhs,itrc)- &
443 & t(i,j,k ,nrhs,itrc))
444 tl_dtdz(i,j,k2)=tl_cff*(t(i,j,k+1,nrhs,itrc)- &
445 & t(i,j,k ,nrhs,itrc))+ &
446 & cff*(tl_t(i,j,k+1,nrhs,itrc)- &
447 & tl_t(i,j,k ,nrhs,itrc))- &
448# ifdef TL_IOMS
449 & dtdz(i,j,k2)
450# endif
451#endif
452 END DO
453 END DO
454 END IF
455 IF (k.gt.0) THEN
456 DO j=jmin,jmax
457 DO i=imin,imax+1
458#ifdef DIFF_3DCOEF
459# ifdef TS_U3ADV_SPLIT
460 cff=0.5_r8*diff3d_u(i,j,k)*on_u(i,j)
461# else
462 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
463 & on_u(i,j)
464# endif
465#else
466 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
467 & on_u(i,j)
468#endif
469 fx(i,j)=cff* &
470 & (hz(i,j,k)+hz(i-1,j,k))* &
471 & (dtdx(i,j,k1)- &
472 & 0.5_r8*(min(dzdx(i,j,k1),0.0_r8)* &
473 & (dtdz(i-1,j,k1)+ &
474 & dtdz(i ,j,k2))+ &
475 & max(dzdx(i,j,k1),0.0_r8)* &
476 & (dtdz(i-1,j,k2)+ &
477 & dtdz(i ,j,k1))))
478 tl_fx(i,j)=cff* &
479 & ((tl_hz(i,j,k)+tl_hz(i-1,j,k))* &
480 & (dtdx(i,j,k1)- &
481 & 0.5_r8*(min(dzdx(i,j,k1),0.0_r8)* &
482 & (dtdz(i-1,j,k1)+ &
483 & dtdz(i ,j,k2))+ &
484 & max(dzdx(i,j,k1),0.0_r8)* &
485 & (dtdz(i-1,j,k2)+ &
486 & dtdz(i ,j,k1))))+ &
487 & (hz(i,j,k)+hz(i-1,j,k))* &
488 & (tl_dtdx(i,j,k1)- &
489 & 0.5_r8*(min(dzdx(i,j,k1),0.0_r8)* &
490 & (tl_dtdz(i-1,j,k1)+ &
491 & tl_dtdz(i ,j,k2))+ &
492 & max(dzdx(i,j,k1),0.0_r8)* &
493 & (tl_dtdz(i-1,j,k2)+ &
494 & tl_dtdz(i ,j,k1)))- &
495 & 0.5_r8*((0.5_r8+ &
496 & sign(0.5_r8,-dzdx(i,j,k1)))* &
497 & tl_dzdx(i,j,k1)* &
498 & (dtdz(i-1,j,k1)+dtdz(i,j,k2))+ &
499 & (0.5_r8+ &
500 & sign(0.5_r8, dzdx(i,j,k1)))* &
501 & tl_dzdx(i,j,k1)* &
502 & (dtdz(i-1,j,k2)+dtdz(i,j,k1)))))- &
503#ifdef TL_IOMS
504 & cff* &
505 & (hz(i,j,k)+hz(i-1,j,k))* &
506 & (dtdx(i,j,k1)- &
507 & (min(dzdx(i,j,k1),0.0_r8)* &
508 & (dtdz(i-1,j,k1)+ &
509 & dtdz(i ,j,k2))+ &
510 & max(dzdx(i,j,k1),0.0_r8)* &
511 & (dtdz(i-1,j,k2)+ &
512 & dtdz(i ,j,k1))))
513#endif
514 END DO
515 END DO
516 DO j=jmin,jmax+1
517 DO i=imin,imax
518#ifdef DIFF_3DCOEF
519# ifdef TS_U3ADV_SPLIT
520 cff=0.5_r8*diff3d_v(i,j,k)*om_v(i,j)
521# else
522 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
523 & om_v(i,j)
524# endif
525#else
526 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
527 & om_v(i,j)
528#endif
529 fe(i,j)=cff* &
530 & (hz(i,j,k)+hz(i,j-1,k))* &
531 & (dtde(i,j,k1)- &
532 & 0.5_r8*(min(dzde(i,j,k1),0.0_r8)* &
533 & (dtdz(i,j-1,k1)+ &
534 & dtdz(i,j ,k2))+ &
535 & max(dzde(i,j,k1),0.0_r8)* &
536 & (dtdz(i,j-1,k2)+ &
537 & dtdz(i,j ,k1))))
538 tl_fe(i,j)=cff* &
539 & ((tl_hz(i,j,k)+tl_hz(i,j-1,k))* &
540 & (dtde(i,j,k1)- &
541 & 0.5_r8*(min(dzde(i,j,k1),0.0_r8)* &
542 & (dtdz(i,j-1,k1)+ &
543 & dtdz(i,j ,k2))+ &
544 & max(dzde(i,j,k1),0.0_r8)* &
545 & (dtdz(i,j-1,k2)+ &
546 & dtdz(i,j ,k1))))+ &
547 & (hz(i,j,k)+hz(i,j-1,k))* &
548 & (tl_dtde(i,j,k1)- &
549 & 0.5_r8*(min(dzde(i,j,k1),0.0_r8)* &
550 & (tl_dtdz(i,j-1,k1)+ &
551 & tl_dtdz(i,j ,k2))+ &
552 & max(dzde(i,j,k1),0.0_r8)* &
553 & (tl_dtdz(i,j-1,k2)+ &
554 & tl_dtdz(i,j ,k1)))- &
555 & 0.5_r8*((0.5_r8+ &
556 & sign(0.5_r8,-dzde(i,j,k1)))* &
557 & tl_dzde(i,j,k1)* &
558 & (dtdz(i,j-1,k1)+dtdz(i,j,k2))+ &
559 & (0.5_r8+ &
560 & sign(0.5_r8, dzde(i,j,k1)))* &
561 & tl_dzde(i,j,k1)* &
562 & (dtdz(i,j-1,k2)+dtdz(i,j,k1)))))- &
563#ifdef TL_IOMS
564 & cff* &
565 & (hz(i,j,k)+hz(i,j-1,k))* &
566 & (dtde(i,j,k1)- &
567 & (min(dzde(i,j,k1),0.0_r8)* &
568 & (dtdz(i,j-1,k1)+ &
569 & dtdz(i,j ,k2))+ &
570 & max(dzde(i,j,k1),0.0_r8)* &
571 & (dtdz(i,j-1,k2)+ &
572 & dtdz(i,j ,k1))))
573#endif
574 END DO
575 END DO
576 IF (k.lt.n(ng)) THEN
577 DO j=jmin,jmax
578 DO i=imin,imax
579#ifdef DIFF_3DCOEF
580# ifdef TS_U3ADV_SPLIT
581 difx=0.125_r8*(diff3d_u(i,j,k )+diff3d_u(i+1,j,k )+ &
582 & diff3d_u(i,j,k+1)+diff3d_u(i+1,j,k+1))
583 dife=0.125_r8*(diff3d_v(i,j,k )+diff3d_v(i,j+1,k )+ &
584 & diff3d_v(i,j,k+1)+diff3d_v(i,j+1,k+1))
585# else
586 difx=0.5_r8*diff3d_r(i,j,k)
587 dife=difx
588# endif
589#else
590 difx=0.5_r8*diff4(i,j,itrc)
591 dife=difx
592#endif
593 cff1=min(dzdx(i ,j,k1),0.0_r8)
594 cff2=min(dzdx(i+1,j,k2),0.0_r8)
595 cff3=max(dzdx(i ,j,k2),0.0_r8)
596 cff4=max(dzdx(i+1,j,k1),0.0_r8)
597 tl_cff1=(0.5_r8+sign(0.5_r8,-dzdx(i ,j,k1)))* &
598 & tl_dzdx(i ,j,k1)
599 tl_cff2=(0.5_r8+sign(0.5_r8,-dzdx(i+1,j,k2)))* &
600 & tl_dzdx(i+1,j,k2)
601 tl_cff3=(0.5_r8+sign(0.5_r8, dzdx(i ,j,k2)))* &
602 & tl_dzdx(i ,j,k2)
603 tl_cff4=(0.5_r8+sign(0.5_r8, dzdx(i+1,j,k1)))* &
604 & tl_dzdx(i+1,j,k1)
605 fs(i,j,k2)=difx* &
606 & (cff1*(cff1*dtdz(i,j,k2)- &
607 & dtdx(i ,j,k1))+ &
608 & cff2*(cff2*dtdz(i,j,k2)- &
609 & dtdx(i+1,j,k2))+ &
610 & cff3*(cff3*dtdz(i,j,k2)- &
611 & dtdx(i ,j,k2))+ &
612 & cff4*(cff4*dtdz(i,j,k2)- &
613 & dtdx(i+1,j,k1)))
614 tl_fs(i,j,k2)=difx* &
615 & (tl_cff1*(cff1*dtdz(i,j,k2)- &
616 & dtdx(i ,j,k1))+ &
617 & tl_cff2*(cff2*dtdz(i,j,k2)- &
618 & dtdx(i+1,j,k2))+ &
619 & tl_cff3*(cff3*dtdz(i,j,k2)- &
620 & dtdx(i ,j,k2))+ &
621 & tl_cff4*(cff4*dtdz(i,j,k2)- &
622 & dtdx(i+1,j,k1))+ &
623 & cff1*(tl_cff1*dtdz(i,j,k2)+ &
624 & cff1*tl_dtdz(i,j,k2)- &
625 & tl_dtdx(i ,j,k1))+ &
626 & cff2*(tl_cff2*dtdz(i,j,k2)+ &
627 & cff2*tl_dtdz(i,j,k2)- &
628 & tl_dtdx(i+1,j,k2))+ &
629 & cff3*(tl_cff3*dtdz(i,j,k2)+ &
630 & cff3*tl_dtdz(i,j,k2)- &
631 & tl_dtdx(i ,j,k2))+ &
632 & cff4*(tl_cff4*dtdz(i,j,k2)+ &
633 & cff4*tl_dtdz(i,j,k2)- &
634 & tl_dtdx(i+1,j,k1)))- &
635#ifdef TL_IOMS
636 & difx* &
637 & (cff1*(2.0_r8*cff1*dtdz(i,j,k2)- &
638 & dtdx(i ,j,k1))+ &
639 & cff2*(2.0_r8*cff2*dtdz(i,j,k2)- &
640 & dtdx(i+1,j,k2))+ &
641 & cff3*(2.0_r8*cff3*dtdz(i,j,k2)- &
642 & dtdx(i ,j,k2))+ &
643 & cff4*(2.0_r8*cff4*dtdz(i,j,k2)- &
644 & dtdx(i+1,j,k1)))
645#endif
646!
647 cff1=min(dzde(i,j ,k1),0.0_r8)
648 cff2=min(dzde(i,j+1,k2),0.0_r8)
649 cff3=max(dzde(i,j ,k2),0.0_r8)
650 cff4=max(dzde(i,j+1,k1),0.0_r8)
651 tl_cff1=(0.5_r8+sign(0.5_r8,-dzde(i,j ,k1)))* &
652 & tl_dzde(i,j ,k1)
653 tl_cff2=(0.5_r8+sign(0.5_r8,-dzde(i,j+1,k2)))* &
654 & tl_dzde(i,j+1,k2)
655 tl_cff3=(0.5_r8+sign(0.5_r8, dzde(i,j ,k2)))* &
656 & tl_dzde(i,j ,k2)
657 tl_cff4=(0.5_r8+sign(0.5_r8, dzde(i,j+1,k1)))* &
658 & tl_dzde(i,j+1,k1)
659 fs(i,j,k2)=fs(i,j,k2)+ &
660 & dife* &
661 & (cff1*(cff1*dtdz(i,j,k2)- &
662 & dtde(i,j ,k1))+ &
663 & cff2*(cff2*dtdz(i,j,k2)- &
664 & dtde(i,j+1,k2))+ &
665 & cff3*(cff3*dtdz(i,j,k2)- &
666 & dtde(i,j ,k2))+ &
667 & cff4*(cff4*dtdz(i,j,k2)- &
668 & dtde(i,j+1,k1)))
669 tl_fs(i,j,k2)=tl_fs(i,j,k2)+ &
670 & dife* &
671 & (tl_cff1*(cff1*dtdz(i,j,k2)- &
672 & dtde(i,j ,k1))+ &
673 & tl_cff2*(cff2*dtdz(i,j,k2)- &
674 & dtde(i,j+1,k2))+ &
675 & tl_cff3*(cff3*dtdz(i,j,k2)- &
676 & dtde(i,j ,k2))+ &
677 & tl_cff4*(cff4*dtdz(i,j,k2)- &
678 & dtde(i,j+1,k1))+ &
679 & cff1*(tl_cff1*dtdz(i,j,k2)+ &
680 & cff1*tl_dtdz(i,j,k2)- &
681 & tl_dtde(i,j ,k1))+ &
682 & cff2*(tl_cff2*dtdz(i,j,k2)+ &
683 & cff2*tl_dtdz(i,j,k2)- &
684 & tl_dtde(i,j+1,k2))+ &
685 & cff3*(tl_cff3*dtdz(i,j,k2)+ &
686 & cff3*tl_dtdz(i,j,k2)- &
687 & tl_dtde(i,j ,k2))+ &
688 & cff4*(tl_cff4*dtdz(i,j,k2)+ &
689 & cff4*tl_dtdz(i,j,k2)- &
690 & tl_dtde(i,j+1,k1)))- &
691#ifdef TL_IOMS
692 & dife* &
693 & (cff1*(2.0_r8*cff1*dtdz(i,j,k2)- &
694 & dtde(i,j ,k1))+ &
695 & cff2*(2.0_r8*cff2*dtdz(i,j,k2)- &
696 & dtde(i,j+1,k2))+ &
697 & cff3*(2.0_r8*cff3*dtdz(i,j,k2)- &
698 & dtde(i,j ,k2))+ &
699 & cff4*(2.0_r8*cff4*dtdz(i,j,k2)- &
700 & dtde(i,j+1,k1)))
701#endif
702 END DO
703 END DO
704 END IF
705!
706! Compute first harmonic operator, without mixing coefficient.
707! Multiply by the metrics of the second harmonic operator. Save
708! into work array "LapT".
709!
710 DO j=jmin,jmax
711 DO i=imin,imax
712 cff=pm(i,j)*pn(i,j)
713 cff1=1.0_r8/hz(i,j,k)
714 tl_cff1=-cff1*cff1*tl_hz(i,j,k)+ &
715#ifdef TL_IOMS
716 & 2.0_r8*cff1
717#endif
718 lapt(i,j,k)=cff1*(cff* &
719 & (fx(i+1,j)-fx(i,j)+ &
720 & fe(i,j+1)-fe(i,j))+ &
721 & (fs(i,j,k2)-fs(i,j,k1)))
722 tl_lapt(i,j,k)=tl_cff1*(cff* &
723 & (fx(i+1,j)-fx(i,j)+ &
724 & fe(i,j+1)-fe(i,j))+ &
725 & (fs(i,j,k2)-fs(i,j,k1)))+ &
726 & cff1*(cff* &
727 & (tl_fx(i+1,j)-tl_fx(i,j)+ &
728 & tl_fe(i,j+1)-tl_fe(i,j))+ &
729 & (tl_fs(i,j,k2)-tl_fs(i,j,k1)))- &
730#ifdef TL_IOMS
731 & lapt(i,j,k)
732#endif
733 END DO
734 END DO
735 END IF
736 END DO k_loop1
737!
738! Apply boundary conditions (except periodic; closed or gradient)
739! to the first harmonic operator.
740!
741 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
742 IF (domain(ng)%Western_Edge(tile)) THEN
743 IF (tl_lbc(iwest,istvar(itrc),ng)%closed) THEN
744 DO k=1,n(ng)
745 DO j=jmin,jmax
746 lapt(istr-1,j,k)=0.0_r8
747 tl_lapt(istr-1,j,k)=0.0_r8
748 END DO
749 END DO
750 ELSE
751 DO k=1,n(ng)
752 DO j=jmin,jmax
753 lapt(istr-1,j,k)=lapt(istr,j,k)
754 tl_lapt(istr-1,j,k)=tl_lapt(istr,j,k)
755 END DO
756 END DO
757 END IF
758 END IF
759 END IF
760!
761 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
762 IF (domain(ng)%Eastern_Edge(tile)) THEN
763 IF (tl_lbc(ieast,istvar(itrc),ng)%closed) THEN
764 DO k=1,n(ng)
765 DO j=jmin,jmax
766 lapt(iend+1,j,k)=0.0_r8
767 tl_lapt(iend+1,j,k)=0.0_r8
768 END DO
769 END DO
770 ELSE
771 DO k=1,n(ng)
772 DO j=jmin,jmax
773 lapt(iend+1,j,k)=lapt(iend,j,k)
774 tl_lapt(iend+1,j,k)=tl_lapt(iend,j,k)
775 END DO
776 END DO
777 END IF
778 END IF
779 END IF
780!
781 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
782 IF (domain(ng)%Southern_Edge(tile)) THEN
783 IF (tl_lbc(isouth,istvar(itrc),ng)%closed) THEN
784 DO k=1,n(ng)
785 DO i=imin,imax
786 lapt(i,jstr-1,k)=0.0_r8
787 tl_lapt(i,jstr-1,k)=0.0_r8
788 END DO
789 END DO
790 ELSE
791 DO k=1,n(ng)
792 DO i=imin,imax
793 lapt(i,jstr-1,k)=lapt(i,jstr,k)
794 tl_lapt(i,jstr-1,k)=tl_lapt(i,jstr,k)
795 END DO
796 END DO
797 END IF
798 END IF
799 END IF
800!
801 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
802 IF (domain(ng)%Northern_Edge(tile)) THEN
803 IF (tl_lbc(inorth,istvar(itrc),ng)%closed) THEN
804 DO k=1,n(ng)
805 DO i=imin,imax
806 lapt(i,jend+1,k)=0.0_r8
807 tl_lapt(i,jend+1,k)=0.0_r8
808 END DO
809 END DO
810 ELSE
811 DO k=1,n(ng)
812 DO i=imin,imax
813 lapt(i,jend+1,k)=lapt(i,jend,k)
814 tl_lapt(i,jend+1,k)=tl_lapt(i,jend,k)
815 END DO
816 END DO
817 END IF
818 END IF
819 END IF
820!
821 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
822 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
823 IF (domain(ng)%SouthWest_Corner(tile)) THEN
824 DO k=1,n(ng)
825 lapt(istr-1,jstr-1,k)=0.5_r8* &
826 & (lapt(istr ,jstr-1,k)+ &
827 & lapt(istr-1,jstr ,k))
828 tl_lapt(istr-1,jstr-1,k)=0.5_r8* &
829 & (tl_lapt(istr ,jstr-1,k)+ &
830 & tl_lapt(istr-1,jstr ,k))
831 END DO
832 END IF
833 END IF
834
835 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
836 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
837 IF (domain(ng)%SouthEast_Corner(tile)) THEN
838 DO k=1,n(ng)
839 lapt(iend+1,jstr-1,k)=0.5_r8* &
840 & (lapt(iend ,jstr-1,k)+ &
841 & lapt(iend+1,jstr ,k))
842 tl_lapt(iend+1,jstr-1,k)=0.5_r8* &
843 & (tl_lapt(iend ,jstr-1,k)+ &
844 & tl_lapt(iend+1,jstr ,k))
845 END DO
846 END IF
847 END IF
848
849 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
850 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
851 IF (domain(ng)%NorthWest_Corner(tile)) THEN
852 DO k=1,n(ng)
853 lapt(istr-1,jend+1,k)=0.5_r8* &
854 & (lapt(istr ,jend+1,k)+ &
855 & lapt(istr-1,jend ,k))
856 tl_lapt(istr-1,jend+1,k)=0.5_r8* &
857 & (tl_lapt(istr ,jend+1,k)+ &
858 & tl_lapt(istr-1,jend ,k))
859 END DO
860 END IF
861 END IF
862
863 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
864 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
865 IF (domain(ng)%NorthEast_Corner(tile)) THEN
866 DO k=1,n(ng)
867 lapt(iend+1,jend+1,k)=0.5_r8* &
868 & (lapt(iend ,jend+1,k)+ &
869 & lapt(iend+1,jend ,k))
870 tl_lapt(iend+1,jend+1,k)=0.5_r8* &
871 & (tl_lapt(iend ,jend+1,k)+ &
872 & tl_lapt(iend+1,jend ,k))
873 END DO
874 END IF
875 END IF
876!
877! Compute tangent linear horizontal and vertical gradients associated
878! with the second rotated harmonic operator.
879!
880 k2=1
881 k_loop2: DO k=0,n(ng)
882 k1=k2
883 k2=3-k1
884 IF (k.lt.n(ng)) THEN
885 DO j=jstr,jend
886 DO i=istr,iend+1
887 cff=0.5_r8*(pm(i,j)+pm(i-1,j))
888#ifdef MASKING
889 cff=cff*umask(i,j)
890#endif
891#ifdef WET_DRY_NOT_YET
892 cff=cff*umask_wet(i,j)
893#endif
894 dzdx(i,j,k2)=cff*(z_r(i ,j,k+1)- &
895 & z_r(i-1,j,k+1))
896 tl_dzdx(i,j,k2)=cff*(tl_z_r(i ,j,k+1)- &
897 & tl_z_r(i-1,j,k+1))
898 dtdx(i,j,k2)=cff*(lapt(i ,j,k+1)- &
899 & lapt(i-1,j,k+1))
900 tl_dtdx(i,j,k2)=cff*(tl_lapt(i ,j,k+1)- &
901 & tl_lapt(i-1,j,k+1))
902 END DO
903 END DO
904 DO j=jstr,jend+1
905 DO i=istr,iend
906 cff=0.5_r8*(pn(i,j)+pn(i,j-1))
907#ifdef MASKING
908 cff=cff*vmask(i,j)
909#endif
910#ifdef WET_DRY_NOT_YET
911 cff=cff*vmask_wet(i,j)
912#endif
913 dzde(i,j,k2)=cff*(z_r(i,j ,k+1)- &
914 & z_r(i,j-1,k+1))
915 tl_dzde(i,j,k2)=cff*(tl_z_r(i,j ,k+1)- &
916 & tl_z_r(i,j-1,k+1))
917 dtde(i,j,k2)=cff*(lapt(i,j ,k+1)- &
918 & lapt(i,j-1,k+1))
919 tl_dtde(i,j,k2)=cff*(tl_lapt(i,j ,k+1)- &
920 & tl_lapt(i,j-1,k+1))
921 END DO
922 END DO
923 END IF
924 IF ((k.eq.0).or.(k.eq.n(ng))) THEN
925 DO j=jstr-1,jend+1
926 DO i=istr-1,iend+1
927 dtdz(i,j,k2)=0.0_r8
928 tl_dtdz(i,j,k2)=0.0_r8
929 fs(i,j,k2)=0.0_r8
930 tl_fs(i,j,k2)=0.0_r8
931 END DO
932 END DO
933 ELSE
934 DO j=jstr-1,jend+1
935 DO i=istr-1,iend+1
936 cff=1.0_r8/(z_r(i,j,k+1)- &
937 & z_r(i,j,k ))
938 tl_cff=-cff*cff*(tl_z_r(i,j,k+1)- &
939 & tl_z_r(i,j,k ))+ &
940#ifdef TL_IOMS
941 & 2.0_r8*cff
942#endif
943 dtdz(i,j,k2)=cff*(lapt(i,j,k+1)- &
944 & lapt(i,j,k ))
945 tl_dtdz(i,j,k2)=tl_cff*(lapt(i,j,k+1)- &
946 & lapt(i,j,k ))+ &
947 & cff*(tl_lapt(i,j,k+1)- &
948 & tl_lapt(i,j,k ))- &
949#ifdef TL_IOMS
950 & dtdz(i,j,k2)
951#endif
952 END DO
953 END DO
954 END IF
955!
956! Compute components of the rotated tracer flux (T m4/s) along
957! geopotential surfaces.
958!
959 IF (k.gt.0) THEN
960 DO j=jstr,jend
961 DO i=istr,iend+1
962#ifdef DIFF_3DCOEF
963# ifdef TS_U3ADV_SPLIT
964 cff=0.5_r8*diff3d_u(i,j,k)*on_u(i,j)
965# else
966 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
967 & on_u(i,j)
968# endif
969#else
970 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
971 & on_u(i,j)
972#endif
973!^ FX(i,j)=cff*
974!^ & (Hz(i,j,k)+Hz(i-1,j,k))* &
975!^ & (dTdx(i ,j,k1)- &
976!^ & 0.5_r8*(MIN(dZdx(i,j,k1),0.0_r8)* &
977!^ & (dTdz(i-1,j,k1)+ &
978!^ & dTdz(i ,j,k2))+ &
979!^ & MAX(dZdx(i,j,k1),0.0_r8)* &
980!^ & (dTdz(i-1,j,k2)+ &
981!^ & dTdz(i ,j,k1))))
982!^
983 tl_fx(i,j)=cff* &
984 & ((tl_hz(i,j,k)+tl_hz(i-1,j,k))* &
985 & (dtdx(i,j,k1)- &
986 & 0.5_r8*(min(dzdx(i,j,k1),0.0_r8)* &
987 & (dtdz(i-1,j,k1)+ &
988 & dtdz(i ,j,k2))+ &
989 & max(dzdx(i,j,k1),0.0_r8)* &
990 & (dtdz(i-1,j,k2)+ &
991 & dtdz(i ,j,k1))))+ &
992 & (hz(i,j,k)+hz(i-1,j,k))* &
993 & (tl_dtdx(i,j,k1)- &
994 & 0.5_r8*(min(dzdx(i,j,k1),0.0_r8)* &
995 & (tl_dtdz(i-1,j,k1)+ &
996 & tl_dtdz(i ,j,k2))+ &
997 & max(dzdx(i,j,k1),0.0_r8)* &
998 & (tl_dtdz(i-1,j,k2)+ &
999 & tl_dtdz(i ,j,k1)))- &
1000 & 0.5_r8*((0.5_r8+ &
1001 & sign(0.5_r8,-dzdx(i,j,k1)))* &
1002 & tl_dzdx(i,j,k1)* &
1003 & (dtdz(i-1,j,k1)+dtdz(i,j,k2))+ &
1004 & (0.5_r8+ &
1005 & sign(0.5_r8, dzdx(i,j,k1)))* &
1006 & tl_dzdx(i,j,k1)* &
1007 & (dtdz(i-1,j,k2)+dtdz(i,j,k1)))))- &
1008#ifdef TL_IOMS
1009 & cff* &
1010 & (hz(i,j,k)+hz(i-1,j,k))* &
1011 & (dtdx(i ,j,k1)- &
1012 & (min(dzdx(i,j,k1),0.0_r8)* &
1013 & (dtdz(i-1,j,k1)+ &
1014 & dtdz(i ,j,k2))+ &
1015 & max(dzdx(i,j,k1),0.0_r8)* &
1016 & (dtdz(i-1,j,k2)+ &
1017 & dtdz(i ,j,k1))))
1018#endif
1019 END DO
1020 END DO
1021 DO j=jstr,jend+1
1022 DO i=istr,iend
1023#ifdef DIFF_3DCOEF
1024# ifdef TS_U3ADV_SPLIT
1025 cff=0.5_r8*diff3d_v(i,j,k)*om_v(i,j)
1026# else
1027 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
1028 & om_v(i,j)
1029# endif
1030#else
1031 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
1032 & om_v(i,j)
1033#endif
1034!^ FE(i,j)=cff* &
1035!^ & (Hz(i,j,k)+Hz(i,j-1,k))* &
1036!^ & (dTde(i,j,k1)- &
1037!^ & 0.5_r8*(MIN(dZde(i,j,k1),0.0_r8)* &
1038!^ & (dTdz(i,j-1,k1)+ &
1039!^ & dTdz(i,j ,k2))+ &
1040!^ & MAX(dZde(i,j,k1),0.0_r8)* &
1041!^ & (dTdz(i,j-1,k2)+ &
1042!^ & dTdz(i,j ,k1))))
1043!^
1044 tl_fe(i,j)=cff* &
1045 & ((tl_hz(i,j,k)+tl_hz(i,j-1,k))* &
1046 & (dtde(i,j,k1)- &
1047 & 0.5_r8*(min(dzde(i,j,k1),0.0_r8)* &
1048 & (dtdz(i,j-1,k1)+ &
1049 & dtdz(i,j ,k2))+ &
1050 & max(dzde(i,j,k1),0.0_r8)* &
1051 & (dtdz(i,j-1,k2)+ &
1052 & dtdz(i,j ,k1))))+ &
1053 & (hz(i,j,k)+hz(i,j-1,k))* &
1054 & (tl_dtde(i,j,k1)- &
1055 & 0.5_r8*(min(dzde(i,j,k1),0.0_r8)* &
1056 & (tl_dtdz(i,j-1,k1)+ &
1057 & tl_dtdz(i,j ,k2))+ &
1058 & max(dzde(i,j,k1),0.0_r8)* &
1059 & (tl_dtdz(i,j-1,k2)+ &
1060 & tl_dtdz(i,j ,k1)))- &
1061 & 0.5_r8*((0.5_r8+ &
1062 & sign(0.5_r8,-dzde(i,j,k1)))* &
1063 & tl_dzde(i,j,k1)* &
1064 & (dtdz(i,j-1,k1)+dtdz(i,j,k2))+ &
1065 & (0.5_r8+ &
1066 & sign(0.5_r8, dzde(i,j,k1)))* &
1067 & tl_dzde(i,j,k1)* &
1068 & (dtdz(i,j-1,k2)+dtdz(i,j,k1)))))- &
1069#ifdef TL_IOMS
1070 & cff* &
1071 & (hz(i,j,k)+hz(i,j-1,k))* &
1072 & (dtde(i,j,k1)- &
1073 & (min(dzde(i,j,k1),0.0_r8)* &
1074 & (dtdz(i,j-1,k1)+ &
1075 & dtdz(i,j ,k2))+ &
1076 & max(dzde(i,j,k1),0.0_r8)* &
1077 & (dtdz(i,j-1,k2)+ &
1078 & dtdz(i,j ,k1))))
1079#endif
1080 END DO
1081 END DO
1082 IF (k.lt.n(ng)) THEN
1083 DO j=jstr,jend
1084 DO i=istr,iend
1085#ifdef DIFF_3DCOEF
1086# ifdef TS_U3ADV_SPLIT
1087 difx=0.125_r8*(diff3d_u(i,j,k )+diff3d_u(i+1,j,k )+ &
1088 & diff3d_u(i,j,k+1)+diff3d_u(i+1,j,k+1))
1089 dife=0.125_r8*(diff3d_v(i,j,k )+diff3d_v(i,j+1,k )+ &
1090 & diff3d_v(i,j,k+1)+diff3d_v(i,j+1,k+1))
1091# else
1092 difx=0.5_r8*diff3d_r(i,j,k)
1093 dife=difx
1094# endif
1095#else
1096 difx=0.5_r8*diff4(i,j,itrc)
1097 dife=difx
1098#endif
1099 cff1=min(dzdx(i ,j,k1),0.0_r8)
1100 cff2=min(dzdx(i+1,j,k2),0.0_r8)
1101 cff3=max(dzdx(i ,j,k2),0.0_r8)
1102 cff4=max(dzdx(i+1,j,k1),0.0_r8)
1103 tl_cff1=(0.5_r8+sign(0.5_r8,-dzdx(i ,j,k1)))* &
1104 & tl_dzdx(i ,j,k1)
1105 tl_cff2=(0.5_r8+sign(0.5_r8,-dzdx(i+1,j,k2)))* &
1106 & tl_dzdx(i+1,j,k2)
1107 tl_cff3=(0.5_r8+sign(0.5_r8, dzdx(i ,j,k2)))* &
1108 & tl_dzdx(i ,j,k2)
1109 tl_cff4=(0.5_r8+sign(0.5_r8, dzdx(i+1,j,k1)))* &
1110 & tl_dzdx(i+1,j,k1)
1111!^
1112!^ FS(i,j,k2)=difx* &
1113!^ & (cff1*(cff1*dTdz(i,j,k2)- &
1114!^ & dTdx(i ,j,k1))+ &
1115!^ & cff2*(cff2*dTdz(i,j,k2)- &
1116!^ & dTdx(i+1,j,k2))+ &
1117!^ & cff3*(cff3*dTdz(i,j,k2)- &
1118!^ & dTdx(i ,j,k2))+ &
1119!^ & cff4*(cff4*dTdz(i,j,k2)- &
1120!^ & dTdx(i+1,j,k1)))
1121!^
1122 tl_fs(i,j,k2)=difx* &
1123 & (tl_cff1*(cff1*dtdz(i,j,k2)- &
1124 & dtdx(i ,j,k1))+ &
1125 & tl_cff2*(cff2*dtdz(i,j,k2)- &
1126 & dtdx(i+1,j,k2))+ &
1127 & tl_cff3*(cff3*dtdz(i,j,k2)- &
1128 & dtdx(i ,j,k2))+ &
1129 & tl_cff4*(cff4*dtdz(i,j,k2)- &
1130 & dtdx(i+1,j,k1))+ &
1131 & cff1*(tl_cff1*dtdz(i,j,k2)+ &
1132 & cff1*tl_dtdz(i,j,k2)- &
1133 & tl_dtdx(i ,j,k1))+ &
1134 & cff2*(tl_cff2*dtdz(i,j,k2)+ &
1135 & cff2*tl_dtdz(i,j,k2)- &
1136 & tl_dtdx(i+1,j,k2))+ &
1137 & cff3*(tl_cff3*dtdz(i,j,k2)+ &
1138 & cff3*tl_dtdz(i,j,k2)- &
1139 & tl_dtdx(i ,j,k2))+ &
1140 & cff4*(tl_cff4*dtdz(i,j,k2)+ &
1141 & cff4*tl_dtdz(i,j,k2)- &
1142 & tl_dtdx(i+1,j,k1)))- &
1143#ifdef TL_IOMS
1144 & difx* &
1145 & (cff1*(2.0_r8*cff1*dtdz(i,j,k2)- &
1146 & dtdx(i ,j,k1))+ &
1147 & cff2*(2.0_r8*cff2*dtdz(i,j,k2)- &
1148 & dtdx(i+1,j,k2))+ &
1149 & cff3*(2.0_r8*cff3*dtdz(i,j,k2)- &
1150 & dtdx(i ,j,k2))+ &
1151 & cff4*(2.0_r8*cff4*dtdz(i,j,k2)- &
1152 & dtdx(i+1,j,k1)))
1153#endif
1154!
1155 cff1=min(dzde(i,j ,k1),0.0_r8)
1156 cff2=min(dzde(i,j+1,k2),0.0_r8)
1157 cff3=max(dzde(i,j ,k2),0.0_r8)
1158 cff4=max(dzde(i,j+1,k1),0.0_r8)
1159 tl_cff1=(0.5_r8+sign(0.5_r8,-dzde(i,j ,k1)))* &
1160 & tl_dzde(i,j ,k1)
1161 tl_cff2=(0.5_r8+sign(0.5_r8,-dzde(i,j+1,k2)))* &
1162 & tl_dzde(i,j+1,k2)
1163 tl_cff3=(0.5_r8+sign(0.5_r8, dzde(i,j ,k2)))* &
1164 & tl_dzde(i,j ,k2)
1165 tl_cff4=(0.5_r8+sign(0.5_r8, dzde(i,j+1,k1)))* &
1166 & tl_dzde(i,j+1,k1)
1167!^ FS(i,j,k2)=FS(i,j,k2)+ &
1168!^ & dife* &
1169!^ & (cff1*(cff1*dTdz(i,j,k2)- &
1170!^ & dTde(i,j ,k1))+ &
1171!^ & cff2*(cff2*dTdz(i,j,k2)- &
1172!^ & dTde(i,j+1,k2))+ &
1173!^ & cff3*(cff3*dTdz(i,j,k2)- &
1174!^ & dTde(i,j ,k2))+ &
1175!^ & cff4*(cff4*dTdz(i,j,k2)- &
1176!^ & dTde(i,j+1,k1)))
1177!^
1178 tl_fs(i,j,k2)=tl_fs(i,j,k2)+ &
1179 & dife* &
1180 & (tl_cff1*(cff1*dtdz(i,j,k2)- &
1181 & dtde(i,j ,k1))+ &
1182 & tl_cff2*(cff2*dtdz(i,j,k2)- &
1183 & dtde(i,j+1,k2))+ &
1184 & tl_cff3*(cff3*dtdz(i,j,k2)- &
1185 & dtde(i,j ,k2))+ &
1186 & tl_cff4*(cff4*dtdz(i,j,k2)- &
1187 & dtde(i,j+1,k1))+ &
1188 & cff1*(tl_cff1*dtdz(i,j,k2)+ &
1189 & cff1*tl_dtdz(i,j,k2)- &
1190 & tl_dtde(i,j ,k1))+ &
1191 & cff2*(tl_cff2*dtdz(i,j,k2)+ &
1192 & cff2*tl_dtdz(i,j,k2)- &
1193 & tl_dtde(i,j+1,k2))+ &
1194 & cff3*(tl_cff3*dtdz(i,j,k2)+ &
1195 & cff3*tl_dtdz(i,j,k2)- &
1196 & tl_dtde(i,j ,k2))+ &
1197 & cff4*(tl_cff4*dtdz(i,j,k2)+ &
1198 & cff4*tl_dtdz(i,j,k2)- &
1199 & tl_dtde(i,j+1,k1)))- &
1200#ifdef TL_IOMS
1201 & dife* &
1202 & (cff1*(2.0_r8*cff1*dtdz(i,j,k2)- &
1203 & dtde(i,j ,k1))+ &
1204 & cff2*(2.0_r8*cff2*dtdz(i,j,k2)- &
1205 & dtde(i,j+1,k2))+ &
1206 & cff3*(2.0_r8*cff3*dtdz(i,j,k2)- &
1207 & dtde(i,j ,k2))+ &
1208 & cff4*(2.0_r8*cff4*dtdz(i,j,k2)- &
1209 & dtde(i,j+1,k1)))
1210#endif
1211 END DO
1212 END DO
1213 END IF
1214!
1215! Time-step biharmonic, geopotential diffusion term (m Tunits).
1216!
1217 DO j=jstr,jend
1218 DO i=istr,iend
1219!^ cff=dt(ng)*pm(i,j)*pn(i,j)* &
1220!^ & (FX(i+1,j )-FX(i,j)+ &
1221!^ & FE(i ,j+1)-FE(i,j))+ &
1222!^ & dt(ng)*(FS(i,j,k2)-FS(i,j,k1))
1223!^
1224 tl_cff=dt(ng)*pm(i,j)*pn(i,j)* &
1225 & (tl_fx(i+1,j)-tl_fx(i,j)+ &
1226 & tl_fe(i,j+1)-tl_fe(i,j))+ &
1227 & dt(ng)*(tl_fs(i,j,k2)-tl_fs(i,j,k1))
1228!^ t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)-cff
1229!^
1230 tl_t(i,j,k,nnew,itrc)=tl_t(i,j,k,nnew,itrc)-tl_cff
1231#ifdef DIAGNOSTICS_TS
1232!! DiaTwrk(i,j,k,itrc,iThdif)=-cff
1233#endif
1234 END DO
1235 END DO
1236 END IF
1237 END DO k_loop2
1238 END DO t_loop
1239!
1240 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 rp_t3dmix4().

Here is the caller graph for this function:

◆ rp_t3dmix4_iso_tile()

subroutine rp_t3dmix4_mod::rp_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 rp_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# ifdef TL_IOMS
427 & cff2
428# endif
429 cff3=max(pden(i,j,k)-pden(i,j,k+1),small)
430 tl_cff3=(0.5_r8+sign(0.5_r8,pden(i,j,k)-pden(i,j,k+1)- &
431 & small))* &
432 & (tl_pden(i,j,k)-tl_pden(i,j,k+1))+ &
433# ifdef TL_IOMS
434 & (0.5_r8-sign(0.5_r8, &
435 & pden(i,j,k)-pden(i,j,k+1)-small))* &
436 & small
437# endif
438 cff4=max(cff2,cff3)
439 tl_cff4=(0.5_r8+sign(0.5_r8,cff2-cff3))*tl_cff2+ &
440 & (0.5_r8-sign(0.5_r8,cff2-cff3))*tl_cff3
441 cff=-1.0_r8/cff4
442 tl_cff=cff*cff*tl_cff4+ &
443# ifdef TL_IOMS
444 & 2.0_r8*cff
445# endif
446#elif defined TS_MIX_MIN_STRAT
447 cff1=max(pden(i,j,k)-pden(i,j,k+1), &
448 & strat_min*(z_r(i,j,k+1)-z_r(i,j,k)))
449 tl_cff1=(0.5_r8+sign(0.5_r8, &
450 & pden(i,j,k)-pden(i,j,k+1)- &
451 & strat_min*(z_r(i,j,k+1)- &
452 & z_r(i,j,k ))))* &
453 & (tl_pden(i,j,k)-tl_pden(i,j,k+1))+ &
454 & (0.5_r8-sign(0.5_r8, &
455 & pden(i,j,k)-pden(i,j,k+1)- &
456 & strat_min*(z_r(i,j,k+1)- &
457 & z_r(i,j,k ))))* &
458 & (strat_min*(tl_z_r(i,j,k+1)-tl_z_r(i,j,k )))
459 cff=-1.0_r8/cff1
460 tl_cff=cff*cff*tl_cff1+ &
461# ifdef TL_IOMS
462 & 2.0_r8*cff
463# endif
464#else
465 cff1=max(pden(i,j,k)-pden(i,j,k+1),eps)
466 tl_cff1=(0.5_r8+sign(0.5_r8, &
467 & pden(i,j,k)-pden(i,j,k+1)-eps))* &
468 & (tl_pden(i,j,k)-tl_pden(i,j,k+1))+ &
469# ifdef TL_IOMS
470 & (0.5_r8-sign(0.5_r8, &
471 & pden(i,j,k)-pden(i,j,k+1)-eps))*eps
472# endif
473 cff=-1.0_r8/cff1
474 tl_cff=cff*cff*tl_cff1+ &
475# ifdef TL_IOMS
476 & 2.0_r8*cff
477# endif
478#endif
479#if defined TS_MIX_STABILITY
480 dtdr(i,j,k2)=cff*(0.75_r8*(t(i,j,k+1,nrhs,itrc)- &
481 & t(i,j,k ,nrhs,itrc))+ &
482 & 0.25_r8*(t(i,j,k+1,nstp,itrc)- &
483 & t(i,j,k ,nstp,itrc)))
484 tl_dtdr(i,j,k2)=tl_cff* &
485 & (0.75_r8*(t(i,j,k+1,nrhs,itrc)- &
486 & t(i,j,k ,nrhs,itrc))+ &
487 & 0.25_r8*(t(i,j,k+1,nstp,itrc)- &
488 & t(i,j,k ,nstp,itrc)))+ &
489 & cff*
490 & (0.75_r8*(tl_t(i,j,k+1,nrhs,itrc)- &
491 & tl_t(i,j,k ,nrhs,itrc))+ &
492 & 0.25_r8*(tl_t(i,j,k+1,nstp,itrc)- &
493 & tl_t(i,j,k ,nstp,itrc)))- &
494# ifdef TL_IOMS
495 & dtdr(i,j,k2)
496# endif
497#elif defined TS_MIX_CLIMA
498 IF (ltracerclm(itrc,ng)) THEN
499 dtdr(i,j,k2)=cff*((t(i,j,k+1,nrhs,itrc)- &
500 & tclm(i,j,k+1,itrc))- &
501 & (t(i,j,k ,nrhs,itrc)- &
502 & tclm(i,j,k ,itrc)))
503 tl_dtdr(i,j,k2)=tl_cff*((t(i,j,k+1,nrhs,itrc)- &
504 & tclm(i,j,k+1,itrc))- &
505 & (t(i,j,k ,nrhs,itrc)- &
506 & tclm(i,j,k ,itrc)))+ &
507 & cff*(tl_t(i,j,k+1,nrhs,itrc)- &
508 & tl_t(i,j,k ,nrhs,itrc))- &
509# ifdef TL_IOMS
510 & dtdr(i,j,k2)
511# endif
512 ELSE
513 dtdr(i,j,k2)=cff*(t(i,j,k+1,nrhs,itrc)- &
514 & t(i,j,k ,nrhs,itrc))
515 tl_dtdr(i,j,k2)=tl_cff*(t(i,j,k+1,nrhs,itrc)- &
516 & t(i,j,k ,nrhs,itrc))+ &
517 & cff*(tl_t(i,j,k+1,nrhs,itrc)- &
518 & tl_t(i,j,k ,nrhs,itrc))- &
519# ifdef TL_IOMS
520 & dtdr(i,j,k2)
521# endif
522 END IF
523#else
524 dtdr(i,j,k2)=cff*(t(i,j,k+1,nrhs,itrc)- &
525 & t(i,j,k ,nrhs,itrc))
526 tl_dtdr(i,j,k2)=tl_cff*(t(i,j,k+1,nrhs,itrc)- &
527 & t(i,j,k ,nrhs,itrc))+ &
528 & cff*(tl_t(i,j,k+1,nrhs,itrc)- &
529 & tl_t(i,j,k ,nrhs,itrc))- &
530# ifdef TL_IOMS
531 & dtdr(i,j,k2)
532# endif
533#endif
534 fs(i,j,k2)=cff*(z_r(i,j,k+1)- &
535 & z_r(i,j,k ))
536 tl_fs(i,j,k2)=tl_cff*(z_r(i,j,k+1)- &
537 & z_r(i,j,k ))+ &
538 & cff*(tl_z_r(i,j,k+1)- &
539 & tl_z_r(i,j,k ))- &
540#ifdef TL_IOMS
541 & fs(i,j,k2)
542#endif
543 END DO
544 END DO
545 END IF
546 IF (k.gt.0) THEN
547 DO j=jmin,jmax
548 DO i=imin,imax+1
549#ifdef DIFF_3DCOEF
550# ifdef TS_U3ADV_SPLIT
551 cff=0.5_r8*diff3d_u(i,j,k)*on_u(i,j)
552# else
553 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
554 & on_u(i,j)
555# endif
556#else
557 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
558 & on_u(i,j)
559#endif
560 fx(i,j)=cff* &
561 & (hz(i,j,k)+hz(i-1,j,k))* &
562 & (dtdx(i,j,k1)- &
563 & 0.5_r8*(max(drdx(i,j,k1),0.0_r8)* &
564 & (dtdr(i-1,j,k1)+ &
565 & dtdr(i ,j,k2))+ &
566 & min(drdx(i,j,k1),0.0_r8)* &
567 & (dtdr(i-1,j,k2)+ &
568 & dtdr(i ,j,k1))))
569 tl_fx(i,j)=cff* &
570 & ((tl_hz(i,j,k)+tl_hz(i-1,j,k))* &
571 & (dtdx(i,j,k1)- &
572 & 0.5_r8*(max(drdx(i,j,k1),0.0_r8)* &
573 & (dtdr(i-1,j,k1)+ &
574 & dtdr(i ,j,k2))+ &
575 & min(drdx(i,j,k1),0.0_r8)* &
576 & (dtdr(i-1,j,k2)+ &
577 & dtdr(i ,j,k1))))+ &
578 & (hz(i,j,k)+hz(i-1,j,k))* &
579 & (tl_dtdx(i,j,k1)- &
580 & 0.5_r8*(max(drdx(i,j,k1),0.0_r8)* &
581 & (tl_dtdr(i-1,j,k1)+ &
582 & tl_dtdr(i ,j,k2))+ &
583 & min(drdx(i,j,k1),0.0_r8)* &
584 & (tl_dtdr(i-1,j,k2)+ &
585 & tl_dtdr(i ,j,k1)))- &
586 & 0.5_r8*((0.5_r8+ &
587 & sign(0.5_r8, drdx(i,j,k1)))* &
588 & tl_drdx(i,j,k1)* &
589 & (dtdr(i-1,j,k1)+dtdr(i,j,k2))+ &
590 & (0.5_r8+ &
591 & sign(0.5_r8,-drdx(i,j,k1)))* &
592 & tl_drdx(i,j,k1)* &
593 & (dtdr(i-1,j,k2)+dtdr(i,j,k1)))))- &
594#ifdef TL_IOMS
595 & cff* &
596 & (hz(i,j,k)+hz(i-1,j,k))* &
597 & (dtdx(i,j,k1)- &
598 & (max(drdx(i,j,k1),0.0_r8)* &
599 & (dtdr(i-1,j,k1)+ &
600 & dtdr(i ,j,k2))+ &
601 & min(drdx(i,j,k1),0.0_r8)* &
602 & (dtdr(i-1,j,k2)+ &
603 & dtdr(i ,j,k1))))
604#endif
605 END DO
606 END DO
607 DO j=jmin,jmax+1
608 DO i=imin,imax
609#ifdef DIFF_3DCOEF
610# ifdef TS_U3ADV_SPLIT
611 cff=0.5_r8*diff3d_v(i,j,k)*om_v(i,j)
612# else
613 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
614 & om_v(i,j)
615# endif
616#else
617 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
618 & om_v(i,j)
619#endif
620 fe(i,j)=cff* &
621 & (hz(i,j,k)+hz(i,j-1,k))* &
622 & (dtde(i,j,k1)- &
623 & 0.5_r8*(max(drde(i,j,k1),0.0_r8)* &
624 & (dtdr(i,j-1,k1)+ &
625 & dtdr(i,j ,k2))+ &
626 & min(drde(i,j,k1),0.0_r8)* &
627 & (dtdr(i,j-1,k2)+ &
628 & dtdr(i,j ,k1))))
629 tl_fe(i,j)=cff* &
630 & ((tl_hz(i,j,k)+tl_hz(i,j-1,k))* &
631 & (dtde(i,j,k1)- &
632 & 0.5_r8*(max(drde(i,j,k1),0.0_r8)* &
633 & (dtdr(i,j-1,k1)+ &
634 & dtdr(i,j ,k2))+ &
635 & min(drde(i,j,k1),0.0_r8)* &
636 & (dtdr(i,j-1,k2)+ &
637 & dtdr(i,j ,k1))))+ &
638 & (hz(i,j,k)+hz(i,j-1,k))* &
639 & (tl_dtde(i,j,k1)- &
640 & 0.5_r8*(max(drde(i,j,k1),0.0_r8)* &
641 & (tl_dtdr(i,j-1,k1)+ &
642 & tl_dtdr(i,j ,k2))+ &
643 & min(drde(i,j,k1),0.0_r8)* &
644 & (tl_dtdr(i,j-1,k2)+ &
645 & tl_dtdr(i,j ,k1)))- &
646 & 0.5_r8*((0.5_r8+ &
647 & sign(0.5_r8, drde(i,j,k1)))* &
648 & tl_drde(i,j,k1)* &
649 & (dtdr(i,j-1,k1)+dtdr(i,j,k2))+ &
650 & (0.5_r8+ &
651 & sign(0.5_r8,-drde(i,j,k1)))* &
652 & tl_drde(i,j,k1)* &
653 & (dtdr(i,j-1,k2)+dtdr(i,j,k1)))))- &
654#ifdef TL_IOMS
655 & cff* &
656 & (hz(i,j,k)+hz(i,j-1,k))* &
657 & (dtde(i,j,k1)- &
658 & (max(drde(i,j,k1),0.0_r8)* &
659 & (dtdr(i,j-1,k1)+ &
660 & dtdr(i,j ,k2))+ &
661 & min(drde(i,j,k1),0.0_r8)* &
662 & (dtdr(i,j-1,k2)+ &
663 & dtdr(i,j ,k1))))
664#endif
665 END DO
666 END DO
667 IF (k.lt.n(ng)) THEN
668 DO j=jmin,jmax
669 DO i=imin,imax
670#ifdef DIFF_3DCOEF
671# ifdef TS_U3ADV_SPLIT
672 difx=0.125_r8*(diff3d_u(i,j,k )+diff3d_u(i+1,j,k )+ &
673 & diff3d_u(i,j,k+1)+diff3d_u(i+1,j,k+1))
674 dife=0.125_r8*(diff3d_v(i,j,k )+diff3d_v(i,j+1,k )+ &
675 & diff3d_v(i,j,k+1)+diff3d_v(i,j+1,k+1))
676# else
677 difx=0.5_r8*diff3d_r(i,j,k)
678 dife=difx
679# endif
680#else
681 difx=0.5_r8*diff4(i,j,itrc)
682 dife=difx
683#endif
684 cff1=max(drdx(i ,j,k1),0.0_r8)
685 cff2=max(drdx(i+1,j,k2),0.0_r8)
686 cff3=min(drdx(i ,j,k2),0.0_r8)
687 cff4=min(drdx(i+1,j,k1),0.0_r8)
688 tl_cff1=(0.5_r8+sign(0.5_r8, drdx(i ,j,k1)))* &
689 & tl_drdx(i ,j,k1)
690 tl_cff2=(0.5_r8+sign(0.5_r8, drdx(i+1,j,k2)))* &
691 & tl_drdx(i+1,j,k2)
692 tl_cff3=(0.5_r8+sign(0.5_r8,-drdx(i ,j,k2)))* &
693 & tl_drdx(i ,j,k2)
694 tl_cff4=(0.5_r8+sign(0.5_r8,-drdx(i+1,j,k1)))* &
695 & tl_drdx(i+1,j,k1)
696 cff=difx* &
697 & (cff1*(cff1*dtdr(i,j,k2)-dtdx(i ,j,k1))+ &
698 & cff2*(cff2*dtdr(i,j,k2)-dtdx(i+1,j,k2))+ &
699 & cff3*(cff3*dtdr(i,j,k2)-dtdx(i ,j,k2))+ &
700 & cff4*(cff4*dtdr(i,j,k2)-dtdx(i+1,j,k1)))
701 tl_cff=difx* &
702 & (tl_cff1*(cff1*dtdr(i ,j,k2)- &
703 & dtdx(i ,j,k1))+ &
704 & tl_cff2*(cff2*dtdr(i,j,k2)- &
705 & dtdx(i+1,j,k2))+ &
706 & tl_cff3*(cff3*dtdr(i,j,k2)- &
707 & dtdx(i ,j,k2))+ &
708 & tl_cff4*(cff4*dtdr(i,j,k2)- &
709 & dtdx(i+1,j,k1))+ &
710 & cff1*(tl_cff1*dtdr(i,j,k2)+ &
711 & cff1*tl_dtdr(i,j,k2)- &
712 & tl_dtdx(i ,j,k1))+ &
713 & cff2*(tl_cff2*dtdr(i,j,k2)+ &
714 & cff2*tl_dtdr(i,j,k2)- &
715 & tl_dtdx(i+1,j,k2))+ &
716 & cff3*(tl_cff3*dtdr(i,j,k2)+ &
717 & cff3*tl_dtdr(i,j,k2)- &
718 & tl_dtdx(i ,j,k2))+ &
719 & cff4*(tl_cff4*dtdr(i,j,k2)+ &
720 & cff4*tl_dtdr(i,j,k2)- &
721 & tl_dtdx(i+1,j,k1)))- &
722#ifdef TL_IOMS
723 & difx* &
724 & (cff1*(2.0_r8*cff1*dtdr(i,j,k2)- &
725 & dtdx(i ,j,k1))- &
726 & cff2*(2.0_r8*cff2*dtdr(i,j,k2)- &
727 & dtdx(i+1,j,k2))- &
728 & cff3*(2.0_r8*cff3*dtdr(i,j,k2)- &
729 & dtdx(i ,j,k2))- &
730 & cff4*(2.0_r8*cff4*dtdr(i,j,k2)- &
731 & dtdx(i+1,j,k1)))
732#endif
733!
734 cff1=max(drde(i,j ,k1),0.0_r8)
735 cff2=max(drde(i,j+1,k2),0.0_r8)
736 cff3=min(drde(i,j ,k2),0.0_r8)
737 cff4=min(drde(i,j+1,k1),0.0_r8)
738 tl_cff1=(0.5_r8+sign(0.5_r8, drde(i,j ,k1)))* &
739 & tl_drde(i,j ,k1)
740 tl_cff2=(0.5_r8+sign(0.5_r8, drde(i,j+1,k2)))* &
741 & tl_drde(i,j+1,k2)
742 tl_cff3=(0.5_r8+sign(0.5_r8,-drde(i,j ,k2)))* &
743 & tl_drde(i,j ,k2)
744 tl_cff4=(0.5_r8+sign(0.5_r8,-drde(i,j+1,k1)))* &
745 & tl_drde(i,j+1,k1)
746 cff=cff+ &
747 & dife* &
748 & (cff1*(cff1*dtdr(i,j,k2)-dtde(i,j ,k1))+ &
749 & cff2*(cff2*dtdr(i,j,k2)-dtde(i,j+1,k2))+ &
750 & cff3*(cff3*dtdr(i,j,k2)-dtde(i,j ,k2))+ &
751 & cff4*(cff4*dtdr(i,j,k2)-dtde(i,j+1,k1)))
752 tl_cff=tl_cff+ &
753 & dife* &
754 & (tl_cff1*(cff1*dtdr(i,j,k2)- &
755 & dtde(i,j ,k1))+ &
756 & tl_cff2*(cff2*dtdr(i,j,k2)- &
757 & dtde(i,j+1,k2))+ &
758 & tl_cff3*(cff3*dtdr(i,j,k2)- &
759 & dtde(i,j ,k2))+ &
760 & tl_cff4*(cff4*dtdr(i,j,k2)- &
761 & dtde(i,j+1,k1))+ &
762 & cff1*(tl_cff1*dtdr(i,j,k2)+ &
763 & cff1*tl_dtdr(i,j,k2)- &
764 & tl_dtde(i,j ,k1))+ &
765 & cff2*(tl_cff2*dtdr(i,j,k2)+ &
766 & cff2*tl_dtdr(i,j,k2)- &
767 & tl_dtde(i,j+1,k2))+ &
768 & cff3*(tl_cff3*dtdr(i,j,k2)+ &
769 & cff3*tl_dtdr(i,j,k2)- &
770 & tl_dtde(i,j ,k2))+ &
771 & cff4*(tl_cff4*dtdr(i,j,k2)+ &
772 & cff4*tl_dtdr(i,j,k2)- &
773 & tl_dtde(i,j+1,k1)))- &
774#ifdef TL_IOMS
775 & dife* &
776 & (cff1*(2.0_r8*cff1*dtdr(i,j,k2)- &
777 & dtde(i,j ,k1))- &
778 & cff2*(2.0_r8*cff2*dtdr(i,j,k2)- &
779 & dtde(i,j+1,k2))- &
780 & cff3*(2.0_r8*cff3*dtdr(i,j,k2)- &
781 & dtde(i,j ,k2))- &
782 & cff4*(2.0_r8*cff4*dtdr(i,j,k2)-
783 & dtde(i,j+1,k1)))
784#endif
785!^ FS(i,j,k2)=cff*FS(i,j,k2) ! recursive
786!^ ! compute after TLM
787!^
788 tl_fs(i,j,k2)=tl_cff*fs(i,j,k2)+ &
789 & cff*tl_fs(i,j,k2)- &
790#ifdef TL_IOMS
791 & cff*fs(i,j,k2)
792#endif
793 fs(i,j,k2)=cff*fs(i,j,k2) ! recursive
794 END DO
795 END DO
796 END IF
797!
798! Compute first harmonic operator, without mixing coefficient.
799! Multiply by the metrics of the second harmonic operator. Save
800! into work array "LapT".
801!
802 DO j=jmin,jmax
803 DO i=imin,imax
804 cff=pm(i,j)*pn(i,j)
805 cff1=1.0_r8/hz(i,j,k)
806 tl_cff1=-cff1*cff1*tl_hz(i,j,k)+ &
807#ifdef TL_IOMS
808 & 2.0_r8*cff1
809#endif
810 lapt(i,j,k)=cff1*(cff* &
811 & (fx(i+1,j)-fx(i,j)+ &
812 & fe(i,j+1)-fe(i,j))+ &
813 & (fs(i,j,k2)-fs(i,j,k1)))
814 tl_lapt(i,j,k)=tl_cff1*(cff* &
815 & (fx(i+1,j)-fx(i,j)+ &
816 & fe(i,j+1)-fe(i,j))+ &
817 & (fs(i,j,k2)-fs(i,j,k1)))+ &
818 & cff1*(cff* &
819 & (tl_fx(i+1,j)-tl_fx(i,j)+ &
820 & tl_fe(i,j+1)-tl_fe(i,j))+ &
821 & (tl_fs(i,j,k2)-tl_fs(i,j,k1)))- &
822#ifdef TL_IOMS
823 & lapt(i,j,k)
824#endif
825 END DO
826 END DO
827 END IF
828 END DO k_loop1
829!
830! Apply boundary conditions (except periodic; closed or gradient)
831! to the first harmonic operator.
832!
833 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
834 IF (domain(ng)%Western_Edge(tile)) THEN
835 IF (tl_lbc(iwest,istvar(itrc),ng)%closed) THEN
836 DO k=1,n(ng)
837 DO j=jmin,jmax
838 lapt(istr-1,j,k)=0.0_r8
839 tl_lapt(istr-1,j,k)=0.0_r8
840 END DO
841 END DO
842 ELSE
843 DO k=1,n(ng)
844 DO j=jmin,jmax
845 lapt(istr-1,j,k)=lapt(istr,j,k)
846 tl_lapt(istr-1,j,k)=tl_lapt(istr,j,k)
847 END DO
848 END DO
849 END IF
850 END IF
851 END IF
852!
853 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
854 IF (domain(ng)%Eastern_Edge(tile)) THEN
855 IF (tl_lbc(ieast,istvar(itrc),ng)%closed) THEN
856 DO k=1,n(ng)
857 DO j=jmin,jmax
858 lapt(iend+1,j,k)=0.0_r8
859 tl_lapt(iend+1,j,k)=0.0_r8
860 END DO
861 END DO
862 ELSE
863 DO k=1,n(ng)
864 DO j=jmin,jmax
865 lapt(iend+1,j,k)=lapt(iend,j,k)
866 tl_lapt(iend+1,j,k)=tl_lapt(iend,j,k)
867 END DO
868 END DO
869 END IF
870 END IF
871 END IF
872!
873 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
874 IF (domain(ng)%Southern_Edge(tile)) THEN
875 IF (tl_lbc(isouth,istvar(itrc),ng)%closed) THEN
876 DO k=1,n(ng)
877 DO i=imin,imax
878 lapt(i,jstr-1,k)=0.0_r8
879 tl_lapt(i,jstr-1,k)=0.0_r8
880 END DO
881 END DO
882 ELSE
883 DO k=1,n(ng)
884 DO i=imin,imax
885 lapt(i,jstr-1,k)=lapt(i,jstr,k)
886 tl_lapt(i,jstr-1,k)=tl_lapt(i,jstr,k)
887 END DO
888 END DO
889 END IF
890 END IF
891 END IF
892!
893 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
894 IF (domain(ng)%Northern_Edge(tile)) THEN
895 IF (tl_lbc(inorth,istvar(itrc),ng)%closed) THEN
896 DO k=1,n(ng)
897 DO i=imin,imax
898 lapt(i,jend+1,k)=0.0_r8
899 tl_lapt(i,jend+1,k)=0.0_r8
900 END DO
901 END DO
902 ELSE
903 DO k=1,n(ng)
904 DO i=imin,imax
905 lapt(i,jend+1,k)=lapt(i,jend,k)
906 tl_lapt(i,jend+1,k)=tl_lapt(i,jend,k)
907 END DO
908 END DO
909 END IF
910 END IF
911 END IF
912!
913 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
914 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
915 IF (domain(ng)%SouthWest_Corner(tile)) THEN
916 DO k=1,n(ng)
917 lapt(istr-1,jstr-1,k)=0.5_r8* &
918 & (lapt(istr ,jstr-1,k)+ &
919 & lapt(istr-1,jstr ,k))
920 tl_lapt(istr-1,jstr-1,k)=0.5_r8* &
921 & (tl_lapt(istr ,jstr-1,k)+ &
922 tl_lapt(istr-1,jstr ,k))
923 END DO
924 END IF
925 END IF
926
927 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
928 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
929 IF (domain(ng)%SouthEast_Corner(tile)) THEN
930 DO k=1,n(ng)
931 lapt(iend+1,jstr-1,k)=0.5_r8* &
932 & (lapt(iend ,jstr-1,k)+ &
933 & lapt(iend+1,jstr ,k))
934 tl_lapt(iend+1,jstr-1,k)=0.5_r8* &
935 & (tl_lapt(iend ,jstr-1,k)+ &
936 & tl_lapt(iend+1,jstr ,k))
937 END DO
938 END IF
939 END IF
940
941 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
942 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
943 IF (domain(ng)%NorthWest_Corner(tile)) THEN
944 DO k=1,n(ng)
945 lapt(istr-1,jend+1,k)=0.5_r8* &
946 & (lapt(istr ,jend+1,k)+ &
947 & lapt(istr-1,jend ,k))
948 tl_lapt(istr-1,jend+1,k)=0.5_r8* &
949 & (tl_lapt(istr ,jend+1,k)+ &
950 & tl_lapt(istr-1,jend ,k))
951 END DO
952 END IF
953 END IF
954
955 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
956 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
957 IF (domain(ng)%NorthEast_Corner(tile)) THEN
958 DO k=1,n(ng)
959 lapt(iend+1,jend+1,k)=0.5_r8* &
960 & (lapt(iend ,jend+1,k)+ &
961 & lapt(iend+1,jend ,k))
962 tl_lapt(iend+1,jend+1,k)=0.5_r8* &
963 & (tl_lapt(iend ,jend+1,k)+ &
964 & tl_lapt(iend+1,jend ,k))
965 END DO
966 END IF
967 END IF
968!
969! Compute horizontal and density gradients associated with the
970! second rotated harmonic operator.
971!
972 k2=1
973 k_loop2: DO k=0,n(ng)
974 k1=k2
975 k2=3-k1
976 IF (k.lt.n(ng)) THEN
977 DO j=jstr,jend
978 DO i=istr,iend+1
979 cff=0.5_r8*(pm(i,j)+pm(i-1,j))
980#ifdef MASKING
981 cff=cff*umask(i,j)
982#endif
983#ifdef WET_DRY_NOT_YET
984 cff=cff*umask_wet(i,j)
985#endif
986 drdx(i,j,k2)=cff*(pden(i ,j,k+1)- &
987 & pden(i-1,j,k+1))
988 tl_drdx(i,j,k2)=cff*(tl_pden(i ,j,k+1)- &
989 & tl_pden(i-1,j,k+1))
990 dtdx(i,j,k2)=cff*(lapt(i ,j,k+1)- &
991 & lapt(i-1,j,k+1))
992 tl_dtdx(i,j,k2)=cff*(tl_lapt(i ,j,k+1)- &
993 & tl_lapt(i-1,j,k+1))
994 END DO
995 END DO
996 DO j=jstr,jend+1
997 DO i=istr,iend
998 cff=0.5_r8*(pn(i,j)+pn(i,j-1))
999#ifdef MASKING
1000 cff=cff*vmask(i,j)
1001#endif
1002#ifdef WET_DRY_NOT_YET
1003 cff=cff*vmask_wet(i,j)
1004#endif
1005 drde(i,j,k2)=cff*(pden(i,j ,k+1)- &
1006 & pden(i,j-1,k+1))
1007 tl_drde(i,j,k2)=cff*(tl_pden(i,j ,k+1)- &
1008 & tl_pden(i,j-1,k+1))
1009 dtde(i,j,k2)=cff*(lapt(i,j ,k+1)- &
1010 & lapt(i,j-1,k+1))
1011 tl_dtde(i,j,k2)=cff*(tl_lapt(i,j ,k+1)- &
1012 & tl_lapt(i,j-1,k+1))
1013 END DO
1014 END DO
1015 END IF
1016 IF ((k.eq.0).or.(k.eq.n(ng))) THEN
1017 DO j=jstr-1,jend+1
1018 DO i=istr-1,iend+1
1019 dtdr(i,j,k2)=0.0_r8
1020 tl_dtdr(i,j,k2)=0.0_r8
1021 fs(i,j,k2)=0.0_r8
1022 tl_fs(i,j,k2)=0.0_r8
1023 END DO
1024 END DO
1025 ELSE
1026 DO j=jstr-1,jend+1
1027 DO i=istr-1,iend+1
1028#if defined TS_MIX_MAX_SLOPE
1029 cff1=sqrt(drdx(i,j,k2)**2+drdx(i+1,j,k2)**2+ &
1030 & drdx(i,j,k1)**2+drdx(i+1,j,k1)**2+ &
1031 & drde(i,j,k2)**2+drde(i,j+1,k2)**2+ &
1032 & drde(i,j,k1)**2+drde(i,j+1,k1)**2)
1033 IF (cff1.ne.0.0_r8) THEN
1034 tl_cff1=(drdx(i ,j,k2)*tl_drdx(i ,j,k2)+ &
1035 & drdx(i+1,j,k2)*tl_drdx(i+1,j,k2)+ &
1036 & drdx(i ,j,k1)*tl_drdx(i ,j,k1)+ &
1037 & drdx(i+1,j,k1)*tl_drdx(i+1,j,k1)+ &
1038 & drde(i,j ,k2)*tl_drde(i,j ,k2)+ &
1039 & drde(i,j+1,k2)*tl_drde(i,j+1,k2)+ &
1040 & drde(i,j ,k1)*tl_drde(i,j ,k1)+ &
1041 & drde(i,j+1,k1)*tl_drde(i,j+1,k1))/cff1
1042 ELSE
1043 tl_cff1=0.0_r8
1044 END IF
1045 cff2=0.25_r8*slope_max* &
1046 & (z_r(i,j,k+1)-z_r(i,j,k))*cff1
1047 tl_cff2=0.25_r8*slope_max* &
1048 & ((tl_z_r(i,j,k+1)-tl_z_r(i,j,k))*cff1+ &
1049 & (z_r(i,j,k+1)-z_r(i,j,k))*tl_cff1)- &
1050# ifdef TL_IOMS
1051 & cff2
1052# endif
1053 cff3=max(pden(i,j,k)-pden(i,j,k+1),small)
1054 tl_cff3=(0.5_r8+sign(0.5_r8,pden(i,j,k)-pden(i,j,k+1)- &
1055 & small))* &
1056 & (tl_pden(i,j,k)-tl_pden(i,j,k+1))+ &
1057# ifdef TL_IOMS
1058 & (0.5_r8-sign(0.5_r8, &
1059 & pden(i,j,k)-pden(i,j,k+1)-small))* &
1060 & small
1061# endif
1062 cff4=max(cff2,cff3)
1063 tl_cff4=(0.5_r8+sign(0.5_r8,cff2-cff3))*tl_cff2+ &
1064 & (0.5_r8-sign(0.5_r8,cff2-cff3))*tl_cff3
1065 cff=-1.0_r8/cff4
1066 tl_cff=cff*cff*tl_cff4+ &
1067# ifdef TL_IOMS
1068 & 2.0_r8*cff
1069# endif
1070#elif defined TS_MIX_MIN_STRAT
1071 cff1=max(pden(i,j,k)-pden(i,j,k+1), &
1072 & strat_min*(z_r(i,j,k+1)-z_r(i,j,k)))
1073 tl_cff1=(0.5_r8+sign(0.5_r8, &
1074 & pden(i,j,k)-pden(i,j,k+1)- &
1075 & strat_min*(z_r(i,j,k+1)- &
1076 & z_r(i,j,k ))))* &
1077 & (tl_pden(i,j,k)-tl_pden(i,j,k+1))+ &
1078 & (0.5_r8-sign(0.5_r8, &
1079 & pden(i,j,k)-pden(i,j,k+1)- &
1080 & strat_min*(z_r(i,j,k+1)- &
1081 & z_r(i,j,k ))))* &
1082 & (strat_min*(tl_z_r(i,j,k+1)-tl_z_r(i,j,k )))
1083 cff=-1.0_r8/cff1
1084 tl_cff=cff*cff*tl_cff1+ &
1085# ifdef TL_IOMS
1086 & 2.0_r8*cff
1087# endif
1088#else
1089 cff1=max(pden(i,j,k)-pden(i,j,k+1),eps)
1090 tl_cff1=(0.5_r8+sign(0.5_r8, &
1091 & pden(i,j,k)-pden(i,j,k+1)-eps))* &
1092 & (tl_pden(i,j,k)-tl_pden(i,j,k+1))+ &
1093# ifdef TL_IOMS
1094 & (0.5_r8-sign(0.5_r8, &
1095 & pden(i,j,k)-pden(i,j,k+1)-eps))*eps
1096# endif
1097 cff=-1.0_r8/cff1
1098 tl_cff=cff*cff*tl_cff1+ &
1099# ifdef TL_IOMS
1100 & 2.0_r8*cff
1101# endif
1102#endif
1103 dtdr(i,j,k2)=cff*(lapt(i,j,k+1)- &
1104 & lapt(i,j,k ))
1105 tl_dtdr(i,j,k2)=tl_cff*(lapt(i,j,k+1)- &
1106 & lapt(i,j,k ))+ &
1107 & cff*(tl_lapt(i,j,k+1)- &
1108 & tl_lapt(i,j,k ))- &
1109#ifdef TL_IOMS
1110 & dtdr(i,j,k2)
1111#endif
1112 fs(i,j,k2)=cff*(z_r(i,j,k+1)- &
1113 & z_r(i,j,k ))
1114 tl_fs(i,j,k2)=tl_cff*(z_r(i,j,k+1)- &
1115 & z_r(i,j,k ))+ &
1116 & cff*(tl_z_r(i,j,k+1)- &
1117 & tl_z_r(i,j,k ))- &
1118#ifdef TL_IOMS
1119 & fs(i,j,k2)
1120#endif
1121 END DO
1122 END DO
1123 END IF
1124!
1125! Compute components of the rotated tracer flux (T m4/s) along
1126! isopycnic surfaces.
1127!
1128 IF (k.gt.0) THEN
1129 DO j=jstr,jend
1130 DO i=istr,iend+1
1131#ifdef DIFF_3DCOEF
1132# ifdef TS_U3ADV_SPLIT
1133 cff=0.5_r8*diff3d_u(i,j,k)*on_u(i,j)
1134# else
1135 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
1136 & on_u(i,j)
1137# endif
1138#else
1139 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
1140 & on_u(i,j)
1141#endif
1142!^ FX(i,j)=cff* &
1143!^ & (Hz(i,j,k)+Hz(i-1,j,k))* &
1144!^ & (dTdx(i,j,k1)- &
1145!^ & 0.5_r8*(MAX(dRdx(i,j,k1),0.0_r8)* &
1146!^ & (dTdr(i-1,j,k1)+ &
1147!^ & dTdr(i ,j,k2))+ &
1148!^ & MIN(dRdx(i,j,k1),0.0_r8)* &
1149!^ & (dTdr(i-1,j,k2)+ &
1150!^ & dTdr(i ,j,k1))))
1151!^
1152 tl_fx(i,j)=cff* &
1153 & ((tl_hz(i,j,k)+tl_hz(i-1,j,k))* &
1154 & (dtdx(i,j,k1)- &
1155 & 0.5_r8*(max(drdx(i,j,k1),0.0_r8)* &
1156 & (dtdr(i-1,j,k1)+ &
1157 & dtdr(i ,j,k2))+ &
1158 & min(drdx(i,j,k1),0.0_r8)* &
1159 & (dtdr(i-1,j,k2)+ &
1160 & dtdr(i ,j,k1))))+ &
1161 & (hz(i,j,k)+hz(i-1,j,k))* &
1162 & (tl_dtdx(i,j,k1)- &
1163 & 0.5_r8*(max(drdx(i,j,k1),0.0_r8)* &
1164 & (tl_dtdr(i-1,j,k1)+ &
1165 & tl_dtdr(i ,j,k2))+ &
1166 & min(drdx(i,j,k1),0.0_r8)* &
1167 & (tl_dtdr(i-1,j,k2)+ &
1168 & tl_dtdr(i ,j,k1)))- &
1169 & 0.5_r8*((0.5_r8+ &
1170 & sign(0.5_r8, drdx(i,j,k1)))* &
1171 & tl_drdx(i,j,k1)* &
1172 & (dtdr(i-1,j,k1)+dtdr(i,j,k2))+ &
1173 & (0.5_r8+ &
1174 & sign(0.5_r8,-drdx(i,j,k1)))* &
1175 & tl_drdx(i,j,k1)* &
1176 & (dtdr(i-1,j,k2)+dtdr(i,j,k1)))))- &
1177#ifdef TL_IOMS
1178 cff* &
1179 & (hz(i,j,k)+hz(i-1,j,k))* &
1180 & (dtdx(i,j,k1)- &
1181 & (max(drdx(i,j,k1),0.0_r8)* &
1182 & (dtdr(i-1,j,k1)+ &
1183 & dtdr(i ,j,k2))+ &
1184 & min(drdx(i,j,k1),0.0_r8)* &
1185 & (dtdr(i-1,j,k2)+ &
1186 & dtdr(i ,j,k1))))
1187#endif
1188 END DO
1189 END DO
1190 DO j=jstr,jend+1
1191 DO i=istr,iend
1192#ifdef DIFF_3DCOEF
1193# ifdef TS_U3ADV_SPLIT
1194 cff=0.5_r8*diff3d_v(i,j,k)*om_v(i,j)
1195# else
1196 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
1197 & om_v(i,j)
1198# endif
1199#else
1200 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
1201 & om_v(i,j)
1202#endif
1203!^ FE(i,j)=cff* &
1204!^ & (Hz(i,j,k)+Hz(i,j-1,k))* &
1205!^ & (dTde(i,j,k1)- &
1206!^ & 0.5_r8*(MAX(dRde(i,j,k1),0.0_r8)* &
1207!^ & (dTdr(i,j-1,k1)+ &
1208!^ & dTdr(i,j ,k2))+ &
1209!^ & MIN(dRde(i,j,k1),0.0_r8)* &
1210!^ & (dTdr(i,j-1,k2)+ &
1211!^ & dTdr(i,j ,k1))))
1212!^
1213 tl_fe(i,j)=cff* &
1214 & ((tl_hz(i,j,k)+tl_hz(i,j-1,k))* &
1215 & (dtde(i,j,k1)- &
1216 & 0.5_r8*(max(drde(i,j,k1),0.0_r8)* &
1217 & (dtdr(i,j-1,k1)+ &
1218 & dtdr(i,j ,k2))+ &
1219 & min(drde(i,j,k1),0.0_r8)* &
1220 & (dtdr(i,j-1,k2)+ &
1221 & dtdr(i,j ,k1))))+ &
1222 & (hz(i,j,k)+hz(i,j-1,k))* &
1223 & (tl_dtde(i,j,k1)- &
1224 & 0.5_r8*(max(drde(i,j,k1),0.0_r8)* &
1225 & (tl_dtdr(i,j-1,k1)+ &
1226 & tl_dtdr(i,j ,k2))+ &
1227 & min(drde(i,j,k1),0.0_r8)* &
1228 & (tl_dtdr(i,j-1,k2)+ &
1229 & tl_dtdr(i,j ,k1)))- &
1230 & 0.5_r8*((0.5_r8+ &
1231 & sign(0.5_r8, drde(i,j,k1)))* &
1232 & tl_drde(i,j,k1)* &
1233 & (dtdr(i,j-1,k1)+dtdr(i,j,k2))+ &
1234 & (0.5_r8+ &
1235 & sign(0.5_r8,-drde(i,j,k1)))* &
1236 & tl_drde(i,j,k1)* &
1237 & (dtdr(i,j-1,k2)+dtdr(i,j,k1)))))- &
1238#ifdef TL_IOMS
1239 & cff* &
1240 & (hz(i,j,k)+hz(i,j-1,k))* &
1241 & (dtde(i,j,k1)- &
1242 & (max(drde(i,j,k1),0.0_r8)* &
1243 & (dtdr(i,j-1,k1)+ &
1244 & dtdr(i,j ,k2))+ &
1245 & min(drde(i,j,k1),0.0_r8)* &
1246 & (dtdr(i,j-1,k2)+ &
1247 & dtdr(i,j ,k1))))
1248#endif
1249 END DO
1250 END DO
1251 IF (k.lt.n(ng)) THEN
1252 DO j=jstr,jend
1253 DO i=istr,iend
1254#ifdef DIFF_3DCOEF
1255# ifdef TS_U3ADV_SPLIT
1256 difx=0.125_r8*(diff3d_u(i,j,k )+diff3d_u(i+1,j,k )+ &
1257 & diff3d_u(i,j,k+1)+diff3d_u(i+1,j,k+1))
1258 dife=0.125_r8*(diff3d_v(i,j,k )+diff3d_v(i,j+1,k )+ &
1259 & diff3d_v(i,j,k+1)+diff3d_v(i,j+1,k+1))
1260# else
1261 difx=0.5_r8*diff3d_r(i,j,k)
1262 dife=difx
1263# endif
1264#else
1265 difx=0.5_r8*diff4(i,j,itrc)
1266 dife=difx
1267#endif
1268 cff1=max(drdx(i ,j,k1),0.0_r8)
1269 cff2=max(drdx(i+1,j,k2),0.0_r8)
1270 cff3=min(drdx(i ,j,k2),0.0_r8)
1271 cff4=min(drdx(i+1,j,k1),0.0_r8)
1272 tl_cff1=(0.5_r8+sign(0.5_r8, drdx(i ,j,k1)))* &
1273 & tl_drdx(i ,j,k1)
1274 tl_cff2=(0.5_r8+sign(0.5_r8, drdx(i+1,j,k2)))* &
1275 & tl_drdx(i+1,j,k2)
1276 tl_cff3=(0.5_r8+sign(0.5_r8,-drdx(i ,j,k2)))* &
1277 & tl_drdx(i ,j,k2)
1278 tl_cff4=(0.5_r8+sign(0.5_r8,-drdx(i+1,j,k1)))* &
1279 & tl_drdx(i+1,j,k1)
1280 cff=difx* &
1281 & (cff1*(cff1*dtdr(i,j,k2)-dtdx(i ,j,k1))+ &
1282 & cff2*(cff2*dtdr(i,j,k2)-dtdx(i+1,j,k2))+ &
1283 & cff3*(cff3*dtdr(i,j,k2)-dtdx(i ,j,k2))+ &
1284 & cff4*(cff4*dtdr(i,j,k2)-dtdx(i+1,j,k1)))
1285 tl_cff=difx* &
1286 & (tl_cff1*(cff1*dtdr(i ,j,k2)- &
1287 & dtdx(i ,j,k1))+ &
1288 & tl_cff2*(cff2*dtdr(i,j,k2)- &
1289 & dtdx(i+1,j,k2))+ &
1290 & tl_cff3*(cff3*dtdr(i,j,k2)- &
1291 & dtdx(i ,j,k2))+ &
1292 & tl_cff4*(cff4*dtdr(i,j,k2)- &
1293 & dtdx(i+1,j,k1))+ &
1294 & cff1*(tl_cff1*dtdr(i,j,k2)+ &
1295 & cff1*tl_dtdr(i,j,k2)- &
1296 & tl_dtdx(i ,j,k1))+ &
1297 & cff2*(tl_cff2*dtdr(i,j,k2)+ &
1298 & cff2*tl_dtdr(i,j,k2)- &
1299 & tl_dtdx(i+1,j,k2))+ &
1300 & cff3*(tl_cff3*dtdr(i,j,k2)+ &
1301 & cff3*tl_dtdr(i,j,k2)- &
1302 & tl_dtdx(i ,j,k2))+ &
1303 & cff4*(tl_cff4*dtdr(i,j,k2)+ &
1304 & cff4*tl_dtdr(i,j,k2)- &
1305 & tl_dtdx(i+1,j,k1)))- &
1306#ifdef TL_IOMS
1307 & difx* &
1308 & (cff1*(2.0_r8*cff1*dtdr(i,j,k2)- &
1309 & dtdx(i ,j,k1))- &
1310 & cff2*(2.0_r8*cff2*dtdr(i,j,k2)- &
1311 & dtdx(i+1,j,k2))- &
1312 & cff3*(2.0_r8*cff3*dtdr(i,j,k2)- &
1313 & dtdx(i ,j,k2))- &
1314 & cff4*(2.0_r8*cff4*dtdr(i,j,k2)- &
1315 & dtdx(i+1,j,k1)))
1316#endif
1317!
1318 cff1=max(drde(i,j ,k1),0.0_r8)
1319 cff2=max(drde(i,j+1,k2),0.0_r8)
1320 cff3=min(drde(i,j ,k2),0.0_r8)
1321 cff4=min(drde(i,j+1,k1),0.0_r8)
1322 tl_cff1=(0.5_r8+sign(0.5_r8, drde(i,j ,k1)))* &
1323 & tl_drde(i,j ,k1)
1324 tl_cff2=(0.5_r8+sign(0.5_r8, drde(i,j+1,k2)))* &
1325 & tl_drde(i,j+1,k2)
1326 tl_cff3=(0.5_r8+sign(0.5_r8,-drde(i,j ,k2)))* &
1327 & tl_drde(i,j ,k2)
1328 tl_cff4=(0.5_r8+sign(0.5_r8,-drde(i,j+1,k1)))* &
1329 & tl_drde(i,j+1,k1)
1330 cff=cff+ &
1331 & dife* &
1332 & (cff1*(cff1*dtdr(i,j,k2)-dtde(i,j ,k1))+ &
1333 & cff2*(cff2*dtdr(i,j,k2)-dtde(i,j+1,k2))+ &
1334 & cff3*(cff3*dtdr(i,j,k2)-dtde(i,j ,k2))+ &
1335 & cff4*(cff4*dtdr(i,j,k2)-dtde(i,j+1,k1)))
1336 tl_cff=tl_cff+ &
1337 & dife* &
1338 & (tl_cff1*(cff1*dtdr(i,j,k2)- &
1339 & dtde(i,j ,k1))+ &
1340 & tl_cff2*(cff2*dtdr(i,j,k2)- &
1341 & dtde(i,j+1,k2))+ &
1342 & tl_cff3*(cff3*dtdr(i,j,k2)- &
1343 & dtde(i,j ,k2))+ &
1344 & tl_cff4*(cff4*dtdr(i,j,k2)- &
1345 & dtde(i,j+1,k1))+ &
1346 & cff1*(tl_cff1*dtdr(i,j,k2)+ &
1347 & cff1*tl_dtdr(i,j,k2)- &
1348 & tl_dtde(i,j ,k1))+ &
1349 & cff2*(tl_cff2*dtdr(i,j,k2)+ &
1350 & cff2*tl_dtdr(i,j,k2)- &
1351 & tl_dtde(i,j+1,k2))+ &
1352 & cff3*(tl_cff3*dtdr(i,j,k2)+ &
1353 & cff3*tl_dtdr(i,j,k2)- &
1354 & tl_dtde(i,j ,k2))+ &
1355 & cff4*(tl_cff4*dtdr(i,j,k2)+ &
1356 & cff4*tl_dtdr(i,j,k2)- &
1357 & tl_dtde(i,j+1,k1)))- &
1358#ifdef TL_IOMS
1359 & dife* &
1360 & (cff1*(2.0_r8*cff1*dtdr(i,j,k2)- &
1361 & dtde(i,j ,k1))- &
1362 & cff2*(2.0_r8*cff2*dtdr(i,j,k2)- &
1363 & dtde(i,j+1,k2))- &
1364 & cff3*(2.0_r8*cff3*dtdr(i,j,k2)- &
1365 & dtde(i,j ,k2))- &
1366 & cff4*(2.0_r8*cff4*dtdr(i,j,k2)- &
1367 & dtde(i,j+1,k1)))
1368#endif
1369!^ FS(i,j,k2)=cff*FS(i,j,k2) ! recursive
1370!^ ! compute after TLM
1371!^
1372 tl_fs(i,j,k2)=tl_cff*fs(i,j,k2)+ &
1373 & cff*tl_fs(i,j,k2)
1374 fs(i,j,k2)=cff*fs(i,j,k2) ! recursive
1375#ifdef TL_IOMS
1376 tl_fs(i,j,k2)=tl_fs(i,j,k2)-fs(i,j,k2)
1377#endif
1378 END DO
1379 END DO
1380 END IF
1381!
1382! Time-step biharmonic, isopycnal diffusion term (m Tunits).
1383!
1384 DO j=jstr,jend
1385 DO i=istr,iend
1386!^ cff=dt(ng)*pm(i,j)*pn(i,j)* &
1387!^ & (FX(i+1,j )-FX(i,j)+ &
1388!^ & FE(i ,j+1)-FE(i,j))+ &
1389!^ & dt(ng)*(FS(i,j,k2)-FS(i,j,k1))
1390!^
1391 tl_cff=dt(ng)*pm(i,j)*pn(i,j)* &
1392 & (tl_fx(i+1,j)-tl_fx(i,j)+ &
1393 & tl_fe(i,j+1)-tl_fe(i,j))+ &
1394 & dt(ng)*(tl_fs(i,j,k2)-tl_fs(i,j,k1))
1395!^ t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)-cff
1396!^
1397 tl_t(i,j,k,nnew,itrc)=tl_t(i,j,k,nnew,itrc)-tl_cff
1398#ifdef DIAGNOSTICS_TS
1399!! DiaTwrk(i,j,k,itrc,iThdif)=-cff
1400#endif
1401 END DO
1402 END DO
1403 END IF
1404 END DO k_loop2
1405 END DO t_loop
1406!
1407 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.

◆ rp_t3dmix4_s_tile()

subroutine rp_t3dmix4_mod::rp_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 rp_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# ifdef TL_IOMS
292 & fx(i,j)
293# endif
294#elif defined TS_MIX_CLIMA
295 IF (ltracerclm(itrc,ng)) THEN
296 fx(i,j)=cff* &
297 & (hz(i,j,k)+hz(i-1,j,k))* &
298 & ((t(i ,j,k,nrhs,itrc)-tclm(i ,j,k,itrc))- &
299 & (t(i-1,j,k,nrhs,itrc)-tclm(i-1,j,k,itrc)))
300 tl_fx(i,j)=cff* &
301 & ((tl_hz(i,j,k)+tl_hz(i-1,j,k))* &
302 & ((t(i ,j,k,nrhs,itrc)- &
303 & tclm(i ,j,k,itrc))- &
304 & (t(i-1,j,k,nrhs,itrc)- &
305 & tclm(i-1,j,k,itrc)))+ &
306 & (hz(i,j,k)+hz(i-1,j,k))* &
307 & (tl_t(i ,j,k,nrhs,itrc)- &
308 & tl_t(i-1,j,k,nrhs,itrc)))- &
309# ifdef TL_IOMS
310 & fx(i,j)
311# endif
312 ELSE
313 fx(i,j)=cff* &
314 & (hz(i,j,k)+hz(i-1,j,k))* &
315 & (t(i ,j,k,nrhs,itrc)- &
316 & t(i-1,j,k,nrhs,itrc))
317 tl_fx(i,j)=cff* &
318 & ((tl_hz(i,j,k)+tl_hz(i-1,j,k))* &
319 & (t(i ,j,k,nrhs,itrc)- &
320 & t(i-1,j,k,nrhs,itrc))+ &
321 & (hz(i,j,k)+hz(i-1,j,k))* &
322 & (tl_t(i ,j,k,nrhs,itrc)- &
323 & tl_t(i-1,j,k,nrhs,itrc)))- &
324# ifdef TL_IOMS
325 & fx(i,j)
326# endif
327 END IF
328#else
329 fx(i,j)=cff* &
330 & (hz(i,j,k)+hz(i-1,j,k))* &
331 & (t(i ,j,k,nrhs,itrc)- &
332 & t(i-1,j,k,nrhs,itrc))
333 tl_fx(i,j)=cff* &
334 & ((tl_hz(i,j,k)+tl_hz(i-1,j,k))* &
335 & (t(i ,j,k,nrhs,itrc)- &
336 & t(i-1,j,k,nrhs,itrc))+ &
337 & (hz(i,j,k)+hz(i-1,j,k))* &
338 & (tl_t(i ,j,k,nrhs,itrc)- &
339 & tl_t(i-1,j,k,nrhs,itrc)))- &
340# ifdef TL_IOMS
341 & fx(i,j)
342# endif
343#endif
344 END DO
345 END DO
346 DO j=jmin,jmax+1
347 DO i=imin,imax
348#ifdef DIFF_3DCOEF
349# ifdef TS_U3ADV_SPLIT_NOT_YET
350 cff=0.5_r8*diff3d_v(i,j,k)*pnom_v(i,j)
351# else
352 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
353 & pnom_v(i,j)
354# endif
355#else
356 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
357 & pnom_v(i,j)
358#endif
359#ifdef MASKING
360 cff=cff*vmask(i,j)
361#endif
362#ifdef WET_DRY_NOT_YET
363 cff=cff*vmask_wet(i,j)
364#endif
365#if defined TS_MIX_STABILITY
366 fe(i,j)=cff* &
367 & (hz(i,j,k)+hz(i,j-1,k))* &
368 & (0.75_r8*(t(i,j ,k,nrhs,itrc)- &
369 & t(i,j-1,k,nrhs,itrc))+ &
370 & 0.25_r8*(t(i,j ,k,nstp,itrc)- &
371 & t(i,j-1,k,nstp,itrc)))
372 tl_fe(i,j)=cff* &
373 & ((tl_hz(i,j,k)+tl_hz(i,j-1,k))* &
374 & (0.75_r8*(t(i,j ,k,nrhs,itrc)- &
375 & t(i,j-1,k,nrhs,itrc))+ &
376 & 0.25_r8*(t(i,j ,k,nstp,itrc)- &
377 & t(i,j-1,k,nstp,itrc)))+ &
378 & (hz(i,j,k)+hz(i,j-1,k))* &
379 & (0.75_r8*(tl_t(i,j ,k,nrhs,itrc)- &
380 & tl_t(i,j-1,k,nrhs,itrc))+ &
381 & 0.25_r8*(tl_t(i,j ,k,nstp,itrc)- &
382 & tl_t(i,j-1,k,nstp,itrc))))- &
383# ifdef TL_IOMS
384 & fe(i,j)
385# endif
386#elif defined TS_MIX_CLIMA
387 IF (ltracerclm(itrc,ng)) THEN
388 fe(i,j)=cff* &
389 & (hz(i,j,k)+hz(i,j-1,k))* &
390 & ((t(i,j ,k,nrhs,itrc)-tclm(i,j ,k,itrc))- &
391 & (t(i,j-1,k,nrhs,itrc)-tclm(i,j-1,k,itrc)))
392 tl_fe(i,j)=cff* &
393 & ((tl_hz(i,j,k)+tl_hz(i,j-1,k))* &
394 & ((t(i,j ,k,nrhs,itrc)- &
395 & tclm(i,j ,k,itrc))- &
396 & (t(i,j-1,k,nrhs,itrc)- &
397 & tclm(i,j-1,k,itrc)))+ &
398 & (hz(i,j,k)+hz(i,j-1,k))* &
399 & (tl_t(i,j ,k,nrhs,itrc)- &
400 & tl_t(i,j-1,k,nrhs,itrc)))- &
401# ifdef TL_IOMS
402 & fe(i,j)
403# endif
404 ELSE
405 fe(i,j)=cff* &
406 & (hz(i,j,k)+hz(i,j-1,k))* &
407 & (t(i,j ,k,nrhs,itrc)- &
408 & t(i,j-1,k,nrhs,itrc))
409 tl_fe(i,j)=cff* &
410 & ((tl_hz(i,j,k)+tl_hz(i,j-1,k))* &
411 & (t(i,j ,k,nrhs,itrc)- &
412 & t(i,j-1,k,nrhs,itrc))+ &
413 & (hz(i,j,k)+hz(i,j-1,k))* &
414 & (tl_t(i,j ,k,nrhs,itrc)- &
415 & tl_t(i,j-1,k,nrhs,itrc)))- &
416# ifdef TL_IOMS
417 & fe(i,j)
418# endif
419 END IF
420#else
421 fe(i,j)=cff* &
422 & (hz(i,j,k)+hz(i,j-1,k))* &
423 & (t(i,j ,k,nrhs,itrc)- &
424 & t(i,j-1,k,nrhs,itrc))
425 tl_fe(i,j)=cff* &
426 & ((tl_hz(i,j,k)+tl_hz(i,j-1,k))* &
427 & (t(i,j ,k,nrhs,itrc)- &
428 & t(i,j-1,k,nrhs,itrc))+ &
429 & (hz(i,j,k)+hz(i,j-1,k))* &
430 & (tl_t(i,j ,k,nrhs,itrc)- &
431 & tl_t(i,j-1,k,nrhs,itrc)))- &
432# ifdef TL_IOMS
433 & fe(i,j)
434# endif
435#endif
436 END DO
437 END DO
438!
439! Compute first harmonic operator and multiply by the metrics of the
440! second harmonic operator.
441!
442 DO j=jmin,jmax
443 DO i=imin,imax
444 cff=1.0_r8/hz(i,j,k)
445 tl_cff=-cff*cff*tl_hz(i,j,k)+ &
446#ifdef TL_IOMS
447 & 2.0_r8*cff
448#endif
449 lapt(i,j)=pm(i,j)*pn(i,j)*cff* &
450 & (fx(i+1,j)-fx(i,j)+ &
451 & fe(i,j+1)-fe(i,j))
452 tl_lapt(i,j)=pm(i,j)*pn(i,j)* &
453 & (tl_cff* &
454 & (fx(i+1,j)-fx(i,j)+ &
455 & fe(i,j+1)-fe(i,j))+ &
456 & cff* &
457 & (tl_fx(i+1,j)-tl_fx(i,j)+ &
458 & tl_fe(i,j+1)-tl_fe(i,j)))- &
459#ifdef TL_IOMS
460 & lapt(i,j)
461#endif
462 END DO
463 END DO
464!
465! Apply boundary conditions (except periodic; closed or gradient)
466! to the first harmonic operator.
467!
468 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
469 IF (domain(ng)%Western_Edge(tile)) THEN
470 IF (tl_lbc(iwest,istvar(itrc),ng)%closed) THEN
471 DO j=jmin,jmax
472 lapt(istr-1,j)=0.0_r8
473 tl_lapt(istr-1,j)=0.0_r8
474 END DO
475 ELSE
476 DO j=jmin,jmax
477 lapt(istr-1,j)=lapt(istr,j)
478 tl_lapt(istr-1,j)=tl_lapt(istr,j)
479 END DO
480 END IF
481 END IF
482 END IF
483!
484 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
485 IF (domain(ng)%Eastern_Edge(tile)) THEN
486 IF (tl_lbc(ieast,istvar(itrc),ng)%closed) THEN
487 DO j=jmin,jmax
488 lapt(iend+1,j)=0.0_r8
489 tl_lapt(iend+1,j)=0.0_r8
490 END DO
491 ELSE
492 DO j=jmin,jmax
493 lapt(iend+1,j)=lapt(iend,j)
494 tl_lapt(iend+1,j)=tl_lapt(iend,j)
495 END DO
496 END IF
497 END IF
498 END IF
499!
500 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
501 IF (domain(ng)%Southern_Edge(tile)) THEN
502 IF (tl_lbc(isouth,istvar(itrc),ng)%closed) THEN
503 DO i=imin,imax
504 lapt(i,jstr-1)=0.0_r8
505 tl_lapt(i,jstr-1)=0.0_r8
506 END DO
507 ELSE
508 DO i=imin,imax
509 lapt(i,jstr-1)=lapt(i,jstr)
510 tl_lapt(i,jstr-1)=tl_lapt(i,jstr)
511 END DO
512 END IF
513 END IF
514 END IF
515!
516 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
517 IF (domain(ng)%Northern_Edge(tile)) THEN
518 IF (tl_lbc(inorth,istvar(itrc),ng)%closed) THEN
519 DO i=imin,imax
520 lapt(i,jend+1)=0.0_r8
521 tl_lapt(i,jend+1)=0.0_r8
522 END DO
523 ELSE
524 DO i=imin,imax
525 lapt(i,jend+1)=lapt(i,jend)
526 tl_lapt(i,jend+1)=tl_lapt(i,jend)
527 END DO
528 END IF
529 END IF
530 END IF
531!
532! Compute FX=d(LapT)/d(xi) and FE=d(LapT)/d(eta) terms.
533!
534 DO j=jstr,jend
535 DO i=istr,iend+1
536#ifdef DIFF_3DCOEF
537# ifdef TS_U3ADV_SPLIT_NOT_YET
538 cff=0.5_r8*diff3d_u(i,j,k)*pmon_u(i,j)
539# else
540 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
541 & pmon_u(i,j)
542# endif
543#else
544 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
545 & pmon_u(i,j)
546#endif
547!^ FX(i,j)=cff* &
548!^ & (Hz(i,j,k)+Hz(i-1,j,k))* &
549!^ & (LapT(i,j)-LapT(i-1,j))
550!^
551 tl_fx(i,j)=cff* &
552 & ((tl_hz(i,j,k)+tl_hz(i-1,j,k))* &
553 & (lapt(i,j)-lapt(i-1,j))+ &
554 & (hz(i,j,k)+hz(i-1,j,k))* &
555 & (tl_lapt(i,j)-tl_lapt(i-1,j)))- &
556# ifdef TL_IOMS
557 & cff* &
558 & (hz(i,j,k)+hz(i-1,j,k))* &
559 & (lapt(i,j)-lapt(i-1,j))
560# endif
561#ifdef MASKING
562!^ FX(i,j)=FX(i,j)*umask(i,j)
563!^
564 tl_fx(i,j)=tl_fx(i,j)*umask(i,j)
565#endif
566#ifdef WET_DRY_NOT_YET
567 fx(i,j)=fx(i,j)*umask_wet(i,j)
568#endif
569 END DO
570 END DO
571 DO j=jstr,jend+1
572 DO i=istr,iend
573#ifdef DIFF_3DCOEF
574# ifdef TS_U3ADV_SPLIT_NOT_YET
575 cff=0.5_r8*diff3d_v(i,j,k)*pnom_v(i,j)
576# else
577 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
578 & pnom_v(i,j)
579# endif
580#else
581 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
582 & pnom_v(i,j)
583#endif
584!^ FE(i,j)=cff* &
585!^ & (Hz(i,j,k)+Hz(i,j-1,k))* &
586!^ & (LapT(i,j)-LapT(i,j-1))
587!^
588 tl_fe(i,j)=cff* &
589 & ((tl_hz(i,j,k)+tl_hz(i,j-1,k))* &
590 & (lapt(i,j)-lapt(i,j-1))+ &
591 & (hz(i,j,k)+hz(i,j-1,k))* &
592 & (tl_lapt(i,j)-tl_lapt(i,j-1)))- &
593#ifdef TL_IOMS
594 & cff* &
595 & (hz(i,j,k)+hz(i,j-1,k))* &
596 & (lapt(i,j)-lapt(i,j-1))
597#endif
598#ifdef MASKING
599!^ FE(i,j)=FE(i,j)*vmask(i,j)
600!^
601 tl_fe(i,j)=tl_fe(i,j)*vmask(i,j)
602#endif
603#ifdef WET_DRY_NOT_YET
604 fe(i,j)=fe(i,j)*vmask_wet(i,j)
605#endif
606 END DO
607 END DO
608!
609! Time-step biharmonic, S-surfaces diffusion term (m Tunits).
610!
611 DO j=jstr,jend
612 DO i=istr,iend
613!^ cff=dt(ng)*pm(i,j)*pn(i,j)* &
614!^ & (FX(i+1,j)-FX(i,j)+ &
615!^ & FE(i,j+1)-FE(i,j))
616!^
617 tl_cff=dt(ng)*pm(i,j)*pn(i,j)* &
618 & (tl_fx(i+1,j)-tl_fx(i,j)+ &
619 & tl_fe(i,j+1)-tl_fe(i,j))
620!^ t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)-cff
621!^
622 tl_t(i,j,k,nnew,itrc)=tl_t(i,j,k,nnew,itrc)-tl_cff
623#ifdef DIAGNOSTICS_TS
624!! DiaTwrk(i,j,k,itrc,iThdif)=-cff
625#endif
626 END DO
627 END DO
628 END DO
629 END DO
630!
631 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.