- Timestamp:
- 05/02/06 14:59:12 (18 years ago)
- Location:
- trunk
- Files:
-
- 3 added
- 1 deleted
- 17 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 -
trunk/ToBeReviewed/CALCULS/depth2level.pro
r23 r25 60 60 ;------------------------------------------------------------ 61 61 in = litchamp(tab) 62 grille,mask,glam,gphi,gdep,nx,ny,nz,premierx,premiery,premierz,dernierx,derniery,dernierz 63 glam = 1 64 gphi = 1 62 grille,mask,-1,-1,gdep,nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz 65 63 ;--------------------------------------------------------------- 66 64 ; verification de la coherence entre la taille du tableau et le domaine definit par domdef … … 70 68 if taille[0] NE 2 then return, report('le champ en entree doit contenir un tableau 2d') 71 69 case 1 of 72 taille[1] eq jpi and taille[2] eq jpj:in=in[ premierx:dernierx, premiery:derniery]70 taille[1] eq jpi and taille[2] eq jpj:in=in[firstx:lastx, firsty:lasty] 73 71 taille[1] eq nx and taille[2] eq ny: 74 72 else:return, report('Probleme d''adequation entre les tailles du domaine et celle du champ.') … … 85 83 ;------------------------------------------------------------ 86 84 ; on passe en tableaux qui ont la taille des tableaux 3d 87 prof=replicate(1,nx*ny)#gdep[ premierz:dernierz]85 prof=replicate(1,nx*ny)#gdep[firstz:lastz] 88 86 in = in[*]#replicate(1, nz) 89 87 ; … … 94 92 if keyword_set(upper) then begin 95 93 levels = levels-1 96 notvalid = where(levels EQ 0)94 notvalid = where(levels EQ -1) 97 95 ENDIF ELSE notvalid = where(levels EQ nz) 98 96 IF notvalid[0] NE -1 THEN levels[notvalid] = !values.f_nan -
trunk/ToBeReviewed/CALCULS/div.pro
r23 r25 56 56 tempsun = systime(1) ; pour key_performance 57 57 @common 58 u = uu 59 v = vv 60 58 ; 59 IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $ 60 return, report(['This version of div is based on Arakawa C-grid.' $ 61 , 'U and V grids must therefore be defined']) 62 ; 63 u = litchamp(uu) 64 v = litchamp(vv) 65 ; 61 66 date1 = time[0] 62 67 if n_elements(jpt) EQ 0 then date2 = date1 ELSE date2 = time[jpt-1] … … 66 71 ; on trouve les points que u et v ont en communs 67 72 ;------------------------------------------------------------ 68 indicexu = (lindgen(jpi))[ premierxu:premierxu+nxu-1]69 indicexv = (lindgen(jpi))[ premierxv:premierxv+nxv-1]73 indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1] 74 indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1] 70 75 indicex = inter(indicexu, indicexv) 71 indiceyu = (lindgen(jpj))[ premieryu:premieryu+nyu-1]72 indiceyv = (lindgen(jpj))[ premieryv:premieryv+nyv-1]76 indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1] 77 indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1] 73 78 indicey = inter(indiceyu, indiceyv) 74 79 nx = n_elements(indicex) … … 86 91 ;------------------------------------------------------------ 87 92 case 1 of 88 (size( u))[0] NE 3 OR (size(v))[0] NE 3: return, -193 (size(v))[0] NE 3: return, -1 89 94 (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $ 90 95 (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 91 96 case 1 of 92 nxu NE nx:if indicex[0] EQ premierxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]93 nxv NE nx:if indicex[0] EQ premierxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]94 nyu NE ny:if indicey[0] EQ premieryu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]95 nyv NE ny:if indicey[0] EQ premieryv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]97 nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 98 nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 99 nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 100 nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 96 101 ELSE : 97 102 endcase … … 114 119 zu = reform(zu, nx, ny, nzt, /over) 115 120 zu = temporary(u)*temporary(zu) $ 116 *(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt]121 *(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 117 122 terreu = where(zu EQ 0) 118 123 if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan … … 121 126 zv = reform(zv, nx, ny, nzt, /over) 122 127 zv = temporary(v)*temporary(zv) $ 123 *(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt]128 *(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 124 129 terrev = where(zv EQ 0) 125 130 if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan … … 129 134 zdiv = reform(zdiv, nx, ny, nzt, /over) 130 135 zdiv = temporary(zdiv)*(zu-shift(zu, 1, 0, 0)+zv-shift(zv, 0, 1, 0)) $ 131 *tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt]136 *tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 132 137 ;------------------------------------------------------------ 133 138 ; mise a !values.f_nan de la bordure 134 139 ;------------------------------------------------------------ 135 if NOT keyword_set(key_periodi que) OR nx NE jpi then begin140 if NOT keyword_set(key_periodic) OR nx NE jpi then begin 136 141 zdiv(0, *, *) = !values.f_nan 137 142 zdiv(nx-1, *, *) = !values.f_nan … … 143 148 ; 144 149 if n_elements(valmask) EQ 0 THEN valmask = 1e20 145 terre = where(tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt] EQ 0)150 terre = where(tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] EQ 0) 146 151 if terre[0] NE -1 then zdiv[temporary(terre)] = valmask 147 152 ;------------------------------------------------------------ … … 151 156 varname = 'div' 152 157 varunits = '1e6*s-1' 153 domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, grille = ['t']154 if keyword_set(direc) then zdiv = moyenne(zdiv,direc,/nan , boite = boite)158 domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t'] 159 if keyword_set(direc) then zdiv = moyenne(zdiv,direc,/nan) 155 160 END 156 161 ;---------------------------------------------------------------------------- … … 168 173 (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 169 174 case 1 of 170 nxu NE nx:if indicex[0] EQ premierxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]171 nxv NE nx:if indicex[0] EQ premierxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]172 nyu NE ny:if indicey[0] EQ premieryu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]173 nyv NE ny:if indicey[0] EQ premieryv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]175 nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 176 nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 177 nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 178 nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 174 179 ELSE : 175 180 endcase … … 185 190 ; calcul de la divergence 186 191 ;------------------------------------------------------------ 187 zu = e2u[indice2d]*(umask())[indice2d+jpi*jpj* premierzt]192 zu = e2u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt] 188 193 terreu = where(zu EQ 0) 189 194 if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan … … 192 197 zu = temporary(u)*temporary(zu) 193 198 ; 194 zv = e1v[indice2d]*(vmask())[indice2d+jpi*jpj* premierzt]199 zv = e1v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt] 195 200 terrev = where(zv EQ 0) 196 201 if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan … … 199 204 zv = temporary(v)*temporary(zv) 200 205 ; 201 zdiv = 1e6*tmask[indice2d+jpi*jpj* premierzt]/(e1t[indice2d]*e2t[indice2d])206 zdiv = 1e6*tmask[indice2d+jpi*jpj*firstzt]/(e1t[indice2d]*e2t[indice2d]) 202 207 zdiv = (zdiv)[*]#replicate(1, jpt) 203 208 zdiv = reform(zdiv, nx, ny, jpt, /over) … … 207 212 ; mise a !values.f_nan de la bordure 208 213 ;------------------------------------------------------------ 209 if NOT keyword_set(key_periodi que) OR nx NE jpi then begin214 if NOT keyword_set(key_periodic) OR nx NE jpi then begin 210 215 zdiv(0, *, *) = !values.f_nan 211 216 zdiv(nx-1, *, *) = !values.f_nan … … 222 227 varname = 'div' 223 228 varunits = '1e6*s-1' 224 domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, grille = ['t']225 if keyword_set(direc) then zdiv = grossemoyenne(zdiv,direc,/nan , boite = boite)229 domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t'] 230 if keyword_set(direc) then zdiv = grossemoyenne(zdiv,direc,/nan) 226 231 END 227 232 ;---------------------------------------------------------------------------- … … 240 245 ELSE:BEGIN ;xy 241 246 indice3d = lindgen(jpi, jpj, jpk) 242 indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, niveau -1]247 indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt] 243 248 ;------------------------------------------------------------ 244 249 ; extraction de u et v sur le domaine qui convient … … 252 257 (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 253 258 case 1 of 254 nxu NE nx:if indicex[0] EQ premierxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *]255 nxv NE nx:if indicex[0] EQ premierxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *]256 nyu NE ny:if indicey[0] EQ premieryu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny]257 nyv NE ny:if indicey[0] EQ premieryv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny]259 nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *] 260 nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *] 261 nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny] 262 nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny] 258 263 ELSE : 259 264 endcase … … 280 285 ; mise a !values.f_nan de la bordure 281 286 ;------------------------------------------------------------ 282 if NOT keyword_set(key_periodi que) OR nx NE jpi then begin287 if NOT keyword_set(key_periodic) OR nx NE jpi then begin 283 288 zdiv(0, *) = !values.f_nan 284 289 zdiv(nx-1, *) = !values.f_nan … … 298 303 varname = 'div' 299 304 varunits = '1e6*s-1' 300 domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, grille = ['t']301 if keyword_set(direc) then zdiv = moyenne(zdiv,direc,/nan , boite = boite)305 domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t'] 306 if keyword_set(direc) then zdiv = moyenne(zdiv,direc,/nan) 302 307 303 308 END -
trunk/ToBeReviewed/CALCULS/grad.pro
r23 r25 34 34 @common 35 35 ;------------------------------------------------------------ 36 ; 37 IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $ 38 return, report(['This version of grad is based on Arakawa C-grid.' $ 39 , 'U and V grids must therefore be defined']) 40 ; 36 41 res = litchamp(field) 37 42 taille=size(res) 38 grille, mask, glam, gphi, gdep, nx, ny,nz, premierx,premiery,premierz,dernierx, derniery, dernierz43 grille, mask, glam, gphi, gdep, nx, ny,nz,firstx,firsty,firstz,lastx, lasty, lastz 39 44 if n_elements(valmask) EQ 0 then valmask = 1e20 40 45 case strupcase(vargrid) of … … 42 47 case direc of 43 48 'x':BEGIN 44 divi = e1u[ premierx:dernierx, premiery:derniery]45 newmask = (umask())[ premierx:dernierx, premiery:derniery, premierz:dernierz]49 divi = e1u[firstx:lastx, firsty:lasty] 50 newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz] 46 51 vargrid = 'U' 47 domdef, glamt[ premierx, 0], glamu[dernierx, 0] $48 , gphit[0, premiery], gphiu[0, derniery], grille = ['T','U']52 domdef, glamt[firstx, 0], glamu[lastx, 0] $ 53 , gphit[0, firsty], gphiu[0, lasty], gridtype = ['T','U'] 49 54 END 50 55 'y':BEGIN 51 divi = e2v[ premierx:dernierx, premiery:derniery]52 newmask = (vmask())[ premierx:dernierx, premiery:derniery, premierz:dernierz]56 divi = e2v[firstx:lastx, firsty:lasty] 57 newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz] 53 58 vargrid = 'V' 54 domdef, glamt[ premierx, 0], glamv[dernierx, 0] $55 , gphit[0, premiery], gphiv[0, derniery], grille = ['T','V']56 END 57 'z':BEGIN 58 divi = e3w[ premierz:dernierz]59 domdef, glamt[firstx, 0], glamv[lastx, 0] $ 60 , gphit[0, firsty], gphiv[0, lasty], gridtype = ['T','V'] 61 END 62 'z':BEGIN 63 divi = e3w[firstz:lastz] 59 64 newmask = mask 60 65 vargrid = 'W' … … 65 70 'W':BEGIN 66 71 case direc of 67 'x':divi = e1u[ premierx:dernierx, premiery:derniery]68 'y':divi = e2v[ premierx:dernierx, premiery:derniery]69 'z':BEGIN 70 divi = e3t[ premierz:dernierz]72 'x':divi = e1u[firstx:lastx, firsty:lasty] 73 'y':divi = e2v[firstx:lastx, firsty:lasty] 74 'z':BEGIN 75 divi = e3t[firstz:lastz] 71 76 newmask = mask 72 77 vargrid = 'T' … … 78 83 case direc of 79 84 'x':BEGIN 80 divi = (shift(e1t, -1, 0))[ premierx:dernierx, premiery:derniery]81 newmask = tmask[ premierx:dernierx, premiery:derniery, premierz:dernierz]85 divi = (shift(e1t, -1, 0))[firstx:lastx, firsty:lasty] 86 newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz] 82 87 vargrid = 'T' 83 domdef, glamt[ premierx, 0], glamu[dernierx] $84 , gphit[0, premiery], gphiu[0, derniery], grille = ['T','U']88 domdef, glamt[firstx, 0], glamu[lastx] $ 89 , gphit[0, firsty], gphiu[0, lasty], gridtype = ['T','U'] 85 90 END 86 91 'y':BEGIN 87 divi = e2f[ premierx:dernierx, premiery:derniery]88 newmask = (fmask())[ premierx:dernierx, premiery:derniery, premierz:dernierz]92 divi = e2f[firstx:lastx, firsty:lasty] 93 newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz] 89 94 vargrid = 'F' 90 domdef, glamu[ premierx, 0], glamf[dernierx, 0] $91 , gphiu[0, premiery], gphif[0, derniery], grille = ['U','F']92 END 93 'z':BEGIN 94 divi = e3w[ premierz:dernierz]95 domdef, glamu[firstx, 0], glamf[lastx, 0] $ 96 , gphiu[0, firsty], gphif[0, lasty], gridtype = ['U','F'] 97 END 98 'z':BEGIN 99 divi = e3w[firstz:lastz] 95 100 newmask = mask 96 101 vargrid = 'W' … … 102 107 case direc of 103 108 'x':BEGIN 104 divi = e1f[ premierx:dernierx, premiery:derniery]105 newmask = (fmask())[ premierx:dernierx, premiery:derniery, premierz:dernierz]109 divi = e1f[firstx:lastx, firsty:lasty] 110 newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz] 106 111 vargrid = 'F' 107 domdef, glamv[ premierx, 0], glamf[dernierx, 0] $108 , gphiv[0, premiery], gphif[0, derniery], grille = ['V','F']112 domdef, glamv[firstx, 0], glamf[lastx, 0] $ 113 , gphiv[0, firsty], gphif[0, lasty], gridtype = ['V','F'] 109 114 END 110 115 'y':BEGIN 111 divi = (shift(e2t, 0, -1))[ premierx:dernierx, premiery:derniery]112 newmask = tmask[ premierx:dernierx, premiery:derniery, premierz:dernierz]116 divi = (shift(e2t, 0, -1))[firstx:lastx, firsty:lasty] 117 newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz] 113 118 vargrid = 'T' 114 domdef, glamt[ premierx, 0], glamv[dernierx, 0] $115 , gphit[0, premiery], gphiv[0, derniery], grille = ['T','V']116 END 117 'z':BEGIN 118 divi = e3w[ premierz:dernierz]119 domdef, glamt[firstx, 0], glamv[lastx, 0] $ 120 , gphit[0, firsty], gphiv[0, lasty], gridtype = ['T','V'] 121 END 122 'z':BEGIN 123 divi = e3w[firstz:lastz] 119 124 newmask = mask 120 125 vargrid = 'W' … … 125 130 ; 'F':BEGIN 126 131 ; case direc of 127 ; 'x':divi = (shift(e1v, -1, 0))[ premierx:dernierx, premiery:derniery]128 ; 'y':divi = (shift(e2u, 0, -1))[ premierx:dernierx, premiery:derniery]129 ; 'z':divi = e3w[ premierz:dernierz]132 ; 'x':divi = (shift(e1v, -1, 0))[firstx:lastx, firsty:lasty] 133 ; 'y':divi = (shift(e2u, 0, -1))[firstx:lastx, firsty:lasty] 134 ; 'z':divi = e3w[firstz:lastz] 130 135 ; ELSE:return, report('Bad definition of direction argument') 131 136 ; endcase … … 139 144 ;---------------------------------------------------------------------------- 140 145 taille[0] EQ 2:BEGIN 141 earth = where(mask[*, *, premierz] EQ 0)146 earth = where(mask[*, *, firstz] EQ 0) 142 147 if earth[0] NE -1 then res[earth] = !values.f_nan 143 148 case direc of 144 149 'x':BEGIN 145 150 res = (shift(res, -1, 0)-res)/divi 146 if key_periodi queEQ 0 OR nx NE jpi THEN res[nx-1, *] = !values.f_nan151 if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *] = !values.f_nan 147 152 if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0) 148 153 END … … 154 159 ELSE:return, report('Bad definition of direction argument for the type of array') 155 160 ENDCASE 156 earth = where(newmask[*, *, premierz] EQ 0)161 earth = where(newmask[*, *, firstz] EQ 0) 157 162 if earth[0] NE -1 then res[earth] = valmask 158 163 END … … 161 166 ;---------------------------------------------------------------------------- 162 167 taille[0] EQ 3 AND jpt NE 1:BEGIN 163 earth = where(mask[*, *, premierz] EQ 0)168 earth = where(mask[*, *, firstz] EQ 0) 164 169 if earth[0] NE -1 then BEGIN 165 170 earth = earth#replicate(1, jpt)+replicate(1, n_elements(earth))#(nx*ny*lindgen(jpt)) … … 170 175 'x':BEGIN 171 176 res = (shift(res, -1, 0, 0)-res)/divi 172 if key_periodi queEQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan177 if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan 173 178 if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0, 0) 174 179 END … … 180 185 ELSE:return, report('Bad definition of direction argument for the type of array') 181 186 ENDCASE 182 earth = where(newmask[*, *, premierz] EQ 0) 183 if earth[0] NE -1 then BEGIN 184 earth = earth#replicate(1, jpt)+replicate(1, n_elements(earth))#(nx*ny*lindgen(jpt)) 185 res[earth] = valmask 186 ENDIF 187 earth = where(newmask[*, *, firstz] EQ 0) 188 if earth[0] NE -1 THEN res[earth] = valmask 187 189 END 188 190 ;---------------------------------------------------------------------------- … … 196 198 divi = divi[*]#replicate(1, nz) 197 199 res = (shift(res, -1, 0, 0)-res)/divi 198 if key_periodi queEQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan200 if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan 199 201 if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0, 0) 200 202 END … … 227 229 if earth[0] NE -1 then res[earth] = !values.f_nan 228 230 case direc OF 231 'x':BEGIN 232 divi = divi[*]#replicate(1, nz*jpt) 233 res = (shift(res, -1, 0, 0, 0)-res)/divi 234 if key_periodic EQ 0 OR nx NE jpi THEN res[nx-1, *, *, *] = !values.f_nan 235 if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0, 0, 0) 236 END 237 'y':BEGIN 238 divi = divi[*]#replicate(1, nz*jpt) 239 res = (shift(res, 0, -1, 0, 0)-res)/divi 240 res[*, ny-1, *, *] = !values.f_nan 241 if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(res, 0, 1, 0, 0) 242 END 229 243 'z':BEGIN 230 244 divi = replicate(1, nx*ny)#divi … … 238 252 res[*, *, nz-1, *] = !values.f_nan 239 253 ENDELSE 240 if earth[0] NE -1 then res[earth] = valmask241 END242 ENDCASE254 END 255 ENDCASE 256 if earth[0] NE -1 then res[earth] = valmask 243 257 END 244 258 ;------------------------------------------------------------ -
trunk/ToBeReviewed/CALCULS/grossemoyenne.pro
r23 r25 11 11 ; CATEGORY: 12 12 ; 13 ; CALLING SEQUENCE: result = grossemoyenne(tab,'direc',BO ITE=boite)13 ; CALLING SEQUENCE: result = grossemoyenne(tab,'direc',BOXZOOM=boxzoom) 14 14 ; 15 15 ; INPUTS: tab = 3 or 4d field … … 19 19 ; KEYWORD PARAMETERS: 20 20 ; 21 ; bo ite= [xmin,xmax,ymin,ymax (,(zmin,)zmax)] pour plus21 ; boxzoom = [xmin,xmax,ymin,ymax (,(zmin,)zmax)] pour plus 22 22 ; de detail cf domdef 23 ; bo itepeut prendre 5 formes:24 ; [ prof2], [prof1, prof2],[lon1, lon2, lat1, lat2],25 ; [lon1, lon2, lat1, lat2, prof2],[lon1, lon2, lat1, lat2, prof1,prof2]23 ; boxzoom peut prendre 5 formes: 24 ; [vert2], [vert1, vert2],[lon1, lon2, lat1, lat2], 25 ; [lon1, lon2, lat1, lat2, vert2],[lon1, lon2, lat1, lat2, vert1,vert2] 26 26 ; 27 27 ; NAN: not a number, a activer si l''on peut faire veut … … 38 38 ; 39 39 ; NODOMDEF: activer si l''on ne veut pas passer ds 40 ; domdef bien que le mot cle bo itesoit present (comme40 ; domdef bien que le mot cle boxzoom soit present (comme 41 41 ; c''est le cas qd grossemoyenne est appelee via 42 42 ; checkfield) 43 43 ; 44 ; INTEGRATION: pour faire une integrale plutot qu''une moyenne 44 ; INTEGRATION: pour faire une integrale plutot qu''une 45 ; moyenne 46 ; 47 ; /SPATIALFIRST, when performing at the same time 48 ; spatial and temporal mean, grossemoyenne is assuming 49 ; that the mask is not changing with the time. In 50 ; consequence, grossemoyenne performs temporal mean 51 ; first and then call moyenne. Activate /SPATIALFIRST if 52 ; you want to perform the spatial mean before the 53 ; temporal mean. Note that if NAN is activated, then 54 ; SPATIALFIRST is activated automatically. 55 ; 56 ; /TEMPORALFIRST: to force to perform first temporal 57 ; mean even if nan is activated (see SPATIALFIRST 58 ; explanations...) 59 ; 60 ; 61 ; /WDEPTH: to specify that the field is at W depth instad of T 62 ; depth (automatically activated if vargrid eq 'W') 45 63 ; 46 64 ; OUTPUTS: … … 90 108 ;------------------------------------------------------------ 91 109 ;------------------------------------------------------------ 92 function grossemoyenne ,tab,direc,BOITE=boite, INTEGRATION=integration $ 93 , NAN = nan, NODOMDEF = nodomdef, _extra = ex 94 @common 95 tempsun = systime(1) ; pour key_performance 110 function grossemoyenne, tab, direc, BOXZOOM = boxzoom, INTEGRATION = integration $ 111 , NAN = nan, NODOMDEF = nodomdef, WDEPTH = wdepth $ 112 , SPATIALFIRST = spatialfirst, TEMPORALFIRST = temporalfirst $ 113 , _extra = ex 114 ;--------------------------------------------------------- 115 @cm_4mesh 116 @cm_4data 117 @cm_4cal 118 IF NOT keyword_set(key_forgetold) THEN BEGIN 119 @updatenew 120 @updatekwd 121 ENDIF 122 ;--------------------------------------------------------- 123 tempsun = systime(1) ; pour key_performance 96 124 ;------------------------------------------------------------ 97 125 ;I) preliminaires 98 126 ;------------------------------------------------------------ 99 100 101 102 103 127 dirt = 0 128 dirx = 0 129 diry = 0 130 dirz = 0 131 dim = 'aa' 104 132 ;------------------------------------------------------------ 105 133 ; I.1) direction(s) suivants lesquelles on integre 106 134 ;------------------------------------------------------------ 107 if ( strpos(direc,'t') ge 0 ) then dirt = 1 108 if ( strpos(direc,'x') ge 0 ) then dirx = 1 109 if ( strpos(direc,'y') ge 0 ) then diry = 1 110 if ( strpos(direc,'z') ge 0 ) then dirz = 1 135 if ( strpos(direc, 't') ge 0 ) then dirt = 1 136 if ( strpos(direc, 'x') ge 0 ) then dirx = 1 137 if ( strpos(direc, 'y') ge 0 ) then diry = 1 138 if ( strpos(direc, 'z') ge 0 ) then dirz = 1 139 IF keyword_set(NAN) AND (dirx EQ 1 OR diry EQ 1 OR dirz EQ 1) $ 140 THEN spatialfirst = 1 141 IF keyword_set(temporalfirst) THEN spatialfirst = 0 111 142 ;------------------------------------------------------------ 112 143 ; I.2) verification de la taille du tableau d'entree 113 144 ;------------------------------------------------------------ 114 115 116 117 118 119 120 121 122 123 124 125 126 127 145 taille = size(tab) 146 case 1 of 147 taille[0] eq 1 :return, report('Le tableau n''a qu''une dimension, cas non traite!') 148 taille[0] eq 2 :return, report('Le tableau n''a qu''deux dimension, cas non traite!') 149 taille[0] eq 3 :BEGIN 150 dim = '3d' 151 if dirx eq 0 and diry eq 0 and dirt eq 0 then return, tab 152 END 153 taille[0] eq 4 :BEGIN 154 dim = '4d' 155 if dirx eq 0 and diry eq 0 and dirz eq 0 and dirt eq 0 then return, tab 156 END 157 else : return, report('Le tableau d entree doit etre a 3 ou 4 dimensions s''il ne contient pas de dim temporelle utiliser moyenne') 158 endcase 128 159 ;------------------------------------------------------------ 129 160 ; I.4) obtention des facteurs d''echelles et du masque sur le sous 130 161 ; domaine concerne par la moyenne 131 ; redefinition du domaine ajuste a bo ite(a 6 elements)162 ; redefinition du domaine ajuste a boxzoom (a 6 elements) 132 163 ; ceci va nous permetre de faire les calcules que sur le sous domaine 133 164 ; comcerne par la moyenne. domdef, suivit de grille nous donne tous 134 165 ; les tableaux de la grille sur le sous domaine 135 166 ;------------------------------------------------------------ 136 if keyword_set(boite) then BEGIN 137 Case 1 Of 138 N_Elements(Boite) Eq 1: bte = [lon1, lon2, lat1, lat2, 0.,boite[0]] 139 N_Elements(Boite) Eq 2: bte = [lon1, lon2, lat1, lat2, boite[0],boite[1]] 140 N_Elements(Boite) Eq 4: bte = [Boite,prof1, prof2] 141 N_Elements(Boite) Eq 5: bte = [Boite[0:3], 0, Boite[4]] 142 N_Elements(Boite) Eq 6: bte = Boite 143 Else: return, report('Mauvaise Definition de Boite') 144 endcase 145 oldboite = [lon1, lon2, lat1, lat2, prof1, prof2] 146 if NOT keyword_set(nodomdef) then domdef, bte,GRILLE=vargrid, _extra = ex 147 ENDIF 167 if keyword_set(boxzoom) then BEGIN 168 Case 1 Of 169 N_Elements(Boxzoom) Eq 1: bte = [lon1, lon2, lat1, lat2, 0., boxzoom[0]] 170 N_Elements(Boxzoom) Eq 2: bte = [lon1, lon2, lat1, lat2, boxzoom[0], boxzoom[1]] 171 N_Elements(Boxzoom) Eq 4: bte = [Boxzoom, vert1, vert2] 172 N_Elements(Boxzoom) Eq 5: bte = [Boxzoom[0:3], 0, Boxzoom[4]] 173 N_Elements(Boxzoom) Eq 6: bte = Boxzoom 174 Else: return, report('Wrong Definition of Boxzoom') 175 endcase 176 if NOT keyword_set(nodomdef) then BEGIN 177 savedbox = 1b 178 saveboxparam, 'boxparam4grmoyenne.dat' 179 domdef, bte, GRIDTYPE = vargrid, _extra = ex 180 ENDIF 181 ENDIF 148 182 ;--------------------------------------------------------------- 149 183 ; attribution du mask et des tableaux de longitude et latitude... 150 184 ;--------------------------------------------------------------- 151 grille, mask, glam, gphi, gdep, nx, ny,nz,premierx,premiery,premierz, dernierx, derniery, dernierz,e1,e2,e3185 grille, mask, glam, gphi, gdep, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz, e1, e2, e3, WDEPTH = wdepth 152 186 ;------------------------------------------------------------ 153 187 ; I.3) si dirt eq 1 on fait la moyenne temporelle et on envoie ds moyenne 154 188 ;------------------------------------------------------------ 155 if dirt EQ 1then begin156 189 if dirt EQ 1 AND NOT keyword_set(spatialfirst) then begin 190 if dim EQ 3d then BEGIN 157 191 case 1 of 158 taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpt: $ 159 res = tab[premierx:premierx+nx-1 $ 160 ,premiery:premiery+ny-1, *] 161 taille[1] EQ nx and taille[2] eq ny and taille[3] eq jpt:res=tab 162 else:return, report('Probleme d''adequation entre les tailles du domaine nx*ny*jpt '+strtrim(nx, 1)+'*'+strtrim(ny, 1)+'*'+strtrim(jpt, 1)+' et du tableau '+strtrim(taille[1], 1)+'*'+strtrim(taille[2], 1)+'*'+strtrim(taille[3], 1)) 192 taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpt: $ 193 res = tab[firstx:firstx+nx-1 $ 194 , firsty:firsty+ny-1, *] 195 taille[1] EQ nx and taille[2] eq ny and taille[3] eq jpt:res = tab 196 else:BEGIN 197 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 198 return, report('Probleme d''adequation entre les tailles du domaine nx*ny*jpt '+strtrim(nx, 1)+'*'+strtrim(ny, 1)+'*'+strtrim(jpt, 1)+' et du tableau '+strtrim(taille[1], 1)+'*'+strtrim(taille[2], 1)+'*'+strtrim(taille[3], 1)) 199 END 163 200 ENDCASE 164 divi = finite(res) 165 divi = total(temporary(divi),3, nan = keyword_set(nan)) 166 notanum = where(divi EQ 0) 167 if keyword_set(integration) then begin 168 res = total(res,3, nan = nan) 169 ENDIF ELSE BEGIN 170 if keyword_set(nan) then BEGIN 171 res = total(res,3, nan = keyword_set(nan))/ (1 > divi) 172 ENDIF ELSE res = total(res,3)/(1.*taille[3]) 173 ENDELSE 201 if keyword_set(integration) then begin 202 res = total(res, 3, nan = nan) 174 203 ENDIF ELSE BEGIN 204 if keyword_set(nan) then BEGIN 205 divi = finite(res) 206 divi = total(temporary(divi), 3) 207 notanum = where(divi EQ 0) 208 res = total(res, 3, nan = keyword_set(nan))/ (1 > divi) 209 if notanum[0] NE -1 then res[temporary(notanum)] = !values.f_nan 210 ENDIF ELSE res = total(res, 3)/(1.*taille[3]) 211 ENDELSE 212 ENDIF ELSE BEGIN 175 213 case 1 of 176 taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpk and taille[4] eq jpt: $ 177 res =tab[premierx:dernierx,premiery:derniery,premierz:dernierz, *] 178 taille[1] eq jpi and taille[2] eq jpj and taille[3] eq nz and taille[4] eq jpt: $ 179 res =tab[premierx:dernierx,premiery:derniery,*, *] 180 taille[1] EQ nx and taille[2] eq ny and taille[3] eq nz and taille[4] eq jpt:res=tab 181 taille[1] EQ nx and taille[2] eq ny and taille[3] eq jpk and taille[4] eq jpt: $ 182 res=tab[*, *, premierz:dernierz, *] 183 else:return, report('Probleme d''adequation entre les tailles du domaine nx*ny*nz*jpt ' $ 184 +strtrim(nx, 1)+'*'+strtrim(ny, 1)+'*'+strtrim(nz, 1)+'*' $ 185 +strtrim(jpt, 1)+' et du tableau '+strtrim(taille[1], 1)+'*' $ 186 +strtrim(taille[2], 1)+'*'+strtrim(taille[3], 1)+'*' $ 187 +strtrim(taille[4], 1)) 188 endcase 189 divi = finite(res) 190 divi = total(temporary(divi),4, nan = keyword_set(nan)) 191 notanum = where(divi EQ 0) 192 if keyword_set(integration) then begin 193 res = total(res,4, nan = nan) 194 ENDIF ELSE BEGIN 195 if keyword_set(nan) then begin 196 res = total(res,4, /nan)/(1 > divi) 197 ENDIF ELSE res = total(res,4)/(1.*taille[4]) 198 ENDELSE 214 taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpk and taille[4] eq jpt: $ 215 res = tab[firstx:lastx, firsty:lasty, firstz:lastz, *] 216 taille[1] eq jpi and taille[2] eq jpj and taille[3] eq nz and taille[4] eq jpt: $ 217 res = tab[firstx:lastx, firsty:lasty, *, *] 218 taille[1] EQ nx and taille[2] eq ny and taille[3] eq nz and taille[4] eq jpt:res = tab 219 taille[1] EQ nx and taille[2] eq ny and taille[3] eq jpk and taille[4] eq jpt: $ 220 res = tab[*, *, firstz:lastz, *] 221 else:BEGIN 222 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 223 return, report('Probleme d''adequation entre les tailles du domaine nx*ny*nz*jpt ' $ 224 +strtrim(nx, 1)+'*'+strtrim(ny, 1)+'*'+strtrim(nz, 1)+'*' $ 225 +strtrim(jpt, 1)+' et du tableau '+strtrim(taille[1], 1)+'*' $ 226 +strtrim(taille[2], 1)+'*'+strtrim(taille[3], 1)+'*' $ 227 +strtrim(taille[4], 1)) 228 END 229 endcase 230 if keyword_set(integration) then begin 231 res = total(res, 4, nan = nan) 232 ENDIF ELSE BEGIN 233 if keyword_set(nan) then begin 234 divi = finite(res) 235 divi = total(temporary(divi), 4) 236 notanum = where(divi EQ 0) 237 res = total(res, 4, /nan)/(1 > divi) 238 if notanum[0] NE -1 then res[temporary(notanum)] = !values.f_nan 239 ENDIF ELSE res = total(res, 4)/(1.*taille[4]) 199 240 ENDELSE 200 if notanum[0] NE -1 then res[temporary(notanum)] = !values.f_nan 201 return, moyenne(temporary(res), direc, BOITE = boite, NAN = nan, INTEGRATION = integration, _extra = ex) 202 endif 241 ENDELSE 242 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 243 return, moyenne(temporary(res), direc, BOXZOOM = boxzoom, NAN = nan, INTEGRATION = integration, NODOMDEF = nodomdef, WDEPTH = wdepth, _extra = ex) 244 ENDIF ELSE res = tab 245 IF jpt EQ 1 THEN BEGIN 246 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 247 return, moyenne(reform(res, /over), direc, BOXZOOM = boxzoom, NAN = nan, INTEGRATION = integration, NODOMDEF = nodomdef, WDEPTH = wdepth, _extra = ex) 248 END 203 249 ;------------------------------------------------------------ 204 250 ;------------------------------------------------------------ … … 206 252 ;------------------------------------------------------------ 207 253 ;------------------------------------------------------------ 208 254 if (dim eq '3d') then begin 209 255 ;--------------------------------------------------------------- 210 256 ; II.1) verification de la coherence de la taille du tableau a … … 215 261 ; (jpi,jpj,jpt) soit celle du domaine reduit (nx,ny,jpt) 216 262 ;--------------------------------------------------------------- 217 case 1 of 218 taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpt: $ 219 res = tab[premierx:premierx+nx-1 $ 220 ,premiery:premiery+ny-1, *] 221 taille[1] EQ nx and taille[2] eq ny and taille[3] eq jpt:res=tab 222 else:return, report('Probleme d''adequation entre les tailles du domaine nx*ny*jpt '+strtrim(nx, 1)+'*'+strtrim(ny, 1)+'*'+strtrim(jpt, 1)+' et du tableau '+strtrim(taille[1], 1)+'*'+strtrim(taille[2], 1)+'*'+strtrim(taille[3], 1)) 223 endcase 224 if keyword_set(nan) NE 0 then BEGIN 225 if nan NE 1 then BEGIN ; si nan n''est pas !values.f_nan 263 case 1 of 264 taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpt: $ 265 res = tab[firstx:firstx+nx-1 $ 266 , firsty:firsty+ny-1, *] 267 taille[1] EQ nx and taille[2] eq ny and taille[3] eq jpt:res = tab 268 else:BEGIN 269 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 270 return, report('Probleme d''adequation entre les tailles du domaine nx*ny*jpt '+strtrim(nx, 1)+'*'+strtrim(ny, 1)+'*'+strtrim(jpt, 1)+' et du tableau '+strtrim(taille[1], 1)+'*'+strtrim(taille[2], 1)+'*'+strtrim(taille[3], 1)) 271 enD 272 endcase 273 if keyword_set(nan) NE 0 then BEGIN 274 if nan NE 1 then BEGIN ; si nan n''est pas !values.f_nan 226 275 ; on le met a !values.f_nan 227 if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $ 228 ELSE notanumber = where(abs(res) GT abs(nan)/10.) 229 if notanumber[0] NE -1 then res[temporary(notanumber)] = !values.f_nan 230 ENDIF 276 if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $ 277 ELSE notanumber = where(abs(res) GT abs(nan)/10.) 278 if notanumber[0] NE -1 then res[temporary(notanumber)] = !values.f_nan 231 279 ENDIF 280 ENDIF 232 281 ;--------------------------------------------------------------- 233 282 ; rq IL FAUT FAIRE ATTENTION AUX CAS OU LA DIM A MOYENNER = 1, ET … … 235 284 ; PEUVENT SEMBLER INUTILE AU DEPART 236 285 ;--------------------------------------------------------------- 237 if nx EQ 1 OR ny EQ 1 OR jpt EQ 1 then res=reform(res,nx,ny,jpt, /over) 238 if nx EQ 1 OR ny EQ 1 then BEGIN 239 mask = reform(mask, nx, ny, nz, /over) 240 e1 = reform(e1, nx, ny, /over) 241 e2 = reform(e2, nx, ny, /over) 242 endif 286 if nx EQ 1 OR ny EQ 1 then BEGIN 287 res = reform(res, nx, ny, jpt, /over) 288 e1 = reform(e1, nx, ny, /over) 289 e2 = reform(e2, nx, ny, /over) 290 endif 291 if nx EQ 1 OR ny EQ 1 OR nz EQ 1 THEN $ 292 mask = reform(mask, nx, ny, nz, /over) 243 293 ;--------------------------------------------------------------- 244 294 ; II.3) differents types de moyennes 245 295 ;--------------------------------------------------------------- 246 if keyword_set(nan) NE 0 then begin 247 msknan = replicate(1., nx, ny, jpt) 248 notanumber = where(finite(res) EQ 0) 249 if notanumber[0] NE -1 then msknan[temporary(notanumber)] = !values.f_nan 250 ENDIF ELSE msknan = 1 251 if nz eq jpk then mask = mask[*, *,niveau-1] ELSE mask = mask[*, *,nz-1] 252 case 1 of 253 (dirx eq 1) and (diry eq 0) : begin 254 e = temporary(e1)*temporary(mask) 255 if nx EQ 1 OR ny EQ 1 then e = reform(e, nx, ny, /over) 256 echelle = temporary(e[*])#replicate(1, jpt) 257 echelle = reform(echelle, nx, ny, jpt, /over) 258 if keyword_set(integration) then divi=1 $ 259 ELSE divi = total(echelle*msknan,1, nan = nan) 260 res=total(temporary(res)*echelle,1, nan = nan)/(divi > 1.) 261 if keyword_set(nan) then begin 262 testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0) 263 testnan = total(testnan,1) 264 endif 265 end 266 (dirx eq 0) and (diry eq 1) : begin 267 e = temporary(e2)*temporary(mask) 268 if nx EQ 1 OR ny EQ 1 then e = reform(e, nx, ny, /over) 269 echelle = temporary(e[*])#replicate(1, jpt) 270 echelle = reform(echelle, nx, ny, jpt, /over) 271 if keyword_set(integration) then divi=1 $ 272 ELSE divi = total(echelle*msknan,2, nan = nan) 273 res=total(temporary(res)*echelle,2, nan = nan)/(divi > 1.) 274 if keyword_set(nan) then begin 275 testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0) 276 testnan = total(testnan,2) 277 endif 278 end 279 (dirx eq 1) and (diry eq 1) : begin 280 echelle=(temporary(e1)*temporary(e2)*temporary(mask))[*]#replicate(1,jpt) 281 echelle=reform(echelle,nx,ny,jpt, /over) 282 if keyword_set(integration) then divi=1 $ 283 ELSE divi = total(temporary(total(echelle*msknan,1, nan = nan)),1, nan = nan) 284 res = total(temporary(total(temporary(res)*echelle,1, nan=nan)),1, nan=nan)/(divi > 1.) 285 if keyword_set(nan) then begin 286 testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0) 287 testnan = total(total(testnan,1),1) 288 endif 289 end 290 endcase 291 endif 296 if keyword_set(nan) NE 0 then msknan = finite(res) ELSE msknan = -1 297 mask = mask[*, *, 0] 298 case 1 of 299 (dirx eq 1) and (diry eq 0) : begin 300 e = temporary(e1)*temporary(mask) 301 echelle = (temporary(e))[*]#replicate(1, jpt) 302 echelle = reform(echelle, nx, ny, jpt, /over) 303 if keyword_set(integration) then divi = 1 ELSE BEGIN 304 IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 1) $ 305 ELSE divi = total(echelle, 1) 306 ENDELSE 307 res = total(temporary(res)*echelle, 1, nan = nan)/(divi > 1.) 308 if msknan[0] NE -1 then BEGIN 309 echelle = temporary(echelle) NE 0 310 testnan = temporary(msknan)*echelle 311 testnan = total(temporary(testnan), 1) $ 312 +(total(temporary(echelle), 1) EQ 0) 313 endif 314 end 315 (dirx eq 0) and (diry eq 1) : begin 316 e = temporary(e2)*temporary(mask) 317 if nx EQ 1 OR ny EQ 1 then e = reform(e, nx, ny, /over) 318 echelle = (temporary(e))[*]#replicate(1, jpt) 319 echelle = reform(echelle, nx, ny, jpt, /over) 320 if keyword_set(integration) then divi = 1 ELSE BEGIN 321 IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 2) $ 322 ELSE divi = total(echelle, 2) 323 ENDELSE 324 res = total(temporary(res)*echelle, 2, nan = nan)/(divi > 1.) 325 if msknan[0] NE -1 then begin 326 echelle = temporary(echelle) NE 0 327 testnan = temporary(msknan)*echelle 328 testnan = total(temporary(testnan), 2) $ 329 +(total(temporary(echelle), 2) EQ 0) 330 endif 331 end 332 (dirx eq 1) and (diry eq 1) : begin 333 echelle = (temporary(e1)*temporary(e2)*temporary(mask))[*]#replicate(1, jpt) 334 echelle = reform(echelle, nx, ny, jpt, /over) 335 if keyword_set(integration) then divi = 1 ELSE BEGIN 336 IF msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 1), 1) $ 337 ELSE divi = total(total(echelle, 1), 1) 338 ENDELSE 339 res = total(temporary(total(temporary(res)*echelle, 1, nan = nan)), 1, nan = nan)/(divi > 1.) 340 if msknan[0] NE -1 then begin 341 echelle = temporary(echelle) NE 0 342 testnan = temporary(msknan)*echelle 343 testnan = total(total(temporary(testnan), 1), 1) $ 344 +(total(total(temporary(echelle), 1), 1) EQ 0) 345 endif 346 end 347 endcase 348 endif 292 349 ;------------------------------------------------------------ 293 350 ;------------------------------------------------------------ 294 351 ; III) Cas serie tableaux 3d (tab4d) 295 352 ;------------------------------------------------------------ 296 353 if (dim eq '4d') then begin 297 354 ;--------------------------------------------------------------- 298 355 ; III.1) verification de la coherence de la taille du tableau a … … 303 360 ; (jpi,jpj,jpk,jpt) soit celle du domaine reduit (nx,ny,ny,jpt) 304 361 ;--------------------------------------------------------------- 305 case 1 of 306 taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpk and taille[4] eq jpt: $ 307 res =tab[premierx:dernierx,premiery:derniery,premierz:dernierz, *] 308 taille[1] eq jpi and taille[2] eq jpj and taille[3] eq nz and taille[4] eq jpt: $ 309 res =tab[premierx:dernierx,premiery:derniery,*, *] 310 taille[1] EQ nx and taille[2] eq ny and taille[3] eq nz and taille[4] eq jpt:res=tab 311 taille[1] EQ nx and taille[2] eq ny and taille[3] eq jpk and taille[4] eq jpt: $ 312 res=tab[*, *, premierz:dernierz, *] 313 else:return, report('Probleme d''adequation entre les tailles du domaine nx*ny*nz*jpt ' $ 314 +strtrim(nx, 1)+'*'+strtrim(ny, 1)+'*'+strtrim(nz, 1)+'*' $ 315 +strtrim(jpt, 1)+' et du tableau '+strtrim(taille[1], 1)+'*' $ 316 +strtrim(taille[2], 1)+'*'+strtrim(taille[3], 1)+'*' $ 317 +strtrim(taille[4], 1)) 318 endcase 319 if nx EQ 1 OR ny EQ 1 OR nz EQ 1 OR jpt EQ 1 then res=reform(res,nx,ny,nz,jpt, /over) 320 if keyword_set(nan) NE 0 then BEGIN 321 if nan NE 1 then BEGIN ; si nan n''est pas !values.f_nan 362 case 1 of 363 taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpk and taille[4] eq jpt: $ 364 res = tab[firstx:lastx, firsty:lasty, firstz:lastz, *] 365 taille[1] eq jpi and taille[2] eq jpj and taille[3] eq nz and taille[4] eq jpt: $ 366 res = tab[firstx:lastx, firsty:lasty, *, *] 367 taille[1] EQ nx and taille[2] eq ny and taille[3] eq nz and taille[4] eq jpt:res = tab 368 taille[1] EQ nx and taille[2] eq ny and taille[3] eq jpk and taille[4] eq jpt: $ 369 res = tab[*, *, firstz:lastz, *] 370 else:BEGIN 371 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 372 return, report('Probleme d''adequation entre les tailles du domaine nx*ny*nz*jpt ' $ 373 +strtrim(nx, 1)+'*'+strtrim(ny, 1)+'*'+strtrim(nz, 1)+'*' $ 374 +strtrim(jpt, 1)+' et du tableau '+strtrim(taille[1], 1)+'*' $ 375 +strtrim(taille[2], 1)+'*'+strtrim(taille[3], 1)+'*' $ 376 +strtrim(taille[4], 1)) 377 END 378 endcase 379 if nx EQ 1 OR ny EQ 1 OR nz EQ 1 OR jpt EQ 1 then res = reform(res, nx, ny, nz, jpt, /over) 380 if keyword_set(nan) NE 0 then BEGIN 381 if nan NE 1 then BEGIN ; si nan n''est pas !values.f_nan 322 382 ; on le met a !values.f_nan 323 if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $ 324 ELSE notanumber = where(abs(res) GT abs(nan)/10.) 325 if notanumber[0] NE -1 then res[temporary(notanumber)] = !values.f_nan 326 ENDIF 383 if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $ 384 ELSE notanumber = where(abs(res) GT abs(nan)/10.) 385 if notanumber[0] NE -1 then res[temporary(notanumber)] = !values.f_nan 327 386 ENDIF 387 ENDIF 328 388 ;--------------------------------------------------------------- 329 389 ; rq IL FAUT FAIRE ATTENTION AUX CAS OU LA DIM A MOYENNER = 1, ET … … 331 391 ; PEUVENT SEMBLER INUTILE AU DEPART 332 392 ;--------------------------------------------------------------- 333 if nx EQ 1 OR ny EQ 1 OR nz EQ 1 OR jpt EQ 1 then res=reform(res,nx,ny,nz,jpt, /over) 334 if nx EQ 1 OR ny EQ 1 OR nz EQ 1 then mask = reform(mask, nx, ny, nz, /over) 393 if nx EQ 1 OR ny EQ 1 OR nz EQ 1 then BEGIN 394 res = reform(res, nx, ny, nz, jpt, /over) 395 mask = reform(mask, nx, ny, nz, /over) 396 ENDIF 397 IF keyword_set(key_partialstep) THEN BEGIN 398 ; the top of the ocean floor is 399 IF vargrid EQ 'T' OR vargrid EQ 'W' THEN bottom = total(mask, 3) $ 400 ELSE bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3) 401 ; we suppress columns with only ocean or land 402 good = where(bottom NE 0 AND bottom NE nz) 403 ; the bottom of the ocean in 3D index is: 404 bottom = lindgen(nx*ny)+(temporary(bottom)-1L)*nx*ny 405 IF good[0] NE -1 THEN bottom = bottom[good] $ 406 ELSE bottom = -1 407 ENDIF ELSE bottom = -1 335 408 ;--------------------------------------------------------------- 336 409 ; III.2) differents types de moyennes 337 410 ;--------------------------------------------------------------- 338 IF keyword_set(nan) NE 0 then begin 339 msknan = replicate(1., nx, ny, nz, jpt) 340 notanumber = where(finite(res) EQ 0) 341 if notanumber[0] NE -1 then msknan[temporary(notanumber)] = !values.f_nan 342 ENDIF ELSE msknan = 1 343 case 1 of 344 (dirx eq 1) and (diry eq 0) and (dirz eq 0) : begin 345 e13 = temporary(e1[*])#replicate(1.,nz) 346 e13 = reform(e13,nx,ny,nz, /over) 347 echelle = (temporary(e13)*temporary(mask))[*]#replicate(1, jpt) 348 echelle = reform(echelle, nx, ny, nz, jpt, /over) 349 if keyword_set(integration) then divi=1 $ 350 ELSE begin 351 if keyword_set(nan) then divi = echelle*msknan ELSE divi = echelle 352 divi = total(temporary(divi),1, nan = nan) 353 endelse 354 res = temporary(res)*echelle 355 res = total(temporary(res),1, nan = nan)/(divi > 1) 356 if keyword_set(nan) then begin 357 testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0) 358 testnan = total(testnan,1) 359 endif 360 end 361 (dirx eq 0) and (diry eq 1) and (dirz eq 0) : begin 362 e23 = temporary(e2[*])#replicate(1.,nz) 363 e23 = reform(e23,nx,ny,nz, /over) 364 echelle = (temporary(e23)*temporary(mask))[*]#replicate(1, jpt) 365 echelle = reform(echelle, nx, ny, nz, jpt, /over) 366 if keyword_set(integration) then divi=1 $ 367 ELSE begin 368 if keyword_set(nan) then divi = echelle*msknan ELSE divi = echelle 369 divi = total(temporary(divi), 2, nan = nan) 370 endelse 371 res = total(temporary(res)*echelle,2, nan = nan)/(divi > 1) 372 if keyword_set(nan) then begin 373 testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0) 374 testnan = total(testnan,2) 375 endif 376 end 377 (dirx eq 0) and (diry eq 0) and (dirz eq 1) : begin 378 e33 = replicate(1,1.*nx*ny)#e3 379 e33 = reform(e33,nx,ny,nz, /over) 380 echelle =(temporary(e33)*temporary(mask))[*]#replicate(1, jpt) 381 echelle = reform(echelle, nx, ny, nz, jpt, /over) 382 if keyword_set(integration) then divi=1 $ 383 ELSE begin 384 if keyword_set(nan) then divi = echelle*msknan ELSE divi = echelle 385 divi = total(temporary(divi), 3, nan = nan) 386 endelse 387 res = total(temporary(res)*echelle,3, nan = nan)/(divi > 1) 388 if keyword_set(nan) then begin 389 testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0) 390 testnan = total(testnan,3) 391 endif 392 end 393 (dirx eq 1) and (diry eq 1) and (dirz eq 0) : begin 394 e13 = e1[*]#replicate(1.,nz) 395 e13 = reform(e13,nx,ny,nz, /over) 396 e23 = e2[*]#replicate(1.,nz) 397 e23 = reform(e23,nx,ny,nz, /over) 398 echelle = (temporary(e13)*temporary(e23)*temporary(mask))[*]#replicate(1, jpt) 399 echelle = reform(echelle,nx,ny,nz,jpt, /over) 400 if keyword_set(integration) then divi=1 $ 401 ELSE begin 402 if keyword_set(nan) then divi = echelle*msknan ELSE divi = echelle 403 divi = total(total(temporary(divi),1, nan = nan),1, nan = nan) 404 endelse 405 res =total(total(temporary(res)*echelle,1, nan = nan),1, nan = nan)/(divi > 1) 406 if keyword_set(nan) then begin 407 testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0) 408 testnan = total(total(testnan,1),1) 409 endif 410 end 411 (dirx eq 1) and (diry eq 0) and (dirz eq 1) : begin 412 e133 = e1[*]#e3 413 echelle = (temporary(e133[*])*temporary(mask[*]))#replicate(1, jpt) 414 echelle = reform(echelle,nx,ny,nz,jpt, /over) 415 if keyword_set(integration) then divi=1 $ 416 ELSE begin 417 if keyword_set(nan) then divi = echelle*msknan ELSE divi = echelle 418 divi = total(total(temporary(divi),1, nan = nan),2, nan = nan) 419 endelse 420 res = total(total(temporary(res)*echelle,1, nan = nan),2, nan = nan)/(divi > 1) 421 if keyword_set(nan) then begin 422 testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0) 423 testnan = total(total(testnan,1),2) 424 endif 425 end 426 (dirx eq 0) and (diry eq 1) and (dirz eq 1) : begin 427 e233 = e2[*]#e3 428 echelle = (temporary(e233[*])*temporary(mask[*]))#replicate(1, jpt) 429 echelle = reform(echelle,nx,ny,nz,jpt, /over) 430 if keyword_set(integration) then divi=1 $ 431 ELSE begin 432 if keyword_set(nan) then divi = echelle*msknan ELSE divi = echelle 433 divi = total(total(temporary(divi),2, nan = nan),2,nan=nan) 434 endelse 435 res =total(total(temporary(res)*echelle,2, nan = nan),2, nan = nan)/(divi > 1) 436 if keyword_set(nan) then begin 437 testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0) 438 testnan = total(total(testnan,2),2) 439 endif 440 end 441 (dirx eq 1) and (diry eq 1) and (dirz eq 1) : begin 442 e1233 = (e1[*]*e2[*])#e3 443 echelle=(temporary(e1233[*])*temporary(mask[*]))#replicate(1,jpt) 444 echelle=reform(echelle,nx,ny,nz,jpt, /over) 445 if keyword_set(integration) then divi=1 $ 446 ELSE begin 447 if keyword_set(nan) then divi = echelle*msknan ELSE divi = echelle 448 divi = total(total(total(temporary(divi), 1,nan=nan), 1,nan=nan), 1,nan=nan) 449 endelse 450 res =total(total(total(temporary(res)*echelle, 1,nan=nan), 1,nan=nan), 1,nan=nan)/(divi > 1) 451 if keyword_set(nan) then begin 452 testnan = finite(temporary(msknan))+(temporary(echelle) EQ 0) 453 testnan = total(total(total(testnan,1),1),1) 454 endif 455 end 456 endcase 457 endif 411 IF keyword_set(nan) NE 0 THEN msknan = finite(res) ELSE msknan = -1 412 case 1 of 413 (dirx eq 1) and (diry eq 0) and (dirz eq 0) : BEGIN 414 e13 = (temporary(e1))[*]#replicate(1., nz) 415 e13 = reform(e13, nx, ny, nz, /over) 416 echelle = (temporary(e13)*temporary(mask))[*]#replicate(1, jpt) 417 echelle = reform(echelle, nx, ny, nz, jpt, /over) 418 IF keyword_set(key_partialstep) AND bottom[0] NE -1 $ 419 AND nx NE 1 THEN BEGIN 420 IF msknan[0] EQ -1 THEN BEGIN 421 msknan = replicate(1b, nx, ny, nz, jpt) 422 nan = 1 423 ENDIF 424 bottom = bottom#replicate(1, jpt) $ ; 4D bottom! 425 + replicate(1, n_elements(bottom))#(nx*ny*nz*lindgen(jpt)) 426 msknan[bottom] = 0 427 res[temporary(bottom)] = !values.f_nan 428 ENDIF 429 if keyword_set(integration) then divi = 1 ELSE begin 430 IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 1) $ 431 ELSE divi = total(echelle, 1) 432 endelse 433 res = temporary(res)*echelle 434 res = total(temporary(res), 1, nan = nan)/(divi > 1) 435 if msknan[0] NE -1 then begin 436 echelle = temporary(echelle) NE 0 437 testnan = temporary(msknan)*echelle 438 testnan = total(temporary(testnan), 1) $ 439 +(total(temporary(echelle), 1) EQ 0) 440 endif 441 end 442 (dirx eq 0) and (diry eq 1) and (dirz eq 0) : begin 443 e23 = temporary(e2[*])#replicate(1., nz) 444 e23 = reform(e23, nx, ny, nz, /over) 445 echelle = (temporary(e23)*temporary(mask))[*]#replicate(1, jpt) 446 echelle = reform(echelle, nx, ny, nz, jpt, /over) 447 IF keyword_set(key_partialstep) AND bottom[0] NE -1 $ 448 AND ny NE 1 THEN BEGIN 449 IF msknan[0] EQ -1 THEN BEGIN 450 msknan = replicate(1b, nx, ny, nz) 451 nan = 1 452 endif 453 bottom = bottom#replicate(1, jpt) $ ; 4D bottom! 454 + replicate(1, n_elements(bottom))#(nx*ny*nz*lindgen(jpt)) 455 msknan[bottom] = 0 456 res[temporary(bottom)] = !values.f_nan 457 ENDIF 458 if keyword_set(integration) then divi = 1 ELSE begin 459 IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 2) $ 460 ELSE divi = total(echelle, 2) 461 endelse 462 res = total(temporary(res)*echelle, 2, nan = nan)/(divi > 1) 463 if msknan[0] NE -1 then begin 464 echelle = temporary(echelle) NE 0 465 testnan = temporary(msknan)*echelle 466 testnan = total(temporary(testnan), 2) $ 467 +(total(temporary(echelle), 2) EQ 0) 468 endif 469 end 470 (dirx eq 0) and (diry eq 0) and (dirz eq 1) : begin 471 e33 = replicate(1, 1.*nx*ny)#e3 472 e33 = reform(e33, nx, ny, nz, /over) 473 IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN 474 IF keyword_set(wdepth) THEN $ 475 e33[bottom] = (e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $ 476 ELSE e33[bottom] = (e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)] 477 ENDIF 478 echelle = (temporary(e33)*temporary(mask))[*]#replicate(1, jpt) 479 echelle = reform(echelle, nx, ny, nz, jpt, /over) 480 if keyword_set(integration) then divi = 1 ELSE begin 481 IF msknan[0] NE -1 THEN divi = total(echelle*msknan, 3) $ 482 ELSE divi = total(echelle, 3) 483 endelse 484 res = total(temporary(res)*echelle, 3, nan = nan)/(divi > 1) 485 if msknan[0] NE -1 then begin 486 echelle = temporary(echelle) NE 0 487 testnan = temporary(msknan)*echelle 488 testnan = total(temporary(testnan), 3) $ 489 +(total(temporary(echelle), 3) EQ 0) 490 endif 491 end 492 (dirx eq 1) and (diry eq 1) and (dirz eq 0) : begin 493 e13 = e1[*]#replicate(1., nz) 494 e13 = reform(e13, nx, ny, nz, /over) 495 e23 = e2[*]#replicate(1., nz) 496 e23 = reform(e23, nx, ny, nz, /over) 497 echelle = (temporary(e13)*temporary(e23)*temporary(mask))[*]#replicate(1, jpt) 498 echelle = reform(echelle, nx, ny, nz, jpt, /over) 499 IF keyword_set(key_partialstep) AND bottom[0] NE -1 $ 500 AND nx*ny NE 1 THEN BEGIN 501 IF msknan[0] EQ -1 THEN BEGIN 502 msknan = replicate(1b, nx, ny, nz) 503 nan = 1 504 endif 505 bottom = bottom#replicate(1, jpt) $ ; 4D bottom! 506 + replicate(1, n_elements(bottom))#(nx*ny*nz*lindgen(jpt)) 507 msknan[bottom] = 0 508 res[temporary(bottom)] = !values.f_nan 509 ENDIF 510 if keyword_set(integration) then divi = 1 ELSE begin 511 IF msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 1), 1) $ 512 ELSE divi = total(total(echelle, 1), 1) 513 endelse 514 res = total(total(temporary(res)*echelle, 1, nan = nan), 1, nan = nan)/(divi > 1) 515 if msknan[0] NE -1 then begin 516 echelle = temporary(echelle) NE 0 517 testnan = temporary(msknan)*echelle 518 testnan = total(total(temporary(testnan), 1), 1) $ 519 +(total(total(temporary(echelle), 1), 1) EQ 0) 520 endif 521 end 522 (dirx eq 1) and (diry eq 0) and (dirz eq 1) : begin 523 e133 = e1[*]#e3 524 IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN 525 IF keyword_set(wdepth) THEN $ 526 e133[bottom] = (e1*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $ 527 ELSE e133[bottom] = (e1*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)] 528 ENDIF 529 echelle = (temporary(e133[*])*temporary(mask[*]))#replicate(1, jpt) 530 echelle = reform(echelle, nx, ny, nz, jpt, /over) 531 if keyword_set(integration) then divi = 1 ELSE begin 532 IF msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 1), 2) $ 533 ELSE divi = total(total(echelle, 1), 2) 534 endelse 535 res = total(total(temporary(res)*echelle, 1, nan = nan), 2, nan = nan)/(divi > 1) 536 if msknan[0] NE -1 then begin 537 echelle = temporary(echelle) NE 0 538 testnan = temporary(msknan)*echelle 539 testnan = total(total(temporary(testnan), 1), 2) $ 540 +(total(total(temporary(echelle), 1), 2) EQ 0) 541 endif 542 end 543 (dirx eq 0) and (diry eq 1) and (dirz eq 1) : begin 544 e233 = e2[*]#e3 545 IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN 546 IF keyword_set(wdepth) THEN $ 547 e233[bottom] = (e2*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $ 548 ELSE e233[bottom] = (e2*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)] 549 ENDIF 550 echelle = (temporary(e233[*])*temporary(mask[*]))#replicate(1, jpt) 551 echelle = reform(echelle, nx, ny, nz, jpt, /over) 552 if keyword_set(integration) then divi = 1 ELSE begin 553 IF msknan[0] NE -1 THEN divi = total(total(echelle*msknan, 2), 2) $ 554 ELSE divi = total(total(echelle, 2), 2) 555 endelse 556 res = total(total(temporary(res)*echelle, 2, nan = nan), 2, nan = nan)/(divi > 1) 557 if msknan[0] NE -1 then begin 558 echelle = temporary(echelle) NE 0 559 testnan = temporary(msknan)*echelle 560 testnan = total(total(temporary(testnan), 2), 2) $ 561 +(total(total(temporary(echelle), 2), 2) EQ 0) 562 endif 563 end 564 (dirx eq 1) and (diry eq 1) and (dirz eq 1) : begin 565 e1233 = (e1[*]*e2[*])#e3 566 IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN 567 IF keyword_set(wdepth) THEN $ 568 e1233[bottom] = (e1*e2*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $ 569 ELSE e1233[bottom] = (e1*e2*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)] 570 ENDIF 571 echelle = (temporary(e1233[*])*temporary(mask[*]))#replicate(1, jpt) 572 echelle = reform(echelle, nx, ny, nz, jpt, /over) 573 if keyword_set(integration) then divi = 1 ELSE begin 574 IF msknan[0] NE -1 THEN divi = total(total(total(echelle*msknan, 1), 1), 1) $ 575 ELSE divi = total(total(total(echelle, 1), 1), 1) 576 endelse 577 res = total(total(total(temporary(res)*echelle, 1, nan = nan), 1, nan = nan), 1, nan = nan)/(divi > 1) 578 if msknan[0] NE -1 then begin 579 echelle = temporary(echelle) NE 0 580 testnan = temporary(msknan)*echelle 581 testnan = total(total(total(temporary(testnan), 1), 1), 1) $ 582 +(total(total(total(temporary(echelle), 1), 1), 1) EQ 0) 583 endif 584 end 585 endcase 586 endif 587 ;------------------------------------------------------------ 588 if dirt EQ 1 AND keyword_set(spatialfirst) then BEGIN 589 IF (reverse(size(res, /dimension)))[0] NE jpt THEN BEGIN 590 print, 'the last dimension of res is not equal to jpt: '+strtrim(jpt, 2) 591 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 592 return, -1 593 ENDIF 594 tdim = size(res, /n_dimensions) 595 if keyword_set(integration) then res = total(res, tdim, nan = nan) $ 596 ELSE BEGIN 597 if keyword_set(nan) then BEGIN 598 testnan = testnan < 1 599 testnan = total(temporary(testnan), tdim) 600 divi = testnan 601 ENDIF ELSE divi = jpt 602 res = total(res, tdim, nan = nan)/(1 > divi) 603 ENDELSE 604 ENDIF 458 605 ;------------------------------------------------------------ 459 606 ;------------------------------------------------------------ … … 463 610 ; IV.1) on masque les terres par une valeur a 1e+20 464 611 ;------------------------------------------------------------ 465 466 467 468 469 612 valmask = 1e+20 613 terre = where(divi EQ 0) 614 IF terre[0] NE -1 THEN BEGIN 615 res(temporary(terre)) = 1e+20 616 ENDIF 470 617 ;------------------------------------------------------------ 471 618 ; IV.2) on remplace, quand nan ne 1, !values.f_nan par nan 472 619 ;------------------------------------------------------------ 473 474 475 476 477 478 479 480 620 if keyword_set(nan) NE 0 then BEGIN 621 puttonan = where(temporary(testnan) EQ 0) 622 if puttonan[0] NE -1 then res[temporary(puttonan)] = !values.f_nan 623 if nan NE 1 then BEGIN 624 notanumber = where(finite(res) eq 0) 625 if notanumber[0] NE -1 then res[temporary(notanumber)] = nan 626 ENDIF 627 ENDIF 481 628 ;------------------------------------------------------------ 482 629 ; IV.3) on se remplace ds le sous domaine qui etait definit a l''entree de 483 630 ; moyenne 484 631 ;------------------------------------------------------------ 485 if keyword_set(boite) AND NOT keyword_set(nodomdef) then domdef, oldboite,GRILLE=vargrid486 ;------------------------------------------------------------ 487 488 632 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grmoyenne.dat' 633 ;------------------------------------------------------------ 634 if keyword_set(key_performance) THEN print, 'temps grossemoyenne', systime(1)-tempsun 635 return, res 489 636 ;------------------------------------------------------------ 490 637 ;------------------------------------------------------------ -
trunk/ToBeReviewed/CALCULS/hdyn.pro
r23 r25 101 101 case 1 of 102 102 tailles[1] eq jpi and tailles[2] eq jpj: BEGIN 103 sn = tabsn[ premierxt:dernierxt, premieryt:dernieryt, *]104 tn = tabtn[ premierxt:dernierxt, premieryt:dernieryt, *]103 sn = tabsn[firstxt:lastxt, firstyt:lastyt, *] 104 tn = tabtn[firstxt:lastxt, firstyt:lastyt, *] 105 105 end 106 106 tailles[1] eq nx and tailles[2] eq ny:BEGIN … … 116 116 e33d = replicate(1, nx*ny)#e3t 117 117 e33d = reform(e33d, nx, ny, jpk, /over) 118 terre = where(tmask[ premierxt:dernierxt, premieryt:dernieryt, *] EQ 0)118 terre = where(tmask[firstxt:lastxt, firstyt:lastyt, *] EQ 0) 119 119 if terre[0] NE -1 then vol[terre] = !values.f_nan 120 120 case level of … … 128 128 case 1 of 129 129 tailles[1] eq jpi and tailles[2] eq jpj AND tailles[4] EQ jpt: BEGIN 130 sn = tabsn[ premierxt:dernierxt, premieryt:dernieryt, *, *]131 tn = tabtn[ premierxt:dernierxt, premieryt:dernieryt, *, *]130 sn = tabsn[firstxt:lastxt, firstyt:lastyt, *, *] 131 tn = tabtn[firstxt:lastxt, firstyt:lastyt, *, *] 132 132 end 133 133 tailles[1] eq nx and tailles[2] eq ny AND tailles[4] EQ jpt:BEGIN … … 144 144 e33d = e33d[*]#replicate(1, jpt) 145 145 e33d = reform(e33d, nx, ny, jpk, jpt, /over) 146 mask = tmask[ premierxt:dernierxt, premieryt:dernieryt, *]146 mask = tmask[firstxt:lastxt, firstyt:lastyt, *] 147 147 mask = mask[*]#replicate(1, jpt) 148 148 terre = where(mask EQ 0) -
trunk/ToBeReviewed/CALCULS/level2depth.pro
r23 r25 45 45 ;------------------------------------------------------------ 46 46 niveaux = litchamp(tab) 47 grille,mask, glam,gphi,gdep,nx,ny,nz,premierx,premiery,premierz,dernierx,derniery,dernierz47 grille,mask, -1, -1,gdep,nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz 48 48 ;--------------------------------------------------------------- 49 49 ; verification de la coherence entre la taille du tableau et le domaine definit par domdef … … 52 52 if taille[0] NE 2 then return, report('le champ en entree doit contenir un tableau 2d') 53 53 case 1 of 54 taille[1] eq jpi and taille[2] eq jpj:niveaux=niveaux[ premierx:dernierx, premiery:derniery]54 taille[1] eq jpi and taille[2] eq jpj:niveaux=niveaux[firstx:lastx, firsty:lasty] 55 55 taille[1] eq nx and taille[2] eq ny: 56 56 else:return, report('Probleme d''adequation entre les tailles du domaine et celle du champ.') -
trunk/ToBeReviewed/CALCULS/level2mask.pro
r23 r25 6 6 ; 7 7 ; PURPOSE: permet de passer d''un tableau 2d de niveau seuil au 8 ; tableau 3d de mask avec des 1 ds les niveaux au dessus (et sur)du9 ; niveau seuil et des 0 en dessous .8 ; tableau 3d de mask avec des 1 ds les niveaux au dessus du 9 ; niveau seuil et des 0 en dessous (et sur). 10 10 ; 11 11 ; CATEGORY: SANS BOUCLE … … 32 32 ; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr) 33 33 ; 17/6/1999 34 ; Setp 2004: boundary level have 0 values and not 1 (as it was 35 ; explained before in the header). see: 36 ; print, array_equal(niveau, total(level2mask(niveau),3)) 37 ; 34 38 ;- 35 39 ;------------------------------------------------------------ … … 45 49 ;------------------------------------------------------------ 46 50 niveaux = litchamp(tab) 47 grille,maskterre, glam,gphi,gdep,nx,ny,nz,premierx,premiery,premierz,dernierx,derniery,dernierz51 grille,maskterre, -1, -1, -1,nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz 48 52 ;--------------------------------------------------------------- 49 53 ; verification de la coherence entre la taille du tableau et le domaine definit par domdef … … 53 57 if taille[0] NE 2 then return, report('le champ en entree doit contenir un tableau 2d') 54 58 case 1 of 55 taille[1] eq jpi and taille[2] eq jpj:niveaux=niveaux[ premierx:dernierx, premiery:derniery]59 taille[1] eq jpi and taille[2] eq jpj:niveaux=niveaux[firstx:lastx, firsty:lasty] 56 60 taille[1] eq nx and taille[2] eq ny: 57 61 else:return, report('Probleme d''adequation entre les tailles du domaine et celle du champ.') -
trunk/ToBeReviewed/CALCULS/moyenne.pro
r23 r25 11 11 ; CATEGORY: 12 12 ; 13 ; CALLING SEQUENCE: result = moyenne(tab,'direc',BO ITE=boite)13 ; CALLING SEQUENCE: result = moyenne(tab,'direc',BOXZOOM=boxzoom) 14 14 ; 15 15 ; INPUTS: tab = 2 or 3d field … … 18 18 ; KEYWORD PARAMETERS: 19 19 ; 20 ; BO ITE= [xmin,xmax,ymin,ymax (,(zmin,)zmax)] pour plus20 ; BOXZOOM = [xmin,xmax,ymin,ymax (,(zmin,)zmax)] pour plus 21 21 ; de detail cf domdef. 22 ; bo itepeut prendre 5 formes:23 ; [ prof2], [prof1, prof2],[lon1, lon2, lat1, lat2],24 ; [lon1, lon2, lat1, lat2, prof2],[lon1, lon2, lat1,25 ; lat2, prof1,prof2]22 ; boxzoom peut prendre 5 formes: 23 ; [vert2], [vert1, vert2],[lon1, lon2, lat1, lat2], 24 ; [lon1, lon2, lat1, lat2, vert2],[lon1, lon2, lat1, 25 ; lat2, vert1,vert2] 26 26 ; 27 27 ; NAN: not a number, a activer si l''on peut faire veut … … 38 38 ; 39 39 ; NODOMDEF: activer si l''on ne veut pas passer ds 40 ; domdef bien que le mot cle bo itesoit present (comme40 ; domdef bien que le mot cle boxzoom soit present (comme 41 41 ; c''est le cas qd moyenne est appelee via checkfield) 42 42 ; 43 43 ; INTEGRATION: pour faire une integrale plutot qu''une moyenne 44 ; 45 ; /WDEPTH: to specify that the field is at W depth instad of T 46 ; depth (automatically activated if vargrid eq 'W') 44 47 ; 45 48 ; OUTPUTS: result:un tableau … … 88 91 ;------------------------------------------------------------ 89 92 ;------------------------------------------------------------ 90 function moyenne ,tab,direc,BOITE=boite, INTEGRATION=integration $ 91 , NAN = nan, NODOMDEF = nodomdef, _extra = ex 92 @common 93 tempsun = systime(1) ; pour key_performance 93 function moyenne, tab, direc, BOXZOOM = boxzoom, INTEGRATION = integration $ 94 , NAN = nan, NODOMDEF = nodomdef, WDEPTH = wdepth $ 95 , _extra = ex 96 ;--------------------------------------------------------- 97 @cm_4mesh 98 @cm_4data 99 @cm_4cal 100 IF NOT keyword_set(key_forgetold) THEN BEGIN 101 @updatenew 102 @updatekwd 103 ENDIF 104 ;--------------------------------------------------------- 105 tempsun = systime(1) ; pour key_performance 94 106 ;------------------------------------------------------------ 95 107 ;I) preliminaires 96 108 ;------------------------------------------------------------ 97 98 99 100 109 dirt = 0 110 dirx = 0 111 diry = 0 112 dirz = 0 101 113 ;------------------------------------------------------------ 102 114 ; I.1) direction(s) suivants lesquelles on integre 103 115 ;------------------------------------------------------------ 104 if ( strpos(direc,'t') ge 0 ) then dirt = 1 105 if ( strpos(direc,'x') ge 0 ) then dirx = 1 106 if ( strpos(direc,'y') ge 0 ) then diry = 1 107 if ( strpos(direc,'z') ge 0 ) then dirz = 1 108 if (dirx eq 0 and diry eq 0 and dirz eq 0) then BEGIN 109 IF keyword_set(bavard) then $ 110 IF dirt NE 1 THEN ras = report('MOYENNE: aucune valeur de direc ne convient, le champ reste inchange') 111 return, tab 112 ENDIF 116 if ( strpos(direc, 't') ge 0 ) then dirt = 1 117 if ( strpos(direc, 'x') ge 0 ) then dirx = 1 118 if ( strpos(direc, 'y') ge 0 ) then diry = 1 119 if ( strpos(direc, 'z') ge 0 ) then dirz = 1 120 if (dirx eq 0 and diry eq 0 and dirz eq 0) then return, tab 113 121 ;------------------------------------------------------------ 114 122 ; I.2) verification de la taille du tableau d'entree 115 123 ;------------------------------------------------------------ 116 117 118 119 120 121 122 123 124 125 126 127 128 124 taille = size(tab) 125 case 1 of 126 taille[0] eq 1 :dim = '1d' 127 taille[0] eq 2 :BEGIN 128 dim = '2d' 129 if dirx eq 0 and diry eq 0 then return, tab 130 END 131 taille[0] eq 3 :BEGIN 132 dim = '3d' 133 if dirx eq 0 and diry eq 0 and dirz eq 0 then return, tab 134 END 135 else : return, report('Le tableau d''entree doit etre a 2 ou 3 dimensions s''il contient une dim temporelle utiliser grossemoyenne') 136 endcase 129 137 ;------------------------------------------------------------ 130 138 ; I.3) obtention des facteurs d''echelles et du masque sur le sous 131 139 ; domaine concerne par la moyenne 132 ; redefinition du domaine ajuste a bo ite(a 6 elements)140 ; redefinition du domaine ajuste a boxzoom (a 6 elements) 133 141 ; ceci va nous permetre de faire les calcules que sur le sous domaine 134 142 ; comcerne par la moyenne. domdef, suivit de grille nous donne tous 135 143 ; les tableaux de la grille sur le sous domaine 136 144 ;------------------------------------------------------------ 137 if keyword_set(boite) then BEGIN 138 Case 1 Of 139 N_Elements(Boite) Eq 1:bte=[lon1, lon2, lat1, lat2, 0.,boite[0]] 140 N_Elements(Boite) Eq 2:bte=[lon1, lon2, lat1, lat2, boite[0],boite[1]] 141 N_Elements(Boite) Eq 4:bte=[Boite, prof1, prof2] 142 N_Elements(Boite) Eq 5:bte=[Boite[0:3], 0, Boite[4]] 143 N_Elements(Boite) Eq 6:bte=Boite 144 Else: return, report('Mauvaise Definition de Boite') 145 endcase 146 oldboite = [lon1, lon2, lat1, lat2, prof1, prof2] 147 if NOT keyword_set(nodomdef) then domdef, bte,GRILLE=vargrid, _extra = ex 148 ENDIF 145 if keyword_set(boxzoom) then BEGIN 146 Case 1 Of 147 N_Elements(Boxzoom) Eq 1:bte = [lon1, lon2, lat1, lat2, 0., boxzoom[0]] 148 N_Elements(Boxzoom) Eq 2:bte = [lon1, lon2, lat1, lat2, boxzoom[0], boxzoom[1]] 149 N_Elements(Boxzoom) Eq 4:bte = [Boxzoom, vert1, vert2] 150 N_Elements(Boxzoom) Eq 5:bte = [Boxzoom[0:3], 0, Boxzoom[4]] 151 N_Elements(Boxzoom) Eq 6:bte = Boxzoom 152 Else: return, report('Mauvaise Definition de Boxzoom') 153 endcase 154 if NOT keyword_set(nodomdef) then BEGIN 155 savedbox = 1b 156 saveboxparam, 'boxparam4moyenne.dat' 157 domdef, bte, GRIDTYPE = vargrid, _extra = ex 158 ENDIF 159 ENDIF 149 160 ;--------------------------------------------------------------- 150 161 ; attribution du mask et des tableaux de longitude et latitude... 151 162 ;--------------------------------------------------------------- 152 grille, mask, glam, gphi, gdep, nx, ny,nz,premierx,premiery,premierz,dernierx, derniery, dernierz,e1,e2,e3 163 IF vargrid EQ 'W' THEN wdepth = 1 164 grille, mask, glam, gphi, gdep, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz, e1, e2, e3, WDEPTH = wdepth 153 165 ;------------------------------------------------------------ 154 166 ;------------------------------------------------------------ … … 156 168 ;------------------------------------------------------------ 157 169 ;------------------------------------------------------------ 158 if dim EQ '1d' then BEGIN 159 if n_elements(tab) NE nx*ny AND n_elements(tab) NE nx*ny*nz then $ 160 return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 161 case 1 of 162 nx EQ 1 AND ny EQ 1:BEGIN ;vecteur suivant z 163 case n_elements(tab) of 164 jpk:res = tab[premierz:dernierz] 165 nz:res = tab 166 ELSE:return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 167 ENDCASE 168 if dirz EQ 1 then BEGIN 169 dim = '3d' 170 taille = size(reform(res, nx, ny, nz)) 171 ENDIF ELSE return, res 172 END 173 ny EQ 1:BEGIN ;vecteur suivant x 174 case n_elements(tab) of 175 jpi:res = tab[premierx:dernierx] 176 nx:res = tab 177 ELSE:return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 178 ENDCASE 179 if dirx EQ 1 then BEGIN 180 dim = '2d' 181 taille = size(reform(res, nx, ny)) 182 ENDIF ELSE return, res 183 END 184 nx EQ 1:BEGIN ;vecteur suivant y 185 case n_elements(tab) of 186 jpj:res = tab[premiery:derniery] 187 ny:res = tab 188 ELSE:return, report('Probleme d''adequation entre les tailles du domaine et de la boite.') 189 ENDCASE 190 if diry EQ 1 then BEGIN 191 dim = '2d' 192 taille = size(reform(res, nx, ny)) 193 ENDIF ELSE return, res 194 END 195 endcase 196 endif 170 if dim EQ '1d' then BEGIN 171 if n_elements(tab) NE nx*ny AND n_elements(tab) NE nx*ny*nz then BEGIN 172 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat' 173 return, report('Probleme d''adequation entre les tailles du domaine et de la boxzoom.') 174 ENDIF 175 case 1 of 176 nx EQ 1 AND ny EQ 1:BEGIN ;vecteur suivant z 177 case n_elements(tab) of 178 jpk:res = tab[firstz:lastz] 179 nz:res = tab 180 ELSE:BEGIN 181 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat' 182 return, report('Probleme d''adequation entre les tailles du domaine et de la boxzoom.') 183 END 184 ENDCASE 185 if dirz EQ 1 then BEGIN 186 dim = '3d' 187 taille = size(reform(res, nx, ny, nz)) 188 ENDIF ELSE BEGIN 189 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat' 190 return, res 191 ENDELSE 192 END 193 ny EQ 1:BEGIN ;vecteur suivant x 194 case n_elements(tab) of 195 jpi:res = tab[firstx:lastx] 196 nx:res = tab 197 ELSE:BEGIN 198 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat' 199 return, report('Probleme d''adequation entre les tailles du domaine et de la boxzoom.') 200 END 201 ENDCASE 202 if dirx EQ 1 then BEGIN 203 dim = '2d' 204 taille = size(reform(res, nx, ny)) 205 ENDIF ELSE BEGIN 206 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat' 207 return, res 208 ENDELSE 209 END 210 nx EQ 1:BEGIN ;vecteur suivant y 211 case n_elements(tab) of 212 jpj:res = tab[firsty:lasty] 213 ny:res = tab 214 ELSE:BEGIN 215 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat' 216 return, report('Probleme d''adequation entre les tailles du domaine et de la boxzoom.') 217 END 218 ENDCASE 219 if diry EQ 1 then BEGIN 220 dim = '2d' 221 taille = size(reform(res, nx, ny)) 222 ENDIF ELSE BEGIN 223 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat' 224 return, res 225 ENDELSE 226 END 227 endcase 228 endif 197 229 ;------------------------------------------------------------ 198 230 ;------------------------------------------------------------ … … 200 232 ;------------------------------------------------------------ 201 233 ;------------------------------------------------------------ 202 234 if (dim eq '2d') then begin 203 235 ;--------------------------------------------------------------- 204 236 ; II.1) verification de la coherence de la taille du tableau a … … 209 241 ; (jpi,jpj) soit celle du domaine reduit (nx,ny) 210 242 ;--------------------------------------------------------------- 211 case 1 of 212 taille[1] eq jpi and taille[2] eq jpj: $ 213 res = tab[premierx:dernierx, premiery:derniery] 214 taille[1] eq nx and taille[2] eq ny:res = tab 215 else:return, report('Probleme d''adequation entre les tailles du domaine nx*ny '+strtrim(nx, 1)+'*'+strtrim(ny, 1)+' et du tableau '+strtrim(taille[1], 1)+'*'+strtrim(taille[2], 1)) 216 ENDCASE 217 if keyword_set(nan) NE 0 then BEGIN 218 if nan NE 1 then BEGIN ; si nan n''est pas !values.f_nan 243 case 1 of 244 taille[1] eq jpi and taille[2] eq jpj: $ 245 res = tab[firstx:lastx, firsty:lasty] 246 taille[1] eq nx and taille[2] eq ny:res = tab 247 else:BEGIN 248 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat' 249 return, report('Probleme d''adequation entre les tailles du domaine nx*ny '+strtrim(nx, 1)+'*'+strtrim(ny, 1)+' et du tableau '+strtrim(taille[1], 1)+'*'+strtrim(taille[2], 1)) 250 END 251 ENDCASE 252 if keyword_set(nan) NE 0 then BEGIN 253 if nan NE 1 then BEGIN ; si nan n''est pas !values.f_nan 219 254 ; on le met a !values.f_nan 220 if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $ 221 ELSE notanumber = where(abs(res) GT abs(nan)/10.) 222 if notanumber[0] NE -1 then res[temporary(notanumber)] = !values.f_nan 223 ENDIF 255 if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $ 256 ELSE notanumber = where(abs(res) GT abs(nan)/10.) 257 if notanumber[0] NE -1 then res[temporary(notanumber)] = !values.f_nan 224 258 ENDIF 259 ENDIF 225 260 ;--------------------------------------------------------------- 226 261 ; rq IL FAUT FAIRE ATTENTION AUX CAS OU LA DIM A MOYENNER = 1, ET … … 228 263 ; PEUVENT SEMBLER INUTILE AU DEPART 229 264 ;--------------------------------------------------------------- 230 if nx EQ 1 OR ny EQ 1 then BEGIN 231 res = reform(res, nx, ny, /over) 232 mask = reform(mask, nx, ny, nz, /over) 233 e1 = reform(e1, nx, ny, /over) 234 e2 = reform(e2, nx, ny, /over) 235 endif 265 if nx EQ 1 OR ny EQ 1 then BEGIN 266 res = reform(res, nx, ny, /over) 267 e1 = reform(e1, nx, ny, /over) 268 e2 = reform(e2, nx, ny, /over) 269 endif 270 if nx EQ 1 OR ny EQ 1 OR nz EQ 1 THEN $ 271 mask = reform(mask, nx, ny, nz, /over) 236 272 ;--------------------------------------------------------------- 237 273 ; II.3) differents types de moyennes 238 274 ;--------------------------------------------------------------- 239 if nz eq jpk then mask = mask[*, *,niveau-1] ELSE mask = mask[*, *,nz-1] 240 if keyword_set(nan) NE 0 then begin 241 msknan = replicate(1., nx, ny) 242 notanumber = where(finite(res) EQ 0) 243 if notanumber[0] NE -1 then msknan[notanumber] = !values.f_nan 244 ENDIF ELSE msknan = 1 245 case 1 of 246 (dirx eq 1) and (diry eq 0) : begin 247 e = e1*mask 248 if keyword_set(integration) then divi=1 $ 249 else begin 250 divi = e*msknan 251 if ny EQ 1 then divi = reform(divi,nx,ny, /over) 252 divi = total(divi,1, nan = nan) 253 endelse 254 res = res*e 255 if ny EQ 1 then res = reform(res,nx,ny, /over) 256 res = total(res,1, nan = nan)/(divi > 1.) 257 if keyword_set(nan) then begin 258 testnan = finite(msknan)+(1-mask) 259 if ny EQ 1 then testnan = reform(testnan,nx,ny, /over) 260 testnan = total(testnan,1) 261 endif 262 end 263 (dirx eq 0) and (diry eq 1) : begin 264 e = e2*mask 265 if keyword_set(integration) then divi=1 $ 266 else begin 267 divi = e*msknan 268 if ny EQ 1 then divi = reform(divi,nx,ny, /over) 269 divi = total(divi,2, nan = nan) 270 endelse 271 res = res*e 272 if ny EQ 1 then res = reform(res,nx,ny, /over) 273 res = total(res,2, nan = nan)/(divi > 1.) 274 if keyword_set(nan) then begin 275 testnan = finite(msknan)+(1-mask) 276 if ny EQ 1 then testnan = reform(testnan,nx,ny, /over) 277 testnan = total(testnan,2) 278 endif 279 end 280 (dirx eq 1) and (diry eq 1) : begin 281 if keyword_set(integration) then divi=1 else divi = total(e1*e2*mask*msknan, nan = nan) 282 res = total(res*e1*e2*mask, nan = nan)/(divi > 1.) 283 if keyword_set(nan) then begin 284 testnan = finite(msknan)+(1-mask) 285 testnan = total(testnan) 286 endif 287 end 288 endcase 289 endif 275 mask = mask[*, *, 0] 276 if keyword_set(nan) NE 0 then msknan = finite(res) ELSE msknan = -1 277 case 1 of 278 (dirx eq 1) and (diry eq 0) : begin 279 e = e1*mask 280 if keyword_set(integration) then divi = 1 $ 281 else begin 282 divi = e 283 IF msknan[0] NE -1 THEN divi = temporary(divi)*msknan 284 if ny EQ 1 then divi = reform(divi, nx, ny, /over) 285 divi = total(divi, 1) 286 endelse 287 res = res*e 288 if ny EQ 1 then res = reform(res, nx, ny, /over) 289 res = total(res, 1, nan = nan)/(divi > 1.) 290 if msknan[0] NE -1 then begin 291 testnan = msknan*mask 292 if ny EQ 1 then testnan = reform(testnan, nx, ny, /over) 293 testnan = total(testnan, 1)+(total(mask, 1) EQ 0) 294 endif 295 end 296 (dirx eq 0) and (diry eq 1) : begin 297 e = e2*mask 298 if keyword_set(integration) then divi = 1 $ 299 else begin 300 divi = e 301 IF msknan[0] NE -1 THEN divi = temporary(divi)*msknan 302 if ny EQ 1 then divi = reform(divi, nx, ny, /over) 303 divi = total(divi, 2) 304 endelse 305 res = res*e 306 if ny EQ 1 then res = reform(res, nx, ny, /over) 307 res = total(res, 2, nan = nan)/(divi > 1.) 308 if msknan[0] NE -1 then begin 309 testnan = msknan*mask 310 if ny EQ 1 then testnan = reform(testnan, nx, ny, /over) 311 testnan = total(testnan, 2)+(total(mask, 2) EQ 0) 312 endif 313 end 314 (dirx eq 1) and (diry eq 1) : begin 315 if keyword_set(integration) then divi = 1 else BEGIN 316 IF msknan[0] NE -1 THEN divi = total(e1*e2*mask*msknan) $ 317 ELSE divi = total(e1*e2*mask) 318 ENDELSE 319 res = total(res*e1*e2*mask, nan = nan)/(divi > 1.) 320 if msknan[0] NE -1 then begin 321 testnan = msknan*mask 322 testnan = total(testnan)+(total(mask) EQ 0) 323 endif 324 end 325 endcase 326 endif 290 327 ;------------------------------------------------------------ 291 328 ;------------------------------------------------------------ … … 293 330 ;------------------------------------------------------------ 294 331 ;------------------------------------------------------------ 295 332 if (dim eq '3d') then begin 296 333 ;--------------------------------------------------------------- 297 334 ; III.1) verification de la coherence de la taille du tableau a … … 302 339 ; (jpi,jpj,jpk) soit celle du domaine reduit (nx,ny,ny) 303 340 ;--------------------------------------------------------------- 304 case 1 of 305 taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpk: $ 306 res=tab[premierx:dernierx, premiery:derniery, premierz:dernierz] 307 taille[1] eq jpi and taille[2] eq jpj and taille[3] eq nz: $ 308 res=tab[premierx:dernierx, premiery:derniery, *] 309 taille[1] EQ nx and taille[2] eq ny and taille[3] eq nz :res=tab 310 taille[1] EQ nx and taille[2] eq ny and taille[3] eq jpk : $ 311 res=tab[*, *, premierz:dernierz] 312 else:return, report('Probleme d''adequation entre les tailles du domaine nx*ny*nz '+strtrim(nx, 1)+'*'+strtrim(ny, 1)+'*'+strtrim(nz, 1)+' et du tableau '+strtrim(taille[1], 1)+'*'+strtrim(taille[2], 1)+'*'+strtrim(taille[3], 1)) 313 endcase 314 if keyword_set(nan) NE 0 then BEGIN 315 if nan NE 1 then BEGIN ; si nan n''est pas !values.f_nan 341 case 1 of 342 taille[1] eq jpi and taille[2] eq jpj and taille[3] eq jpk: $ 343 res = tab[firstx:lastx, firsty:lasty, firstz:lastz] 344 taille[1] eq jpi and taille[2] eq jpj and taille[3] eq nz: $ 345 res = tab[firstx:lastx, firsty:lasty, *] 346 taille[1] EQ nx and taille[2] eq ny and taille[3] eq nz :res = tab 347 taille[1] EQ nx and taille[2] eq ny and taille[3] eq jpk : $ 348 res = tab[*, *, firstz:lastz] 349 else:BEGIN 350 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat' 351 return, report('Probleme d''adequation entre les tailles du domaine nx*ny*nz '+strtrim(nx, 1)+'*'+strtrim(ny, 1)+'*'+strtrim(nz, 1)+' et du tableau '+strtrim(taille[1], 1)+'*'+strtrim(taille[2], 1)+'*'+strtrim(taille[3], 1)) 352 END 353 endcase 354 if keyword_set(nan) NE 0 then BEGIN 355 if nan NE 1 then BEGIN ; si nan n''est pas !values.f_nan 316 356 ; on le met a !values.f_nan 317 if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $ 318 ELSE notanumber = where(abs(res) GT abs(nan)/10.) 319 if notanumber[0] NE -1 then res[temporary(notanumber)] = !values.f_nan 320 ENDIF 357 if abs(nan) LT 1e6 then notanumber = where(res EQ nan) $ 358 ELSE notanumber = where(abs(res) GT abs(nan)/10.) 359 if notanumber[0] NE -1 then res[temporary(notanumber)] = !values.f_nan 321 360 ENDIF 361 ENDIF 322 362 ;--------------------------------------------------------------- 323 363 ; rq IL FAUT FAIRE ATTENTION AUX CAS OU LA DIM A MOYENNER = 1, ET … … 325 365 ; PEUVENT SEMBLER INUTILE AU DEPART 326 366 ;--------------------------------------------------------------- 327 if nx EQ 1 OR ny EQ 1 OR nz EQ 1 then BEGIN 328 res = reform(res, nx, ny, nz, /over) 329 mask = reform(mask, nx, ny, nz, /over) 330 e1 = reform(e1, nx, ny, /over) 331 e2 = reform(e2, nx, ny, /over) 332 endif 367 if nx EQ 1 OR ny EQ 1 OR nz EQ 1 then BEGIN 368 res = reform(res, nx, ny, nz, /over) 369 e1 = reform(e1, nx, ny, /over) 370 e2 = reform(e2, nx, ny, /over) 371 endif 372 if nx EQ 1 OR ny EQ 1 OR nz EQ 1 THEN $ 373 mask = reform(mask, nx, ny, nz, /over) 374 IF keyword_set(key_partialstep) THEN BEGIN 375 ; the top of the ocean floor is 376 IF vargrid EQ 'T' OR vargrid EQ 'W' THEN bottom = total(mask, 3) $ 377 ELSE bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3) 378 ; we suppress columns with only ocean or land 379 good = where(bottom NE 0 AND bottom NE nz) 380 ; the bottom of the ocean in 3D index is: 381 bottom = lindgen(nx*ny)+(temporary(bottom)-1L)*nx*ny 382 IF good[0] NE -1 THEN bottom = bottom[good] $ 383 ELSE bottom = -1 384 ENDIF ELSE bottom = -1 333 385 ;--------------------------------------------------------------- 334 386 ; III.2) differents types de moyennes 335 387 ;--------------------------------------------------------------- 336 if keyword_set(nan) NE 0 then begin 337 msknan = replicate(1., nx, ny, nz) 338 notanumber = where(finite(res) EQ 0) 339 if notanumber[0] NE -1 then msknan[notanumber] = !values.f_nan 340 ENDIF ELSE msknan = 1 341 case 1 of 342 (dirx eq 1) and (diry eq 0) and (dirz eq 0) : begin 343 e13 = e1[*]#replicate(1.,nz) 344 e13 = reform(e13,nx,ny,nz, /over) 345 if keyword_set(integration) then divi = 1 else begin 346 divi = e13*mask*msknan 347 if nz EQ 1 then divi = reform(divi,nx,ny,nz, /over) 348 divi = total(divi,1, nan = nan) 349 endelse 350 res = res*e13*mask 351 if nz EQ 1 then res = reform(res,nx,ny,nz, /over) 352 res = total(res,1, nan = nan)/(divi > 1.) 353 e13 = 1 354 if keyword_set(nan) then begin 355 testnan = finite(msknan)+(1-mask) 356 if nz EQ 1 then testnan = reform(testnan,nx,ny,nz, /over) 357 testnan = total(testnan,1) 358 endif 359 end 360 (dirx eq 0) and (diry eq 1) and (dirz eq 0) : begin 361 e23 = e2[*]#replicate(1.,nz) 362 e23 = reform(e23,nx,ny,nz, /over) 363 if keyword_set(integration) then divi =1 else begin 364 divi = e23*mask*msknan 365 if nz EQ 1 then divi = reform(divi,nx,ny,nz, /over) 366 divi = total(divi,2, nan = nan) 367 endelse 368 res = res*e23*mask 369 if nz EQ 1 then res = reform(res,nx,ny,nz, /over) 370 res = total(res,2, nan = nan)/(divi > 1.) 371 e23 = 1 372 if keyword_set(nan) then begin 373 testnan = finite(msknan)+(1-mask) 374 if nz EQ 1 then testnan = reform(testnan,nx,ny,nz, /over) 375 testnan = total(testnan,2) 376 endif 377 end 378 (dirx eq 0) and (diry eq 0) and (dirz eq 1) : begin 379 e33 = replicate(1,1.*nx*ny)#e3 380 e33 = reform(e33,nx,ny,nz, /over) 381 if keyword_set(integration) then divi =1 else begin 382 divi = e33*mask*msknan 383 if nz EQ 1 then divi = reform(divi,nx,ny,nz, /over) 384 divi = total(divi,3, nan = nan) 385 endelse 386 res = res*e33*mask 387 if nz EQ 1 then res = reform(res,nx,ny,nz, /over) 388 res = total(res,3, nan = nan)/(divi > 1.) 389 e33 = 1 390 if keyword_set(nan) then begin 391 testnan = finite(msknan)+(1-mask) 392 if nz EQ 1 then testnan = reform(testnan,nx,ny,nz, /over) 393 testnan = total(testnan,3) 394 endif 395 end 396 (dirx eq 1) and (diry eq 1) and (dirz eq 0) : begin 397 e123 = (e1*e2)[*]#replicate(1.,nz) 398 e123 = reform(e123,nx,ny,nz, /over) 399 if keyword_set(integration) then divi =1 else $ 400 divi = e123*mask*msknan 401 if nz EQ 1 then divi = reform(divi,nx,ny,nz, /over) 402 divi = total(total(divi,1, nan = nan),1, nan = nan) 403 res = res*e123*mask 404 if nz EQ 1 then res = reform(res,nx,ny,nz, /over) 405 res = total(total(res,1, nan = nan),1, nan = nan) / (divi > 1.) 406 e123 = 1 407 if keyword_set(nan) then begin 408 testnan = finite(msknan)+(1-mask) 409 if nz EQ 1 then testnan = reform(testnan,nx,ny,nz, /over) 410 testnan = total(total(testnan,1),1) 411 endif 412 end 413 (dirx eq 1) and (diry eq 0) and (dirz eq 1) : begin 414 e133 = e1[*]#e3 415 e133 = reform(e133,nx,ny,nz, /over) 416 divi = e133*mask*msknan 417 if nz EQ 1 then divi = reform(divi,nx,ny,nz, /over) 418 if keyword_set(integration) then divi =1 else $ 419 divi = total(total(divi,1, nan = nan),2, nan = nan) 420 res = res*e133*mask 421 if nz EQ 1 then res = reform(res,nx,ny,nz, /over) 422 res = total(total(res,1, nan = nan),2, nan = nan) / (divi > 1.) 423 e133 = 1 424 if keyword_set(nan) then begin 425 testnan = finite(msknan)+(1-mask) 426 if nz EQ 1 then testnan = reform(testnan,nx,ny,nz, /over) 427 testnan = total(total(testnan,1),2) 428 endif 429 end 430 (dirx eq 0) and (diry eq 1) and (dirz eq 1) : begin 431 e233 = e2[*]#e3 432 e233 = reform(e233,nx,ny,nz, /over) 433 divi = e233*mask*msknan 434 if nz EQ 1 then divi = reform(divi,nx,ny,nz, /over) 435 if keyword_set(integration) then divi =1 else $ 436 divi = total(total(divi,2, nan = nan),2, nan = nan) 437 res = res*e233*mask 438 if nz EQ 1 then res = reform(res,nx,ny,nz, /over) 439 res = total(total(res,2, nan = nan),2, nan = nan) / (divi > 1.) 440 e233 = 1 441 if keyword_set(nan) then begin 442 testnan = finite(msknan)+(1-mask) 443 if nz EQ 1 then testnan = reform(testnan,nx,ny,nz, /over) 444 testnan = total(total(testnan,2),2) 445 endif 446 end 447 (dirx eq 1) and (diry eq 1) and (dirz eq 1) : begin 448 e1233 = (e1*e2)[*]#e3 449 e1233 = reform(e1233,nx,ny,nz, /over) 450 if keyword_set(integration) then divi =1 else divi = total(e1233*mask*msknan, nan = nan) 451 res = total(res*e1233*mask, nan = nan) / (divi > 1.) 452 e1233 = 1 453 if keyword_set(nan) then begin 454 testnan = finite(msknan)+(1-mask) 455 testnan = total(testnan) 456 endif 457 end 458 endcase 459 endif 388 if keyword_set(nan) NE 0 then msknan = finite(res) ELSE msknan = -1 389 case 1 of 390 (dirx eq 1) and (diry eq 0) and (dirz eq 0) : begin 391 e13 = e1[*]#replicate(1., nz) 392 e13 = reform(e13, nx, ny, nz, /over) 393 IF keyword_set(key_partialstep) AND bottom[0] NE -1 $ 394 AND nx NE 1 THEN BEGIN 395 IF msknan[0] EQ -1 THEN BEGIN 396 msknan = replicate(1b, nx, ny, nz) 397 nan = 1 398 endif 399 msknan[bottom] = 0 400 res[bottom] = !values.f_nan 401 ENDIF 402 if keyword_set(integration) then divi = 1 else begin 403 divi = e13*mask 404 IF msknan[0] NE -1 THEN divi = temporary(divi)*msknan 405 if nz EQ 1 then divi = reform(divi, nx, ny, nz, /over) 406 divi = total(divi, 1) 407 ENDELSE 408 res = res*e13*mask 409 if nz EQ 1 then res = reform(res, nx, ny, nz, /over) 410 res = total(res, 1, nan = nan)/(divi > 1.) 411 e13 = 1 412 if msknan[0] NE -1 then begin 413 testnan = msknan*mask 414 if nz EQ 1 then testnan = reform(testnan, nx, ny, nz, /over) 415 testnan = total(testnan, 1)+(total(mask, 1) EQ 0) 416 endif 417 end 418 (dirx eq 0) and (diry eq 1) and (dirz eq 0) : begin 419 e23 = e2[*]#replicate(1., nz) 420 e23 = reform(e23, nx, ny, nz, /over) 421 IF keyword_set(key_partialstep) AND bottom[0] NE -1 $ 422 AND ny NE 1 THEN BEGIN 423 IF msknan[0] EQ -1 THEN BEGIN 424 msknan = replicate(1b, nx, ny, nz) 425 nan = 1 426 endif 427 msknan[bottom] = 0 428 res[bottom] = !values.f_nan 429 ENDIF 430 if keyword_set(integration) then divi = 1 else begin 431 divi = e23*mask 432 IF msknan[0] NE -1 THEN divi = temporary(divi)*msknan 433 if nz EQ 1 then divi = reform(divi, nx, ny, nz, /over) 434 divi = total(divi, 2) 435 ENDELSE 436 res = res*e23*mask 437 if nz EQ 1 then res = reform(res, nx, ny, nz, /over) 438 res = total(res, 2, nan = nan)/(divi > 1.) 439 e23 = 1 440 if msknan[0] NE -1 then begin 441 testnan = msknan*mask 442 if nz EQ 1 then testnan = reform(testnan, nx, ny, nz, /over) 443 testnan = total(testnan, 2)+(total(mask, 2) EQ 0) 444 endif 445 end 446 (dirx eq 0) and (diry eq 0) and (dirz eq 1) : begin 447 e33 = replicate(1, 1.*nx*ny)#e3 448 e33 = reform(e33, nx, ny, nz, /over) 449 IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN 450 IF keyword_set(wdepth) THEN $ 451 e33[bottom] = (e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $ 452 ELSE e33[bottom] = (e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)] 453 ENDIF 454 if keyword_set(integration) then divi = 1 else begin 455 divi = e33*mask 456 if msknan[0] NE -1 then divi = temporary(divi)*msknan 457 if nz EQ 1 then divi = reform(divi, nx, ny, nz, /over) 458 divi = total(divi, 3) 459 ENDELSE 460 res = res*e33*mask 461 if nz EQ 1 then res = reform(res, nx, ny, nz, /over) 462 res = total(res, 3, nan = nan)/(divi > 1.) 463 e33 = 1 464 if msknan[0] NE -1 then begin 465 testnan = msknan*mask 466 if nz EQ 1 then testnan = reform(testnan, nx, ny, nz, /over) 467 testnan = total(testnan, 3)+(total(mask, 3) EQ 0) 468 endif 469 end 470 (dirx eq 1) and (diry eq 1) and (dirz eq 0) : begin 471 e123 = (e1*e2)[*]#replicate(1., nz) 472 e123 = reform(e123, nx, ny, nz, /over) 473 IF keyword_set(key_partialstep) AND bottom[0] NE -1 $ 474 AND nx*ny NE 1 THEN BEGIN 475 IF msknan[0] EQ -1 THEN BEGIN 476 msknan = replicate(1b, nx, ny, nz) 477 nan = 1 478 endif 479 msknan[bottom] = 0 480 res[bottom] = !values.f_nan 481 ENDIF 482 if keyword_set(integration) then divi = 1 else BEGIN 483 divi = e123*mask 484 IF msknan[0] NE -1 THEN divi = temporary(divi)*msknan 485 if nz EQ 1 then divi = reform(divi, nx, ny, nz, /over) 486 divi = total(total(divi, 1), 1) 487 ENDELSE 488 res = res*e123*mask 489 if nz EQ 1 then res = reform(res, nx, ny, nz, /over) 490 res = total(total(res, 1, nan = nan), 1, nan = nan) / (divi > 1.) 491 e123 = 1 492 if msknan[0] NE -1 then begin 493 testnan = msknan*mask 494 if nz EQ 1 then testnan = reform(testnan, nx, ny, nz, /over) 495 testnan = total(total(testnan, 1), 1)+(total(total(mask, 1), 1) EQ 0) 496 endif 497 end 498 (dirx eq 1) and (diry eq 0) and (dirz eq 1) : begin 499 e133 = e1[*]#e3 500 e133 = reform(e133, nx, ny, nz, /over) 501 IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN 502 IF keyword_set(wdepth) THEN $ 503 e133[bottom] = (e1*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $ 504 ELSE e133[bottom] = (e1*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)] 505 ENDIF 506 if keyword_set(integration) then divi = 1 else BEGIN 507 divi = e133*mask 508 if msknan[0] NE -1 then divi = temporary(divi)*msknan 509 if nz EQ 1 then divi = reform(divi, nx, ny, nz, /over) 510 divi = total(total(divi, 1), 2) 511 ENDELSE 512 res = res*e133*mask 513 if nz EQ 1 then res = reform(res, nx, ny, nz, /over) 514 res = total(total(res, 1, nan = nan), 2, nan = nan) / (divi > 1.) 515 e133 = 1 516 if msknan[0] NE -1 then begin 517 testnan = msknan*mask 518 if nz EQ 1 then testnan = reform(testnan, nx, ny, nz, /over) 519 testnan = total(total(testnan, 1), 2)+(total(total(mask, 1), 2) EQ 0) 520 endif 521 end 522 (dirx eq 0) and (diry eq 1) and (dirz eq 1) : begin 523 e233 = e2[*]#e3 524 e233 = reform(e233, nx, ny, nz, /over) 525 IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN 526 IF keyword_set(wdepth) THEN $ 527 e233[bottom] = (e2*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $ 528 ELSE e233[bottom] = (e2*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)] 529 ENDIF 530 if keyword_set(integration) then divi = 1 else BEGIN 531 divi = e233*mask 532 if msknan[0] NE -1 then divi = temporary(divi)*msknan 533 if nz EQ 1 then divi = reform(divi, nx, ny, nz, /over) 534 divi = total(total(divi, 2), 2) 535 ENDELSE 536 res = res*e233*mask 537 if nz EQ 1 then res = reform(res, nx, ny, nz, /over) 538 res = total(total(res, 2, nan = nan), 2, nan = nan) / (divi > 1.) 539 e233 = 1 540 if msknan[0] NE -1 then begin 541 testnan = msknan*mask 542 if nz EQ 1 then testnan = reform(testnan, nx, ny, nz, /over) 543 testnan = total(total(testnan, 2), 2)+(total(total(mask, 2), 2) EQ 0) 544 endif 545 end 546 (dirx eq 1) and (diry eq 1) and (dirz eq 1) : begin 547 e1233 = (e1*e2)[*]#e3 548 e1233 = reform(e1233, nx, ny, nz, /over) 549 IF keyword_set(key_partialstep) AND bottom[0] NE -1 THEN BEGIN 550 IF keyword_set(wdepth) THEN $ 551 e1233[bottom] = (e1*e2*e3w_ps[firstx:lastx, firsty:lasty])[temporary(good)] $ 552 ELSE e1233[bottom] = (e1*e2*e3t_ps[firstx:lastx, firsty:lasty])[temporary(good)] 553 ENDIF 554 if keyword_set(integration) then divi = 1 else BEGIN 555 if msknan[0] NE -1 then divi = total(e1233*mask*msknan) $ 556 ELSE divi = total(e1233*mask) 557 ENDELSE 558 res = total(res*e1233*mask, nan = nan) / (divi > 1.) 559 e1233 = 1 560 if msknan[0] NE -1 then begin 561 testnan = msknan*mask 562 testnan = total(testnan)+(total(mask) EQ 0) 563 endif 564 end 565 endcase 566 endif 460 567 ;------------------------------------------------------------ 461 568 ;------------------------------------------------------------ … … 465 572 ; IV.1) on masque les terres par une valeur a 1e+20 466 573 ;------------------------------------------------------------ 467 468 469 470 471 574 valmask = 1e+20 575 terre = where(divi EQ 0) 576 IF terre[0] NE -1 THEN BEGIN 577 res(terre) = 1e+20 578 ENDIF 472 579 ;------------------------------------------------------------ 473 580 ; IV.2) on remplace, quand nan ne 1, !values.f_nan par nan 474 581 ;------------------------------------------------------------ 475 476 477 478 479 480 481 482 582 if keyword_set(nan) NE 0 then BEGIN 583 puttonan = where(testnan EQ 0) 584 if puttonan[0] NE -1 then res[puttonan] = !values.f_nan 585 if nan NE 1 then BEGIN 586 notanumber = where(finite(res) eq 0) 587 if notanumber[0] NE -1 then res[notanumber] = nan 588 ENDIF 589 ENDIF 483 590 ;------------------------------------------------------------ 484 591 ; IV.3) on se remplace ds le sous domaine qui etait definit a l''entree de 485 592 ; moyenne 486 593 ;------------------------------------------------------------ 487 if keyword_set(boite) AND NOT keyword_set(nodomdef) then domdef, oldboite,GRILLE=vargrid488 ;------------------------------------------------------------ 489 490 594 if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4moyenne.dat' 595 ;------------------------------------------------------------ 596 if keyword_set(key_performance) THEN print, 'temps moyenne', systime(1)-tempsun 597 return, res 491 598 ;------------------------------------------------------------ 492 599 ;------------------------------------------------------------ -
trunk/ToBeReviewed/CALCULS/norme.pro
r23 r25 21 21 ; KEYWORD PARAMETERS: 22 22 ; 23 ; BO ITE: boitesur laquelle moyenner (par defaut le domaine23 ; BOXZOOM: boxzoom sur laquelle moyenner (par defaut le domaine 24 24 ; selectionner par le dernier domdef effectue) 25 25 ; … … 60 60 ; pour calculer la moyenne de la norme des courants sur tout le 61 61 ; dommaine entre 0 et 50: 62 ; IDL> res=norme(un,vn,bo ite=[0,50],dir='xyz')62 ; IDL> res=norme(un,vn,boxzoom=[0,50],dir='xyz') 63 63 ; 64 64 ; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr) … … 68 68 ;------------------------------------------------------------ 69 69 ;------------------------------------------------------------ 70 FUNCTION norme,composanteu, composantev, BOITE=boite,DIREC=direc 71 @common 70 FUNCTION norme, composanteu, composantev, BOXZOOM = boxzoom, DIREC = direc, _extra = ex 71 ;--------------------------------------------------------- 72 @cm_4mesh 73 @cm_4data 74 @cm_4cal 75 IF NOT keyword_set(key_forgetold) THEN BEGIN 76 @updatenew 77 @updatekwd 78 ENDIF 79 ;--------------------------------------------------------- 72 80 tempsun = systime(1) ; pour key_performance 81 ; 82 IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $ 83 return, report(['This version of norme is based on Arakawa C-grid.' $ 84 , 'U and V grids must therefore be defined']) 85 ; 86 ;------------------------------------------------------------ 87 if keyword_set(boxzoom) then BEGIN 88 Case 1 Of 89 N_Elements(Boxzoom) Eq 1:bte = [lon1, lon2, lat1, lat2, 0., boxzoom[0]] 90 N_Elements(Boxzoom) Eq 2:bte = [lon1, lon2, lat1, lat2, boxzoom[0], boxzoom[1]] 91 N_Elements(Boxzoom) Eq 4:bte = [Boxzoom, vert1, vert2] 92 N_Elements(Boxzoom) Eq 5:bte = [Boxzoom[0:3], 0, Boxzoom[4]] 93 N_Elements(Boxzoom) Eq 6:bte = Boxzoom 94 Else: return, report('Mauvaise Definition de Boxzoom') 95 ENDCASE 96 domdef, boxzoom 97 ENDIF 73 98 ;------------------------------------------------------------ 74 99 if NOT keyword_set(direc) then direc = 0 … … 104 129 ; on trouve les points que u et v ont en communs 105 130 ;------------------------------------------------------------ 106 indicexu = (lindgen(jpi))[ premierxu:premierxu+nxu-1]107 indicexv = (lindgen(jpi))[ premierxv:premierxv+nxv-1]131 indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1] 132 indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1] 108 133 indicex = inter(indicexu, indicexv) 109 indiceyu = (lindgen(jpj))[ premieryu:premieryu+nyu-1]110 indiceyv = (lindgen(jpj))[ premieryv:premieryv+nyv-1]134 indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1] 135 indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1] 111 136 indicey = inter(indiceyu, indiceyv) 112 137 nx = n_elements(indicex) … … 124 149 indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1] 125 150 indice3d = lindgen(jpi, jpj, jpk) 126 indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt]151 indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,firstzt:lastzt] 127 152 ;------------------------------------------------------------ 128 153 ; extraction de u et v sur le domaine qui convient … … 134 159 nzt:BEGIN 135 160 if nxu NE nx then $ 136 if indicex[0] EQ premierxu then u = u[0:nx-1,*,*] ELSE u = u[1: nx, *,*]161 if indicex[0] EQ firstxu then u = u[0:nx-1,*,*] ELSE u = u[1: nx, *,*] 137 162 IF nxv NE nx THEN $ 138 if indicex[0] EQ premierxv then v = v[0:nx-1,*,*] ELSE v = v[1: nx, *,*]163 if indicex[0] EQ firstxv then v = v[0:nx-1,*,*] ELSE v = v[1: nx, *,*] 139 164 IF nyu NE ny THEN $ 140 if indicey[0] EQ premieryu then u = u[*,0:ny-1,*] ELSE u = u[*, 1: ny,*]165 if indicey[0] EQ firstyu then u = u[*,0:ny-1,*] ELSE u = u[*, 1: ny,*] 141 166 IF nyv NE ny THEN $ 142 if indicey[0] EQ premieryv then v = v[*,0:ny-1,*] ELSE v = v[*, 1: ny,*]167 if indicey[0] EQ firstyv then v = v[*,0:ny-1,*] ELSE v = v[*, 1: ny,*] 143 168 end 144 169 jpk:BEGIN 145 170 if nxu NE nx then $ 146 if indicex[0] EQ premierxu then u = u[0:nx-1, *,premierzt:dernierzt] ELSE u = u[1: nx, *,premierzt:dernierzt]171 if indicex[0] EQ firstxu then u = u[0:nx-1, *,firstzt:lastzt] ELSE u = u[1: nx, *,firstzt:lastzt] 147 172 IF nxv NE nx THEN $ 148 if indicex[0] EQ premierxv then v = v[0:nx-1, *,premierzt:dernierzt] ELSE v = v[1: nx, *,premierzt:dernierzt]173 if indicex[0] EQ firstxv then v = v[0:nx-1, *,firstzt:lastzt] ELSE v = v[1: nx, *,firstzt:lastzt] 149 174 IF nyu NE ny THEN $ 150 if indicey[0] EQ premieryu then u = u[*, 0:ny-1,premierzt:dernierzt] ELSE u = u[*, 1: ny,premierzt:dernierzt]175 if indicey[0] EQ firstyu then u = u[*, 0:ny-1,firstzt:lastzt] ELSE u = u[*, 1: ny,firstzt:lastzt] 151 176 IF nyv NE ny THEN $ 152 if indicey[0] EQ premieryv then v = v[*, 0:ny-1,premierzt:dernierzt] ELSE v = v[*, 1: ny,premierzt:dernierzt]177 if indicey[0] EQ firstyv then v = v[*, 0:ny-1,firstzt:lastzt] ELSE v = v[*, 1: ny,firstzt:lastzt] 153 178 end 154 179 ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs') … … 174 199 a=u(0,*,*) 175 200 u=(u+shift(u,1,0,0))/2. 176 if NOT keyword_set(key_periodi que) OR nx NE jpi then u(0,*,*)=a201 if NOT keyword_set(key_periodic) OR nx NE jpi then u(0,*,*)=a 177 202 a=v(*,0,*) 178 203 v=(v+shift(v,0,1,0))/2. 179 if NOT keyword_set(key_periodi que) OR nx NE jpi then v(*,0,*)=a204 if NOT keyword_set(key_periodic) OR nx NE jpi then v(*,0,*)=a 180 205 ;---------------------------------------------------------------------------- 181 206 ; attribution du mask et des tableau de longitude et latitude … … 190 215 if landv[0] NE -1 then v[landv] = 0 191 216 res=sqrt(u^2+v^2) 192 if NOT keyword_set(key_periodi que) OR nx NE jpi then res(0,*, *)=!values.f_nan217 if NOT keyword_set(key_periodic) OR nx NE jpi then res(0,*, *)=!values.f_nan 193 218 res(*,0, *)=!values.f_nan 194 219 mask = where(mask eq 0) 195 220 IF mask[0] NE -1 THEN res(mask) = valmask 196 221 ; moyennes en tous genres 197 oldboite = [lon1, lon2, lat1, lat2, prof1, prof2] 198 domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, /meme 199 if keyword_set(direc) then res = moyenne(res,direc,/nan, boite = boite) 222 domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme 223 if keyword_set(direc) then res = moyenne(res,direc,/nan, boxzoom = boxzoom, /nodomdef) 200 224 ;----------------------------------------------------------- 201 225 END … … 215 239 (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 216 240 if nxu NE nx then $ 217 if indicex[0] EQ premierxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]241 if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 218 242 IF nxv NE nx THEN $ 219 if indicex[0] EQ premierxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]243 if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 220 244 IF nyu NE ny THEN $ 221 if indicey[0] EQ premieryu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]245 if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 222 246 IF nyv NE ny THEN $ 223 if indicey[0] EQ premieryv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]247 if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 224 248 END 225 249 (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ … … 235 259 a=u(0,*,*) 236 260 u=(u+shift(u,1,0,0))/2. 237 if NOT keyword_set(key_periodi que) OR nx NE jpi then u(0,*,*)=a261 if NOT keyword_set(key_periodic) OR nx NE jpi then u(0,*,*)=a 238 262 a=v(*,0,*) 239 263 v=(v+shift(v,0,1,0))/2. 240 if NOT keyword_set(key_periodi que) OR nx NE jpi then v(*,0,*)=a264 if NOT keyword_set(key_periodic) OR nx NE jpi then v(*,0,*)=a 241 265 ;---------------------------------------------------------------------------- 242 266 ; attribution du mask et des tableau de longitude et latitude … … 245 269 ; compte (faire un petit dessin...) 246 270 ;---------------------------------------------------------------------------- 247 mask = tmask[indice2d+jpi*jpj* (niveau-1)]271 mask = tmask[indice2d+jpi*jpj*firstzt] 248 272 if ny EQ 1 then mask = reform(mask, nx, ny, /over) 249 273 ;----------------------------------------------------------- … … 256 280 if landv[0] NE -1 then v[landv] = 0 257 281 res=sqrt(u^2+v^2) 258 if NOT keyword_set(key_periodi que) OR nx NE jpi then res(0,*, *)=!values.f_nan282 if NOT keyword_set(key_periodic) OR nx NE jpi then res(0,*, *)=!values.f_nan 259 283 res(*,0, *)=!values.f_nan 260 284 mask = where(mask eq 0) … … 267 291 ENDIF 268 292 ; moyennes en tous genres 269 domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, /meme270 if keyword_set(direc) then res = grossemoyenne(res,direc,/nan, bo ite = boite)293 domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme 294 if keyword_set(direc) then res = grossemoyenne(res,direc,/nan, boxzoom = boxzoom, /nodomdef) 271 295 END 272 296 ;---------------------------------------------------------------------------- … … 279 303 indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1] 280 304 indice3d = lindgen(jpi, jpj, jpk) 281 indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt]305 indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,firstzt:lastzt] 282 306 ;------------------------------------------------------------ 283 307 ; extraction de u et v sur le domaine qui convient … … 289 313 nzt:BEGIN 290 314 if nxu NE nx then $ 291 if indicex[0] EQ premierxu then u = u[0:nx-1,*,*,*] ELSE u = u[1: nx, *,*,*]315 if indicex[0] EQ firstxu then u = u[0:nx-1,*,*,*] ELSE u = u[1: nx, *,*,*] 292 316 IF nxv NE nx THEN $ 293 if indicex[0] EQ premierxv then v = v[0:nx-1,*,*,*] ELSE v = v[1: nx, *,*,*]317 if indicex[0] EQ firstxv then v = v[0:nx-1,*,*,*] ELSE v = v[1: nx, *,*,*] 294 318 IF nyu NE ny THEN $ 295 if indicey[0] EQ premieryu then u = u[*,0:ny-1,*,*] ELSE u = u[*, 1: ny,*,*]319 if indicey[0] EQ firstyu then u = u[*,0:ny-1,*,*] ELSE u = u[*, 1: ny,*,*] 296 320 IF nyv NE ny THEN $ 297 if indicey[0] EQ premieryv then v = v[*,0:ny-1,*,*] ELSE v = v[*, 1: ny,*,*]321 if indicey[0] EQ firstyv then v = v[*,0:ny-1,*,*] ELSE v = v[*, 1: ny,*,*] 298 322 end 299 323 jpk:BEGIN 300 324 if nxu NE nx then $ 301 if indicex[0] EQ premierxu then u = u[0:nx-1, *,premierzt:dernierzt,*] ELSE u = u[1: nx, *,premierzt:dernierzt,*]325 if indicex[0] EQ firstxu then u = u[0:nx-1, *,firstzt:lastzt,*] ELSE u = u[1: nx, *,firstzt:lastzt,*] 302 326 IF nxv NE nx THEN $ 303 if indicex[0] EQ premierxv then v = v[0:nx-1, *,premierzt:dernierzt,*] ELSE v = v[1: nx, *,premierzt:dernierzt,*]327 if indicex[0] EQ firstxv then v = v[0:nx-1, *,firstzt:lastzt,*] ELSE v = v[1: nx, *,firstzt:lastzt,*] 304 328 IF nyu NE ny THEN $ 305 if indicey[0] EQ premieryu then u = u[*, 0:ny-1,premierzt:dernierzt,*] ELSE u = u[*, 1: ny,premierzt:dernierzt,*]329 if indicey[0] EQ firstyu then u = u[*, 0:ny-1,firstzt:lastzt,*] ELSE u = u[*, 1: ny,firstzt:lastzt,*] 306 330 IF nyv NE ny THEN $ 307 if indicey[0] EQ premieryv then v = v[*, 0:ny-1,premierzt:dernierzt,*] ELSE v = v[*, 1: ny,premierzt:dernierzt,*]331 if indicey[0] EQ firstyv then v = v[*, 0:ny-1,firstzt:lastzt,*] ELSE v = v[*, 1: ny,firstzt:lastzt,*] 308 332 end 309 333 ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs') … … 312 336 (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND (size(u))[3] EQ jpk AND $ 313 337 (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj AND (size(u))[3] EQ jpk :BEGIN 314 u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt, *]315 v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt, *]338 u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,firstzt:lastzt, *] 339 v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1,firstzt:lastzt, *] 316 340 END 317 341 ELSE: return, report('problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs') … … 322 346 a=u(0,*,*,*) 323 347 u=(u+shift(u,1,0,0,0))/2. 324 if NOT keyword_set(key_periodi que) OR nx NE jpi then u(0,*,*,*)=a348 if NOT keyword_set(key_periodic) OR nx NE jpi then u(0,*,*,*)=a 325 349 a=v(*,0,*,*) 326 350 v=(v+shift(v,0,1,0,0))/2. 327 if NOT keyword_set(key_periodi que) OR nx NE jpi then v(*,0,*,*)=a351 if NOT keyword_set(key_periodic) OR nx NE jpi then v(*,0,*,*)=a 328 352 ;---------------------------------------------------------------------------- 329 353 ; attribution du mask et des tableau de longitude et latitude … … 338 362 if landv[0] NE -1 then v[landv] = 0 339 363 res=sqrt(u^2+v^2) 340 if NOT keyword_set(key_periodi que) OR nx NE jpi then res(0,*, *, *)=!values.f_nan364 if NOT keyword_set(key_periodic) OR nx NE jpi then res(0,*, *, *)=!values.f_nan 341 365 res(*,0, *, *)=!values.f_nan 342 366 mask = where(mask eq 0) … … 349 373 ENDIF 350 374 ; moyennes en tous genres 351 domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, /meme352 if keyword_set(direc) then res = grossemoyenne(res,direc,/nan, bo ite = boite)375 domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme 376 if keyword_set(direc) then res = grossemoyenne(res,direc,/nan, boxzoom = boxzoom, /nodomdef) 353 377 END 354 378 ;---------------------------------------------------------------------------- … … 367 391 (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 368 392 if nxu NE nx then $ 369 if indicex[0] EQ premierxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *]393 if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *] 370 394 IF nxv NE nx THEN $ 371 if indicex[0] EQ premierxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *]395 if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *] 372 396 IF nyu NE ny THEN $ 373 if indicey[0] EQ premieryu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny]397 if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny] 374 398 IF nyv NE ny THEN $ 375 if indicey[0] EQ premieryv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny]399 if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny] 376 400 END 377 401 (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ … … 394 418 a=u(0,*) 395 419 u=(u+shift(u,1,0))/2. 396 if NOT keyword_set(key_periodi que) OR nx NE jpi then u(0,*)=a420 if NOT keyword_set(key_periodic) OR nx NE jpi then u(0,*)=a 397 421 a=v(*,0) 398 422 v=(v+shift(v,0,1))/2. 399 if NOT keyword_set(key_periodi que) OR nx NE jpi then v(*,0)=a423 if NOT keyword_set(key_periodic) OR nx NE jpi then v(*,0)=a 400 424 ;---------------------------------------------------------------------------- 401 425 ; attribution du mask et des tableau de longitude et latitude … … 404 428 ; compte (faire un petit dessin...) 405 429 ;---------------------------------------------------------------------------- 406 mask = tmask[indice2d+jpi*jpj* (niveau-1)]430 mask = tmask[indice2d+jpi*jpj*firstzt] 407 431 if nyt EQ 1 THEN mask = reform(mask, nx, ny, /over) 408 432 ;----------------------------------------------------------- … … 415 439 if landv[0] NE -1 then v[landv] = 0 416 440 res=sqrt(u^2+v^2) 417 if NOT keyword_set(key_periodi que) OR nx NE jpi then res(0,*)=!values.f_nan441 if NOT keyword_set(key_periodic) OR nx NE jpi then res(0,*)=!values.f_nan 418 442 res(*,0)=!values.f_nan 419 443 mask = where(mask eq 0) 420 444 IF mask[0] NE -1 THEN res(mask) = valmask 421 445 ; moyennes en tous genres 422 oldboite = [lon1, lon2, lat1, lat2, prof1, prof2] 423 domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, /meme 424 if keyword_set(direc) then res = moyenne(res,direc,/nan, boite = boite) 446 domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme 447 if keyword_set(direc) then res = moyenne(res,direc,/nan, boxzoom = boxzoom, /nodomdef) 425 448 END 426 449 ;---------------------------------------------------------------------------- -
trunk/ToBeReviewed/CALCULS/projectondepth.pro
r23 r25 43 43 ; ->champ nul a 1e-6 pres 44 44 ; 45 ; verifcation en projettant la temperature sur la profondeur de la 2 o45 ; verifcation en projettant la temperature sur la profondeur de la 20 46 46 ; degres par exemple... 47 47 ; … … 65 65 ; verification de la coherence entre la taille du tableau et le 66 66 ; domaine 67 grille, mask,glam,gphi,gdep,nx,ny,nz,premierx,premiery,premierz,dernierx,derniery,dernierz67 grille, mask, -1, -1, -1,nx,ny,nz,firstx,firsty,firstz,lastx,lasty,lastz 68 68 case 1 of 69 tailledepth[1] eq jpi and tailledepth[2] eq jpj:depth=depth[ premierx:dernierx, premiery:derniery]69 tailledepth[1] eq jpi and tailledepth[2] eq jpj:depth=depth[firstx:lastx, firsty:lasty] 70 70 tailledepth[1] eq nx and tailledepth[2] eq ny: 71 71 else:return, report('Probleme d''adequation entre les tailles du domaine et celle du tableau de profondeur') … … 73 73 case 1 OF 74 74 taillearray[3] NE jpk:return, report('Le tableau 3d doit avoir sa 3eme dimension egale a jpk') 75 taillearray[1] eq jpi and taillearray[2] eq jpj:array=array[ premierx:dernierx, premiery:derniery, *]75 taillearray[1] eq jpi and taillearray[2] eq jpj:array=array[firstx:lastx, firsty:lasty, *] 76 76 taillearray[1] eq nx and taillearray[2] eq ny: 77 77 else:return, report('Probleme d''adequation entre les tailles du domaine et celle du tableau de profondeur') … … 100 100 if out[0] NE -1 then res[out] = !values.f_nan 101 101 ; on masque les points terres a valmask 102 grille,mask103 102 if n_elements(valmask) EQ 0 then valmask = 1e20 104 103 terre = where((temporary(mask))[*, *, 0] EQ 0) -
trunk/ToBeReviewed/CALCULS/remplit.pro
r23 r25 1 FUNCTION remplit, zinput, NAN = nan, NITE = nite, BASIQUE = basique, mask = mask, _extra = ex1 FUNCTION remplit, zinput, NAN = nan, NITER = niter, BASIQUE = basique, mask = mask, FILLXDIR = fillxdir, FILLYDIR = fillydir, FILLVAL = fillval, _extra = ex 2 2 @common 3 tempsun = systime(1); pour key_performance3 tempsun = systime(1) ; pour key_performance 4 4 ;+ 5 6 7 5 ;; 6 ;; Extrapole zinout[jpi,jpj] sur les continents en utilisant les 4 7 ;; plus proches voisins masques oceaniquement et construit un nouveau 8 8 ; masque 9 10 ;; Reitere le processus nitefois.11 12 13 9 ;; contenant l'ancien masque oceanique PLUSles points extrapoles. 10 ;; Reitere le processus niter fois. 11 ;; C'est pas clair, essayez ! 12 ;; 13 ;; 14 14 ; 15 15 ; /Nan: to fill the point which have the value … … 18 18 ; 19 19 ; 20 21 22 23 20 ;- 24 21 ; les points non remplis sont masques a valmask 25 z = zinput 26 if n_elements(key_gridtype) EQ 0 then key_gridtype = 'c' 27 if keyword_set(basique) then begin 28 oldkey_gridtype = key_gridtype 29 key_gridtype = 'c' 30 endif 31 IF n_elements(nite) EQ 0 THEN nite = 1 32 IF nite EQ 0 THEN return, zinput 33 if keyword_set(basique) then begin 34 nx = (size(zinput))[1] 35 ny = (size(zinput))[2] 36 if NOT keyword_set(mask) then mmmask = basique ELSE mmmask = mask 37 if key_gridtype eq 'e' then begin 38 case vargrid of 39 'T':glam = glamt[premierxt:dernierxt, premieryt:dernieryt] 40 'U':glam = glamu[premierxu:dernierxu, premieryu:dernieryu] 41 endcase 42 endif 43 ENDIF ELSE grille, mmmask, glam, gphi, gdep, nx, ny, nz, _extra = ex 44 if keyword_set(mask) then mmmask = mask 45 ;--------------------------------------------------------------- 46 if (size(mmmask))[0] EQ 3 THEN BEGIN 47 if (size(mmmask))[3] EQ jpk THEN mmmask = mmmask[*, *, niveau-1] $ 48 ELSE mmmask = mmmask[*, *, nz-1] 49 ENDIF 50 ; 51 if n_elements(mmmask) EQ 1 then mmmask = replicate(1, nx, ny) 52 if keyword_set(nan) then begin 53 nanpoint = where(finite(z) EQ 0) 54 if nanpoint[0] NE -1 then begin 55 mmmask[nanpoint] = 0 56 z[nanpoint] = 0 57 endif 58 endif 22 IF n_elements(niter) EQ 0 THEN niter = 1 23 IF niter EQ 0 THEN return, zinput 24 z = zinput 25 if n_elements(key_gridtype) EQ 0 then key_gridtype = 'c' 26 if keyword_set(basique) then begin 27 oldkey_gridtype = key_gridtype 28 key_gridtype = 'c' 29 nx = (size(zinput))[1] 30 ny = (size(zinput))[2] 31 if NOT keyword_set(mask) then mmmask = basique ELSE mmmask = mask 32 if key_gridtype eq 'e' then begin 33 case vargrid of 34 'T':glam = glamt[firstxt:lastxt, firstyt:lastyt] 35 'U':glam = glamu[firstxu:lastxu, firstyu:lastyu] 36 endcase 37 endif 38 ENDIF ELSE grille, mmmask, glam, gphi, gdep, nx, ny, nz, _extra = ex 39 if keyword_set(mask) then mmmask = mask 40 ;--------------------------------------------------------------- 41 if (size(mmmask))[0] EQ 3 THEN mmmask = mmmask[*, *, 0] 42 ; 43 if n_elements(mmmask) EQ 1 then mmmask = replicate(1b, nx, ny) 44 if keyword_set(nan) then begin 45 nanpoint = where(finite(z) EQ 0) 46 if nanpoint[0] NE -1 then begin 47 mmmask[nanpoint] = 0b 48 z[nanpoint] = 0 49 endif 50 ENDIF 51 mmmask = byte(mmmask) 59 52 ;--------------------------------------------------------------- 60 53 ; on ajoute un cadre de zero a z, mask, e1, e2 … … 62 55 ; soucier des bords du domaine! 63 56 ;--------------------------------------------------------------- 64 tempdeux = systime(1); pour key_performance =265 66 67 68 zero = fltarr(nx+2, ny+2)69 zero[1:nx, 1:ny] = mmmask70 mmmask = zero71 zero= fltarr(nx+2, ny+2)72 zero[1:nx, 1:ny] = z73 if keyword_set(key_periodique) AND nx EQ jpi then begin74 zero[0, 1:ny] = z[jpi-1, *]75 zero[nx+1, 1:ny] = z[0, *]76 77 z = zero78 79 80 zero = fltarr(nx+2, ny+4)81 zero[1:nx, 2:ny+1] = mmmask82 mmmask = zero83 zero= fltarr(nx+2, ny+4)84 zero[1:nx, 2:ny+1] = z85 if keyword_set(key_periodique) AND nx EQ jpi then begin86 zero[0, 2:ny+1] = z[jpi-1, *]87 zero[nx+1, 2:ny+1] = z[0, *]88 89 z = zero90 91 92 57 tempdeux = systime(1) ; pour key_performance =2 58 nx2 = nx+2 59 case key_gridtype of 60 'c':BEGIN 61 ztmp = bytarr(nx+2, ny+2) 62 ztmp[1:nx, 1:ny] = mmmask 63 mmmask = temporary(ztmp) 64 ztmp = fltarr(nx+2, ny+2) 65 ztmp[1:nx, 1:ny] = z 66 if keyword_set(key_periodic) AND nx EQ jpi then begin 67 ztmp[0, 1:ny] = z[jpi-1, *] 68 ztmp[nx+1, 1:ny] = z[0, *] 69 endif 70 z = temporary(ztmp) 71 END 72 'e':BEGIN 73 ztmp = bytarr(nx+2, ny+4) 74 ztmp[1:nx, 2:ny+1] = mmmask 75 mmmask = temporary(ztmp) 76 ztmp = fltarr(nx+2, ny+4) 77 ztmp[1:nx, 2:ny+1] = z 78 if keyword_set(key_periodic) AND nx EQ jpi then begin 79 ztmp[0, 2:ny+1] = z[jpi-1, *] 80 ztmp[nx+1, 2:ny+1] = z[0, *] 81 endif 82 z = temporary(ztmp) 83 END 84 endcase 85 IF testvar(var = key_performance) EQ 2 THEN $ 93 86 print, 'temps remplit: on ajoute un cadre de zero ', systime(1)-tempdeux 94 87 ;--------------------------------------------------------------- … … 97 90 ;--------------------------------------------------------------- 98 91 ;--------------------------------------------------------------- 99 FOR n = 1, nite DO BEGIN 100 ; on trouve les points cote 101 tempdeux = systime(1) ; pour key_performance =2 102 ; les points du bord du cadre ne doivent pas etre dans la cote 103 case key_gridtype of 104 'c':BEGIN 105 mmmask[0, *] = 1 106 mmmask[nx+1, *] = 1 107 mmmask[*, 0] = 1 108 mmmask[*, ny+1] = 1 109 if keyword_set(key_periodique) AND nx EQ jpi then begin 110 mmmask[0,*] = mmmask[nx-2, *] 111 mmmask[nx+1,*] = mmmask[1, *] 112 endif 113 END 114 'e':BEGIN 115 mmmask[0, *] = 1 116 mmmask[nx+1, *] = 1 117 mmmask[*, 0:1] = 1 118 mmmask[*, ny+2:ny+3] = 1 119 if keyword_set(key_periodique) AND nx EQ jpi then begin 120 mmmask[0,*] = mmmask[nx-2, *] 121 mmmask[nx+1,*] = mmmask[1, *] 122 endif 123 END 124 endcase 92 FOR n = 1, niter DO BEGIN 93 ; on trouve les points coast 94 tempdeux = systime(1) ; pour key_performance =2 95 ; les points du bord du cadre ne doivent pas etre selectionnes comme 96 ; la coast 97 case key_gridtype of 98 'c':BEGIN 99 mmmask[0, *] = 1b 100 mmmask[nx+1, *] = 1b 101 mmmask[*, 0] = 1b 102 mmmask[*, ny+1] = 1b 103 END 104 'e':BEGIN 105 mmmask[0, *] = 1b 106 mmmask[nx+1, *] = 1b 107 mmmask[*, 0:1] = 1b 108 mmmask[*, ny+2:ny+3] = 1b 109 END 110 endcase 125 111 ; liste des points terre restant 126 terre = where(mmmask EQ 0) 127 if terre[0] EQ -1 then GOTO, fini 112 IF keyword_set(fillxdir) THEN BEGIN 113 ; we stop if all the lines, that contains data, have been filled 114 test = total(mmmask[1:nx, *], 1) 115 IF total((test EQ 0)+(test EQ nx)) EQ ny+2 THEN GOTO, fini 116 ENDIF 117 IF keyword_set(fillydir) THEN BEGIN 118 ; we stop if all the columns, that contains data, have been filled 119 test = total(mmmask[*, 1:ny], 2) 120 IF total((test EQ 0)+(test EQ ny)) EQ nx+2 THEN GOTO, fini 121 ENDIF 122 land = where(mmmask EQ 0) 123 if land[0] EQ -1 then GOTO, fini 128 124 ; les points du bord du cadre doivent maintenant etre dans la terre 129 case key_gridtype of 130 'c':BEGIN 131 mmmask[0, *] = 0 132 mmmask[nx+1, *] = 0 133 mmmask[*, 0] = 0 134 mmmask[*, ny+1] = 0 135 END 136 'e':BEGIN 137 mmmask[0, *] = 0 138 mmmask[nx+1, *] = 0 139 mmmask[*, 0:1] = 0 140 mmmask[*, ny+2:ny+3] = 0 141 END 142 endcase 125 case key_gridtype of 126 'c':BEGIN 127 mmmask[0, *] = 0b 128 mmmask[nx+1, *] = 0b 129 mmmask[*, 0] = 0b 130 mmmask[*, ny+1] = 0b 131 END 132 'e':BEGIN 133 mmmask[0, *] = 0b 134 mmmask[nx+1, *] = 0b 135 mmmask[*, 0:1] = 0b 136 mmmask[*, ny+2:ny+3] = 0b 137 END 138 endcase 139 if keyword_set(key_periodic) AND nx EQ jpi then begin 140 mmmask[0, *] = mmmask[nx, *] 141 mmmask[nx+1, *] = mmmask[1, *] 142 endif 143 143 ; liste des voisins mer 144 case key_gridtype of 145 'c':BEGIN 146 voisins = shift(mmmask,-1,0)+shift(mmmask,1,0)+shift(mmmask,0,-1)+shift(mmmask,0,1) $ 147 +1./sqrt(2)*(shift(mmmask,-1,-1)+shift(mmmask,1,1)+shift(mmmask,1,-1)+shift(mmmask,-1,1)) 148 cote = where(voisins[terre] GT 0) 149 ; liste des points cote 150 cote = terre[cote] 151 poids = voisins[cote] 152 END 153 'e':BEGIN 154 shifted = glam[0, 0] LT glam[0, 1] 155 oddeven = (terre/nx2+1-shifted) MOD 2 156 voisins = mmmask[terre+1]+mmmask[terre-1] $ 157 +mmmask[2*nx2+terre]+mmmask[-2*nx2+terre] $ 158 +sqrt(2)*(mmmask[terre-nx2+oddeven]+mmmask[terre+nx2+oddeven] $ 159 +mmmask[terre+(-nx2-1)+oddeven]+mmmask[terre+(nx2-1)+oddeven]) 160 cote = where(voisins GT 0) 161 poids = voisins[cote] 162 cote = terre[cote] 163 END 164 endcase 165 ; 166 IF testvar(var = key_performance) EQ 2 THEN $ 167 print, 'temps remplit: trouver la cote ', systime(1)-tempdeux 168 ;--------------------------------------------------------------- 169 ; remplissage des points cote 170 ;--------------------------------------------------------------- 171 tempdeux = systime(1) ; pour key_performance =2 144 case key_gridtype of 145 'c':BEGIN 146 CASE 1 OF 147 keyword_set(fillxdir):weight = mmmask[1+land]+mmmask[-1+land] 148 keyword_set(fillydir):weight = mmmask[nx2+land]+mmmask[-nx2+land] 149 ELSE:weight = mmmask[1+land]+mmmask[-1+land]+mmmask[nx2+land]+mmmask[-nx2+land] $ 150 +1./sqrt(2)*(mmmask[nx2+1+land]+mmmask[nx2-1+land] $ 151 +mmmask[-nx2+1+land]+mmmask[-nx2-1+land]) 152 ENDCASE 153 END 154 'e':BEGIN 155 shifted = glam[0, 0] LT glam[0, 1] 156 oddeven = (land/nx2+1-shifted) MOD 2 157 weight = mmmask[1+land]+mmmask[-1+land] $ 158 +mmmask[2*nx2+land]+mmmask[-2*nx2+land] $ 159 +sqrt(2)*(mmmask[-nx2+oddeven+land]+mmmask[-nx2-1+oddeven+land] $ 160 +mmmask[nx2+oddeven+land]+mmmask[nx2-1+oddeven+land]) 161 END 162 endcase 163 164 ok = where(weight GT 0) 165 weight = weight[ok] 166 coast = land[temporary(ok)] 167 ; 168 IF testvar(var = key_performance) EQ 2 THEN $ 169 print, 'temps remplit: trouver la coast ', systime(1)-tempdeux 170 ;--------------------------------------------------------------- 171 ; remplissage des points coast 172 ;--------------------------------------------------------------- 173 tempdeux = systime(1) ; pour key_performance =2 172 174 ; on masque z 173 z = z*mmmask 174 ; 175 case key_gridtype of 176 'c':BEGIN 177 zcote = z[1+cote]+z[-1+cote]+z[nx2+cote]+z[-nx2+cote] $ 178 +1./sqrt(2)*(z[nx2+1+cote]+z[nx2-1+cote] $ 179 +z[-nx2+1+cote]+z[-nx2-1+cote]) 180 END 181 'e':BEGIN 182 oddeven = (cote/nx2+1-shifted) MOD 2 183 zcote = z[1+cote]+z[-1+cote]+z[2*nx2+cote]+z[-2*nx2+cote] $ 184 +sqrt(2)*(z[-nx2+oddeven+cote]+z[-nx2-1+oddeven+cote] $ 185 +z[nx2+oddeven+cote]+z[nx2-1+oddeven+cote]) 186 END 187 endcase 188 189 z[cote] = zcote/poids 175 z = temporary(z)*mmmask 176 ; 177 case key_gridtype of 178 'c':BEGIN 179 CASE 1 OF 180 keyword_set(fillxdir):zcoast = z[1+coast]+z[-1+coast] 181 keyword_set(fillydir):zcoast = z[nx2+coast]+z[-nx2+coast] 182 ELSE:zcoast = z[1+coast]+z[-1+coast]+z[nx2+coast]+z[-nx2+coast] $ 183 +1./sqrt(2)*(z[nx2+1+coast]+z[nx2-1+coast] $ 184 +z[-nx2+1+coast]+z[-nx2-1+coast]) 185 ENDCASE 186 END 187 'e':BEGIN 188 oddeven = (coast/nx2+1-shifted) MOD 2 189 zcoast = z[1+coast]+z[-1+coast]+z[2*nx2+coast]+z[-2*nx2+coast] $ 190 +sqrt(2)*(z[-nx2+oddeven+coast]+z[-nx2-1+oddeven+coast] $ 191 +z[nx2+oddeven+coast]+z[nx2-1+oddeven+coast]) 192 END 193 endcase 194 ; 195 z[coast] = temporary(zcoast)/ temporary(weight) 196 ; we update the the boundary conditions of z 197 if keyword_set(key_periodic) AND nx EQ jpi then begin 198 z[0, *] = z[nx, *] 199 z[nx+1, *] = z[1, *] 200 endif 190 201 ;--------------------------------------------------------------- 191 202 ; IV) on reduit le masque 192 203 ;--------------------------------------------------------------- 193 mmmask[cote] = 1194 ; 195 196 197 204 mmmask[ temporary(coast)] = 1 205 ; 206 IF testvar(var = key_performance) EQ 2 THEN $ 207 print, 'temps remplit: une iteration ', systime(1)-tempdeux 208 ENDFOR 198 209 fini: 199 210 ; 200 ; on masque les valeurs sur les terres restantes 201 ; 202 IF n_elements(valmask) EQ 0 then valmask = 1e20 203 z = z*mmmask+valmask*(1-mmmask) 211 ; on masque les valeurs sur les lands restantes 212 ; 213 IF n_elements(valmask) EQ 0 then valmask = 1e20 214 IF n_elements(fillval) EQ 0 THEN fillval = valmask 215 z = temporary(z)*mmmask + fillval*(1b-mmmask) 204 216 ;--------------------------------------------------------------- 205 217 ; on redecoupe le tableau pour retirer le cadre! 206 218 ;--------------------------------------------------------------- 207 208 209 210 211 212 213 214 215 ; 216 217 ;--------------------------------------------------------------- 218 219 219 case key_gridtype of 220 'c':BEGIN 221 z = z[1:nx, 1:ny] 222 END 223 'e':BEGIN 224 z = z[1:nx, 2:ny+1] 225 END 226 endcase 227 ; 228 if keyword_set(basique) then key_gridtype = oldkey_gridtype 229 ;--------------------------------------------------------------- 230 if keyword_set(key_performance) THEN print, 'temps remplit', systime(1)-tempsun 231 return, z 220 232 END 221 233
Note: See TracChangeset
for help on using the changeset viewer.