#include "dims.h" ! ! am_02/02: include general wave (origin Kelvin wave) ! module bndry_module ! ! Lower boundary conditions: diurnal, semi-diurnal, general wave, annual t. ! use cons_module,only: jmax,t0,pi,atm_amu,gask,grav,freq_semidi, | re,dlamda,imaxp4,tbound,tgrad,kmax,kmaxp1,len1,cs,cor,tn, | freq_ann implicit none #include "params.h" ! complex,parameter :: ci=(0.,1.) complex :: | zb(zjmx),zb2(zjmx),zba(zjmx),zbw(zjmx), | tb(zjmx),tb2(zjmx),tba(zjmx),tbw(zjmx), | ub(zjmx),ub2(zjmx),uba(zjmx),ubw(zjmx), | vb(zjmx),vb2(zjmx),vba(zjmx),vbw(zjmx), | bnd(zimxp),bnd2(zimxp),bnda(zimxp),bndw(zimxp) ! contains !----------------------------------------------------------------------- subroutine lowbound ! ! Call the lower boundary routines. This is called once per model run ! from main tgcm.F. ! real :: ttbound,nob,avto,hor COMMON/RUNTIM/TTBOUND,NOB(ZJMX),AVTO,HOR(ZJMX) ! ! Local: integer :: k C **** C **** VALUE OF T AT LOWER BOUNDARY AND GRADIENT C **** ! TBOUND = TTBOUND ! TGRAD = 6. T0(1) = TBOUND T0(2) = TBOUND+dz*TGRAD call bndry_semidiurnal call bndry_diurnal call bndry_wave call bndry_annual ! C **** SET T0=0. do k=1,kmax T0(K)=0. enddo t0(kmaxp1)=0. 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. ! ! am 2001-8-29: modes are (2,2),(2,3),(2,4),(2,5),(2,6) ! same order for the amplitude [cm] and phases [hr] ! ! use input_module,only: tide ! ! Local: integer,parameter :: nalf=19, malf=2 real :: p(zjmx,nalf,malf),hough(zjmx,5,malf),cp(nalf/2+1) complex :: dzb(zjmx) real :: B(5,19),RL(5),BHOUR(5),rlamda,xdot(19),ydot(19), | ptscal,theta,ptjm(2*zjmx+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*JMAX+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,jmax p(j,nm1,mm1) = ptjm(2*(jmax+1-j)) enddo enddo DO J=1,JMAX P(J,nm1,2)=SQRT(FLOAT(n*(n+1)-6))*P(J,nm1,2)-2.*TN(J)* | P(J,nm1,1) enddo enddo ! ! Original code, using old alfpac (was alfpac.f): ! DO L=2,20 ! LM1=L-1 ! DO M=2,3 ! MM1=M-1 ! CALL LFK(L,M,CP,W) ! CALL LFP(L,M,JM,CP,PT) ! DO J=1,JMAX ! P(J,LM1,MM1) = PT(2*(JMAX+1-J)) ! enddo ! enddo ! DO J=1,JMAX ! P(J,LM1,2)=SQRT(FLOAT(L*(L+1)-6))*P(J,LM1,2)-2.*TN(J)* ! | P(J,LM1,1) ! enddo ! enddo ! ! util.F: real function sddot(n,x,y) DO L=1,5 DO LD=1,2 DO J=1,JMAX 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,JMAX TB(J)=0. ZB(J)=0. DZB(J)=0. 5 CONTINUE DO L=1,5 DO J=1,JMAX 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,JMAX 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,JMAX 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,IMAXP4 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. ! use input_module,only: tide2 ! ! Local: integer,parameter :: nalf=19, malf=2 real :: p(zjmx,nalf,malf),hough(zjmx,5,malf),cp(nalf/2+1) complex :: dzb(zjmx) real :: B(1,19),RL (1),BHOUR2(1),rlamda,xdot(19),ydot(19), | ptscal,theta,ptjm(2*zjmx+1) 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 = tide2(2) if (all(tide2==0.)) goto 13 DO 1 N=1,1 ZEE(N)=TIDE2(N)*CEXP(CI*pi*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*JMAX+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,jmax p(j,n,m) = ptjm(2*(jmax+1-j)) enddo enddo DO J=1,JMAX P(J,n,2)=SQRT(FLOAT(n*(n+1)-2))*P(J,n,2)-TN(J)*P(J,n,1) enddo enddo ! ! Original code, using old alfpac (was alfpac.f): ! DO L=1,19 ! DO M=1,2 ! CALL LFK(L,M,CP,W) ! CALL LFP(L,M,JM,CP,PT) ! DO J=1,JMAX ! P(J,L,M)=PT(2*(JMAX+1-J)) ! enddo ! enddo ! DO J=1,JMAX ! P(J,L,2)=SQRT(FLOAT(L*(L+1)-2))*P(J,L,2)-TN(J)*P(J,L,1) ! enddo ! enddo ! ! util.F: real function sddot(n,x,y) DO L=1,1 DO LD=1,2 DO J=1,JMAX 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,JMAX TB2(J)=0. ZB2(J)=0. DZB(J)=0. 5 CONTINUE DO 6 L=1,1 DO 6 J=1,JMAX 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,JMAX 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,JMAX 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,IMAXP4 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. ! C **** TIDAL BOUNDARY CONDITION FOR ANNUAL MODE ! use input_module,only: tideann use init_module,only: iday ! include "strt.h" ! for iday ! ! For 1998 spherepack lib code (sphpac.f) ! (replaces old alfpac.f) ! integer,parameter :: nalf=24, malf=2 ! ! Local: real :: p(zjmx,nalf,malf),hough(zjmx,0:6,2),cp(nalf/2+1) complex :: dzb(zjmx),zzb(zjmx) real :: B(6,24),RL(0:6),scale,rt2,rm,factor,xdot(24),ydot(24), | pi,ptscal,ptjm(2*zjmx+1),theta,w(zjmx) integer :: n,l,m,i,k,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*JMAX+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,jmax p(j,n,mp1) = ptjm(2*(jmax+1-j)) enddo enddo m = 0 rm = float(m) do j=1,jmax 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 ! ! Original code, using old alfpac (was alfpac.f): ! DO 2 L=1,24 ! DO 3 M=0,1 ! MP1=M+1 ! CALL LFK(L,M,CP,W) ! CALL LFP(L,M,JM,CP,PT) ! DO 3 J=1,JMAX ! P(J,L,MP1) = PT(2*(JMAX+1-J)) ! 3 CONTINUE ! M = 0 ! RM = FLOAT(M) ! DO 2 J=1,JMAX ! P(J,L,2)=SQRT(FLOAT(L*(L+1)-M*(M+1)))*P(J,L,2)-RM*TN(J)* ! 1 P(J,L,1) ! 2 CONTINUE C **** C **** NOW EVALUATE HOUGH FUNCTIONS C **** DO 4 L=1,6 DO 4 LD=1,2 DO 4 J=1,JMAX 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,JMAX 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,JMAX 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,JMAX 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,JMAX 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,IMAXP4 BNDA(I) = 1. 8 CONTINUE end subroutine bndry_annual !----------------------------------------------------------------------- SUBROUTINE BNDCMP C **** C **** CALCULATE MATRICES B(ZIMXP,2,2) AND VECTORS FB(ZIMXP,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 **** ! include "vscr.h" #include "cmpbnd.h" ! ! Local: real :: EE(ZIMXP,2,2),FF(ZIMXP,2,2),GG(ZIMXP,2),WM1(ZIMXP,2,2), 1 WM2(ZIMXP,2,2),WM3(ZIMXP,2,2),WV1(ZIMXP,2),WS1(ZIMXP) integer :: l,m,i C **** C **** CALL BNDEF TO DEFINE E, F AND G IN S1, S2 AND S3 C **** C **** CALL BNDEF(ee,ff,gg) C **** C **** WM1 = (E/DS - F/2.) C **** C **** WM1 = (E/DS + F/2.) C **** DO 1 L = 1,2 DO 1 M = 1,2 DO 1 I = 1,LEN1 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. 1 CONTINUE C **** C **** WM3 = WM1**(-1) C **** C **** WS1 = DET(WM1) C **** DO 2 I = 1,LEN1 WS1(I) = WM1(I,1,1)*WM1(I,2,2)-WM1(I,1,2)*WM1(I,2,1) 2 CONTINUE C **** C **** NOW INVERSE OF WM1 IN WM3 C **** DO 3 I = 1,LEN1 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) 3 CONTINUE C **** C **** B = WM3 * WM2 C **** DO 4 L = 1,2 DO 4 M = 1,2 DO 4 I = 1,LEN1 B(I,L,M) = WM3(I,L,1)*WM2(I,1,M)+WM3(I,L,2)*WM2(I,2,M) 4 CONTINUE C **** C **** FB = WM3 * G C **** DO 5 L = 1,2 DO 5 I = 1,LEN1 FB(I,L) = WM3(I,L,1)*GG(I,1)+WM3(I,L,2)*GG(I,2) 5 CONTINUE RETURN 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(ZIMXP,2,2),FF(ZIMXP,2,2),GG(ZIMXP,2) ! ! Local: real :: alfa 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 **** DATA ALFA/0.234/ DO 1 I = 1,LEN1 EE(I,1,1) = 0. EE(I,1,2) = 0. EE(I,2,1) = 0. EE(I,2,2) = 1. FF(I,1,1) = 1. FF(I,1,2) = 1. FF(I,2,1) = 0. FF(I,2,2) = -1. GG(I,1) = -ALFA GG(I,2) = 0. 1 CONTINUE RETURN end subroutine bndef !----------------------------------------------------------------------- subroutine bndry_wave ! ! general wave ! based on Chapman & Lindzen: Atmospheric Tides pp.109 ! Richmond: Atmospheric Tides pp. 467 ! tested with diurnal,semidiurnal & Kelvin waves ! difference equation to code: ! - colatitude -> latitude in code (inner derivative -1) ! - southward velocity -> northward velocity (-1*vbw) ! ! input: wtyp: typ of wave: 1 = Kelvin wave ! wper: period of wave (days) ! wper < 0 -> east, wper > 0 -> west ! wamp: amplitude in cm ! wphase: phase in hrs ! wheqn: equivalent depth in cm ! wzno: zonal wave number (wzno > 0) ! ! output: zbw,tbw,ubw,vbw ! use input_module,only: wave #include "consdyn.h" ! include "consts.h" ! ! Local: integer,parameter :: nalf=19, malf=2 real :: p(zjmx,nalf,malf),hough(zjmx,5,malf),cp(nalf/2+1), | expfac,dexpfac,fac,wper,wamp,wphase,wheqn,wzno, | c23,c23w,fd, | B(1,19),rlamda,xdot(19),ydot(19),pi, | ptscal,theta,ptjm(2*zjmx+1) complex :: dzb(zjmx) integer :: l,m,j,n,jm,ld,i,wtyp,nm1,mm1 complex zee,cl,expdlm ! data B/ ! B(19) for the diurnal tide: 2 0.282710, ! B must be changed for other tides 3 0.0 , ! and caculation of Hough function 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/ ! real,external :: sddot ! in util.F ! ! input data wtyp = int(wave(1)) ! typ of wave: 1 = Kelvin wave wper = wave(2) ! period of wave (days) ! wper < 0 -> east, wper > 0 -> west wamp = wave(3) ! amplitude in cm wphase = wave(4) ! phase in hrs wheqn = wave(5) ! equivalent depth in cm wzno = wave(6) ! zonal wave number (wzno > 0) ! ! ... no wave if (all(wave==0.)) goto 13 ! ... kelvin wave if(wtyp.eq.1) then zee = cmplx(wamp) fac = (gask/(atm_amu*grav*wheqn)* * (T0(1)*2./7.+(T0(2)-T0(1))/dz)-.25) cl = -csqrt(cmplx(fac))-.5*CI ! ! ... instead of hough functions ! set up sinus function and derivative ! from dissipationless momentum equations without background mean winds ! hough = cos(phi)^(-2*s*Omega/omega) Omega angular frequency of earth ! omega angular frequency of 3 days ! s is the zonal wave number ! h=cos(phi)**6.01643 and ! h'=-1*6.01643*cos(phi)**5.01643*sin(phi) ! phi: latitude -> inner derivative ! expfac = 6.01643 dexpfac = expfac -1.0 ! do 4 j=1,jmax hough(j,1,1) = cos(YLATG(j))**expfac hough(j,1,2) = -expfac*cos(YLATG(j))**dexpfac*sin(YLATG(j)) 4 continue ! ... no kelvin wave elseif(wtyp.gt.1) then fac = wper/abs(wper) ! travel direction zee = wamp*cexp(CI*2.*pi*wphase/24.0*wzno) cl = fac*csqrt(cmplx(gask/(atm_amu*grav*wheqn)* | (T0(1)*2./7.+(T0(2)-T0(1))/dz)-.25))-.5*CI jm=2*jmax+1 ! ! ... hough functions for diurnal tides ! Using new (1998) spherepack (sphpac.f): pi = 4.*atan(1.) 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,jmax p(j,n,m) = ptjm(2*(jmax+1-j)) enddo enddo DO J=1,JMAX 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,JMAX xdot(:) = p(j,:,ld) ydot(:) = b(l,:) HOUGH(J,L,LD)=sddot(19,xdot,ydot) enddo enddo enddo ! endif ! ... caculate tbw,zbw,vbw,ubw ! do 6 j=1,jmax zbw(j) = zee*hough(j,1,1) dzb(j) = -1*zee*hough(j,1,2) ! -1 inner derivative colat->lat tbw(j) = CI*atm_amu*grav/gask*zee*cl*hough(j,1,1) 6 continue ! c23 = 0.5*freq_semidi ! c23 = 2Pi/(24*60*60) c23w = c23/wper ! wave frequency ! ... kelvin wave if(wtyp.eq.1) then ! factor for approximating the amplitude for ubw ! from solving momentum eqn u = -g/(sigma*R)*h ! sigma angular frequency of the Kelvin wave fac = -1.0*grav/(re*C23w) do 17 j=1,jmax fd = re*(c23w*c23w-COR(j)*COR(j)) vbw(j) = -CI*grav*(c23w*dzb(j)+wzno*COR(j)/CS(j)*zbw(j))/fd ubw(j) = fac*zbw(j) 17 continue ! ... no kelvin waves elseif(wtyp.gt.1) then do 7 j=1,jmax ubw(j) = re*(c23w*c23w-COR(j)*COR(j)) vbw(j) = -CI*grav*(c23w*dzb(j)+wzno*COR(j)/ | CS(j)*zbw(j))/ubw(j) ubw(j) = -grav*(COR(J)*dzb(J)+c23w*wzno/ | CS(J)*zbw(j))/ubw(j) 7 continue endif ! go to 11 13 continue ! ... no waves defined wave(2) = 1.0 ! period of wave to avoid numerical problems in duv.f and dt.f do 12 j=1,jmax zbw(j)=0. tbw(j)=0. vbw(j)=0. ubw(j)=0. 12 continue ! 11 continue ! ! ... calculate longitudial structure rlamda = -2.*dlamda bndw(1)= cexp(CI*wzno*rlamda) expdlm = cexp(CI*wzno*dlamda) ! do 9 i=2,imaxp4 bndw(i)=bndw(i-1)*expdlm 9 continue ! end subroutine bndry_wave end module bndry_module