From f282a9447e9d3c37e7866a7c9cebb96c0bb5068b Mon Sep 17 00:00:00 2001 From: xincjin-NOAA <122945910+xincjin-NOAA@users.noreply.github.com> Date: Tue, 12 Mar 2024 14:02:12 -0400 Subject: [PATCH 1/6] Assimilate GMI in GSI (#692) --- regression/global_4denvar.sh | 55 +++++----- regression/global_enkf.sh | 19 ++-- regression/regression_namelists.sh | 2 +- regression/regression_var.sh | 2 +- src/gsi/clw_mod.f90 | 2 +- src/gsi/deter_sfc_mod.f90 | 138 +++++++++----------------- src/gsi/radiance_mod.f90 | 6 +- src/gsi/read_gmi.f90 | 6 +- src/gsi/setuprad.f90 | 13 +-- src/gsi/ssmis_spatial_average_mod.f90 | 29 ++---- 10 files changed, 106 insertions(+), 166 deletions(-) diff --git a/regression/global_4denvar.sh b/regression/global_4denvar.sh index 08a62f5eb0..056815228b 100755 --- a/regression/global_4denvar.sh +++ b/regression/global_4denvar.sh @@ -55,14 +55,15 @@ cycg=`echo $gdate | cut -c9-10` dumpobs=gdas prefix_obs=${dumpobs}.t${cyca}z prefix_ges=gdas.t${cycg}z -prefix_ens=gdas.t${cycg}z +prefix_ens=enkfgdas.t${cycg}z suffix=tm00.bufr_d dumpges=gdas COMROOTgfs=$casesdir/gfs/prod -datobs=$COMROOTgfs/$dumpobs.$PDYa/${cyca}/atmos -datges=$COMROOTgfs/$dumpges.$PDYg/${cycg}/atmos -datens=$COMROOTgfs/enkfgdas.$PDYg/${cycg}/atmos +datobs=$COMROOTgfs/$dumpobs.$PDYa/${cyca}/obs +dathis=$COMROOTgfs/$dumpges.$PDYg/${cycg}/model_data/atmos/history +datanl=$COMROOTgfs/gdas.$PDYg/${cycg}/analysis/atmos +datens=$COMROOTgfs/enkfgdas.$PDYg/${cycg} # Set up $tmpdir @@ -265,28 +266,28 @@ $nln $datobs/${prefix_obs}.esatms.${suffix} ./atmsbufrears ## $nln $datobs/${prefix_obs}.amsr2.tm00.bufr_d ./amsr2bufr # Copy bias correction, atmospheric and surface files -$nln $datges/${prefix_ges}.abias ./satbias_in -$nln $datges/${prefix_ges}.abias_pc ./satbias_pc -$nln $datges/${prefix_ges}.abias_air ./aircftbias_in -$nln $datges/${prefix_ges}.radstat ./radstat.gdas - -$nln $datges/${prefix_ges}.sfcf003.nc ./sfcf03 -$nln $datges/${prefix_ges}.sfcf004.nc ./sfcf04 -$nln $datges/${prefix_ges}.sfcf005.nc ./sfcf05 -$nln $datges/${prefix_ges}.sfcf006.nc ./sfcf06 -$nln $datges/${prefix_ges}.sfcf007.nc ./sfcf07 -$nln $datges/${prefix_ges}.sfcf008.nc ./sfcf08 -$nln $datges/${prefix_ges}.sfcf009.nc ./sfcf09 - -$nln $datges/${prefix_ges}.atmf003.nc ./sigf03 -$nln $datges/${prefix_ges}.atmf004.nc ./sigf04 -$nln $datges/${prefix_ges}.atmf005.nc ./sigf05 -$nln $datges/${prefix_ges}.atmf006.nc ./sigf06 -$nln $datges/${prefix_ges}.atmf007.nc ./sigf07 -$nln $datges/${prefix_ges}.atmf008.nc ./sigf08 -$nln $datges/${prefix_ges}.atmf009.nc ./sigf09 - -$nln $datens/${prefix_ens}.sfcf006.ensmean.nc ./sfcf06_anlgrid +$nln $datanl/${prefix_ges}.abias ./satbias_in +$nln $datanl/${prefix_ges}.abias_pc ./satbias_pc +$nln $datanl/${prefix_ges}.abias_air ./aircftbias_in +$nln $datanl/${prefix_ges}.radstat ./radstat.gdas + +$nln $dathis/${prefix_ges}.sfcf003.nc ./sfcf03 +$nln $dathis/${prefix_ges}.sfcf004.nc ./sfcf04 +$nln $dathis/${prefix_ges}.sfcf005.nc ./sfcf05 +$nln $dathis/${prefix_ges}.sfcf006.nc ./sfcf06 +$nln $dathis/${prefix_ges}.sfcf007.nc ./sfcf07 +$nln $dathis/${prefix_ges}.sfcf008.nc ./sfcf08 +$nln $dathis/${prefix_ges}.sfcf009.nc ./sfcf09 + +$nln $dathis/${prefix_ges}.atmf003.nc ./sigf03 +$nln $dathis/${prefix_ges}.atmf004.nc ./sigf04 +$nln $dathis/${prefix_ges}.atmf005.nc ./sigf05 +$nln $dathis/${prefix_ges}.atmf006.nc ./sigf06 +$nln $dathis/${prefix_ges}.atmf007.nc ./sigf07 +$nln $dathis/${prefix_ges}.atmf008.nc ./sigf08 +$nln $dathis/${prefix_ges}.atmf009.nc ./sigf09 + +$nln $datens/ensstat/model_data/atmos/history/${prefix_ens}.sfcf006.ensmean.nc ./sfcf06_anlgrid export ENS_PATH='./ensemble_data/' mkdir -p ${ENS_PATH} @@ -296,7 +297,7 @@ for fh in $flist; do imem=1 while [[ $imem -le $NMEM_ENKF ]]; do member="mem"`printf %03i $imem` - $nln $datens/$member/$sigens ${ENS_PATH}sigf${fh}_ens_${member} + $nln $datens/$member/model_data/atmos/history/$sigens ${ENS_PATH}sigf${fh}_ens_${member} (( imem = $imem + 1 )) done done diff --git a/regression/global_enkf.sh b/regression/global_enkf.sh index e458c5830d..ca40abda52 100755 --- a/regression/global_enkf.sh +++ b/regression/global_enkf.sh @@ -51,17 +51,14 @@ cyca=`echo $global_adate | cut -c9-10` PDYg=`echo $gdate | cut -c1-8` cycg=`echo $gdate | cut -c9-10` -dumpobs=gdas -prefix_obs=${dumpobs}.t${cyca}z -prefix_ges=gdas.t${cycg}z -prefix_ens=gdas.t${cycg}z +prefix_obs=enkfgdas.t${cyca}z +prefix_ens=enkfgdas.t${cycg}z suffix=tm00.bufr_d dumpges=gdas COMROOTgfs=$casesdir/gfs/prod -datobs=$COMROOTgfs/enkfgdas.$PDYa/${cyca}/atmos -datges=$COMROOTgfs/$dumpges.$PDYg/${cycg}/atmos -datens=$COMROOTgfs/enkfgdas.$PDYg/${cycg}/atmos +datobs=$COMROOTgfs/enkfgdas.$PDYa/${cyca}/ensstat/analysis/atmos +datens=$COMROOTgfs/enkfgdas.$PDYg/${cycg} # Set up $tmpdir @@ -166,19 +163,19 @@ nfhrs=`echo $IAUFHRS_ENKF | sed 's/,/ /g'` for fhr in $nfhrs; do for imem in $(seq 1 $NMEM_ENKF); do memchar="mem"$(printf %03i $imem) - $nln $datens/$memchar/${prefix_ens}.atmf00${fhr}.nc sfg_${global_adate}_fhr0${fhr}_${memchar} + $nln $datens/$memchar/model_data/atmos/history/${prefix_ens}.atmf00${fhr}.nc sfg_${global_adate}_fhr0${fhr}_${memchar} if [ $cnvw_option = ".true." ]; then - $nln $datens/$memchar/${prefix_ens}sfcf00${fhr}.nc sfgsfc_${global_adate}_fhr0${fhr}_${memchar} + $nln $datens/$memchar/model_data/atmos/history/${prefix_ens}sfcf00${fhr}.nc sfgsfc_${global_adate}_fhr0${fhr}_${memchar} fi (( imem = $imem + 1 )) done - $nln $datens/${prefix_ens}.atmf00${fhr}.ensmean.nc sfg_${global_adate}_fhr0${fhr}_ensmean + $nln $datens/ensstat/model_data/atmos/history/${prefix_ens}.atmf00${fhr}.ensmean.nc sfg_${global_adate}_fhr0${fhr}_ensmean if [ $cnvw_option = ".true." ]; then $nln $datens/${prefix_ens}.sfcf00${fhr}.ensmean.nc sfgsfc_${global_adate}_fhr0${fhr}_ensmean fi done -$nln $datobs/${prefix_obs}.abias_int ./satbias_in +$nln $datobs/${prefix_obs}.abias_int.ensmean ./satbias_in cd $tmpdir diff --git a/regression/regression_namelists.sh b/regression/regression_namelists.sh index 552bc1ba59..7ca183ef3e 100755 --- a/regression/regression_namelists.sh +++ b/regression/regression_namelists.sh @@ -181,7 +181,7 @@ OBS_INPUT:: sstviirs viirs-m j1 viirs-m_j1 0.0 4 0 abibufr abi g18 abi_g18 0.0 1 0 ahibufr ahi himawari9 ahi_himawari9 0.0 1 0 - atmsbufr atms n21 atms_n21 0.0 1 1 + atmsbufr atms n21 atms_n21 0.0 1 0 crisfsbufr cris-fsr n21 cris-fsr_n21 0.0 1 0 sstviirs viirs-m j2 viirs-m_j2 0.0 4 0 ompsnpbufr ompsnp n21 ompsnp_n21 0.0 0 0 diff --git a/regression/regression_var.sh b/regression/regression_var.sh index 9d91b2c41f..315028675c 100755 --- a/regression/regression_var.sh +++ b/regression/regression_var.sh @@ -190,7 +190,7 @@ export savdir="$ptmp" export JCAP="62" # Case Study analysis dates -export global_adate="2022110900" +export global_adate="2024022300" export rtma_adate="2020022420" export fv3_netcdf_adate="2017030100" export rrfs_3denvar_glbens_adate="2021072518" diff --git a/src/gsi/clw_mod.f90 b/src/gsi/clw_mod.f90 index 512aaded01..49387f05eb 100644 --- a/src/gsi/clw_mod.f90 +++ b/src/gsi/clw_mod.f90 @@ -2019,7 +2019,7 @@ subroutine gmi_37pol_diff(tb37v,tb37h,tsim37v,tsim37h,clw,ierrret) clw = one - (tb37v-tb37h)/(tsim37v-tsim37h) clw=max(zero,clw) - if(tb37h > tb37v) then + if ((tb37h > tb37v) .or. (tb37h > 500_r_kind )) then ierrret = 1 clw= r1000 endif diff --git a/src/gsi/deter_sfc_mod.f90 b/src/gsi/deter_sfc_mod.f90 index 300d36cffb..271e81c5d2 100644 --- a/src/gsi/deter_sfc_mod.f90 +++ b/src/gsi/deter_sfc_mod.f90 @@ -33,7 +33,7 @@ module deter_sfc_mod use satthin, only: sno_full,isli_full,sst_full,soil_moi_full, & soil_temp_full,soil_type_full,veg_frac_full,veg_type_full, & fact10_full,zs_full,sfc_rough_full,zs_full_gfs - use constants, only: zero,one,two,one_tenth,deg2rad,rad2deg + use constants, only: zero,one,two,one_tenth,deg2rad,rad2deg, rearth use gridmod, only: nlat,nlon,regional,tll2xy,nlat_sfc,nlon_sfc,rlats_sfc,rlons_sfc, & rlats,rlons,dx_gfs,txy2ll,lpl_gfs use guess_grids, only: nfldsfc,hrdifsfc,ntguessfc @@ -1331,7 +1331,7 @@ subroutine deter_sfc_amsre_low(dlat_earth,dlon_earth,isflg,sfcpct) end subroutine deter_sfc_amsre_low -subroutine deter_sfc_gmi(dlat_earth,dlon_earth,isflg,sfcpct) +subroutine deter_sfc_gmi(dlat_earth,dlon_earth,isflg) !$$$ subprogram documentation block ! . . . . ! subprogram: deter_sfc_gmi determine land surface type @@ -1354,11 +1354,6 @@ subroutine deter_sfc_gmi(dlat_earth,dlon_earth,isflg,sfcpct) ! 2 sea ice ! 3 snow ! 4 mixed -! sfcpct(0:3)- percentage of 4 surface types -! (0) - sea percentage -! (1) - land percentage -! (2) - sea ice percentage -! (3) - snow percentage ! ! attributes: ! language: f90 @@ -1370,15 +1365,11 @@ subroutine deter_sfc_gmi(dlat_earth,dlon_earth,isflg,sfcpct) real(r_kind) ,intent(in ) :: dlat_earth,dlon_earth integer(i_kind) ,intent( out) :: isflg - real(r_kind),dimension(0:3),intent( out) :: sfcpct - - integer(i_kind) jsli,it - integer(i_kind):: klat1,klon1,klatp1,klonp1 - real(r_kind):: dx,dy,dx1,dy1,w00,w10,w01,w11 - real(r_kind) :: dlat,dlon + integer(i_kind) jsli,it, i, j + integer(i_kind):: klat1,klon1,klatp1,klonp1, ksmall, klarge, n_grid + real(r_kind) :: dlat,dlon, grid_dist + integer(i_kind):: klatn,klonn,klatpn,klonpn logical :: outside - integer(i_kind):: klat2,klon2,klatp2,klonp2 - ! ! For interpolation, we usually use o points (4points for land sea decision) ! In case of lowfreq channel (Large FOV), add the check of x points(8 points) @@ -1407,90 +1398,55 @@ subroutine deter_sfc_gmi(dlat_earth,dlon_earth,isflg,sfcpct) end if klon1=int(dlon); klat1=int(dlat) - dx =dlon-klon1; dy =dlat-klat1 - dx1 =one-dx; dy1 =one-dy - w00=dx1*dy1; w10=dx1*dy; w01=dx*dy1; w11=dx*dy klat1=min(max(1,klat1),nlat_sfc); klon1=min(max(0,klon1),nlon_sfc) if(klon1==0) klon1=nlon_sfc klatp1=min(nlat_sfc,klat1+1); klonp1=klon1+1 - if(klonp1==nlon_sfc+1) klonp1=1 - klonp2 = klonp1+1 - if(klonp2==nlon_sfc+1) klonp2=1 - klon2=klon1-1 - if(klon2==0)klon2=nlon_sfc - klat2=max(1,klat1-1) - klatp2=min(nlat_sfc,klatp1+1) ! Set surface type flag. Begin by assuming obs over ice-free water - sfcpct = zero - - jsli = isli_full(klat1 ,klon1 ) - if(sno_full(klat1 ,klon1 ,it) > one .and. jsli == 1)jsli=3 - sfcpct(jsli)=sfcpct(jsli)+one - - jsli = isli_full(klatp1,klon1 ) - if(sno_full(klatp1 ,klon1 ,it) > one .and. jsli == 1)jsli=3 - sfcpct(jsli)=sfcpct(jsli)+one - - jsli = isli_full(klat1 ,klonp1) - if(sno_full(klat1 ,klonp1 ,it) > one .and. jsli == 1)jsli=3 - sfcpct(jsli)=sfcpct(jsli)+one - - jsli = isli_full(klatp1,klonp1) - if(sno_full(klatp1 ,klonp1 ,it) > one .and. jsli == 1)jsli=3 - sfcpct(jsli)=sfcpct(jsli)+one - - jsli = isli_full(klatp2,klon1) - if(sno_full(klatp2 ,klon1 ,it) > one .and. jsli == 1)jsli=3 - sfcpct(jsli)=sfcpct(jsli)+one - - jsli = isli_full(klatp2,klonp1) - if(sno_full(klatp2 ,klonp1 ,it) > one .and. jsli == 1)jsli=3 - sfcpct(jsli)=sfcpct(jsli)+one - - jsli = isli_full(klatp1,klon2) - if(sno_full(klatp1 ,klon2 ,it) > one .and. jsli == 1)jsli=3 - sfcpct(jsli)=sfcpct(jsli)+one - - jsli = isli_full(klatp1,klonp2) - if(sno_full(klatp1 ,klonp2 ,it) > one .and. jsli == 1)jsli=3 - sfcpct(jsli)=sfcpct(jsli)+one - - jsli = isli_full(klat1,klon2) - if(sno_full(klat1 ,klon2 ,it) > one .and. jsli == 1)jsli=3 - sfcpct(jsli)=sfcpct(jsli)+one - - jsli = isli_full(klat1,klonp2) - if(sno_full(klat1 ,klonp2 ,it) > one .and. jsli == 1)jsli=3 - sfcpct(jsli)=sfcpct(jsli)+one - - jsli = isli_full(klat2,klon1) - if(sno_full(klat2 ,klon1 ,it) > one .and. jsli == 1)jsli=3 - sfcpct(jsli)=sfcpct(jsli)+one - - jsli = isli_full(klat2,klonp1) - if(sno_full(klat2 ,klonp1 ,it) > one .and. jsli == 1)jsli=3 - sfcpct(jsli)=sfcpct(jsli)+one - - sfcpct=sfcpct/12.0_r_kind - -! sfcpct(3)=min(sfcpct(3),sfcpct(1)) -! sfcpct(1)=max(zero,sfcpct(1)-sfcpct(3)) - - if(sfcpct(0) > 0.99_r_kind)then - isflg = 0 - else if(sfcpct(1) > 0.99_r_kind)then - isflg = 1 - else if(sfcpct(2) > 0.99_r_kind)then - isflg = 2 - else if(sfcpct(3) > 0.99_r_kind)then - isflg = 3 - else - isflg = 4 - end if + grid_dist=rearth * (rlats_sfc(klatp1) - rlats_sfc(klat1)) + n_grid=int(40000 / grid_dist) + 1 + klatn = max(klat1 - n_grid, 1) + klonn = klon1 - n_grid + if (klonn < 0) klonn = nlon_sfc - klonn + klatpn = min((klat1 + n_grid), nlat_sfc) + klonpn = klon1 + n_grid + if (klonpn > nlon_sfc) klonpn = klonpn - nlon_sfc + + isflg=0 + outer: do i = klatn, klatpn + ! assume n_grid > 2 + if (0 < klonpn - klonn .and. klonpn - klonn < nlon_sfc / 2) then + do j = klonn, klonpn + if (isli_full(i, j) /= 0) then + isflg = 1 + exit outer + end if + end do + else + if (klonpn < klonn) then + ksmall = klonpn + klarge = klonn + else + ksmall = klonn + klarge = klonpn + end if + do j = 1, ksmall + if (isli_full(i, j) /= 0) then + isflg = 1 + exit outer + endif + end do + do j = klarge, nlon_sfc + if (isli_full(i, j) /= 0) then + isflg = 1 + exit outer + end if + end do + end if + end do outer return end subroutine deter_sfc_gmi diff --git a/src/gsi/radiance_mod.f90 b/src/gsi/radiance_mod.f90 index 60aa0bc3cd..aae7794957 100644 --- a/src/gsi/radiance_mod.f90 +++ b/src/gsi/radiance_mod.f90 @@ -1326,11 +1326,7 @@ subroutine radiance_ex_biascor_gmi(radmod,clw_obs,clw_guess_retrieval,nchanl,cld do i=1,nchanl if (radmod%lcloud4crtm(i)<0) cycle - if (clw_obs <= cclr(i) .and. clw_guess_retrieval <= cclr(i) .and. abs(clw_obs-clw_guess_retrieval) < 0.001_r_kind) then - cld_rbc_idx(i)=one !clear/clear - else - cld_rbc_idx(i)=zero - endif + if ((clw_obs-cclr(i))*(clw_guess_retrieval-cclr(i))=0.005_r_kind) cld_rbc_idx(i)=zero end do return diff --git a/src/gsi/read_gmi.f90 b/src/gsi/read_gmi.f90 index f54263b129..6ad4d829a3 100644 --- a/src/gsi/read_gmi.f90 +++ b/src/gsi/read_gmi.f90 @@ -184,7 +184,7 @@ subroutine read_gmi(mype,val_gmi,ithin,rmesh,jsatid,gstime,& integer(i_kind) :: jc,bufsat,n integer(i_kind),dimension(5):: iobsdate integer(i_kind):: method,iobs,num_obs - integer(i_kind),parameter :: maxobs=4000000 + integer(i_kind),parameter :: maxobs=6000000 !-- integer(i_kind),parameter :: nscan=74 ! after binning ifov, 221/3 + 1 integer(i_kind),parameter :: nscan=221 @@ -414,7 +414,7 @@ subroutine read_gmi(mype,val_gmi,ithin,rmesh,jsatid,gstime,& call ufbrep(lnbufr,var_check1,1,nchanl,iret,'GMICHQ') !call ufbrep(lnbufr,gmirfi,1,nchanl,iret,'GMIRFI') call ufbrep(lnbufr,pixelsaza,1,ngs,iret,'SAZA') - call ufbrep(lnbufr,val_angls,n_angls,ngs,iret,'SAMA SZA SMA SGA') + call ufbrep(lnbufr,val_angls,n_angls,ngs,iret,'BEARAZ SOZA SOLAZI SSGA') call ufbint(lnbufr,pixelloc,2, 1,iret,'CLATH CLONH') if (any(var_check1 < 99999999999_r_double)) then ! 100000000000 seems to be the missing value @@ -696,7 +696,7 @@ subroutine read_gmi(mype,val_gmi,ithin,rmesh,jsatid,gstime,& call deter_sfc(dlat,dlon,dlat_earth,dlon_earth,t4dv,isflg,idomsfc,sfcpct, & ts,tsavg,vty,vfr,sty,stp,sm,sn,zz,ff10,sfcr) - call deter_sfc_gmi(dlat_earth,dlon_earth,isflg,sfcpct) + call deter_sfc_gmi(dlat_earth,dlon_earth,isflg) ! Only keep obs over ocean - ej diff --git a/src/gsi/setuprad.f90 b/src/gsi/setuprad.f90 index ad8ccbc4e6..136568d1f3 100644 --- a/src/gsi/setuprad.f90 +++ b/src/gsi/setuprad.f90 @@ -1099,6 +1099,7 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& endif predbias=zero + !$omp parallel do schedule(dynamic,1) private(i,mm,j,k,tlap,node,bias) do i=1,nchanl mm=ich(i) @@ -1177,7 +1178,6 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& end do end if - do j = 1,npred predbias(j,i) = predchan(j,i)*pred(j,i) end do @@ -1263,8 +1263,8 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& if(amsua.or.atms) then call ret_amsua(tsim_bc,nchanl,tsavg5,zasat,clw_guess_retrieval,ierrret) else if(gmi) then - call gmi_37pol_diff(tsim_bc(6),tsim_bc(7),tsim_clr_bc(6),tsim_clr_bc(7),clw_guess_retrieval,ierrret) - call gmi_37pol_diff(tb_obs(6),tb_obs(7),tsim_clr_bc(6),tsim_clr_bc(7),clw_obs,ierrret) + call gmi_37pol_diff(tsim(6),tsim(7),tsim_clr(6),tsim_clr(7),clw_guess_retrieval,ierrret) + call gmi_37pol_diff(tb_obs(6),tb_obs(7),tsim_clr(6),tsim_clr(7),clw_obs,ierrret) end if if (radmod%ex_obserr=='ex_obserr1') then call radiance_ex_biascor(radmod,nchanl,tsim_bc,tsavg5,zasat, & @@ -1298,11 +1298,12 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& do i=1,nchanl pred(6,i) = zero pred(7,i) = zero - clw_avg = half*(clw_obs+clw_guess_retrieval) +! Need to investigate clw_ave = half*(clw_obs+clw_guess_retrieval) + clw_avg = zero if (i > 3 .and. clw_obs > 0.05_r_kind .and. clw_guess_retrieval > 0.05_r_kind .and. & - abs(clw_obs-clw_guess_retrieval) < 0.005_r_kind .and. clw_avg < 0.5_r_kind) cld_rbc_idx2(i) = one + abs(clw_obs-clw_guess_retrieval) < 0.005_r_kind .and. clw_avg < 0.5_r_kind) cld_rbc_idx2(i) = zero if (i < 5 .and. clw_obs > 0.2_r_kind .and. clw_guess_retrieval > 0.2_r_kind .and. & - abs(clw_obs-clw_guess_retrieval) < 0.005_r_kind .and. clw_avg < 0.5_r_kind) cld_rbc_idx2(i) = one + abs(clw_obs-clw_guess_retrieval) < 0.005_r_kind .and. clw_avg < 0.5_r_kind) cld_rbc_idx2(i) = zero if( i > 3 .and. clw_obs > 0.05_r_kind .and. clw_guess_retrieval > 0.05_r_kind .and. cld_rbc_idx(i) == zero) then pred(6,i) = clw_avg*clw_avg diff --git a/src/gsi/ssmis_spatial_average_mod.f90 b/src/gsi/ssmis_spatial_average_mod.f90 index 64dd5c6cf7..30b69c6223 100644 --- a/src/gsi/ssmis_spatial_average_mod.f90 +++ b/src/gsi/ssmis_spatial_average_mod.f90 @@ -682,26 +682,15 @@ SUBROUTINE SSMIS_Spatial_Average(BufrSat, Method, Num_Obs, NChanl, & ! Define grid box by channel - ! Ch 1-2: 1 scan direction, 1 track direction ! Ch 3-13: 3 scan direction, 3 track direction - if ((ic == 1) .or. (ic == 2)) then - ns1 = iscan - ns2 = iscan - if (ns1 < 1) ns1=1 - if (ns2 > max_scan) ns2=max_scan - np1 = ifov - np2 = ifov - if (np1 < 1) np1=1 - if (np2 > max_fov_gmi) np2=max_fov_gmi - else if ((ic > 2) .and. (ic < 14)) then - ns1 = iscan-1 - ns2 = iscan+1 - if (ns1 < 1) ns1=1 - if (ns2 > max_scan) ns2=max_scan - np1 = ifov-1 - np2 = ifov+1 - if (np1 < 1) np1=1 - if (np2 > max_fov_gmi) np2=max_fov_gmi - endif + ns1 = iscan-4 + ns2 = iscan+4 + if (ns1 < 1) ns1=1 + if (ns2 > max_scan) ns2=max_scan + np1 = ifov-8 + np2 = ifov+8 + if (np1 < 1) np1=1 + if (np2 > max_fov_gmi) np2=max_fov_gmi xnum = 0.0_r_kind mta = 0.0_r_kind if (any(bt_image_orig(np1:np2,ns1:ns2,ic) < btmin .or. & @@ -716,7 +705,7 @@ SUBROUTINE SSMIS_Spatial_Average(BufrSat, Method, Num_Obs, NChanl, & lat2 = latitude(ip,is) lon2 = longitude(ip,is) dist = distance(lat1,lon1,lat2,lon2) - if (dist > 50.0_r_kind) cycle gmi_box_x1 ! outside the box + if (dist > 20.0_r_kind) cycle gmi_box_x1 ! outside the box if (gaussian_wgt) then wgt = exp(-0.5_r_kind*(dist/sigma)*(dist/sigma)) else From a8d670c9f7096b2c7ac81801d6f2146b6c08bf4c Mon Sep 17 00:00:00 2001 From: RussTreadon-NOAA <26926959+RussTreadon-NOAA@users.noreply.github.com> Date: Thu, 14 Mar 2024 11:55:05 -0400 Subject: [PATCH 2/6] Update fix submodule (Rcov_crisn21) & GSI_BINARY_SOURCE_DIR (#718) --- fix | 2 +- modulefiles/gsi_gaea.intel.lua | 2 +- modulefiles/gsi_hera.gnu.lua | 2 +- modulefiles/gsi_hera.intel.lua | 2 +- modulefiles/gsi_hercules.intel.lua | 2 +- modulefiles/gsi_jet.intel.lua | 2 +- modulefiles/gsi_orion.intel.lua | 2 +- modulefiles/gsi_s4.intel.lua | 2 +- modulefiles/gsi_wcoss2.intel.lua | 2 +- 9 files changed, 9 insertions(+), 9 deletions(-) diff --git a/fix b/fix index 298bdc0f56..92de100c4d 160000 --- a/fix +++ b/fix @@ -1 +1 @@ -Subproject commit 298bdc0f56692ad25141a946177e5d88884e367a +Subproject commit 92de100c4d5893e9d6409afbdda6937b0de1cb3b diff --git a/modulefiles/gsi_gaea.intel.lua b/modulefiles/gsi_gaea.intel.lua index ef6b9ddba7..96643202a7 100644 --- a/modulefiles/gsi_gaea.intel.lua +++ b/modulefiles/gsi_gaea.intel.lua @@ -31,7 +31,7 @@ local MKLROOT="/opt/intel/oneapi/mkl/2022.0.2/" prepend_path("LD_LIBRARY_PATH",pathJoin(MKLROOT,"lib/intel64")) pushenv("MKLROOT", MKLROOT) -pushenv("GSI_BINARY_SOURCE_DIR", "/lustre/f2/dev/role.epic/contrib/GSI_data/fix/20230911") +pushenv("GSI_BINARY_SOURCE_DIR", "/lustre/f2/dev/role.epic/contrib/GSI_data/fix/20240208") setenv("CC","cc") setenv("FC","ftn") diff --git a/modulefiles/gsi_hera.gnu.lua b/modulefiles/gsi_hera.gnu.lua index d85e79b005..eab352553f 100644 --- a/modulefiles/gsi_hera.gnu.lua +++ b/modulefiles/gsi_hera.gnu.lua @@ -22,6 +22,6 @@ load("gsi_common") load(pathJoin("prod_util", prod_util_ver)) load(pathJoin("openblas", openblas_ver)) -pushenv("GSI_BINARY_SOURCE_DIR", "/scratch1/NCEPDEV/global/glopara/fix/gsi/20230911") +pushenv("GSI_BINARY_SOURCE_DIR", "/scratch1/NCEPDEV/global/glopara/fix/gsi/20240208") whatis("Description: GSI environment on Hera with GNU Compilers") diff --git a/modulefiles/gsi_hera.intel.lua b/modulefiles/gsi_hera.intel.lua index b967c32f74..2d36f375e0 100644 --- a/modulefiles/gsi_hera.intel.lua +++ b/modulefiles/gsi_hera.intel.lua @@ -20,6 +20,6 @@ load(pathJoin("prod_util", prod_util_ver)) pushenv("CFLAGS", "-xHOST") pushenv("FFLAGS", "-xHOST") -pushenv("GSI_BINARY_SOURCE_DIR", "/scratch1/NCEPDEV/global/glopara/fix/gsi/20230911") +pushenv("GSI_BINARY_SOURCE_DIR", "/scratch1/NCEPDEV/global/glopara/fix/gsi/20240208") whatis("Description: GSI environment on Hera with Intel Compilers") diff --git a/modulefiles/gsi_hercules.intel.lua b/modulefiles/gsi_hercules.intel.lua index c72323cbad..66ec9b03e1 100644 --- a/modulefiles/gsi_hercules.intel.lua +++ b/modulefiles/gsi_hercules.intel.lua @@ -21,6 +21,6 @@ load("intel-oneapi-mkl/2022.2.1") pushenv("CFLAGS", "-xHOST") pushenv("FFLAGS", "-xHOST") -pushenv("GSI_BINARY_SOURCE_DIR", "/work/noaa/global/glopara/fix/gsi/20230911") +pushenv("GSI_BINARY_SOURCE_DIR", "/work/noaa/global/glopara/fix/gsi/20240208") whatis("Description: GSI environment on Hercules with Intel Compilers") diff --git a/modulefiles/gsi_jet.intel.lua b/modulefiles/gsi_jet.intel.lua index b2b98be7ff..cc43e98260 100644 --- a/modulefiles/gsi_jet.intel.lua +++ b/modulefiles/gsi_jet.intel.lua @@ -20,6 +20,6 @@ load(pathJoin("prod_util", prod_util_ver)) pushenv("CFLAGS", "-axSSE4.2,AVX,CORE-AVX2") pushenv("FFLAGS", "-axSSE4.2,AVX,CORE-AVX2") -pushenv("GSI_BINARY_SOURCE_DIR", "/mnt/lfs4/HFIP/hfv3gfs/glopara/git/fv3gfs/fix/gsi/20230911") +pushenv("GSI_BINARY_SOURCE_DIR", "/mnt/lfs4/HFIP/hfv3gfs/glopara/git/fv3gfs/fix/gsi/20240208") whatis("Description: GSI environment on Jet with Intel Compilers") diff --git a/modulefiles/gsi_orion.intel.lua b/modulefiles/gsi_orion.intel.lua index 114e5f8ad1..03ea22018d 100644 --- a/modulefiles/gsi_orion.intel.lua +++ b/modulefiles/gsi_orion.intel.lua @@ -20,6 +20,6 @@ load(pathJoin("prod_util", prod_util_ver)) pushenv("CFLAGS", "-xHOST") pushenv("FFLAGS", "-xHOST") -pushenv("GSI_BINARY_SOURCE_DIR", "/work/noaa/global/glopara/fix/gsi/20230911") +pushenv("GSI_BINARY_SOURCE_DIR", "/work/noaa/global/glopara/fix/gsi/20240208") whatis("Description: GSI environment on Orion with Intel Compilers") diff --git a/modulefiles/gsi_s4.intel.lua b/modulefiles/gsi_s4.intel.lua index 08e3d38516..04945eef3e 100644 --- a/modulefiles/gsi_s4.intel.lua +++ b/modulefiles/gsi_s4.intel.lua @@ -20,6 +20,6 @@ load(pathJoin("prod_util", prod_util_ver)) pushenv("CFLAGS", "-march=ivybridge") pushenv("FFLAGS", "-march=ivybridge") -pushenv("GSI_BINARY_SOURCE_DIR", "/data/prod/glopara/fix/gsi/20230911") +pushenv("GSI_BINARY_SOURCE_DIR", "/data/prod/glopara/fix/gsi/20240208") whatis("Description: GSI environment on S4 with Intel Compilers") diff --git a/modulefiles/gsi_wcoss2.intel.lua b/modulefiles/gsi_wcoss2.intel.lua index 28b2f0d09d..c3bfd1156c 100644 --- a/modulefiles/gsi_wcoss2.intel.lua +++ b/modulefiles/gsi_wcoss2.intel.lua @@ -46,6 +46,6 @@ load(pathJoin("ncio", ncio_ver)) load(pathJoin("crtm", crtm_ver)) load(pathJoin("ncdiag",ncdiag_ver)) -pushenv("GSI_BINARY_SOURCE_DIR", "/lfs/h2/emc/global/noscrub/emc.global/FIX/fix/gsi/20230911") +pushenv("GSI_BINARY_SOURCE_DIR", "/lfs/h2/emc/global/noscrub/emc.global/FIX/fix/gsi/20240208") whatis("Description: GSI environment on WCOSS2") From 8d740a764d1b9d32e11cf1d0b4dd0ca26873871f Mon Sep 17 00:00:00 2001 From: RussTreadon-NOAA <26926959+RussTreadon-NOAA@users.noreply.github.com> Date: Mon, 18 Mar 2024 15:26:36 -0400 Subject: [PATCH 3/6] Update read_ozone.f90 to handle GOME data before and after 20240131 18Z (#721) --- src/gsi/read_ozone.f90 | 61 +++++++++++++++++++++++++++++------------- 1 file changed, 42 insertions(+), 19 deletions(-) diff --git a/src/gsi/read_ozone.f90 b/src/gsi/read_ozone.f90 index 2ad95a858e..e6d3411d97 100644 --- a/src/gsi/read_ozone.f90 +++ b/src/gsi/read_ozone.f90 @@ -138,7 +138,7 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & character(8) subset,subset6,subset8,subset8_ompsnp character(49) ozstr,ozostr character(63) lozstr - character(51) ozgstr + character(51) ozgstr_v1,ozgstr_v2 character(27) ozgstr2 character(42) ozostr2 character(64) mlstr @@ -165,11 +165,12 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & ! maximum number of observations set to real(r_kind),allocatable,dimension(:,:):: ozout - real(r_double) toq,poq + real(r_double) toq,poq,orbn real(r_double),dimension(nloz_v6):: ozone_v6 real(r_double),dimension(29,nloz_v8):: ozone_v8 real(r_double),dimension(10):: hdroz - real(r_double),dimension(10):: hdrozg + integer(i_kind):: nhdrozg + real(r_double),allocatable,dimension(:):: hdrozg real(r_double),dimension(5):: hdrozg2 real(r_double),dimension(10):: hdrozo real(r_double),dimension(8) :: hdrozo2 @@ -195,8 +196,10 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & data lozstr & / 'OSP12 OSP11 OSP10 OSP9 OSP8 OSP7 OSP6 OSP5 OSP4 OSP3 OSP2 OSP1 ' / - data ozgstr & - / 'SAID CLAT CLON YEAR DOYR HOUR MINU SECO SOZA SOLAZI' / + data ozgstr_v1 & + / 'SAID CLAT CLON SOZA SOLAZI YEAR DOYR HOUR MINU SECO' / + data ozgstr_v2 & + / 'SAID CLAT CLON SOZA SOLAZI YEAR MNTH DAYS HOUR MINU SECO' / data ozgstr2 & / 'CLDMNT SNOC ACIDX STKO FOVN' / data ozostr & @@ -482,8 +485,19 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & cycle obsloop endif -! extract header information - call ufbint(lunin,hdrozg,10,1,iret,ozgstr) +! Test for BUFR version using ORBN mnemonic + call ufbint(lunin,orbn,1,1,iret,'ORBN') + if (orbn > 100000000.0_r_kind) then + nhdrozg = 11 + else + nhdrozg = 10 + endif + if (.not.allocated(hdrozg)) allocate(hdrozg(nhdrozg)) + if (nhdrozg == 11) then + call ufbint(lunin,hdrozg,nhdrozg,1,iret,ozgstr_v2) + else + call ufbint(lunin,hdrozg,nhdrozg,1,iret,ozgstr_v1) + endif call ufbint(lunin,hdrozg2,5,1,iret,ozgstr2) rsat = hdrozg(1); ksatid=rsat @@ -494,7 +508,7 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & if (ksatid /= kidsat) cycle obsloop ! NESDIS does not put a flag for high SZA gome-2 data (SZA > 84 degree) - if ( hdrozg(9) > r84 ) cycle obsloop + if ( hdrozg(4) > r84 ) cycle obsloop nmrecs=nmrecs+nloz+1 @@ -520,15 +534,24 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & endif ! Convert observation time to relative time - idate5(1) = hdrozg(4) !year - IDAYYR = hdrozg(5) ! Day of year - JULIAN = -31739 + 1461 * (idate5(1) + 4799) /4 & - -3 * ((idate5(1) + 4899) / 100) / 4 + IDAYYR - call w3fs26(JULIAN,idate5(1),idate5(2),idate5(3),IDAYWK,IDAYYR) -! idate5(2) month -! idate5(3) day - idate5(4) = hdrozg(6) !hour - idate5(5) = hdrozg(7) !minute + if (nhdrozg == 11) then + idate5(1) = hdrozg(6) !year + idate5(2) = hdrozg(7) !month + idate5(3) = hdrozg(8) !day + idate5(4) = hdrozg(9) !hour + idate5(5) = hdrozg(10) !minute + else + idate5(1) = hdrozg(6) !year + IDAYYR = hdrozg(7) ! Day of year + JULIAN = -31739 + 1461 * (idate5(1) + 4799) /4 & + -3 * ((idate5(1) + 4899) / 100) / 4 + IDAYYR + call w3fs26(JULIAN,idate5(1),idate5(2),idate5(3),IDAYWK,IDAYYR) +! idate5(2) month +! idate5(3) day + idate5(4) = hdrozg(8) !hour + idate5(5) = hdrozg(9) !minute + endif + call w3fs21(idate5,nmind) t4dv=real((nmind-iwinbgn),r_kind)*r60inv sstime=real(nmind,r_kind) @@ -574,8 +597,8 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, & ozout(5,itx)=dlon_earth_deg ! earth relative longitude (degrees) ozout(6,itx)=dlat_earth_deg ! earth relative latitude (degrees) ozout(7,itx)=toq ! total ozone error flag - ozout(8,itx)=hdrozg(9) ! solar zenith angle - ozout(9,itx)=hdrozg(10) ! solar azimuth angle + ozout(8,itx)=hdrozg(4) ! solar zenith angle + ozout(9,itx)=hdrozg(5) ! solar azimuth angle ozout(10,itx)=hdrozg2(1) ! CLOUD AMOUNT IN SEGMENT ozout(11,itx)=hdrozg2(2) ! SNOW COVER ozout(12,itx)=hdrozg2(3) ! AEROSOL CONTAMINATION INDEX From dfb958fa9372c10c808d20e64f2955652017ee32 Mon Sep 17 00:00:00 2001 From: RussTreadon-NOAA <26926959+RussTreadon-NOAA@users.noreply.github.com> Date: Mon, 18 Mar 2024 15:27:23 -0400 Subject: [PATCH 4/6] Update Hera intel modulefile to Rocky 8 (#715) --- modulefiles/gsi_hera.intel.lua | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modulefiles/gsi_hera.intel.lua b/modulefiles/gsi_hera.intel.lua index 2d36f375e0..d21b9195c3 100644 --- a/modulefiles/gsi_hera.intel.lua +++ b/modulefiles/gsi_hera.intel.lua @@ -1,7 +1,7 @@ help([[ ]]) -prepend_path("MODULEPATH", "/scratch1/NCEPDEV/nems/role.epic/spack-stack/spack-stack-1.6.0/envs/gsi-addon-dev/install/modulefiles/Core") +prepend_path("MODULEPATH", "/scratch1/NCEPDEV/nems/role.epic/spack-stack/spack-stack-1.6.0/envs/gsi-addon-dev-rocky8/install/modulefiles/Core") local python_ver=os.getenv("python_ver") or "3.11.6" local stack_intel_ver=os.getenv("stack_intel_ver") or "2021.5.0" From 4e8107c0ae054bc9dcc7b9a2684479d97a1e4261 Mon Sep 17 00:00:00 2001 From: "Tseganeh Z. Gichamo" Date: Tue, 19 Mar 2024 12:44:32 -0600 Subject: [PATCH 5/6] Add parallel netcdf read/write from EnKF for sfc files (paranc option) (#709) --- src/enkf/controlvec.f90 | 16 +- src/enkf/gridio_gfs.f90 | 453 +++++++++++++++++++++++++++++++++++++++- 2 files changed, 452 insertions(+), 17 deletions(-) diff --git a/src/enkf/controlvec.f90 b/src/enkf/controlvec.f90 index 808eae2e28..0961549634 100644 --- a/src/enkf/controlvec.f90 +++ b/src/enkf/controlvec.f90 @@ -191,7 +191,7 @@ subroutine read_control() ! read ensemble members on IO tasks implicit none real(r_double) :: t1,t2 -integer(i_kind) :: nb,ne +integer(i_kind) :: nb,nlev,ne integer(i_kind) :: q_ind integer(i_kind) :: ierr @@ -218,19 +218,23 @@ subroutine read_control() if (nproc == 0) t1 = mpi_wtime() call readgriddata_pnc(cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,nbackgrounds, & fgfileprefixes,fgsfcfileprefixes,reducedgrid,grdin,qsat) + if (nproc == 0) then + t2 = mpi_wtime() + print *,'time in readgrid_pnc on root',t2-t1,'secs' + end if end if if (nproc <= ntasks_io-1) then if (.not. paranc) then if (nproc == 0) t1 = mpi_wtime() call readgriddata(nanal1(nproc),nanal2(nproc),cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,nbackgrounds, & fgfileprefixes,fgsfcfileprefixes,reducedgrid,grdin,qsat) + if (nproc == 0) then + t2 = mpi_wtime() + print *,'time in readgrid on root',t2-t1,'secs' + end if end if !print *,'min/max qsat',nanal,'=',minval(qsat),maxval(qsat) q_ind = getindex(cvars3d, 'q') - if (nproc == 0) then - t2 = mpi_wtime() - print *,'time in readgridata on root',t2-t1,'secs' - end if if (pseudo_rh .and. q_ind > 0) then do ne=1,nanals_per_iotask do nb=1,nbackgrounds @@ -357,7 +361,7 @@ subroutine write_control(no_inflate_flag) endif deallocate(grdin_mean) t2 = mpi_wtime() - print *,'time in write_control on root',t2-t1,'secs' + print *,'time in write_control paranc on root',t2-t1,'secs' endif end if diff --git a/src/enkf/gridio_gfs.f90 b/src/enkf/gridio_gfs.f90 index 8afb012fcf..35a0c3fbe4 100644 --- a/src/enkf/gridio_gfs.f90 +++ b/src/enkf/gridio_gfs.f90 @@ -89,15 +89,21 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & real(r_kind), dimension(ndimspec) :: vrtspec,divspec real(r_kind), allocatable, dimension(:) :: psg,pstend,ak,bk real(r_single),allocatable,dimension(:,:,:) :: ug3d,vg3d - type(Dataset) :: dset + type(Dataset) :: dset, dset_sfc type(Dimension) :: londim,latdim,levdim integer(i_kind) :: u_ind, v_ind, tv_ind, q_ind, oz_ind, cw_ind integer(i_kind) :: qr_ind, qs_ind, qg_ind integer(i_kind) :: tsen_ind, ql_ind, qi_ind, prse_ind integer(i_kind) :: ps_ind, pst_ind, sst_ind + ! surface + integer(i_kind) :: tmp2m_ind, spfh2m_ind, soilt1_ind, soilt2_ind, soilt3_ind + integer(i_kind) :: soilt4_ind,slc1_ind, slc2_ind, slc3_ind, slc4_ind integer(i_kind) :: k,iret,nb,i,imem,idvc,nlonsin,nlatsin,nlevsin,ne,nanal + ! surface + integer(i_kind) :: nlonsin_sfc,nlatsin_sfc + logical ice logical use_full_hydro integer(i_kind), allocatable, dimension(:) :: mem_pe, lev_pe1, lev_pe2, iocomms @@ -111,12 +117,6 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & call set_ncio_file_flags(vars3d, n3d, vars2d, n2d, read_sfc_file, read_atm_file) - if (read_sfc_file) then - print *,'paranc not supported for reading surface files' - call mpi_barrier(mpi_comm_world,ierr) - call mpi_finalize(ierr) - endif - ! figure out what member to read and do MPI sub-communicator things allocate(mem_pe(0:numproc-1)) allocate(iocomms(nanals)) @@ -152,6 +152,7 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & displs(i+1) = ((lev_pe1(i)-1)*nlons*nlats) end do + if (read_atm_file) then ! loop through times and do the read ne = 1 @@ -159,7 +160,6 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & write(charnanal,'(a3, i3.3)') 'mem', nanal filename = trim(adjustl(datapath))//trim(adjustl(fileprefixes(nb)))//trim(charnanal) - sfcfilename = trim(adjustl(datapath))//trim(adjustl(filesfcprefixes(nb)))//trim(charnanal) if (use_gfs_ncio) then dset = open_dataset(filename, paropen=.true., mpicomm=iocomms(mem_pe(nproc))) londim = get_dim(dset,'grid_xt'); nlonsin = londim%len @@ -496,6 +496,141 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & end do backgroundloop ! loop over backgrounds to read in + end if !read_atm_file + + if (read_sfc_file) then + ! loop through times and do the read + ne = 1 + sfcbackgroundloop: do nb=1,ntimes + + write(charnanal,'(a3, i3.3)') 'mem', nanal + sfcfilename = trim(adjustl(datapath))//trim(adjustl(filesfcprefixes(nb)))//trim(charnanal) + if (use_gfs_ncio) then + dset_sfc = open_dataset(sfcfilename, paropen=.true., mpicomm=iocomms(mem_pe(nproc))) + else + write(6,*)'READGRIDDATA_PNC sfc: ***FATAL ERROR*** parallel read only supported for netCDF' , ' PROGRAM STOPS' + call mpi_barrier(mpi_comm_world,ierr) + call mpi_finalize(ierr) + end if + if ( reducedgrid ) then + write(6,*) "READGRIDDATA_PNC sfc: reducedgrid=T interpolation not valid for writing sfc files" + call mpi_barrier(mpi_comm_world,ierr) + call mpi_finalize(ierr) + endif + + ! land sfc DA variables + tmp2m_ind = getindex(vars2d, 't2m') + spfh2m_ind = getindex(vars2d, 'q2m') + soilt1_ind = getindex(vars2d, 'st1') + slc1_ind = getindex(vars2d, 'sl1') + soilt2_ind = getindex(vars2d, 'st2') + slc2_ind = getindex(vars2d, 'sl2') + soilt3_ind = getindex(vars2d, 'st3') + slc3_ind = getindex(vars2d, 'sl3') + soilt4_ind = getindex(vars2d, 'st4') + slc4_ind = getindex(vars2d, 'sl4') + + ! read in sfc vars, if requested + if (tmp2m_ind > 0) then + call read_vardata(dset_sfc, 'tmp2m', values_2d, errcode=iret) + if (iret /= 0) then + print *,'READGRIDDATA_PNC: error reading tmp2m' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + if (iope==0) call copytogrdin(ug,grdin(:,levels(n3d) + tmp2m_ind,nb,ne)) + endif + if (spfh2m_ind > 0) then + call read_vardata(dset_sfc, 'spfh2m', values_2d, errcode=iret) + if (iret /= 0) then + print *,'READGRIDDATA_PNC: error reading spfh2m' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + if (iope==0) call copytogrdin(ug,grdin(:,levels(n3d) + spfh2m_ind,nb,ne)) + endif + if (soilt1_ind > 0) then + call read_vardata(dset_sfc, 'soilt1', values_2d, errcode=iret) + if (iret /= 0) then + print *,'READGRIDDATA_PNC: error reading soilt1' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + if (iope==0) call copytogrdin(ug,grdin(:,levels(n3d) + soilt1_ind,nb,ne)) + endif + if (soilt2_ind > 0) then + call read_vardata(dset_sfc, 'soilt2', values_2d, errcode=iret) + if (iret /= 0) then + print *,'READGRIDDATA_PNC: error reading soilt2' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + if (iope==0) call copytogrdin(ug,grdin(:,levels(n3d) + soilt2_ind,nb,ne)) + endif + if (soilt3_ind > 0) then + call read_vardata(dset_sfc, 'soilt3', values_2d, errcode=iret) + if (iret /= 0) then + print *,'READGRIDDATA_PNC: error reading soilt3' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + if (iope==0) call copytogrdin(ug,grdin(:,levels(n3d) + soilt3_ind,nb,ne)) + endif + if (soilt4_ind > 0) then + call read_vardata(dset_sfc, 'soilt4', values_2d, errcode=iret) + if (iret /= 0) then + print *,'READGRIDDATA_PNC: error reading soilt2' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + if (iope==0) call copytogrdin(ug,grdin(:,levels(n3d) + soilt4_ind,nb,ne)) + endif + if (slc1_ind > 0) then + call read_vardata(dset_sfc, 'soill1', values_2d, errcode=iret) + if (iret /= 0) then + print *,'READGRIDDATA_PNC: error reading soill1' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + if (iope==0) call copytogrdin(ug,grdin(:,levels(n3d) + slc1_ind,nb,ne)) + endif + if (slc2_ind > 0) then + call read_vardata(dset_sfc, 'soill2', values_2d, errcode=iret) + if (iret /= 0) then + print *,'READGRIDDATA_PNC: error reading soill2' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + if (iope==0) call copytogrdin(ug,grdin(:,levels(n3d) + slc2_ind,nb,ne)) + endif + if (slc3_ind > 0) then + call read_vardata(dset_sfc, 'soill3', values_2d, errcode=iret) + if (iret /= 0) then + print *,'READGRIDDATA_PNC: error reading soill3' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + if (iope==0) call copytogrdin(ug,grdin(:,levels(n3d) + slc3_ind,nb,ne)) + endif + if (slc4_ind > 0) then + call read_vardata(dset_sfc, 'soill4', values_2d, errcode=iret) + if (iret /= 0) then + print *,'READGRIDDATA_PNC: error reading soill4' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + if (iope==0) call copytogrdin(ug,grdin(:,levels(n3d) + slc4_ind,nb,ne)) + endif + + ! bring all the subdomains back to the main PE + call mpi_barrier(iocomms(mem_pe(nproc)), iret) + if (allocated(values_2d)) deallocate(values_2d) + call close_dataset(dset_sfc) + call mpi_barrier(iocomms(mem_pe(nproc)), iret) + + end do sfcbackgroundloop ! loop over backgrounds to read in + end if !if (read_sfc_file) + ! remove the sub communicators call mpi_barrier(iocomms(mem_pe(nproc)), iret) call mpi_comm_free(iocomms(mem_pe(nproc)), iret) @@ -1322,6 +1457,20 @@ subroutine writegriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate_ integer(i_kind) :: ncstart(4), nccount(4) logical :: nocompress + logical :: write_sfc_file, write_atm_file + character(len=max_varname_length), dimension(n3d) :: no_vars3d + character(len=max_varname_length), dimension(n2d) :: no_vars2d + + call set_ncio_file_flags(vars3d, n3d, vars2d, n2d, write_sfc_file, write_atm_file) + + if (write_sfc_file ) then + ! adding the sfc increments requires adjusting several other variables. + ! This is done is a separate program. + if (nproc == 0) write(6,*) 'gridio/writegriddata_pnc: not coded to write sfc analysis, will write increment for sfc fields' + no_vars3d='' + call writeincrement_pnc(no_vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate_flag) + endif + nocompress = .true. if (nccompress) nocompress = .false. @@ -3616,6 +3765,7 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin, ! soil / snow mask (not fixed) integer(i_kind), dimension(nlons,nlats) :: mask logical :: write_sfc_file, write_atm_file + real(r_double) :: t1,t2 call set_ncio_file_flags(vars3d, n3d, vars2d, n2d, write_sfc_file, write_atm_file) @@ -3627,6 +3777,7 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin, nccount = (/nlons, nlats, nlevs/) if ( write_atm_file) then + if (nproc == 0) t1 = mpi_wtime() ne = 0 ensmemloop: do nanal=nanal1,nanal2 ne = ne + 1 @@ -3954,10 +4105,15 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin, end do backgroundloop ! loop over backgrounds to read in end do ensmemloop ! loop over ens members to read in + if (nproc == 0) then + t2 = mpi_wtime() + print *,'time in writeincrement atm_file on root',t2-t1,'secs' + endif endif ! write_atm_file if (write_sfc_file) then + if (nproc == 0) t1 = mpi_wtime() ne = 0 sfcensmemloop: do nanal=nanal1,nanal2 ne = ne + 1 @@ -4199,6 +4355,10 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin, end do sfcbackgroundloop ! loop over backgrounds to read in end do sfcensmemloop ! loop over ens members to read in + if (nproc == 0) then + t2 = mpi_wtime() + print *,'time in writeincrement sfc_file on root',t2-t1,'secs' + endif endif ! write_sfc_file @@ -4225,7 +4385,8 @@ end subroutine writeincrement subroutine writeincrement_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate_flag) use netcdf use params, only: nbackgrounds,incfileprefixes,fgfileprefixes,reducedgrid,& - datestring,nhr_anal + datestring,nhr_anal, & + incsfcfileprefixes,fgsfcfileprefixes use constants, only: grav,qcmin use mpi use module_ncio, only: Dataset, Variable, Dimension, open_dataset,& @@ -4257,12 +4418,17 @@ subroutine writeincrement_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate integer :: ql_ind, qi_ind, qr_ind, qs_ind, qg_ind ! netcdf things - integer(i_kind) :: dimids3(3),nccount(3),ncstart(3) + integer(i_kind) :: dimids3(3),nccount(3),ncstart(3), dimids2(2) integer(i_kind) :: ncid_out, lon_dimid, lat_dimid, lev_dimid, ilev_dimid integer(i_kind) :: lonvarid, latvarid, levvarid, pfullvarid, ilevvarid, & hyaivarid, hybivarid, uvarid, vvarid, delpvarid, delzvarid, & tvarid, sphumvarid, liqwatvarid, o3varid, icvarid, & - rwmrvarid, snmrvarid, grlevarid + rwmrvarid, snmrvarid, grlevarid, & + tmp2mvarid, spfh2mvarid, soilt1varid, soilt2varid, & + soilt3varid, soilt4varid, slc1varid, slc2varid, & + slc3varid, slc4varid, maskvarid + integer(i_kind) :: tmp2m_ind, spfh2m_ind, soilt1_ind, soilt2_ind,soilt3_ind, & + soilt4_ind,slc1_ind, slc2_ind, slc3_ind, slc4_ind integer(i_kind) :: iadateout ! fixed fields such as lat, lon, levs @@ -4274,11 +4440,20 @@ subroutine writeincrement_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate ! increment real(r_kind), dimension(nlons*nlats) :: psinc, inc, ug, vg, work real(r_single), allocatable, dimension(:,:,:) :: inc3d, inc3d2, inc3dout + real(r_single), allocatable, dimension(:,:) :: inc2d, inc2dout real(r_single), allocatable, dimension(:,:,:) :: tv, tvanl, tmp, tmpanl, q, qanl real(r_single), allocatable, dimension(:,:,:) :: q2, qanl2 real(r_kind), allocatable, dimension(:,:) :: values_2d real(r_kind), allocatable, dimension(:) :: psges, delzb, values_1d + ! soil / snow mask (not fixed) + integer(i_kind), dimension(nlons,nlats) :: mask + + logical :: write_sfc_file, write_atm_file + real(r_double) :: t1,t2 + + call set_ncio_file_flags(vars3d, n3d, vars2d, n2d, write_sfc_file, write_atm_file) + use_full_hydro = .false. clip = tiny_r_kind read(datestring,*) iadateout @@ -4317,6 +4492,9 @@ subroutine writeincrement_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate call mpi_bcast(grdin(1,1,nb,1),npts*ndim, mpi_real4, 0, iocomms(mem_pe(nproc)), iret) enddo + if (write_atm_file ) then + + if (nproc == 0) t1 = mpi_wtime() ! loop through times and do the read ne = 1 backgroundloop: do nb=1,nbackgrounds @@ -4772,8 +4950,261 @@ subroutine writeincrement_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate if (allocated(delzb)) deallocate(delzb) if (allocated(psges)) deallocate(psges) + !closing file + call nccheck_incr(nf90_close(ncid_out)) end do backgroundloop ! loop over backgrounds to write out + if (nproc == 0) then + t2 = mpi_wtime() + print *,'time in writeincrement_pnc atm_file on root',t2-t1,'secs' + endif + end if ! if (write_atm_file) + + if (write_sfc_file ) then + + if (nproc == 0) t1 = mpi_wtime() + + tmp2m_ind = getindex(vars2d, 't2m') !< indices in the state or control var arrays + spfh2m_ind = getindex(vars2d, 'q2m') + soilt1_ind = getindex(vars2d, 'st1') + slc1_ind = getindex(vars2d, 'sl1') + soilt2_ind = getindex(vars2d, 'st2') + slc2_ind = getindex(vars2d, 'sl2') + soilt3_ind = getindex(vars2d, 'st3') + slc3_ind = getindex(vars2d, 'sl3') + soilt4_ind = getindex(vars2d, 'st4') + slc4_ind = getindex(vars2d, 'sl4') + + ! loop through times and do the read + ne = 1 + write(charnanal,'(i3.3)') nanal + sfcbackgroundloop: do nb=1,nbackgrounds + + if(no_inflate_flag) then + filenameout = trim(adjustl(datapath))//trim(adjustl(incsfcfileprefixes(nb)))//"nimem"//charnanal + else + filenameout = trim(adjustl(datapath))//trim(adjustl(incsfcfileprefixes(nb)))//"mem"//charnanal + end if + filenamein = trim(adjustl(datapath))//trim(adjustl(fgsfcfileprefixes(nb)))//"mem"//charnanal + + !! note: only iope=0 is writing the outputs. Having all pes in iocomm write to a file slows it down. + !! + if (iope==0) then + dsfg = open_dataset(filenamein) + ! create the output netCDF increment file + call nccheck_incr(nf90_create(path=trim(filenameout), cmode=nf90_netcdf4,ncid=ncid_out)) + + ! create dimensions based on analysis resolution, not guess + call nccheck_incr(nf90_def_dim(ncid_out, "longitude", nlons, lon_dimid)) + call nccheck_incr(nf90_def_dim(ncid_out, "latitude", nlats, lat_dimid)) + dimids2 = (/ lon_dimid, lat_dimid /) + ! create variables + call nccheck_incr(nf90_def_var(ncid_out, "longitude", nf90_real,(/lon_dimid/), lonvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "latitude", nf90_real,(/lat_dimid/), latvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "tmp2m_inc", nf90_real, dimids2,tmp2mvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "spfh2m_inc", nf90_real, dimids2,spfh2mvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilt1_inc", nf90_real, dimids2,soilt1varid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilt2_inc", nf90_real, dimids2,soilt2varid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilt3_inc", nf90_real, dimids2,soilt3varid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilt4_inc", nf90_real, dimids2,soilt4varid)) + call nccheck_incr(nf90_def_var(ncid_out, "slc1_inc", nf90_real, dimids2,slc1varid)) + call nccheck_incr(nf90_def_var(ncid_out, "slc2_inc", nf90_real, dimids2,slc2varid)) + call nccheck_incr(nf90_def_var(ncid_out, "slc3_inc", nf90_real, dimids2,slc3varid)) + call nccheck_incr(nf90_def_var(ncid_out, "slc4_inc", nf90_real, dimids2,slc4varid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilsnow_mask", nf90_int,dimids2, maskvarid)) + ! place global attributes to serial calc_increment output + call nccheck_incr(nf90_put_att(ncid_out, nf90_global, "source", "GSI EnKF")) + call nccheck_incr(nf90_put_att(ncid_out, nf90_global, "comment", & + "global landsfc anal increment from writeincrement")) + call nccheck_incr(nf90_put_att(ncid_out, nf90_global, "analysis_time",iadateout)) + call nccheck_incr(nf90_put_att(ncid_out, nf90_global,"IAU_hour_from_guess", nhr_anal(nb))) + ! add units to lat/lon because that's what the calc_increment utility has + call nccheck_incr(nf90_put_att(ncid_out, lonvarid, "units","degrees_east")) + call nccheck_incr(nf90_put_att(ncid_out, latvarid, "units","degrees_north")) + ! end the netCDF file definition + call nccheck_incr(nf90_enddef(ncid_out)) + + ! longitudes + call read_vardata(dsfg, 'grid_xt', values_1d, errcode=iret) + deglons(:) = values_1d + call nccheck_incr(nf90_put_var(ncid_out, lonvarid, deglons, & + start = (/1/), count = (/nlons/))) + + call read_vardata(dsfg, 'grid_yt', values_1d, errcode=iret) + ! latitudes + do j=1,nlats + deglats(nlats-j+1) = values_1d(j) + end do + call nccheck_incr(nf90_put_var(ncid_out, latvarid, deglats, & + start = (/1/), count = (/nlats/))) + ! construct mask (1 - soil, 2 - snow, 0 - not snow) + ! note: same logic/threshold used in global_cycle to produce + ! mask on model grid. + call read_vardata(dsfg, 'soill1', values_2d, errcode=iret) + mask = 0 + do j=1,nlats + do i = 1, nlons + if (values_2d(i,j) .LT. 1.0) then + mask(i,nlats-j+1) = 1 + endif + enddo + end do + call read_vardata(dsfg, 'weasd', values_2d, errcode=iret) + do j=1,nlats + do i = 1, nlons + if (values_2d(i,j) .GT. 0.001) then + mask(i,nlats-j+1) = 2 + endif + enddo + end do + call nccheck_incr(nf90_put_var(ncid_out, maskvarid, mask, & + start = ncstart(1:2), count = nccount(1:2))) + + allocate(inc2d(nlons,nlats)) + allocate(inc2dout(nlons,nlats)) + + ! tmp2m increment + inc(:) = zero + if (tmp2m_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d) + tmp2m_ind,nb,ne),inc) + endif + inc2d(:,:) = reshape(inc,(/nlons,nlats/)) + do j=1,nlats + inc2dout(:,nlats-j+1) = inc2d(:,j) + end do + call nccheck_incr(nf90_put_var(ncid_out, tmp2mvarid, sngl(inc2dout), & + start = ncstart(1:2), count = nccount(1:2))) + ! spfh2m increment + inc(:) = zero + if (spfh2m_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d)+spfh2m_ind,nb,ne),inc) + endif + inc2d(:,:) = reshape(inc,(/nlons,nlats/)) + do j=1,nlats + inc2dout(:,nlats-j+1) = inc2d(:,j) + end do + call nccheck_incr(nf90_put_var(ncid_out, spfh2mvarid, sngl(inc2dout), & + start = ncstart(1:2), count = nccount(1:2))) + ! soilt1 increment + inc(:) = zero + if (soilt1_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d)+soilt1_ind,nb,ne),inc) + endif + inc2d(:,:) = reshape(inc,(/nlons,nlats/)) + inc2dout=0. + do j=1,nlats + do i = 1, nlons + if (mask(i,nlats-j+1) .NE. 0) inc2dout(i,nlats-j+1) = inc2d(i,j) + enddo + end do + call nccheck_incr(nf90_put_var(ncid_out, soilt1varid, sngl(inc2dout), & + start = ncstart(1:2), count = nccount(1:2))) + ! soilt2 increment + inc(:) = zero + if (soilt2_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d)+soilt2_ind,nb,ne),inc) + endif + inc2d(:,:) = reshape(inc,(/nlons,nlats/)) + inc2dout=0. + do j=1,nlats + do i = 1, nlons + if (mask(i,nlats-j+1) .NE. 0) inc2dout(i,nlats-j+1) = inc2d(i,j) + enddo + end do + call nccheck_incr(nf90_put_var(ncid_out, soilt2varid, sngl(inc2dout), & + start = ncstart(1:2), count = nccount(1:2))) + ! soilt3 increment + inc(:) = zero + if (soilt3_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d)+soilt3_ind,nb,ne),inc) + endif + inc2d(:,:) = reshape(inc,(/nlons,nlats/)) + inc2dout=0. + do j=1,nlats + do i = 1, nlons + if (mask(i,nlats-j+1) .NE. 0) inc2dout(i,nlats-j+1) = inc2d(i,j) + enddo + end do + call nccheck_incr(nf90_put_var(ncid_out, soilt3varid, sngl(inc2dout), & + start = ncstart(1:2), count = nccount(1:2))) + ! soilt4 increment + inc(:) = zero + if (soilt4_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d)+soilt4_ind,nb,ne),inc) + endif + inc2d(:,:) = reshape(inc,(/nlons,nlats/)) + inc2dout=0. + do j=1,nlats + do i = 1, nlons + if (mask(i,nlats-j+1) .NE. 0) inc2dout(i,nlats-j+1) = inc2d(i,j) + enddo + end do + call nccheck_incr(nf90_put_var(ncid_out, soilt4varid, sngl(inc2dout), & + start = ncstart(1:2), count = nccount(1:2))) + ! slc1 increment + inc(:) = zero + if (slc1_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d)+slc1_ind,nb,ne),inc) + endif + inc2d(:,:) = reshape(inc,(/nlons,nlats/)) + do j=1,nlats + inc2dout(:,nlats-j+1) = inc2d(:,j) + end do + call nccheck_incr(nf90_put_var(ncid_out, slc1varid, sngl(inc2dout), & + start = ncstart(1:2), count = nccount(1:2))) + ! slc2 increment + inc(:) = zero + if (slc2_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d)+slc2_ind,nb,ne),inc) + endif + inc2d(:,:) = reshape(inc,(/nlons,nlats/)) + do j=1,nlats + inc2dout(:,nlats-j+1) = inc2d(:,j) + end do + call nccheck_incr(nf90_put_var(ncid_out, slc2varid, sngl(inc2dout), & + start = ncstart(1:2), count = nccount(1:2))) + ! slc3 increment + inc(:) = zero + if (slc3_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d)+slc3_ind,nb,ne),inc) + endif + inc2d(:,:) = reshape(inc,(/nlons,nlats/)) + do j=1,nlats + inc2dout(:,nlats-j+1) = inc2d(:,j) + end do + call nccheck_incr(nf90_put_var(ncid_out, slc3varid, sngl(inc2dout), & + start = ncstart(1:2), count = nccount(1:2))) + ! slc4 increment + inc(:) = zero + if (slc4_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d)+slc4_ind,nb,ne),inc) + endif + inc2d(:,:) = reshape(inc,(/nlons,nlats/)) + do j=1,nlats + inc2dout(:,nlats-j+1) = inc2d(:,j) + end do + call nccheck_incr(nf90_put_var(ncid_out, slc4varid, sngl(inc2dout), & + start = ncstart(1:2), count = nccount(1:2))) + + call close_dataset(dsfg,errcode=iret) + if (iret/=0) then + write(6,*)'gridio/writeincrement_par: problem closing netcdf sfc fg dataset, iret=',iret + call stop2(23) + endif + ! deallocate things + deallocate(inc2d,inc2dout) + + call nccheck_incr(nf90_close(ncid_out)) + + end if + + end do sfcbackgroundloop ! loop over backgrounds to read in + if (nproc == 0) then + t2 = mpi_wtime() + print *,'time in writeincrement_pnc sfc_file on root',t2-t1,'secs' + endif + endif !write_Sfc + ! remove the sub communicators call mpi_barrier(iocomms(mem_pe(nproc)), iret) call mpi_comm_free(iocomms(mem_pe(nproc)), iret) From 2167bc93fe2b1e84657822a323666ccc61f9a339 Mon Sep 17 00:00:00 2001 From: TingLei-NOAA Date: Wed, 20 Mar 2024 12:16:42 -0400 Subject: [PATCH 6/6] =?UTF-8?q?Issue=20694:=20Upgrade/refactoring=20for=20?= =?UTF-8?q?U=20and=20V=20write-out=20sub=20for=20FV3REG=20GSI=20failure=20?= =?UTF-8?q?=E2=80=A6=20(#698)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit **DUE DATE for merger of this PR into `develop` is 3/27/2024 (six weeks after PR creation).** Resolves #693 (Thanks to @edwardhartnett 's suggestions) Resolves # 694 ( this PR is not able to provide a stable solution, more details would be given on the issue page) Resolves # 697: With larger requested memory for each mpi task, it still showed, for some time, the differences in the analysis files between loproc vs hiproc for the control runs on hercules. whether integrating this with the refactored IO part would provide a stable solution remains to be seen. This PR resolved the newly emerged issue with IO of netcdf files in the continuous storage, with upgraded FV3REG IO for the cold start options. (Co author Ming Hu @hu5970 ) This PR is being worked on in collaboration with Pete Johnson through RDHPCS help desk, @RussTreadon-NOAA @DavidHuber-NOAA and thanks to help from @ed Raghue Reddy through RDHPCS help desk. --------- Co-authored-by: Ting Lei Co-authored-by: Ting.Lei-NOAA --- regression/regression_driver.sh | 2 + regression/regression_param.sh | 48 +-- src/gsi/gsi_rfv3io_mod.f90 | 542 ++++++++++++++++++++------------ ush/sub_hera | 5 +- ush/sub_hercules | 4 +- ush/sub_jet | 2 +- ush/sub_orion | 2 +- 7 files changed, 365 insertions(+), 240 deletions(-) diff --git a/regression/regression_driver.sh b/regression/regression_driver.sh index 821cc7cedb..38329778a4 100755 --- a/regression/regression_driver.sh +++ b/regression/regression_driver.sh @@ -36,10 +36,12 @@ for jn in `seq ${RSTART} ${REND}`; do export scripts=${scripts_updat:-$scripts} export fixgsi=${fixgsi_updat:-$fixgsi} export modulefiles=${modulefiles_updat:-$modulefiles} + export ush=${ush_update:-$ush} else export scripts=${scripts_contrl:-$scripts} export fixgsi=${fixgsi_contrl:-$fixgsi} export modulefiles=${modulefiles_contrl:-$modulefiles} + export ush=${ush_cntrl:-$ush} fi rm -f ${job[$jn]}.out diff --git a/regression/regression_param.sh b/regression/regression_param.sh index 2ac615fc4a..a4f5d7035c 100755 --- a/regression/regression_param.sh +++ b/regression/regression_param.sh @@ -87,17 +87,17 @@ case $regtest in rrfs_3denvar_glbens) if [[ "$machine" = "Hera" ]]; then - topts[1]="0:15:00" ; popts[1]="20/1/" ; ropts[1]="/1" - topts[2]="0:15:00" ; popts[2]="20/2/" ; ropts[2]="/1" + topts[1]="0:15:00" ; popts[1]="5/4/" ; ropts[1]="/1" + topts[2]="0:15:00" ; popts[2]="10/4/" ; ropts[2]="/1" elif [[ "$machine" = "Orion" ]]; then - topts[1]="0:15:00" ; popts[1]="20/1/" ; ropts[1]="/1" - topts[2]="0:15:00" ; popts[2]="20/2/" ; ropts[2]="/2" + topts[1]="0:15:00" ; popts[1]="5/4/" ; ropts[1]="/1" + topts[2]="0:15:00" ; popts[2]="10/4/" ; ropts[2]="/2" elif [[ "$machine" = "Hercules" ]]; then - topts[1]="0:15:00" ; popts[1]="20/1/" ; ropts[1]="/1" - topts[2]="0:15:00" ; popts[2]="20/2/" ; ropts[2]="/2" + topts[1]="0:15:00" ; popts[1]="5/4/" ; ropts[1]="/1" + topts[2]="0:15:00" ; popts[2]="10/4/" ; ropts[2]="/2" elif [[ "$machine" = "Jet" ]]; then - topts[1]="0:15:00" ; popts[1]="20/1/" ; ropts[1]="/1" - topts[2]="0:15:00" ; popts[2]="20/2/" ; ropts[2]="/1" + topts[1]="0:15:00" ; popts[1]="5/4/" ; ropts[1]="/1" + topts[2]="0:15:00" ; popts[2]="10/4/" ; ropts[2]="/1" elif [[ "$machine" = "Gaea" ]]; then topts[1]="0:15:00" ; popts[1]="18/1/" ; ropts[1]="/1" topts[2]="0:15:00" ; popts[2]="18/2/" ; ropts[2]="/1" @@ -117,17 +117,17 @@ case $regtest in hafs_3denvar_hybens) if [[ "$machine" = "Hera" ]]; then - topts[1]="0:15:00" ; popts[1]="20/1/" ; ropts[1]="/1" - topts[2]="0:15:00" ; popts[2]="20/2/" ; ropts[2]="/1" + topts[1]="0:15:00" ; popts[1]="5/4/" ; ropts[1]="/1" + topts[2]="0:15:00" ; popts[2]="10/4/" ; ropts[2]="/1" elif [[ "$machine" = "Orion" ]]; then - topts[1]="0:15:00" ; popts[1]="20/1/" ; ropts[1]="/1" - topts[2]="0:15:00" ; popts[2]="20/2/" ; ropts[2]="/2" + topts[1]="0:15:00" ; popts[1]="5/4/" ; ropts[1]="/1" + topts[2]="0:15:00" ; popts[2]="10/4/" ; ropts[2]="/2" elif [[ "$machine" = "Hercules" ]]; then topts[1]="0:15:00" ; popts[1]="20/1/" ; ropts[1]="/1" - topts[2]="0:15:00" ; popts[2]="20/2/" ; ropts[2]="/2" + topts[2]="0:15:00" ; popts[2]="5/8/" ; ropts[2]="/2" elif [[ "$machine" = "Jet" ]]; then - topts[1]="0:15:00" ; popts[1]="20/1/" ; ropts[1]="/1" - topts[2]="0:15:00" ; popts[2]="20/2/" ; ropts[2]="/1" + topts[1]="0:15:00" ; popts[1]="5/4/" ; ropts[1]="/1" + topts[2]="0:15:00" ; popts[2]="10/4/" ; ropts[2]="/1" elif [[ "$machine" = "Gaea" ]]; then topts[1]="0:15:00" ; popts[1]="18/1/" ; ropts[1]="/1" topts[2]="0:15:00" ; popts[2]="18/2/" ; ropts[2]="/1" @@ -146,17 +146,17 @@ case $regtest in hafs_4denvar_glbens) if [[ "$machine" = "Hera" ]]; then - topts[1]="0:15:00" ; popts[1]="20/1/" ; ropts[1]="/1" - topts[2]="0:15:00" ; popts[2]="20/2/" ; ropts[2]="/1" + topts[1]="0:15:00" ; popts[1]="5/4/" ; ropts[1]="/1" + topts[2]="0:15:00" ; popts[2]="10/4/" ; ropts[2]="/1" elif [[ "$machine" = "Orion" ]]; then - topts[1]="0:20:00" ; popts[1]="20/1/" ; ropts[1]="/1" + topts[1]="0:20:00" ; popts[1]="5/4/" ; ropts[1]="/1" topts[2]="0:20:00" ; popts[2]="20/2/" ; ropts[2]="/1" elif [[ "$machine" = "Hercules" ]]; then - topts[1]="0:20:00" ; popts[1]="20/1/" ; ropts[1]="/1" - topts[2]="0:20:00" ; popts[2]="20/2/" ; ropts[2]="/1" + topts[1]="0:20:00" ; popts[1]="5/4/" ; ropts[1]="/1" + topts[2]="0:20:00" ; popts[2]="10/4/" ; ropts[2]="/1" elif [[ "$machine" = "Jet" ]]; then - topts[1]="0:15:00" ; popts[1]="20/1/" ; ropts[1]="/1" - topts[2]="0:15:00" ; popts[2]="20/2/" ; ropts[2]="/1" + topts[1]="0:15:00" ; popts[1]="5/4/" ; ropts[1]="/1" + topts[2]="0:15:00" ; popts[2]="10/4/" ; ropts[2]="/1" elif [[ "$machine" = "Gaea" ]]; then topts[1]="0:15:00" ; popts[1]="18/1/" ; ropts[1]="/1" topts[2]="0:15:00" ; popts[2]="18/2/" ; ropts[2]="/1" @@ -300,10 +300,10 @@ if [[ "$machine" = "Hera" ]]; then export APRUN="srun" elif [[ "$machine" = "Orion" ]]; then export OMP_STACKSIZE=2048M - export APRUN="srun -n \$ntasks --cpus-per-task=\$threads" + export APRUN="srun -n \$ntasks --mem=0 --cpus-per-task=\$threads" elif [[ "$machine" = "Hercules" ]]; then export OMP_STACKSIZE=2048M - export APRUN="srun -n \$ntasks --cpus-per-task=\$threads" + export APRUN="srun -n \$ntasks --mem=0 --cpus-per-task=\$threads" elif [[ "$machine" = "Jet" ]]; then export OMP_STACKSIZE=1024M export MPI_BUFS_PER_PROC=256 diff --git a/src/gsi/gsi_rfv3io_mod.f90 b/src/gsi/gsi_rfv3io_mod.f90 index b3f3488a70..652bad9a33 100644 --- a/src/gsi/gsi_rfv3io_mod.f90 +++ b/src/gsi/gsi_rfv3io_mod.f90 @@ -547,7 +547,7 @@ subroutine gsi_rfv3io_get_ens_grid_specs(grid_spec,ierr) ! !$$$ end documentation block use netcdf, only: nf90_open,nf90_close,nf90_get_var,nf90_noerr - use netcdf, only: nf90_nowrite,nf90_inquire,nf90_inquire_dimension + use netcdf, only: nf90_nowrite,nf90_mpiio,nf90_inquire,nf90_inquire_dimension use netcdf, only: nf90_inquire_variable use mpimod, only: mype use mod_fv3_lola, only: definecoef_regular_grids @@ -1238,7 +1238,6 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) fv3lam_io_phymetvars3d_nouv(jphyvar)=trim(vartem) else write(6,*)'the metvarname ',vartem,' is not expected, stop' - call flush(6) call stop2(333) endif endif @@ -1253,7 +1252,6 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) endif if(jdynvar /= ndynvario3d.or.jtracer /= ntracerio3d.or.jphyvar /= nphyvario3d ) then write(6,*)'ndynvario3d is not as expected, stop' - call flush(6) call stop2(333) endif if(mype == 0) then @@ -1361,7 +1359,6 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) else if (trim(vartem) /= "pm2_5")then write(6,*)'the chemvarname ',vartem,' is not in aeronames_smoke_fv3 !!!' - call flush(6) endif endif enddo @@ -1598,7 +1595,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) if( fv3sar_bg_opt == 0) then call gsi_fv3ncdf_readuv(grd_fv3lam_uv,ges_u,ges_v,fv3filenamegin(it),.false.) else - call gsi_fv3ncdf_readuv_v1(grd_fv3lam_uv,ges_u,ges_v,fv3filenamegin(it),.false.) + call gsi_fv3ncdf_readuv_v1(grd_fv3lam_uv,ges_u,ges_v,fv3filenamegin(it),.false.) endif if( fv3sar_bg_opt == 0) then @@ -2449,7 +2446,7 @@ subroutine gsi_fv3ncdf_read(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin,ens use mpimod, only: mpi_comm_world,mpi_rtype,mype,npe,setcomm,mpi_integer,mpi_max use mpimod, only: MPI_INFO_NULL use netcdf, only: nf90_open,nf90_close,nf90_get_var,nf90_noerr - use netcdf, only: nf90_nowrite,nf90_inquire,nf90_inquire_dimension + use netcdf, only: nf90_nowrite,nf90_mpiio,nf90_inquire,nf90_inquire_dimension use netcdf, only: nf90_inquire_variable use netcdf, only: nf90_inq_varid use mod_fv3_lola, only: fv3_h_to_ll,fv3_h_to_ll_ens @@ -2469,9 +2466,9 @@ subroutine gsi_fv3ncdf_read(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin,ens character(len=max_varname_length) :: name character(len=max_filename_length) :: filenamein2 real(r_kind),allocatable,dimension(:,:):: uu2d_tmp - integer(i_kind) :: countloc_tmp(3),startloc_tmp(3) + integer(i_kind) :: countloc_tmp(4),startloc_tmp(4) - integer(i_kind) nlatcase,nloncase,nxcase,nycase,countloc(3),startloc(3) + integer(i_kind) nlatcase,nloncase,nxcase,nycase,countloc(4),startloc(4) integer(i_kind) ilev,ilevtot,inative integer(i_kind) kbgn,kend,len logical :: phy_smaller_domain @@ -2528,18 +2525,16 @@ subroutine gsi_fv3ncdf_read(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin,ens allocate(gfile_loc_layout(0:fv3_io_layout_y-1)) do nio=0,fv3_io_layout_y-1 write(filename_layout,'(a,a,I4.4)') trim(filenamein),'.',nio - iret=nf90_open(filename_layout,nf90_nowrite,gfile_loc_layout(nio),comm=mpi_comm_read,info=MPI_INFO_NULL) !clt + iret=nf90_open(filename_layout,ior(nf90_nowrite,nf90_mpiio),gfile_loc_layout(nio),comm=mpi_comm_read,info=MPI_INFO_NULL) !clt if(iret/=nf90_noerr) then write(6,*)' gsi_fv3ncdf_read: problem opening ',trim(filename_layout),gfile_loc_layout(nio),', Status = ',iret - call flush(6) call stop2(333) endif enddo else - iret=nf90_open(filenamein,nf90_nowrite,gfile_loc,comm=mpi_comm_read,info=MPI_INFO_NULL) !clt + iret=nf90_open(filenamein,ior(nf90_nowrite,nf90_mpiio),gfile_loc,comm=mpi_comm_read,info=MPI_INFO_NULL) !clt if(iret/=nf90_noerr) then write(6,*)' gsi_fv3ncdf_read: problem opening ',trim(filenamein),gfile_loc,', Status = ',iret - call flush(6) call stop2(333) endif endif @@ -2554,15 +2549,14 @@ subroutine gsi_fv3ncdf_read(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin,ens name=trim(varname) if(trim(filenamein) /= trim(filenamein2)) then write(6,*)'filenamein and filenamein2 are not the same as expected, stop' - call flush(6) call stop2(333) endif ilev=grd_ionouv%lnames(1,ilevtot) nz=grd_ionouv%nsig nzp1=nz+1 inative=nzp1-ilev - startloc=(/1,1,inative/) - countloc=(/nxcase,nycase,1/) + startloc=(/1,1,inative,1/) + countloc=(/nxcase,nycase,1,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' .or. trim(adjustl(varname)) == 'flash_extent_density' )then @@ -2570,23 +2564,23 @@ subroutine gsi_fv3ncdf_read(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin,ens if(trim(name)=='xaxis_1') nx_phy=len if( nx_phy == nxcase )then allocate(uu2d_tmp(nxcase,nycase)) - countloc_tmp=(/nxcase,nycase,1/) + countloc_tmp=(/nxcase,nycase,1,1/) phy_smaller_domain = .false. else allocate(uu2d_tmp(nxcase-6,nycase-6)) - countloc_tmp=(/nxcase-6,nycase-6,1/) + countloc_tmp=(/nxcase-6,nycase-6,1,1/) phy_smaller_domain = .true. end if - startloc_tmp=(/1,1,ilev/) + startloc_tmp=(/1,1,ilev,1/) end if if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 if (ensgrid) then - countloc=(/nxcase,ny_layout_lenens(nio)+1,1/) + countloc=(/nxcase,ny_layout_lenens(nio)+1,1,1/) allocate(uu2d_layout(nxcase,ny_layout_lenens(nio)+1)) else - countloc=(/nxcase,ny_layout_len(nio),1/) + countloc=(/nxcase,ny_layout_len(nio),1,1/) allocate(uu2d_layout(nxcase,ny_layout_len(nio))) end if iret=nf90_inq_varid(gfile_loc_layout(nio),trim(adjustl(varname)),var_id) @@ -2671,10 +2665,10 @@ subroutine gsi_fv3ncdf_read_v1(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin, use kinds, only: r_kind,i_kind - use mpimod, only: mpi_rtype,mpi_comm_world,mype,MPI_INFO_NULL - use mpimod, only: mpi_comm_world,mpi_rtype,mype + use mpimod, only: npe,mpi_rtype,mpi_comm_world,mype,MPI_INFO_NULL + use mpimod, only: mpi_comm_world,mpi_rtype,mype,setcomm,mpi_integer,mpi_max use netcdf, only: nf90_open,nf90_close,nf90_get_var,nf90_noerr - use netcdf, only: nf90_nowrite,nf90_inquire,nf90_inquire_dimension + use netcdf, only: nf90_nowrite,nf90_mpiio,nf90_inquire,nf90_inquire_dimension use netcdf, only: nf90_inquire_variable use netcdf, only: nf90_inq_varid use mod_fv3_lola, only: fv3_h_to_ll,fv3_h_to_ll_ens @@ -2694,12 +2688,18 @@ subroutine gsi_fv3ncdf_read_v1(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin, character(len=max_varname_length) :: varname,vgsiname - integer(i_kind) nlatcase,nloncase,nxcase,nycase,countloc(3),startloc(3) + integer(i_kind) nlatcase,nloncase,nxcase,nycase,countloc(4),startloc(4) integer(i_kind) kbgn,kend integer(i_kind) var_id integer(i_kind) inative,ilev,ilevtot integer(i_kind) gfile_loc,iret integer(i_kind) nzp1,mm1 + + integer(i_kind):: iworld,iworld_group,nread,mpi_comm_read,i,ierror + integer(i_kind),dimension(npe):: members,members_read,mype_read_rank + logical:: procuse + + mm1=mype+1 @@ -2712,13 +2712,34 @@ subroutine gsi_fv3ncdf_read_v1(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin, nxcase=nx nycase=ny end if + allocate(uu2d(nxcase,nycase)) + kbgn=grd_ionouv%kbegin_loc kend=grd_ionouv%kend_loc - allocate(uu2d(nxcase,nycase)) - iret=nf90_open(filenamein,nf90_nowrite,gfile_loc,comm=mpi_comm_world,info=MPI_INFO_NULL) !clt + procuse = .false. + members=-1 + members_read=-1 + if (kbgn<=kend) then + procuse = .true. + members(mm1) = mype + endif + call mpi_allreduce(members,members_read,npe,mpi_integer,mpi_max,mpi_comm_world,ierror) + + nread=0 + mype_read_rank=-1 + do i=1,npe + if (members_read(i) >= 0) then + nread=nread+1 + mype_read_rank(nread) = members_read(i) + endif + enddo + + call setcomm(iworld,iworld_group,nread,mype_read_rank,mpi_comm_read,ierror) + + if (procuse) then + iret=nf90_open(filenamein,ior(nf90_nowrite,nf90_mpiio),gfile_loc,comm=mpi_comm_read,info=MPI_INFO_NULL) !clt if(iret/=nf90_noerr) then write(6,*)' gsi_fv3ncdf_read_v1: problem opening ',trim(filenamein),gfile_loc,', Status = ',iret - call flush(6) call stop2(333) endif @@ -2728,15 +2749,14 @@ subroutine gsi_fv3ncdf_read_v1(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin, call getfv3lamfilevname(vgsiname,fv3filenamegin,filenamein2,varname) if(trim(filenamein) /= trim(filenamein2)) then write(6,*)'filenamein and filenamein2 are not the same as expected, stop' - call flush(6) call stop2(333) endif ilev=grd_ionouv%lnames(1,ilevtot) nz=grd_ionouv%nsig nzp1=nz+1 inative=nzp1-ilev - startloc=(/1,1,inative+1/) - countloc=(/nxcase,nycase,1/) + startloc=(/1,1,inative+1,1/) + countloc=(/nxcase,nycase,1,1/) iret=nf90_inq_varid(gfile_loc,trim(adjustl(varname)),var_id) if(iret/=nf90_noerr) then write(6,*)' wrong to get var_id ',var_id @@ -2752,8 +2772,9 @@ subroutine gsi_fv3ncdf_read_v1(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin, end if enddo ! i - call general_grid2sub(grd_ionouv,hwork,cstate_nouv%values) iret=nf90_close(gfile_loc) + endif + call general_grid2sub(grd_ionouv,hwork,cstate_nouv%values) deallocate (uu2d) @@ -2785,7 +2806,7 @@ subroutine gsi_fv3ncdf_readuv(grd_uv,ges_u,ges_v,fv3filenamegin,ensgrid) use kinds, only: r_kind,i_kind use mpimod, only: mpi_comm_world,mpi_rtype,mype,mpi_info_null,npe,setcomm,mpi_integer,mpi_max use netcdf, only: nf90_open,nf90_close,nf90_get_var,nf90_noerr - use netcdf, only: nf90_nowrite,nf90_inquire,nf90_inquire_dimension + use netcdf, only: nf90_nowrite,nf90_mpiio,nf90_inquire,nf90_inquire_dimension use netcdf, only: nf90_inquire_variable use netcdf, only: nf90_inq_varid use mod_fv3_lola, only: fv3_h_to_ll,fv3uv2earth,fv3_h_to_ll_ens,fv3uv2earthens @@ -2808,7 +2829,7 @@ subroutine gsi_fv3ncdf_readuv(grd_uv,ges_u,ges_v,fv3filenamegin,ensgrid) integer(i_kind) u_grd_VarId,v_grd_VarId integer(i_kind) nlatcase,nloncase integer(i_kind) nxcase,nycase - integer(i_kind) u_countloc(3),u_startloc(3),v_countloc(3),v_startloc(3) + integer(i_kind) u_countloc(4),u_startloc(4),v_countloc(4),v_startloc(4) integer(i_kind) inative,ilev,ilevtot integer(i_kind) kbgn,kend @@ -2873,15 +2894,13 @@ subroutine gsi_fv3ncdf_readuv(grd_uv,ges_u,ges_v,fv3filenamegin,ensgrid) iret=nf90_open(filename_layout,nf90_nowrite,gfile_loc_layout(nio),comm=mpi_comm_read,info=MPI_INFO_NULL) if(iret/=nf90_noerr) then write(6,*)'problem opening6 ',trim(filename_layout),gfile_loc_layout(nio),', Status = ',iret - call flush(6) call stop2(333) endif enddo else - iret=nf90_open(filenamein,nf90_nowrite,gfile_loc,comm=mpi_comm_read,info=MPI_INFO_NULL) !clt + iret=nf90_open(filenamein,ior(nf90_nowrite,nf90_mpiio),gfile_loc,comm=mpi_comm_read,info=MPI_INFO_NULL) !clt if(iret/=nf90_noerr) then write(6,*)' problem opening6 ',trim(filenamein),', Status = ',iret - call flush(6) call stop2(333) endif endif @@ -2891,24 +2910,23 @@ subroutine gsi_fv3ncdf_readuv(grd_uv,ges_u,ges_v,fv3filenamegin,ensgrid) call getfv3lamfilevname(vgsiname,fv3filenamegin,filenamein2,varname) if(trim(filenamein) /= trim(filenamein2)) then write(6,*)'filenamein and filenamein2 are not the same as expected, stop' - call flush(6) call stop2(333) endif ilev=grd_uv%lnames(1,ilevtot) nz=grd_uv%nsig nzp1=nz+1 inative=nzp1-ilev - u_countloc=(/nxcase,nycase+1,1/) - v_countloc=(/nxcase+1,nycase,1/) - u_startloc=(/1,1,inative/) - v_startloc=(/1,1,inative/) + u_countloc=(/nxcase,nycase+1,1,1/) + v_countloc=(/nxcase+1,nycase,1,1/) + u_startloc=(/1,1,inative,1/) + v_startloc=(/1,1,inative,1/) if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 if (ensgrid) then - u_countloc=(/nxcase,ny_layout_lenens(nio)+1,1/) + u_countloc=(/nxcase,ny_layout_lenens(nio)+1,1,1/) allocate(u2d_layout(nxcase,ny_layout_lenens(nio)+1)) else - u_countloc=(/nxcase,ny_layout_len(nio)+1,1/) + u_countloc=(/nxcase,ny_layout_len(nio)+1,1,1/) allocate(u2d_layout(nxcase,ny_layout_len(nio)+1)) end if call check( nf90_inq_varid(gfile_loc_layout(nio),'u',u_grd_VarId) ) @@ -2917,13 +2935,13 @@ subroutine gsi_fv3ncdf_readuv(grd_uv,ges_u,ges_v,fv3filenamegin,ensgrid) u2d(:,ny_layout_bens(nio):ny_layout_eens(nio))=u2d_layout(:,1:ny_layout_lenens(nio)) if(nio==fv3_io_layout_y-1) u2d(:,ny_layout_eens(nio)+1)=u2d_layout(:,ny_layout_lenens(nio)+1) deallocate(u2d_layout) - v_countloc=(/nxcase+1,ny_layout_lenens(nio),1/) + v_countloc=(/nxcase+1,ny_layout_lenens(nio),1,1/) allocate(v2d_layout(nxcase+1,ny_layout_lenens(nio))) else u2d(:,ny_layout_b(nio):ny_layout_e(nio))=u2d_layout(:,1:ny_layout_len(nio)) if(nio==fv3_io_layout_y-1) u2d(:,ny_layout_e(nio)+1)=u2d_layout(:,ny_layout_len(nio)+1) deallocate(u2d_layout) - v_countloc=(/nxcase+1,ny_layout_len(nio),1/) + v_countloc=(/nxcase+1,ny_layout_len(nio),1,1/) allocate(v2d_layout(nxcase+1,ny_layout_len(nio))) end if call check( nf90_inq_varid(gfile_loc_layout(nio),'v',v_grd_VarId) ) @@ -3019,9 +3037,10 @@ subroutine gsi_fv3ncdf_readuv_v1(grd_uv,ges_u,ges_v,fv3filenamegin,ensgrid) !$$$ end documentation block use constants, only: half use kinds, only: r_kind,i_kind - use mpimod, only: mpi_comm_world,mpi_rtype,mype,mpi_info_null + use mpimod, only: setcomm,mpi_integer,mpi_max, npe,mpi_comm_world,mpi_rtype,mype,mpi_info_null use netcdf, only: nf90_open,nf90_close,nf90_get_var,nf90_noerr - use netcdf, only: nf90_nowrite,nf90_inquire,nf90_inquire_dimension + use netcdf, only: nf90_nowrite,nf90_mpiio,nf90_inquire,nf90_inquire_dimension + use netcdf, only: nf90_var_par_access,nf90_netcdf4 use netcdf, only: nf90_inquire_variable use netcdf, only: nf90_inq_varid use mod_fv3_lola, only: fv3_h_to_ll,fv3_h_to_ll_ens @@ -3051,6 +3070,9 @@ subroutine gsi_fv3ncdf_readuv_v1(grd_uv,ges_u,ges_v,fv3filenamegin,ensgrid) integer(i_kind) nxcase,nycase integer(i_kind) us_countloc(3),us_startloc(3) integer(i_kind) vw_countloc(3),vw_startloc(3) + integer(i_kind):: iworld,iworld_group,nread,mpi_comm_read,i,ierror + integer(i_kind),dimension(npe):: members,members_read,mype_read_rank + logical:: procuse allocate (worksub(2,grd_uv%lat2,grd_uv%lon2,grd_uv%nsig)) mm1=mype+1 @@ -3067,11 +3089,33 @@ subroutine gsi_fv3ncdf_readuv_v1(grd_uv,ges_u,ges_v,fv3filenamegin,ensgrid) kend=grd_uv%kend_loc allocate (us2d(nxcase,nycase+1),vw2d(nxcase+1,nycase)) allocate (uorv2d(nxcase,nycase)) + procuse = .false. + members=-1 + members_read=-1 + if (kbgn<=kend) then + procuse = .true. + members(mm1) = mype + endif + + call mpi_allreduce(members,members_read,npe,mpi_integer,mpi_max,mpi_comm_world,ierror) + + nread=0 + mype_read_rank=-1 + do i=1,npe + if (members_read(i) >= 0) then + nread=nread+1 + mype_read_rank(nread) = members_read(i) + endif + enddo + + call setcomm(iworld,iworld_group,nread,mype_read_rank,mpi_comm_read,ierror) + + if (procuse) then + filenamein=fv3filenamegin%dynvars - iret=nf90_open(filenamein,nf90_nowrite,gfile_loc,comm=mpi_comm_world,info=MPI_INFO_NULL) !clt + iret=nf90_open(filenamein,ior(nf90_netcdf4,ior(nf90_nowrite,nf90_mpiio)),gfile_loc,comm=mpi_comm_read,info=MPI_INFO_NULL) !clt if(iret/=nf90_noerr) then write(6,*)' gsi_fv3ncdf_read_v1: problem opening ',trim(filenamein),gfile_loc,', Status = ',iret - call flush(6) call stop2(333) endif @@ -3080,7 +3124,6 @@ subroutine gsi_fv3ncdf_readuv_v1(grd_uv,ges_u,ges_v,fv3filenamegin,ensgrid) filenamein2=fv3filenamegin%dynvars if(trim(filenamein) /= trim(filenamein2)) then write(6,*)'filenamein and filenamein2 are not the same as expected, stop' - call flush(6) call stop2(333) endif ilev=grd_uv%lnames(1,ilevtot) @@ -3099,9 +3142,9 @@ subroutine gsi_fv3ncdf_readuv_v1(grd_uv,ges_u,ges_v,fv3filenamegin,ensgrid) ! transfor to earth u/v, interpolate to analysis grid, reverse vertical order - iret=nf90_inq_varid(gfile_loc,trim(adjustl("u_s")),var_id) - - iret=nf90_get_var(gfile_loc,var_id,us2d,start=us_startloc,count=us_countloc) + call check(nf90_inq_varid(gfile_loc,trim(adjustl("u_s")),var_id)) + + call check(nf90_get_var(gfile_loc,var_id,us2d,start=us_startloc,count=us_countloc)) iret=nf90_inq_varid(gfile_loc,trim(adjustl("v_w")),var_id) iret=nf90_get_var(gfile_loc,var_id,vw2d,start=vw_startloc,count=vw_countloc) do j=1,ny @@ -3123,10 +3166,11 @@ subroutine gsi_fv3ncdf_readuv_v1(grd_uv,ges_u,ges_v,fv3filenamegin,ensgrid) end if enddo ! iilevtoto + iret=nf90_close(gfile_loc) + endif !procuse call general_grid2sub(grd_uv,hwork,worksub) ges_u=worksub(1,:,:,:) ges_v=worksub(2,:,:,:) - iret=nf90_close(gfile_loc) deallocate (us2d,vw2d,worksub) end subroutine gsi_fv3ncdf_readuv_v1 @@ -3160,7 +3204,7 @@ subroutine gsi_fv3ncdf_read_ens_parallel_over_ens(filenamein,fv3filenamegin, & use mpimod, only: mpi_comm_world,mpi_rtype,mype use mpimod, only: MPI_INFO_NULL use netcdf, only: nf90_open,nf90_close,nf90_get_var,nf90_noerr - use netcdf, only: nf90_nowrite,nf90_inquire,nf90_inquire_dimension + use netcdf, only: nf90_nowrite,nf90_mpiio,nf90_inquire,nf90_inquire_dimension use netcdf, only: nf90_inquire_variable use netcdf, only: nf90_inq_varid use gridmod, only: nsig,nlon,nlat @@ -3179,7 +3223,7 @@ subroutine gsi_fv3ncdf_read_ens_parallel_over_ens(filenamein,fv3filenamegin, & character(len=max_varname_length) :: name character(len=max_filename_length), allocatable,dimension(:) :: varname_files - integer(i_kind) nlatcase,nloncase,nxcase,nycase,countloc(3),startloc(3),countloc_tmp(3),startloc_tmp(3) + integer(i_kind) nlatcase,nloncase,nxcase,nycase,countloc(4),startloc(4),countloc_tmp(4),startloc_tmp(4) integer(i_kind) ilev,ilevtot,inative,ivar integer(i_kind) kbgn,kend integer(i_kind) gfile_loc,iret,var_id @@ -3238,15 +3282,13 @@ subroutine gsi_fv3ncdf_read_ens_parallel_over_ens(filenamein,fv3filenamegin, & iret=nf90_open(filename_layout,nf90_nowrite,gfile_loc_layout(nio),comm=mpi_comm_world,info=MPI_INFO_NULL) if(iret/=nf90_noerr) then write(6,*)' gsi_fv3ncdf_read: problem opening ',trim(filename_layout),gfile_loc_layout(nio),', Status = ',iret - call flush(6) call stop2(333) endif enddo else - iret=nf90_open(filenamein,nf90_nowrite,gfile_loc) + iret=nf90_open(filenamein,ior(nf90_nowrite,nf90_mpiio),gfile_loc) if(iret/=nf90_noerr) then write(6,*)' gsi_fv3ncdf_read: problem opening ',trim(filenamein),gfile_loc,', Status = ',iret - call flush(6) call stop2(333) endif endif @@ -3256,8 +3298,8 @@ subroutine gsi_fv3ncdf_read_ens_parallel_over_ens(filenamein,fv3filenamegin, & nz=nsig nzp1=nz+1 inative=nzp1-ilev - startloc=(/1,1,inative/) - countloc=(/nxcase,nycase,1/) + startloc=(/1,1,inative,1/) + countloc=(/nxcase,nycase,1,1/) 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 @@ -3266,19 +3308,19 @@ subroutine gsi_fv3ncdf_read_ens_parallel_over_ens(filenamein,fv3filenamegin, & if(trim(name)=='xaxis_1') nx_phy=len if( nx_phy == nxcase )then allocate(uu2d_tmp(nxcase,nycase)) - countloc_tmp=(/nxcase,nycase,1/) + countloc_tmp=(/nxcase,nycase,1,1/) phy_smaller_domain = .false. else allocate(uu2d_tmp(nxcase-6,nycase-6)) - countloc_tmp=(/nxcase-6,nycase-6,1/) + countloc_tmp=(/nxcase-6,nycase-6,1,1/) phy_smaller_domain = .true. end if - startloc_tmp=(/1,1,ilev/) + startloc_tmp=(/1,1,ilev,1/) end if if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 - countloc=(/nxcase,ny_layout_len(nio),1/) + countloc=(/nxcase,ny_layout_len(nio),1,1/) allocate(uu2d_layout(nxcase,ny_layout_len(nio))) iret=nf90_inq_varid(gfile_loc_layout(nio),trim(adjustl(varname)),var_id) iret=nf90_get_var(gfile_loc_layout(nio),var_id,uu2d_layout,start=startloc,count=countloc) @@ -3408,7 +3450,7 @@ subroutine gsi_fv3ncdf_readuv_ens_parallel_over_ens(ges_u,ges_v,fv3filenamegin,i integer(i_kind) u_grd_VarId,v_grd_VarId integer(i_kind) nlatcase,nloncase integer(i_kind) nxcase,nycase - integer(i_kind) u_countloc(3),u_startloc(3),v_countloc(3),v_startloc(3) + integer(i_kind) u_countloc(4),u_startloc(4),v_countloc(4),v_startloc(4) integer(i_kind) inative,ilev,ilevtot integer(i_kind) kbgn,kend @@ -3442,7 +3484,6 @@ subroutine gsi_fv3ncdf_readuv_ens_parallel_over_ens(ges_u,ges_v,fv3filenamegin,i iret=nf90_open(filename_layout,nf90_nowrite,gfile_loc_layout(nio),comm=mpi_comm_world,info=MPI_INFO_NULL) if(iret/=nf90_noerr) then write(6,*)'problem opening ',trim(filename_layout),gfile_loc_layout(nio),', Status = ',iret - call flush(6) call stop2(333) endif enddo @@ -3450,7 +3491,6 @@ subroutine gsi_fv3ncdf_readuv_ens_parallel_over_ens(ges_u,ges_v,fv3filenamegin,i iret=nf90_open(filenamein,nf90_nowrite,gfile_loc) if(iret/=nf90_noerr) then write(6,*)' problem opening ',trim(filenamein),', Status = ',iret - call flush(6) call stop2(333) endif endif @@ -3459,14 +3499,14 @@ subroutine gsi_fv3ncdf_readuv_ens_parallel_over_ens(ges_u,ges_v,fv3filenamegin,i nz=nsig nzp1=nz+1 inative=nzp1-ilev - u_countloc=(/nxcase,nycase+1,1/) - v_countloc=(/nxcase+1,nycase,1/) - u_startloc=(/1,1,inative/) - v_startloc=(/1,1,inative/) + u_countloc=(/nxcase,nycase+1,1,1/) + v_countloc=(/nxcase+1,nycase,1,1/) + u_startloc=(/1,1,inative,1/) + v_startloc=(/1,1,inative,1/) if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 - u_countloc=(/nxcase,ny_layout_len(nio)+1,1/) + u_countloc=(/nxcase,ny_layout_len(nio)+1,1,1/) allocate(u2d_layout(nxcase,ny_layout_len(nio)+1)) call check( nf90_inq_varid(gfile_loc_layout(nio),'u',u_grd_VarId) ) iret=nf90_get_var(gfile_loc_layout(nio),u_grd_VarId,u2d_layout,start=u_startloc,count=u_countloc) @@ -3474,7 +3514,7 @@ subroutine gsi_fv3ncdf_readuv_ens_parallel_over_ens(ges_u,ges_v,fv3filenamegin,i if(nio==fv3_io_layout_y-1) u2d(:,ny_layout_e(nio)+1)=u2d_layout(:,ny_layout_len(nio)+1) deallocate(u2d_layout) - v_countloc=(/nxcase+1,ny_layout_len(nio),1/) + v_countloc=(/nxcase+1,ny_layout_len(nio),1,1/) allocate(v2d_layout(nxcase+1,ny_layout_len(nio))) call check( nf90_inq_varid(gfile_loc_layout(nio),'v',v_grd_VarId) ) iret=nf90_get_var(gfile_loc_layout(nio),v_grd_VarId,v2d_layout,start=v_startloc,count=v_countloc) @@ -3965,7 +4005,7 @@ subroutine gsi_fv3ncdf_writeuv(grd_uv,ges_u,ges_v,add_saved,fv3filenamegin) use mod_fv3_lola, only: fv3_ll_to_h,fv3_h_to_ll, & fv3uv2earth,earthuv2fv3 use netcdf, only: nf90_open,nf90_close,nf90_noerr - use netcdf, only: nf90_write,nf90_inq_varid + use netcdf, only: nf90_write,nf90_mpiio,nf90_inq_varid,nf90_var_par_access,nf90_collective use netcdf, only: nf90_put_var,nf90_get_var use general_sub2grid_mod, only: sub2grid_info,general_sub2grid @@ -3984,11 +4024,11 @@ subroutine gsi_fv3ncdf_writeuv(grd_uv,ges_u,ges_v,add_saved,fv3filenamegin) integer(i_kind) inative,ilev,ilevtot integer(i_kind) nlatcase,nloncase integer(i_kind) nxcase,nycase - integer(i_kind) u_countloc(3),u_startloc(3),v_countloc(3),v_startloc(3) + integer(i_kind) u_countloc(4),u_startloc(4),v_countloc(4),v_startloc(4) character(:),allocatable:: filenamein ,varname real(r_kind),allocatable,dimension(:,:,:,:):: worksub real(r_kind),allocatable,dimension(:,:):: work_au,work_av - real(r_kind),allocatable,dimension(:,:):: work_bu,work_bv + real(r_kind),allocatable,dimension(:,:,:):: work_bu,work_bv real(r_kind),allocatable,dimension(:,:):: u2d,v2d,workau2,workav2 real(r_kind),allocatable,dimension(:,:):: workbu2,workbv2 @@ -3997,10 +4037,12 @@ subroutine gsi_fv3ncdf_writeuv(grd_uv,ges_u,ges_v,add_saved,fv3filenamegin) logical:: procuse ! for fv3_io_layout_y > 1 - real(r_kind),allocatable,dimension(:,:):: u2d_layout,v2d_layout + real(r_kind),allocatable,dimension(:,:,:):: u2d_layout,v2d_layout integer(i_kind) :: nio integer(i_kind),allocatable :: gfile_loc_layout(:) character(len=180) :: filename_layout + integer(i_kind):: kend_native,kbgn_native + integer(i_kind):: istat mm1=mype+1 @@ -4012,8 +4054,6 @@ subroutine gsi_fv3ncdf_writeuv(grd_uv,ges_u,ges_v,add_saved,fv3filenamegin) kend=grd_uv%kend_loc allocate( u2d(nlon_regional,nlat_regional+1)) allocate( v2d(nlon_regional+1,nlat_regional)) - allocate( work_bu(nlon_regional,nlat_regional+1)) - allocate( work_bv(nlon_regional+1,nlat_regional)) allocate (worksub(2,grd_uv%lat2,grd_uv%lon2,grd_uv%nsig)) allocate( work_au(nlatcase,nloncase),work_av(nlatcase,nloncase)) do k=1,grd_uv%nsig @@ -4054,59 +4094,70 @@ subroutine gsi_fv3ncdf_writeuv(grd_uv,ges_u,ges_v,add_saved,fv3filenamegin) allocate(gfile_loc_layout(0:fv3_io_layout_y-1)) do nio=0,fv3_io_layout_y-1 write(filename_layout,'(a,a,I4.4)') trim(filenamein),".",nio - call check( nf90_open(filename_layout,nf90_write,gfile_loc_layout(nio),comm=mpi_comm_read,info=MPI_INFO_NULL) ) + call check( nf90_open(filename_layout,ior(nf90_write, nf90_mpiio),gfile_loc_layout(nio),comm=mpi_comm_read,info=MPI_INFO_NULL) ) enddo gfile_loc=gfile_loc_layout(0) else - call check( nf90_open(filenamein,nf90_write,gfile_loc,comm=mpi_comm_read,info=MPI_INFO_NULL) ) + call check( nf90_open(filenamein,ior(nf90_write, nf90_mpiio),gfile_loc,comm=mpi_comm_read,info=MPI_INFO_NULL) ) + endif + nz=grd_uv%nsig + nzp1=nz+1 + kend_native=nzp1-grd_uv%lnames(1,kbgn) + kbgn_native=nzp1-grd_uv%lnames(1,kend) + allocate( work_bu(nlon_regional,nlat_regional+1,kbgn_native:kend_native)) + allocate( work_bv(nlon_regional+1,nlat_regional,kbgn_native:kend_native)) + u_startloc=(/1,1,kbgn_native,1/) + u_countloc=(/nxcase,nycase+1,kend_native-kbgn_native+1,1/) + v_startloc=(/1,1,kbgn_native,1/) + v_countloc=(/nxcase+1,nycase,kend_native-kbgn_native+1,1/) + if(fv3_io_layout_y > 1) then + do nio=0,fv3_io_layout_y-1 + allocate(u2d_layout(nxcase,ny_layout_len(nio)+1,kend_native-kbgn_native+1)) + u_countloc=(/nxcase,ny_layout_len(nio)+1,kend_native-kbgn_native+1,1/) + call check( nf90_get_var(gfile_loc_layout(nio),ugrd_VarId,u2d_layout,start=u_startloc,count=u_countloc) ) + work_bu(:,ny_layout_b(nio):ny_layout_e(nio),:)=u2d_layout(:,1:ny_layout_len(nio),:) + if(nio==fv3_io_layout_y-1) work_bu(:,ny_layout_e(nio)+1,:)=u2d_layout(:,ny_layout_len(nio)+1,:) + deallocate(u2d_layout) + + allocate(v2d_layout(nxcase+1,ny_layout_len(nio),kend_native-kbgn_native+1)) + v_countloc=(/nxcase+1,ny_layout_len(nio),kend_native-kbgn_native+1,1/) + call check( nf90_get_var(gfile_loc_layout(nio),vgrd_VarId,v2d_layout,start=v_startloc,count=v_countloc) ) + work_bv(:,ny_layout_b(nio):ny_layout_e(nio),:)=v2d_layout + deallocate(v2d_layout) + enddo + else + call check( nf90_inq_varid(gfile_loc,'u',ugrd_VarId) ) + call check( nf90_inq_varid(gfile_loc,'v',vgrd_VarId) ) + call check( nf90_var_par_access(gfile_loc, ugrd_VarId, nf90_collective)) + call check( nf90_var_par_access(gfile_loc, vgrd_VarId, nf90_collective)) + call check( nf90_get_var(gfile_loc,ugrd_VarId,work_bu,start=u_startloc,count=u_countloc) ) + call check( nf90_get_var(gfile_loc,vgrd_VarId,work_bv,start=v_startloc,count=v_countloc) ) endif + + do ilevtot=kbgn,kend varname=grd_uv%names(1,ilevtot) ilev=grd_uv%lnames(1,ilevtot) - nz=grd_uv%nsig - nzp1=nz+1 inative=nzp1-ilev - u_countloc=(/nxcase,nycase+1,1/) - v_countloc=(/nxcase+1,nycase,1/) - u_startloc=(/1,1,inative/) - v_startloc=(/1,1,inative/) work_au=hwork(1,:,:,ilevtot) work_av=hwork(2,:,:,ilevtot) - call check( nf90_inq_varid(gfile_loc,'u',ugrd_VarId) ) - call check( nf90_inq_varid(gfile_loc,'v',vgrd_VarId) ) if(add_saved)then allocate( workau2(nlatcase,nloncase),workav2(nlatcase,nloncase)) allocate( workbu2(nlon_regional,nlat_regional+1)) allocate( workbv2(nlon_regional+1,nlat_regional)) !!!!!!!! readin work_b !!!!!!!!!!!!!!!! - if(fv3_io_layout_y > 1) then - do nio=0,fv3_io_layout_y-1 - allocate(u2d_layout(nxcase,ny_layout_len(nio)+1)) - u_countloc=(/nxcase,ny_layout_len(nio)+1,1/) - call check( nf90_get_var(gfile_loc_layout(nio),ugrd_VarId,u2d_layout,start=u_startloc,count=u_countloc) ) - work_bu(:,ny_layout_b(nio):ny_layout_e(nio))=u2d_layout(:,1:ny_layout_len(nio)) - if(nio==fv3_io_layout_y-1) work_bu(:,ny_layout_e(nio)+1)=u2d_layout(:,ny_layout_len(nio)+1) - deallocate(u2d_layout) - - allocate(v2d_layout(nxcase+1,ny_layout_len(nio))) - v_countloc=(/nxcase+1,ny_layout_len(nio),1/) - call check( nf90_get_var(gfile_loc_layout(nio),vgrd_VarId,v2d_layout,start=v_startloc,count=v_countloc) ) - work_bv(:,ny_layout_b(nio):ny_layout_e(nio))=v2d_layout - deallocate(v2d_layout) - enddo - else - call check( nf90_get_var(gfile_loc,ugrd_VarId,work_bu,start=u_startloc,count=u_countloc) ) - call check( nf90_get_var(gfile_loc,vgrd_VarId,work_bv,start=v_startloc,count=v_countloc) ) - endif + +!clt for fv3_io_layout<=1 now the nf90_get_var has been moved outside of this do loop +!to avoid failure on hercules when L_MPI_EXTRA_FILESYSTEM=1 if(.not.grid_reverse_flag) then - call reverse_grid_r_uv(work_bu,nlon_regional,nlat_regional+1,1) - call reverse_grid_r_uv(work_bv,nlon_regional+1,nlat_regional,1) + call reverse_grid_r_uv(work_bu(:,:,inative),nlon_regional,nlat_regional+1,1) + call reverse_grid_r_uv(work_bv(:,:,inative),nlon_regional+1,nlat_regional,1) endif - call fv3uv2earth(work_bu,work_bv,nlon_regional,nlat_regional,u2d,v2d) + call fv3uv2earth(work_bu(:,:,inative),work_bv(:,:,inative),nlon_regional,nlat_regional,u2d,v2d) call fv3_h_to_ll(u2d,workau2,nlon_regional,nlat_regional,nloncase,nlatcase,.true.) call fv3_h_to_ll(v2d,workav2,nlon_regional,nlat_regional,nloncase,nlatcase,.true.) !!!!!!!! find analysis_inc: work_a !!!!!!!!!!!!!!!! @@ -4116,38 +4167,38 @@ subroutine gsi_fv3ncdf_writeuv(grd_uv,ges_u,ges_v,add_saved,fv3filenamegin) call fv3_ll_to_h(work_av(:,:),v2d,nloncase,nlatcase,nlon_regional,nlat_regional,.true.) call earthuv2fv3(u2d,v2d,nlon_regional,nlat_regional,workbu2,workbv2) !!!!!!!! add analysis_inc to readin work_b !!!!!!!!!!!!!!!! - work_bu(:,:)=work_bu(:,:)+workbu2(:,:) - work_bv(:,:)=work_bv(:,:)+workbv2(:,:) + work_bu(:,:,inative)=work_bu(:,:,inative)+workbu2(:,:) + work_bv(:,:,inative)=work_bv(:,:,inative)+workbv2(:,:) deallocate(workau2,workbu2,workav2,workbv2) else call fv3_ll_to_h(work_au(:,:),u2d,nloncase,nlatcase,nlon_regional,nlat_regional,.true.) call fv3_ll_to_h(work_av(:,:),v2d,nloncase,nlatcase,nlon_regional,nlat_regional,.true.) - call earthuv2fv3(u2d,v2d,nlon_regional,nlat_regional,work_bu(:,:),work_bv(:,:)) + call earthuv2fv3(u2d,v2d,nlon_regional,nlat_regional,work_bu(:,:,inative),work_bv(:,:,inative)) endif if(.not.grid_reverse_flag) then - call reverse_grid_r_uv(work_bu,nlon_regional,nlat_regional+1,1) - call reverse_grid_r_uv(work_bv,nlon_regional+1,nlat_regional,1) + call reverse_grid_r_uv(work_bu(:,:,inative),nlon_regional,nlat_regional+1,1) + call reverse_grid_r_uv(work_bv(:,:,inative),nlon_regional+1,nlat_regional,1) endif + enddo !ilevltot - if(fv3_io_layout_y > 1) then + if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 - allocate(u2d_layout(nxcase,ny_layout_len(nio)+1)) - u_countloc=(/nxcase,ny_layout_len(nio)+1,1/) - u2d_layout=work_bu(:,ny_layout_b(nio):ny_layout_e(nio)+1) + allocate(u2d_layout(nxcase,ny_layout_len(nio)+1,kend_native-kbgn_native+1)) + u_countloc=(/nxcase,ny_layout_len(nio)+1,kend_native-kbgn_native+1,1/) + u2d_layout=work_bu(:,ny_layout_b(nio):ny_layout_e(nio)+1,:) call check( nf90_put_var(gfile_loc_layout(nio),ugrd_VarId,u2d_layout,start=u_startloc,count=u_countloc) ) deallocate(u2d_layout) - allocate(v2d_layout(nxcase+1,ny_layout_len(nio))) - v_countloc=(/nxcase+1,ny_layout_len(nio),1/) - v2d_layout=work_bv(:,ny_layout_b(nio):ny_layout_e(nio)) + allocate(v2d_layout(nxcase+1,ny_layout_len(nio),kend_native-kbgn_native+1)) + v_countloc=(/nxcase+1,ny_layout_len(nio),kend_native-kbgn_native+1,1/) + v2d_layout=work_bv(:,ny_layout_b(nio):ny_layout_e(nio),:) call check( nf90_put_var(gfile_loc_layout(nio),vgrd_VarId,v2d_layout,start=v_startloc,count=v_countloc) ) deallocate(v2d_layout) enddo - else + else call check( nf90_put_var(gfile_loc,ugrd_VarId,work_bu,start=u_startloc,count=u_countloc) ) call check( nf90_put_var(gfile_loc,vgrd_VarId,work_bv,start=v_startloc,count=v_countloc) ) - endif - enddo !ilevltot + endif if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 @@ -4157,11 +4208,12 @@ subroutine gsi_fv3ncdf_writeuv(grd_uv,ges_u,ges_v,add_saved,fv3filenamegin) else call check( nf90_close(gfile_loc) ) endif + deallocate(work_bu,work_bv) endif call mpi_barrier(mpi_comm_world,ierror) - deallocate(work_bu,work_bv,u2d,v2d) + deallocate(u2d,v2d) deallocate(work_au,work_av) end subroutine gsi_fv3ncdf_writeuv @@ -4193,12 +4245,12 @@ subroutine gsi_fv3ncdf_writeuv_v1(grd_uv,ges_u,ges_v,add_saved,fv3filenamegin) !$$$ end documentation block use constants, only: half,zero - use mpimod, only: mpi_rtype,mpi_comm_world,mype,mpi_info_null + use mpimod, only: npe, setcomm,mpi_integer,mpi_max,mpi_rtype,mpi_comm_world,mype,mpi_info_null use gridmod, only: nlon_regional,nlat_regional use mod_fv3_lola, only: fv3_ll_to_h,fv3_h_to_ll, & fv3uv2earth,earthuv2fv3 use netcdf, only: nf90_open,nf90_close,nf90_noerr - use netcdf, only: nf90_write,nf90_inq_varid + use netcdf, only: nf90_write, nf90_mpiio,nf90_inq_varid,nf90_var_par_access,nf90_collective use netcdf, only: nf90_put_var,nf90_get_var use general_sub2grid_mod, only: sub2grid_info,general_sub2grid implicit none @@ -4220,14 +4272,20 @@ subroutine gsi_fv3ncdf_writeuv_v1(grd_uv,ges_u,ges_v,add_saved,fv3filenamegin) integer(i_kind) inative,ilev,ilevtot real(r_kind),allocatable,dimension(:,:,:,:):: worksub real(r_kind),allocatable,dimension(:,:):: work_au,work_av - real(r_kind),allocatable,dimension(:,:):: work_bu_s,work_bv_s - real(r_kind),allocatable,dimension(:,:):: work_bu_w,work_bv_w + real(r_kind),allocatable,dimension(:,:,:):: work_bu_s,work_bv_s + real(r_kind),allocatable,dimension(:,:,:):: work_bu_w,work_bv_w real(r_kind),allocatable,dimension(:,:):: u2d,v2d,workau2,workav2 real(r_kind),allocatable,dimension(:,:):: workbu_s2,workbv_s2 real(r_kind),allocatable,dimension(:,:):: workbu_w2,workbv_w2 integer(i_kind) nlatcase,nloncase,nxcase,nycase - integer(i_kind) uw_countloc(3),us_countloc(3),uw_startloc(3),us_startloc(3) - integer(i_kind) vw_countloc(3),vs_countloc(3),vw_startloc(3),vs_startloc(3) + integer(i_kind) uw_countloc(4),us_countloc(4),uw_startloc(4),us_startloc(4) + integer(i_kind) vw_countloc(4),vs_countloc(4),vw_startloc(4),vs_startloc(4) + integer(i_kind):: kend_native,kbgn_native,kdim_native + + + integer(i_kind):: iworld,iworld_group,nread,mpi_comm_read,ierror + integer(i_kind),dimension(npe):: members,members_read,mype_read_rank + logical:: procuse mm1=mype+1 nloncase=grd_uv%nlon @@ -4249,61 +4307,96 @@ subroutine gsi_fv3ncdf_writeuv_v1(grd_uv,ges_u,ges_v,add_saved,fv3filenamegin) allocate( u2d(nlon_regional,nlat_regional)) allocate( v2d(nlon_regional,nlat_regional)) - allocate( work_bu_s(nlon_regional,nlat_regional+1)) - allocate( work_bv_s(nlon_regional,nlat_regional+1)) - allocate( work_bu_w(nlon_regional+1,nlat_regional)) - allocate( work_bv_w(nlon_regional+1,nlat_regional)) allocate( work_au(nlatcase,nloncase),work_av(nlatcase,nloncase)) + if(add_saved) allocate( workau2(nlatcase,nloncase),workav2(nlatcase,nloncase)) - allocate( workbu_w2(nlon_regional+1,nlat_regional)) - allocate( workbv_w2(nlon_regional+1,nlat_regional)) - allocate( workbu_s2(nlon_regional,nlat_regional+1)) - allocate( workbv_s2(nlon_regional,nlat_regional+1)) + allocate( workbu_w2(nlon_regional+1,nlat_regional)) + allocate( workbv_w2(nlon_regional+1,nlat_regional)) + allocate( workbu_s2(nlon_regional,nlat_regional+1)) + allocate( workbv_s2(nlon_regional,nlat_regional+1)) filenamein=fv3filenamegin%dynvars - call check( nf90_open(filenamein,nf90_write,gfile_loc,comm=mpi_comm_world,info=MPI_INFO_NULL) ) - do ilevtot=kbgn,kend - varname=grd_uv%names(1,ilevtot) - ilev=grd_uv%lnames(1,ilevtot) - nz=grd_uv%nsig - nzp1=nz+1 - inative=nzp1-ilev + procuse = .false. + members=-1 + members_read=-1 + if (kbgn<=kend) then + procuse = .true. + members(mm1) = mype + endif + + call mpi_allreduce(members,members_read,npe,mpi_integer,mpi_max,mpi_comm_world,ierror) + + nread=0 + mype_read_rank=-1 + do i=1,npe + if (members_read(i) >= 0) then + nread=nread+1 + mype_read_rank(nread) = members_read(i) + endif + enddo + + call setcomm(iworld,iworld_group,nread,mype_read_rank,mpi_comm_read,ierror) + + if (procuse) then + - uw_countloc= (/nlon_regional+1,nlat_regional,1/) - us_countloc= (/nlon_regional,nlat_regional+1,1/) - vw_countloc= (/nlon_regional+1,nlat_regional,1/) - vs_countloc= (/nlon_regional,nlat_regional+1,1/) + + + call check( nf90_open(filenamein,ior(nf90_write, nf90_mpiio),gfile_loc,comm=mpi_comm_read,info=MPI_INFO_NULL) ) + + call check( nf90_inq_varid(gfile_loc,'u_s',u_sgrd_VarId) ) + call check( nf90_var_par_access(gfile_loc, u_sgrd_VarId, nf90_collective)) + call check( nf90_inq_varid(gfile_loc,'u_w',u_wgrd_VarId) ) + call check( nf90_var_par_access(gfile_loc, u_wgrd_VarId, nf90_collective)) + call check( nf90_inq_varid(gfile_loc,'v_s',v_sgrd_VarId) ) + call check( nf90_var_par_access(gfile_loc, v_sgrd_VarId, nf90_collective)) + call check( nf90_inq_varid(gfile_loc,'v_w',v_wgrd_VarId) ) + call check( nf90_var_par_access(gfile_loc, v_wgrd_VarId, nf90_collective)) + nz=grd_uv%nsig + nzp1=nz+1 + kend_native=nzp1-grd_uv%lnames(1,kbgn) + kbgn_native=nzp1-grd_uv%lnames(1,kend) + kdim_native=kend_native-kbgn_native+1 + + uw_countloc= (/nlon_regional+1,nlat_regional,kdim_native,1/) + us_countloc= (/nlon_regional,nlat_regional+1,kdim_native,1/) + vw_countloc= (/nlon_regional+1,nlat_regional,kdim_native,1/) + vs_countloc= (/nlon_regional,nlat_regional+1,kdim_native,1/) - uw_startloc=(/1,1,inative+1/) - us_startloc=(/1,1,inative+1/) - vw_startloc=(/1,1,inative+1/) - vs_startloc=(/1,1,inative+1/) + uw_startloc=(/1,1,kbgn_native+1,1/) !In the coldstart files, there is an extra top level + us_startloc=(/1,1,kbgn_native+1,1/) + vw_startloc=(/1,1,kbgn_native+1,1/) + vs_startloc=(/1,1,kbgn_native+1,1/) + allocate( work_bu_s(nlon_regional,nlat_regional+1,kbgn_native:kend_native)) + allocate( work_bv_s(nlon_regional,nlat_regional+1,kbgn_native:kend_native)) + allocate( work_bu_w(nlon_regional+1,nlat_regional,kbgn_native:kend_native)) + allocate( work_bv_w(nlon_regional+1,nlat_regional,kbgn_native:kend_native)) +!!!!!!!! readin work_b !!!!!!!!!!!!!!!! + call check( nf90_get_var(gfile_loc,u_sgrd_VarId,work_bu_s,start=us_startloc,count=us_countloc) ) + call check( nf90_get_var(gfile_loc,u_wgrd_VarId,work_bu_w,start=uw_startloc,count=uw_countloc) ) + call check( nf90_get_var(gfile_loc,v_sgrd_VarId,work_bv_s,start=vs_startloc,count=vs_countloc) ) + call check( nf90_get_var(gfile_loc,v_wgrd_VarId,work_bv_w,start=vw_startloc,count=vw_countloc) ) + do ilevtot=kbgn,kend + varname=grd_uv%names(1,ilevtot) + ilev=grd_uv%lnames(1,ilevtot) + inative=nzp1-ilev work_au=hwork(1,:,:,ilevtot) work_av=hwork(2,:,:,ilevtot) - call check( nf90_inq_varid(gfile_loc,'u_s',u_sgrd_VarId) ) - call check( nf90_inq_varid(gfile_loc,'u_w',u_wgrd_VarId) ) - call check( nf90_inq_varid(gfile_loc,'v_s',v_sgrd_VarId) ) - call check( nf90_inq_varid(gfile_loc,'v_w',v_wgrd_VarId) ) -!!!!!!!! readin work_b !!!!!!!!!!!!!!!! - call check( nf90_get_var(gfile_loc,u_sgrd_VarId,work_bu_s,start=us_startloc,count=us_countloc) ) - call check( nf90_get_var(gfile_loc,u_wgrd_VarId,work_bu_w,start=uw_startloc,count=uw_countloc) ) - call check( nf90_get_var(gfile_loc,v_sgrd_VarId,work_bv_s,start=vs_startloc,count=vs_countloc) ) - call check( nf90_get_var(gfile_loc,v_wgrd_VarId,work_bv_w,start=vw_startloc,count=vw_countloc) ) if(add_saved)then do j=1,nlat_regional - u2d(:,j)=half * (work_bu_s(:,j)+ work_bu_s(:,j+1)) + u2d(:,j)=half * (work_bu_s(:,j,inative)+ work_bu_s(:,j+1,inative)) enddo do i=1,nlon_regional - v2d(i,:)=half*(work_bv_w(i,:)+work_bv_w(i+1,:)) + v2d(i,:)=half*(work_bv_w(i,:,inative)+work_bv_w(i+1,:,inative)) enddo call fv3_h_to_ll(u2d,workau2,nlon_regional,nlat_regional,nloncase,nlatcase,grid_reverse_flag) call fv3_h_to_ll(v2d,workav2,nlon_regional,nlat_regional,nloncase,nlatcase,grid_reverse_flag) @@ -4333,44 +4426,46 @@ subroutine gsi_fv3ncdf_writeuv_v1(grd_uv,ges_u,ges_v,add_saved,fv3filenamegin) - work_bu_w(:,:)=work_bu_w(:,:)+workbu_w2(:,:) - work_bu_s(:,:)=work_bu_s(:,:)+workbu_s2(:,:) - work_bv_w(:,:)=work_bv_w(:,:)+workbv_w2(:,:) - work_bv_s(:,:)=work_bv_s(:,:)+workbv_s2(:,:) + work_bu_w(:,:,inative)=work_bu_w(:,:,inative)+workbu_w2(:,:) + work_bu_s(:,:,inative)=work_bu_s(:,:,inative)+workbu_s2(:,:) + work_bv_w(:,:,inative)=work_bv_w(:,:,inative)+workbv_w2(:,:) + work_bv_s(:,:,inative)=work_bv_s(:,:,inative)+workbv_s2(:,:) else call fv3_ll_to_h(work_au(:,:),u2d,nloncase,nlatcase,nlon_regional,nlat_regional,grid_reverse_flag) call fv3_ll_to_h(work_av(:,:),v2d,nloncase,nlatcase,nlon_regional,nlat_regional,grid_reverse_flag) do i=2,nlon_regional - work_bu_w(i,:)=half*(u2d(i-1,:)+u2d(i,:)) - work_bv_w(i,:)=half*(v2d(i-1,:)+v2d(i,:)) + work_bu_w(i,:,inative)=half*(u2d(i-1,:)+u2d(i,:)) + work_bv_w(i,:,inative)=half*(v2d(i-1,:)+v2d(i,:)) enddo - work_bu_w(1,:)=u2d(1,:) - work_bv_w(1,:)=v2d(1,:) - work_bu_w(nlon_regional+1,:)=u2d(nlon_regional,:) - work_bv_w(nlon_regional+1,:)=v2d(nlon_regional,:) + work_bu_w(1,:,inative)=u2d(1,:) + work_bv_w(1,:,inative)=v2d(1,:) + work_bu_w(nlon_regional+1,:,inative)=u2d(nlon_regional,:) + work_bv_w(nlon_regional+1,:,inative)=v2d(nlon_regional,:) do j=2,nlat_regional - work_bu_s(:,j)=half*(u2d(:,j-1)+u2d(:,j)) - work_bv_s(:,j)=half*(v2d(:,j-1)+v2d(:,j)) + work_bu_s(:,j,inative)=half*(u2d(:,j-1)+u2d(:,j)) + work_bv_s(:,j,inative)=half*(v2d(:,j-1)+v2d(:,j)) enddo - work_bu_s(:,1)=u2d(:,1) - work_bv_s(:,1)=v2d(:,1) - work_bu_s(:,nlat_regional+1)=u2d(:,nlat_regional) - work_bv_s(:,nlat_regional+1)=v2d(:,nlat_regional) + work_bu_s(:,1,inative)=u2d(:,1) + work_bv_s(:,1,inative)=v2d(:,1) + work_bu_s(:,nlat_regional+1,inative)=u2d(:,nlat_regional) + work_bv_s(:,nlat_regional+1,inative)=v2d(:,nlat_regional) endif - - call check( nf90_put_var(gfile_loc,u_wgrd_VarId,work_bu_w,start=uw_startloc,count=uw_countloc) ) - call check( nf90_put_var(gfile_loc,u_sgrd_VarId,work_bu_s,start=us_startloc,count=us_countloc) ) - call check( nf90_put_var(gfile_loc,v_wgrd_VarId,work_bv_w,start=vw_startloc,count=vw_countloc) ) - call check( nf90_put_var(gfile_loc,v_sgrd_VarId,work_bv_s,start=vs_startloc,count=vs_countloc) ) enddo ! + + call check( nf90_put_var(gfile_loc,u_wgrd_VarId,work_bu_w,start=uw_startloc,count=uw_countloc) ) + call check( nf90_put_var(gfile_loc,u_sgrd_VarId,work_bu_s,start=us_startloc,count=us_countloc) ) + call check( nf90_put_var(gfile_loc,v_wgrd_VarId,work_bv_w,start=vw_startloc,count=vw_countloc) ) + call check( nf90_put_var(gfile_loc,v_sgrd_VarId,work_bv_s,start=vs_startloc,count=vs_countloc) ) call check( nf90_close(gfile_loc) ) deallocate(work_bu_w,work_bv_w) deallocate(work_bu_s,work_bv_s) + endif !procuse + deallocate(work_au,work_av,u2d,v2d) if(add_saved) deallocate(workau2,workav2) if (allocated(workbu_w2)) then @@ -4541,8 +4636,8 @@ subroutine gsi_fv3ncdf_write(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3file use mod_fv3_lola, only: fv3_ll_to_h use mod_fv3_lola, only: fv3_h_to_ll use netcdf, only: nf90_open,nf90_close - use netcdf, only: nf90_write,nf90_inq_varid - use netcdf, only: nf90_put_var,nf90_get_var + use netcdf, only: nf90_write,nf90_netcdf4, nf90_mpiio,nf90_inq_varid + use netcdf, only: nf90_put_var,nf90_get_var,nf90_independent,nf90_var_par_access use netcdf, only: nf90_inquire_dimension use gsi_bundlemod, only: gsi_bundle use general_sub2grid_mod, only: sub2grid_info,general_sub2grid @@ -4558,8 +4653,8 @@ subroutine gsi_fv3ncdf_write(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3file character(len=max_filename_length) :: filenamein2 character(len=max_varname_length) :: varname,vgsiname,name - integer(i_kind) nlatcase,nloncase,nxcase,nycase,countloc(3),startloc(3) - integer(i_kind) countloc_tmp(3),startloc_tmp(3) + integer(i_kind) nlatcase,nloncase,nxcase,nycase,countloc(4),startloc(4) + integer(i_kind) countloc_tmp(4),startloc_tmp(4) integer(i_kind) kbgn,kend integer(i_kind) inative,ilev,ilevtot integer(i_kind) :: VarId,gfile_loc @@ -4621,11 +4716,11 @@ subroutine gsi_fv3ncdf_write(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3file allocate(gfile_loc_layout(0:fv3_io_layout_y-1)) do nio=0,fv3_io_layout_y-1 write(filename_layout,'(a,a,I4.4)') trim(filenamein),'.',nio - call check( nf90_open(filename_layout,nf90_write,gfile_loc_layout(nio),comm=mpi_comm_read,info=MPI_INFO_NULL) ) + call check( nf90_open(filename_layout,ior(nf90_netcdf4,ior(nf90_write, nf90_mpiio)),gfile_loc_layout(nio),comm=mpi_comm_read,info=MPI_INFO_NULL) ) enddo gfile_loc=gfile_loc_layout(0) else - call check( nf90_open(filenamein,nf90_write,gfile_loc,comm=mpi_comm_read,info=MPI_INFO_NULL) ) + call check( nf90_open(filenamein,ior(nf90_netcdf4,ior(nf90_write, nf90_mpiio)),gfile_loc,comm=mpi_comm_read,info=MPI_INFO_NULL) ) endif do ilevtot=kbgn,kend @@ -4637,15 +4732,14 @@ subroutine gsi_fv3ncdf_write(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3file call getfv3lamfilevname(vgsiname,fv3filenamegin,filenamein2,varname) if(trim(filenamein) /= trim(filenamein2)) then write(6,*)'filenamein and filenamein2 are not the same as expected, stop' - call flush(6) call stop2(333) endif ilev=grd_ionouv%lnames(1,ilevtot) nz=grd_ionouv%nsig nzp1=nz+1 inative=nzp1-ilev - countloc=(/nxcase,nycase,1/) - startloc=(/1,1,inative/) + countloc=(/nxcase,nycase,1,1/) + startloc=(/1,1,inative,1/) work_a=hwork(1,:,:,ilevtot) @@ -4654,23 +4748,24 @@ subroutine gsi_fv3ncdf_write(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3file if(trim(name)=='xaxis_1') nx_phy=len if( nx_phy == nxcase )then allocate(work_b_tmp(nxcase,nycase)) - countloc_tmp=(/nxcase,nycase,1/) + countloc_tmp=(/nxcase,nycase,1,1/) phy_smaller_domain = .false. else allocate(work_b_tmp(nxcase-6,nycase-6)) - countloc_tmp=(/nxcase-6,nycase-6,1/) + countloc_tmp=(/nxcase-6,nycase-6,1,1/) phy_smaller_domain = .true. end if - startloc_tmp=(/1,1,ilev/) + startloc_tmp=(/1,1,ilev,1/) end if call check( nf90_inq_varid(gfile_loc,trim(varname),VarId) ) + call check( nf90_var_par_access(gfile_loc, VarId, nf90_independent)) if(index(vgsiname,"delzinc") > 0) then if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 - countloc=(/nxcase,ny_layout_len(nio),1/) + countloc=(/nxcase,ny_layout_len(nio),1,1/) allocate(work_b_layout(nxcase,ny_layout_len(nio))) call check( nf90_get_var(gfile_loc_layout(nio),VarId,work_b_layout,start = startloc, count = countloc) ) work_b(:,ny_layout_b(nio):ny_layout_e(nio))=work_b_layout @@ -4685,7 +4780,7 @@ subroutine gsi_fv3ncdf_write(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3file if(add_saved)then if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 - countloc=(/nxcase,ny_layout_len(nio),1/) + countloc=(/nxcase,ny_layout_len(nio),1,1/) allocate(work_b_layout(nxcase,ny_layout_len(nio))) call check( nf90_get_var(gfile_loc_layout(nio),VarId,work_b_layout,start = startloc, count = countloc) ) work_b(:,ny_layout_b(nio):ny_layout_e(nio))=work_b_layout @@ -4721,7 +4816,7 @@ subroutine gsi_fv3ncdf_write(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3file endif if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 - countloc=(/nxcase,ny_layout_len(nio),1/) + countloc=(/nxcase,ny_layout_len(nio),1,1/) allocate(work_b_layout(nxcase,ny_layout_len(nio))) work_b_layout=work_b(:,ny_layout_b(nio):ny_layout_e(nio)) call check( nf90_put_var(gfile_loc_layout(nio),VarId,work_b_layout, start = startloc, count = countloc) ) @@ -4795,12 +4890,13 @@ subroutine gsi_fv3ncdf_write_v1(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3f ! !$$$ end documentation block - use mpimod, only: mpi_rtype,mpi_comm_world,mype,mpi_info_null + use mpimod, only: npe, setcomm,mpi_integer,mpi_max,mpi_rtype,mpi_comm_world,mype,mpi_info_null use mod_fv3_lola, only: fv3_ll_to_h use mod_fv3_lola, only: fv3_h_to_ll use netcdf, only: nf90_open,nf90_close - use netcdf, only: nf90_write,nf90_inq_varid + use netcdf, only: nf90_write, nf90_netcdf4,nf90_mpiio,nf90_inq_varid use netcdf, only: nf90_put_var,nf90_get_var + use netcdf, only: nf90_independent,nf90_var_par_access use gsi_bundlemod, only: gsi_bundle use general_sub2grid_mod, only: sub2grid_info,general_sub2grid implicit none @@ -4824,6 +4920,10 @@ subroutine gsi_fv3ncdf_write_v1(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3f character(len=max_varname_length) :: varname,vgsiname integer(i_kind) nlatcase,nloncase,nxcase,nycase,countloc(3),startloc(3) + integer(i_kind):: iworld,iworld_group,nread,mpi_comm_read,i,ierror + integer(i_kind),dimension(npe):: members,members_read,mype_read_rank + logical:: procuse + mm1=mype+1 nloncase=grd_ionouv%nlon @@ -4838,7 +4938,30 @@ subroutine gsi_fv3ncdf_write_v1(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3f allocate( work_b(nlon_regional,nlat_regional)) allocate( workb2(nlon_regional,nlat_regional)) allocate( worka2(nlatcase,nloncase)) - call check ( nf90_open(filenamein,nf90_write,gfile_loc,comm=mpi_comm_world,info=MPI_INFO_NULL)) !clt + + procuse = .false. + members=-1 + members_read=-1 + if (kbgn<=kend) then + procuse = .true. + members(mm1) = mype + endif + + call mpi_allreduce(members,members_read,npe,mpi_integer,mpi_max,mpi_comm_world,ierror) + + nread=0 + mype_read_rank=-1 + do i=1,npe + if (members_read(i) >= 0) then + nread=nread+1 + mype_read_rank(nread) = members_read(i) + endif + enddo + + call setcomm(iworld,iworld_group,nread,mype_read_rank,mpi_comm_read,ierror) + + if (procuse) then + call check ( nf90_open(filenamein,ior(nf90_netcdf4,ior(nf90_write, nf90_mpiio)),gfile_loc,comm=mpi_comm_read,info=MPI_INFO_NULL)) !clt do ilevtot=kbgn,kend vgsiname=grd_ionouv%names(1,ilevtot) if(trim(vgsiname)=='amassi') cycle @@ -4848,7 +4971,6 @@ subroutine gsi_fv3ncdf_write_v1(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3f call getfv3lamfilevname(vgsiname,fv3filenamegin,filenamein2,varname) if(trim(filenamein) /= trim(filenamein2)) then write(6,*)'filenamein and filenamein2 are not the same as expected, stop' - call flush(6) call stop2(333) endif ilev=grd_ionouv%lnames(1,ilevtot) @@ -4862,6 +4984,7 @@ subroutine gsi_fv3ncdf_write_v1(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3f call check( nf90_inq_varid(gfile_loc,trim(varname),VarId) ) + call check( nf90_var_par_access(gfile_loc, VarId, nf90_independent)) call check( nf90_get_var(gfile_loc,VarId,work_b,start=startloc,count=countloc) ) if(index(vgsiname,"delzinc") > 0) then write(6,*)'delz is not in the cold start fiels with this option, incompatible setup , stop' @@ -4885,6 +5008,7 @@ subroutine gsi_fv3ncdf_write_v1(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3f call check( nf90_put_var(gfile_loc,VarId,work_b,start=startloc,count=countloc) ) enddo !ilevtot call check(nf90_close(gfile_loc)) + endif deallocate(work_b,work_a) deallocate(worka2,workb2) @@ -5420,7 +5544,7 @@ subroutine gsi_copy_bundle(bundi,bundo) character(len=max_varname_length),dimension(:),allocatable:: target_name_vars3d character(len=max_varname_length) ::varname real(r_kind),dimension(:,:,:),pointer:: pvar3d=>NULL() - real(r_kind),dimension(:,:,:),pointer:: pvar2d =>NULL() + real(r_kind),dimension(:,:),pointer:: pvar2d =>NULL() integer(i_kind):: src_nc3d,src_nc2d,target_nc3d,target_nc2d integer(i_kind):: ivar,jvar,istatus src_nc3d=bundi%n3d diff --git a/ush/sub_hera b/ush/sub_hera index 610756af00..c94b734596 100755 --- a/ush/sub_hera +++ b/ush/sub_hera @@ -120,10 +120,10 @@ echo "#SBATCH --output=$output" echo "#SBATCH --job-name=$jobname" >> $cfile echo "#SBATCH --qos=$queue" >> $cfile echo "#SBATCH --time=$timew" >> $cfile -echo "#SBATCH --nodes=$nodes --ntasks-per-node=$procs --cpus-per-task=$threads" >> $cfile +echo "#SBATCH --nodes=$nodes --ntasks-per-node=$procs --exclusive" >> $cfile #echo "#SBATCH -j oe" >> $cfile echo "#SBATCH --account=$accnt" >> $cfile -echo "#SBATCH --mem=0" >> $cfile +#cltorg echo "#SBATCH --mem=0" >> $cfile #echo "#SBATCH -V" >> $cfile #echo "#PBS -d" >> $cfile #. $exec >> $cfile @@ -143,7 +143,6 @@ echo "module list" >> $cfile echo "" >>$cfile cat $exec >> $cfile - if [[ $nosub = YES ]];then cat $cfile exit diff --git a/ush/sub_hercules b/ush/sub_hercules index 459b480559..78a0f5daee 100755 --- a/ush/sub_hercules +++ b/ush/sub_hercules @@ -111,7 +111,7 @@ echo "#SBATCH --job-name=$jobname" echo "#SBATCH --qos=$queue" >> $cfile echo "#SBATCH --partition=$partition" >> $cfile echo "#SBATCH --time=$timew" >> $cfile -echo "#SBATCH --nodes=$nodes --ntasks-per-node=$procs --cpus-per-task=$threads" >> $cfile +echo "#SBATCH --nodes=$nodes --ntasks-per-node=$procs --exclusive" >> $cfile echo "#SBATCH --account=$accnt" >> $cfile echo "" >>$cfile @@ -131,7 +131,7 @@ echo "module use $modulefiles" >> $cfile echo "module load gsi_hercules.intel" >> $cfile #TODO reenable I_MPI_EXTRA_FILESYSTEM once regional ctests can properly handle parallel I/O on Hercules echo "unset I_MPI_EXTRA_FILESYSTEM" >> $cfile -echo "" >> $cfile + cat $exec >> $cfile if [[ $nosub = YES ]];then diff --git a/ush/sub_jet b/ush/sub_jet index 9bd60486f6..96f3eae9b2 100755 --- a/ush/sub_jet +++ b/ush/sub_jet @@ -108,7 +108,7 @@ echo "#SBATCH --output=$output" echo "#SBATCH --job-name=$jobname" >> $cfile echo "#SBATCH --qos=$queue" >> $cfile echo "#SBATCH --time=$timew" >> $cfile -echo "#SBATCH --nodes=$nodes --ntasks-per-node=$procs --cpus-per-task=$threads" >> $cfile +echo "#SBATCH --nodes=$nodes --ntasks-per-node=$procs --exclusive" >> $cfile echo "#SBATCH --account=$accnt" >> $cfile echo "#SBATCH --mem=0" >> $cfile echo "#SBATCH --partition=kjet" >> $cfile diff --git a/ush/sub_orion b/ush/sub_orion index 5a13f54845..b810576379 100755 --- a/ush/sub_orion +++ b/ush/sub_orion @@ -111,7 +111,7 @@ echo "#SBATCH --job-name=$jobname" echo "#SBATCH --qos=$queue" >> $cfile echo "#SBATCH --partition=$partition" >> $cfile echo "#SBATCH --time=$timew" >> $cfile -echo "#SBATCH --nodes=$nodes --ntasks-per-node=$procs --cpus-per-task=$threads" >> $cfile +echo "#SBATCH --nodes=$nodes --ntasks-per-node=$procs --exclusive" >> $cfile echo "#SBATCH --account=$accnt" >> $cfile echo "" >>$cfile