;:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: ; ; This file is part of ICTP RegCM. ; ; ICTP RegCM is free software: you can redistribute it and/or modify ; it under the terms of the GNU General Public License as published by ; the Free Software Foundation, either version 3 of the License, or ; (at your option) any later version. ; ; ICTP RegCM is distributed in the hope that it will be useful, ; but WITHOUT ANY WARRANTY; without even the implied warranty of ; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ; GNU General Public License for more details. ; ; You should have received a copy of the GNU General Public License ; along with ICTP RegCM. If not, see . ; ;:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 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/contrib/ut_string.ncl" begin ;************************* ; read in data ;************************* ; a = addfile("SUDD_SRF.annual_eman.nc","r") a = addfile("SUDD_SRF.annual_grell_nnrp2.nc","r") ; b = addfile("SUDD_SRF.annual_grell_nnrp2_land_mod.nc","r") b = addfile("SUDD_SRF.annual_100percent.nc","r") ; prec = a->tpr ; prec1 = b->tpr prec = a->t2m prec1 = b->t2m landmask = a->mask ;print(prec(:,{6},{31})) xlat = b->xlat(0,:,:) xlon = b->xlon(0,:,:) dims = dimsizes(xlat) nlat = dims(1) nlon = dims(0) ; pdims = dimsizes(prec) ; ntimes = pdims(0) ; times = a->time; ;*************************************** ; For demo purposes: create a companion TIME array use thus to replace the ; current time coordinate variable, then use this to access specific months ;*************************************** ; time = a->time ; read in current array ; ntime = dimsizes(time) ; number of months ; TIME = new (ntime,integer) ; delete(TIME@_FillValue) ; TIME@long_name = "TIME" ; TIME@units = "yyyymm" ; TIME!0 = "TIME" ; nmo = 0 ; nyear = 1998 ; do nt = 0,ntime-1 ; nmo = nmo+1 ; TIME(nt) = nyear*100 + nmo ; if (nmo.eq.12) then ; nmo = 0 ; nyear = nyear+1 ; end if ; end do ; TIME&TIME = TIME ; coordinate variable ; delete (prec&time) ; delete the currently associated coordinate variable ; prec&time = TIME ; associate the new coordinate variable ;*************************************************** ; indices for special dates ;********************************* ; ind01 = ind(TIME.eq.199801) ; Jan 1996 ; ind02 = ind(TIME.eq.200412) ; Dec 2004 ;print(ind01+" "+ind02+" ") ;******************************** ; over specified periods ;******************************** ; varave = clmMonTLL ( prec(ind01:ind02,:,:) ) ; varave2 = clmMonTLL ( var2(ind01:ind02,:,:) ) ;*********************** ; plot ;*********************** plot = new (4 , graphic) ; create graphical array ; wks = gsn_open_wks("pdf","Precipitation_diff_100percent") wks = gsn_open_wks("pdf","Temperature_diff_100percent_mmmm") ; gsn_define_colormap(wks,"gui_default") ; gsn_define_colormap(wks,"prcp_1") gsn_define_colormap(wks,"BlWhRe") ; gsn_define_colormap(wks,"BlueYellowRed") res = True ; plot mods desired ; !!!!! any plot of data that is on a native grid, must use the "corners" ; method of zooming in on map. res@mpLimitMode = "Corners" ; choose range of map res@mpLeftCornerLatF = xlat(0,0) res@mpLeftCornerLonF = xlon(0,0) res@mpRightCornerLatF = xlat(nlon-1,nlat-1) res@mpRightCornerLonF = xlon(nlon-1,nlat-1) ; The following 4 pieces of information are REQUIRED to properly display ; data on a native grid. This data should be specified somewhere in the ; model itself. prj = a@projection clon = a@longitude_of_projection_origin if (prj .eq. "LAMCON") then trlats = a@standard_parallel res@mpProjection = "LambertConformal" res@mpLambertParallel1F = trlats(0) res@mpLambertParallel2F = trlats(1) res@mpLambertMeridianF = clon end if if (prj .eq. "NORMER") then res@mpProjection = "Mercator" end if if (prj .eq. "POLSTR") then clat = a@latitude_of_projection_origin res@mpProjection = "Stereographic" res@mpRelativeCenterLon = True res@mpCenterLonF = clon res@mpRelativeCenterLat = True res@mpCenterLatF = clat end if if (prj .eq. "ROTMER") then clat = a@latitude_of_projection_origin res@mpProjection = "Mercator" res@mpCenterLonF = clon res@mpCenterLatF = clat end if ; usually, when data is placed onto a map, it is TRANSFORMED to the specified ; projection. Since this model is already on a native lambert conformal grid, ; we want to turn OFF the tranformation. ; res@lbLabelBarOn = True ;False ; turn off individual lb's res@lbLabelBarOn = False ; No single label bar res@gsnDraw = False res@gsnFrame = False res@tfDoNDCOverlay = True ; do not transform res@mpOutlineBoundarySets = "National" res@pmTickMarkDisplayMode = "Always" res@mpDataSetName = "Earth..4" res@mpDataBaseVersion = "Ncarg4_1" res@mpGeophysicalLineThicknessF = 1.5 res@mpNationalLineThicknessF = 1.5 res@gsnPaperOrientation = "Landscape" ;"Portrait" res@gsnMaximize = True ; must include w/ Paper Orientation res@cnFillOn = True ; color plot desired res@cnLinesOn = False ; no contour lines res@gsnSpreadColors = True res@gsnDraw = False ; don't draw yet res@gsnFrame = False ; don't advance frame yet ; res@mpMinLatF = 0 ; zoom in on a subregion ; res@mpMaxLatF = 12 ; res@mpMinLonF = 28 ; res@mpMaxLonF = 34. res@mpGeophysicalLineColor = "red" ; color of continental outlines res@mpPerimOn = True ; draw box around map res@mpGridLineDashPattern = 2 ; lat/lon lines as dashed res@gsnAddCyclic = False ; regional data don't add res@pmTickMarkDisplayMode = "Always" ; ; res@lbAutoManage = True ; res@cnLevelSelectionMode = "ExplicitLevels" ; set manual contour levels ; res@cnLevels =(/50,100,150,200,250,300,350,400,450,500,600,700,800,1000,1200,1400,1600,1800,2000/) ; res@cnLevels =(/3,5,7,10,12,15,17,20,22,25,27,30,32,35,37,40,50/) res@cnLevels =(/-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,-0.5,0.5,1,2,3,4,5,6,7,8,9,10/) ; prec@_FillValue = -1.e+34 ;var@missing_value = -1.e+30 ; prec = where(landmask.eq.0, prec@_FillValue, prec) res@tiMainFuncCode = "~" ; ; res@gsnLeftString = " Precipitation Bias " ; res@gsnLeftString = " Precipitation DRA-CTL" ; res@gsnRightString = "mm/day " ; i = 0 ; do while (i .lt. 12) ;do i=0,11 ; res@tiMainString = a@title + "~C~" + ut_string(times(i), ""); ; plot = gsn_csm_contour_map(wks,prec1(i,:,:)-prec(i,:,:),res) res@gsnLeftString = "Temperature DRA-CTL" res@gsnRightString = "~S~o~N~C " res@gsnCenterString = " DJF " plot(0) = gsn_csm_contour_map(wks,((prec1(11,0,:,:)-prec(11,0,:,:))+(prec1(0,0,:,:)-prec(0,0,:,:))+(prec1(1,0,:,:)-prec(1,0,:,:)))/3, res) ; create plot res@gsnCenterString = " MAM " plot(1) = gsn_csm_contour_map(wks,((prec1(2,0,:,:)-prec(2,0,:,:))+(prec1(3,0,:,:)-prec(3,0,:,:))+(prec1(4,0,:,:)-prec(4,0,:,:)))/3, res) ; create plot res@gsnCenterString = "JJA " plot(2) = gsn_csm_contour_map(wks,((prec1(5,0,:,:)-prec(5,0,:,:))+(prec1(6,0,:,:)-prec(6,0,:,:))+(prec1(7,0,:,:)-prec(7,0,:,:)))/3, res) ; create plot res@gsnCenterString = " SON " plot(3) = gsn_csm_contour_map(wks,((prec1(8,0,:,:)-prec(8,0,:,:))+(prec1(9,0,:,:)-prec(9,0,:,:))+(prec1(10,0,:,:)-prec(10,0,:,:)))/3, res) ; create plot resP = True ; panel options resP@gsnMaximize = True ; maximize image resP@lbLabelAngleF = 270 resP@gsnPanelLabelBar = True ; Add common label bar gsn_panel(wks,plot,(/2,2/),resP) ; plot = gsn_csm_contour_map(wks,prec(i,:,:)*30,res) ; res@gsnLeftString = " Precipitation CTL" ; res@gsnRightString = "mm/day " ; plot = gsn_csm_contour_map(wks,prec1(7,:,:)-prec(7,:,:),res) ; end do ; draw(plot) ; frame(wks) end