Changeset 25 for trunk/ToBeReviewed/CALCULS/curl.pro
- Timestamp:
- 05/02/06 14:59:12 (18 years ago)
- File:
-
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
trunk/ToBeReviewed/CALCULS/curl.pro
r23 r25 61 61 tempsun = systime(1) ; pour key_performance 62 62 ; 63 IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $ 64 return, report(['This version of curl is based on Arakawa C-grid.' $ 65 , 'U and V grids must therefore be defined']) 66 ; 63 67 u = litchamp(uu) 64 68 v = litchamp(vv) … … 71 75 ; on trouve les points que u et v ont en communs 72 76 ;------------------------------------------------------------ 73 indicexu = (lindgen(jpi))[ premierxu:premierxu+nxu-1]74 indicexv = (lindgen(jpi))[ premierxv:premierxv+nxv-1]77 indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1] 78 indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1] 75 79 indicex = inter(indicexu, indicexv) 76 indiceyu = (lindgen(jpj))[ premieryu:premieryu+nyu-1]77 indiceyv = (lindgen(jpj))[ premieryv:premieryv+nyv-1]80 indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1] 81 indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1] 78 82 indicey = inter(indiceyu, indiceyv) 79 83 nx = n_elements(indicex) … … 95 99 (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 96 100 case 1 of 97 nxu NE nx:if indicex[0] EQ premierxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]98 nxv NE nx:if indicex[0] EQ premierxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]99 nyu NE ny:if indicey[0] EQ premieryu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]100 nyv NE ny:if indicey[0] EQ premieryv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]101 nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 102 nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 103 nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 104 nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 101 105 ELSE : 102 106 endcase … … 114 118 coefu = (e1u[indice2d])[*]#replicate(1, nzt) 115 119 coefu = reform(coefu, nx, ny, nzt, /over) 116 coefu = coefu*(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt]120 coefu = coefu*(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 117 121 terreu = where(coefu EQ 0) 118 122 if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan … … 120 124 coefv = (e2v[indice2d])[*]#replicate(1, nzt) 121 125 coefv = reform(coefv, nx, ny, nzt, /over) 122 coefv = coefv*(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt]126 coefv = coefv*(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 123 127 terrev = where(coefv EQ 0) 124 128 if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan 125 129 126 tabf = (fmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt]130 tabf = (fmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 127 131 div = (e1f[indice2d]*e2f[indice2d])[*]#replicate(1, nzt) 128 132 div = reform(div, nx, ny, nzt, /over) … … 137 141 ; mise a !values.f_nan de la bordure 138 142 ;------------------------------------------------------------ 139 if NOT keyword_set(key_periodi que) OR nx NE jpi then begin143 if NOT keyword_set(key_periodic) OR nx NE jpi then begin 140 144 psi(0, *, *) = !values.f_nan 141 145 psi(nx-1, *, *) = !values.f_nan … … 150 154 ; pour le trace graphique 151 155 ;------------------------------------------------------------ 152 domdef, (gla mt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, grille = ['t', 'f']153 if keyword_set(direc) then psi = moyenne(psi,direc,/nan , boite = boite)156 domdef, (glagmt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f'] 157 if keyword_set(direc) then psi = moyenne(psi,direc,/nan) 154 158 155 159 ;---------------------------------------------------------------------------- … … 170 174 (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 171 175 if nxu NE nx then $ 172 if indicex[0] EQ premierxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]176 if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 173 177 IF nxv NE nx THEN $ 174 if indicex[0] EQ premierxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]178 if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 175 179 IF nyu NE ny THEN $ 176 if indicey[0] EQ premieryu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]180 if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 177 181 IF nyv NE ny THEN $ 178 if indicey[0] EQ premieryv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]182 if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 179 183 END 180 184 (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ … … 191 195 ; calcul du rotationnel 192 196 ;---------------------------------------------------------------------------- 193 coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj* premierzt]197 coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt] 194 198 terreu = where(coefu EQ 0) 195 199 if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan … … 197 201 coefu = reform(coefu, nx, ny, jpt, /over) 198 202 ; 199 coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj* premierzt]203 coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt] 200 204 terrev = where(coefv EQ 0) 201 205 if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan … … 203 207 coefv = reform(coefv, nx, ny, jpt, /over) 204 208 ; 205 tabf = (fmask())[indice2d+jpi*jpj* premierzt]/(e1f[indice2d]*e2f[indice2d])209 tabf = (fmask())[indice2d+jpi*jpj*firstzt]/(e1f[indice2d]*e2f[indice2d]) 206 210 tabf = temporary(tabf[*])#replicate(1, jpt) 207 211 tabf = reform(tabf, nx, ny, jpt, /over) … … 217 221 ; mise a !values.f_nan de la bordure 218 222 ;------------------------------------------------------------ 219 if NOT keyword_set(key_periodi que) OR nx NE jpi then begin223 if NOT keyword_set(key_periodic) OR nx NE jpi then begin 220 224 psi(0, *, *) = !values.f_nan 221 225 psi(nx-1, *, *) = !values.f_nan … … 227 231 if terref[0] NE -1 then psi[temporary(terref)] = valmask 228 232 ;---------------------------------------------------------------------------- 229 domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, grille = ['t', 'f']230 if keyword_set(direc) then psi = grossemoyenne(psi,direc,/nan , boite = boite)233 domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f'] 234 if keyword_set(direc) then psi = grossemoyenne(psi,direc,/nan) 231 235 ;---------------------------------------------------------------------------- 232 236 END … … 255 259 (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 256 260 if nxu NE nx then $ 257 if indicex[0] EQ premierxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *]261 if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *] 258 262 IF nxv NE nx THEN $ 259 if indicex[0] EQ premierxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *]263 if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *] 260 264 IF nyu NE ny THEN $ 261 if indicey[0] EQ premieryu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny]265 if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny] 262 266 IF nyv NE ny THEN $ 263 if indicey[0] EQ premieryv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny]267 if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny] 264 268 END 265 269 (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ … … 273 277 ; calcul du rotationnel 274 278 ;------------------------------------------------------------ 275 coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj* premierzt]279 coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt] 276 280 terreu = where(coefu EQ 0) 277 281 if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan 278 coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj* premierzt]282 coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt] 279 283 terrev = where(coefv EQ 0) 280 284 if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan 281 tabf = (fmask())[indice2d+jpi*jpj* premierzt]/(e1f[indice2d]*e2f[indice2d])285 tabf = (fmask())[indice2d+jpi*jpj*firstzt]/(e1f[indice2d]*e2f[indice2d]) 282 286 ; 283 287 zu = u*temporary(coefu) … … 290 294 ; mise a !values.f_nan de la bordure 291 295 ;------------------------------------------------------------ 292 if NOT keyword_set(key_periodi que) OR nx NE jpi then begin296 if NOT keyword_set(key_periodic) OR nx NE jpi then begin 293 297 psi(0, *) = !values.f_nan 294 298 psi(nx-1, *) = !values.f_nan … … 303 307 ; pour le trace graphique 304 308 ;------------------------------------------------------------ 305 domdef, (glamt[indice2d])[0, 0], (glamf[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphif[indice2d])[0, ny-1], prof1, prof2, grille = ['t', 'f'] 306 if keyword_set(direc) then psi = moyenne(psi,direc,/nan, boite = boite) 307 309 domdef, (glamt[indice2d])[0, 0], (glamf[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphif[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f'] 310 if keyword_set(direc) then psi = moyenne(psi,direc,/nan) 308 311 309 312 END
Note: See TracChangeset
for help on using the changeset viewer.