132
133
137
138
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# ifdef DIAGNOSTICS_TS
176 real(r8), intent(inout) :: DiaTwrk(LBi:,LBj:,:,:,:)
177# endif
178 real(r8), intent(inout) :: ad_Hz(LBi:,LBj:,:)
179 real(r8), intent(inout) :: ad_z_r(LBi:,LBj:,:)
180 real(r8), intent(inout) :: ad_pden(LBi:,LBj:,:)
181 real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:)
182#else
183# ifdef MASKING
184 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
185 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
186# endif
187# ifdef WET_DRY_NOT_YET
188 real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
189 real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
190# endif
191# ifdef DIFF_3DCOEF
192# ifdef TS_U3ADV_SPLIT
193 real(r8), intent(in) :: diff3d_u(LBi:UBi,LBj:UBj,N(ng))
194 real(r8), intent(in) :: diff3d_v(LBi:UBi,LBj:UBj,N(ng))
195# else
196 real(r8), intent(in) :: diff3d_r(LBi:UBi,LBj:UBj,N(ng))
197# endif
198# else
199 real(r8), intent(in) :: diff4(LBi:UBi,LBj:UBj,NT(ng))
200# endif
201 real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
202 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
203 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
204 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
205 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
206 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
207 real(r8), intent(in) :: pden(LBi:UBi,LBj:UBj,N(ng))
208 real(r8), intent(in) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
209# ifdef TS_MIX_CLIMA
210 real(r8), intent(in) :: tclm(LBi:UBi,LBj:UBj,N(ng),NT(ng))
211# endif
212# ifdef DIAGNOSTICS_TS
213
214
215# endif
216 real(r8), intent(inout) :: ad_Hz(LBi:UBi,LBj:UBj,N(ng))
217 real(r8), intent(inout) :: ad_z_r(LBi:UBi,LBj:UBj,N(ng))
218 real(r8), intent(inout) :: ad_pden(LBi:UBi,LBj:UBj,N(ng))
219 real(r8), intent(inout) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
220#endif
221
222
223
224 integer :: Imin, Imax, Jmin, Jmax
225 integer :: i, itrc, j, k, kk, kt, k1, k1b, k2, k2b
226
227 real(r8), parameter :: eps = 0.5_r8
228 real(r8), parameter :: small = 1.0e-14_r8
229 real(r8), parameter :: slope_max = 0.0001_r8
230 real(r8), parameter :: strat_min = 0.1_r8
231
232 real(r8) :: cff, cff1, cff2, cff3, cff4, dife, difx
233 real(r8) :: ad_cff, ad_cff1, ad_cff2, ad_cff3, ad_cff4
234 real(r8) :: adfac, adfac1, adfac2, adfac3, adfac4
235
236 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: LapT
237
238 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: ad_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) :: ad_FE
244 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: ad_FX
245
246 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: FS
247 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: FS1
248 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dRde
249 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dRdx
250 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dTde
251 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dTdr
252 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dTdx
253
254 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: ad_FS
255 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: ad_dRde
256 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: ad_dRdx
257 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: ad_dTde
258 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: ad_dTdr
259 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: ad_dTdx
260
261#include "set_bounds.h"
262
263
264
265
266
267 ad_cff=0.0_r8
268 ad_cff1=0.0_r8
269 ad_cff2=0.0_r8
270 ad_cff3=0.0_r8
271 ad_cff4=0.0_r8
272
273 ad_fe(imins:imaxs,jmins:jmaxs)=0.0_r8
274 ad_fx(imins:imaxs,jmins:jmaxs)=0.0_r8
275
276 ad_fs(imins:imaxs,jmins:jmaxs,1:2)=0.0_r8
277
278 ad_drde(imins:imaxs,jmins:jmaxs,1:2)=0.0_r8
279 ad_drdx(imins:imaxs,jmins:jmaxs,1:2)=0.0_r8
280 ad_dtde(imins:imaxs,jmins:jmaxs,1:2)=0.0_r8
281 ad_dtdr(imins:imaxs,jmins:jmaxs,1:2)=0.0_r8
282 ad_dtdx(imins:imaxs,jmins:jmaxs,1:2)=0.0_r8
283
284 ad_lapt(imins:imaxs,jmins:jmaxs,1:
n(ng))=0.0_r8
285
286
287
288
289
290
291
292
293
295 imin=istr-1
296 imax=iend+1
297 ELSE
298 imin=max(istr-1,1)
299 imax=min(iend+1,
lm(ng))
300 END IF
302 jmin=jstr-1
303 jmax=jend+1
304 ELSE
305 jmin=max(jstr-1,1)
306 jmax=min(jend+1,
mm(ng))
307 END IF
308
309
310
311
312
313
314
315
316
317
318#ifdef TS_MIX_STABILITY
319
320
321
322#endif
323
324 t_loop :
DO itrc=1,
nt(ng)
325 k2=1
326 k_loop1 :
DO k=0,
n(ng)
327 k1=k2
328 k2=3-k1
330 DO j=jmin,jmax
331 DO i=imin,imax+1
332 cff=0.5_r8*(pm(i,j)+pm(i-1,j))
333#ifdef MASKING
334 cff=cff*umask(i,j)
335#endif
336#ifdef WET_DRY_NOT_YET
337 cff=cff*umask_wet(i,j)
338#endif
339 drdx(i,j,k2)=cff*(pden(i ,j,k+1)- &
340 & pden(i-1,j,k+1))
341#if defined TS_MIX_STABILITY
342 dtdx(i,j,k2)=cff*(0.75_r8*(t(i ,j,k+1,nrhs,itrc)- &
343 & t(i-1,j,k+1,nrhs,itrc))+ &
344 & 0.25_r8*(t(i ,j,k+1,nstp,itrc)- &
345 & t(i-1,j,k+1,nstp,itrc)))
346#elif defined TS_MIX_CLIMA
348 dtdx(i,j,k2)=cff*((t(i ,j,k+1,nrhs,itrc)- &
349 & tclm(i ,j,k+1,itrc))- &
350 & (t(i-1,j,k+1,nrhs,itrc)- &
351 & tclm(i-1,j,k+1,itrc)))
352 ELSE
353 dtdx(i,j,k2)=cff*(t(i ,j,k+1,nrhs,itrc)- &
354 & t(i-1,j,k+1,nrhs,itrc))
355 END IF
356#else
357 dtdx(i,j,k2)=cff*(t(i ,j,k+1,nrhs,itrc)- &
358 & t(i-1,j,k+1,nrhs,itrc))
359#endif
360 END DO
361 END DO
362 DO j=jmin,jmax+1
363 DO i=imin,imax
364 cff=0.5_r8*(pn(i,j)+pn(i,j-1))
365#ifdef MASKING
366 cff=cff*vmask(i,j)
367#endif
368#ifdef WET_DRY_NOT_YET
369 cff=cff*vmask_wet(i,j)
370#endif
371 drde(i,j,k2)=cff*(pden(i,j ,k+1)- &
372 & pden(i,j-1,k+1))
373#if defined TS_MIX_STABILITY
374 dtde(i,j,k2)=cff*(0.75_r8*(t(i,j ,k+1,nrhs,itrc)- &
375 & t(i,j-1,k+1,nrhs,itrc))+ &
376 & 0.25_r8*(t(i,j ,k+1,nstp,itrc)- &
377 & t(i,j-1,k+1,nstp,itrc)))
378#elif defined TS_MIX_CLIMA
380 dtde(i,j,k2)=cff*((t(i,j ,k+1,nrhs,itrc)- &
381 & tclm(i,j ,k+1,itrc))- &
382 & (t(i,j-1,k+1,nrhs,itrc)- &
383 & tclm(i,j-1,k+1,itrc)))
384 ELSE
385 dtde(i,j,k2)=cff*(t(i,j ,k+1,nrhs,itrc)- &
386 & t(i,j-1,k+1,nrhs,itrc))
387 END IF
388#else
389 dtde(i,j,k2)=cff*(t(i,j ,k+1,nrhs,itrc)- &
390 & t(i,j-1,k+1,nrhs,itrc))
391#endif
392 END DO
393 END DO
394 END IF
395 IF ((k.eq.0).or.(k.eq.
n(ng)))
THEN
396 DO j=-1+jmin,jmax+1
397 DO i=-1+imin,imax+1
398 dtdr(i,j,k2)=0.0_r8
399 fs(i,j,k2)=0.0_r8
400 END DO
401 END DO
402 ELSE
403 DO j=-1+jmin,jmax+1
404 DO i=-1+imin,imax+1
405#if defined TS_MIX_MAX_SLOPE
406 cff1=sqrt(drdx(i,j,k2)**2+drdx(i+1,j,k2)**2+ &
407 & drdx(i,j,k1)**2+drdx(i+1,j,k1)**2+ &
408 & drde(i,j,k2)**2+drde(i,j+1,k2)**2+ &
409 & drde(i,j,k1)**2+drde(i,j+1,k1)**2)
410 cff2=0.25_r8*slope_max* &
411 & (z_r(i,j,k+1)-z_r(i,j,k))*cff1
412 cff3=max(pden(i,j,k)-pden(i,j,k+1),small)
413 cff4=max(cff2,cff3)
414 cff=-1.0_r8/cff4
415#elif defined TS_MIX_MIN_STRAT
416 cff1=max(pden(i,j,k)-pden(i,j,k+1), &
417 & strat_min*(z_r(i,j,k+1)-z_r(i,j,k)))
418 cff=-1.0_r8/cff1
419#else
420 cff1=max(pden(i,j,k)-pden(i,j,k+1),eps)
421 cff=-1.0_r8/cff1
422#endif
423#if defined TS_MIX_STABILITY
424 dtdr(i,j,k2)=cff*(0.75_r8*(t(i,j,k+1,nrhs,itrc)- &
425 & t(i,j,k ,nrhs,itrc))+ &
426 & 0.25_r8*(t(i,j,k+1,nstp,itrc)- &
427 & t(i,j,k ,nstp,itrc)))
428#elif defined TS_MIX_CLIMA
430 dtdr(i,j,k2)=cff*((t(i,j,k+1,nrhs,itrc)- &
431 & tclm(i,j,k+1,itrc))- &
432 & (t(i,j,k ,nrhs,itrc)- &
433 & tclm(i,j,k ,itrc)))
434 ELSE
435 dtdr(i,j,k2)=cff*(t(i,j,k+1,nrhs,itrc)- &
436 & t(i,j,k ,nrhs,itrc))
437 END IF
438#else
439 dtdr(i,j,k2)=cff*(t(i,j,k+1,nrhs,itrc)- &
440 & t(i,j,k ,nrhs,itrc))
441#endif
442 fs(i,j,k2)=cff*(z_r(i,j,k+1)- &
443 & z_r(i,j,k ))
444 END DO
445 END DO
446 END IF
447 IF (k.gt.0) THEN
448 DO j=jmin,jmax
449 DO i=imin,imax+1
450#ifdef DIFF_3DCOEF
451# ifdef TS_U3ADV_SPLIT
452 cff=0.5_r8*diff3d_u(i,j,k)*on_u(i,j)
453# else
454 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
455 & on_u(i,j)
456# endif
457#else
458 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
459 & on_u(i,j)
460#endif
461 fx(i,j)=cff* &
462 & (hz(i,j,k)+hz(i-1,j,k))* &
463 & (dtdx(i,j,k1)- &
464 & 0.5_r8*(max(drdx(i,j,k1),0.0_r8)* &
465 & (dtdr(i-1,j,k1)+ &
466 & dtdr(i ,j,k2))+ &
467 & min(drdx(i,j,k1),0.0_r8)* &
468 & (dtdr(i-1,j,k2)+ &
469 & dtdr(i ,j,k1))))
470 END DO
471 END DO
472 DO j=jmin,jmax+1
473 DO i=imin,imax
474#ifdef DIFF_3DCOEF
475# ifdef TS_U3ADV_SPLIT
476 cff=0.5_r8*diff3d_v(i,j,k)*om_v(i,j)
477# else
478 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
479 & om_v(i,j)
480# endif
481#else
482 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
483 & om_v(i,j)
484#endif
485 fe(i,j)=cff* &
486 & (hz(i,j,k)+hz(i,j-1,k))* &
487 & (dtde(i,j,k1)- &
488 & 0.5_r8*(max(drde(i,j,k1),0.0_r8)* &
489 & (dtdr(i,j-1,k1)+ &
490 & dtdr(i,j ,k2))+ &
491 & min(drde(i,j,k1),0.0_r8)* &
492 & (dtdr(i,j-1,k2)+ &
493 & dtdr(i,j ,k1))))
494 END DO
495 END DO
497 DO j=jmin,jmax
498 DO i=imin,imax
499#ifdef DIFF_3DCOEF
500# ifdef TS_U3ADV_SPLIT
501 difx=0.125_r8*(diff3d_u(i,j,k )+diff3d_u(i+1,j,k )+ &
502 & diff3d_u(i,j,k+1)+diff3d_u(i+1,j,k+1))
503 dife=0.125_r8*(diff3d_v(i,j,k )+diff3d_v(i,j+1,k )+ &
504 & diff3d_v(i,j,k+1)+diff3d_v(i,j+1,k+1))
505# else
506 difx=0.5_r8*diff3d_r(i,j,k)
507 dife=difx
508# endif
509#else
510 difx=0.5_r8*diff4(i,j,itrc)
511 dife=difx
512#endif
513 cff1=max(drdx(i ,j,k1),0.0_r8)
514 cff2=max(drdx(i+1,j,k2),0.0_r8)
515 cff3=min(drdx(i ,j,k2),0.0_r8)
516 cff4=min(drdx(i+1,j,k1),0.0_r8)
517 cff=difx* &
518 & (cff1*(cff1*dtdr(i,j,k2)-dtdx(i ,j,k1))+ &
519 & cff2*(cff2*dtdr(i,j,k2)-dtdx(i+1,j,k2))+ &
520 & cff3*(cff3*dtdr(i,j,k2)-dtdx(i ,j,k2))+ &
521 & cff4*(cff4*dtdr(i,j,k2)-dtdx(i+1,j,k1)))
522 cff1=max(drde(i,j ,k1),0.0_r8)
523 cff2=max(drde(i,j+1,k2),0.0_r8)
524 cff3=min(drde(i,j ,k2),0.0_r8)
525 cff4=min(drde(i,j+1,k1),0.0_r8)
526 cff=cff+ &
527 & dife* &
528 & (cff1*(cff1*dtdr(i,j,k2)-dtde(i,j ,k1))+ &
529 & cff2*(cff2*dtdr(i,j,k2)-dtde(i,j+1,k2))+ &
530 & cff3*(cff3*dtdr(i,j,k2)-dtde(i,j ,k2))+ &
531 & cff4*(cff4*dtdr(i,j,k2)-dtde(i,j+1,k1)))
532 fs(i,j,k2)=cff*fs(i,j,k2)
533 END DO
534 END DO
535 END IF
536
537
538
539
540
541 DO j=jmin,jmax
542 DO i=imin,imax
543 cff=pm(i,j)*pn(i,j)
544 cff1=1.0_r8/hz(i,j,k)
545 lapt(i,j,k)=cff1*(cff* &
546 & (fx(i+1,j)-fx(i,j)+ &
547 & fe(i,j+1)-fe(i,j))+ &
548 & (fs(i,j,k2)-fs(i,j,k1)))
549 END DO
550 END DO
551 END IF
552 END DO k_loop1
553
554
555
556
558 IF (
domain(ng)%Western_Edge(tile))
THEN
561 DO j=jmin,jmax
562 lapt(istr-1,j,k)=0.0_r8
563 END DO
564 END DO
565 ELSE
567 DO j=jmin,jmax
568 lapt(istr-1,j,k)=lapt(istr,j,k)
569 END DO
570 END DO
571 END IF
572 END IF
573 END IF
574
576 IF (
domain(ng)%Eastern_Edge(tile))
THEN
579 DO j=jmin,jmax
580 lapt(iend+1,j,k)=0.0_r8
581 END DO
582 END DO
583 ELSE
585 DO j=jmin,jmax
586 lapt(iend+1,j,k)=lapt(iend,j,k)
587 END DO
588 END DO
589 END IF
590 END IF
591 END IF
592
594 IF (
domain(ng)%Southern_Edge(tile))
THEN
597 DO i=imin,imax
598 lapt(i,jstr-1,k)=0.0_r8
599 END DO
600 END DO
601 ELSE
603 DO i=imin,imax
604 lapt(i,jstr-1,k)=lapt(i,jstr,k)
605 END DO
606 END DO
607 END IF
608 END IF
609 END IF
610
612 IF (
domain(ng)%Northern_Edge(tile))
THEN
615 DO i=imin,imax
616 lapt(i,jend+1,k)=0.0_r8
617 END DO
618 END DO
619 ELSE
621 DO i=imin,imax
622 lapt(i,jend+1,k)=lapt(i,jend,k)
623 END DO
624 END DO
625 END IF
626 END IF
627 END IF
628
631 IF (
domain(ng)%SouthWest_Corner(tile))
THEN
633 lapt(istr-1,jstr-1,k)=0.5_r8* &
634 & (lapt(istr ,jstr-1,k)+ &
635 & lapt(istr-1,jstr ,k))
636 END DO
637 END IF
638 END IF
639
642 IF (
domain(ng)%SouthEast_Corner(tile))
THEN
644 lapt(iend+1,jstr-1,k)=0.5_r8* &
645 & (lapt(iend ,jstr-1,k)+ &
646 & lapt(iend+1,jstr ,k))
647 END DO
648 END IF
649 END IF
650
653 IF (
domain(ng)%NorthWest_Corner(tile))
THEN
655 lapt(istr-1,jend+1,k)=0.5_r8* &
656 & (lapt(istr ,jend+1,k)+ &
657 & lapt(istr-1,jend ,k))
658 END DO
659 END IF
660 END IF
661
664 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
666 lapt(iend+1,jend+1,k)=0.5_r8* &
667 & (lapt(iend ,jend+1,k)+ &
668 & lapt(iend+1,jend ,k))
669 END DO
670 END IF
671 END IF
672
673
674
675 k1=2
676 k2=1
678
679
680
681
682
683
684
685
686
687 k1=k2
688 k2=3-k1
689 END DO
690
691
692
693
694 k_loop2:
DO k=
n(ng),0,-1
695 k2b=1
696 DO kk=0,k
697 k1b=k2b
698 k2b=3-k1b
699
700
701
702
703 IF (kk.lt.
n(ng))
THEN
704 DO j=jstr,jend
705 DO i=istr,iend+1
706 cff=0.5_r8*(pm(i,j)+pm(i-1,j))
707#ifdef MASKING
708 cff=cff*umask(i,j)
709#endif
710#ifdef WET_DRY_NOT_YET
711 cff=cff*umask_wet(i,j)
712#endif
713 drdx(i,j,k2b)=cff*(pden(i ,j,kk+1)- &
714 & pden(i-1,j,kk+1))
715 dtdx(i,j,k2b)=cff*(lapt(i ,j,kk+1)- &
716 & lapt(i-1,j,kk+1))
717 END DO
718 END DO
719 IF (kk.eq.0) THEN
720 DO j=jstr,jend
721 DO i=istr,iend+1
722 drdx(i,j,k1b)=0.0_r8
723 dtdx(i,j,k1b)=0.0_r8
724 END DO
725 END DO
726 END IF
727 DO j=jstr,jend+1
728 DO i=istr,iend
729 cff=0.5_r8*(pn(i,j)+pn(i,j-1))
730#ifdef MASKING
731 cff=cff*vmask(i,j)
732#endif
733#ifdef WET_DRY_NOT_YET
734 cff=cff*vmask_wet(i,j)
735#endif
736 drde(i,j,k2b)=cff*(pden(i,j ,kk+1)- &
737 & pden(i,j-1,kk+1))
738 dtde(i,j,k2b)=cff*(lapt(i,j ,kk+1)- &
739 & lapt(i,j-1,kk+1))
740 END DO
741 END DO
742 IF (kk.eq.0) THEN
743 DO j=jstr,jend+1
744 DO i=istr,iend
745 drde(i,j,k1b)=0.0_r8
746 dtde(i,j,k1b)=0.0_r8
747 END DO
748 END DO
749 END IF
750 END IF
751 IF ((kk.eq.0).or.(kk.eq.
n(ng)))
THEN
752 DO j=jstr-1,jend+1
753 DO i=istr-1,iend+1
754 dtdr(i,j,k2b)=0.0_r8
755 fs(i,j,k2b)=0.0_r8
756 END DO
757 END DO
758 IF (kk.eq.0) THEN
759 DO j=jstr-1,jend+1
760 DO i=istr-1,iend+1
761 dtdr(i,j,k1b)=0.0_r8
762 fs(i,j,k1b)=0.0_r8
763 END DO
764 END DO
765 END IF
766 ELSE
767 DO j=jstr-1,jend+1
768 DO i=istr-1,iend+1
769#if defined TS_MIX_MAX_SLOPE
770 cff1=sqrt(drdx(i,j,k2b)**2+drdx(i+1,j,k2b)**2+ &
771 & drdx(i,j,k1b)**2+drdx(i+1,j,k1b)**2+ &
772 & drde(i,j,k2b)**2+drde(i,j+1,k2b)**2+ &
773 & drde(i,j,k1b)**2+drde(i,j+1,k1b)**2)
774 cff2=0.25_r8*slope_max* &
775 & (z_r(i,j,kk+1)-z_r(i,j,kk))*cff1
776 cff3=max(pden(i,j,kk)-pden(i,j,kk+1),small)
777 cff4=max(cff2,cff3)
778 cff=-1.0_r8/cff4
779#elif defined TS_MIX_MIN_STRAT
780 cff1=max(pden(i,j,kk)-pden(i,j,kk+1), &
781 & strat_min*(z_r(i,j,kk+1)-z_r(i,j,kk)))
782 cff=-1.0_r8/cff1
783#else
784 cff1=max(pden(i,j,kk)-pden(i,j,kk+1),eps)
785 cff=-1.0_r8/cff1
786#endif
787 dtdr(i,j,k2b)=cff*(lapt(i,j,kk+1)- &
788 & lapt(i,j,kk ))
789 fs(i,j,k2b)=cff*(z_r(i,j,kk+1)- &
790 & z_r(i,j,kk ))
791 END DO
792 END DO
793 END IF
794 END DO
795
796 IF (k.gt.0) THEN
797
798
799
800 DO j=jstr,jend
801 DO i=istr,iend
802#ifdef DIAGNOSTICS_TS
803
804#endif
805
806
807 ad_cff=ad_cff-ad_t(i,j,k,nnew,itrc)
808
809
810
811
812
814 adfac1=adfac*pm(i,j)*pn(i,j)
815 ad_fs(i,j,k1)=ad_fs(i,j,k1)-adfac
816 ad_fs(i,j,k2)=ad_fs(i,j,k2)+adfac
817 ad_fx(i ,j)=ad_fx(i ,j)-adfac1
818 ad_fx(i+1,j)=ad_fx(i+1,j)+adfac1
819 ad_fe(i,j )=ad_fe(i,j )-adfac1
820 ad_fe(i,j+1)=ad_fe(i,j+1)+adfac1
821 ad_cff=0.0_r8
822 END DO
823 END DO
825 DO j=jstr,jend
826 DO i=istr,iend
827#ifdef DIFF_3DCOEF
828# ifdef TS_U3ADV_SPLIT
829 difx=0.125_r8*(diff3d_u(i,j,k )+diff3d_u(i+1,j,k )+ &
830 & diff3d_u(i,j,k+1)+diff3d_u(i+1,j,k+1))
831 dife=0.125_r8*(diff3d_v(i,j,k )+diff3d_v(i,j+1,k )+ &
832 & diff3d_v(i,j,k+1)+diff3d_v(i,j+1,k+1))
833# else
834 difx=0.5_r8*diff3d_r(i,j,k)
835 dife=difx
836# endif
837#else
838 difx=0.5_r8*diff4(i,j,itrc)
839 dife=difx
840#endif
841 cff1=max(drdx(i ,j,k1),0.0_r8)
842 cff2=max(drdx(i+1,j,k2),0.0_r8)
843 cff3=min(drdx(i ,j,k2),0.0_r8)
844 cff4=min(drdx(i+1,j,k1),0.0_r8)
845 cff=difx* &
846 & (cff1*(cff1*dtdr(i,j,k2)-dtdx(i ,j,k1))+ &
847 & cff2*(cff2*dtdr(i,j,k2)-dtdx(i+1,j,k2))+ &
848 & cff3*(cff3*dtdr(i,j,k2)-dtdx(i ,j,k2))+ &
849 & cff4*(cff4*dtdr(i,j,k2)-dtdx(i+1,j,k1)))
850 cff1=max(drde(i,j ,k1),0.0_r8)
851 cff2=max(drde(i,j+1,k2),0.0_r8)
852 cff3=min(drde(i,j ,k2),0.0_r8)
853 cff4=min(drde(i,j+1,k1),0.0_r8)
854 cff=cff+ &
855 & dife* &
856 & (cff1*(cff1*dtdr(i,j,k2)-dtde(i,j ,k1))+ &
857 & cff2*(cff2*dtdr(i,j,k2)-dtde(i,j+1,k2))+ &
858 & cff3*(cff3*dtdr(i,j,k2)-dtde(i,j ,k2))+ &
859 & cff4*(cff4*dtdr(i,j,k2)-dtde(i,j+1,k1)))
860
861
862
863 ad_cff=ad_cff+fs(i,j,k2)*ad_fs(i,j,k2)
864 ad_fs(i,j,k2)=cff*ad_fs(i,j,k2)
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888 adfac=dife*ad_cff
889 ad_cff1=ad_cff1+ &
890 & (2.0_r8*cff1*dtdr(i,j,k2)-dtde(i,j ,k1))* &
891 & adfac
892 ad_cff2=ad_cff2+ &
893 & (2.0_r8*cff2*dtdr(i,j,k2)-dtde(i,j+1,k2))* &
894 & adfac
895 ad_cff3=ad_cff3+ &
896 & (2.0_r8*cff3*dtdr(i,j,k2)-dtde(i,j ,k2))* &
897 & adfac
898 ad_cff4=ad_cff4+ &
899 & (2.0_r8*cff4*dtdr(i,j,k2)-dtde(i,j+1,k1))* &
900 & adfac
901 ad_dtdr(i,j,k2)=ad_dtdr(i,j,k2)+ &
902 & (cff1*cff1+ &
903 & cff2*cff2+ &
904 & cff3*cff3+ &
905 & cff4*cff4)*adfac
906 ad_dtde(i,j ,k1)=ad_dtde(i,j ,k1)-cff1*adfac
907 ad_dtde(i,j+1,k2)=ad_dtde(i,j+1,k2)-cff2*adfac
908 ad_dtde(i,j ,k2)=ad_dtde(i,j ,k2)-cff3*adfac
909 ad_dtde(i,j+1,k1)=ad_dtde(i,j+1,k1)-cff4*adfac
910
911
912
913 ad_drde(i,j+1,k1)=ad_drde(i,j+1,k1)+ &
914 & (0.5_r8+sign(0.5_r8, &
915 & -drde(i,j+1,k1)))* &
916 & ad_cff4
917 ad_cff4=0.0_r8
918
919
920
921 ad_drde(i,j ,k2)=ad_drde(i,j ,k2)+ &
922 & (0.5_r8+sign(0.5_r8, &
923 & -drde(i,j ,k2)))* &
924 & ad_cff3
925 ad_cff3=0.0_r8
926
927
928
929 ad_drde(i,j+1,k2)=ad_drde(i,j+1,k2)+ &
930 & (0.5_r8+sign(0.5_r8, &
931 & drde(i,j+1,k2)))* &
932 & ad_cff2
933 ad_cff2=0.0_r8
934
935
936
937 ad_drde(i,j ,k1)=ad_drde(i,j ,k1)+ &
938 & (0.5_r8+sign(0.5_r8, &
939 & drde(i,j ,k1)))* &
940 & ad_cff1
941 ad_cff1=0.0_r8
942
943 cff1=max(drdx(i ,j,k1),0.0_r8)
944 cff2=max(drdx(i+1,j,k2),0.0_r8)
945 cff3=min(drdx(i ,j,k2),0.0_r8)
946 cff4=min(drdx(i+1,j,k1),0.0_r8)
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969 adfac=difx*ad_cff
970 ad_cff1=ad_cff1+ &
971 & (2.0_r8*cff1*dtdr(i,j,k2)-dtdx(i ,j,k1))* &
972 & adfac
973 ad_cff2=ad_cff2+ &
974 & (2.0_r8*cff2*dtdr(i,j,k2)-dtdx(i+1,j,k2))* &
975 & adfac
976 ad_cff3=ad_cff3+ &
977 & (2.0_r8*cff3*dtdr(i,j,k2)-dtdx(i ,j,k2))* &
978 & adfac
979 ad_cff4=ad_cff4+ &
980 & (2.0_r8*cff4*dtdr(i,j,k2)-dtdx(i+1,j,k1))* &
981 & adfac
982 ad_dtdr(i,j,k2)=ad_dtdr(i,j,k2)+ &
983 & (cff1*cff1+ &
984 & cff2*cff2+ &
985 & cff3*cff3+ &
986 & cff4*cff4)*adfac
987 ad_dtdx(i ,j,k1)=ad_dtdx(i ,j,k1)-cff1*adfac
988 ad_dtdx(i+1,j,k2)=ad_dtdx(i+1,j,k2)-cff2*adfac
989 ad_dtdx(i ,j,k2)=ad_dtdx(i ,j,k2)-cff3*adfac
990 ad_dtdx(i+1,j,k1)=ad_dtdx(i+1,j,k1)-cff4*adfac
991 ad_cff=0.0_r8
992
993
994
995 ad_drdx(i+1,j,k1)=ad_drdx(i+1,j,k1)+ &
996 & (0.5_r8+sign(0.5_r8, &
997 & -drdx(i+1,j,k1)))* &
998 & ad_cff4
999 ad_cff4=0.0_r8
1000
1001
1002
1003 ad_drdx(i ,j,k2)=ad_drdx(i ,j,k2)+ &
1004 & (0.5_r8+sign(0.5_r8, &
1005 & -drdx(i ,j,k2)))* &
1006 & ad_cff3
1007 ad_cff3=0.0_r8
1008
1009
1010
1011 ad_drdx(i+1,j,k2)=ad_drdx(i+1,j,k2)+ &
1012 & (0.5_r8+sign(0.5_r8, &
1013 & drdx(i+1,j,k2)))* &
1014 & ad_cff2
1015 ad_cff2=0.0_r8
1016
1017
1018
1019 ad_drdx(i ,j,k1)=ad_drdx(i ,j,k1)+ &
1020 & (0.5_r8+sign(0.5_r8, &
1021 & drdx(i ,j,k1)))* &
1022 & ad_cff1
1023 ad_cff1=0.0_r8
1024 END DO
1025 END DO
1026 END IF
1027 DO j=jstr,jend+1
1028 DO i=istr,iend
1029#ifdef DIFF_3DCOEF
1030# ifdef TS_U3ADV_SPLIT
1031 cff=0.5_r8*diff3d_v(i,j,k)*om_v(i,j)
1032# else
1033 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
1034 & om_v(i,j)
1035# endif
1036#else
1037 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
1038 & om_v(i,j)
1039#endif
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066 adfac=cff*ad_fe(i,j)
1067 adfac1=adfac*(dtde(i,j,k1)- &
1068 & 0.5_r8*(max(drde(i,j,k1),0.0_r8)* &
1069 & (dtdr(i,j-1,k1)+ &
1070 & dtdr(i,j ,k2))+ &
1071 & min(drde(i,j,k1),0.0_r8)* &
1072 & (dtdr(i,j-1,k2)+ &
1073 & dtdr(i,j ,k1))))
1074 adfac2=adfac*(hz(i,j,k)+hz(i,j-1,k))
1075 adfac3=adfac2*0.5_r8*max(drde(i,j,k1),0.0_r8)
1076 adfac4=adfac2*0.5_r8*min(drde(i,j,k1),0.0_r8)
1077 ad_hz(i,j-1,k)=ad_hz(i,j-1,k)+adfac1
1078 ad_hz(i,j ,k)=ad_hz(i,j ,k)+adfac1
1079 ad_dtde(i,j,k1)=ad_dtde(i,j,k1)+adfac2
1080 ad_dtdr(i,j-1,k1)=ad_dtdr(i,j-1,k1)-adfac3
1081 ad_dtdr(i,j ,k2)=ad_dtdr(i,j ,k2)-adfac3
1082 ad_dtdr(i,j-1,k2)=ad_dtdr(i,j-1,k2)-adfac4
1083 ad_dtdr(i,j ,k1)=ad_dtdr(i,j ,k1)-adfac4
1084 ad_drde(i,j,k1)=ad_drde(i,j,k1)- &
1085 & adfac2*0.5_r8* &
1086 & ((0.5_r8+sign(0.5_r8, drde(i,j,k1)))* &
1087 & (dtdr(i,j-1,k1)+dtdr(i,j,k2))+ &
1088 & (0.5_r8+sign(0.5_r8,-drde(i,j,k1)))* &
1089 & (dtdr(i,j-1,k2)+dtdr(i,j,k1)))
1090 ad_fe(i,j)=0.0_r8
1091 END DO
1092 END DO
1093 DO j=jstr,jend
1094 DO i=istr,iend+1
1095#ifdef DIFF_3DCOEF
1096# ifdef TS_U3ADV_SPLIT
1097 cff=0.5_r8*diff3d_u(i,j,k)*on_u(i,j)
1098# else
1099 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
1100 & on_u(i,j)
1101# endif
1102#else
1103 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
1104 & on_u(i,j)
1105#endif
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132 adfac=cff*ad_fx(i,j)
1133 adfac1=adfac*(dtdx(i ,j,k1)- &
1134 & 0.5_r8*(max(drdx(i,j,k1),0.0_r8)* &
1135 & (dtdr(i-1,j,k1)+ &
1136 & dtdr(i ,j,k2))+ &
1137 & min(drdx(i,j,k1),0.0_r8)* &
1138 & (dtdr(i-1,j,k2)+ &
1139 & dtdr(i ,j,k1))))
1140 adfac2=adfac*(hz(i,j,k)+hz(i-1,j,k))
1141 adfac3=adfac2*0.5_r8*max(drdx(i,j,k1),0.0_r8)
1142 adfac4=adfac2*0.5_r8*min(drdx(i,j,k1),0.0_r8)
1143 ad_hz(i-1,j,k)=ad_hz(i-1,j,k)+adfac1
1144 ad_hz(i ,j,k)=ad_hz(i ,j,k)+adfac1
1145 ad_dtdx(i ,j,k1)=ad_dtdx(i ,j,k1)+adfac2
1146 ad_dtdr(i-1,j,k1)=ad_dtdr(i-1,j,k1)-adfac3
1147 ad_dtdr(i ,j,k2)=ad_dtdr(i ,j,k2)-adfac3
1148 ad_dtdr(i-1,j,k2)=ad_dtdr(i-1,j,k2)-adfac4
1149 ad_dtdr(i ,j,k1)=ad_dtdr(i ,j,k1)-adfac4
1150 ad_drdx(i,j,k1)=ad_drdx(i,j,k1)- &
1151 & adfac2*0.5_r8* &
1152 & ((0.5_r8+sign(0.5_r8, drdx(i,j,k1)))* &
1153 & (dtdr(i-1,j,k1)+dtdr(i,j,k2))+ &
1154 & (0.5_r8+sign(0.5_r8,-drdx(i,j,k1)))* &
1155 & (dtdr(i-1,j,k2)+dtdr(i,j,k1)))
1156 ad_fx(i,j)=0.0_r8
1157 END DO
1158 END DO
1159 END IF
1160 IF ((k.eq.0).or.(k.eq.
n(ng)))
THEN
1161 DO j=jstr-1,jend+1
1162 DO i=istr-1,iend+1
1163
1164
1165 ad_fs(i,j,k2)=0.0_r8
1166
1167
1168 ad_dtdr(i,j,k2)=0.0_r8
1169 END DO
1170 END DO
1171 ELSE
1172 DO j=jstr-1,jend+1
1173 DO i=istr-1,iend+1
1174#if defined TS_MIX_MAX_SLOPE
1175 cff1=sqrt(drdx(i,j,k2)**2+drdx(i+1,j,k2)**2+ &
1176 & drdx(i,j,k1)**2+drdx(i+1,j,k1)**2+ &
1177 & drde(i,j,k2)**2+drde(i,j+1,k2)**2+ &
1178 & drde(i,j,k1)**2+drde(i,j+1,k1)**2)
1179 cff2=0.25_r8*slope_max* &
1180 & (z_r(i,j,k+1)-z_r(i,j,k))*cff1
1181 cff3=max(pden(i,j,k)-pden(i,j,k+1),small)
1182 cff4=max(cff2,cff3)
1183 cff=-1.0_r8/cff4
1184#elif defined TS_MIX_MIN_STRAT
1185 cff1=max(pden(i,j,k)-pden(i,j,k+1), &
1186 & strat_min*(z_r(i,j,k+1)-z_r(i,j,k)))
1187 cff=-1.0_r8/cff1
1188#else
1189 cff1=max(pden(i,j,k)-pden(i,j,k+1),eps)
1190 cff=-1.0_r8/cff1
1191#endif
1192
1193
1194
1195
1196
1197 adfac=cff*ad_fs(i,j,k2)
1198 ad_z_r(i,j,k )=ad_z_r(i,j,k )-adfac
1199 ad_z_r(i,j,k+1)=ad_z_r(i,j,k+1)+adfac
1200 ad_cff=ad_cff+(z_r(i,j,k+1)- &
1201 & z_r(i,j,k ))*ad_fs(i,j,k2)
1202 ad_fs(i,j,k2)=0.0_r8
1203
1204
1205
1206
1207
1208 adfac=cff*ad_dtdr(i,j,k2)
1209 ad_lapt(i,j,k )=ad_lapt(i,j,k )-adfac
1210 ad_lapt(i,j,k+1)=ad_lapt(i,j,k+1)+adfac
1211 ad_cff=ad_cff+(lapt(i,j,k+1)- &
1212 & lapt(i,j,k ))*ad_dtdr(i,j,k2)
1213 ad_dtdr(i,j,k2)=0.0_r8
1214#if defined TS_MIX_MAX_SLOPE
1215
1216
1217 ad_cff4=ad_cff4+cff*cff*ad_cff
1218 ad_cff=0.0_r8
1219
1220
1221
1222 ad_cff3=ad_cff3+ &
1223 & (0.5_r8-sign(0.5_r8,cff2-cff3))*ad_cff4
1224 ad_cff2=ad_cff2+ &
1225 & (0.5_r8+sign(0.5_r8,cff2-cff3))*ad_cff4
1226 ad_cff4=0.0_r8
1227
1228
1229
1230
1231 adfac=(0.5_r8+sign(0.5_r8,pden(i,j,k)-pden(i,j,k+1)- &
1232 & small))*ad_cff3
1233 ad_pden(i,j,k )=ad_pden(i,j,k )+adfac
1234 ad_pden(i,j,k+1)=ad_pden(i,j,k+1)-adfac
1235 ad_cff3=0.0_r8
1236
1237
1238
1239
1240 adfac=0.25_r8*slope_max*ad_cff2
1241 adfac1=adfac*cff1
1242 ad_cff1=ad_cff1+(z_r(i,j,k+1)-z_r(i,j,k))*adfac
1243 ad_z_r(i,j,k )=ad_z_r(i,j,k )-adfac1
1244 ad_z_r(i,j,k+1)=ad_z_r(i,j,k+1)+adfac1
1245 ad_cff2=0.0_r8
1246 IF (cff1.ne.0.0_r8) THEN
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256 adfac=ad_cff1/cff1
1257 ad_drdx(i ,j,k1)=ad_drdx(i ,j,k1)+ &
1258 & drdx(i ,j,k1)*adfac
1259 ad_drdx(i+1,j,k1)=ad_drdx(i+1,j,k1)+ &
1260 & drdx(i+1,j,k1)*adfac
1261 ad_drdx(i ,j,k2)=ad_drdx(i ,j,k2)+ &
1262 & drdx(i ,j,k2)*adfac
1263 ad_drdx(i+1,j,k2)=ad_drdx(i+1,j,k2)+ &
1264 & drdx(i+1,j,k2)*adfac
1265 ad_drde(i,j ,k2)=ad_drde(i,j ,k2)+ &
1266 & drde(i,j ,k2)*adfac
1267 ad_drde(i,j+1,k2)=ad_drde(i,j+1,k2)+ &
1268 & drde(i,j+1,k2)*adfac
1269 ad_drde(i,j ,k1)=ad_drde(i,j ,k1)+ &
1270 & drde(i,j ,k1)*adfac
1271 ad_drde(i,j+1,k1)=ad_drde(i,j+1,k1)+ &
1272 & drde(i,j+1,k1)*adfac
1273 ad_cff1=0.0_r8
1274 ELSE
1275
1276
1277 ad_cff1=0.0_r8
1278 END IF
1279#elif defined TS_MIX_MIN_STRAT
1280
1281
1282 ad_cff1=ad_cff1+cff*cff*ad_cff
1283 ad_cff=0.0_r8
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295 adfac1=(0.5_r8+sign(0.5_r8, &
1296 & pden(i,j,k)-pden(i,j,k+1)- &
1297 & strat_min*(z_r(i,j,k+1)- &
1298 & z_r(i,j,k ))))* &
1299 & ad_cff1
1300 adfac2=(0.5_r8-sign(0.5_r8, &
1301 & pden(i,j,k)-pden(i,j,k+1)- &
1302 & strat_min*(z_r(i,j,k+1)- &
1303 & z_r(i,j,k ))))* &
1304 & strat_min*ad_cff1
1305 ad_pden(i,j,k )=ad_pden(i,j,k )+adfac1
1306 ad_pden(i,j,k+1)=ad_pden(i,j,k+1)-adfac1
1307 ad_z_r(i,j,k )=ad_z_r(i,j,k )-adfac2
1308 ad_z_r(i,j,k+1)=ad_z_r(i,j,k+1)+adfac2
1309 ad_cff1=0.0_r8
1310#else
1311
1312
1313 ad_cff1=ad_cff1+cff*cff*ad_cff
1314 ad_cff=0.0_r8
1315
1316
1317
1318
1319 adfac=(0.5_r8+sign(0.5_r8, &
1320 & pden(i,j,k)-pden(i,j,k+1)-eps))* &
1321 & ad_cff1
1322 ad_pden(i,j,k )=ad_pden(i,j,k )+adfac
1323 ad_pden(i,j,k+1)=ad_pden(i,j,k+1)-adfac
1324 ad_cff1=0.0_r8
1325#endif
1326 END DO
1327 END DO
1328 END IF
1329 IF (k.lt.
n(ng))
THEN
1330 DO j=jstr,jend+1
1331 DO i=istr,iend
1332 cff=0.5_r8*(pn(i,j)+pn(i,j-1))
1333#ifdef MASKING
1334 cff=cff*vmask(i,j)
1335#endif
1336#ifdef WET_DRY_NOT_YET
1337 cff=cff*vmask_wet(i,j)
1338#endif
1339
1340
1341
1342 adfac=cff*ad_dtde(i,j,k2)
1343 ad_lapt(i,j-1,k+1)=ad_lapt(i,j-1,k+1)-adfac
1344 ad_lapt(i,j ,k+1)=ad_lapt(i,j ,k+1)+adfac
1345 ad_dtde(i,j,k2)=0.0_r8
1346
1347
1348
1349 adfac=cff*ad_drde(i,j,k2)
1350 ad_pden(i,j-1,k+1)=ad_pden(i,j-1,k+1)-adfac
1351 ad_pden(i,j ,k+1)=ad_pden(i,j ,k+1)+adfac
1352 ad_drde(i,j,k2)=0.0_r8
1353 END DO
1354 END DO
1355 DO j=jstr,jend
1356 DO i=istr,iend+1
1357 cff=0.5_r8*(pm(i,j)+pm(i-1,j))
1358#ifdef MASKING
1359 cff=cff*umask(i,j)
1360#endif
1361#ifdef WET_DRY_NOT_YET
1362 cff=cff*umask_wet(i,j)
1363#endif
1364
1365
1366
1367 adfac=cff*ad_dtdx(i,j,k2)
1368 ad_lapt(i-1,j,k+1)=ad_lapt(i-1,j,k+1)-adfac
1369 ad_lapt(i ,j,k+1)=ad_lapt(i ,j,k+1)+adfac
1370 ad_dtdx(i,j,k2)=0.0_r8
1371
1372
1373
1374 adfac=cff*ad_drdx(i,j,k2)
1375 ad_pden(i-1,j,k+1)=ad_pden(i-1,j,k+1)-adfac
1376 ad_pden(i ,j,k+1)=ad_pden(i ,j,k+1)+adfac
1377 ad_drdx(i,j,k2)=0.0_r8
1378 END DO
1379 END DO
1380 END IF
1381
1382
1383
1384 kt=k2
1385 k2=k1
1386 k1=kt
1387 END DO k_loop2
1388
1389
1390
1391
1394 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
1396
1397
1398
1399
1400 adfac=0.5_r8*ad_lapt(iend+1,jend+1,k)
1401 ad_lapt(iend+1,jend ,k)=ad_lapt(iend+1,jend ,k)+adfac
1402 ad_lapt(iend ,jend+1,k)=ad_lapt(iend ,jend+1,k)+adfac
1403 ad_lapt(iend+1,jend+1,k)=0.0_r8
1404 END DO
1405 END IF
1406 END IF
1407
1410 IF (
domain(ng)%NorthWest_Corner(tile))
THEN
1412
1413
1414
1415
1416 adfac=0.5_r8*ad_lapt(istr-1,jend+1,k)
1417 ad_lapt(istr-1,jend ,k)=ad_lapt(istr-1,jend ,k)+adfac
1418 ad_lapt(istr ,jend+1,k)=ad_lapt(istr ,jend+1,k)+adfac
1419 ad_lapt(istr-1,jend+1,k)=0.0_r8
1420 END DO
1421 END IF
1422 END IF
1423
1426 IF (
domain(ng)%SouthEast_Corner(tile))
THEN
1428
1429
1430
1431
1432 adfac=0.5_r8*ad_lapt(iend+1,jstr-1,k)
1433 ad_lapt(iend ,jstr-1,k)=ad_lapt(iend ,jstr-1,k)+adfac
1434 ad_lapt(iend+1,jstr ,k)=ad_lapt(iend+1,jstr ,k)+adfac
1435 ad_lapt(iend+1,jstr-1,k)=0.0_r8
1436 END DO
1437 END IF
1438 END IF
1439
1442 IF (
domain(ng)%SouthWest_Corner(tile))
THEN
1444
1445
1446
1447
1448 adfac=0.5_r8*ad_lapt(istr-1,jstr-1,k)
1449 ad_lapt(istr ,jstr-1,k)=ad_lapt(istr ,jstr-1,k)+adfac
1450 ad_lapt(istr-1,jstr ,k)=ad_lapt(istr-1,jstr ,k)+adfac
1451 ad_lapt(istr-1,jstr-1,k)=0.0_r8
1452 END DO
1453 END IF
1454 END IF
1455
1457 IF (
domain(ng)%Northern_Edge(tile))
THEN
1460 DO i=imin,imax
1461
1462
1463 ad_lapt(i,jend+1,k)=0.0_r8
1464 END DO
1465 END DO
1466 ELSE
1468 DO i=imin,imax
1469
1470
1471 ad_lapt(i,jend,k)=ad_lapt(i,jend,k)+ &
1472 & ad_lapt(i,jend+1,k)
1473 ad_lapt(i,jend+1,k)=0.0_r8
1474 END DO
1475 END DO
1476 END IF
1477 END IF
1478 END IF
1479
1481 IF (
domain(ng)%Southern_Edge(tile))
THEN
1484 DO i=imin,imax
1485
1486
1487 ad_lapt(i,jstr-1,k)=0.0_r8
1488 END DO
1489 END DO
1490 ELSE
1492 DO i=imin,imax
1493
1494
1495 ad_lapt(i,jstr,k)=ad_lapt(i,jstr,k)+ &
1496 & ad_lapt(i,jstr-1,k)
1497 ad_lapt(i,jstr-1,k)=0.0_r8
1498 END DO
1499 END DO
1500 END IF
1501 END IF
1502 END IF
1503
1505 IF (
domain(ng)%Eastern_Edge(tile))
THEN
1508 DO j=jmin,jmax
1509
1510
1511 ad_lapt(iend+1,j,k)=0.0_r8
1512 END DO
1513 END DO
1514 ELSE
1516 DO j=jmin,jmax
1517
1518
1519 ad_lapt(iend,j,k)=ad_lapt(iend,j,k)+ &
1520 & ad_lapt(iend+1,j,k)
1521 ad_lapt(iend+1,j,k)=0.0_r8
1522 END DO
1523 END DO
1524 END IF
1525 END IF
1526 END IF
1527
1529 IF (
domain(ng)%Western_Edge(tile))
THEN
1532 DO j=jmin,jmax
1533
1534
1535 ad_lapt(istr-1,j,k)=0.0_r8
1536 END DO
1537 END DO
1538 ELSE
1540 DO j=jmin,jmax
1541
1542
1543 ad_lapt(istr,j,k)=ad_lapt(istr,j,k)+ &
1544 & ad_lapt(istr-1,j,k)
1545 ad_lapt(istr-1,j,k)=0.0_r8
1546 END DO
1547 END DO
1548 END IF
1549 END IF
1550 END IF
1551
1552
1553
1554
1555
1556
1557
1558
1559 k1=2
1560 k2=1
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571 k1=k2
1572 k2=3-k1
1573 END DO
1574
1575
1576
1577
1578 k_loop3:
DO k=
n(ng),0,-1
1579 k2b=1
1580 DO kk=0,k
1581 k1b=k2b
1582 k2b=3-k1b
1583 IF (kk.lt.
n(ng))
THEN
1584 DO j=jmin,jmax
1585 DO i=imin,imax+1
1586 cff=0.5_r8*(pm(i,j)+pm(i-1,j))
1587#ifdef MASKING
1588 cff=cff*umask(i,j)
1589#endif
1590#ifdef WET_DRY_NOT_YET
1591 cff=cff*umask_wet(i,j)
1592#endif
1593 drdx(i,j,k2b)=cff*(pden(i ,j,kk+1)- &
1594 & pden(i-1,j,kk+1))
1595#if defined TS_MIX_STABILITY
1596 dtdx(i,j,k2b)=cff*(0.75_r8*(t(i ,j,kk+1,nrhs,itrc)- &
1597 & t(i-1,j,kk+1,nrhs,itrc))+ &
1598 & 0.25_r8*(t(i ,j,kk+1,nstp,itrc)- &
1599 & t(i-1,j,kk+1,nstp,itrc)))
1600#elif defined TS_MIX_CLIMA
1602 dtdx(i,j,k2b)=cff*((t(i ,j,kk+1,nrhs,itrc)- &
1603 & tclm(i ,j,kk+1,itrc))- &
1604 & (t(i-1,j,kk+1,nrhs,itrc)- &
1605 & tclm(i-1,j,kk+1,itrc)))
1606 ELSE
1607 dtdx(i,j,k2b)=cff*(t(i ,j,kk+1,nrhs,itrc)- &
1608 & t(i-1,j,kk+1,nrhs,itrc))
1609 END IF
1610#else
1611 dtdx(i,j,k2b)=cff*(t(i ,j,kk+1,nrhs,itrc)- &
1612 & t(i-1,j,kk+1,nrhs,itrc))
1613#endif
1614 END DO
1615 END DO
1616 IF (kk.eq.0) THEN
1617 DO j=jmin,jmax
1618 DO i=imin,imax+1
1619 drdx(i,j,k1b)=0.0_r8
1620 dtdx(i,j,k1b)=0.0_r8
1621 END DO
1622 END DO
1623 END IF
1624 DO j=jmin,jmax+1
1625 DO i=imin,imax
1626 cff=0.5_r8*(pn(i,j)+pn(i,j-1))
1627#ifdef MASKING
1628 cff=cff*vmask(i,j)
1629#endif
1630#ifdef WET_DRY_NOT_YET
1631 cff=cff*vmask_wet(i,j)
1632#endif
1633 drde(i,j,k2b)=cff*(pden(i,j ,kk+1)- &
1634 & pden(i,j-1,kk+1))
1635#if defined TS_MIX_STABILITY
1636 dtde(i,j,k2b)=cff*(0.75_r8*(t(i,j ,kk+1,nrhs,itrc)- &
1637 & t(i,j-1,kk+1,nrhs,itrc))+ &
1638 & 0.25_r8*(t(i,j ,kk+1,nstp,itrc)- &
1639 & t(i,j-1,kk+1,nstp,itrc)))
1640#elif defined TS_MIX_CLIMA
1642 dtde(i,j,k2b)=cff*((t(i,j ,kk+1,nrhs,itrc)- &
1643 & tclm(i,j ,kk+1,itrc))- &
1644 & (t(i,j-1,kk+1,nrhs,itrc)- &
1645 & tclm(i,j-1,kk+1,itrc)))
1646 ELSE
1647 dtde(i,j,k2b)=cff*(t(i,j ,kk+1,nrhs,itrc)- &
1648 & t(i,j-1,kk+1,nrhs,itrc))
1649 END IF
1650#else
1651 dtde(i,j,k2b)=cff*(t(i,j ,kk+1,nrhs,itrc)- &
1652 & t(i,j-1,kk+1,nrhs,itrc))
1653#endif
1654 END DO
1655 END DO
1656 IF (kk.eq.0) THEN
1657 DO j=jmin,jmax+1
1658 DO i=imin,imax
1659 drde(i,j,k1b)=0.0_r8
1660 dtde(i,j,k1b)=0.0_r8
1661 END DO
1662 END DO
1663 END IF
1664 END IF
1665 IF ((kk.eq.0).or.(kk.eq.
n(ng)))
THEN
1666 DO j=-1+jmin,jmax+1
1667 DO i=-1+imin,imax+1
1668 dtdr(i,j,k2b)=0.0_r8
1669 fs(i,j,k2b)=0.0_r8
1670 END DO
1671 END DO
1672 IF (kk.eq.0) THEN
1673 DO j=-1+jmin,jmax+1
1674 DO i=-1+imin,imax+1
1675 dtdr(i,j,k1b)=0.0_r8
1676 fs(i,j,k1b)=0.0_r8
1677 END DO
1678 END DO
1679 END IF
1680 ELSE
1681 DO j=-1+jmin,jmax+1
1682 DO i=-1+imin,imax+1
1683#if defined TS_MIX_MAX_SLOPE
1684 cff1=sqrt(drdx(i,j,k2b)**2+drdx(i+1,j,k2b)**2+ &
1685 & drdx(i,j,k1b)**2+drdx(i+1,j,k1b)**2+ &
1686 & drde(i,j,k2b)**2+drde(i,j+1,k2b)**2+ &
1687 & drde(i,j,k1b)**2+drde(i,j+1,k1b)**2)
1688 cff2=0.25_r8*slope_max* &
1689 & (z_r(i,j,kk+1)-z_r(i,j,kk))*cff1
1690 cff3=max(pden(i,j,kk)-pden(i,j,kk+1),small)
1691 cff4=max(cff2,cff3)
1692 cff=-1.0_r8/cff4
1693#elif defined TS_MIX_MIN_STRAT
1694 cff1=max(pden(i,j,kk)-pden(i,j,kk+1), &
1695 & strat_min*(z_r(i,j,kk+1)-z_r(i,j,kk)))
1696 cff=-1.0_r8/cff1
1697#else
1698 cff1=max(pden(i,j,kk)-pden(i,j,kk+1),eps)
1699 cff=-1.0_r8/cff1
1700#endif
1701#if defined TS_MIX_STABILITY
1702 dtdr(i,j,k2b)=cff*(0.75_r8*(t(i,j,kk+1,nrhs,itrc)- &
1703 & t(i,j,kk ,nrhs,itrc))+ &
1704 & 0.25_r8*(t(i,j,kk+1,nstp,itrc)- &
1705 & t(i,j,kk ,nstp,itrc)))
1706#elif defined TS_MIX_CLIMA
1708 dtdr(i,j,k2b)=cff*((t(i,j,kk+1,nrhs,itrc)- &
1709 & tclm(i,j,kk+1,itrc))- &
1710 & (t(i,j,kk ,nrhs,itrc)- &
1711 & tclm(i,j,kk ,itrc)))
1712 ELSE
1713 dtdr(i,j,k2b)=cff*(t(i,j,kk+1,nrhs,itrc)- &
1714 & t(i,j,kk ,nrhs,itrc))
1715 END IF
1716#else
1717 dtdr(i,j,k2b)=cff*(t(i,j,kk+1,nrhs,itrc)- &
1718 & t(i,j,kk ,nrhs,itrc))
1719#endif
1720 fs(i,j,k2b)=cff*(z_r(i,j,kk+1)- &
1721 & z_r(i,j,kk ))
1722 END DO
1723 END DO
1724 END IF
1725 IF (kk.gt.0) THEN
1726 DO j=jmin,jmax
1727 DO i=imin,imax+1
1728#ifdef DIFF_3DCOEF
1729# ifdef TS_U3ADV_SPLIT
1730 cff=0.5_r8*diff3d_u(i,j,k)*on_u(i,j)
1731# else
1732 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
1733 & on_u(i,j)
1734# endif
1735#else
1736 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
1737 & on_u(i,j)
1738#endif
1739 fx(i,j)=cff* &
1740 & (hz(i,j,kk)+hz(i-1,j,kk))* &
1741 & (dtdx(i ,j,k1b)- &
1742 & 0.5_r8*(max(drdx(i,j,k1b),0.0_r8)* &
1743 & (dtdr(i-1,j,k1b)+ &
1744 & dtdr(i ,j,k2b))+ &
1745 & min(drdx(i,j,k1b),0.0_r8)* &
1746 & (dtdr(i-1,j,k2b)+ &
1747 & dtdr(i ,j,k1b))))
1748 END DO
1749 END DO
1750 DO j=jmin,jmax+1
1751 DO i=imin,imax
1752#ifdef DIFF_3DCOEF
1753# ifdef TS_U3ADV_SPLIT
1754 cff=0.5_r8*diff3d_v(i,j,k)*om_v(i,j)
1755# else
1756 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
1757 & om_v(i,j)
1758# endif
1759#else
1760 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
1761 & om_v(i,j)
1762#endif
1763 fe(i,j)=cff* &
1764 & (hz(i,j,kk)+hz(i,j-1,kk))* &
1765 & (dtde(i,j,k1b)- &
1766 & 0.5_r8*(max(drde(i,j,k1b),0.0_r8)* &
1767 & (dtdr(i,j-1,k1b)+ &
1768 & dtdr(i,j ,k2b))+ &
1769 & min(drde(i,j,k1b),0.0_r8)* &
1770 & (dtdr(i,j-1,k2b)+ &
1771 & dtdr(i,j ,k1b))))
1772 END DO
1773 END DO
1774 IF (kk.lt.
n(ng))
THEN
1775 DO j=jmin,jmax
1776 DO i=imin,imax
1777#ifdef DIFF_3DCOEF
1778# ifdef TS_U3ADV_SPLIT
1779 difx=0.125_r8* &
1780 & (diff3d_u(i,j,kk )+diff3d_u(i+1,j,kk )+ &
1781 & diff3d_u(i,j,kk+1)+diff3d_u(i+1,j,kk+1))
1782 dife=0.125_r8* &
1783 & (diff3d_v(i,j,kk )+diff3d_v(i,j+1,kk )+ &
1784 & diff3d_v(i,j,kk+1)+diff3d_v(i,j+1,kk+1))
1785# else
1786 difx=0.5_r8*diff3d_r(i,j,kk)
1787 dife=difx
1788# endif
1789#else
1790 difx=0.5_r8*diff4(i,j,itrc)
1791 dife=difx
1792#endif
1793 cff1=max(drdx(i ,j,k1b),0.0_r8)
1794 cff2=max(drdx(i+1,j,k2b),0.0_r8)
1795 cff3=min(drdx(i ,j,k2b),0.0_r8)
1796 cff4=min(drdx(i+1,j,k1b),0.0_r8)
1797 cff=difx* &
1798 & (cff1*(cff1*dtdr(i,j,k2b)-dtdx(i ,j,k1b))+ &
1799 & cff2*(cff2*dtdr(i,j,k2b)-dtdx(i+1,j,k2b))+ &
1800 & cff3*(cff3*dtdr(i,j,k2b)-dtdx(i ,j,k2b))+ &
1801 & cff4*(cff4*dtdr(i,j,k2b)-dtdx(i+1,j,k1b)))
1802 cff1=max(drde(i,j ,k1b),0.0_r8)
1803 cff2=max(drde(i,j+1,k2b),0.0_r8)
1804 cff3=min(drde(i,j ,k2b),0.0_r8)
1805 cff4=min(drde(i,j+1,k1b),0.0_r8)
1806 cff=cff+ &
1807 & dife* &
1808 & (cff1*(cff1*dtdr(i,j,k2b)-dtde(i,j ,k1b))+ &
1809 & cff2*(cff2*dtdr(i,j,k2b)-dtde(i,j+1,k2b))+ &
1810 & cff3*(cff3*dtdr(i,j,k2b)-dtde(i,j ,k2b))+ &
1811 & cff4*(cff4*dtdr(i,j,k2b)-dtde(i,j+1,k1b)))
1812 fs1(i,j,k2b)=fs(i,j,k2b)
1813 fs(i,j,k2b)=cff*fs(i,j,k2b)
1814 END DO
1815 END DO
1816 IF (kk.eq.0) THEN
1817 DO j=jmin,jmax
1818 DO i=imin,imax
1819 fs1(i,j,k1b)=0.0_r8
1820 fs(i,j,k1b)=0.0_r8
1821 END DO
1822 END DO
1823 END IF
1824 END IF
1825 END IF
1826 END DO
1827
1828
1829
1830
1831
1832 IF (k.gt.0) THEN
1833 DO j=jmin,jmax
1834 DO i=imin,imax
1835 cff=pm(i,j)*pn(i,j)
1836 cff1=1.0_r8/hz(i,j,k)
1837
1838
1839
1840
1841
1842
1843
1844
1845
1846 adfac=cff1*ad_lapt(i,j,k)
1847 adfac1=adfac*cff
1848 ad_fs(i,j,k1)=ad_fs(i,j,k1)-adfac
1849 ad_fs(i,j,k2)=ad_fs(i,j,k2)+adfac
1850 ad_fe(i,j )=ad_fe(i,j )-adfac1
1851 ad_fe(i,j+1)=ad_fe(i,j+1)+adfac1
1852 ad_fx(i ,j)=ad_fx(i ,j)-adfac1
1853 ad_fx(i+1,j)=ad_fx(i+1,j)+adfac1
1854 ad_cff1=ad_cff1+(cff* &
1855 & (fx(i+1,j)-fx(i,j)+ &
1856 & fe(i,j+1)-fe(i,j))+ &
1857 & (fs(i,j,k2)-fs(i,j,k1)))* &
1858 & ad_lapt(i,j,k)
1859 ad_lapt(i,j,k)=0.0_r8
1860
1861
1862 ad_hz(i,j,k)=ad_hz(i,j,k)-cff1*cff1*ad_cff1
1863 ad_cff1=0.0_r8
1864 END DO
1865 END DO
1866 IF (k.lt.
n(ng))
THEN
1867 DO j=jmin,jmax
1868 DO i=imin,imax
1869#ifdef DIFF_3DCOEF
1870# ifdef TS_U3ADV_SPLIT
1871 difx=0.125_r8*(diff3d_u(i,j,k )+diff3d_u(i+1,j,k )+ &
1872 & diff3d_u(i,j,k+1)+diff3d_u(i+1,j,k+1))
1873 dife=0.125_r8*(diff3d_v(i,j,k )+diff3d_v(i,j+1,k )+ &
1874 & diff3d_v(i,j,k+1)+diff3d_v(i,j+1,k+1))
1875# else
1876 difx=0.5_r8*diff3d_r(i,j,k)
1877 dife=difx
1878# endif
1879#else
1880 difx=0.5_r8*diff4(i,j,itrc)
1881 dife=difx
1882#endif
1883 cff1=max(drdx(i ,j,k1),0.0_r8)
1884 cff2=max(drdx(i+1,j,k2),0.0_r8)
1885 cff3=min(drdx(i ,j,k2),0.0_r8)
1886 cff4=min(drdx(i+1,j,k1),0.0_r8)
1887 cff=difx* &
1888 & (cff1*(cff1*dtdr(i,j,k2)-dtdx(i ,j,k1))+ &
1889 & cff2*(cff2*dtdr(i,j,k2)-dtdx(i+1,j,k2))+ &
1890 & cff3*(cff3*dtdr(i,j,k2)-dtdx(i ,j,k2))+ &
1891 & cff4*(cff4*dtdr(i,j,k2)-dtdx(i+1,j,k1)))
1892 cff1=max(drde(i,j ,k1),0.0_r8)
1893 cff2=max(drde(i,j+1,k2),0.0_r8)
1894 cff3=min(drde(i,j ,k2),0.0_r8)
1895 cff4=min(drde(i,j+1,k1),0.0_r8)
1896 cff=cff+ &
1897 & dife* &
1898 & (cff1*(cff1*dtdr(i,j,k2)-dtde(i,j ,k1))+ &
1899 & cff2*(cff2*dtdr(i,j,k2)-dtde(i,j+1,k2))+ &
1900 & cff3*(cff3*dtdr(i,j,k2)-dtde(i,j ,k2))+ &
1901 & cff4*(cff4*dtdr(i,j,k2)-dtde(i,j+1,k1)))
1902
1903
1904
1905
1906 ad_fs(i,j,k2)=cff*ad_fs(i,j,k2)
1907 ad_cff=ad_cff+fs1(i,j,k2)*ad_fs(i,j,k2)
1908
1909
1910
1911
1912
1913
1914
1915
1916
1917
1918
1919
1920
1921
1922
1923
1924
1925
1926
1927
1928
1929
1930
1931 adfac*dife*ad_cff
1932 ad_cff1=ad_cff1+ &
1933 & (2.0_r8*cff1*dtdr(i,j,k2)-dtde(i,j ,k1))* &
1934 & adfac
1935 ad_cff2=ad_cff2+ &
1936 & (2.0_r8*cff2*dtdr(i,j,k2)-dtde(i,j+1,k2))* &
1937 & adfac
1938 ad_cff3=ad_cff3+ &
1939 & (2.0_r8*cff3*dtdr(i,j,k2)-dtde(i,j ,k2))* &
1940 & adfac
1941 ad_cff4=ad_cff4+ &
1942 & (2.0_r8*cff4*dtdr(i,j,k2)-dtde(i,j+1,k1))* &
1943 & adfac
1944 ad_dtdr(i,j,k2)=ad_dtdr(i,j,k2)+ &
1945 & (cff1*cff1+ &
1946 & cff2*cff2+ &
1947 & cff3*cff3+ &
1948 & cff4*cff4)*adfac
1949 ad_dtde(i,j ,k1)=ad_dtde(i,j ,k1)-cff1*adfac
1950 ad_dtde(i,j+1,k2)=ad_dtde(i,j+1,k2)-cff2*adfac
1951 ad_dtde(i,j ,k2)=ad_dtde(i,j ,k2)-cff3*adfac
1952 ad_dtde(i,j+1,k1)=ad_dtde(i,j+1,k1)-cff4*adfac
1953
1954
1955
1956 ad_drde(i,j+1,k1)=ad_drde(i,j+1,k1)+ &
1957 & (0.5_r8+sign(0.5_r8, &
1958 & -drde(i,j+1,k1)))* &
1959 & ad_cff4
1960 ad_cff4=0.0_r8
1961
1962
1963
1964 ad_drde(i,j ,k2)=ad_drde(i,j ,k2)+ &
1965 & (0.5_r8+sign(0.5_r8, &
1966 & -drde(i,j ,k2)))* &
1967 & ad_cff3
1968 ad_cff3=0.0_r8
1969
1970
1971
1972 ad_drde(i,j+1,k2)=ad_drde(i,j+1,k2)+ &
1973 & (0.5_r8+sign(0.5_r8, &
1974 & drde(i,j+1,k2)))* &
1975 & ad_cff2
1976 ad_cff2=0.0_r8
1977
1978
1979
1980 ad_drde(i ,j,k1)=ad_drde(i ,j,k1)+ &
1981 & (0.5_r8+sign(0.5_r8, &
1982 & drde(i ,j,k1)))* &
1983 & ad_cff1
1984 ad_cff1=0.0_r8
1985 cff1=max(drdx(i ,j,k1),0.0_r8)
1986 cff2=max(drdx(i+1,j,k2),0.0_r8)
1987 cff3=min(drdx(i ,j,k2),0.0_r8)
1988 cff4=min(drdx(i+1,j,k1),0.0_r8)
1989
1990
1991
1992
1993
1994
1995
1996
1997
1998
1999
2000
2001
2002
2003
2004
2005
2006
2007
2008
2009
2010
2011 adfac=difx*ad_cff
2012 ad_cff1=ad_cff1+ &
2013 & (2.0_r8*cff1*dtdr(i,j,k2)-dtdx(i ,j,k1))* &
2014 & adfac
2015 ad_cff2=ad_cff2+ &
2016 & (2.0_r8*cff2*dtdr(i,j,k2)-dtdx(i+1,j,k2))* &
2017 & adfac
2018 ad_cff3=ad_cff3+ &
2019 & (2.0_r8*cff3*dtdr(i,j,k2)-dtdx(i ,j,k2))* &
2020 & adfac
2021 ad_cff4=ad_cff4+ &
2022 & (2.0_r8*cff4*dtdr(i,j,k2)-dtdx(i+1,j,k1))* &
2023 & adfac
2024 ad_dtdr(i,j,k2)=ad_dtdr(i,j,k2)+ &
2025 & (cff1*cff1+ &
2026 & cff2*cff2+ &
2027 & cff3*cff3+ &
2028 & cff4*cff4)*adfac
2029 ad_dtdx(i ,j,k1)=ad_dtdx(i ,j,k1)-cff1*adfac
2030 ad_dtdx(i+1,j,k2)=ad_dtdx(i+1,j,k2)-cff2*adfac
2031 ad_dtdx(i ,j,k2)=ad_dtdx(i ,j,k2)-cff3*adfac
2032 ad_dtdx(i+1,j,k1)=ad_dtdx(i+1,j,k1)-cff4*adfac
2033 ad_cff=0.0_r8
2034
2035
2036
2037 ad_drdx(i+1,j,k1)=ad_drdx(i+1,j,k1)+ &
2038 & (0.5_r8+sign(0.5_r8, &
2039 & -drdx(i+1,j,k1)))* &
2040 & ad_cff4
2041 ad_cff4=0.0_r8
2042
2043
2044
2045 ad_drdx(i ,j,k2)=ad_drdx(i ,j,k2)+ &
2046 & (0.5_r8+sign(0.5_r8, &
2047 & -drdx(i ,j,k2)))* &
2048 & ad_cff3
2049 ad_cff3=0.0_r8
2050
2051
2052
2053 ad_drdx(i+1,j,k2)=ad_drdx(i+1,j,k2)+ &
2054 & (0.5_r8+sign(0.5_r8, &
2055 & drdx(i+1,j,k2)))* &
2056 & ad_cff2
2057 ad_cff2=0.0_r8
2058
2059
2060
2061 ad_drdx(i ,j,k1)=ad_drdx(i ,j,k1)+ &
2062 & (0.5_r8+sign(0.5_r8, &
2063 & drdx(i ,j,k1)))* &
2064 & ad_cff1
2065 ad_cff1=0.0_r8
2066 END DO
2067 END DO
2068 END IF
2069 DO j=jmin,jmax+1
2070 DO i=imin,imax
2071#ifdef DIFF_3DCOEF
2072# ifdef TS_U3ADV_SPLIT
2073 cff=0.5_r8*diff3d_v(i,j,k)*om_v(i,j)
2074# else
2075 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
2076 & om_v(i,j)
2077# endif
2078#else
2079 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
2080 & om_v(i,j)
2081#endif
2082
2083
2084
2085
2086
2087
2088
2089
2090
2091
2092
2093
2094
2095
2096
2097
2098
2099
2100
2101
2102
2103
2104
2105
2106
2107
2108 adfac=cff*ad_fe(i,j)
2109 adfac1=adfac*(dtde(i,j,k1)- &
2110 & 0.5_r8*(max(drde(i,j,k1),0.0_r8)* &
2111 & (dtdr(i,j-1,k1)+ &
2112 & dtdr(i ,j,k2))+ &
2113 & min(drde(i,j,k1),0.0_r8)* &
2114 & (dtdr(i,j-1,k2)+ &
2115 & dtdr(i ,j,k1))))
2116 adfac2=adfac*(hz(i,j,k)+hz(i,j-1,k))
2117 adfac3=adfac2*0.5_r8*max(drde(i,j,k1),0.0_r8)
2118 adfac4=adfac2*0.5_r8*min(drde(i,j,k1),0.0_r8)
2119 ad_hz(i,j-1,k)=ad_hz(i,j-1,k)+adfac1
2120 ad_hz(i,j ,k)=ad_hz(i,j ,k)+adfac1
2121 ad_dtde(i,j,k1)=ad_dtde(i,j,k1)+adfac2
2122 ad_dtdr(i,j-1,k1)=ad_dtdr(i,j-1,k1)-adfac3
2123 ad_dtdr(i,j, k2)=ad_dtdr(i,j ,k2)-adfac3
2124 ad_dtdr(i,j-1,k2)=ad_dtdr(i,j-1,k2)-adfac4
2125 ad_dtdr(i,j ,k1)=ad_dtdr(i,j ,k1)-adfac4
2126 ad_drde(i,j,k1)=ad_drde(i,j,k1)- &
2127 & adfac2*0.5_r8* &
2128 & ((0.5_r8+sign(0.5_r8, drde(i,j,k1)))* &
2129 & (dtdr(i,j-1,k1)+dtdr(i,j,k2))+ &
2130 & (0.5_r8+sign(0.5_r8,-drde(i,j,k1)))* &
2131 & (dtdr(i,j-1,k2)+dtdr(i,j,k1)))
2132 ad_fe(i,j)=0.0_r8
2133 END DO
2134 END DO
2135 DO j=jmin,jmax
2136 DO i=imin,imax+1
2137#ifdef DIFF_3DCOEF
2138# ifdef TS_U3ADV_SPLIT
2139 cff=0.5_r8*diff3d_u(i,j,k)*on_u(i,j)
2140# else
2141 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
2142 & on_u(i,j)
2143# endif
2144#else
2145 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
2146 & on_u(i,j)
2147#endif
2148
2149
2150
2151
2152
2153
2154
2155
2156
2157
2158
2159
2160
2161
2162
2163
2164
2165
2166
2167
2168
2169
2170
2171
2172
2173
2174 adfac=cff*ad_fx(i,j)
2175 adfac1=adfac*(dtdx(i,j,k1)- &
2176 & 0.5_r8*(max(drdx(i,j,k1),0.0_r8)* &
2177 & (dtdr(i-1,j,k1)+ &
2178 & dtdr(i ,j,k2))+ &
2179 & min(drdx(i,j,k1),0.0_r8)* &
2180 & (dtdr(i-1,j,k2)+ &
2181 & dtdr(i ,j,k1))))
2182 adfac2=adfac*(hz(i,j,k)+hz(i-1,j,k))
2183 adfac3=adfac2*0.5_r8*max(drdx(i,j,k1),0.0_r8)
2184 adfac4=adfac2*0.5_r8*min(drdx(i,j,k1),0.0_r8)
2185 ad_hz(i-1,j,k)=ad_hz(i-1,j,k)+adfac1
2186 ad_hz(i ,j,k)=ad_hz(i ,j,k)+adfac1
2187 ad_dtdx(i,j,k1)=ad_dtdx(i,j,k1)+adfac2
2188 ad_dtdr(i-1,j,k1)=ad_dtdr(i-1,j,k1)-adfac3
2189 ad_dtdr(i ,j,k2)=ad_dtdr(i ,j,k2)-adfac3
2190 ad_dtdr(i-1,j,k2)=ad_dtdr(i-1,j,k2)-adfac4
2191 ad_dtdr(i ,j,k1)=ad_dtdr(i ,j,k1)-adfac4
2192 ad_drdx(i,j,k1)=ad_drdx(i,j,k1)- &
2193 & adfac2*0.5_r8* &
2194 & ((0.5_r8+sign(0.5_r8, drdx(i,j,k1)))* &
2195 & (dtdr(i-1,j,k1)+dtdr(i,j,k2))+ &
2196 & (0.5_r8+sign(0.5_r8,-drdx(i,j,k1)))* &
2197 & (dtdr(i-1,j,k2)+dtdr(i,j,k1)))
2198 ad_fx(i,j)=0.0_r8
2199 END DO
2200 END DO
2201 END IF
2202 IF ((k.eq.0).or.(k.eq.
n(ng)))
THEN
2203 DO j=-1+jmin,jmax+1
2204 DO i=-1+imin,imax+1
2205
2206
2207 ad_dtdr(i,j,k2)=0.0_r8
2208
2209
2210 ad_fs(i,j,k2)=0.0_r8
2211 END DO
2212 END DO
2213 ELSE
2214 DO j=-1+jmin,jmax+1
2215 DO i=-1+imin,imax+1
2216#if defined TS_MIX_MAX_SLOPE
2217 cff1=sqrt(drdx(i,j,k2)**2+drdx(i+1,j,k2)**2+ &
2218 & drdx(i,j,k1)**2+drdx(i+1,j,k1)**2+ &
2219 & drde(i,j,k2)**2+drde(i,j+1,k2)**2+ &
2220 & drde(i,j,k1)**2+drde(i,j+1,k1)**2)
2221 cff2=0.25_r8*slope_max* &
2222 & (z_r(i,j,k+1)-z_r(i,j,k))*cff1
2223 cff3=max(pden(i,j,k)-pden(i,j,k+1),small)
2224 cff4=max(cff2,cff3)
2225 cff=-1.0_r8/cff4
2226#elif defined TS_MIX_MIN_STRAT
2227 cff1=max(pden(i,j,k)-pden(i,j,k+1), &
2228 & strat_min*(z_r(i,j,k+1)-z_r(i,j,k)))
2229 cff=-1.0_r8/cff1
2230#else
2231 cff1=max(pden(i,j,k)-pden(i,j,k+1),eps)
2232 cff=-1.0_r8/cff1
2233#endif
2234
2235
2236
2237
2238
2239 adfac=cff*ad_fs(i,j,k2)
2240 ad_z_r(i,j,k )=ad_z_r(i,j,k )-adfac
2241 ad_z_r(i,j,k+1)=ad_z_r(i,j,k+1)+adfac
2242 ad_cff=ad_cff+(z_r(i,j,k+1)- &
2243 & z_r(i,j,k ))*ad_fs(i,j,k2)
2244 ad_fs(i,j,k2)=0.0_r8
2245#if defined TS_MIX_STABILITY
2246
2247
2248
2249
2250
2251
2252
2253
2254
2255
2256
2257 adfac=cff*ad_dtdr(i,j,k2)
2258 adfac1=adfac*0.75_r8
2259 adfac2=adfac*0.25_r8
2260 ad_t(i,j,k ,nrhs,itrc)=ad_t(i,j,k ,nrhs,itrc)-adfac1
2261 ad_t(i,j,k+1,nrhs,itrc)=ad_t(i,j,k+1,nrhs,itrc)+adfac1
2262 ad_t(i,j,k ,nstp,itrc)=ad_t(i,j,k ,nstp,itrc)-adfac2
2263 ad_t(i,j,k+1,nstp,itrc)=ad_t(i,j,k+1,nstp,itrc)+adfac2
2264 ad_cff=ad_cff+(0.75_r8*(t(i,j,k+1,nrhs,itrc)- &
2265 & t(i,j,k ,nrhs,itrc))+ &
2266 & 0.25_r8*(t(i,j,k+1,nstp,itrc)- &
2267 & t(i,j,k ,nstp,itrc)))* &
2268 & ad_dtdz(i,j,k2)
2269 ad_dtdz(i,j,k2)=0.0_r8
2270#elif defined TS_MIX_CLIMA
2272
2273
2274
2275
2276
2277
2278
2279 adfac=cff*ad_dtdr(i,j,k2)
2280 ad_t(i,j,k ,nrhs,itrc)=ad_t(i,j,k ,nrhs,itrc)-adfac
2281 ad_t(i,j,k+1,nrhs,itrc)=ad_t(i,j,k+1,nrhs,itrc)+adfac
2282 ad_cff=ad_cff+((t(i,j,k+1,nrhs,itrc)- &
2283 & tclm(i,j,k+1,itrc))- &
2284 & (t(i,j,k ,nrhs,itrc)- &
2285 & tclm(i,j,k ,itrc)))*ad_dtdr(i,j,k2)
2286 ad_dtdr(i,j,k2)=0.0_r8
2287 ELSE
2288
2289
2290
2291
2292
2293 END IF
2294#else
2295
2296
2297
2298
2299
2300 adfac=cff*ad_dtdr(i,j,k2)
2301 ad_t(i,j,k ,nrhs,itrc)=ad_t(i,j,k ,nrhs,itrc)-adfac
2302 ad_t(i,j,k+1,nrhs,itrc)=ad_t(i,j,k+1,nrhs,itrc)+adfac
2303 ad_cff=ad_cff+(t(i,j,k+1,nrhs,itrc)- &
2304 & t(i,j,k ,nrhs,itrc))*ad_dtdr(i,j,k2)
2305 ad_dtdr(i,j,k2)=0.0_r8
2306#endif
2307#if defined TS_MIX_MAX_SLOPE
2308
2309
2310 ad_cff4=ad_cff4+cff*cff*ad_cff
2311 ad_cff=0.0_r8
2312
2313
2314
2315 ad_cff3=ad_cff3+ &
2316 & (0.5_r8-sign(0.5_r8,cff2-cff3))*ad_cff4
2317 ad_cff2=ad_cff2+ &
2318 & (0.5_r8+sign(0.5_r8,cff2-cff3))*ad_cff4
2319 ad_cff4=0.0_r8
2320
2321
2322
2323
2324 adfac=(0.5_r8+sign(0.5_r8,pden(i,j,k)-pden(i,j,k+1)- &
2325 & small))*ad_cff3
2326 ad_pden(i,j,k )=ad_pden(i,j,k )+adfac
2327 ad_pden(i,j,k+1)=ad_pden(i,j,k+1)-adfac
2328 ad_cff3=0.0_r8
2329
2330
2331
2332
2333 adfac=0.25_r8*slope_max*ad_cff2
2334 adfac1=adfac*cff1
2335 ad_cff1=ad_cff1+(z_r(i,j,k+1)-z_r(i,j,k))*adfac
2336 ad_z_r(i,j,k )=ad_z_r(i,j,k )-adfac1
2337 ad_z_r(i,j,k+1)=ad_z_r(i,j,k+1)+adfac1
2338 ad_cff2=0.0_r8
2339 IF (cff1.ne.0.0_r8) THEN
2340
2341
2342
2343
2344
2345
2346
2347
2348
2349 adfac=ad_cff1/cff1
2350 ad_drdx(i ,j,k1)=ad_drdx(i ,j,k1)+ &
2351 & drdx(i ,j,k1)*adfac
2352 ad_drdx(i+1,j,k1)=ad_drdx(i+1,j,k1)+ &
2353 & drdx(i+1,j,k1)*adfac
2354 ad_drdx(i ,j,k2)=ad_drdx(i ,j,k2)+ &
2355 & drdx(i ,j,k2)*adfac
2356 ad_drdx(i+1,j,k2)=ad_drdx(i+1,j,k2)+ &
2357 & drdx(i+1,j,k2)*adfac
2358 ad_drde(i,j ,k2)=ad_drde(i,j ,k2)+ &
2359 & drde(i,j ,k2)*adfac
2360 ad_drde(i,j+1,k2)=ad_drde(i,j+1,k2)+ &
2361 & drde(i,j+1,k2)*adfac
2362 ad_drde(i,j ,k1)=ad_drde(i,j ,k1)+ &
2363 & drde(i,j ,k1)*adfac
2364 ad_drde(i,j+1,k1)=ad_drde(i,j+1,k1)+ &
2365 & drde(i,j+1,k1)*adfac
2366 ad_cff1=0.0_r8
2367 ELSE
2368
2369
2370 ad_cff1=0.0_r8
2371 END IF
2372#elif defined TS_MIX_MIN_STRAT
2373
2374
2375 ad_cff1=ad_cff1+cff*cff*ad_cff
2376 ad_cff=0.0_r8
2377
2378
2379
2380
2381
2382
2383
2384
2385
2386
2387
2388 adfac1=(0.5_r8+sign(0.5_r8, &
2389 & pden(i,j,k)-pden(i,j,k+1)- &
2390 & strat_min*(z_r(i,j,k+1)- &
2391 & z_r(i,j,k ))))* &
2392 & ad_cff1
2393 adfac2=(0.5_r8-sign(0.5_r8, &
2394 & pden(i,j,k)-pden(i,j,k+1)- &
2395 & strat_min*(z_r(i,j,k+1)- &
2396 & z_r(i,j,k ))))* &
2397 & strat_min*ad_cff1
2398 ad_pden(i,j,k )=ad_pden(i,j,k )+adfac1
2399 ad_pden(i,j,k+1)=ad_pden(i,j,k+1)-adfac1
2400 ad_z_r(i,j,k )=ad_z_r(i,j,k )-adfac2
2401 ad_z_r(i,j,k+1)=ad_z_r(i,j,k+1)+adfac2
2402 ad_cff1=0.0_r8
2403#else
2404
2405
2406 ad_cff1=ad_cff1+cff*cff*ad_cff
2407 ad_cff=0.0_r8
2408
2409
2410
2411
2412 adfac=(0.5_r8+sign(0.5_r8, &
2413 & pden(i,j,k)-pden(i,j,k+1)-eps))* &
2414 & ad_cff1
2415 ad_pden(i,j,k )=ad_pden(i,j,k )+adfac
2416 ad_pden(i,j,k+1)=ad_pden(i,j,k+1)-adfac
2417 ad_cff1=0.0_r8
2418#endif
2419 END DO
2420 END DO
2421 END IF
2422 IF (k.lt.
n(ng))
THEN
2423 DO j=jmin,jmax+1
2424 DO i=imin,imax
2425 cff=0.5_r8*(pn(i,j)+pn(i,j-1))
2426#ifdef MASKING
2427 cff=cff*vmask(i,j)
2428#endif
2429#ifdef WET_DRY_NOT_YET
2430 cff=cff*vmask_wet(i,j)
2431#endif
2432#if defined TS_MIX_STABILITY
2433
2434
2435
2436
2437
2438
2439 adfac=cff*ad_dtde(i,j,k2)
2440 adfac1=adfac*0.75_r8
2441 adfac2=adfac*0.25_r8
2442 ad_t(i,j-1,k+1,nrhs,itrc)=ad_t(i,j-1,k+1,nrhs,itrc)- &
2443 & adfac1
2444 ad_t(i,j ,k+1,nrhs,itrc)=ad_t(i,j ,k+1,nrhs,itrc)+ &
2445 & adfac1
2446 ad_t(i,j-1,k+1,nstp,itrc)=ad_t(i,j-1,k+1,nstp,itrc)- &
2447 & adfac2
2448 ad_t(i,j ,k+1,nstp,itrc)=ad_t(i,j ,k+1,nstp,itrc)+ &
2449 & adfac2
2450 ad_dtde(i,j,k2)=0.0_r8
2451#elif defined TS_MIX_CLIMA
2452
2453
2454
2455 adfac=cff*ad_dtde(i,j,k2)
2456 ad_t(i,j-1,k+1,nrhs,itrc)=ad_t(i,j-1,k+1,nrhs,itrc)- &
2457 & adfac
2458 ad_t(i,j ,k+1,nrhs,itrc)=ad_t(i,j ,k+1,nrhs,itrc)+ &
2459 & adfac
2460 ad_dtde(i,j,k2)=0.0_r8
2461#else
2462
2463
2464
2465 adfac=cff*ad_dtde(i,j,k2)
2466 ad_t(i,j-1,k+1,nrhs,itrc)=ad_t(i,j-1,k+1,nrhs,itrc)- &
2467 & adfac
2468 ad_t(i,j ,k+1,nrhs,itrc)=ad_t(i,j ,k+1,nrhs,itrc)+ &
2469 & adfac
2470 ad_dtde(i,j,k2)=0.0_r8
2471#endif
2472
2473
2474
2475 adfac=cff*ad_drde(i,j,k2)
2476 ad_pden(i,j-1,k+1)=ad_pden(i,j-1,k+1)-adfac
2477 ad_pden(i,j ,k+1)=ad_pden(i,j ,k+1)+adfac
2478 ad_drde(i,j,k2)=0.0_r8
2479 END DO
2480 END DO
2481 DO j=jmin,jmax
2482 DO i=imin,imax+1
2483 cff=0.5_r8*(pm(i,j)+pm(i-1,j))
2484#ifdef MASKING
2485 cff=cff*umask(i,j)
2486#endif
2487#ifdef WET_DRY_NOT_YET
2488 cff=cff*umask_wet(i,j)
2489#endif
2490#if defined TS_MIX_STABILITY
2491
2492
2493
2494
2495
2496
2497 adfac=cff*ad_dtdx(i,j,k2)
2498 adfac1=adfac*0.75_r8
2499 adfac2=adfac*0.25_r8
2500 ad_t(i-1,j,k+1,nrhs,itrc)=ad_t(i-1,j,k+1,nrhs,itrc)- &
2501 & adfac1
2502 ad_t(i ,j,k+1,nrhs,itrc)=ad_t(i ,j,k+1,nrhs,itrc)+ &
2503 & adfac1
2504 ad_t(i-1,j,k+1,nstp,itrc)=ad_t(i-1,j,k+1,nstp,itrc)- &
2505 & adfac2
2506 ad_t(i ,j,k+1,nstp,itrc)=ad_t(i ,j,k+1,nstp,itrc)+ &
2507 & adfac2
2508 ad_dtdx(i,j,k2)=0.0_r8
2509#elif defined TS_MIX_CLIMA
2510
2511
2512
2513 adfac=cff*ad_dtdx(i,j,k2)
2514 ad_t(i-1,j,k+1,nrhs,itrc)=ad_t(i-1,j,k+1,nrhs,itrc)- &
2515 & adfac
2516 ad_t(i ,j,k+1,nrhs,itrc)=ad_t(i ,j,k+1,nrhs,itrc)+ &
2517 & adfac
2518 ad_dtdx(i,j,k2)=0.0_r8
2519#else
2520
2521
2522
2523 adfac=cff*ad_dtdx(i,j,k2)
2524 ad_t(i-1,j,k+1,nrhs,itrc)=ad_t(i-1,j,k+1,nrhs,itrc)- &
2525 & adfac
2526 ad_t(i ,j,k+1,nrhs,itrc)=ad_t(i ,j,k+1,nrhs,itrc)+ &
2527 & adfac
2528 ad_dtdx(i,j,k2)=0.0_r8
2529#endif
2530
2531
2532
2533 adfac=cff*ad_drdx(i,j,k2)
2534 ad_pden(i-1,j,k+1)=ad_pden(i-1,j,k+1)-adfac
2535 ad_pden(i ,j,k+1)=ad_pden(i ,j,k+1)+adfac
2536 ad_drdx(i,j,k2)=0.0_r8
2537 END DO
2538 END DO
2539 END IF
2540
2541
2542
2543 kt=k2
2544 k2=k1
2545 k1=kt
2546 END DO k_loop3
2547 END DO t_loop
2548
2549 RETURN