problem with gridpak

[kee]

Hi,

I tried to use Gridpak to generate a grid, and simply used the two given examples, namely Carribean and USwest. The first case was successfully duplicated, while I failed in the USwest case.
The error messages are

Code: Select all

[gridpak]$ gridG < Include_USwest/grid.in 
 rectangularity error in mapped contour at iteration  1 is 7.7235E-02
 rectangularity error in mapped contour at iteration  2 is 5.6629E-03
 rectangularity error in mapped contour at iteration  3 is 2.5965E-03
 rectangularity error in mapped contour at iteration  4 is 1.1331E-03
 rectangularity error in mapped contour at iteration  5 is 4.6255E-04
 rectangularity error in mapped contour at iteration  6 is 1.8686E-04
forrtl: severe (408): fort: (3): Subscript #1 of the array B has value 0 which is less than the lower bound of 1

Image              PC                Routine            Line        Source             
libintlc.so.5      00002B36C8232B2D  Unknown               Unknown  Unknown
libintlc.so.5      00002B36C8231635  Unknown               Unknown  Unknown
libifcore.so.5     00002B36C7B2F985  Unknown               Unknown  Unknown
libifcore.so.5     00002B36C7AA0FDD  Unknown               Unknown  Unknown
libifcore.so.5     00002B36C7AA1420  Unknown               Unknown  Unknown
gridG              000000000043EC39  blktr1_                   376  blktri.f90
gridG              000000000043E0E3  blktri_                   302  blktri.f90
gridG              000000000043B1DA  spelip_                   849  sepeli.f90
gridG              0000000000431ACA  sepeli_                   541  sepeli.f90
gridG              000000000040A1C7  MAIN__                    352  grid.f90
gridG              0000000000405B9C  Unknown               Unknown  Unknown
libc.so.6          00000030AE61ECDD  Unknown               Unknown  Unknown
gridG              0000000000405A99  Unknown               Unknown  Unknown

I compiled the program under DEBUG mode, so it is easy to find the position where the program goes wrong. However, I didnot find anything about array B the error messages implied. There is no array B at all.

zhao

[kee]

the grid.h file

Code: Select all

cat grid.h 
#include "griddefs.h"
#include "bathy.h"
      integer         ITMAX, IBIG
      parameter (     ITMAX=6, IBIG=500    )
!  ITMAX is the number of iterations to perform
!  IBIG is the largest number of points to be read in for one
!  boundary.
!
!  original distribution of x,y points is preserved on boundary kb1
!  and kb2:
      integer         kb1, kb2
      parameter (     kb1 = 1, kb2 = 2   )

      integer         L2, M2, L2big, M2big, nwrk
      integer         N, N1, N2, N3, N4
      parameter (     L2=2*(L-1), M2=2*(M-1)   )
      parameter (     L2big=2*Lm, M2big=2*Mm   )
      parameter (     N1=M2, N2=M2+L2, N3=M2+L2+M2, &
     &                N4=M2+L2+M2+L2,  N=N4    )
      integer         KK
      parameter (     KK = 9   )
      parameter (     nwrk = 2*(KK-2)*(2**(KK+1)) + KK + 10*M2big + &
     &                       12*L2big + 27    )
      BIGREAL         sxi(0:L2big), seta(0:M2big)
      common / xiej / sxi, seta
      BIGREAL         x1spl(IBIG),x2spl(IBIG),x3spl(IBIG),x4spl(IBIG), &
     &                y1spl(IBIG),y2spl(IBIG),y3spl(IBIG),y4spl(IBIG), &
     &                s1spl(IBIG),s2spl(IBIG),s3spl(IBIG),s4spl(IBIG), &
     &                b1spl(IBIG),b2spl(IBIG),b3spl(IBIG),b4spl(IBIG), &
     &                c1spl(IBIG),c2spl(IBIG),c3spl(IBIG),c4spl(IBIG)
      integer         nb1pts,nb2pts,nb3pts,nb4pts
      common / bdata/ x1spl, x2spl, x3spl, x4spl, &
     &                y1spl, y2spl, y3spl, y4spl, &
     &                s1spl, s2spl, s3spl, s4spl, &
     &                b1spl, b2spl, b3spl, b4spl, &
     &                c1spl, c2spl, c3spl, c4spl, &
     &                nb1pts, nb2pts, nb3pts, nb4pts
!  The boundary values are read from stdin for edges which have
!  rbnd true.  For boundaries which are read in, the grid spacing
!  is proportional to distance along the boundary if even? is true.
!  Otherwise, it is proportional to the spacing of the supplied
!  boundary points.
      logical         rbnd1, rbnd2, rbnd3, rbnd4, &
     &                even1, even2, even3, even4
      parameter   (  rbnd1=.true., rbnd2=.true., &
     &               rbnd3=.true., rbnd4=.true., &
     &               even1=.true., even2=.true., &
     &               even3=.true., even4=.true.  )

!  The following are used when you need to fit a boundary with
!  bumps on opposite sides and need to make intermediate partial
!  grids.  Set pleft1 or pbot1 to true to print out the boundaries
!  of a partial left or bottom grid.  Set pleft2 or pbot2 to true
!  to print out the new left or bottom boundary.  Lmiddle or Mmiddle
!  gives the position of the interior boundary for the intermediate
!  grid.  The boundaries are written out to iout1 or iout2.
!
!  Don't forget to adjust the evenx flags, kb1 and kb2 accordingly.
      logical         pleft1, pleft2, pbot1, pbot2
      integer         Lmiddle, Mmiddle, iout1, iout2
      parameter   (   pleft1=.false., pleft2=.false., &
     &                pbot1=.false., pbot2=.false., &
     &                Lmiddle=49, Mmiddle=25, &
     &                iout1=13, iout2=14             )

!  These variables are used for writing out a subset of the psi points
!  to be used in generating a nested domain.
      logical        subset
      integer        Lwest, Least, Msouth, Mnorth, iout3
      parameter    ( subset = .false., Lwest = 10, Least = 20, &
     &               Msouth = 10, Mnorth = 20, iout3 = 15    )

!  xpots unit numbers
      integer         ipot1, ipot2, ipot3, ipot4
      parameter   (   ipot1=41, ipot2=42, ipot3=43, ipot4=44   )

the griddefs.h file

Code: Select all

cat griddefs.h 
! define as 1 for ETOPO5 bathymetry
#define ETOPO5        1
! define as 1 for ETOPO2 bathymetry
#undef ETOPO2
#undef GEBCO

! for 64-bit output
#define DBLEPREC      1

! to draw coastlines on some plots
#define DRAW_COASTS   1

! to keep ellipsoidal terms in Earth's shape
#undef ELLIPSOID

! for averaging bathymetry in gridbox (for EW/NS grids only)
#undef IMG_AVG

#define KEEP_SHALLOW 1

! for NCAR graphics (3.2 or better) */
#define PLOTS	     1
! for X windows rather than metafile */
#define X_WIN        1

#undef SYS_POTS       /* unimplimented system calls */
#undef XPOTS1	      /* read ipot1 file */
#undef XPOTS2	      /* read ipot2 file */
#undef XPOTS3	      /* read ipot3 file */
#undef XPOTS4	      /* read ipot4 file */

#ifdef cray
#undef DCOMPLEX
#define DBLEPREC      1	/* for 64-bit output */
#define BIGREAL real
#define SMALLREAL real
#define BIGCOMPLEX complex
#define FLoaT float
#else
#if DBLEPREC
#define DCOMPLEX      1    /* for compilers which support complex*16 */
#define SMALLREAL real
#define BIGREAL real*8
#define BIGCOMPLEX complex*16
#define FLoaT dfloat
#else
#undef DCOMPLEX
#define BIGREAL real
#define SMALLREAL real
#define BIGCOMPLEX complex
#define FLoaT float
#endif  /* DBLEPREC */
#endif  /* cray */

the proj.h file

Code: Select all

cat proj.h 
      character*2   JPRJ, JLTS
      real          PLAT, PLONG, ROTA, P1, P2, P3, P4, XOFF, YOFF
      integer       JGRD
      parameter  (  JPRJ = 'LC', PLAT = 35.500000,&
     &              PLONG = -125.000000, ROTA = 45.500000,&
     &              JLTS = 'CO', P1 = 28.000000,&
     &              P2 = -135.000000, P3 = 50.000000,&
     &              P4 = -112.000000, JGRD = 5)
      parameter  (  XOFF = 0.,  YOFF = 0.)

the gridparame.h

Code: Select all

cat gridparam.h 
      integer         Lm, Mm
      parameter (     Lm=62      , Mm=254      )

[kate]

Indeed, it fails for me with gfortran and ifort, runs with pgi. It got through the hairy "rect" routine only to die in the solver. I don't have a quick fix, sorry.

[kate]

Try this patch:

Code: Select all

diff --git a/Utility/blktri.F b/Utility/blktri.F
index 1d2a23b..e16be34 100644
--- a/Utility/blktri.F
+++ b/Utility/blktri.F
@@ -951,6 +951,7 @@
       common /cblkt/  npp        ,k          ,eps        ,cnv        , &
      &                nm         ,ncmplx     ,ik
       integer         izh        ,ipl        ,id
+      idx = MAX(i,1)
       idp = 0
       if (ir) 107,101,103
   101 if (i-nm) 102,102,107

[kee]

dear Kata,

Thanks for your help!
The compiler I use is ifort, and it did go wrong. I followed what you suggested, added idx = MAX(i,1) into Utility/blktri.F, and now it goes through.
Now I have other questions:
1) how to use bathsuds and bathsoap to get a good r-Value?
---------
In the Caribbean case, ETOPO5 topography data was used, and Min-depth, Max-depth, filtering times without KEEP-SHALLOW were set to 20, 6000, and 20 respectively, followed by a pass of a regular Shapiro filter. The resulting Bottom MAX-slope and Max-RATIO of r-Value are 0.024 and 0.963.
When I use 100 filtering times without KEEP-SHALLOW, the Bottom MAX-slope and Max-RATIO of r-Value are 0.019 and 0.888.
---------
In the USwest case, also ETOPO5 topography data was used, 500 filtering times with KEEP-SHALLOW, and depth1 is 200 in ./Drivers/bathsuds.F. Then a pass of a regular Shapiro filter was applied. The resulting Bottom MAX-slope and Max-RATIO of r-Value are 0.068 and 0.388.
---------
As you can see above, the r-value is large. Is it OK? can it be used to run the model? if not, how to deal with bathsuds and bathsoap to obtain a good r-Value? is there a generic way that not such experiential to do this?
for instance, if using ETOPO2, how many times should be done?

2) there is something wrong when I changed from ETOPO5 to ETOPO2
the dimension of ETOPO5 is

Code: Select all

	lon = 4320 ;
	lat = 2161 ;

while the ETOPO2 is

Code: Select all

 	lon = 10800 ;
	lat = 5400 ;

so I have to modify one place in Drivers/bathtub.F to make it through

Code: Select all

#elif ETOPO2
      parameter     ( evenflag = .true. )
      parameter     ( ILON=360*30+1, ILAT=180*30 )

3) other problems using ETOPO2

grid and tolat seems ok, while not the bathtub.
this case is Caribbean, so the domain is

Code: Select all

 Latitude range    8.73607540130615        30.2309589385986     
 Longitude range   -98.2675476074219       -59.7324600219727 

what's strange is that I saw topography and coastline about other area overlaped this domain(one can have a look at 1.png in the attachment), it is obviously not what I want but there is no error information at the command prompt.

I then checked the difference between dimensions of ETOPO5 and ETOPO2,
in ETOPO5, the latitude is from -90 to 90, longitude is from 0 to 180, then from -180 to 0.
in ETOPO2, the latitude is from -90 to ~89.97 as in ETOPO5, while the longitude is from -180 to 180.
I just wondered that whether the order in zonal of the dataset matters or not.

zhao

[kate]

Hi Kee,

Perhaps these aren't the best examples. It's been years since I used KEEP_SHALLOW.

1. I wouldn't settle for any r-max greater than maybe 0.42.

How many passes of bathsuds depends on your grid resolution, your minimum depth and sometimes never gets there. The resolution of Caribbean is very, very coarse and will only get a good r-max with minimum depths at 100 m or greater. Yes, it's an art, but I like art. I have been known to mask out a problematic point or manually change the depth at a point. The NEP domain was giving us trouble in the Western Bering Sea so we changed the minimum depth linearly from 10 m over the bulk of the domain to 30 m off Kamchatka.

Our newest domain includes Cook Inlet with wetting and drying. That implies a negative minimum depth for which the default bathsuds r-value checker breaks down. Compromises have to be made.

I'm gathering up all the latest versions of gridpak and friends and will most likely move them to github.

2. Yes, my latest has both ETOPO2 and ETOPO1.

3. Yes, you need the longitude to be monotonically increasing. Sometimes my longitudes go -180 to +180, sometimes they go 0 to 360. Sometimes I switch them back and forth for the various needs of bathymetry and winds. Yes, the codes could be more flexible to manage that internally.

ETA: I've now got this code in bathtub.F:

Code: Select all

!  Longitudes between 0 and 360
      lon(ILON) = lon(1)+360.
# if ETOPO5
      do i=1,ILON
        if (lon(i) .lt. 0.) lon(i) = lon(i) + 360.
      enddo
# endif

For ETOPO2, that's when I'd use the code to force the grid into the range of ETOPO2.

[kee]

dear Kate,

I have used roms_tools to prepare input datasets. these days I try to use roms own matlab scripts to do the thing, and met with many problems that still unresovled. I felt terrible these days. I wish you can help me.

<1> I obtain grid.nc using gridpak, then I turn to matlab scripts to edit masks in order to make sure there is no isolated wet points. I did in this way:
firstly: ijcoast('grid.nc', 'Coastline.mat')

Code: Select all

??? Maximum recursion limit of 500 reached. Use set(0,'RecursionLimit',N)
to change the limit.  Be aware that exceeding your available stack space can
crash MATLAB and/or your computer.

Error in ==> ncinfo

questions (1): 'Coastline.mat' used here is globe. Should it be the global one or the one just match the grid area?
questions (2): I really donot know how to deal with the error-"??? Maximum recursion limit of 500 reached.". each time I set N in set(0,'RecursionLimit',N) to a large number, the matlab crashes as it warns.

second: I fail to get ijcoast.dat. Ok, fine. I just donot use it. Instead I use editmask('grid.nc'), and the same error occurs:

Code: Select all

??? Maximum recursion limit of 500 reached. Use set(0,'RecursionLimit',N)
to change the limit.  Be aware that exceeding your available stack space can
crash MATLAB and/or your computer.

Error in ==> ncinfo

<2> I try to use Fotran codes-initial.tar.gz and forcing.tar.gz. However, I meet with the same error when compile these to packeges.

Code: Select all

[dir]$ make
g77 -c -g -DNO_EXIT initial.F
g77 -c -g -DNO_EXIT bleck.F
g77 -c -g -DNO_EXIT crash.F
g77 -c -g -DNO_EXIT day_code.F
g77 -c -g -DNO_EXIT get_date.F
get_date.F: In subroutine `get_date':
get_date.F:111: 
         mon=1      
         ^
Invalid declaration of or reference to symbol `mon' at (^) [initially seen at (^)]
make: *** [get_date.o] Error 1

According to the errors, I opened the get_date.F, and checked the declaration section

Code: Select all

!
      implicit none
      INTEGER_TYPE
     &        dstat, half, hour, iday, imon, len1, len2, len3, min,
     &        nday, sec, tstat, year
      INTEGER_TYPE
     &        lday(31), lmonth(12)
      INTEGER_TYPE
     &        lenstr
#if defined cray
      INTEGER_TYPE
     &        century
      parameter (century=1900)
      character*8  tstring
#elif defined sun || defined sgi || defined DECALPHA
      character*3  day3, mon
      character*28 fdate, tmpday
#elif defined AIX
      character*3  day3, mon
      character*28 tmpday
#endif
      character*3  ampm(0:1)
      character*9  day(0:6), month(12)
      character*11 ctime
      character*18 today
      character*20 fmt
      character*44 wkday
      character*(*) date_str
      data ampm /' AM',' PM'/
      data day /'Sunday','Monday','Tuesday','Wednesday','Thursday',
     &          'Friday','Saturday'/
      data lmonth, lday /7,8,5,5,3,4,4,6,9,7,8,8,9*1,22*2/
      data month /'January','February','March','April','May','June',
     &            'July','August','September','October','November',
     &            'December'/ 

Meantime, in the cppdefs.h, I donot define sun || defined sgi || defined DECALPHA | cray | AIX,

Code: Select all

#undef DOUBLE_PRECISION
/*
**  Undefine double precision flag for default 64-bit computers.
*/
#if defined cray || defined CM2
# ifdef  DOUBLE_PRECISION
#  undef  DOUBLE_PRECISION
# endif
#endif
/*
**  Set double/single precision for real type variables and associated
**  intrinsic functions.
*/
.........

So the error happened because the variable mon has been used before defined.
I donot know how to fix this.

<3> I also find that these FORTRAN codes maybe have been discarded long before, and there is no new Vtransform to select. is it right?

<4> I also use c_grid.m to create a grid successfully. But there is no value for the grid variables, it is an empty grid file. I can fill the lat_rho and lon_rho as I want. However, there are so many variables, I can not do this the same way. I do wish someone can write a guide to use the matlab scripts, even very simple, the step will be helpful.

zhao

[kate]

I'm afraid I can't help you with the Matlab scripts - I'm not using any of them these days. Nor am I using the Fortran codes you mention. If you can fix the Fortran by declaring mon, do so and move on. The pre- and post-processing world is such that the needs change year to year as we change our inputs from NCEP to CORE to MERRA to whatever's next. You need to become comfortable enough with some language that you can modify these tools as needed. If these bits appear, it's because someone wrote a tool and hopes it helps someone else, but they don't get as polished as ROMS itself.

#3. Yes, the old tools don't have the new Vtransform stuff. Feel free to add it yourself. Part of the evolutionary process.

[kee]

ok, thanks!