Adding I/O for direct analysis of near-surface wind gust for RRFS-bas…

…ed 3DRTMA (NOAA-EMC#730)   **Description**    To improve the analysis of the near-surface wind gust in 3DRTMA, the observations of near-surface wind gust would be analyzed directly in GSI (3DVar and Hybrid 3DEnVar), instead of being a derived product from the near-surface wind analysis. Since the core subroutines for direct variational assimilation of wind gust (e.g., setupgust.f90, intgust.f90, stpgust.f90, etc.) had already been implemented in GSI for 2DRTMA, so in the work the development in GSI mainly focus on adding I/O of 2-D wind gust firstguess and analysis fields for RRFS-based 3DRTMA, and some minor modifications in observation and background error for wind gust, options to control the analysis of wind gust, etc.   This PR is to address the issue NOAA-EMC#726 : Adding I/O for direct analysis of near-surface wind gust for RRFS-based 3DRTMA Fixes NOAA-EMC#726 **Type of change** Please delete options that are not relevant. - [x] New feature (non-breaking change which adds functionality) **How Has This Been Tested?**    **Checklist** - [x] My code follows the style guidelines of this project - [x] I have performed a self-review of my own code - [x] I have commented my code, particularly in hard-to-understand areas - [x] New and existing tests pass with my changes tested with a real case - 2024-02-20_12:00Z, 1. dry-run (using my updated GSI code with wind gust analysis, but actually no wind gust obs is analyzed, so-called dryrun) is compared to control-run (original GSI code running without wind gust obs): the results are identical. This indicates that if without analyzing wind-gust obs, then the updated code generates the analysis identical to the analysis of original/control code. Or say, the added code does not have influence on the other part of code. 2. real case run with updated GSI code to analyze the obs of wind gust: The following figure shows the used observations of near-surface wind gust: ![var_obs_2024022012_gust_used_maprll_datll_reg_ncf](https://github.com/NOAA-EMC/GSI/assets/53267411/ecbe479a-03c6-490f-a179-9e0027291468) the following figure shows the analysis increments: ![GUST_hyb_hwllp90_corptuned_inc_incrintrp_maprll_datrll_reg_grb2](https://github.com/NOAA-EMC/GSI/assets/53267411/a01fca0d-dc1f-438b-b8eb-e624de35a631) - [x] Any dependent changes have been merged and published - [x] Regression tests on WCOSS2 (Cactus) and Hera (Rocky-8) : my updated GSI commit [#f91f247d](GangZhao-NOAA@f91f247)) vs control/original GSI code (commit [#6d9ebbb7](NOAA-EMC@6d9ebbb)) Here is the reports of the regression tests on WCOSS2 (Cactus): ~~~ [gang.zhao@clogin02:build] (feature/windgust_in_3dvar_for_3drtma)$ ctest -j 7 Test project /lfs/h2/emc/da/save/gang.zhao/WorkDir/ProdGSI_Dev/gsi_dev/build Start 1: global_4denvar Start 2: rtma Start 3: rrfs_3denvar_glbens Start 4: netcdf_fv3_regional Start 5: hafs_4denvar_glbens Start 6: hafs_3denvar_hybens Start 7: global_enkf 1/7 Test NOAA-EMC#4: netcdf_fv3_regional .............. Passed 483.15 sec 2/7 Test NOAA-EMC#3: rrfs_3denvar_glbens .............. Passed 486.74 sec 3/7 Test NOAA-EMC#7: global_enkf ...................... Passed 850.98 sec 4/7 Test #2: rtma ............................. Passed 970.28 sec 5/7 Test NOAA-EMC#6: hafs_3denvar_hybens .............. Passed 1152.82 sec 6/7 Test NOAA-EMC#5: hafs_4denvar_glbens .............. Passed 1213.93 sec 7/7 Test #1: global_4denvar ................... Passed 1683.16 sec 100% tests passed, 0 tests failed out of 7 Total Test time (real) = 1683.19 sec ~~~ Here is the reports of the regression tests on Hera (Rocky8): ~~~ (base) [Gang.Zhao@hfe11:build] (feature/windgust_in_3dvar_for_3drtma)$ ctest -j 7 Test project /scratch1/NCEPDEV/da/Gang.Zhao/ProdGSI_dev/gsi_dev/build Start 1: global_4denvar Start 2: rtma Start 3: rrfs_3denvar_glbens Start 4: netcdf_fv3_regional Start 5: hafs_4denvar_glbens Start 6: hafs_3denvar_hybens Start 7: global_enkf 1/7 Test NOAA-EMC#4: netcdf_fv3_regional .............. Passed 491.53 sec 2/7 Test NOAA-EMC#3: rrfs_3denvar_glbens ..............***Failed 495.27 sec 3/7 Test #2: rtma ............................. Passed 982.45 sec 4/7 Test NOAA-EMC#6: hafs_3denvar_hybens .............. Passed 1168.99 sec 5/7 Test NOAA-EMC#7: global_enkf ...................... Passed 1239.77 sec 6/7 Test NOAA-EMC#5: hafs_4denvar_glbens ..............***Failed 1347.87 sec 7/7 Test #1: global_4denvar ................... Passed 1974.45 sec 71% tests passed, 2 tests failed out of 7 Total Test time (real) = 1974.91 sec The following tests FAILED: 3 - rrfs_3denvar_glbens (Failed) 5 - hafs_4denvar_glbens (Failed) Errors while running CTest Output from these tests are in: /scratch1/NCEPDEV/da/Gang.Zhao/ProdGSI_dev/gsi_dev/build/Testing/Temporary/LastTest.log Use "--rerun-failed --output-on-failure" to re-run the failed cases verbosely. (base) [Gang.Zhao@hfe11:build] (feature/windgust_in_3dvar_for_3drtma)$ ctest -R rrfs_3denvar_glbens Test project /scratch1/NCEPDEV/da/Gang.Zhao/ProdGSI_dev/gsi_dev/build Start 3: rrfs_3denvar_glbens 1/1 Test NOAA-EMC#3: rrfs_3denvar_glbens .............. Passed 430.52 sec 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 430.55 sec (base) [Gang.Zhao@hfe11:build] (feature/windgust_in_3dvar_for_3drtma)$ ctest -R hafs_4denvar_glbens Test project /scratch1/NCEPDEV/da/Gang.Zhao/ProdGSI_dev/gsi_dev/build Start 5: hafs_4denvar_glbens 1/1 Test NOAA-EMC#5: hafs_4denvar_glbens .............. Passed 1225.37 sec 100% tests passed, 0 tests failed out of 1 Total Test time (real) = 1225.39 sec ~~~ **Note**: _When I was running the regression tests, GSI code was just updated to the latest commit [#b53740a7](GangZhao-NOAA@f91f247). Considering the frequent update in EMC GSI code recently and saving the time, after this PR has been reviewed and approved by peer-reviewers, I will update the code to latest EMC GSI commit, then re-run the regression tests for final approval.
TingLei-daprediction · Apr 5, 2024 · 29d9d8f · 29d9d8f
1 parent db477e3
commit 29d9d8f
Show file tree

Hide file tree

Showing 5 changed files with 140 additions and 16 deletions.
diff --git a/src/gsi/gsi_rfv3io_mod.f90 b/src/gsi/gsi_rfv3io_mod.f90
@@ -61,7 +61,7 @@ module gsi_rfv3io_mod
   use rapidrefresh_cldsurf_mod, only: i_use_2mq4b,i_use_2mt4b
   use chemmod, only: naero_cmaq_fv3,aeronames_cmaq_fv3,imodes_cmaq_fv3,laeroana_fv3cmaq
   use chemmod, only: naero_smoke_fv3,aeronames_smoke_fv3,laeroana_fv3smoke  
-  use rapidrefresh_cldsurf_mod, only: i_howv_3dda
+  use rapidrefresh_cldsurf_mod, only: i_howv_3dda, i_gust_3dda
 
   implicit none
   public type_fv3regfilenameg
@@ -147,7 +147,7 @@ module gsi_rfv3io_mod
   public :: mype_u,mype_v,mype_t,mype_q,mype_p,mype_oz,mype_ql
   public :: mype_qi,mype_qr,mype_qs,mype_qg,mype_qnr,mype_w
   public :: k_slmsk,k_tsea,k_vfrac,k_vtype,k_stype,k_zorl,k_smc,k_stc
-  public :: k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc,k_t2m,k_q2m,k_howv
+  public :: k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc,k_t2m,k_q2m,k_howv,k_gust
   public :: ijns,ijns2d,displss,displss2d,ijnz,displsz_g
   public :: fv3lam_io_dynmetvars3d_nouv,fv3lam_io_tracermetvars3d_nouv
   public :: fv3lam_io_tracerchemvars3d_nouv,fv3lam_io_tracersmokevars3d_nouv
@@ -158,7 +158,7 @@ module gsi_rfv3io_mod
   integer(i_kind) mype_qi,mype_qr,mype_qs,mype_qg,mype_qnr,mype_w
 
   integer(i_kind) k_slmsk,k_tsea,k_vfrac,k_vtype,k_stype,k_zorl,k_smc,k_stc
-  integer(i_kind) k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc,k_t2m,k_q2m,k_howv
+  integer(i_kind) k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc,k_t2m,k_q2m,k_howv,k_gust
   parameter(                   &  
     k_f10m =1,                  &   !fact10
     k_stype=2,                  &   !soil_type
@@ -174,7 +174,8 @@ module gsi_rfv3io_mod
     k_q2m  =12,                 & ! 2 m Q
     k_orog =13,                 & !terrain
     k_howv =14,                 &   ! significant wave height (aka howv in GSI)
-    n2d=14                   )
+    k_gust =15,                 &   ! wind gust (aka gust in GSI)
+    n2d=15                   )
   logical :: grid_reverse_flag
   character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_io_dynmetvars3d_nouv 
                                     ! copy of cvars3d excluding uv 3-d fields   
@@ -996,6 +997,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin)
     real(r_kind),dimension(:,:),pointer::ges_t2m=>NULL()
     real(r_kind),dimension(:,:),pointer::ges_q2m=>NULL()
     real(r_kind),dimension(:,:),pointer::ges_howv=>NULL()
+    real(r_kind),dimension(:,:),pointer::ges_gust=>NULL()
 
     real(r_kind),dimension(:,:,:),pointer::ges_ql=>NULL()
     real(r_kind),dimension(:,:,:),pointer::ges_qi=>NULL()
@@ -1274,6 +1276,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin)
           else if(trim(vartem)=='t2m') then
           else if(trim(vartem)=='q2m') then
           else if(trim(vartem)=='howv') then
+          else if(trim(vartem)=='gust') then
           else 
             write(6,*)'the metvarname2 ',trim(vartem),' has not been considered yet, stop'
             call stop2(333)
@@ -1294,7 +1297,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin)
         vartem=trim(name_metvars2d(i))
         if(.not.( (trim(vartem)=='ps'.and.fv3sar_bg_opt==0).or.(trim(vartem)=="z") &
                   .or.(trim(vartem)=="t2m").or.(trim(vartem)=="q2m")               &
-                  .or.(trim(vartem)=="howv")))  then ! z is treated separately
+                  .or.(trim(vartem)=="howv").or.(trim(vartem)=="gust")))  then ! z is treated separately
           if (ifindstrloc(vardynvars,trim(vartem)) > 0) then
             jdynvar=jdynvar+1
             fv3lam_io_dynmetvars2d_nouv(jdynvar)=trim(vartem)
@@ -1557,6 +1560,12 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin)
             if (ier/=0) call die(trim(myname),'cannot get pointers for howv, ier=',ier)
          endif
 
+!---     wind gust (gust)
+         if ( i_gust_3dda == 1 ) then
+            call GSI_BundleGetPointer(GSI_MetGuess_Bundle(it),'gust',ges_gust,istatus ); ier=ier+istatus
+            if (ier/=0) call die(trim(myname),'cannot get pointers for gust, ier=',ier)
+         endif
+
          if(mype == 0 ) then
            call check(nf90_open(fv3filenamegin(it)%dynvars,nf90_nowrite,loc_id))
            call check(nf90_inquire(loc_id,formatNum=ncfmt))
@@ -1739,7 +1748,8 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin)
          endif
 
 
-         call gsi_fv3ncdf2d_read(fv3filenamegin(it),it,ges_z,ges_t2m,ges_q2m,ges_howv)
+         call gsi_fv3ncdf2d_read(fv3filenamegin(it),it,ges_z,ges_t2m,ges_q2m,   &
+                                 ges_howv,ges_gust)
 
          if(i_use_2mq4b > 0 .and. i_use_2mt4b > 0 ) then
 ! Convert 2m guess mixing ratio to specific humidity
@@ -1975,7 +1985,8 @@ end subroutine gsi_bundlegetpointer_fv3lam_tracerchem_nouv
 
 end subroutine read_fv3_netcdf_guess
 
-subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv)
+subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,          &
+                              ges_howv,ges_gust)
 !$$$  subprogram documentation block
 !                .      .    .                                       .
 ! subprogram:    gsi_fv3ncdf2d_read       
@@ -2022,6 +2033,7 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv)
     real(r_kind),               intent(in),dimension(:,:),pointer::ges_t2m
     real(r_kind),               intent(in),dimension(:,:),pointer::ges_q2m
     real(r_kind),               intent(in),dimension(:,:),pointer::ges_howv
+    real(r_kind),               intent(in),dimension(:,:),pointer::ges_gust
     type (type_fv3regfilenameg),intent(in) :: fv3filenamegin
 
     character(len=max_varname_length) :: name
@@ -2036,8 +2048,8 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv)
     integer(i_kind) kk,n,ns,j,ii,jj,mm1
       character(len=:),allocatable :: sfcdata   !='fv3_sfcdata'
       character(len=:),allocatable :: dynvars   !='fv3_dynvars'
-! for checking the existence of howv in firstguess file
-    integer(i_kind) id_howv
+! for checking the existence of howv/gust in firstguess file
+    integer(i_kind) id_howv, id_gust
     integer(i_kind) iret_bcast
 
 ! for io_layout > 1
@@ -2057,8 +2069,9 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv)
     allocate(work(itotsub*n2d))
     allocate( sfcn2d(lat2,lon2,n2d))
 
-!-- initialisation of the array for howv
+!-- initialisation of the array for howv/gust
     sfcn2d(:,:,k_howv) = zero
+    sfcn2d(:,:,k_gust) = zero
 
 !-- initialisation of the array for sfc_var_exist 
     sfc_var_exist = .false.
@@ -2107,6 +2120,21 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv)
              trim(sfcdata), ', iret, varid = ',iret, id_howv,' (on pe: ', mype,').'
          end if
        end if
+!---   check the existence of wind gust (gust) in 2D FV3-LAM firstguess file
+!      (similar as done above for howv)
+       if ( i_gust_3dda == 1 ) then
+         iret = nf90_inq_varid(gfile_loc,'gust',id_gust)
+         if ( iret /= nf90_noerr ) then
+           iret = nf90_inq_varid(gfile_loc,'GUST',id_gust) ! double check with name in uppercase
+         end if
+         if ( iret /= nf90_noerr ) then
+           i_gust_3dda = 0                ! gust does not exist in firstguess, then stop GSI run.
+           call die('gsi_fv3ncdf2d_read','Warning: CANNOT find gust in firstguess, aborting..., iret = ', iret)
+         else
+           write(6,'(1x,A,1x,A,1x,A,1x,I4,1x,I4,1x,A,1x,I4.4,A)') 'gsi_fv3ncdf2d_read:: Found gust in firstguess ',  &
+             trim(sfcdata), ', iret, varid = ',iret, id_gust,' (on pe: ', mype,').'
+         end if
+       end if
 
    !!!!!!!!!!!! read in 2d variables !!!!!!!!!!!!!!!!!!!!!!!!!!
        do i=ndimensions+1,nvariables
@@ -2150,6 +2178,9 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv)
           else if( trim(name)=='HOWV'.or.trim(name)=='howv' ) then
              k=k_howv
              sfc_var_exist(k) = .true.
+          else if( trim(name)=='GUST'.or.trim(name)=='gust' ) then
+             k=k_gust
+             sfc_var_exist(k) = .true.
           else
              cycle 
           endif
@@ -2283,8 +2314,9 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv)
        if(allocated(sfc_fulldomain)) deallocate (sfc_fulldomain)
     endif  ! mype
 
-!-- broadcast the updated i_howv_3dda to all tasks (!!!!)
+!-- broadcast the updated i_howv_3dda, i_gust_3dda to all tasks (!!!!)
     call mpi_bcast(i_howv_3dda, 1, mpi_itype, mype_2d, mpi_comm_world, iret_bcast)
+    call mpi_bcast(i_gust_3dda, 1, mpi_itype, mype_2d, mpi_comm_world, iret_bcast)
 
 !-- broadcast the updated sfc_var_exist to all tasks (!!!!)
     call mpi_bcast(sfc_var_exist, n2d, mpi_itype, mype_2d, mpi_comm_world, iret_bcast)
@@ -2313,6 +2345,9 @@ subroutine gsi_fv3ncdf2d_read(fv3filenamegin,it,ges_z,ges_t2m,ges_q2m,ges_howv)
     if ( i_howv_3dda == 1 ) then
        if ( sfc_var_exist(k_howv)  ) ges_howv(:,:)=sfcn2d(:,:,k_howv)
     endif
+    if ( i_gust_3dda == 1 ) then
+       if ( sfc_var_exist(k_gust)  ) ges_gust(:,:)=sfcn2d(:,:,k_gust)
+    endif
     deallocate (sfcn2d,a)
     return
 end subroutine gsi_fv3ncdf2d_read
@@ -3628,6 +3663,7 @@ subroutine wrfv3_netcdf(fv3filenamegin)
     real(r_kind),pointer,dimension(:,:  ):: ges_t2m =>NULL()
     real(r_kind),pointer,dimension(:,:  ):: ges_q2m  =>NULL()
     real(r_kind),pointer,dimension(:,:  ):: ges_howv =>NULL()
+    real(r_kind),pointer,dimension(:,:  ):: ges_gust =>NULL()
 
     integer(i_kind) i,k
 
@@ -3750,6 +3786,9 @@ subroutine wrfv3_netcdf(fv3filenamegin)
     if ( i_howv_3dda == 1 ) then
        call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'howv',ges_howv,istatus); ier=ier+istatus
     endif
+    if ( i_gust_3dda == 1 ) then
+       call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'gust',ges_gust,istatus); ier=ier+istatus
+    endif
     if (ier/=0) call die('wrfv3_netcdf','cannot get pointers for fv3 met-fields, ier =',ier)
 
     if (laeroana_fv3cmaq) then
@@ -3964,6 +4003,10 @@ subroutine wrfv3_netcdf(fv3filenamegin)
     if ( i_howv_3dda == 1 ) then
       call gsi_fv3ncdf_write_sfc(fv3filenamegin,'howv',ges_howv,add_saved)
     endif
+!-- output analysis of gust
+    if ( i_gust_3dda == 1 ) then
+      call gsi_fv3ncdf_write_sfc(fv3filenamegin,'gust',ges_gust,add_saved)
+    endif
 
     if(allocated(g_prsi)) deallocate(g_prsi)
 

diff --git 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_howv, hwllp_howv, corp_gust, hwllp_gust, oerr_gust
   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
@@ -1604,6 +1604,9 @@ module gsimod
 !                           = 0.42 meters (default)
 !      hwllp_howv    - real, background error de-correlation length scale of howv 
 !                           = 170,000.0 meters (default 170 km)
+!      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
 !
   namelist/rapidrefresh_cldsurf/dfi_radar_latent_heat_time_period, &
                                 metar_impact_radius,metar_impact_radius_lowcloud, &
@@ -1625,7 +1628,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_howv, hwllp_howv, corp_gust, hwllp_gust, oerr_gust
 
 ! 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
@@ -313,6 +313,7 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt
       use mpeu_util,only: getindex
       use radiance_mod, only: icloud_cv,n_clouds_fwd,cloud_names_fwd
       use chemmod, only: berror_fv3_cmaq_regional,berror_fv3_sd_regional
+      use rapidrefresh_cldsurf_mod, only: corp_gust, hwllp_gust, l_rtma3d
 
       implicit none
 
@@ -825,11 +826,35 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt
 !       end if
      else if (n==nrf2_gust) then
         do i=1,mlat
-           corp(i,n)=three
+           corp(i,n)=three    ! background error stddev of wind gust = 3 m/s (default: legacy code from 2DRTMA) 
         end do
         do i=0,mlat+1
-           hwllp(i,n)=hwll(i,1,nrf3_q)
+           hwllp(i,n)=hwll(i,1,nrf3_q)  ! de-correlation length of bkgd error of gust is
+                                        ! same as the value of q at bottom level (default: legacy code from 2DRTMA)
+                                        ! for other DA apps, it is recommended to change it 
+                                        ! by setting hwllp_gust in GSI namelist.
         end do
+        if ( l_rtma3d ) then  ! For 3drtma only: allowing to change the stddev and 
+                              ! de-correlation length of bkgd error of gust:
+                              !   corp_gust : set in namelist(if <=0, using default value above (3.0)
+                              !   hwllp_gust: set in namelist(if <=0, using default value above (value of q)
+           if ( corp_gust  .gt. 0.0_r_kind ) then
+              corp(1:mlat,   n) = corp_gust
+              if (mype==0) write(6,'(1x,A,A,I5.5,A,F8.3)') &
+                           myname_,"@pe=",mype," (3drtma) set b_error stddev of gust = ",corp_gust
+           else
+              if (mype==0) write(6,'(1x,A,A,I5.5,A,F8.3)') &
+                           myname_,"@pe=",mype," (3drtma) set b_error stddev of gust (default) = ",three
+           end if
+           if ( hwllp_gust .gt. 0.0_r_kind ) then
+              hwllp(0:mlat+1,n) = hwllp_gust
+              if (mype==0) write(6,'(1x,A,A,I5.5,A,F12.3)') &
+                           myname_,"@pe=",mype," (3drtma) set b_error de-corr length of gust = ",hwllp_gust
+           else
+              if (mype==0) write(6,'(1x,A,A,I5.5,A)') &
+                           myname_,"@pe=",mype," (3drtma) set b_error de-corr length of gust is same as length of q."
+           end if
+        end if
      else if (n==nrf2_vis) then
         do i=1,mlat
            corp(i,n)=3.0_r_kind