From a827064c4ecdcb9eb6c60a45c1a578b8ebd29522 Mon Sep 17 00:00:00 2001 From: KariA-Spire <104657787+KariA-Spire@users.noreply.github.com> Date: Tue, 9 Jul 2024 08:57:22 -0600 Subject: [PATCH 01/11] Gnssrwnd1.0 (#747) New GNSS-R ocean wind speed observation operator and assimilation capability in GSI. --- src/gsi/convinfo.f90 | 16 +- src/gsi/gsi_files.cmake | 6 + src/gsi/gsi_gnssrspdOper.F90 | 161 +++++ src/gsi/gsi_obOperTypeManager.F90 | 12 +- src/gsi/gsimain.f90 | 4 +- src/gsi/intgnssrspd.f90 | 239 ++++++++ src/gsi/intjo.f90 | 5 +- src/gsi/m_gnssrspdNode.F90 | 257 ++++++++ src/gsi/m_obsNodeTypeManager.F90 | 15 +- src/gsi/m_rhs.F90 | 4 +- src/gsi/obsmod.F90 | 11 +- src/gsi/read_gnssrspd.f90 | 463 +++++++++++++++ src/gsi/read_obs.F90 | 22 +- src/gsi/setupgnssrspd.f90 | 940 ++++++++++++++++++++++++++++++ src/gsi/setuprhsall.f90 | 5 +- src/gsi/statsconv.f90 | 55 +- src/gsi/stpgnssrspd.f90 | 175 ++++++ 17 files changed, 2363 insertions(+), 27 deletions(-) create mode 100644 src/gsi/gsi_gnssrspdOper.F90 create mode 100644 src/gsi/intgnssrspd.f90 create mode 100644 src/gsi/m_gnssrspdNode.F90 create mode 100644 src/gsi/read_gnssrspd.f90 create mode 100644 src/gsi/setupgnssrspd.f90 create mode 100644 src/gsi/stpgnssrspd.f90 diff --git a/src/gsi/convinfo.f90 b/src/gsi/convinfo.f90 index f6ba29c2e1..83c84d3f38 100755 --- a/src/gsi/convinfo.f90 +++ b/src/gsi/convinfo.f90 @@ -155,6 +155,8 @@ subroutine convinfo_read ! 2009-01-22 todling - protect against non-initialized destroy call ! 2010-05-29 todling - interface consistent w/ similar routines ! 2014-07-10 carley - add check to bypass blank lines in convinfo file +! 2023-05-09 s.vetra-carvalho - extended iotype to be 8 characters long to +! allow for more descriptive observations names ! ! input argument list: ! mype - mpi task id @@ -171,7 +173,7 @@ subroutine convinfo_read implicit none character(len=1)cflg - character(len=7) iotype + character(len=8) iotype character(len=140) crecord integer(i_kind) lunin,i,nc,ier,istat integer(i_kind) nlines @@ -188,13 +190,13 @@ subroutine convinfo_read nlines=0 read1: do cflg=' ' - iotype=' ' + iotype=' ' read(lunin,1030,iostat=istat,end=1130)cflg,iotype,crecord -1030 format(a1,a7,2x,a140) +1030 format(a1,a8,1x,a140) if (istat /= 0) exit nlines=nlines+1 if(cflg == '!')cycle - if (cflg==' '.and.iotype==' ') then + if (cflg==' '.and.iotype==' ') then if(print_verbose)write(6,*) 'Encountered a blank line in convinfo file at line number: ',nlines,' skipping!' cycle end if @@ -242,9 +244,9 @@ subroutine convinfo_read do i=1,nlines cflg=' ' - iotype=' ' + iotype=' ' read(lunin,1030)cflg,iotype,crecord - if (cflg==' '.and.iotype==' ') then + if (cflg==' '.and.iotype==' ') then if(print_verbose)write(6,*) 'Encountered a blank line in convinfo file at line number: ',i,' skipping!' cycle end if @@ -282,7 +284,7 @@ subroutine convinfo_read if(print_verbose .and. mype == 0)write(6,1031)ioctype(nc),ictype(nc),icsubtype(nc),icuse(nc),ctwind(nc),ncnumgrp(nc), & ncgroup(nc),ncmiter(nc),cgross(nc),cermax(nc),cermin(nc),cvar_b(nc),cvar_pg(nc), & ithin_conv(nc),rmesh_conv(nc),pmesh_conv(nc),idum,pmot_conv(nc),ptime_conv(nc),index_sub(nc),ibeta(nc),ikapa(nc) -1031 format('READ_CONVINFO: ',a7,1x,i3,1x,i4,1x,i2,1x,g13.6,1x,3(I3,1x),5g13.6,i5,2g13.6,i5,2g13.6,3i5) +1031 format('READ_CONVINFO: ',a8,1x,i3,1x,i4,1x,i2,1x,g13.6,1x,3(I3,1x),5g13.6,i5,2g13.6,i5,2g13.6,3i5) enddo close(lunin) diff --git a/src/gsi/gsi_files.cmake b/src/gsi/gsi_files.cmake index 95d885e2ee..97dcc3101f 100644 --- a/src/gsi/gsi_files.cmake +++ b/src/gsi/gsi_files.cmake @@ -218,6 +218,7 @@ gsi_dbzOper.F90 gsi_dwOper.F90 gsi_enscouplermod.f90 gsi_fedOper.F90 +gsi_gnssrspdOper.F90 gsi_gpsbendOper.F90 gsi_gpsrefOper.F90 gsi_gustOper.F90 @@ -276,6 +277,7 @@ intco.f90 intdbz.f90 intfed.f90 intdw.f90 +intgnssrspd.f90 intgps.f90 intgust.f90 inthowv.f90 @@ -342,6 +344,7 @@ m_dwNode.F90 m_extOzone.F90 m_fedNode.F90 m_find.f90 +m_gnssrspdNode.F90 m_gpsNode.F90 m_gpsrhs.F90 m_gsiBiases.f90 @@ -486,6 +489,7 @@ read_files.f90 read_fl_hdob.f90 read_gfs_ozone_for_regional.f90 read_gmi.f90 +read_gnssrspd.f90 read_goesglm.f90 read_goesimg.f90 read_goesimgr_skycover.f90 @@ -537,6 +541,7 @@ setupdbz.f90 setupdbz_lib.f90 setupdw.f90 setupfed.f90 +setupgnssrspd.f90 setupgust.f90 setuphowv.f90 setuplag.f90 @@ -597,6 +602,7 @@ stpco.f90 stpdbz.f90 stpfed.f90 stpdw.f90 +stpgnssrspd.f90 stpgps.f90 stpgust.f90 stphowv.f90 diff --git a/src/gsi/gsi_gnssrspdOper.F90 b/src/gsi/gsi_gnssrspdOper.F90 new file mode 100644 index 0000000000..92f0d0a2a6 --- /dev/null +++ b/src/gsi/gsi_gnssrspdOper.F90 @@ -0,0 +1,161 @@ +module gsi_gnssrspdOper +!$$$ subprogram documentation block +! . . . . +! subprogram: module gsi_gnssrspdOper +! (based on gsi_*Oper.F90 routines by j guo ) +! prgmmr: K. Apodaca +! org: Spire Global, Inc. +! date: 2023-04-28 +! +! abstract: an obOper extension for gnssrspdNode type +! +! program history log: +! 2023-04-21 k apodaca - initial version +! input argument list: see Fortran 90 style document below +! +! output argument list: see Fortran 90 style document below +! +! attributes: +! language: Fortran 90 and/or above +! machine: +! +!$$$ end subprogram documentation block + +! module interface: + + use gsi_obOper, only: obOper + use m_gnssrspdNode , only: gnssrspdNode + implicit none + public:: gnssrspdOper ! data stracture + + type,extends(obOper):: gnssrspdOper + contains + procedure,nopass:: mytype + procedure,nopass:: nodeMold + procedure:: setup_ + procedure:: intjo1_ + procedure:: stpjo1_ + end type gnssrspdOper + +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + character(len=*),parameter :: myname='gsi_gnssrspdOper' + type(gnssrspdNode),save,target:: myNodeMold_ + +contains + function mytype(nodetype) + implicit none + character(len=:),allocatable:: mytype + logical,optional, intent(in):: nodetype + mytype="[gnssrspdOper]" + if(present(nodetype)) then + if(nodetype) mytype=myNodeMold_%mytype() + endif + end function mytype + + function nodeMold() + !> %nodeMold() returns a mold of its corresponding obsNode + use m_obsNode, only: obsNode + implicit none + class(obsNode),pointer:: nodeMold + nodeMold => myNodeMold_ + end function nodeMold + + subroutine setup_(self, lunin, mype, is, nobs, init_pass,last_pass) + use gnssrspd_setup, only: setup + use kinds, only: i_kind + use gsi_obOper, only: len_obstype + use gsi_obOper, only: len_isis + + use m_rhs , only: awork => rhs_awork + use m_rhs , only: bwork => rhs_bwork + use m_rhs , only: iwork => i_gnssrspd + + use obsmod , only: write_diag + use convinfo, only: diag_conv + use jfunc , only: jiter + + use mpeu_util, only: die + implicit none + class(gnssrspdOper ), intent(inout):: self + integer(i_kind), intent(in):: lunin + integer(i_kind), intent(in):: mype + integer(i_kind), intent(in):: is + integer(i_kind), intent(in):: nobs + logical , intent(in):: init_pass ! supporting multi-pass setup() + logical , intent(in):: last_pass ! with incremental backgrounds. + + !---------------------------------------- + character(len=*),parameter:: myname_=myname//"::setup_" + + character(len=len_obstype):: obstype + character(len=len_isis ):: isis + integer(i_kind):: nreal,nchanl,ier,nele + logical:: diagsave + + if(nobs == 0) return + + read(lunin,iostat=ier) obstype,isis,nreal,nchanl + if(ier/=0) call die(myname_,'read(obstype,...), iostat =',ier) + nele = nreal+nchanl + + diagsave = write_diag(jiter) .and. diag_conv + + call setup(self%obsLL(:), self%odiagLL(:), & + lunin,mype,bwork,awork(:,iwork),nele,nobs,is,diagsave) + + end subroutine setup_ + + subroutine intjo1_(self, ibin, rval,sval, qpred,sbias) + use intgnssrspdmod, only: intjo => intgnssrspd + use gsi_bundlemod , only: gsi_bundle + use bias_predictors, only: predictors + use m_obsNode , only: obsNode + use m_obsLList, only: obsLList_headNode + use kinds , only: i_kind, r_quad + implicit none + class(gnssrspdOper ),intent(in ):: self + integer(i_kind ),intent(in ):: ibin + type(gsi_bundle),intent(inout):: rval ! (ibin) + type(gsi_bundle),intent(in ):: sval ! (ibin) + real(r_quad ),target,dimension(:),intent(inout):: qpred ! (ibin) + type(predictors),target, intent(in ):: sbias + + !---------------------------------------- + character(len=*),parameter:: myname_=myname//"::intjo1_" + class(obsNode),pointer:: headNode + + headNode => obsLList_headNode(self%obsLL(ibin)) + call intjo(headNode, rval,sval) + headNode => null() + + end subroutine intjo1_ + + subroutine stpjo1_(self, ibin, dval,xval,pbcjo,sges,nstep,dbias,xbias) + use stpgnssrspdmod, only: stpjo => stpgnssrspd + use gsi_bundlemod, only: gsi_bundle + use bias_predictors, only: predictors + use m_obsNode , only: obsNode + use m_obsLList, only: obsLList_headNode + use kinds, only: r_quad,r_kind,i_kind + implicit none + class(gnssrspdOper ),intent(in):: self + integer(i_kind ),intent(in):: ibin + type(gsi_bundle),intent(in):: dval + type(gsi_bundle),intent(in):: xval + real(r_quad ),dimension(:),intent(inout):: pbcjo ! (1:4) + real(r_kind ),dimension(:),intent(in ):: sges + integer(i_kind),intent(in):: nstep + + type(predictors),target, intent(in):: dbias + type(predictors),target, intent(in):: xbias + + !---------------------------------------- + character(len=*),parameter:: myname_=myname//"::stpjo1_" + class(obsNode),pointer:: headNode + + headNode => obsLList_headNode(self%obsLL(ibin)) + call stpjo(headNode,dval,xval,pbcjo(:),sges,nstep) + headNode => null() + end subroutine stpjo1_ + +end module gsi_gnssrspdOper diff --git a/src/gsi/gsi_obOperTypeManager.F90 b/src/gsi/gsi_obOperTypeManager.F90 index 6db7921905..dfcde80a32 100644 --- a/src/gsi/gsi_obOperTypeManager.F90 +++ b/src/gsi/gsi_obOperTypeManager.F90 @@ -12,7 +12,8 @@ module gsi_obOperTypeManager ! 2018-07-12 j guo - a type-manager for all obOper extensions. ! - an enum mapping of obsinput::dtype(:) to obOper type ! extensions. -! +! 2022-03-15 k apodaca - add GNSS-R L2 ocean wind speed type-manager +! 2023-03-20 k apodaca - add GNSS-R DDM type-manager ! ! input argument list: see Fortran 90 style document below ! ! output argument list: see Fortran 90 style document below @@ -52,6 +53,7 @@ module gsi_obOperTypeManager use gsi_radOper , only: radOper use gsi_rwOper , only: rwOper use gsi_spdOper , only: spdOper + use gsi_gnssrspdOper , only: gnssrspdOper use gsi_sstOper , only: sstOper use gsi_swcpOper , only: swcpOper use gsi_tcamtOper , only: tcamtOper @@ -102,6 +104,7 @@ module gsi_obOperTypeManager public:: iobOper_w public:: iobOper_q public:: iobOper_spd + public:: iobOper_gnssrspd public:: iobOper_rw public:: iobOper_dw public:: iobOper_sst @@ -148,6 +151,7 @@ module gsi_obOperTypeManager enumerator:: iobOper_w enumerator:: iobOper_q enumerator:: iobOper_spd + enumerator:: iobOper_gnssrspd enumerator:: iobOper_rw enumerator:: iobOper_dw enumerator:: iobOper_sst @@ -210,6 +214,7 @@ module gsi_obOperTypeManager type( wOper), target, save:: wOper_mold type( qOper), target, save:: qOper_mold type( spdOper), target, save:: spdOper_mold + type( gnssrspdOper), target, save:: gnssrspdOper_mold type( rwOper), target, save:: rwOper_mold type( dwOper), target, save:: dwOper_mold type( sstOper), target, save:: sstOper_mold @@ -264,6 +269,7 @@ function dtype2index_(dtype) result(index_) case("q" ,"[qoper]" ); index_= iobOper_q case("spd" ,"[spdoper]" ); index_= iobOper_spd + case("gnssrspd" ,"[gnssrspdoper]" ); index_= iobOper_gnssrspd case("rw" ,"[rwoper]" ); index_= iobOper_rw case("dw" ,"[dwoper]" ); index_= iobOper_dw case("sst" ,"[sstoper]" ); index_= iobOper_sst @@ -293,7 +299,7 @@ function dtype2index_(dtype) result(index_) case("ompslp" ); index_= iobOper_o3l case("ompslpuv" ); index_= iobOper_o3l case("ompslpvis"); index_= iobOper_o3l - case("ompslpnc" ); index_= iobOper_o3l + case("ompslpnc" ); index_= iobOper_o3l case("gpsbend","[gpsbendoper]"); index_= iobOper_gpsbend case("gps_bnd"); index_= iobOper_gpsbend @@ -457,6 +463,7 @@ function index2vmold_(iobOper) result(vmold_) case(iobOper_w ); vmold_ => wOper_mold case(iobOper_q ); vmold_ => qOper_mold case(iobOper_spd ); vmold_ => spdOper_mold + case(iobOper_gnssrspd ); vmold_ => gnssrspdOper_mold case(iobOper_rw ); vmold_ => rwOper_mold case(iobOper_dw ); vmold_ => dwOper_mold case(iobOper_sst ); vmold_ => sstOper_mold @@ -573,6 +580,7 @@ subroutine cobstype_config_() cobstype(iobOper_w ) ="wind " ! w_ob_type cobstype(iobOper_q ) ="moisture " ! q_ob_type cobstype(iobOper_spd ) ="wind speed " ! spd_ob_type + cobstype(iobOper_gnssrspd ) ="gnss-r wind speed " ! gnssrspd_ob_type cobstype(iobOper_rw ) ="radial wind " ! rw_ob_type cobstype(iobOper_dw ) ="doppler wind " ! dw_ob_type cobstype(iobOper_sst ) ="sst " ! sst_ob_type diff --git a/src/gsi/gsimain.f90 b/src/gsi/gsimain.f90 index 5eaa38286a..5c56259bc7 100644 --- a/src/gsi/gsimain.f90 +++ b/src/gsi/gsimain.f90 @@ -131,6 +131,8 @@ program gsi ! related to the GOES/GLM lightnig assimilation ! 2019-07-09 todling - add initialization of abstract layer defining use of GFS ensemble ! 2019-08-04 guo - moved ensemble object configuration into module gsi_fixture. +! 2022-03-15 K Apodaca - add CYGNSS and Spire Ocean wind speed observations +! 2023-03-15 K Apodaca - add GNSS-R Doppler Delay Map observations ! ! usage: ! input files: @@ -194,7 +196,7 @@ program gsi ! read_goesglm, read_goesndr, read_gps_ref, read_guess, read_ieeetovs, ! read_lidar, read_obs, read_ozone, read_pcp, read_prepbufr, read_radar, ! read_superwinds, read_wrf_mass_files, read_wrf_mass_guess, -! read_wrf_nmm_files, read_wrf_nmm_guess, rfdpar, rsearch, satthin, +! read_wrf_nmm_files, read_wrf_nmm_guess, read_gnssrspd, rfdpar, rsearch, satthin, ! setupdw, setupoz, setuppcp, setupps, setuppw, setupq, setuprad, ! setupref, setupbend, setuplight, setuprhsall, setuprw, setupspd, setupsst, ! setupt, setupw, simpin1, simpin1_init, smooth121, smoothrf, diff --git a/src/gsi/intgnssrspd.f90 b/src/gsi/intgnssrspd.f90 new file mode 100644 index 0000000000..1310ec7e14 --- /dev/null +++ b/src/gsi/intgnssrspd.f90 @@ -0,0 +1,239 @@ +module intgnssrspdmod +!$$$ module documentation block +! . . . . +! module: intgnssrspdmod module for intgnssrspd and its tangent linear i ntgnssrspd_tl +! (Based on intspd.f90) +! prgmmr: Apodaca, K. - Add CYGNSS TL and AD modules +! email: Karina.Apodaca@Spire.com +! +! abstract: module for intgnssrspd and its tangent linear intgnssrspd_tl +! +! program history log: +! 2023-04-21 k apodaca - initial version +! +! subroutines included: +! sub intgnssrspd_ +! +! variable definitions: +! +! attributes: +! language: f90 +! machine: +! +!$$$ end documentation block + +use m_obsNode, only: obsNode +use m_gnssrspdNode, only: gnssrspdNode +use m_gnssrspdNode, only: gnssrspdNode_typecast +use m_gnssrspdNode, only: gnssrspdNode_nextcast +use m_obsdiagNode, only: obsdiagNode_set +implicit none + +PRIVATE +PUBLIC intgnssrspd + +interface intgnssrspd; module procedure & + intgnssrspd_ +end interface + +contains + +subroutine intgnssrspd_(gnssrspdhead,rval,sval) +!$$$ subprogram documentation block +! . . . . +! subprogram: intgnssrspd apply nonlin qc obs operator for wind speed +! prgmmr: K. Apodaca org: Spire Global, Inc. date: 2023-03-15 +! email: Karina.Apodaca@Spire.com +! Based on intspd.f90 by prgmmr: derber org: np23 date: 1991-02-26 +! +! abstract: apply nonlinear observation operator and adjoint for winds +! +! program history log: +! 2023-04-21 k apodaca - initial version` +! +! input argument list: +! gnssrspdhead - obs type pointer to obs structure +! su - u increment in grid space +! sv - v increment in grid space +! ru +! rv +! +! output argument list: +! gnssrspdhead - obs type pointer to obs structure +! ru - u results from observation operator +! rv - v results from observation operator +! +! attributes: +! language: f90 +! machine: ibm RS/6000 SP +! +!$$$ + use kinds, only: r_kind,i_kind + use obsmod, only: lsaveobsens,l_do_adjoint,luse_obsdiag + use qcmod, only: nlnqc_iter,varqc_iter + use constants, only: zero, half, one, tiny_r_kind,cg_term,r3600 + use gsi_4dvar, only: ltlint + use jfunc, only: jiter + use gsi_bundlemod, only: gsi_bundle + use gsi_bundlemod, only: gsi_bundlegetpointer + use gsi_4dvar, only: ladtest_obs + implicit none + +! Declare passed variables + class(obsNode), pointer, intent(in ) :: gnssrspdhead + type(gsi_bundle), intent(in ) :: sval + type(gsi_bundle), intent(inout) :: rval + +! Declare local variables + integer(i_kind) ier,istatus + integer(i_kind) j1,j2,j3,j4 + real(r_kind) w1,w2,w3,w4,term +! real(r_kind) penalty + real(r_kind) uanl,vanl,gnssrspdanl,gnssrspd,valv,valu + real(r_kind) uatl,vatl,gnssrspdatl,gnssrspdtra,grad + real(r_kind) cg_gnssrspd,p0,wnotgross,wgross,pg_gnssrspd + real(r_kind),pointer,dimension(:) :: su,sv + real(r_kind),pointer,dimension(:) :: ru,rv + type(gnssrspdNode), pointer :: gnssrspdptr + logical :: ltlint_tmp + +! If no gnssrspd data return + if(.not. associated(gnssrspdhead))return + +! Retrieve pointers +! Simply return if any pointer not found + ier=0 + call gsi_bundlegetpointer(sval,'u',su,istatus);ier=istatus+ier + call gsi_bundlegetpointer(sval,'v',sv,istatus);ier=istatus+ier + call gsi_bundlegetpointer(rval,'u',ru,istatus);ier=istatus+ier + call gsi_bundlegetpointer(rval,'v',rv,istatus);ier=istatus+ier + if(ier/=0)return + + if( ladtest_obs) then + ltlint_tmp = ltlint + ltlint = .true. + end if + !gnssrspdptr => gnssrspdhead + gnssrspdptr => gnssrspdNode_typecast(gnssrspdhead) + do while (associated(gnssrspdptr)) + + j1 = gnssrspdptr%ij(1) + j2 = gnssrspdptr%ij(2) + j3 = gnssrspdptr%ij(3) + j4 = gnssrspdptr%ij(4) + w1 = gnssrspdptr%wij(1) + w2 = gnssrspdptr%wij(2) + w3 = gnssrspdptr%wij(3) + w4 = gnssrspdptr%wij(4) + + + valu=zero + valv=zero + gnssrspdtra=sqrt(gnssrspdptr%uges*gnssrspdptr%uges+gnssrspdptr%vges*gnssrspdptr%vges) + + if (ltlint) then + + if (gnssrspdtra>EPSILON(gnssrspdtra)) then +! Forward model + uatl=w1*su(j1)+w2*su(j2)+w3*su(j3)+w4*su(j4) + vatl=w1*sv(j1)+w2*sv(j2)+w3*sv(j3)+w4*sv(j4) + gnssrspdatl=gnssrspdptr%uges*uatl+gnssrspdptr%vges*vatl + gnssrspdatl=gnssrspdatl/gnssrspdtra + + if(luse_obsdiag)then + if (lsaveobsens) then + grad=gnssrspdptr%raterr2*gnssrspdptr%err2*gnssrspdatl + !-- gnssrspdptr%diags%obssen(jiter)=grad + call obsdiagNode_set(gnssrspdptr%diags,jiter=jiter,obssen=grad) + else + !-- if (gnssrspdptr%luse) gnssrspdptr%diags%tldepart(jiter)=gnssrspdatl + if (gnssrspdptr%luse) call obsdiagNode_set(gnssrspdptr%diags,jiter=jiter,tldepart=gnssrspdatl) + endif + endif + + if (l_do_adjoint) then + if (.not. lsaveobsens) then + if( ladtest_obs ) then + grad=gnssrspdatl + else + gnssrspd=gnssrspdatl-gnssrspdptr%res+sqrt(gnssrspdptr%uges*gnssrspdptr%uges+gnssrspdptr%vges*gnssrspdptr%vges) + grad=gnssrspdptr%raterr2*gnssrspdptr%err2*gnssrspd + end if + endif + +! Adjoint + valu=grad*gnssrspdptr%uges/gnssrspdtra + valv=grad*gnssrspdptr%vges/gnssrspdtra + endif + else + if(luse_obsdiag)then + !-- if (gnssrspdptr%luse) gnssrspdptr%diags%tldepart(jiter)=zero + if (gnssrspdptr%luse) call obsdiagNode_set(gnssrspdptr%diags,jiter=jiter,tldepart=zero) + !-- if (lsaveobsens) gnssrspdptr%diags%obssen(jiter)=zero + if (lsaveobsens) call obsdiagNode_set(gnssrspdptr%diags,jiter=jiter,obssen=zero) + end if + endif + + + else ! < ltlint > + +! Forward model + uanl=gnssrspdptr%uges+w1* su(j1)+w2* su(j2)+w3* su(j3)+w4* su(j4) + vanl=gnssrspdptr%vges+w1* sv(j1)+w2* sv(j2)+w3* sv(j3)+w4* sv(j4) + gnssrspdanl=sqrt(uanl*uanl+vanl*vanl) + !-- if (luse_obsdiag .and. gnssrspdptr%luse) gnssrspdptr%diags%tldepart(jiter)=gnssrspdanl-gnssrspdtra + if (luse_obsdiag .and. gnssrspdptr%luse) call obsdiagNode_set(gnssrspdptr%diags,jiter=jiter,tldepart=gnssrspdanl-gnssrspdtra) + + if (l_do_adjoint) then + valu=zero + valv=zero + gnssrspd=gnssrspdanl-gnssrspdptr%res + grad=gnssrspdptr%raterr2*gnssrspdptr%err2*gnssrspd + +! Adjoint +! if(gnssrspdanl > tiny_r_kind*100._r_kind) then + if (gnssrspdanl>EPSILON(gnssrspdanl)) then + !-- if (luse_obsdiag .and. lsaveobsens) gnssrspdptr%diags%obssen(jiter)=grad + if (luse_obsdiag .and. lsaveobsens) call obsdiagNode_set(gnssrspdptr%diags,jiter=jiter,obssen=grad) + valu=uanl/gnssrspdanl + valv=vanl/gnssrspdanl + if (nlnqc_iter .and. gnssrspdptr%pg > tiny_r_kind .and. & + gnssrspdptr%b > tiny_r_kind) then + pg_gnssrspd=gnssrspdptr%pg*varqc_iter + cg_gnssrspd=cg_term/gnssrspdptr%b + wnotgross= one-pg_gnssrspd + wgross = pg_gnssrspd*cg_gnssrspd/wnotgross + p0 = wgross/(wgross+exp(-half*gnssrspdptr%err2*gnssrspd**2)) + term = (one-p0) + grad = grad*term + endif + end if + endif ! < l_do_adjoint > + + valu=valu*grad + valv=valv*grad + + endif ! < ltlint > + + + if (l_do_adjoint) then + ru(j1)=ru(j1)+w1*valu + ru(j2)=ru(j2)+w2*valu + ru(j3)=ru(j3)+w3*valu + ru(j4)=ru(j4)+w4*valu + rv(j1)=rv(j1)+w1*valv + rv(j2)=rv(j2)+w2*valv + rv(j3)=rv(j3)+w3*valv + rv(j4)=rv(j4)+w4*valv + + endif + + !gnssrspdptr => gnssrspdptr%llpoint + gnssrspdptr => gnssrspdNode_nextcast(gnssrspdptr) + + end do + if( ladtest_obs) ltlint = ltlint_tmp + return +end subroutine intgnssrspd_ + +end module intgnssrspdmod diff --git a/src/gsi/intjo.f90 b/src/gsi/intjo.f90 index a68355471b..4ce2952062 100644 --- a/src/gsi/intjo.f90 +++ b/src/gsi/intjo.f90 @@ -14,6 +14,7 @@ module intjomod ! 2016-08-29 J Guo - Separated calls to intozlay() and intozlev() ! 2018-08-10 guo - a new generic intjo() implementation replaced type ! specific intXYZ() calls with polymorphic %intjo(). +! 2022-03-15 K Apodaca - add CYGNSS and Spire ocean wind speed ! ! subroutines included: ! sub intjo_ @@ -31,7 +32,7 @@ module intjomod use gsi_obOperTypeManager, only: & iobOper_t, iobOper_pw, iobOper_q, & iobOper_cldtot, iobOper_w, iobOper_dw, & - iobOper_rw, iobOper_dbz, iobOper_fed, & + iobOper_rw, iobOper_dbz, iobOper_fed, iobOper_gnssrspd, & iobOper_spd, iobOper_oz, iobOper_o3l, iobOper_colvk, & iobOper_pm2_5, iobOper_pm10, iobOper_ps, iobOper_tcp, iobOper_sst, & iobOper_gpsbend, iobOper_gpsref, & @@ -60,7 +61,7 @@ module intjomod integer(i_kind),parameter,dimension(obOper_count):: ix_obtype = (/ & iobOper_t, iobOper_pw, iobOper_q, & iobOper_cldtot, iobOper_w, iobOper_dw, & - iobOper_rw, iobOper_dbz, iobOper_fed, & + iobOper_rw, iobOper_dbz, iobOper_fed, iobOper_gnssrspd, & iobOper_spd, iobOper_oz, iobOper_o3l, iobOper_colvk, & iobOper_pm2_5, iobOper_pm10, iobOper_ps, iobOper_tcp, iobOper_sst, & iobOper_gpsbend, iobOper_gpsref, & diff --git a/src/gsi/m_gnssrspdNode.F90 b/src/gsi/m_gnssrspdNode.F90 new file mode 100644 index 0000000000..1a8400440f --- /dev/null +++ b/src/gsi/m_gnssrspdNode.F90 @@ -0,0 +1,257 @@ +module m_gnssrspdNode +!$$$ subprogram documentation block +! . . . . +! subprogram: module m_gnssrspdNode +! prgmmr: K. Apodaca +! org: Spire Global, Inc. +! date: 2023-03-05 +! (based on m_*Node.F90 modules by: +! prgmmr: j guo +! org: NASA/GSFC, Global Modeling and Assimilation Office, 610.3 +! date: 2016-05-18 +! +! abstract: class-module of obs-type gnssrspdNode (wind speed) +! +! program history log: +! 2016-05-18 j guo - added this document block for the initial polymorphic +! implementation. +! +! input argument list: see Fortran 90 style document below +! +! output argument list: see Fortran 90 style document below +! +! attributes: +! language: Fortran 90 and/or above +! machine: +! +!$$$ end subprogram documentation block + +! module interface: + use m_obsdiagNode, only: obs_diag + use m_obsdiagNode, only: obs_diags + use kinds , only: i_kind,r_kind + use mpeu_util, only: assert_,die,perr,warn,tell + use m_obsNode, only: obsNode + implicit none + private + + public:: gnssrspdNode + + type,extends(obsNode):: gnssrspdNode + !type(gnssrspd_ob_type),pointer :: llpoint => NULL() + type(obs_diag), pointer :: diags => NULL() + real(r_kind) :: res ! speed observation + real(r_kind) :: err2 ! speed error squared + real(r_kind) :: raterr2 ! square of ratio of final obs error + ! to original obs error + !real(r_kind) :: time ! observation time in sec + real(r_kind) :: b ! variational quality control parameter + real(r_kind) :: pg ! variational quality control parameter + real(r_kind) :: wij(4) ! horizontal interpolation weights + real(r_kind) :: uges ! guess u value + real(r_kind) :: vges ! guess v value + integer(i_kind) :: ij(4) ! horizontal locations + !logical :: luse ! flag indicating if ob is used in pen. + + !integer(i_kind) :: idv,iob ! device id and obs index for sorting + !real (r_kind) :: elat, elon ! earth lat-lon for redistribution + !real (r_kind) :: dlat, dlon ! earth lat-lon for redistribution + real (r_kind) :: factw ! factor of 10m wind + contains + procedure,nopass:: mytype + procedure:: setHop => obsNode_setHop_ + procedure:: xread => obsNode_xread_ + procedure:: xwrite => obsNode_xwrite_ + procedure:: isvalid => obsNode_isvalid_ + procedure:: gettlddp => gettlddp_ + + ! procedure, nopass:: headerRead => obsHeader_read_ + ! procedure, nopass:: headerWrite => obsHeader_write_ + ! procedure:: init => obsNode_init_ + ! procedure:: clean => obsNode_clean_ + end type gnssrspdNode + + public:: gnssrspdNode_typecast + public:: gnssrspdNode_nextcast + interface gnssrspdNode_typecast; module procedure typecast_ ; end interface + interface gnssrspdNode_nextcast; module procedure nextcast_ ; end interface + + public:: gnssrspdNode_appendto + interface gnssrspdNode_appendto; module procedure appendto_ ; end interface + + character(len=*),parameter:: MYNAME="m_gnssrspdNode" + +#include "myassert.H" +#include "mytrace.H" +contains +function typecast_(aNode) result(ptr_) +!-- cast a class(obsNode) to a type(gnssrspdNode) + use m_obsNode, only: obsNode + implicit none + type(gnssrspdNode ),pointer:: ptr_ + class(obsNode),pointer,intent(in):: aNode + ptr_ => null() + if(.not.associated(aNode)) return + ! logically, typecast of a null-reference is a null pointer. + select type(aNode) + type is(gnssrspdNode) + ptr_ => aNode + end select +return +end function typecast_ + +function nextcast_(aNode) result(ptr_) +!-- cast an obsNode_next(obsNode) to a type(gnssrspdNode) + use m_obsNode, only: obsNode,obsNode_next + implicit none + type(gnssrspdNode ),pointer:: ptr_ + class(obsNode),target ,intent(in):: aNode + + class(obsNode),pointer:: inode_ + inode_ => obsNode_next(aNode) + ptr_ => typecast_(inode_) +return +end function nextcast_ + +subroutine appendto_(aNode,oll) +!-- append aNode to linked-list oLL + use m_obsNode , only: obsNode + use m_obsLList, only: obsLList,obsLList_appendNode + implicit none + type(gnssrspdNode),pointer,intent(in):: aNode + type(obsLList),intent(inout):: oLL + + class(obsNode),pointer:: inode_ + inode_ => aNode + call obsLList_appendNode(oLL,inode_) + inode_ => null() +end subroutine appendto_ + +! obsNode implementations + +function mytype() + implicit none + character(len=:),allocatable:: mytype + mytype="[gnssrspdNode]" +end function mytype + +subroutine obsNode_xread_(aNode,iunit,istat,diagLookup,skip) + use m_obsdiagNode, only: obsdiagLookup_locate + implicit none + class(gnssrspdNode),intent(inout):: aNode + integer(i_kind),intent(in ):: iunit + integer(i_kind),intent( out):: istat + type(obs_diags),intent(in ):: diagLookup + logical,optional,intent(in ):: skip + + character(len=*),parameter:: myname_=MYNAME//'.obsNode_xread_' + logical:: skip_ +_ENTRY_(myname_) + skip_=.false. + if(present(skip)) skip_=skip + + istat=0 + if(skip_) then + read(iunit,iostat=istat) + if(istat/=0) then + call perr(myname_,'skipping read(%(res,err2,...)), iostat =',istat) + _EXIT_(myname_) + return + endif + + else + read(iunit,iostat=istat) aNode%res , & + aNode%err2 , & + aNode%raterr2, & + aNode%b , & + aNode%pg , & + aNode%uges , & + aNode%vges , & + aNode%factw , & + aNode%wij , & + aNode%ij + if (istat/=0) then + call perr(myname_,'read(%(res,err2,...)), iostat =',istat) + _EXIT_(myname_) + return + end if + + aNode%diags => obsdiagLookup_locate(diagLookup,aNode%idv,aNode%iob,1_i_kind) + if(.not.associated(aNode%diags)) then + call perr(myname_,'obsdiagLookup_locate(), %idv =',aNode%idv) + call perr(myname_,' %iob =',aNode%iob) + call die(myname_) + endif + endif +_EXIT_(myname_) +return +end subroutine obsNode_xread_ + +subroutine obsNode_xwrite_(aNode,junit,jstat) + implicit none + class(gnssrspdNode),intent(in):: aNode + integer(i_kind),intent(in ):: junit + integer(i_kind),intent( out):: jstat + + character(len=*),parameter:: myname_=MYNAME//'.obsNode_xwrite_' +_ENTRY_(myname_) + + jstat=0 + write(junit,iostat=jstat) aNode%res , & + aNode%err2 , & + aNode%raterr2, & + aNode%b , & + aNode%pg , & + aNode%uges , & + aNode%vges , & + aNode%factw , & + aNode%wij , & + aNode%ij + if (jstat/=0) then + call perr(myname_,'write(%(res,err2,...)), iostat =',jstat) + _EXIT_(myname_) + return + end if +_EXIT_(myname_) +return +end subroutine obsNode_xwrite_ + +subroutine obsNode_setHop_(aNode) + use m_cvgridLookup, only: cvgridLookup_getiw + implicit none + class(gnssrspdNode),intent(inout):: aNode + + character(len=*),parameter:: myname_=MYNAME//'::obsNode_setHop_' +_ENTRY_(myname_) + call cvgridLookup_getiw(aNode%elat,aNode%elon,aNode%ij,aNode%wij) + aNode%wij(1:4) = aNode%wij(1:4)*aNode%factw +_EXIT_(myname_) +return +end subroutine obsNode_setHop_ + +function obsNode_isvalid_(aNode) result(isvalid_) + implicit none + logical:: isvalid_ + class(gnssrspdNode),intent(in):: aNode + + character(len=*),parameter:: myname_=MYNAME//'::obsNode_isvalid_' +_ENTRY_(myname_) + isvalid_=associated(aNode%diags) +_EXIT_(myname_) +return +end function obsNode_isvalid_ + +pure subroutine gettlddp_(aNode,jiter,tlddp,nob) + use kinds, only: r_kind + implicit none + class(gnssrspdNode), intent(in):: aNode + integer(kind=i_kind),intent(in):: jiter + real(kind=r_kind),intent(inout):: tlddp + integer(kind=i_kind),optional,intent(inout):: nob + + tlddp = tlddp + aNode%diags%tldepart(jiter)*aNode%diags%tldepart(jiter) + if(present(nob)) nob=nob+1 +return +end subroutine gettlddp_ + +end module m_gnssrspdNode diff --git a/src/gsi/m_obsNodeTypeManager.F90 b/src/gsi/m_obsNodeTypeManager.F90 index 43b42e4bf2..9da356c808 100644 --- a/src/gsi/m_obsNodeTypeManager.F90 +++ b/src/gsi/m_obsNodeTypeManager.F90 @@ -15,6 +15,8 @@ module m_obsNodeTypeManager ! 2018-01-23 k apodaca - add a new observation type i.e. lightning (light) ! suitable for the GOES/GLM instrument ! +! 2024-04-07 k apodaca - add gnssrspdnode +! ! input argument list: see Fortran 90 style document below ! ! output argument list: see Fortran 90 style document below @@ -54,6 +56,8 @@ module m_obsNodeTypeManager use m_wspd10mNode, only: wspd10mNode use m_uwnd10mNode, only: uwnd10mNode use m_vwnd10mNode, only: vwnd10mNode + + use m_gnssrspdNode, only: gnssrspdNode use m_td2mNode , only: td2mNode use m_mxtmNode , only: mxtmNode @@ -122,10 +126,10 @@ module m_obsNodeTypeManager public:: iobsNode_cldch public:: iobsNode_swcp public:: iobsNode_lwcp - public:: iobsNode_light public:: iobsNode_dbz public:: iobsNode_fed + public:: iobsNode_gnssrspd public :: obsNode_typeMold public :: obsNode_typeIndex @@ -166,6 +170,8 @@ module m_obsNodeTypeManager type(wspd10mNode), target, save:: wspd10m_mold type(uwnd10mNode), target, save:: uwnd10m_mold type(vwnd10mNode), target, save:: vwnd10m_mold + + type(gnssrspdNode), target, save:: gnssrspd_mold type( td2mNode), target, save:: td2m_mold type( mxtmNode), target, save:: mxtm_mold @@ -214,6 +220,7 @@ module m_obsNodeTypeManager enumerator:: iobsNode_w enumerator:: iobsNode_q enumerator:: iobsNode_spd + enumerator:: iobsNode_gnssrspd enumerator:: iobsNode_rw enumerator:: iobsNode_dw enumerator:: iobsNode_sst @@ -301,6 +308,8 @@ function vname2index_(vname) result(index_) "[uwnd10mnode]"); index_ = iobsNode_uwnd10m case("vwnd10m", & "[vwnd10mnode]"); index_ = iobsNode_vwnd10m + case("gnssrspd", & + "[gnssrspdnode]"); index_ = iobsNode_gnssrspd case("td2m" , "[td2mnode]"); index_ = iobsNode_td2m case("mxtm" , "[mxtmnode]"); index_ = iobsNode_mxtm @@ -365,6 +374,8 @@ function vmold2index_select_(mold) result(index_) type is(wspd10mNode); index_ = iobsNode_wspd10m type is(uwnd10mNode); index_ = iobsNode_uwnd10m type is(vwnd10mNode); index_ = iobsNode_vwnd10m + + type is(gnssrspdNode); index_ = iobsNode_gnssrspd type is( td2mNode); index_ = iobsNode_td2m type is( mxtmNode); index_ = iobsNode_mxtm @@ -423,6 +434,8 @@ function index2vmold_(i_obType) result(obsmold_) case(iobsNode_wspd10m); obsmold_ => wspd10m_mold case(iobsNode_uwnd10m); obsmold_ => uwnd10m_mold case(iobsNode_vwnd10m); obsmold_ => vwnd10m_mold + + case(iobsNode_gnssrspd ); obsmold_ => gnssrspd_mold case(iobsNode_td2m ); obsmold_ => td2m_mold case(iobsNode_mxtm ); obsmold_ => mxtm_mold diff --git a/src/gsi/m_rhs.F90 b/src/gsi/m_rhs.F90 index aea417fe27..723b0a6918 100644 --- a/src/gsi/m_rhs.F90 +++ b/src/gsi/m_rhs.F90 @@ -17,6 +17,7 @@ module m_rhs ! through an enum block. ! - removed external dimension argument aworkdim2 of ! rhs_alloc(). +! 2022-03-15 K. Apodaca - add GNSS-R L2 ocean wind speed index (i_gnssrspd) ! ! input argument list: see Fortran 90 style document below ! @@ -66,6 +67,7 @@ module m_rhs public:: i_vis public:: i_pblh public:: i_wspd10m + public:: i_gnssrspd public:: i_td2m public:: i_mxtm public:: i_mitm @@ -133,6 +135,7 @@ module m_rhs enumerator:: i_vis enumerator:: i_pblh enumerator:: i_wspd10m + enumerator:: i_gnssrspd enumerator:: i_td2m enumerator:: i_mxtm enumerator:: i_mitm @@ -207,7 +210,6 @@ subroutine rhs_alloc() allocate(rhs_stats_oz(9,jpch_oz)) allocate(rhs_toss_gps(max(1,nprof_gps))) - rhs_awork =zero rhs_bwork =zero rhs_aivals =zero diff --git a/src/gsi/obsmod.F90 b/src/gsi/obsmod.F90 index 1c45c62bc8..19f6a072fb 100644 --- a/src/gsi/obsmod.F90 +++ b/src/gsi/obsmod.F90 @@ -160,6 +160,7 @@ module obsmod ! 2021-11-16 Zhao - add option l_obsprvdiag (if true) to trigger the output of ! observation provider and sub-provider information into ! obsdiags files (used for AutoObsQC) +! 2022-03-15 K. Apodaca - add GNSS-R L2 ocean wind speed observations (CYGNSS, Spire) ! 2023-07-10 Y. Wang, D. Dowell - add variables for flash extent density ! 2023-10-10 H. Wang (GSL) - add variables for flash extent density EnVar DA ! @@ -463,8 +464,8 @@ module obsmod public :: iout_lag,iout_uv,iout_gps,iout_ps,iout_light,mype_light public :: mype_gust,mype_vis,mype_pblh,iout_gust,iout_vis,iout_pblh public :: mype_tcamt,mype_lcbas,iout_tcamt,iout_lcbas - public :: mype_wspd10m,mype_td2m,iout_wspd10m,iout_td2m - public :: mype_uwnd10m,mype_vwnd10m,iout_uwnd10m,iout_vwnd10m + public :: mype_wspd10m,mype_gnssrspd,mype_td2m,iout_wspd10m,iout_gnssrspd,iout_td2m + public :: mype_uwnd10m,mype_vwnd10m,iout_uwnd10m,iout_vwnd10m public :: mype_mxtm,mype_mitm,iout_mxtm,iout_mitm public :: mype_pmsl,mype_howv,iout_pmsl,iout_howv public :: mype_swcp,mype_lwcp,iout_swcp,iout_lwcp @@ -600,12 +601,12 @@ module obsmod integer(i_kind) iout_dw,iout_gps,iout_sst,iout_tcp,iout_lag integer(i_kind) iout_co,iout_gust,iout_vis,iout_pblh,iout_tcamt,iout_lcbas integer(i_kind) iout_cldch - integer(i_kind) iout_wspd10m,iout_td2m,iout_mxtm,iout_mitm,iout_pmsl,iout_howv + integer(i_kind) iout_wspd10m,iout_gnssrspd,iout_td2m,iout_mxtm,iout_mitm,iout_pmsl,iout_howv integer(i_kind) iout_uwnd10m,iout_vwnd10m,iout_fed integer(i_kind) mype_t,mype_q,mype_uv,mype_ps,mype_pw, & mype_rw,mype_dw,mype_gps,mype_sst, & mype_tcp,mype_lag,mype_co,mype_gust,mype_vis,mype_pblh, & - mype_wspd10m,mype_td2m,mype_mxtm,mype_mitm,mype_pmsl,mype_howv,& + mype_wspd10m,mype_gnssrspd,mype_td2m,mype_mxtm,mype_mitm,mype_pmsl,mype_howv,& mype_uwnd10m,mype_vwnd10m, mype_tcamt,mype_lcbas, mype_dbz, mype_fed integer(i_kind) mype_cldch integer(i_kind) iout_swcp, iout_lwcp @@ -887,6 +888,7 @@ subroutine init_obsmod_dflts iout_light=237 ! lightning iout_dbz=238 ! radar reflectivity iout_fed=239 ! flash extent density + iout_gnssrspd=240 ! GNSS-R wind speed mype_ps = npe-1 ! surface pressure mype_t = max(0,npe-2) ! temperature @@ -922,6 +924,7 @@ subroutine init_obsmod_dflts mype_light=max(0,npe-32)! GOES/GLM lightning mype_dbz=max(0,npe-33) ! radar reflectivity mype_fed= max(0,npe-34) ! flash extent density + mype_gnssrspd= max(0,npe-35) ! surface GNSS-R speed ! Initialize arrays used in namelist obs_input diff --git a/src/gsi/read_gnssrspd.f90 b/src/gsi/read_gnssrspd.f90 new file mode 100644 index 0000000000..859345a9cc --- /dev/null +++ b/src/gsi/read_gnssrspd.f90 @@ -0,0 +1,463 @@ +subroutine read_gnssrspd(nread,ndata,nodata,infile,obstype,lunout,twind,sis,& + nobs) + +!$$$ subprogram documentation block +! . . . . +! subprogram: read_gnssrspd read obs from gnssrspd bufr file +! prgmmr: kapodaca org: Spire Global, Inc. date: 2022-03-12 +! Largely based on other read_* routines. For the !Winds --- surface wind speed section +! ~L438, we used the exixting read_fl_hdob.f90 routine +! +! abstract: This routine reads GNSSRSPD L2 wind speed observations +! +! When running the gsi in regional mode, the code only +! retains those observations that fall within the regional +! domain + +! program history log: +! 2022-03-12 k apodaca- initial coding +! +! input argument list: +! infile - unit from which to read BUFR data +! obstype - observation type to process +! lunout - unit to which to write data for further processing +! twind - input group time window (hours) +! sis - satellite/instrument/sensor indicator +! +! output argument list: +! nread - number of type "obstype" observations read +! nodata - number of individual "obstype" observations read +! ndata - number of type "obstype" observations retained for further processing +! nobs - array of observations on each subdomain for each processor +! +! attributes: +! language: f90 +! machine: +! +!$$$ + use kinds, only: r_single,r_kind,r_double,i_kind + use constants, only: zero,one_tenth,one,two,ten,deg2rad,t0c,half,& + three,four,rad2deg,r0_01,& + r60inv,r10,r100,r2000,hvap + use gridmod, only: diagnostic_reg,regional,nlon,nlat,& + tll2xy,txy2ll,rotate_wind_ll2xy,rotate_wind_xy2ll,& + rlats,rlons,twodvar_regional + use convinfo, only: nconvtype, & + icuse,ictype,ioctype + use obsmod, only: ran01dom + use obsmod, only: iadate,bmiss,offtime_data + use gsi_4dvar, only: l4dvar,l4densvar,time_4dvar,winlen + use qcmod, only: errormod + use convthin, only: make3grids,del3grids + use ndfdgrids,only: init_ndfdgrid,destroy_ndfdgrid,relocsfcob,adjust_error + use deter_sfc_mod, only: deter_sfc_type,deter_sfc2 + use mpimod, only: npe + + implicit none + +! Declare passed variables + character(len=*), intent(in ) :: infile,obstype + character(len=20),intent(in ) :: sis + integer(i_kind) , intent(in ) :: lunout + integer(i_kind) , dimension(npe), intent(inout) :: nobs + integer(i_kind) , intent(inout) :: nread,ndata,nodata + real(r_kind) , intent(in ) :: twind + +! Declare local variables +! Logical variables + logical :: outside + logical :: inflate_error + logical :: lspdob + +! Character variables + character(40) :: timestr,locstr,wndstr,oestr + character( 8) :: subset + character( 8) :: c_prvstg,c_sprvstg + character( 8) :: c_station_id + character(10) date +! Integer variables + integer(i_kind), parameter :: mxib = 31 + + integer(i_kind) :: i,k + integer(i_kind) :: ihh,idd,idate,im,iy + integer(i_kind) :: lunin + integer(i_kind) :: ireadmg,ireadsb + integer(i_kind) :: ilat,ilon + integer(i_kind) :: nlv + integer(i_kind) :: nreal,nchanl + integer(i_kind) :: idomsfc,isflg + integer(i_kind) :: iout + integer(i_kind) :: nc,ncsave + integer(i_kind) :: ntmatch,ntb + integer(i_kind) :: nmsg + integer(i_kind) :: maxobs + integer(i_kind) :: itype,iecol + integer(i_kind) :: qcm,lim_qm + integer(i_kind) :: wspd_qm + integer(i_kind) :: ntest,nvtest + integer(i_kind) :: igood + integer(i_kind) :: iuse + + integer(i_kind) :: idate5(5) + integer(i_kind) :: minobs,minan + + integer(i_kind), allocatable,dimension(:) :: isort + +! Real variables + real(r_kind), parameter :: r0_001 = 0.001_r_kind + real(r_kind), parameter :: r1_2 = 1.2_r_kind + real(r_kind), parameter :: r3_0 = 3.0_r_kind + real(r_kind), parameter :: r0_7 = 0.7_r_kind + real(r_kind), parameter :: r6 = 6.0_r_kind + real(r_kind), parameter :: r50 = 50.0_r_kind + real(r_kind), parameter :: r1200 = 1200.0_r_kind + real(r_kind), parameter :: emerr = 0.2_r_kind ! RH + real(r_kind), parameter :: missing = 1.0e+11_r_kind + real(r_kind), parameter :: r180 = 180.0_r_kind + real(r_kind), parameter :: r360 = 360.0_r_kind + + real(r_kind) :: toff,t4dv + real(r_kind) :: usage + real(r_kind) :: woe,obserr + real(r_kind) :: dlat,dlon,dlat_earth,dlon_earth + real(r_kind) :: dlat_earth_deg,dlon_earth_deg + real(r_kind) :: cdist,disterr,disterrmax,rlon00,rlat00 + real(r_kind) :: vdisterrmax + real(r_kind) :: spdob + real(r_kind) :: gob + real(r_kind) :: dlnpsob + real(r_kind) :: tdiff + real(r_kind) :: tsavg,ff10,sfcr,zz + real(r_kind) :: errmin + real(r_kind) :: log100 + + real(r_kind) :: obstime(6,1) + real(r_kind) :: obsloc(2,1) + real(r_kind) :: gnssrw(2,1) + + real(r_double) :: rstation_id + real(r_double) :: r_prvstg(1,1),r_sprvstg(1,1) + + real(r_kind), allocatable,dimension(:,:) :: cdata_all,cdata_out + +! Equivalence to handle character names + equivalence(r_prvstg(1,1),c_prvstg) + equivalence(r_sprvstg(1,1),c_sprvstg) + equivalence(rstation_id,c_station_id) + +! Data + !data timestr / 'YEAR MNTH DAYS HOUR MINU SECO' / + data timestr / 'DHR RPT' / + !data locstr / 'CLAT CLON' / + data locstr / 'XOB YOB' / + !data wndstr / 'WSPD' / !GNSSRSPD Wind speed + data wndstr / 'SOB' / !GNSSRSPD Wind speed + data oestr / 'WSU' / !GNSSRSPD Wind speed uncertainty/error + data lunin / 13 / + +!------------------------------------------------------------------------------------------------ + + write(6,*)'READ_GNSSRSPD: begin to read gnssrspd satellite data ...' + +! Initialize parameters + +! Set common variables + lspdob = obstype == 'gnssrspd' + + nreal = 0 + iecol = 0 + log100=log(100._r_kind) + + + lim_qm = 4 + iecol=0 + if (lspdob) then + nreal = 23 + iecol = 4 + errmin = one + else + write(6,*) ' illegal obs type in read_gnssrspd ' + call stop2(94) + end if + + inflate_error = .true. + +! Check if the obs type specified in the convinfo is in the fl hdob bufr file +! If found, get the index (nc) from the convinfo for the specified type +ntmatch = 0 +ncsave = 0 +do nc = 1, nconvtype + if (trim(ioctype(nc)) == trim(obstype)) then + if (trim(ioctype(nc)) == 'gnssrspd' .and. ictype(nc) == 600 ) then + ntmatch = ntmatch + 1 + ncsave = nc + itype = ictype(nc) + end if + end if +end do +if (ntmatch == 0) then ! Return if not specified in convinfo + write(6,*) ' READ_GNSSRSPD: No matching obstype found in obsinfo ', obstype + return +else + nc = ncsave + write(6,*) ' READ_GNSSRSPD: Processing GNSSRSPD data : ', ntmatch, nc, ioctype(nc), ictype(nc), itype +end if + + +!------------------------------------------------------------------------------------------------ + +! Go through the bufr file to find out how mant subsets to process + nmsg = 0 + maxobs = 0 + call closbf(lunin) + open(lunin,file=trim(infile),form='unformatted') + call openbf(lunin,'IN',lunin) + call datelen(10) + + loop_msg1: do while(ireadmg(lunin,subset,idate) >= 0) + if(nmsg == 0) call time_4dvar(idate,toff) ! time offset (hour) + + nmsg = nmsg+1 + loop_readsb1: do while(ireadsb(lunin) == 0) + maxobs = maxobs+1 + end do loop_readsb1 + end do loop_msg1 + call closbf(lunin) + write(6,*) 'READ_GNSSRSPD: total number of data found in the bufr file ',maxobs,obstype + write(6,*) 'READ_GNSSRSPD: time offset is ',toff,' hours' + +!--------------------------------------------------------------------------------------------------- + +! Allocate array to hold data + allocate(cdata_all(nreal,maxobs)) + allocate(isort(maxobs)) + +! Initialize + cdata_all = zero + isort = 0 + nread = 0 + nchanl = 0 + ntest = 0 + nvtest = 0 + ilon = 2 + ilat = 3 + +! Open bufr file again for reading + call closbf(lunin) + open(lunin,file=trim(infile),form='unformatted') + call openbf(lunin,'IN',lunin) + call datelen(10) + ntb = 0 + igood = 0 +! Loop through BUFR file + loop_msg2: do while(ireadmg(lunin,subset,idate) >= 0) + loop_readsb2: do while(ireadsb(lunin) == 0) + + ntb = ntb+1 + + c_station_id = subset + +! QC mark 9: will be monitored but not assimilated +! QC mark 4: reject - will not be monitored nor assimilated +! QC mark 3: suspect +! QC mark 2: neutral or not checked +! QC mark 1: good +! QC mark 0: keep - will be always assimilated + qcm = 0 + wspd_qm = 0 + + +! Read observation time + call ufbint(lunin,obstime,2,1,nlv,timestr) + +! If date in gnssrspd file does not agree with analysis date, +! print message and stop program execution. + write(date,'( i10)') idate + read (date,'(i4,3i2)') iy,im,idd,ihh + if(offtime_data) then + +! in time correction for observations to account for analysis + idate5(1)=iy + idate5(2)=im + idate5(3)=idd + idate5(4)=ihh + idate5(5)=0 + call w3fs21(idate5,minobs) ! obs ref time in minutes relative to historic date + idate5(1)=iadate(1) + idate5(2)=iadate(2) + idate5(3)=iadate(3) + idate5(4)=iadate(4) + idate5(5)=0 + call w3fs21(idate5,minan) ! analysis ref time in minutes relative to historic date +! Add obs reference time, then subtract analysis time to get obs time relative to analysis + + tdiff=float(minobs-minan)*r60inv + + else + tdiff=zero + end if + + t4dv = toff+obstime(1,1) + + if (l4dvar.or.l4densvar) then + if (t4dv < zero .OR. t4dv > winlen) cycle loop_readsb2 + else + if (abs(tdiff)>twind) cycle loop_readsb2 + endif + nread = nread+1 + + usage = zero ! will be considered for assimilation + ! subject to further QC in setupt subroutine + iuse = icuse(nc) ! assimilation flag + if (iuse <=0) usage = r100 ! will be monitored but not assimilated + +! Read observation location (lat/lon degree) + call ufbint(lunin,obsloc,2,1,nlv,locstr) + + if (obsloc(1,1) == missing .or. abs(obsloc(1,1)) < -180.0_r_kind .or. & + obsloc(1,1) == missing .or. abs(obsloc(1,1)) > 180.0_r_kind .or. & + obsloc(2,1) == missing .or. abs(obsloc(2,1)) < -90.0_r_kind .or. & + obsloc(2,1) == missing .or. abs(obsloc(2,1)) > 90.0_r_kind ) then + write(6,*) 'READ_GNSSRSPD: bad lon/lat values: ', obsloc(1,1),obsloc(2,1) + cycle loop_readsb2 + endif +! GNSSRSPD BUFR longitudes are in the +/- 180 range, need to convert to 0 to 360 +! deg + if (obsloc(1,1) < 0.0_r_kind) obsloc(1,1) = obsloc(1,1) + 360.0_r_kind + !if (obsloc(1,1) >= 0.0_r_kind .or. obsloc(1,1) <= 180.0_r_kind) obsloc(1,1) = obsloc(1,1) + 180.0_r_kind + + dlon_earth_deg = obsloc(1,1) + dlat_earth_deg = obsloc(2,1) + dlon_earth = obsloc(1,1)*deg2rad ! degree to radian + dlat_earth = obsloc(2,1)*deg2rad ! degree to radian + +! Convert obs lat/lon to rotated coordinate and check +! if the obs is outside of the regional domain + if (regional) then + call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) + if (diagnostic_reg) then + call txy2ll(dlon,dlat,rlon00,rlat00) + ntest = ntest+1 + cdist = sin(dlat_earth)*sin(rlat00)+cos(dlat_earth)*cos(rlat00)* & + (sin(dlon_earth)*sin(rlon00)+cos(dlon_earth)*cos(rlon00)) + cdist = max(-one,min(cdist,one)) + disterr = acos(cdist)*rad2deg + disterrmax = max(disterrmax,disterr) + end if + if(outside) cycle loop_readsb2 + else + dlon = dlon_earth + dlat = dlat_earth + call grdcrd1(dlat,rlats,nlat,1) + call grdcrd1(dlon,rlons,nlon,1) + endif + + +! Read surface wind speed [m/s] from GNSSRSPD + if (lspdob) then + !usage=r100 + !usage = zero ! will be considered for assimilation + + ! Get Wind Speed observations from bufr file + call ufbint(lunin,gnssrw,1,1,nlv,wndstr) + spdob = gnssrw(1,1) ! surface wind speed + + ! Don't permit observations with ws <= 1 m/s + if (spdob <= one) cycle loop_readsb2 + if (spdob >= missing) cycle loop_readsb2 + + ! Get observation error from bufr file + call ufbint(lunin,gnssrw,1,1,nlv,oestr) + obserr = max(gnssrw(1,1),1.5_r_kind) ! surface wind speed error + endif + + + if ( .not. twodvar_regional) then + call deter_sfc_type(dlat_earth,dlon_earth,t4dv,isflg,tsavg) + endif + + ! Get information from surface file necessary for conventional data + call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,tsavg,ff10,sfcr,zz) + ! Process data passed quality control + igood = igood + 1 + ndata = ndata + 1 + nodata = nodata + 1 + iout = ndata + +! Read extrapolated surface pressure [pa] and convert to [cb] + dlnpsob = log100 ! default (1000mb) + +!------------------------------------------------------------------------------------------------- + + ! Winds --- surface wind speed + if (lspdob) then + woe = obserr + !if (inflate_error) woe = woe*r3_0 + !if (inflate_error) woe = woe*r1_2 + !if (qcm > lim_qm ) woe = woe*1.0e6_r_kind + cdata_all( 1,iout)=woe ! wind error + cdata_all( 2,iout)=dlon ! grid relative longitude + cdata_all( 3,iout)=dlat ! grid relative latitude + cdata_all( 4,iout)=dlnpsob ! ln(surface pressure in cb) + cdata_all( 5,iout)=spdob*sqrt(two)*half ! u obs + cdata_all( 6,iout)=spdob*sqrt(two)*half ! v obs + cdata_all( 7,iout)=rstation_id ! station id + cdata_all( 8,iout)=t4dv ! time + cdata_all( 9,iout)=nc ! type + cdata_all(10,iout)=r10 ! elevation of observation ! 10-m wind + cdata_all(11,iout)=qcm ! quality mark + cdata_all(12,iout)=obserr ! original obs error + cdata_all(13,iout)=usage ! usage parameter + cdata_all(14,iout)=idomsfc ! dominate surface type + cdata_all(15,iout)=tsavg ! skin temperature + cdata_all(16,iout)=ff10 ! 10 meter wind factor + cdata_all(17,iout)=sfcr ! surface roughness + cdata_all(18,iout)=dlon_earth_deg ! earth relative longitude (degree) + cdata_all(19,iout)=dlat_earth_deg ! earth relative latitude (degree) + cdata_all(20,iout)=gob ! station elevation (m) + cdata_all(21,iout)=zz ! terrain height at ob location + cdata_all(22,iout)=r_prvstg(1,1) ! provider name + cdata_all(23,iout)=r_sprvstg(1,1) ! subprovider name + endif + + end do loop_readsb2 + end do loop_msg2 + +! Close unit to bufr file + call closbf(lunin) + +! Write header record and data to output file for further processing + write(6,*) "READ_GNSSRSPD: nreal=",nreal," ndata=",ndata + allocate(cdata_out(nreal,ndata)) + do i=1,ndata + do k=1,nreal + cdata_out(k,i)=cdata_all(k,i) + end do + end do + deallocate(cdata_all) + + call count_obs(ndata,nreal,ilat,ilon,cdata_out,nobs) + write(lunout) obstype,sis,nreal,nchanl,ilat,ilon + write(lunout) cdata_out + deallocate(cdata_out) +900 continue + if(diagnostic_reg .and. ntest>0) write(6,*)'READ_GNSSRSPD: ',& + 'ntest, disterrmax=', ntest,disterrmax + if(diagnostic_reg .and. nvtest>0) write(6,*)'READ_GNSSRSPD: ',& + 'nvtest,vdisterrmax=',ntest,vdisterrmax + + if (ndata == 0) then + call closbf(lunin) + write(6,*)'READ_GNSSRSPD: no data to process' + endif + write(6,*)'READ_GNSSRSPD: nreal=',nreal + write(6,*)'READ_GNSSRSPD: ntb,nread,ndata,nodata=',ntb,nread,ndata,nodata + + call closbf(lunin) + close(lunin) + +! End of routine + return + +end subroutine read_gnssrspd + diff --git a/src/gsi/read_obs.F90 b/src/gsi/read_obs.F90 index aa0f11b4e3..6f0e32a4b6 100644 --- a/src/gsi/read_obs.F90 +++ b/src/gsi/read_obs.F90 @@ -141,6 +141,9 @@ subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread) ! 2019-08-21 H. Shao - add METOPC-C, COSMIC-2 and PAZ to the GPS check list ! 2020-05-21 H. Shao - add commercial GNSSRO (Spire, PlanetIQ, GeoOptics) and other existing missions to the check list ! 2021-02-20 X.Li - add viirs-m and get_hsst +! 2022-03-12 K. Apodaca - Enable CYGNSS ocean wind speed reading capability +! 2023-03-12 K. Apodaca - Enable GNSS-R L2 ocean wind speed reading (CYGNSS, Spire) +! 2023-03-12 K. Apodaca - Enable GNSS-R DDM reading capability ! ! input argument list: ! lexist - file status @@ -487,6 +490,16 @@ subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread) endif nread = nread + 1 end do loop_hdob + else if(trim(filename) == 'gnssrwndbufr')then + lexist = .false. + gnssrwndloop: do while(ireadmg(lnbufr,subset,idate2) >= 0) + if (trim(dtype)=='gnssrspd') then + lexist = .true. + exit gnssrwndloop + endif + nread = nread + 1 + end do gnssrwndloop + else if(trim(dtype) == 'pm2_5')then if (oneobtest_chem .and. oneob_type_chem=='pm2_5') then lexist=.true. @@ -901,7 +914,7 @@ subroutine read_obs(ndata,mype) if (obstype == 't' .or. obstype == 'uv' .or. & obstype == 'q' .or. obstype == 'ps' .or. & obstype == 'pw' .or. obstype == 'spd'.or. & - obstype == 'sst'.or. & + obstype == 'sst'.or. obstype == 'gnssrspd'.or. & obstype == 'tcp'.or. obstype == "lag".or. & obstype == 'dw' .or. obstype == 'rw' .or. & obstype == 'mta_cld' .or. obstype == 'gos_ctp' .or. & @@ -1449,8 +1462,13 @@ subroutine read_obs(ndata,mype) call read_prepbufr(nread,npuse,nouse,infile,obstype,lunout,twind,sis,& prsl_full,nobs_sub1(1,i),read_rec(i)) string='READ_PREPBUFR' - endif + + else if (obstype == 'gnssrspd' .and. index(infile,'gnssrwndbufr') /=0 ) then + call read_gnssrspd(nread,npuse,nouse,infile,obstype,lunout,twind,sis, & + nobs_sub1(1,i)) + string='READ_GNSSRSPD' + else if(obstype == 'howv') then if ( index(infile,'satmar') /=0) then diff --git a/src/gsi/setupgnssrspd.f90 b/src/gsi/setupgnssrspd.f90 new file mode 100644 index 0000000000..e90b94a1bd --- /dev/null +++ b/src/gsi/setupgnssrspd.f90 @@ -0,0 +1,940 @@ +module gnssrspd_setup + implicit none + private + public:: setup + interface setup; module procedure setupgnssrspd; end interface + +contains +subroutine setupgnssrspd(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsave) +!$$$ subprogram documentation block +! . . . . +! subprogram: setupgnssrspd compute rhs of oi for wind speed obs +! prgmmr: Karina Apodaca org: Spire Global, Inc. date: 2023-03-15 +! Based on: setup10mspd.f90 by: +! prgmmr: parrish org: np22 date: 1990-10-06 +! +! abstract: For wind speed observations, this routine +! a) reads obs assigned to given mpi task (geographic region), +! b) simulates obs from guess, +! c) apply some quality control to obs, +! d) load weight and innovation arrays used in minimization +! e) collects statistics for runtime diagnostic output +! f) writes additional diagnostic information to output file +! +! program history log: +! 2023-03-15 K. Apodaca - add GNSS-R L2 ocean wind speed +! +! input argument list: +! lunin - unit from which to read observations +! mype - mpi task id +! nele - number of data elements per observation +! nobs - number of observations +! +! output argument list: +! bwork - array containing information about obs-ges statistics +! awork - array containing information for data counts and gross checks +! +! attributes: +! language: f90 +! machine: ibm RS/6000 SP +! +!$$$ + use mpeu_util, only: die,perr,getindex + use kinds, only: r_kind,r_single,r_double,i_kind + use m_obsdiagNode, only: obs_diag + use m_obsdiagNode, only: obs_diags + use m_obsdiagNode, only: obsdiagLList_nextNode + use m_obsdiagNode, only: obsdiagNode_set + use m_obsdiagNode, only: obsdiagNode_get + use m_obsdiagNode, only: obsdiagNode_assert + + use obsmod, only: rmiss_single,& + lobsdiagsave,nobskeep,lobsdiag_allocated,time_offset,& + lobsdiag_forenkf + use obsmod, only: netcdf_diag, binary_diag, dirname, ianldate + use nc_diag_write_mod, only: nc_diag_init, nc_diag_header, nc_diag_metadata, & + nc_diag_write, nc_diag_data2d + use nc_diag_read_mod, only: nc_diag_read_init, nc_diag_read_get_dim, nc_diag_read_close + use m_obsNode, only: obsNode + use m_gnssrspdNode, only: gnssrspdNode + use m_gnssrspdNode, only: gnssrspdNode_appendto + use m_obsLList, only: obsLList + use obsmod, only: luse_obsdiag + use gsi_4dvar, only: nobs_bins,hr_obsbin + use guess_grids, only: nfldsig,hrdifsig,ges_lnprsl, & + comp_fact10,sfcmod_gfs,sfcmod_mm5 + use guess_grids, only: geop_hgtl + use gridmod, only: nsig,get_ij,twodvar_regional + use qcmod, only: npres_print,ptop,pbot + use constants, only: one,grav,rd,zero,four,tiny_r_kind, & + half,two,cg_term,huge_single,r1000,wgtlim,ten + use jfunc, only: jiter,last,miter,jiterstart + use state_vectors, only: svars3d, levels + use qcmod, only: dfact,dfact1 + use convinfo, only: nconvtype,cermin,cermax,cgross,cvar_b,cvar_pg,ictype + use convinfo, only: icsubtype + use gsi_bundlemod, only : gsi_bundlegetpointer + use gsi_metguess_mod, only : gsi_metguess_get,gsi_metguess_bundle + use m_dtime, only: dtime_setup, dtime_check + use sparsearr, only: sparr2, new, size, writearray, fullarray + + ! The following variables are the coefficients that describe the + ! linear regression fits that are used to define the dynamic + ! observation error (DOE) specifications for all reconnissance + ! observations collected within hurricanes/tropical cyclones; these + ! apply only to the regional forecast models (e.g., HWRF); Henry + ! R. Winterbottom (henry.winterbottom@noaa.gov). + + + implicit none + +! Declare passed variables + type(obsLList ),target,dimension(:),intent(in):: obsLL + type(obs_diags),target,dimension(:),intent(in):: odiagLL + + logical ,intent(in ) :: conv_diagsave + integer(i_kind) ,intent(in ) :: lunin,mype,nele,nobs + real(r_kind),dimension(100+7*nsig) ,intent(inout) :: awork + real(r_kind),dimension(npres_print,nconvtype,5,3),intent(inout) :: bwork + integer(i_kind) ,intent(in ) :: is ! ndat index + +! Declare local variables + character(len=*),parameter:: myname='setupgnssrspd' + +! Declare external calls for code analysis + external:: tintrp2a1,tintrp2a11 + external:: tintrp31 + external:: grdcrd1 + external:: stop2 + +! Declare local variables + + real(r_double) rstation_id + + real(r_kind) uob,vob,gnssrspdges,gnssrspdob,gnssrspdob0,goverrd,ratio_errors + real(r_kind) presw,factw,dpres,ugesin,vgesin,sfcr,skint + real(r_kind) scale + real(r_kind) val2,ressw,ress,error,ddiff,dx10,rhgh,prsfc,r0_001 + real(r_kind) sfcchk,prsln2,rwgt,tfact + real(r_kind) thirty,rsig,ratio,residual,obserrlm,obserror + real(r_kind) val,valqc,psges,drpx,dlat,dlon,dtime,rlow + real(r_kind) cg_gnssrspd,wgross,wnotgross,wgt,arg,exp_arg,term,rat_err2 + real(r_kind) errinv_input,errinv_adjst,errinv_final + real(r_kind) err_input,err_adjst,err_final + real(r_kind),dimension(nele,nobs):: data + real(r_kind),dimension(nobs):: dup + real(r_kind),dimension(nsig)::prsltmp,tges + real(r_single),allocatable,dimension(:,:)::rdiagbuf + + integer(i_kind) mm1,ibin,ioff,ioff0 + integer(i_kind) ii,jj,i,nchar,nreal,k,j,l,nty,nn,ikxx + integer(i_kind) ier,ilon,ilat,ipres,iuob,ivob,id,itime,ikx + integer(i_kind) ihgt,iqc,ier2,iuse,ilate,ilone,istnelv,izz,iprvd,isprvd + integer(i_kind) idomsfc,iskint,iff10,isfcr,isli + + type(sparr2) :: dhx_dx + integer(i_kind) :: iz, u_ind, v_ind, nnz, nind + real(r_kind) :: delz + + logical,dimension(nobs):: luse,muse + integer(i_kind),dimension(nobs):: ioid ! initial (pre-distribution) obs ID + logical proceed + + character(8) station_id + character(8),allocatable,dimension(:):: cdiagbuf + character(8),allocatable,dimension(:):: cprvstg,csprvstg + character(8) c_prvstg,c_sprvstg + real(r_double) r_prvstg,r_sprvstg + + logical:: in_curbin, in_anybin, save_jacobian + type(gnssrspdNode),pointer:: my_head + type(obs_diag),pointer:: my_diag + type(obs_diags),pointer:: my_diagLL + + logical z_height + real(r_kind) zsges,dstn + real(r_kind),dimension(nsig):: zges + real(r_kind) dz,zob,z1,z2,p1,p2,dz21,dlnp21,pobl + integer(i_kind) k1,k2 + + equivalence(rstation_id,station_id) + equivalence(r_prvstg,c_prvstg) + equivalence(r_sprvstg,c_sprvstg) + + real(r_kind),allocatable,dimension(:,:,: ) :: ges_ps + real(r_kind),allocatable,dimension(:,:,: ) :: ges_z + real(r_kind),allocatable,dimension(:,:,:,:) :: ges_u + real(r_kind),allocatable,dimension(:,:,:,:) :: ges_v + real(r_kind),allocatable,dimension(:,:,:,:) :: ges_tv + + type(obsLList),pointer,dimension(:):: gnssrspdhead + gnssrspdhead => obsLL(:) + +!****************************************************************************** +! Read and reformat observations in work arrays. + read(lunin)data,luse,ioid + +! index information for data array (see reading routine) + ier=1 ! index of obs error + ilon=2 ! index of grid relative obs location (x) + ilat=3 ! index of grid relative obs location (y) + ipres=4 ! index of pressure + iuob=5 ! index of u observation + ivob=6 ! index of v observation + id=7 ! index of station id + itime=8 ! index of observation time in data array + ikxx=9 ! index of ob type + ihgt=10 ! index of observation elevation + iqc=11 ! index of quality mark + ier2=12 ! index of original-original obs error ratio + iuse=13 ! index of use parameter + idomsfc=14 ! index of dominant surface type + iskint=15 ! index of surface skin temperature + iff10=16 ! index of 10 meter wind factor + isfcr=17 ! index of surface roughness + ilone=18 ! index of longitude (degrees) + ilate=19 ! index of latitude (degrees) + istnelv=20 ! index of station elevation (m) + izz=21 ! index of surface height + iprvd=22 ! index of observation provider + isprvd=23 ! index of observation subprovider + + mm1=mype+1 + scale=one + rsig=nsig + thirty = 30.0_r_kind + r0_001=0.001_r_kind + goverrd=grav/rd + + save_jacobian = conv_diagsave .and. jiter==jiterstart .and. lobsdiag_forenkf + +! Check to see if required guess fields are available + call check_vars_(proceed) + if(.not.proceed) return ! not all vars available, simply return + +! If require guess vars available, extract from bundle ... + call init_vars_ + +! If requested, save select data for output to diagnostic file + if(conv_diagsave)then + ii=0 + nchar=1 + ioff0=21 + nreal=ioff0 + if (lobsdiagsave) nreal=nreal+4*miter+1 + if (twodvar_regional) then; nreal=nreal+2; allocate(cprvstg(nobs),csprvstg(nobs)); endif + if (save_jacobian) then + nnz = 4 ! number of non-zero elements in dH(x)/dx profile + nind = 2 + call new(dhx_dx, nnz, nind) + nreal = nreal + size(dhx_dx) + endif + allocate(cdiagbuf(nobs),rdiagbuf(nreal,nobs)) + if(netcdf_diag) call init_netcdf_diag_ + end if + + !!awork(:) = zero + + do i=1,nobs + muse(i)=nint(data(iuse,i)) <= jiter + end do + + dup=one + do k=1,nobs + do l=k+1,nobs + if(data(ilat,k) == data(ilat,l) .and. & + data(ilon,k) == data(ilon,l) .and. & + data(ipres,k)== data(ipres,l) .and. & + data(ier,k) < r1000 .and. data(ier,l) < r1000 .and. & + muse(l) .and. muse(k))then + + tfact=min(one,abs(data(itime,k)-data(itime,l))/dfact1) + dup(k)=dup(k)+one-tfact*tfact*(one-dfact) + dup(l)=dup(l)+one-tfact*tfact*(one-dfact) + end if + end do + end do + + call dtime_setup() + do i=1,nobs + dtime=data(itime,i) + call dtime_check(dtime, in_curbin, in_anybin) + if(.not.in_anybin) cycle + + if(in_curbin) then + dlat=data(ilat,i) + dlon=data(ilon,i) + + dpres=data(ipres,i) + error=data(ier2,i) + ikx=nint(data(ikxx,i)) + endif + +! Link observation to appropriate observation bin + if (nobs_bins>1) then + ibin = NINT( dtime/hr_obsbin ) + 1 + else + ibin = 1 + endif + IF (ibin<1.OR.ibin>nobs_bins) write(6,*)mype,'Error nobs_bins,ibin= ',nobs_bins,ibin + + if (luse_obsdiag) my_diagLL => odiagLL(ibin) + +! Link obs to diagnostics structure + if (luse_obsdiag) then + my_diag => obsdiagLList_nextNode(my_diagLL ,& + create = .not.lobsdiag_allocated ,& + idv = is ,& + iob = ioid(i) ,& + ich = 1 ,& + elat = data(ilate,i) ,& + elon = data(ilone,i) ,& + luse = luse(i) ,& + miter = miter ) + + if(.not.associated(my_diag)) call die(myname, & + 'obsdiagLList_nextNode(), create =', .not.lobsdiag_allocated) + endif + + if(.not.in_curbin) cycle + +! Load obs error and u,v obs + obserror = max(cermin(ikx),min(cermax(ikx),data(ier,i))) + uob = data(iuob,i) + vob = data(ivob,i) + + + gnssrspdob=sqrt(uob*uob+vob*vob) + call tintrp2a1(ges_tv,tges,dlat,dlon,dtime,hrdifsig,& + nsig,mype,nfldsig) + call tintrp2a11(ges_ps,psges,dlat,dlon,dtime,hrdifsig,& + mype,nfldsig) + call tintrp2a1(ges_lnprsl,prsltmp,dlat,dlon,dtime,hrdifsig,& + nsig,mype,nfldsig) + + factw = data(iff10,i) + if(sfcmod_gfs .or. sfcmod_mm5)then + sfcr = data(isfcr,i) + skint = data(iskint,i) + isli=data(idomsfc,i) + call comp_fact10(dlat,dlon,dtime,skint,sfcr,isli,mype,factw) + end if + + nty=ictype(ikx) + +! Initialize logical + z_height = .false. + +! Process observations reported with height differently than those +! reported with pressure. Type 260=nacelle 261=tower wind spd are +! encoded in NCEP prepbufr files with geometric height above +! sea level. + +! Process cygnss observations at mslp + if ( nty == 600 ) then + z_height = .true. + data(ihgt,i) = ten + endif + + if (z_height) then + + drpx = zero + dpres = data(ihgt,i) + dstn = data(istnelv,i) + call tintrp2a11(ges_z,zsges,dlat,dlon,dtime,hrdifsig,& + mype,nfldsig) + +! Get guess surface elevation and geopotential height profile +! at observation location. + call tintrp2a1(geop_hgtl,zges,dlat,dlon,dtime,hrdifsig,& + nsig,mype,nfldsig) + +! Convert observation height (in dpres) from meters to grid relative +! units. Save the observation height in zob for later use. + zob = dpres + call grdcrd1(dpres,zges,nsig,1) + factw=one + rlow=zero + rhgh=zero + +! Compute observation pressure (only used for diagnostics) +! Set indices of model levels below (k1) and above (k2) observation. + if (dpresnsig) then + z1=zges(nsig-1); p1=prsltmp(nsig-1) + z2=zges(nsig); p2=prsltmp(nsig) + drpx = 1.e6_r_kind + else + k=dpres + k1=min(max(1,k),nsig) + k2=max(1,min(k+1,nsig)) + z1=zges(k1); p1=prsltmp(k1) + z2=zges(k2); p2=prsltmp(k2) + endif + + dz21 = z2-z1 + dlnp21 = p2-p1 + dz = zob-z1 + pobl = p1 + (dlnp21/dz21)*dz + presw = ten*exp(pobl) + + prsfc=psges + prsln2=log(exp(prsltmp(1))/prsfc) + sfcchk=log(psges) + if(dpres <= prsln2)then + factw=one + else + dx10=-goverrd*ten/tges(1) + if (dpres < dx10)then + factw=(dpres-dx10+factw*(prsln2-dpres))/(prsln2-dx10) + end if + end if + +! Put obs pressure in correct units to get grid coord. number + dpres=log(exp(dpres)*prsfc) + call grdcrd1(dpres,prsltmp(1),nsig,-1) + +! Get approx k value of sfc by using surface pressure of 1st ob + call grdcrd1(sfcchk,prsltmp(1),nsig,-1) + + +! Check to see if observations is below what is seen to be the surface + rlow=max(sfcchk-dpres,zero) + + rhgh=max(dpres-r0_001-rsig,zero) + + endif ! end of process observations with reported pressure + + if(luse(i))then + awork(1) = awork(1) + one + if(rlow/=zero) awork(2) = awork(2) + one + if(rhgh/=zero) awork(3) = awork(3) + one + end if + + ratio_errors=error/(data(ier,i)+drpx+1.0e6_r_kind*rhgh+four*rlow) + + +! Interpolate guess u and v to observation location and time. + call tintrp31(ges_u,ugesin,dlat,dlon,dpres,dtime, & + hrdifsig,mype,nfldsig) + call tintrp31(ges_v,vgesin,dlat,dlon,dpres,dtime, & + hrdifsig,mype,nfldsig) + + +! Apply 10-meter wind reduction factor to guess winds. Compute +! guess wind speed. + ugesin=factw*ugesin + vgesin=factw*vgesin + gnssrspdges=sqrt(ugesin*ugesin+vgesin*vgesin) + + iz = max(1, min( int(dpres), nsig)) + delz = max(zero, min(dpres - float(iz), one)) + + if (save_jacobian) then + + u_ind = getindex(svars3d, 'u') + if (u_ind < 0) then + print *, 'Error: no variable u in state vector. Exiting.' + call stop2(1300) + endif + v_ind = getindex(svars3d, 'v') + if (v_ind < 0) then + print *, 'Error: no variable v in state vector. Exiting.' + call stop2(1300) + endif + + dhx_dx%st_ind(1) = iz + sum(levels(1:u_ind-1)) + dhx_dx%end_ind(1) = min(iz + 1,nsig) + sum(levels(1:u_ind-1)) + + dhx_dx%val(1) = (one - delz) * two * ugesin + dhx_dx%val(2) = delz * two * ugesin + + dhx_dx%st_ind(2) = iz + sum(levels(1:v_ind-1)) + dhx_dx%end_ind(2) = min(iz + 1,nsig) + sum(levels(1:v_ind-1)) + + dhx_dx%val(3) = (one - delz) * two * vgesin + dhx_dx%val(4) = delz * two * ugesin + endif + + + ddiff = gnssrspdob-gnssrspdges + +! write(6,*) "cygnss gnssrspdob=",gnssrspdob," gnssrspdges=",gnssrspdges," ddiff=",ddiff + + error=one/error + +! Check to see if observations is above the top of the model (regional mode) + if (dpres>rsig) ratio_errors=zero + + +! Gross error checks + obserror = one/max(ratio_errors*error,tiny_r_kind) + obserrlm = max(cermin(ikx),min(cermax(ikx),obserror)) + residual = abs(ddiff) + ratio = residual/obserrlm + if (ratio>cgross(ikx) .or. ratio_errors < tiny_r_kind) then + if (luse(i)) awork(4) = awork(4)+one + error = zero + ratio_errors = zero + muse(i)=.false. + else + ratio_errors=ratio_errors/sqrt(dup(i)) + end if + + if (ratio_errors*error <=tiny_r_kind) muse(i)=.false. + if (nobskeep>0.and.luse_obsdiag) call obsdiagNode_get(my_diag, jiter=nobskeep, muse=muse(i)) + +! Compute penalty terms (linear & nonlinear qc). + val = error*ddiff + val2 = val*val + exp_arg = -half*val2 + rat_err2 = ratio_errors**2 + if (cvar_pg(ikx) > tiny_r_kind .and. error > tiny_r_kind) then + arg = exp(exp_arg) + wnotgross= one-cvar_pg(ikx) + cg_gnssrspd=cvar_b(ikx) + wgross = cg_term*cvar_pg(ikx)/(cg_gnssrspd*wnotgross) + term = log((arg+wgross)/(one+wgross)) + wgt = one-wgross/(arg+wgross) + wgt = wgt/wgtlim + else + term = exp_arg + wgt = wgtlim + rwgt = wgt/wgtlim + endif + valqc = -two*rat_err2*term + + +! Accumulate statistics for obs belonging to this task + if (luse(i) .and. muse(i)) then + if(rwgt < one) awork(61) = awork(61)+one + awork(5)=awork(5) + val2*rat_err2 + awork(6)=awork(6) + one + awork(22)=awork(22) + valqc + end if + +! Loop over pressure level groupings and obs to accumulate statistics +! as a function of observation type. + do k = 1,npres_print + if(luse(i) .and.presw >ptop(k) .and. presw<=pbot(k))then + ress = scale*ddiff + ressw = ress*ress + nn=1 + if (.not. muse(i)) then + nn=2 + if(ratio_errors*error >=tiny_r_kind)nn=3 + end if + bwork(k,ikx,1,nn) = bwork(k,ikx,1,nn)+one ! count + bwork(k,ikx,2,nn) = bwork(k,ikx,2,nn)+ddiff ! bias + bwork(k,ikx,3,nn) = bwork(k,ikx,3,nn)+ressw ! (o-g)**2 + bwork(k,ikx,4,nn) = bwork(k,ikx,4,nn)+val2*rat_err2 ! penalty + bwork(k,ikx,5,nn) = bwork(k,ikx,5,nn)+valqc ! nonlin qc penalty + + end if + end do + + !write(6,*) "cygnss: i,luse(i),muse(i)=",i,luse(i),muse(i) + + if (luse_obsdiag) then + call obsdiagNode_set(my_diag, wgtjo=(error*ratio_errors)**2, & + jiter=jiter, muse=muse(i), nldepart=gnssrspdob-sqrt(ugesin*ugesin+vgesin*vgesin)) + endif + +! If obs is "acceptable", load array with obs info for use +! in inner loop minimization (int* and stp* routines) + if (.not. last .and. muse(i)) then + + allocate(my_head) + call gnssrspdNode_appendto(my_head,gnssrspdhead(ibin)) + + my_head%idv = is + my_head%iob = ioid(i) + my_head%elat= data(ilate,i) + my_head%elon= data(ilone,i) + +! Set (i,j) indices of guess gridpoint that bound obs location + call get_ij(mm1,dlat,dlon,my_head%ij,my_head%wij) + + my_head%factw=factw + do j=1,4 + my_head%wij(j)=factw*my_head%wij(j) + end do + my_head%raterr2= ratio_errors**2 + my_head%res = gnssrspdob + my_head%uges = ugesin + my_head%vges = vgesin + my_head%err2 = error**2 + my_head%time = dtime + my_head%luse = luse(i) + my_head%b = cvar_b(ikx) + my_head%pg = cvar_pg(ikx) + + if (luse_obsdiag) then + call obsdiagNode_assert(my_diag, my_head%idv,my_head%iob,1,myname,'my_diag:my_head') + my_head%diags => my_diag + endif + + my_head => null() + end if +! Save select output for diagnostic file + if(conv_diagsave .and. luse(i))then + ii=ii+1 + gnssrspdob0 = sqrt(data(iuob,i)*data(iuob,i)+data(ivob,i)*data(ivob,i)) + rstation_id = data(id,i) + err_input = data(ier2,i) + err_adjst = data(ier,i) + if (ratio_errors*error>tiny_r_kind) then + err_final = one/(ratio_errors*error) + else + err_final = huge_single + endif + + errinv_input = huge_single + errinv_adjst = huge_single + errinv_final = huge_single + if (err_input>tiny_r_kind) errinv_input = one/err_input + if (err_adjst>tiny_r_kind) errinv_adjst = one/err_adjst + if (err_final>tiny_r_kind) errinv_final = one/err_final + + if(binary_diag) call contents_binary_diag_(my_diag) + if(netcdf_diag) call contents_netcdf_diag_(my_diag) + + end if + +! End of loop over observations + end do + +! Release memory of local guess arrays + call final_vars_ + +! Write information to diagnostic file + if(conv_diagsave) then + if(netcdf_diag) call nc_diag_write + if(binary_diag .and. ii>0)then + write(7)'gnssrspd',nchar,nreal,ii,mype,ioff0 + write(7)cdiagbuf(1:ii),rdiagbuf(:,1:ii) + deallocate(cdiagbuf,rdiagbuf) + + if (twodvar_regional) then + write(7)cprvstg(1:ii),csprvstg(1:ii) + deallocate(cprvstg,csprvstg) + endif + end if + end if + +! End of routine + + return + contains + + subroutine check_vars_ (proceed) + logical,intent(inout) :: proceed + integer(i_kind) ivar, istatus +! Check to see if required guess fields are available + call gsi_metguess_get ('var::ps', ivar, istatus ) + proceed=ivar>0 + call gsi_metguess_get ('var::z' , ivar, istatus ) + proceed=proceed.and.ivar>0 + call gsi_metguess_get ('var::u' , ivar, istatus ) + proceed=proceed.and.ivar>0 + call gsi_metguess_get ('var::v' , ivar, istatus ) + proceed=proceed.and.ivar>0 + call gsi_metguess_get ('var::tv', ivar, istatus ) + proceed=proceed.and.ivar>0 + end subroutine check_vars_ + + subroutine init_vars_ + + real(r_kind),dimension(:,: ),pointer:: rank2=>NULL() + real(r_kind),dimension(:,:,:),pointer:: rank3=>NULL() + character(len=5) :: varname + integer(i_kind) ifld, istatus + +! If require guess vars available, extract from bundle ... + if(size(gsi_metguess_bundle)==nfldsig) then +! get ps ... + varname='ps' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus) + if (istatus==0) then + if(allocated(ges_ps))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif + allocate(ges_ps(size(rank2,1),size(rank2,2),nfldsig)) + ges_ps(:,:,1)=rank2 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank2,istatus) + ges_ps(:,:,ifld)=rank2 + enddo + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + endif +! get z ... + varname='z' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus) + if (istatus==0) then + if(allocated(ges_z))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif + allocate(ges_z(size(rank2,1),size(rank2,2),nfldsig)) + ges_z(:,:,1)=rank2 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank2,istatus) + ges_z(:,:,ifld)=rank2 + enddo + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + endif +! get u ... + varname='u' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) + if (istatus==0) then + if(allocated(ges_u))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif + allocate(ges_u(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) + ges_u(:,:,:,1)=rank3 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) + ges_u(:,:,:,ifld)=rank3 + enddo + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + endif +! get v ... + varname='v' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) + if (istatus==0) then + if(allocated(ges_v))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif + allocate(ges_v(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) + ges_v(:,:,:,1)=rank3 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) + ges_v(:,:,:,ifld)=rank3 + enddo + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + endif +! get tv ... + varname='tv' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) + if (istatus==0) then + if(allocated(ges_tv))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif + allocate(ges_tv(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) + ges_tv(:,:,:,1)=rank3 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) + ges_tv(:,:,:,ifld)=rank3 + enddo + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + endif + else + write(6,*) trim(myname), ': inconsistent vector sizes (nfldsig,size(metguess_bundle) ',& + nfldsig,size(gsi_metguess_bundle) + call stop2(999) + endif + end subroutine init_vars_ + + subroutine init_netcdf_diag_ + character(len=80) string + character(len=128) diag_conv_file + integer(i_kind) ncd_fileid,ncd_nobs + logical append_diag + logical,parameter::verbose=.false. + write(string,900) jiter +900 format('conv_gnssrspd_',i2.2,'.nc4') + diag_conv_file=trim(dirname) // trim(string) + + inquire(file=diag_conv_file, exist=append_diag) + + if (append_diag) then + call nc_diag_read_init(diag_conv_file,ncd_fileid) + ncd_nobs = nc_diag_read_get_dim(ncd_fileid,'nobs') + call nc_diag_read_close(diag_conv_file) + + if (ncd_nobs > 0) then + if(verbose) print *,'file ' // trim(diag_conv_file) // ' exists. Appending. nobs,mype=',ncd_nobs,mype + else + if(verbose) print *,'file ' // trim(diag_conv_file) // ' exists but contains no obs. Not appending. nobs,mype=',ncd_nobs,mype + append_diag = .false. ! if there are no obs in existing file, then do not try to append + endif + end if + + call nc_diag_init(diag_conv_file, append=append_diag) + + if (.not. append_diag) then ! don't write headers on append - the module will break? + call nc_diag_header("date_time",ianldate ) + if (save_jacobian) then + call nc_diag_header("jac_nnz", nnz) + call nc_diag_header("jac_nind", nind) + endif + endif + end subroutine init_netcdf_diag_ + subroutine contents_binary_diag_(odiag) + type(obs_diag),pointer,intent(in):: odiag + cdiagbuf(ii) = station_id ! station id + + rdiagbuf(1,ii) = ictype(ikx) ! observation type + rdiagbuf(2,ii) = icsubtype(ikx) ! observation subtype + + rdiagbuf(3,ii) = data(ilate,i) ! observation latitude (degrees) + rdiagbuf(4,ii) = data(ilone,i) ! observation longitude (degrees) + rdiagbuf(5,ii) = data(istnelv,i) ! station elevation (meters) + rdiagbuf(6,ii) = presw ! observation pressure (hPa) + rdiagbuf(7,ii) = data(ihgt,i) ! observation height (meters) + rdiagbuf(8,ii) = dtime-time_offset ! obs time (hours relative to analysis time) + + rdiagbuf(9,ii) = data(iqc,i) ! input prepbufr qc or event mark + rdiagbuf(10,ii) = rmiss_single ! setup qc or event mark + rdiagbuf(11,ii) = data(iuse,i) ! read_prepbufr data usage flag + if(muse(i)) then + rdiagbuf(12,ii) = one ! analysis usage flag (1=use, -1=not used) + else + rdiagbuf(12,ii) = -one + endif + + rdiagbuf(13,ii) = rwgt ! nonlinear qc relative weight + rdiagbuf(14,ii) = errinv_input ! prepbufr inverse obs error (m/s)**-1 + rdiagbuf(15,ii) = errinv_adjst ! read_prepbufr inverse obs error (m/s)**-1 + rdiagbuf(16,ii) = errinv_final ! final inverse observation error (m/s)**-1 + + rdiagbuf(17,ii) = gnssrspdob ! wind speed observation (m/s) + rdiagbuf(18,ii) = ddiff ! obs-ges used in analysis (m/s) + rdiagbuf(19,ii) = gnssrspdob0-gnssrspdges ! obs-ges w/o bias correction (m/s) (future slot) + + rdiagbuf(20,ii) = factw ! 10m wind reduction factor + + rdiagbuf(21,ii) = 1.e+10_r_single ! ges ensemble spread (filled in by EnKF) + + ioff=ioff0 + if (lobsdiagsave) then + do jj=1,miter + ioff=ioff+1 + if (odiag%muse(jj)) then + rdiagbuf(ioff,ii) = one + else + rdiagbuf(ioff,ii) = -one + endif + enddo + do jj=1,miter+1 + ioff=ioff+1 + rdiagbuf(ioff,ii) = odiag%nldepart(jj) + enddo + do jj=1,miter + ioff=ioff+1 + rdiagbuf(ioff,ii) = odiag%tldepart(jj) + enddo + do jj=1,miter + ioff=ioff+1 + rdiagbuf(ioff,ii) = odiag%obssen(jj) + enddo + endif + + if (twodvar_regional) then + ioff = ioff + 1 + rdiagbuf(ioff,ii) = data(idomsfc,i) ! dominate surface type + ioff = ioff + 1 + rdiagbuf(ioff,ii) = data(izz,i) ! model terrain at ob location + r_prvstg = data(iprvd,i) + cprvstg(ii) = c_prvstg ! provider name + r_sprvstg = data(isprvd,i) + csprvstg(ii) = c_sprvstg ! subprovider name + endif + + if (save_jacobian) then + call writearray(dhx_dx, rdiagbuf(ioff+1:nreal,ii)) + ioff = ioff + size(dhx_dx) + endif + + end subroutine contents_binary_diag_ + subroutine contents_netcdf_diag_(odiag) + type(obs_diag),pointer,intent(in):: odiag +! Observation class + character(8),parameter :: obsclass = 'gnssrspd' + real(r_kind),dimension(miter) :: obsdiag_iuse + call nc_diag_metadata("Station_ID", station_id ) + call nc_diag_metadata("Observation_Class", obsclass ) + call nc_diag_metadata("Observation_Type", ictype(ikx) ) + call nc_diag_metadata("Observation_Subtype", icsubtype(ikx) ) + call nc_diag_metadata("Latitude", sngl(data(ilate,i)) ) + call nc_diag_metadata("Longitude", sngl(data(ilone,i)) ) + call nc_diag_metadata("Station_Elevation", sngl(data(istnelv,i)) ) + call nc_diag_metadata("Pressure", sngl(presw) ) + call nc_diag_metadata("Height", sngl(data(ihgt,i)) ) + call nc_diag_metadata("Time", sngl(dtime-time_offset)) + call nc_diag_metadata("Prep_QC_Mark", sngl(data(iqc,i)) ) + call nc_diag_metadata("Prep_Use_Flag", sngl(data(iuse,i)) ) +! call nc_diag_metadata("Nonlinear_QC_Var_Jb", var_jb ) + call nc_diag_metadata("Nonlinear_QC_Rel_Wgt", sngl(rwgt) ) + if(muse(i)) then + call nc_diag_metadata("Analysis_Use_Flag", sngl(one) ) + else + call nc_diag_metadata("Analysis_Use_Flag", sngl(-one) ) + endif + + call nc_diag_metadata("Errinv_Input", sngl(errinv_input) ) + call nc_diag_metadata("Errinv_Adjust", sngl(errinv_adjst) ) + call nc_diag_metadata("Errinv_Final", sngl(errinv_final) ) + + call nc_diag_metadata("Observation", sngl(gnssrspdob) ) + call nc_diag_metadata("Obs_Minus_Forecast_adjusted", sngl(ddiff) ) + call nc_diag_metadata("Obs_Minus_Forecast_unadjusted", sngl(gnssrspdob0-gnssrspdges) ) + + if (lobsdiagsave) then + do jj=1,miter + if (odiag%muse(jj)) then + obsdiag_iuse(jj) = one + else + obsdiag_iuse(jj) = -one + endif + enddo + + call nc_diag_data2d("ObsDiagSave_iuse", obsdiag_iuse ) + call nc_diag_data2d("ObsDiagSave_nldepart", odiag%nldepart ) + call nc_diag_data2d("ObsDiagSave_tldepart", odiag%tldepart ) + call nc_diag_data2d("ObsDiagSave_obssen", odiag%obssen ) + endif + + if (twodvar_regional) then + call nc_diag_metadata("Dominant_Sfc_Type", data(idomsfc,i) ) + call nc_diag_metadata("Model_Terrain", data(izz,i) ) + r_prvstg = data(iprvd,i) + call nc_diag_metadata("Provider_Name", c_prvstg ) + r_sprvstg = data(isprvd,i) + call nc_diag_metadata("Subprovider_Name", c_sprvstg ) + endif + if (save_jacobian) then + call nc_diag_data2d("Observation_Operator_Jacobian_stind", dhx_dx%st_ind) + call nc_diag_data2d("Observation_Operator_Jacobian_endind", dhx_dx%end_ind) + call nc_diag_data2d("Observation_Operator_Jacobian_val", real(dhx_dx%val,r_single)) + endif + + + end subroutine contents_netcdf_diag_ + + subroutine final_vars_ + if(allocated(ges_tv)) deallocate(ges_tv) + if(allocated(ges_v )) deallocate(ges_v ) + if(allocated(ges_u )) deallocate(ges_u ) + if(allocated(ges_z )) deallocate(ges_z ) + if(allocated(ges_ps)) deallocate(ges_ps) + end subroutine final_vars_ + +end subroutine setupgnssrspd +end module gnssrspd_setup diff --git a/src/gsi/setuprhsall.f90 b/src/gsi/setuprhsall.f90 index 8075956431..6c00de5803 100644 --- a/src/gsi/setuprhsall.f90 +++ b/src/gsi/setuprhsall.f90 @@ -101,6 +101,7 @@ subroutine setuprhsall(ndata,mype,init_pass,last_pass) ! 2019-03-15 Ladwig - add option for cloud analysis in observer ! 2019-03-28 Ladwig - add metar cloud obs as pseudo water vapor in var analysis ! 2020-09-08 CAPS(G. Zhao) - add 'l_use_dbz_directDA' flag not to sort obsdiag +! 2023-03-20 K Apodaca - add GNSS-R L2 Ocean Wind Speed ! ! input argument list: ! ndata(*,1)- number of prefiles retained for further processing @@ -165,7 +166,7 @@ subroutine setuprhsall(ndata,mype,init_pass,last_pass) use m_rhs, only: toss_gps_sub => rhs_toss_gps use m_rhs, only: i_ps,i_uv,i_t,i_q,i_pw,i_rw,i_dw,i_gps,i_sst,i_tcp,i_lag, & - i_gust,i_vis,i_pblh,i_wspd10m,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv, & + i_gust,i_vis,i_pblh,i_wspd10m,i_gnssrspd,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv, & i_tcamt,i_lcbas,i_cldch,i_uwnd10m,i_vwnd10m,i_swcp,i_lwcp use m_rhs, only: i_dbz use m_rhs, only: i_fed @@ -625,7 +626,7 @@ subroutine setuprhsall(ndata,mype,init_pass,last_pass) ! Compute and print statistics for "conventional" data call statsconv(mype,& i_ps,i_uv,i_t,i_q,i_pw,i_rw,i_dw,i_gps,i_sst,i_tcp,i_lag, & - i_gust,i_vis,i_pblh,i_wspd10m,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv, & + i_gust,i_vis,i_pblh,i_wspd10m,i_gnssrspd,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv, & i_tcamt,i_lcbas,i_cldch,i_uwnd10m,i_vwnd10m,i_swcp,i_lwcp,i_fed,i_dbz, & size(awork1,2),bwork1,awork1,ndata) diff --git a/src/gsi/statsconv.f90 b/src/gsi/statsconv.f90 index fc105515ff..b6bc849a37 100644 --- a/src/gsi/statsconv.f90 +++ b/src/gsi/statsconv.f90 @@ -1,6 +1,6 @@ subroutine statsconv(mype,& i_ps,i_uv,i_t,i_q,i_pw,i_rw,i_dw,i_gps,i_sst,i_tcp,i_lag, & - i_gust,i_vis,i_pblh,i_wspd10m,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv, & + i_gust,i_vis,i_pblh,i_wspd10m,i_gnssrspd,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv, & i_tcamt,i_lcbas,i_cldch,i_uwnd10m,i_vwnd10m,& i_swcp,i_lwcp,i_fed,i_dbz,i_ref,bwork,awork,ndata) !$$$ subprogram documentation block @@ -44,6 +44,7 @@ subroutine statsconv(mype,& ! 2015-07-10 pondeca - add cldch ! 2016-05-05 pondeca - add uwnd10m, vwnd10m ! 2017-05-12 Y. Wang and X. Wang - add dbz, POC: xuguang.wang@ou.edu +! 2022-03-15 K. Apodaca - add GNSS-R L2 ocean wind speed ! ! input argument list: ! mype - mpi task number @@ -62,6 +63,7 @@ subroutine statsconv(mype,& ! i_vis - index in awork array holding vis info ! i_pblh - index in awork array holding pblh info ! i_wspd10m- index in awork array holding wspd10m info +! i_gnssrspd - index in awork array holding gnssrspd info ! i_td2m - index in awork array holding td2m info ! i_mxtm - index in awork array holding mxtm info ! i_mitm - index in awork array holding mitm info @@ -94,13 +96,13 @@ subroutine statsconv(mype,& use constants, only: zero,three,five use obsmod, only: iout_sst,iout_pw,iout_t,iout_rw,iout_dw,& iout_uv,iout_gps,iout_ps,iout_q,iout_tcp,iout_lag,& - iout_gust,iout_vis,iout_pblh,iout_wspd10m,iout_td2m,& + iout_gust,iout_vis,iout_pblh,iout_wspd10m,iout_gnssrspd,iout_td2m,& iout_mxtm,iout_mitm,iout_pmsl,iout_howv,iout_tcamt,iout_lcbas,iout_cldch,& iout_uwnd10m,iout_vwnd10m,& iout_fed,iout_dbz,iout_swcp,iout_lwcp,& mype_dw,mype_rw,mype_sst,mype_gps,mype_uv,mype_ps,& mype_t,mype_pw,mype_q,mype_tcp,ndat,dtype,mype_lag,mype_gust,& - mype_vis,mype_pblh,mype_wspd10m,mype_td2m,mype_mxtm,mype_mitm,& + mype_vis,mype_pblh,mype_wspd10m,mype_gnssrspd,mype_td2m,mype_mxtm,mype_mitm,& mype_pmsl,mype_howv,mype_tcamt,mype_lcbas,mype_cldch,mype_uwnd10m,mype_vwnd10m,& mype_fed,mype_dbz,mype_swcp,mype_lwcp use qcmod, only: npres_print,ptop,pbot,ptopq,pbotq @@ -112,7 +114,7 @@ subroutine statsconv(mype,& ! Declare passed variables integer(i_kind) ,intent(in ) :: mype,i_ps,i_uv,& i_t,i_q,i_pw,i_rw,i_dw,i_gps,i_sst,i_tcp,i_lag,i_gust,i_vis,i_pblh,& - i_wspd10m,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv,i_tcamt,i_lcbas,& + i_wspd10m,i_gnssrspd,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv,i_tcamt,i_lcbas,& i_cldch,i_uwnd10m,i_vwnd10m,i_swcp,i_lwcp,i_fed,i_dbz,i_ref real(r_kind),dimension(7*nsig+100,i_ref) ,intent(in ) :: awork real(r_kind),dimension(npres_print,nconvtype,5,3),intent(in ) :: bwork @@ -122,7 +124,7 @@ subroutine statsconv(mype,& character(100) mesage integer(i_kind) numgrspw,numsst,nsuperp,nump,nhitopo,ntoodif - integer(i_kind) numgrsq,numhgh,numgust,numvis,numpblh,numwspd10m,numuwnd10m,numvwnd10m + integer(i_kind) numgrsq,numhgh,numgust,numvis,numpblh,numwspd10m,numgnssrspd,numuwnd10m,numvwnd10m integer(i_kind) numtd2m,nummxtm,nummitm,numpmsl,numhowv,numtcamt,numlcbas,numcldch integer(i_kind) numgrsswcp,numgrslwcp integer(i_kind) ntot,numlow,k,numssm,i,j @@ -677,6 +679,49 @@ subroutine statsconv(mype,& end if end if +! Summary report for conventional gnssrspd +if (mype == mype_gnssrspd) then + nread = 0 + nkeep = 0 + do i = 1, ndat + if (dtype(i) == 'gnssrspd') then + nread = nread + ndata(i, 2) + nkeep = nkeep + ndata(i, 3) + end if + end do + if (nread > 0) then + if (first) then + open(iout_gnssrspd) + else + open(iout_gnssrspd, position = 'append') + end if + + numgnssrspd = nint(awork(5, i_gnssrspd)) + pw = zero + pw3 = zero + if (nkeep > 0) then + mesage = 'current fit of conventional gnssrspd data, ranges in m/s' + do j = 1, nconvtype + pflag(j) = trim(ioctype(j)) == 'gnssrspd' + end do + call dtast(bwork, 1, pbotall, ptopall, mesage, jiter, iout_gnssrspd, pflag) + + numgross = nint(awork(6, i_gnssrspd)) + numfailqc = nint(awork(21, i_gnssrspd)) + if (numgnssrspd > 0) then + pw = awork(4, i_gnssrspd) / numgnssrspd + pw3 = awork(22, i_gnssrspd) / numgnssrspd + end if + write(iout_gnssrspd, 925) 'gnssrspd', numgross, numfailqc + end if + write(iout_gnssrspd, 950) 'gnssrspd', jiter, nread, nkeep, numgnssrspd + write(iout_gnssrspd, 951) 'gnssrspd', awork(4, i_gnssrspd), awork(22, i_gnssrspd), pw, pw3 + + close(iout_gnssrspd) + end if +end if + + ! Summary report for conventional td2m if(mype==mype_td2m) then nread=0 diff --git a/src/gsi/stpgnssrspd.f90 b/src/gsi/stpgnssrspd.f90 new file mode 100644 index 0000000000..3909c40534 --- /dev/null +++ b/src/gsi/stpgnssrspd.f90 @@ -0,0 +1,175 @@ +module stpgnssrspdmod + + +!$$$ module documentation block +! . . . . +! module: stpgnssrspdmod module for stpgnssrspd and its tangent linear stpgnssrspd_tl +! prgmmr: K. Apodaca org: Spire Global, Inc. date: 2022-03-12 +! Largely based on other stp_* routines +! +! abstract: module for stpgnssrspd and its tangent linear stpgnssrspd_tl +! +! program history log: +! 2023-09-21 K. Apodaca - add documentation +! subroutine included: +! sub stpgnssrspd +! +! attributes: +! language: f90 +! machine: +! +!$$$ end documentation block + +implicit none + +PRIVATE +PUBLIC stpgnssrspd + +contains + +subroutine stpgnssrspd(gnssrspdhead,rval,sval,out,sges,nstep) +!$$$ subprogram documentation block +! . . . . +! subprogram: stpgnssrspd calculate penalty and stepsize terms +! for wind speed, with nonlinear qc. +! 2023-03-15 K. Apodaca - add GNSS-R L2 ocean wind speed +! +! abstract: calculate penalty and stepsize terms for wind speed +! +! program history log: + +! +! input argument list: +! gnssrspdhead +! ru - search direction for u +! rv - search direction for v +! su - analysis increment for u +! sv - analysis increment for v +! sges - step size estimates (nstep) +! nstep - number of stepsizes (==0 means use outer iteration values) +! +! output argument list +! out(1:nstep) - contribution to penalty from wind speed sges(1:nstep) +! +! attributes: +! language: f90 +! machine: ibm RS/6000 SP +! +!$$$ + use kinds, only: r_kind,i_kind,r_quad + use qcmod, only: nlnqc_iter,varqc_iter + use constants, only: zero,half,one,two,tiny_r_kind,cg_term,zero_quad,r3600 + use gsi_4dvar, only: ltlint + use gsi_bundlemod, only: gsi_bundle + use gsi_bundlemod, only: gsi_bundlegetpointer + use m_obsNode, only: obsNode + use m_gnssrspdNode, only: gnssrspdNode + use m_gnssrspdNode, only: gnssrspdNode_typecast + use m_gnssrspdNode, only: gnssrspdNode_nextcast + implicit none + +! Declare passed variables + class(obsNode), pointer ,intent(in ) :: gnssrspdhead + integer(i_kind) ,intent(in ) :: nstep + real(r_kind),dimension(max(1,nstep)),intent(in ) :: sges + real(r_quad),dimension(max(1,nstep)),intent(inout) :: out + type(gsi_bundle) ,intent(in ) :: rval,sval + +! Declare local variables + integer(i_kind) j1,j2,j3,j4,kk,ier,istatus + real(r_kind) w1,w2,w3,w4,time_gnssrspd + real(r_kind) valu,valv,ucur,vcur,gnssrspdnl,gnssrspdtl,uu,vv,gnssrspd + real(r_kind),dimension(max(1,nstep)):: pen + real(r_kind) cg_gnssrspd,pencur,wgross,wnotgross + real(r_kind) pg_gnssrspd,pentl + real(r_kind),pointer,dimension(:) :: su,sv + real(r_kind),pointer,dimension(:) :: ru,rv + type(gnssrspdNode), pointer :: gnssrspdptr + + out=zero_quad + +! If no gnssrspd data return + if(.not. associated(gnssrspdhead))return + + time_gnssrspd=zero +! Retrieve pointers +! Simply return if any pointer not found + ier=0 + call gsi_bundlegetpointer(sval,'u',su,istatus);ier=istatus+ier + call gsi_bundlegetpointer(sval,'v',sv,istatus);ier=istatus+ier + call gsi_bundlegetpointer(rval,'u',ru,istatus);ier=istatus+ier + call gsi_bundlegetpointer(rval,'v',rv,istatus);ier=istatus+ier + if(ier/=0)return + + gnssrspdptr => gnssrspdNode_typecast(gnssrspdhead) + do while (associated(gnssrspdptr)) + + if(gnssrspdptr%luse)then + if(nstep > 0)then + j1 = gnssrspdptr%ij(1) + j2 = gnssrspdptr%ij(2) + j3 = gnssrspdptr%ij(3) + j4 = gnssrspdptr%ij(4) + w1 = gnssrspdptr%wij(1) + w2 = gnssrspdptr%wij(2) + w3 = gnssrspdptr%wij(3) + w4 = gnssrspdptr%wij(4) + + valu=w1* ru(j1)+w2* ru(j2)+w3* ru(j3)+w4* ru(j4) + valv=w1* rv(j1)+w2* rv(j2)+w3* rv(j3)+w4* rv(j4) + ucur=w1* su(j1)+w2* su(j2)+w3* su(j3)+w4* su(j4)+gnssrspdptr%uges + vcur=w1* sv(j1)+w2* sv(j2)+w3* sv(j3)+w4* sv(j4)+gnssrspdptr%vges + + if (ltlint) then + gnssrspd=sqrt(ucur*ucur+vcur*vcur)-gnssrspdptr%res + pencur=gnssrspd*gnssrspd*gnssrspdptr%err2 + do kk=1,nstep + gnssrspdnl=sqrt(ucur*ucur+vcur*vcur) + gnssrspdtl=ucur*valu+vcur*valv + if (gnssrspdnl>tiny_r_kind*100._r_kind) then + gnssrspdtl=gnssrspdtl/gnssrspdnl + else + gnssrspdtl=zero + endif + pentl =two*gnssrspdtl*gnssrspd*gnssrspdptr%err2 + pen(kk)=pencur+sges(kk)*pentl + end do + else + do kk=1,nstep + uu=ucur+sges(kk)*valu + vv=vcur+sges(kk)*valv + gnssrspd=sqrt(uu*uu+vv*vv)-gnssrspdptr%res + pen(kk)=gnssrspd*gnssrspd*gnssrspdptr%err2 + end do + end if + else + pen(1)=gnssrspdptr%res*gnssrspdptr%res*gnssrspdptr%err2 + end if + +! Modify penalty term if nonlinear QC + if (nlnqc_iter .and. gnssrspdptr%pg > tiny_r_kind .and. & + gnssrspdptr%b > tiny_r_kind) then + pg_gnssrspd=gnssrspdptr%pg*varqc_iter + cg_gnssrspd=cg_term/gnssrspdptr%b + wnotgross= one-pg_gnssrspd + wgross = pg_gnssrspd*cg_gnssrspd/wnotgross + pencur = -two*log((exp(-half*pencur) + wgross)/(one+wgross)) + do kk=1,max(1,nstep) + pen(kk) = -two*log((exp(-half*pen(kk) ) + wgross)/(one+wgross)) + enddo + endif + + out(1) = out(1)+pen(1)*gnssrspdptr%raterr2 + do kk=2,nstep + out(kk) = out(kk)+(pen(kk)-pen(1))*gnssrspdptr%raterr2 + end do + + end if + + gnssrspdptr => gnssrspdNode_nextcast(gnssrspdptr) + + end do + return +end subroutine stpgnssrspd + +end module stpgnssrspdmod From a5e2a439ee8a8cd87571156d5b14f58af7f64bfb Mon Sep 17 00:00:00 2001 From: Jeff Whitaker Date: Tue, 9 Jul 2024 09:22:42 -0600 Subject: [PATCH 02/11] fix for cris-fsr memory corruption (#767) --- src/gsi/correlated_obsmod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gsi/correlated_obsmod.F90 b/src/gsi/correlated_obsmod.F90 index 17cd94efe1..9c8791b0f3 100644 --- a/src/gsi/correlated_obsmod.F90 +++ b/src/gsi/correlated_obsmod.F90 @@ -422,7 +422,7 @@ subroutine set_(instrument,fname,mask,method,kreq,kmut,ErrorCov) ErrorCov%nch_active = coun if (.not.GMAO_ObsErrorCov) ErrorCov%nctot = nctot call create_(coun,ErrorCov) - allocate(indxRf(nch_active),indxR(nch_active),Rcov(nctot,nctot)) + allocate(indxRf(coun),indxR(nch_active),Rcov(nctot,nctot)) ! Read GSI-like channel numbers used in estimating R for this instrument read(lu,IOSTAT=ioflag) indxR From 3e27bb8ada1c7202ea1eb99957f0b143464410c7 Mon Sep 17 00:00:00 2001 From: RussTreadon-NOAA <26926959+RussTreadon-NOAA@users.noreply.github.com> Date: Mon, 22 Jul 2024 10:53:55 -0400 Subject: [PATCH 03/11] Update global_4denvar and global_enkf ctests to reflect GFS v17 (#774) --- regression/global_enkf.sh | 3 --- regression/regression_namelists.sh | 8 ++++---- regression/regression_namelists_db.sh | 2 +- 3 files changed, 5 insertions(+), 8 deletions(-) diff --git a/regression/global_enkf.sh b/regression/global_enkf.sh index a35f8d109f..f5b4f3d1c6 100755 --- a/regression/global_enkf.sh +++ b/regression/global_enkf.sh @@ -76,14 +76,11 @@ export use_correlated_oberrs=${use_correlated_oberrs:-".false."} if [ $USE_CORRELATED_OBERRS == "YES" ]; then export use_correlated_oberrs=".true." fi -export imp_physics=${imp_physics:-"11"} export lupp=${lupp:-".true."} export corrlength=${corrlength:-1250} export lnsigcutoff=${lnsigcutoff:-2.75} export analpertwt=${analpertwt:-0.85} export readin_localization_enkf=".false." -export readin_localization_enkf=${readin_localization_enkf:-".true."} -export reducedgrid=${reducedgrid:-".true."} export letkf_flag=${letkf_flag:-".true."} export getkf=${getkf:-".true."} export denkf=${denkf:-".true."} diff --git a/regression/regression_namelists.sh b/regression/regression_namelists.sh index a4f283f92b..a6bd78994d 100755 --- a/regression/regression_namelists.sh +++ b/regression/regression_namelists.sh @@ -23,12 +23,12 @@ export gsi_namelist=" crtm_coeffs_path='./crtm_coeffs/', newpc4pred=.true.,adp_anglebc=.true.,angord=4,passive_bc=.true.,use_edges=.false., diag_precon=.true.,step_start=1.e-3,emiss_bc=.true.,thin4d=.true.,cwoption=3, - verbose=.false.,imp_physics=11,lupp=.true., + verbose=.false.,imp_physics=8,lupp=.true., binary_diag=.false.,netcdf_diag=.true., lobsdiag_forenkf=.false., nhr_anal=3,6,9,nhr_obsbin=1, l4densvar=.true.,ens_nstarthr=3,nhr_assimilation=6,lwrite4danl=.true., - optconv=0.05,cao_check=.true.,ta2tb=.false., + optconv=0.05,cao_check=.true.,ta2tb=.true., tzr_qc=1,sfcnst_comb=.true., write_fv3_incr=.true.,incvars_to_zero= 'liq_wat_inc','icmr_inc','rwmr_inc','snmr_inc','grle_inc', incvars_zero_strat='sphum_inc','liq_wat_inc','icmr_inc','rwmr_inc','snmr_inc','grle_inc',incvars_efold=5, use_gfs_ncio=.true., @@ -962,8 +962,8 @@ export gsi_namelist=" sprd_tol=1.e30,paoverpb_thresh=0.98, nlons=$LONA,nlats=$LATA,nlevs=$LEVS,nanals=$NMEM_ENKF, deterministic=.true.,sortinc=.true.,lupd_satbiasc=.false., - reducedgrid=.true.,readin_localization=.true., - use_gfs_nemsio=${use_gfs_nemsio},use_gfs_ncio=${use_gfs_ncio},imp_physics=$imp_physics,lupp=$lupp, + reducedgrid=.false.,readin_localization=.false., + use_gfs_nemsio=${use_gfs_nemsio},use_gfs_ncio=${use_gfs_ncio},imp_physics=8,lupp=$lupp, univaroz=.false.,adp_anglebc=.true.,angord=4,use_edges=.false.,emiss_bc=.true., letkf_flag=${letkf_flag},nobsl_max=${nobsl_max},denkf=${denkf},getkf=${getkf}., nhr_anal=${IAUFHRS_ENKF},nhr_state=${IAUFHRS_ENKF}, diff --git a/regression/regression_namelists_db.sh b/regression/regression_namelists_db.sh index b96c208070..34d80cba96 100755 --- a/regression/regression_namelists_db.sh +++ b/regression/regression_namelists_db.sh @@ -21,7 +21,7 @@ export gsi_namelist=" crtm_coeffs_path='./crtm_coeffs/', newpc4pred=.true.,adp_anglebc=.true.,angord=4,passive_bc=.true.,use_edges=.false., diag_precon=.true.,step_start=1.e-3,emiss_bc=.true.,thin4d=.true.,cwoption=3, - verbose=.false.,imp_physics=11,lupp=.true.,binary_diag=.false.,netcdf_diag=.true., + verbose=.false.,imp_physics=8,lupp=.true.,binary_diag=.false.,netcdf_diag=.true., lobsdiag_forenkf=.false., nhr_anal=3,6,9, l4densvar=.true.,ens_nstarthr=3,nhr_obsbin=1,nhr_assimilation=6,lwrite4danl=.true., tzr_qc=1,sfcnst_comb=.true., From de74bb2a5a14b61b6131a242c8bc600e1936aeb9 Mon Sep 17 00:00:00 2001 From: James Jung Date: Tue, 30 Jul 2024 20:16:10 -0400 Subject: [PATCH 04/11] Format changes for EUMETSAT metop-sg and CADS debug fix (#773) --- src/gsi/gsimod.F90 | 2 +- src/gsi/mrmsmod.f90 | 7 +++- src/gsi/obsmod.F90 | 3 +- src/gsi/qcmod.f90 | 2 +- src/gsi/radinfo.f90 | 4 +- src/gsi/read_cris.f90 | 89 +++++++++++++++++++++++------------------- src/gsi/read_diag.f90 | 8 ++-- src/gsi/read_iasi.f90 | 81 +++++++++++++++++++++----------------- src/gsi/read_obs.F90 | 59 ++++++++++++++-------------- src/gsi/read_viirs.f90 | 2 +- src/gsi/statsrad.f90 | 6 +-- 11 files changed, 144 insertions(+), 119 deletions(-) diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index d7f5667252..8c6977d78c 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -2320,7 +2320,7 @@ subroutine gsimain_initialize endif do i=1,ndat write(6,401)dfile(i),dtype(i),dplat(i),dsis(i),dval(i),dthin(i),dsfcalc(i),time_window(i) - 401 format(1x,a20,1x,a10,1x,a10,1x,a20,1x,f10.2,1x,I3,1x,I3,1x,f10.2) + 401 format(1x,a20,1x,a10,1x,a11,1x,a20,1x,f10.2,1x,I3,1x,I3,1x,f10.2) end do write(6,superob_radar) write(6,lag_data) diff --git a/src/gsi/mrmsmod.f90 b/src/gsi/mrmsmod.f90 index 713b974d7a..1cbad54904 100644 --- a/src/gsi/mrmsmod.f90 +++ b/src/gsi/mrmsmod.f90 @@ -27,7 +27,9 @@ subroutine load_mrms_data_info (mrms_listfile,nrows0,ntot_mrms,nrows_mrms,nrows, integer(i_kind),parameter::nobstype_mrms=2 ! first for vr, and the second for ref integer(i_kind),intent(in):: nrows0,ntot_mrms,nrows_mrms,nrows character(len=*),intent(in),optional ::mrms_listfile, rcname ! input filename - character(10),intent(inout),dimension(nrows):: dtype,ditype,dplat + character(len=*),intent(inout),dimension(nrows):: dtype + character(len=*),intent(inout),dimension(nrows):: ditype + character(len=*),intent(inout),dimension(nrows):: dplat character(20),intent(inout),dimension(nrows):: obsfile_all character(*),intent(inout),dimension(nrows):: dfile character(20),intent(inout),dimension(nrows):: dsis @@ -42,7 +44,8 @@ subroutine load_mrms_data_info (mrms_listfile,nrows0,ntot_mrms,nrows_mrms,nrows, real(r_kind),allocatable,dimension(:):: dmesh_mrms - character(10),allocatable,dimension(:):: dtype_mrms,ditype_mrms,dplat_mrms + character(10),allocatable,dimension(:):: dtype_mrms,ditype_mrms + character(11),allocatable,dimension(:):: dplat_mrms character(120),allocatable,dimension(:):: dfile_mrms character(20),allocatable,dimension(:):: dsis_mrms real(r_kind) ,allocatable,dimension(:):: dval_mrms diff --git a/src/gsi/obsmod.F90 b/src/gsi/obsmod.F90 index 19f6a072fb..eab9c828f7 100644 --- a/src/gsi/obsmod.F90 +++ b/src/gsi/obsmod.F90 @@ -626,7 +626,8 @@ module obsmod character(128) dirname character(128) obs_input_common character(20),allocatable,dimension(:):: obsfile_all - character(10),allocatable,dimension(:):: dtype,ditype,dplat + character(10),allocatable,dimension(:):: dtype,ditype + character(11),allocatable,dimension(:):: dplat character(120),allocatable,dimension(:):: dfile character(20),allocatable,dimension(:):: dsis real(r_kind) ,allocatable,dimension(:):: dval diff --git a/src/gsi/qcmod.f90 b/src/gsi/qcmod.f90 index 9804965573..2c4cbf6317 100644 --- a/src/gsi/qcmod.f90 +++ b/src/gsi/qcmod.f90 @@ -4240,7 +4240,7 @@ subroutine qc_goesimg(nchanl,is,ndat,nsig,ich,dplat,sea,land,ice,snow,luse, & real(r_kind),dimension(nsig,nchanl),intent(in ) :: temp,wmix real(r_kind),dimension(nchanl), intent(in ) :: tb_obs,tb_obs_sdv,tbc,tnoise,emissivity_k,ts real(r_kind),dimension(nchanl), intent(inout) :: errf,varinv - character(10), intent(in ) :: dplat + character(len=*), intent(in ) :: dplat ! Declare local parameters diff --git a/src/gsi/radinfo.f90 b/src/gsi/radinfo.f90 index 4ad17626e6..8bb496c645 100644 --- a/src/gsi/radinfo.f90 +++ b/src/gsi/radinfo.f90 @@ -805,10 +805,10 @@ subroutine radinfo_read end do close(lunin) 100 format(a1,a120) -110 format(i4,1x,a20,' chan= ',i4, & +110 format(i4,1x,a20,' chan= ',i5, & ' var= ',f7.3,' varch_cld=',f7.3,' use= ',i2,' ermax= ',F7.3, & ' b_rad= ',F7.2,' pg_rad=',F7.2,' icld_det=',I2,' icloud=',I2,' iaeros=',I2) -111 format(i4,1x,a20,' chan= ',i4, & +111 format(i4,1x,a20,' chan= ',i5, & ' var= ',f7.3,' varch_cld=',f7.3,' use= ',i2,' ermax= ',F7.3, & ' b_rad= ',F7.2,' pg_rad=',F7.2,' iextra_det=',I2, 'icloud=',I2,'iaeros=', I2) diff --git a/src/gsi/read_cris.f90 b/src/gsi/read_cris.f90 index d5668f8864..84288a7f04 100644 --- a/src/gsi/read_cris.f90 +++ b/src/gsi/read_cris.f90 @@ -112,9 +112,9 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& integer(i_kind) ,intent(in ) :: mype_sub integer(i_kind) ,intent(in ) :: npe_sub integer(i_kind) ,intent(in ) :: mpi_comm_sub - character(len=*), intent(in ) :: infile - character(len=10),intent(in ) :: jsatid - character(len=*), intent(in ) :: obstype + character(len=*) ,intent(in ) :: infile + character(len=*) ,intent(in ) :: jsatid + character(len=*) ,intent(in ) :: obstype character(len=20),intent(in ) :: sis real(r_kind) ,intent(in ) :: twind real(r_kind) ,intent(inout) :: val_cris @@ -200,7 +200,6 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& real(r_kind),dimension(83,7) :: imager_info real(r_kind),dimension(7) :: imager_cluster_size real(r_kind),dimension(2) :: imager_mean, imager_std_dev, imager_conversion - real(r_kind) :: imager_cluster_tot ! bufr error codes ! real(r_kind),dimension(7,3) :: error_codes @@ -848,10 +847,9 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& if ( cris_cads ) then call ufbseq(lnbufr,imager_info,83,7,iret,'CRISCS') if ( iret == 7 .and. imager_info(3,1) <= 100.0_r_kind .and. & - imager_info(3,1) >= zero .and. imager_coeff ) then ! if imager cluster info exists + sum(imager_info(3,:)) > zero .and. imager_coeff ) then ! if imager cluster info exists imager_mean = zero imager_std_dev = zero - imager_cluster_tot = zero imager_cluster_flag = .TRUE. imager_cluster_size = imager_info(3,1:7) imager_cluster_size(:) = imager_cluster_size(:) / sum(imager_cluster_size(:)) @@ -871,51 +869,62 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& imager_cluster_info: do j=1,7 i = imager_cluster_index(j) - data_all(maxinfo+j,itx) = imager_cluster_size(i) ! Imager cluster fraction - imager_cluster_tot = imager_cluster_tot + imager_info(3,i) +! If the cluster size, or radiance values of channel 4 or 5 are zero, do not compute statistics for that cluster - iexponent = -(nint(imager_info(75,i)) -11) ! channel 15 radiance for each cluster - imager_info(76,i) = imager_info(76,i) * imager_conversion(1) * (ten ** iexponent) + if ( imager_cluster_size(i) > zero .and. imager_info(76,i) > zero .and. imager_info(81,i) > zero ) then + data_all(maxinfo+j,itx) = imager_cluster_size(i) ! Imager cluster fraction - iexponent = -(nint(imager_info(77,i)) -11) ! channel 15 radiance std dev for each cluster. - imager_info(78,i) = imager_info(78,i) * imager_conversion(1) * (ten ** iexponent) + iexponent = -(nint(imager_info(75,i)) -11) ! channel 15 radiance for each cluster + imager_info(76,i) = imager_info(76,i) * imager_conversion(1) * (ten ** iexponent) - call crtm_planck_temperature(sensorindex_imager,4,imager_info(76,i),data_all(maxinfo+7+j,itx)) - data_all(maxinfo+7+j,itx) = max(data_all(maxinfo+7+j,itx),zero) + iexponent = -(nint(imager_info(77,i)) -11) ! channel 15 radiance std dev for each cluster. + imager_info(78,i) = imager_info(78,i) * imager_conversion(1) * (ten ** iexponent) - iexponent = -(nint(imager_info(80,i)) -11) ! channel 16 radiance for each cluster - imager_info(81,i) = imager_info(81,i) * imager_conversion(2) * (ten ** iexponent) + iexponent = -(nint(imager_info(80,i)) -11) ! channel 16 radiance for each cluster + imager_info(81,i) = imager_info(81,i) * imager_conversion(2) * (ten ** iexponent) - iexponent = -(nint(imager_info(82,i))-5 ) ! channel 16 radiance std dev for each cluster. - iexponent = -(nint(imager_info(82,i)) -11) ! channel 16 radiance std dev for each cluster. - imager_info(83,i) = imager_info(83,i) * imager_conversion(2) * (ten ** iexponent) + iexponent = -(nint(imager_info(82,i)) -11) ! channel 16 radiance std dev for each cluster. + imager_info(83,i) = imager_info(83,i) * imager_conversion(2) * (ten ** iexponent) - call crtm_planck_temperature(sensorindex_imager,5,imager_info(81,i),data_all(maxinfo+14+j,itx)) - data_all(maxinfo+14+j,itx) = max(data_all(maxinfo+14+j,itx),zero) - - - end do imager_cluster_info + call crtm_planck_temperature(sensorindex_imager,4,imager_info(76,i),data_all(maxinfo+7+j,itx)) + data_all(maxinfo+7+j,itx) = max(data_all(maxinfo+7+j,itx),zero) + call crtm_planck_temperature(sensorindex_imager,5,imager_info(81,i),data_all(maxinfo+14+j,itx)) + data_all(maxinfo+14+j,itx) = max(data_all(maxinfo+14+j,itx),zero) + else + data_all(maxinfo+j,itx) = zero ! something is wrong + data_all(maxinfo+7+j,itx) = zero ! set everything to zero + data_all(maxinfo+14+j,itx) = zero + endif + end do imager_cluster_info ! Compute cluster averages for each channel - imager_mean(1) = sum(imager_cluster_size(:) * imager_info(76,:)) ! Channel 15 radiance cluster average - imager_std_dev(1) = sum(imager_cluster_size(:) * (imager_info(76,:)**2 + imager_info(78,:)**2)) - imager_mean(1)**2 - imager_std_dev(1) = sqrt(max(imager_std_dev(1),zero)) ! Channel 15 radiance RMSE - call crtm_planck_temperature(sensorindex_imager,4,(imager_std_dev(1) + imager_mean(1)),imager_std_dev(1)) - call crtm_planck_temperature(sensorindex_imager,4,imager_mean(1),imager_mean(1)) ! Channel 15 average BT - imager_std_dev(1) = imager_std_dev(1) - imager_mean(1) ! Channel 15 BT std dev - data_all(maxinfo+22,itx) = imager_std_dev(1) - - imager_mean(2) = sum(imager_cluster_size(:) * imager_info(81,:)) ! Channel 16 radiance cluster average - imager_std_dev(2) = sum(imager_cluster_size(:) * (imager_info(81,:)**2 + imager_info(83,:)**2)) - imager_mean(1)**2 - imager_std_dev(2) = sqrt(max(imager_std_dev(1),zero)) ! Channel 16 radiance RMSE - call crtm_planck_temperature(sensorindex_imager,5,(imager_std_dev(2) + imager_mean(2)),imager_std_dev(2)) - call crtm_planck_temperature(sensorindex_imager,5,imager_mean(2),imager_mean(2)) ! Channel 16 average BT - imager_std_dev(2) = imager_std_dev(2) - imager_mean(2) ! Channel 16 BT std dev - data_all(maxinfo+23,itx) = imager_std_dev(2) + imager_mean(1) = sum(imager_cluster_size(:) * imager_info(76,:)) ! Channel 15 radiance cluster average + imager_std_dev(1) = sum(imager_cluster_size(:) * (imager_info(76,:)**2 + imager_info(78,:)**2)) - imager_mean(1)**2 + imager_std_dev(1) = sqrt(max(imager_std_dev(1),zero)) ! Channel 15 radiance RMSE + if ( imager_mean(1) > zero .and. imager_std_dev(1) > zero ) then + call crtm_planck_temperature(sensorindex_imager,4,(imager_std_dev(1) + imager_mean(1)),imager_std_dev(1)) + call crtm_planck_temperature(sensorindex_imager,4,imager_mean(1),imager_mean(1)) ! Channel 15 average BT + imager_std_dev(1) = imager_std_dev(1) - imager_mean(1) ! Channel 15 BT std dev + data_all(maxinfo+22,itx) = imager_std_dev(1) + else + data_all(maxinfo+22,itx) = zero + endif + + imager_mean(2) = sum(imager_cluster_size(:) * imager_info(81,:)) ! Channel 16 radiance cluster average + imager_std_dev(2) = sum(imager_cluster_size(:) * (imager_info(81,:)**2 + imager_info(83,:)**2)) - imager_mean(1)**2 + imager_std_dev(2) = sqrt(max(imager_std_dev(1),zero)) ! Channel 16 radiance RMSE + if ( imager_mean(2) > zero .and. imager_std_dev(2) > zero ) then + call crtm_planck_temperature(sensorindex_imager,5,(imager_std_dev(2) + imager_mean(2)),imager_std_dev(2)) + call crtm_planck_temperature(sensorindex_imager,5,imager_mean(2),imager_mean(2)) ! Channel 16 average BT + imager_std_dev(2) = imager_std_dev(2) - imager_mean(2) ! Channel 16 BT std dev + data_all(maxinfo+23,itx) = imager_std_dev(2) + else + data_all(maxinfo+23,itx) = zero + endif else ! Imager cluster information is missing. Set everything to zero - data_all(maxinfo+1 : maxinfo+25,itx) = zero + data_all(maxinfo+1 : maxinfo+cads_info,itx) = zero endif endif ! cris_cads diff --git a/src/gsi/read_diag.f90 b/src/gsi/read_diag.f90 index e389a708d5..f69f2f5133 100644 --- a/src/gsi/read_diag.f90 +++ b/src/gsi/read_diag.f90 @@ -88,7 +88,7 @@ module read_diag ! Declare structures for radiance diagnostic file information type diag_header_fix_list character(len=20) :: isis ! sat and sensor type - character(len=10) :: id ! sat type + character(len=11) :: id ! sat type character(len=10) :: obstype ! observation type integer(i_kind) :: jiter ! outer loop counter integer(i_kind) :: nchan ! number of channels in the sensor @@ -414,7 +414,8 @@ subroutine read_radiag_header_nc(ftin,header_fix,header_chan,iflag) real(r_kind),allocatable,dimension(:) :: r_var_stor integer(i_kind),allocatable,dimension(:) :: i_var_stor character(20) :: isis - character(10) :: id, obstype + character(11) :: id + character(10) :: obstype ! integer(i_kind),dimension(:),allocatable :: jiter, nchan_diag, npred, idate, & integer(i_kind) :: jiter, nchan_diag, npred, idate, & ireal, ipchan, iextra, jextra, & @@ -515,7 +516,8 @@ subroutine read_radiag_header_bin(ftin,npred_radiag,retrieval,header_fix,header_ ! Declare local variables character(len=2):: string - character(len=10):: satid,sentype + character(len=11):: satid + character(len=10):: sentype character(len=20):: sensat integer(i_kind) :: i,ich integer(i_kind):: jiter,nchanl,npred,ianldate,ireal,ipchan,iextra,jextra diff --git a/src/gsi/read_iasi.f90 b/src/gsi/read_iasi.f90 index edd7a9b50e..75dc50bb76 100644 --- a/src/gsi/read_iasi.f90 +++ b/src/gsi/read_iasi.f90 @@ -142,9 +142,9 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& integer(i_kind) ,intent(in ) :: mype_sub integer(i_kind) ,intent(in ) :: npe_sub integer(i_kind) ,intent(in ) :: mpi_comm_sub - character(len=*), intent(in ) :: infile - character(len=10),intent(in ) :: jsatid - character(len=*), intent(in ) :: obstype + character(len=*) ,intent(in ) :: infile + character(len=*) ,intent(in ) :: jsatid + character(len=*) ,intent(in ) :: obstype character(len=20),intent(in ) :: sis real(r_kind) ,intent(in ) :: twind real(r_kind) ,intent(inout) :: val_iasi @@ -233,7 +233,6 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& real(r_kind),dimension(33,7) :: imager_info real(r_kind),dimension(7) :: imager_cluster_size real(r_kind),dimension(2) :: imager_mean, imager_std_dev - real(r_kind) :: imager_cluster_tot ! Set standard parameters character(8),parameter:: fov_flag="crosstrk" @@ -813,10 +812,9 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& if ( iasi_cads ) then call ufbseq(lnbufr,imager_info,33,7,iret,'IASIL1CS') if (iret == 7 .and. imager_info(3,1) <= 100.0_r_kind .and. & - imager_info(3,1) >= zero .and. imager_coeff ) then ! if imager cluster info exists + sum(imager_info(3,:)) > zero .and. imager_coeff ) then ! if imager cluster info exists imager_mean = zero imager_std_dev = zero - imager_cluster_tot = zero imager_cluster_flag = .TRUE. imager_cluster_size = imager_info(3,1:7) imager_cluster_size(:) = imager_cluster_size(:) / sum(imager_cluster_size(:)) @@ -834,49 +832,62 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& imager_cluster_info: do j=1,7 i = imager_cluster_index(j) - data_all(maxinfo+j,itx) = imager_cluster_size(i) ! Imager cluster fraction - imager_cluster_tot = imager_cluster_tot + imager_info(3,i) +! If the cluster size, or radiance values of channel 4 or 5 are zero, do not compute statistics for that cluster + if ( imager_cluster_size(i) > zero .and. imager_info(26,i) > zero .and. imager_info(31,i) > zero ) then + data_all(maxinfo+j,itx) = imager_cluster_size(i) ! Imager cluster fraction - iexponent = -(nint(imager_info(25,i))-5 ) ! channel 4 radiance for each cluster. - imager_info(26,i) = imager_info(26,i) * (ten ** iexponent) + iexponent = -(nint(imager_info(25,i))-5 ) ! channel 4 radiance for each cluster. + imager_info(26,i) = imager_info(26,i) * (ten ** iexponent) - iexponent = -(nint(imager_info(27,i))-5 ) ! channel 4 radiance std dev for each cluster. - imager_info(28,i) = imager_info(28,i) * (ten ** iexponent) + iexponent = -(nint(imager_info(27,i))-5 ) ! channel 4 radiance std dev for each cluster. + imager_info(28,i) = imager_info(28,i) * (ten ** iexponent) - call crtm_planck_temperature(sensorindex_imager,2,imager_info(26,i),data_all(maxinfo+7+j,itx)) - data_all(maxinfo+7+j,itx) = max(data_all(maxinfo+7+j,itx),zero) + iexponent = -(nint(imager_info(30,i))-5 ) ! channel 5 radiance for each cluster + imager_info(31,i) = imager_info(31,i) * (ten ** iexponent) - iexponent = -(nint(imager_info(30,i))-5 ) ! channel 5 radiance for each cluster - imager_info(31,i) = imager_info(31,i) * (ten ** iexponent) + iexponent = -(nint(imager_info(32,i))-5 ) ! channel 5 radiance std dev for each cluser. + imager_info(33,i) = imager_info(33,i) * (ten ** iexponent) - iexponent = -(nint(imager_info(32,i))-5 ) ! channel 5 radiance std dev for each cluser. - imager_info(33,i) = imager_info(33,i) * (ten ** iexponent) - - call crtm_planck_temperature(sensorindex_imager,3,imager_info(31,i),data_all(maxinfo+14+j,itx)) - data_all(maxinfo+14+j,itx) = max(data_all(maxinfo+14+j,itx),zero) + call crtm_planck_temperature(sensorindex_imager,2,imager_info(26,i),data_all(maxinfo+7+j,itx)) + data_all(maxinfo+7+j,itx) = max(data_all(maxinfo+7+j,itx),zero) + call crtm_planck_temperature(sensorindex_imager,3,imager_info(31,i),data_all(maxinfo+14+j,itx)) + data_all(maxinfo+14+j,itx) = max(data_all(maxinfo+14+j,itx),zero) + else ! something is wrong + data_all(maxinfo+j,itx) = zero ! set everything to zero + data_all(maxinfo+7+j,itx) = zero + data_all(maxinfo+14+j,itx) = zero + endif end do imager_cluster_info ! Compute cluster averages for each channel - imager_mean(1) = sum(imager_cluster_size(:) * imager_info(26,:)) ! Channel 4 radiance cluster average + imager_mean(1) = sum(imager_cluster_size(:) * imager_info(26,:)) ! Channel 4 radiance cluster average imager_std_dev(1) = sum(imager_cluster_size(:) * (imager_info(26,:)**2 + imager_info(28,:)**2)) - imager_mean(1)**2 - imager_std_dev(1) = sqrt(max(imager_std_dev(1),zero)) ! Channel 4 radiance RMSE - call crtm_planck_temperature(sensorindex_imager,2,(imager_std_dev(1) + imager_mean(1)),imager_std_dev(1)) - call crtm_planck_temperature(sensorindex_imager,2,imager_mean(1),imager_mean(1)) ! Channel 4 average BT - imager_std_dev(1) = imager_std_dev(1) - imager_mean(1) ! Channel 4 BT std dev - data_all(maxinfo+22,itx) = imager_std_dev(1) - - imager_mean(2) = sum(imager_cluster_size(:) * imager_info(31,:)) ! Channel 5 radiance cluster average + imager_std_dev(1) = sqrt(max(imager_std_dev(1),zero)) ! Channel 4 radiance RMSE + if ( imager_mean(1) > zero .and. imager_std_dev(1) > zero ) then + call crtm_planck_temperature(sensorindex_imager,2,(imager_std_dev(1) + imager_mean(1)),imager_std_dev(1)) + call crtm_planck_temperature(sensorindex_imager,2,imager_mean(1),imager_mean(1)) ! Channel 4 average BT + imager_std_dev(1) = imager_std_dev(1) - imager_mean(1) ! Channel 4 BT std dev + data_all(maxinfo+22,itx) = imager_std_dev(1) + else + data_all(maxinfo+22,itx) = zero + endif + + imager_mean(2) = sum(imager_cluster_size(:) * imager_info(31,:)) ! Channel 5 radiance cluster average imager_std_dev(2) = sum(imager_cluster_size(:) * (imager_info(31,:)**2 + imager_info(33,:)**2)) - imager_mean(1)**2 - imager_std_dev(2) = sqrt(max(imager_std_dev(1),zero)) ! Channel 5 radiance RMSE - call crtm_planck_temperature(sensorindex_imager,3,(imager_std_dev(2) + imager_mean(2)),imager_std_dev(2)) - call crtm_planck_temperature(sensorindex_imager,3,imager_mean(2),imager_mean(2)) ! Channel 5 average BT - imager_std_dev(2) = imager_std_dev(2) - imager_mean(2) ! Channel 5 BT std dev - data_all(maxinfo+23,itx) = imager_std_dev(2) + imager_std_dev(2) = sqrt(max(imager_std_dev(1),zero)) ! Channel 5 radiance RMSE + if ( imager_mean(2) > zero .and. imager_std_dev(2) > zero ) then + call crtm_planck_temperature(sensorindex_imager,3,(imager_std_dev(2) + imager_mean(2)),imager_std_dev(2)) + call crtm_planck_temperature(sensorindex_imager,3,imager_mean(2),imager_mean(2)) ! Channel 5 average BT + imager_std_dev(2) = imager_std_dev(2) - imager_mean(2) ! Channel 5 BT std dev + data_all(maxinfo+23,itx) = imager_std_dev(2) + else + data_all(maxinfo+23,itx) = zero + endif else ! Imager cluster information is missing. Set everything to zero - data_all(maxinfo+1 : maxinfo+25,itx) = zero + data_all(maxinfo+1 : maxinfo+cads_info,itx) = zero endif endif ! iasi_cads = .true. ! diff --git a/src/gsi/read_obs.F90 b/src/gsi/read_obs.F90 index 6f0e32a4b6..b81bdb5752 100644 --- a/src/gsi/read_obs.F90 +++ b/src/gsi/read_obs.F90 @@ -791,7 +791,7 @@ subroutine read_obs(ndata,mype) logical,dimension(ndat):: belong,parallel_read,ears_possible,db_possible logical :: modis,use_sfc_any logical :: acft_profl_file - character(10):: obstype,platid + character(10):: obstype character(22):: string character(120):: infile character(20):: sis @@ -1415,7 +1415,6 @@ subroutine read_obs(ndata,mype) task_belongs: & if (i > 0 .and. belong(i)) then - platid=dplat(i) ! platid - satellites to read obstype=dtype(i) ! obstype - observation types to process infile=trim(dfile(i)) ! infile - units from which to read data sis=dsis(i) ! sensor/instrument/satellite indicator @@ -1553,12 +1552,12 @@ subroutine read_obs(ndata,mype) ! Process conventional SST (nsstbufr, at this moment) data elseif ( obstype == 'sst' ) then - string="--"//trim(ditype(i))//":sst:"//trim(platid) - if ( platid == 'nsst') then + string="--"//trim(ditype(i))//":sst:"//trim(dplat(i)) + if ( dplat(i) == 'nsst') then call read_nsstbufr(nread,npuse,nouse,gstime,infile,obstype, & lunout,twind,sis,nobs_sub1(1,i)) string='READ_NSSTBUFR' - elseif ( platid == 'mods') then + elseif ( dplat(i) == 'mods') then call read_modsbufr(nread,npuse,nouse,gstime,infile,obstype, & lunout,twind,sis,nobs_sub1(1,i)) string='READ_MODSBUFR' @@ -1716,12 +1715,12 @@ subroutine read_obs(ndata,mype) ! Process TOVS 1b data rad_obstype_select: & - if (platid /= 'aqua' .and. (obstype == 'amsua' .or. & + if (dplat(i) /= 'aqua' .and. (obstype == 'amsua' .or. & obstype == 'amsub' .or. obstype == 'msu' .or. & obstype == 'mhs' .or. obstype == 'hirs4' .or. & obstype == 'hirs3' .or. obstype == 'hirs2' .or. & obstype == 'ssu' )) then - call read_bufrtovs(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,& + call read_bufrtovs(mype,val_dat,ithin,isfcalc,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis, & mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), nobs_sub1(1,i), & read_rec(i),read_ears_rec(i),read_db_rec(i),dval_use,radmod) @@ -1729,7 +1728,7 @@ subroutine read_obs(ndata,mype) ! Process atms data else if (obstype == 'atms') then - call read_atms(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,& + call read_atms(mype,val_dat,ithin,isfcalc,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis, & mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i),nobs_sub1(1,i),& read_rec(i),read_ears_rec(i),read_db_rec(i),dval_use,radmod) @@ -1737,7 +1736,7 @@ subroutine read_obs(ndata,mype) ! Process saphir data else if (obstype == 'saphir') then - call read_saphir(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,& + call read_saphir(mype,val_dat,ithin,isfcalc,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis, & mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),dval_use) @@ -1745,9 +1744,9 @@ subroutine read_obs(ndata,mype) ! Process airs data - else if(platid == 'aqua' .and. (obstype == 'airs' .or. & + else if(dplat(i) == 'aqua' .and. (obstype == 'airs' .or. & obstype == 'amsua' .or. obstype == 'hsb' ))then - call read_airs(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,& + call read_airs(mype,val_dat,ithin,isfcalc,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis,& mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) @@ -1755,7 +1754,7 @@ subroutine read_obs(ndata,mype) ! Process iasi data else if(obstype == 'iasi')then - call read_iasi(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,& + call read_iasi(mype,val_dat,ithin,isfcalc,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis,& mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i),nobs_sub1(1,i), & read_rec(i),read_ears_rec(i),read_db_rec(i),dval_use) @@ -1763,7 +1762,7 @@ subroutine read_obs(ndata,mype) ! Process cris data else if(obstype == 'cris' .or. obstype =='cris-fsr' )then - call read_cris(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,& + call read_cris(mype,val_dat,ithin,isfcalc,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis,& mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i),nobs_sub1(1,i), & read_rec(i),read_ears_rec(i),read_db_rec(i),dval_use) @@ -1774,7 +1773,7 @@ subroutine read_obs(ndata,mype) else if (obstype == 'sndr' .or. & obstype == 'sndrd1' .or. obstype == 'sndrd2' .or. & obstype == 'sndrd3' .or. obstype == 'sndrd4') then - call read_goesndr(mype,val_dat,ithin,rmesh,platid,& + call read_goesndr(mype,val_dat,ithin,rmesh,dplat(i),& infile,lunout,obstype,nread,npuse,nouse,twind,gstime,sis,& mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) @@ -1782,7 +1781,7 @@ subroutine read_obs(ndata,mype) ! Process ssmi data else if (obstype == 'ssmi' ) then - call read_ssmi(mype,val_dat,ithin,rmesh,platid,gstime,& + call read_ssmi(mype,val_dat,ithin,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis,& mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) @@ -1791,7 +1790,7 @@ subroutine read_obs(ndata,mype) ! Process amsre data else if ( obstype == 'amsre_low' .or. obstype == 'amsre_mid' .or. & obstype == 'amsre_hig' ) then - call read_amsre(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,& + call read_amsre(mype,val_dat,ithin,isfcalc,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis,& mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) @@ -1801,7 +1800,7 @@ subroutine read_obs(ndata,mype) else if (obstype == 'ssmis' .or. & obstype == 'ssmis_las' .or. obstype == 'ssmis_uas' .or. & obstype == 'ssmis_img' .or. obstype == 'ssmis_env' ) then - call read_ssmis(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,& + call read_ssmis(mype,val_dat,ithin,isfcalc,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis,& mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) @@ -1809,7 +1808,7 @@ subroutine read_obs(ndata,mype) ! Process AMSR2 data else if(obstype == 'amsr2')then - call read_amsr2(mype,val_dat,ithin,rmesh,platid,gstime,& + call read_amsr2(mype,val_dat,ithin,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis,& mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i)) @@ -1817,7 +1816,7 @@ subroutine read_obs(ndata,mype) ! Process GOES IMAGER RADIANCE data else if(obstype == 'goes_img') then - call read_goesimg(mype,val_dat,ithin,rmesh,platid,gstime,& + call read_goesimg(mype,val_dat,ithin,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis, & mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) @@ -1825,7 +1824,7 @@ subroutine read_obs(ndata,mype) ! Process GMI data else if (obstype == 'gmi') then - call read_gmi(mype,val_dat,ithin,rmesh,platid,gstime,& + call read_gmi(mype,val_dat,ithin,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis,& mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),dval_use) @@ -1833,14 +1832,14 @@ subroutine read_obs(ndata,mype) ! Process Meteosat SEVIRI RADIANCE data else if(obstype == 'seviri') then - call read_seviri(mype,val_dat,ithin,rmesh,platid,gstime,& + call read_seviri(mype,val_dat,ithin,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis, & mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) string='READ_SEVIRI' ! Process GOES-R ABI RADIANCE data else if(obstype == 'abi') then - call read_abi(mype,val_dat,ithin,rmesh,platid,gstime,& + call read_abi(mype,val_dat,ithin,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis, & mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) @@ -1848,7 +1847,7 @@ subroutine read_obs(ndata,mype) ! Process Himawari-8 AHI RADIANCE data else if(obstype == 'ahi') then - call read_ahi(mype,val_dat,ithin,rmesh,platid,gstime,& + call read_ahi(mype,val_dat,ithin,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis, & mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) @@ -1857,7 +1856,7 @@ subroutine read_obs(ndata,mype) ! Process NAVY AVHRR RADIANCE data else if(obstype == 'avhrr_navy') then - call read_avhrr_navy(mype,val_dat,ithin,rmesh,platid,gstime,& + call read_avhrr_navy(mype,val_dat,ithin,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis, & mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) @@ -1865,14 +1864,14 @@ subroutine read_obs(ndata,mype) ! Process NESDIS AVHRR RADIANCE data else if(obstype == 'avhrr') then - call read_avhrr(mype,val_dat,ithin,rmesh,platid,gstime,& + call read_avhrr(mype,val_dat,ithin,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis, & mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) string='READ_AVHRR' ! Process SST VIIRS RADIANCE data else if(obstype == 'viirs-m') then - call read_sst_viirs(mype,val_dat,ithin,rmesh,platid,gstime,& + call read_sst_viirs(mype,val_dat,ithin,rmesh,dplat(i),gstime,& infile,lunout,obstype,nread,npuse,nouse,twind,sis, & mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i),read_rec(i),dval_use) @@ -1885,7 +1884,7 @@ subroutine read_obs(ndata,mype) if (is_extOzone(infile,obstype,dplat(i))) then call extOzone_read(infile,obstype,dplat(i),dsis(i), & - iread,ipuse,iouse, platid,gstime,lunout,twind,ithin,rmesh, & + iread,ipuse,iouse, dplat(i),gstime,lunout,twind,ithin,rmesh, & nobs_sub1(:,i)) string='extOzone_read' @@ -1895,7 +1894,7 @@ subroutine read_obs(ndata,mype) else call read_ozone(nread,npuse,nouse,& - platid,infile,gstime,lunout,obstype,twind,sis,ithin,rmesh, & + dplat(i),infile,gstime,lunout,obstype,twind,sis,ithin,rmesh, & nobs_sub1(1,i)) string='READ_OZONE' endif ozone_obstype_select @@ -1921,7 +1920,7 @@ subroutine read_obs(ndata,mype) ! Process aerosol data else if (ditype(i) == 'aero' )then call read_aerosol(nread,npuse,nouse,& - platid,infile,gstime,lunout,obstype,twind,sis,ithin,rmesh, & + dplat(i),infile,gstime,lunout,obstype,twind,sis,ithin,rmesh, & mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), & nobs_sub1(1,i)) string='READ_AEROSOL' @@ -1947,7 +1946,7 @@ subroutine read_obs(ndata,mype) call warn('read_obs',' infile =',trim(infile)) call warn('read_obs',' ditype(i) =',trim(ditype(i))) call warn('read_obs',' obstype =',trim(obstype)) - call warn('read_obs',' platid =',trim(platid)) + call warn('read_obs',' dplat =',trim(dplat(i))) call warn('read_obs',' string =',trim(string)) endif diff --git a/src/gsi/read_viirs.f90 b/src/gsi/read_viirs.f90 index ab80642f29..3ea50352ec 100644 --- a/src/gsi/read_viirs.f90 +++ b/src/gsi/read_viirs.f90 @@ -475,7 +475,7 @@ subroutine read_sst_viirs(mype,val_viirs,ithin,rmesh,jsatid,& nele,itxmax,nread,ndata_mesh,data_mesh,score_crit,nrec) if ( nread > 0 ) then - write(*,'(a,a10,I3,F6.1,3I10)') 'read_viirs,jsatid,imesh,amesh,itxmax,nread,ndata_mesh :',jsatid,imesh,amesh(imesh),itxmax,nread,ndata_mesh + write(*,'(a,a11,I3,F6.1,3I10)') 'read_viirs,jsatid,imesh,amesh,itxmax,nread,ndata_mesh :',jsatid,imesh,amesh(imesh),itxmax,nread,ndata_mesh endif ! ! get data_all by combining data from all thinning box sizes diff --git a/src/gsi/statsrad.f90 b/src/gsi/statsrad.f90 index d42a53f7d6..dfedcc92ab 100644 --- a/src/gsi/statsrad.f90 +++ b/src/gsi/statsrad.f90 @@ -100,7 +100,7 @@ subroutine statsrad(aivals,stats,ndata) write(iout_rad,4002) 'qcpenalty','qc1','qc2','qc3','qc4','qc5','qc6','qc7' write(iout_rad,4003) aivals(39,is),(nint(aivals(i,is)),i=8,14) 4000 format(2x,a3,7x,a4,14x,a7,1x,8(a7,1x)) -4001 format(1x,a10,1x,a10,f16.8,8(i7,1x)) +4001 format(1x,a11,1x,a10,f16.8,8(i7,1x)) 4002 format(28x,a9,1x,8(a7,1x)) 4003 format(22x,f16.8,8(i7,1x)) @@ -161,11 +161,11 @@ subroutine statsrad(aivals,stats,ndata) 2011 format(8x,f16.8,8(i7,1x)) 2012 format(12x,A7,5x,8(a7,1x)) 2999 format(' Illegal satellite type ') -1102 format(1x,i4,i5,1x,a16,2i7,1x,f10.3,1x,6(f11.7,1x)) +1102 format(1x,i4,i6,1x,a20,2i7,1x,f10.3,1x,6(f11.7,1x)) 1109 format(t5,'it',t13,'satellite',t23,'instrument',t40, & '# read',t53,'# keep',t65,'# assim',& t75,'penalty',t88,'qcpnlty',t104,'cpen',t115,'qccpen') -1115 format('o-g',1x,i2.2,1x,'rad',2x,2A10,2x,3(i11,2x),4(g12.5,1x)) +1115 format('o-g',1x,i2.2,1x,'rad',2x,2A12,2x,3(i11,2x),4(g12.5,1x)) ! Close output unit close(iout_rad) From 04896a7852fa7904851173231f797970d0de52d3 Mon Sep 17 00:00:00 2001 From: Azadeh Gholoubi <51101867+azadeh-gh@users.noreply.github.com> Date: Wed, 31 Jul 2024 05:47:29 -0400 Subject: [PATCH 05/11] Transition to Thompson Microphysics for Microwave All-sky Assimilation (#743) --- fix | 2 +- src/enkf/gridio_gfs.f90 | 6 - src/gsi/crtm_interface.f90 | 359 ++++++++++++++++++++++++++++++-- src/gsi/general_read_gfsatm.f90 | 120 +++++++++-- src/gsi/netcdfgfs_io.f90 | 28 ++- src/gsi/prt_guess.f90 | 64 +++++- src/gsi/set_crtm_cloudmod.f90 | 2 +- 7 files changed, 528 insertions(+), 53 deletions(-) diff --git a/fix b/fix index 15ffa60307..2a86d5b1a0 160000 --- a/fix +++ b/fix @@ -1 +1 @@ -Subproject commit 15ffa60307bbc19746d8caeb41782de0b7604be6 +Subproject commit 2a86d5b1a09a8908618d1c8a376c68e8d9b3d02c diff --git a/src/enkf/gridio_gfs.f90 b/src/enkf/gridio_gfs.f90 index 35a0c3fbe4..6e68f1a347 100644 --- a/src/enkf/gridio_gfs.f90 +++ b/src/enkf/gridio_gfs.f90 @@ -195,9 +195,6 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & sst_ind = getindex(vars2d, 'sst') use_full_hydro = ( ql_ind > 0 .and. qi_ind > 0 .and. & qr_ind > 0 .and. qs_ind > 0 .and. qg_ind > 0 ) - ! Currently, we do not let precipiation to affect the enkf analysis - ! The following line will be removed after testing - use_full_hydro = .false. if (.not. isinitialized) call init_spec_vars(nlons,nlats,ntrunc,4) @@ -459,9 +456,6 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & end if ! cloud derivatives - ! Currently, we do not let precipiation to affect the enkf analysis - ! The following line will be removed after testing - use_full_hydro = .true. if (.not. use_full_hydro .and. iope==0) then if (ql_ind > 0 .or. qi_ind > 0) then do k=1,nlevs diff --git a/src/gsi/crtm_interface.f90 b/src/gsi/crtm_interface.f90 index 4bb1191001..32a2cdb067 100644 --- a/src/gsi/crtm_interface.f90 +++ b/src/gsi/crtm_interface.f90 @@ -169,7 +169,8 @@ module crtm_interface real(r_kind) , save ,allocatable,dimension(:,:) :: cloudefr ! effective radius of cloud type in CRTM real(r_kind) , save ,allocatable,dimension(:,:) :: cloud_cont ! cloud content fed into CRTM real(r_kind) , save ,allocatable,dimension(:,:) :: cloud_efr ! effective radius of cloud type in CRTM - real(r_kind) , save ,allocatable,dimension(:) :: cf ! effective radius of cloud type in CRTM + real(r_kind) , save ,allocatable,dimension(:) :: cf ! cloud fraction + real(r_kind) , save ,allocatable,dimension(:) :: ni,nr ! ice and rain number concentration real(r_kind) , save ,allocatable,dimension(:) :: hwp_guess ! column total for each hydrometeor real(r_kind) , save ,allocatable,dimension(:) :: table,table2,tablew ! GFDL saturation water vapor pressure tables real(r_kind) , save ,allocatable,dimension(:) :: des2,desw ! GFDL saturation water vapor presure @@ -435,13 +436,17 @@ subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmo allocate(cloud(nsig,n_clouds_fwd)) allocate(cloudefr(nsig,n_clouds_fwd)) allocate(icloud(n_actual_clouds)) - allocate(cf(nsig)) + allocate(cf(nsig)) + allocate(ni(nsig)) + allocate(nr(nsig)) allocate(hwp_guess(n_clouds_fwd)) cloud_cont=zero cloud_efr =zero cloud =zero cloudefr =zero - cf =zero + cf =zero + ni =zero + nr =zero hwp_guess =zero call gsi_bundlegetpointer(gsi_metguess_bundle(1),cloud_names,icloud,ier) @@ -859,10 +864,9 @@ subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmo endif ! nvege_type endif ! regional or IGBP -! Initial GFDL saturation water vapor pressure tables +! Initial GFDL saturation water vapor pressure tables if (use_gfdl_qsat) then if (n_actual_aerosols_wk>0 .or. n_clouds_fwd_wk>0 .and. imp_physics==11) then - if (mype==0) write(6,*)myname_,':initial and load GFDL saturation water vapor pressure tables' allocate(table (length)) @@ -965,7 +969,9 @@ subroutine destroy_crtm if(allocated(jcloud)) deallocate(jcloud) if(allocated(cloud_cont)) deallocate(cloud_cont) if(allocated(cloud_efr)) deallocate(cloud_efr) - if(allocated(cf)) deallocate(cf) + if(allocated(cf)) deallocate(cf) + if(allocated(ni)) deallocate(ni) + if(allocated(nr)) deallocate(nr) if(allocated(hwp_guess)) deallocate(hwp_guess) if(allocated(icw)) deallocate(icw) if(allocated(lcloud4crtm_wk)) deallocate(lcloud4crtm_wk) @@ -1076,7 +1082,8 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & use crtm_module, only: limit_exp,o3_id,toa_pressure use obsmod, only: iadate use aeroinfo, only: nsigaerojac - use chemmod, only: lread_ext_aerosol !for separate aerosol input file + use chemmod, only: lread_ext_aerosol !for separate aerosol input file + use crtm_cloudcover_define, only: cloudcover_maximum_overlap, & cloudcover_random_overlap, & cloudcover_maxran_overlap, & @@ -1123,7 +1130,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & reshape((/0.0_r_kind, 1.0_r_kind, 1.0_r_kind, 2.0_r_kind, 1.0_r_kind, & -1.0_r_kind, 1.0_r_kind, -1.0_r_kind/), (/4, 2/)) real(r_kind),parameter:: jac_pert = 1.0_r_kind - + real(r_kind),parameter:: xrc3 = 200.0_r_kind ! Declare local variables integer(i_kind):: iquadrant integer(i_kind):: ier,ii,kk,kk2,i,itype,leap_day,day_of_year @@ -1132,12 +1139,12 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & integer(i_kind):: i_minus, i_plus, j_minus, j_plus integer(i_kind):: itsig,itsigp,itsfc,itsfcp,itaer,itaerp integer(i_kind):: istyp00,istyp01,istyp10,istyp11 - integer(i_kind):: iqs,iozs,icfs + integer(i_kind):: iqs,iozs,icfs + integer(i_kind):: inis,inrs integer(i_kind):: error_status_clr integer(i_kind):: idx700,dprs,dprs_min integer(i_kind),dimension(8)::obs_time,anal_time integer(i_kind),dimension(msig) :: klevel - ! ****************************** ! Constrained indexing for lai ! CRTM 2.1 implementation change @@ -1169,6 +1176,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & real(r_kind),dimension(nsig) :: rho_air ! density of air (kg/m3) real(r_kind),dimension(nsig) :: cf_calc ! GFDL cloud fraction calculation real(r_kind),dimension(nsig) :: qmix ! water vapor mixing ratio + real(r_kind),dimension(nsig) :: qcond real(r_kind),allocatable,dimension(:,:) :: tgas1d real(r_kind),pointer,dimension(:,: )::psges_itsig =>NULL() real(r_kind),pointer,dimension(:,: )::psges_itsigp=>NULL() @@ -1185,8 +1193,12 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & real(r_kind),pointer,dimension(:,:,:)::aeroges_itsig =>NULL() real(r_kind),pointer,dimension(:,:,:)::aeroges_itsigp=>NULL() real(r_kind),pointer,dimension(:,:,:)::cfges_itsig =>NULL() - real(r_kind),pointer,dimension(:,:,:)::cfges_itsigp=>NULL() - + real(r_kind),pointer,dimension(:,:,:)::cfges_itsigp=>NULL() + real(r_kind),pointer,dimension(:,:,:)::niges_itsig =>NULL() + real(r_kind),pointer,dimension(:,:,:)::niges_itsigp=>NULL() + real(r_kind),pointer,dimension(:,:,:)::nrges_itsig =>NULL() + real(r_kind),pointer,dimension(:,:,:)::nrges_itsigp=>NULL() + logical :: lmfdeep2 logical :: sea,icmask integer(i_kind),parameter,dimension(12):: mday=(/0,31,59,90,& @@ -1194,7 +1206,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & real(r_kind) :: lai m1=mype+1 - + if (mype==0) write(6,*) myname_, ' imp_physics = ', imp_physics if (n_clouds_fwd_wk>0) hwp_guess=zero hwp_total=zero theta_700=zero @@ -1210,6 +1222,8 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & cloud_efr = zero cloudefr = zero cf = zero + ni = zero + nr = zero endif dx = data_s(ilat) ! grid relative latitude @@ -1352,10 +1366,19 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & iqs=istatus call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'q',qges_itsigp,istatus) iqs=iqs+istatus - call gsi_bundlegetpointer(gsi_metguess_bundle(itsig ),'cf',cfges_itsig ,icfs) + call gsi_bundlegetpointer(gsi_metguess_bundle(itsig ),'cf',cfges_itsig ,istatus) icfs=istatus - call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'cf',cfges_itsigp,icfs) + call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'cf',cfges_itsigp,istatus) icfs=icfs+istatus + call gsi_bundlegetpointer(gsi_metguess_bundle(itsig ),'ni',niges_itsig,istatus) + inis=istatus + call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'ni',niges_itsigp,istatus) + inis=inis+istatus + call gsi_bundlegetpointer(gsi_metguess_bundle(itsig ),'nr',nrges_itsig,istatus) + inrs=istatus + call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'nr',nrges_itsigp,istatus) + inrs=inrs+istatus + ! Space-time interpolation of temperature (h) and q fields from sigma files !$omp parallel do schedule(dynamic,1) private(k,ii,iii) @@ -1730,8 +1753,10 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & ges_qsat (ixp,iy ,k, itsigp)*w10+ & ges_qsat (ix ,iyp,k, itsigp)*w01+ & ges_qsat (ixp,iyp,k, itsigp)*w11)*dtsigp + c3(k)=r1000/(one-q(k)) qmix(k)=q(k)*c3(k) !convert specific humidity to mixing ratio + rh(k) = q(k)/qs(k) !calculate relative humidity; added for cloud fraction computation ! Space-time interpolation of ozone(poz) if (iozs==0) then poz(k)=((ozges_itsig (ix ,iy ,k)*w00+ & @@ -1758,10 +1783,36 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & cfges_itsigp(ix ,iyp,k)*w01+ & cfges_itsigp(ixp,iyp,k)*w11)*dtsigp) -! Ensure ozone is greater than ozsmall +! Ensure cloud fraction is between 0 and 1 cf(k)=min(max(zero,cf(k)),one) endif ! cf + +! Space-time interpolation of ice number concentration (ni) + if (n_clouds_fwd_wk>0 .and. inis==0) then + ni(k)=((niges_itsig (ix ,iy ,k)*w00+ & + niges_itsig (ixp,iy ,k)*w10+ & + niges_itsig (ix ,iyp,k)*w01+ & + niges_itsig (ixp,iyp,k)*w11)*dtsig + & + (niges_itsigp(ix ,iy ,k)*w00+ & + niges_itsigp(ixp,iy ,k)*w10+ & + niges_itsigp(ix ,iyp,k)*w01+ & + niges_itsigp(ixp,iyp,k)*w11)*dtsigp) + ni(k)=max(zero,ni(k)) + endif ! ni + +! Space-time interpolation of rain number concentration (nr) + if (n_clouds_fwd_wk>0 .and. inrs==0) then + nr(k)=((nrges_itsig (ix ,iy ,k)*w00+ & + nrges_itsig (ixp,iy ,k)*w10+ & + nrges_itsig (ix ,iyp,k)*w01+ & + nrges_itsig (ixp,iyp,k)*w11)*dtsig + & + (nrges_itsigp(ix ,iy ,k)*w00+ & + nrges_itsigp(ixp,iy ,k)*w10+ & + nrges_itsigp(ix ,iyp,k)*w01+ & + nrges_itsigp(ixp,iyp,k)*w11)*dtsigp) + nr(k)=max(zero,nr(k)) + endif ! nr ! Quantities required for MW cloudy radiance calculations if (n_clouds_fwd_wk>0) then @@ -1822,8 +1873,17 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & end do endif + + ! Calculate Thompson effective radius for each hydrometeor + if ( icmask .and. n_clouds_fwd_wk > 0 .and. imp_physics==08 .and. lprecip_wk) then + do ii = 1, n_clouds_fwd_wk + iii=jcloud(ii) + call calc_thompson_reff(rho_air,h,cloud(:,ii),cloud_names(iii),cloudefr(:,ii)) + end do + endif + ! Calculate GFDL cloud fraction (if no cf in metguess table) based on PDF scheme -! if ( icmask .and. n_clouds_fwd_wk > 0 .and. imp_physics==11 .and. lcalc_gfdl_cfrac ) then + ! if ( icmask .and. n_clouds_fwd_wk > 0 .and. imp_physics==11 .and. lcalc_gfdl_cfrac ) then if ( icmask .and. n_clouds_fwd_wk > 0 .and. imp_physics==11 .and. lprecip_wk ) then cf_calc = zero call calc_gfdl_cloudfrac(rho_air,h,q,cloud,hs,garea,qs,cf_calc) @@ -1831,6 +1891,16 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & icfs = 0 ! load cloud fraction into CRTM endif + ! Calculate Thompson cloud fraction + if ( icmask .and. n_clouds_fwd_wk > 0 .and. imp_physics==08 .and. lprecip_wk ) then + lmfdeep2 = .true. + cf_calc = zero + ! sum of the cloud condensate amount for all 5 hydrometeos (ql + qi + qs + qg + qh + qr) + qcond(:) = cloud(:,1) + cloud(:,2) + cloud(:,3) + cloud(:,4) + cloud(:,5) + call calc_thompson_cloudfrac (nsig , lmfdeep2, xrc3, prsl, qcond(:), rh, qs, cf_calc) + cf = cf_calc + icfs = 0 ! load cloud fraction into CRTM + endif ! Interpolate level pressure to observation point for top interface prsi(nsig+1)=(ges_prsi(ix ,iy ,nsig+1,itsig )*w00+ & ges_prsi(ixp,iy ,nsig+1,itsig )*w10+ & @@ -2046,7 +2116,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & do ii=1,n_clouds_fwd_wk !cloud_cont(k,ii)=cloud(kk2,ii)*kgkg_kgm2 cloud_cont(k,ii)=cloud(kk2,ii)*c6(k) - if (imp_physics==11 .and. lprecip_wk .and. cloud_cont(k,ii) > 1.0e-6_r_kind) then + if (lprecip_wk .and. cloud_cont(k,ii) > 1.0e-6_r_kind) then cloud_efr(k,ii)=cloudefr(kk2,ii) else cloud_efr(k,ii)=zero @@ -2597,6 +2667,185 @@ subroutine calc_gfdl_reff(rho_air,tsen,qxmr,cloud_name,reff) end subroutine calc_gfdl_reff + subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) +!$$$ subprogram documentation block +! . . . . +! subprogram: calc_thompson_reff: calculate effective radius based on Thompson scheme +! prgmmr: Azadeh.Gholoubi +! +! program history log: +! +! 2023-01-24 +! +! language: f90 +! +!$$$ +!-------- + use constants, only: zero, pi + implicit none + +! Declare passed variables + character(10) ,intent(in ) :: cloud_name + real(r_kind), dimension(nsig) ,intent(in ) :: rho_air ! [ kg/m3 ] + real(r_kind), dimension(nsig) ,intent(in ) :: tsen ! [ K ] + real(r_kind), dimension(nsig) ,intent(in ) :: qxmr ! [ kg/kg ] + real(r_kind), dimension(nsig) ,intent(inout) :: reff ! [ micron ] + +! Declare local variables + character(len=*), parameter :: myname_ = 'calc_thompson_reff' + integer(i_kind) :: k + integer(i_kind) :: mu_w + real(r_kind) :: qx + real(r_kind) :: reff_min, reff_max + real(r_kind):: lam_i,lam_w,lam_r, lam_g,lam_exp, am_r,am_w,am_i,am_g + + ! Parameters + real(r_kind), parameter :: qmin = 1.0e-12_r_kind ! [kg/kg ] + !.. droplet number concentration. + real, parameter :: Nt_c = 100.E6_r_kind ![m-3] + + ! Parameters for water cloud + real(r_kind), parameter :: ccn = 1.0e8_r_kind + real(r_kind), parameter :: rho_w = 1000.0_r_kind ! [kg/m3 ] + real(r_kind), parameter :: reff_w_min = 2.0_r_kind ! previous value was 5.0_r_kind + real(r_kind), parameter :: reff_w_max = 25.0_r_kind ! previous value was 10.0_r_kind + + ! Parameters for ice cloud (Hemisfield and mcFarquhar 1996) + real(r_kind), parameter :: rho_i = 890.0_r_kind ! [kg/m3 ] + real(r_kind), parameter :: reff_i_min = 2.5_r_kind ! previous value was 10_r_kind + real(r_kind), parameter :: reff_i_max = 250.0_r_kind ! previous value was 150_r_kind + real(r_kind), parameter:: mu_i = 0.0_r_kind + + ! Parameters for rain (Lin 1983) + real(r_kind), parameter :: rho_r = 1000.0_r_kind ! [kg/m3 ] + real(r_kind), parameter :: reff_r_min = 50.0_r_kind ! [micron] ! previous value was 0.0_r_kind + real(r_kind), parameter :: reff_r_max = 1000.0_r_kind ! [micron] ! previous value was 10000.0_r_kind + real(r_kind), parameter:: mu_r = 0.0_r_kind + ! Parameters for snow + real(r_kind), parameter :: reff_s_min = 5.0_r_kind ! [micron] ! previous value was 0.0_r_kind + real(r_kind), parameter :: reff_s_max = 5000.0_r_kind ! [micron] ! previous value was 10000.0_r_kind + +!For snow moments conversions (from Field et al. 2005) + real(r_kind), dimension(10), parameter:: & + sa = (/ 5.065339_r_kind, -0.062659_r_kind, -3.032362_r_kind, 0.029469_r_kind, -0.000285_r_kind, & + & 0.31255_r_kind, 0.000204_r_kind, 0.003199_r_kind, 0.0_r_kind, -0.015952_r_kind/) + real(r_kind), dimension(10), parameter:: & + sb = (/ 0.476221_r_kind, -0.015896_r_kind, 0.165977_r_kind, 0.007468_r_kind, -0.000141_r_kind, & + & 0.060366_r_kind, 0.000079_r_kind, 0.000594_r_kind, 0.0_r_kind, -0.003577_r_kind/) + real(r_kind), parameter :: am_s = 0.069_r_kind + real(r_kind), parameter :: bm_s = 2.0_r_kind + real(r_kind), dimension(1), parameter :: cse = (/ bm_s + 1.0_r_kind /) + real(r_kind) :: tc0, smob, smoc, a_, b_, loga_ + + ! Parameters for graupel (Lin 1983) + real(r_kind), parameter :: rho_g = 500.0_r_kind ! [kg/m3 ] + real(r_kind), parameter :: reff_g_min = 150.0_r_kind ! [micron] ! previous value was 0.0_r_kind + real(r_kind), parameter :: reff_g_max = 5000.0_r_kind ! [micron] ! previous value was 10000.0_r_kind + real(r_kind), parameter:: mu_g = 0.0_r_kind + real(r_kind), parameter :: no_exp = 10.0_r_kind + +!cloud water + if (trim(cloud_name)=='ql') then + am_w = rho_w*pi/6.0_r_kind + reff_min = reff_w_min + reff_max = reff_w_max + do k = 1, nsig + qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) + if (qx > qmin) then + mu_w = MAX(2, MIN((NINT(1000.E6_r_kind/Nt_c) + 2), 15)) + lam_w=exp(1.0_r_kind / 3.0_r_kind * log ((am_w*Nt_c *gamma(mu_w + 3.0_r_kind + 1.0_r_kind))/(qx*gamma(mu_w+1.0_r_kind)))) + + reff(k) = 0.5_r_kind * ((3.0_r_kind+mu_w)/lam_w)*1.0e6_r_kind + reff(k) = max(reff_min, min(reff_max, reff(k))) + else + reff(k) = zero + endif + enddo + + ! Cloud Ice + else if (trim(cloud_name)=='qi') then + am_i = rho_i*pi/6.0_r_kind + reff_min = reff_i_min + reff_max = reff_i_max + do k = 1, nsig + qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) + if (qx > qmin) then + lam_i=exp(1.0_r_kind / 3.0_r_kind * log((am_i*ni(k) *gamma(mu_i + 3.0_r_kind + 1.0_r_kind))/(qx*gamma(mu_i+1.0_r_kind)))) + reff(k) = 0.5_r_kind * (3.0_r_kind /lam_i)*1.0e6_r_kind + reff(k) = max(reff_min, min(reff_max, reff(k))) + else + reff(k) = zero + endif + enddo + !Rain + else if (trim(cloud_name)=='qr') then + am_r = rho_r*pi/6.0_r_kind + reff_min = reff_r_min + reff_max = reff_r_max + do k = 1, nsig + qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) + if (qx > qmin) then + lam_r=exp(1.0_r_kind / 3.0_r_kind * log ((am_r*nr(k) *gamma(mu_r + 3.0_r_kind + 1.0_r_kind))/(qx*gamma(mu_r + 1.0_r_kind)))) + reff(k) = 0.5_r_kind *(3.0_r_kind/lam_r)*1.0e6_r_kind + reff(k) = max(reff_min, min(reff_max, reff(k))) + else + reff(k) = zero + endif + enddo + +! Snow (Field et al. 2005) + + else if (trim(cloud_name)=='qs') then + reff_min = reff_s_min + reff_max = reff_s_max + do k = 1, nsig +!..Calculate bm_s+1 (th) moment, smoc. Useful for diameter calcs. + tc0 = MIN(-0.1_r_kind, tsen(k)-273.15_r_kind) + qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) + if (qx > qmin) then + smob =qx/am_s + loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) & + + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 & + + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) & + + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 & + + sa(10)*cse(1)*cse(1)*cse(1) + a_ = 10.0_r_kind**loga_ + b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) & + + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) & + + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) & + + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1) + smoc = a_ * smob**b_ + reff(k) = MAX(2.51E-6_r_kind, MIN(0.5_r_kind*(smoc/smob), 1999.E-6_r_kind)) + reff(k) = max(reff_min, min(reff_max, reff(k)*1.0e6_r_kind)) + else + reff(k) = zero + endif + enddo + + ! Graupel + else if (trim(cloud_name)=='qg') then + reff_min = reff_g_min + reff_max = reff_g_max + am_g = rho_g*pi/6.0_r_kind + do k = 1, nsig + qx = qxmr(k)*rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) + if (qx > qmin) then + lam_exp = no_exp* exp(1.0_r_kind/(3.0_r_kind+1.0_r_kind) * log ((am_g*gamma(3.0_r_kind +1.0_r_kind))/qx)) + lam_g=lam_exp*exp(1/3.0_r_kind*log(gamma(3.0_r_kind+mu_g+1.0_r_kind)/((3.0_r_kind+mu_g+1.0_r_kind)*(mu_g+1.0_r_kind)))) + reff(k) = 0.5_r_kind *(3.0_r_kind/ lam_g)*1.0e6_r_kind + reff(k) = max(reff_min, min(reff_max, reff(k))) + else + reff(k) = zero + endif + enddo + + ! Mysterious + else + call die(myname_,"cannot recognize cloud name <"//trim(myname_)//">") + endif + + end subroutine calc_thompson_reff + subroutine calc_gfdl_cloudfrac(den,pt1,qv,cloud,hs,area,qsat,cfrac) !$$$ subprogram documentation block ! . . . . @@ -2810,6 +3059,80 @@ subroutine calc_gfdl_cloudfrac(den,pt1,qv,cloud,hs,area,qsat,cfrac) enddo !k-loop end subroutine calc_gfdl_cloudfrac + subroutine calc_thompson_cloudfrac(nsig, lmfdeep2, xrc3, prsl, clwf, rhly, qstl,cldtot) +!$$$ subprogram documentation block +! . . . . +! subprogram: calc_thompson_cloudfrac calculate cloud fractions from cloud using Thompson scheme +! +! prgmmr: Azadeh.Gholoubi +! +! abstract: +! +! program history log: +! +! 2023-01-24 +! +! language: f90 +! +! ==================== definition of variables ==================== ! +! input variables: ! +! prsl (msig) : model layer mean pressure in cb ! +! qstl (nsig) : layer saturate humidity in gm/gm ! +! rhly (nsig) : layer relative humidity (=qlyr/qstl) ! +! nsig : vertical layer/level dimensions ! +! lmfdeep2 : logical - true for mass flux deep convection ! +! clwf : sum of the all 5 Hydrometeors ! +! lmfdeep2 : scale-aware mass-flux deep conv scheme flag ! ! +! output variables: ! +! cldtot (nsig) - layer total cloud fraction ! +! ==================== end of description ===================== ! +!$$$ +!-------- + use constants + implicit none + +! --- inputs: + integer(i_kind), intent(in) :: nsig + real(r_kind) , dimension(nsig), intent(in) :: prsl, clwf, rhly, qstl + logical, intent(in) :: lmfdeep2 + real(r_kind) , intent(in):: xrc3 + +! --- outputs + real(r_kind) , dimension(nsig), intent(inout) :: cldtot + +! --- local variables: + real(r_kind) :: clwmin, clwm, clwt, onemrh, value, tem1, tem2 + integer(i_kind) :: k + +!> - Compute layer cloud fraction. + clwmin = 0.0_r_kind + do k = 1, nsig-1 +! converting prsl from cb to bar + clwt = 1.0e-10_r_kind * (prsl(k)*0.01_r_kind) + if (clwf(k) > clwt) then + if(rhly(k) > 0.99_r_kind) then + cldtot(k) = 1.0_r_kind + else + onemrh= max( 1.e-10_r_kind, 1.0_r_kind-rhly(k) ) + clwm = clwmin / max( 0.01_r_kind, prsl(k)*0.01_r_kind ) + tem1 = min(max((onemrh*qstl(k))**0.49_r_kind,0.0001_r_kind),1.0_r_kind) + if (lmfdeep2) then + tem1 = xrc3 / tem1 + else + tem1 = 100.0_r_kind / tem1 + endif + value = max( min( tem1*(clwf(k)-clwm), 50.0_r_kind ), 0.0_r_kind ) + tem2 = sqrt( sqrt(rhly(k)) ) + cldtot(k) = max( tem2*(1.0_r_kind-exp(-value)), 0.0_r_kind ) + endif + else + cldtot(k) = 0.0_r_kind + endif + enddo + + end subroutine calc_thompson_cloudfrac +!............................................ + subroutine qs_table(n) !$$$ subprogram documentation block ! . . . . diff --git a/src/gsi/general_read_gfsatm.f90 b/src/gsi/general_read_gfsatm.f90 index 2971f27148..c34b269950 100755 --- a/src/gsi/general_read_gfsatm.f90 +++ b/src/gsi/general_read_gfsatm.f90 @@ -189,13 +189,14 @@ subroutine general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, end subroutine general_reload subroutine general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vdflag,g_cf) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vdflag,g_ni,g_nr,g_cf) ! !USES: use kinds, only: r_kind,i_kind use mpimod, only: npe,mpi_comm_world,ierror,mpi_rtype use general_sub2grid_mod, only: sub2grid_info + use ncepnems_io, only: imp_physics implicit none ! !INPUT PARAMETERS: @@ -211,7 +212,8 @@ subroutine general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & real(r_kind),dimension(grd%lat2,grd%lon2), intent( out) :: g_ps real(r_kind),dimension(grd%lat2,grd%lon2), intent(inout) :: g_z real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig),intent( out) :: g_u,g_v,& - g_vor,g_div,g_q,g_oz,g_tv,g_ql,g_qi,g_qr,g_qs,g_qg + g_vor,g_div,g_q,g_oz,g_tv,g_ql,g_qi,g_qr,g_qs,g_qg,g_ni,g_nr + real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig),intent( out),optional :: g_cf @@ -392,7 +394,25 @@ subroutine general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & g_qg(i,j,klev)=sub(ij,k) enddo enddo - elseif ( iflag(k) == 15 .and. present(g_cf) ) then + elseif ( iflag(k) == 15 .and. imp_physics == 8) then + klev=ilev(k) + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_ni(i,j,klev)=sub(ij,k) + enddo + enddo + elseif ( iflag(k) == 16 .and. imp_physics == 8) then + klev=ilev(k) + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_nr(i,j,klev)=sub(ij,k) + enddo + enddo + elseif ( iflag(k) == 17 .and. present(g_cf) ) then klev=ilev(k) ij=0 do j=1,grd%lon2 @@ -2828,7 +2848,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z real(r_kind),pointer,dimension(:,:) :: g_t2m, g_q2m real(r_kind),pointer,dimension(:,:,:) :: g_vor,g_div,& g_q,g_oz,g_tv - real(r_kind),pointer,dimension(:,:,:) :: g_ql,g_qi,g_qr,g_qs,g_qg + real(r_kind),pointer,dimension(:,:,:) :: g_ql,g_qi,g_qr,g_qs,g_qg,g_ni,g_nr real(r_kind),allocatable,dimension(:,:) :: g_z real(r_kind),allocatable,dimension(:,:,:) :: g_u,g_v @@ -3038,6 +3058,8 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z call gsi_bundlegetpointer(gfs_bundle,'qs',g_qs ,ier);istatus1=istatus1+ier call gsi_bundlegetpointer(gfs_bundle,'qg',g_qg ,ier);istatus1=istatus1+ier ! call gsi_bundlegetpointer(gfs_bundle,'cf',g_cf ,ier);istatus1=istatus1+ier + call gsi_bundlegetpointer(gfs_bundle,'ni',g_ni ,ier);istatus1=istatus1+ier + call gsi_bundlegetpointer(gfs_bundle,'nr',g_nr ,ier);istatus1=istatus1+ier if ( istatus1 /= 0 ) then if ( mype == 0 ) then write(6,*) 'general_read_gfsatm_allhydro_nc: ERROR' @@ -3100,7 +3122,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif endif @@ -3135,7 +3157,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif ! Thermodynamic variable: s-->g transform, communicate to all tasks @@ -3172,7 +3194,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif end do @@ -3229,7 +3251,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif end do do k=1,nlevs @@ -3285,7 +3307,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif end do @@ -3321,7 +3343,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif icount=icount+1 @@ -3352,7 +3374,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif end do endif ! if ( uvflag ) @@ -3383,7 +3405,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif end do @@ -3414,7 +3436,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif end do @@ -3445,7 +3467,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif enddo ! do k=1,nlevs @@ -3475,7 +3497,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif enddo ! do k=1,nlevs @@ -3505,7 +3527,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif enddo ! do k=1,nlevs @@ -3535,7 +3557,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm ) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif enddo ! do k=1,nlevs @@ -3565,13 +3587,73 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z endif if ( icount == icm .or. k==nlevs) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) + endif + enddo ! do k=1,nlevs + + do k=1,nlevs + icount=icount+1 + iflag(icount)=15 + ilev(icount)=k + kr = levs+1-k ! netcdf is top to bottom, need to flip + + if (mype==mype_use(icount)) then + call read_vardata(filges, 'nccice', rwork3d0, nslice=kr, slicedim=3) + ! cloud ice water number concentration. + if ( diff_res ) then + grid_b=rwork3d0(:,:,1) + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork3d0(:,:,1) + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm .or. k==nlevs ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) + endif + enddo ! do k=1,nlevs + + do k=1,nlevs + icount=icount+1 + iflag(icount)=16 + ilev(icount)=k + kr = levs+1-k ! netcdf is top to bottom, need to flip + + if (mype==mype_use(icount)) then + call read_vardata(filges, 'nconrd', rwork3d0, nslice=kr, slicedim=3) + ! rain number concentration. + if ( diff_res ) then + grid_b=rwork3d0(:,:,1) + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork3d0(:,:,1) + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm .or. k==nlevs) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) endif enddo ! do k=1,nlevs ! do k=1,nlevs ! icount=icount+1 -! iflag(icount)=15 +! iflag(icount)=17 ! ilev(icount)=k ! kr = levs+1-k ! netcdf is top to bottom, need to flip ! diff --git a/src/gsi/netcdfgfs_io.f90 b/src/gsi/netcdfgfs_io.f90 index 8aa6c4cc6a..7db02636cd 100644 --- a/src/gsi/netcdfgfs_io.f90 +++ b/src/gsi/netcdfgfs_io.f90 @@ -159,17 +159,17 @@ subroutine read_ real(r_kind),pointer,dimension(:,:,:):: ges_qs_it => NULL() real(r_kind),pointer,dimension(:,:,:):: ges_qg_it => NULL() real(r_kind),pointer,dimension(:,:,:):: ges_cf_it => NULL() - + real(r_kind),pointer,dimension(:,:,:):: ges_ni_it => NULL() + real(r_kind),pointer,dimension(:,:,:):: ges_nr_it => NULL() type(sub2grid_info) :: grd_t logical regional logical:: l_cld_derived,zflag,inithead - type(gsi_bundle) :: atm_bundle type(gsi_grid) :: atm_grid integer(i_kind),parameter :: n2d=2 ! integer(i_kind),parameter :: n3d=8 integer(i_kind),parameter :: n2d_2m=4 - integer(i_kind),parameter :: n3d=14 + integer(i_kind),parameter :: n3d=16 character(len=4), parameter :: vars2d(n2d) = (/ 'z ', 'ps ' /) character(len=4), parameter :: vars2d_with2m(n2d_2m) = (/ 'z ', 'ps ','t2m ','q2m ' /) ! character(len=4), parameter :: vars3d(n3d) = (/ 'u ', 'v ', & @@ -182,13 +182,16 @@ subroutine read_ 'cw ', 'oz ', & 'ql ', 'qi ', & 'qr ', 'qs ', & - 'qg ', 'cf ' /) + 'qg ', 'ni ', & + 'nr ', 'cf ' /) + real(r_kind),pointer,dimension(:,:):: ptr2d =>NULL() real(r_kind),pointer,dimension(:,:,:):: ptr3d =>NULL() regional=.false. inner_vars=1 - num_fields=min(14*grd_a%nsig+2,npe) + + num_fields=min(n3d*grd_a%nsig+2,npe) ! Create temporary communication information fore read routines call general_sub2grid_create_info(grd_t,inner_vars,grd_a%nlat,grd_a%nlon, & grd_a%nsig,num_fields,regional) @@ -200,6 +203,7 @@ subroutine read_ else call gsi_bundlecreate(atm_bundle,atm_grid,'aux-atm-read',istatus,names2d=vars2d,names3d=vars3d) endif + if(istatus/=0) then write(6,*) myname_,': trouble creating atm_bundle' call stop2(999) @@ -248,6 +252,8 @@ subroutine read_ if (associated(ges_qr_it)) ges_qr_it(i,j,k) = max(qcmin,ges_qr_it(i,j,k)) if (associated(ges_qs_it)) ges_qs_it(i,j,k) = max(qcmin,ges_qs_it(i,j,k)) if (associated(ges_qg_it)) ges_qg_it(i,j,k) = max(qcmin,ges_qg_it(i,j,k)) + if (associated(ges_ni_it)) ges_ni_it(i,j,k) = max(qcmin,ges_ni_it(i,j,k)) + if (associated(ges_nr_it)) ges_nr_it(i,j,k) = max(qcmin,ges_nr_it(i,j,k)) if (associated(ges_cf_it)) ges_cf_it(i,j,k) = min(max(zero,ges_cf_it(i,j,k)),one) enddo enddo @@ -256,6 +262,8 @@ subroutine read_ l_cld_derived = associated(ges_cwmr_it).and.& associated(ges_q_it) .and.& associated(ges_ql_it) .and.& + associated(ges_ni_it) .and.& + associated(ges_nr_it) .and.& associated(ges_qi_it) .and.& associated(ges_tv_it) ! call set_cloud_lower_bound(ges_cwmr_it) @@ -343,6 +351,16 @@ subroutine set_guess_ call gsi_bundlegetpointer (gsi_metguess_bundle(it),'ql',ges_ql_it,istatus) if(istatus==0) ges_ql_it = ptr3d endif + call gsi_bundlegetpointer (atm_bundle,'ni',ptr3d,istatus) + if (istatus==0) then + call gsi_bundlegetpointer (gsi_metguess_bundle(it),'ni',ges_ni_it,istatus) + if(istatus==0) ges_ni_it = ptr3d + endif + call gsi_bundlegetpointer (atm_bundle,'nr',ptr3d,istatus) + if (istatus==0) then + call gsi_bundlegetpointer (gsi_metguess_bundle(it),'nr',ges_nr_it,istatus) + if(istatus==0) ges_nr_it = ptr3d + endif call gsi_bundlegetpointer (atm_bundle,'qi',ptr3d,istatus) if (istatus==0) then call gsi_bundlegetpointer (gsi_metguess_bundle(it),'qi',ges_qi_it,istatus) diff --git a/src/gsi/prt_guess.f90 b/src/gsi/prt_guess.f90 index 97c1e0a7c1..dbce7e0f79 100644 --- a/src/gsi/prt_guess.f90 +++ b/src/gsi/prt_guess.f90 @@ -287,7 +287,8 @@ subroutine prt_guess2(sgrep) integer(i_kind), parameter :: nvars3=5 integer(i_kind) :: nvars,nvars2,nvarsc,nc integer(i_kind) ii,istatus,ier,ivar - integer(i_kind) iql,iqi,iqr,iqs,iqg,icf + integer(i_kind) iql,iqi,iqr,iqs,iqg,icf + integer(i_kind) ini,inr integer(i_kind) ntsig integer(i_kind) ntsfc integer(i_kind) n_actual_clouds @@ -307,6 +308,8 @@ subroutine prt_guess2(sgrep) real(r_kind),pointer,dimension(:,:,:)::ges_qr_it => NULL() real(r_kind),pointer,dimension(:,:,:)::ges_qs_it => NULL() real(r_kind),pointer,dimension(:,:,:)::ges_qg_it => NULL() + real(r_kind),pointer,dimension(:,:,:)::ges_ni_it => NULL() + real(r_kind),pointer,dimension(:,:,:)::ges_nr_it => NULL() real(r_kind),pointer,dimension(:,:,:)::ges_cf_it => NULL() ! character(len=4) :: cvar(nvars+3) character(len=4),allocatable,dimension(:) :: cvar @@ -338,6 +341,7 @@ subroutine prt_guess2(sgrep) ! get pointer to cloud water condensate ier=0;nvarsc=0 iql=0;iqi=0;iqr=0;iqs=0;iqg=0 + ini=0;inr=0 call gsi_metguess_get('clouds::3d',n_actual_clouds,ier) if (mype==0) write(6,*)'prt_guess2: n_actual_clouds = ', n_actual_clouds if (n_actual_clouds>0) then @@ -345,36 +349,52 @@ subroutine prt_guess2(sgrep) if (ivar > 0) then call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'ql',ges_ql_it,istatus) ier=ier+istatus + if (ier==0) nvarsc=nvarsc+1 endif call gsi_metguess_get ( 'var::qi', ivar, ier ); iqi=ivar if (ivar > 0) then call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'qi',ges_qi_it,istatus) ier=ier+istatus + if (ier==0) nvarsc=nvarsc+1 endif call gsi_metguess_get ( 'var::qr', ivar, ier ); iqr=ivar if (ivar > 0) then call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'qr',ges_qr_it,istatus) ier=ier+istatus + if (ier==0) nvarsc=nvarsc+1 endif call gsi_metguess_get ( 'var::qs', ivar, ier ); iqs=ivar if (ivar > 0) then call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'qs',ges_qs_it,istatus) ier=ier+istatus + if (ier==0) nvarsc=nvarsc+1 endif call gsi_metguess_get ( 'var::qg', ivar, ier ); iqg=ivar if (ivar > 0) then call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'qg',ges_qg_it,istatus) ier=ier+istatus + if (ier==0) nvarsc=nvarsc+1 + endif + call gsi_metguess_get ( 'var::ni', ivar, ier ); ini=ivar + if (ivar > 0) then + call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'ni',ges_ni_it,istatus) + ier=ier+istatus + if (ier==0) nvarsc=nvarsc+1 + endif + call gsi_metguess_get ( 'var::nr', ivar, ier ); inr=ivar + if (ivar > 0) then + call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'nr',ges_nr_it,istatus) + ier=ier+istatus + if (ier==0) nvarsc=nvarsc+1 endif end if - nvarsc=n_actual_clouds call gsi_metguess_get ( 'var::cf', ivar, ier ); icf=ivar if (ivar > 0) then call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'cf',ges_cf_it,istatus) ier=ier+istatus nvarsc=nvarsc+1 endif -! if (ier/=0) return ! this is a fundamental routine, when some not found just ! emily +! if (ier/=0) return ! this is a fundamental routine, when some not found just ! nvars = nvars1+nvarsc+nvars3 nvars2 = nvars1+nvarsc @@ -415,6 +435,14 @@ subroutine prt_guess2(sgrep) nc=nc+1 cvar(nvars1+nc)='QG ' endif + if(ini>0) then + nc=nc+1 + cvar(nvars1+nc)='NI ' + endif + if(inr>0) then + nc=nc+1 + cvar(nvars1+nc)='NR ' + endif if(icf>0) then nc=nc+1 cvar(nvars1+nc)='CF ' @@ -455,6 +483,17 @@ subroutine prt_guess2(sgrep) nc=nc+1 zloc(nvars1+nc) = sum (ges_qg_it (2:lat1+1,2:lon1+1,1:nsig)) endif + + if(ini>0) then + nc=nc+1 + zloc(nvars1+nc) = sum (ges_ni_it (2:lat1+1,2:lon1+1,1:nsig)) + endif + + if(inr>0) then + nc=nc+1 + zloc(nvars1+nc) = sum (ges_nr_it (2:lat1+1,2:lon1+1,1:nsig)) + endif + if(icf>0) then nc=nc+1 zloc(nvars1+nc) = sum (ges_cf_it (2:lat1+1,2:lon1+1,1:nsig)) @@ -491,6 +530,16 @@ subroutine prt_guess2(sgrep) nc=nc+1 zloc(nvars+nvars1+nc) = minval(ges_qg_it (2:lat1+1,2:lon1+1,1:nsig)) endif + + if(ini>0) then + nc=nc+1 + zloc(nvars+nvars1+nc) = minval(ges_ni_it (2:lat1+1,2:lon1+1,1:nsig)) + endif + if(inr>0) then + nc=nc+1 + zloc(nvars+nvars1+nc) = minval(ges_nr_it (2:lat1+1,2:lon1+1,1:nsig)) + endif + if(icf>0) then nc=nc+1 zloc(nvars+nvars1+nc) = minval(ges_cf_it (2:lat1+1,2:lon1+1,1:nsig)) @@ -527,6 +576,15 @@ subroutine prt_guess2(sgrep) nc=nc+1 zloc(2*nvars+nvars1+nc) = maxval(ges_qg_it (2:lat1+1,2:lon1+1,1:nsig)) endif + if(ini>0) then + nc=nc+1 + zloc(2*nvars+nvars1+nc) = maxval(ges_ni_it (2:lat1+1,2:lon1+1,1:nsig)) + endif + if(inr>0) then + nc=nc+1 + zloc(2*nvars+nvars1+nc) = maxval(ges_nr_it (2:lat1+1,2:lon1+1,1:nsig)) + endif + if(icf>0) then nc=nc+1 zloc(2*nvars+nvars1+nc) = maxval(ges_cf_it (2:lat1+1,2:lon1+1,1:nsig)) diff --git a/src/gsi/set_crtm_cloudmod.f90 b/src/gsi/set_crtm_cloudmod.f90 index 9e54bd73ce..b52acad2bd 100644 --- a/src/gsi/set_crtm_cloudmod.f90 +++ b/src/gsi/set_crtm_cloudmod.f90 @@ -187,7 +187,7 @@ subroutine setCloud (cloud_name, icmask, cloud_cont, cloud_efr,jcloud, dp, tp, p cloud(n)%Effective_Radius(:) = cloud_efr(:,n) else !cloud(n)%Effective_Radius(:) = EftSize_(cloud_name(jcloud(n))) - if ( imp_physics==11 .and. lprecip ) then + if (lprecip) then cloud(n)%Effective_Radius(:) = cloud_efr(:,n) else do k = 1, km From e82365d94cac057327de439f01b8aa7793c93100 Mon Sep 17 00:00:00 2001 From: TingLei-NOAA Date: Tue, 6 Aug 2024 11:03:43 -0400 Subject: [PATCH 06/11] Fix incorrect usage of real(i_kind) in mg_input.f90 (#760) --- src/mgbf/mg_input.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mgbf/mg_input.f90 b/src/mgbf/mg_input.f90 index 80b0772c12..ad2a47034d 100644 --- a/src/mgbf/mg_input.f90 +++ b/src/mgbf/mg_input.f90 @@ -57,7 +57,7 @@ subroutine input_2d & integer(i_kind),intent(in):: imax0 integer(i_kind),intent(in):: ampl real(r_kind),dimension(imin:imax,jmin:jmax),intent(out):: V -real(i_kind):: ng,mg,L,m,n +integer(i_kind):: ng,mg,L,m,n !----------------------------------------------------------------------- do m=imin,jmax @@ -134,7 +134,7 @@ subroutine input_3d & integer(i_kind),intent(in):: imax0 integer(i_kind),intent(in):: ampl,incrm real(r_kind),dimension(lmin:lmax,imin:imax,jmin:jmax),intent(out):: V -real(i_kind):: ng,mg,L,m,n +integer(i_kind):: ng,mg,L,m,n !----------------------------------------------------------------------- do l=lmin,lmax From c1eb61c13a21b0b0c0868e42c4d8fd1804922481 Mon Sep 17 00:00:00 2001 From: RussTreadon-NOAA <26926959+RussTreadon-NOAA@users.noreply.github.com> Date: Mon, 12 Aug 2024 08:26:07 -0400 Subject: [PATCH 07/11] add safeguard to thompson_reff (#779) --- src/gsi/crtm_interface.f90 | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/gsi/crtm_interface.f90 b/src/gsi/crtm_interface.f90 index 32a2cdb067..c30a2e06c8 100644 --- a/src/gsi/crtm_interface.f90 +++ b/src/gsi/crtm_interface.f90 @@ -1206,7 +1206,6 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & real(r_kind) :: lai m1=mype+1 - if (mype==0) write(6,*) myname_, ' imp_physics = ', imp_physics if (n_clouds_fwd_wk>0) hwp_guess=zero hwp_total=zero theta_700=zero @@ -2715,15 +2714,17 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) real(r_kind), parameter :: reff_i_min = 2.5_r_kind ! previous value was 10_r_kind real(r_kind), parameter :: reff_i_max = 250.0_r_kind ! previous value was 150_r_kind real(r_kind), parameter:: mu_i = 0.0_r_kind + real(r_kind), parameter:: ni_min = 1.0e-6_r_kind ! minimum number concentration (ccpp-physics) ! Parameters for rain (Lin 1983) real(r_kind), parameter :: rho_r = 1000.0_r_kind ! [kg/m3 ] - real(r_kind), parameter :: reff_r_min = 50.0_r_kind ! [micron] ! previous value was 0.0_r_kind - real(r_kind), parameter :: reff_r_max = 1000.0_r_kind ! [micron] ! previous value was 10000.0_r_kind + real(r_kind), parameter :: reff_r_min = 50.0_r_kind ! [micron] ! previous value was 0.0_r_kind + real(r_kind), parameter :: reff_r_max = 1000.0_r_kind ! [micron] ! previous value was 10000.0_r_kind real(r_kind), parameter:: mu_r = 0.0_r_kind ! Parameters for snow real(r_kind), parameter :: reff_s_min = 5.0_r_kind ! [micron] ! previous value was 0.0_r_kind - real(r_kind), parameter :: reff_s_max = 5000.0_r_kind ! [micron] ! previous value was 10000.0_r_kind + real(r_kind), parameter :: reff_s_max = 5000.0_r_kind ! [micron] ! previous value was 10000.0_r_kind + real(r_kind), parameter:: nr_min = 1.0e-6_r_kind ! minimum number concentration (ccpp-physics) !For snow moments conversions (from Field et al. 2005) real(r_kind), dimension(10), parameter:: & @@ -2769,10 +2770,10 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) reff_max = reff_i_max do k = 1, nsig qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) - if (qx > qmin) then + if (qx > qmin .and. ni(k)>ni_min) then lam_i=exp(1.0_r_kind / 3.0_r_kind * log((am_i*ni(k) *gamma(mu_i + 3.0_r_kind + 1.0_r_kind))/(qx*gamma(mu_i+1.0_r_kind)))) - reff(k) = 0.5_r_kind * (3.0_r_kind /lam_i)*1.0e6_r_kind - reff(k) = max(reff_min, min(reff_max, reff(k))) + reff(k) = 0.5_r_kind * (3.0_r_kind /lam_i)*1.0e6_r_kind + reff(k) = max(reff_min, min(reff_max, reff(k))) else reff(k) = zero endif @@ -2784,7 +2785,7 @@ subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff) reff_max = reff_r_max do k = 1, nsig qx = qxmr(k) * rho_air(k) ! convert mixing ratio (kg/kg) to water content (kg/m3) - if (qx > qmin) then + if (qx > qmin .and. nr(k)>nr_min) then lam_r=exp(1.0_r_kind / 3.0_r_kind * log ((am_r*nr(k) *gamma(mu_r + 3.0_r_kind + 1.0_r_kind))/(qx*gamma(mu_r + 1.0_r_kind)))) reff(k) = 0.5_r_kind *(3.0_r_kind/lam_r)*1.0e6_r_kind reff(k) = max(reff_min, min(reff_max, reff(k))) From 5ff4f3fcdb78f357b5399c89dbdd333b64fe78d9 Mon Sep 17 00:00:00 2001 From: xincjin-NOAA <122945910+xincjin-NOAA@users.noreply.github.com> Date: Thu, 15 Aug 2024 16:47:25 -0400 Subject: [PATCH 08/11] fix a bug in deter_sfc_gmi (#781) --- src/gsi/deter_sfc_mod.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gsi/deter_sfc_mod.f90 b/src/gsi/deter_sfc_mod.f90 index 271e81c5d2..2fa85f8ba5 100644 --- a/src/gsi/deter_sfc_mod.f90 +++ b/src/gsi/deter_sfc_mod.f90 @@ -1410,7 +1410,7 @@ subroutine deter_sfc_gmi(dlat_earth,dlon_earth,isflg) n_grid=int(40000 / grid_dist) + 1 klatn = max(klat1 - n_grid, 1) klonn = klon1 - n_grid - if (klonn < 0) klonn = nlon_sfc - klonn + if (klonn < 1) klonn = nlon_sfc - klonn klatpn = min((klat1 + n_grid), nlat_sfc) klonpn = klon1 + n_grid if (klonpn > nlon_sfc) klonpn = klonpn - nlon_sfc From 8a1c0e725e4a318ea72cd303c7e274f4b46ff52c Mon Sep 17 00:00:00 2001 From: emilyhcliu <36091766+emilyhcliu@users.noreply.github.com> Date: Thu, 22 Aug 2024 14:20:36 -0400 Subject: [PATCH 09/11] Recover the capability of handling model fields from operation gfs.v16.3 (#785) --- src/gsi/general_read_gfsatm.f90 | 177 ++++++++++++++++---------------- 1 file changed, 89 insertions(+), 88 deletions(-) diff --git a/src/gsi/general_read_gfsatm.f90 b/src/gsi/general_read_gfsatm.f90 index c34b269950..e1a5406a13 100755 --- a/src/gsi/general_read_gfsatm.f90 +++ b/src/gsi/general_read_gfsatm.f90 @@ -3588,99 +3588,100 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z if ( icount == icm .or. k==nlevs) then call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) - endif - enddo ! do k=1,nlevs - - do k=1,nlevs - icount=icount+1 - iflag(icount)=15 - ilev(icount)=k - kr = levs+1-k ! netcdf is top to bottom, need to flip + endif + enddo ! do k=1,nlevs - if (mype==mype_use(icount)) then - call read_vardata(filges, 'nccice', rwork3d0, nslice=kr, slicedim=3) - ! cloud ice water number concentration. - if ( diff_res ) then - grid_b=rwork3d0(:,:,1) - vector(1)=.false. - call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) - call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work(kk)=grid2(i,j,1) - enddo - else - grid=rwork3d0(:,:,1) - call general_fill_ns(grd,grid,work) - endif - endif - if ( icount == icm .or. k==nlevs ) then - call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) - endif - enddo ! do k=1,nlevs + ! Read fields specific to Thompson microphysics + if (imp_physics == 8) then + do k=1,nlevs + icount=icount+1 + iflag(icount)=15 + ilev(icount)=k + kr = levs+1-k ! netcdf is top to bottom, need to flip - do k=1,nlevs - icount=icount+1 - iflag(icount)=16 - ilev(icount)=k - kr = levs+1-k ! netcdf is top to bottom, need to flip + if (mype==mype_use(icount)) then + call read_vardata(filges, 'nccice', rwork3d0, nslice=kr, slicedim=3) + ! cloud ice water number concentration. + if ( diff_res ) then + grid_b=rwork3d0(:,:,1) + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork3d0(:,:,1) + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm .or. k==nlevs ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) + endif + enddo ! do k=1,nlevs - if (mype==mype_use(icount)) then - call read_vardata(filges, 'nconrd', rwork3d0, nslice=kr, slicedim=3) - ! rain number concentration. - if ( diff_res ) then - grid_b=rwork3d0(:,:,1) - vector(1)=.false. - call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) - call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) - do kk=1,grd%itotsub - i=grd%ltosi_s(kk) - j=grd%ltosj_s(kk) - work(kk)=grid2(i,j,1) - enddo - else - grid=rwork3d0(:,:,1) - call general_fill_ns(grd,grid,work) - endif - endif - if ( icount == icm .or. k==nlevs) then - call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) - endif - enddo ! do k=1,nlevs + do k=1,nlevs + icount=icount+1 + iflag(icount)=16 + ilev(icount)=k + kr = levs+1-k ! netcdf is top to bottom, need to flip -! do k=1,nlevs -! icount=icount+1 -! iflag(icount)=17 -! ilev(icount)=k -! kr = levs+1-k ! netcdf is top to bottom, need to flip -! -! if (mype==mype_use(icount)) then -! call read_vardata(filges, 'cld_amt', rwork3d0, nslice=kr, slicedim=3) -! ! Cloud amount (cloud fraction). -! if ( diff_res ) then -! grid_b=rwork3d0(:,:,1) -! vector(1)=.false. -! call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) -! call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) -! do kk=1,grd%itotsub -! i=grd%ltosi_s(kk) -! j=grd%ltosj_s(kk) -! work(kk)=grid2(i,j,1) -! enddo -! else -! grid=rwork3d0(:,:,1) -! call general_fill_ns(grd,grid,work) -! endif + if (mype==mype_use(icount)) then + call read_vardata(filges, 'nconrd', rwork3d0, nslice=kr, slicedim=3) + ! rain number concentration. + if ( diff_res ) then + grid_b=rwork3d0(:,:,1) + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork3d0(:,:,1) + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm .or. k==nlevs) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr) + endif + enddo ! do k=1,nlevs + endif ! imp_physics +! do k=1,nlevs +! icount=icount+1 +! iflag(icount)=17 +! ilev(icount)=k +! kr = levs+1-k ! netcdf is top to bottom, need to flip ! -! endif -! if ( icount == icm .or. k==nlevs ) then -! call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & -! g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_cf) -! endif -! enddo ! do k=1,nlevs +! if (mype==mype_use(icount)) then +! call read_vardata(filges, 'cld_amt', rwork3d0, nslice=kr, slicedim=3) +! ! Cloud amount (cloud fraction). +! if ( diff_res ) then +! grid_b=rwork3d0(:,:,1) +! vector(1)=.false. +! call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) +! call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) +! do kk=1,grd%itotsub +! i=grd%ltosi_s(kk) +! j=grd%ltosj_s(kk) +! work(kk)=grid2(i,j,1) +! enddo +! else +! grid=rwork3d0(:,:,1) +! call general_fill_ns(grd,grid,work) +! endif +! endif +! if ( icount == icm .or. k==nlevs ) then +! call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & +! g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_cf) +! endif +! enddo ! do k=1,nlevs else ! read_2m From 8412e3f605a7a60768bc7bddaff3b3b83875fc3d Mon Sep 17 00:00:00 2001 From: TingLei-NOAA Date: Tue, 3 Sep 2024 14:34:00 -0400 Subject: [PATCH 10/11] a quick workaround for increasing the mpi task numbers on orion for ctest :: rrfs_3denvar_rdasens (#788) This PR would provide a quick fix/workaround for rrfs_3denvar_rdasens running on orion by increasing the mpi task numbers running the test. --- regression/regression_param.sh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/regression/regression_param.sh b/regression/regression_param.sh index bfc6f042fc..209762569b 100755 --- a/regression/regression_param.sh +++ b/regression/regression_param.sh @@ -90,8 +90,8 @@ case $regtest in topts[1]="0:05:00" ; popts[1]="40/3/" ; ropts[1]="/1" topts[2]="0:05:00" ; popts[2]="40/5/" ; ropts[2]="/1" elif [[ "$machine" = "Orion" ]]; then - topts[1]="0:15:00" ; popts[1]="5/4/" ; ropts[1]="/1" - topts[2]="0:15:00" ; popts[2]="10/4/" ; ropts[2]="/2" + topts[1]="0:15:00" ; popts[1]="20/6/" ; ropts[1]="/1" + topts[2]="0:15:00" ; popts[2]="20/12/" ; ropts[2]="/2" elif [[ "$machine" = "Hercules" ]]; then topts[1]="0:05:00" ; popts[1]="40/3/" ; ropts[1]="/1" topts[2]="0:05:00" ; popts[2]="40/5/" ; ropts[2]="/2" From 9f44c8798c2087aca06df8f629699632e57df431 Mon Sep 17 00:00:00 2001 From: Innocent Souopgui <162634017+InnocentSouopgui-NOAA@users.noreply.github.com> Date: Fri, 6 Sep 2024 07:47:08 -0500 Subject: [PATCH 11/11] Move to contrib spack-stack on Jet (#787) --- modulefiles/gsi_jet.intel.lua | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/modulefiles/gsi_jet.intel.lua b/modulefiles/gsi_jet.intel.lua index 48189ba241..b210e6efa2 100644 --- a/modulefiles/gsi_jet.intel.lua +++ b/modulefiles/gsi_jet.intel.lua @@ -1,7 +1,7 @@ help([[ ]]) -prepend_path("MODULEPATH", "/mnt/lfs4/HFIP/hfv3gfs/role.epic/spack-stack/spack-stack-1.6.0/envs/gsi-addon-dev-rocky8/install/modulefiles/Core") +prepend_path("MODULEPATH", "/contrib/spack-stack/spack-stack-1.6.0/envs/gsi-addon-intel/install/modulefiles/Core") local python_ver=os.getenv("python_ver") or "3.11.6" local stack_intel_ver=os.getenv("stack_intel_ver") or "2021.5.0" @@ -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", "/mnt/lfs4/HFIP/hfv3gfs/glopara/git/fv3gfs/fix/gsi/20240208") +pushenv("GSI_BINARY_SOURCE_DIR", "/lfs5/HFIP/hfv3gfs/glopara/FIX/fix/gsi/20240208") whatis("Description: GSI environment on Jet with Intel Compilers")