[2] | 1 | ;------------------------------------------------------------ |
---|
| 2 | ;------------------------------------------------------------ |
---|
| 3 | ;------------------------------------------------------------ |
---|
| 4 | ;+ |
---|
| 5 | ; NAME:ncdf_meshlec |
---|
| 6 | ; |
---|
| 7 | ; PURPOSE:lit le fichier meshmask au format netcdf cree par OPA ds sa |
---|
| 8 | ; version tout NETCDF! |
---|
| 9 | ; |
---|
| 10 | ; CATEGORY:lecture de grille |
---|
| 11 | ; |
---|
| 12 | ; CALLING SEQUENCE:ncdf_meshlec [,' nomfich'] |
---|
| 13 | ; |
---|
| 14 | ; INPUTS: |
---|
| 15 | ; |
---|
| 16 | ; nomfich: un string definissant le nom du fichier a lire. Si |
---|
| 17 | ; ce string ne contient pas l''adresse entiere du fichier nomfich est |
---|
| 18 | ; complete par iodir. Par defaut nomfich=iodir+'meshmask.nc' |
---|
| 19 | ; |
---|
| 20 | ; KEYWORD PARAMETERS: |
---|
| 21 | ; |
---|
| 22 | ; GLAMBOUNDARY: un vecteur de 2 elements specifaint le min et le |
---|
| 23 | ; max qui doivent etre imposes en longitude (obligatoire si le |
---|
| 24 | ; tableau depasse 360 degres) |
---|
| 25 | ; |
---|
| 26 | ; /CHECKDAT: cherche si il existe un ficher du meme nom que le |
---|
| 27 | ; fichier meshmask mais terminant par .dat au lieu de .dat. |
---|
| 28 | ; Si ce fichier n''existe pas, ncdf_meshlec lit le fichier |
---|
| 29 | ; NetCdf et sauve tous les tableaux lus ds le fichier |
---|
| 30 | ; ...dat. Rq: Ce fichier est bcp plus petit que le fichier |
---|
| 31 | ; NetCdf d''origine car par ex on ne sauve pas les tableaux |
---|
| 32 | ; umask, vmask, fmask mais simplement certains de leurs bords |
---|
| 33 | ; (cf. umask.pro vmask.pro et fmask.pro) et ces masks sont |
---|
| 34 | ; sauves en bytes et toues les autres tableaux en |
---|
| 35 | ; float. Cependant ce fichier .dat est cree en fonction des |
---|
| 36 | ; parametres rentres ds le init... (ixminmesh.....) et Je ne |
---|
| 37 | ; sais rien sur la portabilite du fichier .dat. |
---|
| 38 | ; Si ce fichier .dat existe on le lit avec restore. |
---|
| 39 | ; |
---|
| 40 | ; OUTPUTS:mise a jours des variables de common glam, common |
---|
| 41 | ; gphi,common e12,common mask,common tab (cf. common.pro) |
---|
| 42 | ; |
---|
| 43 | ; COMMON BLOCKS:common.pro |
---|
| 44 | ; |
---|
| 45 | ; SIDE EFFECTS:prend en compte (si il sont definits) les |
---|
| 46 | ; ixminmesh,...et definit jpiglo,jpi,.... |
---|
| 47 | ; |
---|
| 48 | ; RESTRICTIONS: |
---|
| 49 | ; |
---|
| 50 | ; La definition de ixminmesh,ixmaxmesh,iyminmesh,iymaxmesh |
---|
| 51 | ; ,izminmesh,izmaxmesh doit etre faite avant l''entree dans cette |
---|
| 52 | ; routine. pour attribuer automatiquement ces valeurs au maximum |
---|
| 53 | ; possible les mettre toutes a -1 et ncdf_meshlec les calculera. |
---|
| 54 | ; |
---|
| 55 | ; EXAMPLE: |
---|
| 56 | ; |
---|
| 57 | ; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr) |
---|
| 58 | ; 12/1999 |
---|
| 59 | ;- |
---|
| 60 | ;------------------------------------------------------------ |
---|
| 61 | ;------------------------------------------------------------ |
---|
| 62 | ;------------------------------------------------------------ |
---|
| 63 | PRO ncdf_meshlec, nomfich, GLAMBOUNDARY = glamboundary, GETDIMENSIONS = getdimensions, CHECKDAT = checkdat |
---|
| 64 | @common |
---|
| 65 | tempsun = systime(1) ; pour key_performance |
---|
| 66 | ;---------------------------------------------------------- |
---|
| 67 | ; constitution de l''adresse s_fichier et |
---|
| 68 | ; ouverture du fichier a l''adresse s_fichier |
---|
| 69 | ;------------------------------------------------------- |
---|
| 70 | ; def de nomfich par defaut |
---|
| 71 | IF n_params() EQ 0 then nomfich = 'meshmask.nc' |
---|
| 72 | ;------------------ |
---|
| 73 | if keyword_set(CHECKDAT) then begin |
---|
| 74 | nomfichdat = strmid(nomfich, 0, strlen(nomfich)-2)+'dat' |
---|
| 75 | thisOS = strupcase(strmid(!version.os_family, 0, 3)) |
---|
| 76 | CASE thisOS OF |
---|
| 77 | 'MAC':sep = ':' |
---|
| 78 | 'WIN':sep = '\' |
---|
| 79 | ELSE:sep = '/' |
---|
| 80 | ENDCASE |
---|
| 81 | if strpos(nomfichdat, sep) EQ -1 then begin |
---|
| 82 | ; si iodir ne finit pas par sep on le rajoute |
---|
| 83 | if rstrpos(iodir, sep) NE strlen(iodir)-1 then iodir = iodir+sep |
---|
| 84 | nomfichdat = iodir+nomfichdat |
---|
| 85 | endif |
---|
| 86 | test = findfile(nomfichdat) ; le nom cherche correspond bien a un fichier? |
---|
| 87 | if test[0] NE '' then begin |
---|
| 88 | noticebase=xnotice('Lecture du fichier !C '+nomfichdat+'!C ...') |
---|
| 89 | restore, nomfichdat |
---|
| 90 | widget_control, noticebase, bad_id = nothing, /destroy |
---|
| 91 | if keyword_set(key_performance) THEN print, 'temps ncdf_meshlec', systime(1)-tempsun |
---|
| 92 | return |
---|
| 93 | ENDIF |
---|
| 94 | endif |
---|
| 95 | ;------------------ |
---|
| 96 | s_fichier = isafile(file = nomfich, iodir = iodir) |
---|
| 97 | ; |
---|
| 98 | noticebase=xnotice('Lecture du fichier !C '+s_fichier+'!C ...') |
---|
| 99 | ; |
---|
| 100 | cdfid=ncdf_open(s_fichier) |
---|
| 101 | contient=ncdf_inquire(cdfid) |
---|
| 102 | |
---|
| 103 | ;------------------------------------------------------------ |
---|
| 104 | ; differentes dimensions |
---|
| 105 | ;------------------------------------------------------------ |
---|
| 106 | ncdf_diminq,cdfid,'x',name,jpiglo |
---|
| 107 | ncdf_diminq,cdfid,'y',name,jpjglo |
---|
| 108 | ncdf_diminq,cdfid,'z',name,jpkglo |
---|
| 109 | ; |
---|
| 110 | if keyword_set(getdimensions) then begin |
---|
| 111 | ncdf_close, cdfid |
---|
| 112 | return |
---|
| 113 | endif |
---|
| 114 | ;------------------------------------------------------- |
---|
| 115 | ;------------------------------------------------------- |
---|
| 116 | if n_elements(ixminmesh) EQ 0 THEN ixminmesh = 0 |
---|
| 117 | if n_elements(ixmaxmesh) EQ 0 then ixmaxmesh = jpiglo-1 |
---|
| 118 | if ixminmesh EQ -1 THEN ixminmesh = 0 |
---|
| 119 | IF ixmaxmesh EQ -1 then ixmaxmesh = jpiglo-1 |
---|
| 120 | if n_elements(iyminmesh) EQ 0 THEN iyminmesh = 0 |
---|
| 121 | IF n_elements(iymaxmesh) EQ 0 then iymaxmesh = jpjglo-1 |
---|
| 122 | if iyminmesh EQ -1 THEN iyminmesh = 0 |
---|
| 123 | IF iymaxmesh EQ -1 then iymaxmesh = jpjglo-1 |
---|
| 124 | if n_elements(izminmesh) EQ 0 THEN izminmesh = 0 |
---|
| 125 | IF n_elements(izmaxmesh) EQ 0 then izmaxmesh = jpkglo-1 |
---|
| 126 | if izminmesh EQ -1 THEN izminmesh = 0 |
---|
| 127 | IF izmaxmesh EQ -1 then izmaxmesh = jpkglo-1 |
---|
| 128 | ; |
---|
| 129 | ixmindtasauve = testvar(var = ixmindta) |
---|
| 130 | iymindtasauve = testvar(var = iymindta) |
---|
| 131 | izmindtasauve = testvar(var = izmindta) |
---|
| 132 | ; |
---|
| 133 | ixmindta = 0l |
---|
| 134 | iymindta = 0l |
---|
| 135 | izmindta = 0l |
---|
| 136 | ; |
---|
| 137 | jpi = long(ixmaxmesh-ixminmesh+1) |
---|
| 138 | jpj = long(iymaxmesh-iyminmesh+1) |
---|
| 139 | jpk = long(izmaxmesh-izminmesh+1) |
---|
| 140 | ; |
---|
| 141 | if n_elements(key_stride) LE 2 then key_stride = [1, 1, 1] |
---|
| 142 | key_stride = 1l > long(key_stride) |
---|
| 143 | ; |
---|
| 144 | jpitotal = jpi |
---|
| 145 | key_shift = long(testvar(var = key_shift)) |
---|
| 146 | key = long(key_shift MOD jpitotal) |
---|
| 147 | if key LT 0 then key = key+jpitotal |
---|
| 148 | ixmin = ixminmesh-ixmindta |
---|
| 149 | ; |
---|
| 150 | jpi = jpi/key_stride[0] |
---|
| 151 | jpj = jpj/key_stride[1] |
---|
| 152 | jpk = jpk/key_stride[2] |
---|
| 153 | ; |
---|
| 154 | jpt = 1 |
---|
| 155 | premiertps = 0 |
---|
| 156 | ; |
---|
| 157 | premierx = 0 |
---|
| 158 | dernierx = jpi-1 |
---|
| 159 | premiery = 0 |
---|
| 160 | derniery = jpj-1 |
---|
| 161 | premierz = 0 |
---|
| 162 | dernierz = jpk-1 |
---|
| 163 | nx = jpi |
---|
| 164 | ny = jpj |
---|
| 165 | nz = 1 |
---|
| 166 | izminmeshsauve = izminmesh |
---|
| 167 | izminmesh = 0 |
---|
| 168 | ; |
---|
| 169 | ;------------------------------------------------------- |
---|
| 170 | ; tableaux 2d: |
---|
| 171 | ;------------------------------------------------------- |
---|
| 172 | if total(key_stride) EQ 3 then $ |
---|
| 173 | namevar = ['glamt','glamu','glamv','glamf','gphit','gphiu','gphiv','gphif','e1t','e1u','e1v','e1f','e2t','e2u','e2v','e2f'] $ |
---|
| 174 | ELSE namevar = ['glamt','glamu','glamv','gphit','gphiu','gphiv'] |
---|
| 175 | ; |
---|
| 176 | for i = 0, n_elements(namevar)-1 do begin |
---|
| 177 | varcontient=ncdf_varinq(cdfid, namevar[i]) |
---|
| 178 | name = varcontient.name |
---|
| 179 | @read_ncdf_varget |
---|
| 180 | commande = namevar[i]+'=float(res)' |
---|
| 181 | rien = execute(commande) |
---|
| 182 | ENDFOR |
---|
| 183 | ; |
---|
| 184 | if total(key_stride) NE 3 then BEGIN |
---|
| 185 | glamf = glamt |
---|
| 186 | glamf = [glamf, reform(2.*reform(glamt[jpi-1, *])-reform(glamt[jpi-2, *]), 1, jpj)] |
---|
| 187 | glamf = [ [glamf], [reform(2*glamf[*, jpj-1]-glamf[*, jpj-2], jpi+1, 1)] ] |
---|
| 188 | glamf = glamf+shift(glamf, -1, -1) |
---|
| 189 | glamf = .5*glamf[0:jpi-1, 0:jpj-1] |
---|
| 190 | ; |
---|
| 191 | gphif = gphit |
---|
| 192 | gphif = [gphif, reform(2.*reform(gphit[jpi-1, *])-reform(gphit[jpi-2, *]), 1, jpj)] |
---|
| 193 | gphif = [ [gphif], [reform(2*gphif[*, jpj-1]-gphif[*, jpj-2], jpi+1, 1)] ] |
---|
| 194 | gphif = gphif+shift(gphif, -1, -1) |
---|
| 195 | gphif = .5*gphif[0:jpi-1, 0:jpj-1] |
---|
| 196 | ; |
---|
| 197 | e1t = replicate(1, jpi, jpj) |
---|
| 198 | e2t = e1t |
---|
| 199 | e1u = e1t |
---|
| 200 | e2u = e1t |
---|
| 201 | e1v = e1t |
---|
| 202 | e2v = e1t |
---|
| 203 | e1f = e1t |
---|
| 204 | e2f = e1t |
---|
| 205 | ENDIF |
---|
| 206 | ;------------------------------------------------------- |
---|
| 207 | ; tableaux 3d: |
---|
| 208 | ;------------------------------------------------------- |
---|
| 209 | nz = jpk |
---|
| 210 | izminmesh = izminmeshsauve |
---|
| 211 | ; |
---|
| 212 | varcontient=ncdf_varinq(cdfid,'tmask') |
---|
| 213 | name = varcontient.name |
---|
| 214 | @read_ncdf_varget |
---|
| 215 | tmask = byte(res) |
---|
| 216 | ; |
---|
| 217 | varcontient=ncdf_varinq(cdfid,'umask') |
---|
| 218 | name = varcontient.name |
---|
| 219 | nx = 1 |
---|
| 220 | premierx = jpi-1 |
---|
| 221 | @read_ncdf_varget |
---|
| 222 | umaskred = byte(res) |
---|
| 223 | umaskred = reform(umaskred, /over) |
---|
| 224 | ; |
---|
| 225 | varcontient=ncdf_varinq(cdfid,'vmask') |
---|
| 226 | name = varcontient.name |
---|
| 227 | nx = jpi |
---|
| 228 | premierx = 0 |
---|
| 229 | ny = 1 |
---|
| 230 | premiery = jpj-1 |
---|
| 231 | @read_ncdf_varget |
---|
| 232 | vmaskred = byte(res) |
---|
| 233 | vmaskred = reform(vmaskred, /over) |
---|
| 234 | ; |
---|
| 235 | varcontient=ncdf_varinq(cdfid,'fmask') |
---|
| 236 | name = varcontient.name |
---|
| 237 | nx = 1 |
---|
| 238 | premierx = jpi-1 |
---|
| 239 | ny = jpj |
---|
| 240 | premiery = 0 |
---|
| 241 | @read_ncdf_varget |
---|
| 242 | fmaskredy = byte(res) |
---|
| 243 | coast = where(fmaskredy NE 0 and fmaskredy NE 1) |
---|
| 244 | IF coast[0] NE -1 THEN fmaskredy[coast] = 0b |
---|
| 245 | fmaskredy= reform(fmaskredy, /over) |
---|
| 246 | nx = jpi |
---|
| 247 | premierx = 0 |
---|
| 248 | ny = 1 |
---|
| 249 | premiery = jpj-1 |
---|
| 250 | @read_ncdf_varget |
---|
| 251 | fmaskredx = byte(res) |
---|
| 252 | coast = where(fmaskredx NE 0 and fmaskredx NE 1) |
---|
| 253 | IF coast[0] NE -1 THEN fmaskredx[coast] = 0b |
---|
| 254 | fmaskredx = reform(fmaskredx, /over) |
---|
| 255 | ;------------------------------------------------------- |
---|
| 256 | ; tableaux 1d: |
---|
| 257 | ;------------------------------------------------------- |
---|
| 258 | namevar = ['e3t','e3w','gdept','gdepw'] |
---|
| 259 | for i = 0, n_elements(namevar)-1 do begin |
---|
| 260 | varcontient=ncdf_varinq(cdfid, namevar[i]) |
---|
| 261 | if n_elements(varcontient.dim) EQ 4 THEN BEGIN |
---|
| 262 | commande = 'ncdf_varget,cdfid,namevar[i],'+namevar[i] $ |
---|
| 263 | +',offset = [0,0,izminmesh,0], count = [1,1,jpk,1]' |
---|
| 264 | if key_stride[2] NE 1 then commande = commande+', stride=[1,1,key_stride[2],1]' |
---|
| 265 | ENDIF ELSE BEGIN |
---|
| 266 | commande = 'ncdf_varget,cdfid,namevar[i],'+namevar[i] $ |
---|
| 267 | +',offset = [izminmesh], count = [jpk]' |
---|
| 268 | if key_stride[2] NE 1 then commande = commande+', stride=key_stride[2]' |
---|
| 269 | ENDELSE |
---|
| 270 | rien = execute(commande) |
---|
| 271 | commande = namevar[i]+'=float('+namevar[i]+')' |
---|
| 272 | rien = execute(commande) |
---|
| 273 | commande = namevar[i]+' = reform('+namevar[i]+', /over)' |
---|
| 274 | rien = execute(commande) |
---|
| 275 | ENDFOR |
---|
| 276 | ;------------------------------------------------------- |
---|
| 277 | ncdf_close, cdfid |
---|
| 278 | ;------------------------------------------------------- |
---|
| 279 | ; bornes de glam qui ne doivent pas depasser 360 degres... |
---|
| 280 | ;------------------------------------------------------- |
---|
| 281 | ; minglam = min(glamt, max = maxglam) |
---|
| 282 | ; if maxglam-minglam GE 360 AND NOT keyword_set(glamboundary) then $ |
---|
| 283 | ; nothing = execute('glamboundary = '+xquestion('What are the longitudes boundary?', '[-180,180]',/chkwidget)) |
---|
| 284 | ;------------------------------------------------------- |
---|
| 285 | if keyword_set(glamboundary) then BEGIN |
---|
| 286 | if glamboundary[0] NE glamboundary[1] then BEGIN |
---|
| 287 | glamt = glamt MOD 360 |
---|
| 288 | smaller = where(glamt LT glamboundary[0]) |
---|
| 289 | if smaller[0] NE -1 then glamt[smaller] = glamt[smaller]+360 |
---|
| 290 | bigger = where(glamt GE glamboundary[1]) |
---|
| 291 | if bigger[0] NE -1 then glamt[bigger] = glamt[bigger]-360 |
---|
| 292 | glamu = glamu MOD 360 |
---|
| 293 | smaller = where(glamu LT glamboundary[0]) |
---|
| 294 | if smaller[0] NE -1 then glamu[smaller] = glamu[smaller]+360 |
---|
| 295 | bigger = where(glamu GE glamboundary[1]) |
---|
| 296 | if bigger[0] NE -1 then glamu[bigger] = glamu[bigger]-360 |
---|
| 297 | glamv = glamv MOD 360 |
---|
| 298 | smaller = where(glamv LT glamboundary[0]) |
---|
| 299 | if smaller[0] NE -1 then glamv[smaller] = glamv[smaller]+360 |
---|
| 300 | bigger = where(glamv GE glamboundary[1]) |
---|
| 301 | if bigger[0] NE -1 then glamv[bigger] = glamv[bigger]-360 |
---|
| 302 | glamf = glamf MOD 360 |
---|
| 303 | smaller = where(glamf LT glamboundary[0]) |
---|
| 304 | if smaller[0] NE -1 then glamf[smaller] = glamf[smaller]+360 |
---|
| 305 | bigger = where(glamf GE glamboundary[1]) |
---|
| 306 | if bigger[0] NE -1 then glamf[bigger] = glamf[bigger]-360 |
---|
| 307 | endif |
---|
| 308 | endif |
---|
| 309 | ;------------------------------------------------------- |
---|
| 310 | ixmindta = ixmindtasauve |
---|
| 311 | iymindta = iymindtasauve |
---|
| 312 | izmindta = izmindtasauve |
---|
| 313 | ;------------------------------------------------------- |
---|
| 314 | widget_control, noticebase, bad_id = nothing, /destroy |
---|
| 315 | ; |
---|
| 316 | if keyword_set(checkdat) then BEGIN |
---|
| 317 | noticebase=xnotice('Ecriture du fichier !C '+nomfichdat+'!C ...') |
---|
| 318 | ; |
---|
| 319 | save, jpi, jpj, jpk, ixminmesh, ixmaxmesh, iyminmesh, iymaxmesh, izminmesh, izmaxmesh, glamt,glamu,glamv,glamf,gphit,gphiu,gphiv,gphif,e1t,e1u,e1v,e1f,e2t,e2u,e2v,e2f, tmask, umaskred, vmaskred, fmaskredx, fmaskredy, gdept, gdepw, e3t, e3w, filename = nomfichdat |
---|
| 320 | ; |
---|
| 321 | widget_control, noticebase, bad_id = nothing, /destroy |
---|
| 322 | endif |
---|
| 323 | ; |
---|
| 324 | if keyword_set(key_performance) THEN print, 'temps ncdf_meshlec', systime(1)-tempsun |
---|
| 325 | |
---|
| 326 | ;------------------------------------------------------- |
---|
| 327 | return |
---|
| 328 | end |
---|