Changeset 11258
- Timestamp:
- 2019-07-11T18:21:02+02:00 (5 years ago)
- Location:
- NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdydyn2d.F90
r11224 r11258 199 199 ELSE ; ibeg = idx%nblenrim0(igrd)+1 ; iend = idx%nblenrim(igrd) 200 200 END IF 201 sshdta(:,:) = 0.0201 ! 202 202 DO jb = ibeg, iend 203 203 ii = idx%nbi(jb,igrd) -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyini.F90
r11223 r11258 116 116 ! Open boundaries definition (arrays and masks) 117 117 CALL bdy_def 118 IF( ln_meshmask ) CALL bdy_meshwri() 118 119 ! 119 120 ! Open boundaries initialisation of external data arrays … … 146 147 INTEGER :: ib_bdy, ii, ij, igrd, ib, ir, iseg ! dummy loop indices 147 148 INTEGER :: icount, icountr, icountr0, ibr_max ! local integers 148 INTEGER :: ilen1 , ibm1! - -149 INTEGER :: ilen1 ! - - 149 150 INTEGER :: iwe, ies, iso, ino, inum, id_dummy ! - - 150 151 INTEGER :: jpbdta ! - - … … 153 154 INTEGER :: iibe, ijbe, iibi, ijbi ! - - 154 155 INTEGER :: flagu, flagv ! short cuts 155 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zz_read ! work space for 2D global boundary data 156 REAL(wp), POINTER, DIMENSION(:,:) :: pmask ! pointer to 2D mask fields 157 INTEGER, DIMENSION (4) :: kdimsz 158 INTEGER, DIMENSION(jpbgrd,jp_bdy) :: nblendta ! Length of index arrays 159 INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbidta, nbjdta ! Index arrays: i and j indices of bdy dta 160 INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbrdta ! Discrete distance from rim points 161 CHARACTER(LEN=1),DIMENSION(jpbgrd) :: cgrid 162 REAL(wp), DIMENSION(jpi,jpj) :: zfmask ! temporary fmask array excluding coastal boundary condition (shlat) 163 REAL(wp), DIMENSION(jpi,jpj) :: ztmask, zumask, zvmask ! temporary u/v mask array 164 !! 165 INTEGER :: nbdyind, nbdybeg, nbdyend 156 INTEGER :: nbdyind, nbdybeg, nbdyend 157 INTEGER , DIMENSION(4) :: kdimsz 158 INTEGER , DIMENSION(jpbgrd,jp_bdy) :: nblendta ! Length of index arrays 159 INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbidta, nbjdta ! Index arrays: i and j indices of bdy dta 160 INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: nbrdta ! Discrete distance from rim points 161 CHARACTER(LEN=1) , DIMENSION(jpbgrd) :: cgrid 162 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zz_read ! work space for 2D global boundary data 163 REAL(wp), POINTER , DIMENSION(:,:) :: zmask ! pointer to 2D mask fields 164 REAL(wp) , DIMENSION(jpi,jpj) :: zfmask ! temporary fmask array excluding coastal boundary condition (shlat) 165 REAL(wp) , DIMENSION(jpi,jpj) :: ztmask, zumask, zvmask ! temporary u/v mask array 166 166 !!---------------------------------------------------------------------- 167 167 ! … … 483 483 END DO 484 484 END DO 485 485 ! 486 486 ! Find lenght of boundaries and rim on local mpi domain 487 487 !------------------------------------------------------ … … 502 502 DO ib = 1, nblendta(igrd,ib_bdy) 503 503 ! check that data is in correct order in file 504 ibm1 = MAX(1,ib-1) 505 IF(lwp) THEN ! Since all procs read global data only need to do this check on one proc... 506 IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN 504 IF( ib > 1 ) THEN 505 IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ib-1,igrd,ib_bdy) ) THEN 507 506 CALL ctl_stop('bdy_segs : ERROR : boundary data in file must be defined ', & 508 507 & ' in order of distance from edge nbr A utility for re-ordering ', & … … 636 635 ! For the flagu/flagv calculation below we require a version of fmask without 637 636 ! the land boundary condition (shlat) included: 638 DO ij = 2, jpjm1639 DO ii = 2, jpim1640 zfmask(ii,ij) = ztmask(ii,ij ) * ztmask(ii+1,ij ) &641 & 637 DO ij = 1, jpjm1 638 DO ii = 1, jpim1 639 zfmask(ii,ij) = ztmask(ii,ij ) * ztmask(ii+1,ij ) & 640 & * ztmask(ii,ij+1) * ztmask(ii+1,ij+1) 642 641 END DO 643 642 END DO … … 678 677 ! handle rim1, do as if rim 0 was land 679 678 ! ------------------------------------ 680 679 680 ! z[tuv]mask are now set to zero on rim 0 points: 681 681 DO ib_bdy = 1, nb_bdy 682 682 DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(1) ! extent of rim 0 … … 692 692 693 693 ! Recompute zfmask 694 DO ij = 2, jpjm1695 DO ii = 2, jpim1696 zfmask(ii,ij) = ztmask(ii,ij ) * ztmask(ii+1,ij ) &697 & 694 DO ij = 1, jpjm1 695 DO ii = 1, jpim1 696 zfmask(ii,ij) = ztmask(ii,ij ) * ztmask(ii+1,ij ) & 697 & * ztmask(ii,ij+1) * ztmask(ii+1,ij+1) 698 698 END DO 699 699 END DO … … 714 714 715 715 CALL bdy_rim_treat( zumask, zvmask, zfmask, .false. ) ! compute flagu, flagv, ntreat on rim 1 716 717 716 ! 718 717 ! Check which boundaries might need communication … … 727 726 DO ib_bdy = 1, nb_bdy 728 727 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 728 IF( idx_bdy(ib_bdy)%ntreat(ib,igrd) == -1 ) CYCLE 729 729 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 730 730 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 731 IF( ib .LE. idx_bdy(ib_bdy)%nblenrim0(igrd) ) THEN ; ir = 0 732 ELSE ; ir = 1 733 END IF 731 ir = idx_bdy(ib_bdy)%nbr(ib,igrd) 734 732 flagu = NINT(idx_bdy(ib_bdy)%flagu(ib,igrd)) 735 733 flagv = NINT(idx_bdy(ib_bdy)%flagv(ib,igrd)) … … 784 782 785 783 DO ib_bdy = 1,nb_bdy 786 IF( cn_dyn2d(ib_bdy) == 'orlanski' .OR. cn_dyn2d(ib_bdy) == 'orlanski_npo' .OR. &784 IF( cn_dyn2d(ib_bdy) == 'orlanski' .OR. cn_dyn2d(ib_bdy) == 'orlanski_npo' .OR. & 787 785 & cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' .OR. & 788 786 & cn_tra(ib_bdy) == 'orlanski' .OR. cn_tra(ib_bdy) == 'orlanski_npo' ) THEN … … 793 791 IF( mig(ii) > 2 .AND. mig(ii) < jpiglo-2 .AND. mjg(ij) > 2 .AND. mjg(ij) < jpjglo-2 ) THEN 794 792 WRITE(ctmp1,*) ' Orlanski is not safe when the open boundaries are on the interior of the computational domain' 795 WRITE(ctmp2,*) ' ========== ' 796 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 793 CALL ctl_stop( ctmp1 ) 797 794 END IF 798 795 END DO … … 800 797 END IF 801 798 END DO 802 803 ! 804 ! Tidy up 805 !-------- 799 ! 806 800 DEALLOCATE( nbidta, nbjdta, nbrdta ) 807 801 ! … … 809 803 810 804 811 SUBROUTINE bdy_rim_treat( zumask, zvmask, zfmask, lrim0 )805 SUBROUTINE bdy_rim_treat( pumask, pvmask, pfmask, lrim0 ) 812 806 !!---------------------------------------------------------------------- 813 807 !! *** ROUTINE bdy_rim_treat *** … … 826 820 !! - and look at the ocean neighbours to compute ntreat 827 821 !!---------------------------------------------------------------------- 828 REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in ) :: zfmask ! temporary fmask excluding coastal boundary condition (shlat)829 REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in ) :: zumask, zvmask ! temporary t/u/v mask array822 REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in ) :: pfmask ! temporary fmask excluding coastal boundary condition (shlat) 823 REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in ) :: pumask, pvmask ! temporary t/u/v mask array 830 824 LOGICAL , INTENT (in ) :: lrim0 ! .true. -> rim 0 .false. -> rim 1 831 825 INTEGER :: ib_bdy, ii, ij, igrd, ib, icount ! dummy loop indices … … 833 827 INTEGER :: ibeg, iend ! local integer 834 828 LOGICAL :: llnon, llson, llean, llwen ! local logicals indicating the presence of a ocean neighbour 835 REAL(wp), POINTER, DIMENSION(:,:) :: pmask ! pointer to 2D mask fields829 REAL(wp), POINTER, DIMENSION(:,:) :: zmask ! pointer to 2D mask fields 836 830 REAL(wp) :: zefl, zwfl, znfl, zsfl ! local scalars 837 831 CHARACTER(LEN=1), DIMENSION(jpbgrd) :: cgrid … … 849 843 DO igrd = 1, jpbgrd 850 844 SELECT CASE( igrd ) 851 CASE( 1 ) ; pmask => zumask ; i_offset = 0852 CASE( 2 ) ; pmask => bdytmask ; i_offset = 1853 CASE( 3 ) ; pmask => zfmask ; i_offset = 0845 CASE( 1 ) ; zmask => pumask ; i_offset = 0 846 CASE( 2 ) ; zmask => bdytmask ; i_offset = 1 847 CASE( 3 ) ; zmask => pfmask ; i_offset = 0 854 848 END SELECT 855 849 icount = 0 856 ztmp(:,:) = 0._wp850 ztmp(:,:) = -999._wp 857 851 IF( lrim0 ) THEN ! extent of rim 0 858 852 ibeg = 1 ; iend = idx_bdy(ib_bdy)%nblenrim0(igrd) … … 864 858 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 865 859 IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE 866 zwfl = pmask(ii+i_offset-1,ij)867 zefl = pmask(ii+i_offset ,ij)860 zwfl = zmask(ii+i_offset-1,ij) 861 zefl = zmask(ii+i_offset ,ij) 868 862 ! This error check only works if you are using the bdyXmask arrays 869 863 IF( i_offset == 1 .and. zefl + zwfl == 2. ) THEN … … 875 869 END DO 876 870 IF( icount /= 0 ) THEN 877 WRITE(ctmp1,*) ' E R R O R :Some ',cgrid(igrd),' grid points,', &871 WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,', & 878 872 ' are not boundary points (flagu calculation). Check nbi, nbj, indices for boundary set ',ib_bdy 879 WRITE(ctmp2,*) ' ========== ' 880 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 873 CALL ctl_stop( ctmp1 ) 881 874 ENDIF 882 CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 875 SELECT CASE( igrd ) 876 CASE( 1 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 877 CASE( 2 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'U', 1. ) 878 CASE( 3 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'V', 1. ) 879 END SELECT 883 880 DO ib = ibeg, iend 884 881 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) … … 894 891 DO igrd = 1, jpbgrd 895 892 SELECT CASE( igrd ) 896 CASE( 1 ) ; pmask => zvmask ; j_offset = 0897 CASE( 2 ) ; pmask => zfmask ; j_offset = 0898 CASE( 3 ) ; pmask => bdytmask ; j_offset = 1893 CASE( 1 ) ; zmask => pvmask ; j_offset = 0 894 CASE( 2 ) ; zmask => pfmask ; j_offset = 0 895 CASE( 3 ) ; zmask => bdytmask ; j_offset = 1 899 896 END SELECT 900 897 icount = 0 901 ztmp(:,:) = 0._wp898 ztmp(:,:) = -999._wp 902 899 IF( lrim0 ) THEN ! extent of rim 0 903 900 ibeg = 1 ; iend = idx_bdy(ib_bdy)%nblenrim0(igrd) … … 909 906 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 910 907 IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE 911 zsfl = pmask(ii,ij+j_offset-1)912 znfl = pmask(ii,ij+j_offset )908 zsfl = zmask(ii,ij+j_offset-1) 909 znfl = zmask(ii,ij+j_offset ) 913 910 ! This error check only works if you are using the bdyXmask arrays 914 911 IF( j_offset == 1 .and. znfl + zsfl == 2. ) THEN … … 920 917 END DO 921 918 IF( icount /= 0 ) THEN 922 WRITE(ctmp1,*) ' E R R O R :Some ',cgrid(igrd),' grid points,', &919 WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,', & 923 920 ' are not boundary points (flagv calculation). Check nbi, nbj, indices for boundary set ',ib_bdy 924 WRITE(ctmp2,*) ' ========== ' 925 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 921 CALL ctl_stop( ctmp1 ) 926 922 ENDIF 927 CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 923 SELECT CASE( igrd ) 924 CASE( 1 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 925 CASE( 2 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'U', 1. ) 926 CASE( 3 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'V', 1. ) 927 END SELECT 928 928 DO ib = ibeg, iend 929 929 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) … … 938 938 DO igrd = 1, jpbgrd 939 939 SELECT CASE( igrd ) 940 CASE( 1 ) ; pmask => bdytmask941 CASE( 2 ) ; pmask => bdyumask942 CASE( 3 ) ; pmask => bdyvmask940 CASE( 1 ) ; zmask => bdytmask 941 CASE( 2 ) ; zmask => bdyumask 942 CASE( 3 ) ; zmask => bdyvmask 943 943 END SELECT 944 ztmp(:,:) = 0._wp944 ztmp(:,:) = -999._wp 945 945 IF( lrim0 ) THEN ! extent of rim 0 946 946 ibeg = 1 ; iend = idx_bdy(ib_bdy)%nblenrim0(igrd) … … 952 952 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 953 953 IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE 954 llnon = pmask(ii ,ij+1) == 1.955 llson = pmask(ii ,ij-1) == 1.956 llean = pmask(ii+1,ij ) == 1.957 llwen = pmask(ii-1,ij ) == 1.954 llnon = zmask(ii ,ij+1) == 1. 955 llson = zmask(ii ,ij-1) == 1. 956 llean = zmask(ii+1,ij ) == 1. 957 llwen = zmask(ii-1,ij ) == 1. 958 958 inn = COUNT( (/ llnon, llson, llean, llwen /) ) 959 959 IF( inn == 0 ) THEN ! no neighbours -> interior of a corner or cluster of rim points … … 961 961 ! 1 | o ! 2 o | ! 3 | x ! 4 x | ! | | -> error 962 962 ! |_x_ _ ! _ _x_| ! | o ! o | ! |x_x| 963 IF( pmask(ii+1,ij+1) == 1. ) THEN ; ztmp(ii,ij) = 1.964 ELSEIF( pmask(ii-1,ij+1) == 1. ) THEN ; ztmp(ii,ij) = 2.965 ELSEIF( pmask(ii+1,ij-1) == 1. ) THEN ; ztmp(ii,ij) = 3.966 ELSEIF( pmask(ii-1,ij-1) == 1. ) THEN ; ztmp(ii,ij) = 4.967 ELSE 968 WRITE(ctmp1,*) ' E R R O R :Problem with ',cgrid(igrd) ,' grid point', ii, ij, &963 IF( zmask(ii+1,ij+1) == 1. ) THEN ; ztmp(ii,ij) = 1. 964 ELSEIF( zmask(ii-1,ij+1) == 1. ) THEN ; ztmp(ii,ij) = 2. 965 ELSEIF( zmask(ii+1,ij-1) == 1. ) THEN ; ztmp(ii,ij) = 3. 966 ELSEIF( zmask(ii-1,ij-1) == 1. ) THEN ; ztmp(ii,ij) = 4. 967 ELSE ; ztmp(ii,ij) = -1. 968 WRITE(ctmp1,*) 'Problem with ',cgrid(igrd) ,' grid point', ii, ij, & 969 969 ' on boundary set ', ib_bdy, ' has no free ocean neighbour' 970 970 IF( lrim0 ) THEN … … 973 973 WRITE(ctmp2,*) ' There seems to be a cluster of rim 1 points.' 974 974 END IF 975 WRITE(ctmp3,*) ' ========== ' 976 CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ' ) 975 CALL ctl_warn( ctmp1, ctmp2 ) 977 976 END IF 978 977 END IF … … 1005 1004 END IF 1006 1005 IF( inn == 4 ) THEN 1007 WRITE(ctmp1,*) ' E R R O R :Problem with ',cgrid(igrd) ,' grid point', ii, ij, &1006 WRITE(ctmp1,*) 'Problem with ',cgrid(igrd) ,' grid point', ii, ij, & 1008 1007 ' on boundary set ', ib_bdy, ' have 4 neighbours' 1009 WRITE(ctmp2,*) ' ========== ' 1010 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 1008 CALL ctl_stop( ctmp1 ) 1011 1009 END IF 1012 1010 END DO 1013 CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 1011 SELECT CASE( igrd ) 1012 CASE( 1 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 1013 CASE( 2 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'U', 1. ) 1014 CASE( 3 ) ; CALL lbc_lnk( 'bdyini', ztmp, 'V', 1. ) 1015 END SELECT 1014 1016 DO ib = ibeg, iend 1015 1017 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) … … 1031 1033 !! 1032 1034 !! ** Method : use itreat input to select a case 1033 !! N.B. ntreat is defined for all bdy points in routine bdy_ segs1035 !! N.B. ntreat is defined for all bdy points in routine bdy_rim_treat 1034 1036 !! 1035 1037 !!---------------------------------------------------------------------- … … 1243 1245 icorns(ib2,1) = npckgw(ib1) 1244 1246 ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN 1245 WRITE(ctmp1,*) ' E R R O R :Found an acute open boundary corner at point (i,j)= ', &1247 WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', & 1246 1248 & jpisft(ib2), jpjwft(ib1) 1247 WRITE(ctmp2,*) ' ==========Not allowed yet'1248 WRITE(ctmp3,*) ' 1249 & 1250 CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ')1249 WRITE(ctmp2,*) ' Not allowed yet' 1250 WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), & 1251 & ' and South segment: ',npckgs(ib2) 1252 CALL ctl_stop( ctmp1, ctmp2, ctmp3 ) 1251 1253 ELSE 1252 WRITE(ctmp1,*) ' E R R O R :Check South and West Open boundary indices'1253 WRITE(ctmp2,*) ' ==========Crossing problem with West segment: ',npckgw(ib1) , &1254 & 1255 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ')1254 WRITE(ctmp1,*) ' Check South and West Open boundary indices' 1255 WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1) , & 1256 & ' and South segment: ',npckgs(ib2) 1257 CALL ctl_stop( ctmp1, ctmp2 ) 1256 1258 END IF 1257 1259 END IF … … 1275 1277 icorns(ib2,2) = npckge(ib1) 1276 1278 ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN 1277 WRITE(ctmp1,*) ' E R R O R :Found an acute open boundary corner at point (i,j)= ', &1279 WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', & 1278 1280 & jpisdt(ib2), jpjeft(ib1) 1279 WRITE(ctmp2,*) ' ==========Not allowed yet'1280 WRITE(ctmp3,*) ' 1281 & 1282 CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ')1281 WRITE(ctmp2,*) ' Not allowed yet' 1282 WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), & 1283 & ' and South segment: ',npckgs(ib2) 1284 CALL ctl_stop( ctmp1, ctmp2, ctmp3 ) 1283 1285 ELSE 1284 WRITE(ctmp1,*) ' E R R O R :Check South and East Open boundary indices'1285 WRITE(ctmp2,*) ' ==========Crossing problem with East segment: ',npckge(ib1), &1286 & 1287 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ')1286 WRITE(ctmp1,*) ' Check South and East Open boundary indices' 1287 WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), & 1288 & ' and South segment: ',npckgs(ib2) 1289 CALL ctl_stop( ctmp1, ctmp2 ) 1288 1290 END IF 1289 1291 END IF … … 1307 1309 icornn(ib2,1) = npckgw(ib1) 1308 1310 ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN 1309 WRITE(ctmp1,*) ' E R R O R :Found an acute open boundary corner at point (i,j)= ', &1311 WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', & 1310 1312 & jpinft(ib2), jpjwdt(ib1) 1311 WRITE(ctmp2,*) ' ==========Not allowed yet'1312 WRITE(ctmp3,*) ' 1313 & 1314 CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ')1313 WRITE(ctmp2,*) ' Not allowed yet' 1314 WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), & 1315 & ' and North segment: ',npckgn(ib2) 1316 CALL ctl_stop( ctmp1, ctmp2, ctmp3 ) 1315 1317 ELSE 1316 WRITE(ctmp1,*) ' E R R O R :Check North and West Open boundary indices'1317 WRITE(ctmp2,*) ' ==========Crossing problem with West segment: ',npckgw(ib1), &1318 & 1319 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ')1318 WRITE(ctmp1,*) ' Check North and West Open boundary indices' 1319 WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1), & 1320 & ' and North segment: ',npckgn(ib2) 1321 CALL ctl_stop( ctmp1, ctmp2 ) 1320 1322 END IF 1321 1323 END IF … … 1339 1341 icornn(ib2,2) = npckge(ib1) 1340 1342 ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN 1341 WRITE(ctmp1,*) ' E R R O R :Found an acute open boundary corner at point (i,j)= ', &1343 WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', & 1342 1344 & jpindt(ib2), jpjedt(ib1) 1343 WRITE(ctmp2,*) ' ==========Not allowed yet'1344 WRITE(ctmp3,*) ' 1345 & 1346 CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ')1345 WRITE(ctmp2,*) ' Not allowed yet' 1346 WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), & 1347 & ' and North segment: ',npckgn(ib2) 1348 CALL ctl_stop( ctmp1, ctmp2, ctmp3 ) 1347 1349 ELSE 1348 WRITE(ctmp1,*) ' E R R O R :Check North and East Open boundary indices'1349 WRITE(ctmp2,*) ' ==========Crossing problem with East segment: ',npckge(ib1), &1350 & 1351 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ')1350 WRITE(ctmp1,*) ' Check North and East Open boundary indices' 1351 WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), & 1352 & ' and North segment: ',npckgn(ib2) 1353 CALL ctl_stop( ctmp1, ctmp2 ) 1352 1354 END IF 1353 1355 END IF … … 1375 1377 IF (ztestmask(1)==1) THEN 1376 1378 IF (icornw(ib,1)==0) THEN 1377 WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgw(ib) 1378 WRITE(ctmp2,*) ' ========== does not start on land or on a corner' 1379 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 1379 WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib) 1380 CALL ctl_stop( ctmp1, ' does not start on land or on a corner' ) 1380 1381 ELSE 1381 1382 ! This is a corner … … 1387 1388 IF (ztestmask(2)==1) THEN 1388 1389 IF (icornw(ib,2)==0) THEN 1389 WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgw(ib) 1390 WRITE(ctmp2,*) ' ========== does not end on land or on a corner' 1391 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 1390 WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib) 1391 CALL ctl_stop( ' ', ctmp1, ' does not end on land or on a corner' ) 1392 1392 ELSE 1393 1393 ! This is a corner … … 1415 1415 IF (ztestmask(1)==1) THEN 1416 1416 IF (icorne(ib,1)==0) THEN 1417 WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckge(ib) 1418 WRITE(ctmp2,*) ' ========== does not start on land or on a corner' 1419 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 1417 WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib) 1418 CALL ctl_stop( ctmp1, ' does not start on land or on a corner' ) 1420 1419 ELSE 1421 1420 ! This is a corner … … 1427 1426 IF (ztestmask(2)==1) THEN 1428 1427 IF (icorne(ib,2)==0) THEN 1429 WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckge(ib) 1430 WRITE(ctmp2,*) ' ========== does not end on land or on a corner' 1431 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 1428 WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib) 1429 CALL ctl_stop( ctmp1, ' does not end on land or on a corner' ) 1432 1430 ELSE 1433 1431 ! This is a corner … … 1454 1452 1455 1453 IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN 1456 WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgs(ib) 1457 WRITE(ctmp2,*) ' ========== does not start on land or on a corner' 1458 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 1454 WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib) 1455 CALL ctl_stop( ctmp1, ' does not start on land or on a corner' ) 1459 1456 ENDIF 1460 1457 IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN 1461 WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgs(ib) 1462 WRITE(ctmp2,*) ' ========== does not end on land or on a corner' 1463 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 1458 WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib) 1459 CALL ctl_stop( ctmp1, ' does not end on land or on a corner' ) 1464 1460 ENDIF 1465 1461 END DO … … 1480 1476 1481 1477 IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN 1482 WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgn(ib) 1483 WRITE(ctmp2,*) ' ========== does not start on land' 1484 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 1478 WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib) 1479 CALL ctl_stop( ctmp1, ' does not start on land' ) 1485 1480 ENDIF 1486 1481 IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN 1487 WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgn(ib) 1488 WRITE(ctmp2,*) ' ========== does not end on land' 1489 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 1482 WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib) 1483 CALL ctl_stop( ctmp1, ' does not end on land' ) 1490 1484 ENDIF 1491 1485 END DO … … 1727 1721 ! 1728 1722 IF( itest>0 ) THEN 1729 WRITE(ctmp1,*) ' E R R O R : Segments ', ib1, 'and ', ib2 1730 WRITE(ctmp2,*) ' ========== have different open bdy schemes' 1731 CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' ) 1723 WRITE(ctmp1,*) ' Segments ', ib1, 'and ', ib2 1724 CALL ctl_stop( ctmp1, ' have different open bdy schemes' ) 1732 1725 ENDIF 1733 1726 ! 1734 1727 END SUBROUTINE bdy_ctl_corn 1735 1728 1729 1730 SUBROUTINE bdy_meshwri() 1731 !!---------------------------------------------------------------------- 1732 !! *** ROUTINE bdy_meshwri *** 1733 !! 1734 !! ** Purpose : write netcdf file with nbr, flagu, flagv, ntreat for T, U 1735 !! and V points in 2D arrays for easier visualisation/control 1736 !! 1737 !! ** Method : use iom_rstput as in domwri.F 1738 !!---------------------------------------------------------------------- 1739 INTEGER :: ib_bdy, ii, ij, igrd, ib ! dummy loop indices 1740 INTEGER :: inum ! - - 1741 REAL(wp), POINTER, DIMENSION(:,:) :: zmask ! pointer to 2D mask fields 1742 REAL(wp) , DIMENSION(jpi,jpj) :: ztmp 1743 CHARACTER(LEN=1) , DIMENSION(jpbgrd) :: cgrid 1744 !!---------------------------------------------------------------------- 1745 cgrid = (/'t','u','v'/) 1746 CALL iom_open( 'bdy_mesh', inum, ldwrt = .TRUE. ) 1747 DO igrd = 1, jpbgrd 1748 SELECT CASE( igrd ) 1749 CASE( 1 ) ; zmask => tmask(:,:,1) 1750 CASE( 2 ) ; zmask => umask(:,:,1) 1751 CASE( 3 ) ; zmask => vmask(:,:,1) 1752 END SELECT 1753 ztmp(:,:) = zmask(:,:) 1754 DO ib_bdy = 1, nb_bdy 1755 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) ! nbr deined for all rims 1756 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 1757 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 1758 ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%nbr(ib,igrd), wp) + 10. 1759 IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij) 1760 END DO 1761 END DO 1762 CALL iom_rstput( 0, 0, inum, 'bdy_nbr_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 ) 1763 ztmp(:,:) = zmask(:,:) 1764 DO ib_bdy = 1, nb_bdy 1765 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) ! flagu defined only for rims 0 and 1 1766 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 1767 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 1768 ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagu(ib,igrd), wp) + 10. 1769 IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij) 1770 END DO 1771 END DO 1772 CALL iom_rstput( 0, 0, inum, 'flagu_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 ) 1773 ztmp(:,:) = zmask(:,:) 1774 DO ib_bdy = 1, nb_bdy 1775 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) ! flagv defined only for rims 0 and 1 1776 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 1777 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 1778 ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagv(ib,igrd), wp) + 10. 1779 IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij) 1780 END DO 1781 END DO 1782 CALL iom_rstput( 0, 0, inum, 'flagv_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 ) 1783 ztmp(:,:) = zmask(:,:) 1784 DO ib_bdy = 1, nb_bdy 1785 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) ! ntreat defined only for rims 0 and 1 1786 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 1787 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 1788 ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%ntreat(ib,igrd), wp) + 10. 1789 IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij) 1790 END DO 1791 END DO 1792 CALL iom_rstput( 0, 0, inum, 'ntreat_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 ) 1793 END DO 1794 CALL iom_close( inum ) 1795 1796 END SUBROUTINE bdy_meshwri 1797 1736 1798 !!================================================================================= 1737 1799 END MODULE bdyini -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdylib.F90
r11191 r11258 149 149 REAL(wp) :: zdy_1, zdy_2, zsign_ups 150 150 REAL(wp), PARAMETER :: zepsilon = 1.e-30 ! local small value 151 REAL(wp), POINTER, DIMENSION(:,:) :: pmask ! land/sea mask for field152 REAL(wp), POINTER, DIMENSION(:,:) :: pmask_xdif ! land/sea mask for x-derivatives153 REAL(wp), POINTER, DIMENSION(:,:) :: pmask_ydif ! land/sea mask for y-derivatives151 REAL(wp), POINTER, DIMENSION(:,:) :: zmask ! land/sea mask for field 152 REAL(wp), POINTER, DIMENSION(:,:) :: zmask_xdif ! land/sea mask for x-derivatives 153 REAL(wp), POINTER, DIMENSION(:,:) :: zmask_ydif ! land/sea mask for y-derivatives 154 154 REAL(wp), POINTER, DIMENSION(:,:) :: pe_xdif ! scale factors for x-derivatives 155 155 REAL(wp), POINTER, DIMENSION(:,:) :: pe_ydif ! scale factors for y-derivatives … … 162 162 SELECT CASE(igrd) 163 163 CASE(1) 164 pmask => tmask(:,:,1)165 pmask_xdif => umask(:,:,1)166 pmask_ydif => vmask(:,:,1)164 zmask => tmask(:,:,1) 165 zmask_xdif => umask(:,:,1) 166 zmask_ydif => vmask(:,:,1) 167 167 pe_xdif => e1u(:,:) 168 168 pe_ydif => e2v(:,:) … … 170 170 ij_offset = 0 171 171 CASE(2) 172 pmask => umask(:,:,1)173 pmask_xdif => tmask(:,:,1)174 pmask_ydif => fmask(:,:,1)172 zmask => umask(:,:,1) 173 zmask_xdif => tmask(:,:,1) 174 zmask_ydif => fmask(:,:,1) 175 175 pe_xdif => e1t(:,:) 176 176 pe_ydif => e2f(:,:) … … 178 178 ij_offset = 0 179 179 CASE(3) 180 pmask => vmask(:,:,1)181 pmask_xdif => fmask(:,:,1)182 pmask_ydif => tmask(:,:,1)180 zmask => vmask(:,:,1) 181 zmask_xdif => fmask(:,:,1) 182 zmask_ydif => tmask(:,:,1) 183 183 pe_xdif => e1f(:,:) 184 184 pe_ydif => e2t(:,:) … … 214 214 ! 215 215 ! Calculate scale factors for calculation of spatial derivatives. 216 zex1 = ( abs(iibm1-iibm2) * pe_xdif(iibm1 +ii_offset,ijbm1 )&217 & + abs(ijbm1-ijbm2) * pe_ydif(iibm1 ,ijbm1+ij_offset) )218 zex2 = ( abs(iibm1-iibm2) * pe_xdif(iibm2 +ii_offset,ijbm2 )&219 & + abs(ijbm1-ijbm2) * pe_ydif(iibm2 ,ijbm2+ij_offset) )220 zey1 = ( (iibm1-iibm1jm1) * pe_xdif(iibm1jm1+ii_offset,ijbm1jm1 ) &216 zex1 = ( abs(iibm1-iibm2) * pe_xdif(iibm1 +ii_offset,ijbm1 ) & 217 & + abs(ijbm1-ijbm2) * pe_ydif(iibm1 ,ijbm1 +ij_offset) ) 218 zex2 = ( abs(iibm1-iibm2) * pe_xdif(iibm2 +ii_offset,ijbm2 ) & 219 & + abs(ijbm1-ijbm2) * pe_ydif(iibm2 ,ijbm2 +ij_offset) ) 220 zey1 = ( (iibm1-iibm1jm1) * pe_xdif(iibm1jm1+ii_offset,ijbm1jm1 ) & 221 221 & + (ijbm1-ijbm1jm1) * pe_ydif(iibm1jm1 ,ijbm1jm1+ij_offset) ) 222 zey2 = ( (iibm1jp1-iibm1) * pe_xdif(iibm1 +ii_offset,ijbm1)&223 & + (ijbm1jp1-ijbm1) * pe_ydif(iibm1 ,ijbm1+ij_offset) )222 zey2 = ( (iibm1jp1-iibm1) * pe_xdif(iibm1 +ii_offset,ijbm1 ) & 223 & + (ijbm1jp1-ijbm1) * pe_ydif(iibm1 ,ijbm1 +ij_offset) ) 224 224 ! make sure scale factors are nonzero 225 225 if( zey1 .lt. rsmall ) zey1 = zey2 … … 228 228 zey1 = max(zey1,rsmall); zey2 = max(zey2,rsmall); 229 229 ! 230 ! Calculate masks for calculation of spatial derivatives. 231 zmask_x = ( abs(iibm1-iibm2) * pmask_xdif(iibm2+ii_offset,ijbm2 )&232 & + abs(ijbm1-ijbm2) * pmask_ydif(iibm2 ,ijbm2+ij_offset) )233 zmask_y1 = ( (iibm1-iibm1jm1) * pmask_xdif(iibm1jm1+ii_offset,ijbm1jm1 )&234 & + (ijbm1-ijbm1jm1) * pmask_ydif(iibm1jm1 ,ijbm1jm1+ij_offset) )235 zmask_y2 = ( (iibm1jp1-iibm1) * pmask_xdif(iibm1+ii_offset,ijbm1)&236 & + (ijbm1jp1-ijbm1) * pmask_ydif(iibm1 ,ijbm1+ij_offset) )230 ! Calculate masks for calculation of spatial derivatives. 231 zmask_x = ( abs(iibm1-iibm2) * zmask_xdif(iibm2 +ii_offset,ijbm2 ) & 232 & + abs(ijbm1-ijbm2) * zmask_ydif(iibm2 ,ijbm2 +ij_offset) ) 233 zmask_y1 = ( (iibm1-iibm1jm1) * zmask_xdif(iibm1jm1+ii_offset,ijbm1jm1 ) & 234 & + (ijbm1-ijbm1jm1) * zmask_ydif(iibm1jm1 ,ijbm1jm1+ij_offset) ) 235 zmask_y2 = ( (iibm1jp1-iibm1) * zmask_xdif(iibm1 +ii_offset,ijbm1 ) & 236 & + (ijbm1jp1-ijbm1) * zmask_ydif(iibm1 ,ijbm1 +ij_offset) ) 237 237 238 238 ! Calculation of terms required for both versions of the scheme. … … 242 242 ! Note no rdt factor in expression for zdt because it cancels in the expressions for 243 243 ! zrx and zry. 244 zdt = phia(iibm1,ijbm1) - phib(iibm1,ijbm1)245 zdx = ( ( phia(iibm1,ijbm1) - phia(iibm2,ijbm2) ) / zex2 ) * zmask_x244 zdt = phia(iibm1 ,ijbm1 ) - phib(iibm1 ,ijbm1 ) 245 zdx = ( ( phia(iibm1 ,ijbm1 ) - phia(iibm2 ,ijbm2 ) ) / zex2 ) * zmask_x 246 246 zdy_1 = ( ( phib(iibm1 ,ijbm1 ) - phib(iibm1jm1,ijbm1jm1) ) / zey1 ) * zmask_y1 247 zdy_2 = ( ( phib(iibm1jp1,ijbm1jp1) - phib(iibm1 ,ijbm1 )) / zey2 ) * zmask_y2247 zdy_2 = ( ( phib(iibm1jp1,ijbm1jp1) - phib(iibm1 ,ijbm1 ) ) / zey2 ) * zmask_y2 248 248 zdy_centred = 0.5 * ( zdy_1 + zdy_2 ) 249 249 !!$ zdy_centred = phib(iibm1jp1,ijbm1jp1) - phib(iibm1jm1,ijbm1jm1) … … 276 276 & + zwgt * ( phi_ext(jb) - phib(ii,ij) ) ) / ( 1. + zrx ) 277 277 end if 278 phia(ii,ij) = phia(ii,ij) * pmask(ii,ij)278 phia(ii,ij) = phia(ii,ij) * zmask(ii,ij) 279 279 END DO 280 280 ! … … 314 314 REAL(wp) :: zdy_1, zdy_2, zsign_ups 315 315 REAL(wp), PARAMETER :: zepsilon = 1.e-30 ! local small value 316 REAL(wp), POINTER, DIMENSION(:,:,:) :: pmask ! land/sea mask for field317 REAL(wp), POINTER, DIMENSION(:,:,:) :: pmask_xdif ! land/sea mask for x-derivatives318 REAL(wp), POINTER, DIMENSION(:,:,:) :: pmask_ydif ! land/sea mask for y-derivatives316 REAL(wp), POINTER, DIMENSION(:,:,:) :: zmask ! land/sea mask for field 317 REAL(wp), POINTER, DIMENSION(:,:,:) :: zmask_xdif ! land/sea mask for x-derivatives 318 REAL(wp), POINTER, DIMENSION(:,:,:) :: zmask_ydif ! land/sea mask for y-derivatives 319 319 REAL(wp), POINTER, DIMENSION(:,:) :: pe_xdif ! scale factors for x-derivatives 320 320 REAL(wp), POINTER, DIMENSION(:,:) :: pe_ydif ! scale factors for y-derivatives … … 327 327 SELECT CASE(igrd) 328 328 CASE(1) 329 pmask => tmask(:,:,:)330 pmask_xdif => umask(:,:,:)331 pmask_ydif => vmask(:,:,:)329 zmask => tmask(:,:,:) 330 zmask_xdif => umask(:,:,:) 331 zmask_ydif => vmask(:,:,:) 332 332 pe_xdif => e1u(:,:) 333 333 pe_ydif => e2v(:,:) … … 335 335 ij_offset = 0 336 336 CASE(2) 337 pmask => umask(:,:,:)338 pmask_xdif => tmask(:,:,:)339 pmask_ydif => fmask(:,:,:)337 zmask => umask(:,:,:) 338 zmask_xdif => tmask(:,:,:) 339 zmask_ydif => fmask(:,:,:) 340 340 pe_xdif => e1t(:,:) 341 341 pe_ydif => e2f(:,:) … … 343 343 ij_offset = 0 344 344 CASE(3) 345 pmask => vmask(:,:,:)346 pmask_xdif => fmask(:,:,:)347 pmask_ydif => tmask(:,:,:)345 zmask => vmask(:,:,:) 346 zmask_xdif => fmask(:,:,:) 347 zmask_ydif => tmask(:,:,:) 348 348 pe_xdif => e1f(:,:) 349 349 pe_ydif => e2t(:,:) … … 381 381 ! 382 382 ! Calculate scale factors for calculation of spatial derivatives. 383 zex1 = ( abs(iibm1-iibm2) * pe_xdif(iibm1 +ii_offset,ijbm1 )&384 & + abs(ijbm1-ijbm2) * pe_ydif(iibm1 ,ijbm1+ij_offset) )385 zex2 = ( abs(iibm1-iibm2) * pe_xdif(iibm2 +ii_offset,ijbm2 )&386 & + abs(ijbm1-ijbm2) * pe_ydif(iibm2 ,ijbm2+ij_offset) )387 zey1 = ( (iibm1-iibm1jm1) * pe_xdif(iibm1jm1+ii_offset,ijbm1jm1 ) &383 zex1 = ( abs(iibm1-iibm2) * pe_xdif(iibm1 +ii_offset,ijbm1 ) & 384 & + abs(ijbm1-ijbm2) * pe_ydif(iibm1 ,ijbm1+ij_offset ) ) 385 zex2 = ( abs(iibm1-iibm2) * pe_xdif(iibm2 +ii_offset,ijbm2 ) & 386 & + abs(ijbm1-ijbm2) * pe_ydif(iibm2 ,ijbm2+ij_offset ) ) 387 zey1 = ( (iibm1-iibm1jm1) * pe_xdif(iibm1jm1+ii_offset,ijbm1jm1 ) & 388 388 & + (ijbm1-ijbm1jm1) * pe_ydif(iibm1jm1 ,ijbm1jm1+ij_offset) ) 389 zey2 = ( (iibm1jp1-iibm1) * pe_xdif(iibm1 +ii_offset,ijbm1)&390 & + (ijbm1jp1-ijbm1) * pe_ydif(iibm1 ,ijbm1+ij_offset) )389 zey2 = ( (iibm1jp1-iibm1) * pe_xdif(iibm1 +ii_offset,ijbm1 ) & 390 & + (ijbm1jp1-ijbm1) * pe_ydif(iibm1 ,ijbm1+ij_offset ) ) 391 391 ! make sure scale factors are nonzero 392 392 if( zey1 .lt. rsmall ) zey1 = zey2 … … 396 396 ! 397 397 ! Calculate masks for calculation of spatial derivatives. 398 zmask_x = ( abs(iibm1-iibm2) * pmask_xdif(iibm2+ii_offset,ijbm2 ,jk)&399 & + abs(ijbm1-ijbm2) * pmask_ydif(iibm2 ,ijbm2+ij_offset,jk) )400 zmask_y1 = ( (iibm1-iibm1jm1) * pmask_xdif(iibm1jm1+ii_offset,ijbm1jm1 ,jk) &401 & + (ijbm1-ijbm1jm1) * pmask_ydif(iibm1jm1 ,ijbm1jm1+ij_offset,jk) )402 zmask_y2 = ( (iibm1jp1-iibm1) * pmask_xdif(iibm1+ii_offset,ijbm1 ,jk)&403 & + (ijbm1jp1-ijbm1) * pmask_ydif(iibm1 ,ijbm1+ij_offset,jk) )398 zmask_x = ( abs(iibm1-iibm2) * zmask_xdif(iibm2 +ii_offset,ijbm2 ,jk) & 399 & + abs(ijbm1-ijbm2) * zmask_ydif(iibm2 ,ijbm2 +ij_offset,jk) ) 400 zmask_y1 = ( (iibm1-iibm1jm1) * zmask_xdif(iibm1jm1+ii_offset,ijbm1jm1 ,jk) & 401 & + (ijbm1-ijbm1jm1) * zmask_ydif(iibm1jm1 ,ijbm1jm1+ij_offset,jk) ) 402 zmask_y2 = ( (iibm1jp1-iibm1) * zmask_xdif(iibm1 +ii_offset,ijbm1 ,jk) & 403 & + (ijbm1jp1-ijbm1) * zmask_ydif(iibm1 ,ijbm1 +ij_offset,jk) ) 404 404 ! 405 405 ! Calculate normal (zrx) and tangential (zry) components of radiation velocities. … … 407 407 ! Centred derivative is calculated as average of "left" and "right" derivatives for 408 408 ! this reason. 409 zdt = phia(iibm1,ijbm1,jk) - phib(iibm1,ijbm1,jk)410 zdx = ( ( phia(iibm1,ijbm1,jk) - phia(iibm2,ijbm2,jk) ) / zex2 ) * zmask_x409 zdt = phia(iibm1 ,ijbm1 ,jk) - phib(iibm1 ,ijbm1 ,jk) 410 zdx = ( ( phia(iibm1 ,ijbm1 ,jk) - phia(iibm2 ,ijbm2 ,jk) ) / zex2 ) * zmask_x 411 411 zdy_1 = ( ( phib(iibm1 ,ijbm1 ,jk) - phib(iibm1jm1,ijbm1jm1,jk) ) / zey1 ) * zmask_y1 412 412 zdy_2 = ( ( phib(iibm1jp1,ijbm1jp1,jk) - phib(iibm1 ,ijbm1 ,jk) ) / zey2 ) * zmask_y2 … … 442 442 & + zwgt * ( phi_ext(jb,jk) - phib(ii,ij,jk) ) ) / ( 1. + zrx ) 443 443 end if 444 phia(ii,ij,jk) = phia(ii,ij,jk) * pmask(ii,ij,jk)444 phia(ii,ij,jk) = phia(ii,ij,jk) * zmask(ii,ij,jk) 445 445 END DO 446 446 ! … … 469 469 !! 470 470 REAL(wp) :: zweight 471 REAL(wp), POINTER, DIMENSION(:,:,:) :: pmask ! land/sea mask for field471 REAL(wp), POINTER, DIMENSION(:,:,:) :: zmask ! land/sea mask for field 472 472 INTEGER :: ib, ik ! dummy loop indices 473 473 INTEGER :: ii, ij ! 2D addresses … … 480 480 ! 481 481 SELECT CASE(igrd) 482 CASE(1) ; pmask => tmask(:,:,:)483 CASE(2) ; pmask => umask(:,:,:)484 CASE(3) ; pmask => vmask(:,:,:)482 CASE(1) ; zmask => tmask(:,:,:) 483 CASE(2) ; zmask => umask(:,:,:) 484 CASE(3) ; zmask => vmask(:,:,:) 485 485 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for igrd in bdy_nmn' ) 486 486 END SELECT … … 502 502 IF( ii1 < 1 .OR. ii1 > jpi .OR. ij1 < 1 .OR. ij1 > jpj ) CYCLE 503 503 DO ik = 1, ipkm1 504 IF( pmask(ii1,ij1,ik) /= 0. ) phia(ii,ij,ik) = phia(ii1,ij1,ik)504 IF( zmask(ii1,ij1,ik) /= 0. ) phia(ii,ij,ik) = phia(ii1,ij1,ik) 505 505 END DO 506 506 CASE( 9:12 ) … … 508 508 IF( ii2 < 1 .OR. ii2 > jpi .OR. ij2 < 1 .OR. ij2 > jpj ) CYCLE 509 509 DO ik = 1, ipkm1 510 zweight = pmask(ii1,ij1,ik) + pmask(ii2,ij2,ik)510 zweight = zmask(ii1,ij1,ik) + zmask(ii2,ij2,ik) 511 511 IF( zweight /= 0. ) phia(ii,ij,ik) = ( phia(ii1,ij1,ik) + phia(ii2,ij2,ik) ) / zweight 512 512 END DO … … 516 516 IF( ii3 < 1 .OR. ii3 > jpi .OR. ij3 < 1 .OR. ij3 > jpj ) CYCLE 517 517 DO ik = 1, ipkm1 518 zweight = pmask(ii1,ij1,ik) + pmask(ii2,ij2,ik) + pmask(ii3,ij3,ik)518 zweight = zmask(ii1,ij1,ik) + zmask(ii2,ij2,ik) + zmask(ii3,ij3,ik) 519 519 IF( zweight /= 0. ) phia(ii,ij,ik) = ( phia(ii1,ij1,ik) + phia(ii2,ij2,ik) + phia(ii3,ij3,ik) ) / zweight 520 520 END DO -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DOM/domain.F90
r10425 r11258 101 101 CASE( 0 ) ; WRITE(numout,*) ' (i.e. closed)' 102 102 CASE( 1 ) ; WRITE(numout,*) ' (i.e. cyclic east-west)' 103 CASE( 2 ) ; WRITE(numout,*) ' (i.e. equatorial symmetric)'103 CASE( 2 ) ; WRITE(numout,*) ' (i.e. cyclic north-south)' 104 104 CASE( 3 ) ; WRITE(numout,*) ' (i.e. north fold with T-point pivot)' 105 105 CASE( 4 ) ; WRITE(numout,*) ' (i.e. cyclic east-west and north fold with T-point pivot)' … … 527 527 INTEGER :: inum, ii ! local integer 528 528 REAL(wp) :: zorca_res ! local scalars 529 REAL(wp) :: ziglo, zjglo, zkglo, zperio ! - - 529 REAL(wp) :: zperio ! - - 530 INTEGER, DIMENSION(4) :: idvar, idimsz ! size of dimensions 530 531 !!---------------------------------------------------------------------- 531 532 ! … … 559 560 ! 560 561 ENDIF 561 ! 562 CALL iom_get( inum, 'jpiglo', ziglo ) ; kpi = NINT( ziglo ) 563 CALL iom_get( inum, 'jpjglo', zjglo ) ; kpj = NINT( zjglo ) 564 CALL iom_get( inum, 'jpkglo', zkglo ) ; kpk = NINT( zkglo ) 562 ! 563 idvar = iom_varid( inum, 'e3t_0', kdimsz = idimsz ) ! use e3t_0, that must exist, to get jp(ijk)glo 564 kpi = idimsz(1) 565 kpj = idimsz(2) 566 kpk = idimsz(3) 565 567 CALL iom_get( inum, 'jperio', zperio ) ; kperio = NINT( zperio ) 566 568 CALL iom_close( inum ) -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DYN/dynkeg.F90
r11195 r11258 74 74 INTEGER, INTENT( in ) :: kscheme ! =0/1 type of KEG scheme 75 75 ! 76 INTEGER :: ji, jj, jk, jb ! dummy loop indices 77 INTEGER :: igrd, ib_bdy ! local integers 76 INTEGER :: ji, jj, jk ! dummy loop indices 78 77 REAL(wp) :: zu, zv ! local scalars 79 78 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhke 80 79 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv 81 REAL(wp) :: zweightu, zweightv82 LOGICAL, DIMENSION(4) :: llsend1, llrecv1 ! indicate how bdy communications are to be carried out83 80 !!---------------------------------------------------------------------- 84 81 ! … … 113 110 END DO 114 111 END DO 115 !116 IF (ln_bdy) THEN117 ! Maria Luneva & Fred Wobus: July-2016118 ! compensate for lack of turbulent kinetic energy on liquid bdy points119 DO ib_bdy = 1, nb_bdy120 IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN121 igrd = 1 ! compensating null velocity on the bdy122 DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)123 ji = idx_bdy(ib_bdy)%nbi(jb,igrd) ! maximum extent : from 1 to jpi124 jj = idx_bdy(ib_bdy)%nbj(jb,igrd) ! maximum extent : from 1 to jpj125 IF( ji == 1 .OR. jj == 1 ) CYCLE126 DO jk = 1, jpkm1127 zhke(ji,jj,jk) = 0._wp128 zweightu = umask(ji-1,jj ,jk) + umask(ji,jj,jk)129 zweightv = vmask(ji ,jj-1,jk) + vmask(ji,jj,jk)130 zu = un(ji-1,jj ,jk) * un(ji-1,jj ,jk) + un(ji ,jj ,jk) * un(ji ,jj ,jk)131 zv = vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) + vn(ji ,jj ,jk) * vn(ji ,jj ,jk)132 IF( zweightu > 0._wp ) zhke(ji,jj,jk) = zhke(ji,jj,jk) + zu / (2._wp * zweightu)133 IF( zweightv > 0._wp ) zhke(ji,jj,jk) = zhke(ji,jj,jk) + zv / (2._wp * zweightv)134 END DO135 END DO136 END IF137 END DO138 ! send jpi-1, jpj-1 and receive 1 used in the computation of the speed tendencies139 llsend1(:) = .false.140 llrecv1(:) = .false.141 DO ib_bdy = 1, nb_bdy142 llsend1(2) = llsend1(2) .OR. ANY(lsend_bdy(ib_bdy,igrd,2,:)) ! send east143 llsend1(4) = llsend1(4) .OR. ANY(lsend_bdy(ib_bdy,igrd,4,:)) ! send north144 llrecv1(1) = llrecv1(1) .OR. ANY(lrecv_bdy(ib_bdy,igrd,1,:)) ! receive west145 llrecv1(3) = llrecv1(3) .OR. ANY(lrecv_bdy(ib_bdy,igrd,3,:)) ! receive south146 END DO147 148 IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN ! if need to send/recv in at least one direction149 CALL lbc_lnk( 'bdydyn2d', zhke, 'T', 1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )150 END IF151 END IF152 !153 112 CASE ( nkeg_HW ) !-- Hollingsworth scheme --! 154 113 DO jk = 1, jpkm1 … … 168 127 END DO 169 128 END DO 170 IF (ln_bdy) THEN171 ! Maria Luneva & Fred Wobus: July-2016172 ! compensate for lack of turbulent kinetic energy on liquid bdy points173 DO ib_bdy = 1, nb_bdy174 IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN175 igrd = 1 ! compensation null velocity on land at the bdy176 DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)177 ji = idx_bdy(ib_bdy)%nbi(jb,igrd) ! maximum extent : from 1 to jpi178 jj = idx_bdy(ib_bdy)%nbj(jb,igrd) ! maximum extent : from 1 to jpj179 IF( ji == 1 .OR. ji == jpi .OR. jj == 1 .OR. jj == jpj ) CYCLE180 DO jk = 1, jpkm1181 zhke(ji,jj,jk) = 0._wp182 zweightu = 8._wp * ( umask(ji-1,jj ,jk) + umask(ji ,jj ,jk) ) &183 & + 2._wp * ( umask(ji-1,jj-1,jk) + umask(ji-1,jj+1,jk) + umask(ji ,jj-1,jk) + umask(ji ,jj+1,jk) )184 zweightv = 8._wp * ( vmask(ji ,jj-1,jk) + vmask(ji ,jj-1,jk) ) &185 & + 2._wp * ( vmask(ji-1,jj-1,jk) + vmask(ji+1,jj-1,jk) + vmask(ji-1,jj ,jk) + vmask(ji+1,jj ,jk) )186 zu = 8._wp * ( un(ji-1,jj ,jk) * un(ji-1,jj ,jk) &187 & + un(ji ,jj ,jk) * un(ji ,jj ,jk) ) &188 & + ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) &189 & + ( un(ji ,jj-1,jk) + un(ji ,jj+1,jk) ) * ( un(ji ,jj-1,jk) + un(ji ,jj+1,jk) )190 zv = 8._wp * ( vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) &191 & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) &192 & + ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) &193 & + ( vn(ji-1,jj ,jk) + vn(ji+1,jj ,jk) ) * ( vn(ji-1,jj ,jk) + vn(ji+1,jj ,jk) )194 IF( zweightu > 0._wp ) zhke(ji,jj,jk) = zhke(ji,jj,jk) + zu / ( 2._wp * zweightu )195 IF( zweightv > 0._wp ) zhke(ji,jj,jk) = zhke(ji,jj,jk) + zv / ( 2._wp * zweightv )196 END DO197 END DO198 END IF199 END DO200 END IF201 129 CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. ) 202 130 ! -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/SBC/fldread.F90
r11223 r11258 465 465 ELSE ; it_offset = 0 466 466 ENDIF 467 IF( PRESENT(kt_offset) ) it_offset = kt_offset467 IF( PRESENT(kt_offset) ) it_offset = kt_offset 468 468 IF( PRESENT(kit) ) THEN ; it_offset = ( kit + it_offset ) * NINT( rdt/REAL(nn_baro,wp) ) 469 469 ELSE ; it_offset = it_offset * NINT( rdt ) … … 486 486 ! forcing record : 1 487 487 ! 488 ztmp = REAL( nsec_year, wp ) / ( REAL( nyear_len(1), wp ) * rday ) + 0.5 &489 &+ REAL( it_offset, wp ) / ( REAL( nyear_len(1), wp ) * rday )488 ztmp = REAL( nsec_year, wp ) / ( REAL( nyear_len(1), wp ) * rday ) + 0.5 & 489 & + REAL( it_offset, wp ) / ( REAL( nyear_len(1), wp ) * rday ) 490 490 sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/)) 491 491 ! swap at the middle of the year … … 517 517 ! forcing record : nmonth 518 518 ! 519 ztmp = REAL( nsec_month, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday ) + 0.5 &520 & + REAL(it_offset, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday )519 ztmp = REAL( nsec_month, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday ) + 0.5 & 520 & + REAL( it_offset, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday ) 521 521 imth = nmonth + INT( ztmp ) - COUNT((/llbefore/)) 522 522 IF( sdjf%cltype == 'monthly' ) THEN ; sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/))
Note: See TracChangeset
for help on using the changeset viewer.