8c8,12
< ! Compute Laplacian using Spectral Approximation.
---
> ! Compute scalar surface Laplacian:  Nabla^2 F = 1/sin(theta)*d/d(theta)*
> !  (sin(theta)*dF/d(theta)) + 1/sin^2(theta)*d^2*F/d(phi)^2 using a 
> !  truncated spectral approximation, where theta = colatitude, phi = 
> !  longitude, and F includes the constant factor 1/re^2 that is normally 
> !  present on the rhs.
22c26
<       use he_coefs_sres,only: lpmn,zmn ! he_coefs_sres.F
---
>       use he_coefs_sres,only: pmn,zmn  ! he_coefs_sres.F
24c28
<       use he_coef0_dres,only: lpmn     ! he_coef0_dres.F
---
>       use he_coef0_dres,only: pmn      ! he_coef0_dres.F
42a47,53
> #if (NLAT==36 && NLON==72)
>       integer,parameter :: truncdeg=4 ! Degree of truncation to maintain numerical stability
> #elif (NLAT==72 && NLON==144)
>       integer,parameter :: truncdeg=8 ! Degree of truncation to maintain numerical stability
> #else
>       UNKNOWN NLAT,NLON      ! compilation will stop here if unknown res
> #endif
89c100
<         do n=0,nmax-1
---
>         do n=0,truncdeg ! for no truncation, sum to nmax-1
102,112c113,124
<         ! n=nmax,m=0
<         do lat=1,nlat
<           amn(0,nmax) = amn(0,nmax)+zmn(lat,0,nmax)*fx(1,lat)
<         enddo
<         ! n=nmax,m=2,nmax-1,2
<         do m=2,nmax-1,2
<           do lat=1,nlat
<             amn(m,nmax) = amn(m,nmax)+zmn(lat,m,nmax)*fx(2*m+1,lat)
<             bmn(m,nmax) = bmn(m,nmax)+zmn(lat,m,nmax)*fx(2*m+2,lat)
<           enddo
<         enddo ! m=2,nmax-1,2
---
>         ! for (no) truncation (un)comment the following loops
>         ! n=nmax;m=0
> !        do lat=1,nlat
> !          amn(0,nmax) = amn(0,nmax)+zmn(lat,0,nmax)*fx(1,lat)
> !        enddo
>         ! n=nmax;m=2,nmax-1,2
> !        do m=2,nmax-1,2
> !          do lat=1,nlat
> !            amn(m,nmax) = amn(m,nmax)+zmn(lat,m,nmax)*fx(2*m+1,lat)
> !            bmn(m,nmax) = bmn(m,nmax)+zmn(lat,m,nmax)*fx(2*m+2,lat)
> !          enddo
> !        enddo ! m=2,nmax-1,2
121c133
<         do n=0,8 !12 !nmax-1
---
>         do n=0,truncdeg ! for no truncation, sum to nmax-1
124c136
<             fx(1,lat) = fx(1,lat)+amn(0,n)*lpmn(lat,0,n) ! A(0)
---
>             fx(1,lat) = fx(1,lat)-n*(n+1)*amn(0,n)*pmn(lat,0,n)         ! A(0)
129,130c141,144
<               fx(2*m+1,lat) = fx(2*m+1,lat)+amn(m,n)*lpmn(lat,m,n) ! A(1),A(2),...
<               fx(2*m+2,lat) = fx(2*m+2,lat)+bmn(m,n)*lpmn(lat,m,n) ! B(1),B(2),...
---
>               fx(2*m+1,lat) = fx(2*m+1,lat)-n*(n+1)*amn(m,n)*
>      |                        pmn(lat,m,n)                              ! A(1),A(2),...
>               fx(2*m+2,lat) = fx(2*m+2,lat)-n*(n+1)*bmn(m,n)*
>      |                        pmn(lat,m,n)                              ! B(1),B(2),...
133a148,162
>         ! for (no) truncation (un)comment the following loops
>         ! n=nmax;m=0
> !        do lat=lat0,lat1
> !          fx(1,lat) = fx(1,lat)-nmax*(nmax+1)*amn(0,nmax)*
> !     |                pmn(lat,0,nmax)                                   ! A(0)
> !        enddo ! lat=lat0,lat1
>         ! n=nmax;m=2,nmax-1,2
> !        do m=2,nmax-1,2
> !          do lat=lat0,lat1
> !            fx(2*m+1,lat) = fx(2*m+1,lat)-nmax*(nmax+1)*amn(m,nmax)*
> !     |                      pmn(lat,m,nmax)                             ! A(2),A(4),...
> !            fx(2*m+2,lat) = fx(2*m+2,lat)-nmax*(nmax+1)*bmn(m,nmax)*
> !     |                      pmn(lat,m,nmax)                             ! B(2),B(4),...
> !          enddo ! lat=lat0,lat1
> !        enddo ! m=2,nmax-1,2