Ignore:
Timestamp:
08/30/07 14:44:23 (17 years ago)
Author:
smasson
Message:

bugfix for interpolation from ORCA2 without mask

Location:
trunk/SRC/Interpolation
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/SRC/Interpolation/compute_fromirr_bilinear_weigaddr.pro

    r257 r271  
    4040; @restrictions 
    4141;  -  the input grid must be an "irregular 2D grid", defined as a grid made 
    42 ;  of quadrilateral cells which corners positions are defined with olonin and olat 
     42;  of quadrilateral cells 
    4343;  -  We supposed the data are located on a sphere, with a periodicity along 
    4444;  the longitude 
    … …  
    7171  IF n_elements(omsk) EQ 1 AND omsk[0] EQ -1 THEN BEGIN 
    7272    omsk = replicate(1b, jpio, jpjo) 
    73    IF jpio EQ 182 AND jpjo EQ 149 THEN BEGIN 
    74       lontmp = (olonin[1:180, *] + 3600) MOD 360 
    75       lattmp = olat[1:180, *] 
    76       bad = abs(abs(lontmp - shift(lontmp, 1, 0)) - 180) LT 176 $ 
    77             OR  abs(abs(lontmp - shift(lontmp, 0, 1)) - 180) LT 176 $ 
    78             OR  abs(abs(lontmp - shift(lontmp, 1, 1)) - 180) LT 176 $ 
    79             OR  abs(abs(lattmp - shift(lattmp, 1, 0)) - 180) LT 176 $ 
    80             OR  abs(abs(lattmp - shift(lattmp, 0, 1)) - 180) LT 176 $ 
    81             OR  abs(abs(lattmp - shift(lattmp, 1, 1)) - 180) LT 176 
    82       bad = bad AND lattmp LT 50 
    83       omsk[1:180, 1:148] = 1b - bad[*, 1:148] 
    84       omsk[0, *] = omsk[180, *] 
    85       omsk[181, *] = omsk[1, *] 
    86     ENDIF ELSE omsk = replicate(1b, jpio, jpjo) 
    87   ENDIF 
    88   IF n_elements(amsk) EQ 1 AND amsk[0] EQ -1 THEN amsk = replicate(1b, jpia, jpja) 
    89   IF n_elements(omsk) NE jpio*jpjo THEN BEGIN 
    90     ras = report('input grid mask do not have the good size') 
    91     stop 
    92   ENDIF 
    93   IF n_elements(amsk) NE jpia*jpja THEN BEGIN 
    94     ras = report('output grid mask do not have the good size') 
    95     stop 
    96   ENDIF 
     73; if this is ORCA2 grid... 
     74    IF (jpio EQ 180 OR jpio EQ 182) AND (jpjo EQ 149  OR jpjo EQ 148  OR jpjo EQ 147 ) THEN BEGIN 
     75; we look for ill defined cells. 
     76      IF jpio EQ 182 THEN lontmp = olonin[1:180, *] ELSE lontmp = olonin 
     77; longitudinal size of the cells... 
     78      a = (lontmp-shift(lontmp, 1, 0)) 
     79      d1 = 0.5 * max( [[[a]], [[360+a]]] MOD 360, dimension = 3) 
     80      a = (lontmp-shift(lontmp, 1, 1)) 
     81      d2 = 0.5 * max( [[[a]], [[360+a]]] MOD 360, dimension = 3) 
     82      a = (shift(lontmp, 0, 1)-shift(lontmp, 1, 0)) 
     83      d3 = 0.5 * max( [[[a]], [[360+a]]] MOD 360, dimension = 3) 
     84      a = (shift(lontmp, 0, 1)-shift(lontmp, 1, 1)) 
     85      d4 = 0.5 * max( [[[a]], [[360+a]]] MOD 360, dimension = 3) 
     86      md = [[[d1]], [[d2]], [[d3]], [[d4]]] 
     87      bad = max(md, dimension = 3) GE 1.5 
     88      bad = (bad + shift(bad, -1, -1) + shift(bad, 0, -1) + shift(bad, -1, 0)) < 1 
     89      IF jpio EQ 182 THEN BEGIN 
     90        omsk[1:180, 80:120] = 1b - bad[*, 80:120] 
     91        omsk[0, *] = omsk[180, *] 
     92        omsk[181, *] = omsk[1, *] 
     93      ENDIF ELSE omsk[*, 80:120] = 1b - bad[*, 80:120] 
     94    ENDIF 
     95  ENDIF ELSE BEGIN  
     96    IF n_elements(omsk) NE jpio*jpjo THEN BEGIN 
     97      ras = report('input grid mask do not have the good size') 
     98      stop 
     99    ENDIF 
     100  ENDELSE  
     101  IF n_elements(amsk) EQ 1 AND amsk[0] EQ -1 THEN amsk = replicate(1b, jpia, jpja) ELSE BEGIN  
     102    IF n_elements(amsk) NE jpia*jpja THEN BEGIN 
     103      ras = report('output grid mask do not have the good size') 
     104      stop 
     105    ENDIF 
     106  ENDELSE  
    97107; 
    98108; longitude, between 0 and 360 
    … …  
    175185; ORCA cases : orca grid is irregular only northward of 40N 
    176186        CASE 1 OF 
    177           jpio EQ 92  AND (jpjo EQ 76   OR jpjo EQ 75   OR jpjo EQ 74  ):onsphere = yy GT 40 
    178           jpio EQ 180 AND (jpjo EQ 149  OR jpjo EQ 148  OR jpjo EQ 147 ):onsphere = yy GT 40 
    179           jpio EQ 720  AND (jpjo EQ 522  OR jpjo EQ 521  OR jpjo EQ 520 ):onsphere = yy GT 40 
    180           jpio EQ 1440 AND (jpjo EQ 1021 OR jpjo EQ 1020 OR jpjo EQ 1019):onsphere = yy GT 40 
     187          (jpio EQ 90 OR jpio EQ 92) AND (jpjo EQ 76   OR jpjo EQ 75   OR jpjo EQ 74  ):onsphere = yy GT 40 
     188          (jpio EQ 180 OR jpio EQ 182) AND (jpjo EQ 149  OR jpjo EQ 148  OR jpjo EQ 147 ):onsphere = yy GT 40 
     189          (jpio EQ 720 OR jpio EQ 722)  AND (jpjo EQ 522  OR jpjo EQ 521  OR jpjo EQ 520 ):onsphere = yy GT 40 
     190          (jpio EQ 1440 OR jpio EQ 1442) AND (jpjo EQ 1021 OR jpjo EQ 1020 OR jpjo EQ 1019):onsphere = yy GT 40 
    181191          ELSE:onsphere = 1 
    182192        ENDCASE 
    … …  
    233243          ELSE: BEGIN 
    234244; we keep its address 
    235         foraddr[n] = ind 
     245            foraddr[n] = ind 
    236246; keep the new coordinates 
    237247            forweight[n, 0] = xy[0] 
    … …  
    274284  a = reform(forweight[*, 0], 1, nawater) 
    275285  b = reform(forweight[*, 1], 1, nawater) 
    276   forweight =  -1 ; free memory 
     286  forweight =  -1               ; free memory 
    277287  newaweig = [(1-a)*(1-b), (1-b)*a, a*b, (1-a)*b] 
    278   a = -1 &  b = -1 ; free memory 
     288  a = -1 &  b = -1              ; free memory 
    279289; mask the weight to suppress the corner located on land 
    280290  newaweig = newaweig*((omsk)[newaaddr]) 

  • trunk/SRC/Interpolation/extrapolate.pro

    r262 r271  
    1616; put -1 if input data are not masked 
    1717; 
    18 ; @param nb_iteration {in}{optional}{type=integer scalar}{default=10.E20} 
     18; @param nb_iteration {in}{optional}{type=integer}{default=large enough to fill everything} 
    1919; Maximum number of iterations done in the extrapolation process. If there 
    2020; is no more masked values we exit extrapolate before reaching nb_iteration 
    21 ; (to be sure to fill everything, you can use a very large value) 
    2221; 
    2322; @keyword X_PERIODIC {type=scalar, 0 or 1}{default=0} 
    … …  
    3231; @keyword GE0 {type=scalar 0 or 1}{default=0} 
    3332; put 1 to force the extrapolated values to be larger than 0, same as using minval=0. 
     33; 
     34; @keyword  
     35; _EXTRA to be able to call extrapolate with _extra keyword 
    3436; 
    3537; @returns 
    … …  
    5456; 
    5557FUNCTION extrapolate, zinput, maskinput, nb_iteration, X_PERIODIC = x_periodic $ 
    56                       , MINVAL = minval, MAXVAL = maxval, GE0 = ge0 
     58                      , MINVAL = minval, MAXVAL = maxval, GE0 = ge0, _EXTRA = ex 
    5759; 
    5860  compile_opt idl2, strictarrsubs 
    5961; 
    6062; check the number of iteration used in the extrapolation. 
    61   IF n_elements(nb_iteration) EQ 0 THEN nb_iteration = 10.E20 
     63  szin = size(zinput) 
     64  IF szin[0] NE 2 THEN return, -1. ELSE szin = szin[1:2] 
     65  nx = szin[0] 
     66  ny = szin[1] 
     67  IF n_elements(nb_iteration) EQ 0 THEN nb_iteration = max(szin) 
    6268  IF nb_iteration EQ 0 THEN return, zinput 
    63   nx = (size(zinput))[1] 
    64   ny = (size(zinput))[2] 
    6569; take care of the boundary conditions... 
    6670; 

  • trunk/SRC/Interpolation/extrapsmooth.pro

    r262 r271  
    2727; put 1 to force the extrapolated values to be larger than 0, same as using minval=0. 
    2828; 
     29; @keyword  
     30; _EXTRA to be able to call extrapsmooth with _extra keyword 
     31; 
    2932; @returns 
    3033; the extrapolated array with no more masked values 
    … …  
    4750;- 
    4851; 
    49 FUNCTION extrapsmooth, in, mskin, X_PERIODIC = x_periodic, MINVAL = minval, MAXVAL = maxval, GE0 = ge0 
     52FUNCTION extrapsmooth, in, mskin, X_PERIODIC = x_periodic, MINVAL = minval, MAXVAL = maxval, GE0 = ge0, _EXTRA = ex 
    5053; 
    5154  compile_opt strictarr, strictarrsubs 

  • trunk/SRC/Interpolation/fromirr.pro

    r238 r271  
    4949; lonout, latout are used only to know the output domain size 
    5050; 
     51; @keyword  
     52; _EXTRA to be able to call fromirr with _extra keyword 
     53; 
    5154; @returns 
    5255; 2D array the interpolated data 
    … …  
    8790; 
    8891FUNCTION fromirr, method, datain, lonin, latin, mskin, lonout, latout, mskout $ 
    89                   , WEIG = weig, ADDR = addr 
     92                  , WEIG = weig, ADDR = addr, _EXTRA = ex 
    9093; 
    9194  compile_opt strictarr, strictarrsubs 

  • trunk/SRC/Interpolation/fromreg.pro

    r238 r271  
    5353; of the input data when performing the interpolation. 
    5454; 
     55; @keyword  
     56; _EXTRA to be able to call fromreg with _extra keyword 
     57; 
    5558; @returns 
    5659; 2D array the interpolated data 
    … …  
    9093                  , WEIG = weig, ADDR = addr $ 
    9194                  , NONORTHERNLINE = nonorthernline $ 
    92                   , NOSOUTHERNLINE = nosouthernline 
     95                  , NOSOUTHERNLINE = nosouthernline, _EXTRA = ex 
    9396; 
    9497  compile_opt idl2, strictarrsubs 

  • trunk/SRC/Interpolation/get_gridparams.pro

    r242 r271  
    22; 
    33; @file_comments 
    4 ; 1) extract from a NetCDF file the longitude, latitude, and their dimensions 
     4; Case 1: extract from a NetCDF file longitude and latitude arrays, their dimensions 
    55; and make sure it is 1D or 2D arrays 
    66; 
    7 ; or 
    8 ; 2) given longitude and latitude arrays, get their dimensions and make 
     7; Case 2: given longitude and latitude arrays, get their dimensions and make 
    98; sure they are 1D or 2D arrays 
    109; 
    … …  
    1413; @examples 
    1514; 
    16 ; 1) 
    17 ; IDL> get_gridparams, file, lonname, latname, lon, lat, jpi, jpj, n_dimensions 
    18 ; 
    19 ; or 
    20 ; 
    21 ; 2) 
     15; Case 1: 
     16; IDL> get_gridparams, file name/id, lonname, latname, lon, lat, jpi, jpj, n_dimensions 
     17; 
     18; Case 2: 
    2219; IDL> get_gridparams, lon, lat, jpi, jpj, n_dimensions 
    2320; 
    2421; @param in1 {in}{required} 
    25 ; Case 1: the name of the netcdf file 
    26 ; Case 2: 1d or 2d arrays defining longitudes and latitudes. 
    27 ; Out: the variable that will contain the longitudes 
     22; Case 1: name or the id (returned by ncdf_open) of the netcdf file 
     23; Case 2: 1D or 2D arrays defining longitudes, will be forced to have 
     24;         n_dimensions dimensions when returned 
    2825; 
    2926; @param in2 {in}{required} 
    30 ; Case 1: the name of the variable that contains the longitude in the NetCDF file 
    31 ; Case 2: 1d or 2d arrays defining longitudes and latitudes. 
    32 ;         Note that these arrays are also outputs and can therefore be modified. 
    33 ; Out: the variable that will contain the latitudes 
    34 ; 
    35 ; @param in3 {in}{required} 
    36 ; Case 1: the name of the variable that contains the latitude in the NetCDF file 
    37 ; Case 2: the number of points in the longitudinal direction. 
     27; Case 1: name of the variable containing the longitude in the NetCDF file 
     28; Case 2: 1D or 2D arrays defining latitudes, will be forced to have 
     29;         n_dimensions dimensions when returned 
     30; 
     31; @param in3  
     32; Case 1: name of the variable containing the latitude in the NetCDF file 
     33; Case 2: returned number of points in longitudinal direction. 
    3834; 
    3935; @param in4 {out} 
    40 ; Case 1: the number of points in the longitudinal direction 
    41 ; Case 2: the number of points in the latitudinal direction 
    42 ; 
    43 ; @param in5 {out} 
    44 ; Case 1: the number of points in the latitudinal direction 
    45 ; Case 2: 1 or 2 to specify if lon and lat should be 1D (jpi or jpj) 
    46 ; arrays or 2D arrays (jpi,jpj). Note that of  n_dimensions = 1, then the 
    47 ; grid must be regular (each longitude must be the same for all latitudes 
    48 ; and each latitude should be the same for all longitudes). 
     36; Case 1: returned longitudes array, with n_dimensions dimensions 
     37; Case 2: returned number of points in latitudinal direction 
     38; 
     39; @param in5 
     40; Case 1: returned latitudes array, with n_dimensions dimensions  
     41; Case 2: input scalar (1 or 2) to specify if lon and lat should be returned  
     42;         as 1D or 2D arrays. Note that if n_dimensions = 1, the grid must be 
     43;         regular (longitude and latitudes can be described as 1D arrays). 
    4944; 
    5045; @param in6 {out} 
    51 ; the variable that will contain the longitudes 
     46; Case 1: returned number of points in longitudinal direction. 
    5247; 
    5348; @param in7 {out} 
    54 ; the variable that will contain the latitudes 
    55 ; 
    56 ; @param in8 {out} 
    57 ; 1 or 2 to specify if lon and lat should be 1D (jpi or jpj) 
    58 ; 
    59 ; @keyword DOUBLE 
    60 ; use double precision to perform the computation 
     49; Case 1: returned number of points in latitudinal direction 
     50; 
     51; @param in8 {in} 
     52; Case 1: input scalar (1 or 2) to specify if lon and lat should be returned  
     53;         as 1D or 2D arrays. Note that if n_dimensions = 1, the grid must be 
     54;         regular (longitude and latitudes can be described as 1D arrays). 
     55; 
     56; @keyword DOUBLE {default=0}{type=scalar: 0 or 1} 
     57; activate to force double precision for lon/lat arrays 
    6158; 
    6259; @examples 
    … …  
    8481    8:BEGIN 
    8582; get longitude and latitude 
    86       IF file_test(in1) EQ 0 THEN BEGIN 
    87         ras = report('file ' + in1 + ' does not exist') 
    88         stop 
    89       ENDIF 
    90       cdfido = ncdf_open(in1) 
     83      IF size(in1, /type) EQ 7 THEN BEGIN 
     84        IF file_test(in1) EQ 0 THEN BEGIN 
     85          ras = report('file ' + in1 + ' does not exist') 
     86          stop 
     87        ENDIF 
     88        cdfido = ncdf_open(in1) 
     89      ENDIF ELSE cdfido = in1 
    9190      ncdf_varget, cdfido, in2, lon 
    9291      ncdf_varget, cdfido, in3, lat 
    93       ncdf_close, cdfido 
     92      IF size(in1, /type) EQ 7 THEN ncdf_close, cdfido 
    9493      n_dimensions = in8 
    9594    END 
Note: See TracChangeset for help on using the changeset viewer.