Skip to content

NSSL Microphysics log(0) FPE #1216

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
davegill merged 1 commit into from
Jul 14, 2020
Merged

NSSL Microphysics log(0) FPE #1216

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
89 changes: 62 additions & 27 deletions phys/module_mp_nssl_2mom.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -142,7 +142,7 @@ MODULE module_mp_nssl_2mom
public nssl_2mom_driver
public nssl_2mom_init
private gamma_sp,gamxinf,GAML02, GAML02d300, GAML02d500, fqvs, fqis
private gamma_dp, gamxinfdp
private gamma_dp, gamxinfdp, gamma_dpr
private delbk, delabk
private gammadp

Expand Down Expand Up @@ -675,7 +675,7 @@ MODULE module_mp_nssl_2mom
real :: ciacrratio(0:nqiacrratio,0:nqiacralpha)
real :: qiacrratio(0:nqiacrratio,0:nqiacralpha)
real :: ziacrratio(0:nqiacrratio,0:nqiacralpha)
double precision :: gamxinflu(0:nqiacrratio,0:nqiacralpha,10,2) ! last index for graupel (1) or hail (2)
double precision :: gamxinflu(0:nqiacrratio,0:nqiacralpha,12,2) ! last index for graupel (1) or hail (2)

integer, parameter :: ngdnmm = 9
real :: mmgraupvt(ngdnmm,3) ! Milbrandt and Morrison (2013) fall speed coefficients for graupel/hail
Expand Down Expand Up @@ -885,7 +885,8 @@ SUBROUTINE nssl_2mom_init( &
integer :: isub
real :: bxh,bxhl

real :: alp,ratio,x,y,y7
real :: alp,ratio !,x,y,y7
Copy link
Contributor

@davegill davegill Jun 29, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@MicroTed
Ted,
Instead of commenting these out, maybe just remove them?

double precision :: x,y,y2,y7
logical :: turn_on_ccna


Expand Down Expand Up @@ -1012,45 +1013,48 @@ SUBROUTINE nssl_2mom_init( &

DO j = 0,nqiacralpha
alp = float(j)*dqiacralpha
y = gamma_sp(1.+alp)
y = gamma_dpr(1.+alp)
y2 = gamma_dpr(real(2.+alp))
Copy link
Contributor

@davegill davegill Jun 29, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@MicroTed
Ted,
Why are these two lines different? Why does one need real and the other does not?

      y = gamma_dpr(1.+alp)
      y2 = gamma_dpr(real(2.+alp))

DO i = 0,nqiacrratio
ratio = float(i)*dqiacrratio
x = gamxinf( 1.+alp, ratio )
x = gamxinfdp( 1.+alp, ratio )
! write(0,*) 'i, x/y = ',i, x/y
ciacrratio(i,j) = x/y

! graupel (.,.,.,1)
gamxinflu(i,j,1,1) = x/y
gamxinflu(i,j,2,1) = gamxinf( 2.0+alp, ratio )/y
gamxinflu(i,j,3,1) = gamxinf( 2.5+alp+0.5*bxh, ratio )/y
gamxinflu(i,j,5,1) = (gamma_sp(5.0+alp) - gamxinf( 5.0+alp, ratio ))/y
gamxinflu(i,j,6,1) = (gamma_sp(5.5+alp+0.5*bxh) - gamxinf( 5.5+alp+0.5*bxh, ratio ))/y
gamxinflu(i,j,9,1) = gamxinf( 1.0+alp, ratio )/y
gamxinflu(i,j,10,1)= gamxinf( 4.0+alp, ratio )/y
gamxinflu(i,j,2,1) = gamxinfdp( 2.0+alp, ratio )/y
gamxinflu(i,j,3,1) = gamxinfdp( 2.5+alp+0.5*bxh, ratio )/y
gamxinflu(i,j,5,1) = (gamma_dpr(5.0+alp) - gamxinfdp( 5.0+alp, ratio ))/y
gamxinflu(i,j,6,1) = (gamma_dpr(5.5+alp+0.5*bxh) - gamxinfdp( 5.5+alp+0.5*bxh, ratio ))/y
gamxinflu(i,j,9,1) = gamxinfdp( 1.0+alp, ratio )/y
gamxinflu(i,j,10,1)= gamxinfdp( 4.0+alp, ratio )/y

gamxinflu(i,j,12,1) = gamxinfdp( 2.0+alp, ratio )/y2

! hail (.,.,.,2)
gamxinflu(i,j,1,2) = gamxinflu(i,j,1,1)
gamxinflu(i,j,2,2) = gamxinflu(i,j,2,1)
gamxinflu(i,j,3,2) = gamxinf( 2.5+alp+0.5*bxhl, ratio )/y
gamxinflu(i,j,3,2) = gamxinfdp( 2.5+alp+0.5*bxhl, ratio )/y
gamxinflu(i,j,5,2) = gamxinflu(i,j,5,1)
gamxinflu(i,j,6,2) = (gamma_sp(5.5+alp+0.5*bxhl) - gamxinf( 5.5+alp+0.5*bxhl, ratio ))/y
gamxinflu(i,j,6,2) = (gamma_dpr(5.5+alp+0.5*bxhl) - gamxinfdp( 5.5+alp+0.5*bxhl, ratio ))/y
gamxinflu(i,j,9,2) = gamxinflu(i,j,9,1)
gamxinflu(i,j,10,2)= gamxinflu(i,j,10,1)

IF ( alp > 1.1 ) THEN
! gamxinflu(i,j,7,1) = gamxinf( alp - 1., ratio )/y
gamxinflu(i,j,7,1) = (gamma_sp(alp - 1.) - gamxinf( alp - 1., ratio ))/y
! gamxinflu(i,j,8,1) = gamxinf( alp - 0.5 + 0.5*bxh, ratio )/y
gamxinflu(i,j,8,1) = (gamma_sp(alp - 0.5 + 0.5*bxh) - gamxinf( alp - 0.5 + 0.5*bxh, ratio ))/y
! gamxinflu(i,j,8,2) = gamxinf( alp - 0.5 + 0.5*bxhl, ratio )/y
gamxinflu(i,j,8,2) = (gamma_sp(alp - 0.5 + 0.5*bxhl) - gamxinf( alp - 0.5 + 0.5*bxhl, ratio ))/y
! gamxinflu(i,j,7,1) = gamxinfdp( alp - 1., ratio )/y
gamxinflu(i,j,7,1) = (gamma_dpr(alp - 1.) - gamxinfdp( alp - 1., ratio ))/y
! gamxinflu(i,j,8,1) = gamxinfdp( alp - 0.5 + 0.5*bxh, ratio )/y
gamxinflu(i,j,8,1) = (gamma_dpr(alp - 0.5 + 0.5*bxh) - gamxinfdp( alp - 0.5 + 0.5*bxh, ratio ))/y
! gamxinflu(i,j,8,2) = gamxinfdp( alp - 0.5 + 0.5*bxhl, ratio )/y
gamxinflu(i,j,8,2) = (gamma_dpr(alp - 0.5 + 0.5*bxhl) - gamxinfdp( alp - 0.5 + 0.5*bxhl, ratio ))/y
ELSE
! gamxinflu(i,j,7,1) = gamxinf( .1, ratio )/y
gamxinflu(i,j,7,1) = (gamma_sp(0.1) - gamxinf( 0.1, ratio ) )/y
! gamxinflu(i,j,8,1) = gamxinf( 1.1 - 0.5 + 0.5*bxh, ratio )/y
! gamxinflu(i,j,8,2) = gamxinf( 1.1 - 0.5 + 0.5*bxhl, ratio )/y
gamxinflu(i,j,8,1) = (gamma_sp(1.1 - 0.5 + 0.5*bxh) - gamxinf( 1.1 - 0.5 + 0.5*bxh, ratio ) )/y
gamxinflu(i,j,8,2) = (gamma_sp(1.1 - 0.5 + 0.5*bxhl) - gamxinf( 1.1 - 0.5 + 0.5*bxhl, ratio ) )/y
! gamxinflu(i,j,7,1) = gamxinfdp( .1, ratio )/y
gamxinflu(i,j,7,1) = (gamma_dpr(0.1) - gamxinfdp( 0.1, ratio ) )/y
! gamxinflu(i,j,8,1) = gamxinfdp( 1.1 - 0.5 + 0.5*bxh, ratio )/y
! gamxinflu(i,j,8,2) = gamxinfdp( 1.1 - 0.5 + 0.5*bxhl, ratio )/y
gamxinflu(i,j,8,1) = (gamma_dpr(1.1 - 0.5 + 0.5*bxh) - gamxinfdp( 1.1 - 0.5 + 0.5*bxh, ratio ) )/y
gamxinflu(i,j,8,2) = (gamma_dpr(1.1 - 0.5 + 0.5*bxhl) - gamxinfdp( 1.1 - 0.5 + 0.5*bxhl, ratio ) )/y
ENDIF

gamxinflu(i,j,7,2) = gamxinflu(i,j,7,1)
Expand All @@ -1065,14 +1069,19 @@ SUBROUTINE nssl_2mom_init( &
y7 = gamma_sp(7.+alp)
DO i = 0,nqiacrratio
ratio = float(i)*dqiacrratio
x = gamxinf( 4.+alp, ratio )

! mass fraction
x = gamxinfdp( 4.+alp, ratio )
! write(0,*) 'i, x/y = ',i, x/y
qiacrratio(i,j) = x/y
gamxinflu(i,j,4,1) = x/y
gamxinflu(i,j,4,2) = x/y

x = gamxinf( 7.+alp, ratio )
! reflectivity fraction
x = gamxinfdp( 7.+alp, ratio )
ziacrratio(i,j) = x/y7
gamxinflu(i,j,11,1) = x/y7
gamxinflu(i,j,11,2) = x/y7

ENDDO
ENDDO
Expand Down Expand Up @@ -2430,6 +2439,24 @@ REAL FUNCTION GAMMA_SP(xx)
RETURN
END FUNCTION GAMMA_SP

! #####################################################################

DOUBLE PRECISION FUNCTION GAMMA_DPR(x)
! dp gamma with real input
implicit none
real :: x
double precision :: xx

xx = x

gamma_dpr = gamma_dp(xx)

return
end FUNCTION GAMMA_DPR




! #####################################################################

real function GAMXINF(A1,X1)
Expand All @@ -2452,6 +2479,10 @@ real function GAMXINF(A1,X1)

a = a1
x = x1
IF ( x1 <= 0.0 ) THEN
gamxinf = GAMMA_SP(A1)
return
ENDIF
XAM=-X+A*DLOG(X)
IF (XAM.GT.700.0.OR.A.GT.170.0) THEN
WRITE(*,*)'a and/or x too large'
Expand Down Expand Up @@ -2509,6 +2540,10 @@ double precision function GAMXINFDP(A1,X1)

a = a1
x = x1
IF ( x1 <= 0.0 ) THEN
gamxinfdp = GAMMA_DP(A)
return
ENDIF
XAM=-X+A*DLOG(X)
IF (XAM.GT.700.0.OR.A.GT.170.0) THEN
WRITE(*,*)'a and/or x too large'
Expand Down