diff --git a/bld/namelist_files/use_cases/sd_waccmx_ma_cam6.xml b/bld/namelist_files/use_cases/sd_waccmx_ma_cam6.xml
index 6b56c46b17..f380c36d60 100644
--- a/bld/namelist_files/use_cases/sd_waccmx_ma_cam6.xml
+++ b/bld/namelist_files/use_cases/sd_waccmx_ma_cam6.xml
@@ -135,4 +135,9 @@
.false.
.false.
+
+ 'SolIonRate_Tot = jeuv_1 + jeuv_2 + jeuv_3 + jeuv_4 + jeuv_5 + jeuv_6 + jeuv_7 + jeuv_8 + jeuv_9 + jeuv_10 + jeuv_11 + jeuv_14 + jeuv_15 + jeuv_16 +',
+ 'jeuv_17 + jeuv_18 + jeuv_19 + jeuv_20 + jeuv_21 + jeuv_22 + jeuv_23',
+
+
diff --git a/doc/ChangeLog b/doc/ChangeLog
index 35efde58ef..83409c55c7 100644
--- a/doc/ChangeLog
+++ b/doc/ChangeLog
@@ -81,6 +81,76 @@ izumi/gnu/aux_cam:
===============================================================
+Tag name: cam6_3_156
+Originator(s): fvitt
+Date: 16 Apr 2024
+One-line Summary: Misc code clean up for WACCM
+Github PR URL: https://github.com/ESCOMP/CAM/pull/1001
+
+Purpose of changes (include the issue number and title text for each relevant GitHub issue):
+
+ Use supported lapack library routine to solve a matrix equation in WACCM physics
+ efield module (issue #999)
+
+ Misc code clean up in calculations of effective cross section of O2
+
+ Fix for sd_waccmx_ma_cam6 use case file for waccm_mad_mam5 chemistry
+
+ Minor change to APEX module needed for when NAG compiler '-nan' option is used
+
+Describe any changes made to build system: N/A
+
+Describe any changes made to the namelist: N/A
+
+List any changes to the defaults for the boundary datasets: N/A
+
+Describe any substantial timing or memory changes: N/A
+
+Code reviewed by: cacraigucar nusbaume
+
+List all files eliminated:
+D src/chemistry/mozart/sv_decomp.F90
+ - remove deprecated matrix solve routines -- replaced by LAPACK DGESV routine
+
+List all files added and what they do: N/A
+
+List all existing files that have been modified, and describe the changes:
+M bld/namelist_files/use_cases/sd_waccmx_ma_cam6.xml
+ - fix for waccm_mad_mam5 chem
+
+M src/chemistry/mozart/mo_jshort.F90
+ - code clean up in calculations of effective cross section of O2
+
+M src/chemistry/utils/apex.F90
+ - minor changes for NAG compiler '-nan' option is used
+
+M src/physics/waccm/efield.F90
+ - use LAPACK DGESV routine to solve matrix equation
+
+If there were any failures reported from running test_driver.sh on any test
+platform, and checkin with these failures has been OK'd by the gatekeeper,
+then copy the lines from the td.*.status files for the failed tests to the
+appropriate machine below. All failed tests must be justified.
+
+derecho/intel/aux_cam:
+ PEND ERP_Ln9.C96_C96_mg17.F2000climo.derecho_intel.cam-outfrq9s_mg3
+ FAIL SMS_Ld1.f09_f09_mg17.FCHIST_GC.derecho_intel.cam-outfrq1d
+ - pre-existing failures
+
+ FAIL ERP_Ln9.f09_f09_mg17.FCSD_HCO.derecho_intel.cam-outfrq9s
+ - pre-existing failure -- should be fixed with an external cime update
+
+izumi/nag/aux_cam:
+ FAIL DAE.f45_f45_mg37.FHS94.izumi_nag.cam-dae
+ - pre-existing failure
+
+izumi/gnu/aux_cam: All PASS
+
+Summarize any changes to answers: bit-for-bit unchanged
+
+===============================================================
+===============================================================
+
Tag name: cam6_3_155
Originator(s): katec,vlarson,bstephens82,huebleruwm,zarzycki,JulioTBacmeister
Date: April 11, 2024
diff --git a/src/chemistry/mozart/mo_jshort.F90 b/src/chemistry/mozart/mo_jshort.F90
index aa47dffb31..97ec5f1375 100644
--- a/src/chemistry/mozart/mo_jshort.F90
+++ b/src/chemistry/mozart/mo_jshort.F90
@@ -71,6 +71,9 @@ module mo_jshort
real(r8), allocatable :: xs_o3b(:)
real(r8), allocatable :: xs_wl(:,:)
+ real(r8), parameter :: lno2_llimit = 38._r8 ! ln(NO2) lower limit
+ real(r8), parameter :: lno2_ulimit = 56._r8 ! ln(NO2) upper limit
+
contains
subroutine jshort_init( xs_coef_file, xs_short_file, sht_indexer )
@@ -1492,13 +1495,13 @@ subroutine calc_o2srb( nlev, nid, o2col, tlev, tsrb, xscho2 )
do k = 1,nlev
x = log( o2col(k) )
- if( x >= 38._r8 .and. x <= 56._r8 ) then
+ if( x >= lno2_llimit .and. x <= lno2_ulimit ) then
call effxs( x, tlev(k), xs )
xscho2(k,:) = xs(:)
- else if( x < 38._r8 ) then
+ else if( x < lno2_llimit ) then
ktop1 = k-1
ktop = min( ktop1,ktop )
- else if( x > 56._r8 ) then
+ else if( x > lno2_ulimit ) then
kbot = k
end if
end do
@@ -1601,9 +1604,9 @@ subroutine effxs( x, t, xs )
! method:
! ln(xs) = A(X)[T-220]+B(X)
! X = log of slant column of O2
-! A,B calculated from chebyshev polynomial coeffs
-! AC and BC using NR routine chebev. Assume interval
-! is 38lno2_ulimit) then
+ call endrun('mo_jshort::calc_params of O2 abs xs: x is not in the valid range. ')
+ end if
+
!-------------------------------------------------------------
-! ... call chebyshev evaluation routine to calc a and b from
-! set of 20 coeficients for each wavelength
+! ... evaluate at each wavelength
+! for a set of 20 Chebyshev coeficients
!-------------------------------------------------------------
do i = 1,nsrbtuv
- a(i) = jchebev( 38._r8, 56._r8, ac(1,i), 20, x )
- b(i) = jchebev( 38._r8, 56._r8, bc(1,i), 20, x )
+ a(i) = evalchebpoly( ac(:,i), x )
+ b(i) = evalchebpoly( bc(:,i), x )
end do
contains
- function jchebev( a, b, c, m, x )
-!-------------------------------------------------------------
-! Chebyshev evaluation algorithm
-! See Numerical recipes p193
-!-------------------------------------------------------------
+ ! Use Clenshaw summation algorithm to evaluate Chebyshev polynomial at point
+ ! [pnt - (lno2_ulimit + lno2_llimit)/2]/[(lno2_ulimit - lno2_llimit)/2]
+ ! given coefficients coefs within limits lim1 and lim2
+ pure function evalchebpoly( coefs, pnt ) result(cval)
+ real(r8), intent(in) :: coefs(:)
+ real(r8), intent(in) :: pnt
-!-------------------------------------------------------------
-! ... Dummy arguments
-!-------------------------------------------------------------
- integer, intent(in) :: m
- real(r8), intent(in) :: a, b, x
- real(r8), intent(in) :: c(m)
+ real(r8) :: cval
+ real(r8) :: fac(2)
+ real(r8) :: csum(2) ! Clenshaw summation
+ integer :: ndx
+ integer :: ncoef
- real(r8) :: jchebev
-!-------------------------------------------------------------
-! ... Local variables
-!-------------------------------------------------------------
- integer :: j
- real(r8) :: d, dd, sv, y, y2
+ ncoef = size(coefs)
- if( (x - a)*(x - b) > 0._r8 ) then
- write(iulog,*) 'x not in range in chebev', x
- jchebev = 0._r8
- return
- end if
+ fac(1) = (2._r8*pnt-lno2_llimit-lno2_ulimit)/(lno2_ulimit-lno2_llimit)
+ fac(2) = 2._r8*fac(1)
- d = 0._r8
- dd = 0._r8
- y = (2._r8*x - a - b)/(b - a)
- y2 = 2._r8*y
- do j = m,2,-1
- sv = d
- d = y2*d - dd + c(j)
- dd = sv
- end do
+ ! Clenshaw recurrence summation
+ csum(:) = 0.0_r8
+ do ndx = ncoef, 2, -1
+ cval = csum(1)
+ csum(1) = fac(2)*csum(1) - csum(2) + coefs(ndx)
+ csum(2) = cval
+ end do
- jchebev = y*d - dd + .5_r8*c(1)
+ cval = fac(1)*csum(1) - csum(2) + 0.5_r8*coefs(1)
- end function jchebev
+ end function evalchebpoly
end subroutine calc_params
diff --git a/src/chemistry/mozart/sv_decomp.F90 b/src/chemistry/mozart/sv_decomp.F90
deleted file mode 100644
index 0540f1f575..0000000000
--- a/src/chemistry/mozart/sv_decomp.F90
+++ /dev/null
@@ -1,364 +0,0 @@
-!-------------------------------------------------------------------------
-! purpose: singular value decomposition
-!
-! method:
-! given a matrix a(1:m,1:n), with physical dimensions mp by np,
-! this routine computes its singular value decomposition,
-! the matrix u replaces a on output. the
-! diagonal matrix of singular values w is output as a vector
-! w(1:n). the matrix v (not the transpose v^t) is output as
-! v(1:n,1:n).
-!
-! author: a. maute dec 2003
-! (* copyright (c) 1985 numerical recipes software -- svdcmp *!
-! from numerical recipes 1986 pp. 60 or can be find on web-sites
-!-------------------------------------------------------------------------
-
- module sv_decomp
-
- use shr_kind_mod, only : r8 => shr_kind_r8
-
- implicit none
-
- private
- public :: svdcmp
- public :: svbksb
-
- integer, parameter :: nmax = 1600
-
- contains
-
- subroutine svdcmp( a, m, n, mp, np, w, v )
-!-------------------------------------------------------------------------
-! ... dummy arguments
-!-------------------------------------------------------------------------
- integer, intent(in) :: m
- integer, intent(in) :: n
- integer, intent(in) :: mp
- integer, intent(in) :: np
- real(r8), intent(inout) :: a(mp,np)
- real(r8), intent(out) :: v(np,np)
- real(r8), intent(out) :: w(np)
-
-!-------------------------------------------------------------------------
-! ... local variables
-!-------------------------------------------------------------------------
- integer :: i, its, j, k, l, nm
- real(r8) :: anorm
- real(r8) :: c
- real(r8) :: f
- real(r8) :: g
- real(r8) :: h
- real(r8) :: s
- real(r8) :: scale
- real(r8) :: x, y, z
- real(r8) :: rv1(nmax)
- logical :: cnd1
- logical :: cnd2
-
- g = 0.0_r8
- scale = 0.0_r8
- anorm = 0.0_r8
-
-loop1 : &
- do i = 1,n
- l = i + 1
- rv1(i) = scale*g
- g = 0.0_r8
- s = 0.0_r8
- scale = 0.0_r8
- if( i <= m ) then
- do k = i,m
- scale = scale + abs(a(k,i))
- end do
- if( scale /= 0.0_r8 ) then
- do k = i,m
- a(k,i) = a(k,i)/scale
- s = s + a(k,i)*a(k,i)
- end do
- f = a(i,i)
- g = -sign(sqrt(s),f)
- h = f*g - s
- a(i,i) = f - g
- if( i /= n ) then
- do j = l,n
- s = 0.0_r8
- do k = i,m
- s = s + a(k,i)*a(k,j)
- end do
- f = s/h
- do k = i,m
- a(k,j) = a(k,j) + f*a(k,i)
- end do
- end do
- end if
- do k = i,m
- a(k,i) = scale*a(k,i)
- end do
- endif
- endif
- w(i) = scale *g
- g = 0.0_r8
- s = 0.0_r8
- scale = 0.0_r8
- if( i <= m .and. i /= n ) then
- do k = l,n
- scale = scale + abs(a(i,k))
- end do
- if( scale /= 0.0_r8 ) then
- do k = l,n
- a(i,k) = a(i,k)/scale
- s = s + a(i,k)*a(i,k)
- end do
- f = a(i,l)
- g = -sign(sqrt(s),f)
- h = f*g - s
- a(i,l) = f - g
- do k = l,n
- rv1(k) = a(i,k)/h
- end do
- if( i /= m ) then
- do j = l,m
- s = 0.0_r8
- do k = l,n
- s = s + a(j,k)*a(i,k)
- end do
- do k = l,n
- a(j,k) = a(j,k) + s*rv1(k)
- end do
- end do
- end if
- do k = l,n
- a(i,k) = scale*a(i,k)
- end do
- end if
- end if
- anorm = max( anorm,(abs(w(i)) + abs(rv1(i))) )
- end do loop1
-
- do i = n,1,-1
- if( i < n ) then
- if( g /= 0.0_r8 ) then
- do j = l,n
- v(j,i) = (a(i,j)/a(i,l))/g
- end do
- do j = l,n
- s = 0.0_r8
- do k = l,n
- s = s + a(i,k)*v(k,j)
- end do
- do k = l,n
- v(k,j) = v(k,j) + s*v(k,i)
- end do
- end do
- end if
- do j = l,n
- v(i,j) = 0.0_r8
- v(j,i) = 0.0_r8
- end do
- end if
- v(i,i) = 1.0_r8
- g = rv1(i)
- l = i
- end do
-
- do i = n,1,-1
- l = i + 1
- g = w(i)
- if( i < n ) then
- do j = l,n
- a(i,j) = 0.0_r8
- end do
- end if
- if( g /= 0.0_r8 ) then
- g = 1.0_r8/g
- if( i /= n ) then
- do j = l,n
- s = 0.0_r8
- do k = l,m
- s = s + a(k,i)*a(k,j)
- end do
- f = (s/a(i,i))*g
- do k = i,m
- a(k,j) = a(k,j) + f*a(k,i)
- end do
- end do
- end if
- do j = i,m
- a(j,i) = a(j,i)*g
- end do
- else
- do j = i,m
- a(j,i) = 0.0_r8
- end do
- end if
- a(i,i) = a(i,i) + 1.0_r8
- end do
-
- do k = n,1,-1
-loop2 : do its = 1,30
- do l = k,1,-1
- nm = l - 1
- cnd1 = abs( rv1(l) ) + anorm == anorm
- if( cnd1 ) then
- cnd2 = .false.
- exit
- end if
- cnd2 = abs( w(nm) ) + anorm == anorm
- if( cnd2 ) then
- cnd1 = .true.
- exit
- else if( l == 1 ) then
- cnd1 = .true.
- cnd2 = .true.
- end if
- end do
-
- if( cnd2 ) then
- c = 0.0_r8
- s = 1.0_r8
- do i = l,k
- f = s*rv1(i)
- if( (abs(f) + anorm) /= anorm ) then
- g = w(i)
- h = sqrt(f*f + g*g)
- w(i) = h
- h = 1.0_r8/h
- c = (g*h)
- s = -(f*h)
- do j = 1,m
- y = a(j,nm)
- z = a(j,i)
- a(j,nm) = (y*c) + (z*s)
- a(j,i) = -(y*s) + (z*c)
- end do
- end if
- end do
- end if
-
- if( cnd1 ) then
- z = w(k)
- if( l == k ) then
- if( z < 0.0_r8 ) then
- w(k) = -z
- do j = 1,n
- v(j,k) = -v(j,k)
- end do
- end if
- exit loop2
- end if
- end if
-
- x = w(l)
- nm = k - 1
- y = w(nm)
- g = rv1(nm)
- h = rv1(k)
- f = ((y - z)*(y + z) + (g - h)*(g + h))/(2.0_r8*h*y)
- g = sqrt( f*f + 1.0_r8 )
- f = ((x - z)*(x + z) + h*((y/(f + sign(g,f))) - h))/x
- c = 1.0_r8
- s = 1.0_r8
- do j = l,nm
- i = j + 1
- g = rv1(i)
- y = w(i)
- h = s*g
- g = c*g
- z = sqrt( f*f + h*h )
- rv1(j) = z
- c = f/z
- s = h/z
- f = (x*c)+(g*s)
- g = -(x*s)+(g*c)
- h = y*s
- y = y*c
- do nm = 1,n
- x = v(nm,j)
- z = v(nm,i)
- v(nm,j) = (x*c)+(z*s)
- v(nm,i) = -(x*s)+(z*c)
- end do
- z = sqrt( f*f + h*h )
- w(j) = z
- if( z /= 0.0_r8 ) then
- z = 1.0_r8/z
- c = f*z
- s = h*z
- end if
- f = (c*g)+(s*y)
- x = -(s*g)+(c*y)
- do nm = 1,m
- y = a(nm,j)
- z = a(nm,i)
- a(nm,j) = (y*c)+(z*s)
- a(nm,i) = -(y*s)+(z*c)
- end do
- end do
- rv1(l) = 0.0_r8
- rv1(k) = f
- w(k) = x
- end do loop2
- end do
-
- end subroutine svdcmp
-
-!-------------------------------------------------------------------------
-! purpose: solves a*x = b
-!
-! method:
-! solves a*x = b for a vector x, where a is specified by the arrays
-! u,w,v as returned by svdcmp. m and n
-! are the logical dimensions of a, and will be equal for square matrices.
-! mp and np are the physical dimensions of a. b(1:m) is the input right-hand
-! side. x(1:n) is the output solution vector. no input quantities are
-! destroyed, so the routine may be called sequentially with different b
-!
-! author: a. maute dec 2002
-! (* copyright (c) 1985 numerical recipes software -- svbksb *!
-! from numerical recipes 1986 pp. 57 or can be find on web-sites
-!-------------------------------------------------------------------------
-
- subroutine svbksb( u, w, v, m, n, mp, np, b, x )
-!-------------------------------------------------------------------------
-! ... dummy arguments
-!-------------------------------------------------------------------------
- integer, intent(in) :: m
- integer, intent(in) :: n
- integer, intent(in) :: mp
- integer, intent(in) :: np
- real(r8), intent(in) :: u(mp,np)
- real(r8), intent(in) :: w(np)
- real(r8), intent(in) :: v(np,np)
- real(r8), intent(in) :: b(mp)
- real(r8), intent(out) :: x(np)
-
-!-------------------------------------------------------------------------
-! ... local variables
-!-------------------------------------------------------------------------
- integer :: i, j, jj
- real(r8) :: s
- real(r8) :: tmp(nmax)
-
- do j = 1,n
- s = 0._r8
- if( w(j) /= 0._r8 ) then
- do i = 1,m
- s = s + u(i,j)*b(i)
- end do
- s = s/w(j)
- endif
- tmp(j) = s
- end do
-
- do j = 1,n
- s = 0._r8
- do jj = 1,n
- s = s + v(j,jj)*tmp(jj)
- end do
- x(j) = s
- end do
-
- end subroutine svbksb
-
- end module sv_decomp
diff --git a/src/chemistry/utils/apex.F90 b/src/chemistry/utils/apex.F90
index d4b60af9b1..bb690a8b42 100644
--- a/src/chemistry/utils/apex.F90
+++ b/src/chemistry/utils/apex.F90
@@ -2015,8 +2015,8 @@ subroutine cofrm(date)
! Set outputs gb(ncoef) and gv(ncoef)
! These are module data above.
!
- gb(1) = 0._r8
- gv(1) = 0._r8
+ gb(:) = 0._r8
+ gv(:) = 0._r8
f0 = -1.e-5_r8
do k=2,kmx
if (n < m) then
diff --git a/src/physics/waccm/efield.F90 b/src/physics/waccm/efield.F90
index 3ad30a970a..90508549b2 100644
--- a/src/physics/waccm/efield.F90
+++ b/src/physics/waccm/efield.F90
@@ -81,7 +81,7 @@ module efield
integer, parameter :: &
nmlon1f = nmlon/4, & ! 1 fourth mlon
nmlon2f = nmlon/2, & ! 2 fourths mlon
- nmlon3f = 3*nmlon/4 ! 3 fourths mlon
+ nmlon3f = 3*nmlon/4 ! 3 fourths mlon
real(r8) :: &
ylatm(0:nmlat), & ! magnetic latitudes (deg)
@@ -1194,7 +1194,7 @@ subroutine bnd_sinus( ihlat_bnd, itrans_width )
! Author: A. Maute Nov 2003 am 11/20/03
!----------------------------------------------------------------------------
- use sv_decomp, only : svdcmp, svbksb
+ external DGESV ! LAPACK routine to solve matrix eq
!----------------------------------------------------------------------------
! ... dummy arguments
@@ -1216,6 +1216,11 @@ subroutine bnd_sinus( ihlat_bnd, itrans_width )
real(r8) :: w(nmax_a,nmax_a)
real(r8) :: f(-nmax_sin:nmax_sin,0:nmlon)
+ real(r8) :: x(nmax_a)
+ integer :: ipiv(nmax_a), info
+
+ character(len=120) :: msg
+
!----------------------------------------------------------------------------
! Sinusoidal Boundary calculation
!----------------------------------------------------------------------------
@@ -1224,6 +1229,7 @@ subroutine bnd_sinus( ihlat_bnd, itrans_width )
u(:,:) = 0._r8
v(:,:) = 0._r8
w(:,:) = 0._r8
+ ipiv(:) = 0
do ilon = 0,nmlon ! long.
bnd = nmlath - ihlat_bnd(ilon) ! switch from pole=0 to pole =90
@@ -1238,19 +1244,18 @@ subroutine bnd_sinus( ihlat_bnd, itrans_width )
end do
end do
end do
-
-! if (debug) write(iulog,*) ' Single Value Decomposition'
- call svdcmp( u, nmax_a, nmax_a, nmax_a, nmax_a, w, v )
-
-! if (debug) write(iulog,*) ' Solving'
- call svbksb( u, w, v, nmax_a, nmax_a, nmax_a, nmax_a, rhs, lsg )
+!
+ x(:) = rhs(:)
+ call DGESV( nmax_a, 1, u, nmax_a, ipiv, x, nmax_a, info)
+ if (info/=0) then
+ write(msg,'(a,i4)') 'bnd_sinus -- LAPACK DGESV return error code: ',info
+ if (masterproc) write(iulog,*) trim(msg)
+ call endrun(trim(msg))
+ end if
+ lsg(:) = x(:)
!
do ilon = 0,nmlon ! long.
-! sum = 0._r8
sum = dot_product( lsg(-nmax_sin+ishf:nmax_sin+ishf),f(-nmax_sin:nmax_sin,ilon) )
-! do i = -nmax_sin,nmax_sin
-! sum = sum + lsg(i+ishf)*f(i,ilon)
-! end do
ihlat_bnd(ilon) = nmlath - int( sum + .5_r8 ) ! closest point
itrans_width(ilon) = int( 8._r8 - 2._r8*cos( ylonm(ilon)*dtr ) + .5_r8 )/dlatm ! 6 to 10 deg.
end do