module projmethod_support !$$$ module documentation block ! . . . . ! module: projmethod_support ! ! abstract: ! ! program history log: ! 2007-10-01 pondeca ! ! subroutines included: ! sub init_mgram_schmidt ! sub destroy_mgram_schmidt ! sub mgram_schmidt ! ! functions included: ! dplev_mask ! fast_dplev ! dplev5 ! ! variable definition: ! ! attributes: ! language: f90 ! machine: ! !$$$ end documentation block ! Uses: use kinds, only: r_kind use control_vectors implicit none PRIVATE PUBLIC init_mgram_schmidt, mgram_schmidt, & destroy_mgram_schmidt ! Declare local variables real(r_kind),allocatable,dimension(:,:)::gx,gy contains subroutine init_mgram_schmidt !$$$ subprogram documentation block ! . . . . ! subprogram: init_mgram_schmidt ! prgmmr: ! ! abstract: ! ! program history log: ! 2009-09-22 lueken - added subprogram doc block ! ! input argument list: ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ! !$$$ end documentation block use jfunc, only: nclen,jiter,niter implicit none ! Declare passed variables ! Declare local variables allocate(gx(nclen,0:niter(jiter))) allocate(gy(nclen,0:niter(jiter))) end subroutine init_mgram_schmidt subroutine destroy_mgram_schmidt !$$$ subprogram documentation block ! . . . . ! subprogram: destroy_mgram_schmidt ! prgmmr: ! ! abstract: ! ! program history log: ! 2009-09-22 lueken - added subprogram doc block ! ! input argument list: ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ! !$$$ end documentation block implicit none ! Declare passed variables deallocate(gx) deallocate(gy) end subroutine destroy_mgram_schmidt subroutine mgram_schmidt(gradx,grady) !$$$ subprogram documentation block ! . . . . ! subprogram: mgram_schmidt ! prgmmr: pondeca org: np22 date: 2007-08-01 ! ! abstract: apply modified gram-schmidt orthoginalization procedure ! to set of bi-orthogonal vectors ! ! program history log: ! 2007-08-01 pondeca ! ! input argument list: ! gradx - gradient of the cost function w.r.t control variable ! grady - B*(gradx) where B is background error covariance matrix ! ! ! output argument list: ! gradx - modified gradient of the cost function ! grady - modified B*(gradx) ! ! ! remarks: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: i_kind use jfunc, only: iter,nclen use constants, only: izero,ione,tiny_r_kind use mpimod, only: mype implicit none ! Declare passed variables type(control_vector),intent(inout) :: gradx,grady ! Declare local variables integer(i_kind) i,k real(r_kind) prd0 !********************************************************************** gx(1:nclen,iter)=gradx%values(1:nclen) gy(1:nclen,iter)=grady%values(1:nclen) !==> orthogonalization + renormalization do k=1,2 do i=izero,iter-ione prd0=dplev_mask(gy(1,iter),gx(1,i),mype) gx(1:nclen,iter)=gx(1:nclen,iter)-gx(1:nclen,i)*prd0 gy(1:nclen,iter)=gy(1:nclen,iter)-gy(1:nclen,i)*prd0 enddo prd0=dplev_mask(gx(1,iter),gy(1,iter),mype) if (prd0 <= tiny_r_kind) then if (mype==izero) then print*,'in mgram_schmidt: unable to bi-orthogonalize due to round-off error for iter,k=',iter,k print*,'in mgram_schmidt: likely to happen when using fast version of inner product' print*,'in mgram_schmidt: iter,k,prd0=',iter,k,prd0 endif goto 100 endif gx(1:nclen,iter)=gx(1:nclen,iter)/sqrt(prd0) gy(1:nclen,iter)=gy(1:nclen,iter)/sqrt(prd0) enddo !==> update gradx and grady and put correct B-norm back prd0=dplev_mask(gradx%values,grady%values,mype) gradx%values(1:nclen)=gx(1:nclen,iter)*sqrt(prd0) grady%values(1:nclen)=gy(1:nclen,iter)*sqrt(prd0) 100 continue contains real(r_kind) function dplev_mask(dx,dy,mype) !$$$ subprogram documentation block ! . . . . ! subprogram: dplev_mask ! prgmmr: ! ! abstract: ! ! program history log: ! 2009-09-22 lueken - added subprogram doc block ! ! input argument list: ! mype ! dx,dy ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ! !$$$ end documentation block use jfunc, only: nval_levs use gridmod, only: lat2,lon2,twodvar_regional implicit none ! Declar passed variables real(r_kind),dimension(lat2,lon2,nval_levs),intent(in ) :: dx,dy integer(i_kind) ,intent(in ) :: mype ! Declare local variables logical mask(nval_levs) logical fast fast=.false. mask=.true. ! set fast to .true. for twodvar_regional, ! substantially faster, but no roundoff error reduction and ! results differ for different number of processors. if(twodvar_regional) then ! fast=.true. mask(5)=.false. mask(6)=.false. mask(8)=.false. end if if(fast) then dplev_mask=fast_dplev(dx,dy,mask) else dplev_mask=dplev5(dx,dy,mype,mask) end if end function dplev_mask real(r_kind) function fast_dplev(dx,dy,mask) !$$$ subprogram documentation block ! . . . . ! subprogram: fast_dplev ! prgmmr: ! ! abstract: ! ! program history log: ! 2009-09-22 lueken - added subprogram doc block ! ! input argument list: ! dx,dy ! mask ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ! !$$$ end documentation block use jfunc, only: nval_levs use gridmod, only: lat2,lon2 use mpimod, only: npe,mpi_comm_world,ierror,mpi_rtype use constants, only: zero implicit none ! Declar passed variables real(r_kind),dimension(lat2,lon2,nval_levs),intent(in ) :: dx,dy logical ,intent(in ) :: mask(nval_levs) ! Declare local variables real(r_kind),dimension(npe):: sumall integer(i_kind) i,j,k real(r_kind) sum sum=zero do k=1,nval_levs if(.not.mask(k)) cycle do j=2,lon2-ione do i=2,lat2-ione sum=sum+dx(i,j,k)*dy(i,j,k) end do end do end do call mpi_allgather(sum,ione,mpi_rtype,sumall,ione,mpi_rtype,mpi_comm_world,ierror) fast_dplev=zero do i=1,npe fast_dplev=fast_dplev+sumall(i) end do end function fast_dplev real(r_kind) function dplev5(dx,dy,mype,mask) !$$$ subprogram documentation block ! . . . . ! subprogram: dplev calculates dot product for data on subdomain ! prgmmr: derber org: np23 date: 2004-05-13 ! ! abstract: calculates dot product for data on subdomain. Note loops over ! streamfunction, velocity potential, temperature, etc. Also, only ! over interior of subdomain. ! ! program history log: ! 2004-05-13 derber, document ! 2004-07-28 treadon - add only on use declarations; add intent in/out ! 2010-04-01 treadon - move strip to grimod ! ! input argument list: ! dx - input vector 1 ! dy - input vector 2 ! mype ! mask ! ! output argument list ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use jfunc, only: nval_levs use gridmod, only: nlat,nlon,lat2,lon2,lat1,lon1,& ltosi,ltosj,iglobal,ijn,displs_g,strip use mpimod, only: mpi_comm_world,ierror,mpi_rtype use constants, only: zero implicit none ! Declare passed variables real(r_kind),dimension(lat2,lon2,nval_levs),intent(in ) :: dx,dy integer(i_kind) ,intent(in ) :: mype logical ,intent(in ) :: mask(nval_levs) ! Declare local variables real(r_kind),dimension(lat1*lon1):: zsm real(r_kind),dimension(iglobal):: work1 real(r_kind),dimension(lat2,lon2):: sum real(r_kind),dimension(nlat,nlon):: sumall integer(i_kind) i,j,k,mm1 real(r_kind) e,y,temp mm1=mype+ione sum=zero do k=1,nval_levs if(.not.mask(k)) cycle do j=1,lon2 do i=1,lat2 sum(i,j)=sum(i,j)+dx(i,j,k)*dy(i,j,k) end do end do end do do j=1,lon1*lat1 zsm(j)=zero end do call strip(sum,zsm,ione) call mpi_allgatherv(zsm,ijn(mm1),mpi_rtype,& work1,ijn,displs_g,mpi_rtype,& mpi_comm_world,ierror) do k=1,iglobal i=ltosi(k) ; j=ltosj(k) sumall(i,j)=work1(k) end do dplev5=zero e=zero do j=1,nlon do i=1,nlat ! Compensated summation version of sum temp=dplev5 y=sumall(i,j)+e dplev5=temp+y e=(temp-dplev5)+y ! dplev=dplev+sumall(i,j) end do end do end function dplev5 end subroutine mgram_schmidt end module projmethod_support subroutine writeout_gradients(dx,dy,nv,alpha,gamma,mype) !$$$ subprogram documentation block ! . . . . ! subprogram: writeout_gradients ! prgmmr: pondeca org: np23 date: 2006-10-17 ! ! abstract: writes out the x and y gradidents of the costfunction ! also writes out B*y ! ! program history log: ! 2006-10-17 pondeca, document ! 2010-03-29 zhu - make changes for generalizing control variables ! 2010-04-01 treadon - move strip to grimod ! ! input argument list: ! nv ! mype ! dx - input vector 1 ! dy - input vector 2 ! alpha ! gamma ! ! output argument list ! !attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ !************************************************************************* use kinds, only: r_kind,i_kind use mpimod, only: mpi_comm_world,ierror,mpi_rtype use gridmod, only: nsig,lat1,lon1,lat2,lon2,nlon,nlat,itotsub,iglobal, & ltosi,ltosj,latlon11,ijn,displs_g,strip use radinfo, only: npred,jpch_rad use pcpinfo, only: npredp,npcptype use jfunc, only: iter,jiter use control_vectors use constants, only: izero,ione implicit none ! Declare passed variables integer(i_kind) ,intent(in ) :: nv,mype type(control_vector),intent(in ) :: dx,dy real(r_kind) ,intent(in ) :: alpha,gamma ! Declare local variables integer(i_kind) i,k,k1,k2,lun,ifield,icase,ii real(r_kind),allocatable,dimension(:)::tempa real(r_kind),allocatable,dimension(:,:)::slab real(r_kind),allocatable,dimension(:)::strp real(r_kind),allocatable,dimension(:)::field type(control_vector)::dz character(2) clun1 character(3) clun2 !************************************************************************* call allocate_cv(dz) allocate(tempa(itotsub)) allocate(slab(nlon,nlat)) allocate(strp(lat1*lon1)) allocate(field(lat2*lon2*nsig)) write (clun1(1:2),'(i2.2)') jiter write (clun2(1:3),'(i3.3)') iter lun=19_i_kind do icase=1,2 if (icase==ione) then dz=dx open (lun,file='gradx.dat_'//clun1//'_'//clun2,form='unformatted') else if (icase==2_i_kind) then dz=dy open (lun,file='grady.dat_'//clun1//'_'//clun2,form='unformatted') endif write (lun) nlon,nlat,nsig,jpch_rad,npred,npcptype,npredp,jiter,nv,alpha,gamma ii=ione do ifield=1,nrf3 if (ifield==nrf3_sf) field=dz%step(ii)%st if (ifield==nrf3_vp) field=dz%step(ii)%vp if (ifield==nrf3_t) field=dz%step(ii)%t if (ifield==nrf3_q) field=dz%step(ii)%rh if (ifield==nrf3_oz) field=dz%step(ii)%oz if (ifield==nrf3_cw) field=dz%step(ii)%cw do k=1,nsig k1=ione+(k-ione)*latlon11 k2=k1+latlon11-ione call strip(field(k1:k2),strp,ione) call mpi_gatherv(strp,ijn(mype+ione),mpi_rtype, & tempa,ijn,displs_g,mpi_rtype,izero,mpi_comm_world,ierror) if(mype == izero) then do i=1,iglobal slab(ltosj(i),ltosi(i))=tempa(i) end do write(lun) slab endif end do enddo !ifield ! gradient wrt sfcp call strip(dz%step(ii)%p,strp,ione) call mpi_gatherv(strp,ijn(mype+ione),mpi_rtype, & tempa,ijn,displs_g,mpi_rtype,izero,mpi_comm_world,ierror) if(mype == izero) then do i=1,iglobal slab(ltosj(i),ltosi(i))=tempa(i) end do write(lun) slab endif ! gradient wrt sfct if (nrf2_sst>izero) then call strip(dz%step(ii)%sst,strp,ione) call mpi_gatherv(strp,ijn(mype+ione),mpi_rtype, & tempa,ijn,displs_g,mpi_rtype,izero,mpi_comm_world,ierror) if(mype == izero) then do i=1,iglobal slab(ltosj(i),ltosi(i))=tempa(i) end do write(lun) slab endif endif ! gradient wrt satellite radiance bias correction coefficients if (mype==izero) write(lun) dz%predr ! gradient wrt precipitation bias correction coefficients if (mype==izero)write(lun) dz%predp close(lun) end do ! icase call deallocate_cv(dz) deallocate(tempa) deallocate(slab) deallocate(strp) deallocate(field) return end subroutine writeout_gradients