Changeset 169
- Timestamp:
- 11/16/09 15:39:10 (15 years ago)
- Location:
- trunk
- Files:
-
- 151 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/procs/ajoutvectz.pro
r165 r169 49 49 , TYPE = type $ 50 50 , BOITE = boite 51 ; 52 compile_opt idl2, strictarrsubs 51 53 ; 52 54 @common … … 129 131 130 132 terre = where(u GE valmask/10.) 131 IF terre[0] NE -1 THEN u (terre)= !VALUES.F_NAN133 IF terre[0] NE -1 THEN u[terre] = !VALUES.F_NAN 132 134 terre = where(w GE valmask/10.) 133 IF terre[0] NE -1 THEN w (terre)= !VALUES.F_NAN134 135 a=u (0,*)135 IF terre[0] NE -1 THEN w[terre] = !VALUES.F_NAN 136 137 a=u[0,*] 136 138 u=(u+shift(u,1,0))/2. 137 u (0,*)=a138 139 a=w (*,0)139 u[0,*]=a 140 141 a=w[*,0] 140 142 w=(w+shift(w,0,1))/2. 141 w (*,0)= a143 w[*,0] = a 142 144 143 145 vargrid='T' … … 171 173 172 174 CASE type OF 173 'yz': xgrid = gphit ((firstxt+lastxt)/2, indicex[0]:indicex[nx-1])174 'xz': xgrid = glamt (indicex[0]:indicex[nx-1], (firstyt+lastyt)/2)175 'yz': xgrid = gphit[(firstxt+lastxt)/2, indicex[0]:indicex[nx-1]] 176 'xz': xgrid = glamt[indicex[0]:indicex[nx-1], (firstyt+lastyt)/2] 175 177 ENDCASE 176 178 177 179 x0 = xgrid[*]#replicate(1, nz) 178 y0 = replicate(1, nx)#gdept (indicez[0]:indicez[nz-1])180 y0 = replicate(1, nx)#gdept[indicez[0]:indicez[nz-1]] 179 181 180 182 -
trunk/procs/alpha.pro
r165 r169 9 9 FUNCTION alpha, t, s, p 10 10 ; 11 compile_opt idl2, strictarrsubs 12 ; 11 13 dt = 0.05 12 14 siga = eos(t, s)-1000. -
trunk/procs/betar.pro
r165 r169 9 9 FUNCTION betar, t, s, p 10 10 ; 11 compile_opt idl2, strictarrsubs 12 ; 11 13 ds = 0.01 12 14 siga = eos(t, s)-1000. -
trunk/procs/bin_sigma.pro
r168 r169 24 24 PRO bin_sigma, cmd, pfildi 25 25 ; 26 26 compile_opt idl2, strictarrsubs 27 ; 27 28 @common 28 29 @com_eg -
trunk/procs/bining2.pro
r168 r169 88 88 , E3W=e3w $ 89 89 , TMASK=tmask 90 90 ; 91 compile_opt idl2, strictarrsubs 92 ; 91 93 size3d = size(density) 92 jpi = size3d (1)93 jpj = size3d (2)94 jpk = size3d (3)94 jpi = size3d[1] 95 jpj = size3d[2] 96 jpk = size3d[3] 95 97 ;print, 'size3d=', size3d 96 98 ;print, 'size(tmask)=', size(tmask) … … 194 196 IF i_ocean[0] NE -1 THEN BEGIN ; on n entre que si il y a des points ocean 195 197 ; print, 'ocean point' 196 i_bottom = i_ocean (n_elements(i_ocean)-1)198 i_bottom = i_ocean[n_elements[i_ocean]-1] 197 199 ;print, 'i_bottom=', i_bottom 198 z_s(N_s) = z_zw (i_bottom)199 c1_s(N_s) = x1_content(i, j, jpk-1 )200 z_s(N_s) = z_zw[i_bottom] 201 c1_s(N_s) = x1_content(i, j, jpk-1] 200 202 201 s_z (*) = density(i, j, *)202 c1_z (*) = x1_content(i, j, *)203 s_z[*] = density[i, j, *] 204 c1_z[*] = x1_content[i, j, *] 203 205 ; print, 'density profile s_z', s_z 204 206 ; print, 'field profile c1_z', c1_z … … 207 209 208 210 ; extraction d'un sous profil strictement croissant 209 ; print, 's_z (i_ocean)=', s_z(i_ocean)210 mini = min(s_z (i_ocean))211 maxi = max(s_z (i_ocean))212 i_min = where(s_z (i_ocean)EQ mini)213 i_max = where(s_z (i_ocean)EQ maxi)211 ; print, 's_z[i_ocean]=', s_z[i_ocean] 212 mini = min(s_z[i_ocean]) 213 maxi = max(s_z[i_ocean]) 214 i_min = where(s_z[i_ocean] EQ mini) 215 i_max = where(s_z[i_ocean] EQ maxi) 214 216 ; print, 'mini, maxi', mini, maxi 215 217 ; print, 'i_min, i_max', i_min, i_max … … 228 230 ; l isopycne est mise en surface (z_s=0) 229 231 ; 230 ; print, 's_z (i_min)', s_z(i_min)232 ; print, 's_z[i_min]', s_z[i_min] 231 233 ;print, 's_s=', s_s 232 ind = where(s_s LT s_z (i_min))234 ind = where(s_s LT s_z[i_min]) 233 235 IF ind[0] NE -1 THEN BEGIN 234 236 ; IF i_min GT i_max THEN print, 'min reached at sigma indices', ind 235 z_s (ind)= 0236 c1_s (ind)= !VALUES.F_NAN237 z_s[ind] = 0 238 c1_s[ind] = !VALUES.F_NAN 237 239 ENDIF 238 240 … … 240 242 ; l isopycne est mise au fond (z_s=z_zw(i_bottom)) 241 243 ; 242 ; print, 's_z (i_max)', s_z(i_max)243 ind = where(s_s GT s_z (i_max))244 ; print, 's_z[i_max]', s_z[i_max] 245 ind = where(s_s GT s_z[i_max]) 244 246 245 247 IF ind[0] NE -1 THEN BEGIN 246 248 ; IF i_min GT i_max THEN print, 'max reached at sigma indices', ind 247 z_s (ind)= z_s(N_s)248 c1_s (ind)= c1_s(N_s)249 c1_s (ind)= !VALUES.F_NAN249 z_s[ind] = z_s(N_s) 250 c1_s[ind] = c1_s(N_s) 251 c1_s[ind] = !VALUES.F_NAN 250 252 ENDIF 251 253 ; cas general 252 ind = where( (s_s GE s_z (i_min)) AND (s_s LE s_z(i_max)) )254 ind = where( (s_s GE s_z[i_min]) AND (s_s LE s_z[i_max]) ) 253 255 ; IF i_min GT i_max THEN print, 's_s indices inside min/max', ind 254 256 IF ind[0] NE -1 THEN BEGIN 255 i_profil = i_ocean (i_min:i_max); problem line256 257 z_s (ind) = interpol(z_zt(i_profil), s_z(i_profil), s_s(ind)) ; original258 ; z_s (ind) = interpol(z_zw(i_profil), s_z(i_profil), s_s(ind)) ; changed to z_zw 1/7/04257 i_profil = i_ocean[i_min:i_max] ; problem line 258 259 z_s[ind] = interpol(z_zt[i_profil], s_z[i_profil], s_s[ind]) ; original 260 ; z_s[ind] = interpol(z_zw[i_profil], s_z[i_profil], s_s[ind]) ; changed to z_zw 1/7/04 259 261 260 262 ; … … 263 265 ; me semble moins propre. Je la donne en remarque. 264 266 ; 265 IF n_elements (i_profil)GT 2 THEN BEGIN266 ; c1_s (ind) = spline(z_zw(i_profil), c1_z(i_profil), z_s(ind), 1) ; cubic spline267 c1_s (ind) = interpol(c1_z(i_profil), z_zt(i_profil), z_s(ind)) ; linear interpolation268 ; SPLINE_P, z_zw (i_profil), c1_z(i_profil), z_s(ind), c1_s(ind); generalization of spline(*)269 ;print, z_zw (i_profil), c1_z(i_profil), z_s(ind), c1_s(ind)267 IF n_elements[i_profil] GT 2 THEN BEGIN 268 ; c1_s[ind] = spline(z_zw[i_profil], c1_z[i_profil], z_s[ind], 1) ; cubic spline 269 c1_s[ind] = interpol(c1_z[i_profil], z_zt[i_profil], z_s[ind]) ; linear interpolation 270 ; SPLINE_P, z_zw[i_profil], c1_z[i_profil], z_s[ind], c1_s[ind] ; generalization of spline(*) 271 ;print, z_zw[i_profil], c1_z[i_profil], z_s[ind], c1_s[ind] 270 272 ENDIF ELSE BEGIN 271 c1_s (ind) = interpol(c1_z(i_profil), z_zt(i_profil), z_s(ind))273 c1_s[ind] = interpol(c1_z[i_profil], z_zt[i_profil], z_s[ind]) 272 274 ENDELSE 273 275 ; … … 275 277 276 278 IF sig_bowl EQ 1 THEN BEGIN 277 bowl_s = interpol(s_z (i_profil), z_zt(i_profil), sobwlmax(i, j))279 bowl_s = interpol(s_z[i_profil], z_zt[i_profil], sobwlmax(i, j)) 278 280 ENDIF ELSE bowl_s = 0 279 281 … … 281 283 ; 282 284 ; 283 ; x1_s (ind) = (c1_s(ind+1)-c1_s(ind))/(z_s(ind+1)-z_s(ind))285 ; x1_s[ind] = (c1_s(ind+1)-c1_s[ind])/(z_s(ind+1)-z_s[ind]) 284 286 ; x1_s = c1_s 285 287 -
trunk/procs/bining3.pro
r165 r169 67 67 , E3T=e3t $ 68 68 , TMASK=tmask 69 69 ; 70 compile_opt idl2, strictarrsubs 71 ; 70 72 size3d = size(density) 71 jpi = size3d (1)72 jpj = size3d (2)73 jpk = size3d (3)73 jpi = size3d[1] 74 jpj = size3d[2] 75 jpk = size3d[3] 74 76 75 77 s_s = sigma … … 83 85 vol_bin2 = fltarr(N_s) 84 86 85 x1_bin (*, *, *)= !VALUES.F_NAN86 depth_bin (*, *, *)= !VALUES.F_NAN87 thick_bin (*, *, *)= !VALUES.F_NAN88 bowl_bin (*, *)= !VALUES.F_NAN89 vol_bin (*, *, *)=090 vol_bin2 (*)=087 x1_bin[*, *, *] = !VALUES.F_NAN 88 depth_bin[*, *, *] = !VALUES.F_NAN 89 thick_bin[*, *, *] = !VALUES.F_NAN 90 bowl_bin[*, *] = !VALUES.F_NAN 91 vol_bin[*, *, *] =0 92 vol_bin2[*] =0 91 93 92 delta_sig = s_s (1)-s_s(0)94 delta_sig = s_s[1]-s_s[0] 93 95 94 96 vol_tot = 0. 95 97 count_bin = fltarr(jpi, jpj, N_s) 96 98 count_bin2 = fltarr(N_s) 97 count_bin (*, *, *)=098 count_bin2 (*)=099 count_bin[*, *, *] =0 100 count_bin2[*] =0 99 101 100 102 ; perform bining - loop over i,j,k … … 102 104 FOR i = 0, (jpi-1) DO BEGIN 103 105 FOR j = 0, (jpj-1) DO BEGIN 104 i_ocean = where(tmask (i, j, *)EQ 1)106 i_ocean = where(tmask[i, j, *] EQ 1) 105 107 IF i_ocean[0] NE -1 THEN BEGIN ; on n entre que si il y a des points ocean 106 i_bottom = i_ocean (n_elements(i_ocean)-1)108 i_bottom = i_ocean[n_elements(i_ocean)-1] 107 109 FOR k = 0, i_bottom-1 DO BEGIN 108 bin_index = floor((density (i, j, k)-s_s(0))/delta_sig+0.5)110 bin_index = floor((density[i, j, k]-s_s[0])/delta_sig+0.5) 109 111 IF bin_index LT 0 OR bin_index GT (N_s-1) THEN BEGIN 110 ; print, ' WARNING: bin index out of bounds - density (i, j, k), bin_index, N_s = ', density(i, j, k), bin_index, N_s112 ; print, ' WARNING: bin index out of bounds - density[i, j, k], bin_index, N_s = ', density[i, j, k], bin_index, N_s 111 113 ENDIF ELSE BEGIN 112 vol_bin (i, j, bin_index) = vol_bin(i, j, bin_index) + e1t(i,j)*e2t(i,j)*e3t(k)113 vol_bin2 (bin_index) = vol_bin2(bin_index) + e1t(i,j)*e2t(i,j)*e3t(k)114 vol_tot = vol_tot + e1t (i,j)*e2t(i,j)*e3t(k)115 count_bin (i, j, bin_index) = count_bin(i, j, bin_index)+1116 count_bin2 (bin_index) = count_bin2(bin_index)+1114 vol_bin[i, j, bin_index] = vol_bin[i, j, bin_index] + e1t[i,j]*e2t[i,j]*e3t[k] 115 vol_bin2[bin_index] = vol_bin2[bin_index] + e1t[i,j]*e2t[i,j]*e3t[k] 116 vol_tot = vol_tot + e1t[i,j]*e2t[i,j]*e3t[k] 117 count_bin[i, j, bin_index] = count_bin[i, j, bin_index]+1 118 count_bin2[bin_index] = count_bin2[bin_index]+1 117 119 ENDELSE 118 120 ENDFOR … … 122 124 123 125 mask_index = where(count_bin EQ 0) 124 vol_bin (mask_index)= !VALUES.F_NAN126 vol_bin[mask_index] = !VALUES.F_NAN 125 127 mask_index2 = where(count_bin2 EQ 0) 126 ; vol_bin2 (mask_index2)= !VALUES.F_NAN128 ; vol_bin2[mask_index2] = !VALUES.F_NAN 127 129 128 130 print, ' total volume of box (m3) = ', vol_tot -
trunk/procs/compute_time.pro
r165 r169 15 15 FUNCTION compute_time, timeave, date1, date2 $ 16 16 , NEXT=next 17 ; 18 compile_opt idl2, strictarrsubs 17 19 ; 18 20 @common -
trunk/procs/cpsw.pro
r165 r169 8 8 ;- 9 9 FUNCTION cpsw, t, s, p 10 ; 11 compile_opt idl2, strictarrsubs 10 12 ; 11 13 CP1 = 0. -
trunk/procs/data_read.pro
r165 r169 12 12 , _EXTRA=extra 13 13 ; 14 compile_opt idl2, strictarrsubs 15 ; 14 16 @common 15 17 @com_eg 16 17 18 ; 18 19 ; same data as previous read ? 19 20 … … 330 331 FOR i = 0, timedim-1 DO BEGIN 331 332 print, 'data_read: Performing monthly bining. Indices: first, last, current = 0,', timedim-1, i 332 sfild2 = {name: sfild.name, data: sfild.data (*, *, *, i), legend: sfild.legend, units: sfild.units, origin: sfild.origin, dim: sfild.dim-1}333 sfild2 = {name: sfild.name, data: sfild.data[*, *, *, i], legend: sfild.legend, units: sfild.units, origin: sfild.origin, dim: sfild.dim-1} 333 334 bin_sigma, cmd, sfild2, BOXZOOM = box_plot, ALL_DATA = all_data 334 fldrd (*, *, *, i)= sfild2.data335 fldrd[*, *, *, i]= sfild2.data 335 336 ENDFOR 336 337 fldr = {name: sfild.name, data:fldrd, legend: sfild.legend, units: sfild.units, origin: sfild.origin, dim: sfild.dim} … … 407 408 print, 'Be careful some data are null. They will be masked ! ' 408 409 idx0 = where(field2.data EQ 0.0) 409 field2.data (idx0)= valmask410 fldr.data (idx0)= valmask411 div (idx0)= valmask410 field2.data[idx0] = valmask 411 fldr.data[idx0] = valmask 412 div[idx0] = valmask 412 413 ENDIF 413 414 414 415 idx = where(field2.data NE valmask) 415 div (idx) = fldr.data(idx) / field2.data(idx)416 div[idx] = fldr.data[idx] / field2.data[idx] 416 417 417 418 fldr.data = div -
trunk/procs/daypm.pro
r165 r169 9 9 FUNCTION daypm, month_i, year_i 10 10 ; 11 compile_opt idl2, strictarrsubs 12 ; 11 13 @common 12 14 @com_eg 13 ;14 15 ; 15 16 IF calendar_type LE 1 THEN BEGIN -
trunk/procs/decode_cmd.pro
r165 r169 8 8 ;- 9 9 FUNCTION decode_cmd, cmdline, index 10 ; 11 compile_opt idl2, strictarrsubs 10 12 ; 11 13 @common -
trunk/procs/def_box.pro
r165 r169 9 9 FUNCTION def_box, plt, dimplot, boxdef, stride 10 10 ; 11 compile_opt idl2, strictarrsubs 12 ; 11 13 @common 12 14 @com_eg 13 ;14 15 ; 15 16 IF debug_w THEN print, ' ENTER def_box...' -
trunk/procs/def_dbase.pro
r165 r169 8 8 ;- 9 9 FUNCTION def_dbase, expid 10 ; 11 compile_opt idl2, strictarrsubs 10 12 ; 11 13 @common -
trunk/procs/def_dptyp.pro
r165 r169 8 8 ;- 9 9 FUNCTION def_dptyp, cmd 10 ; 11 compile_opt idl2, strictarrsubs 10 12 ; 11 13 @common -
trunk/procs/def_file_name.pro
r165 r169 9 9 ;- 10 10 PRO def_file_name, cmd, ncdf_db, file_name, delta_t1 11 ; 12 compile_opt idl2, strictarrsubs 11 13 ; 12 14 @common -
trunk/procs/def_grid.pro
r165 r169 8 8 ;- 9 9 PRO def_grid, cmd 10 ; 11 compile_opt idl2, strictarrsubs 10 12 ; 11 13 @common -
trunk/procs/def_month.pro
r165 r169 8 8 ;- 9 9 FUNCTION def_month, timave, date 10 ; 11 compile_opt idl2, strictarrsubs 10 12 ; 11 13 @common -
trunk/procs/def_time_origin.pro
r165 r169 11 11 ;- 12 12 FUNCTION def_time_origin, timeave, date1 13 13 ; 14 compile_opt idl2, strictarrsubs 15 ; 14 16 len1 = strlen(date1) 15 17 -
trunk/procs/def_work.pro
r165 r169 14 14 ;- 15 15 PRO def_work, data_base_list, out_ps, cmdline, out_all, other_file, spec_base_list 16 ; 17 compile_opt idl2, strictarrsubs 16 18 ; 17 19 @common -
trunk/procs/densit_pltmap_plot.pro
r165 r169 9 9 10 10 sizefld=size(fld) 11 fld=fld (*,*,round((sig_surface-sig_min)/sig_del),*)12 fld = reform(fld, sizefld (1), sizefld(2), sizefld(4))13 fld2 = fltarr(sizefld (4))14 ;FOR count = 0, sizefld (4)-1 DO fld2(count) = mean(fld(*,*,count), /nan)15 FOR count = 0, sizefld (4)-1 DO fld2(count) = stddev(fld(*,*,count), /nan)11 fld=fld[*,*,round((sig_surface-sig_min)/sig_del),*] 12 fld = reform(fld, sizefld[1], sizefld[2], sizefld[4]) 13 fld2 = fltarr(sizefld[4]) 14 ;FOR count = 0, sizefld[4]-1 DO fld2[count] = mean(fld[*,*,count], /nan) 15 FOR count = 0, sizefld[4]-1 DO fld2[count] = stddev(fld[*,*,count], /nan) 16 16 fld = fld2 17 17 mask_z, fld, cmd, boite_plt1d, dimplot, legz -
trunk/procs/densit_pltmap_read.pro
r165 r169 35 35 cumulative = cumulative + fields.data 36 36 ENDFOR 37 output (*, *, *, year_s-uint(start_year))= cumulative/12.37 output[*, *, *, year_s-uint(start_year)] = cumulative/12. 38 38 ENDFOR 39 39 ; redefine cmd.* to be their initial values before the monthly loop -
trunk/procs/eos.pro
r165 r169 8 8 ;- 9 9 FUNCTION eos, t, s 10 10 ; 11 compile_opt idl2, strictarrsubs 12 ; 11 13 sr=sqrt(abs(s)) 12 14 r1=((((6.536332E-9*t-1.120083E-6)*t+1.001685E-4)*t $ -
trunk/procs/extract_str.pro
r165 r169 8 8 ;- 9 9 FUNCTION extract_str, line, char1, char2 10 ; 11 compile_opt idl2, strictarrsubs 10 12 ; 11 13 pos_char1 = strpos(line,char1) -
trunk/procs/find_jptmax.pro
r165 r169 8 8 , IODIRECTORY=iodirectory $ 9 9 , _EXTRA=extra 10 ; 11 compile_opt idl2, strictarrsubs 10 12 ; 11 13 @common -
trunk/procs/fld_pltext.pro
r165 r169 24 24 FUNCTION fld_pltext, var_name, plttyp, dimplot, hotyp 25 25 ; 26 compile_opt idl2, strictarrsubs 27 ; 26 28 @common 27 29 @com_eg … … 46 48 print, ' *** no min & max found for field ', var_name, $ 47 49 ' in file fld_glo_mmx.def' 48 fldextr.min = min(field.data (where(field.data NE valmask)))49 fldextr.max = max(field.data (where(field.data NE valmask)))50 fldextr.min = min(field.data[where(field.data NE valmask)]) 51 fldextr.max = max(field.data[where(field.data NE valmask)]) 50 52 fldextr.homin = fldextr.min 51 53 fldextr.homax = fldextr.max -
trunk/procs/fld_pltint.pro
r165 r169 24 24 ;- 25 25 FUNCTION fld_pltint, var_name, plttyp, dimplot, hotyp 26 ; 27 compile_opt idl2, strictarrsubs 26 28 ; 27 29 @common -
trunk/procs/get_file.pro
r165 r169 8 8 ;- 9 9 PRO get_file, file_name, ncdf_db 10 ; 11 compile_opt idl2, strictarrsubs 10 12 ; 11 13 @common -
trunk/procs/gfunct.pro
r165 r169 23 23 , pder 24 24 25 ; 26 compile_opt idl2, strictarrsubs 27 ; 25 28 F = A[0] * EXP(A[1] * X) + A[2] 26 29 pder = [[EXP(A[1] * X)], [A[0] * X * EXP(A[1] * X)], [replicate(1.0, N_ELEMENTS(X))]] -
trunk/procs/grep.pro
r165 r169 9 9 FUNCTION grep, command, separator, index 10 10 ; 11 compile_opt idl2, strictarrsubs 12 ; 11 13 spawn, command, line 12 14 line = strcompress(strtrim(line[0])) -
trunk/procs/grepsed.pro
r164 r169 9 9 10 10 FUNCTION grepsed, command1,command2 11 11 ; 12 compile_opt idl2, strictarrsubs 13 ; 12 14 spawn, command1+'|'+command2, line 13 15 -
trunk/procs/lec_pal_gmt.pro
r165 r169 10 10 ;- 11 11 PRO lec_pal_gmt, var_name, c_anot_str, fmt, found, readpal 12 ; 13 compile_opt idl2, strictarrsubs 12 14 ; 13 15 @common … … 70 72 71 73 spawn, 'grep piso '+file_cpt, line 72 IF n_elements(line) NE 0 AND (strpos(line, 'piso')) (0)NE -1 THEN BEGIN74 IF n_elements(line) NE 0 AND (strpos(line, 'piso'))[0] NE -1 THEN BEGIN 73 75 iso = abs(float(line[0])) 74 76 IF iso EQ 0.1 THEN BEGIN -
trunk/procs/macro_read.pro
r165 r169 13 13 , ALL_DATA=all_data $ 14 14 , _EXTRA=extra 15 ; 16 compile_opt idl2, strictarrsubs 15 17 ; 16 18 @common -
trunk/procs/macros/make_anomaly.pro
r168 r169 14 14 , ALL_DATA=all_data $ 15 15 , ZMTYP=zmtyp 16 ; 17 compile_opt idl2, strictarrsubs 16 18 ; 17 19 @common -
trunk/procs/macros/make_bsf.pro
r165 r169 13 13 , ALL_DATA=all_data $ 14 14 , ZMTYP=zmtyp 15 ; 16 compile_opt idl2, strictarrsubs 15 17 ; 16 18 @common … … 74 76 FOR t = 0, jpt DO BEGIN 75 77 76 bsf= diag_bsf(u.data (*, *, *, t), v.data(*, *, *, t))78 bsf= diag_bsf(u.data[*, *, *, t], v.data[*, *, *, t]) 77 79 bsft = bsf[1:jpi-2, 0:jpj-2] 78 bsfr (*, *, t)= shift(bsft, key_shift_old, 0)80 bsfr[*, *, t] = shift(bsft, key_shift_old, 0) 79 81 ENDFOR 80 82 … … 100 102 101 103 isign=where(glamt gt 380.) 102 IF isign (0) NE -1 THEN glamt(isign) = glamt(isign)-360.104 IF isign[0] NE -1 THEN glamt[isign] = glamt[isign]-360. 103 105 isign=where(glamu gt 380.) 104 IF isign (0) NE -1 THEN glamu(isign) = glamu(isign)-360.106 IF isign[0] NE -1 THEN glamu[isign] = glamu[isign]-360. 105 107 isign=where(glamv gt 380.) 106 IF isign (0) NE -1 THEN glamv(isign) = glamv(isign)-360.108 IF isign[0] NE -1 THEN glamv[isign] = glamv[isign]-360. 107 109 isign=where(glamf gt 380.) 108 IF isign (0) NE -1 THEN glamf(isign) = glamf(isign)-360.110 IF isign[0] NE -1 THEN glamf[isign] = glamf[isign]-360. 109 111 ; domdef, old_boite 110 112 -
trunk/procs/macros/make_correlation.pro
r165 r169 12 12 , TIME_2=time_2 $ 13 13 , ALL_DATA=all_data 14 ; 15 compile_opt idl2, strictarrsubs 14 16 ; 15 17 @common … … 51 53 IF debug_w THEN print, ' month idx/value: ', imth, strd(imth) 52 54 53 data1 = (var1.data)[*, *, reform(idxm (imth,*), njpt)]54 data2 = (var2.data)[*, *, reform(idxm (imth,*), njpt)]55 data1 = (var1.data)[*, *, reform(idxm[imth,*], njpt)] 56 data2 = (var2.data)[*, *, reform(idxm[imth,*], njpt)] 55 57 nval = njpt 56 58 … … 65 67 pt_correl = pt_correl/float(nmth) 66 68 67 IF idm[0] NE -1 THEN pt_correl (idm)= valmask69 IF idm[0] NE -1 THEN pt_correl[idm] = valmask 68 70 69 71 varname = varname+' '+ntxt -
trunk/procs/macros/make_correldomain.pro
r165 r169 13 13 , TIME_2=time_2 $ 14 14 , ALL_DATA=all_data 15 ; 16 compile_opt idl2, strictarrsubs 15 17 ; 16 18 @common … … 62 64 IF debug_w THEN print, ' month idx/value: ', imth, strd(imth) 63 65 64 data = (var.data)[*, *, reform(idxm (imth,*), njpt)]65 var_domc = (var_dom)[reform(idxm (imth,*), njpt)]66 data = (var.data)[*, *, reform(idxm[imth,*], njpt)] 67 var_domc = (var_dom)[reform(idxm[imth,*], njpt)] 66 68 var_domc = var_domc-mean(var_domc) 67 69 nval = njpt … … 70 72 FOR idx = 0, nxa-1 DO FOR idy = 0, nya-1 DO BEGIN 71 73 72 IF data (idx, idy, 0)NE valmask THEN BEGIN73 x1i = reform(data (idx, idy, *), nval)74 x1 = x1i-mean (x1i)75 zcor = c_timecorrelate (x1, var_domc)76 correl (idx, idy)= correl(idx, idy)+ zcor74 IF data[idx, idy, 0] NE valmask THEN BEGIN 75 x1i = reform(data[idx, idy, *], nval) 76 x1 = x1i-mean[x1i] 77 zcor = c_timecorrelate[x1, var_domc] 78 correl[idx, idy]= correl[idx, idy] + zcor 77 79 ENDIF ELSE BEGIN 78 correl (idx, idy)= valmask80 correl[idx, idy] = valmask 79 81 ENDELSE 80 82 ENDFOR … … 85 87 correl = correl/float(nmth) 86 88 87 IF idm[0] NE -1 THEN correl (idm)= valmask89 IF idm[0] NE -1 THEN correl[idm] = valmask 88 90 89 91 varname = varname+' '+ntxt -
trunk/procs/macros/make_crtm.pro
r168 r169 14 14 , ALL_DATA=all_data $ 15 15 , ZMTYP=zmtyp 16 ; 17 compile_opt idl2, strictarrsubs 16 18 ; 17 19 @common -
trunk/procs/macros/make_curltau.pro
r165 r169 14 14 , ALL_DATA=all_data $ 15 15 , ZMTYP=zmtyp 16 ; 17 compile_opt idl2, strictarrsubs 16 18 ; 17 19 @common -
trunk/procs/macros/make_denflw.pro
r165 r169 11 11 , ALL_DATA=all_data $ 12 12 , ZMTYP=zmtyp 13 ; 14 compile_opt idl2, strictarrsubs 13 15 ; 14 16 @common -
trunk/procs/macros/make_denflx.pro
r168 r169 8 8 , TIME_1=time_1 $ 9 9 , TIME_2=time_2 10 ; 11 compile_opt idl2, strictarrsubs 10 12 ; 11 13 @common -
trunk/procs/macros/make_depth.pro
r165 r169 14 14 , ZMTYP=zmtyp $ 15 15 , _EXTRA=extra 16 ; 17 compile_opt idl2, strictarrsubs 16 18 ; 17 19 @common -
trunk/procs/macros/make_energetics.pro
r168 r169 13 13 , ALL_DATA=all_data $ 14 14 , ZMTYP=zmtyp 15 ; 16 compile_opt idl2, strictarrsubs 15 17 ; 16 18 @common … … 126 128 valmask = 1.e20 127 129 idxw = where (tmaskr EQ 0) 128 IF idxw[0] NE -1 THEN w (idxw)= valmask130 IF idxw[0] NE -1 THEN w[idxw] = valmask 129 131 idxt = where (tmaskr EQ 0) 130 IF idxt[0] NE -1 THEN t (idxt)= 0.131 IF idxt[0] NE -1 THEN s (idxt)= 0.132 IF idxt[0] NE -1 THEN t [idxt] = 0. 133 IF idxt[0] NE -1 THEN s [idxt] = 0. 132 134 tmaskr = 0 ; free memory 133 135 ENDIF ELSE BEGIN 134 136 idxt = where (t EQ valmask) 135 IF idxt (0) NE -1 THEN t (idxt)= 0.136 IF idxt (0) NE -1 THEN s (idxt)= 0.137 IF idxt[0] NE -1 THEN t [idxt] = 0. 138 IF idxt[0] NE -1 THEN s [idxt] = 0. 137 139 ENDELSE 138 140 … … 140 142 maskw = w LT valmask/10. 141 143 w_T = 0.5*( w*maskw + shift(w, 0, 0, -1, 0)*shift(maskw, 0, 0, -1, 0) ) 142 w_T (*, *, (size(w))(3)-1, *) = w_T(*, *, (size(w))(3)-2, *)144 w_T[*, *, (size(w))[3]-1, *] = w_T[*, *, (size(w))[3]-2, *] 143 145 w = 0 144 146 maskw = 0 ; free memory … … 152 154 maskw = w LT valmask/10. 153 155 w_T = 0.5*( w*maskw + shift(w, 0, 0, -1, 0)*shift(maskw, 0, 0, -1, 0) ) 154 w_T (*, *, (size(w))(3)-1, *) = w_T(*, *, (size(w))(3)-2, *)156 w_T[*, *, (size(w))[3]-1, *] = w_T[*, *, (size(w))[3]-2, *] 155 157 w = 0 156 w_T (idx_t)= valmask158 w_T[idx_t] = valmask 157 159 maskw = 0 ; free memory 158 160 ENDIF ELSE w_T = w 159 161 160 idx_2d = where (u (*, *, 0, 0)GT valmask/10.)161 IF idx_2d (0) NE -1 THEN tx(idx_2d)= valmask162 idx_2d = where (v (*, *, 0, 0)GT valmask/10.)163 IF idx_2d (0) NE -1 THEN ty(idx_2d)= valmask162 idx_2d = where (u[*, *, 0, 0] GT valmask/10.) 163 IF idx_2d[0] NE -1 THEN tx[idx_2d] = valmask 164 idx_2d = where (v[*, *, 0, 0] GT valmask/10.) 165 IF idx_2d[0] NE -1 THEN ty[idx_2d] = valmask 164 166 idx = where (t LT valmask/10.) 165 167 IF t_unit NE "C" THEN BEGIN 166 IF idx (0) NE -1 THEN t(idx) = t(idx)-273.15168 IF idx[0] NE -1 THEN t[idx] = t[idx]-273.15 167 169 ENDIF 168 170 idx1 = n_elements(idx_t) … … 172 174 print, ' *** WARNING wo and thetao do not have same masks !', idx1, idx2, min(idx_t -idx_w), max(idx_t- idx_w) 173 175 print, ' *** Setting wo to zero on all masked points' 174 w_T (idx_w)= 0.176 w_T[idx_w] = 0. 175 177 ENDIF ELSE print, ' check these are 0 (W,T): ', min(idx_t -idx_w), max(idx_t- idx_w) 176 178 w = 0 & idx = 0 & idx_t = 0 & idx_w = 0 ; free memory … … 178 180 idxs=where(s GT valmask/10.) 179 181 IF source_model EQ 'ipcc' THEN print, ' check these are 0 (T,S): ', min(idxt -idxs), max(idxt-idxs) 180 IF idxt[0] NE -1 THEN t (idxt)=0.181 IF idxs[0] NE -1 THEN s (idxs)=0.182 IF idxt[0] NE -1 THEN t[idxt]=0. 183 IF idxs[0] NE -1 THEN s[idxs]=0. 182 184 idxs = 0 ; free memory 183 185 … … 198 200 rhop = ( ( 4.8314E-4*s + r3*sr +r2)*s +r1) 199 201 print, ' rhop min/max', min(rhop), max(rhop) 200 IF idxt[0] NE -1 THEN rhop (idxt)= valmask202 IF idxt[0] NE -1 THEN rhop[idxt] = valmask 201 203 sr = 0 & r1 = 0 & r2 = 0 & r3 = 0 & t = 0 & s = 0 ; free memory 202 204 … … 216 218 rho_diff = (rho_s-shift(rho_s,-1))/shift(e3w, -1) 217 219 rho_diff = shift(temporary(rho_diff), 1) 218 rho_diff (0)= 0.220 rho_diff[0] = 0. 219 221 220 222 ; transform onto T grid 221 223 222 224 rho_diff_T = 0.5*(rho_diff+shift(rho_diff, -1)) 223 rho_diff_T ((size(rho_diff))(1)-1) = rho_diff((size(rho_diff))(1)-2)225 rho_diff_T[(size(rho_diff))[1]-1] = rho_diff[(size(rho_diff))[1]-2] 224 226 225 227 stab_inv = ABS(1./rho_diff_T) … … 241 243 242 244 int_val2 = ((rhop-rho_s4d)^2)*stab_inv 243 IF idxt[0] NE -1 THEN int_val2 (idxt)= 0.245 IF idxt[0] NE -1 THEN int_val2[idxt] = 0. 244 246 245 247 print, ' compute APE...' … … 257 259 258 260 int_val = (rhop-rho_s4d)*(w_T) 259 IF idxt[0] NE -1 THEN int_val (idxt)= 0.261 IF idxt[0] NE -1 THEN int_val[idxt] = 0. 260 262 261 263 ; remove first 2 levels (MXL too unstable) … … 284 286 idyv = where(vmean GT valmask/10.) 285 287 286 IF idx[0] NE -1 THEN tx (idx)= 0.287 IF idy[0] NE -1 THEN ty (idy)= 0.288 IF idxu[0] NE -1 THEN umean (idxu)= 0.289 IF idyv[0] NE -1 THEN vmean (idyv)= 0.288 IF idx[0] NE -1 THEN tx[idx] = 0. 289 IF idy[0] NE -1 THEN ty[idy] = 0. 290 IF idxu[0] NE -1 THEN umean[idxu] = 0. 291 IF idyv[0] NE -1 THEN vmean[idyv] = 0. 290 292 291 293 dot_prodx = tx*umean … … 358 360 359 361 maxdepth = 10 360 IF gdept (0) GT maxdepth THEN maxdepth = gdept(0)362 IF gdept[0] GT maxdepth THEN maxdepth = gdept[0] 361 363 362 364 ;domdef, [210., 280., -5., 5., 0., maxdepth] -
trunk/procs/macros/make_eos.pro
r168 r169 62 62 , ZMTYP=zmtyp 63 63 ; 64 compile_opt idl2, strictarrsubs 65 ; 64 66 @common 65 67 @com_eg … … 77 79 idxt=where(t eq valmask) 78 80 idxs=where(s eq valmask) 79 IF idxt[0] NE -1 THEN t (idxt)=0.80 IF idxs[0] NE -1 THEN s (idxs)=0.81 IF idxt[0] NE -1 THEN t[idxt]=0. 82 IF idxs[0] NE -1 THEN s[idxs]=0. 81 83 ; 82 84 ; potential volumic mass … … 89 91 rhopn = ( ( 4.8314E-4*s + r3*sr +r2)*s +r1) 90 92 91 IF idxs[0] NE -1 THEN rhopn (idxt)= valmask93 IF idxs[0] NE -1 THEN rhopn[idxt] = valmask 92 94 93 95 fdirec = tn.direc … … 108 110 remove_dim = 1 109 111 fdirec = 'xyt' 110 name_suff = ' averaged in '+vert_type+'['+strtrim(string(vert_mean (0)), 2)+','+strtrim(string(vert_mean(1)), 2)+']'112 name_suff = ' averaged in '+vert_type+'['+strtrim(string(vert_mean[0]), 2)+','+strtrim(string(vert_mean[1]), 2)+']' 111 113 flegend = name_suff 112 114 ENDIF -
trunk/procs/macros/make_eos_sal.pro
r168 r169 52 52 , ZMTYP=zmtyp 53 53 ; 54 compile_opt idl2, strictarrsubs 55 ; 54 56 @common 55 57 @com_eg … … 67 69 idxt=where(t eq valmask) 68 70 idxs=where(s eq valmask) 69 IF idxt[0] NE -1 THEN t (idxt)=0.70 IF idxs[0] NE -1 THEN s (idxs)=0.71 IF idxt[0] NE -1 THEN t[idxt]=0. 72 IF idxs[0] NE -1 THEN s[idxs]=0. 71 73 ; 72 74 ; potential volumic mass … … 79 81 rhopn = ( ( 0.0E-4*s + 0.0*r3*sr +r2)*s +0.0*r1) 80 82 81 IF idxs[0] NE -1 THEN rhopn (idxt)= valmask83 IF idxs[0] NE -1 THEN rhopn[idxt] = valmask 82 84 83 85 fdirec = tn.direc … … 98 100 remove_dim = 1 99 101 fdirec = 'xyt' 100 name_suff = ' averaged in '+vert_type+'['+strtrim(string(vert_mean (0)), 2)+','+strtrim(string(vert_mean(1)), 2)+']'102 name_suff = ' averaged in '+vert_type+'['+strtrim(string(vert_mean[0]), 2)+','+strtrim(string(vert_mean[1]), 2)+']' 101 103 flegend = name_suff 102 104 ENDIF -
trunk/procs/macros/make_eos_tem.pro
r168 r169 49 49 , ZMTYP=zmtyp 50 50 ; 51 compile_opt idl2, strictarrsubs 52 ; 51 53 @common 52 54 @com_eg … … 64 66 idxt=where(t eq valmask) 65 67 idxs=where(s eq valmask) 66 IF idxt[0] NE -1 THEN t (idxt)=0.67 IF idxs[0] NE -1 THEN s (idxs)=0.68 IF idxt[0] NE -1 THEN t[idxt]=0. 69 IF idxs[0] NE -1 THEN s[idxs]=0. 68 70 ; 69 71 ; potential volumic mass … … 76 78 rhopn = ( ( 4.8314E-4*s + r3*sr +r2)*s +r1) 77 79 78 IF idxs[0] NE -1 THEN rhopn (idxt)= valmask80 IF idxs[0] NE -1 THEN rhopn[idxt] = valmask 79 81 80 82 fdirec = tn.direc … … 95 97 remove_dim = 1 96 98 fdirec = 'xyt' 97 name_suff = ' averaged in '+vert_type+'['+strtrim(string(vert_mean (0)), 2)+','+strtrim(string(vert_mean(1)), 2)+']'99 name_suff = ' averaged in '+vert_type+'['+strtrim(string(vert_mean[0]), 2)+','+strtrim(string(vert_mean[1]), 2)+']' 98 100 flegend = name_suff 99 101 ENDIF -
trunk/procs/macros/make_grad.pro
r165 r169 13 13 , ALL_DATA=all_data $ 14 14 , ZMTYP=zmtyp 15 ; 16 compile_opt idl2, strictarrsubs 15 17 ; 16 18 @common -
trunk/procs/macros/make_interp.pro
r165 r169 16 16 , ALL_DATA=all_data $ 17 17 , ZMTYP=zmtyp 18 ; 19 compile_opt idl2, strictarrsubs 18 20 ; 19 21 @common … … 56 58 IF( idx[0] EQ -1 ) THEN STOP, 'STOP. Something wrong in the data' 57 59 58 interp (idx) = ( nb EQ 6 ) ? coef1*var1.data(idx) + coef2*var2.data(idx) + coef3*var3.data(idx) : coef1*var1.data(idx) + coef2*var2.data(idx)60 interp[idx] = ( nb EQ 6 ) ? coef1*var1.data[idx] + coef2*var2.data[idx] + coef3*var3.data[idx] : coef1*var1.data[idx] + coef2*var2.data[idx] 59 61 60 62 legend = ( nb EQ 6 ) ? ' : ('+strtrim(string(coef1), 2)+')*'+varname1+'+('+ strtrim(string(coef2), 2)+')*'+varname2+'+('+ strtrim(string(coef3), 2)+')*'+varname3 : ' : ('+strtrim(string(coef1), 2)+')*'+varname1+'+('+ strtrim(string(coef2), 2)+')*'+varname2 … … 72 74 73 75 END 74 75 76 -
trunk/procs/macros/make_linfit.pro
r168 r169 12 12 , TIME_2=time_2 $ 13 13 , ALL_DATA=all_data 14 ; 15 compile_opt idl2, strictarrsubs 14 16 ; 15 17 @common … … 55 57 IF debug_w THEN print, ' month idx/value: ', imth, strd(imth) 56 58 57 data1 = (var1.data)[*, *, reform(idxm (imth,*), njpt)]58 data2 = (var2.data)[*, *, reform(idxm (imth,*), njpt)]59 data1 = (var1.data)[*, *, reform(idxm[imth,*], njpt)] 60 data2 = (var2.data)[*, *, reform(idxm[imth,*], njpt)] 59 61 nval = njpt 60 62 … … 66 68 IF data1(idx, idy, 0) LT valmask/10. THEN BEGIN 67 69 68 x1i = reform(data1 (idx, idy, *), nval)69 x2i = reform(data2 (idx, idy, *), nval)70 x1 = x1i-mean (x1i)71 x2 = x2i-mean (x2i)70 x1i = reform(data1[idx, idy, *], nval) 71 x2i = reform(data2[idx, idy, *], nval) 72 x1 = x1i-mean[x1i] 73 x2 = x2i-mean[x2i] 72 74 73 75 no_compute = 0 … … 81 83 IF (size(idt))[1] GE 3 THEN BEGIN 82 84 IF idt[0] NE -1 THEN BEGIN 83 x1 = x1 (idt)84 x2 = x2 (idt)85 x1 = x1[idt] 86 x2 = x2[idt] 85 87 ENDIF 86 88 ENDIF ELSE BEGIN … … 98 100 99 101 coeffl = linfit(x2, x1, CHISQ = linerrl, PROB = proberrl, SIGMA = sigmaerrl) 100 correl = c_timecorrelate (x2,x1)101 pt_linfit (idx, idy) = pt_linfit(idx, idy)+coeffl(1)102 pt_err (idx, idy) = min(pt_err(idx, idy), proberrl)103 pt_corr (idx, idy) = pt_corr(idx, idy)+ correl102 correl = c_timecorrelate[x2,x1] 103 pt_linfit[idx, idy] = pt_linfit[idx, idy]+coeffl[1] 104 pt_err[idx, idy] = min(pt_err[idx, idy], proberrl) 105 pt_corr[idx, idy] = pt_corr[idx, idy] + correl 104 106 105 107 IF debug_w THEN BEGIN … … 107 109 print, ' idx, idy, linfit_map ', idx, idy, linfit_map 108 110 print, ' size(idt)',size(idt) 109 print, ' pt_linfit, correl = ',coeffl (1), correl111 print, ' pt_linfit, correl = ',coeffl[1], correl 110 112 ENDIF 111 113 ENDIF 112 114 ENDIF ELSE BEGIN 113 pt_linfit (idx, idy)= valmask114 pt_err (idx, idy)= 0.115 pt_corr (idx, idy)= 0.115 pt_linfit[idx, idy] = valmask 116 pt_err[idx, idy] = 0. 117 pt_corr[idx, idy] = 0. 116 118 ENDELSE 117 119 118 120 ENDIF ELSE BEGIN 119 pt_linfit (idx, idy)= valmask120 pt_err (idx, idy)= 0.121 pt_corr (idx, idy)= 0.121 pt_linfit[idx, idy] = valmask 122 pt_err[idx, idy] = 0. 123 pt_corr[idx, idy] = 0. 122 124 ENDELSE 123 125 … … 130 132 pt_linfit = pt_linfit/float(nmth) 131 133 pt_corr = pt_corr/float(nmth) 132 IF idm[0] NE -1 THEN pt_linfit (idm)= valmask134 IF idm[0] NE -1 THEN pt_linfit[idm] = valmask 133 135 134 136 ; only take points where correlation larger than 0.2 135 137 136 138 idc = where(abs(pt_corr) LT corr_thres) 137 IF idc[0] NE -1 THEN pt_linfit (idc)= !values.f_nan139 IF idc[0] NE -1 THEN pt_linfit[idc] = !values.f_nan 138 140 139 141 corr_txt = '[r>'+strtrim(string(corr_thres, format = '(f5.2)'), 2)+']' -
trunk/procs/macros/make_linfitdom.pro
r168 r169 12 12 , TIME_2=time_2 $ 13 13 , ALL_DATA=all_data 14 ; 15 compile_opt idl2, strictarrsubs 14 16 ; 15 17 @common … … 62 64 IF debug_w THEN print, ' month idx/value: ', imth, strd(imth) 63 65 64 data1 = (var1.data)[*, *, reform(idxm (imth,*), njpt)]65 var_domc = (var_dom)[reform(idxm (imth,*), njpt)]66 data1 = (var1.data)[*, *, reform(idxm[imth,*], njpt)] 67 var_domc = (var_dom)[reform(idxm[imth,*], njpt)] 66 68 var_domc = var_domc-mean(var_domc) 67 69 nval = njpt … … 71 73 FOR idx = 0, nxa-1 DO FOR idy = 0, nya-1 DO BEGIN 72 74 73 IF data1 (idx, idy, 0)NE valmask THEN BEGIN75 IF data1[idx, idy, 0] NE valmask THEN BEGIN 74 76 75 x1i = reform(data1 (idx, idy, *), nval)76 x1 = x1i-mean (x1i)77 x1i = reform(data1[idx, idy, *], nval) 78 x1 = x1i-mean[x1i] 77 79 78 80 IF linfit_map NE '' THEN BEGIN … … 82 84 ELSE: idt = where (var_domc GE linfit_sep) 83 85 ENDCASE 84 x1 = x1 (idt)85 var_domc = var_domc (idt)86 x1 = x1[idt] 87 var_domc = var_domc[idt] 86 88 ENDIF 87 89 88 90 coeffl = linfit(var_domc, x1, CHISQ = linerrl, PROB = proberrl, SIGMA = sigmaerrl) 89 correl = c_timecorrelate (var_domc,x1)90 pt_linfit (idx, idy) = pt_linfit(idx, idy)+coeffl(1)91 pt_err (idx, idy) = min(pt_err(idx, idy), proberrl)92 pt_corr (idx, idy) = pt_corr(idx, idy)+ correl91 correl = c_timecorrelate[var_domc,x1] 92 pt_linfit[idx, idy] = pt_linfit[idx, idy]+coeffl[1] 93 pt_err[idx, idy] = min[pt_err[idx, idy], proberrl] 94 pt_corr[idx, idy] = pt_corr[idx, idy] + correl 93 95 94 96 IF debug_w THEN BEGIN … … 96 98 print, ' idx, idy, linfit_map ', idx, idy, linfit_map 97 99 print, ' size(idt)',size(idt) 98 print, ' pt_linfit, correl = ',coeffl (1), correl100 print, ' pt_linfit, correl = ',coeffl[1], correl 99 101 stop 100 102 ENDIF … … 102 104 103 105 ENDIF ELSE BEGIN 104 pt_linfit (idx, idy)= valmask105 pt_err (idx, idy)= 0.106 pt_corr (idx, idy)= 0.106 pt_linfit[idx, idy] = valmask 107 pt_err[idx, idy] = 0. 108 pt_corr[idx, idy] = 0. 107 109 ENDELSE 108 110 … … 115 117 pt_linfit = pt_linfit/float(nmth) 116 118 pt_corr = pt_corr/float(nmth) 117 IF idm[0] NE -1 THEN pt_linfit (idm)= valmask119 IF idm[0] NE -1 THEN pt_linfit[idm] = valmask 118 120 119 121 ; only take points where ABS(correlation) larger than a value 120 122 121 123 idc = where(abs(pt_corr) LT corr_thres) 122 pt_linfit (idc)= !values.f_nan124 pt_linfit[idc] = !values.f_nan 123 125 corr_txt = '[r>'+strtrim(string(corr_thres, format = '(f5.2)'), 2)+']' 124 126 -
trunk/procs/macros/make_msf.pro
r168 r169 14 14 , ZMTYP=zmtyp 15 15 ; 16 compile_opt idl2, strictarrsubs 17 ; 16 18 @common 17 19 @com_eg … … 27 29 28 30 idx = where(v.data EQ valmask) 29 IF idx (0) NE -1 THEN v.data(idx)= 0.031 IF idx[0] NE -1 THEN v.data[idx] = 0.0 30 32 ; 31 33 ; Mask bathymetry ? … … 49 51 IF full_name EQ 'atlantic' OR full_name EQ 'pacific' $ 50 52 OR full_name EQ 'indopacific' OR full_name EQ 'indian' THEN BEGIN 51 idx = where (gphit (diaznl_idx+key_shift, *)LE -32.)53 idx = where (gphit[diaznl_idx+key_shift, *] LE -32.) 52 54 msf[idx, *] = valmask 53 55 ENDIF ELSE BEGIN -
trunk/procs/macros/make_msf2.pro
r165 r169 13 13 , ZMTYP=zmtyp 14 14 ; 15 compile_opt idl2, strictarrsubs 16 ; 15 17 @common 16 18 @com_eg … … 26 28 27 29 idx = where(v.data EQ valmask) 28 IF idx (0) NE -1 THEN v.data(idx)= 0.030 IF idx[0] NE -1 THEN v.data[idx] = 0.0 29 31 ; 30 32 ; Mask bathymetry ? … … 39 41 full_name = strmid(full_name, 5, 50) 40 42 maskzm = maskread(full_name, 'orca', /d3) 41 maskzm = reform((reform(maskzm (*, *, *), jpi*jpj*jpk, /overwrite)#replicate(1, jpt)),jpi, jpj, jpk, jpt, /overwrite)43 maskzm = reform((reform(maskzm[*, *, *], jpi*jpj*jpk, /overwrite)#replicate(1, jpt)),jpi, jpj, jpk, jpt, /overwrite) 42 44 v.data = v.data*maskzm 43 45 ; reduction du champ de vitesse meridienne a la bathymetrique. … … 56 58 IF full_name EQ 'atlantic' OR full_name EQ 'pacific' $ 57 59 OR full_name EQ 'indopacific' OR full_name EQ 'indian' THEN BEGIN 58 idx = where (gphit (diaznl_idx+key_shift, *)LE -32.)60 idx = where (gphit[diaznl_idx+key_shift, *] LE -32.) 59 61 msf[idx, *,*] = valmask 60 62 ENDIF ELSE BEGIN … … 82 84 ;pas de condition : equivalent a msf_mean=0 ---> plot yz a la date T 83 85 IF msf_mean EQ 1 THEN BEGIN 84 field = {name: '', data: msf (*,*,0), legend: '', units: 'm3/s', origin: '', dim: 2, direc:'yz'}86 field = {name: '', data: msf[*,*,0], legend: '', units: 'm3/s', origin: '', dim: 2, direc:'yz'} 85 87 ENDIF 86 88 -
trunk/procs/macros/make_ratio.pro
r165 r169 13 13 , ALL_DATA=all_data $ 14 14 , ZMTYP=zmtyp 15 ; 16 compile_opt idl2, strictarrsubs 15 17 ; 16 18 @common … … 53 55 IF( idx[0] EQ -1 ) THEN STOP, 'STOP. Division cannot be computed' 54 56 55 div (idx) = ( var_a.data(idx) - var_b.data(idx) )/( diff(idx))57 div[idx] = ( var_a.data[idx] - var_b.data[idx] )/( diff[idx] ) 56 58 legend = ' : ( '+a+' - '+b+' ) / ( '+c+' - '+d+' )' 57 59 … … 60 62 IF ( (size(tab_var))[1] EQ 2 ) THEN BEGIN 61 63 62 a = strtrim(tab_var (0), 2)63 b = strtrim(tab_var (1), 2)64 a = strtrim(tab_var[0], 2) 65 b = strtrim(tab_var[1], 2) 64 66 65 67 print, 'a : ', a … … 75 77 IF( idx[0] EQ -1 ) THEN STOP, 'STOP. Division cannot be computed' 76 78 77 div (idx) = ( var_a.data(idx) )/( var_b.data(idx))79 div[idx] = ( var_a.data[idx] )/( var_b.data[idx] ) 78 80 legend = ' : '+a+' / '+b 79 81 … … 92 94 93 95 END 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 -
trunk/procs/macros/make_skewness.pro
r168 r169 13 13 , ALL_DATA=all_data $ 14 14 , ZMTYP=zmtyp 15 ; 16 compile_opt idl2, strictarrsubs 15 17 ; 16 18 @common -
trunk/procs/macros/make_sobwlmax.pro
r165 r169 13 13 , ALL_DATA=all_data $ 14 14 , ZMTYP=zmtyp 15 ; 16 compile_opt idl2, strictarrsubs 15 17 ; 16 18 @common -
trunk/procs/macros/make_stddev.pro
r168 r169 13 13 , ALL_DATA=all_data $ 14 14 , ZMTYP=zmtyp 15 ; 16 compile_opt idl2, strictarrsubs 15 17 ; 16 18 @common … … 44 46 IF debug_w THEN print, ' month idx/value: ', imth, strd(imth) 45 47 46 data = (mfld.data)[*, *, reform(idxm (imth,*), njpt)]48 data = (mfld.data)[*, *, reform(idxm[imth,*], njpt)] 47 49 nval = njpt 48 50 ; compute std dev … … 56 58 stdw = stdw/float(nmth) 57 59 58 IF idm[0] NE -1 THEN stdw (idm)= valmask60 IF idm[0] NE -1 THEN stdw[idm] = valmask 59 61 60 62 varname = varname+' '+ntxt … … 64 66 std_tot = a_timecorrelate(mfld.data, 0,/covariance) 65 67 std_tot = sqrt( jpt/(jpt-1)*std_tot ) 66 IF idm[0] NE -1 THEN std_tot (idm)= valmask68 IF idm[0] NE -1 THEN std_tot[idm] = valmask 67 69 idx_diff = WHERE (stdw NE valmask) 68 stdw (idx_diff) = stdw(idx_diff) - std_tot(idx_diff)70 stdw[idx_diff] = stdw[idx_diff] - std_tot[idx_diff] 69 71 ENDIF 70 72 -
trunk/procs/macros/make_thdepth.pro
r168 r169 12 12 , TIME_2=time_2 $ 13 13 , ALL_DATA=all_data 14 ; 15 compile_opt idl2, strictarrsubs 14 16 ; 15 17 @common … … 43 45 ; depth is a linear interpolation 44 46 45 hd20 = gdept (ik20c)47 hd20 = gdept[ik20c] 46 48 IF debug_w THEN help, hd20 47 49 IF debug_w THEN print, 'min(hd20), max(hd20) : ', min(hd20), max(hd20) 48 50 49 51 IF valmask EQ 0. THEN BEGIN 50 tempm = tmask (firstxt:lastxt, firstyt:lastyt, firstzt:lastzt)52 tempm = tmask[firstxt:lastxt, firstyt:lastyt, firstzt:lastzt] 51 53 idm = where (tempm EQ 0) 52 tn.data (idm)= 1.e2054 tn.data[idm] = 1.e20 53 55 valmask = 1.e20 54 56 tempm = 0 … … 57 59 FOR jj = 0, (size(tn.data))[2]-1 DO BEGIN 58 60 FOR ji = 0, (size(tn.data))[1]-1 DO BEGIN 59 IF tn.data (ji, jj, 0)LT valmask/10. THEN BEGIN60 hd20[ji, jj] = gdept (ik20c[ji, jj])+ $61 ( gdept (ik20c[ji, jj]+1)-gdept(ik20c[ji, jj])) $62 * tmask (ji, jj, ik20c[ji, jj]+1) * ( 20. - tn.data(ji, jj, ik20c[ji, jj])) $63 / ( tn.data (ji, jj, ik20c[ji, jj]+1) - tn.data(ji, jj, ik20c[ji, jj]))61 IF tn.data[ji, jj, 0] LT valmask/10. THEN BEGIN 62 hd20[ji, jj] = gdept[ik20c[ji, jj]] + $ 63 ( gdept[ik20c[ji, jj]+1]-gdept[ik20c[ji, jj]] ) $ 64 * tmask[ji, jj, ik20c[ji, jj]+1] * ( 20. - tn.data[ji, jj, ik20c[ji, jj]] ) $ 65 / ( tn.data[ji, jj, ik20c[ji, jj]+1] - tn.data[ji, jj, ik20c[ji, jj]] ) 64 66 ENDIF ELSE hd20[ji, jj] = !values.f_nan 65 67 ENDFOR … … 72 74 73 75 ; bound by the ocean depth, minimum value, first T-point depth 74 ; hd20 (ji,jj)= min( zd20, fsdepw(ji,jj,mbathy(ji,jj)) )76 ; hd20[ji,jj] = min( zd20, fsdepw(ji,jj,mbathy(ji,jj)) ) 75 77 76 78 field = {name: '', data: hd20, legend: '', units: '', origin: '', dim: 2, direc:'xy'} -
trunk/procs/macros/make_transfo.pro
r168 r169 21 21 , ZMTYP=zmtyp 22 22 ; 23 compile_opt idl2, strictarrsubs 24 ; 23 25 @common 24 26 @com_eg … … 38 40 idxt=where(sst eq valmask) 39 41 idxs=where(sss eq valmask) 40 IF idxt[0] NE -1 THEN sst (idxt)=0.41 IF idxs[0] NE -1 THEN sss (idxs)=0.42 IF idxt[0] NE -1 THEN sst[idxt]=0. 43 IF idxs[0] NE -1 THEN sss[idxs]=0. 42 44 43 45 ; Compute potential density (sigma_0) 44 46 45 rho0 = eos(sst, sss)-1000. 46 IF idxs[0] NE -1 THEN rho0(idxt) = valmask 47 47 rho0 = eos[sst, sss]-1000. 48 IF idxs[0] NE -1 THEN rho0[idxt] = valmask 48 49 49 50 ; Define density grid … … 89 90 90 91 trf_bin = fltarr(n_sig+1) 91 trf_bin (*)= 0.92 trf_bin[*] = 0. 92 93 ; 93 94 ; Loop on 2D domain (and time as '*' dimension) … … 96 97 FOR jj = 0, (jpj-1) DO BEGIN 97 98 98 IF rho0 (ji, jj, 0)NE valmask THEN BEGIN99 IF rho0[ji, jj, 0] NE valmask THEN BEGIN 99 100 ; 100 101 ; input fields 101 t = sst (ji, jj, *)102 s = sss (ji, jj, *)103 hf = hef (ji, jj, *)104 ep = emp (ji, jj, *)105 sig = rho0 (ji, jj, *)102 t = sst[ji, jj, *] 103 s = sss[ji, jj, *] 104 hf = hef[ji, jj, *] 105 ep = emp[ji, jj, *] 106 sig = rho0[ji, jj, *] 106 107 107 108 ; compte alpha, beta, Cp 108 al = alpha (t, s, P)109 be = betar (t, s, P)110 Cp = Cpsw (t, s, P)109 al = alpha[t, s, P] 110 be = betar[t, s, P] 111 Cp = Cpsw[t, s, P] 111 112 112 113 ; find density bin … … 118 119 ; compute area and integral of fluxes 119 120 120 dA = e1t (ji, jj)*e2t(ji, jj)121 dA = e1t[ji, jj]*e2t[ji, jj] 121 122 ; conv = m2 -> 1.e6 km2 122 123 conv = 1.e-6*1.e-6 … … 132 133 133 134 FOR jt = 0, (jpt-1) DO BEGIN 134 area (ib(jt), jt) = area(ib(jt), jt)+ dA135 nval (ib(jt), jt) = nval(ib(jt), jt)+ 1135 area[ib[jt], jt] = area[ib[jt], jt] + dA 136 nval[ib[jt], jt] = nval[ib[jt], jt] + 1 136 137 ENDFOR 137 138 … … 139 140 ; convwf : kg/m2/s=mm/s -> m/s 140 141 convwf = 1.e-3 141 fheat (ji, jj, *)= ( -al/Cp)*hf142 fwafl (ji, jj, *)= (1000.+sig)*be*s*ep*convwf142 fheat[ji, jj, *] = ( -al/Cp)*hf 143 fwafl[ji, jj, *] = (1000.+sig)*be*s*ep*convwf 143 144 144 145 ; volume transport per density class (m3/s) 145 146 ; period averaged density flux 146 147 FOR jt = 0, (jpt-1) DO BEGIN 147 fht = fheat (ji, jj, jt)*dA148 fwf = fwafl (ji, jj, jt)*dA149 hist (ib(jt), jt) = hist(ib(jt), jt)+ (fht+fwf)/drho150 histwf (ib(jt), jt) = histwf(ib(jt), jt)+ fwf/drho148 fht = fheat[ji, jj, jt]*dA 149 fwf = fwafl[ji, jj, jt]*dA 150 hist[ib[jt], jt] = hist[ib[jt], jt] + (fht+fwf)/drho 151 histwf[ib[jt], jt] = histwf[ib[jt], jt] + fwf/drho 151 152 fluxbar = fluxbar + (fht+fwf)*dt 152 153 fluxsal = fluxsal + (fwf)*dt -
trunk/procs/macros/make_wcurl.pro
r168 r169 13 13 , ZMTYP=zmtyp $ 14 14 , ALL_DATA=all_data 15 ; 16 compile_opt idl2, strictarrsubs 15 17 ; 16 18 @common -
trunk/procs/macros/make_wdiv.pro
r165 r169 13 13 , ZMTYP=zmtyp $ 14 14 , ALL_DATA=all_data 15 ; 16 compile_opt idl2, strictarrsubs 15 17 ; 16 18 @common -
trunk/procs/macros/make_wm.pro
r168 r169 13 13 , ALL_DATA=all_data $ 14 14 , ZMTYP=zmtyp 15 ; 16 compile_opt idl2, strictarrsubs 15 17 ; 16 18 @common -
trunk/procs/macros/make_work.pro
r165 r169 11 11 , TIME_2=time_2 $ 12 12 , ZMTYP=zmtyp 13 ; 14 compile_opt idl2, strictarrsubs 13 15 ; 14 16 @common … … 42 44 43 45 work = glamt 44 work (*, *)= 0.46 work[*, *] = 0. 45 47 alpha = glamt 46 alpha (*, *)= 0.48 alpha[*, *] = 0. 47 49 48 50 A = ((glamu-x0)^2)/(a1^2) + ((gphiu-y0)^2)/b1^2 … … 50 52 51 53 idx = where(A LE 1.) 52 work (idx)= 1.54 work[idx] = 1. 53 55 idx = where((A0-A) NE 0) 54 alpha (idx) = (A0(idx)-1)*A(idx)/(A0(idx)-A(idx))56 alpha[idx] = (A0[idx]-1)*A[idx]/(A0[idx]-A[idx]) 55 57 idx = where(A GT 1 AND A0 LE 1) 56 work (idx) = alpha(idx)58 work[idx] = alpha[idx] 57 59 58 VALMASK= 1.e2060 valmask = 1.e20 59 61 60 62 field = {name: '', data: work, legend: '', units: '', origin: '', dim: 0, direc:''} … … 65 67 66 68 return, field 67 END 69 END -
trunk/procs/macros/make_wwv.pro
r165 r169 12 12 , TIME_2=time_2 $ 13 13 , ALL_DATA=all_data 14 ; 15 compile_opt idl2, strictarrsubs 14 16 ; 15 17 @common … … 69 71 FOR jj = 0, nyt-1 DO BEGIN 70 72 FOR ji = 0, nxt-1 DO BEGIN 71 IF tn.data (ji, jj, 0, jt)LT valmask/10. THEN BEGIN72 idx20 = ik20c (ji, jj, jt)73 IF tn.data[ji, jj, 0, jt] LT valmask/10. THEN BEGIN 74 idx20 = ik20c[ji, jj, jt] 73 75 IF idx20 GE 0 THEN BEGIN 74 hd20[ji, jj, jt] = gdept (idx20)+$76 hd20[ji, jj, jt] = gdept[idx20] +$ 75 77 76 ( gdept (idx20+1)- gdept(idx20)) * $ ;DELTAZ entre les niveaux i et i+177 mask (ji, jj, idx20+1) * ( thermo - tn.data(ji, jj, idx20, jt)) $ ;DELTAT entre la hd20 et le niveau i78 / ( tn.data (ji, jj, idx20+1, jt) - tn.data(ji, jj, idx20, jt)) ;variation de T entre les niveau i et i+1(dT/dz)78 ( gdept[idx20+1]- gdept[idx20] ) * $ ;DELTAZ entre les niveaux i et i+1 79 mask[ji, jj, idx20+1] * ( thermo - tn.data[ji, jj, idx20, jt] ) $ ;DELTAT entre la hd20 et le niveau i 80 / ( tn.data[ji, jj, idx20+1, jt] - tn.data[ji, jj, idx20, jt] ) ;variation de T entre les niveau i et i+1(dT/dz) 79 81 ENDIF 80 82 ENDIF ELSE BEGIN hd20[ji, jj, *] = -1. … … 95 97 FOR ji = 0, nxt-1 DO BEGIN 96 98 IF hd20[ji, jj, jt] NE -1. THEN BEGIN 97 wwv (jt) = wwv(jt)+ hd20[ji, jj,jt]*e1t[firstx+ji, firsty+jj]*e2t[firstx+ji, firsty+jj]99 wwv[jt] = wwv[jt] + hd20[ji, jj,jt]*e1t[firstx+ji, firsty+jj]*e2t[firstx+ji, firsty+jj] 98 100 ENDIF 99 101 ENDFOR … … 105 107 ;-------------------------------------------------------------------------------- 106 108 ; bound by the ocean depth, minimum value, first T-point depth 107 ; hd20 (ji,jj) = min( zd20, fsdepw(ji,jj,mbathy(ji,jj)))109 ; hd20[ji,jj] = min( zd20, fsdepw[ji,jj,mbathy[ji,jj]] ) 108 110 field = {name: '', data: wwv, legend: '', units: '', origin: '', dim: 1, direc:'t'} 109 111 -
trunk/procs/macros/make_ytran.pro
r168 r169 14 14 , ZMTYP=zmtyp 15 15 ; 16 compile_opt idl2, strictarrsubs 17 ; 16 18 @common 17 19 @com_eg … … 28 30 29 31 idx = where(v.data EQ valmask) 30 IF idx (0) NE -1 THEN v.data(idx)= 0.032 IF idx[0] NE -1 THEN v.data[idx] = 0.0 31 33 32 34 ; Transport along Y : try = v*e1t*e3t (Sv) -
trunk/procs/macros/make_ztran.pro
r168 r169 13 13 , TIME_2=time_2 $ 14 14 , ZMTYP=zmtyp 15 ; 16 compile_opt idl2, strictarrsubs 15 17 ; 16 18 @common -
trunk/procs/mask_z.pro
r165 r169 10 10 ;- 11 11 PRO mask_z, fld, cmd, boite_pltz, dimplot, legz 12 ; 13 compile_opt idl2, strictarrsubs 12 14 ; 13 15 @common … … 37 39 ; vertical_domain] 38 40 idxmsk=where(fld eq valmask) 39 IF idxmsk (0)ne -1 THEN BEGIN41 IF idxmsk[0] ne -1 THEN BEGIN 40 42 print, ' mask_z: building mask from fld valmask' 41 IF (size(fld)) (0)EQ 3 THEN BEGIN42 zmsk = fld (*,*,1)LT valmask/10.43 IF (size(fld))[0] EQ 3 THEN BEGIN 44 zmsk = fld[*,*,1] LT valmask/10. 43 45 ENDIF ELSE BEGIN 44 46 zmsk = fld LT valmask/10. … … 69 71 ; mask MSF as field=valmask 70 72 idx=where (fld eq valmask) 71 fld (idx)= !values.f_nan73 fld[idx] = !values.f_nan 72 74 boite_pltz = [box_h, vertical_domain] 73 75 ENDELSE … … 100 102 ; mask MSF as field=valmask 101 103 idx=where (fld eq valmask) 102 fld (idx)= !values.f_nan104 fld[idx] = !values.f_nan 103 105 boite_pltz = [box_h, vertical_domain] 104 106 ENDELSE … … 140 142 boite_pltz = [idx[0], idx[0], boite_pltz[2:5]] 141 143 domdef, boite_pltz, /xindex ; indice pour x 142 IF (size(fld)) (0)EQ 3 THEN BEGIN143 tmask[idx[0], firstyt:lastyt, *] = fld (*,*,1)LT valmask/10.144 IF (size(fld))[0] EQ 3 THEN BEGIN 145 tmask[idx[0], firstyt:lastyt, *] = fld[*,*,1] LT valmask/10. 144 146 ENDIF ELSE BEGIN 145 147 tmask[idx[0], firstyt:lastyt, *] = fld LT valmask/10. … … 190 192 zbox = def_box(cmd.plt, dimplot, legz, time_stride) 191 193 boite_pltz = zbox 192 prof1 = gdept (0)193 prof2 = gdept (0)194 prof1 = gdept[0] 195 prof2 = gdept[0] 194 196 IF dimplot EQ 2 THEN BEGIN 195 197 CASE strmid(cmd.plt, 1, 1) OF -
trunk/procs/maskread.pro
r165 r169 24 24 FUNCTION maskread, bathy, model $ 25 25 , D3=d3 26 ; 27 compile_opt idl2, strictarrsubs 26 28 ; 27 29 @common … … 103 105 iformat = string('(', il3+2, 'i3)') 104 106 FOR jj = jpjglo-1, 0, -1 DO BEGIN 105 zbati (0:ifreq-1)= 0107 zbati[0:ifreq-1] = 0 106 108 readf, iunit, FORMAT = iformat, ij, zbati 107 zbat (il1:il2, jj) = zbati(0:il3)109 zbat[il1:il2, jj] = zbati[0:il3] 108 110 END 109 111 il1 = il1 + ifreq … … 118 120 iformat = string('(', il3+2, 'i3)') 119 121 FOR jj = jpjglo-1, 0, -1 DO BEGIN 120 zbati2 (0:irest-1)= 0122 zbati2[0:irest-1] = 0 121 123 readf, iunit, FORMAT = iformat, ij, zbati2 122 zbat (il1:il2, jj) = zbati2(0:il3)124 zbat[il1:il2, jj] = zbati2[0:il3] 123 125 END 124 126 ;; -
trunk/procs/meshes/build_mesh_glosea.pro
r168 r169 13 13 , DELTA_K=delta_k $ 14 14 , ZONAL=zonal 15 15 ; 16 compile_opt idl2, strictarrsubs 17 ; 16 18 @common 17 19 @com_eg … … 188 190 189 191 isign=where(glamt gt 380.) 190 IF isign (0) NE -1 THEN glamt(isign) = glamt(isign)-360.192 IF isign[0] NE -1 THEN glamt[isign] = glamt[isign]-360. 191 193 isign=where(glamu gt 380.) 192 IF isign (0) NE -1 THEN glamu(isign) = glamu(isign)-360.194 IF isign[0] NE -1 THEN glamu[isign] = glamu[isign]-360. 193 195 isign=where(glamv gt 380.) 194 IF isign (0) NE -1 THEN glamv(isign) = glamv(isign)-360.196 IF isign[0] NE -1 THEN glamv[isign] = glamv[isign]-360. 195 197 isign=where(glamf gt 380.) 196 IF isign (0) NE -1 THEN glamf(isign) = glamf(isign)-360.198 IF isign[0] NE -1 THEN glamf[isign] = glamf[isign]-360. 197 199 198 200 ; reduce grid to zonal mean -
trunk/procs/meshes/mesh_from_file.pro
r165 r169 7 7 PRO mesh_from_file, model, file_name, ncdf_db, var_name $ 8 8 , ALL_DATA=all_data 9 ; 10 compile_opt idl2, strictarrsubs 9 11 ; 10 12 @common … … 39 41 idx = where(tmask EQ valmask) 40 42 41 IF idx (0) NE -1 THEN tmask(idx)= 0.43 IF idx[0] NE -1 THEN tmask[idx] = 0. 42 44 idx = where(tmask LE 50.) 43 tmask (idx)= 0.45 tmask[idx] = 0. 44 46 tmask = tmask < 1 45 47 tmask = 1-tmask -
trunk/procs/meshes/mesh_gaussian.pro
r165 r169 14 14 , WHOLE_ARRAYS=whole_arrays $ 15 15 , GLAMBOUNDARY=glamboundary 16 ; 17 compile_opt idl2, strictarrsubs 16 18 ; 17 19 @common -
trunk/procs/meshes/mesh_glosea.pro
r165 r169 15 15 , NO_SHIFT = no_shift $ 16 16 , WHOLE_ARRAYS=whole_arrays 17 17 ; 18 compile_opt idl2, strictarrsubs 19 ; 18 20 @common 19 21 @com_eg … … 137 139 umask = shift(umask,key_shift, 0, 0) 138 140 139 isign=where(glamt LE glamt (0, 0))140 IF isign(0) NE -1 THEN glamt (isign) = glamt(isign)+360.141 isign=where(glamu LE glamu (0, 0))142 IF isign(0) NE -1 THEN glamu (isign) = glamu(isign)+360.141 isign=where(glamt LE glamt[0,0]) 142 IF isign(0) NE -1 THEN glamt[isign] = glamt[isign]+360. 143 isign=where(glamu LE glamu[0,0]) 144 IF isign(0) NE -1 THEN glamu[isign] = glamu[isign]+360. 143 145 144 146 glamv = glamu -
trunk/procs/meshes/mesh_micom.pro
r165 r169 17 17 , H_CONFIG=h_config $ 18 18 , V_CONFIG=V_config 19 19 ; 20 compile_opt idl2, strictarrsubs 21 ; 20 22 @common 21 23 @com_eg -
trunk/procs/meshes/mesh_orca.pro
r165 r169 17 17 , H_CONFIG=h_config $ 18 18 , V_CONFIG=V_config 19 19 ; 20 compile_opt idl2, strictarrsubs 21 ; 20 22 @common 21 23 @com_eg -
trunk/procs/meshes/mesh_pcmdi.pro
r165 r169 8 8 , NO_SHIFT=no_shift $ 9 9 , WHOLE_ARRAYS=whole_arrays 10 ; 11 compile_opt idl2, strictarrsubs 10 12 ; 11 13 @common … … 60 62 61 63 idx = where(tmask EQ valmask) 62 IF idx (0) NE -1 THEN tmask(idx)= 0.64 IF idx[0] NE -1 THEN tmask[idx] = 0. 63 65 64 66 idx = where(tmask LE 50.) 65 tmask (idx)= 0.67 tmask[idx] = 0. 66 68 tmask = tmask < 1 67 69 tmask = 1-tmask … … 89 91 90 92 return 91 END 92 93 END -
trunk/procs/meshes/mesh_regular.pro
r165 r169 18 18 , REVERSE_Y=reverse_y $ 19 19 , GLAMBOUNDARY=glamboundary 20 ; 21 compile_opt idl2, strictarrsubs 20 22 ; 21 23 @common -
trunk/procs/meshes/mesh_tao.pro
r165 r169 8 8 , NO_SHIFT=no_shift $ 9 9 , WHOLE_ARRAYS=whole_arrays 10 ; 11 compile_opt idl2, strictarrsubs 10 12 ; 11 13 @common … … 57 59 ybot = (shift(gphit, 0, +1)+gphit)*0.5*!pi/180.0 58 60 59 ytop (*, jpj-1)= !pi/260 ybot (*, 0)= -!pi/261 ytop[*, jpj-1] = !pi/2 62 ybot[*, 0] = -!pi/2 61 63 62 64 e2t = zrad*(ytop-ybot) … … 98 100 glamf = 0.5*(shift(glamt, -1, 0)+glamt) 99 101 gphif = 0.5*(shift(gphit, 0, -1)+gphit) 100 glamf (jpi-1, *) = glamf(jpi-2, *) + (glamf(jpi-2, *)-glamf(jpi-3, *))101 gphif (*, jpj-1) = gphif(*, jpj-2) + (gphif(*, jpj-2)-gphif(*, jpj-3))102 glamf[jpi-1, *] = glamf[jpi-2, *] + (glamf[jpi-2, *]-glamf[jpi-3, *]) 103 gphif[*, jpj-1] = gphif[*, jpj-2] + (gphif[*, jpj-2]-gphif[*, jpj-3]) 102 104 e1f = e1t 103 105 e2f = e2t -
trunk/procs/meshes/mesh_um.pro
r165 r169 9 9 , WHOLE_ARRAYS=whole_arrays $ 10 10 , REVERSE_Y=reverse_y 11 ; 12 compile_opt idl2, strictarrsubs 11 13 ; 12 14 @common … … 87 89 ybot = (shift(gphit, 0, +1)+gphit)*0.5*!pi/180.0 88 90 89 ytop (*, jpj-1)= !pi/290 ybot (*, 0)= -!pi/291 ytop[*, jpj-1] = !pi/2 92 ybot[*, 0] = -!pi/2 91 93 92 94 e2t = zrad*(-ytop+ybot) … … 140 142 e1t = nc_get(s_file, 'umat.srf') 141 143 e2t = e1t 142 e2t (*)= 1.144 e2t[*] = 1. 143 145 s_file = hom_idl+'grids/masks_um_N144_nofrac.nc' 144 146 tmask = nc_get(s_file, 'umat.msk') … … 151 153 e1t = nc_get(s_file, 'umau.srf') 152 154 e2t = e1t 153 e2t (*)= 1.155 e2t[*] = 1. 154 156 s_file = hom_idl+'grids/masks_um_N144_nofrac.nc' 155 157 tmask = nc_get(s_file, 'umau.msk') -
trunk/procs/msf_mult.pro
r165 r169 14 14 PRO msf_mult, v, msf, msfmsk $ 15 15 , MSK=msk 16 ; 17 compile_opt idl2, strictarrsubs 16 18 ; 17 19 @common … … 36 38 ; donc cette procedure est ineffective. 37 39 zvmask = boundperio( (msk +shift(msk, 0, -1) ) < 1 ) 38 zvmask = zvmask (*)#vert40 zvmask = zvmask[*]#vert 39 41 zv = zv*zvmask 40 42 ENDIF … … 46 48 ; calcul du flux 47 49 ; 48 FOR i = 0, jpi-1 DO BEGIN for j= 0, jpj-1 do begin ze1v (i,j,*,*) = replicate(e1v(i,j),jpk*jpt) & endfor &endfor49 FOR k = 0, jpk-1 DO BEGIN ze3v (*,*,k,*)=replicate(e3t(k),jpi*jpj*jpt) & endfor50 FOR i = 0, jpi-1 DO BEGIN for j= 0, jpj-1 do begin ze1v[i,j,*,*] = replicate(e1v[i,j],jpk*jpt) & endfor &endfor 51 FOR k = 0, jpk-1 DO BEGIN ze3v[*,*,k,*]=replicate(e3t[k],jpi*jpj*jpt) & endfor 50 52 ; 51 53 z= -v*ze1v*ze3v … … 58 60 ; calcul de la msf en integrant depuis le fond 59 61 ; 60 FOR k = jpk-2, 0, -1 DO begin msf[*, k,*] = msf[*, k+1,*]-fm (*, k,*)& endfor62 FOR k = jpk-2, 0, -1 DO begin msf[*, k,*] = msf[*, k+1,*]-fm[*, k,*] & endfor 61 63 ; 62 64 ; msfmsk est le masque associe a msf (utilise pour les graphiques) -
trunk/procs/msf_simple.pro
r165 r169 15 15 PRO msf_simple, v, msf, msfmsk $ 16 16 , MSK=msk 17 ; 18 compile_opt idl2, strictarrsubs 17 19 ; 18 20 @common … … 35 37 IF n_elements(msk) NE 0 THEN BEGIN 36 38 zvmask = boundperio( (msk +shift(msk, 0, -1) ) < 1 ) 37 zvmask = zvmask (*)#vert39 zvmask = zvmask[*]#vert 38 40 zv = zv*zvmask 39 41 ENDIF … … 45 47 ; calcul du flux 46 48 ; 47 FOR i = 0, jpi-1 DO BEGIN for j= 0, jpj-1 do begin ze1v (i,j,*) = replicate(e1v(i,j),jpk) & endfor &endfor48 FOR k = 0, jpk-1 DO BEGIN ze3v (*,*,k)=replicate(e3t(k),jpi*jpj) & endfor49 FOR i = 0, jpi-1 DO BEGIN for j= 0, jpj-1 do begin ze1v[i,j,*] = replicate(e1v[i,j],jpk) & endfor &endfor 50 FOR k = 0, jpk-1 DO BEGIN ze3v[*,*,k]=replicate(e3t[k],jpi*jpj) & endfor 49 51 ; 50 52 z= -v*ze1v*ze3v … … 56 58 ; calcul de la msf en integrant depuis le fond 57 59 ; 58 FOR k = jpk-2, 0, -1 DO begin msf[*, k] = msf[*, k+1]-fm (*, k)& endfor60 FOR k = jpk-2, 0, -1 DO begin msf[*, k] = msf[*, k+1]-fm[*, k] & endfor 59 61 ; 60 62 ; msfmsk est le masque associe a msf (utilise pour les graphiques) -
trunk/procs/mth_decode.pro
r165 r169 40 40 41 41 FOR imth = 0, nmth-1 DO BEGIN 42 idxm (imth, *)= (lindgen(jpt))[strd(imth)-1:*:12]42 idxm[imth, *]= (lindgen(jpt))[strd[imth]-1:*:12] 43 43 ENDFOR 44 44 -
trunk/procs/nc_build.pro
r165 r169 6 6 ;- 7 7 FUNCTION nc_build, cmd, fld, direc, grid_type 8 ; 9 compile_opt idl2, strictarrsubs 8 10 ; 9 11 @com_eg -
trunk/procs/nc_get.pro
r165 r169 6 6 ;- 7 7 FUNCTION nc_get, file_name, var_name 8 ; 9 compile_opt idl2, strictarrsubs 8 10 ; 9 11 ; ouverture fichier netCDF + contenu ? -
trunk/procs/nc_put.pro
r165 r169 12 12 , T=t $ 13 13 , GLOBAL=global 14 14 ; 15 compile_opt idl2, strictarrsubs 16 ; 15 17 @common 16 18 @com_eg … … 178 180 NCDF_ATTPUT, cdfid, fldid, 'short_name', varn 179 181 NCDF_ATTPUT, cdfid, fldid, 'missing_value', fld.missing_value 180 NCDF_ATTPUT, cdfid, fldid, 'valid_min', min(fld.data (where (fld.data NE fld.missing_value))), /float181 NCDF_ATTPUT, cdfid, fldid, 'valid_max', max(fld.data (where (fld.data NE fld.missing_value))), /float182 NCDF_ATTPUT, cdfid, fldid, 'valid_min', min(fld.data[where (fld.data NE fld.missing_value)]), /float 183 NCDF_ATTPUT, cdfid, fldid, 'valid_max', max(fld.data[where (fld.data NE fld.missing_value)]), /float 182 184 183 185 … … 223 225 print, ' Field legend = ',fld.long_name 224 226 print, ' Field dimensions = ',size(fld.data) 225 print, ' Field min/max/unit = ',min(fld.data (where (fld.data NE fld.missing_value))), max(fld.data(where (fld.data NE fld.missing_value))), ' ', fld.units227 print, ' Field min/max/unit = ',min(fld.data[where (fld.data NE fld.missing_value)]), max(fld.data[where (fld.data NE fld.missing_value)]), ' ', fld.units 226 228 print, ' Global attributes= Conventions: ', global.conventions 227 229 print, ' Title: ', global.title -
trunk/procs/nc_read.pro
r165 r169 24 24 , ALL_DATA=all_data $ 25 25 , NO_MEAN=no_mean 26 ; 27 compile_opt idl2, strictarrsubs 26 28 ; 27 29 @common … … 116 118 print, ' *** Warning valmask is negative - changing sign: ', valmask 117 119 idx_t = where (field.data EQ valmask) 118 IF idx_t[0] NE -1 THEN field.data (idx_t)= -valmask120 IF idx_t[0] NE -1 THEN field.data[idx_t] = -valmask 119 121 valmask = -valmask 120 122 idx_t=0 ; free memory … … 136 138 idx_valmsk = (where (field.data EQ valmask)) 137 139 idx_novalmsk = (where (field.data NE valmask)) 138 minf = idx_valmsk[0] EQ -1 ? min(field.data) : min(field.data (idx_novalmsk))139 maxf = idx_valmsk[0] EQ -1 ? max(field.data) : max(field.data (idx_novalmsk))140 minf = idx_valmsk[0] EQ -1 ? min(field.data) : min(field.data[idx_novalmsk]) 141 maxf = idx_valmsk[0] EQ -1 ? max(field.data) : max(field.data[idx_novalmsk]) 140 142 141 143 print, ' = ',field.legend, ' [min/max = ',minf , maxf,' ', field.units,' - masked values = ',valmask, chardim, ']' -
trunk/procs/new_filename.pro
r165 r169 8 8 ;- 9 9 FUNCTION new_filename, file_in, oldgrid, newgrid 10 10 ; 11 compile_opt idl2, strictarrsubs 12 ; 11 13 @common 12 14 @com_eg -
trunk/procs/npc.pro
r165 r169 25 25 FUNCTION npc, s3d 26 26 ; 27 compile_opt idl2, strictarrsubs 28 ; 27 29 @common 28 30 ;; … … 48 50 49 51 ; Indices des points T dans l ocean 50 i_ocean = where(tmask (i,j,*)EQ 1)52 i_ocean = where(tmask[i,j,*] EQ 1) 51 53 52 54 IF (i_ocean[0] NE -1) THEN BEGIN ; on n'entre que si il y a des points ocean 53 55 ; 54 56 ; density profil 55 s_z (*)= s3d(i,j,*)56 ;s_z (*) = rho(i,j,*)57 s_z[*)= s3d[i,j,*) 58 ;s_z[*) = rho[i,j,*) 57 59 ; 58 60 ; 1. Static instability pointer 59 61 ; ----------------------------- 60 62 ; 61 ds =(shift(s_z,-1)-s_z) (i_ocean(0:n_elements(i_ocean)-2))63 ds =(shift(s_z,-1)-s_z)[i_ocean[0:n_elements(i_ocean)-2]] 62 64 ind_c = where(ds LT 0.) 63 65 ; 64 66 ; 2. Vertical mixing for each instable portion of the density profil 65 67 ; 66 IF ( ind_c (0)NE -1 ) THEN BEGIN68 IF ( ind_c[0] NE -1 ) THEN BEGIN 67 69 ncompt=ncompt+1 68 70 ; print, 'static instability at i,j=', i,j … … 73 75 ; vertical iteration 74 76 jiter = 0 75 WHILE ( (ind_c (0)NE -1) AND (jiter LT jpk-1) ) DO BEGIN77 WHILE ( (ind_c[0] NE -1) AND (jiter LT jpk-1) ) DO BEGIN 76 78 jiter = jiter+1 77 79 ; ikup : the first static instability from the sea surface 78 ikup = ind_c (0)80 ikup = ind_c[0] 79 81 ; the density profil is instable below ikup 80 82 ; ikdown : bottom of the instable portion of the density profil 81 83 ; search of ikdown and vertical mixing from ikup to ikdown 82 ze3tot= e3t (ikup)83 zraua = s_z (ikup)84 ze3tot= e3t[ikup] 85 zraua = s_z[ikup] 84 86 jkdown = ikup+1 85 87 ; 86 WHILE (jkdown LE ikbot AND zraua GT s_z (jkdown)) DO BEGIN87 ze3dwn = e3t (jkdown)88 WHILE (jkdown LE ikbot AND zraua GT s_z[jkdown] ) DO BEGIN 89 ze3dwn = e3t[jkdown] 88 90 ze3tot = ze3tot+ze3dwn 89 zraua = ( zraua*(ze3tot-ze3dwn) + s_z (jkdown)*ze3dwn )/ze3tot91 zraua = ( zraua*(ze3tot-ze3dwn) + s_z[jkdown]*ze3dwn )/ze3tot 90 92 jkdown = jkdown + 1 91 93 ENDWHILE 92 94 93 95 FOR jkp = ikup,jkdown-1 DO BEGIN 94 s_z (jkp)= zraua96 s_z[jpk] = zraua 95 97 ENDFOR 96 ds =(shift(s_z, -1)-s_z) (i_ocean(0:n_elements(i_ocean)-2))98 ds =(shift(s_z, -1)-s_z)[i_ocean[0:n_elements(i_ocean)-2]] 97 99 ind_c = where(ds LT 0.) 98 100 ENDWHILE 99 101 ENDIF 100 102 ; save the modifications 101 rhos (i,j,*) = s_z(*)103 rhos[i,j,*] = s_z[*] 102 104 ; 103 105 ; <<-- no more static instability on slab jj -
trunk/procs/overlay_type.pro
r165 r169 8 8 ;- 9 9 FUNCTION overlay_type, iover, dimplot 10 ; 11 compile_opt idl2, strictarrsubs 10 12 ; 11 13 @com_eg -
trunk/procs/plt_map.pro
r165 r169 20 20 ;- 21 21 PRO plt_map, cmd, iplot, win, iover, landscape 22 ; 23 compile_opt idl2, strictarrsubs 22 24 ; 23 25 @common … … 213 215 idx = where(field.data EQ valmask) 214 216 field.data = field.data > (-1.8) 215 IF idx[0] NE -1 THEN field.data (idx)= valmask217 IF idx[0] NE -1 THEN field.data[idx] = valmask 216 218 ENDIF 217 219 … … 828 830 IF strpos(cmd.plt, '@z') GT 1 THEN BEGIN 829 831 fldzm = moyenne(fld, 'x', boite=[20,380,-90,90], /NAN) 830 fldzm = replicate(1, (size(fld)) (1))#fldzm832 fldzm = replicate(1, (size(fld))[1])#fldzm 831 833 fld = fld-fldzm 832 834 edd_txt = ' Eddy ' … … 1147 1149 ENDIF ELSE BEGIN 1148 1150 time_int = time-shift(time, 1) 1149 time_int (0)= 0.1151 time_int[0] = 0. 1150 1152 divi = total(time_int^2) 1151 1153 rms_spec = sqrt( total( ((fld-spec_prev)*time_int)^2) / (divi > 0.)) -
trunk/procs/saxo_mods/pltsc.pro
r165 r169 15 15 , FRACTION=fraction $ 16 16 , _EXTRA=extra 17 ; 18 compile_opt idl2, strictarrsubs 17 19 ; 18 20 @common -
trunk/procs/saxo_mods/plttg.pro
r165 r169 252 252 , _EXTRA=extra 253 253 ; 254 compile_opt idl2, strictarrsubs 255 ; 254 256 @common 255 257 -
trunk/procs/search_time_file.pro
r165 r169 10 10 ;- 11 11 PRO search_time_file, cmd, ncdf_db, suffix, date1, date2, delta_t1, timavef 12 ; 13 compile_opt idl2, strictarrsubs 12 14 ; 13 15 @common -
trunk/procs/set_ps_devices.pro
r165 r169 9 9 PRO set_ps_devices 10 10 ; 11 compile_opt idl2, strictarrsubs 11 12 ; font 12 13 !p.font = 0 -
trunk/procs/set_x_devices.pro
r165 r169 9 9 PRO set_x_devices 10 10 ; 11 compile_opt idl2, strictarrsubs 12 ; 11 13 ; font 12 14 !p.font = -1 -
trunk/procs/spectra.pro
r164 r169 58 58 FUNCTION spectra, datin, nyears, tukey, running, full, bootwin 59 59 ; 60 compile_opt idl2, strictarrsubs 61 ; 60 62 ; 0. Initializations 61 63 ; ------------------ … … 100 102 ; 101 103 FOR t = 0, 12-1 DO BEGIN 102 amean (t)= mean(datin(long(findgen(ntime/12)*12+t)))104 amean[t] = mean(datin(long(findgen(ntime/12)*12+t))) 103 105 ENDFOR 104 106 ; … … 126 128 temperr = stat_error(anom, bootwin) 127 129 ; err= +- sqrt(sum(temperr(i)-sdev)**2)/N) 128 err_std (0)= sqrt( total((temperr-sdev)^2)/n_elements(temperr) )129 ; err_std (0)= sqrt((moment(temperr))[1])130 err_std (1)= mean(temperr)130 err_std[0] = sqrt( total((temperr-sdev)^2)/n_elements(temperr) ) 131 ; err_std[0] = sqrt((moment(temperr))[1]) 132 err_std[1] = mean(temperr) 131 133 132 134 ; … … 134 136 ; 135 137 FOR t = 0, ntime-(running*12) DO BEGIN 136 run_std (t)= sqrt((moment(datad[t:t+(running*12)-1]))[1])138 run_std[t] = sqrt((moment(datad[t:t+(running*12)-1]))[1]) 137 139 ; dattmp = datin[t:t+(running*12)-1] 138 140 ; ntsc = running*12 139 141 ; FOR tsc = 0, 12-1 DO BEGIN 140 ; ameansc (tsc)= mean(dattmp(long(findgen(ntsc/12)*12+tsc)))142 ; ameansc[tsc] = mean(dattmp(long(findgen(ntsc/12)*12+tsc))) 141 143 ; ENDFOR 142 ; sc_run (t)= sqrt((moment(ameansc))[1])143 ; sc_run (t)= max(ameansc)-min(ameansc)144 ; sc_run (t)= mean(dattmp(long(findgen(running/12)*12+ (t MOD 12))))144 ; sc_run[t] = sqrt((moment(ameansc))[1]) 145 ; sc_run[t] = max(ameansc)-min(ameansc) 146 ; sc_run[t] = mean(dattmp(long(findgen(running/12)*12+ (t MOD 12)))) 145 147 ENDFOR 146 148 147 149 ; anomaly time serie wrt time varying SC 148 150 149 ; anom_stdsc = datin - [sc_run (0)#replicate(1, long(running*12/2)), sc_run, sc_run(ntime-(running*12))#replicate(1, long(running*12/2-1))]151 ; anom_stdsc = datin - [sc_run[0]#replicate(1, long(running*12/2)), sc_run, sc_run(ntime-(running*12))#replicate(1, long(running*12/2-1))] 150 152 151 153 ; run_std = smooth(run_std, running*12) … … 156 158 ; 157 159 FOR t = 0, 12-1 DO BEGIN 158 stdsc (t)= sqrt((moment(datin(long(findgen(ntime/12)*12+t))))[1])160 stdsc[t] = sqrt((moment(datin(long(findgen(ntime/12)*12+t))))[1]) 159 161 ENDFOR 160 162 ; … … 171 173 ; 172 174 FOR k = 0, ntim1 DO BEGIN 173 cx (k) = (total(datad(0:ntime-k-1)*(shift(datad, -k))(0:ntime-k-1)))/float(ntime)175 cx[k] = (total(datad[0:ntime-k-1]*(shift(datad, -k))[0:ntime-k-1]))/float(ntime) 174 176 ENDFOR 175 177 ; 176 178 ; Compute autocorrelations for plotting 177 179 ; 178 rx = cx/cx (0)180 rx = cx/cx[0] 179 181 ; 180 182 ; Compute the spectrum with no windowing … … 182 184 FOR i=0,lfreq-2 DO BEGIN 183 185 omega=twopin*float(i+1) 184 powerux (i) = (cx(0)+total((2.0*shift(cx, -1)*cos(omega*float(findgen(ntim1)+1)))(0:ntim1-1)))/pi186 powerux[i] = (cx[0]+total((2.0*shift(cx, -1)*cos(omega*float(findgen(ntim1)+1)))[0:ntim1-1]))/pi 185 187 ENDFOR 186 188 ; … … 198 200 FOR i=0,lfreq-2 DO BEGIN 199 201 omega=twopin*float(i+1) 200 powerwx (i) = (weight(0)*cx(0)+total((2.0*shift(weight, -1)*shift(cx, -1)*cos(omega*float(findgen(mm)+1)))(0:mm-1)))/pi202 powerwx[i] = (weight(0)*cx(0)+total((2.0*shift(weight, -1)*shift(cx, -1)*cos(omega*float(findgen(mm)+1)))[0:mm-1]))/pi 201 203 ENDFOR 202 204 … … 219 221 return, spectra 220 222 221 END 223 END -
trunk/procs/spectrum.pro
r165 r169 10 10 FUNCTION spectrum, ts, tinc 11 11 ; 12 compile_opt idl2, strictarrsubs 13 ; 12 14 N=n_elements(ts) 13 15 N21 = N/2 14 16 ps=ABS(FFT(ts, -1)) 15 ps=ps (0:N21-1)17 ps=ps[0:N21-1] 16 18 17 19 period=tinc*N/(findgen(N21)+0.0000000000001) 18 period (0)=99999999999.20 period[0]=99999999999. 19 21 20 22 ret=findgen(2,N21) 21 ret (1,*)=ps22 ret (0,*)=period23 ret[1,*]=ps 24 ret[0,*]=period 23 25 24 26 return, ret -
trunk/procs/stat_error.pro
r165 r169 27 27 , BAVARD=bavard $ 28 28 , MEAN=mean 29 ; 30 compile_opt idl2, strictarrsubs 29 31 ; 30 32 … … 96 98 ; ---------------- 97 99 IF keyword_set(mean) THEN BEGIN 98 errstats (it)= mean(testarr)100 errstats[it] = mean(testarr) 99 101 ENDIF ELSE BEGIN 100 ; errstats (it)= sqrt((moment(testarr))[1])101 errstats (it)= sqrt(total(testarr^2)/n_elements(testarr))102 ; errstats[it] = sqrt((moment(testarr))[1]) 103 errstats[it] = sqrt(total(testarr^2)/n_elements(testarr)) 102 104 ENDELSE 103 105 ; method 1 keep mean i.e = sqrt(sum(testarr(i)**2)/N) 104 ; method 2 errstats (it)= sqrt((moment(testarr))[1])106 ; method 2 errstats[it] = sqrt((moment(testarr))[1]) 105 107 106 108 it = it + 1 -
trunk/procs/tools_wavelets/wave_signif.pro
r165 r169 127 127 , CDELTA=cdelta $ 128 128 , PSI0=psi0 129 129 ; 130 compile_opt idl2, strictarrsubs 131 ; 130 132 ON_ERROR,2 131 133 IF (N_ELEMENTS(y) LT 1) THEN MESSAGE,'Time series Y must be input' … … 133 135 IF (N_ELEMENTS(scale) LT 1) THEN MESSAGE,'Scales must be input' 134 136 IF (N_PARAMS() LT 4) THEN sigtest = 0 ; the default 135 IF (N_ELEMENTS(y) EQ 1) THEN variance=y ELSE variance=(MOMENT(y)) (1)137 IF (N_ELEMENTS(y) EQ 1) THEN variance=y ELSE variance=(MOMENT(y))[1] 136 138 137 139 ;....check keywords & optional inputs … … 153 155 fourier_factor = (4*!PI)/(k0 + SQRT(2+k0^2)) ; [Sec.3h] 154 156 empir = [2.,-1,-1,-1] 155 IF (k0 EQ 6) THEN empir (1:3)=[0.776,2.32,0.60]157 IF (k0 EQ 6) THEN empir[1:3]=[0.776,2.32,0.60] 156 158 END 157 159 'PAUL': BEGIN ;****************** PAUL … … 159 161 fourier_factor = 4*!PI/(2*m+1) 160 162 empir = [2.,-1,-1,-1] 161 IF (m EQ 4) THEN empir (1:3)=[1.132,1.17,1.5]163 IF (m EQ 4) THEN empir[1:3]=[1.132,1.17,1.5] 162 164 END 163 165 'DOG': BEGIN ;******************* DOG … … 165 167 fourier_factor = 2*!PI*SQRT(2./(2*m+1)) 166 168 empir = [1.,-1,-1,-1] 167 IF (m EQ 2) THEN empir (1:3)= [3.541,1.43,1.4]168 IF (m EQ 6) THEN empir (1:3)= [1.966,1.37,0.97]169 IF (m EQ 2) THEN empir[1:3] = [3.541,1.43,1.4] 170 IF (m EQ 6) THEN empir[1:3] = [1.966,1.37,0.97] 169 171 END 170 172 ENDCASE 171 173 172 174 period = scale*fourier_factor 173 dofmin = empir (0); Degrees of freedom with no smoothing174 Cdelta = empir (1); reconstruction factor175 gamma = empir (2); time-decorrelation factor176 dj0 = empir (3); scale-decorrelation factor175 dofmin = empir[0] ; Degrees of freedom with no smoothing 176 Cdelta = empir[1] ; reconstruction factor 177 gamma = empir[2] ; time-decorrelation factor 178 dj0 = empir[3] ; scale-decorrelation factor 177 179 178 180 ;....significance levels [Sec.4] … … 206 208 IF (NOT confidence) THEN BEGIN 207 209 FOR a1=0,J DO BEGIN 208 chisqr = CHISQR_CVF(1. - siglvl,dof (a1))/dof(a1)209 signif (a1) = fft_theor(a1)*chisqr210 chisqr = CHISQR_CVF(1. - siglvl,dof[a1])/dof[a1] 211 signif[a1] = fft_theor[a1]*chisqr 210 212 ENDFOR 211 213 ENDIF ELSE BEGIN … … 213 215 sig = (1. - siglvl)/2. 214 216 FOR a1=0,J DO BEGIN 215 chisqr = dof (a1)/ $216 [CHISQR_CVF(sig,dof (a1)),CHISQR_CVF(1.-sig,dof(a1))]217 signif (a1,*) = fft_theor(a1)*chisqr217 chisqr = dof[a1]/ $ 218 [CHISQR_CVF(sig,dof[a1]),CHISQR_CVF(1.-sig,dof[a1])] 219 signif[a1,*] = fft_theor[a1]*chisqr 218 220 ENDFOR 219 221 ENDELSE -
trunk/procs/tools_wavelets/wavelet.pro
r165 r169 71 71 FUNCTION morlet $ 72 72 , k0, scale, k, period, coi, dofmin, cdelta, psi0 73 73 ; 74 compile_opt idl2, strictarrsubs 75 ; 74 76 IF (k0 EQ -1) THEN k0 = 6d 75 77 n = N_ELEMENTS(k) 76 78 expnt = -(scale*k - k0)^2/2d*(k GT 0.) 77 dt = 2*!PI/(n*k (1))79 dt = 2*!PI/(n*k[1]) 78 80 norm = SQRT(2*!PI*scale/dt)*(!PI^(-0.25)) ; total energy=N [Eqn(7)] 79 81 morlet = norm*EXP(expnt > (-100d)) … … 84 86 coi = fourier_factor/SQRT(2) ; Cone-of-influence [Sec.3g] 85 87 dofmin = 2 ; Degrees of freedom with no smoothing 86 Cdelta = -187 IF (k0 EQ 6) THEN Cdelta = 0.776 ; reconstruction factor88 cdelta = -1 89 IF (k0 EQ 6) THEN cdelta = 0.776 ; reconstruction factor 88 90 psi0 = !PI^(-0.25) 89 91 ; PRINT,scale,n,SQRT(TOTAL(ABS(morlet)^2,/DOUBLE)) … … 96 98 FUNCTION paul $ 97 99 , m, scale, k, period, coi, dofmin, cdelta, psi0 100 ; 101 compile_opt idl2, strictarrsubs 102 ; 98 103 99 104 IF (m EQ -1) THEN m = 4d 100 105 n = N_ELEMENTS(k) 101 106 expnt = -(scale*k)*(k GT 0.) 102 dt = 2d*!PI/(n*k (1))107 dt = 2d*!PI/(n*k[1]) 103 108 norm = SQRT(2*!PI*scale/dt)*(2^m/SQRT(m*FACTORIAL(2*m-1))) 104 109 paul = norm*((scale*k)^m)*EXP(expnt > (-100d))*(expnt GT -100) … … 108 113 coi = fourier_factor*SQRT(2) 109 114 dofmin = 2 ; Degrees of freedom with no smoothing 110 Cdelta = -1111 IF (m EQ 4) THEN Cdelta = 1.132 ; reconstruction factor115 cdelta = -1 116 IF (m EQ 4) THEN cdelta = 1.132 ; reconstruction factor 112 117 psi0 = 2.^m*FACTORIAL(m)/SQRT(!PI*FACTORIAL(2*m)) 113 118 ; PRINT,scale,n,norm,SQRT(TOTAL(paul^2,/DOUBLE))*SQRT(n) … … 120 125 FUNCTION dog $ 121 126 , m, scale, k, period, coi, dofmin, cdelta, psi0 127 ; 128 compile_opt idl2, strictarrsubs 129 ; 122 130 123 131 IF (m EQ -1) THEN m = 2 124 132 n = N_ELEMENTS(k) 125 133 expnt = -(scale*k)^2/2d 126 dt = 2d*!PI/(n*k (1))134 dt = 2d*!PI/(n*k[1]) 127 135 norm = SQRT(2*!PI*scale/dt)*SQRT(1d/GAMMA(m+0.5)) 128 136 I = DCOMPLEX(0,1) … … 132 140 coi = fourier_factor/SQRT(2) 133 141 dofmin = 1 ; Degrees of freedom with no smoothing 134 Cdelta = -1142 cdelta = -1 135 143 psi0 = -1 136 144 IF (m EQ 2) THEN BEGIN 137 Cdelta = 3.541 ; reconstruction factor145 cdelta = 3.541 ; reconstruction factor 138 146 psi0 = 0.867325 139 147 ENDIF 140 148 IF (m EQ 6) THEN BEGIN 141 Cdelta = 1.966 ; reconstruction factor149 cdelta = 1.966 ; reconstruction factor 142 150 psi0 = 0.88406 143 151 ENDIF … … 257 265 , VOICE=voice 258 266 ; 267 compile_opt idl2, strictarrsubs 268 ; 259 269 ON_ERROR,2 260 270 r = CHECK_MATH(0,1) … … 298 308 ;....construct wavenumber array used in transform [Eqn(5)] 299 309 k = (DINDGEN(n/2) + 1)*(2*!PI)/(DOUBLE(n)*dt) 300 k = [0d,k,-REVERSE(k (0:(n-1)/2 - 1))]310 k = [0d,k,-REVERSE(k[0:(n-1)/2 - 1])] 301 311 302 312 ;....compute FFT of the (padded) time series … … 318 328 FOR a1=0,na-1 DO BEGIN ;scale 319 329 psi_fft=CALL_FUNCTION(mother, $ 320 param,scale (a1),k,period1,coi,dofmin,Cdelta,psi0)330 param,scale[a1],k,period1,coi,dofmin,cdelta,psi0) 321 331 IF (do_wave) THEN $ 322 wave (*,a1)= FFT(yfft*psi_fft,1,/DOUBLE) ;wavelet transform[Eqn(4)]323 period (a1)= period1 ; save period324 fft_theor (a1)= TOTAL((ABS(psi_fft)^2)*fft_theor_k)/n332 wave[*,a1] = FFT(yfft*psi_fft,1,/DOUBLE) ;wavelet transform[Eqn(4)] 333 period[a1] = period1 ; save period 334 fft_theor[a1] = TOTAL((ABS(psi_fft)^2)*fft_theor_k)/n 325 335 IF (do_daughter) THEN $ 326 daughter (*,a1)= FFT(psi_fft,1,/DOUBLE) ; save daughter327 IF (verbose) THEN PRINT,a1,scale (a1),period(a1), $328 TOTAL(ABS(wave (*,a1))^2),CHECK_MATH(0), $336 daughter[*,a1] = FFT(psi_fft,1,/DOUBLE) ; save daughter 337 IF (verbose) THEN PRINT,a1,scale[a1],period[a1], $ 338 TOTAL(ABS(wave[*,a1])^2),CHECK_MATH(0), $ 329 339 FORMAT='(I3,3F11.3,I6)' 330 340 ENDFOR ;scale … … 333 343 334 344 IF (do_daughter) THEN $ ; shift so DAUGHTERs are in middle of array 335 daughter = [daughter (n-n1/2:*,*),daughter(0:n1/2-1,*)]345 daughter = [daughter[n-n1/2:*,*],daughter[0:n1/2-1,*]] 336 346 337 347 ;....significance levels [Sec.4] 338 sdev = (MOMENT(y1)) (1)348 sdev = (MOMENT(y1))[1] 339 349 fft_theor = sdev*fft_theor ; include time-series variance 340 350 dof = dofmin … … 342 352 343 353 IF (recon) THEN BEGIN ; Reconstruction [Eqn(11)] 344 IF ( Cdelta EQ -1) THEN BEGIN354 IF (cdelta EQ -1) THEN BEGIN 345 355 y1 = -1 346 356 MESSAGE,/INFO, $ 347 ' Cdelta undefined, cannot reconstruct with this wavelet'357 'cdelta undefined, cannot reconstruct with this wavelet' 348 358 ENDIF ELSE BEGIN 349 y1=dj*SQRT(dt)/( Cdelta*psi0)*(FLOAT(wave) # (1./SQRT(scale)))359 y1=dj*SQRT(dt)/(cdelta*psi0)*(FLOAT(wave) # (1./SQRT(scale))) 350 360 y1 = y1[0:n1-1] 351 361 ENDELSE 352 362 ENDIF 353 363 354 RETURN,wave (0:n1-1,*); get rid of padding before returning364 RETURN,wave[0:n1-1,*] ; get rid of padding before returning 355 365 356 366 END -
trunk/procs/trans_month.pro
r164 r169 6 6 ;- 7 7 FUNCTION trans_month, month 8 8 ; 9 compile_opt idl2, strictarrsubs 10 ; 9 11 CASE month OF 10 12 01: mn = 'JAN' -
trunk/procs/trends.pro
r165 r169 7 7 FUNCTION trends, fld, trend_int, type $ 8 8 , SILENT=silent 9 ; 10 compile_opt idl2, strictarrsubs 9 11 ; 10 12 @common … … 114 116 idx = where(fld EQ valmask) 115 117 fld = fldrem-fld 116 IF idx (0) NE -1 THEN fld(idx)= valmask118 IF idx[0] NE -1 THEN fld[idx] = valmask 117 119 ENDELSE 118 120 END … … 136 138 fldrem = fltarr(running) 137 139 FOR t = 0, running-1 DO BEGIN 138 fldrem (t)= mean(fld(long(findgen(lenght/running)*running+t)))140 fldrem[t] = mean(fld(long(findgen(lenght/running)*running+t))) 139 141 ENDFOR 140 142 IF fld_flag EQ 1 THEN BEGIN … … 179 181 mean_sc = fldrem 180 182 181 IF idx (0)NE -1 THEN BEGIN182 fld (idx)= valmask183 IF idx[0] NE -1 THEN BEGIN 184 fld[idx] = valmask 183 185 ENDIF 184 186 185 187 ; compute time lag correlation 186 188 IF ioverchk EQ 2 AND c_correl EQ 1 THEN BEGIN 187 IF (size(fld)) (1) EQ (size(fld_prev_t))(1)THEN BEGIN189 IF (size(fld))[1] EQ (size(fld_prev_t))[1] THEN BEGIN 188 190 lag_array = findgen(2*lag_correl+1)-lag_correl 189 191 lag_correlation = C_CORRELATE(fld, fld_prev_t, lag_array) … … 191 193 indx = where (lag_correlation EQ max(lag_correlation)) 192 194 IF indx[0] NE -1 THEN BEGIN 193 print, ' Lag of max correlation =', lag_array (indx),lag_correlation(indx)195 print, ' Lag of max correlation =', lag_array[indx],lag_correlation[indx] 194 196 print, ' Lag correlation limit [-N,N] =', lag_correl 195 197 ENDIF … … 261 263 ; keep mean seasonal cycle in common 262 264 mean_sc = fldrem 263 IF idx (0) NE -1 THEN fld(idx)= valmask265 IF idx[0] NE -1 THEN fld[idx] = valmask 264 266 ENDELSE 265 267 … … 304 306 IF int_win EQ 0 THEN BEGIN ; integral from start 305 307 FOR t = 0, lenght-1 DO BEGIN 306 fldint (t)= total(fld[0:t])308 fldint[t] = total(fld[0:t]) 307 309 ENDFOR 308 310 ENDIF ELSE BEGIN ; running integral 309 311 FOR t = 0, int_win-1 DO BEGIN 310 fldint (t)= total(fld[0:t])312 fldint[t] = total(fld[0:t]) 311 313 ENDFOR 312 314 FOR t = int_win,lenght-1 DO BEGIN 313 fldint (t)= total(fld[t-int_win+1:t])315 fldint[t] = total(fld[t-int_win+1:t]) 314 316 ENDFOR 315 317 ENDELSE … … 323 325 lenght = (size(fld))[2] 324 326 FOR i = 0, ny-1 DO BEGIN 325 zwrk = fld (i, *)327 zwrk = fld[i, *] 326 328 zint = zwrk 327 329 IF int_win EQ 0 THEN BEGIN ; integral from start 328 330 FOR t = 0, lenght-1 DO BEGIN 329 zint (t)= total(zwrk[0:t])331 zint[t] = total(zwrk[0:t]) 330 332 ENDFOR 331 333 ENDIF ELSE BEGIN ; running integral 332 334 FOR t = 0, int_win-1 DO BEGIN 333 zint (t)= total(zwrk[0:t])335 zint[t] = total(zwrk[0:t]) 334 336 ENDFOR 335 337 FOR t = int_win,lenght-1 DO BEGIN 336 zint (t)= total(zwrk[t-int_win+1:t])338 zint[t] = total(zwrk[t-int_win+1:t]) 337 339 ENDFOR 338 340 ENDELSE 339 fld (i, *)= zint341 fld[i, *] = zint 340 342 ENDFOR 341 343 END … … 344 346 lenght = (size(fld))[2] 345 347 FOR i = 0, nx-1 DO BEGIN 346 zwrk = fld (i, *)348 zwrk = fld[i, *] 347 349 zint = zwrk 348 350 IF int_win EQ 0 THEN BEGIN ; integral from start 349 351 FOR t = 0, lenght-1 DO BEGIN 350 zint (t)= total(zwrk[0:t])352 zint[t] = total(zwrk[0:t]) 351 353 ENDFOR 352 354 ENDIF ELSE BEGIN ; running integral 353 355 FOR t = 0, int_win-1 DO BEGIN 354 zint (t)= total(zwrk[0:t])356 zint[t] = total(zwrk[0:t]) 355 357 ENDFOR 356 358 FOR t = int_win,lenght-1 DO BEGIN 357 zint (t)= total(zwrk[t-int_win+1:t])359 zint[t] = total(zwrk[t-int_win+1:t]) 358 360 ENDFOR 359 361 ENDELSE 360 fld (i, *)= zint362 fld[i, *] = zint 361 363 ENDFOR 362 364 END -
trunk/procs/tvnplot.pro
r165 r169 18 18 , PS=ps $ 19 19 , _EXTRA=extra 20 20 ; 21 compile_opt idl2, strictarrsubs 22 ; 21 23 IF keyword_set(min) EQ 0 THEN min = min(image) 22 24 IF keyword_set(max) EQ 0 THEN max = max(image) … … 55 57 image2 = ((min > image < max) -min)*(!d.n_colors-2)/(max-min)+1 56 58 IF n_elements(msk) NE 0 THEN BEGIN 57 image2 (where(msk EQ 0))= 059 image2[where(msk EQ 0)] = 0 58 60 ENDIF 59 61 image2 = rebin(image2, gamma*ni, gamma*nj, sample = sample) … … 98 100 xticklen = 1, yticklen = 1, color = 0, _EXTRA=extra ; cadre 99 101 100 x0 = !x.window (0)101 y0 = !y.window (0)102 dx = !x.window (1)-!x.window(0)103 dy = !y.window (1)-!y.window(0)102 x0 = !x.window[0] 103 y0 = !y.window[0] 104 dx = !x.window[1]-!x.window[0] 105 dy = !y.window[1]-!y.window[0] 104 106 105 107 deltax = (!x.window[1]-!x.window[0]) … … 122 124 CASE button OF 123 125 1: BEGIN 124 print, 'i=', i+ibase, ' / j=', j+jbase, ' / z (i,j)=', image(i, j)126 print, 'i=', i+ibase, ' / j=', j+jbase, ' / z[i,j]=', image[i, j] 125 127 wait, .2 126 128 END … … 133 135 IF i1 LE i0 THEN i1 = (i0+10) < (ni-1) 134 136 IF j1 LE j0 THEN j1 = (j0+10) < (nj-1) 135 tvnplot, image (i0:i1, j0:j1), _EXTRA=extra, /erase, $137 tvnplot, image[i0:i1, j0:j1], _EXTRA=extra, /erase, $ 136 138 ibase = i0+ibase, jbase = j0+jbase, MIN = min, Max = max 137 139 END -
trunk/procs/update_data.pro
r165 r169 17 17 , NO_MEAN=no_mean $ 18 18 , _EXTRA=extra 19 ; 20 compile_opt idl2, strictarrsubs 19 21 ; 20 22 @common … … 44 46 'z': BEGIN 45 47 IF vert_mean[0] EQ vert_mean[1] THEN BEGIN 46 suffix = ' at '+strtrim(string(long(vert_mean (0))), 2)+zunits+' '47 ENDIF ELSE suffix = ' at '+strtrim(string(long(vert_mean (0))), 2)+'-'+strtrim(string(long(vert_mean(1))), 2)+zunits+' '48 suffix = ' at '+strtrim(string(long(vert_mean[0])), 2)+zunits+' ' 49 ENDIF ELSE suffix = ' at '+strtrim(string(long(vert_mean[0])), 2)+'-'+strtrim(string(long(vert_mean[1])), 2)+zunits+' ' 48 50 END 49 51 'level':BEGIN ; case level (zindex) 50 suffix = ' at '+vert_type+' '+strtrim(string(long(gdept (vert_mean(0))+.1)), 2)+' '52 suffix = ' at '+vert_type+' '+strtrim(string(long(gdept[vert_mean[0]]+.1)), 2)+' ' 51 53 END 52 54 ENDCASE … … 66 68 'z': BEGIN 67 69 zmean = moyenne(temporary(tab), 'z', boite = boxzoom, /NAN) 68 suffix = ' averaged in ['+strtrim(string(long(vert_mean (0))), 2)+zunits+', '+strtrim(string(long(vert_mean(1))), 2)+zunits+']'70 suffix = ' averaged in ['+strtrim(string(long(vert_mean[0])), 2)+zunits+', '+strtrim(string(long(vert_mean[1])), 2)+zunits+']' 69 71 END 70 72 'level':BEGIN ; case level (zindex) 71 73 zmean = moyenne(temporary(tab), 'z', boite = boxzoom, /zindex, /NAN) 72 suffix = ' averaged in '+vert_type+' ['+strtrim(string(long(vert_mean (0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']'74 suffix = ' averaged in '+vert_type+' ['+strtrim(string(long(vert_mean[0])), 2)+','+strtrim(string(long(vert_mean[1])), 2)+']' 73 75 END 74 76 ENDCASE … … 93 95 'z': BEGIN 94 96 zmean = grossemoyenne(tab, 'z', boite = boxzoom, /NAN) 95 suffix = ' averaged in ['+strtrim(string(long(vert_mean (0))), 2)+zunits+','+strtrim(string(long(vert_mean(1))), 2)+zunits+']'97 suffix = ' averaged in ['+strtrim(string(long(vert_mean[0])), 2)+zunits+','+strtrim(string(long(vert_mean[1])), 2)+zunits+']' 96 98 END 97 99 'level': BEGIN ; case level (zindex) 98 100 zmean = grossemoyenne(temporary(tab), 'z', boite = boxzoom, /zindex, /NAN) 99 suffix = ' averaged in '+vert_type+' ['+strtrim(string(long(vert_mean (0))), 2)+','+strtrim(string(long(vert_mean(1))), 2)+']'101 suffix = ' averaged in '+vert_type+' ['+strtrim(string(long(vert_mean[0])), 2)+','+strtrim(string(long(vert_mean[1])), 2)+']' 100 102 END 101 103 ENDCASE -
trunk/procs/ybinx.pro
r165 r169 16 16 ; fld2 is the field being bined 17 17 18 IF (size(bin_interval)) (0)EQ 1 THEN BEGIN18 IF (size(bin_interval))[0] EQ 1 THEN BEGIN 19 19 ; build... 20 20 ENDIF 21 21 22 nbins = (size(bin_interval)) (1)22 nbins = (size(bin_interval))[1] 23 23 24 24 bin_plt = (bin_interval+shift(bin_interval, -1))/2 25 bin_plt = [2*bin_plt (0)-bin_plt(1), bin_plt[0:nbins-2], 2*bin_plt(nbins-2)-bin_plt(nbins-3)]25 bin_plt = [2*bin_plt[0]-bin_plt[1], bin_plt[0:nbins-2], 2*bin_plt[nbins-2]-bin_plt[nbins-3]] 26 26 27 27 IF debug_w THEN print, ' bin_plt= ', bin_plt … … 36 36 IF sw3 THEN BEGIN 37 37 idxmsk = where(fld GT valmask/10.) 38 IF idxmsk (0) NE -1 THEN fld(idxmsk)= valmask39 IF idxmsk (0) NE -1 THEN fld3(idxmsk)= valmask38 IF idxmsk[0] NE -1 THEN fld[idxmsk] = valmask 39 IF idxmsk[0] NE -1 THEN fld3[idxmsk] = valmask 40 40 ENDIF 41 41 … … 44 44 idxmskpn = where(finite(fld, /nan)) 45 45 idxmskpn2 = where(finite(fld2, /nan)) 46 IF idxmskpn (0) NE -1 THEN fld(idxmskpn)= valmask47 IF idxmskpn2 (0) NE -1 THEN fld2(idxmskpn2)= valmask46 IF idxmskpn[0] NE -1 THEN fld[idxmskpn] = valmask 47 IF idxmskpn2[0] NE -1 THEN fld2[idxmskpn2] = valmask 48 48 49 49 idxmskpn = where(fld GT valmask/10.) 50 50 idxmskpn2 = where(fld2 GT valmask/10.) 51 IF idxmskpn (0) NE -1 THEN fld2(idxmskpn)= valmask52 IF idxmskpn2 (0) NE -1 THEN fld(idxmskpn2)= valmask51 IF idxmskpn[0] NE -1 THEN fld2[idxmskpn] = valmask 52 IF idxmskpn2[0] NE -1 THEN fld[idxmskpn2] = valmask 53 53 54 54 ; print min/max of field for debug 55 55 idxmskp = where(fld LE valmask/10.) 56 IF debug_w THEN print, ' Min/max fld= ', min(fld (idxmskp)), max(fld(idxmskp))57 IF debug_w THEN print, ' Min/max fld2= ', min(fld2 (idxmskp)), max(fld2(idxmskp))56 IF debug_w THEN print, ' Min/max fld= ', min(fld[idxmskp]), max(fld[idxmskp]) 57 IF debug_w THEN print, ' Min/max fld2= ', min(fld2[idxmskp]), max(fld2[idxmskp]) 58 58 59 59 ; remove mean seasonal cycle if required … … 77 77 @mth_decode 78 78 79 fldn = (fld)[*, *, reform(idxm (0,*), njpt)]80 fld2n = (fld2)[*, *, reform(idxm (0,*), njpt)]81 IF sw3 THEN fld3n = (fld3)[*, *, reform(idxm (0,*), njpt)]79 fldn = (fld)[*, *, reform(idxm[0,*], njpt)] 80 fld2n = (fld2)[*, *, reform(idxm[0,*], njpt)] 81 IF sw3 THEN fld3n = (fld3)[*, *, reform(idxm[0,*], njpt)] 82 82 83 83 FOR imth = 1, nmth-1 DO BEGIN 84 84 85 fldn = [fldn, (fld)[*, *, reform(idxm (imth,*), njpt)]]86 fld2n = [fld2n, (fld2)[*, *, reform(idxm (imth,*), njpt)]]87 IF sw3 THEN fld3n = [fld3n, (fld3)[*, *, reform(idxm (imth,*), njpt)]]85 fldn = [fldn, (fld)[*, *, reform(idxm[imth,*], njpt)]] 86 fld2n = [fld2n, (fld2)[*, *, reform(idxm[imth,*], njpt)]] 87 IF sw3 THEN fld3n = [fld3n, (fld3)[*, *, reform(idxm[imth,*], njpt)]] 88 88 89 89 ENDFOR … … 110 110 WHILE ib LE nbins-1 DO BEGIN 111 111 112 indices = where(fld2 GT bin_interval (ib-1) AND fld2 LE bin_interval(ib))113 binpop (ib)= n_elements(indices)114 idxb[ib, 0:binpop (ib)-1] = where(fld2 GT bin_interval(ib-1) AND fld2 LE bin_interval(ib))115 IF debug_w THEN print, ' Size of bin (i) ', ib, binpop(ib)112 indices = where(fld2 GT bin_interval[ib-1] AND fld2 LE bin_interval[ib]) 113 binpop[ib] = n_elements(indices) 114 idxb[ib, 0:binpop[ib]-1] = where(fld2 GT bin_interval[ib-1] AND fld2 LE bin_interval[ib]) 115 IF debug_w THEN print, ' Size of bin[i] ', ib, binpop[ib] 116 116 117 117 ib = ib + 1 118 118 ENDWHILE 119 119 120 index0 = where(fld2 LE bin_interval (0))121 binpop (0)= n_elements(index0)122 idxb (0, 0:binpop(0)-1)= index0123 124 indexm = n_elements(where(fld2 GT bin_interval (nbins-1)))125 binpop (nbins)= n_elements(indexm)126 idxb (nbins, 0:binpop(nbins)-1)= indexm120 index0 = where(fld2 LE bin_interval[0]) 121 binpop[0] = n_elements(index0) 122 idxb[0, 0:binpop[0]-1] = index0 123 124 indexm = n_elements(where(fld2 GT bin_interval[nbins-1])) 125 binpop[nbins] = n_elements(indexm) 126 idxb[nbins, 0:binpop[nbins]-1] = indexm 127 127 128 128 ; compute maximum number of indices for one bin … … 148 148 ; first build fld*e1te2t for later average 149 149 150 surf = reform(e1t (firstxt:lastxt,firstyt:lastyt)*e2t(firstxt:lastxt,firstyt:lastyt), nxt*nyt)150 surf = reform(e1t[firstxt:lastxt,firstyt:lastyt]*e2t[firstxt:lastxt,firstyt:lastyt], nxt*nyt) 151 151 surf = reform(surf#replicate(1, jpt), nxt, nyt, jpt) 152 152 … … 157 157 158 158 IF debug_w THEN print, 'bin = ', ib 159 binsz = binpop (ib)159 binsz = binpop[ib] 160 160 IF binsz GT 1 THEN BEGIN 161 fldy (ib, 0:binsz-1) = fld(idxb(ib, 0:binsz-1))162 IF debug_w THEN print, 'fld (idxb(ib, 0:binsz-1)) =',fld(idxb(ib, 0:binsz-1))163 fldys (ib, 0:binsz-1) = flds(idxb(ib, 0:binsz-1))164 surfb (ib, 0:binsz-1) = surf(idxb(ib, 0:binsz-1))165 IF sw3 THEN fldy2 (ib, 0:binsz-1) = fld3(idxb(ib, 0:binsz-1))161 fldy[ib, 0:binsz-1] = fld[idxb[ib, 0:binsz-1]] 162 IF debug_w THEN print, 'fld[idxb[ib, 0:binsz-1]] =',fld[idxb[ib, 0:binsz-1]] 163 fldys[ib, 0:binsz-1] = flds[idxb[ib, 0:binsz-1]] 164 surfb[ib, 0:binsz-1] = surf[idxb[ib, 0:binsz-1]] 165 IF sw3 THEN fldy2[ib, 0:binsz-1] = fld3[idxb[ib, 0:binsz-1]] 166 166 ENDIF 167 167 ib = ib + 1 … … 169 169 170 170 171 yplt = fltarr (nbins+1)172 yerr = fltarr (nbins+1)171 yplt = fltarr[nbins+1] 172 yerr = fltarr[nbins+1] 173 173 174 174 IF NOT sw3 THEN BEGIN … … 182 182 WHILE ib LE nbins DO BEGIN 183 183 184 binsz = binpop (ib)184 binsz = binpop[ib] 185 185 186 186 IF binsz GT 1 THEN BEGIN 187 187 188 sfc_tot = total(surfb (ib, 0:binsz-1))189 yplt (ib) = total(fldys(ib, 0:binsz-1))/sfc_tot190 yerr (ib) = sqrt((moment(fldy(ib, 0:binsz-1)))[1])191 192 mean_fld = mean_fld + yplt (ib)*float(binpop(ib))188 sfc_tot = total(surfb[ib, 0:binsz-1]) 189 yplt[ib] = total(fldys[ib, 0:binsz-1])/sfc_tot 190 yerr[ib] = sqrt((moment(fldy[ib, 0:binsz-1]))[1]) 191 192 mean_fld = mean_fld + yplt[ib]*float(binpop[ib]) 193 193 194 194 ; print bin info 195 IF ib GT 0 AND ib LT nbins THEN print, ' Bin size, occurence, average: ', bin_interval (ib-1), bin_interval(ib), binpop(ib), (binpop(ib)/total(binpop))*100., yplt(ib)196 IF ib EQ 0 THEN print, ' Bin size, occurence, average: min' , bin_interval (ib), binpop(ib), (binpop(ib)/total(binpop))*100. , yplt(ib)197 IF ib EQ nbins THEN print, ' Bin size, occurence, average: ', bin_interval (ib-1),' max ', binpop(ib), (binpop(ib)/total(binpop))*100. , yplt(ib)198 ENDIF ELSE yplt (ib)= !values.f_nan195 IF ib GT 0 AND ib LT nbins THEN print, ' Bin size, occurence, average: ', bin_interval[ib-1], bin_interval[ib], binpop[ib], (binpop[ib]/total(binpop))*100., yplt[ib] 196 IF ib EQ 0 THEN print, ' Bin size, occurence, average: min' , bin_interval[ib], binpop[ib], (binpop[ib]/total(binpop))*100. , yplt[ib] 197 IF ib EQ nbins THEN print, ' Bin size, occurence, average: ', bin_interval[ib-1],' max ', binpop[ib], (binpop[ib]/total(binpop))*100. , yplt[ib] 198 ENDIF ELSE yplt[ib] = !values.f_nan 199 199 200 200 ib = ib + 1 … … 210 210 211 211 ENDIF ELSE BEGIN 212 ; if bining fld1=f (fld3), compute regression in each bin and plot212 ; if bining fld1=f[fld3], compute regression in each bin and plot 213 213 214 214 ; compute regression for each bin … … 219 219 WHILE ib LE nbins DO BEGIN 220 220 221 binsz = binpop (ib)221 binsz = binpop[ib] 222 222 223 223 IF binsz GT 1 THEN BEGIN 224 224 225 idx1 = where(fldy (ib, 0:binsz-1)NE valmask)226 idx2 = where(fldy2 (ib, 0:binsz-1)NE valmask)227 tab1 = (fldy (ib, 0:binsz-1))(idx1)228 tab2 = (fldy2 (ib, 0:binsz-1))(idx2)225 idx1 = where(fldy[ib, 0:binsz-1] NE valmask) 226 idx2 = where(fldy2[ib, 0:binsz-1] NE valmask) 227 tab1 = (fldy[ib, 0:binsz-1])[idx1] 228 tab2 = (fldy2[ib, 0:binsz-1])[idx2] 229 229 coeff = linfit(tab2, tab1, CHISQ = linerr, PROB = proberr, SIGMA = sigmaerr) 230 yplt (ib) = coeff(1)231 yerr (ib) = sigmaerr(1)232 233 mean_fld = mean_fld + yplt (ib)*float(binpop(ib))230 yplt[ib] = coeff[1] 231 yerr[ib] = sigmaerr[1] 232 233 mean_fld = mean_fld + yplt[ib]*float(binpop[ib]) 234 234 235 235 ; print bin info 236 IF ib GT 0 AND ib LT nbins THEN print, ' Bin size, occurence, regress.: ', bin_interval (ib-1), bin_interval(ib), binpop(ib), (binpop(ib)/total(binpop))*100., yplt(ib)237 IF ib EQ 0 THEN print, ' Bin size, occurence, regress.: min' , bin_interval (ib), binpop(ib), (binpop(ib)/total(binpop))*100. , yplt(ib)238 IF ib EQ nbins THEN print, ' Bin size, occurence, regress.: ', bin_interval (ib-1),' max ', binpop(ib), (binpop(ib)/total(binpop))*100. , yplt(ib)239 240 ENDIF ELSE yplt (ib)= !values.f_nan236 IF ib GT 0 AND ib LT nbins THEN print, ' Bin size, occurence, regress.: ', bin_interval[ib-1], bin_interval[ib], binpop[ib], (binpop[ib]/total(binpop))*100., yplt[ib] 237 IF ib EQ 0 THEN print, ' Bin size, occurence, regress.: min' , bin_interval[ib], binpop[ib], (binpop[ib]/total(binpop))*100. , yplt[ib] 238 IF ib EQ nbins THEN print, ' Bin size, occurence, regress.: ', bin_interval[ib-1],' max ', binpop[ib], (binpop[ib]/total(binpop))*100. , yplt[ib] 239 240 ENDIF ELSE yplt[ib] = !values.f_nan 241 241 242 242 ib = ib + 1 … … 261 261 varunit = ' ['+ntxt+' in box '+leg_name+']' 262 262 263 overc = overlay_type (iover, dimplot)263 overc = overlay_type[iover, dimplot] 264 264 265 265 pltcmd = 'pltsc,yplt,bin_plt,minc,maxc,minc2,maxc2,varlegend'+com_strplt+overc+',STY1D=-6,subtitle=""' -
trunk/procs/yfx.pro
r165 r169 95 95 ; whole serie 96 96 coeff = linfit(fld2, fld, CHISQ = linerr, PROB = proberr, SIGMA = sigmaerr) 97 print, ' Linfit coef (a+bx) = ', coeff (0), coeff(1), ' errbar = ',sigmaerr97 print, ' Linfit coef (a+bx) = ', coeff[0], coeff[1], ' errbar = ',sigmaerr 98 98 ; text for plot info in legend_overlay 99 coeff_txt = 'Linfit = '+strtrim(string(coeff (1), format = fmt_linfit), 2)99 coeff_txt = 'Linfit = '+strtrim(string(coeff[1], format = fmt_linfit), 2) 100 100 IF proberr NE 1.0 THEN print, ' WARNING: proba = ',proberr 101 101 ; divide serie in two domains (on x axis) seperated by linfit_sep (plt_def) … … 108 108 coefflu_txt = ' [' 109 109 110 IF idx_low (0)NE -1 THEN BEGIN110 IF idx_low[0] NE -1 THEN BEGIN 111 111 coeffl = linfit(fld2(idx_low), fld(idx_low), CHISQ = linerrl, PROB = proberrl, SIGMA = sigmaerrl) 112 print, ' Linfit coef below ',strtrim(string(linfit_sep), 2),' = ', coeffl (0), coeffl(1), ' errbar = ',sigmaerrl112 print, ' Linfit coef below ',strtrim(string(linfit_sep), 2),' = ', coeffl[0], coeffl[1], ' errbar = ',sigmaerrl 113 113 IF proberrl NE 1.0 THEN print, ' *** WARNING: proba lower = ',proberrl 114 coefflu_txt = coefflu_txt+strtrim(string(coeffl (1), format = fmt_linfit), 2)+'<'+strtrim(string(linfit_sep, format = '(F4.2)'), 2)114 coefflu_txt = coefflu_txt+strtrim(string(coeffl[1], format = fmt_linfit), 2)+'<'+strtrim(string(linfit_sep, format = '(F4.2)'), 2) 115 115 ENDIF 116 116 117 IF idx_up (0)NE -1 THEN BEGIN117 IF idx_up[0] NE -1 THEN BEGIN 118 118 coeffu = linfit(fld2(idx_up), fld(idx_up), CHISQ = linerru, PROB = proberru, SIGMA = sigmaerru) 119 print, ' Linfit coef above ',strtrim(string(linfit_sep), 2),' = ', coeffu (0), coeffu(1), ' errbar = ',sigmaerru119 print, ' Linfit coef above ',strtrim(string(linfit_sep), 2),' = ', coeffu[0], coeffu[1], ' errbar = ',sigmaerru 120 120 IF proberru NE 1.0 THEN print, ' *** WARNING: proba upper = ',proberru 121 coefflu_txt = coefflu_txt+'>'+strtrim(string(coeffu (1), format = fmt_linfit), 2)121 coefflu_txt = coefflu_txt+'>'+strtrim(string(coeffu[1], format = fmt_linfit), 2) 122 122 ENDIF 123 123 coefflu_txt = coefflu_txt+']' … … 132 132 cm3 = linfit(shift(fld2, -3), fld, CHISQ = linerrm3, PROB = proberrm3, SIGMA = sigmaerrm3) 133 133 print, ' ' 134 print, ' -3 Linfit coef (a+bx) = ', cm3 (0), cm3(1), ' errbar = ',sigmaerrm3135 print, ' -2 Linfit coef (a+bx) = ', cm2 (0), cm2(1), ' errbar = ',sigmaerrm2136 print, ' -1 Linfit coef (a+bx) = ', cm1 (0), cm1(1), ' errbar = ',sigmaerrm1137 print, ' +1 Linfit coef (a+bx) = ', cp1 (0), cp1(1), ' errbar = ',sigmaerrp1138 print, ' +2 Linfit coef (a+bx) = ', cp2 (0), cp2(1), ' errbar = ',sigmaerrp2139 print, ' +3 Linfit coef (a+bx) = ', cp3 (0), cp3(1), ' errbar = ',sigmaerrp3134 print, ' -3 Linfit coef (a+bx) = ', cm3[0], cm3[1], ' errbar = ',sigmaerrm3 135 print, ' -2 Linfit coef (a+bx) = ', cm2[0], cm2[1], ' errbar = ',sigmaerrm2 136 print, ' -1 Linfit coef (a+bx) = ', cm1[0], cm1[1], ' errbar = ',sigmaerrm1 137 print, ' +1 Linfit coef (a+bx) = ', cp1[0], cp1[1], ' errbar = ',sigmaerrp1 138 print, ' +2 Linfit coef (a+bx) = ', cp2[0], cp2[1], ' errbar = ',sigmaerrp2 139 print, ' +3 Linfit coef (a+bx) = ', cp3[0], cp3[1], ' errbar = ',sigmaerrp3 140 140 141 141 ; plot domain … … 152 152 IF strpos(symbol_families, 'x') NE -1 THEN BEGIN 153 153 coding = strsplit(symbol_families, 'x', /EXTRACT) 154 ncolors = long(coding (0))155 ntimes = long(coding (1))154 ncolors = long(coding[0]) 155 ntimes = long(coding[1]) 156 156 ENDIF ELSE BEGIN 157 157 ncolors = long(symbol_families) … … 175 175 jl = jl+1 176 176 ENDWHILE 177 fldp = fld (idx)178 fldp2 = fld2 (idx)177 fldp = fld[idx] 178 fldp2 = fld2[idx] 179 179 ; coeff = linfit(fldp2(where (fldp2(*) GT 0.)), fldp(where (fldp2(*) GT 0.)), CHISQ = errlin, PROB = prblin) 180 180 coeff = linfit(fldp2, fldp, CHISQ = errlin, PROB = prblin, SIGMA = sigmaerr) 181 linerro (ic) = sigmaerr(1)182 linprob (ic)= prblin183 lincoef (ic) = coeff(1)184 linerr0 (ic) = sigmaerr(0)185 print, ' Period '+strtrim(string(ic+1), 2)+' Linfit coef =', coeff (0),coeff(1), ' errbar = ',sigmaerr186 coeffm_txt = coeffm_txt+strtrim(string(coeff (1), format = fmt_linfit), 2)+'/'181 linerro[ic] = sigmaerr[1] 182 linprob[ic] = prblin 183 lincoef[ic] = coeff[1] 184 linerr0[ic] = sigmaerr[0] 185 print, ' Period '+strtrim(string(ic+1), 2)+' Linfit coef =', coeff[0],coeff[1], ' errbar = ',sigmaerr 186 coeffm_txt = coeffm_txt+strtrim(string(coeff[1], format = fmt_linfit), 2)+'/' 187 187 valmask = ABS(valmask) 188 188 IF mean_sc_only EQ 0 OR mean_sc_only EQ 1 THEN BEGIN 189 pltcmd = 'pltsc,fldp,fldp2,minc,maxc,minc2,maxc2,varlegend2, boite=boxyfx'+com_strplt+',ov1d='+string (ic)+',COLOR='+string(symbol_color(ic))+', STY1D='+string(symbol_style(ic)-1)+', SYMSIZE='+string(symbol_size)189 pltcmd = 'pltsc,fldp,fldp2,minc,maxc,minc2,maxc2,varlegend2, boite=boxyfx'+com_strplt+',ov1d='+string[ic]+',COLOR='+string(symbol_color[ic])+', STY1D='+string(symbol_style[ic]-1)+', SYMSIZE='+string(symbol_size) 190 190 printf, nulhis, strcompress(pltcmd) 191 191 IF debug_w THEN print, ' ',pltcmd … … 193 193 ENDIF 194 194 IF mean_sc_only EQ 4 AND ic EQ ncolors-1 THEN BEGIN 195 print, ' Linfit slope + error (July-Dec)=', mean(lincoef (5:10)), mean(linerro(5:10))196 print, ' Linfit slope + error (Jan-June)=', mean([lincoef (0:4), lincoef(11)]), mean([linerro(0:4), linerro(11)])195 print, ' Linfit slope + error (July-Dec)=', mean(lincoef[5:10]), mean(linerro[5:10]) 196 print, ' Linfit slope + error (Jan-June)=', mean([lincoef[0:4], lincoef[11]]), mean([linerro[0:4], linerro[11]]) 197 197 vardate = 'toto' 198 198 jpt = ncolors … … 226 226 ENDIF ELSE BEGIN ; one color 227 227 IF debug_w THEN print, ' ',pltcmdstd 228 res = execute(pltcmdstd (0))228 res = execute(pltcmdstd[0]) 229 229 ENDELSE 230 230 ENDIF 231 231 ENDIF ELSE BEGIN ; standard scatter plot 232 232 IF debug_w THEN print, ' ',pltcmdstd 233 res = execute(pltcmdstd (0))233 res = execute(pltcmdstd[0]) 234 234 ENDELSE 235 235 IF cmd.trend EQ 412 THEN coeffm_txt = '' … … 250 250 fldrem_t2 = fltarr(running) 251 251 FOR t1 = 0, running-1 DO BEGIN 252 fldrem_t1 (t1)= mean(fld(long(findgen(lenght/running)*running+t1)))253 fldrem_t2 (t1)= mean(fld2(long(findgen(lenght/running)*running+t1)))252 fldrem_t1[t1] = mean(fld(long(findgen(lenght/running)*running+t1))) 253 fldrem_t2[t1] = mean(fld2(long(findgen(lenght/running)*running+t1))) 254 254 ENDFOR 255 255 ENDELSE … … 258 258 print, ' Seasonal cycle var/stddev fld2 ', (moment(fldrem_t2))[1], sqrt((moment(fldrem_t2))[1]) 259 259 260 fldrem_t1 = [fldrem_t1, fldrem_t1 (0)]261 fldrem_t2 = [fldrem_t2, fldrem_t2 (0)]260 fldrem_t1 = [fldrem_t1, fldrem_t1[0]] 261 fldrem_t2 = [fldrem_t2, fldrem_t2[0]] 262 262 sw_ov1d = mean_sc_only 263 263 IF mean_sc_only EQ 3 AND cmd.trend EQ '412' THEN BEGIN … … 267 267 print, ' Stddev of monthly stddevs fld', sqrt((moment(stat.sc_std-(moment(stat.sc_std))[0]))[1]) 268 268 print, ' Stddev of monthly stddevs fld2', sqrt((moment(stat2.sc_std-(moment(stat2.sc_std))[0]))[1]) 269 stdsc = [stat.sc_std, stat.sc_std (0)]270 stdsc2 = [stat2.sc_std, stat2.sc_std (0)]269 stdsc = [stat.sc_std, stat.sc_std[0]] 270 stdsc2 = [stat2.sc_std, stat2.sc_std[0]] 271 271 pltcmd = 'pltsc,stdsc,stdsc2,0. ,maxc,0.,maxc2,varlegend2, boite=boxyfx'+com_strplt+',ov1d=1-min(1,sw_ov1d), STY1D=-1, THICKN=3, SYMSIZE='+string(symbol_size)+', FRACTION=fraction' 272 272 printf, nulhis, strcompress(pltcmd) 273 273 IF debug_w THEN print, pltcmd 274 res = execute(pltcmd (0))275 pltcmd = 'xyouts,stdsc2-(maxc2)/(.3*win (0)*win(1)),stdsc,string(long(findgen(12))+1),charsize=1.3,alignment=0.5'274 res = execute(pltcmd[0]) 275 pltcmd = 'xyouts,stdsc2-(maxc2)/(.3*win[0]*win[1]),stdsc,string(long(findgen(12))+1),charsize=1.3,alignment=0.5' 276 276 printf, nulhis, strcompress(pltcmd) 277 277 IF debug_w THEN print, pltcmd 278 res = execute(pltcmd (0))278 res = execute(pltcmd[0]) 279 279 ENDIF ELSE BEGIN 280 280 ov1d_val = 0 … … 283 283 printf, nulhis, strcompress(pltcmd) 284 284 IF debug_w THEN print, pltcmd 285 IF mean_sc_only NE 0 THEN res = execute(pltcmd (0))285 IF mean_sc_only NE 0 THEN res = execute(pltcmd[0]) 286 286 IF mean_sc_only EQ 2 THEN BEGIN 287 ; add month info win (0)*win(1)288 pltcmd = 'xyouts,fldrem_t2-(maxc2-minc2)/(0.5*win (0)*win(1)),fldrem_t1,string(long(findgen(12))+1),charsize=1.3,alignment=0.5'287 ; add month info win[0]*win[1] 288 pltcmd = 'xyouts,fldrem_t2-(maxc2-minc2)/(0.5*win[0]*win[1]),fldrem_t1,string(long(findgen(12))+1),charsize=1.3,alignment=0.5' 289 289 printf, nulhis, strcompress(pltcmd) 290 290 IF debug_w THEN print, pltcmd 291 res = execute(pltcmd (0))291 res = execute(pltcmd[0]) 292 292 ENDIF 293 293 ENDELSE … … 322 322 levels = findgen((max_sig-min_sig)/delta_sig+1)*delta_sig+min_sig 323 323 labels = levels 324 labels (*)= 1324 labels[*] = 1 325 325 contour, rhopn, s, t, /overplot, levels = levels, c_labels = labels, c_linestyle = 0 326 326 t = t_ -
trunk/tools/bsf/boundperio.pro
r165 r169 11 11 , VV=vv $ 12 12 , FF1=ff1 13 ; 14 compile_opt idl2, strictarrsubs 13 15 ; 14 16 @common … … 25 27 CASE (size(z))[0] OF 26 28 2 : BEGIN 27 zz (0, *) = zz(jpi-2, *)28 zz (jpi-1, *) = zz(1, *)29 zz[0, *] = zz[jpi-2, *] 30 zz[jpi-1, *] = zz[1, *] 29 31 CASE type OF 30 32 0:BEGIN 31 zz (*, jpj-1) = shift(reverse( zz(*, jpj-3)), 1)32 zz (milieu:jpi-1, jpj-2) = reverse( zz(1:milieu, jpj-2))33 zz[*, jpj-1] = shift(reverse( zz[*, jpj-3]), 1) 34 zz[milieu:jpi-1, jpj-2] = reverse( zz[1:milieu, jpj-2] ) 33 35 END 34 36 ; Point u 35 37 1:BEGIN 36 zz (*, jpj-1) = -reverse( zz(*, jpj-3))37 zz (milieu-1:jpi-1, jpj-2) = -reverse( zz(0:milieu, jpj-2))38 zz[*, jpj-1] = -reverse( zz[*, jpj-3]) 39 zz[milieu-1:jpi-1, jpj-2] = -reverse( zz[0:milieu, jpj-2] ) 38 40 END 39 41 ; Point v 40 42 2:BEGIN 41 zz (*, jpj-2) = -shift(reverse( zz(*, jpj-3)), 1)42 zz (*, jpj-1) = -shift(reverse( zz(*, jpj-4)), 1)43 zz[*, jpj-2] = -shift(reverse( zz[*, jpj-3]), 1) 44 zz[*, jpj-1] = -shift(reverse( zz[*, jpj-4]), 1) 43 45 END 44 46 ; Point f 45 47 3:BEGIN 46 zz (*, jpj-2) = reverse( zz(*, jpj-3))47 zz (*, jpj-1) = reverse( zz(*, jpj-4))48 zz[*, jpj-2] = reverse( zz[*, jpj-3]) 49 zz[*, jpj-1] = reverse( zz[*, jpj-4]) 48 50 END 49 51 ENDCASE 50 52 END 51 53 3 : BEGIN 52 zz (0, *, *) = zz(jpi-2, *, *)53 zz (jpi-1, *, *) = zz(1, *, *)54 zz[0, *, *] = zz[jpi-2, *, *] 55 zz[jpi-1, *, *] = zz[1, *, *] 54 56 CASE type OF 55 57 0:BEGIN 56 zz (*, jpj-1, *) = shift(reverse( zz(*, jpj-3, *)), 1, 0, 0)57 zz (milieu:jpi-1, jpj-2, *) = reverse( zz(1:milieu, jpj-2, *))58 zz[*, jpj-1, *] = shift(reverse( zz[*, jpj-3, *]), 1, 0, 0) 59 zz[milieu:jpi-1, jpj-2, *] = reverse( zz[1:milieu, jpj-2, *] ) 58 60 END 59 61 ; Point u 60 62 1:BEGIN 61 zz (*, jpj-1, *) = reverse( zz(*, jpj-3, *))62 zz (milieu-1:jpi-1, jpj-2, *) = -reverse( zz(0:milieu, jpj-2, *))63 zz[*, jpj-1, *] = reverse( zz[*, jpj-3, *]) 64 zz[milieu-1:jpi-1, jpj-2, *] = -reverse( zz[0:milieu, jpj-2, *] ) 63 65 END 64 66 ; Point v 65 67 2:BEGIN 66 zz (*, jpj-2, *) = -shift(reverse( zz(*, jpj-3, *)), 1, 0, 0)67 zz (*, jpj-1, *) = -shift(reverse( zz(*, jpj-4, *)), 1, 0, 0)68 zz[*, jpj-2, *] = -shift(reverse( zz[*, jpj-3, *]), 1, 0, 0) 69 zz[*, jpj-1, *] = -shift(reverse( zz[*, jpj-4, *]), 1, 0, 0) 68 70 END 69 71 ; Point f 70 72 3:BEGIN 71 zz (*, jpj-2, *) = reverse( zz(*, jpj-3, *))72 zz (*, jpj-1, *) = reverse( zz(*, jpj-4, *))73 zz[*, jpj-2, *] = reverse( zz[*, jpj-3, *]) 74 zz[*, jpj-1, *] = reverse( zz[*, jpj-4, *]) 73 75 END 74 76 ENDCASE -
trunk/tools/bsf/build_laplacien_f.pro
r165 r169 23 23 PRO build_laplacien_f 24 24 ; 25 compile_opt idl2, strictarrsubs 26 ; 25 27 @common 26 28 @com_eg … … 32 34 33 35 34 fm = fmaskr (*, *, 0)/(e1f*e2f)36 fm = fmaskr[*, *, 0]/(e1f*e2f) 35 37 36 38 ; 37 39 ; On masque le domaine non physique (points en double) 38 40 ; 39 fm (0, *)= 040 fm (jpi-1, *)= 041 fm[0, *] = 0 42 fm[jpi-1, *] = 0 41 43 ;vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv 42 44 ; specificite ORCA 43 fm (*, jpj-1)= 044 fm (*, jpj-2)= 045 fm[*, jpj-1] = 0 46 fm[*, jpj-2] = 0 45 47 ;^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 46 um = umaskr (*, *, 0)*e1u/e2u47 vm = vmaskr (*, *, 0)*e2v/e1v48 um = umaskr[*, *, 0]*e1u/e2u 49 vm = vmaskr[*, *, 0]*e2v/e1v 48 50 ; 49 51 ; construction des 5 bandes … … 51 53 52 54 ; Sud 53 bande ( 0:N2-1 , 1)= fm*shift(um, 0, 0)55 bande[ 0:N2-1 , 1] = fm*shift(um, 0, 0) 54 56 ; Ouest 55 bande ( 0:N2-1 , 2)= fm*shift(vm, 0, 0)57 bande[ 0:N2-1 , 2] = fm*shift(vm, 0, 0) 56 58 ; Est 57 bande ( 0:N2-1 , 3)= fm*shift(vm,-1, 0)59 bande[ 0:N2-1 , 3] = fm*shift(vm,-1, 0) 58 60 ; Nord 59 bande ( 0:N2-1 , 4)= fm*shift(um, 0,-1)61 bande[ 0:N2-1 , 4] = fm*shift(um, 0,-1) 60 62 ; 61 63 ; 62 64 ; 63 bande (0:N-1 , 1)= 0.65 bande[0:N-1 , 1] = 0. 64 66 ; 65 67 ; elements ci dessous non nuls du fait des CL periodiques E-W 66 68 ; decommenter les 2 lignes si votre domaine est ferme a l''Est et a l''Ouest 67 ; bande (0 , 2)= 0.68 ; bande (N2-1 , 3)= 0.69 ; bande[0 , 2] = 0. 70 ; bande[N2-1 , 3] = 0. 69 71 70 bande (N2-N:N2-1 , 4)= 0.72 bande[N2-N:N2-1 , 4] = 0. 71 73 72 74 ; Centre 73 bande ( * , 0) = - bande(*, 1) - bande(*, 2) - bande(*, 3) - bande(*, 4)75 bande[ * , 0] = - bande[*, 1] - bande[*, 2] - bande[*, 3] - bande[*, 4] 74 76 75 77 ; … … 78 80 79 81 ;; nb d'elements de la matrice creuse 80 Nb = N2+n_elements(where(bande (*, 1:4)NE 0))+182 Nb = N2+n_elements(where(bande[*, 1:4] NE 0))+1 81 83 82 84 ;; definition de la matrice creuse … … 84 86 85 87 ;; elements diagonaux 86 lapla_f.sa[0:N2-1] = bande (0:N2-1, 0)88 lapla_f.sa[0:N2-1] = bande[0:N2-1, 0] 87 89 lapla_f.sa[N2] = 0 88 90 lapla_f.ija[N2] = Nb+1 … … 105 107 ;; creuse ^^^^^^^^ 106 108 107 i = long (0)108 k = long (N2)109 compteur = long (0)109 i = long[0] 110 k = long[N2] 111 compteur = long[0] 110 112 REPEAT BEGIN 111 113 ligne = i/4 … … 142 144 ENDIF 143 145 144 colonne = ligne+decalage (ordre(i MOD 4)-1)146 colonne = ligne+decalage[ordre(i MOD 4)-1] 145 147 ; 146 148 IF (colonne GE 0) AND (colonne LT N2) THEN BEGIN 147 x = bande (ligne, ordre(i MOD 4) )149 x = bande[ligne, ordre[i MOD 4] ] 148 150 IF x NE 0. THEN BEGIN 149 151 k = k+1 … … 160 162 161 163 END 162 -
trunk/tools/bsf/curl2.pro
r165 r169 12 12 FUNCTION curl2, u, v 13 13 ; 14 compile_opt idl2, strictarrsubs 15 ; 14 16 @common 15 17 @com_eg 16 18 17 zu = u*e1u*umaskr (*, *, 0)18 zv = v*e2v*vmaskr (*, *, 0)19 zu = u*e1u*umaskr[*, *, 0] 20 zv = v*e2v*vmaskr[*, *, 0] 19 21 20 22 zcurl = shift(zv, -1, 0)-zv+zu-shift(zu, 0, -1) 21 zcurl = zcurl*fmaskr (*, *, 0)/(e1f*e2f)23 zcurl = zcurl*fmaskr[*, *, 0]/(e1f*e2f) 22 24 23 25 zcurl = boundperio(zcurl, /ff1) -
trunk/tools/bsf/diag_bsf.pro
r165 r169 13 13 PRO diag_bsf_mu, hu, hv, mu 14 14 ; 15 compile_opt idl2, strictarrsubs 16 ; 15 17 @common 16 18 @com_eg … … 20 22 mu = fltarr(14) 21 23 ; 0/ Antartic 22 ; mu (0) = total( (hu*e2u*umaskr(*, *, 0))[149, 5:35])23 mu (0) = total( (hu*e2u*umaskr(*, *, 0))[151, 14:48])24 ; mu[0] = total( (hu*e2u*umaskr[*, *, 0])[149, 5:35]) 25 mu[0] = total( (hu*e2u*umaskr[*, *, 0])[151, 14:48]) 24 26 ; 1/ America 25 ; mu (1) = -total( (hv*e1v*vmaskr(*, *, 0))[110:150, 60])26 mu (1) = -total( (hv*e1v*vmaskr(*, *, 0))[110:150, 60])27 ; mu[1] = -total( (hv*e1v*vmaskr[*, *, 0])[110:150, 60]) 28 mu[1] = -total( (hv*e1v*vmaskr[*, *, 0])[110:150, 60]) 27 29 ; 3/ Australia 28 ; mu (3) = total( (hu*e2u*umaskr(*, *, 0))[17, 40:49])29 mu (3) = total( (hu*e2u*umaskr(*, *, 0))[20, 50:61])30 ; mu[3] = total( (hu*e2u*umaskr[*, *, 0])[17, 40:49]) 31 mu[3] = total( (hu*e2u*umaskr[*, *, 0])[20, 50:61]) 30 32 ; 2/ New Zealand 31 ; mu (2) = mu(3)+ total( (hv*e1v*vmaskr(*, *, 0))[30:40, 25])32 mu (2) = mu(3)+ total( (hv*e1v*vmaskr(*, *, 0))[35:49, 44])33 ; mu[2] = mu[3]+ total( (hv*e1v*vmaskr[*, *, 0])[30:40, 25]) 34 mu[2] = mu[3]+ total( (hv*e1v*vmaskr[*, *, 0])[35:49, 44]) 33 35 ; 4/ Madagascar 34 ; mu (4) = total( (hv*e1v*vmaskr(*, *, 0))[150:159, 40])35 mu (4) = total( (hv*e1v*vmaskr(*, *, 0))[160:164, 56])36 ; mu[4] = total( (hv*e1v*vmaskr[*, *, 0])[150:159, 40]) 37 mu[4] = total( (hv*e1v*vmaskr[*, *, 0])[160:164, 56]) 36 38 ; 5/ New Guinea 37 ; mu (5) = - total( (hu*e2u*umaskr(*, *, 0))[16, 50:55])38 mu (5) = total( (hv*e2v*vmaskr(*, *, 0))[23:33, 61])39 ; mu[5] = - total( (hu*e2u*umaskr[*, *, 0])[16, 50:55]) 40 mu[5] = total( (hv*e2v*vmaskr[*, *, 0])[23:33, 61]) 39 41 ; 6/ Celebes 40 mu (6) = total( (hv*e1v*vmaskr(*, *, 0))[13:22, 64])42 mu[6] = total( (hv*e1v*vmaskr[*, *, 0])[13:22, 64]) 41 43 ; 7/ Borneo 42 ; mu (6) = total( (hv*e1v*vmaskr(*, *, 0))[7:12, 60])43 mu (7) = total( (hv*e1v*vmaskr(*, *, 0))[14:17, 67])44 ; mu[6] = total( (hv*e1v*vmaskr[*, *, 0])[7:12, 60]) 45 mu[7] = total( (hv*e1v*vmaskr[*, *, 0])[14:17, 67]) 44 46 ; 8/ Philippines 45 ; mu (7) = total( (hv*e1v*vmaskr(*, *, 0))[10:16, 90])46 mu (8) = total( (hv*e1v*vmaskr(*, *, 0))[14:22, 86])47 ; mu[7] = total( (hv*e1v*vmaskr[*, *, 0])[10:16, 90]) 48 mu[8] = total( (hv*e1v*vmaskr[*, *, 0])[14:22, 86]) 47 49 ; 9/ Tasmania 48 mu (9) = mu(3)+ total( (hu*e1u*umaskr(*, *, 0))[34, 41:45])50 mu[9] = mu[3]+ total( (hu*e1u*umaskr[*, *, 0])[34, 41:45]) 49 51 ; 10/ Cuba 50 ; mu (8) = - total( (hv*e1v*vmaskr(*, *, 0))[101:127, 86])51 mu (10) = - total( (hv*e1v*vmaskr(*, *, 0))[105:134, 90])52 ; mu[8] = - total( (hv*e1v*vmaskr[*, *, 0])[101:127, 86]) 53 mu[10] = - total( (hv*e1v*vmaskr[*, *, 0])[105:134, 90]) 52 54 ; 11/ Island 53 ; mu (9) = - total( (hv*e1v*vmaskr(*, *, 0))[117:130, 124])54 mu (11) = - total( (hv*e1v*vmaskr(*, *, 0))[132:143, 120])55 ; mu[9] = - total( (hv*e1v*vmaskr[*, *, 0])[117:130, 124]) 56 mu[11] = - total( (hv*e1v*vmaskr[*, *, 0])[132:143, 120]) 55 57 ; 12/ Spitzberg 56 ; mu (10) = - total( (hv*e1v*vmaskr(*, *, 0))[110:120, 137])57 mu (12) = - total( (hv*e1v*vmaskr(*, *, 0))[140:160, 135])58 ; mu[10] = - total( (hv*e1v*vmaskr[*, *, 0])[110:120, 137]) 59 mu[12] = - total( (hv*e1v*vmaskr[*, *, 0])[140:160, 135]) 58 60 ; 13/ Japan 59 mu (13) = total( (hv*e1v*vmaskr(*, *, 0))[25:131, 100])61 mu[13] = total( (hv*e1v*vmaskr[*, *, 0])[25:131, 100]) 60 62 end 61 63 ;+ … … 70 72 ;- 71 73 PRO diag_bsf_lec, bsf 74 ; 72 75 @common 73 76 @com_eg … … 80 83 ile = 'island'+strtrim(string(k), 2) 81 84 z2d = nc_get(fichier, ile, NCDF_DB =hom_idl+'grids/') 82 bsf (*, *, k-1)= shift(z2d, 0, 0)85 bsf[*, *, k-1] = shift(z2d, 0, 0) 83 86 ENDFOR 84 87 end … … 90 93 , FLUX = flux $ 91 94 , _EXTRA=extra 95 ; 96 compile_opt idl2, strictarrsubs 92 97 ; 93 98 @common … … 121 126 IF n_elements(bsfislands) EQ 0 THEN diag_bsf_lec, bsfislands 122 127 123 FOR k = 0, jpisl-1 DO bsf = bsf+mu (k)*bsfislands(*, *, k)128 FOR k = 0, jpisl-1 DO bsf = bsf+mu[k]*bsfislands[*, *, k] 124 129 125 130 return, bsf 126 131 127 END 132 END -
trunk/tools/bsf/inversion.pro
r165 r169 3 3 ; Resolution de l'equation A.X=b ou A est une matrice creuse 4 4 ; et (x,b) sont des champs horizontaux 5 ; 6 ; @todo 7 ; missing build_laplacien_t module 5 8 ; 6 9 ; @history … … 16 19 , TOL=tol $ 17 20 , _EXTRA=extra 21 ; 22 compile_opt idl2, strictarrsubs 18 23 ; 19 24 @common … … 40 45 zx = double(guess) 41 46 42 zr = linbcg(A, zb (*), zx(*), tol = tol, iter = iter, /double, _EXTRA=extra)47 zr = linbcg(A, zb[*], zx[*], tol = tol, iter = iter, /double, _EXTRA=extra) 43 48 print, ' Number of iterations :', iter 44 49 -
trunk/tools/bsf/moyz.pro
r165 r169 15 15 , INT=int 16 16 ; 17 compile_opt idl2, strictarrsubs 18 ; 17 19 @common 18 20 … … 24 26 0 : tt = tt 25 27 1 : tt = tt*(replicate(1, jpi*jpj)#msk) 26 2 : tt = tt*(msk (*)#replicate(1, jpk))28 2 : tt = tt*(msk[*]#replicate(1, jpk)) 27 29 3 : tt = tt*msk 28 30 ENDCASE 29 31 30 32 height = reform( replicate(1, jpi*jpj)#e3t, jpi, jpj, jpk)*tt 31 ; IF s_exp EQ 'GF1' THEN height (*, *, 0) = height(*, *, 0)+ssh33 ; IF s_exp EQ 'GF1' THEN height[*, *, 0] = height[*, *, 0]+ssh 32 34 totz = total( z*height, 3) 33 35 totheight = total( height, 3) … … 37 39 IF n_elements(int) EQ 0 THEN BEGIN 38 40 ind = where(totheight NE 0) 39 moy (ind) = totz(ind) / totheight(ind)41 moy[ind] = totz[ind] / totheight[ind] 40 42 ENDIF 41 43 -
trunk/tools/density_binning/bining.pro
r165 r169 68 68 , TMASK=tmask 69 69 ; 70 compile_opt idl2, strictarrsubs 71 ; 70 72 size3d = size(density) 71 jpi = size3d (1)72 jpj = size3d (2)73 jpk = size3d (3)73 jpi = size3d[1] 74 jpj = size3d[2] 75 jpk = size3d[3] 74 76 75 77 indm = where(tmask eq 0) 76 78 77 x (indm)= 0.78 density (indm)= !VALUES.F_NAN79 x[indm] = 0. 80 density[indm] = !VALUES.F_NAN 79 81 80 82 N_z = jpk … … 117 119 cont_s = fltarr(jpi, jpj, N_s+2) 118 120 119 x_bin (*, *, *)= !VALUES.F_NAN120 depth_bin (*, *, *)= !VALUES.F_NAN121 x_bin[*, *, *] = !VALUES.F_NAN 122 depth_bin[*, *, *] = !VALUES.F_NAN 121 123 122 124 ; … … 124 126 ; 125 127 x_content = fltarr(jpi, jpj, jpk) 126 x_content (*, *, 0)= 0.128 x_content[*, *, 0] = 0. 127 129 ; FOR k = 0, jpk-1 DO BEGIN 128 ; print, k, z_zw (k), e3t(k), x(4, 30, k)130 ; print, k, z_zw[k], e3t[k], x[4, 30, k] 129 131 ; ENDFOR 130 132 FOR k = 1, jpk-1 DO BEGIN 131 x_content (*, *, k) = x_content(*, *, k-1)$132 ; + e3t (k-1) * x(*, *, k-1)133 + ( z_zw(k)-z_zw(k-1) ) * x(*, *, k)134 ; print, k, z_zw (k), e3t(k-1), ( z_zw(k)-z_zw(k-1) ), ( z_zw(k)-z_zw(k-1) )-e3t(k-1)133 x_content[*, *, k] = x_content[*, *, k-1] $ 134 ; + e3t[k-1] * x[*, *, k-1] 135 + [ z_zw[k]-z_zw[k-1] ] * x[*, *, k] 136 ; print, k, z_zw[k], e3t[k-1], [ z_zw[k]-z_zw[k-1] ], [ z_zw[k]-z_zw[k-1] ]-e3t[k-1] 135 137 ENDFOR 136 138 ; … … 144 146 ; Indices des points T dans l ocean 145 147 ; 146 i_ocean = where(tmask (i, j, *)EQ 1)147 148 z_s (*)= !VALUES.F_NAN148 i_ocean = where(tmask[i, j, *] EQ 1) 149 150 z_s[*] = !VALUES.F_NAN 149 151 c_s = c_s*0. 150 x_s (*)= !VALUES.F_NAN152 x_s[*] = !VALUES.F_NAN 151 153 152 154 IF i_ocean[0] NE -1 THEN BEGIN ; on n entre que si il y a des points ocean 153 155 154 i_bottom = i_ocean (n_elements(i_ocean)-1)155 156 z_s (0)= 0157 z_s (N_s+1) = z_zw(i_bottom)158 159 c_s (0)= 0160 c_s (N_s+1) = x_content(i, j, jpk-1)156 i_bottom = i_ocean[n_elements(i_ocean)-1] 157 158 z_s[0] = 0 159 z_s[N_s+1] = z_zw[i_bottom] 160 161 c_s[0] = 0 162 c_s[N_s+1] = x_content[i, j, jpk-1] 161 163 162 s_z (*) = density(i, j, *)163 c_z (*) = x_content(i, j, *)164 s_z[*] = density[i, j, *] 165 c_z[*] = x_content[i, j, *] 164 166 165 167 ; exctraction d'un sous profil strictement croissant 166 mini = min(s_z (i_ocean))167 maxi = max(s_z (i_ocean))168 i_min = where(s_z (i_ocean)EQ mini)169 i_max = where(s_z (i_ocean)EQ maxi)168 mini = min(s_z[i_ocean]) 169 maxi = max(s_z[i_ocean]) 170 i_min = where(s_z[i_ocean] EQ mini) 171 i_max = where(s_z[i_ocean] EQ maxi) 170 172 ; on prend le plus grand des indices min 171 173 ; et le plus petit des indices max 172 i_min = i_min (n_elements(i_min)-1)174 i_min = i_min[n_elements(i_min)-1] 173 175 i_max = i_max[0] 174 176 … … 176 178 ; l isopycne est mise en surface (z_s=0) 177 179 ; 178 ind = where(s_s LT s_z (i_min))180 ind = where(s_s LT s_z[i_min]) 179 181 IF ind[0] NE -1 THEN BEGIN 180 z_s (ind+1)= 0181 c_s (ind+1)= 0182 z_s[ind+1] = 0 183 c_s[ind+1] = 0 182 184 ENDIF 183 185 184 186 ; Si la valeur du niveau (s_s) est plus elevee que la densite du fond, 185 ; l isopycne est mise au fond (z_s=z_zw (i_bottom))186 ; 187 ind = where(s_s GT s_z (i_max))187 ; l isopycne est mise au fond (z_s=z_zw[i_bottom]) 188 ; 189 ind = where(s_s GT s_z[i_max]) 188 190 IF ind[0] NE -1 THEN BEGIN 189 z_s (ind+1) = z_s(N_s+1)190 c_s (ind+1) = c_s(N_s+1)191 z_s[ind+1] = z_s[N_s+1] 192 c_s[ind+1] = c_s[N_s+1] 191 193 ENDIF 192 194 193 ind = where( (s_s GE s_z (i_min)) AND (s_s LE s_z(i_max)) )195 ind = where( (s_s GE s_z[i_min]) AND (s_s LE s_z[i_max]) ) 194 196 IF ind[0] NE -1 THEN BEGIN 195 197 196 i_profil = i_ocean (i_min:i_max)198 i_profil = i_ocean[i_min:i_max] 197 199 198 200 ; print, i, j, n_elements(i_profil) 199 201 200 z_s (ind+1) = interpol(z_zt(i_profil), s_z(i_profil), s_s(ind))202 z_s[ind+1] = interpol(z_zt[i_profil], s_z[i_profil], s_s[ind]) 201 203 202 204 ; … … 204 206 ; 205 207 IF n_elements(i_profil) GT 2 THEN BEGIN 206 c_s (ind+1) = spline(z_zw(i_profil), c_z(i_profil), z_s(ind+1), 1)208 c_s[ind+1] = spline(z_zw[i_profil], c_z[i_profil], z_s[ind+1], 1) 207 209 ENDIF ELSE BEGIN 208 c_s (ind+1) = interpol(c_z(i_profil), z_zw(i_profil), z_s(ind+1))210 c_s[ind+1] = interpol(c_z[i_profil], z_zw[i_profil], z_s[ind+1]) 209 211 ENDELSE 210 212 ; … … 212 214 ; me semble moins propre. Je la donne en remarque. 213 215 ; 214 ; c_s (ind+1) = interpol(c_z(i_profil), z_zw(i_profil), z_s(ind+1))215 216 x_s (ind+1) = (c_s(ind+2)-c_s(ind+1))/(z_s(ind+2)-z_s(ind+1))216 ; c_s[ind+1] = interpol(c_z[i_profil], z_zw[i_profil], z_s[ind+1]) 217 218 x_s[ind+1] = (c_s[ind+2]-c_s[ind+1])/(z_s[ind+2]-z_s[ind+1]) 217 219 218 220 ENDIF 219 221 220 ; x_s (0) = c_s(1)/z_s(1)222 ; x_s[0] = c_s[1]/z_s[1] 221 223 222 224 ENDIF 223 225 224 depth_bin (i, j, *)= z_s225 x_bin (i, j, *)= x_s226 depth_bin[i, j, *] = z_s 227 x_bin [i, j, *] = x_s 226 228 ENDFOR 227 229 ENDFOR 228 230 229 231 ; mask depth_bin with undefined values of x_bin 230 depth_bin (where(finite(x_bin, /nan) EQ 1))= !VALUES.F_NAN232 depth_bin[where(finite(x_bin, /nan) EQ 1)] = !VALUES.F_NAN 231 233 232 234 ; print, min(x_bin(where (abs(x_bin) LT 1.e19))), max(x_bin(where (abs(x_bin) LT 1.e19))) 233 235 ; print, min(depth_bin(where (abs(depth_bin) LT 1.e19))), max(depth_bin(where (abs(depth_bin) LT 1.e19))) 234 236 235 END 237 END -
trunk/tools/density_binning/bining2_bak.pro
r165 r169 63 63 ; 64 64 ;- 65 PRO bining2 $65 PRO bining2_bak $ 66 66 , density $ 67 67 , x1 $ … … 76 76 , TMASK=tmask 77 77 ; 78 compile_opt idl2, strictarrsubs 79 ; 78 80 size3d = size(density) 79 jpi = size3d (1)80 jpj = size3d (2)81 jpk = size3d (3)81 jpi = size3d[1] 82 jpj = size3d[2] 83 jpk = size3d[3] 82 84 83 85 indm = where(tmask eq 0) 84 86 85 x1 (indm)= 0.86 density (indm)= !VALUES.F_NAN87 x1[indm] = 0. 88 density[indm] = !VALUES.F_NAN 87 89 88 90 N_z = jpk … … 129 131 x1_bin = fltarr(jpi, jpj, N_s+1) 130 132 131 x1_bin (*, *, *)= !VALUES.F_NAN132 depth_bin (*, *, *)= !VALUES.F_NAN133 thick_bin (*, *, *)= !VALUES.F_NAN133 x1_bin[*, *, *] = !VALUES.F_NAN 134 depth_bin[*, *, *] = !VALUES.F_NAN 135 thick_bin[*, *, *] = !VALUES.F_NAN 134 136 135 137 ; … … 137 139 ; 138 140 x1_content = fltarr(jpi, jpj, jpk) 139 x1_content (*, *, 0) = e3t(0) * x1(*, *, 0) * tmask (*, *, 0)140 FOR k = 1, jpk-1 DO x1_content (*, *, k) = x1_content(*, *, k-1)$141 + e3t (k) * x1(*, *, k)* tmask (*, *, k)141 x1_content[*, *, 0] = e3t[0] * x1[*, *, 0] * tmask [*, *, 0] 142 FOR k = 1, jpk-1 DO x1_content[*, *, k] = x1_content[*, *, k-1] $ 143 + e3t[k] * x1[*, *, k]* tmask [*, *, k] 142 144 143 145 … … 150 152 ; Indices des points T dans l ocean 151 153 ; 152 i_ocean = where(tmask (i, j, *)EQ 1)154 i_ocean = where(tmask[i, j, *] EQ 1) 153 155 154 156 z_s = z_s*0. … … 159 161 IF i_ocean[0] NE -1 THEN BEGIN ; on n entre que si il y a des points ocean 160 162 161 i_bottom = i_ocean (n_elements(i_ocean)-1)162 163 z_s (N_s) = z_zw(i_bottom)164 c1_s (N_s) = x1_content(i, j, jpk-1)163 i_bottom = i_ocean[n_elements(i_ocean)-1] 164 165 z_s[N_s] = z_zw[i_bottom] 166 c1_s[N_s] = x1_content[i, j, jpk-1] 165 167 166 s_z (*) = density(i, j, *)167 c1_z (*) = x1_content(i, j, *)168 s_z[*] = density[i, j, *] 169 c1_z[*] = x1_content[i, j, *] 168 170 169 171 ; extraction d'un sous profil strictement croissant 170 mini = min(s_z (i_ocean))171 maxi = max(s_z (i_ocean))172 i_min = where(s_z (i_ocean)EQ mini)173 i_max = where(s_z (i_ocean)EQ maxi)172 mini = min(s_z[i_ocean]) 173 maxi = max(s_z[i_ocean]) 174 i_min = where(s_z[i_ocean] EQ mini) 175 i_max = where(s_z[i_ocean] EQ maxi) 174 176 ; on prend le plus grand des indices min 175 177 ; et le plus petit des indices max 176 i_min = i_min (n_elements(i_min)-1)178 i_min = i_min[n_elements(i_min)-1] 177 179 i_max = i_max[0] 178 180 … … 183 185 ; l isopycne est mise en surface (z_s=0) 184 186 ; 185 ind = where(s_s LT s_z (i_min))187 ind = where(s_s LT s_z[i_min]) 186 188 IF ind[0] NE -1 THEN BEGIN 187 z_s (ind)= 0188 c1_s (ind)= 0189 z_s[ind] = 0 190 c1_s[ind] = 0 189 191 ENDIF 190 192 191 193 ; Si la valeur du niveau (s_s) est plus elevee que la densite du fond, 192 ; l isopycne est mise au fond (z_s=z_zw (i_bottom))193 ; 194 ind = where(s_s GT s_z (i_max))194 ; l isopycne est mise au fond (z_s=z_zw[i_bottom]) 195 ; 196 ind = where(s_s GT s_z[i_max]) 195 197 IF ind[0] NE -1 THEN BEGIN 196 z_s (ind) = z_s(N_s)197 c1_s (ind) = c1_s(N_s)198 z_s[ind] = z_s[N_s] 199 c1_s[ind] = c1_s[N_s] 198 200 ENDIF 199 201 ; cas general 200 ind = where( (s_s GE s_z (i_min)) AND (s_s LE s_z(i_max)) )202 ind = where( (s_s GE s_z[i_min]) AND (s_s LE s_z[i_max]) ) 201 203 IF ind[0] NE -1 THEN BEGIN 202 204 203 i_profil = i_ocean (i_min:i_max)204 205 z_s (ind) = interpol(z_zt(i_profil), s_z(i_profil), s_s(ind))205 i_profil = i_ocean[i_min:i_max] 206 207 z_s[ind] = interpol(z_zt[i_profil], s_z[i_profil], s_s[ind]) 206 208 207 209 ; … … 211 213 ; 212 214 IF n_elements(i_profil) GT 2 THEN BEGIN 213 c1_s (ind) = spline(z_zw(i_profil), c1_z(i_profil), z_s(ind), 1)214 ; c1_s (ind) = interpol(c1_z(i_profil), z_zw(i_profil), z_s(ind))215 c1_s[ind] = spline(z_zw[i_profil], c1_z[i_profil], z_s[ind], 1) 216 ; c1_s[ind] = interpol(c1_z[i_profil], z_zw[i_profil], z_s[ind]) 215 217 ENDIF ELSE BEGIN 216 c1_s (ind) = interpol(c1_z(i_profil), z_zw(i_profil), z_s(ind))218 c1_s[ind] = interpol(c1_z[i_profil], z_zw[i_profil], z_s[ind]) 217 219 ENDELSE 218 220 ; 219 221 ; 220 x1_s (ind) = (c1_s(ind+1)-c1_s(ind))/(z_s(ind+1)-z_s(ind))222 x1_s[ind] = (c1_s[ind+1]-c1_s[ind])/(z_s[ind+1]-z_s[ind]) 221 223 222 224 ENDIF 223 225 ENDIF 224 226 ; 225 depth_bin (i, j, *)= z_s226 thick_bin (i, j, 0) = z_s(0)227 thick_bin (i, j, 1:N_s) = z_s(1:N_s)-z_s(0:N_s-1)228 x1_bin (i, j, *)= x1_s227 depth_bin [i, j, *] = z_s 228 thick_bin [i, j, 0] = z_s[0] 229 thick_bin [i, j, 1:N_s] = z_s[1:N_s]-z_s[0:N_s-1] 230 x1_bin [i, j, *] = x1_s 229 231 ; 230 232 ENDFOR 231 233 ENDFOR 232 234 ; mask depth_bin with undefined values of x_bin 233 depth_bin (where(finite(x1_bin, /nan) EQ 1))= !VALUES.F_NAN234 thick_bin (where(finite(x1_bin, /nan) EQ 1))= !VALUES.F_NAN235 depth_bin[where(finite(x1_bin, /nan) EQ 1)] = !VALUES.F_NAN 236 thick_bin[where(finite(x1_bin, /nan) EQ 1)] = !VALUES.F_NAN 235 237 236 238 237 END 239 END -
trunk/tools/density_binning/binning_neutral_and_co/bin_velocity_new.pro
r166 r169 9 9 , VV=vv $ 10 10 , TT=tt 11 ; 12 compile_opt idl2, strictarrsubs 11 13 ; 12 14 @common … … 42 44 ; z_zw = gdepw (old) 43 45 z_zw=shift(gdepw,-1) ; W (k=0 -> 10m) 44 z_zw (jpk-1)=gdepw(jpk-1)46 z_zw[jpk-1]=gdepw[jpk-1] 45 47 46 48 ; s or partial step case: … … 48 50 fich_msh = 'mesh_mask_NDV1.nc' 49 51 z3d_zt = read_ncdf('e3v_ps',0,/timestep,iodir=ioMESH,/nostruct,/tout,filename=fich_msh) 50 z3d_zt (where(z3d_zt eq 1.e20))=500.52 z3d_zt[where(z3d_zt eq 1.e20)]=500. 51 53 z3d_zw = total(z3d_zt,3,/cumulative) ; W (k=0 -> 10m) 52 54 z3d_zt = (shift(z3d_zw,0,0,1)+z3d_zw)*0.5 53 z3d_zt (*,*,0) = z3d_zw(*,*,0)*0.5 ; T (k=0 -> 5m)55 z3d_zt[*,*,0] = z3d_zw[*,*,0]*0.5 ; T (k=0 -> 5m) 54 56 55 57 ; Profiles along density level (suffixe _s) … … 62 64 ; vertical content of x1 (ensure the integrale conservation) 63 65 x1_content = fltarr(jpi, jpj, jpk) 64 x1_content (*, *, 0) = e3t(0) * x1(*, *, 0) * xmask(*,*,0)65 FOR k = 1, jpk-1 DO x1_content (*, *, k) = x1_content(*, *, k-1)$66 + e3t (k) * x1(*, *, k) * xmask(*,*,k)66 x1_content[*, *, 0] = e3t[0] * x1[*, *, 0] * xmask[*,*,0] 67 FOR k = 1, jpk-1 DO x1_content[*, *, k] = x1_content[*, *, k-1] $ 68 + e3t[k] * x1[*, *, k] * xmask[*,*,k] 67 69 ; 68 70 ; Loop over the horizontal domain 2D … … 72 74 73 75 ; depth at T and W points 74 z_zt (*) = z3d_zt(i,j,*)75 z_zw (*) = z3d_zw(i,j,*)76 z_zt[*] = z3d_zt[i,j,*] 77 z_zw[*] = z3d_zw[i,j,*] 76 78 77 79 ; 78 80 ; Indices des points T dans l ocean 79 81 ; 80 i_ocean = where(xmask (i, j, *)EQ 1)82 i_ocean = where(xmask[i, j, *] EQ 1) 81 83 z_s = z_s*0. 82 84 c1_s = c1_s*0. … … 84 86 IF (i_ocean[0] NE -1) THEN BEGIN ; on n'entre que si il y a des points ocean 85 87 86 s_z (*) = density(i, j, *)87 c1_z (*) = x1_content(i, j, *)88 s_z[*] = density[i, j, *] 89 c1_z[*] = x1_content[i, j, *] 88 90 ; critere supplementaire a imposer sur le profil pour eviter les cas 89 91 ; pathologiques en attendant d'ecrire une vraie commande d'extraction 90 92 ; de progils stritement croissant. Il s'agit donc d'un test adhoc. 91 IF s_z (0) LT s_z(i_ocean(n_elements(i_ocean)-1))THEN BEGIN93 IF s_z[0] LT s_z[i_ocean[n_elements(i_ocean)-1]] THEN BEGIN 92 94 ;------------------------------------------------------------------------ 93 95 ; controle si le profil est bien strictement croissant (pour l'instant 94 ; inutilis (avis aux amateurs)96 ; inutilisé (avis aux amateurs) 95 97 ; 96 ; ds = (shift(s_z, -1)-s_z)(i_ocean (0:n_elements(i_ocean)-2))97 ; ds (0)= 098 ; ds = (shift(s_z, -1)-s_z)(i_ocean[0:n_elements(i_ocean)-2]) 99 ; ds[0] = 0 98 100 ; ind = where(ds LT 0.) 99 ; croissant = ind (0)EQ -1101 ; croissant = ind[0] EQ -1 100 102 ;------------------------------------------------------------------------ 101 i_bottom = i_ocean (n_elements(i_ocean)-1)102 z_s (N_s) = z_zw(i_bottom)103 c1_s (N_s) = x1_content(i, j, jpk-1)103 i_bottom = i_ocean[n_elements(i_ocean)-1] 104 z_s[N_s] = z_zw[i_bottom] 105 c1_s[N_s] = x1_content[i, j, jpk-1] 104 106 105 107 ; extraction d'un sous profil strictement croissant 106 mini = min(s_z (i_ocean))107 maxi = max(s_z (i_ocean))108 i_min = where(s_z (i_ocean)EQ mini)109 i_max = where(s_z (i_ocean)EQ maxi)108 mini = min(s_z[i_ocean]) 109 maxi = max(s_z[i_ocean]) 110 i_min = where(s_z[i_ocean] EQ mini) 111 i_max = where(s_z[i_ocean] EQ maxi) 110 112 ; on prend le plus grand des indices min 111 113 ; et le plus petit des indices max 112 ;gr i_min = i_min (n_elements(i_min)-1)114 ;gr i_min = i_min[n_elements(i_min)-1] 113 115 ;gr i_max = i_max[0] 114 116 i_min = i_min[0] 115 i_max = i_max (n_elements(i_min)-1)117 i_max = i_max[n_elements(i_min)-1] 116 118 ; IF i_max GE jpk-1 THEN print, i, j, i_max 117 119 ; IF i_min LE 1 THEN print, i, j, i_min … … 120 122 ; l isopycne est mise en surface (z_s=0) 121 123 ; 122 ind = where(s_s LT s_z (i_min))124 ind = where(s_s LT s_z[i_min]) 123 125 IF ind[0] NE -1 THEN BEGIN &$ 124 z_s (ind)= 0 &$125 c1_s (ind)= 0 &$126 z_s[ind] = 0 &$ 127 c1_s[ind] = 0 &$ 126 128 ENDIF 127 129 ; Si la valeur du niveau (s_s) est plus elevee que la densite du fond, 128 ; l isopycne est mise au fond (z_s=z_zw (i_bottom))130 ; l isopycne est mise au fond (z_s=z_zw[i_bottom]) 129 131 ; 130 ind = where(s_s GT s_z (i_max))132 ind = where(s_s GT s_z[i_max]) 131 133 IF ind[0] NE -1 THEN BEGIN &$ 132 z_s (ind) = z_s(N_s)&$133 c1_s (ind) = c1_s(N_s)&$134 z_s[ind] = z_s[N_s] &$ 135 c1_s[ind] = c1_s[N_s] &$ 134 136 ENDIF 135 137 ; 136 ds =(shift(s_z, -1)-s_z) (i_ocean(0:n_elements(i_ocean)-2))138 ds =(shift(s_z, -1)-s_z)[i_ocean[0:n_elements(i_ocean)-2]] 137 139 ind_c = where(ds LT 0.) 138 ; croissant = ind_c (0)EQ -1140 ; croissant = ind_c[0] EQ -1 139 141 ; IF ( i_min GT i_max ) then begin 140 142 ; print, 'bug ',i_min, i_max,' en i,j=', i,j 141 143 ; endif 142 IF ( ind_c (0)NE -1 ) then begin144 IF ( ind_c[0] NE -1 ) then begin 143 145 print, 'bug en i,j=', i,j,' ind_c = ', ind_c 144 146 stop 145 147 endif 146 148 ; IF ( i_min LE i_max ) then begin 147 IF ( ind_c (0)EQ -1 ) then begin148 ind = where( (s_s GE s_z (i_min)) AND (s_s LT s_z(i_max)) )149 IF ( ind_c[0] EQ -1 ) then begin 150 ind = where( (s_s GE s_z[i_min]) AND (s_s LT s_z[i_max]) ) 149 151 IF ind[0] NE -1 THEN BEGIN 150 152 IF ( n_elements(ind) NE 1 ) THEN BEGIN 151 153 152 i_profil = i_ocean (i_min:i_max)153 z_s (ind) = interpol(z_zt(i_profil), s_z(i_profil), s_s(ind))154 i_profil = i_ocean[i_min:i_max] 155 z_s[ind] = interpol(z_zt[i_profil], s_z[i_profil], s_s[ind]) 154 156 155 157 ; 156 158 ; j'utilise la fonction spline pour interpoler le contenu 157 159 ; 158 c1_s (ind) = spline(z_zw(i_profil), c1_z(i_profil), z_s(ind), 1)160 c1_s[ind] = spline(z_zw[i_profil], c1_z[i_profil], z_s[ind], 1) 159 161 ; 160 162 ; l'interpolation lineaire marche aussi. Elle donne des resultats differents. Elle 161 163 ; me semble moins propre. Je la donne en remarque. 162 164 ; 163 ; c_s (ind+1) = interpol(c_z(i_profil), z_zw(i_profil), z_s(ind+1))165 ; c_s[ind+1] = interpol(c_z[i_profil], z_zw[i_profil], z_s[ind+1]) 164 166 ; 165 167 ; Remarque : on ne divise pas par l'epaisseur de la couche … … 168 170 IF ( n_elements(ind) EQ 1 ) THEN BEGIN 169 171 ; print, 'one density bin en i,j=', i,j,' ind = ', ind 170 c1_s (ind) = c1_z(jpk-1)172 c1_s[ind] = c1_z[jpk-1] 171 173 ENDIF 172 174 ; 173 x1_s(ind[0]-1) = (c1_s(ind[0])-c1_s(ind[0]-1));/(z_s(ind+1)-z_s(ind)) 174 x1_s(ind) = (c1_s(ind+1)-c1_s(ind));/(z_s(ind+1)-z_s(ind)) 175 ; x1_s(ind[0]-1) = (c1_s(ind[0])-c1_s(ind[0]-1))/(z_s(ind+1)-z_s(ind)) 176 ; x1_s(ind) = (c1_s(ind+1)-c1_s(ind))/(z_s(ind+1)-z_s(ind)) 175 x1_s[ind[0]-1] = (c1_s[ind[0]]-c1_s[ind[0]-1]);/(z_s[ind+1]-z_s[ind]) 176 x1_s[ind] = (c1_s[ind+1]-c1_s[ind]);/(z_s[ind+1)-z_s[ind]) 177 177 ENDIF 178 178 endif 179 179 ENDIF 180 180 ENDIF 181 x1_bin (i, j, *)= x1_s181 x1_bin[i, j, *] = x1_s 182 182 183 183 ENDFOR -
trunk/tools/density_binning/binning_neutral_and_co/bn2.pro
r168 r169 46 46 ;- 47 47 FUNCTION bn2, tn, sn, neos 48 ; 49 compile_opt idl2, strictarrsubs 48 50 ; 49 51 @common … … 136 138 ; ------------------------ 137 139 ; bn^2=0. at first vertical point and jk=jpk 138 z (*, *, 0)= 0139 IF (gdep (nz-1) eq gdepw(jpk-1)) THEN z(*, *, nz-1)= 0140 z[*, *, 0] = 0 141 IF (gdep[nz-1] eq gdepw[jpk-1]) THEN z[*, *, nz-1] = 0 140 142 sortie: 141 143 return, z -
trunk/tools/density_binning/binning_neutral_and_co/eos.pro
r168 r169 53 53 ;- 54 54 FUNCTION eos, t, s, sigma 55 ; 56 compile_opt idl2, strictarrsubs 55 57 ; 56 58 @common -
trunk/tools/density_binning/binning_neutral_and_co/eos_neutral.pro
r168 r169 38 38 FUNCTION eos_neutral, t, s 39 39 ; 40 compile_opt idl2, strictarrsubs 41 ; 40 42 @common 41 43 ; -
trunk/tools/density_binning/binning_neutral_and_co/fig_msfs_new.pro
r168 r169 8 8 ; 9 9 ;- 10 PRO fig_msf_sigma 11 10 PRO fig_msfs_new 11 ; 12 compile_opt idl2, strictarrsubs 13 ; 12 14 @common 13 15 ; 14 16 ;meshmask 15 17 @0_initorca2_NDV1 … … 46 48 ; 47 49 ; define vmask for each sub basin (mask in the ACC: MSF computation) 48 vglomsk = shift(tmask (*,*,0), 0, -1) *tmask(*,*,0)50 vglomsk = shift(tmask[*,*,0], 0, -1) *tmask[*,*,0] 49 51 atlmsk=glo_atl*glo_noACC 50 52 vatlmsk = shift(atlmsk, 0, -1) * atlmsk … … 104 106 105 107 s_v=sig_min+findgen(n_sig)*sig_del 106 s_v (0)= 0.107 s_v (n_sig-1)= 100.108 s_v[0] = 0. 109 s_v[n_sig-1] = 100. 108 110 109 111 … … 223 225 ;@initorca2_o 224 226 @0_initorca2_NDV1 227 ; 228 END -
trunk/tools/density_binning/binning_neutral_and_co/fsalbt.pro
r168 r169 31 31 FUNCTION fsalbt, pft, pfs, pfh 32 32 ; 33 compile_opt idl2, strictarrsubs 34 ; 33 35 ; ratio alpha/beta 34 36 ; ================ -
trunk/tools/density_binning/binning_neutral_and_co/fsbeta.pro
r168 r169 31 31 FUNCTION fsbeta, pft, pfs, pfh 32 32 ; 33 compile_opt idl2, strictarrsubs 34 ; 33 35 z = $ 34 36 ( ( -0.415613e-09 * pft + 0.555579e-07 ) * pft $ -
trunk/tools/density_binning/binning_neutral_and_co/msfs.pro
r165 r169 24 24 , ORCA=orca 25 25 ; 26 compile_opt idl2, strictarrsubs 27 ; 26 28 @common 27 29 ;------------------------- … … 31 33 ; density = vmask()*sig 32 34 vvmask = vmask() 33 density (*,jpj-1,*)=0.e0 ; set to zero the northern j-line35 density[*,jpj-1,*]=0.e0 ; set to zero the northern j-line 34 36 density=density*vmask() 35 37 ;------------------------- … … 46 48 ; level 20 2.5 Sv 47 49 ; level 21 1.0 Sv 48 vvf (167,100, 0:13) = vvf(167,100, 0:13)+0.8e+6/14.*vvmask(167,100, 0:13)49 vvf (167,100,14:19) = vvf(167,100,14:19)+0.7e+6/6. *vvmask(167,100,14:19)50 vvf (167,100, 20 ) = vvf(167,100, 20 )+2.5e+6 *vvmask(167,100, 20 )51 vvf (167,100, 21 ) = vvf(167,100, 21 )+0.8e+6/14.*vvmask(167,100, 21 )50 vvf[167,100, 0:13] = vvf[167,100, 0:13]+0.8e+6/14.*vvmask[167,100, 0:13] 51 vvf[167,100,14:19] = vvf[167,100,14:19]+0.7e+6/6. *vvmask[167,100,14:19] 52 vvf[167,100, 20 ] = vvf[167,100, 20 ]+2.5e+6 *vvmask[167,100, 20 ] 53 vvf[167,100, 21 ] = vvf[167,100, 21 ]+0.8e+6/14.*vvmask[167,100, 21 ] 52 54 ;------------------------- 53 55 if keyword_set(mask) then begin &$ … … 58 60 ; crude summation 59 61 ;------------------------- 60 N_s = n_elements(s_s) ; number of density levels61 x1_bin = fltarr(jpj, N_s) ; output array (transport function of density)62 x1_bin (*,*)= 0.e063 FOR is = N_s-1,0,-1 DO begin &$64 dmsk= density GT s_s (is)&$65 x1_bin (*,is)= total( total(vvf*dmsk,1),2) &$62 n_s = n_elements(s_s) ; number of density levels 63 x1_bin = fltarr(jpj, n_s) ; output array (transport function of density) 64 x1_bin[*,*] = 0.e0 65 FOR is = n_s-1,0,-1 DO begin &$ 66 dmsk= density GT s_s[is] &$ 67 x1_bin[*,is] = total( total(vvf*dmsk,1),2) &$ 66 68 ENDFOR &$ 67 69 msfs=x1_bin *1.e-6 -
trunk/tools/density_binning/binning_neutral_and_co/npc.pro
r168 r169 30 30 FUNCTION npc, s3d 31 31 ; 32 compile_opt idl2, strictarrsubs 33 ; 32 34 @common 33 35 ;; … … 52 54 ; 53 55 ; Indices des points T dans l ocean 54 i_ocean = where(tmask (i,j,*)EQ 1)56 i_ocean = where(tmask[i,j,*] EQ 1) 55 57 ; 56 58 IF (i_ocean[0] NE -1) THEN BEGIN ; on n'entre que si il y a des points ocean 57 59 ; 58 60 ; density profil 59 s_z (*)= s3d(i,j,*)60 ;s_z (*) = rho(i,j,*)61 s_z[*]= s3d[i,j,*] 62 ;s_z[*] = rho[i,j,*] 61 63 ; 62 64 ; 1. Static instability pointer 63 65 ; ----------------------------- 64 66 ; 65 ds =(shift(s_z,-1)-s_z) (i_ocean(0:n_elements(i_ocean)-2))67 ds =(shift(s_z,-1)-s_z)[i_ocean[0:n_elements(i_ocean)-2]] 66 68 ind_c = where(ds LT 0.) 67 69 ; 68 70 ; 2. Vertical mixing for each instable portion of the density profil 69 71 ; 70 IF ( ind_c (0)NE -1 ) THEN BEGIN72 IF ( ind_c[0] NE -1 ) THEN BEGIN 71 73 ncompt=ncompt+1 72 74 ; print, 'static instability at i,j=', i,j … … 77 79 ; vertical iteration 78 80 jiter = 0 79 WHILE ( (ind_c (0)NE -1) AND (jiter LT jpk-1) ) DO BEGIN &$81 WHILE ( (ind_c[0] NE -1) AND (jiter LT jpk-1) ) DO BEGIN &$ 80 82 jiter = jiter+1 &$ 81 83 ; ikup : the first static instability from the sea surface 82 ikup = ind_c (0)&$84 ikup = ind_c[0] &$ 83 85 ; the density profil is instable below ikup 84 86 ; ikdown : bottom of the instable portion of the density profil 85 87 ; search of ikdown and vertical mixing from ikup to ikdown 86 ze3tot= e3t (ikup)&$87 zraua = s_z (ikup)&$88 ze3tot= e3t[ikup] &$ 89 zraua = s_z[ikup] &$ 88 90 jkdown = ikup+1 &$ 89 91 ; 90 WHILE (jkdown LE ikbot AND zraua GT s_z (jkdown)) DO BEGIN &$91 ze3dwn = e3t (jkdown)&$92 WHILE (jkdown LE ikbot AND zraua GT s_z[jkdown] ) DO BEGIN &$ 93 ze3dwn = e3t[jkdown] &$ 92 94 ze3tot = ze3tot+ze3dwn &$ 93 zraua = ( zraua*(ze3tot-ze3dwn) + s_z (jkdown)*ze3dwn )/ze3tot &$95 zraua = ( zraua*(ze3tot-ze3dwn) + s_z[jkdown]*ze3dwn )/ze3tot &$ 94 96 jkdown=jkdown+1 &$ 95 97 ; print, jkdown, zraua &$ … … 97 99 ; 98 100 FOR jkp = ikup,jkdown-1 DO BEGIN &$ 99 s_z (jkp)= zraua &$101 s_z[jkp] = zraua &$ 100 102 ENDFOR &$ 101 ds =(shift(s_z, -1)-s_z) (i_ocean(0:n_elements(i_ocean)-2))&$103 ds =(shift(s_z, -1)-s_z)[i_ocean[0:n_elements(i_ocean)-2]] &$ 102 104 ind_c = where(ds LT 0.) &$ 103 105 ENDWHILE 104 106 ENDIF 105 107 ; save the modifications 106 rhos (i,j,*) = s_z(*)108 rhos[i,j,*] = s_z[*] 107 109 ; 108 110 ; <<-- no more static instability on slab jj -
trunk/tools/density_binning/density_bin_IDL_gm/bin_sigma.pro
r168 r169 25 25 , BOXZOOM=boxzoom $ 26 26 , ALL_DATA=all_data 27 ; 28 compile_opt idl2, strictarrsubs 27 29 ; 28 30 @common … … 88 90 89 91 CASE pfildi.name OF 90 '@@vodeptht': nfild = depth_bin (*, *, 0:n_sig-1)91 '@@vosigthi': nfild = thick_bin (*, *, 0:n_sig-1)92 '@@vodeptht': nfild = depth_bin[*, *, 0:n_sig-1] 93 '@@vosigthi': nfild = thick_bin[*, *, 0:n_sig-1] 92 94 '@@vosigvol': BEGIN 93 95 nfild = vol_bin 94 96 final_dim = 1 95 97 ; e12t = e1t*e2t 96 ; e12t3 = reform((reform(e12t (premierx:dernierx, premiery:derniery), nx*ny))#replicate(1,n_sig +2), nx, ny, n_sig+2)97 ; nfild = (thick_bin*e12t3) (*, *, 0:n_sig-1)98 ; e12t3 = reform((reform(e12t[premierx:dernierx, premiery:derniery], nx*ny))#replicate(1,n_sig +2), nx, ny, n_sig+2) 99 ; nfild = (thick_bin*e12t3)[*, *, 0:n_sig-1] 98 100 99 101 END 100 ELSE: nfild = f_bin (*, *, 0:n_sig-1)102 ELSE: nfild = f_bin[*, *, 0:n_sig-1] 101 103 ENDCASE 102 104 … … 108 110 109 111 CASE pfildi.name OF 110 '@@vodeptht': nfild = depth_bin (*, *, index)111 '@@vosigthi': nfild = thick_bin (*, *, index)112 '@@vodeptht': nfild = depth_bin[*, *, index] 113 '@@vosigthi': nfild = thick_bin[*, *, index] 112 114 '@@vosigvol': BEGIN 113 115 nfild = vol_bin 114 116 final_dim = 1 115 117 ; e12t = e1t*e2t 116 ; e12t3 = reform((reform(e12t (premierx:dernierx, premiery:derniery), nx*ny))#replicate(1,n_sig +2), nx, ny, n_sig+2)117 ; nfild = ((thick_bin*e12t3) (*, *, index))118 ; e12t3 = reform((reform(e12t[premierx:dernierx, premiery:derniery], nx*ny))#replicate(1,n_sig +2), nx, ny, n_sig+2) 119 ; nfild = ((thick_bin*e12t3)[*, *, index]) 118 120 END 119 ELSE: nfild = f_bin (*, *, index)121 ELSE: nfild = f_bin[*, *, index] 120 122 ENDCASE 121 123 … … 124 126 IF final_dim eq 3 THEN BEGIN 125 127 indx = where(finite(nfild, /NAN)) 126 nfild (indx)= -1.e20127 nfild (indx)= !VALUES.F_NAN128 indx = where(tmask (*, *, 0)EQ 0)129 nfild (indx)= !VALUES.F_NAN128 nfild[indx] = -1.e20 129 nfild[indx] = !VALUES.F_NAN 130 indx = where(tmask[*, *, 0] EQ 0) 131 nfild[indx] = !VALUES.F_NAN 130 132 ENDIF ELSE BEGIN 131 133 -
trunk/tools/density_binning/density_bin_IDL_gm/bin_velocity.pro
r165 r169 9 9 , VV=vv $ 10 10 , TT=tt 11 ; 12 compile_opt idl2, strictarrsubs 11 13 ; 12 14 @common … … 33 35 ; initializations 34 36 s_s = sigma_values ; density range over which we bin 35 N_s = n_elements(s_s) ; number of density levels36 N_z = jpk ; number of z levels37 n_s = n_elements(s_s) ; number of density levels 38 n_z = jpk ; number of z levels 37 39 ; Profiles along z level ( suffixe _z) 38 c1_z = fltarr( N_z) ; profil du contenu vertical de x39 s_z = fltarr( N_z) ; profil de la densite40 c1_z = fltarr(n_z) ; profil du contenu vertical de x 41 s_z = fltarr(n_z) ; profil de la densite 40 42 z_zt = gdept ; profondeur au point T (k=0 -> 5m) 41 43 ; z_zw = gdepw (old) 42 44 z_zw=shift(gdepw,-1) ; W (k=0 -> 10m) 43 z_zw (jpk-1)=gdepw(jpk-1)45 z_zw[jpk-1]=gdepw[jpk-1] 44 46 ; Profiles along density level (suffixe _s) 45 z_s = fltarr( N_s+1)46 c1_s = fltarr( N_s+1)47 x1_s = fltarr( N_s+1)47 z_s = fltarr(n_s+1) 48 c1_s = fltarr(n_s+1) 49 x1_s = fltarr(n_s+1) 48 50 ; x1 binned on density (output array) 49 x1_bin = fltarr(jpi, jpj, N_s+1)51 x1_bin = fltarr(jpi, jpj, n_s+1) 50 52 ; 51 53 ; vertical content of x1 (ensure the integrale conservation) 52 54 x1_content = fltarr(jpi, jpj, jpk) 53 x1_content (*, *, 0) = e3t(0) * x1(*, *, 0) * xmask(*,*,0)54 FOR k = 1, jpk-1 DO x1_content (*, *, k) = x1_content(*, *, k-1)$55 + e3t (k) * x1(*, *, k) * xmask(*,*,k)55 x1_content[*, *, 0] = e3t[0] * x1[*, *, 0] * xmask[*,*,0] 56 FOR k = 1, jpk-1 DO x1_content[*, *, k] = x1_content[*, *, k-1] $ 57 + e3t[k] * x1[*, *, k] * xmask[*,*,k] 56 58 ; 57 59 ; Loop over the horizontal domain 2D … … 62 64 ; Indices des points T dans l ocean 63 65 ; 64 i_ocean = where(xmask (i, j, *)EQ 1)66 i_ocean = where(xmask[i, j, *] EQ 1) 65 67 z_s = z_s*0. 66 68 c1_s = c1_s*0. … … 68 70 IF (i_ocean[0] NE -1) THEN BEGIN ; on n'entre que si il y a des points ocean 69 71 70 s_z (*) = density(i, j, *)71 c1_z (*) = x1_content(i, j, *)72 s_z[*] = density[i, j, *] 73 c1_z[*] = x1_content[i, j, *] 72 74 ; critere supplementaire a imposer sur le profil pour eviter les cas 73 75 ; pathologiques en attendant d'ecrire une vraie commande d'extraction 74 76 ; de progils stritement croissant. Il s'agit donc d'un test adhoc. 75 IF s_z (0) LT s_z(i_ocean(n_elements(i_ocean)-1))THEN BEGIN77 IF s_z[0] LT s_z[i_ocean[n_elements(i_ocean)-1]] THEN BEGIN 76 78 ;------------------------------------------------------------------------ 77 79 ; controle si le profil est bien strictement croissant (pour l'instant 78 80 ; inutilis (avis aux amateurs) 79 81 ; 80 ; ds = (shift(s_z, -1)-s_z) (i_ocean(0:n_elements(i_ocean)-2))82 ; ds = (shift(s_z, -1)-s_z)[i_ocean[0:n_elements(i_ocean)-2]] 81 83 ; ds(0) = 0 82 84 ; ind = where(ds LT 0.) 83 85 ; croissant = ind(0) EQ -1 84 86 ;------------------------------------------------------------------------ 85 i_bottom = i_ocean (n_elements(i_ocean)-1)86 z_s (N_s) = z_zw(i_bottom)87 c1_s (N_s) = x1_content(i, j, jpk-1)87 i_bottom = i_ocean[n_elements(i_ocean)-1] 88 z_s[n_s] = z_zw[i_bottom] 89 c1_s[n_s] = x1_content[i, j, jpk-1] 88 90 89 91 ; extraction d'un sous profil strictement croissant 90 mini = min(s_z (i_ocean))91 maxi = max(s_z (i_ocean))92 i_min = where(s_z (i_ocean)EQ mini)93 i_max = where(s_z (i_ocean)EQ maxi)92 mini = min(s_z[i_ocean]) 93 maxi = max(s_z[i_ocean]) 94 i_min = where(s_z[i_ocean] EQ mini) 95 i_max = where(s_z[i_ocean] EQ maxi) 94 96 ; on prend le plus grand des indices min 95 97 ; et le plus petit des indices max … … 97 99 ;gr i_max = i_max[0] 98 100 i_min = i_min[0] 99 i_max = i_max (n_elements(i_min)-1)101 i_max = i_max[n_elements(i_min)-1] 100 102 ; IF i_max GE jpk-1 THEN print, i, j, i_max 101 103 ; IF i_min LE 1 THEN print, i, j, i_min … … 104 106 ; l isopycne est mise en surface (z_s=0) 105 107 ; 106 ind = where(s_s LT s_z (i_min))108 ind = where(s_s LT s_z[i_min]) 107 109 IF ind[0] NE -1 THEN BEGIN &$ 108 z_s (ind)= 0 &$109 c1_s (ind)= 0 &$110 z_s[ind] = 0 &$ 111 c1_s[ind] = 0 &$ 110 112 ENDIF 111 113 ; Si la valeur du niveau (s_s) est plus elevee que la densite du fond, 112 114 ; l isopycne est mise au fond (z_s=z_zw(i_bottom)) 113 115 ; 114 ind = where(s_s GT s_z (i_max))116 ind = where(s_s GT s_z[i_max]) 115 117 IF ind[0] NE -1 THEN BEGIN &$ 116 z_s (ind) = z_s(N_s)&$117 c1_s (ind) = c1_s(N_s)&$118 z_s[ind] = z_s[n_s] &$ 119 c1_s[ind] = c1_s[n_s] &$ 118 120 ENDIF 119 121 ; 120 ds =(shift(s_z, -1)-s_z) (i_ocean(0:n_elements(i_ocean)-2))122 ds =(shift(s_z, -1)-s_z)[i_ocean[0:n_elements(i_ocean)-2]] 121 123 ind_c = where(ds LT 0.) 122 ; croissant = ind_c (0)EQ -1124 ; croissant = ind_c[0] EQ -1 123 125 ; IF ( i_min GT i_max ) then begin 124 126 ; print, 'bug ',i_min, i_max,' en i,j=', i,j 125 127 ; endif 126 IF ( ind_c (0)NE -1 ) then begin128 IF ( ind_c[0] NE -1 ) then begin 127 129 print, 'bug en i,j=', i,j,' ind_c = ', ind_c 128 130 stop 129 131 endif 130 132 ; IF ( i_min LE i_max ) then begin 131 IF ( ind_c (0)EQ -1 ) then begin132 ind = where( (s_s GE s_z (i_min)) AND (s_s LT s_z(i_max)) )133 IF ( ind_c[0] EQ -1 ) then begin 134 ind = where( (s_s GE s_z[i_min]) AND (s_s LT s_z[i_max]) ) 133 135 IF ind[0] NE -1 THEN BEGIN 134 136 IF ( n_elements(ind) NE 1 ) THEN BEGIN 135 137 136 i_profil = i_ocean (i_min:i_max)137 z_s (ind) = interpol(z_zt(i_profil), s_z(i_profil), s_s(ind))138 i_profil = i_ocean[i_min:i_max] 139 z_s[ind] = interpol(z_zt[i_profil], s_z[i_profil], s_s[ind]) 138 140 139 141 ; 140 142 ; j'utilise la fonction spline pour interpoler le contenu 141 143 ; 142 c1_s (ind) = spline(z_zw(i_profil), c1_z(i_profil), z_s(ind), 1)144 c1_s[ind] = spline(z_zw[i_profil], c1_z[i_profil], z_s[ind], 1) 143 145 ; 144 146 ; l'interpolation lineaire marche aussi. Elle donne des resultats differents. Elle … … 152 154 IF ( n_elements(ind) EQ 1 ) THEN BEGIN 153 155 ; print, 'one density bin en i,j=', i,j,' ind = ', ind 154 c1_s (ind) = c1_z(jpk-1)156 c1_s[ind] = c1_z[jpk-1] 155 157 ENDIF 156 158 ; 157 x1_s (ind[0]-1) = (c1_s(ind[0])-c1_s(ind[0]-1));/(z_s(ind+1)-z_s(ind))158 x1_s (ind) = (c1_s(ind+1)-c1_s(ind));/(z_s(ind+1)-z_s(ind))159 ; x1_s (ind[0]-1) = (c1_s(ind[0])-c1_s(ind[0]-1))/(z_s(ind+1)-z_s(ind))160 ; x1_s (ind) = (c1_s(ind+1)-c1_s(ind))/(z_s(ind+1)-z_s(ind))159 x1_s[ind[0]-1] = (c1_s[ind[0]]-c1_s[ind[0]-1]);/(z_s[ind+1]-z_s[ind]) 160 x1_s[ind] = (c1_s[ind+1]-c1_s[ind]);/(z_s[ind+1]-z_s[ind]) 161 ; x1_s[ind[0]-1] = (c1_s[ind[0]]-c1_s[ind[0]-1])/(z_s[ind+1]-z_s[ind]) 162 ; x1_s[ind] = (c1_s[ind+1]-c1_s[ind])/(z_s[ind+1]-z_s[ind]) 161 163 ENDIF 162 164 endif 163 165 ENDIF 164 166 ENDIF 165 x1_bin (i, j, *)= x1_s167 x1_bin[i, j, *] = x1_s 166 168 167 169 ENDFOR -
trunk/tools/density_binning/density_bin_IDL_gm/bining2.pro
r168 r169 1 1 ;+ 2 ; 2 3 ; PRINCIPE DU BINING 3 4 ; … … 87 88 , TMASK=tmask 88 89 ; 90 compile_opt idl2, strictarrsubs 91 ; 89 92 size3d = size(density) 90 jpi = size3d (1)91 jpj = size3d (2)92 jpk = size3d (3)93 jpi = size3d[1] 94 jpj = size3d[2] 95 jpk = size3d[3] 93 96 ;print, 'size3d=', size3d 94 97 ;print, 'size(tmask)=', size(tmask) … … 96 99 indm = where(tmask eq 0) 97 100 98 x1 (indm)= !VALUES.F_NAN99 density (indm)= !VALUES.F_NAN100 101 N_z = jpk101 x1[indm] = !VALUES.F_NAN 102 density[indm] = !VALUES.F_NAN 103 104 n_z = jpk 102 105 103 106 ; adjust profiles to avoid static instabilities (there are almost none IF using pac_40) 104 107 IF jpi EQ 180 THEN BEGIN 105 density = npc (density)108 density = npc[density] 106 109 ENDIF 107 110 108 ; N_s, defini ci dessous est le nombre d'interface111 ; n_s, defini ci dessous est le nombre d'interface 109 112 110 113 IF keyword_set(sigma) THEN BEGIN 111 114 s_s = sigma 112 N_s = n_elements(s_s)115 n_s = n_elements(s_s) 113 116 ENDIF ELSE BEGIN 114 117 ; … … 116 119 ; les eaux intermediaires) 117 120 ; 118 N_s = 20119 s_s = 27.6+0.2*findgen( N_s)/(N_s-1)120 N_s = 3121 n_s = 20 122 s_s = 27.6+0.2*findgen(n_s)/(n_s-1) 123 n_s = 3 121 124 s_s = [26.8, 27.6, 27.8] 122 125 s_s = [22.8, 27.6, 27.8] … … 129 132 ; Profils selon les niveaux du modele (suffixe _z) 130 133 ; 131 c1_z = fltarr( N_z) ; profil du contenu vertical de x132 s_z = fltarr( N_z) ; profil de la densite134 c1_z = fltarr(n_z) ; profil du contenu vertical de x 135 s_z = fltarr(n_z) ; profil de la densite 133 136 z_zt = depth_t ; profondeur au point T (k=0 -> 5m) 134 137 z_zw = depth_w ; W (k=0 -> 10m) … … 136 139 ; Profils selon les couches de densite (suffixe _s) 137 140 ; 138 z_s = fltarr( N_s+1)139 c1_s = fltarr( N_s+1)140 x1_s = fltarr( N_s+1)141 z_s = fltarr(n_s+1) 142 c1_s = fltarr(n_s+1) 143 x1_s = fltarr(n_s+1) 141 144 bowl_s = 0. 142 145 143 146 ; Tableaux de sorties 144 147 ; 145 depth_bin = fltarr(jpi, jpj, N_s+1)146 thick_bin = fltarr(jpi, jpj, N_s+1)147 x1_bin = fltarr(jpi, jpj, N_s+1)148 depth_bin = fltarr(jpi, jpj, n_s+1) 149 thick_bin = fltarr(jpi, jpj, n_s+1) 150 x1_bin = fltarr(jpi, jpj, n_s+1) 148 151 bowl_bin = fltarr(jpi, jpj) 149 152 150 x1_bin (*, *, *)= !VALUES.F_NAN151 depth_bin (*, *, *)= !VALUES.F_NAN152 thick_bin (*, *, *)= !VALUES.F_NAN153 bowl_bin (*, *)= !VALUES.F_NAN153 x1_bin[*, *, *] = !VALUES.F_NAN 154 depth_bin[*, *, *] = !VALUES.F_NAN 155 thick_bin[*, *, *] = !VALUES.F_NAN 156 bowl_bin[*, *] = !VALUES.F_NAN 154 157 155 158 … … 159 162 160 163 161 ; x1_content (*, *, 0) = e3t(0) * x1(*, *, 0) * tmask (*, *, 0)162 ; FOR k = 1, jpk-1 DO x1_content (*, *, k) = x1_content(*, *, k-1)$163 ; + e3t (k) * x1(*, *, k)* tmask (*, *, k)164 ; x1_content[*, *, 0] = e3t[0] * x1[*, *, 0] * tmask [*, *, 0] 165 ; FOR k = 1, jpk-1 DO x1_content[*, *, k] = x1_content[*, *, k-1] $ 166 ; + e3t[k] * x1[*, *, k]* tmask [*, *, k] 164 167 x1_content = x1 165 168 … … 181 184 ; 182 185 ;stop 183 i_ocean = where(tmask (i, j, *)EQ 1)184 ;print, 'tmask (i, j, *)=', tmask(i, j, *)186 i_ocean = where(tmask[i, j, *] EQ 1) 187 ;print, 'tmask[i, j, *]=', tmask[i, j, *] 185 188 ;print, 'i_ocean=', i_ocean 186 189 z_s = z_s*0. 187 c1_s (*)= !VALUES.F_NAN190 c1_s[*] = !VALUES.F_NAN 188 191 bowl_s = !VALUES.F_NAN 189 ; x1_s (*)= !VALUES.F_NAN192 ; x1_s[*] = !VALUES.F_NAN 190 193 191 194 192 195 IF i_ocean[0] NE -1 THEN BEGIN ; on n entre que si il y a des points ocean 193 196 ; print, 'ocean point' 194 i_bottom = i_ocean (n_elements(i_ocean)-1)197 i_bottom = i_ocean[n_elements(i_ocean)-1] 195 198 ;print, 'i_bottom=', i_bottom 196 z_s (N_s) = z_zw(i_bottom)197 c1_s (N_s) = x1_content(i, j, jpk-1)199 z_s[n_s] = z_zw[i_bottom] 200 c1_s[n_s] = x1_content[i, j, jpk-1] 198 201 199 s_z (*) = density(i, j, *)200 c1_z (*) = x1_content(i, j, *)202 s_z[*] = density[i, j, *] 203 c1_z[*] = x1_content[i, j, *] 201 204 ; print, 'density profile s_z', s_z 202 205 ; print, 'field profile c1_z', c1_z … … 205 208 206 209 ; extraction d'un sous profil strictement croissant 207 ; print, 's_z (i_ocean)=', s_z(i_ocean)208 mini = min(s_z (i_ocean))209 maxi = max(s_z (i_ocean))210 i_min = where(s_z (i_ocean)EQ mini)211 i_max = where(s_z (i_ocean)EQ maxi)210 ; print, 's_z[i_ocean]=', s_z[i_ocean] 211 mini = min(s_z[i_ocean]) 212 maxi = max(s_z[i_ocean]) 213 i_min = where(s_z[i_ocean] EQ mini) 214 i_max = where(s_z[i_ocean] EQ maxi) 212 215 ; print, 'mini, maxi', mini, maxi 213 216 ; print, 'i_min, i_max', i_min, i_max 214 217 ; on prend le plus grand des indices min 215 218 ; et le plus petit des indices max 216 i_min = i_min (n_elements(i_min)-1)219 i_min = i_min[n_elements(i_min)-1] 217 220 i_max = i_max[0] 218 221 ; print, 'i_min, i_max', i_min, i_max … … 226 229 ; l isopycne est mise en surface (z_s=0) 227 230 ; 228 ; print, 's_z (i_min)', s_z(i_min)231 ; print, 's_z[i_min]', s_z[i_min] 229 232 ;print, 's_s=', s_s 230 ind = where(s_s LT s_z (i_min))233 ind = where(s_s LT s_z[i_min]) 231 234 IF ind[0] NE -1 THEN BEGIN 232 235 ; IF i_min GT i_max THEN print, 'min reached at sigma indices', ind 233 z_s (ind)= 0234 c1_s (ind)= !VALUES.F_NAN236 z_s[ind] = 0 237 c1_s[ind] = !VALUES.F_NAN 235 238 ENDIF 236 239 … … 238 241 ; l isopycne est mise au fond (z_s=z_zw(i_bottom)) 239 242 ; 240 ; print, 's_z (i_max)', s_z(i_max)241 ind = where(s_s GT s_z (i_max))243 ; print, 's_z[i_max]', s_z[i_max] 244 ind = where(s_s GT s_z[i_max]) 242 245 243 246 IF ind[0] NE -1 THEN BEGIN 244 247 ; IF i_min GT i_max THEN print, 'max reached at sigma indices', ind 245 z_s (ind) = z_s(N_s)246 c1_s (ind) = c1_s(N_s)247 c1_s (ind)= !VALUES.F_NAN248 z_s[ind] = z_s[n_s] 249 c1_s[ind] = c1_s[n_s] 250 c1_s[ind] = !VALUES.F_NAN 248 251 ENDIF 249 252 ; cas general 250 ind = where( (s_s GE s_z (i_min)) AND (s_s LE s_z(i_max)) )253 ind = where( (s_s GE s_z[i_min]) AND (s_s LE s_z[i_max]) ) 251 254 ; IF i_min GT i_max THEN print, 's_s indices inside min/max', ind 252 255 IF ind[0] NE -1 THEN BEGIN 253 i_profil = i_ocean (i_min:i_max); problem line254 255 z_s (ind) = interpol(z_zt(i_profil), s_z(i_profil), s_s(ind)) ; original256 ; z_s (ind) = interpol(z_zw(i_profil), s_z(i_profil), s_s(ind)) ; changed to z_zw 1/7/04256 i_profil = i_ocean[i_min:i_max] ; problem line 257 258 z_s[ind] = interpol(z_zt[i_profil], s_z[i_profil], s_s[ind]) ; original 259 ; z_s[ind] = interpol(z_zw[i_profil], s_z[i_profil], s_s[ind]) ; changed to z_zw 1/7/04 257 260 258 261 ; … … 262 265 ; 263 266 IF n_elements(i_profil) GT 2 THEN BEGIN 264 ; c1_s (ind) = spline(z_zw(i_profil), c1_z(i_profil), z_s(ind), 1) ; cubic spline265 c1_s (ind) = interpol(c1_z(i_profil), z_zt(i_profil), z_s(ind)) ; linear interpolation266 ; SPLINE_P, z_zw (i_profil), c1_z(i_profil), z_s(ind), c1_s(ind) ; generalization of spline(*)267 ;print, z_zw (i_profil), c1_z(i_profil), z_s(ind), c1_s(ind)267 ; c1_s[ind] = spline(z_zw[i_profil], c1_z[i_profil], z_s[ind], 1) ; cubic spline 268 c1_s[ind] = interpol(c1_z[i_profil], z_zt[i_profil], z_s[ind]) ; linear interpolation 269 ; SPLINE_P, z_zw[i_profil], c1_z[i_profil], z_s[ind], c1_s[ind] ; generalization of spline[*] 270 ;print, z_zw[i_profil], c1_z[i_profil], z_s[ind], c1_s[ind] 268 271 ENDIF ELSE BEGIN 269 c1_s (ind) = interpol(c1_z(i_profil), z_zt(i_profil), z_s(ind))272 c1_s[ind] = interpol(c1_z[i_profil], z_zt[i_profil], z_s[ind]) 270 273 ENDELSE 271 274 ; … … 273 276 274 277 IF sig_bowl EQ 1 THEN BEGIN 275 bowl_s = interpol(s_z (i_profil), z_zt(i_profil), sobwlmax(i, j))278 bowl_s = interpol(s_z[i_profil], z_zt[i_profil], sobwlmax[i, j]) 276 279 ENDIF ELSE bowl_s = 0 277 280 … … 279 282 ; 280 283 ; 281 ; x1_s (ind) = (c1_s(ind+1)-c1_s(ind))/(z_s(ind+1)-z_s(ind))284 ; x1_s[ind] = (c1_s[ind+1]-c1_s[ind])/(z_s[ind+1]-z_s[ind]) 282 285 ; x1_s = c1_s 283 286 … … 285 288 ;ELSE print, 'ind[0] = -1', ind, i_bottom 286 289 ENDIF 287 ;ELSE print, ' land ', tmask (i, j, 1)288 depth_bin (i, j, *)= z_s289 thick_bin (i, j, 0) = z_s(0)290 thick_bin (i, j, 1:N_s) = z_s(1:N_s)-z_s(0:N_s-1)291 IF total(thick_bin (i, j, 1:N_s)LT 0) GT 0 THEN BEGIN292 ;print, 'WARNING: negative thick_bin (i, j, 1:N_s)at i,j= ', i, j293 ;print, 'depth_bin (i, j, 1:N_s)= ', depth_bin (i, j, 1:N_s)294 ;print, 'thick_bin (i, j, 1:N_s)= ', thick_bin (i, j, 1:N_s)290 ;ELSE print, ' land ', tmask[i, j, 1] 291 depth_bin[i, j, *] = z_s 292 thick_bin[i, j, 0] = z_s[0] 293 thick_bin[i, j, 1:n_s] = z_s[1:n_s]-z_s[0:n_s-1] 294 IF total(thick_bin[i, j, 1:n_s] LT 0) GT 0 THEN BEGIN 295 ;print, 'WARNING: negative thick_bin[i, j, 1:n_s] at i,j= ', i, j 296 ;print, 'depth_bin[i, j, 1:n_s]= ', depth_bin[i, j, 1:n_s] 297 ;print, 'thick_bin[i, j, 1:n_s]= ', thick_bin[i, j, 1:n_s] 295 298 ;print, 's_z= ', s_z 296 299 ;stop 297 300 ENDIF 298 x1_bin (i, j, *)= c1_s299 bowl_bin (i, j)= bowl_s301 x1_bin[i, j, *] = c1_s 302 bowl_bin[i, j] = bowl_s 300 303 ; print, 'c1_s', c1_s 301 304 ; … … 303 306 ENDFOR 304 307 ; mask depth_bin with undefined values of x_bin 305 depth_bin (where(finite(x1_bin, /nan) EQ 1))= !VALUES.F_NAN306 thick_bin (where(finite(x1_bin, /nan) EQ 1))= !VALUES.F_NAN308 depth_bin[where(finite(x1_bin, /nan) EQ 1)] = !VALUES.F_NAN 309 thick_bin[where(finite(x1_bin, /nan) EQ 1)] = !VALUES.F_NAN 307 310 308 311 ; OPTIONAL: compute spiciness by removing isopycnal mean for each sigma (T or S) … … 310 313 print, ' Computing spiciness...' 311 314 FOR i = 0, (size(x1_bin))[3]-1 DO BEGIN 312 x1_bin (*,*,i) = x1_bin(*,*,i) - mean(x1_bin(*,*,i),/NAN)315 x1_bin[*,*,i] = x1_bin[*,*,i] - mean(x1_bin[*,*,i],/NAN) 313 316 ENDFOR 314 317 ENDIF … … 317 320 IF 0 THEN BEGIN 318 321 print, ' Removing domain mean...' 319 x1_bin (*,*,*) = x1_bin(*,*,*) - mean(x1_bin(*,*,*),/NAN)322 x1_bin[*,*,*] = x1_bin[*,*,*] - mean(x1_bin[*,*,*],/NAN) 320 323 ENDIF 321 324 … … 323 326 IF 1 THEN BEGIN 324 327 print, ' Removing 34.6psu...' 325 x1_bin (*,*,*) = x1_bin(*,*,*)- 34.6328 x1_bin[*,*,*] = x1_bin[*,*,*] - 34.6 326 329 ENDIF 327 330 -
trunk/tools/density_binning/density_bin_IDL_gm/eos.pro
r168 r169 52 52 ;- 53 53 FUNCTION eos, t, s, sigma 54 ; 55 compile_opt idl2, strictarrsubs 54 56 ; 55 57 @common -
trunk/tools/density_binning/density_bin_IDL_gm/msf.pro
r165 r169 14 14 PRO msf, v, sf $ 15 15 , MSK=msk 16 ; 17 compile_opt idl2, strictarrsubs 16 18 ; 17 19 @common … … 34 36 IF n_elements(msk) NE 0 THEN BEGIN 35 37 zvmask = boundperio( (msk +shift(msk, 0, -1) ) < 1 ) 36 zvmask = zvmask (*)#vert38 zvmask = zvmask[*]#vert 37 39 zv = zv*zvmask 38 40 ENDIF … … 44 46 ; calcul du flux 45 47 ; 46 FOR i = 0, jpi-1 DO BEGIN for j= 0, jpj-1 do begin ze1v (i,j,*) = replicate(e1v(i,j),jpk) & endfor &endfor47 FOR k = 0, jpk-1 DO BEGIN ze3v (*,*,k)=replicate(e3t(k),jpi*jpj) & endfor48 FOR i = 0, jpi-1 DO BEGIN for j= 0, jpj-1 do begin ze1v[i,j,*] = replicate(e1v[i,j],jpk) & endfor &endfor 49 FOR k = 0, jpk-1 DO BEGIN ze3v[*,*,k]=replicate(e3t[k],jpi*jpj) & endfor 48 50 ; 49 51 z= -v*ze1v*ze3v … … 55 57 ; calcul de la msf en integrant depuis le fond 56 58 ; 57 FOR k = jpk-2, 0, -1 DO begin sf (*, k) = sf(*, k+1)+fm(*, k)& endfor59 FOR k = jpk-2, 0, -1 DO begin sf[*, k] = sf[*, k+1]+fm[*, k] & endfor 58 60 ; 59 61 ; msfmsk est le masque associe a msf (utilise pour les graphiques) -
trunk/tools/density_binning/density_bin_IDL_gm/msfs.pro
r165 r169 8 8 , MASK=mask $ 9 9 , ORCA=orca 10 ; 11 compile_opt idl2, strictarrsubs 10 12 ; 11 13 @common -
trunk/tools/density_binning/density_bin_IDL_gm/msfz.pro
r165 r169 9 9 , ORCA=orca $ 10 10 , _EXTRA=extra 11 ; 12 compile_opt idl2, strictarrsubs 11 13 ; 12 14 @common … … 40 42 for k=0, jpk-1 do begin 41 43 for j=0,jpj-1 do begin 42 msf (j,k) = msf(j,k) - msf(j,jpk-1)44 msf[j,k] = msf[j,k] - msf[j,jpk-1] 43 45 endfor 44 46 endfor -
trunk/tools/density_binning/density_bin_IDL_gm/npc.pro
r168 r169 29 29 ;- 30 30 FUNCTION npc, s3d 31 ; 32 compile_opt idl2, strictarrsubs 31 33 ; 32 34 @common … … 63 65 ; ----------------------------- 64 66 ; 65 ds =(shift(s_z,-1)-s_z) (i_ocean(0:n_elements(i_ocean)-2))67 ds =(shift(s_z,-1)-s_z)[i_ocean[0:n_elements(i_ocean)-2]] 66 68 ind_c = where(ds LT 0.) 67 69 ; 68 70 ; 2. Vertical mixing for each instable portion of the density profil 69 71 ; 70 IF ( ind_c (0)NE -1 ) THEN BEGIN72 IF ( ind_c[0] NE -1 ) THEN BEGIN 71 73 ncompt=ncompt+1 72 74 ; print, 'static instability at i,j=', i,j … … 77 79 ; vertical iteration 78 80 jiter = 0 79 WHILE ( (ind_c (0)NE -1) AND (jiter LT jpk-1) ) DO BEGIN &$81 WHILE ( (ind_c[0] NE -1) AND (jiter LT jpk-1) ) DO BEGIN &$ 80 82 jiter = jiter+1 &$ 81 83 ; ikup : the first static instability from the sea surface 82 ikup = ind_c (0)&$84 ikup = ind_c[0] &$ 83 85 ; the density profil is instable below ikup 84 86 ; ikdown : bottom of the instable portion of the density profil 85 87 ; search of ikdown and vertical mixing from ikup to ikdown 86 ze3tot= e3t (ikup)&$87 zraua = s_z (ikup)&$88 ze3tot= e3t[ikup] &$ 89 zraua = s_z[ikup] &$ 88 90 jkdown = ikup+1 &$ 89 91 ; 90 WHILE (jkdown LE ikbot AND zraua GT s_z (jkdown)) DO BEGIN &$91 ze3dwn = e3t (jkdown)&$92 WHILE (jkdown LE ikbot AND zraua GT s_z[jkdown] ) DO BEGIN &$ 93 ze3dwn = e3t[jkdown] &$ 92 94 ze3tot = ze3tot+ze3dwn &$ 93 zraua = ( zraua*(ze3tot-ze3dwn) + s_z (jkdown)*ze3dwn )/ze3tot &$95 zraua = ( zraua*(ze3tot-ze3dwn) + s_z[jkdown]*ze3dwn )/ze3tot &$ 94 96 jkdown=jkdown+1 &$ 95 97 ; print, jkdown, zraua &$ … … 99 101 s_z(jkp) = zraua &$ 100 102 ENDFOR &$ 101 ds =(shift(s_z, -1)-s_z) (i_ocean(0:n_elements(i_ocean)-2))&$103 ds =(shift(s_z, -1)-s_z)[i_ocean[0:n_elements(i_ocean)-2]] &$ 102 104 ind_c = where(ds LT 0.) &$ 103 105 ENDWHILE 104 106 ENDIF 105 107 ; save the modifications 106 rhos (i,j,*) = s_z(*)108 rhos[i,j,*] = s_z[*] 107 109 ; 108 110 ; <<-- no more static instability on slab jj -
trunk/tools/idl_netcdf/ncdf_ajoute.pro
r165 r169 35 35 ;- 36 36 PRO ncdf_ajoute, nomfich, nom, tab, dimension, attribut 37 ; 38 compile_opt idl2, strictarrsubs 37 39 ; 38 40 ;----------------------------------------------------------------------------------------- -
trunk/tools/idl_netcdf/ncdf_ajoute2.pro
r165 r169 34 34 ; 35 35 ;- 36 PRO ncdf_ajoute, nomfich, nom, tab, dimension, attribut 36 PRO ncdf_ajoute2, nomfich, nom, tab, dimension, attribut 37 ; 38 compile_opt idl2, strictarrsubs 37 39 ; 38 40 ;----------------------------------------------------------------------------------------- -
trunk/tools/idl_netcdf/ncdf_colle.pro
r165 r169 75 75 , GARDE=garde $ 76 76 , EXCLU=exclu 77 ; 78 compile_opt idl2, strictarrsubs 77 79 ; 78 80 nomdirec = strlowcase(nomdirec) -
trunk/tools/idl_netcdf/ncdf_colle2.pro
r165 r169 72 72 ; 73 73 ;- 74 PRO ncdf_colle , listin, nomfichout, nomdirec $74 PRO ncdf_colle2, listin, nomfichout, nomdirec $ 75 75 , GARDE=garde $ 76 76 , EXCLU=exclu 77 ; 78 compile_opt idl2, strictarrsubs 77 79 ; 78 80 nomdirec = strlowcase(nomdirec) -
trunk/tools/idl_netcdf/ncdf_extrait.pro
r165 r169 56 56 , GARDE=garde $ 57 57 , EXCLU=exclu 58 ; 59 compile_opt idl2, strictarrsubs 58 60 ; 59 61 ;------------------------------------------------------------ -
trunk/tools/lma/post_lma_patt_am.pro
r166 r169 6 6 ;- 7 7 PRO post_lma_patt_am 8 ; 9 compile_opt idl2, strictarrsubs 8 10 ; 9 11 ; 1.Declaration des variables … … 31 33 FOR i = 0,xdim-1 DO BEGIN 32 34 READF, 10, FORMAT = '(2E12.5)',re, im 33 patt_sum1r (i, j)= re34 patt_sum1i (i, j)= im35 patt_sum1r[i, j] = re 36 patt_sum1i[i, j] = im 35 37 ENDFOR 36 38 ENDFOR … … 39 41 FOR i = 0,xdim-1 DO BEGIN 40 42 READF, 10, FORMAT = '(2E12.5)',re, im 41 patt_sum2r (i, j)= re42 patt_sum2i (i, j)= im43 patt_sum2r[i, j] = re 44 patt_sum2i[i, j] = im 43 45 ENDFOR 44 46 ENDFOR … … 47 49 FOR i = 0,xdim-1 DO BEGIN 48 50 READF, 10, FORMAT = '(2E12.5)',re, im 49 patt_win1r (i, j)= re50 patt_win1i (i, j)= im51 patt_win1r[i, j] = re 52 patt_win1i[i, j] = im 51 53 ENDFOR 52 54 ENDFOR … … 55 57 FOR i = 0,xdim-1 DO BEGIN 56 58 READF, 10, FORMAT = '(2E12.5)',re, im 57 patt_win2r (i, j)= re58 patt_win2i (i, j)= im59 patt_win2r[i, j] = re 60 patt_win2i[i, j] = im 59 61 ENDFOR 60 62 ENDFOR … … 70 72 FOR t = 0, NCMAX-1 DO BEGIN 71 73 READF, 20, FORMAT = '(2E12.5,12I5)',res1, res2, res3, res4, res5, res6, res7, res8, res9 72 selec (t)= res874 selec[t] = res8 73 75 ENDFOR 74 76 close, 20 … … 79 81 FOR i = 0,xdim-1 DO BEGIN 80 82 READF, 30, FORMAT = '(3E12.5)', re, im, var 81 vari (i, j, ta)= var83 vari[i, j, ta] = var 82 84 ENDFOR 83 85 ENDFOR … … 348 350 return 349 351 end 350 -
trunk/tools/lma/postlma_distpm_am.pro
r166 r169 9 9 PRO postlma_distpm_am $ 10 10 , POST=post 11 ; 12 compile_opt idl2, strictarrsubs 11 13 ; 12 14 @common … … 41 43 FOR t = 0, NCMAX-1 DO BEGIN 42 44 READF, 10, FORMAT = '(2E12.5,12I5)',res1, res2, res3, res4, res5, res6, res7, res8, res9 43 xver (t)= res144 dto (t)= res245 an (t)= res346 mois (t)= res447 jour (t)= res548 selec (t)= res849 time (t)= julday(res4, res5, res3, 0, 0, 0)45 xver[t] = res1 46 dto[t] = res2 47 an[t] = res3 48 mois[t] = res4 49 jour[t] = res5 50 selec[t] = res8 51 time[t] = julday(res4, res5, res3, 0, 0, 0) 50 52 ENDFOR 51 53 close, 10 … … 58 60 FOR t = 0, NCMAX-1 DO BEGIN 59 61 READF, 20, FORMAT = '(2E12.5)', res 60 DPM_sum1 (t)= res62 DPM_sum1[t] = res 61 63 ENDFOR 62 64 READF, 20, FORMAT = '(A6)', u 63 65 FOR t = 0, NCMAX-1 DO BEGIN 64 66 READF, 20, FORMAT = '(2E12.5)', res 65 DPM_sum2 (t)= res67 DPM_sum2[t] = res 66 68 ENDFOR 67 69 READF, 20, FORMAT = '(A6)', u 68 70 FOR t = 0, NCMAX-1 DO BEGIN 69 71 READF, 20, FORMAT = '(2E12.5)', res 70 DPM_win1 (t)= res72 DPM_win1[t] = res 71 73 ENDFOR 72 74 READF, 20, FORMAT = '(A6)', u 73 75 FOR t = 0, NCMAX-1 DO BEGIN 74 76 READF, 20, FORMAT = '(2E12.5)', res 75 DPM_win2 (t)= res77 DPM_win2[t] = res 76 78 ENDFOR 77 79 READF, 20, u … … 91 93 print, res13 92 94 print, res4, res5, res6 93 pourvar (t)= res1395 pourvar[t] = res13 94 96 ENDFOR 95 97 close, 30 … … 117 119 FOR i = 0, NCMAX-1 DO BEGIN 118 120 g = 0 119 IF (pourvar (i) GT 0.4) and (pourvar(i)LE 0.5) THEN g = 3120 IF (pourvar (i) GT 0.5) and (pourvar(i)LE 0.6) THEN g = 6121 IF (pourvar (i) GT 0.6) and (pourvar(i)LE 0.7) THEN g = 9122 IF (pourvar (i) GT 0.7) and (pourvar(i)LE 0.8) THEN g = 12123 IF (pourvar (i) GT 0.8) and (pourvar(i)LE 0.9) THEN g = 15124 IF selec (i) EQ summer THEN xyouts, time(i),DPM_sum1(i), '.',charthick=5, charsize = g, alignment = align121 IF (pourvar[i] GT 0.4) and (pourvar[i] LE 0.5) THEN g = 3 122 IF (pourvar[i] GT 0.5) and (pourvar[i] LE 0.6) THEN g = 6 123 IF (pourvar[i] GT 0.6) and (pourvar[i] LE 0.7) THEN g = 9 124 IF (pourvar[i] GT 0.7) and (pourvar[i] LE 0.8) THEN g = 12 125 IF (pourvar[i] GT 0.8) and (pourvar[i] LE 0.9) THEN g = 15 126 IF selec[i] EQ summer THEN xyouts, time[i],DPM_sum1[i], '.',charthick=5, charsize = g, alignment = align 125 127 ENDFOR 126 128 … … 129 131 FOR i = 0, NCMAX-1 DO BEGIN 130 132 g = 0 131 IF (pourvar (i) GT 0.4) and (pourvar(i)LE 0.5) THEN g = 3132 IF (pourvar (i) GT 0.5) and (pourvar(i)LE 0.6) THEN g = 6133 IF (pourvar (i) GT 0.6) and (pourvar(i)LE 0.7) THEN g = 9134 IF (pourvar (i) GT 0.7) and (pourvar(i)LE 0.8) THEN g = 12135 IF (pourvar (i) GT 0.8) and (pourvar(i)LE 0.9) THEN g = 15136 IF selec (i) EQ summer THEN xyouts, time(i),DPM_sum2(i), '.',charthick=5, charsize = g, alignment = align, color = 30133 IF (pourvar[i] GT 0.4) and (pourvar[i] LE 0.5) THEN g = 3 134 IF (pourvar[i] GT 0.5) and (pourvar[i] LE 0.6) THEN g = 6 135 IF (pourvar[i] GT 0.6) and (pourvar[i] LE 0.7) THEN g = 9 136 IF (pourvar[i] GT 0.7) and (pourvar[i] LE 0.8) THEN g = 12 137 IF (pourvar[i] GT 0.8) and (pourvar[i] LE 0.9) THEN g = 15 138 IF selec[i] EQ summer THEN xyouts, time[i],DPM_sum2[i], '.',charthick=5, charsize = g, alignment = align, color = 30 137 139 ENDFOR 138 140 … … 143 145 FOR i = 0, NCMAX-1 DO BEGIN 144 146 g = 0 145 IF (pourvar (i) GT 0.4) and (pourvar(i)LE 0.5) THEN g = 3146 IF (pourvar (i) GT 0.5) and (pourvar(i)LE 0.6) THEN g = 6147 IF (pourvar (i) GT 0.6) and (pourvar(i)LE 0.7) THEN g = 9148 IF (pourvar (i) GT 0.7) and (pourvar(i)LE 0.8) THEN g = 12149 IF (pourvar (i) GT 0.8) and (pourvar(i)LE 0.9) THEN g = 15150 IF selec (i) EQ win THEN xyouts, time(i),DPM_win1(i), '.',charthick=5, charsize = g, alignment = align147 IF (pourvar[i] GT 0.4) and (pourvar[i] LE 0.5) THEN g = 3 148 IF (pourvar[i] GT 0.5) and (pourvar[i] LE 0.6) THEN g = 6 149 IF (pourvar[i] GT 0.6) and (pourvar[i] LE 0.7) THEN g = 9 150 IF (pourvar[i] GT 0.7) and (pourvar[i] LE 0.8) THEN g = 12 151 IF (pourvar[i] GT 0.8) and (pourvar[i] LE 0.9) THEN g = 15 152 IF selec[i] EQ win THEN xyouts, time[i],DPM_win1[i], '.',charthick=5, charsize = g, alignment = align 151 153 ENDFOR 152 154 … … 154 156 FOR i = 0, NCMAX-1 DO BEGIN 155 157 g = 0 156 IF (pourvar (i) GT 0.4) and (pourvar(i)LE 0.5) THEN g = 3157 IF (pourvar (i) GT 0.5) and (pourvar(i)LE 0.6) THEN g = 6158 IF (pourvar (i) GT 0.6) and (pourvar(i)LE 0.7) THEN g = 9159 IF (pourvar (i) GT 0.7) and (pourvar(i)LE 0.8) THEN g = 12160 IF (pourvar (i) GT 0.8) and (pourvar(i)LE 0.9) THEN g = 15161 IF selec (i) EQ win THEN xyouts, time(i),DPM_win2(i), '.',charthick=5, charsize = g, alignment = align, color = 30158 IF (pourvar[i] GT 0.4) and (pourvar[i] LE 0.5) THEN g = 3 159 IF (pourvar[i] GT 0.5) and (pourvar[i] LE 0.6) THEN g = 6 160 IF (pourvar[i] GT 0.6) and (pourvar[i] LE 0.7) THEN g = 9 161 IF (pourvar[i] GT 0.7) and (pourvar[i] LE 0.8) THEN g = 12 162 IF (pourvar[i] GT 0.8) and (pourvar[i] LE 0.9) THEN g = 15 163 IF selec[i] EQ win THEN xyouts, time[i],DPM_win2[i], '.',charthick=5, charsize = g, alignment = align, color = 30 162 164 ENDFOR 163 165 if keyword_set(post) THEN BEGIN -
trunk/tools/lma/postlma_var.pro
r166 r169 9 9 PRO postlma_var $ 10 10 , POST=post 11 ; 12 compile_opt idl2, strictarrsubs 11 13 ; 12 14 @common … … 42 44 FOR t = 0, NCMAX-1 DO BEGIN 43 45 READF, 10, FORMAT = '(2E12.5,12I5)',res1, res2, res3, res4, res5, res6, res7, res8, res9 44 xver (t)= res145 dto (t)= res246 an (t)= res347 mois (t)= res448 jour (t)= res549 selec (t)= res850 time (t)= julday(res4, res5, res3, 0, 0, 0)46 xver[t] = res1 47 dto[t] = res2 48 an[t] = res3 49 mois[t] = res4 50 jour[t] = res5 51 selec[t] = res8 52 time[t] = julday(res4, res5, res3, 0, 0, 0) 51 53 ENDFOR 52 54 close, 10 53 55 54 print, 'Initial date dd/mm/yy = ', jour (0), mois(0), an(0)55 print, 'Final date dd/mm/yy = ', jour (NCMAX-1), mois(NCMAX-1), an(NCMAX-1)56 print, 'Initial date dd/mm/yy = ', jour[0], mois[0], an[0] 57 print, 'Final date dd/mm/yy = ', jour[NCMAX-1], mois[NCMAX-1], an[NCMAX-1] 56 58 57 inidat = (an (0)-1)*12+mois(0)58 findat = (an (NCMAX-1)-1)*12+mois(NCMAX-1)59 inidat = (an[0]-1)*12+mois[0] 60 findat = (an[NCMAX-1]-1)*12+mois[NCMAX-1] 59 61 60 62 ; read % var, variance, period, … … 68 70 FOR t = 0, NCMAX-1 DO BEGIN 69 71 READF, 20, res1, res2, res3, res4, res5, res6, res7, res8, res9, res10, res11, res12, res13 70 DPM_pouvar (t)= res971 DPM_varian (t)= sqrt(res12)*res1372 DPM_period (t)= res1072 DPM_pouvar[t] = res9 73 DPM_varian[t] = sqrt(res12)*res13 74 DPM_period[t] = res10 73 75 ENDFOR 74 76 … … 120 122 121 123 IF res1 GE inidat AND res1 LE findat THEN BEGIN 122 nino3_ssta (idx)= res2123 year (idx)= an_t124 month (idx)= mon_t125 time (idx) = julday(month(idx), 15, year(idx), 0, 0, 0)124 nino3_ssta[idx] = res2 125 year[idx] = an_t 126 month[idx] = mon_t 127 time[idx] = julday(month[idx], 15, year[idx], 0, 0, 0) 126 128 idx = idx + 1 127 129 ENDIF … … 134 136 ; take from inidat to findat 135 137 136 print, 'Initial date dd/mm/yy = ', 15, month (0), year(0)137 print, 'Final date dd/mm/yy = ', 15, month (idx), year(idx)138 print, 'Initial date dd/mm/yy = ', 15, month[0], year[0] 139 print, 'Final date dd/mm/yy = ', 15, month[idx], year[idx] 138 140 nino3_ssta = nino3_ssta[0:idx] 139 141 time = time[0:idx] -
trunk/tools/lma/slec.pro
r168 r169 10 10 , _EXTRA=extra 11 11 ; 12 compile_opt idl2, strictarrsubs 13 ; 12 14 @common 13 15 if strlowcase(name) EQ 'un' then name = 'AM_win1i' -
trunk/tools/merid_trans.pro
r165 r169 16 16 ; 17 17 ;- 18 PRO merid $18 PRO merid_trans $ 19 19 , SALT=salt $ 20 20 , PS=ps 21 ; 22 compile_opt idl2, strictarrsubs 21 23 ; 22 24 @common … … 141 143 142 144 ; 2nd window - Atlantic 143 plot,latitude (index),(field7+field8)(index)*ymult,/noerase, linestyle = 0, thick = 2144 plot,latitude (index),field5(index)*ymult,/noerase , linestyle = 1, thick = 2145 plot,latitude (index),field6(index)*ymult,/noerase, linestyle = 2 , thick = 2146 plot,latitude (index),field7(index)*ymult,/noerase, linestyle = 3 , thick = 2145 plot,latitude[index],(field7+field8)[index]*ymult,/noerase, linestyle = 0, thick = 2 146 plot,latitude[index],field5[index]*ymult,/noerase , linestyle = 1, thick = 2 147 plot,latitude[index],field6[index]*ymult,/noerase, linestyle = 2 , thick = 2 148 plot,latitude[index],field7[index]*ymult,/noerase, linestyle = 3 , thick = 2 147 149 plot,!x.crange, [0., 0.], /noerase , thick = 0.5 148 150 plot,[0., 0.],!y.crange, /noerase , thick = 0.5 , title = 'atlantic' … … 154 156 indo3 = (field2-field6) 155 157 indo4 = (field3-field7) 156 plot,latitude (index),indo1(index)*ymult,/noerase, linestyle = 0, thick = 2157 plot,latitude (index),indo2(index)*ymult,/noerase , linestyle = 1, thick = 2158 plot,latitude (index),indo3(index)*ymult,/noerase, linestyle = 2 , thick = 2159 plot,latitude (index),indo4(index)*ymult,/noerase, linestyle = 3 , thick = 2158 plot,latitude[index],indo1[index]*ymult,/noerase, linestyle = 0, thick = 2 159 plot,latitude[index],indo2[index]*ymult,/noerase , linestyle = 1, thick = 2 160 plot,latitude[index],indo3[index]*ymult,/noerase, linestyle = 2 , thick = 2 161 plot,latitude[index],indo4[index]*ymult,/noerase, linestyle = 3 , thick = 2 160 162 plot,!x.crange, [0., 0.], /noerase , thick = 0.5 161 163 plot,[0., 0.],!y.crange, /noerase , thick = 0.5 , title = 'indopac' -
trunk/tools/msf_bin/boundperio.pro
r165 r169 11 11 , VV=vv $ 12 12 , FF=ff 13 ; 14 compile_opt idl2, strictarrsubs 13 15 ; 14 16 @common -
trunk/tools/msf_bin/iniflx.pro
r165 r169 18 18 ;- 19 19 PRO iniflx 20 ; 21 compile_opt idl2, strictarrsubs 20 22 ; 21 23 @common … … 45 47 46 48 FOR l = 0, jpl-1 DO BEGIN 47 lat = gphiv (90, l)49 lat = gphiv[90, l] 48 50 uvmsk = fltarr(jpi, jpj) 49 51 ; uvmsk = uvmsk*0. 50 uvmsk (where(gphit GT lat))= 151 zuleng = boundcut( e2u*(shift(uvmsk, -1, 0)-uvmsk)*(umask) (*, *, 0))52 uvmsk[where(gphit GT lat)] = 1 53 zuleng = boundcut( e2u*(shift(uvmsk, -1, 0)-uvmsk)*(umask)[*, *, 0] ) 52 54 ; ^ 53 55 ; il manquait le signe moins ci dessous/dessus, corrige le 10/04/99 ... 54 56 ; v 55 zvleng = boundcut( e1v*(shift(uvmsk, 0, -1)-uvmsk)*(vmask) (*, *, 0))57 zvleng = boundcut( e1v*(shift(uvmsk, 0, -1)-uvmsk)*(vmask)[*, *, 0] ) 56 58 indu = where( (zuleng NE 0) ) 57 59 IF indu[0] EQ -1 THEN indu = [0] … … 60 62 Nu = n_elements(indu) 61 63 Nv = n_elements(indv) 62 indxu (0:Nu-1, l)= indu63 indxv (0:Nv-1, l)= indv64 zxuleng (0:Nu-1, l) = zuleng(indu)65 zxvleng (0:Nv-1, l) = zvleng(indv)64 indxu [0:Nu-1, l] = indu 65 indxv [0:Nv-1, l] = indv 66 zxuleng[0:Nu-1, l] = zuleng[indu] 67 zxvleng[0:Nv-1, l] = zvleng[indv] 66 68 ENDFOR 67 69 … … 73 75 74 76 FOR l = 0, jpi-3 DO BEGIN 75 lon = glamu (l, 0)77 lon = glamu[l, 0] 76 78 uvmsk = fltarr(jpi, jpj) 77 79 uvmsk = uvmsk*0. 78 uvmsk (where(glamt GT lon))= 179 zuleng = boundcut( e2u*(shift(uvmsk, 1, 0)-uvmsk)*(umask) (*, *, 0))80 zvleng = boundcut( e1v*(shift(uvmsk, 0, 1)-uvmsk)*(vmask) (*, *, 0))80 uvmsk[where(glamt GT lon)] = 1 81 zuleng = boundcut( e2u*(shift(uvmsk, 1, 0)-uvmsk)*(umask)[*, *, 0] ) 82 zvleng = boundcut( e1v*(shift(uvmsk, 0, 1)-uvmsk)*(vmask)[*, *, 0] ) 81 83 indu = where( (zuleng NE 0) ) 82 84 IF indu[0] EQ -1 THEN indu = [0] … … 85 87 Nu = n_elements(indu) 86 88 Nv = n_elements(indv) 87 indyu (l+1, 0:Nu-1)= indu88 indyv (l+1, 0:Nv-1)= indv89 zyuleng (l+1, 0:Nu-1) = zuleng(indu)90 zyvleng (l+1, 0:Nv-1) = zvleng(indv)89 indyu [l+1, 0:Nu-1] = indu 90 indyv [l+1, 0:Nv-1] = indv 91 zyuleng[l+1, 0:Nu-1] = zuleng[indu] 92 zyvleng[l+1, 0:Nv-1] = zvleng[indv] 91 93 ENDFOR 92 END 94 END -
trunk/tools/msf_bin/msf_sigma.pro
r165 r169 29 29 , MSK=msk 30 30 ; 31 compile_opt idl2, strictarrsubs 32 ; 31 33 print, 'Bining U ...' 32 34 u_s = bin_velocity(density, sigma_values, u, /uu) … … 47 49 , UU = uu $ 48 50 , VV = vv 51 ; 52 compile_opt idl2, strictarrsubs 49 53 ; 50 54 @common … … 71 75 c1_z = fltarr(N_z) ; profil du contenu vertical de x 72 76 s_z = fltarr(N_z) ; profil de la densite 73 z_zt = zdept (1:jpk); profondeur au point T (k=0 -> 5m)74 z_zw = zdepw (1:jpk); W (k=0 -> 10m)77 z_zt = zdept[1:jpk] ; profondeur au point T (k=0 -> 5m) 78 z_zw = zdepw[1:jpk] ; W (k=0 -> 10m) 75 79 76 80 ; Profils selon les couches de densite (suffixe _s) … … 88 92 ; 89 93 x1_content = fltarr(jpi, jpj, jpk) 90 x1_content (*, *, 0) = e3t(0) * x1(*, *, 0)91 FOR k = 1, jpk-1 DO x1_content (*, *, k) = x1_content(*, *, k-1)$92 + e3t (k) * x1(*, *, k)94 x1_content[*, *, 0] = e3t[0] * x1[*, *, 0] 95 FOR k = 1, jpk-1 DO x1_content[*, *, k] = x1_content[*, *, k-1] $ 96 + e3t[k] * x1[*, *, k] 93 97 94 98 ; … … 100 104 ; Indices des points T dans l ocean 101 105 ; 102 i_ocean = where(xmask (i, j, *)EQ 1)106 i_ocean = where(xmask[i, j, *] EQ 1) 103 107 104 108 z_s = z_s*0. … … 108 112 IF (i_ocean[0] NE -1) THEN BEGIN ; on n'entre que si il y a des points ocean 109 113 110 s_z (*) = density(i, j, *)111 c1_z (*) = x1_content(i, j, *)114 s_z[*] = density[i, j, *] 115 c1_z[*] = x1_content[i, j, *] 112 116 113 117 ; critere supplementaire a imposer sur le profil pour eviter les cas 114 118 ; pathologiques en attendant d'ecrire une vraie commande d'extraction 115 119 ; de progils stritement croissant. Il s'agit donc d'un test adhoc. 116 IF s_z (0) LT s_z(i_ocean(n_elements(i_ocean)-1))THEN BEGIN120 IF s_z[0] LT s_z[i_ocean[n_elements(i_ocean)-1]] THEN BEGIN 117 121 118 122 ;------------------------------------------------------------------------ … … 120 124 ; inutilisé (avis aux amateurs) 121 125 ; 122 ; ds = (shift(s_z, -1)-s_z) (i_ocean(0:n_elements(i_ocean)-2))123 ; ds (0)= 0126 ; ds = (shift(s_z, -1)-s_z)[i_ocean[0:n_elements(i_ocean)-2]] 127 ; ds[0] = 0 124 128 ; ind = where(ds LT 0.) 125 ; croissant = ind (0)EQ -1129 ; croissant = ind[0] EQ -1 126 130 ;------------------------------------------------------------------------ 127 131 128 i_bottom = i_ocean (n_elements(i_ocean)-1)129 130 z_s (N_s) = z_zw(i_bottom)131 c1_s (N_s) = x1_content(i, j, jpk-1)132 i_bottom = i_ocean[n_elements(i_ocean)-1] 133 134 z_s[N_s) = z_zw(i_bottom] 135 c1_s[N_s] = x1_content[i, j, jpk-1] 132 136 133 137 ; extraction d'un sous profil strictement croissant 134 mini = min(s_z (i_ocean))135 maxi = max(s_z (i_ocean))136 137 i_min = where(s_z (i_ocean)EQ mini)138 i_max = where(s_z (i_ocean)EQ maxi)138 mini = min(s_z[i_ocean]) 139 maxi = max(s_z[i_ocean]) 140 141 i_min = where(s_z[i_ocean] EQ mini) 142 i_max = where(s_z[i_ocean] EQ maxi) 139 143 ; on prend le plus grand des indices min 140 144 ; et le plus petit des indices max … … 148 152 ; l isopycne est mise en surface (z_s=0) 149 153 ; 150 ind = where(s_s LT s_z (i_min))154 ind = where(s_s LT s_z[i_min]) 151 155 IF ind[0] NE -1 THEN BEGIN 152 z_s (ind)= 0153 c1_s (ind)= 0156 z_s[ind] = 0 157 c1_s[ind] = 0 154 158 ENDIF 155 159 … … 157 161 ; l isopycne est mise au fond (z_s=z_zw(i_bottom)) 158 162 ; 159 ind = where(s_s GT s_z (i_max))163 ind = where(s_s GT s_z[i_max]) 160 164 IF ind[0] NE -1 THEN BEGIN 161 z_s (ind) = z_s(N_s)162 c1_s (ind) = c1_s(N_s)165 z_s[ind] = z_s[N_s] 166 c1_s[ind] = c1_s[N_s] 163 167 ENDIF 164 168 165 ind = where( (s_s GE s_z (i_min)) AND (s_s LT s_z(i_max)) )169 ind = where( (s_s GE s_z[i_min]) AND (s_s LT s_z[i_max]) ) 166 170 IF ind[0] NE -1 THEN BEGIN 167 171 168 i_profil = i_ocean (i_min:i_max)169 170 z_s (ind) = interpol(z_zt(i_profil), s_z(i_profil), s_s(ind))172 i_profil = i_ocean[i_min:i_max] 173 174 z_s[ind] = interpol(z_zt[i_profil], s_z[i_profil], s_s[ind]) 171 175 172 176 ; 173 177 ; j'utilise la fonction spline pour interpoler le contenu 174 178 ; 175 c1_s (ind) = spline(z_zw(i_profil), c1_z(i_profil), z_s(ind), 1)179 c1_s[ind] = spline(z_zw[i_profil], c1_z[i_profil], z_s[ind], 1) 176 180 ; 177 181 ; l'interpolation lineaire marche aussi. Elle donne des resultats differents. Elle 178 182 ; me semble moins propre. Je la donne en remarque. 179 183 ; 180 ; c_s (ind+1) = interpol(c_z(i_profil), z_zw(i_profil), z_s(ind+1))184 ; c_s[ind+1] = interpol(c_z[i_profil], z_zw[i_profil], z_s[ind+1]) 181 185 ; 182 186 ; Remarque : on ne divise pas par l'epaisseur de la couche 183 187 ; 184 x1_s (ind) = (c1_s(ind+1)-c1_s(ind));/(z_s(ind+1)-z_s(ind))188 x1_s[ind] = (c1_s[ind+1]-c1_s[ind]);/(z_s[ind+1]-z_s[ind]) 185 189 186 190 ENDIF … … 188 192 ENDIF 189 193 190 x1_bin (i, j, *)= x1_s194 x1_bin[i, j, *] = x1_s 191 195 192 196 ENDFOR … … 196 200 197 201 END 198 199 200 202 ;+ 201 203 ; … … 203 205 PRO compute_msf, u, v, msf $ 204 206 , MSK = msk 207 ; 208 compile_opt idl2, strictarrsubs 205 209 ; 206 210 ; … … 231 235 zumask = boundperio( (msk +shift(msk, -1, 0) ) < 1 ) 232 236 zvmask = boundperio( (msk +shift(msk, 0, -1) ) < 1 ) 233 zumask = zumask (*)#vert234 zvmask = zvmask (*)#vert237 zumask = zumask[*]#vert 238 zvmask = zvmask[*]#vert 235 239 zu = zu*zumask 236 240 zv = zv*zvmask … … 248 252 ; 249 253 offset = (fltarr(180, jpl)+long(jpi)*jpj)(*)#findgen(jpk2) 250 indu = indxu (*)#vert+offset251 indv = indxv (*)#vert+offset254 indu = indxu[*]#vert+offset 255 indv = indxv[*]#vert+offset 252 256 ; 253 257 ; extension sur la verticale 254 258 ; 255 zuleng = zxuleng (*)#vert256 zvleng = zxvleng (*)#vert259 zuleng = zxuleng[*]#vert 260 zvleng = zxvleng[*]#vert 257 261 ; 258 262 ; calcul du flux 259 263 ; 260 z = reform(zu (indu)*zuleng+zv(indv)*zvleng, 180, jpl, jpk2)264 z = reform(zu[indu]*zuleng+zv[indv]*zvleng, 180, jpl, jpk2) 261 265 ; 262 266 ; integration zonale du flux ! … … 270 274 ; couche). 271 275 ; 272 FOR k = jpk2-2, 0, -1 DO msf (*, k) = msf(*, k+1)-fm(*, k); *e3t(k)273 274 END 276 FOR k = jpk2-2, 0, -1 DO msf[*, k] = msf[*, k+1]-fm[*, k]; *e3t[k] 277 278 END -
trunk/tools/nadw.pro
r165 r169 10 10 PRO nadw, msf, nadw, lat, prof 11 11 ; 12 compile_opt idl2, strictarrsubs 13 ; 12 14 @common 13 15 ; … … 20 22 21 23 FOR k = 0, N-1 DO BEGIN 22 z = msf (*, *, k)23 nadw (k)= max(z) ; prend le max de msf24 z = msf[*, *, k] 25 nadw[k] = max(z) ; prend le max de msf 24 26 ; 25 27 ; 26 28 ; 27 29 indice = (where(z EQ flow(k)))[0] 28 lat(k) = gphit(90, indice MOD jpl) ; recupere la latitude 29 prof(k) = zdept(indice/jpl) ; et la profondeur 30 ; recupere la latitude 31 lat[k] = gphit[90, indice MOD jpl] 32 ; et la profondeur 33 prof[k] = zdept[indice/jpl] 30 34 ENDFOR 31 35 -
trunk/tools/nino_composite.pro
r165 r169 1 1 ;+ 2 ; 3 ; @todo 4 ; missing initorca0 module 2 5 ; 3 6 ; @version … … 7 10 PRO nino_composite $ 8 11 , FILENAME=filename 12 ; 13 compile_opt idl2, strictarrsubs 9 14 ; 10 15 dbase='/home/dynamite/mkolas/database/' … … 37 42 sstaC_temp = total(sstaC_sc, 2)*12/jpt 38 43 tab_sst_sc = fltarr(jpt) 39 FOR i = 0, jpt/12.-1 DO tab_sst_sc (i*12:i*12+11) = sstaC_temp(*)44 FOR i = 0, jpt/12.-1 DO tab_sst_sc[i*12:i*12+11] = sstaC_temp[*] 40 45 41 46 ;Moyenne glissante Trenberth 3 ou 5 mois ?? 42 sstaC = ts_smooth (sstaC-tab_sst_sc, 5)47 sstaC = ts_smooth[sstaC-tab_sst_sc, 5] 43 48 ; 44 49 resp=fltarr(jpt) … … 51 56 Deb_Nino = a(where(a-shift(a,1) GT 1)) 52 57 Fin_Nino = a(where(a-shift(a,1) GT 1)-1) 53 IF (Fin_Nino (0) LT Deb_Nino(0)) THEN Fin_Nino = Fin_Nino(1:n_elements(Fin_Nino)-1)54 IF(n_elements(Deb_Nino) GT n_elements(Fin_Nino)) THEN Deb_Nino = Deb_Nino (0:n_elements(Fin_Nino)-1)58 IF (Fin_Nino[0] LT Deb_Nino[0]) THEN Fin_Nino = Fin_Nino[1:n_elements(Fin_Nino)-1] 59 IF(n_elements(Deb_Nino) GT n_elements(Fin_Nino)) THEN Deb_Nino = Deb_Nino[0:n_elements(Fin_Nino)-1] 55 60 Length_Nino = Fin_Nino-Deb_Nino+1 56 Deb_Nino = Deb_Nino (where(Length_Nino GT 5))57 Fin_Nino = Fin_Nino (where(Length_Nino GT 5))58 Length_Nino = Length_Nino (where(Length_Nino GT 5))61 Deb_Nino = Deb_Nino[where(Length_Nino GT 5)] 62 Fin_Nino = Fin_Nino[where(Length_Nino GT 5)] 63 Length_Nino = Length_Nino[where(Length_Nino GT 5)] 59 64 60 65 ;;Max_Nino = intarr(n_elements(Deb_Nino)) 61 66 ;;FOR i = 0, n_elements(Deb_Nino)-1 DO BEGIN 62 ;;a = max(sstaC(Deb_Nino (i):Fin_Nino(i)),ind)63 ;;Max_Nino (i) = ind+Deb_Nino(i)67 ;;a = max(sstaC(Deb_Nino[i]:Fin_Nino[i]),ind) 68 ;;Max_Nino[i] = ind+Deb_Nino[i] 64 69 ;;ENDFOR 65 70 ;;Deb_Nino_sel = Deb_Nino(where((Max_Nino mod 12) GE 8 OR (Max_Nino mod 12) LE 1)) … … 73 78 74 79 FOR i=0,n_elements(Year_Nino)-1 DO BEGIN 75 sstaC_Nino (i, *) = sstaC(Year_Nino(i)*12: (Year_Nino(i)+2)*12-1)80 sstaC_Nino[i, *] = sstaC[Year_Nino[i]*12: (Year_Nino[i]+2)*12-1] 76 81 ENDFOR 77 82 78 83 jpt = 24 79 84 time = julday(01, 15, 01)+findgen(jpt)*30 80 plt1d,reform(sstaC_Nino (0,*)),'t', min = -stddev(sstaC)*3, max = stddev(sstaC)*3, petit = [2, 2, 1], /rempli,xtitle ='', ytitle = 'SSTA (C)', title = model+' '+exp+' '+strtrim(string(n_elements(Year_Nino_Weak)), 1)+' Weak Nino (PrecMAM<2.): Nino34 SSTA',subtitle = '', /port, win = 285 plt1d,reform(sstaC_Nino[0,*]),'t', min = -stddev(sstaC)*3, max = stddev(sstaC)*3, petit = [2, 2, 1], /rempli,xtitle ='', ytitle = 'SSTA (C)', title = model+' '+exp+' '+strtrim(string(n_elements(Year_Nino_Weak)), 1)+' Weak Nino (PrecMAM<2.): Nino34 SSTA',subtitle = '', /port, win = 2 81 86 FOR i=1,n_elements(Year_Nino)-1 DO BEGIN 82 plt1d,reform(sstaC_Nino (i,*)),'t', min = -stddev(sstaC)*3, max = stddev(sstaC)*3, /ov1d, /noerase, petit = [2, 2, 1], /rempli,xtitle ='', ytitle = '', title = '',subtitle = '', /port87 plt1d,reform(sstaC_Nino[i,*]),'t', min = -stddev(sstaC)*3, max = stddev(sstaC)*3, /ov1d, /noerase, petit = [2, 2, 1], /rempli,xtitle ='', ytitle = '', title = '',subtitle = '', /port 83 88 ENDFOR 84 89 -
trunk/tools/retouche.pro
r165 r169 39 39 FUNCTION retouche, tab $ 40 40 , TAILLEPOINT=taillepoint 41 ; 42 compile_opt idl2, strictarrsubs 41 43 ; 42 44 nouvmask=-1 -
trunk/tools/transfo/grads2cdf.pro
r165 r169 12 12 ;- 13 13 PRO grads2cdf 14 ; 15 compile_opt idl2, strictarrsubs 14 16 ; 15 17 @common … … 203 205 ; print, longname[0] 204 206 ; print, extract_str(line, '[', ']') 205 fid (in-1)= NCDF_VARDEF(id,varname[in-1], [xid, yid, tid], /FLOAT)206 NCDF_ATTPUT, id, fid (in-1), 'units', extract_str(line, '[', ']')207 NCDF_ATTPUT, id, fid (in-1), 'long_name', longname[0]208 NCDF_ATTPUT, id, fid (in-1), 'short_name', varname[in-1]207 fid[in-1] = NCDF_VARDEF(id,varname[in-1], [xid, yid, tid], /FLOAT) 208 NCDF_ATTPUT, id, fid[in-1], 'units', extract_str(line, '[', ']') 209 NCDF_ATTPUT, id, fid[in-1], 'long_name', longname[0] 210 NCDF_ATTPUT, id, fid[in-1], 'short_name', varname[in-1] 209 211 ENDFOR 210 212 … … 227 229 ; print, longname[0] 228 230 ; print, extract_str(line, '[', ']') 229 fid (idx+in-1)= NCDF_VARDEF(id,varname[idx+in-1], [xid, yid, tid], /FLOAT)231 fid[idx+in-1] = NCDF_VARDEF(id,varname[idx+in-1], [xid, yid, tid], /FLOAT) 230 232 NCDF_ATTPUT, id, fid(idx+in-1), 'units', extract_str(line, '[', ']') 231 233 NCDF_ATTPUT, id, fid(idx+in-1), 'long_name', longname[0] … … 252 254 ; print, longname[0] 253 255 ; print, extract_str(line, '[', ']') 254 fid (idx+in-1)= NCDF_VARDEF(id,varname[idx+in-1], [xid, yid, zid, tid], /FLOAT)256 fid[idx+in-1] = NCDF_VARDEF(id,varname[idx+in-1], [xid, yid, zid, tid], /FLOAT) 255 257 NCDF_ATTPUT, id, fid(idx+in-1), 'units', extract_str(line, '[', ']') 256 258 NCDF_ATTPUT, id, fid(idx+in-1), 'long_name', longname[0] … … 289 291 field = grads_read(filename, varname[in-1], time_2 = it, NCDF_DB = iodir, /all_data, /no_key_shift) 290 292 tidx = ntime*(month-1)+it-1 291 fldata (in-1, *, tidx)= field.data293 fldata[in-1, *, tidx] = field.data 292 294 ENDFOR 293 295 ENDFOR … … 303 305 field = grads_read(filename, varname[in-1], time_2 = it, NCDF_DB = iodir, /all_data, /no_key_shift) 304 306 tidx = ntime*(month-1)+it-1 305 fldata (in-1, *, tidx)= field.data307 fldata[in-1, *, tidx] = field.data 306 308 ENDFOR 307 309 ENDFOR … … 320 322 field = grads_read(filename, varname[idx+in-1], time_2 = it, NCDF_DB = iodir, /all_data, /no_key_shift, level = il) 321 323 tidx = ntime*(month-1)+it-1 322 fldata3 (in-1, *, il-1, tidx)= field.data324 fldata3[in-1, *, il-1, tidx] = field.data 323 325 ENDFOR 324 326 ENDFOR … … 333 335 FOR in = 1, nfields DO BEGIN 334 336 335 NCDF_ATTPUT, id, fid (in-1), 'valid_min', min(fldata(in-1, *, *)), /float336 NCDF_ATTPUT, id, fid (in-1), 'valid_max', max(fldata(in-1, *, *)), /float337 NCDF_ATTPUT, id, fid (in-1), 'missing_value',1.e20, /float337 NCDF_ATTPUT, id, fid[in-1], 'valid_min', min(fldata[in-1, *, *]), /float 338 NCDF_ATTPUT, id, fid[in-1], 'valid_max', max(fldata[in-1, *, *]), /float 339 NCDF_ATTPUT, id, fid[in-1], 'missing_value',1.e20, /float 338 340 339 341 ENDFOR … … 341 343 FOR in = 1, nfields3 DO BEGIN 342 344 343 NCDF_ATTPUT, id, fid(idx+in-1), 'valid_min', min(fldata3 (in-1, *, *, *)), /float344 NCDF_ATTPUT, id, fid(idx+in-1), 'valid_max', max(fldata3 (in-1, *, *, *)), /float345 NCDF_ATTPUT, id, fid(idx+in-1), 'valid_min', min(fldata3[in-1, *, *, *]), /float 346 NCDF_ATTPUT, id, fid(idx+in-1), 'valid_max', max(fldata3[in-1, *, *, *]), /float 345 347 NCDF_ATTPUT, id, fid(idx+in-1), 'missing_value',1.e20, /float 346 348 … … 365 367 FOR in = 1, nfields DO BEGIN 366 368 367 fieldm = reform(fldata (in-1, *, *), jpi, jpj, ndates)369 fieldm = reform(fldata[in-1, *, *], jpi, jpj, ndates) 368 370 print, 'write ', varname[in-1] 369 371 370 NCDF_VARPUT, id, fid (in-1), fieldm372 NCDF_VARPUT, id, fid[in-1], fieldm 371 373 372 374 ENDFOR 373 375 FOR in = 1, nfields3 DO BEGIN 374 376 375 fieldm = reform(fldata3 (in-1, *, *, *), jpi, jpj, nlevs, ndates)377 fieldm = reform(fldata3[in-1, *, *, *], jpi, jpj, nlevs, ndates) 376 378 print, 'write ', varname[nfields2d + nfields_sf+in-1] 377 379 -
trunk/tools/transfo/grads_read.pro
r168 r169 17 17 , NO_KEY_SHIFT=no_key_shift $ 18 18 , LEVEL=level 19 ; 20 compile_opt idl2, strictarrsubs 19 21 ; 20 22 @common -
trunk/usr/correl.pro
r165 r169 6 6 ;- 7 7 PRO correl 8 ; 9 compile_opt idl2, strictarrsubs 8 10 ; 9 11 @common -
trunk/usr/p.pro
r165 r169 7 7 PRO p 8 8 ; 9 compile_opt idl2, strictarrsubs 10 ; 9 11 resolve_routine, 'post_it' 10 12 post_it -
trunk/usr/plt_def.pro
r165 r169 8 8 ;- 9 9 PRO plt_def 10 ; 11 compile_opt idl2, strictarrsubs 10 12 ; 11 13 @com_eg -
trunk/usr/post_it.pro
r165 r169 11 11 ;- 12 12 PRO post_it 13 ; 14 compile_opt idl2, strictarrsubs 13 15 ; 14 16 @common -
trunk/usr/post_it2.pro
r165 r169 6 6 ;- 7 7 PRO post_it2 8 ; 9 compile_opt idl2, strictarrsubs 8 10 ; 9 11 @common
Note: See TracChangeset
for help on using the changeset viewer.