module scanslt !----------------------------------------------------------------------- ! ! Module to handle Semi-Lagrangian transport in the context of ! Semi-Lagrangian Dynamics. ! !----------------------------------------------------------------------- ! ! $Id$ ! !----------------------------------------------------------------------- use shr_kind_mod, only: r8 => shr_kind_r8 use pmgrid, only: plon, plev, plevp, plat, beglat, endlat use constituents, only: pcnst use abortutils, only: endrun use hycoef, only: hyai, hybi, hypm, hypi, ps0, hybm, nprlev, hybd use perf_mod use cam_logfile, only: iulog !----------------------------------------------------------------------- implicit none !----------------------------------------------------------------------- private ! ! Public interfaces ! public slt_alloc ! Initialization public slt_initial ! Initialization public slt_run_setup ! setup internal data for run public slt_run ! run method public slt_final ! Finalization ! ! Public data ! public plond, platd, i1, j1, nxpt, beglatex, endlatex, nlonex, jintmx integer, parameter :: nxpt = 1 ! no.of pts outside active domain of interpolnt integer, parameter :: jintmx = 2 ! # of extra lats in polar region integer, parameter :: plond = plon + 1 + 2*nxpt ! slt extended domain longitude integer, parameter :: platd = plat + 2*nxpt + 2*jintmx ! slt extended domain lat. integer, parameter :: i1 = 1 + nxpt ! model starting longitude index integer, parameter :: j1 = 1 + nxpt + jintmx ! model starting latitude index integer, parameter :: numbnd = nxpt + jintmx ! no.of lats passed N/S of forecast lat integer, parameter :: plndlv = plond*plev ! Length of multilevel 3-d field slice integer :: numlatsex ! number of ext lats owned by a given proc integer :: beglatex ! extended grid beglat integer :: endlatex ! extended grid endlat #if ( ! defined SPMD ) parameter (numlatsex = platd) parameter (beglatex = 1) parameter (endlatex = platd) #endif integer :: nlonex(platd) ! Number of longitudes for each latitude on ! extended grid integer, public, parameter :: plondfft = plon + 2 public numbnd, plndlv ! ! SPMD decomposition information ! integer, allocatable :: cutex(:,:) ! extended partition integer :: neighs ! number of south neighbors to comm guardcells integer, allocatable :: neighs_proc(:) ! sorted south process neighbors integer :: neighn ! number of north neighbors to comm guardcells integer, allocatable :: neighn_proc(:) ! sorted north process neighbors ! ! Public types ! public advection_state type advection_state real(r8), pointer :: u3(:,:,:) ! u wind real(r8), pointer :: v3(:,:,:) ! v wind real(r8), pointer :: t3(:,:,:) ! Temperature real(r8), pointer :: q3(:,:,:,:) ! Specific humidity real(r8), pointer :: u3sld(:,:,:) ! u-wind (time n) (used for advection) real(r8), pointer :: v3sld(:,:,:) ! v-wind (time n) (used for advection) real(r8), pointer :: lnpssld(:,:,:) ! RHS Ps term for SLD real(r8), pointer :: prhssld(:,:,:) ! RHS Ps term for SLD real(r8), pointer :: etadot(:,:,:) ! vertical motion (on previous time-step) end type advection_state public phigs, levkntl ! Can become private as soon as slttraj is part of scanslt public engy1lat, hw1lat ! Used in scan2, restart, and linemsdyn public epssld public qfcst ! Used in restart real(r8) :: hw1lat (pcnst,plat) ! lat contribution to const. mass integral real(r8) :: engy1lat(plat) ! lat contribution to total energy integral real(r8), parameter :: epssld = 0.2_r8 ! "epsilon" for SLD decentering algorithm real(r8), allocatable :: qfcst(:,:,:,:) ! slt forecast of moisture and constituents ! ! Global integrals ! public hw1, hw2, hw3, alpha ! Needed for restart real(r8) :: hw1 (pcnst) ! Pre-SLT global integral of constitent real(r8) :: hw2 (pcnst) ! Post-SLT global integral of const. real(r8) :: hw3 (pcnst) ! Global integral for denom. of expr. for alpha real(r8) :: alpha (pcnst) ! alpha(m) = ( hw1(m) - hw2(m) )/hw3(m) ! ! Private module data !----------------------------------------------------------------------- ! !------------------------------Parameters------------------------------- ! integer, parameter :: pmap = 20000 ! max dimension of evenly spaced vert. grid used ! ! by SLT code to map the departure pts into true ! ! model levels. ! integer :: kdpmpf (pmap) ! artificial full vert grid indices integer :: kdpmph (pmap) ! artificial half vert grid indices real(r8) :: lam (plond,platd) ! longitude coords of extended grid real(r8) :: phi (platd) ! latitude coords of extended grid real(r8) :: dphi (platd) ! latitude intervals (radians) real(r8) :: sinlam (plond,platd) ! sin(lam) model domain only real(r8) :: coslam (plond,platd) ! cos(lam) model domain only real(r8) :: lbasdy (4,2,platd) ! latitude derivative weights real(r8) :: lbasdz (4,2,plev) ! vert (full levels) deriv wghts real(r8) :: lbassd (4,2,plevp) ! vert (half levels) deriv wghts real(r8) :: lbasiy (4,2,platd) ! Lagrange cubic interp wghts (lat.) real(r8) :: lbasiz (4,2,plev) ! Lagrange cubic interp wghts (vert) real(r8) :: lbassi (4,2,plevp) ! Lagrange cubic interp wghts (vert) real(r8) :: detai (plevp) ! intervals between vert half levs. real(r8) :: dlam (platd) ! longitudinal grid interval (radians) real(r8) :: etaint (plevp) ! vertical coords at interfaces integer :: kstep ! Keep track of number of steps ran real(r8) :: gamma (plev,plev) ! SLD coefficient real(r8) :: levknt (plev,plev) ! Counter for departure point binning statistics real(r8) :: levkntl(plev,plev,plat) ! Counter for departure point binning statistics real(r8) :: phigs ! Latitude cutoff for using local geodesic algorithm contains ! !----------------------------------------------------------------------- ! #if (defined SPMD) subroutine slt_spmd() !----------------------------------------------------------------------- ! ! Purpose: ! Setup SPMD information regarding SLT halos will be used in bndexch. ! !----------------------------------------------------------------------- use spmd_utils, only: iam, masterproc, npes use spmd_dyn, only: cut, dyn_npes, dyn_npes_stride ! ! Local variables ! integer :: lat, procid ! Loop indices integer :: procid_s ! strided process id integer :: neighn_minlat(plat) ! minimum latitude in north neighbor integer :: neighs_maxlat(plat) ! maximum latitude in south neighbor ! ! The extended regions are simply "numbnd" wider at each ! side. The extended region do not go beyond 1 and plat, though ! allocate (cutex(2,0:npes-1)) cutex(1,0:npes-1) = 1 cutex(2,0:npes-1) = 0 ! ! Initialization for inactive processes ! beglatex = 1 endlatex = 0 numlatsex = 0 ! do procid=0,dyn_npes-1 procid_s = dyn_npes_stride*procid cutex(1,procid_s) = cut(1,procid_s) - numbnd cutex(2,procid_s) = cut(2,procid_s) + numbnd if (iam == procid_s) then beglatex = cutex(1,procid_s) + numbnd endlatex = cutex(2,procid_s) + numbnd numlatsex = endlatex - beglatex + 1 end if end do ! ! Determine neighbor processors needed for boundary communication. ! North first. ! neighn = 0 neighn_minlat(:) = -1 do procid=0,dyn_npes-1 procid_s = dyn_npes_stride*procid if (procid_s /= iam) then if ((cut(1,procid_s) > cut(2,iam)) .and. & (cut(1,procid_s) <= cut(2,iam)+numbnd)) then neighn_minlat(cut(1,procid_s)) = procid_s neighn = neighn + 1 endif endif enddo ! ! Sort north processes by increasing latitude ! allocate (neighn_proc (neighn)) neighn = 0 do lat=1,plat if (neighn_minlat(lat) /= -1) then neighn = neighn + 1 neighn_proc(neighn) = neighn_minlat(lat) endif enddo ! ! South next. ! neighs = 0 neighs_maxlat(:) = -1 do procid=0,dyn_npes-1 procid_s = dyn_npes_stride*procid if (procid_s /= iam) then if ((cut(2,procid_s) < cut(1,iam)) .and. & (cut(2,procid_s) >= cut(1,iam)-numbnd)) then neighs_maxlat(cut(2,procid_s)) = procid_s neighs = neighs + 1 endif endif enddo ! ! Sort south processes by decreasing latitude ! allocate (neighs_proc (neighs)) neighs = 0 do lat=plat,1,-1 if (neighs_maxlat(lat) /= -1) then neighs = neighs + 1 neighs_proc(neighs) = neighs_maxlat(lat) endif enddo ! if (masterproc) then write(iulog,*)'-----------------------------------------' write(iulog,*)'Number of lats passed north & south = ',numbnd write(iulog,*)'Node Extended Partition' write(iulog,*)'-----------------------------------------' do procid=0,npes-1 write(iulog,200) procid,cutex(1,procid), cutex(2,procid) 200 format(i3,4x,i3,'-',i3) end do end if ! write(iulog,*)'iam=',iam,'Number of south neighbors needed for bndry exchange = ',neighs ! write(iulog,*)'iam=',iam,'Number of north neighbors needed for bndry exchange = ',neighn end subroutine slt_spmd #endif ! !----------------------------------------------------------------------- ! subroutine slt_alloc() !----------------------------------------------------------------------- ! ! Purpose: ! slt initialization ! !----------------------------------------------------------------------- use infnan, only: bigint, nan nlonex(:) = bigint allocate (qfcst(plon,plev,pcnst,beglat:endlat)) qfcst (:,:,:,:) = nan end subroutine slt_alloc ! !----------------------------------------------------------------------- ! subroutine slt_initial( detam, gw, cwava, etamid, adv_state ) !----------------------------------------------------------------------- ! ! Purpose: ! slt initialization ! !----------------------------------------------------------------------- use pmgrid, only: plev, plevp use commap, only: t0, clat use rgrid, only: nlon use physconst, only: cappa, gravit use prognostics, only: u3, v3, n3, div, dpsl, dpsm, ps, ed1 use time_manager, only: is_first_step use sld_control_mod, only: pdela ! ! Arguments ! real(r8), intent(in) :: etamid (plev) ! vertical coords at midpoints real(r8), intent(out) :: gw (plat) ! Gaussian weights real(r8), intent(out) :: detam (plev) ! intervals between vert full levs. real(r8), intent(out) :: cwava (plat) ! weight applied to global integrals type(advection_state), intent(out) :: adv_state ! Advection state ! ! Local variables ! real(r8) :: coslat(plon) ! cosine of latitude real(r8) :: rcoslat(plon) ! 1 over cosine of latitude real(r8) :: hyad (plev) ! del (A) integer :: i, j, k, kk, l, c, lat ! Indices real(r8) :: tmp1 ! temp variable real(r8) :: pdel(plon,plev) ! Pressure difference over layer real(r8) :: pint(plon,plevp) ! Pressure at interfaces real(r8) :: pmid(plon,plev) ! Pressure at layer mid-points !----------------------------------------------------------------------- #if ( defined SPMD ) call slt_spmd() #endif call alloc_advstate( adv_state ) ! ! Eta coordinates on interfaces: Used for calculation etadot vertical velocity ! do k=1,plevp etaint(k) = hyai(k) + hybi(k) end do ! ! Initialize matrix gamma (used in sld T computation) ! gamma(:,:) = 0._r8 do k = 1,plev tmp1 = cappa*t0(k)*hypi(plevp)/hypm(k) gamma(k,k) = 0.5_r8*tmp1 do l=1,k-1 gamma(l,k) = tmp1 end do end do ! ! Set slt common block variables ! call grdini(pmap ,etamid ,etaint ,gravit ,dlam , & lam ,phi ,dphi ,gw ,sinlam , & coslam ,lbasdy ,lbasdz ,lbassd ,lbasiy , & lbasiz ,lbassi ,detam ,detai ,kdpmpf , & kdpmph ,cwava ,phigs ) ! ! Compute pdel from "A" portion of hybrid vertical grid ! do k=1,plev hyad(k) = hyai(k+1) - hyai(k) end do do k=1,plev do i=1,plon pdela(i,k) = hyad(k)*ps0 end do end do if (is_first_step()) then do lat=beglat,endlat do i=1,nlon(lat) coslat(i) = cos(clat(lat)) rcoslat(i) = 1._r8/coslat(i) end do ! ! Set current time pressure arrays for model levels etc. ! call plevs0(nlon(lat), plon, plev, ps(1,lat,n3), pint, pmid, pdel) ! ! Compute appropriate (1/ps)etadot(dp/deta) ! call etadt0 (lat, nlon(lat), & rcoslat ,div(1,1,lat,n3), u3(1,1,lat,n3), v3(1,1,lat,n3), dpsl(1,lat), & dpsm(1,lat), pdel, ps(1,lat,n3), ed1(1,1,lat)) end do end if ! !---------------------------------------------------------- ! Initialize departure point bin arrays to 0. !---------------------------------------------------------- ! kstep = 0 do lat = beglat,endlat do k = 1,plev do kk = 1,plev levkntl(kk,k,lat) = 0._r8 end do end do end do end subroutine slt_initial ! !----------------------------------------------------------------------- ! subroutine slt_run( ztodt ,detam ,cwava ,etamid , & grlps1 ,grlps2 ,grt1 ,grt2 ,grq1 , & grq2 ,grfu1 ,grfu2 ,grfv1 ,grfv2 , & grfu ,grfv ,lnpssld ,prhssld ,adv_state ) !----------------------------------------------------------------------- ! ! Purpose: ! slt run method ! !----------------------------------------------------------------------- use pmgrid, only: plat, plon use comspe, only: maxm use constituents, only: pcnst use rgrid, only: nlon use prognostics, only: parrsld, tarrsld, u3, v3, etadot, n3 #if ( defined SPMD ) use mpishorthand, only: mpicom #endif ! !------------------------------Arguments-------------------------------- ! real(r8), intent(in) :: ztodt ! twice the timestep unless nstep=0 real(r8), intent(in) :: detam (plev) ! delta eta at levels real(r8), intent(in) :: cwava (plat) ! weight for global water vapor int. real(r8), intent(in) :: etamid (plev) ! eta at levels real(r8), intent(out) :: grlps1(2*maxm,plat/2) ! ------------------------------ real(r8), intent(out) :: grlps2(2*maxm,plat/2) ! | real(r8), intent(out) :: grt1 (2*maxm,plev,plat/2) ! | real(r8), intent(out) :: grt2 (2*maxm,plev,plat/2) ! | real(r8), intent(out) :: grq1 (2*maxm,plev,plat/2) ! | real(r8), intent(out) :: grq2 (2*maxm,plev,plat/2) ! |- see quad for definitions real(r8), intent(out) :: grfu1 (2*maxm,plev,plat/2) ! | real(r8), intent(out) :: grfu2 (2*maxm,plev,plat/2) ! | real(r8), intent(out) :: grfv1 (2*maxm,plev,plat/2) ! | real(r8), intent(out) :: grfv2 (2*maxm,plev,plat/2) ! ------------------------------ real(r8), intent(inout) :: grfu (plon,plev,beglat:endlat) ! nonlinear term - u momentum eqn real(r8), intent(inout) :: grfv (plon,plev,beglat:endlat) ! nonlinear term - v momentum eqn real(r8), intent(in) :: lnpssld(plon,plev,beglat:endlat) real(r8), intent(in) :: prhssld(plon,plev,beglat:endlat) type(advection_state), intent(inout) :: adv_state !---------------------------Local workspace----------------------------- ! integer, parameter :: iter=1 ! number of iterations to be used in dep pt calc. integer :: c ! constituent index integer :: irow ! latitude pair index integer :: lat,latn,lats ! latitude indices integer :: nlon_fft_in ! FFT work array inner dimension integer :: nlon_fft_out ! FFT work array inner dimension real(r8) :: onepeps ! 1 + epssld real(r8) :: dtr ! 1/dt ! ! FFT buffers ! real(r8), allocatable:: fftbuf_in(:,:,:,:) ! fftbuf_in(nlon_fft_in,plev,5,beglat:endlat) real(r8), allocatable:: fftbuf_out(:,:,:,:) ! fftbuf_out(nlon_fft_out,plev,5,plat) call t_startf ('da_coupling') call da_coupling( adv_state, cwava, lnpssld, prhssld ) call t_stopf ('da_coupling') #if ( defined SPMD ) ! ! Communicate boundary information ! call t_barrierf ('sync_bndexch', mpicom) call t_startf ('bndexch') call bndexch( adv_state ) call t_stopf ('bndexch') #endif onepeps = 1._r8 + epssld dtr = 2._r8/(ztodt*onepeps) ! ! Fill latitude/longitude extensions of constituents and dynamics terms ! if (beglatex .le. endlatex) then !$OMP PARALLEL DO PRIVATE (C) do c = 1, pcnst call extys(1 ,plev ,adv_state%q3(1,1,beglatex,c), 1) end do call extyv(1 ,plev ,coslam ,sinlam ,adv_state%u3sld (1,1,beglatex ), & adv_state%v3sld (1,1,beglatex )) call extys(1 ,plevp ,adv_state%etadot (1,1,beglatex), 1) ! call extyv(1 ,plev ,coslam ,sinlam ,adv_state%u3 (1,1,beglatex), & adv_state%v3 (1,1,beglatex)) call extys(1 ,plev ,adv_state%t3 (1,1,beglatex), 1) call extys(1 ,plev ,adv_state%lnpssld(1,1,beglatex ), 1) call extys(1 ,plev ,adv_state%prhssld(1,1,beglatex ), 1) ! !$OMP PARALLEL DO PRIVATE (C) do C = 1, pcnst call extx (1 ,plev ,adv_state%q3 (1,1,beglatex,c), 1) end do call extx (1 ,plev ,adv_state%u3sld (1,1,beglatex ), 1) call extx (1 ,plev ,adv_state%v3sld (1,1,beglatex ), 1) call extx (1 ,plevp ,adv_state%etadot (1,1,beglatex), 1) ! call extx (1 ,plev ,adv_state%u3 (1,1,beglatex), 1) call extx (1 ,plev ,adv_state%v3 (1,1,beglatex), 1) call extx (1 ,plev ,adv_state%t3 (1,1,beglatex), 1) call extx (1 ,plev ,adv_state%lnpssld(1,1,beglatex ), 1) call extx (1 ,plev ,adv_state%prhssld(1,1,beglatex ), 1) endif ! ! Begin SLT interpolation ! ! call t_startf ('scanslt_alloc') nlon_fft_in = plondfft allocate(fftbuf_in(nlon_fft_in,plev,5,beglat:endlat)) #if ( defined SPMD ) #ifdef NEC_SX nlon_fft_out = 2*maxm + 1 #else nlon_fft_out = 2*maxm #endif allocate(fftbuf_out(nlon_fft_out,plev,5,plat)) #else nlon_fft_out = 1 allocate(fftbuf_out(1,1,1,1)) #endif call t_stopf ('scanslt_alloc') call t_startf ('scanslt_bft') !$OMP PARALLEL DO PRIVATE(LAT) do lat = beglat,endlat call scanslt_bft(ztodt ,lat ,dtr ,iter , & detam ,etamid ,grfu ,grfv , & tarrsld ,parrsld ,u3(1,1,lat,n3) , & v3(1,1,lat,n3) ,etadot(1,1,lat) , & nlon(lat) ,nlon_fft_in , & fftbuf_in(1,1,1,lat), adv_state ) end do ! end lat-loop call t_stopf ('scanslt_bft') ! call t_startf ('scanslt_fft') call scanslt_fft(nlon_fft_in,nlon_fft_out,fftbuf_in,fftbuf_out) call t_stopf ('scanslt_fft') ! call t_startf ('scanslt_aft') !$OMP PARALLEL DO PRIVATE (IROW, LATN, LATS) do irow=1,plat/2 lats = irow latn = plat - irow + 1 #if ( defined SPMD ) call scanslt_aft(irow, nlon_fft_out, & fftbuf_out(1,1,1,lats), fftbuf_out(1,1,1,latn), & grlps1(1,irow) ,grlps2(1,irow) , & grt1(1,1,irow) ,grt2(1,1,irow) , & grq1(1,1,irow) ,grq2(1,1,irow) , & grfu1(1,1,irow) ,grfu2(1,1,irow), & grfv1(1,1,irow) ,grfv2(1,1,irow) ) #else call scanslt_aft(irow, nlon_fft_in, & fftbuf_in(1,1,1,lats), fftbuf_in(1,1,1,latn), & grlps1(1,irow) ,grlps2(1,irow) , & grt1(1,1,irow) ,grt2(1,1,irow) , & grq1(1,1,irow) ,grq2(1,1,irow) , & grfu1(1,1,irow) ,grfu2(1,1,irow), & grfv1(1,1,irow) ,grfv2(1,1,irow) ) #endif end do ! end irow-loop call t_stopf ('scanslt_aft') call t_startf ('scanslt_dealloc') deallocate(fftbuf_in) deallocate(fftbuf_out) call t_stopf ('scanslt_dealloc') call t_startf ('ad_coupling') call ad_coupling( adv_state ) call t_stopf ('ad_coupling') end subroutine slt_run ! !----------------------------------------------------------------------- ! subroutine alloc_advstate( adv_state) !----------------------------------------------------------------------- ! ! Purpose: ! Allocate advection_state structure ! !----------------------------------------------------------------------- use infnan, only: nan type(advection_state), intent(out) :: adv_state character(len=*), parameter :: sub = 'alloc_advstate' !pw if ( (beglatex == 0) .or. (endlatex == 0) ) then !pw call endrun(sub//': extended grid begining and ending latitude have NOT been set yet!') !pw end if allocate (adv_state%u3 (plond,plev, beglatex:endlatex) ) allocate (adv_state%v3 (plond,plev, beglatex:endlatex) ) allocate (adv_state%t3 (plond,plev, beglatex:endlatex) ) allocate (adv_state%q3 (plond,plev, beglatex:endlatex,pcnst)) allocate (adv_state%u3sld (plond,plev, beglatex:endlatex) ) allocate (adv_state%v3sld (plond,plev, beglatex:endlatex) ) allocate (adv_state%lnpssld(plond,plev, beglatex:endlatex) ) allocate (adv_state%prhssld(plond,plev, beglatex:endlatex) ) allocate (adv_state%etadot (plond,plevp,beglatex:endlatex) ) if (beglatex .le. endlatex) then adv_state%u3(:,:,:) = nan adv_state%v3(:,:,:) = nan adv_state%t3(:,:,:) = nan adv_state%q3(:,:,:,:) = nan adv_state%u3sld(:,:,:) = nan adv_state%v3sld(:,:,:) = nan adv_state%lnpssld(:,:,:) = nan adv_state%prhssld(:,:,:) = nan adv_state%etadot(:,:,:) = nan end if end subroutine alloc_advstate ! !----------------------------------------------------------------------- ! subroutine dealloc_advstate( adv_state ) !----------------------------------------------------------------------- ! ! Purpose: ! De-allocate advection_state structure ! !----------------------------------------------------------------------- type(advection_state), intent(inout) :: adv_state deallocate (adv_state%u3) deallocate (adv_state%v3) deallocate (adv_state%t3) deallocate (adv_state%q3) deallocate (adv_state%u3sld) deallocate (adv_state%v3sld) deallocate (adv_state%lnpssld) deallocate (adv_state%prhssld) deallocate (adv_state%etadot) end subroutine dealloc_advstate ! !----------------------------------------------------------------------- ! subroutine da_coupling ( adv_state, cwava, lnpssld, prhssld ) !----------------------------------------------------------------------- ! ! Purpose: ! Dynamics to advection coupling ! !----------------------------------------------------------------------- use prognostics, only: u3, v3, t3, q3, etadot, n3, n3m1, ps use time_manager, only: is_first_step use rgrid, only: nlon use commap, only: w use qmassa, only: qmassarun type(advection_state), intent(inout) :: adv_state real(r8), intent(in) :: cwava (plat) ! weight applied to global integrals real(r8), intent(in) :: lnpssld(plon,plev,beglat:endlat) real(r8), intent(in) :: prhssld(plon,plev,beglat:endlat) real(r8) pmid (plon,plev) ! pressure at model levels real(r8) pint (plon,plevp) ! pressure at interfaces real(r8) pdel (plon,plev) ! pressure difference between integer :: i, j, k, c, jcen, icen ! Indices integer :: irow ! latitude pair index if(is_first_step()) then !$OMP PARALLEL DO PRIVATE (I,K,J,JCEN,ICEN) do j = beglat, endlat jcen = j1 - 1 + j do k = 1, plevp do i = 1, nlon(j) icen = i1 + i - 1 adv_state%etadot(icen,k,jcen) = etadot(i,k,j) end do end do end do end if !$OMP PARALLEL DO PRIVATE (I,K,J,ICEN,JCEN,IROW,PINT,PMID,PDEL) do j = beglat, endlat jcen = j1 - 1 + j if( j <= plat/2 ) then irow = j else irow = plat + 1 - j end if do k = 1, plev do i = 1, nlon(j) icen = i1 + i - 1 adv_state%u3(icen,k,jcen) = u3(i,k,j,n3m1) adv_state%v3(icen,k,jcen) = v3(i,k,j,n3m1) adv_state%t3(icen,k,jcen) = t3(i,k,j,n3m1) adv_state%lnpssld(icen,k,jcen) = lnpssld(i,k,j) adv_state%prhssld(icen,k,jcen) = prhssld(i,k,j) end do end do ! ! Calculate mass of moisture in field being advected by slt. ! call plevs0(nlon(j), plon, plev, ps(1,j,n3), pint ,pmid ,pdel) call qmassarun(cwava(j) , w(irow) , q3(1,1,1,j,n3), pdel ,& hw1lat(1,j), nlon(j), q3(1,1,1,j,n3m1), j) ! ! The modified etadot @n3m1 will be used later for trajectory calculation in SCANSLT ! do k = 1, plevp do i = 1, nlon(j) icen = i1 + i - 1 adv_state%etadot(icen,k,jcen) = 2._r8*etadot(i,k,j) - & adv_state%etadot(icen,k,jcen) end do end do end do !$OMP PARALLEL DO PRIVATE (I,K,J,C,ICEN,JCEN) do c = 1, pcnst do j = beglat, endlat jcen = j1 - 1 + j do k = 1, plev do i = 1, nlon(j) icen = i1 + i - 1 adv_state%q3(icen,k,jcen,c) = q3(i,k,c,j,n3) end do end do end do end do end subroutine da_coupling ! !----------------------------------------------------------------------- ! subroutine slt_run_setup( adv_state ) !----------------------------------------------------------------------- ! ! Purpose: ! Setup internal advection state data for running scanslt. ! !----------------------------------------------------------------------- use prognostics, only: u3, v3, n3, n3m1, etadot use rgrid, only: nlon use time_manager, only: is_first_restart_step type(advection_state), intent(out) :: adv_state integer :: i, j, k, icen, jcen ! Indices ! ! Compute U/V for time n + 1/2 ! This will be used later for trajectory calculation ! !$OMP PARALLEL DO PRIVATE (I,K,J,JCEN,ICEN) do j = beglat, endlat jcen = j1 - 1 + j do k = 1, plev do i = 1, nlon(j) icen = i1 + i - 1 adv_state%u3sld(icen,k,jcen) = 2._r8*u3(i,k,j,n3) - u3(i,k,j,n3m1) adv_state%v3sld(icen,k,jcen) = 2._r8*v3(i,k,j,n3) - v3(i,k,j,n3m1) end do end do end do if(is_first_restart_step()) then !$OMP PARALLEL DO PRIVATE (I,K,J,ICEN,JCEN) do j = beglat, endlat jcen = j1 - 1 + j do k = 1, plevp do i = 1, nlon(j) icen = i1 + i - 1 adv_state%etadot(icen,k,jcen) = etadot(i,k,j) end do end do end do end if end subroutine slt_run_setup ! !----------------------------------------------------------------------- ! subroutine ad_coupling ( adv_state ) !----------------------------------------------------------------------- ! ! Purpose: ! Advection to dynamics coupling ! !----------------------------------------------------------------------- use prognostics, only: u3, v3, t3, q3, n3, n3m1, etadot use rgrid, only: nlon type(advection_state), intent(inout) :: adv_state integer :: i, j, k, c, icen, jcen !$OMP PARALLEL DO PRIVATE (I,K,J,ICEN,JCEN) do j = beglat, endlat jcen = j1 - 1 + j do k = 1, plev do i = 1, nlon(j) icen = i1 + i - 1 u3(i,k,j,n3m1) = adv_state%u3(icen,k,jcen) v3(i,k,j,n3m1) = adv_state%v3(icen,k,jcen) t3(i,k,j,n3m1) = adv_state%t3(icen,k,jcen) end do end do do k = 1, plevp do i = 1, nlon(j) icen = i1 + i - 1 adv_state%etadot(icen,k,jcen) = etadot(i,k,j) end do end do end do !$OMP PARALLEL DO PRIVATE (I,K,J,C,ICEN,JCEN) do c = 1, pcnst do j = beglat, endlat jcen = j1 - 1 + j do k = 1, plev do i = 1, nlon(j) icen = i1 + i - 1 q3(i,k,c,j,n3) = adv_state%q3(icen,k,jcen,c) end do end do end do end do ! ! Increment kstep ! kstep = kstep + 1 end subroutine ad_coupling ! !----------------------------------------------------------------------- ! subroutine scanslt_bft (ztodt ,lat ,dtr ,iter ,detam , & etamid ,grfu ,grfv ,tarrsld ,parrsld , & u3 ,v3 ,etadot ,nlon ,nlon_fft, & fftbuf ,adv_state ) !----------------------------------------------------------------------- ! ! Purpose: ! Interpolate terms for semi-lagrangian transport and SLD dynamics. ! One latitude slice only ! ! Author: J. Olson ! !----------------------------------------------------------------------- use rgrid, only: nmmax use pspect, only: pmmax use commap, only: w, clat, t0 use physconst, only: cappa, ra !------------------------------Arguments-------------------------------- ! real(r8), intent(in) :: ztodt ! twice the time step unless nstep = 0 integer , intent(in) :: lat ! latitude index real(r8), intent(in) :: dtr ! 1/dt integer , intent(in) :: iter ! number of iterations for trajectory real(r8), intent(in) :: detam (plev) ! delta eta at levels real(r8), intent(in) :: etamid (plev) ! eta at levels real(r8), intent(in) :: grfu (plon,plev,beglat:endlat) ! nonlinear term - u momentum eqn real(r8), intent(in) :: grfv (plon,plev,beglat:endlat) ! nonlinear term - v momentum eqn real(r8), intent(in) :: tarrsld (plon,plev,beglat:endlat) ! T at arr. pt.(SLD) real(r8), intent(inout):: parrsld (plon,plev,beglat:endlat) ! Ps at arr. pt.(SLD) real(r8), intent(in) :: u3 (plon,plev) ! u-wind at current time real(r8), intent(in) :: v3 (plon,plev) ! v-wind at current time real(r8), intent(in) :: etadot (plon,plevp) ! vertical motion at current time integer , intent(in) :: nlon ! # of longitudes integer, intent(in) :: nlon_fft ! first dimension of FFT work array real(r8), intent(out) :: fftbuf(nlon_fft,plev,5) ! buffer used for in-place FFTs type(advection_state), intent(in) :: adv_state ! !---------------------------Local workspace----------------------------- ! integer i ! index integer k ! index integer l ! index integer m ! constituent index integer irow ! N/S latitude pair index integer jcen ! lat index (extended grid) ! ! of forecast real(r8) qtmp (plond,plev,beglatex:endlatex) ! Temporary for q array real(r8) fdp (plon,plev,2) ! interpolant real(r8) lamdp(plon,plev) ! x-coord of dep pt real(r8) phidp(plon,plev) ! y-coord of dep pt real(r8) sigdp(plon,plev) ! z-coord of dep pt integer idp (plon,plev,4) ! zonal dep point index integer jdp (plon,plev) ! meridional dep point index integer kdp (plon,plev) ! vertical dep point index integer kkdp (plon,plev) ! index of z-coordinate of dep pt (alt) real(r8) xl (plon,plev,4) ! weight for x-interpolants (left) real(r8) xr (plon,plev,4) ! weight for x-interpolants (right) real(r8) wgt1x(plon,plev,4) ! weight for x-interpolants (Lag Cubic) real(r8) wgt2x(plon,plev,4) ! weight for x-interpolants (Lag Cubic) real(r8) wgt3x(plon,plev,4) ! weight for x-interpolants (Lag Cubic) real(r8) wgt4x(plon,plev,4) ! weight for x-interpolants (Lag Cubic) real(r8) hl (plon,plev,4) ! weight for x-interpolants (Hermite) real(r8) hr (plon,plev,4) ! weight for x-interpolants (Hermite) real(r8) dhl (plon,plev,4) ! weight for x-interpolants (Hermite) real(r8) dhr (plon,plev,4) ! weight for x-interpolants (Hermite) real(r8) ys (plon,plev) ! weight for y-interpolants (south) real(r8) yn (plon,plev) ! weight for y-interpolants (north) real(r8) wgt1y(plon,plev) ! weight for y-interpolants (Lag Cubic) real(r8) wgt2y(plon,plev) ! weight for y-interpolants (Lag Cubic) real(r8) wgt3y(plon,plev) ! weight for y-interpolants (Lag Cubic) real(r8) wgt4y(plon,plev) ! weight for y-interpolants (Lag Cubic) real(r8) hs (plon,plev) ! weight for y-interpolants (Hermite) real(r8) hn (plon,plev) ! weight for y-interpolants (Hermite) real(r8) dhs (plon,plev) ! weight for y-interpolants (Hermite) real(r8) dhn (plon,plev) ! weight for y-interpolants (Hermite) real(r8) rdphi(plon,plev) ! reciprocal of y-interval real(r8) wgt1z(plon,plev) ! weight for z-interpolants (Lag Cubic) real(r8) wgt2z(plon,plev) ! weight for z-interpolants (Lag Cubic) real(r8) wgt3z(plon,plev) ! weight for z-interpolants (Lag Cubic) real(r8) wgt4z(plon,plev) ! weight for z-interpolants (Lag Cubic) real(r8) hb (plon,plev) ! weight for z-interpolants (Hermite) real(r8) ht (plon,plev) ! weight for z-interpolants (Hermite) real(r8) dhb (plon,plev) ! weight for z-interpolants (Hermite) real(r8) dht (plon,plev) ! weight for z-interpolants (Hermite) real(r8) rdz (plon,plev) ! reciprocal of z-interval real(r8) zt (plon,plev) ! top vertical interpolation weight real(r8) zb (plon,plev) ! bot vertical interpolation weight real(r8) lampr(plon,plev) ! trajectory increment (x-direction) real(r8) phipr(plon,plev) ! trajectory increment (y-direction) real(r8) upr (plon,plev) ! interpolated u field (local geodesic) real(r8) vpr (plon,plev) ! interpolated v field (local geodesic) ! real(r8) pd (plon) ! RHS term for Ps and (1/ps)etadot(dp/deta) real(r8) pdsum(plon) ! RHS term for Ps and (1/ps)etadot(dp/deta) real(r8) pd1 (plon) ! RHS term for Ps and (1/ps)etadot(dp/deta) real(r8) pdsm1(plon) ! RHS term for Ps and (1/ps)etadot(dp/deta) real(r8) pa (plon) ! RHS term for Ps and (1/ps)etadot(dp/deta) real(r8) pasum(plon) ! RHS term for Ps and (1/ps)etadot(dp/deta) real(r8) coslat ! cos(latitude) real(r8) tmp1 ! temp space ! logical limdrh ! horizontal derivative limiter flag logical limdrv ! vertical derivative limiter flag logical lhrzint ! horizontal interp flag logical lvrtint ! vertical interp flag logical lhrzwgt ! flag to compute horizontal weights logical lvrtwgt ! flag to compute vertical weights ! real(r8) grfulat(plon,plev) ! non-linear terms for u-momentum real(r8) grfvlat(plon,plev) ! non-linear terms for u-momentum real(r8) grtlat (plon,plev) ! RHS of T-eqn real(r8) grqlat (plon,plev) ! q real(r8) grpslat(plon) ! RHS of Ps-eqn ! integer, parameter :: fudex = 1 ! indices into fftbuf integer, parameter :: fvdex = 2 integer, parameter :: tdex = 3 integer, parameter :: qdex = 4 integer, parameter :: psdex = 5 ! !----------------------------------------------------------------------- ! if(lat.le.plat/2) then irow = lat else irow = plat + 1 - lat end if jcen = j1 - 1 + lat coslat = cos(clat(lat)) ! ! Initial guess for trajectory midpoints in spherical coords. ! Use arrival points as initial guess for trajectory midpoints. ! do k=1,plev do i=1,nlon phidp(i,k) = clat(lat) sigdp(i,k) = etamid(k) end do end do ! ! Offset bottom level departure point first guess by epsilon ! do i = 1,nlon sigdp(i,plev) = sigdp(i,plev)*(1._r8 - 10._r8*epsilon(sigdp)) end do ! ! Loop through latitudes producing departure point calculation ! call slttraj(pmap ,jcen ,lat ,ztodt ,ra , & iter ,lam ,phi ,dphi ,etamid , & etaint ,detam ,detai ,lbasiy ,lbasiz , & lbassi ,kdpmpf ,kdpmph ,idp ,jdp , & kdp ,kkdp ,xl ,xr ,wgt1x , & wgt2x ,wgt3x ,wgt4x ,hl ,hr , & dhl ,dhr ,ys ,yn ,wgt1y , & wgt2y ,wgt3y ,wgt4y ,hs ,hn , & dhs ,dhn ,rdphi ,wgt1z ,wgt2z , & wgt3z ,wgt4z ,hb ,ht ,dhb , & dht ,rdz ,lampr ,phipr ,upr , & vpr ,lamdp ,phidp ,sigdp ,adv_state%u3 , & adv_state%v3 ,adv_state%u3sld ,adv_state%v3sld , & adv_state%etadot ,u3 ,v3 ,etadot , & dlam ,nlon ) ! ! Compute constituent forecast ! lhrzwgt = .true. lvrtwgt = .true. lhrzint = .true. lvrtint = .true. limdrh = .true. limdrv = .true. call bandij (dlam ,phi ,lamdp ,phidp ,idp , & jdp ,nlon ) call kdpfnd (plev ,pmap ,etamid ,sigdp ,kdpmpf , & kdp ,nlon ) call sltwgts(limdrh ,limdrv ,lhrzwgt ,lvrtwgt ,plev , & idp ,jdp ,kdp ,lam ,phi , & etamid ,dphi ,detam ,lamdp ,phidp , & sigdp ,lbasiy ,lbasiz ,kkdp ,xl , & xr ,wgt1x ,wgt2x ,wgt3x ,wgt4x , & hl ,hr ,dhl ,dhr ,ys , & yn ,wgt1y ,wgt2y ,wgt3y ,wgt4y , & hs ,hn ,dhs ,dhn ,rdphi , & wgt1z ,wgt2z ,wgt3z ,wgt4z ,hb , & ht ,dhb ,dht ,rdz ,zt , & zb ,nlon ) do m = 1,pcnst call sltint (jcen ,adv_state%q3(1,1,beglatex,m),lam ,rdphi , & rdz ,lbasdy ,lbasdz ,xl ,xr , & wgt1x ,wgt2x ,wgt3x ,wgt4x ,hl , & hr ,dhl ,dhr ,ys ,yn , & wgt1y ,wgt2y ,wgt3y ,wgt4y ,hs , & hn ,dhs ,dhn ,wgt1z ,wgt2z , & wgt3z ,wgt4z ,hb ,ht ,dhb , & dht ,idp ,jdp ,kdp ,kkdp , & lhrzint ,lvrtint ,limdrh ,limdrv ,qfcst(1,1,m,lat), & nlon ) if ( m == 1 ) grqlat(:nlon,:plev) = qfcst(:nlon,:plev,m,lat) end do ! ! Accumulate P-interpolants into T equation ! do i = 1,nlon grpslat(i) = 0._r8 pasum (i) = 0._r8 pa (i) = 0._r8 end do do k = 1,plev do i = 1,nlon grtlat(i,k) = tarrsld (i,k,lat) end do end do do k = 1,plev do l = 1,k do i = 1,nlon grtlat(i,k) = grtlat(i,k) + parrsld(i,l,lat)*gamma(l,k) end do end do end do ! ! Accumulate Ps interpolants in 3-D array ! do k = 1,plev do i = 1,nlon grpslat(i) = grpslat(i) + parrsld(i,k,lat) pasum (i) = pasum (i) + parrsld(i,k,lat) end do end do ! ! Compute first part of (1/ps)etadot(dp/deta) ! do k = 1,plev-1 do i = 1,nlon pa(i) = pa(i) + parrsld(i,k,lat) parrsld(i,k,lat) = pa(i) end do ! if(k.ge.nprlev) then do i = 1,nlon parrsld(i,k,lat) = parrsld(i,k,lat) - hybi(k+1)*( pasum(i) ) end do end if end do ! ! Compute U, V interpolants: Non-monotonic, 3-D interpolation ! limdrh = .false. limdrv = .true. lhrzint = .true. lvrtint = .true. call sltint (jcen ,adv_state%u3,lam , & rdphi ,rdz ,lbasdy ,lbasdz ,xl , & xr ,wgt1x ,wgt2x ,wgt3x ,wgt4x , & hl ,hr ,dhl ,dhr ,ys , & yn ,wgt1y ,wgt2y ,wgt3y ,wgt4y , & hs ,hn ,dhs ,dhn ,wgt1z , & wgt2z ,wgt3z ,wgt4z ,hb ,ht , & dhb ,dht ,idp ,jdp ,kdp , & kkdp ,lhrzint ,lvrtint ,limdrh ,limdrv , & fdp ,nlon ) call sltint (jcen ,adv_state%v3,lam , & rdphi ,rdz ,lbasdy ,lbasdz ,xl , & xr ,wgt1x ,wgt2x ,wgt3x ,wgt4x , & hl ,hr ,dhl ,dhr ,ys , & yn ,wgt1y ,wgt2y ,wgt3y ,wgt4y , & hs ,hn ,dhs ,dhn ,wgt1z , & wgt2z ,wgt3z ,wgt4z ,hb ,ht , & dhb ,dht ,idp ,jdp ,kdp , & kkdp ,lhrzint ,lvrtint ,limdrh ,limdrv , & fdp(1,1,2),nlon ) ! ! Evaluate last half of grfu and grfv (Nu,Nv) ! call nunv1(lam(i1,jcen) ,phi(jcen),lamdp ,phidp ,fdp(1,1,1), & fdp(1,1,2) ,coslat ,grfu(1,1,lat),grfv(1,1,lat) ,grfulat , & grfvlat ,nlon) ! ! Compute T interpolants: Non-monotonic, 3-D interpolation ! limdrh = .false. limdrv = .true. lhrzint = .true. lvrtint = .true. call sltint (jcen ,adv_state%t3,lam , & rdphi ,rdz ,lbasdy ,lbasdz ,xl , & xr ,wgt1x ,wgt2x ,wgt3x ,wgt4x , & hl ,hr ,dhl ,dhr ,ys , & yn ,wgt1y ,wgt2y ,wgt3y ,wgt4y , & hs ,hn ,dhs ,dhn ,wgt1z , & wgt2z ,wgt3z ,wgt4z ,hb ,ht , & dhb ,dht ,idp ,jdp ,kdp , & kkdp ,lhrzint ,lvrtint ,limdrh ,limdrv , & fdp(1,1,1),nlon ) ! ! Accumulate T interpolants in 3-D array ! do k = 1,plev do i = 1,nlon grtlat(i,k) = grtlat(i,k) + fdp(i,k,1) end do end do ! ! Reset "kdp" to arrival indices everywhere so that we can do true ! horizontal interpolation rather than "vertical non-interpolation" ! do k = 1,plev do i = 1,nlon kdp(i,k) = k end do end do ! ! Compute Ps and remaining T interpolants: ! Non-monotonic, 2-D interpolation ! limdrh = .false. limdrv = .false. lhrzint = .true. lvrtint = .false. call sltint (jcen ,adv_state%lnpssld ,lam ,rdphi , & rdz ,lbasdy ,lbasdz ,xl ,xr , & wgt1x ,wgt2x ,wgt3x ,wgt4x ,hl , & hr ,dhl ,dhr ,ys ,yn , & wgt1y ,wgt2y ,wgt3y ,wgt4y ,hs , & hn ,dhs ,dhn ,wgt1z ,wgt2z , & wgt3z ,wgt4z ,hb ,ht ,dhb , & dht ,idp ,jdp ,kdp ,kkdp , & lhrzint ,lvrtint ,limdrh ,limdrv ,fdp(1,1,1), & nlon ) call sltint (jcen ,adv_state%prhssld ,lam ,rdphi , & rdz ,lbasdy ,lbasdz ,xl ,xr , & wgt1x ,wgt2x ,wgt3x ,wgt4x ,hl , & hr ,dhl ,dhr ,ys ,yn , & wgt1y ,wgt2y ,wgt3y ,wgt4y ,hs , & hn ,dhs ,dhn ,wgt1z ,wgt2z , & wgt3z ,wgt4z ,hb ,ht ,dhb , & dht ,idp ,jdp ,kdp ,kkdp , & lhrzint ,lvrtint ,limdrh ,limdrv ,fdp(1,1,2), & nlon ) do i = 1,nlon pdsum(i) = 0._r8 pd (i) = 0._r8 pdsm1(i) = 0._r8 pd1 (i) = 0._r8 end do ! ! Accumulate P-interpolants into T equation ! do k = nprlev,plev tmp1 = cappa*t0(k)*hypi(plevp)/hypm(k) do i = 1,nlon grtlat(i,k) = grtlat(i,k) - fdp(i,k,1)*tmp1*hybm(k) end do end do do k = nprlev,plev do l = nprlev,k do i = 1,nlon grtlat(i,k) = grtlat(i,k) + fdp(i,l,1)*hybd(l)*gamma(l,k) end do end do end do do k = 1,plev do l = 1,k do i = 1,nlon grtlat(i,k) = grtlat(i,k) + fdp(i,l,2)*gamma(l,k) end do end do end do ! ! Accumulate Ps interpolants in 3-D array ! do k = 1,plev do i = 1,nlon grpslat(i) = grpslat(i) + fdp(i,k,2) pdsum(i) = pdsum(i) + fdp(i,k,2) end do end do do k = nprlev,plev do i = 1,nlon grpslat(i) = grpslat(i) + fdp(i,k,1)*hybd(k) pdsm1(i) = pdsm1(i) + fdp(i,k,1)*hybd(k) end do end do ! ! Compute remainder of (1/ps)etadot(dp/deta) ! do k = 1,plev-1 do i = 1,nlon pd (i) = pd (i) + fdp(i,k,2) parrsld(i,k,lat) = parrsld(i,k,lat) + pd (i) end do ! if(k.ge.nprlev) then do i=1,nlon pd1(i) = pd1(i) + fdp(i,k,1)*hybd(k) parrsld(i,k,lat) = parrsld(i,k,lat) + pd1(i) - & hybi(k+1)*(pdsum(i) + pdsm1(i)) end do end if do i = 1,nlon parrsld(i,k,lat) = parrsld(i,k,lat)*dtr end do end do ! ! Copy fu,fv,T,q,Ps into FFT buffer ! do k = 1,plev do i = 1,nlon fftbuf(i,k,fudex) = grfulat(i,k) fftbuf(i,k,fvdex) = grfvlat(i,k) fftbuf(i,k,tdex) = grtlat(i,k) fftbuf(i,k,qdex) = grqlat(i,k) enddo enddo do i = 1,nlon fftbuf(i,1,psdex) = grpslat(i) enddo return end subroutine scanslt_bft ! !----------------------------------------------------------------------- ! subroutine scanslt_fft (nlon_fft,nlon_fft2,fftbuf,fftbuf2) !----------------------------------------------------------------------- ! ! Purpose: ! Compute FFT of non-linear dynamical terms ! in preparation for Fourier -> spectral quadrature. ! ! Author: J. Olson ! Modified: P. Worley ! !----------------------------------------------------------------------- use pmgrid, only: plon, plat use rgrid, only: nlon #if (defined SPMD) use mpishorthand, only: mpicom use comspe #endif use sld_control_mod, only : pcray, trig, ifax ! ! Input arguments ! integer, intent(in) :: nlon_fft ! first dimension of first FFT work array integer, intent(in) :: nlon_fft2 ! first dimension of second FFT work array ! ! Input/Output arguments ! real(r8), intent(inout) :: fftbuf(nlon_fft,plev,5,beglat:endlat) ! buffer used for in-place FFTs ! ! Output arguments ! #if (defined SPMD) real(r8), intent(out) :: fftbuf2(nlon_fft2,plev,5,plat) ! buffer for returning reorderd Fourier coefficients #else real(r8), intent(in) :: fftbuf2(1) ! buffer unused #endif ! !---------------------------Local workspace----------------------------- ! ! The "work" array has a different size requirement depending upon whether ! the proprietary Cray assembly language version of the FFT library ! routines, or the all-Fortran version, is being used. ! #if ( ! defined USEFFTLIB ) real(r8) work((plon+1)*5*plev) #else real(r8) work((plon+1)*pcray) ! workspace array for fft991 #endif integer lat ! latitude index integer inc ! increment for fft991 integer isign ! flag indicates transform direction integer ntr ! number of transforms to perform ! inc = 1 isign = -1 ntr = 4*plev + 1 !$OMP PARALLEL DO PRIVATE (LAT,WORK) do lat=beglat,endlat fftbuf(nlon(lat)+1:nlon_fft,:,:,lat) = 0.0_r8 call fft991(fftbuf(1,1,1,lat) ,work ,trig(1,lat),ifax(1,lat),inc ,& nlon_fft ,nlon(lat) ,ntr ,isign ) enddo ! #if ( defined SPMD ) ! ! reorder Fourier coefficients ! call t_barrierf ('sync_realloc4a', mpicom) call t_startf('realloc4a') call realloc4a(nlon_fft, nlon_fft2, fftbuf, fftbuf2) call t_stopf('realloc4a') #endif return end subroutine scanslt_fft ! !----------------------------------------------------------------------- ! subroutine scanslt_aft (irow ,nlon_fft,fftbufs ,fftbufn , & grlps1 ,grlps2 ,grt1 ,grt2 , & grq1 ,grq2 ,grfu1 ,grfu2 , & grfv1 ,grfv2 ) !----------------------------------------------------------------------- ! ! Purpose: ! Combine terms in preparation for Fourier -> spectral quadrature. ! ! Author: J. Olson ! Modified: P. Worley ! !----------------------------------------------------------------------- use spmd_utils, only: iam use pspect, only: pmmax #if (defined SPMD) use comspe, only: numm, maxm #else use comspe, only: maxm use rgrid, only: nmmax #endif ! ! Input arguments ! integer, intent(in) :: irow ! latitude pair index integer, intent(in) :: nlon_fft ! first dimension of FFT work arrays real(r8), intent(in) :: fftbufs(nlon_fft,plev,5) ! southern latitude Fourier coefficients real(r8), intent(in) :: fftbufn(nlon_fft,plev,5) ! northern latitude Fourier coefficients real(r8), intent(out) :: grlps1(2*maxm) ! ------------------------------ real(r8), intent(out) :: grlps2(2*maxm) ! | real(r8), intent(out) :: grt1 (2*maxm,plev) ! | real(r8), intent(out) :: grt2 (2*maxm,plev) ! | real(r8), intent(out) :: grq1 (2*maxm,plev) ! |- see quad for definitions real(r8), intent(out) :: grq2 (2*maxm,plev) ! | real(r8), intent(out) :: grfu1 (2*maxm,plev) ! | real(r8), intent(out) :: grfu2 (2*maxm,plev) ! | real(r8), intent(out) :: grfv1 (2*maxm,plev) ! | real(r8), intent(out) :: grfv2 (2*maxm,plev) ! ------------------------------ ! !---------------------------Local workspace----------------------------- ! integer i ! index integer k ! index integer mlength, mstrt ! number of wavenumbers and index offset ! integer, parameter :: fudex = 1 ! indices into fftbuf integer, parameter :: fvdex = 2 integer, parameter :: tdex = 3 integer, parameter :: qdex = 4 integer, parameter :: psdex = 5 ! !----------------------------------------------------------------------- ! #if (defined SPMD) mlength = numm(iam) #else mlength = nmmax(irow) #endif do k = 1,plev do i = 1,2*mlength grfu1(i,k) = 0.5_r8*(fftbufn(i,k,fudex)+fftbufs(i,k,fudex)) grfu2(i,k) = 0.5_r8*(fftbufn(i,k,fudex)-fftbufs(i,k,fudex)) grfv1(i,k) = 0.5_r8*(fftbufn(i,k,fvdex)+fftbufs(i,k,fvdex)) grfv2(i,k) = 0.5_r8*(fftbufn(i,k,fvdex)-fftbufs(i,k,fvdex)) grt1 (i,k) = 0.5_r8*(fftbufn(i,k,tdex)+fftbufs(i,k,tdex)) grt2 (i,k) = 0.5_r8*(fftbufn(i,k,tdex)-fftbufs(i,k,tdex)) grq1 (i,k) = 0.5_r8*(fftbufn(i,k,qdex)+fftbufs(i,k,qdex)) grq2 (i,k) = 0.5_r8*(fftbufn(i,k,qdex)-fftbufs(i,k,qdex)) end do end do do i = 1,2*mlength grlps1(i) = 0.5_r8*(fftbufn(i,1,psdex)+fftbufs(i,1,psdex)) grlps2(i) = 0.5_r8*(fftbufn(i,1,psdex)-fftbufs(i,1,psdex)) end do ! return end subroutine scanslt_aft ! !----------------------------------------------------------------------- ! subroutine slt_final( adv_state ) !----------------------------------------------------------------------- ! ! Purpose: ! slt finalization ! !----------------------------------------------------------------------- use rgrid, only: nlon use spmd_utils, only: masterproc, npes #if ( defined SPMD ) use mpishorthand, only: mpicom, mpir8 use spmd_dyn, only: compute_gsfactors #endif type(advection_state), intent(inout) :: adv_state ! ! Local variables ! integer :: i, lat, k ! | integer :: kk ! | integer :: kk1 ! | integer :: kk2 ! | - indices integer :: l ! | integer :: lvsum ! counter integer :: lsum ! counter real(r8) :: vsum ! accumulator for SLD binning statistics logical :: vflag ! logical flag to indicate that binning has been done ! ! correctly #ifdef SPMD integer :: numsend ! number of items to be sent integer :: numrecv(0:npes-1) ! number of items to be received integer :: displs(0:npes-1) ! displacement array integer :: numperlat ! number of items per latitude band #endif !----------------------------------------------------------------------- call dealloc_advstate( adv_state ) ! ! Compute some statistical info ! #if ( defined SPMD ) numperlat = plev*plev call compute_gsfactors (numperlat, numsend, numrecv, displs) call mpiallgatherv ( levkntl(:,:,beglat:endlat), numsend, mpir8, & levkntl(:,:,:), numrecv, displs, mpir8, mpicom) #endif if (masterproc) then do k = 1,plev do kk = 1,plev levknt(kk,k) = 0._r8 do lat = 1,plat levknt(kk,k) = levknt(kk,k) + levkntl(kk,k,lat) end do end do end do lsum = 0 do lat = 1,plat lsum = lsum + nlon(lat) end do vflag = .false. do k = 1,plev vsum = 0._r8 do kk = 1,plev vsum = vsum + levknt(kk,k) end do lvsum = vsum + 0.01_r8 if( lvsum .ne. lsum*kstep ) vflag = .true. end do do k = 1,plev do kk = 1,plev levknt(kk,k) = levknt(kk,k)/real(lsum*kstep,r8) if (levknt(kk,k) .eq. 0._r8) levknt(kk,k) = 1.e+36_r8 end do end do if(vflag) then write(iulog,*) '********************************************' write(iulog,1000) write(iulog,*) '********************************************' else write(iulog,2000) k = 1 write(iulog,3001) hypm(k)/100._r8,(levknt(kk,k),kk = 1, 6) k = 2 write(iulog,3002) hypm(k)/100._r8,(levknt(kk,k),kk = 1, 7) k = 3 write(iulog,3003) hypm(k)/100._r8,(levknt(kk,k),kk = 1, 8) k = 4 write(iulog,3004) hypm(k)/100._r8,(levknt(kk,k),kk = 1, 9) k = 5 write(iulog,3005) hypm(k)/100._r8,(levknt(kk,k),kk = 1,10) do k = 6,plev kk1 = k-5 kk2 = k+5 if(kk2 .gt. plev) kk2 = plev write(iulog,3000) hypm(k)/100._r8,(levknt(kk,k),kk = kk1,kk2) end do end if end if 1000 format(' ERROR in binning departure points; sums were done'// & ' incorrectly. not printing results.') 2000 format(/40x,' BINNING STATISTICS FOR VERTICAL DEPARTURE POINTS'// & ' level -5 -4 -3 -2 ', & ' -1 0 (arr pt) 1 2 3 4 ', & ' 5'/) 3000 format(' ',f9.4, 4x,11f10.5) 3001 format(' ',f9.4,54x,11f10.5) 3002 format(' ',f9.4,44x,11f10.5) 3003 format(' ',f9.4,34x,11f10.5) 3004 format(' ',f9.4,24x,11f10.5) 3005 format(' ',f9.4,14x,11f10.5) end subroutine slt_final ! !----------------------------------------------------------------------- ! subroutine bndexch( adv_state ) !----------------------------------------------------------------------- ! ! Purpose: ! Pack and Exchange initial prognostic information among all the ! processors ! !----------------------------------------------------------------------- #ifdef SPMD use spmd_utils,only: iam use spmd_dyn, only: cut, dyn_npes #endif ! ! Arguments ! type(advection_state), intent(inout) :: adv_state #ifdef SPMD ! ! Local workspace ! integer ns, nn integer inreg ( 2 ) integer outreg( 2 ) integer others,othern ! Other node ! ! Return if number of processors is less than 2 ! if (dyn_npes .lt. 2) return ! ! For each partition (south and north) communicate boundaries ! on each side of partition among however many neighbors necessary ! ! send south, receive north ! ns = 1 nn = 1 do while (ns .le. neighs .or. nn .le. neighn) if (ns .le. neighs) then others = neighs_proc(ns) ! ! Intersection of my cuts and neighbor processor's extended ! cuts tells if this node needs to send data to neighbor ! call intersct(cut(1,iam),cutex(1,others),outreg ) ns = ns + 1 else others = -1 outreg(1) = 0 outreg(2) = 0 end if if (nn .le. neighn) then othern = neighn_proc(nn) ! ! Intersection of neighbor cuts and this node's extended ! cut tells if this node receives data from neighbor ! call intersct(cut(1,othern),cutex(1,iam),inreg ) nn = nn + 1 else othern = -1 inreg(1) = 0 inreg(2) = 0 end if call bndexch_mpi(others ,outreg ,othern ,inreg, adv_state ) end do ! ! send north, receive south ! ns = 1 nn = 1 do while (ns .le. neighs .or. nn .le. neighn) if (nn .le. neighn) then othern = neighn_proc(nn) ! ! Intersection of my cuts and neighbor processor's extended ! cuts tells if this node needs to send data to neighbor ! call intersct(cut(1,iam),cutex(1,othern),outreg ) nn = nn + 1 else othern = -1 outreg(1) = 0 outreg(2) = 0 end if if (ns .le. neighs) then others = neighs_proc(ns) ! ! Intersection of neighbor cuts and this node's extended ! cut tells if this node receives data from neighbor ! call intersct(cut(1,others),cutex(1,iam),inreg ) ns = ns + 1 else others = -1 inreg(1) = 0 inreg(2) = 0 end if call bndexch_mpi(othern ,outreg ,others ,inreg, adv_state ) end do #endif return end subroutine bndexch ! !----------------------------------------------------------------------- ! #ifdef SPMD subroutine bndexch_mpi(othero ,outreg ,otheri ,inreg, adv_state ) !----------------------------------------------------------------------- ! Send initial prognostic information to my peer process !----------------------------------------------------------------------- use mpishorthand !---------------------------- Parameters ------------------------------- integer, parameter :: msgtype = 6000 ! id for message passing integer, parameter :: j1m = j1 - 1 ! lat index just before first "real" model lat !------------------------------Arguments-------------------------------- ! integer , intent(in) :: othero integer , intent(in) :: outreg(2) integer , intent(in) :: otheri integer , intent(in) :: inreg(2) type(advection_state), intent(inout) :: adv_state ! !---------------------------Local workspace----------------------------- ! integer num integer msg integer, parameter :: numvars = 8+pcnst integer reqs(numvars*(plat+1)) integer stats(MPI_STATUS_SIZE, numvars*(plat+1)) integer reqr(numvars*(plat+1)) integer statr(MPI_STATUS_SIZE, numvars*(plat+1)) integer i,j,c integer reqs_i,reqr_i ! !----------------------------------------------------------------------- ! reqs_i = 0 if (othero .ne. -1) then do i = outreg(1), outreg(2) j = numvars*(i-outreg(1)) msg = msgtype + j reqs_i = reqs_i + 1 call mpiisend (adv_state%u3(1,1,j1m+i),plndlv,mpir8, & othero,msg,mpicom,reqs(reqs_i)) msg = msgtype + j + 1 reqs_i = reqs_i + 1 call mpiisend (adv_state%v3(1,1,j1m+i),plndlv,mpir8, & othero,msg,mpicom,reqs(reqs_i)) msg = msgtype + j + 2 reqs_i = reqs_i + 1 call mpiisend (adv_state%t3(1,1,j1m+i),plndlv,mpir8, & othero,msg,mpicom,reqs(reqs_i)) num = plndlv do c = 1, pcnst msg = msgtype + j + 2 + c reqs_i = reqs_i + 1 call mpiisend (adv_state%q3(1,1,j1m+i,c),num,mpir8, & othero,msg,mpicom,reqs(reqs_i)) end do msg = msgtype + j + 3 + pcnst reqs_i = reqs_i + 1 call mpiisend (adv_state%lnpssld(1,1,j1m+i),plndlv,mpir8, & othero,msg,mpicom,reqs(reqs_i)) msg = msgtype + j + 4 + pcnst reqs_i = reqs_i + 1 call mpiisend (adv_state%prhssld(1,1,j1m+i),plndlv,mpir8, & othero,msg,mpicom,reqs(reqs_i)) msg = msgtype + j + 5 + pcnst reqs_i = reqs_i + 1 call mpiisend (adv_state%u3sld(1,1,j1m+i),plndlv,mpir8, & othero,msg,mpicom,reqs(reqs_i)) msg = msgtype + j + 6 + pcnst reqs_i = reqs_i + 1 call mpiisend (adv_state%v3sld(1,1,j1m+i),plndlv,mpir8, & othero,msg,mpicom,reqs(reqs_i)) msg = msgtype + j + 7 + pcnst reqs_i = reqs_i + 1 call mpiisend (adv_state%etadot(1,1,j1m+i),plond*plevp,mpir8, & othero,msg,mpicom,reqs(reqs_i)) end do end if reqr_i = 0 if (otheri .ne. -1) then do i = inreg(1), inreg(2) j = numvars*(i-inreg(1)) msg = msgtype + j reqr_i = reqr_i + 1 call mpiirecv (adv_state%u3(1,1,j1m+i),plndlv,mpir8, & otheri,msg,mpicom,reqr(reqr_i)) msg = msgtype + j + 1 reqr_i = reqr_i + 1 call mpiirecv (adv_state%v3(1,1,j1m+i),plndlv,mpir8, & otheri,msg,mpicom,reqr(reqr_i)) msg = msgtype + j + 2 reqr_i = reqr_i + 1 call mpiirecv (adv_state%t3(1,1,j1m+i),plndlv,mpir8, & otheri,msg,mpicom,reqr(reqr_i)) num = plndlv do c = 1, pcnst msg = msgtype + j + 2 + c reqr_i = reqr_i + 1 call mpiirecv (adv_state%q3(1,1,j1m+i,c),num,mpir8, & otheri,msg,mpicom,reqr(reqr_i)) end do msg = msgtype + j + 3 + pcnst reqr_i = reqr_i + 1 call mpiirecv (adv_state%lnpssld(1,1,j1m+i),plndlv,mpir8, & otheri,msg,mpicom,reqr(reqr_i)) msg = msgtype + j + 4 + pcnst reqr_i = reqr_i + 1 call mpiirecv (adv_state%prhssld(1,1,j1m+i),plndlv,mpir8, & otheri,msg,mpicom,reqr(reqr_i)) msg = msgtype + j + 5 + pcnst reqr_i = reqr_i + 1 call mpiirecv (adv_state%u3sld(1,1,j1m+i),plndlv,mpir8, & otheri,msg,mpicom,reqr(reqr_i)) msg = msgtype + j + 6 + pcnst reqr_i = reqr_i + 1 call mpiirecv (adv_state%v3sld(1,1,j1m+i),plndlv,mpir8, & otheri,msg,mpicom,reqr(reqr_i)) msg = msgtype + j + 7 + pcnst reqr_i = reqr_i + 1 call mpiirecv (adv_state%etadot(1,1,j1m+i),plond*plevp,mpir8, & otheri,msg,mpicom,reqr(reqr_i)) end do end if if (reqs_i .ne. 0) then call mpiwaitall(reqs_i,reqs,stats) end if if (reqr_i .ne. 0) then call mpiwaitall(reqr_i,reqr,statr) end if return end subroutine bndexch_mpi ! !----------------------------------------------------------------------- ! subroutine intersct (regiona ,regionb ,regionc ) ! !----------------------------------------------------------------------- ! ! Purpose: ! Given two regions (a,b) output the intersection (common latitudes) ! of these two sets. The routine is used in bndexch to determine which ! latitudes need to be communicated to neighboring processors. Typically ! this routine is invoked as the intersection of the set of resident ! latitudes on processor A with the set of extended latitudes (needed for ! the SLT) of processor B. Any common latitudes will need to be ! communicated to B to complete SLT processing. ! ! Original version: J. Rosinski ! !------------------------------Arguments-------------------------------- ! integer , intent(in) :: regiona( 2 ) integer , intent(in) :: regionb( 2 ) integer , intent(out) :: regionc( 2 ) ! !----------------------------------------------------------------------- ! regionc( 1 ) = max( regiona( 1 ), regionb( 1 ) ) regionc( 2 ) = min( regiona( 2 ), regionb( 2 ) ) return end subroutine intersct #endif ! !----------------------------------------------------------------------- ! subroutine grdini(pmap ,etamid ,etaint ,gravit ,dlam , & lam ,phi ,dphi ,gw ,sinlam , & coslam ,lbasdy ,lbasdz ,lbassd ,lbasiy , & lbasiz ,lbassi ,detam ,detai ,kdpmpf , & kdpmph ,cwava ,phigs ) !----------------------------------------------------------------------- ! ! Purpose: ! Initialize model and extended grid parameters ! Initialize weights for Lagrange cubic derivative estimates ! Initialize weights for Lagrange cubic interpolant ! ! Author: J. Olson ! !----------------------------------------------------------------------- use rgrid, only: nlon !------------------------------Arguments-------------------------------- ! integer , intent(in) :: pmap ! dimension of artificial vert. grid real(r8), intent(in) :: etamid(plev) ! full-level model vertical grid real(r8), intent(in) :: etaint(plevp) ! half-level model vertical grid real(r8), intent(in) :: gravit ! gravitational constant real(r8), intent(out) :: dlam (platd) ! longitudinal grid interval (radians) real(r8), intent(out) :: lam (plond,platd) ! longitudinal coords of extended grid real(r8), intent(out) :: phi (platd) ! latitudinal coords of extended grid real(r8), intent(out) :: dphi (platd) ! latitude intervals (radians) real(r8), intent(out) :: gw (plat) ! Gaussian weights real(r8), intent(out) :: sinlam(plond,platd) ! sin(lam) model domain only real(r8), intent(out) :: coslam(plond,platd) ! cos(lam) model domain only real(r8), intent(out) :: lbasdy(4,2,platd) ! latitude derivative weights real(r8), intent(out) :: lbasdz(4,2,plev) ! vertical (full levels) deriv weights real(r8), intent(out) :: lbassd(4,2,plevp) ! vertical (half levels) deriv weights real(r8), intent(out) :: lbasiy(4,2,platd) ! Lagrange cubic interp weights (lat.) real(r8), intent(out) :: lbasiz(4,2,plev) ! Lagrange cubic interp wghts (full lev) real(r8), intent(out) :: lbassi(4,2,plevp) ! Lagrange cubic interp wghts (half lev) real(r8), intent(out) :: detam (plev) ! intervals between vertical full levs. real(r8), intent(out) :: detai (plevp) ! intervals between vertical half levs. integer , intent(out) :: kdpmpf(pmap) ! artificial full vertical grid indices integer , intent(out) :: kdpmph(pmap) ! artificial half vertical grid indices real(r8), intent(out) :: cwava(plat) ! weight applied to global integrals real(r8), intent(out) :: phigs ! Cutoff latitude for local geodesic ! !---------------------------Local variables----------------------------- ! integer j ! index integer k ! index real(r8) etamln (plev) ! log(etamid) real(r8) etailn (plevp) ! log(etaint) real(r8) detamln(plev) ! dlog(etamid) real(r8) detailn(plevp) ! dlog(etaint) ! !----------------------------------------------------------------------- ! ! Initialize extended horizontal grid coordinates. ! call grdxy(dlam ,lam ,phi ,gw ,sinlam , & coslam ) ! ! Basis functions for computing Lagrangian cubic derivatives ! on unequally spaced latitude and vertical grids. ! call basdy(phi ,lbasdy ) call basdz(plev ,etamid ,lbasdz ) call basdz(plevp ,etaint ,lbassd ) ! ! Basis functions for computing weights for Lagrangian cubic ! interpolation on unequally spaced latitude and height grids. ! call basiy(phi ,lbasiy ) call basiz(plev ,etamid ,lbasiz ) call basiz(plevp ,etaint ,lbassi ) ! ! Compute interval lengths in latitudinal grid ! do j = 1,platd-1 dphi(j) = phi(j+1) - phi(j) end do ! ! Compute interval lengths in vertical grids. ! do k = 1,plev etamln(k) = log(etamid(k)) end do do k = 1,plevp etailn(k) = log(etaint(k)) end do do k = 1,plev-1 detam (k) = etamid(k+1) - etamid(k) detamln(k) = etamln(k+1) - etamln(k) end do do k = 1,plev detai (k) = etaint(k+1) - etaint(k) detailn(k) = etailn(k+1) - etailn(k) end do ! ! Build artificial evenly spaced vertical grid for use in determining ! vertical position of departure point. ! Build one grid for full model levels and one for half levels. ! call vrtmap(plev ,pmap ,etamln ,detamln ,kdpmpf ) call vrtmap(plevp ,pmap ,etailn ,detailn ,kdpmph ) ! ! Compute tracer integral constant ! do j=1,plat cwava(j) = 1._r8/(nlon(j)*gravit) end do ! ! Initialize cutoff latitude (poleward of which local geodesic ! coordinates will be used); radians ! phigs = 1.221730_r8 ! return end subroutine grdini ! !----------------------------------------------------------------------- ! subroutine grdxy(dlam ,lam ,phi ,w ,sinlam , & coslam ) !----------------------------------------------------------------------- ! ! Purpose: ! Define the "extended" grid used in the semi-Lagrangian transport ! scheme. The longitudes are equally spaced and the latitudes are ! Gaussian. The global grid is extended to include "wraparound" points ! on all sides. ! ! Method: ! ! Author: J. Olson ! !----------------------------------------------------------------------- use gauaw_mod, only: gauaw !------------------------------Parameters------------------------------- integer, parameter :: istart = nxpt+1 ! index for first model long. integer, parameter :: jstart = nxpt+jintmx+1 ! index for first model lat. integer, parameter :: jstop = jstart-1+plat ! index for last model lat. !----------------------------------------------------------------------- !------------------------------Arguments-------------------------------- real(r8), intent(out) :: dlam(platd) ! longitudinal increment real(r8), intent(out) :: lam (plond,platd) ! long. coords. in extended grid real(r8), intent(out) :: phi (platd) ! lat. coords. in extended grid real(r8), intent(out) :: w (plat) ! Gaussian weights real(r8), intent(out) :: sinlam(plond,platd) ! sin(lam) real(r8), intent(out) :: coslam(plond,platd) ! cos(lam) ! ! dlam Length of increment in longitude grid. ! lam Longitude values in the extended grid. ! phi Latitude values in the extended grid. ! w Gauss weights for latitudes in the global grid. (These sum ! to 2.0 like the ones in CCM1.) ! sinlam Sine of longitudes in global grid (no extension points). ! coslam Cosine of longitudes in global grid (no extension points). !----------------------------------------------------------------------- !---------------------------Local variables----------------------------- integer i,j,ig ! indices integer nlond ! extended long dim real(r8) lam0 ! lamda = 0 real(r8) pi ! 3.14... real(r8) wrk(platd) ! work space !----------------------------------------------------------------------- ! lam0 = 0.0_r8 pi = 4._r8*atan(1._r8) ! ! Interval length in equally spaced longitude grid. ! do j=1,platd dlam(j) = 2._r8*pi/real(nlonex(j),r8) ! ! Longitude values on extended grid. ! nlond = nlonex(j) + 1 + 2*nxpt do i = 1,nlond lam(i,j) = real(i-istart,r8)*dlam(j) + lam0 end do end do ! ! Compute Gauss latitudes and weights. On return; phi contains the ! sine of the latitudes starting closest to the north pole and going ! toward the south; w contains the corresponding Gauss weights. ! call gauaw(phi ,w ,plat ) ! ! Reorder and compute latitude values. ! do j = jstart,jstop wrk(j) = asin( phi(jstop-j+1) ) end do phi(jstart:jstop) = wrk(jstart:jstop) ! ! North and south poles. ! phi(jstart-1) = -pi/2.0_r8 phi(jstop +1) = pi/2.0_r8 ! ! Extend Gauss latitudes below south pole so that the spacing above ! the pole is symmetric, and phi is decreasing, i.e., phi < -pi/2 ! if( jstart > 2 )then do j = 1,jstart-2 phi(j) = -pi - phi(2*jstart-2-j) end do end if ! ! Analogously for Northern Hemisphere ! if( platd > jstop+1 )then do j = jstop+2,platd phi(j) = pi - phi(2*jstop+2-j) end do end if ! ! Sine and cosine of longitude. ! do j=1,platd ig = 0 do i = istart,nlonex(j)+nxpt ig = ig + 1 sinlam(ig,j) = sin( lam(i,j) ) coslam(ig,j) = cos( lam(i,j) ) end do end do return end subroutine grdxy ! !----------------------------------------------------------------------- ! subroutine sltlinint(kdim ,ikcnst ,fb ,xl ,xr , & ys ,yn ,zt ,zb ,idp , & jdp ,kdp ,fdp ,nlon ) ! !----------------------------------------------------------------------- ! ! Purpose: ! Interpolate field to departure points using tri-linear interpolation ! ! Author: J. Olson ! !----------------------------------------------------------------------- !------------------------------Arguments-------------------------------- ! integer , intent(in) :: kdim ! vertical dimension integer , intent(in) :: ikcnst ! constituent index real(r8), intent(in) :: fb (plond,kdim*ikcnst,beglatex:endlatex) ! input field real(r8), intent(in) :: xl (plon,plev,4) ! weight for x-interpolants (left) real(r8), intent(in) :: xr (plon,plev,4) ! weight for x-interpolants (right) real(r8), intent(in) :: ys (plon,plev) ! weight for y-interpolants (south) real(r8), intent(in) :: yn (plon,plev) ! weight for y-interpolants (north) real(r8), intent(in) :: zt (plon,plev) ! top vertical interpolation weight real(r8), intent(in) :: zb (plon,plev) ! bot vertical interpolation weight integer , intent(in) :: idp(plon,plev,4) ! index of x-coordinate of dep pt integer , intent(in) :: jdp(plon,plev) ! index of y-coordinate of dep pt integer , intent(in) :: kdp(plon,plev) ! index of z-coordinate of dep pt real(r8), intent(out) :: fdp(plon,plev) ! interpolant integer , intent(in) :: nlon ! number of longitudes for this latitude ! !---------------------------Local workspace----------------------------- ! integer i, k, ii2, ii3, jj, kk ! indices ! !----------------------------------------------------------------------- ! do k=1,plev do i = 1,nlon ii2 = idp(i,k,2) ii3 = idp(i,k,3) jj = jdp(i,k) kk = kdp(i,k) fdp(i,k) = (( fb (ii2 ,kk ,jj )*xl (i,k,2) & + fb (ii2+1,kk ,jj )*xr (i,k,2) )*ys(i,k) & + ( fb (ii3 ,kk ,jj+1)*xl (i,k,3) & + fb (ii3+1,kk ,jj+1)*xr (i,k,3) )*yn(i,k) )*zt(i,k) & + (( fb (ii2 ,kk+1,jj )*xl (i,k,2) & + fb (ii2+1,kk+1,jj )*xr (i,k,2) )*ys(i,k) & + ( fb (ii3 ,kk+1,jj+1)*xl (i,k,3) & + fb (ii3+1,kk+1,jj+1)*xr (i,k,3) )*yn(i,k) )*zb(i,k) end do end do return end subroutine sltlinint ! !----------------------------------------------------------------------- ! subroutine slttraj(pmap ,jcen ,lat ,ztodt ,ra , & iter ,lam ,phi ,dphi ,etamid , & etaint ,detam ,detai ,lbasiy ,lbasiz , & lbassi ,kdpmpf ,kdpmph ,idp ,jdp , & kdp ,kkdp ,xl ,xr ,wgt1x , & wgt2x ,wgt3x ,wgt4x ,hl ,hr , & dhl ,dhr ,ys ,yn ,wgt1y , & wgt2y ,wgt3y ,wgt4y ,hs ,hn , & dhs ,dhn ,rdphi ,wgt1z ,wgt2z , & wgt3z ,wgt4z ,hb ,ht ,dhb , & dht ,rdz ,lampr ,phipr ,upr , & vpr ,lamdp ,phidp ,sigdp ,u3 , & v3 ,u3sld ,v3sld ,etadot ,u3p , & v3p ,etadotp ,dlam ,nlon ) ! !----------------------------------------------------------------------- ! ! Purpose: ! Determine trajectory departure point for each arrival point. ! Assume that the first guess for the dep point coordinates are the ! arrival point coordinates ! ! Author: J. Olson ! !----------------------------------------------------------------------- #if (!defined UNICOSMP) use srchutil #endif !------------------------------Arguments-------------------------------- ! integer , intent(in) :: pmap ! artificial vert grid dim. integer , intent(in) :: jcen ! index of lat slice(extend) integer , intent(in) :: lat ! index of lat slice (model) real(r8), intent(in) :: ztodt ! time step (seconds) real(r8), intent(in) :: ra ! 1./(radius of earth) integer , intent(in) :: iter ! iteration count real(r8), intent(in) :: lam (plond,platd) ! long. coord of model grid real(r8), intent(in) :: phi (platd) ! lat. coord of model grid real(r8), intent(in) :: dphi (platd) ! increment between lats. real(r8), intent(in) :: etamid(plev) ! vertical full levels real(r8), intent(in) :: etaint(plevp) ! vertical half levels real(r8), intent(in) :: detam (plev) ! increment between full levs real(r8), intent(in) :: detai (plevp) ! increment between half levs real(r8), intent(in) :: lbasiy(4,2,platd) ! lat interp wts(lagrng) real(r8), intent(in) :: lbasiz(4,2,plev) ! vert interp wghts (lagrng) real(r8), intent(in) :: lbassi(4,2,plevp) ! vert interp wghts (lagrng) integer , intent(in) :: kdpmpf(pmap) ! artificial vert grid indices integer , intent(in) :: kdpmph(pmap) ! artificial vert grid indices integer , intent(out) :: idp (plon,plev,4) ! zonal dep point index integer , intent(out) :: jdp (plon,plev) ! meridional dep point index integer , intent(out) :: kdp (plon,plev) ! vertical dep point index integer , intent(in) :: kkdp (plon,plev) ! index of z-coordinate of dep pt (alt) real(r8), intent(out) :: xl (plon,plev,4) ! weight for x-interpolants (left) real(r8), intent(out) :: xr (plon,plev,4) ! weight for x-interpolants (right) real(r8), intent(in) :: wgt1x (plon,plev,4) ! weight for x-interpolants (Lag Cubic) real(r8), intent(in) :: wgt2x (plon,plev,4) ! weight for x-interpolants (Lag Cubic) real(r8), intent(in) :: wgt3x (plon,plev,4) ! weight for x-interpolants (Lag Cubic) real(r8), intent(in) :: wgt4x (plon,plev,4) ! weight for x-interpolants (Lag Cubic) real(r8), intent(in) :: hl (plon,plev,4) ! weight for x-interpolants (Hermite) real(r8), intent(in) :: hr (plon,plev,4) ! weight for x-interpolants (Hermite) real(r8), intent(in) :: dhl (plon,plev,4) ! weight for x-interpolants (Hermite) real(r8), intent(in) :: dhr (plon,plev,4) ! weight for x-interpolants (Hermite) real(r8), intent(out) :: ys (plon,plev) ! weight for y-interpolants (south) real(r8), intent(out) :: yn (plon,plev) ! weight for y-interpolants (north) real(r8), intent(in) :: wgt1y (plon,plev) ! weight for y-interpolants (Lag Cubic) real(r8), intent(in) :: wgt2y (plon,plev) ! weight for y-interpolants (Lag Cubic) real(r8), intent(in) :: wgt3y (plon,plev) ! weight for y-interpolants (Lag Cubic) real(r8), intent(in) :: wgt4y (plon,plev) ! weight for y-interpolants (Lag Cubic) real(r8), intent(in) :: hs (plon,plev) ! weight for y-interpolants (Hermite) real(r8), intent(in) :: hn (plon,plev) ! weight for y-interpolants (Hermite) real(r8), intent(in) :: dhs (plon,plev) ! weight for y-interpolants (Hermite) real(r8), intent(in) :: dhn (plon,plev) ! weight for y-interpolants (Hermite) real(r8), intent(in) :: rdphi (plon,plev) ! reciprocal of y-interval real(r8), intent(in) :: wgt1z (plon,plev) ! weight for z-interpolants (Lag Cubic) real(r8), intent(in) :: wgt2z (plon,plev) ! weight for z-interpolants (Lag Cubic) real(r8), intent(in) :: wgt3z (plon,plev) ! weight for z-interpolants (Lag Cubic) real(r8), intent(in) :: wgt4z (plon,plev) ! weight for z-interpolants (Lag Cubic) real(r8), intent(in) :: hb (plon,plev) ! weight for z-interpolants (Hermite) real(r8), intent(in) :: ht (plon,plev) ! weight for z-interpolants (Hermite) real(r8), intent(in) :: dhb (plon,plev) ! weight for z-interpolants (Hermite) real(r8), intent(in) :: dht (plon,plev) ! weight for z-interpolants (Hermite) real(r8), intent(in) :: rdz (plon,plev) ! reciprocal of z-interval real(r8), intent(out) :: lampr (plon,plev) ! trajectory increment (x-direction) real(r8), intent(out) :: phipr (plon,plev) ! trajectory increment (y-direction) real(r8), intent(out) :: upr (plon,plev) ! interpolated u field (local geodesic) real(r8), intent(out) :: vpr (plon,plev) ! interpolated v field (local geodesic) real(r8), intent(out) :: lamdp (plon,plev) ! zonal departure pt. coord. real(r8), intent(inout):: phidp (plon,plev) ! meridional departure pt. coord. real(r8), intent(inout):: sigdp (plon,plev) ! vertical departure pt. coord. real(r8), intent(in) :: u3 (plond, plev, beglatex:endlatex) ! u wind component real(r8), intent(in) :: v3 (plond, plev, beglatex:endlatex) ! v wind component real(r8), intent(in) :: u3sld (plond,plev ,beglatex:endlatex) ! u3 inpt for SLD int real(r8), intent(in) :: v3sld (plond,plev ,beglatex:endlatex) ! v3 inpt for SLD int real(r8), intent(in) :: etadot(plond,plevp,beglatex:endlatex) ! Vertical motion real(r8), intent(in) :: u3p (plon,plev) ! u3 wind at current time real(r8), intent(in) :: v3p (plon,plev) ! v3 wind at current time real(r8), intent(in) :: etadotp(plon,plevp) ! Vertical motion at current time real(r8), intent(in) :: dlam (platd) ! d-long for each lat integer , intent(in) :: nlon ! # of longitudes #if (defined UNICOSMP) integer,external :: ismin #endif ! !---------------------------Local workspace----------------------------- ! real(r8) fac ! 1-eps real(r8) tmp ! temp variable real(r8) ump (plon,plev) ! interpolated u field real(r8) vmp (plon,plev) ! interpolated v field real(r8) wmp (plon,plev) ! interpolated w field real(r8) wmpa (plon,plev) ! interpolated w field (arrival level) real(r8) sigpr (plon,plev) ! trajectory increment (vertical) real(r8) etamin ! minimum eta real(r8) etamax ! maximum eta real(r8) delsig (plev) ! dist between dep point and model levs real(r8) zt (plon,plev) ! top vertical interpolation weight real(r8) zb (plon,plev) ! bot vertical interpolation weight real(r8) zttmp ! top vertical interpolation weight real(r8) zbtmp ! bot vertical interpolation weight real(r8) rdx (platd) ! reciprocal of del-x logical locgeo ! local geodesic flag logical lvert ! flag to compute vertical trajectory logical larrival(plon) ! flag to indicate whether or not to ! ! place the alternative dep point at the ! ! arrival point. integer kount ! index counter integer kkmin ! index of minimum "delsig" integer i ! | integer ii ! | integer j ! | integer jj2 ! | - indices integer jj3 ! | integer k ! | integer kk ! | integer n ! | ! !----------------------------------------------------------------------- ! fac = 1._r8 - 10._r8*epsilon (fac) locgeo = .false. if(abs(phi(jcen)) .ge. phigs) locgeo = .true. ! ! Interpolate arrival and departure vertical velocities ! do k = 1,plev do i = 1,nlon zttmp = ( etaint(k+1) - sigdp(i,k) )/detai(k) zbtmp = 1._r8 - zttmp wmpa(i,k) = etadotp(i,k )*zttmp & + etadotp(i,k+1 )*zbtmp wmp (i,k) = etadot(i1+i-1,k ,jcen)*zttmp & + etadot(i1+i-1,k+1,jcen)*zbtmp end do end do ! ! Set up computation of trajectory ! do k = 1,plev do i = 1,nlon ii = i + i1 - 1 ! ! Place u/v on unit sphere ! ump(i,k) = 0.5_r8*(u3sld(ii,k,jcen) + u3p(i,k))*ra vmp(i,k) = 0.5_r8*(v3sld(ii,k,jcen) + v3p(i,k))*ra wmp(i,k) = 0.5_r8*(wmp(i,k) + wmpa(i,k) ) ! ! Estimate departure point of parcel trajectory. ! lampr(i,k) = -ztodt*ump(i,k) phipr(i,k) = -ztodt*vmp(i,k) sigpr(i,k) = -ztodt*wmp(i,k) if (.not. locgeo) lampr(i,k) = lampr(i,k)/cos( phidp(i,k) ) ! ! Initialize winds for use in g.c. calculations. ! if (locgeo) then upr(i,k) = ump(i,k) vpr(i,k) = vmp(i,k) end if end do end do ! ! Estimate initial trajectory departure points ! lvert = .true. call trajdp(lat ,jcen ,ztodt ,locgeo ,lvert , & lam(1,jcen),phi ,etamid ,upr ,vpr , & lampr ,phipr ,sigpr ,lamdp ,phidp , & sigdp ,nlon ) ! ! Loop over departure point iterates. ! do n = 1,iter ! ! Determine departure point indicies and interpolate u,v,w ! call bandij (dlam ,phi ,lamdp ,phidp ,idp , & jdp ,nlon ) call kdpfnd (plev ,pmap ,etamid ,sigdp ,kdpmpf , & kdp ,nlon ) ! ! Compute weights for x,y,z dimensions ! do j = 1,platd rdx(j) = 1._r8/(lam(nxpt+2,j) - lam(nxpt+1,j)) end do do k = 1,plev do i = 1,nlon jj2 = jdp(i,k) jj3 = jdp(i,k) + 1 xl(i,k,2) = (lam(idp(i,k,2)+1,jj2) - lamdp(i,k))*rdx(jj2) xl(i,k,3) = (lam(idp(i,k,3)+1,jj3) - lamdp(i,k))*rdx(jj3) xr(i,k,2) = 1._r8 - xl(i,k,2) xr(i,k,3) = 1._r8 - xl(i,k,3) ys(i,k) = ( phi(jdp(i,k)+1) - phidp(i,k) )/dphi(jdp(i,k)) yn(i,k) = 1._r8 - ys(i,k) zt(i,k) = (etamid(kdp(i,k)+1)-sigdp(i,k))/detam(kdp(i,k)) zb(i,k) = 1._r8 - zt(i,k) end do end do ! call sltlinint(plev ,1 ,u3sld(1,1,beglatex),xl ,xr , & ys ,yn ,zt ,zb ,idp , & jdp ,kdp ,ump ,nlon ) call sltlinint(plev ,1 ,v3sld(1,1,beglatex),xl ,xr , & ys ,yn ,zt ,zb ,idp , & jdp ,kdp ,vmp ,nlon ) call kdpfnd (plevp ,pmap ,etaint ,sigdp ,kdpmph , & kdp ,nlon ) do k = 1,plev do i = 1,nlon zt(i,k) = (etaint(kdp(i,k)+1)-sigdp(i,k))/detai(kdp(i,k)) zb(i,k) = 1._r8 - zt(i,k) end do end do call sltlinint(plevp ,1 ,etadot(1,1,beglatex),xl ,xr , & ys ,yn ,zt ,zb ,idp , & jdp ,kdp ,wmp ,nlon ) ! ! Add arrival point velocities (n) to interpolated velocities ! do k = 1,plev do i = 1,nlon ii = i + i1 - 1 ump(i,k) = 0.5_r8*(ump(i,k) + u3p(i,k)) vmp(i,k) = 0.5_r8*(vmp(i,k) + v3p(i,k)) wmp(i,k) = 0.5_r8*(wmp(i,k) + wmpa(i,k) ) end do end do ! ! Compute new trajectory departure points ! lvert = .true. call depinc (jcen ,ztodt ,ra ,locgeo ,lvert , & lam(1,jcen),phi ,ump ,vmp ,wmp , & upr ,vpr ,lamdp ,phidp ,lampr , & phipr ,sigpr ,nlon ) call trajdp (lat ,jcen ,ztodt ,locgeo ,lvert , & lam(1,jcen),phi ,etamid ,upr ,vpr , & lampr ,phipr ,sigpr ,lamdp ,phidp , & sigdp ,nlon ) end do ! ! Compile departure point binning statistics ! do k = 1,plev if(k .eq. 1) then etamin = etamid(k ) etamax = etamid(k+1) elseif(k .eq. plev) then etamin = etamid(k-1) etamax = etamid(k ) else etamin = etamid(k-1) etamax = etamid(k+1) end if do i = 1,nlon larrival(i) = sigdp(i,k) .ge. etamin .and. sigdp(i,k) .le. etamax end do ! ! If: Departure point is within one grid interval from the arrival ! point, bin it in the ARRIVAL point bin ! tmp = 0._r8 do i = 1,nlon if (larrival(i)) then tmp = tmp + 1._r8 end if end do levkntl(k,k,lat) = levkntl(k,k,lat) + tmp ! ! Else: find departure point bin ! kount = 0 do i=1,nlon if (larrival(i)) kount = kount + 1 end do if (kount .ne. nlon) then do i = 1,nlon if (.not. larrival(i)) then delsig(1) = 1.e+35_r8 do kk = 2,plev delsig(kk) = abs( sigdp(i,k) - etaint(kk) ) end do kkmin = ismin( plev,delsig,1 ) if(kkmin .gt. k) kkmin = kkmin-1 levkntl(kkmin,k,lat) = levkntl(kkmin,k,lat) + 1._r8 end if end do end if end do ! return end subroutine slttraj ! !----------------------------------------------------------------------- ! subroutine sltwgts(limdrh ,limdrv ,lhrzwgt ,lvrtwgt ,kdim , & idp ,jdp ,kdp ,lam ,phi , & z ,dphi ,dz ,lamdp ,phidp , & sigdp ,lbasiy ,wiz ,kkdp ,xl , & xr ,wgt1x ,wgt2x ,wgt3x ,wgt4x , & hl ,hr ,dhl ,dhr ,ys , & yn ,wgt1y ,wgt2y ,wgt3y ,wgt4y , & hs ,hn ,dhs ,dhn ,rdphi , & wgt1z ,wgt2z ,wgt3z ,wgt4z ,hb , & ht ,dhb ,dht ,rdz ,zt , & zb ,nlon ) !----------------------------------------------------------------------- ! ! Purpose: ! Compute weights for SLT interpolation ! ! Author: J. Olson ! !----------------------------------------------------------------------- #if (!defined UNICOSMP) use srchutil, only: whenieq #endif !------------------------------Arguments-------------------------------- ! logical , intent(in) :: limdrh ! horizontal derivative limiter flag logical , intent(in) :: limdrv ! vertical derivative limiter flag logical , intent(in) :: lhrzwgt ! flag to compute horizontal weights logical , intent(in) :: lvrtwgt ! flag to compute vertical weights integer , intent(in) :: kdim ! vertical coordinate integer , intent(in) :: idp (plon,plev,4) ! index of x-coordinate of dep pt integer , intent(in) :: jdp (plon,plev) ! index of y-coordinate of dep pt integer , intent(in) :: kdp (plon,plev) ! index of z-coordinate of dep pt real(r8), intent(in) :: lam (plond,platd) ! longitude coordinates of model grid real(r8), intent(in) :: phi (platd) ! latitude coordinates of model grid real(r8), intent(in) :: z (kdim) ! vertical coordinates of model grid real(r8), intent(in) :: dphi (platd) ! latitudinal grid increments real(r8), intent(in) :: dz (kdim) ! vertical grid increments real(r8), intent(in) :: lamdp (plon,plev) ! x-coordinates of dep pts. real(r8), intent(in) :: phidp (plon,plev) ! y-coordinates of dep pt real(r8), intent(in) :: sigdp (plon,plev) ! z-coordinates of dep pt real(r8), intent(in) :: lbasiy(4,2,platd) ! y-interpolation weights for Lag. cubic real(r8), intent(in) :: wiz (4,2,kdim) ! z-interpolation weights for Lag. cubic integer , intent(out) :: kkdp (plon,plev) ! index of z-coordinate of dep pt (alt) real(r8), intent(out) :: xl (plon,plev,4) ! weight for x-interpolants (left) real(r8), intent(out) :: xr (plon,plev,4) ! weight for x-interpolants (right) real(r8), intent(out) :: wgt1x (plon,plev,4) ! weight for x-interpolants (Lag Cubic) real(r8), intent(out) :: wgt2x (plon,plev,4) ! weight for x-interpolants (Lag Cubic) real(r8), intent(out) :: wgt3x (plon,plev,4) ! weight for x-interpolants (Lag Cubic) real(r8), intent(out) :: wgt4x (plon,plev,4) ! weight for x-interpolants (Lag Cubic) real(r8), intent(out) :: hl (plon,plev,4) ! weight for x-interpolants (Hermite) real(r8), intent(out) :: hr (plon,plev,4) ! weight for x-interpolants (Hermite) real(r8), intent(out) :: dhl (plon,plev,4) ! weight for x-interpolants (Hermite) real(r8), intent(out) :: dhr (plon,plev,4) ! weight for x-interpolants (Hermite) real(r8), intent(out) :: ys (plon,plev) ! weight for y-interpolants (south) real(r8), intent(out) :: yn (plon,plev) ! weight for y-interpolants (north) real(r8), intent(out) :: wgt1y (plon,plev) ! weight for y-interpolants (Lag Cubic) real(r8), intent(out) :: wgt2y (plon,plev) ! weight for y-interpolants (Lag Cubic) real(r8), intent(out) :: wgt3y (plon,plev) ! weight for y-interpolants (Lag Cubic) real(r8), intent(out) :: wgt4y (plon,plev) ! weight for y-interpolants (Lag Cubic) real(r8), intent(out) :: hs (plon,plev) ! weight for y-interpolants (Hermite) real(r8), intent(out) :: hn (plon,plev) ! weight for y-interpolants (Hermite) real(r8), intent(out) :: dhs (plon,plev) ! weight for y-interpolants (Hermite) real(r8), intent(out) :: dhn (plon,plev) ! weight for y-interpolants (Hermite) real(r8), intent(out) :: rdphi (plon,plev) ! reciprocal of y-interval real(r8), intent(out) :: wgt1z (plon,plev) ! weight for z-interpolants (Lag Cubic) real(r8), intent(out) :: wgt2z (plon,plev) ! weight for z-interpolants (Lag Cubic) real(r8), intent(out) :: wgt3z (plon,plev) ! weight for z-interpolants (Lag Cubic) real(r8), intent(out) :: wgt4z (plon,plev) ! weight for z-interpolants (Lag Cubic) real(r8), intent(out) :: hb (plon,plev) ! weight for z-interpolants (Hermite) real(r8), intent(out) :: ht (plon,plev) ! weight for z-interpolants (Hermite) real(r8), intent(out) :: dhb (plon,plev) ! weight for z-interpolants (Hermite) real(r8), intent(out) :: dht (plon,plev) ! weight for z-interpolants (Hermite) real(r8), intent(out) :: rdz (plon,plev) ! reciprocal of z-interval real(r8), intent(out) :: zt (plon,plev) ! linear interpolation weight real(r8), intent(out) :: zb (plon,plev) ! linear interpolation weight integer , intent(in) :: nlon ! number of longitudes for this latitude ! !---------------------------Local workspace----------------------------- ! integer i ! | integer ii,jj(plon,plev,4)! | integer k,n,j ! | integer icount ! | integer jdpval ! | integer kdpval ! | -- indices integer jmin ! | integer jmax ! | integer kdimm2 ! | integer nval ! | integer indx (plond) ! | ! real(r8) dx (platd) ! | real(r8) rdx (platd) ! | real(r8) dyj ! | real(r8) tmp1 ! | real(r8) tmp2 ! | real(r8) tmp3 ! | real(r8) tmp4 ! | -- tmp variables real(r8) dzk ! | real(r8) denom1 ! | real(r8) denom2 ! | real(r8) denom3 ! | real(r8) denom4 ! | real(r8) coef12 ! | real(r8) coef34 ! | ! !----------------------------------------------------------------------- ! denom1 = -1._r8/6._r8 denom2 = 0.5_r8 denom3 = -0.5_r8 denom4 = 1._r8/6._r8 ! ! HORIZONTAL weights ! if (lhrzwgt) then ! ! Determine N/S extent of all departure points in this latitude slice ! jmin = 1000000 jmax = -1000000 do k=1,plev do i=1,nlon if(jdp(i,k) .lt. jmin) jmin = jdp(i,k) if(jdp(i,k) .gt. jmax) jmax = jdp(i,k) end do end do ! ! Compute weights for x-direction ! do j = 1,platd dx (j) = lam(nxpt+2,j) - lam(nxpt+1,j) rdx(j) = 1._r8/dx(j) end do ! do n=1,4 do k=1,plev do i=1,nlon jj(i,k,n) = jdp(i,k) - 2 + n xl(i,k,n) = (lam(idp(i,k,n)+1,jj(i,k,n)) - lamdp(i,k))*rdx(jj(i,k,n)) xr(i,k,n) = 1._r8 - xl(i,k,n) end do end do end do do n=2,3 do k=1,plev do i=1,nlon hl (i,k,n) = ( 3.0_r8 - 2.0_r8*xl(i,k,n) )*xl(i,k,n)**2 hr (i,k,n) = ( 3.0_r8 - 2.0_r8*xr(i,k,n) )*xr(i,k,n)**2 dhl(i,k,n) = -dx(jj(i,k,n))*( xl(i,k,n) - 1._r8 )*xl(i,k,n)**2 dhr(i,k,n) = dx(jj(i,k,n))*( xr(i,k,n) - 1._r8 )*xr(i,k,n)**2 end do end do do k=1,plev do i=1,nlon tmp1 = xr(i,k,n) + 1._r8 tmp4 = xr(i,k,n) - 2._r8 coef12 = -xl(i,k,n)*tmp4 coef34 = xr(i,k,n)*tmp1 wgt1x(i,k,n) = denom1*coef12*xr(i,k,n) wgt2x(i,k,n) = denom2*coef12*tmp1 wgt3x(i,k,n) = denom3*coef34*tmp4 wgt4x(i,k,n) = -denom4*coef34*xl(i,k,n) end do end do end do ! ! Compute weights for y-direction ! icount = 0 do jdpval=jmin,jmax do k=1,plev call whenieq(nlon ,jdp(1,k),1 ,jdpval ,indx ,nval ) icount = icount + nval dyj = dphi(jdpval) do ii = 1,nval i = indx(ii) ys(i,k) = ( phi(jdpval+1) - phidp(i,k) )/dyj yn(i,k) = 1._r8 - ys(i,k) end do do ii = 1,nval i = indx(ii) rdphi(i,k) = 1._r8/dyj hs (i,k) = ( 3.0_r8 - 2.0_r8*ys(i,k) )*ys(i,k)**2 hn (i,k) = ( 3.0_r8 - 2.0_r8*yn(i,k) )*yn(i,k)**2 dhs (i,k) = -dyj*( ys(i,k) - 1._r8 )*ys(i,k)**2 dhn (i,k) = dyj*( yn(i,k) - 1._r8 )*yn(i,k)**2 end do do ii = 1,nval i = indx(ii) tmp1 = phidp(i,k) - lbasiy(1,1,jdpval) tmp2 = phidp(i,k) - lbasiy(2,1,jdpval) tmp3 = phidp(i,k) - lbasiy(3,1,jdpval) tmp4 = phidp(i,k) - lbasiy(4,1,jdpval) coef12 = tmp3*tmp4 coef34 = tmp1*tmp2 wgt1y(i,k) = coef12*tmp2*lbasiy(1,2,jdpval) wgt2y(i,k) = coef12*tmp1*lbasiy(2,2,jdpval) wgt3y(i,k) = coef34*tmp4*lbasiy(3,2,jdpval) wgt4y(i,k) = coef34*tmp3*lbasiy(4,2,jdpval) end do end do end do if (icount.ne.nlon*plev) then call endrun ('SLTWGTS: Did not complete computations for all departure points') end if end if ! ! VERTICAL weights ! if (lvrtwgt) then ! ! Limit kdp to between "2" and "kdim-2" when computing weights and ! derivatives. ! kdimm2 = kdim - 2 do k=1,plev do i=1,nlon kkdp(i,k) = min0( kdimm2,max0( 2,kdp(i,k) ) ) dzk = dz(kdp(i,k)) rdz(i,k) = 1._r8/dzk zt(i,k) = ( z (kdp(i,k)+1) - sigdp(i,k) )/dzk zb(i,k) = 1._r8 - zt(i,k) ht (i,k) = ( 3.0_r8 - 2.0_r8*zt(i,k) )*zt(i,k)**2 hb (i,k) = ( 3.0_r8 - 2.0_r8*zb(i,k) )*zb(i,k)**2 dht(i,k) = -dzk*( zt(i,k) - 1._r8 )*zt(i,k)**2 dhb(i,k) = dzk*( zb(i,k) - 1._r8 )*zb(i,k)**2 end do end do ! if(.not. limdrv) then icount = 0 do kdpval=2,kdimm2 do k=1,plev call whenieq(nlon ,kkdp(1,k),1 ,kdpval ,indx ,nval ) icount = icount + nval do ii = 1,nval i = indx(ii) tmp1 = sigdp(i,k) - wiz(1,1,kdpval) tmp2 = sigdp(i,k) - wiz(2,1,kdpval) tmp3 = sigdp(i,k) - wiz(3,1,kdpval) tmp4 = sigdp(i,k) - wiz(4,1,kdpval) coef12 = tmp3*tmp4 coef34 = tmp1*tmp2 wgt1z(i,k) = coef12*tmp2*wiz(1,2,kdpval) wgt2z(i,k) = coef12*tmp1*wiz(2,2,kdpval) wgt3z(i,k) = coef34*tmp4*wiz(3,2,kdpval) wgt4z(i,k) = coef34*tmp3*wiz(4,2,kdpval) end do end do end do if (icount.ne.nlon*plev) then call endrun ('SLTWGTS: Did not complete computations for all departure points') end if end if end if ! return end subroutine sltwgts ! !----------------------------------------------------------------------- ! subroutine trajdp(lat ,jcen ,dt ,locgeo ,lvert , & lam ,phi ,etamid ,upr ,vpr , & lampr ,phipr ,sigpr ,lamdp ,phidp , & sigdp ,nlon ) ! !----------------------------------------------------------------------- ! ! Purpose: ! Determine trajectory departure point for each arrival point based upon ! departure point increments ! ! Author: J. Olson ! !----------------------------------------------------------------------- !------------------------------Arguments-------------------------------- ! integer , intent(in) :: lat ! index of lat slice (model) integer , intent(in) :: jcen ! index of lat slice(extend) real(r8), intent(in) :: dt ! time step (seconds) logical , intent(in) :: locgeo ! local geodesic flag logical , intent(in) :: lvert ! flag to compute vertical trajectory real(r8), intent(in) :: lam (plond) ! long. coord of model grid real(r8), intent(in) :: phi (platd) ! lat. coord of model grid real(r8), intent(in) :: etamid(plev) ! vertical full levels real(r8), intent(in) :: upr (plon,plev) ! interpolated u field (local geodesic) real(r8), intent(in) :: vpr (plon,plev) ! interpolated v field (local geodesic) real(r8), intent(in) :: lampr (plon,plev) ! trajectory increment (x-direction) real(r8), intent(in) :: phipr (plon,plev) ! trajectory increment (y-direction) real(r8), intent(in) :: sigpr (plon,plev) ! trajectory increment (vertical) real(r8), intent(out) :: lamdp (plon,plev) ! zonal departure pt. coord. real(r8), intent(out) :: phidp (plon,plev) ! meridional departure pt. coord. real(r8), intent(out) :: sigdp (plon,plev) ! vertical departure pt. coord. integer , intent(in) :: nlon ! number of longitudes for this latitude ! !---------------------------Local workspace----------------------------- ! real(r8) fac ! 1 - eps real(r8) botlim ! bottom limit of the trajectory real(r8) phi0 ! Current latitude (radians) real(r8) cphi0 ! cos latitude real(r8) sphi0 ! sin latitude real(r8) dist ! computed distance**2 real(r8) distmx ! max distance**2 real(r8) pi2 ! pi/2 real(r8) phipi2 ! pi/2 real(r8) sgnphi0 ! holds sign of phi0 real(r8) pi ! pi real(r8) twopi ! 2*pi ! real(r8) clamgc ! | real(r8) cphigc ! | real(r8) dlamsc ! | real(r8) slamgc ! | -- temporary variables real(r8) slam2 ! | real(r8) sphigc ! | ! integer i ! | integer ii ! | integer k ! | -- indices integer ibad ! | integer kbad ! | ! logical lstop ! flag to stop run if departure point is over the pole ! !----------------------------------------------------------------------- ! fac = 1._r8 - 10._r8*epsilon(fac) pi = 4._r8*atan(1._r8) twopi = 2._r8*pi pi2 = pi/2._r8 phi0 = phi(jcen) cphi0 = cos( phi0 ) sphi0 = sin( phi0 ) sgnphi0 = sign( 1._r8, phi0 ) distmx = (sign(pi2,phi0) - phi0)/(1.1_r8*dt) distmx = distmx*distmx ! ! Compute coordinates of departure point. ! lstop = .false. do k = 1,plev ! ! if near pole, convert from g.c. to spherical ! if (locgeo) then do i = 1,nlon ii = i + i1 - 1 sphigc = sin( phipr(i,k) ) cphigc = cos( phipr(i,k) ) slamgc = sin( lampr(i,k) ) clamgc = cos( lampr(i,k) ) phidp(i,k) = asin((sphigc*cphi0 + cphigc*sphi0*clamgc)*fac) if ( abs(phidp(i,k)) .ge. phi(j1+plat)*fac ) & phidp(i,k) = sign( phi(j1+plat),phidp(i,k) )*fac dlamsc = asin((slamgc*cphigc/cos(phidp(i,k)))*fac) ! ! If traj is over pole, check for proper branch of arcsin ! dist = upr(i,k)**2 + vpr(i,k)**2 if( dist .gt. distmx) then slam2 = slamgc**2 phipi2 = asin((sqrt((slam2 - 1._r8)/(slam2 - 1._r8/cphi0**2)))*fac) if(sgnphi0*phipr(i,k) .gt. phipi2) dlamsc = sign(pi,lampr(i,k)) - dlamsc end if lamdp(i,k) = lam(ii) + dlamsc end do else do i = 1,nlon ii = i + i1 - 1 lamdp(i,k) = lam(ii) + lampr(i,k) phidp(i,k) = phi0 + phipr(i,k) ! ! If traj is over pole, STOP ! if (phidp(i,k) >= phi(j1+plat)) then lstop = .true. ibad = i kbad = k end if if (phidp(i,k) < phi(j1-1 )) then lstop = .true. ibad = i kbad = k end if end do end if #if ( defined SPMD ) ! ! If traj goes off-processor (out-of-bounds), STOP ! do i = 1,nlon if (phidp(i,k) >= phi(endlatex-nxpt) ) then lstop = .true. ibad = i kbad = k end if if (phidp(i,k) < phi(beglatex+nxpt) ) then lstop = .true. ibad = i kbad = k end if end do #endif ! ! Apply appropriate limiters ! do i = 1,nlon if(lamdp(i,k) .ge. twopi) lamdp(i,k) = lamdp(i,k) - twopi + 10._r8*epsilon(lamdp) if(lamdp(i,k) .lt. 0._r8) lamdp(i,k) = lamdp(i,k) + twopi - 10._r8*epsilon(lamdp) end do ! ! Compute vertical departure points and apply appropriate limiters ! if (lvert) then botlim = etamid(plev)*fac do i = 1,nlon sigdp(i,k) = etamid(k) + sigpr(i,k) if(sigdp(i,k) .lt. etamid( 1)) sigdp(i,k) = etamid(1) if(sigdp(i,k) .ge. etamid(plev)) sigdp(i,k) = botlim end do end if end do ! ! Test that the latitudinal extent of trajectory is NOT over the poles ! if (lstop) then write(iulog,*)'ERROR IN TRAJDP: ****** MODEL IS BLOWING UP *********' write(iulog,9000) ibad,kbad,lat call endrun end if ! return ! ! Formats ! 9000 format(//'Departure point out of bounds. Parcel associated with longitude ' & ,i5,', level ',i5, ' and latitude ',i5/' is outside the model domain. ') end subroutine trajdp ! !----------------------------------------------------------------------- ! end module scanslt