New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 2351 for branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90 – NEMO

Ignore:
Timestamp:
2010-11-02T08:09:00+01:00 (13 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.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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.