Fortran 中逻辑变量的分段错误
Segmentation fault on logical variable in Fortran
我正在尝试 运行 在这里找到的子例程 'condensation' (http://mesonh.aero.obs-mip.fr/chaboureau/PUB/NCL/)，因此我写了一个 'main' 程序来初始化数组并调用子程序。
我在 'condensation' 子程序中做了 'print' 语句来查找错误，我发现问题（分段错误）发生在每次提到逻辑变量 'LUSERI'.
但是，我不知道为什么。
在主程序中，我写了：
program main
logical :: luseri
...
luseri = .true.
...
call condensation(...,luseri)
end program main
（'LUSERI' 是子例程中的最后一个参数）
似乎一切正常：主程序中的变量声明和赋值，子程序中的声明'condensation'及其提及。
这是我写的'main'程序(main_cst.f90):
program main_cst
implicit none
integer, parameter :: klon = 8 ! horizontal dimension
integer, parameter :: klev = 28 ! vertical dimension
integer, parameter :: kidia = 1 ! value of the first point in x; default=1
integer, parameter :: kfdia = 8 ! value of the last point in x; default=KLON
integer, parameter :: kbdia = 1 ! vert. comp. start at KBDIA that is at least 1
integer, parameter :: ktdia = 1 ! vert. comp. can be limited to KLEV + 1 - KTDIA
! default=1
logical :: luseri ! logical switch to compute both liquid
! and solid condensate (LUSERI=.TRUE.)
! or only liquid condensate (LUSERI=.FALSE.)
real, dimension(klon,klev) :: ppabs ! pressure (Pa)
real, dimension(klon,klev) :: pzz ! height of model levels (m)
real, dimension(klon,klev) :: pt ! grid scale T (K)
real, dimension(klon,klev) :: prv ! grid scale water vapor mixing ratio (kg/kg)
real, dimension(klon,klev) :: pmflx ! convective mass flux (kg/(s m^2))
real, dimension(klon,klev) :: prc ! grid scale r_c mixing ratio (kg/kg)
real, dimension(klon,klev) :: pri ! grid scale r_i (kg/kg)
integer :: kx ! horizontal loop counter
do kx = kidia, kfdia
ppabs(kx,:)=(/1000, 975, 950, 925, 900, 875, 850, 825, 800, 775, 750, 725, 700, &
675, 650, 600, 500, 300, 250, 200, 150, 100, 70, 50, 30, 20, 10, 3/)
pzz(kx,:)=(/111.31, 333.80, 561.25, 794.10, 1032.35, 1276.02, 1525.42, &
1780.62, 2041.78, 2309.21, 2583.33, 2864.61, 3153.70, 3451.53, &
3758.06, 4397.20, 5770.03, 9485.69, 10710.62, 12141.42, 13897.52, &
16283.89, 18361.25, 20367.34, 23523.77, 26158.47, 30748.48, 38877.28/)
pt(kx,:)=(/300.81, 299.97, 299.63, 298.31, 296.93, 295.81, 294.07, 292.04, &
290.02, 287.85, 285.50, 283.34, 281.39, 279.55, 278.16, 275.40, &
267.42, 239.07, 229.57, 219.98, 208.53, 197.10, 193.60, 200.84, &
213.77, 221.35, 230.21, 231.10/)
prv(kx,:)=(/0.012570000, 0.012460000, 0.011830000, 0.011390000, 0.010560000, &
0.009922000, 0.009529000, 0.009226000, 0.009008000, 0.008739000, &
0.008411000, 0.007836000, 0.007077000, 0.006153000, 0.004703000, &
0.002451000, 0.000774100, 0.000107400, 0.000056510, 0.000031610, &
0.000008088, 0.000004136, 0.000002686, 0.000002901, 0.000003876, &
0.000004562, 0.000004388, 0.000007886/)
end do
pmflx(:,:) = 0.0
prc(:,:) = 0.0
pri(:,:) = 0.0
luseri = .true.
call condensation(klon, klev, kidia, kfdia, kbdia, ktdia, ppabs, pzz, pt, prv, &
pmflx, prc, pri, luseri)
end program main_cst
冷凝子例程 (condensation.f90) 可在 http://mesonh.aero.obs-mip.fr/chaboureau/PUB/NCL/:
处找到
! ######spl
SUBROUTINE CONDENSATION( KLON, KLEV, KIDIA, KFDIA, KBDIA, KTDIA, &
PPABS, PZZ, PT, PRV, PRC, PRI, PMFLX, PCLDFR, LUSERI )
! ############################################################################
!
!!
!! PURPOSE
!! -------
!!** Routine to diagnose cloud fraction and liquid and ice condensate mixing ratios
!!
!!
!!** METHOD
!! ------
!! Based on the large-scale fields of temperature, water vapor, and possibly
!! liquid and solid condensate, the conserved quantities r_t and h_l are constructed
!! and then fractional cloudiness, liquid and solid condensate is diagnosed.
!!
!! The total variance is parameterized as the sum of stratiform/turbulent variance
!! and a convective variance.
!! The turbulent variance is parameterized as a function of first-order moments, and
!! the convective variance is modelled as a function of the convective mass flux (units kg/s m^2)
!! as provided by a mass flux convection scheme.
!!
!! Nota: if the host model does not use prognostic values for liquid and solid condensate
!! or does not provide a convective mass flux, put all these values to zero.
!! Also, it is supposed that vertical model levels are numbered from
!! 1 to KLEV, where 1 is the first model level above the surface
!!
!!
!!
!! EXTERNAL
!! --------
!! INI_CST
!!
!! IMPLICIT ARGUMENTS
!! ------------------
!! Module MODD_CST : contains physical constants
!!
!! REFERENCE
!! ---------
!! Chaboureau J.P. and P. Bechtold (J. Atmos. Sci. 2002)
!!
!! AUTHOR
!! ------
!! P. BECHTOLD * Laboratoire d'Aerologie *
!!
!! MODIFICATIONS
!! -------------
!! Original 13/06/2001
!! modified 20/03/2002 : add convective Sigma_s and improve turbulent
!! length-scale in boundary-layer and near tropopause
!!
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
USE MODD_CST
!
IMPLICIT NONE
!
!* 0.1 Declarations of dummy arguments :
!
!
INTEGER, INTENT(IN) :: KLON ! horizontal dimension
INTEGER, INTENT(IN) :: KLEV ! vertical dimension
INTEGER, INTENT(IN) :: KIDIA ! value of the first point in x
! default=1
INTEGER, INTENT(IN) :: KFDIA ! value of the last point in x
! default=KLON
INTEGER, INTENT(IN) :: KBDIA ! vertical computations start at
! ! KBDIA that is at least 1
INTEGER, INTENT(IN) :: KTDIA ! vertical computations can be
! limited to KLEV + 1 - KTDIA
! default=1
REAL, DIMENSION(KLON,KLEV), INTENT(IN) :: PPABS ! pressure (Pa)
REAL, DIMENSION(KLON,KLEV), INTENT(IN) :: PZZ ! height of model levels (m)
REAL, DIMENSION(KLON,KLEV), INTENT(IN) :: PT ! grid scale T (K)
REAL, DIMENSION(KLON,KLEV), INTENT(IN) :: PRV ! grid scale water vapor mixing ratio (kg/kg)
LOGICAL :: LUSERI ! logical switch to compute both
! liquid and solid condensate (LUSERI=.TRUE.)
! or only liquid condensate (LUSERI=.FALSE.)
REAL, DIMENSION(KLON,KLEV), INTENT(IN) :: PMFLX ! convective mass flux (kg/(s m^2))
REAL, DIMENSION(KLON,KLEV), INTENT(INOUT) :: PRC ! grid scale r_c mixing ratio (kg/kg)
REAL, DIMENSION(KLON,KLEV), INTENT(INOUT) :: PRI ! grid scale r_i (kg/kg)
REAL, DIMENSION(KLON,KLEV), INTENT(OUT) :: PCLDFR ! fractional cloudiness (between 0 and 1)
!
!
!* 0.2 Declarations of local variables :
!
INTEGER :: JI, JK, JKT, JKP, JKM ! loop index
REAL, DIMENSION(KLON,KLEV) :: ZTLK, ZRT ! work arrays for T_l, r_t
REAL, DIMENSION(KLON,KLEV) :: ZL ! length-scale
INTEGER, DIMENSION(KLON) :: ITPL ! top levels of tropopause/highest inversion
REAL, DIMENSION(KLON) :: ZTMIN ! min Temp. related to ITPL
!
REAL :: ZTEMP, ZLV, ZLS, ZTL, ZPV, ZQSL, ZPIV, ZQSI, ZFRAC, ZCOND, ZCPD ! thermodynamics
REAL :: ZLL, DZZ, ZZZ ! length scales
REAL :: ZAH, ZA, ZB, ZSBAR, ZQ1, ZSIGMA, ZDRW, ZDTL ! related to computation of Sig_s
REAL :: ZSIG_CONV ! convective part of Sig_s
!
!* 0.3 Definition of constants :
!
!-------------------------------------------------------------------------------
!
REAL :: ZL0 = 600. ! tropospheric length scale
! changed to 600 m instead of 900 m to give a consistent
! value (linear increase) in general 500 m deep oceanic
! mixed layer - but could be put back to 900 m if wished
REAL :: ZCSIGMA = 0.2 ! constant in sigma_s parameterization
REAL :: ZCSIG_CONV = 0.30E-2 ! scaling factor for ZSIG_CONV as function of mass flux
!
!-------------------------------------------------------------------------------
!
CALL INI_CST ! Initialize thermodynamic constants in module MODD_CST
PCLDFR(:,:) = 0. ! Initialize values
JKT = KLEV+1-KTDIA
DO JK=KBDIA,JKT
DO JI=KIDIA,KFDIA
ZTEMP = PT(JI,JK)
!latent heat of vaporisation/sublimation
ZLV = XLVTT + ( XCPV - XCL ) * ( ZTEMP - XTT )
ZLS = XLSTT + ( XCPV - XCI ) * ( ZTEMP - XTT )
!store temperature at saturation and total water mixing ratio
ZRT(JI,JK) = PRV(JI,JK) + PRC(JI,JK) + PRI(JI,JK)
ZCPD = XCPD + ZRT(JI,JK) * XCPV
ZTLK(JI,JK) = ZTEMP - ZLV*PRC(JI,JK)/ZCPD - ZLS*PRI(JI,JK)/ZCPD
END DO
END DO
!-------------------------------------------------------------------------------
! Determine tropopause/inversion height from minimum temperature
ITPL(:) = KBDIA+1
ZTMIN(:) = 400.
DO JK = KBDIA+1,JKT-1
DO JI=KIDIA,KFDIA
IF ( PT(JI,JK) < ZTMIN(JI) ) THEN
ZTMIN(JI) = PT(JI,JK)
ITPL(JI) = JK
END IF
END DO
END DO
! Set the mixing length scale - used for computing the "turbulent part" of Sigma_s
ZL(:,KBDIA) = 20.
DO JK = KBDIA+1,JKT
DO JI=KIDIA,KFDIA
! free troposphere
ZL(JI,JK) = ZL0
JKP = ITPL(JI)
ZZZ = PZZ(JI,JK) - PZZ(JI,KBDIA)
! approximate length for boundary-layer : linear increase
IF ( ZL0 > ZZZ ) ZL(JI,JK) = ZZZ
! gradual decrease of length-scale near and above tropopause/top inversion
IF ( ZZZ > 0.9*(PZZ(JI,JKP)-PZZ(JI,KBDIA)) ) &
ZL(JI,JK) = .6 * ZL(JI,JK-1)
END DO
END DO
!-------------------------------------------------------------------------------
DO JK=KBDIA+1,JKT-1
JKP=JK+1
JKM=JK-1
DO JI=KIDIA,KFDIA
ZTEMP = PT(JI,JK)
!latent heat of vaporisation/sublimation
ZLV = XLVTT + ( XCPV - XCL ) * ( ZTEMP - XTT )
ZLS = XLSTT + ( XCPV - XCI ) * ( ZTEMP - XTT )
ZCPD = XCPD + ZRT(JI,JK) * XCPV
!temperature at saturation
ZTL = ZTEMP - ZLV*PRC(JI,JK)/ZCPD - ZLS*PRI(JI,JK)/ZCPD
!saturated water vapor mixing ratio over liquid water
ZPV = EXP( XALPW - XBETAW / ZTL - XGAMW * LOG( ZTL ) )
ZQSL = XRD / XRV * ZPV / ( PPABS(JI,JK) - ZPV )
!saturated water vapor mixing ratio over ice
ZPIV = EXP( XALPI - XBETAI / ZTL - XGAMI * LOG( ZTL ) )
ZQSI = XRD / XRV * ZPIV / ( PPABS(JI,JK) - ZPIV )
!interpolate between liquid and solid as function of temperature
! glaciation interval is specified here to 20 K
ZFRAC = ( ZTL - 250.16 ) / ( XTT - 250.16 ) ! liquid/solid fraction
ZFRAC = MAX( 0., MIN(1., ZFRAC ) )
ZFRAC = ZFRAC * ZFRAC
IF(.NOT. LUSERI) ZFRAC=1.
ZQSL = ( 1. - ZFRAC ) * ZQSI + ZFRAC * ZQSL
ZLV = ( 1. - ZFRAC ) * ZLS + ZFRAC * ZLV
!coefficients a and b
ZAH = ZLV * ZQSL / ( XRV * ZTL**2 )
ZA = 1. / ( 1. + ZLV/ZCPD * ZAH )
ZB = ZAH * ZA
!parameterize Sigma_s with first_order closure
DZZ = PZZ(JI,JKP) - PZZ(JI,JKM)
ZDRW = ZRT(JI,JKP) - ZRT(JI,JKM)
ZDTL = ZTLK(JI,JKP) - ZTLK(JI,JKM) + XG/ZCPD * DZZ
ZLL = ZL(JI,JK)
ZSIG_CONV = ZCSIG_CONV * PMFLX(JI,JK) / ZA ! standard deviation due to convection
ZSIGMA = SQRT( MAX( 1.E-25, ZCSIGMA*ZCSIGMA* ZLL*ZLL/(DZZ*DZZ) * ( &
ZA*ZA*ZDRW*ZDRW - 2.*ZA*ZB*ZDRW*ZDTL + ZB*ZB*ZDTL*ZDTL ) &
+ ZSIG_CONV * ZSIG_CONV ) )
!zsigma should be of order 4.e-4 in lowest 5 km of atmosphere
ZSIGMA = MAX( ZSIGMA, 1.E-12 )
!normalized saturation deficit
ZSBAR = ZA * ( ZRT(JI,JK) - ZQSL )
ZQ1 = ZSBAR / ZSIGMA
!cloud fraction
PCLDFR(JI,JK) = MAX( 0., MIN(1.,0.5+0.36*ATAN(1.55*ZQ1)) )
!total condensate
IF (ZQ1 > 0. .AND. ZQ1 <= 2. ) THEN
ZCOND = EXP(-1.)+.66*ZQ1+.086*ZQ1*ZQ1
ELSE IF (ZQ1 > 2.) THEN
ZCOND = ZQ1
ELSE
ZCOND = EXP( 1.2*ZQ1-1 )
END IF
ZCOND = ZCOND * ZSIGMA
if ( zcond<1.e-6) then
zcond = 0.
pcldfr(ji,jk) = 0.
end if
PRC(JI,JK) = ZFRAC * ZCOND ! liquid condensate
IF (LUSERI) THEN
PRI(JI,JK) = (1.-ZFRAC) * ZCOND ! solid condensate
END IF
! compute s'rl'/Sigs^2
! used in w'rl'= w's' * s'rl'/Sigs^2
! PSIGRC(JI,JK) = PCLDFR(JI,JK) ! Gaussian relation
END DO
END DO
!
END SUBROUTINE CONDENSATION
还有两个文件，也在 http://mesonh.aero.obs-mip.fr/chaboureau/PUB/NCL/ 中找到。
ini_cst.f90:
! ######spl
SUBROUTINE INI_CST
! ##################
!
!!**** *INI_CST * - routine to initialize the constants modules
!!
!! PURPOSE
!! -------
! The purpose of this routine is to initialize the constants
! stored in modules MODD_CST
!
!
!!** METHOD
!! ------
!! The thermodynamic constants are set to their numerical values
!!
!!
!! EXTERNAL
!! --------
!!
!! IMPLICIT ARGUMENTS
!! ------------------
!! Module MODD_CST : contains physical constants
!!
!! REFERENCE
!! ---------
!! Chaboureau J.P. and P. Bechtold (J. Atmos. Sci. 2002)
!!
!!
!! AUTHOR
!! ------
!! P. BECHTOLD * Laboratoire d'Aerologie *
!!
!! MODIFICATIONS
!! -------------
!! Original 13/06/2001
!! modified 20/03/2002 : add convective Sigma_s and improve turbulent
!! length-scale in boundary-layer and near tropopause
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
USE MODD_CST
!
IMPLICIT NONE
!
!-------------------------------------------------------------------------------
!
!* 1. Set the fundamental thermodynamical constants
! these have the same values (not names) as in ARPEGE IFS
! -------------------------------------------------------
!
!
XP00 = 1.E5 ! reference pressure
XPI = 3.141592654 ! Pi
XG = 9.80665 ! gravity constant
XMD = 28.9644E-3 ! molecular weight of dry air
XMV = 18.0153E-3 ! molecular weight of water vapor
XRD = 287.05967 ! gaz constant for dry air
XRV = 461.524993 ! gaz constant for water vapor
XCPD = 1004.708845 ! specific heat of dry air
XCPV = 1846.1 ! specific heat of water vapor
XRHOLW = 1000. ! density of liquid water
XCL = 4218. ! specific heat of liquid water
XCI = 2106. ! specific heat of ice
XTT = 273.16 ! triple point temperature
XLVTT = 2.5008E6 ! latent heat of vaporisation at XTT
XLSTT = 2.8345E6 ! latent heat of sublimation at XTT
XLMTT = 0.3337E6 ! latent heat of melting at XTT
XESTT = 611.14 ! saturation pressure at XTT
XALPW = 60.22416 ! constants in saturation pressure over liquid water
XBETAW = 6822.459384
XGAMW = 5.13948
XALPI = 32.62116 ! constants in saturation pressure over ice
XBETAI = 6295.421
XGAMI = 0.56313
!
!
END SUBROUTINE INI_CST
! ######spl
MODULE MODD_CST
! ###############
!
IMPLICIT NONE
!
REAL, SAVE :: XP00 ! reference pressure
REAL, SAVE :: XPI ! Pi
REAL, SAVE :: XG ! gravity constant
REAL, SAVE :: XMD ! molecular weight of dry air
REAL, SAVE :: XMV ! molecular weight of water vapor
REAL, SAVE :: XRD ! gaz constant for dry air
REAL, SAVE :: XRV ! gaz constant for water vapor
REAL, SAVE :: XCPD ! specific heat of dry air
REAL, SAVE :: XCPV ! specific heat of water vapor
REAL, SAVE :: XRHOLW ! density of liquid water
REAL, SAVE :: XCL ! specific heat of liquid water
REAL, SAVE :: XCI ! specific heat of ice
REAL, SAVE :: XTT ! triple point temperature
REAL, SAVE :: XLVTT ! latent heat of vaporisation at XTT
REAL, SAVE :: XLSTT ! latent heat of sublimation at XTT
REAL, SAVE :: XLMTT ! latent heat of melting at XTT
REAL, SAVE :: XESTT ! saturation pressure at XTT
REAL, SAVE :: XALPW ! constants in saturation pressure over liquid water
REAL, SAVE :: XBETAW
REAL, SAVE :: XGAMW
REAL, SAVE :: XALPI ! constants in saturation pressure over ice
REAL, SAVE :: XBETAI
REAL, SAVE :: XGAMI
!
END MODULE MODD_CST
modd_cst.f90:
! ######spl
MODULE MODD_CST
! ###############
!
IMPLICIT NONE
!
REAL, SAVE :: XP00 ! reference pressure
REAL, SAVE :: XPI ! Pi
REAL, SAVE :: XG ! gravity constant
REAL, SAVE :: XMD ! molecular weight of dry air
REAL, SAVE :: XMV ! molecular weight of water vapor
REAL, SAVE :: XRD ! gaz constant for dry air
REAL, SAVE :: XRV ! gaz constant for water vapor
REAL, SAVE :: XCPD ! specific heat of dry air
REAL, SAVE :: XCPV ! specific heat of water vapor
REAL, SAVE :: XRHOLW ! density of liquid water
REAL, SAVE :: XCL ! specific heat of liquid water
REAL, SAVE :: XCI ! specific heat of ice
REAL, SAVE :: XTT ! triple point temperature
REAL, SAVE :: XLVTT ! latent heat of vaporisation at XTT
REAL, SAVE :: XLSTT ! latent heat of sublimation at XTT
REAL, SAVE :: XLMTT ! latent heat of melting at XTT
REAL, SAVE :: XESTT ! saturation pressure at XTT
REAL, SAVE :: XALPW ! constants in saturation pressure over liquid water
REAL, SAVE :: XBETAW
REAL, SAVE :: XGAMW
REAL, SAVE :: XALPI ! constants in saturation pressure over ice
REAL, SAVE :: XBETAI
REAL, SAVE :: XGAMI
!
END MODULE MODD_CST
这是我正在使用的所有代码。
我使用命令编译：
g95 modd_cst.f90 ini_cst.f90 condensation.f90 main_cst.f90 -o cst.g95
您在调用中输入的内容与子例程获得的内容之间似乎不匹配
call condensation(klon, klev, kidia, kfdia, kbdia, ktdia, ppabs, pzz, pt, prv,
pmflx, prc, pri, luseri)
子程序需要
SUBROUTINE CONDENSATION( KLON, KLEV, KIDIA, KFDIA, KBDIA, KTDIA, &
PPABS, PZZ, PT, PRV, PRC, PRI, PMFLX, PCLDFR, LUSERI )
所以在调用它的 prv pmflx prc pri
在子程序中它的 prv prc pri pmflx。此外，您还缺少一个传递给子例程的变量。 Fortran 对您传递的变量的顺序很敏感。在不检查其余代码的情况下，我认为应该可以解决问题。
我正在尝试 运行 在这里找到的子例程 'condensation' (http://mesonh.aero.obs-mip.fr/chaboureau/PUB/NCL/)，因此我写了一个 'main' 程序来初始化数组并调用子程序。
我在 'condensation' 子程序中做了 'print' 语句来查找错误，我发现问题（分段错误）发生在每次提到逻辑变量 'LUSERI'.
但是，我不知道为什么。
在主程序中，我写了：
program main
logical :: luseri
...
luseri = .true.
...
call condensation(...,luseri)
end program main
（'LUSERI' 是子例程中的最后一个参数）
似乎一切正常：主程序中的变量声明和赋值，子程序中的声明'condensation'及其提及。
这是我写的'main'程序(main_cst.f90):
program main_cst
implicit none
integer, parameter :: klon = 8 ! horizontal dimension
integer, parameter :: klev = 28 ! vertical dimension
integer, parameter :: kidia = 1 ! value of the first point in x; default=1
integer, parameter :: kfdia = 8 ! value of the last point in x; default=KLON
integer, parameter :: kbdia = 1 ! vert. comp. start at KBDIA that is at least 1
integer, parameter :: ktdia = 1 ! vert. comp. can be limited to KLEV + 1 - KTDIA
! default=1
logical :: luseri ! logical switch to compute both liquid
! and solid condensate (LUSERI=.TRUE.)
! or only liquid condensate (LUSERI=.FALSE.)
real, dimension(klon,klev) :: ppabs ! pressure (Pa)
real, dimension(klon,klev) :: pzz ! height of model levels (m)
real, dimension(klon,klev) :: pt ! grid scale T (K)
real, dimension(klon,klev) :: prv ! grid scale water vapor mixing ratio (kg/kg)
real, dimension(klon,klev) :: pmflx ! convective mass flux (kg/(s m^2))
real, dimension(klon,klev) :: prc ! grid scale r_c mixing ratio (kg/kg)
real, dimension(klon,klev) :: pri ! grid scale r_i (kg/kg)
integer :: kx ! horizontal loop counter
do kx = kidia, kfdia
ppabs(kx,:)=(/1000, 975, 950, 925, 900, 875, 850, 825, 800, 775, 750, 725, 700, &
675, 650, 600, 500, 300, 250, 200, 150, 100, 70, 50, 30, 20, 10, 3/)
pzz(kx,:)=(/111.31, 333.80, 561.25, 794.10, 1032.35, 1276.02, 1525.42, &
1780.62, 2041.78, 2309.21, 2583.33, 2864.61, 3153.70, 3451.53, &
3758.06, 4397.20, 5770.03, 9485.69, 10710.62, 12141.42, 13897.52, &
16283.89, 18361.25, 20367.34, 23523.77, 26158.47, 30748.48, 38877.28/)
pt(kx,:)=(/300.81, 299.97, 299.63, 298.31, 296.93, 295.81, 294.07, 292.04, &
290.02, 287.85, 285.50, 283.34, 281.39, 279.55, 278.16, 275.40, &
267.42, 239.07, 229.57, 219.98, 208.53, 197.10, 193.60, 200.84, &
213.77, 221.35, 230.21, 231.10/)
prv(kx,:)=(/0.012570000, 0.012460000, 0.011830000, 0.011390000, 0.010560000, &
0.009922000, 0.009529000, 0.009226000, 0.009008000, 0.008739000, &
0.008411000, 0.007836000, 0.007077000, 0.006153000, 0.004703000, &
0.002451000, 0.000774100, 0.000107400, 0.000056510, 0.000031610, &
0.000008088, 0.000004136, 0.000002686, 0.000002901, 0.000003876, &
0.000004562, 0.000004388, 0.000007886/)
end do
pmflx(:,:) = 0.0
prc(:,:) = 0.0
pri(:,:) = 0.0
luseri = .true.
call condensation(klon, klev, kidia, kfdia, kbdia, ktdia, ppabs, pzz, pt, prv, &
pmflx, prc, pri, luseri)
end program main_cst
冷凝子例程 (condensation.f90) 可在 http://mesonh.aero.obs-mip.fr/chaboureau/PUB/NCL/:
处找到! ######spl
SUBROUTINE CONDENSATION( KLON, KLEV, KIDIA, KFDIA, KBDIA, KTDIA, &
PPABS, PZZ, PT, PRV, PRC, PRI, PMFLX, PCLDFR, LUSERI )
! ############################################################################
!
!!
!! PURPOSE
!! -------
!!** Routine to diagnose cloud fraction and liquid and ice condensate mixing ratios
!!
!!
!!** METHOD
!! ------
!! Based on the large-scale fields of temperature, water vapor, and possibly
!! liquid and solid condensate, the conserved quantities r_t and h_l are constructed
!! and then fractional cloudiness, liquid and solid condensate is diagnosed.
!!
!! The total variance is parameterized as the sum of stratiform/turbulent variance
!! and a convective variance.
!! The turbulent variance is parameterized as a function of first-order moments, and
!! the convective variance is modelled as a function of the convective mass flux (units kg/s m^2)
!! as provided by a mass flux convection scheme.
!!
!! Nota: if the host model does not use prognostic values for liquid and solid condensate
!! or does not provide a convective mass flux, put all these values to zero.
!! Also, it is supposed that vertical model levels are numbered from
!! 1 to KLEV, where 1 is the first model level above the surface
!!
!!
!!
!! EXTERNAL
!! --------
!! INI_CST
!!
!! IMPLICIT ARGUMENTS
!! ------------------
!! Module MODD_CST : contains physical constants
!!
!! REFERENCE
!! ---------
!! Chaboureau J.P. and P. Bechtold (J. Atmos. Sci. 2002)
!!
!! AUTHOR
!! ------
!! P. BECHTOLD * Laboratoire d'Aerologie *
!!
!! MODIFICATIONS
!! -------------
!! Original 13/06/2001
!! modified 20/03/2002 : add convective Sigma_s and improve turbulent
!! length-scale in boundary-layer and near tropopause
!!
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
USE MODD_CST
!
IMPLICIT NONE
!
!* 0.1 Declarations of dummy arguments :
!
!
INTEGER, INTENT(IN) :: KLON ! horizontal dimension
INTEGER, INTENT(IN) :: KLEV ! vertical dimension
INTEGER, INTENT(IN) :: KIDIA ! value of the first point in x
! default=1
INTEGER, INTENT(IN) :: KFDIA ! value of the last point in x
! default=KLON
INTEGER, INTENT(IN) :: KBDIA ! vertical computations start at
! ! KBDIA that is at least 1
INTEGER, INTENT(IN) :: KTDIA ! vertical computations can be
! limited to KLEV + 1 - KTDIA
! default=1
REAL, DIMENSION(KLON,KLEV), INTENT(IN) :: PPABS ! pressure (Pa)
REAL, DIMENSION(KLON,KLEV), INTENT(IN) :: PZZ ! height of model levels (m)
REAL, DIMENSION(KLON,KLEV), INTENT(IN) :: PT ! grid scale T (K)
REAL, DIMENSION(KLON,KLEV), INTENT(IN) :: PRV ! grid scale water vapor mixing ratio (kg/kg)
LOGICAL :: LUSERI ! logical switch to compute both
! liquid and solid condensate (LUSERI=.TRUE.)
! or only liquid condensate (LUSERI=.FALSE.)
REAL, DIMENSION(KLON,KLEV), INTENT(IN) :: PMFLX ! convective mass flux (kg/(s m^2))
REAL, DIMENSION(KLON,KLEV), INTENT(INOUT) :: PRC ! grid scale r_c mixing ratio (kg/kg)
REAL, DIMENSION(KLON,KLEV), INTENT(INOUT) :: PRI ! grid scale r_i (kg/kg)
REAL, DIMENSION(KLON,KLEV), INTENT(OUT) :: PCLDFR ! fractional cloudiness (between 0 and 1)
!
!
!* 0.2 Declarations of local variables :
!
INTEGER :: JI, JK, JKT, JKP, JKM ! loop index
REAL, DIMENSION(KLON,KLEV) :: ZTLK, ZRT ! work arrays for T_l, r_t
REAL, DIMENSION(KLON,KLEV) :: ZL ! length-scale
INTEGER, DIMENSION(KLON) :: ITPL ! top levels of tropopause/highest inversion
REAL, DIMENSION(KLON) :: ZTMIN ! min Temp. related to ITPL
!
REAL :: ZTEMP, ZLV, ZLS, ZTL, ZPV, ZQSL, ZPIV, ZQSI, ZFRAC, ZCOND, ZCPD ! thermodynamics
REAL :: ZLL, DZZ, ZZZ ! length scales
REAL :: ZAH, ZA, ZB, ZSBAR, ZQ1, ZSIGMA, ZDRW, ZDTL ! related to computation of Sig_s
REAL :: ZSIG_CONV ! convective part of Sig_s
!
!* 0.3 Definition of constants :
!
!-------------------------------------------------------------------------------
!
REAL :: ZL0 = 600. ! tropospheric length scale
! changed to 600 m instead of 900 m to give a consistent
! value (linear increase) in general 500 m deep oceanic
! mixed layer - but could be put back to 900 m if wished
REAL :: ZCSIGMA = 0.2 ! constant in sigma_s parameterization
REAL :: ZCSIG_CONV = 0.30E-2 ! scaling factor for ZSIG_CONV as function of mass flux
!
!-------------------------------------------------------------------------------
!
CALL INI_CST ! Initialize thermodynamic constants in module MODD_CST
PCLDFR(:,:) = 0. ! Initialize values
JKT = KLEV+1-KTDIA
DO JK=KBDIA,JKT
DO JI=KIDIA,KFDIA
ZTEMP = PT(JI,JK)
!latent heat of vaporisation/sublimation
ZLV = XLVTT + ( XCPV - XCL ) * ( ZTEMP - XTT )
ZLS = XLSTT + ( XCPV - XCI ) * ( ZTEMP - XTT )
!store temperature at saturation and total water mixing ratio
ZRT(JI,JK) = PRV(JI,JK) + PRC(JI,JK) + PRI(JI,JK)
ZCPD = XCPD + ZRT(JI,JK) * XCPV
ZTLK(JI,JK) = ZTEMP - ZLV*PRC(JI,JK)/ZCPD - ZLS*PRI(JI,JK)/ZCPD
END DO
END DO
!-------------------------------------------------------------------------------
! Determine tropopause/inversion height from minimum temperature
ITPL(:) = KBDIA+1
ZTMIN(:) = 400.
DO JK = KBDIA+1,JKT-1
DO JI=KIDIA,KFDIA
IF ( PT(JI,JK) < ZTMIN(JI) ) THEN
ZTMIN(JI) = PT(JI,JK)
ITPL(JI) = JK
END IF
END DO
END DO
! Set the mixing length scale - used for computing the "turbulent part" of Sigma_s
ZL(:,KBDIA) = 20.
DO JK = KBDIA+1,JKT
DO JI=KIDIA,KFDIA
! free troposphere
ZL(JI,JK) = ZL0
JKP = ITPL(JI)
ZZZ = PZZ(JI,JK) - PZZ(JI,KBDIA)
! approximate length for boundary-layer : linear increase
IF ( ZL0 > ZZZ ) ZL(JI,JK) = ZZZ
! gradual decrease of length-scale near and above tropopause/top inversion
IF ( ZZZ > 0.9*(PZZ(JI,JKP)-PZZ(JI,KBDIA)) ) &
ZL(JI,JK) = .6 * ZL(JI,JK-1)
END DO
END DO
!-------------------------------------------------------------------------------
DO JK=KBDIA+1,JKT-1
JKP=JK+1
JKM=JK-1
DO JI=KIDIA,KFDIA
ZTEMP = PT(JI,JK)
!latent heat of vaporisation/sublimation
ZLV = XLVTT + ( XCPV - XCL ) * ( ZTEMP - XTT )
ZLS = XLSTT + ( XCPV - XCI ) * ( ZTEMP - XTT )
ZCPD = XCPD + ZRT(JI,JK) * XCPV
!temperature at saturation
ZTL = ZTEMP - ZLV*PRC(JI,JK)/ZCPD - ZLS*PRI(JI,JK)/ZCPD
!saturated water vapor mixing ratio over liquid water
ZPV = EXP( XALPW - XBETAW / ZTL - XGAMW * LOG( ZTL ) )
ZQSL = XRD / XRV * ZPV / ( PPABS(JI,JK) - ZPV )
!saturated water vapor mixing ratio over ice
ZPIV = EXP( XALPI - XBETAI / ZTL - XGAMI * LOG( ZTL ) )
ZQSI = XRD / XRV * ZPIV / ( PPABS(JI,JK) - ZPIV )
!interpolate between liquid and solid as function of temperature
! glaciation interval is specified here to 20 K
ZFRAC = ( ZTL - 250.16 ) / ( XTT - 250.16 ) ! liquid/solid fraction
ZFRAC = MAX( 0., MIN(1., ZFRAC ) )
ZFRAC = ZFRAC * ZFRAC
IF(.NOT. LUSERI) ZFRAC=1.
ZQSL = ( 1. - ZFRAC ) * ZQSI + ZFRAC * ZQSL
ZLV = ( 1. - ZFRAC ) * ZLS + ZFRAC * ZLV
!coefficients a and b
ZAH = ZLV * ZQSL / ( XRV * ZTL**2 )
ZA = 1. / ( 1. + ZLV/ZCPD * ZAH )
ZB = ZAH * ZA
!parameterize Sigma_s with first_order closure
DZZ = PZZ(JI,JKP) - PZZ(JI,JKM)
ZDRW = ZRT(JI,JKP) - ZRT(JI,JKM)
ZDTL = ZTLK(JI,JKP) - ZTLK(JI,JKM) + XG/ZCPD * DZZ
ZLL = ZL(JI,JK)
ZSIG_CONV = ZCSIG_CONV * PMFLX(JI,JK) / ZA ! standard deviation due to convection
ZSIGMA = SQRT( MAX( 1.E-25, ZCSIGMA*ZCSIGMA* ZLL*ZLL/(DZZ*DZZ) * ( &
ZA*ZA*ZDRW*ZDRW - 2.*ZA*ZB*ZDRW*ZDTL + ZB*ZB*ZDTL*ZDTL ) &
+ ZSIG_CONV * ZSIG_CONV ) )
!zsigma should be of order 4.e-4 in lowest 5 km of atmosphere
ZSIGMA = MAX( ZSIGMA, 1.E-12 )
!normalized saturation deficit
ZSBAR = ZA * ( ZRT(JI,JK) - ZQSL )
ZQ1 = ZSBAR / ZSIGMA
!cloud fraction
PCLDFR(JI,JK) = MAX( 0., MIN(1.,0.5+0.36*ATAN(1.55*ZQ1)) )
!total condensate
IF (ZQ1 > 0. .AND. ZQ1 <= 2. ) THEN
ZCOND = EXP(-1.)+.66*ZQ1+.086*ZQ1*ZQ1
ELSE IF (ZQ1 > 2.) THEN
ZCOND = ZQ1
ELSE
ZCOND = EXP( 1.2*ZQ1-1 )
END IF
ZCOND = ZCOND * ZSIGMA
if ( zcond<1.e-6) then
zcond = 0.
pcldfr(ji,jk) = 0.
end if
PRC(JI,JK) = ZFRAC * ZCOND ! liquid condensate
IF (LUSERI) THEN
PRI(JI,JK) = (1.-ZFRAC) * ZCOND ! solid condensate
END IF
! compute s'rl'/Sigs^2
! used in w'rl'= w's' * s'rl'/Sigs^2
! PSIGRC(JI,JK) = PCLDFR(JI,JK) ! Gaussian relation
END DO
END DO
!
END SUBROUTINE CONDENSATION
还有两个文件，也在 http://mesonh.aero.obs-mip.fr/chaboureau/PUB/NCL/ 中找到。
ini_cst.f90:
! ######spl
SUBROUTINE INI_CST
! ##################
!
!!**** *INI_CST * - routine to initialize the constants modules
!!
!! PURPOSE
!! -------
! The purpose of this routine is to initialize the constants
! stored in modules MODD_CST
!
!
!!** METHOD
!! ------
!! The thermodynamic constants are set to their numerical values
!!
!!
!! EXTERNAL
!! --------
!!
!! IMPLICIT ARGUMENTS
!! ------------------
!! Module MODD_CST : contains physical constants
!!
!! REFERENCE
!! ---------
!! Chaboureau J.P. and P. Bechtold (J. Atmos. Sci. 2002)
!!
!!
!! AUTHOR
!! ------
!! P. BECHTOLD * Laboratoire d'Aerologie *
!!
!! MODIFICATIONS
!! -------------
!! Original 13/06/2001
!! modified 20/03/2002 : add convective Sigma_s and improve turbulent
!! length-scale in boundary-layer and near tropopause
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
USE MODD_CST
!
IMPLICIT NONE
!
!-------------------------------------------------------------------------------
!
!* 1. Set the fundamental thermodynamical constants
! these have the same values (not names) as in ARPEGE IFS
! -------------------------------------------------------
!
!
XP00 = 1.E5 ! reference pressure
XPI = 3.141592654 ! Pi
XG = 9.80665 ! gravity constant
XMD = 28.9644E-3 ! molecular weight of dry air
XMV = 18.0153E-3 ! molecular weight of water vapor
XRD = 287.05967 ! gaz constant for dry air
XRV = 461.524993 ! gaz constant for water vapor
XCPD = 1004.708845 ! specific heat of dry air
XCPV = 1846.1 ! specific heat of water vapor
XRHOLW = 1000. ! density of liquid water
XCL = 4218. ! specific heat of liquid water
XCI = 2106. ! specific heat of ice
XTT = 273.16 ! triple point temperature
XLVTT = 2.5008E6 ! latent heat of vaporisation at XTT
XLSTT = 2.8345E6 ! latent heat of sublimation at XTT
XLMTT = 0.3337E6 ! latent heat of melting at XTT
XESTT = 611.14 ! saturation pressure at XTT
XALPW = 60.22416 ! constants in saturation pressure over liquid water
XBETAW = 6822.459384
XGAMW = 5.13948
XALPI = 32.62116 ! constants in saturation pressure over ice
XBETAI = 6295.421
XGAMI = 0.56313
!
!
END SUBROUTINE INI_CST
! ######spl
MODULE MODD_CST
! ###############
!
IMPLICIT NONE
!
REAL, SAVE :: XP00 ! reference pressure
REAL, SAVE :: XPI ! Pi
REAL, SAVE :: XG ! gravity constant
REAL, SAVE :: XMD ! molecular weight of dry air
REAL, SAVE :: XMV ! molecular weight of water vapor
REAL, SAVE :: XRD ! gaz constant for dry air
REAL, SAVE :: XRV ! gaz constant for water vapor
REAL, SAVE :: XCPD ! specific heat of dry air
REAL, SAVE :: XCPV ! specific heat of water vapor
REAL, SAVE :: XRHOLW ! density of liquid water
REAL, SAVE :: XCL ! specific heat of liquid water
REAL, SAVE :: XCI ! specific heat of ice
REAL, SAVE :: XTT ! triple point temperature
REAL, SAVE :: XLVTT ! latent heat of vaporisation at XTT
REAL, SAVE :: XLSTT ! latent heat of sublimation at XTT
REAL, SAVE :: XLMTT ! latent heat of melting at XTT
REAL, SAVE :: XESTT ! saturation pressure at XTT
REAL, SAVE :: XALPW ! constants in saturation pressure over liquid water
REAL, SAVE :: XBETAW
REAL, SAVE :: XGAMW
REAL, SAVE :: XALPI ! constants in saturation pressure over ice
REAL, SAVE :: XBETAI
REAL, SAVE :: XGAMI
!
END MODULE MODD_CST
modd_cst.f90:
! ######spl
MODULE MODD_CST
! ###############
!
IMPLICIT NONE
!
REAL, SAVE :: XP00 ! reference pressure
REAL, SAVE :: XPI ! Pi
REAL, SAVE :: XG ! gravity constant
REAL, SAVE :: XMD ! molecular weight of dry air
REAL, SAVE :: XMV ! molecular weight of water vapor
REAL, SAVE :: XRD ! gaz constant for dry air
REAL, SAVE :: XRV ! gaz constant for water vapor
REAL, SAVE :: XCPD ! specific heat of dry air
REAL, SAVE :: XCPV ! specific heat of water vapor
REAL, SAVE :: XRHOLW ! density of liquid water
REAL, SAVE :: XCL ! specific heat of liquid water
REAL, SAVE :: XCI ! specific heat of ice
REAL, SAVE :: XTT ! triple point temperature
REAL, SAVE :: XLVTT ! latent heat of vaporisation at XTT
REAL, SAVE :: XLSTT ! latent heat of sublimation at XTT
REAL, SAVE :: XLMTT ! latent heat of melting at XTT
REAL, SAVE :: XESTT ! saturation pressure at XTT
REAL, SAVE :: XALPW ! constants in saturation pressure over liquid water
REAL, SAVE :: XBETAW
REAL, SAVE :: XGAMW
REAL, SAVE :: XALPI ! constants in saturation pressure over ice
REAL, SAVE :: XBETAI
REAL, SAVE :: XGAMI
!
END MODULE MODD_CST
这是我正在使用的所有代码。
我使用命令编译：
g95 modd_cst.f90 ini_cst.f90 condensation.f90 main_cst.f90 -o cst.g95
您在调用中输入的内容与子例程获得的内容之间似乎不匹配
call condensation(klon, klev, kidia, kfdia, kbdia, ktdia, ppabs, pzz, pt, prv,
pmflx, prc, pri, luseri)
子程序需要
SUBROUTINE CONDENSATION( KLON, KLEV, KIDIA, KFDIA, KBDIA, KTDIA, &
PPABS, PZZ, PT, PRV, PRC, PRI, PMFLX, PCLDFR, LUSERI )
所以在调用它的 prv pmflx prc pri 在子程序中它的 prv prc pri pmflx。此外，您还缺少一个传递给子例程的变量。 Fortran 对您传递的变量的顺序很敏感。在不检查其余代码的情况下，我认为应该可以解决问题。