source: trunk/UTIL/NCL/visu-utils.ncl @ 2236

Last change on this file since 2236 was 854, checked in by slebonnois, 12 years ago

SL: correction d'un bug pour Titan (zlsm1) + quelques petits ajustements
File size: 16.5 KB
RevLine 
[497]1;;;;;;;;;;;; positionDims ;;;;;;;;;;;;;;;;
2
3undef("positionDims")
4function positionDims(var,nbdim,names)
5;============================
6
7; var   : variable
8; nbdim : number of dimensions for var
9; names : names of the 4 dimensions: time,levs,lat,lon
10
11begin
12
13  flag=(/10,10,10,10/)
14  do i=0,nbdim-1
15     if (var!i.eq.names(0)) then
16        flag(0)=i
17     end if
18     if (var!i.eq.names(1)) then
19        flag(1)=i
20     end if
21     if (var!i.eq.names(2)) then
22        flag(2)=i
23     end if
24     if (var!i.eq.names(3)) then
25        flag(3)=i
26     end if
27  end do
28
29  return(flag)
30
31end
32
33;;;;;;;;;;;; trimVar ;;;;;;;;;;;;;;;;
34
35undef("trimVar")
36function trimVar(var,nbdim,dimname,mindimval,maxdimval)
37;============================
38
39; var      : variable
40; nbdim    : number of dimensions for var
41; dimname  : names of the dimensions for this var
42; mindimval: min values for the dimensions
43; maxdimval: max values for the dimensions
44
45begin
46
47  if (nbdim.eq.1) then
48     vartrim=var({$dimname(0)$|mindimval(0):maxdimval(0)})
49  end if
50  if (nbdim.eq.2) then
51       vartrim=var({$dimname(0)$|mindimval(0):maxdimval(0)},{$dimname(1)$|mindimval(1):maxdimval(1)})
52  end if
53  if (nbdim.eq.3) then
54       vartrim=var({$dimname(0)$|mindimval(0):maxdimval(0)},{$dimname(1)$|mindimval(1):maxdimval(1)},{$dimname(2)$|mindimval(2):maxdimval(2)})
55  end if
56  if (nbdim.eq.4) then
57       vartrim=var({$dimname(0)$|mindimval(0):maxdimval(0)},{$dimname(1)$|mindimval(1):maxdimval(1)},{$dimname(2)$|mindimval(2):maxdimval(2)},{$dimname(3)$|mindimval(3):maxdimval(3)})
[808]58;       vartrim=var($dimname(0)$|0:0,{$dimname(1)$|mindimval(1):maxdimval(1)},{$dimname(2)$|mindimval(2):maxdimval(2)},{$dimname(3)$|mindimval(3):maxdimval(3)})
[497]59  end if
60
61  return(vartrim)
62
63end
64
65;;;;;;;;;;;; avgVar ;;;;;;;;;;;;;;;;
66
67undef("avgVar")
68function avgVar(vartrim,nbdim,dimavg)
69;============================
70
71; vartrim : variable
72; nbdim   : number of dimensions for var
73; dimavg  : logicals for averaging the dimensions
74
75begin
76
77  if ((nbdim.eq.4).and.(dimavg(3))) then
78       vartrim_3=dim_avg_n_Wrap(vartrim,3)
79  else
80       vartrim_3=vartrim
81  end if
82  if ((nbdim.ge.3).and.(dimavg(2))) then
83       vartrim_2=dim_avg_n_Wrap(vartrim_3,2)
84  else
85       vartrim_2=vartrim_3
86  end if
87  if ((nbdim.ge.2).and.(dimavg(1))) then
88       vartrim_1=dim_avg_n_Wrap(vartrim_2,1)
89  else
90       vartrim_1=vartrim_2
91  end if
92  if ((nbdim.ge.1).and.(dimavg(0))) then
93       varavg=dim_avg_n_Wrap(vartrim_1,0)
94  else
95       varavg=vartrim_1
96  end if
97
98  return(varavg)
99
100end
101
102;;;;;;;;;;;; deflimits_2d ;;;;;;;;;;;;;;;;
103
104undef("deflimits_2d")
105function deflimits_2d(var,res)
106;============================
107
108; var : 2D variable
109; res : ressource
110
111begin
112
113  dims=dimsizes(var)
114  dimname=getvardims(var)
115  res@trXMinF = tofloat(var&$dimname(1)$(0))
116  res@trXMaxF = tofloat(var&$dimname(1)$(dims(1)-1))
117  res@trYMinF = tofloat(var&$dimname(0)$(0))
118  res@trYMaxF = tofloat(var&$dimname(0)$(dims(0)-1))
119
120  return(res)
121
122end
123
124;;;;;;;;;;;; contours_2d ;;;;;;;;;;;;;;;;
125
126undef("contours_2d")
127function contours_2d(res,limits)
128;============================
129
130; res    : modified ressource for the plot
131; limits : min, max and values for contours (-888 => automatic)
132
133begin
134
135  if (limits(0).ne.-888) then
136     res@cnLevelSelectionMode = "ManualLevels"
137     res@cnMinLevelValF       = limits(0)
138     res@cnMaxLevelValF       = limits(1)
139     if (limits(2).ne.-888) then
140        res@cnLevelSpacingF   = limits(2)
141     end if
142  end if
143
144  return(res)
145
146end
147
148;;;;;;;;;;;; optionsOrtho ;;;;;;;;;;;;;;;;
149
150undef("optionsOrtho")
151function optionsOrtho(res,center,grid,spacing,gridcol)
152;============================
153
154; res    : modified ressource for the plot
155; center : lon,lat coordinates for center of projection
156; grid   : logical for printing grid
157; spacing: lon,lat spacing for grid
158; gridcol: color for grid
159
160begin
161
162  res@mpCenterLonF      = center(0)  ; choose center lon
163  res@mpCenterLatF      = center(1)  ; choose center lat
164  res@mpGridAndLimbOn   = grid
165  res@mpGridLonSpacingF = spacing(0)
166  res@mpGridLatSpacingF = spacing(1)
167  res@mpGridLineColor   = gridcol
168
169  return(res)
170
171end
172
173;;;;;;;;;;;; optionsStereo ;;;;;;;;;;;;;;;;
174
175undef("optionsStereo")
176function optionsStereo(res,hm,grid,spacing,gridcol)
177;============================
178
179; res    : modified ressource for the plot
180; hm     : hemisphere to be plotted
181; grid   : logical for printing grid
182; spacing: lon,lat spacing for grid
183; gridcol: color for grid
184
185begin
186
187  res@gsnPolar          = hm
188
189  res@mpGridAndLimbOn   = grid
190  res@mpGridLonSpacingF = spacing(0)
191  res@mpGridLatSpacingF = spacing(1)
192  res@mpGridLineColor   = gridcol
193
194  return(res)
195
196end
197
198;;;;;;;;;;;; optionsOverPlot_2d ;;;;;;;;;;;;;;;;
199
200undef("optionsOverPlot_2d")
201function optionsOverPlot_2d(res)
202;============================
203
204; res    : modified ressource for the plot
205
206begin
207
208  res@cnLineLabelFormat            = "@5.1f"
209  res@cnLineLabelsOn               = True
210  res@cnLabelMasking               = True               
211  res@cnLineLabelBackgroundColor   = -1
212  res@cnInfoLabelOn                = False
213  res@gsnContourNegLineDashPattern = 1 
214
215  return(res)
216
217end
218
219;;;;;;;;;;;; axePress_2d ;;;;;;;;;;;;;;;;
220
221undef("axePress_2d")
222function axePress_2d(res)
223;============================
224
225; res : modified ressource for the plot
226
227begin
228
229  res@tiYAxisString  = "Pressure"
230  res@trYReverse     = True
231  res@gsnYAxisIrregular2Log = True
232
233  return(res)
234
235end
236
237;;;;;;;;;;;; axePress_1d ;;;;;;;;;;;;;;;;
238
239undef("axePress_1d")
240function axePress_1d(res,axe,limits)
241;============================
242
243; res    : modified ressource for the plot
244; axe    : Y axis for plot (pressure)
245; limits : min and max values for variable axis
246
247begin
248
249  res@tiYAxisString  = "Pressure"
250  res@trYReverse     = True
251  res@xyYStyle       = "Log"
252  res@trYMinF         = axe(0)
253  res@trYMaxF         = axe(dimsizes(axe)-1)
254  if (limits(0).ne.-888) then
255    res@trXMinF       = limits(0)
256    res@trXMaxF       = limits(1)
257  end if
258
259  return(res)
260
261end
262
263
264undef("axeNoPress_1d")
265function axeNoPress_1d(res,axe,limits)
266;============================
267
268; res    : modified ressource for the plot
269; axe    : X axis for plot (not pressure)
270; limits : min and max values for variable axis
271
272begin
273
274  res@trXMinF         = axe(0)
275  res@trXMaxF         = axe(dimsizes(axe)-1)
276  if (limits(0).ne.-888) then
277    res@trYMinF       = limits(0)
278    res@trYMaxF       = limits(1)
279  end if
280
281  return(res)
282
283end
284
285;;;;;;;;;;;; axeTimeX_2d ;;;;;;;;;;;;;;;;
286
287undef("axeTimeX_2d")
288function axeTimeX_2d(res,axe)
289;============================
290
291; res : modified ressource for the plot
292; axe : X axis is time
293
294begin
295
296  res@tiXAxisString        = "Time (local days)"
297  res@sfXArray             = axe
298  res@tmXBLabelFontHeightF = 0.015
299  res@trXMinF              = axe(0)
300  res@trXMaxF              = axe(dimsizes(axe)-1)
301
302  return(res)
303
304end
305
306;;;;;;;;;;;; axeTimeXls_2d ;;;;;;;;;;;;;;;;
307
308undef("axeTimeXls_2d")
309function axeTimeXls_2d(res,axe)
310;============================
311
312; res : modified ressource for the plot
313; axe : X axis is Ls
314
315begin
316
317  res@tiXAxisString            = "Solar longitude"
318  res@sfXArray                 = axe
319  res@gsnXAxisIrregular2Linear = True
320  res@tmXBMode                 = "Manual"       
321  res@tmXBTickStartF           = 0
322  res@tmXBTickEndF             = 360
323  res@tmXBTickSpacingF         = 30
324  res@tmXBMinorPerMajor        = 2
325  res@tmXBLabelFontHeightF     = 0.015
326  res@trXMinF                  = axe(0)
327  res@trXMaxF                  = axe(dimsizes(axe)-1)
328
329  return(res)
330
331end
332
333;;;;;;;;;;;; axeTimeY_2d ;;;;;;;;;;;;;;;;
334
335undef("axeTimeY_2d")
336function axeTimeY_2d(res,axe)
337;============================
338
339; res : modified ressource for the plot
340; axe : Y axis is time
341
342begin
343
344  res@tiYAxisString        = "Time (local days)"
345  res@sfYArray             = axe
346  res@tmYLLabelFontHeightF = 0.015
347  res@trYMinF              = axe(0)
348  res@trYMaxF              = axe(dimsizes(axe)-1)
349
350  return(res)
351
352end
353
354;;;;;;;;;;;; axeTimeYls_2d ;;;;;;;;;;;;;;;;
355
356undef("axeTimeYls_2d")
357function axeTimeYls_2d(res,axe)
358;============================
359
360; res : modified ressource for the plot
361; axe : Y axis is Ls
362
363begin
364
365  res@tiYAxisString            = "Solar longitude"
366  res@sfYArray                 = axe
367  res@gsnYAxisIrregular2Linear = True
368  res@tmYLMode                 = "Manual"       
369  res@tmYLTickStartF           = 0
370  res@tmYLTickEndF             = 360
371  res@tmYLTickSpacingF         = 30
372  res@tmYLMinorPerMajor        = 2
373  res@tmYLLabelFontHeightF     = 0.015
374  res@trYMinF                  = axe(0)
375  res@trYMaxF                  = axe(dimsizes(axe)-1)
376
377  return(res)
378
379end
380
381;;;;;;;;;;;; axeTime_1d ;;;;;;;;;;;;;;;;
382
383undef("axeTime_1d")
384function axeTime_1d(res,axe)
385;============================
386
387; res : modified ressource for the plot
388; axe : X axis is time
389
390begin
391
392  res@tiXAxisString        = "Time (local days)"
393  res@tmXBLabelFontHeightF = 0.015
394
395  return(res)
396
397end
398
399;;;;;;;;;;;; axeTimels_1d ;;;;;;;;;;;;;;;;
400
401undef("axeTimels_1d")
402function axeTimels_1d(res,axe)
403;============================
404
405; res : modified ressource for the plot
406; axe : X axis is Ls
407
408begin
409
410  res@tiXAxisString            = "Solar longitude"
411  res@tmXBMode                 = "Manual"       
412  res@tmXBTickStartF           = 0
413  res@tmXBTickEndF             = 360
414  res@tmXBTickSpacingF         = 30
415  res@tmXBMinorPerMajor        = 2
416  res@tmXBLabelFontHeightF     = 0.015
417
418  return(res)
419
420end
421
422;;;;;;;;;;;; inversDim_2d ;;;;;;;;;;;;;;;;
423
424undef("inversDim_2d")
425function inversDim_2d(var)
426;============================
427
428; var is a 2d variable
429
430begin
431
432  dimname=getvardims(var)
433  tmpvar = var($dimname(1)$|:,$dimname(0)$|:)
434
435  return(tmpvar)
436
437end
438
439;;;;;;;;;;;; revertLon_2dmap ;;;;;;;;;;;;;;;;
440
441undef("revertLon_2dmap")
442function revertLon_2dmap(var)
443;============================
444
445; var is a 2d lat-lon variable
446
447begin
448
449  tmpvar  = var
450  dimname = getvardims(var)
451  Lname   = dimname(1)
452  nbL     = dimsizes(var&$Lname$)
453  tmpaxL  = new(nbL,float)
454
455;revert lon axis:
456  do i=0,nbL-1
457    tmpaxL(i)=var&$Lname$(nbL-1-i)
458  end do
459  tmpvar&$Lname$ = tmpaxL
460
461;revert values along lon axis:
462  do i=0,nbL-1
463    tmpvar(:,i) = var(:,nbL-1-i)
464  end do
465
466  return(tmpvar)
467
468end
469
470;;;;;;;;;;;; Map_2d ;;;;;;;;;;;;;;;;
471
472undef("Map_2d")
473function Map_2d(var,wks,res)
474;============================
475
476; var : 2D variable
477; wks : workstation
478; res : ressource
479
480begin
481
482; main plot characteristics
483
484  res@gsnSpreadColors = True
485  res@tmXTOn          = False
486  res@cnFillOn        = True
487  res@cnLineLabelsOn  = True
488  res@cnLabelMasking  = True
489  res@cnLineLabelBackgroundColor = "white"
490  res@lbOrientation   = "vertical"
491
492; Create plot
493
494  plotvar  = gsn_csm_contour(wks, var, res )
495
496  return(plotvar)
497
498end
499
500;;;;;;;;;;;; OPlot_2d ;;;;;;;;;;;;;;;;
501
502undef("OPlot_2d")
503function OPlot_2d(var,wks,res)
504;============================
505
506; var : 2D variable
507; wks : workstation
508; res : ressource
509
510begin
511
512; main plot characteristics
513
514  res@tmXTOn          = False
515  res@tmYROn          = False
516  res@cnFillOn        = False
517
518  res@gsnLeftString   = ""
519  res@gsnRightString  = ""
520
521; Create plot
522
523  plotvar  = gsn_csm_contour(wks, var, res )
524
525  return(plotvar)
526
527end
528
529;;;;;;;;;;;; Plot_1d ;;;;;;;;;;;;;;;;
530
531undef("Plot_1d")
532function Plot_1d(var,wks,res,axe,axep,limits)
533;============================
534
535; var            : 1D variable
536; wks            : workstation
537; res            : ressource
538; axe            : axis for plot
539; axep           : logical indication for pressure axis
540; limits         : min and max values for variable axis
541
542begin
543
544; main plot characteristics
545
546; Create plot
547; pressure axis ?
548  if (axep) then
549    res=axePress_1d(res,axe,limits)
550    plotvar  = gsn_csm_xy(wks, var, axe, res )
551  else
552    res1d=axeNoPress_1d(res,axe,limits)
553    plotvar  = gsn_csm_xy(wks, axe, var, res )
554  end if
555
556  return(plotvar)
557
558end
559
560;;;;;;;;;;;; Ortho_2d ;;;;;;;;;;;;;;;;
561
562undef("Ortho_2d")
563function Ortho_2d(var,wks,res)
564;============================
565
566; var : 2D variable
567; wks : workstation
568; res : ressource
569
570begin
571
572; main plot characteristics
573
574  res@gsnSpreadColors = True
575  res@tmXTOn          = False
576  res@cnFillOn        = True
577  res@cnLineLabelsOn  = True
578  res@cnLabelMasking  = True
579  res@cnLineLabelBackgroundColor = "white"
580  res@lbOrientation   = "vertical"
581
582  res@mpProjection      = "Orthographic"       ; choose projection
583  res@mpOutlineOn       = False
584  res@mpPerimOn         = False             ; turn off box around plot
585  res@mpFillOn          = False
586  res@mpGridLineDashPattern = 2
587 
588; Create plot
589
590  plotvar  = gsn_csm_contour_map(wks, var, res )
591
592  return(plotvar)
593
594end
595
596;;;;;;;;;;;; Stereo_2d ;;;;;;;;;;;;;;;;
597
598undef("Stereo_2d")
599function Stereo_2d(var,wks,res)
600;============================
601
602; var : 2D variable
603; wks : workstation
604; res : ressource
605
606begin
607
608; main plot characteristics
609
610  res@gsnSpreadColors = True
611  res@tmXTOn          = False
612  res@cnFillOn        = True
613  res@cnLineLabelsOn  = True
614  res@cnLabelMasking  = True
615  res@cnLineLabelBackgroundColor = "white"
616  res@lbOrientation   = "vertical"
617
618  res@mpProjection      = "Stereographic"       ; choose projection
619  res@mpOutlineOn       = False
620  res@mpPerimOn         = False             ; turn off box around plot
621  res@mpFillOn          = False
622 
623; Create plot
624
625  plotvar  = gsn_csm_contour_map_polar(wks, var, res )
626
627  return(plotvar)
628
629end
630
[808]631;;;;;;;;;;;; customVar ;;;;;;;;;;;;;;;;
632
633undef("customVar")
634function customVar(infile,labelvar)
635;============================
636
637; infile   : input file
638; labelvar : input label for variable
639
640; THIS IS TO BE CUSTOMIZED EACH TIME ACCORDING TO YOUR NEEDS
641; because it is impossible to make it automatic...
642
643begin
644
645prepared = "dummy"    ; DEFAULT
646
[854]647if (labelvar.eq."ustar") then
648  if (isfilevar(infile,"vitu")) then
649    u = infile->vitu
650  else
651    u = infile->U
652  end if
653  myvar = dim_rmvmean_n_Wrap(u,3)
654  myvar@long_name="Deviation from zonal-mean zonal wind"
655  prepared = labelvar
656end if
657
658if (labelvar.eq."uprime") then
659  if (isfilevar(infile,"vitu")) then
660    u = infile->vitu
661  else
662    u = infile->U
663  end if
664  myvar = dim_rmvmean_n_Wrap(u,0)
665  myvar@long_name="Deviation from time-mean zonal wind"
666  prepared = labelvar
667end if
668
[808]669if (labelvar.eq."dTdyn") then
670  if (isfilevar(infile,"vitu")) then
671    dtdyn = infile->dtdyn
672  else
673    dtdyn = infile->DTCORE
674  end if
675  myvar = dtdyn*1e7
676  copy_VarMeta(dtdyn,myvar)
677  myvar@units    ="K/Vd"
678  myvar@long_name="Dynamics heating rate"
679  prepared = labelvar
680end if
681
682if (labelvar.eq."dTadj") then
683  if (isfilevar(infile,"vitu")) then
684    dtadj = infile->dtajs
685  else
686    dtadj = infile->DADJ
687  end if
688  myvar = dtadj*1e7
689  copy_VarMeta(dtadj,myvar)
690  myvar@units    ="K/Vd"
691  myvar@long_name="Dry convection heating rate"
692  prepared = labelvar
693end if
694
695if (labelvar.eq."dTvdf") then
696  if (isfilevar(infile,"vitu")) then
697    dtvdf = infile->dtvdf
698  else
699    dtvdf = infile->DTV
700  end if
701  myvar = dtvdf*1e7
702  copy_VarMeta(dtvdf,myvar)
703  myvar@units    ="K/Vd"
704  myvar@long_name="PBL heating rate"
705  prepared = labelvar
706end if
707
708if (labelvar.eq."dTrad") then
709  if (isfilevar(infile,"vitu")) then
710    dtlwr = infile->dtlwr
711    dtswr = infile->dtswr
712  else
713    dtlwr = infile->QRL
714    dtswr = infile->QRS
715  end if
716  myvar = (dtlwr+dtswr)*1e7
717  copy_VarMeta(dtlwr,myvar)
718  myvar@units    ="K/Vd"
719  myvar@long_name="Radiative balance"
720  prepared = labelvar
721end if
722
723if (labelvar.eq."Tbal") then
724  if (isfilevar(infile,"vitu")) then
725    dtlwr = infile->dtlwr
726    dtswr = infile->dtswr
727    dtajs = infile->dtajs
728    dtdyn = infile->dtdyn
729    dtvdf = infile->dtvdf
730  else
731    dtlwr = infile->QRL
732    dtswr = infile->QRS
733    dtajs = infile->DADJ
734    dtdyn = infile->DTCORE
735    dtvdf = infile->DTV
736  end if
737  myvar = (dtlwr+dtswr+dtajs+dtdyn+dtvdf)*1e7
738  copy_VarMeta(dtlwr,myvar)
739  myvar@units    ="K/Vd"
740  myvar@long_name="Thermal balance"
741  prepared = labelvar
742end if
743
744if (labelvar.eq."Tforc") then
745  if (isfilevar(infile,"vitu")) then
746    dtlwr = infile->dtlwr
747    dtswr = infile->dtswr
748    dtajs = infile->dtajs
749    if (isfilevar(infile,"dtvdf")) then
750          dtvdf = infile->dtvdf
751    else
752          dtvdf = dtajs*0.
753    end if
754  else
755    dtlwr = infile->QRL
756    dtswr = infile->QRS
757    dtajs = infile->DADJ
758    if (isfilevar(infile,"DTV")) then
759          dtvdf = infile->DTV
760    else
761          dtvdf = dtajs*0.
762    end if
763  end if
764  myvar = (dtlwr+dtswr+dtajs+dtvdf)*1e7
765  copy_VarMeta(dtlwr,myvar)
766  myvar@units    ="K/Vd"
767  myvar@long_name="Thermal physical forcing"
768  prepared = labelvar
769end if
770
771if (labelvar.eq."ske") then
772  if (isfilevar(infile,"vitu")) then
773    u = infile->vitu
774    v = infile->vitv
775  else
776    u = infile->U
777    v = infile->V
778  end if
779  myvar = u*u+v*v
780  copy_VarMeta(u,myvar)
781  myvar@units    ="m2/s2"
782  myvar@long_name="Specific kinetic energy"
783  prepared = labelvar
784end if
785
786if (prepared.eq."dummy") then   ; DEFAULT: LEAVE AS IS
787  print("You chose a customized variable !")
788  print("Modify the function customVar in visu-utils.ncl to prepare this variable")
789  exit
790end if
791
792  return(myvar)
793
794end
795
Note: See TracBrowser for help on using the repository browser.