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
COMMON/PASQLF/ITPRO,PASQBUOY,BUOSSG,MOLEN
LOGICAL PASQBUOY,BUOSSG
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.
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)
IF(ITPRO==0.OR.ITPRO==4) THEN ! uniform reference density
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)
ELSE ! reference density varies with height
CALL GXBLIN
ENDIF
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(ITPRO>0 .AND. ITPRO/= 4) THEN ! Pasquill Reference temperature varies with height
CALL GXBLIN
ELSE ! uniform reference density
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('tem1',itorh)
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
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