[NARCCAP-discuss] gfdl timeslice pressure coordinates

Jeanne Thibeault jeanne.thibeault at uconn.edu
Sun Apr 15 20:01:09 MDT 2012 

Hi Kelly,
Thanks for sharing this.  I appreciate the help.
Jeanne
On Apr 15, 2012, at 3:16 PM, Kelly Mahoney wrote:

> Hi Jeanne,
> 
> My colleague, Jamie Scott, developed a code to convert GFDL timeslice 
> experiment data to pressure coordinates.
> 
> It will need some adapting, as we were using it with files that we had 
> subset-ed from the larger GFDL data, but I've copied it below just in 
> case it is helpful to you at all.
> 
> Good luck,
> Kelly Mahoney
> 
> 
> 
> 
>> 
>> 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
>> setvalues NhlGetWorkspaceObjectId()
>>        "wsMaximumSize" : 100000000000
>> end setvalues
>> ;input directory
>> diri="/nas/kmahoney/WRF_ICs_From_NARCCAP/ICs/GFDL/Timeslices/Future/Rank_10_20610812/"
>> 
>> ;these variables are used to loop through and build the dates that are 
>> in the input file names.
>> ;each array position in years,months,days,hours creates one date to 
>> read.  Loop through all dates.
>> 
>> years=(/"2061","2061","2061","2061","2061","2061","2061","2061","2061"/)  
>> ;;input years to process
>> months=(/"08","08","08","08","08","08","08","08","08"/)       ;;input 
>> months
>> days=(/"12","12","13","13","13","13","14","14","14"/)         ;;input 
>> days of month
>> hours=(/"12","18","00","06","12","18","00","06","12"/)        ;;input 
>> hours
>> 
>> 
>> do it=0,dimsizes(years)-1
>> yyyy=years(it)
>> mm=months(it)
>> dd=days(it)
>> hh=hours(it)
>> 
>> inname=diri+"ua_gfdl_"+yyyy+mm+dd+hh+".nc"
>> infile=addfile(inname,"r")
>> 
>> inname2="/Projects/IPCC/NARCCAP/GFDL/Timeslices/orog_gfdl.nc"
>> infile2=addfile(inname2,"r")
>> 
>> ;inname3=diri+"hus_gfdl_1986071706.nc"
>> ;infile3=addfile(inname3,"r")
>> 
>> if (it.eq.0) then
>> a_bnds=infile->a_bnds
>> b_bnds=infile->b_bnds
>> a=dim_avg(a_bnds) ;;a @ mid point of level
>> a=a/100000.  ;; normalize a
>> b=dim_avg(b_bnds) ;;b @ mid point of level
>> ;store a & b and level interfaces
>> hyai=new((/dimsizes(a)+1/),typeof(a))
>> hybi=new((/dimsizes(b)+1/),typeof(b))
>> hyai(0:dimsizes(a)-1)=a_bnds(:,0)
>> hyai(dimsizes(a))=0.
>> hybi(0:dimsizes(b)-1)=b_bnds(:,0)
>> hybi(dimsizes(b))=1.
>> hyai=hyai/100000. ;; normalize by P0
>> 
>> print(hyai+" "+hybi)
>> plevel = 
>> (/1000.,950.,900.,850.,800.,750.,700.,650.,600.,550.,500.,450.,400.,350.,300.,250.,200.,150.,100.,50./)
>> plevel!0="level"
>> plevel at units = "millibar" ;
>> plevel at long_name = "Level" ;
>> plevel at positive = "down" ;
>> plevel&level=plevel
>> end if ;do once
>> time=infile->time  ;;
>> print("time coord value: "+time)
>> 
>> ua=infile->ua
>> ;hus=infile3->hus
>> ;tv=ta*(1.0+0.61*hus)
>> lat=infile->lat
>> lon=infile->lon
>> lat1=min(lat)
>> lat2=max(lat)
>> lon1=min(lon)
>> lon2=max(lon)
>> if (it.eq.0)
>> z=infile2->orog({lat1:lat2},{lon1:lon2})
>> gz=z*9.81 ;;geopotential hgt @ sfc
>> end if ;; do once
>> tbot=ua(:,{996.1084},:,:)
>> ps=infile->ps
>> 
>> 
>> print("interpolating")
>> tplev=vinth2p_ecmwf(ua,a,b,plevel,ps,2,1000.,1,True,0,tbot,gz)
>> 
>> outname=diri+"ua.plev."+yyyy+mm+dd+hh+".nc"
>> system("rm "+outname)
>> outfile=addfile(outname,"c")
>> ;;;output netcdf efficiently
>> print("reading globatts")
>> globatts=getvaratts(infile) ; get global attributes
>> ;print(dimsizes(globatts))
>> tplevatts=True    ;transfer attributes from input file to output file
>> do i=0,dimsizes(globatts)-1
>> ; print(i)
>> ; print(file1@$globatts(i)$)
>> tplevatts@$globatts(i)$=infile@$globatts(i)$
>> end do
>> 
>> print("defining dims")
>> ;tplev at coordinates = "level lon lat time"
>> dimNames=(/"time","level","lat","lon"/)  ;predefine coordinate vars
>> dimSizes=(/-1,dimsizes(plevel),dimsizes(lat),dimsizes(lon)/)
>> dimUnlim=(/True,False,False,False/)
>> filedimdef(outfile,dimNames,dimSizes,dimUnlim)
>> 
>> print("defining coords")
>> filevardef(outfile,"time","double","time") ;predefine var names,type,dims
>> filevardef(outfile,"level","float","level")
>> filevardef(outfile,"lat","double","lat")
>> filevardef(outfile,"lon","double","lon")
>> print("defining tplev")
>> filevardef(outfile,"ua","float",(/"time","level","lat","lon"/))
>> print("writing varatts")
>> 
>> filevarattdef(outfile,"ua",ua)         ;copy variable atts to file
>> filevarattdef(outfile,"lon",lon)
>> filevarattdef(outfile,"lat",lat)
>> filevarattdef(outfile,"level",plevel)
>> filevarattdef(outfile,"time",time)
>> 
>> fileattdef(outfile,tplevatts)     ;create glabal atts
>> ;output data
>> print("writing coord data")
>> outfile->time=(/time/)
>> outfile->level=(/plevel/)
>> outfile->lat =(/lat/)
>> outfile->lon =(/lon/)
>> 
>> 
>> 
>> print("writing netcdf")
>> outfile->ua=(/tplev/)
>> delete(outfile)
>> 
>> delete(ps)
>> delete(ua)
>> delete(tplev)
>> 
>> end do ;;loop through time
>> 
>> end
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> On 4/13/2012 6:33 PM, Jeanne Thibeault wrote:
>> Thank you, Seth.
>> Jeanne
>> On Apr 13, 2012, at 2:08 PM, Seth McGinnis wrote:
>> 
>>> Hi Jeanne,
>>> 
>>> The NARCCAP modelers at GFDL are Isaac Held and Bruce Wyman.  You can find
>>> contact info for them (and all the modelers) in the "Participants Directory" on
>>> the NARCCAP website: http://www.narccap.ucar.edu/about/narccap-directory.html
>>> 
>>> Cheers,
>>> 
>>> --Seth
>>> 
>>> 
>>> On Thu, 12 Apr 2012 16:13:30 -0400
>>> "Thibeault, Jeanne"<jeanne.thibeault at uconn.edu>  wrote:
>>>> Hello NARCCAP users,
>>>> 
>>>> I am analyzing the GFDL timeslice experiments and need to convert the pressure
>>>> coordinates from hybrid sigma coordinates to pressure levels (hPa).  I have a
>>>> fortran code that I downloaded from GFDL
>>>> (http://www.gfdl.noaa.gov/narccap-am2--data, see attachment), but I received
>>>> an error message trying to run the program.
>>>> 
>>>> NETCDF ERROR
>>>> Invalid argument
>>>> STOP 111
>>>> 
>>>> The fortran code seems to fail at this line:
>>>> istat = NF90_OPEN (trim(file_in), NF90_NOWRITE, ncid_in)
>>>> if (istat /= NF90_NOERR) call error_handler (ncode=istat)
>>>> 
>>>> Does anyone know if there is a contact at GFDL that can help with this?  Is
>>>> there another way I can interpolate the data to pressure levels?
>>>> 
>>>> Thank you for your help,
>>>> Jeanne Thibeault
>>> _______________________________________________
>>> narccap-discuss mailing list
>>> narccap-discuss at mailman.ucar.edu
>>> http://mailman.ucar.edu/mailman/listinfo/narccap-discuss
>> _______________________________________________
>> narccap-discuss mailing list
>> narccap-discuss at mailman.ucar.edu
>> http://mailman.ucar.edu/mailman/listinfo/narccap-discuss
> _______________________________________________
> narccap-discuss mailing list
> narccap-discuss at mailman.ucar.edu
> http://mailman.ucar.edu/mailman/listinfo/narccap-discuss

More information about the narccap-discuss mailing list