From 94706a711b084cbd2589c14db5002ac3e11131ef Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Sat, 9 Jul 2022 16:23:17 -0500 Subject: [PATCH 01/12] Added endif statement to controlvec that had been inadvertently removed. --- src/enkf/controlvec.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/enkf/controlvec.f90 b/src/enkf/controlvec.f90 index 6a825fe278..4f8c2c7255 100644 --- a/src/enkf/controlvec.f90 +++ b/src/enkf/controlvec.f90 @@ -374,6 +374,7 @@ subroutine write_control(no_inflate_flag) enddo endif end if + end if if (.not. paranc) then if (write_fv3_incr) then if (global_2mDA) then From 45d3aebc2022f08843e2cac693e63997578c8aa6 Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Sat, 9 Jul 2022 16:27:33 -0500 Subject: [PATCH 02/12] Added end if statement back to controlvec that had inavertently been removed. --- src/enkf/controlvec.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/enkf/controlvec.f90 b/src/enkf/controlvec.f90 index 6a825fe278..4f8c2c7255 100644 --- a/src/enkf/controlvec.f90 +++ b/src/enkf/controlvec.f90 @@ -374,6 +374,7 @@ subroutine write_control(no_inflate_flag) enddo endif end if + end if if (.not. paranc) then if (write_fv3_incr) then if (global_2mDA) then From 51d9b4be7fce24dbc6b1a77bc1e9fd09b1467452 Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Wed, 13 Jul 2022 12:50:56 -0500 Subject: [PATCH 03/12] Separated reading in 2m obs from doing sfc update. Changes observer to be able to read 2m obs together with other obs. --- src/gsi/buddycheck_mod.f90 | 8 +++--- src/gsi/general_read_gfsatm.f90 | 28 +++++++++++---------- src/gsi/gsimod.F90 | 9 ++----- src/gsi/jfunc.f90 | 6 ++--- src/gsi/netcdfgfs_io.f90 | 12 ++++----- src/gsi/read_prepbufr.f90 | 10 ++++---- src/gsi/setupt.f90 | 44 ++++++++++++++++----------------- 7 files changed, 57 insertions(+), 60 deletions(-) diff --git a/src/gsi/buddycheck_mod.f90 b/src/gsi/buddycheck_mod.f90 index 4ff46b6f58..d23198a3ed 100644 --- a/src/gsi/buddycheck_mod.f90 +++ b/src/gsi/buddycheck_mod.f90 @@ -60,7 +60,7 @@ subroutine buddy_check_t(is,data,luse,mype,nele,nobs,muse,buddyuse) use aircraftinfo, only: aircraft_t_bc_pof,aircraft_t_bc, & aircraft_t_bc_ext use convinfo, only: ictype - use jfunc, only: do_global_2mDA + use jfunc, only: use_2m_obs implicit none ! !INPUT PARAMETERS: @@ -281,8 +281,8 @@ subroutine buddy_check_t(is,data,luse,mype,nele,nobs,muse,buddyuse) f10ges,u10ges,v10ges, t2ges, q2ges, regime, iqtflg) tges = t2ges - elseif(sfctype .and. do_global_2mDA) then -! SCENARIO 2: obs is sfctype, and do_global_2mDA scheme is on. + elseif(sfctype .and. use_2m_obs) then +! SCENARIO 2: obs is sfctype, and use_2m_obs scheme is on. ! 2m forecast has been read from the sfc guess files ! note: this block should match that in setupt. @@ -467,7 +467,7 @@ subroutine init_vars_ write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus call stop2(999) endif - if( do_global_2mDA) then + if( use_2m_obs) then ! get t2m ... varname='t2m' call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus) diff --git a/src/gsi/general_read_gfsatm.f90 b/src/gsi/general_read_gfsatm.f90 index ed3cee326f..576e502f76 100755 --- a/src/gsi/general_read_gfsatm.f90 +++ b/src/gsi/general_read_gfsatm.f90 @@ -1687,7 +1687,7 @@ subroutine general_read_gfsatm_nems(grd,sp_a,filename,uvflag,vordivflag,zflag, & end subroutine general_read_gfsatm_nems -subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag,do_global_2mDA, & +subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag,use_2m_obs, & gfs_bundle,init_head,iret_read) !$$$ subprogram documentation block ! . . . . @@ -1748,7 +1748,7 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag,do_g type(spec_vars) ,intent(in ) :: sp_a character(*) ,intent(in ) :: filename logical ,intent(in ) :: uvflag,zflag,vordivflag,init_head - logical ,intent(in ) :: do_global_2mDA + logical ,intent(in ) :: use_2m_obs integer(i_kind) ,intent( out) :: iret_read type(gsi_bundle) ,intent(inout) :: gfs_bundle @@ -1810,21 +1810,21 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag,do_g read_z = .false. if ( mype == 0 ) write(6,* ) & trim(my_name), ': reading 2m variables from ', trim(filename) - if (.not. do_global_2mDA) then - write(6,*) 'error, requesting sfc file read, but not do_global_2mDA. Check options' + if (.not. use_2m_obs) then + write(6,*) 'error, requesting sfc file read, but not use_2m_obs. Check options' call stop2(101) endif else read_2m = .false. - if (do_global_2mDA) then - read_z =.true. ! over-ride, as this is the only var we want. - if ( mype == 0 ) write(6,* ) & - trim(my_name), ': reading z only from ', trim(filename) - else + !if (use_2m_obs) then + ! read_z =.true. ! over-ride, as this is the only var we want. + ! if ( mype == 0 ) write(6,* ) & + ! trim(my_name), ': reading z only from ', trim(filename) + !else read_z = zflag if ( mype == 0 ) write(6,* ) & trim(my_name), ': reading atmos variables from ', trim(filename) - endif + !endif endif @@ -2053,13 +2053,15 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag,do_g endif endif ! CSD for 2 m case, icm = 2. this doesn't get triggered. - if ( do_global_2mDA .or. ( (.not. do_global_2mDA ) .and. ( icount == icm)) ) then + !if ( use_2m_obs .or. ( (.not. use_2m_obs ) .and. ( icount == icm)) ) then + if ( read_2m .or. ( (.not. read_2m ) .and. ( icount == icm)) ) then call general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, & icount,iflag,ilev,work,uvflag,vordivflag) endif endif - if ((.not. read_2m) .and. (.not. do_global_2mDA)) then + !if ((.not. read_2m) .and. (.not. use_2m_obs)) then + if (.not. read_2m) then icount=icount+1 iflag(icount)=2 @@ -2500,7 +2502,7 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag,do_g !enddo ! Load u->div and v->vor slot when uv are used instead - if ( (.not. read_2m) .and. (.not. do_global_2mDA)) then + if ( .not. read_2m ) then if ( uvflag ) then call gsi_bundlegetpointer(gfs_bundle,'u' ,ptr3d,ier) if ( ier == 0 ) then diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index 8bfde3dad6..bf75561fb0 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -98,7 +98,7 @@ module gsimod factv,factl,factp,factg,factw10m,facthowv,factcldch,niter,niter_no_qc,biascor,& init_jfunc,qoption,cwoption,switch_on_derivatives,tendsflag,jiterstart,jiterend,R_option,& bcoption,diurnalbc,print_diag_pcg,tsensible,diag_precon,step_start,pseudo_q2,& - clip_supersaturation,cnvw_option,do_global_2mDA + clip_supersaturation,cnvw_option,use_2m_obs use state_vectors, only: init_anasv,final_anasv use control_vectors, only: init_anacv,final_anacv,nrf,nvars,nrf_3d,cvars3d,cvars2d,& nrf_var,lcalc_gfdl_cfrac,incvars_to_zero,incvars_zero_strat,incvars_efold @@ -1020,7 +1020,7 @@ module gsimod ! l_foreaft_thin - separate TDR fore/aft scan for thinning namelist/obs_input/dmesh,time_window_max,time_window_rad, & - ext_sonde,l_foreaft_thin + ext_sonde,l_foreaft_thin,use_2m_obs ! SINGLEOB_TEST (one observation test case setup): ! maginnov - magnitude of innovation for one ob @@ -1533,7 +1533,6 @@ module gsimod ! fac_dtl - index to apply diurnal thermocline layer or not: 0 = no; 1 = yes. ! fac_tsl - index to apply thermal skin layer or not: 0 = no; 1 = yes. namelist/nst/nst_gsi,nstinfo,zsea1,zsea2,fac_dtl,fac_tsl - namelist/global_2mDA/do_global_2mDA !EOC @@ -1650,7 +1649,6 @@ subroutine gsimain_initialize read(5,rapidrefresh_cldsurf) read(5,chem) read(5,nst) - read(5,global_2mDA) #else ! Initialize table of instruments and data types call obsmod_init_instr_table(nhr_assimilation,ndat,rcname='gsiparm.anl') @@ -1702,9 +1700,6 @@ subroutine gsimain_initialize read(11,nst,iostat=ios) if(ios/=0) call die(myname_,'read(nst)',ios) - read(11,global_2mDA,iostat=ios) - if(ios/=0) call die(myname_,'read(global_2mDA)',ios) - close(11) #endif if(jcap > jcap_cut)then diff --git a/src/gsi/jfunc.f90 b/src/gsi/jfunc.f90 index 58bd19c1ec..bafa59b85a 100644 --- a/src/gsi/jfunc.f90 +++ b/src/gsi/jfunc.f90 @@ -136,12 +136,12 @@ module jfunc public :: pseudo_q2 public :: varq public :: cnvw_option - public :: do_global_2mDA + public :: use_2m_obs logical first,last,switch_on_derivatives,tendsflag,print_diag_pcg,tsensible,diag_precon logical clip_supersaturation,R_option logical pseudo_q2,limitqobs - logical do_global_2mDA + logical use_2m_obs logical cnvw_option integer(i_kind) iout_iter,miter,iguess,nclen,qoption,cwoption integer(i_kind) jiter,jiterstart,jiterend,iter @@ -251,7 +251,7 @@ subroutine init_jfunc ! option for including convective clouds in the all-sky assimilation cnvw_option=.false. - do_global_2mDA=.false. + use_2m_obs=.false. return end subroutine init_jfunc diff --git a/src/gsi/netcdfgfs_io.f90 b/src/gsi/netcdfgfs_io.f90 index 7198989ead..12a8b8f779 100644 --- a/src/gsi/netcdfgfs_io.f90 +++ b/src/gsi/netcdfgfs_io.f90 @@ -128,7 +128,7 @@ subroutine read_ use general_sub2grid_mod, only: sub2grid_info,general_sub2grid_create_info,general_sub2grid_destroy_info use mpimod, only: npe,mype use cloud_efr_mod, only: cloud_calc_gfs,set_cloud_lower_bound - use jfunc, only: do_global_2mDA + use jfunc, only: use_2m_obs implicit none character(len=*),parameter::myname_=myname//'*read_' @@ -178,7 +178,7 @@ subroutine read_ ! Allocate bundle used for reading members call gsi_gridcreate(atm_grid,lat2,lon2,nsig) - if (do_global_2mDA) then + if (use_2m_obs) then call gsi_bundlecreate(atm_bundle,atm_grid,'aux-atm-read',istatus,names2d=vars2d_with2m,names3d=vars3d) else call gsi_bundlecreate(atm_bundle,atm_grid,'aux-atm-read',istatus,names2d=vars2d,names3d=vars3d) @@ -191,20 +191,20 @@ subroutine read_ do it=1,nfldsig ! Read background fields into bundle - if (do_global_2mDA) then ! current read_sfc routines called from differebt part of + if (use_2m_obs) then ! current read_sfc routines called from differebt part of ! of the code, can't easily read into the met-bundle ! wrote a new routine here write(filename,'(''sfcf'',i2.2)') ifilesig(it) call general_read_gfsatm_nc(grd_t,sp_a,filename,.true.,.true.,.true.,& - do_global_2mDA,atm_bundle,.true.,istatus) + use_2m_obs,atm_bundle,.true.,istatus) write(filename,'(''sigf'',i2.2)') ifilesig(it) call general_read_gfsatm_nc(grd_t,sp_a,filename,.true.,.true.,.true.,& - do_global_2mDA,atm_bundle,.true.,istatus) + use_2m_obs,atm_bundle,.true.,istatus) else write(filename,'(''sigf'',i2.2)') ifilesig(it) call general_read_gfsatm_nc(grd_t,sp_a,filename,.true.,.true.,.true.,& - do_global_2mDA,atm_bundle,.true.,istatus) + use_2m_obs,atm_bundle,.true.,istatus) endif inithead=.false. diff --git a/src/gsi/read_prepbufr.f90 b/src/gsi/read_prepbufr.f90 index 0a2b027440..190cbf7ea0 100644 --- a/src/gsi/read_prepbufr.f90 +++ b/src/gsi/read_prepbufr.f90 @@ -211,7 +211,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& use hilbertcurve,only: init_hilbertcurve, accum_hilbertcurve, & apply_hilbertcurve,destroy_hilbertcurve use ndfdgrids,only: init_ndfdgrid,destroy_ndfdgrid,relocsfcob,adjust_error - use jfunc, only: tsensible, do_global_2mDA + use jfunc, only: tsensible, use_2m_obs use deter_sfc_mod, only: deter_sfc_type,deter_sfc2 use gsi_nstcouplermod, only: nst_gsi,nstinfo use gsi_nstcouplermod, only: gsi_nstcoupler_deter @@ -1617,7 +1617,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& ! versus sensible temperature if(tob) then ! use tvirtual if tsensible flag not set, and not in either 2Dregional or global_2m DA mode - if ( (.not. tsensible) .and. .not. (twodvar_regional .or. do_global_2mDA) ) then + if ( (.not. tsensible) .and. .not. (twodvar_regional .or. use_2m_obs) ) then do k=1,levs tvflg(k)=one ! initialize as sensible @@ -1917,7 +1917,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& ! CSD - temporary hack ( move to prepbufr pre-processing) ! Over-ride QM=9 for sfc obs - if (sfctype .and. do_global_2mDA ) then + if (sfctype .and. use_2m_obs ) then if (tob .and. qm == 9 ) qm = 2 ! 2=not checked if (qob .and. qm == 9 ) qm = 2 ! 2=not checked endif @@ -2079,7 +2079,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& if (aircraftobs .and. aircraft_t_bc .and. acft_profl_file) then call errormod_aircraft(pqm,tqm,levs,plevs,errout,k,presl,dpres,nsig,lim_qm,hdr3) else - if (.not. (sfctype .and. do_global_2mDA)) then + if (.not. (sfctype .and. use_2m_obs)) then call errormod(pqm,tqm,levs,plevs,errout,k,presl,dpres,nsig,lim_qm) endif end if @@ -2336,7 +2336,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& end if qobcon=obsdat(2,k)*convert tdry=r999 - if (sfctype .and. do_global_2mDA) then + if (sfctype .and. use_2m_obs) then if (tqm(k) < 10) tdry=(obsdat(3,k)+t0c)/(one+fv*qobcon) else if (tqm(k)179.and.itype<190).or.(itype>=192.and.itype<=199) do l=k+1,nobs - if (twodvar_regional .or. (do_global_2mDA .and. sfctype) ) then + if (twodvar_regional .or. (use_2m_obs .and. sfctype) ) then duplogic=data(ilat,k) == data(ilat,l) .and. & data(ilon,k) == data(ilon,l) .and. & data(ier,k) < r1000 .and. data(ier,l) < r1000 .and. & @@ -473,19 +473,19 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav end do end do -! Run a buddy-check (does not work for do_global_2mDA, +! Run a buddy-check (does not work for use_2m_obs, ! pre-existing routine crashes with an MPI problem ) ! Note: current params have buddy radius of 108 km, max diff of 8 K. ! gross error check removes O-F > 7., so this is probably removing ! most obs that fail the buddy check already ! GSD uses the buddy check to expand gross error check for obs that pass ! the buddy check. Not sure if we want to do this globally. Turn off for now. - !if ( (twodvar_regional .or. do_global_2mDA) .and. buddycheck_t) & + !if ( (twodvar_regional .or. use_2m_obs) .and. buddycheck_t) & if ( (twodvar_regional) .and. buddycheck_t) & call buddy_check_t(is,data,luse,mype,nele,nobs,muse,buddyuse) -! for 2mDA, remove obs that failed the buddcheck - ! if (do_global_2mDA) then +! for 2m obs, remove obs that failed the buddcheck + ! if (use_2m_obs) then ! do i = 1,nobs ! ikx=nint(data(ikxx,k)) ! itype=ictype(ikx) @@ -674,17 +674,17 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav end if ! Interpolate log(ps) & log(pres) at mid-layers to obs locations/times - ! for global_2mDA case, skipping reading the pressure levels into the met guess (to speed up the observer) + ! for 2m obs case, skipping reading the pressure levels into the met guess (to speed up the observer) ! these are not currently used here, so skip for now. - if (.not. do_global_2mDA) then + if (sfctype .and. .not. use_2m_obs ) then call tintrp2a11(ges_ps,psges,dlat,dlon,dtime,hrdifsig,& mype,nfldsig) call tintrp2a1(ges_lnprsl,prsltmp,dlat,dlon,dtime,hrdifsig,& nsig,mype,nfldsig) drpx=zero - !if(sfctype .and. .not.twodvar_regional .and. .not. do_global_2mDA) then + !if(sfctype .and. .not.twodvar_regional .and. .not. use_2m_obs ) then if(sfctype .and. .not.twodvar_regional ) then drpx=abs(one-((one/exp(dpres-log(psges))))**rd_over_cp)*t0c end if @@ -729,8 +729,8 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav prsltmp2(2), tvtmp(2), qtmp(2), hsges(1), roges, msges, & f10ges,u10ges,v10ges, t2ges, q2ges, regime, iqtflg) tges = t2ges - elseif (sfctype .and. do_global_2mDA) then -! SCENARIO 2: obs is sfctype, and do_global_2mDA scheme is on. + elseif (sfctype .and. use_2m_obs ) then +! SCENARIO 2: obs is sfctype, and use_2m_obs scheme is on. ! 2m forecast has been read from the sfc guess files ! note: any changes to ges_t2 and tob in this block should also be made in buddy_check_t @@ -754,7 +754,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav else -! SCENARIO 3: obs is sfctype, and neither sfcmodel nor do_global_2mDA is chosen +! SCENARIO 3: obs is sfctype, and neither sfcmodel nor use_2m_obs is chosen ! .or. obs is not sfctype ! ! SCENARIO 3a: obs is a virtual temp. if(iqtflg)then @@ -840,13 +840,13 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav call grdcrd1(sfcchk,prsltmp(1),nsig,-1) ! Check to see if observations is above the top of the model (regional mode) - if(sfctype .and. .not. do_global_2mDA )then + if(sfctype .and. .not. use_2m_obs )then if(abs(dpres)>four) drpx=1.0e10_r_kind pres_diff=prest-r10*psges if (twodvar_regional .and. abs(pres_diff)>=r1000) drpx=1.0e10_r_kind end if - if (.not. (do_global_2mDA .and. sfctype) ) then + if (.not. (use_2m_obs .and. sfctype) ) then rlow=max(sfcchk-dpres,zero) ! sfcchk = k of sfc, dpres k of obs ! linear variation of observation ramp [between grid points 1(~3mb) and 15(~45mb) below the surface] if(l_sfcobserror_ramp_t) then @@ -870,7 +870,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ! inflate error for uncertainty in the terrain adjustment lapse_error = 0. - if ( do_global_2mDA .and. sfctype) then + if ( use_2m_obs .and. sfctype) then if (abs(delta_z)300 m. do not assim. ! inflate obs error to account for error in lapse_rate ! also include some representativity error here (assuming @@ -881,7 +881,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav endif endif -! note: for global_2mDA option, three middle terms in denominator will be zero +! note: for use_2m_obs tion, three middle terms in denominator will be zero ! (elevation mistmatch between obs and model dealt with elsewhere) ratio_errors=error/(data(ier,i)+drpx+1.0e6_r_kind*rhgh+r8*ramp + lapse_error) @@ -889,7 +889,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ! Compute innovation !if(i_use_2mt4b>0 .and. sfctype) then - if(sfctype .and. (( i_use_2mt4b>0) .or. do_global_2mDA) ) then + if(sfctype .and. (( i_use_2mt4b>0) .or. use_2m_obs ) ) then ddiff = tob-tges2m tges = tges2m ! not necessary? else @@ -1518,7 +1518,7 @@ subroutine init_vars_ call stop2(999) endif endif - if( do_global_2mDA) then + if( use_2m_obs ) then ! get t2m ... varname='t2m' call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus) @@ -1538,7 +1538,7 @@ subroutine init_vars_ call stop2(999) endif endif - if(i_use_2mt4b>0 .or. do_global_2mDA) then + if(i_use_2mt4b>0 .or. use_2m_obs ) then ! get q2m ... varname='q2m' call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus) @@ -1753,7 +1753,7 @@ subroutine contents_netcdf_diag_(odiag) call nc_diag_metadata("Height", sngl(data(iobshgt,i)) ) ! this is the obs height (= stn elevation, before being corrected) call nc_diag_metadata("Time", sngl(dtime-time_offset)) call nc_diag_metadata("Prep_QC_Mark", sngl(data(iqc,i)) ) - if (do_global_2mDA == .true. ) then + if (use_2m_obs == .true. ) then idg = 1 else idg = 0 @@ -1777,7 +1777,7 @@ subroutine contents_netcdf_diag_(odiag) call nc_diag_metadata("Errinv_Input", sngl(errinv_input) ) call nc_diag_metadata("Errinv_Adjust", sngl(errinv_adjst) ) call nc_diag_metadata("Errinv_Final", sngl(errinv_final) ) - if (do_global_2mDA) then + if (use_2m_obs ) then call nc_diag_metadata("Observation", sngl(tob) ) call nc_diag_metadata("Observation_Before_Elev_Correction", sngl(data(itob,i)) ) else @@ -1821,7 +1821,7 @@ subroutine contents_netcdf_diag_(odiag) call nc_diag_data2d("ObsDiagSave_obssen", odiag%obssen ) endif - if (twodvar_regional .or. do_global_2mDA ) then + if (twodvar_regional .or. use_2m_obs ) then call nc_diag_metadata("Dominant_Sfc_Type", data(idomsfc,i) ) ! this is the model height interpolated to the obs location in read_prepbufr call nc_diag_metadata("Model_Terrain", data(izz,i) ) From d5fae407745b1fdeb5dde02720053db011b36c73 Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Wed, 13 Jul 2022 13:18:26 -0500 Subject: [PATCH 04/12] Changed setupt output QC marks back to original values. --- src/gsi/setupt.f90 | 23 +++++------------------ 1 file changed, 5 insertions(+), 18 deletions(-) diff --git a/src/gsi/setupt.f90 b/src/gsi/setupt.f90 index a4e07eeadd..4342005c77 100644 --- a/src/gsi/setupt.f90 +++ b/src/gsi/setupt.f90 @@ -30,7 +30,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav use m_obsdiagNode, only: obsdiagNode_assert use obsmod, only: sfcmodel,perturb_obs,oberror_tune,lobsdiag_forenkf,ianldate,& - lobsdiagsave,nobskeep,lobsdiag_allocated,time_offset,aircraft_recon + lobsdiagsave,nobskeep,lobsdiag_allocated,time_offset,aircraft_recon use m_obsNode, only: obsNode use m_tNode, only: tNode use m_tNode, only: tNode_appendto @@ -110,11 +110,11 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ! !INPUT/OUTPUT PARAMETERS: + ! array containing information ... real(r_kind),dimension(npres_print,nconvtype,5,3), intent(inout) :: bwork ! about o-g stats real(r_kind),dimension(100+7*nsig) , intent(inout) :: awork ! for data counts and gross checks - ! !DESCRIPTION: For temperature observations, this routine ! \begin{enumerate} ! \item reads obs assigned to given mpi task (geographic region), @@ -365,7 +365,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ! call GSD terrain match for surface temperature observation if(l_gsd_terrain_match_surftobs) then - call gsd_terrain_match_surfTobs(mype,nele,nobs,data) + call gsd_terrain_match_surfTobs(mype,nele,nobs,data) endif ! index information for data array (see reading routine) @@ -984,7 +984,6 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav end if endif - if (sfctype .and. i_sfct_gross==1) then ! extend the threshold for surface T if(i_use_2mt4b<=0) tges2m=tges @@ -1017,7 +1016,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ratio_errors = ratio_errors/sqrt(dup(i)) end if endif - + if (ratio_errors*error <=tiny_r_kind) muse(i)=.false. if (nobskeep>0 .and. luse_obsdiag) call obsdiagNode_get(my_diag, jiter=nobskeep, muse=muse(i)) @@ -1753,19 +1752,7 @@ subroutine contents_netcdf_diag_(odiag) call nc_diag_metadata("Height", sngl(data(iobshgt,i)) ) ! this is the obs height (= stn elevation, before being corrected) call nc_diag_metadata("Time", sngl(dtime-time_offset)) call nc_diag_metadata("Prep_QC_Mark", sngl(data(iqc,i)) ) - if (use_2m_obs == .true. ) then - idg = 1 - else - idg = 0 - endif - if (sfctype == .true. ) then - ist = 1 - else - ist = 0 - endif - write(flag_char,"(3I1)") ist, int(data(iqt,i)), idg - call nc_diag_metadata("Setup_Flags", flag_char ) - !call nc_diag_metadata("Setup_Flag", sngl(data(iqt,i)) ) + call nc_diag_metadata("Setup_QC_Mark", sngl(data(iqt,i)) ) ! this is the virtual temp flag call nc_diag_metadata("Prep_Use_Flag", sngl(data(iuse,i)) ) if(muse(i)) then call nc_diag_metadata("Analysis_Use_Flag", sngl(one) ) From 2c6d9e59a44da768c9f9946166f45259c0b2643c Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Thu, 14 Jul 2022 10:19:39 -0500 Subject: [PATCH 05/12] Reverted drpx and dpres code for use_2m_obs option, and changes to compile. --- src/gsi/setupt.f90 | 37 +++++++++++++++++-------------------- src/gsi/update_guess.f90 | 6 +++--- 2 files changed, 20 insertions(+), 23 deletions(-) diff --git a/src/gsi/setupt.f90 b/src/gsi/setupt.f90 index c8180fef8a..e8325d5b84 100644 --- a/src/gsi/setupt.f90 +++ b/src/gsi/setupt.f90 @@ -58,7 +58,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav use guess_grids, only: nfldsig, hrdifsig,ges_lnprsl,& geop_hgtl,ges_tsen,pbl_height - use state_vectors, only: svars3d, levels + use state_vectors, only: svars3d, levels, ns3d, svars2d use constants, only: zero, one, four,t0c,rd_over_cp,three,rd_over_cp_mass,ten use constants, only: tiny_r_kind,half,two @@ -335,9 +335,6 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav type(obsLList),pointer,dimension(:):: thead - character(len=3) :: flag_char - integer :: idg, ist, itm - real(r_kind) :: delta_z, lapse_error real(r_kind), parameter :: T_lapse = -0.0045 ! standard lapse rate, K/m ! should match value in buddy_check_t @@ -677,18 +674,18 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav call tintrp2a1(ges_lnprsl,prsltmp,dlat,dlon,dtime,hrdifsig,& nsig,mype,nfldsig) - drpx=zero - if(sfctype ) then - if ( use_2m_obs) then - dpres = one ! put obs on surface (not sure about this) - elseif ( .not.twodvar_regional) then - drpx=abs(one-((one/exp(dpres-log(psges))))**rd_over_cp)*t0c - endif - end if - + if ( .not. use_2m_obs) then + drpx=zero + if(sfctype .and. .not.twodvar_regional) then + drpx=abs(one-((one/exp(dpres-log(psges))))**rd_over_cp)*t0c + end if -! Put obs pressure in correct units to get grid coord. number - call grdcrd1(dpres,prsltmp(1),nsig,-1) + ! Put obs pressure in correct units to get grid coord. number + call grdcrd1(dpres,prsltmp(1),nsig,-1) + else + drpx = zero + dpres = one ! put obs at surface + endif ! Implementation of forward model ---------- @@ -747,13 +744,13 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav data(istnelv,i) = data(izz,i) if(save_jacobian) then - t2m_ind = getindex(svars2d, 't2m') - if (t2m_ind < 0) then + t_ind = getindex(svars2d, 't2m') + if (t_ind < 0) then print *, 'Error: no variable t2m in state vector.Exiting.' call stop2(1300) endif - dhx_dx%st_ind(1) = sum(levels(1:ns3d)) + t2m_ind - dhx_dx%end_ind(1) = sum(levels(1:ns3d)) + t2m_ind + dhx_dx%st_ind(1) = sum(levels(1:ns3d)) + t_ind + dhx_dx%end_ind(1) = sum(levels(1:ns3d)) + t_ind dhx_dx%val(1) = one endif else @@ -883,7 +880,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav endif endif -! note: for use_2m_obs tion, three middle terms in denominator will be zero +! note: for use_2m_obs, three middle terms in denominator will be zero ! (elevation mistmatch between obs and model dealt with elsewhere) ratio_errors=error/(data(ier,i)+drpx+1.0e6_r_kind*rhgh+r8*ramp + lapse_error) diff --git a/src/gsi/update_guess.f90 b/src/gsi/update_guess.f90 index 189d1afece..c8656053b5 100644 --- a/src/gsi/update_guess.f90 +++ b/src/gsi/update_guess.f90 @@ -113,7 +113,7 @@ subroutine update_guess(sval,sbias) use mpimod, only: mype use constants, only: zero,one,fv,max_varname_length,qmin,qcmin,tgmin,& r100,one_tenth,tiny_r_kind - use jfunc, only: iout_iter,bcoption,tsensible,clip_supersaturation,superfact,do_global_2mDA + use jfunc, only: iout_iter,bcoption,tsensible,clip_supersaturation,superfact,use_2m_obs use gridmod, only: lat2,lon2,nsig,& regional,twodvar_regional,regional_ozone,& l_reg_update_hydro_delz @@ -454,7 +454,7 @@ subroutine update_guess(sval,sbias) endif call gsd_update_soil_tq(tinc_1st,is_t,qinc_1st,is_q,it) endif ! l_gsd_soilTQ_nudge - if ( (i_use_2mt4b > 0.or.do_global_2mDA) .and. is_t>0) then + if ( (i_use_2mt4b > 0.or. use_2m_obs) .and. is_t>0) then do j=1,lon2 do i=1,lat2 tinc_1st(i,j)=p_tv(i,j,1) @@ -462,7 +462,7 @@ subroutine update_guess(sval,sbias) end do call gsd_update_t2m(tinc_1st,it) endif ! l_gsd_t2m_adjust - if ( (i_use_2mq4b > 0.or.do_global_2mDA) .and. is_q>0) then + if ( (i_use_2mq4b > 0.or. use_2m_obs) .and. is_q>0) then do j=1,lon2 do i=1,lat2 qinc_1st(i,j)=p_q(i,j,1) From 16f9d2c6d1730df4536bad001d6511f5c5b7151a Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Tue, 26 Jul 2022 09:45:01 -0500 Subject: [PATCH 06/12] Midway through generalising the 2m DA case: reverted nlev to always be number of atmos levels merged read routines to use single routine to read from atm and/or sfc file works for reading control variables. Can read state variable, but perturbations are not being counted properly. --- src/enkf/controlvec.f90 | 39 +- src/enkf/gridinfo_gfs.f90 | 29 +- src/enkf/gridio_gfs.f90 | 1158 ++++++++++++++++++------------------- src/enkf/mpi_readobs.f90 | 2 +- src/enkf/params.f90 | 15 +- src/enkf/readconvobs.f90 | 18 +- src/enkf/statevec.f90 | 13 +- src/gsi/setupt.f90 | 4 +- 8 files changed, 606 insertions(+), 672 deletions(-) diff --git a/src/enkf/controlvec.f90 b/src/enkf/controlvec.f90 index 4f8c2c7255..4789bdff94 100644 --- a/src/enkf/controlvec.f90 +++ b/src/enkf/controlvec.f90 @@ -47,7 +47,7 @@ module controlvec mpi_integer,mpi_wtime,mpi_status,mpi_real8 use gridio, only: readgriddata, readgriddata_pnc, writegriddata, writegriddata_pnc, & - writeincrement, writeincrement_pnc, readgriddata_2mDA, & + writeincrement, writeincrement_pnc, & writegriddata_2mDA, writeincrement_2mDA use gridinfo, only: getgridinfo, gridinfo_cleanup, & npts, vars3d_supported, vars2d_supported @@ -132,7 +132,7 @@ subroutine init_controlvec() cvars3d(nc3d) = trim(adjustl(var)) clevels(nc3d) = ilev + clevels(nc3d-1) else - if (nproc .eq. 0) print *,'Error: only ', nlevs, ' and ', nlevs+1,' number of levels is supported in current version, got ',ilev + if (nproc .eq. 0) print *,'Error: controlvec only ', nlevs, ' and ', nlevs+1,' number of levels is supported in current version, got ',ilev call stop2(503) endif enddo @@ -213,9 +213,10 @@ subroutine read_control() ! read in whole control vector on i/o procs - keep in memory ! (needed in write_ensemble) allocate(grdin(npts,ncdim,nbackgrounds,nanals_per_iotask)) -if (.not. global_2mDA) allocate(qsat(npts,nlevs,nbackgrounds,nanals_per_iotask)) +q_ind = getindex(cvars3d, 'q') +if (q_ind > 0) allocate(qsat(npts,nlevs,nbackgrounds,nanals_per_iotask)) if (paranc) then - if (global_2mDA) then + if (global_2mDA) then ! flag for parallel read support print *,'paranc not supported for global_2mDA' call mpi_barrier(mpi_comm_world,ierr) call mpi_finalize(ierr) @@ -227,16 +228,12 @@ subroutine read_control() if (nproc <= ntasks_io-1) then if (.not. paranc) then if (nproc == 0) t1 = mpi_wtime() - if (global_2mDA) then - call readgriddata_2mDA(nanal1(nproc),nanal2(nproc),cvars2d,nc2d,ncdim,nbackgrounds, & - fgsfcfileprefixes,reducedgrid,grdin) - else - call readgriddata(nanal1(nproc),nanal2(nproc),cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,nbackgrounds, & - fgfileprefixes,fgsfcfileprefixes,reducedgrid,grdin,qsat) - endif + call readgriddata(nanal1(nproc),nanal2(nproc),cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,nbackgrounds, & + fgfileprefixes,fgsfcfileprefixes,reducedgrid,grdin,qsat) end if !print *,'min/max qsat',nanal,'=',minval(qsat),maxval(qsat) - if (use_qsatensmean .and. .not. global_2mDA) then + q_ind = getindex(cvars3d, 'q') + if (use_qsatensmean .and. q_ind>0 ) then ! flag for if q is in control vector allocate(qsatmean(npts,nlevs,nbackgrounds)) allocate(qsat_tmp(npts)) ! compute ensemble mean qsat @@ -268,8 +265,6 @@ subroutine read_control() ! print *,'min/max qsatmean proc',nproc,'=',& ! minval(qsatmean(:,:,nbackgrounds/2+1)),maxval(qsatmean(:,:,nbackgrounds/2+1)) !endif - if (.not. global_2mDA) then - q_ind = getindex(cvars3d, 'q') if (pseudo_rh .and. q_ind > 0) then if (use_qsatensmean) then do ne=1,nanals_per_iotask @@ -289,7 +284,6 @@ subroutine read_control() enddo endif end if - end if endif @@ -345,7 +339,6 @@ subroutine write_control(no_inflate_flag) 100 format('ens. mean anal. increment min/max ',a,2x,g19.12,2x,g19.12) deallocate(grdin_mean_tmp) - if (.not. global_2mDA) then q_ind = getindex(cvars3d, 'q') if (pseudo_rh .and. q_ind > 0) then if (use_qsatensmean) then @@ -374,16 +367,16 @@ subroutine write_control(no_inflate_flag) enddo endif end if - end if + if (.not. paranc) then if (write_fv3_incr) then - if (global_2mDA) then + if (global_2mDA) then ! flag for writing sfc increments call writeincrement_2mDA(nanal1(nproc),nanal2(nproc),cvars2d,nc2d,ncdim,grdin,no_inflate_flag) else call writeincrement(nanal1(nproc),nanal2(nproc),cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin,no_inflate_flag) endif else - if (global_2mDA) then + if (global_2mDA) then ! flag for witing sfc analyis call writegriddata_2mDA(nanal1(nproc),nanal2(nproc),cvars2d,nc2d,ncdim,grdin,no_inflate_flag) else call writegriddata(nanal1(nproc),nanal2(nproc),cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin,no_inflate_flag) @@ -393,13 +386,13 @@ subroutine write_control(no_inflate_flag) if (write_ensmean) then ! also write out ens mean on root task. if (write_fv3_incr) then - if (global_2mDA) then + if (global_2mDA) then ! flag for writing sfc inctements call writeincrement_2mDA(0,0,cvars2d,nc2d,ncdim,grdin_mean,no_inflate_flag) - else + else call writeincrement(0,0,cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin_mean,no_inflate_flag) endif else - if (global_2mDA) then + if (global_2mDA) then ! flag for writing sfc analysis call writegriddata_2mDA(0,0,cvars2d,nc2d,ncdim,grdin_mean,no_inflate_flag) else call writegriddata(0,0,cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin_mean,no_inflate_flag) @@ -415,7 +408,7 @@ subroutine write_control(no_inflate_flag) end if ! io task if (paranc) then - if (global_2mDA) then + if (global_2mDA) then ! flag for parralell write print *,'paranc not supported for global_2mDA' call mpi_barrier(mpi_comm_world,ierr) call mpi_finalize(ierr) diff --git a/src/enkf/gridinfo_gfs.f90 b/src/enkf/gridinfo_gfs.f90 index cce5b6993d..10e8369f8e 100644 --- a/src/enkf/gridinfo_gfs.f90 +++ b/src/enkf/gridinfo_gfs.f90 @@ -45,7 +45,7 @@ module gridinfo use mpisetup, only: nproc, mpi_integer, mpi_real4 use mpimod, only: mpi_comm_world -use params, only: datapath,nlevs,nlons,nlats,use_gfs_nemsio,use_gfs_ncio,fgfileprefixes, global_2mDA +use params, only: datapath,nlevs,nlons,nlats,use_gfs_nemsio,use_gfs_ncio,fgfileprefixes use kinds, only: r_kind, i_kind, r_double, r_single use constants, only: one,zero,pi,cp,rd,grav,rearth,max_varname_length use specmod, only: sptezv_s, sptez_s, init_spec_vars, isinitialized, asin_gaulats, & @@ -138,11 +138,7 @@ subroutine getgridinfo(fileprefix, reducedgrid) dset = open_dataset(filename) londim = get_dim(dset,'grid_xt'); nlonsin = londim%len latdim = get_dim(dset,'grid_yt'); nlatsin = latdim%len - if (.not. global_2mDA) then - levdim = get_dim(dset,'pfull'); nlevsin = levdim%len - else - nlevsin = 1 - endif + levdim = get_dim(dset,'pfull'); nlevsin = levdim%len idvc = 2; ntrunc = nlatsin-2 if (nlons /= nlonsin .or. nlats /= nlatsin .or. nlevs /= nlevsin) then print *,'incorrect dims in netcdf file' @@ -233,13 +229,8 @@ subroutine getgridinfo(fileprefix, reducedgrid) endif ! convert to 1d array, units to millibars, flip so lats go N to S. spressmn = 0.01_r_kind*reshape(values_2d,(/nlons*nlats/)) - if (.not. global_2mDA) then - call read_attribute(dset, 'ak', ak) - call read_attribute(dset, 'bk', bk) - else - allocate(ak(nlevs+1),bk(nlevs+1)) - ak=0; bk=1 - endif + call read_attribute(dset, 'ak', ak) + call read_attribute(dset, 'bk', bk) call close_dataset(dset) ! pressure at interfaces do k=1,nlevs+1 @@ -311,15 +302,11 @@ subroutine getgridinfo(fileprefix, reducedgrid) enddo enddo endif - if (global_2mDA) then - presslmn(:,nlevs) = spressmn - else - do k=1,nlevs - ! layer pressure from Phillips vertical interpolation. - presslmn(:,k) = ((pressimn(:,k)**kap1-pressimn(:,k+1)**kap1)/& + do k=1,nlevs + ! layer pressure from Phillips vertical interpolation. + presslmn(:,k) = ((pressimn(:,k)**kap1-pressimn(:,k+1)**kap1)/& (kap1*(pressimn(:,k)-pressimn(:,k+1))))**kapr - end do - endif + end do print *,'ensemble mean first guess surface pressure:' print *,minval(spressmn),maxval(spressmn) !do k=1,nlevs diff --git a/src/enkf/gridio_gfs.f90 b/src/enkf/gridio_gfs.f90 index c5d66a7642..1b882e272f 100644 --- a/src/enkf/gridio_gfs.f90 +++ b/src/enkf/gridio_gfs.f90 @@ -35,6 +35,7 @@ module gridio ! a required input for EFSO calculations ! 2019-03-13 Add precipitation components ! 2019-07-10 Add convective clouds +! 2022-07-21 Draper: added reading vars from sfc file for nc case. ! ! attributes: ! language: f95 @@ -55,7 +56,7 @@ module gridio implicit none private public :: readgriddata, readgriddata_pnc, writegriddata, writegriddata_pnc - public :: readgriddata_2mDA,writeincrement, writeincrement_pnc, & + public :: writeincrement, writeincrement_pnc, & writegriddata_2mDA, writeincrement_2mDA contains @@ -416,187 +417,6 @@ end subroutine copytogrdin end subroutine readgriddata_pnc - subroutine readgriddata_2mDA(nanal1,nanal2,vars2d,n2d,ndim,ntimes, & - fileprefixes,reducedgrid,grdin) - use module_fv3gfs_ncio, only: Dataset, Variable, Dimension, open_dataset,& - quantize_data, read_attribute, close_dataset, get_dim, read_vardata - implicit none - - integer, intent(in) :: nanal1,nanal2 - character(len=max_varname_length), dimension(n2d), intent(in) :: vars2d - integer, intent(in) :: n2d - integer, intent(in) :: ndim, ntimes - character(len=120), dimension(7), intent(in) :: fileprefixes - logical, intent(in) :: reducedgrid - real(r_single), dimension(npts,ndim,ntimes,nanal2-nanal1+1), intent(out) :: grdin - - character(len=500) :: filename - character(len=7) charnanal - - real(r_single), allocatable, dimension(:,:) :: values_2d - real(r_kind), dimension(nlons*nlats) :: ug - type(Dataset) :: dset - type(Dimension) :: londim,latdim - - 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) :: iunitsig,iret,nb,nlonsin,nlatsin,ne,nanal - - if (reducedgrid) then - print *,"reducedgrid=T not valid for global_2mDA=T" - call stop2(22) - endif - - ne = 0 - ensmemloop: do nanal=nanal1,nanal2 - ne = ne + 1 - backgroundloop: do nb=1,ntimes - - if (nanal > 0) then - write(charnanal,'(a3, i3.3)') 'mem', nanal - else - charnanal = 'ensmean' - endif - iunitsig = 77 - filename = trim(adjustl(datapath))//trim(adjustl(fileprefixes(nb)))//trim(charnanal) - dset = open_dataset(filename) - londim = get_dim(dset,'grid_xt'); nlonsin = londim%len - latdim = get_dim(dset,'grid_yt'); nlatsin = latdim%len - - tmp2m_ind = getindex(vars2d, 't2m') !< indices in the state or control var arrays - spfh2m_ind = getindex(vars2d, 'q2m') - soilt1_ind = getindex(vars2d, 'soilt1') - slc1_ind = getindex(vars2d, 'slc1') - soilt2_ind = getindex(vars2d, 'soilt2') - slc2_ind = getindex(vars2d, 'slc2') - soilt3_ind = getindex(vars2d, 'soilt3') - slc3_ind = getindex(vars2d, 'slc3') - soilt4_ind = getindex(vars2d, 'soilt4') - slc4_ind = getindex(vars2d, 'slc4') - -! if (nproc == 0) then -! print *, 'indices: ' -! print *, 't2m_ind: ', t2m_ind, ', q2m_ind: ', q2m_ind, ', soilt_ind: ', soilt_ind, ', soilm_ind: ',soilm_ind -! endif - - if (tmp2m_ind > 0) then - call read_vardata(dset, 'tmp2m', values_2d, errcode=iret) - if (iret /= 0) then - print *,'error reading tmp2m' - call stop2(22) - endif - ug = reshape(values_2d,(/nlons*nlats/)) - call copytogrdin(ug,grdin(:,tmp2m_ind,nb,ne)) - endif - if (spfh2m_ind > 0) then - call read_vardata(dset, 'spfh2m', values_2d, errcode=iret) - if (iret /= 0) then - print *,'error reading spfh2m' - call stop2(22) - endif - ug = reshape(values_2d,(/nlons*nlats/)) - call copytogrdin(ug,grdin(:,spfh2m_ind,nb,ne)) - endif - if (soilt1_ind > 0) then - call read_vardata(dset, 'soilt1', values_2d, errcode=iret) - if (iret /= 0) then - print *,'error reading soilt1' - call stop2(22) - endif - ug = reshape(values_2d,(/nlons*nlats/)) - call copytogrdin(ug,grdin(:,soilt1_ind,nb,ne)) - endif - if (soilt2_ind > 0) then - call read_vardata(dset, 'soilt2', values_2d, errcode=iret) - if (iret /= 0) then - print *,'error reading soilt2' - call stop2(22) - endif - ug = reshape(values_2d,(/nlons*nlats/)) - call copytogrdin(ug,grdin(:,soilt2_ind,nb,ne)) - endif - if (soilt3_ind > 0) then - call read_vardata(dset, 'soilt3', values_2d, errcode=iret) - if (iret /= 0) then - print *,'error reading soilt3' - call stop2(22) - endif - ug = reshape(values_2d,(/nlons*nlats/)) - call copytogrdin(ug,grdin(:,soilt3_ind,nb,ne)) - endif - if (soilt4_ind > 0) then - call read_vardata(dset, 'soilt4', values_2d, errcode=iret) - if (iret /= 0) then - print *,'error reading soilt2' - call stop2(22) - endif - ug = reshape(values_2d,(/nlons*nlats/)) - call copytogrdin(ug,grdin(:,soilt4_ind,nb,ne)) - endif - if (slc1_ind > 0) then - call read_vardata(dset, 'slc1', values_2d, errcode=iret) - if (iret /= 0) then - print *,'error reading slc1' - call stop2(22) - endif - ug = reshape(values_2d,(/nlons*nlats/)) - call copytogrdin(ug,grdin(:,slc1_ind,nb,ne)) - endif - if (slc2_ind > 0) then - call read_vardata(dset, 'slc2', values_2d, errcode=iret) - if (iret /= 0) then - print *,'error reading slc2' - call stop2(22) - endif - ug = reshape(values_2d,(/nlons*nlats/)) - call copytogrdin(ug,grdin(:,slc2_ind,nb,ne)) - endif - if (slc3_ind > 0) then - call read_vardata(dset, 'slc3', values_2d, errcode=iret) - if (iret /= 0) then - print *,'error reading slc3' - call stop2(22) - endif - ug = reshape(values_2d,(/nlons*nlats/)) - call copytogrdin(ug,grdin(:,slc3_ind,nb,ne)) - endif - if (slc4_ind > 0) then - call read_vardata(dset, 'slc4', values_2d, errcode=iret) - if (iret /= 0) then - print *,'error reading slc2' - call stop2(22) - endif - ug = reshape(values_2d,(/nlons*nlats/)) - call copytogrdin(ug,grdin(:,slc4_ind,nb,ne)) - endif - - call close_dataset(dset) - - end do backgroundloop ! loop over backgrounds to read in - end do ensmemloop ! loop over ens members to read in - - return - - contains - ! copying to grdin (calling regtoreduced if reduced grid) - subroutine copytogrdin(field, grdin) - implicit none - - real(r_kind), dimension(:), intent(in) :: field - real(r_single), dimension(:), intent(inout) :: grdin - - if (reducedgrid) then - call regtoreduced(field, grdin) - else - grdin = field - endif - - end subroutine copytogrdin - - end subroutine readgriddata_2mDA - - subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & fileprefixes,filesfcprefixes,reducedgrid,grdin,qsat) use sigio_module, only: sigio_head, sigio_data, sigio_sclose, sigio_sropen, & @@ -638,7 +458,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, type(sigio_head) :: sighead type(sigio_data) :: sigdata type(nemsio_gfile) :: gfile - type(Dataset) :: dset + type(Dataset) :: dset, dset_sfc type(Dimension) :: londim,latdim,levdim type(nemsio_gfile) :: gfilesfc @@ -646,14 +466,68 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, 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 + 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,iunitsig,iret,nb,i,idvc,nlonsin,nlatsin,nlevsin,ne,nanal integer(i_kind) :: nlonsin_sfc,nlatsin_sfc logical ice logical use_full_hydro + logical read_sfc_file, read_atm_file use_full_hydro = .false. + u_ind = getindex(vars3d, 'u') !< indices in the state or control var arrays + v_ind = getindex(vars3d, 'v') ! U and V (3D) + tv_ind = getindex(vars3d, 'tv') ! Tv (3D) + q_ind = getindex(vars3d, 'q') ! Q (3D) + oz_ind = getindex(vars3d, 'oz') ! Oz (3D) + cw_ind = getindex(vars3d, 'cw') ! CW (3D) + tsen_ind = getindex(vars3d, 'tsen') !sensible T (3D) + ql_ind = getindex(vars3d, 'ql') ! QL (3D) + qi_ind = getindex(vars3d, 'qi') ! QI (3D) + prse_ind = getindex(vars3d, 'prse') + qr_ind = getindex(vars3d, 'qr') ! QR (3D) + qs_ind = getindex(vars3d, 'qs') ! QS (3D) + qg_ind = getindex(vars3d, 'qg') ! QG (3D) + ps_ind = getindex(vars2d, 'ps') ! Ps (2D) + pst_ind = getindex(vars2d, 'pst') ! Ps tendency (2D) // equivalent of + ! old logical massbal_adjust, if non-zero + sst_ind = getindex(vars2d, 'sst') + + read_atm_file = ( u_ind>0 .or. v_ind>0 .or. tv_ind>0 .or. q_ind>0 .or. & + oz_ind>0 .or. cw_ind>0 .or. tsen_ind>0 .or. ql_ind>0 .or. & + qi_ind>0 .or. prse_ind>0 .or. qr_ind>0 .or. qs_ind>0 .or. qg_ind>0 ) + + use_full_hydro = ( ql_ind > 0 .and. qi_ind > 0 .and. & + qr_ind > 0 .and. qs_ind > 0 .and. qg_ind > 0 ) + +! if (nproc == 0) then +! print *, 'indices: ' +! print *, 'u: ', u_ind, ', v: ', v_ind, ', tv: ', tv_ind, ', tsen: ', tsen_ind +! print *, 'q: ', q_ind, ', oz: ', oz_ind, ', cw: ', cw_ind, ', qi: ', qi_ind +! print *, 'ql: ', ql_ind, ', prse: ', prse_ind +! print *, 'ps: ', ps_ind, ', pst: ', pst_ind, ', sst: ', sst_ind +! endif + + ! land sfc DA variables + tmp2m_ind = getindex(vars2d, 't2m') + spfh2m_ind = getindex(vars2d, 'q2m') + soilt1_ind = getindex(vars2d, 'soilt1') + slc1_ind = getindex(vars2d, 'slc1') + soilt2_ind = getindex(vars2d, 'soilt2') + slc2_ind = getindex(vars2d, 'slc2') + soilt3_ind = getindex(vars2d, 'soilt3') + slc3_ind = getindex(vars2d, 'slc3') + soilt4_ind = getindex(vars2d, 'soilt4') + slc4_ind = getindex(vars2d, 'slc4') + + read_sfc_file = ( tmp2m_ind > 0 .or. spfh2m_ind > 0 .or. soilt1_ind > 0 .or. & + slc1_ind > 0 .or. soilt2_ind > 0 .or. slc2_ind > 0 .or. & + soilt3_ind > 0 .or. slc3_ind > 0 .or. soilt4_ind > 0 .or. & + slc4_ind > 0 ) + + ne = 0 ensmemloop: do nanal=nanal1,nanal2 ne = ne + 1 @@ -701,6 +575,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, endif endif else if (use_gfs_ncio) then + ! for read_sfc_file assuming atmos file available dset = open_dataset(filename) londim = get_dim(dset,'grid_xt'); nlonsin = londim%len latdim = get_dim(dset,'grid_yt'); nlatsin = latdim%len @@ -714,480 +589,569 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, call stop2(23) end if endif - ice = .false. ! calculate qsat w/resp to ice? - kap = rd/cp - kapr = cp/rd - kap1 = kap+one - - u_ind = getindex(vars3d, 'u') !< indices in the state or control var arrays - v_ind = getindex(vars3d, 'v') ! U and V (3D) - tv_ind = getindex(vars3d, 'tv') ! Tv (3D) - q_ind = getindex(vars3d, 'q') ! Q (3D) - oz_ind = getindex(vars3d, 'oz') ! Oz (3D) - cw_ind = getindex(vars3d, 'cw') ! CW (3D) - tsen_ind = getindex(vars3d, 'tsen') !sensible T (3D) - ql_ind = getindex(vars3d, 'ql') ! QL (3D) - qi_ind = getindex(vars3d, 'qi') ! QI (3D) - prse_ind = getindex(vars3d, 'prse') - qr_ind = getindex(vars3d, 'qr') ! QR (3D) - qs_ind = getindex(vars3d, 'qs') ! QS (3D) - qg_ind = getindex(vars3d, 'qg') ! QG (3D) - ps_ind = getindex(vars2d, 'ps') ! Ps (2D) - pst_ind = getindex(vars2d, 'pst') ! Ps tendency (2D) // equivalent of - ! old logical massbal_adjust, if non-zero - sst_ind = getindex(vars2d, 'sst') - use_full_hydro = ( ql_ind > 0 .and. qi_ind > 0 .and. & - qr_ind > 0 .and. qs_ind > 0 .and. qg_ind > 0 ) -! if (nproc == 0) then -! print *, 'indices: ' -! print *, 'u: ', u_ind, ', v: ', v_ind, ', tv: ', tv_ind, ', tsen: ', tsen_ind -! print *, 'q: ', q_ind, ', oz: ', oz_ind, ', cw: ', cw_ind, ', qi: ', qi_ind -! print *, 'ql: ', ql_ind, ', prse: ', prse_ind -! print *, 'ps: ', ps_ind, ', pst: ', pst_ind, ', sst: ', sst_ind -! endif + if (read_atm_file) then - if (.not. isinitialized) call init_spec_vars(nlons,nlats,ntrunc,4) + ice = .false. ! calculate qsat w/resp to ice? + kap = rd/cp + kapr = cp/rd + kap1 = kap+one - allocate(pressi(nlons*nlats,nlevs+1)) - allocate(pslg(npts,nlevs)) - allocate(psg(nlons*nlats)) - if (pst_ind > 0) allocate(vmassdiv(nlons*nlats,nlevs),pstend(nlons*nlats)) - if (use_gfs_nemsio) then - call nemsio_readrecv(gfile,'pres','sfc',1,nems_wrk,iret=iret) - if (iret/=0) then - write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(ps), iret=',iret - call stop2(23) - endif - psg = 0.01_r_kind*nems_wrk ! convert ps to millibars. - - if (allocated(nems_vcoord)) deallocate(nems_vcoord) - allocate(nems_vcoord(nlevs+1,3,2)) - call nemsio_getfilehead(gfile,iret=iret,vcoord=nems_vcoord) - if ( iret /= 0 ) then - write(6,*)' gridio: ***ERROR*** problem reading header ', & - 'vcoord, Status = ',iret - call stop2(99) - endif + if (.not. isinitialized) call init_spec_vars(nlons,nlats,ntrunc,4) - allocate(ak(nlevs+1),bk(nlevs+1)) - - if ( idvc == 0 ) then ! sigma coordinate, old file format. - ak = zero - bk = nems_vcoord(1:nlevs+1,1,1) - elseif ( idvc == 1 ) then ! sigma coordinate - ak = zero - bk = nems_vcoord(1:nlevs+1,2,1) - elseif ( idvc == 2 .or. idvc == 3 ) then ! hybrid coordinate - ak = 0.01_r_kind*nems_vcoord(1:nlevs+1,1,1) ! convert to mb - bk = nems_vcoord(1:nlevs+1,2,1) - else - write(6,*)'gridio: ***ERROR*** INVALID value for idvc=',idvc - call stop2(85) - endif - if (nanal .eq. 1) then - print *,'time level ',nb - print *,'---------------' - endif - ! pressure at interfaces - do k=1,nlevs+1 - pressi(:,k)=ak(k)+bk(k)*psg - if (nanal .eq. 1) print *,'nemsio, min/max pressi',k,minval(pressi(:,k)),maxval(pressi(:,k)) - enddo - deallocate(ak,bk) - else if (use_gfs_ncio) then - call read_vardata(dset, 'pressfc', values_2d,errcode=iret) - if (iret /= 0) then - print *,'error reading ps' - call stop2(31) - endif - psg = 0.01_r_kind*reshape(values_2d,(/nlons*nlats/)) - call read_attribute(dset, 'ak', ak) - call read_attribute(dset, 'bk', bk) - if (nanal .eq. 1) then - print *,'time level ',nb - print *,'---------------' - endif - ! pressure at interfaces - do k=1,nlevs+1 - ! k=1 in ak,bk is model top - pressi(:,k) = 0.01_r_kind*ak(nlevs-k+2)+bk(nlevs-k+2)*psg - if (nanal .eq. 1) print *,'netcdf, min/max pressi',k,minval(pressi(:,k)),maxval(pressi(:,k)) - enddo - deallocate(ak,bk,values_2d) - else - vrtspec = sigdata%ps - call sptez_s(vrtspec,psg,1) - !==> input psg is ln(ps) in centibars - convert to ps in millibars. - psg = 10._r_kind*exp(psg) - allocate(ak(nlevs+1),bk(nlevs+1)) - if (sighead%idvc .eq. 0) then ! sigma coordinate, old file format. - ak = zero - bk = sighead%si(1:nlevs+1) - else if (sighead%idvc == 1) then ! sigma coordinate - ak = zero - bk = sighead%vcoord(1:nlevs+1,2) - else if (sighead%idvc == 2 .or. sighead%idvc == 3) then ! hybrid coordinate - bk = sighead%vcoord(1:nlevs+1,2) - ak = 0.01_r_kind*sighead%vcoord(1:nlevs+1,1) ! convert to mb - else - print *,'unknown vertical coordinate type',sighead%idvc - call stop2(32) - end if - !==> pressure at interfaces. - if (nanal .eq. 1) then - print *,'time level ',nb - print *,'---------------' - endif - do k=1,nlevs+1 - pressi(:,k)=ak(k)+bk(k)*psg - if (nanal .eq. 1) print *,'sigio, min/max pressi',k,minval(pressi(:,k)),maxval(pressi(:,k)) - enddo - deallocate(ak,bk) - endif + allocate(pressi(nlons*nlats,nlevs+1)) + allocate(pslg(npts,nlevs)) + allocate(psg(nlons*nlats)) + if (pst_ind > 0) allocate(vmassdiv(nlons*nlats,nlevs),pstend(nlons*nlats)) - !==> get U,V,temp,q,ps on gaussian grid. - ! u is first nlevs, v is second, t is third, then tracers. - if (use_gfs_nemsio) then - clip=tiny_r_kind - do k=1,nlevs - call nemsio_readrecv(gfile,'ugrd','mid layer',k,nems_wrk,iret=iret) + if (use_gfs_nemsio) then + call nemsio_readrecv(gfile,'pres','sfc',1,nems_wrk,iret=iret) if (iret/=0) then - write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(ugrd), iret=',iret + write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(ps), iret=',iret call stop2(23) endif - ug = nems_wrk - call nemsio_readrecv(gfile,'vgrd','mid layer',k,nems_wrk,iret=iret) - if (iret/=0) then - write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(vgrd), iret=',iret - call stop2(23) + psg = 0.01_r_kind*nems_wrk ! convert ps to millibars. + + if (allocated(nems_vcoord)) deallocate(nems_vcoord) + allocate(nems_vcoord(nlevs+1,3,2)) + call nemsio_getfilehead(gfile,iret=iret,vcoord=nems_vcoord) + if ( iret /= 0 ) then + write(6,*)' gridio: ***ERROR*** problem reading header ', & + 'vcoord, Status = ',iret + call stop2(99) endif - vg = nems_wrk - if (u_ind > 0) call copytogrdin(ug,grdin(:,levels(u_ind-1) + k,nb,ne)) - if (v_ind > 0) call copytogrdin(vg,grdin(:,levels(v_ind-1) + k,nb,ne)) - ! calculate vertical integral of mass flux div (ps tendency) - ! this variable is analyzed in order to enforce mass balance in the analysis - if (pst_ind > 0) then - ug = ug*(pressi(:,k)-pressi(:,k+1)) - vg = vg*(pressi(:,k)-pressi(:,k+1)) - call sptezv_s(divspec,vrtspec,ug,vg,-1) ! u,v to div,vrt - call sptez_s(divspec,vmassdiv(:,k),1) ! divspec to divgrd + + allocate(ak(nlevs+1),bk(nlevs+1)) + + if ( idvc == 0 ) then ! sigma coordinate, old file format. + ak = zero + bk = nems_vcoord(1:nlevs+1,1,1) + elseif ( idvc == 1 ) then ! sigma coordinate + ak = zero + bk = nems_vcoord(1:nlevs+1,2,1) + elseif ( idvc == 2 .or. idvc == 3 ) then ! hybrid coordinate + ak = 0.01_r_kind*nems_vcoord(1:nlevs+1,1,1) ! convert to mb + bk = nems_vcoord(1:nlevs+1,2,1) + else + write(6,*)'gridio: ***ERROR*** INVALID value for idvc=',idvc + call stop2(85) endif - call nemsio_readrecv(gfile,'tmp','mid layer',k,nems_wrk,iret=iret) - if (iret/=0) then - write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(tmp), iret=',iret - call stop2(23) + if (nanal .eq. 1) then + print *,'time level ',nb + print *,'---------------' endif - call nemsio_readrecv(gfile,'spfh','mid layer',k,nems_wrk2,iret=iret) - if (iret/=0) then - write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(spfh), iret=',iret - call stop2(23) + ! pressure at interfaces + do k=1,nlevs+1 + pressi(:,k)=ak(k)+bk(k)*psg + if (nanal .eq. 1) print *,'nemsio, min/max pressi',k,minval(pressi(:,k)),maxval(pressi(:,k)) + enddo + deallocate(ak,bk) + else if (use_gfs_ncio) then + call read_vardata(dset, 'pressfc', values_2d,errcode=iret) + if (iret /= 0) then + print *,'error reading ps' + call stop2(31) endif - if (cliptracers) where (nems_wrk2 < clip) nems_wrk2 = clip - ug = nems_wrk - if (tsen_ind > 0) call copytogrdin(ug,grdin(:,levels(tsen_ind-1)+k,nb,ne)) - nems_wrk = nems_wrk * ( 1.0 + fv*nems_wrk2 ) ! convert T to Tv - ug = nems_wrk - vg = nems_wrk2 - call copytogrdin(ug,tv(:,k)) - call copytogrdin(vg, q(:,k)) - if (tv_ind > 0) grdin(:,levels(tv_ind-1)+k,nb,ne) = tv(:,k) - if (q_ind > 0) grdin(:,levels( q_ind-1)+k,nb,ne) = q(:,k) - if (oz_ind > 0) then - call nemsio_readrecv(gfile,'o3mr','mid layer',k,nems_wrk2,iret=iret) - if (iret/=0) then - write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(o3mr), iret=',iret - call stop2(23) - endif - if (cliptracers) where (nems_wrk2 < clip) nems_wrk2 = clip - ug = nems_wrk2 - call copytogrdin(ug,grdin(:,levels(oz_ind-1)+k,nb,ne)) + psg = 0.01_r_kind*reshape(values_2d,(/nlons*nlats/)) + call read_attribute(dset, 'ak', ak) + call read_attribute(dset, 'bk', bk) + if (nanal .eq. 1) then + print *,'time level ',nb + print *,'---------------' endif - if (.not. use_full_hydro) then - if (cw_ind > 0 .or. ql_ind > 0 .or. qi_ind > 0) then - call nemsio_readrecv(gfile,'clwmr','mid layer',k,nems_wrk2,iret=iret) - if (iret/=0) then - write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(clwmr), iret=',iret + ! pressure at interfaces + do k=1,nlevs+1 + ! k=1 in ak,bk is model top + pressi(:,k) = 0.01_r_kind*ak(nlevs-k+2)+bk(nlevs-k+2)*psg + if (nanal .eq. 1) print *,'netcdf, min/max pressi',k,minval(pressi(:,k)),maxval(pressi(:,k)) + enddo + deallocate(ak,bk) + else + vrtspec = sigdata%ps + call sptez_s(vrtspec,psg,1) + !==> input psg is ln(ps) in centibars - convert to ps in millibars. + psg = 10._r_kind*exp(psg) + allocate(ak(nlevs+1),bk(nlevs+1)) + if (sighead%idvc .eq. 0) then ! sigma coordinate, old file format. + ak = zero + bk = sighead%si(1:nlevs+1) + else if (sighead%idvc == 1) then ! sigma coordinate + ak = zero + bk = sighead%vcoord(1:nlevs+1,2) + else if (sighead%idvc == 2 .or. sighead%idvc == 3) then ! hybrid coordinate + bk = sighead%vcoord(1:nlevs+1,2) + ak = 0.01_r_kind*sighead%vcoord(1:nlevs+1,1) ! convert to mb + else + print *,'unknown vertical coordinate type',sighead%idvc + call stop2(32) + end if + !==> pressure at interfaces. + if (nanal .eq. 1) then + print *,'time level ',nb + print *,'---------------' + endif + do k=1,nlevs+1 + pressi(:,k)=ak(k)+bk(k)*psg + if (nanal .eq. 1) print *,'sigio, min/max pressi',k,minval(pressi(:,k)),maxval(pressi(:,k)) + enddo + deallocate(ak,bk) + endif + + !==> get U,V,temp,q,ps on gaussian grid. + ! u is first nlevs, v is second, t is third, then tracers. + if (use_gfs_nemsio) then + clip=tiny_r_kind + do k=1,nlevs + call nemsio_readrecv(gfile,'ugrd','mid layer',k,nems_wrk,iret=iret) + if (iret/=0) then + write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(ugrd), iret=',iret + call stop2(23) + endif + ug = nems_wrk + call nemsio_readrecv(gfile,'vgrd','mid layer',k,nems_wrk,iret=iret) + if (iret/=0) then + write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(vgrd), iret=',iret + call stop2(23) + endif + vg = nems_wrk + if (u_ind > 0) call copytogrdin(ug,grdin(:,levels(u_ind-1) + k,nb,ne)) + if (v_ind > 0) call copytogrdin(vg,grdin(:,levels(v_ind-1) + k,nb,ne)) + ! calculate vertical integral of mass flux div (ps tendency) + ! this variable is analyzed in order to enforce mass balance in the analysis + if (pst_ind > 0) then + ug = ug*(pressi(:,k)-pressi(:,k+1)) + vg = vg*(pressi(:,k)-pressi(:,k+1)) + call sptezv_s(divspec,vrtspec,ug,vg,-1) ! u,v to div,vrt + call sptez_s(divspec,vmassdiv(:,k),1) ! divspec to divgrd + endif + call nemsio_readrecv(gfile,'tmp','mid layer',k,nems_wrk,iret=iret) + if (iret/=0) then + write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(tmp), iret=',iret + call stop2(23) + endif + call nemsio_readrecv(gfile,'spfh','mid layer',k,nems_wrk2,iret=iret) + if (iret/=0) then + write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(spfh), iret=',iret + call stop2(23) + endif + if (cliptracers) where (nems_wrk2 < clip) nems_wrk2 = clip + ug = nems_wrk + if (tsen_ind > 0) call copytogrdin(ug,grdin(:,levels(tsen_ind-1)+k,nb,ne)) + nems_wrk = nems_wrk * ( 1.0 + fv*nems_wrk2 ) ! convert T to Tv + ug = nems_wrk + vg = nems_wrk2 + call copytogrdin(ug,tv(:,k)) + call copytogrdin(vg, q(:,k)) + if (tv_ind > 0) grdin(:,levels(tv_ind-1)+k,nb,ne) = tv(:,k) + if (q_ind > 0) grdin(:,levels( q_ind-1)+k,nb,ne) = q(:,k) + if (oz_ind > 0) then + call nemsio_readrecv(gfile,'o3mr','mid layer',k,nems_wrk2,iret=iret) + if (iret/=0) then + write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(o3mr), iret=',iret call stop2(23) endif - if (imp_physics == 11) then - call nemsio_readrecv(gfile,'icmr','mid layer',k,nems_wrk,iret=iret) + if (cliptracers) where (nems_wrk2 < clip) nems_wrk2 = clip + ug = nems_wrk2 + call copytogrdin(ug,grdin(:,levels(oz_ind-1)+k,nb,ne)) + endif + if (.not. use_full_hydro) then + if (cw_ind > 0 .or. ql_ind > 0 .or. qi_ind > 0) then + call nemsio_readrecv(gfile,'clwmr','mid layer',k,nems_wrk2,iret=iret) if (iret/=0) then - write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(icmr), iret=',iret + write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(clwmr), iret=',iret call stop2(23) - else - nems_wrk2 = nems_wrk2 + nems_wrk endif + if (imp_physics == 11) then + call nemsio_readrecv(gfile,'icmr','mid layer',k,nems_wrk,iret=iret) + if (iret/=0) then + write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(icmr), iret=',iret + call stop2(23) + else + nems_wrk2 = nems_wrk2 + nems_wrk + endif + endif + if (cnvw_option) then + call nemsio_readrecv(gfilesfc,'cnvcldwat','mid layer',k,nems_wrk,iret=iret) + if (iret/=0) then + write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(cnvw), iret=',iret + call stop2(23) + else + nems_wrk2 = nems_wrk2 + nems_wrk + end if + end if + if (cliptracers) where (nems_wrk2 < clip) nems_wrk2 = clip + ug = nems_wrk2 + call copytogrdin(ug,cw(:,k)) + if (cw_ind > 0) grdin(:,levels(cw_ind-1)+k,nb,ne) = cw(:,k) endif - if (cnvw_option) then - call nemsio_readrecv(gfilesfc,'cnvcldwat','mid layer',k,nems_wrk,iret=iret) + else + if (ql_ind > 0) then + call nemsio_readrecv(gfile,'clwmr','mid layer',k,nems_wrk,iret=iret) if (iret/=0) then - write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(cnvw), iret=',iret + write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(clwmr), iret=',iret call stop2(23) - else - nems_wrk2 = nems_wrk2 + nems_wrk - end if - end if - if (cliptracers) where (nems_wrk2 < clip) nems_wrk2 = clip - ug = nems_wrk2 - call copytogrdin(ug,cw(:,k)) - if (cw_ind > 0) grdin(:,levels(cw_ind-1)+k,nb,ne) = cw(:,k) - endif - else - if (ql_ind > 0) then - call nemsio_readrecv(gfile,'clwmr','mid layer',k,nems_wrk,iret=iret) - if (iret/=0) then - write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(clwmr), iret=',iret - call stop2(23) + endif + if (cliptracers) where (nems_wrk < clip) nems_wrk = clip + ug = nems_wrk + call copytogrdin(ug,ql(:,k)) + grdin(:,levels(ql_ind-1)+k,nb,ne) = ql(:,k) endif - if (cliptracers) where (nems_wrk < clip) nems_wrk = clip - ug = nems_wrk - call copytogrdin(ug,ql(:,k)) - grdin(:,levels(ql_ind-1)+k,nb,ne) = ql(:,k) - endif - if (qi_ind > 0) then - call nemsio_readrecv(gfile,'icmr','mid layer',k,nems_wrk,iret=iret) - if (iret/=0) then - write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(icmr), iret=',iret - call stop2(23) + if (qi_ind > 0) then + call nemsio_readrecv(gfile,'icmr','mid layer',k,nems_wrk,iret=iret) + if (iret/=0) then + write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(icmr), iret=',iret + call stop2(23) + endif + if (cliptracers) where (nems_wrk < clip) nems_wrk = clip + ug = nems_wrk + call copytogrdin(ug,qi(:,k)) + grdin(:,levels(qi_ind-1)+k,nb,ne) = qi(:,k) endif - if (cliptracers) where (nems_wrk < clip) nems_wrk = clip - ug = nems_wrk - call copytogrdin(ug,qi(:,k)) - grdin(:,levels(qi_ind-1)+k,nb,ne) = qi(:,k) - endif - if (qr_ind > 0) then - call nemsio_readrecv(gfile,'rwmr','mid layer',k,nems_wrk,iret=iret) - if (iret/=0) then - write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(rwmr), iret=',iret - call stop2(23) + if (qr_ind > 0) then + call nemsio_readrecv(gfile,'rwmr','mid layer',k,nems_wrk,iret=iret) + if (iret/=0) then + write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(rwmr), iret=',iret + call stop2(23) + endif + if (cliptracers) where (nems_wrk < clip) nems_wrk = clip + ug = nems_wrk + call copytogrdin(ug,qr(:,k)) + grdin(:,levels(qr_ind-1)+k,nb,ne) = qr(:,k) endif - if (cliptracers) where (nems_wrk < clip) nems_wrk = clip - ug = nems_wrk - call copytogrdin(ug,qr(:,k)) - grdin(:,levels(qr_ind-1)+k,nb,ne) = qr(:,k) - endif - if (qs_ind > 0) then - call nemsio_readrecv(gfile,'snmr','mid layer',k,nems_wrk,iret=iret) - if (iret/=0) then - write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(snmr), iret=',iret - call stop2(23) + if (qs_ind > 0) then + call nemsio_readrecv(gfile,'snmr','mid layer',k,nems_wrk,iret=iret) + if (iret/=0) then + write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(snmr), iret=',iret + call stop2(23) + endif + if (cliptracers) where (nems_wrk < clip) nems_wrk = clip + ug = nems_wrk + call copytogrdin(ug,qs(:,k)) + grdin(:,levels(qs_ind-1)+k,nb,ne) = qs(:,k) endif - if (cliptracers) where (nems_wrk < clip) nems_wrk = clip - ug = nems_wrk - call copytogrdin(ug,qs(:,k)) - grdin(:,levels(qs_ind-1)+k,nb,ne) = qs(:,k) - endif - if (qg_ind > 0) then - call nemsio_readrecv(gfile,'grle','mid layer',k,nems_wrk,iret=iret) - if (iret/=0) then - write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(grle), iret=',iret - call stop2(23) + if (qg_ind > 0) then + call nemsio_readrecv(gfile,'grle','mid layer',k,nems_wrk,iret=iret) + if (iret/=0) then + write(6,*)'gridio/readgriddata: gfs model: problem with nemsio_readrecv(grle), iret=',iret + call stop2(23) + endif + if (cliptracers) where (nems_wrk < clip) nems_wrk = clip + ug = nems_wrk + call copytogrdin(ug,qg(:,k)) + grdin(:,levels(qg_ind-1)+k,nb,ne) = qg(:,k) endif - if (cliptracers) where (nems_wrk < clip) nems_wrk = clip - ug = nems_wrk - call copytogrdin(ug,qg(:,k)) - grdin(:,levels(qg_ind-1)+k,nb,ne) = qg(:,k) - endif - endif ! use_full_hydro - enddo - else if (use_gfs_ncio) then - call read_vardata(dset, 'ugrd', ug3d,errcode=iret) - if (iret /= 0) then - print *,'error reading ugrd' - call stop2(22) - endif - call read_vardata(dset, 'vgrd', vg3d,errcode=iret) - if (iret /= 0) then - print *,'error reading vgrd' - call stop2(23) - endif - do k=1,nlevs - ug = reshape(ug3d(:,:,nlevs-k+1),(/nlons*nlats/)) - vg = reshape(vg3d(:,:,nlevs-k+1),(/nlons*nlats/)) - if (u_ind > 0) call copytogrdin(ug,grdin(:,levels(u_ind-1) + k,nb,ne)) - if (v_ind > 0) call copytogrdin(vg,grdin(:,levels(v_ind-1) + k,nb,ne)) - ! calculate vertical integral of mass flux div (ps tendency) - ! this variable is analyzed in order to enforce mass balance in the analysis - if (pst_ind > 0) then - ug = ug*(pressi(:,k)-pressi(:,k+1)) - vg = vg*(pressi(:,k)-pressi(:,k+1)) - call sptezv_s(divspec,vrtspec,ug,vg,-1) ! u,v to div,vrt - call sptez_s(divspec,vmassdiv(:,k),1) ! divspec to divgrd + endif ! use_full_hydro + enddo + else if (use_gfs_ncio) then + call read_vardata(dset, 'ugrd', ug3d,errcode=iret) + if (iret /= 0) then + print *,'error reading ugrd' + call stop2(22) endif - enddo - call read_vardata(dset,'tmp', ug3d,errcode=iret) - if (iret /= 0) then - print *,'error reading tmp' - call stop2(24) - endif - call read_vardata(dset,'spfh', vg3d,errcode=iret) - if (iret /= 0) then - print *,'error reading spfh' - call stop2(25) - endif - do k=1,nlevs - ug = reshape(ug3d(:,:,nlevs-k+1),(/nlons*nlats/)) - vg = reshape(vg3d(:,:,nlevs-k+1),(/nlons*nlats/)) - if (tsen_ind > 0) call copytogrdin(ug,grdin(:,levels(tsen_ind-1)+k,nb,ne)) - call copytogrdin(vg, q(:,k)) - ug = ug * ( 1.0 + fv*vg ) ! convert T to Tv - call copytogrdin(ug,tv(:,k)) - if (tv_ind > 0) grdin(:,levels(tv_ind-1)+k,nb,ne) = tv(:,k) - if (q_ind > 0) grdin(:,levels( q_ind-1)+k,nb,ne) = q(:,k) - enddo - if (oz_ind > 0) then - call read_vardata(dset, 'o3mr', ug3d,errcode=iret) + call read_vardata(dset, 'vgrd', vg3d,errcode=iret) if (iret /= 0) then - print *,'error reading o3mr' - call stop2(26) + print *,'error reading vgrd' + call stop2(23) endif - if (cliptracers) where (ug3d < clip) ug3d = clip do k=1,nlevs ug = reshape(ug3d(:,:,nlevs-k+1),(/nlons*nlats/)) - call copytogrdin(ug,grdin(:,levels(oz_ind-1)+k,nb,ne)) + vg = reshape(vg3d(:,:,nlevs-k+1),(/nlons*nlats/)) + if (u_ind > 0) call copytogrdin(ug,grdin(:,levels(u_ind-1) + k,nb,ne)) + if (v_ind > 0) call copytogrdin(vg,grdin(:,levels(v_ind-1) + k,nb,ne)) + ! calculate vertical integral of mass flux div (ps tendency) + ! this variable is analyzed in order to enforce mass balance in the analysis + if (pst_ind > 0) then + ug = ug*(pressi(:,k)-pressi(:,k+1)) + vg = vg*(pressi(:,k)-pressi(:,k+1)) + call sptezv_s(divspec,vrtspec,ug,vg,-1) ! u,v to div,vrt + call sptez_s(divspec,vmassdiv(:,k),1) ! divspec to divgrd + endif enddo - endif - if (cw_ind > 0 .or. ql_ind > 0 .or. qi_ind > 0) then - call read_vardata(dset, 'clwmr', ug3d,errcode=iret) + call read_vardata(dset,'tmp', ug3d,errcode=iret) if (iret /= 0) then - print *,'error reading clwmr' - call stop2(27) + print *,'error reading tmp' + call stop2(24) endif - if (imp_physics == 11) then - call read_vardata(dset, 'icmr', vg3d,errcode=iret) - if (iret /= 0) then - print *,'error reading icmr' - call stop2(28) - endif - ug3d = ug3d + vg3d + call read_vardata(dset,'spfh', vg3d,errcode=iret) + if (iret /= 0) then + print *,'error reading spfh' + call stop2(25) endif - if (cliptracers) where (ug3d < clip) ug3d = clip do k=1,nlevs ug = reshape(ug3d(:,:,nlevs-k+1),(/nlons*nlats/)) - call copytogrdin(ug,cw(:,k)) - if (cw_ind > 0) grdin(:,levels(cw_ind-1)+k,nb,ne) = cw(:,k) + vg = reshape(vg3d(:,:,nlevs-k+1),(/nlons*nlats/)) + if (tsen_ind > 0) call copytogrdin(ug,grdin(:,levels(tsen_ind-1)+k,nb,ne)) + call copytogrdin(vg, q(:,k)) + ug = ug * ( 1.0 + fv*vg ) ! convert T to Tv + call copytogrdin(ug,tv(:,k)) + if (tv_ind > 0) grdin(:,levels(tv_ind-1)+k,nb,ne) = tv(:,k) + if (q_ind > 0) grdin(:,levels( q_ind-1)+k,nb,ne) = q(:,k) enddo - endif - deallocate(ug3d,vg3d) - else -!$omp parallel do private(k,ug,vg,divspec,vrtspec) shared(sigdata,pressi,vmassdiv,grdin,tv,q,cw,u_ind,v_ind,pst_ind,q_ind,tsen_ind,cw_ind,qi_ind,ql_ind) - do k=1,nlevs - - vrtspec = sigdata%z(:,k); divspec = sigdata%d(:,k) - call sptezv_s(divspec,vrtspec,ug,vg,1) - if (u_ind > 0) then - call copytogrdin(ug,grdin(:,levels(u_ind-1)+k,nb,ne)) + if (oz_ind > 0) then + call read_vardata(dset, 'o3mr', ug3d,errcode=iret) + if (iret /= 0) then + print *,'error reading o3mr' + call stop2(26) + endif + if (cliptracers) where (ug3d < clip) ug3d = clip + do k=1,nlevs + ug = reshape(ug3d(:,:,nlevs-k+1),(/nlons*nlats/)) + call copytogrdin(ug,grdin(:,levels(oz_ind-1)+k,nb,ne)) + enddo endif - if (v_ind > 0) then - call copytogrdin(vg,grdin(:,levels(v_ind-1)+k,nb,ne)) + if (cw_ind > 0 .or. ql_ind > 0 .or. qi_ind > 0) then + call read_vardata(dset, 'clwmr', ug3d,errcode=iret) + if (iret /= 0) then + print *,'error reading clwmr' + call stop2(27) + endif + if (imp_physics == 11) then + call read_vardata(dset, 'icmr', vg3d,errcode=iret) + if (iret /= 0) then + print *,'error reading icmr' + call stop2(28) + endif + ug3d = ug3d + vg3d + endif + if (cliptracers) where (ug3d < clip) ug3d = clip + do k=1,nlevs + ug = reshape(ug3d(:,:,nlevs-k+1),(/nlons*nlats/)) + call copytogrdin(ug,cw(:,k)) + if (cw_ind > 0) grdin(:,levels(cw_ind-1)+k,nb,ne) = cw(:,k) + enddo endif + deallocate(ug3d,vg3d) + + else + !$omp parallel do private(k,ug,vg,divspec,vrtspec) shared(sigdata,pressi,vmassdiv,grdin,tv,q,cw,u_ind,v_ind,pst_ind,q_ind,tsen_ind,cw_ind,qi_ind,ql_ind) + do k=1,nlevs -! calculate vertical integral of mass flux div (ps tendency) -! this variable is analyzed in order to enforce mass balance in the analysis - if (pst_ind > 0) then - ug = ug*(pressi(:,k)-pressi(:,k+1)) - vg = vg*(pressi(:,k)-pressi(:,k+1)) - call sptezv_s(divspec,vrtspec,ug,vg,-1) ! u,v to div,vrt - call sptez_s(divspec,vmassdiv(:,k),1) ! divspec to divgrd - endif + vrtspec = sigdata%z(:,k); divspec = sigdata%d(:,k) + call sptezv_s(divspec,vrtspec,ug,vg,1) + if (u_ind > 0) then + call copytogrdin(ug,grdin(:,levels(u_ind-1)+k,nb,ne)) + endif + if (v_ind > 0) then + call copytogrdin(vg,grdin(:,levels(v_ind-1)+k,nb,ne)) + endif - divspec = sigdata%t(:,k) - call sptez_s(divspec,ug,1) - call copytogrdin(ug,tv(:,k)) - if (tv_ind > 0) grdin(:,levels(tv_ind-1)+k,nb,ne) = tv(:,k) + ! calculate vertical integral of mass flux div (ps tendency) + ! this variable is analyzed in order to enforce mass balance in the analysis + if (pst_ind > 0) then + ug = ug*(pressi(:,k)-pressi(:,k+1)) + vg = vg*(pressi(:,k)-pressi(:,k+1)) + call sptezv_s(divspec,vrtspec,ug,vg,-1) ! u,v to div,vrt + call sptez_s(divspec,vmassdiv(:,k),1) ! divspec to divgrd + endif - divspec = sigdata%q(:,k,1) - call sptez_s(divspec,vg,1) - call copytogrdin(vg,q(:,k)) - if (q_ind > 0) grdin(:,levels( q_ind-1)+k,nb,ne) = q(:,k) + divspec = sigdata%t(:,k) + call sptez_s(divspec,ug,1) + call copytogrdin(ug,tv(:,k)) + if (tv_ind > 0) grdin(:,levels(tv_ind-1)+k,nb,ne) = tv(:,k) - if (tsen_ind > 0) grdin(:,levels(tsen_ind-1)+k,nb,ne) = tv(:,k) / (one + fv*max(0._r_kind,q(:,k))) + divspec = sigdata%q(:,k,1) + call sptez_s(divspec,vg,1) + call copytogrdin(vg,q(:,k)) + if (q_ind > 0) grdin(:,levels( q_ind-1)+k,nb,ne) = q(:,k) - if (oz_ind > 0) then - divspec = sigdata%q(:,k,2) - call sptez_s(divspec,ug,1) - call copytogrdin(ug,grdin(:,levels(oz_ind-1)+k,nb,ne)) - endif + if (tsen_ind > 0) grdin(:,levels(tsen_ind-1)+k,nb,ne) = tv(:,k) / (one + fv*max(0._r_kind,q(:,k))) - if (cw_ind > 0 .or. ql_ind > 0 .or. qi_ind > 0) then - divspec = sigdata%q(:,k,3) - call sptez_s(divspec,ug,1) - call copytogrdin(ug,cw(:,k)) - if (cw_ind > 0) grdin(:,levels(cw_ind-1)+k,nb,ne) = cw(:,k) - endif + if (oz_ind > 0) then + divspec = sigdata%q(:,k,2) + call sptez_s(divspec,ug,1) + call copytogrdin(ug,grdin(:,levels(oz_ind-1)+k,nb,ne)) + endif - enddo -!$omp end parallel do - endif + if (cw_ind > 0 .or. ql_ind > 0 .or. qi_ind > 0) then + divspec = sigdata%q(:,k,3) + call sptez_s(divspec,ug,1) + call copytogrdin(ug,cw(:,k)) + if (cw_ind > 0) grdin(:,levels(cw_ind-1)+k,nb,ne) = cw(:,k) + endif - ! surface pressure - if (ps_ind > 0) then - call copytogrdin(psg,grdin(:,levels(n3d) + ps_ind,nb,ne)) - endif - if (.not. use_gfs_nemsio) call sigio_axdata(sigdata,iret) + enddo + !$omp end parallel do + endif - ! surface pressure tendency - if (pst_ind > 0) then - pstend = sum(vmassdiv,2) - if (nanal .eq. 1) & - print *,nanal,'min/max first-guess ps tend',minval(pstend),maxval(pstend) - call copytogrdin(pstend,grdin(:,levels(n3d) + pst_ind,nb,ne)) - endif + ! surface pressure + if (ps_ind > 0) then + call copytogrdin(psg,grdin(:,levels(n3d) + ps_ind,nb,ne)) + endif + if (.not. use_gfs_nemsio) call sigio_axdata(sigdata,iret) - ! compute saturation q. - do k=1,nlevs - ! pressure at bottom of layer interface (for gps jacobian, see prsltmp in setupbend.f90) - if (prse_ind > 0) then - ug(:) = pressi(:,k) + ! surface pressure tendency + if (pst_ind > 0) then + pstend = sum(vmassdiv,2) + if (nanal .eq. 1) & + print *,nanal,'min/max first-guess ps tend',minval(pstend),maxval(pstend) + call copytogrdin(pstend,grdin(:,levels(n3d) + pst_ind,nb,ne)) + endif + + ! compute saturation q. + do k=1,nlevs + ! pressure at bottom of layer interface (for gps jacobian, see prsltmp in setupbend.f90) + if (prse_ind > 0) then + ug(:) = pressi(:,k) + call copytogrdin(ug,pslg(:,k)) + ! Jacobian for gps in pressure is saved in different units in GSI; need to + ! multiply pressure by 0.1 + grdin(:,levels(prse_ind-1)+k,nb,ne) = 0.1*pslg(:,k) + endif + ! layer pressure from phillips vertical interolation (used for qsat + ! calculation) + ug(:) = ((pressi(:,k)**kap1-pressi(:,k+1)**kap1)/& + (kap1*(pressi(:,k)-pressi(:,k+1))))**kapr call copytogrdin(ug,pslg(:,k)) - ! Jacobian for gps in pressure is saved in different units in GSI; need to - ! multiply pressure by 0.1 - grdin(:,levels(prse_ind-1)+k,nb,ne) = 0.1*pslg(:,k) - endif - ! layer pressure from phillips vertical interolation (used for qsat - ! calculation) - ug(:) = ((pressi(:,k)**kap1-pressi(:,k+1)**kap1)/& - (kap1*(pressi(:,k)-pressi(:,k+1))))**kapr - call copytogrdin(ug,pslg(:,k)) - end do - if (pseudo_rh) then - call genqsat1(q,qsat(:,:,nb,ne),pslg,tv,ice,npts,nlevs) - else - qsat(:,:,nb,ne) = 1._r_double - end if + end do + if (pseudo_rh) then + call genqsat1(q,qsat(:,:,nb,ne),pslg,tv,ice,npts,nlevs) + else + qsat(:,:,nb,ne) = 1._r_double + end if - ! cloud derivatives - if (.not. use_full_hydro) then - if (ql_ind > 0 .or. qi_ind > 0) then - do k = 1, nlevs - do i = 1, npts - qi_coef = -r0_05*(tv(i,k)/(one+fv*q(i,k))-t0c) - qi_coef = max(zero,qi_coef) - qi_coef = min(one,qi_coef) ! 0<=qi_coef<=1 - if (ql_ind > 0) then - grdin(i,levels(ql_ind-1)+k,nb,ne) = cw(i,k)*(one-qi_coef) - endif - if (qi_ind > 0) then - grdin(i,levels(qi_ind-1)+k,nb,ne) = cw(i,k)*qi_coef - endif + ! cloud derivatives + if (.not. use_full_hydro) then + if (ql_ind > 0 .or. qi_ind > 0) then + do k = 1, nlevs + do i = 1, npts + qi_coef = -r0_05*(tv(i,k)/(one+fv*q(i,k))-t0c) + qi_coef = max(zero,qi_coef) + qi_coef = min(one,qi_coef) ! 0<=qi_coef<=1 + if (ql_ind > 0) then + grdin(i,levels(ql_ind-1)+k,nb,ne) = cw(i,k)*(one-qi_coef) + endif + if (qi_ind > 0) then + grdin(i,levels(qi_ind-1)+k,nb,ne) = cw(i,k)*qi_coef + endif + enddo enddo - enddo - endif - endif + endif + endif - if (sst_ind > 0) then - grdin(:,levels(n3d)+sst_ind, nb,ne) = zero - endif + if (sst_ind > 0) then + grdin(:,levels(n3d)+sst_ind, nb,ne) = zero + endif + + deallocate(pressi,pslg) + deallocate(psg) + if (pst_ind > 0) deallocate(vmassdiv,pstend) + endif ! read_atm_file - deallocate(pressi,pslg) - deallocate(psg) - if (pst_ind > 0) deallocate(vmassdiv,pstend) if (use_gfs_nemsio) call nemsio_close(gfile,iret=iret) if (use_gfs_ncio) call close_dataset(dset) if (use_gfs_nemsio) call nemsio_close(gfilesfc,iret=iret) + if ( read_sfc_file ) then + + if ( .not. use_gfs_ncio ) then + write(6,*) 'griddio/griddata for sfc update vars only coded for nc io' + call stop2(23) + endif + if ( reducedgrid ) then + write(6,*) "reducedgrid=T not valid for global_2mDA=T" + call stop2(22) + endif + + dset_sfc = open_dataset(filenamesfc) + ! 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 *,'error reading tmp2m' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + call copytogrdin(ug,grdin(:,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 *,'error reading spfh2m' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + call copytogrdin(ug,grdin(:,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 *,'error reading soilt1' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + call copytogrdin(ug,grdin(:,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 *,'error reading soilt2' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + call copytogrdin(ug,grdin(:,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 *,'error reading soilt3' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + call copytogrdin(ug,grdin(:,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 *,'error reading soilt2' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + call copytogrdin(ug,grdin(:,soilt4_ind,nb,ne)) + endif + if (slc1_ind > 0) then + call read_vardata(dset_sfc, 'slc1', values_2d, errcode=iret) + if (iret /= 0) then + print *,'error reading slc1' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + call copytogrdin(ug,grdin(:,slc1_ind,nb,ne)) + endif + if (slc2_ind > 0) then + call read_vardata(dset_sfc, 'slc2', values_2d, errcode=iret) + if (iret /= 0) then + print *,'error reading slc2' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + call copytogrdin(ug,grdin(:,slc2_ind,nb,ne)) + endif + if (slc3_ind > 0) then + call read_vardata(dset_sfc, 'slc3', values_2d, errcode=iret) + if (iret /= 0) then + print *,'error reading slc3' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + call copytogrdin(ug,grdin(:,slc3_ind,nb,ne)) + endif + if (slc4_ind > 0) then + call read_vardata(dset_sfc, 'slc4', values_2d, errcode=iret) + if (iret /= 0) then + print *,'error reading slc2' + call stop2(22) + endif + ug = reshape(values_2d,(/nlons*nlats/)) + call copytogrdin(ug,grdin(:,slc4_ind,nb,ne)) + endif + + call close_dataset(dset_sfc) + + endif ! sfc read + + if ( allocated(values_2d) ) deallocate(values_2d) + end do backgroundloop ! loop over backgrounds to read in end do ensmemloop ! loop over ens members to read in diff --git a/src/enkf/mpi_readobs.f90 b/src/enkf/mpi_readobs.f90 index e332609d30..db9c276f7d 100644 --- a/src/enkf/mpi_readobs.f90 +++ b/src/enkf/mpi_readobs.f90 @@ -89,7 +89,7 @@ subroutine mpi_getobs(obspath, datestring, nobs_conv, nobs_oz, nobs_sat, nobs_to ! get total number of conventional and sat obs for ensmean. id = 'ensmean' if(nproc == 0)call get_num_convobs(obspath,datestring,nobs_conv,nobs_convdiag,id) - if (global_2mDA) then + if (global_2mDA) then ! flag for not reading other conventional obs if(nproc == iozproc) nobs_oz=0; nobs_ozdiag=0 if(nproc == isatproc) nobs_sat=0; nobs_satdiag=0 else diff --git a/src/enkf/params.f90 b/src/enkf/params.f90 index 1f804e89b3..6416ceca35 100644 --- a/src/enkf/params.f90 +++ b/src/enkf/params.f90 @@ -484,9 +484,6 @@ subroutine read_namelist() endif close(912) -! all fields are 2d if global_2mDA -if (global_2mDA) nlevs=1 - ! find number of satellite files nsats_rad=0 do i=1,nsatmax_rad @@ -717,11 +714,7 @@ subroutine read_namelist() fgfileprefixes(nbackgrounds+1)="firstguess." endif else ! global - if (global_2mDA) then - fgfileprefixes(nbackgrounds+1)="bfg_"//datestring//"_fhr"//charfhr_anal(nbackgrounds+1)//"_" - else fgfileprefixes(nbackgrounds+1)="sfg_"//datestring//"_fhr"//charfhr_anal(nbackgrounds+1)//"_" - endif endif endif if (trim(fgsfcfileprefixes(nbackgrounds+1)) .eq. "") then @@ -743,11 +736,11 @@ subroutine read_namelist() statefileprefixes(nstatefields+1)="firstguess." endif else ! global - if (global_2mda) then - statefileprefixes(nstatefields+1)="bfg_"//datestring//"_fhr"//charfhr_state(nstatefields+1)//"_" - else + !if (global_2mda) then + !statefileprefixes(nstatefields+1)="bfg_"//datestring//"_fhr"//charfhr_state(nstatefields+1)//"_" + !else statefileprefixes(nstatefields+1)="sfg_"//datestring//"_fhr"//charfhr_state(nstatefields+1)//"_" - endif + !endif endif endif if (trim(statesfcfileprefixes(nstatefields+1)) .eq. "") then diff --git a/src/enkf/readconvobs.f90 b/src/enkf/readconvobs.f90 index 2db8b1c77e..6b47d58ce4 100644 --- a/src/enkf/readconvobs.f90 +++ b/src/enkf/readconvobs.f90 @@ -282,7 +282,7 @@ subroutine get_num_convobs_nc(obspath,datestring,num_obs_tot,num_obs_totdiag,id) obtype = obtypes(itype) ! only read t and q obs for global_2mDA - if (global_2mDA .and. (obtype .ne. ' t' .and. obtype .ne. ' q')) cycle obtypeloop + if (global_2mDA .and. (obtype .ne. ' t' .and. obtype .ne. ' q')) cycle obtypeloop ! flag for nor reading other convention obs type peloop: do ipe=0,npefiles write(pe_name,'(i4.4)') ipe @@ -338,10 +338,10 @@ subroutine get_num_convobs_nc(obspath,datestring,num_obs_tot,num_obs_totdiag,id) num_obs_totdiag = num_obs_totdiag + nobs_curr do i = 1, nobs_curr - ! for global_2mDA skip if not 2m (surface) ob + ! for global_2mDA skip if not 2m (surface) ob ityp = Observation_Type(i) sfctype=(ityp>179.and.ityp<190).or.(ityp>=192.and.ityp<=199) - if (global_2mDA .and. .not. sfctype) cycle + if (global_2mDA .and. .not. sfctype) cycle ! flag for not reading other conventional obs types errorlimit2=errorlimit2_obs @@ -350,7 +350,7 @@ subroutine get_num_convobs_nc(obspath,datestring,num_obs_tot,num_obs_totdiag,id) endif ! for q, normalize by qsatges - if (obtype == ' q' .and. .not. global_2mDA) then + if (obtype == ' q' .and. .not. global_2mDA) then ! flag for normalzing q by ens mean (if q read in) obmax = abs(Observation(i) / Forecast_Saturation_Spec_Hum(i)) error = Errinv_Final(i) * Forecast_Saturation_Spec_Hum(i) else @@ -553,7 +553,7 @@ subroutine get_convobs_data_nc(obspath, datestring, nobs_max, nobs_maxdiag, & obtype = obtypes(itype) ! only read t and q obs for global_2mDA - if (global_2mDA .and. (obtype .ne. ' t' .and. obtype .ne. ' q')) cycle obtypeloop + if (global_2mDA .and. (obtype .ne. ' t' .and. obtype .ne. ' q')) cycle obtypeloop ! flag for not reading other obs types peloop: do ipe=0,npefiles write(pe_name,'(i4.4)') ipe @@ -678,7 +678,7 @@ subroutine get_convobs_data_nc(obspath, datestring, nobs_max, nobs_maxdiag, & ! for global_2mDA skip if not 2m (surface) ob ityp = Observation_Type(i) sfctype=(ityp>179.and.ityp<190).or.(ityp>=192.and.ityp<=199) - if (global_2mDA .and. .not. sfctype) cycle + if (global_2mDA .and. .not. sfctype) cycle ! flag for not reading other conventional obs types nobdiag = nobdiag + 1 ! special handling for error limits for GPS bend angle if (obtype == 'gps') then @@ -686,7 +686,7 @@ subroutine get_convobs_data_nc(obspath, datestring, nobs_max, nobs_maxdiag, & endif ! for q, normalize by qsatges - if (.not. global_2mDA .and. obtype == ' q') then + if (.not. global_2mDA .and. obtype == ' q') then ! flag for normazling q obmax = abs(real(Observation(i),r_single) / real(Forecast_Saturation_Spec_Hum(i),r_single)) errororig = real(Errinv_Input(i),r_single) * real(Forecast_Saturation_Spec_Hum(i),r_single) error = real(Errinv_Final(i),r_single) * real(Forecast_Saturation_Spec_Hum(i),r_single) @@ -795,13 +795,13 @@ subroutine get_convobs_data_nc(obspath, datestring, nobs_max, nobs_maxdiag, & endif ! normalize q by qsatges - if (obtype == ' q' .and. .not. global_2mDA) then + if (obtype == ' q' .and. .not. global_2mDA) then ! flag for normalizing q hx(nob) = hx(nob) / Forecast_Saturation_Spec_Hum(i) endif endif ! normalize q by qsatges - if (obtype == ' q' .and. .not. global_2mDA) then + if (obtype == ' q' .and. .not. global_2mDA) then ! flag for normalizing q x_obs(nob) = x_obs(nob) /Forecast_Saturation_Spec_Hum(i) hx_mean(nob) = hx_mean(nob) /Forecast_Saturation_Spec_Hum(i) hx_mean_nobc(nob) = hx_mean_nobc(nob) /Forecast_Saturation_Spec_Hum(i) diff --git a/src/enkf/statevec.f90 b/src/enkf/statevec.f90 index d2641327cc..4469550e86 100644 --- a/src/enkf/statevec.f90 +++ b/src/enkf/statevec.f90 @@ -14,7 +14,7 @@ module statevec ! ! Public Variables: ! nanals: (integer scalar) number of ensemble members (from module params) -! nlevs: number of analysis vertical levels (from module params). +! nlevs: number of analysis atmos vertical levels (from module params). ! ns3d: number of 3D variables ! ns2d: number of 2D variables ! svars3d: names of 3D variables @@ -39,7 +39,7 @@ module statevec ! !$$$ -use gridio, only: readgriddata, readgriddata_pnc, readgriddata_2mDA +use gridio, only: readgriddata, readgriddata_pnc use mpisetup, only: mpi_real4,mpi_sum,mpi_comm_io,mpi_in_place,numproc,nproc use mpimod, only: mpi_comm_world use gridinfo, only: getgridinfo, gridinfo_cleanup, & @@ -120,7 +120,7 @@ subroutine init_statevec() svars3d(ns3d)=trim(adjustl(var)) slevels(ns3d)=ilev + slevels(ns3d-1) else - if (nproc .eq. 0) print *,'Error: only ', nlevs, ' and ', nlevs+1,' number of levels is supported in current version, got ',ilev + if (nproc .eq. 0) print *,'Error: statevec - only ', nlevs, ' and ', nlevs+1,' number of levels is supported in current version, got ',ilev call stop2(503) endif enddo @@ -188,7 +188,7 @@ subroutine read_state() allocate(state_d(npts,nsdim,nstatefields,nanals_per_iotask)) allocate(qsat(npts,nlevs,nstatefields,nanals_per_iotask)) if (paranc) then - if (global_2mDA) then + if (global_2mDA) then ! flag for parallel rad print *,'paranc not supported for global_2mDA' call mpi_barrier(mpi_comm_world,ierr) call mpi_finalize(ierr) @@ -200,13 +200,8 @@ subroutine read_state() if (nproc <= ntasks_io-1) then nanal = nproc + 1 if ( .not. paranc) then - if (global_2mDA) then - call readgriddata_2mDA(nanal1(nproc),nanal2(nproc),svars2d,ns2d,nsdim,nstatefields, & - statesfcfileprefixes,.false.,state_d) - else call readgriddata(nanal1(nproc),nanal2(nproc),svars3d,svars2d,ns3d,ns2d,slevels,nsdim,nstatefields, & statefileprefixes,statesfcfileprefixes,.false.,state_d,qsat) - endif end if ! subtract the mean diff --git a/src/gsi/setupt.f90 b/src/gsi/setupt.f90 index e8325d5b84..bee321ba2a 100644 --- a/src/gsi/setupt.f90 +++ b/src/gsi/setupt.f90 @@ -503,7 +503,9 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav if (lobsdiagsave) nreal=nreal+4*miter+1 if (twodvar_regional) then; nreal=nreal+2; allocate(cprvstg(nobs),csprvstg(nobs)); endif if (save_jacobian) then - nnz = 2 ! number of non-zero elements in dH(x)/dx profile + write(6,*) 'CSD - setting nnz to 1' + !nnz = 2 ! number of non-zero elements in dH(x)/dx profile + nnz = 1 ! for T2m. nind = 1 call new(dhx_dx, nnz, nind) nreal = nreal + size(dhx_dx) From e7e87547bc4291814b8494a3df106f2e34ec48b3 Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Tue, 26 Jul 2022 14:19:16 -0500 Subject: [PATCH 07/12] Fixed bug in sfc variable indexing into grdin --- src/enkf/controlvec.f90 | 8 ++--- src/enkf/gridio_gfs.f90 | 62 ++++++++++++++++++++------------------- src/enkf/observer_gfs.f90 | 1 + 3 files changed, 37 insertions(+), 34 deletions(-) diff --git a/src/enkf/controlvec.f90 b/src/enkf/controlvec.f90 index 4789bdff94..e57fba9395 100644 --- a/src/enkf/controlvec.f90 +++ b/src/enkf/controlvec.f90 @@ -371,13 +371,13 @@ subroutine write_control(no_inflate_flag) if (.not. paranc) then if (write_fv3_incr) then if (global_2mDA) then ! flag for writing sfc increments - call writeincrement_2mDA(nanal1(nproc),nanal2(nproc),cvars2d,nc2d,ncdim,grdin,no_inflate_flag) + call writeincrement_2mDA(nanal1(nproc),nanal2(nproc),cvars2d,nc2d,grdin(:,clevels(nc3d)+1:ncdim, :, :),no_inflate_flag) else call writeincrement(nanal1(nproc),nanal2(nproc),cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin,no_inflate_flag) endif else if (global_2mDA) then ! flag for witing sfc analyis - call writegriddata_2mDA(nanal1(nproc),nanal2(nproc),cvars2d,nc2d,ncdim,grdin,no_inflate_flag) + call writegriddata_2mDA(nanal1(nproc),nanal2(nproc),cvars2d,nc2d,grdin(:,clevels(nc3d)+1:ncdim,:,:),no_inflate_flag) else call writegriddata(nanal1(nproc),nanal2(nproc),cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin,no_inflate_flag) endif @@ -387,13 +387,13 @@ subroutine write_control(no_inflate_flag) ! also write out ens mean on root task. if (write_fv3_incr) then if (global_2mDA) then ! flag for writing sfc inctements - call writeincrement_2mDA(0,0,cvars2d,nc2d,ncdim,grdin_mean,no_inflate_flag) + call writeincrement_2mDA(0,0,cvars2d,nc2d,grdin_mean(:,clevels(nc3d)+1:ncdim,:,:),no_inflate_flag) else call writeincrement(0,0,cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin_mean,no_inflate_flag) endif else if (global_2mDA) then ! flag for writing sfc analysis - call writegriddata_2mDA(0,0,cvars2d,nc2d,ncdim,grdin_mean,no_inflate_flag) + call writegriddata_2mDA(0,0,cvars2d,nc2d,grdin_mean(:,clevels(nc3d)+1:ncdim,:,:),no_inflate_flag) else call writegriddata(0,0,cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin_mean,no_inflate_flag) endif diff --git a/src/enkf/gridio_gfs.f90 b/src/enkf/gridio_gfs.f90 index 1b882e272f..ac08671f45 100644 --- a/src/enkf/gridio_gfs.f90 +++ b/src/enkf/gridio_gfs.f90 @@ -1043,7 +1043,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, if (use_gfs_nemsio) call nemsio_close(gfilesfc,iret=iret) if ( read_sfc_file ) then - + if ( .not. use_gfs_ncio ) then write(6,*) 'griddio/griddata for sfc update vars only coded for nc io' call stop2(23) @@ -1062,7 +1062,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, call stop2(22) endif ug = reshape(values_2d,(/nlons*nlats/)) - call copytogrdin(ug,grdin(:,tmp2m_ind,nb,ne)) + 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) @@ -1071,7 +1071,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, call stop2(22) endif ug = reshape(values_2d,(/nlons*nlats/)) - call copytogrdin(ug,grdin(:,spfh2m_ind,nb,ne)) + 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) @@ -1080,7 +1080,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, call stop2(22) endif ug = reshape(values_2d,(/nlons*nlats/)) - call copytogrdin(ug,grdin(:,soilt1_ind,nb,ne)) + 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) @@ -1089,7 +1089,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, call stop2(22) endif ug = reshape(values_2d,(/nlons*nlats/)) - call copytogrdin(ug,grdin(:,soilt2_ind,nb,ne)) + 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) @@ -1098,7 +1098,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, call stop2(22) endif ug = reshape(values_2d,(/nlons*nlats/)) - call copytogrdin(ug,grdin(:,soilt3_ind,nb,ne)) + 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) @@ -1107,7 +1107,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, call stop2(22) endif ug = reshape(values_2d,(/nlons*nlats/)) - call copytogrdin(ug,grdin(:,soilt4_ind,nb,ne)) + call copytogrdin(ug,grdin(:,levels(n3d) + soilt4_ind,nb,ne)) endif if (slc1_ind > 0) then call read_vardata(dset_sfc, 'slc1', values_2d, errcode=iret) @@ -1116,7 +1116,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, call stop2(22) endif ug = reshape(values_2d,(/nlons*nlats/)) - call copytogrdin(ug,grdin(:,slc1_ind,nb,ne)) + call copytogrdin(ug,grdin(:,levels(n3d) + slc1_ind,nb,ne)) endif if (slc2_ind > 0) then call read_vardata(dset_sfc, 'slc2', values_2d, errcode=iret) @@ -1125,7 +1125,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, call stop2(22) endif ug = reshape(values_2d,(/nlons*nlats/)) - call copytogrdin(ug,grdin(:,slc2_ind,nb,ne)) + call copytogrdin(ug,grdin(:,levels(n3d) + slc2_ind,nb,ne)) endif if (slc3_ind > 0) then call read_vardata(dset_sfc, 'slc3', values_2d, errcode=iret) @@ -1134,7 +1134,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, call stop2(22) endif ug = reshape(values_2d,(/nlons*nlats/)) - call copytogrdin(ug,grdin(:,slc3_ind,nb,ne)) + call copytogrdin(ug,grdin(:,levels(n3d) + slc3_ind,nb,ne)) endif if (slc4_ind > 0) then call read_vardata(dset_sfc, 'slc4', values_2d, errcode=iret) @@ -1143,7 +1143,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, call stop2(22) endif ug = reshape(values_2d,(/nlons*nlats/)) - call copytogrdin(ug,grdin(:,slc4_ind,nb,ne)) + call copytogrdin(ug,grdin(:,levels(n3d) + slc4_ind,nb,ne)) endif call close_dataset(dset_sfc) @@ -1976,7 +1976,7 @@ end subroutine copyfromgrdin end subroutine writegriddata_pnc - subroutine writegriddata_2mDA(nanal1,nanal2,vars2d,n2d,ndim,grdin,no_inflate_flag) + subroutine writegriddata_2mDA(nanal1,nanal2,vars2d,n2d,grdin_2d,no_inflate_flag) use netcdf use module_fv3gfs_ncio, only: Dataset, Variable, Dimension, open_dataset,& read_attribute, close_dataset, get_dim, read_vardata,& @@ -1989,8 +1989,9 @@ subroutine writegriddata_2mDA(nanal1,nanal2,vars2d,n2d,ndim,grdin,no_inflate_fla integer, intent(in) :: nanal1,nanal2 character(len=max_varname_length), dimension(n2d), intent(in) :: vars2d - integer, intent(in) :: n2d,ndim - real(r_single), dimension(npts,ndim,nbackgrounds,nanal2-nanal1+1), intent(inout) :: grdin + integer, intent(in) :: n2d + ! only 2d segment of grdin should be input here + real(r_single), dimension(npts,n2d,nbackgrounds,nanal2-nanal1+1), intent(inout) :: grdin_2d logical, intent(in) :: no_inflate_flag character(len=500):: filenamein, filenameout real(r_kind), dimension(nlons*nlats) :: ug,vg @@ -2126,7 +2127,7 @@ subroutine writegriddata_2mDA(nanal1,nanal2,vars2d,n2d,ndim,grdin,no_inflate_fla call stop2(29) endif vg = reshape(values_2d,(/nlons*nlats/)) - call copyfromgrdin(grdin(:,tmp2m_ind,nb,ne),ug) + call copyfromgrdin(grdin_2d(:,tmp2m_ind,nb,ne),ug) ! add increment to background. values_2d = reshape(ug+vg,(/nlons,nlats/)) ! quantize data if lossy compression desired. @@ -2155,7 +2156,7 @@ subroutine writegriddata_2mDA(nanal1,nanal2,vars2d,n2d,ndim,grdin,no_inflate_fla call stop2(29) endif vg = reshape(values_2d,(/nlons*nlats/)) - call copyfromgrdin(grdin(:,spfh2m_ind,nb,ne),ug) + call copyfromgrdin(grdin_2d(:,spfh2m_ind,nb,ne),ug) ! add increment to background. values_2d = reshape(ug+vg,(/nlons,nlats/)) ! quantize data if lossy compression desired. @@ -2184,7 +2185,7 @@ subroutine writegriddata_2mDA(nanal1,nanal2,vars2d,n2d,ndim,grdin,no_inflate_fla call stop2(29) endif vg = reshape(values_2d,(/nlons*nlats/)) - call copyfromgrdin(grdin(:,soilt1_ind,nb,ne),ug) + call copyfromgrdin(grdin_2d(:,soilt1_ind,nb,ne),ug) ! add increment to background. values_2d = reshape(ug+vg,(/nlons,nlats/)) ! quantize data if lossy compression desired. @@ -3708,7 +3709,7 @@ end subroutine copyfromgrdin end subroutine writegriddata - subroutine writeincrement_2mDA(nanal1,nanal2,vars2d,n2d,ndim,grdin,no_inflate_flag) + subroutine writeincrement_2mDA(nanal1,nanal2,vars2d,n2d,grdin_2d,no_inflate_flag) use netcdf use params, only: nbackgrounds,fgsfcfileprefixes,incsfcfileprefixes,& datestring,nhr_anal,write_ensmean,reducedgrid @@ -3722,8 +3723,9 @@ subroutine writeincrement_2mDA(nanal1,nanal2,vars2d,n2d,ndim,grdin,no_inflate_fl integer, intent(in) :: nanal1,nanal2 character(len=max_varname_length), dimension(n2d), intent(in) :: vars2d - integer, intent(in) :: n2d,ndim - real(r_single), dimension(npts,ndim,nbackgrounds,1), intent(inout) :: grdin + integer, intent(in) :: n2d + ! only the 2d segment of grdin is input + real(r_single), dimension(npts,n2d,nbackgrounds,1), intent(inout) :: grdin_2d logical, intent(in) :: no_inflate_flag character(len=500):: filenamein, filenameout integer(i_kind) :: i, j, nb, ne, nanal @@ -3874,7 +3876,7 @@ subroutine writeincrement_2mDA(nanal1,nanal2,vars2d,n2d,ndim,grdin,no_inflate_fl ! tmp2m increment inc(:) = zero if (tmp2m_ind > 0) then - call copyfromgrdin(grdin(:,tmp2m_ind,nb,ne),inc) + call copyfromgrdin(grdin_2d(:,tmp2m_ind,nb,ne),inc) endif inc2d(:,:) = reshape(inc,(/nlons,nlats/)) do j=1,nlats @@ -3885,7 +3887,7 @@ subroutine writeincrement_2mDA(nanal1,nanal2,vars2d,n2d,ndim,grdin,no_inflate_fl ! spfh2m increment inc(:) = zero if (spfh2m_ind > 0) then - call copyfromgrdin(grdin(:,spfh2m_ind,nb,ne),inc) + call copyfromgrdin(grdin_2d(:,spfh2m_ind,nb,ne),inc) endif inc2d(:,:) = reshape(inc,(/nlons,nlats/)) do j=1,nlats @@ -3896,7 +3898,7 @@ subroutine writeincrement_2mDA(nanal1,nanal2,vars2d,n2d,ndim,grdin,no_inflate_fl ! soilt1 increment inc(:) = zero if (soilt1_ind > 0) then - call copyfromgrdin(grdin(:,soilt1_ind,nb,ne),inc) + call copyfromgrdin(grdin_2d(:,soilt1_ind,nb,ne),inc) endif inc2d(:,:) = reshape(inc,(/nlons,nlats/)) inc2dout=0. @@ -3910,7 +3912,7 @@ subroutine writeincrement_2mDA(nanal1,nanal2,vars2d,n2d,ndim,grdin,no_inflate_fl ! soilt2 increment inc(:) = zero if (soilt2_ind > 0) then - call copyfromgrdin(grdin(:,soilt2_ind,nb,ne),inc) + call copyfromgrdin(grdin_2d(:,soilt2_ind,nb,ne),inc) endif inc2d(:,:) = reshape(inc,(/nlons,nlats/)) inc2dout=0. @@ -3924,7 +3926,7 @@ subroutine writeincrement_2mDA(nanal1,nanal2,vars2d,n2d,ndim,grdin,no_inflate_fl ! soilt3 increment inc(:) = zero if (soilt3_ind > 0) then - call copyfromgrdin(grdin(:,soilt3_ind,nb,ne),inc) + call copyfromgrdin(grdin_2d(:,soilt3_ind,nb,ne),inc) endif inc2d(:,:) = reshape(inc,(/nlons,nlats/)) inc2dout=0. @@ -3938,7 +3940,7 @@ subroutine writeincrement_2mDA(nanal1,nanal2,vars2d,n2d,ndim,grdin,no_inflate_fl ! soilt4 increment inc(:) = zero if (soilt4_ind > 0) then - call copyfromgrdin(grdin(:,soilt4_ind,nb,ne),inc) + call copyfromgrdin(grdin_2d(:,soilt4_ind,nb,ne),inc) endif inc2d(:,:) = reshape(inc,(/nlons,nlats/)) inc2dout=0. @@ -3952,7 +3954,7 @@ subroutine writeincrement_2mDA(nanal1,nanal2,vars2d,n2d,ndim,grdin,no_inflate_fl ! slc1 increment inc(:) = zero if (slc1_ind > 0) then - call copyfromgrdin(grdin(:,slc1_ind,nb,ne),inc) + call copyfromgrdin(grdin_2d(:,slc1_ind,nb,ne),inc) endif inc2d(:,:) = reshape(inc,(/nlons,nlats/)) do j=1,nlats @@ -3963,7 +3965,7 @@ subroutine writeincrement_2mDA(nanal1,nanal2,vars2d,n2d,ndim,grdin,no_inflate_fl ! slc2 increment inc(:) = zero if (slc2_ind > 0) then - call copyfromgrdin(grdin(:,slc2_ind,nb,ne),inc) + call copyfromgrdin(grdin_2d(:,slc2_ind,nb,ne),inc) endif inc2d(:,:) = reshape(inc,(/nlons,nlats/)) do j=1,nlats @@ -3974,7 +3976,7 @@ subroutine writeincrement_2mDA(nanal1,nanal2,vars2d,n2d,ndim,grdin,no_inflate_fl ! slc3 increment inc(:) = zero if (slc3_ind > 0) then - call copyfromgrdin(grdin(:,slc3_ind,nb,ne),inc) + call copyfromgrdin(grdin_2d(:,slc3_ind,nb,ne),inc) endif inc2d(:,:) = reshape(inc,(/nlons,nlats/)) do j=1,nlats @@ -3985,7 +3987,7 @@ subroutine writeincrement_2mDA(nanal1,nanal2,vars2d,n2d,ndim,grdin,no_inflate_fl ! slc4 increment inc(:) = zero if (slc4_ind > 0) then - call copyfromgrdin(grdin(:,slc4_ind,nb,ne),inc) + call copyfromgrdin(grdin_2d(:,slc4_ind,nb,ne),inc) endif inc2d(:,:) = reshape(inc,(/nlons,nlats/)) do j=1,nlats diff --git a/src/enkf/observer_gfs.f90 b/src/enkf/observer_gfs.f90 index 983b25f959..6104e20b9a 100644 --- a/src/enkf/observer_gfs.f90 +++ b/src/enkf/observer_gfs.f90 @@ -178,6 +178,7 @@ subroutine calc_linhx(hx, dens, dhx_dx, hxpert, hx_ens, & ! interpolate state horizontally and in time and do dot product with dHx/dx profile ! saves from calculating interpolated x_ens for each state variable hx_ens = hx + do i = 1, dhx_dx%nnz j = dhx_dx%ind(i) hxpert%val(i) = (( dens( ix*nlons + iy , j, it) *delxp*delyp & From 53c3ebb52c84c0d334e85963d99e16629ca4a730 Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Tue, 26 Jul 2022 15:52:45 -0500 Subject: [PATCH 08/12] changed dhx_dx dimension in setupt to cover both profile and single layer obs. --- src/gsi/setupt.f90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/gsi/setupt.f90 b/src/gsi/setupt.f90 index bee321ba2a..a5e511a4f6 100644 --- a/src/gsi/setupt.f90 +++ b/src/gsi/setupt.f90 @@ -503,9 +503,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav if (lobsdiagsave) nreal=nreal+4*miter+1 if (twodvar_regional) then; nreal=nreal+2; allocate(cprvstg(nobs),csprvstg(nobs)); endif if (save_jacobian) then - write(6,*) 'CSD - setting nnz to 1' - !nnz = 2 ! number of non-zero elements in dH(x)/dx profile - nnz = 1 ! for T2m. + nnz = 2 ! number of non-zero elements in dH(x)/dx profile nind = 1 call new(dhx_dx, nnz, nind) nreal = nreal + size(dhx_dx) @@ -725,7 +723,6 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav elseif (sfctype .and. use_2m_obs ) then ! SCENARIO 2: obs is sfctype, and use_2m_obs scheme is on. ! 2m forecast has been read from the sfc guess files -! note: any changes to ges_t2m and tob in this block should also be made in buddy_check_t ! mask: 0 - sea, 1 - land, 2-ice (val + 3 = interpolated from at least one different land cover ! for now, use only pure land (interpolation from mixed grid cells is pretty sketchy, but @@ -754,6 +751,11 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav dhx_dx%st_ind(1) = sum(levels(1:ns3d)) + t_ind dhx_dx%end_ind(1) = sum(levels(1:ns3d)) + t_ind dhx_dx%val(1) = one + dhx_dx%val(2) = zero ! in this case, there is no vertical interp + ! and nnz (=dim(dhx_dx%val)) should be one, + ! but nnz is a file attribute, so need to use + ! same value as for vertical profile obs. Get + ! around this by setting val(2) to zero. endif else ! SCENARIO 3: obs is sfctype, and neither sfcmodel nor use_2m_obs is chosen From 965ca0fa741d96f3acf3f77dd3ee3c54927d600e Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Fri, 12 Aug 2022 15:12:50 -0500 Subject: [PATCH 09/12] Combined write routines to single routine for sfc and/or atmos fields. --- src/enkf/controlvec.f90 | 15 +- src/enkf/gridio_gfs.f90 | 1298 ++++++++++++++++++++------------------- 2 files changed, 655 insertions(+), 658 deletions(-) diff --git a/src/enkf/controlvec.f90 b/src/enkf/controlvec.f90 index e57fba9395..3eebe2ef4a 100644 --- a/src/enkf/controlvec.f90 +++ b/src/enkf/controlvec.f90 @@ -48,7 +48,7 @@ module controlvec use gridio, only: readgriddata, readgriddata_pnc, writegriddata, writegriddata_pnc, & writeincrement, writeincrement_pnc, & - writegriddata_2mDA, writeincrement_2mDA + writegriddata_2mDA use gridinfo, only: getgridinfo, gridinfo_cleanup, & npts, vars3d_supported, vars2d_supported use params, only: nlevs, nbackgrounds, fgfileprefixes, reducedgrid, & @@ -370,27 +370,15 @@ subroutine write_control(no_inflate_flag) if (.not. paranc) then if (write_fv3_incr) then - if (global_2mDA) then ! flag for writing sfc increments - call writeincrement_2mDA(nanal1(nproc),nanal2(nproc),cvars2d,nc2d,grdin(:,clevels(nc3d)+1:ncdim, :, :),no_inflate_flag) - else call writeincrement(nanal1(nproc),nanal2(nproc),cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin,no_inflate_flag) - endif else - if (global_2mDA) then ! flag for witing sfc analyis - call writegriddata_2mDA(nanal1(nproc),nanal2(nproc),cvars2d,nc2d,grdin(:,clevels(nc3d)+1:ncdim,:,:),no_inflate_flag) - else call writegriddata(nanal1(nproc),nanal2(nproc),cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin,no_inflate_flag) - endif end if if (nproc == 0) then if (write_ensmean) then ! also write out ens mean on root task. if (write_fv3_incr) then - if (global_2mDA) then ! flag for writing sfc inctements - call writeincrement_2mDA(0,0,cvars2d,nc2d,grdin_mean(:,clevels(nc3d)+1:ncdim,:,:),no_inflate_flag) - else call writeincrement(0,0,cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin_mean,no_inflate_flag) - endif else if (global_2mDA) then ! flag for writing sfc analysis call writegriddata_2mDA(0,0,cvars2d,nc2d,grdin_mean(:,clevels(nc3d)+1:ncdim,:,:),no_inflate_flag) @@ -448,3 +436,4 @@ subroutine controlvec_cleanup() end subroutine controlvec_cleanup end module controlvec + \ No newline at end of file diff --git a/src/enkf/gridio_gfs.f90 b/src/enkf/gridio_gfs.f90 index ac08671f45..55873a451d 100644 --- a/src/enkf/gridio_gfs.f90 +++ b/src/enkf/gridio_gfs.f90 @@ -57,7 +57,7 @@ module gridio private public :: readgriddata, readgriddata_pnc, writegriddata, writegriddata_pnc public :: writeincrement, writeincrement_pnc, & - writegriddata_2mDA, writeincrement_2mDA + writegriddata_2mDA contains @@ -477,6 +477,9 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, use_full_hydro = .false. + ! determine which files will be read in + call set_ncio_file_flags(vars3d, n3d, vars2d, n2d, read_sfc_file, read_atm_file) + u_ind = getindex(vars3d, 'u') !< indices in the state or control var arrays v_ind = getindex(vars3d, 'v') ! U and V (3D) tv_ind = getindex(vars3d, 'tv') ! Tv (3D) @@ -495,10 +498,6 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, ! old logical massbal_adjust, if non-zero sst_ind = getindex(vars2d, 'sst') - read_atm_file = ( u_ind>0 .or. v_ind>0 .or. tv_ind>0 .or. q_ind>0 .or. & - oz_ind>0 .or. cw_ind>0 .or. tsen_ind>0 .or. ql_ind>0 .or. & - qi_ind>0 .or. prse_ind>0 .or. qr_ind>0 .or. qs_ind>0 .or. qg_ind>0 ) - use_full_hydro = ( ql_ind > 0 .and. qi_ind > 0 .and. & qr_ind > 0 .and. qs_ind > 0 .and. qg_ind > 0 ) @@ -510,24 +509,6 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, ! print *, 'ps: ', ps_ind, ', pst: ', pst_ind, ', sst: ', sst_ind ! endif - ! land sfc DA variables - tmp2m_ind = getindex(vars2d, 't2m') - spfh2m_ind = getindex(vars2d, 'q2m') - soilt1_ind = getindex(vars2d, 'soilt1') - slc1_ind = getindex(vars2d, 'slc1') - soilt2_ind = getindex(vars2d, 'soilt2') - slc2_ind = getindex(vars2d, 'slc2') - soilt3_ind = getindex(vars2d, 'soilt3') - slc3_ind = getindex(vars2d, 'slc3') - soilt4_ind = getindex(vars2d, 'soilt4') - slc4_ind = getindex(vars2d, 'slc4') - - read_sfc_file = ( tmp2m_ind > 0 .or. spfh2m_ind > 0 .or. soilt1_ind > 0 .or. & - slc1_ind > 0 .or. soilt2_ind > 0 .or. slc2_ind > 0 .or. & - soilt3_ind > 0 .or. slc3_ind > 0 .or. soilt4_ind > 0 .or. & - slc4_ind > 0 ) - - ne = 0 ensmemloop: do nanal=nanal1,nanal2 ne = ne + 1 @@ -1044,14 +1025,26 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, if ( read_sfc_file ) then - if ( .not. use_gfs_ncio ) then - write(6,*) 'griddio/griddata for sfc update vars only coded for nc io' - call stop2(23) - endif - if ( reducedgrid ) then - write(6,*) "reducedgrid=T not valid for global_2mDA=T" - call stop2(22) - endif + if ( .not. use_gfs_ncio ) then + write(6,*) 'griddio/griddata for sfc update vars only coded for nc io' + call stop2(23) + endif + if ( reducedgrid ) then + write(6,*) "reducedgrid=T not valid for global_2mDA=T" + call stop2(22) + endif + + ! land sfc DA variables + tmp2m_ind = getindex(vars2d, 't2m') + spfh2m_ind = getindex(vars2d, 'q2m') + soilt1_ind = getindex(vars2d, 'soilt1') + slc1_ind = getindex(vars2d, 'slc1') + soilt2_ind = getindex(vars2d, 'soilt2') + slc2_ind = getindex(vars2d, 'slc2') + soilt3_ind = getindex(vars2d, 'soilt3') + slc3_ind = getindex(vars2d, 'slc3') + soilt4_ind = getindex(vars2d, 'soilt4') + slc4_ind = getindex(vars2d, 'slc4') dset_sfc = open_dataset(filenamesfc) ! read in sfc vars, if requested @@ -3709,323 +3702,10 @@ end subroutine copyfromgrdin end subroutine writegriddata - subroutine writeincrement_2mDA(nanal1,nanal2,vars2d,n2d,grdin_2d,no_inflate_flag) - use netcdf - use params, only: nbackgrounds,fgsfcfileprefixes,incsfcfileprefixes,& - datestring,nhr_anal,write_ensmean,reducedgrid - use mpi - use module_fv3gfs_ncio, only: Dataset, Variable, Dimension, open_dataset,& - read_attribute, close_dataset, get_dim, read_vardata,& - create_dataset, get_idate_from_time_units, & - get_time_units_from_idate, write_vardata, & - write_attribute, quantize_data, has_var, has_attr - implicit none - - integer, intent(in) :: nanal1,nanal2 - character(len=max_varname_length), dimension(n2d), intent(in) :: vars2d - integer, intent(in) :: n2d - ! only the 2d segment of grdin is input - real(r_single), dimension(npts,n2d,nbackgrounds,1), intent(inout) :: grdin_2d - logical, intent(in) :: no_inflate_flag - character(len=500):: filenamein, filenameout - integer(i_kind) :: i, j, nb, ne, nanal - character(len=3) charnanal - type(Dataset) :: dsfg - - integer(i_kind) :: iret - - integer(i_kind) :: tmp2m_ind, spfh2m_ind, soilt1_ind, soilt2_ind, soilt3_ind, & - soilt4_ind,slc1_ind, slc2_ind, slc3_ind, slc4_ind - - ! netcdf things - integer(i_kind) :: dimids(2), ncstart(2), nccount(2) - integer(i_kind) :: ncid_out, lon_dimid, lat_dimid - integer(i_kind) :: lonvarid, latvarid, & - tmp2mvarid, spfh2mvarid, & - soilt1varid, soilt2varid, soilt3varid, soilt4varid, & - slc1varid, slc2varid, slc3varid, slc4varid, & - maskvarid - integer(i_kind) :: iadateout - - ! fixed fields such as lat, lon, levs - real(r_kind),dimension(nlons) :: deglons - real(r_kind),dimension(nlats) :: deglats - - ! soil / snow mask (not fixed) - integer(i_kind), dimension(nlons,nlats) :: mask - - ! increment - real(r_kind), dimension(nlons*nlats) :: inc - real(r_single), allocatable, dimension(:,:) :: inc2d, inc2dout - real(r_kind), allocatable, dimension(:) :: values_1d - real(r_kind), allocatable, dimension(:,:) :: values_2d - - read(datestring,*) iadateout - - ncstart = (/1, 1/) - nccount = (/nlons, nlats/) - - ne = 0 - ensmemloop: do nanal=nanal1,nanal2 - ne = ne + 1 - write(charnanal,'(i3.3)') nanal - backgroundloop: do nb=1,nbackgrounds - - if (nanal == 0 .and. write_ensmean) then - filenamein = trim(adjustl(datapath))//trim(adjustl(fgsfcfileprefixes(nb)))//"ensmean" - filenameout = trim(adjustl(datapath))//trim(adjustl(incsfcfileprefixes(nb)))//"ensmean" - else - 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 - endif - - ! 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)) - dimids = (/ 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, dimids, tmp2mvarid)) - call nccheck_incr(nf90_def_var(ncid_out, "spfh2m_inc", nf90_real, dimids, spfh2mvarid)) - call nccheck_incr(nf90_def_var(ncid_out, "soilt1_inc", nf90_real, dimids, soilt1varid)) - call nccheck_incr(nf90_def_var(ncid_out, "soilt2_inc", nf90_real, dimids, soilt2varid)) - call nccheck_incr(nf90_def_var(ncid_out, "soilt3_inc", nf90_real, dimids, soilt3varid)) - call nccheck_incr(nf90_def_var(ncid_out, "soilt4_inc", nf90_real, dimids, soilt4varid)) - call nccheck_incr(nf90_def_var(ncid_out, "slc1_inc", nf90_real, dimids, slc1varid)) - call nccheck_incr(nf90_def_var(ncid_out, "slc2_inc", nf90_real, dimids, slc2varid)) - call nccheck_incr(nf90_def_var(ncid_out, "slc3_inc", nf90_real, dimids, slc3varid)) - call nccheck_incr(nf90_def_var(ncid_out, "slc4_inc", nf90_real, dimids, slc4varid)) - call nccheck_incr(nf90_def_var(ncid_out, "soilsnow_mask", nf90_int, dimids, 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_2mDA")) - 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)) - - tmp2m_ind = getindex(vars2d, 't2m') !< indices in the state or control var arrays - spfh2m_ind = getindex(vars2d, 'q2m') - soilt1_ind = getindex(vars2d, 'soilt1') - slc1_ind = getindex(vars2d, 'slc1') - soilt2_ind = getindex(vars2d, 'soilt2') - slc2_ind = getindex(vars2d, 'slc2') - soilt3_ind = getindex(vars2d, 'soilt3') - slc3_ind = getindex(vars2d, 'slc3') - soilt4_ind = getindex(vars2d, 'soilt4') - slc4_ind = getindex(vars2d, 'slc4') - - dsfg = open_dataset(filenamein) - - ! 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, 'slc1', 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, count = nccount)) - - allocate(inc2d(nlons,nlats)) - allocate(inc2dout(nlons,nlats)) - - ! tmp2m increment - inc(:) = zero - if (tmp2m_ind > 0) then - call copyfromgrdin(grdin_2d(:,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, count = nccount)) - ! spfh2m increment - inc(:) = zero - if (spfh2m_ind > 0) then - call copyfromgrdin(grdin_2d(:,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, count = nccount)) - ! soilt1 increment - inc(:) = zero - if (soilt1_ind > 0) then - call copyfromgrdin(grdin_2d(:,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, count = nccount)) - ! soilt2 increment - inc(:) = zero - if (soilt2_ind > 0) then - call copyfromgrdin(grdin_2d(:,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, count = nccount)) - ! soilt3 increment - inc(:) = zero - if (soilt3_ind > 0) then - call copyfromgrdin(grdin_2d(:,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, count = nccount)) - ! soilt4 increment - inc(:) = zero - if (soilt4_ind > 0) then - call copyfromgrdin(grdin_2d(:,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, count = nccount)) - ! slc1 increment - inc(:) = zero - if (slc1_ind > 0) then - call copyfromgrdin(grdin_2d(:,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, count = nccount)) - ! slc2 increment - inc(:) = zero - if (slc2_ind > 0) then - call copyfromgrdin(grdin_2d(:,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, count = nccount)) - ! slc3 increment - inc(:) = zero - if (slc3_ind > 0) then - call copyfromgrdin(grdin_2d(:,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, count = nccount)) - ! slc4 increment - inc(:) = zero - if (slc4_ind > 0) then - call copyfromgrdin(grdin_2d(:,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, count = nccount)) - - ! deallocate things - deallocate(inc2d,inc2dout) - - end do backgroundloop ! loop over backgrounds to read in - end do ensmemloop ! loop over ens members to read in - - return - - contains -! copying to grdin (calling regtoreduced if reduced grid) - subroutine copyfromgrdin(grdin, field) - implicit none - - real(r_single), dimension(:), intent(in) :: grdin - real(r_kind), dimension(:), intent(inout) :: field - - if (reducedgrid) then - call reducedtoreg(grdin, field) - else - field = grdin - endif - - end subroutine copyfromgrdin - - end subroutine writeincrement_2mDA - subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate_flag) use netcdf use params, only: nbackgrounds,incfileprefixes,fgfileprefixes,reducedgrid,& - datestring,nhr_anal,write_ensmean + datestring,nhr_anal,write_ensmean, fgsfcfileprefixes,incsfcfileprefixes use constants, only: grav use mpi use module_fv3gfs_ncio, only: Dataset, Variable, Dimension, open_dataset,& @@ -4061,7 +3741,12 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin, 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 + tvarid, sphumvarid, liqwatvarid, o3varid, icvarid, & + 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 @@ -4073,337 +3758,598 @@ subroutine writeincrement(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin, ! 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_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 - use_full_hydro = .false. - clip = tiny_r_kind - read(datestring,*) iadateout + call set_ncio_file_flags(vars3d, n3d, vars2d, n2d, write_sfc_file, write_atm_file) + + if ( write_atm_file) then + use_full_hydro = .false. + clip = tiny_r_kind + read(datestring,*) iadateout + + ncstart = (/1, 1, 1/) + nccount = (/nlons, nlats, nlevs/) + + ne = 0 + ensmemloop: do nanal=nanal1,nanal2 + ne = ne + 1 + write(charnanal,'(i3.3)') nanal + backgroundloop: do nb=1,nbackgrounds + + if (nanal == 0 .and. write_ensmean) then + filenameout = trim(adjustl(datapath))//trim(adjustl(incfileprefixes(nb)))//"ensmean" + filenamein = trim(adjustl(datapath))//trim(adjustl(fgfileprefixes(nb)))//"ensmean" + else + if(no_inflate_flag) then + filenameout = trim(adjustl(datapath))//trim(adjustl(incfileprefixes(nb)))//"nimem"//charnanal + else + filenameout = trim(adjustl(datapath))//trim(adjustl(incfileprefixes(nb)))//"mem"//charnanal + end if + filenamein = trim(adjustl(datapath))//trim(adjustl(fgfileprefixes(nb)))//"mem"//charnanal + endif - ncstart = (/1, 1, 1/) - nccount = (/nlons, nlats, nlevs/) + ! 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, "lon", nlons, lon_dimid)) + call nccheck_incr(nf90_def_dim(ncid_out, "lat", nlats, lat_dimid)) + call nccheck_incr(nf90_def_dim(ncid_out, "lev", nlevs, lev_dimid)) + call nccheck_incr(nf90_def_dim(ncid_out, "ilev", nlevs+1, ilev_dimid)) + dimids3 = (/ lon_dimid, lat_dimid, lev_dimid /) + ! create variables + call nccheck_incr(nf90_def_var(ncid_out, "lon", nf90_real, (/lon_dimid/), lonvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "lat", nf90_real, (/lat_dimid/), latvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "lev", nf90_real, (/lev_dimid/), levvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "pfull", nf90_real, (/lev_dimid/), pfullvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "ilev", nf90_real, (/ilev_dimid/), ilevvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "hyai", nf90_real, (/ilev_dimid/), hyaivarid)) + call nccheck_incr(nf90_def_var(ncid_out, "hybi", nf90_real, (/ilev_dimid/), hybivarid)) + call nccheck_incr(nf90_def_var(ncid_out, "u_inc", nf90_real, dimids3, uvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "v_inc", nf90_real, dimids3, vvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "delp_inc", nf90_real, dimids3, delpvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "delz_inc", nf90_real, dimids3, delzvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "T_inc", nf90_real, dimids3, tvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "sphum_inc", nf90_real, dimids3, sphumvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "liq_wat_inc", nf90_real, dimids3, liqwatvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "o3mr_inc", nf90_real, dimids3, o3varid)) + call nccheck_incr(nf90_def_var(ncid_out, "icmr_inc", nf90_real, dimids3, icvarid)) + ! 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 analysis 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)) + + u_ind = getindex(vars3d, 'u') !< indices in the control var arrays + v_ind = getindex(vars3d, 'v') ! U and V (3D) + tv_ind = getindex(vars3d, 'tv') ! Tv (3D) + q_ind = getindex(vars3d, 'q') ! Q (3D) + ps_ind = getindex(vars2d, 'ps') ! Ps (2D) + oz_ind = getindex(vars3d, 'oz') ! Oz (3D) + cw_ind = getindex(vars3d, 'cw') ! CW (3D) + ql_ind = getindex(vars3d, 'ql') ! QL (3D) + qi_ind = getindex(vars3d, 'qi') ! QI (3D) + qr_ind = getindex(vars3d, 'qr') ! QR (3D) + qs_ind = getindex(vars3d, 'qs') ! QS (3D) + qg_ind = getindex(vars3d, 'qg') ! QG (3D) + pst_ind = getindex(vars2d, 'pst') ! Ps tendency (2D) // equivalent of + ! old logical massbal_adjust, if non-zero + use_full_hydro = ( ql_ind > 0 .and. qi_ind > 0 .and. & + qr_ind > 0 .and. qs_ind > 0 .and. qg_ind > 0 ) - ne = 0 - ensmemloop: do nanal=nanal1,nanal2 - ne = ne + 1 - write(charnanal,'(i3.3)') nanal - backgroundloop: do nb=1,nbackgrounds + dsfg = open_dataset(filenamein) + call read_attribute(dsfg, 'ak', values_1d,errcode=iret) + if (iret /= 0) then + print *,'error reading ak' + call stop2(29) + endif + do k=1,nlevs+1 + ! k=1 in values_1d is model top, flip so k=1 in ak is bottom + ak(nlevs-k+2) = 0.01_r_kind*values_1d(k) + enddo + call read_attribute(dsfg, 'bk', values_1d,errcode=iret) + if (iret /= 0) then + print *,'error reading bk' + call stop2(29) + endif + do k=1,nlevs+1 + ! k=1 in values_1d is model top, flip so k=1 in ak is bottom + bk(nlevs-k+2) = values_1d(k) + enddo - if (nanal == 0 .and. write_ensmean) then - filenameout = trim(adjustl(datapath))//trim(adjustl(incfileprefixes(nb)))//"ensmean" - filenamein = trim(adjustl(datapath))//trim(adjustl(fgfileprefixes(nb)))//"ensmean" - else - if(no_inflate_flag) then - filenameout = trim(adjustl(datapath))//trim(adjustl(incfileprefixes(nb)))//"nimem"//charnanal - else - filenameout = trim(adjustl(datapath))//trim(adjustl(incfileprefixes(nb)))//"mem"//charnanal - end if - filenamein = trim(adjustl(datapath))//trim(adjustl(fgfileprefixes(nb)))//"mem"//charnanal - endif + ! levels + do k=1,nlevs + levsout(k) = float(k) + ilevsout(k) = float(k) + end do + ilevsout(nlevs+1) = float(nlevs+1) + + ! 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 - ! create the output netCDF increment file - call nccheck_incr(nf90_create(path=trim(filenameout), cmode=nf90_netcdf4, ncid=ncid_out)) + call nccheck_incr(nf90_put_var(ncid_out, latvarid, deglats, & + start = (/1/), count = (/nlats/))) + + ! write to file + call nccheck_incr(nf90_put_var(ncid_out, levvarid, sngl(levsout), & + start = (/1/), count = (/nlevs/))) + ! pfull + call nccheck_incr(nf90_put_var(ncid_out, pfullvarid, sngl(levsout), & + start = (/1/), count = (/nlevs/))) + ! ilev + call nccheck_incr(nf90_put_var(ncid_out, ilevvarid, sngl(ilevsout), & + start = (/1/), count = (/nlevs+1/))) + ! hyai + call nccheck_incr(nf90_put_var(ncid_out, hyaivarid, sngl(ilevsout), & + start = (/1/), count = (/nlevs+1/))) + ! hybi + call nccheck_incr(nf90_put_var(ncid_out, hybivarid, sngl(ilevsout), & + start = (/1/), count = (/nlevs+1/))) + + allocate(inc3d(nlons,nlats,nccount(3))) + allocate(inc3d2(nlons,nlats,nccount(3))) + allocate(inc3dout(nlons,nlats,nccount(3))) + ! u increment + do k=1,nlevs + krev = nlevs-k+1 + inc(:) = zero + if (u_ind > 0) then + call copyfromgrdin(grdin(:,levels(u_ind-1) + krev,nb,ne),inc) + endif + inc3d(:,:,k) = reshape(inc,(/nlons,nlats/)) + end do + do j=1,nlats + inc3dout(:,nlats-j+1,:) = inc3d(:,j,:) + end do + if (should_zero_increments_for('u_inc')) inc3dout = zero + call nccheck_incr(nf90_put_var(ncid_out, uvarid, sngl(inc3dout), & + start = ncstart, count = nccount)) + ! v increment + do k=1,nlevs + krev = nlevs-k+1 + inc(:) = zero + if (u_ind > 0) then + call copyfromgrdin(grdin(:,levels(v_ind-1) + krev,nb,ne),inc) + endif + inc3d(:,:,k) = reshape(inc,(/nlons,nlats/)) + end do + do j=1,nlats + inc3dout(:,nlats-j+1,:) = inc3d(:,j,:) + end do + if (should_zero_increments_for('v_inc')) inc3dout = zero + call nccheck_incr(nf90_put_var(ncid_out, vvarid, sngl(inc3dout), & + start = ncstart, count = nccount)) - ! create dimensions based on analysis resolution, not guess - call nccheck_incr(nf90_def_dim(ncid_out, "lon", nlons, lon_dimid)) - call nccheck_incr(nf90_def_dim(ncid_out, "lat", nlats, lat_dimid)) - call nccheck_incr(nf90_def_dim(ncid_out, "lev", nlevs, lev_dimid)) - call nccheck_incr(nf90_def_dim(ncid_out, "ilev", nlevs+1, ilev_dimid)) - dimids3 = (/ lon_dimid, lat_dimid, lev_dimid /) - ! create variables - call nccheck_incr(nf90_def_var(ncid_out, "lon", nf90_real, (/lon_dimid/), lonvarid)) - call nccheck_incr(nf90_def_var(ncid_out, "lat", nf90_real, (/lat_dimid/), latvarid)) - call nccheck_incr(nf90_def_var(ncid_out, "lev", nf90_real, (/lev_dimid/), levvarid)) - call nccheck_incr(nf90_def_var(ncid_out, "pfull", nf90_real, (/lev_dimid/), pfullvarid)) - call nccheck_incr(nf90_def_var(ncid_out, "ilev", nf90_real, (/ilev_dimid/), ilevvarid)) - call nccheck_incr(nf90_def_var(ncid_out, "hyai", nf90_real, (/ilev_dimid/), hyaivarid)) - call nccheck_incr(nf90_def_var(ncid_out, "hybi", nf90_real, (/ilev_dimid/), hybivarid)) - call nccheck_incr(nf90_def_var(ncid_out, "u_inc", nf90_real, dimids3, uvarid)) - call nccheck_incr(nf90_def_var(ncid_out, "v_inc", nf90_real, dimids3, vvarid)) - call nccheck_incr(nf90_def_var(ncid_out, "delp_inc", nf90_real, dimids3, delpvarid)) - call nccheck_incr(nf90_def_var(ncid_out, "delz_inc", nf90_real, dimids3, delzvarid)) - call nccheck_incr(nf90_def_var(ncid_out, "T_inc", nf90_real, dimids3, tvarid)) - call nccheck_incr(nf90_def_var(ncid_out, "sphum_inc", nf90_real, dimids3, sphumvarid)) - call nccheck_incr(nf90_def_var(ncid_out, "liq_wat_inc", nf90_real, dimids3, liqwatvarid)) - call nccheck_incr(nf90_def_var(ncid_out, "o3mr_inc", nf90_real, dimids3, o3varid)) - call nccheck_incr(nf90_def_var(ncid_out, "icmr_inc", nf90_real, dimids3, icvarid)) - ! 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 analysis 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)) + ! delp increment + psinc(:) = zero + if (ps_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d) + ps_ind,nb,ne),psinc) + endif + do k=1,nlevs + krev = nlevs-k+1 + inc(:) = zero + inc = psinc*(bk(krev)-bk(krev+1))*100_r_kind + inc3d(:,:,k) = reshape(inc,(/nlons,nlats/)) + end do + do j=1,nlats + inc3dout(:,nlats-j+1,:) = inc3d(:,j,:) + end do + if (should_zero_increments_for('delp_inc')) inc3dout = zero + call nccheck_incr(nf90_put_var(ncid_out, delpvarid, sngl(inc3dout), & + start = ncstart, count = nccount)) + + ! sphum increment + allocate(tmp(nlons,nlats,nccount(3)),tv(nlons,nlats,nccount(3)),q(nlons,nlats,nccount(3))) + allocate(tvanl(nlons,nlats,nccount(3)),tmpanl(nlons,nlats,nccount(3)),qanl(nlons,nlats,nccount(3))) + call read_vardata(dsfg, 'spfh', q, ncstart=ncstart, nccount=nccount, errcode=iret) + if (iret /= 0) then + print *,'error reading spfh' + call stop2(29) + endif + do k=1,nlevs + krev = nlevs-k+1 + inc(:) = zero + if (q_ind > 0) then + call copyfromgrdin(grdin(:,levels(q_ind-1) + krev,nb,ne),inc) + endif + inc3d(:,:,k) = reshape(inc,(/nlons,nlats/)) + qanl(:,:,k) = q(:,:,k) + inc3d(:,:,k) + end do + if (cliptracers) where (qanl < clip) qanl = clip + inc3d = qanl - q + do j=1,nlats + inc3dout(:,nlats-j+1,:) = inc3d(:,j,:) + end do + if (should_zero_increments_for('sphum_inc')) inc3dout = zero + call nccheck_incr(nf90_put_var(ncid_out, sphumvarid, sngl(inc3dout), & + start = ncstart, count = nccount)) - u_ind = getindex(vars3d, 'u') !< indices in the control var arrays - v_ind = getindex(vars3d, 'v') ! U and V (3D) - tv_ind = getindex(vars3d, 'tv') ! Tv (3D) - q_ind = getindex(vars3d, 'q') ! Q (3D) - ps_ind = getindex(vars2d, 'ps') ! Ps (2D) - oz_ind = getindex(vars3d, 'oz') ! Oz (3D) - cw_ind = getindex(vars3d, 'cw') ! CW (3D) - ql_ind = getindex(vars3d, 'ql') ! QL (3D) - qi_ind = getindex(vars3d, 'qi') ! QI (3D) - qr_ind = getindex(vars3d, 'qr') ! QR (3D) - qs_ind = getindex(vars3d, 'qs') ! QS (3D) - qg_ind = getindex(vars3d, 'qg') ! QG (3D) - pst_ind = getindex(vars2d, 'pst') ! Ps tendency (2D) // equivalent of - ! old logical massbal_adjust, if non-zero - use_full_hydro = ( ql_ind > 0 .and. qi_ind > 0 .and. & - qr_ind > 0 .and. qs_ind > 0 .and. qg_ind > 0 ) + ! t increment + call read_vardata(dsfg, 'tmp', tmp, ncstart=ncstart, nccount=nccount, errcode=iret) + if (iret /= 0) then + print *,'error reading tmp' + call stop2(29) + endif + tv = tmp * ( 1.0 + fv*q) + do k=1,nlevs + krev = nlevs-k+1 + inc(:) = zero + if (tv_ind > 0) then + call copyfromgrdin(grdin(:,levels(tv_ind-1) + krev,nb,ne),inc) + endif + inc3d(:,:,k) = reshape(inc,(/nlons,nlats/)) + tvanl(:,:,k) = tv(:,:,k) + inc3d(:,:,k) + tmpanl(:,:,k) = tvanl(:,:,k)/(1. + fv*qanl(:,:,k)) + end do + inc3d = tmpanl - tmp + do j=1,nlats + inc3dout(:,nlats-j+1,:) = inc3d(:,j,:) + end do + if (should_zero_increments_for('T_inc')) inc3dout = zero + call nccheck_incr(nf90_put_var(ncid_out, tvarid, sngl(inc3dout), & + start = ncstart, count = nccount)) - dsfg = open_dataset(filenamein) - call read_attribute(dsfg, 'ak', values_1d,errcode=iret) - if (iret /= 0) then - print *,'error reading ak' - call stop2(29) - endif - do k=1,nlevs+1 - ! k=1 in values_1d is model top, flip so k=1 in ak is bottom - ak(nlevs-k+2) = 0.01_r_kind*values_1d(k) - enddo - call read_attribute(dsfg, 'bk', values_1d,errcode=iret) - if (iret /= 0) then - print *,'error reading bk' - call stop2(29) - endif - do k=1,nlevs+1 - ! k=1 in values_1d is model top, flip so k=1 in ak is bottom - bk(nlevs-k+2) = values_1d(k) - enddo + ! delz increment + inc3d(:,:,:) = zero + if (has_var(dsfg,'delz')) then + allocate(delzb(nlons*nlats)) + call read_vardata(dsfg,'pressfc',values_2d,errcode=iret) + if (allocated(psges)) deallocate(psges) + allocate(psges(nlons*nlats)) + psges = reshape(values_2d,(/nlons*nlats/)) + vg = psges + (psinc*100_r_kind) + do k=1,nlevs + krev = nlevs-k+1 ! k=1 is model top + ug=(rd/grav)*reshape(tvanl(:,:,k),(/nlons*nlats/)) + ! ps in Pa here, need to multiply ak by 100. + ! calculate ug (analysis delz) so it is negative. + ! (note that ak,bk are already reversed to go from bottom to top) + ug=ug*log((100_r_kind*ak(krev+1)+bk(krev+1)*vg)/(100_r_kind*ak(krev)+bk(krev)*vg)) + ! ug is hydrostatic analysis delz inferred from analysis ps,Tv + ! delzb is hydrostatic background delz inferred from background ps,Tv + delzb=(rd/grav)*reshape(tv(:,:,k),(/nlons*nlats/)) + delzb=delzb*log((100_r_kind*ak(krev+1)+bk(krev+1)*psges)/(100_r_kind*ak(krev)+bk(krev)*psges)) + inc3d(:,:,k)=reshape(ug-delzb,(/nlons,nlats/)) + end do + end if + do j=1,nlats + inc3dout(:,nlats-j+1,:) = inc3d(:,j,:) + end do + if (should_zero_increments_for('delz_inc')) inc3dout = zero + call nccheck_incr(nf90_put_var(ncid_out, delzvarid, sngl(inc3dout), & + start = ncstart, count = nccount)) - ! levels - do k=1,nlevs - levsout(k) = float(k) - ilevsout(k) = float(k) - end do - ilevsout(nlevs+1) = float(nlevs+1) + ! o3mr increment + do k=1,nlevs + krev = nlevs-k+1 + inc(:) = zero + if (oz_ind > 0) then + call copyfromgrdin(grdin(:,levels(oz_ind-1) + krev,nb,ne),inc) + endif + inc3d(:,:,k) = reshape(inc,(/nlons,nlats/)) + end do + do j=1,nlats + inc3dout(:,nlats-j+1,:) = inc3d(:,j,:) + end do + if (should_zero_increments_for('o3mr_inc')) inc3dout = zero + call nccheck_incr(nf90_put_var(ncid_out, o3varid, sngl(inc3dout), & + start = ncstart, count = nccount)) - ! 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/))) + ! liq wat increment + ! icmr increment + do k=1,nlevs + krev = nlevs-k+1 + ug = zero + if (cw_ind > 0) then + call copyfromgrdin(grdin(:,levels(cw_ind-1)+krev,nb,ne),ug) + end if + if (imp_physics == 11) then + work = -r0_05 * (reshape(tmpanl(:,:,k),(/nlons*nlats/)) - t0c) + do i=1,nlons*nlats + work(i) = max(zero,work(i)) + work(i) = min(one,work(i)) + enddo + vg = ug * work ! cloud ice + ug = ug * (one - work) ! cloud water + inc3d2(:,:,k) = reshape(vg,(/nlons,nlats/)) + endif + inc3d(:,:,k) = reshape(ug,(/nlons,nlats/)) + enddo + do j=1,nlats + inc3dout(:,nlats-j+1,:) = inc3d(:,j,:) + end do + if (should_zero_increments_for('liq_wat_inc')) inc3dout = zero + call nccheck_incr(nf90_put_var(ncid_out, liqwatvarid, sngl(inc3dout), & + start = ncstart, count = nccount)) + do j=1,nlats + inc3dout(:,nlats-j+1,:) = inc3d2(:,j,:) + end do + if (should_zero_increments_for('icmr_inc')) inc3dout = zero + call nccheck_incr(nf90_put_var(ncid_out, icvarid, sngl(inc3dout), & + start = ncstart, count = nccount)) - 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 close_dataset(dsfg,errcode=iret) + if (iret/=0) then + write(6,*)'gridio/writeincrement gfs model: problem closing netcdf fg dataset, iret=',iret + call stop2(23) + endif + ! deallocate things + deallocate(inc3d,inc3d2,inc3dout) + deallocate(tmp,tv,q,tmpanl,tvanl,qanl) + if (allocated(delzb)) deallocate(delzb) + if (allocated(psges)) deallocate(psges) - call nccheck_incr(nf90_put_var(ncid_out, latvarid, deglats, & - start = (/1/), count = (/nlats/))) - - ! write to file - call nccheck_incr(nf90_put_var(ncid_out, levvarid, sngl(levsout), & - start = (/1/), count = (/nlevs/))) - ! pfull - call nccheck_incr(nf90_put_var(ncid_out, pfullvarid, sngl(levsout), & - start = (/1/), count = (/nlevs/))) - ! ilev - call nccheck_incr(nf90_put_var(ncid_out, ilevvarid, sngl(ilevsout), & - start = (/1/), count = (/nlevs+1/))) - ! hyai - call nccheck_incr(nf90_put_var(ncid_out, hyaivarid, sngl(ilevsout), & - start = (/1/), count = (/nlevs+1/))) - ! hybi - call nccheck_incr(nf90_put_var(ncid_out, hybivarid, sngl(ilevsout), & - start = (/1/), count = (/nlevs+1/))) + end do backgroundloop ! loop over backgrounds to read in + end do ensmemloop ! loop over ens members to read in - allocate(inc3d(nlons,nlats,nccount(3))) - allocate(inc3d2(nlons,nlats,nccount(3))) - allocate(inc3dout(nlons,nlats,nccount(3))) - ! u increment - do k=1,nlevs - krev = nlevs-k+1 + endif ! write_atm_file + + if (write_sfc_file) then + + ne = 0 + sfcensmemloop: do nanal=nanal1,nanal2 + ne = ne + 1 + write(charnanal,'(i3.3)') nanal + sfcbackgroundloop: do nb=1,nbackgrounds + + if (nanal == 0 .and. write_ensmean) then + filenamein = trim(adjustl(datapath))//trim(adjustl(fgsfcfileprefixes(nb)))//"ensmean" + filenameout = trim(adjustl(datapath))//trim(adjustl(incsfcfileprefixes(nb)))//"ensmean" + else + 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 + endif + + ! 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)) + ! 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, dimids3(1:2), tmp2mvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "spfh2m_inc", nf90_real, dimids3(1:2), spfh2mvarid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilt1_inc", nf90_real, dimids3(1:2), soilt1varid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilt2_inc", nf90_real, dimids3(1:2), soilt2varid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilt3_inc", nf90_real, dimids3(1:2), soilt3varid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilt4_inc", nf90_real, dimids3(1:2), soilt4varid)) + call nccheck_incr(nf90_def_var(ncid_out, "slc1_inc", nf90_real, dimids3(1:2), slc1varid)) + call nccheck_incr(nf90_def_var(ncid_out, "slc2_inc", nf90_real, dimids3(1:2), slc2varid)) + call nccheck_incr(nf90_def_var(ncid_out, "slc3_inc", nf90_real, dimids3(1:2), slc3varid)) + call nccheck_incr(nf90_def_var(ncid_out, "slc4_inc", nf90_real, dimids3(1:2), slc4varid)) + call nccheck_incr(nf90_def_var(ncid_out, "soilsnow_mask", nf90_int, dimids3(1:2), 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)) + + tmp2m_ind = getindex(vars2d, 't2m') !< indices in the state or control var arrays + spfh2m_ind = getindex(vars2d, 'q2m') + soilt1_ind = getindex(vars2d, 'soilt1') + slc1_ind = getindex(vars2d, 'slc1') + soilt2_ind = getindex(vars2d, 'soilt2') + slc2_ind = getindex(vars2d, 'slc2') + soilt3_ind = getindex(vars2d, 'soilt3') + slc3_ind = getindex(vars2d, 'slc3') + soilt4_ind = getindex(vars2d, 'soilt4') + slc4_ind = getindex(vars2d, 'slc4') + + dsfg = open_dataset(filenamein) + + ! 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, 'slc1', 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 (u_ind > 0) then - call copyfromgrdin(grdin(:,levels(u_ind-1) + krev,nb,ne),inc) + if (tmp2m_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d) + tmp2m_ind,nb,ne),inc) endif - inc3d(:,:,k) = reshape(inc,(/nlons,nlats/)) - end do - do j=1,nlats - inc3dout(:,nlats-j+1,:) = inc3d(:,j,:) - end do - if (should_zero_increments_for('u_inc')) inc3dout = zero - call nccheck_incr(nf90_put_var(ncid_out, uvarid, sngl(inc3dout), & - start = ncstart, count = nccount)) - ! v increment - do k=1,nlevs - krev = nlevs-k+1 + 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 (u_ind > 0) then - call copyfromgrdin(grdin(:,levels(v_ind-1) + krev,nb,ne),inc) + if (spfh2m_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d)+spfh2m_ind,nb,ne),inc) endif - inc3d(:,:,k) = reshape(inc,(/nlons,nlats/)) - end do - do j=1,nlats - inc3dout(:,nlats-j+1,:) = inc3d(:,j,:) - end do - if (should_zero_increments_for('v_inc')) inc3dout = zero - call nccheck_incr(nf90_put_var(ncid_out, vvarid, sngl(inc3dout), & - start = ncstart, count = nccount)) - - ! delp increment - psinc(:) = zero - if (ps_ind > 0) then - call copyfromgrdin(grdin(:,levels(n3d) + ps_ind,nb,ne),psinc) - endif - do k=1,nlevs - krev = nlevs-k+1 + 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 - inc = psinc*(bk(krev)-bk(krev+1))*100_r_kind - inc3d(:,:,k) = reshape(inc,(/nlons,nlats/)) - end do - do j=1,nlats - inc3dout(:,nlats-j+1,:) = inc3d(:,j,:) - end do - if (should_zero_increments_for('delp_inc')) inc3dout = zero - call nccheck_incr(nf90_put_var(ncid_out, delpvarid, sngl(inc3dout), & - start = ncstart, count = nccount)) - - ! sphum increment - allocate(tmp(nlons,nlats,nccount(3)),tv(nlons,nlats,nccount(3)),q(nlons,nlats,nccount(3))) - allocate(tvanl(nlons,nlats,nccount(3)),tmpanl(nlons,nlats,nccount(3)),qanl(nlons,nlats,nccount(3))) - call read_vardata(dsfg, 'spfh', q, ncstart=ncstart, nccount=nccount, errcode=iret) - if (iret /= 0) then - print *,'error reading spfh' - call stop2(29) - endif - do k=1,nlevs - krev = nlevs-k+1 + 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 (q_ind > 0) then - call copyfromgrdin(grdin(:,levels(q_ind-1) + krev,nb,ne),inc) + if (soilt2_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d)+soilt2_ind,nb,ne),inc) endif - inc3d(:,:,k) = reshape(inc,(/nlons,nlats/)) - qanl(:,:,k) = q(:,:,k) + inc3d(:,:,k) - end do - if (cliptracers) where (qanl < clip) qanl = clip - inc3d = qanl - q - do j=1,nlats - inc3dout(:,nlats-j+1,:) = inc3d(:,j,:) - end do - if (should_zero_increments_for('sphum_inc')) inc3dout = zero - call nccheck_incr(nf90_put_var(ncid_out, sphumvarid, sngl(inc3dout), & - start = ncstart, count = nccount)) - - ! t increment - call read_vardata(dsfg, 'tmp', tmp, ncstart=ncstart, nccount=nccount, errcode=iret) - if (iret /= 0) then - print *,'error reading tmp' - call stop2(29) - endif - tv = tmp * ( 1.0 + fv*q) - do k=1,nlevs - krev = nlevs-k+1 + 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 (tv_ind > 0) then - call copyfromgrdin(grdin(:,levels(tv_ind-1) + krev,nb,ne),inc) + if (soilt3_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d)+soilt3_ind,nb,ne),inc) endif - inc3d(:,:,k) = reshape(inc,(/nlons,nlats/)) - tvanl(:,:,k) = tv(:,:,k) + inc3d(:,:,k) - tmpanl(:,:,k) = tvanl(:,:,k)/(1. + fv*qanl(:,:,k)) - end do - inc3d = tmpanl - tmp - do j=1,nlats - inc3dout(:,nlats-j+1,:) = inc3d(:,j,:) - end do - if (should_zero_increments_for('T_inc')) inc3dout = zero - call nccheck_incr(nf90_put_var(ncid_out, tvarid, sngl(inc3dout), & - start = ncstart, count = nccount)) - - ! delz increment - inc3d(:,:,:) = zero - if (has_var(dsfg,'delz')) then - allocate(delzb(nlons*nlats)) - call read_vardata(dsfg,'pressfc',values_2d,errcode=iret) - if (allocated(psges)) deallocate(psges) - allocate(psges(nlons*nlats)) - psges = reshape(values_2d,(/nlons*nlats/)) - vg = psges + (psinc*100_r_kind) - do k=1,nlevs - krev = nlevs-k+1 ! k=1 is model top - ug=(rd/grav)*reshape(tvanl(:,:,k),(/nlons*nlats/)) - ! ps in Pa here, need to multiply ak by 100. - ! calculate ug (analysis delz) so it is negative. - ! (note that ak,bk are already reversed to go from bottom to top) - ug=ug*log((100_r_kind*ak(krev+1)+bk(krev+1)*vg)/(100_r_kind*ak(krev)+bk(krev)*vg)) - ! ug is hydrostatic analysis delz inferred from analysis ps,Tv - ! delzb is hydrostatic background delz inferred from background ps,Tv - delzb=(rd/grav)*reshape(tv(:,:,k),(/nlons*nlats/)) - delzb=delzb*log((100_r_kind*ak(krev+1)+bk(krev+1)*psges)/(100_r_kind*ak(krev)+bk(krev)*psges)) - inc3d(:,:,k)=reshape(ug-delzb,(/nlons,nlats/)) + 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 - end if - do j=1,nlats - inc3dout(:,nlats-j+1,:) = inc3d(:,j,:) - end do - if (should_zero_increments_for('delz_inc')) inc3dout = zero - call nccheck_incr(nf90_put_var(ncid_out, delzvarid, sngl(inc3dout), & - start = ncstart, count = nccount)) - - ! o3mr increment - do k=1,nlevs - krev = nlevs-k+1 + call nccheck_incr(nf90_put_var(ncid_out, soilt3varid, sngl(inc2dout), & + start = ncstart(1:2), count = nccount(1:2))) + ! soilt4 increment inc(:) = zero - if (oz_ind > 0) then - call copyfromgrdin(grdin(:,levels(oz_ind-1) + krev,nb,ne),inc) + if (soilt4_ind > 0) then + call copyfromgrdin(grdin(:,levels(n3d)+soilt4_ind,nb,ne),inc) endif - inc3d(:,:,k) = reshape(inc,(/nlons,nlats/)) - end do - do j=1,nlats - inc3dout(:,nlats-j+1,:) = inc3d(:,j,:) - end do - if (should_zero_increments_for('o3mr_inc')) inc3dout = zero - call nccheck_incr(nf90_put_var(ncid_out, o3varid, sngl(inc3dout), & - start = ncstart, count = nccount)) + 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))) - ! liq wat increment - ! icmr increment - do k=1,nlevs - krev = nlevs-k+1 - ug = zero - if (cw_ind > 0) then - call copyfromgrdin(grdin(:,levels(cw_ind-1)+krev,nb,ne),ug) - end if - if (imp_physics == 11) then - work = -r0_05 * (reshape(tmpanl(:,:,k),(/nlons*nlats/)) - t0c) - do i=1,nlons*nlats - work(i) = max(zero,work(i)) - work(i) = min(one,work(i)) - enddo - vg = ug * work ! cloud ice - ug = ug * (one - work) ! cloud water - inc3d2(:,:,k) = reshape(vg,(/nlons,nlats/)) + call close_dataset(dsfg,errcode=iret) + if (iret/=0) then + write(6,*)'gridio/writeincrement gfs model: problem closing netcdf sfc fg dataset, iret=',iret + call stop2(23) endif - inc3d(:,:,k) = reshape(ug,(/nlons,nlats/)) - enddo - do j=1,nlats - inc3dout(:,nlats-j+1,:) = inc3d(:,j,:) - end do - if (should_zero_increments_for('liq_wat_inc')) inc3dout = zero - call nccheck_incr(nf90_put_var(ncid_out, liqwatvarid, sngl(inc3dout), & - start = ncstart, count = nccount)) - do j=1,nlats - inc3dout(:,nlats-j+1,:) = inc3d2(:,j,:) - end do - if (should_zero_increments_for('icmr_inc')) inc3dout = zero - call nccheck_incr(nf90_put_var(ncid_out, icvarid, sngl(inc3dout), & - start = ncstart, count = nccount)) + ! deallocate things + deallocate(inc2d,inc2dout) - ! deallocate things - deallocate(inc3d,inc3d2,inc3dout) - deallocate(tmp,tv,q,tmpanl,tvanl,qanl) - deallocate(delzb,psges) + end do sfcbackgroundloop ! loop over backgrounds to read in + end do sfcensmemloop ! loop over ens members to read in - end do backgroundloop ! loop over backgrounds to read in - end do ensmemloop ! loop over ens members to read in + + endif ! write_sfc_file return @@ -4883,6 +4829,68 @@ end subroutine copyfromgrdin end subroutine writeincrement_pnc + subroutine set_ncio_file_flags(vars3d, n3d, vars2d, n2d, sfc_file, atm_file) + ! determine if variables are in sfc and/or atm file, for ncio case. + character(len=max_varname_length), dimension(n2d), intent(in) :: vars2d + character(len=max_varname_length), dimension(n3d), intent(in) :: vars3d + integer, intent(in) :: n2d, n3d + logical, intent(out) :: sfc_file, atm_file + + 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 + 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 + +!do i = 1, nc2d + !if (getindex(vars2d_supported, cvars2d(i))<0) then + !if (nproc .eq. 0) then + + ! atmos file variables + u_ind = getindex(vars3d, 'u') !< indices in the state or control var arrays + v_ind = getindex(vars3d, 'v') ! U and V (3D) + tv_ind = getindex(vars3d, 'tv') ! Tv (3D) + q_ind = getindex(vars3d, 'q') ! Q (3D) + oz_ind = getindex(vars3d, 'oz') ! Oz (3D) + cw_ind = getindex(vars3d, 'cw') ! CW (3D) + tsen_ind = getindex(vars3d, 'tsen') !sensible T (3D) + ql_ind = getindex(vars3d, 'ql') ! QL (3D) + qi_ind = getindex(vars3d, 'qi') ! QI (3D) + prse_ind = getindex(vars3d, 'prse') + qr_ind = getindex(vars3d, 'qr') ! QR (3D) + qs_ind = getindex(vars3d, 'qs') ! QS (3D) + qg_ind = getindex(vars3d, 'qg') ! QG (3D) + ps_ind = getindex(vars2d, 'ps') ! Ps (2D) + pst_ind = getindex(vars2d, 'pst') ! Ps tendency (2D) // equivalent of + ! old logical massbal_adjust, if non-zero + sst_ind = getindex(vars2d, 'sst') ! is this really in the atmos file? + + ! for nc gfs io determine if requested variables are in sfc and/or atmos file + atm_file = ( u_ind>0 .or. v_ind>0 .or. tv_ind>0 .or. q_ind>0 .or. sst_ind>0 .or. & + oz_ind>0 .or. cw_ind>0 .or. tsen_ind>0 .or. ql_ind>0 .or. & + qi_ind>0 .or. prse_ind>0 .or. qr_ind>0 .or. qs_ind>0 .or. qg_ind>0 ) + + ! sfc file variables + tmp2m_ind = getindex(vars2d, 't2m') + spfh2m_ind = getindex(vars2d, 'q2m') + soilt1_ind = getindex(vars2d, 'soilt1') + slc1_ind = getindex(vars2d, 'slc1') + soilt2_ind = getindex(vars2d, 'soilt2') + slc2_ind = getindex(vars2d, 'slc2') + soilt3_ind = getindex(vars2d, 'soilt3') + slc3_ind = getindex(vars2d, 'slc3') + soilt4_ind = getindex(vars2d, 'soilt4') + slc4_ind = getindex(vars2d, 'slc4') + + sfc_file = ( tmp2m_ind > 0 .or. spfh2m_ind > 0 .or. soilt1_ind > 0 .or. & + slc1_ind > 0 .or. soilt2_ind > 0 .or. slc2_ind > 0 .or. & + soilt3_ind > 0 .or. slc3_ind > 0 .or. soilt4_ind > 0 .or. & + slc4_ind > 0 ) + + end subroutine set_ncio_file_flags + + logical function checkfield(field,fields,nrec) result(hasfield) use nemsio_module, only: nemsio_charkind integer, intent(in) :: nrec From 299e368736d4236d7ae091a9ccece739f70ab0cf Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Fri, 12 Aug 2022 15:52:54 -0500 Subject: [PATCH 10/12] Tidy up of read/write for sfc files. --- src/enkf/controlvec.f90 | 9 +- src/enkf/gridio_gfs.f90 | 289 +++------------------------------------- src/enkf/statevec.f90 | 7 +- 3 files changed, 24 insertions(+), 281 deletions(-) diff --git a/src/enkf/controlvec.f90 b/src/enkf/controlvec.f90 index 3eebe2ef4a..4a9ad2129b 100644 --- a/src/enkf/controlvec.f90 +++ b/src/enkf/controlvec.f90 @@ -47,8 +47,7 @@ module controlvec mpi_integer,mpi_wtime,mpi_status,mpi_real8 use gridio, only: readgriddata, readgriddata_pnc, writegriddata, writegriddata_pnc, & - writeincrement, writeincrement_pnc, & - writegriddata_2mDA + writeincrement, writeincrement_pnc use gridinfo, only: getgridinfo, gridinfo_cleanup, & npts, vars3d_supported, vars2d_supported use params, only: nlevs, nbackgrounds, fgfileprefixes, reducedgrid, & @@ -380,11 +379,7 @@ subroutine write_control(no_inflate_flag) if (write_fv3_incr) then call writeincrement(0,0,cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin_mean,no_inflate_flag) else - if (global_2mDA) then ! flag for writing sfc analysis - call writegriddata_2mDA(0,0,cvars2d,nc2d,grdin_mean(:,clevels(nc3d)+1:ncdim,:,:),no_inflate_flag) - else call writegriddata(0,0,cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin_mean,no_inflate_flag) - endif end if endif deallocate(grdin_mean) @@ -436,4 +431,4 @@ subroutine controlvec_cleanup() end subroutine controlvec_cleanup end module controlvec - \ No newline at end of file + diff --git a/src/enkf/gridio_gfs.f90 b/src/enkf/gridio_gfs.f90 index 55873a451d..b897bb776e 100644 --- a/src/enkf/gridio_gfs.f90 +++ b/src/enkf/gridio_gfs.f90 @@ -56,8 +56,7 @@ module gridio implicit none private public :: readgriddata, readgriddata_pnc, writegriddata, writegriddata_pnc - public :: writeincrement, writeincrement_pnc, & - writegriddata_2mDA + public :: writeincrement, writeincrement_pnc contains @@ -102,12 +101,21 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, & logical ice logical use_full_hydro integer(i_kind), allocatable, dimension(:) :: mem_pe, lev_pe1, lev_pe2, iocomms - integer(i_kind) :: iope, ionumproc, iolevs, krev + integer(i_kind) :: iope, ionumproc, iolevs, krev, ierr integer(i_kind) :: ncstart(3), nccount(3) ! mpi gatherv things integer(i_kind), allocatable, dimension(:) :: recvcounts, displs real(r_single), dimension(nlons,nlats,nlevs) :: ug3d_0, vg3d_0 + logical :: read_sfc_file, read_atm_file + + call set_ncio_file_flags(vars3d, n3d, vars2d, n2d, read_sfc_file, read_atm_file) + + if (read_sfc_file) then ! flag for parallel rad + 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)) @@ -1030,7 +1038,7 @@ subroutine readgriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, call stop2(23) endif if ( reducedgrid ) then - write(6,*) "reducedgrid=T not valid for global_2mDA=T" + write(6,*) "reducedgrid=T not valid for writing sfc files" call stop2(22) endif @@ -1969,270 +1977,6 @@ end subroutine copyfromgrdin end subroutine writegriddata_pnc - subroutine writegriddata_2mDA(nanal1,nanal2,vars2d,n2d,grdin_2d,no_inflate_flag) - use netcdf - use module_fv3gfs_ncio, only: Dataset, Variable, Dimension, open_dataset,& - read_attribute, close_dataset, get_dim, read_vardata,& - create_dataset, get_idate_from_time_units, & - get_time_units_from_idate, write_vardata, & - write_attribute, quantize_data, has_var, has_attr - use params, only: nbackgrounds,anlsfcfileprefixes,fgsfcfileprefixes,& - nccompress,write_ensmean - implicit none - - integer, intent(in) :: nanal1,nanal2 - character(len=max_varname_length), dimension(n2d), intent(in) :: vars2d - integer, intent(in) :: n2d - ! only 2d segment of grdin should be input here - real(r_single), dimension(npts,n2d,nbackgrounds,nanal2-nanal1+1), intent(inout) :: grdin_2d - logical, intent(in) :: no_inflate_flag - character(len=500):: filenamein, filenameout - real(r_kind), dimension(nlons*nlats) :: ug,vg - real(r_kind), allocatable, dimension(:) :: values_1d - real(r_single), allocatable, dimension(:,:) :: values_2d, values_2d_tmp - integer iadate(4),idate(4),nfhour,idat(7),iret,jdat(6),nbits - integer,dimension(8):: ida,jda - real(r_double),dimension(5):: fha - real(r_kind) fhour - type(Dataset) :: dsfg, dsanl - character(len=3) charnanal - character(len=nf90_max_name) :: time_units - - real(r_single) compress_err - - integer(i_kind) :: tmp2m_ind, spfh2m_ind, soilt1_ind, soilt2_ind, soilt3_ind, & - soilt4_ind,slc1_ind, slc2_ind, slc3_ind, slc4_ind - - - integer nb,ne,nanal - - logical :: nocompress - - nocompress = .true. - if (nccompress) nocompress = .false. - - ne = 0 - ensmemloop: do nanal=nanal1,nanal2 - ne = ne + 1 - write(charnanal,'(i3.3)') nanal - backgroundloop: do nb=1,nbackgrounds - - if (nanal == 0 .and. write_ensmean) then - filenamein = trim(adjustl(datapath))//trim(adjustl(fgsfcfileprefixes(nb)))//"ensmean" - filenameout = trim(adjustl(datapath))//trim(adjustl(anlsfcfileprefixes(nb)))//"ensmean" - else - if(no_inflate_flag) then - filenameout = trim(adjustl(datapath))//trim(adjustl(anlsfcfileprefixes(nb)))//"nimem"//charnanal - else - filenameout = trim(adjustl(datapath))//trim(adjustl(anlsfcfileprefixes(nb)))//"mem"//charnanal - end if - filenamein = trim(adjustl(datapath))//trim(adjustl(fgsfcfileprefixes(nb)))//"mem"//charnanal - endif - - dsfg = open_dataset(filenamein) - jdat = get_idate_from_time_units(dsfg) - idat(4) = jdat(1) ! yr - idat(2) = jdat(2) ! mon - idat(3) = jdat(3) ! day - idat(1) = jdat(4) ! hr - call read_vardata(dsfg,'time',values_1d,errcode=iret) - if (iret /= 0) then - print *,'error reading time' - call stop2(29) - endif - nfhour = int(values_1d(1)) - - tmp2m_ind = getindex(vars2d, 't2m') !< indices in the state or control var arrays - spfh2m_ind = getindex(vars2d, 'q2m') - soilt1_ind = getindex(vars2d, 'soilt1') - slc1_ind = getindex(vars2d, 'slc1') - soilt2_ind = getindex(vars2d, 'soilt2') - slc2_ind = getindex(vars2d, 'slc2') - soilt3_ind = getindex(vars2d, 'soilt3') - slc3_ind = getindex(vars2d, 'slc3') - soilt4_ind = getindex(vars2d, 'soilt4') - slc4_ind = getindex(vars2d, 'slc4') - - idate(3)=idat(3) !day - idate(2)=idat(2) !mon - idate(4)=idat(4) !yr - idate(1)=idat(1) !hr - fhour = nfhour - fha=zero; ida=0; jda=0 - fha(2)=fhour ! relative time interval in hours - ida(1)=idate(4) ! year - ida(2)=idate(2) ! month - ida(3)=idate(3) ! day - ida(4)=0 ! time zone - ida(5)=idate(1) ! hour - call w3movdat(fha,ida,jda) -! -! INPUT VARIABLES: -! RINC REAL (5) NCEP RELATIVE TIME INTERVAL -! (DAYS, HOURS, MINUTES, SECONDS, MILLISECONDS) -! IDAT INTEGER (8) NCEP ABSOLUTE DATE AND TIME -! (YEAR, MONTH, DAY, TIME ZONE, -! HOUR, MINUTE, SECOND, MILLISECOND) -! -! OUTPUT VARIABLES: -! JDAT INTEGER (8) NCEP ABSOLUTE DATE AND TIME -! (YEAR, MONTH, DAY, TIME ZONE, -! HOUR, MINUTE, SECOND, MILLISECOND) -! (JDAT IS LATER THAN IDAT IF TIME INTERVAL IS POSITIVE.) - iadate(1)=jda(5) ! hour - iadate(2)=jda(2) ! mon - iadate(3)=jda(3) ! day - iadate(4)=jda(1) ! year - if (nproc .eq. 0) then - print *,'idate = ',idate - print *,'iadate = ',iadate - end if - - dsanl = create_dataset(filenameout, dsfg, copy_vardata=.true., nocompress=nocompress, errcode=iret) - if (iret /= 0) then - print *,'error creating netcdf file' - call stop2(29) - endif - deallocate(values_1d) - allocate(values_1d(1)) - values_1d(1)=zero - call write_vardata(dsanl,'time',values_1d,errcode=iret) - if (iret /= 0) then - print *,'error writing time' - call stop2(29) - endif - jdat(1) = iadate(4) - jdat(2) = iadate(2) - jdat(3) = iadate(3) - jdat(4) = iadate(1) - jdat(5) = jda(6); jdat(6) = jda(7) - time_units = get_time_units_from_idate(jdat) - call write_attribute(dsanl,'units',time_units,'time',errcode=iret) - if (iret /= 0) then - print *,'error writing time units attribute' - call stop2(29) - endif - - if (tmp2m_ind > 0) then - call read_vardata(dsfg, 'tmp2m', values_2d, errcode=iret) - if (iret /= 0) then - print *,'error reading tmp2m' - call stop2(29) - endif - vg = reshape(values_2d,(/nlons*nlats/)) - call copyfromgrdin(grdin_2d(:,tmp2m_ind,nb,ne),ug) - ! add increment to background. - values_2d = reshape(ug+vg,(/nlons,nlats/)) - ! quantize data if lossy compression desired. - if (has_attr(dsfg, 'nbits', 'tmp2m') .and. .not. nocompress) then - call read_attribute(dsfg, 'nbits', nbits, 'tmp2m') - if (.not. allocated(values_2d_tmp)) allocate(values_2d_tmp(nlons,nlats)) - values_2d_tmp = values_2d - call quantize_data(values_2d_tmp, values_2d, nbits, compress_err) - call write_attribute(dsanl,& - 'max_abs_compression_error',compress_err,'tmp2m',errcode=iret) - if (iret /= 0) then - print *,'error writing tmp2m attribute' - call stop2(29) - endif - endif - call write_vardata(dsanl,'tmp2m',values_2d,errcode=iret) - if (iret /= 0) then - print *,'error writing tmp2m' - call stop2(29) - endif - endif - if (spfh2m_ind > 0) then - call read_vardata(dsfg, 'spfh2m', values_2d, errcode=iret) - if (iret /= 0) then - print *,'error reading spfh2m' - call stop2(29) - endif - vg = reshape(values_2d,(/nlons*nlats/)) - call copyfromgrdin(grdin_2d(:,spfh2m_ind,nb,ne),ug) - ! add increment to background. - values_2d = reshape(ug+vg,(/nlons,nlats/)) - ! quantize data if lossy compression desired. - if (has_attr(dsfg, 'nbits', 'spfh2m') .and. .not. nocompress) then - call read_attribute(dsfg, 'nbits', nbits, 'spfh2m') - if (.not. allocated(values_2d_tmp)) allocate(values_2d_tmp(nlons,nlats)) - values_2d_tmp = values_2d - call quantize_data(values_2d_tmp, values_2d, nbits, compress_err) - call write_attribute(dsanl,& - 'max_abs_compression_error',compress_err,'spfh2m',errcode=iret) - if (iret /= 0) then - print *,'error writing spfh2m attribute' - call stop2(29) - endif - endif - call write_vardata(dsanl,'spfh2m',values_2d,errcode=iret) - if (iret /= 0) then - print *,'error writing spfh2m' - call stop2(29) - endif - endif - if (soilt1_ind > 0) then - call read_vardata(dsfg, 'soilt1', values_2d, errcode=iret) - if (iret /= 0) then - print *,'error reading soilt1' - call stop2(29) - endif - vg = reshape(values_2d,(/nlons*nlats/)) - call copyfromgrdin(grdin_2d(:,soilt1_ind,nb,ne),ug) - ! add increment to background. - values_2d = reshape(ug+vg,(/nlons,nlats/)) - ! quantize data if lossy compression desired. - if (has_attr(dsfg, 'nbits', 'soilt1') .and. .not. nocompress) then - call read_attribute(dsfg, 'nbits', nbits, 'soilt1') - if (.not. allocated(values_2d_tmp)) allocate(values_2d_tmp(nlons,nlats)) - values_2d_tmp = values_2d - call quantize_data(values_2d_tmp, values_2d, nbits, compress_err) - call write_attribute(dsanl,& - 'max_abs_compression_error',compress_err,'soilt1',errcode=iret) - if (iret /= 0) then - print *,'error writing soilt1 attribute' - call stop2(29) - endif - endif - call write_vardata(dsanl,'soilt1',values_2d,errcode=iret) - if (iret /= 0) then - print *,'error writing soilt1' - call stop2(29) - endif - endif - - if (allocated(values_2d)) deallocate(values_2d) - if (allocated(values_2d_tmp)) deallocate(values_2d_tmp) - if (allocated(values_1d)) deallocate(values_1d) - - call close_dataset(dsfg,errcode=iret) - if (iret/=0) then - write(6,*)'gridio/writegriddata: gfs model: problem closing netcdf fg dataset, iret=',iret - call stop2(23) - endif - call close_dataset(dsanl,errcode=iret) - if (iret/=0) then - write(6,*)'gridio/writegriddata: gfs model: problem closing netcdf anal dataset, iret=',iret - call stop2(23) - endif - - end do backgroundloop ! loop over backgrounds to write out - end do ensmemloop ! loop over ens members to write out - - contains -! copying to grdin (calling regtoreduced if reduced grid) - subroutine copyfromgrdin(grdin, field) - implicit none - - real(r_single), dimension(:), intent(in) :: grdin - real(r_kind), dimension(:), intent(inout) :: field - - field = grdin - - end subroutine copyfromgrdin - - end subroutine writegriddata_2mDA - subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,no_inflate_flag) use netcdf use sigio_module, only: sigio_head, sigio_data, sigio_sclose, sigio_sropen, & @@ -2301,6 +2045,15 @@ subroutine writegriddata(nanal1,nanal2,vars3d,vars2d,n3d,n2d,levels,ndim,grdin,n integer k,krev,nt,ierr,iunitsig,nb,i,ne,nanal logical :: nocompress + logical :: write_sfc_file, write_atm_file + + call set_ncio_file_flags(vars3d, n3d, vars2d, n2d, write_sfc_file, write_atm_file) + + if (write_sfc_file .and. nproc==0 ) then + ! adding the sfc increments requires adjusting several other variables. This is done is a separate + ! program. + write(6,*)'gridio/writegriddata: not coded to write sfc analysis, use separate add_incr program instead' + endif nocompress = .true. if (nccompress) nocompress = .false. diff --git a/src/enkf/statevec.f90 b/src/enkf/statevec.f90 index 4469550e86..1d18a01316 100644 --- a/src/enkf/statevec.f90 +++ b/src/enkf/statevec.f90 @@ -46,7 +46,7 @@ module statevec npts, vars3d_supported, vars2d_supported use params, only: nlevs,nstatefields,nanals,statefileprefixes,& ntasks_io,nanals_per_iotask,nanal1,nanal2, & - statesfcfileprefixes, paranc, global_2mDA + statesfcfileprefixes, paranc use kinds, only: r_kind, i_kind, r_double, r_single use mpeu_util, only: gettablesize, gettable, getindex use constants, only : max_varname_length @@ -188,11 +188,6 @@ subroutine read_state() allocate(state_d(npts,nsdim,nstatefields,nanals_per_iotask)) allocate(qsat(npts,nlevs,nstatefields,nanals_per_iotask)) if (paranc) then - if (global_2mDA) then ! flag for parallel rad - print *,'paranc not supported for global_2mDA' - call mpi_barrier(mpi_comm_world,ierr) - call mpi_finalize(ierr) - endif call readgriddata_pnc(svars3d,svars2d,ns3d,ns2d,slevels,nsdim,nstatefields, & statefileprefixes,statesfcfileprefixes,.false.,state_d,qsat) end if From 7046544a8b69fb963f4a55cd3297a01f061ebc05 Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Fri, 12 Aug 2022 16:15:20 -0500 Subject: [PATCH 11/12] Removed global_2mDA input flag. --- src/enkf/controlvec.f90 | 12 +----------- src/enkf/mpi_readobs.f90 | 11 +++-------- src/enkf/params.f90 | 4 +--- src/enkf/readconvobs.f90 | 19 +++++-------------- 4 files changed, 10 insertions(+), 36 deletions(-) diff --git a/src/enkf/controlvec.f90 b/src/enkf/controlvec.f90 index 4a9ad2129b..b6cd0d9eaa 100644 --- a/src/enkf/controlvec.f90 +++ b/src/enkf/controlvec.f90 @@ -52,7 +52,7 @@ module controlvec npts, vars3d_supported, vars2d_supported use params, only: nlevs, nbackgrounds, fgfileprefixes, reducedgrid, & nanals, pseudo_rh, use_qsatensmean, nlons, nlats,& - nanals_per_iotask, ntasks_io, nanal1, nanal2, global_2mDA, & + nanals_per_iotask, ntasks_io, nanal1, nanal2,& fgsfcfileprefixes, paranc, write_fv3_incr, write_ensmean use kinds, only: r_kind, i_kind, r_double, r_single use mpeu_util, only: gettablesize, gettable, getindex @@ -215,11 +215,6 @@ subroutine read_control() q_ind = getindex(cvars3d, 'q') if (q_ind > 0) allocate(qsat(npts,nlevs,nbackgrounds,nanals_per_iotask)) if (paranc) then - if (global_2mDA) then ! flag for parallel read support - print *,'paranc not supported for global_2mDA' - call mpi_barrier(mpi_comm_world,ierr) - call mpi_finalize(ierr) - endif if (nproc == 0) t1 = mpi_wtime() call readgriddata_pnc(cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,nbackgrounds, & fgfileprefixes,fgsfcfileprefixes,reducedgrid,grdin,qsat) @@ -391,11 +386,6 @@ subroutine write_control(no_inflate_flag) end if ! io task if (paranc) then - if (global_2mDA) then ! flag for parralell write - print *,'paranc not supported for global_2mDA' - call mpi_barrier(mpi_comm_world,ierr) - call mpi_finalize(ierr) - endif if (write_fv3_incr) then call writeincrement_pnc(cvars3d,cvars2d,nc3d,nc2d,clevels,ncdim,grdin,no_inflate_flag) else diff --git a/src/enkf/mpi_readobs.f90 b/src/enkf/mpi_readobs.f90 index db9c276f7d..d789dfcfff 100644 --- a/src/enkf/mpi_readobs.f90 +++ b/src/enkf/mpi_readobs.f90 @@ -33,7 +33,7 @@ module mpi_readobs !$$$ use kinds, only: r_double,i_kind,r_kind,r_single,num_bytes_for_r_single -use params, only: ntasks_io, nanals_per_iotask, nanal1, nanal2, global_2mDA +use params, only: ntasks_io, nanals_per_iotask, nanal1, nanal2 use radinfo, only: npred use readconvobs use readsatobs @@ -89,13 +89,8 @@ subroutine mpi_getobs(obspath, datestring, nobs_conv, nobs_oz, nobs_sat, nobs_to ! get total number of conventional and sat obs for ensmean. id = 'ensmean' if(nproc == 0)call get_num_convobs(obspath,datestring,nobs_conv,nobs_convdiag,id) - if (global_2mDA) then ! flag for not reading other conventional obs - if(nproc == iozproc) nobs_oz=0; nobs_ozdiag=0 - if(nproc == isatproc) nobs_sat=0; nobs_satdiag=0 - else - if(nproc == iozproc)call get_num_ozobs(obspath,datestring,nobs_oz,nobs_ozdiag,id) - if(nproc == isatproc)call get_num_satobs(obspath,datestring,nobs_sat,nobs_satdiag,id) - endif + if(nproc == iozproc)call get_num_ozobs(obspath,datestring,nobs_oz,nobs_ozdiag,id) + if(nproc == isatproc)call get_num_satobs(obspath,datestring,nobs_sat,nobs_satdiag,id) call mpi_bcast(nobs_conv,1,mpi_integer,0,mpi_comm_world,ierr) call mpi_bcast(nobs_convdiag,1,mpi_integer,0,mpi_comm_world,ierr) call mpi_bcast(nobs_oz,1,mpi_integer,iozproc,mpi_comm_world,ierr) diff --git a/src/enkf/params.f90 b/src/enkf/params.f90 index 6416ceca35..6ac3d2580b 100644 --- a/src/enkf/params.f90 +++ b/src/enkf/params.f90 @@ -157,8 +157,6 @@ module params ! matrix are read from a file called 'vlocal_eig.dat' ! (created by an external python utility). logical,public :: modelspace_vloc=.false. -! if global_2mDA true, update land states only using 2m obs in diag_t, diag_q -logical, public :: global_2mDA=.false. ! use correlated obs errors ! (implies letkf_flag=T, modelspace_vloc=T and lobsdiag_forenkf=T) ! if T, extra fields read from diag file and innovation stats @@ -269,7 +267,7 @@ module params lnsigcutoffsatnh,lnsigcutoffsattr,lnsigcutoffsatsh,& lnsigcutoffpsnh,lnsigcutoffpstr,lnsigcutoffpssh,& fgfileprefixes,fgsfcfileprefixes,anlfileprefixes, & - anlsfcfileprefixes,incfileprefixes,incsfcfileprefixes,global_2mDA,& + anlsfcfileprefixes,incfileprefixes,incsfcfileprefixes,& statefileprefixes,statesfcfileprefixes, & covl_minfact,covl_efold,lupd_obspace_serial,letkf_novlocal,& analpertwtnh,analpertwtsh,analpertwttr,sprd_tol,& diff --git a/src/enkf/readconvobs.f90 b/src/enkf/readconvobs.f90 index 6b47d58ce4..5901bfdc1d 100644 --- a/src/enkf/readconvobs.f90 +++ b/src/enkf/readconvobs.f90 @@ -32,7 +32,7 @@ module readconvobs use kinds, only: r_kind,i_kind,r_single,r_double use constants, only: one,zero,deg2rad -use params, only: npefiles, netcdf_diag, modelspace_vloc, global_2mDA, & +use params, only: npefiles, netcdf_diag, modelspace_vloc,& l_use_enkf_directZDA,diagprefix implicit none @@ -281,8 +281,6 @@ subroutine get_num_convobs_nc(obspath,datestring,num_obs_tot,num_obs_totdiag,id) obtypeloop: do itype=1, nobtype obtype = obtypes(itype) - ! only read t and q obs for global_2mDA - if (global_2mDA .and. (obtype .ne. ' t' .and. obtype .ne. ' q')) cycle obtypeloop ! flag for nor reading other convention obs type peloop: do ipe=0,npefiles write(pe_name,'(i4.4)') ipe @@ -334,14 +332,11 @@ subroutine get_num_convobs_nc(obspath,datestring,num_obs_tot,num_obs_totdiag,id) call nc_diag_read_close(obsfile) - num_obs_totdiag = num_obs_totdiag + nobs_curr do i = 1, nobs_curr - ! for global_2mDA skip if not 2m (surface) ob ityp = Observation_Type(i) sfctype=(ityp>179.and.ityp<190).or.(ityp>=192.and.ityp<=199) - if (global_2mDA .and. .not. sfctype) cycle ! flag for not reading other conventional obs types errorlimit2=errorlimit2_obs @@ -350,7 +345,7 @@ subroutine get_num_convobs_nc(obspath,datestring,num_obs_tot,num_obs_totdiag,id) endif ! for q, normalize by qsatges - if (obtype == ' q' .and. .not. global_2mDA) then ! flag for normalzing q by ens mean (if q read in) + if (obtype == ' q' .and. .not. sfctype) then ! not normalizing for sfc q (should we?) obmax = abs(Observation(i) / Forecast_Saturation_Spec_Hum(i)) error = Errinv_Final(i) * Forecast_Saturation_Spec_Hum(i) else @@ -552,8 +547,6 @@ subroutine get_convobs_data_nc(obspath, datestring, nobs_max, nobs_maxdiag, & obtypeloop: do itype=1, nobtype obtype = obtypes(itype) - ! only read t and q obs for global_2mDA - if (global_2mDA .and. (obtype .ne. ' t' .and. obtype .ne. ' q')) cycle obtypeloop ! flag for not reading other obs types peloop: do ipe=0,npefiles write(pe_name,'(i4.4)') ipe @@ -675,10 +668,8 @@ subroutine get_convobs_data_nc(obspath, datestring, nobs_max, nobs_maxdiag, & errorlimit2=errorlimit2_obs do i = 1, nobs - ! for global_2mDA skip if not 2m (surface) ob ityp = Observation_Type(i) sfctype=(ityp>179.and.ityp<190).or.(ityp>=192.and.ityp<=199) - if (global_2mDA .and. .not. sfctype) cycle ! flag for not reading other conventional obs types nobdiag = nobdiag + 1 ! special handling for error limits for GPS bend angle if (obtype == 'gps') then @@ -686,7 +677,7 @@ subroutine get_convobs_data_nc(obspath, datestring, nobs_max, nobs_maxdiag, & endif ! for q, normalize by qsatges - if (.not. global_2mDA .and. obtype == ' q') then ! flag for normazling q + if (.not. sfctype .and. obtype == ' q') then ! not normalizing for sfc q (should we?) obmax = abs(real(Observation(i),r_single) / real(Forecast_Saturation_Spec_Hum(i),r_single)) errororig = real(Errinv_Input(i),r_single) * real(Forecast_Saturation_Spec_Hum(i),r_single) error = real(Errinv_Final(i),r_single) * real(Forecast_Saturation_Spec_Hum(i),r_single) @@ -795,13 +786,13 @@ subroutine get_convobs_data_nc(obspath, datestring, nobs_max, nobs_maxdiag, & endif ! normalize q by qsatges - if (obtype == ' q' .and. .not. global_2mDA) then ! flag for normalizing q + if (obtype == ' q' .and. .not. sfctype ) then hx(nob) = hx(nob) / Forecast_Saturation_Spec_Hum(i) endif endif ! normalize q by qsatges - if (obtype == ' q' .and. .not. global_2mDA) then ! flag for normalizing q + if (obtype == ' q' .and. .not. sfctype ) then x_obs(nob) = x_obs(nob) /Forecast_Saturation_Spec_Hum(i) hx_mean(nob) = hx_mean(nob) /Forecast_Saturation_Spec_Hum(i) hx_mean_nobc(nob) = hx_mean_nobc(nob) /Forecast_Saturation_Spec_Hum(i) From cf6849796447e91812ad94c319e1ce490d40dc01 Mon Sep 17 00:00:00 2001 From: ClaraDraper-NOAA Date: Fri, 12 Aug 2022 17:34:53 -0500 Subject: [PATCH 12/12] Updates ot gsi: 1. change variable from use_2m_obs to hofx_2m_sfcfile 2. removed unnecessary use of above variable 3. reverted to original buddycheck routine (doesn't work for hofx_2m_sfcfile option anyway). 4. removed soil temp and moisture nuding with hofx_2m_sfcfile option. --- src/gsi/buddycheck_mod.f90 | 39 --------------------------- src/gsi/general_read_gfsatm.f90 | 22 +++------------ src/gsi/gsimod.F90 | 4 +-- src/gsi/jfunc.f90 | 7 ++--- src/gsi/netcdfgfs_io.f90 | 12 ++++----- src/gsi/read_prepbufr.f90 | 10 +++---- src/gsi/setupt.f90 | 47 +++++++++++++-------------------- src/gsi/update_guess.f90 | 6 ++--- 8 files changed, 43 insertions(+), 104 deletions(-) diff --git a/src/gsi/buddycheck_mod.f90 b/src/gsi/buddycheck_mod.f90 index d23198a3ed..ae4a29cc63 100644 --- a/src/gsi/buddycheck_mod.f90 +++ b/src/gsi/buddycheck_mod.f90 @@ -60,7 +60,6 @@ subroutine buddy_check_t(is,data,luse,mype,nele,nobs,muse,buddyuse) use aircraftinfo, only: aircraft_t_bc_pof,aircraft_t_bc, & aircraft_t_bc_ext use convinfo, only: ictype - use jfunc, only: use_2m_obs implicit none ! !INPUT PARAMETERS: @@ -154,15 +153,11 @@ subroutine buddy_check_t(is,data,luse,mype,nele,nobs,muse,buddyuse) real(r_kind),allocatable,dimension(:,:,: ) :: ges_ps real(r_kind),allocatable,dimension(:,:,: ) :: ges_z - real(r_kind),allocatable,dimension(:,:,: ) :: ges_t2 real(r_kind),allocatable,dimension(:,:,:,:) :: ges_u real(r_kind),allocatable,dimension(:,:,:,:) :: ges_v real(r_kind),allocatable,dimension(:,:,:,:) :: ges_tv real(r_kind),allocatable,dimension(:,:,:,:) :: ges_q - real(r_kind) :: delta_z - real(r_kind), parameter :: T_lapse = -0.0045 ! must match setupt - ! Check to see if required guess fields are available call check_vars_(proceed) if(.not.proceed) return ! not all vars available, simply return @@ -281,19 +276,6 @@ subroutine buddy_check_t(is,data,luse,mype,nele,nobs,muse,buddyuse) f10ges,u10ges,v10ges, t2ges, q2ges, regime, iqtflg) tges = t2ges - elseif(sfctype .and. use_2m_obs) then -! SCENARIO 2: obs is sfctype, and use_2m_obs scheme is on. -! 2m forecast has been read from the sfc guess files -! note: this block should match that in setupt. - - call tintrp2a11(ges_t2,tges,dlat,dlon,dtime,hrdifsig,& - mype,nfldsig) - -! correct obs to model terrain height - equivalent to gsd_terrain_match - - delta_z = data(izz,i) - data(istnelv,i) - tob = tob + delta_z*T_lapse - else if(iqtflg)then ! Interpolate guess tv to observation location and time @@ -467,26 +449,6 @@ subroutine init_vars_ write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus call stop2(999) endif - if( use_2m_obs) then -! get t2m ... - varname='t2m' - call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus) - if (istatus==0) then - if(allocated(ges_t2))then - write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc ' - call stop2(999) - endif - allocate(ges_t2(size(rank2,1),size(rank2,2),nfldsig)) - ges_t2(:,:,1)=rank2 - do ifld=2,nfldsig - call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank2,istatus) - ges_t2(:,:,ifld)=rank2 - enddo - else - write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus - call stop2(999) - endif - endif else write(6,*) trim(myname), ': inconsistent vector sizes (nfldsig,size(metguess_bundle) ',& nfldsig,size(gsi_metguess_bundle) @@ -500,7 +462,6 @@ subroutine final_vars_ if(allocated(ges_v )) deallocate(ges_v ) if(allocated(ges_u )) deallocate(ges_u ) if(allocated(ges_ps)) deallocate(ges_ps) - if(allocated(ges_t2)) deallocate(ges_t2) end subroutine final_vars_ end subroutine buddy_check_t diff --git a/src/gsi/general_read_gfsatm.f90 b/src/gsi/general_read_gfsatm.f90 index 576e502f76..d7059dc84a 100755 --- a/src/gsi/general_read_gfsatm.f90 +++ b/src/gsi/general_read_gfsatm.f90 @@ -1687,7 +1687,7 @@ subroutine general_read_gfsatm_nems(grd,sp_a,filename,uvflag,vordivflag,zflag, & end subroutine general_read_gfsatm_nems -subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag,use_2m_obs, & +subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & gfs_bundle,init_head,iret_read) !$$$ subprogram documentation block ! . . . . @@ -1748,7 +1748,6 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag,use_ type(spec_vars) ,intent(in ) :: sp_a character(*) ,intent(in ) :: filename logical ,intent(in ) :: uvflag,zflag,vordivflag,init_head - logical ,intent(in ) :: use_2m_obs integer(i_kind) ,intent( out) :: iret_read type(gsi_bundle) ,intent(inout) :: gfs_bundle @@ -1810,21 +1809,11 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag,use_ read_z = .false. if ( mype == 0 ) write(6,* ) & trim(my_name), ': reading 2m variables from ', trim(filename) - if (.not. use_2m_obs) then - write(6,*) 'error, requesting sfc file read, but not use_2m_obs. Check options' - call stop2(101) - endif else read_2m = .false. - !if (use_2m_obs) then - ! read_z =.true. ! over-ride, as this is the only var we want. - ! if ( mype == 0 ) write(6,* ) & - ! trim(my_name), ': reading z only from ', trim(filename) - !else - read_z = zflag - if ( mype == 0 ) write(6,* ) & - trim(my_name), ': reading atmos variables from ', trim(filename) - !endif + read_z = zflag + if ( mype == 0 ) write(6,* ) & + trim(my_name), ': reading atmos variables from ', trim(filename) endif @@ -2052,15 +2041,12 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag,use_ call general_fill_ns(grd,grid,work) endif endif - ! CSD for 2 m case, icm = 2. this doesn't get triggered. - !if ( use_2m_obs .or. ( (.not. use_2m_obs ) .and. ( icount == icm)) ) then if ( read_2m .or. ( (.not. read_2m ) .and. ( icount == icm)) ) then call general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, & icount,iflag,ilev,work,uvflag,vordivflag) endif endif - !if ((.not. read_2m) .and. (.not. use_2m_obs)) then if (.not. read_2m) then icount=icount+1 diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index bf75561fb0..d03685677a 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -98,7 +98,7 @@ module gsimod factv,factl,factp,factg,factw10m,facthowv,factcldch,niter,niter_no_qc,biascor,& init_jfunc,qoption,cwoption,switch_on_derivatives,tendsflag,jiterstart,jiterend,R_option,& bcoption,diurnalbc,print_diag_pcg,tsensible,diag_precon,step_start,pseudo_q2,& - clip_supersaturation,cnvw_option,use_2m_obs + clip_supersaturation,cnvw_option,hofx_2m_sfcfile use state_vectors, only: init_anasv,final_anasv use control_vectors, only: init_anacv,final_anacv,nrf,nvars,nrf_3d,cvars3d,cvars2d,& nrf_var,lcalc_gfdl_cfrac,incvars_to_zero,incvars_zero_strat,incvars_efold @@ -1020,7 +1020,7 @@ module gsimod ! l_foreaft_thin - separate TDR fore/aft scan for thinning namelist/obs_input/dmesh,time_window_max,time_window_rad, & - ext_sonde,l_foreaft_thin,use_2m_obs + ext_sonde,l_foreaft_thin,hofx_2m_sfcfile ! SINGLEOB_TEST (one observation test case setup): ! maginnov - magnitude of innovation for one ob diff --git a/src/gsi/jfunc.f90 b/src/gsi/jfunc.f90 index bafa59b85a..1f5924102e 100644 --- a/src/gsi/jfunc.f90 +++ b/src/gsi/jfunc.f90 @@ -136,12 +136,12 @@ module jfunc public :: pseudo_q2 public :: varq public :: cnvw_option - public :: use_2m_obs + public :: hofx_2m_sfcfile logical first,last,switch_on_derivatives,tendsflag,print_diag_pcg,tsensible,diag_precon logical clip_supersaturation,R_option logical pseudo_q2,limitqobs - logical use_2m_obs + logical hofx_2m_sfcfile logical cnvw_option integer(i_kind) iout_iter,miter,iguess,nclen,qoption,cwoption integer(i_kind) jiter,jiterstart,jiterend,iter @@ -251,7 +251,8 @@ subroutine init_jfunc ! option for including convective clouds in the all-sky assimilation cnvw_option=.false. - use_2m_obs=.false. +! to calculate hofx for 2m obs, interpolate from gfs sfc file. + hofx_2m_sfcfile=.false. return end subroutine init_jfunc diff --git a/src/gsi/netcdfgfs_io.f90 b/src/gsi/netcdfgfs_io.f90 index 12a8b8f779..f687ba779a 100644 --- a/src/gsi/netcdfgfs_io.f90 +++ b/src/gsi/netcdfgfs_io.f90 @@ -128,7 +128,7 @@ subroutine read_ use general_sub2grid_mod, only: sub2grid_info,general_sub2grid_create_info,general_sub2grid_destroy_info use mpimod, only: npe,mype use cloud_efr_mod, only: cloud_calc_gfs,set_cloud_lower_bound - use jfunc, only: use_2m_obs + use jfunc, only: hofx_2m_sfcfile implicit none character(len=*),parameter::myname_=myname//'*read_' @@ -178,7 +178,7 @@ subroutine read_ ! Allocate bundle used for reading members call gsi_gridcreate(atm_grid,lat2,lon2,nsig) - if (use_2m_obs) then + if (hofx_2m_sfcfile) then call gsi_bundlecreate(atm_bundle,atm_grid,'aux-atm-read',istatus,names2d=vars2d_with2m,names3d=vars3d) else call gsi_bundlecreate(atm_bundle,atm_grid,'aux-atm-read',istatus,names2d=vars2d,names3d=vars3d) @@ -191,20 +191,20 @@ subroutine read_ do it=1,nfldsig ! Read background fields into bundle - if (use_2m_obs) then ! current read_sfc routines called from differebt part of + if (hofx_2m_sfcfile) then ! current read_sfc routines called from different part of ! of the code, can't easily read into the met-bundle ! wrote a new routine here write(filename,'(''sfcf'',i2.2)') ifilesig(it) call general_read_gfsatm_nc(grd_t,sp_a,filename,.true.,.true.,.true.,& - use_2m_obs,atm_bundle,.true.,istatus) + atm_bundle,.true.,istatus) write(filename,'(''sigf'',i2.2)') ifilesig(it) call general_read_gfsatm_nc(grd_t,sp_a,filename,.true.,.true.,.true.,& - use_2m_obs,atm_bundle,.true.,istatus) + atm_bundle,.true.,istatus) else write(filename,'(''sigf'',i2.2)') ifilesig(it) call general_read_gfsatm_nc(grd_t,sp_a,filename,.true.,.true.,.true.,& - use_2m_obs,atm_bundle,.true.,istatus) + atm_bundle,.true.,istatus) endif inithead=.false. diff --git a/src/gsi/read_prepbufr.f90 b/src/gsi/read_prepbufr.f90 index 190cbf7ea0..24572877a7 100644 --- a/src/gsi/read_prepbufr.f90 +++ b/src/gsi/read_prepbufr.f90 @@ -211,7 +211,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& use hilbertcurve,only: init_hilbertcurve, accum_hilbertcurve, & apply_hilbertcurve,destroy_hilbertcurve use ndfdgrids,only: init_ndfdgrid,destroy_ndfdgrid,relocsfcob,adjust_error - use jfunc, only: tsensible, use_2m_obs + use jfunc, only: tsensible, hofx_2m_sfcfile use deter_sfc_mod, only: deter_sfc_type,deter_sfc2 use gsi_nstcouplermod, only: nst_gsi,nstinfo use gsi_nstcouplermod, only: gsi_nstcoupler_deter @@ -1617,7 +1617,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& ! versus sensible temperature if(tob) then ! use tvirtual if tsensible flag not set, and not in either 2Dregional or global_2m DA mode - if ( (.not. tsensible) .and. .not. (twodvar_regional .or. use_2m_obs) ) then + if ( (.not. tsensible) .and. .not. (twodvar_regional .or. hofx_2m_sfcfile) ) then do k=1,levs tvflg(k)=one ! initialize as sensible @@ -1917,7 +1917,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& ! CSD - temporary hack ( move to prepbufr pre-processing) ! Over-ride QM=9 for sfc obs - if (sfctype .and. use_2m_obs ) then + if (sfctype .and. hofx_2m_sfcfile ) then if (tob .and. qm == 9 ) qm = 2 ! 2=not checked if (qob .and. qm == 9 ) qm = 2 ! 2=not checked endif @@ -2079,7 +2079,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& if (aircraftobs .and. aircraft_t_bc .and. acft_profl_file) then call errormod_aircraft(pqm,tqm,levs,plevs,errout,k,presl,dpres,nsig,lim_qm,hdr3) else - if (.not. (sfctype .and. use_2m_obs)) then + if (.not. (sfctype .and. hofx_2m_sfcfile)) then call errormod(pqm,tqm,levs,plevs,errout,k,presl,dpres,nsig,lim_qm) endif end if @@ -2336,7 +2336,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,& end if qobcon=obsdat(2,k)*convert tdry=r999 - if (sfctype .and. use_2m_obs) then + if (sfctype .and. hofx_2m_sfcfile) then if (tqm(k) < 10) tdry=(obsdat(3,k)+t0c)/(one+fv*qobcon) else if (tqm(k)0) hofx_2m_sfcfile=.false. + thead => obsLL(:) save_jacobian = conv_diagsave .and. jiter==jiterstart .and. lobsdiag_forenkf @@ -437,7 +440,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav itype=ictype(ikx) sfctype=(itype>179.and.itype<190).or.(itype>=192.and.itype<=199) do l=k+1,nobs - if (twodvar_regional .or. (use_2m_obs .and. sfctype) ) then + if (twodvar_regional .or. (hofx_2m_sfcfile .and. sfctype) ) then duplogic=data(ilat,k) == data(ilat,l) .and. & data(ilon,k) == data(ilon,l) .and. & data(ier,k) < r1000 .and. data(ier,l) < r1000 .and. & @@ -468,29 +471,17 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav end do end do -! Run a buddy-check (does not work for use_2m_obs, +! Run a buddy-check (does not work for hofx_2m_sfcfile, ! pre-existing routine crashes with an MPI problem ) ! Note: current params have buddy radius of 108 km, max diff of 8 K. ! gross error check removes O-F > 7., so this is probably removing ! most obs that fail the buddy check already ! GSD uses the buddy check to expand gross error check for obs that pass ! the buddy check. Not sure if we want to do this globally. Turn off for now. - !if ( (twodvar_regional .or. use_2m_obs) .and. buddycheck_t) & + !if ( (twodvar_regional .or. hofx_2m_sfcfile) .and. buddycheck_t) & if ( (twodvar_regional) .and. buddycheck_t) & call buddy_check_t(is,data,luse,mype,nele,nobs,muse,buddyuse) -! for 2m obs, remove obs that failed the buddcheck - ! if (use_2m_obs) then - ! do i = 1,nobs - ! ikx=nint(data(ikxx,k)) - ! itype=ictype(ikx) - ! sfctype=(itype>179.and.itype<190).or.(itype>=192.and.itype<=199) - ! if (sfctype) then - ! if ( buddycheck_t .and. buddyuse(i) == -1) muse(i) = .false. - ! endif - ! enddo - ! endif - ! If requested, save select data for output to diagnostic file if(conv_diagsave)then ii=0 @@ -674,7 +665,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav call tintrp2a1(ges_lnprsl,prsltmp,dlat,dlon,dtime,hrdifsig,& nsig,mype,nfldsig) - if ( .not. use_2m_obs) then + if ( .not. hofx_2m_sfcfile) then drpx=zero if(sfctype .and. .not.twodvar_regional) then drpx=abs(one-((one/exp(dpres-log(psges))))**rd_over_cp)*t0c @@ -720,8 +711,8 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav prsltmp2(2), tvtmp(2), qtmp(2), hsges(1), roges, msges, & f10ges,u10ges,v10ges, t2ges, q2ges, regime, iqtflg) tges = t2ges - elseif (sfctype .and. use_2m_obs ) then -! SCENARIO 2: obs is sfctype, and use_2m_obs scheme is on. + elseif (sfctype .and. hofx_2m_sfcfile ) then +! SCENARIO 2: obs is sfctype, and hofx_2m_sfcfile scheme is on. ! 2m forecast has been read from the sfc guess files ! mask: 0 - sea, 1 - land, 2-ice (val + 3 = interpolated from at least one different land cover @@ -758,7 +749,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ! around this by setting val(2) to zero. endif else -! SCENARIO 3: obs is sfctype, and neither sfcmodel nor use_2m_obs is chosen +! SCENARIO 3: obs is sfctype, and neither sfcmodel nor hofx_2m_sfcfile is chosen ! .or. obs is not sfctype ! ! SCENARIO 3a: obs is a virtual temp. if(iqtflg)then @@ -842,13 +833,13 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav call grdcrd1(sfcchk,prsltmp(1),nsig,-1) ! Check to see if observations is above the top of the model (regional mode) - if(sfctype .and. .not. use_2m_obs )then + if(sfctype .and. .not. hofx_2m_sfcfile )then if(abs(dpres)>four) drpx=1.0e10_r_kind pres_diff=prest-r10*psges if (twodvar_regional .and. abs(pres_diff)>=r1000) drpx=1.0e10_r_kind end if - if (.not. (use_2m_obs .and. sfctype) ) then + if (.not. (hofx_2m_sfcfile .and. sfctype) ) then rlow=max(sfcchk-dpres,zero) ! sfcchk = k of sfc, dpres k of obs ! linear variation of observation ramp [between grid points 1(~3mb) and 15(~45mb) below the surface] if(l_sfcobserror_ramp_t) then @@ -873,7 +864,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ! inflate error for uncertainty in the terrain adjustment lapse_error = 0. - if ( use_2m_obs .and. sfctype) then + if ( hofx_2m_sfcfile .and. sfctype) then if (abs(delta_z)300 m. do not assim. ! inflate obs error to account for error in lapse_rate ! also include some representativity error here (assuming @@ -884,7 +875,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav endif endif -! note: for use_2m_obs, three middle terms in denominator will be zero +! note: for hofx_2m_sfcfile, three middle terms in denominator will be zero ! (elevation mistmatch between obs and model dealt with elsewhere) ratio_errors=error/(data(ier,i)+drpx+1.0e6_r_kind*rhgh+r8*ramp + lapse_error) @@ -892,7 +883,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ! Compute innovation !if(i_use_2mt4b>0 .and. sfctype) then - if(sfctype .and. (( i_use_2mt4b>0) .or. use_2m_obs ) ) then + if(sfctype .and. (( i_use_2mt4b>0) .or. hofx_2m_sfcfile ) ) then ddiff = tob-tges2m else ddiff = tob-tges @@ -1500,7 +1491,7 @@ subroutine init_vars_ write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus call stop2(999) endif - if(i_use_2mt4b>0 .or. use_2m_obs) then + if(i_use_2mt4b>0 .or. hofx_2m_sfcfile) then ! get t2m ... varname='t2m' call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus) @@ -1746,7 +1737,7 @@ subroutine contents_netcdf_diag_(odiag) call nc_diag_metadata("Errinv_Input", sngl(errinv_input) ) call nc_diag_metadata("Errinv_Adjust", sngl(errinv_adjst) ) call nc_diag_metadata("Errinv_Final", sngl(errinv_final) ) - if (use_2m_obs ) then + if (hofx_2m_sfcfile ) then call nc_diag_metadata("Observation", sngl(tob) ) call nc_diag_metadata("Observation_Before_Elev_Correction", sngl(data(itob,i)) ) else @@ -1790,7 +1781,7 @@ subroutine contents_netcdf_diag_(odiag) call nc_diag_data2d("ObsDiagSave_obssen", odiag%obssen ) endif - if (twodvar_regional .or. use_2m_obs ) then + if (twodvar_regional .or. hofx_2m_sfcfile ) then call nc_diag_metadata("Dominant_Sfc_Type", data(idomsfc,i) ) ! this is the model height interpolated to the obs location in read_prepbufr call nc_diag_metadata("Model_Terrain", data(izz,i) ) diff --git a/src/gsi/update_guess.f90 b/src/gsi/update_guess.f90 index c8656053b5..8c7e40459a 100644 --- a/src/gsi/update_guess.f90 +++ b/src/gsi/update_guess.f90 @@ -113,7 +113,7 @@ subroutine update_guess(sval,sbias) use mpimod, only: mype use constants, only: zero,one,fv,max_varname_length,qmin,qcmin,tgmin,& r100,one_tenth,tiny_r_kind - use jfunc, only: iout_iter,bcoption,tsensible,clip_supersaturation,superfact,use_2m_obs + use jfunc, only: iout_iter,bcoption,tsensible,clip_supersaturation,superfact use gridmod, only: lat2,lon2,nsig,& regional,twodvar_regional,regional_ozone,& l_reg_update_hydro_delz @@ -454,7 +454,7 @@ subroutine update_guess(sval,sbias) endif call gsd_update_soil_tq(tinc_1st,is_t,qinc_1st,is_q,it) endif ! l_gsd_soilTQ_nudge - if ( (i_use_2mt4b > 0.or. use_2m_obs) .and. is_t>0) then + if ( (i_use_2mt4b > 0) .and. is_t>0) then do j=1,lon2 do i=1,lat2 tinc_1st(i,j)=p_tv(i,j,1) @@ -462,7 +462,7 @@ subroutine update_guess(sval,sbias) end do call gsd_update_t2m(tinc_1st,it) endif ! l_gsd_t2m_adjust - if ( (i_use_2mq4b > 0.or. use_2m_obs) .and. is_q>0) then + if ( (i_use_2mq4b > 0) .and. is_q>0) then do j=1,lon2 do i=1,lat2 qinc_1st(i,j)=p_q(i,j,1)