module netcdfgfs_io
!$$$ module documentation block
! . . . .
! module: netcdfgfs_io
! prgmmr: Martin org: NCEP/EMC date: 2019-09-24
!
! abstract: This module contains routines which handle input/output
! operations for NCEP FV3 GFS netCDF global atmospheric and surface files.
!
! program history log:
! 2019-09-24 Martin Initial version. Based on ncepnems_io
!
! Subroutines Included:
! sub read_gfsnc - driver to read fv3gfs netcdf atmospheric and surface
! sub read_gfsnc_chem
! sub read_gfsncatm - read fv3gfs netcdf atmospheric file, scatter
! on grid to analysis subdomains
! sub read_gfsncsfc - read fv3gfs netcdf surface file, scatter on grid to
! analysis subdomains
! sub read_gfsncsfc_anl- read ncep EnKF fv3gfs netcdf surface file, scatter on grid to
! analysis subdomains
! sub write_gfsncsfc - gather/write on grid ncep surface analysis file
! sub read_gfsncnst - read ncep nst file, scatter on grid to analysis subdomains
! sub write_gfsnc_sfc_nst - gather/write on grid ncep surface & nst analysis file
! sub intrp22 - interpolate from one grid to another grid (2D)
! sub read_gfsnc_sfcnst - read sfc hist file, including sfc and nst vars, scatter on grid to analysis subdomains
!
! Variable Definitions:
!
! language: f90
! machine:
!
! NOTE: This module adds capability to read netCDF FV3 first guess files
! and to write netCDF FV3 analysis files using the ncio interface
! Using this is controled by a namelist argument "use_gfs_ncio"
!
!
!$$$ end documentation block
use constants, only: zero,one,fv,r60,r3600
implicit none
private
public read_gfsnc
public read_gfsnc_chem
public read_gfsncatm
public read_gfsncsfc
public read_gfsncsfc_anl
public write_gfsncsfc
public read_gfsncnst
public write_gfsnc_sfc_nst
public intrp22
public tran_gfsncsfc
public write_gfsncatm
interface read_gfsnc
module procedure read_
end interface
interface read_gfsnc_chem
module procedure read_chem_
end interface
interface read_gfsncatm
module procedure read_atm_
end interface
interface read_gfsncsfc
module procedure read_gfsncsfc_
end interface
interface read_gfsncsfc_anl
module procedure read_gfsncsfc_anl_
end interface
interface read_gfsncnst
module procedure read_gfsncnst_
end interface
interface write_gfsncsfc
module procedure write_sfc_
end interface
interface write_gfsnc_sfc_nst
module procedure write_sfc_nst_
end interface
interface write_gfsncatm
module procedure write_atm_
end interface
character(len=*),parameter::myname='netcdfgfs_io'
contains
subroutine read_
!$$$ subprogram documentation block
! . . .
! subprogram: read_gfsnc
!
! prgrmmr: Martin
!
! abstract:
!
! program history log:
! 2019-09-24 Martin - create routine based on read_nems
!
! input argument list:
!
! output argument list:
!
! attributes:
! language: f90
!
!$$$ end documentation block
use kinds, only: i_kind,r_kind
use constants, only: qcmin
use gridmod, only: sp_a,grd_a,lat2,lon2,nsig
use guess_grids, only: ifilesig,nfldsig,ntguessig
use gsi_metguess_mod, only: gsi_metguess_bundle
use gsi_bundlemod, only: gsi_bundlegetpointer
use gsi_bundlemod, only: gsi_bundlecreate
use gsi_bundlemod, only: gsi_grid
use gsi_bundlemod, only: gsi_gridcreate
use gsi_bundlemod, only: gsi_bundle
use gsi_bundlemod, only: gsi_bundledestroy
use general_sub2grid_mod, only: sub2grid_info,general_sub2grid_create_info,general_sub2grid_destroy_info
use mpimod, only: npe,mype
use cloud_efr_mod, only: cloud_calc_gfs,set_cloud_lower_bound
use gridmod, only: fv3_full_hydro
implicit none
character(len=*),parameter::myname_=myname//'*read_'
character(24) filename
integer(i_kind):: it, istatus, inner_vars, num_fields
integer(i_kind):: i,j,k
real(r_kind),pointer,dimension(:,: ):: ges_ps_it =>NULL()
real(r_kind),pointer,dimension(:,: ):: ges_z_it =>NULL()
real(r_kind),pointer,dimension(:,:,:):: ges_u_it =>NULL()
real(r_kind),pointer,dimension(:,:,:):: ges_v_it =>NULL()
real(r_kind),pointer,dimension(:,:,:):: ges_div_it =>NULL()
real(r_kind),pointer,dimension(:,:,:):: ges_vor_it =>NULL()
real(r_kind),pointer,dimension(:,:,:):: ges_tv_it =>NULL()
real(r_kind),pointer,dimension(:,:,:):: ges_q_it =>NULL()
real(r_kind),pointer,dimension(:,:,:):: ges_oz_it =>NULL()
real(r_kind),pointer,dimension(:,:,:):: ges_cwmr_it=>NULL()
real(r_kind),pointer,dimension(:,:,:):: ges_ql_it => NULL()
real(r_kind),pointer,dimension(:,:,:):: ges_qi_it => NULL()
real(r_kind),pointer,dimension(:,:,:):: ges_qr_it => NULL()
real(r_kind),pointer,dimension(:,:,:):: ges_qs_it => NULL()
real(r_kind),pointer,dimension(:,:,:):: ges_qg_it => NULL()
real(r_kind),pointer,dimension(:,:,:):: ges_cf_it => NULL()
type(sub2grid_info) :: grd_t
logical regional
logical:: l_cld_derived,zflag,inithead
type(gsi_bundle) :: atm_bundle
type(gsi_grid) :: atm_grid
integer(i_kind),parameter :: n2d=2
! integer(i_kind),parameter :: n3d=8
integer(i_kind),parameter :: n3d=14
character(len=4), parameter :: vars2d(n2d) = (/ 'z ', 'ps ' /)
! character(len=4), parameter :: vars3d(n3d) = (/ 'u ', 'v ', &
! 'vor ', 'div ', &
! 'tv ', 'q ', &
! 'cw ', 'oz ' /)
character(len=4), parameter :: vars3d(n3d) = (/ 'u ', 'v ', &
'vor ', 'div ', &
'tv ', 'q ', &
'cw ', 'oz ', &
'ql ', 'qi ', &
'qr ', 'qs ', &
'qg ', 'cf ' /)
real(r_kind),pointer,dimension(:,:):: ptr2d =>NULL()
real(r_kind),pointer,dimension(:,:,:):: ptr3d =>NULL()
regional=.false.
inner_vars=1
num_fields=min(14*grd_a%nsig+2,npe)
! Create temporary communication information fore read routines
call general_sub2grid_create_info(grd_t,inner_vars,grd_a%nlat,grd_a%nlon, &
grd_a%nsig,num_fields,regional)
! Allocate bundle used for reading members
call gsi_gridcreate(atm_grid,lat2,lon2,nsig)
call gsi_bundlecreate(atm_bundle,atm_grid,'aux-atm-read',istatus,names2d=vars2d,names3d=vars3d)
if(istatus/=0) then
write(6,*) myname_,': trouble creating atm_bundle'
call stop2(999)
endif
do it=1,nfldsig
write(filename,'(''sigf'',i2.2)') ifilesig(it)
! Read background fields into bundle
if (fv3_full_hydro) then
if (mype==0) write(*,*) 'calling general_read_gfsatm_allhydro_nc', it
call general_read_gfsatm_allhydro_nc(grd_t,sp_a,filename,.true.,.true.,.true.,&
atm_bundle,istatus) ! this loads cloud and precip
if (mype==0) write(*,*) 'done with general_read_gfsatm_allhydro_nc', it
else
if (mype==0) write(*,*) 'calling general_read_gfsatm_nc'
call general_read_gfsatm_nc(grd_t,sp_a,filename,.true.,.true.,.true.,&
atm_bundle,.true.,istatus)
if (mype==0) write(*,*) 'done with general_read_gfsatm_nc'
end if
inithead=.false.
zflag=.false.
! Set values to actual MetGuess fields
call set_guess_
if (it==ntguessig) then
if (mype==0) write(6,*)'Print guess field ... after set_guess'
call prt_guess('guess')
endif
if (fv3_full_hydro) then
do k=1, nsig
do j=1, lon2
do i=1, lat2
! set lower bound to hydrometeors
if (associated(ges_ql_it)) ges_ql_it(i,j,k) = max(qcmin,ges_ql_it(i,j,k))
if (associated(ges_qi_it)) ges_qi_it(i,j,k) = max(qcmin,ges_qi_it(i,j,k))
if (associated(ges_qr_it)) ges_qr_it(i,j,k) = max(qcmin,ges_qr_it(i,j,k))
if (associated(ges_qs_it)) ges_qs_it(i,j,k) = max(qcmin,ges_qs_it(i,j,k))
if (associated(ges_qg_it)) ges_qg_it(i,j,k) = max(qcmin,ges_qg_it(i,j,k))
if (associated(ges_cf_it)) ges_cf_it(i,j,k) = min(max(zero,ges_cf_it(i,j,k)),one)
enddo
enddo
enddo
else
l_cld_derived = associated(ges_cwmr_it).and.&
associated(ges_q_it) .and.&
associated(ges_ql_it) .and.&
associated(ges_qi_it) .and.&
associated(ges_tv_it)
! call set_cloud_lower_bound(ges_cwmr_it)
if (mype==0) write(6,*)'READ_GFS_NETCDF: l_cld_derived = ', l_cld_derived
if (l_cld_derived) then
call cloud_calc_gfs(ges_ql_it,ges_qi_it,ges_cwmr_it,ges_q_it,ges_tv_it,.true.)
end if
end if
if (it==ntguessig) then
if (mype==0) write(6,*)'Print guess field ... after set lower bound for hydrometers'
call prt_guess('guess')
endif
end do
call general_sub2grid_destroy_info(grd_t)
call gsi_bundledestroy(atm_bundle,istatus)
contains
subroutine set_guess_
call gsi_bundlegetpointer (atm_bundle,'ps',ptr2d,istatus)
if (istatus==0) then
call gsi_bundlegetpointer (gsi_metguess_bundle(it),'ps',ges_ps_it ,istatus)
if(istatus==0) ges_ps_it = ptr2d
endif
call gsi_bundlegetpointer (atm_bundle,'z',ptr2d,istatus)
if (istatus==0) then
call gsi_bundlegetpointer (gsi_metguess_bundle(it),'z' ,ges_z_it ,istatus)
if(istatus==0) ges_z_it = ptr2d
endif
call gsi_bundlegetpointer (atm_bundle,'u',ptr3d,istatus)
if (istatus==0) then
call gsi_bundlegetpointer (gsi_metguess_bundle(it),'u' ,ges_u_it ,istatus)
if(istatus==0) ges_u_it = ptr3d
endif
call gsi_bundlegetpointer (atm_bundle,'v',ptr3d,istatus)
if (istatus==0) then
call gsi_bundlegetpointer (gsi_metguess_bundle(it),'v' ,ges_v_it ,istatus)
if(istatus==0) ges_v_it = ptr3d
endif
call gsi_bundlegetpointer (atm_bundle,'vor',ptr3d,istatus)
if (istatus==0) then
call gsi_bundlegetpointer (gsi_metguess_bundle(it),'vor',ges_vor_it,istatus)
if(istatus==0) ges_vor_it = ptr3d
endif
call gsi_bundlegetpointer (atm_bundle,'div',ptr3d,istatus)
if (istatus==0) then
call gsi_bundlegetpointer (gsi_metguess_bundle(it),'div',ges_div_it,istatus)
if(istatus==0) ges_div_it = ptr3d
endif
call gsi_bundlegetpointer (atm_bundle,'tv',ptr3d,istatus)
if (istatus==0) then
call gsi_bundlegetpointer (gsi_metguess_bundle(it),'tv',ges_tv_it ,istatus)
if(istatus==0) ges_tv_it = ptr3d
endif
call gsi_bundlegetpointer (atm_bundle,'q',ptr3d,istatus)
if (istatus==0) then
call gsi_bundlegetpointer (gsi_metguess_bundle(it),'q' ,ges_q_it ,istatus)
if(istatus==0) ges_q_it = ptr3d
endif
call gsi_bundlegetpointer (atm_bundle,'oz',ptr3d,istatus)
if (istatus==0) then
call gsi_bundlegetpointer (gsi_metguess_bundle(it),'oz',ges_oz_it ,istatus)
if(istatus==0) ges_oz_it = ptr3d
endif
call gsi_bundlegetpointer (atm_bundle,'cw',ptr3d,istatus)
if (istatus==0) then
call gsi_bundlegetpointer (gsi_metguess_bundle(it),'cw',ges_cwmr_it,istatus)
if(istatus==0) ges_cwmr_it = ptr3d
endif
call gsi_bundlegetpointer (atm_bundle,'ql',ptr3d,istatus)
if (istatus==0) then
call gsi_bundlegetpointer (gsi_metguess_bundle(it),'ql',ges_ql_it,istatus)
if(istatus==0) ges_ql_it = ptr3d
endif
call gsi_bundlegetpointer (atm_bundle,'qi',ptr3d,istatus)
if (istatus==0) then
call gsi_bundlegetpointer (gsi_metguess_bundle(it),'qi',ges_qi_it,istatus)
if(istatus==0) ges_qi_it = ptr3d
endif
call gsi_bundlegetpointer (atm_bundle,'qr',ptr3d,istatus)
if (istatus==0) then
call gsi_bundlegetpointer (gsi_metguess_bundle(it),'qr',ges_qr_it,istatus)
if(istatus==0) ges_qr_it = ptr3d
endif
call gsi_bundlegetpointer (atm_bundle,'qs',ptr3d,istatus)
if (istatus==0) then
call gsi_bundlegetpointer (gsi_metguess_bundle(it),'qs',ges_qs_it,istatus)
if(istatus==0) ges_qs_it = ptr3d
endif
call gsi_bundlegetpointer (atm_bundle,'qg',ptr3d,istatus)
if (istatus==0) then
call gsi_bundlegetpointer (gsi_metguess_bundle(it),'qg',ges_qg_it,istatus)
if(istatus==0) ges_qg_it = ptr3d
endif
call gsi_bundlegetpointer (atm_bundle,'cf',ptr3d,istatus)
if (istatus==0) then
call gsi_bundlegetpointer (gsi_metguess_bundle(it),'cf',ges_cf_it,istatus)
if(istatus==0) ges_cf_it = ptr3d
endif
end subroutine set_guess_
end subroutine read_
subroutine read_chem_ ( iyear, month,idd )
!$$$ subprogram documentation block
! . . .
! subprogram: read_gfsnc_chem
!
! prgrmmr: martin
!
! abstract: fills chemguess_bundle with GFS chemistry.
!
! remarks:
! 1. Right now, only CO2 is done and even this is treated
! as constant througout the assimialation window.
! 2. iyear and month could come from obsmod, but logically
! this program should never depend on obsmod
!
!
! program history log:
! 2019-09-24 martin - initial code, based on read_nems_chem
!
! input argument list:
!
! output argument list:
!
! attributes:
! language: f90
!
!$$$ end documentation block
use kinds, only: i_kind, r_kind
use mpimod, only: mype
use gridmod, only: lat2,lon2,nsig,nlat,rlats,istart
use ncepgfs_ghg, only: read_gfsco2
use guess_grids, only: nfldsig
use gsi_bundlemod, only: gsi_bundlegetpointer
use gsi_chemguess_mod, only: gsi_chemguess_bundle
use gsi_chemguess_mod, only: gsi_chemguess_get
implicit none
! Declared argument list
integer(i_kind), intent(in):: iyear
integer(i_kind), intent(in):: month
integer(i_kind), intent(in):: idd
! Declare local variables
integer(i_kind) :: igfsco2, i, j, n, iret
real(r_kind),dimension(lat2):: xlats
real(r_kind),pointer,dimension(:,:,:)::p_co2=>NULL()
real(r_kind),pointer,dimension(:,:,:)::ptr3d=>NULL()
if(.not.associated(gsi_chemguess_bundle)) return
call gsi_bundlegetpointer(gsi_chemguess_bundle(1),'co2',p_co2,iret)
if(iret /= 0) return
! Get subdomain latitude array
j = mype + 1
do i = 1, lat2
n = min(max(1, istart(j)+i-2), nlat)
xlats(i) = rlats(n)
enddo
! Read in CO2
call gsi_chemguess_get ( 'i4crtm::co2', igfsco2, iret )
call read_gfsco2 ( iyear,month,idd,igfsco2,xlats,&
lat2,lon2,nsig,mype, p_co2 )
! Approximation: setting all times co2 values equal to the daily co2 values
do n = 2, nfldsig
call gsi_bundlegetpointer(gsi_chemguess_bundle(n),'co2',ptr3d,iret)
ptr3d = p_co2
enddo
end subroutine read_chem_
subroutine read_atm_ (grd,filename,sp_a,uvflag,vordivflag,zflag, &
g_z,g_ps,g_vor,g_div,g_u,g_v,&
g_tv,g_q,g_cwmr,g_oz)
!$$$ subprogram documentation block
! . . . .
! subprogram: read_gfsncatm read GFS netCDF atm and send to all mpi tasks
! prgmmr: Martin org: NCEP/EMC date: 2019-09-23
!
! abstract: read ncep netcdf/gfs atmospheric guess field and
! scatter to subdomains
!
! program history log:
! 2019-09-23 Martin Initial version. Based on sub read_nemsatm
!
! 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)
! uvflag - logical to use u,v (.true.) or st,vp (.false.) perturbations
! vordivflag - logical to determine if routine should output vorticity and
! divergence
! zflag - logical to determine if surface height field should be output
!
! output argument list:
! g_* - guess fields
!
! attributes:
! language: f90
!
!$$$
use kinds, only: r_kind,i_kind, r_single
use gridmod, only: ntracer,ncloud,reload,itotsub
use general_commvars_mod, only: fill_ns,filluv_ns,fill2_ns,filluv2_ns,ltosj_s,ltosi_s
use general_specmod, only: spec_vars
use general_sub2grid_mod, only: sub2grid_info
use mpimod, only: npe,mpi_comm_world,ierror,mpi_rtype,mype
use module_ncio, only: Dataset, Variable, Dimension, open_dataset,&
quantize_data,close_dataset, get_dim, read_vardata, get_idate_from_time_units
use egrid2agrid_mod,only: g_egrid2agrid,g_create_egrid2agrid,egrid2agrid_parm,destroy_egrid2agrid
use constants, only: two,pi,half,deg2rad
implicit none
! Declare local parameters
real(r_kind),parameter:: r0_001 = 0.001_r_kind
! Declare passed variables
type(sub2grid_info) ,intent(in ) :: grd
character(len=24) ,intent(in ) :: filename
logical ,intent(in ) :: uvflag,vordivflag,zflag
real(r_kind),dimension(grd%lat2,grd%lon2) ,intent( out) :: g_z,g_ps
real(r_kind),dimension(grd%lat2,grd%lon2,grd%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
! Declare local variables
character(len=120) :: my_name = 'READ_GFSNCATM'
integer(i_kind),dimension(6):: idate
integer(i_kind),dimension(4):: odate
integer(i_kind) :: nlatm2,nflds
integer(i_kind) :: k,icount,icount_prev,mm1,i,j,kk,kr
integer(i_kind) :: mype_hs, mype_ps,nord_int
integer(i_kind) :: latb, lonb, levs
real(r_kind),allocatable,dimension(:,:) :: grid, grid_v, &
grid_vor, grid_div, grid_b, grid_b2
real(r_kind),allocatable,dimension(:,:,:) :: grid_c, grid2, grid_c2
real(r_kind),allocatable,dimension(:) :: work, work_vor, work_div, &
work_v
real(r_kind),allocatable,dimension(:,:) :: sub, sub_vor, sub_div, &
sub_v
real(r_kind),dimension(sp_a%nc):: spec_vor,spec_div
real(r_kind),allocatable,dimension(:) :: rlats,rlons,clons,slons
real(r_kind),allocatable,dimension(:) :: rlats_tmp,rlons_tmp
real(r_single),allocatable, dimension(:,:) :: rwork2d, rwork2d1
real(r_single),allocatable, dimension(:,:,:) :: rwork3d, rwork3d1
real(r_single),allocatable, dimension(:) :: fhour
logical diff_res,eqspace
logical,dimension(1) :: vector
type(egrid2agrid_parm) :: p_high
type(Dataset) :: atmges
type(Dimension) :: ncdim
!******************************************************************************
! Initialize variables used below
mm1=mype+1
mype_hs=min(1,npe-1)
mype_ps=0
nlatm2=grd%nlat-2
nflds=5*grd%nsig+1
if(zflag) nflds=nflds+1
if(vordivflag .or. .not. uvflag)nflds=nflds+2*grd%nsig
! nflds=npe
nflds=grd%nsig
levs=grd%nsig
allocate( work(grd%itotsub),work_v(grd%itotsub) )
work=zero
work_v=zero
allocate( sub(grd%lat2*grd%lon2,max(grd%nsig,npe)),sub_v(grd%lat2*grd%lon2,max(grd%nsig,npe)) )
allocate( sub_div(grd%lat2*grd%lon2,max(grd%nsig,npe)),sub_vor(grd%lat2*grd%lon2,max(grd%nsig,npe)) )
if(mype < nflds)then
! open the netCDF file
atmges = open_dataset(filename)
! get dimension sizes
ncdim = get_dim(atmges, 'grid_xt'); lonb = ncdim%len
ncdim = get_dim(atmges, 'grid_yt'); latb = ncdim%len
ncdim = get_dim(atmges, 'pfull'); levs = ncdim%len
! get time information
idate = get_idate_from_time_units(atmges)
odate(1) = idate(4) !hour
odate(2) = idate(2) !month
odate(3) = idate(3) !day
odate(4) = idate(1) !year
call read_vardata(atmges, 'time', fhour) ! might need to change this to attribute later
! depends on model changes from Jeff Whitaker
!
! 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
! get lat/lon coordinates
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))
allocate( grid_b2(lonb,latb),grid_c2(latb+2,lonb,1))
end if
allocate( rlats(latb+2),rlons(lonb),clons(lonb),slons(lonb))
call read_vardata(atmges, 'grid_xt', rlons_tmp)
call read_vardata(atmges, 'grid_yt', rlats_tmp)
do j=1,latb
rlats(j+1)=deg2rad*rlats_tmp(j)
end do
do j=1,lonb
rlons(j)=deg2rad*rlons_tmp(j)
end do
deallocate(rlats_tmp,rlons_tmp)
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)
deallocate(rlats,rlons)
end if
!
! Load values into rows for south and north pole before scattering
!
! Terrain: scatter to all mpi tasks
!
if(zflag)then
if (mype==mype_hs) then
call read_vardata(atmges, 'hgtsfc', rwork2d)
if(diff_res)then
grid_b=rwork2d
vector(1)=.false.
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,itotsub
i=ltosi_s(kk)
j=ltosj_s(kk)
work(kk)=grid2(i,j,1)
end do
else
grid=rwork2d
call fill_ns(grid,work)
end if
endif
call mpi_scatterv(work,grd%ijn_s,grd%displs_s,mpi_rtype,&
g_z,grd%ijn_s(mm1),mpi_rtype,mype_hs,mpi_comm_world,ierror)
end if
! Surface pressure: same procedure as terrain, but handled by task mype_ps
!
if (mype==mype_ps) then
call read_vardata(atmges, 'pressfc', rwork2d)
rwork2d = r0_001*rwork2d
if(diff_res)then
vector(1)=.false.
grid_b=rwork2d
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,itotsub
i=ltosi_s(kk)
j=ltosj_s(kk)
work(kk)=grid2(i,j,1)
end do
else
grid=rwork2d
call fill_ns(grid,work)
endif
endif
call mpi_scatterv(work,grd%ijn_s,grd%displs_s,mpi_rtype,&
g_ps,grd%ijn_s(mm1),mpi_rtype,mype_ps,mpi_comm_world,ierror)
! Divergence and voriticity. Compute u and v from div and vor
sub_vor=zero
sub_div=zero
sub =zero
sub_v =zero
icount =0
icount_prev=1
allocate( work_vor(grd%itotsub),work_div(grd%itotsub) )
call read_vardata(atmges, 'ugrd', rwork3d)
call read_vardata(atmges, 'vgrd', rwork3d1)
! TODO CRM, would above be faster if only done on one PE and then distributed?
do k=1,levs
kr = levs+1-k ! netcdf is top to bottom need to flip
icount=icount+1
if (mype==mod(icount-1,npe)) then
! Convert grid u,v to div and vor
if(diff_res)then
grid_b=rwork3d(:,:,kr)
grid_b2=rwork3d1(:,:,kr)
vector(1)=.true.
call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons)
call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector)
do kk=1,itotsub
i=ltosi_s(kk)
j=ltosj_s(kk)
work(kk)=grid2(i,j,1)
end do
do j=1,grd%nlon
do i=2,grd%nlat-1
grid(j,grd%nlat-i)=grid2(i,j,1)
end do
end do
call g_egrid2agrid(p_high,grid_c2,grid2,1,1,vector)
do kk=1,itotsub
i=ltosi_s(kk)
j=ltosj_s(kk)
work_v(kk)=grid2(i,j,1)
end do
do j=1,grd%nlon
do i=2,grd%nlat-1
grid_v(j,grd%nlat-i)=grid2(i,j,1)
end do
end do
else
grid=rwork3d(:,:,kr)
grid_v=rwork3d1(:,:,kr)
call filluv_ns(grid,grid_v,work,work_v)
end if
if(vordivflag .or. .not. uvflag)then
allocate( grid_vor(grd%nlon,nlatm2), grid_div(grd%nlon,nlatm2) )
call general_sptez_v(sp_a,spec_div,spec_vor,grid,grid_v,-1)
call general_sptez_s_b(sp_a,sp_a,spec_div,grid_div,1)
call general_sptez_s_b(sp_a,sp_a,spec_vor,grid_vor,1)
! Load values into rows for south and north pole
call fill_ns(grid_div,work_div)
call fill_ns(grid_vor,work_vor)
deallocate(grid_vor,grid_div)
end if
endif
! Scatter to sub
if (mod(icount,npe)==0 .or. icount==levs) then
if(vordivflag .or. .not. uvflag)then
call mpi_alltoallv(work_vor,grd%ijn_s,grd%displs_s,mpi_rtype,&
sub_vor(1,icount_prev),grd%irc_s,grd%ird_s,mpi_rtype,&
mpi_comm_world,ierror)
call mpi_alltoallv(work_div,grd%ijn_s,grd%displs_s,mpi_rtype,&
sub_div(1,icount_prev),grd%irc_s,grd%ird_s,mpi_rtype,&
mpi_comm_world,ierror)
end if
if(uvflag)then
call mpi_alltoallv(work,grd%ijn_s,grd%displs_s,mpi_rtype,&
sub(1,icount_prev),grd%irc_s,grd%ird_s,mpi_rtype,&
mpi_comm_world,ierror)
call mpi_alltoallv(work_v,grd%ijn_s,grd%displs_s,mpi_rtype,&
sub_v(1,icount_prev),grd%irc_s,grd%ird_s,mpi_rtype,&
mpi_comm_world,ierror)
end if
icount_prev=icount+1
endif
end do
deallocate(work_vor,work_div)
! 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,g_u)
call reload(sub_v,g_v)
deallocate(sub_vor,sub_div)
! Thermodynamic variable and Specific humidity: communicate to all tasks
!
sub=zero
icount=0
icount_prev=1
call read_vardata(atmges, 'spfh', rwork3d)
call read_vardata(atmges, 'tmp', rwork3d1)
allocate(rwork2d1(lonb,latb))
do k=1,levs
kr = levs+1-k ! netcdf is top to bottom need to flip
icount=icount+1
if (mype==mod(icount-1,npe)) then
if(diff_res)then
grid_b=rwork3d(:,:,kr)
vector(1)=.false.
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,itotsub
i=ltosi_s(kk)
j=ltosj_s(kk)
work(kk)=grid2(i,j,1)
end do
else
grid=rwork3d(:,:,kr)
call fill_ns(grid,work)
end if
rwork2d1 = rwork3d1(:,:,kr)*(one+fv*rwork3d(:,:,kr))
if(diff_res)then
grid_b=rwork2d1
vector(1)=.false.
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,itotsub
i=ltosi_s(kk)
j=ltosj_s(kk)
work_v(kk)=grid2(i,j,1)
end do
else
grid_v=rwork2d1
call fill_ns(grid_v,work_v)
end if
endif
if (mod(icount,npe)==0 .or. icount==levs) then
call mpi_alltoallv(work_v,grd%ijn_s,grd%displs_s,mpi_rtype,&
sub_v(1,icount_prev),grd%irc_s,grd%ird_s,mpi_rtype,&
mpi_comm_world,ierror)
call mpi_alltoallv(work,grd%ijn_s,grd%displs_s,mpi_rtype,&
sub(1,icount_prev),grd%irc_s,grd%ird_s,mpi_rtype,&
mpi_comm_world,ierror)
icount_prev=icount+1
endif
end do
deallocate(rwork2d1)
call reload(sub_v,g_tv)
call reload(sub,g_q)
deallocate(sub_v,work_v)
sub=zero
icount=0
icount_prev=1
call read_vardata(atmges, 'o3mr', rwork3d) ! need k thrown in here somewhere
do k=1,levs
kr = levs+1-k ! netcdf is top to bottom need to flip
icount=icount+1
if (mype==mod(icount-1,npe)) then
if(diff_res)then
grid_b=rwork3d(:,:,kr)
vector(1)=.false.
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,itotsub
i=ltosi_s(kk)
j=ltosj_s(kk)
work(kk)=grid2(i,j,1)
end do
else
grid=rwork3d(:,:,kr)
call fill_ns(grid,work)
end if
endif
if (mod(icount,npe)==0 .or. icount==levs) then
call mpi_alltoallv(work,grd%ijn_s,grd%displs_s,mpi_rtype,&
sub(1,icount_prev),grd%irc_s,grd%ird_s,mpi_rtype,&
mpi_comm_world,ierror)
icount_prev=icount+1
endif
end do
call reload(sub,g_oz)
! Cloud condensate mixing ratio.
if (ntracer>2 .or. ncloud>=1) then
sub=zero
icount=0
icount_prev=1
call read_vardata(atmges, 'clwmr', rwork3d)
call read_vardata(atmges, 'icmr', rwork3d1)
do k=1,levs
kr = levs+1-k ! netcdf is top to bottom need to flip
icount=icount+1
if (mype==mod(icount-1,npe)) then
rwork2d = rwork3d(:,:,kr) + rwork3d1(:,:,kr)
if(diff_res)then
grid_b=rwork2d
vector(1)=.false.
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,itotsub
i=ltosi_s(kk)
j=ltosj_s(kk)
work(kk)=grid2(i,j,1)
end do
else
grid=rwork2d
call fill_ns(grid,work)
end if
endif
if (mod(icount,npe)==0 .or. icount==levs) then
call mpi_alltoallv(work,grd%ijn_s,grd%displs_s,mpi_rtype,&
sub(1,icount_prev),grd%irc_s,grd%ird_s,mpi_rtype,&
mpi_comm_world,ierror)
icount_prev=icount+1
endif
end do
call reload(sub,g_cwmr)
else
g_cwmr = zero
endif
if(mype < nflds)then
if(diff_res) deallocate(grid_b,grid_b2,grid_c,grid_c2,grid2)
call destroy_egrid2agrid(p_high)
deallocate(clons,slons)
deallocate(grid,grid_v)
call close_dataset(atmges)
end if
deallocate(work,sub)
! Print date/time stamp
if ( mype == 0 ) write(6, &
'(a,'': ges read/scatter,lonb,latb,levs= '',3i6,'',hour= '',f4.1,'',idate= '',4i5)') &
trim(my_name),lonb,latb,levs,fhour,odate
end subroutine read_atm_
subroutine read_sfc_(sfct,soil_moi,sno,soil_temp,veg_frac,fact10,sfc_rough, &
veg_type,soil_type,terrain,isli,use_sfc_any, &
tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d)
!$$$ subprogram documentation block
! . . . .
! subprogram: read_sfc_ read netCDF sfc hist file
! prgmmr: Martin org: NCEP/EMC date: 2019-09-23
!
! abstract: read nems sfc & nst combined file
!
! program history log:
! 2019-09-23 Martin Initial version.
!
! input argument list:
! use_sfc_any - true if any processor uses extra surface fields
!
! output argument list:
! sfct - surface temperature (skin temp)
! soil_moi - soil moisture of first layer
! sno - snow depth
! soil_temp - soil temperature of first layer
! veg_frac - vegetation fraction
! fact10 - 10 meter wind factor
! sfc_rough - surface roughness
! veg_type - vegetation type
! soil_type - soil type
! terrain - terrain height
! isli - sea/land/ice mask
! tref - optional, oceanic foundation temperature
! dt_cool - optional, sub-layer cooling amount at sub-skin layer
! z_c - optional, depth of sub-layer cooling layer
! dt_warm - optional, diurnal warming amount at sea surface
! z_w - optional, depth of diurnal warming layer
! c_0 - optional, coefficient to calculate d(Tz)/d(tr) in dimensionless
! c_d - optional, coefficient to calculate d(Tz)/d(tr) in m^-1
! w_0 - optional, coefficient to calculate d(Tz)/d(tr) in dimensionless
! w_d - optional, coefficient to calculate d(Tz)/d(tr) in m^-1
!
! attributes:
! language: f90
!
!$$$
use mpimod, only: mype
use kinds, only: r_kind,i_kind,r_single
use gridmod, only: nlat_sfc,nlon_sfc
use guess_grids, only: nfldsfc,ifilesfc
use constants, only: zero,two
use module_ncio, only: Dataset, Variable, Dimension, open_dataset,&
close_dataset, get_dim, read_vardata, get_idate_from_time_units
implicit none
! Declare passed variables
logical, intent(in) :: use_sfc_any
real(r_single), dimension(nlat_sfc,nlon_sfc,nfldsfc), intent(out) :: sfct,soil_moi,sno,soil_temp,veg_frac,fact10,sfc_rough
real(r_single), dimension(nlat_sfc,nlon_sfc), intent(out) :: veg_type,soil_type,terrain
integer(i_kind), dimension(nlat_sfc,nlon_sfc), intent(out) :: isli
real(r_single), optional, dimension(nlat_sfc,nlon_sfc,nfldsfc), intent(out) :: tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d
! Declare local parameters
integer(i_kind), parameter :: nsfc_all=11
integer(i_kind),dimension(6):: idate
integer(i_kind),dimension(4):: odate
! Declare local variables
real(r_single), dimension(nlat_sfc,nlon_sfc,nfldsfc) :: xt
character(len=24) :: filename
character(len=120) :: my_name = 'READ_SFCNST'
integer(i_kind) :: i,j,it,n,nsfc,iret
integer(i_kind) :: lonb, latb
real(r_single),allocatable, dimension(:) :: fhour
real(r_single), allocatable, dimension(:,:) :: work,outtmp
type(Dataset) :: sfcges
type(Dimension) :: ncdim
!-----------------------------------------------------------------------------
do it = 1, nfldsfc
! read a surface history file on the task
write(filename,200)ifilesfc(it)
200 format('sfcf',i2.2)
! open the netCDF file
sfcges = open_dataset(filename,errcode=iret)
if (iret/=0) then
write(6,*) trim(my_name),': ***ERROR*** ',trim(filename),' NOT AVAILABLE: PROGRAM STOPS, later'
endif
! get dimension sizes
ncdim = get_dim(sfcges, 'grid_xt'); lonb = ncdim%len
ncdim = get_dim(sfcges, 'grid_yt'); latb = ncdim%len
! get time information
idate = get_idate_from_time_units(sfcges)
odate(1) = idate(4) !hour
odate(2) = idate(2) !month
odate(3) = idate(3) !day
odate(4) = idate(1) !year
call read_vardata(sfcges, 'time', fhour) ! might need to change this to attribute later
! depends on model changes from Jeff Whitaker
if ( (latb /= nlat_sfc-2) .or. (lonb /= nlon_sfc) ) then
if ( mype == 0 ) write(6, &
'(a,'': inconsistent spatial dimension '',''nlon,nlatm2 = '',2(i4,tr1),''-vs- sfc file lonb,latb = '',i4)') &
trim(my_name),nlon_sfc,nlat_sfc-2,lonb,latb
call stop2(102)
endif
!
! Read the surface records (lonb, latb) and convert to GSI array pattern (nlat_sfc,nlon_sfc)
! Follow the read order sfcio in ncepgfs_io
!
allocate(work(lonb,latb))
work = zero
if(it == 1)then
nsfc=nsfc_all
else
nsfc=nsfc_all-4
end if
do n = 1, nsfc
if (n == 1) then ! skin temperature
! Tsea
call read_vardata(sfcges, 'tmpsfc', work)
call tran_gfsncsfc(work,sfct(1,1,it),lonb,latb)
elseif(n == 2 .and. use_sfc_any) then ! soil moisture
! smc/soilw
call read_vardata(sfcges, 'soilw1', work)
call tran_gfsncsfc(work,soil_moi(1,1,it),lonb,latb)
elseif(n == 3) then ! snow depth
call read_vardata(sfcges, 'weasd', work)
call tran_gfsncsfc(work,sno(1,1,it),lonb,latb)
elseif(n == 4 .and. use_sfc_any) then ! soil temperature
! stc/tmp
call read_vardata(sfcges, 'soilt1', work)
call tran_gfsncsfc(work,soil_temp(1,1,it),lonb,latb)
elseif(n == 5 .and. use_sfc_any) then ! vegetation cover
! vfrac
call read_vardata(sfcges, 'veg', work)
call tran_gfsncsfc(work,veg_frac(1,1,it),lonb,latb)
elseif(n == 6) then ! 10m wind factor
! f10m
call read_vardata(sfcges, 'f10m', work)
call tran_gfsncsfc(work,fact10(1,1,it),lonb,latb)
elseif(n == 7) then ! suface roughness
! zorl
call read_vardata(sfcges, 'sfcr', work)
call tran_gfsncsfc(work,sfc_rough(1,1,it),lonb,latb)
elseif(n == 8 .and. use_sfc_any) then ! vegetation type
! vtype
call read_vardata(sfcges, 'vtype', work)
call tran_gfsncsfc(work,veg_type,lonb,latb)
elseif(n == 9 .and. use_sfc_any) then ! soil type
! stype
call read_vardata(sfcges, 'sotyp', work)
call tran_gfsncsfc(work,soil_type,lonb,latb)
elseif(n == 10) then ! terrain
! orog
call read_vardata(sfcges, 'orog', work)
call tran_gfsncsfc(work,terrain,lonb,latb)
elseif(n == 11) then ! sea/land/ice flag
! slmsk
call read_vardata(sfcges, 'land', work)
allocate(outtmp(latb+2,lonb))
call tran_gfsncsfc(work,outtmp,lonb,latb)
do j=1,lonb
do i=1,latb+2
isli(i,j) = nint(outtmp(i,j))
end do
end do
deallocate(outtmp)
endif
! End of loop over data records
enddo
if( present(tref) ) then
if ( mype == 0 ) write(6,*) ' read 9 optional NSST variables '
call read_vardata(sfcges, 'tref', work)
call tran_gfsncsfc(work,tref(1,1,it),lonb,latb)
call read_vardata(sfcges, 'dtcool', work)
call tran_gfsncsfc(work,dt_cool(1,1,it),lonb,latb)
call read_vardata(sfcges, 'zc', work)
call tran_gfsncsfc(work,z_c(1,1,it),lonb,latb)
call read_vardata(sfcges, 'xt', work)
call tran_gfsncsfc(work,xt(1,1,it),lonb,latb)
call read_vardata(sfcges, 'xz', work)
call tran_gfsncsfc(work,z_w(1,1,it),lonb,latb)
call read_vardata(sfcges, 'c0', work)
call tran_gfsncsfc(work,c_0(1,1,it),lonb,latb)
call read_vardata(sfcges, 'cd', work)
call tran_gfsncsfc(work,c_d(1,1,it),lonb,latb)
call read_vardata(sfcges, 'w0', work)
call tran_gfsncsfc(work,w_0(1,1,it),lonb,latb)
call read_vardata(sfcges, 'wd', work)
call tran_gfsncsfc(work,w_d(1,1,it),lonb,latb)
!
! Get diurnal warming amout at z=0
!
do j = 1,nlon_sfc
do i = 1,nlat_sfc
if (z_w(i,j,it) > zero) then
dt_warm(i,j,it) = two*xt(i,j,it)/z_w(i,j,it)
end if
end do
end do
endif
! Deallocate local work arrays
deallocate(work)
call close_dataset(sfcges)
!
! Print date/time stamp
if ( mype == 0 ) write(6, &
'(a,'': sfc read,nlon,nlat= '',2i6,'',hour= '',f4.1,'',idate= '',4i5)') &
trim(my_name),lonb,latb,fhour,odate
! End of loop over time levels
end do
end subroutine read_sfc_
subroutine read_gfsncsfc_(iope,sfct,soil_moi,sno,soil_temp,veg_frac,fact10,sfc_rough, &
veg_type,soil_type,terrain,isli,use_sfc_any, &
tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d)
!$$$ subprogram documentation block
! . . . .
! subprogram: read_gfsncsfc_ read netcdf sfc hist file
! prgmmr: martin org: NCEP/EMC date: 2019-09-23
!
! abstract: read netcdf surface file
!
! program history log:
! 2019-09-23 Martin Initial version.
!
! input argument list:
! iope - mpi task handling i/o
! use_sfc_any - true if any processor uses extra surface fields
!
! output argument list:
! sfct - surface temperature (skin temp)
! soil_moi - soil moisture of first layer
! sno - snow depth
! soil_temp - soil temperature of first layer
! veg_frac - vegetation fraction
! fact10 - 10 meter wind factor
! sfc_rough - surface roughness
! veg_type - vegetation type
! soil_type - soil type
! terrain - terrain height
! isli - sea/land/ice mask
! tref - oceanic foundation temperature
! dt_cool - optional, sub-layer cooling amount at sub-skin layer
! z_c - optional, depth of sub-layer cooling layer
! dt_warm - optional, diurnal warming amount at sea surface
! z_w - optional, depth of diurnal warming layer
! c_0 - optional, coefficient to calculate d(Tz)/d(tf)
! c_d - optional, coefficient to calculate d(Tz)/d(tf)
! w_0 - optional, coefficient to calculate d(Tz)/d(tf)
! w_d - optional, coefficient to calculate d(Tz)/d(tf)
!
! attributes:
! language: f90
!
!$$$
use kinds, only: r_kind,i_kind,r_single
use gridmod, only: nlat_sfc,nlon_sfc
use guess_grids, only: nfldsfc,sfcmod_mm5,sfcmod_gfs
use mpimod, only: mpi_itype,mpi_rtype4,mpi_comm_world,mype
use constants, only: zero
implicit none
! Declare passed variables
integer(i_kind), intent(in) :: iope
logical, intent(in) :: use_sfc_any
real(r_single), dimension(nlat_sfc,nlon_sfc,nfldsfc), intent(out) :: sfct,soil_moi,sno,soil_temp,veg_frac,fact10,sfc_rough
real(r_single), dimension(nlat_sfc,nlon_sfc), intent(out) :: veg_type,soil_type,terrain
integer(i_kind), dimension(nlat_sfc,nlon_sfc), intent(out) :: isli
real(r_single), optional, dimension(nlat_sfc,nlon_sfc,nfldsfc), intent(out) :: tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d
! Declare local variables
integer(i_kind):: iret,npts,nptsall
!-----------------------------------------------------------------------------
! Read surface history file on processor iope
if(mype == iope)then
if ( present(tref) ) then
call read_sfc_(sfct,soil_moi,sno,soil_temp,veg_frac,fact10,sfc_rough, &
veg_type,soil_type,terrain,isli,use_sfc_any, &
tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d)
write(*,*) 'read_sfc netcdf, with NSST variables'
else
call read_sfc_(sfct,soil_moi,sno,soil_temp,veg_frac,fact10,sfc_rough, &
veg_type,soil_type,terrain,isli,use_sfc_any)
write(*,*) 'read_sfc netcdf, without NSST variables'
endif
endif
! Load onto all processors
npts=nlat_sfc*nlon_sfc
nptsall=npts*nfldsfc
call mpi_bcast(sfct, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
call mpi_bcast(fact10, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
call mpi_bcast(sno, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
if(sfcmod_mm5 .or. sfcmod_gfs)then
call mpi_bcast(sfc_rough, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
else
sfc_rough = zero
endif
call mpi_bcast(terrain, npts, mpi_rtype4,iope,mpi_comm_world,iret)
call mpi_bcast(isli, npts, mpi_itype, iope,mpi_comm_world,iret)
if(use_sfc_any)then
call mpi_bcast(veg_frac, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
call mpi_bcast(soil_temp,nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
call mpi_bcast(soil_moi, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
call mpi_bcast(veg_type, npts, mpi_rtype4,iope,mpi_comm_world,iret)
call mpi_bcast(soil_type,npts, mpi_rtype4,iope,mpi_comm_world,iret)
endif
if ( present(tref) ) then
call mpi_bcast(tref, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
call mpi_bcast(dt_cool, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
call mpi_bcast(z_c, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
call mpi_bcast(dt_warm, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
call mpi_bcast(z_w, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
call mpi_bcast(c_0, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
call mpi_bcast(c_d, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
call mpi_bcast(w_0, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
call mpi_bcast(w_d, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
endif
end subroutine read_gfsncsfc_
subroutine read_sfc_anl_(isli_anl)
!$$$ subprogram documentation block
! . . . .
! subprogram: read_sfc_anl_ read netcdf surface file with analysis resolution
!
! prgmmr: martin org: NCEP/EMC date: 2019-09-23
!
! abstract: read netcdf surface file at analysis grids when nlon /= nlon_sfc or nlat /= nlat_sfc
!
! program history log:
! 2019-09-23 Martin Initial version.
!
! input argument list:
!
! output argument list:
! isli - sea/land/ice mask
!
! attributes:
! language: f90
!
!$$$
use mpimod, only: mype
use kinds, only: r_kind,i_kind,r_single
use gridmod, only: nlat,nlon
use constants, only: zero
use module_ncio, only: Dataset, Variable, Dimension, open_dataset,&
close_dataset, get_dim, read_vardata, get_idate_from_time_units
implicit none
! Declare passed variables
integer(i_kind), dimension(nlat,nlon), intent( out) :: isli_anl
! Declare local parameters
integer(i_kind),dimension(6):: idate
integer(i_kind),dimension(4):: odate
! Declare local variables
character(len=24) :: filename
character(len=120) :: my_name = 'READ_GFSNCSFC_ANL'
integer(i_kind) :: i,j,iret
integer(i_kind) :: lonb, latb
real(r_single),allocatable, dimension(:) :: fhour
real(r_single), allocatable, dimension(:,:) :: work,outtmp
type(Dataset) :: sfcges
type(Dimension) :: ncdim
!-----------------------------------------------------------------------------
filename='sfcf06_anlgrid'
! open the netCDF file
sfcges = open_dataset(filename,errcode=iret)
if (iret/=0) then
write(6,*) trim(my_name),': ***FATAL ERROR*** ',trim(filename),' NOT AVAILABLE: PROGRAM STOPS'
call stop2(999)
endif
! get dimension sizes
ncdim = get_dim(sfcges, 'grid_xt'); lonb = ncdim%len
ncdim = get_dim(sfcges, 'grid_yt'); latb = ncdim%len
! get time information
idate = get_idate_from_time_units(sfcges)
odate(1) = idate(4) !hour
odate(2) = idate(2) !month
odate(3) = idate(3) !day
odate(4) = idate(1) !year
call read_vardata(sfcges, 'time', fhour) ! might need to change this to attribute later
! depends on model changes from Jeff Whitaker
if ( (latb /= nlat-2) .or. (lonb /= nlon) ) then
if ( mype == 0 ) write(6, &
'(a,'': inconsistent spatial dimension '',''nlon,nlatm2 = '',2(i4,tr1),''-vs- sfc file lonb,latb = '',i4)') &
trim(my_name),nlon,nlat-2,lonb,latb
call stop2(102)
endif
!
! Read the surface records (lonb, latb) and convert to GSI array pattern (nlat,nlon)
! Follow the read order sfcio in ncepgfs_io
!
allocate(work(lonb,latb))
work = zero
! slmsk
call read_vardata(sfcges, 'land', work)
allocate(outtmp(latb+2,lonb))
call tran_gfsncsfc(work,outtmp,lonb,latb)
do j=1,lonb
do i=1,latb+2
isli_anl(i,j) = nint(outtmp(i,j))
end do
end do
deallocate(outtmp)
! Deallocate local work arrays
deallocate(work)
call close_dataset(sfcges)
!
! Print date/time stamp
if ( mype == 0 ) write(6, &
'(a,'': read_sfc_anl_ ,nlon,nlat= '',2i6,'',hour= '',f4.1,'',idate= '',4i5)') &
trim(my_name),lonb,latb,fhour,odate
end subroutine read_sfc_anl_
subroutine read_gfsncsfc_anl_(iope,isli_anl)
!$$$ subprogram documentation block
! . . . .
! subprogram: read_gfsncsfc_anl read netcdf surface guess file with analysis resolution
!
! prgmmr: Martin org: NCEP/EMC date: 2019-09-23
!
! abstract: read netcdf surface file at analysis grids
!
! program history log:
! 2019-09-23 Martin Initial version.
!
! input argument list:
! iope - mpi task handling i/o
!
! output argument list:
! isli - sea/land/ice mask
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind,r_single
use gridmod, only: nlat,nlon
use mpimod, only: mpi_itype,mpi_comm_world,mype
implicit none
! Declare passed variables
integer(i_kind), intent(in ) :: iope
integer(i_kind), dimension(nlat,nlon), intent( out) :: isli_anl
! Declare local variables
integer(i_kind):: iret,npts
!-----------------------------------------------------------------------------
! Read surface file on processor iope
if(mype == iope)then
call read_sfc_anl_(isli_anl)
write(*,*) 'read_sfc netcdf'
end if
! Load onto all processors
npts=nlat*nlon
call mpi_bcast(isli_anl,npts,mpi_itype,iope,mpi_comm_world,iret)
end subroutine read_gfsncsfc_anl_
subroutine read_nst_ (tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d)
!$$$ subprogram documentation block
! . . . .
! subprogram: read_nst_ read netcdf nst surface guess file (quadratic
! Gaussin grids) without scattering to tasks
! prgmmr: Martin org: NCEP/EMC date: 2019-09-23
!
! abstract: read netcdf surface NST file
!
! program history log:
! 2019-09-23 Martin Initial version based on sub read_nemsnst
!
! input argument list:
!
! output argument list:
! tref (:,:) ! oceanic foundation temperature
! dt_cool (:,:) ! sub-layer cooling amount at sub-skin layer
! z_c (:,:) ! depth of sub-layer cooling layer
! dt_warm (:,:) ! diurnal warming amount at sea surface (skin layer)
! z_w (:,:) ! depth of diurnal warming layer
! c_0 (:,:) ! coefficient to calculate d(Tz)/d(tr)
! c_d (:,:) ! coefficient to calculate d(Tz)/d(tr)
! w_0 (:,:) ! coefficient to calculate d(Tz)/d(tr)
! w_d (:,:) ! coefficient to calculate d(Tz)/d(tr)
!
! attributes:
! language: f90
!
!$$$
use kinds, only: r_kind,i_kind,r_single
use mpimod, only: mype
use gridmod, only: nlat_sfc,nlon_sfc
use constants, only: zero,two
use guess_grids, only: nfldnst,ifilenst
use module_ncio, only: Dataset, Variable, Dimension, open_dataset,&
close_dataset, get_dim, read_vardata, get_idate_from_time_units
implicit none
! Declare passed variables
real(r_single) , dimension(nlat_sfc,nlon_sfc,nfldnst), intent( out) :: &
tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d
! Declare local parameters
integer(i_kind),parameter :: n_nst=9
integer(i_kind),dimension(6) :: idate
integer(i_kind),dimension(4) :: odate
! Declare local variables
character(len=6) :: filename
character(len=120) :: my_name = 'READ_GFSNCNST'
integer(i_kind) :: i,j,it,latb,lonb
real(r_single),allocatable, dimension(:) :: fhour
real(r_single), dimension(nlat_sfc,nlon_sfc,nfldnst) :: xt
real(r_single), allocatable, dimension(:,:) :: work
type(Dataset) :: sfcges
type(Dimension) :: ncdim
!-----------------------------------------------------------------------------
do it=1,nfldnst
! read a nst file on the task
write(filename,200)ifilenst(it)
200 format('nstf',i2.2)
! open the netCDF file
sfcges = open_dataset(filename)
! get dimension sizes
ncdim = get_dim(sfcges, 'grid_xt'); lonb = ncdim%len
ncdim = get_dim(sfcges, 'grid_yt'); latb = ncdim%len
! get time information
idate = get_idate_from_time_units(sfcges)
odate(1) = idate(4) !hour
odate(2) = idate(2) !month
odate(3) = idate(3) !day
odate(4) = idate(1) !year
call read_vardata(sfcges, 'time', fhour) ! might need to change this to attribute later
! depends on model changes from Jeff Whitaker
if ( (latb /= nlat_sfc-2) .or. (lonb /= nlon_sfc) ) then
if ( mype == 0 ) &
write(6,'(a,'': inconsistent spatial dimension nlon,nlatm2 = '',2(i4,tr1),''-vs- sfc file lonb,latb = '',i4)') &
trim(my_name),nlon_sfc,nlat_sfc-2,lonb,latb
call stop2(80)
endif
!
! Load surface fields into local work array
!
allocate(work(lonb,latb))
work = zero
! Tref
call read_vardata(sfcges, 'tref', work)
call tran_gfsncsfc(work,tref(1,1,it),lonb,latb)
! dt_cool
call read_vardata(sfcges, 'dtcool', work)
call tran_gfsncsfc(work,dt_cool(1,1,it),lonb,latb)
! z_c
call read_vardata(sfcges, 'zc', work)
call tran_gfsncsfc(work,z_c(1,1,it),lonb,latb)
! xt
call read_vardata(sfcges, 'xt', work)
call tran_gfsncsfc(work,xt(1,1,it),lonb,latb)
! xz
call read_vardata(sfcges, 'xz', work)
call tran_gfsncsfc(work,z_w(1,1,it),lonb,latb)
!
! c_0
call read_vardata(sfcges, 'c0', work)
call tran_gfsncsfc(work,c_0(1,1,it),lonb,latb)
! c_d
call read_vardata(sfcges, 'cd', work)
call tran_gfsncsfc(work,c_d(1,1,it),lonb,latb)
! w_0
call read_vardata(sfcges, 'w0', work)
call tran_gfsncsfc(work,w_0(1,1,it),lonb,latb)
! w_d
call read_vardata(sfcges, 'wd', work)
call tran_gfsncsfc(work,w_d(1,1,it),lonb,latb)
!
! Get diurnal warming amout at z=0
!
do j = 1,nlon_sfc
do i = 1,nlat_sfc
if (z_w(i,j,it) > zero) then
dt_warm(i,j,it) = two*xt(i,j,it)/z_w(i,j,it)
end if
end do
end do
! Deallocate local work arrays
deallocate(work)
call close_dataset(sfcges)
! End of loop over time levels
end do
end subroutine read_nst_
subroutine read_gfsncnst_ (iope,tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d)
!$$$ subprogram documentation block
! . . . .
! subprogram: read_gfsnc_nst
! prgmmr: Martin org: NCEP/EMC date: 2019-09-23
!
! abstract: read netcdf nst fields from a specific task and then broadcast to others
!
! input argument list:
! iope - mpi task handling i/o
!
! output argument list:
! tref (:,:) ! oceanic foundation temperature
! dt_cool (:,:) ! sub-layer cooling amount at sub-skin layer
! z_c (:,:) ! depth of sub-layer cooling layer
! dt_warm (:,:) ! diurnal warming amount at sea surface
! z_w (:,:) ! depth of diurnal warming layer
! c_0 (:,:) ! coefficient to calculate d(Tz)/d(tr) in dimensionless
! c_d (:,:) ! coefficient to calculate d(Tz)/d(tr) in m^-1
! w_0 (:,:) ! coefficient to calculate d(Tz)/d(tr) in dimensionless
! w_d (:,:) ! coefficient to calculate d(Tz)/d(tr) in m^-1
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind,r_single
use gridmod, only: nlat_sfc,nlon_sfc
use guess_grids, only: nfldnst
use mpimod, only: mpi_itype,mpi_rtype4,mpi_comm_world
use mpimod, only: mype
use constants, only: zero
implicit none
! Declare passed variables
integer(i_kind), intent(in ) :: iope
real(r_single), dimension(nlat_sfc,nlon_sfc,nfldnst), intent( out) :: &
tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d
! Declare local variables
integer(i_kind):: iret,npts,nptsall
!-----------------------------------------------------------------------------
! Read nst file on processor iope
if(mype == iope)then
write(*,*) 'read_nst netcdf'
call read_nst_(tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d)
end if
! Load onto all processors
npts=nlat_sfc*nlon_sfc
nptsall=npts*nfldnst
call mpi_bcast(tref, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
call mpi_bcast(dt_cool, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
call mpi_bcast(z_c, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
call mpi_bcast(dt_warm, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
call mpi_bcast(z_w, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
call mpi_bcast(c_0, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
call mpi_bcast(c_d, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
call mpi_bcast(w_0, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
call mpi_bcast(w_d, nptsall,mpi_rtype4,iope,mpi_comm_world,iret)
end subroutine read_gfsncnst_
subroutine write_atm_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin)
!$$$ subprogram documentation block
! . . .
! subprogram: write_ncatm --- Gather, transform, and write out to netcdf
!
! prgmmr: whitaker org: oar date: 2019-10-03
!
! abstract: This routine gathers fields needed for the GSI analysis
! file from subdomains and then transforms the fields from
! analysis grid to model guess grid, then written to an
! netcdf atmospheric analysis file.
!
! program history log:
! 2019-10-03 whitaker initial version
!
! input argument list:
! filename - file to open and write to
! mype_out - mpi task to write output file
! gfs_bundle - bundle containing fields on subdomains
! ibin - time bin
!
! 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
use constants, only: r1000,fv,one,zero,qcmin,r0_05,t0c
use mpimod, only: mpi_rtype
use mpimod, only: mpi_comm_world
use mpimod, only: ierror
use mpimod, only: mype
use guess_grids, only: ifilesig
use guess_grids, only: load_geop_hgt,geop_hgti,ges_geopi
use gridmod, only: ntracer
use gridmod, only: ncloud
use gridmod, only: strip,jcap_b,bk5
use general_commvars_mod, only: load_grid,fill2_ns,filluv2_ns
use general_specmod, only: spec_vars
use obsmod, only: iadate
use gsi_4dvar, only: ibdate,nhr_obsbin,lwrite4danl
use general_sub2grid_mod, only: sub2grid_info
use egrid2agrid_mod,only: g_egrid2agrid,g_create_egrid2agrid,egrid2agrid_parm,destroy_egrid2agrid
use constants, only: two,pi,half,deg2rad
use gsi_bundlemod, only: gsi_bundle
use gsi_bundlemod, only: gsi_bundlegetpointer
use cloud_efr_mod, only: cloud_calc_gfs
use netcdf, only: nf90_max_name
use module_ncio, only: open_dataset, close_dataset, Dimension, Dataset,&
read_attribute, write_attribute,get_dim, create_dataset, write_vardata, read_vardata,&
get_idate_from_time_units,quantize_data,get_time_units_from_idate,has_attr,has_var
use ncepnems_io, only: error_msg
implicit none
! !INPUT PARAMETERS:
type(sub2grid_info), intent(in) :: grd
type(spec_vars), intent(in) :: sp_a
character(len=24), intent(in) :: filename ! file to open and write to
integer(i_kind), intent(in) :: mype_out ! mpi task to write output file
type(gsi_bundle), intent(in) :: gfs_bundle
integer(i_kind), intent(in) :: ibin ! time bin
!-------------------------------------------------------------------------
real(r_kind),parameter:: r0_001 = 0.001_r_kind
character(6):: fname_ges
character(len=120) :: my_name = 'WRITE_GFSNCATM'
character(len=1) :: null = ' '
integer(i_kind),dimension(6):: idate,jdate
integer(i_kind) :: k, mm1, nlatm2, nord_int, i, j, kk, kr, nbits
integer(i_kind) :: iret, lonb, latb, levs, istatus
integer(i_kind) :: nfhour
integer(i_kind) :: istop = 104
integer(i_kind),dimension(5):: mydate
integer(i_kind),dimension(8) :: ida,jda
real(r_kind),dimension(5) :: fha
real(r_kind), allocatable, dimension(:) :: fhour
real(r_kind),allocatable,dimension(:) :: rlats_tmp,rlons_tmp
real(r_kind),pointer,dimension(:,:) :: sub_ps
real(r_kind),pointer,dimension(:,:,:) :: sub_u,sub_v,sub_tv
real(r_kind),pointer,dimension(:,:,:) :: sub_q,sub_oz,sub_cwmr
real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig) :: sub_dzb,sub_dza
real(r_kind),dimension(grd%lat1*grd%lon1) :: psm
real(r_kind),dimension(grd%lat1*grd%lon1,grd%nsig):: tvsm, usm, vsm
real(r_kind),dimension(grd%lat1*grd%lon1,grd%nsig):: qsm, ozsm
real(r_kind),dimension(grd%lat1*grd%lon1,grd%nsig):: cwsm, dzsm
real(r_kind),dimension(max(grd%iglobal,grd%itotsub)) :: work1,work2
real(r_kind),dimension(grd%nlon,grd%nlat-2):: grid
real(r_kind),allocatable,dimension(:) :: rlats,rlons,clons,slons
real(r_kind),allocatable,dimension(:,:) :: grid_b,grid_b2
real(r_kind),allocatable,dimension(:,:,:) :: grid_c, grid3, grid_c2, grid3b
real(4), allocatable, dimension(:,:) :: values_2d,values_2d_tmp
real(4), allocatable, dimension(:,:,:) :: values_3d,values_3d_tmp,ug3d,vg3d
real(4) compress_err
type(Dataset) :: atmges, atmanl
type(Dimension) :: ncdim
character(len=nf90_max_name) :: time_units
logical diff_res,eqspace
logical,dimension(1) :: vector
type(egrid2agrid_parm) :: p_low,p_high
!*************************************************************************
! Initialize local variables
mm1=mype+1
nlatm2=grd%nlat-2
diff_res=.false.
istatus=0
call gsi_bundlegetpointer(gfs_bundle,'ps', sub_ps, iret); istatus=istatus+iret
call gsi_bundlegetpointer(gfs_bundle,'u', sub_u, iret); istatus=istatus+iret
call gsi_bundlegetpointer(gfs_bundle,'v', sub_v, 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,*) 'write_atm_: ERROR'
write(6,*) 'Missing some of the required fields'
write(6,*) 'Aborting ... '
endif
call stop2(999)
endif
if ( sp_a%jcap /= jcap_b ) then
if ( mype == 0 ) write(6, &
'('' dual resolution for nems sp_a%jcap,jcap_b = '',2i6)') &
sp_a%jcap,jcap_b
diff_res = .true.
endif
! Single task writes analysis data to analysis file
if ( mype == mype_out ) then
write(fname_ges,'(''sigf'',i2.2)') ifilesig(ibin)
! open the netCDF file
atmges = open_dataset(fname_ges,errcode=iret)
if ( iret /= 0 ) call error_msg(trim(my_name),null,null,'open',istop,iret)
! get dimension sizes
ncdim = get_dim(atmges, 'grid_xt'); lonb = ncdim%len
ncdim = get_dim(atmges, 'grid_yt'); latb = ncdim%len
ncdim = get_dim(atmges, 'pfull'); levs = ncdim%len
if ( levs /= grd%nsig ) then
write(6,*) trim(my_name),': problem in data dimension background levs = ',levs,' nsig = ',grd%nsig
call stop2(103)
endif
! get time information
idate = get_idate_from_time_units(atmges)
call read_vardata(atmges, 'time', fhour) ! might need to change this to attribute later
! depends on model changes from Jeff Whitaker
nfhour = fhour(1)
atmanl = create_dataset(filename, atmges, &
copy_vardata=.true., errcode=iret)
if ( iret /= 0 ) call error_msg(trim(my_name),trim(filename),null,'open',istop,iret)
! Update time information (with ibdate) and write it to analysis file
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)
jdate(1) = jda(1) ! analysis year
jdate(2) = jda(2) ! analysis month
jdate(3) = jda(3) ! analysis day
jdate(4) = jda(5) ! analysis hour
jdate(5) = iadate(5) ! analysis minute
jdate(6) = 0 ! analysis second
fhour = zero
call write_vardata(atmanl, 'time', fhour)
time_units = get_time_units_from_idate(jdate)
call write_attribute(atmanl, 'units', time_units, 'time')
! Allocate structure arrays to hold data
allocate(values_3d_tmp(lonb,latb,levs),values_2d_tmp(lonb,latb))
allocate(grid3b(grd%nlat,grd%nlon,1))
allocate( grid_b(lonb,latb),grid_c(latb+2,lonb,1),grid3(grd%nlat,grd%nlon,1))
allocate( grid_b2(lonb,latb),grid_c2(latb+2,lonb,1))
allocate( rlats(latb+2),rlons(lonb),clons(lonb),slons(lonb))
call read_vardata(atmges, 'grid_xt', rlons_tmp, errcode=iret)
call read_vardata(atmges, 'grid_yt', rlats_tmp, errcode=iret)
do j=1,latb
rlats(latb+2-j)=deg2rad*rlats_tmp(j)
enddo
rlats(1)=-half*pi
rlats(latb+2)=half*pi
do j=1,lonb
rlons(j)=deg2rad*rlons_tmp(j)
enddo
deallocate(rlons_tmp, rlats_tmp)
do j=1,lonb
clons(j)=cos(rlons(j))
slons(j)=sin(rlons(j))
enddo
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_low,.false.,eqspace=eqspace)
call g_create_egrid2agrid(latb+2,rlats,lonb,rlons, &
grd%nlat,sp_a%rlats,grd%nlon,sp_a%rlons,&
nord_int,p_high,.false.,eqspace=eqspace)
deallocate(rlats,rlons)
endif ! if ( mype == mype_out )
! compute delz (so that delz < 0 as model expects)
do k=1,grd%nsig
sub_dzb(:,:,k) = ges_geopi(:,:,k,ibin) - ges_geopi(:,:,k+1,ibin)
enddo
if ((.not. lwrite4danl) .or. ibin == 1) call load_geop_hgt
do k=1,grd%nsig
sub_dza(:,:,k) = geop_hgti(:,:,k,ibin) - geop_hgti(:,:,k+1,ibin)
enddo
sub_dza = sub_dza - sub_dzb !sub_dza is increment
! Strip off boundary points from subdomains
call strip(sub_ps ,psm)
call strip(sub_tv ,tvsm ,grd%nsig)
call strip(sub_q ,qsm ,grd%nsig)
call strip(sub_oz ,ozsm ,grd%nsig)
call strip(sub_cwmr,cwsm ,grd%nsig)
call strip(sub_u ,usm ,grd%nsig)
call strip(sub_v ,vsm ,grd%nsig)
call strip(sub_dza ,dzsm ,grd%nsig)
! Thermodynamic variable
! The GSI analysis variable is virtual temperature (Tv). For NEMSIO
! output we need the sensible temperature.
! Convert Tv to T
tvsm = tvsm/(one+fv*qsm)
! Generate and write analysis fields
! Surface pressure and delp.
call mpi_gatherv(psm,grd%ijn(mm1),mpi_rtype,&
work1,grd%ijn,grd%displs_g,mpi_rtype,&
mype_out,mpi_comm_world,ierror)
if (mype==mype_out) then
call read_vardata(atmges,'pressfc',values_2d,errcode=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),'pres','read',istop,iret)
grid_b = r0_001*values_2d
vector(1)=.false.
call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector)
do kk=1,grd%iglobal
i=grd%ltosi(kk)
j=grd%ltosj(kk)
grid3(i,j,1)=work1(kk)-grid3(i,j,1)
work1(kk)=grid3(i,j,1)
end do
if (has_var(atmges,'dpres')) then ! skip this if delp not in guess file.
call read_vardata(atmges,'dpres',values_3d,errcode=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),'dpres','read',istop,iret)
do k=1,grd%nsig
kr = grd%nsig-k+1
do kk=1,grd%iglobal
i=grd%ltosi(kk)
j=grd%ltosj(kk)
grid3(i,j,1)=work1(kk)*(bk5(k)-bk5(k+1))
enddo
call g_egrid2agrid(p_high,grid3,grid_c2,1,1,vector)
do j=1,latb
do i=1,lonb
values_3d(i,j,kr)=values_3d(i,j,kr)+r1000*(grid_c2(latb-j+2,i,1))
enddo
enddo
enddo
if (has_attr(atmges, 'nbits', 'dpres')) then
call read_attribute(atmges, 'nbits', nbits, 'dpres')
values_3d_tmp = values_3d
call quantize_data(values_3d_tmp, values_3d, nbits, compress_err)
call write_attribute(atmanl,&
'max_abs_compression_error',compress_err,'dpres')
endif
call write_vardata(atmanl,'dpres',values_3d,errcode=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),'dpres','write',istop,iret)
endif
do kk=1,grd%iglobal
i=grd%ltosi(kk)
j=grd%ltosj(kk)
grid3(i,j,1)=work1(kk)
enddo
call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector)
do j=1,latb
do i=1,lonb
grid_b(i,j)=r1000*(grid_b(i,j)+grid_c(latb-j+2,i,1))
end do
end do
values_2d = grid_b
if (has_attr(atmges, 'nbits', 'pressfc')) then
call read_attribute(atmges, 'nbits', nbits, 'pressfc')
values_2d_tmp = values_2d
call quantize_data(values_2d_tmp, values_2d, nbits, compress_err)
call write_attribute(atmanl,&
'max_abs_compression_error',compress_err,'pressfc')
endif
call write_vardata(atmanl,'pressfc',values_2d,errcode=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),'pressfc','write',istop,iret)
endif
! u, v
if (mype==mype_out) then
if (allocated(values_3d)) deallocate(values_3d)
call read_vardata(atmges, 'ugrd', ug3d, errcode=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ugrd','read',istop,iret)
call read_vardata(atmges, 'vgrd', vg3d, errcode=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),'vgrd','read',istop,iret)
endif
do k=1,grd%nsig
kr = grd%nsig-k+1
call mpi_gatherv(usm(1,k),grd%ijn(mm1),mpi_rtype,&
work1,grd%ijn,grd%displs_g,mpi_rtype,&
mype_out,mpi_comm_world,ierror)
call mpi_gatherv(vsm(1,k),grd%ijn(mm1),mpi_rtype,&
work2,grd%ijn,grd%displs_g,mpi_rtype,&
mype_out,mpi_comm_world,ierror)
if (mype==mype_out) then
if(diff_res)then
grid_b = ug3d(:,:,kr)
grid_b2 = vg3d(:,:,kr)
vector(1)=.true.
call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons)
call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector)
do kk=1,grd%iglobal
i=grd%ltosi(kk)
j=grd%ltosj(kk)
grid3(i,j,1)=work1(kk)-grid3(i,j,1)
end do
call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector)
do j=1,latb
do i=1,lonb
grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1)
end do
end do
call g_egrid2agrid(p_low,grid_c2,grid3,1,1,vector)
do kk=1,grd%iglobal
i=grd%ltosi(kk)
j=grd%ltosj(kk)
grid3(i,j,1)=work2(kk)-grid3(i,j,1)
end do
call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector)
do j=1,latb
do i=1,lonb
grid_b2(i,j)=grid_b2(i,j)+grid_c(latb-j+2,i,1)
end do
end do
ug3d(:,:,kr) = grid_b
vg3d(:,:,kr) = grid_b2
else
call load_grid(work1,grid)
ug3d(:,:,kr) = grid
call load_grid(work2,grid)
vg3d(:,:,kr) = grid
end if
endif ! mype_out
end do
! Zonal wind
if (mype==mype_out) then
if (has_attr(atmges, 'nbits', 'ugrd')) then
call read_attribute(atmges, 'nbits', nbits, 'ugrd')
values_3d_tmp = ug3d
call quantize_data(values_3d_tmp, ug3d, nbits, compress_err)
call write_attribute(atmanl,&
'max_abs_compression_error',compress_err,'ugrd')
endif
call write_vardata(atmanl,'ugrd',ug3d,errcode=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),'ugrd','write',istop,iret)
! Meridional wind
if (has_attr(atmges, 'nbits', 'vgrd')) then
call read_attribute(atmges, 'nbits', nbits, 'vgrd')
values_3d_tmp = vg3d
call quantize_data(values_3d_tmp, vg3d, nbits, compress_err)
call write_attribute(atmanl,&
'max_abs_compression_error',compress_err,'vgrd')
endif
call write_vardata(atmanl,'vgrd',vg3d,errcode=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),'vgrd','write',istop,iret)
deallocate(ug3d, vg3d)
endif
! Thermodynamic variable
if (mype==mype_out) then
call read_vardata(atmges, 'tmp', values_3d, errcode=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),'tmp','read',istop,iret)
endif
do k=1,grd%nsig
kr = grd%nsig-k+1
call mpi_gatherv(tvsm(1,k),grd%ijn(mm1),mpi_rtype,&
work1,grd%ijn,grd%displs_g,mpi_rtype,&
mype_out,mpi_comm_world,ierror)
if (mype == mype_out) then
if(diff_res)then
grid_b=values_3d(:,:,kr)
vector(1)=.false.
call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector)
do kk=1,grd%iglobal
i=grd%ltosi(kk)
j=grd%ltosj(kk)
grid3(i,j,1)=work1(kk)-grid3(i,j,1)
end do
call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector)
do j=1,latb
do i=1,lonb
grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1)
end do
end do
values_3d(:,:,kr) = grid_b
else
call load_grid(work1,grid)
values_3d(:,:,kr) = grid
end if
endif
end do
if (mype==mype_out) then
if (has_attr(atmges, 'nbits', 'tmp')) then
call read_attribute(atmges, 'nbits', nbits, 'tmp')
values_3d_tmp = values_3d
call quantize_data(values_3d_tmp, values_3d, nbits, compress_err)
call write_attribute(atmanl,&
'max_abs_compression_error',compress_err,'tmp')
endif
call write_vardata(atmanl,'tmp',values_3d,errcode=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),'tmp','write',istop,iret)
endif
! Specific humidity
if (mype==mype_out) then
call read_vardata(atmges, 'spfh', values_3d, errcode=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),'spfh','read',istop,iret)
endif
do k=1,grd%nsig
kr = grd%nsig-k+1
call mpi_gatherv(qsm(1,k),grd%ijn(mm1),mpi_rtype,&
work1,grd%ijn,grd%displs_g,mpi_rtype,&
mype_out,mpi_comm_world,ierror)
if (mype == mype_out) then
if(diff_res)then
grid_b=values_3d(:,:,kr)
vector(1)=.false.
call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector)
do kk=1,grd%iglobal
i=grd%ltosi(kk)
j=grd%ltosj(kk)
grid3(i,j,1)=work1(kk)-grid3(i,j,1)
end do
call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector)
do j=1,latb
do i=1,lonb
grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1)
end do
end do
values_3d(:,:,kr) = grid_b
else
call load_grid(work1,grid)
values_3d(:,:,kr) = grid
end if
endif
end do
if (mype==mype_out) then
if (has_attr(atmges, 'nbits', 'spfh')) then
call read_attribute(atmges, 'nbits', nbits, 'spfh')
values_3d_tmp = values_3d
call quantize_data(values_3d_tmp, values_3d, nbits, compress_err)
call write_attribute(atmanl,&
'max_abs_compression_error',compress_err,'spfh')
endif
call write_vardata(atmanl,'spfh',values_3d,errcode=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),'spfh','write',istop,iret)
endif
! Ozone
if (mype==mype_out) then
call read_vardata(atmges, 'o3mr', values_3d, errcode=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),'o3mr','read',istop,iret)
endif
do k=1,grd%nsig
kr = grd%nsig-k+1
call mpi_gatherv(ozsm(1,k),grd%ijn(mm1),mpi_rtype,&
work1,grd%ijn,grd%displs_g,mpi_rtype,&
mype_out,mpi_comm_world,ierror)
if (mype == mype_out) then
if(diff_res)then
grid_b=values_3d(:,:,kr)
vector(1)=.false.
call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector)
do kk=1,grd%iglobal
i=grd%ltosi(kk)
j=grd%ltosj(kk)
grid3(i,j,1)=work1(kk)-grid3(i,j,1)
end do
call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector)
do j=1,latb
do i=1,lonb
grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1)
end do
end do
values_3d(:,:,kr) = grid_b
else
call load_grid(work1,grid)
values_3d(:,:,kr) = grid
end if
endif
end do
if (mype==mype_out) then
if (has_attr(atmges, 'nbits', 'o3mr')) then
call read_attribute(atmges, 'nbits', nbits, 'o3mr')
values_3d_tmp = values_3d
call quantize_data(values_3d_tmp, values_3d, nbits, compress_err)
call write_attribute(atmanl,&
'max_abs_compression_error',compress_err,'o3mr')
endif
call write_vardata(atmanl,'o3mr',values_3d,errcode=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),'o3mr','write',istop,iret)
endif
! Cloud condensate mixing ratio
if (ntracer>2 .or. ncloud>=1) then
if (mype==mype_out) then
if (allocated(values_3d)) deallocate(values_3d)
call read_vardata(atmges, 'clwmr', ug3d, errcode=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),'clwmr','read',istop,iret)
call read_vardata(atmges, 'icmr', vg3d, errcode=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),'icmr','read',istop,iret)
endif
do k=1,grd%nsig
kr = grd%nsig-k+1
call mpi_gatherv(cwsm(1,k),grd%ijn(mm1),mpi_rtype,&
work1,grd%ijn,grd%displs_g,mpi_rtype,&
mype_out,mpi_comm_world,ierror)
call mpi_gatherv(tvsm(1,k),grd%ijn(mm1),mpi_rtype,&
work2,grd%ijn,grd%displs_g,mpi_rtype,&
mype_out,mpi_comm_world,ierror)
if (mype == mype_out) then
grid_b = ug3d(:,:,kr)
grid_b2=vg3d(:,:,kr)
grid_b = grid_b + grid_b2
vector(1)=.false.
call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
call g_egrid2agrid(p_low,grid_c,grid3,1,1,vector)
do kk=1,grd%iglobal
i=grd%ltosi(kk)
j=grd%ltosj(kk)
grid3(i,j,1)=work1(kk)-max(grid3(i,j,1),qcmin)
work2(kk) = -r0_05*(work2(kk) - t0c)
work2(kk) = max(zero,work2(kk))
work2(kk) = min(one,work2(kk))
grid3b(i,j,1)=grid3(i,j,1)
grid3(i,j,1)=grid3b(i,j,1)*(one - work2(kk))
end do
call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector)
grid_b = grid_b - grid_b2
do j=1,latb
do i=1,lonb
grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1)
end do
end do
ug3d(:,:,kr) = grid_b
do kk=1,grd%iglobal
i=grd%ltosi(kk)
j=grd%ltosj(kk)
grid3(i,j,1)=grid3b(i,j,1)*work2(kk)
end do
call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector)
do j=1,latb
do i=1,lonb
grid_b2(i,j)=grid_b2(i,j)+grid_c(latb-j+2,i,1)
end do
end do
vg3d(:,:,kr) = grid_b2
endif !mype == mype_out
end do
if (mype==mype_out) then
if (has_attr(atmges, 'nbits', 'clwmr')) then
call read_attribute(atmges, 'nbits', nbits, 'clwmr')
values_3d_tmp = ug3d
call quantize_data(values_3d_tmp, ug3d, nbits, compress_err)
call write_attribute(atmanl,&
'max_abs_compression_error',compress_err,'clwmr')
endif
call write_vardata(atmanl,'clwmr',ug3d,errcode=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),'clwmr','write',istop,iret)
if (has_attr(atmges, 'nbits', 'icmr')) then
call read_attribute(atmges, 'nbits', nbits, 'icmr')
values_3d_tmp = vg3d
call quantize_data(values_3d_tmp, vg3d, nbits, compress_err)
call write_attribute(atmanl,&
'max_abs_compression_error',compress_err,'icmr')
endif
call write_vardata(atmanl,'icmr',vg3d,errcode=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),'icmr','write',istop,iret)
deallocate(ug3d, vg3d)
endif
endif !ntracer
! Variables needed by the Unified Post Processor (dzdt, delz, delp)
if (mype==mype_out) then
if (has_var(atmges,'delz')) then ! if delz in guess file, read it.
call read_vardata(atmges, 'delz', values_3d, errcode=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),'delz','read',istop,iret)
endif
endif
do k=1,grd%nsig
kr = grd%nsig-k+1
call mpi_gatherv(dzsm(1,k),grd%ijn(mm1),mpi_rtype,&
work1,grd%ijn,grd%displs_g,mpi_rtype,&
mype_out,mpi_comm_world,ierror)
if (mype == mype_out) then
if (has_var(atmges,'delz')) then
if(diff_res)then
grid_b=values_3d(:,:,kr)
do kk=1,grd%iglobal
i=grd%ltosi(kk)
j=grd%ltosj(kk)
grid3(i,j,1)=work1(kk)
end do
call g_egrid2agrid(p_high,grid3,grid_c,1,1,vector)
do j=1,latb
do i=1,lonb
grid_b(i,j)=grid_b(i,j)+grid_c(latb-j+2,i,1)
end do
end do
values_3d(:,:,kr) = grid_b
else
call load_grid(work1,grid)
values_3d(:,:,kr) = values_3d(:,:,kr) + grid
end if
end if
endif
end do
if (mype==mype_out) then
! if delz in guess file, write to analysis file.
if (has_var(atmges,'delz')) then
if (has_attr(atmges, 'nbits', 'delz')) then
call read_attribute(atmges, 'nbits', nbits, 'delz')
values_3d_tmp = values_3d
call quantize_data(values_3d_tmp, values_3d, nbits, compress_err)
call write_attribute(atmanl,&
'max_abs_compression_error',compress_err,'delz')
endif
call write_vardata(atmanl,'delz',values_3d,errcode=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),'delz','write',istop,iret)
endif
endif
!
! Deallocate local array
!
if (mype==mype_out) then
deallocate(grid_b,grid_b2,grid_c,grid_c2,grid3,clons,slons)
deallocate(grid3b)
call close_dataset(atmges, errcode=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(fname_ges),null,'close',istop,iret)
call close_dataset(atmanl, errcode=iret)
if (iret /= 0) call error_msg(trim(my_name),trim(filename),null,'close',istop,iret)
!
! Deallocate local array
!
if (allocated(values_2d)) deallocate(values_2d)
if (allocated(values_3d)) deallocate(values_3d)
if (allocated(values_2d_tmp)) deallocate(values_2d_tmp)
if (allocated(values_3d_tmp)) deallocate(values_3d_tmp)
!
write(6,'(a,'': atm anal written for lonb,latb,levs= '',3i6,'',valid hour= '',f4.1,'',idate= '',4i5)') &
trim(my_name),lonb,latb,levs,fhour(1),jdate(1:4)
endif
end subroutine write_atm_
subroutine write_sfc_ (filename,mype_sfc,dsfct)
!$$$ subprogram documentation block
! . . .
! subprogram: write_gfsncsfc --- Write surface analysis to file
!
! prgmmr: Martin org: NCEP/EMC date: 2019-09-24
!
! abstract: This routine writes the updated surface analysis.
!
! 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.
!
! program history log:
! 2019-09-24 Martin Initial version. Based on write_nemssfc
!
! input argument list:
! filename - file to open and write to
! dsfct - delta skin temperature
! mype_sfc - mpi task to write output file
!
! output argument list:
!
! attributes:
! language: f90
!
!$$$ end documentation block
! !USES:
use kinds, only: r_kind,i_kind,r_single
use mpimod, only: mpi_rtype
use mpimod, only: mpi_comm_world
use mpimod, only: ierror
use mpimod, only: mype
use gridmod, only: nlat,nlon
use gridmod, only: lat1,lon1
use gridmod, only: lat2,lon2
use gridmod, only: iglobal
use gridmod, only: ijn
use gridmod, only: displs_g
use gridmod, only: itotsub
use gridmod, only: rlats,rlons,rlats_sfc,rlons_sfc
use general_commvars_mod, only: ltosi,ltosj
use obsmod, only: iadate
use constants, only: zero
use netcdf, only: nf90_max_name
use module_ncio, only: open_dataset, close_dataset, Dimension, Dataset,&
get_dim, create_dataset, write_vardata, read_vardata,&
get_time_units_from_idate, write_attribute
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_sfc ! mpi task to write output file
! !OUTPUT PARAMETERS:
!-------------------------------------------------------------------------
! Declare local parameters
character( 6),parameter:: fname_ges='sfcf06'
! Declare local variables
character(len=120) :: my_name = 'WRITE_GFSNCSFC'
integer(i_kind),dimension(6):: jdate
integer(i_kind),dimension(4):: odate
integer(i_kind) :: i, j, ip1, jp1, ilat, ilon, jj, mm1
integer(i_kind) :: nlatm2, lonb, latb
real(r_kind),allocatable,dimension(:) :: fhour
real(r_kind),dimension(lat1,lon1):: sfcsub
real(r_kind),dimension(nlon,nlat):: grid
real(r_kind),dimension(max(iglobal,itotsub)):: sfcall
real(r_kind),allocatable,dimension(:,:) :: tsea
real(r_single),dimension(nlon,nlat):: buffer
real(r_single),allocatable,dimension(:,:) :: buffer2,grid2
type(Dataset) :: sfcges,sfcanl
type(Dimension) :: ncdim
character(len=nf90_max_name) :: time_units
!*****************************************************************************
! Initialize local variables
mm1=mype+1
nlatm2=nlat-2
! Gather skin temperature information from all tasks.
do j=1,lon1
jp1 = j+1
do i=1,lat1
ip1 = i+1
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+1
do i=1,nlon
buffer(i,j)=grid(i,jj)
end do
end do
! Read surface guess file
! open the netCDF file
sfcges = open_dataset(fname_ges)
! get dimension sizes
ncdim = get_dim(sfcges, 'grid_xt'); lonb = ncdim%len
ncdim = get_dim(sfcges, 'grid_yt'); latb = ncdim%len
!
! Start to write output sfc file : filename
!
! First copy entire data from fname_ges to filename, then do selective update
!
sfcanl = create_dataset(filename, sfcges, copy_vardata=.true.)
!
! Replace header record date with analysis time from iadate
!
jdate(1) = iadate(1) ! analysis year
jdate(2) = iadate(2) ! analysis month
jdate(3) = iadate(3) ! analysis day
jdate(4) = iadate(4) ! analysis hour
jdate(5) = iadate(5) ! analysis minute
jdate(5) = 0 ! analysis minute
jdate(6) = 0 ! analysis scaled seconds
fhour = zero
odate(1) = jdate(4) !hour
odate(2) = jdate(2) !month
odate(3) = jdate(3) !day
odate(4) = jdate(1) !year
call write_vardata(sfcanl, 'time', fhour)
time_units = get_time_units_from_idate(jdate)
call write_attribute(sfcanl, 'units', time_units, 'time')
allocate(buffer2(lonb,latb))
allocate(grid2(lonb,latb))
allocate(tsea(lonb,latb))
!
! Only sea surface temperature will be updated in the SFC files
!
call read_vardata(sfcges, 'tmpsfc', tsea)
if ( (latb /= nlatm2) .or. (lonb /= nlon) ) then
write(6,*)trim(my_name),': different grid dimensions analysis', &
' vs sfc. interpolating sfc temperature nlon,nlat-2=',nlon, &
nlatm2,' -vs- sfc file lonb,latb=',lonb,latb
call intrp22(buffer, rlons,rlats,nlon,nlat, &
buffer2,rlons_sfc,rlats_sfc,lonb,latb)
else
do j=1,latb
do i=1,lonb
buffer2(i,j)=buffer(i,j+1)
end do
end do
endif
grid2 = tsea + buffer2
deallocate(buffer2)
! update tsea record
call write_vardata(sfcanl, 'tmpsfc', grid2)
call close_dataset(sfcges)
call close_dataset(sfcanl)
write(6,'(a,'': sfc anal written for lonb,latb= '',2i6,'',valid hour= '',f4.1,'',idate= '',4i5)') &
trim(my_name),lonb,latb,fhour,odate
endif
end subroutine write_sfc_
subroutine write_sfc_nst_ (mype_so,dsfct)
!$$$ subprogram documentation block
! . . .
! subprogram: write_sfc_nst --- Write both sfc and nst surface analysis to file
!
! prgmmr: Martin org: NCEP/EMC date: 2019-09-24
!
! abstract: This routine writes the sfc & nst analysis files and is nst_gsi dependent.
! Tr (foundation temperature), instead of skin temperature, is the analysis variable.
! nst_gsi > 2: Tr analysis is on
! nst_gsi <= 2: Tr analysis is off
!
! The routine gathers Tr field from subdomains,
! reformats the data records, and then writes each record
! to the output files.
!
! Since the gsi only update the Tr temperature, all
! other fields in surface are simply read from the guess
! files and written to the analysis file.
!
! program history log:
! 2019-09-24 Martin initial version based on routine write_nems_sfc_nst
!
! input argument list:
! dsfct - delta skin temperature
! mype_so - mpi task to write output file
!
! output argument list:
!
! attributes:
! language: f90
!
!$$$ end documentation block
! !USES:
use kinds, only: r_kind,i_kind,r_single
use mpimod, only: mpi_rtype,mpi_itype
use mpimod, only: mpi_comm_world
use mpimod, only: ierror
use mpimod, only: mype
use gridmod, only: nlat,nlon
use gridmod, only: lat1,lon1
use gridmod, only: lat2,lon2
use gridmod, only: nlat_sfc,nlon_sfc
use gridmod, only: iglobal
use gridmod, only: ijn
use gridmod, only: displs_g
use gridmod, only: itotsub
use general_commvars_mod, only: ltosi,ltosj
use obsmod, only: iadate
use constants, only: zero,two,tfrozen,z_w_max
use constants, only: zero_single
use guess_grids, only: isli2
use gsi_nstcouplermod, only: nst_gsi,zsea1,zsea2
use gridmod, only: rlats,rlons,rlats_sfc,rlons_sfc
use module_ncio, only: open_dataset, close_dataset, Dimension, Dataset,&
get_dim, create_dataset, write_vardata, read_vardata,&
get_time_units_from_idate, write_attribute
use netcdf, only: nf90_max_name
implicit none
! !INPUT PARAMETERS:
real(r_kind),dimension(lat2,lon2),intent(in ) :: dsfct ! delta skin temperature
integer(i_kind) ,intent(in ) :: mype_so ! mpi task to write output file
! !OUTPUT PARAMETERS:
!-------------------------------------------------------------------------
! Declare local parameters
character(6), parameter:: fname_sfcges = 'sfcf06'
character(6), parameter:: fname_sfcgcy = 'sfcgcy'
character(6), parameter:: fname_sfctsk = 'sfctsk'
character(6), parameter:: fname_sfcanl = 'sfcanl'
character(6), parameter:: fname_nstges = 'nstf06'
character(6), parameter:: fname_nstanl = 'nstanl'
character(6), parameter:: fname_dtfanl = 'dtfanl'
! Declare local variables
integer(i_kind), parameter:: io_dtfanl = 54
integer(i_kind), parameter:: nprep=15
real(r_kind),parameter :: houra = zero_single
character(len=120) :: my_name = 'WRITE_SFC_NST'
integer(i_kind),dimension(7):: jdate
integer(i_kind),dimension(4):: odate
integer(i_kind) :: i, j, ip1, jp1, ilat, ilon, mm1
integer(i_kind) :: lonb, latb, nlatm2
integer(i_kind) :: lonb_nst, latb_nst
real(r_kind),allocatable,dimension(:) :: fhour
real(r_single) :: r_zsea1,r_zsea2
real(r_kind), dimension(lat1,lon1):: dsfct_sub
integer(i_kind), dimension(lat1,lon1):: isli_sub
real(r_kind), dimension(max(iglobal,itotsub)):: dsfct_all
integer(i_kind), dimension(max(iglobal,itotsub)):: isli_all
real(r_kind), dimension(nlat,nlon):: dsfct_glb,dsfct_tmp
integer(i_kind), dimension(nlat,nlon):: isli_glb,isli_tmp
real(r_kind), dimension(nlat_sfc,nlon_sfc) :: dsfct_gsi
integer(i_kind), dimension(nlat_sfc,nlon_sfc) :: isli_gsi
real(r_kind), dimension(nlon_sfc,nlat_sfc-2):: dsfct_anl
real(r_single), dimension(nlon_sfc,nlat_sfc-2):: dtzm
real(r_single), dimension(nlat_sfc,nlon_sfc) :: work
real(r_single), allocatable, dimension(:,:) :: tsea,xt,xs,xu,xv,xz,zm,xtts,xzts,dt_cool,z_c, &
c_0,c_d,w_0,w_d,d_conv,ifd,tref,qrain
real(r_single), allocatable, dimension(:,:) :: slmsk_ges,slmsk_anl
type(Dataset) :: sfcges,sfcgcy,nstges,sfctsk,sfcanl,nstanl
type(Dimension) :: ncdim
character(len=nf90_max_name) :: time_units
!*****************************************************************************
! Initialize local variables
mm1=mype+1
nlatm2=nlat-2
!
! Extract the analysis increment and surface mask in subdomain without the buffer
!
do j=1,lon1
jp1 = j+1
do i=1,lat1
ip1 = i+1
dsfct_sub(i,j) = dsfct(ip1,jp1)
isli_sub (i,j) = isli2(ip1,jp1)
end do
end do
!
! Gather global analysis increment and surface mask info from subdomains
!
call mpi_gatherv(dsfct_sub,ijn(mm1),mpi_rtype,&
dsfct_all,ijn,displs_g,mpi_rtype,mype_so ,&
mpi_comm_world,ierror)
call mpi_gatherv(isli_sub,ijn(mm1),mpi_itype,&
isli_all,ijn,displs_g,mpi_itype,mype_so ,&
mpi_comm_world,ierror)
! Only MPI task mype_so writes the surface file.
if (mype==mype_so ) then
write(*,'(a,5(1x,a6))') 'write_gfsnc_sfc_nst:',fname_sfcges,fname_nstges,fname_sfctsk,fname_sfcanl,fname_nstanl
!
! get Tf analysis increment and surface mask at analysis (lower resolution) grids
!
do i=1,iglobal
ilon=ltosj(i)
ilat=ltosi(i)
dsfct_glb(ilat,ilon) = dsfct_all(i)
isli_glb (ilat,ilon) = isli_all (i)
end do
!
! write dsfct_anl to a data file for later use (at eupd step at present)
!
open(io_dtfanl,file=fname_dtfanl,form='unformatted')
write(io_dtfanl) nlon,nlat
write(io_dtfanl) dsfct_glb
write(io_dtfanl) isli_glb
! open nsst guess file
nstges = open_dataset(fname_nstges)
! open surface guess file
sfcges = open_dataset(fname_sfcges)
! open surface gcycle file
sfcgcy = open_dataset(fname_sfcgcy)
! read a few surface guess file header records
! get dimension sizes
ncdim = get_dim(sfcges, 'grid_xt'); lonb = ncdim%len
ncdim = get_dim(sfcges, 'grid_yt'); latb = ncdim%len
! read some nsst guess file header records (dimensions)
! get dimension sizes
ncdim = get_dim(nstges, 'grid_xt'); lonb_nst = ncdim%len
ncdim = get_dim(nstges, 'grid_yt'); latb_nst = ncdim%len
! check the dimensions consistency in sfc, nst files and the used.
if ( latb /= latb_nst .or. lonb /= lonb_nst ) then
write(6,*) 'Inconsistent dimension for sfc & nst files. latb,lonb : ',latb,lonb, &
'latb_nst,lonb_nst : ',latb_nst,lonb_nst
call stop2(80)
endif
if ( nlat_sfc /= latb+2 .or. nlon_sfc /= lonb ) then
write(6,*) 'Inconsistent dimension for used and read. nlat_sfc,nlon_sfc : ',nlat_sfc,nlon_sfc, &
'latb+2,lonb :',latb+2,lonb
call stop2(81)
endif
!
allocate(slmsk_ges(lonb,latb),slmsk_anl(lonb,latb))
! read slmsk in fname_sfcges to get slmsk_ges
call read_vardata(sfcges, 'land', slmsk_ges)
! read slmsk in fname_sfcgcy to get slmsk_anl
call read_vardata(sfcgcy, 'land', slmsk_anl)
!
! Replace header record date with analysis time from iadate
!
jdate(1) = iadate(1) ! analysis year
jdate(2) = iadate(2) ! analysis month
jdate(3) = iadate(3) ! analysis day
jdate(4) = iadate(4) ! analysis hour
jdate(5) = iadate(5) ! analysis minute
jdate(5) = 0 ! analysis minute
jdate(6) = 0 ! analysis scaled seconds
fhour = zero
odate(1) = jdate(4) !hour
odate(2) = jdate(2) !month
odate(3) = jdate(3) !day
odate(4) = jdate(1) !year
if ( (latb /= nlatm2) .or. (lonb /= nlon) ) then
write(6,*)'WRITE_GFSNC_SFC_NST: different grid dimensions analysis vs sfc. interpolating sfc temperature ',&
', nlon,nlat-2=',nlon,nlatm2,' -vs- sfc file lonb,latb=',lonb,latb
write(6,*) ' WRITE_GFSNC_SFC_NST, nlon_sfc,nlat_sfc : ', nlon_sfc,nlat_sfc
!
! Get the expanded values for a surface type (0 = water now) and the new mask
!
call int2_msk_glb_prep(dsfct_glb,isli_glb,dsfct_tmp,isli_tmp,nlat,nlon,0,nprep)
!
! Get updated/analysis surface mask info from sfcgcy file
!
call tran_gfsncsfc(slmsk_anl,work,lonb,latb)
do j=1,lonb
do i=1,latb+2
isli_gsi(i,j) = nint(work(i,j))
end do
end do
!
! Interpolate dsfct_tmp(nlat,nlon) to dsfct_gsi(nlat_sfc,nlon_sfc) with surface mask accounted
!
call int22_msk_glb(dsfct_tmp,isli_tmp,rlats,rlons,nlat,nlon, &
dsfct_gsi,isli_gsi,rlats_sfc,rlons_sfc,nlat_sfc,nlon_sfc,0)
!
! transform the dsfct_gsi(latb+2,lonb) to dsfct_anl(lonb,latb) for sfc file format
!
do j = 1, latb
do i = 1, lonb
dsfct_anl(i,j) = dsfct_gsi(latb+2-j,i)
end do
end do
else
!
! transform the dsfct_glb(nlat,nlon) to dsfct_anl(lonb,latb) for sfc file
! format when nlat == latb-2 & nlon = lonb
!
do j=1,latb
do i=1,lonb
dsfct_anl(i,j)=dsfct_glb(latb+1-j,i)
end do
end do
endif ! if ( (latb /= nlatm2) .or. (lonb /= nlon) ) then
!
! Start to write output sfc file : fname_sfcanl & fname_nstanl
! open new output file with new header gfile_sfcanl and gfile_nstanl with "write" access.
! Use this call to update header as well
!
! copy input header info to output header info for sfcanl, need to do this before nemsio_close(gfile)
!
sfcanl = create_dataset(fname_sfcanl, sfcgcy, copy_vardata=.true.)
call write_vardata(sfcanl, 'time', fhour)
time_units = get_time_units_from_idate(jdate)
call write_attribute(sfcanl, 'units', time_units, 'time')
sfctsk = create_dataset(fname_sfctsk, sfcgcy, copy_vardata=.true.)
call write_vardata(sfctsk, 'time', fhour)
time_units = get_time_units_from_idate(jdate)
call write_attribute(sfctsk, 'units', time_units, 'time')
!
! copy input header info to output header info for nstanl, need to do this before nemsio_close(gfile)
!
nstanl = create_dataset(fname_nstanl, nstges, copy_vardata=.true.)
call write_vardata(nstanl, 'time', fhour)
time_units = get_time_units_from_idate(jdate)
call write_attribute(nstanl, 'units', time_units, 'time')
!
! For sfcanl, Only tsea (sea surface temperature) will be updated in the SFC
! Need values from nstges for tref update
! read tsea from sfcges
call read_vardata(sfcges, 'tmpsfc', tsea)
! For nstanl, Only tref (foundation temperature) is updated by analysis
! others are updated for snow melting case
! read 18 nsst variables from nstges
! xt
call read_vardata(nstges, 'xt', xt)
! xs
call read_vardata(nstges, 'xs', xs)
! xu
call read_vardata(nstges, 'xu', xu)
! xv
call read_vardata(nstges, 'xv', xv)
! xz
call read_vardata(nstges, 'xz', xz)
! zm
call read_vardata(nstges, 'zm', zm)
! xtts
call read_vardata(nstges, 'xtts', xtts)
! xzts
call read_vardata(nstges, 'xzts', xzts)
! dt_cool
call read_vardata(nstges, 'dtcool', dt_cool)
! z_c
call read_vardata(nstges, 'zc', z_c)
! c_0
call read_vardata(nstges, 'c0', c_0)
! c_d
call read_vardata(nstges, 'cd', c_d)
! w_0
call read_vardata(nstges, 'w0', w_0)
! w_d
call read_vardata(nstges, 'wd', w_d)
! tref
call read_vardata(nstges, 'tref', tref)
! d_conv
call read_vardata(nstges, 'dconv', d_conv)
! ifd
! CRM - does this exist? what is it's name??
!call read_vardata(nstges, 'ifd', ifd)
! qrain
call read_vardata(nstges, 'qrain', qrain)
!
! update tref (in nst file) & tsea (in the surface file) when Tr analysis is on
! reset NSSTM variables for new open water grids
!
if ( nst_gsi > 2 ) then
!
! For the new open water (sea ice just melted) grids, (1) set dsfct_anl = zero; (2) reset the NSSTM variables
!
! Notes: slmsk_ges is the mask of the background
! slmsk_anl is the mask of the analysis
!
where ( slmsk_anl(:,:) == zero .and. slmsk_ges(:,:) == two )
dsfct_anl(:,:) = zero
xt(:,:) = zero
xs(:,:) = zero
xu(:,:) = zero
xv(:,:) = zero
xz(:,:) = z_w_max
zm(:,:) = zero
xtts(:,:) = zero
xzts(:,:) = zero
dt_cool(:,:) = zero
z_c(:,:) = zero
c_0(:,:) = zero
c_d(:,:) = zero
w_0(:,:) = zero
w_d(:,:) = zero
d_conv(:,:) = zero
ifd(:,:) = zero
tref(:,:) = tfrozen
qrain(:,:) = zero
end where
!
! update analysis variable: Tref (foundation temperature) for nst file
!
where ( slmsk_anl(:,:) == zero )
tref(:,:) = max(tref(:,:) + dsfct_anl(:,:),tfrozen)
elsewhere
tref(:,:) = tsea(:,:)
end where
!
! update SST: tsea for sfc file with NSST profile
!
r_zsea1 = 0.001_r_single*real(zsea1)
r_zsea2 = 0.001_r_single*real(zsea2)
call dtzm_2d(xt,xz,dt_cool,z_c,slmsk_anl,r_zsea1,r_zsea2,lonb,latb,dtzm)
where ( slmsk_anl(:,:) == zero )
tsea(:,:) = max(tref(:,:) + dtzm(:,:), tfrozen)
end where
else ! when (nst_gsi <= 2)
do j=1,latb
do i=1,lonb
tref(i,j) = tsea(i,j) ! keep tref as tsea before analysis
end do
end do
!
! For the new open water (sea ice just melted) grids, reset the NSSTM variables
!
where ( slmsk_anl(:,:) == zero .and. slmsk_ges(:,:) == two )
xt(:,:) = zero
xs(:,:) = zero
xu(:,:) = zero
xv(:,:) = zero
xz(:,:) = z_w_max
zm(:,:) = zero
xtts(:,:) = zero
xzts(:,:) = zero
dt_cool(:,:) = zero
z_c(:,:) = zero
c_0(:,:) = zero
c_d(:,:) = zero
w_0(:,:) = zero
w_d(:,:) = zero
d_conv(:,:) = zero
ifd(:,:) = zero
tref(:,:) = tfrozen
qrain(:,:) = zero
end where
!
! update tsea when NO Tf analysis
!
do j=1,latb
do i=1,lonb
tsea(i,j) = max(tsea(i,j) + dsfct_anl(i,j),tfrozen)
end do
end do
endif ! if ( nst_gsi > 2 ) then
!
! update tsea record in sfcanl
!
call write_vardata(sfcanl, 'tmpsfc', tsea)
write(6,100) fname_sfcanl,lonb,latb,houra,iadate(1:4)
100 format(' WRITE_GFSNCIO_SFC_NST: update tsea in ',a6,2i6,1x,f4.1,4(i4,1x))
!
! update tsea record in sfctsk
!
call write_vardata(sfctsk, 'tmpsfc', tsea)
write(6,101) fname_sfctsk,lonb,latb,houra,iadate(1:4)
101 format(' WRITE_GFSNCIO_SFC_NST: update tsea in ',a6,2i6,1x,f4.1,4(i4,1x))
!
! update nsst records in nstanl
!
! slmsk
call write_vardata(nstanl, 'land', slmsk_anl)
! xt
call write_vardata(nstanl, 'xt', xt)
! xs
call write_vardata(nstanl, 'xs', xs)
! xu
call write_vardata(nstanl, 'xu', xu)
! xv
call write_vardata(nstanl, 'xv', xv)
! xz
call write_vardata(nstanl, 'xz', xz)
! zm
call write_vardata(nstanl, 'zm', zm)
! xtts
call write_vardata(nstanl, 'xtts', xtts)
! xzts
call write_vardata(nstanl, 'xzts', xzts)
! z_0
call write_vardata(nstanl, 'dtcool', dt_cool)
! z_c
call write_vardata(nstanl, 'zc', z_c)
! c_0
call write_vardata(nstanl, 'c0', c_0)
! c_d
call write_vardata(nstanl, 'cd', c_d)
! w_0
call write_vardata(nstanl, 'w0', w_0)
! w_d
call write_vardata(nstanl, 'wd', w_d)
! d_conv
call write_vardata(nstanl, 'dconv', d_conv)
! ifd
! CRM See above ifd issue/comment
!call write_vardata(nstanl, 'ifd', ifd)
! tref
call write_vardata(nstanl, 'tref', tref)
! qrain
call write_vardata(nstanl, 'qrain', qrain)
write(6,200) fname_nstanl,lonb,latb,houra,iadate(1:4)
200 format(' WRITE_GFSNCIO_SFC_NST: update variables in ',a6,2i6,1x,f4.1,4(i4,1x))
deallocate(xt,xs,xu,xv,xz,zm,xtts,xzts,dt_cool,z_c,c_0,c_d,w_0,w_d,d_conv,ifd,tref,qrain)
call close_dataset(sfcges)
call close_dataset(sfcgcy)
call close_dataset(nstges)
call close_dataset(sfcanl)
call close_dataset(nstanl)
call close_dataset(sfctsk)
write(6,'(a,'': gfsncio sfc_nst anal written for lonb,latb= '',2i6,'',valid hour= '',f4.1,'',idate= '',4i5)') &
trim(my_name),lonb,latb,fhour,odate
endif
end subroutine write_sfc_nst_
subroutine intrp22(a,rlons_a,rlats_a,nlon_a,nlat_a, &
b,rlons_b,rlats_b,nlon_b,nlat_b)
!$$$ subprogram documentation block
! . . .
! subprogram: intrp22 --- interpolates from one 2-d grid to another 2-d grid
! like analysis to surface grid or vice versa
! prgrmmr: li - initial version; org: np2
!
! abstract: This routine interpolates a grid to b grid
!
! program history log:
!
! input argument list:
! rlons_a - longitudes of input array
! rlats_a - latitudes of input array
! nlon_a - number of longitude of input array
! nlat_a - number of latitude of input array
! rlons_b - longitudes of output array
! rlats_b - latitudes of output array
! nlon_b - number of longitude of output array
! nlat_b - number of latitude of output array
! a - input values
!
! output argument list:
! b - output 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: zero,one
implicit none
! !INPUT PARAMETERS:
integer(i_kind) ,intent(in ) :: nlon_a,nlat_a,nlon_b,nlat_b
real(r_kind), dimension(nlon_a) ,intent(in ) :: rlons_a
real(r_kind), dimension(nlat_a) ,intent(in ) :: rlats_a
real(r_kind), dimension(nlon_b) ,intent(in ) :: rlons_b
real(r_kind), dimension(nlat_b) ,intent(in ) :: rlats_b
real(r_single), dimension(nlon_a,nlat_a),intent(in ) :: a
! !OUTPUT PARAMETERS:
real(r_single), dimension(nlon_b,nlat_b),intent( out) :: b
! 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,nlat_b
dlat=rlats_b(j)
call grdcrd1(dlat,rlats_a,nlat_a,1)
iy=int(dlat)
iy=min(max(1,iy),nlat_a)
dy =dlat-iy
dy1 =one-dy
iyp=min(nlat_a,iy+1)
do i=1,nlon_b
dlon=rlons_b(i)
call grdcrd1(dlon,rlons_a,nlon_a,1)
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(0,ix),nlon_a)
ixp=ix+1
if(ix==0) ix=nlon_a
if(ixp==nlon_a+1) ixp=1
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 intrp22
subroutine tran_gfsncsfc(ain,aout,lonb,latb)
!$$$ subprogram documentation block
! . . . .
! subprogram: tran_gfsncsfc transform gfs surface file to analysis grid
! prgmmr: derber org: np2 date: 2003-04-10
!
! abstract: transform gfs surface file to analysis grid
!
! program history log:
! 2012-31-38 derber - initial routine
!
! input argument list:
! ain - input surface record on processor iope
! lonb - input number of longitudes
! latb - input number of latitudes
!
! output argument list:
! aout - output transposed surface record
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind,r_single
use constants, only: zero
use sfcio_module, only: sfcio_realkind
implicit none
! Declare passed variables
integer(i_kind) ,intent(in ) :: lonb,latb
real(sfcio_realkind),dimension(lonb,latb),intent(in ) :: ain
real(r_single),dimension(latb+2,lonb),intent(out) :: aout
! Declare local variables
integer(i_kind) i,j
real(r_kind) sumn,sums
! of surface guess array
sumn = zero
sums = zero
do i=1,lonb
sumn = ain(i,1) + sumn
sums = ain(i,latb) + sums
end do
sumn = sumn/float(lonb)
sums = sums/float(lonb)
! Transfer from local work array to surface guess array
do j = 1,lonb
aout(1,j)=sums
do i=2,latb+1
aout(i,j) = ain(j,latb+2-i)
end do
aout(latb+2,j)=sumn
end do
return
end subroutine tran_gfsncsfc
end module netcdfgfs_io