[NARCCAP-discuss] gfdl timeslice pressure coordinates

Sun Apr 15 13:16:05 MDT 2012

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