Changeset 2351
- Timestamp:
- 2010-11-02T08:09:00+01:00 (13 years ago)
- Location:
- branches/nemo_v3_3_beta
- Files:
-
- 4 edited
-
DOC/TexFiles/Chapters/Chap_SBC.tex (modified) (3 diffs)
-
NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90 (modified) (3 diffs)
-
NEMOGCM/NEMO/OPA_SRC/IOM/iom_nf90.F90 (modified) (2 diffs)
-
NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_SBC.tex
r2349 r2351 641 641 if this is left as an empty string no action is taken. 642 642 In the model, weights files are read in and stored in a structured type (WGT) in the fldread module, as and when they are first required. 643 This initialisation procedure tries to determine whether the input data grid should be treated as cyclical or not. 643 This initialisation procedure determines whether the input data grid should be treated as cyclical or not by inspecting a global attribute stored in the weights input file. 644 This attribute must be called "ew_wrap" and be of integer type. 645 If it is negative, the input non-model grid is assumed not to be cyclic. 646 If zero or greater, then the value represents the number of columns that overlap. 647 E.g. if the input grid has columns at longitudes 0, 1, 2, .... , 359, then ew_wrap should be set to 0; 648 if longitudes are 0.5, 2.5, .... , 358.5, 360.5, 362.5, ew_wrap should be 2. 649 If the model does not find attribute ew_wrap, then a value of -999 is assumed. 650 In this case the fld_read routine defaults ew_wrap to value 0 and therefore the grid is assumed to be cyclic 651 with no overlapping columns. 644 652 (In fact this only matters when bicubic interpolation is required.) 645 To do this the model looks in the input data file (i.e. the data to which the weights are to be applied) for a variable with name "nav\_lon" or "lon". 646 If found, it checks the difference between the first and last values of longitude along a single row. 647 If the absolute value of this difference is close to 360 degrees or less than twice the maximum spacing from 360 degrees, the grid is assumed to be cyclical, and the difference determines whether the first column is a repeat of the last one or not. 648 If neither "nav\_lon" or "lon" can be found, the model resorts to looking at the first and last columns of data. 649 If the sum of the absolute values of the differences between the columns is very small, then the grid is assumed to be cyclical with coincident first and last columns. 650 If both of these tests fail, the grid is assumed not to be cyclical. 653 Note that no testing is done to check the validity in the model, since there is no way of knowing the name used for the longitude variable, 654 so it is up to the user to make sure his or her data is correctly represented. 651 655 652 656 Next the routine reads in the weights. … … 655 659 The size of the input data array is determined by examining the values in the "src" arrays to find the minimum and maximum i and j values required. 656 660 Since bicubic interpolation requires the calculation of gradients at each point on the grid, the corresponding arrays are dimensioned with a halo of width one grid point all the way around. 657 When the array of points from the data file is adjacent to an edge of the data grid, the halo is either a copy of the row/column next to it (non-cyclical case), or is a copy of one from the first two rows/columns on the opposite side of the grid (cyclical case with coincident end rows/columns, or cyclical case with non-coincident end rows/columns).661 When the array of points from the data file is adjacent to an edge of the data grid, the halo is either a copy of the row/column next to it (non-cyclical case), or is a copy of one from the first few columns on the opposite side of the grid (cyclical case). 658 662 659 663 \subsection{Limitations} … … 662 666 \begin{description} 663 667 \item 664 Input data grids must be logically rectangular.668 The case where input data grids are not logically rectangular has not been tested. 665 669 \item 666 670 This code is not guaranteed to produce positive definite answers from positive definite inputs. -
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90
r2287 r2351 44 44 #endif 45 45 PUBLIC iom_init, iom_swap, iom_open, iom_close, iom_setkt, iom_varid, iom_get, iom_gettime, iom_rstput, iom_put 46 PUBLIC iom_getatt 46 47 47 48 PRIVATE iom_rp0d, iom_rp1d, iom_rp2d, iom_rp3d … … 54 55 INTERFACE iom_get 55 56 MODULE PROCEDURE iom_g0d, iom_g1d, iom_g2d, iom_g3d 57 END INTERFACE 58 INTERFACE iom_getatt 59 MODULE PROCEDURE iom_g0d_intatt 56 60 END INTERFACE 57 61 INTERFACE iom_rstput … … 824 828 ! 825 829 END SUBROUTINE iom_gettime 830 831 832 !!---------------------------------------------------------------------- 833 !! INTERFACE iom_getatt 834 !!---------------------------------------------------------------------- 835 SUBROUTINE iom_g0d_intatt( kiomid, cdatt, pvar ) 836 INTEGER , INTENT(in ) :: kiomid ! Identifier of the file 837 CHARACTER(len=*), INTENT(in ) :: cdatt ! Name of the attribute 838 INTEGER , INTENT( out) :: pvar ! read field 839 ! 840 IF( kiomid > 0 ) THEN 841 IF( iom_file(kiomid)%nfid > 0 ) THEN 842 SELECT CASE (iom_file(kiomid)%iolib) 843 CASE (jpioipsl ) ; CALL ctl_stop('iom_getatt: only nf90 available') 844 CASE (jpnf90 ) ; CALL iom_nf90_getatt( kiomid, cdatt, pvar ) 845 CASE (jprstdimg) ; CALL ctl_stop('iom_getatt: only nf90 available') 846 CASE DEFAULT 847 CALL ctl_stop( 'iom_g0d_att: accepted IO library are only jpioipsl, jpnf90 and jprstdimg' ) 848 END SELECT 849 ENDIF 850 ENDIF 851 END SUBROUTINE iom_g0d_intatt 826 852 827 853 -
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/IOM/iom_nf90.F90
r2287 r2351 28 28 29 29 PUBLIC iom_nf90_open, iom_nf90_close, iom_nf90_varid, iom_nf90_get, iom_nf90_gettime, iom_nf90_rstput 30 PUBLIC iom_nf90_getatt 30 31 31 32 INTERFACE iom_nf90_get 32 33 MODULE PROCEDURE iom_nf90_g0d, iom_nf90_g123d 34 END INTERFACE 35 INTERFACE iom_nf90_getatt 36 MODULE PROCEDURE iom_nf90_intatt 33 37 END INTERFACE 34 38 INTERFACE iom_nf90_rstput … … 288 292 289 293 294 SUBROUTINE iom_nf90_intatt( kiomid, cdatt, pvar ) 295 !!----------------------------------------------------------------------- 296 !! *** ROUTINE iom_nf90_intatt *** 297 !! 298 !! ** Purpose : read an integer attribute with NF90 299 !!----------------------------------------------------------------------- 300 INTEGER , INTENT(in ) :: kiomid ! Identifier of the file 301 CHARACTER(len=*), INTENT(in ) :: cdatt ! attribute name 302 INTEGER , INTENT( out) :: pvar ! read field 303 ! 304 INTEGER :: if90id ! temporary integer 305 LOGICAL :: llok ! temporary logical 306 CHARACTER(LEN=100) :: clinfo ! info character 307 !--------------------------------------------------------------------- 308 ! 309 if90id = iom_file(kiomid)%nfid 310 llok = NF90_Inquire_attribute(if90id, NF90_GLOBAL, cdatt) == nf90_noerr 311 IF( llok) THEN 312 clinfo = 'iom_nf90_getatt, file: '//TRIM(iom_file(kiomid)%name)//', att: '//TRIM(cdatt) 313 CALL iom_nf90_check(NF90_GET_ATT(if90id, NF90_GLOBAL, cdatt, values=pvar), clinfo) 314 ELSE 315 CALL ctl_warn('iom_nf90_getatt: no attribute '//cdatt//' found') 316 pvar = -999 317 ENDIF 318 ! 319 END SUBROUTINE iom_nf90_intatt 320 321 290 322 SUBROUTINE iom_nf90_gettime( kiomid, kvid, ptime, cdunits, cdcalendar ) 291 323 !!-------------------------------------------------------------------- -
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.
