module lbc ! ! This software is part of the NCAR TIE-GCM. Use is governed by the ! Open Source Academic Research License Agreement contained in the file ! tiegcmlicense.txt. ! ! Calculate lower boundary conditions for T,U,V,Z ! ! 06/13 bae: Added for zm_file use params_module,only: nlon,nlonp2,nlevp1,dlev,zibot use mpi_module,only: lon0,lon1,lat0,lat1 use params_module,only: nlonp4,nlat,nlev,dz use cons_module,only: pi,atm_amu,gask,grav,freq_semidi, | dt,re,dlamda,tgrad,cs,cor,tn use cons_module,only: tbound,zbound use addfld_module,only: addfld implicit none private calc_tzm ! calc_tzm from dt.F in TIMEGCM ! Percentage of tide+zm data to nudge with from lbc to gswm09_nlev or ctmt_nlev real,allocatable,save :: nudge(:) ! dimensioned later in rd_zm integer,save :: nudge_nlev ! Set in input.F to 0, ctmt_nlev, gswm09_nlev ! Use max dimension of 13 for ZM_FILE from 90:2.5:120km (max number) real,dimension(13) :: alt_zm,globtn_zm,globden_zm,globun_zm, | globvn_zm real,dimension(72,13) :: zm_t,zm_d,zm_u,zm_v,zm_dpc,zm_tpert, | zm_zpert ! Arrays on constant p levels for tides plus zonal means (zm) (for lbc,dt,duv.F) real,allocatable,save,dimension(:,:,:) :: ! (lon0:lon1,lat0:lat1,nudge_nlev) | t_zm, u_zm, v_zm, z_zm ! 9/14 add monthly ZM lbc (~95km) from Mack Jones Jr WINDII-HRDI Un,Vn and SABER Tn,gph real,dimension(12,72) :: zm_t12,zm_gph12,zm_u12,zm_v12 ! zonal mean tn,un,vn,z for TIEGCM (was not calc correctly in lbc.F!) ! real,dimension(nlev,nlat) :: tzm_tie,uzm_tie,vzm_tie,zzm_tie ! ! Total lower boundary conditions returned by this module. ! (dimensioned at full global grid, but defined at subdomains only) ! real,dimension(nlonp4,nlat) :: t_lbc, u_lbc, v_lbc, z_lbc, d_lbc ! ! Diurnal and semi-diurnal tidal perturbations using Hough functions: ! complex,dimension(nlat) :: | t_di , u_di , v_di , z_di, t_sdi, u_sdi, v_sdi, z_sdi complex,parameter :: ci=(0.,1.), expta=1. complex :: bnd_sdi(nlonp4), bnd_di(nlonp4) ! ! For bndcmp: real :: b(nlonp4,2,2),fb(nlonp4,2) ! ! This t0 is different than the t0 in cons. real :: t0(nlev+1) ! contains !----------------------------------------------------------------------- subroutine tuvz_lbc(istep,itp) ! HME/CTMT: 2 extra arguements istep,itp in tuvz_lbc called from advance.F ! Model reads itp data and updates its itc data, so get global mean of z(..itp) ! For initial lbc, use zbound=96.3723km (in cm) (and add from ZM and z_pert) ! ! Update lower boundary subdomains of T,U,V,Z ! This is called every timestep from advance, after getgswm and ! before addiag. ! use mpi_module,only: lon0,lon1,lat0,lat1 use init_module,only: iter use input_module,only: step, | gswm_mi_di_ncfile, ! gswm migrating diurnal data file | gswm_mi_sdi_ncfile, ! gswm migrating semi-diurnal data file | gswm_nm_di_ncfile, ! gswm non-migrating diurnal data file | gswm_nm_sdi_ncfile ! gswm non-migrating semi-diurnal data file use gswm_module,only: ! (nlonp4,nlat) | 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 use input_module,only: saber_ncfile, tidi_ncfile use saber_tidi,only: get_saber_tidi,saber_t,saber_z, | tidi_u,tidi_v use hist_module,only: modeltime use input_module,only: gswm09_di_file,gswm09_nlev use gswm09_module,only: getgswm09, nalt9 use input_module,only: ctmt_di_ncfile,ctmt_nlev use ctmt_module,only: getctmt,nalt5 implicit none ! Args: integer,intent(in) :: istep,itp ! Local: integer :: i,j,j2,k,ki,nalt,nalts real :: rstep,fac1 complex :: t_expt_sdi, t_expt_di, uvz_expt_sdi, uvz_expt_di ! ! Calculate exponentials rstep = float(step) t_expt_sdi = cexp(ci*freq_semidi*rstep*iter) t_expt_di = cexp(ci*.5*freq_semidi*rstep*iter) uvz_expt_sdi = cexp(ci*freq_semidi*dt*iter) uvz_expt_di = cexp(ci*.5*freq_semidi*dt*iter) ! ! Set background constants (see cons module): (zbound=96.3723km in cm, T=181.0K) t_lbc(lon0:lon1,lat0:lat1) = tbound u_lbc(lon0:lon1,lat0:lat1) = 0. v_lbc(lon0:lon1,lat0:lat1) = 0. z_lbc(lon0:lon1,lat0:lat1) = zbound d_lbc(lon0:lon1,lat0:lat1) = 0. ! Set nudge_nlev=-1 initially if (istep .eq. 1) nudge_nlev = -1 if (len_trim(saber_ncfile) > 0 .or. | len_trim(tidi_ncfile) > 0) goto 100 if (len_trim(gswm09_di_file) > 0) go to 125 if (len_trim(ctmt_di_ncfile) > 0) go to 150 ! ! Add gswm or Hough mode perturbations: ! GSWM 12076: 141-212K,-96 to +79 m/s Un, -115 to +119 m/s Vn, z 94.10-97.53km ! ! GSWM migrating diurnal: if (len_trim(gswm_mi_di_ncfile) > 0) then do j=lat0,lat1 do i=lon0,lon1 t_lbc(i,j) = t_lbc(i,j)+gswm_mi_di_t(i,j) u_lbc(i,j) = u_lbc(i,j)+gswm_mi_di_u(i,j) v_lbc(i,j) = v_lbc(i,j)+gswm_mi_di_v(i,j) z_lbc(i,j) = z_lbc(i,j)+gswm_mi_di_z(i,j) enddo enddo ! ! Hough mode diurnal: else ! use Hough functions for diurnal tide do j=lat0,lat1 do i=lon0,lon1 t_lbc(i,j) = t_lbc(i,j)+real(t_di(j)*bnd_di(i)*t_expt_di) u_lbc(i,j) = u_lbc(i,j)+real(u_di(j)*bnd_di(i)*uvz_expt_di) v_lbc(i,j) = v_lbc(i,j)+real(v_di(j)*bnd_di(i)*uvz_expt_di) z_lbc(i,j) = z_lbc(i,j)+real(z_di(j)*bnd_di(i)*uvz_expt_di) enddo enddo endif ! ! GSWM migrating semi-diurnal: if (len_trim(gswm_mi_sdi_ncfile) > 0) then do j=lat0,lat1 do i=lon0,lon1 t_lbc(i,j) = t_lbc(i,j)+gswm_mi_sdi_t(i,j) u_lbc(i,j) = u_lbc(i,j)+gswm_mi_sdi_u(i,j) v_lbc(i,j) = v_lbc(i,j)+gswm_mi_sdi_v(i,j) z_lbc(i,j) = z_lbc(i,j)+gswm_mi_sdi_z(i,j) enddo enddo ! ! Hough mode semi-diurnal: else ! use Hough functions for semi-diurnal tide do j=lat0,lat1 do i=lon0,lon1 t_lbc(i,j)=t_lbc(i,j)+real(t_sdi(j)*bnd_sdi(i)*t_expt_sdi) u_lbc(i,j)=u_lbc(i,j)+real(u_sdi(j)*bnd_sdi(i)*uvz_expt_sdi) v_lbc(i,j)=v_lbc(i,j)+real(v_sdi(j)*bnd_sdi(i)*uvz_expt_sdi) z_lbc(i,j)=z_lbc(i,j)+real(z_sdi(j)*bnd_sdi(i)*uvz_expt_sdi) enddo enddo endif ! ! GSWM non-migrating diurnal: if (len_trim(gswm_nm_di_ncfile) > 0) then do j=lat0,lat1 do i=lon0,lon1 t_lbc(i,j) = t_lbc(i,j)+gswm_nm_di_t(i,j) u_lbc(i,j) = u_lbc(i,j)+gswm_nm_di_u(i,j) v_lbc(i,j) = v_lbc(i,j)+gswm_nm_di_v(i,j) z_lbc(i,j) = z_lbc(i,j)+gswm_nm_di_z(i,j) enddo enddo endif ! ! GSWM non-migrating semi-diurnal: if (len_trim(gswm_nm_sdi_ncfile) > 0) then do j=lat0,lat1 do i=lon0,lon1 t_lbc(i,j) = t_lbc(i,j)+gswm_nm_sdi_t(i,j) u_lbc(i,j) = u_lbc(i,j)+gswm_nm_sdi_u(i,j) v_lbc(i,j) = v_lbc(i,j)+gswm_nm_sdi_v(i,j) z_lbc(i,j) = z_lbc(i,j)+gswm_nm_sdi_z(i,j) enddo enddo endif ! ! Add SABER and/or TIDI perturbations: 100 continue if (len_trim(saber_ncfile) > 0 .or. | len_trim(tidi_ncfile) > 0) | call get_saber_tidi(modeltime) if (len_trim(saber_ncfile) > 0) then do j=lat0,lat1 do i=lon0,lon1 t_lbc(i,j) = t_lbc(i,j)+saber_t(i,j) z_lbc(i,j) = z_lbc(i,j)+saber_z(i,j) enddo enddo ! write(6,"('tuvz_lbc: added saber_t,z to t_lbc,z_lbc')") ! write(6,"('tuvz_lbc: saber_z min,max=',2e12.4)") ! | minval(saber_z),maxval(saber_z) ! write(6,"('tuvz_lbc: z_lbc min,max=',2e12.4)") ! | minval(z_lbc),maxval(z_lbc) endif if (len_trim(tidi_ncfile) > 0) then do j=lat0,lat1 do i=lon0,lon1 u_lbc(i,j) = u_lbc(i,j)+tidi_u(i,j) v_lbc(i,j) = v_lbc(i,j)+tidi_v(i,j) enddo enddo ! write(6,"('tuvz_lbc: added tidi_u,v to u_lbc,v_lbc')") endif ! Add GSWM09 (use k=4 for 97.5km for lbc initially) 125 continue if (len_trim(gswm09_di_file) > 0) then ! Must read files first for istep=1 in getgswm09 call getgswm09 (itp) ! 06/13 bae: Read ZM_FILE for istep=1 if gswm09_nlev>0 or set alt_zm if gswm09_nlev=0 if (istep .eq. 1) then nudge_nlev = gswm09_nlev call rd_zm (nalt9) endif call nudge_zm (nalt9,itp,istep) endif ! if (len_trim(gswm09_di_file) > 0) then ! ! Add CTMT: 150 continue if (len_trim(ctmt_di_ncfile) > 0) then ! Must read files first for istep=1 in getctmt call getctmt (itp) ! 06/13 bae: Read ZM_FILE for istep=1 if ctmt_nlev>0 ! alt_zm set in ctmt.F if ctmt_nlev=0 or set alt_zm if ctmt_nlev=0 if (istep .eq. 1) then nudge_nlev = ctmt_nlev call rd_zm (nalt5) ! 9/14 bae: Test reading 95km lbc from Mack Jones Jr WINDII-HRDI Un,Vn and SABER Tn,gph ! nudge_nlev = 1 ! call rd_zm_mack endif call nudge_zm (nalt5,itp,istep) ! do j=lat0,lat1 ! do i=lon0,lon1 ! t_lbc(i,j) = t_lbc(i,j)+gswm_nm_di_t(i,j) ! u_lbc(i,j) = u_lbc(i,j)+gswm_nm_di_u(i,j) ! v_lbc(i,j) = v_lbc(i,j)+gswm_nm_di_v(i,j) ! z_lbc(i,j) = z_lbc(i,j)+gswm_nm_di_z(i,j) ! enddo ! enddo endif ! for CTMT ! Print out min/max values of lbc if (istep<3) write(6,"('tuvz_lbc: t_lbc min,max=',2e12.4, | ' u_lbc=',2e12.4,' v_lbc=',2e12.4,' z_lbc=',2e12.4)") | minval(t_lbc),maxval(t_lbc), | minval(u_lbc),maxval(u_lbc), | minval(v_lbc),maxval(v_lbc), | minval(z_lbc),maxval(z_lbc) ! ! Save 2d boundaries to secondary histories: call addfld('T_LBC','T_LBC',' ',t_lbc(lon0:lon1,lat0:lat1), | 'lon',lon0,lon1,'lat',lat0,lat1,0) call addfld('U_LBC','U_LBC',' ',u_lbc(lon0:lon1,lat0:lat1), | 'lon',lon0,lon1,'lat',lat0,lat1,0) call addfld('V_LBC','V_LBC',' ',v_lbc(lon0:lon1,lat0:lat1), | 'lon',lon0,lon1,'lat',lat0,lat1,0) call addfld('Z_LBC','Z_LBC',' ',z_lbc(lon0:lon1,lat0:lat1), | 'lon',lon0,lon1,'lat',lat0,lat1,0) end subroutine tuvz_lbc !----------------------------------------------------------------------- subroutine bndry_semidiurnal ! ! Lower boundary conditions for semi-diurnal tide, using Hough functions. ! This is called once per run from init, and returns t_sdi, u_sdi, v_sdi, ! z_sdi at nlat latitudes. ! use input_module,only: tide implicit none ! ! 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/ | 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/ data rl/7.8519E5, 3.6665E5, 2.1098E5, 1.3671E5, 0.9565E5/ real,external :: sddot ! util.F ! ! t0 is local to this module (different than the t0 in cons) t0(:) = 0. t0(1) = tbound t0(2) = tbound+dz*tgrad ! ! Longitudinal structure: rlamda = -2.*dlamda bnd_sdi(1)=cexp(ci*2.*rlamda) expdlm=cexp(ci*2.*dlamda) do i=2,nlonp4 bnd_sdi(i)=bnd_sdi(i-1)*expdlm enddo ! ! Zero out if user did not provide amp/phase: if (all(tide==0.)) then t_sdi = 0. u_sdi = 0. v_sdi = 0. z_sdi = 0. return endif ! bhour = tide(6:10) 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 ! ! Set up hough functions (see 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 ! ! Define module data: do j=1,nlat t_sdi(j)=0. z_sdi(j)=0. dzb(j)=0. enddo do l=1,5 do j=1,nlat z_sdi(j)=z_sdi(j)+zee(l)*hough(j,l,1) dzb(j)=dzb(j)+zee(l)*hough(j,l,2) t_sdi(j)=t_sdi(j)+ci*atm_amu*grav/gask*zee(l)*cl(l)* | hough(j,l,1) enddo enddo do j=1,nlat u_sdi(j)=freq_semidi*re*(1.-(cor(j)/freq_semidi)**2) v_sdi(j)=ci*grav*(dzb(j)-2.*cor(j)/(freq_semidi*cs(j))* | z_sdi(j))/u_sdi(j) u_sdi(j)=grav*(cor(j)/freq_semidi*dzb(j)-2./cs(j)* | z_sdi(j))/u_sdi(j) enddo ! write(6,"('bndry_semidiurnal: t_sdi min,max=',2e12.4)") ! | minval(real(t_sdi)),maxval(real(t_sdi)) ! write(6,"('bndry_semidiurnal: u_sdi min,max=',2e12.4)") ! | minval(real(u_sdi)),maxval(real(u_sdi)) ! write(6,"('bndry_semidiurnal: v_sdi min,max=',2e12.4)") ! | minval(real(v_sdi)),maxval(real(v_sdi)) ! write(6,"('bndry_semidiurnal: z_sdi min,max=',2e12.4)") ! | minval(real(z_sdi)),maxval(real(z_sdi)) end subroutine bndry_semidiurnal !----------------------------------------------------------------------- subroutine bndry_diurnal ! ! Lower boundary conditions for diurnal tide, using Hough functions. ! This is called once per run from init, and returns t_di, u_di, v_di, ! z_di at nlat latitudes. ! 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),bhour(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/ | 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/ data rl/0.6909E5/ real,external :: sddot ! in util.F ! ! t0 is local to this module (different than the t0 in cons) t0(:) = 0. t0(1) = tbound t0(2) = tbound+dz*tgrad ! ! Calculate longitudinal structure rlamda = -2.*dlamda bnd_di(1)=cexp(ci*rlamda) expdlm=cexp(ci*dlamda) do i=2,nlonp4 bnd_di(i)=bnd_di(i-1)*expdlm enddo ! ! Zero out if user did not provide amp/phase: if (all(tide2==0.)) then t_di = 0. u_di = 0. v_di = 0. z_di = 0. return endif bhour(1) = tide2(2) pik = 3.14159265358979312 do n=1,1 zee(n)=tide2(n)*cexp(ci*pik*bhour(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 ! ! Set up hough functions: ! 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 ! ! Generate t_di, u_di, v_di, z_di: do j=1,nlat t_di(j)=0. z_di(j)=0. dzb(j)=0. enddo do l=1,1 do j=1,nlat z_di(j)=z_di(j)+zee(l)*hough(j,l,1) dzb(j)=dzb(j)+zee(l)*hough(j,l,2) t_di(j)=t_di(j)+ci*atm_amu*grav/gask*zee(l)*cl(l)*hough(j,l,1) enddo enddo do j=1,nlat u_di(j)=.5*freq_semidi*re*(1.-(cor(j)/(.5*freq_semidi))**2) v_di(j)=ci*grav*(dzb(j)-cor(j)/(.5*freq_semidi*cs(j))*z_di(j))/ | u_di(j) u_di(j)=grav*(cor(j)/(.5*freq_semidi)*dzb(j)-1./cs(j)*z_di(j))/ | u_di(j) enddo ! write(6,"('bndry_diurnal: t_di min,max=',2e12.4)") ! | minval(real(t_di)),maxval(real(t_di)) ! write(6,"('bndry_diurnal: u_di min,max=',2e12.4)") ! | minval(real(u_di)),maxval(real(u_di)) ! write(6,"('bndry_diurnal: v_di min,max=',2e12.4)") ! | minval(real(v_di)),maxval(real(v_di)) ! write(6,"('bndry_diurnal: z_di min,max=',2e12.4)") ! | minval(real(z_di)),maxval(real(z_di)) end subroutine bndry_diurnal !----------------------------------------------------------------------- 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 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 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. enddo end subroutine bndef !----------------------------------------------------------------------- subroutine rd_zm(naltzm) ! 06/13 bae: Read zonal mean file for nudging CTMT or GSWM09 tides C or set alt_zm if nudge_nlev=0 use input_module,only: zm_file integer,external :: nextlu ! Args: integer,intent(in) :: naltzm ! Local: integer :: nalt real,dimension(nlon,nlat) :: utid integer :: il,jl,nn,lu_zm,nd_zm character(len=80) :: zm_char real,dimension(72) :: glat72 integer :: istat,k,i real :: zp(nlevp1),bot_alt ! Set alt_zm if nudge_nlev=0 if (nudge_nlev==0) then bot_alt = 95. if (naltzm==13) bot_alt = 90. do i=1,naltzm alt_zm(i) = bot_alt + (i-1)*2.5 enddo write (6,"(1x,'alt_zm=',13f6.1)") (alt_zm(i),i=1,naltzm) else ! If nudge_nlev>0, read zonal mean (zm) file (=1 means lb only, but still use zm) allocate(nudge(nudge_nlev),stat=istat) lu_zm = nextlu() write (6,"(1x,'lu_zm, zm_file =',i3,1x,a)") lu_zm, | trim(zm_file) open(lu_zm,file=trim(zm_file),status='old') do nalt=1,naltzm write (6,"(1x,'nalt=',i2)") nalt read(lu_zm,"(f5.1,i4)") alt_zm(nalt),nd_zm write (6,"(1x,'alt=',f5.1)") alt_zm(nalt) read(lu_zm,"(f7.2,e14.6,2f8.2)") globtn_zm(nalt), | globden_zm(nalt),globun_zm(nalt),globvn_zm(nalt) write (6,"(1x,'globtn,den =',f7.2,e14.7)") globtn_zm(nalt), | globden_zm(nalt) read(lu_zm,"(a)") zm_char ! write (6,"(1x,'zm_char =',a)") zm_char ! Read double resolution zonal means of glat,Tn(K),Den(g/cm3),Un(m/s),Vn(m/s) do jl=1,72 read(lu_zm,"(f6.2,f7.1,e13.5,2f7.1)") glat72(jl), | zm_t(jl,nalt),zm_d(jl,nalt),zm_u(jl,nalt),zm_v(jl,nalt) ! 08/14: Convert Un and Vn to cm/s from m/s zm_u(jl,nalt) = zm_u(jl,nalt)*100. zm_v(jl,nalt) = zm_v(jl,nalt)*100. ! 10/14: Calculate zm_tpert and zm_dpc (want pert from 1.) zm_tpert(jl,nalt) = zm_t(jl,nalt) - globtn_zm(nalt) zm_dpc(jl,nalt) = 1. - zm_d(jl,nalt)/globden_zm(nalt) if (jl==36) | write (6,"(1x,'jl glat72 zm_t =',i2,f6.2,f7.1)") jl, | glat72(jl),zm_t(jl,nalt) enddo ! Print out min/max values of zm(1:nalt9) write(6,"('zm min,max den',2e12.4, | ' un=',2e12.4,' vn=',2e12.4,' tn=',2e12.4)") | minval(zm_d(:,nalt)),maxval(zm_d(:,nalt)), | minval(zm_u(:,nalt)),maxval(zm_u(:,nalt)), | minval(zm_v(:,nalt)),maxval(zm_v(:,nalt)), | minval(zm_t(:,nalt)),maxval(zm_t(:,nalt)) enddo close(lu_zm) ! Calculate nudging factor once per run (nudge_nlev is from namelist C as gswm09_nlev or ctmt_nlev): do k=1,nlevp1 zp(k) = zibot+(k-1)*dlev enddo do k=1,nudge_nlev if (nudge_nlev==1) then ! 100% lb only, no nudging nudge(k) = 0. else nudge(k) = | cos(pi/2.*((zp(k)-zp(1))/(zp(nudge_nlev)-zp(1))))**2 endif enddo write(6,"('rd_zm: nudge_nlev=',i4,' nudge=', | /,(5f15.5))") nudge_nlev,nudge endif ! if (nudge_nlev==0) then,else end subroutine rd_zm !----------------------------------------------------------------------- subroutine nudge_zm (naltzm,itp,istep) use init_module,only: glon,glat ! added for output to fort.40 for nproc=1 ! interface p levels for Z,DEN and midpoint p levels for TN,UN,VN use fields_module,only: z, ! z(nlev,nlon,nlat,2) usually allocated | tn,un,vn use input_module,only: ictmt,igswm09 use ctmt_module,only: ctmt_u,ctmt_v,ctmt_t,ctmt_d use gswm09_module,only: gswm09_z, gswm09_t, gswm09_u, gswm09_v ! interface p levels for Z,DEN and midpoint p levels for TN,UN,VN use fields_module,only: z, ! z(nlev,nlon,nlat,2) usually allocated | tn,un,vn #ifdef MPI use mpi_module,only: mp_periodic_f2d #endif ! Args: integer,intent(in) :: naltzm,itp,istep ! Local: real :: xkm,fac1,xkmn,gphglob integer :: i,j,j2,k,ki,nalt,nalts,ktop,n97p5 ! NOTE: May not have right dimensions for altitude/p level w calc_tzm subroutine! ! real,dimension(1:nalt5+1,lat0:lat1) :: tzm_tie,uzm_tie,vzm_tie real,dimension(14,lat0:lat1) :: tzm_tie,uzm_tie,vzm_tie real,dimension(lon0:lon1,lat0:lat1) :: t_pert,z_pert,d_pert, | u_pert,v_pert real :: fctmt(lon0:lon1,lat0:lat1,3) ! for mp calls ! TEMP 10/14: real,dimension(72) :: zmzpert,zmtpert,zmdpc ! For nudging, try to put zm in u,v,den as well for lb and calc for alts>lbc ! Only nudge z at lbc - so only need to use den and cal_ctmt_z for ki=1 ! Need to interpolate first to z_pert+zbound to interpolate between 95-97.5km ! Could find z_pert again ! Nudge and revise for actual z of ctmt_nlev p levels, going from top to bottom (ki) ! TEMP: The zonal means are not used ! Calculate TIEGCM zonal means (zm) for TN, UN and VN ! calc_tzm is not working! istep==1, minvals ALL zero, but tn(1,1,17,itp) OK!!? ! call calc_tzm(tn(:,lon0:lon1,lat0:lat1,itp),tzm_tie, ! | 1,nudge_nlev,lon0,lon1,lat0,lat1) ! call calc_tzm(un(:,lon0:lon1,lat0:lat1,itp),uzm_tie, ! | 1,nudge_nlev,lon0,lon1,lat0,lat1) ! call calc_tzm(vn(:,lon0:lon1,lat0:lat1,itp),vzm_tie, ! | 1,nudge_nlev,lon0,lon1,lat0,lat1) ! if (istep<3) then ! do i=1,naltzm+1 ! write (6,"('zm at nalt min t,u,v =',i2,8e12.5)") i, ! | minval(tzm_tie(i,:)),minval(uzm_tie(i,:)), ! | minval(vzm_tie(i,:)) ! enddo ! endif ! Want to add z_pert to zbound=136.291*1.e5/sqrt(2.)~96.37km ! CTMT 12076: 155-205K,-35 to +43 m/s Un, -64 to +47 m/s Vn, z 95.24-97.94km, rho +/-15% ! CTMT 12080 w t_zm: 151-200K or -25to+24K, -34to+43m/s Un, -63to+45m/s Vn, ! z 95.3-97.9km or -1.065to+1.545km, -14.98to+15.68% Nden ! For 13 alts from 90-120km, #4=97.5km at the approx lbc n97p5 = 4 if (istep<2 .and. ictmt==1) | write(6,"('tuvz_lbc 97.5km: ctmt_rho,un,vn,tn ', | 'pert =',4e12.4,4f8.0,2f7.1)") | minval(ctmt_d(:,:,n97p5)),maxval(ctmt_d(:,:,n97p5)), | minval(ctmt_u(:,:,n97p5)),maxval(ctmt_u(:,:,n97p5)), | minval(ctmt_v(:,:,n97p5)),maxval(ctmt_v(:,:,n97p5)), | minval(ctmt_t(:,:,n97p5)),maxval(ctmt_t(:,:,n97p5)) if (istep<2 .and. igswm09==1) | write(6,"('tuvz_lbc 97.5km: gswm09_z,un,vn,tn ', | 'pert =',4e12.4,4f8.0,2f7.1)") | minval(gswm09_z(:,:,n97p5)),maxval(gswm09_z(:,:,n97p5)), | minval(gswm09_u(:,:,n97p5)),maxval(gswm09_u(:,:,n97p5)), | minval(gswm09_v(:,:,n97p5)),maxval(gswm09_v(:,:,n97p5)), | minval(gswm09_t(:,:,n97p5)),maxval(gswm09_t(:,:,n97p5)) ! Reminder that z_lbc (and t_lbc) only defined normally as constants (in tuvz_lbc) ! Set background constants (see cons module): (zbound=96.3723km in cm, T=181.0K) t_lbc(lon0:lon1,lat0:lat1) = tbound u_lbc(lon0:lon1,lat0:lat1) = 0. v_lbc(lon0:lon1,lat0:lat1) = 0. z_lbc(lon0:lon1,lat0:lat1) = zbound d_lbc(lon0:lon1,lat0:lat1) = 0. ! 10/14 bae: Find geopotential ht from ZM and save to fort.41 for nproc=1 write (41,"(1x,10f8.1)") glat do k=1,naltzm write (41,"(1x,i4,2f8.1,e14.5)") k,alt_zm(k), | globtn_zm(k),globden_zm(k) ! Calculate z_pert for ZM do j=lat0,lat1 zmtpert(j) = zm_tpert(j,k) zmdpc(j) = zm_dpc(j,k) enddo call cal_ctmt_z(istep,itp,t_lbc,z_pert,t_pert,d_pert, | 1,globtn_zm(k),zmzpert,zmtpert,zmdpc,k) ! k not correct for zm 90-120km do j=lat0,lat1 zm_zpert(j,k) = zmzpert(j) enddo enddo xkmn = 0. xkm = 0. ! Define ktop to be at least 1 so when nudge_nlev=0, k=1 ktop = max(nudge_nlev,1) do k=1,ktop ki = ktop + 1 - k ! Use .._lbc for higher altitudes and then revise to .._lbc at lbc (bottom) ! 9/14: Find global mean of previous z(ki,i,j,itp) to interpolate ZM or _pert to do j=lat0,lat1 j2 = 2*j do i=lon0,lon1 ! NOTE: alt_zm should be the same altitude range as alt_ctmt,gswm of 13 hts from 90-105km ! Interpolate alt_zm to z potential ht nalts = 0 do nalt=2,naltzm if (z(k,i,j,itp)*1.e-5>=alt_zm(nalt-1) .and. | z(k,i,j,itp)*1.e-5<=alt_zm(nalt)) nalts=nalt enddo ! Extrapolate if zalt_zm if (nalts==0) then if (alt_zm(1)>z(k,i,j,itp)*1.e-5) nalts = 2 if (alt_zm(naltzm)