! module bndry_module ! ! Lower boundary conditions: diurnal, semi-diurnal, annual are ! lower boundary is only set in these subroutines if no input file ! from the GSWM is provided. Therefore a check is included for ! the tidal contribution (diurnal, semi-diurnal) ! ! ----------------------------------------------------------------- ! Titan TGCM modifications. Latest version: 9/15/04 S. W. BOUGHER ! (BNDEF mods mostly) ! ----------------------------------------------------------------- ! use cons_module,only: t0,pi,atm_amu,gask,grav,freq_semidi, | re,dlamda,tbound,tgrad,cs,cor,tn,freq_ann use params_module,only: nlat,nlonp4,nlevp1,dz implicit none ! complex,parameter :: ci=(0.,1.) complex :: | zb(nlat),zb2(nlat),zba(nlat), | tb(nlat),tb2(nlat),tba(nlat), | ub(nlat),ub2(nlat),uba(nlat), | vb(nlat),vb2(nlat),vba(nlat), | bnd(nlonp4),bnd2(nlonp4),bnda(nlonp4) ! real :: b(nlonp4,2,2),fb(nlonp4,2) ! contains !----------------------------------------------------------------------- subroutine lowbound ! use init_module,only: igetgswmdi,igetgswmsdi ! GSWM input flags ! Call the lower boundary routines. This is called once per model run ! from main tgcm.F. ! ! Local: integer :: k ! ! tbound = value of T at lower boundary, and tgrad = gradient ! write(6,"(/,72('-'))") write(6,"('Set Lower Boundary Conditions:')") t0(1) = tbound t0(2) = tbound+dz*tgrad if(igetgswmsdi == 0) call bndry_semidiurnal ! check for GSWM semidiurnal if(igetgswmdi == 0) call bndry_diurnal ! check for GSWM diurnal call bndry_annual write(6,"(72('-'),/)") ! ! t0==0: do k=1,nlevp1 t0(K)=0. enddo end subroutine lowbound !----------------------------------------------------------------------- subroutine bndry_semidiurnal ! ! Lower boundary conditions for semi-diurnal tide: ! 2/00: 1998 spherepack lib code (sphpac.f) replaces old lib ! alfpac.f for legendre polynomials and Hough functions. ! This routine calculates complex ZB,TB,UB,VB. ! ! This will be MUCH different for Titan. Leave code for now. ! Set to zeros for lexical reads. ! use input_module,only: tide ! ! Local: integer,parameter :: nalf=19, malf=2 real :: p(nlat,nalf,malf),hough(nlat,5,malf),cp(nalf/2+1) complex :: dzb(nlat) real :: B(5,19),RL(5),BHOUR(5),rlamda,xdot(19),ydot(19), | ptscal,theta,ptjm(2*nlat+1) integer :: n,jm,l,lm1,m,mm1,j,ld,i,nm1 ! COMPLEX ZEE(5),CL(5),EXPDLM DATA B/ 2 0.969152, 0.0 , 0.216046, 0.0 , 0.093838, 3 0.0 , 0.909763, 0.0 , 0.342113, 0.0 , 4-0.245226, 0.0 , 0.798445, 0.0 , 0.421218, 5 0.0 ,-0.408934, 0.0 , 0.645517, 0.0 , 6 0.024633, 0.0 ,-0.543993, 0.0 , 0.464159, 7 0.0 , 0.071127, 0.0 ,-0.643189, 0.0 , 8-0.001292, 0.0 , 0.139613, 0.0 ,-0.699495, 9 0.0 ,-0.006673, 0.0 , 0.225090, 0.0 , 1 0.000042, 0.0 ,-0.019654, 0.0 , 0.320141, 1 0.0 , 0.000394, 0.0 ,-0.043345, 0.0 , 2-0.000001, 0.0 , 0.001772, 0.0 ,-0.079831, 3 0.0 ,-0.000016, 0.0 , 0.005401, 0.0 , 4 0.0 , 0.0 ,-0.000112, 0.0 , 0.012932, 5 0.0 , 0.0 , 0.0 ,-0.000476, 0.0 , 6 0.0 , 0.0 , 0.000005, 0.0 ,-0.001490, 7 0.0 , 0.0 , 0.0 , 0.000031, 0.0 , 8 0.0 , 0.0 , 0.0 , 0.0 , 0.000129, 9 0.0 , 0.0 , 0.0 ,-0.000002, 0.0 , 1 0.0 , 0.0 , 0.0 , 0.0 ,-0.000009/ DATA RL/7.8519E5, 3.6665E5, 2.1098E5, 1.3671E5, 0.9565E5/ real,external :: sddot ! in util.F ! bhour = tide(6:10) if (all(tide==0.)) goto 13 DO N=1,5 ZEE(N)=TIDE(N)*CEXP(CI*pi*BHOUR(N)/6.) CL(N)=CSQRT(CMPLX(gask/(atm_amu*grav*RL(N))* | (T0(1)*2./7.+(T0(2)-T0(1))/dz)-.25))-.5*CI enddo JM=2*nlat+1 ! C **** SET UP HOUGH FUNCTIONS ! ! Using new (1998) spherepack (sphpac.f): do n=2,nalf+1 nm1 = n-1 do m=2,malf+1 mm1=m-1 call alfk(n,m,cp) do j=1,jm theta = float(j-1)*pi/float(jm-1) call lfpt(n,m,theta,cp,ptscal) ptjm(j) = ptscal enddo do j=1,nlat p(j,nm1,mm1) = ptjm(2*(nlat+1-j)) enddo enddo DO J=1,nlat P(J,nm1,2)=SQRT(FLOAT(n*(n+1)-6))*P(J,nm1,2)-2.*TN(J)* | P(J,nm1,1) enddo enddo ! ! util.F: real function sddot(n,x,y) DO L=1,5 DO LD=1,2 DO J=1,nlat xdot(:) = p(j,:,ld) ydot(:) = b(l,:) HOUGH(J,L,LD)=sddot(19,xdot,ydot) enddo enddo enddo C **** GENERATE ZB, UB, VB AND TB DO 5 J=1,nlat TB(J)=0. ZB(J)=0. DZB(J)=0. 5 CONTINUE DO L=1,5 DO J=1,nlat ZB(J)=ZB(J)+ZEE(L)*HOUGH(J,L,1) DZB(J)=DZB(J)+ZEE(L)*HOUGH(J,L,2) TB(J)=TB(J)+CI*atm_amu*grav/gask*ZEE(L)*CL(L)*HOUGH(J,L,1) enddo enddo DO 7 J=1,nlat UB(J)=freq_semidi*re*(1.-(COR(J)/freq_semidi)**2) VB(J)=CI*grav*(DZB(J)-2.*COR(J)/(freq_semidi*CS(J))*ZB(J))/UB(J) UB(J)=grav*(COR(J)/freq_semidi*DZB(J)-2./CS(J)*ZB(J))/UB(J) 7 CONTINUE GO TO 11 13 CONTINUE C **** ZERO BOUNDARY CONDITION DO 12 J=1,nlat ZB(J)=0. TB(J)=0. UB(J)=0. VB(J)=0. 12 CONTINUE 11 CONTINUE C **** CALCULATE LONGITUDINAL STRUCTURE RLAMDA = -2.*dlamda BND(1)=CEXP(CI*2.*RLAMDA) EXPDLM=CEXP(CI*2.*dlamda) DO 9 I=2,nlonp4 BND(I)=BND(I-1)*EXPDLM 9 CONTINUE end subroutine bndry_semidiurnal !----------------------------------------------------------------------- subroutine bndry_diurnal ! C **** TIDAL BOUNDARY CONDITION FOR DIURNAL GRAVITATIONAL MODE C **** (1,1) ! 2/00: 1998 spherepack lib code (sphpac.f) replaces old lib ! alfpac.f for legendre polynomials and Hough functions. ! This routine calculates complex ZB2,TB2,UB2,VB2. ! ! This will be MUCH different for Titan. Leave code for now. ! Set to zeros for lexical reads. ! use input_module,only: tide2 ! ! Local: integer,parameter :: nalf=19, malf=2 real :: p(nlat,nalf,malf),hough(nlat,5,malf),cp(nalf/2+1) complex :: dzb(nlat) real :: B(1,19),RL (1),BHOUR2(1),rlamda,xdot(19),ydot(19), | ptscal,theta,ptjm(2*nlat+1),pik integer :: l,m,j,n,jm,ld,i COMPLEX ZEE(1),CL(1),EXPDLM ! DATA B/ 2 0.282710, 3 0.0 , 4-0.638229, 5 0.0 , 6 0.620521, 7 0.0 , 8-0.336408, 9 0.0 , 1 0.117021, 1 0.0 , 2-0.028332, 3 0.0 , 4 0.005042, 5 0.0 , 6-0.000686, 7 0.0 , 8 0.000074, 9 0.0 , 1-0.000006/ DATA RL/0.6909E5/ real,external :: sddot ! in util.F ! bhour2(1) = tide2(2) ! ! pik is as per tgcm15 (4*atan(1)), without -lmass (pi with and without ! -lmass differs in 16-17 decimal places, apparently causing diffs in zee, ! and subsequently zb2). This does not seem to be a problem with semidiurnal. ! pik = 3.14159265358979312 if (all(tide2==0.)) goto 13 DO 1 N=1,1 ZEE(N)=TIDE2(N)*CEXP(CI*pik*BHOUR2(N)/12.) CL(N)=CSQRT(CMPLX(gask/(atm_amu*grav*RL(N))* | (T0(1)*2./7.+(T0(2)-T0(1))/dz)-.25))-.5*CI 1 CONTINUE JM=2*nlat+1 ! C **** SET UP HOUGH FUNCTIONS ! ! Using new (1998) spherepack (sphpac.f): do n=1,19 do m=1,2 call alfk(n,m,cp) do j=1,jm theta = float(j-1)*pi/float(jm-1) call lfpt(n,m,theta,cp,ptscal) ptjm(j)=ptscal enddo do j=1,nlat p(j,n,m) = ptjm(2*(nlat+1-j)) enddo enddo DO J=1,nlat P(J,n,2)=SQRT(FLOAT(n*(n+1)-2))*P(J,n,2)-TN(J)*P(J,n,1) enddo enddo ! ! util.F: real function sddot(n,x,y) DO L=1,1 DO LD=1,2 DO J=1,nlat xdot(:) = p(j,:,ld) ydot(:) = b(l,:) HOUGH(J,L,LD)=sddot(19,xdot,ydot) enddo enddo enddo C **** GENERATE ZB2, UB2, VB2 AND TB2 DO 5 J=1,nlat TB2(J)=0. ZB2(J)=0. DZB(J)=0. 5 CONTINUE DO 6 L=1,1 DO 6 J=1,nlat ZB2(J)=ZB2(J)+ZEE(L)*HOUGH(J,L,1) DZB(J)=DZB(J)+ZEE(L)*HOUGH(J,L,2) TB2(J)=TB2(J)+CI*atm_amu*grav/gask*ZEE(L)*CL(L)*HOUGH(J,L,1) 6 CONTINUE DO 7 J=1,nlat UB2(J)=.5*freq_semidi*re*(1.-(COR(J)/(.5*freq_semidi))**2) VB2(J)=CI*grav*(DZB(J)-COR(J)/(.5*freq_semidi*CS(J))*ZB2(J))/ | UB2(J) UB2(J)=grav*(COR(J)/(.5*freq_semidi)*DZB(J)-1./CS(J)*ZB2(J))/ | UB2(J) 7 CONTINUE GO TO 11 13 CONTINUE C **** ZERO BOUNDARY CONDITION DO 12 J=1,nlat ZB2(J)=0. TB2(J)=0. UB2(J)=0. VB2(J)=0. 12 CONTINUE 11 CONTINUE C **** CALCULATE LONGITUDINAL STRUCTURE RLAMDA = -2.*dlamda BND2(1)=CEXP(CI*RLAMDA) EXPDLM=CEXP(CI*dlamda) DO 9 I=2,nlonp4 BND2(I)=BND2(I-1)*EXPDLM 9 CONTINUE end subroutine bndry_diurnal !----------------------------------------------------------------------- subroutine bndry_annual ! ! 2/00: 1998 spherepack lib code (sphpac.f) replaces old lib ! alfpac.f for legendre polynomials and Hough functions. ! This routine calculates complex ZBA,TBA,UBA,VBA. ! ! This will be MUCH different for Titan. Leave code for now. ! Set to zeros for lexical reads. ! C **** TIDAL BOUNDARY CONDITION FOR ANNUAL MODE ! use input_module,only: tideann use init_module,only: iday ! ! For 1998 spherepack lib code (sphpac.f) ! (replaces old alfpac.f) ! integer,parameter :: nalf=24, malf=2 ! ! Local: real :: p(nlat,nalf,malf),hough(nlat,0:6,2),cp(nalf/2+1) complex :: dzb(nlat),zzb(nlat) real :: B(6,24),RL(0:6),scale,rt2,rm,factor,xdot(24),ydot(24), | pi,ptscal,ptjm(2*nlat+1),theta,w(nlat) integer :: n,l,m,k,i,jm,mp1,ld,j COMPLEX CC(0:6,0:6),CL(0:6),EXPT ! DATA ((B(I,J),I = 1,6),J = 1,12)/ 1-0.882922, 0.000000,-0.345087, 0.000000,-0.202228, 0.000000, 2 0.000000,-0.930826, 0.000000,-0.301357, 0.000000,-0.152720, 3-0.466226, 0.000000, 0.567457, 0.000000, 0.407114, 0.000000, 4 0.000000,-0.362673, 0.000000, 0.694431, 0.000000, 0.438014, 5-0.055436, 0.000000, 0.711847, 0.000000,-0.163050, 0.000000, 6 0.000000,-0.044983, 0.000000, 0.625545, 0.000000,-0.325772, 7-0.002909, 0.000000, 0.225723, 0.000000,-0.749160, 0.000000, 8 0.000000,-0.002773, 0.000000, 0.186467, 0.000000,-0.723674, 9-0.000086, 0.000000, 0.034940, 0.000000,-0.435919, 0.000000, * 0.000000,-0.000103, 0.000000, 0.029425, 0.000000,-0.379254, 1-0.000002, 0.000000, 0.003267, 0.000000,-0.122687, 0.000000, 2 0.000000,-0.000003, 0.000000, 0.002928, 0.000000,-0.104008/ DATA ((B(I,J),I = 1,6),J = 13,24)/ 3 0.0 , 0.000000, 0.000206, 0.000000,-0.021267, 0.000000, 4 0.0 , 0.0 , 0.000000, 0.000202, 0.000000,-0.018228, 5 0.0 , 0.0 , 0.000009, 0.000000,-0.002540, 0.000000, 6 0.0 , 0.0 , 0.000000, 0.000010, 0.000000,-0.002252, 7 0.0 , 0.0 , 0.0 , 0.000000,-0.000223, 0.000000, 8 0.0 , 0.0 , 0.0 , 0.0 , 0.000000,-0.000208, 9 0.0 , 0.0 , 0.0 , 0.0 ,-0.000015, 0.000000, * 0.0 , 0.0 , 0.0 , 0.0 , 0.000000,-0.000015, 1 0.0 , 0.0 , 0.0 , 0.0 ,-0.000001, 0.000000, 2 0.0 , 0.0 , 0.0 , 0.0 , 0.000000,-0.000001, 3 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.000000, 4 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 / DATA (RL(N),N=1,6)/ 1 -10.8409E5,-7.0243E5,-2.4874E5,-1.9696E5,-1.0694E5, 2 -0.9119E5/ C **** C **** ANNUAL BOUNDARY COEFFICIENTS FOR LOWER BOUNDARY AT C **** Z = -7. C **** DATA((CC(K,N),N=0,6),K=0,3)/ 1 ( 0.136291E+03, 0.000000E+00),(-0.200536E-01, 0.000000E+00), 2 ( 0.423456E+00, 0.000000E+00),(-0.143623E-02, 0.000000E+00), 3 ( 0.262889E+00, 0.000000E+00),( 0.365122E-02, 0.000000E+00), 4 ( 0.102716E+00, 0.000000E+00),( 0.826209E+00,-0.116536E-01), 5 (-0.455993E+00,-0.311380E-02),( 0.271258E-01,-0.373420E-01), 6 (-0.433095E-01, 0.165546E-01),( 0.713386E-02,-0.112538E-01), 7 (-0.153849E-01, 0.103063E-01),( 0.214366E-02,-0.570878E-02), 8 (-0.316298E+00, 0.115053E+00),(-0.159072E-01,-0.245495E-02), 9 ( 0.302211E+00,-0.132446E-01),( 0.230750E-02, 0.170566E-03), * ( 0.100434E+00,-0.299227E-02),( 0.264555E-02,-0.137723E-03), 1 ( 0.499098E-01,-0.110255E-02),(-0.520584E-03,-0.114124E-02), 2 ( 0.178599E-01, 0.561092E-02),( 0.557591E-03, 0.176165E-02), 3 (-0.148151E-02, 0.749397E-03),(-0.617325E-03, 0.778743E-03), 4 (-0.530835E-03, 0.641767E-03),(-0.964206E-03, 0.551394E-03)/ DATA((CC(K,N),N=0,6),K=4,6)/ 1 (-0.137927E-02, 0.866386E-03),(-0.242825E-02, 0.441184E-03), 2 ( 0.120715E-02,-0.136729E-02),( 0.122657E-03, 0.316213E-04), 3 ( 0.390769E-03,-0.162978E-03),( 0.378377E-03,-0.195668E-04), 4 ( 0.366912E-03,-0.681579E-04),(-0.470068E-03,-0.118650E-05), 5 ( 0.120025E-02,-0.797459E-03),( 0.622700E-03,-0.424648E-04), 6 (-0.537275E-03, 0.101658E-03),( 0.222407E-03,-0.828812E-05), 7 (-0.209097E-03, 0.828365E-04),( 0.945940E-04, 0.317248E-04), 8 ( 0.341903E-03, 0.192246E-04),( 0.129833E-03, 0.247156E-04), 9 (-0.610206E-03, 0.591081E-06),(-0.102160E-03,-0.434110E-04), * ( 0.196672E-04,-0.305687E-04),(-0.905354E-04,-0.813929E-04), 1 ( 0.569460E-05,-0.116661E-03)/ C****************************** DATA SCALE/1./ C****************************** real,external :: sddot ! in util.F ! RT2 = SQRT(2.) if (tideann==0) goto 13 JM = 2*nlat+1 C **** C **** HEIGHT VARIATION C **** DO 1 N = 1,6 CL(N) = -CSQRT(CMPLX(gask/(atm_amu*grav*RL(N))* 1 (T0(1)*2./7.+(T0(2)-T0(1))/dz)-.25))-.5*CI 1 CONTINUE CL(0) = 0. C **** C **** SET UP LEGENDRE POLYNOMIALS C **** ! ! Using new (1998) spherepack (sphpac.f): ! (nalf=24, see cbndrya.h) ! pi = 4.*atan(1.) do n=1,nalf do m=0,1 mp1=m+1 call alfk(n,m,cp) do j=1,jm theta = float(j-1)*pi/float(jm-1) call lfpt(n,m,theta,cp,ptscal) ptjm(j) = ptscal enddo do j=1,nlat p(j,n,mp1) = ptjm(2*(nlat+1-j)) enddo enddo m = 0 rm = float(m) do j=1,nlat p(j,n,2)=sqrt(float(n*(n+1)-m*(m+1)))*p(j,n,2)-rm*tn(j)* | p(j,n,1) enddo enddo C **** C **** NOW EVALUATE HOUGH FUNCTIONS C **** DO 4 L=1,6 DO 4 LD=1,2 DO 4 J=1,nlat xdot(:) = p(j,:,ld) ydot(:) = b(l,:) HOUGH(J,L,LD)=sddot(24,xdot,ydot) 4 CONTINUE C **** C **** HOUGH FUNCTION OF ORDE ZERO` C **** DO 5 J =1,nlat HOUGH(J,0,1) = 1./RT2 HOUGH(J,0,2) = 0. 5 CONTINUE C **** C **** GENERATE ZBA, TBA, UBA, VBA C **** DO 6 J = 1,nlat ZBA(J) = 0. TBA(J) = 0. UBA(J) = 0. VBA(J) = 0. 6 CONTINUE C **** C **** SUMMATION OVER FREQUENCY, K, CALCULATION OF PHASE FACTOR C **** DO 7 K = 0,6 EXPT = CEXP(CI*FLOAT(K*(IDAY-1))*86400.*freq_ann) C **** C **** SUMMATION OVER ORDER, N C **** DO 7 N = 0,6 FACTOR = SCALE IF(K.EQ.0.AND.N.EQ.0)FACTOR = 1. DO 7 J = 1,nlat ZZB(J) = CC(K,N)*HOUGH(J,N,1)*EXPT*1.E5*FACTOR DZB(J) = CC(K,N)*HOUGH(J,N,2)*EXPT*1.E5*FACTOR W(J) = grav/(re*((FLOAT(K)*freq_ann)**2-COR(J)**2)) ZBA(J) = ZBA(J)+ZZB(J) TBA(J) = TBA(J)+CI*atm_amu*grav/gask*CL(N)*ZZB(J) UBA(J) = UBA(J)+W(J)*(COR(J)*DZB(J)-RM*FLOAT(K)*freq_ann/ 1 CS(J)*ZZB(J)) VBA(J) = VBA(J)+CI*W(J)*(FLOAT(K)*freq_ann*DZB(J)-RM*COR(J)/ 1 CS(J)*ZZB(J)) 7 CONTINUE GO TO 11 13 CONTINUE C **** C **** ZERO BOUNDARY CONDITION C **** DO 12 J = 1,nlat ZBA(J) = CC(0,0)*1./RT2*1.E5 TBA(J) = 0. UBA(J) = 0. VBA(J) = 0. 12 CONTINUE 11 CONTINUE C **** C **** LONGITUDINAL STRUCTURE C **** DO 8 I = 1,nlonp4 BNDA(I) = 1. 8 CONTINUE end subroutine bndry_annual !----------------------------------------------------------------------- subroutine bndcmp ! BNDEF mods mostly) C **** C **** CALCULATE MATRICES B(nlonp4,2,2) AND VECTORS FB(nlonp4,2) C **** REPRESENTING THE LOWER BOUBNDARY CONDITION IN COMP, C **** WHERE PSI1 AND PSI2 ARE CALCULATED: C **** C **** PSI(K=-1/2) = B * PSI(K=1/2) + FB C **** C **** BNDCMP CALLS THE SUBROUTINE BNDEF TO DEFINE THE 2 X 2 C **** MATRICES E, F AND THE 2 VECTOR G IN THE GENERAL C **** LOWER BOUNDARY CONDITION: C **** C **** E * D(PSI)/DS + F * PSI + G = 0. C **** C **** WHERE: C **** PSI = |PSI1| AND THE BOUNDARY CONDITION IS APPLIED C **** | | C **** |PSI2| C **** C **** AT LEVEL ZERO C **** C **** THIS SUBROUTINE THEN EVALUATES B AND FB FROM: C **** C **** B = (E/DS - F/2.)**(-1) * (E/DS + F/2.) C **** C **** FB = (E/DS - F/2.)**(-1) * G C **** ! ! Local: real :: ee(nlonp4,2,2),ff(nlonp4,2,2),gg(nlonp4,2), | wm1(nlonp4,2,2),wm2(nlonp4,2,2),wm3(nlonp4,2,2), | ws1(nlonp4) integer :: l,m,i C **** C **** CALL BNDEF TO DEFINE E, F AND G IN S1, S2 AND S3 C **** call bndef(ee,ff,gg) C **** C **** WM1 = (E/DS - F/2.) C **** C **** WM1 = (E/DS + F/2.) C **** do l = 1,2 do m = 1,2 do i = 1,nlonp4 wm1(i,l,m) = ee(i,l,m)/dz-ff(i,l,m)/2. wm2(i,l,m) = ee(i,l,m)/dz+ff(i,l,m)/2. enddo enddo enddo C **** C **** WM3 = WM1**(-1) C **** C **** WS1 = DET(WM1) C **** do i = 1,nlonp4 ws1(i) = wm1(i,1,1)*wm1(i,2,2)-wm1(i,1,2)*wm1(i,2,1) enddo C **** C **** NOW INVERSE OF WM1 IN WM3 C **** do i = 1,nlonp4 wm3(i,1,1) = wm1(i,2,2)/ws1(i) wm3(i,1,2) = -wm1(i,1,2)/ws1(i) wm3(i,2,1) = -wm1(i,2,1)/ws1(i) wm3(i,2,2) = wm1(i,1,1)/ws1(i) enddo C **** C **** B = WM3 * WM2 C **** ! b and fb are bndry_module module data. do l = 1,2 do m = 1,2 do i = 1,nlonp4 b(i,l,m) = wm3(i,l,1)*wm2(i,1,m)+wm3(i,l,2)*wm2(i,2,m) enddo enddo enddo C **** C **** FB = WM3 * G C **** do l = 1,2 do i = 1,nlonp4 fb(i,l) = wm3(i,l,1)*gg(i,1)+wm3(i,l,2)*gg(i,2) enddo enddo end subroutine bndcmp !----------------------------------------------------------------------- SUBROUTINE BNDEF(ee,ff,gg) C **** C **** BNDEF DEFINES THE LOWER BOUNDARY CONDITION FOR THIS C **** VERSION OF THE MODEL C **** C **** THE LOWER BOUNDARY CONDITION FOR COMP IS: C **** C **** E * D(PSI)/DS + F * PSI +G = 0. C **** C **** WHERE: C **** PSI = VECTOR(PSI1,PSI2) C **** E AND F ARE 2 X 2 MATRICES C **** G = VECTOR(G1,G2) C **** C **** E, F AND G MAY BE FUNCTIONS OF LATITUDE & LONGITUDE C **** C **** THIS SUBROUTINE DEFINES E, F AND G for BNDCMP C **** ! ! Args: real,intent(out) :: EE(nlonp4,2,2),FF(nlonp4,2,2),GG(nlonp4,2) ! ! Local: real,parameter :: alfa1=9.676E-03, alfa2=7.1E-05 integer :: i C **** C **** IN TIGCM AND TIEGCM: C **** C **** E = |0. 0.| C **** | | C **** |0. 1.| C **** C **** F = |1. 1.| C **** | | C **** |0. -1.| C **** C **** G = |-ALFA| C **** | | C **** | 0.| C **** C **** WHERE: C **** ALFA = 0.22 + 0.014 = 0.234 C **** C **** C **** IN VTGCM, MTGCM, TTGCM: NO/LITTLE DATA AT L.BOUNDARY C **** C **** E = |0. 0.| C **** | | C **** |0. 0.| C **** C **** F = |1. 0.| C **** | | C **** |0. 1.| C **** C **** G = |-ALFA1| C **** | | C **** |-ALFA2| C **** C **** WHERE GLOBAL MEAN DENSITY LBCS: MASS MIXING C **** Titan values for CH4 and H2 from Yung et al (1984) C **** MIXCH4 = 1.7%; MIXH2 = 0.1% C **** ALFA1 = 9.767E-03 PHICH4 C **** ALFA2 = 7.100E-05 PHIH2 C **** DO 1 I = 1,nlonp4 EE(I,1,1) = 0. EE(I,1,2) = 0. EE(I,2,1) = 0. EE(I,2,2) = 0. FF(I,1,1) = 1. FF(I,1,2) = 0. FF(I,2,1) = 0. FF(I,2,2) = 1. GG(I,1) = -ALFA1 GG(I,2) = -ALFA2 1 CONTINUE RETURN end subroutine bndef end module bndry_module