Changeset 11191 for NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyini.F90
- Timestamp:
- 2019-06-27T10:14:39+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyini.F90
r11183 r11191 33 33 PRIVATE 34 34 35 PUBLIC bdy_init ! routine called in nemo_init35 PUBLIC bdy_init ! routine called in nemo_init 36 36 PUBLIC find_neib ! routine called in bdy_nmn 37 37 … … 126 126 !! ** Input : bdy_init.nc, input file for unstructured open boundaries 127 127 !!---------------------------------------------------------------------- 128 INTEGER :: ib_bdy, ii, ij, ik, igrd, ib, ir, iseg ! dummy loop indices 129 INTEGER :: icount, icountr, ibr_max, ilen1, ibm1 ! local integers 128 INTEGER :: ib_bdy, ii, ij, igrd, ib, ir, iseg ! dummy loop indices 129 INTEGER :: icount, icountr, icountr0, ibr_max ! local integers 130 INTEGER :: ilen1, ibm1 ! - - 130 131 INTEGER :: iwe, ies, iso, ino, inum, id_dummy ! - - 131 INTEGER :: igrd_start, igrd_end, jpbdta ! - - 132 INTEGER :: jpbdtau, jpbdtas ! - - 132 INTEGER :: jpbdta, jpbdtau, jpbdtas ! - - 133 133 INTEGER :: ib_bdy1, ib_bdy2, ib1, ib2 ! - - 134 INTEGER :: i_offset, j_offset, inn ! - -135 134 INTEGER :: ii1, ii2, ii3, ij1, ij2, ij3 ! - - 136 135 INTEGER :: iibe, ijbe, iibi, ijbi ! - - 137 136 INTEGER :: flagu, flagv ! short cuts 138 INTEGER , POINTER :: nbi, nbj, nbr ! short cuts139 137 REAL(wp), POINTER, DIMENSION(:,:) :: pmask ! pointer to 2D mask fields 140 138 REAL(wp) :: zefl, zwfl, znfl, zsfl ! local scalars … … 145 143 CHARACTER(LEN=1),DIMENSION(jpbgrd) :: cgrid 146 144 INTEGER :: jpk_max 147 LOGICAL :: llsend_ea, llsend_we, llsend_so, llsend_no ! Flags for boundaries sending 148 LOGICAL :: llrecv_ea, llrecv_we, llrecv_so, llrecv_no ! Flags for boundaries receiving 149 REAL(wp), TARGET, DIMENSION(jpi,jpj) :: zfmask ! temporary fmask array excluding coastal boundary condition (shlat) 150 REAL(wp) , DIMENSION(jpi,jpj) :: ztmp 151 LOGICAL :: llnon, llson, llean, llwen ! local logicals 145 REAL(wp), DIMENSION(jpi,jpj) :: zfmask ! temporary fmask array excluding coastal boundary condition (shlat) 146 REAL(wp), DIMENSION(jpi,jpj) :: ztmask, zumask, zvmask ! temporary u/v mask array 152 147 !! 153 148 CHARACTER(LEN=1) :: ctypebdy ! - - … … 798 793 ! 799 794 iwe = mig(1) 800 ies = mig( nlci)795 ies = mig(jpi) 801 796 iso = mjg(1) 802 ino = mjg( nlcj)797 ino = mjg(jpj) 803 798 ! 804 799 DO ib_bdy = 1, nb_bdy 805 800 DO igrd = 1, jpbgrd 806 icount = 0 ! initialization of local bdy length 807 icountr = 0 ! initialization of local rim bdy length 808 idx_bdy(ib_bdy)%nblen(igrd) = 0 809 idx_bdy(ib_bdy)%nblenrim(igrd) = 0 801 icount = 0 ! initialization of local bdy length 802 icountr = 0 ! initialization of local rim 0 and rim 1 bdy length 803 icountr0 = 0 ! initialization of local rim 0 bdy length 804 idx_bdy(ib_bdy)%nblen(igrd) = 0 805 idx_bdy(ib_bdy)%nblenrim(igrd) = 0 806 idx_bdy(ib_bdy)%nblenrim0(igrd) = 0 810 807 DO ib = 1, nblendta(igrd,ib_bdy) 811 808 ! check that data is in correct order in file … … 823 820 ! 824 821 icount = icount + 1 825 IF( nbrdta(ib,igrd,ib_bdy) == 1 ) icountr = icountr + 1 822 IF( nbrdta(ib,igrd,ib_bdy) == 1 .OR. nbrdta(ib,igrd,ib_bdy) == 0 ) icountr = icountr + 1 823 IF( nbrdta(ib,igrd,ib_bdy) == 0 ) icountr0 = icountr0 + 1 826 824 ENDIF 827 825 END DO 828 idx_bdy(ib_bdy)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc 829 idx_bdy(ib_bdy)%nblen (igrd) = icount !: length of boundary data on each proc 826 idx_bdy(ib_bdy)%nblen (igrd) = icount !: length of boundary data on each proc 827 idx_bdy(ib_bdy)%nblenrim (igrd) = icountr !: length of rim 0 and rim 1 boundary data on each proc 828 idx_bdy(ib_bdy)%nblenrim0(igrd) = icountr0 !: length of rim 0 boundary data on each proc 830 829 END DO ! igrd 831 830 … … 849 848 icount = 0 850 849 ! Outer loop on rimwidth to ensure outermost points come first in the local arrays. 851 DO ir =1, nn_rimwidth(ib_bdy)850 DO ir = 0, nn_rimwidth(ib_bdy) 852 851 DO ib = 1, nblendta(igrd,ib_bdy) 853 852 ! check if point is in local domain and equals ir … … 870 869 ! Initialize array indicating communications in bdy 871 870 ! ------------------------------------------------- 872 ALLOCATE( lsend_bdy(nb_bdy,jpbgrd,4 ), lrecv_bdy(nb_bdy,jpbgrd,4) )873 lsend_bdy(:,:,: ) = .false.874 lrecv_bdy(:,:,: ) = .false.871 ALLOCATE( lsend_bdy(nb_bdy,jpbgrd,4,0:1), lrecv_bdy(nb_bdy,jpbgrd,4,0:1) ) 872 lsend_bdy(:,:,:,:) = .false. 873 lrecv_bdy(:,:,:,:) = .false. 875 874 876 875 DO ib_bdy = 1, nb_bdy … … 879 878 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 880 879 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 880 IF( ib .LE. idx_bdy(ib_bdy)%nblenrim0(igrd) ) THEN ; ir = 0 881 ELSE ; ir = 1 882 END IF 881 883 ! 882 884 ! check if point has to be sent to a neighbour 883 885 ! W neighbour and on the inner left side 884 IF( ii == 2 .and. (nbondi == 0 .or. nbondi == 1) ) lsend_bdy(ib_bdy,igrd,1 ) = .true.886 IF( ii == 2 .and. (nbondi == 0 .or. nbondi == 1) ) lsend_bdy(ib_bdy,igrd,1,ir) = .true. 885 887 ! E neighbour and on the inner right side 886 IF( ii == jpi-1 .and. (nbondi == 0 .or. nbondi == -1) ) lsend_bdy(ib_bdy,igrd,2 ) = .true.888 IF( ii == jpi-1 .and. (nbondi == 0 .or. nbondi == -1) ) lsend_bdy(ib_bdy,igrd,2,ir) = .true. 887 889 ! S neighbour and on the inner down side 888 IF( ij == 2 .and. (nbondj == 0 .or. nbondj == 1) ) lsend_bdy(ib_bdy,igrd,3 ) = .true.890 IF( ij == 2 .and. (nbondj == 0 .or. nbondj == 1) ) lsend_bdy(ib_bdy,igrd,3,ir) = .true. 889 891 ! N neighbour and on the inner up side 890 IF( ij == jpj-1 .and. (nbondj == 0 .or. nbondj == -1) ) lsend_bdy(ib_bdy,igrd,4 ) = .true.892 IF( ij == jpj-1 .and. (nbondj == 0 .or. nbondj == -1) ) lsend_bdy(ib_bdy,igrd,4,ir) = .true. 891 893 ! 892 894 ! check if point has to be received from a neighbour 893 895 ! W neighbour and on the outter left side 894 IF( ii == 1 .and. (nbondi == 0 .or. nbondi == 1) ) lrecv_bdy(ib_bdy,igrd,1) = .true.896 IF( ii == 1 .and. (nbondi == 0 .or. nbondi == 1) ) lrecv_bdy(ib_bdy,igrd,1,ir) = .true. 895 897 ! E neighbour and on the outter right side 896 IF( ii == jpi .and. (nbondi == 0 .or. nbondi == -1) ) lrecv_bdy(ib_bdy,igrd,2) = .true.898 IF( ii == jpi .and. (nbondi == 0 .or. nbondi == -1) ) lrecv_bdy(ib_bdy,igrd,2,ir) = .true. 897 899 ! S neighbour and on the outter down side 898 IF( ij == 1 .and. (nbondj == 0 .or. nbondj == 1) ) lrecv_bdy(ib_bdy,igrd,3) = .true.900 IF( ij == 1 .and. (nbondj == 0 .or. nbondj == 1) ) lrecv_bdy(ib_bdy,igrd,3,ir) = .true. 899 901 ! N neighbour and on the outter up side 900 IF( ij == jpj .and. (nbondj == 0 .or. nbondj == -1) ) lrecv_bdy(ib_bdy,igrd,4) = .true.902 IF( ij == jpj .and. (nbondj == 0 .or. nbondj == -1) ) lrecv_bdy(ib_bdy,igrd,4,ir) = .true. 901 903 ! 902 904 END DO … … 907 909 DO igrd = 1, jpbgrd 908 910 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 909 nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)910 idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( REAL( nbr - 1 ) *0.5 ) ! tanh formulation911 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = (REAL(nn_rimwidth(ib_bdy)+1- nbr)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic912 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = REAL(nn_rimwidth(ib_bdy)+1- nbr)/REAL(nn_rimwidth(ib_bdy)) ! linear911 ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) ) ! both rim 0 and rim 1 have the same weights 912 idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( REAL( ir - 1 ) *0.5 ) ! tanh formulation 913 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = (REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic 914 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)) ! linear 913 915 END DO 914 916 END DO … … 918 920 DO igrd = 1, jpbgrd 919 921 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 920 nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)922 ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) ) ! both rim 0 and rim 1 have the same damping coefficients 921 923 idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) & 922 & *(REAL(nn_rimwidth(ib_bdy)+1- nbr)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic924 & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic 923 925 idx_bdy(ib_bdy)%nbdout(ib,igrd) = 1. / ( rn_time_dmp_out(ib_bdy) * rday ) & 924 & *(REAL(nn_rimwidth(ib_bdy)+1- nbr)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic926 & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2. ! quadratic 925 927 END DO 926 928 END DO … … 931 933 ! Initialise masks and find normal/tangential directions 932 934 ! ------------------------------------------------------ 935 936 ! ------------------------------------------ 937 ! handel rim0, do as if rim 1 was free ocean 938 ! ------------------------------------------ 939 940 ztmask(:,:) = tmask(:,:,1) ; zumask(:,:) = umask(:,:,1) ; zvmask(:,:) = vmask(:,:,1) 941 ! For the flagu/flagv calculation below we require a version of fmask without 942 ! the land boundary condition (shlat) included: 943 DO ij = 2, jpjm1 944 DO ii = 2, jpim1 945 zfmask(ii,ij) = ztmask(ii,ij ) * ztmask(ii+1,ij ) & 946 & * ztmask(ii,ij+1) * ztmask(ii+1,ij+1) 947 END DO 948 END DO 949 CALL lbc_lnk( 'bdyini', zfmask, 'F', 1. ) 933 950 934 951 ! Read global 2D mask at T-points: bdytmask … … 940 957 941 958 ! Derive mask on U and V grid from mask on T grid 942 943 bdyumask(:,:) = 0._wp944 bdyvmask(:,:) = 0._wp945 959 DO ij = 1, jpjm1 946 960 DO ii = 1, jpim1 947 bdyumask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii+1, ij)961 bdyumask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii+1,ij ) 948 962 bdyvmask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii ,ij+1) 949 963 END DO 950 964 END DO 951 CALL lbc_lnk_multi( 'bdyini', bdyumask, 'U', 1. , bdyvmask, 'V', 1. ) ! Lateral boundary cond. 952 953 ! bdy masks are now set to zero on boundary points: 954 ! 955 igrd = 1 965 CALL lbc_lnk_multi( 'bdyini', bdyumask, 'U', 1., bdyvmask, 'V', 1. ) ! Lateral boundary cond. 966 967 ! bdy masks are now set to zero on rim 0 points: 956 968 DO ib_bdy = 1, nb_bdy 957 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 958 bdytmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0._wp 959 END DO 969 DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(1) ! extent of rim 0 970 bdytmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp 971 END DO 972 DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(2) ! extent of rim 0 973 bdyumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp 974 END DO 975 DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(3) ! extent of rim 0 976 bdyvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp 977 END DO 960 978 END DO 961 ! 962 igrd = 2 979 980 CALL bdy_rim_treat( zumask, zvmask, zfmask, .true. ) ! compute flagu, flagv, ntreat on rim 0 981 982 ! ------------------------------------ 983 ! handel rim1, do as if rim 0 was land 984 ! ------------------------------------ 985 963 986 DO ib_bdy = 1, nb_bdy 964 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 965 bdyumask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0._wp 966 END DO 987 DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(1) ! extent of rim 0 988 ztmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp 989 END DO 990 DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(2) ! extent of rim 0 991 zumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp 992 END DO 993 DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(3) ! extent of rim 0 994 zvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp 995 END DO 967 996 END DO 968 ! 969 igrd = 3 970 DO ib_bdy = 1, nb_bdy 971 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 972 bdyvmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0._wp 973 END DO 974 END DO 975 976 ! For the flagu/flagv calculation below we require a version of fmask without 977 ! the land boundary condition (shlat) included: 978 zfmask(:,:) = 0 997 998 ! Recompute zfmask 979 999 DO ij = 2, jpjm1 980 1000 DO ii = 2, jpim1 981 zfmask(ii,ij) = tmask(ii,ij ,1) * tmask(ii+1,ij ,1) &982 & * tmask(ii,ij+1,1) * tmask(ii+1,ij+1,1)1001 zfmask(ii,ij) = ztmask(ii,ij ) * ztmask(ii+1,ij ) & 1002 & * ztmask(ii,ij+1) * ztmask(ii+1,ij+1) 983 1003 END DO 984 1004 END DO 985 986 ! Lateral boundary conditions 987 CALL lbc_lnk( 'bdyini', zfmask, 'F', 1. ) 988 CALL lbc_lnk_multi( 'bdyini', bdyumask, 'U', 1. , bdyvmask, 'V', 1., bdytmask, 'T', 1. ) 1005 CALL lbc_lnk( 'bdyini', zfmask, 'F', 1. ) 1006 1007 ! bdy masks are now set to zero on rim1 points: 1008 DO ib_bdy = 1, nb_bdy 1009 DO ib = idx_bdy(ib_bdy)%nblenrim0(1) + 1, idx_bdy(ib_bdy)%nblenrim(1) ! extent of rim 1 1010 bdytmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp 1011 END DO 1012 DO ib = idx_bdy(ib_bdy)%nblenrim0(2) + 1, idx_bdy(ib_bdy)%nblenrim(2) ! extent of rim 1 1013 bdyumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp 1014 END DO 1015 DO ib = idx_bdy(ib_bdy)%nblenrim0(3) + 1, idx_bdy(ib_bdy)%nblenrim(3) ! extent of rim 1 1016 bdyvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp 1017 END DO 1018 END DO 1019 1020 CALL bdy_rim_treat( zumask, zvmask, zfmask, .false. ) ! compute flagu, flagv, ntreat on rim 1 1021 1022 1023 ! 1024 ! Check which boundaries might need communication 1025 ALLOCATE( lsend_bdyint(nb_bdy,jpbgrd,4,0:1), lrecv_bdyint(nb_bdy,jpbgrd,4,0:1) ) 1026 lsend_bdyint(:,:,:,:) = .false. 1027 lrecv_bdyint(:,:,:,:) = .false. 1028 ALLOCATE( lsend_bdyext(nb_bdy,jpbgrd,4,0:1), lrecv_bdyext(nb_bdy,jpbgrd,4,0:1) ) 1029 lsend_bdyext(:,:,:,:) = .false. 1030 lrecv_bdyext(:,:,:,:) = .false. 1031 ! 1032 DO igrd = 1, jpbgrd 1033 DO ib_bdy = 1, nb_bdy 1034 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 1035 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 1036 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 1037 IF( ib .LE. idx_bdy(ib_bdy)%nblenrim0(igrd) ) THEN ; ir = 0 1038 ELSE ; ir = 1 1039 END IF 1040 flagu = NINT(idx_bdy(ib_bdy)%flagu(ib,igrd)) 1041 flagv = NINT(idx_bdy(ib_bdy)%flagv(ib,igrd)) 1042 iibe = ii - flagu ! neighbouring point towards the exterior of the computational domain 1043 ijbe = ij - flagv 1044 iibi = ii + flagu ! neighbouring point towards the interior of the computational domain 1045 ijbi = ij + flagv 1046 CALL find_neib( ii, ij, idx_bdy(ib_bdy)%ntreat(ib,igrd), ii1, ij1, ii2, ij2, ii3, ij3 ) ! free ocean neighbours 1047 ! 1048 ! search neighbour in the west/east direction 1049 ! Rim is on the halo and computed ocean is towards exterior of mpi domain 1050 ! <-- (o exterior) --> 1051 ! (1) o|x OR (2) x|o 1052 ! |___ ___| 1053 IF( iibi == 0 .OR. ii1 == 0 .OR. ii2 == 0 .OR. ii3 == 0 ) lrecv_bdyint(ib_bdy,igrd,1,ir) = .true. 1054 IF( iibi == jpi+1 .OR. ii1 == jpi+1 .OR. ii2 == jpi+1 .OR. ii3 == jpi+1 ) lrecv_bdyint(ib_bdy,igrd,2,ir) = .true. 1055 IF( iibe == 0 ) lrecv_bdyext(ib_bdy,igrd,1,ir) = .true. 1056 IF( iibe == jpi+1 ) lrecv_bdyext(ib_bdy,igrd,2,ir) = .true. 1057 ! Check if neighbour has its rim parallel to its mpi subdomain border and located next to its halo 1058 ! :¨¨¨¨¨|¨¨--> | | <--¨¨|¨¨¨¨¨: 1059 ! : | x:o | neighbour limited by ... would need o | o:x | : 1060 ! :.....|_._:_____| (1) W neighbour E neighbour (2) |_____:_._|.....: 1061 IF( ii == 2 .AND. ( nbondi == 1 .OR. nbondi == 0 ) .AND. & 1062 & ( iibi == 3 .OR. ii1 == 3 .OR. ii2 == 3 .OR. ii3 == 3 ) ) lsend_bdyint(ib_bdy,igrd,1,ir) = .true. 1063 IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. & 1064 & ( iibi == jpi-2 .OR. ii1 == jpi-2 .OR. ii2 == jpi-2 .OR. ii3 == jpi-2) ) lsend_bdyint(ib_bdy,igrd,2,ir) = .true. 1065 IF( ii == 2 .AND. ( nbondi == 1 .OR. nbondi == 0 ) .AND. iibe == 3 ) lsend_bdyext(ib_bdy,igrd,1,ir) = .true. 1066 IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. iibe == jpi-2) lsend_bdyext(ib_bdy,igrd,2,ir) = .true. 1067 ! 1068 ! search neighbour in the north/south direction 1069 ! Rim is on the halo and computed ocean is towards exterior of mpi domain 1070 !(3) | | ^ ___o___ 1071 ! | |___x___| OR | | x | 1072 ! v o (4) | | 1073 IF( ijbi == 0 .OR. ij1 == 0 .OR. ij2 == 0 .OR. ij3 == 0 ) lrecv_bdyint(ib_bdy,igrd,3,ir) = .true. 1074 IF( ijbi == jpj+1 .OR. ij1 == jpj+1 .OR. ij2 == jpj+1 .OR. ij3 == jpj+1 ) lrecv_bdyint(ib_bdy,igrd,4,ir) = .true. 1075 IF( ijbe == 0 ) lrecv_bdyext(ib_bdy,igrd,3,ir) = .true. 1076 IF( ijbe == jpj+1 ) lrecv_bdyext(ib_bdy,igrd,4,ir) = .true. 1077 ! Check if neighbour has its rim parallel to its mpi subdomain _________ border and next to its halo 1078 ! ^ | o | : : 1079 ! | |¨¨¨¨x¨¨¨¨| neighbour limited by ... would need o | |....x....| 1080 ! :_________: (3) S neighbour N neighbour (4) v | o | 1081 IF( ij == 2 .AND. ( nbondj == 1 .OR. nbondj == 0 ) .AND. & 1082 & ( ijbi == 3 .OR. ij1 == 3 .OR. ij2 == 3 .OR. ij3 == 3 ) ) lsend_bdyint(ib_bdy,igrd,3,ir) = .true. 1083 IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. & 1084 & ( ijbi == jpj-2 .OR. ij1 == jpj-2 .OR. ij2 == jpj-2 .OR. ij3 == jpj-2) ) lsend_bdyint(ib_bdy,igrd,4,ir) = .true. 1085 IF( ij == 2 .AND. ( nbondj == 1 .OR. nbondj == 0 ) .AND. ijbe == 3 ) lsend_bdyext(ib_bdy,igrd,3,ir) = .true. 1086 IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. ijbe == jpj-2) lsend_bdyext(ib_bdy,igrd,4,ir) = .true. 1087 END DO 1088 END DO 1089 END DO 1090 1091 DO ib_bdy = 1,nb_bdy 1092 IF( cn_dyn2d(ib_bdy) == 'orlanski' .OR. cn_dyn2d(ib_bdy) == 'orlanski_npo' .OR. & 1093 & cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' .OR. & 1094 & cn_tra(ib_bdy) == 'orlanski' .OR. cn_tra(ib_bdy) == 'orlanski_npo' ) THEN 1095 DO igrd = 1, jpbgrd 1096 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 1097 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 1098 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 1099 IF( mig(ii) > 2 .AND. mig(ii) < jpiglo-2 .AND. mjg(ij) > 2 .AND. mjg(ij) < jpjglo-2 ) THEN 1100 WRITE(ctmp1,*) ' Orlanski is not safe when the open boundaries are on the interior of the computational domain' 1101 WRITE(ctmp2,*) ' ========== ' 1102 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 1103 END IF 1104 END DO 1105 END DO 1106 END IF 1107 END DO 1108 1109 ! 1110 ! Tidy up 1111 !-------- 1112 IF( nb_bdy>0 ) DEALLOCATE( nbidta, nbjdta, nbrdta ) 1113 ! 1114 END SUBROUTINE bdy_segs 1115 1116 1117 SUBROUTINE bdy_rim_treat( zumask, zvmask, zfmask, lrim0 ) 1118 !!---------------------------------------------------------------------- 1119 !! *** ROUTINE bdy_rim_treat *** 1120 !! 1121 !! ** Purpose : Initialize structures ( flagu, flagv, ntreat ) indicating how rim points 1122 !! are to be handeled in the boundary condition treatment 1123 !! 1124 !! ** Method : - to handel rim 0 zmasks must indicate ocean points (set at one on rim 0 and rim 1 and interior) 1125 !! and bdymasks must be set at 0 on rim 0 (set at one on rim 1 and interior) 1126 !! (as if rim 1 was free ocean) 1127 !! - to handel rim 1 zmasks must be set at 0 on rim 0 (set at one on rim 1 and interior) 1128 !! and bdymasks must indicate free ocean points (set at one on interior) 1129 !! (as if rim 0 was land) 1130 !! - we can then check in which direction the interior of the computational domain is with the difference 1131 !! mask array values on both sides to compute flagu and flagv 1132 !! - and look at the ocean neighbours to compute ntreat 1133 !!---------------------------------------------------------------------- 1134 REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in ) :: zfmask ! temporary fmask excluding coastal boundary condition (shlat) 1135 REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in ) :: zumask, zvmask ! temporary t/u/v mask array 1136 LOGICAL , INTENT (in ) :: lrim0 ! .true. -> rim 0 .false. -> rim 1 1137 INTEGER :: ib_bdy, ii, ij, igrd, ib, icount ! dummy loop indices 1138 INTEGER :: i_offset, j_offset, inn ! local integer 1139 INTEGER :: ibeg, iend ! local integer 1140 LOGICAL :: llnon, llson, llean, llwen ! local logicals indicating the presence of a ocean neighbour 1141 REAL(wp), POINTER, DIMENSION(:,:) :: pmask ! pointer to 2D mask fields 1142 REAL(wp) :: zefl, zwfl, znfl, zsfl ! local scalars 1143 CHARACTER(LEN=1), DIMENSION(jpbgrd) :: cgrid 1144 REAL(wp) , DIMENSION(jpi,jpj) :: ztmp 1145 !!---------------------------------------------------------------------- 1146 1147 cgrid = (/'t','u','v'/) 1148 989 1149 DO ib_bdy = 1, nb_bdy ! Indices and directions of rim velocity components 990 icount = 0991 1150 992 1151 ! Calculate relationship of U direction to the local orientation of the boundary … … 994 1153 ! flagu = 0 : u is tangential 995 1154 ! flagu = 1 : u is normal to the boundary and is direction is inward 996 997 1155 DO igrd = 1, jpbgrd 998 1156 SELECT CASE( igrd ) 999 CASE( 1 ) ; pmask => umask (:,:,1); i_offset = 01000 CASE( 2 ) ; pmask => bdytmask (:,:); i_offset = 11001 CASE( 3 ) ; pmask => zfmask (:,:); i_offset = 01157 CASE( 1 ) ; pmask => zumask ; i_offset = 0 1158 CASE( 2 ) ; pmask => bdytmask ; i_offset = 1 1159 CASE( 3 ) ; pmask => zfmask ; i_offset = 0 1002 1160 END SELECT 1003 1161 icount = 0 1004 1162 ztmp(:,:) = 0._wp 1005 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 1163 IF( lrim0 ) THEN ! extent of rim 0 1164 ibeg = 1 ; iend = idx_bdy(ib_bdy)%nblenrim0(igrd) 1165 ELSE ! extent of rim 1 1166 ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1 ; iend = idx_bdy(ib_bdy)%nblenrim(igrd) 1167 END IF 1168 DO ib = ibeg, iend 1006 1169 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 1007 1170 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 1008 1171 IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE 1009 z efl = pmask(ii+i_offset-1,ij)1010 z wfl = pmask(ii+i_offset ,ij)1172 zwfl = pmask(ii+i_offset-1,ij) 1173 zefl = pmask(ii+i_offset ,ij) 1011 1174 ! This error check only works if you are using the bdyXmask arrays 1012 1175 IF( i_offset == 1 .and. zefl + zwfl == 2. ) THEN … … 1014 1177 IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij) 1015 1178 ELSE 1016 ztmp(ii,ij) = -z efl + zwfl1179 ztmp(ii,ij) = -zwfl + zefl 1017 1180 ENDIF 1018 1181 END DO … … 1024 1187 ENDIF 1025 1188 CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 1026 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)1189 DO ib = ibeg, iend 1027 1190 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 1028 1191 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) … … 1035 1198 ! flagv = 0 : v is tangential 1036 1199 ! flagv = 1 : v is normal to the boundary and is direction is inward 1037 1038 1200 DO igrd = 1, jpbgrd 1039 1201 SELECT CASE( igrd ) 1040 CASE( 1 ) ; pmask => vmask (:,:,1); j_offset = 01041 CASE( 2 ) ; pmask => zfmask (:,:); j_offset = 01042 CASE( 3 ) ; pmask => bdytmask 1202 CASE( 1 ) ; pmask => zvmask ; j_offset = 0 1203 CASE( 2 ) ; pmask => zfmask ; j_offset = 0 1204 CASE( 3 ) ; pmask => bdytmask ; j_offset = 1 1043 1205 END SELECT 1044 1206 icount = 0 1045 1207 ztmp(:,:) = 0._wp 1046 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 1208 IF( lrim0 ) THEN ! extent of rim 0 1209 ibeg = 1 ; iend = idx_bdy(ib_bdy)%nblenrim0(igrd) 1210 ELSE ! extent of rim 1 1211 ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1 ; iend = idx_bdy(ib_bdy)%nblenrim(igrd) 1212 END IF 1213 DO ib = ibeg, iend 1047 1214 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 1048 1215 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 1049 1216 IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE 1050 z nfl = pmask(ii,ij+j_offset-1)1051 z sfl = pmask(ii,ij+j_offset )1217 zsfl = pmask(ii,ij+j_offset-1) 1218 znfl = pmask(ii,ij+j_offset ) 1052 1219 ! This error check only works if you are using the bdyXmask arrays 1053 1220 IF( j_offset == 1 .and. znfl + zsfl == 2. ) THEN … … 1055 1222 icount = icount + 1 1056 1223 ELSE 1057 ztmp(ii,ij) = -z nfl + zsfl1224 ztmp(ii,ij) = -zsfl + znfl 1058 1225 END IF 1059 1226 END DO … … 1065 1232 ENDIF 1066 1233 CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 1067 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)1234 DO ib = ibeg, iend 1068 1235 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 1069 1236 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) … … 1072 1239 END DO 1073 1240 ! 1074 END DO 1241 END DO ! ib_bdy 1075 1242 1076 1243 DO ib_bdy = 1, nb_bdy … … 1082 1249 END SELECT 1083 1250 ztmp(:,:) = 0._wp 1084 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 1251 IF( lrim0 ) THEN ! extent of rim 0 1252 ibeg = 1 ; iend = idx_bdy(ib_bdy)%nblenrim0(igrd) 1253 ELSE ! extent of rim 1 1254 ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1 ; iend = idx_bdy(ib_bdy)%nblenrim(igrd) 1255 END IF 1256 DO ib = ibeg, iend 1085 1257 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 1086 1258 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 1087 IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE 1088 ll ean = pmask(ii+1,ij ) == 1.1089 ll wen = pmask(ii-1,ij ) == 1.1090 ll non = pmask(ii ,ij+1) == 1.1091 ll son = pmask(ii ,ij-1) == 1.1259 IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE 1260 llnon = pmask(ii ,ij+1) == 1. 1261 llson = pmask(ii ,ij-1) == 1. 1262 llean = pmask(ii+1,ij ) == 1. 1263 llwen = pmask(ii-1,ij ) == 1. 1092 1264 inn = COUNT( (/ llnon, llson, llean, llwen /) ) 1093 1265 IF( inn == 0 ) THEN ! no neighbours -> interior of a corner or cluster of rim points … … 1101 1273 ELSE 1102 1274 WRITE(ctmp1,*) ' E R R O R : Problem with ',cgrid(igrd) ,' grid point', ii, ij, & 1103 ' on boundary set ', ib_bdy, ' has no free ocean neighbour' 1104 WRITE(ctmp2,*) ' There seems to be a cluster of rim points.' 1275 ' on boundary set ', ib_bdy, ' has no free ocean neighbour' 1276 IF( lrim0 ) THEN 1277 WRITE(ctmp2,*) ' There seems to be a cluster of rim 0 points.' 1278 ELSE 1279 WRITE(ctmp2,*) ' There seems to be a cluster of rim 1 points.' 1280 END IF 1105 1281 WRITE(ctmp3,*) ' ========== ' 1106 1282 CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ' ) … … 1142 1318 END DO 1143 1319 CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 1144 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)1320 DO ib = ibeg, iend 1145 1321 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 1146 1322 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) … … 1150 1326 END DO 1151 1327 1152 1153 ALLOCATE( lsend_bdyint(nb_bdy,jpbgrd,4), lrecv_bdyint(nb_bdy,jpbgrd,4) ) 1154 lsend_bdyint(:,:,:) = .false. 1155 lrecv_bdyint(:,:,:) = .false. 1156 ALLOCATE( lsend_bdyext(nb_bdy,jpbgrd,4), lrecv_bdyext(nb_bdy,jpbgrd,4) ) 1157 lsend_bdyext(:,:,:) = .false. 1158 lrecv_bdyext(:,:,:) = .false. 1159 ! 1160 ! Check which boundaries might need communication 1161 DO igrd = 1, jpbgrd 1162 DO ib_bdy = 1, nb_bdy 1163 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 1164 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 1165 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 1166 flagu = NINT(idx_bdy(ib_bdy)%flagu(ib,igrd)) 1167 flagv = NINT(idx_bdy(ib_bdy)%flagv(ib,igrd)) 1168 iibe = ii - flagu ! neighbouring point towards the exterior of the computational domain 1169 ijbe = ij - flagv 1170 iibi = ii + flagu ! neighbouring point towards the interior of the computational domain 1171 ijbi = ij + flagv 1172 CALL find_neib( ii, ij, idx_bdy(ib_bdy)%ntreat(ib,igrd), ii1, ij1, ii2, ij2, ii3, ij3 ) 1173 ! 1174 ! search neighbour in the west/east direction 1175 ! Rim is on the halo and computed ocean is towards exterior of mpi domain 1176 ! <-- (o exterior) --> 1177 ! (1) o|x OR (2) x|o 1178 ! |___ ___| 1179 IF( iibi == 0 .OR. ii1 == 0 .OR. ii2 == 0 .OR. ii3 == 0 ) lrecv_bdyint(ib_bdy,igrd,1) = .true. 1180 IF( iibi == jpi+1 .OR. ii1 == jpi+1 .OR. ii2 == jpi+1 .OR. ii3 == jpi+1 ) lrecv_bdyint(ib_bdy,igrd,2) = .true. 1181 IF( iibe == 0 ) lrecv_bdyext(ib_bdy,igrd,1) = .true. 1182 IF( iibe == jpi+1 ) lrecv_bdyext(ib_bdy,igrd,2) = .true. 1183 ! Check if neighbour has its rim parallel to its mpi subdomain border and located next to its halo 1184 ! :¨¨¨¨¨|¨¨--> | | <--¨¨|¨¨¨¨¨: 1185 ! : | x:o | neighbour limited by ... would need o | o:x | : 1186 ! :.....|_._:_____| (1) W neighbour E neighbour (2) |_____:_._|.....: 1187 IF( ii == 2 .AND. ( nbondi == 1 .OR. nbondi == 0 ) .AND. & 1188 & ( iibi == 3 .OR. ii1 == 3 .OR. ii2 == 3 .OR. ii3 == 3 ) ) lsend_bdyint(ib_bdy,igrd,1) = .true. 1189 IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. & 1190 & ( iibi == jpi-2 .OR. ii1 == jpi-2 .OR. ii2 == jpi-2 .OR. ii3 == jpi-2) ) lsend_bdyint(ib_bdy,igrd,2) = .true. 1191 IF( ii == 2 .AND. ( nbondi == 1 .OR. nbondi == 0 ) .AND. iibe == 3 ) lsend_bdyext(ib_bdy,igrd,1) = .true. 1192 IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. iibe == jpi-2 ) lsend_bdyext(ib_bdy,igrd,2) = .true. 1193 ! 1194 ! search neighbour in the north/south direction 1195 ! Rim is on the halo and computed ocean is towards exterior of mpi domain 1196 !(3) | | ^ ___o___ 1197 ! | |___x___| OR | | x | 1198 ! v o (4) | | 1199 IF( ijbi == 0 .OR. ij1 == 0 .OR. ij2 == 0 .OR. ij3 == 0 ) lrecv_bdyint(ib_bdy,igrd,3) = .true. 1200 IF( ijbi == jpj+1 .OR. ij1 == jpj+1 .OR. ij2 == jpj+1 .OR. ij3 == jpj+1 ) lrecv_bdyint(ib_bdy,igrd,4) = .true. 1201 IF( ijbe == 0 ) lrecv_bdyext(ib_bdy,igrd,3) = .true. 1202 IF( ijbe == jpj+1 ) lrecv_bdyext(ib_bdy,igrd,4) = .true. 1203 ! Check if neighbour has its rim parallel to its mpi subdomain _________ border and next to its halo 1204 ! ^ | o | : : 1205 ! | |¨¨¨¨x¨¨¨¨| neighbour limited by ... would need o | |....x....| 1206 ! :_________: (3) S neighbour N neighbour (4) v | o | 1207 IF( ij == 2 .AND. ( nbondj == 1 .OR. nbondj == 0 ) .AND. & 1208 & ( ijbi == 3 .OR. ij1 == 3 .OR. ij2 == 3 .OR. ij3 == 3 ) ) lsend_bdyint(ib_bdy,igrd,3) = .true. 1209 IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. & 1210 & ( ijbi == jpj-2 .OR. ij1 == jpj-2 .OR. ij2 == jpj-2 .OR. ij3 == jpj-2) ) lsend_bdyint(ib_bdy,igrd,4) = .true. 1211 IF( ij == 2 .AND. ( nbondj == 1 .OR. nbondj == 0 ) .AND. ijbe == 3 ) lsend_bdyext(ib_bdy,igrd,3) = .true. 1212 IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. ijbe == jpj-2 ) lsend_bdyext(ib_bdy,igrd,4) = .true. 1213 END DO 1214 END DO 1215 END DO 1216 1217 DO ib_bdy = 1,nb_bdy 1218 IF( cn_dyn2d(ib_bdy) == 'orlanski' .OR. cn_dyn2d(ib_bdy) == 'orlanski_npo' .OR. & 1219 & cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' .OR. & 1220 & cn_tra(ib_bdy) == 'orlanski' .OR. cn_tra(ib_bdy) == 'orlanski_npo' ) THEN 1221 DO igrd = 1, jpbgrd 1222 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 1223 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 1224 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 1225 IF( mig(ii) > 2 .AND. mig(ii) < jpiglo-2 .AND. mjg(ij) > 2 .AND. mjg(ij) < jpjglo-2 ) THEN 1226 WRITE(ctmp1,*) ' Orlanski can not be used when the open boundaries are on the interior of the computational domain' 1227 WRITE(ctmp2,*) ' ========== ' 1228 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 1229 END IF 1230 END DO 1231 END DO 1232 END IF 1233 END DO 1234 1235 ! 1236 ! Tidy up 1237 !-------- 1238 IF( nb_bdy>0 ) DEALLOCATE( nbidta, nbjdta, nbrdta ) 1239 ! 1240 END SUBROUTINE bdy_segs 1241 1242 1243 SUBROUTINE find_neib( ii, ij, itreat, ii1, ij1, ii2, ij2, ii3, ij3 ) 1244 !!---------------------------------------------------------------------- 1245 !! *** ROUTINE find_neib *** 1246 !! 1247 !! ** Purpose : get ii1, ij1, ii2, ij2, ii3, ij3, the indices of 1248 !! the free ocean neighbours of (ii,ij) for bdy treatment 1249 !! 1250 !!---------------------------------------------------------------------- 1251 INTEGER, INTENT(in ) :: ii, ij, itreat 1252 INTEGER, INTENT( out) :: ii1, ij1, ii2, ij2, ii3, ij3 1253 !!---------------------------------------------------------------------- 1254 SELECT CASE( itreat ) ! points that will be used by bdy routines, -1 will be discarded 1255 ! ! ! _____ ! _____ 1256 ! 1 | o ! 2 o | ! 3 | x ! 4 x | 1257 ! |_x_ _ ! _ _x_| ! | o ! o | 1258 CASE( 1 ) ; ii1 = ii+1 ; ij1 = ij+1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1 1259 CASE( 2 ) ; ii1 = ii-1 ; ij1 = ij+1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1 1260 CASE( 3 ) ; ii1 = ii+1 ; ij1 = ij-1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1 1261 CASE( 4 ) ; ii1 = ii-1 ; ij1 = ij-1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1 1262 ! | ! | ! o ! ______ ! or incomplete corner 1263 ! 5 | x o ! 6 o x | ! 7 __x__ ! 8 x ! 7 ____ o 1264 ! | ! | ! ! o ! |x___ 1265 CASE( 5 ) ; ii1 = ii+1 ; ij1 = ij ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1 1266 CASE( 6 ) ; ii1 = ii-1 ; ij1 = ij ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1 1267 CASE( 7 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1 1268 CASE( 8 ) ; ii1 = ii ; ij1 = ij-1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1 1269 ! o ! o ! _____| ! |_____ 1270 ! 9 ____x o ! 10 o x___ ! 11 x o ! 12 o x 1271 ! | ! | ! o ! o 1272 CASE( 9 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = ii+1 ; ij2 = ij ; ii3 = -1 ; ij3 = -1 1273 CASE( 10 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = ii-1 ; ij2 = ij ; ii3 = -1 ; ij3 = -1 1274 CASE( 11 ) ; ii1 = ii ; ij1 = ij-1 ; ii2 = ii+1 ; ij2 = ij ; ii3 = -1 ; ij3 = -1 1275 CASE( 12 ) ; ii1 = ii ; ij1 = ij-1 ; ii2 = ii-1 ; ij2 = ij ; ii3 = -1 ; ij3 = -1 1276 ! |_ o ! o _| ! ¨¨|_|¨¨ ! o 1277 ! 13 _| x o ! 14 o x |_ ! 15 o x o ! 16 o x o 1278 ! | o ! o | ! o ! __|¨|__ 1279 CASE( 13 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = ii+1 ; ij2 = ij ; ii3 = ii ; ij3 = ij-1 1280 CASE( 14 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = ii-1 ; ij2 = ij ; ii3 = ii ; ij3 = ij-1 1281 CASE( 15 ) ; ii1 = ii-1 ; ij1 = ij ; ii2 = ii ; ij2 = ij-1 ; ii3 = ii+1 ; ij3 = ij 1282 CASE( 16 ) ; ii1 = ii-1 ; ij1 = ij ; ii2 = ii ; ij2 = ij+1 ; ii3 = ii+1 ; ij3 = ij 1283 END SELECT 1284 END SUBROUTINE find_neib 1285 1328 END SUBROUTINE bdy_rim_treat 1329 1330 1331 SUBROUTINE find_neib( ii, ij, itreat, ii1, ij1, ii2, ij2, ii3, ij3 ) 1332 !!---------------------------------------------------------------------- 1333 !! *** ROUTINE find_neib *** 1334 !! 1335 !! ** Purpose : get ii1, ij1, ii2, ij2, ii3, ij3, the indices of 1336 !! the free ocean neighbours of (ii,ij) for bdy treatment 1337 !! 1338 !! ** Method : use itreat input to select a case 1339 !! N.B. ntreat is defined for all bdy points in routine bdy_segs 1340 !! 1341 !!---------------------------------------------------------------------- 1342 INTEGER, INTENT(in ) :: ii, ij, itreat 1343 INTEGER, INTENT( out) :: ii1, ij1, ii2, ij2, ii3, ij3 1344 !!---------------------------------------------------------------------- 1345 SELECT CASE( itreat ) ! points that will be used by bdy routines, -1 will be discarded 1346 ! ! ! _____ ! _____ 1347 ! 1 | o ! 2 o | ! 3 | x ! 4 x | 1348 ! |_x_ _ ! _ _x_| ! | o ! o | 1349 CASE( 1 ) ; ii1 = ii+1 ; ij1 = ij+1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1 1350 CASE( 2 ) ; ii1 = ii-1 ; ij1 = ij+1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1 1351 CASE( 3 ) ; ii1 = ii+1 ; ij1 = ij-1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1 1352 CASE( 4 ) ; ii1 = ii-1 ; ij1 = ij-1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1 1353 ! | ! | ! o ! ______ ! or incomplete corner 1354 ! 5 | x o ! 6 o x | ! 7 __x__ ! 8 x ! 7 ____ o 1355 ! | ! | ! ! o ! |x___ 1356 CASE( 5 ) ; ii1 = ii+1 ; ij1 = ij ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1 1357 CASE( 6 ) ; ii1 = ii-1 ; ij1 = ij ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1 1358 CASE( 7 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1 1359 CASE( 8 ) ; ii1 = ii ; ij1 = ij-1 ; ii2 = -1 ; ij2 = -1 ; ii3 = -1 ; ij3 = -1 1360 ! o ! o ! _____| ! |_____ 1361 ! 9 ____x o ! 10 o x___ ! 11 x o ! 12 o x 1362 ! | ! | ! o ! o 1363 CASE( 9 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = ii+1 ; ij2 = ij ; ii3 = -1 ; ij3 = -1 1364 CASE( 10 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = ii-1 ; ij2 = ij ; ii3 = -1 ; ij3 = -1 1365 CASE( 11 ) ; ii1 = ii ; ij1 = ij-1 ; ii2 = ii+1 ; ij2 = ij ; ii3 = -1 ; ij3 = -1 1366 CASE( 12 ) ; ii1 = ii ; ij1 = ij-1 ; ii2 = ii-1 ; ij2 = ij ; ii3 = -1 ; ij3 = -1 1367 ! |_ o ! o _| ! ¨¨|_|¨¨ ! o 1368 ! 13 _| x o ! 14 o x |_ ! 15 o x o ! 16 o x o 1369 ! | o ! o | ! o ! __|¨|__ 1370 CASE( 13 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = ii+1 ; ij2 = ij ; ii3 = ii ; ij3 = ij-1 1371 CASE( 14 ) ; ii1 = ii ; ij1 = ij+1 ; ii2 = ii-1 ; ij2 = ij ; ii3 = ii ; ij3 = ij-1 1372 CASE( 15 ) ; ii1 = ii-1 ; ij1 = ij ; ii2 = ii ; ij2 = ij-1 ; ii3 = ii+1 ; ij3 = ij 1373 CASE( 16 ) ; ii1 = ii-1 ; ij1 = ij ; ii2 = ii ; ij2 = ij+1 ; ii3 = ii+1 ; ij3 = ij 1374 END SELECT 1375 END SUBROUTINE find_neib 1376 1286 1377 1287 1378 SUBROUTINE bdy_ctl_seg
Note: See TracChangeset
for help on using the changeset viewer.
