152
153
157# if defined NESTING && !defined ONE_WAY
159# endif
162
164# ifdef DISTRIBUTE
167# endif
168 USE mpdata_adiff_mod
169# ifdef NESTING
171# endif
173
174
175
176 integer, intent(in) :: ng, tile
177 integer, intent(in) :: LBi, UBi, LBj, UBj
178 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
179 integer, intent(in) :: nrhs, nstp, nnew
180
181# ifdef ASSUMED_SHAPE
182# ifdef MASKING
183 real(r8), intent(in) :: rmask(LBi:,LBj:)
184 real(r8), intent(in) :: umask(LBi:,LBj:)
185 real(r8), intent(in) :: vmask(LBi:,LBj:)
186# endif
187# ifdef WET_DRY
188 real(r8), intent(in) :: rmask_wet(LBi:,LBj:)
189 real(r8), intent(in) :: umask_wet(LBi:,LBj:)
190 real(r8), intent(in) :: vmask_wet(LBi:,LBj:)
191# endif
192 real(r8), intent(in) :: omn(LBi:,LBj:)
193 real(r8), intent(in) :: om_u(LBi:,LBj:)
194 real(r8), intent(in) :: om_v(LBi:,LBj:)
195 real(r8), intent(in) :: on_u(LBi:,LBj:)
196 real(r8), intent(in) :: on_v(LBi:,LBj:)
197 real(r8), intent(in) :: pm(LBi:,LBj:)
198 real(r8), intent(in) :: pn(LBi:,LBj:)
199 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
200 real(r8), intent(in) :: Huon(LBi:,LBj:,:)
201 real(r8), intent(in) :: Hvom(LBi:,LBj:,:)
202 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
203# ifdef SUN
204 real(r8), intent(in) :: Akt(LBi:UBi,LBj:UBj,0:N(ng),NAT)
205# else
206 real(r8), intent(in) :: Akt(LBi:,LBj:,0:,:)
207# endif
208 real(r8), intent(in) :: W(LBi:,LBj:,0:)
209# ifdef OMEGA_IMPLICIT
210 real(r8), intent(in) :: Wi(LBi:,LBj:,0:)
211# endif
212# ifdef WEC_VF
213 real(r8), intent(in) :: W_stokes(LBi:,LBj:,0:)
214# endif
215# if defined SEDIMENT && defined SED_MORPH
216 real(r8), intent(in) :: bed_thick(LBi:,LBj:,:)
217# endif
218# ifdef DIAGNOSTICS_TS
219 real(r8), intent(inout) :: DiaTwrk(LBi:,LBj:,:,:,:)
220# endif
221# ifdef SUN
222 real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
223# else
224 real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:)
225# endif
226# if defined FLOATS && defined FLOAT_VWALK
227 real(r8), intent(out) :: dAktdz(LBi:,LBj:,:)
228# endif
229
230# else
231
232# ifdef MASKING
233 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
234 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
235 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
236# endif
237# ifdef WET_DRY
238 real(r8), intent(in) :: rmask_wet(LBi:UBi,LBj:UBj)
239 real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
240 real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
241# endif
242 real(r8), intent(in) :: omn(LBi:UBi,LBj:UBj)
243 real(r8), intent(in) :: om_u(LBi:UBi,LBj:UBj)
244 real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
245 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
246 real(r8), intent(in) :: on_v(LBi:UBi,LBj:UBj)
247 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
248 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
249 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
250 real(r8), intent(in) :: Huon(LBi:UBi,LBj:UBj,N(ng))
251 real(r8), intent(in) :: Hvom(LBi:UBi,LBj:UBj,N(ng))
252 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
253 real(r8), intent(in) :: Akt(LBi:UBi,LBj:UBj,0:N(ng),NAT)
254 real(r8), intent(in) :: W(LBi:UBi,LBj:UBj,0:N(ng))
255# ifdef OMEGA_IMPLICIT
256 real(r8), intent(in) :: Wi(LBi:UBi,LBj:UBj,0:N(ng))
257# endif
258# ifdef WEC_VF
259 real(r8), intent(in) :: W_stokes(LBi:UBi,LBj:UBj,0:N(ng))
260# endif
261# if defined SEDIMENT && defined SED_MORPH
262 real(r8), intent(in) :: bed_thick(LBi:UBi,LBj:UBj,3)
263# endif
264# ifdef DIAGNOSTICS_TS
265 real(r8), intent(inout) :: DiaTwrk(LBi:UBi,LBj:UBj,N(ng),NT(ng), &
266 & NDT)
267# endif
268 real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
269
270# if defined FLOATS && defined FLOAT_VWALK
271 real(r8), intent(out) :: dAktdz(LBi:UBi,LBj:UBj,N(ng))
272# endif
273# endif
274
275
276
277 logical :: LapplySrc, Lhsimt, Lmpdata
278
279# ifdef NESTING
280 integer :: ILB, IUB, JLB, JUB
281 integer :: dg, cr, rg
282# endif
283 integer :: IminT, ImaxT, JminT, JmaxT
284 integer :: Isrc, Jsrc
285 integer :: i, ic, ii, is, itrc, j, jj, k, ltrc
286# if defined AGE_MEAN && defined T_PASSIVE
287 integer :: iage
288# endif
289# ifdef DIAGNOSTICS_TS
290 integer :: idiag
291# endif
292 real(r8) :: eps = 1.0e-16_r8
293 real(r8) :: eps1 = 1.0e-12_r8
294
295 real(r8) :: cff, cff1, cff2, cff3
296 real(r8) :: betaL, betaR, betaD, betaU
297 real(r8) :: rL, rR, rD, rU, rkaL, rkaR, rkaD, rkaU
298 real(r8) :: a1, b1, sw, sw_eta, sw_xi
299
300 real(r8), dimension(IminS:ImaxS) :: gradX, KaX, oKaX
301 real(r8), dimension(JminS:JmaxS) :: gradE, KaE, oKaE
302 real(r8), dimension(0:N(ng)) :: gradZ, KaZ, oKaZ
303
304 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: CF
305 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: BC
306 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: DC
307 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
308# ifdef OMEGA_IMPLICIT
309 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FCmin
310 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FCmax
311# endif
312
313 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FE
314 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FX
315 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: curv
316 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: grad
317
318 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: oHz
319
320# ifdef DIAGNOSTICS_TS
321 real(r8), allocatable :: Dhadv(:,:,:)
322 real(r8), allocatable :: Dvadv(:,:,:,:)
323# endif
324 real(r8), allocatable :: Ta(:,:,:,:)
325 real(r8), allocatable :: Ua(:,:,:)
326 real(r8), allocatable :: Va(:,:,:)
327 real(r8), allocatable :: Wa(:,:,:)
328
329# include "set_bounds.h"
330
331# ifdef NESTING
332
333
334
335
340# endif
341
342
343
344
345
350
351
352
353 IF (lmpdata) THEN
354# ifdef DIAGNOSTICS_TS
355 IF (.not.allocated(dhadv)) THEN
356 allocate ( dhadv(imins:imaxs,jmins:jmaxs,3) )
357 dhadv=0.0_r8
358 END IF
359 IF (.not.allocated(dvadv)) THEN
360 allocate ( dvadv(imins:imaxs,jmins:jmaxs,
n(ng),
nt(ng)) )
361 dvadv=0.0_r8
362 END IF
363# endif
364 IF (.not.allocated(ta)) THEN
365 allocate ( ta(imins:imaxs,jmins:jmaxs,
n(ng),
nt(ng)) )
366 ta=0.0_r8
367 END IF
368 IF (.not.allocated(ua)) THEN
369 allocate ( ua(imins:imaxs,jmins:jmaxs,
n(ng)) )
370 ua=0.0_r8
371 END IF
372 IF (.not.allocated(va)) THEN
373 allocate ( va(imins:imaxs,jmins:jmaxs,
n(ng)) )
374 va=0.0_r8
375 END IF
376 IF (.not.allocated(wa)) THEN
377 allocate ( wa(imins:imaxs,jmins:jmaxs,0:
n(ng)) )
378 wa=0.0_r8
379 END IF
380 END IF
381
382
383
384 IF (lmpdata.or.lhsimt) THEN
386 DO j=jstrm2,jendp2
387 DO i=istrm2,iendp2
388 ohz(i,j,k)=1.0_r8/hz(i,j,k)
389 END DO
390 END DO
391 END DO
392 ELSE
394 DO j=jstr,jend
395 DO i=istr,iend
396 ohz(i,j,k)=1.0_r8/hz(i,j,k)
397 END DO
398 END DO
399 END DO
400 END IF
401
402
403
404
405 t_loop1 :
DO itrc=1,
nt(ng)
406
407
408
409
410
415 & lbi, ubi, lbj, ubj, 1,
n(ng), &
416 & t(:,:,:,nnew,itrc))
417 END IF
418
419# ifdef DISTRIBUTE
421 & lbi, ubi, lbj, ubj, 1,
n(ng), &
424 & t(:,:,:,nnew,itrc))
425# endif
426 END IF
427
428
429
430 k_loop :
DO k=1,
n(ng)
431
432 hadv_flux :
IF (
hadvection(itrc,ng)%CENTERED2)
THEN
433
434
435
436 DO j=jstr,jend
437 DO i=istr,iend+1
438 fx(i,j)=huon(i,j,k)* &
439 & 0.5_r8*(t(i-1,j,k,3,itrc)+ &
440 & t(i ,j,k,3,itrc))
441 END DO
442 END DO
443 DO j=jstr,jend+1
444 DO i=istr,iend
445 fe(i,j)=hvom(i,j,k)* &
446 & 0.5_r8*(t(i,j-1,k,3,itrc)+ &
447 & t(i,j ,k,3,itrc))
448 END DO
449 END DO
450
452
453
454
455 DO j=jstrvm2,jendp2i
456 DO i=istrum2,iendp3
457 cff1=max(huon(i,j,k),0.0_r8)
458 cff2=min(huon(i,j,k),0.0_r8)
459 fx(i,j)=cff1*t(i-1,j,k,3,itrc)+ &
460 & cff2*t(i ,j,k,3,itrc)
461 END DO
462 END DO
463 DO j=jstrvm2,jendp3
464 DO i=istrum2,iendp2i
465 cff1=max(hvom(i,j,k),0.0_r8)
466 cff2=min(hvom(i,j,k),0.0_r8)
467 fe(i,j)=cff1*t(i,j-1,k,3,itrc)+ &
468 & cff2*t(i,j ,k,3,itrc)
469 END DO
470 END DO
471
473
474
475
476
477
478
479
480
481
482
483 DO j=jstr,jend
484 DO i=istru-1,iendp2
485 cff=0.125_r8*(pm(i-1,j)+pm(i,j))*(pn(i-1,j)+pn(i,j))* &
487 cff1=cff*(ohz(i-1,j,k)+ohz(i,j,k))
488 gradx(i)=t(i,j,k,3,itrc)-t(i-1,j,k,3,itrc)
489 kax(i)=1.0_r8-abs(huon(i,j,k)*cff1)
490# ifdef MASKING
491 gradx(i)=gradx(i)*umask(i,j)
492 kax(i)=kax(i)*umask(i,j)
493# endif
494 END DO
496 IF (
domain(ng)%Western_Edge(tile))
THEN
497 IF (huon(istr,j,k).ge.0.0_r8) THEN
498 gradx(istr-1)=0.0_r8
499 kax(istr-1)=0.0_r8
500 END IF
501 END IF
502 IF (
domain(ng)%Eastern_Edge(tile))
THEN
503 IF (huon(iend+1,j,k).lt.0.0_r8) THEN
504 gradx(iend+2)=0.0_r8
505 kax(iend+2)=0.0_r8
506 END IF
507 END IF
508 END IF
509 DO i=istr,iend+1
510 IF (kax(i).le.eps1) THEN
511 okax(i)=0.0_r8
512 ELSE
513 okax(i)=1.0_r8/max(kax(i),eps1)
514 END IF
515 IF (huon(i,j,k).ge.0.0_r8) THEN
516 IF (abs(gradx(i)).le.eps1) THEN
517 rl=0.0_r8
518 rkal=0.0_r8
519 ELSE
520 rl=gradx(i-1)/gradx(i)
521 rkal=kax(i-1)*okax(i)
522 END IF
525 betal=a1+b1*rl
526 cff=0.5_r8*max(0.0_r8, &
527 & min(2.0_r8, 2.0_r8*rl*rkal, betal))* &
528 & gradx(i)*kax(i)
529# ifdef MASKING
530 ii=max(i-2,0)
531 cff=cff*rmask(ii,j)
532# endif
533 sw_xi=t(i-1,j,k,3,itrc)+cff
534 ELSE
535 IF (abs(gradx(i)).le.eps1) THEN
536 rr=0.0_r8
537 rkar=0.0_r8
538 ELSE
539 rr=gradx(i+1)/gradx(i)
540 rkar=kax(i+1)*okax(i)
541 END IF
544 betar=a1+b1*rr
545 cff=0.5_r8*max(0.0_r8, &
546 & min(2.0_r8, 2.0_r8*rr*rkar, betar))* &
547 & gradx(i)*kax(i)
548# ifdef MASKING
550 cff=cff*rmask(ii,j)
551# endif
552 sw_xi=t(i,j,k,3,itrc)-cff
553 END IF
554 fx(i,j)=sw_xi*huon(i,j,k)
555 END DO
556 END DO
557
558 DO i=istr,iend
559 DO j=jstrv-1,jendp2
560 cff=0.125_r8*(pn(i,j)+pn(i,j-1))*(pm(i,j)+pm(i,j-1))* &
562 cff1=cff*(ohz(i,j,k)+ohz(i,j-1,k))
563 grade(j)=t(i,j,k,3,itrc)-t(i,j-1,k,3,itrc)
564 kae(j)=1.0_r8-abs(hvom(i,j,k)*cff1)
565# ifdef MASKING
566 grade(j)=grade(j)*vmask(i,j)
567 kae(j)=kae(j)*vmask(i,j)
568# endif
569 END DO
571 IF (
domain(ng)%Southern_Edge(tile))
THEN
572 IF (hvom(i,jstr,k).ge.0.0_r8) THEN
573 grade(jstr-1)=0.0_r8
574 kae(jstr-1)=0.0_r8
575 END IF
576 END IF
577 IF (
domain(ng)%Northern_Edge(tile))
THEN
578 IF (hvom(i,jend+1,k).lt.0.0_r8) THEN
579 grade(jend+2)=0.0_r8
580 kae(jend+2)=0.0_r8
581 END IF
582 END IF
583 END IF
584 DO j=jstr,jend+1
585 IF (kae(j).le.eps1) THEN
586 okae(j)=0.0_r8
587 ELSE
588 okae(j)=1.0_r8/max(kae(j),eps1)
589 END IF
590 IF (hvom(i,j,k).ge.0.0_r8) THEN
591 IF (abs(grade(j)).le.eps1) THEN
592 rd=0.0_r8
593 rkad=0.0_r8
594 ELSE
595 rd=grade(j-1)/grade(j)
596 rkad=kae(j-1)*okae(j)
597 END IF
600 betad=a1+b1*rd
601 cff=0.5_r8*max(0.0_r8, &
602 & min(2.0_r8, 2.0_r8*rd*rkad, betad))* &
603 & grade(j)*kae(j)
604# ifdef MASKING
605 jj=max(j-2,0)
606 cff=cff*rmask(i,jj)
607# endif
608 sw_eta=t(i,j-1,k,3,itrc)+cff
609 ELSE
610 IF (abs(grade(j)).le.eps1) THEN
611 ru=0.0_r8
612 rkau=0.0_r8
613 ELSE
614 ru=grade(j+1)/grade(j)
615 rkau=kae(j+1)*okae(j)
616 END IF
619 betau=a1+b1*ru
620 cff=0.5*max(0.0_r8, &
621 & min(2.0_r8, 2.0_r8*ru*rkau, betau))* &
622 & grade(j)*kae(j)
623# ifdef MASKING
625 cff=cff*rmask(i,jj)
626# endif
627 sw_eta=t(i,j,k,3,itrc)-cff
628 END IF
629 fe(i,j)=sw_eta*hvom(i,j,k)
630 END DO
631 END DO
632
637
638
639
640
641 DO j=jstr,jend
642 DO i=istrm1,iendp2
643 fx(i,j)=t(i ,j,k,3,itrc)- &
644 & t(i-1,j,k,3,itrc)
645# ifdef MASKING
646 fx(i,j)=fx(i,j)*umask(i,j)
647# endif
648 END DO
649 END DO
651 IF (
domain(ng)%Western_Edge(tile))
THEN
652 DO j=jstr,jend
653 fx(istr-1,j)=fx(istr,j)
654 END DO
655 END IF
656 END IF
658 IF (
domain(ng)%Eastern_Edge(tile))
THEN
659 DO j=jstr,jend
660 fx(iend+2,j)=fx(iend+1,j)
661 END DO
662 END IF
663 END IF
664
665 DO j=jstr,jend
666 DO i=istr-1,iend+1
668 curv(i,j)=fx(i+1,j)-fx(i,j)
670 cff=2.0_r8*fx(i+1,j)*fx(i,j)
671 IF (cff.gt.eps) THEN
672 grad(i,j)=cff/(fx(i+1,j)+fx(i,j))
673 ELSE
674 grad(i,j)=0.0_r8
675 END IF
676 ELSE IF ((
hadvection(itrc,ng)%CENTERED4).or. &
678 grad(i,j)=0.5_r8*(fx(i+1,j)+fx(i,j))
679 END IF
680 END DO
681 END DO
682
683 cff1=1.0_r8/6.0_r8
684 cff2=1.0_r8/3.0_r8
685 DO j=jstr,jend
686 DO i=istr,iend+1
688 fx(i,j)=huon(i,j,k)*0.5_r8* &
689 & (t(i-1,j,k,3,itrc)+ &
690 & t(i ,j,k,3,itrc))- &
691 & cff1*(curv(i-1,j)*max(huon(i,j,k),0.0_r8)+ &
692 & curv(i ,j)*min(huon(i,j,k),0.0_r8))
696 fx(i,j)=huon(i,j,k)*0.5_r8* &
697 & (t(i-1,j,k,3,itrc)+ &
698 & t(i ,j,k,3,itrc)- &
699 & cff2*(grad(i ,j)- &
700 & grad(i-1,j)))
701 END IF
702 END DO
703 END DO
704
705 DO j=jstrm1,jendp2
706 DO i=istr,iend
707 fe(i,j)=t(i,j ,k,3,itrc)- &
708 & t(i,j-1,k,3,itrc)
709# ifdef MASKING
710 fe(i,j)=fe(i,j)*vmask(i,j)
711# endif
712 END DO
713 END DO
715 IF (
domain(ng)%Southern_Edge(tile))
THEN
716 DO i=istr,iend
717 fe(i,jstr-1)=fe(i,jstr)
718 END DO
719 END IF
720 END IF
722 IF (
domain(ng)%Northern_Edge(tile))
THEN
723 DO i=istr,iend
724 fe(i,jend+2)=fe(i,jend+1)
725 END DO
726 END IF
727 END IF
728
729 DO j=jstr-1,jend+1
730 DO i=istr,iend
732 curv(i,j)=fe(i,j+1)-fe(i,j)
734 cff=2.0_r8*fe(i,j+1)*fe(i,j)
735 IF (cff.gt.eps) THEN
736 grad(i,j)=cff/(fe(i,j+1)+fe(i,j))
737 ELSE
738 grad(i,j)=0.0_r8
739 END IF
740 ELSE IF ((
hadvection(itrc,ng)%CENTERED4).or. &
742 grad(i,j)=0.5_r8*(fe(i,j+1)+fe(i,j))
743 END IF
744 END DO
745 END DO
746
747 cff1=1.0_r8/6.0_r8
748 cff2=1.0_r8/3.0_r8
749 DO j=jstr,jend+1
750 DO i=istr,iend
752 fe(i,j)=hvom(i,j,k)*0.5_r8* &
753 & (t(i,j-1,k,3,itrc)+ &
754 & t(i,j ,k,3,itrc))- &
755 & cff1*(curv(i,j-1)*max(hvom(i,j,k),0.0_r8)+ &
756 & curv(i,j )*min(hvom(i,j,k),0.0_r8))
760 fe(i,j)=hvom(i,j,k)*0.5_r8* &
761 & (t(i,j-1,k,3,itrc)+ &
762 & t(i,j ,k,3,itrc)- &
763 & cff2*(grad(i,j )- &
764 & grad(i,j-1)))
765 END IF
766 END DO
767 END DO
768 END IF hadv_flux
769
770
771
772
773
774
775
780 IF (int(
sources(ng)%Dsrc(is)).eq.0)
THEN
783 lapplysrc=(istrum2.le.isrc).and. &
784 & (isrc.le.iendp3).and. &
785 & (jstrvm2.le.jsrc).and. &
786 & (jsrc.le.jendp2i)
787 ELSE
788 lapplysrc=(istr.le.isrc).and. &
789 & (isrc.le.iend+1).and. &
790 & (jstr.le.jsrc).and. &
791 & (jsrc.le.jend)
792 END IF
793 IF (lapplysrc) THEN
795 fx(isrc,jsrc)=huon(isrc,jsrc,k)* &
797# ifdef MASKING
798 ELSE
799 IF ((rmask(isrc ,jsrc).eq.0.0_r8).and. &
800 & (rmask(isrc-1,jsrc).eq.1.0_r8)) THEN
801 fx(isrc,jsrc)=huon(isrc,jsrc,k)* &
802 & t(isrc-1,jsrc,k,3,itrc)
803 ELSE IF ((rmask(isrc ,jsrc).eq.1.0_r8).and. &
804 & (rmask(isrc-1,jsrc).eq.0.0_r8)) THEN
805 fx(isrc,jsrc)=huon(isrc,jsrc,k)* &
806 & t(isrc ,jsrc,k,3,itrc)
807 END IF
808# endif
809 END IF
810 END IF
811 ELSE IF (int(
sources(ng)%Dsrc(is)).eq.1)
THEN
814 lapplysrc=(istrum2.le.isrc).and. &
815 & (isrc.le.iendp2i).and. &
816 & (jstrvm2.le.jsrc).and. &
817 & (jsrc.le.jendp3)
818 ELSE
819 lapplysrc=(istr.le.isrc).and. &
820 & (isrc.le.iend).and. &
821 & (jstr.le.jsrc).and. &
822 & (jsrc.le.jend+1)
823 END IF
824 IF (lapplysrc) THEN
826 fe(isrc,jsrc)=hvom(isrc,jsrc,k)* &
828# ifdef MASKING
829 ELSE
830 IF ((rmask(isrc,jsrc ).eq.0.0_r8).and. &
831 & (rmask(isrc,jsrc-1).eq.1.0_r8)) THEN
832 fe(isrc,jsrc)=hvom(isrc,jsrc,k)* &
833 & t(isrc,jsrc-1,k,3,itrc)
834 ELSE IF ((rmask(isrc,jsrc ).eq.1.0_r8).and. &
835 & (rmask(isrc,jsrc-1).eq.0.0_r8)) THEN
836 fe(isrc,jsrc)=hvom(isrc,jsrc,k)* &
837 & t(isrc,jsrc ,k,3,itrc)
838 END IF
839# endif
840 END IF
841 END IF
842 END IF
843 END DO
844 END IF
845
846# if defined NESTING && !defined ONE_WAY
847
848
849
850
851
856 IF (ng.eq.rg) THEN
858 & imins, imaxs, jmins, jmaxs, &
859 & ilb, iub, jlb, jub, &
865 END IF
866 END DO
867 END IF
868# endif
869
870
871
872
873 hadv_stepping :
IF (
hadvection(itrc,ng)%MPDATA)
THEN
874 DO j=jstrvm2,jendp2i
875 DO i=istrum2,iendp2i
876 cff=
dt(ng)*pm(i,j)*pn(i,j)
877 cff1=cff*(fx(i+1,j)-fx(i,j))
878 cff2=cff*(fe(i,j+1)-fe(i,j))
879 cff3=cff1+cff2
880 ta(i,j,k,itrc)=t(i,j,k,nnew,itrc)-cff3
881# ifdef DIAGNOSTICS_TS
885# endif
886 END DO
887 END DO
888# ifdef DIAGNOSTICS_TS
889
890 DO j=jstr,jend
891 DO i=istr,iend
895 END DO
896 END DO
897# endif
898
899
900
901
902 ELSE
903 DO j=jstr,jend
904 DO i=istr,iend
905 cff=
dt(ng)*pm(i,j)*pn(i,j)
906 cff1=cff*(fx(i+1,j)-fx(i,j))
907 cff2=cff*(fe(i,j+1)-fe(i,j))
908 cff3=cff1+cff2
909 t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)-cff3
910# ifdef DIAGNOSTICS_TS
911 diatwrk(i,j,k,itrc,
itxadv)=-cff1
912 diatwrk(i,j,k,itrc,
ityadv)=-cff2
913 diatwrk(i,j,k,itrc,
ithadv)=-cff3
914# endif
915 END DO
916 END DO
917 END IF hadv_stepping
918 END DO k_loop
919 END DO t_loop1
920
921
922
923
924
925 t_loop2 :
DO itrc=1,
nt(ng)
927 jmint=jstrvm2
928 jmaxt=jendp2i
929 ELSE
930 jmint=jstr
931 jmaxt=jend
932 END IF
933
934 j_loop1 : DO j=jmint,jmaxt
935
936 vadv_flux :
IF (
vadvection(itrc,ng)%SPLINES)
THEN
937
938
939
940
941
942 DO i=istr,iend
943# ifdef NEUMANN
944 fc(i,0)=1.5_r8*t(i,j,1,3,itrc)
945 cf(i,1)=0.5_r8
946# else
947 fc(i,0)=2.0_r8*t(i,j,1,3,itrc)
948 cf(i,1)=1.0_r8
949# endif
950 END DO
952 DO i=istr,iend
953 cff=1.0_r8/(2.0_r8*hz(i,j,k)+ &
954 & hz(i,j,k+1)*(2.0_r8-cf(i,k)))
955 cf(i,k+1)=cff*hz(i,j,k)
956 fc(i,k)=cff*(3.0_r8*(hz(i,j,k )*t(i,j,k+1,3,itrc)+ &
957 & hz(i,j,k+1)*t(i,j,k ,3,itrc))- &
958 & hz(i,j,k+1)*fc(i,k-1))
959 END DO
960 END DO
961 DO i=istr,iend
962# ifdef NEUMANN
963 fc(i,
n(ng))=(3.0_r8*t(i,j,
n(ng),3,itrc)-fc(i,
n(ng)-1))/ &
964 & (2.0_r8-cf(i,
n(ng)))
965# else
966 fc(i,
n(ng))=(2.0_r8*t(i,j,
n(ng),3,itrc)-fc(i,
n(ng)-1))/ &
967 & (1.0_r8-cf(i,
n(ng)))
968# endif
969 END DO
971 DO i=istr,iend
972 fc(i,k)=fc(i,k)-cf(i,k+1)*fc(i,k+1)
973# ifdef WEC_VF
974 fc(i,k+1)=(w(i,j,k+1)+w_stokes(i,j,k+1))*fc(i,k+1)
975# else
976 fc(i,k+1)=w(i,j,k+1)*fc(i,k+1)
977# endif
978 END DO
979 END DO
980 DO i=istr,iend
982 fc(i,0)=0.0_r8
983 END DO
984
986
987
988
990 DO i=istr,iend
991 fc(i,k)=t(i,j,k+1,3,itrc)- &
992 & t(i,j,k ,3,itrc)
993 END DO
994 END DO
995 DO i=istr,iend
996 fc(i,0)=fc(i,1)
997 fc(i,
n(ng))=fc(i,
n(ng)-1)
998 END DO
1000 DO i=istr,iend
1001 cff=2.0_r8*fc(i,k)*fc(i,k-1)
1002 IF (cff.gt.eps) THEN
1003 cf(i,k)=cff/(fc(i,k)+fc(i,k-1))
1004 ELSE
1005 cf(i,k)=0.0_r8
1006 END IF
1007 END DO
1008 END DO
1009 cff1=1.0_r8/3.0_r8
1011 DO i=istr,iend
1012# ifdef WEC_VF
1013 fc(i,k)=(w(i,j,k)+w_stokes(i,j,k))* &
1014# else
1015 fc(i,k)=w(i,j,k)* &
1016# endif
1017 & 0.5_r8*(t(i,j,k ,3,itrc)+ &
1018 & t(i,j,k+1,3,itrc)- &
1019 & cff1*(cf(i,k+1)-cf(i,k)))
1020 END DO
1021 END DO
1022 DO i=istr,iend
1023 fc(i,0)=0.0_r8
1025 END DO
1026
1028
1029
1030
1031
1033 DO i=istr,iend
1034# ifdef WEC_VF
1035 fc(i,k)=(w(i,j,k)+w_stokes(i,j,k))* &
1036# else
1037 fc(i,k)=w(i,j,k)* &
1038# endif
1039 & 0.5_r8*(t(i,j,k ,3,itrc)+ &
1040 & t(i,j,k+1,3,itrc))
1041 END DO
1042 END DO
1043 DO i=istr,iend
1044 fc(i,0)=0.0_r8
1046 END DO
1047
1049
1050
1051
1052
1053 DO i=istrum2,iendp2i
1055# ifdef WEC_VF
1056 cff1=max(w(i,j,k)+w_stokes(i,j,k),0.0_r8)
1057 cff2=min(w(i,j,k)+w_stokes(i,j,k),0.0_r8)
1058# else
1059 cff1=max(w(i,j,k),0.0_r8)
1060 cff2=min(w(i,j,k),0.0_r8)
1061# endif
1062 fc(i,k)=cff1*t(i,j,k ,3,itrc)+ &
1063 & cff2*t(i,j,k+1,3,itrc)
1064 END DO
1065 fc(i,0)=0.0_r8
1067 END DO
1068
1070
1071
1072
1073
1074
1075 DO i=istr,iend
1076 kaz(0)=0.0_r8
1077 okaz(0)=0.0_r8
1078 gradz(0)=0.0_r8
1080 cff=pm(i,j)*pn(i,j)*
dt(ng)
1081# ifdef WEC_VF
1082 kaz(k)=1.0_r8-abs(cff*(w(i,j,k)+w_stokes(i,j,k))/ &
1083 & (z_r(i,j,k+1)-z_r(i,j,k)))
1084# else
1085 kaz(k)=1.0_r8-abs(cff*w(i,j,k)/ &
1086 & (z_r(i,j,k+1)-z_r(i,j,k)))
1087# endif
1088 okaz(k)=1.0_r8/kaz(k)
1089 gradz(k)=t(i,j,k+1,3,itrc)-t(i,j,k,3,itrc)
1090 END DO
1094
1096# ifdef WEC_VF
1097 cff1=w(i,j,k)+w_stokes(i,j,k)
1098# else
1099 cff1=w(i,j,k)
1100# endif
1101 IF ((k.eq.1).and.(cff1.ge.0.0_r8)) THEN
1102 fc(i,k)=cff1*t(i,j,k,3,itrc)
1103 ELSE IF ((k.eq.
n(ng)-1).and.(cff1.lt.0.0_r8))
THEN
1104 fc(i,k)=cff1*t(i,j,k+1,3,itrc)
1105 ELSE
1106 IF (cff1.ge.0) THEN
1107 IF (abs(gradz(k)).le.eps1) THEN
1108 rd=0.0_r8
1109 rkad=0.0_r8
1110 ELSE
1111 rd=gradz(k-1)/gradz(k)
1112 rkad=kaz(k-1)*okaz(k)
1113 END IF
1116 betad=a1+b1*rd
1117 cff=0.5_r8*max(0.0_r8, &
1118 & min(2.0_r8, 2.0_r8*rd*rkad, betad))* &
1119 & gradz(k)*kaz(k)
1120 sw=t(i,j,k,3,itrc)+cff
1121 ELSE
1122 IF (abs(gradz(k)).le.eps1) THEN
1123 ru=0.0_r8
1124 rkau=0.0_r8
1125 ELSE
1126 ru=gradz(k+1)/gradz(k)
1127 rkau=kaz(k+1)*okaz(k)
1128 END IF
1131 betau=a1+b1*ru
1132 cff=0.5_r8*max(0.0_r8, &
1133 & min(2.0_r8, 2.0_r8*ru*rkau, betau))* &
1134 & gradz(k)*kaz(k)
1135 sw=t(i,j,k+1,3,itrc)-cff
1136 END IF
1137 fc(i,k)=cff1*sw
1138 END IF
1139 END DO
1140 fc(i,0)=0.0_r8
1142 END DO
1143
1144 ELSE IF ((
vadvection(itrc,ng)%CENTERED4).or. &
1146
1147
1148
1149
1150 cff1=0.5_r8
1151 cff2=7.0_r8/12.0_r8
1152 cff3=1.0_r8/12.0_r8
1154 DO i=istr,iend
1155# ifdef WEC_VF
1156 fc(i,k)=(w(i,j,k)+w_stokes(i,j,k))* &
1157# else
1158 fc(i,k)=w(i,j,k)* &
1159# endif
1160 & (cff2*(t(i,j,k ,3,itrc)+ &
1161 & t(i,j,k+1,3,itrc))- &
1162 & cff3*(t(i,j,k-1,3,itrc)+ &
1163 & t(i,j,k+2,3,itrc)))
1164 END DO
1165 END DO
1166 DO i=istr,iend
1167 fc(i,0)=0.0_r8
1168# ifdef WEC_VF
1169 fc(i,1)=(w(i,j,1)+w_stokes(i,j,1))* &
1170# else
1171 fc(i,1)=w(i,j,1)* &
1172# endif
1173 & (cff1*t(i,j,1,3,itrc)+ &
1174 & cff2*t(i,j,2,3,itrc)- &
1175 & cff3*t(i,j,3,3,itrc))
1176# ifdef WEC_VF
1177 fc(i,
n(ng)-1)=(w(i,j,
n(ng)-1)+w_stokes(i,j,
n(ng)-1))* &
1178# else
1179 fc(i,
n(ng)-1)=w(i,j,
n(ng)-1)* &
1180# endif
1181 & (cff1*t(i,j,
n(ng) ,3,itrc)+ &
1182 & cff2*t(i,j,
n(ng)-1,3,itrc)- &
1183 & cff3*t(i,j,
n(ng)-2,3,itrc))
1185 END DO
1186 END IF vadv_flux
1187
1188
1189
1190
1191
1192
1193
1194
1198 IF (int(
sources(ng)%Dsrc(is)).eq.2)
THEN
1201 IF (((istr.le.isrc).and.(isrc.le.iend+1)).and. &
1202 & ((jstr.le.jsrc).and.(jsrc.le.jend+1)).and. &
1203 & (j.eq.jsrc)) THEN
1205 cff=
dt(ng)*pm(isrc,jsrc)*pn(isrc,jsrc)
1207 cff3=
sources(ng)%Tsrc(is,k,itrc)
1208 ELSE
1209 cff3=t(isrc,jsrc,k,3,itrc)
1210 END IF
1211 ta(isrc,jsrc,k,itrc)=ta(isrc,jsrc,k,itrc)+ &
1212 & cff*
sources(ng)%Qsrc(is,k)* &
1213 & cff3
1214 END DO
1215 END IF
1216 END IF
1217 END DO
1218 END IF
1219 END IF
1220
1221# ifdef SED_MORPH
1222
1223
1224
1225
1226
1227
1231 DO i=istrum2,iendp2i
1232 cff=cff1*(bed_thick(i,j,nstp)-bed_thick(i,j,3))
1233 cff3=t(i,j,k,3,itrc)
1234 ta(i,j,k,itrc)=ta(i,j,k,itrc)-cff*cff3
1235 END DO
1236 END DO
1237 END IF
1238# endif
1239
1240
1241
1242# ifdef DIAGNOSTICS_TS
1243
1244# endif
1245
1246 vadv_stepping :
IF (
vadvection(itrc,ng)%MPDATA)
THEN
1247 DO i=istrum2,iendp2i
1248 cf(i,0)=
dt(ng)*pm(i,j)*pn(i,j)
1249 END DO
1251 DO i=istrum2,iendp2i
1252 cff1=cf(i,0)*(fc(i,k)-fc(i,k-1))
1253 ta(i,j,k,itrc)=(ta(i,j,k,itrc)-cff1)*ohz(i,j,k)
1254# ifdef DIAGNOSTICS_TS
1255 dvadv(i,j,k,itrc)=-cff1
1256# endif
1257 END DO
1258 END DO
1259
1260# ifdef OMEGA_IMPLICIT
1261
1262
1263
1264
1265
1266
1267
1268
1271 DO i=istr,iend
1272 cff1=cff*pm(i,j)*pn(i,j)
1273 fcmax(i,k)=max(wi(i,j,k),0.0_r8)*cff1
1274 fcmin(i,k)=min(wi(i,j,k),0.0_r8)*cff1
1275 END DO
1276 END DO
1277 DO i=istr,iend
1278 fcmax(i,0)=0.0_r8
1279 fcmin(i,0)=0.0_r8
1280 fcmax(i,
n(ng))=0.0_r8
1281 fcmin(i,
n(ng))=0.0_r8
1282 END DO
1283
1284
1285
1286
1288 DO i=istr,iend
1289 bc(i,k)=hz(i,j,k)+fcmax(i,k)-fcmin(i,k-1)
1290 dc(i,k)=ta(i,j,k,itrc)*hz(i,j,k)
1291 END DO
1292 END DO
1293
1294
1295
1296 DO i=istr,iend
1297 cff=1.0_r8/bc(i,1)
1298 cf(i,1)=cff*fcmin(i,1)
1299 dc(i,1)=cff*dc(i,1)
1300 END DO
1302 DO i=istr,iend
1303 cff=1.0_r8/(bc(i,k)+fcmax(i,k-1)*cf(i,k-1))
1304 cf(i,k)=cff*fcmin(i,k)
1305 dc(i,k)=cff*(dc(i,k)+fcmax(i,k-1)*dc(i,k-1))
1306 END DO
1307 END DO
1308
1309
1310
1311 DO i=istr,iend
1312# ifdef DIAGNOSTICS_TS
1313 cff1=ta(i,j,
n(ng),itrc)*ohz(i,j,
n(ng))
1314# endif
1315 cff=1.0_r8/(bc(i,
n(ng))+fcmax(i,
n(ng)-1)*cf(i,
n(ng)-1))
1316 dc(i,
n(ng))=cff*(dc(i,
n(ng))+ &
1317 & fcmax(i,
n(ng)-1)*dc(i,
n(ng)-1))
1318 ta(i,j,
n(ng),itrc)=dc(i,
n(ng))
1319# ifdef DIAGNOSTICS_TS
1320 diatwrk(i,j,
n(ng),itrc,
itvadv)= &
1321 & diatwrk(i,j,
n(ng),itrc,
itvadv)+ &
1322 & ta(i,j,
n(ng),itrc)-cff1
1323# endif
1324 END DO
1326 DO i=istr,iend
1327# ifdef DIAGNOSTICS_TS
1328 cff1=ta(i,j,k,itrc)*ohz(i,j,k)
1329# endif
1330 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
1331 ta(i,j,k,itrc)=dc(i,k)
1332# ifdef DIAGNOSTICS_TS
1333 diatwrk(i,j,k,itrc,
itvadv)=diatwrk(i,j,k,itrc,
itvadv)+&
1334 & ta(i,j,k,itrc)-cff1
1335# endif
1336 END DO
1337 END DO
1338# endif
1339
1340
1341
1342# ifdef DIAGNOSTICS_TS
1343
1344# endif
1345
1346 ELSE
1347 DO i=istr,iend
1348 cf(i,0)=
dt(ng)*pm(i,j)*pn(i,j)
1349 END DO
1351 DO i=istr,iend
1352 cff1=cf(i,0)*(fc(i,k)-fc(i,k-1))
1353 t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)-cff1
1354# ifdef SPLINES_VDIFF
1355 t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)*ohz(i,j,k)
1356# endif
1357# ifdef DIAGNOSTICS_TS
1358 diatwrk(i,j,k,itrc,
itvadv)=-cff1
1360 diatwrk(i,j,k,itrc,idiag)=diatwrk(i,j,k,itrc,idiag)* &
1361 & ohz(i,j,k)
1362 END DO
1363# endif
1364 END DO
1365 END DO
1366 END IF vadv_stepping
1367 END DO j_loop1
1368 END DO t_loop2
1369
1370
1371
1372
1373
1374
1375 t_loop3 :
DO itrc=1,
nt(ng)
1376 mpdata :
IF ((
hadvection(itrc,ng)%MPDATA).and. &
1379 & lbi, ubi, lbj, ubj, &
1380 & imins, imaxs, jmins, jmaxs, &
1381# ifdef MASKING
1382 & rmask, umask, vmask, &
1383# endif
1384# ifdef WET_DRY
1385 & rmask_wet, umask_wet, vmask_wet, &
1386# endif
1387 & pm, pn, omn, om_u, on_v, &
1388 & z_r, ohz, &
1389 & huon, hvom, w, &
1390# ifdef WEC_VF
1391 & w_stokes, &
1392# endif
1393# ifdef OMEGA_IMPLICIT
1394 & wi, &
1395# endif
1396 & t(:,:,:,3,itrc), &
1397 & ta(:,:,:,itrc), ua, va, wa)
1398
1399
1400
1402 DO j=jstr,jend
1403 DO i=istr,iend+1
1404 cff1=max(ua(i,j,k),0.0_r8)
1405 cff2=min(ua(i,j,k),0.0_r8)
1406 fx(i,j)=(cff1*ta(i-1,j,k,itrc)+ &
1407 & cff2*ta(i ,j,k,itrc))* &
1408 & 0.5_r8*(hz(i,j,k)+hz(i-1,j,k))*on_u(i,j)
1409 END DO
1410 END DO
1411 DO j=jstr,jend+1
1412 DO i=istr,iend
1413 cff1=max(va(i,j,k),0.0_r8)
1414 cff2=min(va(i,j,k),0.0_r8)
1415 fe(i,j)=(cff1*ta(i,j-1,k,itrc)+ &
1416 & cff2*ta(i,j ,k,itrc))* &
1417 & 0.5_r8*(hz(i,j,k)+hz(i,j-1,k))*om_v(i,j)
1418 END DO
1419 END DO
1420
1421
1422
1423 DO j=jstr,jend
1424 DO i=istr,iend
1425 cff=
dt(ng)*pm(i,j)*pn(i,j)
1426 cff1=cff*(fx(i+1,j)-fx(i,j))
1427 cff2=cff*(fe(i,j+1)-fe(i,j))
1428 cff3=cff1+cff2
1429 t(i,j,k,nnew,itrc)=ta(i,j,k,itrc)*hz(i,j,k)-cff3
1430# ifdef DIAGNOSTICS_TS
1431 diatwrk(i,j,k,itrc,
itxadv)=diatwrk(i,j,k,itrc,
itxadv)- &
1432 & cff1
1433 diatwrk(i,j,k,itrc,
ityadv)=diatwrk(i,j,k,itrc,
ityadv)- &
1434 & cff2
1435 diatwrk(i,j,k,itrc,
ithadv)=diatwrk(i,j,k,itrc,
ithadv)- &
1436 & cff3
1437# endif
1438 END DO
1439 END DO
1440 END DO
1441
1442
1443
1444 DO j=jstr,jend
1446 DO i=istr,iend
1447 cff1=max(wa(i,j,k),0.0_r8)
1448 cff2=min(wa(i,j,k),0.0_r8)
1449 fc(i,k)=cff1*ta(i,j,k ,itrc)+ &
1450 & cff2*ta(i,j,k+1,itrc)
1451 END DO
1452 END DO
1453 DO i=istr,iend
1454 fc(i,0)=0.0_r8
1456 END DO
1457
1458
1459# ifdef DIAGNOSTICS_TS
1460
1461# endif
1462
1463 DO i=istr,iend
1464 cf(i,0)=
dt(ng)*pm(i,j)*pn(i,j)
1465 END DO
1467 DO i=istr,iend
1468 cff1=cf(i,0)*(fc(i,k)-fc(i,k-1))
1469 t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)-cff1
1470# ifdef DIAGNOSTICS_TS
1471 diatwrk(i,j,k,itrc,
itvadv)=dvadv(i,j,k,itrc)- &
1472 & cff1
1474 diatwrk(i,j,k,itrc,idiag)=diatwrk(i,j,k,itrc,idiag)* &
1475 & ohz(i,j,k)
1476 END DO
1477# endif
1478 END DO
1479 END DO
1480 END DO
1481 END IF mpdata
1482 END DO t_loop3
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1497 IF (.not.((
hadvection(itrc,ng)%MPDATA).and. &
1500 IF (int(
sources(ng)%Dsrc(is)).eq.2)
THEN
1503 IF (((istr.le.isrc).and.(isrc.le.iend+1)).and. &
1504 & ((jstr.le.jsrc).and.(jsrc.le.jend+1))) THEN
1506 cff=
dt(ng)*pm(isrc,jsrc)*pn(isrc,jsrc)
1507# ifdef SPLINES_VDIFF
1508 cff=cff*ohz(isrc,jsrc,k)
1509# endif
1511 cff3=
sources(ng)%Tsrc(is,k,itrc)
1512 ELSE
1513 cff3=t(isrc,jsrc,k,3,itrc)
1514 END IF
1515 t(isrc,jsrc,k,nnew,itrc)=t(isrc,jsrc,k,nnew,itrc)+ &
1516 & cff*
sources(ng)%Qsrc(is,k)*&
1517 & cff3
1518 END DO
1519 END IF
1520 END IF
1521 END DO
1522 END IF
1523 END DO
1524 END IF
1525
1526# if defined SEDIMENT && defined SED_MORPH
1527
1528
1529
1530
1531
1534 IF (.not.((
hadvection(itrc,ng)%MPDATA).and. &
1536 DO j=jstr,jend
1538 DO i=istr,iend
1539 cff3=t(i,j,k,3,itrc)
1540# ifdef SPLINES_VDIFF
1541 cff3=cff3*ohz(i,j,k)
1542# endif
1543 cff=cff1*(bed_thick(i,j,nstp)-bed_thick(i,j,3))
1544 t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)- &
1545 & cff*cff3
1546 END DO
1547 END DO
1548 END DO
1549 END IF
1550 END DO
1551# endif
1552
1553# ifdef OMEGA_IMPLICIT
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565 DO j=jstr,jend
1570 DO i=istr,iend
1571 cff1=cff*pm(i,j)*pn(i,j)
1572 fcmax(i,k)=max(wi(i,j,k),0.0_r8)*cff1
1573 fcmin(i,k)=min(wi(i,j,k),0.0_r8)*cff1
1574 END DO
1575 END DO
1576 DO i=istr,iend
1577 fcmax(i,0)=0.0_r8
1578 fcmin(i,0)=0.0_r8
1579 fcmax(i,
n(ng))=0.0_r8
1580 fcmin(i,
n(ng))=0.0_r8
1581 END DO
1582
1583
1584
1585
1587 DO i=istr,iend
1588 bc(i,k)=hz(i,j,k)+fcmax(i,k)-fcmin(i,k-1)
1589# ifdef SPLINES_VDIFF
1590 dc(i,k)=t(i,j,k,nnew,itrc)*hz(i,j,k)
1591# else
1592 dc(i,k)=t(i,j,k,nnew,itrc)
1593# endif
1594 END DO
1595 END DO
1596
1597
1598
1599 DO i=istr,iend
1600 cff=1.0_r8/bc(i,1)
1601 cf(i,1)=cff*fcmin(i,1)
1602 dc(i,1)=cff*dc(i,1)
1603 END DO
1605 DO i=istr,iend
1606 cff=1.0_r8/(bc(i,k)+fcmax(i,k-1)*cf(i,k-1))
1607 cf(i,k)=cff*fcmin(i,k)
1608 dc(i,k)=cff*(dc(i,k)+fcmax(i,k-1)*dc(i,k-1))
1609 END DO
1610 END DO
1611
1612
1613
1614 DO i=istr,iend
1615# ifdef DIAGNOSTICS_TS
1616 cff1=t(i,j,
n(ng),nnew,itrc)*ohz(i,j,
n(ng))
1617# endif
1618 cff=1.0_r8/(bc(i,
n(ng))+ &
1619 & fcmax(i,
n(ng)-1)*cf(i,
n(ng)-1))
1620 dc(i,
n(ng))=cff*(dc(i,
n(ng))+ &
1621 & fcmax(i,
n(ng)-1)*dc(i,
n(ng)-1))
1622# ifdef SPLINES_VDIFF
1623 t(i,j,
n(ng),nnew,itrc)=dc(i,
n(ng))
1624# else
1625 t(i,j,
n(ng),nnew,itrc)=dc(i,
n(ng))*hz(i,j,
n(ng))
1626# endif
1627# ifdef DIAGNOSTICS_TS
1628 diatwrk(i,j,
n(ng),itrc,
itvadv)=diatwrk(i,j,
n(ng),itrc, &
1631# endif
1632 END DO
1633
1635 DO i=istr,iend
1636# ifdef DIAGNOSTICS_TS
1637 cff1=t(i,j,k,nnew,itrc)*ohz(i,j,k)
1638# endif
1639 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
1640# ifdef SPLINES_VDIFF
1641 t(i,j,k,nnew,itrc)=dc(i,k)
1642# else
1643 t(i,j,k,nnew,itrc)=dc(i,k)*hz(i,j,k)
1644# endif
1645# ifdef DIAGNOSTICS_TS
1646 diatwrk(i,j,k,itrc,
itvadv)=diatwrk(i,j,k,itrc,
itvadv)+ &
1647 & dc(i,k)-cff1
1648# endif
1649 END DO
1650 END DO
1651 END IF
1652 END DO
1653 END DO
1654# endif
1655
1656
1657
1658
1659
1660 j_loop2 : DO j=jstr,jend
1663
1664# ifdef SPLINES_VDIFF
1665 IF (.not.((
hadvection(itrc,ng)%MPDATA).and. &
1667
1668
1669
1670
1671
1672 cff1=1.0_r8/6.0_r8
1674 DO i=istr,iend
1675 fc(i,k)=cff1*hz(i,j,k )- &
1676 &
dt(ng)*akt(i,j,k-1,ltrc)*ohz(i,j,k )
1677 cf(i,k)=cff1*hz(i,j,k+1)- &
1678 &
dt(ng)*akt(i,j,k+1,ltrc)*ohz(i,j,k+1)
1679 END DO
1680 END DO
1681 DO i=istr,iend
1682 cf(i,0)=0.0_r8
1683 dc(i,0)=0.0_r8
1684 END DO
1685
1686
1687
1688 cff1=1.0_r8/3.0_r8
1690 DO i=istr,iend
1691 bc(i,k)=cff1*(hz(i,j,k)+hz(i,j,k+1))+ &
1692 &
dt(ng)*akt(i,j,k,ltrc)*(ohz(i,j,k)+ohz(i,j,k+1))
1693 cff=1.0_r8/(bc(i,k)-fc(i,k)*cf(i,k-1))
1694 cf(i,k)=cff*cf(i,k)
1695 dc(i,k)=cff*(t(i,j,k+1,nnew,itrc)-t(i,j,k,nnew,itrc)- &
1696 & fc(i,k)*dc(i,k-1))
1697 END DO
1698 END DO
1699
1700
1701
1702 DO i=istr,iend
1704 END DO
1706 DO i=istr,iend
1707 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
1708 END DO
1709 END DO
1710
1712 DO i=istr,iend
1713 dc(i,k)=dc(i,k)*akt(i,j,k,ltrc)
1714 cff1=
dt(ng)*ohz(i,j,k)*(dc(i,k)-dc(i,k-1))
1715 t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)+cff1
1716# ifdef DIAGNOSTICS_TS
1717 diatwrk(i,j,k,itrc,
itvdif)=diatwrk(i,j,k,itrc,
itvdif)+ &
1718 & cff1
1719# endif
1720 END DO
1721 END DO
1722 ELSE
1723# endif
1724
1725
1726
1727
1728
1729
1732 DO i=istr,iend
1733 cff1=1.0_r8/(z_r(i,j,k+1)-z_r(i,j,k))
1734 fc(i,k)=cff*cff1*akt(i,j,k,ltrc)
1735 END DO
1736 END DO
1737 DO i=istr,iend
1738 fc(i,0)=0.0_r8
1740 END DO
1741
1742
1743
1744
1746 DO i=istr,iend
1747 bc(i,k)=hz(i,j,k)-fc(i,k)-fc(i,k-1)
1748 dc(i,k)=t(i,j,k,nnew,itrc)
1749 END DO
1750 END DO
1751
1752
1753
1754 DO i=istr,iend
1755 cff=1.0_r8/bc(i,1)
1756 cf(i,1)=cff*fc(i,1)
1757 dc(i,1)=cff*dc(i,1)
1758 END DO
1760 DO i=istr,iend
1761 cff=1.0_r8/(bc(i,k)-fc(i,k-1)*cf(i,k-1))
1762 cf(i,k)=cff*fc(i,k)
1763 dc(i,k)=cff*(dc(i,k)-fc(i,k-1)*dc(i,k-1))
1764 END DO
1765 END DO
1766
1767
1768
1769 DO i=istr,iend
1770# ifdef DIAGNOSTICS_TS
1771 cff1=t(i,j,
n(ng),nnew,itrc)*ohz(i,j,
n(ng))
1772# endif
1773 dc(i,
n(ng))=(dc(i,
n(ng))-fc(i,
n(ng)-1)*dc(i,
n(ng)-1))/ &
1774 & (bc(i,
n(ng))-fc(i,
n(ng)-1)*cf(i,
n(ng)-1))
1775 t(i,j,
n(ng),nnew,itrc)=dc(i,
n(ng))
1776# ifdef DIAGNOSTICS_TS
1777 diatwrk(i,j,
n(ng),itrc,
itvdif)= &
1778 & diatwrk(i,j,
n(ng),itrc,
itvdif)+ &
1779 & t(i,j,
n(ng),nnew,itrc)-cff1
1780# endif
1781 END DO
1783 DO i=istr,iend
1784# ifdef DIAGNOSTICS_TS
1785 cff1=t(i,j,k,nnew,itrc)*ohz(i,j,k)
1786# endif
1787 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
1788 t(i,j,k,nnew,itrc)=dc(i,k)
1789# ifdef DIAGNOSTICS_TS
1790 diatwrk(i,j,k,itrc,
itvdif)=diatwrk(i,j,k,itrc,
itvdif)+ &
1791 & t(i,j,k,nnew,itrc)-cff1
1792# endif
1793 END DO
1794 END DO
1795# ifdef SPLINES_VDIFF
1796 END IF
1797# endif
1798 END DO
1799 END DO j_loop2
1800
1801# if defined AGE_MEAN && defined T_PASSIVE
1802
1803
1804
1805
1806
1807
1808
1809
1810
1811
1812
1813
1814
1815
1816
1820 DO j=jstr,jend
1821 DO i=istr,iend
1824 t(i,j,k,nnew,iage)=t(i,j,k,nnew,iage)+ &
1825 &
dt(ng)*t(i,j,k,nnew,
inert(itrc))
1826 ELSE
1827 t(i,j,k,nnew,iage)=t(i,j,k,nnew,iage)+ &
1828 &
dt(ng)*t(i,j,k,3,
inert(itrc))
1829 END IF
1830 END DO
1831 END DO
1832 END DO
1833 END DO
1834# endif
1835
1836
1837
1838
1839
1840
1841
1842
1843
1844
1845 ic=0
1846
1848
1849
1850
1851
1853 ic=ic+1
1854 END IF
1855
1856
1857
1859 & lbi, ubi, lbj, ubj,
n(ng),
nt(ng), &
1860 & imins, imaxs, jmins, jmaxs, &
1861 & nstp, nnew, &
1862 & t)
1863
1864
1865
1868 DO j=jstrr,jendr
1869 DO i=istrr,iendr
1870 t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)+ &
1872 &
clima(ng)%Tnudgcof(i,j,k,ic)* &
1873 & (
clima(ng)%tclm(i,j,k,ic)- &
1874 & t(i,j,k,nnew,itrc))
1875 END DO
1876 END DO
1877 END DO
1878 END IF
1879
1880# ifdef MASKING
1881
1882
1883
1885 DO j=jstrr,jendr
1886 DO i=istrr,iendr
1887 t(i,j,k,nnew,itrc)=t(i,j,k,nnew,itrc)*rmask(i,j)
1888 END DO
1889 END DO
1890 END DO
1891# endif
1892# ifdef DIAGNOSTICS_TS
1893
1894
1895
1897 DO j=jstrr,jendr
1898 DO i=istrr,iendr
1899 diatwrk(i,j,k,itrc,
itrate)=t(i,j,k,nnew,itrc)- &
1900 & diatwrk(i,j,k,itrc,
itrate)
1901
1902
1903 END DO
1904 END DO
1905 END DO
1906# endif
1907
1908
1909
1912 & lbi, ubi, lbj, ubj, 1,
n(ng), &
1913 & t(:,:,:,nnew,itrc))
1914 END IF
1915 END DO
1916# ifdef DISTRIBUTE
1917
1918
1919
1921 & lbi, ubi, lbj, ubj, 1,
n(ng), 1,
nt(ng), &
1924 & t(:,:,:,nnew,:))
1925# endif
1926# if defined FLOATS && defined FLOAT_VWALK
1927
1928
1929
1930
1931
1932
1933 DO j=jstrr,jendr
1934 DO i=istrr,iendr
1936 daktdz(i,j,k)=(akt(i,j,k,1)-akt(i,j,k-1,1))/hz(i,j,k)
1937 END DO
1938 END DO
1939 END DO
1940
1941
1942
1945 & lbi, ubi, lbj, ubj, 1,
n(ng), &
1946 & daktdz)
1947 END IF
1948
1949# ifdef DISTRIBUTE
1951 & lbi, ubi, lbj, ubj, 1,
n(ng), &
1954 & daktdz)
1955# endif
1956# endif
1957
1958
1959
1960
1961
1962 IF (lmpdata) THEN
1963# ifdef DIAGNOSTICS_TS
1964 IF (allocated(dhadv)) deallocate (dhadv)
1965 IF (allocated(dvadv)) deallocate (dvadv)
1966# endif
1967 IF (allocated(ta)) deallocate (ta)
1968 IF (allocated(ua)) deallocate (ua)
1969 IF (allocated(va)) deallocate (va)
1970 IF (allocated(wa)) deallocate (wa)
1971 END IF
1972
1973 RETURN
subroutine mpdata_adiff_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, rmask, umask, vmask, rmask_wet, umask_wet, vmask_wet, pm, pn, omn, om_u, on_v, z_r, ohz, huon, hvom, w, t, ta, ua, va, wa)
subroutine exchange_r3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
type(t_clima), dimension(:), allocatable clima
type(t_bcp), dimension(:,:), allocatable bry_contact
type(t_ngc), dimension(:), allocatable rcontact
type(t_adv), dimension(:,:), allocatable hadvection
integer, dimension(:), allocatable n
type(t_bounds), dimension(:), allocatable bounds
type(t_domain), dimension(:), allocatable domain
integer, dimension(:), allocatable lm
type(t_adv), dimension(:,:), allocatable vadvection
integer, dimension(:), allocatable nt
integer, dimension(:), allocatable mm
logical, dimension(:), allocatable luvsrc
logical, dimension(:,:), allocatable ltracersrc
real(dp), dimension(:), allocatable dt
logical, dimension(:), allocatable ewperiodic
logical, dimension(:), allocatable nsperiodic
logical, dimension(:), allocatable lwsrc
logical, dimension(:,:), allocatable compositegrid
integer, parameter isouth
integer, dimension(:), pointer inert
logical, dimension(:), allocatable refinedgrid
logical, dimension(:,:), allocatable ltracerclm
integer, parameter inorth
logical, dimension(:,:), allocatable lnudgetclm
type(t_sources), dimension(:), allocatable sources
integer, dimension(:), allocatable nsrc
subroutine mp_exchange4d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, lbk, ubk, lbt, ubt, nghost, ew_periodic, ns_periodic, a, b, c)
subroutine mp_exchange3d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, lbk, ubk, nghost, ew_periodic, ns_periodic, a, b, c, d)
subroutine, public bry_fluxes(dg, rg, cr, model, tile, imins, imaxs, jmins, jmaxs, ilb, iub, jlb, jub, scale, fx, fe, f_west, f_east, f_south, f_north)
subroutine, public t3dbc_tile(ng, tile, itrc, ic, lbi, ubi, lbj, ubj, ubk, ubt, imins, imaxs, jmins, jmaxs, nstp, nout, t)