注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

Mr.Right

不顾一切的去想,于是我们有了梦想。脚踏实地的去做,于是梦想成了现实。

 
 
 

日志

 
 
关于我

人生一年又一年,只要每年都有所积累,有所成长,都有那么一次自己认为满意的花开时刻就好。即使一时不顺,也要敞开胸怀。生命的荣枯并不是简单的重复,一时的得失不是成败的尺度。花开不是荣耀,而是一个美丽的结束,花谢也不是耻辱,而是一个低调的开始。

网易考拉推荐

阿英讲gfortran编译运行GMF(Global Mapping Functions)的Fortran 77程序  

2016-12-19 21:34:52|  分类: 编程 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
缘起:最近在研究PPP的方法,需要对流层的GMF函数。但是是Fortran的,自己之前也没弄过。摸索了两个小时才调试通过,贴出来,方便后来人。阿弥陀佛!

1) 运行环境: MSYS + gfortran

2) 源代码及运行结果

gfortran helloworld.f -o helloworld
等价地
gcc helloworld.f -o helloworld -lgfortran -lgfortranbegin
库文件 libgfortranbegin.a (通过命令行选项 -lgfortranbegin 被调用) 包含运行和终止一个 Fortran 程序所必须的开始和退出代码。
库文件 libgfortran.a 包含 Fortran 底层的输入输出等所需要的运行函数。当运行 gfortran 时,会自动链接这两个库。



C     gfortran gmf.f -o gmf
C     gcc gmf.f -o gmf -lgfortran -lgfortranbegin
      PROGRAM MAIN
      DOUBLE PRECISION DMJD,DLAT,DLON,DHGT,ZD,GMFH,GMFW
      DMJD= 55055D0
      DLAT=0.6708665767D0
      DLON=-1.393397187D0
      DHGT=844.715D0
      ZD=1.278564131D0
      CALL GMF (DMJD,DLAT,DLON,DHGT,ZD,GMFH,GMFW)
      WRITE(*,*)GMFH
      WRITE(*,*)GMFW
      PRINT *, 'GMFH=', GMFH
      PRINT *, 'GMFW=', GMFW
      END PROGRAM MAIN
      SUBROUTINE GMF (DMJD,DLAT,DLON,DHGT,ZD,GMFH,GMFW)
*+
*  - - - - - - - - -
*   G M F
*  - - - - - - - - -
*
*  This routine is part of the International Earth Rotation and
*  Reference Systems Service (IERS) Conventions software collection.
*
*  This subroutine determines the Global Mapping Functions GMF (Boehm et al. 2006).
*
*  In general, Class 1, 2, and 3 models represent physical effects that
*  act on geodetic parameters while canonical models provide lower-level
*  representations or basic computations that are used by Class 1, 2, or
*  3 models.
*
*  Status: Class 1 model 
*
*     Class 1 models are those recommended to be used a priori in the
*     reduction of raw space geodetic data in order to determine
*     geodetic parameter estimates.
*     Class 2 models are those that eliminate an observational
*     singularity and are purely conventional in nature.
*     Class 3 models are those that are not required as either Class
*     1 or 2.
*     Canonical models are accepted as is and cannot be classified as a
*     Class 1, 2, or 3 model.
*
*  Given:
*     DMJD           d      Modified Julian Date
*     DLAT           d      Latitude given in radians (North Latitude)
*     DLON           d      Longitude given in radians (East Longitude)
*     DHGT           d      Height in meters (mean sea level)
*     ZD             d      Zenith distance in radians
*
*  Returned:
*     GMFH           d      Hydrostatic mapping function (Note 1)
*     GMFW           d      Wet mapping function (Note 1)
*
*  Notes:
*
*  1) The mapping functions are dimensionless scale factors.
*
*  2) This is from a 9x9 Earth Gravitational Model (EGM).
*
*  Test case:
*     given input: DMJD = 55055D0
*                  DLAT = 0.6708665767D0 radians (NRAO, Green Bank, WV)
*                  DLON = -1.393397187D0 radians
*                  DHGT = 844.715D0 meters
*                  ZD   = 1.278564131D0 radians
*
*     expected output: GMFH = 3.425245519339138678D0
*                      GMFW = 3.449589116182419257D0
*                    
*  References:
*
*     Boehm, J., Niell, A., Tregoning, P. and Schuh, H., (2006),
*     "Global Mapping Functions (GMF): A new empirical mapping
*     function based on numerical weather model data",
*     Geophy. Res. Lett., Vol. 33, L07304, doi:10.1029/2005GL025545.
*
*     Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
*     IERS Technical Note No. 36, BKG (2010)
*
*  Revisions:
*  2005 August 30 J. Boehm    Original code
*  2009 August 11 B.E. Stetzler Added header and copyright
*  2009 August 12 B.E. Stetzler More modifications and defined twopi
*  2009 August 12 B.E. Stetzler Provided test case
*  2009 August 12 B.E. Stetzler Capitalized all variables for FORTRAN 77
*                              compatibility and corrected test case
*                              latitude and longitude coordinates
*-----------------------------------------------------------------------
      IMPLICIT NONE
      DOUBLE PRECISION DMJD, DLAT, DLON, DHGT, ZD, GMFH, GMFW
      DOUBLE PRECISION DFAC(20),P(10,10),AP(55),BP(55),
     .          AH_MEAN(55),BH_MEAN(55),AH_AMP(55),BH_AMP(55),
     .          AW_MEAN(55),BW_MEAN(55),AW_AMP(55),BW_AMP(55)
      INTEGER I, J, K, M, N, IR  
      DOUBLE PRECISION TWOPI, DOY, T, SUM1, BH, C0H, PHH, C11H, C10H,
     .                 CH, AHM, AHA, AH, SINE, BETA, GAMMA, TOPCON,
     .                 A_HT, B_HT, C_HT, HS_KM, HT_CORR_COEF, HT_CORR,
     .                 BW, CW, AWM, AWA, AW, PI
      PARAMETER ( PI = 3.1415926535897932384626433D0 )
      PARAMETER (TWOPI = 6.283185307179586476925287D0)
*+---------------------------------------------------------------------
*     Reference day is 28 January 1980
*     This is taken from Niell (1996) to be consistent
*----------------------------------------------------------------------
      DOY = DMJD  - 44239D0 + 1 - 28
      DATA (AH_MEAN(I),I=1,55)/
     .+1.2517D+02, +8.503D-01, +6.936D-02, -6.760D+00, +1.771D-01,
     . +1.130D-02, +5.963D-01, +1.808D-02, +2.801D-03, -1.414D-03,
     . -1.212D+00, +9.300D-02, +3.683D-03, +1.095D-03, +4.671D-05,
     . +3.959D-01, -3.867D-02, +5.413D-03, -5.289D-04, +3.229D-04,
     . +2.067D-05, +3.000D-01, +2.031D-02, +5.900D-03, +4.573D-04,
     . -7.619D-05, +2.327D-06, +3.845D-06, +1.182D-01, +1.158D-02,
     . +5.445D-03, +6.219D-05, +4.204D-06, -2.093D-06, +1.540D-07,
     . -4.280D-08, -4.751D-01, -3.490D-02, +1.758D-03, +4.019D-04,
     . -2.799D-06, -1.287D-06, +5.468D-07, +7.580D-08, -6.300D-09,
     . -1.160D-01, +8.301D-03, +8.771D-04, +9.955D-05, -1.718D-06,
     . -2.012D-06, +1.170D-08, +1.790D-08, -1.300D-09, +1.000D-10/
      DATA (BH_MEAN(I),I=1,55)/
     . +0.000D+00, +0.000D+00, +3.249D-02, +0.000D+00, +3.324D-02,
     . +1.850D-02, +0.000D+00, -1.115D-01, +2.519D-02, +4.923D-03,
     . +0.000D+00, +2.737D-02, +1.595D-02, -7.332D-04, +1.933D-04,
     . +0.000D+00, -4.796D-02, +6.381D-03, -1.599D-04, -3.685D-04,
     . +1.815D-05, +0.000D+00, +7.033D-02, +2.426D-03, -1.111D-03,
     . -1.357D-04, -7.828D-06, +2.547D-06, +0.000D+00, +5.779D-03,
     . +3.133D-03, -5.312D-04, -2.028D-05, +2.323D-07, -9.100D-08,
     . -1.650D-08, +0.000D+00, +3.688D-02, -8.638D-04, -8.514D-05,
     . -2.828D-05, +5.403D-07, +4.390D-07, +1.350D-08, +1.800D-09,
     . +0.000D+00, -2.736D-02, -2.977D-04, +8.113D-05, +2.329D-07,
     . +8.451D-07, +4.490D-08, -8.100D-09, -1.500D-09, +2.000D-10/
      
      DATA (AH_AMP(I),I=1,55)/
     . -2.738D-01, -2.837D+00, +1.298D-02, -3.588D-01, +2.413D-02,
     . +3.427D-02, -7.624D-01, +7.272D-02, +2.160D-02, -3.385D-03,
     . +4.424D-01, +3.722D-02, +2.195D-02, -1.503D-03, +2.426D-04,
     . +3.013D-01, +5.762D-02, +1.019D-02, -4.476D-04, +6.790D-05,
     . +3.227D-05, +3.123D-01, -3.535D-02, +4.840D-03, +3.025D-06,
     . -4.363D-05, +2.854D-07, -1.286D-06, -6.725D-01, -3.730D-02,
     . +8.964D-04, +1.399D-04, -3.990D-06, +7.431D-06, -2.796D-07,
     . -1.601D-07, +4.068D-02, -1.352D-02, +7.282D-04, +9.594D-05,
     . +2.070D-06, -9.620D-08, -2.742D-07, -6.370D-08, -6.300D-09,
     . +8.625D-02, -5.971D-03, +4.705D-04, +2.335D-05, +4.226D-06,
     . +2.475D-07, -8.850D-08, -3.600D-08, -2.900D-09, +0.000D+00/
      
      DATA (BH_AMP(I),I=1,55)/
     . +0.000D+00, +0.000D+00, -1.136D-01, +0.000D+00, -1.868D-01,
     . -1.399D-02, +0.000D+00, -1.043D-01, +1.175D-02, -2.240D-03,
     . +0.000D+00, -3.222D-02, +1.333D-02, -2.647D-03, -2.316D-05,
     . +0.000D+00, +5.339D-02, +1.107D-02, -3.116D-03, -1.079D-04,
     . -1.299D-05, +0.000D+00, +4.861D-03, +8.891D-03, -6.448D-04,
     . -1.279D-05, +6.358D-06, -1.417D-07, +0.000D+00, +3.041D-02,
     . +1.150D-03, -8.743D-04, -2.781D-05, +6.367D-07, -1.140D-08,
     . -4.200D-08, +0.000D+00, -2.982D-02, -3.000D-03, +1.394D-05,
     . -3.290D-05, -1.705D-07, +7.440D-08, +2.720D-08, -6.600D-09,
     . +0.000D+00, +1.236D-02, -9.981D-04, -3.792D-05, -1.355D-05,
     . +1.162D-06, -1.789D-07, +1.470D-08, -2.400D-09, -4.000D-10/
      
      DATA (AW_MEAN(I),I=1,55)/
     . +5.640D+01, +1.555D+00, -1.011D+00, -3.975D+00, +3.171D-02,
     . +1.065D-01, +6.175D-01, +1.376D-01, +4.229D-02, +3.028D-03,
     . +1.688D+00, -1.692D-01, +5.478D-02, +2.473D-02, +6.059D-04,
     . +2.278D+00, +6.614D-03, -3.505D-04, -6.697D-03, +8.402D-04,
     . +7.033D-04, -3.236D+00, +2.184D-01, -4.611D-02, -1.613D-02,
     . -1.604D-03, +5.420D-05, +7.922D-05, -2.711D-01, -4.406D-01,
     . -3.376D-02, -2.801D-03, -4.090D-04, -2.056D-05, +6.894D-06,
     . +2.317D-06, +1.941D+00, -2.562D-01, +1.598D-02, +5.449D-03,
     . +3.544D-04, +1.148D-05, +7.503D-06, -5.667D-07, -3.660D-08,
     . +8.683D-01, -5.931D-02, -1.864D-03, -1.277D-04, +2.029D-04,
     . +1.269D-05, +1.629D-06, +9.660D-08, -1.015D-07, -5.000D-10/
      
      DATA (BW_MEAN(I),I=1,55)/
     . +0.000D+00, +0.000D+00, +2.592D-01, +0.000D+00, +2.974D-02,
     . -5.471D-01, +0.000D+00, -5.926D-01, -1.030D-01, -1.567D-02,
     . +0.000D+00, +1.710D-01, +9.025D-02, +2.689D-02, +2.243D-03,
     . +0.000D+00, +3.439D-01, +2.402D-02, +5.410D-03, +1.601D-03,
     . +9.669D-05, +0.000D+00, +9.502D-02, -3.063D-02, -1.055D-03,
     . -1.067D-04, -1.130D-04, +2.124D-05, +0.000D+00, -3.129D-01,
     . +8.463D-03, +2.253D-04, +7.413D-05, -9.376D-05, -1.606D-06,
     . +2.060D-06, +0.000D+00, +2.739D-01, +1.167D-03, -2.246D-05,
     . -1.287D-04, -2.438D-05, -7.561D-07, +1.158D-06, +4.950D-08,
     . +0.000D+00, -1.344D-01, +5.342D-03, +3.775D-04, -6.756D-05,
     . -1.686D-06, -1.184D-06, +2.768D-07, +2.730D-08, +5.700D-09/
      
      DATA (AW_AMP(I),I=1,55)/
     . +1.023D-01, -2.695D+00, +3.417D-01, -1.405D-01, +3.175D-01,
     . +2.116D-01, +3.536D+00, -1.505D-01, -1.660D-02, +2.967D-02,
     . +3.819D-01, -1.695D-01, -7.444D-02, +7.409D-03, -6.262D-03,
     . -1.836D+00, -1.759D-02, -6.256D-02, -2.371D-03, +7.947D-04,
     . +1.501D-04, -8.603D-01, -1.360D-01, -3.629D-02, -3.706D-03,
     . -2.976D-04, +1.857D-05, +3.021D-05, +2.248D+00, -1.178D-01,
     . +1.255D-02, +1.134D-03, -2.161D-04, -5.817D-06, +8.836D-07,
     . -1.769D-07, +7.313D-01, -1.188D-01, +1.145D-02, +1.011D-03,
     . +1.083D-04, +2.570D-06, -2.140D-06, -5.710D-08, +2.000D-08,
     . -1.632D+00, -6.948D-03, -3.893D-03, +8.592D-04, +7.577D-05,
     . +4.539D-06, -3.852D-07, -2.213D-07, -1.370D-08, +5.800D-09/
      
      DATA (BW_AMP(I),I=1,55)/
     . +0.000D+00, +0.000D+00, -8.865D-02, +0.000D+00, -4.309D-01,
     . +6.340D-02, +0.000D+00, +1.162D-01, +6.176D-02, -4.234D-03,
     . +0.000D+00, +2.530D-01, +4.017D-02, -6.204D-03, +4.977D-03,
     . +0.000D+00, -1.737D-01, -5.638D-03, +1.488D-04, +4.857D-04,
     . -1.809D-04, +0.000D+00, -1.514D-01, -1.685D-02, +5.333D-03,
     . -7.611D-05, +2.394D-05, +8.195D-06, +0.000D+00, +9.326D-02,
     . -1.275D-02, -3.071D-04, +5.374D-05, -3.391D-05, -7.436D-06,
     . +6.747D-07, +0.000D+00, -8.637D-02, -3.807D-03, -6.833D-04,
     . -3.861D-05, -2.268D-05, +1.454D-06, +3.860D-07, -1.068D-07,
     . +0.000D+00, -2.658D-02, -1.947D-03, +7.131D-04, -3.506D-05,
     . +1.885D-07, +5.792D-07, +3.990D-08, +2.000D-08, -5.700D-09/
*     Define a parameter t
      T = DSIN(DLAT)
*     Define degree n and order m EGM
      N = 9
      M = 9
*     Determine n!  (factorial)  moved by 1
      DFAC(1) = 1
      DO I = 1,(2*N + 1)
        DFAC(I+1) = DFAC(I)*I
      END DO
*     Determine Legendre functions (Heiskanen and Moritz,
*     Physical Geodesy, 1967, eq. 1-62)
      DO I = 0,N
        DO J = 0,MIN(I,M)
          IR = INT((I - J)/2)
          SUM1 = 0
          DO K = 0,IR
            SUM1 = SUM1 + (-1)**K*DFAC(2*I - 2*K + 1)/DFAC(K + 1)/
     .       DFAC(I - K + 1)/DFAC(I - J - 2*K + 1)*T**(I - J - 2*K)
          END DO
*         Legendre functions moved by 1
          P(I + 1,J + 1) = 1.D0/2**I*DSQRT((1 - T**2)**(J))*SUM1
        END DO
      END DO
*     Calculate spherical harmonics
      I = 0
      DO N = 0,9
        DO M = 0,N
          I = I + 1
          AP(I) = P(N+1,M+1)*DCOS(M*DLON)
          BP(I) = P(N+1,M+1)*DSIN(M*DLON)
        END DO
      END DO
*     Compute hydrostatic mapping function
      BH = 0.0029D0
      C0H = 0.062D0
      IF (DLAT.LT.0D0) THEN ! SOUTHERN HEMISPHERE
        PHH  = PI
        C11H = 0.007D0
        C10H = 0.002D0
      ELSE                  ! NORTHERN HEMISPHERE
        PHH  = 0D0
        C11H = 0.005D0
        C10H = 0.001D0
      END IF
      CH = C0H + ((DCOS(DOY/365.25D0*TWOPI + PHH)+1D0)*C11H/2D0 + C10H)*
     .           (1D0-DCOS(DLAT))
      AHM = 0D0
      AHA = 0D0
      DO I = 1,55
        AHM = AHM + (AH_MEAN(I)*AP(I) + BH_MEAN(I)*BP(I))*1D-5
        AHA = AHA + (AH_AMP(I) *AP(I) + BH_AMP(I) *BP(I))*1D-5
      END DO
      AH  = AHM + AHA*DCOS(DOY/365.25D0*TWOPI)
      SINE   = DSIN(PI/2D0 - ZD)
      BETA   = BH/( SINE + CH  )
      GAMMA  = AH/( SINE + BETA)
      TOPCON = (1D0 + AH/(1D0 + BH/(1D0 + CH)))
      GMFH   = TOPCON/(SINE+GAMMA)
*     Height correction for hydrostatic mapping function from Niell (1996)
      A_HT = 2.53D-5
      B_HT = 5.49D-3
      C_HT = 1.14D-3
      HS_KM  = DHGT/1000D0
      BETA   = B_HT/( SINE + C_HT )
      GAMMA  = A_HT/( SINE + BETA)
      TOPCON = (1D0 + A_HT/(1D0 + B_HT/(1D0 + C_HT)))
      HT_CORR_COEF = 1D0/SINE - TOPCON/(SINE + GAMMA)
      HT_CORR      = HT_CORR_COEF * HS_KM
      GMFH         = GMFH + HT_CORR
*     Compute wet mapping function
      BW = 0.00146D0
      CW = 0.04391D0
      AWM = 0D0
      AWA = 0D0
      DO I = 1,55
        AWM = AWM + (AW_MEAN(I)*AP(I) + BW_MEAN(I)*BP(I))*1D-5
        AWA = AWA + (AW_AMP(I) *AP(I) + BW_AMP(I) *BP(I))*1D-5
      END DO
      AW =  AWM + AWA*DCOS(DOY/365.25D0*TWOPI)
      BETA   = BW/( SINE + CW )
      GAMMA  = AW/( SINE + BETA)
      TOPCON = (1D0 + AW/(1D0 + BW/(1D0 + CW)))
      GMFW   = TOPCON/(SINE+GAMMA)
     
* Finished.

      END     

阿英讲gfortran编译运行GMF(Global Mapping Functions)的Fortran 77程序 - 阿英 - Mr.Right
 
  评论这张
 
阅读(215)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017