subroutine general_write_gfsatm(grd,sp_a,sp_b,filename,mype_out,&
gfs_bundle,ibin,inithead,iret_write)
!$$$ subprogram documentation block
! . . . .
! subprogram: general_write_gfsatm adaptation of write_gfsatm for general resolutions
! prgmmr: parrish org: np22 date: 1990-10-10
!
! abstract: copied from write_gfsatm, primarily for writing in gefs sigma files, where the
! input resolution and the grid that variables are reconstructed on can be
! different from the analysis grid/resolution.
!
! program history log:
! 2010-02-25 parrish
! 2010-03-29 todling - add prologue; load_grid now in commvars
! 2014-12-03 derber - simplify if structure and use guess surface height
! directly
! 2016-05-06 thomas - recalculate cw increment to account for qcmin
!
! input argument list:
!
! inithead - logical to read header record. Usually .true. unless
! repeatedly reading similar files(e.g., ensembles)
! output argument list:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind,r_single
use sigio_r_module, only: sigio_dbti,sigio_rropen,sigio_rrhead,sigio_rwhead,&
sigio_rrdbti,sigio_rwdbti,sigio_rwopen,sigio_rclose,sigio_aldbti
use sigio_module, only: sigio_head,sigio_alhead
use general_sub2grid_mod, only: sub2grid_info
use guess_grids, only: ifilesig
use obsmod, only: iadate
use mpimod, only: npe,mype
use general_specmod, only: spec_vars
use gridmod, only: ntracer,ncepgfs_head,idpsfc5,idthrm5,cp5,idvc5,idvm5
use general_commvars_mod, only: load_grid
use ncepgfs_io, only: sigio_cnvtdv8,sighead
use constants, only: zero,zero_single,one,fv,qcmin
use gsi_4dvar, only: ibdate,nhr_obsbin,lwrite4danl
use gsi_bundlemod, only: gsi_bundle
use gsi_bundlemod, only: gsi_bundlegetpointer
implicit none
! INPUT PARAMETERS:
character(*), intent(in ) :: filename ! file to open and write to
integer(i_kind), intent(in ) :: mype_out ! mpi task number
type(sub2grid_info), intent(in ) :: grd
type(spec_vars), intent(in ) :: sp_a,sp_b
type(gsi_bundle), intent(in ) :: gfs_bundle
logical, intent(in ) :: inithead
integer(i_kind), intent(in ) :: ibin
integer(i_kind), intent( out) :: iret_write
! LOCAL VARIABLES
integer(i_kind),parameter:: lunges = 11
integer(i_kind),parameter:: lunanl = 51
character(6):: fname_ges
real(r_kind),pointer,dimension(:,:) :: sub_ps
real(r_kind),pointer,dimension(:,:,:) :: sub_vor,sub_div,sub_tv
real(r_kind),pointer,dimension(:,:,:) :: sub_q,sub_oz,sub_cwmr
real(r_kind),dimension(grd%itotsub):: work
real(r_kind),dimension(grd%nlon,grd%nlat-2):: grid,grid2
real(r_kind),dimension(grd%lat2,grd%lon2):: work_ps
real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig):: work_tv
real(r_kind),dimension(sp_b%nc):: spec_work
real(r_kind),dimension(sp_a%nc):: spec_work_sm
real(r_kind),dimension(sp_b%nc),target :: specges_4
integer(i_kind) :: nlatm2,icount,itotflds,i,j,iret,kvar,klev,k,l,n,ks1,ks2,istatus
integer(i_kind),dimension(npe)::ilev,ivar
integer(i_kind),dimension(5):: mydate
integer(i_kind),dimension(8) :: ida,jda
real(r_kind),dimension(5) :: fha
type(sigio_dbti):: sigdati
logical :: lloop
!*************************************************************************
! Handle case of NCEP SIGIO
! Initialize local variables
iret_write=0
nlatm2=grd%nlat-2
itotflds=6*grd%nsig+2 ! Hardwired for now! vor,div,tv,q,oz,cwmr,ps,z
lloop=.true.
istatus=0
call gsi_bundlegetpointer(gfs_bundle,'ps', sub_ps, iret); istatus=istatus+iret
call gsi_bundlegetpointer(gfs_bundle,'vor',sub_vor, iret); istatus=istatus+iret
call gsi_bundlegetpointer(gfs_bundle,'div',sub_div, iret); istatus=istatus+iret
call gsi_bundlegetpointer(gfs_bundle,'tv', sub_tv, iret); istatus=istatus+iret
call gsi_bundlegetpointer(gfs_bundle,'q', sub_q, iret); istatus=istatus+iret
call gsi_bundlegetpointer(gfs_bundle,'oz', sub_oz, iret); istatus=istatus+iret
call gsi_bundlegetpointer(gfs_bundle,'cw', sub_cwmr,iret); istatus=istatus+iret
if ( istatus /= 0 ) then
if ( mype == 0 ) then
write(6,*) 'general_write_gfsatm: ERROR'
write(6,*) 'Missing some of the required fields'
write(6,*) 'Aborting ... '
endif
call stop2(999)
endif
! Set guess file name
write(fname_ges,100) ifilesig(ibin)
100 format('sigf',i2.2)
! Have all files open ges and read header for now with RanRead
if ( mype < itotflds .or. mype == mype_out ) then
call sigio_rropen(lunges,fname_ges,iret)
if ( inithead ) call sigio_rrhead(lunges,sighead,iret)
! All tasks should also open output file for random write
call sigio_rwopen(lunanl,filename,iret_write)
if ( iret_write /= 0 ) then
write(6,*)'GENERAL_WRITE_GFSATM: ***ERROR*** writing ',&
trim(filename),' mype,iret_write=',mype,iret_write
return
end if
endif
! Load date and write header
if ( mype == mype_out ) then
if ( lwrite4danl ) then
! increment mydate
mydate=ibdate
fha(:)=zero ; ida=0; jda=0
fha(2)=real(nhr_obsbin*(ibin-1)) ! relative time interval in hours
ida(1)=mydate(1) ! year
ida(2)=mydate(2) ! month
ida(3)=mydate(3) ! day
ida(4)=0 ! time zone
ida(5)=mydate(4) ! hour
! Move date-time forward by nhr_assimilation hours
call w3movdat(fha,ida,jda)
mydate(1)=jda(1)
mydate(2)=jda(2)
mydate(3)=jda(3)
mydate(4)=jda(5)
else
mydate=iadate
endif
! Replace header record date with analysis time
sighead%fhour = zero_single
sighead%idate(1) = mydate(4) !hour
sighead%idate(2) = mydate(2) !month
sighead%idate(3) = mydate(3) !day
sighead%idate(4) = mydate(1) !year
! Load grid dimension and other variables used below
! into local header structure
! Write header to analysis file
call sigio_rwhead(lunanl,sighead,iret)
iret_write=iret_write+iret
endif
! Surface pressure.
! NCEP SIGIO has two options for surface pressure. Variable idpsfc5
! indicates the type:
! idpsfc5= 0,1 for ln(psfc)
! idpsfc5= 2 for psfc
work_ps=sub_ps
! If output ln(ps), take log of ps in cb
if ( idpsfc5 /= 2 ) then
do j=1,grd%lon2
do i=1,grd%lat2
if ( work_ps(i,j) <= zero ) &
work_ps(i,j)=one
work_ps(i,j)=log(work_ps(i,j))
enddo
enddo
endif
! Thermodynamic variable
! The GSI analysis variable is virtual temperature (Tv). For SIGIO
! we have three possibilities: Tv, sensible temperature (T), or
! enthalpy (h=CpT). Variable idthrm5 indicates the type
! idthrm5 = 0,1 = virtual temperature (Tv)
! idthrm5 = 2 = sensible (dry) temperature (T)
! idthrm5 = 3 = enthalpy (h=CpT)
work_tv=sub_tv
if ( idthrm5 == 2 .or. idthrm5 == 3 ) then
! Convert virtual temperature to dry temperature
do k=1,grd%nsig
do j=1,grd%lon2
do i=1,grd%lat2
work_tv(i,j,k)=work_tv(i,j,k)/(one+fv*sub_q(i,j,k))
enddo
enddo
enddo
! If output is enthalpy, convert dry temperature to CpT
if ( idthrm5 == 3 ) call sigio_cnvtdv8(grd%lat2*grd%lon2,&
grd%lat2*grd%lon2,grd%nsig,idvc5,idvm5,ntracer,&
iret,work_tv,sub_q,cp5,-1)
endif
! Do loop until total fields have been processed. Stop condition on itotflds
icount=0
gfsfields: do while (lloop)
! First, perform sub2grid for up to npe
call general_gather(grd,work_ps,work_tv,sub_vor,sub_div,sub_q,sub_oz,&
sub_cwmr,icount,ivar,ilev,work)
pe_loop: do k=1,npe ! loop over pe distributed data
kvar=ivar(k)
if ( mype == k-1 .and. kvar > 0 ) then
klev=ilev(k)
if ( kvar == 1 ) then ! hs
sigdati%i = 1
else if ( kvar == 2 ) then ! ps
sigdati%i = 2
else if ( kvar == 3 ) then ! temperature
sigdati%i = 2+klev
else if ( kvar == 4 ) then ! vorticity
sigdati%i = sighead%levs + 2 + (klev-1) * 2 + 2
else if ( kvar == 5 ) then ! divergence
sigdati%i = sighead%levs + 2 + (klev-1) * 2 + 1
else if ( kvar == 6 ) then ! q
sigdati%i = sighead%levs * (2+1) + 2 + klev
else if ( kvar == 7 ) then ! oz
sigdati%i = sighead%levs * (2+2) + 2 + klev
else if ( kvar == 8 ) then ! cw, 3rd tracer
sigdati%i = sighead%levs * (2+3) + 2 + klev
endif
if ( klev > 0 ) then
sigdati%f => specges_4
! Read in full resolution guess spectral coefficients
call sigio_rrdbti(lunges,sighead,sigdati,iret)
do i=1,sp_b%nc
spec_work(i) = specges_4(i)
enddo
! Ensure coefficients that must be zero are zero
do i=1,sp_b%nc
if ( sp_b%factsml(i) ) spec_work(i) = zero
enddo
if ( kvar /= 1 ) then ! if sfc elevation field just write out
! Put current analysis on 2d (full level) grid
call load_grid(work,grid)
! Convert full resolution guess to analysis grid
call general_sptez_s_b(sp_a,sp_b,spec_work,grid2,1)
! Calculation grid increment on analysis grid
if ( kvar == 8 ) then
do i=1,sp_a%imax
do j=1,sp_a%jmax
grid(i,j) = grid(i,j) - max(grid2(i,j),qcmin)
enddo
enddo
else
grid = grid - grid2
endif
! Convert grid increment to spectral space
call general_sptez_s(sp_a,spec_work_sm,grid,-1)
! Add increment in spectral space (possibly lower resolution)
! to guess (taken from sppad)
do l=0,min(sp_b%jcap,sp_a%jcap)
do n=l,min(sp_b%jcap,sp_a%jcap)
ks2=l*(2*sp_b%jcap+1-l)+2*n
ks1=l*(2*sp_a%jcap+1-l)+2*n
specges_4(ks2+1)=specges_4(ks2+1)+spec_work_sm(ks1+1)
specges_4(ks2+2)=specges_4(ks2+2)+spec_work_sm(ks1+2)
enddo
enddo
if ( kvar /= 4 .and. kvar /= 5 ) then
do i=1,sp_b%nc
if ( sp_b%factsml(i) ) specges_4(i) = zero_single
enddo
else
do i=1,sp_b%nc
if ( sp_b%factvml(i) ) specges_4(i) = zero_single
enddo
endif
endif ! if ( kvar /= 1 )
! Write out using RanWrite
call sigio_rwdbti(lunanl,sighead,sigdati,iret)
iret_write=iret_write+iret
endif ! if ( klev > 0 )
endif ! if ( mype == k-1 .and. kvar > 0 )
enddo pe_loop
if ( icount > itotflds ) then
lloop=.false.
exit gfsfields
endif
enddo gfsfields
! Print date/time stamp
if ( mype == mype_out ) then
write(6,700) sighead%jcap,grd%nlon,nlatm2,sighead%levs,&
sighead%fhour,sighead%idate
700 format('GENERAL_WRITE_GFSATM: anl write, jcap,lonb,latb,levs=',&
4i6,', hour=',f10.1,', idate=',4i5)
endif
if ( mype < itotflds .or. mype == mype_out ) then
call sigio_rclose(lunges,iret)
call sigio_rclose(lunanl,iret)
iret_write=iret_write+iret
if ( iret_write /= 0 ) then
write(6,*)'GENERAL_WRITE_GFSATM: ***ERROR*** writing ',&
trim(filename),' mype,iret_write=',mype,iret_write
end if
endif
return
end subroutine general_write_gfsatm
subroutine general_gather(grd,g_ps,g_tv,g_vor,g_div,g_q,g_oz,g_cwmr, &
icountx,ivar,ilev,work)
! !USES:
use kinds, only: r_kind,i_kind
use mpimod, only: npe,mpi_comm_world,ierror,mpi_rtype
use general_sub2grid_mod, only: sub2grid_info
use gridmod, only: strip
use constants, only: zero
implicit none
! !INPUT PARAMETERS:
type(sub2grid_info) ,intent(in ) :: grd
integer(i_kind),intent(inout) :: icountx
integer(i_kind),dimension(npe),intent(inout):: ivar,ilev
real(r_kind),dimension(grd%itotsub),intent(out) :: work
! !OUTPUT PARAMETERS:
real(r_kind),dimension(grd%lat2,grd%lon2) ,intent( in) :: g_ps
real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig),intent( in) :: g_tv,&
g_vor,g_div,g_q,g_oz,g_cwmr
! !DESCRIPTION: Transfer contents of 3d subdomains to 2d work arrays over pes
!
! !REVISION HISTORY:
! 2013-06-19 treadon
! 2013-10-24 todling update interface to strip
!
! !REMARKS:
!
! language: f90
! machine: ibm rs/6000 sp; sgi origin 2000; compaq/hp
!
! !AUTHOR:
! kleist org: np23 date: 2013-06-19
!
!EOP
!-------------------------------------------------------------------------
integer(i_kind) klev,k,icount
real(r_kind),dimension(grd%lat1*grd%lon1,npe):: sub
!$omp parallel do schedule(dynamic,1) private(k,klev,icount)
do k=1,npe
icount=icountx+k
if(icount == 1)then
ivar(k)=1
ilev(k)=1
sub(:,k)=zero
else if(icount == 2)then
ivar(k)=2
ilev(k)=1
call strip(g_ps ,sub(:,k))
else if( icount>= 3 .and. icount<=(grd%nsig+2) )then
ivar(k)=3
klev=icount-2
ilev(k)=klev
call strip(g_tv(:,:,klev) ,sub(:,k))
else if( icount>=(grd%nsig)+3 .and. icount<=2*(grd%nsig)+2 )then
ivar(k)=4
klev=icount-2-(grd%nsig)
ilev(k)=klev
call strip(g_vor(:,:,klev) ,sub(:,k))
else if( icount>=2*(grd%nsig)+3 .and. icount<=3*(grd%nsig)+2 )then
ivar(k)=5
klev=icount-2-2*(grd%nsig)
ilev(k)=klev
call strip(g_div(:,:,klev) ,sub(:,k))
else if( icount>=3*(grd%nsig)+3 .and. icount<=4*(grd%nsig)+2 )then
ivar(k)=6
klev=icount-2-3*(grd%nsig)
ilev(k)=klev
call strip(g_q(:,:,klev) ,sub(:,k))
else if( icount>=4*(grd%nsig)+3 .and. icount<=5*(grd%nsig)+2 )then
ivar(k)=7
klev=icount-2-4*(grd%nsig)
ilev(k)=klev
call strip(g_oz(:,:,klev) ,sub(:,k))
else if( icount>=5*(grd%nsig)+3 .and. icount<=6*(grd%nsig)+2 )then
ivar(k)=8
klev=icount-2-5*(grd%nsig)
ilev(k)=klev
call strip(g_cwmr(:,:,klev) ,sub(:,k))
else
! NULL, No work to be done for this pe
ivar(k)=-1
ilev(k)=-1
end if
end do
icountx=icountx+npe
call mpi_alltoallv(sub,grd%isc_g,grd%isd_g,mpi_rtype,&
work,grd%ijn,grd%displs_g,mpi_rtype,&
mpi_comm_world,ierror)
return
end subroutine general_gather