- Timestamp:
- 2010-11-02T08:09:00+01:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
r2330 r2351 72 72 INTEGER :: numwgt ! number of weights (4=bilinear, 16=bicubic) 73 73 INTEGER :: nestid ! for agrif, keep track of nest we're in 74 INTEGER :: o ffset ! =0 when cyclic grid has coincident first/last columns,75 ! = 1 when they assumed to be one grid spacing apart76 ! =-1 otherwise74 INTEGER :: overlap ! =0 when cyclic grid has no overlapping EW columns 75 ! =>1 when they have one or more overlapping columns 76 ! =-1 not cyclic 77 77 LOGICAL :: cyclic ! east-west cyclic or not 78 78 INTEGER, DIMENSION(:,:,:), POINTER :: data_jpi ! array of source integers 79 79 INTEGER, DIMENSION(:,:,:), POINTER :: data_jpj ! array of source integers 80 80 REAL(wp), DIMENSION(:,:,:), POINTER :: data_wgt ! array of weights on model grid 81 REAL(wp), DIMENSION(:,:,:), POINTER 82 REAL(wp), DIMENSION(:,:,:), POINTER :: col2! temporary array for reading in columns81 REAL(wp), DIMENSION(:,:,:), POINTER :: fly_dta ! array of values on input grid 82 REAL(wp), DIMENSION(:,:,:), POINTER :: col ! temporary array for reading in columns 83 83 END TYPE WGT 84 84 … … 807 807 IF( ref_wgts(kw)%cyclic ) THEN 808 808 WRITE(numout,*) ' cyclical' 809 IF( ref_wgts(kw)%o ffset > 0 ) WRITE(numout,*) ' with offset'809 IF( ref_wgts(kw)%overlap > 0 ) WRITE(numout,*) ' with overlap of ', ref_wgts(kw)%overlap 810 810 ELSE 811 811 WRITE(numout,*) ' not cyclical' … … 835 835 INTEGER , DIMENSION(jpi, jpj) :: data_src 836 836 REAL(wp), DIMENSION(jpi, jpj) :: data_tmp 837 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: line2, lines ! temporary array to read 2 lineumns838 CHARACTER (len=34) :: lonvar839 837 LOGICAL :: cyclical 840 REAL(wp) :: resid, dlon ! temporary array to read 2 lineumns841 INTEGER :: o ffset! temporary integer838 INTEGER :: zwrap ! temporary integer 839 INTEGER :: overlap ! temporary integer 842 840 !!---------------------------------------------------------------------- 843 841 ! … … 857 855 id = iom_varid( inum, sd%clvar, ddims ) 858 856 859 !! check for an east-west cyclic grid860 !! try to guess name of longitude variable861 862 lonvar = 'nav_lon'863 id = iom_varid(inum, TRIM(lonvar), ldstop=.FALSE.)864 IF( id <= 0 ) THEN865 lonvar = 'lon'866 id = iom_varid(inum, TRIM(lonvar), ldstop=.FALSE.)867 ENDIF868 869 offset = -1870 cyclical = .FALSE.871 IF( id > 0 ) THEN872 !! found a longitude variable873 !! now going to assume that grid is regular so we can read a single row874 875 !! because input array is 2d, have to present iom with 2d array even though we only need 1d slice876 !! worse, we cant pass line2(:,1) to iom_get since this is treated as a 1d array which doesnt match input file877 ALLOCATE( lines(ddims(1),2) )878 CALL iom_get(inum, jpdom_unknown, lonvar, lines(:,:), 1, kstart=(/1,1/), kcount=(/ddims(1),2/) )879 880 !! find largest grid spacing881 lines(1:ddims(1)-1,2) = lines(2:ddims(1),1) - lines(1:ddims(1)-1,1)882 dlon = MAXVAL( lines(1:ddims(1)-1,2) )883 884 resid = ABS(ABS(lines(ddims(1),1)-lines(1,1))-360.0)885 IF( resid < rsmall ) THEN886 !! end rows overlap in longitude887 offset = 0888 cyclical = .TRUE.889 ELSEIF( resid < 2.0*dlon ) THEN890 !! also call it cyclic if difference between end points is less than twice dlon from 360891 offset = 1892 cyclical = .TRUE.893 ENDIF894 895 DEALLOCATE( lines )896 897 ELSE898 !! guessing failed899 !! read in first and last columns of data variable900 !! since we dont know the name of the longitude variable (or even if there is one)901 !! we assume that if these two columns are equal, file is cyclic east-west902 903 !! because input array is 2d, have to present iom with 2d array even though we only need 1d slice904 !! worse, we cant pass line2(1,:) to iom_get since this is treated as a 1d array which doesnt match input file905 ALLOCATE( lines(2,ddims(2)), line2(2,ddims(2)) )906 CALL iom_get(inum, jpdom_unknown, sd%clvar, line2(:,:), 1, kstart=(/1,1/), kcount=(/2,ddims(2)/) )907 lines(2,:) = line2(1,:)908 909 CALL iom_get(inum, jpdom_unknown, sd%clvar, line2(:,:), 1, kstart=(/ddims(1)-1,1/), kcount=(/2,ddims(2)/) )910 lines(1,:) = line2(2,:)911 912 resid = SUM( ABS(lines(1,:) - lines(2,:)) )913 IF( resid < ddims(2)*rsmall ) THEN914 offset = 0915 cyclical = .TRUE.916 ENDIF917 918 DEALLOCATE( lines, line2 )919 ENDIF920 921 857 !! close it 922 858 CALL iom_close( inum ) … … 926 862 CALL iom_open ( sd%wgtname, inum ) ! interpolation weights 927 863 IF ( inum > 0 ) THEN 864 865 !! determine whether we have an east-west cyclic grid 866 !! from global attribute called "ew_wrap" in the weights file 867 !! note that if not found, iom_getatt returns -999 and cyclic with no overlap is assumed 868 !! since this is the most common forcing configuration 869 870 CALL iom_getatt(inum, 'ew_wrap', zwrap) 871 IF( zwrap >= 0 ) THEN 872 cyclical = .TRUE. 873 ELSE IF( zwrap == -999 ) THEN 874 cyclical = .TRUE. 875 zwrap = 0 876 ELSE 877 cyclical = .FALSE. 878 ENDIF 928 879 929 880 ref_wgts(nxt_wgt)%ddims(1) = ddims(1) 930 881 ref_wgts(nxt_wgt)%ddims(2) = ddims(2) 931 882 ref_wgts(nxt_wgt)%wgtname = sd%wgtname 932 ref_wgts(nxt_wgt)%offset = -1 933 ref_wgts(nxt_wgt)%cyclic = .FALSE. 934 IF( cyclical ) THEN 935 ref_wgts(nxt_wgt)%offset = offset 936 ref_wgts(nxt_wgt)%cyclic = .TRUE. 937 ENDIF 883 ref_wgts(nxt_wgt)%overlap = zwrap 884 ref_wgts(nxt_wgt)%cyclic = cyclical 938 885 ref_wgts(nxt_wgt)%nestid = 0 939 886 #if defined key_agrif … … 995 942 ipk = SIZE(sd%fnow,3) 996 943 ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3 ,ipk) ) 997 IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col 2(2,ref_wgts(nxt_wgt)%jpjwgt+3,ipk) )944 IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col(1,ref_wgts(nxt_wgt)%jpjwgt+3,ipk) ) 998 945 999 946 nxt_wgt = nxt_wgt + 1 … … 1035 982 !! so we need to have a 4 by 4 subgrid surrounding each model point to cover both cases 1036 983 1037 !! sub grid where we already have weights984 !! sub grid from non-model input grid which encloses all grid points in this nemo process 1038 985 jpimin = ref_wgts(kw)%botleft(1) 1039 986 jpjmin = ref_wgts(kw)%botleft(2) … … 1041 988 jpjwid = ref_wgts(kw)%jpjwgt 1042 989 1043 !! wh at we need to read into sub grid in order to calculategradients990 !! when reading in, expand this sub-grid by one halo point all the way round for calculating gradients 1044 991 rec1(1) = MAX( jpimin-1, 1 ) 1045 992 rec1(2) = MAX( jpjmin-1, 1 ) … … 1049 996 recn(3) = kk 1050 997 1051 !! where we need to read it to 998 !! where we need to put it in the non-nemo grid fly_dta 999 !! note that jpi1 and jpj1 only differ from 1 when jpimin and jpjmin are 1 1000 !! (ie at the extreme west or south of the whole input grid) and similarly for jpi2 and jpj2 1052 1001 jpi1 = 2 + rec1(1) - jpimin 1053 1002 jpj1 = 2 + rec1(2) - jpjmin … … 1064 1013 1065 1014 !! first four weights common to both bilinear and bicubic 1015 !! data_jpi, data_jpj have already been shifted to (1,1) corresponding to botleft 1066 1016 !! note that we have to offset by 1 into fly_dta array because of halo 1067 1017 dta(:,:,:) = 0.0 1068 1018 DO jk = 1,4 1069 DO jn = 1, nlcj1070 DO jm = 1, nlci1019 DO jn = 1, jpj 1020 DO jm = 1,jpi 1071 1021 ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 1072 1022 nj = ref_wgts(kw)%data_jpj(jm,jn,jk) … … 1096 1046 IF( ref_wgts(kw)%cyclic ) THEN 1097 1047 rec1(2) = MAX( jpjmin-1, 1 ) 1098 recn(1) = 21048 recn(1) = 1 1099 1049 recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 ) 1100 1050 jpj1 = 2 + rec1(2) - jpjmin 1101 1051 jpj2 = jpj1 + recn(2) - 1 1102 1052 IF( jpi1 == 2 ) THEN 1103 rec1(1) = ref_wgts(kw)%ddims(1) - 11104 SELECT CASE( SIZE( ref_wgts(kw)%col 2(:,jpj1:jpj2,:),3) )1053 rec1(1) = ref_wgts(kw)%ddims(1) - ref_wgts(kw)%overlap 1054 SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) ) 1105 1055 CASE(1) 1106 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col 2(:,jpj1:jpj2,1), nrec, rec1, recn)1056 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn) 1107 1057 CASE(jpk) 1108 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col 2(:,jpj1:jpj2,:), nrec, rec1, recn)1058 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn) 1109 1059 END SELECT 1110 ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col 2(ref_wgts(kw)%offset+1,jpj1:jpj2,:)1060 ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:) 1111 1061 ENDIF 1112 1062 IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN 1113 rec1(1) = 1 1114 SELECT CASE( SIZE( ref_wgts(kw)%col 2(:,jpj1:jpj2,:),3) )1063 rec1(1) = 1 + ref_wgts(kw)%overlap 1064 SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) ) 1115 1065 CASE(1) 1116 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col 2(:,jpj1:jpj2,1), nrec, rec1, recn)1066 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn) 1117 1067 CASE(jpk) 1118 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col 2(:,jpj1:jpj2,:), nrec, rec1, recn)1068 CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn) 1119 1069 END SELECT 1120 ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col 2(2-ref_wgts(kw)%offset,jpj1:jpj2,:)1070 ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:) 1121 1071 ENDIF 1122 1072 ENDIF … … 1124 1074 ! gradient in the i direction 1125 1075 DO jk = 1,4 1126 DO jn = 1, nlcj1127 DO jm = 1, nlci1076 DO jn = 1, jpj 1077 DO jm = 1,jpi 1128 1078 ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 1129 1079 nj = ref_wgts(kw)%data_jpj(jm,jn,jk) … … 1136 1086 ! gradient in the j direction 1137 1087 DO jk = 1,4 1138 DO jn = 1, nlcj1139 DO jm = 1, nlci1088 DO jn = 1, jpj 1089 DO jm = 1,jpi 1140 1090 ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 1141 1091 nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
Note: See TracChangeset
for help on using the changeset viewer.