Changeset 30
- Timestamp:
- 07/30/07 15:27:55 (17 years ago)
- Location:
- trunk
- Files:
-
- 1 added
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/divfred.pro
r12 r30 1 ;------------------------------------------------------------2 ;------------------------------------------------------------3 ;------------------------------------------------------------4 1 ;+ 5 ; NAME:div 6 ; 7 ; PURPOSE:calcule la divergence d'un champ 2D 8 ; 9 ; CATEGORY:calcule sur les matrices 10 ; 11 ; CALLING SEQUENCE:res=div(u,v) 12 ; 13 ; INPUTS: 14 ; u et v deux matrices representant les coordonnes d''un 15 ; champ de vecteur 16 ; 17 ; KEYWORD PARAMETERS: 18 ; 19 ; OUTPUTS:res: une matrice 2d 20 ; 21 ; COMMON BLOCKS: 22 ; common.pro 23 ; 24 ; SIDE EFFECTS: 25 ; 26 ; RESTRICTIONS: 27 ; les matrices u et v peuvent de 2 a 4 dimensions. 28 ; attention pour distinger les differents configurations de u et v 29 ; (xy, xyz, xyt, xyzt), on regarde la variable du common 2 ; @files_comments 3 ; calcule la divergence d'un champ 2D 4 ; 5 ; @categories 6 ; matrix 7 ; 8 ; @examples 9 ; IDL> res=Divfred(u,v) 10 ; 11 ; @param u {in}{required}{type=2D array} 12 ; Matrix representing the zonal coordinates (U point) 13 ; 14 ; @param v {in}{required}{type=2D array} 15 ; Matrix representing the zonal coordinates (V point) 16 ; 17 ; @returns 18 ; the divergence of the input data (with the same size) 19 ; 20 ; @uses 21 ; common.pro 22 ; 23 ; @restrictions 24 ; - Requires SAXO tools 25 ; - les matrices u et v peuvent de 2 a 4 dimensions. 26 ; attention pour distinguer les différentes configurations de u et v 27 ; (xy, xyz, xyt, xyzt), on regarde la variable du common 30 28 ; -time qui contient le calendrier en jour julien d''IDL auquel 31 ; se rapportent u et v ansi que la variable 32 ; -jpt qui est le nombre de pas de temps a considerer ds time. 33 ; les tableaux u et v sont decoupes sur le meme domaine 34 ; geographique. A cause du decalage des grilles T, U, V et F il est 35 ; possiible que ces 2 tableaux n''aient pas la meme taille et se 36 ; repportent a des indices differents. Si tel est le cas les tableaux 37 ; sont redecoupes sur les indices qu'ils ont en commun et le dommaine 38 ; est redefinit pour qu'il colle a ces indices communs. 39 ; pour eviter ces redecoupes utiliser le mot cles /memeindice ds 40 ; domdef.pro 41 ; 42 ; les points sur le bord du dessin sont mis a !values.f_nan 43 ; 44 ; EXAMPLE: 45 ; 46 ; MODIFICATION HISTORY:Guillaume Roullet (grlod\@ipsl.jussieu.fr) 29 ; se rapportent u et v ansi que la variable 30 ; -jpt qui est le nombre de pas de temps à considérer dans time. 31 ; les tableaux u et v sont découpés sur le même domaine 32 ; géographique. A cause du décalage des grilles T, U, V et F il est 33 ; possible que ces 2 tableaux n''aient pas la même taille et se 34 ; reportent à des indices différents. Si tel est le cas les tableaux 35 ; sont redécoupés sur les indices qu'ils ont en commun et le dommaine 36 ; est redéfini pour qu'il colle à ces indices communs. 37 ; pour éviter ces redécoupes utiliser le mot clé /memeindice dans 38 ; <pro>domdef</pro> 39 ; 40 ; les points sur le bord du dessin sont mis à !values.f_nan 41 ; 42 ; @todo 43 ; ++ pas fini de comprendre, tester (compatibilité saxo), adapter, commenter 44 ; nettoyer 45 ; ++ a comparer et merger avec SAXO_DIR/ToBeReviewed/CALCULS/div.pro 46 ; Written by: fplod 2006-03-27T13:45:14Z aedon.lodyc.jussieu.fr (Darwin) 47 ; 48 ; @history 49 ; fplod 2007-07-30T10:11:33Z aedon.locean-ipsl.upmc.fr (Darwin) 50 ; add header for idldoc 51 ; created from 52 ; /usr/work/sur/fvi/OPA/geomag/divfred.pro 53 ; Guillaume Roullet (grlod\@ipsl.jussieu.fr) 47 54 ; Creation : printemps 1998 48 55 ; Sebastien Masson (smasson\@lodyc.jussieu.fr) 49 56 ; adaptation pour marcher avec un domaine reduit 50 57 ; 12/1/2000 58 ; @version 59 ; $Id$ 60 ; 51 61 ;- 52 ;------------------------------------------------------------ 53 ;------------------------------------------------------------ 54 ;------------------------------------------------------------ 62 ; 55 63 FUNCTION divfred, uu, vv 56 64 tempsun = systime(1) ; pour key_performance … … 66 74 ; on trouve les points que u et v ont en communs 67 75 ;------------------------------------------------------------ 68 indicexu = (lindgen(jpi))[ premierxu:premierxu+nxu-1]69 indicexv = (lindgen(jpi))[ premierxv:premierxv+nxv-1]76 indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1] 77 indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1] 70 78 indicex = inter(indicexu, indicexv) 71 indiceyu = (lindgen(jpj))[ premieryu:premieryu+nyu-1]72 indiceyv = (lindgen(jpj))[ premieryv:premieryv+nyv-1]79 indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1] 80 indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1] 73 81 indicey = inter(indiceyu, indiceyv) 74 nx = n_elements(indicex) 82 nx = n_elements(indicex) 75 83 ny = n_elements(indicey) 76 84 indice2d = lindgen(jpi, jpj) … … 81 89 ;xyz 82 90 ;---------------------------------------------------------------------------- 83 (size(u))[0] EQ 3 AND date1 EQ date2 :BEGIN 91 (size(u))[0] EQ 3 AND date1 EQ date2 :BEGIN 84 92 ;------------------------------------------------------------ 85 93 ; extraction de u et v sur le domaine qui convient … … 88 96 (size(u))[0] NE 3 OR (size(v))[0] NE 3: return, -1 89 97 (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $ 90 (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 98 (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 91 99 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, *]100 nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 101 nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 102 nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 103 nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 96 104 ELSE : 97 105 endcase 98 106 END 99 107 (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ 100 (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 108 (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 101 109 u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *] 102 110 v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *] 103 111 END 104 ELSE:BEGIN 112 ELSE:BEGIN 105 113 zdiv = -1 106 114 GOTO, sortie 107 115 end 108 116 109 117 endcase 110 118 ;------------------------------------------------------------ … … 114 122 zu = reform(zu, nx, ny, nzt, /over) 115 123 zu = temporary(u)*temporary(zu) $ 116 *(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt]124 *(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 117 125 terreu = where(zu EQ 0) 118 126 ;;; if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan … … 121 129 zv = reform(zv, nx, ny, nzt, /over) 122 130 zv = temporary(v)*temporary(zv) $ 123 *(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt]131 *(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 124 132 terrev = where(zv EQ 0) 125 133 ;;; if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan … … 129 137 zdiv = reform(zdiv, nx, ny, nzt, /over) 130 138 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]139 *tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] 132 140 ;------------------------------------------------------------ 133 141 ; mise a !values.f_nan de la bordure 134 142 ;------------------------------------------------------------ 135 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 136 144 zdiv(0, *, *) = !values.f_nan 137 145 zdiv(nx-1, *, *) = !values.f_nan … … 143 151 ; 144 152 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)153 terre = where(tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt] EQ 0) 146 154 if terre[0] NE -1 then zdiv[temporary(terre)] = valmask 147 155 ;------------------------------------------------------------ … … 151 159 varname = 'div' 152 160 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']161 domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, grille = ['t'] 154 162 if keyword_set(direc) then zdiv = moyenne(zdiv,direc,/nan, boite = boite) 155 163 END … … 159 167 ;---------------------------------------------------------------------------- 160 168 ;---------------------------------------------------------------------------- 161 date1 NE date2 AND (size(u))[0] EQ 3 :BEGIN 169 date1 NE date2 AND (size(u))[0] EQ 3 :BEGIN 162 170 ;------------------------------------------------------------ 163 171 ; extraction de u et v sur le domaine qui convient … … 166 174 (size(u))[0] NE 3 OR (size(v))[0] NE 3: return, -1 167 175 (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $ 168 (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 176 (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 169 177 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, *]178 nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 179 nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *] 180 nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 181 nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *] 174 182 ELSE : 175 183 endcase 176 184 END 177 185 (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ 178 (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 186 (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 179 187 u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *] 180 188 v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *] … … 185 193 ; calcul de la divergence 186 194 ;------------------------------------------------------------ 187 zu = e2u[indice2d]*(umask())[indice2d+jpi*jpj* premierzt]195 zu = e2u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt] 188 196 terreu = where(zu EQ 0) 189 197 ;;; if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan … … 192 200 zu = temporary(u)*temporary(zu) 193 201 ; 194 zv = e1v[indice2d]*(vmask())[indice2d+jpi*jpj* premierzt]202 zv = e1v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt] 195 203 terrev = where(zv EQ 0) 196 204 ;;; if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan … … 199 207 zv = temporary(v)*temporary(zv) 200 208 ; 201 zdiv = 1e6*tmask[indice2d+jpi*jpj* premierzt]/(e1t[indice2d]*e2t[indice2d])209 zdiv = 1e6*tmask[indice2d+jpi*jpj*firstzt]/(e1t[indice2d]*e2t[indice2d]) 202 210 zdiv = (zdiv)[*]#replicate(1, jpt) 203 211 zdiv = reform(zdiv, nx, ny, jpt, /over) … … 207 215 ; mise a !values.f_nan de la bordure 208 216 ;------------------------------------------------------------ 209 if NOT keyword_set(key_periodi que) OR nx NE jpi then begin217 if NOT keyword_set(key_periodic) OR nx NE jpi then begin 210 218 zdiv(0, *, *) = !values.f_nan 211 219 zdiv(nx-1, *, *) = !values.f_nan … … 222 230 varname = 'div' 223 231 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']232 domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, grille = ['t'] 225 233 if keyword_set(direc) then zdiv = grossemoyenne(zdiv,direc,/nan, boite = boite) 226 234 END … … 230 238 ;---------------------------------------------------------------------------- 231 239 ;---------------------------------------------------------------------------- 232 date1 NE date2 AND (size(u))[0] EQ 4:BEGIN 240 date1 NE date2 AND (size(u))[0] EQ 4:BEGIN 233 241 return, report('non code!') 234 242 END … … 240 248 ELSE:BEGIN ;xy 241 249 indice3d = lindgen(jpi, jpj, jpk) 242 indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, niveau -1]250 indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt] 243 251 ;------------------------------------------------------------ 244 252 ; extraction de u et v sur le domaine qui convient 245 253 ;------------------------------------------------------------ 246 254 case 1 of 247 (size(u))[0] NE 2 OR (size(v))[0] NE 2: BEGIN 255 (size(u))[0] NE 2 OR (size(v))[0] NE 2: BEGIN 248 256 zdiv = -1 249 257 GOTO, sortie 250 258 end 251 259 (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $ 252 (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 260 (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 253 261 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]262 nxu NE nx:if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *] 263 nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *] 264 nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny] 265 nyv NE ny:if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny] 258 266 ELSE : 259 267 endcase 260 268 END 261 269 (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $ 262 (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 270 (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 263 271 u = u[indice2d] 264 272 v = v[indice2d] … … 280 288 ; mise a !values.f_nan de la bordure 281 289 ;------------------------------------------------------------ 282 if NOT keyword_set(key_periodi que) OR nx NE jpi then begin290 if NOT keyword_set(key_periodic) OR nx NE jpi then begin 283 291 zdiv(0, *) = !values.f_nan 284 292 zdiv(nx-1, *) = !values.f_nan … … 298 306 varname = 'div' 299 307 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']308 domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, grille = ['t'] 301 309 if keyword_set(direc) then zdiv = moyenne(zdiv,direc,/nan, boite = boite) 302 310 303 311 END 304 312 endcase 305 313 306 314 ;------------------------------------------------------------ 307 315 sortie: 308 if keyword_set(key_performance) THEN print, 'temps div', systime(1)-tempsun 309 316 if keyword_set(key_performance) THEN print, 'temps div', systime(1)-tempsun 317 310 318 return, zdiv 311 319 end -
trunk/forcagequimarche.pro
r29 r30 1 ; 1 ;+ 2 2 ; 3 3 ; @file_comments … … 31 31 ; provide tools to plot output files 32 32 ; produce a NetCDF GDT or CF compliant 33 ; add divchoice parameter to choose either div or divfred 33 34 ; 34 35 ; @pre … … 46 47 ; 47 48 ; @history 49 ; fplod 2007-07-30T10:12:38Z aedon.locean-ipsl.upmc.fr (Darwin) 50 ; restore divfred call instead of div ("official" saxo divergence module) 48 51 ; fplod 2007-07-27T10:29:18Z cerbere.locean-ipsl.upmc.fr (Linux) 49 52 ; add orcares in the name of output files … … 232 235 ; 233 236 ; Divergence du champ 234 Diver=div (Bu_star,Bv_star)237 Diver=divfred(Bu_star,Bv_star) 235 238 236 239 Diver=Diver*1e-6
Note: See TracChangeset
for help on using the changeset viewer.