From 6fd4bd568d69f2e80c62f8db205fae2aadf5d506 Mon Sep 17 00:00:00 2001 From: RussTreadon-NOAA <26926959+RussTreadon-NOAA@users.noreply.github.com> Date: Tue, 10 Dec 2024 15:42:41 -0500 Subject: [PATCH 1/2] Add option to use netcdf global_berror files (#810) --- fix | 2 +- 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 +- src/gsi/glbsoi.f90 | 3 - src/gsi/gsi_files.cmake | 1 + src/gsi/m_berror_stats.f90 | 509 +++++++++++++++----- src/gsi/m_nc_berror.f90 | 699 ++++++++++++++++++++++++++++ src/gsi/ncepnems_io.f90 | 165 ++++--- 17 files changed, 1213 insertions(+), 188 deletions(-) create mode 100644 src/gsi/m_nc_berror.f90 diff --git a/fix b/fix index 2a86d5b1a0..1567130f4e 160000 --- a/fix +++ b/fix @@ -1 +1 @@ -Subproject commit 2a86d5b1a09a8908618d1c8a376c68e8d9b3d02c +Subproject commit 1567130f4eba20b301903e30114e72a734fd9a15 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") 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/gsi_files.cmake b/src/gsi/gsi_files.cmake index 04e7926300..5c5c4ec5a5 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 diff --git a/src/gsi/m_berror_stats.f90 b/src/gsi/m_berror_stats.f90 index 088a7619fe..0ed00bb435 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,28 @@ 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 + public :: bin_berror + ! 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 +81,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 +89,11 @@ 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. contains @@ -100,6 +109,7 @@ module m_berror_stats subroutine get_dims(msig,mlat,mlon,lunit) + use mpimod, only: mype implicit none integer(i_kind) ,intent( out) :: msig ! dimension of levels @@ -114,7 +124,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 +157,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 +188,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 +244,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 +325,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 +349,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 +381,31 @@ 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 ( 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 + wgvin = bvars%pscon + pputin=zero + call nc_berror_vars_final(bvars) + end subroutine nc_ + end subroutine read_bal !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -313,7 +420,7 @@ end subroutine read_bal subroutine read_wgt(corz,corp,hwll,hwllp,vz,corsst,hsst,varq,qoption,varcw,cwoption,mype,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,10 +438,12 @@ 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 ! !REVISION HISTORY: @@ -363,19 +472,102 @@ 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 + call bin_() + else + 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 +578,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 +598,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 @@ -468,6 +656,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 +687,116 @@ 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 + end subroutine bin_ - 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 - - ! 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 ( 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 + +! 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) 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) 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) 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) call die(myname_," in cw, failed to find hwllp ", 99) + 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) 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) 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) + 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 + 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 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -620,7 +852,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) @@ -694,7 +926,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 @@ -719,13 +951,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 - 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) @@ -737,7 +969,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_ @@ -754,7 +986,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 +1020,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 +1041,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 ___________________________________________________________________ @@ -892,7 +1126,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 +1144,60 @@ 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) 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 +! 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..0f3d5b329b --- /dev/null +++ b/src/gsi/m_nc_berror.f90 @@ -0,0 +1,699 @@ +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_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_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 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 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..f5ae00ecaf 100755 --- a/src/gsi/ncepnems_io.f90 +++ b/src/gsi/ncepnems_io.f90 @@ -1805,8 +1805,17 @@ 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 m_berror_stats, only: bin_berror + use mpimod, only: mype + implicit none integer(i_kind) , intent(in ) :: cwoption @@ -1830,72 +1839,98 @@ 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' + integer(i_kind) :: nv,n + type(nc_berror_vars) bvars + real(r_single), pointer :: ptr2d(:,:) - write(6,*)'PREBAL: berror nsigstat,nlatstat=', nsigstat,nlatstat, & - ' -vs- GSI nsig,nlat=',nsig,nlat - call stop2(101) - 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 (.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) + 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 +! 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 - 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 092125196a49e856ad760730edac01e398a84d8f Mon Sep 17 00:00:00 2001 From: jianjunj Date: Thu, 12 Dec 2024 14:00:37 -0500 Subject: [PATCH 2/2] Correct scattering index and enable error inflation in ATMS data, (#812) --- src/gsi/clw_mod.f90 | 16 +++++++++++++--- src/gsi/qcmod.f90 | 2 +- 2 files changed, 14 insertions(+), 4 deletions(-) diff --git a/src/gsi/clw_mod.f90 b/src/gsi/clw_mod.f90 index 49387f05eb..9d16ad64e1 100644 --- a/src/gsi/clw_mod.f90 +++ b/src/gsi/clw_mod.f90 @@ -1961,6 +1961,7 @@ subroutine ret_amsua(tb_obs,nchanl,tsavg5,zasat,clwp_amsua,ierrret,scat) ! Declare local variables real(r_kind) :: d0, d1, d2, coszat ! real(r_kind) :: c0, c1, c2 + real(r_kind) :: tb890 = zero coszat=cos(zasat) @@ -1979,9 +1980,18 @@ subroutine ret_amsua(tb_obs,nchanl,tsavg5,zasat,clwp_amsua,ierrret,scat) endif if (present(scat)) then - scat=-113.2_r_kind+(2.41_r_kind-0.0049_r_kind*tb_obs(1))*tb_obs(1) & - +0.454_r_kind*tb_obs(2)-tb_obs(15) - scat=max(zero,scat) + if (nchanl == 15) then +! AMSU-A + tb890 = tb_obs(15) + else if (nchanl == 22) then +! ATMS + tb890 = tb_obs(16) + endif + if (tb890 > zero) then + scat=-113.2_r_kind+(2.41_r_kind-0.0049_r_kind*tb_obs(1))*tb_obs(1) & + +0.454_r_kind*tb_obs(2)-tb890 + endif + scat=max(zero,scat) end if end subroutine ret_amsua diff --git a/src/gsi/qcmod.f90 b/src/gsi/qcmod.f90 index 86207a49fb..abd75108c3 100644 --- a/src/gsi/qcmod.f90 +++ b/src/gsi/qcmod.f90 @@ -3670,7 +3670,7 @@ subroutine qc_amsua(nchanl,is,ndat,nsig,npred,sea,land,ice,snow,mixed,luse, & ework = ework+min(0.002_r_kind*sfc_speed**2*error0(i), 0.5_r_kind*error0(i)) clwtmp=min(abs(clwp_amsua-clw_guess_retrieval), one) ework = ework+min(13.0_r_kind*clwtmp*error0(i), 3.5_r_kind*error0(i)) - if (scatp>9.0_r_kind .and. nchanl==15) then + if (scatp>9.0_r_kind) then ework = ework+min(1.5_r_kind*(scatp-9.0_r_kind)*error0(i), 2.5_r_kind*error0(i)) end if ework=ework**2