subroutine general_read_nemsaero(grd,sp_a,filename,mype,gfschem_bundle, &
naero,aeroname,init_head,iret_read)
!$$$ subprogram documentation block
! . . . .
! subprogram: general_read_nemsaero adaptation of general_read_gfsatm
! for reading in aerosols from NEMSI/O
!
! abstract: copied from general_read_gfsatm, primarily for reading in aerosol
! tracer variables from NEMS GFS I/O files
!
! program history log:
! 2019-04-19 Wei/Martin - copied and modified to read in aerosol arrays
! from either FV3-Chem or NEMS
!
! input argument list:
! grd - structure variable containing information about grid
! (initialized by general_sub2grid_create_info,
! located in general_sub2grid_mod.f90)
! sp_a - structure variable containing spectral information for
! analysis
! (initialized by general_init_spec_vars, located in
! general_specmod.f90)
! filename - input sigma file name
! mype - mpi task id
! naero - number of aerosol tracers to read
! aeroname - len(naero) character strings of aerosol tracers to read
! init_head- flag to read header record. Usually .true. unless
! repeatedly
! reading similar files (ensembles)
!
! input/output list:
! gfschem_bundle - GSI bundle containing chem/aerosol arrays
!
! output argument list:
! iret_read - return code, 0 for successful read.
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,r_single,i_kind
use gridmod, only: use_fv3_aero
use general_commvars_mod, only: fill_ns,fill2_ns
use general_sub2grid_mod, only: sub2grid_info
use general_specmod, only: spec_vars
use mpimod, only: npe
use constants, only: zero,one,r0_01
use nemsio_module, only: nemsio_init,nemsio_open,nemsio_close
use ncepnems_io, only: error_msg
use nemsio_module, only: nemsio_gfile,nemsio_getfilehead,nemsio_readrecv
use egrid2agrid_mod, only: g_egrid2agrid,g_create_egrid2agrid,egrid2agrid_parm,destroy_egrid2agrid
use constants, only: two,pi,half,deg2rad,r60,r3600
use gsi_bundlemod, only: gsi_bundle, gsi_bundlegetpointer
implicit none
! Declare local parameters
real(r_kind),parameter:: r0_001 = 0.001_r_kind
! Declare passed variables
type(sub2grid_info) ,intent(in ) :: grd
type(spec_vars) ,intent(in ) :: sp_a
character(*) ,intent(in ) :: filename
integer(i_kind) ,intent(in ) :: mype
integer(i_kind) ,intent(in ) :: naero
character(*),dimension(naero) ,intent(in ) :: aeroname
logical ,intent(in ) :: init_head
integer(i_kind) ,intent( out) :: iret_read
type(gsi_bundle) ,intent(inout) :: gfschem_bundle
! Declare local variables
character(len=120) :: my_name = 'general_read_nemsaero'
character(len=1) :: null = ' '
character(len=20),dimension(npe) :: ch_aero
integer(i_kind):: iret,nlatm2,nlevs,icm,nord_int
integer(i_kind):: i,j,k,l,icount,kk,istatus,ier
integer(i_kind) :: latb, lonb, levs, nframe
integer(i_kind) :: nfhour, nfminute, nfsecondn, nfsecondd
integer(i_kind) :: istop = 101
integer(i_kind),dimension(npe)::ilev,iflag,mype_use
integer(i_kind),dimension(7):: idate
integer(i_kind),dimension(4):: odate
real(r_kind) :: fhour
real(r_kind),dimension(:,:,:),pointer :: &
ae_d1,ae_d2,ae_d3,ae_d4,ae_d5,&
ae_s1,ae_s2,ae_s3,ae_s4,ae_so4,&
ae_ocpho,ae_ocphi,ae_bcpho,ae_bcphi
real(r_kind),allocatable,dimension(:,:) :: grid, grid_v,grid_b
real(r_kind),allocatable,dimension(:,:,:) :: grid_c, grid2
real(r_kind),allocatable,dimension(:) :: work
real(r_kind),allocatable,dimension(:) :: rwork1d0, rwork1d1, rwork1d2
real(r_kind),allocatable,dimension(:) :: rlats,rlons,clons,slons
real(r_single),allocatable,dimension(:) :: r4lats,r4lons
logical :: procuse,diff_res,eqspace
type(nemsio_gfile) :: gfile
type(egrid2agrid_parm) :: p_high
logical,dimension(1) :: vector
!******************************************************************************
if(mype==0) write(6,*) trim(my_name)," start and filename is ",trim(filename)
! Initialize variables used below
iret_read=0
iret=0
nlatm2=grd%nlat-2
iflag = 0
ilev = 0
nlevs=grd%nsig
mype_use=-1
icount=0
procuse=.false.
if(mype == 0)procuse = .true.
do i=1,npe
if(grd%recvcounts_s(i-1) > 0)then
icount = icount+1
mype_use(icount)=i-1
if(i-1 == mype) procuse=.true.
end if
end do
icm=icount
allocate( work(grd%itotsub))
work=zero
if(procuse)then
if (init_head) call nemsio_init(iret=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'init',istop,iret)
call nemsio_open(gfile,filename,'READ',iret=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'open',istop+1,iret)
call nemsio_getfilehead(gfile,iret=iret, nframe=nframe, &
nfhour=nfhour, nfminute=nfminute, nfsecondn=nfsecondn, nfsecondd=nfsecondd, &
idate=idate, dimx=lonb, dimy=latb,dimz=levs)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'getfilehead',istop+1,iret)
if( nframe /= 0 ) then
if ( mype == 0 ) &
write(6,*)trim(my_name),': ***ERROR*** nframe /= 0 for global model read, nframe = ', nframe
call stop2(101)
end if
fhour = float(nfhour) + float(nfminute)/r60 + float(nfsecondn)/float(nfsecondd)/r3600
odate(1) = idate(4) !hour
odate(2) = idate(2) !month
odate(3) = idate(3) !day
odate(4) = idate(1) !year
if ( iret == 0 .and. mype == 0 ) then
write(6,'(''Aerosol file time='',i4.4,i2.2,i2.2,i2.2)') odate(4),odate(2),odate(3),odate(1)
end if
!
! g_* array already pre-allocate as (lat2,lon2,<nsig>) => 2D and <3D>
! array
!
diff_res=.false.
if(latb /= nlatm2) then
diff_res=.true.
if ( mype == 0 ) write(6, &
'(a,'': different spatial dimension nlatm2 = '',i4,tr1,''latb = '',i4)') &
trim(my_name),nlatm2,latb
end if
if(lonb /= grd%nlon) then
diff_res=.true.
if ( mype == 0 ) write(6, &
'(a,'': different spatial dimension nlon = '',i4,tr1,''lonb = '',i4)') &
trim(my_name),grd%nlon,lonb
end if
if(levs /= grd%nsig)then
if ( mype == 0 ) write(6, &
'(a,'': inconsistent spatial dimension nsig = '',i4,tr1,''levs = '',i4)') &
trim(my_name),grd%nsig,levs
call stop2(101)
end if
allocate( grid(grd%nlon,nlatm2), grid_v(grd%nlon,nlatm2) )
if(diff_res)then
allocate( grid_b(lonb,latb),grid_c(latb+2,lonb,1),grid2(grd%nlat,grd%nlon,1))
end if
allocate( rwork1d0(latb*lonb) )
allocate(rlats(latb+2),rlons(lonb),clons(lonb),slons(lonb),r4lats(lonb*latb),r4lons(lonb*latb))
allocate(rwork1d1(latb*lonb),rwork1d2(latb*lonb))
call nemsio_getfilehead(gfile,lat=r4lats,iret=iret)
call nemsio_getfilehead(gfile,lon=r4lons,iret=iret)
do j=1,latb
rlats(latb+2-j)=deg2rad*r4lats(lonb/2+(j-1)*lonb)
end do
do j=1,lonb
rlons(j)=deg2rad*r4lons(j)
end do
deallocate(r4lats,r4lons)
rlats(1)=-half*pi
rlats(latb+2)=half*pi
do j=1,lonb
clons(j)=cos(rlons(j))
slons(j)=sin(rlons(j))
end do
nord_int=4
eqspace=.false.
call g_create_egrid2agrid(grd%nlat,sp_a%rlats,grd%nlon,sp_a%rlons,&
latb+2,rlats,lonb,rlons,&
nord_int,p_high,.true.,eqspace=eqspace)
deallocate(rlats,rlons)
end if
istatus=0
do l=1,naero
select case(trim(aeroname(l)))
case ('sulf')
call gsi_bundlegetpointer(gfschem_bundle,'sulf' ,ae_so4 ,ier); istatus=istatus+ier
case ('oc1')
call gsi_bundlegetpointer(gfschem_bundle,'oc1' ,ae_ocpho,ier); istatus=istatus+ier
case ('oc2')
call gsi_bundlegetpointer(gfschem_bundle,'oc2' ,ae_ocphi,ier); istatus=istatus+ier
case ('bc1')
call gsi_bundlegetpointer(gfschem_bundle,'bc1' ,ae_bcpho,ier); istatus=istatus+ier
case ('bc2')
call gsi_bundlegetpointer(gfschem_bundle,'bc2' ,ae_bcphi,ier); istatus=istatus+ier
case ('dust1')
call gsi_bundlegetpointer(gfschem_bundle,'dust1',ae_d1 ,ier); istatus=istatus+ier
case ('dust2')
call gsi_bundlegetpointer(gfschem_bundle,'dust2',ae_d2 ,ier); istatus=istatus+ier
case ('dust3')
call gsi_bundlegetpointer(gfschem_bundle,'dust3',ae_d3 ,ier); istatus=istatus+ier
case ('dust4')
call gsi_bundlegetpointer(gfschem_bundle,'dust4',ae_d4 ,ier); istatus=istatus+ier
case ('dust5')
call gsi_bundlegetpointer(gfschem_bundle,'dust5',ae_d5 ,ier); istatus=istatus+ier
case ('seas1')
call gsi_bundlegetpointer(gfschem_bundle,'seas1',ae_s1 ,ier); istatus=istatus+ier
case ('seas2')
call gsi_bundlegetpointer(gfschem_bundle,'seas2',ae_s2 ,ier); istatus=istatus+ier
case ('seas3')
call gsi_bundlegetpointer(gfschem_bundle,'seas3',ae_s3 ,ier); istatus=istatus+ier
case ('seas4')
call gsi_bundlegetpointer(gfschem_bundle,'seas4',ae_s4 ,ier); istatus=istatus+ier
end select
end do
if ( istatus /= 0 ) then
if ( mype == 0 ) then
write(6,*) 'general_read_nemsaero: ERROR'
write(6,*) 'Missing some of the required fields'
write(6,*) 'Aborting ... '
endif
call stop2(999)
endif
icount=0
! Process guess fields according to type of input file. NCEP_SIGIO
! files
! are spectral coefficient files and need to be transformed to the
! grid.
! Once on the grid, fields need to be scattered from the full domain
! to
! sub-domains.
do l=1,naero
do k=1,nlevs
icount=icount+1
ilev(icount)=k
ch_aero(icount)=trim(aeroname(l))
vector(1)=.false.
if (mype==mype_use(icount)) then
if (use_fv3_aero) then
! variable names in FV3GFS-GSDChem
if ( aeroname(l)(1:4) == 'seas') then
select case ( trim(aeroname(l)) )
case ('seas1')
call nemsio_readrecv(gfile,'seas1','mid layer',k,rwork1d1,iret=iret)
call nemsio_readrecv(gfile,'seas2','mid layer',k,rwork1d2,iret=iret)
rwork1d0=rwork1d1+rwork1d2
case ('seas2')
call nemsio_readrecv(gfile,'seas3','mid layer',k,rwork1d0,iret=iret)
case ('seas3')
call nemsio_readrecv(gfile,'seas4','mid layer',k,rwork1d0,iret=iret)
case ('seas4')
call nemsio_readrecv(gfile,'seas5','mid layer',k,rwork1d0,iret=iret)
end select
else
! many of the names are the same in the GSI bundle as they are
! in the FV3GFS-GSDChem NEMSIO files
call nemsio_readrecv(gfile,trim(aeroname(l)),'mid layer',k,rwork1d0,iret=iret)
end if
else
! variable names in NGACv2
select case ( trim(aeroname(l)) )
case ('sulf')
call nemsio_readrecv(gfile,'so4','mid layer',k,rwork1d0,iret=iret)
case ('bc1')
call nemsio_readrecv(gfile,'bcphobic','mid layer',k,rwork1d0,iret=iret)
case ('bc2')
call nemsio_readrecv(gfile,'bcphilic','mid layer',k,rwork1d0,iret=iret)
case ('oc1')
call nemsio_readrecv(gfile,'ocphobic','mid layer',k,rwork1d0,iret=iret)
case ('oc2')
call nemsio_readrecv(gfile,'ocphilic','mid layer',k,rwork1d0,iret=iret)
case ('dust1')
call nemsio_readrecv(gfile,'du001','mid layer',k,rwork1d0,iret=iret)
case ('dust2')
call nemsio_readrecv(gfile,'du002','mid layer',k,rwork1d0,iret=iret)
case ('dust3')
call nemsio_readrecv(gfile,'du003','mid layer',k,rwork1d0,iret=iret)
case ('dust4')
call nemsio_readrecv(gfile,'du004','mid layer',k,rwork1d0,iret=iret)
case ('dust5')
call nemsio_readrecv(gfile,'du005','mid layer',k,rwork1d0,iret=iret)
case ('seas1')
call nemsio_readrecv(gfile,'ss001','mid layer',k,rwork1d1,iret=iret)
call nemsio_readrecv(gfile,'ss002','mid layer',k,rwork1d2,iret=iret)
rwork1d0=rwork1d1+rwork1d2
case ('seas2')
call nemsio_readrecv(gfile,'ss003','mid layer',k,rwork1d0,iret=iret)
case ('seas3')
call nemsio_readrecv(gfile,'ss004','mid layer',k,rwork1d0,iret=iret)
case ('seas4')
call nemsio_readrecv(gfile,'ss005','mid layer',k,rwork1d0,iret=iret)
end select
! Convert NGAC mixing ratio unit from kg/kg( 10^3 g/kg ) to ug/kg( 10^-6 g/kg )
rwork1d0=rwork1d0*1.0e+9_r_kind
end if ! NGAC vs FV3-Chem
if (iret /= 0) call error_msg(trim(my_name),trim(filename),'tmp','read',istop+7,iret)
if(diff_res)then
grid_b=reshape(rwork1d0,(/size(grid_b,1),size(grid_b,2)/))
call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector)
do kk=1,grd%itotsub
i=grd%ltosi_s(kk)
j=grd%ltosj_s(kk)
work(kk)=grid2(i,j,1)
end do
else
grid=reshape(rwork1d0,(/size(grid,1),size(grid,2)/))
call general_fill_ns(grd,grid,work)
end if
end if
if(icount == icm)then
call aerosol_reload(grd,ae_d1,ae_d2,ae_d3,ae_d4,ae_d5, &
ae_s1,ae_s2,ae_s3,ae_s4,ae_so4,&
ae_ocpho,ae_ocphi,ae_bcpho,ae_bcphi, &
icount,ilev,ch_aero,work)
end if
end do
end do
if(procuse)then
if(diff_res) deallocate(grid_b,grid_c,grid2)
call destroy_egrid2agrid(p_high)
deallocate(rwork1d1,clons,slons)
deallocate(rwork1d0)
deallocate(grid)
call nemsio_close(gfile,iret=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'close',istop+9,iret)
end if
deallocate(work)
! Print date/time stamp
if(mype==0) then
write(6,700) lonb,latb,nlevs,grd%nlon,nlatm2,&
fhour,odate
700 format('READ_GLOBAL_AEROSOL: ges read/scatter, lonb,latb,levs=',&
3i6,', nlon,nlat=',2i6,', hour=',f10.1,', idate=',4i5)
end if
return
! ERROR detected while reading file
1000 continue
write(6,*)'GENERAL_READ_GFSATM: ***ERROR*** reading ',&
trim(filename),' mype,iret_read=',mype,iret_read,grd%nsig,nlevs
return
! End of routine. Return
return
end subroutine general_read_nemsaero
!
subroutine aerosol_reload(grd,ae_d1,ae_d2,ae_d3,ae_d4,ae_d5, &
ae_s1,ae_s2,ae_s3,ae_s4,ae_so4, &
ae_ocpho,ae_ocphi,ae_bcpho,ae_bcphi, &
icount,ilev,chaero,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
implicit none
! !INPUT PARAMETERS:
type(sub2grid_info) ,intent(in ) :: grd
integer(i_kind),intent(inout) ::icount
integer(i_kind),dimension(npe),intent(inout):: ilev!,iflag
real(r_kind),dimension(grd%itotsub),intent(in) :: work
character(*),dimension(npe) , intent(in) :: chaero
! !OUTPUT PARAMETERS:
real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig),intent( out) :: &
ae_d1,ae_d2,ae_d3,ae_d4,ae_d5, &
ae_s1,ae_s2,ae_s3,ae_s4,ae_so4, &
ae_ocpho,ae_ocphi,ae_bcpho,ae_bcphi
! !DESCRIPTION: Transfer contents of 2-d array to 3-d array
!
! !REVISION HISTORY:
! 2004-05-14 treadon
! 2004-07-15 todling, protex-compliant prologue
!
! !REMARKS:
!
! language: f90
! machine: ibm rs/6000 sp; sgi origin 2000; compaq/hp
!
! !AUTHOR:
! treadon org: np23 date: 2004-05-14
!
!EOP
!-------------------------------------------------------------------------
integer(i_kind) i,j,k,ij,klev
real(r_kind),dimension(grd%lat2*grd%lon2,npe):: sub
call mpi_alltoallv(work,grd%sendcounts_s,grd%sdispls_s,mpi_rtype,&
sub,grd%recvcounts_s,grd%rdispls_s,mpi_rtype,&
mpi_comm_world,ierror)
!$omp parallel do schedule(dynamic,1) private(k,i,j,ij,klev)
do k=1,icount
klev=ilev(k)
ij=0
select case ( chaero(k) )
case ('sulf')
do j=1,grd%lon2
do i=1,grd%lat2
ij=ij+1
ae_so4(i,j,klev)=sub(ij,k)
end do
end do
case ('bc1')
do j=1,grd%lon2
do i=1,grd%lat2
ij=ij+1
ae_bcpho(i,j,klev)=sub(ij,k)
end do
end do
case ('bc2')
do j=1,grd%lon2
do i=1,grd%lat2
ij=ij+1
ae_bcphi(i,j,klev)=sub(ij,k)
end do
end do
case ('oc1')
do j=1,grd%lon2
do i=1,grd%lat2
ij=ij+1
ae_ocpho(i,j,klev)=sub(ij,k)
end do
end do
case ('oc2')
do j=1,grd%lon2
do i=1,grd%lat2
ij=ij+1
ae_ocphi(i,j,klev)=sub(ij,k)
end do
end do
case ('dust1')
do j=1,grd%lon2
do i=1,grd%lat2
ij=ij+1
ae_d1(i,j,klev)=sub(ij,k)
end do
end do
case ('dust2')
do j=1,grd%lon2
do i=1,grd%lat2
ij=ij+1
ae_d2(i,j,klev)=sub(ij,k)
end do
end do
case ('dust3')
do j=1,grd%lon2
do i=1,grd%lat2
ij=ij+1
ae_d3(i,j,klev)=sub(ij,k)
end do
end do
case ('dust4')
do j=1,grd%lon2
do i=1,grd%lat2
ij=ij+1
ae_d4(i,j,klev)=sub(ij,k)
end do
end do
case ('dust5')
do j=1,grd%lon2
do i=1,grd%lat2
ij=ij+1
ae_d5(i,j,klev)=sub(ij,k)
end do
end do
case ('seas1')
do j=1,grd%lon2
do i=1,grd%lat2
ij=ij+1
ae_s1(i,j,klev)=sub(ij,k)
end do
end do
case ('seas2')
do j=1,grd%lon2
do i=1,grd%lat2
ij=ij+1
ae_s2(i,j,klev)=sub(ij,k)
end do
end do
case ('seas3')
do j=1,grd%lon2
do i=1,grd%lat2
ij=ij+1
ae_s3(i,j,klev)=sub(ij,k)
end do
end do
case ('seas4')
do j=1,grd%lon2
do i=1,grd%lat2
ij=ij+1
ae_s4(i,j,klev)=sub(ij,k)
end do
end do
end select
end do
!$omp end parallel do
icount=0
ilev=0
return
end subroutine aerosol_reload