Changeset 2351


Ignore:
Timestamp:
2010-11-02T08:09:00+01:00 (10 years ago)
Author:
sga
Message:

NEMO branch nemo_v3_3_beta
modify interpolation on the fly scheme within fldread.F90.
Instead of attempting to decide on cyclicity of input non-model grid by inspecting a longitude variable,
fld_read now expects to find an attribute in the input weights file to tell it how many columns overlap
at the east-west edges of the grid.
iom and iom_nf90.F90 have new routines to expose netcdf attributes to the model.

Location:
branches/nemo_v3_3_beta
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_SBC.tex

    r2349 r2351  
    641641if this is left as an empty string no action is taken. 
    642642In 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. 
     643This 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. 
     644This attribute must be called "ew_wrap" and be of integer type. 
     645If it is negative, the input non-model grid is assumed not to be cyclic. 
     646If zero or greater, then the value represents the number of columns that overlap. 
     647E.g. if the input grid has columns at longitudes 0, 1, 2, .... , 359, then ew_wrap should be set to 0; 
     648if longitudes are 0.5, 2.5, .... , 358.5, 360.5, 362.5, ew_wrap should be 2. 
     649If the model does not find attribute ew_wrap, then a value of -999 is assumed. 
     650In this case the fld_read routine defaults ew_wrap to value 0 and therefore the grid is assumed to be cyclic 
     651with no overlapping columns. 
    644652(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. 
     653Note 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, 
     654so it is up to the user to make sure his or her data is correctly represented. 
    651655 
    652656Next the routine reads in the weights. 
     
    655659The 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. 
    656660Since 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). 
     661When 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). 
    658662 
    659663\subsection{Limitations} 
     
    662666\begin{description} 
    663667\item 
    664 Input data grids must be logically rectangular. 
     668The case where input data grids are not logically rectangular has not been tested. 
    665669\item 
    666670This 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  
    4444#endif 
    4545   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 
    4647 
    4748   PRIVATE iom_rp0d, iom_rp1d, iom_rp2d, iom_rp3d 
     
    5455   INTERFACE iom_get 
    5556      MODULE PROCEDURE iom_g0d, iom_g1d, iom_g2d, iom_g3d 
     57   END INTERFACE 
     58   INTERFACE iom_getatt 
     59      MODULE PROCEDURE iom_g0d_intatt 
    5660   END INTERFACE 
    5761   INTERFACE iom_rstput 
     
    824828      ! 
    825829   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 
    826852 
    827853 
  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/IOM/iom_nf90.F90

    r2287 r2351  
    2828 
    2929   PUBLIC iom_nf90_open, iom_nf90_close, iom_nf90_varid, iom_nf90_get, iom_nf90_gettime, iom_nf90_rstput 
     30   PUBLIC iom_nf90_getatt 
    3031 
    3132   INTERFACE iom_nf90_get 
    3233      MODULE PROCEDURE iom_nf90_g0d, iom_nf90_g123d 
     34   END INTERFACE 
     35   INTERFACE iom_nf90_getatt 
     36      MODULE PROCEDURE iom_nf90_intatt 
    3337   END INTERFACE 
    3438   INTERFACE iom_nf90_rstput 
     
    288292 
    289293 
     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 
    290322   SUBROUTINE iom_nf90_gettime( kiomid, kvid, ptime, cdunits, cdcalendar ) 
    291323      !!-------------------------------------------------------------------- 
  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90

    r2330 r2351  
    7272      INTEGER                                 ::   numwgt       ! number of weights (4=bilinear, 16=bicubic) 
    7373      INTEGER                                 ::   nestid       ! for agrif, keep track of nest we're in 
    74       INTEGER                                 ::   offset       ! =0 when cyclic grid has coincident first/last columns,  
    75                                                                 ! =1 when they assumed to be one grid spacing apart       
    76                                                                 ! =-1 otherwise 
     74      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 
    7777      LOGICAL                                 ::   cyclic       ! east-west cyclic or not 
    7878      INTEGER, DIMENSION(:,:,:), POINTER      ::   data_jpi     ! array of source integers 
    7979      INTEGER, DIMENSION(:,:,:), POINTER      ::   data_jpj     ! array of source integers 
    8080      REAL(wp), DIMENSION(:,:,:), POINTER     ::   data_wgt     ! array of weights on model grid 
    81       REAL(wp), DIMENSION(:,:,:), POINTER       ::   fly_dta      ! array of values on input grid 
    82       REAL(wp), DIMENSION(:,:,:), POINTER       ::   col2         ! temporary array for reading in columns 
     81      REAL(wp), DIMENSION(:,:,:), POINTER     ::   fly_dta      ! array of values on input grid 
     82      REAL(wp), DIMENSION(:,:,:), POINTER     ::   col          ! temporary array for reading in columns 
    8383   END TYPE WGT 
    8484 
     
    807807         IF( ref_wgts(kw)%cyclic ) THEN 
    808808            WRITE(numout,*) '       cyclical' 
    809             IF( ref_wgts(kw)%offset > 0 ) WRITE(numout,*) '                 with offset' 
     809            IF( ref_wgts(kw)%overlap > 0 ) WRITE(numout,*) '              with overlap of ', ref_wgts(kw)%overlap 
    810810         ELSE 
    811811            WRITE(numout,*) '       not cyclical' 
     
    835835      INTEGER , DIMENSION(jpi, jpj)           ::   data_src 
    836836      REAL(wp), DIMENSION(jpi, jpj)           ::   data_tmp 
    837       REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::   line2, lines  ! temporary array to read 2 lineumns 
    838       CHARACTER (len=34)                      ::   lonvar 
    839837      LOGICAL                                 ::   cyclical 
    840       REAL(wp)                                ::   resid, dlon   ! temporary array to read 2 lineumns 
    841       INTEGER                                 ::   offset        ! temporary integer 
     838      INTEGER                                 ::   zwrap         ! temporary integer 
     839      INTEGER                                 ::   overlap        ! temporary integer 
    842840      !!---------------------------------------------------------------------- 
    843841      ! 
     
    857855      id = iom_varid( inum, sd%clvar, ddims ) 
    858856 
    859       !! check for an east-west cyclic grid 
    860       !! try to guess name of longitude variable 
    861  
    862       lonvar = 'nav_lon' 
    863       id = iom_varid(inum, TRIM(lonvar), ldstop=.FALSE.) 
    864       IF( id <= 0 ) THEN 
    865          lonvar = 'lon' 
    866          id = iom_varid(inum, TRIM(lonvar), ldstop=.FALSE.) 
    867       ENDIF 
    868  
    869       offset = -1 
    870       cyclical = .FALSE. 
    871       IF( id > 0 ) THEN 
    872          !! found a longitude variable 
    873          !! now going to assume that grid is regular so we can read a single row 
    874  
    875          !! because input array is 2d, have to present iom with 2d array even though we only need 1d slice 
    876          !! worse, we cant pass line2(:,1) to iom_get since this is treated as a 1d array which doesnt match input file 
    877          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 spacing 
    881          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 ) THEN 
    886             !! end rows overlap in longitude 
    887             offset = 0 
    888             cyclical = .TRUE. 
    889          ELSEIF( resid < 2.0*dlon ) THEN 
    890             !! also call it cyclic if difference between end points is less than twice dlon from 360 
    891             offset = 1 
    892             cyclical = .TRUE. 
    893          ENDIF 
    894  
    895          DEALLOCATE( lines ) 
    896  
    897       ELSE 
    898          !! guessing failed 
    899          !! read in first and last columns of data variable 
    900          !! 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-west 
    902  
    903          !! because input array is 2d, have to present iom with 2d array even though we only need 1d slice 
    904          !! worse, we cant pass line2(1,:) to iom_get since this is treated as a 1d array which doesnt match input file 
    905          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 ) THEN 
    914             offset = 0 
    915             cyclical = .TRUE. 
    916          ENDIF 
    917  
    918          DEALLOCATE( lines, line2 ) 
    919       ENDIF 
    920  
    921857      !! close it 
    922858      CALL iom_close( inum ) 
     
    926862      CALL iom_open ( sd%wgtname, inum )   ! interpolation weights 
    927863      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 
    928879 
    929880         ref_wgts(nxt_wgt)%ddims(1) = ddims(1) 
    930881         ref_wgts(nxt_wgt)%ddims(2) = ddims(2) 
    931882         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 
    938885         ref_wgts(nxt_wgt)%nestid = 0 
    939886#if defined key_agrif 
     
    995942         ipk =  SIZE(sd%fnow,3) 
    996943         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)%col2(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) ) 
    998945 
    999946         nxt_wgt = nxt_wgt + 1 
     
    1035982      !! so we need to have a 4 by 4 subgrid surrounding each model point to cover both cases 
    1036983 
    1037       !! sub grid where we already have weights 
     984      !! sub grid from non-model input grid which encloses all grid points in this nemo process 
    1038985      jpimin = ref_wgts(kw)%botleft(1) 
    1039986      jpjmin = ref_wgts(kw)%botleft(2) 
     
    1041988      jpjwid = ref_wgts(kw)%jpjwgt 
    1042989 
    1043       !! what we need to read into sub grid in order to calculate gradients 
     990      !! when reading in, expand this sub-grid by one halo point all the way round for calculating gradients 
    1044991      rec1(1) = MAX( jpimin-1, 1 ) 
    1045992      rec1(2) = MAX( jpjmin-1, 1 ) 
     
    1049996      recn(3) = kk 
    1050997 
    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 
    10521001      jpi1 = 2 + rec1(1) - jpimin 
    10531002      jpj1 = 2 + rec1(2) - jpjmin 
     
    10641013 
    10651014      !! first four weights common to both bilinear and bicubic 
     1015      !! data_jpi, data_jpj have already been shifted to (1,1) corresponding to botleft 
    10661016      !! note that we have to offset by 1 into fly_dta array because of halo 
    10671017      dta(:,:,:) = 0.0 
    10681018      DO jk = 1,4 
    1069         DO jn = 1, nlcj 
    1070           DO jm = 1,nlci 
     1019        DO jn = 1, jpj 
     1020          DO jm = 1,jpi 
    10711021            ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 
    10721022            nj = ref_wgts(kw)%data_jpj(jm,jn,jk) 
     
    10961046        IF( ref_wgts(kw)%cyclic ) THEN 
    10971047           rec1(2) = MAX( jpjmin-1, 1 ) 
    1098            recn(1) = 2 
     1048           recn(1) = 1 
    10991049           recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 ) 
    11001050           jpj1 = 2 + rec1(2) - jpjmin 
    11011051           jpj2 = jpj1 + recn(2) - 1 
    11021052           IF( jpi1 == 2 ) THEN 
    1103               rec1(1) = ref_wgts(kw)%ddims(1) - 1 
    1104               SELECT CASE( SIZE( ref_wgts(kw)%col2(:,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) ) 
    11051055              CASE(1) 
    1106                    CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2,1), nrec, rec1, recn) 
     1056                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn) 
    11071057              CASE(jpk)          
    1108                    CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2,:), nrec, rec1, recn) 
     1058                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn) 
    11091059              END SELECT       
    1110               ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col2(ref_wgts(kw)%offset+1,jpj1:jpj2,:) 
     1060              ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:) 
    11111061           ENDIF 
    11121062           IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN 
    1113               rec1(1) = 1 
    1114               SELECT CASE( SIZE( ref_wgts(kw)%col2(:,jpj1:jpj2,:),3) ) 
     1063              rec1(1) = 1 + ref_wgts(kw)%overlap 
     1064              SELECT CASE( SIZE( ref_wgts(kw)%col(:,jpj1:jpj2,:),3) ) 
    11151065              CASE(1) 
    1116                    CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2,1), nrec, rec1, recn) 
     1066                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,1), nrec, rec1, recn) 
    11171067              CASE(jpk) 
    1118                    CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2,:), nrec, rec1, recn) 
     1068                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn) 
    11191069              END SELECT 
    1120               ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col2(2-ref_wgts(kw)%offset,jpj1:jpj2,:) 
     1070              ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:) 
    11211071           ENDIF 
    11221072        ENDIF 
     
    11241074        ! gradient in the i direction 
    11251075        DO jk = 1,4 
    1126           DO jn = 1, nlcj 
    1127             DO jm = 1,nlci 
     1076          DO jn = 1, jpj 
     1077            DO jm = 1,jpi 
    11281078              ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 
    11291079              nj = ref_wgts(kw)%data_jpj(jm,jn,jk) 
     
    11361086        ! gradient in the j direction 
    11371087        DO jk = 1,4 
    1138           DO jn = 1, nlcj 
    1139             DO jm = 1,nlci 
     1088          DO jn = 1, jpj 
     1089            DO jm = 1,jpi 
    11401090              ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 
    11411091              nj = ref_wgts(kw)%data_jpj(jm,jn,jk) 
Note: See TracChangeset for help on using the changeset viewer.