From 20d8a743f27a2897a9d6ba8810fa0257a089c361 Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Tue, 29 Dec 2015 16:54:51 -0700 Subject: [PATCH] Added smoothing fortran code. Tested and fixed interpolation routines. Fixed bug with SLP constant. --- wrf_open/var/ncl_reference/wrf_vinterp.f | 404 ++++++ wrf_open/var/ncl_reference/wrf_vinterpW.c | 1147 +++++++++++++++++ wrf_open/var/src/python/wrf/var/decorators.py | 47 +- wrf_open/var/src/python/wrf/var/extension.py | 2 +- wrf_open/var/src/python/wrf/var/interp.py | 67 +- wrf_open/var/src/python/wrf/var/wrfext.f90 | 59 +- wrf_open/var/test/ncl_get_var.ncl | 53 +- wrf_open/var/test/utests.py | 156 ++- 8 files changed, 1901 insertions(+), 34 deletions(-) create mode 100644 wrf_open/var/ncl_reference/wrf_vinterp.f create mode 100644 wrf_open/var/ncl_reference/wrf_vinterpW.c diff --git a/wrf_open/var/ncl_reference/wrf_vinterp.f b/wrf_open/var/ncl_reference/wrf_vinterp.f new file mode 100644 index 0000000..8be4844 --- /dev/null +++ b/wrf_open/var/ncl_reference/wrf_vinterp.f @@ -0,0 +1,404 @@ +CThe subroutines in this file were taken directly from RIP code written +C by Dr. Mark Stoelinga. They were modified by Sherrie +C Fredrick(NCAR/MMM) to work with NCL February 2015. +C NCLFORTSTART + subroutine wrf_monotonic(out,in,lvprs,cor,idir,delta, + & ew,ns,nz,icorsw) + implicit none + integer idir,ew,ns,nz,icorsw + double precision delta + double precision in(ew,ns,nz),out(ew,ns,nz) + double precision lvprs(ew,ns,nz),cor(ew,ns) +C NCLEND + + integer i,j,k,ripk,k300 + + do j=1,ns + do i=1,ew + if (icorsw.eq.1.and.cor(i,j).lt.0.) then + do k=1,nz + in(i,j,k)=-in(i,j,k) + enddo + endif + + +c +c First find k index that is at or below (height-wise) the 300 hPa +c level. +c + do k = 1,nz + ripk = nz-k+1 + if (lvprs(i,j,k) .le. 300.d0) then + k300=k + goto 40 + endif + enddo +c + 40 continue + + do k = k300, 1,-1 + if (idir.eq.1) then + out(i,j,k)=min(in(i,j,k),in(i,j,k+1)+delta) + elseif (idir.eq.-1) then + out(i,j,k)=max(in(i,j,k),in(i,j,k+1)-delta) + endif + enddo + + + do k = k300+1, nz + if (idir.eq.1) then + out(i,j,k)=max(in(i,j,k),in(i,j,k-1)-delta) + elseif (idir.eq.-1) then + out(i,j,k)=min(in(i,j,k),in(i,j,k-1)+delta) + endif + enddo + + end do + end do + + return + end + +c-------------------------------------------------------------------- + +C NCLFORTSTART + FUNCTION wrf_intrp_value (wvalp0,wvalp1,vlev,vcp0,vcp1,icase) + implicit none + + integer icase + double precision wvalp0,wvalp1,vlev,vcp0,vcp1 +C NCLEND + double precision valp0,valp1,rvalue,rgas,ussalr,sclht + + double precision wrf_intrp_value,chkdiff + + rgas = 287.04d0 !J/K/kg + ussalr = 0.0065d0 ! deg C per m + sclht = rgas*256.d0/9.81d0 + + valp0 = wvalp0 + valp1 = wvalp1 + if ( icase .eq. 2) then !GHT + valp0=exp(-wvalp0/sclht) + valp1=exp(-wvalp1/sclht) + end if + + chkdiff = vcp1 - vcp0 + if(chkdiff .eq. 0) then + print *,"bad difference in vcp's" + stop + end if + + rvalue = (vlev-vcp0)*(valp1-valp0)/(vcp1-vcp0)+valp0 + if (icase .eq. 2) then !GHT + wrf_intrp_value = -sclht*log(rvalue) + else + wrf_intrp_value = rvalue + endif + + return + end +c------------------------------------------------------------ +C NOTES: +c vcarray is the array holding the values for the vertical +c coordinate. +c It will always come in with the dimensions of +c the staggered U and V grid. +C NCLFORTSTART + + subroutine wrf_vintrp(datain,dataout,pres,tk,qvp,ght,terrain, + & sfp,smsfp,vcarray,interp_levels,numlevels, + & icase,ew,ns,nz,extrap,vcor,logp,rmsg) + + + implicit none + integer ew,ns,nz,icase,extrap + integer vcor,numlevels,logp + double precision datain(ew,ns,nz),pres(ew,ns,nz),tk(ew,ns,nz) + double precision ght(ew,ns,nz) + double precision terrain(ew,ns),sfp(ew,ns),smsfp(ew,ns) + double precision dataout(ew,ns,numlevels),qvp(ew,ns,nz) + double precision vcarray(ew,ns,nz) + double precision interp_levels(numlevels),rmsg +C NCLEND + integer njx,niy,nreqlvs,ripk + integer i,j,k,itriv,kupper + integer ifound,miy,mjx,isign + double precision rlevel,vlev,diff + double precision tempout(ew,ns),tmpvlev + double precision vcp1,vcp0,valp0,valp1 + double precision rgas,rgasmd,sclht,ussalr,cvc,eps + double precision qvlhsl,ttlhsl,vclhsl,vctophsl + double precision wrf_intrp_value + double precision plhsl,zlhsl,ezlhsl,tlhsl,psurf,pratio,tlev + double precision ezsurf,psurfsm,zsurf,qvapor,vt + double precision rconst,expon,exponi + double precision ezlev,plev,zlev,ptarget,dpmin,dp + double precision pbot,zbot,tbotextrap,e + double precision tlclc1,tlclc2,tlclc3,tlclc4 + double precision thtecon1,thtecon2,thtecon3 + double precision tlcl,gamma,cp,cpmd,gammamd,gammam + character cvcord*1 + + rgas = 287.04d0 !J/K/kg + rgasmd = .608d0 + ussalr = .0065d0 ! deg C per m + sclht = rgas*256.d0/9.81d0 + eps = 0.622d0 + rconst = -9.81d0/(rgas * ussalr) + expon = rgas*ussalr/9.81d0 + exponi = 1./expon + tlclc1 = 2840.d0 + tlclc2 = 3.5d0 + tlclc3 = 4.805d0 + tlclc4 = 55.d0 + thtecon1 = 3376.d0 ! K + thtecon2 = 2.54d0 + thtecon3 = 0.81d0 + cp = 1004.d0 + cpmd = 0.887d0 + gamma = rgas/cp + gammamd = rgasmd-cpmd + + if(vcor .eq. 1) then + cvcord = 'p' + else if((vcor .eq. 2) .or. (vcor .eq. 3)) then + cvcord = 'z' + else if((vcor .eq. 4) .or. (vcor .eq. 5)) then + cvcord = 't' + end if + + + miy = ns + mjx = ew + njx = ew + niy = ns + + + do j = 1,mjx + do i = 1,miy + tempout(j,i) = rmsg + end do + end do + + + + do nreqlvs = 1,numlevels + if(cvcord .eq. 'z') then +!Convert rlevel to meters from km + + rlevel = interp_levels(nreqlvs) * 1000.d0 + vlev = exp(-rlevel/sclht) + else if(cvcord .eq. 'p') then + vlev = interp_levels(nreqlvs) + else if(cvcord .eq. 't') then + vlev = interp_levels(nreqlvs) + end if + + + do j=1,mjx + do i=1,miy +cGet the interpolated value that is within the model domain + ifound = 0 + do k = 1,nz-1 + ripk = nz-k+1 + vcp1 = vcarray(j,i,ripk-1) + vcp0 = vcarray(j,i,ripk) + valp0 = datain(j,i,ripk) + valp1 = datain(j,i,ripk-1) + if ((vlev.ge.vcp0.and.vlev.le.vcp1) .or. + & (vlev.le.vcp0.and.vlev.ge.vcp1)) then +c print *,i,j,valp0,valp1 + if((valp0 .eq. rmsg).or.(valp1 .eq. rmsg)) then + tempout(j,i) = rmsg + ifound=1 + else + if(logp .eq. 1) then + vcp1 = log(vcp1) + vcp0 = log(vcp0) + if(vlev .eq. 0.0d0) then + print *,"Pressure value = 0" + print *,"Unable to take log of 0" + stop + end if + tmpvlev = log(vlev) + else + tmpvlev = vlev + end if + tempout(j,i) = wrf_intrp_value(valp0,valp1, + & tmpvlev,vcp0,vcp1,icase) +c print *,"one ",i,j,tempout(j,i) + ifound=1 + end if + goto 115 + end if + end do !end for the k loop + 115 continue + + + if (ifound.eq.1) then !Grid point is in the model domain + goto 333 + end if + +cIf the user has requested no extrapolatin then just assign +call values above or below the model level to rmsg. + if(extrap .eq. 0) then + tempout(j,i) = rmsg + goto 333 + end if + + +c The grid point is either above or below the model domain +c +c First we will check to see if the grid point is above the +c model domain. + vclhsl = vcarray(j,i,1) !lowest model level + vctophsl = vcarray(j,i,nz)!highest model level + diff = vctophsl-vclhsl + isign = nint(diff/abs(diff)) +C + if(isign*vlev.ge.isign*vctophsl) then +C Assign the highest model level to the out array + tempout(j,i)=datain(j,i,nz) +C print *,"at warn",j,i,tempout(j,i) + goto 333 + endif + + +c +c Only remaining possibility is that the specified level is below +c lowest model level. If lowest model level value is missing, +c set interpolated value to missing. +c + if (datain(i,j,1) .eq. rmsg) then + tempout(j,i) = rmsg + goto 333 + endif + +c +c If the field comming in is not a pressure,temperature or height +C field we can set the output array to the value at the lowest +c model level. +c + tempout(j,i) = datain(j,i,1) +c +c For the special cases of pressure on height levels or height on +c pressure levels, or temperature-related variables on pressure or +c height levels, perform a special extrapolation based on +c US Standard Atmosphere. Here we calcualate the surface pressure +c with the altimeter equation. This is how RIP calculates the +c surface pressure. +c + if (icase.gt.0) then + plhsl = pres(j,i,1) * 0.01d0 !pressure at lowest model level + zlhsl = ght(j,i,1) !grid point height a lowest model level + ezlhsl = exp(-zlhsl/sclht) + tlhsl = tk(j,i,1) !temperature in K at lowest model level + zsurf = terrain(j,i) + qvapor = max((qvp(j,i,1)*.001d0),1.e-15) +c virtual temperature +c vt = tlhsl * (eps + qvapor)/(eps*(1.0 + qvapor)) +c psurf = plhsl * (vt/(vt+ussalr * (zlhsl-zsurf)))**rconst + psurf = sfp(j,i) + psurfsm = smsfp(j,i) + ezsurf = exp(-zsurf/sclht) + +cThe if for checking above ground + if ((cvcord.eq.'z'.and.vlev.lt.ezsurf).or. + & (cvcord.eq.'p'.and.vlev.lt.psurf)) then +c +c We are below the lowest data level but above the ground. +c Use linear interpolation (linear in prs and exp-height). +c + if (cvcord.eq.'p') then + plev=vlev + ezlev=((plev-plhsl)*ezsurf+(psurf-plev)*ezlhsl)/ + & (psurf-plhsl) + zlev=-sclht*log(ezlev) + if (icase .eq. 2) then + tempout(j,i)=zlev + goto 333 + endif + + elseif (cvcord.eq.'z') then + ezlev=vlev + zlev=-sclht*log(ezlev) + plev=((ezlev-ezlhsl)*psurf+(ezsurf-ezlev)*plhsl)/ + & (ezsurf-ezlhsl) + if (icase .eq. 1) then + tempout(j,i)=plev + goto 333 + endif + endif + + else !else for checking above ground + ptarget=psurfsm-150.d0 + dpmin=1.e4 + do k=1,nz + ripk = nz-k+1 + dp=abs((pres(j,i,ripk) * 0.01d0)-ptarget) + if (dp.gt.dpmin) goto 334 + dpmin=min(dpmin,dp) + enddo + 334 kupper=k-1 + + ripk = nz - kupper + 1 + pbot = max(plhsl,psurf) + zbot = min(zlhsl,zsurf) + pratio = pbot/(pres(j,i,ripk) * 0.01d0) + tbotextrap = tk(j,i,ripk)*(pratio)**expon +c virtual temperature + vt = tbotextrap * (eps + qvapor)/(eps*(1.0d0+qvapor)) + if (cvcord.eq.'p') then + plev=vlev + zlev=zbot+vt/ussalr*(1.-(vlev/pbot)**expon) + if(icase .eq. 2) then + tempout(j,i)=zlev + goto 333 + endif + elseif (cvcord.eq.'z') then + zlev=-sclht*log(vlev) + plev=pbot*(1.+ussalr/vt*(zbot-zlev))**exponi + if (icase .eq. 1) then + tempout(j,i)=plev + goto 333 + endif + endif + end if !end if for checking above ground + end if !for icase gt 0 + + + if(icase .gt. 2) then !extrapolation for temperature + tlev=tlhsl+(zlhsl-zlev)*ussalr + qvapor = max(qvp(j,i,1),1.e-15) + gammam = gamma*(1.+gammamd*qvapor) + if(icase .eq. 3) then + tempout(j,i) = tlev - 273.16d0 + else if(icase .eq. 4) then + tempout(j,i) = tlev +C Potential temperature - theta + else if (icase. eq. 5) then + tempout(j,i)=tlev*(1000.d0/plev)**gammam +C extraolation for equivalent potential temperature + else if (icase .eq. 6) then + e = qvapor*plev/(eps+qvapor) + tlcl = tlclc1/(log(tlev**tlclc2/e)-tlclc3)+tlclc4 + tempout(j,i)=tlev*(1000.d0/plev)**(gammam)* + & exp((thtecon1/tlcl-thtecon2)*qvapor* + & (1.+thtecon3*qvapor)) + end if + end if + + 333 continue + + end do + end do + ! print *,"----done----",interp_levels(nreqlvs) + do i = 1,njx + do j = 1,niy + dataout(i,j,nreqlvs) = tempout(i,j) + end do + end do + end do !end for the nreqlvs + return + end !wrf_vinterp diff --git a/wrf_open/var/ncl_reference/wrf_vinterpW.c b/wrf_open/var/ncl_reference/wrf_vinterpW.c new file mode 100644 index 0000000..5f64b9b --- /dev/null +++ b/wrf_open/var/ncl_reference/wrf_vinterpW.c @@ -0,0 +1,1147 @@ +#include +#include "wrapper.h" + +extern void NGCALLF(wrf_vintrp,WRF_VINTRP)(double *, double *, double *, + double *, double *, double *, + double *, double *, double *, + double *, double *, int *, + int *, int *, int *, int *, + int *, int *, int *, double *); + +extern void NGCALLF(wrf_monotonic,WRF_MONOTONIC)(double *, double *, double *, + double *, int *, double *, + int *, int *, int *, int *); + + +NhlErrorTypes wrf_vintrp_W( void ) +{ + +/* + * Input variables + */ +/* + * Argument # 0 + */ + void *field; + double *tmp_field; + int ndims_field; + ng_size_t dsizes_field[NCL_MAX_DIMENSIONS]; + int has_missing_field; + NclScalar missing_field, missing_flt_field, missing_dbl_field; + NclBasicDataTypes type_field; + +/* + * Argument # 1 + */ + void *pres; + double *tmp_pres; + int ndims_pres; + ng_size_t dsizes_pres[NCL_MAX_DIMENSIONS]; + NclBasicDataTypes type_pres; + +/* + * Argument # 2 + */ + void *tk; + double *tmp_tk; + int ndims_tk; + ng_size_t dsizes_tk[NCL_MAX_DIMENSIONS]; + NclBasicDataTypes type_tk; + +/* + * Argument # 3 + */ + void *qvp; + double *tmp_qvp; + int ndims_qvp; + ng_size_t dsizes_qvp[NCL_MAX_DIMENSIONS]; + NclBasicDataTypes type_qvp; + +/* + * Argument # 4 + */ + void *ght; + double *tmp_ght; + int ndims_ght; + ng_size_t dsizes_ght[NCL_MAX_DIMENSIONS]; + NclBasicDataTypes type_ght; + +/* + * Argument # 5 + */ + void *ter; + double *tmp_ter; + ng_size_t dsizes_ter[2]; + NclBasicDataTypes type_ter; + +/* + * Argument # 6 + */ + void *sfp; + double *tmp_sfp; + int ndims_sfp; + ng_size_t dsizes_sfp[NCL_MAX_DIMENSIONS]; + NclBasicDataTypes type_sfp; + +/* + * Argument # 7 + */ + void *smsfp; + double *tmp_smsfp; + int ndims_smsfp; + ng_size_t dsizes_smsfp[NCL_MAX_DIMENSIONS]; + NclBasicDataTypes type_smsfp; + +/* + * Argument # 8 + */ + void *vcarray; + double *tmp_vcarray; + int ndims_vcarray; + ng_size_t dsizes_vcarray[NCL_MAX_DIMENSIONS]; + NclBasicDataTypes type_vcarray; + +/* + * Argument # 9 + */ + void *intrp_levels; + double *tmp_intrp_levels; + ng_size_t dsizes_intrp_levels[1]; + NclBasicDataTypes type_intrp_levels; + +/* + * Arguments # 10-13 + */ + int *icase, *extrap, *vcor, *logp; + +/* + * Return variable + */ + void *field_out; + double *tmp_field_out; + ng_size_t *dsizes_field_out; + NclScalar missing_field_out; + NclBasicDataTypes type_field_out; + +/* + * Various + */ + ng_size_t ninlev, nlat, nlon, ninlevlatlon, nlatlon; + ng_size_t noutlev, noutlevlatlon; + ng_size_t index_field, index_sfp, index_field_out; + ng_size_t i, size_leftmost, size_output; + int ininlev, inlat, inlon, inoutlev, ret; + +/* + * Retrieve parameters. + * + * Note any of the pointer parameters can be set to NULL, which + * implies you don't care about its value. + */ +/* + * Get argument # 0 + */ + field = (void*)NclGetArgValue( + 0, + 14, + &ndims_field, + dsizes_field, + &missing_field, + &has_missing_field, + &type_field, + DONT_CARE); + +/* + * Check dimension sizes. + */ + if(ndims_field < 3 || ndims_field > 4) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: The field array must be 3D or 4D"); + return(NhlFATAL); + } + +/* + * Coerce missing value to double if necessary. + */ + coerce_missing(type_field,has_missing_field,&missing_field, + &missing_dbl_field,&missing_flt_field); + + ninlev = dsizes_field[ndims_field-3]; + nlat = dsizes_field[ndims_field-2]; + nlon = dsizes_field[ndims_field-1]; + +/* + * Test dimension sizes. + */ + if(ninlev > INT_MAX || nlat > INT_MAX || nlon > INT_MAX) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: one of bottom_top, south_north, or west_east is greater than INT_MAX"); + return(NhlFATAL); + } + ininlev = (int) ninlev; + inlat = (int) nlat; + inlon = (int) nlon; + +/* + * Get argument # 1 + */ + pres = (void*)NclGetArgValue( + 1, + 14, + &ndims_pres, + dsizes_pres, + NULL, + NULL, + &type_pres, + DONT_CARE); + +/* + * Check dimension sizes. + */ + if(ndims_pres != ndims_field) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: The pres and field arrays must have the same dimensionality"); + return(NhlFATAL); + } + else { + for(i = 0; i < ndims_field; i++) { + if(dsizes_pres[i] != dsizes_field[i]) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: The pres and field arrays must have the same dimensionality"); + return(NhlFATAL); + } + } + } + +/* + * Get argument # 2 + */ + tk = (void*)NclGetArgValue( + 2, + 14, + &ndims_tk, + dsizes_tk, + NULL, + NULL, + &type_tk, + DONT_CARE); + +/* + * Check dimension sizes. + */ + if(ndims_tk != ndims_field) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: The tk and field arrays must have the same dimensionality"); + return(NhlFATAL); + } + else { + for(i = 0; i < ndims_field; i++) { + if(dsizes_tk[i] != dsizes_field[i]) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: The tk and field arrays must have the same dimensionality"); + return(NhlFATAL); + } + } + } + + +/* + * Get argument # 3 + */ + qvp = (void*)NclGetArgValue( + 3, + 14, + &ndims_qvp, + dsizes_qvp, + NULL, + NULL, + &type_qvp, + DONT_CARE); + +/* + * Check dimension sizes. + */ + if(ndims_qvp != ndims_field) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: The qvp and field arrays must have the same dimensionality"); + return(NhlFATAL); + } + else { + for(i = 0; i < ndims_field; i++) { + if(dsizes_qvp[i] != dsizes_field[i]) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: The qvp and field arrays must have the same dimensionality"); + return(NhlFATAL); + } + } + } + +/* + * Get argument # 4 + */ + ght = (void*)NclGetArgValue( + 4, + 14, + &ndims_ght, + dsizes_ght, + NULL, + NULL, + &type_ght, + DONT_CARE); + +/* + * Check dimension sizes. + */ + if(ndims_ght != ndims_field) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: The ght and field arrays must have the same dimensionality"); + return(NhlFATAL); + } + else { + for(i = 0; i < ndims_field; i++) { + if(dsizes_ght[i] != dsizes_field[i]) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: The ght and field arrays must have the same dimensionality"); + return(NhlFATAL); + } + } + } + +/* + * Get argument # 5 + */ + ter = (void*)NclGetArgValue( + 5, + 14, + NULL, + dsizes_ter, + NULL, + NULL, + &type_ter, + DONT_CARE); +/* + * Check dimension sizes for ter. It can either be 2D, or one fewer + * dimensions than field. + */ + if(dsizes_ter[0] != nlat || dsizes_ter[1] != nlon) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: The dimensions of ter must be south_north x west_east"); + return(NhlFATAL); + } + +/* + * Get argument # 6 + */ + sfp = (void*)NclGetArgValue( + 6, + 14, + &ndims_sfp, + dsizes_sfp, + NULL, + NULL, + &type_sfp, + DONT_CARE); + +/* + * Check dimension sizes. + */ + if(ndims_sfp != (ndims_field-1)) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: The sfp array must have the same dimensionality as the field array, minus the level dimension"); + return(NhlFATAL); + } + else { + for(i = 0; i < ndims_field-3; i++) { + if(dsizes_sfp[i] != dsizes_field[i]) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: The sfp array must have the same dimensionality as the field array, minus the level dimension"); + return(NhlFATAL); + } + } + if(dsizes_sfp[ndims_sfp-2] != nlat || dsizes_sfp[ndims_sfp-1] != nlon) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: The sfp array must have the same dimensionality as the field array, minus the level dimension"); + return(NhlFATAL); + } + } + + +/* + * Get argument # 7 + */ + smsfp = (void*)NclGetArgValue( + 7, + 14, + &ndims_smsfp, + dsizes_smsfp, + NULL, + NULL, + &type_smsfp, + DONT_CARE); + +/* + * Check dimension sizes. + */ + if(ndims_smsfp != (ndims_field-1)) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: The smsfp array must have the same dimensionality as the field array, minus the level dimension"); + return(NhlFATAL); + } + else { + for(i = 0; i < ndims_field-3; i++) { + if(dsizes_smsfp[i] != dsizes_field[i]) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: The smsfp array must have the same dimensionality as the field array, minus the level dimension"); + return(NhlFATAL); + } + } + if(dsizes_smsfp[ndims_smsfp-2] != nlat || + dsizes_smsfp[ndims_smsfp-1] != nlon) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: The smsfp array must have the same dimensionality as the field array, minus the level dimension"); + return(NhlFATAL); + } + } + + +/* + * Get argument # 8 + */ + vcarray = (void*)NclGetArgValue( + 8, + 14, + &ndims_vcarray, + dsizes_vcarray, + NULL, + NULL, + &type_vcarray, + DONT_CARE); + +/* + * Check dimension sizes. + */ + if(ndims_vcarray != ndims_field) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: The vcarray and field arrays must have the same dimensionality"); + return(NhlFATAL); + } + else { + for(i = 0; i < ndims_field; i++) { + if(dsizes_vcarray[i] != dsizes_field[i]) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: The vcarray and field arrays must have the same dimensionality"); + return(NhlFATAL); + } + } + } + +/* + * Get argument # 9 + */ + intrp_levels = (void*)NclGetArgValue( + 9, + 14, + NULL, + dsizes_intrp_levels, + NULL, + NULL, + &type_intrp_levels, + DONT_CARE); + + noutlev = dsizes_intrp_levels[0]; + +/* + * Test dimension sizes. + */ + if(noutlev > INT_MAX) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: the # of output levels is greater than INT_MAX"); + return(NhlFATAL); + } + inoutlev = (int) noutlev; + + +/* + * Get argument # 10 + */ + icase = (int*)NclGetArgValue( + 10, + 14, + NULL, + NULL, + NULL, + NULL, + NULL, + DONT_CARE); +/* + * Get argument # 11 + */ + extrap = (int*)NclGetArgValue( + 11, + 14, + NULL, + NULL, + NULL, + NULL, + NULL, + DONT_CARE); +/* + * Get argument # 12 + */ + vcor = (int*)NclGetArgValue( + 12, + 14, + NULL, + NULL, + NULL, + NULL, + NULL, + DONT_CARE); +/* + * Get argument # 13 + */ + logp = (int*)NclGetArgValue( + 13, + 14, + NULL, + NULL, + NULL, + NULL, + NULL, + DONT_CARE); + +/* + * Calculate size of leftmost dimensions. + */ + size_leftmost = 1; + for(i = 0; i < ndims_field-3; i++) size_leftmost *= dsizes_field[i]; + +/* + * Allocate space for coercing input arrays. If any of the input + * is already double, then we don't need to allocate space for + * temporary arrays, because we'll just change the pointer into + * the void array appropriately. + */ +/* + * Allocate space for tmp_field. + */ + nlatlon = nlat * nlon; + ninlevlatlon = ninlev * nlatlon; + noutlevlatlon = noutlev * nlatlon; + + if(type_field != NCL_double) { + tmp_field = (double *)calloc(ninlevlatlon,sizeof(double)); + if(tmp_field == NULL) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: Unable to allocate memory for coercing field array to double"); + return(NhlFATAL); + } + } + +/* + * Allocate space for tmp_pres. + */ + if(type_pres != NCL_double) { + tmp_pres = (double *)calloc(ninlevlatlon,sizeof(double)); + if(tmp_pres == NULL) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: Unable to allocate memory for coercing pressure array to double"); + return(NhlFATAL); + } + } + +/* + * Allocate space for tmp_tk. + */ + if(type_tk != NCL_double) { + tmp_tk = (double *)calloc(ninlevlatlon,sizeof(double)); + if(tmp_tk == NULL) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: Unable to allocate memory for coercing tk array to double"); + return(NhlFATAL); + } + } + +/* + * Allocate space for tmp_qvp. + */ + if(type_qvp != NCL_double) { + tmp_qvp = (double *)calloc(ninlevlatlon,sizeof(double)); + if(tmp_qvp == NULL) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: Unable to allocate memory for coercing qvp array to double"); + return(NhlFATAL); + } + } + +/* + * Allocate space for tmp_ght. + */ + if(type_ght != NCL_double) { + tmp_ght = (double *)calloc(ninlevlatlon,sizeof(double)); + if(tmp_ght == NULL) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: Unable to allocate memory for coercing ght array to double"); + return(NhlFATAL); + } + } + +/* + * Coerce ter to double, if necessary. + */ + tmp_ter = coerce_input_double(ter,type_ter,nlatlon,0,NULL,NULL); + if(tmp_ter == NULL) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: Unable to coerce ter to double precision"); + return(NhlFATAL); + } + +/* + * Allocate space for tmp_sfp. + */ + if(type_sfp != NCL_double) { + tmp_sfp = (double *)calloc(nlatlon,sizeof(double)); + if(tmp_sfp == NULL) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: Unable to allocate memory for coercing sfp array to double"); + return(NhlFATAL); + } + } + +/* + * Allocate space for tmp_smsfp. + */ + if(type_smsfp != NCL_double) { + tmp_smsfp = (double *)calloc(nlatlon,sizeof(double)); + if(tmp_smsfp == NULL) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: Unable to allocate memory for coercing smsfp array to double"); + return(NhlFATAL); + } + } + +/* + * Allocate space for tmp_vcarray. + */ + if(type_vcarray != NCL_double) { + tmp_vcarray = (double *)calloc(ninlevlatlon,sizeof(double)); + if(tmp_vcarray == NULL) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: Unable to allocate memory for coercing vcarray to double"); + return(NhlFATAL); + } + } +/* + * Coerce intrp_levels to double, if necessary. + */ + tmp_intrp_levels = coerce_input_double(intrp_levels,type_intrp_levels, + noutlev,0,NULL,NULL); + if(tmp_intrp_levels == NULL) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: Unable to allocate memory for coercing interp_levels array to double"); + return(NhlFATAL); + } + +/* + * The output type defaults to float, unless one or more input + * arrays are double. + */ + if(type_field == NCL_double || type_pres == NCL_double || + type_tk == NCL_double || type_qvp == NCL_double || + type_ght == NCL_double || type_ter == NCL_double || + type_sfp == NCL_double || type_smsfp == NCL_double || + type_vcarray == NCL_double) { + type_field_out = NCL_double; + } + else { + type_field_out = NCL_float; + } + +/* + * Allocate space for output array and set a missing value. + * Note: a missing value is returned even if the input doesn't + * have a missing value, because the output may contain missing + * values after interpolation is done. + */ + size_output = size_leftmost * noutlevlatlon; + if(type_field_out != NCL_double) { + field_out = (void *)calloc(size_output, sizeof(float)); + tmp_field_out = (double *)calloc(noutlevlatlon,sizeof(double)); + if(field_out == NULL || tmp_field_out == NULL) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: Unable to allocate memory for temporary output array"); + return(NhlFATAL); + } + missing_field_out = missing_flt_field; + } + else { + field_out = (void *)calloc(size_output, sizeof(double)); + if(field_out == NULL) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: Unable to allocate memory for output array"); + return(NhlFATAL); + } + missing_field_out = missing_dbl_field; + } + +/* + * Allocate space for output dimension sizes and set them. + */ + dsizes_field_out = (ng_size_t*)calloc(ndims_field,sizeof(ng_size_t)); + if( dsizes_field_out == NULL ) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: Unable to allocate memory for holding dimension sizes"); + return(NhlFATAL); + } + for(i = 0; i < ndims_field-3; i++) dsizes_field_out[i] = dsizes_field[i]; + dsizes_field_out[ndims_field-3] = noutlev; + dsizes_field_out[ndims_field-2] = nlat; + dsizes_field_out[ndims_field-1] = nlon; + +/* + * Loop across leftmost dimensions and call the Fortran routine for each + * subsection of the input arrays. + */ + index_field = index_sfp = index_field_out = 0; + + for(i = 0; i < size_leftmost; i++) { +/* + * Coerce subsection of field (tmp_field) to double if necessary. + */ + if(type_field != NCL_double) { + coerce_subset_input_double(field,tmp_field,index_field,type_field, + ninlevlatlon,0,NULL,NULL); + } + else { + tmp_field = &((double*)field)[index_field]; + } + +/* + * Coerce subsection of pres (tmp_pres) to double if necessary. + */ + if(type_pres != NCL_double) { + coerce_subset_input_double(pres,tmp_pres,index_field,type_pres,ninlevlatlon,0,NULL,NULL); + } + else { + tmp_pres = &((double*)pres)[index_field]; + } + +/* + * Coerce subsection of tk (tmp_tk) to double if necessary. + */ + if(type_tk != NCL_double) { + coerce_subset_input_double(tk,tmp_tk,index_field,type_tk, + ninlevlatlon,0,NULL,NULL); + } + else { + tmp_tk = &((double*)tk)[index_field]; + } + +/* + * Coerce subsection of qvp (tmp_qvp) to double if necessary. + */ + if(type_qvp != NCL_double) { + coerce_subset_input_double(qvp,tmp_qvp,index_field,type_qvp, + ninlevlatlon,0,NULL,NULL); + } + else { + tmp_qvp = &((double*)qvp)[index_field]; + } + +/* + * Coerce subsection of ght (tmp_ght) to double if necessary. + */ + if(type_ght != NCL_double) { + coerce_subset_input_double(ght,tmp_ght,index_field,type_ght, + ninlevlatlon,0,NULL,NULL); + } + else { + tmp_ght = &((double*)ght)[index_field]; + } + +/* + * Coerce subsection of sfp (tmp_sfp) to double if necessary. + */ + if(type_sfp != NCL_double) { + coerce_subset_input_double(sfp,tmp_sfp,index_sfp,type_sfp, + nlatlon,0,NULL,NULL); + } + else { + tmp_sfp = &((double*)sfp)[index_sfp]; + } + +/* + * Coerce subsection of smsfp (tmp_smsfp) to double if necessary. + */ + if(type_smsfp != NCL_double) { + coerce_subset_input_double(smsfp,tmp_smsfp,index_sfp,type_smsfp, + nlatlon,0,NULL,NULL); + } + else { + tmp_smsfp = &((double*)smsfp)[index_sfp]; + } + +/* + * Coerce subsection of vcarray (tmp_vcarray) to double if necessary. + */ + if(type_vcarray != NCL_double) { + coerce_subset_input_double(vcarray,tmp_vcarray,index_field,type_vcarray, + ninlevlatlon,0,NULL,NULL); + } + else { + tmp_vcarray = &((double*)vcarray)[index_field]; + } + + +/* + * Point temporary output array to void output array if appropriate. + */ + if(type_field_out == NCL_double) { + tmp_field_out = &((double*)field_out)[index_field_out]; + } + +/* + * Call the Fortran routine. + */ + NGCALLF(wrf_vintrp,WRF_VINTRP)(tmp_field, tmp_field_out, tmp_pres, + tmp_tk, tmp_qvp, tmp_ght, tmp_ter, + tmp_sfp, tmp_smsfp, tmp_vcarray, + tmp_intrp_levels, &inoutlev, icase, + &inlon, &inlat, &ininlev, extrap, vcor, + logp, &missing_dbl_field.doubleval); + +/* + * Coerce output back to float if necessary. + */ + if(type_field_out == NCL_float) { + coerce_output_float_only(field_out,tmp_field_out,noutlevlatlon, + index_field_out); + } + index_field += ninlevlatlon; + index_field_out += noutlevlatlon; + index_sfp += nlatlon; + } + +/* + * Free unneeded memory. + */ + + if(type_field != NCL_double) NclFree(tmp_field); + if(type_pres != NCL_double) NclFree(tmp_pres); + if(type_tk != NCL_double) NclFree(tmp_tk); + if(type_qvp != NCL_double) NclFree(tmp_qvp); + if(type_ght != NCL_double) NclFree(tmp_ght); + if(type_ter != NCL_double) NclFree(tmp_ter); + if(type_sfp != NCL_double) NclFree(tmp_sfp); + if(type_smsfp != NCL_double) NclFree(tmp_smsfp); + if(type_vcarray != NCL_double) NclFree(tmp_vcarray); + if(type_intrp_levels != NCL_double) NclFree(tmp_intrp_levels); + if(type_field_out != NCL_double) NclFree(tmp_field_out); + +/* + * Return value back to NCL script. + */ + ret = NclReturnValue(field_out,ndims_field,dsizes_field_out, + &missing_field_out,type_field_out,0); + + NclFree(dsizes_field_out); + return(ret); +} + +NhlErrorTypes wrf_monotonic_W( void ) +{ + +/* + * Input variables + */ +/* + * Argument # 0 + */ + void *x; + double *tmp_x; + int ndims_x; + ng_size_t dsizes_x[NCL_MAX_DIMENSIONS]; + NclBasicDataTypes type_x; + +/* + * Argument # 1 + */ + void *pres; + double *tmp_pres; + int ndims_pres; + ng_size_t dsizes_pres[NCL_MAX_DIMENSIONS]; + NclBasicDataTypes type_pres; + +/* + * Argument # 2 + */ + void *cor; + double *tmp_cor; + ng_size_t dsizes_cor[2]; + NclBasicDataTypes type_cor; + +/* + * Argument # 3 + */ + int *idir; +/* + * Argument # 4 + */ + void *delta; + double *tmp_delta; + NclBasicDataTypes type_delta; + +/* + * Argument # 5 + */ + int *icorsw; +/* + * Return variable + */ + void *xout; + double *tmp_xout; + NclBasicDataTypes type_xout; + + +/* + * Various + */ + ng_size_t nlev, nlat, nlon, nlevlatlon, nlatlon; + ng_size_t i, index_x, size_leftmost, size_output; + int inlev, inlat, inlon, ret; +/* + * Retrieve parameters. + * + * Note any of the pointer parameters can be set to NULL, which + * implies you don't care about its value. + */ +/* + * Get argument # 0 + */ + x = (void*)NclGetArgValue( + 0, + 6, + &ndims_x, + dsizes_x, + NULL, + NULL, + &type_x, + DONT_CARE); + +/* + * Check dimension sizes. + */ + if(ndims_x < 3) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_monotonic: The x array must have at least three dimensions"); + return(NhlFATAL); + } + nlev = dsizes_x[ndims_x-3]; + nlat = dsizes_x[ndims_x-2]; + nlon = dsizes_x[ndims_x-1]; + +/* + * Test dimension sizes. + */ + if(nlev > INT_MAX || nlat > INT_MAX || nlon > INT_MAX) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_monotonic: one of bottom_top, south_north, or west_east is greater than INT_MAX"); + return(NhlFATAL); + } + inlev = (int) nlev; + inlat = (int) nlat; + inlon = (int) nlon; + + +/* + * Get argument # 1 + */ + pres = (void*)NclGetArgValue( + 1, + 6, + &ndims_pres, + dsizes_pres, + NULL, + NULL, + &type_pres, + DONT_CARE); + +/* + * Check dimension sizes. + */ + if(ndims_pres != ndims_x) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_monotonic: The pres and x arrays must have the same dimensionality"); + return(NhlFATAL); + } + else { + for(i = 0; i < ndims_x; i++) { + if(dsizes_pres[i] != dsizes_x[i]) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_monotonic: The pres and x arrays must have the same dimensionality"); + return(NhlFATAL); + } + } + } + +/* + * Get argument # 2 + */ + cor = (void*)NclGetArgValue( + 2, + 6, + NULL, + dsizes_cor, + NULL, + NULL, + &type_cor, + DONT_CARE); + + if(dsizes_cor[0] != nlat || dsizes_cor[1] != nlon) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_vintrp: The dimensions of cor must be south_north x west_east"); + return(NhlFATAL); + } + +/* + * Get argument # 3 + */ + idir = (int*)NclGetArgValue( + 3, + 6, + NULL, + NULL, + NULL, + NULL, + NULL, + DONT_CARE); +/* + * Get argument # 4 + */ + delta = (void*)NclGetArgValue( + 4, + 6, + NULL, + NULL, + NULL, + NULL, + &type_delta, + DONT_CARE); +/* + * Get argument # 5 + */ + icorsw = (int*)NclGetArgValue( + 5, + 6, + NULL, + NULL, + NULL, + NULL, + NULL, + DONT_CARE); + +/* + * Calculate size of leftmost dimensions. + */ + size_leftmost = 1; + for(i = 0; i < ndims_x-3; i++) size_leftmost *= dsizes_x[i]; + +/* + * Allocate space for coercing input arrays. If any of the input + * is already double, then we don't need to allocate space for + * temporary arrays, because we'll just change the pointer into + * the void array appropriately. + */ + +/* + * Allocate space for tmp_x. + */ + nlatlon = nlat * nlon; + nlevlatlon = nlev * nlatlon; + + if(type_x != NCL_double) { + type_xout = NCL_float; + tmp_x = (double *)calloc(nlevlatlon,sizeof(double)); + if(tmp_x == NULL) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_monotonic: Unable to allocate memory for coercing input array to double"); + return(NhlFATAL); + } + } + else { + type_xout = NCL_double; + } +/* + * Allocate space for tmp_pres. + */ + if(type_pres != NCL_double) { + tmp_pres = (double *)calloc(nlevlatlon,sizeof(double)); + if(tmp_pres == NULL) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_monotonic: Unable to allocate memory for coercing input array to double"); + return(NhlFATAL); + } + } +/* + * Coerce cor to double, if necessary. + */ + tmp_cor = coerce_input_double(cor,type_cor,nlatlon,0,NULL,NULL); + if(tmp_cor == NULL) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_monotonic: Unable to coerce cor to double precision"); + return(NhlFATAL); + } + +/* + * Allocate space for tmp_delta. + */ + tmp_delta = coerce_input_double(delta,type_delta,1,0,NULL,NULL); + if(tmp_delta == NULL) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_monotonic: Unable to allocate memory for coercing input array to double"); + return(NhlFATAL); + } + +/* + * Calculate size of output array, which is same as x. + */ + size_output = size_leftmost * nlevlatlon; + +/* + * Allocate space for output array. + */ + if(type_xout != NCL_double) { + xout = (void *)calloc(size_output, sizeof(float)); + tmp_xout = (double *)calloc(nlevlatlon,sizeof(double)); + if(xout == NULL || tmp_xout == NULL) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_monotonic: Unable to allocate memory for temporary output array"); + return(NhlFATAL); + } + } + else { + xout = (void *)calloc(size_output, sizeof(double)); + if(xout == NULL) { + NhlPError(NhlFATAL,NhlEUNKNOWN,"wrf_monotonic: Unable to allocate memory for output array"); + return(NhlFATAL); + } + } + +/* + * Loop across leftmost dimensions and call the Fortran routine for each + * subsection of the input arrays. + */ + index_x = 0; + for(i = 0; i < size_leftmost; i++) { +/* + * Coerce subsection of x (tmp_x) to double if necessary. + */ + if(type_x != NCL_double) { + coerce_subset_input_double(x,tmp_x,index_x,type_x,nlevlatlon,0,NULL,NULL); + } + else { + tmp_x = &((double*)x)[index_x]; + } + +/* + * Coerce subsection of pres (tmp_pres) to double if necessary. + */ + if(type_pres != NCL_double) { + coerce_subset_input_double(pres,tmp_pres,index_x,type_pres,nlevlatlon,0,NULL,NULL); + } + else { + tmp_pres = &((double*)pres)[index_x]; + } + +/* + * Point temporary output array to void output array if appropriate. + */ + if(type_xout == NCL_double) tmp_xout = &((double*)xout)[index_x]; + +/* + * Call the Fortran routine. + */ + NGCALLF(wrf_monotonic,WRF_MONOTONIC)(tmp_xout, tmp_x, tmp_pres, + tmp_cor, idir, tmp_delta, &inlon, + &inlat, &inlev, icorsw); +/* + * Coerce output back to float if necessary. + */ + if(type_xout == NCL_float) { + coerce_output_float_only(xout,tmp_xout,nlevlatlon,index_x); + } + index_x += nlevlatlon; + } + +/* + * Free unneeded memory. + */ + if(type_x != NCL_double) NclFree(tmp_x); + if(type_pres != NCL_double) NclFree(tmp_pres); + if(type_cor != NCL_double) NclFree(tmp_cor); + if(type_delta != NCL_double) NclFree(tmp_delta); + if(type_xout != NCL_double) NclFree(tmp_xout); + +/* + * Return value back to NCL script. + */ + ret = NclReturnValue(xout,ndims_x,dsizes_x,NULL,type_xout,0); + return(ret); +} diff --git a/wrf_open/var/src/python/wrf/var/decorators.py b/wrf_open/var/src/python/wrf/var/decorators.py index 124f538..824f610 100644 --- a/wrf_open/var/src/python/wrf/var/decorators.py +++ b/wrf_open/var/src/python/wrf/var/decorators.py @@ -3,6 +3,7 @@ from inspect import getargspec from itertools import product import numpy as n +import numpy.ma as ma from wrf.var.units import do_conversion, check_units from wrf.var.destag import destagger @@ -66,7 +67,7 @@ def _calc_out_dims(outvar, left_dims): def handle_left_iter(ref_var_expected_dims, ref_var_idx=-1, ref_var_name=None, - ignore_args=(), ignore_kargs={}): + ignore_args=(), ignore_kargs=()): """Decorator to handle iterating over leftmost dimensions when using multiple files and/or multiple times. @@ -118,8 +119,7 @@ def handle_left_iter(ref_var_expected_dims, ref_var_idx=-1, new_args = [arg[left_and_slice_idxs] if i not in ignore_args else arg for i,arg in enumerate(args)] - - + # Slice the kargs if applicable new_kargs = {key:(val[left_and_slice_idxs] if key not in ignore_kargs else val) @@ -132,20 +132,47 @@ def handle_left_iter(ref_var_expected_dims, ref_var_idx=-1, # Output array if not out_inited: outdims = _calc_out_dims(res, extra_dims) - output = n.zeros(outdims, ref_var.dtype) - out_inited = True + if not isinstance(res, ma.MaskedArray): + output = n.zeros(outdims, ref_var.dtype) + masked = False + else: + output = ma.MaskedArray( + n.zeros(outdims, ref_var.dtype), + mask=n.zeros(outdims, n.bool_), + fill_value=res.fill_value) + masked = True + + out_inited = True - output[left_and_slice_idxs] = res[:] + if not masked: + output[left_and_slice_idxs] = res[:] + else: + output.data[left_and_slice_idxs] = res.data[:] + output.mask[left_and_slice_idxs] = res.mask[:] else: # This should be a list or a tuple (cape) if not out_inited: outdims = _calc_out_dims(res[0], extra_dims) - output = [n.zeros(outdims, ref_var.dtype) - for i in xrange(len(res))] + if not isinstance(res[0], ma.MaskedArray): + output = [n.zeros(outdims, ref_var.dtype) + for i in xrange(len(res))] + masked = False + else: + output = [ma.MaskedArray( + n.zeros(outdims, ref_var.dtype), + mask=n.zeros(outdims, n.bool_), + fill_value=res[0].fill_value) + for i in xrange(len(res))] + masked = True + out_inited = True for i,outarr in enumerate(res): - (output[i])[left_and_slice_idxs] = outarr[:] + if not masked: + output[i][left_and_slice_idxs] = outarr[:] + else: + output[i].data[left_and_slice_idxs] = outarr.data[:] + output[i].mask[left_and_slice_idxs] = outarr.mask[:] return output @@ -235,7 +262,7 @@ def uvmet_left_iter(): return indexing_decorator -def handle_casting(ref_idx=0, arg_idxs=(), karg_names=None,dtype=n.float64): +def handle_casting(ref_idx=0, arg_idxs=(), karg_names=(),dtype=n.float64): """Decorator to handle casting to/from required dtype used in numerical routine. diff --git a/wrf_open/var/src/python/wrf/var/extension.py b/wrf_open/var/src/python/wrf/var/extension.py index d11cda1..acc7e5a 100755 --- a/wrf_open/var/src/python/wrf/var/extension.py +++ b/wrf_open/var/src/python/wrf/var/extension.py @@ -32,7 +32,7 @@ def interpz3d(data3d,zdata,desiredloc,missingval): missingval) return res.T -@handle_left_iter(2,0) +@handle_left_iter(3,0) @handle_casting(arg_idxs=(0,1)) def interp2dxy(data3d,xy): res = f_interp2dxy(data3d.T, diff --git a/wrf_open/var/src/python/wrf/var/interp.py b/wrf_open/var/src/python/wrf/var/interp.py index a963000..76632dd 100755 --- a/wrf_open/var/src/python/wrf/var/interp.py +++ b/wrf_open/var/src/python/wrf/var/interp.py @@ -23,6 +23,11 @@ def interplevel(data3d,zdata,desiredloc,missingval=-99999): masked_r1 = ma.masked_values (r1, missingval) return masked_r1 +def _to_positive_idxs(shape, coord): + if (coord[-2] >= 0 and coord[-1] >=0): + return coord + return [x if (x >= 0) else shape[i]+x for (i,x) in enumerate(coord) ] + def _get_xy(xdim, ydim, pivot_point=None, angle=None, start_point=None, end_point=None): """Returns the x,y points for the horizontal cross section line. @@ -31,8 +36,8 @@ def _get_xy(xdim, ydim, pivot_point=None, angle=None, ydim - maximum y-dimension pivot_point - a pivot point of (south_north, west_east) (must be used with angle) angle - the angle through the pivot point in degrees - start_point - a start_point tuple of (south_north1, west_east1) - end_point - an end point tuple of (south_north2, west_east2) + start_point - a start_point sequence of [south_north1, west_east1] + end_point - an end point sequence of [south_north2, west_east2] """ @@ -131,9 +136,9 @@ def _get_xy(xdim, ydim, pivot_point=None, angle=None, return xy -# TODO: Need a decorator that can handle right dimensions based on what -# is returned. In this case, the result of the xy calculation determines -# output dimensions on right. +@handle_left_iter(3, 0, ignore_args=(2,3,4,5,6), + ignore_kargs=("missingval", "pivot_point", "angle", + "start_point", "end_point")) @handle_casting(arg_idxs=(0,1)) def vertcross(data3d, z, missingval=-99999, pivot_point=None,angle=None, @@ -151,20 +156,35 @@ def vertcross(data3d, z, missingval=-99999, """ + if pivot_point is not None: + pos_pivot = _to_positive_idxs(z.shape[-2:], pivot_point) + else: + pos_pivot = pivot_point + + if start_point is not None: + pos_start = _to_positive_idxs(z.shape[-2:], start_point) + else: + pos_start = start_point + + if end_point is not None: + pos_end = _to_positive_idxs(z.shape[-2:], end_point) + else: + pos_end = start_point + xdim = z.shape[-1] ydim = z.shape[-2] - xy = _get_xy(xdim, ydim, pivot_point, angle, start_point, end_point) + xy = _get_xy(xdim, ydim, pos_pivot, angle, pos_start, pos_end) # Interp z var2dz = interp2dxy(z, xy) # interp to constant z grid - if(var2dz[0,0] > var2dz[1,0]): # monotonically decreasing coordinate - z_max = floor(n.amax(z)/10)*10 # bottom value - z_min = ceil(n.amin(z)/10)*10 # top value + if(var2dz[0,0] > var2dz[-1,0]): # monotonically decreasing coordinate + z_max = floor(n.amax(z) / 10) * 10 # bottom value + z_min = ceil(n.amin(z) / 10) * 10 # top value dz = 10 - nlevels = int( (z_max-z_min)/dz) + nlevels = int((z_max-z_min) / dz) z_var2d = n.zeros((nlevels), dtype=z.dtype) z_var2d[0] = z_max dz = -dz @@ -172,7 +192,7 @@ def vertcross(data3d, z, missingval=-99999, z_max = n.amax(z) z_min = 0. dz = 0.01 * z_max - nlevels = int( z_max/dz ) + nlevels = int(z_max / dz) z_var2d = n.zeros((nlevels), dtype=z.dtype) z_var2d[0] = z_min @@ -189,9 +209,9 @@ def vertcross(data3d, z, missingval=-99999, return ma.masked_values(var2d, missingval) -# TODO: Need a decorator that can handle right dimensions based on what -# is returned. In this case, the result of the xy calculation determines -# output dimensions on right. +@handle_left_iter(2, 0, ignore_args=(1,2,3,4), + ignore_kargs=("pivot_point", "angle", + "start_point", "end_point")) @handle_casting(arg_idxs=(0,)) def interpline(data2d, pivot_point=None, angle=None, start_point=None, @@ -207,14 +227,29 @@ def interpline(data2d, pivot_point=None, """ - tmp_shape = [0] + [x for x in data2d.shape] + if pivot_point is not None: + pos_pivot = _to_positive_idxs(data2d.shape[-2:], pivot_point) + else: + pos_pivot = pivot_point + + if start_point is not None: + pos_start = _to_positive_idxs(data2d.shape[-2:], start_point) + else: + pos_start = start_point + + if end_point is not None: + pos_end = _to_positive_idxs(data2d.shape[-2:], end_point) + else: + pos_end = start_point + + tmp_shape = [1] + [x for x in data2d.shape] var2dtmp = n.zeros(tmp_shape, data2d.dtype) xdim = data2d.shape[-1] ydim = data2d.shape[-2] - xy = _get_xy(xdim, ydim, pivot_point, angle, start_point, end_point) + xy = _get_xy(xdim, ydim, pos_pivot, angle, pos_start, pos_end) var2dtmp[0,:,:] = data2d[:,:] diff --git a/wrf_open/var/src/python/wrf/var/wrfext.f90 b/wrf_open/var/src/python/wrf/var/wrfext.f90 index 7d5580c..3b23a82 100755 --- a/wrf_open/var/src/python/wrf/var/wrfext.f90 +++ b/wrf_open/var/src/python/wrf/var/wrfext.f90 @@ -149,7 +149,7 @@ SUBROUTINE f_computeslp(z,t,p,q,t_sea_level,t_surf,level,throw_exception,& ! Specific constants for assumptions made in this routine: - REAL(KIND=8), PARAMETER :: TC=273.16D0, PCONST=10000 + REAL(KIND=8), PARAMETER :: TC=273.16D0+17.5D0, PCONST=10000 LOGICAL, PARAMETER :: ridiculous_mm5_test=.TRUE. ! PARAMETER (ridiculous_mm5_test = .FALSE.) @@ -1819,5 +1819,62 @@ SUBROUTINE f_computectt(prs,tk,qci,qcw,qvp,ght,ter,haveqci,ctt,ew,ns,nz) RETURN END SUBROUTINE f_computectt +SUBROUTINE f_filter2d(a, b, missing, it, nx, ny) + IMPLICIT NONE + + INTEGER, INTENT(IN) :: nx, ny, it + REAL(KIND=8), DIMENSION(nx, ny), INTENT(IN) :: a + REAL(KIND=8), INTENT(IN) :: missing + REAL(KIND=8), DIMENSION(nx, ny), INTENT(OUT) :: b + + REAL(KIND=8), PARAMETER :: COEF=0.25D0 + + INTEGER :: i, j, iter + + DO iter=1,it + DO j=1,ny + DO i = 1,nx + b(i,j) = a(i,j) + END DO + END DO + + DO j=2,ny-1 + DO i=1,nx + IF (b(i,j-1).EQ.missing .OR. b(i,j).EQ.missing .OR. b(i,j+1).EQ.missing) THEN + a(i,j) = a(i,j) + ELSE + a(i,j) = a(i,j) + COEF*(b(i,j-1)-2*b(i,j) + b(i,j+1)) + END IF + END DO + END DO + + DO j=1,ny + DO i=2,nx-1 + IF (b(i-1,j).EQ.missing .OR. b(i,j).EQ.missing .OR. b(i+1,j).EQ.missing) THEN + a(i,j) = a(i,j) + ELSE + a(i,j) = a(i,j) + COEF*(b(i-1,j)-2*b(i,j)+b(i+1,j)) + END IF + END DO + END DO + c do j=1,ny + c do i=1,nx + c b(i,j) = a(i,j) + c enddo + c enddo + c do j=2,ny-1 + c do i=1,nx + c a(i,j) = a(i,j) - .99*coef*(b(i,j-1)-2*b(i,j)+b(i,j+1)) + c enddo + c enddo + c do j=1,ny + c do i=2,nx-1 + c a(i,j) = a(i,j) - .99*coef*(b(i-1,j)-2*b(i,j)+b(i+1,j)) + c enddo + c enddo + END DO + RETURN +END + diff --git a/wrf_open/var/test/ncl_get_var.ncl b/wrf_open/var/test/ncl_get_var.ncl index e386a56..e5037aa 100644 --- a/wrf_open/var/test/ncl_get_var.ncl +++ b/wrf_open/var/test/ncl_get_var.ncl @@ -94,7 +94,7 @@ end do end do - + setfileoption(fout,"DefineMode",True) ; Set global attributes @@ -126,6 +126,57 @@ fout->$wrf_vars[i]$ = (/d/) end do + ; Do the interpolation routines manually + ; 3D vertical cross section + z = wrf_user_getvar(input_file, "z", 0) ; grid point height + p = wrf_user_getvar(input_file, "pressure", 0) ; total pressure + + dimsz = dimsizes(z) + pivot = (/ dimsz(2)/2, dimsz(1)/2 /) ; pivot point is center of domain + + ht_cross = wrf_user_intrp3d(z,p,"v",pivot,90.0,False) + p_cross = wrf_user_intrp3d(p,z,"v",pivot,90.0,False) + + ht_cross_dims = dimsizes(ht_cross) + p_cross_dims = dimsizes(p_cross) + + ; 3D horizontal interpolation + plev = 500. ; 500 MB + z_500 = wrf_user_intrp3d(z,p,"h",plev,0.,False) + z_500_dims = dimsizes(z_500) + + ; 2D interpolation along line + t2 = wrf_user_getvar(input_file, "T2", 0) + dimst2 = dimsizes(t2) + pivot = (/ dimst2(1)/2, dimst2(0)/2 /) + + t2_line = wrf_user_intrp2d(t2, pivot, 90.0, False) + t2_line_dims = dimsizes(t2_line) + + filedimdef(fout, (/"ht_cross_vert", "ht_cross_horiz", "p_cross_vert", "p_cross_horiz"/), \ + (/ht_cross_dims(0), ht_cross_dims(1), p_cross_dims(0), p_cross_dims(1)/), \ + (/False,False,False,False/)) + filedimdef(fout, (/"z500_vert", "z500_horiz"/), \ + (/z_500_dims(0), z_500_dims(1) /), \ + (/False, False/)) + filedimdef(fout, (/"t2_line_horiz"/), \ + (/t2_line_dims(0) /), \ + (/False/)) + + filevardef(fout, "ht_cross", typeof(ht_cross), (/"ht_cross_vert", "ht_cross_horiz"/)) + filevardef(fout, "p_cross", typeof(p_cross), (/"p_cross_vert", "p_cross_horiz"/)) + filevardef(fout, "z_500", typeof(z_500), (/"z500_vert", "z500_horiz"/)) + filevardef(fout, "t2_line", typeof(t2_line), (/"t2_line_horiz"/)) + + filevarattdef(fout,"ht_cross", ht_cross) + filevarattdef(fout,"p_cross", p_cross) + filevarattdef(fout, "z_500", z_500) + filevarattdef(fout, "t2_line", t2_line) + fout->ht_cross = (/ht_cross/) + fout->p_cross = (/p_cross/) + fout->z_500 = (/z_500/) + fout->t2_line = (/t2_line/) + delete(fout) diff --git a/wrf_open/var/test/utests.py b/wrf_open/var/test/utests.py index cc7cfbc..77205cb 100644 --- a/wrf_open/var/test/utests.py +++ b/wrf_open/var/test/utests.py @@ -5,7 +5,7 @@ import numpy.ma as ma import os, sys import subprocess -from wrf.var import getvar +from wrf.var import getvar, interplevel, interpline, vertcross NCL_EXE = "/Users/ladwig/nclbuild/6.3.0/bin/ncl" TEST_FILE = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2010-06-13_21:00:00" @@ -76,7 +76,7 @@ def make_test(varname, wrf_in, referent, multi=False, repeat=3, pynio=False): in_wrfnc = [nc for i in xrange(repeat)] refnc = NetCDF(referent) - + if not multi: ref_vals = refnc.variables[varname][:] else: @@ -124,14 +124,14 @@ def make_test(varname, wrf_in, referent, multi=False, repeat=3, pynio=False): nt.assert_allclose(my_vals, ref_vals, tol, atol) elif (varname == "slp"): my_vals = getvar(in_wrfnc, varname, units="hpa") - tol = 2/100. + tol = 1/100. atol = 0 nt.assert_allclose(my_vals, ref_vals, tol, atol) elif (varname == "uvmet"): my_vals = getvar(in_wrfnc, varname) - tol = 1/100. - atol = .5 + tol = 0/100. + atol = 0.0001 nt.assert_allclose(my_vals, ref_vals, tol, atol) elif (varname == "uvmet10"): my_vals = getvar(in_wrfnc, varname) @@ -170,8 +170,133 @@ def make_test(varname, wrf_in, referent, multi=False, repeat=3, pynio=False): return test +def _get_refvals(referent, varname, repeat, multi): + try: + from netCDF4 import Dataset as NetCDF + except: + pass + + refnc = NetCDF(referent) + + if not multi: + ref_vals = refnc.variables[varname][:] + else: + data = refnc.variables[varname][:] + new_dims = [repeat] + [x for x in data.shape] + masked=False + if (isinstance(data, ma.core.MaskedArray)): + masked=True + + if not masked: + ref_vals = n.zeros(new_dims, data.dtype) + else: + ref_vals = ma.asarray(n.zeros(new_dims, data.dtype)) + + for i in xrange(repeat): + ref_vals[i,:] = data[:] + + if masked: + ref_vals.mask[i,:] = data.mask[:] + + return ref_vals + +def make_interp_test(varname, wrf_in, referent, multi=False, + repeat=3, pynio=False): + def test(self): + try: + from netCDF4 import Dataset as NetCDF + except: + pass + + try: + from PyNIO import Nio + except: + pass + + if not multi: + if not pynio: + in_wrfnc = NetCDF(wrf_in) + else: + # Note: Python doesn't like it if you reassign an outer scoped + # variable (wrf_in in this case) + if not wrf_in.endswith(".nc"): + wrf_file = wrf_in + ".nc" + else: + wrf_file = wrf_in + in_wrfnc = Nio.open_file(wrf_file) + else: + if not pynio: + nc = NetCDF(wrf_in) + in_wrfnc = [nc for i in xrange(repeat)] + else: + if not wrf_in.endswith(".nc"): + wrf_file = wrf_in + ".nc" + else: + wrf_file = wrf_in + nc = Nio.open_file(wrf_file) + in_wrfnc = [nc for i in xrange(repeat)] + + if (varname == "interplevel"): + ref_ht_500 = _get_refvals(referent, "z_500", repeat, multi) + hts = getvar(in_wrfnc, "z") + p = getvar(in_wrfnc, "pressure", units="hpa") + hts_500 = interplevel(hts, p, 500) + + nt.assert_allclose(hts_500, ref_ht_500) + + elif (varname == "vertcross"): + ref_ht_cross = _get_refvals(referent, "ht_cross", repeat, multi) + ref_p_cross = _get_refvals(referent, "p_cross", repeat, multi) + + hts = getvar(in_wrfnc, "z") + p = getvar(in_wrfnc, "pressure", units="hpa") + + pivot_point = (hts.shape[-2] / 2, hts.shape[-1] / 2) + ht_cross = vertcross(hts, p, pivot_point=pivot_point,angle=90.) + + nt.assert_allclose(ht_cross, ref_ht_cross, rtol=.01) + + # Test opposite + p_cross1 = vertcross(p,hts,pivot_point=pivot_point, angle=90.0) + + nt.assert_allclose(p_cross1, + ref_p_cross, + rtol=.01) + # Test point to point + start_point = (hts.shape[-2]/2, 0) + end_point = (hts.shape[-2]/2, -1) + + p_cross2 = vertcross(p,hts,start_point=start_point, + end_point=end_point) + + nt.assert_allclose(p_cross1, p_cross2) + + elif (varname == "interpline"): + ref_t2_line = _get_refvals(referent, "t2_line", repeat, multi) + + t2 = getvar(in_wrfnc, "T2") + pivot_point = (t2.shape[-2] / 2, t2.shape[-1] / 2) + + t2_line1 = interpline(t2, pivot_point=pivot_point, angle=90.0) + + nt.assert_allclose(t2_line1, ref_t2_line) + + # Test point to point + start_point = (t2.shape[-2]/2, 0) + end_point = (t2.shape[-2]/2, -1) + + t2_line2 = interpline(t2, start_point=start_point, + end_point=end_point) + + nt.assert_allclose(t2_line1, t2_line2) + + return test + class WRFVarsTest(ut.TestCase): longMessage = True + +class WRFInterpTest(ut.TestCase): + longMessage = True if __name__ == "__main__": @@ -181,6 +306,7 @@ if __name__ == "__main__": "pvo", "pw", "rh2", "rh", "slp", "ter", "td2", "td", "tc", "theta", "tk", "tv", "twb", "updraft_helicity", "ua", "va", "wa", "uvmet10", "uvmet", "z", "ctt", "cape_2d", "cape_3d"] + interp_methods = ["interplevel", "vertcross", "interpline"] try: import netCDF4 @@ -195,6 +321,16 @@ if __name__ == "__main__": test_func2 = make_test(var, TEST_FILE, OUT_NC_FILE, multi=True) setattr(WRFVarsTest, 'test_{0}'.format(var), test_func1) setattr(WRFVarsTest, 'test_multi_{0}'.format(var), test_func2) + + for method in interp_methods: + test_interp_func1 = make_interp_test(method, TEST_FILE, + OUT_NC_FILE) + test_interp_func2 = make_interp_test(method, TEST_FILE, + OUT_NC_FILE, multi=True) + setattr(WRFInterpTest, 'test_{0}'.format(method), + test_interp_func1) + setattr(WRFInterpTest, 'test_multi_{0}'.format(method), + test_interp_func2) try: import PyNIO @@ -211,6 +347,16 @@ if __name__ == "__main__": setattr(WRFVarsTest, 'test_pynio_{0}'.format(var), test_func1) setattr(WRFVarsTest, 'test_pynio_multi_{0}'.format(var), test_func2) + + for method in interp_methods: + test_interp_func1 = make_interp_test(method, TEST_FILE, + OUT_NC_FILE) + test_interp_func2 = make_interp_test(method, TEST_FILE, + OUT_NC_FILE, multi=True) + setattr(WRFInterpTest, 'test_pynio_{0}'.format(method), + test_interp_func1) + setattr(WRFInterpTest, 'test_pynio_multi_{0}'.format(method), + test_interp_func2) ut.main() \ No newline at end of file