subroutine calctends_model_tl(u,v,t,q,oz,cw,pri, & phi_x,phi_y,u_x,u_y,v_x,v_y,t_x,t_y,pri_x,pri_y,q_x,q_y,oz_x,oz_y,cw_x,cw_y, & mype,u_t,v_t,t_t,q_t,oz_t,cw_t,ps_t) !$$$ subprogram documentation block ! . . . . ! subprogram: calctends_model_tl tlm of calctends_model ! prgmmr: kleist org: np20 date: 2005-09-29 ! ! abstract: TLM of routine that compute tendencies for pressure, Tv, u, v ! ! program history log: ! 2005-09-29 kleist ! 2005-10-17 kleist - changes to improve computational efficiency ! 2005-11-21 kleist - add tracer tendencies, use new module ! 2006-04-12 treadon - replace sigi with bk5 ! 2006-04-21 kleist - add divergence tendency parts ! 2006-07-31 kleist - changes to use ps instead of ln(ps) ! 2006-09-21 kleist - add rescaling to divergence tendency formulation ! 2006-10-04 rancic - correct bug in tracer advection terms ! 2007-04-16 kleist - move constraint specific items elsewhere ! 2007-05-08 kleist - add bits for fully generalized vertical coordinate ! 2007-06-21 rancic - add pbl ! 2007-07-26 cucurull - add 3d pressure pri in argument list ; ! move getprs_tl outside calctends_tl; ! call getprs_horiz_tl; ! remove ps from argument list ! 2007-08-08 derber - optimize ! 2008-06-05 safford - rm unused uses ! 2010-02-24 rancic - adjust for use in 4dvar perturbation model ! ! usage: ! input argument list: ! u - zonal wind on subdomain ! v - meridional wind on subdomain ! t - virtual temperature on subdomain ! q - q on subdomain ! oz - ozone on subdomain ! cw - cloud water mixing ratio on subdomain ! pri - 3d pressure ! u_x - zonal derivative of u ! u_y - meridional derivative of u ! v_x - zonal derivative of v ! v_y - meridional derivative of v ! t_x - zonal derivative of t ! t_y - meridional derivative of t !!! ps_x - zonal derivative of ps !!! ps_y - meridional derivative of ps ! pri_x - zonal gradient of pri ! pri_y - meridional gradient of por ! q_x - zonal derivative of q ! q_y - meridional derivative of q ! oz_x - zonal derivative of ozone ! oz_y - meridional derivative of ozone ! cw_x - zonal derivative of cloud water ! cw_y - meridional derivative of cloud water ! mype - task id ! tracer - logical flag if true tracer time derivatives calculated ! ! output argument list: ! u_t - time tendency of u ! v_t - time tendency of v ! t_t - time tendency of t !!! p_t - time tendency of 3d pressure ! p_t - time tendency of surface pressure ! q_t - time tendency of q ! oz_t - time tendency of ozone ! cw_t - time tendency of cloud water ! ! notes: ! TLM check performed & verified on 2005-09-29 by d. kleist ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds,only: r_kind,i_kind use gridmod, only: lat2,lon2,nsig,istart,nlon,nlat,rlats,idvc5,bk5,& wrf_nmm_regional,eta2_ll,regional,region_lat,region_lon,jstart,corlats use constants, only: ione,zero,half,one,two,rd,rcp,grav,rearth,pi,omega use tendsmod, only: coriolis,ctph0,stph0,tlm0 use tends4pertmod, only: curvfct,time_step use nonlinmod, only: bck_u,bck_v,bck_tv,bck_q,bck_oz,bck_cw, & bck_u_lon,bck_u_lat,bck_v_lon,bck_v_lat,bck_tvlat,bck_tvlon, & bck_qlon,bck_qlat,bck_ozlon,bck_ozlat,bck_cwlon,bck_cwlat, & what_bck,prsth_bck,prsum_bck,r_prsum_bck,prdif_bck,r_prdif_bck,pr_xsum_bck,& pr_ysum_bck,rdtop_bck implicit none ! Declare passed variables real(r_kind),dimension(lat2,lon2,nsig) ,intent(in ) :: u,v,t,u_x,u_y,v_x,v_y,& t_x,t_y,q,oz,cw,q_x,q_y,oz_x,oz_y,cw_x,cw_y,phi_x,phi_y real(r_kind),dimension(lat2,lon2,nsig+ione),intent(in ) :: pri,pri_x,pri_y integer(i_kind) ,intent(in ) :: mype real(r_kind),dimension(lat2,lon2,nsig) ,intent( out) :: u_t,v_t,t_t,q_t,oz_t,cw_t real(r_kind),dimension(lat2,lon2) ,intent( out) :: ps_t ! Declare local variables real(r_kind),dimension(lat2,lon2,nsig+ione):: prsth,what real(r_kind),dimension(lat2,lon2,nsig):: prsum,prdif,pr_xsum,pr_ysum real(r_kind),dimension(lat2,lon2,nsig):: t_thor9 !m-----------------------------------------------------------------------------B real(r_kind),dimension(lat2,lon2,nsig):: prdifu,prdifv real(r_kind),dimension(lat2,lon2,nsig):: pgf_x,pgf_y,pgf_xx,pgf_yy real(r_kind),dimension(lat2,lon2,nsig):: pgf1_x,pgf1_y,pgf1_xx,pgf1_yy real(r_kind),dimension(lat2,lon2,nsig):: rdtop,div,prdifu_x,prdifv_y2 real(r_kind),dimension(lat2,lon2):: rdrdx2,rdrdy2 real(r_kind) psum,psum9,psump1,psum9p1 real(r_kind) rr,dlam,dphi,wpdar,pdif real(r_kind) relm,crlm,aph,sph,cph,cc,tph !m-----------------------------------------------------------------------------E real(r_kind) tmp,tmp2,tmp9,tmp9u,tmp9v,tmp9t,tmp9q,tmp9oz,tmp9cw integer(i_kind) i,j,k,ix,jx integer(i_kind) :: jjstart,jjstop integer(i_kind) :: nth,tid,omp_get_num_threads,omp_get_thread_num !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ real(r_kind), parameter:: dts=600.,eps_damp=0.2,fac_nk_damp=2.6 integer(i_kind) k_top,nk_damp,k_damp real(r_kind) dampwt,rdampwt,pihalf,arg,rnk_damp !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! linearized about guess solution, so set it flag accordingly jjstart=1 jjstop=lon2 ! preliminaries: do k=1,nsig do j=jjstart,jjstop do i=1,lat2 prsum(i,j,k)=pri(i,j,k)+pri(i,j,k+ione) prdif(i,j,k)=pri(i,j,k)-pri(i,j,k+ione) pr_xsum(i,j,k)=pri_x(i,j,k)+pri_x(i,j,k+ione) pr_ysum(i,j,k)=pri_y(i,j,k)+pri_y(i,j,k+ione) end do end do end do ! constants if(wrf_nmm_regional) then do j=1,lon2 jx=j+jstart(mype+ione)-2_i_kind jx=min(max(1,jx),nlon) do i=1,lat2 ix=istart(mype+ione)+i-2_i_kind ix=min(max(ix,1),nlat) coriolis(i,j)=two*omega*sin(region_lat(ix,jx)) relm=region_lon(ix,jx)-tlm0 crlm=cos(relm) aph=region_lat(ix,jx) sph=sin(aph) cph=cos(aph) cc=cph*crlm tph=asin(ctph0*sph-stph0*cc) curvfct(i,j)=tan(tph)/rearth end do end do else do j=1,lon2 do i=1,lat2 ix=istart(mype+ione)+i-2_i_kind ix=min(max(ix,2),nlat-1) coriolis(i,j)=corlats(ix) curvfct(i,j)=tan(rlats(ix))/rearth end do end do end if ! Preliminaries for 0) rr=rd/rearth**2 dlam=two*pi/nlon dphi=pi/nlat !!MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM !! wpdar= grav*time_step * 0.09 !m>>>> test !! wpdar= grav*time_step * 0.0 !m>>>> test !!MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM do j=1,lon2 do i=1,lat2 ix=istart(mype+ione)+i-2_i_kind ix=min(max(ix,2),nlat-ione) rdrdx2(i,j)=rr/(cos(rlats(ix))*dlam)**2 rdrdy2(i,j)=rr/dphi**2 end do end do ! 0) TLM of A-grid modification for divergence do k=1,nsig do j=1,lon2 do i=1,lat2 rdtop(i,j,k)=-rd*t(i,j,k)*r_prsum_bck(i,j,k) & -rdtop_bck(i,j,k)*r_prsum_bck(i,j,k)*prsum(i,j,k) pgf_x(i,j,k)=-rdtop_bck(i,j,k)*pr_xsum(i,j,k) & -rdtop(i,j,k)*pr_xsum_bck(i,j,k)+phi_x(i,j,k) pgf_y(i,j,k)=-rdtop_bck(i,j,k)*pr_ysum(i,j,k) & -rdtop(i,j,k)*pr_ysum_bck(i,j,k)+phi_y(i,j,k) end do end do end do !!MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM !! pgf1_x(:,:,:)=zero !! pgf1_y(:,:,:)=zero !! call mp_compact_dlon2(pgf_x,pgf_xx,.false.,nsig,mype) !! call mp_compact_dlat2(pgf_y,pgf_yy,.false.,nsig,mype) !! do k=1,nsig !! do j=1,lon2-ione !! do i=1,lat2 !! psum9=prsum_bck(i,j,k) !! psum9p1=prsum_bck(i,j+ione,k) !! psum=prsum(i,j,k) !! psump1=prsum(i,j+ione,k) !! tmp9=one/(psum9+psum9p1) !! pgf1_x(i,j,k)=rdrdx2(i,j)*tmp9*( & !! (psum9p1-psum9)*(t(i,j+ione,k)+t(i,j,k))+& !! two*tmp9*(bck_tv(i,j+ione,k)+bck_tv(i,j,k))* & !! (psum9*psump1-psum9p1*psum) ) !! end do !! end do !! end do !! do k=1,nsig !! do j=1,lon2 !! do i=1,lat2-ione !! psum9=prsum_bck(i,j,k) !! psum9p1=prsum_bck(i+ione,j,k) !! psum=prsum(i,j,k) !! psump1=prsum(i+ione,j,k) !! tmp9=one/(psum9+psum9p1) !! pgf1_y(i,j,k)=rdrdy2(i,j)*tmp9*( & !! (psum9p1-psum9)*(t(i+ione,j,k)+t(i,j,k))+& !! two*tmp9*(bck_tv(i+ione,j,k)+bck_tv(i,j,k))* & !! (psum9*psump1-psum9p1*psum) ) !! end do !! end do !! end do !! do k=1,nsig !! do j=2,lon2-ione !! do i=1,lat2 !! pgf1_xx(i,j,k)=pgf1_x(i,j,k)-pgf1_x(i,j-ione,k) !! end do !! end do !! end do !! do k=1,nsig !! do j=1,lon2 !! do i=2,lat2-ione !! pgf1_yy(i,j,k)=pgf1_y(i,j,k)-pgf1_y(i-ione,j,k) !! end do !! end do !! end do !!MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM do k=1,nsig do j=1,lon2 do i=1,lat2 prdifu(i,j,k)=prdif_bck(i,j,k)*u(i,j,k)+prdif(i,j,k)*bck_u(i,j,k) prdifv(i,j,k)=prdif_bck(i,j,k)*v(i,j,k)+prdif(i,j,k)*bck_v(i,j,k) end do end do end do call mp_compact_dlon2(prdifu,prdifu_x,.false.,nsig,mype) call mp_compact_dlat2(prdifv,prdifv_y2,.true.,nsig,mype) div(:,:,:)=zero do k=1,nsig do j=1,lon2 do i=1,lat2 !! do j=2,lon2-ione !! do i=2,lat2-ione !!MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM !! div(i,j,k)=wpdar*( pgf1_xx(i,j,k) + pgf1_yy(i,j,k)- & !! pgf_xx(i,j,k) - pgf_yy(i,j,k) ) !!MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM div(i,j,k)=div(i,j,k)+prdifu_x(i,j,k)+prdifv_y2(i,j,k) end do end do end do ! 1) Compute horizontal part of tendency for 3d pressure do j=jjstart,jjstop do i=1,lat2 prsth(i,j,nsig+ione)=zero end do end do do k=nsig,1,-ione do j=jjstart,jjstop do i=1,lat2 prsth(i,j,k)=prsth(i,j,k+ione) - div(i,j,k) end do end do end do ! 1.1) Get horizontal part of temperature tendency for vertical velocity term do k=1,nsig do j=jjstart,jjstop do i=1,lat2 tmp9=-rdtop_bck(i,j,k) tmp =-rdtop(i,j,k) t_t(i,j,k)=-u(i,j,k)*bck_tvlon(i,j,k)-bck_u(i,j,k)*t_x(i,j,k) - & v(i,j,k)*bck_tvlat(i,j,k)-bck_v(i,j,k)*t_y(i,j,k) + & tmp9*rcp*(bck_u(i,j,k)*pr_xsum(i,j,k)+u(i,j,k)*pr_xsum_bck(i,j,k) + & bck_v(i,j,k)*pr_ysum(i,j,k)+v(i,j,k)*pr_ysum_bck(i,j,k) + & prsth(i,j,k) + prsth(i,j,k+ione) ) + & tmp *rcp*(bck_u(i,j,k)*pr_xsum_bck(i,j,k)+bck_v(i,j,k)*pr_ysum_bck(i,j,k)+& prsth_bck(i,j,k)+prsth_bck(i,j,k+ione) ) end do end do end do ! 2) calculate vertical velocity term: z(dp/dz) (zero at top/bottom interfaces) ! if running global, and there is a c(k) coefficient, we call the 'getvvel' ! subroutine if ( (.not.regional) .AND. (idvc5.eq.3)) then ! 1.1) Get horizontal part of temperature tendency for vertical velocity term do k=1,nsig do j=jjstart,jjstop do i=1,lat2 tmp=-rdtop_bck(i,j,k) t_thor9(i,j,k)=-bck_u(i,j,k)*bck_tvlon(i,j,k) - & bck_v(i,j,k)*bck_tvlat(i,j,k) t_thor9(i,j,k)=t_thor9(i,j,k) - & tmp*rcp * ( bck_u(i,j,k)*pr_xsum_bck(i,j,k) + & bck_v(i,j,k)*pr_ysum_bck(i,j,k) + & prsth_bck(i,j,k) + prsth_bck(i,j,k+ione) ) end do end do end do call getvvel_tl(t,t_t,t_thor9,prsth,prdif,what,jjstart,jjstop) else do k=2,nsig do j=jjstart,jjstop do i=1,lat2 if(wrf_nmm_regional) then what(i,j,k)=prsth(i,j,k)-eta2_ll(k)*prsth(i,j,1) else what(i,j,k)=prsth(i,j,k)-bk5(k)*prsth(i,j,1) end if end do end do end do end if ! top/bottom boundary condition: do j=jjstart,jjstop do i=1,lat2 what(i,j,1)=zero what(i,j,nsig+ione)=zero enddo enddo !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! Shiqiu Peng added implicit Rayleigh damping of veritical velocity ! pihalf=pi*half ! k_top=nsig ! nk_damp=k_top/fac_nk_damp+ione ! k_damp=k_top-nk_damp ! rnk_damp=one/nk_damp ! do k=2,nsig ! do j=1,lon2 ! do i=1,lat2 ! if(k.ge.k_damp) then ! arg=pihalf*(k-k_damp)*rnk_damp ! dampwt=eps_damp*dts**sin(arg)**2 ! rdampwt=one/(one+dampwt) ! else ! rdampwt=one ! end if ! what(i,j,k)=what(i,j,k)*rdampwt ! end do ! end do ! end do !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! 3) load actual dp/dt !! do k=1,nsig+ione do j=jjstart,jjstop do i=1,lat2 !! p_t(i,j,k)=prsth(i,j,k)-what(i,j,k) ps_t(i,j)=prsth(i,j,1)-what(i,j,1) end do end do !! end do ! 4) Compute terms for tendencies of wind components & Temperature do k=1,nsig do j=jjstart,jjstop do i=1,lat2 ! horizontal part of momentum equations u_t(i,j,k)= & +coriolis(i,j)*v(i,j,k) & -pgf_x(i,j,k) & +curvfct(i,j)*(bck_u(i,j,k)*v(i,j,k) + u(i,j,k)*bck_v(i,j,k)) & -u(i,j,k)*bck_u_lon(i,j,k)-bck_u(i,j,k)*u_x(i,j,k) & -v(i,j,k)*bck_u_lat(i,j,k)-bck_v(i,j,k)*u_y(i,j,k) v_t(i,j,k)= & -coriolis(i,j)*u(i,j,k) & -pgf_y(i,j,k) & -two*curvfct(i,j)*bck_u(i,j,k)*u(i,j,k) & -u(i,j,k)*bck_v_lon(i,j,k)-bck_u(i,j,k)*v_x(i,j,k) & -v(i,j,k)*bck_v_lat(i,j,k)-bck_v(i,j,k)*v_y(i,j,k) end do !end do i end do !end do j end do !end do k do k=1,nsig do j=jjstart,jjstop do i=1,lat2 ! horizontal advection of "tracer" quantities q_t(i,j,k)=-u(i,j,k)*bck_qlon(i,j,k)-bck_u(i,j,k)*q_x(i,j,k) & -v(i,j,k)*bck_qlat(i,j,k)-bck_v(i,j,k)*q_y(i,j,k) oz_t(i,j,k)=-u(i,j,k)*bck_ozlon(i,j,k)-bck_u(i,j,k)*oz_x(i,j,k) & -v(i,j,k)*bck_ozlat(i,j,k)-bck_v(i,j,k)*oz_y(i,j,k) cw_t(i,j,k)=-u(i,j,k)*bck_cwlon(i,j,k)-bck_u(i,j,k)*cw_x(i,j,k) & -v(i,j,k)*bck_cwlat(i,j,k)-bck_v(i,j,k)*cw_y(i,j,k) end do !end do i end do !end do j end do !end do k ! vertical flux terms do k=1,nsig do j=jjstart,jjstop do i=1,lat2 pdif=prdif(i,j,k)*r_prdif_bck(i,j,k) if (k > 1) then tmp=half*what(i,j,k)*r_prdif_bck(i,j,k) tmp2=half*what_bck(i,j,k)*r_prdif_bck(i,j,k) tmp9u=bck_u(i,j,k-ione)-bck_u(i,j,k) tmp9v=bck_v(i,j,k-ione)-bck_v(i,j,k) tmp9t=bck_tv(i,j,k-ione)-bck_tv(i,j,k) tmp9q=bck_q(i,j,k-ione)-bck_q(i,j,k) tmp9oz=bck_oz(i,j,k-ione)-bck_oz(i,j,k) tmp9cw=bck_cw(i,j,k-ione)-bck_cw(i,j,k) u_t(i,j,k) = u_t(i,j,k) - tmp*tmp9u - & tmp2*( (u(i,j,k-ione)-u(i,j,k))-tmp9u*pdif ) v_t(i,j,k) = v_t(i,j,k) - tmp*tmp9v - & tmp2*( (v(i,j,k-ione)-v(i,j,k))-tmp9v*pdif ) t_t(i,j,k) = t_t(i,j,k)-tmp*tmp9t - & tmp2*( (t(i,j,k-ione)-t(i,j,k))-tmp9t*pdif ) q_t(i,j,k)= q_t(i,j,k) - tmp*tmp9q - & tmp2*( (q(i,j,k-ione)-q(i,j,k))-tmp9q*pdif ) oz_t(i,j,k)= oz_t(i,j,k)-tmp*tmp9oz - & tmp2*( (oz(i,j,k-ione)-oz(i,j,k))-tmp9oz*pdif ) cw_t(i,j,k)= cw_t(i,j,k)-tmp*tmp9cw - & tmp2*( (cw(i,j,k-ione)-cw(i,j,k))-tmp9cw*pdif ) end if if (k < nsig) then tmp=half*what(i,j,k+ione)*r_prdif_bck(i,j,k) tmp2=half*what_bck(i,j,k+ione)*r_prdif_bck(i,j,k) tmp9u=bck_u(i,j,k)-bck_u(i,j,k+ione) tmp9v=bck_v(i,j,k)-bck_v(i,j,k+ione) tmp9t=bck_tv(i,j,k)-bck_tv(i,j,k+ione) tmp9q=bck_q(i,j,k)-bck_q(i,j,k+ione) tmp9oz=bck_oz(i,j,k)-bck_oz(i,j,k+ione) tmp9cw=bck_cw(i,j,k)-bck_cw(i,j,k+ione) u_t(i,j,k)= u_t(i,j,k) - tmp*tmp9u - & tmp2*( (u(i,j,k)-u(i,j,k+ione))-tmp9u*pdif ) v_t(i,j,k)= v_t(i,j,k) - tmp*tmp9v - & tmp2*( (v(i,j,k)-v(i,j,k+ione))-tmp9v*pdif ) t_t(i,j,k)= t_t(i,j,k)- tmp*tmp9t - & tmp2*( (t(i,j,k)-t(i,j,k+ione))-tmp9t*pdif ) q_t(i,j,k)= q_t(i,j,k) - tmp*tmp9q - & tmp2*( (q(i,j,k)-q(i,j,k+ione)) - tmp9q*pdif ) oz_t(i,j,k)= oz_t(i,j,k) - tmp*tmp9oz - & tmp2*( (oz(i,j,k)-oz(i,j,k+ione))-tmp9oz*pdif ) cw_t(i,j,k)= cw_t(i,j,k)-tmp*tmp9cw - & tmp2*( (cw(i,j,k)-cw(i,j,k+ione))-tmp9cw*pdif ) end if end do !end do i end do !end do j end do !end do k call hdiff(u_x,u_y,v_x,v_y,t_x,t_y,u_t,v_t,t_t,mype) call sfcdrag_tl(u,v,t,pri,u_t,v_t) if(.not.wrf_nmm_regional)then do k=1,nsig do j=1,lon2 do i=1,lat2 ix=istart(mype+ione)+i-2_i_kind if (ix == 1 .OR. ix == nlat) then u_t(i,j,k)=zero v_t(i,j,k)=zero end if end do end do end do !end do k end if return end subroutine calctends_model_tl