247
248
253#if defined SEDIMENT && defined SED_MORPH
255#endif
257
259#ifdef DISTRIBUTE
261#endif
266#ifdef WET_DRY
268#endif
269
270
271
272 integer, intent(in ) :: ng, tile
273 integer, intent(in ) :: LBi, UBi, LBj, UBj, UBk
274 integer, intent(in ) :: IminS, ImaxS, JminS, JmaxS
275 integer, intent(in ) :: krhs, kstp, knew
276#ifdef SOLVE3D
277 integer, intent(in ) :: nstp, nnew
278#endif
279
280#ifdef ASSUMED_SHAPE
281# ifdef MASKING
282 real(r8), intent(in ) :: pmask(LBi:,LBj:)
283 real(r8), intent(in ) :: rmask(LBi:,LBj:)
284 real(r8), intent(in ) :: umask(LBi:,LBj:)
285 real(r8), intent(in ) :: vmask(LBi:,LBj:)
286# endif
287 real(r8), intent(in ) :: fomn(LBi:,LBj:)
288# if defined SEDIMENT && defined SED_MORPH
289 real(r8), intent(inout) :: h(LBi:,LBj:)
290# else
291 real(r8), intent(in ) :: h(LBi:,LBj:)
292# endif
293 real(r8), intent(in ) :: om_u(LBi:,LBj:)
294 real(r8), intent(in ) :: om_v(LBi:,LBj:)
295 real(r8), intent(in ) :: on_u(LBi:,LBj:)
296 real(r8), intent(in ) :: on_v(LBi:,LBj:)
297 real(r8), intent(in ) :: omn(LBi:,LBj:)
298 real(r8), intent(in ) :: pm(LBi:,LBj:)
299 real(r8), intent(in ) :: pn(LBi:,LBj:)
300# if defined CURVGRID && defined UV_ADV
301 real(r8), intent(in ) :: dndx(LBi:,LBj:)
302 real(r8), intent(in ) :: dmde(LBi:,LBj:)
303# endif
304# if defined UV_VIS2 || defined UV_VIS4
305 real(r8), intent(in ) :: pmon_r(LBi:,LBj:)
306 real(r8), intent(in ) :: pnom_r(LBi:,LBj:)
307 real(r8), intent(in ) :: pmon_p(LBi:,LBj:)
308 real(r8), intent(in ) :: pnom_p(LBi:,LBj:)
309 real(r8), intent(in ) :: om_r(LBi:,LBj:)
310 real(r8), intent(in ) :: on_r(LBi:,LBj:)
311 real(r8), intent(in ) :: om_p(LBi:,LBj:)
312 real(r8), intent(in ) :: on_p(LBi:,LBj:)
313# ifdef UV_VIS2
314 real(r8), intent(in ) :: visc2_p(LBi:,LBj:)
315 real(r8), intent(in ) :: visc2_r(LBi:,LBj:)
316# endif
317# ifdef UV_VIS4
318 real(r8), intent(in ) :: visc4_p(LBi:,LBj:)
319 real(r8), intent(in ) :: visc4_r(LBi:,LBj:)
320# endif
321# endif
322# if defined SEDIMENT && defined SED_MORPH
323 real(r8), intent(in ) :: bed_thick(LBi:,LBj:,:)
324# endif
325# ifdef WEC
326# ifdef WEC_VF
327# ifdef WEC_ROLLER
328 real(r8), intent(in) :: rurol2d(LBi:,LBj:)
329 real(r8), intent(in) :: rvrol2d(LBi:,LBj:)
330# endif
331# ifdef BOTTOM_STREAMING
332 real(r8), intent(in) :: rubst2d(LBi:,LBj:)
333 real(r8), intent(in) :: rvbst2d(LBi:,LBj:)
334# endif
335# ifdef SURFACE_STREAMING
336 real(r8), intent(in) :: russt2d(LBi:,LBj:)
337 real(r8), intent(in) :: rvsst2d(LBi:,LBj:)
338# endif
339 real(r8), intent(in) :: rubrk2d(LBi:,LBj:)
340 real(r8), intent(in) :: rvbrk2d(LBi:,LBj:)
341 real(r8), intent(in) :: rukvf2d(LBi:,LBj:)
342 real(r8), intent(in) :: rvkvf2d(LBi:,LBj:)
343 real(r8), intent(in) :: bh(LBi:,LBj:)
344 real(r8), intent(in) :: qsp(LBi:,LBj:)
345 real(r8), intent(in) :: zetaw(LBi:,LBj:)
346# endif
347 real(r8), intent(in) :: ubar_stokes(LBi:,LBj:)
348 real(r8), intent(in) :: vbar_stokes(LBi:,LBj:)
349# endif
350# if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
351 real(r8), intent(in ) :: eq_tide(LBi:,LBj:)
352# endif
353# ifndef SOLVE3D
354 real(r8), intent(in ) :: sustr(LBi:,LBj:)
355 real(r8), intent(in ) :: svstr(LBi:,LBj:)
356 real(r8), intent(in ) :: bustr(LBi:,LBj:)
357 real(r8), intent(in ) :: bvstr(LBi:,LBj:)
358# ifdef ATM_PRESS
359 real(r8), intent(in ) :: Pair(LBi:,LBj:)
360# endif
361# else
362# ifdef VAR_RHO_2D
363 real(r8), intent(in ) :: rhoA(LBi:,LBj:)
364 real(r8), intent(in ) :: rhoS(LBi:,LBj:)
365# endif
366 real(r8), intent(inout) :: DU_avg1(LBi:,LBj:)
367 real(r8), intent(inout) :: DU_avg2(LBi:,LBj:)
368 real(r8), intent(inout) :: DV_avg1(LBi:,LBj:)
369 real(r8), intent(inout) :: DV_avg2(LBi:,LBj:)
370 real(r8), intent(inout) :: Zt_avg1(LBi:,LBj:)
371 real(r8), intent(inout) :: rufrc(LBi:,LBj:)
372 real(r8), intent(inout) :: rvfrc(LBi:,LBj:)
373 real(r8), intent(inout) :: ru(LBi:,LBj:,0:,:)
374 real(r8), intent(inout) :: rv(LBi:,LBj:,0:,:)
375# endif
376# ifdef WET_DRY
377 real(r8), intent(inout) :: pmask_full(LBi:,LBj:)
378 real(r8), intent(inout) :: rmask_full(LBi:,LBj:)
379 real(r8), intent(inout) :: umask_full(LBi:,LBj:)
380 real(r8), intent(inout) :: vmask_full(LBi:,LBj:)
381
382 real(r8), intent(inout) :: pmask_wet(LBi:,LBj:)
383 real(r8), intent(inout) :: rmask_wet(LBi:,LBj:)
384 real(r8), intent(inout) :: umask_wet(LBi:,LBj:)
385 real(r8), intent(inout) :: vmask_wet(LBi:,LBj:)
386# ifdef SOLVE3D
387 real(r8), intent(inout) :: rmask_wet_avg(LBi:,LBj:)
388# endif
389# endif
390# ifdef DIAGNOSTICS_UV
391 real(r8), intent(inout) :: DiaU2wrk(LBi:,LBj:,:)
392 real(r8), intent(inout) :: DiaV2wrk(LBi:,LBj:,:)
393 real(r8), intent(inout) :: DiaRUbar(LBi:,LBj:,:,:)
394 real(r8), intent(inout) :: DiaRVbar(LBi:,LBj:,:,:)
395# ifdef SOLVE3D
396 real(r8), intent(inout) :: DiaU2int(LBi:,LBj:,:)
397 real(r8), intent(inout) :: DiaV2int(LBi:,LBj:,:)
398 real(r8), intent(inout) :: DiaRUfrc(LBi:,LBj:,:,:)
399 real(r8), intent(inout) :: DiaRVfrc(LBi:,LBj:,:,:)
400# endif
401# endif
402 real(r8), intent(inout) :: rubar(LBi:,LBj:,:)
403 real(r8), intent(inout) :: rvbar(LBi:,LBj:,:)
404 real(r8), intent(inout) :: rzeta(LBi:,LBj:,:)
405 real(r8), intent(inout) :: ubar(LBi:,LBj:,:)
406 real(r8), intent(inout) :: vbar(LBi:,LBj:,:)
407 real(r8), intent(inout) :: zeta(LBi:,LBj:,:)
408# if defined NESTING && !defined SOLVE3D
409 real(r8), intent(out ) :: DU_flux(LBi:,LBj:)
410 real(r8), intent(out ) :: DV_flux(LBi:,LBj:)
411# endif
412
413#else
414
415# ifdef MASKING
416 real(r8), intent(in ) :: pmask(LBi:UBi,LBj:UBj)
417 real(r8), intent(in ) :: rmask(LBi:UBi,LBj:UBj)
418 real(r8), intent(in ) :: umask(LBi:UBi,LBj:UBj)
419 real(r8), intent(in ) :: vmask(LBi:UBi,LBj:UBj)
420# endif
421 real(r8), intent(in ) :: fomn(LBi:UBi,LBj:UBj)
422# if defined SEDIMENT && defined SED_MORPH
423 real(r8), intent(inout) :: h(LBi:UBi,LBj:UBj)
424# else
425 real(r8), intent(in ) :: h(LBi:UBi,LBj:UBj)
426# endif
427 real(r8), intent(in ) :: om_u(LBi:UBi,LBj:UBj)
428 real(r8), intent(in ) :: om_v(LBi:UBi,LBj:UBj)
429 real(r8), intent(in ) :: on_u(LBi:UBi,LBj:UBj)
430 real(r8), intent(in ) :: on_v(LBi:UBi,LBj:UBj)
431 real(r8), intent(in ) :: omn(LBi:UBi,LBj:UBj)
432 real(r8), intent(in ) :: pm(LBi:UBi,LBj:UBj)
433 real(r8), intent(in ) :: pn(LBi:UBi,LBj:UBj)
434# if defined CURVGRID && defined UV_ADV
435 real(r8), intent(in ) :: dndx(LBi:UBi,LBj:UBj)
436 real(r8), intent(in ) :: dmde(LBi:UBi,LBj:UBj)
437# endif
438# if defined UV_VIS2 || defined UV_VIS4
439 real(r8), intent(in ) :: pmon_r(LBi:UBi,LBj:UBj)
440 real(r8), intent(in ) :: pnom_r(LBi:UBi,LBj:UBj)
441 real(r8), intent(in ) :: pmon_p(LBi:UBi,LBj:UBj)
442 real(r8), intent(in ) :: pnom_p(LBi:UBi,LBj:UBj)
443 real(r8), intent(in ) :: om_r(LBi:UBi,LBj:UBj)
444 real(r8), intent(in ) :: on_r(LBi:UBi,LBj:UBj)
445 real(r8), intent(in ) :: om_p(LBi:UBi,LBj:UBj)
446 real(r8), intent(in ) :: on_p(LBi:UBi,LBj:UBj)
447# ifdef UV_VIS2
448 real(r8), intent(in ) :: visc2_p(LBi:UBi,LBj:UBj)
449 real(r8), intent(in ) :: visc2_r(LBi:UBi,LBj:UBj)
450# endif
451# ifdef UV_VIS4
452 real(r8), intent(in ) :: visc4_p(LBi:UBi,LBj:UBj)
453 real(r8), intent(in ) :: visc4_r(LBi:UBi,LBj:UBj)
454# endif
455# endif
456# if defined SEDIMENT && defined SED_MORPH
457 real(r8), intent(in ) :: bed_thick(LBi:UBi,LBj:UBj,1:3)
458# endif
459# ifdef WEC
460# ifdef WEC_VF
461# ifdef WEC_ROLLER
462 real(r8), intent(in) :: rurol2d(LBi:UBi,LBj:UBj)
463 real(r8), intent(in) :: rvrol2d(LBi:UBi,LBj:UBj)
464# endif
465# ifdef BOTTOM_STREAMING
466 real(r8), intent(in) :: rubst2d(LBi:UBi,LBj:UBj)
467 real(r8), intent(in) :: rvbst2d(LBi:UBi,LBj:UBj)
468# endif
469# ifdef SURFACE_STREAMING
470 real(r8), intent(in) :: russt2d(LBi:UBi,LBj:UBj)
471 real(r8), intent(in) :: rvsst2d(LBi:UBi,LBj:UBj)
472# endif
473 real(r8), intent(in) :: rubrk2d(LBi:UBi,LBj:UBj)
474 real(r8), intent(in) :: rvbrk2d(LBi:UBi,LBj:UBj)
475 real(r8), intent(in) :: rukvf2d(LBi:UBi,LBj:UBj)
476 real(r8), intent(in) :: rvkvf2d(LBi:UBi,LBj:UBj)
477 real(r8), intent(in) :: bh(LBi:UBi,LBj:UBj)
478 real(r8), intent(in) :: qsp(LBi:UBi,LBj:UBj)
479 real(r8), intent(in) :: zetaw(LBi:UBi,LBj:UBj)
480# endif
481 real(r8), intent(in) :: ubar_stokes(LBi:UBi,LBj:UBj)
482 real(r8), intent(in) :: vbar_stokes(LBi:UBi,LBj:UBj)
483# endif
484# if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
485 real(r8), intent(in ) :: eq_tide(LBi:UBi,LBj:UBj)
486# endif
487# ifndef SOLVE3D
488 real(r8), intent(in ) :: sustr(LBi:UBi,LBj:UBj)
489 real(r8), intent(in ) :: svstr(LBi:UBi,LBj:UBj)
490 real(r8), intent(in ) :: bustr(LBi:UBi,LBj:UBj)
491 real(r8), intent(in ) :: bvstr(LBi:UBi,LBj:UBj)
492# ifdef ATM_PRESS
493 real(r8), intent(in ) :: Pair(LBi:UBi,LBj:UBj)
494# endif
495# else
496# ifdef VAR_RHO_2D
497 real(r8), intent(in ) :: rhoA(LBi:UBi,LBj:UBj)
498 real(r8), intent(in ) :: rhoS(LBi:UBi,LBj:UBj)
499# endif
500 real(r8), intent(inout) :: DU_avg1(LBi:UBi,LBj:UBj)
501 real(r8), intent(inout) :: DU_avg2(LBi:UBi,LBj:UBj)
502 real(r8), intent(inout) :: DV_avg1(LBi:UBi,LBj:UBj)
503 real(r8), intent(inout) :: DV_avg2(LBi:UBi,LBj:UBj)
504 real(r8), intent(inout) :: Zt_avg1(LBi:UBi,LBj:UBj)
505 real(r8), intent(inout) :: rufrc(LBi:UBi,LBj:UBj)
506 real(r8), intent(inout) :: rvfrc(LBi:UBi,LBj:UBj)
507 real(r8), intent(inout) :: ru(LBi:UBi,LBj:UBj,0:UBk,2)
508 real(r8), intent(inout) :: rv(LBi:UBi,LBj:UBj,0:UBk,2)
509# endif
510# ifdef WET_DRY
511 real(r8), intent(inout) :: pmask_full(LBi:UBi,LBj:UBj)
512 real(r8), intent(inout) :: rmask_full(LBi:UBi,LBj:UBj)
513 real(r8), intent(inout) :: umask_full(LBi:UBi,LBj:UBj)
514 real(r8), intent(inout) :: vmask_full(LBi:UBi,LBj:UBj)
515
516 real(r8), intent(inout) :: pmask_wet(LBi:UBi,LBj:UBj)
517 real(r8), intent(inout) :: rmask_wet(LBi:UBi,LBj:UBj)
518 real(r8), intent(inout) :: umask_wet(LBi:UBi,LBj:UBj)
519 real(r8), intent(inout) :: vmask_wet(LBi:UBi,LBj:UBj)
520# ifdef SOLVE3D
521 real(r8), intent(inout) :: rmask_wet_avg(LBi:UBi,LBj:UBj)
522# endif
523# endif
524# ifdef DIAGNOSTICS_UV
525 real(r8), intent(inout) :: DiaU2wrk(LBi:UBi,LBj:UBj,NDM2d)
526 real(r8), intent(inout) :: DiaV2wrk(LBi:UBi,LBj:UBj,NDM2d)
527 real(r8), intent(inout) :: DiaRUbar(LBi:UBi,LBj:UBj,2,NDM2d-1)
528 real(r8), intent(inout) :: DiaRVbar(LBi:UBi,LBj:UBj,2,NDM2d-1)
529# ifdef SOLVE3D
530 real(r8), intent(inout) :: DiaU2int(LBi:UBi,LBj:UBj,NDM2d)
531 real(r8), intent(inout) :: DiaV2int(LBi:UBi,LBj:UBj,NDM2d)
532 real(r8), intent(inout) :: DiaRUfrc(LBi:UBi,LBj:UBj,3,NDM2d-1)
533 real(r8), intent(inout) :: DiaRVfrc(LBi:UBi,LBj:UBj,3,NDM2d-1)
534# endif
535# endif
536 real(r8), intent(inout) :: rubar(LBi:UBi,LBj:UBj,2)
537 real(r8), intent(inout) :: rvbar(LBi:UBi,LBj:UBj,2)
538 real(r8), intent(inout) :: rzeta(LBi:UBi,LBj:UBj,2)
539 real(r8), intent(inout) :: ubar(LBi:UBi,LBj:UBj,:)
540 real(r8), intent(inout) :: vbar(LBi:UBi,LBj:UBj,:)
541 real(r8), intent(inout) :: zeta(LBi:UBi,LBj:UBj,:)
542# if defined NESTING && !defined SOLVE3D
543 real(r8), intent(out ) :: DU_flux(LBi:UBi,LBj:UBj)
544 real(r8), intent(out ) :: DV_flux(LBi:UBi,LBj:UBj)
545# endif
546#endif
547
548
549
550 logical :: CORRECTOR_2D_STEP
551
552 integer :: i, is, j, ptsk
553#ifdef DIAGNOSTICS_UV
554 integer :: idiag
555#endif
556
557 real(r8) :: cff, cff1, cff2, cff3, cff4, cff5, cff6, cff7
558 real(r8) :: fac, fac1, fac2, fac3
559
560 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Dgrad
561 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Dnew
562 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Drhs
563 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Drhs_p
564 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Dstp
565 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: DUon
566 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: DVom
567#ifdef WEC
568 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: DUSon
569 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: DVSom
570#endif
571#ifdef WEC_VF
572 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: DUSom
573 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: DVSon
574#endif
575#ifdef UV_VIS4
576 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: LapU
577 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: LapV
578#endif
579 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: UFe
580 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: UFx
581 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: VFe
582 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: VFx
583 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: grad
584 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: gzeta
585 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: gzeta2
586#if defined VAR_RHO_2D && defined SOLVE3D
587 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: gzetaSA
588#endif
589 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rhs_ubar
590 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rhs_vbar
591 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rhs_zeta
592 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: zeta_new
593 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: zwrk
594#ifdef WET_DRY
595 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: wetdry
596#endif
597#ifdef DIAGNOSTICS_UV
598 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Uwrk
599 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Vwrk
600 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,NDM2d-1) :: DiaU2rhs
601 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,NDM2d-1) :: DiaV2rhs
602#endif
603
604#include "set_bounds.h"
605
606 ptsk=3-kstp
608
609
610
611
612
613#if defined DISTRIBUTE && !defined NESTING
614
615
616
617
618
619
620
621
622 DO j=jstrv-2,jendp2
623 DO i=istru-2,iendp2
624 drhs(i,j)=zeta(i,j,krhs)+h(i,j)
625 END DO
626 END DO
627 DO j=jstrv-2,jendp2
628 DO i=istru-1,iendp2
629 cff=0.5_r8*on_u(i,j)
630 cff1=cff*(drhs(i,j)+drhs(i-1,j))
631 duon(i,j)=ubar(i,j,krhs)*cff1
632# ifdef WEC
633# ifdef WET_DRY
634 cff5=abs(abs(umask_wet(i,j))-1.0_r8)
635 cff6=0.5_r8+dsign(0.5_r8,ubar_stokes(i,j))*umask_wet(i,j)
636 cff7=0.5_r8*umask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
637 cff1=cff1*cff7
638# endif
639 duson(i,j)=ubar_stokes(i,j)*cff1
640 duon(i,j)=duon(i,j)+duson(i,j)
641# endif
642 END DO
643 END DO
644 DO j=jstrv-1,jendp2
645 DO i=istru-2,iendp2
646 cff=0.5_r8*om_v(i,j)
647 cff1=cff*(drhs(i,j)+drhs(i,j-1))
648 dvom(i,j)=vbar(i,j,krhs)*cff1
649# ifdef WEC
650# ifdef WET_DRY
651 cff5=abs(abs(vmask_wet(i,j))-1.0_r8)
652 cff6=0.5_r8+dsign(0.5_r8,vbar_stokes(i,j))*vmask_wet(i,j)
653 cff7=0.5_r8*vmask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
654 cff1=cff1*cff7
655# endif
656 dvsom(i,j)=vbar_stokes(i,j)*cff1
657 dvom(i,j)=dvom(i,j)+dvsom(i,j)
658# endif
659 END DO
660 END DO
661
662#else
663
664 DO j=jstrvm2-1,jendp2
665 DO i=istrum2-1,iendp2
666 drhs(i,j)=zeta(i,j,krhs)+h(i,j)
667 END DO
668 END DO
669 DO j=jstrvm2-1,jendp2
670 DO i=istrum2,iendp2
671 cff=0.5_r8*on_u(i,j)
672 cff1=cff*(drhs(i,j)+drhs(i-1,j))
673 duon(i,j)=ubar(i,j,krhs)*cff1
674# ifdef WEC
675# ifdef WET_DRY
676 cff5=abs(abs(umask_wet(i,j))-1.0_r8)
677 cff6=0.5_r8+dsign(0.5_r8,ubar_stokes(i,j))*umask_wet(i,j)
678 cff7=0.5_r8*umask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
679 cff1=cff1*cff7
680# endif
681 duson(i,j)=ubar_stokes(i,j)*cff1
682 duon(i,j)=duon(i,j)+duson(i,j)
683# endif
684 END DO
685 END DO
686 DO j=jstrvm2,jendp2
687 DO i=istrum2-1,iendp2
688 cff=0.5_r8*om_v(i,j)
689 cff1=cff*(drhs(i,j)+drhs(i,j-1))
690 dvom(i,j)=vbar(i,j,krhs)*cff1
691# ifdef WEC
692# ifdef WET_DRY
693 cff5=abs(abs(vmask_wet(i,j))-1.0_r8)
694 cff6=0.5_r8+dsign(0.5_r8,vbar_stokes(i,j))*vmask_wet(i,j)
695 cff7=0.5_r8*vmask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
696 cff1=cff1*cff7
697# endif
698 dvsom(i,j)=vbar_stokes(i,j)*cff1
699 dvom(i,j)=dvom(i,j)+dvsom(i,j)
700# endif
701 END DO
702 END DO
703#endif
704#ifdef DISTRIBUTE
705
708 & imins, imaxs, jmins, jmaxs, &
709 & duon)
711 & imins, imaxs, jmins, jmaxs, &
712 & dvom)
713 END IF
715 & imins, imaxs, jmins, jmaxs, &
718 & duon, dvom)
719#endif
720
721
722
723
726 & lbi, ubi, lbj, ubj, &
727 & imins, imaxs, jmins, jmaxs, &
728 & krhs, &
729#ifdef MASKING
730 & umask, vmask, &
731#endif
732 & om_v, on_u, &
733 & ubar, vbar, &
734 & drhs, duon, dvom)
735 END IF
736#ifdef SOLVE3D
737
738
739
740
741
743 IF (first_2d_step) THEN
744
745
746
747 cff2=(-1.0_r8/12.0_r8)*
weight(2,
iif(ng)+1,ng)
748 DO j=jstrr,jendr
749 DO i=istrr,iendr
750 zt_avg1(i,j)=0.0_r8
751 END DO
752 DO i=istr,iendr
753 du_avg1(i,j)=0.0_r8
754 du_avg2(i,j)=cff2*duon(i,j)
755 END DO
756 END DO
757 DO j=jstr,jendr
758 DO i=istrr,iendr
759 dv_avg1(i,j)=0.0_r8
760 dv_avg2(i,j)=cff2*dvom(i,j)
761 END DO
762 END DO
763 ELSE
764
765
766
767
768
770 cff2=(8.0_r8/12.0_r8)*
weight(2,
iif(ng) ,ng)- &
771 & (1.0_r8/12.0_r8)*
weight(2,
iif(ng)+1,ng)
772 DO j=jstrr,jendr
773 DO i=istrr,iendr
774 zt_avg1(i,j)=zt_avg1(i,j)+cff1*zeta(i,j,krhs)
775 END DO
776 DO i=istr,iendr
777 du_avg1(i,j)=du_avg1(i,j)+cff1*duon(i,j)
778# ifdef WEC
779 du_avg1(i,j)=du_avg1(i,j)-cff1*duson(i,j)
780# endif
781 du_avg2(i,j)=du_avg2(i,j)+cff2*duon(i,j)
782 END DO
783 END DO
784 DO j=jstr,jendr
785 DO i=istrr,iendr
786 dv_avg1(i,j)=dv_avg1(i,j)+cff1*dvom(i,j)
787# ifdef WEC
788 dv_avg1(i,j)=dv_avg1(i,j)-cff1*dvsom(i,j)
789# endif
790 dv_avg2(i,j)=dv_avg2(i,j)+cff2*dvom(i,j)
791 END DO
792 END DO
793 END IF
794 ELSE
795 IF (first_2d_step) THEN
797 ELSE
798 cff2=(5.0_r8/12.0_r8)*
weight(2,
iif(ng),ng)
799 END IF
800 DO j=jstrr,jendr
801 DO i=istr,iendr
802 du_avg2(i,j)=du_avg2(i,j)+cff2*duon(i,j)
803 END DO
804 END DO
805 DO j=jstr,jendr
806 DO i=istrr,iendr
807 dv_avg2(i,j)=dv_avg2(i,j)+cff2*dvom(i,j)
808 END DO
809 END DO
810 END IF
811
812
813
814# ifdef NESTING
815
816
817
818
819# endif
820
824 & lbi, ubi, lbj, ubj, &
825 & zt_avg1)
827 & lbi, ubi, lbj, ubj, &
828 & du_avg1)
830 & lbi, ubi, lbj, ubj, &
831 & dv_avg1)
832# ifdef NESTING
834 & lbi, ubi, lbj, ubj, &
835 & du_avg2)
837 & lbi, ubi, lbj, ubj, &
838 & dv_avg2)
839# endif
840 END IF
841# ifdef DISTRIBUTE
843 & lbi, ubi, lbj, ubj, &
846 & zt_avg1, du_avg1, dv_avg1)
847# ifdef NESTING
849 & lbi, ubi, lbj, ubj, &
852 & du_avg2, dv_avg2)
853# endif
854# endif
855 END IF
856#endif
857#ifdef WET_DRY
858
859
860
861
862
864 & lbi, ubi, lbj, ubj, &
865 & imins, imaxs, jmins, jmaxs, &
866# ifdef MASKING
867 & pmask, rmask, umask, vmask, &
868# endif
869 & h, zeta(:,:,kstp), &
870# ifdef SOLVE3D
871 & du_avg1, dv_avg1, &
872 & rmask_wet_avg, &
873# endif
874 & pmask_wet, pmask_full, &
875 & rmask_wet, rmask_full, &
876 & umask_wet, umask_full, &
877 & vmask_wet, vmask_full)
878#endif
879
880
881
882
884
885
886
887
888
889
890
891
892#if defined VAR_RHO_2D && defined SOLVE3D
893
894
895
897#endif
898
899 IF (first_2d_step) THEN
901 DO j=jstrv-1,jend
902 DO i=istru-1,iend
903 rhs_zeta(i,j)=(duon(i,j)-duon(i+1,j))+ &
904 & (dvom(i,j)-dvom(i,j+1))
905 zeta_new(i,j)=zeta(i,j,kstp)+ &
906 & pm(i,j)*pn(i,j)*cff1*rhs_zeta(i,j)
907#ifdef MASKING
908 zeta_new(i,j)=zeta_new(i,j)*rmask(i,j)
909#endif
910 dnew(i,j)=zeta_new(i,j)+h(i,j)
911
912 zwrk(i,j)=0.5_r8*(zeta(i,j,kstp)+zeta_new(i,j))
913#if defined VAR_RHO_2D && defined SOLVE3D
914 gzeta(i,j)=(fac+rhos(i,j))*zwrk(i,j)
915 gzeta2(i,j)=gzeta(i,j)*zwrk(i,j)
916 gzetasa(i,j)=zwrk(i,j)*(rhos(i,j)-rhoa(i,j))
917#else
918 gzeta(i,j)=zwrk(i,j)
919 gzeta2(i,j)=zwrk(i,j)*zwrk(i,j)
920#endif
921 END DO
922 END DO
925 cff4=4.0_r8/25.0_r8
926 cff5=1.0_r8-2.0_r8*cff4
927 DO j=jstrv-1,jend
928 DO i=istru-1,iend
929 rhs_zeta(i,j)=(duon(i,j)-duon(i+1,j))+ &
930 & (dvom(i,j)-dvom(i,j+1))
931 zeta_new(i,j)=zeta(i,j,kstp)+ &
932 & pm(i,j)*pn(i,j)*cff1*rhs_zeta(i,j)
933#ifdef MASKING
934 zeta_new(i,j)=zeta_new(i,j)*rmask(i,j)
935#endif
936 dnew(i,j)=zeta_new(i,j)+h(i,j)
937
938 zwrk(i,j)=cff5*zeta(i,j,krhs)+ &
939 & cff4*(zeta(i,j,kstp)+zeta_new(i,j))
940#if defined VAR_RHO_2D && defined SOLVE3D
941 gzeta(i,j)=(fac+rhos(i,j))*zwrk(i,j)
942 gzeta2(i,j)=gzeta(i,j)*zwrk(i,j)
943 gzetasa(i,j)=zwrk(i,j)*(rhos(i,j)-rhoa(i,j))
944#else
945 gzeta(i,j)=zwrk(i,j)
946 gzeta2(i,j)=zwrk(i,j)*zwrk(i,j)
947#endif
948 END DO
949 END DO
950 ELSE IF (corrector_2d_step) THEN
951 cff1=
dtfast(ng)*5.0_r8/12.0_r8
952 cff2=
dtfast(ng)*8.0_r8/12.0_r8
953 cff3=
dtfast(ng)*1.0_r8/12.0_r8
954 cff4=2.0_r8/5.0_r8
955 cff5=1.0_r8-cff4
956 DO j=jstrv-1,jend
957 DO i=istru-1,iend
958 cff=cff1*((duon(i,j)-duon(i+1,j))+ &
959 & (dvom(i,j)-dvom(i,j+1)))
960 zeta_new(i,j)=zeta(i,j,kstp)+ &
961 & pm(i,j)*pn(i,j)*(cff+ &
962 & cff2*rzeta(i,j,kstp)- &
963 & cff3*rzeta(i,j,ptsk))
964#ifdef MASKING
965 zeta_new(i,j)=zeta_new(i,j)*rmask(i,j)
966#endif
967 dnew(i,j)=zeta_new(i,j)+h(i,j)
968
969 zwrk(i,j)=cff5*zeta_new(i,j)+cff4*zeta(i,j,krhs)
970#if defined VAR_RHO_2D && defined SOLVE3D
971 gzeta(i,j)=(fac+rhos(i,j))*zwrk(i,j)
972 gzeta2(i,j)=gzeta(i,j)*zwrk(i,j)
973 gzetasa(i,j)=zwrk(i,j)*(rhos(i,j)-rhoa(i,j))
974#else
975 gzeta(i,j)=zwrk(i,j)
976 gzeta2(i,j)=zwrk(i,j)*zwrk(i,j)
977#endif
978 END DO
979 END DO
980 END IF
981
982
983
984#ifdef WET_DRY
985
986
987#endif
988
989 DO j=jstr,jend
990 DO i=istr,iend
991 zeta(i,j,knew)=zeta_new(i,j)
992#if defined WET_DRY && defined MASKING
993 zeta(i,j,knew)=zeta(i,j,knew)+ &
994 & (
dcrit(ng)-h(i,j))*(1.0_r8-rmask(i,j))
995#endif
996 END DO
997 END DO
998
999
1000
1002 DO j=jstr,jend
1003 DO i=istr,iend
1004 rzeta(i,j,krhs)=rhs_zeta(i,j)
1005 END DO
1006 END DO
1009 & lbi, ubi, lbj, ubj, &
1010 & rzeta(:,:,krhs))
1011 END IF
1012#ifdef DISTRIBUTE
1014 & lbi, ubi, lbj, ubj, &
1017 & rzeta(:,:,krhs))
1018#endif
1019 END IF
1020
1021
1022
1023
1024
1027 IF (int(
sources(ng)%Dsrc(is)).eq.2)
THEN
1030 IF (((istrr.le.i).and.(i.le.iendr)).and. &
1031 & ((jstrr.le.j).and.(j.le.jendr))) THEN
1032 zeta(i,j,knew)=zeta(i,j,knew)+ &
1034 & pm(i,j)*pn(i,j)*
dtfast(ng)
1035 END IF
1036 END IF
1037 END DO
1038 END IF
1039
1040#if defined SEDIMENT && defined SED_MORPH
1041
1042
1043
1044
1045
1047 DO j=jstr,jend
1048 DO i=istr,iend
1049 cff=fac*(bed_thick(i,j,nstp)-bed_thick(i,j,nnew))
1050 h(i,j)=h(i,j)-cff
1051 END DO
1052 END DO
1053#endif
1054
1055
1056
1058 & lbi, ubi, lbj, ubj, &
1059 & imins, imaxs, jmins, jmaxs, &
1060 & krhs, kstp, knew, &
1061 & zeta)
1064 & lbi, ubi, lbj, ubj, &
1065 & zeta(:,:,knew))
1066 END IF
1067#ifdef DISTRIBUTE
1069 & lbi, ubi, lbj, ubj, &
1072 & zeta(:,:,knew))
1073#endif
1074
1075
1076
1077
1078
1079
1080
1081
1082
1084 cff2=1.0_r8/3.0_r8
1085#if !defined SOLVE3D && defined ATM_PRESS
1086 fac3=0.5_r8*100.0_r8/
rho0
1087#endif
1088 DO j=jstr,jend
1089 DO i=istru,iend
1090 rhs_ubar(i,j)=cff1*on_u(i,j)* &
1091 & ((h(i-1,j)+ &
1092 & h(i ,j))* &
1093 & (gzeta(i-1,j)- &
1094 & gzeta(i ,j))+ &
1095#if defined VAR_RHO_2D && defined SOLVE3D
1096 & (h(i-1,j)- &
1097 & h(i ,j))* &
1098 & (gzetasa(i-1,j)+ &
1099 & gzetasa(i ,j)+ &
1100 & cff2*(rhoa(i-1,j)- &
1101 & rhoa(i ,j))* &
1102 & (zwrk(i-1,j)- &
1103 & zwrk(i ,j)))+ &
1104#endif
1105 & (gzeta2(i-1,j)- &
1106 & gzeta2(i ,j)))
1107#if defined ATM_PRESS && !defined SOLVE3D
1108 rhs_ubar(i,j)=rhs_ubar(i,j)- &
1109 & fac3*on_u(i,j)* &
1110 & (h(i-1,j)+h(i,j)+ &
1111 & gzeta(i-1,j)+gzeta(i,j))* &
1112 & (pair(i,j)-pair(i-1,j))
1113#endif
1114#if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
1115 rhs_ubar(i,j)=rhs_ubar(i,j)- &
1116 & cff1*on_u(i,j)* &
1117 & (h(i-1,j)+h(i,j)+ &
1118 & gzeta(i-1,j)+gzeta(i,j))* &
1119 & (eq_tide(i,j)-eq_tide(i-1,j))
1120#endif
1121#ifdef DIAGNOSTICS_UV
1122 diau2rhs(i,j,
m2pgrd)=rhs_ubar(i,j)
1123#endif
1124#if defined WEC_VF
1125 cff3=0.5_r8*on_u(i,j)* &
1126 & (h(i-1,j)+h(i,j)+ &
1127 & gzeta(i-1,j)+gzeta(i,j))
1128 cff4=cff3*
g*(zetaw(i-1,j)-zetaw(i,j))
1129 cff5=cff3*
g*(qsp(i-1,j)-qsp(i,j))
1130 cff6=cff3*(bh(i-1,j)-bh(i,j))
1131 cff7=rukvf2d(i,j)
1132 rhs_ubar(i,j)=rhs_ubar(i,j)-cff4-cff5+cff6+cff7
1133# ifdef DIAGNOSTICS_UV
1134 diau2rhs(i,j,m2zeta)=diau2rhs(i,j,
m2pgrd)
1135 diau2rhs(i,j,
m2pgrd)=diau2rhs(i,j,
m2pgrd)-cff4-cff5+cff6
1136 diau2rhs(i,j,m2zetw)=-cff4
1137 diau2rhs(i,j,m2zqsp)=-cff5
1138 diau2rhs(i,j,m2zbeh)=cff6
1139 diau2rhs(i,j,m2kvrf)=cff7
1140# ifndef UV_ADV
1141 diau2rhs(i,j,m2hjvf)=0.0_r8
1142# endif
1143# endif
1144#endif
1145 END DO
1146 IF (j.ge.jstrv) THEN
1147 DO i=istr,iend
1148 rhs_vbar(i,j)=cff1*om_v(i,j)* &
1149 & ((h(i,j-1)+ &
1150 & h(i,j ))* &
1151 & (gzeta(i,j-1)- &
1152 & gzeta(i,j ))+ &
1153#if defined VAR_RHO_2D && defined SOLVE3D
1154 & (h(i,j-1)- &
1155 & h(i,j ))* &
1156 & (gzetasa(i,j-1)+ &
1157 & gzetasa(i,j )+ &
1158 & cff2*(rhoa(i,j-1)- &
1159 & rhoa(i,j ))* &
1160 & (zwrk(i,j-1)- &
1161 & zwrk(i,j )))+ &
1162#endif
1163 & (gzeta2(i,j-1)- &
1164 & gzeta2(i,j )))
1165#if defined ATM_PRESS && !defined SOLVE3D
1166 rhs_vbar(i,j)=rhs_vbar(i,j)- &
1167 & fac3*om_v(i,j)* &
1168 & (h(i,j-1)+h(i,j)+ &
1169 & gzeta(i,j-1)+gzeta(i,j))* &
1170 & (pair(i,j)-pair(i,j-1))
1171#endif
1172#if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
1173 rhs_vbar(i,j)=rhs_vbar(i,j)- &
1174 & cff1*om_v(i,j)* &
1175 & (h(i,j-1)+h(i,j)+ &
1176 & gzeta(i,j-1)+gzeta(i,j))* &
1177 & (eq_tide(i,j)-eq_tide(i,j-1))
1178#endif
1179#ifdef DIAGNOSTICS_UV
1180 diav2rhs(i,j,
m2pgrd)=rhs_vbar(i,j)
1181#endif
1182#if defined WEC_VF
1183 cff3=0.5_r8*om_v(i,j)* &
1184 & (h(i,j-1)+h(i,j)+ &
1185 & gzeta(i,j-1)+gzeta(i,j))
1186 cff4=cff3*
g*(zetaw(i,j-1)-zetaw(i,j))
1187 cff5=cff3*
g*(qsp(i,j-1)-qsp(i,j))
1188 cff6=cff3*(bh(i,j-1)-bh(i,j))
1189 cff7=rvkvf2d(i,j)
1190 rhs_vbar(i,j)=rhs_vbar(i,j)-cff4-cff5+cff6+cff7
1191# ifdef DIAGNOSTICS_UV
1192 diav2rhs(i,j,m2zeta)=diav2rhs(i,j,
m2pgrd)
1193 diav2rhs(i,j,
m2pgrd)=diav2rhs(i,j,
m2pgrd)-cff4-cff5+cff6
1194 diav2rhs(i,j,m2zetw)=-cff4
1195 diav2rhs(i,j,m2zqsp)=-cff5
1196 diav2rhs(i,j,m2zbeh)=cff6
1197 diav2rhs(i,j,m2kvrf)=cff7
1198# ifndef UV_ADV
1199 diav2rhs(i,j,m2hjvf)=0.0_r8
1200# endif
1201# endif
1202#endif
1203 END DO
1204 END IF
1205 END DO
1206#ifdef UV_ADV
1207
1208
1209
1210
1211
1212# ifdef UV_C2ADVECTION
1213
1214
1215
1216 DO j=jstr,jend
1217 DO i=istru-1,iend
1218 ufx(i,j)=0.25_r8*(duon(i,j)+duon(i+1,j))* &
1219 & (ubar(i ,j,krhs)+ &
1220 & ubar(i+1,j,krhs))
1221 END DO
1222 END DO
1223
1224 DO j=jstr,jend+1
1225 DO i=istru,iend
1226 ufe(i,j)=0.25_r8*(dvom(i,j)+dvom(i-1,j))* &
1227 & (ubar(i,j ,krhs)+ &
1228 & ubar(i,j-1,krhs))
1229 END DO
1230 END DO
1231
1232 DO j=jstrv,jend
1233 DO i=istr,iend+1
1234 vfx(i,j)=0.25_r8*(duon(i,j)+duon(i,j-1))* &
1235 & (vbar(i ,j,krhs)+ &
1236 & vbar(i-1,j,krhs))
1237 END DO
1238 END DO
1239
1240 DO j=jstrv-1,jend
1241 DO i=istr,iend
1242 vfe(i,j)=0.25_r8*(dvom(i,j)+dvom(i,j+1))* &
1243 & (vbar(i,j ,krhs)+ &
1244 & vbar(i,j+1,krhs))
1245 END DO
1246 END DO
1247# else
1248
1249
1250
1251 DO j=jstr,jend
1252 DO i=istrum1,iendp1
1253 grad(i,j)=ubar(i-1,j,krhs)-2.0_r8*ubar(i,j,krhs)+ &
1254 & ubar(i+1,j,krhs)
1255 dgrad(i,j)=duon(i-1,j)-2.0_r8*duon(i,j)+duon(i+1,j)
1256 END DO
1257 END DO
1259 IF (
domain(ng)%Western_Edge(tile))
THEN
1260 DO j=jstr,jend
1261 grad(istr,j)=grad(istr+1,j)
1262 dgrad(istr,j)=dgrad(istr+1,j)
1263 END DO
1264 END IF
1265 END IF
1267 IF (
domain(ng)%Eastern_Edge(tile))
THEN
1268 DO j=jstr,jend
1269 grad(iend+1,j)=grad(iend,j)
1270 dgrad(iend+1,j)=dgrad(iend,j)
1271 END DO
1272 END IF
1273 END IF
1274
1275 cff=1.0_r8/6.0_r8
1276 DO j=jstr,jend
1277 DO i=istru-1,iend
1278 ufx(i,j)=0.25_r8*(ubar(i ,j,krhs)+ &
1279 & ubar(i+1,j,krhs)- &
1280 & cff*(grad(i,j)+grad(i+1,j)))* &
1281 & (duon(i,j)+duon(i+1,j)- &
1282 & cff*(dgrad(i,j)+dgrad(i+1,j)))
1283 END DO
1284 END DO
1285
1286 DO j=jstrm1,jendp1
1287 DO i=istru,iend
1288 grad(i,j)=ubar(i,j-1,krhs)-2.0_r8*ubar(i,j,krhs)+ &
1289 & ubar(i,j+1,krhs)
1290 END DO
1291 END DO
1293 IF (
domain(ng)%Southern_Edge(tile))
THEN
1294 DO i=istru,iend
1295 grad(i,jstr-1)=grad(i,jstr)
1296 END DO
1297 END IF
1298 END IF
1300 IF (
domain(ng)%Northern_Edge(tile))
THEN
1301 DO i=istru,iend
1302 grad(i,jend+1)=grad(i,jend)
1303 END DO
1304 END IF
1305 END IF
1306 DO j=jstr,jend+1
1307 DO i=istru-1,iend
1308 dgrad(i,j)=dvom(i-1,j)-2.0_r8*dvom(i,j)+dvom(i+1,j)
1309 END DO
1310 END DO
1311
1312 cff=1.0_r8/6.0_r8
1313 DO j=jstr,jend+1
1314 DO i=istru,iend
1315 ufe(i,j)=0.25_r8*(ubar(i,j ,krhs)+ &
1316 & ubar(i,j-1,krhs)- &
1317 & cff*(grad(i,j)+grad(i,j-1)))* &
1318 & (dvom(i,j)+dvom(i-1,j)- &
1319 & cff*(dgrad(i,j)+dgrad(i-1,j)))
1320 END DO
1321 END DO
1322
1323 DO j=jstrv,jend
1324 DO i=istrm1,iendp1
1325 grad(i,j)=vbar(i-1,j,krhs)-2.0_r8*vbar(i,j,krhs)+ &
1326 & vbar(i+1,j,krhs)
1327 END DO
1328 END DO
1330 IF (
domain(ng)%Western_Edge(tile))
THEN
1331 DO j=jstrv,jend
1332 grad(istr-1,j)=grad(istr,j)
1333 END DO
1334 END IF
1335 END IF
1337 IF (
domain(ng)%Eastern_Edge(tile))
THEN
1338 DO j=jstrv,jend
1339 grad(iend+1,j)=grad(iend,j)
1340 END DO
1341 END IF
1342 END IF
1343 DO j=jstrv-1,jend
1344 DO i=istr,iend+1
1345 dgrad(i,j)=duon(i,j-1)-2.0_r8*duon(i,j)+duon(i,j+1)
1346 END DO
1347 END DO
1348
1349 cff=1.0_r8/6.0_r8
1350 DO j=jstrv,jend
1351 DO i=istr,iend+1
1352 vfx(i,j)=0.25_r8*(vbar(i ,j,krhs)+ &
1353 & vbar(i-1,j,krhs)- &
1354 & cff*(grad(i,j)+grad(i-1,j)))* &
1355 & (duon(i,j)+duon(i,j-1)- &
1356 & cff*(dgrad(i,j)+dgrad(i,j-1)))
1357 END DO
1358 END DO
1359
1360 DO j=jstrvm1,jendp1
1361 DO i=istr,iend
1362 grad(i,j)=vbar(i,j-1,krhs)-2.0_r8*vbar(i,j,krhs)+ &
1363 & vbar(i,j+1,krhs)
1364 dgrad(i,j)=dvom(i,j-1)-2.0_r8*dvom(i,j)+dvom(i,j+1)
1365 END DO
1366 END DO
1368 IF (
domain(ng)%Southern_Edge(tile))
THEN
1369 DO i=istr,iend
1370 grad(i,jstr)=grad(i,jstr+1)
1371 dgrad(i,jstr)=dgrad(i,jstr+1)
1372 END DO
1373 END IF
1374 END IF
1376 IF (
domain(ng)%Northern_Edge(tile))
THEN
1377 DO i=istr,iend
1378 grad(i,jend+1)=grad(i,jend)
1379 dgrad(i,jend+1)=dgrad(i,jend)
1380 END DO
1381 END IF
1382 END IF
1383
1384 cff=1.0_r8/6.0_r8
1385 DO j=jstrv-1,jend
1386 DO i=istr,iend
1387 vfe(i,j)=0.25_r8*(vbar(i,j ,krhs)+ &
1388 & vbar(i,j+1,krhs)- &
1389 & cff*(grad(i,j)+grad(i,j+1)))* &
1390 & (dvom(i,j)+dvom(i,j+1)- &
1391 & cff*(dgrad(i,j)+dgrad(i,j+1)))
1392 END DO
1393 END DO
1394# endif
1395
1396
1397
1398 DO j=jstr,jend
1399 DO i=istru,iend
1400 cff1=ufx(i,j)-ufx(i-1,j)
1401 cff2=ufe(i,j+1)-ufe(i,j)
1402 fac=cff1+cff2
1403 rhs_ubar(i,j)=rhs_ubar(i,j)-fac
1404# if defined DIAGNOSTICS_UV
1405 diau2rhs(i,j,
m2xadv)=-cff1
1406 diau2rhs(i,j,
m2yadv)=-cff2
1407 diau2rhs(i,j,
m2hadv)=-fac
1408# endif
1409 END DO
1410 END DO
1411 DO j=jstrv,jend
1412 DO i=istr,iend
1413 cff1=vfx(i+1,j)-vfx(i,j)
1414 cff2=vfe(i,j)-vfe(i,j-1)
1415 fac=cff1+cff2
1416 rhs_vbar(i,j)=rhs_vbar(i,j)-fac
1417# if defined DIAGNOSTICS_UV
1418 diav2rhs(i,j,
m2xadv)=-cff1
1419 diav2rhs(i,j,
m2yadv)=-cff2
1420 diav2rhs(i,j,
m2hadv)=-fac
1421# endif
1422 END DO
1423 END DO
1424#endif
1425
1426#ifdef UV_COR
1427
1428
1429
1430
1431
1432 DO j=jstrv-1,jend
1433 DO i=istru-1,iend
1434 cff=0.5_r8*drhs(i,j)*fomn(i,j)
1435 ufx(i,j)=cff*(vbar(i,j ,krhs)+ &
1436 & vbar(i,j+1,krhs))
1437 vfe(i,j)=cff*(ubar(i ,j,krhs)+ &
1438 & ubar(i+1,j,krhs))
1439 END DO
1440 END DO
1441 DO j=jstr,jend
1442 DO i=istru,iend
1443 fac1=0.5_r8*(ufx(i,j)+ufx(i-1,j))
1444 rhs_ubar(i,j)=rhs_ubar(i,j)+fac1
1445# if defined DIAGNOSTICS_UV
1446 diau2rhs(i,j,
m2fcor)=fac1
1447# endif
1448 END DO
1449 END DO
1450 DO j=jstrv,jend
1451 DO i=istr,iend
1452 fac1=0.5_r8*(vfe(i,j)+vfe(i,j-1))
1453 rhs_vbar(i,j)=rhs_vbar(i,j)-fac1
1454# if defined DIAGNOSTICS_UV
1455 diav2rhs(i,j,
m2fcor)=-fac1
1456# endif
1457 END DO
1458 END DO
1459
1460# ifdef WEC
1461 DO j=jstrv-1,jend
1462 DO i=istru-1,iend
1463 cff=0.5_r8*drhs(i,j)*fomn(i,j)
1464 ufx(i,j)=cff*(vbar_stokes(i,j )+ &
1465 & vbar_stokes(i,j+1))
1466 vfe(i,j)=cff*(ubar_stokes(i ,j)+ &
1467 & ubar_stokes(i+1,j))
1468 END DO
1469 END DO
1470 DO j=jstr,jend
1471 DO i=istru,iend
1472 fac1=0.5_r8*(ufx(i,j)+ufx(i-1,j))
1473 rhs_ubar(i,j)=rhs_ubar(i,j)+fac1
1474# if defined DIAGNOSTICS_UV
1475 diau2rhs(i,j,m2fsco)=fac1
1476# endif
1477 END DO
1478 END DO
1479 DO j=jstrv,jend
1480 DO i=istr,iend
1481 fac1=0.5_r8*(vfe(i,j)+vfe(i,j-1))
1482 rhs_vbar(i,j)=rhs_vbar(i,j)-fac1
1483# if defined DIAGNOSTICS_UV
1484 diav2rhs(i,j,m2fsco)=-fac1
1485# endif
1486 END DO
1487 END DO
1488# endif
1489#endif
1490
1491#if defined CURVGRID && defined UV_ADV
1492
1493
1494
1495
1496
1497 DO j=jstrv-1,jend
1498 DO i=istru-1,iend
1499 cff1=0.5_r8*(vbar(i,j ,krhs)+ &
1500# ifdef WEC
1501 & vbar_stokes(i,j )+ &
1502 & vbar_stokes(i,j+1)+ &
1503# endif
1504 & vbar(i,j+1,krhs))
1505 cff2=0.5_r8*(ubar(i ,j,krhs)+ &
1506# ifdef WEC
1507 & ubar_stokes(i ,j)+ &
1508 & ubar_stokes(i+1,j)+ &
1509# endif
1510 & ubar(i+1,j,krhs))
1511 cff3=cff1*dndx(i,j)
1512 cff4=cff2*dmde(i,j)
1513# ifdef WEC_VF
1514 cff5=0.5_r8*(vbar_stokes(i,j )+ &
1515 & vbar_stokes(i,j+1))
1516 cff6=0.5_r8*(ubar_stokes(i ,j)+ &
1517 & ubar_stokes(i+1,j))
1518 cff7=cff5*dndx(i,j)
1519 cff8=cff6*dmde(i,j)
1520# endif
1521 cff=drhs(i,j)*(cff3-cff4)
1522 ufx(i,j)=cff*cff1
1523 vfe(i,j)=cff*cff2
1524# ifdef WEC_VF
1525 ufx(i,j)=ufx(i,j)-(cff5*drhs(i,j)*(cff7-cff8))
1526 vfe(i,j)=vfe(i,j)-(cff6*drhs(i,j)*(cff7-cff8))
1527# endif
1528# if defined DIAGNOSTICS_UV
1529 cff=drhs(i,j)*cff4
1530 uwrk(i,j)=-cff*cff1
1531 vwrk(i,j)=-cff*cff2
1532# ifdef WEC_VF
1533 uwrk(i,j)=uwrk(i,j)+drhs(i,j)*cff5*cff8
1534 vwrk(i,j)=vwrk(i,j)-drhs(i,j)*cff6*cff8
1535# endif
1536# endif
1537 END DO
1538 END DO
1539 DO j=jstr,jend
1540 DO i=istru,iend
1541 fac1=0.5_r8*(ufx(i,j)+ufx(i-1,j))
1542 rhs_ubar(i,j)=rhs_ubar(i,j)+fac1
1543# if defined DIAGNOSTICS_UV
1544 fac2=0.5_r8*(uwrk(i,j)+uwrk(i-1,j))
1548# endif
1549 END DO
1550 END DO
1551 DO j=jstrv,jend
1552 DO i=istr,iend
1553 fac1=0.5_r8*(vfe(i,j)+vfe(i,j-1))
1554 rhs_vbar(i,j)=rhs_vbar(i,j)-fac1
1555# if defined DIAGNOSTICS_UV
1556 fac2=0.5_r8*(vwrk(i,j)+vwrk(i,j-1))
1560# endif
1561 END DO
1562 END DO
1563#endif
1564#if defined UV_VIS2 || defined UV_VIS4
1565
1566
1567
1568
1569
1570# ifdef UV_VIS4
1571 DO j=jstrm1,jendp2
1572 DO i=istrm1,iendp2
1573# else
1574 DO j=jstr,jend+1
1575 DO i=istr,iend+1
1576# endif
1577 drhs_p(i,j)=0.25_r8*(drhs(i,j )+drhs(i-1,j )+ &
1578 & drhs(i,j-1)+drhs(i-1,j-1))
1579 END DO
1580 END DO
1581#endif
1582#ifdef UV_VIS2
1583
1584
1585
1586
1587
1588
1589
1590
1591 DO j=jstrv-1,jend
1592 DO i=istru-1,iend
1593 cff=visc2_r(i,j)*drhs(i,j)*0.5_r8* &
1594 & (pmon_r(i,j)* &
1595 & ((pn(i ,j)+pn(i+1,j))*ubar(i+1,j,krhs)- &
1596 & (pn(i-1,j)+pn(i ,j))*ubar(i ,j,krhs))- &
1597 & pnom_r(i,j)* &
1598 & ((pm(i,j )+pm(i,j+1))*vbar(i,j+1,krhs)- &
1599 & (pm(i,j-1)+pm(i,j ))*vbar(i,j ,krhs)))
1600 ufx(i,j)=on_r(i,j)*on_r(i,j)*cff
1601 vfe(i,j)=om_r(i,j)*om_r(i,j)*cff
1602 END DO
1603 END DO
1604 DO j=jstr,jend+1
1605 DO i=istr,iend+1
1606 cff=visc2_p(i,j)*drhs_p(i,j)*0.5_r8* &
1607 & (pmon_p(i,j)* &
1608 & ((pn(i ,j-1)+pn(i ,j))*vbar(i ,j,krhs)- &
1609 & (pn(i-1,j-1)+pn(i-1,j))*vbar(i-1,j,krhs))+ &
1610 & pnom_p(i,j)* &
1611 & ((pm(i-1,j )+pm(i,j ))*ubar(i,j ,krhs)- &
1612 & (pm(i-1,j-1)+pm(i,j-1))*ubar(i,j-1,krhs)))
1613# ifdef MASKING
1614 cff=cff*pmask(i,j)
1615# endif
1616# ifdef WET_DRY
1617 cff=cff*pmask_wet(i,j)
1618# endif
1619 ufe(i,j)=om_p(i,j)*om_p(i,j)*cff
1620 vfx(i,j)=on_p(i,j)*on_p(i,j)*cff
1621 END DO
1622 END DO
1623
1624
1625
1626 DO j=jstr,jend
1627 DO i=istru,iend
1628 cff1=0.5_r8*(pn(i-1,j)+pn(i,j))*(ufx(i,j )-ufx(i-1,j))
1629 cff2=0.5_r8*(pm(i-1,j)+pm(i,j))*(ufe(i,j+1)-ufe(i ,j))
1630 fac=cff1+cff2
1631 rhs_ubar(i,j)=rhs_ubar(i,j)+fac
1632# if defined DIAGNOSTICS_UV
1634 diau2rhs(i,j,
m2xvis)=cff1
1635 diau2rhs(i,j,
m2yvis)=cff2
1636# endif
1637 END DO
1638 END DO
1639 DO j=jstrv,jend
1640 DO i=istr,iend
1641 cff1=0.5_r8*(pn(i,j-1)+pn(i,j))*(vfx(i+1,j)-vfx(i,j ))
1642 cff2=0.5_r8*(pm(i,j-1)+pm(i,j))*(vfe(i ,j)-vfe(i,j-1))
1643 fac=cff1-cff2
1644 rhs_vbar(i,j)=rhs_vbar(i,j)+fac
1645# if defined DIAGNOSTICS_UV
1647 diav2rhs(i,j,
m2xvis)= cff1
1648 diav2rhs(i,j,
m2yvis)=-cff2
1649# endif
1650 END DO
1651 END DO
1652#endif
1653#ifdef UV_VIS4
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667 DO j=jstrvm2,jendp1
1668 DO i=istrum2,iendp1
1669 cff=visc4_r(i,j)*0.5_r8* &
1670 & (pmon_r(i,j)* &
1671 & ((pn(i ,j)+pn(i+1,j))*ubar(i+1,j,krhs)- &
1672 & (pn(i-1,j)+pn(i ,j))*ubar(i ,j,krhs))- &
1673 & pnom_r(i,j)* &
1674 & ((pm(i,j )+pm(i,j+1))*vbar(i,j+1,krhs)- &
1675 & (pm(i,j-1)+pm(i,j ))*vbar(i,j ,krhs)))
1676 ufx(i,j)=on_r(i,j)*on_r(i,j)*cff
1677 vfe(i,j)=om_r(i,j)*om_r(i,j)*cff
1678 END DO
1679 END DO
1680 DO j=jstrm1,jendp2
1681 DO i=istrm1,iendp2
1682 cff=visc4_p(i,j)*0.5_r8* &
1683 & (pmon_p(i,j)* &
1684 & ((pn(i ,j-1)+pn(i ,j))*vbar(i ,j,krhs)- &
1685 & (pn(i-1,j-1)+pn(i-1,j))*vbar(i-1,j,krhs))+ &
1686 & pnom_p(i,j)* &
1687 & ((pm(i-1,j )+pm(i,j ))*ubar(i,j ,krhs)- &
1688 & (pm(i-1,j-1)+pm(i,j-1))*ubar(i,j-1,krhs)))
1689# ifdef MASKING
1690 cff=cff*pmask(i,j)
1691# endif
1692# ifdef WET_DRY
1693 cff=cff*pmask_wet(i,j)
1694# endif
1695 ufe(i,j)=om_p(i,j)*om_p(i,j)*cff
1696 vfx(i,j)=on_p(i,j)*on_p(i,j)*cff
1697 END DO
1698 END DO
1699
1700
1701
1702 DO j=jstrm1,jendp1
1703 DO i=istrum1,iendp1
1704 lapu(i,j)=0.125_r8* &
1705 & (pm(i-1,j)+pm(i,j))*(pn(i-1,j)+pn(i,j))* &
1706 & ((pn(i-1,j)+pn(i,j))* &
1707 & (ufx(i,j )-ufx(i-1,j))+ &
1708 & (pm(i-1,j)+pm(i,j))* &
1709 & (ufe(i,j+1)-ufe(i ,j)))
1710 END DO
1711 END DO
1712 DO j=jstrvm1,jendp1
1713 DO i=istrm1,iendp1
1714 lapv(i,j)=0.125_r8* &
1715 & (pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))* &
1716 & ((pn(i,j-1)+pn(i,j))* &
1717 & (vfx(i+1,j)-vfx(i,j ))- &
1718 & (pm(i,j-1)+pm(i,j))* &
1719 & (vfe(i ,j)-vfe(i,j-1)))
1720 END DO
1721 END DO
1722
1723
1724
1725
1726
1728 IF (
domain(ng)%Western_Edge(tile))
THEN
1730 DO j=jstrm1,jendp1
1731 lapu(istru-1,j)=0.0_r8
1732 END DO
1733 ELSE
1734 DO j=jstrm1,jendp1
1735 lapu(istru-1,j)=lapu(istru,j)
1736 END DO
1737 END IF
1739 DO j=jstrvm1,jendp1
1740 lapv(istr-1,j)=
gamma2(ng)*lapv(istr,j)
1741 END DO
1742 ELSE
1743 DO j=jstrvm1,jendp1
1744 lapv(istr-1,j)=0.0_r8
1745 END DO
1746 END IF
1747 END IF
1748 END IF
1749
1751 IF (
domain(ng)%Eastern_Edge(tile))
THEN
1753 DO j=jstrm1,jendp1
1754 lapu(iend+1,j)=0.0_r8
1755 END DO
1756 ELSE
1757 DO j=jstrm1,jendp1
1758 lapu(iend+1,j)=lapu(iend,j)
1759 END DO
1760 END IF
1762 DO j=jstrvm1,jendp1
1763 lapv(iend+1,j)=
gamma2(ng)*lapv(iend,j)
1764 END DO
1765 ELSE
1766 DO j=jstrvm1,jendp1
1767 lapv(iend+1,j)=0.0_r8
1768 END DO
1769 END IF
1770 END IF
1771 END IF
1772
1774 IF (
domain(ng)%Southern_Edge(tile))
THEN
1776 DO i=istrum1,iendp1
1777 lapu(i,jstr-1)=
gamma2(ng)*lapu(i,jstr)
1778 END DO
1779 ELSE
1780 DO i=istrum1,iendp1
1781 lapu(i,jstr-1)=0.0_r8
1782 END DO
1783 END IF
1785 DO i=istrm1,iendp1
1786 lapv(i,jstrv-1)=0.0_r8
1787 END DO
1788 ELSE
1789 DO i=istrm1,iendp1
1790 lapv(i,jstrv-1)=lapv(i,jstrv)
1791 END DO
1792 END IF
1793 END IF
1794 END IF
1795
1797 IF (
domain(ng)%Northern_Edge(tile))
THEN
1799 DO i=istrum1,iendp1
1800 lapu(i,jend+1)=
gamma2(ng)*lapu(i,jend)
1801 END DO
1802 ELSE
1803 DO i=istrum1,iendp1
1804 lapu(i,jend+1)=0.0_r8
1805 END DO
1806 END IF
1808 DO i=istrm1,iendp1
1809 lapv(i,jend+1)=0.0_r8
1810 END DO
1811 ELSE
1812 DO i=istrm1,iendp1
1813 lapv(i,jend+1)=lapv(i,jend)
1814 END DO
1815 END IF
1816 END IF
1817 END IF
1818
1821 IF (
domain(ng)%SouthWest_Corner(tile))
THEN
1822 lapu(istr ,jstr-1)=0.5_r8*(lapu(istr+1,jstr-1)+ &
1823 & lapu(istr ,jstr ))
1824 lapv(istr-1,jstr )=0.5_r8*(lapv(istr-1,jstr+1)+ &
1825 & lapv(istr ,jstr ))
1826 END IF
1827 END IF
1828
1831 IF (
domain(ng)%SouthEast_Corner(tile))
THEN
1832 lapu(iend+1,jstr-1)=0.5_r8*(lapu(iend ,jstr-1)+ &
1833 & lapu(iend+1,jstr ))
1834 lapv(iend+1,jstr )=0.5_r8*(lapv(iend ,jstr )+ &
1835 & lapv(iend+1,jstr+1))
1836 END IF
1837 END IF
1838
1841 IF (
domain(ng)%NorthWest_Corner(tile))
THEN
1842 lapu(istr ,jend+1)=0.5_r8*(lapu(istr+1,jend+1)+ &
1843 & lapu(istr ,jend ))
1844 lapv(istr-1,jend+1)=0.5_r8*(lapv(istr ,jend+1)+ &
1845 & lapv(istr-1,jend ))
1846 END IF
1847 END IF
1848
1851 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
1852 lapu(iend+1,jend+1)=0.5_r8*(lapu(iend ,jend+1)+ &
1853 & lapu(iend+1,jend ))
1854 lapv(iend+1,jend+1)=0.5_r8*(lapv(iend ,jend+1)+ &
1855 & lapv(iend+1,jend ))
1856 END IF
1857 END IF
1858
1859
1860
1861
1862 DO j=jstrv-1,jend
1863 DO i=istru-1,iend
1864 cff=visc4_r(i,j)*drhs(i,j)*0.5_r8* &
1865 & (pmon_r(i,j)* &
1866 & ((pn(i ,j)+pn(i+1,j))*lapu(i+1,j)- &
1867 & (pn(i-1,j)+pn(i ,j))*lapu(i ,j))- &
1868 & pnom_r(i,j)* &
1869 & ((pm(i,j )+pm(i,j+1))*lapv(i,j+1)- &
1870 & (pm(i,j-1)+pm(i,j ))*lapv(i,j )))
1871 ufx(i,j)=on_r(i,j)*on_r(i,j)*cff
1872 vfe(i,j)=om_r(i,j)*om_r(i,j)*cff
1873 END DO
1874 END DO
1875 DO j=jstr,jend+1
1876 DO i=istr,iend+1
1877 cff=visc4_p(i,j)*drhs_p(i,j)*0.5_r8* &
1878 & (pmon_p(i,j)* &
1879 & ((pn(i ,j-1)+pn(i ,j))*lapv(i ,j)- &
1880 & (pn(i-1,j-1)+pn(i-1,j))*lapv(i-1,j))+ &
1881 & pnom_p(i,j)* &
1882 & ((pm(i-1,j )+pm(i,j ))*lapu(i,j )- &
1883 & (pm(i-1,j-1)+pm(i,j-1))*lapu(i,j-1)))
1884# ifdef MASKING
1885 cff=cff*pmask(i,j)
1886# endif
1887# ifdef WET_DRY
1888 cff=cff*pmask_wet(i,j)
1889# endif
1890 ufe(i,j)=om_p(i,j)*om_p(i,j)*cff
1891 vfx(i,j)=on_p(i,j)*on_p(i,j)*cff
1892 END DO
1893 END DO
1894
1895
1896
1897 DO j=jstr,jend
1898 DO i=istru,iend
1899 cff1=0.5_r8*(pn(i-1,j)+pn(i,j))*(ufx(i,j )-ufx(i-1,j))
1900 cff2=0.5_r8*(pm(i-1,j)+pm(i,j))*(ufe(i,j+1)-ufe(i ,j))
1901 fac=cff1+cff2
1902 rhs_ubar(i,j)=rhs_ubar(i,j)-fac
1903# if defined DIAGNOSTICS_UV
1904 diau2rhs(i,j,
m2hvis)=-fac
1905 diau2rhs(i,j,
m2xvis)=-cff1
1906 diau2rhs(i,j,
m2yvis)=-cff2
1907# endif
1908 END DO
1909 END DO
1910 DO j=jstrv,jend
1911 DO i=istr,iend
1912 cff1=0.5_r8*(pn(i,j-1)+pn(i,j))*(vfx(i+1,j)-vfx(i,j ))
1913 cff2=0.5_r8*(pm(i,j-1)+pm(i,j))*(vfe(i ,j)-vfe(i,j-1))
1914 fac=cff1-cff2
1915 rhs_vbar(i,j)=rhs_vbar(i,j)-fac
1916# if defined DIAGNOSTICS_UV
1917 diav2rhs(i,j,
m2hvis)=-fac
1918 diav2rhs(i,j,
m2xvis)=-cff1
1919 diav2rhs(i,j,
m2yvis)= cff2
1920# endif
1921 END DO
1922 END DO
1923#endif
1924#if defined WEC_VF && defined SOLVE3D
1925
1926
1927
1928
1929
1930 DO j=jstr,jend
1931 DO i=istru,iend
1932 cff1=rubrk2d(i,j)
1933 rhs_ubar(i,j)=rhs_ubar(i,j)+cff1
1934# ifdef DIAGNOSTICS_UV
1935 diau2rhs(i,j,m2wbrk)=cff1
1936# endif
1937# ifdef WEC_ROLLER
1938 cff1=rurol2d(i,j)
1939 rhs_ubar(i,j)=rhs_ubar(i,j)+cff1
1940# ifdef DIAGNOSTICS_UV
1941 diau2rhs(i,j,m2wrol)=cff1
1942# endif
1943# else
1944# ifdef DIAGNOSTICS_UV
1945 diau2rhs(i,j,m2wrol)=0.0_r8
1946# endif
1947# endif
1948# ifdef BOTTOM_STREAMING
1949 cff1=rubst2d(i,j)
1950 rhs_ubar(i,j)=rhs_ubar(i,j)+cff1
1951# ifdef DIAGNOSTICS_UV
1952 diau2rhs(i,j,m2bstm)=cff1
1953# endif
1954# endif
1955# ifdef SURFACE_STREAMING
1956 cff1=russt2d(i,j)
1957 rhs_ubar(i,j)=rhs_ubar(i,j)+cff1
1958# ifdef DIAGNOSTICS_UV
1959 diau2rhs(i,j,m2sstm)=cff1
1960# endif
1961# endif
1962 END DO
1963 IF (j.ge.jstrv) THEN
1964 DO i=istr,iend
1965 cff1=rvbrk2d(i,j)
1966 rhs_vbar(i,j)=rhs_vbar(i,j)+cff1
1967# ifdef DIAGNOSTICS_UV
1968 diav2rhs(i,j,m2wbrk)=cff1
1969# endif
1970# ifdef WEC_ROLLER
1971 cff1=rvrol2d(i,j)
1972 rhs_vbar(i,j)=rhs_vbar(i,j)+cff1
1973# ifdef DIAGNOSTICS_UV
1974 diav2rhs(i,j,m2wrol)=cff1
1975# endif
1976# else
1977# ifdef DIAGNOSTICS_UV
1978 diav2rhs(i,j,m2wrol)=0.0_r8
1979# endif
1980# endif
1981# ifdef BOTTOM_STREAMING
1982 cff1=rvbst2d(i,j)
1983 rhs_vbar(i,j)=rhs_vbar(i,j)+cff1
1984# ifdef DIAGNOSTICS_UV
1985 diav2rhs(i,j,m2bstm)=cff1
1986# endif
1987# endif
1988# ifdef SURFACE_STREAMING
1989 cff1=rvsst2d(i,j)
1990 rhs_vbar(i,j)=rhs_vbar(i,j)+cff1
1991# ifdef DIAGNOSTICS_UV
1992 diav2rhs(i,j,m2sstm)=cff1
1993# endif
1994# endif
1995 END DO
1996 END IF
1997 END DO
1998# ifdef UV_ADV
1999# ifdef DIAGNOSTICS_UV
2000
2001
2002
2003
2004
2005
2006 DO j=jstr,jend
2007 DO i=istru,iend
2008 cff=0.5_r8*(drhs(i-1,j)+drhs(i,j))
2009 dvsom(i,j)=0.25_r8*cff*om_u(i,j)* &
2010 & (vbar_stokes(i ,j )+ &
2011 & vbar_stokes(i ,j+1)+ &
2012 & vbar_stokes(i-1,j )+ &
2013 & vbar_stokes(i-1,j+1))
2014 END DO
2015 END DO
2016 DO j=jstr,jend+1
2017 DO i=istru,iend
2018 ufx(i,j)=0.5_r8*(ubar(i ,j-1,krhs)+ &
2019 ubar(i ,j ,krhs))
2020 END DO
2021 END DO
2022 DO j=jstr,jend
2023 DO i=istru,iend
2024 cff1=ufx(i,j+1)-ufx(i,j)
2025 cff=cff1*dvsom(i,j)
2028 diau2rhs(i,j,m2hjvf)=-cff
2029 END DO
2030 END DO
2031 DO j=jstrv,jend
2032 DO i=istr,iend
2033 cff=0.5_r8*(drhs(i,j)+drhs(i,j-1))
2034 duson(i,j)=cff*0.25_r8*on_v(i,j)* &
2035 & (ubar_stokes(i ,j )+ &
2036 & ubar_stokes(i+1,j )+ &
2037 & ubar_stokes(i ,j-1)+ &
2038 & ubar_stokes(i+1,j-1))
2039 END DO
2040 END DO
2041 DO j=jstrv,jend
2042 DO i=istr,iend+1
2043 vfe(i,j)=0.5_r8*(vbar(i-1,j ,krhs)+ &
2044 & vbar(i ,j ,krhs))
2045 END DO
2046 END DO
2047 DO i=istr,iend
2048 DO j=jstrv,jend
2049 cff2=vfe(i+1,j)-vfe(i,j)
2050 cff=cff2*duson(i,j)
2053 diav2rhs(i,j,m2hjvf)=-cff
2054 END DO
2055 END DO
2056# endif
2057
2058
2059
2060
2061
2062
2063
2064 DO j=jstr,jend
2065 DO i=istru,iend
2066 cff=0.5_r8*(drhs(i-1,j)+drhs(i,j))
2067 duson(i,j)=cff*on_u(i,j)*ubar_stokes(i,j)
2068 dvson(i,j)=0.25_r8*cff*on_u(i,j)* &
2069 & (vbar_stokes(i ,j )+ &
2070 & vbar_stokes(i ,j+1)+ &
2071 & vbar_stokes(i-1,j )+ &
2072 & vbar_stokes(i-1,j+1))
2073 END DO
2074 DO i=istru-1,iend
2075 ufx(i,j)=0.5_r8*(ubar(i ,j ,krhs)+ &
2076 ubar(i+1,j ,krhs))
2077 vfx(i,j)=0.5_r8*(vbar(i ,j ,krhs)+ &
2078 & vbar(i ,j+1,krhs))
2079 END DO
2080 END DO
2081 DO j=jstrv,jend
2082 DO i=istr,iend
2083 cff=0.5_r8*(drhs(i,j)+drhs(i,j-1))
2084 dusom(i,j)=cff*0.25_r8*om_v(i,j)* &
2085 & (ubar_stokes(i ,j )+ &
2086 & ubar_stokes(i+1,j )+ &
2087 & ubar_stokes(i ,j-1)+ &
2088 & ubar_stokes(i+1,j-1))
2089 dvsom(i,j)=cff*om_v(i,j)*vbar_stokes(i,j)
2090 END DO
2091 END DO
2092 DO j=jstrv-1,jend
2093 DO i=istr,iend
2094 cff=0.5_r8*(drhs(i,j)+drhs(i,j-1))
2095 ufe(i,j)=0.5_r8*(ubar(i+1,j ,krhs)+ &
2096 & ubar(i ,j ,krhs))
2097 vfe(i,j)=0.5_r8*(vbar(i ,j ,krhs)+ &
2098 & vbar(i ,j+1,krhs))
2099 END DO
2100 END DO
2101 DO j=jstr,jend
2102 DO i=istru,iend
2103 cff1=ufx(i,j)-ufx(i-1,j)
2104 cff2=vfx(i,j)-vfx(i-1,j)
2105 cff3=duson(i,j)*cff1
2106 cff4=dvson(i,j)*cff2
2107 rhs_ubar(i,j)=rhs_ubar(i,j)+cff3+cff4
2108
2109# ifdef DIAGNOSTICS_UV
2112 diau2rhs(i,j,m2hjvf)=diau2rhs(i,j,m2hjvf)+cff4
2113# endif
2114 END DO
2115 END DO
2116 DO i=istr,iend
2117 DO j=jstrv,jend
2118 cff1=ufe(i,j)-ufe(i,j-1)
2119 cff2=vfe(i,j)-vfe(i,j-1)
2120 cff3=dusom(i,j)*cff1
2121 cff4=dvsom(i,j)*cff2
2122 rhs_vbar(i,j)=rhs_vbar(i,j)+cff3+cff4
2123
2124# ifdef DIAGNOSTICS_UV
2127 diav2rhs(i,j,m2hjvf)=diav2rhs(i,j,m2hjvf)+cff3
2128# endif
2129 END DO
2130 END DO
2131# endif
2132#endif
2133#ifndef SOLVE3D
2134
2135
2136
2137
2138
2139 DO j=jstr,jend
2140 DO i=istru,iend
2141 fac=bustr(i,j)*om_u(i,j)*on_u(i,j)
2142 rhs_ubar(i,j)=rhs_ubar(i,j)-fac
2143# ifdef DIAGNOSTICS_UV
2144 diau2rhs(i,j,
m2bstr)=-fac
2145# endif
2146 END DO
2147 END DO
2148 DO j=jstrv,jend
2149 DO i=istr,iend
2150 fac=bvstr(i,j)*om_v(i,j)*on_v(i,j)
2151 rhs_vbar(i,j)=rhs_vbar(i,j)-fac
2152# ifdef DIAGNOSTICS_UV
2153 diav2rhs(i,j,
m2bstr)=-fac
2154# endif
2155 END DO
2156 END DO
2157#else
2158# ifdef DIAGNOSTICS_UV
2159
2160
2161
2162 DO j=jstr,jend
2163 DO i=istru,iend
2164 diau2rhs(i,j,
m2bstr)=0.0_r8
2165 END DO
2166 END DO
2167 DO j=jstrv,jend
2168 DO i=istr,iend
2169 diav2rhs(i,j,
m2bstr)=0.0_r8
2170 END DO
2171 END DO
2172# endif
2173#endif
2174
2175
2176
2177
2178
2180 DO j=jstr,jend
2181 DO i=istru,iend
2182 cff=0.25_r8*(
clima(ng)%M2nudgcof(i-1,j)+ &
2183 &
clima(ng)%M2nudgcof(i ,j))* &
2184 & om_u(i,j)*on_u(i,j)
2185 rhs_ubar(i,j)=rhs_ubar(i,j)+ &
2186 & cff*(drhs(i-1,j)+drhs(i,j))* &
2187 & (
clima(ng)%ubarclm(i,j)- &
2188 & ubar(i,j,krhs))
2189 END DO
2190 END DO
2191 DO j=jstrv,jend
2192 DO i=istr,iend
2193 cff=0.25_r8*(
clima(ng)%M2nudgcof(i,j-1)+ &
2194 &
clima(ng)%M2nudgcof(i,j ))* &
2195 & om_v(i,j)*on_v(i,j)
2196 rhs_vbar(i,j)=rhs_vbar(i,j)+ &
2197 & cff*(drhs(i,j-1)+drhs(i,j))* &
2198 & (
clima(ng)%vbarclm(i,j)- &
2199 & vbar(i,j,krhs))
2200 END DO
2201 END DO
2202 END IF
2203
2204#ifdef SOLVE3D
2205# ifdef WET_DRY
2206 DO j=jstr,jend
2207 DO i=istru,iend
2208 cff5=abs(abs(umask_wet(i,j))-1.0_r8)
2209 cff6=0.5_r8+dsign(0.5_r8,rhs_ubar(i,j))*umask_wet(i,j)
2210 cff7=0.5_r8*umask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
2211 rhs_ubar(i,j)=rhs_ubar(i,j)*cff7
2212 END DO
2213 END DO
2214 DO j=jstrv,jend
2215 DO i=istr,iend
2216 cff5=abs(abs(vmask_wet(i,j))-1.0_r8)
2217 cff6=0.5_r8+dsign(0.5_r8,rhs_vbar(i,j))*vmask_wet(i,j)
2218 cff7=0.5_r8*vmask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
2219 rhs_vbar(i,j)=rhs_vbar(i,j)*cff7
2220 END DO
2221 END DO
2222# endif
2223
2224
2225
2226
2227
2228
2229
2230
2231
2232
2233
2234
2235
2236
2237
2238
2239
2240
2243 DO j=jstr,jend
2244 DO i=istru,iend
2245 rufrc(i,j)=rufrc(i,j)-rhs_ubar(i,j)
2246 rhs_ubar(i,j)=rhs_ubar(i,j)+rufrc(i,j)
2247 ru(i,j,0,nstp)=rufrc(i,j)
2248# ifdef DIAGNOSTICS_UV
2250 diarufrc(i,j,3,idiag)=diarufrc(i,j,3,idiag)- &
2251 & diau2rhs(i,j,idiag)
2252 diau2rhs(i,j,idiag)=diau2rhs(i,j,idiag)+ &
2253 & diarufrc(i,j,3,idiag)
2254 diarufrc(i,j,nstp,idiag)=diarufrc(i,j,3,idiag)
2255 END DO
2260# ifdef WEC_VF
2261 diau2rhs(i,j,m2zeta)=diau2rhs(i,j,m2zeta)+ &
2263# endif
2264# endif
2265 END DO
2266 END DO
2267 DO j=jstrv,jend
2268 DO i=istr,iend
2269 rvfrc(i,j)=rvfrc(i,j)-rhs_vbar(i,j)
2270 rhs_vbar(i,j)=rhs_vbar(i,j)+rvfrc(i,j)
2271 rv(i,j,0,nstp)=rvfrc(i,j)
2272# ifdef DIAGNOSTICS_UV
2274 diarvfrc(i,j,3,idiag)=diarvfrc(i,j,3,idiag)- &
2275 & diav2rhs(i,j,idiag)
2276 diav2rhs(i,j,idiag)=diav2rhs(i,j,idiag)+ &
2277 & diarvfrc(i,j,3,idiag)
2278 diarvfrc(i,j,nstp,idiag)=diarvfrc(i,j,3,idiag)
2279 END DO
2284# ifdef WEC_VF
2285 diav2rhs(i,j,m2zeta)=diav2rhs(i,j,m2zeta)+ &
2287# endif
2288# endif
2289 END DO
2290 END DO
2292 DO j=jstr,jend
2293 DO i=istru,iend
2294 rufrc(i,j)=rufrc(i,j)-rhs_ubar(i,j)
2295 rhs_ubar(i,j)=rhs_ubar(i,j)+ &
2296 & 1.5_r8*rufrc(i,j)-0.5_r8*ru(i,j,0,nnew)
2297 ru(i,j,0,nstp)=rufrc(i,j)
2298# ifdef DIAGNOSTICS_UV
2300 diarufrc(i,j,3,idiag)=diarufrc(i,j,3,idiag)- &
2301 & diau2rhs(i,j,idiag)
2302 diau2rhs(i,j,idiag)=diau2rhs(i,j,idiag)+ &
2303 & 1.5_r8*diarufrc(i,j,3,idiag)- &
2304 & 0.5_r8*diarufrc(i,j,nnew,idiag)
2305 diarufrc(i,j,nstp,idiag)=diarufrc(i,j,3,idiag)
2306 END DO
2308 & 0.5_r8*diarufrc(i,j,nnew,
m2sstr)
2311 & 0.5_r8*diarufrc(i,j,nnew,
m2bstr)
2313# ifdef WEC_VF
2314 diau2rhs(i,j,m2zeta)=diau2rhs(i,j,m2zeta)+ &
2315 & 1.5_r8*diarufrc(i,j,3,
m2pgrd)- &
2316 & 0.5_r8*diarufrc(i,j,nnew,
m2pgrd)
2317# endif
2318# endif
2319 END DO
2320 END DO
2321 DO j=jstrv,jend
2322 DO i=istr,iend
2323 rvfrc(i,j)=rvfrc(i,j)-rhs_vbar(i,j)
2324 rhs_vbar(i,j)=rhs_vbar(i,j)+ &
2325 & 1.5_r8*rvfrc(i,j)-0.5_r8*rv(i,j,0,nnew)
2326 rv(i,j,0,nstp)=rvfrc(i,j)
2327# ifdef DIAGNOSTICS_UV
2329 diarvfrc(i,j,3,idiag)=diarvfrc(i,j,3,idiag)- &
2330 & diav2rhs(i,j,idiag)
2331 diav2rhs(i,j,idiag)=diav2rhs(i,j,idiag)+ &
2332 & 1.5_r8*diarvfrc(i,j,3,idiag)- &
2333 & 0.5_r8*diarvfrc(i,j,nnew,idiag)
2334 diarvfrc(i,j,nstp,idiag)=diarvfrc(i,j,3,idiag)
2335 END DO
2337 & 0.5_r8*diarvfrc(i,j,nnew,
m2sstr)
2340 & 0.5_r8*diarvfrc(i,j,nnew,
m2bstr)
2342# ifdef WEC_VF
2343 diav2rhs(i,j,m2zeta)=diav2rhs(i,j,m2zeta)+ &
2344 & 1.5_r8*diarvfrc(i,j,3,
m2pgrd)- &
2345 & 0.5_r8*diarvfrc(i,j,nnew,
m2pgrd)
2346# endif
2347# endif
2348 END DO
2349 END DO
2350 ELSE
2351 cff1=23.0_r8/12.0_r8
2352 cff2=16.0_r8/12.0_r8
2353 cff3= 5.0_r8/12.0_r8
2354 DO j=jstr,jend
2355 DO i=istru,iend
2356 rufrc(i,j)=rufrc(i,j)-rhs_ubar(i,j)
2357 rhs_ubar(i,j)=rhs_ubar(i,j)+ &
2358 & cff1*rufrc(i,j)- &
2359 & cff2*ru(i,j,0,nnew)+ &
2360 & cff3*ru(i,j,0,nstp)
2361 ru(i,j,0,nstp)=rufrc(i,j)
2362# ifdef DIAGNOSTICS_UV
2364 diarufrc(i,j,3,idiag)=diarufrc(i,j,3,idiag)- &
2365 & diau2rhs(i,j,idiag)
2366 diau2rhs(i,j,idiag)=diau2rhs(i,j,idiag)+ &
2367 & cff1*diarufrc(i,j,3,idiag)- &
2368 & cff2*diarufrc(i,j,nnew,idiag)+ &
2369 & cff3*diarufrc(i,j,nstp,idiag)
2370 diarufrc(i,j,nstp,idiag)=diarufrc(i,j,3,idiag)
2371 END DO
2373 & cff2*diarufrc(i,j,nnew,
m2sstr)+ &
2374 & cff3*diarufrc(i,j,nstp,
m2sstr)
2377 & cff2*diarufrc(i,j,nnew,
m2bstr)+ &
2378 & cff3*diarufrc(i,j,nstp,
m2bstr)
2380# ifdef WEC_VF
2381 diau2rhs(i,j,m2zeta)=diau2rhs(i,j,m2zeta)+ &
2382 & cff1*diarufrc(i,j,3,
m2pgrd)- &
2383 & cff2*diarufrc(i,j,nnew,
m2pgrd)+ &
2384 & cff3*diarufrc(i,j,nstp,
m2pgrd)
2385# endif
2386# endif
2387 END DO
2388 END DO
2389 DO j=jstrv,jend
2390 DO i=istr,iend
2391 rvfrc(i,j)=rvfrc(i,j)-rhs_vbar(i,j)
2392 rhs_vbar(i,j)=rhs_vbar(i,j)+ &
2393 & cff1*rvfrc(i,j)- &
2394 & cff2*rv(i,j,0,nnew)+ &
2395 & cff3*rv(i,j,0,nstp)
2396 rv(i,j,0,nstp)=rvfrc(i,j)
2397# ifdef DIAGNOSTICS_UV
2399 diarvfrc(i,j,3,idiag)=diarvfrc(i,j,3,idiag)- &
2400 & diav2rhs(i,j,idiag)
2401 diav2rhs(i,j,idiag)=diav2rhs(i,j,idiag)+ &
2402 & cff1*diarvfrc(i,j,3,idiag)- &
2403 & cff2*diarvfrc(i,j,nnew,idiag)+ &
2404 & cff3*diarvfrc(i,j,nstp,idiag)
2405 diarvfrc(i,j,nstp,idiag)=diarvfrc(i,j,3,idiag)
2406 END DO
2408 & cff2*diarvfrc(i,j,nnew,
m2sstr)+ &
2409 & cff3*diarvfrc(i,j,nstp,
m2sstr)
2412 & cff2*diarvfrc(i,j,nnew,
m2bstr)+ &
2413 & cff3*diarvfrc(i,j,nstp,
m2bstr)
2415# ifdef WEC_VF
2416 diav2rhs(i,j,m2zeta)=diav2rhs(i,j,m2zeta)+ &
2417 & cff1*diarvfrc(i,j,3,
m2pgrd)- &
2418 & cff2*diarvfrc(i,j,nnew,
m2pgrd)+ &
2419 & cff3*diarvfrc(i,j,nstp,
m2pgrd)
2420# endif
2421# endif
2422 END DO
2423 END DO
2424 END IF
2425 ELSE
2426 DO j=jstr,jend
2427 DO i=istru,iend
2428 rhs_ubar(i,j)=rhs_ubar(i,j)+rufrc(i,j)
2429# ifdef DIAGNOSTICS_UV
2431 diau2rhs(i,j,idiag)=diau2rhs(i,j,idiag)+ &
2432 & diarufrc(i,j,3,idiag)
2433 END DO
2436# ifdef WEC_VF
2437 diau2rhs(i,j,m2zeta)=diau2rhs(i,j,m2zeta)+ &
2439# endif
2440# endif
2441 END DO
2442 END DO
2443 DO j=jstrv,jend
2444 DO i=istr,iend
2445 rhs_vbar(i,j)=rhs_vbar(i,j)+rvfrc(i,j)
2446# ifdef DIAGNOSTICS_UV
2448 diav2rhs(i,j,idiag)=diav2rhs(i,j,idiag)+ &
2449 & diarvfrc(i,j,3,idiag)
2450 END DO
2453# ifdef WEC_VF
2454 diav2rhs(i,j,m2zeta)=diav2rhs(i,j,m2zeta)+ &
2456# endif
2457# endif
2458 END DO
2459 END DO
2460 END IF
2461#else
2462
2463
2464
2465
2466
2467 DO j=jstr,jend
2468 DO i=istru,iend
2469 fac=sustr(i,j)*om_u(i,j)*on_u(i,j)
2470 rhs_ubar(i,j)=rhs_ubar(i,j)+fac
2471# ifdef DIAGNOSTICS_UV
2473# endif
2474 END DO
2475 END DO
2476 DO j=jstrv,jend
2477 DO i=istr,iend
2478 fac=svstr(i,j)*om_v(i,j)*on_v(i,j)
2479 rhs_vbar(i,j)=rhs_vbar(i,j)+fac
2480# ifdef DIAGNOSTICS_UV
2482# endif
2483 END DO
2484 END DO
2485#endif
2486
2487
2488
2489
2490
2491
2492
2493 DO j=jstrv-1,jend
2494 DO i=istru-1,iend
2495 dstp(i,j)=zeta(i,j,kstp)+h(i,j)
2496 END DO
2497 END DO
2498
2499
2500
2501
2502
2503 IF (first_2d_step) THEN
2505#ifdef WET_DRY
2506 cff2=1.0_r8/cff1
2507#endif
2508 DO j=jstr,jend
2509 DO i=istru,iend
2510 cff=(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))
2511 fac=1.0_r8/(dnew(i,j)+dnew(i-1,j))
2512 ubar(i,j,knew)=(ubar(i,j,kstp)* &
2513 & (dstp(i,j)+dstp(i-1,j))+ &
2514 & cff*cff1*rhs_ubar(i,j))*fac
2515#ifdef MASKING
2516 ubar(i,j,knew)=ubar(i,j,knew)*umask(i,j)
2517#endif
2518#ifdef WET_DRY
2519 cff5=abs(abs(umask_wet(i,j))-1.0_r8)
2520 cff6=0.5_r8+dsign(0.5_r8,ubar(i,j,knew))*umask_wet(i,j)
2521 cff7=0.5_r8*umask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
2522 ubar(i,j,knew)=ubar(i,j,knew)*cff7
2523 rhs_ubar(i,j)=rhs_ubar(i,j)*cff7
2524# ifdef SOLVE3D
2526 rufrc(i,j)=rufrc(i,j)*cff7
2527 ru(i,j,0,nstp)=rufrc(i,j)
2528 END IF
2529# endif
2530#endif
2531#if defined NESTING && !defined SOLVE3D
2532 du_flux(i,j)=ubar(i,j,knew)* &
2533 & 0.5_r8*(dnew(i,j)+dnew(i-1,j))*on_u(i,j)
2534#endif
2535 END DO
2536 END DO
2537 DO j=jstrv,jend
2538 DO i=istr,iend
2539 cff=(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
2540 fac=1.0_r8/(dnew(i,j)+dnew(i,j-1))
2541 vbar(i,j,knew)=(vbar(i,j,kstp)* &
2542 & (dstp(i,j)+dstp(i,j-1))+ &
2543 & cff*cff1*rhs_vbar(i,j))*fac
2544#ifdef MASKING
2545 vbar(i,j,knew)=vbar(i,j,knew)*vmask(i,j)
2546#endif
2547#ifdef WET_DRY
2548 cff5=abs(abs(vmask_wet(i,j))-1.0_r8)
2549 cff6=0.5_r8+dsign(0.5_r8,vbar(i,j,knew))*vmask_wet(i,j)
2550 cff7=0.5_r8*vmask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
2551 vbar(i,j,knew)=vbar(i,j,knew)*cff7
2552 rhs_vbar(i,j)=rhs_vbar(i,j)*cff7
2553# ifdef SOLVE3D
2555 rvfrc(i,j)=rvfrc(i,j)*cff7
2556 rv(i,j,0,nstp)=rvfrc(i,j)
2557 END IF
2558# endif
2559#endif
2560#if defined NESTING && !defined SOLVE3D
2561 dv_flux(i,j)=vbar(i,j,knew)* &
2562 & 0.5_r8*(dnew(i,j)+dnew(i,j-1))*om_v(i,j)
2563#endif
2564 END DO
2565 END DO
2568#ifdef WET_DRY
2569 cff2=1.0_r8/cff1
2570#endif
2571 DO j=jstr,jend
2572 DO i=istru,iend
2573 cff=(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))
2574 fac=1.0_r8/(dnew(i,j)+dnew(i-1,j))
2575 ubar(i,j,knew)=(ubar(i,j,kstp)* &
2576 & (dstp(i,j)+dstp(i-1,j))+ &
2577 & cff*cff1*rhs_ubar(i,j))*fac
2578#ifdef MASKING
2579 ubar(i,j,knew)=ubar(i,j,knew)*umask(i,j)
2580#endif
2581#ifdef WET_DRY
2582 cff5=abs(abs(umask_wet(i,j))-1.0_r8)
2583 cff6=0.5_r8+dsign(0.5_r8,ubar(i,j,knew))*umask_wet(i,j)
2584 cff7=0.5_r8*umask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
2585 ubar(i,j,knew)=ubar(i,j,knew)*cff7
2586 rhs_ubar(i,j)=rhs_ubar(i,j)*cff7
2587#endif
2588#if defined NESTING && !defined SOLVE3D
2589 du_flux(i,j)=ubar(i,j,knew)* &
2590 & 0.5_r8*(dnew(i,j)+dnew(i-1,j))*on_u(i,j)
2591#endif
2592 END DO
2593 END DO
2594 DO j=jstrv,jend
2595 DO i=istr,iend
2596 cff=(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
2597 fac=1.0_r8/(dnew(i,j)+dnew(i,j-1))
2598 vbar(i,j,knew)=(vbar(i,j,kstp)* &
2599 & (dstp(i,j)+dstp(i,j-1))+ &
2600 & cff*cff1*rhs_vbar(i,j))*fac
2601#ifdef MASKING
2602 vbar(i,j,knew)=vbar(i,j,knew)*vmask(i,j)
2603#endif
2604#ifdef WET_DRY
2605 cff5=abs(abs(vmask_wet(i,j))-1.0_r8)
2606 cff6=0.5_r8+dsign(0.5_r8,vbar(i,j,knew))*vmask_wet(i,j)
2607 cff7=0.5_r8*vmask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
2608 vbar(i,j,knew)=vbar(i,j,knew)*cff7
2609 rhs_vbar(i,j)=rhs_vbar(i,j)*cff7
2610#endif
2611#if defined NESTING && !defined SOLVE3D
2612 dv_flux(i,j)=vbar(i,j,knew)* &
2613 & 0.5_r8*(dnew(i,j)+dnew(i,j-1))*om_v(i,j)
2614#endif
2615 END DO
2616 END DO
2617 ELSE IF (corrector_2d_step) THEN
2618 cff1=0.5_r8*
dtfast(ng)*5.0_r8/12.0_r8
2619 cff2=0.5_r8*
dtfast(ng)*8.0_r8/12.0_r8
2620 cff3=0.5_r8*
dtfast(ng)*1.0_r8/12.0_r8
2621#ifdef WET_DRY
2622 cff4=1.0_r8/cff1
2623#endif
2624 DO j=jstr,jend
2625 DO i=istru,iend
2626 cff=(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))
2627 fac=1.0_r8/(dnew(i,j)+dnew(i-1,j))
2628 ubar(i,j,knew)=(ubar(i,j,kstp)* &
2629 & (dstp(i,j)+dstp(i-1,j))+ &
2630 & cff*(cff1*rhs_ubar(i,j)+ &
2631 & cff2*rubar(i,j,kstp)- &
2632 & cff3*rubar(i,j,ptsk)))*fac
2633#ifdef MASKING
2634 ubar(i,j,knew)=ubar(i,j,knew)*umask(i,j)
2635#endif
2636#ifdef WET_DRY
2637 cff5=abs(abs(umask_wet(i,j))-1.0_r8)
2638 cff6=0.5_r8+dsign(0.5_r8,ubar(i,j,knew))*umask_wet(i,j)
2639 cff7=0.5_r8*umask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
2640 ubar(i,j,knew)=ubar(i,j,knew)*cff7
2641 rhs_ubar(i,j)=rhs_ubar(i,j)*cff7
2642#endif
2643#if defined NESTING && !defined SOLVE3D
2644 du_flux(i,j)=ubar(i,j,knew)* &
2645 & 0.5_r8*(dnew(i,j)+dnew(i-1,j))*on_u(i,j)
2646#endif
2647 END DO
2648 END DO
2649 DO j=jstrv,jend
2650 DO i=istr,iend
2651 cff=(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
2652 fac=1.0_r8/(dnew(i,j)+dnew(i,j-1))
2653 vbar(i,j,knew)=(vbar(i,j,kstp)* &
2654 & (dstp(i,j)+dstp(i,j-1))+ &
2655 & cff*(cff1*rhs_vbar(i,j)+ &
2656 & cff2*rvbar(i,j,kstp)- &
2657 & cff3*rvbar(i,j,ptsk)))*fac
2658#ifdef MASKING
2659 vbar(i,j,knew)=vbar(i,j,knew)*vmask(i,j)
2660#endif
2661#ifdef WET_DRY
2662 cff5=abs(abs(vmask_wet(i,j))-1.0_r8)
2663 cff6=0.5_r8+dsign(0.5_r8,vbar(i,j,knew))*vmask_wet(i,j)
2664 cff7=0.5_r8*vmask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
2665 vbar(i,j,knew)=vbar(i,j,knew)*cff7
2666 rhs_vbar(i,j)=rhs_vbar(i,j)*cff7
2667#endif
2668#if defined NESTING && !defined SOLVE3D
2669 dv_flux(i,j)=vbar(i,j,knew)* &
2670 & 0.5_r8*(dnew(i,j)+dnew(i,j-1))*om_v(i,j)
2671#endif
2672 END DO
2673 END DO
2674 END IF
2675
2676#ifdef DIAGNOSTICS_UV
2677
2678
2679
2680
2681
2682# ifdef MASKING
2683
2684
2685
2687 DO j=jstr,jend
2688 DO i=istru,iend
2689 diau2rhs(i,j,idiag)=diau2rhs(i,j,idiag)*umask(i,j)
2690 END DO
2691 END DO
2692 DO j=jstrv,jend
2693 DO i=istr,iend
2694 diav2rhs(i,j,idiag)=diav2rhs(i,j,idiag)*vmask(i,j)
2695 END DO
2696 END DO
2697 END DO
2698# endif
2699# ifdef SOLVE3D
2700
2701
2702
2703
2704
2705
2706
2708 IF (first_2d_step.and.corrector_2d_step) THEN
2711 DO j=jstr,jend
2712 DO i=istru,iend
2713 diau2int(i,j,idiag)=cff1*diau2rhs(i,j,idiag)
2714 diau2wrk(i,j,idiag)=diau2int(i,j,idiag)* &
2715 & (pm(i-1,j)+pm(i,j))*fac
2716 END DO
2717 END DO
2718 DO j=jstrv,jend
2719 DO i=istr,iend
2720 diav2int(i,j,idiag)=cff1*diav2rhs(i,j,idiag)
2721 diav2wrk(i,j,idiag)=diav2int(i,j,idiag)* &
2722 & (pn(i,j)+pn(i,j-1))*fac
2723 END DO
2724 END DO
2725 END DO
2726 ELSE IF (corrector_2d_step) THEN
2727 cff1=0.5_r8*
dtfast(ng)*5.0_r8/12.0_r8
2728 cff2=0.5_r8*
dtfast(ng)*8.0_r8/12.0_r8
2729 cff3=0.5_r8*
dtfast(ng)*1.0_r8/12.0_r8
2731 DO j=jstr,jend
2732 DO i=istru,iend
2733 diau2int(i,j,idiag)=diau2int(i,j,idiag)+ &
2734 & (cff1*diau2rhs(i,j,idiag)+ &
2735 & cff2*diarubar(i,j,kstp,idiag)- &
2736 & cff3*diarubar(i,j,ptsk,idiag))
2737 diau2wrk(i,j,idiag)=diau2wrk(i,j,idiag)+ &
2738 & diau2int(i,j,idiag)* &
2739 & (pm(i-1,j)+pm(i,j))*fac
2740 END DO
2741 END DO
2742 DO j=jstrv,jend
2743 DO i=istr,iend
2744 diav2int(i,j,idiag)=diav2int(i,j,idiag)+ &
2745 & (cff1*diav2rhs(i,j,idiag)+ &
2746 & cff2*diarvbar(i,j,kstp,idiag)- &
2747 & cff3*diarvbar(i,j,ptsk,idiag))
2748 diav2wrk(i,j,idiag)=diav2wrk(i,j,idiag)+ &
2749 & diav2int(i,j,idiag)* &
2750 & (pn(i,j)+pn(i,j-1))*fac
2751 END DO
2752 END DO
2753 END DO
2754 END IF
2755# else
2756
2757
2758
2759 IF (first_2d_step.and.corrector_2d_step) THEN
2762 DO j=jstr,jend
2763 DO i=istru,iend
2764 cff=(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))
2765 fac=1.0_r8/(dnew(i,j)+dnew(i-1,j))
2766 diau2wrk(i,j,idiag)=cff*cff1*diau2rhs(i,j,idiag)*fac
2767 END DO
2768 END DO
2769 DO j=jstrv,jend
2770 DO i=istr,iend
2771 cff=(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
2772 fac=1.0_r8/(dnew(i,j)+dnew(i,j-1))
2773 diav2wrk(i,j,idiag)=cff*cff1*diav2rhs(i,j,idiag)*fac
2774 END DO
2775 END DO
2776 END DO
2777 DO j=jstr,jend
2778 DO i=istru,iend
2779 fac=1.0_r8/(dnew(i,j)+dnew(i-1,j))
2780 diau2wrk(i,j,
m2rate)=ubar(i,j,knew)-ubar(i,j,kstp)* &
2781 & (dstp(i,j)+dstp(i-1,j))*fac
2782 END DO
2783 END DO
2784 DO j=jstrv,jend
2785 DO i=istr,iend
2786 fac=1.0_r8/(dnew(i,j)+dnew(i,j-1))
2787 diav2wrk(i,j,
m2rate)=vbar(i,j,knew)-vbar(i,j,kstp)* &
2788 & (dstp(i,j)+dstp(i,j-1))*fac
2789 END DO
2790 END DO
2791 ELSE IF (corrector_2d_step) THEN
2792 cff1=0.5_r8*
dtfast(ng)*5.0_r8/12.0_r8
2793 cff2=0.5_r8*
dtfast(ng)*8.0_r8/12.0_r8
2794 cff3=0.5_r8*
dtfast(ng)*1.0_r8/12.0_r8
2796 DO j=jstr,jend
2797 DO i=istru,iend
2798 cff=(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))
2799 fac=1.0_r8/(dnew(i,j)+dnew(i-1,j))
2800 diau2wrk(i,j,idiag)=cff*(cff1*diau2rhs(i,j,idiag)+ &
2801 & cff2*diarubar(i,j,kstp,idiag)- &
2802 & cff3*diarubar(i,j,ptsk,idiag))* &
2803 & fac
2804 END DO
2805 END DO
2806 DO j=jstrv,jend
2807 DO i=istr,iend
2808 cff=(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
2809 fac=1.0_r8/(dnew(i,j)+dnew(i,j-1))
2810 diav2wrk(i,j,idiag)=cff*(cff1*diav2rhs(i,j,idiag)+ &
2811 & cff2*diarvbar(i,j,kstp,idiag)- &
2812 & cff3*diarvbar(i,j,ptsk,idiag))* &
2813 & fac
2814 END DO
2815 END DO
2816 END DO
2817 DO j=jstr,jend
2818 DO i=istru,iend
2819 fac=1.0_r8/(dnew(i,j)+dnew(i-1,j))
2820 diau2wrk(i,j,
m2rate)=ubar(i,j,knew)- &
2821 & ubar(i,j,kstp)* &
2822 & (dstp(i,j)+dstp(i-1,j))*fac
2823 END DO
2824 END DO
2825 DO j=jstrv,jend
2826 DO i=istr,iend
2827 fac=1.0_r8/(dnew(i,j)+dnew(i,j-1))
2828 diav2wrk(i,j,
m2rate)=vbar(i,j,knew)- &
2829 & vbar(i,j,kstp)* &
2830 & (dstp(i,j)+dstp(i,j-1))*fac
2831 END DO
2832 END DO
2833 END IF
2834# endif
2835#endif
2836
2837
2838
2839
2841 DO j=jstr,jend
2842 DO i=istru,iend
2843 rubar(i,j,krhs)=rhs_ubar(i,j)
2844 END DO
2845 END DO
2846 DO j=jstrv,jend
2847 DO i=istr,iend
2848 rvbar(i,j,krhs)=rhs_vbar(i,j)
2849 END DO
2850 END DO
2851#ifdef DIAGNOSTICS_UV
2853 DO j=jstr,jend
2854 DO i=istru,iend
2855 diarubar(i,j,krhs,idiag)=diau2rhs(i,j,idiag)
2856 END DO
2857 END DO
2858 DO j=jstrv,jend
2859 DO i=istr,iend
2860 diarvbar(i,j,krhs,idiag)=diav2rhs(i,j,idiag)
2861 END DO
2862 END DO
2863 END DO
2864#endif
2865 END IF
2866
2867
2868
2869
2870
2872 & lbi, ubi, lbj, ubj, &
2873 & imins, imaxs, jmins, jmaxs, &
2874 & krhs, kstp, knew, &
2875 & ubar, vbar, zeta)
2877 & lbi, ubi, lbj, ubj, &
2878 & imins, imaxs, jmins, jmaxs, &
2879 & krhs, kstp, knew, &
2880 & ubar, vbar, zeta)
2881
2882
2883
2884
2887 & lbi, ubi, lbj, ubj, &
2888 & imins, imaxs, jmins, jmaxs, &
2889 & knew, &
2890#ifdef MASKING
2891 & umask, vmask, &
2892#endif
2893 & h, om_v, on_u, &
2894 & ubar, vbar, zeta)
2895 END IF
2896
2897#if defined NESTING && !defined SOLVE3D
2898
2899
2900
2901
2902
2904 IF (
domain(ng)%Western_Edge(tile))
THEN
2905 DO j=jstr-1,jendr
2906 dnew(istr-1,j)=h(istr-1,j)+zeta_new(istr-1,j)
2907 END DO
2908 END IF
2909 END IF
2911 IF (
domain(ng)%Eastern_Edge(tile))
THEN
2912 DO j=jstr-1,jendr
2913 dnew(iend+1,j)=h(iend+1,j)+zeta_new(iend+1,j)
2914 END DO
2915 END IF
2916 END IF
2918 IF (
domain(ng)%Southern_Edge(tile))
THEN
2919 DO i=istr-1,iendr
2920 dnew(i,jstr-1)=h(i,jstr-1)+zeta_new(i,jstr-1)
2921 END DO
2922 END IF
2923 END IF
2925 IF (
domain(ng)%Northern_Edge(tile))
THEN
2926 DO i=istr-1,iendr
2927 dnew(i,jend+1)=h(i,jend+1)+zeta_new(i,jend+1)
2928 END DO
2929 END IF
2930 END IF
2931
2933 IF (
domain(ng)%Western_Edge(tile))
THEN
2934 DO j=jstrr,jendr
2935 du_flux(istru-1,j)=ubar(istru-1,j,knew)*on_u(istru-1,j)* &
2936 & 0.5_r8*(dnew(istru-1,j)+dnew(istru-2,j))
2937 END DO
2938 DO j=jstrv,jend
2939 dv_flux(istr-1,j)=vbar(istr-1,j,knew)*om_v(istr-1,j)* &
2940 & 0.5_r8*(dnew(istr-1,j)+dnew(istr-1,j-1))
2941 END DO
2942 END IF
2943 END IF
2945 IF (
domain(ng)%Eastern_Edge(tile))
THEN
2946 DO j=jstrr,jendr
2947 du_flux(iend+1,j)=ubar(iend+1,j,knew)*on_u(iend+1,j)* &
2948 & 0.5_r8*(dnew(iend+1,j)+dnew(iend,j))
2949 END DO
2950 DO j=jstrv,jend
2951 dv_flux(iend+1,j)=vbar(iend+1,j,knew)*om_v(iend+1,j)* &
2952 & 0.5_r8*(dnew(iend+1,j)+dnew(iend+1,j-1))
2953 END DO
2954 END IF
2955 END IF
2957 IF (
domain(ng)%Southern_Edge(tile))
THEN
2958 DO i=istru,iend
2959 du_flux(i,jstr-1)=ubar(i,jstr-1,knew)*on_u(i,jstr-1)* &
2960 & 0.5_r8*(dnew(i,jstr-1)+dnew(i-1,jstr-1))
2961 END DO
2962 DO i=istrr,iendr
2963 dv_flux(i,jstrv-1)=vbar(i,jstrv-1,knew)*om_v(i,jstrv-1)* &
2964 & 0.5_r8*(dnew(i,jstrv-1)+dnew(i,jstrv-2))
2965 END DO
2966 END IF
2967 END IF
2969 IF (
domain(ng)%Northern_Edge(tile))
THEN
2970 DO i=istru,iend
2971 du_flux(i,jend+1)=ubar(i,jend+1,knew)*on_u(i,jend+1)* &
2972 & 0.5_r8*(dnew(i,jend+1)+dnew(i-1,jend+1))
2973 END DO
2974 DO i=istrr,iendr
2975 dv_flux(i,jend+1)=vbar(i,jend+1,knew)*om_v(i,jend+1)* &
2976 & 0.5_r8*(dnew(i,jend+1)+dnew(i,jend))
2977 END DO
2978 END IF
2979 END IF
2980#endif
2981
2982
2983
2984
2985
2986
2987
2988
2993 IF (((istrr.le.i).and.(i.le.iendr)).and. &
2994 & ((jstrr.le.j).and.(j.le.jendr))) THEN
2995 IF (int(
sources(ng)%Dsrc(is)).eq.0)
THEN
2996 cff=1.0_r8/(on_u(i,j)* &
2997 & 0.5_r8*(zeta(i-1,j,knew)+h(i-1,j)+ &
2998 & zeta(i ,j,knew)+h(i ,j)))
2999 ubar(i,j,knew)=
sources(ng)%Qbar(is)*cff
3000#if defined NESTING && !defined SOLVE3D
3001 du_flux(i,j)=
sources(ng)%Qbar(is)
3002#endif
3003 ELSE IF (int(
sources(ng)%Dsrc(is)).eq.1)
THEN
3004 cff=1.0_r8/(om_v(i,j)* &
3005 & 0.5_r8*(zeta(i,j-1,knew)+h(i,j-1)+ &
3006 & zeta(i,j ,knew)+h(i,j )))
3007 vbar(i,j,knew)=
sources(ng)%Qbar(is)*cff
3008#if defined NESTING && !defined SOLVE3D
3009 dv_flux(i,j)=
sources(ng)%Qbar(is)
3010#endif
3011 END IF
3012 END IF
3013 END DO
3014 END IF
3015
3016
3017
3018
3019
3022 & lbi, ubi, lbj, ubj, &
3023 & ubar(:,:,knew))
3025 & lbi, ubi, lbj, ubj, &
3026 & vbar(:,:,knew))
3027
3028#if defined NESTING && !defined SOLVE3D
3030 & lbi, ubi, lbj, ubj, &
3031 & du_flux)
3033 & lbi, ubi, lbj, ubj, &
3034 & dv_flux)
3035#endif
3036 END IF
3037
3038#ifdef DISTRIBUTE
3039
3040# if defined NESTING && !defined SOLVE3D
3042# else
3044# endif
3045 & lbi, ubi, lbj, ubj, &
3048# if defined NESTING && !defined SOLVE3D
3049 & du_flux, dv_flux, &
3050# endif
3051 & ubar(:,:,knew), &
3052 & vbar(:,:,knew))
3053#endif
3054
3055 RETURN
subroutine exchange_r2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
subroutine exchange_u2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
subroutine exchange_v2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
type(t_clima), dimension(:), allocatable clima
type(t_lbc), dimension(:,:,:), allocatable lbc
type(t_domain), dimension(:), allocatable domain
logical, dimension(:), allocatable luvsrc
logical, dimension(:), allocatable lnudgem2clm
integer, dimension(:), allocatable iic
real(dp), dimension(:), allocatable dt
real(r8), dimension(:), allocatable dcrit
logical, dimension(:), allocatable ewperiodic
logical, dimension(:), allocatable nsperiodic
logical, dimension(:), allocatable predictor_2d_step
real(dp), dimension(:,:,:), allocatable weight
logical, dimension(:,:), allocatable volcons
integer, dimension(:), allocatable nfast
integer, dimension(:), allocatable ndtfast
logical, dimension(:), allocatable lwsrc
real(r8), dimension(:), allocatable gamma2
logical, dimension(:,:), allocatable compositegrid
integer, parameter isouth
real(dp), dimension(:), allocatable dtfast
integer, dimension(:), allocatable ntfirst
integer, parameter inorth
integer, dimension(:), allocatable iif
type(t_sources), dimension(:), allocatable sources
integer, dimension(:), allocatable nsrc
subroutine mp_exchange2d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, nghost, ew_periodic, ns_periodic, a, b, c, d)
subroutine, public set_duv_bc_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, kinp, umask, vmask, om_v, on_u, ubar, vbar, drhs, duon, dvom)
subroutine, public obc_flux_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, kinp, umask, vmask, h, om_v, on_u, ubar, vbar, zeta)
subroutine, public u2dbc_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, krhs, kstp, kout, ubar, vbar, zeta)
subroutine, public v2dbc_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, krhs, kstp, kout, ubar, vbar, zeta)
subroutine wetdry_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, pmask, rmask, umask, vmask, h, zeta, du_avg1, dv_avg1, rmask_wet_avg, pmask_wet, pmask_full, rmask_wet, rmask_full, umask_wet, umask_full, vmask_wet, vmask_full)
subroutine, public zetabc_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, krhs, kstp, kout, zeta)