Changeset 226 for trunk/SRC/Computation/curl.pro
- Timestamp:
- 03/16/07 10:22:26 (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/SRC/Computation/curl.pro
r168 r226 9 9 ; @categories 10 10 ; Calculation 11 ; 11 ; 12 12 ; @param UU 13 13 ; Matrix representing the zonal coordinates (U point) of a field of vectors 14 14 ; A 2D (xy), 3D (xyz or yt) or a structure readable by litchamp and containing 15 15 ; a 2D (xy), 3D (xyz or yt) array (4d case is not coded yet). 16 ; note that the dimension of the arr y must suit the domain dimension.17 ; 18 ; @param VV 16 ; note that the dimension of the array must suit the domain dimension. 17 ; 18 ; @param VV 19 19 ; Matrix representing the meridional coordinates (V point) of a field of vectors 20 20 ; A 2D (xy), 3D (xyz or yt) or a structure readable by litchamp and containing 21 21 ; a 2D (xy), 3D (xyz or yt) array (4d case is not coded yet). 22 ; note that the dimension of the arr y must suit the domain dimension.22 ; note that the dimension of the array must suit the domain dimension. 23 23 ; 24 24 ; @keyword DIREC {type=scalar string} … … 30 30 ; 31 31 ; @uses 32 ; cm_4cal, cm_4data, cm_4mmesh 32 ; cm_4cal, cm_4data, cm_4mmesh 33 33 ; 34 34 ; @restrictions 35 35 ; 36 ; - Works only for Arakawa C-grid. 36 ; - Works only for Arakawa C-grid. 37 37 ; - UU must be on U grid, VV must be on V grid 38 38 ; - 4d case is not coded yet 39 39 ; - the common variable jpt is used to differ xyz (jpt=1) and xyt (jpt\=1) cases. 40 ; - U and V arrays are cut in the same geographic domain. Because of the shift between 41 ; T, U, V and F grids, it is possible that these two arrays do not have the same 42 ; size and refer to different indexes. In this case, arrays are re-cut on 40 ; - U and V arrays are cut in the same geographic domain. Because of the shift between 41 ; T, U, V and F grids, it is possible that these two arrays do not have the same 42 ; size and refer to different indexes. In this case, arrays are re-cut on 43 43 ; common indexes. To avoid these re-cuts, use the keyword /memeindice in domdef.pro 44 44 ; - When computing the divergence, we update, vargrid, varname, varunits and the 45 45 ; grid position parameters (firstxf, lastxf, nxf, firstyf, lastyf, nyf). 46 ; - points that cannot be computed (domain bo ndaries, coastline) are set to NaN46 ; - points that cannot be computed (domain boundaries, coastline) are set to NaN 47 47 ; 48 48 ; @examples … … 50 50 ; IDL> plt, curl(dist(jpi,jpj), dist(jpi,jpj)) 51 51 ; 52 ; @history 52 ; @history 53 53 ; Guillaume Roullet (grlod\@ipsl.jussieu.fr) 54 54 ; Sebastien Masson (smasson\@lodyc.jussieu.fr) 55 55 ; adaptation to work with a reduce domain 56 56 ; 21/5/1999: missing values at !values.f_nan 57 ; 57 ; 58 58 ; @version 59 59 ; $Id$ 60 60 ; 61 ; @todo 61 ; @todo 62 62 ; code the 4d case 63 63 ;- … … 71 71 @cm_4cal ; for jpt 72 72 @cm_4data ; for varname, vargrid, vardate, varunit, valmask 73 @cm_4mesh 73 @cm_4mesh 74 74 ; 75 75 tempsun = systime(1) ; To key_performance … … 84 84 szu = size(u) 85 85 szv = size(v) 86 86 87 87 if szu[0] NE szv[0] then return, report('U and V input data must have the same number of dimensions') 88 88 … … 96 96 indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1] 97 97 indicey = inter(indiceyu, indiceyv) 98 nx = n_elements(indicex) 98 nx = n_elements(indicex) 99 99 ny = n_elements(indicey) 100 100 indice2d = lindgen(jpi, jpj) … … 105 105 varunits = 's-1' 106 106 if n_elements(valmask) EQ 0 THEN valmask = 1e20 107 firstxt = indicex[0] & lastxt = indicex[0]+nx-1 & nxt = nx 107 firstxt = indicex[0] & lastxt = indicex[0]+nx-1 & nxt = nx 108 108 firstyt = indicey[0] & lastyt = indicey[0]+ny-1 & nyt = ny 109 109 ;---------------------------------------------------------------------------- … … 113 113 ;xyz 114 114 ;---------------------------------------------------------------------------- 115 szu[0] EQ 3 AND jpt EQ 1:BEGIN 115 szu[0] EQ 3 AND jpt EQ 1:BEGIN 116 116 ;------------------------------------------------------------ 117 117 ; extraction of U and V on the appropriated domain … … 119 119 case 1 of 120 120 szu[1] EQ nxu AND szu[2] EQ nyu AND $ 121 szv[1] EQ nxv AND szv[2] EQ nyv:BEGIN 121 szv[1] EQ nxv AND szv[2] EQ nyv:BEGIN 122 122 case 1 of 123 123 nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] … … 129 129 END 130 130 szu[1] EQ jpi AND szu[2] EQ jpj AND $ 131 szv[1] EQ jpi AND szv[2] EQ jpj:BEGIN 131 szv[1] EQ jpi AND szv[2] EQ jpj:BEGIN 132 132 u = u[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, *] 133 133 v = v[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, *] … … 142 142 landu = where(coefu EQ 0) 143 143 if landu[0] NE -1 then coefu[temporary(landu)] = !values.f_nan 144 144 145 145 coefv = ((e2v[indice2d])[*]#replicate(1., nzt)) $ 146 146 *(vmask())[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt] … … 175 175 ;---------------------------------------------------------------------------- 176 176 ;---------------------------------------------------------------------------- 177 szu[0] EQ 3 AND jpt GT 1:BEGIN 177 szu[0] EQ 3 AND jpt GT 1:BEGIN 178 178 ;------------------------------------------------------------ 179 179 ; extraction of U and V on the appropriated domain … … 181 181 case 1 of 182 182 szu[1] EQ nxu AND szu[2] EQ nyu AND $ 183 szv[1] EQ nxv AND szv[2] EQ nyv:BEGIN 183 szv[1] EQ nxv AND szv[2] EQ nyv:BEGIN 184 184 if nxu NE nx then $ 185 if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 185 if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 186 186 IF nxv NE nx THEN $ 187 187 if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 188 188 IF nyu NE ny THEN $ 189 if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 189 if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 190 190 IF nyv NE ny THEN $ 191 191 if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 192 192 END 193 193 szu[1] EQ jpi AND szu[2] EQ jpj AND $ 194 szv[1] EQ jpi AND szv[2] EQ jpj:BEGIN 194 szv[1] EQ jpi AND szv[2] EQ jpj:BEGIN 195 195 u = u[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, *] 196 196 v = v[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, *] 197 197 END 198 ELSE:BEGIN 198 ELSE:BEGIN 199 199 print, 'problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs' 200 200 return, -1 … … 240 240 ;---------------------------------------------------------------------------- 241 241 ;---------------------------------------------------------------------------- 242 szu[0] EQ 4:BEGIN 242 szu[0] EQ 4:BEGIN 243 243 return, report('Case not coded contact saxo team or make a do loop!') 244 244 END … … 248 248 ;---------------------------------------------------------------------------- 249 249 ;---------------------------------------------------------------------------- 250 szu[0] EQ 2:BEGIN 250 szu[0] EQ 2:BEGIN 251 251 ;------------------------------------------------------------ 252 252 ;------------------------------------------------------------ 253 253 case 1 of 254 254 szu[1] EQ nxu AND szu[2] EQ nyu AND $ 255 szv[1] EQ nxv AND szv[2] EQ nyv:BEGIN 255 szv[1] EQ nxv AND szv[2] EQ nyv:BEGIN 256 256 if nxu NE nx then $ 257 if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *] 257 if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *] 258 258 IF nxv NE nx THEN $ 259 259 if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *] 260 260 IF nyu NE ny THEN $ 261 if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny] 261 if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny] 262 262 IF nyv NE ny THEN $ 263 263 if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny] 264 264 END 265 265 szu[1] EQ jpi AND szu[2] EQ jpj AND $ 266 szv[1] EQ jpi AND szv[2] EQ jpj:BEGIN 266 szv[1] EQ jpi AND szv[2] EQ jpj:BEGIN 267 267 u = u[indice2d] 268 268 v = v[indice2d] … … 305 305 ENDCASE 306 306 ;------------------------------------------------------------ 307 if keyword_set(key_performance) THEN print, 'temps curl', systime(1)-tempsun 307 if keyword_set(key_performance) THEN print, 'temps curl', systime(1)-tempsun 308 308 309 309 return, psi
Note: See TracChangeset
for help on using the changeset viewer.