- Timestamp:
- 2020-07-02T17:50:26+02:00 (4 years ago)
- Location:
- NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src
- Files:
-
- 3 edited
-
NST/agrif_oce_interp.F90 (modified) (43 diffs)
-
NST/agrif_oce_sponge.F90 (modified) (9 diffs)
-
OCE/DYN/sshwzv.F90 (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/NST/agrif_oce_interp.F90
r13229 r13230 34 34 USE lib_mpp 35 35 USE vremap 36 USE lbclnk 36 37 37 38 IMPLICIT NONE … … 44 45 PUBLIC interpunb, interpvnb , interpub2b, interpvb2b 45 46 PUBLIC interpe3t, interpglamt, interpgphit 46 #if defined key_vertical47 47 PUBLIC interpht0, interpmbkt 48 # endif 48 PUBLIC agrif_initts, agrif_initssh 49 49 50 INTEGER :: bdy_tinterp = 0 50 51 … … 89 90 Agrif_UseSpecialValue = ln_spc_dyn 90 91 ! 92 use_sign_north = .TRUE. 93 sign_north = -1. 91 94 CALL Agrif_Bc_variable( un_interp_id, procname=interpun ) 92 95 CALL Agrif_Bc_variable( vn_interp_id, procname=interpvn ) 96 use_sign_north = .FALSE. 93 97 ! 94 98 Agrif_UseSpecialValue = .FALSE. … … 99 103 ibdy2 = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells 100 104 ! 101 IF( .NOT.ln_dynspg_ts ) THEN ! Store t angential transport105 IF( .NOT.ln_dynspg_ts ) THEN ! Store transport 102 106 DO ji = mi0(ibdy1), mi1(ibdy2) 103 107 uu_b(ji,:,Krhs_a) = 0._wp … … 116 120 ! 117 121 DO ji = mi0(ibdy1), mi1(ibdy2) 118 zub(ji,:) = 0._wp ! Correct t angential transport122 zub(ji,:) = 0._wp ! Correct transport 119 123 DO jk = 1, jpkm1 120 124 DO jj = 1, jpj … … 126 130 zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 127 131 END DO 128 132 129 133 DO jk = 1, jpkm1 130 134 DO jj = 1, jpj … … 133 137 END DO 134 138 END DO 135 139 136 140 IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate 137 141 DO ji = mi0(ibdy1), mi1(ibdy2) … … 155 159 156 160 ! --- East --- ! 157 IF( lk_east 161 IF( lk_east) THEN 158 162 ibdy1 = jpiglo - ( nn_hls + nbghostcells + 1) ! halo + land + nbghostcells 159 163 ibdy2 = jpiglo - ( nn_hls + 2 ) ! halo + land + 1 … … 185 189 zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 186 190 END DO 187 191 188 192 DO jk = 1, jpkm1 189 193 DO jj = 1, jpj … … 193 197 END DO 194 198 END DO 195 199 196 200 IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate 197 ibdy1 = jpiglo - ( nn_hls + nbghostcells ) ! halo + land + nbghostcells - 1198 ibdy2 = jpiglo - ( nn_hls + 1 ) ! halo + land + 1 - 1201 ibdy1 = jpiglo-nbghostcells 202 ibdy2 = jpiglo-1 199 203 DO ji = mi0(ibdy1), mi1(ibdy2) 200 204 zvb(ji,:) = 0._wp … … 257 261 END DO 258 262 END DO 259 263 260 264 IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate 261 265 DO jj = mj0(jbdy1), mj1(jbdy2) … … 270 274 zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 271 275 END DO 272 276 273 277 DO jk = 1, jpkm1 274 278 DO ji = 1, jpi 275 279 uu(ji,jj,jk,Krhs_a) = ( uu(ji,jj,jk,Krhs_a) & 276 280 & + uu_b(ji,jj,Krhs_a) - zub(ji,jj) ) * umask(ji,jj,jk) 277 281 END DO 278 282 END DO … … 316 320 DO ji = 1, jpi 317 321 vv(ji,jj,jk,Krhs_a) = ( vv(ji,jj,jk,Krhs_a) & 318 319 END DO 320 END DO 321 END DO 322 322 & + vv_b(ji,jj,Krhs_a) - zvb(ji,jj) ) * vmask(ji,jj,jk) 323 END DO 324 END DO 325 END DO 326 323 327 IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate 324 jbdy1 = jpjglo - ( nn_hls + nbghostcells ) ! halo + land + nbghostcells - 1325 jbdy2 = jpjglo - ( nn_hls + 1 ) ! halo + land + 1 -1328 jbdy1 = jpjglo-nbghostcells 329 jbdy2 = jpjglo-1 326 330 DO jj = mj0(jbdy1), mj1(jbdy2) 327 331 zub(:,jj) = 0._wp … … 335 339 zub(ji,jj) = zub(ji,jj) * r1_hu(ji,jj,Krhs_a) 336 340 END DO 337 341 338 342 DO jk = 1, jpkm1 339 343 DO ji = 1, jpi … … 378 382 iend = jpiglo - ( nn_hls + 1 ) ! halo + land + 1 - 1 379 383 DO ji = mi0(istart), mi1(iend) 384 380 385 DO jj=1,jpj 381 386 va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj) … … 389 394 END DO 390 395 END DO 391 ENDIF 396 ENDIF 392 397 ! 393 398 !--- South ---! … … 395 400 jstart = nn_hls + 2 ! halo + land + 1 396 401 jend = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells 402 jstart = 2 403 jend = nbghostcells+1 397 404 DO jj = mj0(jstart), mj1(jend) 405 398 406 DO ji=1,jpi 399 407 ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj) … … 401 409 END DO 402 410 END DO 403 ENDIF 411 ENDIF 404 412 ! 405 413 !--- North ---! … … 419 427 END DO 420 428 END DO 421 ENDIF 429 ENDIF 422 430 ! 423 431 END SUBROUTINE Agrif_dyn_ts 424 432 433 425 434 SUBROUTINE Agrif_dyn_ts_flux( jn, zu, zv ) 426 435 !!---------------------------------------------------------------------- … … 498 507 END SUBROUTINE Agrif_dyn_ts_flux 499 508 509 500 510 SUBROUTINE Agrif_dta_ts( kt ) 501 511 !!---------------------------------------------------------------------- … … 513 523 ! 514 524 ! Enforce volume conservation if no time refinement: 515 IF 525 IF ( Agrif_rhot()==1 ) ll_int_cons=.TRUE. 516 526 ! 517 527 ! Interpolate barotropic fluxes 518 528 Agrif_SpecialValue = 0._wp 519 529 Agrif_UseSpecialValue = ln_spc_dyn 530 531 use_sign_north = .TRUE. 532 sign_north = -1. 533 520 534 ! 521 535 ! Set bdy time interpolation stage to 0 (latter incremented locally do deal with corners) … … 542 556 ENDIF 543 557 Agrif_UseSpecialValue = .FALSE. 558 use_sign_north = .FALSE. 544 559 ! 545 560 END SUBROUTINE Agrif_dta_ts … … 566 581 ! 567 582 ! --- West --- ! 568 IF( lk_west) THEN569 istart = nn_hls + 2 ! halo + land + 1583 IF(lk_west) THEN 584 istart = nn_hls + 2 ! halo + land + 1 570 585 iend = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells 571 586 DO ji = mi0(istart), mi1(iend) 572 587 DO jj = 1, jpj 573 588 ssh(ji,jj,Krhs_a) = hbdy(ji,jj) 574 END DO575 END DO589 END DO 590 END DO 576 591 ENDIF 577 592 ! 578 593 ! --- East --- ! 579 IF( lk_east) THEN594 IF(lk_east) THEN 580 595 istart = jpiglo - ( nn_hls + nbghostcells ) ! halo + land + nbghostcells - 1 581 596 iend = jpiglo - ( nn_hls + 1 ) ! halo + land + 1 - 1 … … 583 598 DO jj = 1, jpj 584 599 ssh(ji,jj,Krhs_a) = hbdy(ji,jj) 585 END DO586 END DO600 END DO 601 END DO 587 602 ENDIF 588 603 ! 589 604 ! --- South --- ! 590 IF( lk_south) THEN605 IF(lk_south) THEN 591 606 jstart = nn_hls + 2 ! halo + land + 1 592 607 jend = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells … … 594 609 DO ji = 1, jpi 595 610 ssh(ji,jj,Krhs_a) = hbdy(ji,jj) 596 END DO597 END DO611 END DO 612 END DO 598 613 ENDIF 599 614 ! 600 615 ! --- North --- ! 601 IF( lk_north) THEN616 IF(lk_north) THEN 602 617 jstart = jpjglo - ( nn_hls + nbghostcells ) ! halo + land + nbghostcells - 1 603 618 jend = jpjglo - ( nn_hls + 1 ) ! halo + land + 1 - 1 … … 605 620 DO ji = 1, jpi 606 621 ssh(ji,jj,Krhs_a) = hbdy(ji,jj) 607 END DO608 END DO622 END DO 623 END DO 609 624 ENDIF 610 625 ! … … 625 640 ! 626 641 ! --- West --- ! 627 IF( lk_west) THEN642 IF(lk_west) THEN 628 643 istart = nn_hls + 2 ! halo + land + 1 629 644 iend = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells … … 631 646 DO jj = 1, jpj 632 647 ssha_e(ji,jj) = hbdy(ji,jj) 633 END DO634 END DO648 END DO 649 END DO 635 650 ENDIF 636 651 ! 637 652 ! --- East --- ! 638 IF( lk_east) THEN653 IF(lk_east) THEN 639 654 istart = jpiglo - ( nn_hls + nbghostcells ) ! halo + land + nbghostcells - 1 640 655 iend = jpiglo - ( nn_hls + 1 ) ! halo + land + 1 - 1 … … 642 657 DO jj = 1, jpj 643 658 ssha_e(ji,jj) = hbdy(ji,jj) 644 END DO645 END DO659 END DO 660 END DO 646 661 ENDIF 647 662 ! 648 663 ! --- South --- ! 649 IF( lk_south) THEN664 IF(lk_south) THEN 650 665 jstart = nn_hls + 2 ! halo + land + 1 651 666 jend = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells … … 653 668 DO ji = 1, jpi 654 669 ssha_e(ji,jj) = hbdy(ji,jj) 655 END DO656 END DO670 END DO 671 END DO 657 672 ENDIF 658 673 ! 659 674 ! --- North --- ! 660 IF( lk_north) THEN675 IF(lk_north) THEN 661 676 jstart = jpjglo - ( nn_hls + nbghostcells ) ! halo + land + nbghostcells - 1 662 677 jend = jpjglo - ( nn_hls + 1 ) ! halo + land + 1 - 1 … … 664 679 DO ji = 1, jpi 665 680 ssha_e(ji,jj) = hbdy(ji,jj) 666 END DO667 END DO681 END DO 682 END DO 668 683 ENDIF 669 684 ! 670 685 END SUBROUTINE Agrif_ssh_ts 671 686 687 672 688 SUBROUTINE Agrif_avm 673 689 !!---------------------------------------------------------------------- … … 690 706 ! 691 707 END SUBROUTINE Agrif_avm 692 708 693 709 694 710 SUBROUTINE interptsn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before ) … … 702 718 INTEGER :: ji, jj, jk, jn ! dummy loop indices 703 719 INTEGER :: N_in, N_out 720 INTEGER :: item 704 721 ! vertical interpolation: 705 722 REAL(wp) :: zhtot 706 723 REAL(wp), DIMENSION(k1:k2,1:jpts) :: tabin 707 REAL(wp), DIMENSION(k1:k2) :: h_in 708 REAL(wp), DIMENSION(1:jpk) :: h_out 709 !!---------------------------------------------------------------------- 710 711 IF( before ) THEN 724 REAL(wp), DIMENSION(k1:k2) :: h_in, z_in 725 REAL(wp), DIMENSION(1:jpk) :: h_out, z_out 726 !!---------------------------------------------------------------------- 727 728 IF( before ) THEN 729 730 item = Kmm_a 731 IF( l_ini_child ) Kmm_a = Kbb_a 732 712 733 DO jn = 1,jpts 713 734 DO jk=k1,k2 … … 718 739 END DO 719 740 END DO 720 END DO 721 722 # if defined key_vertical 723 ! Interpolate thicknesses 724 ! Warning: these are masked, hence extrapolated prior interpolation. 725 DO jk=k1,k2 726 DO jj=j1,j2 727 DO ji=i1,i2 728 ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 729 END DO 730 END DO 731 END DO 732 733 ! Extrapolate thicknesses in partial bottom cells: 734 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 735 IF (ln_zps) THEN 736 DO jj=j1,j2 737 DO ji=i1,i2 738 jk = mbkt(ji,jj) 739 ptab(ji,jj,jk,jpts+1) = 0._wp 740 END DO 741 END DO 742 END IF 743 744 ! Save ssh at last level: 745 IF (.NOT.ln_linssh) THEN 746 ptab(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 747 ELSE 748 ptab(i1:i2,j1:j2,k2,jpts+1) = 0._wp 749 END IF 750 # endif 741 END DO 742 743 IF( l_vremap .OR. l_ini_child) THEN 744 ! Interpolate thicknesses 745 ! Warning: these are masked, hence extrapolated prior interpolation. 746 DO jk=k1,k2 747 DO jj=j1,j2 748 DO ji=i1,i2 749 ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 750 751 END DO 752 END DO 753 END DO 754 755 ! Extrapolate thicknesses in partial bottom cells: 756 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 757 IF (ln_zps) THEN 758 DO jj=j1,j2 759 DO ji=i1,i2 760 jk = mbkt(ji,jj) 761 ptab(ji,jj,jk,jpts+1) = 0._wp 762 END DO 763 END DO 764 END IF 765 766 ! Save ssh at last level: 767 IF (.NOT.ln_linssh) THEN 768 ptab(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 769 ELSE 770 ptab(i1:i2,j1:j2,k2,jpts+1) = 0._wp 771 END IF 772 ENDIF 773 Kmm_a = item 774 751 775 ELSE 752 753 # if defined key_vertical 754 IF (ln_linssh) ptab(i1:i2,j1:j2,k2,n2) = 0._wp 755 756 DO jj=j1,j2 757 DO ji=i1,i2 758 ts(ji,jj,:,:,Krhs_a) = 0._wp 759 N_in = mbkt_parent(ji,jj) 760 zhtot = 0._wp 761 DO jk=1,N_in !k2 = jpk of parent grid 762 IF (jk==N_in) THEN 763 h_in(jk) = ht0_parent(ji,jj) + ptab(ji,jj,k2,n2) - zhtot 764 ELSE 765 h_in(jk) = ptab(ji,jj,jk,n2) 776 item = Krhs_a 777 IF( l_ini_child ) Krhs_a = Kbb_a 778 779 IF( l_vremap .OR. l_ini_child ) THEN 780 IF (ln_linssh) ptab(i1:i2,j1:j2,k2,n2) = 0._wp 781 782 DO jj=j1,j2 783 DO ji=i1,i2 784 ts(ji,jj,:,:,Krhs_a) = 0. 785 ! IF( l_ini_child) ts(ji,jj,:,:,Krhs_a) = ptab(ji,jj,:,1:jpts) 786 N_in = mbkt_parent(ji,jj) 787 zhtot = 0._wp 788 DO jk=1,N_in !k2 = jpk of parent grid 789 IF (jk==N_in) THEN 790 h_in(jk) = ht0_parent(ji,jj) + ptab(ji,jj,k2,n2) - zhtot 791 ELSE 792 h_in(jk) = ptab(ji,jj,jk,n2) 793 ENDIF 794 zhtot = zhtot + h_in(jk) 795 tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1) 796 END DO 797 z_in(1) = 0.5_wp * h_in(1) - zhtot + ht0_parent(ji,jj) 798 DO jk=2,N_in 799 z_in(jk) = z_in(jk-1) + 0.5_wp * h_in(jk) 800 END DO 801 802 N_out = 0 803 DO jk=1,jpk ! jpk of child grid 804 IF (tmask(ji,jj,jk) == 0._wp) EXIT 805 N_out = N_out + 1 806 h_out(jk) = e3t(ji,jj,jk,Krhs_a) 807 END DO 808 809 z_out(1) = 0.5_wp * h_out(1) - SUM(h_out(1:N_out)) + ht_0(ji,jj) 810 DO jk=2,N_out 811 z_out(jk) = z_out(jk-1) + 0.5_wp * h_out(jk) 812 END DO 813 814 IF (N_in*N_out > 0) THEN 815 IF( l_ini_child ) THEN 816 CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a), & 817 & z_out(1:N_out),N_in,N_out,jpts) 818 ELSE 819 CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a), & 820 & h_out(1:N_out),N_in,N_out,jpts) 821 ENDIF 766 822 ENDIF 767 zhtot = zhtot + h_in(jk) 768 tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1) 769 END DO 770 N_out = 0 771 DO jk=1,jpk ! jpk of child grid 772 IF (tmask(ji,jj,jk) == 0._wp) EXIT 773 N_out = N_out + 1 774 h_out(jk) = e3t(ji,jj,jk,Krhs_a) 775 ENDDO 776 IF (N_in*N_out > 0) THEN 777 CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),ts(ji,jj,1:N_out,1:jpts,Krhs_a),h_out(1:N_out),N_in,N_out,jpts) 778 ENDIF 779 ENDDO 780 ENDDO 781 # else 782 ! 783 DO jn=1, jpts 784 ts(i1:i2,j1:j2,1:jpk,jn,Krhs_a)=ptab(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 785 END DO 786 # endif 823 END DO 824 END DO 825 Krhs_a = item 826 827 ELSE 828 829 DO jn=1, jpts 830 ts(i1:i2,j1:j2,1:jpk,jn,Krhs_a)=ptab(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 831 END DO 832 ENDIF 787 833 788 834 ENDIF … … 790 836 END SUBROUTINE interptsn 791 837 838 792 839 SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before ) 793 840 !!---------------------------------------------------------------------- … … 808 855 END SUBROUTINE interpsshn 809 856 857 810 858 SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before ) 811 859 !!---------------------------------------------------------------------- … … 820 868 REAL(wp) :: zrhoy, zhtot 821 869 ! vertical interpolation: 822 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 823 REAL(wp), DIMENSION(1:jpk) :: h_out 824 INTEGER :: N_in, N_out 870 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in, z_in 871 REAL(wp), DIMENSION(1:jpk) :: h_out, z_out 872 INTEGER :: N_in, N_out,item 825 873 REAL(wp) :: h_diff 826 874 !!--------------------------------------------- 827 875 ! 828 876 IF (before) THEN 877 878 item = Kmm_a 879 IF( l_ini_child ) Kmm_a = Kbb_a 880 829 881 DO jk=1,jpk 830 882 DO jj=j1,j2 831 883 DO ji=i1,i2 832 884 ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a)*umask(ji,jj,jk)) 833 # if defined key_vertical 834 ! Interpolate thicknesses (masked for subsequent extrapolation) 835 ptab(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) 836 # endif 837 END DO 838 END DO 839 END DO 840 # if defined key_vertical 885 IF( l_vremap .OR. l_ini_child) THEN 886 ! Interpolate thicknesses (masked for subsequent extrapolation) 887 ptab(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) 888 ENDIF 889 END DO 890 END DO 891 END DO 892 893 IF( l_vremap .OR. l_ini_child) THEN 841 894 ! Extrapolate thicknesses in partial bottom cells: 842 895 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 843 IF (ln_zps) THEN 844 DO jj=j1,j2 845 DO ji=i1,i2 846 jk = mbku(ji,jj) 847 ptab(ji,jj,jk,2) = 0._wp 848 END DO 849 END DO 850 END IF 851 ! Save ssh at last level: 852 ptab(i1:i2,j1:j2,k2,2) = 0._wp 853 IF (.NOT.ln_linssh) THEN 854 ! This vertical sum below should be replaced by the sea-level at U-points (optimization): 855 DO jk=1,jpk 856 ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3u(i1:i2,j1:j2,jk,Kmm_a) * umask(i1:i2,j1:j2,jk) 857 END DO 858 ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hu_0(i1:i2,j1:j2) 859 END IF 860 # endif 896 IF (ln_zps) THEN 897 DO jj=j1,j2 898 DO ji=i1,i2 899 jk = mbku(ji,jj) 900 ptab(ji,jj,jk,2) = 0._wp 901 END DO 902 END DO 903 END IF 904 905 ! Save ssh at last level: 906 ptab(i1:i2,j1:j2,k2,2) = 0._wp 907 IF (.NOT.ln_linssh) THEN 908 ! This vertical sum below should be replaced by the sea-level at U-points (optimization): 909 DO jk=1,jpk 910 ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3u(i1:i2,j1:j2,jk,Kmm_a) * umask(i1:i2,j1:j2,jk) 911 END DO 912 ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hu_0(i1:i2,j1:j2) 913 END IF 914 ENDIF 915 916 Kmm_a = item 861 917 ! 862 918 ELSE 863 919 zrhoy = Agrif_rhoy() 864 # if defined key_vertical 920 921 IF( l_vremap .OR. l_ini_child) THEN 865 922 ! VERTICAL REFINEMENT BEGIN 866 923 867 IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 868 869 DO ji=i1,i2 870 DO jj=j1,j2 871 uu(ji,jj,:,Krhs_a) = 0._wp 872 N_in = mbku_parent(ji,jj) 873 zhtot = 0._wp 874 DO jk=1,N_in 875 IF (jk==N_in) THEN 876 h_in(jk) = hu0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot 877 ELSE 878 h_in(jk) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 879 ENDIF 880 zhtot = zhtot + h_in(jk) 881 tabin(jk) = ptab(ji,jj,jk,1)/(e2u(ji,jj)*zrhoy*h_in(jk)) 882 ENDDO 883 884 N_out = 0 885 DO jk=1,jpk 886 if (umask(ji,jj,jk) == 0) EXIT 887 N_out = N_out + 1 888 h_out(N_out) = e3u(ji,jj,jk,Krhs_a) 889 ENDDO 890 IF (N_in*N_out > 0) THEN 891 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),uu(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out,1) 892 ENDIF 893 ENDDO 894 ENDDO 895 896 # else 897 DO jk = 1, jpkm1 898 DO jj=j1,j2 899 uu(i1:i2,jj,jk,Krhs_a) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u(i1:i2,jj,jk,Krhs_a) ) 900 END DO 901 END DO 902 # endif 924 IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 925 926 DO ji=i1,i2 927 DO jj=j1,j2 928 uu(ji,jj,:,Krhs_a) = 0._wp 929 N_in = mbku_parent(ji,jj) 930 zhtot = 0._wp 931 DO jk=1,N_in 932 IF (jk==N_in) THEN 933 h_in(jk) = hu0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot 934 ELSE 935 h_in(jk) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 936 ENDIF 937 zhtot = zhtot + h_in(jk) 938 IF( h_in(jk) .GT. 0. ) THEN 939 tabin(jk) = ptab(ji,jj,jk,1)/(e2u(ji,jj)*zrhoy*h_in(jk)) 940 ELSE 941 tabin(jk) = 0. 942 ENDIF 943 END DO 944 z_in(1) = 0.5_wp * h_in(1) - zhtot + hu0_parent(ji,jj) 945 DO jk=2,N_in 946 z_in(jk) = z_in(jk-1) + 0.5_wp * h_in(jk) 947 END DO 948 949 N_out = 0 950 DO jk=1,jpk 951 IF (umask(ji,jj,jk) == 0) EXIT 952 N_out = N_out + 1 953 h_out(N_out) = e3u(ji,jj,jk,Krhs_a) 954 END DO 955 956 z_out(1) = 0.5_wp * h_out(1) - SUM(h_out(1:N_out)) + hu_0(ji,jj) 957 DO jk=2,N_out 958 z_out(jk) = z_out(jk-1) + 0.5_wp * h_out(jk) 959 END DO 960 961 IF (N_in*N_out > 0) THEN 962 IF( l_ini_child ) THEN 963 CALL remap_linear (tabin(1:N_in),z_in(1:N_in),uu(ji,jj,1:N_out,Krhs_a),z_out(1:N_out),N_in,N_out,1) 964 ELSE 965 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),uu(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out,1) 966 ENDIF 967 ENDIF 968 END DO 969 END DO 970 ELSE 971 DO jk = 1, jpkm1 972 DO jj=j1,j2 973 uu(i1:i2,jj,jk,Krhs_a) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u(i1:i2,jj,jk,Krhs_a) ) 974 END DO 975 END DO 976 ENDIF 903 977 904 978 ENDIF … … 906 980 END SUBROUTINE interpun 907 981 982 908 983 SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before ) 909 984 !!---------------------------------------------------------------------- … … 918 993 REAL(wp) :: zrhox 919 994 ! vertical interpolation: 920 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 921 REAL(wp), DIMENSION(1:jpk) :: h_out 922 INTEGER :: N_in, N_out 995 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in, z_in 996 REAL(wp), DIMENSION(1:jpk) :: h_out, z_out 997 INTEGER :: N_in, N_out, item 923 998 REAL(wp) :: h_diff, zhtot 924 999 !!--------------------------------------------- 925 1000 ! 926 IF (before) THEN 1001 IF (before) THEN 1002 1003 item = Kmm_a 1004 IF( l_ini_child ) Kmm_a = Kbb_a 1005 927 1006 DO jk=k1,k2 928 1007 DO jj=j1,j2 929 1008 DO ji=i1,i2 930 1009 ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)*vmask(ji,jj,jk)) 931 # if defined key_vertical 932 ! Interpolate thicknesses (masked for subsequent extrapolation) 933 ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) 934 # endif 935 END DO 936 END DO 937 END DO 938 # if defined key_vertical 1010 IF( l_vremap .OR. l_ini_child) THEN 1011 ! Interpolate thicknesses (masked for subsequent extrapolation) 1012 ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) 1013 ENDIF 1014 END DO 1015 END DO 1016 END DO 1017 1018 IF( l_vremap .OR. l_ini_child) THEN 939 1019 ! Extrapolate thicknesses in partial bottom cells: 940 1020 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 941 IF (ln_zps) THEN 1021 IF (ln_zps) THEN 1022 DO jj=j1,j2 1023 DO ji=i1,i2 1024 jk = mbkv(ji,jj) 1025 ptab(ji,jj,jk,2) = 0._wp 1026 END DO 1027 END DO 1028 END IF 1029 ! Save ssh at last level: 1030 ptab(i1:i2,j1:j2,k2,2) = 0._wp 1031 IF (.NOT.ln_linssh) THEN 1032 ! This vertical sum below should be replaced by the sea-level at V-points (optimization): 1033 DO jk=1,jpk 1034 ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3v(i1:i2,j1:j2,jk,Kmm_a) * vmask(i1:i2,j1:j2,jk) 1035 END DO 1036 ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hv_0(i1:i2,j1:j2) 1037 END IF 1038 ENDIF 1039 item = Kmm_a 1040 1041 ELSE 1042 zrhox = Agrif_rhox() 1043 1044 IF( l_vremap .OR. l_ini_child ) THEN 1045 1046 IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 1047 942 1048 DO jj=j1,j2 943 1049 DO ji=i1,i2 944 jk = mbkv(ji,jj) 945 ptab(ji,jj,jk,2) = 0._wp 946 END DO 947 END DO 948 END IF 949 ! Save ssh at last level: 950 ptab(i1:i2,j1:j2,k2,2) = 0._wp 951 IF (.NOT.ln_linssh) THEN 952 ! This vertical sum below should be replaced by the sea-level at V-points (optimization): 953 DO jk=1,jpk 954 ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3v(i1:i2,j1:j2,jk,Kmm_a) * vmask(i1:i2,j1:j2,jk) 955 END DO 956 ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hv_0(i1:i2,j1:j2) 957 END IF 958 # endif 959 ELSE 960 zrhox = Agrif_rhox() 961 # if defined key_vertical 962 963 IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 964 965 DO jj=j1,j2 966 DO ji=i1,i2 967 vv(ji,jj,:,Krhs_a) = 0._wp 968 N_in = mbkv_parent(ji,jj) 969 zhtot = 0._wp 970 DO jk=1,N_in 971 IF (jk==N_in) THEN 972 h_in(jk) = hv0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot 973 ELSE 974 h_in(jk) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox) 1050 vv(ji,jj,:,Krhs_a) = 0._wp 1051 N_in = mbkv_parent(ji,jj) 1052 zhtot = 0._wp 1053 DO jk=1,N_in 1054 IF (jk==N_in) THEN 1055 h_in(jk) = hv0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot 1056 ELSE 1057 h_in(jk) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox) 1058 ENDIF 1059 zhtot = zhtot + h_in(jk) 1060 IF( h_in(jk) .GT. 0. ) THEN 1061 tabin(jk) = ptab(ji,jj,jk,1)/(e1v(ji,jj)*zrhox*h_in(jk)) 1062 ELSE 1063 tabin(jk) = 0. 1064 ENDIF 1065 END DO 1066 1067 z_in(1) = 0.5_wp * h_in(1) - zhtot + hv0_parent(ji,jj) 1068 DO jk=2,N_in 1069 z_in(jk) = z_in(jk-1) + 0.5_wp * h_in(jk) 1070 END DO 1071 1072 N_out = 0 1073 DO jk=1,jpk 1074 IF (vmask(ji,jj,jk) == 0) EXIT 1075 N_out = N_out + 1 1076 h_out(N_out) = e3v(ji,jj,jk,Krhs_a) 1077 END DO 1078 1079 z_out(1) = 0.5_wp * h_out(1) - SUM(h_out(1:N_out)) + hv_0(ji,jj) 1080 DO jk=2,N_out 1081 z_out(jk) = z_out(jk-1) + 0.5_wp * h_out(jk) 1082 END DO 1083 1084 IF (N_in*N_out > 0) THEN 1085 IF( l_ini_child ) THEN 1086 CALL remap_linear (tabin(1:N_in),z_in(1:N_in),vv(ji,jj,1:N_out,Krhs_a),z_out(1:N_out),N_in,N_out,1) 1087 ELSE 1088 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),vv(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out,1) 1089 ENDIF 975 1090 ENDIF 976 zhtot = zhtot + h_in(jk) 977 tabin(jk) = ptab(ji,jj,jk,1)/(e1v(ji,jj)*zrhox*h_in(jk)) 978 ENDDO 979 980 N_out = 0 981 DO jk=1,jpk 982 if (vmask(ji,jj,jk) == 0) EXIT 983 N_out = N_out + 1 984 h_out(N_out) = e3v(ji,jj,jk,Krhs_a) 985 END DO 986 IF (N_in*N_out > 0) THEN 987 call reconstructandremap(tabin(1:N_in),h_in(1:N_in),vv(ji,jj,1:N_out,Krhs_a),h_out(1:N_out),N_in,N_out,1) 988 ENDIF 989 END DO 990 END DO 991 # else 992 DO jk = 1, jpkm1 993 vv(i1:i2,j1:j2,jk,Krhs_a) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v(i1:i2,j1:j2,jk,Krhs_a) ) 994 END DO 995 # endif 1091 END DO 1092 END DO 1093 ELSE 1094 DO jk = 1, jpkm1 1095 vv(i1:i2,j1:j2,jk,Krhs_a) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v(i1:i2,j1:j2,jk,Krhs_a) ) 1096 END DO 1097 ENDIF 996 1098 ENDIF 997 1099 ! … … 1203 1305 END SUBROUTINE interpe3t 1204 1306 1205 1206 1307 SUBROUTINE interpglamt( ptab, i1, i2, j1, j2, before ) 1207 1308 !!---------------------------------------------------------------------- … … 1283 1384 END DO 1284 1385 END DO 1285 END DO 1286 1287 # if defined key_vertical 1288 ! Interpolate thicknesses 1289 ! Warning: these are masked, hence extrapolated prior interpolation. 1290 DO jk=k1,k2 1291 DO jj=j1,j2 1292 DO ji=i1,i2 1293 ptab(ji,jj,jk,2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 1294 END DO 1295 END DO 1296 END DO 1297 1298 ! Extrapolate thicknesses in partial bottom cells: 1299 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 1300 IF (ln_zps) THEN 1301 DO jj=j1,j2 1302 DO ji=i1,i2 1303 jk = mbkt(ji,jj) 1304 ptab(ji,jj,jk,2) = 0._wp 1305 END DO 1306 END DO 1307 END IF 1308 1309 ! Save ssh at last level: 1310 IF (.NOT.ln_linssh) THEN 1311 ptab(i1:i2,j1:j2,k2,2) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 1312 ELSE 1313 ptab(i1:i2,j1:j2,k2,2) = 0._wp 1314 END IF 1315 # endif 1386 END DO 1387 1388 IF( l_vremap ) THEN 1389 ! Interpolate thicknesses 1390 ! Warning: these are masked, hence extrapolated prior interpolation. 1391 DO jk=k1,k2 1392 DO jj=j1,j2 1393 DO ji=i1,i2 1394 ptab(ji,jj,jk,2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 1395 END DO 1396 END DO 1397 END DO 1398 1399 ! Extrapolate thicknesses in partial bottom cells: 1400 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 1401 IF (ln_zps) THEN 1402 DO jj=j1,j2 1403 DO ji=i1,i2 1404 jk = mbkt(ji,jj) 1405 ptab(ji,jj,jk,2) = 0._wp 1406 END DO 1407 END DO 1408 END IF 1409 1410 ! Save ssh at last level: 1411 IF (.NOT.ln_linssh) THEN 1412 ptab(i1:i2,j1:j2,k2,2) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 1413 ELSE 1414 ptab(i1:i2,j1:j2,k2,2) = 0._wp 1415 END IF 1416 ENDIF 1417 1316 1418 ELSE 1317 #ifdef key_vertical 1318 IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 1319 avm_k(i1:i2,j1:j2,k1:k2) = 0._wp 1320 1321 DO jj = j1, j2 1322 DO ji =i1, i2 1323 N_in = mbkt_parent(ji,jj) 1324 IF ( tmask(ji,jj,1) == 0._wp) N_in = 0 1325 z_in(N_in+1) = ht0_parent(ji,jj) + ptab(ji,jj,k2,2) 1326 DO jk = N_in, 1, -1 ! Parent vertical grid 1327 z_in(jk) = z_in(jk+1) - ptab(ji,jj,jk,2) 1328 tabin(jk) = ptab(ji,jj,jk,1) 1329 END DO 1330 N_out = mbkt(ji,jj) 1331 DO jk = 1, N_out ! Child vertical grid 1332 z_out(jk) = gdepw(ji,jj,jk,Kmm_a) 1333 ENDDO 1334 IF (N_in*N_out > 0) THEN 1335 CALL remap_linear(tabin(1:N_in),z_in(1:N_in),avm_k(ji,jj,1:N_out),z_out(1:N_out),N_in,N_out,1) 1336 ENDIF 1337 ENDDO 1338 ENDDO 1339 #else 1340 avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2,1) 1341 #endif 1419 1420 IF( l_vremap ) THEN 1421 IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 1422 avm_k(i1:i2,j1:j2,k1:k2) = 0._wp 1423 1424 DO jj = j1, j2 1425 DO ji =i1, i2 1426 N_in = mbkt_parent(ji,jj) 1427 IF ( tmask(ji,jj,1) == 0._wp) N_in = 0 1428 z_in(N_in+1) = ht0_parent(ji,jj) + ptab(ji,jj,k2,2) 1429 DO jk = N_in, 1, -1 ! Parent vertical grid 1430 z_in(jk) = z_in(jk+1) - ptab(ji,jj,jk,2) 1431 tabin(jk) = ptab(ji,jj,jk,1) 1432 END DO 1433 N_out = mbkt(ji,jj) 1434 DO jk = 1, N_out ! Child vertical grid 1435 z_out(jk) = gdepw(ji,jj,jk,Kmm_a) 1436 END DO 1437 IF (N_in*N_out > 0) THEN 1438 CALL remap_linear(tabin(1:N_in),z_in(1:N_in),avm_k(ji,jj,1:N_out),z_out(1:N_out),N_in,N_out,1) 1439 ENDIF 1440 END DO 1441 END DO 1442 ELSE 1443 avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2,1) 1444 ENDIF 1342 1445 ENDIF 1343 1446 ! 1344 1447 END SUBROUTINE interpavm 1345 1448 1346 # if defined key_vertical 1449 1347 1450 SUBROUTINE interpmbkt( ptab, i1, i2, j1, j2, before ) 1348 1451 !!---------------------------------------------------------------------- … … 1363 1466 END SUBROUTINE interpmbkt 1364 1467 1468 1365 1469 SUBROUTINE interpht0( ptab, i1, i2, j1, j2, before ) 1366 1470 !!---------------------------------------------------------------------- … … 1380 1484 ! 1381 1485 END SUBROUTINE interpht0 1382 #endif 1383 1486 1487 1488 SUBROUTINE agrif_initts(tabres,i1,i2,j1,j2,k1,k2,m1,m2,before) 1489 INTEGER :: i1, i2, j1, j2, k1, k2, m1, m2 1490 REAL(wp):: tabres(i1:i2,j1:j2,k1:k2,m1:m2) 1491 LOGICAL :: before 1492 1493 INTEGER :: jm 1494 1495 IF (before) THEN 1496 DO jm=1,jpts 1497 tabres(i1:i2,j1:j2,k1:k2,jm) = ts(i1:i2,j1:j2,k1:k2,jm,Kbb_a) 1498 END DO 1499 ELSE 1500 DO jm=1,jpts 1501 ts(i1:i2,j1:j2,k1:k2,jm,Kbb_a)=tabres(i1:i2,j1:j2,k1:k2,jm) 1502 END DO 1503 ENDIF 1504 END SUBROUTINE agrif_initts 1505 1506 1507 SUBROUTINE agrif_initssh( ptab, i1, i2, j1, j2, before ) 1508 !!---------------------------------------------------------------------- 1509 !! *** ROUTINE interpsshn *** 1510 !!---------------------------------------------------------------------- 1511 INTEGER , INTENT(in ) :: i1, i2, j1, j2 1512 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab 1513 LOGICAL , INTENT(in ) :: before 1514 ! 1515 !!---------------------------------------------------------------------- 1516 ! 1517 IF( before) THEN 1518 ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Kbb_a) 1519 ELSE 1520 ssh(i1:i2,j1:j2,Kbb_a) = ptab(i1:i2,j1:j2)*tmask(i1:i2,j1:j2,1) 1521 ENDIF 1522 ! 1523 END SUBROUTINE agrif_initssh 1524 1384 1525 #else 1385 1526 !!---------------------------------------------------------------------- -
NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/NST/agrif_oce_sponge.F90
r13229 r13230 32 32 33 33 !! * Substitutions 34 # include "do_loop_substitute.h90"35 !!----------------------------------------------------------------------36 !! NEMO/NST 4.0 , NEMO Consortium (2018)37 !! $Id$38 !! Software governed by the CeCILL license (see ./LICENSE)39 !!----------------------------------------------------------------------40 CONTAINS41 42 SUBROUTINE Agrif_Sponge_Tra43 !!----------------------------------------------------------------------44 !! *** ROUTINE Agrif_Sponge_Tra ***45 !!----------------------------------------------------------------------46 REAL(wp) :: zcoef ! local scalar47 48 !!----------------------------------------------------------------------49 !50 #if defined SPONGE51 !! Assume persistence:52 zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())53 54 CALL Agrif_Sponge55 Agrif_SpecialValue = 0._wp56 Agrif_UseSpecialValue = .TRUE.57 tabspongedone_tsn = .FALSE.58 !59 CALL Agrif_Bc_Variable( tsn_sponge_id, calledweight=zcoef, procname=interptsn_sponge )60 !61 Agrif_UseSpecialValue = .FALSE.62 #endif63 !64 CALL iom_put( 'agrif_spu', fspu(:,:))65 CALL iom_put( 'agrif_spv', fspv(:,:))66 !67 END SUBROUTINE Agrif_Sponge_Tra68 69 70 SUBROUTINE Agrif_Sponge_dyn71 !!----------------------------------------------------------------------72 !! *** ROUTINE Agrif_Sponge_dyn ***73 !!----------------------------------------------------------------------74 REAL(wp) :: zcoef ! local scalar75 !!----------------------------------------------------------------------76 !77 #if defined SPONGE78 zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())79 80 Agrif_SpecialValue=0.81 Agrif_UseSpecialValue = ln_spc_dyn82 !83 tabspongedone_u = .FALSE.84 tabspongedone_v = .FALSE.85 CALL Agrif_Bc_Variable( un_sponge_id, calledweight=zcoef, procname=interpun_sponge )86 !87 tabspongedone_u = .FALSE.88 tabspongedone_v = .FALSE.89 CALL Agrif_Bc_Variable( vn_sponge_id, calledweight=zcoef, procname=interpvn_sponge )90 !91 Agrif_UseSpecialValue = .FALSE.92 #endif93 !94 CALL iom_put( 'agrif_spt', fspt(:,:))95 CALL iom_put( 'agrif_spf', fspf(:,:))96 !97 END SUBROUTINE Agrif_Sponge_dyn98 99 100 SUBROUTINE Agrif_Sponge101 !!----------------------------------------------------------------------102 !! *** ROUTINE Agrif_Sponge ***103 !!----------------------------------------------------------------------104 INTEGER :: ji, jj, ind1, ind2105 INTEGER :: ispongearea, jspongearea106 REAL(wp) :: z1_ispongearea, z1_jspongearea107 REAL(wp), DIMENSION(jpi,jpj) :: ztabramp108 #if defined key_vertical109 REAL(wp), DIMENSION(jpi,jpj) :: ztabrampu110 REAL(wp), DIMENSION(jpi,jpj) :: ztabrampv111 #endif112 REAL(wp), DIMENSION(jpjmax) :: zmskwest, zmskeast113 REAL(wp), DIMENSION(jpimax) :: zmsknorth, zmsksouth114 !!----------------------------------------------------------------------115 !116 ! Sponge 1d example with:117 ! iraf = 3 ; nbghost = 3 ; nn_sponge_len = 2118 !119 !coarse : U T U T U T U120 !| | | | |121 !fine : t u t u t u t u t u t u t u t u t u t u t122 !sponge val:0 0 0 1 5/6 4/6 3/6 2/6 1/6 0 0123 ! | ghost | <-- sponge area -- > |124 ! | points | |125 ! |--> dynamical interface126 127 #if defined SPONGE || defined SPONGE_TOP128 IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN129 !130 ! Retrieve masks at open boundaries:131 132 ! --- West --- !133 IF( lk_west ) THEN134 ztabramp(:,:) = 0._wp135 ind1 = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells136 DO ji = mi0(ind1), mi1(ind1)137 ztabramp(ji,:) = ssumask(ji,:)138 END DO139 !140 zmskwest(1:jpj) = MAXVAL(ztabramp(:,:), dim=1)141 zmskwest(jpj+1:jpjmax) = 0._wp142 ENDIF143 144 ! --- East --- !145 IF( lk_east ) THEN146 ztabramp(:,:) = 0._wp147 ind1 = jpiglo - ( nn_hls + nbghostcells + 1) ! halo + land + nbghostcells148 DO ji = mi0(ind1), mi1(ind1)149 ztabramp(ji,:) = ssumask(ji,:)150 END DO151 !152 zmskeast(1:jpj) = MAXVAL(ztabramp(:,:), dim=1)153 zmskeast(jpj+1:jpjmax) = 0._wp154 ENDIF155 156 ! --- South --- !157 IF( lk_south ) THEN158 ztabramp(:,:) = 0._wp159 ind1 = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells160 DO jj = mj0(ind1), mj1(ind1)161 ztabramp(:,jj) = ssvmask(:,jj)162 END DO163 !164 zmsksouth(1:jpi) = MAXVAL(ztabramp(:,:), dim=2)165 zmsksouth(jpi+1:jpimax) = 0._wp166 ENDIF167 168 ! --- North --- !169 IF( lk_north ) THEN170 ztabramp(:,:) = 0._wp171 ind1 = jpjglo - ( nn_hls + nbghostcells + 1) ! halo + land + nbghostcells172 DO jj = mj0(ind1), mj1(ind1)173 ztabramp(:,jj) = ssvmask(:,jj)174 END DO175 !176 zmsknorth(1:jpi) = MAXVAL(ztabramp(:,:), dim=2)177 zmsknorth(jpi+1:jpimax) = 0._wp178 ENDIF179 180 ! JC: SPONGE MASKING TO BE SORTED OUT:181 zmskwest(:) = 1._wp182 zmskeast(:) = 1._wp183 zmsksouth(:) = 1._wp184 zmsknorth(:) = 1._wp185 #if defined key_mpp_mpi186 ! CALL mpp_max( 'AGRIF_sponge', zmskwest(:) , jpjmax )187 ! CALL mpp_max( 'AGRIF_Sponge', zmskeast(:) , jpjmax )188 ! CALL mpp_max( 'AGRIF_Sponge', zmsksouth(:), jpimax )189 ! CALL mpp_max( 'AGRIF_Sponge', zmsknorth(:), jpimax )190 #endif191 192 ! Define ramp from boundaries towards domain interior at T-points193 ! Store it in ztabramp194 195 ispongearea = nn_sponge_len * Agrif_irhox()196 z1_ispongearea = 1._wp / REAL( ispongearea )197 jspongearea = nn_sponge_len * Agrif_irhoy()198 z1_jspongearea = 1._wp / REAL( jspongearea )199 200 ztabramp(:,:) = 0._wp201 202 ! Trick to remove sponge in 2DV domains:203 IF ( nbcellsx <= 3 ) ispongearea = -1204 IF ( nbcellsy <= 3 ) jspongearea = -1205 206 ! --- West --- !207 IF( lk_west ) THEN208 ind1 = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells209 ind2 = nn_hls + 1 + nbghostcells + ispongearea210 DO ji = mi0(ind1), mi1(ind2)211 DO jj = 1, jpj212 ztabramp(ji,jj) = REAL( ind2 - mig(ji) ) * z1_ispongearea * zmskwest(jj)213 END DO214 END DO215 216 ! ghost cells:217 ind1 = 1218 ind2 = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells219 DO ji = mi0(ind1), mi1(ind2)220 DO jj = 1, jpj221 ztabramp(ji,jj) = zmskwest(jj)222 END DO223 END DO224 ENDIF225 226 ! --- East --- !227 IF( lk_east ) THEN228 ind1 = jpiglo - ( nn_hls + nbghostcells ) - ispongearea229 ind2 = jpiglo - ( nn_hls + nbghostcells ) ! halo + land + nbghostcells - 1230 DO ji = mi0(ind1), mi1(ind2)231 DO jj = 1, jpj232 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mig(ji) - ind1 ) * z1_ispongearea) * zmskeast(jj)233 ENDDO234 END DO235 236 ! ghost cells:237 ind1 = jpiglo - ( nn_hls + nbghostcells ) ! halo + land + nbghostcells - 1238 ind2 = jpiglo239 DO ji = mi0(ind1), mi1(ind2)240 DO jj = 1, jpj241 ztabramp(ji,jj) = zmskeast(jj)242 ENDDO243 END DO244 ENDIF245 246 ! --- South --- !247 IF( lk_south ) THEN248 ind1 = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells249 ind2 = nn_hls + 1 + nbghostcells + jspongearea250 DO jj = mj0(ind1), mj1(ind2)251 DO ji = 1, jpi252 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ind2 - mjg(jj) ) * z1_jspongearea) * zmsksouth(ji)253 END DO254 END DO255 256 ! ghost cells:257 ind1 = 1258 ind2 = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells259 DO jj = mj0(ind1), mj1(ind2)260 DO ji = 1, jpi261 ztabramp(ji,jj) = zmsksouth(ji)262 END DO263 END DO264 ENDIF265 266 ! --- North --- !267 IF( lk_north ) THEN268 ind1 = jpjglo - ( nn_hls + nbghostcells ) - jspongearea269 ind2 = jpjglo - ( nn_hls + nbghostcells ) ! halo + land + nbghostcells - 1270 DO jj = mj0(ind1), mj1(ind2)271 DO ji = 1, jpi272 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mjg(jj) - ind1 ) * z1_jspongearea) * zmsknorth(ji)273 END DO274 END DO275 276 ! ghost cells:277 ind1 = jpjglo - ( nn_hls + nbghostcells ) ! halo + land + nbghostcells - 1278 ind2 = jpjglo279 DO jj = mj0(ind1), mj1(ind2)280 DO ji = 1, jpi281 ztabramp(ji,jj) = zmsknorth(ji)282 END DO283 END DO284 ENDIF285 286 ENDIF287 288 ! Tracers289 IF( .NOT. spongedoneT ) THEN290 fspu(:,:) = 0._wp291 fspv(:,:) = 0._wp292 DO_2D_00_00293 fspu(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji+1,jj ) ) * ssumask(ji,jj)294 fspv(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji ,jj+1) ) * ssvmask(ji,jj)295 END_2D296 ENDIF297 298 ! Dynamics299 IF( .NOT. spongedoneU ) THEN300 fspt(:,:) = 0._wp301 fspf(:,:) = 0._wp302 DO_2D_00_00303 fspt(ji,jj) = ztabramp(ji,jj) * ssmask(ji,jj)304 fspf(ji,jj) = 0.25_wp * ( ztabramp(ji ,jj ) + ztabramp(ji ,jj+1) &305 & +ztabramp(ji+1,jj+1) + ztabramp(ji+1,jj ) ) &306 & * ssvmask(ji,jj) * ssvmask(ji,jj+1)307 END_2D308 ENDIF309 310 IF( .NOT. spongedoneT .AND. .NOT. spongedoneU ) THEN311 CALL lbc_lnk_multi( 'agrif_Sponge', fspu, 'U', 1., fspv, 'V', 1., fspt, 'T', 1., fspf, 'F', 1. )312 spongedoneT = .TRUE.313 spongedoneU = .TRUE.314 ENDIF315 IF( .NOT. spongedoneT ) THEN316 CALL lbc_lnk_multi( 'agrif_Sponge', fspu, 'U', 1., fspv, 'V', 1. )317 spongedoneT = .TRUE.318 ENDIF319 IF( .NOT. spongedoneT .AND. .NOT. spongedoneU ) THEN320 CALL lbc_lnk_multi( 'agrif_Sponge', fspt, 'T', 1., fspf, 'F', 1. )321 spongedoneU = .TRUE.322 ENDIF323 324 #if defined key_vertical325 ! Remove vertical interpolation where not needed:326 DO_2D_00_00327 IF ((fspu(ji-1,jj)==0._wp).AND.(fspu(ji,jj)==0._wp).AND. &328 & (fspv(ji,jj-1)==0._wp).AND.(fspv(ji,jj)==0._wp)) mbkt_parent(ji,jj) = 0329 !330 IF ((fspt(ji+1,jj)==0._wp).AND.(fspt(ji,jj)==0._wp).AND. &331 & (fspf(ji,jj-1)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbku_parent(ji,jj) = 0332 !333 IF ((fspt(ji,jj+1)==0._wp).AND.(fspt(ji,jj)==0._wp).AND. &334 & (fspf(ji-1,jj)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbkv_parent(ji,jj) = 0335 !336 IF ( ssmask(ji,jj) == 0._wp) mbkt_parent(ji,jj) = 0337 IF (ssumask(ji,jj) == 0._wp) mbku_parent(ji,jj) = 0338 IF (ssvmask(ji,jj) == 0._wp) mbkv_parent(ji,jj) = 0339 END_2D340 !341 ztabramp (:,:) = REAL( mbkt_parent (:,:), wp )342 ztabrampu(:,:) = REAL( mbku_parentu(:,:), wp )343 ztabrampv(:,:) = REAL( mbkv_parentv(:,:), wp )344 CALL lbc_lnk_multi( 'Agrif_Sponge', ztabramp, 'T', 1., ztabrampu, 'U', 1., ztabrampv, 'V', 1. )345 mbkt_parent(:,:) = NINT( ztabramp (:,:) )346 mbku_parent(:,:) = NINT( ztabrampu(:,:) )347 mbkv_parent(:,:) = NINT( ztabrampv(:,:) )348 #endif349 !350 #endif351 !352 END SUBROUTINE Agrif_Sponge353 354 SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )355 !!----------------------------------------------------------------------356 !! *** ROUTINE interptsn_sponge ***357 !!----------------------------------------------------------------------358 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2359 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres360 LOGICAL , INTENT(in ) :: before361 !362 INTEGER :: ji, jj, jk, jn ! dummy loop indices363 INTEGER :: iku, ikv364 REAL(wp) :: ztsa, zabe1, zabe2, zbtr, zhtot, ztrelax365 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv366 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff367 ! vertical interpolation:368 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child369 REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin370 REAL(wp), DIMENSION(k1:k2) :: h_in371 REAL(wp), DIMENSION(1:jpk) :: h_out372 INTEGER :: N_in, N_out373 !!----------------------------------------------------------------------374 !375 IF( before ) THEN376 DO jn = 1, jpts377 DO jk=k1,k2378 DO jj=j1,j2379 DO ji=i1,i2380 tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kbb_a)381 END DO382 END DO383 END DO384 END DO385 386 # if defined key_vertical387 ! Interpolate thicknesses388 ! Warning: these are masked, hence extrapolated prior interpolation.389 DO jk=k1,k2390 DO jj=j1,j2391 DO ji=i1,i2392 tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kbb_a)393 END DO394 END DO395 END DO396 397 ! Extrapolate thicknesses in partial bottom cells:398 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on399 IF (ln_zps) THEN400 DO jj=j1,j2401 DO ji=i1,i2402 jk = mbkt(ji,jj)403 tabres(ji,jj,jk,jpts+1) = 0._wp404 END DO405 END DO406 END IF407 408 ! Save ssh at last level:409 IF (.NOT.ln_linssh) THEN410 tabres(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kbb_a)*tmask(i1:i2,j1:j2,1)411 ELSE412 tabres(i1:i2,j1:j2,k2,jpts+1) = 0._wp413 END IF414 # endif415 416 ELSE417 !418 # if defined key_vertical419 420 IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp421 422 DO jj=j1,j2423 DO ji=i1,i2424 tabres_child(ji,jj,:,:) = 0._wp425 N_in = mbkt_parent(ji,jj)426 zhtot = 0._wp427 DO jk=1,N_in !k2 = jpk of parent grid428 IF (jk==N_in) THEN429 h_in(jk) = ht0_parent(ji,jj) + tabres(ji,jj,k2,n2) - zhtot430 ELSE431 h_in(jk) = tabres(ji,jj,jk,n2)432 ENDIF433 zhtot = zhtot + h_in(jk)434 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)435 END DO436 N_out = 0437 DO jk=1,jpk ! jpk of child grid438 IF (tmask(ji,jj,jk) == 0) EXIT439 N_out = N_out + 1440 h_out(jk) = e3t(ji,jj,jk,Kbb_a) !Child grid scale factors. Could multiply by e1e2t here instead of division above441 ENDDO442 443 ! Account for small differences in free-surface444 IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN445 h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )446 ELSE447 h_in(1) = h_in(1) - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) )448 ENDIF449 IF (N_in*N_out > 0) THEN450 CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts)451 ENDIF452 ENDDO453 ENDDO454 # endif455 456 DO jj=j1,j2457 DO ji=i1,i2458 DO jk=1,jpkm1459 # if defined key_vertical460 tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres_child(ji,jj,jk,1:jpts)) * tmask(ji,jj,jk)461 # else462 tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres(ji,jj,jk,1:jpts))*tmask(ji,jj,jk)463 # endif464 ENDDO465 ENDDO466 ENDDO467 468 !* set relaxation time scale469 IF( l_1st_euler .AND. lk_agrif_fstep ) THEN ; ztrelax = rn_trelax_tra / ( rn_Dt )470 ELSE ; ztrelax = rn_trelax_tra / (2._wp * rn_Dt )471 ENDIF472 473 DO jn = 1, jpts474 DO jk = 1, jpkm1475 ztu(i1:i2,j1:j2,jk) = 0._wp476 DO jj = j1,j2477 DO ji = i1,i2-1478 zabe1 = rn_sponge_tra * fspu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm_a)479 ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj ,jk,jn) - tsbdiff(ji,jj,jk,jn) )480 END DO481 END DO482 ztv(i1:i2,j1:j2,jk) = 0._wp483 DO ji = i1,i2484 DO jj = j1,j2-1485 zabe2 = rn_sponge_tra * fspv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a)486 ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )487 END DO488 END DO489 !490 IF( ln_zps ) THEN ! set gradient at partial step level491 DO jj = j1,j2492 DO ji = i1,i2493 ! last level494 iku = mbku(ji,jj)495 ikv = mbkv(ji,jj)496 IF( iku == jk ) ztu(ji,jj,jk) = 0._wp497 IF( ikv == jk ) ztv(ji,jj,jk) = 0._wp498 END DO499 END DO500 ENDIF501 END DO502 !503 DO jk = 1, jpkm1504 DO jj = j1+1,j2-1505 DO ji = i1+1,i2-1506 IF (.NOT. tabspongedone_tsn(ji,jj)) THEN507 zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a)508 ! horizontal diffusive trends509 ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) &510 & - ztrelax * fspt(ji,jj) * tsbdiff(ji,jj,jk,jn)511 ! add it to the general tracer trends512 ts(ji,jj,jk,jn,Krhs_a) = ts(ji,jj,jk,jn,Krhs_a) + ztsa513 ENDIF514 END DO515 END DO516 END DO517 !518 END DO519 !520 tabspongedone_tsn(i1+1:i2-1,j1+1:j2-1) = .TRUE.521 !522 ENDIF523 !524 END SUBROUTINE interptsn_sponge525 526 SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before)527 !!---------------------------------------------528 !! *** ROUTINE interpun_sponge ***529 !!---------------------------------------------530 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2531 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres532 LOGICAL, INTENT(in) :: before533 534 INTEGER :: ji,jj,jk,jmax535 INTEGER :: ind1536 ! sponge parameters537 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot, ztrelax538 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff539 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff540 ! vertical interpolation:541 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child542 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in543 REAL(wp), DIMENSION(1:jpk) :: h_out544 INTEGER ::N_in, N_out545 !!---------------------------------------------546 !547 IF( before ) THEN548 DO jk=k1,k2549 DO jj=j1,j2550 DO ji=i1,i2551 tabres(ji,jj,jk,m1) = uu(ji,jj,jk,Kbb_a)552 # if defined key_vertical553 tabres(ji,jj,jk,m2) = e3u(ji,jj,jk,Kbb_a)*umask(ji,jj,jk)554 # endif555 END DO556 END DO557 END DO558 559 # if defined key_vertical560 ! Extrapolate thicknesses in partial bottom cells:561 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on562 IF (ln_zps) THEN563 DO jj=j1,j2564 DO ji=i1,i2565 jk = mbku(ji,jj)566 tabres(ji,jj,jk,m2) = 0._wp567 END DO568 END DO569 END IF570 ! Save ssh at last level:571 tabres(i1:i2,j1:j2,k2,m2) = 0._wp572 IF (.NOT.ln_linssh) THEN573 ! This vertical sum below should be replaced by the sea-level at U-points (optimization):574 DO jk=1,jpk575 tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3u(i1:i2,j1:j2,jk,Kbb_a) * umask(i1:i2,j1:j2,jk)576 END DO577 tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hu_0(i1:i2,j1:j2)578 END IF579 # endif580 581 ELSE582 583 # if defined key_vertical584 IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp585 586 DO jj=j1,j2587 DO ji=i1,i2588 tabres_child(ji,jj,:) = 0._wp589 N_in = mbku_parent(ji,jj)590 zhtot = 0._wp591 DO jk=1,N_in592 IF (jk==N_in) THEN593 h_in(jk) = hu0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot594 ELSE595 h_in(jk) = tabres(ji,jj,jk,m2)596 ENDIF597 zhtot = zhtot + h_in(jk)598 tabin(jk) = tabres(ji,jj,jk,m1)599 ENDDO600 !601 N_out = 0602 DO jk=1,jpk603 IF (umask(ji,jj,jk) == 0) EXIT604 N_out = N_out + 1605 h_out(N_out) = e3u(ji,jj,jk,Kbb_a)606 ENDDO607 608 ! Account for small differences in free-surface609 IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN610 h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )611 ELSE612 h_in(1) = h_in(1) - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) )613 ENDIF614 615 IF (N_in * N_out > 0) THEN616 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1)617 ENDIF618 ENDDO619 ENDDO620 621 ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:)622 #else623 ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:)624 #endif625 !* set relaxation time scale626 IF( l_1st_euler .AND. lk_agrif_fstep ) THEN ; ztrelax = rn_trelax_dyn / ( rn_Dt )627 ELSE ; ztrelax = rn_trelax_dyn / (2._wp * rn_Dt )628 ENDIF629 !630 DO jk = 1, jpkm1 ! Horizontal slab631 ! ! ===============632 633 ! ! --------634 ! Horizontal divergence ! div635 ! ! --------636 DO jj = j1,j2637 DO ji = i1+1,i2 ! vector opt.638 zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb_a) * rn_sponge_dyn * fspt(ji,jj)639 hdivdiff(ji,jj,jk) = ( e2u(ji ,jj)*e3u(ji ,jj,jk,Kbb_a) * ubdiff(ji ,jj,jk) &640 & -e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb_a) * ubdiff(ji-1,jj,jk) ) * zbtr641 END DO642 END DO643 644 DO jj = j1,j2-1645 DO ji = i1,i2 ! vector opt.646 zbtr = r1_e1e2f(ji,jj) * e3f(ji,jj,jk) * rn_sponge_dyn * fspf(ji,jj)647 rotdiff(ji,jj,jk) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1,jk) &648 & +e1u(ji,jj ) * ubdiff(ji,jj ,jk) ) * fmask(ji,jj,jk) * zbtr649 END DO650 END DO651 END DO652 !653 DO jj = j1+1, j2-1654 DO ji = i1+1, i2-1 ! vector opt.655 656 IF (.NOT. tabspongedone_u(ji,jj)) THEN657 DO jk = 1, jpkm1 ! Horizontal slab658 ze2u = rotdiff (ji,jj,jk)659 ze1v = hdivdiff(ji,jj,jk)660 ! horizontal diffusive trends661 zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) ) &662 & + ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj) &663 & - ztrelax * fspu(ji,jj) * ubdiff(ji,jj,jk)664 665 ! add it to the general momentum trends666 uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a) + zua667 END DO668 ENDIF669 670 END DO671 END DO672 673 tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE.674 675 jmax = j2-1676 ind1 = jpjglo - ( nn_hls + nbghostcells + 2 ) ! North677 DO jj = mj0(ind1), mj1(ind1)678 jmax = MIN(jmax,jj)679 END DO680 681 DO jj = j1+1, jmax682 DO ji = i1+1, i2 ! vector opt.683 684 IF (.NOT. tabspongedone_v(ji,jj)) THEN685 DO jk = 1, jpkm1 ! Horizontal slab686 ze2u = rotdiff (ji,jj,jk)687 ze1v = hdivdiff(ji,jj,jk)688 689 ! horizontal diffusive trends690 zva = + ( ze2u - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) ) &691 + ( hdivdiff(ji,jj+1,jk) - ze1v ) * r1_e2v(ji,jj)692 693 ! add it to the general momentum trends694 vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a) + zva695 END DO696 ENDIF697 !698 END DO699 END DO700 !701 tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.702 !703 ENDIF704 !705 END SUBROUTINE interpun_sponge706 707 SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before,nb,ndir)708 !!---------------------------------------------709 !! *** ROUTINE interpvn_sponge ***710 !!---------------------------------------------711 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2712 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres713 LOGICAL, INTENT(in) :: before714 INTEGER, INTENT(in) :: nb , ndir715 !716 INTEGER :: ji, jj, jk, imax717 INTEGER :: ind1718 ! sponge parameters719 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot, ztrelax720 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff721 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff722 ! vertical interpolation:723 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child724 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in725 REAL(wp), DIMENSION(1:jpk) :: h_out726 INTEGER :: N_in, N_out727 !!---------------------------------------------728 729 IF( before ) THEN730 DO jk=k1,k2731 DO jj=j1,j2732 34 DO ji=i1,i2 733 35 tabres(ji,jj,jk,m1) = vv(ji,jj,jk,Kbb_a) … … 778 80 zhtot = zhtot + h_in(jk) 779 81 tabin(jk) = tabres(ji,jj,jk,m1) 780 END DO82 END DO 781 83 ! 782 84 N_out = 0 … … 785 87 N_out = N_out + 1 786 88 h_out(N_out) = e3v(ji,jj,jk,Kbb_a) 787 END DO89 END DO 788 90 789 91 ! Account for small differences in free-surface … … 797 99 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 798 100 ENDIF 799 END DO800 END DO101 END DO 102 END DO 801 103 802 104 vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:) … … 804 106 vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:) 805 107 # endif 806 !* set relaxation time scale807 IF( l_1st_euler .AND. lk_agrif_fstep ) THEN ; ztrelax = rn_trelax_dyn / ( rn_Dt )808 ELSE ; ztrelax = rn_trelax_dyn / (2._wp * rn_Dt )809 ENDIF810 108 ! 811 109 DO jk = 1, jpkm1 ! Horizontal slab … … 817 115 DO jj = j1+1,j2 818 116 DO ji = i1,i2 ! vector opt. 819 zbtr = r 1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb_a) * rn_sponge_dyn * fspt(ji,jj)117 zbtr = rn_sponge_dyn * r1_Dt * fspt(ji,jj) / e3t(ji,jj,jk,Kbb_a) 820 118 hdivdiff(ji,jj,jk) = ( e1v(ji,jj ) * e3v(ji,jj ,jk,Kbb_a) * vbdiff(ji,jj ,jk) & 821 119 & -e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kbb_a) * vbdiff(ji,jj-1,jk) ) * zbtr … … 824 122 DO jj = j1,j2 825 123 DO ji = i1,i2-1 ! vector opt. 826 zbtr = r 1_e1e2f(ji,jj) * e3f(ji,jj,jk) * rn_sponge_dyn * fspf(ji,jj)124 zbtr = rn_sponge_dyn * r1_Dt * fspf(ji,jj) * e3f(ji,jj,jk) 827 125 rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) & 828 126 & -e2v(ji ,jj) * vbdiff(ji ,jj,jk) ) * fmask(ji,jj,jk) * zbtr … … 844 142 IF( .NOT. tabspongedone_u(ji,jj) ) THEN 845 143 DO jk = 1, jpkm1 846 uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a) 144 uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a) & 847 145 & - ( rotdiff (ji ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) ) & 848 146 & + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj ,jk)) * r1_e1u(ji,jj) … … 858 156 IF( .NOT. tabspongedone_v(ji,jj) ) THEN 859 157 DO jk = 1, jpkm1 860 vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a) 158 vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a) & 861 159 & + ( rotdiff (ji,jj ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) ) & 862 & + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji ,jj,jk) ) * r1_e2v(ji,jj) &863 & - ztrelax* fspv(ji,jj) * vbdiff(ji,jj,jk)160 & + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji ,jj,jk) ) * r1_e2v(ji,jj) & 161 & - rn_trelax_dyn * r1_Dt * fspv(ji,jj) * vbdiff(ji,jj,jk) 864 162 END DO 865 163 ENDIF -
NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/DYN/sshwzv.F90
r13229 r13230 200 200 ENDIF 201 201 ! 202 #if defined key_agrif203 202 IF( .NOT. AGRIF_Root() ) THEN 204 203 ! … … 207 206 DO jk = 1, jpkm1 208 207 ! --- West --- ! 209 IF( lk_west 208 IF( lk_west) THEN 210 209 DO ji = mi0(2+nn_hls), mi1(2+nn_hls) 211 210 DO jj = 1, jpj … … 216 215 ! 217 216 ! --- East --- ! 218 IF( lk_east 217 IF( lk_east) THEN 219 218 DO ji = mi0(jpiglo-1-nn_hls), mi1(jpiglo-1-nn_hls) 220 219 DO jj = 1, jpj … … 225 224 ! 226 225 ! --- South --- ! 227 IF( lk_south 226 IF( lk_south) THEN 228 227 DO jj = mj0(2+nn_hls), mj1(2+nn_hls) 229 228 DO ji = 1, jpi … … 234 233 ! 235 234 ! --- North --- ! 236 IF( lk_north 235 IF( lk_north) THEN 237 236 DO jj = mj0(jpjglo-1-nn_hls), mj1(jpjglo-1-nn_hls) 238 237 DO ji = 1, jpi … … 241 240 END DO 242 241 ENDIF 242 ! 243 243 END DO 244 244 ! 245 245 ENDIF 246 #endif247 246 ! 248 247 IF( ln_timing ) CALL timing_stop('wzv')
Note: See TracChangeset
for help on using the changeset viewer.
