module gsi_bias !$$$ module documentation block ! . . . . ! module: gsi_bias ! prgmmr: derber org: np23 date: 2017-06-07 ! ! abstract: This module contains routines which handle bias ! operations for GSI atmospheric and surface files. ! ! program history log: ! 2017-06-07 derber - copy from gsi_io to separate out bias functions ! ! Subroutines Included: ! sub read_bias - read gsi guess bias file from binary file, scatter ! from full grid to subdomains ! sub write_bias - gather gsi guess bias from subdomains to full ! grid, write to binary file ! sub reorder21s_ - ! sub reorder21d_ - ! sub reorder12s_ - ! sub reorder12d_ - ! ! Variable Definitions: ! ! attributes: ! language: f90 ! machine: ! !$$$ end documentation block use kinds, only: i_kind implicit none private public read_bias public write_bias public reorder21 public reorder12 interface reorder21; module procedure & reorder21s_, & reorder21d_ end interface interface reorder12; module procedure & reorder12s_, & reorder12d_ end interface character(len=*), parameter :: myname='gsi_bias' contains subroutine read_bias(filename,mype,nbc,sub_z,bundle,istatus) !$$$ subprogram documentation block ! . . . . ! subprogram: read_bias read bias, convert to grid and ! send to all mpi tasks ! prgmmr: treadon org: np23 date: 2006-04-15 ! ! abstract: read bias, convert to grid, and ! scatter to subdomains ! ! program history log: ! 2006-04-15 treadon ! 2006-12-04 todling - add nbc and loop over nbc ! 2007-06-01 todling - bug fix: loops were only copying to (1,1) element ! 2014-10-05 todling - bias estimates now kept in bundle ! ! input argument list: ! filename - name of local file from which to read bias ! mype - mpi task id ! ! output argument list: ! sub_z - terrain ! bundle - bundle with background bias estimates ! istatus - read status indicator ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_kind,r_single use gridmod, only: itotsub,nlon,nlat,lat2,lon2,nsig,displs_s,ijn_s use constants, only: zero use mpimod, only: mpi_rtype,ierror,mpi_comm_world use gsi_bundlemod, only: gsi_bundle use gsi_bundlemod, only: gsi_bundlegetpointer use m_gsiBiases, only: bvars2d,bvars3d use mpeu_util, only: die implicit none ! Declare local parameters integer(i_kind):: lunin=11 integer(i_kind):: nsize=4 ! Declare passed variables character(24) ,intent(in ) :: filename integer(i_kind) ,intent(in ) :: mype integer(i_kind) ,intent(in ) :: nbc integer(i_kind) ,intent( out) :: istatus type(gsi_bundle) ,intent(inout) :: bundle(nbc) real(r_kind),dimension(lat2,lon2,nbc) ,intent( out) :: sub_z ! Declare local variables integer(i_kind) i,j,k,mm1,nv integer(i_kind) mype_in,iret integer(i_kind):: ib,nb,ka,n real(r_kind),dimension(itotsub):: work real(r_single),dimension(nlon,nlat):: grid4 real(r_kind),dimension(lat2,lon2,nsig)::work3d real(r_kind),pointer,dimension(:,:) :: ptr2d=>NULL() real(r_kind),pointer,dimension(:,:,:):: ptr3d=>NULL() !****************************************************************************** ! Initialize variables used below mype_in=0 mm1=mype+1 ib=-1 nb=nsize*nlon*nlat ! Open file to read bias fields istatus=0 call baopenr(lunin,filename,iret) if (iret/=0) then if (mype==mype_in) write(6,*) & 'READ_BIAS: ***ERROR*** opening output file, iret=',iret,lunin,filename istatus=istatus+iret return endif ! Loop over all coefficients of bias model do n=1,nbc ! Terrain: spectral --> grid transform, scatter to all mpi tasks if (mype==mype_in) then call baread(lunin,ib,nb,ka,grid4) call reorder21(grid4,work) endif call mpi_scatterv(work,ijn_s,displs_s,mpi_rtype,& sub_z(1,1,n),ijn_s(mm1),mpi_rtype,mype_in,mpi_comm_world,ierror) ! 2d fields: do nv=1,size(bvars2d) if (mype==mype_in) then call baread(lunin,ib,nb,ka,grid4) call reorder21(grid4,work) endif call gsi_bundlegetpointer(bundle(n),bvars2d(nv),ptr2d,ierror) if(ierror/=0) call die('trouble in reading 2d bias estimates') call mpi_scatterv(work,ijn_s,displs_s,mpi_rtype,& ptr2d,ijn_s(mm1),mpi_rtype,mype_in,mpi_comm_world,ierror) enddo ! 3d fields: do nv=1,size(bvars3d) call gsi_bundlegetpointer(bundle(n),bvars3d(nv),ptr3d,ierror) if(ierror/=0) call die('trouble in reading 3d bias estimates') do k=1,nsig if (mype==mype_in) then call baread(lunin,ib,nb,ka,grid4) call reorder21(grid4,work) endif call mpi_scatterv(work,ijn_s,displs_s,mpi_rtype,& work3d(1,1,k),ijn_s(mm1),mpi_rtype,mype_in,mpi_comm_world,ierror) do j=1,lon2 do i=1,lat2 ptr3d(i,j,k) = work3d(i,j,k) end do end do end do end do end do ! End loop over coefficients ! Close input file call baclose(lunin,iret) if (iret==0) then write(6,*)'READ_BIAS: read in previous estimate of bkg bias' else write(6,*)'READ_BIAS: ***ERROR*** closing input file, iret=',iret endif istatus=istatus+iret ! End of routine. Return return end subroutine read_bias !------------------------------------------------------------------------- ! NOAA/NCEP, National Centers for Environmental Prediction GSI ! !------------------------------------------------------------------------- !BOP ! ! !ROUTINE: write_bias --- Gather, transform, and write out spectal coefficients ! ! !INTERFACE: ! subroutine write_bias(filename,mype_out,nbc,sub_z,bundle,istatus) ! ! !USES: ! use kinds, only: r_kind,r_single use mpimod, only: mpi_rtype use mpimod, only: mpi_comm_world use mpimod, only: ierror use mpimod, only: mype use gridmod, only: nlat, nlon ! no. lat/lon use gridmod, only: lat1, lon1 ! no. lat/lon on subdomain (no buffer) use gridmod, only: lat2, lon2 ! no. lat/lon on subdomain (buffer pnts on ends) use gridmod, only: nsig ! no. levels use gridmod, only: iglobal ! no. of horizontal points on global grid use gridmod, only: ijn ! no. of horiz. pnts for each subdomain (no buffer) use gridmod, only: displs_g ! comm. array, displacement for receive on global grid use gridmod, only: itotsub ! no. of horizontal points of all subdomains combined use gridmod, only: strip use gsi_bundlemod, only: gsi_bundle use gsi_bundlemod, only: gsi_bundlegetpointer use gsi_bundlemod, only: gsi_bundlecreate use gsi_bundlemod, only: gsi_bundledestroy use m_gsiBiases, only: bvars2d,bvars3d use m_gsiBiases, only: bkg_bias_model use gsi_4dcouplermod, only: gsi_4dcoupler_putpert use obsmod, only: iadate use mpeu_util, only: die implicit none ! ! !LOCAL PARAMETER: ! ! ! !INPUT PARAMETERS: ! character(24) ,intent(in):: filename ! file to open and write to integer(i_kind) ,intent(in ) :: mype_out ! mpi task to write output file integer(i_kind) ,intent(in ) :: nbc ! number of bias coefficients in bias model integer(i_kind) ,intent( out) :: istatus ! write status real(r_kind),dimension(lat2,lon2,nbc),intent(in ) :: sub_z type(gsi_bundle),intent(inout) :: bundle(nbc) ! holds all bkg bias estimates ! ! !OUTPUT PARAMETERS: ! ! !DESCRIPTION: This routine gathers fields needed for the GSI analysis ! file from subdomains and then transforms the fields from ! grid to spectral space. The spectral coefficients are ! then written to an atmospheric analysis file. ! ! !REVISION HISTORY: ! ! 2006-12-04 todling - add nbc and loop over nbc ! 2010-04-01 treadon - move strip to gridmod ! 2013-10-24 todling - revisit strip interface ! 2014-10-05 todling - bias estimates now kept in bundle ! - rename bkg-bias interface prog for clarity ! ! !REMARKS: ! ! language: f90 ! machines: ibm RS/6000 SP; SGI Origin 2000; Compaq HP ! ! !AUTHOR: ! ! 2006-04-15 treadon ! !EOP !------------------------------------------------------------------------- character(len=*), parameter :: myname_=myname//'*write_bias_' integer(i_kind),parameter:: lunout = 51 integer(i_kind),parameter:: nsize=4 integer(i_kind) k,mm1 integer(i_kind):: iret integer(i_kind):: nb,n,nv integer(i_kind):: nymd,nhms real(r_kind),dimension(lat1*lon1):: zsm,work2dm real(r_kind),dimension(lat1*lon1,nsig):: work3dm real(r_kind),dimension(max(iglobal,itotsub)):: work real(r_single),dimension(nlon,nlat):: grid4 real(r_kind),pointer,dimension(:,:) :: ptr2d=>NULL() real(r_kind),pointer,dimension(:,:,:):: ptr3d=>NULL() type(gsi_bundle) xbundle !************************************************************************* ! Initialize local variables mm1=mype+1 nb=nsize*nlon*nlat ! Open file to receive bias fields istatus=0 if (mype==mype_out) then call baopenwt(lunout,filename,iret) if (iret/=0) then write(6,*)'WRITE_BIAS: ***ERROR*** opening output file, iret=',iret endif istatus=istatus+iret endif ! Loop over number of coefficients in bias model do n=1,nbc ! Terrain call strip(sub_z(:,:,n),zsm) call mpi_gatherv(zsm,ijn(mm1),mpi_rtype,& work,ijn,displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) if (mype==mype_out) then call reorder12(work,grid4) call wryte(lunout,nb,grid4) endif ! 2d fields: do nv=1,size(bvars2d) call gsi_bundlegetpointer(bundle(n),bvars2d(nv),ptr2d,ierror) if(ierror/=0) call die('trouble in writing 2d bias estimates') ! Strip off boundary points from subdomains call strip(ptr2d,work2dm) ! Create global grid by gathering from subdomains call mpi_gatherv(work2dm,ijn(mm1),mpi_rtype,& work,ijn,displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) ! Write full grid field to output file if (mype==mype_out) then call reorder12(work,grid4) call wryte(lunout,nb,grid4) endif enddo ! 3d fields: do nv=1,size(bvars3d) call gsi_bundlegetpointer(bundle(n),bvars3d(nv),ptr3d,ierror) if(ierror/=0) call die('trouble in writing 3d bias estimates') ! Strip off boundary points from subdomains call strip(ptr3d,work3dm,nsig) ! For each level ... do k=1,nsig ! Create global grid by gathering from subdomains call mpi_gatherv(work3dm(1,k),ijn(mm1),mpi_rtype,& work,ijn,displs_g,mpi_rtype,& mype_out,mpi_comm_world,ierror) ! Write slice of 3d field to output file if (mype==mype_out) then call reorder12(work,grid4) call wryte(lunout,nb,grid4) endif end do enddo end do ! End loop over nbc ! Single task writes message to stdout if (mype==mype_out) then write(6,*) 'WRITE_BIAS: bias file written to ',& trim(filename) call baclose(lunout,iret) if (iret/=0) then write(6,*)'WRITE_BIAS: ***ERROR*** closing output file, iret=',iret endif istatus=istatus+iret endif nymd = 10000*iadate(1)+iadate(2)*100+iadate(3) nhms = 10000*iadate(4) if(mype==0) write(6,'(2a,i8.8,2x,i6.6)')trim(myname_),': writing out bias on ',& nymd, nhms call gsi_bundlecreate(xbundle,bundle(1),'Bias Estimate',iret) call bkg_bias_model(xbundle,iadate(4)) call gsi_4dcoupler_putpert (xbundle,nymd,nhms,'tlm','bbias') call gsi_bundledestroy(xbundle,iret) ! return end subroutine write_bias !------------------------------------------------------------------------- ! NOAA/NCEP, National Centers for Environmental Prediction GSI ! !------------------------------------------------------------------------- !BOP ! ! !IROUTINE: reorder21s_ --- reorder 2d array to 1d order ! ! !INTERFACE: ! subroutine reorder21s_(grid_in,grid_out) ! !USES: use kinds, only: r_kind,r_single use gridmod, only: itotsub,nlat,nlon use general_commvars_mod, only: ltosi_s,ltosj_s implicit none ! !INPUT PARAMETERS: real(r_single),dimension(nlon,nlat),intent(in ) :: grid_in ! input grid real(r_kind) ,dimension(itotsub) ,intent( out) :: grid_out ! output grid ! !DESCRIPTION: This routine transfers the contents of a two-diemnsional, ! type r_single array into a one-dimension, type r_kind ! array. ! ! ! !REVISION HISTORY: ! 2004-08-27 treadon ! 2013-10-25 todling - repositioned ltosi and others to commvars ! ! !REMARKS: ! language: f90 ! machine: ibm rs/6000 ! ! !AUTHOR: ! treadon org: np23 date: 2004-08-27 ! !EOP !------------------------------------------------------------------------- ! Declare local variables integer(i_kind) i,j,k ! Transfer input 2d array to output 1d array do k=1,itotsub i=ltosi_s(k) j=ltosj_s(k) grid_out(k)=grid_in(j,i) end do return end subroutine reorder21s_ !------------------------------------------------------------------------- ! NOAA/NCEP, National Centers for Environmental Prediction GSI ! !------------------------------------------------------------------------- !BOP ! ! !IROUTINE: reorder21d_ --- reorder 2d array to 1d order ! ! !INTERFACE: ! subroutine reorder21d_(grid_in,grid_out) ! !USES: use kinds, only: r_kind,r_double use gridmod, only: itotsub,nlat,nlon use general_commvars_mod, only: ltosi_s,ltosj_s implicit none ! !INPUT PARAMETERS: real(r_double),dimension(nlon,nlat),intent(in ) :: grid_in ! input grid real(r_kind),dimension(itotsub) ,intent( out) :: grid_out ! output grid ! !DESCRIPTION: This routine transfers the contents of a two-diemnsional, ! type r_single array into a one-dimension, type r_kind ! array. ! ! ! !REVISION HISTORY: ! 2004-08-27 treadon ! 2007-05-27 todling - add double precision version ! 2011-07-03 todling - true double prec interface ! 2013-10-25 todling - repositioned ltosi and others to commvars ! ! !REMARKS: ! language: f90 ! machine: ibm rs/6000 ! ! !AUTHOR: ! treadon org: np23 date: 2004-08-27 ! !EOP !------------------------------------------------------------------------- ! Declare local variables integer(i_kind) i,j,k ! Transfer input 2d array to output 1d array do k=1,itotsub i=ltosi_s(k) j=ltosj_s(k) grid_out(k)=grid_in(j,i) end do return end subroutine reorder21d_ !------------------------------------------------------------------------- ! NOAA/NCEP, National Centers for Environmental Prediction GSI ! !------------------------------------------------------------------------- !BOP ! ! !IROUTINE: reorder12s_ --- reorder 1d array to 2d order ! ! !INTERFACE: ! subroutine reorder12s_(grid_in,grid_out) ! !USES: use kinds, only: r_kind,r_single use gridmod, only: itotsub,iglobal,nlat,nlon use general_commvars_mod, only: ltosi,ltosj implicit none ! !INPUT PARAMETERS: real(r_kind) ,dimension(max(iglobal,itotsub)),intent(in ) :: grid_in ! input grid real(r_single),dimension(nlon,nlat) ,intent( out) :: grid_out ! input grid ! !DESCRIPTION: This routine transfers the contents of a one-diemnsional, ! type r_kind array into a two-dimensional, type r_single ! array. ! ! ! !REVISION HISTORY: ! 2004-08-27 treadon ! 2013-10-25 todling - repositioned ltosi and others to commvars ! ! !REMARKS: ! language: f90 ! machine: ibm rs/6000 ! ! !AUTHOR: ! treadon org: np23 date: 2004-08-27 ! !EOP !------------------------------------------------------------------------- ! Declare local variables integer(i_kind) i,j,k ! Transfer input 1d array to output 2d array do k=1,iglobal i=ltosi(k) j=ltosj(k) grid_out(j,i) = grid_in(k) end do return end subroutine reorder12s_ !------------------------------------------------------------------------- ! NOAA/NCEP, National Centers for Environmental Prediction GSI ! !------------------------------------------------------------------------- !BOP ! ! !IROUTINE: reorder12d_ --- reorder 1d array to 2d order ! ! !INTERFACE: ! subroutine reorder12d_(grid_in,grid_out) ! !USES: use kinds, only: r_kind,r_double use gridmod, only: itotsub,iglobal,nlat,nlon use general_commvars_mod, only: ltosi,ltosj implicit none ! !INPUT PARAMETERS: real(r_kind),dimension(max(iglobal,itotsub)),intent(in ) :: grid_in ! input grid real(r_double),dimension(nlon,nlat) ,intent( out) :: grid_out ! input grid ! !DESCRIPTION: This routine transfers the contents of a one-diemnsional, ! type r_kind array into a two-dimensional, type r_single ! array. ! ! ! !REVISION HISTORY: ! 2004-08-27 treadon ! 2007-05-27 todling - add double precision version ! 2011-07-03 todling - true double prec interface ! 2013-10-25 todling - repositioned ltosi and others to commvars ! ! !REMARKS: ! language: f90 ! machine: ibm rs/6000 ! ! !AUTHOR: ! treadon org: np23 date: 2004-08-27 ! !EOP !------------------------------------------------------------------------- ! Declare local variables integer(i_kind) i,j,k ! Transfer input 1d array to output 2d array do k=1,iglobal i=ltosi(k) j=ltosj(k) grid_out(j,i) = grid_in(k) end do return end subroutine reorder12d_ end module gsi_bias