diff --git a/src/gsi/cplr_read_wrf_mass_guess.f90 b/src/gsi/cplr_read_wrf_mass_guess.f90 index 6afc077d6..54db338ac 100644 --- a/src/gsi/cplr_read_wrf_mass_guess.f90 +++ b/src/gsi/cplr_read_wrf_mass_guess.f90 @@ -1394,6 +1394,7 @@ subroutine read_wrf_mass_netcdf_guess_wrf(this,mype) aerotot_guess,init_aerotot_guess,wrf_pm2_5,aero_ratios use rapidrefresh_cldsurf_mod, only: l_hydrometeor_bkio,l_gsd_soiltq_nudge use rapidrefresh_cldsurf_mod, only: i_use_2mq4b,i_use_2mt4b + use rapidrefresh_cldsurf_mod, only: i_howv_3dda, i_gust_3dda use wrf_mass_guess_mod, only: soil_temp_cld,isli_cld,ges_xlon,ges_xlat,ges_tten,create_cld_grids use gsi_bundlemod, only: GSI_BundleGetPointer use gsi_metguess_mod, only: gsi_metguess_get,GSI_MetGuess_Bundle @@ -1446,6 +1447,7 @@ subroutine read_wrf_mass_netcdf_guess_wrf(this,mype) integer(i_kind) i_qc,i_qi,i_qr,i_qs,i_qg,i_qnr,i_qni,i_qnc,i_w,i_dbz integer(i_kind) kqc,kqi,kqr,kqs,kqg,kqnr,kqni,kqnc,i_xlon,i_xlat,i_tt,ktt integer(i_kind) i_th2,i_q2,i_soilt1,ksmois,ktslb + integer(i_kind) i_howv, i_gust integer(i_kind) ier, istatus integer(i_kind) n_actual_clouds integer(i_kind) iv,n_gocart_var @@ -1469,6 +1471,8 @@ subroutine read_wrf_mass_netcdf_guess_wrf(this,mype) real(r_kind), pointer :: ges_tv_it (:,:,:)=>NULL() real(r_kind), pointer :: ges_q_it (:,:,:)=>NULL() real(r_kind), pointer :: ges_w_it (:,:,:)=>NULL() + real(r_kind), pointer :: ges_howv_it (:,:)=>NULL() + real(r_kind), pointer :: ges_gust_it (:,:)=>NULL() real(r_kind), pointer :: ges_qc (:,:,:)=>NULL() real(r_kind), pointer :: ges_qi (:,:,:)=>NULL() @@ -1560,6 +1564,8 @@ subroutine read_wrf_mass_netcdf_guess_wrf(this,mype) if(l_gsd_soilTQ_nudge) num_mass_fields=num_mass_fields+2*(nsig_soil-1)+1 if(i_use_2mt4b > 0 ) num_mass_fields=num_mass_fields + 2 if(i_use_2mq4b > 0 .and. i_use_2mt4b <=0 ) num_mass_fields=num_mass_fields + 1 + if( i_howv_3dda > 0 ) num_mass_fields = num_mass_fields + 1 + if( i_gust_3dda > 0 ) num_mass_fields = num_mass_fields + 1 if (laeroana_gocart .and. wrf_pm2_5 ) then if(mype==0) write(6,*)'laeroana_gocart canoot be both true' @@ -1740,6 +1746,18 @@ subroutine read_wrf_mass_netcdf_guess_wrf(this,mype) write(identity(i),'("record ",i3,"--th2(",i2,")")')i,k jsig_skip(i)=0 ; igtype(i)=1 endif + ! for wave height (howv is after tsk/q2/soilt1/th2, and before cloud hydrometers) + if ( i_howv_3dda >0 ) then + i=i+1 ; i_howv=i ! howv + write(identity(i),'("record ",i3,"--howv")')i + jsig_skip(i)=0 ; igtype(i)=1 + end if + ! for wind gust (gust is after tsk/q2/soilt1/th2, and before cloud hydrometers) + if ( i_gust_3dda >0 ) then + i=i+1 ; i_gust=i ! gust + write(identity(i),'("record ",i3,"--gust")')i + jsig_skip(i)=0 ; igtype(i)=1 + end if ! for cloud array if(l_hydrometeor_bkio .and. n_actual_clouds>0) then i_qc=i+1 @@ -1944,6 +1962,14 @@ subroutine read_wrf_mass_netcdf_guess_wrf(this,mype) call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 't2m',ges_t2m_it, istatus );ier=ier+istatus if (ier/=0) call die(trim(myname),'cannot get pointers for t2m,ier =',ier) endif + if ( i_howv_3dda >0 ) then + call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'howv',ges_howv_it,istatus );ier=ier+istatus + if (ier/=0) call die(trim(myname),'cannot get pointers for met-field: howv, ier =',ier) + end if + if ( i_gust_3dda >0 ) then + call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'gust',ges_gust_it,istatus );ier=ier+istatus + if (ier/=0) call die(trim(myname),'cannot get pointers for met-field: gust, ier =',ier) + end if if (l_gsd_soilTQ_nudge) then call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'tskn',ges_tsk_it, istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'tsoil',ges_soilt1_it,istatus);ier=ier+istatus @@ -2273,6 +2299,14 @@ subroutine read_wrf_mass_netcdf_guess_wrf(this,mype) ges_q2_it(j,i)=real(all_loc(j,i,i_0+i_q2),r_kind) ges_q2_it(j,i)=ges_q2_it(j,i)/(one+ges_q2_it(j,i)) endif + ! wave height (howv) + if ( i_howv_3dda >0 ) then + ges_howv_it(j,i) = real(all_loc(j,i,i_0+i_howv),r_kind) + end if + ! wind gust (gust) + if ( i_gust_3dda >0 ) then + ges_gust_it(j,i) = real(all_loc(j,i,i_0+i_gust),r_kind) + end if ! for cloud analysis if(l_hydrometeor_bkio .and. n_actual_clouds>0) then soil_temp_cld(j,i,it)=soil_temp(j,i,it) diff --git a/src/gsi/cplr_regional_io.f90 b/src/gsi/cplr_regional_io.f90 index 43d07f86d..62a27b8f9 100644 --- a/src/gsi/cplr_regional_io.f90 +++ b/src/gsi/cplr_regional_io.f90 @@ -76,7 +76,7 @@ subroutine convert_regional_guess_wrf(this,mype,ctph0,stph0,tlm0) !$$$ end documentation block use kinds, only: i_kind,r_kind - use mpimod, only: mpi_integer4,mpi_rtype + use mpimod, only: mpi_integer4,mpi_rtype,mpi_itype use gridmod, only: wrf_mass_regional,wrf_nmm_regional,& nems_nmmb_regional,cmaq_regional,& twodvar_regional,netcdf @@ -88,6 +88,8 @@ subroutine convert_regional_guess_wrf(this,mype,ctph0,stph0,tlm0) use mpimod, only: mpi_comm_world,ierror use wrf_params_mod, only: update_pint,cold_start use gsi_io, only: verbose + use rapidrefresh_cldsurf_mod, only: i_howv_3dda, i_gust_3dda + use rapidrefresh_cldsurf_mod, only: i_howv_mask implicit none @@ -145,7 +147,15 @@ subroutine convert_regional_guess_wrf(this,mype,ctph0,stph0,tlm0) end if call mpi_barrier(mpi_comm_world,ierror) call mpi_bcast(byte_swap,1,mpi_integer4,0,mpi_comm_world,ierror) + call mpi_bcast(i_howv_3dda, 1, mpi_itype, 0, mpi_comm_world, ierror) + call mpi_bcast(i_gust_3dda, 1, mpi_itype, 0, mpi_comm_world, ierror) + call mpi_bcast(i_howv_mask, 1, mpi_itype, 0, mpi_comm_world, ierror) if(print_verbose)write(6,*)' in convert_regional_guess, for wrf arw binary input, byte_swap=',byte_swap + if (mype <= 1 .and. print_verbose) then + write(6,'(1x,A,3(2x,I4),2x,A6,I6.6,A2)') & + ' in convert_regional_guess, i_howv_3dda i_gust_3dda i_howv_mask =', & + i_howv_3dda,i_gust_3dda,i_howv_mask,' (pe=',mype,').' + end if elseif (cmaq_regional) then if (mype==0) then diff --git a/src/gsi/cplr_wrf_netcdf_interface.f90 b/src/gsi/cplr_wrf_netcdf_interface.f90 index 7c1bc400e..7cce1bd36 100644 --- a/src/gsi/cplr_wrf_netcdf_interface.f90 +++ b/src/gsi/cplr_wrf_netcdf_interface.f90 @@ -66,11 +66,14 @@ subroutine convert_netcdf_mass_wrf(this) use gsi_4dvar, only: nhr_assimilation use rapidrefresh_cldsurf_mod, only: l_hydrometeor_bkio,l_gsd_soilTQ_nudge use rapidrefresh_cldsurf_mod, only: i_use_2mt4b,i_use_2mq4b + use rapidrefresh_cldsurf_mod, only: i_howv_3dda, i_gust_3dda + use rapidrefresh_cldsurf_mod, only: i_howv_mask use gsi_metguess_mod, only: gsi_metguess_get use chemmod, only: laeroana_gocart, ppmv_conv,wrf_pm2_5 use gsi_chemguess_mod, only: gsi_chemguess_get use gridmod, only: wrf_mass_hybridcord use netcdf_mod, only: nc_check + use mpeu_util, only: die, perr use wrf_vars_mod, only : w_exist, dbz_exist use constants, only: zero @@ -91,6 +94,12 @@ subroutine convert_netcdf_mass_wrf(this) integer(i_kind) :: dh1 integer(i_kind) :: iunit + + character(len=120) :: filename_lakemask ! lakemask.bin + integer(i_kind) :: iunit_lakemask + logical :: l_iunit_opened + character(len=200) :: derr_msg ! error message used in die + character(len=24),parameter :: myname_ = 'convert_netcdf_mass_wrf' integer(i_kind) :: i,j,k integer(i_kind) :: ndim1 @@ -1022,6 +1031,116 @@ subroutine convert_netcdf_mass_wrf(this) write(iunit)field2 !TH2 endif +! Reading Significant Wave Height (HOWV) in firstguess + if(i_howv_3dda >= 1) then ! if howv is found in anavinfo + rmse_var='HOWV' + call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, & + start_index,end_index, WrfType, ierr ) + if(print_verbose)then + write(6,*)' rmse_var=',trim(rmse_var) + write(6,*)' ordering=',ordering + write(6,*)' WrfType,ierr=',WrfType,ierr + write(6,*)' ndim1=',ndim1 + write(6,*)' staggering=',staggering + write(6,*)' start_index=',start_index + write(6,*)' end_index=',end_index + end if + if(ierr == 0) then + call ext_ncd_read_field(dh1,DateStr1,TRIM(rmse_var), & + field2,WRF_REAL,0,0,0,ordering, & + staggering, dimnames , & + start_index,end_index, & !dom + start_index,end_index, & !mem + start_index,end_index, & !pat + ierr ) + if(print_verbose)write(6,*)' max,min ',trim(adjustl(rmse_var)),'=',maxval(field2),minval(field2) + write(iunit)field2 !HOWV (2D wave height) + i_howv_3dda = 2 ! howv was found both in anavinfo and firstguess + else + i_howv_3dda = 0 ! skipping analysis of this variable + derr_msg='Warning: Error when get info in firstguess for var : '//trim(rmse_var)// & + ' ==> Stop GSI analysis ... ierr=' +! write(6,'(1x,A,A2,A,1x,I6)')trim(adjustl(myname_)),'::',trim(adjustl(derr_msg)),ierr +! call stop2(ierr) + call die(trim(adjustl(myname_)),trim(adjustl(derr_msg)),ierr) + end if + +! reading lakemask for screening off wave height over lake (and/or land) area + if(i_howv_mask .eq. 2) then + rmse_var='LAKEMASK' ! see the code for rmse_var='XLAND' + call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, & + start_index,end_index, WrfType, ierr ) + if(ierr .eq. 0) then + if(print_verbose)then + write(6,*)' rmse_var = ',trim(rmse_var),' ndim1=',ndim1 + write(6,*)' WrfType = ',WrfType,' WRF_REAL=',WRF_REAL,'ierr = ',ierr !DEDE + write(6,*)' ordering = ',trim(ordering),' staggering = ',trim(staggering) + write(6,*)' start_index = ',start_index,' end_index = ',end_index + end if + call ext_ncd_read_field(dh1,DateStr1,TRIM(rmse_var), & + field2,WRF_REAL,0,0,0,ordering, & + staggering, dimnames , & + start_index,end_index, & !dom + start_index,end_index, & !mem + start_index,end_index, & !pat + ierr ) + if(print_verbose)then + write(6,*)' max,min lakemask=',maxval(field2),minval(field2) + write(6,*)' lakemask(1,1),lakemask(nlon,1)=',field2(1,1),field2(nlon_regional,1) + write(6,*)' lakemask(1,nlat),lakemask(nlon,nlat)=', & + field2(1,nlat_regional),field2(nlon_regional,nlat_regional) + end if +! dumping out lakemask into an temporary binary file +! (which would be used for howv mask when updating the analysis of howv) + filename_lakemask="lakemask.bin"; iunit_lakemask=115; l_iunit_opened=.true. + inquire(unit=iunit_lakemask, opened=l_iunit_opened) + if(l_iunit_opened) iunit_lakemask=iunit_lakemask+5000 + if(print_verbose) write(6,*)' I/O unit for writing lakemask: ', iunit_lakemask + open(iunit_lakemask,file=trim(adjustl(filename_lakemask)),form='unformatted') + write(iunit_lakemask)field2 !LAKEMASK (1=lake, 0=water) + close(iunit_lakemask) + else + write(6,'(1x,A)')' No LAKEMASK data in firstguess data. No lakemask is used for HOWV mask. reset i_howv_mask=1' + i_howv_mask=1 ! i_howv_mask is changed on pe 0, so it must be broadcast to all PEs. + endif ! if error in reading lakemask + endif ! i_howv_mask = 2 + endif ! i_howv_3dda (reading 2D wave height from netcdf-format background) + +! Reading 10-m Wind Gust (GUST) in firstguess + if(i_gust_3dda >= 1) then ! if gust is found in anavinfo + rmse_var='GUST' + call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, & + start_index,end_index, WrfType, ierr ) + if(print_verbose)then + write(6,*)' rmse_var=',trim(rmse_var) + write(6,*)' ordering=',ordering + write(6,*)' WrfType,ierr=',WrfType,ierr + write(6,*)' ndim1=',ndim1 + write(6,*)' staggering=',staggering + write(6,*)' start_index=',start_index + write(6,*)' end_index=',end_index + end if + if(ierr == 0) then + call ext_ncd_read_field(dh1,DateStr1,TRIM(rmse_var), & + field2,WRF_REAL,0,0,0,ordering, & + staggering, dimnames , & + start_index,end_index, & !dom + start_index,end_index, & !mem + start_index,end_index, & !pat + ierr ) + if(print_verbose)write(6,*)'convert_netcdf_mass_wrf:: max,min ',trim(adjustl(rmse_var)),'=',maxval(field2),minval(field2) + write(iunit)field2 !GUST (2D Wind Gust) + i_gust_3dda = 2 ! gust was found both in anavinfo and firstguess + else + i_gust_3dda = 0 ! skipping analysis of this variable + derr_msg='Warning: Error when get info in firstguess for var : '//trim(rmse_var)// & + ' ==> Stop GSI analysis ... ierr=' +! write(6,'(1x,A,A2,A,1x,I6)')trim(adjustl(myname_)),'::',trim(adjustl(derr_msg)),ierr +! call stop2(ierr) + call die(trim(adjustl(myname_)),trim(adjustl(derr_msg)),ierr) + end if + endif ! i_gust_3dda (reading 2D 10-m wind gust from netcdf-format background) + if(l_hydrometeor_bkio .and. n_actual_clouds>0) then rmse_var='QCLOUD' call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, & @@ -2437,6 +2556,7 @@ subroutine update_netcdf_mass_wrf(this) use constants, only: h300,tiny_single use rapidrefresh_cldsurf_mod, only: l_hydrometeor_bkio,l_gsd_soilTQ_nudge use rapidrefresh_cldsurf_mod, only: i_gsdcldanal_type + use rapidrefresh_cldsurf_mod, only: i_howv_3dda, i_gust_3dda use gsi_metguess_mod, only: gsi_metguess_get,GSI_MetGuess_Bundle use rapidrefresh_cldsurf_mod, only: i_use_2mt4b,i_use_2mq4b use gsi_bundlemod, only: GSI_BundleGetPointer @@ -3002,6 +3122,61 @@ subroutine update_netcdf_mass_wrf(this) ierr ) endif +! writing analysis of significantg wave height (howv) into firstguess (analysis) file + if(i_howv_3dda==2) then + read(iunit) field2 !HOWV (wave height) + rmse_var='HOWV' + if(print_verbose)write(6,'(1x,4A,2(1x,F15.7))') & + trim(adjustl(myname_)),'::max,min ',trim(rmse_var),'=', & + maxval(field2),minval(field2) + call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, & + start_index,end_index1, WrfType, ierr ) + if(print_verbose)then + write(6,*)' rmse_var=',trim(rmse_var) + write(6,*)' ordering=',ordering + write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL + write(6,*)' ndim1=',ndim1 + write(6,*)' staggering=',staggering + write(6,*)' start_index=',start_index + write(6,*)' end_index1=',end_index1 + end if + call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var), & + field2,WRF_REAL,0,0,0,ordering, & + staggering, dimnames , & + start_index,end_index1, & !dom + start_index,end_index1, & !mem + start_index,end_index1, & !pat + ierr ) + endif ! i_howv_3dda (wave height) + +! writing analysis of 10-m wind gust (gust) into firstguess (analysis) file + if(i_gust_3dda==2) then + read(iunit) field2 !GUST (wind gust) + rmse_var='GUST' + if(print_verbose)write(6,'(1x,4A,2(1x,F15.7))') & + trim(adjustl(myname_)),'::max,min ',trim(rmse_var),'=', & + maxval(field2),minval(field2) + call ext_ncd_get_var_info (dh1,trim(rmse_var),ndim1,ordering,staggering, & + start_index,end_index1, WrfType, ierr ) + if(print_verbose)then + write(6,*)' rmse_var=',trim(rmse_var) + write(6,*)' ordering=',ordering + write(6,*)' WrfType,WRF_REAL=',WrfType,WRF_REAL + write(6,*)' ndim1=',ndim1 + write(6,*)' staggering=',staggering + write(6,*)' start_index=',start_index + write(6,*)' end_index1=',end_index1 + end if + call ext_ncd_write_field(dh1,DateStr1,TRIM(rmse_var), & + field2,WRF_REAL,0,0,0,ordering, & + staggering, dimnames , & + start_index,end_index1, & !dom + start_index,end_index1, & !mem + start_index,end_index1, & !pat + ierr ) + endif ! i_howv_3dda (wave height) + + if (l_hydrometeor_bkio .and. n_actual_clouds>0) then do k=1,nsig_regional read(iunit)((field3(i,j,k),i=1,nlon_regional),j=1,nlat_regional) ! Qc diff --git a/src/gsi/cplr_wrwrfmassa.f90 b/src/gsi/cplr_wrwrfmassa.f90 index 6f416c6bd..0bebf1947 100644 --- a/src/gsi/cplr_wrwrfmassa.f90 +++ b/src/gsi/cplr_wrwrfmassa.f90 @@ -1873,8 +1873,12 @@ subroutine wrwrfmassa_netcdf_wrf(this,mype) use constants, only: one,zero_single,rd_over_cp_mass,one_tenth,r10,r100 use constants, only: soilmoistmin use gsi_io, only: lendian_in, lendian_out + use gsi_io, only: verbose use rapidrefresh_cldsurf_mod, only: l_hydrometeor_bkio,l_gsd_soilTQ_nudge,& i_use_2mq4b,i_use_2mt4b + use rapidrefresh_cldsurf_mod, only: i_howv_3dda, i_gust_3dda + use rapidrefresh_cldsurf_mod, only: i_howv_mask + use chemmod, only: laeroana_gocart,wrf_pm2_5 use gsi_bundlemod, only: GSI_BundleGetPointer use gsi_metguess_mod, only: gsi_metguess_get,GSI_MetGuess_Bundle @@ -1913,6 +1917,7 @@ subroutine wrwrfmassa_netcdf_wrf(this,mype) integer(i_kind) i_qc,i_qi,i_qr,i_qs,i_qg,i_qnr,i_qni,i_qnc integer(i_kind) kqc,kqi,kqr,kqs,kqg,kqnr,kqni,kqnc,i_tt,ktt,kw,kdbz integer(i_kind) i_sst,i_skt,i_th2,i_q2,i_soilt1,i_tslb,i_smois,ktslb,ksmois + integer(i_kind) i_howv, i_gust integer(i_kind) :: iv, n_gocart_var,i_snowT_check integer(i_kind),allocatable :: i_chem(:), kchem(:) integer(i_kind) num_mass_fields,num_all_fields,num_all_pad,num_mass_fields_base @@ -1926,11 +1931,19 @@ subroutine wrwrfmassa_netcdf_wrf(this,mype) real(r_single) aeta10(nsig),eta10(nsig+1),aeta20(nsig),eta20(nsig+1) real(r_single) glon0(nlon_regional,nlat_regional),glat0(nlon_regional,nlat_regional) real(r_single) dx_mc0(nlon_regional,nlat_regional),dy_mc0(nlon_regional,nlat_regional) + +!--- used to read lakemask for screen off analysis over lake area (e.g., for wave height HOWV) + character(len=20) :: filename_lakemask ! lakemask.bin + integer(i_kind) :: iunit_lakemask + logical :: l_iunit_opened + real(r_single),allocatable::lakemask(:) real(r_kind), pointer :: ges_ps(:,: )=>NULL() real(r_kind), pointer :: ges_tsk(:,:)=>NULL() real(r_kind), pointer :: ges_t2m(:,:)=>NULL() real(r_kind), pointer :: ges_q2(:,:)=>NULL() + real(r_kind), pointer :: ges_howv_it(:,:)=>NULL() ! wave height + real(r_kind), pointer :: ges_gust_it(:,:)=>NULL() ! wind gust real(r_kind), pointer :: ges_soilt1(:,:)=>NULL() real(r_kind), pointer :: ges_tslb_it(:,:,:)=>NULL() real(r_kind), pointer :: ges_smois_it(:,:,:)=>NULL() @@ -1965,6 +1978,11 @@ subroutine wrwrfmassa_netcdf_wrf(this,mype) real(r_kind), pointer :: ges_p25 (:,:,:)=>NULL() real(r_kind), pointer :: ges_pm2_5 (:,:,:)=>NULL() + logical :: print_verbose + + print_verbose = .false. + if(verbose .and. mype == 0)print_verbose=.true. + associate( this => this ) ! eliminates warning for unused dummy argument needed for binding end associate it=ntguessig @@ -2003,6 +2021,8 @@ subroutine wrwrfmassa_netcdf_wrf(this,mype) if(l_gsd_soilTQ_nudge) num_mass_fields=num_mass_fields+2*nsig_soil+1 if(i_use_2mt4b > 0 ) num_mass_fields=num_mass_fields+2 if(i_use_2mt4b <= 0 .and. i_use_2mq4b > 0) num_mass_fields=num_mass_fields+1 + if(i_howv_3dda > 0) num_mass_fields = num_mass_fields + 1 + if(i_gust_3dda > 0) num_mass_fields = num_mass_fields + 1 if ( laeroana_gocart ) then call gsi_chemguess_get ( 'aerosols::3d', n_gocart_var, ier ) if ( n_gocart_var > 0 ) then @@ -2060,10 +2080,22 @@ subroutine wrwrfmassa_netcdf_wrf(this,mype) else i_q2=i_skt endif +! howv is after q2, and before gust + if ( i_howv_3dda > 0 ) then + i_howv=i_q2+1 + else + i_howv=i_q2 + end if +! gust is after howv, and before cloud hydrometers + if ( i_gust_3dda > 0 ) then + i_gust=i_howv+1 + else + i_gust=i_howv + end if ! for hydrometeors if(l_hydrometeor_bkio .and. n_actual_clouds>0) then - i_qc=i_q2+1 + i_qc=i_gust+1 i_qr=i_qc+lm i_qs=i_qr+lm i_qi=i_qs+lm @@ -2097,13 +2129,13 @@ subroutine wrwrfmassa_netcdf_wrf(this,mype) else if ( laeroana_gocart) then do iv = 1, n_gocart_var - i_chem(iv)=i_skt+(iv-1)*lm+1 + i_chem(iv)=i_gust+(iv-1)*lm+1 end do endif if ( wrf_pm2_5 ) then iv=1 - i_chem(iv)=i_skt+(iv-1)*lm+1 + i_chem(iv)=i_gust+(iv-1)*lm+1 endif @@ -2594,6 +2626,32 @@ subroutine wrwrfmassa_netcdf_wrf(this,mype) end do end do endif + ! Wave height (howv) + if ( i_howv_3dda > 0 ) then + call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'howv', ges_howv_it, istatus );ier=ier+istatus + if (ier/=0) then ! doesn't have to die - code can be generalized to bypass missing vars + write(6,*)'wrwrfmassa_netcdf_wrf: getpointer for howv failed, cannot retrieve howv, ier =',ier + call stop2(999) + end if + do i=1,lon2 + do j=1,lat2 + all_loc(j,i,i_howv)=ges_howv_it(j,i) + end do + end do + end if + ! Wind gust (gust) + if ( i_gust_3dda > 0 ) then + call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'gust', ges_gust_it, istatus );ier=ier+istatus + if (ier/=0) then ! doesn't have to die - code can be generalized to bypass missing vars + write(6,*)'wrwrfmassa_netcdf_wrf: getpointer for gust failed, cannot retrieve gust, ier =',ier + call stop2(999) + end if + do i=1,lon2 + do j=1,lat2 + all_loc(j,i,i_gust)=ges_gust_it(j,i) + end do + end do + end if if(mype == 0) then ! SM This is landmask @@ -2803,6 +2861,116 @@ subroutine wrwrfmassa_netcdf_wrf(this,mype) write(lendian_out)temp1 end if !endif mype==0 endif ! i_use_2mt4b>0 + + ! update HOWV (wave height): after TSK/Q2/SOILT1/TH2 + if ( i_howv_3dda > 0 ) then + + !--- reading lakemask if existing + if(mype == 0) then + allocate(lakemask(im*jm)) + lakemask(:) = -1.0_r_single + if(i_howv_mask == 2) then + filename_lakemask="lakemask.bin" + iunit_lakemask=115 + l_iunit_opened=.true. + inquire(unit=iunit_lakemask, opened=l_iunit_opened) + if(l_iunit_opened) iunit_lakemask=iunit_lakemask+5000 + if(print_verbose) write(6,*)' I/O unit for reading lakemask: ', iunit_lakemask + open(iunit_lakemask,file=trim(adjustl(filename_lakemask)),form='unformatted',status='old') + read(iunit_lakemask)temp1 !LAKEMASK (1=lake, 0=water) + lakemask = temp1 + close(iunit_lakemask) + end if + !--- checking the land and lake mask + if ( i_howv_mask > 0 ) then + if (maxval(landmask) > 1.01_r_single .or. minval(landmask) < -0.01_r_single) then + i_howv_mask = 0 + if (print_verbose) then + write(6,*) ' bad landmask, do not use land and lake mask to filter howv field.' + endif + else + if (maxval(lakemask) > 1.01_r_single .or. minval(lakemask) < -0.01_r_single) then + i_howv_mask = 1 + if (print_verbose) then + write(6,*) ' bad lakemask, do not use it, only use land mask to filter howv field.' + endif + endif + endif + else + i_howv_mask = 0 + if (print_verbose) then + write(6,*) ' No land and lake mask would be used to filter howv field.' + endif + endif + endif ! mype = 0 (reading mask data only on pe=0) + + if(mype == 0) read(lendian_in)temp1 + if(mype == 0) write(6,*)' at 10.15 (howv: fgs) in wrwrfmassa,max,min(temp1)=', & + maxval(temp1),minval(temp1) + call strip(all_loc(:,:,i_howv),strp) + tempa=zero_single + call mpi_gatherv(strp,ijn(mype+1),mpi_real4, & + tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror) + if(mype == 0) then + write(6,*)' at 10.16 (howv: anl on analysis grid) in wrwrfmassa,max,min(tempa)=', & + maxval(tempa),minval(tempa) + call fill_mass_grid2t(temp1,im,jm,tempb,2) + do i=1,iglobal + tempa(i)=tempa(i)-tempb(i) + end do + write(6,*)' at 10.17 (howv: inc on analysis grid) in wrwrfmassa,max,min(tempa)=', & + maxval(tempa),minval(tempa) + + if ( i_howv_mask > 0 ) then + ! --- mask/screen off the howv analysis over land and lake area + if (print_verbose) write(6,*) & + ' Analysis of howv is masked/modified by landmask/lakemask.' + call unfill_mass_grid2t_ldlkmask(tempa,im,jm,temp1,landmask, & + lakemask,i_howv_mask) + else + if (print_verbose) write(6,*) & + ' Analysis of howv is NOT masked by landmask/lakemask.' + call unfill_mass_grid2t(tempa,im,jm,temp1) + end if + write(6,*)' at 10.18 (howv: anl on model grid) in wrwrfmassa,max,min(temp1)=', & + maxval(temp1),minval(temp1) + write(lendian_out)temp1 + + deallocate(lakemask) + + end if !endif mype==0 + + end if ! i_howv_3dda > 0 + + ! update GUST (wind gust): after howv + if ( i_gust_3dda > 0 ) then + + if(mype == 0) read(lendian_in)temp1 + if(mype == 0) write(6,*)' at 10.19 (gust: fgs) in wrwrfmassa,max,min(temp1)=', & + maxval(temp1),minval(temp1) + call strip(all_loc(:,:,i_gust),strp) + tempa=zero_single + call mpi_gatherv(strp,ijn(mype+1),mpi_real4, & + tempa,ijn,displs_g,mpi_real4,0,mpi_comm_world,ierror) + if(mype == 0) then + write(6,*)' at 10.20 (gust: anl on analysis grid) in wrwrfmassa,max,min(tempa)=', & + maxval(tempa),minval(tempa) + call fill_mass_grid2t(temp1,im,jm,tempb,2) + do i=1,iglobal + tempa(i)=tempa(i)-tempb(i) + end do + write(6,*)' at 10.21 (gust: inc on analysis grid) in wrwrfmassa,max,min(tempa)=', & + maxval(tempa),minval(tempa) + + call unfill_mass_grid2t(tempa,im,jm,temp1) + write(6,*)' at 10.22 (gust: anl on model grid) in wrwrfmassa,max,min(temp1)=', & + maxval(temp1),minval(temp1) + write(lendian_out)temp1 + + end if !endif mype==0 + + end if ! i_gust_3dda > 0 + ! ! for saving cloud analysis results if(l_hydrometeor_bkio .and. n_actual_clouds>0) then diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index 8bae3ca61..754d14499 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -184,7 +184,7 @@ module gsimod cld_bld_coverage,cld_clr_coverage,& i_cloud_q_innovation,i_ens_mean,DTsTmax,& i_T_Q_adjust,l_saturate_bkCloud,l_rtma3d,i_precip_vertical_check, & - corp_howv, hwllp_howv, corp_gust, hwllp_gust, oerr_gust + corp_howv, hwllp_howv, corp_gust, hwllp_gust, oerr_gust, i_howv_mask use gsi_metguess_mod, only: gsi_metguess_init,gsi_metguess_final use gsi_chemguess_mod, only: gsi_chemguess_init,gsi_chemguess_final use tcv_mod, only: init_tcps_errvals,tcp_refps,tcp_width,tcp_ermin,tcp_ermax @@ -1605,6 +1605,10 @@ module gsimod ! = 0.42 meters (default) ! hwllp_howv - real, background error de-correlation length scale of howv ! = 170,000.0 meters (default 170 km) +! i_howv_mask - integer, option to control the mask of the wave height (howv) over the land/lake area +! = 0: do not mask (default) +! = 1: mask the value over the land area only (using land mask data) +! = 2: mask the value over the land & lake area (using land mask and lake mask data) ! corp_gust - real, static background error of gust (stddev error) ! hwllp_gust - real, background error de-correlation length scale of gust ! oerr_gust - real, observation error of gust @@ -1629,7 +1633,7 @@ module gsimod cld_bld_coverage,cld_clr_coverage,& i_cloud_q_innovation,i_ens_mean,DTsTmax, & i_T_Q_adjust,l_saturate_bkCloud,l_rtma3d,i_precip_vertical_check, & - corp_howv, hwllp_howv, corp_gust, hwllp_gust, oerr_gust + corp_howv, hwllp_howv, corp_gust, hwllp_gust, oerr_gust, i_howv_mask ! chem(options for gsi chem analysis) : ! berror_chem - .true. when background for chemical species that require diff --git a/src/gsi/m_berror_stats_reg.f90 b/src/gsi/m_berror_stats_reg.f90 index d7a30808e..f87180e65 100644 --- a/src/gsi/m_berror_stats_reg.f90 +++ b/src/gsi/m_berror_stats_reg.f90 @@ -1167,7 +1167,7 @@ subroutine read_howv_stats(nlat,nlon,npar,arrout,mype) arrout(:,:,2)=50000.0_r_kind ! values were specified by Manuel and Stelio for 2DRTMA else arrout(:,:,1) = corp_howv ! 0.42_r_kind used in 3dvar (default) if not set in namelist - arrout(:,:,2) = hwllp_howv ! 17000.0_r_kind used in 3dvar (default) if not set in namelist + arrout(:,:,2) = hwllp_howv ! 170,000.0_r_kind used in 3dvar (default) if not set in namelist end if reclength=nlat*r_kind diff --git a/src/gsi/rapidrefresh_cldsurf_mod.f90 b/src/gsi/rapidrefresh_cldsurf_mod.f90 index 475f44a9d..d0c0a0b8e 100644 --- a/src/gsi/rapidrefresh_cldsurf_mod.f90 +++ b/src/gsi/rapidrefresh_cldsurf_mod.f90 @@ -185,18 +185,18 @@ module rapidrefresh_cldsurf_mod ! = 2(clean Qg as in 1, and adjustment to the retrieved Qr/Qs/Qnr throughout the whole profile) ! = 3(similar to 2, but adjustment to Qr/Qs/Qnr only below maximum reflectivity level ! and where the dbz_obs is missing); +! ! corp_howv - namelist real, static BE of howv (standard error deviation) ! hwllp_howv - namelist real, static BE de-correlation length scale of howv ! i_howv_3dda - integer, control the analysis of howv in 3D analysis (either var or hybrid) ! = 0 (howv-off: default) : no analysis of howv in 3D analysis. ! = 1 (howv-on) : if variable name "howv" is found in anavinfo, ! set it to be 1 to turn on analysis of howv; -! note: in hybrid envar run, the static BE is redueced by beta_s (<1.0), -! since there is no ensemble of howv currently yet, then no ensemble -! contribution to the total BE of howv, so the total BE of howv is actually -! just the reduced static BE of howv. If to make the analysis of howv -! in hyrbid run is as similar as the analysis of howv in pure 3dvar run, -! the static BE of howv used in hybrid run needs to be tuned (inflated actually). +! i_howv_mask - integer, control the screen of wave height (howv) over the land and lake area +! with the land mask and lake mask data +! = 0 (default): do not screen off the howv data +! = 1 : screen of the howv data over the land area only (with land mask) +! = 2 : screen of the howv data over the land and lake area (with land mask and lake mask) ! corp_gust - namelist real, static BE of gust (standard error deviation) ! note: 1. initialised to be an arbitary negative value, in order to skip this ! negative value, instead to use value (3.0 m/s) set in subroutine @@ -294,6 +294,7 @@ module rapidrefresh_cldsurf_mod public :: i_precip_vertical_check public :: corp_howv, hwllp_howv public :: i_howv_3dda + public :: i_howv_mask public :: corp_gust, hwllp_gust, oerr_gust public :: i_gust_3dda @@ -356,6 +357,7 @@ module rapidrefresh_cldsurf_mod integer(i_kind) i_precip_vertical_check real(r_kind) :: corp_howv, hwllp_howv integer(i_kind) :: i_howv_3dda + integer(i_kind) :: i_howv_mask real(r_kind) :: corp_gust, hwllp_gust, oerr_gust integer(i_kind) :: i_gust_3dda @@ -475,6 +477,7 @@ subroutine init_rapidrefresh_cldsurf corp_howv = 0.42_r_kind ! 0.42 meters (default) hwllp_howv = 170000.0_r_kind ! 170,000.0 meters (170km as default for 3DRTMA, 50km is used in 2DRTMA) i_howv_3dda = 0 ! no analysis of significant wave height (howv) in 3D analysis (default) + i_howv_mask = 0 ! do NOT mask significant wave height (howv) over the land/lake area (default) corp_gust = -1.50_r_kind ! initialised as negative & void to be skipped, in order to use ! the value (3.0 m/s) set in sub berror_read_wgt_reg (as default). ! If user-specified value is preferred, set it in session diff --git a/src/gsi/setuphowv.f90 b/src/gsi/setuphowv.f90 index 3ecb05c8f..053ff85c6 100644 --- a/src/gsi/setuphowv.f90 +++ b/src/gsi/setuphowv.f90 @@ -68,6 +68,7 @@ subroutine setuphowv(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diag lobsdiagsave,nobskeep,lobsdiag_allocated, & time_offset,bmiss,luse_obsdiag,ianldate use obsmod, only: netcdf_diag, binary_diag, dirname + use obsmod, only: l_obsprvdiag 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 @@ -230,7 +231,7 @@ subroutine setuphowv(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diag ioff0=19 nreal=ioff0 if (lobsdiagsave) nreal=nreal+4*miter+1 - if (twodvar_regional) then; nreal=nreal+2; allocate(cprvstg(nobs),csprvstg(nobs)); endif + if (twodvar_regional .or. l_obsprvdiag) then; nreal=nreal+2; allocate(cprvstg(nobs),csprvstg(nobs)); endif allocate(cdiagbuf(nobs),rdiagbuf(nreal,nobs)) if (netcdf_diag) call init_netcdf_diag_ end if @@ -443,7 +444,7 @@ subroutine setuphowv(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diag write(7)cdiagbuf(1:ii),rdiagbuf(:,1:ii) deallocate(cdiagbuf,rdiagbuf) - if (twodvar_regional) then + if (twodvar_regional .or. l_obsprvdiag) then write(7)cprvstg(1:ii),csprvstg(1:ii) deallocate(cprvstg,csprvstg) endif @@ -623,7 +624,7 @@ subroutine contents_binary_diag_(odiag) enddo endif - if (twodvar_regional) then + if (twodvar_regional .or. l_obsprvdiag) then rdiagbuf(ioff+1,ii) = data(idomsfc,i) ! dominate surface type rdiagbuf(ioff+2,ii) = data(izz,i) ! model terrain at observation location r_prvstg = data(iprvd,i) @@ -680,7 +681,7 @@ subroutine contents_netcdf_diag_(odiag) call nc_diag_data2d("ObsDiagSave_obssen", odiag%obssen ) endif - if (twodvar_regional) then + if (twodvar_regional .or. l_obsprvdiag) 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) diff --git a/src/gsi/unfill_mass_grid2.f90 b/src/gsi/unfill_mass_grid2.f90 index 317bf24d0..dca2ea99c 100644 --- a/src/gsi/unfill_mass_grid2.f90 +++ b/src/gsi/unfill_mass_grid2.f90 @@ -381,3 +381,101 @@ subroutine unfill_mass_grid2t_drycheck(gout,nx,ny,gin,qs) end do end subroutine unfill_mass_grid2t_drycheck + +subroutine unfill_mass_grid2t_ldlkmask(gout,nx,ny,gin,landmask, & + lakemask,i_howv_mask) +!$$$ subprogram documentation block +! . . . . +! subprogram: unfill_mass_grid2t_ldlkmask +! prgmmr: zhao org: saic@ncep/emc date: 2025-01-28 +! +! abstract: This subroutine is based on subroutine unfill_mass_grid2t andi +! unfill_mass_grid2t_ldmk, but screen off analysis increment over +! land and/or lake area with land-mask and lake-mask. +! +! +! program history log: +! 2025-01-28 zhao - apply land and lake mask for howv (wave height) fields +! +! input argument list: +! gout - input A-grid (reorganized for distibution to local domains) +! gin - preexisting input values to be added to on C-grid +! nx,ny - input grid dimensions +! landmask - land mask (0: water, 1: land) +! lakemask - lake mask (0: sea and land (including Great Lakes), 1: lake) +! i_howv_mask - input option for HOWV adjustment +! = 0: no adjustment +! = 1: mask off data over land +! = 2: mask off data over land and lake +! +! output argument list: +! gin - output result on C grid +! +! attributes: +! language: f90 +! machine: ibm RS/6000 SP +! +!$$$ + use kinds, only: r_single,i_kind + use gridmod, only: itotsub,iglobal + use general_commvars_mod, only: ltosi,ltosj + use mod_wrfmass_to_a, only: wrfmass_a_to_h4 + use gridmod, only: nlon, nlat + + implicit none + + integer(i_kind), intent(in ) :: nx,ny + real(r_single) , intent(in ) :: gout(itotsub) + real(r_single) , intent(inout) :: gin(nx,ny) + real(r_single) , intent(in) :: landmask(nx,ny) + real(r_single) , intent(in) :: lakemask(nx,ny) + integer(i_kind), intent(in ) :: i_howv_mask + + real(r_single) ba(nlon,nlat) + real(r_single) b(nx,ny) + integer(i_kind) i,j + + do i=1,iglobal + ba(ltosj(i),ltosi(i))=gout(i) + end do + + if(nlon == nx .and. nlat == ny) then + b=ba + else + call wrfmass_a_to_h4(ba,b) + endif + +! Only add analysis increment over sea water (no land and lake) + select case (i_howv_mask) + case(2) + write(6,*) 'unfill_mass_grid2t_ldlkmask: using landlakemask to screen off HOWV analysis increment over land and lake.' + do j=1,ny + do i=1,nx + if (landmask(i,j) > 0.00001_r_single .or. lakemask(i,j) > 0.00001_r_single) then + b(i,j) = 0.0_r_single + end if + end do + end do + case(1) + write(6,*) 'unfill_mass_grid2t_ldlkmask: using landmask to screen off HOWV analysis increment over land.' + do j=1,ny + do i=1,nx + if (landmask(i,j) > 0.00001_r_single) then + b(i,j) = 0.0_r_single + end if + end do + end do + case default + write(6,*) 'unfill_mass_grid2t_ldlkmask: no land iand lake mask used to screen off HOWV analysis increment.' + end select + +! Mass grids--just copy + do j=1,ny + do i=1,nx + gin(i,j)=b(i,j)+gin(i,j) + end do + end do + + return + +end subroutine unfill_mass_grid2t_ldlkmask