#include "dims.h"
      SUBROUTINE HELIUM
      use cons_module,only: jmax,imax,cs,pi
      implicit none
C     ****
C     ****     SOLVE ELLIPTIC EQUATION FOR UPPER BOUNDARY CONDITION OF
C     ****     HE
C     ****
#include "params.h"
!
! This is called from advnce after lat loop and before dynamo.
! For mtask: make f and deltf arrays local, and add wt to helium
! common (see cmphe.f). 
!
      real F(ZJMXP1,ZIMXP1),DELTF(ZIMX),W(1200)
#include "helium.h"
      real :: UNITY(ZIMX),rlamda,folds,foldn,alfas,betas,pshts,alfan,
     |  betan,pshtn,bdts,bdps,bdpf,rnum,denom,dev,pertrb,bdtf,
     |  xdot(zimx),ydot(zimx)
      integer :: i,j,n,ierror
      real,external :: sddot
!     DATA UNITY/ZIMX*1./
      unity(:) = 1.
C     ****
C     ****     IF IHE=0, CALCULATE TOTAL WEIGHT
C     ****
      IF(IHE.EQ.0)THEN
	  WT=sddot(JMAX,CS,UNITY)*FLOAT(IMAX)
      ENDIF
C     ****
C     ****     CALCULATE RLAMDA=MEAN(ALFA)
C     ****
      RLAMDA=0.
      DO 1 J=1,JMAX
      RLAMDA=RLAMDA+CS(J)*sddot(IMAX,ALFA(1,J),UNITY)
    1 CONTINUE
      RLAMDA=RLAMDA/WT
C     ****
C     ****     IF IHE = 0, FOLD = CURRENT VALUE OF MBAR*T**(3./2.)*PSH
C     ****     INTERPOLATE TO POLAR GRID
C     ****
      IF(IHE.EQ.0)THEN
C     ****     CALCULATE VALUES AT POLES
	  FOLDS=sddot(IMAX,FOLD,UNITY)/FLOAT(IMAX)
	  FOLDN=sddot(IMAX,FOLD(1,JMAX),UNITY)/FLOAT(IMAX)
	  DO 2 J=JMAX,2,-1
	  DO 2 I=1,IMAX+1
	  FOLD(I,J)=.5*(FOLD(I,J)+FOLD(I,J-1))
    2     CONTINUE
	  DO 12 I=1,IMAX+1
	  FOLD(I,1)=FOLDS
	  FOLD(I,JMAX+1)=FOLDN
   12     CONTINUE
      ENDIF
C     ****
C     ****     SUBTRACT RLAMDA FROM ALFA
C     ****
      DO 3 I=1,(IMAX+1)*JMAX
      ALFA(I,1)=ALFA(I,1)-RLAMDA
    3 CONTINUE
C     ****
C     ****     INTERPOLATE ALFA AND BETA TO POLAR GRID
C     ****
C     ****     VALUES AT POLES
      ALFAS=sddot(IMAX,ALFA,UNITY)/FLOAT(IMAX)
      BETAS=sddot(IMAX,BETA,UNITY)/FLOAT(IMAX)
      PSHTS=sddot(IMAX,PSHT,UNITY)/FLOAT(IMAX)
      ALFAN=sddot(IMAX,ALFA(1,JMAX),UNITY)/FLOAT(IMAX)
      BETAN=sddot(IMAX,BETA(1,JMAX),UNITY)/FLOAT(IMAX)
      PSHTN=sddot(IMAX,PSHT(1,JMAX),UNITY)/FLOAT(IMAX)
C     ****     NORTH POLE
      DO 4 I=1,IMAX+1
      ALFA(I,JMAX+1)=ALFAN
      BETA(I,JMAX+1)=BETAN
      PSHT(I,JMAX+1)=PSHTN
    4 CONTINUE
C     ****     J = JMAX,2,-1
      DO 5 J=JMAX,2,-1
      DO 5 I=1,IMAX+1
      ALFA(I,J)=.5*(ALFA(I,J)+ALFA(I,J-1))
      BETA(I,J)=.5*(BETA(I,J)+BETA(I,J-1))
      PSHT(I,J)=.5*(PSHT(I,J)+PSHT(I,J-1))
    5 CONTINUE
C     ****     SOUTH POLE
      DO 6 I=1,IMAX+1
      ALFA(I,1)=ALFAS
      BETA(I,1)=BETAS
      PSHT(I,1)=PSHTS
    6 CONTINUE
C     ****
C     ****     ITERATION OF ELLIPTIC SOLVER
C     ****
C     ****     ITERATION COUNT
      DO 7 N=1,10
C     ****     RHS FOR SOLVER
      DO 8 J=1,JMAX+1
      DO 8 I=1,IMAX+1
      F(J,I)=BETA(I,J)-ALFA(I,J)*FOLD(I,J)
    8 CONTINUE
C     ****     CALL SOLVER
!
! hwsssp is in fishpak library:
! "Subroutine for solving a five-point finite difference approximation
!  to the Helmholtz equation in spherical coordinates and on the surface
!  of the unit sphere"
!     SUBROUTINE HWSSSP (TS,TF,M,MBDCND,BDTS,BDTF,PS,PF,N,NBDCND,BDPS,
!    1                   BDPF,ELMBDA,F,IDIMF,PERTRB,IERROR,W)
! 12/00: using source hwsssp.f (fishpak not available on ibm)
!        F array (local) is intent(inout) by hwsssp.
!        FOLD array (set from F) is in global common (helium.h)
!
      CALL HWSSSP(0.,pi,JMAX,9,BDTS,BDTF,0.,2.*pi,IMAX,
     |  0,BDPS,BDPF,RLAMDA,F,JMAX+1,PERTRB,IERROR,W)
C     ****     CHECK FOR SUCCESS
      IF(IERROR.NE.0)THEN
          write(6,"('>>> helium: HWSSP failed: ierror=',i5)") ierror
!         STOP
      ENDIF
C     ****     CALCULATE FRACTIONAL RMS DEVIATION
      RNUM=0.
      DENOM=0.
      DO 9 J=1,JMAX+1
        xdot(:) = f(j,1:imax)
        ydot(:) = f(j,1:imax)
        DO 10 I=1,IMAX
          DELTF(I)=F(J,I)-FOLD(I,J)
   10   CONTINUE
        IF(J.EQ.1)THEN
          RNUM=RNUM+.5*CS(J)*sddot(IMAX,DELTF,DELTF)
          DENOM=DENOM+.5*CS(J)*sddot(IMAX,xdot,ydot)
        ELSE IF(J.EQ.JMAX+1)THEN
          RNUM=RNUM+.5*CS(J-1)*sddot(IMAX,DELTF,DELTF)
          DENOM=DENOM+.5*CS(J-1)*sddot(IMAX,xdot,ydot)
        ELSE
          RNUM=RNUM+.5*(CS(J)+CS(J-1))*sddot(IMAX,DELTF,DELTF)
          DENOM=DENOM+.5*(CS(J)+CS(J-1))*sddot(IMAX,xdot,ydot)
        ENDIF
    9 CONTINUE
      DEV=SQRT(RNUM/DENOM)
C     ****     TRANSFER F TO FOLD
      DO 11 J=1,JMAX+1
      DO 11 I=1,IMAX+1
      FOLD(I,J)=F(J,I)
   11 CONTINUE
    7 CONTINUE
      IHE=1
      RETURN
      END