#include "dims.h" ! module bndry_module ! ! Lower boundary conditions diurnal, semi-diurnal, annual, composition, ! 2-day wave, and planetary. ! implicit none #include "params.h" ! complex,parameter :: ci=(0.,1.) complex :: | zb(zjmx),zb2(zjmx),zba(zjmx), | tb(zjmx),tb2(zjmx),tba(zjmx), | ub(zjmx),ub2(zjmx),uba(zjmx), | vb(zjmx),vb2(zjmx),vba(zjmx), | bnd(zimxp),bnd2(zimxp),bnda(zimxp) ! contains !----------------------------------------------------------------------- subroutine lowbound use init_module,only: iday,iter use input_module,only: step,ncep,nmc use cons_module,only: tbound,dz,t0,tgrad #include "crates_const.h" real :: dday ! t0(1) = tbound t0(2) = tbound+dz*tgrad if (ncep <= 0 .and. nmc <= 0) call zatmos_import(iday) ! ! Call the lower boundary routines. This is called once per model run ! from main tgcm.F. ! call bndry_semidiurnal call bndry_diurnal dday = float(iday)+amod(float(iter)*step,86400.)/86400. call bndry_annual(dday) call cco2gr(co2mix) call bndry_2day ! t0(:) = 0. ! end subroutine lowbound !----------------------------------------------------------------------- subroutine bndry_semidiurnal use input_module,only: tide use cons_module,only: pi,jmax,imaxp4,cs,cor,tn,t0,atm_amu,gask,re, | grav,freq_semidi,dz,dlamda ! ! Local: integer,parameter :: nalf=19, malf=2 real :: p(zjmx,nalf,malf),hough(zjmx,5,malf), | cp(nalf/2+1),w(60),xdot(nalf),ydot(nalf),ptscal,theta, | ptjm(2*zjmx+1) complex :: dzb(zjmx) COMPLEX ZEE(5),CL(5),EXPDLM real :: B(5,nalf) = reshape(source = | (/0.969152, 0.0 , 0.216046, 0.0 , 0.093838, | 0.0 , 0.909763, 0.0 , 0.342113, 0.0 , | -0.245226, 0.0 , 0.798445, 0.0 , 0.421218, | 0.0 ,-0.408934, 0.0 , 0.645517, 0.0 , | 0.024633, 0.0 ,-0.543993, 0.0 , 0.464159, | 0.0 , 0.071127, 0.0 ,-0.643189, 0.0 , | -0.001292, 0.0 , 0.139613, 0.0 ,-0.699495, | 0.0 ,-0.006673, 0.0 , 0.225090, 0.0 , | 0.000042, 0.0 ,-0.019654, 0.0 , 0.320141, | 0.0 , 0.000394, 0.0 ,-0.043345, 0.0 , | -0.000001, 0.0 , 0.001772, 0.0 ,-0.079831, | 0.0 ,-0.000016, 0.0 , 0.005401, 0.0 , | 0.0 , 0.0 ,-0.000112, 0.0 , 0.012932, | 0.0 , 0.0 , 0.0 ,-0.000476, 0.0 , | 0.0 , 0.0 , 0.000005, 0.0 ,-0.001490, | 0.0 , 0.0 , 0.0 , 0.000031, 0.0 , | 0.0 , 0.0 , 0.0 , 0.0 , 0.000129, | 0.0 , 0.0 , 0.0 ,-0.000002, 0.0 , | 0.0 , 0.0 , 0.0 , 0.0 ,-0.000009/), | shape = (/5,19/)) real :: RL(5) = | (/7.8519E5, 3.6665E5, 2.1098E5, 1.3671E5, 0.9565E5/) real :: BHOUR(5),rlamda integer :: i,n,jm,l,m,j,ld,lm1,mm1,nm1 real,external :: sddot ! 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 ! **** 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 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) ! HOUGH(J,L,LD)=sddot(19,P(J,1,LD),JMAX,B(L,1),5) enddo enddo enddo C **** GENERATE ZB, UB, VB AND TB DO J=1,JMAX TB(J)=0. ZB(J)=0. DZB(J)=0. enddo 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 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) enddo GO TO 11 13 CONTINUE C **** ZERO BOUNDARY CONDITION DO J=1,JMAX ZB(J)=0. TB(J)=0. UB(J)=0. VB(J)=0. enddo 11 CONTINUE C **** CALCULATE LONGITUDINAL STRUCTURE RLAMDA = -2.*dlamda BND(1)=CEXP(CI*2.*RLAMDA) EXPDLM=CEXP(CI*2.*dlamda) DO I=2,IMAXP4 BND(I)=BND(I-1)*EXPDLM enddo RETURN end subroutine bndry_semidiurnal !----------------------------------------------------------------------- subroutine bndry_diurnal use input_module,only: tide2 use cons_module,only: pi,jmax,imaxp4,cs,cor,tn,t0,dlamda,dz,grav, | gask,freq_semidi,re,atm_amu C **** C **** TIDAL BOUNDARY CONDITION FOR DIURNAL GRAVITATIONAL MODE C **** (1,1) C **** integer,parameter :: nalf=19, malf=2 real :: p(zjmx,nalf,malf),hough(zjmx,5,malf), | cp(nalf/2+1),w(60),xdot(nalf),ydot(nalf),ptscal,theta, | ptjm(2*zjmx+1) complex :: dzb(zjmx) COMPLEX ZEE(1),CL(1),EXPDLM real :: RL(1) = (/0.6909E5/) real :: BHOUR2(1),rlamda integer :: n,jm,l,m,j,ld,i real :: B(1,19) = reshape(source = | (/0.282710, | 0.0 , | -0.638229, | 0.0 , | 0.620521, | 0.0 , | -0.336408, | 0.0 , | 0.117021, | 0.0 , | -0.028332, | 0.0 , | 0.005042, | 0.0 , | -0.000686, | 0.0 , | 0.000074, | 0.0 , | -0.000006/), | shape=(/1,19/)) real,external :: sddot ! bhour2 = tide2(2) if (tide2(1)==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))* 1(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 2 L=1,19 ! DO 3 M=1,2 ! CALL LFK(L,M,CP,W) ! CALL LFP(L,M,JM,CP,PT) ! DO 3 J=1,JMAX ! P(J,L,M)=PT(2*(JMAX+1-J)) ! 3 CONTINUE ! DO 2 J=1,JMAX ! P(J,L,2)=SQRT(FLOAT(L*(L+1)-2))*P(J,L,2)-TN(J)*P(J,L,1) ! 2 CONTINUE DO 4 L=1,1 DO 4 LD=1,2 DO 4 J=1,JMAX xdot(:) = p(j,:,ld) ydot(:) = b(l,:) HOUGH(J,L,LD)=sddot(19,xdot,ydot) ! HOUGH(J,L,LD)=SDOT(19,P(J,1,LD),JMAX,B(L,1),1) 4 CONTINUE 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 RETURN end subroutine bndry_diurnal !----------------------------------------------------------------------- subroutine bndry_annual(dday) use input_module,only: tideann use cons_module,only: pi,jmax,imaxp4,cs,cor,tn,t0,re,grav,gask, | freq_ann,dz,atm_amu C **** C **** TIDAL BOUNDARY CONDITION FOR ANNUAL MODE ! ! 8/20/98: this routine called initially from lowbound (via tgcm.F), ! then from advance at every time step only if calendar day ! is being advanced. Dday is current decimal day of year. ! ! Args: real,intent(in) :: dday ! integer,intent(in) :: idayin ! ! For 1998 spherepack lib code (sphpac.f) ! (replaces old alfpac.f) ! integer,parameter :: nalf=24, malf=2 ! ! Local: real,parameter :: SCALE=1. integer :: n,l,m,mp1,j,jm,ld,i,k real :: rm,rt2,factor 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),xdot(nalf),ydot(nalf),ptscal, | ptjm(2*zjmx+1),theta,w(zjmx) COMPLEX CC(0:6,0:6),CL(0:6),EXPT real,external :: sddot ! 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 = -17. C **** DATA ((CC(K,N),N=0,6),K=0,3)/ 1 ( 0.421595E+02, 0.000000E+00),(-0.226810E+00, 0.000000E+00), 2 ( 0.415692E+00, 0.000000E+00),( 0.547413E-01, 0.000000E+00), 3 (-0.312575E-01, 0.000000E+00),( 0.280035E-01, 0.000000E+00), 4 (-0.185304E-01, 0.000000E+00),( 0.300941E-01,-0.734916E-03), 5 ( 0.100061E+01,-0.168070E+00),(-0.212277E+00, 0.133133E+00), 6 ( 0.226725E-01,-0.553615E-01),(-0.426632E-01, 0.941169E-03), 7 (-0.109421E-02,-0.877230E-02),(-0.576746E-02,-0.113662E-01), 8 ( 0.215525E-01, 0.269012E-01),( 0.212126E-01, 0.377161E-01), 9 (-0.329359E-01, 0.112051E-01),(-0.117842E-01, 0.145664E-01), * (-0.813882E-02, 0.424029E-01),(-0.705682E-02, 0.100935E-01), 1 (-0.589241E-02, 0.224780E-01),(-0.956909E-03, 0.169400E-02), 2 ( 0.283082E-02,-0.579849E-02),( 0.141201E-02,-0.351546E-02), 3 ( 0.288114E-02, 0.600793E-04),(-0.896708E-04, 0.917448E-04), 4 ( 0.146016E-02, 0.152191E-02),( 0.911123E-04,-0.141497E-03)/ DATA ((CC(K,N),N=0,6),K=4,6)/ 1 (-2.200025E-02,-0.159748E-02),( 0.122976E-03,-0.180022E-03), 2 ( 0.218866E-02, 0.217849E-02),( 0.361138E-03, 0.264706E-03), 3 ( 0.803136E-04, 0.595681E-04),(-0.711489E-04, 0.138662E-04), 4 (-0.205981E-04, 0.237149E-05),( 0.211351E-04, 0.146409E-04), 5 (-0.195730E-02,-0.810895E-03),(-0.442840E-04, 0.292211E-04), 6 (-0.123450E-04,-0.538276E-04),(-0.253900E-04,-0.621475E-04), 7 ( 0.767875E-04,-0.422957E-05),( 0.206095E-04, 0.178031E-05), 8 (-0.191415E-03, 0.722780E-04),(-0.275580E-04,-0.427594E-04), 9 ( 0.132229E-03,-0.182378E-03),( 0.341214E-04, 0.355103E-05), * ( 0.326348E-04, 0.543373E-04),(-0.109405E-04, 0.189210E-04), 1 ( 0.489606E-04, 0.396533E-05)/ ! ! 8/20/98: This routine now called at every time step from advnce if ! the run is advancing calendar day. The single argument is current ! decimal model day. ! 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) ! 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)*(dday-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 RETURN end subroutine bndry_annual !----------------------------------------------------------------------- subroutine bndry_comp use cons_module,only: len1 C **** C **** CALCULATE MATRICES B(ZIMXP,2,2) AND VECTORS FB(ZIMXP,2) C **** REPRESENTING BOUNDARY CONDITION IN COMP, WHERE PSI1 C **** AND PSI2 ARE EVALUATED. C **** #include "fcom.h" #include "vscr.h" #include "index.h" #include "buff.h" #include "cmpbnd.h" #include "lowbnd.h" #include "phys.h" #include "mwt.h" C **** C **** CURRENT BOUNDARY CONDITION IS: C **** C **** 0.5*(PSI1(-1/2)+PSI1(1/2)) = PSO2LB C **** 0.5*(PSI2(-1/2)+PSI2(1/2)) = XOXLB*RMOX(REAL)/MBAR C **** WHERE: C **** XOXLB = OX NUMBER DENSITY MIXING RATIO AT LOWER C **** BOUNDARY C **** PSO2LB= O2 MASS MIXING RATIO AT LOWER BOUNDARY C **** C **** THIS GIVES: C **** PSI1(-1/2) = B(1,1)*PSI1(1/2)+B(1,2)*PSI2(1/2)+FB(1) C **** PSI2(-1/2) = B(2,1)*PSI1(1/2)+B(2,2)*PSI2(1/2)+FB(2) C **** WHERE: C **** B(1,1) = B(2,2) = -1. C **** B(1,2) = B(2,1) = 0. C **** FB(1) = 2.*PSO2LB C **** FB(2) = 2.*XOXLB*RMOX(TRUE)/MBAR C **** ! ! Local: ! PSO2LB adjusted from 0.22 to 0.24, as per modsrc.kibo: real,parameter :: PSO2LB=0.24 integer :: i,k,kk,nmsk C **** C **** SET MATRIX B C **** DO 1 I = 1,LEN1 B(I,1,1) = -1. B(I,2,2) = -1. B(I,1,2) = 0. B(I,2,1) = 0. 1 CONTINUE C **** C **** T1, T2 = RMOX(TRUE), K = 1/2, 3/2 C **** DO K = 1,2 KK = (K-1)*LEN1 DO I = 1,LEN1 T1(I+KK) = RMO3 enddo enddo C **** C **** T1 = RMOX(TRUE), K=0 (EXTRAPOLATION) C **** SET FB(1) AND FB(2) C **** NMSK = NJ+NMS DO 3 I = 1,LEN1 T1(I) = 1.5*T1(I)-0.5*T2(I) FB(I,1) = 2.*PSO2LB FB(I,2) = 2.*XOXLB(J)*T1(I)/F(I,NMSK) 3 CONTINUE RETURN end subroutine bndry_comp !----------------------------------------------------------------------- subroutine bndry_2day use input_module,only: tide3m3 use cons_module,only: pi,jmax,imaxp4,cs,cor,tn,t0,re,dz,dlamda, | gask,freq_3m3,grav,atm_amu C **** C **** Tidal boundary condition for the two day mixed C **** Rossby-Gravity wave. This is denoted by indices (3,-3) C **** and is asymmetric. The charactistic depth is 10.5km. C **** See Forbes, Jeffrey M. Tidal and Planetary Waves. In: C **** Upper Mesosphere and Lower Thermosphere: A Review of C **** Experiment and Theory. Geophysical Monograph 87, C **** Copyright 1995 by American Geophysical Union. C **** C **** The following common blocks are needed: C **** #include "cons3m3.h" C **** C **** Local variables C **** integer,parameter :: nalf=19, malf=2 real :: p(zjmx,nalf,malf),hough(zjmx,5,malf), | cp(nalf/2+1),xdot(10),ydot(10),ptjm(2*zjmx+1),ptscal, | theta,w(zjmx) complex :: dzb(zjmx) COMPLEX ZEE,CL,EXPDLM C **** C **** Where: C **** C **** B(10) = Coefficients of normalised associated Legendre C **** polynomials needed to build Hough(3,-3). C **** i.e.: P(l,3), where l = 4,6,8,10,12 (higher C **** terms are too small to merit consideration) C **** C **** RL = equivalent depth for wave C **** = 10.5km = 0.105E7 cm C **** C **** BHOUR = phase of wave measured in hours with constraint: C **** 0 .le. BHOUR .lt. p C **** ( p = period of wave in hours = 2.07412*24) C **** C **** TIDE = amplitude of wave in cm C **** C **** ZEE = complex amplitude incorporating phase (BHOUR) and C **** amplitude (TIDE) C **** C **** EXPDLM = exp(2*pi*i*dlamda), where dlamda is C **** longitudinal grid spacing (radians) C **** real :: B(10) = (/ | 0.0 , | -0.997138, | 0.0 , | 0.075553, | 0.0 , | -0.002905, | 0.0 , | 0.000069, | 0.0 , | -0.000001/) real,parameter :: RL=0.105E7 real :: bhour,tide,rlamda integer :: jm,n,nm2,m,mm2,j,ld,i real,external :: sddot integer :: itide3m3 ! tide = tide3m3(1) bhour = tide3m3(2) itide3m3 = 0 if (tide3m3(1) /= 0..or.tide3m3(2) /= 0.) itide3m3 = 1 C **** C **** Check presence of two-day wave C **** C******************************************************** C WRITE (6,*) ITIDE3M3, TIDE, BHOUR, IPROC C******************************************************** IF(ITIDE3M3.EQ.1)THEN C **** C **** Two-day tide is present C **** C **** ZEE = complex amplitude of wave at time zero. C **** = amp * exp(2*pi*i*BHOUR/P2DAY) C **** C **** CL = vertical wave length C **** = ((K*H + dH/ds)/hn - 1/4)**(1/2) -i/2 C **** C **** where: C **** K = 2/7 C **** H = scale height (cm) C **** hn = equivalent depth for mode = 1.05E7 cm C **** s = vertical pressure coordinate C **** ZEE = TIDE*CEXP(CI*freq_3m3*BHOUR*60.*60.) CL = CSQRT(CMPLX(gask/(atm_amu*grav*RL)* 1 (T0(1)*2./7.+(T0(2)-T0(1))/dz)-.25))-.5*CI C*********************************** C WRITE(6,*)'CL = ', CL C*********************************** JM=2*JMAX+1 C **** C **** Set up Hough functions C **** ! ! Using new (1998) spherepack (sphpac.f): do n=3,12 nm2 = n-2 do m=3,4 mm2=m-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,nm2,mm2) = ptjm(2*(jmax+1-j)) enddo enddo m = 3 DO J=1,JMAX P(J,nm2,2)=SQRT(FLOAT(n*(n+1)-m*(m+1)))*P(J,nm2,2)- | float(m)*tn(j)*P(J,nm2,1) enddo enddo ! ! DO L = 3,12 ! LM2 = L-2 ! DO M = 3,4 ! MM2 = M-2 ! CALL LFK(L,M,CP,W) ! CALL LFP(L,M,JM,CP,PT) ! DO J = 1,JMAX ! P(J,LM2,MM2) = PT(2*(JMAX+1-J)) ! ENDDO ! ENDDO ! M = 3 ! DO J = 1,JMAX ! P(J,LM2,2) = SQRT(FLOAT(L*(L+1)-M*(M+1)))*P(J,LM2,2)- ! 1 FLOAT(M)*TN(J)*P(J,LM2,1) ! ENDDO ! ENDDO ! DO LD = 1,2 DO J = 1,JMAX xdot(:) = p(j,1:10,ld) ydot(:) = b(:) HOUGH(J,1,LD) = sddot(10,xdot,ydot) ! HOUGH(J,1,LD) = sddot(10,P(J,1,LD),JMAX,B,1) ENDDO ENDDO C **** C **** Generate ZB3M3, UB3M3, VB3M3 and TB3M3 C **** DO J = 1,JMAX ZB3M3(J) = ZEE*HOUGH(J,1,1) DZB(J) = ZEE*HOUGH(J,1,2) TB3M3(J) = CI*atm_amu*grav/gask*ZEE*CL*HOUGH(J,1,1) ENDDO DO J = 1,JMAX UB3M3(J) = freq_3m3*re*(1.-(COR(J)/freq_3m3)**2) VB3M3(J) = CI*grav*(DZB(J)-FLOAT(M)*COR(J)/(freq_3m3*CS(J)) 1 *ZB3M3(J))/UB3M3(J) UB3M3(J) = grav*(COR(J)/freq_3m3*DZB(J)-FLOAT(M)/CS(J)* 1 ZB3M3(J))/UB3M3(J) ENDDO C **** C **** Calculate longitudinal structure C **** RLAMDA = -2.*dlamda BND3M3(1)=CEXP(CI*FLOAT(M)*RLAMDA) EXPDLM=CEXP(CI*FLOAT(M)*dlamda) DO I=2,IMAXP4 BND3M3(I)=BND3M3(I-1)*EXPDLM ENDDO ! IF(ITIDE3M3.EQ.1.AND.IPROC.LE.0)CALL PLT3M3 ELSE C **** C **** Zero boundary condition, (3,-3) wave is absent. C **** C **** ZERO BOUNDARY CONDITION DO J=1,JMAX ZB3M3(J)=0. TB3M3(J)=0. UB3M3(J)=0. VB3M3(J)=0. ENDDO ENDIF RETURN end subroutine bndry_2day !----------------------------------------------------------------------- subroutine bndry_planetary use init_module,only: iter use cons_module,only: len1,jmax,pi,dlamda,dt #include "bndz.h" C **** C **** Calculate contribution to ZB from planetary waves C **** real,parameter :: amplan=5.0E4, toplan=2.5E+5 ! DATA AMPLAN,T0PLAN/ 5.0E4,2.5E+5/ C DATA AMPLAN,T0PLAN/ 11.0E4,2.5E+5/ ! ! Local: real :: time integer :: j,i ! TIME = ITER*dt DO J = 1,24 DO I = 1,LEN1 ZBND(I,J) = 0. ENDDO ENDDO DO J = 25,35 DO I = 1,LEN1 C ZBND(I,J) = AMPLAN*(SIN((-87.5+(J-1)*5.-30.)*pi/(60.- C 1 (-87.5+(J-1)*5.-60.)/3.)))**2* C 2 (1.-EXP(-TIME/T0PLAN))*SIN(1.*(-pi+(I-3)*dlamda)) C ZBND(I,J) = AMPLAN*(SIN((-87.5+(J-1)*5.-30.)*pi/(60.- C 1 (-87.5+(J-1)*5.-60.)/3.)))**2* C 2 SIN(1.*(-pi+(I-3)*dlamda)) ZBND(I,J) = 0. ENDDO ENDDO J = JMAX DO I = 1,LEN1 ZBND(I,J) = 0. ENDDO RETURN end subroutine bndry_planetary end module bndry_module