diff --git a/src/gsi/cplr_get_fv3_regional_ensperts.f90 b/src/gsi/cplr_get_fv3_regional_ensperts.f90 index 2645723c62..2382ff1286 100644 --- a/src/gsi/cplr_get_fv3_regional_ensperts.f90 +++ b/src/gsi/cplr_get_fv3_regional_ensperts.f90 @@ -74,7 +74,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) use netcdf , only: nf90_open, nf90_close,nf90_nowrite,nf90_inquire,nf90_format_netcdf4 use netcdf_mod , only: nc_check use gsi_rfv3io_mod, only: fv3lam_io_phymetvars3d_nouv - use obsmod, only: if_model_dbz + use obsmod, only: if_model_dbz,if_model_fed implicit none @@ -85,10 +85,10 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig):: u,v,tv,oz,rh real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2):: ps - real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig)::w,ql,qi,qr,qg,qs,qnr,dbz + real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig)::w,ql,qi,qr,qg,qs,qnr,dbz,fed real(r_kind),dimension(:,:,:),allocatable :: gg_u,gg_v,gg_tv,gg_rh real(r_kind),dimension(:,:,:),allocatable :: gg_w,gg_dbz,gg_qr,gg_qs, & - gg_qi,gg_qg,gg_oz,gg_cwmr + gg_qi,gg_qg,gg_oz,gg_cwmr,gg_fed real(r_kind),dimension(:,:),allocatable :: gg_ps real(r_single),pointer,dimension(:,:,:):: w3 =>NULL() @@ -117,6 +117,8 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) type(type_fv3regfilenameg)::fv3_filename integer(i_kind):: imem_start,n_fv3sar + integer(i_kind):: i_caseflag + if(n_ens/=(n_ens_gfs+n_ens_fv3sar)) then write(6,*)'wrong, the sum of n_ens_gfs and n_ens_fv3sar not equal n_ens, stop' write(6,*)"n_ens, n_ens_gfs and n_ens_fv3sar are",n_ens, n_ens_gfs , n_ens_fv3sar @@ -317,7 +319,6 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) fv3_filename%sfcdata=trim(ensfilenam_str)//"-fv3_sfcdata" fv3_filename%couplerres=trim(ensfilenam_str)//"-coupler.res" - if( mype==iope) then allocate(gg_u(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) allocate(gg_v(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) @@ -325,20 +326,35 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) allocate(gg_rh(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) allocate(gg_oz(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) allocate(gg_ps(grd_ens%nlat,grd_ens%nlon)) - if ( .not. if_model_dbz ) then - call this%general_read_fv3_regional_parallel_over_ens(iope,fv3_filename,gg_ps,gg_u,gg_v,gg_tv,gg_rh,gg_oz) - else + + if ( if_model_dbz .or. if_model_fed ) then allocate(gg_w(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) - allocate(gg_dbz(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) allocate(gg_qr(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) allocate(gg_qs(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) allocate(gg_qi(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) allocate(gg_qg(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) allocate(gg_cwmr(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) - call this%general_read_fv3_regional_parallel_over_ens(iope,fv3_filename,gg_ps,gg_u,gg_v,gg_tv,gg_rh,gg_oz, & - g_ql=gg_cwmr,g_qi=gg_qi,g_qr=gg_qr,g_qs=gg_qs,g_qg=gg_qg,g_w=gg_w,g_dbz=gg_dbz) end if - end if + if ( if_model_dbz) then + allocate(gg_dbz(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) + end if + if ( if_model_fed) then + allocate(gg_fed(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig)) + end if + + if ( if_model_dbz .and. if_model_fed) then + call this%general_read_fv3_regional_parallel_over_ens(iope,fv3_filename,gg_ps,gg_u,gg_v,gg_tv,gg_rh,gg_oz,& + g_ql=gg_cwmr,g_qi=gg_qi,g_qr=gg_qr,g_qs=gg_qs,g_qg=gg_qg,g_w=gg_w,g_dbz=gg_dbz,g_fed=gg_fed) + elseif ( if_model_dbz ) then + call this%general_read_fv3_regional_parallel_over_ens(iope,fv3_filename,gg_ps,gg_u,gg_v,gg_tv,gg_rh,gg_oz,& + g_ql=gg_cwmr,g_qi=gg_qi,g_qr=gg_qr,g_qs=gg_qs,g_qg=gg_qg,g_w=gg_w,g_dbz=gg_dbz) + elseif ( if_model_fed ) then + call this%general_read_fv3_regional_parallel_over_ens(iope,fv3_filename,gg_ps,gg_u,gg_v,gg_tv,gg_rh,gg_oz,& + g_ql=gg_cwmr,g_qi=gg_qi,g_qr=gg_qr,g_qs=gg_qs,g_qg=gg_qg,g_w=gg_w,g_fed=gg_fed) + else + call this%general_read_fv3_regional_parallel_over_ens(iope,fv3_filename,gg_ps,gg_u,gg_v,gg_tv,gg_rh,gg_oz) + end if + end if !mype end do if(mype==0) then write(6,'(I0,A)') mype,': reading ensemble data in parallel is done (parallelization_over_ensmembers=.true.)' @@ -390,47 +406,129 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) endif ! ! READ ENEMBLE MEMBERS DATA + ! + ! There are three options to control the list of variables that + ! will be read in along with the basic variables, ps,u,v,tv,rh,oz. + + ! parallelization_over_ensmembers=.True. only works for cases when l_use_dbz_directDA=.False. + ! Noted that l_use_dbz_directDA and if_modle_dbz couldn't be true at the same time + + ! + ! I_CASEFLAG defination + ! + + ! default: all the three options ( l_use_dbz_directDA, if_model_dbz, if_model_fed) are turned off .i.e., + ! if(.not. (l_use_dbz_directDA .or. if_model_dbz .or. if_model_fed )) + ! read in ps,u,v,tv,rh,oz + i_caseflag=0 + + ! only l_use_dbz_directDA is true + if (l_use_dbz_directDA .and. .not.if_model_dbz .and. .not.if_model_fed) i_caseflag=1 + + ! only if_model_dbz is true + if(.not.l_use_dbz_directDA .and. if_model_dbz .and. .not.if_model_fed) i_caseflag=2 + + ! only if_model_fed is true + if(.not.l_use_dbz_directDA .and. .not.if_model_dbz .and. .not.if_model_fed) i_caseflag=3 + + ! l_use_dbz_directDA=.true. and if_model_fed=.true. + if(l_use_dbz_directDA .and. .not.if_model_dbz .and. if_model_fed) i_caseflag=4 + + ! if_model_dbz=.true. and if_model_fed=.true. + if(.not. l_use_dbz_directDA.and. if_model_dbz .and. if_model_fed) i_caseflag=5 + + + !-------------------------------------------------- + ! When .not. parallelization_over_ensmembers=.True. + ! All the above 6 cases (i_caseflag=0,1,2,3,4,5) are valid in + ! the current applications as of Oct 20 2023. + + !-------------------------------------------- + ! When parallelization_over_ensmembers=.True. + ! Only i_flagcase=0,2,3,5 are vaild choices. + + if( .not. parallelization_over_ensmembers )then if (mype == 0) write(6,'(a,a)') & 'CALL READ_FV3_REGIONAL_ENSPERTS FOR ENS DATA with the filename str : ',trim(ensfilenam_str) - if (.not. (l_use_dbz_directDA .or. if_model_dbz) ) then ! Read additional hydrometers and w for dirZDA - call this%general_read_fv3_regional(fv3_filename,ps,u,v,tv,rh,oz) - else - if( l_use_dbz_directDA ) then - call this%general_read_fv3_regional(fv3_filename,ps,u,v,tv,rh,oz, & - g_ql=ql,g_qi=qi,g_qr=qr,g_qs=qs,g_qg=qg,g_qnr=qnr,g_w=w) - else if( if_model_dbz )then - call this%general_read_fv3_regional(fv3_filename,ps,u,v,tv,rh,oz, & - g_ql=ql,g_qi=qi,g_qr=qr,g_qs=qs,g_qg=qg,g_qnr=qnr,g_w=w,g_dbz=dbz) - end if - end if + + select case (i_caseflag) + case (0) + call this%general_read_fv3_regional(fv3_filename,ps,u,v,tv,rh,oz) + case (1) + call this%general_read_fv3_regional(fv3_filename,ps,u,v,tv,rh,oz, & + g_ql=ql,g_qi=qi,g_qr=qr,g_qs=qs,g_qg=qg,g_qnr=qnr,g_w=w) + case (2) + call this%general_read_fv3_regional(fv3_filename,ps,u,v,tv,rh,oz, & + g_ql=ql,g_qi=qi,g_qr=qr,g_qs=qs,g_qg=qg,g_qnr=qnr,g_w=w,g_dbz=dbz) + case (3) + call this%general_read_fv3_regional(fv3_filename,ps,u,v,tv,rh,oz, & + g_ql=ql,g_qi=qi,g_qr=qr,g_qs=qs,g_qg=qg,g_w=w,g_fed=fed) + case (4) + call this%general_read_fv3_regional(fv3_filename,ps,u,v,tv,rh,oz, & + g_ql=ql,g_qi=qi,g_qr=qr,g_qs=qs,g_qg=qg,g_qnr=qnr,g_w=w,g_fed=fed) + case (5) + call this%general_read_fv3_regional(fv3_filename,ps,u,v,tv,rh,oz, & + g_ql=ql,g_qi=qi,g_qr=qr,g_qs=qs,g_qg=qg,g_qnr=qnr,g_w=w,g_dbz=dbz,g_fed=fed) + end select end if - + if( parallelization_over_ensmembers )then iope=(n_fv3sar-1)*npe/n_ens_fv3sar if(mype==iope) then write(0,'(I0,A,I0,A)') mype,': scatter member ',n_fv3sar,' to other ranks...' - if( if_model_dbz )then - call this%parallel_read_fv3_step2(mype,iope,& + select case (i_caseflag) + case (0) + call this%parallel_read_fv3_step2(mype,iope,& + g_ps=ps,g_u=u,g_v=v,g_tv=tv,g_rh=rh,g_ql=ql,g_oz=oz, & + gg_ps=gg_ps,gg_tv=gg_tv,gg_u=gg_u,gg_v=gg_v,gg_rh=gg_rh) + case (2) + call this%parallel_read_fv3_step2(mype,iope,& g_ps=ps,g_u=u,g_v=v,g_tv=tv,g_rh=rh,g_ql=ql,& g_oz=oz,g_w=w,g_qr=qr,g_qs=qs,g_qi=qi,g_qg=qg,g_dbz=dbz,& gg_ps=gg_ps,gg_tv=gg_tv,gg_u=gg_u,gg_v=gg_v,& gg_rh=gg_rh,gg_w=gg_w,gg_dbz=gg_dbz,gg_qr=gg_qr,& gg_qs=gg_qs,gg_qi=gg_qi,gg_qg=gg_qg,gg_ql=gg_cwmr) - else - call this%parallel_read_fv3_step2(mype,iope,& - g_ps=ps,g_u=u,g_v=v,g_tv=tv,g_rh=rh,g_ql=ql,g_oz=oz, & - gg_ps=gg_ps,gg_tv=gg_tv,gg_u=gg_u,gg_v=gg_v,gg_rh=gg_rh) - end if + case (3) + call this%parallel_read_fv3_step2(mype,iope,& + g_ps=ps,g_u=u,g_v=v,g_tv=tv,g_rh=rh,g_ql=ql,& + g_oz=oz,g_w=w,g_qr=qr,g_qs=qs,g_qi=qi,g_qg=qg,g_fed=fed,& + gg_ps=gg_ps,gg_tv=gg_tv,gg_u=gg_u,gg_v=gg_v,& + gg_rh=gg_rh,gg_w=gg_w,gg_fed=gg_fed,gg_qr=gg_qr,& + gg_qs=gg_qs,gg_qi=gg_qi,gg_qg=gg_qg,gg_ql=gg_cwmr) + case (5) + call this%parallel_read_fv3_step2(mype,iope,& + g_ps=ps,g_u=u,g_v=v,g_tv=tv,g_rh=rh,g_ql=ql,& + g_oz=oz,g_w=w,g_qr=qr,g_qs=qs,g_qi=qi,g_qg=qg,g_dbz=dbz,g_fed=fed,& + gg_ps=gg_ps,gg_tv=gg_tv,gg_u=gg_u,gg_v=gg_v,& + gg_rh=gg_rh,gg_w=gg_w,gg_dbz=gg_dbz,gg_fed=gg_fed,gg_qr=gg_qr,& + gg_qs=gg_qs,gg_qi=gg_qi,gg_qg=gg_qg,gg_ql=gg_cwmr) + case (1,4) + write(6,*)'i_case_flag=1 or 4 is not a valid choice for parallelization_over_ensmembers=.T. Stop(8880) ' + call stop2(8880) + end select else - if( if_model_dbz )then - call this%parallel_read_fv3_step2(mype,iope,& + select case (i_caseflag) + case (0) + call this%parallel_read_fv3_step2(mype,iope,& + g_ps=ps,g_u=u,g_v=v,g_tv=tv,g_rh=rh,g_ql=ql,g_oz=oz) + case (2) + call this%parallel_read_fv3_step2(mype,iope,& g_ps=ps,g_u=u,g_v=v,g_tv=tv,g_rh=rh,g_ql=ql,& g_oz=oz,g_w=w,g_qr=qr,g_qs=qs,g_qi=qi,g_qg=qg,g_dbz=dbz) - else - call this%parallel_read_fv3_step2(mype,iope,& - g_ps=ps,g_u=u,g_v=v,g_tv=tv,g_rh=rh,g_ql=ql,g_oz=oz) - endif + case (3) + call this%parallel_read_fv3_step2(mype,iope,& + g_ps=ps,g_u=u,g_v=v,g_tv=tv,g_rh=rh,g_ql=ql,& + g_oz=oz,g_w=w,g_qr=qr,g_qs=qs,g_qi=qi,g_qg=qg,g_fed=fed) + case (5) + call this%parallel_read_fv3_step2(mype,iope,& + g_ps=ps,g_u=u,g_v=v,g_tv=tv,g_rh=rh,g_ql=ql,& + g_oz=oz,g_w=w,g_qr=qr,g_qs=qs,g_qi=qi,g_qg=qg,g_dbz=dbz,g_fed=fed) + case (1,4) + write(6,*)'i_case_flag=1 or 4 is not a valid choice for parallelization_over_ensmembers=.T. Stop(8880) ' + call stop2(8880) + end select + endif call MPI_Barrier(mpi_comm_world,ierror) @@ -601,6 +699,16 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) end do end do end do + + case('fed','FED') + do k=1,grd_ens%nsig + do i=1,grd_ens%lon2 + do j=1,grd_ens%lat2 + w3(j,i,k) = fed(j,i,k) + x3(j,i,k)=x3(j,i,k)+fed(j,i,k) + end do + end do + end do end select @@ -709,7 +817,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) end subroutine get_fv3_regional_ensperts_run subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g_rh,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,g_qnr,g_w,g_dbz) + g_ql,g_qi,g_qr,g_qs,g_qg,g_qnr,g_w,g_dbz,g_fed) !$$$ subprogram documentation block ! first compied from general_read_arw_regional . . . . ! subprogram: general_read_fv3_regional read fv3sar model ensemble members @@ -760,7 +868,7 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g use hybrid_ensemble_parameters, only: grd_ens use directDA_radaruse_mod, only: l_use_cvpqx, cvpqx_pval, cld_nt_updt use directDA_radaruse_mod, only: l_cvpnr, cvpnr_pval - use obsmod, only:if_model_dbz + use obsmod, only:if_model_dbz,if_model_fed implicit none @@ -769,7 +877,7 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g class(get_fv3_regional_ensperts_class), intent(inout) :: this type (type_fv3regfilenameg) , intent (in) :: fv3_filenameginput real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),intent(out)::g_u,g_v,g_tv,g_rh,g_oz - real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),optional,intent(out)::g_ql,g_qi,g_qr,g_dbz + real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),optional,intent(out)::g_ql,g_qi,g_qr,g_dbz,g_fed real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),optional,intent(out)::g_qs,g_qg,g_qnr,g_w real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2),intent(out):: g_ps @@ -857,7 +965,7 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g fv3_filenameginput%dynvars,fv3_filenameginput) call gsi_fv3ncdf_read(grd_fv3lam_ens_tracer_io_nouv,gsibundle_fv3lam_ens_tracer_nouv,& fv3_filenameginput%tracers,fv3_filenameginput) - if( if_model_dbz ) then + if( if_model_dbz .or. if_model_fed ) then call gsi_fv3ncdf_read(grd_fv3lam_ens_phyvar_io_nouv,gsibundle_fv3lam_ens_phyvar_nouv,& fv3_filenameginput%phyvars,fv3_filenameginput) end if @@ -869,9 +977,9 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g endif ier=0 call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_dynvar_nouv, 'tsen' ,g_tsen ,istatus );ier=ier+istatus - call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'q' ,g_q ,istatus );ier=ier+istatus - call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'oz' ,g_oz ,istatus );ier=ier+istatus - if (l_use_dbz_directDA .or. if_model_dbz) then + call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'q' ,g_q ,istatus );ier=ier+istatus + call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'oz' ,g_oz ,istatus );ier=ier+istatus + if (l_use_dbz_directDA .or. if_model_dbz .or. if_model_fed) then call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'ql' ,g_ql ,istatus );ier=ier+istatus call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'qi' ,g_qi ,istatus );ier=ier+istatus call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'qr' ,g_qr ,istatus );ier=ier+istatus @@ -882,6 +990,8 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_dynvar_nouv, 'w' , g_w ,istatus );ier=ier+istatus if( if_model_dbz )& call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_phyvar_nouv, 'dbz' , g_dbz ,istatus );ier=ier+istatus + if( if_model_fed )& + call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_phyvar_nouv, 'fed' , g_fed, istatus );ier=ier+istatus end if @@ -990,7 +1100,7 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g end subroutine general_read_fv3_regional subroutine general_read_fv3_regional_parallel_over_ens(this,iope,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g_rh,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,g_qnr,g_w,g_dbz) + g_ql,g_qi,g_qr,g_qs,g_qg,g_qnr,g_w,g_dbz,g_fed) !$$$ subprogram documentation block ! first compied from general_read_arw_regional . . . . ! subprogram: general_read_fv3_regional read fv3sar model ensemble members @@ -1039,7 +1149,7 @@ subroutine general_read_fv3_regional_parallel_over_ens(this,iope,fv3_filenamegin use gsi_bundlemod, only: gsi_grid use gsi_bundlemod, only: gsi_bundlecreate,gsi_bundledestroy use gsi_bundlemod, only: gsi_bundlegetvar - use obsmod, only: if_model_dbz + use obsmod, only: if_model_dbz,if_model_fed use gsi_rfv3io_mod, only: gsi_fv3ncdf_read_ens_parallel_over_ens,gsi_fv3ncdf_readuv_ens_parallel_over_ens @@ -1051,7 +1161,7 @@ subroutine general_read_fv3_regional_parallel_over_ens(this,iope,fv3_filenamegin integer(i_kind), intent (in) :: iope type (type_fv3regfilenameg) , intent (in) :: fv3_filenameginput real(r_kind),dimension(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig),intent(out)::g_u,g_v,g_tv,g_rh,g_oz - real(r_kind),dimension(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig),optional,intent(out)::g_ql,g_qi,g_qr,g_dbz + real(r_kind),dimension(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig),optional,intent(out)::g_ql,g_qi,g_qr,g_dbz,g_fed real(r_kind),dimension(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig),optional,intent(out)::g_qs,g_qg,g_qnr,g_w real(r_kind),dimension(grd_ens%nlat,grd_ens%nlon),intent(out):: g_ps @@ -1102,11 +1212,17 @@ subroutine general_read_fv3_regional_parallel_over_ens(this,iope,fv3_filenamegin endif if(fv3sar_ensemble_opt == 0) then - if (if_model_dbz) then + if (if_model_dbz .or. if_model_fed) then call gsi_fv3ncdf_read_ens_parallel_over_ens(fv3_filenameginput%dynvars,fv3_filenameginput,delp=g_delp,tsen=g_tsen,w=g_w,iope=iope) call gsi_fv3ncdf_read_ens_parallel_over_ens(fv3_filenameginput%tracers,fv3_filenameginput,q=g_q,oz=g_oz,ql=g_ql,qr=g_qr,& qs=g_qs,qi=g_qi,qg=g_qg,iope=iope) - call gsi_fv3ncdf_read_ens_parallel_over_ens(fv3_filenameginput%phyvars,fv3_filenameginput,dbz=g_dbz,iope=iope) + if(if_model_dbz .and. if_model_fed) then + call gsi_fv3ncdf_read_ens_parallel_over_ens(fv3_filenameginput%phyvars,fv3_filenameginput,dbz=g_dbz,fed=g_fed,iope=iope) + elseif(if_model_dbz) then + call gsi_fv3ncdf_read_ens_parallel_over_ens(fv3_filenameginput%phyvars,fv3_filenameginput,dbz=g_dbz,iope=iope) + elseif(if_model_fed) then + call gsi_fv3ncdf_read_ens_parallel_over_ens(fv3_filenameginput%phyvars,fv3_filenameginput,fed=g_fed,iope=iope) + end if else call gsi_fv3ncdf_read_ens_parallel_over_ens(fv3_filenameginput%dynvars,fv3_filenameginput,delp=g_delp,tsen=g_tsen,iope=iope) call gsi_fv3ncdf_read_ens_parallel_over_ens(fv3_filenameginput%tracers,fv3_filenameginput,q=g_q,oz=g_oz,iope=iope) @@ -1169,8 +1285,8 @@ end subroutine general_read_fv3_regional_parallel_over_ens subroutine parallel_read_fv3_step2(this,mype,iope, & g_ps,g_u,g_v,g_tv,g_rh,g_ql,g_oz,g_w,g_qr,g_qs,g_qi,& - g_qg,g_dbz, & - gg_ps,gg_tv,gg_u,gg_v,gg_rh,gg_w,gg_dbz,gg_qr,& + g_qg,g_dbz,g_fed, & + gg_ps,gg_tv,gg_u,gg_v,gg_rh,gg_w,gg_dbz,gg_fed,gg_qr,& gg_qs,gg_qi,gg_qg,gg_ql) !$$$ subprogram documentation block @@ -1210,7 +1326,7 @@ subroutine parallel_read_fv3_step2(this,mype,iope, & real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),intent(out):: & g_u,g_v,g_tv,g_rh,g_ql,g_oz real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),intent(out),optional::& - g_w,g_qr,g_qs,g_qi,g_qg,g_dbz + g_w,g_qr,g_qs,g_qi,g_qg,g_dbz,g_fed integer(i_kind), intent(in) :: mype, iope real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2),intent(out):: g_ps @@ -1219,7 +1335,7 @@ subroutine parallel_read_fv3_step2(this,mype,iope, & gg_u,gg_v,gg_tv,gg_rh real(r_kind),optional,dimension(grd_ens%nlat,grd_ens%nlon,grd_ens%nsig) :: & - gg_w,gg_dbz,gg_qr,gg_qs,gg_qi,gg_qg,gg_ql + gg_w,gg_dbz,gg_fed,gg_qr,gg_qs,gg_qi,gg_qg,gg_ql real(r_kind),optional,dimension(grd_ens%nlat,grd_ens%nlon):: gg_ps ! Declare local variables @@ -1250,13 +1366,10 @@ subroutine parallel_read_fv3_step2(this,mype,iope, & if (mype==iope) call this%fill_regional_2d(gg_rh(1,1,k),wrk_send_2d) call mpi_scatterv(wrk_send_2d,grd_ens%ijn_s,grd_ens%displs_s,mpi_rtype, & g_rh(1,1,k),grd_ens%ijn_s(mype+1),mpi_rtype,iope,mpi_comm_world,ierror) - if( present(g_dbz) )then + if( present(g_dbz) .or. present(g_fed) )then if (mype==iope) call this%fill_regional_2d(gg_w(1,1,k),wrk_send_2d) call mpi_scatterv(wrk_send_2d,grd_ens%ijn_s,grd_ens%displs_s,mpi_rtype, & g_w(1,1,k),grd_ens%ijn_s(mype+1),mpi_rtype,iope,mpi_comm_world,ierror) - if (mype==iope) call this%fill_regional_2d(gg_dbz(1,1,k),wrk_send_2d) - call mpi_scatterv(wrk_send_2d,grd_ens%ijn_s,grd_ens%displs_s,mpi_rtype,& - g_dbz(1,1,k),grd_ens%ijn_s(mype+1),mpi_rtype,iope,mpi_comm_world,ierror) if (mype==iope) call this%fill_regional_2d(gg_qr(1,1,k),wrk_send_2d) call mpi_scatterv(wrk_send_2d,grd_ens%ijn_s,grd_ens%displs_s,mpi_rtype,& g_qr(1,1,k),grd_ens%ijn_s(mype+1),mpi_rtype,iope,mpi_comm_world,ierror) @@ -1272,6 +1385,16 @@ subroutine parallel_read_fv3_step2(this,mype,iope, & if (mype==iope) call this%fill_regional_2d(gg_ql(1,1,k),wrk_send_2d) call mpi_scatterv(wrk_send_2d,grd_ens%ijn_s,grd_ens%displs_s,mpi_rtype,& g_ql(1,1,k),grd_ens%ijn_s(mype+1),mpi_rtype,iope,mpi_comm_world,ierror) + if( present(g_dbz)) then + if (mype==iope) call this%fill_regional_2d(gg_dbz(1,1,k),wrk_send_2d) + call mpi_scatterv(wrk_send_2d,grd_ens%ijn_s,grd_ens%displs_s,mpi_rtype,& + g_dbz(1,1,k),grd_ens%ijn_s(mype+1),mpi_rtype,iope,mpi_comm_world,ierror) + end if + if( present(g_fed)) then + if (mype==iope) call this%fill_regional_2d(gg_fed(1,1,k),wrk_send_2d) + call mpi_scatterv(wrk_send_2d,grd_ens%ijn_s,grd_ens%displs_s,mpi_rtype,& + g_fed(1,1,k),grd_ens%ijn_s(mype+1),mpi_rtype,iope,mpi_comm_world,ierror) + end if end if enddo deallocate(wrk_send_2d) diff --git a/src/gsi/gsi_fedOper.F90 b/src/gsi/gsi_fedOper.F90 index b2b2400ff0..e704a1c056 100644 --- a/src/gsi/gsi_fedOper.F90 +++ b/src/gsi/gsi_fedOper.F90 @@ -9,6 +9,7 @@ module gsi_fedOper ! 2023-07-10 D. Dowell - created new module for FED (flash extent ! density); gsi_dbzOper.F90 code used as a ! starting point for developing this new module +! 2023-08-24 H. Wang - Turned on intfed and stpfed ! ! input argument list: see Fortran 90 style document below ! @@ -128,6 +129,7 @@ subroutine setup_(self, lunin, mype, is, nobs, init_pass, last_pass) end subroutine setup_ subroutine intjo1_(self, ibin, rval,sval, qpred,sbias) + use intfedmod, only: intjo => intfed use gsi_bundlemod , only: gsi_bundle use bias_predictors, only: predictors use m_obsNode , only: obsNode @@ -145,9 +147,14 @@ subroutine intjo1_(self, ibin, rval,sval, qpred,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 stpfedmod, only: stpjo => stpfed use gsi_bundlemod, only: gsi_bundle use bias_predictors, only: predictors use m_obsNode , only: obsNode @@ -169,6 +176,9 @@ subroutine stpjo1_(self, ibin, dval,xval,pbcjo,sges,nstep,dbias,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_fedOper diff --git a/src/gsi/gsi_files.cmake b/src/gsi/gsi_files.cmake index b514e11c1e..5a7d29c208 100644 --- a/src/gsi/gsi_files.cmake +++ b/src/gsi/gsi_files.cmake @@ -274,6 +274,7 @@ intaod.f90 intcldch.f90 intco.f90 intdbz.f90 +intfed.f90 intdw.f90 intgps.f90 intgust.f90 @@ -594,6 +595,7 @@ stpcalc.f90 stpcldch.f90 stpco.f90 stpdbz.f90 +stpfed.f90 stpdw.f90 stpgps.f90 stpgust.f90 diff --git a/src/gsi/gsi_rfv3io_mod.f90 b/src/gsi/gsi_rfv3io_mod.f90 index d0cbd3afbd..e62cc06f2b 100644 --- a/src/gsi/gsi_rfv3io_mod.f90 +++ b/src/gsi/gsi_rfv3io_mod.f90 @@ -97,7 +97,7 @@ module gsi_rfv3io_mod type(sub2grid_info) :: grd_fv3lam_tracersmoke_ionouv type(sub2grid_info) :: grd_fv3lam_phyvar_ionouv type(sub2grid_info) :: grd_fv3lam_uv - integer(i_kind) ,parameter:: ndynvarslist=13, ntracerslist=8, nphyvarslist=1 + integer(i_kind) ,parameter:: ndynvarslist=13, ntracerslist=8, nphyvarslist=2 character(len=max_varname_length), dimension(ndynvarslist), parameter :: & vardynvars = [character(len=max_varname_length) :: & @@ -106,13 +106,13 @@ module gsi_rfv3io_mod vartracers = [character(len=max_varname_length) :: & 'q','oz','ql','qi','qr','qs','qg','qnr',aeronames_cmaq_fv3,'pm25at','pm25ac','pm25co','pm2_5','amassi','amassj','amassk',aeronames_smoke_fv3] character(len=max_varname_length), dimension(nphyvarslist), parameter :: & - varphyvars = [character(len=max_varname_length) :: 'dbz'] - character(len=max_varname_length), dimension(16+naero_cmaq_fv3+7+naero_smoke_fv3), parameter :: & + varphyvars = [character(len=max_varname_length) :: 'dbz','fed'] + character(len=max_varname_length), dimension(16+naero_cmaq_fv3+7+naero_smoke_fv3+1), parameter :: & varfv3name = [character(len=max_varname_length) :: & - 'u','v','W','T','delp','sphum','o3mr','liq_wat','ice_wat','rainwat','snowwat','graupel','rain_nc','ref_f3d','ps','DZ', & + 'u','v','W','T','delp','sphum','o3mr','liq_wat','ice_wat','rainwat','snowwat','graupel','rain_nc','ref_f3d','flash_extent_density','ps','DZ', & aeronames_cmaq_fv3,'pm25at','pm25ac','pm25co','pm2_5','amassi','amassj','amassk',aeronames_smoke_fv3], & vgsiname = [character(len=max_varname_length) :: & - 'u','v','w','tsen','delp','q','oz','ql','qi','qr','qs','qg','qnr','dbz','ps','delzinc', & + 'u','v','w','tsen','delp','q','oz','ql','qi','qr','qs','qg','qnr','dbz','fed','ps','delzinc', & aeronames_cmaq_fv3,'pm25at','pm25ac','pm25co','pm2_5','amassi','amassj','amassk',aeronames_smoke_fv3] character(len=max_varname_length),dimension(:),allocatable:: name_metvars2d character(len=max_varname_length),dimension(:),allocatable:: name_metvars3d @@ -801,7 +801,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) use gsi_metguess_mod, only: gsi_metguess_get use netcdf, only:nf90_open,nf90_close,nf90_inquire,nf90_nowrite, nf90_format_netcdf4 use gsi_chemguess_mod, only: gsi_chemguess_get - use obsmod, only: if_model_dbz + use obsmod, only: if_model_dbz,if_model_fed implicit none @@ -809,7 +809,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) integer(i_kind) :: it character(len=24),parameter :: myname = 'read_fv3_netcdf_guess' integer(i_kind) k,i,j - integer(i_kind) ier,istatus + integer(i_kind) ier,istatus,ivar real(r_kind),dimension(:,:),pointer::ges_ps=>NULL() real(r_kind),dimension(:,:),pointer::ges_ps_readin=>NULL() real(r_kind),dimension(:,:),pointer::ges_z=>NULL() @@ -833,6 +833,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) real(r_kind),dimension(:,:,:),pointer::ges_qnr=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_w=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_dbz=>NULL() + real(r_kind),dimension(:,:,:),pointer::ges_fed=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_aalj=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_acaj=>NULL() @@ -1016,12 +1017,13 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) write(6,*)"the set up for met variable is not as expected, abort" call stop2(222) endif - if ( if_model_dbz ) then - if( nphyvario3d<=0 ) then - write(6,*)"the set up for met variable (phyvar) is not as expected, abort" - call stop2(223) - end if - endif + + ivar=0 ; if (if_model_dbz) ivar=ivar+1; if(if_model_fed) ivar=ivar+1 + if ( ivar > nphyvario3d ) then + write(6,*)"the set up for met variable (dbz and fed in phyvar) is not as expected,abort" + call stop2(223) + end if + if (fv3sar_bg_opt == 0.and.ifindstrloc(name_metvars3d,'delp') <= 0)then ndynvario3d=ndynvario3d+1 ! for delp endif @@ -1217,7 +1219,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) ntracerio2d=0 endif - if( if_model_dbz )then + if( allocated(fv3lam_io_phymetvars3d_nouv) )then call gsi_bundlecreate(gsibundle_fv3lam_phyvar_nouv,GSI_MetGuess_Bundle(it)%grid,'gsibundle_fv3lam_phyvar_nouv',istatus, & names3d=fv3lam_io_phymetvars3d_nouv) end if @@ -1311,7 +1313,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) endif - if ( if_model_dbz )then + if ( if_model_dbz .or. if_model_fed )then inner_vars=1 numfields=inner_vars*(nphyvario3d*grd_a%nsig) deallocate(lnames,names) @@ -1352,7 +1354,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'tv' ,ges_tv ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'q' ,ges_q ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'oz' ,ges_oz ,istatus );ier=ier+istatus - if (l_use_dbz_directDA .or. if_model_dbz) then + if (l_use_dbz_directDA .or. nphyvario3d > 0) then call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'ql' ,ges_ql ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qi' ,ges_qi ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qr' ,ges_qr ,istatus );ier=ier+istatus @@ -1365,6 +1367,8 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) end if if(if_model_dbz) & call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'dbz' , ges_dbz ,istatus );ier=ier+istatus + if(if_model_fed) & + call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'fed' , ges_fed ,istatus );ier=ier+istatus end if if (ier/=0) call die(trim(myname),'cannot get pointers for fv3 met-fields, ier =',ier) @@ -1427,7 +1431,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) & ,fv3filenamegin(it)%dynvars,fv3filenamegin(it)) call gsi_fv3ncdf_read(grd_fv3lam_tracer_ionouv,gsibundle_fv3lam_tracer_nouv & & ,fv3filenamegin(it)%tracers,fv3filenamegin(it)) - if( if_model_dbz )then + if( nphyvario3d > 0 )then call gsi_fv3ncdf_read(grd_fv3lam_phyvar_ionouv,gsibundle_fv3lam_phyvar_nouv & & ,fv3filenamegin(it)%phyvars,fv3filenamegin(it)) end if @@ -1526,8 +1530,9 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) if (laeroana_fv3smoke) then call gsi_copy_bundle(gsibundle_fv3lam_tracersmoke_nouv,GSI_ChemGuess_Bundle(it)) endif - - if(if_model_dbz) call gsi_copy_bundle(gsibundle_fv3lam_phyvar_nouv,GSI_MetGuess_Bundle(it)) + if ( nphyvario3d > 0 ) then + call gsi_copy_bundle(gsibundle_fv3lam_phyvar_nouv,GSI_MetGuess_Bundle(it)) + end if call GSI_BundleGetPointer ( gsibundle_fv3lam_dynvar_nouv, 'tsen' ,ges_tsen_readin ,istatus );ier=ier+istatus !! tsen2tv !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do k=1,nsig @@ -2354,7 +2359,7 @@ subroutine gsi_fv3ncdf_read(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin) countloc=(/nxcase,nycase,1/) ! Variable ref_f3d in phy_data.nc has a smaller domain size than ! dynvariables and tracers as well as a reversed order in vertical - if ( trim(adjustl(varname)) == 'ref_f3d' )then + if ( trim(adjustl(varname)) == 'ref_f3d' .or. trim(adjustl(varname)) == 'flash_extent_density' )then iret=nf90_inquire_dimension(gfile_loc,1,name,len) if(trim(name)=='xaxis_1') nx_phy=len if( nx_phy == nxcase )then @@ -2380,7 +2385,7 @@ subroutine gsi_fv3ncdf_read(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin) enddo else iret=nf90_inq_varid(gfile_loc,trim(adjustl(varname)),var_id) - if ( trim(adjustl(varname)) == 'ref_f3d' )then + if ( trim(adjustl(varname)) == 'ref_f3d'.or. trim(adjustl(varname)) == 'flash_extent_density' )then uu2d = 0.0_r_kind iret=nf90_get_var(gfile_loc,var_id,uu2d_tmp,start=startloc_tmp,count=countloc_tmp) where(uu2d_tmp < 0.0_r_kind) @@ -2845,7 +2850,7 @@ subroutine gsi_fv3ncdf_readuv_v1(grd_uv,ges_u,ges_v,fv3filenamegin) end subroutine gsi_fv3ncdf_readuv_v1 subroutine gsi_fv3ncdf_read_ens_parallel_over_ens(filenamein,fv3filenamegin, & - delp,tsen,w,q,oz,ql,qr,qs,qi,qg,dbz,iope) + delp,tsen,w,q,oz,ql,qr,qs,qi,qg,dbz,fed,iope) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_fv3ncdf_read_ens_parallel_over_ens @@ -2887,7 +2892,7 @@ subroutine gsi_fv3ncdf_read_ens_parallel_over_ens(filenamein,fv3filenamegin, & integer(i_kind) ,intent(in ) :: iope real(r_kind),allocatable,dimension(:,:):: uu2d, uu2d_tmp real(r_kind),dimension(nlat,nlon,nsig):: hwork - real(r_kind),dimension(nlat,nlon,nsig),intent(out),optional:: delp,tsen,w,q,oz,ql,qr,qs,qi,qg,dbz + real(r_kind),dimension(nlat,nlon,nsig),intent(out),optional:: delp,tsen,w,q,oz,ql,qr,qs,qi,qg,dbz,fed character(len=max_varname_length) :: varname character(len=max_varname_length) :: name character(len=max_filename_length), allocatable,dimension(:) :: varname_files @@ -2932,11 +2937,18 @@ subroutine gsi_fv3ncdf_read_ens_parallel_over_ens(filenamein,fv3filenamegin, & varname_files = (/'sphum',' o3mr'/) end if end if - if( present(dbz) )then ! phyvars: dbz + if( present(dbz) .and. present(fed) )then ! phyvars: dbz, fed + allocate(varname_files(2)) + varname_files = (/'ref_f3d ','flash_extent_density'/) + elseif( present(dbz) )then ! phyvars: dbz allocate(varname_files(1)) varname_files = (/'ref_f3d'/) + elseif( present(fed) )then ! phyvars: fed + allocate(varname_files(1)) + varname_files = (/'flash_extent_density'/) end if + if(fv3_io_layout_y > 1) then allocate(gfile_loc_layout(0:fv3_io_layout_y-1)) do nio=0,fv3_io_layout_y-1 @@ -2967,7 +2979,7 @@ subroutine gsi_fv3ncdf_read_ens_parallel_over_ens(filenamein,fv3filenamegin, & varname = trim(varname_files(ivar)) ! Variable ref_f3d in phy_data.nc has a smaller domain size than ! dynvariables and tracers as well as a reversed order in vertical - if ( trim(adjustl(varname)) == 'ref_f3d' )then + if ( trim(adjustl(varname)) == 'ref_f3d' .or. trim(adjustl(varname)) == 'flash_extent_density' )then iret=nf90_inquire_dimension(gfile_loc,1,name,len) if(trim(name)=='xaxis_1') nx_phy=len if( nx_phy == nxcase )then @@ -2993,7 +3005,7 @@ subroutine gsi_fv3ncdf_read_ens_parallel_over_ens(filenamein,fv3filenamegin, & enddo else iret=nf90_inq_varid(gfile_loc,trim(adjustl(varname)),var_id) - if ( trim(adjustl(varname)) == 'ref_f3d' )then + if ( trim(adjustl(varname)) == 'ref_f3d' .or. trim(adjustl(varname)) == 'flash_extent_density' )then uu2d = 0.0_r_kind iret=nf90_get_var(gfile_loc,var_id,uu2d_tmp,start=startloc_tmp,count=countloc_tmp) where(uu2d_tmp < 0.0_r_kind) @@ -3041,8 +3053,13 @@ subroutine gsi_fv3ncdf_read_ens_parallel_over_ens(filenamein,fv3filenamegin, & end if end if end if - if( present(dbz) )then ! phyvars: dbz + if( present(dbz) .and. present(fed) )then ! phyvars: dbz,fed + if(ivar == 1) dbz = hwork + if(ivar == 2) fed = hwork + elseif( present(dbz) )then ! phyvars: dbz dbz = hwork + elseif( present(fed) )then ! phyvars: fed + fed = hwork end if end do @@ -3271,7 +3288,7 @@ subroutine wrfv3_netcdf(fv3filenamegin) use directDA_radaruse_mod, only: l_cvpnr, cvpnr_pval use gridmod, only: eta1_ll,eta2_ll use constants, only: one - use obsmod, only: if_model_dbz + use obsmod, only: if_model_dbz,if_model_fed implicit none @@ -3299,6 +3316,7 @@ subroutine wrfv3_netcdf(fv3filenamegin) real(r_kind),pointer,dimension(:,:,:):: ges_qnr =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_w =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_dbz =>NULL() + real(r_kind),pointer,dimension(:,:,:):: ges_fed =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_delzinc =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_delp =>NULL() real(r_kind),dimension(:,: ),allocatable:: ges_ps_write @@ -3388,7 +3406,7 @@ subroutine wrfv3_netcdf(fv3filenamegin) call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'u' , ges_u ,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'v' , ges_v ,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'q' ,ges_q ,istatus);ier=ier+istatus - if (l_use_dbz_directDA .or. if_model_dbz) then + if (l_use_dbz_directDA .or. if_model_dbz .or.if_model_fed) then call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'ql' ,ges_ql ,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qi' ,ges_qi ,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qr' ,ges_qr ,istatus);ier=ier+istatus @@ -3399,6 +3417,8 @@ subroutine wrfv3_netcdf(fv3filenamegin) call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'w' , ges_w ,istatus);ier=ier+istatus if( if_model_dbz )& call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'dbz' , ges_dbz ,istatus);ier=ier+istatus + if( if_model_fed )& + call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'fed' , ges_fed ,istatus);ier=ier+istatus end if if(i_use_2mq4b > 0 .and. i_use_2mt4b > 0 ) then call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'q2m',ges_q2m,istatus); ier=ier+istatus @@ -3532,7 +3552,7 @@ subroutine wrfv3_netcdf(fv3filenamegin) call gsi_copy_bundle(GSI_MetGuess_Bundle(it),gsibundle_fv3lam_dynvar_nouv) call gsi_copy_bundle(GSI_MetGuess_Bundle(it),gsibundle_fv3lam_tracer_nouv) - if( if_model_dbz ) call gsi_copy_bundle(GSI_MetGuess_Bundle(it),gsibundle_fv3lam_phyvar_nouv) + if( if_model_dbz .or. if_model_fed ) call gsi_copy_bundle(GSI_MetGuess_Bundle(it),gsibundle_fv3lam_phyvar_nouv) if (laeroana_fv3cmaq) then call gsi_copy_bundle(GSI_ChemGuess_Bundle(it),gsibundle_fv3lam_tracerchem_nouv) end if @@ -3581,10 +3601,11 @@ subroutine wrfv3_netcdf(fv3filenamegin) add_saved,fv3filenamegin%dynvars,fv3filenamegin) call gsi_fv3ncdf_write(grd_fv3lam_tracer_ionouv,gsibundle_fv3lam_tracer_nouv, & add_saved,fv3filenamegin%tracers,fv3filenamegin) - if( if_model_dbz ) then + if( if_model_dbz .or. if_model_fed ) then call gsi_fv3ncdf_write(grd_fv3lam_phyvar_ionouv,gsibundle_fv3lam_phyvar_nouv,& add_saved,fv3filenamegin%phyvars,fv3filenamegin) end if + call gsi_fv3ncdf_writeuv(grd_fv3lam_uv,ges_u,ges_v,add_saved,fv3filenamegin) if (laeroana_fv3cmaq) then call gsi_fv3ncdf_write(grd_fv3lam_tracerchem_ionouv,gsibundle_fv3lam_tracerchem_nouv, & @@ -4343,7 +4364,7 @@ subroutine gsi_fv3ncdf_write(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3file work_a=hwork(1,:,:,ilevtot) - if( trim(varname) == 'ref_f3d' )then + if( trim(varname) == 'ref_f3d' .or. trim(adjustl(varname)) == 'flash_extent_density' )then iret=nf90_inquire_dimension(gfile_loc,1,name,len) if(trim(name)=='xaxis_1') nx_phy=len if( nx_phy == nxcase )then @@ -4386,7 +4407,7 @@ subroutine gsi_fv3ncdf_write(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3file deallocate(work_b_layout) enddo else - if( trim(varname) == 'ref_f3d' )then + if( trim(varname) == 'ref_f3d' .or. trim(varname) == 'flash_extent_density' )then work_b = 0.0_r_kind call check( nf90_get_var(gfile_loc,VarId,work_b_tmp,start = startloc_tmp, count = countloc_tmp) ) where(work_b_tmp < 0.0_r_kind) @@ -4419,7 +4440,7 @@ subroutine gsi_fv3ncdf_write(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3file deallocate(work_b_layout) enddo else - if( trim(varname) == 'ref_f3d' )then + if( trim(varname) == 'ref_f3d' .or. trim(varname) == 'flash_extent_density' )then if(phy_smaller_domain)then work_b_tmp = work_b(4:nxcase-3,4:nycase-3) else diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index d0ca1c0fbf..5a06eff27a 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -23,12 +23,13 @@ module gsimod use gsi_dbzOper, only: diag_radardbz use gsi_fedOper, only: diag_fed - use obsmod, only: doradaroneob,oneoblat,oneoblon,oneobheight,oneobvalue,oneobddiff,oneobradid,& + use obsmod, only: doradaroneob,dofedoneob,oneoblat,oneoblon,oneobheight,oneobvalue,oneobddiff,oneobradid,& radar_no_thinning,ens_hx_dbz_cut,static_gsi_nopcp_dbz,rmesh_dbz,& - rmesh_vr,zmesh_dbz,zmesh_vr,if_vterminal, if_model_dbz,if_vrobs_raw,if_use_w_vr,& + rmesh_vr,zmesh_dbz,zmesh_vr,if_vterminal, if_model_dbz,if_model_fed,innov_use_model_fed,if_vrobs_raw,if_use_w_vr,& minobrangedbz,maxobrangedbz,maxobrangevr,maxtiltvr,missing_to_nopcp,& ntilt_radarfiles,whichradar,& - minobrangevr,maxtiltdbz,mintiltvr,mintiltdbz,l2rwthin,hurricane_radar + minobrangevr,maxtiltdbz,mintiltvr,mintiltdbz,l2rwthin,hurricane_radar,& + r_hgt_fed use obsmod, only: lwrite_predterms, & lwrite_peakwt,use_limit,lrun_subdirs,l_foreaft_thin,lobsdiag_forenkf,& @@ -202,7 +203,7 @@ module gsimod use gsi_nstcouplermod, only: gsi_nstcoupler_init_nml use gsi_nstcouplermod, only: nst_gsi,nstinfo,zsea1,zsea2,fac_dtl,fac_tsl use ncepnems_io, only: init_nems,imp_physics,lupp - use wrf_vars_mod, only: init_wrf_vars + use wrf_vars_mod, only: init_wrf_vars,fed_exist,dbz_exist use gsi_rfv3io_mod,only : fv3sar_bg_opt use radarz_cst, only: mphyopt, MFflg use radarz_iface, only: init_mphyopt @@ -510,6 +511,15 @@ module gsimod ! 2023-07-30 Zhao - added namelist options for analysis of significant wave height ! (aka howv in GSI code): corp_howv, hwllp_howv ! (in namelist session rapidrefresh_cldsurf) +! +! 2023-09-14 H. Wang - add namelist option for FED EnVar DA. +! - if_model_fed=.true. : FED in background and ens. If +! perform FED DA, this has to be true along with fed in +! control/analysis and metguess vectors. If only run GSI observer, +! it can be false. +! - innov_use_model_fed=.true. : Use FED from BG to calculate innovation. +! this requires if_model_fed=.true. +! it works either an EnVar DA run or a GSI observer run. ! !EOP !------------------------------------------------------------------------- @@ -770,16 +780,17 @@ module gsimod use_sp_eqspace,lnested_loops,lsingleradob,thin4d,use_readin_anl_sfcmask,& luse_obsdiag,id_drifter,id_ship,verbose,print_obs_para,lsingleradar,singleradar,lnobalance, & missing_to_nopcp,minobrangedbz,minobrangedbz,maxobrangedbz,& - maxobrangevr,maxtiltvr,whichradar,doradaroneob,oneoblat,& + maxobrangevr,maxtiltvr,whichradar,doradaroneob,dofedoneob,oneoblat,& oneoblon,oneobheight,oneobvalue,oneobddiff,oneobradid,& rmesh_vr,zmesh_dbz,zmesh_vr, ntilt_radarfiles, whichradar,& radar_no_thinning,ens_hx_dbz_cut,static_gsi_nopcp_dbz,rmesh_dbz,& minobrangevr, maxtiltdbz, mintiltvr,mintiltdbz,if_vterminal,if_vrobs_raw,if_use_w_vr,& - if_model_dbz,imp_physics,lupp,netcdf_diag,binary_diag,l_wcp_cwm,aircraft_recon,diag_version,& + if_model_dbz,if_model_fed,innov_use_model_fed,imp_physics,lupp,netcdf_diag,binary_diag,l_wcp_cwm,aircraft_recon,diag_version,& write_fv3_incr,incvars_to_zero,incvars_zero_strat,incvars_efold,diag_version,& cao_check,lcalc_gfdl_cfrac,tau_fcst,efsoi_order,lupdqc,lqcoef,cnvw_option,l2rwthin,hurricane_radar,& l_reg_update_hydro_delz, l_obsprvdiag,& - l_use_dbz_directDA, l_use_rw_columntilt, ta2tb, optconv + l_use_dbz_directDA, l_use_rw_columntilt, ta2tb, optconv, & + r_hgt_fed ! GRIDOPTS (grid setup variables,including regional specific variables): ! jcap - spectral resolution @@ -1978,6 +1989,20 @@ subroutine gsimain_initialize endif endif + if (innov_use_model_fed .and. .not.if_model_fed) then + if(mype==0) write(6,*)' GSIMOD: invalid innov_use_model_fed=.true. but if_model_fed=.false.' + call die(myname_,'invalid innov_use_model_fed,if_model_fed, check namelist settings',330) + end if + + if (.not. (miter == 0 .or. lobserver) .and. if_model_fed .and. .not. fed_exist) then + if(mype==0) write(6,*)' GSIMOD: .not. (miter == 0 .or. lobserver) and if_model_fed=.true. but fed is not in anavinfo file' + call die(myname_,'Please check namelist parameters and/or add fed in anavinfo (contro/state_vector and met_guess) when miter > 0 and if_model_fed=.true.',332) + end if + + if (.not. (miter == 0 .or. lobserver) .and. if_model_dbz .and. .not. dbz_exist) then + if(mype==0) write(6,*)' GSIMOD: .not. (miter == 0 .or. lobserver) and if_model_dbz=.true. but dbz is not in anavinfo file' + call die(myname_,'Please check namelist parameters and/or add dbz in anavinfo (contro/state_vector and met_guess) when miter > 0 and if_model_fed=.true.',334) + end if ! Ensure valid number of horizontal scales if (nhscrf<0 .or. nhscrf>3) then diff --git a/src/gsi/intfed.f90 b/src/gsi/intfed.f90 new file mode 100644 index 0000000000..8cb16eba10 --- /dev/null +++ b/src/gsi/intfed.f90 @@ -0,0 +1,187 @@ +module intfedmod +!$$$ module documentation block +! . . . . +! module: intfedmod module for intfed and its tangent linear intfed_tl +! prgmmr: +! +! abstract: module for intfed and its tangent linear intfed_tl +! +! program history log: +! 2023-08-24 H. Wang - add tangent linear of fed operator to directly assimilate FED +! +! subroutines included: +! sub intfed_ +! +! variable definitions: +! +! attributes: +! language: f90 +! machine: +! +!$$$ end documentation block + +use m_obsNode, only: obsNode +use m_fedNode, only: fedNode +use m_fedNode, only: fedNode_typecast +use m_fedNode, only: fedNode_nextcast +use m_obsdiagNode, only: obsdiagNode_set +implicit none + +PRIVATE +PUBLIC intfed + +interface intfed; module procedure & + intfed_ +end interface + +contains + +subroutine intfed_(fedhead,rval,sval) +!$$$ subprogram documentation block +! . . . . +! subprogram: intfed apply nonlin qc operator for GLM FED +! +! abstract: apply observation operator for radar winds +! with nonlinear qc operator +! +! program history log: +! 2023-08-24 H.Wang - modified based on intdbz.f90 +! - using tangent linear fed operator + +! +! input argument list: +! fedhead - obs type pointer to obs structure +! sfed - current fed solution increment +! +! output argument list: +! rfed - fed results from fed observation operator +! +! attributes: +! language: f90 +! machine: ibm RS/6000 SP +! +!$$$ + use kinds, only: r_kind,i_kind + use constants, only: half,one,tiny_r_kind,cg_term,r3600 + use obsmod, only: lsaveobsens,l_do_adjoint,luse_obsdiag + use qcmod, only: nlnqc_iter,varqc_iter + use jfunc, only: jiter + use gsi_bundlemod, only: gsi_bundle + use gsi_bundlemod, only: gsi_bundlegetpointer + use gsi_4dvar, only: ladtest_obs + use wrf_vars_mod, only : fed_exist + implicit none + +! Declare passed variables + class(obsNode), pointer, intent(in ) :: fedhead + type(gsi_bundle), intent(in ) :: sval + type(gsi_bundle), intent(inout) :: rval + +! Declare local variables + integer(i_kind) j1,j2,j3,j4,j5,j6,j7,j8,ier,istatus +! real(r_kind) penalty + real(r_kind) val,w1,w2,w3,w4,w5,w6,w7,w8,valfed + real(r_kind) cg_fed,p0,grad,wnotgross,wgross,pg_fed + real(r_kind),pointer,dimension(:) :: sfed + real(r_kind),pointer,dimension(:) :: rfed + type(fedNode), pointer :: fedptr + +! If no fed obs type data return + if(.not. associated(fedhead))return + +! Retrieve pointers +! Simply return if any pointer not found + ier=0 + if(fed_exist)then + call gsi_bundlegetpointer(sval,'fed',sfed,istatus);ier=istatus+ier + call gsi_bundlegetpointer(rval,'fed',rfed,istatus);ier=istatus+ier + else + return + end if + + if(ier/=0)return + + + fedptr => fedNode_typecast(fedhead) + do while (associated(fedptr)) + j1=fedptr%ij(1) + j2=fedptr%ij(2) + j3=fedptr%ij(3) + j4=fedptr%ij(4) + j5=fedptr%ij(5) + j6=fedptr%ij(6) + j7=fedptr%ij(7) + j8=fedptr%ij(8) + w1=fedptr%wij(1) + w2=fedptr%wij(2) + w3=fedptr%wij(3) + w4=fedptr%wij(4) + w5=fedptr%wij(5) + w6=fedptr%wij(6) + w7=fedptr%wij(7) + w8=fedptr%wij(8) + + +! Forward model + if( fed_exist )then + val = w1* sfed(j1)+w2* sfed(j2)+w3* sfed(j3)+w4* sfed(j4)+ & + w5* sfed(j5)+w6* sfed(j6)+w7* sfed(j7)+w8* sfed(j8) + end if + + if(luse_obsdiag)then + if (lsaveobsens) then + grad = val*fedptr%raterr2*fedptr%err2 + !-- fedptr%diags%obssen(jiter) = grad + call obsdiagNode_set(fedptr%diags,jiter=jiter,obssen=grad) + + else + !-- if (fedptr%luse) fedptr%diags%tldepart(jiter)=val + if (fedptr%luse) call obsdiagNode_set(fedptr%diags,jiter=jiter,tldepart=val) + endif + endif + + if (l_do_adjoint) then + if (.not. lsaveobsens) then + if( .not. ladtest_obs ) val=val-fedptr%res + +! gradient of nonlinear operator + if (nlnqc_iter .and. fedptr%pg > tiny_r_kind .and. & + fedptr%b > tiny_r_kind) then + pg_fed=fedptr%pg*varqc_iter + cg_fed=cg_term/fedptr%b + wnotgross= one-pg_fed + wgross = pg_fed*cg_fed/wnotgross + p0 = wgross/(wgross+exp(-half*fedptr%err2*val**2)) + val = val*(one-p0) + endif + + if( ladtest_obs) then + grad = val + else + grad = val*fedptr%raterr2*fedptr%err2 + end if + + endif + +! Adjoint + if(fed_exist)then + valfed = grad + rfed(j1)=rfed(j1)+w1*valfed + rfed(j2)=rfed(j2)+w2*valfed + rfed(j3)=rfed(j3)+w3*valfed + rfed(j4)=rfed(j4)+w4*valfed + rfed(j5)=rfed(j5)+w5*valfed + rfed(j6)=rfed(j6)+w6*valfed + rfed(j7)=rfed(j7)+w7*valfed + rfed(j8)=rfed(j8)+w8*valfed + end if + + endif + + !fedptr => fedptr%llpoint + fedptr => fedNode_nextcast(fedptr) + end do + return +end subroutine intfed_ + +end module intfedmod diff --git a/src/gsi/m_berror_stats_reg.f90 b/src/gsi/m_berror_stats_reg.f90 index bf9fb20674..8730e56c3b 100644 --- a/src/gsi/m_berror_stats_reg.f90 +++ b/src/gsi/m_berror_stats_reg.f90 @@ -400,7 +400,7 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt integer(i_kind) :: nrf2_td2m,nrf2_mxtm,nrf2_mitm,nrf2_pmsl,nrf2_howv,nrf2_tcamt,nrf2_lcbas,nrf2_cldch integer(i_kind) :: nrf2_uwnd10m,nrf2_vwnd10m integer(i_kind) :: nrf3_sfwter,nrf3_vpwter - integer(i_kind) :: nrf3_dbz + integer(i_kind) :: nrf3_dbz,nrf3_fed integer(i_kind) :: nrf3_ql,nrf3_qi,nrf3_qr,nrf3_qs,nrf3_qg,nrf3_qnr,nrf3_w integer(i_kind) :: inerr,istat integer(i_kind) :: nsigstat,nlatstat,isig @@ -624,6 +624,7 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt nrf3_sf =getindex(cvars3d,'sf') nrf3_vp =getindex(cvars3d,'vp') nrf3_dbz=getindex(cvars3d,'dbz') + nrf3_fed=getindex(cvars3d,'fed') nrf2_sst=getindex(cvars2d,'sst') nrf2_gust=getindex(cvars2d,'gust') nrf2_vis=getindex(cvars2d,'vis') @@ -671,6 +672,16 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt vz(:,:,nrf3_dbz)=vz(:,:,nrf3_t) endif + if( nrf3_fed>0 )then + if(.not. nrf3_t>0) then + write(6,*)'not as expect,stop' + stop + endif + corz(:,:,nrf3_fed)=10.0_r_kind + hwll(:,:,nrf3_fed)=hwll(:,:,nrf3_t) + vz(:,:,nrf3_fed)=vz(:,:,nrf3_t) + endif + if (nrf3_oz>0) then factoz = 0.0002_r_kind*r25 corz(:,:,nrf3_oz)=factoz diff --git a/src/gsi/obsmod.F90 b/src/gsi/obsmod.F90 index 633bde91ab..c43a23c1e6 100644 --- a/src/gsi/obsmod.F90 +++ b/src/gsi/obsmod.F90 @@ -161,6 +161,7 @@ module obsmod ! observation provider and sub-provider information into ! obsdiags files (used for AutoObsQC) ! 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 ! ! Subroutines Included: ! sub init_obsmod_dflts - initialize obs related variables to default values @@ -188,6 +189,12 @@ module obsmod ! def diag_radardbz- namelist logical to compute/write (=true) radar ! reflectiivty diag files ! def diag_fed - namelist logical to compute/write (=true) flash extent density diag files +! def innov_use_model_fed - namelist logical. True: use (the FEB in background to calculate innovation +! False: calculate innvation use +! the obs operator in GSI +! def if_model_fed - namelist logical. True: Read in FED from background +! including from ensemble. +! def r_hgt_fed - height of fed observations ! def reduce_diag - namelist logical to produce reduced radiance diagnostic files ! def use_limit - parameter set equal to -1 if diag files produced or 0 if not diag files or reduce_diag ! def obs_setup - prefix for files passing pe relative obs data to setup routines @@ -487,7 +494,12 @@ module obsmod public :: iout_dbz, mype_dbz ! --- DBZ DA --- + ! ==== FED DA === + public :: if_model_fed, innov_use_model_fed + public :: r_hgt_fed public :: iout_fed, mype_fed + public :: dofedoneob + ! --- FED DA --- public :: obsmod_init_instr_table public :: obsmod_final_instr_table @@ -577,7 +589,7 @@ module obsmod real(r_kind) perturb_fact,time_window_max,time_offset,time_window_rad real(r_kind),dimension(50):: dmesh - + real(r_kind) r_hgt_fed integer(i_kind) nchan_total,ianldate integer(i_kind) ndat,ndat_types,ndat_times,nprof_gps integer(i_kind) lunobs_obs,nloz_v6,nloz_v8,nobskeep,nloz_omi @@ -621,8 +633,8 @@ module obsmod integer(i_kind) ntilt_radarfiles,tcp_posmatch,tcp_box logical :: ta2tb - logical :: doradaroneob - logical :: vr_dealisingopt, if_vterminal, if_model_dbz, inflate_obserr, if_vrobs_raw, if_use_w_vr, l2rwthin + logical :: doradaroneob,dofedoneob + logical :: vr_dealisingopt, if_vterminal, if_model_dbz,if_model_fed, innov_use_model_fed,inflate_obserr, if_vrobs_raw, if_use_w_vr, l2rwthin character(4) :: whichradar,oneobradid real(r_kind) :: oneoblat,oneoblon,oneobddiff,oneobvalue,oneobheight logical :: radar_no_thinning @@ -755,11 +767,15 @@ subroutine init_obsmod_dflts if_vrobs_raw=.false. if_use_w_vr=.true. if_model_dbz=.false. + if_model_fed=.false. + innov_use_model_fed=.false. inflate_obserr=.false. whichradar="KKKK" oneobradid="KKKK" doradaroneob=.false. + r_hgt_fed=6500_r_kind + dofedoneob=.false. oneoblat=-999_r_kind oneoblon=-999_r_kind oneobddiff=-999_r_kind diff --git a/src/gsi/read_dbz_netcdf.f90 b/src/gsi/read_dbz_netcdf.f90 index ecfb9169c4..845660168a 100644 --- a/src/gsi/read_dbz_netcdf.f90 +++ b/src/gsi/read_dbz_netcdf.f90 @@ -526,7 +526,8 @@ subroutine read_dbz_mrms_netcdf(nread,ndata,nodata,infile,obstype,lunout,sis,nob if(thislon>=r360) thislon=thislon-r360 if(thislon=r360 .or. thislat >90.0_r_kind) cycle + !-Convert back to radians thislat = thislat*deg2rad diff --git a/src/gsi/read_fed.f90 b/src/gsi/read_fed.f90 index c478b3d93f..3d3d098b08 100644 --- a/src/gsi/read_fed.f90 +++ b/src/gsi/read_fed.f90 @@ -9,6 +9,12 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs) ! 2019-09-20 Yaping Wang (CIMMS/OU) ! 2021-07-01 David Dowell (DCD; NOAA GSL) - added maximum flashes/min for observed FED ! +! 2023-10-18 Hongli Wang (NOAA GSL) +! - cleanup code, removed hardcoded obs height (6500m) +! - use height fron obs file if they are avaiable, otherwise +! use default value or value from namelist variable r_hgt_fed +! - return if NetCDF file open status /= nf90_noerror +! ! input argument list: ! infile - unit from which to read observation information file ! obstype - observation type to process @@ -29,13 +35,13 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs) !_____________________________________________________________________ ! use kinds, only: r_kind,r_double,i_kind - use constants, only: zero,one,deg2rad + use constants, only: zero,one,deg2rad,r60inv use convinfo, only: nconvtype,ctwind,icuse,ioctype - use gsi_4dvar, only: l4dvar,l4densvar,winlen + use gsi_4dvar, only: iwinbgn use gridmod, only: tll2xy use mod_wrfmass_to_a, only: wrfmass_obs_to_a8 use mpimod, only: npe - use obsmod, only: perturb_obs,iadatemn + use obsmod, only: perturb_obs,iadatemn,dofedoneob,oneoblat,oneoblon,r_hgt_fed use netcdf implicit none @@ -72,7 +78,6 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs) real(r_kind),allocatable,dimension(:,:):: cdata_out real(r_kind) :: federr, thiserr real(r_kind) :: hgt_fed(1) - data hgt_fed / 6500.0 / real(r_kind) :: i_maxloc,j_maxloc,k_maxloc integer(i_kind) :: kint_maxloc @@ -80,70 +85,55 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs) integer(i_kind) :: ndata2 integer(i_kind) :: ppp -! -! for read in bufr -! - real(r_kind) :: hdr(5),obs(1,3) - character(80):: hdrstr='SID XOB YOB DHR TYP' - character(80):: obsstr='FED' - - character(8) subset - character(8) station_id - real(r_double) :: rstation_id - equivalence(rstation_id,station_id) - integer(i_kind) :: lunin,idate - integer(i_kind) :: ireadmg,ireadsb + character(8) station_id + real(r_double) :: rstation_id + equivalence(rstation_id,station_id) - integer(i_kind) :: maxlvl - integer(i_kind) :: numlvl,numfed,nmsgmax,maxobs - integer(i_kind) :: k,iret - integer(i_kind) :: nmsg,ntb + integer(i_kind) :: maxlvl + integer(i_kind) :: numlvl,numfed,nmsgmax,maxobs + integer(i_kind) :: k,iret - real(r_kind),allocatable,dimension(:,:) :: fed3d_column ! 3D fed in column - real(r_kind),allocatable,dimension(:) :: utime ! time + real(r_kind),allocatable,dimension(:,:) :: fed3d_column ! 3D fed in column + real(r_kind),allocatable,dimension(:) :: fed3d_hgt ! fed height + real(r_kind),allocatable,dimension(:) :: utime ! time - integer(i_kind) :: ikx - real(r_kind) :: timeo,t4dv + integer(i_kind) :: ikx - character*128 :: myname='read_fed' + character*128 :: myname='read_fed' - real(r_kind) :: dlat, dlon ! rotated corrdinate - real(r_kind) :: dlat_earth, dlon_earth ! in unit of degree - real(r_kind) :: rlat00, rlon00 ! in unit of rad + real(r_kind) :: dlat, dlon ! rotated corrdinate + real(r_kind) :: dlat_earth, dlon_earth ! in unit of degree + real(r_kind) :: rlat00, rlon00 ! in unit of rad - logical :: l_psot_fed - logical :: l_latlon_fedobs - logical :: outside - integer :: unit_table + logical :: l_psot_fed + logical :: l_latlon_fedobs + logical :: outside ! for read netcdf - integer(i_kind) :: sec70,mins_an - integer(i_kind) :: varID, ncdfID, status - real(r_kind) :: timeb,twindm,rmins_an,rmins_ob - + integer(i_kind) :: sec70,mins_an + integer(i_kind) :: varID, ncdfID, status + real(r_kind) :: timeb,twindm,rmins_an,rmins_ob - unit_table = 23 -!********************************************************************** -! -! END OF DECLARATIONS....start of program -! - write(6,*) "r_kind=",r_kind - l_psot_fed = .FALSE. - l_latlon_fedobs = .TRUE. - - fedob = obstype == 'fed' - if(fedob) then - nreal=25 - else - write(6,*) ' illegal obs type in read_fed : obstype=',obstype - call stop2(94) - end if - if(perturb_obs .and. fedob)nreal=nreal+1 - write(6,*)'read_fed: nreal=',nreal - fedobs = .false. - ikx=0 - do i=1,nconvtype + hgt_fed = r_hgt_fed + + write(6,*) "r_kind=",r_kind + l_psot_fed = .FALSE. + l_latlon_fedobs = .TRUE. + + fedob = obstype == 'fed' + if (fedob) then + nreal=25 + else + write(6,*) ' illegal obs type in read_fed : obstype=',obstype + call stop2(94) + end if + if(perturb_obs .and. fedob)nreal=nreal+1 + write(6,*)'read_fed: nreal=',nreal + + fedobs = .false. + ikx=0 + do i=1,nconvtype if(trim(obstype) == trim(ioctype(i)) .and. abs(icuse(i))== 1) then fedobs=.true. ikx=i @@ -156,150 +146,24 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs) write(6,*) 'read_fed: abort read_fed !' return endif - end do - write(6,'(1x,A,A30,I4,A15,F7.3,A7)') & + end do + write(6,'(1x,A,A30,I4,A15,F7.3,A7)') & trim(myname),': fed in convinfo-->ikx=',ikx,' fed ob err:',thiserr," (fed)" - nread=0 - ndata=0 - nchanl=0 - ifn = 15 + nread=0 + ndata=0 + nchanl=0 + ifn = 15 - if(fedobs) then - maxlvl= 1 ! fed only has one level + if (fedobs) then + maxlvl= 1 ! fed only has one level - if(trim(infile) .eq. "fedbufr") then ! prebufr or netcdf format - !! get message and subset counts - ! nmsgmax and maxobs are read in from BUFR data file, not pre-set. - call getcount_bufr(infile,nmsgmax,maxobs) - write(6,*)'read_fed: nmsgmax=',nmsgmax,' maxobs=',maxobs - -! read in fed obs in bufr code format - lunin = 10 - allocate(fed3d_column(maxlvl+2+2,maxobs)) - - open ( unit = lunin, file = trim(infile),form='unformatted',err=200) - call openbf ( lunin, 'IN', lunin ) - open(unit_table,file='prepobs_kr.bufrtable') !temporily dump the bufr table, which is already saved in file - call dxdump(lunin,unit_table) - call datelen ( 10 ) - - nmsg=0 - ntb = 0 - - ndata =0 - ppp = 0 - msg_report: do while (ireadmg(lunin,subset,idate) == 0) - nmsg=nmsg+1 - if (nmsg>nmsgmax) then - write(6,*)'read_fed: messages exceed maximum ',nmsgmax - call stop2(50) - endif - loop_report: do while (ireadsb(lunin) == 0) - ntb = ntb+1 - if (ntb>maxobs) then - write(6,*)'read_fed: reports exceed maximum ',maxobs - call stop2(50) - endif - - ! Extract type, date, and location information from BUFR file - call ufbint(lunin,hdr,5,1,iret,hdrstr) - if(hdr(3) .gt. r90 ) write(6,*) "Inside read_fed.f90, hdr(2)=",hdr(2),"hdr(3)=",hdr(3) - if ( l_latlon_fedobs ) then - if(abs(hdr(3))>r90 .or. abs(hdr(2))>r360) cycle loop_report - if(hdr(2)== r360)hdr(2)=hdr(2)-r360 - if(hdr(2) < zero)hdr(2)=hdr(2)+r360 - end if - -! check time window in subset - if (l4dvar.or.l4densvar) then - t4dv=hdr(4) - if (t4dvwinlen) then - write(6,*)'read_fed: time outside window ',& - t4dv,' skip this report' - cycle loop_report - endif - else - timeo=hdr(4) - if (abs(timeo)>ctwind(ikx) .or. abs(timeo) > twind) then - write(6,*)'read_fed: time outside window ',& - timeo,' skip this report' - cycle loop_report - endif - endif -! read in observations - call ufbint(lunin,obs,1,3,iret,obsstr) !Single level bufr data, Rong Kong - if(obs(1,1) .gt. 5 ) write(6,*) "Inside read_fed.f90, obs(1,1)=",obs(1,1) - numlvl=min(iret,maxlvl) - if (numlvl .ne. maxlvl) then - write(6,*)' read_fed: numlvl is not equalt to maxlvl:',numlvl,maxlvl - end if - if(hdr(3) .gt. 90) write(6,*) "hdr(3)=",hdr(3) - if ( l_latlon_fedobs ) then - if(hdr(2)>= r360)hdr(2)=hdr(2)-r360 - if(hdr(2) < zero)hdr(2)=hdr(2)+r360 - fed3d_column(1,ntb)=hdr(2) ! observation location, earth lon - fed3d_column(2,ntb)=hdr(3) ! observation location, earth lat -! write(6,*) "Inside read_fed.f90, fed3d_column(1,ntb)=",fed3d_column(1,ntb),"fed3d_column(2,ntb)=",fed3d_column(2,ntb) - else - fed3d_column(1,ntb)=hdr(2)*10.0_r_kind ! observation location, grid index i - fed3d_column(2,ntb)=hdr(3)*10.0_r_kind ! observation location, grid index j - end if - - if (l_psot_fed .and. .NOT. l_latlon_fedobs ) then - do k=1,numlvl - if (NINT(fed3d_column(1,ntb)) .eq. 175 .and. NINT(fed3d_column(2,ntb)) .eq. 105 .and. & - NINT(hgt_fed(k)) .ge. 100 ) then - write(6,*) 'read_fed: single point/column obs run on grid: 175 105' - write(6,*) 'read_fed: found the pseudo single(column) fed obs:',fed3d_column(1:2,ntb),hgt_fed(k) - else - obs(1,1) = -999.0 - end if - end do - end if - - fed3d_column(3,ntb)=obs(1,1) - fed3d_column(4,ntb)=obs(1,2) - fed3d_column(5,ntb)=obs(1,3) - if (obs(1,1) == fed_lowbnd .or. obs(1,1) >= fed_lowbnd2 ) then - if (obs(1,1) == 0.0) then - ppp = ppp + 1 - endif - ndata = ndata + 1 - endif - - enddo loop_report - enddo msg_report - - write(6,*)'read_fed: messages/reports = ',nmsg,'/',ntb - print*,'number of Z that is less than 0 is ppp = ', ppp - numfed=ntb - -! - Finished reading fed observations from BUFR format data file -! - call closbf(lunin) - close(lunin) - - else ! NETCDF format !!!! Start reading fed observations from NETCDF format data file - ! CHECK IF DATA FILE EXISTS - + ! CHECK IF DATA FILE EXISTS ! OPEN NETCDF FILE status = nf90_open(TRIM(infile), NF90_NOWRITE, ncdfID) print*, '*** OPENING GOES FED OBS NETCDF FILE: ', infile, status - - - !------------------------ - ! Get date information - !------------------------- - ! status = nf90_get_att( ncdfID, nf90_global, 'year', idate5s(1) ) - ! print*, 'year ',status - ! status = nf90_get_att( ncdfID, nf90_global, 'month', idate5s(2) ) - ! status = nf90_get_att( ncdfID, nf90_global, 'day', idate5s(3) ) - ! status = nf90_get_att( ncdfID, nf90_global, 'hour', idate5s(4) ) - ! status = nf90_get_att( ncdfID, nf90_global, 'minute', idate5s(5) ) - ! read(idate5s(:) , *) idate5(:) - ! print*, idate5 + if(status/=nf90_noerr)return !------------------------ ! Get Dimension Info (1-D) @@ -307,72 +171,75 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs) status = nf90_inq_varid( ncdfID, 'numobs', varID ) status = nf90_get_var( ncdfID, varID, maxobs ) - !------------------------ - ! Allocate data arrays - !------------------------- - ALLOCATE( fed3d_column( 5, maxobs ) ) - ALLOCATE( utime( 1 ) ) ! seconds since from 2000-01-01 12:00 - - !------------------------ - ! Get useful data arrays - !------------------------- - ! LON - status = nf90_inq_varid( ncdfID, 'lon', varID ) - status = nf90_get_var( ncdfID, varID, fed3d_column(1, :) ) - ! LAT - status = nf90_inq_varid( ncdfID, 'lat', varID ) - status = nf90_get_var( ncdfID, varID, fed3d_column(2, :) ) - ! FED value - status = nf90_inq_varid( ncdfID, 'value', varID ) - status = nf90_get_var( ncdfID, varID, fed3d_column(3, :) ) - ! TIME - status = nf90_inq_varid( ncdfID, 'time', varID ) - status = nf90_get_var( ncdfID, varID, utime ) - - ! CLOSE NETCDF FILE - status = nf90_close( ncdfID ) - - - !-Obtain analysis time in minutes since reference date - sec70 = 694267200.0 ! seconds since from 1978-01-01 00:00 to 2000-01-01 12:00 - ! because the official GOES prescribed epoch time for GLM data is 2000-01-01 12:00:00 - - call w3fs21(iadatemn,mins_an) !mins_an -integer number of mins snce 01/01/1978 - rmins_an=mins_an !convert to real number - - ! SINCE ALL OBS WILL HAVE THE SAME TIME, CHECK TIME HERE: - rmins_ob = ( utime(1) + sec70 )/60 !Convert to Minutes from seconds - twindm = twind*60. !Convert to Minutes from hours - timeb = rmins_ob-rmins_an - - if(abs(timeb) > abs(twindm)) then + !------------------------ + ! Allocate data arrays + !------------------------- + ALLOCATE( fed3d_column( 5, maxobs ) ) + allocate( fed3d_hgt(maxobs) ) + ALLOCATE( utime( 1 ) ) ! seconds since from 2000-01-01 12:00 + fed3d_hgt = -999.0_r_kind + + !------------------------ + ! Get useful data arrays + !------------------------- + ! LON + status = nf90_inq_varid( ncdfID, 'lon', varID ) + status = nf90_get_var( ncdfID, varID, fed3d_column(1, :) ) + ! LAT + status = nf90_inq_varid( ncdfID, 'lat', varID ) + status = nf90_get_var( ncdfID, varID, fed3d_column(2, :) ) + ! FED value + status = nf90_inq_varid( ncdfID, 'value', varID ) + status = nf90_get_var( ncdfID, varID, fed3d_column(3, :) ) + ! TIME + status = nf90_inq_varid( ncdfID, 'time', varID ) + status = nf90_get_var( ncdfID, varID, utime ) + + ! FED height, optional variable + status = nf90_inq_varid( ncdfID, 'height', varID ) + if(status==nf90_noerr)& + status = nf90_get_var( ncdfID, varID, fed3d_hgt ) + + ! CLOSE NETCDF FILE + status = nf90_close( ncdfID ) + + + !-Obtain analysis time in minutes since reference date + sec70 = 694267200.0 ! seconds since from 1978-01-01 00:00 to 2000-01-01 12:00 + ! because the official GOES prescribed epoch time for GLM data is 2000-01-01 12:00:00 + + call w3fs21(iadatemn,mins_an) !mins_an -integer number of mins snce 01/01/1978 + rmins_an=mins_an !convert to real number + + ! SINCE ALL OBS WILL HAVE THE SAME TIME, CHECK TIME HERE: + rmins_ob = ( utime(1) + sec70 )/60 !Convert to Minutes from seconds + twindm = twind*60. !Convert to Minutes from hours + timeb = rmins_ob-rmins_an + + if(abs(timeb) > abs(twindm)) then print*, 'WARNING: ALL FED OBSERVATIONS OUTSIDE ASSIMILATION TIME WINDOW: ', timeb, twindm - ! goto 314 - endif - numfed = maxobs - do i=1,numfed + endif + + !time relative to the beginning of the da time window + timeb=real(rmins_ob-iwinbgn,r_kind) + + numfed = maxobs + do i=1,numfed if (fed3d_column( 3, i ) >= fed_lowbnd2 .or. fed3d_column( 3, i ) == fed_lowbnd ) then ndata = ndata + 1 end if - end do - end if ! end if prebufr or netcdf format + end do write(6,*)'read_fed: total no. of obs = ',ndata nread=ndata nodata=ndata !!! - Finished reading fed observations from NETCDF format data file - - allocate(cdata_out(nreal,ndata)) -! ! do i=1,numfed do k=1,maxlvl - -! DCD 1 July 2021 if (fed3d_column(k+2,i) .gt. fed_highbnd) fed3d_column(k+2,i) = fed_highbnd - end do end do @@ -382,38 +249,56 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs) write(6,*) k,maxval(fed3d_column(k+2,:)),minval(fed3d_column(k+2,:)) end do - i_maxloc=-1.0 j_maxloc=-1.0 k_maxloc=-1.0 kint_maxloc=-1 fed_max=-999.99 ndata2=0 + + ILOOP : & do i=1,numfed + if(fed3d_hgt(i) > 0.0_r_kind)then + hgt_fed=fed3d_hgt(i) + else + hgt_fed = r_hgt_fed + end if do k=1,maxlvl - if( fed3d_column(k+2,i) >= fed_lowbnd2 .or. fed3d_column(k+2,i) == fed_lowbnd) then !Rong Kong + if( fed3d_column(k+2,i) >= fed_lowbnd2 .or. fed3d_column(k+2,i) == fed_lowbnd .or. dofedoneob) then !Rong Kong dlon_earth = fed3d_column(1,i) ! longitude (degrees) of observation ! ilone=18 ! index of longitude (degrees) dlat_earth = fed3d_column(2,i) ! latitude (degrees) of observation ! ilate=19 ! index of latitude (degrees) - !-Check format of longitude and correct if necessary - if(dlon_earth>=r360) dlon_earth=dlon_earth-r360 - if(dlon_earth=r360 .or. dlat_earth >90.0_r_kind) cycle - !-Convert back to radians + if (dofedoneob) then + dlat_earth=oneoblat + dlon_earth=oneoblon + endif + + !-Check format of longitude and correct if necessary + if(dlon_earth>=r360) dlon_earth=dlon_earth-r360 + if(dlon_earth=r360 .or. dlat_earth >90.0_r_kind) cycle + + !-Convert back to radians rlon00 = dlon_earth*deg2rad rlat00 = dlat_earth*deg2rad call tll2xy(rlon00,rlat00,dlon,dlat,outside) + if (dofedoneob) then + if (outside) then + write(6,*)'READ_FED: ONE OB OUTSIDE; STOP2(61) ',dlat_earth,dlon_earth + call stop2(61) + end if + end if if (outside) cycle - !If observation is outside the domain - ! then cycle, but don't increase - ! range right away. - ! Domain could be rectangular, so ob - ! may be out of - ! range at one end, but not the - ! other. + !If observation is outside the domain + ! then cycle, but don't increase + ! range right away. + ! Domain could be rectangular, so ob + ! may be out of + ! range at one end, but not the + ! other. ndata2=ndata2+1 cdata_out( 1,ndata2) = thiserr ! obs error (flashes/min) - inflated/adjusted @@ -425,11 +310,11 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs) cdata_out( 4,ndata2) = hgt_fed(k) ! obs absolute height (m) above MSL ! ipres=4 ! index of pressure cdata_out( 5,ndata2) = fed3d_column(k+2,i) ! FED value - ! idbzob=5 ! index of dbz observation + cdata_out( 6,ndata2) = rstation_id ! station id (charstring equivalent to real double) ! id=6 ! index of station id - cdata_out( 7,ndata2) = 0.0_r_kind ! observation time in data array + cdata_out( 7,ndata2) = timeb*r60inv ! observation time in data array ! itime=7 ! index of observation time in data array cdata_out( 8,ndata2) = ikx ! ob type ! ikxx=8 ! index of ob type @@ -472,7 +357,9 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs) cdata_out(26,ndata2) = 1.0_r_kind ! obs perturbation ! iptrb=26 ! index of q perturbation end if -! print*,'cdata_out(:,ndata2)=',cdata_out(:,ndata2) + + if( dofedoneob ) exit ILOOP + if(fed3d_column(k+2,i) > fed_max)then kint_maxloc=k k_maxloc=real(k,r_kind) @@ -480,29 +367,27 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs) i_maxloc=fed3d_column(1,i) fed_max =fed3d_column(k+2,i) end if + endif - enddo - enddo + enddo ! k + enddo ILOOP ! i !---all looping done now print diagnostic output - write(6,*)'READ_FED: Reached eof on FED file' - write(6,*)'READ_FED: # read in obs. number =',nread - write(6,*)'READ_FED: # read in obs. number for further processing =',ndata2 - ! write(6,*)'READ_FED: dlon_earth', cdata_out(18,10:15) + write(6,*)'READ_FED: Reached eof on FED file' + write(6,*)'READ_FED: # read in obs. number =',nread + write(6,*)'READ_FED: # read in obs. number for further processing =',ndata2 +! write(6,*)'READ_FED: dlon_earth', cdata_out(18,10:15) ilon=2 ! array index for longitude ilat=3 ! array index for latitude in obs information array ndata=ndata2 nodata=ndata2 - !---Write observations to scratch file---! +!---Write observations to scratch file---! -! if(ndata > 0 ) then - call count_obs(ndata,nreal,ilat,ilon,cdata_out,nobs) - write(lunout) obstype,sis,nreal,nchanl,ilat,ilon,ndata - write(lunout) ((cdata_out(k,i),k=1,nreal),i=1,ndata) - ! print*,'cdata_out',cdata_out -! endif + call count_obs(ndata,nreal,ilat,ilon,cdata_out,nobs) + write(lunout) obstype,sis,nreal,nchanl,ilat,ilon,ndata + write(lunout) ((cdata_out(k,i),k=1,nreal),i=1,ndata) deallocate(cdata_out) if (allocated(fed3d_column)) deallocate(fed3d_column) @@ -511,15 +396,9 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs) 'read_fed: max fed =',fed_max, '@ i j k =', & i_maxloc,j_maxloc,k_maxloc,kint_maxloc - end if -! close(lunout) ! ???? - return - -200 continue - write(6,*) 'read_fed, Warning : cannot find or open bufr fed data file: ', trim(infile) + end if + return -314 continue -print* ,'FINISHED WITH READ_FED' end subroutine read_fed ! ! diff --git a/src/gsi/setupfed.f90 b/src/gsi/setupfed.f90 index cf6334e567..dbb2f56111 100644 --- a/src/gsi/setupfed.f90 +++ b/src/gsi/setupfed.f90 @@ -17,7 +17,9 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa ! - added a second option (tanh) for observation operator, based on the ! work of Sebok and Back (2021, unpublished) ! - capped maximum model FED -! +! Hongli Wang NOAA GSL 2023-09-14 +! - Add option to use fed from background file to calculate fed innov +! - cleanup code, removed hardcoded obs height (6500m) ! use mpeu_util, only: die,perr use kinds, only: r_kind,r_single,r_double,i_kind @@ -30,6 +32,7 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa use obsmod, only: rmiss_single,& lobsdiagsave,nobskeep,lobsdiag_allocated,time_offset use obsmod, only: oberror_tune + use obsmod, only: if_model_fed,innov_use_model_fed,dofedoneob,oneobddiff,oneobvalue use m_obsNode, only: obsNode use m_fedNode, only: fedNode use m_fedNode, only: fedNode_appendto @@ -93,7 +96,6 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa real(r_kind),parameter:: d_coeff = -1.215e9_r_kind ! d (kg) in tanh operator real(r_kind),parameter:: fed_highbnd = 8.0_r_kind ! DCD: Back (2023, unpublished) for FV3 - real(r_kind),parameter:: fed_height = 6500.0_r_kind ! assumed height (m) of FED observations real(r_kind),parameter:: r0_001 = 0.001_r_kind real(r_kind),parameter:: r8 = 8.0_r_kind real(r_kind),parameter:: ten = 10.0_r_kind @@ -135,7 +137,7 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa real(r_kind),allocatable,dimension(:,:,: ) :: ges_ps real(r_kind),allocatable,dimension(:,:,:,:) :: ges_q real(r_kind),allocatable,dimension(:,:,:,:) :: ges_qg,ges_qg_mask - + real(r_kind),allocatable,dimension(:,:,:,:) :: ges_fed real(r_kind) :: presq real(r_kind) :: T1D,RHO real(r_kind) :: glmcoeff = 2.088_r_kind*10.0**(-8.0) ! Allen et al. (2016,MWR) @@ -230,6 +232,10 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa muse(i)=nint(data(iuse,i)) <= jiter end do + if (dofedoneob) then + muse=.true. + end if + numequal=0 numnotequal=0 @@ -293,10 +299,10 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa print*, 'nsig = ', nsig print*, 'lon2 = ', lon2 print*, 'lat2 = ', lat2 - + if (.not. innov_use_model_fed .or. .not. if_model_fed) then ! compute graupel mass, in kg per 15 km x 15 km column - do jj=1,nfldsig - do k=1,nsig + do jj=1,nfldsig + do k=1,nsig do i=1,lon2 do j=1,lat2 !How to handle MPI???? do igx=1,2*ngx+1 !horizontal integration of qg around point to calculate FED @@ -312,12 +318,12 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa end do !jgy end do !j end do !i - end do !k - end do !jj + end do !k + end do !jj ! compute FED, in flashes/min - do jj=1,nfldsig - do i=1,lon2 + do jj=1,nfldsig + do i=1,lon2 do j=1,lat2 if (fed_obs_ob_shape .eq. 1) then rp(j,i,jj) = CM * glmcoeff * rp(j,i,jj) @@ -331,9 +337,8 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa end if if (rp(j,i,jj) .gt. fed_highbnd) rp(j,i,jj) = fed_highbnd end do !j - end do !i - end do !jj - + end do !i + end do !jj !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! write(6,*) 'fed_obs_ob_shape=',fed_obs_ob_shape if (fed_obs_ob_shape .eq. 2) then @@ -344,7 +349,7 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa end if write(6,*) 'fed_highbnd=',fed_highbnd write(6,*) 'maxval(ges_qg)=',maxval(ges_qg),'pe=',mype - + end if ! .not. innov_use_model_fed .or. .not. if_model_fed !============================================================================================ @@ -352,28 +357,21 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa nlon_ll=size(ges_qg,2) nsig_ll=size(ges_qg,3) nfld_ll=size(ges_qg,4) - -! - Observation times are checked in read routine - comment out for now - -! call dtime_setup() - -!print*,"maxval(data(ifedob,:)),mmaxval(data(ilat,:))=",minval(data(ifedob,:)),maxval(data(ifedob,:)),maxval(data(ilat,:)) -!write(6,*) "OKOKOKOKOK, nobs=", nobs do i=1,nobs - dtime=data(itime,i) - dlat=data(ilat,i) - dlon=data(ilon,i) - - dlon8km=data(iprvd,i) !iprvd=23 - dlat8km=data(isprvd,i) !isprvd=24 - - dpres=data(ipres,i) ! from rdararef_mosaic2: this height abv MSL - ikx = nint(data(ikxx,i)) - error=data(ier2,i) - slat=data(ilate,i)*deg2rad ! needed when converting geophgt to - dlon_earth = data(ilone,i) !the lontitude and latitude on the obs pts. - dlat_earth = data(ilate,i) + dtime=data(itime,i) + dlat=data(ilat,i) + dlon=data(ilon,i) + + dlon8km=data(iprvd,i) !iprvd=23 + dlat8km=data(isprvd,i) !isprvd=24 + + dpres=data(ipres,i) ! from rdararef_mosaic2: this height abv MSL + ikx = nint(data(ikxx,i)) + error=data(ier2,i) + slat=data(ilate,i)*deg2rad ! needed when converting geophgt to + dlon_earth = data(ilone,i) !the lontitude and latitude on the obs pts. + dlat_earth = data(ilate,i) ! geometric hgh (hges --> zges below) if (nobs_bins>1) then @@ -382,7 +380,7 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa ibin = 1 end if - IF (ibin<1.OR.ibin>nobs_bins) write(6,*)mype,'Error nobs_bins,ibin= ',nobs_bins,ibin + if (ibin<1.or.ibin>nobs_bins) write(6,*)mype,'Error nobs_bins,ibin= ',nobs_bins,ibin if (luse_obsdiag) my_diagLL => odiagLL(ibin) @@ -402,28 +400,28 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa end if ! Interpolate terrain height(model elevation) to obs location. - call tintrp2a11(ges_z,zsges,dlat,dlon,dtime,hrdifsig,& - mype,nfldsig) + call tintrp2a11(ges_z,zsges,dlat,dlon,dtime,hrdifsig,& + mype,nfldsig) ! print*,'i,after tintrp2all',i,mype,dlat,zsges ! 1. dpres (MRMS obs height is height above MSL) is adjusted by zsges, so it ! is changed to height relative to model elevation (terrain). ! because in GSI, geop_hgtl is the height relative to terrain (ges_z) ! (subroutine guess_grids) - dpres=dpres-zsges - if (dpres 0_r_kind) then + data(ifedob,i) = oneobvalue + ddiff = data(ifedob,i) - FEDMdiag(i) + else + ddiff = oneobddiff + data(ifedob,i) = FEDMdiag(i)+ddiff + oneobvalue = data(ifedob,i) + endif + write(6,*)"FED_ONEOB: O_Val,B_Val= ",data(ifedob,i),FEDMdiag(i) + write(6,*)"FED_ONEOB: Innov,Error= ",ddiff,magoberr + else + data(ifedob,i) = oneobvalue + ddiff = data(ifedob,i) - FEDMdiag(i) + end if + end if !oneob ! Gross error checks obserror = one/max(ratio_errors*error,tiny_r_kind) obserrlm = max(cermin(ikx),min(cermax(ikx),obserror)) - residual = abs(ddiff) != y-H(xb) ratio = residual/obserrlm != y-H(xb)/sqrt(R) @@ -782,7 +788,7 @@ end subroutine check_vars_ subroutine init_vars_ ! use radaremul_cst, only: mphyopt - + use obsmod, only: if_model_fed real(r_kind),dimension(:,: ),pointer:: rank2=>NULL() real(r_kind),dimension(:,:,:),pointer:: rank3=>NULL() character(len=5) :: varname @@ -822,6 +828,28 @@ subroutine init_vars_ call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank2,istatus) ges_z(:,:,ifld)=rank2 end do + + if(if_model_fed)then + ! get fed .... + varname='fed' + call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus) + if (istatus==0) then + if(allocated(ges_fed))then + write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' + call stop2(999) + endif + allocate(ges_fed(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig)) + ges_fed(:,:,:,1)=rank3 + do ifld=2,nfldsig + call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus) + ges_fed(:,:,:,ifld)=rank3 + enddo + else + write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus + call stop2(999) + endif + endif + else write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus call stop2(999) @@ -938,7 +966,7 @@ subroutine contents_binary_diag_(odiag) rdiagbuf(5,ii) = data(istnelv,i) ! station elevation (meters) rdiagbuf(6,ii) = presq ! observation pressure (hPa) - rdiagbuf(7,ii) = fed_height ! observation height (meters) + rdiagbuf(7,ii) = data(iobshgt,i) ! observation height (meters) rdiagbuf(8,ii) = (dtime*r60)-time_offset ! obs time (sec 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 @@ -1006,7 +1034,7 @@ subroutine contents_netcdf_diag_(odiag) 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(presq) ) - call nc_diag_metadata("Height", sngl(fed_height) ) + call nc_diag_metadata("Height", sngl(data(iobshgt,i)) ) call nc_diag_metadata("Time", sngl(dtime*r60-time_offset)) call nc_diag_metadata("Prep_QC_Mark", sngl(data(iqc,i)) ) call nc_diag_metadata("Prep_Use_Flag", sngl(data(iuse,i)) ) @@ -1048,6 +1076,7 @@ subroutine final_vars_ ! if(allocated(ges_tv)) deallocate(ges_tv) if(allocated(ges_ps)) deallocate(ges_ps) if(allocated(ges_qg)) deallocate(ges_qg) + if(allocated(ges_fed)) deallocate(ges_fed) end subroutine final_vars_ subroutine init_qcld(t_cld, qxmin_cld, icat_cld, t_dpnd) diff --git a/src/gsi/stpfed.f90 b/src/gsi/stpfed.f90 new file mode 100644 index 0000000000..2a69dd08ec --- /dev/null +++ b/src/gsi/stpfed.f90 @@ -0,0 +1,171 @@ +module stpfedmod + +!$$$ module documentation block +! . . . . +! module: stpfedmod module for stpfed and its tangent linear stpfed_tl +! prgmmr: +! +! abstract: module for stpfed and its tangent linear stpfed_tl +! +! program history log: +! 2023-08-23 H. Wang - Modified based on sftdbzmod +! - add adjoint of fed operator +! +! subroutines included: +! sub stpfed +! +! attributes: +! language: f90 +! machine: +! +!$$$ end documentation block + +implicit none + +PRIVATE +PUBLIC stpfed + +contains + +subroutine stpfed(fedhead,rval,sval,out,sges,nstep) +!$$$ subprogram documentation block +! . . . . +! subprogram: stpfed calculate penalty and contribution to +! stepsize with nonlinear qc added. +! prgmmr: derber org: np23 date: 1991-02-26 +! +! +! program history log: +! 2019-08-23 H.Wang - added for FED DA +! input argument list: +! fedhead +! sges - step size estimates (nstep) +! nstep - number of step sizes (== 0 means use outer iteration value) +! +! output argument list - output for step size calculation +! out(1:nstep) - penalty from fed 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: half,one,two,tiny_r_kind,cg_term,zero_quad,r3600 + use gsi_bundlemod, only: gsi_bundle + use gsi_bundlemod, only: gsi_bundlegetpointer + use gridmod, only: wrf_mass_regional, fv3_regional + use wrf_vars_mod, only : fed_exist + use m_obsNode, only: obsNode + use m_fedNode , only: fedNode + use m_fedNode , only: fedNode_typecast + use m_fedNode , only: fedNode_nextcast +! use directDA_radaruse_mod, only: l_use_fed_directDA + use radarz_cst, only: mphyopt + + implicit none + +! Declare passed variables + class(obsNode), pointer ,intent(in ) :: fedhead + integer(i_kind) ,intent(in ) :: nstep + real(r_quad),dimension(max(1,nstep)),intent(inout) :: out + type(gsi_bundle) ,intent(in ) :: rval,sval + real(r_kind),dimension(max(1,nstep)),intent(in ) :: sges + +! Declare local variables + integer(i_kind) ier,istatus + integer(i_kind) j1,j2,j3,j4,j5,j6,j7,j8,kk + real(r_kind) w1,w2,w3,w4,w5,w6,w7,w8 + real(r_kind) valfed + real(r_kind) fedcur + real(r_kind) cg_fed,fed,wgross,wnotgross + real(r_kind),dimension(max(1,nstep))::pen + real(r_kind) pg_fed + real(r_kind),pointer,dimension(:) :: sfed + real(r_kind),pointer,dimension(:) :: rfed + type(fedNode), pointer :: fedptr + + out=zero_quad + +! If no fed data return + if(.not. associated(fedhead))return + +! Retrieve pointers +! Simply return if any pointer not found + ier=0 + if(fed_exist)then + call gsi_bundlegetpointer(sval,'fed',sfed,istatus);ier=istatus+ier + call gsi_bundlegetpointer(rval,'fed',rfed,istatus);ier=istatus+ier + else + return + end if + + if(ier/=0)return + + fedptr => fedNode_typecast(fedhead) + do while (associated(fedptr)) + if(fedptr%luse)then + if(nstep > 0)then + j1=fedptr%ij(1) + j2=fedptr%ij(2) + j3=fedptr%ij(3) + j4=fedptr%ij(4) + j5=fedptr%ij(5) + j6=fedptr%ij(6) + j7=fedptr%ij(7) + j8=fedptr%ij(8) + w1=fedptr%wij(1) + w2=fedptr%wij(2) + w3=fedptr%wij(3) + w4=fedptr%wij(4) + w5=fedptr%wij(5) + w6=fedptr%wij(6) + w7=fedptr%wij(7) + w8=fedptr%wij(8) + + if( fed_exist )then + valfed= w1* rfed(j1)+w2*rfed(j2)+w3*rfed(j3)+w4*rfed(j4)+ & + w5* rfed(j5)+w6*rfed(j6)+w7*rfed(j7)+w8*rfed(j8) + + fedcur= w1* sfed(j1)+w2* sfed(j2)+w3* sfed(j3)+w4*sfed(j4)+ & + w5* sfed(j5)+w6* sfed(j6)+w7* sfed(j7)+w8* sfed(j8)- & + fedptr%res + end if + + + do kk=1,nstep + fed=fedcur+sges(kk)*valfed + pen(kk)=fed*fed*fedptr%err2 + end do + else + pen(1)=fedptr%res*fedptr%res*fedptr%err2 + end if + +! Modify penalty term if nonlinear QC + if (nlnqc_iter .and. fedptr%pg > tiny_r_kind .and. & + fedptr%b > tiny_r_kind) then + + pg_fed=fedptr%pg*varqc_iter + cg_fed=cg_term/fedptr%b + wnotgross= one-pg_fed + wgross = pg_fed*cg_fed/wnotgross + do kk=1,max(1,nstep) + pen(kk)= -two*log((exp(-half*pen(kk)) + wgross)/(one+wgross)) + end do + end if + + out(1) = out(1)+pen(1)*fedptr%raterr2 + kk=1 + do kk=2,nstep + out(kk) = out(kk)+(pen(kk)-pen(1))*fedptr%raterr2 + end do + end if + + fedptr => fedNode_nextcast(fedptr) + + end do + return +end subroutine stpfed + +end module stpfedmod diff --git a/src/gsi/wrf_vars_mod.f90 b/src/gsi/wrf_vars_mod.f90 index 97c36c43cf..f7a5e6c83d 100644 --- a/src/gsi/wrf_vars_mod.f90 +++ b/src/gsi/wrf_vars_mod.f90 @@ -39,6 +39,8 @@ module wrf_vars_mod use mpimod, only: mype use control_vectors, only: nc3d,cvars3d use kinds, only: i_kind +use gsi_metguess_mod, only: gsi_metguess_get +use constants, only: max_varname_length implicit none private ! public methods @@ -46,21 +48,54 @@ module wrf_vars_mod ! common block variables public :: w_exist public :: dbz_exist +public :: fed_exist -logical,save :: w_exist, dbz_exist +logical,save :: w_exist, dbz_exist, fed_exist contains subroutine init_wrf_vars -integer(i_kind) ii +integer(i_kind) ii,istatus +character(max_varname_length),allocatable,dimension(:) :: cloud +integer(i_kind) ncloud +logical :: dbz_cloud_exist,fed_cloud_exist w_exist=.false. dbz_exist=.false. +fed_exist=.false. +dbz_cloud_exist=.false. +fed_cloud_exist=.false. + do ii=1,nc3d if(mype == 0 ) write(6,*)"anacv cvars3d is ",cvars3d(ii) if(trim(cvars3d(ii)) == 'w'.or.trim(cvars3d(ii))=='W') w_exist=.true. if(trim(cvars3d(ii))=='dbz'.or.trim(cvars3d(ii))=='DBZ') dbz_exist=.true. + if(trim(cvars3d(ii))=='fed'.or.trim(cvars3d(ii))=='FED') fed_exist=.true. enddo +! Inquire about clouds + +call gsi_metguess_get('clouds::3d',ncloud,istatus) +if (ncloud>0) then + allocate(cloud(ncloud)) + call gsi_metguess_get('clouds::3d',cloud,istatus) +endif + +do ii=1,ncloud + if(mype == 0 ) write(6,*)"metguess cloud3d is ",cloud(ii) + if(trim(cloud(ii))=='fed'.or.trim(cloud(ii))=='FED')fed_cloud_exist=.true. + if(trim(cloud(ii))=='dbz'.or.trim(cloud(ii))=='DBZ')dbz_cloud_exist=.true. +end do + +if(.not.fed_exist .or. .not.fed_cloud_exist )then + fed_exist=.false. +endif + +if(.not.dbz_exist .or. .not.dbz_cloud_exist )then + dbz_exist=.false. +endif + +if(ncloud>0) deallocate(cloud) + end subroutine init_wrf_vars end module wrf_vars_mod