关于模式中fortran读.a和.b文件不匹配的问题

taotao0718 · 发表于 2014-12-1 21:14:32

本帖最后由 taotao0718 于 2014-12-2 09:18 编辑

首先，先说一下，我是在学习一个模式，叫HYCOM，如果论坛里有HYCOM达人，那真是最好不过了。

这里 .a和.b文件是网格构造生成的两个地形文件，.a文件里是ASCII码，.b戳开就能看，

.b文件的内容是：
1440 'idm '
  960 'jdm '
-1 'mapflg'
plon:  min,max =    -179.99986    179.99980
plat:  min,max =    -79.00232       89.91841
qlon:  min,max =    -179.99979    179.99986
qlat:  min,max =    -79.02532       89.95042
ulon:  min,max =    -179.99977    179.99992
ulat:  min,max =    -79.00232       89.92188
vlon:  min,max =    -179.99983    179.99965
vlat:  min,max =    -79.02532       89.94513
pang:  min,max =    -179.95934    179.86636
scpx:  min,max =    5120.85742    39728.21484
scpy:  min,max =    5120.86133    39728.18359
scqx:  min,max =    5110.66846    39728.27734
scqy:  min,max =    5110.67236    39728.24609
scux:  min,max =    5120.85791    39728.21875
scuy:  min,max =    5120.86133    39728.18750
scvx:  min,max =    5110.66846    39728.27734
scvy:  min,max =    5110.67236    39728.24609
cori:  min,max = -0.0001427841 0.0001454440
pasp:  min,max =       0.99998       1.00004

.a文件里就是详细的经纬度的分布，可以用matlab的程序读出。

我在构造网格后，检查了.a文件的最大最小值，和.b文件匹配，
但是在最后一步提交作业时，出现了一个错误，如下

reading grid file from regional.grid.[ab]
-1 'mapflg'
plon:  min,max =    -179.99986    179.99980       %%（标注） q1

error - .a and .b files not consistent:
.a,.b min =  -3.400760E+38 -1.799999E+02 -3.400760E+38 %%（标注） q2
.a,.b max = 3.401695E+38  1.799998E+02  3.401695E+38

就是说.a .b的最大最小值不匹配。
因此我查了这一部分的相关F文件
if    (k.eq.1) then
call zaiord(plon, ip,.false., hmina,hmaxa, 9)          %%  这里对应 q1 即 hmina=-179.999986  hmaxa=179.99980  (hmina hmaxa即 .a文件中plon的最大最小值)
......

if    (abs(hmina-hminb).gt.abs(hminb)*1.e-4 .or.
&       abs(hmaxa-hmaxb).gt.abs(hmaxb)*1.e-4    ) then
   if    (mnproc.eq.1) then
         write(lp,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)')
   &    'error - .a and .b files not consistent:',
   &    '.a,.b min = ',hmina,hminb,hmina-hminb,    %% 这里对应q2 即hmina= -3.400760E+38， 可是之前读出来是和q1行对应的-179.999986啊，这里真的搞不明白了
   &    '.a,.b max = ',hmaxa,hmaxb,hmaxa-hmaxb
         endif
      call xcstop('(geopar)')
      stop '(geopar)'
endif

红字部分就是我的疑问了，我对模式还不是很熟悉，想不太明白为什么会出现这种问题，真心请教各位，谢谢！

vvt · 发表于 2014-12-2 08:57:44

第一，我不知道什么 HYCOM，也不知道什么模式。
第二，世界上就2种文件：ASCII文件（也叫文本文件），二进制文件。
没有“十进制”文件的说法。从你的内容来看，b 文件也是ASCII文件。
第三，你的问题比较诡异，需要检查
if (k.eq.1) then
call zaiord(plon, ip,.false., hmina,hmaxa, 9)

这里至

'.a,.b min = ',hmina,hminb,hmina-hminb,

这里之间的代码过程，例如函数传递是否正确、中间是否改变了 hmina 和 hminb。

如果你检查起来有困难，请给出全部代码（如有特殊输入文件，请一并给出，如果HYCOM需要特殊的环境，恐怕就很难办了。这个环境允许单步Debug吗？）

taotao0718 · 发表于 2014-12-2 09:16:10

对很多术语都不是很清楚，.b文件反正就是戳开就能看了的。
大概了解您说的单步debug说的是啥意思了，一定要这么做一定是有办法做到，就是比较麻烦，因为里面还调用了很多别的F文件里的函数
这一段的程序如下：

[Fortran] 纯文本查看 复制代码

 implicit none
 c
       include 'common_blocks.h'
 c
       real      dp0kf,dpm,dpms,ds0kf,dsm,dsms
       real      hmina,hminb,hmaxa,hmaxb
       integer   i,ios,j,k,ktr,l
       character preambl(5)*79,cline*80
#if defined(USE_CCSM3)
       real      plinei(itdm),plinej(jtdm)
       save      plinei,plinej
 #endif
 c
       real       aspmax
       parameter (aspmax=2.0)  ! maximum grid aspect ratio for diffusion
 *     parameter (aspmax=1.0)  ! ignore  grid aspect ratio in  diffusion
 c
 c --- read grid location,spacing,coriolis arrays
 c
       if     (mnproc.eq.1) then  ! .b file from 1st tile only
         write (lp,'(3a)') ' reading grid file from ',
      &                         trim(flnmgrd),'.[ab]'
         open (unit=uoff+9,file=trim(flnmgrd)//'.b',
      &        status='old')
       endif
       call xcsync(flush_lp)
       call zagetc(cline,ios, uoff+9)
       if     (ios.ne.0) then
         if     (mnproc.eq.1) then
           write(lp,'(/ a,i4,i9 /)')
      &      'geopar: I/O error from zagetc, iunit,ios = ',uoff+9,ios
         endif !1st tile
         call xcstop('(geopar)')
                stop '(geopar)'
       endif
       read(cline,*) i
 c
       call zagetc(cline,ios, uoff+9)
       if     (ios.ne.0) then
         if     (mnproc.eq.1) then
           write(lp,'(/ a,i4,i9 /)')
      &      'geopar: I/O error from zagetc, iunit,ios = ',uoff+9,ios
         endif !1st tile
call xcstop('(geopar)')
                stop '(geopar)'
       endif
       read (cline,*) j
 c
       if     (i.ne.itdm .or. j.ne.jtdm) then
         if     (mnproc.eq.1) then
         write(lp,'(/ a /)')
      &    'error - wrong array size in grid file'
         endif
         call xcstop('(geopar)')
                stop '(geopar)'
       endif
       call zagetc(cline,ios, uoff+9)
       if     (ios.ne.0) then
         if     (mnproc.eq.1) then
           write(lp,'(/ a,i4,i9 /)')
      &      'geopar: I/O error from zagetc, iunit,ios = ',uoff+9,ios
         endif !1st tile
         call xcstop('(geopar)')
                stop '(geopar)'
       endif
       if     (mnproc.eq.1) then
       write (lp,'(a)') trim(cline)
       endif
read (cline,*) mapflg
 c
       call zaiopf(trim(flnmgrd)//'.a','old', 9)
 c
       do k= 1,15
         call zagetc(cline,ios, uoff+9)
         if     (ios.ne.0) then
           if     (mnproc.eq.1) then
             write(lp,'(/ a,i4,i9 /)')
      &        'geopar: I/O error from zagetc, iunit,ios = ',uoff+9,ios
           endif !1st tile
           call xcstop('(geopar)')
                  stop '(geopar)'
         endif
         i = index(cline,'=')
         read (cline(i+1:),*) hminb,hmaxb
         if     (mnproc.eq.1) then
         write (lp,'(a)') trim(cline)
         endif
         call xcsync(flush_lp)
 c
         if     (k.eq.1) then
           call zaiord(plon, ip,.false., hmina,hmaxa, 9)
         elseif (k.eq.2) then
           call zaiord(plat, ip,.false., hmina,hmaxa, 9)
           do i= 1,2  !skip qlon,qlat
             call zagetc(cline,ios, uoff+9)
             if     (ios.ne.0) then
               if     (mnproc.eq.1) then
                 write(lp,'(/ a,i4,i9 /)')
      &            'geopar: I/O error from zagetc, iunit,ios = ',
      &            uoff+9,ios
               endif !1st tile
               call xcstop('(geopar)')
                      stop '(geopar)'
             endif
             call zaiosk(9)
           enddo
elseif (k.eq.3) then
           call zaiord(ulon, ip,.false., hmina,hmaxa, 9)
         elseif (k.eq.4) then
           call zaiord(ulat, ip,.false., hmina,hmaxa, 9)
         elseif (k.eq.5) then
           call zaiord(vlon, ip,.false., hmina,hmaxa, 9)
         elseif (k.eq.6) then
           call zaiord(vlat, ip,.false., hmina,hmaxa, 9)
           call zagetc(cline,ios, uoff+9)
           if     (ios.ne.0) then
             if     (mnproc.eq.1) then
               write(lp,'(/ a,i4,i9 /)')
      &          'geopar: I/O error from zagetc, iunit,ios = ',uoff+9,ios
             endif !1st tile
             call xcstop('(geopar)')
                    stop '(geopar)'
           endif
 #if defined(USE_CCSM3)
 c         pang in ANGLET
           i = index(cline,'=')
           read (cline(i+1:),*) hminb,hmaxb
           if     (mnproc.eq.1) then
           write (lp,'(a)') trim(cline)
           endif
           call xcsync(flush_lp)
           call zaiord(ANGLET, ip,.false., hmina,hmaxa, 9)
 #else
 c         skip pang
           call zaiosk(9)
 #endif
         elseif (k.eq.7) then
           call zaiord(scpx, ip,.false., hmina,hmaxa, 9)
         elseif (k.eq.8) then
           call zaiord(scpy, ip,.false., hmina,hmaxa, 9)
         elseif (k.eq.9) then
           call zaiord(scqx, iq,.false., hmina,hmaxa, 9)
         elseif (k.eq.10) then
           call zaiord(scqy, iq,.false., hmina,hmaxa, 9)
         elseif (k.eq.11) then
           call zaiord(scux, iu,.false., hmina,hmaxa, 9)
         elseif (k.eq.12) then
           call zaiord(scuy, iu,.false., hmina,hmaxa, 9)
         elseif (k.eq.13) then
           call zaiord(scvx, iv,.false., hmina,hmaxa, 9)
         elseif (k.eq.14) then
           call zaiord(scvy, iv,.false., hmina,hmaxa, 9)
         else
           call zaiord(corio,iq,.false., hmina,hmaxa, 9)
         endif
 c
         if     (abs(hmina-hminb).gt.abs(hminb)*1.e-4 .or.
      &          abs(hmaxa-hmaxb).gt.abs(hmaxb)*1.e-4     ) then
           if     (mnproc.eq.1) then
           write(lp,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)')
      &      'error - .a and .b files not consistent:',
      &      '.a,.b min = ',hmina,hminb,hmina-hminb,
      &      '.a,.b max = ',hmaxa,hmaxb,hmaxa-hmaxb
           endif
           call xcstop('(geopar)')
                  stop '(geopar)'
         endif
       enddo
 c
       call zaiocl(9)
       if     (mnproc.eq.1) then  ! .b file from 1st tile only
         close(unit=uoff+9)
       endif
 c
       if (itest.gt.0 .and. jtest.gt.0) then
         i=itest
         j=jtest
         write (lp,'(/ a,2i5,a,f8.3,a,f12.9,2f10.2/)')
      &   ' i,j=',i+i0,j+j0,
      &   ' plat=',plat(i,j),
      &   ' corio,scux,vy=',corio(i,j),scux(i,j),scvy(i,j)
       endif
       call xcsync(flush_lp)
 #if defined(USE_CCSM3)
 c --- printout similar to ccsm ice model
       call xclget(plinei,itdm, plon, 1,1, +1, 0, 1)
       call xclget(plinej,jtdm, plat, 1,1,  0,+1, 1)
       if     (mnproc.eq.1) then
         write (lp,*)
         write (lp,'(a,4f9.3,a,4f9.3)')
      &    '(domain) plon(:,1): ',
      &    plinei(1:4),' ...', plinei(itdm-2:itdm)
         write (lp,'(a,4f9.3,a,4f9.3)')
      &    '(domain) plat(1,:): ',
      &    plinej(1:4),' ...', plinej(jtdm-2:jtdm)
         write (lp,*)
       endif
       call xcsync(flush_lp)
 #endif

taotao0718 · 发表于 2014-12-2 09:20:45

vvt 发表于 2014-12-2 08:57
第一，我不知道什么 HYCOM，也不知道什么模式。
第二，世界上就2种文件：ASCII文件（也叫文本文件），二进 ...

全部代码贴上了，但是没选在回复里，请您看看，谢谢

vvt · 发表于 2014-12-2 10:51:31

凡事打开能看懂的，就是 ASCII 文件，也叫文本文件。二进制文件一般是看不懂的。（不要说十进制文件，没有这个说法）

对于你的问题，单步调试是最简单、有效的办法。

求助往往意义不大。因为我没有你的包含文件 include 'common_blocks.h' 等等，没有你的输入数据，甚至我不确定 HYCOM 是否需要特定的环境。

从代码来看，两端之间还有很多次的 call zaiord，这个函数也可能会改变 hmina 和 hminb 的值。

单步debug并不麻烦，点几下鼠标的事儿。如果你确实有困难，可以尝试一步一步把 hmina 和 hminb 都 write 出来，看是哪一步改变了，然后再研究这一步做了什么，该不该做

taotao0718 · 发表于 2014-12-2 11:18:47

vvt 发表于 2014-12-2 10:51
凡事打开能看懂的，就是 ASCII 文件，也叫文本文件。二进制文件一般是看不懂的。（不要说十进制文件，没有 ...

谢谢，认识到这种问题求助也没啥用了，也谢谢您对于文件说法的纠正。

[文件读写] 关于模式中fortran读.a和.b文件不匹配的问题

规矩勋章

元老勋章

新人勋章

水王勋章

热心勋章