Skip to content

Commit

Permalink
Merge dom's rap GF b4b fix
Browse files Browse the repository at this point in the history
  • Loading branch information
@SamuelTrahanNOAA
SamuelTrahanNOAA committed Jun 13, 2022
2 parents 42184e7 + 77bcfb1 commit 5544dab
Show file tree
Hide file tree
Showing 2 changed files with 110 additions and 106 deletions.
214 changes: 109 additions & 105 deletions physics/cu_gf_deep.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ module cu_gf_deep
!> flag to turn off or modify mom transport by downdrafts
real(kind=kind_phys), parameter :: pgcd = 0.1
!
!> aerosol awareness, do not user yet!
!> aerosol awareness, do not use yet!
integer, parameter :: autoconv=2
integer, parameter :: aeroevap=3
real(kind=kind_phys), parameter :: scav_factor = 0.5
Expand All @@ -47,11 +47,11 @@ module cu_gf_deep

contains

integer function my_maxloc1d(A,N,dir)
integer function my_maxloc1d(A,N)
!$acc routine vector
implicit none
real(kind_phys), intent(in) :: A(:)
integer, intent(in) :: N,dir
integer, intent(in) :: N

real(kind_phys) :: imaxval
integer :: i
Expand All @@ -71,7 +71,7 @@ end function my_maxloc1d
!>\ingroup cu_gf_deep_group
!> \section general_gf_deep GF Deep Convection General Algorithm
!> @{
subroutine cu_gf_deep_run( &
subroutine cu_gf_deep_run( &
itf,ktf,its,ite, kts,kte &
,dicycle & ! diurnal cycle flag
,ichoice & ! choice of closure, use "0" for ensemble average
Expand Down Expand Up @@ -337,27 +337,28 @@ subroutine cu_gf_deep_run( &
real(kind=kind_phys), dimension (its:ite) :: &
axx,edtmax,edtmin,entr_rate
integer, dimension (its:ite) :: &
kzdown,kdet,k22,jmin,kstabi,kstabm,k22x,xland1, &
kzdown,kdet,k22,jmin,kstabi,kstabm,k22x,xland1, &
ktopdby,kbconx,ierr2,ierr3,kbmax
!$acc declare create(edt,edto,edtm,aa1,aa0,xaa0,hkb, &
!$acc hkbo,xhkb, &
!$acc xmb,pwavo,ccnloss, &
!$acc pwevo,bu,bud,cap_max, &
!$acc cap_max_increment,closure_n,psum,psumh,sig,sigd, &
!$acc axx,edtmax,edtmin,entr_rate, &
!$acc kzdown,kdet,k22,jmin,kstabi,kstabm,k22x,xland1, &
!$acc kzdown,kdet,k22,jmin,kstabi,kstabm,k22x,xland1, &
!$acc ktopdby,kbconx,ierr2,ierr3,kbmax)

integer, dimension (its:ite), intent(inout) :: ierr
integer, dimension (its:ite), intent(in) :: csum
!$acc declare copy(ierr) copyin(csum)
integer :: &
iloop,nens3,ki,kk,i,k
real(kind=kind_phys) :: &
dz,dzo,mbdt,radius,pefc, &
real(kind=kind_phys) :: &
dz,dzo,mbdt,radius, &
zcutdown,depth_min,zkbmax,z_detr,zktop, &
dh,cap_maxs,trash,trash2,frh,sig_thresh
real(kind=kind_phys) entdo,dp,subin,detdo,entup, &
real(kind=kind_phys), dimension (its:ite) :: pefc
real(kind=kind_phys) entdo,dp,subin,detdo,entup, &
detup,subdown,entdoj,entupk,detupk,totmas

real(kind=kind_phys), dimension (its:ite) :: lambau,flux_tun,zws,ztexec,zqexec
Expand Down Expand Up @@ -2269,7 +2270,7 @@ subroutine cu_gf_deep_run( &
if(ierr(i).eq.0) then
if(aeroevap.gt.1)then
! aerosol scavagening
ccnloss(i)=ccn(i)*pefc*xmb(i) ! HCB
ccnloss(i)=ccn(i)*pefc(i)*xmb(i) ! HCB
ccn(i) = ccn(i) - ccnloss(i)*scav_factor
endif
endif
Expand Down Expand Up @@ -2605,7 +2606,8 @@ subroutine cup_dd_edt(ierr,us,vs,z,ktop,kbcon,edt,p,pwav, &
real(kind=kind_phys), dimension (its:ite,1) &
,intent (out ) :: &
edtc
real(kind=kind_phys), intent (out ) :: &
real(kind=kind_phys), dimension (its:ite) &
,intent(out) :: &
pefc
real(kind=kind_phys), dimension (its:ite) &
,intent (out ) :: &
Expand Down Expand Up @@ -2639,7 +2641,7 @@ subroutine cup_dd_edt(ierr,us,vs,z,ktop,kbcon,edt,p,pwav, &
prop_c=0. !10.386
alpha3 = 0.75
beta3 = -0.15
pefc=0.
pefc(:)=0.
pefb=0.
pef=0.

Expand Down Expand Up @@ -2702,12 +2704,12 @@ subroutine cup_dd_edt(ierr,us,vs,z,ktop,kbcon,edt,p,pwav, &
prop_c=.5*(pefb+pef)/aeroadd
aeroadd=((1.e-2*ccn(i))**beta3)*(psum2(i)**(alpha3-1))
aeroadd=prop_c*aeroadd
pefc=aeroadd
pefc(i)=aeroadd

if(pefc.gt.0.9)pefc=0.9
if(pefc.lt.0.1)pefc=0.1
edt(i)=1.-pefc
if(aeroevap.eq.2)edt(i)=1.-.25*(pefb+pef+2.*pefc)
if(pefc(i).gt.0.9)pefc(i)=0.9
if(pefc(i).lt.0.1)pefc(i)=0.1
edt(i)=1.-pefc(i)
if(aeroevap.eq.2)edt(i)=1.-.25*(pefb+pef+2.*pefc(i))
endif
endif

Expand Down Expand Up @@ -4905,7 +4907,7 @@ subroutine get_zu_zd_pdf_fim(kklev,p,rand_vmas,zubeg,ipr,xland,zuh2,draft,ierr,k

if(zu(kpbli).gt.0.) &
zu(kts:min(ktf,kt-1))= zu(kts:min(ktf,kt-1))/zu(kpbli)
do k=my_maxloc1d(zu(:),kte,1),1,-1
do k=my_maxloc1d(zu(:),kte),1,-1
if(zu(k).lt.1.e-6)then
kb_adj=k+1
exit
Expand Down Expand Up @@ -4964,7 +4966,7 @@ subroutine get_zu_zd_pdf_fim(kklev,p,rand_vmas,zubeg,ipr,xland,zuh2,draft,ierr,k
! zu(kts:min(ktf,kt+1))= zu(kts:min(ktf,kt+1))/maxval(zu(kts:min(ktf,kt+1)))
if(zu(kpbli).gt.0.) &
zu(kts:min(ktf,kt-1))= zu(kts:min(ktf,kt-1))/zu(kpbli)
do k=my_maxloc1d(zu(:),kte,1),1,-1
do k=my_maxloc1d(zu(:),kte),1,-1
if(zu(k).lt.1.e-6)then
kb_adj=k+1
exit
Expand Down Expand Up @@ -5013,7 +5015,7 @@ subroutine get_zu_zd_pdf_fim(kklev,p,rand_vmas,zubeg,ipr,xland,zuh2,draft,ierr,k

if(zu(kpbli).gt.0.) &
zu(kts:min(ktf,kt-1))= zu(kts:min(ktf,kt-1))/zu(kpbli)
do k=my_maxloc1d(zu(:),kte,1),1,-1
do k=my_maxloc1d(zu(:),kte),1,-1
if(zu(k).lt.1.e-6)then
kb_adj=k+1
exit
Expand Down Expand Up @@ -5254,91 +5256,93 @@ subroutine get_inversion_layers(ierr,p_cup,t_cup,z_cup,qo_cup,qeso_cup,k_inv_lay
!$acc end parallel

end subroutine get_inversion_layers
!-----------------------------------------------------------------------------------
!>\ingroup cu_gf_deep_group
!> This function calcualtes
function deriv3(xx, xi, yi, ni, m)
!$acc routine vector
!============================================================================*/
! evaluate first- or second-order derivatives
! using three-point lagrange interpolation
! written by: alex godunov (october 2009)
! input ...
! xx - the abscissa at which the interpolation is to be evaluated
! xi() - the arrays of data abscissas
! yi() - the arrays of data ordinates
! ni - size of the arrays xi() and yi()
! m - order of a derivative (1 or 2)
! output ...
! deriv3 - interpolated value
!============================================================================*/

implicit none
integer, parameter :: n=3
integer ni, m,i, j, k, ix
real(kind=kind_phys):: deriv3, xx
real(kind=kind_phys):: xi(ni), yi(ni), x(n), f(n)

! exit if too high-order derivative was needed,
if (m > 2) then
deriv3 = 0.0
return
end if

! if x is ouside the xi(1)-xi(ni) interval set deriv3=0.0
if (xx < xi(1) .or. xx > xi(ni)) then
deriv3 = 0.0
#ifndef _OPENACC
stop "problems with finding the 2nd derivative"
#else
return
#endif
end if

! a binary (bisectional) search to find i so that xi(i-1) < x < xi(i)
i = 1
j = ni
do while (j > i+1)
k = (i+j)/2
if (xx < xi(k)) then
j = k
else
i = k
end if
end do

! shift i that will correspond to n-th order of interpolation
! the search point will be in the middle in x_i, x_i+1, x_i+2 ...
i = i + 1 - n/2

! check boundaries: if i is ouside of the range [1, ... n] -> shift i
if (i < 1) i=1
if (i + n > ni) i=ni-n+1

! old output to test i
! write(*,100) xx, i
! 100 format (f10.5, i5)

! just wanted to use index i
ix = i
! initialization of f(n) and x(n)
do i=1,n
f(i) = yi(ix+i-1)
x(i) = xi(ix+i-1)
end do

! calculate the first-order derivative using lagrange interpolation
if (m == 1) then
deriv3 = (2.0*xx - (x(2)+x(3)))*f(1)/((x(1)-x(2))*(x(1)-x(3)))
deriv3 = deriv3 + (2.0*xx - (x(1)+x(3)))*f(2)/((x(2)-x(1))*(x(2)-x(3)))
deriv3 = deriv3 + (2.0*xx - (x(1)+x(2)))*f(3)/((x(3)-x(1))*(x(3)-x(2)))
! calculate the second-order derivative using lagrange interpolation
else
deriv3 = 2.0*f(1)/((x(1)-x(2))*(x(1)-x(3)))
deriv3 = deriv3 + 2.0*f(2)/((x(2)-x(1))*(x(2)-x(3)))
deriv3 = deriv3 + 2.0*f(3)/((x(3)-x(1))*(x(3)-x(2)))
end if
end function deriv3
! DH* 20220604 - this isn't used at all
!!!!-----------------------------------------------------------------------------------
!!!!>\ingroup cu_gf_deep_group
!!!!> This function calcualtes
!!! function deriv3(xx, xi, yi, ni, m)
!!!!$acc routine vector
!!! !============================================================================*/
!!! ! evaluate first- or second-order derivatives
!!! ! using three-point lagrange interpolation
!!! ! written by: alex godunov (october 2009)
!!! ! input ...
!!! ! xx - the abscissa at which the interpolation is to be evaluated
!!! ! xi() - the arrays of data abscissas
!!! ! yi() - the arrays of data ordinates
!!! ! ni - size of the arrays xi() and yi()
!!! ! m - order of a derivative (1 or 2)
!!! ! output ...
!!! ! deriv3 - interpolated value
!!! !============================================================================*/
!!!
!!! implicit none
!!! integer, parameter :: n=3
!!! integer ni, m,i, j, k, ix
!!! real(kind=kind_phys):: deriv3, xx
!!! real(kind=kind_phys):: xi(ni), yi(ni), x(n), f(n)
!!!
!!! ! exit if too high-order derivative was needed,
!!! if (m > 2) then
!!! deriv3 = 0.0
!!! return
!!! end if
!!!
!!! ! if x is ouside the xi(1)-xi(ni) interval set deriv3=0.0
!!! if (xx < xi(1) .or. xx > xi(ni)) then
!!! deriv3 = 0.0
!!!#ifndef _OPENACC
!!! stop "problems with finding the 2nd derivative"
!!!#else
!!! return
!!!#endif
!!! end if
!!!
!!! ! a binary (bisectional) search to find i so that xi(i-1) < x < xi(i)
!!! i = 1
!!! j = ni
!!! do while (j > i+1)
!!! k = (i+j)/2
!!! if (xx < xi(k)) then
!!! j = k
!!! else
!!! i = k
!!! end if
!!! end do
!!!
!!! ! shift i that will correspond to n-th order of interpolation
!!! ! the search point will be in the middle in x_i, x_i+1, x_i+2 ...
!!! i = i + 1 - n/2
!!!
!!! ! check boundaries: if i is ouside of the range [1, ... n] -> shift i
!!! if (i < 1) i=1
!!! if (i + n > ni) i=ni-n+1
!!!
!!! ! old output to test i
!!! ! write(*,100) xx, i
!!! ! 100 format (f10.5, i5)
!!!
!!! ! just wanted to use index i
!!! ix = i
!!! ! initialization of f(n) and x(n)
!!! do i=1,n
!!! f(i) = yi(ix+i-1)
!!! x(i) = xi(ix+i-1)
!!! end do
!!!
!!! ! calculate the first-order derivative using lagrange interpolation
!!! if (m == 1) then
!!! deriv3 = (2.0*xx - (x(2)+x(3)))*f(1)/((x(1)-x(2))*(x(1)-x(3)))
!!! deriv3 = deriv3 + (2.0*xx - (x(1)+x(3)))*f(2)/((x(2)-x(1))*(x(2)-x(3)))
!!! deriv3 = deriv3 + (2.0*xx - (x(1)+x(2)))*f(3)/((x(3)-x(1))*(x(3)-x(2)))
!!! ! calculate the second-order derivative using lagrange interpolation
!!! else
!!! deriv3 = 2.0*f(1)/((x(1)-x(2))*(x(1)-x(3)))
!!! deriv3 = deriv3 + 2.0*f(2)/((x(2)-x(1))*(x(2)-x(3)))
!!! deriv3 = deriv3 + 2.0*f(3)/((x(3)-x(1))*(x(3)-x(2)))
!!! end if
!!! end function deriv3
! *DH 20220604
!=============================================================================================
!>\ingroup cu_gf_deep_group
subroutine get_lateral_massflux(itf,ktf, its,ite, kts,kte &
Expand Down
2 changes: 1 addition & 1 deletion physics/cu_gf_driver.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ module cu_gf_driver
! DH* TODO: replace constants with arguments to cu_gf_driver_run
!use physcons , g => con_g, cp => con_cp, xlv => con_hvap, r_v => con_rv
use machine , only: kind_phys
use cu_gf_deep, only: cu_gf_deep_run,neg_check,autoconv,aeroevap,fct1d3
use cu_gf_deep, only: cu_gf_deep_run,neg_check,fct1d3
use cu_gf_sh , only: cu_gf_sh_run

implicit none
Expand Down

0 comments on commit 5544dab

Please sign in to comment.