Re: Error allocation

From: Mary Haley <haley_at_nyahnyahspammersnyahnyah>
Date: Thu Jul 12 2012 - 11:23:01 MDT

Cecille,

This definitely looks like a memory issue. To help free up memory so you can create more variables later in your script, you should delete variables as soon as you no longer need them.

For example, after you call the "average" procedure, you can delete several of those variables with:

   delete([/suma,sumb,sumc,sumd,amor,.../])

You have some additional arrays that I think you don't need. For example, you have both mavg_wsm6 and avg_wsm6. avg_wsm6 is just mavg_wsm6 multiplied by a million, correct? If so, you can get rid of avg_wsm6 and just do:

    mavg_wsm6 = mavg_wsm6* 1000000

Of course, you need to replace all references to "avg_wsm6" to "mavg_wsm6(i,:)". Ditto for avg_tho, avg_mil, and avg_mor. I think you can just use the "m" version of these arrays.

Here are some other improvements you can make. For one, you don't need to do:

  awsm6 = (d->MELTOUT3D(:,:,:,:))
  tho = (b->MELTOUT3D(:,:,:,:))
  mil = (c->MELTOUT3D(:,:,:,:))
  amor = (a->MELTOUT3D(:,:,:,:))
  ph = (a->PH(:,:,:,:))/9.81 ;eliminate geopotential
  phb = (a->PHB(:,:,:,:))/9.81

This can cause NCL to go into unnecessary subscripting mode. Just use this:

  awsm6 = d->MELTOUT3D
  tho = b->MELTOUT3D
  mil = c->MELTOUT3D
  amor = a->MELTOUT3D
  ph = (a->PH)/9.81 ;eliminate geopotential
  phb = (a->PHB)/9.81

Here, it looks like you are creating a set of numbers from 0 to num_vert in steps of 1. Instead of:

  heightplot = new(num_vert, float)

  heightplot(0) = 0
  do j=1,num_vert-1
    heightplot(j) = heightplot(j-1)+1
  end do

you can use the "ispan" function:

  heightplot = ispan(0,num_vert-1,1)

Whenver you put things in a do loop, this can really slow things down. You might be able to improve a little on memory,
and a lot on making your code run faster if you can get rid of do loops, or move code outside of do loops.

For example, here's a double do loop where you can use "where" instead:

  do i = 0, dimsizes(mavg_wsm6(:,0))-1
    do k = 0, dimsizes(mavg_wsm6(0,:))-1
      if (mavg_wsm6(i,k) .lt. 0.) then
        mavg_wsm6(i,k) = 0.
      end if
    end do
  end do

Replacement code:

  mavg_wsm6 = where(mavg_wsm6.lt.0,0, mavg_wsm6)

or even simpler:

  mavg_wsm6 = mavg_wsm6 > 0

You can use NC's array arithmetic to get rid of other do loops. For example, instead of:

  do k=0,num_vert-1
    avg_wsm6(k) = mavg_wsm6(i,k)* 1000*1000
    avg_tho(k) = mavg_tho(i,k)* 1000*1000
    avg_mil(k) = mavg_mil(i,k)* 1000*1000
    avg_mor(k) = mavg_mor(i,k)* 1000*1000
    allq(0,k) = avg_mor(k)
    allq(1,k) = avg_tho(k)
    allq(2,k) = avg_mil(k)
    allq(3,k) = avg_wsm6(k)
  end do

use:

    avg_wsm6(k) = mavg_wsm6(i,:)* 1000000
    avg_tho(k) = mavg_tho(i,:)* 1000000
    avg_mil(k) = mavg_mil(i,:)* 1000000
    avg_mor(k) = mavg_mor(i,:)* 1000000
    allq(0,:) = avg_mor
    allq(1,:) = avg_tho
    allq(2,:) = avg_mil
    allq(3,:) = avg_wsm6

However, see my note above about not needing both avg_wsm6 and mavg_wsm6.

You have a do loop across "num_times", and then you have code inside this loop that shouldn't be there. For example:

  do i = 1,num_times -1
    . . .
  leg_label = (/"MOR","THOMP","MILB","WSM6"/)

  res = True

  res@tiMainString = timelabl(i)
  res@tiXAxisString = type
  res@tiYAxisString = "Height (km)"
  res@tiXAxisFontHeightF = 0.02
  res@tiYAxisFontHeightF = 0.02

  res@pmLegendDisplayMode = "Always"
  res@pmLegendZone = 0
  res@pmLegendSide = "Top"
  res@pmLegendWidthF = 0.15
  res@pmLegendHeightF = 0.15
  res@lgLabelFontHeightF = .40
  res@xyExplicitLegendLabels = leg_label
  res@pmLegendParallelPosF = 0.32
  res@pmLegendOrthogonalPosF = 0.32
  . . .
  end do

With the exception of:

  res@tiMainString = timelabl(i)

the rest of this code should appear above the line:

  do i = 1,num_times -1

Please try these suggestions and let me know if this doesn't help.

--Mary

On Jul 12, 2012, at 9:55 AM, Cecille M. Villanueva-Birriel wrote:

> ;This program will produce vertical profiles of total
> ;hydrometeor species at a given time. This plot can be used
> ;to see relationships such as the melting of graupel / freezing
> ;of rain couplet near the freezing level.
> ;
> ;All previous comments consolidated with general script cleanup.
> ;*****************************************************
>
> ;Standard loading scripts for functions and plotting?
> 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 "./contributed.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
> external AVERAGE "./average_rates_meltout.so"
>
> begin
>
> ;;EDITABLE CODE BEGINS HERE
>
> ; RUN_NUM = "Hour-1"
> ; WRF_TIME = 8 ; This is the time index, i.e...index *dt = time.
> ; timelabl = "0110"
> ; mp_label = "GCE"
> ; dif_labl = "WSM6gro2gacs - WSM6defgro2"
> ; di2_labl = "WSM6gro2gacs-WSM6defgro2"
> dif_labl = "WSM6cv"
> di2_labl = "WSM6cv"
> MULTI_FILE_FLAG = 0 ; 1 for yes, 0 for no.
> RESTART_NEWTIME_INDEX = 24 ; i.e. dt=5 & restart time=120 min,
> ; 120/5=24+1(time zero)-1(duplicate time)->Index=24
>
>
> type = "Ice/Snow/Gaupel/Hail Melting (mg kg-1 s-1)"
> pdftype = "meltout" ;name of the pdf file
>
> timelabl = (/"0000","0005","0010","0015","0020","0025","0030","0035","0040","0045",\
> "0050","0055","0100","0105","0110","0115","0120","0125","0130"/)
>
> ; timelabl = (/"0000","0030","0100"/)
> ; WRF_TIME = ispan(0,7,1)
>
> ;;END EDITABLE CODE
>
> a = addfile("../wrfout_d01_0001-01-01_00:00:00_MOR_250m.nc", "r")
> b = addfile("../wrfout_d01_0001-01-01_00:00:00_THO_250m.nc", "r")
> c = addfile("../wrfout_d01_0001-01-01_00:00:00_MY_250m.nc", "r")
> d = addfile("../wrfout_d01_0001-01-01_00:00:00_WSM6_250m.nc", "r")
> ;f = addfile("wrfout_d01_0001-01-01_00:00:00_nssl.nc", "r")
> ; a = addfile("wrfout-wsm6cvfut.nc", "r")
> ; b = addfile("wrfout-wsm6defgro2.nc", "r")
>
> awsm6 = (d->MELTOUT3D(:,:,:,:))
> tho = (b->MELTOUT3D(:,:,:,:))
> mil = (c->MELTOUT3D(:,:,:,:))
> amor = (a->MELTOUT3D(:,:,:,:))
>
> ph = (a->PH(:,:,:,:))/9.81 ;eliminate geopotential
> phb = (a->PHB(:,:,:,:))/9.81
> height = (ph + phb)/1000
> ; height = wrf_user_getvar(a,"z",-1)/1000
>
> delete(ph)
> delete(phb)
>
> ; times = (a->Times(:,:))
> times = wrf_user_list_times(a)
>
> delete(a)
> ; delete(b)
>
> ;I will need the number of vertical/horiz levels for the height labels later
>
> num_times = dimsizes(times)
> num_vert = dimsizes(tho(0,:,1,1))
> num_horiza = dimsizes(tho(0,0,:,1))
> num_horizb = dimsizes(tho(0,0,1,:))
>
> nt = num_times
> nz = num_vert
> ny = num_horiza
> nx = num_horizb
>
> print (num_horiza)
> print (num_horizb)
>
>
> timid = ispan(0,18,1)
>
> ;The following makes an array with the model height (in km)
> ;to be used for the plot labels, in actuallity, we need to make a
> ;code that plots intergers vs. qmass
>
> tmp1 = floattoint((num_vert/5)+1)
>
> height_array = new(tmp1, float)
>
> height_array(0) = avg(height(:,0,:,:)) ;because (j*5)-1 for 0 gives a negative index
>
> do j=1,tmp1-1
> tmp2 = (j*5)-1
> height_array(j) = floattoint(avg(height(:,tmp2,:,:)))
> end do
>
> ;Now generate the intergers for plotting
>
> heightplot = new(num_vert, float)
>
> heightplot(0) = 0
> do j=1,num_vert-1
> heightplot(j) = heightplot(j-1)+1
> end do
>
> ;*****************************************************
> ;These are all of my array declarations, which will
> ;be the length of the num_vert variable from above
> ;*****************************************************
>
> ; qavg_vapor = new(num_vert, float)
> avg_wsm6 = new(num_vert, float)
> avg_tho = new(num_vert, float)
> avg_mil = new(num_vert, float)
> avg_mor = new(num_vert, float)
> ; avg_nssl = new(num_vert, float)
>
> ; qavg_vapob = new(num_vert, float)
> ; qavg_cloub = new(num_vert, float)
> ; qavg_raib = new(num_vert, float)
> ; qavg_icb = new(num_vert, float)
> ; qavg_snob = new(num_vert, float)
> ; qavg_graub = new(num_vert, float)
>
> ; qdif_vapor = new(num_vert, float)
> ; qdif_cloud = new(num_vert, float)
> ; qdif_rain = new(num_vert, float)
> ; qdif_ice = new(num_vert, float)
> ; qdif_snow = new(num_vert, float)
> ; qdif_graup = new(num_vert, float)
>
> qds_vapor = new((/num_times/), float)
> qds_cloud = new((/num_times/), float)
> qds_rain = new((/num_times/), float)
> qds_ice = new((/num_times/), float)
> qds_snow = new((/num_times/), float)
> qds_graup = new((/num_times/), float)
>
> mavg_wsm6 = new((/num_times,num_vert/), float)
> mavg_tho = new((/num_times,num_vert/), float)
> mavg_mil = new((/num_times,num_vert/), float)
> mavg_mor = new((/num_times,num_vert/), float)
> ; mavg_nssl = new((/num_times,num_vert/), float)
>
>
> suma = new((/num_times,num_vert/), float)
> sumb = new((/num_times,num_vert/), float)
> sumc = new((/num_times,num_vert/), float)
> sumd = new((/num_times,num_vert/), float)
> sume = new((/num_times,num_vert/), float)
>
> suma = 0.
> sumb = 0.
> sumc = 0.
> sumd = 0.
> sume = 0.
>
> li = new((/num_times,num_vert/), float)
> th = new((/num_times,num_vert/), float)
> mi = new((/num_times,num_vert/), float)
> mo = new((/num_times,num_vert/), float)
> ms = new((/num_times,num_vert/), float)
>
> li = 0.
> th = 0.
> mi = 0.
> mo = 0.
> ms = 0.
>
>
> ; allq = new((/6, num_vert/), float)
> allq = new((/4, num_vert/), float)
> ;*****************************************************
> ;This code will obtain the maximum values of some
> ;hydrometeor anywhere in the 3D domain for a fixed time.
> ;This may need to be updated to correct for domains
> ;that don't include the cells which develop on the
> ;outflow of my supercell. NOTE: units now in g/kg
> ;*****************************************************
>
> ; suma = dim_sum_n(mor,(/2,3/))
> ; sumb = dim_sum_n(tho,(/2,3/))
> ; sumc = dim_sum_n(mil,(/2,3/))
> ; sumd = dim_sum_n(wsm6,(/2,3/))
> ;print (suma*1000)
>
> ; li = dim_num_n(wsm6.ne.0,(/2,3/))
> ; th = dim_num_n(tho.ne.0,(/2,3/))
> ; mi = dim_num_n(mi.ne.0,(/2,3/))
> ; mo = dim_num_n(mo.ne.0,(/2,3/))
>
>
> AVERAGE::average(suma,sumb,sumc,sumd,amor,tho,mil,awsm6,li,th,mi,mo,mavg_wsm6,mavg_tho,mavg_mil,mavg_mor,nx,ny,nz,nt)
>
> print (suma*1000)
> print (mo)
>
>
>
> ;mavg_ice = smth9(mavg_ice, 0.50, 0.25, False)
> do i = 0, dimsizes(mavg_wsm6(:,0))-1
> do k = 0, dimsizes(mavg_wsm6(0,:))-1
> if (mavg_wsm6(i,k) .lt. 0.) then
> mavg_wsm6(i,k) = 0.
> end if
> end do
> end do
> do i = 0, dimsizes(mavg_tho(:,0))-1
> do k = 0, dimsizes(mavg_tho(0,:))-1
> if (mavg_tho(i,k) .lt. 0.) then
> mavg_tho(i,k) = 0.
> end if
> end do
> end do
> do i = 0, dimsizes(mavg_mil(:,0))-1
> do k = 0, dimsizes(mavg_mil(0,:))-1
> if (mavg_mil(i,k) .lt. 0.) then
> mavg_mil(i,k) = 0.
> end if
> end do
> end do
> do i = 0, dimsizes(mavg_mor(:,0))-1
> do k = 0, dimsizes(mavg_mor(0,:))-1
> if (mavg_mor(i,k) .lt. 0.) then
> mavg_mor(i,k) = 0.
> end if
> end do
> end do
> ; do i = 0, dimsizes(mavg_nssl(:,0))-1
> ; do k = 0, dimsizes(mavg_nssl(0,:))-1
> ; if (mavg_nssl(i,k) .lt. 0.) then
> ; mavg_nssl(i,k) = 0.
> ; end if
> ; end do
> ; end do
>
> print (mavg_mor)
>
> ;Open the workstation
> wks = gsn_open_wks ("pdf", pdftype+"_scheme_avg_height")
>
> do i = 1,num_times -1
>
> ;fub = i
>
> do k=0,num_vert-1
> ; avg_vapor(k) = (sum(qvapor(fub,k,:,:))) * 1000
> avg_wsm6(k) = mavg_wsm6(i,k)* 1000*1000
> avg_tho(k) = mavg_tho(i,k)* 1000*1000
> avg_mil(k) = mavg_mil(i,k)* 1000*1000
> avg_mor(k) = mavg_mor(i,k)* 1000*1000
> ; avg_nssl(k) = mavg_nssl(i,k)* 1000
> ; qavg_vapob(k) = (sum(qvaporb(fub,k,:,:))) * 1000
> ; qavg_cloub(k) = (sum(qcloudb(fub,k,:,:))) * 1000
> ; qavg_raib(k) = (sum(qrainb(fub,k,:,:))) * 1000
> ; qavg_icb(k) = (sum(qiceb(fub,k,:,:))) * 1000
> ; qavg_snob(k) = (sum(qsnowb(fub,k,:,:))) * 1000
> ; qavg_graub(k) = (sum(qgraupb(fub,k,:,:))) * 1000
> ; qdif_vapor(k) = qavg_vapor(k) - qavg_vapob(k)
> ; qdif_cloud(k) = qavg_cloud(k) - qavg_cloub(k)
> ; qdif_rain(k) = qavg_rain(k) - qavg_raib(k)
> ; qdif_ice(k) = qavg_ice(k) - qavg_icb(k)
> ; qdif_snow(k) = qavg_snow(k) - qavg_snob(k)
> ; qdif_graup(k) = qavg_graup(k) - qavg_graub(k)
>
> ;print (avg_mor)
> ;allq(0,k) = avg_vapor(k)
> allq(0,k) = avg_mor(k)
> allq(1,k) = avg_tho(k)
> allq(2,k) = avg_mil(k)
> allq(3,k) = avg_wsm6(k)
> ; allq(4,i,k) = avg_nssl(k)
> end do
>
> ; qds_vapor(i) = sum(avg_vapor)
> ;qds_cloud(i) = sum(avg_cloud)
> ;qds_rain(i) = sum(avg_rain)
> ;qds_ice(i) = sum(avg_ice)
> ;qds_snow(i) = sum(avg_snow)
> ;qds_graup(i) = sum(avg_graup)
>
> ;*****************************************************
> ;Time to make some plots, yo
> ;*****************************************************
>
> ;Open the workstation
> ; wks = gsn_open_wks ("pdf", "Qdiftot_vALT"\
> ; +di2_labl+"")
>
>
> ; allq(0,:) = qavg_vapor(i)
>
> ; leg_label = (/"QVapor","QCloud","QRain","QIce","QSnow","QGraupel"/)
> leg_label = (/"MOR","THOMP","MILB","WSM6"/)
>
> res = True
>
> res@tiMainString = timelabl(i)
> res@tiXAxisString = type
> res@tiYAxisString = "Height (km)"
> res@tiXAxisFontHeightF = 0.02
> res@tiYAxisFontHeightF = 0.02
>
> res@pmLegendDisplayMode = "Always"
> res@pmLegendZone = 0
> res@pmLegendSide = "Top"
> res@pmLegendWidthF = 0.15
> res@pmLegendHeightF = 0.15
> res@lgLabelFontHeightF = .40
> res@xyExplicitLegendLabels = leg_label
> res@pmLegendParallelPosF = 0.32
> res@pmLegendOrthogonalPosF = 0.32
>
> res@tmYLMode = "Explicit"
> res@tmYLMinorOn = True
> res@trYMaxF = num_vert-1
> res@tmYLTickSpacingF = 5
> res@tmYLMinorPerMajor = 4
> res@tmYLLabelFontHeightF = 0.015
> res@tmYLValues = ispan(0,num_vert-1,5)
> res@tmYLMinorValues = ispan(0,num_vert-1,1)
> res@tmYLLabels = height_array
> res@tmYLPrecision = 1
>
> ; res@xyMarkLineMode = "MarkLines"
> ; res@xyMarkers = (/9,4/)
> ;res@xyLineThicknessF = 2.0
> res@xyLineThicknesses = (/2.0,3.0,2.0,3.0/)
> ;res@xyLineColors = (/"black","green","blue","brown","orange","purple"/)
> res@xyLineColors = (/"black","gray","black","gray55"/)
>
> ;res@TimeLabel = times(i)
>
> plot = gsn_xy(wks,allq,heightplot,res)
>
>
> end do
>
> ; grandtotnovap = qds_cloud+qds_rain+qds_ice+qds_snow+qds_graup
>
> ; print (avg_mor)
>
> ; print(qds_vapor)
>
> ; print(qds_cloud)
> ; print(qds_rain)
> ; print(qds_ice)
> ; print(qds_snow)
> ; print(qds_graup)
> ; print(grandtotnovap)
>
> end
>
> --
> ********************************************************
> Cecille M. Villanueva Birriel
> Graduate Research Assistant
> Cloud Microphysics Research Group
> Purdue University
> Department of Earth and Atmospheric Sciences
> 550 Stadium Mall Drive, West Lafayette, IN 47907-2051
> email: cvillanu@purdue.edu
> ********************************************************
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Jul 12 11:23:12 2012

This archive was generated by hypermail 2.1.8 : Wed Jul 18 2012 - 14:33:00 MDT