module ncepgfs_io
!$$$ module documentation block
! . . . .
! module: ncepgfs_io
! prgmmr: treadon org: np23 date: 2006-01-10
!
! abstract: This module contains routines which handle input/output
! operations for NCEP GFS atmospheric and surface files.
!
! program history log:
! 2006-01-10 treadon
! 2010-02-20 parrish - make sigio_cnvtdv8 public so can be accessed by general_read_gfsatm, when
! reading in gefs sigma files at resolution different from analysis.
! 2010-03-31 treadon - add read_gfs, use sp_a and sp_b
!
! Subroutines Included:
! sub read_gfs - driver to read ncep gfs atmospheric ("sigma") files
! sub read_gfsatm - read ncep gfs atmospheric ("sigma") file, scatter
! on grid to analysis subdomains
! sub read_gfssfc - read ncep gfs surface file, scatter on grid to
! analysis subdomains
! sub sfc_interpolate - interpolate from gfs atm grid to gfs sfc grid
! sub write_gfs - driver to write ncep gfs atmospheric and surface
! analysis files
! sub write_gfsatm - gather on grid, transform to spectral, write ncep
! gfs atmospheric analysis file
! sub write_gfssfc - gather/write on grid ncep surface analysis file
!
! Variable Definitions:
! none
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
private
public read_sigma
public read_sfc
public read_gfs
public read_gfssfc
public write_gfs
public sfc_interpolate
public sigio_cnvtdv8
contains
subroutine read_gfs(mype)
!$$$ subprogram documentation block
! . . .
! subprogram: read_gfs
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2010-03-31 treadon - create routine
!
! input argument list:
! mype - mpi task id
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind
use gridmod, only: hires_b,sp_a,sp_b
use guess_grids, only: ges_z,ges_ps,ges_vor,ges_div,&
ges_u,ges_v,ges_tv,ges_q,ges_cwmr,ges_oz,&
ifilesig,nfldsig
use gsi_io, only: mype_io
implicit none
integer(i_kind),intent(in ) :: mype
character(24) filename
integer(i_kind):: it,iret
! If hires_b, spectral to grid transform for background
! uses double FFT. Need to pass in sp_a and sp_b
if (hires_b) then
do it=1,nfldsig
write(filename,100) ifilesig(it)
100 format('sigf',i2.2)
call read_gfsatm(filename,mype_io,mype,sp_a,sp_b,&
ges_z(1,1,it),ges_ps(1,1,it),&
ges_vor(1,1,1,it),ges_div(1,1,1,it),&
ges_u(1,1,1,it),ges_v(1,1,1,it),&
ges_tv(1,1,1,it),ges_q(1,1,1,it),&
ges_cwmr(1,1,1,it),ges_oz(1,1,1,it),iret)
end do
! Otherwise, use standard transform. Use sp_a in place of sp_b.
else
do it=1,nfldsig
write(filename,100) ifilesig(it)
call read_gfsatm(filename,mype_io,mype,sp_a,sp_a,&
ges_z(1,1,it),ges_ps(1,1,it),&
ges_vor(1,1,1,it),ges_div(1,1,1,it),&
ges_u(1,1,1,it),ges_v(1,1,1,it),&
ges_tv(1,1,1,it),ges_q(1,1,1,it),&
ges_cwmr(1,1,1,it),ges_oz(1,1,1,it),iret)
end do
endif
end subroutine read_gfs
subroutine read_sigma(lunges,filename,gfshead,sigdata,iope,mype,iret)
!$$$ subprogram documentation block
! . . .
! subprogram: read_sigma
!
! prgrmmr: whitaker
!
! abstract: read a ncep GFS spectral sigma file on a specified task,
! broadcast data to other tasks.
!
! program history log:
! 2012-01-24 whitaker - create routine
!
! input argument list:
! lunges - unit number to use for IO
! mype - mpi task id
! filename - gfs spectral file to read
! iope - mpi task to perform IO
!
! output argument list:
! sigdata (inout) - sigio data structure to hold data
! gfshead (inout) - gfs header structure to hole metadata
! iret - return code (0 for sucessful completion)
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use sigio_module, only: sigio_srohdc,sigio_head,sigio_data
use kinds, only: i_kind,r_single,r_kind
use gridmod, only: ncepgfs_head
use mpimod, only: mpi_integer4,mpi_real4,mpi_comm_world
character(*),intent(in) :: filename
type(sigio_head):: sighead
type(sigio_data), intent(inout):: sigdata
integer(i_kind), intent(inout) :: iret
integer(i_kind), intent(in) :: iope,mype,lunges
type(ncepgfs_head), intent(inout):: gfshead
integer(i_kind) nc,idate(4),levs,ntrac,ncldt,latb,lonb
real(r_single) fhour
! read a file on a specified task, broadcast data to other tasks.
! iope is task that does IO for this file.
if (mype == iope) then
call sigio_srohdc(lunges,filename,sighead,sigdata,iret)
if (iret /= 0) print *,'error in read_sigma',trim(filename),iret
nc = (sighead%jcap+1)*(sighead%jcap+2)
levs = sighead%levs
idate = sighead%idate
ntrac = sighead%ntrac
ncldt = sighead%ncldt
fhour = sighead%fhour
lonb = sighead%lonb
latb = sighead%latb
endif
call mpi_bcast(nc,1,mpi_integer4,iope,mpi_comm_world,iret)
call mpi_bcast(levs,1,mpi_integer4,iope,mpi_comm_world,iret)
call mpi_bcast(idate,4,mpi_integer4,iope,mpi_comm_world,iret)
call mpi_bcast(ntrac,1,mpi_integer4,iope,mpi_comm_world,iret)
call mpi_bcast(ncldt,1,mpi_integer4,iope,mpi_comm_world,iret)
call mpi_bcast(fhour,1,mpi_real4,iope,mpi_comm_world,iret)
call mpi_bcast(lonb,1,mpi_integer4,iope,mpi_comm_world,iret)
call mpi_bcast(latb,1,mpi_integer4,iope,mpi_comm_world,iret)
gfshead%fhour = fhour
gfshead%idate = idate
gfshead%lonb = lonb
gfshead%latb = latb
gfshead%levs = levs
gfshead%ntrac = ntrac
gfshead%ncldt = ncldt
if (mype /= iope) then
! allocate data structure for non-IO tasks.
allocate(sigdata%hs(nc),sigdata%ps(nc),&
sigdata%t(nc,levs),sigdata%d(nc,levs),sigdata%z(nc,levs),&
sigdata%q(nc,levs,ntrac))
endif
call mpi_bcast(sigdata%ps(1),nc,mpi_real4,iope,mpi_comm_world,iret)
call mpi_bcast(sigdata%hs(1),nc,mpi_real4,iope,mpi_comm_world,iret)
call mpi_bcast(sigdata%t(1,1),nc*levs,mpi_real4,iope,mpi_comm_world,iret)
call mpi_bcast(sigdata%z(1,1),nc*levs,mpi_real4,iope,mpi_comm_world,iret)
call mpi_bcast(sigdata%d(1,1),nc*levs,mpi_real4,iope,mpi_comm_world,iret)
call mpi_bcast(sigdata%q(1,1,1),nc*levs*ntrac,mpi_real4,iope,mpi_comm_world,iret)
end subroutine read_sigma
subroutine read_sfc(lunges,filename,sfchead,sfcdata,iope,mype,iret)
!$$$ subprogram documentation block
! . . .
! subprogram: read_sfc
!
! prgrmmr: whitaker
!
! abstract: read a ncep GFS surface file on a specified task,
! broadcast data to other tasks.
!
! program history log:
! 2012-01-24 whitaker - create routine
!
! input argument list:
! lunges - unit number to use for IO
! mype - mpi task id
! filename - gfs surface file to read
! iope - mpi task to perform IO
!
! output argument list:
! sfcdata (inout) - sfc data structure to hold data
! sfchead (inout) - sfc header structure to hold metadata
! iret - return code (0 for sucessful completion)
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
! read data from sfc file on a single task, bcast data to other tasks.
use sfcio_module, only: sfcio_srohdc,sfcio_head,sfcio_data
use kinds, only: i_kind,r_single,r_kind
use mpimod, only: mpi_integer4,mpi_real4,mpi_comm_world
character(*),intent(in) :: filename
type(sfcio_head), intent(inout) :: sfchead
type(sfcio_data), intent(inout):: sfcdata
integer(i_kind), intent(inout) :: iret
integer(i_kind), intent(in) :: iope,mype,lunges
integer(i_kind) idate(4),latb,lonb
real(r_single) fhour
! read a file on a specified task, broadcast data to other tasks.
! iope is task that does IO for this file.
if (mype == iope) then
call sfcio_srohdc(lunges,filename,sfchead,sfcdata,iret)
if (iret /= 0) print *,'error in read_sfc',trim(filename),iret
idate = sfchead%idate
lonb = sfchead%lonb
latb = sfchead%latb
fhour = sfchead%fhour
endif
call mpi_bcast(idate,4,mpi_integer4,iope,mpi_comm_world,iret)
call mpi_bcast(fhour,1,mpi_real4,iope,mpi_comm_world,iret)
call mpi_bcast(lonb,1,mpi_integer4,iope,mpi_comm_world,iret)
call mpi_bcast(latb,1,mpi_integer4,iope,mpi_comm_world,iret)
sfchead%fhour = fhour
sfchead%idate = idate
sfchead%latb = latb
sfchead%lonb = lonb
if (mype /= iope) then
allocate(&
sfcdata%tsea(lonb,latb),&
sfcdata%smc(lonb,latb,1),&
sfcdata%sheleg(lonb,latb),&
sfcdata%stc(lonb,latb,1),&
sfcdata%slmsk(lonb,latb),&
sfcdata%zorl(lonb,latb),&
sfcdata%vfrac(lonb,latb),&
sfcdata%f10m(lonb,latb),&
sfcdata%vtype(lonb,latb),&
sfcdata%stype(lonb,latb),&
sfcdata%orog(lonb,latb))
endif
call mpi_bcast(sfcdata%tsea(1,1),lonb*latb,mpi_real4,iope,mpi_comm_world,iret)
call mpi_bcast(sfcdata%smc(1,1,1),lonb*latb,mpi_real4,iope,mpi_comm_world,iret)
call mpi_bcast(sfcdata%stc(1,1,1),lonb*latb,mpi_real4,iope,mpi_comm_world,iret)
call mpi_bcast(sfcdata%sheleg(1,1),lonb*latb,mpi_real4,iope,mpi_comm_world,iret)
call mpi_bcast(sfcdata%zorl(1,1),lonb*latb,mpi_real4,iope,mpi_comm_world,iret)
call mpi_bcast(sfcdata%vfrac(1,1),lonb*latb,mpi_real4,iope,mpi_comm_world,iret)
call mpi_bcast(sfcdata%slmsk(1,1),lonb*latb,mpi_real4,iope,mpi_comm_world,iret)
call mpi_bcast(sfcdata%f10m(1,1),lonb*latb,mpi_real4,iope,mpi_comm_world,iret)
call mpi_bcast(sfcdata%vtype(1,1),lonb*latb,mpi_real4,iope,mpi_comm_world,iret)
call mpi_bcast(sfcdata%stype(1,1),lonb*latb,mpi_real4,iope,mpi_comm_world,iret)
call mpi_bcast(sfcdata%orog(1,1),lonb*latb,mpi_real4,iope,mpi_comm_world,iret)
end subroutine read_sfc
subroutine read_gfsatm(filename,iope,mype,sp_a,sp_b,g_z,g_ps,g_vor,g_div,g_u,g_v,&
g_tv,g_q,g_cwmr,g_oz,iret_read)
!$$$ subprogram documentation block
! . . . .
! subprogram: read_gfsatm read gfs atm, convert to grid and
! send to all mpi tasks
! prgmmr: parrish org: np22 date: 1990-10-10
!
! abstract: read ncep gfs atmospheric guess, convert to grid, and
! scatter to subdomains
!
! program history log:
! 1990-10-10 parrish
! 1997-09-23 weiyu yang
! 1998-05-15 weiyu yang mpp version
! 1999-08-24 derber, j., treadon, r., yang, w., first frozen mpp version
! 2004-05-18 kleist, documentation
! 2004-05-15 treadon - transform spectral coef to grid,
! communicate grids to all tasks
! 2004-06-17 treadon - update documentation
! 2004-08-03 treadon - add only to module use, add intent in/out
! 2004-08-23 treadon - declare tracers,vtid,pdryini,xncld as real(single)
! 2004-08-27 treadon - use splib routines for grid <---> spectral transforms
! 2005-03-07 dee - support gmao model interface
! 2005-03-30 treadon - clean up formatting of write statement
! 2005-12-09 guo - removed special GMAO spectral input format
! 2006-01-09 treadon - use sigio
! 2006-03-13 treadon - increase filename to 24 characters
! 2006-09-18 treadon - replace lnps with ps
! 2007-05-07 treadon - add gfsio
! 2007-05-08 kleist - add option for lnps or ps
! 2008-05-28 safford - rm unused vars
! 2010-03-10 sela,iredell,lueken - remove hires_b
! 2010-03-18 treadon - remove zonal mean q check
! 2010-03-31 treadon - add sp_a and sp_b
!
! input argument list:
! inges - unit number of guess coefs
! iope - mpi task handling i/o
! mype - mpi task id
!
! output argument list:
! hourg - guess forecast hour
! idateg - initial date of guess
! g_* - guess fields
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,r_single,i_kind
use gridmod, only: displs_s,irc_s,ijn_s,&
ird_s,nsig,nlat,nlon,lat2,lon2,&
itotsub,fill_ns,filluv_ns,ncep_sigio,ncepgfs_head,idpsfc5,idthrm5,&
ntracer,idvc5,cp5,idvm5,reload
use general_specmod, only: spec_vars
use mpimod, only: npe,mpi_comm_world,ierror,mpi_rtype
use constants, only: izero,ione,zero,one,fv
use sigio_module, only: sigio_intkind,sigio_head,sigio_data,&
sigio_srohdc,sigio_axdata
use gfsio_module, only: gfsio_gfile,gfsio_open,gfsio_close,&
gfsio_init,gfsio_getfilehead,gfsio_readrecv
implicit none
! Declare local parameters
integer(sigio_intkind):: lunges = 11
real(r_kind),parameter:: r0_01 = 0.01_r_kind
real(r_kind),parameter:: r0_001 = 0.001_r_kind
real(r_kind),parameter:: qsmall = 1.e-11_r_kind
! Declare passed variables
character(24) ,intent(in ) :: filename
integer(i_kind) , intent(in ) :: iope
integer(i_kind) ,intent(in ) :: mype
integer(i_kind) ,intent( out) :: iret_read
real(r_kind),dimension(lat2,lon2) ,intent( out) :: g_z,g_ps
real(r_kind),dimension(lat2,lon2,nsig),intent( out) :: g_u,g_v,&
g_vor,g_div,g_cwmr,g_q,g_oz,g_tv
type(spec_vars) ,intent(in ) :: sp_a,sp_b
! Declare local variables
integer(i_kind):: iret,nlatm2,ij,n,ii1,l,m
integer(i_kind) i,j,k,icount,icount_prev,mm1
integer(i_kind) mype_hs,mype_ps
real(r_kind),dimension(nlon,nlat-2_i_kind):: grid,grid_u,grid_v,&
grid_vor,grid_div,grid2
real(r_kind),dimension(nlon,nlat-2_i_kind,ntracer):: grid_q
real(r_kind),dimension(sp_b%nc):: spec_work,spec_vor,spec_div
real(r_kind),dimension(sp_a%nc):: spec2_work,spec2_vor,spec2_div
real(r_kind),dimension(itotsub):: work,work_vor,work_div,&
work_u,work_v
real(r_kind),dimension(lat2*lon2,max(2*nsig,npe)):: sub,sub_div,sub_vor,&
sub_u,sub_v
type(sigio_head):: sighead
type(sigio_data):: sigdata
type(gfsio_gfile) :: gfile
type(ncepgfs_head):: gfshead
type:: gfsio_data
real(r_single),allocatable:: hs(:) !surface height, m
real(r_single),allocatable:: ps(:) !surface pressure, pa
real(r_single),allocatable:: t(:,:) !layer temperature, k
real(r_single),allocatable:: u(:,:) !layer zonal wind, m/s
real(r_single),allocatable:: v(:,:) !layer meridional wind, m/s
real(r_single),allocatable:: q(:,:,:) !tracers, 1-spfh; 2-O3; 3-CLW , kg/kg
end type gfsio_data
type(gfsio_data):: gfsdata
!******************************************************************************
! Initialize variables used below
mm1=mype+ione
mype_hs=izero
mype_ps=npe-ione
iret_read=izero
nlatm2=nlat-2_i_kind
! Read NCEP gfs guess file using appropriate io module
! Do IO on task iope, bcast data to other tasks.
call read_sigma(lunges,filename,gfshead,sigdata,iope,mype,iret)
if (iret /= 0) goto 1000
! Process guess fields according to type of input file. NCEP_SIGIO files
! are spectral coefficient files and need to be transformed to the grid.
! NCEP_GFSIO files are grid files. Both formats need to be scattered
! from the full domain to sub-domains.
! Terrain: spectral --> grid transform, scatter to all mpi tasks
if (mype==mype_hs) then
if (ncep_sigio) then
do i=1,sp_b%nc
spec_work(i)=sigdata%hs(i)
if(sp_b%factsml(i))spec_work(i)=zero
end do
call general_sptez_s_b(sp_a,sp_b,spec_work,grid,ione)
else
ij=izero
do j=1,gfshead%latb
do i=1,gfshead%lonb
ij=ij+ione
grid(i,j)=gfsdata%hs(ij)
end do
end do
endif
call fill_ns(grid,work)
endif
call mpi_scatterv(work,ijn_s,displs_s,mpi_rtype,&
g_z,ijn_s(mm1),mpi_rtype,mype_hs,mpi_comm_world,ierror)
! Surface pressure: same procedure as terrain, but handled by task mype_ps
! NCEP SIGIO has two options for surface pressure. Variable idpsfc5
! indicates the type:
! idpsfc5= 0,1 for ln(psfc)
! idpsfc5= 2 for psfc
!
if (mype==mype_ps) then
if (ncep_sigio) then
do i=1,sp_b%nc
spec_work(i)=sigdata%ps(i)
if(sp_b%factsml(i))spec_work(i)=zero
end do
call general_sptez_s_b(sp_a,sp_b,spec_work,grid,ione)
call fill_ns(grid,work)
! If ln(ps), take exponential to convert to ps in cb
if (idpsfc5 /= 2_i_kind) then
do i=1,itotsub
work(i)=exp(work(i))
end do
endif
else
ij=izero
do j=1,gfshead%latb
do i=1,gfshead%lonb
ij=ij+ione
grid(i,j) = r0_001*gfsdata%ps(ij) ! convert Pa to cb
end do
end do
call fill_ns(grid,work)
endif
endif
call mpi_scatterv(work,ijn_s,displs_s,mpi_rtype,&
g_ps,ijn_s(mm1),mpi_rtype,mype_ps,mpi_comm_world,ierror)
! Thermodynamic variable: s-->g transform, communicate to all tasks
! For multilevel fields, each task handles a given level. Periodic
! mpi_alltoallv calls communicate the grids to all mpi tasks.
! Finally, the grids are loaded into guess arrays used later in the
! code.
sub=zero
icount=izero
icount_prev=ione
do k=1,gfshead%levs
icount=icount+ione
if (mype==mod(icount-ione,npe)) then
if (ncep_sigio) then
do i=1,sp_b%nc
spec_work(i)=sigdata%t(i,k)
if(sp_b%factsml(i))spec_work(i)=zero
end do
call general_sptez_s_b(sp_a,sp_b,spec_work,grid,ione)
! SIGIO has three possible thermodynamic variables
! Variable idthrm5 indicates the type
! idthrm5 = 0,1 = virtual temperature (Tv)
! idthrm5 = 2 = sensible (dry) temperature (T)
! idthrm5 = 3 = enthalpy (h=CpT)
! The GSI analysis variable is Tv
! If needed, convert T or h to Tv
if (idthrm5==2_i_kind .or. idthrm5==3_i_kind) then
! Convert tracers from spectral coefficients to grid
do n=1,ntracer
do i=1,sp_b%nc
spec_work(i)=sigdata%q(i,k,n)
if(sp_b%factsml(i))spec_work(i)=zero
end do
call general_sptez_s_b(sp_a,sp_b,spec_work,grid_q(1,1,n),ione)
end do
! Convert input thermodynamic variable to dry temperature
call sigio_cnvtdv8(nlon*nlatm2,nlon*nlatm2,ione,idvc5,&
idvm5,ntracer,iret,grid,grid_q,cp5,ione)
! Convert dry temperature to virtual
do j=1,nlatm2
do i=1,nlon
grid(i,j) = grid(i,j)*(one+fv*grid_q(i,j,1))
end do
end do
endif
else
! GFSIO thermodynamic variable is sensible temperature (T).
! The GSI analysis variable is virtual temperature (Tv).
! Convert T to Tv in the loop below.
ij=izero
do j=1,gfshead%latb
do i=1,gfshead%lonb
ij=ij+ione
grid(i,j)=gfsdata%t(ij,k)*(one+fv*gfsdata%q(ij,k,1))
end do
end do
endif
! Load values into rows for south and north pole
call fill_ns(grid,work)
endif
if (mod(icount,npe)==izero .or. icount==gfshead%levs) then
call mpi_alltoallv(work,ijn_s,displs_s,mpi_rtype,&
sub(1,icount_prev),irc_s,ird_s,mpi_rtype,&
mpi_comm_world,ierror)
icount_prev=icount+ione
endif
end do
call reload(sub,g_tv)
! Divergence and voriticity. Compute u and v from div and vor
sub_vor=zero
sub_div=zero
sub_u=zero
sub_v=zero
icount=izero
icount_prev=ione
do k=1,gfshead%levs
icount=icount+ione
! The work in the loop below is spread over all mpi tasks
if (mype==mod(icount-ione,npe)) then
! Convert spectral coefficients of div and vor to grid space
if (ncep_sigio) then
do i=1,sp_b%nc
spec_div(i)=sigdata%d(i,k) !div
spec_vor(i)=sigdata%z(i,k) !vor
if(sp_b%factvml(i))then
spec_div(i)=zero
spec_vor(i)=zero
end if
end do
call general_sptez_s_b(sp_a,sp_b,spec_div,grid_div,ione)
call general_sptez_s_b(sp_a,sp_b,spec_vor,grid_vor,ione)
call general_sptez_v_b(sp_a,sp_b,spec_div,spec_vor,grid_u,grid_v,ione)
! Convert grid u,v to div and vor
else
ij=izero
do j=1,gfshead%latb
do i=1,gfshead%lonb
ij=ij+ione
grid_u(i,j)=gfsdata%u(ij,k)
grid_v(i,j)=gfsdata%v(ij,k)
end do
end do
call general_sptez_v(sp_b,spec2_div,spec2_vor,grid_u,grid_v,-ione)
call general_sptez_s(sp_b,spec2_div,grid_div,ione)
call general_sptez_s(sp_b,spec2_vor,grid_vor,ione)
endif
call fill_ns(grid_div,work_div)
call fill_ns(grid_vor,work_vor)
call filluv_ns(grid_u,grid_v,work_u,work_v)
endif
! Periodically exchange vor,div,u,v between all mpi tasks.
if (mod(icount,npe)==izero .or. icount==gfshead%levs) then
call mpi_alltoallv(work_vor,ijn_s,displs_s,mpi_rtype,&
sub_vor(1,icount_prev),irc_s,ird_s,mpi_rtype,&
mpi_comm_world,ierror)
call mpi_alltoallv(work_div,ijn_s,displs_s,mpi_rtype,&
sub_div(1,icount_prev),irc_s,ird_s,mpi_rtype,&
mpi_comm_world,ierror)
call mpi_alltoallv(work_u,ijn_s,displs_s,mpi_rtype,&
sub_u(1,icount_prev),irc_s,ird_s,mpi_rtype,&
mpi_comm_world,ierror)
call mpi_alltoallv(work_v,ijn_s,displs_s,mpi_rtype,&
sub_v(1,icount_prev),irc_s,ird_s,mpi_rtype,&
mpi_comm_world,ierror)
icount_prev=icount+ione
endif
end do
! Transfer vor,div,u,v into real(r_kind) guess arrays
call reload(sub_vor,g_vor)
call reload(sub_div,g_div)
call reload(sub_u,g_u)
call reload(sub_v,g_v)
! Specific humidity
sub=zero
icount=izero
icount_prev=ione
do k=1,gfshead%levs
icount=icount+ione
if (mype==mod(icount-ione,npe)) then
if (ncep_sigio) then
do i=1,sp_b%nc
spec_work(i)=sigdata%q(i,k,1)
if(sp_b%factsml(i))spec_work(i)=zero
end do
call general_sptez_s_b(sp_a,sp_b,spec_work,grid,ione)
else
ij=izero
do j=1,gfshead%latb
do i=1,gfshead%lonb
ij=ij+ione
grid(i,j)=gfsdata%q(ij,k,1)
end do
end do
endif
call fill_ns(grid,work)
endif
if (mod(icount,npe)==izero .or. icount==gfshead%levs) then
call mpi_alltoallv(work,ijn_s,displs_s,mpi_rtype,&
sub(1,icount_prev),irc_s,ird_s,mpi_rtype,&
mpi_comm_world,ierror)
icount_prev=icount+ione
endif
end do
call reload(sub,g_q)
! Ozone mixing ratio
sub=zero
icount=izero
icount_prev=ione
do k=1,gfshead%levs
icount=icount+ione
if (mype==mod(icount-ione,npe)) then
if (ncep_sigio) then
do i=1,sp_b%nc
spec_work(i)=sigdata%q(i,k,2)
if(sp_b%factsml(i))spec_work(i)=zero
end do
call general_sptez_s_b(sp_a,sp_b,spec_work,grid,ione)
else
ij=izero
do j=1,gfshead%latb
do i=1,gfshead%lonb
ij=ij+ione
grid(i,j)=gfsdata%q(ij,k,2)
end do
end do
endif
call fill_ns(grid,work)
endif
if (mod(icount,npe)==izero .or. icount==gfshead%levs) then
call mpi_alltoallv(work,ijn_s,displs_s,mpi_rtype,&
sub(1,icount_prev),irc_s,ird_s,mpi_rtype,&
mpi_comm_world,ierror)
icount_prev=icount+ione
endif
end do
call reload(sub,g_oz)
! Cloud condensate mixing ratio.
if (gfshead%ntrac>2_i_kind .or. gfshead%ncldt>=ione) then
sub=zero
icount=izero
icount_prev=ione
do k=1,gfshead%levs
icount=icount+ione
if (mype==mod(icount-ione,npe)) then
if (ncep_sigio) then
do i=1,sp_b%nc
spec_work(i)=sigdata%q(i,k,3)
if(sp_b%factsml(i))spec_work(i)=zero
end do
call general_sptez_s_b(sp_a,sp_b,spec_work,grid,ione)
else
ij=izero
do j=1,gfshead%latb
do i=1,gfshead%lonb
ij=ij+ione
grid(i,j)=gfsdata%q(ij,k,3)
end do
end do
endif
call fill_ns(grid,work)
endif
if (mod(icount,npe)==izero .or. icount==gfshead%levs) then
call mpi_alltoallv(work,ijn_s,displs_s,mpi_rtype,&
sub(1,icount_prev),irc_s,ird_s,mpi_rtype,&
mpi_comm_world,ierror)
icount_prev=icount+ione
endif
end do
call reload(sub,g_cwmr)
else
do k=1,gfshead%levs
do j=1,lon2
do i=1,lat2
g_cwmr(i,j,k)=zero
end do
end do
end do
endif
! Deallocate sigio data array
if (ncep_sigio) call sigio_axdata(sigdata,iret)
iret_read=iret_read+iret
! Print date/time stamp
if(mype==iope) then
write(6,700) gfshead%lonb,gfshead%latb,gfshead%levs,&
gfshead%fhour,gfshead%idate
700 format('READ_GFSATM: ges read/scatter, lonb,latb,levs=',&
3i6,', hour=',f10.1,', idate=',4i5)
end if
return
! ERROR detected while reading file
1000 continue
if (mype==iope) write(6,*)'READ_GFSATM: ***ERROR*** while reading ',&
filename,' from unit ',lunges,'. iret=',iret
call sigio_axdata(sigdata,iret)
iret_read=iret_read+iret
! End of routine. Return
return
end subroutine read_gfsatm
subroutine read_gfssfc(filename,iope,mype,fact10,sfct,sno,veg_type,&
veg_frac,soil_type,soil_temp,soil_moi,isli,sfc_rough,terrain)
!$$$ subprogram documentation block
! . . . .
! subprogram: read_gfssfc read gfs surface file
! prgmmr: treadon org: np23 date: 2003-04-10
!
! abstract: read gfs surface file
!
! program history log:
! 2003-04-10 treadon
! 2004-05-18 kleist, add global isli & documentation
! 2004-09-07 treadon fix mpi bug when npe > nsfc
! 2005-01-27 treadon - rewrite to make use of sfcio module
! 2005-03-07 todling - die gracefully when return error from sfcio
! 2006-09-28 treadon - pull out surface roughness
! 2008-05-28 safford - rm unused vars
! 2009-01-12 gayno - add read of terrain height
!
! input argument list:
! filename - name of surface guess file
! iope - mpi task handling i/o
! mype - mpi task id
!
! output argument list:
! fact10 - 10 meter wind factor
! sfct - surface temperature (skin temp)
! sno - snow depth
! veg_type - vegetation type
! veg_frac - vegetation fraction
! soil_type - soil type
! soil_temp - soil temperature of first layer
! soil_moi - soil moisture of first layer
! isli - sea/land/ice mask (subdomain)
! isli_g - global sea/land/ice mask
! sfc_rough - surface roughness
! terrain - terrain height
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind
use gridmod, only: nlat_sfc,nlon_sfc
use sfcio_module, only: sfcio_intkind,sfcio_head,sfcio_data,&
sfcio_srohdc,sfcio_axdata
use constants, only: izero,ione,zero
implicit none
! Declare passed variables
character(24) ,intent(in ) :: filename
integer(i_kind) ,intent(in ) :: iope
integer(i_kind) ,intent(in ) :: mype
integer(i_kind),dimension(nlat_sfc,nlon_sfc),intent( out) :: isli
real(r_kind) ,dimension(nlat_sfc,nlon_sfc),intent( out) :: fact10,sfct,sno,&
veg_type,veg_frac,soil_type,soil_temp,soil_moi,sfc_rough,terrain
! Declare local parameters
integer(sfcio_intkind):: lunges = 11
integer(i_kind),parameter:: nsfc=11_i_kind
! Declare local variables
integer(i_kind) i,j,k,latb,lonb,mm1
integer(sfcio_intkind):: irets
real(r_kind) sumn,sums
real(r_kind),allocatable,dimension(:,:,:):: work,sfcges
type(sfcio_head):: sfc_head
type(sfcio_data):: sfc_data
mm1=mype+ione
!-----------------------------------------------------------------------------
! Read surface file
call read_sfc(lunges,filename,sfc_head,sfc_data,iope,mype,irets)
! Check for possible problems
if (irets /= izero) then
write(6,*)'READ_GFSSFC: ***ERROR*** problem reading ',filename,&
', irets=',irets
call sfcio_axdata(sfc_data,irets)
call stop2(80)
endif
latb=sfc_head%latb
lonb=sfc_head%lonb
if ( (latb /= nlat_sfc-2_i_kind) .or. &
(lonb /= nlon_sfc) ) then
write(6,*)'READ_GFSSFC: ***ERROR*** inconsistent grid dimensions. ',&
', nlon,nlat-2=',nlon_sfc,nlat_sfc-2_i_kind,' -vs- sfc file lonb,latb=',&
lonb,latb
call sfcio_axdata(sfc_data,irets)
call stop2(80)
endif
! Load surface fields into local work array
allocate(work(lonb,latb,nsfc),sfcges(latb+2_i_kind,lonb,nsfc))
do k=1,nsfc
do j=1,latb
do i=1,lonb
work(i,j,k) = zero
end do
end do
end do
do j=1,latb
do i=1,lonb
work(i,j,1) = sfc_data%tsea (i,j) ! skin temperature
work(i,j,2) = sfc_data%smc (i,j,1) ! soil moisture
work(i,j,3) = sfc_data%sheleg(i,j) ! snow depth
work(i,j,4) = sfc_data%stc (i,j,1) ! soil temperature
work(i,j,5) = sfc_data%slmsk (i,j) ! sea/land/ice mask
work(i,j,6) = sfc_data%vfrac (i,j) ! vegetation cover
work(i,j,7) = sfc_data%f10m (i,j) ! 10m wind factor
work(i,j,8) = sfc_data%vtype (i,j) ! vegetation type
work(i,j,9) = sfc_data%stype (i,j) ! soil type
work(i,j,10) = sfc_data%zorl (i,j) ! surface roughness length (cm)
work(i,j,11) = sfc_data%orog (i,j) ! terrain
end do
end do
! Fill surface guess array
do k=1,nsfc
! Compute mean for southern- and northern-most rows
! of surface guess array
sumn = zero
sums = zero
do i=1,lonb
sumn = work(i,1,k) + sumn
sums = work(i,latb,k) + sums
end do
sumn = sumn/float(lonb)
sums = sums/float(lonb)
! Transfer from local work array to surface guess array
do j = 1,lonb
sfcges(1,j,k)=sums
sfcges(latb+2_i_kind,j,k)=sumn
do i=2,latb+ione
sfcges(i,j,k) = work(j,latb+2_i_kind-i,k)
end do
end do
! End of loop over data records
end do
! Deallocate local work arrays
deallocate(work)
! Load data into output arrays
do j=1,lonb
do i=1,latb+2_i_kind
sfct(i,j) = sfcges(i,j,1)
soil_moi(i,j) = sfcges(i,j,2)
sno(i,j) = sfcges(i,j,3)
soil_temp(i,j) = sfcges(i,j,4)
isli(i,j) = nint(sfcges(i,j,5)+0.0000001_r_kind)
veg_frac(i,j) = sfcges(i,j,6)
fact10(i,j) = sfcges(i,j,7)
veg_type(i,j) = sfcges(i,j,8)
soil_type(i,j) = sfcges(i,j,9)
sfc_rough(i,j) = sfcges(i,j,10)
terrain(i,j) = sfcges(i,j,11)
end do
end do
deallocate(sfcges)
! Print date/time stamp
if(mype==iope) then
write(6,700) latb,lonb,sfc_head%fhour,sfc_head%idate
700 format('READ_GFSSFC: ges read/scatter, nlat,nlon=',&
2i6,', hour=',f10.1,', idate=',4i5)
end if
return
end subroutine read_gfssfc
subroutine write_gfs(increment,mype,mype_atm,mype_sfc)
!$$$ subprogram documentation block
! . . .
! subprogram: write_gfs
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2006-07-31 kleist - pass ges_ps instead of ges_lnps
! 2006-10-11 treadon - update 10m wind factor in sfc file
! 2008-05-28 safford - rm unused vars, add doc block
! 2008-12-05 todling - adjustment for dsfct time dimension addition
! 2009-11-28 todling - add increment option (hook-only for now)
! 2010-03-31 treadon - add hires_b, sp_a, and sp_b
!
! input argument list:
! increment - when .t. will name files as increment files
! mype - mpi task id
! mype_atm,mype_sfc -
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind
use constants, only: izero
use guess_grids, only: ges_z,ges_ps,ges_vor,ges_div,&
ges_tv,ges_q,ges_oz,ges_cwmr,ges_prsl,&
ges_u,ges_v,ges_prsi,dsfct
use guess_grids, only: ntguessig,ntguessfc
use gridmod, only: hires_b,sp_a,sp_b
implicit none
logical ,intent(in ) :: increment
integer(i_kind),intent(in ) :: mype,mype_atm,mype_sfc
character(24):: filename
! Write atmospheric analysis file
if (increment) then
filename='siginc'
if(mype==izero) write(6,'(a)') 'WRITE_GFS: sorry, I do not know how to write increments yet'
return
else
filename='siganl'
endif
! If hires_b, spectral to grid transform for background
! uses double FFT. Need to pass in sp_a and sp_b
if (hires_b) then
call write_gfsatm(filename,mype,mype_atm,&
sp_a,sp_b,&
ges_z(1,1,ntguessig),ges_ps(1,1,ntguessig),&
ges_vor(1,1,1,ntguessig),ges_div(1,1,1,ntguessig),&
ges_tv(1,1,1,ntguessig),ges_q(1,1,1,ntguessig),&
ges_oz(1,1,1,ntguessig),ges_cwmr(1,1,1,ntguessig),&
ges_prsl(1,1,1,ntguessig),ges_u(1,1,1,ntguessig),&
ges_v(1,1,1,ntguessig),ges_prsi(1,1,1,ntguessig))
! Otherwise, use standard transform. Use sp_a in place of sp_b.
else
call write_gfsatm(filename,mype,mype_atm,&
sp_a,sp_a,&
ges_z(1,1,ntguessig),ges_ps(1,1,ntguessig),&
ges_vor(1,1,1,ntguessig),ges_div(1,1,1,ntguessig),&
ges_tv(1,1,1,ntguessig),ges_q(1,1,1,ntguessig),&
ges_oz(1,1,1,ntguessig),ges_cwmr(1,1,1,ntguessig),&
ges_prsl(1,1,1,ntguessig),ges_u(1,1,1,ntguessig),&
ges_v(1,1,1,ntguessig),ges_prsi(1,1,1,ntguessig))
endif
! Write surface analysis file
if (increment) then
filename='sfcinc.gsi'
else
filename='sfcanl.gsi'
endif
call write_gfssfc(filename,mype,mype_sfc,dsfct(1,1,ntguessfc))
end subroutine write_gfs
subroutine write_gfsatm(filename,mype,mype_out,sp_a,sp_b,sub_z,sub_ps,&
sub_vor,sub_div,sub_tv,sub_q,sub_oz,sub_cwmr,sub_prsl,&
sub_u,sub_v,sub_prsi)
!$$$ subprogram documentation block
! . . .
! subprogram: write_gfsatm --- Gather, transform, and write out
! spectal coefficients
! prgrmmr: parrish - author; org: np22
!
! abstract: 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.
!
! program history log:
! 1998-07-10 weiyu yang
! 1999-08-24 derber, j., treadon, r., yang, w., first frozen mpp version
! 2003-10-31 kleist, d. - add capability to generate output file for
! either hybrid or sigma vertical coordinate
! 2004-06-15 treadon - update documentation
! 2004-07-15 todling - protex-compliant prologue; added intent/only's
! 2004-08-27 treadon - use splib routine for grid <---> spectral transforms
! 2005-03-07 dee - support gmao model interface
! 2005-03-10 treadon - remove iadate from calling list, access via obsmod
! 2005-04-05 wgu - bug fix: modified iadate not properly merge w/ gmao_intfc case
! 2005-10-13 treadon - properly specify vcid4 in NCEP sigma file header
! 2005-12-09 guo - removed special GMAO spectral output format
! 2006-01-09 treadon - use sigio
! 2006-09-18 treadon - convert ps to lnps
! 2007-05-07 treadon - add gfsio
! 2007-05-08 kleist - add options for ps or lnps
! 2008-05-28 safford - rm unused vars and uses
! 2009-06-11 kleist - add sppad for multiple spectral resolutions
! 2010-03-10 sela,iredell,lueken - remove hires_b
! 2010-03-31 treadon - add sp_a and sp_b
!
! input argument list:
! filename - file to open and write to
! mype - mpi task number
! mype_out - mpi task to write output file
! sub_z - GFS terrain field on subdomains
! sub_ps - surface pressure on subdomains
! sub_vor - vorticity on subdomains
! sub_div - divergence on subdomains
! sub_tv - virtual temperature on subdomains
! sub_q - specific humidity on subdomains
! sub_oz - ozone on subdomains
! sub_cwmr - cloud condensate mixing ratio on subdomains
! sub_prsl - layer midpoint pressure
! sub_u - zonal wind
! sub_v - meridional wind
! sub_prsi - interface pressure
!
! output argument list:
!
! attributes:
! language: f90
! machines: ibm RS/6000 SP; SGI Origin 2000; Compaq HP
!
!$$$ end documentation block
! !USES:
use kinds, only: r_kind,i_kind,r_single
use constants, only: izero,ione,zero_single,r1000,fv,one,zero
use mpimod, only: mpi_rtype,mpi_rtype4
use mpimod, only: mpi_comm_world
use mpimod, only: ierror
use mpimod, only: mpi_status_size
use mpimod, only: npe
use guess_grids, only: ntguessig,ifilesig
use gridmod, only: nlat, nlon ! no. lat/lon on analysis grid
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: load_grid
use gridmod, only: idpsfc5 ! identifier for ps vs ln(ps)
use gridmod, only: idthrm5 ! identifier for thermodynamic variable
use gridmod, only: cp5
use gridmod, only: idvc5
use gridmod, only: idvm5
use gridmod, only: ntracer
use gridmod, only: ncloud
use gridmod, only: ncep_sigio
use gridmod, only: ncepgfs_head
use gridmod, only: strip
use obsmod, only: iadate
use general_specmod, only: spec_vars
use sigio_module, only: sigio_intkind,sigio_head,sigio_data,&
sigio_swopen,sigio_swhead,sigio_swdata,sigio_axdata,&
sigio_srohdc,sigio_realkind
use gfsio_module, only: gfsio_gfile,gfsio_open,gfsio_init,&
gfsio_getfilehead,gfsio_close,gfsio_writerecv
implicit none
! !INPUT PARAMETERS:
character(24) ,intent(in ) :: filename ! file to open and write to
integer(i_kind) ,intent(in ) :: mype ! mpi task number
integer(i_kind) ,intent(in ) :: mype_out ! mpi task to write output file
real(r_kind),dimension(lat2,lon2) ,intent(in ) :: sub_z ! GFS terrain field on subdomains
real(r_kind),dimension(lat2,lon2) ,intent(in ) :: sub_ps ! surface pressure on subdomains
real(r_kind),dimension(lat2,lon2,nsig) ,intent(in ) :: sub_vor ! vorticity on subdomains
real(r_kind),dimension(lat2,lon2,nsig) ,intent(in ) :: sub_div ! divergence on subdomains
real(r_kind),dimension(lat2,lon2,nsig) ,intent(in ) :: sub_tv ! virtual temperature on subdomains
real(r_kind),dimension(lat2,lon2,nsig) ,intent(in ) :: sub_q ! specific humidity on subdomains
real(r_kind),dimension(lat2,lon2,nsig) ,intent(in ) :: sub_oz ! ozone on subdomains
real(r_kind),dimension(lat2,lon2,nsig) ,intent(in ) :: sub_cwmr ! cloud condensate mixing ratio on subdomains
real(r_kind),dimension(lat2,lon2,nsig) ,intent(in ) :: sub_prsl ! layer midpoint pressure
real(r_kind),dimension(lat2,lon2,nsig) ,intent(in ) :: sub_u ! zonal wind
real(r_kind),dimension(lat2,lon2,nsig) ,intent(in ) :: sub_v ! meridional wind
real(r_kind),dimension(lat2,lon2,nsig+ione),intent(in ) :: sub_prsi ! interface pressure
type(spec_vars) ,intent(in ) :: sp_a
type(spec_vars) ,intent(in ) :: sp_b
!-------------------------------------------------------------------------
integer(i_kind),parameter:: lunges = 11_i_kind
integer(i_kind),parameter:: lunanl = 51_i_kind
character(5):: string
character(6):: fname_ges
integer(i_kind) i,j,ij,k,mm1,nlatm2
integer(sigio_intkind):: iret
real(r_kind),dimension(lat1*lon1):: hsm,psm
real(r_kind),dimension(lat2,lon2,nsig):: sub_dp
real(r_kind),dimension(lat1*lon1,nsig):: tvsm,vorsm,divsm
real(r_kind),dimension(lat1*lon1,nsig):: prslm,usm,vsm,dpsm
real(r_kind),dimension(lat1*lon1,nsig,ntracer):: qsm
real(r_kind),dimension(max(iglobal,itotsub)):: work1
real(r_kind),dimension(max(iglobal,itotsub),nsig):: work1_k
real(r_kind),dimension(nlon,nlat-2_i_kind):: grid,grid2
real(r_kind),dimension(sp_b%nc):: spec_work
real(r_kind),dimension(sp_a%nc):: spec_work_sm
type(sigio_head):: sighead
type(sigio_data):: sigdata
type(gfsio_gfile) :: gfilei,gfileo
type(ncepgfs_head) :: gfshead
type:: gfsio_data
real(r_single),allocatable:: hs(:) !surface height, m
real(r_single),allocatable:: ps(:) !surface pressure, pa
real(r_single),allocatable:: p(:,:) !layer pressure, pa
real(r_single),allocatable:: dp(:,:) !layer pressure thickness, pa
real(r_single),allocatable:: t(:,:) !layer temperature, k
real(r_single),allocatable:: u(:,:) !layer zonal wind, m/s
real(r_single),allocatable:: v(:,:) !layer meridional wind, m/s
real(r_single),allocatable:: q(:,:,:) !tracers, 1-spfh; 2-O3; 3-CLW , kg/kg
end type gfsio_data
type(gfsio_data) :: gfsdata
integer(i_kind) :: mype_th,mype_sh,mype_oz,mype_clc,mype_div,mype_vort
integer(i_kind) :: itag_th,itag_sh,itag_oz,itag_clc,itag_div,itag_vort
integer(i_kind) :: status(mpi_status_size),istat,pe_stride
real(kind=sigio_realkind),allocatable :: temp(:,:)
!*************************************************************************
! Initialize local variables
mm1=mype+ione
nlatm2=nlat-2_i_kind
! Set mpi tasks and tags for PE's to do grid transformations
pe_stride = max(izero,(npe-2_i_kind)/6)
mype_th = min(2_i_kind,npe-ione)
itag_th = 10000_i_kind
mype_sh = mype_th+pe_stride
itag_sh = 10001_i_kind
mype_oz = mype_sh+pe_stride
itag_oz = 10002_i_kind
mype_clc = mype_oz+pe_stride
itag_clc = 10003_i_kind
mype_div = mype_clc+pe_stride
itag_div = 10004_i_kind
mype_vort = mype_div+pe_stride
itag_vort = 10005_i_kind
! Strip off boundary points from subdomains
call strip(sub_z ,hsm ,ione)
call strip(sub_ps ,psm ,ione)
call strip(sub_vor ,vorsm ,nsig)
call strip(sub_div ,divsm ,nsig)
call strip(sub_tv ,tvsm ,nsig)
call strip(sub_q ,qsm(1,1,1),nsig)
call strip(sub_oz ,qsm(1,1,2),nsig)
call strip(sub_cwmr,qsm(1,1,3),nsig)
! Read guess file. Pull out header information.
! Update header and write out to analysis file.
! These operations only need to be done on the
! analysis output task, mype_out
if (ncep_sigio) then
if (mype==mype_out) then
! Set guess file name
write(fname_ges,100) ifilesig(ntguessig)
100 format('sigf',i2.2)
! Handle case of NCEP SIGIO
! Read header and spectral coefficients from guess
call sigio_srohdc(lunges,fname_ges,sighead,sigdata,iret)
!send data to compute pes
call mpi_send(sigdata%t,sp_b%nc*nsig,mpi_rtype4,mype_th,&
itag_th,mpi_comm_world,ierror)
call mpi_send(sigdata%q(1,1,1),sp_b%nc*nsig,mpi_rtype4,mype_sh,&
itag_sh,mpi_comm_world,ierror)
call mpi_send(sigdata%q(1,1,2),sp_b%nc*nsig,mpi_rtype4,mype_oz,&
itag_oz,mpi_comm_world,ierror)
if (ntracer>2_i_kind .or. ncloud>=ione) then
call mpi_send(sigdata%q(1,1,3),sp_b%nc*nsig,mpi_rtype4,mype_clc,&
itag_clc,mpi_comm_world,ierror)
endif
call mpi_send(sigdata%d,sp_b%nc*nsig,mpi_rtype4,mype_div,&
itag_div,mpi_comm_world,ierror)
call mpi_send(sigdata%z,sp_b%nc*nsig,mpi_rtype4,mype_vort,&
itag_vort,mpi_comm_world,ierror)
!
! Replace header record date with analysis time
sighead%fhour = zero_single
sighead%idate(1) = iadate(4) !hour
sighead%idate(2) = iadate(2) !month
sighead%idate(3) = iadate(3) !day
sighead%idate(4) = iadate(1) !year
! Load grid dimension and other variables used below
! into local header structure
gfshead%fhour = sighead%fhour
gfshead%idate = sighead%idate
gfshead%levs = sighead%levs
gfshead%ntrac = sighead%ntrac
gfshead%ncldt = sighead%ncldt
gfshead%jcap = sighead%jcap
gfshead%lonb = nlon
gfshead%latb = nlatm2
gfshead%idrt = 4_i_kind
! Write header to analysis file
call sigio_swopen(lunanl,filename,iret)
call sigio_swhead(lunanl,sighead,iret)
else
if (mype==mype_th) then
allocate (temp(sp_b%nc,nsig),stat=istat)
call mpi_recv(temp,sp_b%nc*nsig,mpi_rtype4,mype_out,&
itag_th,mpi_comm_world,status,ierror)
endif
if (mype==mype_sh) then
allocate (temp(sp_b%nc,nsig),stat=istat)
call mpi_recv(temp,sp_b%nc*nsig,mpi_rtype4,mype_out,&
itag_sh,mpi_comm_world,status,ierror)
endif
if (mype==mype_oz) then
allocate (temp(sp_b%nc,nsig),stat=istat)
call mpi_recv(temp,sp_b%nc*nsig,mpi_rtype4,mype_out,&
itag_oz,mpi_comm_world,status,ierror)
endif
if (mype==mype_clc) then
if (ntracer>2_i_kind .or. ncloud>=ione) then
allocate (temp(sp_b%nc,nsig),stat=istat)
call mpi_recv(temp,sp_b%nc*nsig,mpi_rtype4,mype_out,&
itag_clc,mpi_comm_world,status,ierror)
endif
endif
if (mype==mype_div) then
allocate (temp(sp_b%nc,nsig),stat=istat)
call mpi_recv(temp,sp_b%nc*nsig,mpi_rtype4,mype_out,&
itag_div,mpi_comm_world,status,ierror)
endif
if (mype==mype_vort) then
allocate (temp(sp_b%nc,nsig),stat=istat)
call mpi_recv(temp,sp_b%nc*nsig,mpi_rtype4,mype_out,&
itag_vort,mpi_comm_world,status,ierror)
endif
endif
! Handle case of NCEP GFSIO
else
if (mype==mype_out) then
do k=1,nsig
do j=1,lon2
do i=1,lat2
sub_dp(i,j,k) = sub_prsi(i,j,k)-sub_prsi(i,j,k+1)
end do
end do
end do
call strip(sub_dp,dpsm,nsig)
call strip(sub_prsl,prslm,nsig)
call strip(sub_u,usm,nsig)
call strip(sub_v,vsm,nsig)
! Read header record from guess
call gfsio_init(iret)
call gfsio_open(gfilei,fname_ges,'read',iret)
call gfsio_getfilehead(gfilei,&
version=gfshead%version,&
fhour=gfshead%fhour,&
idate=gfshead%idate,&
levs=gfshead%levs, &
ntrac=gfshead%ntrac, &
ncldt=gfshead%ncldt, &
lonb=gfshead%lonb, &
latb=gfshead%latb, &
jcap=gfshead%jcap, &
idrt=gfshead%idrt, &
iret=iret)
gfileo=gfilei
call gfsio_close(gfilei,iret)
! Replace header record date with analysis time
gfshead%fhour = zero_single
gfshead%idate(1) = iadate(4) !hour
gfshead%idate(3) = iadate(3) !month
gfshead%idate(2) = iadate(2) !day
gfshead%idate(4) = iadate(1) !year
! Write header record to gfsio output file
call gfsio_init(iret)
call gfsio_open(gfileo,filename,'write',&
fhour=gfshead%fhour,&
idate=gfshead%idate,&
iret=iret)
if (iret/=izero) then
write(6,*)'WRITE_GFSATM: ***ERROR*** with write to gfsio_open ',&
filename,iret
call stop2(99)
endif
! Allocate structure arrays to hold data
allocate(gfsdata%hs(gfshead%latb*gfshead%lonb))
allocate(gfsdata%ps(gfshead%latb*gfshead%lonb))
allocate(gfsdata%p(gfshead%latb*gfshead%lonb,gfshead%levs))
allocate(gfsdata%dp(gfshead%latb*gfshead%lonb,gfshead%levs))
allocate(gfsdata%t(gfshead%latb*gfshead%lonb,gfshead%levs))
allocate(gfsdata%u(gfshead%latb*gfshead%lonb,gfshead%levs))
allocate(gfsdata%v(gfshead%latb*gfshead%lonb,gfshead%levs))
allocate(gfsdata%q(gfshead%latb*gfshead%lonb,gfshead%levs,gfshead%ntrac))
endif
end if
!gather the fields on the processors that will perform grid transforms
! Thermodynamic variable
! The GSI analysis variable is virtual temperature (Tv). For GFSIO
! output we need the sensible temperature. For SIGIO we have three
! possibilities: Tv, sensible temperature (T), or enthalpy (h=CpT)
if (ncep_sigio) then
if (idthrm5==2_i_kind .or. idthrm5==3_i_kind) then
! Convert Tv to T
do k=1,nsig
do i=1,lat1*lon1
tvsm(i,k)=tvsm(i,k)/(one+fv*qsm(i,k,1))
end do
end do
! If CpT is requested, call function to make conversion
if (idthrm5==3_i_kind) call sigio_cnvtdv8(lat1*lon1,lat1*lon1,&
nsig,idvc5,idvm5,ntracer,iret,tvsm,qsm,cp5,-ione)
endif
else
! Handle the case of GFSIO. Convert Tv to T
do k=1,nsig
do i=1,lat1*lon1
tvsm(i,k) = tvsm(i,k)/(one+fv*qsm(i,k,1))
end do
end do
endif
! create global grid by gathering from subdomains
if (ncep_sigio) then
do k=1,nsig
call mpi_gatherv(tvsm(1,k),ijn(mm1),mpi_rtype,&
work1_k(1,k),ijn,displs_g,mpi_rtype,&
mype_th,mpi_comm_world,ierror)
end do
! Specific humidity
do k=1,nsig
call mpi_gatherv(qsm(1,k,1),ijn(mm1),mpi_rtype,&
work1_k(1,k),ijn,displs_g,mpi_rtype,&
mype_sh,mpi_comm_world,ierror)
end do
! Ozone
do k=1,nsig
call mpi_gatherv(qsm(1,k,2),ijn(mm1),mpi_rtype,&
work1_k(1,k),ijn,displs_g,mpi_rtype,&
mype_oz,mpi_comm_world,ierror)
end do
! Cloud condensate mixing ratio
if (ntracer>2_i_kind .or. ncloud>=ione) then
do k=1,nsig
call mpi_gatherv(qsm(1,k,3),ijn(mm1),mpi_rtype,&
work1_k(1,k),ijn,displs_g,mpi_rtype,&
mype_clc,mpi_comm_world,ierror)
end do
endif
! Horizontal divergence and voriticy
do k=1,nsig
call mpi_gatherv(divsm(1,k),ijn(mm1),mpi_rtype,&
work1_k(1,k),ijn,displs_g,mpi_rtype,&
mype_div,mpi_comm_world,ierror)
end do
do k=1,nsig
call mpi_gatherv(vorsm(1,k),ijn(mm1),mpi_rtype,&
work1_k(1,k),ijn,displs_g,mpi_rtype,&
mype_vort,mpi_comm_world,ierror)
end do
endif
! Generate and write analysis fields
! For each output grid, the following steps are repeated
! 1) create global grid by gathering from subdomains
! 2) transfrom from grid space representation to spectral coefficients
! 3) apply factor to ensure certain coefficients are zero
! 4) write spectral coefficients to output file
! Note that steps 2-4 are done on a single task (here mpi task 0)
! Terrain
call mpi_gatherv(hsm,ijn(mm1),mpi_rtype,&
work1,ijn,displs_g,mpi_rtype,&
mype_out,mpi_comm_world,ierror)
if (mype==mype_out) then
call load_grid(work1,grid)
if (ncep_sigio) then
do i=1,sp_b%nc
spec_work(i) = sigdata%hs(i)
if(sp_b%factsml(i))spec_work(i)=zero
end do
call general_sptez_s_b(sp_a,sp_b,spec_work,grid2,ione)
grid=grid-grid2
call general_sptez_s(sp_a,spec_work_sm,grid,-ione)
call sppad(izero,sp_a%jcap,spec_work_sm,izero,sp_b%jcap,spec_work)
do i=1,sp_b%nc
sigdata%hs(i)=sigdata%hs(i)+spec_work(i)
if(sp_b%factsml(i))sigdata%hs(i)=zero_single
end do
else
ij=izero
do j=1,nlatm2
do i=1,nlon
ij=ij+ione
gfsdata%hs(ij)=grid(i,j)
end do
end do
call gfsio_writerecv(gfileo,'hgt','sfc',ione,gfsdata%hs,iret,idrt=gfshead%idrt)
endif
endif
! Surface pressure.
! NCEP SIGIO outputs surface pressure or ln(surface pressure)
if (ncep_sigio) then
if (idpsfc5 /= 2_i_kind) then
do i=1,lat1*lon1
psm(i)=log(psm(i))
end do
endif
endif
call mpi_gatherv(psm,ijn(mm1),mpi_rtype,&
work1,ijn,displs_g,mpi_rtype,&
mype_out,mpi_comm_world,ierror)
if (mype==mype_out) then
call load_grid(work1,grid)
if (ncep_sigio) then
do i=1,sp_b%nc
spec_work(i) = sigdata%ps(i)
if(sp_b%factsml(i))spec_work(i)=zero
end do
call general_sptez_s_b(sp_a,sp_b,spec_work,grid2,ione)
grid=grid-grid2
call general_sptez_s(sp_a,spec_work_sm,grid,-ione)
call sppad(izero,sp_a%jcap,spec_work_sm,izero,sp_b%jcap,spec_work)
do i=1,sp_b%nc
sigdata%ps(i)=sigdata%ps(i)+spec_work(i)
if(sp_b%factsml(i))sigdata%ps(i)=zero_single
end do
else
ij=izero
do j=1,nlatm2
do i=1,nlon
ij=ij+ione
gfsdata%ps(ij)=grid(i,j)*r1000
end do
end do
call gfsio_writerecv(gfileo,'pres','sfc',ione,gfsdata%ps,iret,idrt=gfshead%idrt)
endif
endif
! Thermodynamic variable
if (ncep_sigio) then
if (mype==mype_th) then
!!$omp parallel do private(k,grid,i,spec_work,grid2,spec_work_sm)
do k=1,nsig
call load_grid(work1_k(1,k),grid)
do i=1,sp_b%nc
spec_work(i) = temp(i,k)
if(sp_b%factsml(i))spec_work(i)=zero
end do
call general_sptez_s_b(sp_a,sp_b,spec_work,grid2,ione)
grid=grid-grid2
call general_sptez_s(sp_a,spec_work_sm,grid,-ione)
call sppad(izero,sp_a%jcap,spec_work_sm,izero,sp_b%jcap,spec_work)
do i=1,sp_b%nc
temp(i,k)=temp(i,k)+spec_work(i)
if(sp_b%factsml(i))temp(i,k)=zero_single
end do
end do
!!$omp end parallel do
!send temperature back to mype_out
call mpi_send(temp,sp_b%nc*nsig,mpi_rtype4,mype_out,&
itag_th,mpi_comm_world,ierror)
endif
else !GFS I/O
do k=1,nsig
call mpi_gatherv(tvsm(1,k),ijn(mm1),mpi_rtype,&
work1,ijn,displs_g,mpi_rtype,&
mype_out,mpi_comm_world,ierror)
if (mype == mype_out) then
call load_grid(work1,grid)
ij=izero
do j=1,nlatm2
do i=1,nlon
ij=ij+ione
gfsdata%t(ij,k)=grid(i,j)
end do
end do
call gfsio_writerecv(gfileo,'tmp','layer',k,gfsdata%t(:,k),iret,&
idrt=gfshead%idrt)
endif
end do
endif
! Specific humidity
if (ncep_sigio) then
if (mype==mype_sh) then
!!$omp parallel do private(k,grid,i,spec_work,grid2,spec_work_sm)
do k=1,nsig
call load_grid(work1_k(1,k),grid)
do i=1,sp_b%nc
spec_work(i) = temp(i,k)
if(sp_b%factsml(i))spec_work(i)=zero
end do
call general_sptez_s_b(sp_a,sp_b,spec_work,grid2,ione)
grid=grid-grid2
call general_sptez_s(sp_a,spec_work_sm,grid,-ione)
call sppad(izero,sp_a%jcap,spec_work_sm,izero,sp_b%jcap,spec_work)
do i=1,sp_b%nc
temp(i,k) =temp(i,k)+spec_work(i)
if(sp_b%factsml(i))temp(i,k)=zero_single
end do
end do
!!$omp end parallel do
!send sh back to mype_out
call mpi_send(temp,sp_b%nc*nsig,mpi_rtype4,mype_out,&
itag_sh,mpi_comm_world,ierror)
endif
else !GFS I/O
do k=1,nsig
call mpi_gatherv(qsm(1,k,1),ijn(mm1),mpi_rtype,&
work1,ijn,displs_g,mpi_rtype,&
mype_out,mpi_comm_world,ierror)
if (mype == mype_out) then
call load_grid(work1,grid)
ij=izero
do j=1,nlatm2
do i=1,nlon
ij=ij+ione
gfsdata%q(ij,k,1)=grid(i,j)
end do
end do
call gfsio_writerecv(gfileo,'spfh','layer',k,gfsdata%q(:,k,1),iret,&
idrt=gfshead%idrt)
endif
end do
endif
! Ozone
if (ncep_sigio) then
if (mype==mype_oz) then
!!$omp parallel do private(k,grid,i,spec_work,grid2,spec_work_sm)
do k=1,nsig
call load_grid(work1_k(1,k),grid)
do i=1,sp_b%nc
spec_work(i) = temp(i,k)
if(sp_b%factsml(i))spec_work(i)=zero
end do
call general_sptez_s_b(sp_a,sp_b,spec_work,grid2,ione)
grid=grid-grid2
call general_sptez_s(sp_a,spec_work_sm,grid,-ione)
call sppad(izero,sp_a%jcap,spec_work_sm,izero,sp_b%jcap,spec_work)
do i=1,sp_b%nc
temp(i,k) =temp(i,k) + spec_work(i)
if(sp_b%factsml(i))temp(i,k)=zero_single
end do
end do
!!$omp end parallel do
!send sh back to mype_out
call mpi_send(temp,sp_b%nc*nsig,mpi_rtype4,mype_out,&
itag_oz,mpi_comm_world,ierror)
endif
else !GFS I/O
do k=1,nsig
call mpi_gatherv(qsm(1,k,2),ijn(mm1),mpi_rtype,&
work1,ijn,displs_g,mpi_rtype,&
mype_out,mpi_comm_world,ierror)
if (mype == mype_out) then
call load_grid(work1,grid)
ij=izero
do j=1,nlatm2
do i=1,nlon
ij=ij+ione
gfsdata%q(ij,k,2)=grid(i,j)
end do
end do
call gfsio_writerecv(gfileo,'o3mr','layer',k,gfsdata%q(:,k,2),iret,&
idrt=gfshead%idrt)
endif
end do
endif
! Cloud condensate mixing ratio
if (ntracer>2_i_kind .or. ncloud>=ione) then
if (ncep_sigio) then
if (mype==mype_clc) then
!!$omp parallel do private(k,grid,i,spec_work,grid2,spec_work_sm)
do k=1,nsig
call load_grid(work1_k(1,k),grid)
do i=1,sp_b%nc
spec_work(i) = temp(i,k)
if(sp_b%factsml(i))spec_work(i)=zero
end do
call general_sptez_s_b(sp_a,sp_b,spec_work,grid2,ione)
call load_grid(work1_k(1,k),grid)
grid=grid-grid2
call general_sptez_s(sp_a,spec_work_sm,grid,-ione)
call sppad(izero,sp_a%jcap,spec_work_sm,izero,sp_b%jcap,spec_work)
do i=1,sp_b%nc
temp(i,k) =temp(i,k)+spec_work(i)
if(sp_b%factsml(i))temp(i,k)=zero_single
end do
end do
!!$omp end parallel do
!send sh back to mype_out
call mpi_send(temp,sp_b%nc*nsig,mpi_rtype4,mype_out,&
itag_clc,mpi_comm_world,ierror)
endif
else !GFS I/O
do k=1,nsig
call mpi_gatherv(qsm(1,k,3),ijn(mm1),mpi_rtype,&
work1,ijn,displs_g,mpi_rtype,&
mype_out,mpi_comm_world,ierror)
if (mype == mype_out) then
ij=izero
do j=1,nlatm2
do i=1,nlon
ij=ij+ione
gfsdata%q(ij,k,3)=grid(i,j)
end do
end do
call gfsio_writerecv(gfileo,'clwmr','layer',k,gfsdata%q(:,k,3),iret,&
idrt=gfshead%idrt)
endif
end do
endif
endif
! NCEP_SIGIO specific output
if (ncep_sigio) then
! Horizontal divergence and voriticy
if (mype==mype_div) then
!!$omp parallel do private(k,grid,i,spec_work,grid2,spec_work_sm)
do k=1,nsig
do i=1,sp_b%nc
spec_work(i) = temp(i,k)
if(sp_b%factsml(i))spec_work(i)=zero
end do
call general_sptez_s_b(sp_a,sp_b,spec_work,grid2,ione)
call load_grid(work1_k(1,k),grid)
grid=grid-grid2
call general_sptez_s(sp_a,spec_work_sm,grid,-ione)
call sppad(izero,sp_a%jcap,spec_work_sm,izero,sp_b%jcap,spec_work)
do i=1,sp_b%nc
temp(i,k) = temp(i,k) + spec_work(i)
if(sp_b%factvml(i))temp(i,k)=zero_single
end do
end do
!!$omp end parallel do
!send sh back to mype_out
call mpi_send(temp,sp_b%nc*nsig,mpi_rtype4,mype_out,&
itag_div,mpi_comm_world,ierror)
endif
if (mype==mype_vort) then
!!$omp parallel do private(k,grid,i,spec_work,grid2,spec_work_sm)
do k=1,nsig
do i=1,sp_b%nc
spec_work(i) = temp(i,k)
if(sp_b%factsml(i))spec_work(i)=zero
end do
call general_sptez_s_b(sp_a,sp_b,spec_work,grid2,ione)
call load_grid(work1_k(1,k),grid)
grid=grid-grid2
call general_sptez_s(sp_a,spec_work_sm,grid,-ione)
call sppad(izero,sp_a%jcap,spec_work_sm,izero,sp_b%jcap,spec_work)
do i=1,sp_b%nc
temp(i,k) =temp(i,k) + spec_work(i)
if(sp_b%factvml(i))temp(i,k)=zero_single
end do
end do
!!$omp end parallel do
call mpi_send(temp,sp_b%nc*nsig,mpi_rtype4,mype_out,&
itag_vort,mpi_comm_world,ierror)
endif
endif
! GFSIO specific output
if (.not.ncep_sigio) then
! Pressure depth
do k=1,nsig
call mpi_gatherv(dpsm(1,k),ijn(mm1),mpi_rtype,&
work1,ijn,displs_g,mpi_rtype,&
mype_out,mpi_comm_world,ierror)
if (mype==mype_out) then
call load_grid(work1,grid)
ij=izero
do j=1,nlatm2
do i=1,nlon
ij=ij+ione
gfsdata%dp(ij,k)=grid(i,j)*r1000
end do
end do
call gfsio_writerecv(gfileo,'dpres','layer',k,gfsdata%dp(:,k),iret,&
idrt=gfshead%idrt)
endif
end do
! Layer mean pressure
do k=1,nsig
call mpi_gatherv(prslm(1,k),ijn(mm1),mpi_rtype,&
work1,ijn,displs_g,mpi_rtype,&
mype_out,mpi_comm_world,ierror)
if (mype==mype_out) then
call load_grid(work1,grid)
ij=izero
do j=1,nlatm2
do i=1,nlon
ij=ij+ione
gfsdata%p(ij,k)=grid(i,j)*r1000
end do
end do
call gfsio_writerecv(gfileo,'pres','layer',k,gfsdata%p(:,k),iret,&
idrt=gfshead%idrt)
endif
end do
! Zonal wind
do k=1,nsig
call mpi_gatherv(usm(1,k),ijn(mm1),mpi_rtype,&
work1,ijn,displs_g,mpi_rtype,&
mype_out,mpi_comm_world,ierror)
if (mype==mype_out) then
call load_grid(work1,grid)
ij=izero
do j=1,nlatm2
do i=1,nlon
ij=ij+ione
gfsdata%u(ij,k)=grid(i,j)
end do
end do
call gfsio_writerecv(gfileo,'ugrd','layer',k,gfsdata%u(:,k),iret,&
idrt=gfshead%idrt)
endif
end do
! Meridional wind
do k=1,nsig
call mpi_gatherv(vsm(1,k),ijn(mm1),mpi_rtype,&
work1,ijn,displs_g,mpi_rtype,&
mype_out,mpi_comm_world,ierror)
if (mype==mype_out) then
call load_grid(work1,grid)
ij=izero
do j=1,nlatm2
do i=1,nlon
ij=ij+ione
gfsdata%v(ij,k)=grid(i,j)
end do
end do
call gfsio_writerecv(gfileo,'vgrd','layer',k,gfsdata%v(:,k),iret,&
idrt=gfshead%idrt)
endif
end do
endif
! Single task writes analysis data to analysis file
if (mype==mype_out) then
if (ncep_sigio) then
!receive temperature from mype_th
call mpi_recv(sigdata%t,sp_b%nc*nsig,mpi_rtype4,mype_th,&
itag_th,mpi_comm_world,status,ierror)
!receive specific humidity from mype_sh
call mpi_recv(sigdata%q(1,1,1),sp_b%nc*nsig,mpi_rtype4,mype_sh,&
itag_sh,mpi_comm_world,status,ierror)
!receive ozone from mype_oz
call mpi_recv(sigdata%q(1,1,2),sp_b%nc*nsig,mpi_rtype4,mype_oz,&
itag_oz,mpi_comm_world,status,ierror)
!receive cloud condensate mixing ratio from mype_clc
if (ntracer>2_i_kind .or. ncloud>=ione) then
call mpi_recv(sigdata%q(1,1,3),sp_b%nc*nsig,mpi_rtype4,mype_clc,&
itag_clc,mpi_comm_world,status,ierror)
endif
!receive divergence from mype_div
call mpi_recv(sigdata%d,sp_b%nc*nsig,mpi_rtype4,mype_div,&
itag_div,mpi_comm_world,status,ierror)
!receive vorticity from mype_vort
call mpi_recv(sigdata%z,sp_b%nc*nsig,mpi_rtype4,mype_vort,&
itag_vort,mpi_comm_world,status,ierror)
!
string='sigio'
call sigio_swdata(lunanl,sighead,sigdata,iret)
call sigio_axdata(sigdata,iret)
else !GFS IO
string='gfsio'
call gfsio_close(gfileo,iret)
endif
write(6,110) string,gfshead%jcap,gfshead%latb,gfshead%lonb,&
gfshead%levs,gfshead%fhour,gfshead%idate
110 format('WRITE_GFSATM: NCEP ',a5,&
' atm anal written for jcap,latb,lonb,levs= ',4i6,&
' valid hour,idate= ',f3.1,4(i4,1x))
endif
return
end subroutine write_gfsatm
subroutine write_gfssfc(filename,mype,mype_sfc,dsfct)
!$$$ subprogram documentation block
! . . .
! subprogram: write_gfssfc --- Write surface analysis to file
!
! prgrmmr: treadon - initial version; org: np22
!
! abstract: This routine writes the updated surface analysis. At
! this point (20040615) the only surface field update by
! the gsi is the skin temperature. The current (20040615)
! GDAS setup does use the updated surface file. Rather,
! the output from surface cycle is used as the surface
! analysis for subsequent GFS runs.
!
! The routine gathers surface fields from subdomains,
! reformats the data records, and then writes each record
! to the output file.
!
! Since the gsi only update the skin temperature, all
! other surface fields are simply read from the guess
! surface file and written to the analysis file.
!
! Structure of GFS surface file
! data record 1 label
! data record 2 date, dimension, version, lons/lat record
! data record 3 tsf
! data record 4 soilm(two layers)
! data record 5 snow
! data record 6 soilt(two layers)
! data record 7 tg3
! data record 8 zor
! data record 9 cv
! data record 10 cvb
! data record 11 cvt
! data record 12 albedo (four types)
! data record 13 slimsk
! data record 14 vegetation cover
! data record 15 plantr
! data record 16 f10m
! data record 17 canopy water content (cnpanl)
! data record 18 vegetation type
! data record 19 soil type
! data record 20 zenith angle dependent vegetation fraction (two types)
!
! program history log:
! 2004-06-15 treadon - updated documentation
! 2004-07-15 todling - protex-compliant prologue; added intent/only's
! 2004-12-03 treadon - replace mpe_igatherv (IBM extension) with
! standard mpi_gatherv
! 2005-01-27 treadon - rewrite to make use of sfcio module
! 2005-02-09 kleist - clean up unit number and filename for updated surface file
! 2005-03-07 todling - die gracefully when return error from sfcio
! 2005-03-10 treadon - remove iadate from calling list, access via obsmod
! 2006-10-11 treadon - update 10m wind factor in sfc file
! 2008-05-28 safford - rm unused vars
!
! input argument list:
! filename - file to open and write to
! dsfct - delta skin temperature
! mype - mpi task number
! mype_sfc - mpi task to write output file
!
! output argument list:
!
! attributes:
! language: f90
! machines: ibm RS/6000 SP; SGI Origin 2000; Compaq HP
!
!$$$ end documentation block
! !USES:
use kinds, only: r_kind,r_single,i_kind
use mpimod, only: mpi_rtype
use mpimod, only: mpi_comm_world
use mpimod, only: ierror
use gridmod, only: nlat,nlon
use gridmod, only: lat1,lon1
use gridmod, only: lat2,lon2
use gridmod, only: iglobal
use gridmod, only: ijn
use gridmod, only: ltosi,ltosj
use gridmod, only: displs_g
use gridmod, only: itotsub
use obsmod, only: iadate
use constants, only: izero,ione,zero_single
use sfcio_module, only: sfcio_intkind,sfcio_head,sfcio_data,&
sfcio_srohdc,sfcio_swohdc,sfcio_axdata
implicit none
! !INPUT PARAMETERS:
character(24) ,intent(in ) :: filename ! file to open and write to
real(r_kind),dimension(lat2,lon2),intent(in ) :: dsfct ! delta skin temperature
integer(i_kind) ,intent(in ) :: mype ! mpi task number
integer(i_kind) ,intent(in ) :: mype_sfc ! mpi task to write output file
! !OUTPUT PARAMETERS:
!-------------------------------------------------------------------------
! Declare local parameters
character( 6),parameter:: fname_ges='sfcf06'
integer(sfcio_intkind),parameter:: ioges = 12
integer(sfcio_intkind),parameter:: ioanl = 52
real(r_kind),parameter :: houra = zero_single
! Declare local variables
integer(sfcio_intkind):: iret
integer(i_kind) latb,lonb,nlatm2
integer(i_kind) i,j,ip1,jp1,ilat,ilon,jj,mm1
real(r_kind),dimension(nlon,nlat):: buffer
real(r_kind),dimension(lat1,lon1):: sfcsub
real(r_kind),dimension(nlon,nlat):: grid
real(r_kind),dimension(max(iglobal,itotsub)):: sfcall
real(r_single),allocatable,dimension(:,:):: buffer2
type(sfcio_head):: head
type(sfcio_data):: data
!*****************************************************************************
! Initialize local variables
mm1=mype+ione
nlatm2=nlat-2_i_kind
! Gather skin temperature information from all tasks.
do j=1,lon1
jp1 = j+ione
do i=1,lat1
ip1 = i+ione
sfcsub(i,j)=dsfct(ip1,jp1)
end do
end do
call mpi_gatherv(sfcsub,ijn(mm1),mpi_rtype,&
sfcall,ijn,displs_g,mpi_rtype,mype_sfc,&
mpi_comm_world,ierror)
! Only MPI task mype_sfc writes the surface file.
if (mype==mype_sfc) then
! Reorder updated skin temperature to output format
do i=1,iglobal
ilon=ltosj(i)
ilat=ltosi(i)
grid(ilon,ilat)=sfcall(i)
end do
do j=1,nlat
jj=nlat-j+ione
do i=1,nlon
buffer(i,j)=grid(i,jj)
end do
end do
! For now, rather than carry around all the surface fields in memory from
! the read in ingesfc, just read fields from surface file. Also, for
! now, only update the 6-hour forecast surface guess file.
! Read surface guess file
call sfcio_srohdc(ioges,fname_ges,head,data,iret)
if (iret /= izero) then
write(6,*)'WRITE_GFSSFC: ***ERROR*** problem reading ',fname_ges,&
', iret=',iret
call sfcio_axdata(data,iret)
call stop2(80)
endif
latb=head%latb
lonb=head%lonb
allocate(buffer2(lonb,latb))
if ( (latb /= nlatm2) .or. &
(lonb /= nlon) ) then
write(6,*)'WRITE_GFSSFC: different grid dimensions analysis vs sfc. interpolating sfc temperature ',&
', nlon,nlat-2=',nlon,nlatm2,' -vs- sfc file lonb,latb=',&
lonb,latb
call sfc_interpolate(buffer,nlon,nlat,buffer2,lonb,latb)
else
do j=1,latb
do i=1,lonb
buffer2(i,j)=buffer(i,j+ione)
end do
end do
endif
! Update guess date/time to analysis date/time
head%fhour = houra ! forecast hour
head%idate(1)=iadate(4) ! hour
head%idate(2)=iadate(2) ! month
head%idate(3)=iadate(3) ! day
head%idate(4)=iadate(1) ! year
do j=1,latb
do i=1,lonb
data%tsea(i,j) = data%tsea(i,j)+buffer2(i,j)
end do
end do
deallocate(buffer2)
! Write updated information to surface analysis file
call sfcio_swohdc(ioanl,filename,head,data,iret)
! Deallocate local work arrays
call sfcio_axdata(data,iret)
write(6,100) lonb,latb,houra,iadate(1:4),iret
100 format(' WRITE_GFSSFC: sfc analysis written for ',&
2i6,1x,f3.1,4(i4,1x),' with iret=',i2)
endif
! End of routine
return
end subroutine write_gfssfc
subroutine reorder_gfsgrib(nx,ny,grid_1d,grid_2d)
!$$$ subprogram documentation block
! . . .
! subprogram: reorder_gfsgrib --- transfer gfsgrib data 1d <--> 2d
!
! prgrmmr: treadon - initial version; org: np23
!
! abstract: This routine transfers the contents of gfs grib arrays
! between 1d and 2d.
!
! program history log:
! 2007-04-30 treadon - original routine
! 2008-05-28 safford -- add subprogram doc block
!
! input argument list:
! nx - number of grid points in zonal direction
! ny - number of grid points in meridional direction
! grid_1d - 1d array
!
! output argument list:
! grid_2d - 2d array
!
! attributes:
! language: f90
! machines: ibm RS/6000 SP; SGI Origin 2000; Compaq HP
!
!$$$ end documentation block
! !USES:
use kinds, only: r_kind,i_kind
use constants, only: izero,ione
implicit none
! !INPUT PARAMETERS:
integer(i_kind) ,intent(in ) :: nx ! number of grid points in zonal direction
integer(i_kind) ,intent(in ) :: ny ! number of grid points in meridional direction
real(r_kind),dimension(nx*ny),intent(in ) :: grid_1d ! 1d array
! !OUTPUT PARAMETERS:
real(r_kind),dimension(nx,ny),intent( out) :: grid_2d ! 2d array
!-------------------------------------------------------------------------
! Declare local variables
integer(i_kind) i,j,ij
!*****************************************************************************
! Loop to transfer array contents
ij=izero
do j=1,ny
do i=1,nx
ij=ij+ione
grid_2d(i,j)=grid_1d(ij)
end do
end do
! End of routine
return
end subroutine reorder_gfsgrib
subroutine sfc_interpolate(a,na_lon,na_lat,b,ns_lon,ns_lat)
!$$$ subprogram documentation block
! . . .
! subprogram: sfc_interpolate --- interpolates from analysis grid to
! surface grid
! prgrmmr: derber - initial version; org: np2
!
! abstract: This routine interpolates a on analysis grid to b on
! surface grid
!
! program history log:
! 2008-02-26 derber - original routine
! 2008-05-28 safford - add subprogram doc block, rm unused uses
!
! input argument list:
! na_lon - number of longitude grid analysis
! na_lat - number of latitude grid analysis
! ns_lon - number of longitude grid sfc
! ns_lat - number of latitude grid sfc
! a - analysis values
!
! output argument list:
! b - surface values
!
! attributes:
! language: f90
! machines: ibm RS/6000 SP; SGI Origin 2000; Compaq HP
!
!$$$ end documentation block
! !USES:
use kinds, only: r_kind,i_kind,r_single
use constants, only: izero,ione,zero,one
use gridmod, only: rlats,rlons,rlats_sfc,rlons_sfc
implicit none
! !INPUT PARAMETERS:
integer(i_kind) ,intent(in ) :: na_lon ! number of longitude grid analysis
integer(i_kind) ,intent(in ) :: na_lat ! number of latitude grid analysis
integer(i_kind) ,intent(in ) :: ns_lon ! number of longitude grid sfc
integer(i_kind) ,intent(in ) :: ns_lat ! number of latitude grid sfc
real(r_kind) ,dimension(na_lon,na_lat),intent(in ) :: a ! analysis values
! !OUTPUT PARAMETERS:
real(r_single),dimension(ns_lon,ns_lat),intent( out) :: b ! surface values
! Declare local variables
integer(i_kind) i,j,ix,iy,ixp,iyp
real(r_kind) dx1,dy1,dx,dy,w00,w01,w10,w11,bout,dlat,dlon
!*****************************************************************************
b=zero
! Loop over all points to get interpolated value
do j=1,ns_lat
dlat=rlats_sfc(j)
call grdcrd(dlat,ione,rlats,na_lat,ione)
iy=int(dlat)
iy=min(max(ione,iy),na_lat)
dy =dlat-iy
dy1 =one-dy
iyp=min(na_lat,iy+ione)
do i=1,ns_lon
dlon=rlons_sfc(i)
call grdcrd(dlon,ione,rlons,na_lon,ione)
ix=int(dlon)
dx =dlon-ix
dx=max(zero,min(dx,one))
dx1 =one-dx
w00=dx1*dy1; w10=dx1*dy; w01=dx*dy1; w11=dx*dy
ix=min(max(izero,ix),na_lon)
ixp=ix+ione
if(ix==izero) ix=na_lon
if(ixp==na_lon+ione) ixp=ione
bout=w00*a(ix,iy)+w01*a(ix,iyp)+w10*a(ixp,iy)+w11*a(ixp,iyp)
b(i,j)=bout
end do
end do
! End of routine
return
end subroutine sfc_interpolate
!-------------------------------------------------------------------------------
subroutine sigio_cnvtdv8(im,ix,km,idvc,idvm,ntrac,iret,t,q,cpi,cnflg)
!$$$ subprogram documentation block
! . . .
! subprogram: sigio_cnvtdv8
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-05-28 safford -- add subprogram doc block
!
! input argument list:
! im,ix,km,idvc,idvm,ntrac,cnflg
! q, cpi
! t
!
! output argument list:
! iret
! t
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind,r_kind
use constants, only: izero,zero,one,fv
implicit none
integer(i_kind),intent(in ) :: im,ix,km,idvc,idvm,ntrac,cnflg
integer(i_kind),intent( out) :: iret
real(r_kind) ,intent(in ) :: q(ix,km,ntrac), cpi(0:ntrac)
real(r_kind) ,intent(inout) :: t(ix,km)
integer(i_kind) :: thermodyn_id, n
real(r_kind) :: xcp(ix,km), sumq(ix,km)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
iret=izero
thermodyn_id = mod(IDVM/10,10_i_kind)
!
if (thermodyn_id == 3_i_kind .and. idvc == 3_i_kind) then
xcp(1:im,:) = zero
sumq(1:im,:) = zero
do n=1,NTRAC
if( cpi(n) /= zero) then
xcp(1:im,:) = xcp(1:im,:) + q(1:im,:,n) * cpi(n)
sumq(1:im,:) = sumq(1:im,:) + q(1:im,:,n)
endif
enddo
xcp(1:im,:) = (one-sumq(1:im,:))*cpi(0) + xcp(1:im,:) ! Mean Cp
!
else
xcp(1:im,:) = one + fv*Q(1:im,:,1) ! Virt factor
endif
if (cnflg > izero) then
t(1:im,:) = t(1:im,:) / xcp(1:im,:)
else
t(1:im,:) = t(1:im,:) * xcp(1:im,:)
endif
!
return
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
end subroutine sigio_cnvtdv8
end module ncepgfs_io