module bndry_module ! ! Lower boundary conditions diurnal, semi-diurnal, annual, composition, ! 2-day wave, and planetary. ! Lower boundary conditions: diurnal, semi-diurnal are 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 params_module,only: nlat,nlonp4,nlevp1,dlat ! ! 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), | zb3m3(nlat),tb3m3(nlat),ub3m3(nlat),vb3m3(nlat),bnd3m3(nlonp4) real :: b(nlonp4,2,2),fb(nlonp4,2,nlat) real :: zbplanet(nlonp4,nlat) ! planetary wave (sub bndry_planetary) real :: zbkelvin(nlonp4,nlat) ! kelvin wave (sub bndry_kelvin) ! contains !----------------------------------------------------------------------- subroutine lowbound use init_module,only: iday,iter,igswm_mi_di,igswm_mi_sdi use input_module,only: step,ncep_ncfile use cons_module,only: tbound,dz,t0,tgrad,dt use chemrates_module,only: co2mix use zatmos_module,only: zatmos_bndry use cco2gr_module,only: cco2gr,set_cco2_data real :: dday,secs ! t0(1) = tbound t0(2) = tbound+dz*tgrad if (len_trim(ncep_ncfile) <= 0) then secs = amod(float(iter)*dt,86400.) call zatmos_bndry(iday,int(secs)) endif ! ! Call the lower boundary routines. This is called once per model run ! from main tgcm.F. ! 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. call bndry_annual(dday) ! tba,uba,vba,zba call set_cco2_data call cco2gr(co2mix) call bndry_2day ! tb3m3,ub3m3,vb3m3,bnd3m3 call bndry_planetary ! zbplanet call bndry_kelvin ! zbkelvin ! t0(:) = 0. ! end subroutine lowbound !----------------------------------------------------------------------- subroutine bndry_semidiurnal use input_module,only: tide use cons_module,only: pi,cs,cor,tn,t0,atm_amu,gask,re, | grav,freq_semidi,dz,dlamda ! ! Local: integer,parameter :: nalf=19, malf=2 real :: p(nlat,nalf,malf),hough(nlat,5,malf), | cp(nalf/2+1),w(60),xdot(nalf),ydot(nalf),ptscal,theta, | ptjm(2*nlat+1) complex :: dzb(nlat) 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,j,n,jm,l,m,ld,lm1,mm1,nm1 real,external :: sddot ! bhour = tide(6:10) ! ! Zero boundary condition if user input TIDE==0. ! if (all(tide==0.)) then do j=1,nlat zb(j)=0. tb(j)=0. ub(j)=0. vb(j)=0. enddo rlamda = -2.*dlamda bnd(1)=cexp(ci*2.*rlamda) expdlm=cexp(ci*2.*dlamda) do i=2,nlonp4 bnd(i)=bnd(i-1)*expdlm enddo return endif ! tide==0 ! ! Set up hough functions ! 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 ! ! 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 ! 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) ! hough(j,l,ld)=sddot(19,p(j,1,ld),nlat,b(l,1),5) enddo enddo enddo ! ! Generate zb, ub, vb and tb ! do j=1,nlat tb(j)=0. zb(j)=0. dzb(j)=0. enddo 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 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) enddo ! ! Calculate longitudinal structure ! rlamda = -2.*dlamda bnd(1)=cexp(ci*2.*rlamda) expdlm=cexp(ci*2.*dlamda) do i=2,nlonp4 bnd(i)=bnd(i-1)*expdlm enddo end subroutine bndry_semidiurnal !----------------------------------------------------------------------- subroutine bndry_diurnal use input_module,only: tide2 use cons_module,only: pi,cs,cor,tn,t0,dlamda,dz,grav, | gask,freq_semidi,re,atm_amu ! ! Tidal boundary condition for diurnal gravitational mode (1,1) ! integer,parameter :: nalf=19, malf=2 real :: p(nlat,nalf,malf),hough(nlat,5,malf), | cp(nalf/2+1),w(60),xdot(nalf),ydot(nalf),ptscal,theta, | ptjm(2*nlat+1) complex :: dzb(nlat) 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) ! ! Zero boundary condition ! if (tide2(1)==0.) then do j=1,nlat zb2(j)=0. tb2(j)=0. ub2(j)=0. vb2(j)=0. enddo rlamda = -2.*dlamda bnd2(1)=cexp(ci*rlamda) expdlm=cexp(ci*dlamda) do i=2,nlonp4 bnd2(i)=bnd2(i-1)*expdlm enddo return endif ! ! Set up hough functions ! do 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 enddo jm=2*nlat+1 ! ! 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 ! 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) ! hough(j,l,ld)=sdot(19,p(j,1,ld),nlat,b(l,1),1) enddo enddo enddo ! ! Generate zb2, ub2, vb2 and tb2 ! do j=1,nlat tb2(j)=0. zb2(j)=0. dzb(j)=0. enddo do l=1,1 do 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) enddo enddo do 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) enddo ! ! Calculate longitudinal structure ! rlamda = -2.*dlamda bnd2(1)=cexp(ci*rlamda) expdlm=cexp(ci*dlamda) do i=2,nlonp4 bnd2(i)=bnd2(i-1)*expdlm enddo end subroutine bndry_diurnal !----------------------------------------------------------------------- subroutine bndry_annual(dday) use input_module,only: tideann use cons_module,only: pi,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(nlat,nalf,malf),hough(nlat,0:6,2),cp(nalf/2+1) complex dzb(nlat),zzb(nlat) real :: b(6,24),rl(0:6),xdot(nalf),ydot(nalf),ptscal, | ptjm(2*nlat+1),theta,w(nlat) complex cc(0:6,0:6),cl(0:6),expt real,external :: sddot ! data ((b(i,j),i = 1,6),j = 1,12)/ | -0.882922, 0.000000,-0.345087, 0.000000,-0.202228, 0.000000, | 0.000000,-0.930826, 0.000000,-0.301357, 0.000000,-0.152720, | -0.466226, 0.000000, 0.567457, 0.000000, 0.407114, 0.000000, | 0.000000,-0.362673, 0.000000, 0.694431, 0.000000, 0.438014, | -0.055436, 0.000000, 0.711847, 0.000000,-0.163050, 0.000000, | 0.000000,-0.044983, 0.000000, 0.625545, 0.000000,-0.325772, | -0.002909, 0.000000, 0.225723, 0.000000,-0.749160, 0.000000, | 0.000000,-0.002773, 0.000000, 0.186467, 0.000000,-0.723674, | -0.000086, 0.000000, 0.034940, 0.000000,-0.435919, 0.000000, | 0.000000,-0.000103, 0.000000, 0.029425, 0.000000,-0.379254, | -0.000002, 0.000000, 0.003267, 0.000000,-0.122687, 0.000000, | 0.000000,-0.000003, 0.000000, 0.002928, 0.000000,-0.104008/ data ((b(i,J),i = 1,6),j = 13,24)/ | 0.0 , 0.000000, 0.000206, 0.000000,-0.021267, 0.000000, | 0.0 , 0.0 , 0.000000, 0.000202, 0.000000,-0.018228, | 0.0 , 0.0 , 0.000009, 0.000000,-0.002540, 0.000000, | 0.0 , 0.0 , 0.000000, 0.000010, 0.000000,-0.002252, | 0.0 , 0.0 , 0.0 , 0.000000,-0.000223, 0.000000, | 0.0 , 0.0 , 0.0 , 0.0 , 0.000000,-0.000208, | 0.0 , 0.0 , 0.0 , 0.0 ,-0.000015, 0.000000, | 0.0 , 0.0 , 0.0 , 0.0 , 0.000000,-0.000015, | 0.0 , 0.0 , 0.0 , 0.0 ,-0.000001, 0.000000, | 0.0 , 0.0 , 0.0 , 0.0 , 0.000000,-0.000001, | 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.000000, | 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 / data (rl(n),n=1,6)/ | -10.8409E5,-7.0243E5,-2.4874E5,-1.9696E5,-1.0694E5, | -0.9119E5/ ! ! Annual boundary coefficients for lower boundary at z = -17. ! data ((cc(k,n),n=0,6),k=0,3)/ | ( 0.421595E+02, 0.000000E+00),(-0.226810E+00, 0.000000E+00), | ( 0.415692E+00, 0.000000E+00),( 0.547413E-01, 0.000000E+00), | (-0.312575E-01, 0.000000E+00),( 0.280035E-01, 0.000000E+00), | (-0.185304E-01, 0.000000E+00),( 0.300941E-01,-0.734916E-03), | ( 0.100061E+01,-0.168070E+00),(-0.212277E+00, 0.133133E+00), | ( 0.226725E-01,-0.553615E-01),(-0.426632E-01, 0.941169E-03), | (-0.109421E-02,-0.877230E-02),(-0.576746E-02,-0.113662E-01), | ( 0.215525E-01, 0.269012E-01),( 0.212126E-01, 0.377161E-01), | (-0.329359E-01, 0.112051E-01),(-0.117842E-01, 0.145664E-01), | (-0.813882E-02, 0.424029E-01),(-0.705682E-02, 0.100935E-01), | (-0.589241E-02, 0.224780E-01),(-0.956909E-03, 0.169400E-02), | ( 0.283082E-02,-0.579849E-02),( 0.141201E-02,-0.351546E-02), | ( 0.288114E-02, 0.600793E-04),(-0.896708E-04, 0.917448E-04), | ( 0.146016E-02, 0.152191E-02),( 0.911123E-04,-0.141497E-03)/ data ((cc(k,n),n=0,6),k=4,6)/ | (-2.200025E-02,-0.159748E-02),( 0.122976E-03,-0.180022E-03), | ( 0.218866E-02, 0.217849E-02),( 0.361138E-03, 0.264706E-03), | ( 0.803136E-04, 0.595681E-04),(-0.711489E-04, 0.138662E-04), | (-0.205981E-04, 0.237149E-05),( 0.211351E-04, 0.146409E-04), | (-0.195730E-02,-0.810895E-03),(-0.442840E-04, 0.292211E-04), | (-0.123450E-04,-0.538276E-04),(-0.253900E-04,-0.621475E-04), | ( 0.767875E-04,-0.422957E-05),( 0.206095E-04, 0.178031E-05), | (-0.191415E-03, 0.722780E-04),(-0.275580E-04,-0.427594E-04), | ( 0.132229E-03,-0.182378E-03),( 0.341214E-04, 0.355103E-05), | ( 0.326348E-04, 0.543373E-04),(-0.109405E-04, 0.189210E-04), | ( 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) then ! ! Zero boundary condition, except for geopotential: ! do j = 1,nlat zba(j) = cc(0,0)*1./rt2*1.e5 tba(j) = 0. uba(j) = 0. vba(j) = 0. enddo bnda = 1. ! whole array op return endif jm = 2*nlat+1 ! ! Height variation ! do n = 1,6 cl(n) = -csqrt(cmplx(gask/(atm_amu*grav*rl(n))* | (t0(1)*2./7.+(t0(2)-t0(1))/dz)-.25))-.5*ci enddo cl(0) = 0. ! ! Set up legendre polynomials ! ! 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,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 ! ! Now evaluate hough functions ! do l=1,6 do ld=1,2 do j=1,nlat xdot(:) = p(j,:,ld) ydot(:) = b(l,:) hough(j,l,ld)=sddot(24,xdot,ydot) enddo enddo enddo ! ! Hough function of order zero ! do j=1,nlat hough(j,0,1) = 1./rt2 hough(j,0,2) = 0. enddo ! ! Generate zba, tba, uba, vba ! do j=1,nlat zba(j) = 0. tba(j) = 0. uba(j) = 0. vba(j) = 0. enddo ! ! Summation over frequency, k, calculation of phase factor ! do k = 0,6 expt = cexp(ci*float(k)*(dday-1.)*86400.*freq_ann) ! ! Summation over order, n ! do n = 0,6 factor = scale if (k.eq.0.and.n.eq.0) factor = 1. do 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/ | cs(j)*zzb(j)) vba(j) = vba(j)+ci*w(j)*(float(k)*freq_ann*dzb(j)-rm*cor(j)/ | cs(j)*zzb(j)) enddo enddo enddo ! ! Longitudinal structure ! bnda = 1. ! whole array op end subroutine bndry_annual !----------------------------------------------------------------------- subroutine bndry_comp(barm,lev0,lev1,lon0,lon1,lat) use cons_module,only: rmass_o3 use solgar_module,only: xoxlb ! ! This is called from subdomain latitude loop in dynamics. ! Note fb is calculated at longitude subdomains only. ! ! 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 ! ! Args: integer,intent(in) :: | lev0,lev1, ! first and last level indices, this task | lon0,lon1, ! first and last longitude indices, this task | lat ! latitude index real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) :: | barm ! mbar (from addiag) ! ! Local: real,parameter :: pso2lb=0.24 integer :: i,k,kk,nmsk real :: rmo3a(nlonp4),rmo3b(nlonp4) ! ! 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. enddo ! ! rmo3a, rmo3b = RMOX(TRUE), K = 1/2, 3/2 ! rmo3a(:) = rmass_o3 rmo3b(:) = rmass_o3 ! ! rmo3a = RMOX(TRUE), K=0 (EXTRAPOLATION) ! Set fb(1) AND fb(2). Note fb is defined at subdomains only. ! do i = lon0,lon1 rmo3a(i) = 1.5*rmo3a(i)-0.5*rmo3b(i) fb(i,1,lat) = 2.*pso2lb fb(i,2,lat) = 2.*xoxlb(lat)*rmo3a(i)/barm(1,i) ! write(6,"('bndry_comp: i=',i3,' lat=',i3,' xoxlb(lat)=',e12.4, ! | ' rmo3a(i)=',e12.4,' barm(1,i)=',e12.4,' fb(i,2)=', ! | e12.4)") i,lat,xoxlb(lat),rmo3a(i),barm(1,i), ! | fb(i,2) enddo ! write(6,"('bndry_comp: lat=',i3,' fb(:,1,lat)=',/,(6e12.4))") ! | lat,fb(:,1) ! write(6,"('bndry_comp: lat=',i3,' fb(:,2,lat)=',/,(6e12.4))") ! | lat,fb(:,2) end subroutine bndry_comp !----------------------------------------------------------------------- subroutine bndry_2day use input_module,only: tide3m3 use cons_module,only: pi,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 **** Local variables C **** integer,parameter :: nalf=19, malf=2 real :: p(nlat,nalf,malf),hough(nlat,5,malf), | cp(nalf/2+1),xdot(10),ydot(10),ptjm(2*nlat+1),ptscal, | theta,w(nlat) complex :: dzb(nlat) complex zee,cl,expdlm ! ! Where: ! ! B(10) = Coefficients of normalised associated Legendre ! polynomials needed to build Hough(3,-3). ! i.e.: P(l,3), where l = 4,6,8,10,12 (higher ! terms are too small to merit consideration) ! ! RL = equivalent depth for wave ! = 10.5km = 0.105E7 cm ! ! BHOUR = phase of wave measured in hours with constraint: ! 0 .le. BHOUR .lt. p ! ( p = period of wave in hours = 2.07412*24) ! ! TIDE = amplitude of wave in cm ! ! ZEE = complex amplitude incorporating phase (BHOUR) and ! amplitude (TIDE) ! ! EXPDLM = exp(2*pi*i*dlamda), where dlamda is ! longitudinal grid spacing (radians) ! 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 ! ! Check presence of two-day wave ! if (itide3m3.eq.1) then ! ! Two-day tide is present ! ! zee = complex amplitude of wave at time zero. ! = amp * exp(2*pi*i*BHOUR/P2DAY) ! ! cl = vertical wave length ! = ((K*H + dH/ds)/hn - 1/4)**(1/2) -i/2 ! ! where: ! K = 2/7 ! H = scale height (cm) ! hn = equivalent depth for mode = 1.05E7 cm ! s = vertical pressure coordinate ! zee = tide*cexp(ci*freq_3m3*bhour*60.*60.) cl = csqrt(cmplx(gask/(atm_amu*grav*rl)* | (t0(1)*2./7.+(t0(2)-t0(1))/dz)-.25))-.5*ci jm=2*nlat+1 ! ! Set up Hough functions ! ! 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,nlat p(j,nm2,mm2) = ptjm(2*(nlat+1-j)) enddo enddo m = 3 do j=1,nlat 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 ld = 1,2 do j = 1,nlat 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),nlat,b,1) enddo enddo ! ! Generate zb3m3, ub3m3, vb3m3 AND tb3m3 ! do j = 1,nlat 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,nlat 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)) | *zb3m3(j))/ub3m3(j) ub3m3(j) = grav*(cor(j)/freq_3m3*dzb(j)-float(m)/cs(j)* | zb3m3(j))/ub3m3(j) ENDDO ! ! Calculate longitudinal structure ! rlamda = -2.*dlamda bnd3m3(1)=cexp(ci*float(m)*rlamda) expdlm=cexp(ci*float(m)*dlamda) do i=2,nlonp4 bnd3m3(i)=bnd3m3(i-1)*expdlm enddo else ! ! Zero boundary condition, (3,-3) wave is absent. ! do j=1,nlat zb3m3(j)=0. tb3m3(j)=0. ub3m3(j)=0. vb3m3(j)=0. enddo endif end subroutine bndry_2day !----------------------------------------------------------------------- subroutine bndry_planetary use init_module,only: iter,glat use cons_module,only: pi,dlamda,dt ! ! Calculate contribution to ZB from planetary waves ! real,parameter :: amplan=5.0E4, t0plan=2.5E+5 ! real,parameter :: amplan=11.0E4, t0plan=2.5E+5 ! data amplan,t0plan/ 5.0E4,2.5E+5/ ! data amplan,t0plan/ 11.0E4,2.5E+5/ ! ! real,parameter :: frqplan=1.212e-5 ! 6 day real,parameter :: frqplan=3.63e-5 ! 6 day ! ! Local: real :: time,fac1 integer :: j,i,istartlat,iendlat,iterstart real,parameter :: startlat=-82.5 real,parameter :: endlat=-32.5 ! iterstart = (150.*86400.+.1)/dt ! starting iter at day 255 istartlat = 0 iendlat = 0 do j=1,nlat-1 if (startlat >= glat(j) .and. startlat < glat(j+1)) istartlat=j if (endlat > glat(j) .and. endlat <= glat(j+1)) iendlat=j enddo if (istartlat==0 .or. iendlat==0) then write(6,"(/,'>>> bndry_planetary: could not find index to', | ' startlat=',f8.2,': glat=',/,(6f8.2))") startlat,glat call shutdown('startlat') endif zbplanet = 0. ! init whole-array time = (iter-iterstart)*dt fac1 = 1.-exp(-time/t0plan) ! ! 10/8/04 btf: To implement planetary waves, uncomment below code. ! ! do j = istartlat,iendlat ! do i = 1,nlonp4 ! **** PLANETARY WAVE 1 ! zbplanet(i,j) = amplan*(sin((glat(j)-30.)*pi/(60.- ! | (glat(j)-60.)/3.)))**2* ! | (1.-exp(-time/t0plan))*sin(1.*(-pi+(i-3)*dlamda)) ! zbplanet(i,j) = amplan*(sin((glat(j)+30.)*pi/(-60.- ! | (glat(j)+60.)/3.)))**2* ! | fac1*sin(-pi+(i-3)*dlamda+frqplan*time) ! zbplanet(i,j) = amplan*(sin((glat(j)-30.)*pi/(60.- ! | (glat(j)-60.)/3.)))**2* ! | sin(1.*(-pi+(i-3)*dlamda)) ! zbplanet(i,j) = 0. ! ! **** PLANETARY WAVE 2 ! zbplanet(i,j) = amplan*(sin((glat(j)-30.)*pi/(60.- ! | (glat(j)-60.)/3.)))**2* ! | (1.-exp(-time/t0plan))*sin(2.*(-pi+(i-3)*dlamda)) ! zbplanet(i,j) = amplan*(sin((glat(j)-30.)*pi/(60.- ! | (glat(j)-60.)/3.)))**2* ! | sin(2.*(-pi+(i-3)*dlamda)) ! enddo ! enddo ! ! Top north latitude is zero: zbplanet(:,nlat) = 0. zbplanet(:,1) = 0. ! write(6,"('bndry_planet: zbplanet min,max=',2e12.4)") ! | minval(zbplanet),maxval(zbplanet) end subroutine bndry_planetary !----------------------------------------------------------------------- subroutine bndry_kelvin ! ! Calculate contribution to ZB from ultra-fast Kelvin Waves ! use init_module,only: iter,glat use cons_module,only: pi,dt use hist_module,only: nstep,modeltime ! for print only ! ! Local: ! 10/7/04 btf: amplan=1.0e5 resulted in NaN's in most fields after ~1/2 hour. ! 10/7/04 btf: amplan=1.0e4 resulted in NaN's in most fields after ~14 hours. ! 2/1/05 btf: changed t0plan from 2.5e4 to 2.5e5 for timegcm1.2 ! ! real,parameter :: amplan=1.0E5, t0plan=2.5E+4 ! real,parameter :: amplan=1.0E4, t0plan=2.5E+4 ! timegcm1.1 real,parameter :: amplan=1.0E4, t0plan=2.5E+5 ! timegcm1.2 real :: omg3d,wvx1,time integer :: j,i ! omg3d = 2.*pi/86400./3. ! radian frequency of 3 days (rad/sec) wvx1 = 2.*pi/360. ! wavenumber 1 (rad/degree) zbkelvin = 0. ! init whole-array time = iter*dt ! ! 10/8/04 btf: To implement planetary waves, uncomment below code. ! ! do j = 2,nlat-1 ! do i = 1,nlonp4 ! zbkelvin(i,j) = zbkelvin(i,j)+ ! | amplan*exp(-(glat(j)/30.)**2)* ! | sin(omg3d*time-wvx1*(i-3)*5.)* ! | (1.-exp(-time/t0plan)) ! enddo ! ! write(6,"('bndry_kelvin: mtime=',3i4,' nstep=',i4,' dt=',f8.2, ! | ' iter=',i5,' time=',e12.4,' j=',i3,' zbkelvin(:,j)=',/, ! | (6e12.4))") modeltime(1:3),nstep,dt,iter,time,j,zbkelvin(:,j) ! ! enddo zbkelvin(:,1) = 0. ! highest south latitude zbkelvin(:,nlat) = 0. ! highest north latitude ! write(6,"('bndry_kelvin: zbkelvin min,max=',2e12.4)") ! | minval(zbkelvin),maxval(zbkelvin) end subroutine bndry_kelvin ! !----------------------------------------------------------------------- ! Begin GSWM subroutines (formerly in tgcm24/bndry_gswm.F): ! ! Get gswm tn lbc perturbations for sub dt. ! ! Note tbound is not added here, as in tiegcm1. In timegcm, ! zatmos boundary is added in dt intead of tbound. ! !----------------------------------------------------------------------- ! subroutine lbc_gswm_dt(tnlbc,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 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) ! semidiurnal tide tnlbc(i,lat) = tnlbc(i,lat)+ gswm_mi_di_t(i,lat) ! diurnal tide ! Double gswm perturbations as in timegcm1.2: ! tnlbc(i,lat) = gswm_mi_sdi_t(i,lat)*2.0 ! semidiurnal tide ! tnlbc(i,lat) = tnlbc(i,lat)+gswm_mi_di_t(i,lat)*2.0 ! 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) ! 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) ! 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) ! 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 ! 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. ! write(6,"('lbc_gswm_addiag: igswm_mi_di=',i2,' igswm_mi_sdi=', ! | i2,' igswm_nm_di=',i2,' igswm_nm_sdi=',i2)") ! | igswm_mi_sdi,igswm_mi_di,igswm_nm_di,igswm_nm_sdi ! ! 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 ! 5/5/06 btf: double gswm tides, as per Roble (8/04): z(i,j) = gswm_mi_sdi_z(i,j) ! semidiurnal tide ! z(i,j) = gswm_mi_sdi_z(i,j)*2.0 ! semidiurnal tide z(i,j) = z(i,j)+ gswm_mi_di_z(i,j) ! diurnal tide ! z(i,j) = z(i,j)+ gswm_mi_di_z(i,j)*2.0 ! diurnal tide ! 5/5/06 btf: do not add annual tide, as per timegcm1.2 ! z(i,j) = z(i,j)+real(zba(j)*bnda(i)*expta) ! annual tide 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 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 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 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