[ncl-talk] Vertical interpolation to a 2d array instead of a constant pressure/height
Bill Ladwig
ladwig at ucar.edu
Thu Oct 13 09:18:37 MDT 2016
Hi Michael,
The thread is up on ncl-talk now. Glad it's working for you.
Regards,
Bill
On Thu, Oct 13, 2016 at 6:43 AM, Michael Weston <mjweston at masdar.ac.ae>
wrote:
> Hi Bill,
>
> Success!!!
> Indeed it works. Thanks very much for your help!
> Please confirm if I can send our correspondence as a reply to my original
> post on ncl-talk so others can see the solution for posterity.
>
> Thanks and regards
> Michael
>
>
> On 12/10/2016 20:38, Bill Ladwig wrote:
>
> Hi Michael,
>
> You have several problems. First, you removed the output argument from
> your SUBROUTINE declaration.
>
> This:
>
> SUBROUTINE MYINTERP(V3D,Z,SFC,NX,NY,NZ)
>
> Should be this:
>
> SUBROUTINE MYINTERP(V3D,V2D,Z,SFC,NX,NY,NZ)
>
> Otherwise, your data is only getting written to a local variable with no
> way of getting it back to the caller.
>
> Second, in Fortran, SUBROUTINES and FUNCTIONS are not the same thing.
> SUBROUTINES do not return values, FUNCTIONS do. In your NCL code, you are
> trying to use MYINTERP as a FUNCTION, not a SUBROUTINE, so this is why NCL
> is claiming it can't find the FUNCTION.
>
> To fix this, you need to first allocate space for the output by using the
> NCL 'new' function, then call the MYINTERP subroutine. In other words, you
> should make your last three lines this:
>
> result = new((/ny, nx/), float) ; Note, NCL uses C (row major)
> ordering, Fortran uses column major
> IMPORT::MYINTERP(tk, result, p, psf, nx, ny, klev)
> print(result)
>
> This should work now, at least it does on my machine. Good luck.
>
> Bill
>
>
>
> On Wed, Oct 12, 2016 at 2:16 AM, Michael Weston <mjweston at masdar.ac.ae>
> wrote:
>
>> Hi Bill
>>
>> Thanks, yes that is much simpler!
>> My error message now is:
>>
>> *fatal:syntax error: MYINTERP is not a function in package IMPORT*
>>
>> I am not sure what the problem is now as I have checked the examples on
>> lthe WRAPIT page and have tried to follow them exaclty...
>>
>> I do not get a message saying
>>
>> opening: /usr/local/ncl/source_code/ncl_ncarg-6.3.0/ni/src/lib/nfpfort/myinterp.so
>>
>> But I think it is opening the external statement because I recevied
>> warning messages about the type of data not matching that of the fortran
>> code. So, myinterp.f looks like this: C NCLFORTSTART SUBROUTINE
>> MYINTERP(V3D,Z,SFC,NX,NY,NZ) IMPLICIT NONE INTEGER NX,NY,NZ
>> REAL V3D(NX,NY,NZ),V2D(NX,NY) ! This used to be DOUBLE PRECISION
>> REAL Z(NX,NY,NZ) REAL SFC(NX,NY) C NCLEND INTEGER
>> I,J,KP,IP,IM LOGICAL INTERP DOUBLE PRECISION HEIGHT,W1,W2 c
>> does vertical coordinate increase or decrease with increasing k? c set
>> offset appropriately IP = 0 IM = 1 IF
>> (Z(1,1,1).GT.Z(1,1,NZ)) THEN IP = 1 IM = 0 END IF
>> DO I = 1,NX DO J = 1,NY C Initialize to missing. Was
>> initially hard-coded to -999999. INTERP = .false.
>> KP = NZ DO WHILE ((.NOT.INTERP) .AND.
>> (KP.GE.2)) HEIGHT=SFC(I,J) IF
>> (((Z(I,J,KP-IM).LE.HEIGHT).AND. (Z(I,J, +
>> KP-IP).GT.HEIGHT))) THEN W2 = (HEIGHT-Z(I,J,KP-IM))/
>> + (Z(I,J,KP-IP)-Z(I,J,KP-IM))
>> W1 = 1.D0 - W2 V2D(I,J) =
>> W1*V3D(I,J,KP-IM) + W2*V3D(I,J,KP-IP) INTERP = .true.
>> END IF KP = KP - 1 END DO
>> END DO END DO RETURN END And my ncl script
>> looks like this: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" load
>> "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" external IMPORT
>> "/usr/local/ncl/source_code/ncl_ncarg-6.3.0/ni/src/lib/nfpfort/myinterp.so"
>> a = addfile("../wrfout_d02_2016-01-15_06:00:00","r") tk =
>> wrf_user_getvar(a,"tk",1) ; T in Kelvin p =
>> wrf_user_getvar(a,"pressure",1) ; T in Kelvin psfc=
>> wrf_user_getvar(a,"PSFC",1)/100. ; PSFC is in Pa, I want SFC-35mb, so
>> used 3500 Pa. psfc= psfc-35. dimensions = dimsizes(tk)
>> print(dimensions) printVarSummary(psfc) printMinMax(psfc,True) klev =
>> dimensions(0) ny = dimensions(1) nx = dimensions(2) print(typeof(tk))
>> print(typeof(p)) print(typeof(psfc)) t1 =
>> IMPORT::MYINTERP(tk,p,psfc,nx,ny,klev) print(t1) Any assistance would
>> be appreciated. Regards Michael
>> On 11/10/2016 19:39, Bill Ladwig wrote:
>>
>> Hi Michael,
>> I should have explained this better. Take your version of DINTERP3DZ and
>> change the name of the routine to something unique, like MYINTERP, and
>> place this in its own .f file, like myinterp.f. Then, only use WRAPIT on
>> this new fortran file. You're trying to make a new routine here, not
>> replace the one in NCL. WRAPIT will generate a new shared library file for
>> it. You can then use this new routine by adding
>> external MYINTERP "/path/to/myinterp.so"
>> to the top of your script.
>> If you don't assign a new name to the routine, you're going to have name
>> conflict problems with the one included in NCL. I would try to get the
>> WRAPIT solution working first, because building NCL from source can be very
>> difficult.
>> Bill
>> On Tue, Oct 11, 2016 at 2:01 AM, Michael Weston <mjweston at masdar.ac.ae>
>> wrote:
>>>
>>> Hi Bill Thank you, this is exactly what I was looking for. Modifying the
>>> source subroutine for DINTERP3DZ was simple enough (I edited it in
>>> wrf_user.f but attached as dinterp3dz.f fyi). I followed the WRAPIT
>>> instructions and then tried to run with the precompiled binary NCL. i.e.
>>> edited ni/src/lib/nfp/wrfW.c to include a new 2d double point precision
>>> variable as input to the function. diff wrfW.bak wrfW.c 17c17
>>> < double *,int *,int *,
>>> int*, --- > double *, double
>>> *,int *,int *, int*, 1918c1918 < NGCALLF(dinterp3dz,DINTERP3DZ)
>>> (tmp_v3d,tmp_v2d,tmp_z,tmp_loc, --- > NGCALLF(dinterp3dz,DINTERP3DZ)
>>> (tmp_v3d,tmp_v2d,tmp_z,tmp_loc,tmp_sfc, I then ran WRAPIT (from
>>> precompiled binary) WRAPIT $NCL_SOURCE/ni/src/lib/nfp/wrfW.c
>>> $NCL_SOURCE/ni/src/lib/nfpfortran/wrf_user.f to create wrfW.so I then
>>> tried to call it as suggested in an ncl script using (wrf_t1.ncl attached)
>>> summary below: external INTERP "/usr/local/ncl/source_code/nc
>>> l_ncarg-6.3.0/ni/src/lib/nfp/wrfW.so" begin
>>> INTERP::dinterp3dz(tk,p,psfc) end gives a warning warning:Could not find
>>> Init() in external file...file not loaded fatal:syntax error: line 16 in
>>> file ./wrf_t1.ncl before or near : EX01: ------^ fatal:error in statement
>>> So I tried using the environment route. export
>>> NCL_DEF_LIB_DIR=/path/to/shared_objects and removed the double EX01::
>>> fatal:syntax error: line 16 in file ./wrf_t1.ncl before or near \n
>>> dinterp3dz(tk,p,psfc) -----------------------^ fatal:syntax error: possibly
>>> an undefined procedure The next trick is to compile NCL form source...and
>>> retry the above. Regards Michael
>>> On 10/10/2016 20:08, Bill Ladwig wrote:
>>>
>>> Hi Michael,
>>> You have found the NCL and C wrappers, but the Fortran function is
>>> located at: ni/src/lib/nfpfort/wrf_user.f and the routine name is
>>> DINTERP1D. You can copy and modify this routine to add the X, Y loops and
>>> then wrap it using the wrapit script (see http://www.ncl.ucar.edu/Document/Tools/WRAPIT.shtml).
>>> Or, better yet, you could copy and modify DINTERP3DZ to take your surface
>>> pressure field, as it's mostly set up to do what you want. You just need
>>> to change "HEIGHT" to grab from the 2D array inside the I,J loop instead of
>>> using the constant "LOC" value.
>>> Best of luck,
>>> Bill
>>> On Mon, Oct 10, 2016 at 3:39 AM, Michael Weston <mjweston at masdar.ac.ae>
>>> wrote:
>>>>
>>>> Hi Bill, Thanks for the reply. I hadn't considered your approach yet.
>>>> Unfortunately the time required in looping through all points is prhibitive
>>>> and I would rather modify code to use an array approach. Do you know where
>>>> the built in functions are defined in ncl? e.g. wrf_interp_1d is called
>>>> from WRFUserARW.ncl, but is not defined in the script. I downloaded the
>>>> source code and did a grep on ALL the files; grep -i -H wrf_interp_1d
>>>> source_files. It is not defined anywhere except for error handling in
>>>> ./ni/src/lib/nfp/wrfW.c and ./ni/src/lib/nfp/wrapper.c. Perhaps I am
>>>> missing something here? Regards Michael
>>>> On 05/10/2016 21:01, Bill Ladwig wrote:
>>>>
>>>> Hi Michael,
>>>> There is nothing in the WRF routines that will do this out of the box,
>>>> but you can manually construct this yourself using wrf_interp_1d. See
>>>> http://www.ncl.ucar.edu/Document/Functions/Built-in/wrf_interp_1d.shtml.
>>>> Assuming I'm completely understanding what you are trying to do, you would
>>>> have to iterate over the nx x ny horizontal grid, then supply the vertical
>>>> column values for each grid point and run the wrf_interp_1d routine to get
>>>> the interpolated value. So, at each horizontal grid point, v_in would be
>>>> the vertical column for your variable of interest, z_in would be the
>>>> vertical coordinates for the v_in column, and z_out would be the desired
>>>> surface pressure value (as an array with one element).
>>>> This will probably take a while to run, so plan accordingly.
>>>> Hope this helps,
>>>> Bill
>>>> On Tue, Oct 4, 2016 at 11:59 PM, Michael Weston <mjweston at masdar.ac.ae>
>>>> wrote:
>>>>>
>>>>> Folks, I am looking for an ncl function that can do vertical
>>>>> interpolation in 3d, but instead of using a constant level to interpolate
>>>>> to, I want to use an array variable to interpolate to. e.g. from existing
>>>>> function wrf_interp_3d_z
>>>>>
>>>>> function wrf_interp_3d_z (
>>>>> v3d : numeric,
>>>>> vert : numeric,
>>>>> loc [1] : numeric
>>>>> )
>>>>>
>>>>> loc is a constant pressure or constant height. I want to use "surface
>>>>> pressure -x" to get variables at the top of the surface layer, so loc would
>>>>> be a 2d array with ny x nx. Does such a function exist? Thanks
>>>>> --
>>>>> *Michael* *Weston*
>>>>> _______________________________________________ ncl-talk mailing list
>>>>> ncl-talk at ucar.edu List instructions, subscriber options, unsubscribe:
>>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>>
>>>> --
>>>> *Michael* *Weston* Research Engineer
>>>> PO Box 54224, Abu Dhabi, United Arab Emirates Office +971 2 810 9510
>>>> <%2B971%202%20810%209510>
>>>> Email mjweston at masdar.ac.ae
>>>> <http://www.masdar.ac.ae/>http://www.masdar.ac.ae
>>>>
>>>>
>>>> P *Please consider the environment before printing this email* This
>>>> transmission is confidential and intended solely for the person or
>>>> organization to whom it is addressed. It may contain privileged and
>>>> confidential information. If you are not the intended recipient, you should
>>>> not copy, distribute or take any action in reliance on it. If you have
>>>> received this transmission in error, please notify us immediately by e-mail
>>>> at *info at masdar.ae <info at masdar.ae>**.*
>>>>
>>> --
>>> *Michael* *Weston* Research Engineer
>>> PO Box 54224, Abu Dhabi, United Arab Emirates Office +971 2 810 9510
>>> <%2B971%202%20810%209510>
>>> Email mjweston at masdar.ac.ae
>>> <http://www.masdar.ac.ae/>http://www.masdar.ac.ae
>>>
>>>
>>> P *Please consider the environment before printing this email* This
>>> transmission is confidential and intended solely for the person or
>>> organization to whom it is addressed. It may contain privileged and
>>> confidential information. If you are not the intended recipient, you should
>>> not copy, distribute or take any action in reliance on it. If you have
>>> received this transmission in error, please notify us immediately by e-mail
>>> at *info at masdar.ae <info at masdar.ae>**.*
>>>
>> --
>> *Michael* *Weston* Research Engineer
>> PO Box 54224, Abu Dhabi, United Arab Emirates Office +971 2 810 9510
>> Email mjweston at masdar.ac.ae
>> <http://www.masdar.ac.ae/>http://www.masdar.ac.ae
>>
>>
>> P *Please consider the environment before printing this email* This
>> transmission is confidential and intended solely for the person or
>> organization to whom it is addressed. It may contain privileged and
>> confidential information. If you are not the intended recipient, you should
>> not copy, distribute or take any action in reliance on it. If you have
>> received this transmission in error, please notify us immediately by e-mail
>> at *info at masdar.ae <info at masdar.ae>**.*
>>
> --
> *Michael* *Weston* Research Engineer
> PO Box 54224, Abu Dhabi, United Arab Emirates Office +971 2 810 9510
> Email mjweston at masdar.ac.ae
> http://www.masdar.ac.ae
>
>
> P *Please consider the environment before printing this email* This
> transmission is confidential and intended solely for the person or
> organization to whom it is addressed. It may contain privileged and
> confidential information. If you are not the intended recipient, you should
> not copy, distribute or take any action in reliance on it. If you have
> received this transmission in error, please notify us immediately by e-mail
> at *info at masdar.ae <info at masdar.ae>**.*
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20161013/9be628de/attachment.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: image/jpeg
Size: 5411 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20161013/9be628de/attachment.jpe
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: image/jpeg
Size: 5411 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20161013/9be628de/attachment-0001.jpe
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: image/jpeg
Size: 5411 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20161013/9be628de/attachment-0002.jpe
-------------- next part --------------
A non-text attachment was scrubbed...
Name: masdarlogo.jpg
Type: image/jpeg
Size: 5411 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20161013/9be628de/attachment.jpg
More information about the ncl-talk
mailing list