c

C.... file name GXBUOYSO.F              140206
C**** SUBROUTINE GXBUOY  is called from group 13 of GREX3, and is
C     entered only when the patch name begins with the characters 'BUOY'.
C
C.... The coding provided here presumes that CO in the relevant COVAL
C     statement is FIXFLUX, and that the VAL is:
C     GRND1, or (what is equivalent) DENSITY,     (ISC = 13)
C     GRND2, or (what is equivalent) DENDIFF,     (ISC = 14)
C     GRND3, or (what is equivalent) BOUSS,       (ISC = 15)
C     GRND4, or (what is equivalent) LINBC.       (ISC = 16)
C
C     The arguments BFC, CARTES and PARAB are logicals in the SATELLITE,
C     signifying, when set to T, the body-fitted coordinate, rectangular
C     cartesian coordinate and parabolic flow respectively.
C
C.... The following library cases make use of it:
C     213,P108-P112,740-743(bfctst)                 (VAL=GRND1)
C     122,125,136,137,278,251-254,256,279           (VAL=GRND3),
C     975                                           (VAL=GRND4).
C
      SUBROUTINE GXBUOY(BFC,CARTES,PARAB,DEN1,DEN2)
      INCLUDE 'farray'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      INCLUDE 'grdear'
      INCLUDE 'grdbfc'
      INCLUDE 'prpcmn'
      COMMON/IDATA/NX,NY,IFL2(118) /LDATA/LDFL1(7),XCYCLE,LDFL2(76)
     1      /RDATA/RDAT1(10),XULAST,RDAT2(67),DVO1DT,DVO2DT,RDAT3(5)
     1      /GENL/ LGN1(51),SOLIDM,SOLIDN,LGN3(7)
     1      /GENI/ NXNY,IGNI1(8),NFM,IGFIL(39),ITEM1,ITEM2,IGFIL2(4),
     1             IPRPS,IFIL1(4)
     1      /FNDRI/IDF(11),IJPW(6),IJPS(6),IDSG(6),IDR(6),JDR(6),KDR(6)
     1      /NAMFN/NAMFUN,NAMSUB
      INTEGER DEN1,DEN2
      LOGICAL BFC,CARTES,PARAB,FLUID1,QNE,XCYCLE,LDFL1,LDFL2,SLD,
     1        TESTSOL,LGN1,SOLIDM,SOLIDN,LGN3,GRN,dbbuoy
      CHARACTER DIRCTN, NAMFUN*6, NAMSUB*6
C
      dbbuoy=.false.
cccc      dbbuoy=.true.
      if(dbbuoy) then
        write(14,*) 'gxbuoy entered for '
        call writ2i('indvar  ',indvar,'isc     ',isc)
        call writ3r('buoya   ',buoya,'buoyb   ',buoyb,'buoyc   ',buoyc)
        call writ2r('buoyd   ',buoyd,'buoye   ',buoye)
      endif
      NAMSUB = 'GXBUOY'
      IF(INDVAR.LT.5) THEN
c.... which velocity?
c                                             U1 or U2
        DIRCTN = 'E'
        BUOABC=BUOYA
        IPLPHI=NY
        IPLSLD=NY
        TESTSOL=SOLIDM
      ELSEIF(INDVAR.LT.7) THEN
c                                             V1 or V2
        DIRCTN = 'N'
        BUOABC=BUOYB
        IPLPHI=1
        IPLSLD=1
        TESTSOL=SOLIDM
      ELSE
c                                             W1 or W2
        DIRCTN = 'H'
        BUOABC=BUOYC
        IPLPHI=NFM
        IPLSLD=NXNY
        TESTSOL=SOLIDM.OR.SOLIDN
      ENDIF
      if(dbbuoy) then
        call writ40('direction is '//dirctn)
      endif
c
      IF(ISC.LT.16) THEN
        if(dbbuoy) then
          call writ40('in gxbuoyso for value')
          call writ2l('bfc     ',bfc,'cartes  ',cartes)
        endif
C---------------------------------------- VAL=DENSITY, DENSDIFF or BOUSS
C.... Resolve gravity vector into local resolute direction; BUOYA,
C     BUOYB and BUOYC signify the resolutes of the acceleration
C     due to gravity in the cartesian coordinate direction XC,
C     YC and ZC.
!!MH-start
        IF(BFC.AND..NOT.GCV) THEN
C---------------------------------------------------------------- BFC=T
          CALL BFCBUO
C....          BFCBUO sets:
C              val = easp2*buoya + easp3*buoyb + easp4*buoyb ,
C                    where the easp's are the relevant direction cosines
       ELSEIF(GCV.AND.BFC)THEN
        CALL BUOPAT
C
         RETURN
!!MH-end
        ELSEIF(CARTES) THEN
C----------------------------------------------------- BFC=F & CARTES=T
c.... fn1(y,a)                 y = a
          CALL FN1(VAL,BUOABC)
C------------------------------------------- BFC=F & CARTES=F, ie POLAR
        ELSEIF(INDVAR.GT.6) THEN
C          BUOYC is the acceleration due to gravity in the z-direction
          CALL FN1(VAL,BUOYC)
        ELSEIF(QNE(BUOYB,0.0).OR.QNE(BUOYA,0.0)) THEN
C.... The following sequence assumes that gravity is aligned with x=0,
C     so that only BUOYB ( the y-direction cartesian component of
C     the gravity vector ) is relevant.
          IF(INDVAR.LT.5) THEN
C.... The u-direction source contains:  buoya*cos(xg) - buoyb*sin(xg)
            L0XG  = L0F(LXXG)
            L0VAL = L0F(VAL)
            IXLAST= NX-1
            IF(XCYCLE) IXLAST= NX
            PI=4*ATAN(1.0)
            DO 10 IX= 1,IXLAST
              ANGLE1= F(L0XG+IX)
              IF(IX.EQ.NX) THEN
                ANGLE2= F(L0XG+1)
                IF(ANGLE1.GT.PI) ANGLE2=ANGLE2+2*PI
              ELSE
                ANGLE2= F(L0XG+IX+1)
              ENDIF
              ANGLE = 0.5*(ANGLE1+ANGLE2)
              SOURCE= COS(ANGLE)*BUOYA - SIN(ANGLE)*BUOYB
              IPLUS = (IX-1)*NY
              DO 10 IY= 1+IPLUS,NY+IPLUS
                F(L0VAL+IY)= SOURCE
  10        CONTINUE
          ELSE
C.... The v-direction source contains:  buoya*sin(xg) + buoyb*cos(xg)
            L0XG2= L0F(XG2D)
            L0VAL = L0F(VAL)
            DO 20 IX= 1,NX
              DO 20 IY= 1,NY
                I= IY+(IX-1)*NY
                ANGLE= F(L0XG2+I)
                F(L0VAL+I)= BUOYA*SIN(ANGLE) + BUOYB*COS(ANGLE)
  20        CONTINUE
          ENDIF
        ELSE
c....                       set val to zero if buoya and buoyb are zero
          CALL FN1(VAL,0.0)
        ENDIF
        if(dbbuoy) call prn('val1',val)
        if(dbbuoy) call prn('ld9 ',ld9)
        IF(ISC.EQ.14) THEN
C.... 'Phasem' force = (rho - refrho) * grav ---- VAL=DENSDIFF (ie GRND2)
C     FN69(Y,X,A,B)             Y = Y * (A + B/X)
c.... LD9 is the block-location index of the correctly-staggered density,
c     therefore no further staggering is needed below
c
c....        Note that, apart from zeroing in solids, all has now been
c            done for VAL = DENSITY and VAL = DENSDIFF,
c            The other two options require staggering.
c.... note that the do loop provides better accuracy than FN69
          L0VAL= L0F(VAL)
          L0D9 = L0F(LD9)
          DO I=1,NXNY
            F(L0VAL+I) = F(L0VAL+I)*( F(L0D9+I) - BUOYD ) / F(L0D9+I)
          ENDDO
cccc          CALL FN69(VAL,LD9,1.0,-BUOYD)
          if(dbbuoy) call prn('val2',val)
        endif
      ELSE
C....                                    ISC=16: for VAL=LINBC (ie GRND4)
C.... Value set to a linear function of variables IBUOYB and IBUOYC
        CALL FNAV(EASP2,IBUOYB,DIRCTN)
        CALL FNAV(EASP3,IBUOYC,DIRCTN)
        CALL FN10(VAL,EASP2,EASP3,BUOYA,BUOYB,BUOYC)
        if(dbbuoy) call prn('val3',val)
      ENDIF
C----------------------------------------------------------------------
C---- VAL= BOUSS  (ie GRND3)
      IF(ISC.EQ.15) THEN
C.... 'Phasem' force = rho * grav * (exco   ) * (Tref - T1)
C....             OR   rho * grav * (exco/cp) * (href - h1)
C.... Preliminaries
C....                                fn27: y = y / x
        if(dbbuoy) call writ40('boussinesq option')
        IF(FLUID1(INDVAR)) THEN
          DVDT12= DVO1DT
          ITHMEX= THRME1
          IF(SOLVE(ITEM1)) THEN
            ITORH= ITEM1
          ELSE
            ITORH= H1
C.... Divide Val by Cp when enthalpy is solved
            CALL FN27(VAL,LCP1)
          ENDIF
        ELSE
          DVDT12= DVO2DT
          ITHMEX= THRME2
          IF(SOLVE(ITEM2)) THEN
            ITORH= ITEM2
          ELSE
            ITORH= H2
C.... Divide Val by Cp when enthalpy is solved
            CALL FN27(VAL,LCP2)
          ENDIF
        ENDIF
C.... Put staggered T or H in EASP2
	IF((INDVAR.LT.7.OR..NOT.PARAB)) THEN
          CALL FNAV(EASP2,ITORH,DIRCTN)
        ELSE
          CALL FN0(EASP2,ITORH)
        ENDIF
        if(dbbuoy) then
          call writ40('staggered t or h')
          call prn('easp2',easp2)
        endif
C.... Multiply gravity by Exco * (Tref - T1) or Exco * (Href - H)
        IF(IPRPS.NE.0 .OR. GRN(DVDT12)) THEN
C.... Variable expansion coefficient is set in INDPRP and put into
C     the LDVO1 storage
          CALL INDPRP(ITHMEX,0)
          if(dbbuoy) call writ1i('idvo1   ',idvo1)
          if(dbbuoy) call prn('dvo1',idvo1)
C....                                fn85: y = y * x1 * (A - x2)
          CALL FN85(VAL, IDVO1, EASP2, BUOYE)
 
          if(dbbuoy) call prn('vall',val)
        ELSE
C....                                fn46: y = y * (A + B * x1)
          CALL FN46(VAL, EASP2, BUOYE*DVDT12, -DVDT12)
        ENDIF
      ENDIF
C.... Set VAL to zero in solids
      IF(TESTSOL) THEN
        L0VAL= L0F(VAL)
        DO 30 IX= IXF,IXL
          IADD= (IX-1)*NY
          DO 30 IY= IYF,IYL
            I= IY+IADD
            IF(SLD(I)) THEN
              F(L0VAL+I)= 0.0
            ELSE
              IF(SLD(I+IPLSLD)) F(L0VAL+I)= 0.0
            ENDIF
  30    CONTINUE
      ENDIF
      NAMSUB= 'gxbuoy'
      if(dbbuoy) write(14,*) 'leaving gxbuoy '
      END