! 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) ! 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,zibot ! ! GSWM boundary perturbation (see gswm.F): use gswm_module,only: | gswm_mi_di_z, gswm_mi_sdi_z, gswm_nm_di_z, gswm_nm_sdi_z, | gswm_mi_di_t, gswm_mi_sdi_t, gswm_nm_di_t, gswm_nm_sdi_t, | gswm_mi_di_u, gswm_mi_sdi_u, gswm_nm_di_u, gswm_nm_sdi_u, | gswm_mi_di_v, gswm_mi_sdi_v, gswm_nm_di_v, gswm_nm_sdi_v 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) ! for bndcmp (old) real :: bnew(nlonp4,2,2),fbnew(nlonp4,2,nlat) ! for bndry_comp (new) ! contains !----------------------------------------------------------------------- subroutine lowbound ! use init_module,only: igswm_mi_di,igswm_mi_sdi, ! GSWM input flags | iday,iter ! for dday calculation use input_module,only: step ! Call the lower boundary routines. This is called once per model run ! from main tgcm.F. ! ! Local: integer :: k real :: dday ! ! 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(igswm_mi_sdi == 0) call bndry_semidiurnal ! check for GSWM semidiurnal if(igswm_mi_di == 0) call bndry_diurnal ! check for GSWM diurnal dday = float(iday)+amod(float(iter)*float(step),86400.)/86400. ! am 4/06 annual tide -> from MSIS ! call bndry_annual(dday) 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. ! 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. ! 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(dday) ! ! 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 ! ! am 1/06 add dday to advance the day for annual tides, only done ! if calendar day is also advanced ! Args: real,intent(in) :: dday ! ! 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)*(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,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 ! am 10/20/05 only temporary to get the geopotential height to approx. 80 km ! btf 1/17/06: Check zibot to determine zba bottom boundary: if (zibot==-7.) then ZBA(J) = CC(0,0)*1./RT2*1.E5 elseif (zibot==-10.) then ZBA(J) = 80./96.37*CC(0,0)*1./RT2*1.E5 else write(6,"('>>> bndry_annual: unknown zibot=',e12.4)") zibot call shutdown('bad zibot in bndry_annual') endif ! 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 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 bndry_comp(o2,o1,o3,no,o1d,barm,xnmbar,xolb, | lev0,lev1,lon0,lon1,lat) ! ! 6/16/06 btf: ! Calculate lower boundary matrices for ox and o2, similiar to timegcm. ! This is for new OX and partitioning to O3 and O in tiegcm-hist-lbc, ! and will replace sub bndcmp above. Note this is called per latitude ! from dynamics, before qrj. ! ! Also note that xoxlb is defined by solving quadratic (using o,o3,ox from ! previous timestep), rather than prescribing xoxlb from Solomon-Garcia, ! as in timegcm. ! ! This currently defines module data bnew and fbnew. When this is working, ! and the old bndcmp is removed, these can be renamed b and fb, as before. ! ! The following comments are from timegcm1.2 (sub bndry_comp in bndry.F): ! ! CALCULATE MATRICES B(nlonp4,2,2) AND VECTORS FB(nlonp4,2) ! REPRESENTING BOUNDARY CONDITION IN COMP, WHERE PSI1 ! AND PSI2 ARE EVALUATED. ! ! CURRENT BOUNDARY CONDITION IS: ! ! 0.5*(PSI1(-1/2)+PSI1(1/2)) = PSO2LB ! 0.5*(PSI2(-1/2)+PSI2(1/2)) = XOXLB*RMOX(REAL)/MBAR ! WHERE: ! XOXLB = OX NUMBER DENSITY MIXING RATIO AT LOWER BOUNDARY ! PSO2LB= O2 MASS MIXING RATIO AT LOWER BOUNDARY ! ! THIS GIVES: ! PSI1(-1/2) = B(1,1)*PSI1(1/2)+B(1,2)*PSI2(1/2)+FB(1) ! PSI2(-1/2) = B(2,1)*PSI1(1/2)+B(2,2)*PSI2(1/2)+FB(2) ! WHERE: ! B(1,1) = B(2,2) = -1. ! B(1,2) = B(2,1) = 0. ! FB(1) = 2.*PSO2LB ! FB(2) = 2.*XOXLB*RMOX(TRUE)/MBAR ! use cons_module,only: rmass_o3,rmassinv_o3,rmassinv_o1, | rmassinv_o2,rmassinv_no use addfld_module,only: addfld use chemrates_module,only: rkm21,rkm24,rkm7a,rkm7b,beta9,rk20 ! ! Args: integer,intent(in) :: | lev0,lev1, ! subdomain level indices | lon0,lon1, ! subdomain longitude indices | lat ! current latitude index real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) :: | o2, ! o2 | o1, ! o | o3, ! o3 | no, ! no | o1d, ! o1d | barm, ! mbar | xnmbar ! p0*e(-z)*barm/kT real,intent(in) :: xolb(nlonp4) ! ! Local: integer :: i,k real,parameter :: pso2lb=0.24 real :: rmo3a(nlonp4),rmo3b(nlonp4) real :: xoxlb real,dimension(lev0:lev1,lon0:lon1) :: | r6, ! ratio o3/o (cm3) | r7, ! ratio o/ox (cm3) | b2,b3,b4 ! quadratic terms ! ! Save inputs to secondary history: call addfld('BND_O3' ,' ',' ',o3(lev0:lev1-1,lon0:lon1), | 'lev',lev0,lev1-1,'lon',lon0,lon1,lat) call addfld('BND_O1' ,' ',' ',o1(lev0:lev1-1,lon0:lon1), | 'lev',lev0,lev1-1,'lon',lon0,lon1,lat) call addfld('BND_NO' ,' ',' ',no(lev0:lev1-1,lon0:lon1), | 'lev',lev0,lev1-1,'lon',lon0,lon1,lat) call addfld('BND_O1D',' ',' ',o1d(lev0:lev1-1,lon0:lon1), | 'lev',lev0,lev1-1,'lon',lon0,lon1,lat) call addfld('BND_BARM',' ',' ',barm(lev0:lev1-1,lon0:lon1), | 'lev',lev0,lev1-1,'lon',lon0,lon1,lat) call addfld('BND_XNMBAR',' ',' ',xnmbar(lev0:lev1-1,lon0:lon1), | 'lev',lev0,lev1-1,'lon',lon0,lon1,lat) ! ! Set matrix b: ! do i = 1,nlonp4 b(i,1,1) = -1. b(i,2,2) = -1. b(i,1,2) = 0. b(i,2,1) = 0. ! fb(i,1) = 2.*pso2lb fb(i,2) = 2.*xolb(i) enddo ! ! Define ratios: ! do i=lon0,lon1 ! do k=lev0,lev1-1 ! o3/o: ! r6(k,i) = xnmbar(k,i)*rkm21(k,i,lat)*o2(k,i)*rmassinv_o2 / ! | (xnmbar(k,i)* ! | (rkm24(k,i,lat)*o1 (k,i)*rmassinv_o1+ ! | (rkm7a+rkm7b) *o1d(k,i)*rmassinv_o1+ ! | beta9(k,i,lat) *no (k,i)*rmassinv_no)) ! o/ox: ! r7(k,i) = 1./(1.+r6(k,i)) ! quadratic terms: ! b2(k,i) = 2.*rkm24*r6(k,i)*r7(k,i)**2 + ! | 2.*rk20 ! ! enddo ! k=lev0,lev1-1 ! enddo ! i=lon0,lon1 ! call addfld('BND_R6',' ',' ',r6(lev0:lev1-1,:), ! | 'lev',lev0,lev1-1,'lon',lon0,lon1,lat) ! ! rmo3a, rmo3b = RMOX(TRUE), K = 1/2, 3/2 ! ! rmo3a(:) = rmass_o3 ! rmo3b(:) = rmass_o3 ! end subroutine bndry_comp !----------------------------------------------------------------------- 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 :: 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,nlonp4 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 ! !----------------------------------------------------------------------- ! Begin GSWM subroutines (formerly in bndry_gswm.F): ! ! Get gswm tn lbc perturbations for sub dt. !----------------------------------------------------------------------- ! subroutine lbc_gswm_dt(tnlbc,lon0,lon1,lat0,lat1) ! tbound is replaced by NRL-MSISOO temperatures use input_module,only: step use init_module,only: iter,igswm_mi_di,igswm_mi_sdi, | igswm_nm_di,igswm_nm_sdi use cons_module,only: freq_semidi,tbound implicit none ! Args: integer,intent(in) :: lon0,lon1,lat0,lat1 ! dimension for tnlbc real,intent(out) :: tnlbc(lon0:lon1,lat0:lat1) ! lower boundary condition ! Local: integer :: i,lat real :: rstep complex :: expt,expt2,expta ! ! Calculate exponentials ! rstep = float(step) expt = cexp(ci*freq_semidi*rstep*iter) expt2 = cexp(ci*.5*freq_semidi*rstep*iter) expta = 1. ! ! GSWM migrating diurnal and semi-diurnal: if(igswm_mi_di == 1.and.igswm_mi_sdi == 1) then do lat=lat0,lat1 do i=lon0,lon1 ! tnlbc(i,lat) = gswm_mi_sdi_t(i,lat)+tbound ! semidiurnal tide tnlbc(i,lat) = gswm_mi_sdi_t(i,lat) ! semidiurnal tide tnlbc(i,lat) = tnlbc(i,lat)+ gswm_mi_di_t(i,lat) ! diurnal tide ! don't use annual tides (see timegcm1.2) ! tnlbc(i,lat) = tnlbc(i,lat)+ ! | real(tba(lat)*bnda(i)*expta) ! annual tide enddo ! i=lon0,lon1 enddo ! j=lat0,lat1 ! ! GSWM migrating semi-diurnal: elseif(igswm_mi_di == 0.and.igswm_mi_sdi == 1) then do lat=lat0,lat1 do i=lon0,lon1 ! tnlbc(i,lat) = gswm_mi_sdi_t(i,lat)+tbound ! semidiurnal tide tnlbc(i,lat) = gswm_mi_sdi_t(i,lat) ! semidiurnal tide tnlbc(i,lat) = tnlbc(i,lat)+ | real(tb2(lat)*bnd2(i)*expt2) ! diurnal tide ! tnlbc(i,lat) = tnlbc(i,lat)+ ! | real(tba(lat)*bnda(i)*expta) ! annual tide enddo ! i=lon0,lon1 enddo ! j=lat0,lat1 ! ! GSWM migrating diurnal: elseif(igswm_mi_di == 1.and.igswm_mi_sdi == 0) then do lat=lat0,lat1 do i=lon0,lon1 ! tnlbc(i,lat) = real(tb(lat)*bnd(i)*expt)+tbound ! semidiurnal tide tnlbc(i,lat) = real(tb(lat)*bnd(i)*expt) ! semidiurnal tide tnlbc(i,lat) = tnlbc(i,lat)+ gswm_mi_di_t(i,lat) ! diurnal tide ! tnlbc(i,lat) = tnlbc(i,lat)+ ! | real(tba(lat)*bnda(i)*expta) ! annual tide enddo ! i=lon0,lon1 enddo ! j=lat0,lat1 else ! ! No gswm: do lat=lat0,lat1 do i=lon0,lon1 ! tnlbc(i,lat) = real(tb(lat)*bnd(i)*expt)+tbound ! semidiurnal tide tnlbc(i,lat) = real(tb(lat)*bnd(i)*expt) ! semidiurnal tide tnlbc(i,lat) = tnlbc(i,lat)+ | real(tb2(lat)*bnd2(i)*expt2) ! diurnal tide ! tnlbc(i,lat) = tnlbc(i,lat)+ ! | real(tba(lat)*bnda(i)*expta) ! annual tide enddo ! i=lon0,lon1 enddo ! j=lat0,lat1 endif ! ! GSWM non-migrating diurnal: if(igswm_nm_di == 1) then ! nonmigrating diurnal tide do lat=lat0,lat1 do i=lon0,lon1 tnlbc(i,lat) = tnlbc(i,lat)+ gswm_nm_di_t(i,lat) enddo ! i=lon0,lon1 enddo ! j=lat0,lat1 endif ! ! GSWM non-migrating semi-diurnal: if(igswm_nm_sdi == 1) then ! nonmigrating semidiurnal tide do lat=lat0,lat1 do i=lon0,lon1 tnlbc(i,lat) = tnlbc(i,lat)+ gswm_nm_sdi_t(i,lat) enddo ! i=lon0,lon1 enddo ! j=lat0,lat1 endif ! end subroutine lbc_gswm_dt ! !----------------------------------------------------------------------- ! calculate geopotential height ! z(1) = zb ! subroutine lbc_gswm_addiag(z,lon0,lon1,lat0,lat1) ! use input_module,only: step use init_module,only: iter,igswm_mi_di,igswm_mi_sdi, | igswm_nm_di,igswm_nm_sdi use cons_module,only: freq_semidi,dt use msis_module, only: msis_z ! implicit none ! ! Args: integer,intent(in) :: lon0,lon1,lat0,lat1 ! dimension for tnlbc real,intent(out) :: z(lon0:lon1,lat0:lat1) ! lower boundary condition ! Local: integer :: i,j real :: rstep complex :: expt,expt2,expta ! ! Calculate exponentials ! expt = cexp(ci*freq_semidi*dt*iter) expt2 = cexp(ci*.5*freq_semidi*dt*iter) expta = 1. ! ! am 4/06 annual tides from msis ! GSWM migrating diurnal and semi-diurnal: if(igswm_mi_di == 1.and.igswm_mi_sdi == 1) then do j = lat0,lat1 do i=lon0,lon1 z(i,j) = gswm_mi_sdi_z(i,j) ! semidiurnal tide z(i,j) = z(i,j)+ gswm_mi_di_z(i,j) ! diurnal tide ! z(i,j) = z(i,j)+real(zba(j)*bnda(i)*expta) ! annual tide z(i,j) = z(i,j)+msis_z(i,j) enddo enddo ! ! GSWM migrating semi-diurnal: elseif(igswm_mi_di == 0.and.igswm_mi_sdi == 1) then do j = lat0,lat1 do i=lon0,lon1 z(i,j) = gswm_mi_sdi_z(i,j) ! semidiurnal tide z(i,j) = z(i,j)+real(zb2(j)*bnd2(i)*expt2)! diurnal tide ! z(i,j) = z(i,j)+real(zba(j)*bnda(i)*expta) ! annual tide z(i,j) = z(i,j)+msis_z(i,j) enddo enddo ! ! GSWM migrating diurnal: elseif(igswm_mi_di == 1.and.igswm_mi_sdi == 0) then do j = lat0,lat1 do i=lon0,lon1 z(i,j) = real(zb(j)*bnd(i)*expt) ! semidiurnal tide z(i,j) = z(i,j)+ gswm_mi_di_z(i,j) ! diurnal tide ! z(i,j) = z(i,j)+real(zba(j)*bnda(i)*expta) ! annual tide z(i,j) = z(i,j)+msis_z(i,j) enddo enddo ! ! No gswm: else do j = lat0,lat1 do i=lon0,lon1 z(i,j) = real(zb(j)*bnd(i)*expt) ! semidiurnal tide z(i,j) = z(i,j)+real(zb2(j)*bnd2(i)*expt2) ! diurnal tide ! z(i,j) = z(i,j)+real(zba(j)*bnda(i)*expta) ! annual tide z(i,j) = z(i,j)+msis_z(i,j) enddo enddo endif ! ! GSWM non-migrating diurnal: if(igswm_nm_di == 1) then ! nonmigrating diurnal tide do j = lat0,lat1 do i=lon0,lon1 z(i,j) = z(i,j)+ gswm_nm_di_z(i,j) enddo enddo endif ! ! GSWM non-migrating semi-diurnal: if(igswm_nm_sdi == 1) then ! nonmigrating semidiurnal tide do j = lat0,lat1 do i=lon0,lon1 z(i,j) = z(i,j)+ gswm_nm_sdi_z(i,j) enddo enddo endif end subroutine lbc_gswm_addiag !----------------------------------------------------------------------- subroutine lbc_gswm_duv(unlbc,vnlbc,unlbc_diag,vnlbc_diag, | lon0,lon1,lat,lev0,lev1,expt,expt2,expta) use init_module,only: igswm_mi_di,igswm_mi_sdi, | igswm_nm_di,igswm_nm_sdi ! implicit none ! ! Args: integer,intent(in) :: lon0,lon1,lat,lev0,lev1 ! dimension for lbc complex,intent(in) :: expt,expt2,expta real,intent(out) :: unlbc(lon0:lon1), ! lower boundary condition | vnlbc(lon0:lon1), | unlbc_diag(lev0:lev1,lon0:lon1), ! for output | vnlbc_diag(lev0:lev1,lon0:lon1) ! Local: integer :: i ! ! GSWM migrating diurnal and semi-diurnal: if(igswm_mi_di == 1.and.igswm_mi_sdi == 1) then do i=lon0,lon1 unlbc(i) = gswm_mi_sdi_u(i,lat) ! semidiurnal tide vnlbc(i) = gswm_mi_sdi_v(i,lat) unlbc(i) = unlbc(i) + gswm_mi_di_u(i,lat) ! diurnal tide vnlbc(i) = vnlbc(i) + gswm_mi_di_v(i,lat) unlbc(i) = unlbc(i) + real(uba(lat)*bnda(i)*expta)! annual tide vnlbc(i) = vnlbc(i) + real(vba(lat)*bnda(i)*expta) unlbc_diag(:,i) = unlbc(i) ! for output vnlbc_diag(:,i) = vnlbc(i) enddo ! i=lon0,lon1 ! ! GSWM migrating semi-diurnal: elseif(igswm_mi_di == 0.and.igswm_mi_sdi == 1) then do i=lon0,lon1 unlbc(i) = gswm_mi_sdi_u(i,lat) ! semidiurnal tide vnlbc(i) = gswm_mi_sdi_v(i,lat) unlbc(i) = unlbc(i) + real(ub2(lat)*bnd2(i)*expt2) ! diurnal tide vnlbc(i) = vnlbc(i) + real(vb2(lat)*bnd2(i)*expt2) unlbc(i) = unlbc(i) + real(uba(lat)*bnda(i)*expta) ! annual tide vnlbc(i) = vnlbc(i) + real(vba(lat)*bnda(i)*expta) unlbc_diag(:,i) = unlbc(i) ! for output vnlbc_diag(:,i) = vnlbc(i) enddo ! i=lon0,lon1 ! ! GSWM migrating diurnal: elseif(igswm_mi_di == 1.and.igswm_mi_sdi == 0) then do i=lon0,lon1 unlbc(i) = real(ub(lat)*bnd(i)*expt) ! semidiurnal tide vnlbc(i) = real(vb(lat)*bnd(i)*expt) unlbc(i) = unlbc(i) + gswm_mi_di_u(i,lat) ! diurnal tide vnlbc(i) = vnlbc(i) + gswm_mi_di_v(i,lat) unlbc(i) = unlbc(i) + real(uba(lat)*bnda(i)*expta)! annual tide vnlbc(i) = vnlbc(i) + real(vba(lat)*bnda(i)*expta) unlbc_diag(:,i) = unlbc(i) ! for output vnlbc_diag(:,i) = vnlbc(i) enddo ! i=lon0,lon1 ! ! No gswm: else do i=lon0,lon1 unlbc(i) = real(ub(lat)*bnd(i)*expt) ! Semidiurnal tide vnlbc(i) = real(vb(lat)*bnd(i)*expt) unlbc(i) = unlbc(i) + real(ub2(lat)*bnd2(i)*expt2)! diurnal tide vnlbc(i) = vnlbc(i) + real(vb2(lat)*bnd2(i)*expt2) unlbc(i) = unlbc(i) + real(uba(lat)*bnda(i)*expta)! annual tide vnlbc(i) = vnlbc(i) + real(vba(lat)*bnda(i)*expta) unlbc_diag(:,i) = unlbc(i) ! for output vnlbc_diag(:,i) = vnlbc(i) enddo ! i=lon0,lon1 endif ! ! GSWM non-migrating diurnal: if(igswm_nm_di == 1) then ! nonmigrating diurnal tide do i=lon0,lon1 unlbc(i) = unlbc(i) + gswm_nm_di_u(i,lat) vnlbc(i) = vnlbc(i) + gswm_nm_di_v(i,lat) unlbc_diag(:,i) = unlbc(i) ! for output vnlbc_diag(:,i) = vnlbc(i) enddo ! i=lon0,lon1 endif ! ! GSWM non-migrating semi-diurnal: if(igswm_nm_sdi == 1) then ! nonmigrating semidiurnal tide do i=lon0,lon1 unlbc(i) = unlbc(i) + gswm_nm_sdi_u(i,lat) vnlbc(i) = vnlbc(i) + gswm_nm_sdi_v(i,lat) unlbc_diag(:,i) = unlbc(i) ! for output vnlbc_diag(:,i) = vnlbc(i) enddo ! i=lon0,lon1 endif end subroutine lbc_gswm_duv !----------------------------------------------------------------------- end module bndry_module