[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