Not sure of the problem. Your wrapped code and the NCL script look good.
Some assorted comments:
[1]
The fortran documentation states:
F must be dimensioned at least (M+1)*(N+1).
Are the arguments the correct sizes?
--- [2] I know that NCL is distributed with LAPACK. However, I did not know NCL was distributed with FISHPACK. The LAPACK is double precision so if NCL has FISHPACK, it is likely double precision also. You are passing real values and not double. This could cause the type of error you are getting. Please try double precision. +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ The following is text extracted from 2 power point slides discussed in the NCL Workshop ================================================ Accessing LAPACK (1 of 2) double precision LAPACK (BLAS) => distributed with NCL explicitly link LAPACK lib via fortran interface: WRAPIT eg: subroutine dgels solves [over/under]determined real linear systems C NCLFORTSTART SUBROUTINE DGELSI( M, N, NRHS, A, B, LWORK, WORK ) IMPLICIT NONE INTEGER M, N, NRHS, LWORK DOUBLE PRECISION A( M, N ), B( M, NRHS), WORK(LWORK) C NCLEND C declare local variables INTEGER INFO CHARACTER*1 TRANS TRANS = "N” CALL DGELS(TRANS, M,N,NRHS,A,LDA,B,LDB,WORK,LWORK,INFO) RETURN END ================================================ Accessing LAPACK (2 of 2) external DGELS "./dgels_interface.so” ; NAG example: http://www.nag.com/lapack-ex/node45.html ; These are transposed from the fortran example A = (/ (/ -0.57, -1.93, 2.30, -1.93, 0.15, -0.02 /), \ ; (4,6) (/ -1.28, 1.08, 0.24, 0.64, 0.30, 1.03 /), \ (/ -0.39, -0.31, 0.40, -0.66, 0.15, -1.43 /), \ (/ 0.25, -2.14,-0.35, 0.08,-2.13, 0.50 /) /)*1d0 ; must be double dimA = dimsizes(A) N = dimA(0) ; 4 M = dimA(1) ; 6 B = (/-2.67 ,-0.55 ,3.34, -0.77, 0.48, 4.10/)*1d0 ; must be double ; LAPACK wants 2D nrhs = 1 B2 = conform_dims ( (/nrhs,M/), B, 1 ) ; (1,6) B2(0,:) =B lwork=500 ; allocate space work = new ( lwork, "double", "No_FillValue") ; must be double DGELS::dgelsiw(M, N, nrhs, A , B2, lwork, work ) print(B2(0,0:N-1)) --- Hope this helps. On 10/23/13 9:39 AM, Leung YU Ting wrote: > HI > > When I want to use the HWSSSP of fishpack Fortran subroutine in NCL. I also add a subroutine hstub.f to define the dimension of variables as shown below; > > C NCLFORTSTART > SUBROUTINE HWSSSPDRV (TS,TF,M,MBDCND,BDTS,BDTF,PS,PF,N,NBDCND, > 1 BDPS,BDPF,ELMBDA,F,NLON,NLAT,PERTRB, > 2 IERROR,W,LWORK) > DIMENSION F(NLON,NLAT) ,BDTS(NLON),BDTF(NLON),BDPS(NLAT), > 1 BDPF(NLAT),W(LWORK) > C NCLEND > print*, "1" > CALL HWSSSP(TS,TF,M,MBDCND,BDTS,BDTF,PS,PF,N,NBDCND,BDPS, > 1 BDPF,ELMBDA,F,NLON,PERTRB,IERROR,W) > print*, "2" > RETURN > print*, "3" > END > > The ncl script used is shown below; > > rawtime=hdata->time > time=ut_calendar(rawtime, -3) > lat=hdata->lat(::-1) > lon=hdata->lon > > wantedtime=2000120100 > > timeind=ind(wantedtime .eq. time) > > wantedhgt1=short2flt(hdata->hgt(timeind-1, {500}, ::-1, :)) > wantedhgt2=short2flt(hdata->hgt(timeind+1, {500}, ::-1, :)) > > guv1=z2geouv(wantedhgt1, lat, lon, 1) > gu1=guv1(0, :, :) > gv1=guv1(1, :, :) > gvort1=uv2vr_cfd(gu1, gv1, lat, lon, 3) > copy_VarCoords(wantedhgt1, gvort1) > > guv2=z2geouv(wantedhgt2, lat, lon, 1) > gu2=guv2(0, :, :) > gv2=guv2(1, :, :) > gvort2=uv2vr_cfd(gu2, gv2, lat, lon, 3) > copy_VarCoords(wantedhgt2, gvort2) > > gvort_tend=gvort2-gvort1 > copy_VarCoords(gvort1, gvort_tend) > > PI = 3.14159265358979 > NLAT = dimsizes(lat) > NLON = dimsizes(lon) > PERTRB = 0. > IERROR = 1 > TS= 0. > TF= PI > M=NLAT-1 > MBDCND=9 > BDTS=new(NLON, float) > BDTF=new(NLON, float) > PS=0. > PF=2.*PI*357.5/360. > N=NLON-1 > M=NLAT-1 > NBDCND=0 > BDPS=new(NLAT, float) > BDPF=new(NLAT, float) > ELMBDA=0. > IDIMF= NLON > LWORK = 11*(M+1)+6*(N+1) > W=new(LWORK,float) > F=gvort_tend > > HSTUB::HWSSSPDRV(TS,TF,M,MBDCND,BDTS,BDTF,PS,PF,N,NBDCND,BDPS,BDPF,ELMBDA,F,NLON,NLAT,PERTRB,IERROR,W,LWORK) > > The output is; > 1 > 2 > Segmentation fault (core dumped) > > It seems the script stop at return the variables. So I use gdb to check the script. It is stop at _NclGetObj. What can I do to solve this problem? > > Program received signal SIGSEGV, Segmentation fault. > 0x000000000061d5f2 in _NclGetObj () > (gdb) where > #0 0x000000000061d5f2 in _NclGetObj () > #1 0x000000000064828f in VarValRep () > #2 0x000000000064de18 in _NclGetVarRepValue () > #3 0x0000000000649689 in VarVarWrite () > #4 0x000000000064e2d4 in _NclAssignVarToVar () > #5 0x00000000006cd4f9 in _NclRemapIntrParameters () > #6 0x00000000006c2070 in CallINTRINSIC_PROC_CALL () > #7 0x00000000006c91ad in _NclExecute () > #8 0x00000000005af1b1 in yyparse () > #9 0x00000000005aa592 in main () > (gdb) > > Regards > Marco > > > > > _______________________________________________ > ncl-talk mailing list > List instructions, subscriber options, unsubscribe: > http://mailman.ucar.edu/mailman/listinfo/ncl-talk > _______________________________________________ ncl-talk mailing list List instructions, subscriber options, unsubscribe: http://mailman.ucar.edu/mailman/listinfo/ncl-talkReceived on Thu Oct 24 21:45:49 2013
This archive was generated by hypermail 2.1.8 : Fri Nov 01 2013 - 08:58:14 MDT