I quickly looked at your code. It looks like you are programming in C or
fortran? Your are prealloccating memory (malloc, allocate) and the looping
(for, do). All NCL functions automatically allocate array space and there
are functions to average and sum. Use of built-in functions makes for
cleaner and more efficient code. Please read the mini-Language manual at:
http://www.ncl.ucar.edu/Document/Manuals/
====
-->[1]<============
Klimaklasse = new(dimsizes(Niederschlag),float)
Klimaklasse = 99999.9
can be done in one line
Klimaklasse = new(dimsizes(Niederschlag),float, 99999.9)
-->[2]<============
temp = new(12,typeof(Temperatur))
temp = Temperatur(:,0,loc_lat,loc_lon)
nied = new(12,typeof(Niederschlag))
nied = Niederschlag(:,loc_lat,loc_lon)
summe_nied = 0.0
mean_temp = 0.0
do i = 0,11
summe_nied = summe_nied + nied(i)
mean_temp = mean_temp + temp(i)
end do
mean_temp = mean_temp / 12.0
can be replaced by:
For example
mean_temp = dim_avg_n(Temperatur(:,0,loc_lat,loc_lon) , 0)
summe_nied = dim_sum_n( Niederschlag(:,loc_lat,loc_lon), 0)
or, if you stored the variables in local arrays
mean_temp = dim_avg_n(temp , 0)
summe_nied = dim_sum_n(nied , 0)
or, if meta data are desired
mean_temp = dim_avg_n_Wrap(Temperatur(:,0,loc_lat,loc_lon), 0)
summe_nied = dim_sum_n_Wrap(Niederschlag(:,loc_lat,loc_lon), 0)
=======
...
do k=0,11
if(temp(k) .gt. 10.0) then
warm = warm + 1.0
end if
end do
can be replaced with
warm = dim_num_n(temp.gt.10.0, 0)
...
if(summe_nied / 10.0 .gt. mean_temp + regen) then
klasse = klasse + "S" ; Steppenklima
else
klasse = klasse + "W" ; Wuestenklima
end if
if(mean_temp .gt. 18.0) then
klasse = klasse + "h" ; heiß
else if(mean_temp .le. 18.0) then
klasse = klasse + "k" ; kalt
end if
replaced by
klasse = where(summe_nied / 10.0 .gt. mean_temp + regen \
,klasse + "S", klasse + "W")
klasse = where(mean_temp .gt. 18.0, klasse + "h", klasse + "k")
etc.
++++++++
Please look at the NCL function list.
Good Luck
On Sat, Jun 21, 2014 at 3:58 PM, Karin Meier-Fleischer <
meier-fleischer@dkrz.de> wrote:
> Hi Stephan,
>
> the do loops with loc_lat and loc_lon are wrong. You have to change it to
>
> do loc_lat = 0,dimsizes(lat2d(:,0))-1
> do loc_lon = 0,dimsizes(lon2d(0,:))-1
>
> To define the variable Klimaklasse with all attributes and dimensions of
> the variable Niederschlag you shouldn't use 'new' but copy it:
>
> Klimaklasse = Niederschlag
>
> The changes can be found in your attached script.
>
> Bye,
> Karin
>
> Am 21.06.14 19:15, schrieb Stephan Herrmann:
>
> Hello everybody,
>
>
>
> I want to calculate and plot the climate zones after the Köppen & Ge=
iger
> climate classification with the following NCL code:
>
>
>
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
>
> begin
>
> f =
> addfile("MPIM-REMO_EH5OMA1Br3_HR_044_BOL_1961-1990_167_minus_1_hour_Celsi=
ustt_ymonmean.nc","r")
> Temperatur = f->var167(:,:,:,:)
>
> latsize = getfilevardimsizes(f,"rlat")
> lonsize = getfilevardimsizes(f,"rlon")
> rotpole = f->rotated_pole ; need attributes for correct remap
>
> ll = asciiread("ll.out",(/latsize,lonsize,2/),"float")
> lat2d = ll(::,:,0)
> lon2d = ll(::,:,1)
>
> Temperatur@lat2d = lat2d
> Temperatur@lon2d = lon2d
>
> f =
> addfile("MPIM-REMO_EH5OMA1Br3_HR_044_BOL_1961-1990_142_143_minus_1_hour_y=
monsum_mean_30_years.nc","r")
> Niederschlag = f->var142(:,:,:)
>
> latsize = getfilevardimsizes(f,"rlat")
> lonsize = getfilevardimsizes(f,"rlon")
> rotpole = f->rotated_pole ; need attributes for correct remap
>
> ll = asciiread("ll.out",(/latsize,lonsize,2/),"float")
> lat2d = ll(::,:,0)
> lon2d = ll(::,:,1)
>
> Niederschlag@lat2d = lat2d
> Niederschlag@lon2d = lon2d
>
> Klimaklasse = new(dimsizes(Niederschlag),float)
>
> Klimaklasse = 99999.9
>
> do loc_lat = 0,max(dimsizes(lat2d))-1
> do loc_lon = 0,max(dimsizes(lon2d))-1
>
> temp = new(12,typeof(Temperatur))
> temp = Temperatur(:,0,loc_lat,loc_lon)
>
> nied = new(12,typeof(Niederschlag))
> nied = Niederschlag(:,loc_lat,loc_lon)
>
> summe_nied = 0.0
> mean_temp = 0.0
> do i = 0,11
> summe_nied = summe_nied + nied(i)
> mean_temp = mean_temp + temp(i)
> end do
> mean_temp = mean_temp / 12.0
>
> ;**** Klimaklassifikation nach Koeppen-Geiger
>
> hk = 6
>
> sommer_nied = 0.0
>
> do k = min((/mod(3 + hk, 12), mod(8 + hk, 12)/)), max((/mod(3 + hk, 12=
),
> mod(8 + hk, 12)/)) ; oder hydrologisch Mai bis Oktober?
> sommer_nied = sommer_nied + nied(k)
> end do
>
> if(summe_nied * 0.7 .lt. sommer_nied) then
> regen = 14.0 ; Sommerregen
> else if(summe_nied * 0.3 .gt. sommer_nied) then
> regen = 0.0 ; Winterregen
> else
> regen = 7.0 ; gleichmäßige Niederschlagsver=
teilung
> end if
> end if
>
> warm = 0.0
>
> do k=0,11
> if(temp(k) .gt. 10.0) then
> warm = warm + 1.0
> end if
> end do
>
> if((min(temp(:)) .gt. 18.0) .and. (summe_nied / 10.0 .gt. 2.0 *
> (mean_temp + regen))) then
> klasse = "A" ; Tropisches Regenwald- oder Savannenkli=
ma
> ohne Winter
>
> if(min(nied(:)) .gt. 60.0) then
> klasse = klasse + "f" ; Tropisches Regenwaldklima
> else if (min(nied(:)) .gt. 0.04 * (2400.0 - summe_nied)) then ;aender=
n
> siehe Zettel
> klasse = klasse + "m" ; Mittelform, Regenwaldklima trotz einer
> Trockenzeit
> else
> klasse = klasse + "w" ; Savannenklima
> end if
> end if
>
> else if(summe_nied / 10.0 .lt. 2.0 * (mean_temp + regen)) then
> klasse = "B" ; Trockenklima
>
> if(summe_nied / 10.0 .gt. mean_temp + regen) then
> klasse = klasse + "S" ; Steppenklima
> else
> klasse = klasse + "W" ; Wuestenklima
> end if
>
> if(mean_temp .gt. 18.0) then
> klasse = klasse + "h" ; heiß
> else if(mean_temp .le. 18.0) then
> klasse = klasse + "k" ; kalt
> end if
> end if
>
> else if((min(temp(:)) .gt. -3.0) .and. (min(temp(:)) .lt. 18.0) .and.
> (summe_nied / 10.0 .gt. 2.0 * (mean_temp + regen))) then
> klasse = "C" ; Warm-gemaeßigtes Klima
>
> if(min((/nied(mod(11 + hk, 12)),nied(mod(0 + hk, 12)),nied(mod(1 + hk=
,
> 12))/)) .lt. 0.1 * max((/nied(mod(5 + hk, 12)),nied(mod(6 + hk,
> 12)),nied(mod(7 + hk, 12))/))) then
>
>
> klasse = klasse + "w" ; warmes wintertrockenes Klima
>
>
>
> else if((min((/nied(mod(5 + hk, 12)),nied(mod(6 + hk, 12)),nied(mod(7
> + hk, 12))/)) .lt. 40.0) .and. (max((/nied(mod(11 + hk, 12)),nied(mod(0 +
> hk, 12)),nied(mod(1 + hk, 12))/) / 3.0) .gt. min((/nied(mod(5 + hk,
> 12)),nied(mod(6 + hk, 12)),nied(mod(7 + hk, 12))/)))) then
> klasse = klasse + "s" ; warmes sommertrockenes Klima
> else
> klasse = klasse + "f" ; feuchtgemaessigtes Klima
> end if
> end if
>
> if(max(temp(:)) .gt. 22.0) then
> klasse = klasse + "a" ;heisse Sommer
> else if((max(temp(:)) .lt. 22.0) .and. (warm .ge. 4.0)) then ;ab hier
> neu einfuegen
> klasse = klasse + "b" ; warme Sommer
> else if((max(temp(:)) .lt. 22.0) .and. (warm .le. 3.0) .and. (warm
> .ge. 1.0)) then
> klasse = klasse + "c" ; kuehle Sommer
> end if
> end if
> end if
>
> else if((min(temp(:)) .lt. -3.0) .and. (max(temp(:)) .gt. 10.0)) then
> klasse = "D" ; Boreales oder Schnee-Wald-Klima
>
> if(min((/nied(mod(11 + hk, 12)),nied(mod(0 + hk, 12)),nied(mod(1 + hk=
,
> 12))/)) .lt. 0.1 * max((/nied(mod(5 + hk, 12)),nied(mod(6 + hk,
> 12)),nied(mod(7 + hk, 12))/))) then
> klasse = klasse + "w" ; wintertrockenkaltes Klima
> else
> klasse = klasse + "f" ; winterfeuchtkaltes Klima
> end if
>
> if(max(temp(:)) .gt. 22.0) then
> klasse = klasse + "a" ; heiße Sommer
> else if((max(temp(:)) .lt. 22.0) .and. (warm .ge. 4.0)) then
> klasse = klasse + "b" ; warme Sommer
> else if((max(temp(:)) .lt. 22.0) .and. (warm .le. 3.0) .and. (warm
> .ge. 1.0)) then
> klasse = klasse + "c" ; kuehle Sommer
> else if(min(temp(:)) .lt. -38.0) then
> klasse = klasse + "d" ; strenge Winter
> end if
> end if
> end if
> end if
>
> else if(max(temp(:)) .lt. 10.0) then
> klasse = "E" ; Schneeklima
>
> if(max(temp(:)) .lt. 0.0) then
> klasse = klasse + "F" ; Tundrenklima
> else
> klasse = klasse + "T" ; Klima ewigen Frostes
> end if
>
> end if
> end if
> end if
> end if
> end if
>
> if (klasse .eq. "Af") then
> Klimaklasse(0,loc_lat,loc_lon) = 30.1
> else if (klasse .eq. "Am") then
> Klimaklasse(0,loc_lat,loc_lon) = 29.1
> else if (klasse .eq. "As") then
> Klimaklasse(0,loc_lat,loc_lon) = 28.1
> else if (klasse .eq. "Aw") then
> Klimaklasse(0,loc_lat,loc_lon) = 27.1
> else if (klasse .eq. "BWk") then
> Klimaklasse(0,loc_lat,loc_lon) = 26.1
> else if (klasse .eq. "BWh") then
> Klimaklasse(0,loc_lat,loc_lon) = 25.1
> else if (klasse .eq. "BSk") then
> Klimaklasse(0,loc_lat,loc_lon) = 24.1
> else if (klasse .eq. "BSh") then
> Klimaklasse(0,loc_lat,loc_lon) = 23.1
> else if (klasse .eq. "Cfa") then
> Klimaklasse(0,loc_lat,loc_lon) = 22.1
> else if (klasse .eq. "Cfb") then
> Klimaklasse(0,loc_lat,loc_lon) = 21.1
> else if (klasse .eq. "Cfc") then
> Klimaklasse(0,loc_lat,loc_lon) = 20.1
> else if (klasse .eq. "Csa") then
> Klimaklasse(0,loc_lat,loc_lon) = 19.1
> else if (klasse .eq. "Csb") then
> Klimaklasse(0,loc_lat,loc_lon) = 18.1
> else if (klasse .eq. "Csc") then
> Klimaklasse(0,loc_lat,loc_lon) = 17.1
> else if (klasse .eq. "Cwa") then
> Klimaklasse(0,loc_lat,loc_lon) = 16.1
> else if (klasse .eq. "Cwb") then
> Klimaklasse(0,loc_lat,loc_lon) = 15.1
> else if (klasse .eq. "Cwc") then
> Klimaklasse(0,loc_lat,loc_lon) = 14.1
> else if (klasse .eq. "Dfa") then
> Klimaklasse(0,loc_lat,loc_lon) = 13.1
> else if (klasse .eq. "Dfb") then
> Klimaklasse(0,loc_lat,loc_lon) = 12.1
> else if (klasse .eq. "Dfc") then
> Klimaklasse(0,loc_lat,loc_lon) = 11.1
> else if (klasse .eq. "Dfd") then
> Klimaklasse(0,loc_lat,loc_lon) = 10.1
> else if (klasse .eq. "Dsa") then
> Klimaklasse(0,loc_lat,loc_lon) = 9.1
> else if (klasse .eq. "Dsb") then
> Klimaklasse(0,loc_lat,loc_lon) = 8.1
> else if (klasse .eq. "Dsc") then
> Klimaklasse(0,loc_lat,loc_lon) = 7.1
> else if (klasse .eq. "Dsd") then
> Klimaklasse(0,loc_lat,loc_lon) = 6.1
> else if (klasse .eq. "Dwa") then
> Klimaklasse(0,loc_lat,loc_lon) = 5.1
> else if (klasse .eq. "Dwb") then
> Klimaklasse(0,loc_lat,loc_lon) = 4.1
> else if (klasse .eq. "Dwc") then
> Klimaklasse(0,loc_lat,loc_lon) = 3.1
> else if (klasse .eq. "Dwd") then
> Klimaklasse(0,loc_lat,loc_lon) = 2.1
> else if (klasse .eq. "EF") then
> Klimaklasse(0,loc_lat,loc_lon) = 1.1
> else if (klasse .eq. "ET") then
> Klimaklasse(0,loc_lat,loc_lon) = 0.1
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
> end if
>
>
> end do
> end do
>
> wks = gsn_open_wks("ps","Koeppen") ; open a ps file
> ;plot = new(4,graphic) ; create a plot array
>
> gsn_define_colormap(wks,"cb_9step")
>
> res = True
> res@gsnMaximize = True
> res@gsnAddCyclic = False
> res@mpDataBaseVersion = "mediumres"
> res@mpOutlineBoundarySets = "national"
>
> res@mpCenterLonF = rotpole@grid_north_pole_longitude
> res@mpCenterLatF = rotpole@grid_north_pole_latitude -90
> res@cnFillOn = True
> res@cnLinesOn = False
> res@cnLineLabelsOn = False
>
> res@cnExplicitLegendLabelsOn = True
> ;res@lbLabelAlignment = "BoxCenters" ;ab hier neu einfuegen
>
> res@cnLevelSelectionMode = "ExplicitLevels"
> res@cnLevels = ispan(16,27,1)
> res@lbLabelStrings = (/"ET", "EF", "Dwd", "Dwc", "Dwb", "Dwa", "Dsd",
> "Dsc", "Dsb", "Dsa", "Dfd", "Dfc", "Dfb", "Dfa", "Cwc", "Cwb", "Cwa",
> "Csc", "Csb", "Csa", "Cfc", "Cfb", "Cfa", "BSh", "BSk", "BWh", "BWk", "Aw=
",
> "As", "Am", "Af"/)
> res@cnMinLevelValF = min(Klimaklasse(0,:,:)) - 0.5
> res@cnMaxLevelValF = max(Klimaklasse(0,:,:)) + 0.5
> res@cnLevelSpacingF = 1.0
>
> res@pmTickMarkDisplayMode = "Always"; use NCL default lat/lon labels
> res@mpMinLatF = -23
> res@mpMaxLatF = -10
> res@mpMinLonF = -72
> res@mpMaxLonF = -55
> res@mpLimitMode = "latlon"
> res@mpLeftCornerLatF = lat2d(0,0)
> res@mpLeftCornerLonF = lon2d(0,0)
> res@mpRightCornerLatF = lat2d(latsize-1,lonsize-1)
> res@mpRightCornerLonF = lon2d(latsize-1,lonsize-1)
>
> res@gsnCenterString = "Klimaklassifikation"
> plot = gsn_csm_contour_map(wks,Klimaklasse(0,:,:),res)
>
> ;res@gsnCenterString = "Mitteltemperatur 1991-2020"
> ;plot(1) = gsn_csm_contour_map(wks,N,res)
>
> ;res@gsnCenterString = "Mitteltemperatur 2021-2050"
> ;plot(2) = gsn_csm_contour_map(wks,O,res)
>
> ;res@gsnCenterString = "Mitteltemperatur 2051-2080"
> ;plot(3) = gsn_csm_contour_map(wks,P,res)
>
> ;gsn_panel(wks,plot,(/2,2/),False) ; now draw as one plot
>
> end
>
> I get the following error message:
>
>
>
> fatal:Subscript out of range, error in subscript #3
> fatal:An error occurred reading Temperatur
> fatal:An error occurred reading Temperatur
> fatal:["Execute.c":8126]:Execute: Error occurred at or near line 43
>
>
>
> Line 43 is the following line: temp = Temperatur(:,0,loc_lat,loc_lon)
>
>
>
> How can I solve this problem? I hope you can help me.
>
>
>
> Attachment:
>
>
>
> - file:
> MPIM-REMO_EH5OMA1Br3_HR_044_BOL_1961-1990_167_minus_1_hour_Celsiustt_ymon=
mean.nc
>
> - file:
> MPIM-REMO_EH5OMA1Br3_HR_044_BOL_1961-1990_142_143_minus_1_hour_ymonsum_me=
an_30_years.nc
>
> - file: Koeppen_Geiger_neu.ncl (my NCL code)
>
> - file: ll.out
>
>
>
>
>
> Best regards
>
>
>
> Stephan
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:http://mailman.ucar.ed=
u/mailman/listinfo/ncl-talk
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
Received on Mon Jun 23 02:09:05 2014
This archive was generated by hypermail 2.1.8 : Wed Jul 23 2014 - 15:33:46 MDT