! following contains library for version 3 of rtlnmc (regional tangent linear normal mode constraint. ! ! module helmholtz_fmg_mod ! multigrid solver on fixed domain with 5 point stencil and dirichlet bc. ! TO DO: ! coded tested (NOTE: pay special attention to boundary treatment--likely several errors to fix) ! 20090215 20090215 grid heirarchy generator ! 20090216 20090216 transfer operators ! 20090216 20090216 bilinear interpolation operator ! 20090216 20090216 bicubic interpolation operaror ! 20090216 20090217 forward operator ( f = delsqr(v) - r*v ) ! 20090216 20090220 relaxation operator ! 20090216 20090225 fmg algorithm ! 20090227 20090309 fmg_ad use kinds, only: r_kind,i_kind implicit none private public :: mgvar public :: generate_grids public :: check_xbdy_err public :: check_ybdy_err public :: f2c_full_weight public :: c2f_bilinear public :: c2f_bicubic public :: init_mgvars public :: fmg public :: fmg_ad type mgvar sequence integer(i_kind) nx ! x dimension is 0:nx east-west integer(i_kind) ny ! y dimension is 0:ny north-south real(r_kind) phi0 ! scale geopotential real(r_kind) a_east,b_east,a_east_inv ! for dirichlet boundary conditions, interpolation real(r_kind) a_north,b_north,a_north_inv ! boundary weights on east and north bdy real(r_kind) b_over_a_east,b_over_a_east_neg ! boundary values (note, with fmg these all = 0 except real(r_kind) b_over_a_north,b_over_a_north_neg ! on finest grid, where solution is desired. real(r_kind) half_b_over_a_east_neg ! but this code is specific to zero bc. real(r_kind) half_b_over_a_north_neg ! real(r_kind) half_b_over_a_east_north_neg ! boundary formula (weights are for lin interp): ! ! a_east*v(j,nx)+b_east*v(j,nx-1) = 0 ! for j=1,ny-1 ! ! a_north*v(ny,i)+b_north*v(ny-1,i) = 0 ! for i=1,nx-1 real(r_kind),pointer:: f(:,:) => NULL() ! residual (forcing on finest grid) real(r_kind),pointer:: v(:,:) => NULL() ! correction (solution on finest grid) real(r_kind),pointer:: w(:,:) => NULL() ! scratch space real(r_kind),pointer:: w2(:,:) => NULL() ! scratch space real(r_kind),pointer:: y(:,:) => NULL() ! scratch space real(r_kind),pointer:: z(:,:) => NULL() ! scratch space real(r_kind),pointer:: yy(:,:) => NULL() ! scratch space real(r_kind),pointer:: zz(:,:) => NULL() ! scratch space real(r_kind),pointer:: x00(:,:) => NULL() ! increment in x direction (meters) real(r_kind),pointer:: y00(:,:) => NULL() ! increment in y direction (meters) real(r_kind),pointer:: h00(:,:) => NULL() ! all forward real(r_kind),pointer:: r00(:,:) => NULL() ! operator constants real(r_kind),pointer:: h0m(:,:) => NULL() ! all forward real(r_kind),pointer:: h0p(:,:) => NULL() ! operator constants real(r_kind),pointer:: hm0(:,:) => NULL() ! all forward real(r_kind),pointer:: hp0(:,:) => NULL() ! operator constants real(r_kind),pointer:: dbarx(:,:) => NULL() ! x direction (for poisson problem) real(r_kind),pointer:: ubarx(:,:) => NULL() ! adi real(r_kind),pointer:: obarx(:,:) => NULL() ! coefs real(r_kind),pointer:: dbary(:,:) => NULL() ! y direction real(r_kind),pointer:: ubary(:,:) => NULL() ! adi real(r_kind),pointer:: obary(:,:) => NULL() ! coefs real(r_kind),pointer:: dbarxh(:,:) => NULL() ! x direction (for helmholtz problem) real(r_kind),pointer:: ubarxh(:,:) => NULL() ! adi real(r_kind),pointer:: obarxh(:,:) => NULL() ! coefs real(r_kind),pointer:: dbaryh(:,:) => NULL() ! y direction real(r_kind),pointer:: ubaryh(:,:) => NULL() ! adi real(r_kind),pointer:: obaryh(:,:) => NULL() ! coefs logical,pointer:: bicubic_mask(:,:) => NULL() ! flags points that are covered by bicubic interpolation end type mgvar integer(i_kind) mgridx,mgridy ! number of possible grids in each direction integer(i_kind) mgrid ! number of grids used in multigrid solution ! ( = min(mgridx,mgridy) integer(i_kind) nxall(100),nyall(100) ! dimensions of each grid integer(i_kind) nxa,nya ! dimensions of analysis grid (nya,nxa) ! nxall(1) type(mgvar),save:: mg(100) contains subroutine generate_grids(nxa_in,nya_in) use constants, only: half,one integer(i_kind),intent(in):: nxa_in,nya_in ! grids have south and west row in common, but allow north and east row to go one grid point beyond ! north and east boundary of finest grid. by doing this, it is not necessary to have any particular ! number of points, ie 2**p+1, so that all boundaries are part of every grid. except for the north ! and east points outside the analysis domain, all interior points of all grids coincide with points ! of the target analysis grid. ! w e ! ..................................................................... ! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ! . . . . . . . . . . . . . . . . . . ! . . . . . . . . . x . ! . . . . . x . ! . . . x . integer(i_kind) k real(r_kind) a_end(100),b_end(100) nxa=nxa_in nya=nya_in call generate_1d_grid(nxa-1,nxall,mgridx,a_end,b_end) do k=1,mgridx mg(k)%nx=nxall(k) mg(k)%a_east=a_end(k) mg(k)%a_east_inv=one/a_end(k) mg(k)%b_east=b_end(k) mg(k)%b_over_a_east=b_end(k)/a_end(k) mg(k)%b_over_a_east_neg=-b_end(k)/a_end(k) mg(k)%half_b_over_a_east_neg=-half*b_end(k)/a_end(k) end do call generate_1d_grid(nya-1,nyall,mgridy,a_end,b_end) do k=1,mgridy mg(k)%ny=nyall(k) mg(k)%a_north=a_end(k) mg(k)%a_north_inv=one/a_end(k) mg(k)%b_north=b_end(k) mg(k)%b_over_a_north=b_end(k)/a_end(k) mg(k)%b_over_a_north_neg=-b_end(k)/a_end(k) mg(k)%half_b_over_a_north_neg=-half*b_end(k)/a_end(k) end do mgrid=min(mgridx,mgridy) do k=1,mgrid mg(k)%half_b_over_a_east_north_neg=mg(k)%half_b_over_a_east_neg+mg(k)%half_b_over_a_north_neg end do write(6,*)' in generate_grids, mgridx,mgridy,mgrid=',mgridx,mgridy,mgrid do k=1,mgrid write(6,'(" k,nx,a_east,a_east_inv,b_east=",i3,i5,3f14.7)') & k,mg(k)%nx,mg(k)%a_east,mg(k)%a_east_inv,mg(k)%b_east end do do k=1,mgrid write(6,'(" k,ny,a_north,a_north_inv,b_north=",i3,i5,3f14.7)') & k,mg(k)%ny,mg(k)%a_north,mg(k)%a_north_inv,mg(k)%b_north end do end subroutine generate_grids subroutine generate_1d_grid(nxf,nxall,mgridx,a_east,b_east) use constants, only: zero,one integer(i_kind),intent(in):: nxf integer(i_kind),intent(out):: mgridx integer(i_kind),intent(out):: nxall(100) real(r_kind),intent(out):: a_east(100),b_east(100) integer(i_kind) ifactor,k,next real(r_kind) end,thisend,factor nxall(1)=nxf end=nxall(1) k=1 ifactor=1 a_east(1)=one b_east(1)=zero do next=nxall(k)/2 ifactor=2*ifactor if(next*ifactor < nxall(1)) next=next+1 if(next < 3) exit factor=ifactor thisend=next*factor k=k+1 nxall(k)=next b_east(k)=(thisend-end)/factor a_east(k)=one-b_east(k) end do mgridx=k end subroutine generate_1d_grid subroutine check_xbdy_err use constants, only: zero,one,two real(r_kind) factor,xend,errmax,xend_interp integer(i_kind) k factor=one xend=mg(1)%nx errmax=zero do k=1,mgridx xend_interp=factor*mg(k)%nx*mg(k)%a_east+factor*(mg(k)%nx-one)*mg(k)%b_east errmax=max(abs(xend_interp-xend),errmax) factor=two*factor end do write(6,*)' for nx=',mg(1)%nx,', mgridx,mgrid,east errmax=',mgridx,mgrid,errmax end subroutine check_xbdy_err subroutine check_ybdy_err use constants, only: zero,one,two real(r_kind) factor,yend,errmax,yend_interp integer(i_kind) k factor=one yend=mg(1)%ny errmax=zero do k=1,mgridy yend_interp=factor*mg(k)%ny*mg(k)%a_north+factor*(mg(k)%ny-one)*mg(k)%b_north errmax=max(abs(yend_interp-yend),errmax) factor=two*factor end do write(6,*)' for ny=',mg(1)%ny,', mgridy,mgrid,north errmax=',mgridy,mgrid,errmax end subroutine check_ybdy_err subroutine init_mgvars(region_dx,region_dy,region_lat,phi0) use constants, only: zero,half,one,two,omega real(r_kind),dimension(0:nya-1,0:nxa-1),intent(in):: region_dx,region_dy,region_lat real(r_kind),intent(in):: phi0 integer(i_kind) i,ia,j,ja,nx,ny integer(i_kind) is1,ie1,js1,je1,is2,ie2,js2,je2 integer(i_kind) k,jump integer(i_kind) ii,im,ip,jj,jm,jp real(r_kind) x00,y00,x0p,x0m,xp0,xm0,y0p,y0m,yp0,ym0 real(r_kind),dimension(0:nya-1,0:nxa-1):: dx_fine,dy_fine,coriolis_fine real(r_kind),allocatable:: d(:),dh(:),u(:),o(:) real(r_kind),allocatable:: dbar(:),ubar(:),obar(:) real(r_kind),allocatable:: dbarh(:),ubarh(:),obarh(:) real(r_kind),allocatable::vtest(:),ftest(:),vtest0(:),ftest0(:) real(r_kind) errmaxv,errmaxvh,errmaxf,errmaxfh integer(i_kind) ibad,jbad integer(i_kind) ic,if,jc,jf nx=mg(1)%nx ; ny=mg(1)%ny if(nx /= nxa-1.or.ny /= nya-1) then write(6,*)' inconsistent args in init_mgvars, nx,nxa-1,ny,nya-1=',nx,nxa-1,ny,nya-1 stop end if do i=0,nx do j=0,ny dx_fine(j,i)=region_dx(j,i) dy_fine(j,i)=region_dy(j,i) coriolis_fine(j,i)=two*omega*sin(region_lat(j,i)) end do end do ! allocate grids at all levels do k=1,mgrid nx=mg(k)%nx ; ny=mg(k)%ny mg(k)%phi0=phi0 allocate(mg(k)%f(0:ny,0:nx),mg(k)%v(0:ny,0:nx),mg(k)%w(0:ny,0:nx)) allocate(mg(k)%w2(0:ny,0:nx)) allocate(mg(k)%y(0:ny,0:nx),mg(k)%z(0:ny,0:nx)) allocate(mg(k)%yy(0:ny,0:nx),mg(k)%zz(0:ny,0:nx)) allocate(mg(k)%x00(0:ny,0:nx),mg(k)%y00(0:ny,0:nx)) allocate(mg(k)%h00(0:ny,0:nx),mg(k)%r00(0:ny,0:nx)) allocate(mg(k)%h0m(0:ny,0:nx),mg(k)%h0p(0:ny,0:nx)) allocate(mg(k)%hm0(0:ny,0:nx),mg(k)%hp0(0:ny,0:nx)) allocate(mg(k)%dbarx(0:ny,0:nx),mg(k)%dbary(0:ny,0:nx)) allocate(mg(k)%ubarx(0:ny,0:nx),mg(k)%ubary(0:ny,0:nx)) allocate(mg(k)%obarx(0:ny,0:nx),mg(k)%obary(0:ny,0:nx)) allocate(mg(k)%dbarxh(0:ny,0:nx),mg(k)%dbaryh(0:ny,0:nx)) allocate(mg(k)%ubarxh(0:ny,0:nx),mg(k)%ubaryh(0:ny,0:nx)) allocate(mg(k)%obarxh(0:ny,0:nx),mg(k)%obaryh(0:ny,0:nx)) allocate(mg(k)%bicubic_mask(0:ny,0:nx)) mg(k)%f=zero ; mg(k)%v=zero ; mg(k)%w=zero ; mg(k)%y=zero ; mg(k)%z=zero mg(k)%x00=zero ; mg(k)%y00=zero mg(k)%h00=zero ; mg(k)%r00=zero mg(k)%h0m=zero ; mg(k)%h0p=zero ; mg(k)%hm0=zero ; mg(k)%hp0=zero mg(k)%dbarx=zero ; mg(k)%dbary=zero mg(k)%ubarx=zero ; mg(k)%ubary=zero mg(k)%obarx=zero ; mg(k)%obary=zero mg(k)%dbarxh=zero ; mg(k)%dbaryh=zero mg(k)%ubarxh=zero ; mg(k)%ubaryh=zero mg(k)%obarxh=zero ; mg(k)%obaryh=zero mg(k)%bicubic_mask=.false. end do ! create x00, y00 and r00 on each grid do k=1,mgrid nx=nxall(k) ; ny=nyall(k) jump=2**(k-1) ibad=0 jbad=0 i=-1 do ii=0,nxall(1),jump i=i+1 j=-1 do jj=0,nyall(1),jump j=j+1 mg(k)%x00(j,i)=dx_fine(jj,ii)*jump mg(k)%y00(j,i)=dy_fine(jj,ii)*jump mg(k)%r00(j,i)=coriolis_fine(jj,ii)**2/mg(k)%phi0 end do if(j < nyall(k)-1.or.j > nyall(k)) jbad=jbad+1 if(j == nyall(k)-1) then j=nyall(k) jj=nyall(1) mg(k)%x00(j,i)=dx_fine(jj,ii)*jump mg(k)%y00(j,i)=dy_fine(jj,ii)*jump mg(k)%r00(j,i)=coriolis_fine(jj,ii)**2/mg(k)%phi0 end if end do if(i < nxall(k)-1.or.i > nxall(k)) ibad=ibad+1 if(i == nxall(k)-1) then i=nxall(k) ii=nxall(1) j=-1 do jj=0,nyall(1),jump j=j+1 mg(k)%x00(j,i)=dx_fine(jj,ii)*jump mg(k)%y00(j,i)=dy_fine(jj,ii)*jump mg(k)%r00(j,i)=coriolis_fine(jj,ii)**2/mg(k)%phi0 end do if(j < nyall(k)-1) jbad=jbad+1 if(j > nyall(k)) jbad=jbad+1 if(j == nyall(k)-1) then j=nyall(k) jj=nyall(1) mg(k)%x00(j,i)=dx_fine(jj,ii)*jump mg(k)%y00(j,i)=dy_fine(jj,ii)*jump mg(k)%r00(j,i)=coriolis_fine(jj,ii)**2/mg(k)%phi0 end if end if write(6,*)' in init_mgvars, k,ibad,jbad=',k,ibad,jbad end do ! next create forward operators do k=1,mgrid nx=nxall(k) ; ny=nyall(k) do i=0,nx ip=min(i+1,nx) ; im=max(i-1,0) do j=0,ny jp=min(j+1,ny) ; jm=max(j-1,0) x00=mg(k)%x00(j,i) y00=mg(k)%y00(j,i) x0p=half*(x00+mg(k)%x00(j ,ip)) x0m=half*(x00+mg(k)%x00(j ,im)) xp0=half*(x00+mg(k)%x00(jp,i )) xm0=half*(x00+mg(k)%x00(jm,i )) y0p=half*(y00+mg(k)%y00(j ,ip)) y0m=half*(y00+mg(k)%y00(j ,im)) yp0=half*(y00+mg(k)%y00(jp,i )) ym0=half*(y00+mg(k)%y00(jm,i )) mg(k)%h0m(j,i)=y0m/(x0m*x00*y00) mg(k)%h0p(j,i)=y0p/(x0p*x00*y00) mg(k)%hm0(j,i)=xm0/(ym0*x00*y00) mg(k)%hp0(j,i)=xp0/(yp0*x00*y00) mg(k)%h00(j,i)=mg(k)%h0m(j,i)+mg(k)%h0p(j,i)+mg(k)%hm0(j,i)+mg(k)%hp0(j,i) end do end do end do ! next create implicit x operators do k=1,mgrid nx=nxall(k) ; ny=nyall(k) allocate(d(nx-1),dh(nx-1),u(nx-2),o(nx-2)) allocate(dbar(nx-1),ubar(nx-2),obar(nx-2)) allocate(dbarh(nx-1),ubarh(nx-2),obarh(nx-2)) ! dirichlet boundary conditions do j=1,ny-1 ii=0 do i=1,nx-1 ii=ii+1 d(ii)=mg(k)%h00(j,i) if(i == nx-1) d(ii)=d(ii)+mg(k)%h0p(j,i)*mg(k)%b_over_a_east dh(ii)=d(ii)+mg(k)%r00(j,i) if(i < nx-1) o(ii)=-mg(k)%h0p(j,i) if(i > 1) u(ii-1)=-mg(k)%h0m(j,i) end do call multi_factorize(d,u,o,nx-1,dbar,ubar,obar) call multi_factorize(dh,u,o,nx-1,dbarh,ubarh,obarh) ii=0 do i=1,nx-1 ii=ii+1 mg(k)%dbarx(j,i)=dbar(ii) mg(k)%dbarxh(j,i)=dbarh(ii) end do ii=0 do i=1,nx-2 ii=ii+1 mg(k)%ubarx(j,i)=ubar(ii) mg(k)%obarx(j,i)=obar(ii) mg(k)%ubarxh(j,i)=ubarh(ii) mg(k)%obarxh(j,i)=obarh(ii) end do end do deallocate(d,dh,u,o,dbar,ubar,obar,dbarh,ubarh,obarh) end do ! next create implicit y operators do k=1,mgrid nx=nxall(k) ; ny=nyall(k) allocate(d(ny-1),dh(ny-1),u(ny-2),o(ny-2)) allocate(dbar(ny-1),ubar(ny-2),obar(ny-2)) allocate(dbarh(ny-1),ubarh(ny-2),obarh(ny-2)) ! dirichlet boundary conditions do i=1,nx-1 jj=0 do j=1,ny-1 jj=jj+1 d(jj)=mg(k)%h00(j,i) if(j == ny-1) d(jj)=d(jj)+mg(k)%hp0(j,i)*mg(k)%b_north/mg(k)%a_north dh(jj)=d(jj)+mg(k)%r00(j,i) if(j < ny-1) o(jj)=-mg(k)%hp0(j,i) if(j > 1) u(jj-1)=-mg(k)%hm0(j,i) end do call multi_factorize(d,u,o,ny-1,dbar,ubar,obar) call multi_factorize(dh,u,o,ny-1,dbarh,ubarh,obarh) jj=0 do j=1,ny-1 jj=jj+1 mg(k)%dbary(j,i)=dbar(jj) mg(k)%dbaryh(j,i)=dbarh(jj) end do jj=0 do j=1,ny-2 jj=jj+1 mg(k)%ubary(j,i)=ubar(jj) mg(k)%obary(j,i)=obar(jj) mg(k)%ubaryh(j,i)=ubarh(jj) mg(k)%obaryh(j,i)=obarh(jj) end do end do deallocate(d,dh,u,o,dbar,ubar,obar,dbarh,ubarh,obarh) end do ! set bicubic_mask do k=2,mgrid do ic=1,mg(k)%nx-1 if=2*ic do jf=3,mg(k-1)%ny-3,2 mg(k-1)%bicubic_mask(jf,if)=.true. end do end do do if=3,mg(k-1)%nx-3,2 do jc=1,mg(k)%ny-1 jf=2*jc mg(k-1)%bicubic_mask(jf,if)=.true. end do end do do if=3,mg(k-1)%nx-3,2 do jf=3,mg(k-1)%ny-3,2 mg(k-1)%bicubic_mask(jf,if)=.true. end do end do end do end subroutine init_mgvars subroutine eval_bdys(f,mgk,nx,ny) ! evaluate boundary point values for field f use constants, only: zero,half,one integer(i_kind),intent(in):: nx,ny real(r_kind),intent(inout):: f(0:ny,0:nx) type(mgvar):: mgk integer(i_kind) i,j if(mgk%nx /= nx.or.mgk%ny /= ny) then write(6,*)' argument error in call to eval_ne_bdys, mgk%nx,nx,mgk%ny,ny=',mgk%nx,nx,mgk%ny,ny stop end if ! north and south boundaries do i=1,nx-1 f(ny,i)=mgk%b_over_a_north_neg*f(ny-1,i) f(0,i)=zero end do ! east and west boundaries do j=1,ny-1 f(j,nx)=mgk%b_over_a_east_neg*f(j,nx-1) f(j,0) =zero end do ! corners--just average of adjacent points on each boundary f( 0, 0)=zero f(ny, 0)=mgk%half_b_over_a_north_neg*f(ny-1,1) f( 0,nx)=mgk%half_b_over_a_east_neg*f(1,nx-1) f(ny,nx)=f(ny-1,nx-1)*mgk%half_b_over_a_east_north_neg end subroutine eval_bdys subroutine eval_bdys_ad(f,mgk,nx,ny) ! evaluate boundary point values for field f use constants, only: zero,half,one integer(i_kind),intent(in):: nx,ny real(r_kind),intent(inout):: f(0:ny,0:nx) type(mgvar):: mgk integer(i_kind) i,j if(mgk%nx /= nx.or.mgk%ny /= ny) then write(6,*)' argument error in call to eval_ne_bdys, mgk%nx,nx,mgk%ny,ny=',mgk%nx,nx,mgk%ny,ny stop end if ! adjoint of corners--just average of adjacent points on each boundary f(ny-1,nx-1)=f(ny-1,nx-1)+f(ny,nx)*mgk%half_b_over_a_east_north_neg f(ny,nx)=zero f(1,nx-1)=f(1,nx-1)+f(0,nx)*mgk%half_b_over_a_east_neg f(0,nx)=zero f(ny-1,1)=f(ny-1,1)+f(ny,0)*mgk%half_b_over_a_north_neg f(ny,0)=zero f(0,0)=zero ! adjoint of east and west boundaries do j=1,ny-1 f(j,0) =zero f(j,nx-1)=f(j,nx-1)+f(j,nx)*mgk%b_over_a_east_neg f(j,nx)=zero end do ! adjoint of north and south boundaries do i=1,nx-1 f(0,i)=zero f(ny-1,i)=f(ny-1,i)+f(ny,i)*mgk%b_over_a_north_neg f(ny,i)=zero end do end subroutine eval_bdys_ad !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! following is section containing transfer operators! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! NOTES ABOUT TRANSFER OPERATORS: ! ! c ! I transfer operator: transfer from fine grid to coarse grid using full weighting ! f ! implemented by subroutine f2c_full_weight ! ! the full weighting stencil is ! ! 1/16 1/8 1/16 ! ! 1/8 1/4 1/8 ! ! 1/16 1/8 1/16 ! ! where the central point coincides with a coarse grid point, and the surrounding ! points are fine grid points in between the coarse points. ! ! f ! I transfer operator: transfer from coarse grid to fine grid using bilinear interpolation ! c ! implemented by subroutine c2f_bilinear ! ! f ! II transfer operator: transfer from coarse grid to fine grid using bicubic interpolation ! c ! implemented by subroutine c2f_bicubic ! ! ! NOTE 1: all grids in this multigrid system have common point (0,0). ! the grid dimensions are f(0:nyf,0:nxf), c(0:nyc,0:nxc) ! ! ! NOTE 2: the coarse grid increment is exactly twice the fine grid increment, so ! every other point of the fine grid moving away from 0, coincides with ! a coarse grid point. ! ! NOTE 3: the east (nxc) and north (nyc) edges of the coarse grid are defined using ! boundary interpolation weights, since in general they do not coincide with ! the boundary of the finest grid (this is because there is no condition on ! the number of points, so grid coarsening can result in requiring an extra ! row beyond east/north edges of the domain. ! ! NOTE 4: weights used here correspond to adjoint of bilinear interpolation from coarse to fine ! with the exception that the weights are normalized to preserve area integral. ! ! ! subroutine f2c_full_weight(f,c,nxf,nyf,nxc,nyc) ! full weighting fine to coarse transfer similar to adjoint of coarse to fine bilinear interpolation ! except coefficients are normalized by their sum and only points interior to fine grid ! are updated. use kinds, only: r_kind,i_kind use constants, only: zero implicit none integer(i_kind),intent(in):: nxc,nyc,nxf,nyf real(r_kind),intent(in):: f(0:nyf,0:nxf) real(r_kind),intent(out):: c(0:nyc,0:nxc) real(r_kind),parameter:: one_quarter=1._8/4._8,one_eighth=1._8/8._8,one_sixteenth=1._8/16._8 integer(i_kind) if,ic,jf,jc,jfp,jfm,ifp,ifm c=zero do ic=1,nxc-1 if=2*ic ifp=if+1 ifm=if-1 do jc=1,nyc-1 jf=2*jc jfp=jf+1 jfm=jf-1 c(jc,ic)=one_quarter* f(jf ,if ) & +one_eighth*(f(jf ,ifp)+f(jf ,ifm)+f(jfp, if)+f(jfm, if)) & +one_sixteenth*(f(jfm,ifm)+f(jfm,ifp)+f(jfp,ifm)+f(jfp,ifp)) end do end do end subroutine f2c_full_weight subroutine f2c_full_weight_ad(f,c,nxf,nyf,nxc,nyc) ! full weighting fine to coarse transfer similar to adjoint of coarse to fine bilinear interpolation ! except coefficients are normalized by their sum and only points interior to fine grid ! are updated. use kinds, only: r_kind,i_kind use constants, only: zero implicit none integer(i_kind),intent(in):: nxc,nyc,nxf,nyf real(r_kind),intent(out):: f(0:nyf,0:nxf) real(r_kind),intent(in):: c(0:nyc,0:nxc) real(r_kind),parameter:: one_quarter=1._8/4._8,one_eighth=1._8/8._8,one_sixteenth=1._8/16._8 integer(i_kind) if,ic,jf,jc,jfp,jfm,ifp,ifm f=zero do ic=1,nxc-1 if=2*ic ifp=if+1 ifm=if-1 do jc=1,nyc-1 jf=2*jc jfp=jf+1 jfm=jf-1 f(jf ,if )=f(jf ,if ) + one_quarter *c(jc,ic) f(jf ,ifp)=f(jf ,ifp) + one_eighth *c(jc,ic) f(jf ,ifm)=f(jf ,ifm) + one_eighth *c(jc,ic) f(jfp,if )=f(jfp,if ) + one_eighth *c(jc,ic) f(jfm,if )=f(jfm,if ) + one_eighth *c(jc,ic) f(jfm,ifm)=f(jfm,ifm) + one_sixteenth*c(jc,ic) f(jfm,ifp)=f(jfm,ifp) + one_sixteenth*c(jc,ic) f(jfp,ifm)=f(jfp,ifm) + one_sixteenth*c(jc,ic) f(jfp,ifp)=f(jfp,ifp) + one_sixteenth*c(jc,ic) end do end do end subroutine f2c_full_weight_ad subroutine c2f_bilinear(c,f,nxc,nyc,nxf,nyf) ! bilinear interpolation from coarse to fine grid ! note: c boundary should be filled in to get correct results on 1st interior row of f use kinds, only: r_kind,i_kind use constants, only: zero,quarter,half implicit none integer(i_kind),intent(in):: nxc,nyc,nxf,nyf real(r_kind),intent(in):: c(0:nyc,0:nxc) real(r_kind),intent(out):: f(0:nyf,0:nxf) integer(i_kind) if,ic,jf,jc,icm,icp,jcm,jcp ! fill in coincident points: f=zero do ic=1,nxc-1 if=2*ic do jc=1,nyc-1 jf=2*jc f(jf,if)=c(jc,ic) end do end do ! fill in y blanks on coincident x lines do ic=1,nxc-1 if=2*ic do jf=1,nyf-1,2 jcm=(jf-1)/2 jcp=(jf+1)/2 f(jf,if)=half*(c(jcm,ic)+c(jcp,ic)) end do end do ! fill in x blanks on coincident y lines do if=1,nxf-1,2 icm=(if-1)/2 icp=(if+1)/2 do jc=1,nyc-1 jf=2*jc f(jf,if)=half*(c(jc,icm)+c(jc,icp)) end do end do ! fill in remaining holes do if=1,nxf-1,2 icm=(if-1)/2 icp=(if+1)/2 do jf=1,nyf-1,2 jcm=(jf-1)/2 jcp=(jf+1)/2 f(jf,if)=quarter*(c(jcm,icm)+c(jcm,icp)+c(jcp,icm)+c(jcp,icp)) end do end do end subroutine c2f_bilinear subroutine c2f_bilinear_ad(c,f,nxc,nyc,nxf,nyf) ! adjoint of c2f_bilinear use kinds, only: r_kind,i_kind use constants, only: zero,quarter,half integer(i_kind),intent(in):: nxc,nyc,nxf,nyf real(r_kind),intent(out):: c(0:nyc,0:nxc) real(r_kind),intent(in):: f(0:nyf,0:nxf) integer(i_kind) if,ic,jf,jc,icm,icp,jcm,jcp c=zero ! adjoint of fill in remaining holes do if=1,nxf-1,2 icm=(if-1)/2 icp=(if+1)/2 do jf=1,nyf-1,2 jcm=(jf-1)/2 jcp=(jf+1)/2 c(jcm,icm)=c(jcm,icm)+quarter*f(jf,if) c(jcm,icp)=c(jcm,icp)+quarter*f(jf,if) c(jcp,icm)=c(jcp,icm)+quarter*f(jf,if) c(jcp,icp)=c(jcp,icp)+quarter*f(jf,if) end do end do ! adjoint of fill in x blanks on coincident y lines do if=1,nxf-1,2 icm=(if-1)/2 icp=(if+1)/2 do jc=1,nyc-1 jf=2*jc c(jc,icm)=c(jc,icm)+half*f(jf,if) c(jc,icp)=c(jc,icp)+half*f(jf,if) end do end do ! adjoint of fill in y blanks on coincident x lines do ic=1,nxc-1 if=2*ic do jf=1,nyf-1,2 jcm=(jf-1)/2 jcp=(jf+1)/2 c(jcm,ic)=c(jcm,ic)+half*f(jf,if) c(jcp,ic)=c(jcp,ic)+half*f(jf,if) end do end do ! adjoint of fill in coincident points: do ic=1,nxc-1 if=2*ic do jc=1,nyc-1 jf=2*jc c(jc,ic)=c(jc,ic)+f(jf,if) end do end do end subroutine c2f_bilinear_ad subroutine c2f_bicubic(c,f,nxc,nyc,nxf,nyf,bicubic_mask) ! bicubic interpolation from coarse to fine grid use kinds, only: r_kind,i_kind use constants, only: zero,quarter,half implicit none integer(i_kind),intent(in):: nxc,nyc,nxf,nyf real(r_kind),intent(in):: c(0:nyc,0:nxc) real(r_kind),intent(out):: f(0:nyf,0:nxf) logical,intent(in):: bicubic_mask(0:nyf,0:nxf) real(r_kind),parameter:: one_sixteenth=1._8/16._8,nine_sixteenths=9._8/16._8 real(r_kind),parameter:: one_sixt2=one_sixteenth**2,onenine_sixt2=one_sixteenth*nine_sixteenths real(r_kind),parameter:: nine_sixt2=nine_sixteenths**2 integer(i_kind) if,ic,jf,jc integer(i_kind) icm2,icp,icp2,icm,jcm2,jcm,jcp2,jcp,nxflim,nyflim ! start with coincident points and points closest to boundary that are linearly interpolated: ! fill in coincident points: f=zero do ic=1,nxc-1 if=2*ic do jc=1,nyc-1 jf=2*jc if(.not.bicubic_mask(jf,if)) f(jf,if)=c(jc,ic) end do end do ! fill in y blanks with linear interpolation on coincident x lines adjacent to y boundarys do ic=1,nxc-1 if=2*ic do jf=1,nyf-1,2 jcm=(jf-1)/2 jcp=(jf+1)/2 if(.not.bicubic_mask(jf,if)) f(jf,if)=half*(c(jcm,ic)+c(jcp,ic)) end do end do ! fill in x blanks on coincident y lines do if=1,nxf-1,2 icm=(if-1)/2 icp=(if+1)/2 do jc=1,nyc-1 jf=2*jc if(.not.bicubic_mask(jf,if)) f(jf,if)=half*(c(jc,icm)+c(jc,icp)) end do end do ! fill in remaining holes do if=1,nxf-1,2 icm=(if-1)/2 icp=(if+1)/2 do jf=1,nyf-1,2 jcm=(jf-1)/2 jcp=(jf+1)/2 if(.not.bicubic_mask(jf,if)) f(jf,if)=quarter*(c(jcm,icm)+c(jcm,icp)+c(jcp,icm)+c(jcp,icp)) end do end do ! fill in y blanks on coincident interior x lines do ic=1,nxc-1 if=2*ic do jf=3,nyf-3,2 jcm=(jf-1)/2 ; jcm2=(jf-3)/2 ; jcp=(jf+1)/2 ; jcp2=(jf+3)/2 f(jf,if)=nine_sixteenths*(c(jcm,ic)+c(jcp,ic))-one_sixteenth*(c(jcm2,ic)+c(jcp2,ic)) end do end do ! fill in x blanks on coincident interior y lines do if=3,nxf-3,2 icm=(if-1)/2 ; icm2=(if-3)/2 ; icp=(if+1)/2 ; icp2=(if+3)/2 do jc=1,nyc-1 jf=2*jc f(jf,if)=nine_sixteenths*(c(jc,icm)+c(jc,icp))-one_sixteenth*(c(jc,icm2)+c(jc,icp2)) end do end do ! fill in remaining interior holes do if=3,nxf-3,2 icm=(if-1)/2 ; icm2=(if-3)/2 ; icp=(if+1)/2 ; icp2=(if+3)/2 do jf=3,nyf-3,2 jcm=(jf-1)/2 ; jcm2=(jf-3)/2 ; jcp=(jf+1)/2 ; jcp2=(jf+3)/2 f(jf,if)=nine_sixt2*(c(jcm ,icm )+c(jcm ,icp )+c(jcp ,icm )+c(jcp ,icp )) & -onenine_sixt2*(c(jcm ,icm2)+c(jcm ,icp2)+c(jcp ,icm2)+c(jcp ,icp2) & +c(jcm2,icm )+c(jcp2,icm )+c(jcm2,icp )+c(jcp2,icp )) & +one_sixt2*(c(jcm2,icm2)+c(jcm2,icp2)+c(jcp2,icm2)+c(jcp2,icp2)) end do end do end subroutine c2f_bicubic subroutine c2f_bicubic_ad(c,f,nxc,nyc,nxf,nyf,bicubic_mask) ! bicubic interpolation from coarse to fine grid use kinds, only: r_kind,i_kind use constants, only: zero,quarter,half implicit none integer(i_kind),intent(in):: nxc,nyc,nxf,nyf real(r_kind),intent(out):: c(0:nyc,0:nxc) real(r_kind),intent(in):: f(0:nyf,0:nxf) logical,intent(in):: bicubic_mask(0:nyf,0:nxf) real(r_kind),parameter:: one_sixteenth=1._8/16._8,nine_sixteenths=9._8/16._8 real(r_kind),parameter:: one_sixt2=one_sixteenth**2,onenine_sixt2=one_sixteenth*nine_sixteenths real(r_kind),parameter:: nine_sixt2=nine_sixteenths**2 integer(i_kind) if,ic,jf,jc integer(i_kind) icm2,icp,icp2,icm,jcm2,jcm,jcp2,jcp,nxflim,nyflim c=zero ! adjoint of fill in remaining interior holes do if=3,nxf-3,2 icm=(if-1)/2 ; icm2=(if-3)/2 ; icp=(if+1)/2 ; icp2=(if+3)/2 do jf=3,nyf-3,2 jcm=(jf-1)/2 ; jcm2=(jf-3)/2 ; jcp=(jf+1)/2 ; jcp2=(jf+3)/2 c(jcm ,icm )=c(jcm ,icm )+f(jf,if)*nine_sixt2 c(jcm ,icp )=c(jcm ,icp )+f(jf,if)*nine_sixt2 c(jcp ,icm )=c(jcp ,icm )+f(jf,if)*nine_sixt2 c(jcp ,icp )=c(jcp ,icp )+f(jf,if)*nine_sixt2 c(jcm ,icm2)=c(jcm ,icm2)-f(jf,if)*onenine_sixt2 c(jcm ,icp2)=c(jcm ,icp2)-f(jf,if)*onenine_sixt2 c(jcp ,icm2)=c(jcp ,icm2)-f(jf,if)*onenine_sixt2 c(jcp ,icp2)=c(jcp ,icp2)-f(jf,if)*onenine_sixt2 c(jcm2,icm )=c(jcm2,icm )-f(jf,if)*onenine_sixt2 c(jcp2,icm )=c(jcp2,icm )-f(jf,if)*onenine_sixt2 c(jcm2,icp )=c(jcm2,icp )-f(jf,if)*onenine_sixt2 c(jcp2,icp )=c(jcp2,icp )-f(jf,if)*onenine_sixt2 c(jcm2,icm2)=c(jcm2,icm2)+f(jf,if)*one_sixt2 c(jcm2,icp2)=c(jcm2,icp2)+f(jf,if)*one_sixt2 c(jcp2,icm2)=c(jcp2,icm2)+f(jf,if)*one_sixt2 c(jcp2,icp2)=c(jcp2,icp2)+f(jf,if)*one_sixt2 end do end do ! adjoint of fill in x blanks on coincident interior y lines do if=3,nxf-3,2 icm=(if-1)/2 ; icm2=(if-3)/2 ; icp=(if+1)/2 ; icp2=(if+3)/2 do jc=1,nyc-1 jf=2*jc c(jc,icm )=c(jc,icm )+f(jf,if)*nine_sixteenths c(jc,icp )=c(jc,icp )+f(jf,if)*nine_sixteenths c(jc,icm2)=c(jc,icm2)-f(jf,if)*one_sixteenth c(jc,icp2)=c(jc,icp2)-f(jf,if)*one_sixteenth end do end do ! adjoint of fill in y blanks on coincident interior x lines do ic=1,nxc-1 if=2*ic do jf=3,nyf-3,2 jcm=(jf-1)/2 ; jcm2=(jf-3)/2 ; jcp=(jf+1)/2 ; jcp2=(jf+3)/2 c(jcm ,ic)=c(jcm ,ic)+f(jf,if)*nine_sixteenths c(jcp ,ic)=c(jcp ,ic)+f(jf,if)*nine_sixteenths c(jcm2,ic)=c(jcm2,ic)-f(jf,if)*one_sixteenth c(jcp2,ic)=c(jcp2,ic)-f(jf,if)*one_sixteenth end do end do ! adjoint of start with coincident points and points closest to boundary that are linearly interpolated: ! adjoint of fill in remaining holes do if=1,nxf-1,2 icm=(if-1)/2 icp=(if+1)/2 do jf=1,nyf-1,2 jcm=(jf-1)/2 jcp=(jf+1)/2 if(.not.bicubic_mask(jf,if)) then c(jcm,icm)=c(jcm,icm)+f(jf,if)*quarter c(jcm,icp)=c(jcm,icp)+f(jf,if)*quarter c(jcp,icm)=c(jcp,icm)+f(jf,if)*quarter c(jcp,icp)=c(jcp,icp)+f(jf,if)*quarter end if end do end do ! adjoint of fill in x blanks on coincident y lines do if=1,nxf-1,2 icm=(if-1)/2 icp=(if+1)/2 do jc=1,nyc-1 jf=2*jc if(.not.bicubic_mask(jf,if)) then c(jc,icm)=c(jc,icm)+f(jf,if)*half c(jc,icp)=c(jc,icp)+f(jf,if)*half end if end do end do ! adjoint of fill in y blanks with linear interpolation on coincident x lines adjacent to y boundarys do ic=1,nxc-1 if=2*ic do jf=1,nyf-1,2 jcm=(jf-1)/2 jcp=(jf+1)/2 if(.not.bicubic_mask(jf,if)) then c(jcm,ic)=c(jcm,ic)+f(jf,if)*half c(jcp,ic)=c(jcp,ic)+f(jf,if)*half end if end do end do ! adjoint of fill in coincident points: do ic=1,nxc-1 if=2*ic do jc=1,nyc-1 jf=2*jc if(.not.bicubic_mask(jf,if)) then c(jc,ic)=c(jc,ic)+f(jf,if) end if end do end do end subroutine c2f_bicubic_ad subroutine multi_factorize(d,u,o,nx,dbar,ubar,obar) ! factorize special input matrix into lower times upper triangular matrices (illustrated for nx=5): ! _ _ _ ! |d(1) o(1) 0 0 0 | |d(1) 0 0 0 0 | |d(1) o(1) 0 0 0 | ! | | |_ _ | | _ _ | ! |u(1) d(2) o(2) 0 0 | |u(1) d(2) 0 0 0 | |0 d(2) o(2) 0 0 | ! | | | _ _ | | _ _ | ! |0 u(2) d(3) o(3) 0 | = |0 u(2) d(2) 0 0 |*|0 0 d(3) o(3) 0 | ! | | | _ _ | | _ _ | ! |0 0 u(3) d(4) o(4)| |0 0 u(2) d(4) 0 | |0 0 0 d(4) o(4)| ! | | | _ _ | | _ | ! |0 0 0 u(4) d(5)| |0 0 0 u(4) d(5)| |0 0 0 0 d(5)| integer(i_kind),intent(in):: nx real(r_kind),intent(in):: d(nx),u(nx-1),o(nx-1) real(r_kind),intent(out):: dbar(nx),ubar(nx-1),obar(nx-1) integer(i_kind) i dbar(1)=sqrt(d(1)) obar(1)=o(1)/dbar(1) ubar(1)=u(1)/dbar(1) do i=2,nx-1 dbar(i)=sqrt(d(i)-ubar(i-1)*obar(i-1)) obar(i)=o(i)/dbar(i) ubar(i)=u(i)/dbar(i) end do i=nx dbar(i)=sqrt(d(i)-ubar(i-1)*obar(i-1)) end subroutine multi_factorize subroutine forward_op_x1(v,f,d,u,o,nx) use constants, only: zero integer(i_kind),intent(in)::nx real(r_kind),intent(in):: v(nx),d(nx),u(nx-1),o(nx-1) real(r_kind),intent(out)::f(nx) integer(i_kind) i f(1)=d(1)*v(1)+o(1)*v(2) do i=2,nx-1 f(i)=u(i-1)*v(i-1)+d(i)*v(i)+o(i)*v(i+1) end do f(nx)=u(nx-1)*v(nx-1)+d(nx)*v(nx) end subroutine forward_op_x1 subroutine inverse_op_x1(f,v,dbar,ubar,obar,nx) use constants, only: zero integer(i_kind),intent(in)::nx real(r_kind),intent(in):: f(nx),dbar(nx),ubar(nx-1),obar(nx-1) real(r_kind),intent(out)::v(nx) integer(i_kind) i real(r_kind) temp(nx) temp(1)=f(1)/dbar(1) do i=2,nx temp(i)=(f(i)-ubar(i-1)*temp(i-1))/dbar(i) end do v(nx)=temp(nx)/dbar(nx) do i=nx-1,1,-1 v(i)=(temp(i)-obar(i)*v(i+1))/dbar(i) end do end subroutine inverse_op_x1 subroutine forward_op_x(v,f,nx,ny,j,k,helmholtz) use constants, only: zero integer(i_kind),intent(in)::nx,ny,j,k real(r_kind),intent(in)::helmholtz real(r_kind),intent(in):: v(0:ny,0:nx) real(r_kind),intent(out)::f(0:ny,0:nx) integer(i_kind) i,im1,ip1 do i=1,nx-1 im1=i-1 ; ip1=i+1 f(j,i)=mg(k)%h0m(j,i)*v(j, im1)+mg(k)%h0p(j,i)*v(j,ip1) & -(mg(k)%h00(j,i)+helmholtz*mg(k)%r00(j,i))*v(j,i) end do end subroutine forward_op_x subroutine inverse_luop_x_1(rhs,v,dbarx,ubarx,obarx,nx,ny,j) integer(i_kind),intent(in):: nx,ny,j real(r_kind),dimension(0:ny,0:nx),intent(in):: rhs,dbarx,ubarx,obarx real(r_kind),dimension(0:ny,0:nx),intent(inout):: v integer(i_kind) i,im1,ip1 real(r_kind) temp(0:ny,0:nx) ! inverse of lower triangle first: temp(j,1)=rhs(j,1)/dbarx(j,1) do i=2,nx-1 im1=i-1 temp(j,i)=(rhs(j,i)-ubarx(j,im1)*temp(j,im1))/dbarx(j,i) end do ! finish with inverse of upper triangle: v(j,nx-1)=temp(j,nx-1)/dbarx(j,nx-1) do i=nx-2,1,-1 ip1=i+1 v(j,i)=(temp(j,i)-obarx(j,i)*v(j,ip1))/dbarx(j,i) end do end subroutine inverse_luop_x_1 subroutine forward_op_y(v,f,nx,ny,i,k,helmholtz) use constants, only: zero integer(i_kind),intent(in)::nx,ny,i,k real(r_kind),intent(in)::helmholtz real(r_kind),intent(in):: v(0:ny,0:nx) real(r_kind),intent(out)::f(0:ny,0:nx) integer(i_kind) j,jm1,jp1 do j=1,ny-1 jm1=j-1 ; jp1=j+1 f(j,i)=mg(k)%hm0(j,i)*v(jm1,i )+mg(k)%hp0(j,i)*v(jp1,i) & -(mg(k)%h00(j,i)+helmholtz*mg(k)%r00(j,i))*v(j,i) end do end subroutine forward_op_y subroutine inverse_luop_y_1(rhs,v,dbary,ubary,obary,nx,ny,i) integer(i_kind),intent(in):: nx,ny,i real(r_kind),dimension(0:ny,0:nx),intent(in):: rhs,dbary,ubary,obary real(r_kind),dimension(0:ny,0:nx),intent(inout):: v integer(i_kind) j real(r_kind) temp(0:ny) ! inverse of lower triangle first: temp(1)=rhs(1,i)/dbary(1,i) do j=2,ny-1 temp(j)=(rhs(j,i)-ubary(j-1,i)*temp(j-1))/dbary(j,i) end do ! finish with inverse of upper triangle: v(ny-1,i)=temp(ny-1)/dbary(ny-1,i) do j=ny-2,1,-1 v(j,i)=(temp(j)-obary(j,i)*v(j+1,i))/dbary(j,i) end do end subroutine inverse_luop_y_1 subroutine forward_op(v,f,nx,ny,k,helmholtz) use constants, only: zero integer(i_kind),intent(in)::nx,ny,k real(r_kind),intent(in)::helmholtz real(r_kind),intent(in):: v(0:ny,0:nx) real(r_kind),intent(out)::f(0:ny,0:nx) integer(i_kind) i,im1,ip1,j ! check that nx,ny are correct values for level k: if(mg(k)%nx /= nx.or.mg(k)%ny /= ny) then write(6,*)' inconsistent input nx,ny for fmg_forward_op, program stops' stop end if f=zero do i=1,nx-1 im1=i-1 ; ip1=i+1 do j=1,ny-1 f(j,i)=mg(k)%h0m(j,i)*v(j, im1)+mg(k)%h0p(j,i)*v(j,ip1) & +mg(k)%hm0(j,i)*v(j-1,i )+mg(k)%hp0(j,i)*v(j+1,i) & -(mg(k)%h00(j,i)+helmholtz*mg(k)%r00(j,i))*v(j,i) end do end do end subroutine forward_op subroutine forward_op_ad(v,f,nx,ny,k,helmholtz) use constants, only: zero integer(i_kind),intent(in)::nx,ny,k real(r_kind),intent(in)::helmholtz real(r_kind),intent(out):: v(0:ny,0:nx) real(r_kind),intent(in)::f(0:ny,0:nx) integer(i_kind) i,im1,ip1,j ! check that nx,ny are correct values for level k: if(mg(k)%nx /= nx.or.mg(k)%ny /= ny) then write(6,*)' inconsistent input nx,ny for fmg_forward_op_ad, program stops' stop end if v=zero do i=1,nx-1 im1=i-1 ; ip1=i+1 do j=1,ny-1 v(j ,im1)=v(j ,im1)+mg(k)%h0m(j,i)*f(j,i) v(j ,ip1)=v(j ,ip1)+mg(k)%h0p(j,i)*f(j,i) v(j-1,i )=v(j-1,i )+mg(k)%hm0(j,i)*f(j,i) v(j+1,i )=v(j+1,i )+mg(k)%hp0(j,i)*f(j,i) v(j ,i )=v(j ,i )-(mg(k)%h00(j,i)+helmholtz*mg(k)%r00(j,i))*f(j,i) end do end do end subroutine forward_op_ad subroutine relax(v,f,dbarx,ubarx,obarx,dbary,ubary,obary,h0p,h0m,hp0,hm0,nx,ny, & b_over_a_east_neg,b_over_a_north_neg) use constants, only: zero ! do one iteration of adi relaxation, first in x direction, then in y direction integer(i_kind),intent(in):: nx,ny real(r_kind),dimension(0:ny,0:nx),intent(inout):: v real(r_kind),dimension(0:ny,0:nx),intent(in):: f real(r_kind),dimension(0:ny,0:nx),intent(in):: dbarx,ubarx,obarx real(r_kind),dimension(0:ny,0:nx),intent(in):: dbary,ubary,obary real(r_kind),dimension(0:ny,0:nx),intent(in):: h0p,h0m,hp0,hm0 real(r_kind),intent(in):: b_over_a_east_neg,b_over_a_north_neg integer(i_kind) i,ii,im1,ip1,j,jj real(r_kind) rhs(0:ny,0:nx) ! do x lines first: ii=1 do jj=1,2 do i=1,nx-1 do j=jj,ny-1,2 rhs(j,i)=-f(j,i)+hp0(j,i)*v(j+1,i)+hm0(j,i)*v(j-1,i) end do end do call inverse_luop_x(rhs,v,dbarx,ubarx,obarx,nx,ny,jj,ii) do j=jj,ny-1,2 v(j,0)=zero v(j,nx)=b_over_a_east_neg*v(j,nx-1) end do end do ! now do y lines: jj=1 do ii=1,2 do i=ii,nx-1,2 ip1=i+1 ; im1=i-1 do j=1,ny-1 rhs(j,i)=-f(j,i)+h0p(j,i)*v(j,ip1)+h0m(j,i)*v(j,im1) end do end do call inverse_luop_y(rhs,v,dbary,ubary,obary,nx,ny,ii,jj) do i=ii,nx-1,2 v(0,i)=zero v(ny,i)=b_over_a_north_neg*v(ny-1,i) end do end do end subroutine relax subroutine relax_ad(v,f,dbarx,ubarx,obarx,dbary,ubary,obary,h0p,h0m,hp0,hm0,nx,ny, & b_over_a_east_neg,b_over_a_north_neg) use constants, only: zero ! do one iteration of adi relaxation, first in x direction, then in y direction integer(i_kind),intent(in):: nx,ny real(r_kind),dimension(0:ny,0:nx),intent(inout):: v real(r_kind),dimension(0:ny,0:nx),intent(inout):: f real(r_kind),dimension(0:ny,0:nx),intent(in):: dbarx,ubarx,obarx real(r_kind),dimension(0:ny,0:nx),intent(in):: dbary,ubary,obary real(r_kind),dimension(0:ny,0:nx),intent(in):: h0p,h0m,hp0,hm0 real(r_kind),intent(in):: b_over_a_east_neg,b_over_a_north_neg integer(i_kind) i,ii,im1,ip1,j,jj real(r_kind) rhs(0:ny,0:nx) ! adjoint of now do y lines: jj=1 do ii=2,1,-1 do i=ii,nx-1,2 v(ny-1,i)=v(ny-1,i)+b_over_a_north_neg*v(ny,i) v(ny,i)=zero v(0,i)=zero end do call inverse_luop_y_ad(rhs,v,dbary,ubary,obary,nx,ny,ii,jj) do i=ii,nx-1,2 ip1=i+1 ; im1=i-1 do j=ny-1,1,-1 f(j,i)=f(j,i)-rhs(j,i) v(j,ip1)=v(j,ip1)+h0p(j,i)*rhs(j,i) v(j,im1)=v(j,im1)+h0m(j,i)*rhs(j,i) v(j,i)=zero end do end do end do ! adjoint of do x lines first: ii=1 do jj=2,1,-1 do j=jj,ny-1,2 v(j,nx-1)=v(j,nx-1)+b_over_a_east_neg*v(j,nx) v(j,nx)=zero v(j,0)=zero end do call inverse_luop_x_ad(rhs,v,dbarx,ubarx,obarx,nx,ny,jj,ii) do i=nx-1,1,-1 do j=jj,ny-1,2 f(j,i)=f(j,i)-rhs(j,i) v(j+1,i)=v(j+1,i)+hp0(j,i)*rhs(j,i) v(j-1,i)=v(j-1,i)+hm0(j,i)*rhs(j,i) v(j,i)=zero end do end do end do end subroutine relax_ad subroutine inverse_luop_x(rhs,v,dbarx,ubarx,obarx,nx,ny,jj,ii) integer(i_kind),intent(in):: nx,ny,jj,ii real(r_kind),dimension(0:ny,0:nx),intent(in):: rhs,dbarx,ubarx,obarx real(r_kind),dimension(0:ny,0:nx),intent(inout):: v integer(i_kind) i,im1,ip1,j real(r_kind) temp(0:ny,0:nx) ! inverse of lower triangle first: do j=jj,ny-1,2 temp(j,ii)=rhs(j,ii)/dbarx(j,ii) end do do i=ii+1,nx-ii im1=i-1 do j=jj,ny-1,2 temp(j,i)=(rhs(j,i)-ubarx(j,im1)*temp(j,im1))/dbarx(j,i) end do end do ! finish with inverse of upper triangle: do j=jj,ny-1,2 v(j,nx-ii)=temp(j,nx-ii)/dbarx(j,nx-ii) end do do i=nx-ii-1,ii,-1 ip1=i+1 do j=jj,ny-1,2 v(j,i)=(temp(j,i)-obarx(j,i)*v(j,ip1))/dbarx(j,i) end do end do end subroutine inverse_luop_x subroutine inverse_luop_x_ad(rhs,v,dbarx,ubarx,obarx,nx,ny,jj,ii) use constants, only: zero integer(i_kind),intent(in):: nx,ny,jj,ii real(r_kind),dimension(0:ny,0:nx),intent(in):: dbarx,ubarx,obarx real(r_kind),dimension(0:ny,0:nx),intent(inout):: rhs real(r_kind),dimension(0:ny,0:nx),intent(inout):: v integer(i_kind) i,im1,ip1,j real(r_kind) temp(0:ny,0:nx) ! adjoint of finish with inverse of upper triangle: do i=ii,nx-ii-1 do j=jj,ny-1,2 temp(j,i)=zero end do end do do j=jj,ny-1,2 temp(j,nx-ii)=zero end do do i=ii,nx-ii-1 ip1=i+1 do j=jj,ny-1,2 temp(j,i)=temp(j,i)+v(j,i)/dbarx(j,i) v(j,ip1)=v(j,ip1)-obarx(j,i)*v(j,i)/dbarx(j,i) v(j,i)=zero !??? end do end do do j=jj,ny-1,2 temp(j,nx-ii)=temp(j,nx-ii)+v(j,nx-ii)/dbarx(j,nx-ii) end do ! adjoint of inverse of lower triangle first: do i=nx-ii,ii+1,-1 do j=jj,ny-1,2 rhs(j,i)=zero end do end do do j=jj,ny-1,2 rhs(j,ii)=zero end do do i=nx-ii,ii+1,-1 im1=i-1 do j=jj,ny-1,2 rhs(j,i)=rhs(j,i)+temp(j,i)/dbarx(j,i) temp(j,im1)=temp(j,im1)-ubarx(j,im1)*temp(j,i)/dbarx(j,i) end do end do do j=jj,ny-1,2 rhs(j,ii)=rhs(j,ii)+temp(j,ii)/dbarx(j,ii) end do end subroutine inverse_luop_x_ad subroutine inverse_luop_y(rhs,v,dbary,ubary,obary,nx,ny,ii,jj) integer(i_kind),intent(in):: nx,ny,ii,jj real(r_kind),dimension(0:ny,0:nx),intent(in):: rhs,dbary,ubary,obary real(r_kind),dimension(0:ny,0:nx),intent(inout):: v integer(i_kind) i,j real(r_kind) temp(0:ny) do i=ii,nx-1,2 ! inverse of lower triangle first: temp(jj)=rhs(jj,i)/dbary(jj,i) do j=jj+1,ny-jj temp(j)=(rhs(j,i)-ubary(j-1,i)*temp(j-1))/dbary(j,i) end do ! finish with inverse of upper triangle: v(ny-jj,i)=temp(ny-jj)/dbary(ny-jj,i) do j=ny-jj-1,jj,-1 v(j,i)=(temp(j)-obary(j,i)*v(j+1,i))/dbary(j,i) end do end do end subroutine inverse_luop_y subroutine inverse_luop_y_ad(rhs,v,dbary,ubary,obary,nx,ny,ii,jj) use constants, only: zero integer(i_kind),intent(in):: nx,ny,ii,jj real(r_kind),dimension(0:ny,0:nx),intent(in):: dbary,ubary,obary real(r_kind),dimension(0:ny,0:nx),intent(inout):: rhs real(r_kind),dimension(0:ny,0:nx),intent(inout):: v integer(i_kind) i,j real(r_kind) temp(0:ny) do i=ii,nx-1,2 ! adjoint of finish with inverse of upper triangle: do j=jj,ny-jj-1 temp(j)=zero end do temp(ny-jj)=zero do j=jj,ny-jj-1 temp(j)=temp(j)+v(j,i)/dbary(j,i) v(j+1,i)=v(j+1,i)-obary(j,i)*v(j,i)/dbary(j,i) v(j,i)=zero end do temp(ny-jj)=temp(ny-jj)+v(ny-jj,i)/dbary(ny-jj,i) ! adjoint of inverse of lower triangle first: do j=ny-jj,jj+1,-1 rhs(j,i)=zero end do rhs(jj,i)=zero do j=ny-jj,jj+1,-1 rhs(j,i)=rhs(j,i)+temp(j)/dbary(j,i) temp(j-1)=temp(j-1)-ubary(j-1,i)*temp(j)/dbary(j,i) end do rhs(jj,i)=rhs(jj,i)+temp(jj)/dbary(jj,i) end do end subroutine inverse_luop_y_ad subroutine fmg(v,f,helmholtz,nrelax1,nrelax2,nrelax_solve,nx,ny) ! full multigrid solution on fixed domain with homogeneous dirichlet boundary conditions ! solve delsqr(v) = f helmholtz=.false. ! (delsqr-R)(v) = f helmholtz=.true. use constants, only: zero,one integer(i_kind),intent(in)::nx,ny,nrelax1,nrelax2,nrelax_solve logical,intent(in):: helmholtz real(r_kind),intent(in):: f(0:ny,0:nx) ! ! in actual use, f is in, and v is out, but real(r_kind),intent(out):: v(0:ny,0:nx) ! ! make both inout for adjoint test purposes integer(i_kind) irelax,k,m real(r_kind) helmholtz_factor real(r_kind) fmax,fmean if(nxall(1) /= nx.or.nyall(1) /= ny) then write(6,*)' inconsistent input nx,ny in fmg, program stops' stop end if helmholtz_factor=zero if(helmholtz) helmholtz_factor=one !write(6,'(" at 1 in fmg, min,max f =",2e15.4)')minval(f),maxval(f) mg(1)%f=f mg(1)%v=zero do k=2,mgrid call f2c_full_weight(mg(k-1)%f,mg(k)%f,mg(k-1)%nx,mg(k-1)%ny,mg(k)%nx,mg(k)%ny) mg(k)%v=zero end do !write(6,'(" at 2 in fmg, min,max f =",2e15.4)')minval(f),maxval(f) ! main fmg loop: do m=mgrid,1,-1 ! exact solution at level mgrid k=mgrid do irelax=1,nrelax_solve if(helmholtz) then call relax(mg(k)%v,mg(k)%f, & mg(k)%dbarxh,mg(k)%ubarxh,mg(k)%obarxh, & mg(k)%dbaryh,mg(k)%ubaryh,mg(k)%obaryh, & mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,mg(k)%nx,mg(k)%ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) else call relax(mg(k)%v,mg(k)%f, & mg(k)%dbarx,mg(k)%ubarx,mg(k)%obarx, & mg(k)%dbary,mg(k)%ubary,mg(k)%obary, & mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,mg(k)%nx,mg(k)%ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) end if end do do k=mgrid-1,m,-1 ! update v(k) with v(k+1) call eval_bdys(mg(k+1)%v,mg(k+1),mg(k+1)%nx,mg(k+1)%ny) call c2f_bilinear(mg(k+1)%v,mg(k)%w,mg(k+1)%nx,mg(k+1)%ny,mg(k)%nx,mg(k)%ny) call eval_bdys(mg(k)%w,mg(k),mg(k)%nx,mg(k)%ny) mg(k)%v=mg(k)%v+mg(k)%w ! relax do irelax=1,nrelax2 if(helmholtz) then call relax(mg(k)%v,mg(k)%f, & mg(k)%dbarxh,mg(k)%ubarxh,mg(k)%obarxh, & mg(k)%dbaryh,mg(k)%ubaryh,mg(k)%obaryh, & mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,mg(k)%nx,mg(k)%ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) else call relax(mg(k)%v,mg(k)%f, & mg(k)%dbarx,mg(k)%ubarx,mg(k)%obarx, & mg(k)%dbary,mg(k)%ubary,mg(k)%obary, & mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,mg(k)%nx,mg(k)%ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) end if end do end do if(m > 1) then call eval_bdys(mg(m)%v,mg(m),mg(m)%nx,mg(m)%ny) call c2f_bicubic(mg(m)%v,mg(m-1)%v,mg(m)%nx,mg(m)%ny,mg(m-1)%nx,mg(m-1)%ny,mg(m-1)%bicubic_mask) call eval_bdys(mg(m-1)%v,mg(m-1),mg(m-1)%nx,mg(m-1)%ny) end if if(m == 1) exit do k=m-1,mgrid-1 ! relax do irelax=1,nrelax1 if(helmholtz) then call relax(mg(k)%v,mg(k)%f, & mg(k)%dbarxh,mg(k)%ubarxh,mg(k)%obarxh, & mg(k)%dbaryh,mg(k)%ubaryh,mg(k)%obaryh, & mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,mg(k)%nx,mg(k)%ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) else call relax(mg(k)%v,mg(k)%f, & mg(k)%dbarx,mg(k)%ubarx,mg(k)%obarx, & mg(k)%dbary,mg(k)%ubary,mg(k)%obary, & mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,mg(k)%nx,mg(k)%ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) end if end do ! compute residual call eval_bdys(mg(k)%v,mg(k),mg(k)%nx,mg(k)%ny) mg(k)%w=zero call forward_op(mg(k)%v,mg(k)%w,mg(k)%nx,mg(k)%ny,k,helmholtz_factor) mg(k)%w=mg(k)%f-mg(k)%w ! now have residual in mg(k)%w ! transfer residual to next coarser grid !write(6,'(" k, max val new residual=",i4,e14.5)')k, & ! maxval(abs(mg(k)%w(1:mg(k)%ny-1,1:mg(k)%nx-1))) !write(6,'(" k+1, max val prev residual=",i4,e14.5)')k+1, & ! maxval(abs(mg(k+1)%f(1:mg(k+1)%ny-1,1:mg(k+1)%nx-1))) call f2c_full_weight(mg(k)%w,mg(k+1)%f,mg(k)%nx,mg(k)%ny,mg(k+1)%nx,mg(k+1)%ny) !write(6,'(" k+1, max val transfered residual=",i4,e14.5)')k+1, & ! maxval(abs(mg(k+1)%f(1:mg(k+1)%ny-1,1:mg(k+1)%nx-1))) mg(k+1)%v=zero end do end do !write(6,'(" at 3 in fmg, min,max f =",2e15.4)')minval(f),maxval(f) v=mg(1)%v !write(6,'(" at 4 in fmg, min,max f =",2e15.4)')minval(f),maxval(f) end subroutine fmg subroutine fmg_ad(v,f,helmholtz,nrelax1,nrelax2,nrelax_solve,nx,ny) ! adjoint of full multigrid solution on fixed domain with homogeneous dirichlet boundary conditions ! solve delsqr(v) = f helmholtz=.false. ! (delsqr-R)(v) = f helmholtz=.true. use constants, only: zero,one integer(i_kind),intent(in)::nx,ny,nrelax1,nrelax2,nrelax_solve logical,intent(in):: helmholtz real(r_kind),intent(inout):: f(0:ny,0:nx) ! ! in actual practice, v is in and real(r_kind),intent(in):: v(0:ny,0:nx) ! ! f is inout, but make both inout for adjoint test integer(i_kind) irelax,k,m real(r_kind) helmholtz_factor real(r_kind) fmax,fmean if(nxall(1) /= nx.or.nyall(1) /= ny) then write(6,*)' inconsistent input nx,ny in fmg, program stops' stop end if helmholtz_factor=zero if(helmholtz) helmholtz_factor=one mg(1)%v=v mg(1)%f=f ! adjoint of main fmg loop: do m=1,mgrid if(m > 1) then do k=mgrid-1,m-1,-1 ! adjoint of transfer residual to next coarser grid call f2c_full_weight_ad(mg(k)%w,mg(k+1)%f,mg(k)%nx,mg(k)%ny,mg(k+1)%nx,mg(k+1)%ny) ! adjoint of compute residual mg(k)%f=mg(k)%w+mg(k)%f call forward_op_ad(mg(k)%w2,mg(k)%w,mg(k)%nx,mg(k)%ny,k,helmholtz_factor) mg(k)%v=mg(k)%v-mg(k)%w2 call eval_bdys_ad(mg(k)%v,mg(k),mg(k)%nx,mg(k)%ny) ! adjoint of relax do irelax=1,nrelax1 if(helmholtz) then call relax_ad(mg(k)%v,mg(k)%f, & mg(k)%dbarxh,mg(k)%ubarxh,mg(k)%obarxh, & mg(k)%dbaryh,mg(k)%ubaryh,mg(k)%obaryh, & mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,mg(k)%nx,mg(k)%ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) else call relax_ad(mg(k)%v,mg(k)%f, & mg(k)%dbarx,mg(k)%ubarx,mg(k)%obarx, & mg(k)%dbary,mg(k)%ubary,mg(k)%obary, & mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,mg(k)%nx,mg(k)%ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) end if end do end do call eval_bdys_ad(mg(m-1)%v,mg(m-1),mg(m-1)%nx,mg(m-1)%ny) call c2f_bicubic_ad(mg(m)%v,mg(m-1)%v,mg(m)%nx,mg(m)%ny,mg(m-1)%nx,mg(m-1)%ny,mg(m-1)%bicubic_mask) call eval_bdys_ad(mg(m)%v,mg(m),mg(m)%nx,mg(m)%ny) mg(m)%f=zero end if do k=m,mgrid-1 ! adjoint of relax do irelax=1,nrelax2 if(helmholtz) then call relax_ad(mg(k)%v,mg(k)%f, & mg(k)%dbarxh,mg(k)%ubarxh,mg(k)%obarxh, & mg(k)%dbaryh,mg(k)%ubaryh,mg(k)%obaryh, & mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,mg(k)%nx,mg(k)%ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) else call relax_ad(mg(k)%v,mg(k)%f, & mg(k)%dbarx,mg(k)%ubarx,mg(k)%obarx, & mg(k)%dbary,mg(k)%ubary,mg(k)%obary, & mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,mg(k)%nx,mg(k)%ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) end if end do ! adjoint of update v(k) with v(k+1) mg(k)%w=mg(k)%v call eval_bdys_ad(mg(k)%w,mg(k),mg(k)%nx,mg(k)%ny) call c2f_bilinear_ad(mg(k+1)%v,mg(k)%w,mg(k+1)%nx,mg(k+1)%ny,mg(k)%nx,mg(k)%ny) call eval_bdys_ad(mg(k+1)%v,mg(k+1),mg(k+1)%nx,mg(k+1)%ny) mg(k+1)%f=zero end do ! adjoint of exact solution at level mgrid k=mgrid do irelax=1,nrelax_solve if(helmholtz) then call relax_ad(mg(k)%v,mg(k)%f, & mg(k)%dbarxh,mg(k)%ubarxh,mg(k)%obarxh, & mg(k)%dbaryh,mg(k)%ubaryh,mg(k)%obaryh, & mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,mg(k)%nx,mg(k)%ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) else call relax_ad(mg(k)%v,mg(k)%f, & mg(k)%dbarx,mg(k)%ubarx,mg(k)%obarx, & mg(k)%dbary,mg(k)%ubary,mg(k)%obary, & mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,mg(k)%nx,mg(k)%ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) end if end do end do do k=mgrid,2,-1 call f2c_full_weight_ad(mg(k-1)%w,mg(k)%f,mg(k-1)%nx,mg(k-1)%ny,mg(k)%nx,mg(k)%ny) mg(k-1)%f=mg(k-1)%f+mg(k-1)%w end do f=mg(1)%f end subroutine fmg_ad subroutine vcycle(f,v,nx,ny,mgrid2,nrelax1,nrelax2,nrelax_solve,helmholtz) ! vcycle over mgrid2 grid levels, using nrelax1 relaxation cycles going from fine to coarse ! with residuals, nrelax_solve relaxation cycles on coarsest grid, and nrelax2 relaxation ! cycles going back up from coarse to fine with corrections. use constants, only: zero,one integer(i_kind),intent(in)::nx,ny,mgrid2,nrelax1,nrelax2,nrelax_solve logical,intent(in):: helmholtz real(r_kind),intent(in):: f(0:ny,0:nx) real(r_kind),intent(inout):: v(0:ny,0:nx) integer(i_kind) irelax,k real(r_kind) helmholtz_factor helmholtz_factor=zero if(helmholtz) helmholtz_factor=one mg(1)%f=f ! mg(1)%v=zero mg(1)%v=v ! fine to coarse residuals (descend left side of V cycle) do k=1,mgrid2-1 ! relax do irelax=1,nrelax1 if(helmholtz) then call relax(mg(k)%v,mg(k)%f, & mg(k)%dbarxh,mg(k)%ubarxh,mg(k)%obarxh, & mg(k)%dbaryh,mg(k)%ubaryh,mg(k)%obaryh, & mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,mg(k)%nx,mg(k)%ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) else call relax(mg(k)%v,mg(k)%f, & mg(k)%dbarx,mg(k)%ubarx,mg(k)%obarx, & mg(k)%dbary,mg(k)%ubary,mg(k)%obary, & mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,mg(k)%nx,mg(k)%ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) end if end do ! compute residual call eval_bdys(mg(k)%v,mg(k),mg(k)%nx,mg(k)%ny) mg(k)%w=zero call forward_op(mg(k)%v,mg(k)%w,mg(k)%nx,mg(k)%ny,k,helmholtz_factor) mg(k)%w=mg(k)%f-mg(k)%w ! now have residual in mg(k)%w ! transfer residual to next coarser grid call f2c_full_weight(mg(k)%w,mg(k+1)%f,mg(k)%nx,mg(k)%ny,mg(k+1)%nx,mg(k+1)%ny) mg(k+1)%v=zero end do ! solve (bottom of V cycle) k=mgrid2 do irelax=1,nrelax_solve if(helmholtz) then call relax(mg(k)%v,mg(k)%f, & mg(k)%dbarxh,mg(k)%ubarxh,mg(k)%obarxh, & mg(k)%dbaryh,mg(k)%ubaryh,mg(k)%obaryh, & mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,mg(k)%nx,mg(k)%ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) else call relax(mg(k)%v,mg(k)%f, & mg(k)%dbarx,mg(k)%ubarx,mg(k)%obarx, & mg(k)%dbary,mg(k)%ubary,mg(k)%obary, & mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,mg(k)%nx,mg(k)%ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) end if end do ! coarse to fine corrections (ascend right side of V cycle) do k=mgrid2-1,1,-1 ! update v(k) with v(k+1) call eval_bdys(mg(k+1)%v,mg(k+1),mg(k+1)%nx,mg(k+1)%ny) call c2f_bilinear(mg(k+1)%v,mg(k)%w,mg(k+1)%nx,mg(k+1)%ny,mg(k)%nx,mg(k)%ny) call eval_bdys(mg(k)%w,mg(k),mg(k)%nx,mg(k)%ny) mg(k)%v=mg(k)%v+mg(k)%w ! relax do irelax=1,nrelax2 if(helmholtz) then call relax(mg(k)%v,mg(k)%f, & mg(k)%dbarxh,mg(k)%ubarxh,mg(k)%obarxh, & mg(k)%dbaryh,mg(k)%ubaryh,mg(k)%obaryh, & mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,mg(k)%nx,mg(k)%ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) else call relax(mg(k)%v,mg(k)%f, & mg(k)%dbarx,mg(k)%ubarx,mg(k)%obarx, & mg(k)%dbary,mg(k)%ubary,mg(k)%obary, & mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,mg(k)%nx,mg(k)%ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) end if end do end do v=mg(1)%v end subroutine vcycle subroutine vcycle_ad(f,v,nx,ny,mgrid2,nrelax1,nrelax2,nrelax_solve,helmholtz) ! vcycle over mgrid2 grid levels, using nrelax1 relaxation cycles going from fine to coarse ! with residuals, nrelax_solve relaxation cycles on coarsest grid, and nrelax2 relaxation ! cycles going back up from coarse to fine with corrections. use constants, only: zero,one integer(i_kind),intent(in)::nx,ny,mgrid2,nrelax1,nrelax2,nrelax_solve logical,intent(in):: helmholtz real(r_kind),intent(inout):: f(0:ny,0:nx) real(r_kind),intent(inout):: v(0:ny,0:nx) integer(i_kind) irelax,k real(r_kind) helmholtz_factor helmholtz_factor=zero if(helmholtz) helmholtz_factor=one mg(1)%v=v mg(1)%f=f ! adjoint of coarse to fine corrections (ascend right side of V cycle) do k=1,mgrid2-1 ! adjoint of relax do irelax=1,nrelax2 if(helmholtz) then call relax_ad(mg(k)%v,mg(k)%f, & mg(k)%dbarxh,mg(k)%ubarxh,mg(k)%obarxh, & mg(k)%dbaryh,mg(k)%ubaryh,mg(k)%obaryh, & mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,mg(k)%nx,mg(k)%ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) else call relax_ad(mg(k)%v,mg(k)%f, & mg(k)%dbarx,mg(k)%ubarx,mg(k)%obarx, & mg(k)%dbary,mg(k)%ubary,mg(k)%obary, & mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,mg(k)%nx,mg(k)%ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) end if end do ! adjoint of update v(k) with v(k+1) mg(k)%w=mg(k)%v call eval_bdys_ad(mg(k)%w,mg(k),mg(k)%nx,mg(k)%ny) call c2f_bilinear_ad(mg(k+1)%v,mg(k)%w,mg(k+1)%nx,mg(k+1)%ny,mg(k)%nx,mg(k)%ny) call eval_bdys_ad(mg(k+1)%v,mg(k+1),mg(k+1)%nx,mg(k+1)%ny) mg(k+1)%f=zero end do ! adjoint of solve (bottom of V cycle) k=mgrid2 do irelax=1,nrelax_solve if(helmholtz) then call relax_ad(mg(k)%v,mg(k)%f, & mg(k)%dbarxh,mg(k)%ubarxh,mg(k)%obarxh, & mg(k)%dbaryh,mg(k)%ubaryh,mg(k)%obaryh, & mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,mg(k)%nx,mg(k)%ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) else call relax_ad(mg(k)%v,mg(k)%f, & mg(k)%dbarx,mg(k)%ubarx,mg(k)%obarx, & mg(k)%dbary,mg(k)%ubary,mg(k)%obary, & mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,mg(k)%nx,mg(k)%ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) end if end do ! adjoint of fine to coarse residuals (descend left side of V cycle) do k=mgrid2-1,1,-1 ! adjoint of transfer residual to next coarser grid call f2c_full_weight_ad(mg(k)%w,mg(k+1)%f,mg(k)%nx,mg(k)%ny,mg(k+1)%nx,mg(k+1)%ny) ! adjoint of compute residual call forward_op_ad(mg(k)%w2,mg(k)%w,mg(k)%nx,mg(k)%ny,k,helmholtz_factor) mg(k)%v=mg(k)%v-mg(k)%w2 mg(k)%f=mg(k)%f+mg(k)%w call eval_bdys_ad(mg(k)%v,mg(k),mg(k)%nx,mg(k)%ny) ! adjoint of relax do irelax=1,nrelax1 if(helmholtz) then call relax_ad(mg(k)%v,mg(k)%f, & mg(k)%dbarxh,mg(k)%ubarxh,mg(k)%obarxh, & mg(k)%dbaryh,mg(k)%ubaryh,mg(k)%obaryh, & mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,mg(k)%nx,mg(k)%ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) else call relax_ad(mg(k)%v,mg(k)%f, & mg(k)%dbarx,mg(k)%ubarx,mg(k)%obarx, & mg(k)%dbary,mg(k)%ubary,mg(k)%obary, & mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,mg(k)%nx,mg(k)%ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) end if end do end do v=mg(1)%v f=mg(1)%f end subroutine vcycle_ad subroutine test_hop(test_div,nx1,ny1,helmholtz_on) use constants, only: zero,one,two,rearth integer(i_kind),intent(in):: nx1,ny1 real(r_kind),intent(in):: test_div(0:ny1,0:nx1) logical,intent(in):: helmholtz_on integer(i_kind) irelax,k,i,j real(r_kind) helmholtz character(120) fname real(r_kind),allocatable:: f(:,:),v(:,:),w(:,:) real(r_kind) time0,timef ! run tests with test_div_full allocate(f(0:ny1,0:nx1)) allocate(v(0:ny1,0:nx1)) allocate(w(0:ny1,0:nx1)) f=test_div time0=timef() write(6,'(" min,max f before call fmg=",2e15.4)')minval(f),maxval(f) call fmg(v,f,helmholtz_on,1,1,20,mg(1)%nx,mg(1)%ny) write(6,'(" time in fmg =",f15.6," seconds")') .001*(timef()-time0) w=zero helmholtz=zero if(helmholtz_on) helmholtz=one call forward_op(v,w,mg(1)%nx,mg(1)%ny,1,helmholtz) w=f-w write(6,'(" min,max f after call fmg=",2e15.4)') & minval(f(1:ny1-1,1:nx1-1)), & maxval(f(1:ny1-1,1:nx1-1)) write(6,'(" min,max r after call fmg=",2e15.4)') & minval(w(1:ny1-1,1:nx1-1)), & maxval(w(1:ny1-1,1:nx1-1)) if(helmholtz_on) then call outgrad1(v,'vhs_test20x',mg(1)%nx+1,mg(1)%ny+1) call outgrad1(f,'fhs_test20x',mg(1)%nx+1,mg(1)%ny+1) call outgrad1(w,'rhs_test20x',mg(1)%nx+1,mg(1)%ny+1) else call outgrad1(v,'v_test20x',mg(1)%nx+1,mg(1)%ny+1) call outgrad1(f,'f_test20x',mg(1)%nx+1,mg(1)%ny+1) call outgrad1(w,'r_test20x',mg(1)%nx+1,mg(1)%ny+1) end if end subroutine test_hop subroutine test_relax use constants, only: zero,one integer(i_kind) i,j,k,lwave,nx,ny character(50) string !do k=1,mgrid do k=1,1 nx=mg(k)%nx ny=mg(k)%ny write(6,*)' at 1 in test_relax, k,nx,ny=',k,nx,ny call random_number(mg(k)%v) call eval_bdys(mg(k)%v,mg(k),mg(k)%nx,mg(k)%ny) mg(k)%z=mg(k)%v ! mg(k)%v=zero ! mg(k)%v=mg(k)%r00 ! mg(k)%z=mg(k)%r00 ! do j=-ny+1,ny-1 !! mg(k)%v=mg(k)%r00 !! mg(k)%z=mg(k)%r00 ! call forward_op_x(mg(k)%v,mg(k)%f,nx,ny,j,k,zero) ! end do ! mg(k)%w=-mg(k)%f ! do j=1,ny-1 ! mg(k)%w(j,1 )=mg(k)%w(j, 1)+mg(k)%h0m(j, 1)*mg(k)%v(j, 0) ! mg(k)%w(j,nx-1)=mg(k)%w(j, nx-1)+mg(k)%h0p(j, nx-1)*mg(k)%v(j, nx) ! call inverse_luop_x_1(mg(k)%w,mg(k)%v,mg(k)%dbarx,mg(k)%ubarx,mg(k)%obarx,nx,ny,j) ! end do ! call outgrad1(mg(k)%z,'exact_v_xtest',nx+1,ny+1) ! call outgrad1(mg(k)%v,'solve_v_xtest',nx+1,ny+1) ! if(1 /= 0) stop call eval_bdys(mg(k)%v,mg(k),mg(k)%nx,mg(k)%ny) mg(k)%f=zero call forward_op(mg(k)%v,mg(k)%f,nx,ny,k,zero) mg(k)%v(1:ny-1,1:nx-1)=zero !do i=1,1000 do i=1,60 call relax(mg(k)%v,mg(k)%f,mg(k)%dbarx,mg(k)%ubarx,mg(k)%obarx,mg(k)%dbary,mg(k)%ubary,mg(k)%obary,& mg(k)%h0p,mg(k)%h0m,mg(k)%hp0,mg(k)%hm0,nx,ny, & mg(k)%b_over_a_east_neg,mg(k)%b_over_a_north_neg) call eval_bdys(mg(k)%v,mg(k),mg(k)%nx,mg(k)%ny) mg(k)%w=zero call forward_op(mg(k)%v,mg(k)%w,nx,ny,k,zero) mg(k)%w=mg(k)%f-mg(k)%w write(6,'(" iter,max(|r|)/max(|f| = ",i5,e15.4)') i,maxval(abs(mg(k)%w(1:ny-1,1:nx-1)))/ & maxval(abs(mg(k)%f(1:ny-1,1:nx-1))) write(6,'(" iter,max(|v-v0|)/max(|v0| = ",i5,e15.4)') & i,maxval(abs(mg(k)%v-mg(k)%z))/maxval(abs(mg(k)%z)) if(i > 20) cycle write(string,'("f",i2.2)')i call outgrad1(mg(k)%f,trim(string),nx+1,ny+1) write(string,'("r_odev_xy",i2.2)')i call outgrad1(mg(k)%w,trim(string),nx+1,ny+1) end do call outgrad1(mg(k)%z,'exact_solution',nx+1,ny+1) call outgrad1(mg(k)%v,'relax60_solution',nx+1,ny+1) end do end subroutine test_relax subroutine test_adjoint use constants, only: zero,two real(r_kind) xtz,yty,helmholtz_factor,errmax integer(i_kind) i,ii,j,jj,kc,kf,ihem,nrelax1,nrelax2,nrelax_solve logical helmholtz integer(i_kind) mgrid2,ibicubic,ivar logical bicubic ! c2f_bilinear errmax=zero do kf=1,mgrid-1 kc=kf+1 call random_number(mg(kc)%v) mg(kc)%w=mg(kc)%v call c2f_bilinear(mg(kc)%v,mg(kf)%f,mg(kc)%nx,mg(kc)%ny,mg(kf)%nx,mg(kf)%ny) yty=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny yty=yty+mg(kf)%f(j,i)**2 end do end do call c2f_bilinear_ad(mg(kc)%v,mg(kf)%f,mg(kc)%nx,mg(kc)%ny,mg(kf)%nx,mg(kf)%ny) xtz=zero do i=0,mg(kc)%nx do j=0,mg(kc)%ny xtz=xtz+mg(kc)%v(j,i)*mg(kc)%w(j,i) end do end do write(6,'(" c2f_bilinear_ad (aT*a), kc,kf,xtz,yty=",2i3,2e22.15,e10.3)') & kc,kf,xtz,yty,abs(two*(xtz-yty)/(abs(xtz)+abs(yty))) if(abs(xtz)+abs(yty) /= zero) errmax=max(abs(two*(xtz-yty)/(abs(xtz)+abs(yty))),errmax) end do do kf=1,mgrid-1 kc=kf+1 call random_number(mg(kf)%v) mg(kf)%w=mg(kf)%v call c2f_bilinear_ad(mg(kc)%f,mg(kf)%v,mg(kc)%nx,mg(kc)%ny,mg(kf)%nx,mg(kf)%ny) yty=zero do i=0,mg(kc)%nx do j=0,mg(kc)%ny yty=yty+mg(kc)%f(j,i)**2 end do end do call c2f_bilinear(mg(kc)%f,mg(kf)%v,mg(kc)%nx,mg(kc)%ny,mg(kf)%nx,mg(kf)%ny) xtz=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny xtz=xtz+mg(kf)%v(j,i)*mg(kf)%w(j,i) end do end do write(6,'(" c2f_bilinear_ad (a*aT), kc,kf,xtz,yty=",2i3,2e22.15,e10.3)') & kc,kf,xtz,yty,abs(two*(xtz-yty)/(abs(xtz)+abs(yty))) if(abs(xtz)+abs(yty) /= zero) errmax=max(abs(two*(xtz-yty)/(abs(xtz)+abs(yty))),errmax) end do write(6,'(" max error for c2f_bilinear_ad=",e10.3)') errmax ! c2f_bicubic errmax=zero do kf=1,mgrid-1 kc=kf+1 call random_number(mg(kc)%v) mg(kc)%w=mg(kc)%v call c2f_bicubic(mg(kc)%v,mg(kf)%f,mg(kc)%nx,mg(kc)%ny,mg(kf)%nx,mg(kf)%ny,mg(kf)%bicubic_mask) yty=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny yty=yty+mg(kf)%f(j,i)**2 end do end do call c2f_bicubic_ad(mg(kc)%v,mg(kf)%f,mg(kc)%nx,mg(kc)%ny,mg(kf)%nx,mg(kf)%ny,mg(kf)%bicubic_mask) xtz=zero do i=0,mg(kc)%nx do j=0,mg(kc)%ny xtz=xtz+mg(kc)%v(j,i)*mg(kc)%w(j,i) end do end do write(6,'(" c2f_bicubic_ad (aT*a), kc,kf,xtz,yty=",2i3,2e22.15,e10.3)') & kc,kf,xtz,yty,abs(two*(xtz-yty)/(abs(xtz)+abs(yty))) if(abs(xtz)+abs(yty) /= zero) errmax=max(abs(two*(xtz-yty)/(abs(xtz)+abs(yty))),errmax) end do do kf=1,mgrid-1 kc=kf+1 call random_number(mg(kf)%v) mg(kf)%w=mg(kf)%v call c2f_bicubic_ad(mg(kc)%f,mg(kf)%v,mg(kc)%nx,mg(kc)%ny,mg(kf)%nx,mg(kf)%ny,mg(kf)%bicubic_mask) yty=zero do i=0,mg(kc)%nx do j=0,mg(kc)%ny yty=yty+mg(kc)%f(j,i)**2 end do end do call c2f_bicubic(mg(kc)%f,mg(kf)%v,mg(kc)%nx,mg(kc)%ny,mg(kf)%nx,mg(kf)%ny,mg(kf)%bicubic_mask) xtz=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny xtz=xtz+mg(kf)%v(j,i)*mg(kf)%w(j,i) end do end do write(6,'(" c2f_bicubic_ad (a*aT), kc,kf,xtz,yty=",2i3,2e22.15,e10.3)') & kc,kf,xtz,yty,abs(two*(xtz-yty)/(abs(xtz)+abs(yty))) if(abs(xtz)+abs(yty) /= zero) errmax=max(abs(two*(xtz-yty)/(abs(xtz)+abs(yty))),errmax) end do write(6,'(" max error for c2f_bicubic_ad=",e10.3)') errmax ! eval_bdys errmax=zero do kf=1,mgrid call random_number(mg(kf)%v) mg(kf)%w=mg(kf)%v call eval_bdys(mg(kf)%v,mg(kf),mg(kf)%nx,mg(kf)%ny) yty=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny yty=yty+mg(kf)%v(j,i)**2 end do end do call eval_bdys_ad(mg(kf)%v,mg(kf),mg(kf)%nx,mg(kf)%ny) xtz=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny xtz=xtz+mg(kf)%v(j,i)*mg(kf)%w(j,i) end do end do write(6,'(" eval_bdys_ad (aT*a), k,xtz,yty=",i3,2e22.15,e10.3)') & kf,xtz,yty,abs(two*(xtz-yty)/(abs(xtz)+abs(yty))) if(abs(xtz)+abs(yty) /= zero) errmax=max(abs(two*(xtz-yty)/(abs(xtz)+abs(yty))),errmax) end do do kf=1,mgrid call random_number(mg(kf)%v) mg(kf)%w=mg(kf)%v call eval_bdys_ad(mg(kf)%v,mg(kf),mg(kf)%nx,mg(kf)%ny) yty=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny yty=yty+mg(kf)%v(j,i)**2 end do end do call eval_bdys(mg(kf)%v,mg(kf),mg(kf)%nx,mg(kf)%ny) xtz=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny xtz=xtz+mg(kf)%v(j,i)*mg(kf)%w(j,i) end do end do write(6,'(" eval_bdys_ad (a*aT), k,xtz,yty=",i3,2e22.15,e10.3)') & kf,xtz,yty,abs(two*(xtz-yty)/(abs(xtz)+abs(yty))) if(abs(xtz)+abs(yty) /= zero) errmax=max(abs(two*(xtz-yty)/(abs(xtz)+abs(yty))),errmax) end do write(6,'(" max error for eval_bdys_ad=",e10.3)') errmax ! f2c_full_weight errmax=zero do kf=1,mgrid-1 kc=kf+1 call random_number(mg(kf)%v) mg(kf)%w=mg(kf)%v call f2c_full_weight(mg(kf)%v,mg(kc)%f,mg(kf)%nx,mg(kf)%ny,mg(kc)%nx,mg(kc)%ny) yty=zero do i=0,mg(kc)%nx do j=0,mg(kc)%ny yty=yty+mg(kc)%f(j,i)**2 end do end do call f2c_full_weight_ad(mg(kf)%v,mg(kc)%f,mg(kf)%nx,mg(kf)%ny,mg(kc)%nx,mg(kc)%ny) xtz=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny xtz=xtz+mg(kf)%v(j,i)*mg(kf)%w(j,i) end do end do write(6,'(" f2c_full_weight_ad (aT*a), kc,kf,xtz,yty=",2i3,2e22.15,e10.3)') & kc,kf,xtz,yty,abs(two*(xtz-yty)/(abs(xtz)+abs(yty))) if(abs(xtz)+abs(yty) /= zero) errmax=max(abs(two*(xtz-yty)/(abs(xtz)+abs(yty))),errmax) end do do kf=1,mgrid-1 kc=kf+1 call random_number(mg(kc)%v) mg(kc)%w=mg(kc)%v call f2c_full_weight_ad(mg(kf)%f,mg(kc)%v,mg(kf)%nx,mg(kf)%ny,mg(kc)%nx,mg(kc)%ny) yty=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny yty=yty+mg(kf)%f(j,i)**2 end do end do call f2c_full_weight(mg(kf)%f,mg(kc)%v,mg(kf)%nx,mg(kf)%ny,mg(kc)%nx,mg(kc)%ny) xtz=zero do i=0,mg(kc)%nx do j=0,mg(kc)%ny xtz=xtz+mg(kc)%v(j,i)*mg(kc)%w(j,i) end do end do write(6,'(" f2c_full_weight_ad (a*aT), kc,kf,xtz,yty=",2i3,2e22.15,e10.3)') & kc,kf,xtz,yty,abs(two*(xtz-yty)/(abs(xtz)+abs(yty))) if(abs(xtz)+abs(yty) /= zero) errmax=max(abs(two*(xtz-yty)/(abs(xtz)+abs(yty))),errmax) end do write(6,'(" max error for f2c_full_weight_ad=",e10.3)') errmax ! forward_op errmax=zero do ihem=0,1 helmholtz_factor=ihem do kf=1,mgrid call random_number(mg(kf)%v) mg(kf)%w=mg(kf)%v call forward_op(mg(kf)%v,mg(kf)%f,mg(kf)%nx,mg(kf)%ny,kf,helmholtz_factor) yty=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny yty=yty+mg(kf)%f(j,i)**2 end do end do call forward_op_ad(mg(kf)%v,mg(kf)%f,mg(kf)%nx,mg(kf)%ny,kf,helmholtz_factor) xtz=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny xtz=xtz+mg(kf)%v(j,i)*mg(kf)%w(j,i) end do end do write(6,'(" forward_op_ad (aT*a), h,k,xtz,yty=",i2,i3,2e22.15,e10.3)') & ihem,kf,xtz,yty,abs(two*(xtz-yty)/(abs(xtz)+abs(yty))) if(abs(xtz)+abs(yty) /= zero) errmax=max(abs(two*(xtz-yty)/(abs(xtz)+abs(yty))),errmax) end do do kf=1,mgrid call random_number(mg(kf)%f) mg(kf)%w=mg(kf)%f call forward_op_ad(mg(kf)%v,mg(kf)%f,mg(kf)%nx,mg(kf)%ny,kf,helmholtz_factor) yty=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny yty=yty+mg(kf)%v(j,i)**2 end do end do call forward_op(mg(kf)%v,mg(kf)%f,mg(kf)%nx,mg(kf)%ny,kf,helmholtz_factor) xtz=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny xtz=xtz+mg(kf)%f(j,i)*mg(kf)%w(j,i) end do end do write(6,'(" forward_op_ad (a*aT), h,k,xtz,yty=",i2,i3,2e22.15,e10.3)') & ihem,kf,xtz,yty,abs(two*(xtz-yty)/(abs(xtz)+abs(yty))) if(abs(xtz)+abs(yty) /= zero) errmax=max(abs(two*(xtz-yty)/(abs(xtz)+abs(yty))),errmax) end do end do write(6,'(" max error for forward_op_ad=",e10.3)') errmax ! inverse_luop_x errmax=zero do ii=1,1 do jj=1,2 do ihem=0,1 helmholtz_factor=ihem do kf=1,mgrid call random_number(mg(kf)%v) mg(kf)%w=mg(kf)%v mg(kf)%f=zero if(ihem == 0) then call inverse_luop_x(mg(kf)%v,mg(kf)%f,mg(kf)%dbarx,mg(kf)%ubarx,mg(kf)%obarx, & mg(kf)%nx,mg(kf)%ny,jj,ii) else call inverse_luop_x(mg(kf)%v,mg(kf)%f,mg(kf)%dbarxh,mg(kf)%ubarxh,mg(kf)%obarxh, & mg(kf)%nx,mg(kf)%ny,jj,ii) end if yty=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny yty=yty+mg(kf)%f(j,i)**2 end do end do mg(kf)%v=zero if(ihem == 0) then call inverse_luop_x_ad(mg(kf)%v,mg(kf)%f,mg(kf)%dbarx,mg(kf)%ubarx,mg(kf)%obarx, & mg(kf)%nx,mg(kf)%ny,jj,ii) else call inverse_luop_x_ad(mg(kf)%v,mg(kf)%f,mg(kf)%dbarxh,mg(kf)%ubarxh,mg(kf)%obarxh, & mg(kf)%nx,mg(kf)%ny,jj,ii) end if xtz=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny xtz=xtz+mg(kf)%v(j,i)*mg(kf)%w(j,i) end do end do write(6,'(" inverse_luop_x_ad (aT*a), jj,ii,h,k,xtz,yty=",3i2,i3,2e22.15,e10.3)') & jj,ii,ihem,kf,xtz,yty,abs(two*(xtz-yty)/(abs(xtz)+abs(yty))) if(abs(xtz)+abs(yty) /= zero) errmax=max(abs(two*(xtz-yty)/(abs(xtz)+abs(yty))),errmax) end do do kf=1,mgrid call random_number(mg(kf)%f) mg(kf)%w=mg(kf)%f mg(kf)%v=zero if(ihem == 0) then call inverse_luop_x_ad(mg(kf)%v,mg(kf)%f,mg(kf)%dbarx,mg(kf)%ubarx,mg(kf)%obarx, & mg(kf)%nx,mg(kf)%ny,jj,ii) else call inverse_luop_x_ad(mg(kf)%v,mg(kf)%f,mg(kf)%dbarxh,mg(kf)%ubarxh,mg(kf)%obarxh, & mg(kf)%nx,mg(kf)%ny,jj,ii) end if yty=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny yty=yty+mg(kf)%v(j,i)**2 end do end do mg(kf)%f=zero if(ihem == 0) then call inverse_luop_x(mg(kf)%v,mg(kf)%f,mg(kf)%dbarx,mg(kf)%ubarx,mg(kf)%obarx, & mg(kf)%nx,mg(kf)%ny,jj,ii) else call inverse_luop_x(mg(kf)%v,mg(kf)%f,mg(kf)%dbarxh,mg(kf)%ubarxh,mg(kf)%obarxh, & mg(kf)%nx,mg(kf)%ny,jj,ii) end if xtz=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny xtz=xtz+mg(kf)%f(j,i)*mg(kf)%w(j,i) end do end do write(6,'(" inverse_luop_x_ad (a*aT), jj,ii,h,k,xtz,yty=",3i2,i3,2e22.15,e10.3)') & jj,ii,ihem,kf,xtz,yty,abs(two*(xtz-yty)/(abs(xtz)+abs(yty))) if(abs(xtz)+abs(yty) /= zero) errmax=max(abs(two*(xtz-yty)/(abs(xtz)+abs(yty))),errmax) end do end do end do end do write(6,'(" max error for inverse_luop_x_ad=",e10.3)') errmax ! inverse_luop_y errmax=zero do ii=1,2 do jj=1,1 do ihem=0,1 helmholtz_factor=ihem do kf=1,mgrid call random_number(mg(kf)%v) mg(kf)%w=mg(kf)%v mg(kf)%f=zero if(ihem == 0) then call inverse_luop_y(mg(kf)%v,mg(kf)%f,mg(kf)%dbary,mg(kf)%ubary,mg(kf)%obary, & mg(kf)%nx,mg(kf)%ny,ii,jj) else call inverse_luop_y(mg(kf)%v,mg(kf)%f,mg(kf)%dbaryh,mg(kf)%ubaryh,mg(kf)%obaryh, & mg(kf)%nx,mg(kf)%ny,ii,jj) end if yty=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny yty=yty+mg(kf)%f(j,i)**2 end do end do mg(kf)%v=zero if(ihem == 0) then call inverse_luop_y_ad(mg(kf)%v,mg(kf)%f,mg(kf)%dbary,mg(kf)%ubary,mg(kf)%obary, & mg(kf)%nx,mg(kf)%ny,ii,jj) else call inverse_luop_y_ad(mg(kf)%v,mg(kf)%f,mg(kf)%dbaryh,mg(kf)%ubaryh,mg(kf)%obaryh, & mg(kf)%nx,mg(kf)%ny,ii,jj) end if xtz=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny xtz=xtz+mg(kf)%v(j,i)*mg(kf)%w(j,i) end do end do write(6,'(" inverse_luop_y_ad (aT*a), ii,jj,h,k,xtz,yty=",3i2,i3,2e22.15,e10.3)') & ii,jj,ihem,kf,xtz,yty,abs(two*(xtz-yty)/(abs(xtz)+abs(yty))) if(abs(xtz)+abs(yty) /= zero) errmax=max(abs(two*(xtz-yty)/(abs(xtz)+abs(yty))),errmax) end do do kf=1,mgrid call random_number(mg(kf)%f) mg(kf)%w=mg(kf)%f mg(kf)%v=zero if(ihem == 0) then call inverse_luop_y_ad(mg(kf)%v,mg(kf)%f,mg(kf)%dbary,mg(kf)%ubary,mg(kf)%obary, & mg(kf)%nx,mg(kf)%ny,ii,jj) else call inverse_luop_y_ad(mg(kf)%v,mg(kf)%f,mg(kf)%dbaryh,mg(kf)%ubaryh,mg(kf)%obaryh, & mg(kf)%nx,mg(kf)%ny,ii,jj) end if yty=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny yty=yty+mg(kf)%v(j,i)**2 end do end do mg(kf)%f=zero if(ihem == 0) then call inverse_luop_y(mg(kf)%v,mg(kf)%f,mg(kf)%dbary,mg(kf)%ubary,mg(kf)%obary, & mg(kf)%nx,mg(kf)%ny,ii,jj) else call inverse_luop_y(mg(kf)%v,mg(kf)%f,mg(kf)%dbaryh,mg(kf)%ubaryh,mg(kf)%obaryh, & mg(kf)%nx,mg(kf)%ny,ii,jj) end if xtz=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny xtz=xtz+mg(kf)%f(j,i)*mg(kf)%w(j,i) end do end do write(6,'(" inverse_luop_y_ad (a*aT), ii,jj,h,k,xtz,yty=",3i2,i3,2e22.15,e10.3)') & ii,jj,ihem,kf,xtz,yty,abs(two*(xtz-yty)/(abs(xtz)+abs(yty))) if(abs(xtz)+abs(yty) /= zero) errmax=max(abs(two*(xtz-yty)/(abs(xtz)+abs(yty))),errmax) end do end do end do end do write(6,'(" max error for inverse_luop_y_ad=",e10.3)') errmax ! relax_ad errmax=zero do ihem=0,1 do kf=1,mgrid call random_number(mg(kf)%v) call random_number(mg(kf)%f) mg(kf)%w=mg(kf)%v mg(kf)%y=mg(kf)%f if(ihem == 1) then call relax(mg(kf)%v,mg(kf)%f, & mg(kf)%dbarxh,mg(kf)%ubarxh,mg(kf)%obarxh, & mg(kf)%dbaryh,mg(kf)%ubaryh,mg(kf)%obaryh, & mg(kf)%h0p,mg(kf)%h0m,mg(kf)%hp0,mg(kf)%hm0, & mg(kf)%nx,mg(kf)%ny,mg(kf)%b_over_a_east_neg,mg(kf)%b_over_a_north_neg) else call relax(mg(kf)%v,mg(kf)%f, & mg(kf)%dbarx,mg(kf)%ubarx,mg(kf)%obarx, & mg(kf)%dbary,mg(kf)%ubary,mg(kf)%obary, & mg(kf)%h0p,mg(kf)%h0m,mg(kf)%hp0,mg(kf)%hm0, & mg(kf)%nx,mg(kf)%ny,mg(kf)%b_over_a_east_neg,mg(kf)%b_over_a_north_neg) end if yty=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny yty=yty+mg(kf)%f(j,i)**2+mg(kf)%v(j,i)**2 end do end do if(ihem == 1) then call relax_ad(mg(kf)%v,mg(kf)%f, & mg(kf)%dbarxh,mg(kf)%ubarxh,mg(kf)%obarxh, & mg(kf)%dbaryh,mg(kf)%ubaryh,mg(kf)%obaryh, & mg(kf)%h0p,mg(kf)%h0m,mg(kf)%hp0,mg(kf)%hm0, & mg(kf)%nx,mg(kf)%ny,mg(kf)%b_over_a_east_neg,mg(kf)%b_over_a_north_neg) else call relax_ad(mg(kf)%v,mg(kf)%f, & mg(kf)%dbarx,mg(kf)%ubarx,mg(kf)%obarx, & mg(kf)%dbary,mg(kf)%ubary,mg(kf)%obary, & mg(kf)%h0p,mg(kf)%h0m,mg(kf)%hp0,mg(kf)%hm0, & mg(kf)%nx,mg(kf)%ny,mg(kf)%b_over_a_east_neg,mg(kf)%b_over_a_north_neg) end if xtz=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny xtz=xtz+mg(kf)%v(j,i)*mg(kf)%w(j,i)+mg(kf)%f(j,i)*mg(kf)%y(j,i) end do end do write(6,'(" relax_ad (aT*a), h,k,xtz,yty=",i2,i3,2e22.15,e10.3)') & ihem,kf,xtz,yty,abs(two*(xtz-yty)/(abs(xtz)+abs(yty))) if(abs(xtz)+abs(yty) /= zero) errmax=max(abs(two*(xtz-yty)/(abs(xtz)+abs(yty))),errmax) end do do kf=1,mgrid call random_number(mg(kf)%v) call random_number(mg(kf)%f) mg(kf)%w=mg(kf)%v mg(kf)%y=mg(kf)%f if(ihem == 1) then call relax_ad(mg(kf)%v,mg(kf)%f, & mg(kf)%dbarxh,mg(kf)%ubarxh,mg(kf)%obarxh, & mg(kf)%dbaryh,mg(kf)%ubaryh,mg(kf)%obaryh, & mg(kf)%h0p,mg(kf)%h0m,mg(kf)%hp0,mg(kf)%hm0, & mg(kf)%nx,mg(kf)%ny,mg(kf)%b_over_a_east_neg,mg(kf)%b_over_a_north_neg) else call relax_ad(mg(kf)%v,mg(kf)%f, & mg(kf)%dbarx,mg(kf)%ubarx,mg(kf)%obarx, & mg(kf)%dbary,mg(kf)%ubary,mg(kf)%obary, & mg(kf)%h0p,mg(kf)%h0m,mg(kf)%hp0,mg(kf)%hm0, & mg(kf)%nx,mg(kf)%ny,mg(kf)%b_over_a_east_neg,mg(kf)%b_over_a_north_neg) end if yty=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny yty=yty+mg(kf)%f(j,i)**2+mg(kf)%v(j,i)**2 end do end do if(ihem == 1) then call relax(mg(kf)%v,mg(kf)%f, & mg(kf)%dbarxh,mg(kf)%ubarxh,mg(kf)%obarxh, & mg(kf)%dbaryh,mg(kf)%ubaryh,mg(kf)%obaryh, & mg(kf)%h0p,mg(kf)%h0m,mg(kf)%hp0,mg(kf)%hm0, & mg(kf)%nx,mg(kf)%ny,mg(kf)%b_over_a_east_neg,mg(kf)%b_over_a_north_neg) else call relax(mg(kf)%v,mg(kf)%f, & mg(kf)%dbarx,mg(kf)%ubarx,mg(kf)%obarx, & mg(kf)%dbary,mg(kf)%ubary,mg(kf)%obary, & mg(kf)%h0p,mg(kf)%h0m,mg(kf)%hp0,mg(kf)%hm0, & mg(kf)%nx,mg(kf)%ny,mg(kf)%b_over_a_east_neg,mg(kf)%b_over_a_north_neg) end if xtz=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny xtz=xtz+mg(kf)%v(j,i)*mg(kf)%w(j,i)+mg(kf)%f(j,i)*mg(kf)%y(j,i) end do end do write(6,'(" relax_ad (a*aT), h,k,xtz,yty=",i2,i3,2e22.15,e10.3)') & ihem,kf,xtz,yty,abs(two*(xtz-yty)/(abs(xtz)+abs(yty))) if(abs(xtz)+abs(yty) /= zero) errmax=max(abs(two*(xtz-yty)/(abs(xtz)+abs(yty))),errmax) end do end do write(6,'(" max error for relax_ad=",e10.3)') errmax ! vcycle_ad errmax=zero do ivar=1,2 do ihem=0,1 helmholtz=ihem == 1 do nrelax1=1,1 do nrelax2=1,1 do nrelax_solve=1,1 do mgrid2=2,mgrid ! do mgrid2=2,2 kf=1 mg(kf)%y=zero ; mg(kf)%z=zero if(ivar == 1) call random_number(mg(kf)%y) if(ivar == 2) call random_number(mg(kf)%z) mg(kf)%yy=mg(kf)%y mg(kf)%zz=mg(kf)%z call vcycle(mg(kf)%y,mg(kf)%z,mg(kf)%nx,mg(kf)%ny,mgrid2,nrelax1,nrelax2,nrelax_solve,helmholtz) yty=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny yty=yty+mg(kf)%y(j,i)**2+mg(kf)%z(j,i)**2 end do end do call vcycle_ad(mg(kf)%y,mg(kf)%z,mg(kf)%nx,mg(kf)%ny,mgrid2,nrelax1,nrelax2,nrelax_solve,helmholtz) xtz=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny xtz=xtz+mg(kf)%y(j,i)*mg(kf)%yy(j,i)+mg(kf)%z(j,i)*mg(kf)%zz(j,i) end do end do write(6,'(" vcycle_ad (aT*a), vhn12sm2,xz,yy=",6i2,2e22.15,e10.3)') & ivar,ihem,nrelax1,nrelax2,nrelax_solve,mgrid2,xtz,yty,abs(two*(xtz-yty)/(abs(xtz)+abs(yty))) if(abs(xtz)+abs(yty) /= zero) errmax=max(abs(two*(xtz-yty)/(abs(xtz)+abs(yty))),errmax) mg(kf)%y=zero ; mg(kf)%z=zero if(ivar == 1) call random_number(mg(kf)%y) if(ivar == 2) call random_number(mg(kf)%z) mg(kf)%yy=mg(kf)%y mg(kf)%zz=mg(kf)%z call vcycle_ad(mg(kf)%y,mg(kf)%z,mg(kf)%nx,mg(kf)%ny,mgrid2,nrelax1,nrelax2,nrelax_solve,helmholtz) yty=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny yty=yty+mg(kf)%y(j,i)**2+mg(kf)%z(j,i)**2 end do end do call vcycle(mg(kf)%y,mg(kf)%z,mg(kf)%nx,mg(kf)%ny,mgrid2,nrelax1,nrelax2,nrelax_solve,helmholtz) xtz=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny xtz=xtz+mg(kf)%y(j,i)*mg(kf)%yy(j,i)+mg(kf)%z(j,i)*mg(kf)%zz(j,i) end do end do write(6,'(" vcycle_ad (a*aT), vhn12sm2,xz,yy=",6i2,2e22.15,e10.3)') & ivar,ihem,nrelax1,nrelax2,nrelax_solve,mgrid2,xtz,yty,abs(two*(xtz-yty)/(abs(xtz)+abs(yty))) if(abs(xtz)+abs(yty) /= zero) errmax=max(abs(two*(xtz-yty)/(abs(xtz)+abs(yty))),errmax) end do end do end do end do end do end do write(6,'(" max error for vcycle_ad=",e10.3)') errmax ! fmg_ad errmax=zero do ivar=1,2 do ihem=0,1 helmholtz=ihem == 1 do nrelax1=1,1 do nrelax2=1,1 do nrelax_solve=25,25 kf=1 mg(kf)%y=zero ; mg(kf)%z=zero if(ivar == 1) call random_number(mg(kf)%y) if(ivar == 2) call random_number(mg(kf)%z) mg(kf)%yy=mg(kf)%y mg(kf)%zz=mg(kf)%z call fmg(mg(kf)%y,mg(kf)%z,helmholtz,nrelax1,nrelax2,nrelax_solve,mg(kf)%nx,mg(kf)%ny) yty=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny yty=yty+mg(kf)%y(j,i)**2+mg(kf)%z(j,i)**2 end do end do call fmg_ad(mg(kf)%y,mg(kf)%z,helmholtz,nrelax1,nrelax2,nrelax_solve,mg(kf)%nx,mg(kf)%ny) xtz=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny xtz=xtz+mg(kf)%y(j,i)*mg(kf)%yy(j,i)+mg(kf)%z(j,i)*mg(kf)%zz(j,i) end do end do write(6,'(" fmg_ad (aT*a), vhn12s,xz,yy=",5i2,2e22.15,e10.3)') & ivar,ihem,nrelax1,nrelax2,nrelax_solve,xtz,yty,abs(two*(xtz-yty)/(abs(xtz)+abs(yty))) if(abs(xtz)+abs(yty) /= zero) errmax=max(abs(two*(xtz-yty)/(abs(xtz)+abs(yty))),errmax) mg(kf)%y=zero ; mg(kf)%z=zero if(ivar == 1) call random_number(mg(kf)%y) if(ivar == 2) call random_number(mg(kf)%z) mg(kf)%yy=mg(kf)%y mg(kf)%zz=mg(kf)%z call fmg_ad(mg(kf)%y,mg(kf)%z,helmholtz,nrelax1,nrelax2,nrelax_solve,mg(kf)%nx,mg(kf)%ny) yty=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny yty=yty+mg(kf)%y(j,i)**2+mg(kf)%z(j,i)**2 end do end do call fmg(mg(kf)%y,mg(kf)%z,helmholtz,nrelax1,nrelax2,nrelax_solve,mg(kf)%nx,mg(kf)%ny) xtz=zero do i=0,mg(kf)%nx do j=0,mg(kf)%ny xtz=xtz+mg(kf)%y(j,i)*mg(kf)%yy(j,i)+mg(kf)%z(j,i)*mg(kf)%zz(j,i) end do end do write(6,'(" fmg_ad (a*aT), vhn12s,xz,yy=",5i2,2e22.15,e10.3)') & ivar,ihem,nrelax1,nrelax2,nrelax_solve,xtz,yty,abs(two*(xtz-yty)/(abs(xtz)+abs(yty))) if(abs(xtz)+abs(yty) /= zero) errmax=max(abs(two*(xtz-yty)/(abs(xtz)+abs(yty))),errmax) end do end do end do end do end do write(6,'(" max error for fmg_ad=",e10.3)') errmax end subroutine test_adjoint end module helmholtz_fmg_mod subroutine fmg_initialize(mype) use kinds, only: r_kind,i_kind,r_single use constants, only: deg2rad,rad2deg,zero,half,one,rearth use gridmod, only: region_lat,region_dx,region_dy,nlon,nlat use helmholtz_fmg_mod, only: generate_grids,init_mgvars use mod_vtrans, only: nvmodes_keep,depths implicit none integer(i_kind),intent(in):: mype integer(i_kind) mode_number real(r_kind) phi0 ! for this version, make parallel by vertical modes, so only have nvmodes_keep copies ! this is most simpleminded approach--later modify to take advantage of all processors ! available. if(mype+1 > 2*nvmodes_keep) return mode_number=mod(mype,nvmodes_keep)+1 phi0=depths(mode_number) write(6,*)' in fmg_initialize, mype,mode_number,phi0=',mype,mode_number,phi0 call generate_grids(nlon,nlat) call init_mgvars(region_dx,region_dy,region_lat,phi0) end subroutine fmg_initialize subroutine fmg_initialize_e(mype) use kinds, only: r_kind,i_kind,r_single use constants, only: deg2rad,rad2deg,zero,half,one,rearth use gridmod, only: region_lat,region_dx,region_dy,nlon,nlat,txy2ll use helmholtz_fmg_mod, only: generate_grids,init_mgvars use mod_vtrans, only: nvmodes_keep,depths use zrnmi_mod, only: zrnmi_initialize implicit none integer(i_kind),intent(in):: mype integer(i_kind) mode_number real(r_kind) phi0 real(r_kind) region_lat_e(0:nlat+1,0:nlon+1) real(r_kind) region_dx_e(0:nlat+1,0:nlon+1),region_dy_e(0:nlat+1,0:nlon+1) integer(i_kind) i,j real(r_kind) rx,ry,rlon,rlat ! for this version, make parallel by vertical modes, so only have nvmodes_keep copies ! this is most simpleminded approach--later modify to take advantage of all processors ! available. call zrnmi_initialize(mype) if(mype+1 > 2*nvmodes_keep) return mode_number=mod(mype,nvmodes_keep)+1 phi0=depths(mode_number) write(6,*)' in fmg_initialize_e, mype,mode_number,phi0=',mype,mode_number,phi0 do j=1,nlon do i=1,nlat region_lat_e(i,j)=region_lat(i,j) region_dx_e(i,j)=region_dx(i,j) region_dy_e(i,j)=region_dy(i,j) end do end do do j=1,nlon rx=j ry=0 call txy2ll(rx,ry,rlon,region_lat_e(0,j)) region_dx_e(0,j)=region_dx_e(1,j) region_dy_e(0,j)=region_dy_e(1,j) ry=nlat+1 call txy2ll(rx,ry,rlon,region_lat_e(nlat+1,j)) region_dx_e(nlat+1,j)=region_dx_e(nlat,j) region_dy_e(nlat+1,j)=region_dy_e(nlat,j) end do do i=0,nlat+1 ry=i rx=0 call txy2ll(rx,ry,rlon,region_lat_e(i,0)) region_dx_e(i,0)=region_dx_e(i,1) region_dy_e(i,0)=region_dy_e(i,1) rx=nlon+1 call txy2ll(rx,ry,rlon,region_lat_e(i,nlon+1)) region_dx_e(i,nlon+1)=region_dx_e(i,nlon) region_dy_e(i,nlon+1)=region_dy_e(i,nlon) end do if(mype == 0) call outgrad1(region_lat_e,'region_lat_e',nlon+2,nlat+2) if(mype == 0) call outgrad1(region_dx_e,'region_dx_e',nlon+2,nlat+2) if(mype == 0) call outgrad1(region_dy_e,'region_dy_e',nlon+2,nlat+2) call generate_grids(nlon+2,nlat+2) call init_mgvars(region_dx_e,region_dy_e,region_lat_e,phi0) end subroutine fmg_initialize_e subroutine fmg_strong_bal_correction(u_t,v_t,t_t,ps_t,psi,chi,t,ps,bal_diagnostic,fullfield,update,mype) use kinds, only: r_kind,i_kind use constants, only: zero,two,omega use gridmod, only: lat2,lon2,nsig,nlat,nlon,region_lat use gridmod, only: istart,jstart use mod_vtrans, only: vtrans,vtrans_inv,nvmodes_keep,depths use mpimod, only: mpi_comm_world,ierror,mpi_sum,mpi_rtype use zrnmi_mod, only: zrnmi_filter_uvm2 use jfunc,only: jiter use hybrid_ensemble_parameters, only: uv_hyb_ens implicit none integer(i_kind),intent(in)::mype logical,intent(in)::bal_diagnostic,update,fullfield real(r_kind),dimension(lat2,lon2,nsig),intent(in)::u_t,v_t,t_t real(r_kind),dimension(lat2,lon2),intent(in)::ps_t real(r_kind),dimension(lat2,lon2,nsig),intent(inout)::psi,chi,t real(r_kind),dimension(lat2,lon2),intent(inout)::ps real(r_kind),dimension(lat2,lon2,nvmodes_keep)::utilde,vtilde,mtilde real(r_kind),dimension(lat2,lon2,nvmodes_keep)::delpsitilde,delchitilde,delmtilde,dummytilde real(r_kind),dimension(nlat,nlon)::u0t,v0t,m0t,div0t,vor0t,rhs,mtg,delm,divtg,vortg,mtg_x,mtg_y real(r_kind),dimension(nlat,nlon)::delvor,deldiv,delpsi,delchi,u_psi,v_psi,u_chi,v_chi,delf1,delf2 real(r_kind),dimension(lat2,lon2,nsig)::dpsi,dchi,dt real(r_kind),dimension(lat2,lon2)::dps integer(i_kind) mode_number_a,mode_number_b real(r_kind) bal_a (nvmodes_keep),bal_b (nvmodes_keep) real(r_kind) bal_a0(nvmodes_keep),bal_b0(nvmodes_keep) real(r_kind) bal(nvmodes_keep) real(r_kind),dimension(lat2,lon2,nvmodes_keep)::delupsitilde,delvpsitilde,deluchitilde,delvchitilde real(r_kind),dimension(lat2,lon2,nvmodes_keep)::delutilde,delvtilde character(1) bourke_mcgregor_scheme !??????????????????????????????????? real(r_kind),dimension(lat2,lon2,nvmodes_keep)::utild2,vtild2,mtild2 real(r_kind) errmaxu,errmaxv,errmaxm real(r_kind),dimension(nlat,nlon)::uwork,vwork,twork,pwork integer(i_kind) iy,jx,i,j,k character(100) fname real(r_kind) errmaxpsi,errmaxchi,errmaxdummy !??????????????????????????????????? if(uv_hyb_ens) then write(0,*)' STOP IN fmg_strong_bal_correction--NOT ABLE TO ACCEPT uv_hyb_ens=.true. YET' call stop2(998) end if !bourke_mcgregor_scheme='A' bourke_mcgregor_scheme='B' mode_number_a=0 ; mode_number_b=0 if(mype+1 <= 2*nvmodes_keep) then if(mype+1 <= nvmodes_keep) then mode_number_a=mype+1 else mode_number_b=mype+1-nvmodes_keep end if end if !write(6,*)' in fmg_strong_bal_correction, mode_number_a,mode_number_b=',mode_number_a,mode_number_b ! if(mype > -10) then ! call mpi_finalize(i) ! stop ! end if ! vertical mode transform utilde=zero ; vtilde=zero ; mtilde=zero !!! redundant--remove eventually!!!! call vtrans(u_t,v_t,t_t,ps_t,utilde,vtilde,mtilde) ! filter low frequency scales out of tendencies call zrnmi_filter_uvm2(utilde,vtilde,mtilde,mype) ! layout for new rtlnmc forward code: !1. u0t,v0t,m0t input on subdomains -- output will be delchi, delpsi, delm-->delt,delps !2. set up slab areas A: 0 <= mype <= nvmodes_keep-1 ! B: nvmodes_keep <= mype <= 2*nvmodes_keep-1 ! require npe >= 2*nvmodes_keep !3. subdomains to slabs: u0t,v0t,m0t --> A and B u0t=zero ; v0t=zero ; m0t=zero ! may be redundant--consider removing!! call special_for_llfmg_sub2grid3(utilde,vtilde,mtilde,utilde,vtilde,mtilde,u0t,v0t,m0t,mype) !4. compute div0t on A from u0t,v0t and vort0t on B from u0t,v0t if(mode_number_a /= 0) then call get_div_reg(u0t,v0t,div0t) end if if(mode_number_b /= 0) then call get_vor_reg(u0t,v0t,vor0t) end if !5. solve HOP*delm = div0t on A ! HOP*mtg = delsqr*m0t-f*vor0t on B if(mode_number_a /= 0) then call fmg_wrapper_e(delm,div0t,.true.,nlon,nlat) end if if(mode_number_b /= 0) then call get_delsqr_reg(m0t,rhs) do j=1,nlon do i=1,nlat rhs(i,j)=rhs(i,j)-two*omega*sin(region_lat(i,j))*vor0t(i,j) end do end do call fmg_wrapper_e(mtg,rhs,.true.,nlon,nlat) end if !6. diagnostic BAL computation: ! BAL = integral( mtg_x**2 + mtg_y**2 + phi0*(vortg**2+divtg**2) ) ! bal_a = integral( phi0*divtg**2 ) ! bal_b = integral( mtg_x**2 + mtg_y**2 + phi0*vortg**2 ) ! BAL = bal_a + bal_b if(bal_diagnostic) then bal_a=zero ; bal_b=zero if(mode_number_a /= 0) then divtg=div0t bal_a=zero do j=1,nlon do i=1,nlat bal_a(mode_number_a)=bal_a(mode_number_a)+depths(mode_number_a)*divtg(i,j)**2 end do end do end if if(mode_number_b /= 0) then ! vortg=f*mtg/phi0 ; compute grad(mtg) do j=1,nlon do i=1,nlat vortg(i,j)=two*omega*sin(region_lat(i,j))*mtg(i,j)/depths(mode_number_b) end do end do call delx_reg(mtg,mtg_x,.false.) call dely_reg(mtg,mtg_y,.false.) bal_b=zero do j=1,nlon do i=1,nlat bal_b(mode_number_b)=bal_b(mode_number_b)+depths(mode_number_b)*vortg(i,j)**2 & +mtg_x(i,j)**2+mtg_y(i,j)**2 end do end do end if call mpi_allreduce(bal_a,bal_a0,nvmodes_keep,mpi_rtype,mpi_sum,mpi_comm_world,ierror) call mpi_allreduce(bal_b,bal_b0,nvmodes_keep,mpi_rtype,mpi_sum,mpi_comm_world,ierror) bal=bal_a0+bal_b0 if(mype == 0) then if(fullfield) then write(6,*)'FMG_STRONG_BAL: FULL FIELD BALANCE DIAGNOSTICS -- ' else write(6,*) 'FMG_STRONG_BAL: INCREMENTAL BALANCE DIAGNOSTICS -- ' end if do i=1,nvmodes_keep write(6,'(" jiter,mode,bal_a,b,bal=",2i3,3e20.9)')jiter,i,bal_a0(i),bal_b0(i),bal(i) end do end if end if !7. delvor = f*delm/phi0 --> solve delsqr(delpsi) = delvor on A ! deldiv = mtg/phi0 --> solve delsqr(delchi) = deldiv on B (for bourke_mcgregor_scheme = B ) ! deldiv = m0t/phi0 (for bourke_mcgregor_scheme = A ) if(mode_number_a /= 0) then do j=1,nlon do i=1,nlat delvor(i,j)=two*omega*sin(region_lat(i,j))*delm(i,j)/depths(mode_number_a) end do end do call fmg_wrapper_e(delpsi,delvor,.false.,nlon,nlat) end if if(mode_number_b /= 0) then if(bourke_mcgregor_scheme == 'A') then !write(6,*)' using bourke_mcgregor_scheme A, so deldiv = m0t/phi0' do j=1,nlon do i=1,nlat deldiv(i,j)=m0t(i,j)/depths(mode_number_b) end do end do elseif(bourke_mcgregor_scheme == 'B') then !write(6,*)' using bourke_mcgregor_scheme B, so deldiv = mtg/phi0' do j=1,nlon do i=1,nlat deldiv(i,j)=mtg(i,j)/depths(mode_number_b) end do end do end if call fmg_wrapper_e(delchi,deldiv,.false.,nlon,nlat) end if !8. delpsi,delm A slab --> subdomains ! delchi,dummy B slab --> subdomains ! if update_uv, u_psi,v_psi A slab --> subdomains ! if update_uv, u_chi,v_chi B slab --> subdomains delf1=zero ; delf2=zero if(mode_number_a /= 0) then do j=1,nlon do i=1,nlat delf1(i,j)= delpsi(i,j) delf2 (i,j)=delm(i,j) end do end do end if if(mode_number_b /= 0) then do j=1,nlon do i=1,nlat delf1(i,j)= delchi(i,j) delf2 (i,j)=delchi(i,j) end do end do end if call special_for_llfmg_grid2sub2(delf1,delf2,delpsitilde,delmtilde,delchitilde,dummytilde,mype) call vtrans_inv(delpsitilde,delchitilde,delmtilde,dpsi,dchi,dt,dps) if(update) then psi=psi+dpsi ; chi=chi+dchi end if if(update) then t=t+dt ; ps=ps+dps end if end subroutine fmg_strong_bal_correction subroutine fmg_strong_bal_correction_ad(u_t,v_t,t_t,ps_t,psi,chi,t,ps,update,mype) use kinds, only: r_kind,i_kind use constants, only: zero,two,omega use gridmod, only: lat2,lon2,nsig,nlat,nlon,region_lat use gridmod, only: istart,jstart use mod_vtrans, only: vtrans_ad,vtrans_inv_ad,nvmodes_keep,depths ! use helmholtz_fmg_mod, only: mg use mpimod, only: mpi_comm_world,ierror,mpi_sum,mpi_rtype use zrnmi_mod, only: zrnmi_filter_uvm2_ad use hybrid_ensemble_parameters, only: uv_hyb_ens implicit none integer(i_kind),intent(in)::mype logical,intent(in)::update real(r_kind),dimension(lat2,lon2,nsig),intent(inout)::u_t,v_t,t_t real(r_kind),dimension(lat2,lon2),intent(inout)::ps_t real(r_kind),dimension(lat2,lon2,nsig),intent(in)::psi,chi,t real(r_kind),dimension(lat2,lon2),intent(in)::ps real(r_kind),dimension(lat2,lon2,nvmodes_keep)::utilde,vtilde,mtilde,utdum,vtdum,mtdum real(r_kind),dimension(lat2,lon2,nvmodes_keep)::delpsitilde,delchitilde,delmtilde,dummytilde real(r_kind),dimension(nlat,nlon)::u0t,v0t,m0t,div0t,vor0t,rhs,mtg,delm,divtg,vortg,mtg_x,mtg_y real(r_kind),dimension(nlat,nlon)::delvor,deldiv,delpsi,delchi,u_psi,v_psi,u_chi,v_chi,delf1,delf2 real(r_kind),dimension(lat2,lon2,nsig)::dpsi,dchi,dt real(r_kind),dimension(lat2,lon2)::dps integer(i_kind) i,j,mode_number_a,mode_number_b real(r_kind) bal_a (nvmodes_keep),bal_b (nvmodes_keep) real(r_kind) bal_a0(nvmodes_keep),bal_b0(nvmodes_keep) real(r_kind) bal(nvmodes_keep) real(r_kind),dimension(lat2,lon2,nvmodes_keep)::delupsitilde,delvpsitilde,deluchitilde,delvchitilde real(r_kind),dimension(lat2,lon2,nvmodes_keep)::delutilde,delvtilde character(1) bourke_mcgregor_scheme if(uv_hyb_ens) then write(0,*)' STOP IN fmg_strong_bal_correction--NOT ABLE TO ACCEPT uv_hyb_ens=.true. YET' call stop2(998) end if !bourke_mcgregor_scheme='A' bourke_mcgregor_scheme='B' mode_number_a=0 ; mode_number_b=0 if(mype+1 <= 2*nvmodes_keep) then if(mype+1 <= nvmodes_keep) then mode_number_a=mype+1 else mode_number_b=mype+1-nvmodes_keep end if end if !write(6,*)' in fmg_strong_bal_correction_ad, mode_number_a,mode_number_b=',mode_number_a,mode_number_b dpsi=zero ; dchi=zero ; dt=zero ; dps=zero if(update) then dpsi=psi ; dchi=chi dt=t ; dps=ps end if delpsitilde=zero ; delchitilde=zero ; delmtilde=zero call vtrans_inv_ad(delpsitilde,delchitilde,delmtilde,dpsi,dchi,dt,dps) !8. delpsi,delm A slab --> subdomains ! delchi,dummy B slab --> subdomains ! if update_uv, u_psi,v_psi A slab --> subdomains ! if update_uv, u_chi,v_chi B slab --> subdomains call special_for_llfmg_sub2grid2(delpsitilde,delmtilde,delchitilde,dummytilde,delf1,delf2,mype) delpsi=zero ; delchi=zero ; delm=zero if(mode_number_a /= 0) then do j=1,nlon do i=1,nlat delpsi(i,j)=delf1(i,j) delm(i,j)=delf2(i,j) end do end do end if if(mode_number_b /= 0) then do j=1,nlon do i=1,nlat delchi(i,j)=delf1(i,j) end do end do end if !7. delvor = f*delm/phi0 --> solve delsqr(delpsi) = delvor on A ! deldiv = mtg/phi0 --> solve delsqr(delchi) = deldiv on B (for bourke_mcgregor_scheme = B ) ! deldiv = m0t/phi0 (for bourke_mcgregor_scheme = A ) if(mode_number_a /= 0) then delvor=zero call fmg_wrapper_e_ad(delpsi,delvor,.false.,nlon,nlat) do j=1,nlon do i=1,nlat delm(i,j)=delm(i,j)+two*omega*sin(region_lat(i,j))*delvor(i,j)/depths(mode_number_a) end do end do end if if(mode_number_b /= 0) then deldiv=zero call fmg_wrapper_e_ad(delchi,deldiv,.false.,nlon,nlat) m0t=zero ; mtg=zero if(bourke_mcgregor_scheme == 'A') then !write(6,*)' using bourke_mcgregor_scheme A, so deldiv = m0t/phi0' do j=1,nlon do i=1,nlat m0t(i,j)=deldiv(i,j)/depths(mode_number_b) end do end do elseif(bourke_mcgregor_scheme == 'B') then !write(6,*)' using bourke_mcgregor_scheme B, so deldiv = mtg/phi0' do j=1,nlon do i=1,nlat mtg(i,j)=deldiv(i,j)/depths(mode_number_b) end do end do end if end if !5. solve HOP*delm = div0t on A ! HOP*mtg = delsqr*m0t-f*vor0t on B if(mode_number_a /= 0) then div0t=zero call fmg_wrapper_e_ad(delm,div0t,.true.,nlon,nlat) end if if(mode_number_b /= 0) then rhs=zero call fmg_wrapper_e_ad(mtg,rhs,.true.,nlon,nlat) do j=1,nlon do i=1,nlat vor0t(i,j)= -two*omega*sin(region_lat(i,j))*rhs(i,j) end do end do call tget_delsqr_reg(rhs,m0t) end if !4. compute div0t on A from u0t,v0t and vort0t on B from u0t,v0t if(mode_number_a /= 0) then call get_div_reg_ad(u0t,v0t,div0t) end if if(mode_number_b /= 0) then call get_vor_reg_ad(u0t,v0t,vor0t) end if !3. subdomains to slabs: u0t,v0t,m0t --> A and B utilde=zero ; vtilde=zero ; mtilde=zero call special_for_llfmg_grid2sub3(u0t,v0t,m0t,utilde,vtdum,mtdum,utdum,vtilde,mtilde,mype) utilde=utilde+utdum ; vtilde=vtilde+vtdum ! adjoint of filter low frequency scales out of tendencies call zrnmi_filter_uvm2_ad(utilde,vtilde,mtilde,mype) ! adjoint of vertical mode transform call vtrans_ad(u_t,v_t,t_t,ps_t,utilde,vtilde,mtilde) end subroutine fmg_strong_bal_correction_ad subroutine fmg_strong_bal_correction_ad_test(u_t,v_t,t_t,ps_t,psi,chi,t,ps,mype) !<><><><<><><><>test of fmg_strong_bal_correction_ad><>><<><><<<<<<<<<<<<<<<<>>>>>>>>>>>>> use constants, only: zero,two use mpimod, only: mpi_rtype,mpi_sum,mpi_comm_world,ierror use kinds, only: r_kind,i_kind use gridmod, only: lat2,lon2,nsig implicit none integer(i_kind),intent(in)::mype real(r_kind),dimension(lat2,lon2,nsig),intent(inout)::u_t,v_t,t_t real(r_kind),dimension(lat2,lon2),intent(inout)::ps_t real(r_kind),dimension(lat2,lon2,nsig),intent(inout)::psi,chi,t real(r_kind),dimension(lat2,lon2),intent(inout)::ps real(r_kind),dimension(lat2,lon2,nsig)::uu_t,vv_t,tt_t real(r_kind),dimension(lat2,lon2)::pps_t real(r_kind),dimension(lat2,lon2,nsig)::psi2,chi2,t2 real(r_kind),dimension(lat2,lon2)::ps2 integer(i_kind) ivar,i,j,k real(r_kind) yty,yty0,xtz,xtz0,errmax errmax=zero do ivar=0,4 if(ivar > 0) then u_t=zero ; v_t=zero ; t_t=zero ; ps_t=zero if(ivar == 1) call random_number(u_t) if(ivar == 2) call random_number(v_t) if(ivar == 3) call random_number(t_t) if(ivar == 4) call random_number(ps_t) end if uu_t=u_t ; vv_t=v_t ; tt_t=t_t ; pps_t=ps_t psi=zero ; chi=zero ; t=zero ; ps=zero call fmg_strong_bal_correction(u_t,v_t,t_t,ps_t,psi,chi,t,ps,.false.,.false.,.true.,mype) yty=zero do k=1,nsig do i=2,lon2-1 do j=2,lat2-1 yty=yty+psi(j,i,k)**2+chi(j,i,k)**2+t(j,i,k)**2 end do end do end do do i=2,lon2-1 do j=2,lat2-1 yty=yty+ps(j,i)**2 end do end do call mpi_allreduce(yty,yty0,1,mpi_rtype,mpi_sum,mpi_comm_world,ierror) u_t=zero ; v_t=zero ; t_t=zero ; ps_t=zero call fmg_strong_bal_correction_ad(u_t,v_t,t_t,ps_t,psi,chi,t,ps,.true.,mype) xtz=zero do k=1,nsig do i=2,lon2-1 do j=2,lat2-1 xtz=xtz+uu_t(j,i,k)*u_t(j,i,k)+vv_t(j,i,k)*v_t(j,i,k)+tt_t(j,i,k)*t_t(j,i,k) end do end do end do do i=2,lon2-1 do j=2,lat2-1 xtz=xtz+pps_t(j,i)*ps_t(j,i) end do end do call mpi_allreduce(xtz,xtz0,1,mpi_rtype,mpi_sum,mpi_comm_world,ierror) if(mype == 0) & write(6,'(" fmg_strong_bal_correction_ad (aT*a), ivar,xx,yy=",i2,2e22.15,e10.3)') & ivar,xtz0,yty0,abs(two*(xtz0-yty0)/(abs(xtz0)+abs(yty0))) if(abs(xtz0)+abs(yty0) /= zero) errmax=max(abs(two*(xtz0-yty0)/(abs(xtz0)+abs(yty0))),errmax) if(ivar == 0) cycle psi=zero ; chi=zero ; t=zero ; ps=zero if(ivar == 1) call random_number(psi) if(ivar == 2) call random_number(chi) if(ivar == 3) call random_number(t) if(ivar == 4) call random_number(ps) psi2=psi ; chi2=chi ; t2=t ; ps2=ps u_t=zero ; v_t=zero ; t_t=zero ; ps_t=zero call fmg_strong_bal_correction_ad(u_t,v_t,t_t,ps_t,psi,chi,t,ps,.true.,mype) yty=zero do k=1,nsig do i=2,lon2-1 do j=2,lat2-1 yty=yty+u_t(j,i,k)*u_t(j,i,k)+v_t(j,i,k)*v_t(j,i,k)+t_t(j,i,k)*t_t(j,i,k) end do end do end do do i=2,lon2-1 do j=2,lat2-1 yty=yty+ps_t(j,i)*ps_t(j,i) end do end do call mpi_allreduce(yty,yty0,1,mpi_rtype,mpi_sum,mpi_comm_world,ierror) psi=zero ; chi=zero ; t=zero ; ps=zero call fmg_strong_bal_correction(u_t,v_t,t_t,ps_t,psi,chi,t,ps,.false.,.false.,.true.,mype) xtz=zero do k=1,nsig do i=2,lon2-1 do j=2,lat2-1 xtz=xtz+psi2(j,i,k)*psi(j,i,k)+chi2(j,i,k)*chi(j,i,k)+t2(j,i,k)*t(j,i,k) end do end do end do do i=2,lon2-1 do j=2,lat2-1 xtz=xtz+ps2(j,i)*ps(j,i) end do end do call mpi_allreduce(xtz,xtz0,1,mpi_rtype,mpi_sum,mpi_comm_world,ierror) if(mype == 0) & write(6,'(" fmg_strong_bal_correction_ad (a*aT), ivar,xx,yy=",i2,2e22.15,e10.3)') & ivar,xtz0,yty0,abs(two*(xtz0-yty0)/(abs(xtz0)+abs(yty0))) if(abs(xtz0)+abs(yty0) /= zero) errmax=max(abs(two*(xtz0-yty0)/(abs(xtz0)+abs(yty0))),errmax) end do if(mype == 0) write(6,'(" max error for fmg_strong_bal_correction_ad=",e10.3)') errmax if(mype > -1000) then call mpi_finalize(ierror) stop end if !<><><><<><><><><<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>> end subroutine fmg_strong_bal_correction_ad_test subroutine zrnmi_filter_uvm_ad_test(mype) !<><><><<><><><>test of zrnmi_filter_uvm_ad_test><>><<><><<<<<<<<<<<<<<<<>>>>>>>>>>>>> use constants, only: zero,two use mpimod, only: mpi_rtype,mpi_sum,mpi_comm_world,ierror use kinds, only: r_kind,i_kind use gridmod, only: lat2,lon2,nsig use mod_vtrans, only: nvmodes_keep use zrnmi_mod, only: zrnmi_filter_uvm2,zrnmi_filter_uvm2_ad implicit none integer(i_kind),intent(in)::mype real(r_kind),dimension(lat2,lon2,nvmodes_keep)::u1,u2,v1,v2,m1,m2 integer(i_kind) ivar,i,j,k real(r_kind) yty,yty0,xtz,xtz0,errmax errmax=zero do ivar=1,3 u1=zero ; v1=zero ; m1=zero if(ivar == 1) call random_number(u1) if(ivar == 2) call random_number(v1) if(ivar == 3) call random_number(m1) u2=u1 ; v2=v1 ; m2=m1 call zrnmi_filter_uvm2(u1,v1,m1,mype) yty=zero do k=1,nvmodes_keep do i=2,lon2-1 do j=2,lat2-1 yty=yty+u1(j,i,k)**2+v1(j,i,k)**2+m1(j,i,k)**2 end do end do end do call mpi_allreduce(yty,yty0,1,mpi_rtype,mpi_sum,mpi_comm_world,ierror) call zrnmi_filter_uvm2_ad(u1,v1,m1,mype) xtz=zero do k=1,nvmodes_keep do i=2,lon2-1 do j=2,lat2-1 xtz=xtz+u1(j,i,k)*u2(j,i,k)+v1(j,i,k)*v2(j,i,k)+m1(j,i,k)*m2(j,i,k) end do end do end do call mpi_allreduce(xtz,xtz0,1,mpi_rtype,mpi_sum,mpi_comm_world,ierror) if(mype == 0) & write(6,'(" zrnmi_filter_uvm2_ad (aT*a), ivar,xx,yy=",i2,2e22.15,e10.3)') & ivar,xtz0,yty0,abs(two*(xtz0-yty0)/(abs(xtz0)+abs(yty0))) if(abs(xtz0)+abs(yty0) /= zero) errmax=max(abs(two*(xtz0-yty0)/(abs(xtz0)+abs(yty0))),errmax) u1=zero ; v1=zero ; m1=zero if(ivar == 1) call random_number(u1) if(ivar == 2) call random_number(v1) if(ivar == 3) call random_number(m1) u2=u1 ; v2=v1 ; m2=m1 call zrnmi_filter_uvm2_ad(u1,v1,m1,mype) yty=zero do k=1,nvmodes_keep do i=2,lon2-1 do j=2,lat2-1 yty=yty+u1(j,i,k)**2+v1(j,i,k)**2+m1(j,i,k)**2 end do end do end do call mpi_allreduce(yty,yty0,1,mpi_rtype,mpi_sum,mpi_comm_world,ierror) call zrnmi_filter_uvm2(u1,v1,m1,mype) xtz=zero do k=1,nvmodes_keep do i=2,lon2-1 do j=2,lat2-1 xtz=xtz+u1(j,i,k)*u2(j,i,k)+v1(j,i,k)*v2(j,i,k)+m1(j,i,k)*m2(j,i,k) end do end do end do call mpi_allreduce(xtz,xtz0,1,mpi_rtype,mpi_sum,mpi_comm_world,ierror) if(mype == 0) & write(6,'(" zrnmi_filter_uvm2_ad (a*aT), ivar,xx,yy=",i2,2e22.15,e10.3)') & ivar,xtz0,yty0,abs(two*(xtz0-yty0)/(abs(xtz0)+abs(yty0))) if(abs(xtz0)+abs(yty0) /= zero) errmax=max(abs(two*(xtz0-yty0)/(abs(xtz0)+abs(yty0))),errmax) end do if(mype == 0) write(6,'(" max error for zrnmi_filter_uvm2_ad=",e10.3)') errmax if(mype > -1000) then call mpi_finalize(ierror) stop end if !<><><><<><><><><<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>> end subroutine zrnmi_filter_uvm_ad_test subroutine fmg_wrapper(v,f,helmholtz,nlon,nlat) use kinds, only: r_kind,i_kind use helmholtz_fmg_mod,only: fmg implicit none integer(i_kind),intent(in):: nlon,nlat real(r_kind),intent(in):: f(0:nlat-1,0:nlon-1) real(r_kind),intent(out):: v(0:nlat-1,0:nlon-1) logical,intent(in):: helmholtz call fmg(v,f,helmholtz,1,1,60,nlon-1,nlat-1) end subroutine fmg_wrapper subroutine fmg_wrapper_e(v,f,helmholtz,nlon,nlat) use kinds, only: r_kind,i_kind use constants, only: zero use helmholtz_fmg_mod,only: fmg implicit none integer(i_kind),intent(in):: nlon,nlat real(r_kind),intent(in):: f(nlat,nlon) real(r_kind),intent(out):: v(nlat,nlon) logical,intent(in):: helmholtz real(r_kind),dimension(0:nlat+1,0:nlon+1)::fe,ve integer(i_kind) i,j fe=zero do j=1,nlon do i=1,nlat fe(i,j)=f(i,j) end do end do call fmg(ve,fe,helmholtz,1,1,60,nlon+1,nlat+1) do j=1,nlon do i=1,nlat v(i,j)=ve(i,j) end do end do end subroutine fmg_wrapper_e subroutine fmg_wrapper_e_ad(v,f,helmholtz,nlon,nlat) use kinds, only: r_kind,i_kind use constants, only: zero use helmholtz_fmg_mod,only: fmg_ad implicit none integer(i_kind),intent(in):: nlon,nlat real(r_kind),intent(inout):: f(nlat,nlon) real(r_kind),intent(in):: v(nlat,nlon) logical,intent(in):: helmholtz real(r_kind),dimension(0:nlat+1,0:nlon+1)::fe,ve integer(i_kind) i,j ve=zero do j=1,nlon do i=1,nlat ve(i,j)=v(i,j) fe(i,j)=f(i,j) end do end do call fmg_ad(ve,fe,helmholtz,1,1,60,nlon+1,nlat+1) do j=1,nlon do i=1,nlat f(i,j)=fe(i,j) end do end do end subroutine fmg_wrapper_e_ad subroutine get_delsqr_reg(work1,work2) !$$$ subprogram documentation block ! . . . . ! subprogram: get_delsqr_reg compute laplacian ! prgmmr: parrish org: np23 date: 2006-02-13 ! ! abstract: compute laplacian ! ! program history log: ! 2006-02-13 parrish ! ! input argument list: ! work1 - array to be differentiated ! ! output argument list: ! work2 - array containing laplacian ! ! attributes: ! language: f90 ! machine: ibm rs/6000 sp !$$$ use kinds, only: r_kind,i_kind use gridmod, only: nlat,nlon implicit none integer(i_kind) i,j real(r_kind),dimension(nlat,nlon):: work1,work2 real(r_kind),dimension(nlat,nlon):: grid1,grid2,grid3,grid4 call delx_reg(work1,grid1,(.false.)) call delx_reg(grid1,grid2,(.true.)) call dely_reg(work1,grid3,(.false.)) call dely_reg(grid3,grid4,(.true.)) do j=1,nlon do i=1,nlat work2(i,j)=grid2(i,j)+grid4(i,j) end do end do return end subroutine get_delsqr_reg subroutine tget_delsqr_reg(work2,work1) !$$$ subprogram documentation block ! . . . . ! subprogram: tget_delsqr_reg adjoint of get_delsqr_reg ! prgmmr: parrish org: np23 date: 2006-02-13 ! ! abstract: adjoint of get_delsqr_reg ! ! program history log: ! 2006-02-13 parrish ! ! input argument list: ! work2 - array containing delsqr variable ! work1 - previous contents to be accumulated to ! ! output argument list: ! work1 - array containing accumulation of adjoint delsqr ! ! attributes: ! language: f90 ! machine: ibm rs/6000 sp !$$$ use kinds, only: r_kind,i_kind use gridmod, only: nlat,nlon use constants, only: zero implicit none integer(i_kind) i,j real(r_kind),dimension(nlat,nlon):: work1,work2 real(r_kind),dimension(nlat,nlon):: grid1,grid2,grid3,grid4 do j=1,nlon do i=1,nlat grid2(i,j)=zero grid4(i,j)=zero grid1(i,j)=zero grid3(i,j)=zero end do end do do j=1,nlon do i=1,nlat grid2(i,j)=grid2(i,j)+work2(i,j) grid4(i,j)=grid4(i,j)+work2(i,j) end do end do call tdely_reg(grid4,grid3,(.true.)) call tdely_reg(grid3,work1,(.false.)) call tdelx_reg(grid2,grid1,(.true.)) call tdelx_reg(grid1,work1,(.false.)) return end subroutine tget_delsqr_reg subroutine outgrad1(f,label,nx,ny) character(*) label integer(4) nx,ny real(8) f(ny,nx) character(80) dsdes,dsdat character(80) datdes(1000) character(1) blank data blank/' '/ data undef/-9.99e33/ real(4) out(nx,ny) np=1 ioutcor=10 ioutdat=11 write(dsdes,'(a,".des")')trim(label) write(dsdat,'(a,".dat")')trim(label) open(unit=ioutcor,file=dsdes,form='formatted') open(unit=ioutdat,file=dsdat,form='unformatted') ntime=1 rlonmap0=1. dlonmap=1. rlatmap0=1. dlatmap=1. startp=1. pinc=1. koutmax=1 do i=1,1000 write(datdes(i),'(80a1)')(blank,k=1,80) end do write(datdes(1),'("DSET ",a)')trim(dsdat) write(datdes(2),'("options big_endian sequential")') write(datdes(3),'("TITLE ",a)')trim(label) write(datdes(4),'("UNDEF ",e11.2)')undef write(datdes(5),'("XDEF ",i5," LINEAR ",f7.2,f7.2)')nx,rlonmap0,dlonmap write(datdes(6),'("YDEF ",i5," LINEAR ",f7.2,f7.2)')ny,rlatmap0,dlatmap next=7 write(datdes(next),'("ZDEF ",i5," LINEAR ",f7.2,f7.2)')np,startp,pinc next=next+1 write(datdes(next),'("TDEF ",i5," LINEAR 0Z23may1992 24hr")')koutmax next=next+1 write(datdes(next),'("VARS 1")') next=next+1 write(datdes(next),'("f ",i5," 99 f ")')np next=next+1 write(datdes(next),'("ENDVARS")') last=next write(ioutcor,'(a80)')(datdes(i),i=1,last) close(ioutcor) do i=1,nx do j=1,ny out(i,j)=f(j,i) end do end do write(ioutdat)((out(i,j),i=1,nx),j=1,ny) close(ioutdat) return end subroutine outgrad1 subroutine special_for_llfmg_sub2grid3(a1,a2,a3,b1,b2,b3,f1,f2,f3,mype) use kinds, only: i_kind,r_kind use mpimod, only: npe use gridmod, only: lat1,lon1,lat2,lon2,nlat,nlon,itotsub,iglobal,ltosi,ltosj use mod_vtrans, only: nvmodes_keep implicit none integer(i_kind),intent(in):: mype real(r_kind),dimension(lat2,lon2,nvmodes_keep),intent(in)::a1,a2,a3,b1,b2,b3 real(r_kind),dimension(nlat,nlon),intent(out):: f1,f2,f3 real(r_kind) all_loc(lat1,lon1,3,2*nvmodes_keep),tempa(itotsub,3) integer(i_kind) i,ip1,j,jp1,k,ka,kb integer(i_kind) kbegin(0:npe-1),kend(0:npe-1),numlevs(0:npe-1) do k=1,nvmodes_keep ka=k kb=ka+nvmodes_keep do i=1,lon1 ip1=i+1 do j=1,lat1 jp1=j+1 all_loc(j,i,1,ka)=a1(jp1,ip1,k) all_loc(j,i,2,ka)=a2(jp1,ip1,k) all_loc(j,i,3,ka)=a3(jp1,ip1,k) all_loc(j,i,1,kb)=b1(jp1,ip1,k) all_loc(j,i,2,kb)=b2(jp1,ip1,k) all_loc(j,i,3,kb)=b3(jp1,ip1,k) end do end do end do do k=0,2*nvmodes_keep-1 numlevs(k)=3 end do do k=2*nvmodes_keep,npe-1 numlevs(k)=0 end do kbegin(0)=1 do k=1,npe-1 kbegin(k)=kbegin(k-1)+numlevs(k-1) end do do k=0,npe-1 kend(k)=kbegin(k)+numlevs(k)-1 end do call generic_sub2grid8(all_loc,tempa,kbegin(mype),kend(mype),kbegin,kend,mype,6*nvmodes_keep) do k=1,kend(mype)-kbegin(mype)+1 if(k == 1) then do i=1,iglobal f1(ltosi(i),ltosj(i))=tempa(i,k) end do elseif(k == 2) then do i=1,iglobal f2(ltosi(i),ltosj(i))=tempa(i,k) end do elseif(k == 3) then do i=1,iglobal f3(ltosi(i),ltosj(i))=tempa(i,k) end do end if end do end subroutine special_for_llfmg_sub2grid3 subroutine special_for_llfmg_sub2grid2(a1,a2,b1,b2,f1,f2,mype) use kinds, only: i_kind,r_kind use mpimod, only: npe use gridmod, only: lat1,lon1,lat2,lon2,nlat,nlon,itotsub,iglobal,ltosi,ltosj use mod_vtrans, only: nvmodes_keep implicit none integer(i_kind),intent(in):: mype real(r_kind),dimension(lat2,lon2,nvmodes_keep),intent(in)::a1,a2,b1,b2 real(r_kind),dimension(nlat,nlon),intent(out):: f1,f2 real(r_kind) all_loc(lat1,lon1,2,2*nvmodes_keep),tempa(itotsub,2) integer(i_kind) i,ip1,j,jp1,k,ka,kb integer(i_kind) kbegin(0:npe-1),kend(0:npe-1),numlevs(0:npe-1) do k=1,nvmodes_keep ka=k kb=ka+nvmodes_keep do i=1,lon1 ip1=i+1 do j=1,lat1 jp1=j+1 all_loc(j,i,1,ka)=a1(jp1,ip1,k) all_loc(j,i,2,ka)=a2(jp1,ip1,k) all_loc(j,i,1,kb)=b1(jp1,ip1,k) all_loc(j,i,2,kb)=b2(jp1,ip1,k) end do end do end do do k=0,2*nvmodes_keep-1 numlevs(k)=2 end do do k=2*nvmodes_keep,npe-1 numlevs(k)=0 end do kbegin(0)=1 do k=1,npe-1 kbegin(k)=kbegin(k-1)+numlevs(k-1) end do do k=0,npe-1 kend(k)=kbegin(k)+numlevs(k)-1 end do call generic_sub2grid8(all_loc,tempa,kbegin(mype),kend(mype),kbegin,kend,mype,4*nvmodes_keep) do k=1,kend(mype)-kbegin(mype)+1 if(k == 1) then do i=1,iglobal f1(ltosi(i),ltosj(i))=tempa(i,k) end do elseif(k == 2) then do i=1,iglobal f2(ltosi(i),ltosj(i))=tempa(i,k) end do end if end do end subroutine special_for_llfmg_sub2grid2 subroutine special_for_llfmg_grid2sub2(f1,f2,a1,a2,b1,b2,mype) use kinds, only: i_kind,r_kind use mpimod, only: npe use gridmod, only: lat1,lon1,lat2,lon2,nlat,nlon,itotsub,ltosi_s,ltosj_s use mod_vtrans, only: nvmodes_keep implicit none integer(i_kind),intent(in):: mype real(r_kind),dimension(nlat,nlon),intent(in):: f1,f2 real(r_kind),dimension(lat2,lon2,nvmodes_keep),intent(out)::a1,a2,b1,b2 real(r_kind) all_loc(lat2,lon2,2,2*nvmodes_keep),tempa(itotsub,2) integer(i_kind) i,j,k,ka,kb integer(i_kind) kbegin(0:npe-1),kend(0:npe-1),numlevs(0:npe-1) do k=0,2*nvmodes_keep-1 numlevs(k)=2 end do do k=2*nvmodes_keep,npe-1 numlevs(k)=0 end do kbegin(0)=1 do k=1,npe-1 kbegin(k)=kbegin(k-1)+numlevs(k-1) end do do k=0,npe-1 kend(k)=kbegin(k)+numlevs(k)-1 end do do k=1,kend(mype)-kbegin(mype)+1 if(k == 1) then do i=1,itotsub tempa(i,k)=f1(ltosi_s(i),ltosj_s(i)) end do elseif(k == 2) then do i=1,itotsub tempa(i,k)=f2(ltosi_s(i),ltosj_s(i)) end do end if end do call generic_grid2sub8(tempa,all_loc,kbegin(mype),kend(mype),kbegin,kend,mype,4*nvmodes_keep) do k=1,nvmodes_keep ka=k kb=ka+nvmodes_keep do i=1,lon2 do j=1,lat2 a1(j,i,k)=all_loc(j,i,1,ka) a2(j,i,k)=all_loc(j,i,2,ka) b1(j,i,k)=all_loc(j,i,1,kb) b2(j,i,k)=all_loc(j,i,2,kb) end do end do end do end subroutine special_for_llfmg_grid2sub2 subroutine special_for_llfmg_grid2sub3(f1,f2,f3,a1,a2,a3,b1,b2,b3,mype) use kinds, only: i_kind,r_kind use mpimod, only: npe use gridmod, only: lat1,lon1,lat2,lon2,nlat,nlon,itotsub,ltosi_s,ltosj_s use mod_vtrans, only: nvmodes_keep implicit none integer(i_kind),intent(in):: mype real(r_kind),dimension(nlat,nlon),intent(in):: f1,f2,f3 real(r_kind),dimension(lat2,lon2,nvmodes_keep),intent(out)::a1,a2,a3,b1,b2,b3 real(r_kind) all_loc(lat2,lon2,3,2*nvmodes_keep),tempa(itotsub,3) integer(i_kind) i,j,k,ka,kb integer(i_kind) kbegin(0:npe-1),kend(0:npe-1),numlevs(0:npe-1) do k=0,2*nvmodes_keep-1 numlevs(k)=3 end do do k=2*nvmodes_keep,npe-1 numlevs(k)=0 end do kbegin(0)=1 do k=1,npe-1 kbegin(k)=kbegin(k-1)+numlevs(k-1) end do do k=0,npe-1 kend(k)=kbegin(k)+numlevs(k)-1 end do do k=1,kend(mype)-kbegin(mype)+1 if(k == 1) then do i=1,itotsub tempa(i,k)=f1(ltosi_s(i),ltosj_s(i)) end do elseif(k == 2) then do i=1,itotsub tempa(i,k)=f2(ltosi_s(i),ltosj_s(i)) end do elseif(k == 3) then do i=1,itotsub tempa(i,k)=f3(ltosi_s(i),ltosj_s(i)) end do end if end do call generic_grid2sub8(tempa,all_loc,kbegin(mype),kend(mype),kbegin,kend,mype,6*nvmodes_keep) do k=1,nvmodes_keep ka=k kb=ka+nvmodes_keep do i=1,lon2 do j=1,lat2 a1(j,i,k)=all_loc(j,i,1,ka) a2(j,i,k)=all_loc(j,i,2,ka) a3(j,i,k)=all_loc(j,i,3,ka) b1(j,i,k)=all_loc(j,i,1,kb) b2(j,i,k)=all_loc(j,i,2,kb) b3(j,i,k)=all_loc(j,i,3,kb) end do end do end do end subroutine special_for_llfmg_grid2sub3 subroutine generic_sub2grid8(all_loc,tempa,kbegin_loc,kend_loc,kbegin,kend,mype,num_fields) !$$$ subprogram documentation block ! . . . . ! subprogram: generic_sub2grid converts from subdomains to full horizontal grid ! prgmmr: parrish org: np22 date: 2004-11-29 ! ! abstract: variation on subroutine sub2grid, with more general distribution of variables ! along the k index. ! ! program history log: ! 2004-02-03 kleist, new mpi strategy ! 2004-05-06 derber ! 2004-07-15 treadon - handle periodic subdomains ! 2004-07-28 treadon - add only on use declarations; add intent in/out ! 2004-10-26 kleist - u,v removed; periodicity accounted for only in ! sub2grid routine if necessary ! 2004-11-29 parrish - adapt sub2grid for related use with mpi io. ! ! input argument list: ! all_loc - input grid values in vertical subdomain mode ! kbegin_loc - starting k index for tempa on local processor ! kend_loc - ending k index for tempa on local processor ! kbegin - starting k indices for tempa for all processors ! kend - ending k indices for tempa for all processors ! mype - local processor number ! num_fields - total range of k index (1 <= k <= num_fields) ! ! output argument list: ! tempa - output grid values in horizontal slab mode. ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use mpimod, only: ierror,mpi_comm_world,mpi_rtype,npe use gridmod, only: ijn,itotsub,lat1,lon1 use kinds, only: r_kind,i_kind use constants, only: zero implicit none integer(i_kind) kbegin_loc,kend_loc,mype,num_fields integer(i_kind) kbegin(0:npe),kend(0:npe-1) real(r_kind) tempa(itotsub,kbegin_loc:max(kbegin_loc,kend_loc)) real(r_kind) all_loc(lat1*lon1*num_fields) integer(i_kind) k integer(i_kind) sendcounts(0:npe-1),sdispls(0:npe),recvcounts(0:npe-1),rdispls(0:npe) ! first get alltoallv indices sdispls(0)=0 do k=0,npe-1 sendcounts(k)=ijn(k+1)*(kend_loc-kbegin_loc+1) sdispls(k+1)=sdispls(k)+sendcounts(k) end do rdispls(0)=0 do k=0,npe-1 recvcounts(k)=ijn(mype+1)*(kend(k)-kbegin(k)+1) rdispls(k+1)=rdispls(k)+recvcounts(k) end do call mpi_alltoallv(all_loc,recvcounts,rdispls,mpi_rtype, & tempa,sendcounts,sdispls,mpi_rtype,mpi_comm_world,ierror) if(kbegin_loc <= kend_loc) then call reorder_s8(tempa,kend_loc-kbegin_loc+1) else tempa=zero end if end subroutine generic_sub2grid8 subroutine reorder_s8(work,k_in) ! !USES: use mpimod, only: npe use kinds, only: r_kind,i_kind use constants, only: zero use gridmod, only: ijn,itotsub implicit none ! !INPUT PARAMETERS: integer(i_kind), intent(in) :: k_in ! number of levs in work array ! !INPUT/OUTPUT PARAMETERS: real(r_kind),dimension(itotsub*k_in),intent(inout):: work ! array to reorder ! !OUTPUT PARAMETERS: ! !DESCRIPTION: adapt reorder to work with single precision ! ! !REVISION HISTORY: ! ! 2004-01-25 kleist ! 2004-05-14 kleist, documentation ! 2004-07-15 todling, protex-complaint prologue ! 2004-11-29 adapt reorder to work with single precision ! ! !REMAKRS: ! ! language: f90 ! machine: ibm rs/6000 sp; sgi origin 2000; compaq/hp ! ! !AUTHOR: ! kleist org: np20 date: 2004-01-25 ! !EOP !------------------------------------------------------------------------- integer(i_kind) iloc,iskip,i,j,k,n real(r_kind),dimension(itotsub,k_in):: temp ! Zero out temp array do k=1,k_in do i=1,itotsub temp(i,k)=zero end do end do ! Load temp array in desired order do k=1,k_in iskip=0 iloc=0 do n=1,npe if (n/=1) then iskip=iskip+ijn(n-1)*k_in end if do i=1,ijn(n) iloc=iloc+1 temp(iloc,k)=work(i + iskip + & (k-1)*ijn(n)) end do end do end do ! Load the temp array back into work iloc=0 do k=1,k_in do i=1,itotsub iloc=iloc+1 work(iloc)=temp(i,k) end do end do return end subroutine reorder_s8 subroutine generic_grid2sub8(tempa,all_loc,kbegin_loc,kend_loc,kbegin,kend,mype,num_fields) !$$$ subprogram documentation block ! . . . . ! subprogram: generic_grid2sub converts from full horizontal grid to subdomains ! prgmmr: parrish org: np22 date: 2004-11-29 ! ! abstract: variation on subroutine grid2sub, with more general distribution of variables ! along the k index. ! ! program history log: ! 2004-02-03 kleist, new mpi strategy ! 2004-05-06 derber ! 2004-07-15 treadon - handle periodic subdomains ! 2004-07-28 treadon - add only on use declarations; add intent in/out ! 2004-10-26 kleist - u,v removed; periodicity accounted for only in ! sub2grid routine if necessary ! 2004-11-29 parrish - adapt grid2sub for related use with mpi io. ! ! input argument list: ! tempa - input grid values in horizontal slab mode. ! kbegin_loc - starting k index for tempa on local processor ! kend_loc - ending k index for tempa on local processor ! kbegin - starting k indices for tempa for all processors ! kend - ending k indices for tempa for all processors ! mype - local processor number ! num_fields - total range of k index (1 <= k <= num_fields) ! ! output argument list: ! all_loc - output grid values in vertical subdomain mode ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use mpimod, only: ierror,mpi_comm_world,mpi_rtype,npe use gridmod, only: ijn_s,itotsub,lat2,lon2 use kinds, only: r_kind,i_kind use constants, only: zero implicit none integer(i_kind) kbegin_loc,kend_loc,mype,num_fields integer(i_kind) kbegin(0:npe),kend(0:npe-1) real(r_kind) tempa(itotsub,kbegin_loc:max(kbegin_loc,kend_loc)) real(r_kind) all_loc(lat2*lon2*num_fields) integer(i_kind) k integer(i_kind) sendcounts(0:npe-1),sdispls(0:npe),recvcounts(0:npe-1),rdispls(0:npe) ! first get alltoallv indices sdispls(0)=0 do k=0,npe-1 sendcounts(k)=ijn_s(k+1)*(kend_loc-kbegin_loc+1) sdispls(k+1)=sdispls(k)+sendcounts(k) end do rdispls(0)=0 do k=0,npe-1 recvcounts(k)=ijn_s(mype+1)*(kend(k)-kbegin(k)+1) rdispls(k+1)=rdispls(k)+recvcounts(k) end do ! then call reorder2 if(kbegin_loc <= kend_loc) then call reorder2_s8(tempa,kend_loc-kbegin_loc+1) else tempa=zero end if ! then alltoallv and i think we are done?? call mpi_alltoallv(tempa,sendcounts,sdispls,mpi_rtype, & all_loc,recvcounts,rdispls,mpi_rtype,mpi_comm_world,ierror) end subroutine generic_grid2sub8 subroutine reorder2_s8(work,k_in) ! !USES: use constants, only: zero use mpimod, only: npe use gridmod, only: ijn_s,itotsub use kinds, only: r_kind,i_kind implicit none ! !INPUT PARAMETERS: integer(i_kind), intent(in) :: k_in ! number of levs in work array ! !INPUT/OUTPUT PARAMETERS: real(r_kind),dimension(itotsub,k_in),intent(inout):: work ! !OUTPUT PARAMETERS: ! !DESCRIPTION: adapt reorder2 to single precision ! ! !REVISION HISTORY: ! ! 2004-01-25 kleist ! 2004-05-14 kleist, documentation ! 2004-07-15 todling, protex-complaint prologue ! 2004-11-29 parrish, adapt reorder2 to single precision ! ! !REMARKS: ! language: f90 ! machine: ibm rs/6000 sp; sgi origin 2000; compaq/hp ! ! !AUTHOR: ! kleist org: np20 date: 2004-01-25 ! !EOP !------------------------------------------------------------------------- integer(i_kind) iloc,iskip,i,k,n real(r_kind),dimension(itotsub*k_in):: temp ! Zero out temp array do k=1,itotsub*k_in temp(k)=zero end do ! Load temp array in order of subdomains iloc=0 iskip=0 do n=1,npe if (n/=1) then iskip=iskip+ijn_s(n-1) end if do k=1,k_in do i=1,ijn_s(n) iloc=iloc+1 temp(iloc)=work(iskip+i,k) end do end do end do ! Now load the tmp array back into work iloc=0 do k=1,k_in do i=1,itotsub iloc=iloc+1 work(i,k)=temp(iloc) end do end do return end subroutine reorder2_s8 subroutine get_div_reg(u,v,div) use kinds, only: r_kind,i_kind use gridmod, only: nlat,nlon implicit none real(r_kind),dimension(nlat,nlon),intent(in):: u,v real(r_kind),dimension(nlat,nlon),intent(out)::div integer(i_kind) i,j real(r_kind),dimension(nlat,nlon)::ux,vy call delx_reg(u,ux,.true.) call dely_reg(v,vy,.true.) do j=1,nlon do i=1,nlat div(i,j)=ux(i,j)+vy(i,j) end do end do end subroutine get_div_reg subroutine get_div_reg_ad(u,v,div) use kinds, only: r_kind,i_kind use gridmod, only: nlat,nlon use constants, only: zero implicit none real(r_kind),dimension(nlat,nlon),intent(out):: u,v real(r_kind),dimension(nlat,nlon),intent(in)::div integer(i_kind) i,j real(r_kind),dimension(nlat,nlon)::ux,vy do j=1,nlon do i=1,nlat ux(i,j)=div(i,j) vy(i,j)=div(i,j) end do end do u=zero ; v=zero call tdelx_reg(ux,u,.true.) call tdely_reg(vy,v,.true.) end subroutine get_div_reg_ad subroutine get_vor_reg(u,v,vor) use kinds, only: r_kind,i_kind use gridmod, only: nlat,nlon implicit none real(r_kind),dimension(nlat,nlon),intent(in):: u,v real(r_kind),dimension(nlat,nlon),intent(out)::vor integer(i_kind) i,j real(r_kind),dimension(nlat,nlon)::uy,vx call delx_reg(v,vx,.true.) call dely_reg(u,uy,.true.) do j=1,nlon do i=1,nlat vor(i,j)=vx(i,j)-uy(i,j) end do end do end subroutine get_vor_reg subroutine get_vor_reg_ad(u,v,vor) use kinds, only: r_kind,i_kind use gridmod, only: nlat,nlon use constants, only: zero implicit none real(r_kind),dimension(nlat,nlon),intent(out):: u,v real(r_kind),dimension(nlat,nlon),intent(in)::vor integer(i_kind) i,j real(r_kind),dimension(nlat,nlon)::uy,vx do j=1,nlon do i=1,nlat vx(i,j)=vor(i,j) uy(i,j)=-vor(i,j) end do end do u=zero ; v=zero call tdelx_reg(vx,v,.true.) call tdely_reg(uy,u,.true.) end subroutine get_vor_reg_ad