From 9456a0bfe5b1200f7ce17b15134348106e9060a6 Mon Sep 17 00:00:00 2001 From: "russ.treadon" Date: Mon, 25 Nov 2024 19:55:19 +0000 Subject: [PATCH 1/9] modify GSI to read netcdf format global_berror files (#808) --- src/gsi/gsi_files.cmake | 2 +- src/gsi/m_berror_stats.f90 | 504 +++++++++++++++----- src/gsi/m_nc_berror.f90 | 908 +++++++++++++++++++++++++++++++++++++ src/gsi/ncepnems_io.f90 | 175 ++++--- 4 files changed, 1416 insertions(+), 173 deletions(-) create mode 100644 src/gsi/m_nc_berror.f90 diff --git a/src/gsi/gsi_files.cmake b/src/gsi/gsi_files.cmake index 04e7926300..d4b6e6b3ed 100644 --- a/src/gsi/gsi_files.cmake +++ b/src/gsi/gsi_files.cmake @@ -357,6 +357,7 @@ m_lightNode.F90 m_lwcpNode.F90 m_mitmNode.F90 m_mxtmNode.F90 +m_nc_berror.f90 m_o3lNode.F90 m_obsLList.F90 m_obsNode.F90 @@ -497,7 +498,6 @@ read_goesndr.f90 read_gps.f90 read_guess.F90 read_iasi.f90 -read_iasing.f90 read_l2bufr_mod.f90 read_lag.f90 read_lidar.f90 diff --git a/src/gsi/m_berror_stats.f90 b/src/gsi/m_berror_stats.f90 index 088a7619fe..8ca7b44aea 100644 --- a/src/gsi/m_berror_stats.f90 +++ b/src/gsi/m_berror_stats.f90 @@ -11,6 +11,7 @@ module m_berror_stats ! program history log: ! 2010-03-24 j guo - added this document block ! 2011-08-01 lueken - changed F90 to f90 (no machine logic) and fix indentation +! 2014-04-01 weir - added some chem support ! ! input argument list: see Fortran 90 style document below ! @@ -34,22 +35,27 @@ module m_berror_stats ! ! !INTERFACE: - use kinds, only : i_kind + use kinds, only: i_kind,r_kind use constants, only: one,zero use control_vectors,only: cvars2d,cvars3d - use mpeu_util, only: getindex,check_iostat + use mpeu_util, only: getindex,check_iostat,die + + use m_nc_berror, only: nc_berror_dims + use m_nc_berror, only: nc_berror_read + use m_nc_berror, only: nc_berror_vars + use m_nc_berror, only: nc_berror_vars_final + use m_nc_berror, only: nc_berror_getpointer + use module_ncio, only: open_dataset, Dataset, Dimension, & + close_dataset, read_vardata, get_dim implicit none private ! except -! def usenewgfsberror - use modified gfs berror stats for global and regional. -! for global skips extra record -! for regional properly defines array sizes, etc. -! def berror_stats - reconfigurable filename via NAMELIST/setup/ - + ! reconfigurable parameters, via NAMELIST/setup/ public :: usenewgfsberror - public :: berror_stats,inquire_berror + public :: berror_stats,inquire_berror ! reconfigurable filename + ! interfaces to file berror_stats. public :: berror_get_dims ! get dimensions, jfunc::createj_func() public :: berror_read_bal ! get cross-cov.stats., balmod::prebal() @@ -74,6 +80,7 @@ module m_berror_stats ! 18Dec15 - Rahul.Mahajan ! - replace die calls with check_iostat to clean code ! - fix haphazard indentation and add return to all sub-routines +! 30Aug21 - Todling - introduce netcdf capability !EOP ___________________________________________________________________ character(len=*),parameter :: myname='m_berror_stats' @@ -81,10 +88,14 @@ module m_berror_stats integer(i_kind),parameter :: default_unit_ = 22 integer(i_kind),parameter :: default_rc_ = 2 + character(len=256),save :: berror_stats = "berror_stats" ! filename logical,save :: cwcoveqqcov_ - logical usenewgfsberror + logical usenewgfsberror - character(len=256),save :: berror_stats = "berror_stats" ! filename + logical,save :: bin_berror=.false. + +!! real(r_kind),allocatable,dimension(:,:):: varq +!! real(r_kind),allocatable,dimension(:,:):: varcw contains @@ -98,10 +109,13 @@ module m_berror_stats ! ! !INTERFACE: +!!subroutine get_dims(mype,msig,mlat,mlon,lunit) subroutine get_dims(msig,mlat,mlon,lunit) + use mpimod, only: mype implicit none +!! integer(i_kind) ,intent( in) :: mype ! proc identifier integer(i_kind) ,intent( out) :: msig ! dimension of levels integer(i_kind) ,intent( out) :: mlat ! dimension of latitudes integer(i_kind),optional,intent( out) :: mlon ! dimension of longitudes @@ -114,7 +128,31 @@ subroutine get_dims(msig,mlat,mlon,lunit) !EOP ___________________________________________________________________ character(len=*),parameter :: myname_=myname//'::get_dims' - integer(i_kind) :: inerr,mlon_,ier + integer(i_kind) :: inerr,mlon_,status + +! Try reading as NetCDF ... + if(mype==0) print*, myname_, ": Try reading berror from NetCDF file" + call nc_(status) + if (status==0) then + if(mype==0) print*, myname_, ": Reading berror from NetCDF file" + bin_berror = .false. + else ! if failed, read as NetCDF + call bin_(status) + if (status==0) then + bin_berror = .true. + if(mype==0) print*, myname_, ": Read berror from Binary file" + else + if(mype==0) call die(myname_,'Failed reading Berror file', 99) + endif + endif + if ( present(mlon) ) mlon = mlon_ + + return + + contains + + subroutine bin_(ier) + integer,intent(inout) :: ier ! Read dimension of stats file inerr = default_unit_ @@ -123,12 +161,26 @@ subroutine get_dims(msig,mlat,mlon,lunit) call check_iostat(ier,myname_,'open('//trim(berror_stats)//')') rewind inerr read(inerr,iostat=ier) msig,mlat,mlon_ + if (ier/=0) return call check_iostat(ier,myname_,'read header') close(inerr,iostat=ier) call check_iostat(ier,myname_,'close('//trim(berror_stats)//')') - if ( present(mlon) ) mlon = mlon_ return + end subroutine bin_ + + subroutine nc_(ier) + integer,intent(inout) :: ier + + integer :: mlat_, mlev_ + + call nc_berror_dims (berror_stats,mlat_,mlon_,mlev_,ier, mype,0) + + mlat = mlat_ + msig = mlev_ + + return + end subroutine nc_ end subroutine get_dims subroutine inquire_berror(lunit,mype) @@ -140,30 +192,55 @@ subroutine inquire_berror(lunit,mype) integer(i_kind),intent(in ) :: lunit ! logical unit [22] integer(i_kind),intent(in ) :: mype + logical :: ncio character(len=*),parameter :: myname_=myname//'::inquire_berror' integer(i_kind) :: inerr,msig,mlat,mlon_,ier,errtot,i real(r_single),dimension(:),allocatable:: clat_avn,sigma_avn + type(Dataset) :: dset + type(Dimension) :: levdim + + ! Determine type of berror_stats file: binary or netcdf + dset = open_dataset(berror_stats,errcode=inerr) + if (inerr==0) then + ! this is a netcdf file + ncio = .true. + else + ! this is a binary file + ncio = .false. + endif + ! Read dimension of stats file - inerr = lunit - open(inerr,file=berror_stats,form='unformatted',status='old',iostat=ier) - call check_iostat(ier,myname_,'open('//trim(berror_stats)//')') - rewind inerr - read(inerr,iostat=ier) msig,mlat,mlon_ - call check_iostat(ier,myname_,'read header') - errtot=ier + if (ncio) then + levdim = get_dim(dset,'lev'); msig = levdim%len + allocate(sigma_avn(1:msig)) + call read_vardata(dset,'lev',sigma_avn) + call close_dataset(dset) + else + inerr = lunit + open(inerr,file=berror_stats,form='unformatted',status='old',iostat=ier) + call check_iostat(ier,myname_,'open('//trim(berror_stats)//')') + rewind inerr + read(inerr,iostat=ier) msig,mlat,mlon_ + call check_iostat(ier,myname_,'read header') + errtot=ier + + allocate ( clat_avn(mlat) ) + allocate ( sigma_avn(1:msig) ) + read(inerr,iostat=ier)clat_avn,sigma_avn + close(inerr,iostat=ier) + call check_iostat(ier,myname_,'close('//trim(berror_stats)//')') + deallocate(clat_avn) + endif - errtot=0 - allocate ( clat_avn(mlat) ) - allocate ( sigma_avn(1:msig) ) - read(inerr,iostat=ier)clat_avn,sigma_avn ! Checking to see if sigma_avn fits required format if so newgfsberror file. + errtot=0 do i=1,msig-1 if(sigma_avn(i) < sigma_avn(i+1) .or. sigma_avn(i) < zero .or. sigma_avn(i) > one)then errtot=1 end if end do - deallocate(clat_avn,sigma_avn) + deallocate(sigma_avn) if(errtot /= 0)then usenewgfsberror=.false. if(mype == 0)write(6,*) 'usenewgfsberror = .false. old format file ' @@ -171,8 +248,6 @@ subroutine inquire_berror(lunit,mype) usenewgfsberror=.true. if(mype == 0)write(6,*) 'usenewgfsberror = .true. new format file ' end if - close(inerr,iostat=ier) - call check_iostat(ier,myname_,'close('//trim(berror_stats)//')') return end subroutine inquire_berror @@ -254,6 +329,19 @@ subroutine read_bal(agvin,bvin,wgvin,pputin,fut2ps,mype,lunit) integer(i_kind) :: nsigstat,nlatstat integer(i_kind) :: inerr,ier + ier=0 + if (bin_berror) then + call bin_ + else + call nc_(mype) + endif + + return + + contains + + subroutine bin_ + ! Open background error statistics file inerr = default_unit_ if ( present(lunit) ) inerr = lunit @@ -265,10 +353,9 @@ subroutine read_bal(agvin,bvin,wgvin,pputin,fut2ps,mype,lunit) rewind inerr read(inerr,iostat=ier) nsigstat,nlatstat + if (ier/=0) return call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (nsigstat,nlatstat)') -! dummy read to skip lats,sigma - if(usenewgfsberror)read(inerr) if ( mype==0 ) then if ( nsig/=nsigstat .or. nlat/=nlatstat ) then write(6,*) myname_,'(PREBAL): ***ERROR*** resolution of ', & @@ -298,7 +385,25 @@ subroutine read_bal(agvin,bvin,wgvin,pputin,fut2ps,mype,lunit) close(inerr,iostat=ier) call check_iostat(ier,myname_,'close('//trim(berror_stats)//')') - return + end subroutine bin_ + + subroutine nc_(myid) + integer, intent(in) :: myid + type(nc_berror_vars) bvars + if ( fut2ps ) then + call die(myname_," fut2ps not available in this form "//trim(berror_stats), 99) + endif + call nc_berror_read (berror_stats,bvars,ier, myid=myid,root=0) + if (nlat/=bvars%nlat .or. nsig/=bvars%nsig ) then + call die(myname_," inconsistent dims in "//trim(berror_stats), 99) + endif + agvin = bvars%tcon + bvin = bvars%vpcon + wgvin = bvars%pscon + pputin=zero + call nc_berror_vars_final(bvars) + end subroutine nc_ + end subroutine read_bal !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -312,8 +417,9 @@ end subroutine read_bal ! !INTERFACE: subroutine read_wgt(corz,corp,hwll,hwllp,vz,corsst,hsst,varq,qoption,varcw,cwoption,mype,lunit) +!! n_clouds_fwd,cloud_names_fwd,lunit) - use kinds,only : r_single,r_kind + use kinds, only : r_single,r_kind use gridmod,only : nlat,nlon,nsig use radiance_mod, only: n_clouds_fwd,cloud_names_fwd @@ -331,11 +437,15 @@ subroutine read_wgt(corz,corp,hwll,hwllp,vz,corsst,hsst,varq,qoption,varcw,cwopt real(r_kind), dimension(:,:) ,intent(out ) :: varq real(r_kind), dimension(:,:) ,intent(out ) :: varcw - + integer(i_kind) ,intent(in ) :: qoption integer(i_kind) ,intent(in ) :: cwoption integer(i_kind) ,intent(in ) :: mype ! "my" processor ID + +! Optionals integer(i_kind),optional ,intent(in ) :: lunit ! an alternative unit +!! integer(i_kind), optional ,intent(in ) :: n_clouds_fwd +!! character(len=*),optional ,intent(in ) :: cloud_names_fwd(:) ! !REVISION HISTORY: ! 30Jul08 - Jing Guo @@ -363,19 +473,104 @@ subroutine read_wgt(corz,corp,hwll,hwllp,vz,corsst,hsst,varq,qoption,varcw,cwopt integer(i_kind) :: i,n,k,iq,icw,ivar,ic integer(i_kind) :: inerr,istat,ier - integer(i_kind) :: nsigstat,nlatstat,mlon_ + integer(i_kind) :: msig,mlat + integer(i_kind) :: nsigstat,nlatstat integer(i_kind) :: isig real(r_kind) :: corq2x character(len=5) :: var logical,allocatable,dimension(:) :: found3d logical,allocatable,dimension(:) :: found2d -! real(r_single),allocatable,dimension(:) :: clat,sigma real(r_single),allocatable,dimension(:,:):: hwllin real(r_single),allocatable,dimension(:,:):: corzin real(r_single),allocatable,dimension(:,:):: corq2 real(r_single),allocatable,dimension(:,:):: vscalesin + allocate(found3d(size(cvars3d)),found2d(size(cvars2d))) + found3d=.false. + found2d=.false. + + call berror_get_dims(msig,mlat) + if ( bin_berror ) then + write(6,*)'call binary berror read' + call bin_() + else + write(6,*)'call netcdf berror read' + call nc_(mype) + endif + + ! corz, hwll & vz for undefined 3d variables + do n=1,size(cvars3d) + if ( .not.found3d(n) ) then + if ( n>0 ) then + if ( trim(cvars3d(n))=='oz' ) then + call setcoroz_(corz(:,:,n),mype) + call sethwlloz_(hwll(:,:,n),mype) + else + call setcorchem_(cvars3d(n),corz(:,:,n),ier) + call sethwllchem_(hwll(:,:,n),mype) + if(ier/=0) cycle ! if this happens, code will crash later + endif + call setvscalesoz_(vz(:,:,n)) + endif + if ( mype==0 ) write(6,*) myname_, ': WARNING, using general Berror template for ', cvars3d(n) + endif + enddo + +! if so, overwrite cw-cov with q-cov + iq=-1;icw=-1 + do n=1,size(cvars3d) + if(trim(cvars3d(n))=='q' ) iq =n + if(trim(cvars3d(n))=='cw') icw=n + enddo + if (cwcoveqqcov_) then + if(iq>0.and.icw>0) then + hwll(:,:,icw)=hwll(:,:,iq) + vz (:,:,icw)=vz (:,:,iq) + end if + end if + if (cwoption==1 .or. cwoption==3) then + if (iq>0.and.icw>0) then + do k=1,nsig + do i=1,nlat + corz(i,k,icw)=one + end do + end do + hwll(:,:,icw)=0.5_r_kind*hwll(:,:,iq) + vz (:,:,icw)=0.5_r_kind*vz (:,:,iq) + end if + +!! if (present(n_clouds_fwd) .and. present(cloud_names_fwd)) then + if (n_clouds_fwd>0 .and. icw<=0) then + do n=1,size(cvars3d) + do ic=1,n_clouds_fwd + if(trim(cvars3d(n))==trim(cloud_names_fwd(ic))) then + ivar=n + do k=1,nsig + do i=1,nlat + corz(i,k,ivar)=one + end do + end do + hwll(:,:,ivar)=0.5_r_kind*hwll(:,:,iq) + vz (:,:,ivar)=0.5_r_kind*vz (:,:,iq) + exit + end if + end do + end do + end if +!! end if + endif + + ! need simliar general template for undefined 2d variables ... + + deallocate(found3d,found2d) + + return + + contains + + subroutine bin_ + ! Open background error statistics file inerr = default_unit_ if ( present(lunit) ) inerr = lunit @@ -386,9 +581,8 @@ subroutine read_wgt(corz,corp,hwll,hwllp,vz,corsst,hsst,varq,qoption,varcw,cwopt ! with that specified via the user namelist rewind inerr - read(inerr,iostat=ier)nsigstat,nlatstat,mlon_ + read(inerr,iostat=ier)nsigstat,nlatstat call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (nsigstat,nlatstat)') - if(usenewgfsberror)read(inerr,iostat=ier) if ( mype==0 ) then if ( nsigstat/=nsig .or. nlatstat/=nlat ) then @@ -407,9 +601,6 @@ subroutine read_wgt(corz,corp,hwll,hwllp,vz,corsst,hsst,varq,qoption,varcw,cwopt call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (agvin,bvin,wgvin)') ! Read amplitudes - allocate(found3d(size(cvars3d)),found2d(size(cvars2d))) - found3d=.false. - found2d=.false. readloop: do read(inerr,iostat=istat) var, isig if ( istat/=0 ) exit @@ -455,6 +646,7 @@ subroutine read_wgt(corz,corp,hwll,hwllp,vz,corsst,hsst,varq,qoption,varcw,cwopt do i=1,nlat corq2x=corq2(i,k) varq(i,k)=min(max(corq2x,0.00015_r_kind),one) +!! varq_out(i,k) = varq(i,k) enddo enddo do k=1,isig @@ -468,6 +660,7 @@ subroutine read_wgt(corz,corp,hwll,hwllp,vz,corsst,hsst,varq,qoption,varcw,cwopt do i=1,nlat corq2x=corq2(i,k) varcw(i,k)=max(corq2x,zero) +!! varcw_out(i,k) = varcw(i,k) enddo enddo do k=1,isig @@ -498,73 +691,107 @@ subroutine read_wgt(corz,corp,hwll,hwllp,vz,corsst,hsst,varq,qoption,varcw,cwopt deallocate(corzin,hwllin) if ( isig>1 ) deallocate(vscalesin) if ( var=='q' .or. var=='cw' ) deallocate(corq2) - enddo readloop + enddo readloop close(inerr) - ! corz, hwll & vz for undefined 3d variables - do n=1,size(cvars3d) - if ( .not.found3d(n) ) then - if ( n>0 ) then - if ( cvars3d(n)=='oz' ) then - call setcoroz_(corz(:,:,n),mype) - else - call setcorchem_(cvars3d(n),corz(:,:,n),ier) - if ( ier/=0 ) cycle ! if this happens, code will crash later - endif - call sethwlloz_(hwll(:,:,n),mype) - call setvscalesoz_(vz(:,:,n)) - endif - if ( mype==0 ) write(6,*) myname_, ': WARNING, using general Berror template for ', cvars3d(n) - endif - enddo - -! if so, overwrite cw-cov with q-cov - iq=-1;icw=-1 - do n=1,size(cvars3d) - if(trim(cvars3d(n))=='q' ) iq =n - if(trim(cvars3d(n))=='cw') icw=n - enddo - if (cwcoveqqcov_) then - if(iq>0.and.icw>0) then - hwll(:,:,icw)=hwll(:,:,iq) - vz (:,:,icw)=vz (:,:,iq) - end if - end if - if (cwoption==1 .or. cwoption==3) then - if (iq>0.and.icw>0) then - do k=1,nsig - do i=1,nlat - corz(i,k,icw)=one - end do - end do - hwll(:,:,icw)=0.5_r_kind*hwll(:,:,iq) - vz (:,:,icw)=0.5_r_kind*vz (:,:,iq) - end if - if (n_clouds_fwd>0 .and. icw<=0) then - do n=1,size(cvars3d) - do ic=1,n_clouds_fwd - if(trim(cvars3d(n))==trim(cloud_names_fwd(ic))) then - ivar=n - do k=1,nsig - do i=1,nlat - corz(i,k,ivar)=one - end do - end do - hwll(:,:,ivar)=0.5_r_kind*hwll(:,:,iq) - vz (:,:,ivar)=0.5_r_kind*vz (:,:,iq) - exit - end if - end do - end do - end if - endif + end subroutine bin_ - ! need simliar general template for undefined 2d variables ... + subroutine nc_(myid) - deallocate(found3d,found2d) + integer,intent(in) :: myid + type(nc_berror_vars) bvars + real(r_single), pointer :: ptr1d(:) + real(r_single), pointer :: ptr2d(:,:) + integer :: nv + call nc_berror_read (berror_stats,bvars,ier, myid=myid,root=0) + if (nlat/=bvars%nlat .or. nlon/=bvars%nlon .or. nsig/=bvars%nsig ) then + call die(myname_," inconsistent dims in "//trim(berror_stats), 99) + endif + isig=bvars%nsig + +! RTodling: the following is bad since it wires all naming conventions ... to be revised + do nv=1,size(cvars2d) + if (trim(cvars2d(nv))=='sst') then + n = getindex(cvars2d,'sst') + found2d(n)=.true. + call nc_berror_getpointer (cvars2d(nv),bvars,ptr2d,ier) + if(ier==0) corsst=ptr2d + call nc_berror_getpointer ('h'//cvars2d(nv),bvars,ptr2d,ier) + if(ier==0) hsst=ptr2d + endif + if (trim(cvars2d(nv))=='ps') then + n = getindex(cvars2d,'ps') + found2d(n)=.true. + call nc_berror_getpointer (cvars2d(nv),bvars,ptr1d,ier) + if(ier==0) corp(:,n)=ptr1d + call nc_berror_getpointer ('h'//cvars2d(nv),bvars,ptr1d,ier) + if(ier==0) hwllp(:,n)=ptr1d + endif + enddo + do nv=1,size(cvars3d) + call nc_berror_getpointer (cvars3d(nv),bvars,ptr2d,ier) + if (ier==0) then + n = getindex(cvars3d,cvars3d(nv)) + found3d(n)=.true. + corz(:,:,n)=ptr2d + call nc_berror_getpointer ('h'//trim(cvars3d(nv)),bvars,ptr2d,ier) + if(ier==0) hwll(:,:,n)=ptr2d + call nc_berror_getpointer ('v'//trim(cvars3d(nv)),bvars,ptr2d,ier) + if(ier==0) vz(:,:,n)=transpose(ptr2d) + if (trim(cvars3d(nv))=='cw' .and. cwoption==2) then + allocate(corq2(bvars%nlat,bvars%nsig)) + call nc_berror_getpointer ('nrh',bvars,ptr2d,ier) + if (ier==0) then + corq2=ptr2d + do k=1,bvars%nsig + do i=1,bvars%nlat + corq2x=corq2(i,k) + varcw(i,k)=max(corq2x,zero) + enddo + enddo + else + call die(myname_," in cw, failed to find bvar nrh ", 99) + endif + corz(:,:,n)=one + deallocate(corq2) + endif + if (trim(cvars3d(nv))=='q' .and. qoption==2) then + allocate(corq2(bvars%nlat,bvars%nsig)) + call nc_berror_getpointer ('nrh',bvars,ptr2d,ier) + if (ier==0) then + corq2=ptr2d +! corq2=max(0.0_r_kind,corq2) ! hack 1 +! corq2=min(1.0_r_kind,2*corq2) ! hack 2 +! print *, 'DEBUG (berr): ', minval(corq2),maxval(corq2) + do k=1,bvars%nsig + do i=1,bvars%nlat + corq2x=corq2(i,k) + varq(i,k)=min(max(corq2x,0.00015_r_kind),one) + enddo + enddo + else + call die(myname_," in q, failed to find bvar nrh ", 99) + endif + corz(:,:,n)=one + deallocate(corq2) + endif + cycle + endif + if (trim(cvars3d(nv))=='q') then + n = getindex(cvars3d,'q') + found3d(n)=.true. + corz(:,:,n)=bvars%qvar + if (qoption == 2) then + endif + hwll(:,:,n)=bvars%qhln + vz(:,:,n)=transpose(bvars%qvln) + cycle + endif + enddo + call nc_berror_vars_final(bvars) + end subroutine nc_ - return end subroutine read_wgt !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -646,7 +873,7 @@ subroutine setcoroz_(coroz,mype) enddo enddo enddo - work_oz(nsig+1,mm1)=real(lon1*lat1,r_kind) + work_oz(nsig+1,mm1)=float(lon1*lat1) call mpi_allreduce(work_oz,work_oz1,(nsig+1)*npe,mpi_rtype,mpi_sum,& mpi_comm_world,ierror) @@ -694,7 +921,7 @@ end subroutine setcoroz_ subroutine sethwlloz_(hwlloz,mype) - use kinds, only: r_single,r_kind + use kinds, only: r_single,r_kind use mpimod, only: levs_id use gridmod, only: nnnn1o,nsig,nlon,nlat use constants,only: two,three,pi,rearth_equator @@ -725,7 +952,7 @@ subroutine sethwlloz_(hwlloz,mype) do k=1,nnnn1o k1=levs_id(k) if ( k1>0 ) then - write(6,*) myname_,'(PREWGT): mype = ',mype, k1 + if(mype==0) write(6,*) myname_,'(PREWGT): mype = ',mype, k1 if ( k1<=nsig*3/4 ) then ! fact=1./hwl fact=r40000/(r400*nlon) @@ -754,7 +981,7 @@ end subroutine sethwlloz_ subroutine setvscalesoz_(vscalesoz) - use kinds, only: r_single,r_kind + use kinds,only : r_single,r_kind use gridmod,only: nlat,nsig implicit none @@ -788,7 +1015,7 @@ end subroutine setvscalesoz_ subroutine setcorchem_(cname,corchem,rc) - use kinds, only: r_single,r_kind + use kinds, only: r_single,r_kind use mpimod, only: mype use constants,only: zero,one use mpimod, only: npe,mpi_rtype,mpi_sum,mpi_comm_world @@ -809,7 +1036,9 @@ subroutine setcorchem_(cname,corchem,rc) integer(i_kind) ,intent( out) :: rc ! return error code ! !REVISION HISTORY: -! 15Jul20010 - Todling - created from Guo's OZ routine +! 15Jul2010 - Todling - created from Guo's OZ routine +! 20Apr2015 - Weir - relaced chemz with a constant, old approach +! wasn't appropriate for CO ! !EOP ___________________________________________________________________ @@ -869,7 +1098,7 @@ subroutine setcorchem_(cname,corchem,rc) enddo enddo enddo - work_chem(nsig+1,mm1)=real(lon1*lat1,r_kind) + work_chem(nsig+1,mm1)=float(lon1*lat1) call mpi_allreduce(work_chem,work_chem1,(nsig+1)*npe,mpi_rtype,mpi_sum,& mpi_comm_world,ierror) @@ -892,7 +1121,10 @@ subroutine setcorchem_(cname,corchem,rc) do n=1,npe asum=asum+work_chem1(k,n) enddo - if ( bsum>zero ) chemz(k)=asum/bsum +! Not appropriate for co, just replacing with a constant, will revisit +! later (bweir) +! if (bsum>zero) chemz(k)=asum/bsum + chemz(k) = 1._r_kind enddo ! now this part is taken from prewgt(). @@ -907,4 +1139,62 @@ subroutine setcorchem_(cname,corchem,rc) return end subroutine setcorchem_ +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +! NASA/GSFC, Global Modeling and Assimilation Office, 900.3, GEOS/DAS ! +!BOP ------------------------------------------------------------------- +! +! !IROUTINE: sethwllchem_ - a modeled hwll of chem +! +! !DESCRIPTION: +! +! !INTERFACE: + + subroutine sethwllchem_(hwll, mype) + + use kinds, only: r_single, r_kind + use mpimod , only: levs_id + use gridmod, only: nnnn1o, nsig, nlon, nlat + use constants, only: two, three, pi, rearth_equator + + implicit none + + real(r_single), intent( out) :: hwll(nlat,nsig) + integer(i_kind), intent(in ) :: mype + +! !REVISION HISTORY: +! 20May14 - Weir - Initial code, based on sethwlloz_ +!EOP ___________________________________________________________________ + + character(len=*), parameter :: myname_ = myname//'::sethwllchem_' + + real(r_kind), parameter :: r400 = 400._r_kind + real(r_kind), parameter :: r800 = 800._r_kind + real(r_kind), parameter :: r40000 = 40000._r_kind + + integer(i_kind) :: k, k1 + real(r_kind) :: fact + real(r_kind) :: s2u + + if (mype == 0) then + write(6,*) myname_, '(PREWGT): mype = ', mype + end if + + s2u = (two*pi*rearth_equator)/nlon + do k = 1,nnnn1o + k1 = levs_id(k) + if (k1 > 0) then + if (mype == 0) write(6,*) myname_, '(PREWGT): mype = ', mype, k1 +! make everything constant +! fact = real(k1,r_kind)**2._r_kind + fact = 1._r_kind + fact = r40000/(r400*nlon*fact) + hwll(:,k1) = s2u/fact + end if + end do + + if (mype == 0) then + write(6,*) myname_, '(PREWGT): mype = ', mype, 'finish sethwllchem_' + end if + + end subroutine sethwllchem_ end module m_berror_stats diff --git a/src/gsi/m_nc_berror.f90 b/src/gsi/m_nc_berror.f90 new file mode 100644 index 0000000000..559690112d --- /dev/null +++ b/src/gsi/m_nc_berror.f90 @@ -0,0 +1,908 @@ +module m_nc_berror +use netcdf +implicit none +private + +public :: nc_berror_vars_init +public :: nc_berror_vars_final +public :: nc_berror_vars_comp +public :: nc_berror_vars_copy +public :: nc_berror_vars +public :: nc_berror_dims +public :: nc_berror_read +public :: nc_berror_write +public :: nc_berror_getpointer + +type nc_berror_vars + logical :: initialized=.false. + integer :: nlon,nlat,nsig + real(4),pointer,dimension(:,:,:):: tcon + real(4),pointer,dimension(:,:) :: sfvar,vpvar,tvar,qvar,cvar,nrhvar,ozvar + real(4),pointer,dimension(:,:) :: qivar,qlvar,qrvar,qsvar + real(4),pointer,dimension(:,:) :: sfhln,vphln,thln,qhln,chln,ozhln + real(4),pointer,dimension(:,:) :: qihln,qlhln,qrhln,qshln + real(4),pointer,dimension(:,:) :: sfvln,vpvln,tvln,qvln,cvln,ozvln + real(4),pointer,dimension(:,:) :: qivln,qlvln,qrvln,qsvln + real(4),pointer,dimension(:,:) :: vpcon,pscon,varsst,corlsst + real(4),pointer,dimension(:) :: psvar,pshln + real(4),pointer,dimension(:) :: v1d + real(4),pointer,dimension(:,:) :: v2d + real(4),pointer,dimension(:,:,:):: v3d +end type nc_berror_vars + +character(len=*), parameter :: myname = 'm_nc_berror' + +integer, parameter :: nv1d = 2 +character(len=4),parameter :: cvars1d(nv1d) = (/ 'ps ', 'hps ' /) + +integer, parameter :: nv2d = 33 +character(len=5),parameter :: cvars2d(nv2d) = (/ & + 'sf ', 'hsf ', 'vsf ', & + 'vp ', 'hvp ', 'vvp ', & + 't ', 'ht ', 'vt ', & + 'q ', 'hq ', 'vq ', & + 'qi ', 'hqi ', 'vqi ', & + 'ql ', 'hql ', 'vql ', & + 'qr ', 'hqr ', 'vqr ', & + 'qs ', 'hqs ', 'vqs ', & + 'oz ', 'hoz ', 'voz ', & + 'cw ', 'hcw ', 'vcw ', & + 'pscon', 'vpcon', 'nrh ' & + /) + +integer, parameter :: nvmll = 1 ! meriodional, level, level +character(len=4),parameter :: cvarsMLL(nvmll) = (/ 'tcon' /) + +integer, parameter :: nv2dx = 2 +character(len=4),parameter :: cvars2dx(nv2dx) = (/ 'sst ', 'hsst' /) + +interface nc_berror_dims; module procedure & + read_dims_ ; end interface +interface nc_berror_read; module procedure & + read_berror_ ; end interface +interface nc_berror_write; module procedure & + write_berror_ ; end interface +interface nc_berror_vars_init; module procedure & + init_berror_vars_ ; end interface +interface nc_berror_vars_final; module procedure & + final_berror_vars_ ; end interface +interface nc_berror_vars_comp; module procedure & + comp_berror_vars_ ; end interface +interface nc_berror_vars_copy; module procedure & + copy_ ; end interface +interface nc_berror_getpointer + module procedure get_pointer_1d_ + module procedure get_pointer_2d_ +end interface + +contains + +subroutine read_dims_ (fname,nlat,nlon,nlev,rc, myid,root) + implicit none + character(len=*), intent(in) :: fname ! input filename + integer, intent(out) :: rc + integer, intent(out) :: nlat,nlon,nlev + integer, intent(in), optional :: myid, root + +! This will be the netCDF ID for the file and data variable. + integer :: ncid, varid, ier + integer :: mype_,root_ + +! Local variables + character(len=*), parameter :: myname_ = myname//"::dims_" + +! Return code (status) + rc=0; mype_=0; root_=0 + if(present(myid) .and. present(root) ) then + mype_ = myid + root_ = root + endif + +! Open the file. NF90_NOWRITE tells netCDF we want read-only access to +! the file. + call check_( nf90_open(fname, NF90_NOWRITE, ncid), rc, mype_, root_ ) + if(rc/=0) return + +! Read global attributes + call check_( nf90_inq_dimid(ncid, "lon", varid), rc, mype_, root_) + call check_( nf90_inquire_dimension(ncid, varid, len=nlon), rc, mype_, root_ ) + call check_( nf90_inq_dimid(ncid, "lat", varid), rc, mype_, root_ ) + call check_( nf90_inquire_dimension(ncid, varid, len=nlat), rc, mype_, root_ ) + call check_( nf90_inq_dimid(ncid, "lev", varid), rc, mype_, root_ ) + call check_( nf90_inquire_dimension(ncid, varid, len=nlev), rc, mype_, root_ ) + +! Close the file, freeing all resources. + call check_( nf90_close(ncid), rc, mype_, root_ ) + + return + +end subroutine read_dims_ + +subroutine read_berror_ (fname,bvars,rc, myid,root) + implicit none + character(len=*), intent(in) :: fname ! input filename + type(nc_berror_vars),intent(inout) :: bvars ! background error variables + integer, intent(out) :: rc + integer, intent(in), optional :: myid,root ! accommodate MPI calling programs + +! This will be the netCDF ID for the file and data variable. + integer :: ncid, varid + +! Local variables + character(len=*), parameter :: myname_ = myname//"::read_" + character(len=4) :: cindx + integer :: nv,nl,nlat,nlon,nlev + integer :: ndims_, nvars_, ngatts_, unlimdimid_ + integer :: nlat_,nlon_,nlev_ + integer :: mype_,root_ + real(4), allocatable :: data_in(:,:,:) + logical :: verbose + logical :: init_ + + +! Return code (status) + rc=0; mype_=0; root_=0 + verbose=.true. + init_=.false. + if(present(myid).and.present(root) )then + if(myid/=root) verbose=.false. + mype_ = myid + root_ = root + endif + +! Get dimensions + call read_dims_ (fname,nlat_,nlon_,nlev_,rc, mype_,root_) + + init_ = bvars%initialized + if ( init_ ) then +! Set dims + nlat=bvars%nlat + nlon=bvars%nlon + nlev=bvars%nsig + +! Consistency check + if (nlon_ /= nlon .or. nlat_ /=nlat .or. nlev_/=nlev ) then + rc=1 + if(myid==root) then + print *, 'nlat(file) = ', nlat_, 'nlat(required) = ', nlat + print *, 'nlon(file) = ', nlon_, 'nlon(required) = ', nlon + print *, 'nlev(file) = ', nlev_, 'nlev(required) = ', nlev + print *, myname_, 'Inconsistent dimensions, aborting ... ' + endif + return + endif + else +! Set dims + nlat=nlat_ + nlon=nlon_ + nlev=nlev_ + call init_berror_vars_(bvars,nlon,nlat,nlev) + endif + +! Open the file. NF90_NOWRITE tells netCDF we want read-only access to +! the file. + + call check_( nf90_open(fname, NF90_NOWRITE, ncid), rc, mype_, root_ ) + if(rc/=0) return + +! Read global attributes +! call check_( nf90_inquire(ncid, ndims_, nvars_, ngatts_, unlimdimid_), rc, mype_, root_ ) +! call check_( nf90_inq_dimid(ncid, "lon", varid), rc, mype_, root_ ) +! call check_( nf90_inquire_dimension(ncid, varid, len=nlon_), rc, mype_, root_ ) +! call check_( nf90_inq_dimid(ncid, "lat", varid), rc, mype_, root_ ) +! call check_( nf90_inquire_dimension(ncid, varid, len=nlat_), rc, mype_, root_ ) +! call check_( nf90_inq_dimid(ncid, "lev", varid), rc, mype_, root_ ) +! call check_( nf90_inquire_dimension(ncid, varid, len=nlev_), rc, mype_, root_ ) + +! Read data to file + allocate(data_in(1,nlat,1)) + do nv = 1, nv1d + call check_( nf90_inq_varid(ncid, trim(cvars1d(nv)), varid), rc, mype_, root_ ) + call check_( nf90_get_var(ncid, varid, data_in(1,:,1)), rc, mype_, root_ ) + if(trim(cvars1d(nv))=="ps" ) bvars%psvar = data_in(1,:,1) + if(trim(cvars1d(nv))=="hps" ) bvars%pshln = data_in(1,:,1) + enddo + deallocate(data_in) + allocate(data_in(1,nlat,nlev)) + do nv = 1, nv2d + call check_( nf90_inq_varid(ncid, trim(cvars2d(nv)), varid), rc, mype_, root_ ) + call check_( nf90_get_var(ncid, varid, data_in(1,:,:)), rc, mype_, root_ ) + + if(trim(cvars2d(nv))=="sf" ) bvars%sfvar = data_in(1,:,:) + if(trim(cvars2d(nv))=="hsf") bvars%sfhln = data_in(1,:,:) + if(trim(cvars2d(nv))=="vsf") bvars%sfvln = data_in(1,:,:) +! + if(trim(cvars2d(nv))=="vp" ) bvars%vpvar = data_in(1,:,:) + if(trim(cvars2d(nv))=="hvp") bvars%vphln = data_in(1,:,:) + if(trim(cvars2d(nv))=="vvp") bvars%vpvln = data_in(1,:,:) +! + if(trim(cvars2d(nv))=="t" ) bvars%tvar = data_in(1,:,:) + if(trim(cvars2d(nv))=="ht" ) bvars%thln = data_in(1,:,:) + if(trim(cvars2d(nv))=="vt" ) bvars%tvln = data_in(1,:,:) +! + if(trim(cvars2d(nv))=="q" ) bvars%qvar = data_in(1,:,:) + if(trim(cvars2d(nv))=="hq" ) bvars%qhln = data_in(1,:,:) + if(trim(cvars2d(nv))=="vq" ) bvars%qvln = data_in(1,:,:) +! + if(trim(cvars2d(nv))=="qi" ) bvars%qivar = data_in(1,:,:) + if(trim(cvars2d(nv))=="hqi") bvars%qihln = data_in(1,:,:) + if(trim(cvars2d(nv))=="vqi") bvars%qivln = data_in(1,:,:) +! + if(trim(cvars2d(nv))=="ql" ) bvars%qlvar = data_in(1,:,:) + if(trim(cvars2d(nv))=="hql") bvars%qlhln = data_in(1,:,:) + if(trim(cvars2d(nv))=="vql") bvars%qlvln = data_in(1,:,:) +! + if(trim(cvars2d(nv))=="qr" ) bvars%qrvar = data_in(1,:,:) + if(trim(cvars2d(nv))=="hqr") bvars%qrhln = data_in(1,:,:) + if(trim(cvars2d(nv))=="vqr") bvars%qrvln = data_in(1,:,:) +! + if(trim(cvars2d(nv))=="nrh") bvars%nrhvar = data_in(1,:,:) + if(trim(cvars2d(nv))=="qs" ) bvars%qsvar = data_in(1,:,:) + if(trim(cvars2d(nv))=="hqs") bvars%qshln = data_in(1,:,:) + if(trim(cvars2d(nv))=="vqs") bvars%qsvln = data_in(1,:,:) +! + if(trim(cvars2d(nv))=="cw" ) bvars%cvar = data_in(1,:,:) + if(trim(cvars2d(nv))=="hcw") bvars%chln = data_in(1,:,:) + if(trim(cvars2d(nv))=="vcw") bvars%cvln = data_in(1,:,:) +! + if(trim(cvars2d(nv))=="oz" ) bvars%ozvar = data_in(1,:,:) + if(trim(cvars2d(nv))=="hoz") bvars%ozhln = data_in(1,:,:) + if(trim(cvars2d(nv))=="voz") bvars%ozvln = data_in(1,:,:) +! + if(trim(cvars2d(nv))=="pscon") bvars%pscon = data_in(1,:,:) + if(trim(cvars2d(nv))=="vpcon") bvars%vpcon = data_in(1,:,:) +! + enddo + +! Get matrix NLATxNLEVxNLEV that has been written as NLEV 2d-fields + do nv = 1, nvmll + do nl = 1, nlev + write(cindx,'(i4.4)') nl + if(trim(cvarsMLL(nv))=="tcon") then + call check_( nf90_inq_varid(ncid, trim(cvarsMLL(nv))//cindx, varid), rc, mype_, root_ ) + call check_( nf90_get_var(ncid, varid, data_in(1,:,:)), rc, mype_, root_ ) + bvars%tcon(:,:,nl) = data_in(1,:,:) + endif + enddo + enddo + deallocate(data_in) + +! Read in lat/lon fields + allocate(data_in(nlon,nlat,1)) + do nv = 1, nv2dx + call check_( nf90_inq_varid(ncid, trim(cvars2dx(nv)), varid), rc, mype_, root_ ) + call check_( nf90_get_var(ncid, varid, data_in(:,:,1)), rc, mype_, root_ ) + if(trim(cvars2dx(nv))=="sst" ) then + bvars%varsst = transpose(data_in(:,:,1)) + endif + if(trim(cvars2dx(nv))=="hsst" ) then + bvars%corlsst = transpose(data_in(:,:,1)) + endif + enddo + deallocate(data_in) + +! Close the file, freeing all resources. + call check_( nf90_close(ncid), rc, mype_, root_ ) + + if(verbose) print *,"*** Finish reading file: ", trim(fname) + + return + +end subroutine read_berror_ + +subroutine write_berror_ (fname,bvars,plevs,lats,lons,rc, myid,root) + implicit none + character(len=*), intent(in) :: fname ! input filename + type(nc_berror_vars),intent(in) :: bvars ! background error variables + real(4), intent(in) :: lats(:) ! latitudes per GSI: increase index from South to North Pole + real(4), intent(in) :: lons(:) ! longitudea per GSI: increase index from East to West + real(4), intent(in) :: plevs(:) + integer, intent(out) :: rc + integer, intent(in), optional :: myid,root ! accommodate MPI calling programs + + character(len=*), parameter :: myname_ = myname//"::write_" + integer, parameter :: NDIMS = 3 + +! When we create netCDF files, variables and dimensions, we get back +! an ID for each one. + character(len=4) :: cindx + integer :: ncid, dimids(NDIMS) + integer :: x_dimid, y_dimid, z_dimid + integer :: lon_varid, lat_varid, lev_varid + integer :: ii,jj,nl,nv,nn,nlat,nlon,nlev + integer :: mype_,root_ + integer, allocatable :: varid1d(:), varid2d(:), varid2dx(:), varidMLL(:) + logical :: verbose + +! This is the data array we will write. It will just be filled with +! a progression of integers for this example. + real(4), allocatable :: data_out(:,:,:) + +! Return code (status) + rc=0; mype_=0; root_=0 + verbose=.true. + if(present(myid).and.present(root) )then + if(myid/=root) verbose=.false. + mype_ = myid + root_ = root + endif + +! Set dims + nlat=bvars%nlat + nlon=bvars%nlon + nlev=bvars%nsig + +! Always check the return code of every netCDF function call. In +! this example program, wrapping netCDF calls with "call check()" +! makes sure that any return which is not equal to nf90_noerr (0) +! will print a netCDF error message and exit. + +! Create the netCDF file. The nf90_clobber parameter tells netCDF to +! overwrite this file, if it already exists. + call check_( nf90_create(fname, NF90_CLOBBER, ncid), rc, mype_, root_ ) + if(rc/=0) return + +! Define the dimensions. NetCDF will hand back an ID for each. + call check_( nf90_def_dim(ncid, "lon", nlon, x_dimid), rc, mype_, root_ ) + call check_( nf90_def_dim(ncid, "lat", nlat, y_dimid), rc, mype_, root_ ) + call check_( nf90_def_dim(ncid, "lev", nlev, z_dimid), rc, mype_, root_ ) + + call check_( nf90_def_var(ncid, "lon", NF90_REAL, x_dimid, lon_varid), rc, mype_, root_ ) + call check_( nf90_def_var(ncid, "lat", NF90_REAL, y_dimid, lat_varid), rc, mype_, root_ ) + call check_( nf90_def_var(ncid, "lev", NF90_REAL, z_dimid, lev_varid), rc, mype_, root_ ) + + call check_( nf90_put_att(ncid, lon_varid, "units", "degress"), rc, mype_, root_ ) + call check_( nf90_put_att(ncid, lat_varid, "units", "degress"), rc, mype_, root_ ) + call check_( nf90_put_att(ncid, lev_varid, "units", "hPa"), rc, mype_, root_ ) + +! The dimids array is used to pass the IDs of the dimensions of +! the variables. Note that in fortran arrays are stored in +! column-major format. + dimids = (/ x_dimid, y_dimid, z_dimid /) + +! Define variables. + allocate(varid1d(nv1d)) + do nv = 1, nv1d + call check_( nf90_def_var(ncid, trim(cvars1d(nv)), NF90_REAL, (/ y_dimid /), varid1d(nv)), rc, mype_, root_ ) + enddo + allocate(varid2d(nv2d)) + do nv = 1, nv2d + call check_( nf90_def_var(ncid, trim(cvars2d(nv)), NF90_REAL, (/ y_dimid, z_dimid /), varid2d(nv)), rc, mype_, root_ ) + enddo + allocate(varidMLL(nlev*nvmll)) + nn=0 + do nv = 1, nvmll + do nl = 1, nlev + nn=nn+1 + write(cindx,'(i4.4)') nl + call check_( nf90_def_var(ncid, trim(cvarsMLL(nv))//cindx, NF90_REAL, (/ y_dimid, z_dimid /), varidMLL(nn)), rc, & + mype_, root_ ) + enddo + enddo + allocate(varid2dx(nv2dx)) + do nv = 1, nv2dx + call check_( nf90_def_var(ncid, trim(cvars2dx(nv)), NF90_REAL, (/ x_dimid, y_dimid /), varid2dx(nv)), rc, mype_, root_ ) + enddo + +! End define mode. This tells netCDF we are done defining metadata. + call check_( nf90_enddef(ncid), rc, mype_, root_ ) + +! Write coordinate variables data + call check_( nf90_put_var(ncid, lon_varid, lons ), rc, mype_, root_ ) + call check_( nf90_put_var(ncid, lat_varid, lats ), rc, mype_, root_ ) + call check_( nf90_put_var(ncid, lev_varid, plevs), rc, mype_, root_ ) + +! Write data to file + allocate(data_out(1,nlat,1)) + do nv = 1, nv1d + if(trim(cvars1d(nv))=="ps" ) data_out(1,:,1) = bvars%psvar + if(trim(cvars1d(nv))=="hps" ) data_out(1,:,1) = bvars%pshln + call check_( nf90_put_var(ncid, varid1d(nv), data_out(1,:,1)), rc, mype_, root_) + enddo + deallocate(data_out) + allocate(data_out(1,nlat,nlev)) + do nv = 1, nv2d + if(trim(cvars2d(nv))=="sf" ) data_out(1,:,:) = bvars%sfvar + if(trim(cvars2d(nv))=="hsf") data_out(1,:,:) = bvars%sfhln + if(trim(cvars2d(nv))=="vsf") data_out(1,:,:) = bvars%sfvln +! + if(trim(cvars2d(nv))=="vp" ) data_out(1,:,:) = bvars%vpvar + if(trim(cvars2d(nv))=="hvp") data_out(1,:,:) = bvars%vphln + if(trim(cvars2d(nv))=="vvp") data_out(1,:,:) = bvars%vpvln +! + if(trim(cvars2d(nv))=="t" ) data_out(1,:,:) = bvars%tvar + if(trim(cvars2d(nv))=="ht" ) data_out(1,:,:) = bvars%thln + if(trim(cvars2d(nv))=="vt" ) data_out(1,:,:) = bvars%tvln +! + if(trim(cvars2d(nv))=="q" ) data_out(1,:,:) = bvars%qvar + if(trim(cvars2d(nv))=="hq" ) data_out(1,:,:) = bvars%qhln + if(trim(cvars2d(nv))=="vq" ) data_out(1,:,:) = bvars%qvln +! + if(trim(cvars2d(nv))=="qi" ) data_out(1,:,:) = bvars%qivar + if(trim(cvars2d(nv))=="hqi") data_out(1,:,:) = bvars%qihln + if(trim(cvars2d(nv))=="vqi") data_out(1,:,:) = bvars%qivln +! + if(trim(cvars2d(nv))=="ql" ) data_out(1,:,:) = bvars%qlvar + if(trim(cvars2d(nv))=="hql") data_out(1,:,:) = bvars%qlhln + if(trim(cvars2d(nv))=="vql") data_out(1,:,:) = bvars%qlvln +! + if(trim(cvars2d(nv))=="qr" ) data_out(1,:,:) = bvars%qrvar + if(trim(cvars2d(nv))=="hqr") data_out(1,:,:) = bvars%qrhln + if(trim(cvars2d(nv))=="vqr") data_out(1,:,:) = bvars%qrvln +! + if(trim(cvars2d(nv))=="nrh") data_out(1,:,:) = bvars%nrhvar + if(trim(cvars2d(nv))=="qs" ) data_out(1,:,:) = bvars%qsvar + if(trim(cvars2d(nv))=="hqs") data_out(1,:,:) = bvars%qshln + if(trim(cvars2d(nv))=="vqs") data_out(1,:,:) = bvars%qsvln +! + if(trim(cvars2d(nv))=="cw" ) data_out(1,:,:) = bvars%cvar + if(trim(cvars2d(nv))=="hcw") data_out(1,:,:) = bvars%chln + if(trim(cvars2d(nv))=="vcw") data_out(1,:,:) = bvars%cvln +! + if(trim(cvars2d(nv))=="oz" ) data_out(1,:,:) = bvars%ozvar + if(trim(cvars2d(nv))=="hoz") data_out(1,:,:) = bvars%ozhln + if(trim(cvars2d(nv))=="voz") data_out(1,:,:) = bvars%ozvln +! + if(trim(cvars2d(nv))=="pscon") data_out(1,:,:) = bvars%pscon + if(trim(cvars2d(nv))=="vpcon") data_out(1,:,:) = bvars%vpcon +! + call check_( nf90_put_var(ncid, varid2d(nv), data_out(1,:,:)), rc, mype_, root_ ) + enddo + +! Choose to write out NLATxNLEVxNLEV vars as to facilitate visualization + nn=0 + do nv = 1, nvmll + do nl = 1, nlev + nn = nn + 1 + write(cindx,'(i4.4)') nl + if(trim(cvarsMLL(nv))=="tcon") data_out(1,:,:) = bvars%tcon(:,:,nl) + call check_( nf90_put_var(ncid, varidMLL(nn), data_out(1,:,:)), rc, mype_, root_ ) + enddo + enddo + deallocate(data_out) + +! Write out lat/lon fields + allocate(data_out(nlon,nlat,1)) + do nv = 1, nv2dx + if(trim(cvars2dx(nv))=="sst" ) then + data_out(:,:,1) = transpose(bvars%varsst) + endif + if(trim(cvars2dx(nv))=="hsst" ) then + data_out(:,:,1) = transpose(bvars%corlsst) + endif + call check_( nf90_put_var(ncid, varid2dx(nv), data_out(:,:,1)), rc, mype_, root_ ) + enddo + deallocate(data_out) + +! Close file + call check_( nf90_close(ncid), rc, mype_, root_ ) + + deallocate(varidMLL) + deallocate(varid2d) + deallocate(varid1d) + + if(verbose) print *,"*** Finish writing file: ", trim(fname) + + return + +end subroutine write_berror_ + +subroutine init_berror_vars_(vr,nlon,nlat,nsig) + + integer,intent(in) :: nlon,nlat,nsig + type(nc_berror_vars) vr + + if(vr%initialized) return + + vr%nlon=nlon + vr%nlat=nlat + vr%nsig=nsig + +! allocate single precision arrays + allocate(vr%sfvar(nlat,nsig),vr%vpvar(nlat,nsig),vr%tvar(nlat,nsig),vr%qvar(nlat,nsig), & + vr%qivar(nlat,nsig),vr%qlvar(nlat,nsig),vr%qrvar(nlat,nsig),vr%qsvar(nlat,nsig),& + vr%cvar(nlat,nsig),vr%nrhvar(nlat,nsig),vr%ozvar(nlat,nsig)) + allocate(vr%sfhln(nlat,nsig),vr%vphln(nlat,nsig),vr%thln(nlat,nsig),vr%qhln(nlat,nsig), & + vr%qihln(nlat,nsig),vr%qlhln(nlat,nsig),vr%qrhln(nlat,nsig),vr%qshln(nlat,nsig),& + vr%chln(nlat,nsig), vr%ozhln(nlat,nsig)) + allocate(vr%sfvln(nlat,nsig),vr%vpvln(nlat,nsig),vr%tvln(nlat,nsig),vr%qvln(nlat,nsig), & + vr%qivln(nlat,nsig),vr%qlvln(nlat,nsig),vr%qrvln(nlat,nsig),vr%qsvln(nlat,nsig),& + vr%cvln(nlat,nsig), vr%ozvln(nlat,nsig)) + allocate(vr%pscon(nlat,nsig),vr%vpcon(nlat,nsig)) + allocate(vr%varsst(nlat,nlon),vr%corlsst(nlat,nlon)) + allocate(vr%tcon(nlat,nsig,nsig)) + allocate(vr%psvar(nlat),vr%pshln(nlat)) + vr%initialized=.true. + end subroutine init_berror_vars_ + + subroutine final_berror_vars_(vr) + type(nc_berror_vars) vr +! deallocate arrays + if(.not. vr%initialized) return + deallocate(vr%sfvar,vr%vpvar,vr%tvar,vr%qvar, & + vr%qivar,vr%qlvar,vr%qrvar,vr%qsvar,& + vr%cvar,vr%nrhvar,vr%ozvar) + deallocate(vr%sfhln,vr%vphln,vr%thln,vr%qhln, & + vr%qihln,vr%qlhln,vr%qrhln,vr%qshln,& + vr%chln, vr%ozhln) + deallocate(vr%sfvln,vr%vpvln,vr%tvln,vr%qvln, & + vr%qivln,vr%qlvln,vr%qrvln,vr%qsvln,& + vr%cvln, vr%ozvln) + deallocate(vr%pscon,vr%vpcon) + deallocate(vr%varsst,vr%corlsst) + deallocate(vr%tcon) + deallocate(vr%psvar,vr%pshln) + vr%initialized=.false. +end subroutine final_berror_vars_ + +subroutine comp_berror_vars_(va,vb,rc, myid,root) + type(nc_berror_vars) va + type(nc_berror_vars) vb + integer, intent(out) :: rc + integer, intent(in), optional :: myid,root ! accommodate MPI calling programs + character(len=*), parameter :: myname_ = myname//"::comp_berror_vars_" + integer :: ii,jj,ier(50) + logical :: verbose, failed + real :: tolerance = 10.e-10 +! + rc=0 + verbose=.true. + if(present(myid).and.present(root) )then + if(myid/=root) verbose=.false. + endif +! Consistency check + if (va%nlon/=vb%nlon .or. va%nlat/=vb%nlat .or. va%nsig/=vb%nsig ) then + rc=1 + if(myid==root) then + print *, 'nlat(va) = ', va%nlat, 'nlat(vb) = ', vb%nlat + print *, 'nlon(va) = ', va%nlon, 'nlon(vb) = ', vb%nlon + print *, 'nlev(va) = ', va%nsig, 'nlev(vb) = ', vb%nsig + print *, myname_, 'Inconsistent dimensions, aborting ... ' + endif + return + endif + + ii=0;ier=0 + ii=ii+1; if(abs(sum(va%sfvar - vb%sfvar)) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%vpvar - vb%vpvar)) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%tvar - vb%tvar)) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%qvar - vb%qvar )) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%qivar - vb%qivar)) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%qlvar - vb%qlvar)) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%qrvar - vb%qrvar)) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%qsvar - vb%qsvar) )>tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%cvar - vb%cvar )) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%nrhvar- vb%nrhvar))>tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%ozvar - vb%ozvar)) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%sfhln - vb%sfhln)) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%vphln - vb%vphln ))>tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%thln - vb%thln)) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%qhln - vb%qhln) ) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%qihln - vb%qihln)) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%qlhln - vb%qlhln)) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%qrhln - vb%qrhln) )>tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%qshln - vb%qshln ))>tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%chln - vb%chln )) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%ozhln - vb%ozhln)) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%sfvln - vb%sfvln)) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%vpvln - vb%vpvln)) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%tvln - vb%tvln)) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%qvln - vb%qvln )) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%qivln - vb%qivln)) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%qlvln - vb%qlvln)) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%qrvln - vb%qrvln)) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%qsvln - vb%qsvln) )>tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%cvln - vb%cvln )) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%ozvln - vb%ozvln)) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%pscon - vb%pscon)) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%vpcon - vb%vpcon)) >tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%varsst- vb%varsst))>tolerance) ier(ii)=ii + ii=ii+1; if(abs(sum(va%corlsst-vb%corlsst))>tolerance)ier(ii)=ii + ii=ii+1; if(abs(sum(va%tcon - vb%tcon)) >tolerance) ier(ii)=ii + failed=.false. + do jj=1,ii + if(ier(jj)/=0.and.verbose) then + print *, 'Found field ', jj, ' not to match' + failed=.true. + endif + enddo + if (.not.failed) then + if(verbose) print *, 'Comp finds all fields to match' + endif +end subroutine comp_berror_vars_ + +subroutine copy_(ivars,ovars,hydro) + type(nc_berror_vars) ivars + type(nc_berror_vars) ovars + logical, intent(in), optional :: hydro + + logical wrtall,hydro_ + + hydro_=.true. + wrtall=.true. + if (ovars%nlon/=ivars%nlon .or. & + ovars%nlat/=ivars%nlat ) then + print*, 'copy_berror_vars_: Trying to copy inconsistent vectors, aborting ...' + call exit(1) + endif + if ( ovars%nsig/=ivars%nsig ) then + wrtall=.false. + endif + if(present(hydro)) then + hydro_ = hydro + endif + + if (wrtall) then + ovars%tcon = ivars%tcon + ovars%vpcon = ivars%vpcon + ovars%pscon = ivars%pscon + ovars%sfvar = ivars%sfvar + ovars%sfhln = ivars%sfhln + ovars%sfvln = ivars%sfvln + ovars%vpvar = ivars%vpvar + ovars%vphln = ivars%vphln + ovars%vpvln = ivars%vpvln + ovars%tvar = ivars%tvar + ovars%thln = ivars%thln + ovars%tvln = ivars%tvln + ovars%qvar = ivars%qvar + ovars%nrhvar = ivars%nrhvar + ovars%qhln = ivars%qhln + ovars%qvln = ivars%qvln + if(hydro_) then + ovars%qivar = ivars%qivar + ovars%qihln = ivars%qihln + ovars%qivln = ivars%qivln + ovars%qlvar = ivars%qlvar + ovars%qlhln = ivars%qlhln + ovars%qlvln = ivars%qlvln + ovars%qrvar = ivars%qrvar + ovars%qrhln = ivars%qrhln + ovars%qrvln = ivars%qrvln + ovars%qsvar = ivars%qsvar + ovars%qshln = ivars%qshln + ovars%qsvln = ivars%qsvln + endif + ovars%ozvar = ivars%ozvar + ovars%ozhln = ivars%ozhln + ovars%ozvln = ivars%ozvln + ovars%cvar = ivars%cvar + ovars%chln = ivars%chln + ovars%cvln = ivars%cvln + endif + + ovars%psvar = ivars%psvar + ovars%pshln = ivars%pshln + ovars%varsst = ivars%varsst + ovars%corlsst = ivars%corlsst + +end subroutine copy_ + +subroutine get_pointer_1d_ (vname, bvars, ptr, rc ) +implicit none +character(len=*), intent(in) :: vname +type(nc_berror_vars) bvars +real(4),pointer,intent(inout) :: ptr(:) +integer,intent(out) :: rc +rc=-1 +if(trim(vname)=='ps') then + ptr => bvars%psvar + rc=0 +endif +if(trim(vname)=='hps') then + ptr => bvars%pshln + rc=0 +endif +end subroutine get_pointer_1d_ + +subroutine get_pointer_2d_ (vname, bvars, ptr, rc ) +implicit none +character(len=*), intent(in) :: vname +type(nc_berror_vars) bvars +real(4),pointer,intent(inout) :: ptr(:,:) +integer,intent(out) :: rc +character(len=5) :: var +rc=-1 +! +var='sst' +if(trim(vname)==trim(var)) then + ptr => bvars%varsst + rc=0 + return +endif +if(trim(vname)=='h'//trim(var)) then + ptr => bvars%corlsst + rc=0 + return +endif +! +var='sf' +if(trim(vname)==trim(var)) then + ptr => bvars%sfvar + rc=0 + return +endif +if(trim(vname)=='h'//trim(var)) then + ptr => bvars%sfhln + rc=0 + return +endif +if(trim(vname)=='v'//trim(var)) then + ptr => bvars%sfvln + rc=0 + return +endif +! +var='vp' +if(trim(vname)==trim(var)) then + ptr => bvars%vpvar + rc=0 + return +endif +if(trim(vname)=='h'//trim(var)) then + ptr => bvars%vphln + rc=0 + return +endif +if(trim(vname)=='v'//trim(var)) then + ptr => bvars%vpvln + rc=0 + return +endif +! +var='t' +if(trim(vname)==trim(var)) then + ptr => bvars%tvar + rc=0 + return +endif +if(trim(vname)=='h'//trim(var)) then + ptr => bvars%thln + rc=0 + return +endif +if(trim(vname)=='v'//trim(var)) then + ptr => bvars%tvln + rc=0 + return +endif +! +var='q' +if(trim(vname)==trim(var)) then + ptr => bvars%qvar + rc=0 + return +endif +if(trim(vname)=='h'//trim(var)) then + ptr => bvars%qhln + rc=0 + return +endif +if(trim(vname)=='v'//trim(var)) then + ptr => bvars%qvln + rc=0 + return +endif +! +var='cw' +if(trim(vname)==trim(var)) then + ptr => bvars%cvar + rc=0 + return +endif +if(trim(vname)=='h'//trim(var)) then + ptr => bvars%chln + rc=0 + return +endif +if(trim(vname)=='v'//trim(var)) then + ptr => bvars%cvln + rc=0 + return +endif +! +var='qi' +if(trim(vname)==trim(var)) then + ptr => bvars%qivar + rc=0 + return +endif +if(trim(vname)=='h'//trim(var)) then + ptr => bvars%qihln + rc=0 + return +endif +if(trim(vname)=='v'//trim(var)) then + ptr => bvars%qivln + rc=0 + return +endif +! +var='ql' +if(trim(vname)==trim(var)) then + ptr => bvars%qlvar + rc=0 + return +endif +if(trim(vname)=='h'//trim(var)) then + ptr => bvars%qlhln + rc=0 + return +endif +if(trim(vname)=='v'//trim(var)) then + ptr => bvars%qlvln + rc=0 + return +endif +! +var='qr' +if(trim(vname)==trim(var)) then + ptr => bvars%qrvar + rc=0 + return +endif +if(trim(vname)=='h'//trim(var)) then + ptr => bvars%qrhln + rc=0 + return +endif +if(trim(vname)=='v'//trim(var)) then + ptr => bvars%qrvln + rc=0 + return +endif +! +var='qs' +if(trim(vname)==trim(var)) then + ptr => bvars%qsvar + rc=0 + return +endif +if(trim(vname)=='h'//trim(var)) then + ptr => bvars%qshln + rc=0 + return +endif +if(trim(vname)=='v'//trim(var)) then + ptr => bvars%qsvln + rc=0 + return +endif +! +var='oz' +if(trim(vname)==trim(var)) then + ptr => bvars%ozvar + rc=0 + return +endif +if(trim(vname)=='h'//trim(var)) then + ptr => bvars%ozhln + rc=0 + return +endif +if(trim(vname)=='v'//trim(var)) then + ptr => bvars%ozvln + rc=0 + return +endif +! +var='nrh' +if(trim(vname)==trim(var)) then + ptr => bvars%nrhvar + rc=0 + return +endif +end subroutine get_pointer_2d_ + +subroutine check_(status,rc, myid, root) + integer, intent ( in) :: status + integer, intent (out) :: rc + integer, intent ( in) :: myid, root + rc=0 + if(status /= nf90_noerr) then + if(myid==root) print *, trim(nf90_strerror(status)) + rc=999 + end if +end subroutine check_ + +end module m_nc_berror diff --git a/src/gsi/ncepnems_io.f90 b/src/gsi/ncepnems_io.f90 index 595f07e152..5a8b80d71b 100755 --- a/src/gsi/ncepnems_io.f90 +++ b/src/gsi/ncepnems_io.f90 @@ -1805,8 +1805,16 @@ subroutine read_hsst_(hsst,cwoption) ! use kinds, only : i_kind,r_single,r_kind use gridmod, only : nlat,nlon,nsig - use mpeu_util, only : check_iostat - + use mpeu_util, only : check_iostat,die + use module_ncio, only: open_dataset, Dataset, Dimension, & + close_dataset, read_vardata, get_dim + use control_vectors,only: cvars2d + use m_nc_berror, only: nc_berror_read + use m_nc_berror, only: nc_berror_vars + use m_nc_berror, only: nc_berror_getpointer + use m_nc_berror, only: nc_berror_vars_final + use mpimod, only: mype + implicit none integer(i_kind) , intent(in ) :: cwoption @@ -1830,72 +1838,109 @@ subroutine read_hsst_(hsst,cwoption) character(len=256) :: berror_stats = "berror_stats" ! filename - ! Open background error statistics file - inerr = 22 - open(inerr,file=berror_stats,form='unformatted',status='old',iostat=ier) - call check_iostat(ier,myname_,'open('//trim(berror_stats)//')') - - ! Read header. Ensure that vertical resolution is consistent - ! with that specified via the user namelist - - rewind inerr - read(inerr,iostat=ier)nsigstat,nlatstat - call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (nsigstat,nlatstat)') - - if ( nsigstat/=nsig .or. nlatstat/=nlat ) then - write(6,*)'PREBAL: **ERROR** resolution of berror_stats incompatiable with GSI' - - write(6,*)'PREBAL: berror nsigstat,nlatstat=', nsigstat,nlatstat, & - ' -vs- GSI nsig,nlat=',nsig,nlat - call stop2(101) + logical :: ncio + type(Dataset) :: dset + type(Dimension) :: londim,latdim,levdim + integer(i_kind) :: nv,n + type(nc_berror_vars) bvars + real(r_single), pointer :: ptr2d(:,:) + + dset = open_dataset(berror_stats,errcode=ier) + if (ier==0) then + ! this is a netcdf file + ncio = .true. + else + ! this is a binary file + ncio = .false. endif - - write(6,*) myname_,'(read_hsst): read error amplitudes ', & - '"',trim(berror_stats),'". ', & - 'nsigstat,nlatstat =',nsigstat,nlatstat - - read(inerr,iostat=ier) agvin,bvin,wgvin - call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (agvin,bvin,wgvin)') - - write(*,*) 'in read_hsst_, cwoption = ',cwoption - - readloop: do - read(inerr,iostat=istat) var, isig - if ( istat/=0 ) exit - write(*,*) 'var, isig, istat : ',var, isig, istat - allocate(corzin(nlat,isig)) - if ( var=='q' .or. var=='cw' ) allocate(corq2(nlat,isig)) - allocate(hwllin(nlat,isig)) - if ( isig>1 ) allocate(vscalesin(nlat,isig)) - if ( var/='sst' ) then - if ( var=='q' .or. var=='Q' .or. (var=='cw' .and. cwoption==2) ) then - read(inerr,iostat=ier) corzin,corq2 - call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (corzin,corq2)') - else - read(inerr,iostat=ier) corzin - call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (corzin)') - endif - read(inerr,iostat=ier) hwllin - call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (hwllin)') - if ( isig>1 ) then - read(inerr,iostat=ier) vscalesin - call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (vscalein)') + + if (ncio) then + call nc_berror_read (berror_stats,bvars,ier, myid=mype,root=0) + if (nlat/=bvars%nlat .or. nlon/=bvars%nlon .or. nsig/=bvars%nsig ) then + call die(myname_," inconsistent dims in "//trim(berror_stats), 99) + endif + isig=bvars%nsig + + ! RTodling: the following is bad since it wires all naming conventions ... to be revised + do nv=1,size(cvars2d) + if (trim(cvars2d(nv))=='sst') then +!! n = getindex(cvars2d,'sst') +!! found2d(n)=.true. +!! call nc_berror_getpointer (cvars2d(nv),bvars,ptr2d,ier) +!! if(ier==0) corsst=ptr2d + call nc_berror_getpointer ('h'//cvars2d(nv),bvars,ptr2d,ier) + if(ier==0) hsst=ptr2d endif - else - read(inerr,iostat=ier) corsst - call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (corsst)') - read(inerr,iostat=ier) hsst - call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (hsst)') + enddo + call nc_berror_vars_final(bvars) + else + + ! Open background error statistics file + inerr = 22 + open(inerr,file=berror_stats,form='unformatted',status='old',iostat=ier) + call check_iostat(ier,myname_,'open('//trim(berror_stats)//')') + + ! Read header. Ensure that vertical resolution is consistent + ! with that specified via the user namelist + + rewind inerr + read(inerr,iostat=ier)nsigstat,nlatstat + call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (nsigstat,nlatstat)') + + if ( nsigstat/=nsig .or. nlatstat/=nlat ) then + write(6,*)'PREBAL: **ERROR** resolution of berror_stats incompatiable with GSI' + + write(6,*)'PREBAL: berror nsigstat,nlatstat=', nsigstat,nlatstat, & + ' -vs- GSI nsig,nlat=',nsig,nlat + call stop2(101) endif - - - deallocate(corzin,hwllin) - if (isig>1) deallocate(vscalesin) - if ( var=='q' .or. var=='cw' ) deallocate(corq2) - - enddo readloop - close(inerr) - + + write(6,*) myname_,'(read_hsst): read error amplitudes ', & + '"',trim(berror_stats),'". ', & + 'nsigstat,nlatstat =',nsigstat,nlatstat + + read(inerr,iostat=ier) agvin,bvin,wgvin + call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (agvin,bvin,wgvin)') + + write(*,*) 'in read_hsst_, cwoption = ',cwoption + + readloop: do + read(inerr,iostat=istat) var, isig + if ( istat/=0 ) exit + write(*,*) 'var, isig, istat : ',var, isig, istat + allocate(corzin(nlat,isig)) + if ( var=='q' .or. var=='cw' ) allocate(corq2(nlat,isig)) + allocate(hwllin(nlat,isig)) + if ( isig>1 ) allocate(vscalesin(nlat,isig)) + if ( var/='sst' ) then + if ( var=='q' .or. var=='Q' .or. (var=='cw' .and. cwoption==2) ) then + read(inerr,iostat=ier) corzin,corq2 + call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (corzin,corq2)') + else + read(inerr,iostat=ier) corzin + call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (corzin)') + endif + read(inerr,iostat=ier) hwllin + call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (hwllin)') + if ( isig>1 ) then + read(inerr,iostat=ier) vscalesin + call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (vscalein)') + endif + else + read(inerr,iostat=ier) corsst + call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (corsst)') + read(inerr,iostat=ier) hsst + call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (hsst)') + endif + + + deallocate(corzin,hwllin) + if (isig>1) deallocate(vscalesin) + if ( var=='q' .or. var=='cw' ) deallocate(corq2) + + enddo readloop + close(inerr) + endif return end subroutine read_hsst_ From 5f348de35e8c4393ac0a7db75f48df067d0d20e8 Mon Sep 17 00:00:00 2001 From: "russ.treadon" Date: Mon, 25 Nov 2024 20:37:22 +0000 Subject: [PATCH 2/9] add missing file to gsi_files.cmake (#808) --- src/gsi/gsi_files.cmake | 1 + 1 file changed, 1 insertion(+) diff --git a/src/gsi/gsi_files.cmake b/src/gsi/gsi_files.cmake index d4b6e6b3ed..5c5c4ec5a5 100644 --- a/src/gsi/gsi_files.cmake +++ b/src/gsi/gsi_files.cmake @@ -498,6 +498,7 @@ read_goesndr.f90 read_gps.f90 read_guess.F90 read_iasi.f90 +read_iasing.f90 read_l2bufr_mod.f90 read_lag.f90 read_lidar.f90 From c484d4c88ae4b42c295dcca8b6f97d894e4ff320 Mon Sep 17 00:00:00 2001 From: "russ.treadon" Date: Tue, 26 Nov 2024 20:20:53 +0000 Subject: [PATCH 3/9] remove extraneous prints from m_berror_stats.f90 (#808) --- src/gsi/m_berror_stats.f90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/gsi/m_berror_stats.f90 b/src/gsi/m_berror_stats.f90 index 8ca7b44aea..cd66b570f2 100644 --- a/src/gsi/m_berror_stats.f90 +++ b/src/gsi/m_berror_stats.f90 @@ -492,10 +492,8 @@ subroutine read_wgt(corz,corp,hwll,hwllp,vz,corsst,hsst,varq,qoption,varcw,cwopt call berror_get_dims(msig,mlat) if ( bin_berror ) then - write(6,*)'call binary berror read' call bin_() else - write(6,*)'call netcdf berror read' call nc_(mype) endif From 0a3bb05d27cf3cbfd07d368828e3a57735078403 Mon Sep 17 00:00:00 2001 From: "russ.treadon" Date: Mon, 2 Dec 2024 15:12:14 +0000 Subject: [PATCH 4/9] update fix submodule hash (#808) --- fix | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fix b/fix index 2a86d5b1a0..1567130f4e 160000 --- a/fix +++ b/fix @@ -1 +1 @@ -Subproject commit 2a86d5b1a09a8908618d1c8a376c68e8d9b3d02c +Subproject commit 1567130f4eba20b301903e30114e72a734fd9a15 From 68d52ac65d08b51e46cdcfac8fb59939ebad28b5 Mon Sep 17 00:00:00 2001 From: "russ.treadon" Date: Mon, 2 Dec 2024 16:28:32 +0000 Subject: [PATCH 5/9] code cleanup (#808) --- src/gsi/m_berror_stats.f90 | 23 ++-- src/gsi/m_nc_berror.f90 | 209 ------------------------------------- 2 files changed, 7 insertions(+), 225 deletions(-) diff --git a/src/gsi/m_berror_stats.f90 b/src/gsi/m_berror_stats.f90 index cd66b570f2..9a97a199ca 100644 --- a/src/gsi/m_berror_stats.f90 +++ b/src/gsi/m_berror_stats.f90 @@ -94,9 +94,6 @@ module m_berror_stats logical,save :: bin_berror=.false. -!! real(r_kind),allocatable,dimension(:,:):: varq -!! real(r_kind),allocatable,dimension(:,:):: varcw - contains !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -109,13 +106,11 @@ module m_berror_stats ! ! !INTERFACE: -!!subroutine get_dims(mype,msig,mlat,mlon,lunit) subroutine get_dims(msig,mlat,mlon,lunit) use mpimod, only: mype implicit none -!! integer(i_kind) ,intent( in) :: mype ! proc identifier integer(i_kind) ,intent( out) :: msig ! dimension of levels integer(i_kind) ,intent( out) :: mlat ! dimension of latitudes integer(i_kind),optional,intent( out) :: mlon ! dimension of longitudes @@ -417,7 +412,6 @@ end subroutine read_bal ! !INTERFACE: subroutine read_wgt(corz,corp,hwll,hwllp,vz,corsst,hsst,varq,qoption,varcw,cwoption,mype,lunit) -!! n_clouds_fwd,cloud_names_fwd,lunit) use kinds, only : r_single,r_kind use gridmod,only : nlat,nlon,nsig @@ -444,8 +438,6 @@ subroutine read_wgt(corz,corp,hwll,hwllp,vz,corsst,hsst,varq,qoption,varcw,cwopt ! Optionals integer(i_kind),optional ,intent(in ) :: lunit ! an alternative unit -!! integer(i_kind), optional ,intent(in ) :: n_clouds_fwd -!! character(len=*),optional ,intent(in ) :: cloud_names_fwd(:) ! !REVISION HISTORY: ! 30Jul08 - Jing Guo @@ -644,7 +636,6 @@ subroutine bin_ do i=1,nlat corq2x=corq2(i,k) varq(i,k)=min(max(corq2x,0.00015_r_kind),one) -!! varq_out(i,k) = varq(i,k) enddo enddo do k=1,isig @@ -845,7 +836,7 @@ subroutine setcoroz_(coroz,mype) if ( ierror/=0 ) return ! nothing to do ! sanity check - if ( mype==0 ) write(6,*) myname_,'(PREWGT): mype = ',mype + if ( mype==0 ) write(6,*) myname_,'(PREWGT): enter routine' mlat=size(coroz,1) msig=size(coroz,2) @@ -871,7 +862,7 @@ subroutine setcoroz_(coroz,mype) enddo enddo enddo - work_oz(nsig+1,mm1)=float(lon1*lat1) + work_oz(nsig+1,mm1)=real(lon1*lat1,r_kind) call mpi_allreduce(work_oz,work_oz1,(nsig+1)*npe,mpi_rtype,mpi_sum,& mpi_comm_world,ierror) @@ -944,13 +935,13 @@ subroutine sethwlloz_(hwlloz,mype) real(r_kind) :: fact real(r_kind) :: s2u - if ( mype==0 ) write(6,*) myname_,'(PREWGT): mype = ',mype + if ( mype==0 ) write(6,*) myname_,'(PREWGT): enter routine' s2u=(two*pi*rearth_equator)/nlon do k=1,nnnn1o k1=levs_id(k) if ( k1>0 ) then - if(mype==0) write(6,*) myname_,'(PREWGT): mype = ',mype, k1 + if(mype==0) write(6,*) myname_,'(PREWGT): k1 = ',k1 if ( k1<=nsig*3/4 ) then ! fact=1./hwl fact=r40000/(r400*nlon) @@ -962,7 +953,7 @@ subroutine sethwlloz_(hwlloz,mype) endif enddo - if ( mype==0 ) write(6,*) myname_,'(PREWGT): mype = ',mype, 'finish sethwlloz_' + if ( mype==0 ) write(6,*) myname_,'(PREWGT): finish sethwlloz_' return end subroutine sethwlloz_ @@ -1055,7 +1046,7 @@ subroutine setcorchem_(cname,corchem,rc) rc=0 ! sanity check - if ( mype==0 ) write(6,*) myname_,'(PREWGT): mype = ',mype + if ( mype==0 ) write(6,*) myname_,'(PREWGT): enter routine' ! Get information for how to use CO2 iptr=-1 @@ -1096,7 +1087,7 @@ subroutine setcorchem_(cname,corchem,rc) enddo enddo enddo - work_chem(nsig+1,mm1)=float(lon1*lat1) + work_chem(nsig+1,mm1)=real(lon1*lat1,r_kind) call mpi_allreduce(work_chem,work_chem1,(nsig+1)*npe,mpi_rtype,mpi_sum,& mpi_comm_world,ierror) diff --git a/src/gsi/m_nc_berror.f90 b/src/gsi/m_nc_berror.f90 index 559690112d..0f3d5b329b 100644 --- a/src/gsi/m_nc_berror.f90 +++ b/src/gsi/m_nc_berror.f90 @@ -10,7 +10,6 @@ module m_nc_berror public :: nc_berror_vars public :: nc_berror_dims public :: nc_berror_read -public :: nc_berror_write public :: nc_berror_getpointer type nc_berror_vars @@ -60,8 +59,6 @@ module m_nc_berror read_dims_ ; end interface interface nc_berror_read; module procedure & read_berror_ ; end interface -interface nc_berror_write; module procedure & - write_berror_ ; end interface interface nc_berror_vars_init; module procedure & init_berror_vars_ ; end interface interface nc_berror_vars_final; module procedure & @@ -185,15 +182,6 @@ subroutine read_berror_ (fname,bvars,rc, myid,root) call check_( nf90_open(fname, NF90_NOWRITE, ncid), rc, mype_, root_ ) if(rc/=0) return -! Read global attributes -! call check_( nf90_inquire(ncid, ndims_, nvars_, ngatts_, unlimdimid_), rc, mype_, root_ ) -! call check_( nf90_inq_dimid(ncid, "lon", varid), rc, mype_, root_ ) -! call check_( nf90_inquire_dimension(ncid, varid, len=nlon_), rc, mype_, root_ ) -! call check_( nf90_inq_dimid(ncid, "lat", varid), rc, mype_, root_ ) -! call check_( nf90_inquire_dimension(ncid, varid, len=nlat_), rc, mype_, root_ ) -! call check_( nf90_inq_dimid(ncid, "lev", varid), rc, mype_, root_ ) -! call check_( nf90_inquire_dimension(ncid, varid, len=nlev_), rc, mype_, root_ ) - ! Read data to file allocate(data_in(1,nlat,1)) do nv = 1, nv1d @@ -290,203 +278,6 @@ subroutine read_berror_ (fname,bvars,rc, myid,root) end subroutine read_berror_ -subroutine write_berror_ (fname,bvars,plevs,lats,lons,rc, myid,root) - implicit none - character(len=*), intent(in) :: fname ! input filename - type(nc_berror_vars),intent(in) :: bvars ! background error variables - real(4), intent(in) :: lats(:) ! latitudes per GSI: increase index from South to North Pole - real(4), intent(in) :: lons(:) ! longitudea per GSI: increase index from East to West - real(4), intent(in) :: plevs(:) - integer, intent(out) :: rc - integer, intent(in), optional :: myid,root ! accommodate MPI calling programs - - character(len=*), parameter :: myname_ = myname//"::write_" - integer, parameter :: NDIMS = 3 - -! When we create netCDF files, variables and dimensions, we get back -! an ID for each one. - character(len=4) :: cindx - integer :: ncid, dimids(NDIMS) - integer :: x_dimid, y_dimid, z_dimid - integer :: lon_varid, lat_varid, lev_varid - integer :: ii,jj,nl,nv,nn,nlat,nlon,nlev - integer :: mype_,root_ - integer, allocatable :: varid1d(:), varid2d(:), varid2dx(:), varidMLL(:) - logical :: verbose - -! This is the data array we will write. It will just be filled with -! a progression of integers for this example. - real(4), allocatable :: data_out(:,:,:) - -! Return code (status) - rc=0; mype_=0; root_=0 - verbose=.true. - if(present(myid).and.present(root) )then - if(myid/=root) verbose=.false. - mype_ = myid - root_ = root - endif - -! Set dims - nlat=bvars%nlat - nlon=bvars%nlon - nlev=bvars%nsig - -! Always check the return code of every netCDF function call. In -! this example program, wrapping netCDF calls with "call check()" -! makes sure that any return which is not equal to nf90_noerr (0) -! will print a netCDF error message and exit. - -! Create the netCDF file. The nf90_clobber parameter tells netCDF to -! overwrite this file, if it already exists. - call check_( nf90_create(fname, NF90_CLOBBER, ncid), rc, mype_, root_ ) - if(rc/=0) return - -! Define the dimensions. NetCDF will hand back an ID for each. - call check_( nf90_def_dim(ncid, "lon", nlon, x_dimid), rc, mype_, root_ ) - call check_( nf90_def_dim(ncid, "lat", nlat, y_dimid), rc, mype_, root_ ) - call check_( nf90_def_dim(ncid, "lev", nlev, z_dimid), rc, mype_, root_ ) - - call check_( nf90_def_var(ncid, "lon", NF90_REAL, x_dimid, lon_varid), rc, mype_, root_ ) - call check_( nf90_def_var(ncid, "lat", NF90_REAL, y_dimid, lat_varid), rc, mype_, root_ ) - call check_( nf90_def_var(ncid, "lev", NF90_REAL, z_dimid, lev_varid), rc, mype_, root_ ) - - call check_( nf90_put_att(ncid, lon_varid, "units", "degress"), rc, mype_, root_ ) - call check_( nf90_put_att(ncid, lat_varid, "units", "degress"), rc, mype_, root_ ) - call check_( nf90_put_att(ncid, lev_varid, "units", "hPa"), rc, mype_, root_ ) - -! The dimids array is used to pass the IDs of the dimensions of -! the variables. Note that in fortran arrays are stored in -! column-major format. - dimids = (/ x_dimid, y_dimid, z_dimid /) - -! Define variables. - allocate(varid1d(nv1d)) - do nv = 1, nv1d - call check_( nf90_def_var(ncid, trim(cvars1d(nv)), NF90_REAL, (/ y_dimid /), varid1d(nv)), rc, mype_, root_ ) - enddo - allocate(varid2d(nv2d)) - do nv = 1, nv2d - call check_( nf90_def_var(ncid, trim(cvars2d(nv)), NF90_REAL, (/ y_dimid, z_dimid /), varid2d(nv)), rc, mype_, root_ ) - enddo - allocate(varidMLL(nlev*nvmll)) - nn=0 - do nv = 1, nvmll - do nl = 1, nlev - nn=nn+1 - write(cindx,'(i4.4)') nl - call check_( nf90_def_var(ncid, trim(cvarsMLL(nv))//cindx, NF90_REAL, (/ y_dimid, z_dimid /), varidMLL(nn)), rc, & - mype_, root_ ) - enddo - enddo - allocate(varid2dx(nv2dx)) - do nv = 1, nv2dx - call check_( nf90_def_var(ncid, trim(cvars2dx(nv)), NF90_REAL, (/ x_dimid, y_dimid /), varid2dx(nv)), rc, mype_, root_ ) - enddo - -! End define mode. This tells netCDF we are done defining metadata. - call check_( nf90_enddef(ncid), rc, mype_, root_ ) - -! Write coordinate variables data - call check_( nf90_put_var(ncid, lon_varid, lons ), rc, mype_, root_ ) - call check_( nf90_put_var(ncid, lat_varid, lats ), rc, mype_, root_ ) - call check_( nf90_put_var(ncid, lev_varid, plevs), rc, mype_, root_ ) - -! Write data to file - allocate(data_out(1,nlat,1)) - do nv = 1, nv1d - if(trim(cvars1d(nv))=="ps" ) data_out(1,:,1) = bvars%psvar - if(trim(cvars1d(nv))=="hps" ) data_out(1,:,1) = bvars%pshln - call check_( nf90_put_var(ncid, varid1d(nv), data_out(1,:,1)), rc, mype_, root_) - enddo - deallocate(data_out) - allocate(data_out(1,nlat,nlev)) - do nv = 1, nv2d - if(trim(cvars2d(nv))=="sf" ) data_out(1,:,:) = bvars%sfvar - if(trim(cvars2d(nv))=="hsf") data_out(1,:,:) = bvars%sfhln - if(trim(cvars2d(nv))=="vsf") data_out(1,:,:) = bvars%sfvln -! - if(trim(cvars2d(nv))=="vp" ) data_out(1,:,:) = bvars%vpvar - if(trim(cvars2d(nv))=="hvp") data_out(1,:,:) = bvars%vphln - if(trim(cvars2d(nv))=="vvp") data_out(1,:,:) = bvars%vpvln -! - if(trim(cvars2d(nv))=="t" ) data_out(1,:,:) = bvars%tvar - if(trim(cvars2d(nv))=="ht" ) data_out(1,:,:) = bvars%thln - if(trim(cvars2d(nv))=="vt" ) data_out(1,:,:) = bvars%tvln -! - if(trim(cvars2d(nv))=="q" ) data_out(1,:,:) = bvars%qvar - if(trim(cvars2d(nv))=="hq" ) data_out(1,:,:) = bvars%qhln - if(trim(cvars2d(nv))=="vq" ) data_out(1,:,:) = bvars%qvln -! - if(trim(cvars2d(nv))=="qi" ) data_out(1,:,:) = bvars%qivar - if(trim(cvars2d(nv))=="hqi") data_out(1,:,:) = bvars%qihln - if(trim(cvars2d(nv))=="vqi") data_out(1,:,:) = bvars%qivln -! - if(trim(cvars2d(nv))=="ql" ) data_out(1,:,:) = bvars%qlvar - if(trim(cvars2d(nv))=="hql") data_out(1,:,:) = bvars%qlhln - if(trim(cvars2d(nv))=="vql") data_out(1,:,:) = bvars%qlvln -! - if(trim(cvars2d(nv))=="qr" ) data_out(1,:,:) = bvars%qrvar - if(trim(cvars2d(nv))=="hqr") data_out(1,:,:) = bvars%qrhln - if(trim(cvars2d(nv))=="vqr") data_out(1,:,:) = bvars%qrvln -! - if(trim(cvars2d(nv))=="nrh") data_out(1,:,:) = bvars%nrhvar - if(trim(cvars2d(nv))=="qs" ) data_out(1,:,:) = bvars%qsvar - if(trim(cvars2d(nv))=="hqs") data_out(1,:,:) = bvars%qshln - if(trim(cvars2d(nv))=="vqs") data_out(1,:,:) = bvars%qsvln -! - if(trim(cvars2d(nv))=="cw" ) data_out(1,:,:) = bvars%cvar - if(trim(cvars2d(nv))=="hcw") data_out(1,:,:) = bvars%chln - if(trim(cvars2d(nv))=="vcw") data_out(1,:,:) = bvars%cvln -! - if(trim(cvars2d(nv))=="oz" ) data_out(1,:,:) = bvars%ozvar - if(trim(cvars2d(nv))=="hoz") data_out(1,:,:) = bvars%ozhln - if(trim(cvars2d(nv))=="voz") data_out(1,:,:) = bvars%ozvln -! - if(trim(cvars2d(nv))=="pscon") data_out(1,:,:) = bvars%pscon - if(trim(cvars2d(nv))=="vpcon") data_out(1,:,:) = bvars%vpcon -! - call check_( nf90_put_var(ncid, varid2d(nv), data_out(1,:,:)), rc, mype_, root_ ) - enddo - -! Choose to write out NLATxNLEVxNLEV vars as to facilitate visualization - nn=0 - do nv = 1, nvmll - do nl = 1, nlev - nn = nn + 1 - write(cindx,'(i4.4)') nl - if(trim(cvarsMLL(nv))=="tcon") data_out(1,:,:) = bvars%tcon(:,:,nl) - call check_( nf90_put_var(ncid, varidMLL(nn), data_out(1,:,:)), rc, mype_, root_ ) - enddo - enddo - deallocate(data_out) - -! Write out lat/lon fields - allocate(data_out(nlon,nlat,1)) - do nv = 1, nv2dx - if(trim(cvars2dx(nv))=="sst" ) then - data_out(:,:,1) = transpose(bvars%varsst) - endif - if(trim(cvars2dx(nv))=="hsst" ) then - data_out(:,:,1) = transpose(bvars%corlsst) - endif - call check_( nf90_put_var(ncid, varid2dx(nv), data_out(:,:,1)), rc, mype_, root_ ) - enddo - deallocate(data_out) - -! Close file - call check_( nf90_close(ncid), rc, mype_, root_ ) - - deallocate(varidMLL) - deallocate(varid2d) - deallocate(varid1d) - - if(verbose) print *,"*** Finish writing file: ", trim(fname) - - return - -end subroutine write_berror_ - subroutine init_berror_vars_(vr,nlon,nlat,nsig) integer,intent(in) :: nlon,nlat,nsig From 849b3b19cc2160a7c32dfb3fe04c650cf78f3abc Mon Sep 17 00:00:00 2001 From: "russ.treadon" Date: Mon, 2 Dec 2024 16:56:55 +0000 Subject: [PATCH 6/9] additional code cleanup, enhance error checking (#808) --- src/gsi/m_berror_stats.f90 | 25 +++++++++++++++---------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/src/gsi/m_berror_stats.f90 b/src/gsi/m_berror_stats.f90 index 9a97a199ca..b4735889cd 100644 --- a/src/gsi/m_berror_stats.f90 +++ b/src/gsi/m_berror_stats.f90 @@ -705,17 +705,21 @@ subroutine nc_(myid) n = getindex(cvars2d,'sst') found2d(n)=.true. call nc_berror_getpointer (cvars2d(nv),bvars,ptr2d,ier) - if(ier==0) corsst=ptr2d + if (ier/=0) call die(myname_," in cw, failed to find corsst ", 99) + corsst=ptr2d call nc_berror_getpointer ('h'//cvars2d(nv),bvars,ptr2d,ier) - if(ier==0) hsst=ptr2d + if(ier/=0) call die(myname_," in cw, failed to find hsst ", 99) + hsst=ptr2d endif if (trim(cvars2d(nv))=='ps') then n = getindex(cvars2d,'ps') found2d(n)=.true. call nc_berror_getpointer (cvars2d(nv),bvars,ptr1d,ier) - if(ier==0) corp(:,n)=ptr1d + if(ier/=0) call die(myname_," in cw, failed to find corp ", 99) + corp(:,n)=ptr1d call nc_berror_getpointer ('h'//cvars2d(nv),bvars,ptr1d,ier) - if(ier==0) hwllp(:,n)=ptr1d + if(ier/=0) call die(myname_," in cw, failed to find hwllp ", 99) + hwllp(:,n)=ptr1d endif enddo do nv=1,size(cvars3d) @@ -725,9 +729,11 @@ subroutine nc_(myid) found3d(n)=.true. corz(:,:,n)=ptr2d call nc_berror_getpointer ('h'//trim(cvars3d(nv)),bvars,ptr2d,ier) - if(ier==0) hwll(:,:,n)=ptr2d + if(ier/=0) call die(myname_," in cw, failed to find hwll ", 99) + hwll(:,:,n)=ptr2d call nc_berror_getpointer ('v'//trim(cvars3d(nv)),bvars,ptr2d,ier) - if(ier==0) vz(:,:,n)=transpose(ptr2d) + if(ier/=0) call die(myname_," in cw, failed to find vz ", 99) + vz(:,:,n)=transpose(ptr2d) if (trim(cvars3d(nv))=='cw' .and. cwoption==2) then allocate(corq2(bvars%nlat,bvars%nsig)) call nc_berror_getpointer ('nrh',bvars,ptr2d,ier) @@ -750,9 +756,6 @@ subroutine nc_(myid) call nc_berror_getpointer ('nrh',bvars,ptr2d,ier) if (ier==0) then corq2=ptr2d -! corq2=max(0.0_r_kind,corq2) ! hack 1 -! corq2=min(1.0_r_kind,2*corq2) ! hack 2 -! print *, 'DEBUG (berr): ', minval(corq2),maxval(corq2) do k=1,bvars%nsig do i=1,bvars%nlat corq2x=corq2(i,k) @@ -766,7 +769,9 @@ subroutine nc_(myid) deallocate(corq2) endif cycle - endif + else + call die(myname_," in cw, failed to find cvars3d bvars ", 99) + endif if (trim(cvars3d(nv))=='q') then n = getindex(cvars3d,'q') found3d(n)=.true. From 712551f8b18c6d5aeb861868f9d9b85a7aaf321e Mon Sep 17 00:00:00 2001 From: "russ.treadon" Date: Mon, 2 Dec 2024 18:30:47 +0000 Subject: [PATCH 7/9] remove erroneously placed error trapping line (#808) --- src/gsi/m_berror_stats.f90 | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/gsi/m_berror_stats.f90 b/src/gsi/m_berror_stats.f90 index b4735889cd..b72fbdb3eb 100644 --- a/src/gsi/m_berror_stats.f90 +++ b/src/gsi/m_berror_stats.f90 @@ -769,9 +769,7 @@ subroutine nc_(myid) deallocate(corq2) endif cycle - else - call die(myname_," in cw, failed to find cvars3d bvars ", 99) - endif + endif if (trim(cvars3d(nv))=='q') then n = getindex(cvars3d,'q') found3d(n)=.true. From 0e07d71fdd846c5590e44285f14ad4b9f23f94c9 Mon Sep 17 00:00:00 2001 From: "russ.treadon" Date: Mon, 2 Dec 2024 18:39:03 +0000 Subject: [PATCH 8/9] update GSI_BINARY_SOURCE_DIR in machine specific modulefiles (#808) --- modulefiles/gsi_acorn.intel.lua | 2 +- modulefiles/gsi_gaeac5.intel.lua | 2 +- modulefiles/gsi_gaeac6.intel.lua | 2 +- modulefiles/gsi_hera.gnu.lua | 2 +- modulefiles/gsi_hera.intel.lua | 2 +- modulefiles/gsi_hercules.intel.lua | 2 +- modulefiles/gsi_jet.intel.lua | 2 +- modulefiles/gsi_noaacloud.intel.lua | 2 +- modulefiles/gsi_orion.intel.lua | 2 +- modulefiles/gsi_s4.intel.lua | 2 +- modulefiles/gsi_wcoss2.intel.lua | 2 +- 11 files changed, 11 insertions(+), 11 deletions(-) diff --git a/modulefiles/gsi_acorn.intel.lua b/modulefiles/gsi_acorn.intel.lua index 03928be822..cde59c7ea4 100644 --- a/modulefiles/gsi_acorn.intel.lua +++ b/modulefiles/gsi_acorn.intel.lua @@ -49,6 +49,6 @@ load(pathJoin("ncdiag",ncdiag_ver)) append_path("MODULEPATH", "/lfs/h1/emc/nceplibs/noscrub/hpc-stack/libs/hpc-stack/modulefiles/compiler/intel/19.1.3.304") load(pathJoin("crtm", crtm_ver)) -pushenv("GSI_BINARY_SOURCE_DIR", "/lfs/h2/emc/global/noscrub/emc.global/FIX/fix/gsi/20240208") +pushenv("GSI_BINARY_SOURCE_DIR", "/lfs/h2/emc/global/noscrub/emc.global/FIX/fix/gsi/20241022") whatis("Description: GSI environment on WCOSS2 Acorn") diff --git a/modulefiles/gsi_gaeac5.intel.lua b/modulefiles/gsi_gaeac5.intel.lua index de259f1741..f660a742b1 100644 --- a/modulefiles/gsi_gaeac5.intel.lua +++ b/modulefiles/gsi_gaeac5.intel.lua @@ -17,7 +17,7 @@ load(pathJoin("cmake", cmake_ver)) load("gsi_common") load(pathJoin("prod_util", prod_util_ver)) -pushenv("GSI_BINARY_SOURCE_DIR", "/gpfs/f5/ufs-ard/world-shared/GSI_data/fix/gsi/20240208") +pushenv("GSI_BINARY_SOURCE_DIR", "/gpfs/f5/ufs-ard/world-shared/GSI_data/fix/gsi/20241022") setenv("CC","cc") setenv("FC","ftn") diff --git a/modulefiles/gsi_gaeac6.intel.lua b/modulefiles/gsi_gaeac6.intel.lua index 883bd02a72..bbf1a49181 100644 --- a/modulefiles/gsi_gaeac6.intel.lua +++ b/modulefiles/gsi_gaeac6.intel.lua @@ -18,7 +18,7 @@ load(pathJoin("cmake", cmake_ver)) load("gsi_common") load(pathJoin("prod_util", prod_util_ver)) -pushenv("GSI_BINARY_SOURCE_DIR", "/gpfs/f6/bil-fire8/world-shared/GSI_data/fix/gsi/20240208") +pushenv("GSI_BINARY_SOURCE_DIR", "/gpfs/f6/bil-fire8/world-shared/GSI_data/fix/gsi/20241022") setenv("CC","cc") setenv("FC","ftn") diff --git a/modulefiles/gsi_hera.gnu.lua b/modulefiles/gsi_hera.gnu.lua index eab352553f..815005855d 100644 --- a/modulefiles/gsi_hera.gnu.lua +++ b/modulefiles/gsi_hera.gnu.lua @@ -22,6 +22,6 @@ load("gsi_common") load(pathJoin("prod_util", prod_util_ver)) load(pathJoin("openblas", openblas_ver)) -pushenv("GSI_BINARY_SOURCE_DIR", "/scratch1/NCEPDEV/global/glopara/fix/gsi/20240208") +pushenv("GSI_BINARY_SOURCE_DIR", "/scratch1/NCEPDEV/global/glopara/fix/gsi/20241022") whatis("Description: GSI environment on Hera with GNU Compilers") diff --git a/modulefiles/gsi_hera.intel.lua b/modulefiles/gsi_hera.intel.lua index d21b9195c3..0bf49f4224 100644 --- a/modulefiles/gsi_hera.intel.lua +++ b/modulefiles/gsi_hera.intel.lua @@ -20,6 +20,6 @@ load(pathJoin("prod_util", prod_util_ver)) pushenv("CFLAGS", "-xHOST") pushenv("FFLAGS", "-xHOST") -pushenv("GSI_BINARY_SOURCE_DIR", "/scratch1/NCEPDEV/global/glopara/fix/gsi/20240208") +pushenv("GSI_BINARY_SOURCE_DIR", "/scratch1/NCEPDEV/global/glopara/fix/gsi/20241022") whatis("Description: GSI environment on Hera with Intel Compilers") diff --git a/modulefiles/gsi_hercules.intel.lua b/modulefiles/gsi_hercules.intel.lua index 66ec9b03e1..4426170fc9 100644 --- a/modulefiles/gsi_hercules.intel.lua +++ b/modulefiles/gsi_hercules.intel.lua @@ -21,6 +21,6 @@ load("intel-oneapi-mkl/2022.2.1") pushenv("CFLAGS", "-xHOST") pushenv("FFLAGS", "-xHOST") -pushenv("GSI_BINARY_SOURCE_DIR", "/work/noaa/global/glopara/fix/gsi/20240208") +pushenv("GSI_BINARY_SOURCE_DIR", "/work/noaa/global/glopara/fix/gsi/20241022") whatis("Description: GSI environment on Hercules with Intel Compilers") diff --git a/modulefiles/gsi_jet.intel.lua b/modulefiles/gsi_jet.intel.lua index b210e6efa2..ea9dd38209 100644 --- a/modulefiles/gsi_jet.intel.lua +++ b/modulefiles/gsi_jet.intel.lua @@ -20,6 +20,6 @@ load(pathJoin("prod_util", prod_util_ver)) pushenv("CFLAGS", "-axSSE4.2,AVX,CORE-AVX2") pushenv("FFLAGS", "-axSSE4.2,AVX,CORE-AVX2") -pushenv("GSI_BINARY_SOURCE_DIR", "/lfs5/HFIP/hfv3gfs/glopara/FIX/fix/gsi/20240208") +pushenv("GSI_BINARY_SOURCE_DIR", "/lfs5/HFIP/hfv3gfs/glopara/FIX/fix/gsi/20241022") whatis("Description: GSI environment on Jet with Intel Compilers") diff --git a/modulefiles/gsi_noaacloud.intel.lua b/modulefiles/gsi_noaacloud.intel.lua index e2e019628e..05226bed89 100644 --- a/modulefiles/gsi_noaacloud.intel.lua +++ b/modulefiles/gsi_noaacloud.intel.lua @@ -20,6 +20,6 @@ load(pathJoin("prod_util", prod_util_ver)) pushenv("CFLAGS", "-xHOST") pushenv("FFLAGS", "-xHOST") -pushenv("GSI_BINARY_SOURCE_DIR", "/contrib/Wei.Huang/data/hack-orion/fix/gsi/20240208") +pushenv("GSI_BINARY_SOURCE_DIR", "/contrib/Wei.Huang/data/hack-orion/fix/gsi/20241022") whatis("Description: GSI environment on NOAA Cloud with Intel Compilers") diff --git a/modulefiles/gsi_orion.intel.lua b/modulefiles/gsi_orion.intel.lua index d05bda5b2e..05cf13b7fc 100644 --- a/modulefiles/gsi_orion.intel.lua +++ b/modulefiles/gsi_orion.intel.lua @@ -21,6 +21,6 @@ load("intel-oneapi-mkl/2022.2.1") pushenv("CFLAGS", "-xHOST") pushenv("FFLAGS", "-xHOST") -pushenv("GSI_BINARY_SOURCE_DIR", "/work/noaa/global/glopara/fix/gsi/20240208") +pushenv("GSI_BINARY_SOURCE_DIR", "/work/noaa/global/glopara/fix/gsi/20241022") whatis("Description: GSI environment on Orion with Intel Compilers") diff --git a/modulefiles/gsi_s4.intel.lua b/modulefiles/gsi_s4.intel.lua index 04945eef3e..069b6f90c4 100644 --- a/modulefiles/gsi_s4.intel.lua +++ b/modulefiles/gsi_s4.intel.lua @@ -20,6 +20,6 @@ load(pathJoin("prod_util", prod_util_ver)) pushenv("CFLAGS", "-march=ivybridge") pushenv("FFLAGS", "-march=ivybridge") -pushenv("GSI_BINARY_SOURCE_DIR", "/data/prod/glopara/fix/gsi/20240208") +pushenv("GSI_BINARY_SOURCE_DIR", "/data/prod/glopara/fix/gsi/20241022") whatis("Description: GSI environment on S4 with Intel Compilers") diff --git a/modulefiles/gsi_wcoss2.intel.lua b/modulefiles/gsi_wcoss2.intel.lua index c3bfd1156c..e4f367ab95 100644 --- a/modulefiles/gsi_wcoss2.intel.lua +++ b/modulefiles/gsi_wcoss2.intel.lua @@ -46,6 +46,6 @@ load(pathJoin("ncio", ncio_ver)) load(pathJoin("crtm", crtm_ver)) load(pathJoin("ncdiag",ncdiag_ver)) -pushenv("GSI_BINARY_SOURCE_DIR", "/lfs/h2/emc/global/noscrub/emc.global/FIX/fix/gsi/20240208") +pushenv("GSI_BINARY_SOURCE_DIR", "/lfs/h2/emc/global/noscrub/emc.global/FIX/fix/gsi/20241022") whatis("Description: GSI environment on WCOSS2") From 63df6ff4230305b497525c07de55da8502e4112c Mon Sep 17 00:00:00 2001 From: "russ.treadon" Date: Fri, 6 Dec 2024 19:58:36 +0000 Subject: [PATCH 9/9] clean up printout, simplify global_berror file format logic (#808) --- src/gsi/glbsoi.f90 | 3 --- src/gsi/m_berror_stats.f90 | 35 +++++++++++++++++++++++------------ src/gsi/ncepnems_io.f90 | 24 +++++++----------------- 3 files changed, 30 insertions(+), 32 deletions(-) diff --git a/src/gsi/glbsoi.f90 b/src/gsi/glbsoi.f90 index a210ba3258..d9dbee7b3b 100644 --- a/src/gsi/glbsoi.f90 +++ b/src/gsi/glbsoi.f90 @@ -161,7 +161,6 @@ subroutine glbsoi use m_prad, only: prad_updatePredx ! was -- prad_bias() use m_obsdiags, only: obsdiags_write use gsi_io,only: verbose - use m_berror_stats,only: inquire_berror implicit none @@ -257,8 +256,6 @@ subroutine glbsoi end if end if else - lunit=22 - call inquire_berror(lunit,mype) call create_balance_vars if(anisotropic) then call create_anberror_vars(mype) diff --git a/src/gsi/m_berror_stats.f90 b/src/gsi/m_berror_stats.f90 index b72fbdb3eb..0ed00bb435 100644 --- a/src/gsi/m_berror_stats.f90 +++ b/src/gsi/m_berror_stats.f90 @@ -55,6 +55,7 @@ module m_berror_stats ! reconfigurable parameters, via NAMELIST/setup/ public :: usenewgfsberror public :: berror_stats,inquire_berror ! reconfigurable filename + public :: bin_berror ! interfaces to file berror_stats. public :: berror_get_dims ! get dimensions, jfunc::createj_func() @@ -389,8 +390,14 @@ subroutine nc_(myid) call die(myname_," fut2ps not available in this form "//trim(berror_stats), 99) endif call nc_berror_read (berror_stats,bvars,ier, myid=myid,root=0) - if (nlat/=bvars%nlat .or. nsig/=bvars%nsig ) then - call die(myname_," inconsistent dims in "//trim(berror_stats), 99) + if ( mype == 0 ) then + if (nlat/=bvars%nlat .or. nsig/=bvars%nsig ) then + call die(myname_," inconsistent dims in "//trim(berror_stats), 99) + endif + write(6,*) myname_,'(PREBAL): get balance variables', & + '"',trim(berror_stats),'". ', & + 'mype,nsigstat,nlatstat =', & + mype,bvars%nsig,bvars%nlat endif agvin = bvars%tcon bvin = bvars%vpcon @@ -694,8 +701,14 @@ subroutine nc_(myid) real(r_single), pointer :: ptr2d(:,:) integer :: nv call nc_berror_read (berror_stats,bvars,ier, myid=myid,root=0) - if (nlat/=bvars%nlat .or. nlon/=bvars%nlon .or. nsig/=bvars%nsig ) then - call die(myname_," inconsistent dims in "//trim(berror_stats), 99) + if ( mype==0 ) then + if (nlat/=bvars%nlat .or. nlon/=bvars%nlon .or. nsig/=bvars%nsig ) then + call die(myname_," inconsistent dims in "//trim(berror_stats), 99) + endif + write(6,*) myname_,'(PREWGT): read error amplitudes ', & + '"',trim(berror_stats),'". ', & + 'mype,nsigstat,nlatstat =', & + mype,bvars%nsig,bvars%nlat endif isig=bvars%nsig @@ -1049,7 +1062,7 @@ subroutine setcorchem_(cname,corchem,rc) rc=0 ! sanity check - if ( mype==0 ) write(6,*) myname_,'(PREWGT): enter routine' + if ( mype==0 ) write(6,*) myname_,'(PREWGT): mype = ',mype ! Get information for how to use CO2 iptr=-1 @@ -1167,15 +1180,13 @@ subroutine sethwllchem_(hwll, mype) real(r_kind) :: fact real(r_kind) :: s2u - if (mype == 0) then - write(6,*) myname_, '(PREWGT): mype = ', mype - end if + if (mype == 0) write(6,*) myname_, '(PREWGT): mype = ',mype s2u = (two*pi*rearth_equator)/nlon do k = 1,nnnn1o k1 = levs_id(k) if (k1 > 0) then - if (mype == 0) write(6,*) myname_, '(PREWGT): mype = ', mype, k1 +! if (mype == 0) write(6,*) myname_, '(PREWGT): mype = ', mype, k1 ! make everything constant ! fact = real(k1,r_kind)**2._r_kind fact = 1._r_kind @@ -1184,9 +1195,9 @@ subroutine sethwllchem_(hwll, mype) end if end do - if (mype == 0) then - write(6,*) myname_, '(PREWGT): mype = ', mype, 'finish sethwllchem_' - end if +! if (mype == 0) then +! write(6,*) myname_, '(PREWGT): mype = ', mype, 'finish sethwllchem_' +! end if end subroutine sethwllchem_ end module m_berror_stats diff --git a/src/gsi/ncepnems_io.f90 b/src/gsi/ncepnems_io.f90 index 5a8b80d71b..f5ae00ecaf 100755 --- a/src/gsi/ncepnems_io.f90 +++ b/src/gsi/ncepnems_io.f90 @@ -1813,6 +1813,7 @@ subroutine read_hsst_(hsst,cwoption) use m_nc_berror, only: nc_berror_vars use m_nc_berror, only: nc_berror_getpointer use m_nc_berror, only: nc_berror_vars_final + use m_berror_stats, only: bin_berror use mpimod, only: mype implicit none @@ -1838,23 +1839,11 @@ subroutine read_hsst_(hsst,cwoption) character(len=256) :: berror_stats = "berror_stats" ! filename - logical :: ncio - type(Dataset) :: dset - type(Dimension) :: londim,latdim,levdim integer(i_kind) :: nv,n type(nc_berror_vars) bvars real(r_single), pointer :: ptr2d(:,:) - dset = open_dataset(berror_stats,errcode=ier) - if (ier==0) then - ! this is a netcdf file - ncio = .true. - else - ! this is a binary file - ncio = .false. - endif - - if (ncio) then + if (.not.bin_berror) then call nc_berror_read (berror_stats,bvars,ier, myid=mype,root=0) if (nlat/=bvars%nlat .or. nlon/=bvars%nlon .or. nsig/=bvars%nsig ) then call die(myname_," inconsistent dims in "//trim(berror_stats), 99) @@ -1864,10 +1853,11 @@ subroutine read_hsst_(hsst,cwoption) ! RTodling: the following is bad since it wires all naming conventions ... to be revised do nv=1,size(cvars2d) if (trim(cvars2d(nv))=='sst') then -!! n = getindex(cvars2d,'sst') -!! found2d(n)=.true. -!! call nc_berror_getpointer (cvars2d(nv),bvars,ptr2d,ier) -!! if(ier==0) corsst=ptr2d +! Do not need corsst in this routine so comment out code. Only need hsst +! n = getindex(cvars2d,'sst') +! found2d(n)=.true. +! call nc_berror_getpointer (cvars2d(nv),bvars,ptr2d,ier) +! if(ier==0) corsst=ptr2d call nc_berror_getpointer ('h'//cvars2d(nv),bvars,ptr2d,ier) if(ier==0) hsst=ptr2d endif