Changeset 167 for trunk/SRC/Computation/grad.pro
- Timestamp:
- 09/05/06 14:24:07 (18 years ago)
- Location:
- trunk/SRC/Computation
- Files:
-
- 1 added
- 1 moved
-
. (added)
-
grad.pro (moved) (moved from trunk/SRC/ToBeReviewed/CALCULS/grad.pro) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/SRC/Computation/grad.pro
r164 r167 4 4 ;+ 5 5 ; @file_comments 6 ; 6 ; compute the gradient of a variable 7 7 ; 8 8 ; @categories 9 ; 9 ; Calculation 10 10 ; 11 11 ; @param FIELD 12 ; 13 ; 14 ; @param DIREC 15 ; 16 ; 17 ; @returns 18 ; 12 ; The field for which we want to compute the gradient. A 2D (xy), 13 ; 3D (xyz or yt) or 4D (xyzt) array or a structure readable by litchamp 14 ; and containing a 2D (xy), 3D (xyz or yt) or 4D (xyzt) array. 15 ; note that the dimension of the arry must suit the domain dimension. 16 ; 17 ; @param DIREC {type=scalar string} 18 ; the gradient direction: 'x', 'y', 'z' 19 ; 20 ; @returns RES {type=2D, 3D or 4D array} 21 ; the gradient of the input data (with the same size) 19 22 ; 20 23 ; @uses 21 ; 24 ; cm_4cal, cm_4data, cm_4mmesh 22 25 ; 23 26 ; @restrictions 24 ; 27 ; - Works only for Arakawa C-grid. 28 ; - When computing the gradient, the result is not on the same grid point 29 ; than the input data. In consequence, we update, vargrid and the grid position 30 ; parameters (firstx[tuvf], lastx[tuvf], nx[tuvf], firsty[tuvf], lasty[tuvf], 31 ; ny[tuvf], firstz[tw], lastz[tw], nz[tw]). 32 ; - points that cannot be computed (domain bondaries, coastline) are set to NaN 33 ; - the common variable jpt is used to differ xyz (jpt=1) and xyt (jpt\=1) cases. 25 34 ; 26 35 ; @examples 27 ; 36 ; IDL> @tst_initorca2 37 ; IDL> plt, grad({arr:gphit,g:'T'}, 'x') 38 ; IDL> plt, grad({arr:gphit,g:'T'}, 'y') 28 39 ; 29 40 ; @history … … 33 44 ; $Id$ 34 45 ; 35 ; @todo36 ; seb: remplir les truc!37 46 ;- 38 47 ;------------------------------------------------------------ … … 43 52 compile_opt idl2, strictarrsubs 44 53 ; 45 @common 46 ;------------------------------------------------------------ 47 ; 48 IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $ 54 @cm_4cal ; for jpt 55 @cm_4data ; for varname, vargrid, vardate, varunit, valmask 56 @cm_4mesh 57 ;------------------------------------------------------------ 58 ; 59 IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $ 49 60 return, report(['This version of grad is based on Arakawa C-grid.' $ 50 61 , 'U and V grids must therefore be defined']) 51 62 ; 52 res = litchamp(field) 53 taille = size(res) 54 grille, mask, glam, gphi, gdep, nx, ny, nz $ 55 , firstx, firsty, firstz, lastx, lasty, lastz 56 if n_elements(valmask) EQ 0 then valmask = 1e20 57 case strupcase(vargrid) of 58 'T':BEGIN 59 case direc of 60 'x':BEGIN 61 divi = e1u[firstx:lastx, firsty:lasty] 62 newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz] 63 vargrid = 'U' 64 firstxu = firstxt & lastxu = lastxt & nxu = nxt 65 firstyu = firstyt & lastyu = lastyt & nyu = nyt 66 END 67 'y':BEGIN 68 divi = e2v[firstx:lastx, firsty:lasty] 69 newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz] 70 vargrid = 'V' 71 firstxv = firstxt & lastxv = lastxt & nxv = nxt 72 firstyv = firstyt & lastyv = lastyt & nyv = nyt 73 END 74 'z':BEGIN 75 divi = e3w[firstz:lastz] 76 newmask = mask 77 vargrid = 'W' 78 firstzw = firstzt & lastzw = lastzt & nzw = nzt 79 END 80 ELSE:return, report('Bad definition of direction argument') 81 ENDCASE 82 END 83 'W':BEGIN 84 case direc of 85 'x':BEGIN 86 divi = e1u[firstx:lastx, firsty:lasty] 87 newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz] 88 vargrid = 'U' 89 firstxu = firstxt & lastxu = lastxt & nxu = nxt 90 firstyu = firstyt & lastyu = lastyt & nyu = nyt 91 END 92 'y':BEGIN 93 divi = e2v[firstx:lastx, firsty:lasty] 94 newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz] 95 vargrid = 'V' 96 firstxv = firstxt & lastxv = lastxt & nxv = nxt 97 firstyv = firstyt & lastyv = lastyt & nyv = nyt 98 END 99 'z':BEGIN 100 divi = e3t[firstz:lastz] 101 newmask = mask 102 vargrid = 'T' 103 firstzt = firstzw & lastzt = lastzw & nzt = nzw 104 END 105 ELSE:return, report('Bad definition of direction argument') 106 endcase 107 END 108 'U':BEGIN 109 case direc of 110 'x':BEGIN 111 divi = (shift(e1t, -1, 0))[firstx:lastx, firsty:lasty] 112 newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz] 113 vargrid = 'T' 114 firstxt = firstxu & lastxt = lastxu & nxt = nxu 115 firstyt = firstyu & lastyt = lastyu & nyt = nyu 116 END 117 'y':BEGIN 118 divi = e2f[firstx:lastx, firsty:lasty] 119 newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz] 120 vargrid = 'F' 121 firstxf = firstxu & lastxf = lastxu & nxf = nxu 122 firstyf = firstyu & lastyf = lastyu & nyf = nyu 123 END 124 'z':BEGIN 125 divi = e3w[firstz:lastz] 126 newmask = mask 127 vargrid = 'W' 128 firstzw = firstzt & lastzw = lastzt & nzw = nzt 129 END 130 ELSE:return, report('Bad definition of direction argument') 131 endcase 132 END 133 'V':BEGIN 134 case direc of 135 'x':BEGIN 136 divi = e1f[firstx:lastx, firsty:lasty] 137 newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz] 138 vargrid = 'F' 139 firstxf = firstxv & lastxf = lastxv & nxf = nxv 140 firstyf = firstyv & lastyf = lastyv & nyf = nyv 141 END 142 'y':BEGIN 143 divi = (shift(e2t, 0, -1))[firstx:lastx, firsty:lasty] 144 newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz] 145 vargrid = 'T' 146 firstxt = firstxv & lastxt = lastxv & nxt = nxv 147 firstyt = firstyv & lastyt = lastyv & nyt = nyv 148 END 149 'z':BEGIN 150 divi = e3w[firstz:lastz] 151 newmask = mask 152 vargrid = 'W' 153 firstzw = firstzt & lastzw = lastzt & nzw = nzt 154 END 155 ELSE:return, report('Bad definition of direction argument') 156 endcase 157 END 158 ; 'F':BEGIN 63 res = litchamp(field) 64 szres = size(res) 65 grille, mask, glam, gphi, gdep, nx, ny, nz $ 66 , firstx, firsty, firstz, lastx, lasty, lastz 67 ; 68 if n_elements(valmask) EQ 0 then valmask = 1.e20 69 varname = 'grad of '+varname 70 varunit = varunit+'/m' 71 case strupcase(vargrid) of 72 'T':BEGIN 73 case direc of 74 'x':BEGIN 75 divi = e1u[firstx:lastx, firsty:lasty] 76 newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz] 77 vargrid = 'U' 78 firstxu = firstxt & lastxu = lastxt & nxu = nxt 79 firstyu = firstyt & lastyu = lastyt & nyu = nyt 80 END 81 'y':BEGIN 82 divi = e2v[firstx:lastx, firsty:lasty] 83 newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz] 84 vargrid = 'V' 85 firstxv = firstxt & lastxv = lastxt & nxv = nxt 86 firstyv = firstyt & lastyv = lastyt & nyv = nyt 87 END 88 'z':BEGIN 89 divi = e3w[firstz:lastz] 90 newmask = mask 91 vargrid = 'W' 92 firstzw = firstzt & lastzw = lastzt & nzw = nzt 93 END 94 ELSE:return, report('Bad definition of direction argument') 95 ENDCASE 96 END 97 'W':BEGIN 98 case direc of 99 'x':BEGIN 100 divi = e1u[firstx:lastx, firsty:lasty] 101 newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz] 102 vargrid = 'U' 103 firstxu = firstxt & lastxu = lastxt & nxu = nxt 104 firstyu = firstyt & lastyu = lastyt & nyu = nyt 105 END 106 'y':BEGIN 107 divi = e2v[firstx:lastx, firsty:lasty] 108 newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz] 109 vargrid = 'V' 110 firstxv = firstxt & lastxv = lastxt & nxv = nxt 111 firstyv = firstyt & lastyv = lastyt & nyv = nyt 112 END 113 'z':BEGIN 114 divi = e3t[firstz:lastz] 115 newmask = mask 116 vargrid = 'T' 117 firstzt = firstzw & lastzt = lastzw & nzt = nzw 118 END 119 ELSE:return, report('Bad definition of direction argument') 120 endcase 121 END 122 'U':BEGIN 123 case direc of 124 'x':BEGIN 125 divi = (shift(e1t, -1, 0))[firstx:lastx, firsty:lasty] 126 newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz] 127 vargrid = 'T' 128 firstxt = firstxu & lastxt = lastxu & nxt = nxu 129 firstyt = firstyu & lastyt = lastyu & nyt = nyu 130 END 131 'y':BEGIN 132 divi = e2f[firstx:lastx, firsty:lasty] 133 newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz] 134 vargrid = 'F' 135 firstxf = firstxu & lastxf = lastxu & nxf = nxu 136 firstyf = firstyu & lastyf = lastyu & nyf = nyu 137 END 138 'z':BEGIN 139 divi = e3w[firstz:lastz] 140 newmask = mask 141 vargrid = 'W' 142 firstzw = firstzt & lastzw = lastzt & nzw = nzt 143 END 144 ELSE:return, report('Bad definition of direction argument') 145 endcase 146 END 147 'V':BEGIN 148 case direc of 149 'x':BEGIN 150 divi = e1f[firstx:lastx, firsty:lasty] 151 newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz] 152 vargrid = 'F' 153 firstxf = firstxv & lastxf = lastxv & nxf = nxv 154 firstyf = firstyv & lastyf = lastyv & nyf = nyv 155 END 156 'y':BEGIN 157 divi = (shift(e2t, 0, -1))[firstx:lastx, firsty:lasty] 158 newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz] 159 vargrid = 'T' 160 firstxt = firstxv & lastxt = lastxv & nxt = nxv 161 firstyt = firstyv & lastyt = lastyv & nyt = nyv 162 END 163 'z':BEGIN 164 divi = e3w[firstz:lastz] 165 newmask = mask 166 vargrid = 'W' 167 firstzw = firstzt & lastzw = lastzt & nzw = nzt 168 END 169 ELSE:return, report('Bad definition of direction argument') 170 endcase 171 END 172 'F':BEGIN 159 173 ; case direc of 160 174 ; 'x':divi = (shift(e1v, -1, 0))[firstx:lastx, firsty:lasty] … … 163 177 ; ELSE:return, report('Bad definition of direction argument') 164 178 ; endcase 165 ; END 166 ELSE:return, report('Bad definition of vargrid') 167 ENDCASE 168 res = fitintobox(res) 169 IF n_elements(res) EQ 1 AND res[0] EQ -1 THEN return, res 170 case 1 of 179 return, report('F grid: case not coded, please contact SAXO team...') 180 END 181 ELSE:return, report('Bad definition of vargrid') 182 ENDCASE 183 res = fitintobox(temporary(res)) 184 IF n_elements(res) EQ 1 AND res[0] EQ -1 THEN return, res 185 case 1 of 171 186 ;---------------------------------------------------------------------------- 172 187 ;xy 173 188 ;---------------------------------------------------------------------------- 174 taille[0] EQ 2:BEGIN175 earth = where(mask[*, *, firstz] EQ 0)176 if earth[0] NE -1 then res[earth] = !values.f_nan177 178 179 res = (shift(res, -1, 0)-res)/divi180 181 if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0)182 183 184 res = (shift(res, 0, -1)-res)/divi185 186 if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(res, 0, 1)187 188 189 190 earth = where(newmask[*, *, firstz] EQ 0)191 if earth[0] NE -1 then res[earth] = valmask192 189 szres[0] EQ 2:BEGIN 190 land = where((temporary(mask))[*, *, firstz] EQ 0) 191 if land[0] NE -1 then res[temporary(land)] = !values.f_nan 192 case direc of 193 'x':BEGIN 194 res = (shift(res, -1, 0)-res)/temporary(divi) 195 if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *] = !values.f_nan 196 if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(temporary(res), 1, 0) 197 END 198 'y':BEGIN 199 res = (shift(res, 0, -1)-res)/temporary(divi) 200 res[*, ny-1] = !values.f_nan 201 if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1) 202 END 203 ELSE:return, report('Bad definition of direction argument for the type of array') 204 ENDCASE 205 land = where((temporary(newmask))[*, *, firstz] EQ 0) 206 if land[0] NE -1 then res[temporary(land)] = valmask 207 END 193 208 ;---------------------------------------------------------------------------- 194 209 ;xyt 195 210 ;---------------------------------------------------------------------------- 196 taille[0] EQ 3 AND jpt NE 1:BEGIN 197 earth = where(mask[*, *, firstz] EQ 0) 198 if earth[0] NE -1 then BEGIN 199 earth = earth#replicate(1, jpt)+replicate(1, n_elements(earth))#(nx*ny*lindgen(jpt)) 200 res[earth[*]] = !values.f_nan 201 ENDIF 202 divi = divi[*]#replicate(1, jpt) 203 case direc of 204 'x':BEGIN 205 res = (shift(res, -1, 0, 0)-res)/divi 206 if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan 207 if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0, 0) 208 END 209 'y':BEGIN 210 res = (shift(res, 0, -1, 0)-res)/divi 211 res[*, ny-1, *] = !values.f_nan 212 if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(res, 0, 1, 0) 213 END 214 ELSE:return, report('Bad definition of direction argument for the type of array') 215 ENDCASE 216 earth = where(newmask[*, *, firstz] EQ 0) 217 if earth[0] NE -1 THEN res[earth] = valmask 218 END 211 szres[0] EQ 3 AND jpt NE 1:BEGIN 212 land = where((temporary(mask))[*, *, firstz] EQ 0, cnt) 213 if land[0] NE -1 then BEGIN 214 land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*lindgen(jpt)) 215 res[temporary(land)] = !values.f_nan 216 ENDIF 217 divi = (temporary(divi))[*]#replicate(1., jpt) 218 case direc of 219 'x':BEGIN 220 res = (shift(res, -1, 0, 0)-res)/temporary(divi) 221 if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan 222 if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(temporary(res), 1, 0, 0) 223 END 224 'y':BEGIN 225 res = (shift(res, 0, -1, 0)-res)/temporary(divi) 226 res[*, ny-1, *] = !values.f_nan 227 if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1, 0) 228 END 229 ELSE:return, report('Bad definition of direction argument for the type of array') 230 ENDCASE 231 land = where((temporary(newmask))[*, *, firstz] EQ 0, cnt) 232 if land[0] NE -1 then BEGIN 233 land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*lindgen(jpt)) 234 res[temporary(land)] = valmask 235 ENDIF 236 END 219 237 ;---------------------------------------------------------------------------- 220 238 ;xyz 221 239 ;---------------------------------------------------------------------------- 222 taille[0] EQ 3 AND jpt EQ 1:BEGIN223 earth= where(mask EQ 0)224 if earth[0] NE -1 then res[earth] = !values.f_nan225 226 227 divi = divi[*]#replicate(1, nz)228 res = (shift(res, -1, 0, 0)-res)/divi229 230 if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0, 0)231 232 233 divi = divi[*]#replicate(1, nz)234 res = (shift(res, 0, -1, 0)-res)/divi235 236 if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(res, 0, 1, 0)237 238 239 divi = reform(replicate(1, nx*ny)#divi, nx, ny, nz)240 if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz)241 242 res = (shift(res, 0, 0, 1)-res)/divi243 244 245 res = (res-shift(res, 0, 0, -1))/divi246 247 248 if earth[0] NE -1 then res[earth] = valmask249 END250 ENDCASE251 END252 ;------------------------------------------------------------ 240 szres[0] EQ 3 AND jpt EQ 1:BEGIN 241 land = where(mask EQ 0) 242 if land[0] NE -1 then res[temporary(land)] = !values.f_nan 243 case direc OF 244 'x':BEGIN 245 divi = (temporary(divi))[*]#replicate(1., nz) 246 res = (shift(res, -1, 0, 0)-res)/temporary(divi) 247 if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan 248 if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(temporary(res), 1, 0, 0) 249 END 250 'y':BEGIN 251 divi = (temporary(divi))[*]#replicate(1., nz) 252 res = (shift(res, 0, -1, 0)-res)/temporary(divi) 253 res[*, ny-1, *] = !values.f_nan 254 if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1, 0) 255 END 256 'z':BEGIN 257 divi = replicate(1., nx*ny)#(temporary(divi))[*] 258 if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz, /overwrite) 259 if vargrid EQ 'W' THEN BEGIN 260 res = (shift(res, 0, 0, 1)-res)/temporary(divi) 261 res[*, *, 0] = !values.f_nan 262 ENDIF ELSE BEGIN 263 res = (res-shift(res, 0, 0, -1))/temporary(divi) 264 res[*, *, nz-1] = !values.f_nan 265 ENDELSE 266 END 267 ENDCASE 268 land = where(temporary(newmask) EQ 0) 269 if land[0] NE -1 then res[temporary(land)] = valmask 270 END 253 271 ;---------------------------------------------------------------------------- 254 272 ;xyzt 255 273 ;---------------------------------------------------------------------------- 256 taille[0] EQ 4:BEGIN 257 earth = where((mask[*]#replicate(1, jpt)) EQ 0) 258 if earth[0] NE -1 then res[earth] = !values.f_nan 259 case direc OF 260 'x':BEGIN 261 divi = divi[*]#replicate(1, nz*jpt) 262 res = (shift(res, -1, 0, 0, 0)-res)/divi 263 if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *, *] = !values.f_nan 264 if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0, 0, 0) 265 END 266 'y':BEGIN 267 divi = divi[*]#replicate(1, nz*jpt) 268 res = (shift(res, 0, -1, 0, 0)-res)/divi 269 res[*, ny-1, *, *] = !values.f_nan 270 if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(res, 0, 1, 0, 0) 271 END 272 'z':BEGIN 273 divi = replicate(1, nx*ny)#divi 274 divi = reform(divi[*]#replicate(1, jpt), nx, ny, nz, jpt, /over) 275 if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz, jpt) 276 if vargrid EQ 'W' THEN BEGIN 277 res = (shift(res, 0, 0, 1, 0)-res)/divi 278 res[*, *, 0, *] = !values.f_nan 279 ENDIF ELSE BEGIN 280 res = (res-shift(res, 0, 0, -1, 0))/divi 281 res[*, *, nz-1, *] = !values.f_nan 282 ENDELSE 283 END 284 ENDCASE 285 if earth[0] NE -1 then res[earth] = valmask 286 END 287 ;------------------------------------------------------------ 288 ;------------------------------------------------------------ 289 endcase 290 varname = 'grad of '+varname 291 varunit = varunit+'/m' 292 293 294 ;------------------------------------------------------------ 295 return, res 296 end 297 298 299 300 301 302 274 szres[0] EQ 4:BEGIN 275 land = where(temporary(mask) EQ 0, cnt) 276 if land[0] NE -1 then BEGIN 277 land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*nz*lindgen(jpt)) 278 res[temporary(land)] = !values.f_nan 279 ENDIF 280 case direc OF 281 'x':BEGIN 282 divi = (temporary(divi))[*]#replicate(1., nz*jpt) 283 res = (shift(res, -1, 0, 0, 0)-res)/temporary(divi) 284 if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *, *] = !values.f_nan 285 if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(temporary(res), 1, 0, 0, 0) 286 END 287 'y':BEGIN 288 divi = (temporary(divi))[*]#replicate(1., nz*jpt) 289 res = (shift(res, 0, -1, 0, 0)-res)/temporary(divi) 290 res[*, ny-1, *, *] = !values.f_nan 291 if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1, 0, 0) 292 END 293 'z':BEGIN 294 divi = replicate(1., nx*ny)#(temporary(divi))[*] 295 divi = (temporary(divi))[*]#replicate(1L, jpt) 296 if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz, jpt, /overwrite) 297 if vargrid EQ 'W' THEN BEGIN 298 res = (shift(res, 0, 0, 1, 0)-res)/temporary(divi) 299 res[*, *, 0, *] = !values.f_nan 300 ENDIF ELSE BEGIN 301 res = (res-shift(res, 0, 0, -1, 0))/temporary(divi) 302 res[*, *, nz-1, *] = !values.f_nan 303 ENDELSE 304 END 305 ENDCASE 306 land = where(newmask EQ 0, cnt) 307 if land[0] NE -1 then BEGIN 308 land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*nz*lindgen(jpt)) 309 res[temporary(land)] = valmask 310 ENDIF 311 END 312 ;------------------------------------------------------------ 313 ;------------------------------------------------------------ 314 ELSE:return, report('input array must have 2, 3 or 4 dimensions') 315 ENDCASE 316 ;------------------------------------------------------------ 317 318 319 ;------------------------------------------------------------ 320 return, res 321 END 322 323 324 325 326 327
Note: See TracChangeset
for help on using the changeset viewer.
