From 3c361da3b3cc86513586db39dcd8c7adf15c0030 Mon Sep 17 00:00:00 2001 From: Emily Liu Date: Thu, 31 Mar 2022 02:17:58 +0000 Subject: [PATCH 1/3] GitHub Issue NOAA-EMC/GSI#341 This update includes the following: (1) Add the capability of assimilating precipitation-affected radiances from AMSU-A and ATMS (2) Enable the all-sky radiance assimilation for GMI (3) Add the switch to use brightness temperature (SDR) instead of antenna temperature (TDR) for AMSU-A, ATMS and MHS --- scripts/exglobal_atmos_analysis.sh | 12 +- src/gsi/antcorr_application.f90 | 104 +++ src/gsi/cloud_efr_mod.f90 | 2 +- src/gsi/cplr_gfs_ensmod.f90 | 56 +- src/gsi/crtm_interface.f90 | 192 +++-- src/gsi/general_read_gfsatm.f90 | 1086 ++++++++++++++++++++++++++++ src/gsi/gsimod.F90 | 6 +- src/gsi/netcdfgfs_io.f90 | 119 ++- src/gsi/obsmod.F90 | 4 + src/gsi/prt_guess.f90 | 2 +- src/gsi/qcmod.f90 | 46 +- src/gsi/read_atms.f90 | 8 +- src/gsi/read_bufrtovs.f90 | 195 +++-- src/gsi/read_gmi.f90 | 50 +- src/gsi/setuprad.f90 | 34 +- src/gsi/write_incr.f90 | 145 +++- 16 files changed, 1845 insertions(+), 216 deletions(-) diff --git a/scripts/exglobal_atmos_analysis.sh b/scripts/exglobal_atmos_analysis.sh index 87a0ab4f53..5aa43fc5a6 100755 --- a/scripts/exglobal_atmos_analysis.sh +++ b/scripts/exglobal_atmos_analysis.sh @@ -81,6 +81,10 @@ export imp_physics=${imp_physics:-99} lupp=${lupp:-".true."} cnvw_option=${cnvw_option:-".false."} +# Observation usage options +cao_check=${cao_check:-".false."} +ta2tb=${ta2tb:-".false."} + # Diagnostic files options lobsdiag_forenkf=${lobsdiag_forenkf:-".false."} netcdf_diag=${netcdf_diag:-".true."} @@ -461,6 +465,7 @@ for file in $(awk '{if($1!~"!"){print $1}}' satinfo | sort | uniq); do $NLN $RTMFIX/${file}.SpcCoeff.bin ./crtm_coeffs/ $NLN $RTMFIX/${file}.TauCoeff.bin ./crtm_coeffs/ done +$NLN $RTMFIX/amsua_metop-a_v2.SpcCoeff.bin ./crtm_coeffs/amsua_metop-a_v2.SpcCoeff.bin $NLN $RTMFIX/Nalli.IRwater.EmisCoeff.bin ./crtm_coeffs/Nalli.IRwater.EmisCoeff.bin $NLN $RTMFIX/NPOESS.IRice.EmisCoeff.bin ./crtm_coeffs/NPOESS.IRice.EmisCoeff.bin @@ -473,6 +478,8 @@ $NLN $RTMFIX/NPOESS.VISwater.EmisCoeff.bin ./crtm_coeffs/NPOESS.VISwater.EmisCoe $NLN $RTMFIX/FASTEM6.MWwater.EmisCoeff.bin ./crtm_coeffs/FASTEM6.MWwater.EmisCoeff.bin $NLN $RTMFIX/AerosolCoeff.bin ./crtm_coeffs/AerosolCoeff.bin $NLN $RTMFIX/CloudCoeff.bin ./crtm_coeffs/CloudCoeff.bin +#$NLN $RTMFIX/CloudCoeff.GFDLFV3.-109z-1.bin ./crtm_coeffs/CloudCoeff.bin + ############################################################## # Observational data @@ -759,11 +766,12 @@ cat > gsiparm.anl << EOF crtm_coeffs_path='./crtm_coeffs/', newpc4pred=.true.,adp_anglebc=.true.,angord=4,passive_bc=.true.,use_edges=.false., diag_precon=.true.,step_start=1.e-3,emiss_bc=.true.,nhr_obsbin=${nhr_obsbin:-3}, - cwoption=3,imp_physics=$imp_physics,lupp=$lupp,cnvw_option=$cnvw_option, + cwoption=3,imp_physics=$imp_physics,lupp=$lupp,cnvw_option=$cnvw_option,cao_check=${cao_check}, netcdf_diag=$netcdf_diag,binary_diag=$binary_diag, lobsdiag_forenkf=$lobsdiag_forenkf, write_fv3_incr=$write_fv3_increment, nhr_anal=${IAUFHRS}, + ta2tb=${ta2tb}, $WRITE_INCR_ZERO $WRITE_ZERO_STRAT $WRITE_STRAT_EFOLD @@ -901,7 +909,7 @@ OBS_INPUT:: avhambufr avhrr metop-c avhrr3_metop-c 0.0 4 0 avhpmbufr avhrr n19 avhrr3_n19 0.0 4 0 amsr2bufr amsr2 gcom-w1 amsr2_gcom-w1 0.0 3 0 - gmibufr gmi gpm gmi_gpm 0.0 3 0 + gmibufr gmi gpm gmi_gpm 0.0 1 0 saphirbufr saphir meghat saphir_meghat 0.0 3 0 ahibufr ahi himawari8 ahi_himawari8 0.0 1 0 abibufr abi g16 abi_g16 0.0 1 0 diff --git a/src/gsi/antcorr_application.f90 b/src/gsi/antcorr_application.f90 index db2808d08b..7b625a1096 100644 --- a/src/gsi/antcorr_application.f90 +++ b/src/gsi/antcorr_application.f90 @@ -10,6 +10,8 @@ ! paul.vandelst@noaa.gov ! ! 2011-04-25 A.Collard Modified to be consistent with CRTM 2.1 +! 2020-04-21 S.Sieron Added Apply_Antcorr so that all assimlated +! observations be brightness temperatures ! MODULE AntCorr_Application @@ -35,6 +37,7 @@ MODULE AntCorr_Application ! The AntCorr structure definition PUBLIC :: ACCoeff_type PUBLIC :: Remove_AntCorr + PUBLIC :: Apply_AntCorr ! ----------------- ! Module parameters @@ -153,5 +156,106 @@ SUBROUTINE Remove_AntCorr( AC , & ! Input T(l) = AC%A_earth(iFOV,l)*T(l) + AC%A_platform(iFOV,l)*T(l) + AC%A_space(iFOV,l)*TSPACE END DO END SUBROUTINE Remove_AntCorr + + +!-------------------------------------------------------------------------------- +! +! NAME: +! Apply_AntCorr +! +! PURPOSE: +! Subroutine to apply an antenna correction to microwave instrument +! antenna temperatures, Ta, to produce brightness temperatures, Tb. +! +! CALLING SEQUENCE: +! CALL Apply_AntCorr( AC , & ! Input +! iFOV, & ! Input +! T ) ! In/Output +! +! INPUT ARGUMENTS: +! AC: Structure containing the antenna correction coefficients +! for the sensor of interest. +! UNITS: N/A +! TYPE: TYPE(ACCoeff_type) +! DIMENSION: Scalar +! ATTRIBUTES: INTENT(IN) +! +! iFOV: The FOV index for a scanline of the sensor of interest. +! UNITS: N/A +! TYPE: INTEGER +! DIMENSION: Scalar +! ATTRIBUTES: INTENT(IN) +! +! T: On input, this argument contains the antenna +! temperatures for the sensor channels. +! UNITS: Kelvin +! TYPE: REAL(fp) +! DIMENSION: Rank-1 (n_Channels) +! ATTRIBUTES: INTENT(IN OUT) +! +! OUTPUT ARGUMENTS: +! T: On output, this argument contains the brightness +! temperatures for the sensor channels. +! If an error occurs, the return values are all -1. +! UNITS: Kelvin +! TYPE: REAL(fp) +! DIMENSION: Rank-1 (n_Channels) +! ATTRIBUTES: INTENT(IN OUT) +! +! SIDE EFFECTS: +! The temperature array argument, T, is modified. +! +! PROCEDURE: +! For every FOV and channel, the antenna temperature, Ta, is converted +! to brightness temperature, Tb, using, +! +! Tb = (Ta - As.Ts) / (Ae + Ap) +! +! where Ae == antenna efficiency for the Earth view +! Ap == antenna efficiency for satellite platform view +! As == antenna efficiency for cold space view +! Ts == cosmic background temperature. +! +! Note that the observed earth view brightness temperature is used as a +! proxy for the platform temperature for the (Ap.Tb) component since +! there is no measurement of the platform temperature in-flight. +! +!-------------------------------------------------------------------------------- + + SUBROUTINE Apply_AntCorr( AC , & ! Input + iFOV, & ! Input + T ) ! In/Output + + use kinds,only: i_kind, r_kind + + implicit none + + ! Arguments + TYPE(ACCoeff_type), INTENT(IN) :: AC + INTEGER(i_kind) , INTENT(IN) :: iFOV + REAL(r_kind) , INTENT(IN OUT) :: T(:) + ! Local parameters + CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'Apply_AntCorr' + ! Local variables + INTEGER(i_kind) :: l + + ! Check input + IF ( iFOV < 1 .OR. iFOV > AC%n_FOVS ) THEN + CALL Display_Message( ROUTINE_NAME, 'Input iFOV inconsistent with AC data', FAILURE ) + T = INVALID + RETURN + END IF + IF ( SIZE(T) /= AC%n_Channels ) THEN + CALL Display_Message( ROUTINE_NAME, 'Size of T() inconsistent with AC data', FAILURE ) + T = INVALID + RETURN + END IF + ! Compute the brightness temperature + DO l = 1, AC%n_Channels + T(l) = (T(l) - AC%A_space(iFOV,l)*TSPACE) / (AC%A_earth(iFOV,l) + & + AC%A_platform(iFOV,l)) + END DO + END SUBROUTINE Apply_AntCorr + END MODULE AntCorr_Application diff --git a/src/gsi/cloud_efr_mod.f90 b/src/gsi/cloud_efr_mod.f90 index e965d29e38..6fe17fb903 100644 --- a/src/gsi/cloud_efr_mod.f90 +++ b/src/gsi/cloud_efr_mod.f90 @@ -285,7 +285,7 @@ subroutine cloud_calc_gfs(g_ql,g_qi,g_cwmr,g_q,g_tv,lower_bound,g_cf) real(r_kind),dimension(lat2,lon2,nsig),intent(inout):: g_cwmr ! mixing ratio of total condensates [Kg/Kg] real(r_kind),dimension(lat2,lon2,nsig),intent(in ):: g_q ! specific humidity [Kg/Kg] real(r_kind),dimension(lat2,lon2,nsig),intent(in ):: g_tv ! virtual temperature [K] - real(r_kind),dimension(lat2,lon2,nsig),intent(inout), optional:: g_cf ! cloud fractio + real(r_kind),dimension(lat2,lon2,nsig),intent(inout), optional:: g_cf ! cloud fraction logical,intent(in):: lower_bound ! If .true., set lower bound to cloud ! Declare local variables diff --git a/src/gsi/cplr_gfs_ensmod.f90 b/src/gsi/cplr_gfs_ensmod.f90 index ac00db4c41..c9aba6f9f0 100644 --- a/src/gsi/cplr_gfs_ensmod.f90 +++ b/src/gsi/cplr_gfs_ensmod.f90 @@ -641,7 +641,7 @@ subroutine parallel_read_nemsio_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsi ! Declare local variables integer(i_kind) i,ii,j,jj,k,lonb,latb,levs,latb2,lonb2 integer(i_kind) k2,k3,k3u,k3v,k3t,k3q,k3cw,k3oz,kf - integer(i_kind) k3ql,k3qi,k3qr,k3qs,k3qg + integer(i_kind) k3ql,k3qi,k3qr,k3qs,k3qg integer(i_kind) iret integer(i_kind) :: istop = 101 integer(i_kind),dimension(7):: idate @@ -916,6 +916,7 @@ subroutine parallel_read_gfsnc_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsig ! Declare local variables integer(i_kind) i,ii,j,jj,k,lonb,latb,levs,kr,ierror integer(i_kind) k2,k3,k3u,k3v,k3t,k3q,k3cw,k3oz,kf + integer(i_kind) k3ql,k3qi,k3qr,k3qs,k3qg character(len=120) :: myname_ = 'parallel_read_gfsnc_state_' real(r_single),allocatable,dimension(:,:,:) :: rwork3d1, rwork3d2 real(r_single),allocatable,dimension(:,:) :: temp2,rwork2d @@ -959,6 +960,7 @@ subroutine parallel_read_gfsnc_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsig allocate(rwork3d1(nlon,(nlat-2),nsig)) allocate(temp3(nlat,nlon,nsig,nc3d)) k3u=0 ; k3v=0 ; k3t=0 ; k3q=0 ; k3cw=0 ; k3oz=0 + k3ql=0; k3qi=0; k3qr=0; k3qs=0; k3qg=0 do k3=1,nc3d if (trim(cvars3d(k3))=='cw') then k3cw=k3 @@ -973,6 +975,46 @@ subroutine parallel_read_gfsnc_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsig call move1_(rwork3d1(:,:,kr),temp3(:,:,k,k3),nlon,nlat) call fillpoles_ss_(temp3(:,:,k,k3),nlon,nlat) end do + else if(trim(cvars3d(k3))=='ql') then + k3ql=k3 + call read_vardata(atmges, 'clwmr', rwork3d1) + do k=1,nsig + kr = levs+1-k + call move1_(rwork3d1(:,:,kr),temp3(:,:,k,k3),nlon,nlat) + call fillpoles_ss_(temp3(:,:,k,k3),nlon,nlat) + end do + else if(trim(cvars3d(k3))=='qi') then + k3qi=k3 + call read_vardata(atmges, 'icmr', rwork3d1) + do k=1,nsig + kr = levs+1-k + call move1_(rwork3d1(:,:,kr),temp3(:,:,k,k3),nlon,nlat) + call fillpoles_ss_(temp3(:,:,k,k3),nlon,nlat) + end do + else if(trim(cvars3d(k3))=='qr') then + k3qr=k3 + call read_vardata(atmges, 'rwmr', rwork3d1) + do k=1,nsig + kr = levs+1-k + call move1_(rwork3d1(:,:,kr),temp3(:,:,k,k3),nlon,nlat) + call fillpoles_ss_(temp3(:,:,k,k3),nlon,nlat) + end do + else if(trim(cvars3d(k3))=='qs') then + k3qs=k3 + call read_vardata(atmges, 'snmr', rwork3d1) + do k=1,nsig + kr = levs+1-k + call move1_(rwork3d1(:,:,kr),temp3(:,:,k,k3),nlon,nlat) + call fillpoles_ss_(temp3(:,:,k,k3),nlon,nlat) + end do + else if(trim(cvars3d(k3))=='qg') then + k3qg=k3 + call read_vardata(atmges, 'grle', rwork3d1) + do k=1,nsig + kr = levs+1-k + call move1_(rwork3d1(:,:,kr),temp3(:,:,k,k3),nlon,nlat) + call fillpoles_ss_(temp3(:,:,k,k3),nlon,nlat) + end do else if(trim(cvars3d(k3))=='oz') then k3oz=k3 call read_vardata(atmges, 'o3mr', rwork3d1) @@ -1015,7 +1057,8 @@ subroutine parallel_read_gfsnc_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsig enddo deallocate(rwork3d1) - if (k3u==0.or.k3v==0.or.k3t==0.or.k3q==0.or.k3cw==0.or.k3oz==0) & +! if (k3u==0.or.k3v==0.or.k3t==0.or.k3q==0.or.k3cw==0.or.k3oz==0) & + if (k3u==0.or.k3v==0.or.k3t==0.or.k3q==0.or.k3oz==0) & write(6,'(" WARNING, problem with one of k3-")') do k=1,nsig @@ -1316,8 +1359,13 @@ subroutine get_gfs_ens(this,grd,member,ntindex,atm_bundle,iret) endif else if ( use_gfs_ncio ) then - call general_read_gfsatm_nc(grd,sp_ens,filename,uv_hyb_ens,.false., & - zflag,atm_bundle,.true.,iret) + if (fv3_full_hydro) then + call general_read_gfsatm_allhydro_nc(grd,sp_ens,filename,uv_hyb_ens,.false., & + zflag,atm_bundle,iret) + else + call general_read_gfsatm_nc(grd,sp_ens,filename,uv_hyb_ens,.false., & + zflag,atm_bundle,.true.,iret) + endif else call general_read_gfsatm(grd,sp_ens,sp_ens,filename,uv_hyb_ens,.false., & zflag,atm_bundle,inithead,iret) diff --git a/src/gsi/crtm_interface.f90 b/src/gsi/crtm_interface.f90 index 3de679f7f4..e1b8e73d8b 100644 --- a/src/gsi/crtm_interface.f90 +++ b/src/gsi/crtm_interface.f90 @@ -200,6 +200,7 @@ module crtm_interface logical ,save :: cld_sea_only_wk logical ,save :: lprecip_wk logical ,save :: mixed_use + logical ,save :: use_gfdl_qsat integer(i_kind), parameter :: min_n_absorbers = 2 integer(i_kind),save :: iedge_log @@ -341,7 +342,7 @@ subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmo ! local parameters character(len=*), parameter :: myname_=myname//'*init_crtm' - integer(i_kind), parameter :: length = 2621 ! lenth of GFL qsat table + integer(i_kind), parameter :: length = 2621 ! length of GFL qsat table ! local variables integer(i_kind) :: ier,ii,error_status @@ -356,7 +357,7 @@ subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmo logical quiet logical print_verbose - + use_gfdl_qsat=.false. print_verbose=.false. if(verbose)print_verbose=.true. isst=-1 @@ -838,27 +839,29 @@ subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmo endif ! regional or IGBP ! Initial GFDL saturation water vapor pressure tables - if (n_actual_aerosols_wk>0 .or. n_clouds_fwd_wk>0 .and. imp_physics==11) then + if (use_gfdl_qsat) then + if (n_actual_aerosols_wk>0 .or. n_clouds_fwd_wk>0 .and. imp_physics==11) then - if (mype==0) write(6,*)myname_,':initial and load GFDL saturation water vapor pressure tables' + if (mype==0) write(6,*)myname_,':initial and load GFDL saturation water vapor pressure tables' - allocate(table (length)) - allocate(table2(length)) - allocate(tablew(length)) - allocate(des2 (length)) - allocate(desw (length)) - - call qs_table (length) - call qs_table2(length) - call qs_tablew(length) - - do ii = 1, length - 1 - des2 (ii) = max (zero, table2 (ii + 1) - table2 (ii)) - desw (ii) = max (zero, tablew (ii + 1) - tablew (ii)) - enddo - des2 (length) = des2 (length - 1) - desw (length) = desw (length - 1) + allocate(table (length)) + allocate(table2(length)) + allocate(tablew(length)) + allocate(des2 (length)) + allocate(desw (length)) + + call qs_table (length) + call qs_table2(length) + call qs_tablew(length) + + do ii = 1, length - 1 + des2 (ii) = max (zero, table2 (ii + 1) - table2 (ii)) + desw (ii) = max (zero, tablew (ii + 1) - tablew (ii)) + enddo + des2 (length) = des2 (length - 1) + desw (length) = desw (length - 1) + endif endif return @@ -891,14 +894,12 @@ subroutine destroy_crtm error_status = crtm_destroy(channelinfo) if (error_status /= success) & write(6,*)myname_,': ***ERROR*** error_status=',error_status - if (n_actual_aerosols_wk>0 .or. n_clouds_fwd_wk>0) then - if (imp_physics==11) then - deallocate(table) - deallocate(table2) - deallocate(tablew) - deallocate(des2) - deallocate(desw) - endif + if (use_gfdl_qsat) then + deallocate(table) + deallocate(table2) + deallocate(tablew) + deallocate(des2) + deallocate(desw) endif call crtm_atmosphere_destroy(atmosphere(1)) call crtm_surface_destroy(surface(1)) @@ -1054,6 +1055,11 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & use obsmod, only: iadate use aeroinfo, only: nsigaerojac use chemmod, only: lread_ext_aerosol !for separate aerosol input file + use crtm_cloudcover_define, only: cloudcover_maximum_overlap, & + cloudcover_random_overlap, & + cloudcover_maxran_overlap, & + cloudcover_average_overlap, & !default + cloudcover_overcast_overlap implicit none @@ -1106,6 +1112,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & integer(i_kind):: idx700,dprs,dprs_min integer(i_kind),dimension(8)::obs_time,anal_time integer(i_kind),dimension(msig) :: klevel + real(r_kind),dimension(nsig) :: qsat ! ****************************** ! Constrained indexing for lai @@ -1172,6 +1179,14 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & if (present(tcwv)) tcwv=zero if (present(tcc)) tcc=zero + if (n_clouds_fwd_wk>0) then + cloud = zero + cloud_cont = zero + cloud_efr = zero + cloudefr = zero + cf = zero + endif + dx = data_s(ilat) ! grid relative latitude dy = data_s(ilon) ! grid relative longitude hs = data_s(izz) ! surface height @@ -1682,6 +1697,16 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & else q(k) = qsmall endif + if (n_clouds_fwd_wk>0) then + qsat(k)=(ges_qsat (ix ,iy ,k, itsig)*w00+ & + ges_qsat (ixp,iy ,k, itsig)*w10+ & + ges_qsat (ix ,iyp,k, itsig)*w01+ & + ges_qsat (ixp,iyp,k, itsig)*w11)*dtsig + & + (ges_qsat (ix ,iy ,k, itsigp)*w00+ & + ges_qsat (ixp,iy ,k, itsigp)*w10+ & + ges_qsat (ix ,iyp,k, itsigp)*w01+ & + ges_qsat (ixp,iyp,k, itsigp)*w11)*dtsigp + end if c2(k)=one/(one+fv*q(k)) c3(k)=one/(one-q(k)) c4(k)=fv*h(k)*c2(k) @@ -1778,9 +1803,10 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & endif ! Calculate GFDL cloud fraction (if no cf in metguess table) based on PDF scheme - if ( icmask .and. n_clouds_fwd_wk > 0 .and. imp_physics==11 .and. lcalc_gfdl_cfrac ) then +! if ( icmask .and. n_clouds_fwd_wk > 0 .and. imp_physics==11 .and. lcalc_gfdl_cfrac ) then + if ( icmask .and. n_clouds_fwd_wk > 0 .and. imp_physics==11 .and. lprecip_wk ) then cf_calc = zero - call calc_gfdl_cloudfrac(rho_air,h,qmix,cloud,hs,garea,cf_calc) + call calc_gfdl_cloudfrac(rho_air,h,q,cloud,hs,garea,qsat,cf_calc) cf = cf_calc icfs = 0 ! load cloud fraction into CRTM endif @@ -1972,7 +1998,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & if (n_clouds_fwd_wk>0) then kgkg_kgm2=(atmosphere(1)%level_pressure(k)-atmosphere(1)%level_pressure(k-1))*r100/grav - if (cw_cv.or.ql_cv) then + if ((cw_cv.or.ql_cv).and.(.not. lprecip_wk)) then if (icmask) then c6(k) = kgkg_kgm2 auxdp(k)=abs(prsi_rtm(kk+1)-prsi_rtm(kk))*r10 @@ -2000,7 +2026,6 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & if (ii==2 .and. atmosphere(1)%temperature(k) 1.0e-6_r_kind) then - cloud_efr (k,ii)=cloudefr(kk2,ii) + cloud_efr(k,ii)=cloudefr(kk2,ii) else - cloud_efr (k,ii)=zero + cloud_efr(k,ii)=zero endif enddo - + +! In CRTM, if cloud fraction of the layer < 1.0E-12, set cloud content and +! effective radius of all hydrometer types in that layer to zero +! CRTM minimum thresholds: cloud content=1.0E-6 and cloud fraction=1.E-12 + do ii=1,n_clouds_fwd_wk + if(cloud_cont(k,ii) > 1.000_r_kind*1.0E-6_r_kind .and. cf(kk2) < 1.000_r_kind*1.0E-12_r_kind) then + cf(kk2)=1.001_r_kind*1.0E-12_r_kind + end if + end do + if (cloud_cont(k,1) >= 1.0e-6_r_kind) clw_guess = clw_guess + cloud_cont(k,1) - tcwv = tcwv + (atmosphere(1)%absorber(k,1)*0.001_r_kind)*c6(k) + if (present(tcwv)) tcwv = tcwv + (atmosphere(1)%absorber(k,1)*0.001_r_kind)*c6(k) do ii=1,n_clouds_fwd_wk if (cloud_cont(k,ii) >= 1.0e-6_r_kind) hwp_guess(ii) = hwp_guess(ii) + cloud_cont(k,ii) enddo @@ -2024,18 +2058,36 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & !Add lower bound to all hydrometers !note: may want to add lower bound value for effective radius do ii=1,n_clouds_fwd_wk - if (trim(cloud_names_fwd(ii))=='ql' .and. atmosphere(1)%temperature(k)-t0c>-20.0_r_kind) & + if (trim(cloud_names_fwd(ii))=='ql' .and. atmosphere(1)%temperature(k)-t0c>-20.0_r_kind) then cloud_cont(k,ii)=max(1.001_r_kind*1.0E-6_r_kind, cloud_cont(k,ii)) - if (trim(cloud_names_fwd(ii))=='qi' .and. atmosphere(1)%temperature(k)-20.0_r_kind) & + cloud_efr(k,ii)=max(5.001_r_kind, cloud_efr(k,ii)) + endif + if (trim(cloud_names_fwd(ii))=='qr' .and. atmosphere(1)%temperature(k)-t0c>-20.0_r_kind) then cloud_cont(k,ii)=max(1.001_r_kind*1.0E-6_r_kind, cloud_cont(k,ii)) - if (trim(cloud_names_fwd(ii))=='qs' .and. atmosphere(1)%temperature(k) 1.000_r_kind*1.0E-6_r_kind .and. atmosphere(1)%cloud_fraction(k) < 1.001_r_kind*1.0E-12_r_kind) then + atmosphere(1)%cloud_fraction(k)=1.001_r_kind*1.0E-12_r_kind + end if end do -!crtm2.3.x if (.not. regional .and. icfs==0 ) atmosphere(1)%cloud_fraction(k) = cf(kk2) end if endif endif @@ -2056,7 +2108,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & end do if (n_clouds_fwd_wk>0 .and. icmask) then - if ((hwp_guess(1)+hwp_guess(2))>=1.0e-06_r_kind) hwp_ratio = hwp_guess(1)/(hwp_guess(1)+hwp_guess(2)) + if ((hwp_guess(1)+hwp_guess(2))>=1.0e-06_r_kind .and. present(hwp_ratio)) hwp_ratio = hwp_guess(1)/(hwp_guess(1)+hwp_guess(2)) hwp_total = sum(hwp_guess(:)) theta_700 = atmosphere(1)%temperature(idx700)*(r1000/atmosphere(1)%pressure(idx700))**rd_over_cp theta_sfc = data_s(itsavg)*(r100/ps)**rd_over_cp @@ -2158,7 +2210,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & ! Simulated brightness temperatures tsim(i)=rtsolution(i,1)%brightness_temperature -! if (present(tcc)) tcc(i)=rtsolution(i,1)%total_cloud_cover !crtm2.3.x + if (present(tcc)) tcc(i)=rtsolution(i,1)%total_cloud_cover !crtm2.3.x if (n_clouds_fwd_wk>0 .and. present(tsim_clr)) then if (mixed_use) then @@ -2509,7 +2561,7 @@ subroutine calc_gfdl_reff(rho_air,tsen,qxmr,cloud_name,reff) end subroutine calc_gfdl_reff - subroutine calc_gfdl_cloudfrac(den,pt1,qv,cloud,hs,area,cfrac) + subroutine calc_gfdl_cloudfrac(den,pt1,qv,cloud,hs,area,qsat,cfrac) !$$$ subprogram documentation block ! . . . . ! subprogram: calc_gfdl_cloudfrac calculate GFDL cloud fraction @@ -2552,6 +2604,7 @@ subroutine calc_gfdl_cloudfrac(den,pt1,qv,cloud,hs,area,cfrac) real(r_kind), dimension(nsig,n_clouds_fwd_wk) ,intent(in ) :: cloud ! hydroeteor mixing ratio real(r_kind), intent(in ) :: hs ! surface elevation [ m ] real(r_kind), intent(in ) :: area ! analysis grid area [ m2 ] + real(r_kind), dimension(nsig) ,intent(in ) :: qsat ! saturation specific humidity real(r_kind), dimension(nsig) ,intent(inout) :: cfrac ! cloud fraction ! ! Declare local variables @@ -2572,6 +2625,7 @@ subroutine calc_gfdl_cloudfrac(den,pt1,qv,cloud,hs,area,cfrac) real(r_kind) :: dw,dw_land,dw_ocean real(r_kind) :: q_sol,q_liq real(r_kind) :: qa + real(r_kind) :: phi real(r_kind), dimension(nsig) :: ql,qi,qr,qs,qg logical :: hydrostatic @@ -2604,7 +2658,9 @@ subroutine calc_gfdl_cloudfrac(den,pt1,qv,cloud,hs,area,cfrac) ! default area dependent form: use dx ~ 100 km as the base ! ----------------------------------------------------------------------- ! higher than 10 m is considered "land" and will have higher subgrid variability - dw = dw_ocean + (dw_land - dw_ocean) * min (one, abs(hs) / (ten * grav)) + phi = hs * grav ! surface potential + dw = dw_ocean + (dw_land - dw_ocean) * min (one, abs(phi) / (ten * grav)) +! dw = dw_ocean + (dw_land - dw_ocean) * min (one, abs(hs) / (ten * grav)) ! "scale - aware" subgrid variability: 100 - km as the base hvar = min (0.2_r_kind, max (0.01_r_kind, dw * sqrt (sqrt (area) / 100.e3_r_kind))) @@ -2641,25 +2697,29 @@ subroutine calc_gfdl_cloudfrac(den,pt1,qv,cloud,hs,area,cfrac) ! ----------------------------------------------------------------------- ! determine saturated specific humidity ! ----------------------------------------------------------------------- - if (tin <= t_wfr) then - ! ice phase: - qstar = iqs1 (tin, den(k)) - elseif (tin >= tice) then - ! liquid phase: - qstar = wqs1 (tin, den(k)) - else - ! mixed phase: - qsi = iqs1 (tin, den(k)) - qsw = wqs1 (tin, den(k)) - if (q_cond > 1.e-6_r_kind) then - rqi = q_sol / q_cond + if (use_gfdl_qsat) then + if (tin <= t_wfr) then + ! ice phase: + qstar = iqs1 (tin, den(k)) + elseif (tin >= tice) then + ! liquid phase: + qstar = wqs1 (tin, den(k)) else - ! -------------------------------------------------------------- - ! mostly liquid water q_cond at initial cloud development stage - ! -------------------------------------------------------------- - rqi = (tice - tin) / (tice - t_wfr) + ! mixed phase: + qsi = iqs1 (tin, den(k)) + qsw = wqs1 (tin, den(k)) + if (q_cond > 1.e-6_r_kind) then + rqi = q_sol / q_cond + else + ! -------------------------------------------------------------- + ! mostly liquid water q_cond at initial cloud development stage + ! -------------------------------------------------------------- + rqi = (tice - tin) / (tice - t_wfr) + endif + qstar = rqi * qsi + (one - rqi) * qsw endif - qstar = rqi * qsi + (one - rqi) * qsw + else + qstar = qsat(k) endif ! ----------------------------------------------------------------------- @@ -2674,14 +2734,14 @@ subroutine calc_gfdl_cloudfrac(den,pt1,qv,cloud,hs,area,cfrac) ! icloud_f = 1: old fvgfs gfdl) mp implementation ! icloud_f = 2: binary cloud scheme (0 / 1) ! ----------------------------------------------------------------------- - if (rh > 0.75_r_kind .and. qpz > 1.e-6_r_kind) then + if (rh > 0.75_r_kind .and. qpz > 1.e-8_r_kind) then dq = hvar * qpz q_plus = qpz + dq q_minus = qpz - dq if (icloud_f == 2) then if (qpz > qstar) then qa = one - elseif (qstar < q_plus .and. q_cond > 1.e-6_r_kind) then + elseif (qstar < q_plus .and. q_cond > 1.e-8_r_kind) then qa = ((q_plus - qstar) / dq) ** 2 qa = min (one, qa) else @@ -2701,7 +2761,7 @@ subroutine calc_gfdl_cloudfrac(den,pt1,qv,cloud,hs,area,cfrac) qa = zero endif ! impose minimum cloudiness if substantial q_cond exist - if (q_cond > 1.e-6_r_kind) then + if (q_cond > 1.e-8_r_kind) then qa = max (cld_min, qa) endif qa = min (one, qa) diff --git a/src/gsi/general_read_gfsatm.f90 b/src/gsi/general_read_gfsatm.f90 index fb7c568b33..57b6c6033c 100755 --- a/src/gsi/general_read_gfsatm.f90 +++ b/src/gsi/general_read_gfsatm.f90 @@ -188,7 +188,229 @@ subroutine general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr, return end subroutine general_reload +subroutine general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vdflag,g_cf) +! !USES: + + use kinds, only: r_kind,i_kind + use mpimod, only: npe,mpi_comm_world,ierror,mpi_rtype + use general_sub2grid_mod, only: sub2grid_info + implicit none + +! !INPUT PARAMETERS: + + type(sub2grid_info), intent(in ) :: grd + integer(i_kind), intent(inout) :: icount + integer(i_kind),dimension(npe), intent(inout) :: ilev,iflag + real(r_kind),dimension(grd%itotsub),intent(in ) :: work + logical, intent(in ) :: uvflag,vdflag + +! !OUTPUT PARAMETERS: + + real(r_kind),dimension(grd%lat2,grd%lon2), intent( out) :: g_ps + real(r_kind),dimension(grd%lat2,grd%lon2), intent(inout) :: g_z + real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig),intent( out) :: g_u,g_v,& + g_vor,g_div,g_q,g_oz,g_tv,g_ql,g_qi,g_qr,g_qs,g_qg + real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig),intent( out),optional :: g_cf + + +! !DESCRIPTION: Transfer contents of 2-d array global to 3-d subdomain array +! +! !REVISION HISTORY: +! 2004-05-14 treadon +! 2004-07-15 todling, protex-compliant prologue +! 2014-12-03 derber - introduce vdflag and optimize routines +! +! !REMARKS: +! +! language: f90 +! machine: ibm rs/6000 sp; sgi origin 2000; compaq/hp +! +! !AUTHOR: +! treadon org: np23 date: 2004-05-14 +! +!EOP +!------------------------------------------------------------------------- + + integer(i_kind) i,j,k,ij,klev + real(r_kind),dimension(grd%lat2*grd%lon2,npe):: sub + + call mpi_alltoallv(work,grd%sendcounts_s,grd%sdispls_s,mpi_rtype,& + sub,grd%recvcounts_s,grd%rdispls_s,mpi_rtype,& + mpi_comm_world,ierror) + +!$omp parallel do schedule(dynamic,1) private(k,i,j,ij,klev) + do k=1,icount + if ( iflag(k) == 1 ) then + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_z(i,j)=sub(ij,k) + enddo + enddo + elseif ( iflag(k) == 2 ) then + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_ps(i,j)=sub(ij,k) + enddo + enddo + elseif ( iflag(k) == 3 ) then + klev=ilev(k) + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_tv(i,j,klev)=sub(ij,k) + enddo + enddo + elseif ( iflag(k) == 4 ) then + klev=ilev(k) + if ( vdflag ) then + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_vor(i,j,klev)=sub(ij,k) + enddo + enddo + endif + if ( .not. uvflag ) then + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_u(i,j,klev)=sub(ij,k) + enddo + enddo + endif + elseif ( iflag(k) == 5 ) then + klev=ilev(k) + if ( vdflag ) then + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_div(i,j,klev)=sub(ij,k) + enddo + enddo + endif + if ( .not. uvflag ) then + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_v(i,j,klev)=sub(ij,k) + enddo + enddo + endif + elseif ( iflag(k) == 6 ) then + if ( .not. uvflag) then + endif + klev=ilev(k) + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_u(i,j,klev)=sub(ij,k) + enddo + enddo + elseif ( iflag(k) == 7 ) then + if ( .not. uvflag) then + endif + klev=ilev(k) + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_v(i,j,klev)=sub(ij,k) + enddo + enddo + elseif ( iflag(k) == 8 ) then + klev=ilev(k) + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_q(i,j,klev)=sub(ij,k) + enddo + enddo + elseif ( iflag(k) == 9 ) then + klev=ilev(k) + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_oz(i,j,klev)=sub(ij,k) + enddo + enddo + elseif ( iflag(k) == 10 ) then + klev=ilev(k) + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_ql(i,j,klev)=sub(ij,k) + enddo + enddo + elseif ( iflag(k) == 11 ) then + klev=ilev(k) + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_qi(i,j,klev)=sub(ij,k) + enddo + enddo + elseif ( iflag(k) == 12 ) then + klev=ilev(k) + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_qr(i,j,klev)=sub(ij,k) + enddo + enddo + elseif ( iflag(k) == 13 ) then + klev=ilev(k) + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_qs(i,j,klev)=sub(ij,k) + enddo + enddo + elseif ( iflag(k) == 14 ) then + klev=ilev(k) + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_qg(i,j,klev)=sub(ij,k) + enddo + enddo + elseif ( iflag(k) == 15 .and. present(g_cf) ) then + klev=ilev(k) + ij=0 + do j=1,grd%lon2 + do i=1,grd%lat2 + ij=ij+1 + g_cf(i,j,klev)=sub(ij,k) + enddo + enddo + endif + enddo ! do k=1,icount + + icount=0 + ilev=0 + iflag=0 + + return + +end subroutine general_reload2 end module gfsreadmod subroutine general_read_gfsatm(grd,sp_a,sp_b,filename,uvflag,vordivflag,zflag, & @@ -2348,7 +2570,871 @@ subroutine general_read_gfsatm_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & return end subroutine general_read_gfsatm_nc +subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,zflag, & + gfs_bundle,iret_read) +!$$$ subprogram documentation block +! . . . . +! subprogram: general_read_gfsatm_nc adaptation of read_gfsatm for general resolutions +! prgmmr: martin org: NCEP/EMC date: 2019-09-25 +! +! abstract: copied from read_gfsatm_nems +! +! program history log: +! 2019-09-25 martin - original; copied from general_read_gfsatm_nems +! +! input argument list: +! grd - structure variable containing information about grid +! (initialized by general_sub2grid_create_info, located in general_sub2grid_mod.f90) +! sp_a - structure variable containing spectral information for analysis +! (initialized by general_init_spec_vars, located in general_specmod.f90) +! sp_b - structure variable containing spectral information for input +! fields +! (initialized by general_init_spec_vars, located in general_specmod.f90) +! filename - input netcdf file name +! uvflag - logical to use u,v (.true.) or st,vp (.false.) perturbations +! vordivflag - logical to determine if routine should output vorticity and +! divergence +! zflag - logical to determine if surface height field should be output +! +! output argument list: +! gfs_bundle - bundle carrying guess fields +! iret_read - return code, 0 for successful read. +! +! attributes: +! language: f90 +! machine: ibm RS/6000 SP +! +!$$$ + use kinds, only: r_kind,r_single,i_kind + use mpimod, only: mype + use general_sub2grid_mod, only: sub2grid_info + use general_specmod, only: spec_vars + use mpimod, only: npe + use constants, only: zero,one,fv,r0_01 + use egrid2agrid_mod,only: g_egrid2agrid,g_create_egrid2agrid,egrid2agrid_parm,destroy_egrid2agrid + use general_commvars_mod, only: fill2_ns,filluv2_ns + use constants, only: two,pi,half,deg2rad,r60,r3600 + use gsi_bundlemod, only: gsi_bundle + use gsi_bundlemod, only: gsi_bundlegetpointer + use module_fv3gfs_ncio, only: Dataset, Variable, Dimension, open_dataset,& + close_dataset, get_dim, read_vardata,get_idate_from_time_units + use gfsreadmod, only: general_reload2 + use ncepnems_io, only: imp_physics + + implicit none + + ! Declare local parameters + real(r_kind),parameter:: r0_001 = 0.001_r_kind + + ! Declare passed variables + type(sub2grid_info) ,intent(in ) :: grd + type(spec_vars) ,intent(in ) :: sp_a + character(*) ,intent(in ) :: filename + logical ,intent(in ) :: uvflag,zflag,vordivflag + integer(i_kind) ,intent( out) :: iret_read + type(gsi_bundle) ,intent(inout) :: gfs_bundle + + real(r_kind),pointer,dimension(:,:) :: ptr2d + real(r_kind),pointer,dimension(:,:,:) :: ptr3d + real(r_kind),pointer,dimension(:,:) :: g_ps + real(r_kind),pointer,dimension(:,:,:) :: g_vor,g_div,& + g_q,g_oz,g_tv + real(r_kind),pointer,dimension(:,:,:) :: g_ql,g_qi,g_qr,g_qs,g_qg + + real(r_kind),allocatable,dimension(:,:) :: g_z + real(r_kind),allocatable,dimension(:,:,:) :: g_u,g_v + + ! Declare local variables + character(len=120) :: my_name = 'GENERAL_READ_GFSATM_ALLHYDRO_NC' + integer(i_kind):: iret,nlatm2,nlevs,icm,nord_int + integer(i_kind):: i,j,k,icount,kk,kr + integer(i_kind) :: ier,istatus,istatus1,iredundant + integer(i_kind) :: latb, lonb, levs + integer(i_kind),dimension(npe)::ilev,iflag,mype_use + integer(i_kind),dimension(6):: idate + integer(i_kind),dimension(4):: odate + real(r_kind),allocatable,dimension(:) :: fhour + + real(r_kind),allocatable,dimension(:):: spec_div,spec_vor + real(r_kind),allocatable,dimension(:,:) :: grid, grid_v, & + grid_vor, grid_div, grid_b, grid_b2 + real(r_kind),allocatable,dimension(:,:,:) :: grid_c, grid2, grid_c2 + real(r_single),allocatable,dimension(:,:,:) :: rwork3d0, rwork3d1 + real(r_single),allocatable,dimension(:,:) :: rwork2d + real(r_kind),allocatable,dimension(:) :: work, work_v + real(r_kind),allocatable,dimension(:) :: rlats,rlons,clons,slons + real(r_kind),allocatable,dimension(:) :: rlats_tmp,rlons_tmp + + logical :: procuse,diff_res,eqspace + type(egrid2agrid_parm) :: p_high + logical,dimension(1) :: vector + type(Dataset) :: atmges + type(Dimension) :: ncdim + + + + !****************************************************************************** + ! Initialize variables used below + iret_read=0 + iret=0 + nlatm2=grd%nlat-2 + iflag = 0 + ilev = 0 + + nlevs=grd%nsig + mype_use=-1 + icount=0 + procuse=.false. + if ( mype == 0 ) procuse = .true. + do i=1,npe + if ( grd%recvcounts_s(i-1) > 0 ) then + icount = icount+1 + mype_use(icount)=i-1 + if ( i-1 == mype ) procuse=.true. + endif + enddo + icm=icount + allocate( work(grd%itotsub),work_v(grd%itotsub) ) + work=zero + work_v=zero + + if ( procuse ) then + + atmges = open_dataset(filename, paropen=.true.) + + ! get dimension sizes + ncdim = get_dim(atmges, 'grid_xt'); lonb = ncdim%len + ncdim = get_dim(atmges, 'grid_yt'); latb = ncdim%len + ncdim = get_dim(atmges, 'pfull'); levs = ncdim%len + + ! get time information + idate = get_idate_from_time_units(atmges) + odate(1) = idate(4) !hour + odate(2) = idate(2) !month + odate(3) = idate(3) !day + odate(4) = idate(1) !year + call read_vardata(atmges, 'time', fhour) ! might need to change this to attribute later + ! depends on model changes from + ! Jeff Whitaker + fhour = float(nint(fhour)) + + odate(1) = idate(4) !hour + odate(2) = idate(2) !month + odate(3) = idate(3) !day + odate(4) = idate(1) !year + + diff_res=.false. + if ( latb /= nlatm2 ) then + diff_res=.true. + if ( mype == 0 ) write(6, & + '(a,'': different spatial dimension nlatm2 = '',i4,tr1,''latb = '',i4)') & + trim(my_name),nlatm2,latb + !call stop2(101) + endif + if ( lonb /= grd%nlon ) then + diff_res=.true. + if ( mype == 0 ) write(6, & + '(a,'': different spatial dimension nlon = '',i4,tr1,''lonb = '',i4)') & + trim(my_name),grd%nlon,lonb + !call stop2(101) + endif + if ( levs /= grd%nsig ) then + if ( mype == 0 ) write(6, & + '(a,'': inconsistent spatial dimension nsig = '',i4,tr1,''levs = '',i4)') & + trim(my_name),grd%nsig,levs + call stop2(101) + endif + + allocate( spec_vor(sp_a%nc), spec_div(sp_a%nc) ) + allocate( grid(grd%nlon,nlatm2), grid_v(grd%nlon,nlatm2) ) + if ( diff_res ) then + allocate(grid_b(lonb,latb),grid_c(latb+2,lonb,1),grid2(grd%nlat,grd%nlon,1)) + allocate(grid_b2(lonb,latb),grid_c2(latb+2,lonb,1)) + endif + allocate(rwork3d0(lonb,latb,1)) + allocate(rwork3d1(lonb,latb,1)) + allocate(rwork2d(lonb,latb)) + allocate(rlats(latb+2),rlons(lonb),clons(lonb),slons(lonb)) + call read_vardata(atmges, 'grid_xt', rlons_tmp) + call read_vardata(atmges, 'grid_yt', rlats_tmp) + do j=1,latb + rlats(latb+2-j)=deg2rad*rlats_tmp(j) + end do + do j=1,lonb + rlons(j)=deg2rad*rlons_tmp(j) + end do + deallocate(rlats_tmp,rlons_tmp) + rlats(1)=-half*pi + rlats(latb+2)=half*pi + do j=1,lonb + clons(j)=cos(rlons(j)) + slons(j)=sin(rlons(j)) + enddo + + nord_int=4 + eqspace=.false. + call g_create_egrid2agrid(grd%nlat,sp_a%rlats,grd%nlon,sp_a%rlons, & + latb+2,rlats,lonb,rlons,& + nord_int,p_high,.true.,eqspace=eqspace) + deallocate(rlats,rlons) + + endif ! if ( procuse ) + + ! Get pointer to relevant variables (this should be made flexible and general) + iredundant=0 + call gsi_bundlegetpointer(gfs_bundle,'sf',g_div ,ier) + if ( ier == 0 ) iredundant = iredundant + 1 + call gsi_bundlegetpointer(gfs_bundle,'div',g_div ,ier) + if ( ier == 0 ) iredundant = iredundant + 1 + if ( iredundant==2 ) then + if ( mype == 0 ) then + write(6,*) 'general_read_gfsatm_allhydro_nc: ERROR' + write(6,*) 'cannot handle having both sf and div' + write(6,*) 'Aborting ... ' + endif + call stop2(999) + endif + iredundant=0 + call gsi_bundlegetpointer(gfs_bundle,'vp',g_vor ,ier) + if ( ier == 0 ) iredundant = iredundant + 1 + call gsi_bundlegetpointer(gfs_bundle,'vor',g_vor ,ier) + if ( ier == 0 ) iredundant = iredundant + 1 + if ( iredundant==2 ) then + if ( mype == 0 ) then + write(6,*) 'general_read_gfsatm_allhydro_nc: ERROR' + write(6,*) 'cannot handle having both vp and vor' + write(6,*) 'Aborting ... ' + endif + call stop2(999) + endif + iredundant=0 + call gsi_bundlegetpointer(gfs_bundle,'t' ,g_tv ,ier) + if ( ier == 0 ) iredundant = iredundant + 1 + call gsi_bundlegetpointer(gfs_bundle,'tv',g_tv ,ier) + if ( ier == 0 ) iredundant = iredundant + 1 + if ( iredundant==2 ) then + if ( mype == 0 ) then + write(6,*) 'general_read_gfsatm_allhydro_nc: ERROR' + write(6,*) 'cannot handle having both t and tv' + write(6,*) 'Aborting ... ' + endif + call stop2(999) + endif + istatus=0 + call gsi_bundlegetpointer(gfs_bundle,'ps',g_ps ,ier);istatus=istatus+ier + call gsi_bundlegetpointer(gfs_bundle,'q' ,g_q ,ier);istatus=istatus+ier + call gsi_bundlegetpointer(gfs_bundle,'oz',g_oz ,ier);istatus=istatus+ier +! call gsi_bundlegetpointer(gfs_bundle,'cw',g_cwmr,ier);istatus=istatus+ier + call gsi_bundlegetpointer(gfs_bundle,'ql',g_ql ,ier);istatus1=istatus1+ier + call gsi_bundlegetpointer(gfs_bundle,'qi',g_qi ,ier);istatus1=istatus1+ier + call gsi_bundlegetpointer(gfs_bundle,'qr',g_qr ,ier);istatus1=istatus1+ier + call gsi_bundlegetpointer(gfs_bundle,'qs',g_qs ,ier);istatus1=istatus1+ier + call gsi_bundlegetpointer(gfs_bundle,'qg',g_qg ,ier);istatus1=istatus1+ier +! call gsi_bundlegetpointer(gfs_bundle,'cf',g_cf ,ier);istatus1=istatus1+ier + if ( istatus1 /= 0 ) then + if ( mype == 0 ) then + write(6,*) 'general_read_gfsatm_allhydro_nc: ERROR' + write(6,*) 'Missing some of the required hydrometeor fields for imp_physics = ', imp_physics + write(6,*) 'Aborting ... ' + endif + call stop2(999) + endif + + allocate(g_u(grd%lat2,grd%lon2,grd%nsig),g_v(grd%lat2,grd%lon2,grd%nsig)) + allocate(g_z(grd%lat2,grd%lon2)) + + icount=0 + + ! Process guess fields according to type of input file. NCEP_SIGIO files + ! are spectral coefficient files and need to be transformed to the grid. + ! Once on the grid, fields need to be scattered from the full domain to + ! sub-domains. + + ! Only read Terrain when zflag is true. + + if ( zflag ) then + + icount=icount+1 + iflag(icount)=1 + ilev(icount)=1 + + ! Terrain: spectral --> grid transform, scatter to all mpi tasks + if (mype==mype_use(icount)) then + ! read hs + call read_vardata(atmges, 'hgtsfc', rwork2d) + if ( diff_res ) then + grid_b=rwork2d + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork2d + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + endif + + icount=icount+1 + iflag(icount)=2 + ilev(icount)=1 + + ! Surface pressure: same procedure as terrain + if (mype==mype_use(icount)) then + ! read ps + call read_vardata(atmges, 'pressfc', rwork2d) + rwork2d = r0_001*rwork2d ! convert Pa to cb + if ( diff_res ) then + vector(1)=.false. + grid_b=rwork2d + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork2d + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + + ! Thermodynamic variable: s-->g transform, communicate to all tasks + ! For multilevel fields, each task handles a given level. Periodic + ! mpi_alltoallv calls communicate the grids to all mpi tasks. + ! Finally, the grids are loaded into guess arrays used later in the + ! code. + + do k=1,nlevs + + icount=icount+1 + iflag(icount)=3 + ilev(icount)=k + kr = levs+1-k ! netcdf is top to bottom, need to flip + + if (mype==mype_use(icount)) then + call read_vardata(atmges, 'spfh', rwork3d1, nslice=kr, slicedim=3) + call read_vardata(atmges, 'tmp', rwork3d0, nslice=kr, slicedim=3) + rwork2d = rwork3d0(:,:,1) * (one+fv*rwork3d1(:,:,1)) + if ( diff_res ) then + grid_b=rwork2d + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork2d + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + end do + + if ( vordivflag .or. .not. uvflag ) then + do k=1,nlevs + kr = levs+1-k ! netcdf is top to bottom, need to flip + icount=icount+1 + iflag(icount)=4 + ilev(icount)=k + + if (mype==mype_use(icount)) then + call read_vardata(atmges, 'ugrd', rwork3d0, nslice=kr, slicedim=3) + call read_vardata(atmges, 'vgrd', rwork3d1, nslice=kr, slicedim=3) + ! Vorticity + ! Convert grid u,v to div and vor + if ( diff_res ) then + grid_b = rwork3d0(:,:,1) + grid_b2 = rwork3d1(:,:,1) + vector(1)=.true. + call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + do j=1,grd%nlon + do i=2,grd%nlat-1 + grid(j,grd%nlat-i)=grid2(i,j,1) + enddo + enddo + call g_egrid2agrid(p_high,grid_c2,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work_v(kk)=grid2(i,j,1) + enddo + do j=1,grd%nlon + do i=2,grd%nlat-1 + grid_v(j,grd%nlat-i)=grid2(i,j,1) + enddo + enddo + else + grid = rwork3d0(:,:,1) + grid_v = rwork3d1(:,:,1) + call general_filluv_ns(grd,slons,clons,grid,grid_v,work,work_v) + endif + allocate( grid_vor(grd%nlon,nlatm2)) + call general_sptez_v(sp_a,spec_div,spec_vor,grid,grid_v,-1) + call general_sptez_s_b(sp_a,sp_a,spec_vor,grid_vor,1) + ! Load values into rows for south and north pole + call general_fill_ns(grd,grid_vor,work) + deallocate(grid_vor) + endif + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + end do + do k=1,nlevs + kr = levs+1-k ! netcdf is top to bottom, need to flip + + icount=icount+1 + iflag(icount)=5 + ilev(icount)=k + + if (mype==mype_use(icount)) then + call read_vardata(atmges, 'ugrd', rwork3d0, nslice=kr, slicedim=3) + call read_vardata(atmges, 'vgrd', rwork3d1, nslice=kr, slicedim=3) + ! Divergence + ! Convert grid u,v to div and vor + if ( diff_res ) then + grid_b = rwork3d0(:,:,1) + grid_b2 = rwork3d1(:,:,1) + vector(1)=.true. + call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + do j=1,grd%nlon + do i=2,grd%nlat-1 + grid(j,grd%nlat-i)=grid2(i,j,1) + enddo + enddo + call g_egrid2agrid(p_high,grid_c2,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work_v(kk)=grid2(i,j,1) + enddo + do j=1,grd%nlon + do i=2,grd%nlat-1 + grid_v(j,grd%nlat-i)=grid2(i,j,1) + enddo + enddo + else + grid = rwork3d0(:,:,1) + grid_v = rwork3d1(:,:,1) + call general_filluv_ns(grd,slons,clons,grid,grid_v,work,work_v) + endif + allocate( grid_div(grd%nlon,nlatm2) ) + call general_sptez_v(sp_a,spec_div,spec_vor,grid,grid_v,-1) + call general_sptez_s_b(sp_a,sp_a,spec_div,grid_div,1) + ! Load values into rows for south and north pole + call general_fill_ns(grd,grid_div,work) + deallocate(grid_div) + endif + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + + end do + endif ! if ( vordivflag .or. .not. uvflag ) + + if ( uvflag ) then + do k=1,nlevs + kr = levs+1-k ! netcdf is top to bottom, need to flip + icount=icount+1 + iflag(icount)=6 + ilev(icount)=k + + if (mype==mype_use(icount)) then + call read_vardata(atmges, 'ugrd', rwork3d0, nslice=kr, slicedim=3) + call read_vardata(atmges, 'vgrd', rwork3d1, nslice=kr, slicedim=3) + + if ( diff_res ) then + grid_b = rwork3d0(:,:,1) + grid_b2 = rwork3d1(:,:,1) + vector(1)=.true. + call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid = rwork3d0(:,:,1) + grid_v = rwork3d1(:,:,1) + call general_filluv_ns(grd,slons,clons,grid,grid_v,work,work_v) + endif + endif + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + + icount=icount+1 + iflag(icount)=7 + ilev(icount)=k + + if (mype==mype_use(icount)) then + ! V + call read_vardata(atmges, 'ugrd', rwork3d0, nslice=kr, slicedim=3) + call read_vardata(atmges, 'vgrd', rwork3d1, nslice=kr, slicedim=3) + if ( diff_res ) then + grid_b = rwork3d0(:,:,1) + grid_b2 = rwork3d1(:,:,1) + vector(1)=.true. + call filluv2_ns(grid_b,grid_b2,grid_c(:,:,1),grid_c2(:,:,1),latb+2,lonb,slons,clons) + call g_egrid2agrid(p_high,grid_c2,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid = rwork3d0(:,:,1) + grid_v = rwork3d1(:,:,1) + ! Note work_v and work are switched because output must be in work. + call general_filluv_ns(grd,slons,clons,grid,grid_v,work_v,work) + endif + endif + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + end do + endif ! if ( uvflag ) + + do k=1,nlevs + kr = levs+1-k ! netcdf is top to bottom, need to flip + icount=icount+1 + iflag(icount)=8 + ilev(icount)=k + + if (mype==mype_use(icount)) then + ! Specific humidity + call read_vardata(atmges, 'spfh', rwork3d0, nslice=kr, slicedim=3) + if ( diff_res ) then + grid_b=rwork3d0(:,:,1) + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid = rwork3d0(:,:,1) + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + end do + + do k=1,nlevs + kr = levs+1-k ! netcdf is top to bottom, need to flip + + icount=icount+1 + iflag(icount)=9 + ilev(icount)=k + + if (mype==mype_use(icount)) then + call read_vardata(atmges, 'o3mr', rwork3d0, nslice=kr, slicedim=3) + ! Ozone mixing ratio + if ( diff_res ) then + grid_b=rwork3d0(:,:,1) + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork3d0(:,:,1) + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + end do + + do k=1,nlevs + icount=icount+1 + iflag(icount)=10 + ilev(icount)=k + kr = levs+1-k ! netcdf is top to bottom, need to flip + + if (mype==mype_use(icount)) then + call read_vardata(atmges, 'clwmr', rwork3d0, nslice=kr, slicedim=3) + ! Cloud liquid water mixing ratio. + if ( diff_res ) then + grid_b=rwork3d0(:,:,1) + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork3d0(:,:,1) + call general_fill_ns(grd,grid,work) + endif + endif + + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + enddo ! do k=1,nlevs + + do k=1,nlevs + icount=icount+1 + iflag(icount)=11 + ilev(icount)=k + kr = levs+1-k ! netcdf is top to bottom, need to flip + + if (mype==mype_use(icount)) then + call read_vardata(atmges, 'icmr', rwork3d0, nslice=kr, slicedim=3) + ! Cloud ice water mixing ratio. + if ( diff_res ) then + grid_b=rwork3d0(:,:,1) + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork3d0(:,:,1) + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + enddo ! do k=1,nlevs + + do k=1,nlevs + icount=icount+1 + iflag(icount)=12 + ilev(icount)=k + kr = levs+1-k ! netcdf is top to bottom, need to flip + + if (mype==mype_use(icount)) then + call read_vardata(atmges, 'rwmr', rwork3d0, nslice=kr, slicedim=3) + ! Rain water mixing ratio. + if ( diff_res ) then + grid_b=rwork3d0(:,:,1) + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork3d0(:,:,1) + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + enddo ! do k=1,nlevs + + do k=1,nlevs + icount=icount+1 + iflag(icount)=13 + ilev(icount)=k + kr = levs+1-k ! netcdf is top to bottom, need to flip + + if (mype==mype_use(icount)) then + call read_vardata(atmges, 'snmr', rwork3d0, nslice=kr, slicedim=3) + ! Snow water mixing ratio. + if ( diff_res ) then + grid_b=rwork3d0(:,:,1) + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork3d0(:,:,1) + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm ) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + enddo ! do k=1,nlevs + + do k=1,nlevs + icount=icount+1 + iflag(icount)=14 + ilev(icount)=k + kr = levs+1-k ! netcdf is top to bottom, need to flip + + if (mype==mype_use(icount)) then + call read_vardata(atmges, 'grle', rwork3d0, nslice=kr, slicedim=3) + ! Graupel mixing ratio. + if ( diff_res ) then + grid_b=rwork3d0(:,:,1) + vector(1)=.false. + call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) + call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) + do kk=1,grd%itotsub + i=grd%ltosi_s(kk) + j=grd%ltosj_s(kk) + work(kk)=grid2(i,j,1) + enddo + else + grid=rwork3d0(:,:,1) + call general_fill_ns(grd,grid,work) + endif + endif + if ( icount == icm .or. k==nlevs) then + call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & + g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag) + endif + enddo ! do k=1,nlevs + +! do k=1,nlevs +! icount=icount+1 +! iflag(icount)=15 +! ilev(icount)=k +! kr = levs+1-k ! netcdf is top to bottom, need to flip +! +! if (mype==mype_use(icount)) then +! call read_vardata(atmges, 'cld_amt', rwork3d0, nslice=kr, slicedim=3) +! ! Cloud amount (cloud fraction). +! if ( diff_res ) then +! grid_b=rwork3d0(:,:,1) +! vector(1)=.false. +! call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb) +! call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector) +! do kk=1,grd%itotsub +! i=grd%ltosi_s(kk) +! j=grd%ltosj_s(kk) +! work(kk)=grid2(i,j,1) +! enddo +! else +! grid=rwork3d0(:,:,1) +! call general_fill_ns(grd,grid,work) +! endif +! +! endif +! if ( icount == icm .or. k==nlevs ) then +! call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, & +! g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_cf) +! endif +! enddo ! do k=1,nlevs + + if ( procuse ) then + if ( diff_res) deallocate(grid_b,grid_b2,grid_c,grid_c2,grid2) + call destroy_egrid2agrid(p_high) + deallocate(spec_div,spec_vor) + deallocate(rwork3d1,rwork3d0,clons,slons) + deallocate(rwork2d) + deallocate(grid,grid_v) + call close_dataset(atmges) + endif + deallocate(work, work_v) + + ! Convert dry temperature to virtual temperature + !do k=1,grd%nsig + ! do j=1,grd%lon2 + ! do i=1,grd%lat2 + ! g_tv(i,j,k) = g_tv(i,j,k)*(one+fv*g_q(i,j,k)) + ! enddo + ! enddo + !enddo + + ! Load u->div and v->vor slot when uv are used instead + if ( uvflag ) then + call gsi_bundlegetpointer(gfs_bundle,'u' ,ptr3d,ier) + if ( ier == 0 ) then + ptr3d=g_u + call gsi_bundlegetpointer(gfs_bundle,'v' ,ptr3d,ier) + if ( ier == 0 ) ptr3d=g_v + else ! in this case, overload: return u/v in sf/vp slot + call gsi_bundlegetpointer(gfs_bundle,'sf' ,ptr3d,ier) + if ( ier == 0 ) then + ptr3d=g_u + call gsi_bundlegetpointer(gfs_bundle,'vp' ,ptr3d,ier) + if ( ier == 0 ) ptr3d=g_v + endif + endif + else ! in this case, overload: return u/v in sf/vp slot + call gsi_bundlegetpointer(gfs_bundle,'sf' ,ptr3d,ier) + if ( ier == 0 ) ptr3d=g_u + call gsi_bundlegetpointer(gfs_bundle,'vp' ,ptr3d,ier) + if ( ier == 0 ) ptr3d=g_v + endif + if (zflag) then + call gsi_bundlegetpointer(gfs_bundle,'z' ,ptr2d,ier) + if ( ier == 0 ) ptr2d=g_z + endif + + ! Clean up + deallocate(g_z) + deallocate(g_u,g_v) + + ! Print date/time stamp + if ( mype == 0 ) then + write(6,700) lonb,latb,nlevs,grd%nlon,nlatm2,& + fhour,odate,trim(filename) +700 format('GENERAL_READ_GFSATM_NC: read lonb,latb,levs=',& + 3i6,', scatter nlon,nlat=',2i6,', hour=',f6.1,', idate=',4i5,1x,a) + endif + + return +end subroutine general_read_gfsatm_allhydro_nc subroutine general_fill_ns(grd,grid_in,grid_out) ! !USES: diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index 7177719fde..bdf340801f 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -19,7 +19,7 @@ module gsimod dtbduv_on,time_window_max,offtime_data,init_directories,oberror_tune,ext_sonde, & blacklst,init_obsmod_vars,lobsdiagsave,lobskeep,lobserver,hilbert_curve,& lread_obs_save,lread_obs_skip,time_window_rad,tcp_posmatch,tcp_box, & - neutral_stability_windfact_2dvar,use_similarity_2dvar + neutral_stability_windfact_2dvar,use_similarity_2dvar,ta2tb use gsi_dbzOper, only: diag_radardbz use obsmod, only: doradaroneob,oneoblat,oneoblon,oneobheight,oneobvalue,oneobddiff,oneobradid,& @@ -699,6 +699,8 @@ module gsimod ! efsoi_order - sets order of EFSOI-like calculation ! lupdqc - logical to replace the obs errors from satinfo with diag of est(R) in the case of correlated obs ! lqcoef - logical to combine the inflation coefficients generated by qc with est(R) +! ta2tb - logical to use brightness temperature (SDR) instead of antenna +! temperature (TDR) for assimilation ! l_use_rw_columntilt - option to assimilate radar column-tilt radial wind obs in GSI ! (.TRUE.: on; .FALSE.: off) / Inputfile: l2rwbufr_cltl (bufr format) ! l_use_dbz_directDA - option to assimilate radar reflectivity obs directly in GSI @@ -750,7 +752,7 @@ module gsimod write_fv3_incr,incvars_to_zero,incvars_zero_strat,incvars_efold,diag_version,& cao_check,lcalc_gfdl_cfrac,tau_fcst,efsoi_order,lupdqc,lqcoef,cnvw_option,l2rwthin,hurricane_radar,& l_reg_update_hydro_delz, l_obsprvdiag,& - l_use_dbz_directDA, l_use_rw_columntilt + l_use_dbz_directDA, l_use_rw_columntilt, ta2tb ! GRIDOPTS (grid setup variables,including regional specific variables): ! jcap - spectral resolution diff --git a/src/gsi/netcdfgfs_io.f90 b/src/gsi/netcdfgfs_io.f90 index dd3f991da5..5ebf7296bd 100644 --- a/src/gsi/netcdfgfs_io.f90 +++ b/src/gsi/netcdfgfs_io.f90 @@ -116,8 +116,9 @@ subroutine read_ !$$$ end documentation block use kinds, only: i_kind,r_kind + use constants, only: qcmin use gridmod, only: sp_a,grd_a,lat2,lon2,nsig - use guess_grids, only: ifilesig,nfldsig + use guess_grids, only: ifilesig,nfldsig,ntguessig use gsi_metguess_mod, only: gsi_metguess_bundle use gsi_bundlemod, only: gsi_bundlegetpointer use gsi_bundlemod, only: gsi_bundlecreate @@ -128,12 +129,15 @@ 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 gridmod, only: fv3_full_hydro + implicit none character(len=*),parameter::myname_=myname//'*read_' character(24) filename integer(i_kind):: it, istatus, inner_vars, num_fields - integer(i_kind):: iret_ql,iret_qi + integer(i_kind):: i,j,k + real(r_kind),pointer,dimension(:,: ):: ges_ps_it =>NULL() real(r_kind),pointer,dimension(:,: ):: ges_z_it =>NULL() @@ -147,6 +151,10 @@ subroutine read_ real(r_kind),pointer,dimension(:,:,:):: ges_cwmr_it=>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_ql_it => NULL() real(r_kind),pointer,dimension(:,:,:):: ges_qi_it => NULL() + real(r_kind),pointer,dimension(:,:,:):: ges_qr_it => NULL() + real(r_kind),pointer,dimension(:,:,:):: ges_qs_it => NULL() + real(r_kind),pointer,dimension(:,:,:):: ges_qg_it => NULL() + real(r_kind),pointer,dimension(:,:,:):: ges_cf_it => NULL() type(sub2grid_info) :: grd_t logical regional @@ -155,24 +163,33 @@ subroutine read_ type(gsi_bundle) :: atm_bundle type(gsi_grid) :: atm_grid integer(i_kind),parameter :: n2d=2 - integer(i_kind),parameter :: n3d=8 + ! integer(i_kind),parameter :: n3d=8 + integer(i_kind),parameter :: n3d=14 character(len=4), parameter :: vars2d(n2d) = (/ 'z ', 'ps ' /) + ! character(len=4), parameter :: vars3d(n3d) = (/ 'u ', 'v ', & + ! 'vor ', 'div ', & + ! 'tv ', 'q ', & + ! 'cw ', 'oz ' /) character(len=4), parameter :: vars3d(n3d) = (/ 'u ', 'v ', & 'vor ', 'div ', & 'tv ', 'q ', & - 'cw ', 'oz ' /) + 'cw ', 'oz ', & + 'ql ', 'qi ', & + 'qr ', 'qs ', & + 'qg ', 'cf ' /) real(r_kind),pointer,dimension(:,:):: ptr2d =>NULL() real(r_kind),pointer,dimension(:,:,:):: ptr3d =>NULL() regional=.false. inner_vars=1 - num_fields=min(8*grd_a%nsig+2,npe) + num_fields=min(14*grd_a%nsig+2,npe) ! Create temporary communication information fore read routines call general_sub2grid_create_info(grd_t,inner_vars,grd_a%nlat,grd_a%nlon, & grd_a%nsig,num_fields,regional) ! Allocate bundle used for reading members call gsi_gridcreate(atm_grid,lat2,lon2,nsig) + call gsi_bundlecreate(atm_bundle,atm_grid,'aux-atm-read',istatus,names2d=vars2d,names3d=vars3d) if(istatus/=0) then write(6,*) myname_,': trouble creating atm_bundle' @@ -184,8 +201,17 @@ subroutine read_ write(filename,'(''sigf'',i2.2)') ifilesig(it) ! Read background fields into bundle - call general_read_gfsatm_nc(grd_t,sp_a,filename,.true.,.true.,.true.,& - atm_bundle,.true.,istatus) + if (fv3_full_hydro) then + if (mype==0) write(*,*) 'calling general_read_gfsatm_allhydro_nc', it + call general_read_gfsatm_allhydro_nc(grd_t,sp_a,filename,.true.,.true.,.true.,& + atm_bundle,istatus) ! this loads cloud and precip + if (mype==0) write(*,*) 'done with general_read_gfsatm_allhydro_nc', it + else + if (mype==0) write(*,*) 'calling general_read_gfsatm_nc' + call general_read_gfsatm_nc(grd_t,sp_a,filename,.true.,.true.,.true.,& + atm_bundle,.true.,istatus) + if (mype==0) write(*,*) 'done with general_read_gfsatm_nc' + end if inithead=.false. zflag=.false. @@ -193,17 +219,41 @@ subroutine read_ ! Set values to actual MetGuess fields call set_guess_ - l_cld_derived = associated(ges_cwmr_it).and.& - associated(ges_q_it) .and.& - associated(ges_ql_it) .and.& - associated(ges_qi_it) .and.& - associated(ges_tv_it) -! call set_cloud_lower_bound(ges_cwmr_it) - if (mype==0) write(6,*)'READ_GFS_NETCDF: l_cld_derived = ', l_cld_derived - - if (l_cld_derived) then - call cloud_calc_gfs(ges_ql_it,ges_qi_it,ges_cwmr_it,ges_q_it,ges_tv_it,.true.) + if (it==ntguessig) then + if (mype==0) write(6,*)'Print guess field ... after set_guess' + call prt_guess('guess') + endif + if (fv3_full_hydro) then + do k=1, nsig + do j=1, lon2 + do i=1, lat2 + ! set lower bound to hydrometeors + if (associated(ges_ql_it)) ges_ql_it(i,j,k) = max(qcmin,ges_ql_it(i,j,k)) + if (associated(ges_qi_it)) ges_qi_it(i,j,k) = max(qcmin,ges_qi_it(i,j,k)) + if (associated(ges_qr_it)) ges_qr_it(i,j,k) = max(qcmin,ges_qr_it(i,j,k)) + if (associated(ges_qs_it)) ges_qs_it(i,j,k) = max(qcmin,ges_qs_it(i,j,k)) + if (associated(ges_qg_it)) ges_qg_it(i,j,k) = max(qcmin,ges_qg_it(i,j,k)) + if (associated(ges_cf_it)) ges_cf_it(i,j,k) = min(max(zero,ges_cf_it(i,j,k)),one) + enddo + enddo + enddo + else + l_cld_derived = associated(ges_cwmr_it).and.& + associated(ges_q_it) .and.& + associated(ges_ql_it) .and.& + associated(ges_qi_it) .and.& + associated(ges_tv_it) +! call set_cloud_lower_bound(ges_cwmr_it) + if (mype==0) write(6,*)'READ_GFS_NETCDF: l_cld_derived = ', l_cld_derived + + if (l_cld_derived) then + call cloud_calc_gfs(ges_ql_it,ges_qi_it,ges_cwmr_it,ges_q_it,ges_tv_it,.true.) + end if end if + if (it==ntguessig) then + if (mype==0) write(6,*)'Print guess field ... after set lower bound for hydrometers' + call prt_guess('guess') + endif end do call general_sub2grid_destroy_info(grd_t) @@ -263,15 +313,36 @@ subroutine set_guess_ call gsi_bundlegetpointer (gsi_metguess_bundle(it),'cw',ges_cwmr_it,istatus) if(istatus==0) ges_cwmr_it = ptr3d endif - call gsi_bundlegetpointer (gsi_metguess_bundle(it),'ql',ges_ql_it, iret_ql) - call gsi_bundlegetpointer (gsi_metguess_bundle(it),'qi',ges_qi_it, iret_qi) - if (iret_ql/=0) then - if (mype==0) write(6,*)'READ_ NETCDF: cannot get pointer to ql,iret_ql=',iret_ql + call gsi_bundlegetpointer (atm_bundle,'ql',ptr3d,istatus) + if (istatus==0) then + call gsi_bundlegetpointer (gsi_metguess_bundle(it),'ql',ges_ql_it,istatus) + if(istatus==0) ges_ql_it = ptr3d + endif + call gsi_bundlegetpointer (atm_bundle,'qi',ptr3d,istatus) + if (istatus==0) then + call gsi_bundlegetpointer (gsi_metguess_bundle(it),'qi',ges_qi_it,istatus) + if(istatus==0) ges_qi_it = ptr3d endif - if (iret_qi/=0) then - if (mype==0) write(6,*)'READ_ NETCDF: cannot get pointer to qi,iret_qi=',iret_qi + call gsi_bundlegetpointer (atm_bundle,'qr',ptr3d,istatus) + if (istatus==0) then + call gsi_bundlegetpointer (gsi_metguess_bundle(it),'qr',ges_qr_it,istatus) + if(istatus==0) ges_qr_it = ptr3d + endif + call gsi_bundlegetpointer (atm_bundle,'qs',ptr3d,istatus) + if (istatus==0) then + call gsi_bundlegetpointer (gsi_metguess_bundle(it),'qs',ges_qs_it,istatus) + if(istatus==0) ges_qs_it = ptr3d + endif + call gsi_bundlegetpointer (atm_bundle,'qg',ptr3d,istatus) + if (istatus==0) then + call gsi_bundlegetpointer (gsi_metguess_bundle(it),'qg',ges_qg_it,istatus) + if(istatus==0) ges_qg_it = ptr3d + endif + call gsi_bundlegetpointer (atm_bundle,'cf',ptr3d,istatus) + if (istatus==0) then + call gsi_bundlegetpointer (gsi_metguess_bundle(it),'cf',ges_cf_it,istatus) + if(istatus==0) ges_cf_it = ptr3d endif - end subroutine set_guess_ end subroutine read_ diff --git a/src/gsi/obsmod.F90 b/src/gsi/obsmod.F90 index 40dc3cc0dc..1fb03d0940 100644 --- a/src/gsi/obsmod.F90 +++ b/src/gsi/obsmod.F90 @@ -444,6 +444,7 @@ module obsmod public :: use_similarity_2dvar public :: time_window_rad public :: perturb_fact,dtbduv_on,nsat1,obs_sub_comm,mype_diaghdr + public :: ta2tb public :: lobsdiag_allocated public :: nloz_v8,nloz_v6,nloz_omi,nlco,nobskeep public :: rmiss_single,nchan_total,mype_sst,mype_gps @@ -614,6 +615,7 @@ module obsmod integer(i_kind) ntilt_radarfiles,tcp_posmatch,tcp_box + logical :: ta2tb logical :: doradaroneob logical :: vr_dealisingopt, if_vterminal, if_model_dbz, inflate_obserr, if_vrobs_raw, l2rwthin character(4) :: whichradar,oneobradid @@ -795,6 +797,8 @@ subroutine init_obsmod_dflts nobskeep=0 lsaveobsens=.false. l_do_adjoint=.true. ! .true. = apply H^T when in int routines + ta2tb=.false. ! .true. = assimilation antenna temperature for + ! AMSU-A, ATMS and MHS oberrflg = .false. bflag = .false. ! sfcmodel = .false. ! .false. = do not use boundary layer model diff --git a/src/gsi/prt_guess.f90 b/src/gsi/prt_guess.f90 index e1e251195a..97c1e0a7c1 100644 --- a/src/gsi/prt_guess.f90 +++ b/src/gsi/prt_guess.f90 @@ -374,7 +374,7 @@ subroutine prt_guess2(sgrep) ier=ier+istatus nvarsc=nvarsc+1 endif - if (ier/=0) return ! this is a fundamental routine, when some not found just +! if (ier/=0) return ! this is a fundamental routine, when some not found just ! emily nvars = nvars1+nvarsc+nvars3 nvars2 = nvars1+nvarsc diff --git a/src/gsi/qcmod.f90 b/src/gsi/qcmod.f90 index 525da81979..824226e133 100644 --- a/src/gsi/qcmod.f90 +++ b/src/gsi/qcmod.f90 @@ -2807,7 +2807,7 @@ subroutine qc_avhrr(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse, & end subroutine qc_avhrr subroutine qc_amsua(nchanl,is,ndat,nsig,npred,sea,land,ice,snow,mixed,luse, & - zsges,cenlat,tb_obsbc1,si_mean,cosza,clw,tbc,ptau5,emissivity_k,ts, & + zsges,cenlat,tb_obsbc1,cosza,clw,tbc,ptau5,emissivity_k,ts, & pred,predchan,id_qc,aivals,errf,errf0,clwp_amsua,varinv,cldeff_obs,cldeff_fg,factch6, & cld_rbc_idx,sfc_speed,error0,clw_guess_retrieval,scatp,radmod) @@ -2890,7 +2890,7 @@ subroutine qc_amsua(nchanl,is,ndat,nsig,npred,sea,land,ice,snow,mixed,luse, & logical, intent(in ) :: sea,land,ice,snow,mixed,luse integer(i_kind), intent(in ) :: ndat,nsig,npred,nchanl,is integer(i_kind),dimension(nchanl), intent(inout) :: id_qc - real(r_kind), intent(in ) :: zsges,cenlat,tb_obsbc1,si_mean + real(r_kind), intent(in ) :: zsges,cenlat,tb_obsbc1 real(r_kind),dimension(nchanl), intent(in ) :: cldeff_obs,cldeff_fg real(r_kind), intent(in ) :: cosza,clw,clwp_amsua,clw_guess_retrieval real(r_kind), intent(in ) :: sfc_speed,scatp @@ -3121,9 +3121,7 @@ subroutine qc_amsua(nchanl,is,ndat,nsig,npred,sea,land,ice,snow,mixed,luse, & enddo endif else if (latms) then - if (si_mean >= 20.0_r_kind) then - efactmc=zero - vfactmc=zero + if (abs(cldeff_obs(16)-cldeff_obs(17))>10.0_r_kind) then if(id_qc(ich890) == igood_qc)id_qc(ich890)=ifail_factch1617_qc errf(ich890) = zero varinv(ich890) = zero @@ -3132,12 +3130,34 @@ subroutine qc_amsua(nchanl,is,ndat,nsig,npred,sea,land,ice,snow,mixed,luse, & errf(i) = zero varinv(i) = zero enddo - errf(1:ich544)=zero - varinv(1:ich544)=zero - do i=1,ich544 - if(id_qc(i) == igood_qc)id_qc(i)=ifail_factch1617_qc - end do + if (abs(cldeff_obs(16)-cldeff_obs(17))>15.0_r_kind) then + efactmc=zero + vfactmc=zero + errf(1:ich544)=zero + varinv(1:ich544)=zero + do i=1,ich544 + if(id_qc(i) == igood_qc)id_qc(i)=ifail_factch1617_qc + end do + end if end if + ! test in the future + ! if (si_mean >= 20.0_r_kind) then + ! efactmc=zero + ! vfactmc=zero + ! if(id_qc(ich890) == igood_qc)id_qc(ich890)=ifail_factch1617_qc + ! errf(ich890) = zero + ! varinv(ich890) = zero + ! do i=17,22 ! AMSU-B/MHS like channels + ! if(id_qc(i) == igood_qc)id_qc(i)=ifail_factch1617_qc + ! errf(i) = zero + ! varinv(i) = zero + ! enddo + ! errf(1:ich544)=zero + ! varinv(1:ich544)=zero + ! do i=1,ich544 + ! if(id_qc(i) == igood_qc)id_qc(i)=ifail_factch1617_qc + ! end do + ! end if else ! QC based on the sensitivity of Tb to the surface emissivity ! de1,de2,de3,de15 become smaller as the observation is more cloudy -- ! i.e., less affected by the surface emissivity quality control check @@ -3594,7 +3614,7 @@ subroutine qc_mhs(nchanl,ndat,nsig,is,sea,land,ice,snow,mhs,luse, & end subroutine qc_mhs subroutine qc_atms(nchanl,is,ndat,nsig,npred,sea,land,ice,snow,mixed,luse, & - zsges,cenlat,tb_obsbc1,si_mean,cosza,clw,tbc,ptau5,emissivity_k,ts, & + zsges,cenlat,tb_obsbc1,cosza,clw,tbc,ptau5,emissivity_k,ts, & pred,predchan,id_qc,aivals,errf,errf0,clwp_amsua,varinv,cldeff_obs,cldeff_fg,factch6, & cld_rbc_idx,sfc_speed,error0,clw_guess_retrieval,scatp,radmod) @@ -3667,7 +3687,7 @@ subroutine qc_atms(nchanl,is,ndat,nsig,npred,sea,land,ice,snow,mixed,luse, & logical, intent(in ) :: sea,land,ice,snow,mixed,luse integer(i_kind), intent(in ) :: nchanl,is,ndat,nsig,npred integer(i_kind),dimension(nchanl), intent(inout) :: id_qc - real(r_kind), intent(in ) :: zsges,cenlat,tb_obsbc1,si_mean + real(r_kind), intent(in ) :: zsges,cenlat,tb_obsbc1 real(r_kind),dimension(nchanl), intent(in ) :: cldeff_obs, cldeff_fg real(r_kind), intent(in ) :: cosza,clw,clwp_amsua,clw_guess_retrieval real(r_kind), intent(in ) :: sfc_speed,scatp @@ -3683,7 +3703,7 @@ subroutine qc_atms(nchanl,is,ndat,nsig,npred,sea,land,ice,snow,mixed,luse, & ! For now, just pass all channels to qc_amsua call qc_amsua (nchanl,is,ndat,nsig,npred,sea,land,ice,snow,mixed,luse, & - zsges,cenlat,tb_obsbc1,si_mean,cosza,clw,tbc,ptau5,emissivity_k,ts, & + zsges,cenlat,tb_obsbc1,cosza,clw,tbc,ptau5,emissivity_k,ts, & pred,predchan,id_qc,aivals,errf,errf0,clwp_amsua,varinv,cldeff_obs,cldeff_fg,factch6, & cld_rbc_idx,sfc_speed,error0,clw_guess_retrieval,scatp,radmod) diff --git a/src/gsi/read_atms.f90 b/src/gsi/read_atms.f90 index 6fc0006c93..9f5efb5301 100644 --- a/src/gsi/read_atms.f90 +++ b/src/gsi/read_atms.f90 @@ -76,7 +76,7 @@ subroutine read_atms(mype,val_tovs,ithin,isfcalc,& use satthin, only: super_val,itxmax,makegrids,destroygrids,checkob, & finalcheck,map2tgrid,score_crit use satthin, only: radthin_time_info,tdiff2crit - use obsmod, only: time_window_max + use obsmod, only: time_window_max, ta2tb use radinfo, only: iuse_rad,newchn,cbias,nusis,jpch_rad,air_rad,ang_rad, & use_edges,radedge1,radedge2,nusis,radstart,radstep,newpc4pred,maxscan use radinfo, only: adp_anglebc @@ -482,7 +482,11 @@ subroutine read_atms(mype,val_tovs,ithin,isfcalc,& ! TMBR is actually the antenna temperature for most microwave sounders but for ! ATMS it is stored in TMANT. ! ATMS is assumed not to come via EARS - call ufbrep(lnbufr,data1b8,1,nchanl,iret,'TMANT') + if (ta2tb) then + call ufbrep(lnbufr,data1b8,1,nchanl,iret,'TMBR') + else + call ufbrep(lnbufr,data1b8,1,nchanl,iret,'TMANT') + endif bt_save(1:nchanl,iob) = data1b8(1:nchanl) diff --git a/src/gsi/read_bufrtovs.f90 b/src/gsi/read_bufrtovs.f90 index cf3e7ea17c..4027384bf3 100644 --- a/src/gsi/read_bufrtovs.f90 +++ b/src/gsi/read_bufrtovs.f90 @@ -91,6 +91,8 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,& ! channels are missing. ! 2018-04-19 eliu - allow data selection for precipitation-affected data ! 2018-05-21 j.jin - added time-thinning. Moved the checking of thin4d into satthin.F90. +! 2020-04021 s.sieron - change converting brightness temperatures to antenna temperautres to the +! other way around. added support for multiple SpcCoeff files and ACCoeffs ! ! input argument list: ! mype - mpi task id @@ -131,7 +133,7 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,& use satthin, only: super_val,itxmax,makegrids,destroygrids,checkob, & finalcheck,map2tgrid,score_crit use satthin, only: radthin_time_info,tdiff2crit - use obsmod, only: time_window_max + use obsmod, only: time_window_max, ta2tb use radinfo, only: iuse_rad,newchn,cbias,predx,nusis,jpch_rad,air_rad,ang_rad, & use_edges,radedge1, radedge2, radstart,radstep,newpc4pred use radinfo, only: crtm_coeffs_path,adp_anglebc @@ -141,9 +143,10 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,& crtm_kind => fp, & MAX_SENSOR_ZENITH_ANGLE use crtm_spccoeff, only: sc,crtm_spccoeff_load,crtm_spccoeff_destroy + use ACCoeff_Define, only: ACCoeff_type use calc_fov_crosstrk, only : instrument_init, fov_cleanup, fov_check use gsi_4dvar, only: l4dvar,l4densvar,iwinbgn,winlen - use antcorr_application, only: remove_antcorr + use antcorr_application, only: remove_antcorr, apply_antcorr use mpeu_util, only: getindex use deter_sfc_mod, only: deter_sfc_fov,deter_sfc use gsi_nstcouplermod, only: nst_gsi,nstinfo @@ -174,7 +177,8 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,& ! Declare local parameters character(8),parameter:: fov_flag="crosstrk" - integer(i_kind),parameter:: n1bhdr=13 + ! change from 13 to 14 for sacv + integer(i_kind),parameter:: n1bhdr=14 integer(i_kind),parameter:: n2bhdr=4 real(r_kind),parameter:: r360=360.0_r_kind real(r_kind),parameter:: tbmin=50.0_r_kind @@ -190,6 +194,7 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,& integer(i_kind) ireadsb,ireadmg,irec,next,nrec_startx integer(i_kind) i,j,k,ifov,ntest,llll + integer(i_kind) sacv integer(i_kind) iret,idate,nchanl,n,idomsfc(1) integer(i_kind) ich1,ich2,ich8,ich15,ich16,ich17 integer(i_kind) kidsat,instrument,maxinfo @@ -234,6 +239,11 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,& logical :: critical_channels_missing,quiet logical :: print_verbose + logical :: spc_coeff_found + integer(i_kind) :: spc_coeff_versions + character(len=80) :: spc_filename + type(ACCoeff_type),dimension(3) :: accoeff_sets + !************************************************************************** ! Initialize variables @@ -476,7 +486,7 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,& if(dval_use) maxinfo=maxinfo+2 nreal = maxinfo + nstinfo nele = nreal + nchanl - hdr1b ='SAID FOVN YEAR MNTH DAYS HOUR MINU SECO CLAT CLON CLATH CLONH HOLS' + hdr1b ='SAID FOVN YEAR MNTH DAYS HOUR MINU SECO CLAT CLON CLATH CLONH HOLS SACV' hdr2b ='SAZA SOZA BEARAZ SOLAZI' allocate(data_all(nele,itxmax),data1b8(nchanl),data1b4(nchanl),nrec(itxmax)) @@ -506,35 +516,72 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,& call openbf(lnbufr,'IN',lnbufr) - if(llll >= 2 .and. (amsua .or. amsub .or. mhs))then + ! support multiple spc coefficient files for any given sensor + if(amsua .or. amsub .or. mhs)then quiet=.not.verbose allocate(data1b8x(nchanl)) - sensorlist(1)=sis - if( crtm_coeffs_path /= "" ) then - if(mype_sub==mype_root .and. print_verbose) write(6,*)'READ_BUFRTOVS: crtm_spccoeff_load() on path "'//trim(crtm_coeffs_path)//'"' - error_status = crtm_spccoeff_load(sensorlist,& - File_Path = crtm_coeffs_path, quiet=quiet ) + spc_coeff_versions = 0 + spc_coeff_found = .true. + do while (spc_coeff_found) + if (spc_coeff_versions == 0) then + sensorlist(1)=sis else - error_status = crtm_spccoeff_load(sensorlist,quiet=quiet) - endif - if (error_status /= success) then - write(6,*)'READ_BUFRTOVS: ***ERROR*** crtm_spccoeff_load error_status=',error_status,& - ' TERMINATE PROGRAM EXECUTION' - call stop2(71) - endif - ninstruments = size(sc) - instrument_loop: do n=1,ninstruments - if(sis == sc(n)%sensor_id)then - instrument=n - exit instrument_loop + i = spc_coeff_versions+1 + write(sensorlist(1),'(a,a,i1)') trim(sis),'_v',i end if - end do instrument_loop - if(instrument == 0)then - write(6,*)' failure to find instrument in read_bufrtovs ',sis - end if + + if( crtm_coeffs_path /= "" ) then + if(mype_sub==mype_root .and. print_verbose) write(6,*)'READ_BUFRTOVS: crtm_spccoeff_load() on path "'//trim(crtm_coeffs_path)//'"' + end if + + spc_filename = trim(crtm_coeffs_path) // trim(sensorlist(1)) // '.SpcCoeff.bin' + INQUIRE(FILE=trim(spc_filename), EXIST=spc_coeff_found) + + if (.NOT. spc_coeff_found) then + if (spc_coeff_versions == 0) then + write(6,*)'READ_BUFRTOVS: ***ERROR*** crtm_spccoeff_load error_status=',error_status,& + ' TERMINATE PROGRAM EXECUTION' + call stop2(71) + else + write(6,*)'READ_BUFRTOVS: ', spc_coeff_versions, ' versions of SpcCoeff found for ', trim(sis) + end if + else + spc_coeff_versions = spc_coeff_versions+1 + + if( crtm_coeffs_path /= "" ) then + error_status = crtm_spccoeff_load(sensorlist,& + File_Path = crtm_coeffs_path, quiet=quiet ) + else + error_status = crtm_spccoeff_load(sensorlist,quiet=quiet) + endif + if (error_status /= success) then + write(6,*)'READ_BUFRTOVS: ***ERROR*** crtm_spccoeff_load error_status=',error_status,& + ' despite file ',trim(spc_filename),' existing, TERMINATE PROGRAM EXECUTION' + call stop2(71) + endif + + ninstruments = size(sc) + instrument_loop: do n=1,ninstruments + if(sis == sc(n)%sensor_id)then + instrument=n + exit instrument_loop + end if + end do instrument_loop + if(instrument == 0)then + write(6,*)' failure to find instrument in read_bufrtovs ',sis + end if + + accoeff_sets(spc_coeff_versions) = sc(instrument)%ac + + ! deallocate crtm info + error_status = crtm_spccoeff_destroy() + if (error_status /= success) & + write(6,*)'OBSERVER: ***ERROR*** crtm_spccoeff_destroy error_status=',error_status + end if + end do + end if - ! Loop to read bufr file irecx=0 read_subset: do while(ireadmg(lnbufr,subset,idate)>=0) @@ -633,6 +680,17 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,& call map2tgrid(dlat_earth,dlon_earth,dist1,crit1,itx,ithin,itt,iuse,sis,it_mesh=it_mesh) if(.not. iuse)cycle read_loop +! Extract satellite antenna corrections version number + if (llll > 1) then + sacv = nint(bfr1bhdr(14)) + if (sacv > spc_coeff_versions) then + write(6,*) 'READ_BUFRTOVS WARNING sacv greater than spc_coeff_versions' + end if + else ! normal feed doesn't have antenna correction, so set sacv to 0 + sacv = 0 + end if + + call ufbint(lnbufr,bfr2bhdr,n2bhdr,1,iret,hdr2b) ! Set common predictor parameters @@ -674,22 +732,66 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,& endif ! Read data record. Increment data counter -! TMBR is actually the antenna temperature for most microwave -! sounders. - if (llll == 1) then - call ufbrep(lnbufr,data1b8,1,nchanl,iret,'TMBR') - else ! EARS / DB - call ufbrep(lnbufr,data1b8,1,nchanl,iret,'TMBRST') - if ( amsua .or. amsub .or. mhs )then - data1b8x=data1b8 - data1b4=data1b8 - call remove_antcorr(sc(instrument)%ac,ifov,data1b4) - data1b8=data1b4 - do j=1,nchanl - if(data1b8x(j) > r1000)data1b8(j) = 1000000._r_kind - end do +! TMBR is actually the antenna temperature for most microwave sounders. +! Changed from accepting TMBR and converting TMBRST to +! antenna temperature (remove_antcorr), to converting TMBR to +! brightness temperature (add_antcorr) and accepting TMBRST +! + + if (ta2tb) then + + if (llll == 1) then + call ufbrep(lnbufr,data1b8,1,nchanl,iret,'TMBR') + if ( amsua .or. amsub .or. mhs )then + ! convert antenna temperature to brightness temperature + data1b8x=data1b8 + data1b4=data1b8 + !call apply_antcorr(accoeff_sets(spc_coeff_versions),ifov,data1b4) + call apply_antcorr(accoeff_sets(1),ifov,data1b4) + data1b8=data1b4 + do j=1,nchanl + if(data1b8x(j) > r1000) data1b8(j) = 1000000._r_kind + end do + end if + else ! EARS / DB + call ufbrep(lnbufr,data1b8,1,nchanl,iret,'TMBRST') + !if ( amsua .or. amsub .or. mhs .AND. sacv .ne. spc_coeff_versions)then + if ( amsua .or. amsub .or. mhs .AND. spc_coeff_versions /= 1 .AND. sacv /= 1)then + ! convert brightness temperature to antenna temperature using + ! the satellite antenna correction version (sacv) used by the + ! data originator, + ! then convert back to brightness temperature using the version + ! of parameters used by the CRTM + data1b8x=data1b8 + data1b4=data1b8 + call remove_antcorr(accoeff_sets(sacv),ifov,data1b4) + !call apply_antcorr(accoeff_sets(spc_coeff_versions),ifov,data1b4) + call apply_antcorr(accoeff_sets(1),ifov,data1b4) + data1b8=data1b4 + do j=1,nchanl + if(data1b8x(j) > r1000) data1b8(j) = 1000000._r_kind + end do + end if end if - end if + + else + + if (llll == 1) then + call ufbrep(lnbufr,data1b8,1,nchanl,iret,'TMBR') + else ! EARS / DB + call ufbrep(lnbufr,data1b8,1,nchanl,iret,'TMBRST') + if ( amsua .or. amsub .or. mhs )then + data1b8x=data1b8 + data1b4=data1b8 + call remove_antcorr(accoeff_sets(1),ifov,data1b4) + data1b8=data1b4 + do j=1,nchanl + if(data1b8x(j) > r1000)data1b8(j) = 1000000._r_kind + end do + end if + end if + + endif ! Transfer observed brightness temperature to work array. If any ! temperature exceeds limits, reset observation to "bad" value @@ -947,14 +1049,7 @@ subroutine read_bufrtovs(mype,val_tovs,ithin,isfcalc,& call closbf(lnbufr) close(lnbufr) - if(llll > 1 .and. (amsua .or. amsub .or. mhs))then - deallocate(data1b8x) - -! deallocate crtm info - error_status = crtm_spccoeff_destroy() - if (error_status /= success) & - write(6,*)'OBSERVER: ***ERROR*** crtm_spccoeff_destroy error_status=',error_status - end if + if (allocated(data1b8x)) deallocate(data1b8x) end do ears_db_loop deallocate(data1b8,data1b4) diff --git a/src/gsi/read_gmi.f90 b/src/gsi/read_gmi.f90 index 4361e7144a..0f45aa7e28 100644 --- a/src/gsi/read_gmi.f90 +++ b/src/gsi/read_gmi.f90 @@ -131,7 +131,10 @@ subroutine read_gmi(mype,val_gmi,ithin,rmesh,jsatid,gstime,& integer(i_kind),parameter :: maxchanl=13, ngs=2 integer(i_kind),dimension(maxchanl) :: tbmin ! different tbmin for the channels. - real(r_double),dimension(maxchanl) :: mirad,gmichq,gmirfi ! TBB from strtmbr + !real(r_double),dimension(maxchanl) :: mirad,gmichq,gmirfi ! TBB from strtmbr + real(r_double),dimension(maxchanl) :: mirad,var_check1,var_check2 ! TBB from strtmbr + logical :: use_gmichq + real(r_double) :: fovn,slnm ! FOVN real(r_kind),parameter :: r360=360.0_r_kind character(80),parameter :: satinfo='SAID SIID OGCE GSES SACV' !use for ufbint() @@ -408,12 +411,21 @@ subroutine read_gmi(mype,val_gmi,ithin,rmesh,jsatid,gstime,& ! ----- Read header record to extract obs location information call ufbint(lnbufr,midat,nloc,1,iret,'SCLAT SCLON HMSL') - call ufbrep(lnbufr,gmichq,1,nchanl,iret,'GMICHQ') - call ufbrep(lnbufr,gmirfi,1,nchanl,iret,'GMIRFI') + call ufbrep(lnbufr,var_check1,1,nchanl,iret,'GMICHQ') + !call ufbrep(lnbufr,gmirfi,1,nchanl,iret,'GMIRFI') call ufbrep(lnbufr,pixelsaza,1,ngs,iret,'SAZA') call ufbrep(lnbufr,val_angls,n_angls,ngs,iret,'SAMA SZA SMA SGA') call ufbint(lnbufr,pixelloc,2, 1,iret,'CLATH CLONH') + if (any(var_check1 < 99999999999_r_double)) then ! 100000000000 seems to be the missing value + use_gmichq = .true. + call ufbrep(lnbufr,var_check2,1,nchanl,iret,'GMIRFI') + else + use_gmichq = .false. + call ufbrep(lnbufr,var_check1,1,nchanl,iret,'VIIRSQ') + call ufbrep(lnbufr,var_check2,1,nchanl,iret,'TPQC2') + endif + !--- Extract brightness temperature data. Apply gross check to data. ! If obs fails gross check, reset to missing obs value. call ufbrep(lnbufr,mirad,1,nchanl,iret,strtmbr) @@ -479,13 +491,32 @@ subroutine read_gmi(mype,val_gmi,ithin,rmesh,jsatid,gstime,& iskip=0 do jc=1, nchanla ! only does such check the first 9 channels for GMI 1C-R data - if( mirad(jc)tbmax .or. & - gmichq(jc) < -0.5_r_kind .or. gmichq(jc) > 1.5_r_kind .or. & - gmirfi(jc)>0.0_r_kind) then ! & - iskip = iskip + 1 + if( mirad(jc)tbmax ) then + iskip = iskip+1 + !endif + !if(use_gmichq) then + elseif (use_gmichq) then + if (var_check1(jc) < -0.5_r_kind .or. var_check1(jc) > 1.5_r_kind .or. & + var_check2(jc) > 0.0_r_kind) then + iskip = iskip+1 + else + nread=nread+1 + endif else - nread=nread+1 - end if + if (var_check1(jc) < -0.5_r_kind .or. var_check1(jc) > 0.5_r_kind .or. & + var_check2(jc) > 0.0_r_kind) then + iskip = iskip+1 + else + nread=nread+1 + endif + endif + !if( mirad(jc)tbmax .or. & + ! gmichq(jc) < -0.5_r_kind .or. gmichq(jc) > 1.5_r_kind .or. & + ! gmirfi(jc)>0.0_r_kind) then ! & + ! iskip = iskip + 1 + !else + ! nread=nread+1 + !end if enddo if(iskip == nchanla) then @@ -564,6 +595,7 @@ subroutine read_gmi(mype,val_gmi,ithin,rmesh,jsatid,gstime,& allocate(data_all(nele,itxmax),nrec(itxmax)) nrec=999999 + if (mype==0) write(*,*) 'read_gmi num_obs before thinning: ', num_obs obsloop: do iobs = 1, num_obs t4dv => t4dv_save(iobs) diff --git a/src/gsi/setuprad.f90 b/src/gsi/setuprad.f90 index 66774c605a..491afb402d 100644 --- a/src/gsi/setuprad.f90 +++ b/src/gsi/setuprad.f90 @@ -280,7 +280,7 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& use sst_retrieval, only: setup_sst_retrieval,avhrr_sst_retrieval,& finish_sst_retrieval,spline_cub use m_dtime, only: dtime_setup, dtime_check - use crtm_interface, only: init_crtm,call_crtm,destroy_crtm,sensorindex,surface, & + use crtm_interface, only: init_crtm,call_crtm,destroy_crtm,sensorindex,surface,& itime,ilon,ilat,ilzen_ang,ilazi_ang,iscan_ang,iscan_pos,iszen_ang,isazi_ang, & ifrac_sea,ifrac_lnd,ifrac_ice,ifrac_sno,itsavg, & izz,idomsfc,isfcr,iff10,ilone,ilate, & @@ -357,8 +357,10 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& real(r_kind) bias real(r_kind) factch6 real(r_kind) stability,tcwv,hwp_ratio - real(r_kind) si_obs,si_fg,si_mean - + real(r_kind) si_obs,si_fg +! real(r_kind) si_mean + real(r_kind) total_cloud_cover + logical cao_flag logical hirs2,msu,goessndr,hirs3,hirs4,hirs,amsua,amsub,airs,hsb,goes_img,ahi,mhs,abi type(sparr2) :: dhx_dx @@ -889,6 +891,7 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& ! Output both tsim and tsim_clr for allsky tsim_clr=zero tcc=zero + total_cloud_cover=zero if (radmod%lcloud_fwd) then call call_crtm(obstype,dtime,data_s(:,n),nchanl,nreal,ich, & tvp,qvp,clw_guess,ciw_guess,rain_guess,snow_guess,prsltmp,prsitmp, & @@ -905,7 +908,8 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& tvp,qvp,clw_guess,ciw_guess,rain_guess,snow_guess,prsltmp,prsitmp, & trop5,tzbgr,dtsavg,sfc_speed, & tsim2,emissivity2,ptau52,ts2,emissivity_k2, & - temp2,wmix2,jacobian2,error_status,tsim_clr2) + temp2,wmix2,jacobian2,error_status,tsim_clr=tsim_clr2,tcc=tcc,& + tcwv=tcwv,hwp_ratio=hwp_ratio,stability=stability) ! merge emissivity(10:13) = emissivity2(10:13) ts(10:13) = ts2(10:13) @@ -921,6 +925,8 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& tsim_clr(10:13) = tsim_clr2(10:13) cosza2 = cos(data_s(ilzen_ang2,n)) endif + total_cloud_cover = tcc(1) + cld = total_cloud_cover else call call_crtm(obstype,dtime,data_s(:,n),nchanl,nreal,ich, & tvp,qvp,clw_guess,ciw_guess,rain_guess,snow_guess,prsltmp,prsitmp, & @@ -1047,7 +1053,7 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& endif endif ! Screening for cold-air outbreak area (only applied to MW for now) - if (cao_check) then + if (cao_check .and. radmod%lprecip) then if(microwave .and. sea) then if(radmod%lcloud_fwd) then cao_flag = (stability < 12.0_r_kind) .and. (hwp_ratio < half) .and. (tcwv < 8.0_r_kind) @@ -1368,7 +1374,7 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& end if call qc_amsua(nchanl,is,ndat,nsig,npred,sea,land,ice,snow,mixed,luse(n), & - zsges,cenlat,tb_obsbc1,si_mean,cosza,clw_obs,tbc,ptau5,emissivity_k,ts, & + zsges,cenlat,tb_obsbc1,cosza,clw_obs,tbc,ptau5,emissivity_k,ts, & pred,predchan,id_qc,aivals,errf,errf0,clw_obs,varinv,cldeff_obs,cldeff_fg,factch6, & cld_rbc_idx,sfc_speed,error0,clw_guess_retrieval,scatp,radmod) @@ -1406,10 +1412,10 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& end if si_obs = (tb_obsbc16-tb_obsbc17) - (tsim_clr(16)-tsim_clr(17)) si_fg = (tsim(16)-tsim(17)) - (tsim_clr(16)-tsim_clr(17)) - si_mean= half*(si_obs+si_fg) +! si_mean= half*(si_obs+si_fg) call qc_atms(nchanl,is,ndat,nsig,npred,sea,land,ice,snow,mixed,luse(n), & - zsges,cenlat,tb_obsbc1,si_mean,cosza,clw_obs,tbc,ptau5,emissivity_k,ts, & + zsges,cenlat,tb_obsbc1,cosza,clw_obs,tbc,ptau5,emissivity_k,ts, & pred,predchan,id_qc,aivals,errf,errf0,clw_obs,varinv,cldeff_obs,cldeff_fg,factch6, & cld_rbc_idx,sfc_speed,error0,clw_guess_retrieval,scatp,radmod) @@ -1606,9 +1612,17 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& m=ich(i) if(radmod%lcloud_fwd .and. eff_area) then if(radmod%rtype == 'amsua' .and. (i <=5 .or. i==15) ) then - errf(i) = three*errf(i) + if (radmod%lprecip) then + errf(i) = 2.5_r_kind*errf(i) + else + errf(i) = three*errf(i) + endif else if(radmod%rtype == 'atms' .and. (i <= 6 .or. i>=16) ) then - errf(i) = min(three*errf(i),10.0_r_kind) + if (radmod%lprecip) then + errf(i) = min(2.5_r_kind*errf(i),10.0_r_kind) + else + errf(i) = min(three*errf(i),10.0_r_kind) + endif else if(radmod%rtype == 'gmi') then errf(i) = min(2.0_r_kind*errf(i),ermax_rad(m)) else if (radmod%rtype/='amsua' .and. radmod%rtype/='atms' .and. radmod%rtype/='gmi' .and. radmod%lcloud4crtm(i)>=0) then diff --git a/src/gsi/write_incr.f90 b/src/gsi/write_incr.f90 index 6563775662..69ad96e281 100644 --- a/src/gsi/write_incr.f90 +++ b/src/gsi/write_incr.f90 @@ -91,6 +91,9 @@ subroutine write_fv3_inc_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) use obsmod, only: ianldate use state_vectors, only: allocate_state, deallocate_state + use state_vectors, only: svars3d + use mpeu_util, only: getindex + implicit none ! !INPUT PARAMETERS: @@ -107,6 +110,7 @@ subroutine write_fv3_inc_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) real(r_kind),pointer,dimension(:,:,:) :: sub_u,sub_v real(r_kind),pointer,dimension(:,:,:) :: sub_qanl,sub_oz real(r_kind),pointer,dimension(:,:,:) :: sub_ql, sub_qi + real(r_kind),pointer,dimension(:,:,:) :: sub_qr, sub_qs, sub_qg real(r_kind),pointer,dimension(:,:) :: sub_ps real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig) :: sub_dzb,sub_dza, sub_tsen, sub_q @@ -114,7 +118,7 @@ subroutine write_fv3_inc_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) real(r_kind),dimension(grd%lat1,grd%lon1) :: pssm real(r_kind),dimension(grd%lat1,grd%lon1,grd%nsig):: tsensm, usm, vsm real(r_kind),dimension(grd%lat1,grd%lon1,grd%nsig):: qsm, ozsm - real(r_kind),dimension(grd%lat1,grd%lon1,grd%nsig):: qism, qlsm + real(r_kind),dimension(grd%lat1,grd%lon1,grd%nsig):: qism, qlsm, qrsm, qssm, qgsm real(r_kind),dimension(grd%lat1,grd%lon1,grd%nsig):: dzsm real(r_kind),dimension(grd%lat1,grd%lon1,grd%nsig):: delp real(r_kind),dimension(grd%nlon) :: deglons @@ -128,7 +132,9 @@ subroutine write_fv3_inc_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) 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, & + qrvarid, qsvarid, qgvarid + integer(i_kind) :: iql,iqi,iqr,iqs,iqg integer(i_kind) :: dimids3(3),nccount(3),ncstart(3), cnksize(3), j1, j2 type(gsi_bundle) :: svalinc(nobs_bins) @@ -188,10 +194,20 @@ subroutine write_fv3_inc_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) end if end if + ! Check hydrometeors in control variables + iql = getindex(svars3d,'ql') + iqi = getindex(svars3d,'qi') + iqr = getindex(svars3d,'qr') + iqs = getindex(svars3d,'qs') + iqg = getindex(svars3d,'qg') + istatus=0 call gsi_bundlegetpointer(gfs_bundle,'q', sub_qanl, iret); istatus=istatus+iret - call gsi_bundlegetpointer(svalinc(ibin2),'ql', sub_ql, iret); istatus=istatus+iret - call gsi_bundlegetpointer(svalinc(ibin2),'qi', sub_qi, iret); istatus=istatus+iret + if (iql>0) call gsi_bundlegetpointer(svalinc(ibin2),'ql', sub_ql, iret); istatus=istatus+iret + if (iqi>0) call gsi_bundlegetpointer(svalinc(ibin2),'qi', sub_qi, iret); istatus=istatus+iret + if (iqr>0) call gsi_bundlegetpointer(svalinc(ibin2),'qr', sub_qr, iret); istatus=istatus+iret + if (iqs>0) call gsi_bundlegetpointer(svalinc(ibin2),'qs', sub_qs, iret); istatus=istatus+iret + if (iqg>0) call gsi_bundlegetpointer(svalinc(ibin2),'qg', sub_qg, iret); istatus=istatus+iret call gsi_bundlegetpointer(svalinc(ibin2),'oz', sub_oz, iret); istatus=istatus+iret call gsi_bundlegetpointer(svalinc(ibin2),'u', sub_u, iret); istatus=istatus+iret call gsi_bundlegetpointer(svalinc(ibin2),'v', sub_v, iret); istatus=istatus+iret @@ -235,12 +251,28 @@ subroutine write_fv3_inc_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) call nccheck_incr(nf90_var_par_access(ncid_out, tvarid, nf90_collective)) call nccheck_incr(nf90_def_var(ncid_out, "sphum_inc", nf90_real, dimids3, sphumvarid)) call nccheck_incr(nf90_var_par_access(ncid_out, sphumvarid, nf90_collective)) - call nccheck_incr(nf90_def_var(ncid_out, "liq_wat_inc", nf90_real, dimids3, liqwatvarid)) - call nccheck_incr(nf90_var_par_access(ncid_out, liqwatvarid, nf90_collective)) + if (iql>0) then + call nccheck_incr(nf90_def_var(ncid_out, "liq_wat_inc", nf90_real, dimids3, liqwatvarid)) + call nccheck_incr(nf90_var_par_access(ncid_out, liqwatvarid, nf90_collective)) + endif call nccheck_incr(nf90_def_var(ncid_out, "o3mr_inc", nf90_real, dimids3, o3varid)) call nccheck_incr(nf90_var_par_access(ncid_out, o3varid, nf90_collective)) - call nccheck_incr(nf90_def_var(ncid_out, "icmr_inc", nf90_real, dimids3, icvarid)) - call nccheck_incr(nf90_var_par_access(ncid_out, icvarid, nf90_collective)) + if (iqi>0) then + call nccheck_incr(nf90_def_var(ncid_out, "icmr_inc", nf90_real, dimids3, icvarid)) + call nccheck_incr(nf90_var_par_access(ncid_out, icvarid, nf90_collective)) + endif + if (iqr>0) then + call nccheck_incr(nf90_def_var(ncid_out, "rwmr_inc", nf90_real, dimids3, qrvarid)) + call nccheck_incr(nf90_var_par_access(ncid_out, qrvarid, nf90_collective)) + endif + if (iqs>0) then + call nccheck_incr(nf90_def_var(ncid_out, "snmr_inc", nf90_real, dimids3, qsvarid)) + call nccheck_incr(nf90_var_par_access(ncid_out, qsvarid, nf90_collective)) + endif + if (iqg>0) then + call nccheck_incr(nf90_def_var(ncid_out, "grle_inc", nf90_real, dimids3, qgvarid)) + call nccheck_incr(nf90_var_par_access(ncid_out, qgvarid, nf90_collective)) + endif ! place global attributes to parallel calc_increment output call nccheck_incr(nf90_put_att(ncid_out, nf90_global, "source", "GSI")) call nccheck_incr(nf90_put_att(ncid_out, nf90_global, "comment", & @@ -274,8 +306,11 @@ subroutine write_fv3_inc_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) ! Strip off boundary points from subdomains call strip(sub_tsen ,tsensm ,grd%nsig) call strip(sub_q ,qsm ,grd%nsig) - call strip(sub_ql ,qlsm ,grd%nsig) - call strip(sub_qi ,qism ,grd%nsig) + if (iql>0) call strip(sub_ql ,qlsm ,grd%nsig) + if (iqi>0) call strip(sub_qi ,qism ,grd%nsig) + if (iqr>0) call strip(sub_qr ,qrsm ,grd%nsig) + if (iqs>0) call strip(sub_qs ,qssm ,grd%nsig) + if (iqg>0) call strip(sub_qg ,qgsm ,grd%nsig) call strip(sub_oz ,ozsm ,grd%nsig) call strip(sub_ps ,pssm ) call strip(sub_u ,usm ,grd%nsig) @@ -410,17 +445,19 @@ subroutine write_fv3_inc_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) start = ncstart, count = nccount)) call mpi_barrier(mpi_comm_world,ierror) ! liquid water increment - do k=1,grd%nsig - krev = grd%nsig+1-k - if (zero_increment_strat('liq_wat_inc')) then - call zero_inc_strat(qlsm(:,:,k), k, troplev) - end if - if (should_zero_increments_for('liq_wat_inc')) qlsm(:,:,k) = 0.0_r_kind - out3d(:,:,krev) = transpose(qlsm(j1:j2,:,k)) - end do - call nccheck_incr(nf90_put_var(ncid_out, liqwatvarid, sngl(out3d), & - start = ncstart, count = nccount)) - call mpi_barrier(mpi_comm_world,ierror) + if (iql>0) then + do k=1,grd%nsig + krev = grd%nsig+1-k + if (zero_increment_strat('liq_wat_inc')) then + call zero_inc_strat(qlsm(:,:,k), k, troplev) + end if + if (should_zero_increments_for('liq_wat_inc')) qlsm(:,:,k) = 0.0_r_kind + out3d(:,:,krev) = transpose(qlsm(j1:j2,:,k)) + end do + call nccheck_incr(nf90_put_var(ncid_out, liqwatvarid, sngl(out3d), & + start = ncstart, count = nccount)) + call mpi_barrier(mpi_comm_world,ierror) + endif ! ozone increment do k=1,grd%nsig krev = grd%nsig+1-k @@ -434,17 +471,61 @@ subroutine write_fv3_inc_ (grd,sp_a,filename,mype_out,gfs_bundle,ibin) start = ncstart, count = nccount)) call mpi_barrier(mpi_comm_world,ierror) ! ice mixing ratio increment - do k=1,grd%nsig - krev = grd%nsig+1-k - if (zero_increment_strat('icmr_inc')) then - call zero_inc_strat(qism(:,:,k), k, troplev) - end if - if (should_zero_increments_for('icmr_inc')) qism(:,:,k) = 0.0_r_kind - out3d(:,:,krev) = transpose(qism(j1:j2,:,k)) - end do - call nccheck_incr(nf90_put_var(ncid_out, icvarid, sngl(out3d), & - start = ncstart, count = nccount)) - call mpi_barrier(mpi_comm_world,ierror) + if (iqi>0) then + do k=1,grd%nsig + krev = grd%nsig+1-k + if (zero_increment_strat('icmr_inc')) then + call zero_inc_strat(qism(:,:,k), k, troplev) + end if + if (should_zero_increments_for('icmr_inc')) qism(:,:,k) = 0.0_r_kind + out3d(:,:,krev) = transpose(qism(j1:j2,:,k)) + end do + call nccheck_incr(nf90_put_var(ncid_out, icvarid, sngl(out3d), & + start = ncstart, count = nccount)) + call mpi_barrier(mpi_comm_world,ierror) + endif + ! rain water mixing ratio increment + if (iqr>0) then + do k=1,grd%nsig + krev = grd%nsig+1-k + if (zero_increment_strat('rwmr_inc')) then + call zero_inc_strat(qrsm(:,:,k), k, troplev) + end if + if (should_zero_increments_for('rwmr_inc')) qrsm(:,:,k) = 0.0_r_kind + out3d(:,:,krev) = transpose(qrsm(j1:j2,:,k)) + end do + call nccheck_incr(nf90_put_var(ncid_out, qrvarid, sngl(out3d), & + start = ncstart, count = nccount)) + call mpi_barrier(mpi_comm_world,ierror) + endif + ! snow water mixing ratio increment + if (iqs>0) then + do k=1,grd%nsig + krev = grd%nsig+1-k + if (zero_increment_strat('snmr_inc')) then + call zero_inc_strat(qssm(:,:,k), k, troplev) + end if + if (should_zero_increments_for('snmr_inc')) qssm(:,:,k) = 0.0_r_kind + out3d(:,:,krev) = transpose(qssm(j1:j2,:,k)) + end do + call nccheck_incr(nf90_put_var(ncid_out, qsvarid, sngl(out3d), & + start = ncstart, count = nccount)) + call mpi_barrier(mpi_comm_world,ierror) + endif + ! graupel mixing ratio increment + if (iqg>0) then + do k=1,grd%nsig + krev = grd%nsig+1-k + if (zero_increment_strat('grle_inc')) then + call zero_inc_strat(qgsm(:,:,k), k, troplev) + end if + if (should_zero_increments_for('grle_inc')) qgsm(:,:,k) = 0.0_r_kind + out3d(:,:,krev) = transpose(qgsm(j1:j2,:,k)) + end do + call nccheck_incr(nf90_put_var(ncid_out, qgvarid, sngl(out3d), & + start = ncstart, count = nccount)) + call mpi_barrier(mpi_comm_world,ierror) + endif ! ! cleanup and exit call nccheck_incr(nf90_close(ncid_out)) if ( mype == mype_out ) then From f16e852643570f7d030dd467e83ec5ba7849260e Mon Sep 17 00:00:00 2001 From: "edward.safford" Date: Fri, 29 Apr 2022 14:39:22 +0000 Subject: [PATCH 2/3] Github issue #310. Correct mean correction time series plots. Fixes #310. --- .../image_gen/html/Install_html.sh | 37 ++-- .../image_gen/html/comp.html.glb | 90 --------- .../image_gen/html/comp.html.rgn | 90 --------- .../image_gen/html/install_glb.sh | 189 ++++++++---------- .../image_gen/html/install_rgn.sh | 88 ++------ .../image_gen/html/menu.html.glb | 53 ----- .../image_gen/html/menu.html.rgn | 53 ----- .../image_gen/html/mk_intro.sh | 51 +---- .../image_gen/html/sum_thumbs.tar | Bin 0 -> 204800 bytes .../src/radmon_ig_angle.fd/angle.f90 | 36 +--- util/Radiance_Monitor/parm/RadMon_config | 10 +- .../parm/RadMon_user_settings | 2 +- 12 files changed, 142 insertions(+), 557 deletions(-) delete mode 100755 util/Radiance_Monitor/image_gen/html/comp.html.glb delete mode 100644 util/Radiance_Monitor/image_gen/html/comp.html.rgn delete mode 100755 util/Radiance_Monitor/image_gen/html/menu.html.glb delete mode 100644 util/Radiance_Monitor/image_gen/html/menu.html.rgn create mode 100644 util/Radiance_Monitor/image_gen/html/sum_thumbs.tar diff --git a/util/Radiance_Monitor/image_gen/html/Install_html.sh b/util/Radiance_Monitor/image_gen/html/Install_html.sh index 8ecb71d1da..0d968be994 100755 --- a/util/Radiance_Monitor/image_gen/html/Install_html.sh +++ b/util/Radiance_Monitor/image_gen/html/Install_html.sh @@ -1,5 +1,4 @@ #!/bin/sh -set -ax #-------------------------------------------------------------------- #-------------------------------------------------------------------- @@ -27,7 +26,7 @@ if [[ $nargs -ne 2 ]]; then exit 2 fi -RADMON_SUFFIX=$1 +export RADMON_SUFFIX=$1 echo RADMON_SUFFIX = $RADMON_SUFFIX export RAD_AREA=$2 @@ -36,28 +35,38 @@ this_dir=`dirname $0` top_parm=${this_dir}/../../parm -if [[ -s ${top_parm}/RadMon_config ]]; then - . ${top_parm}/RadMon_config -else - echo "ERROR: Unable to source ${top_parm}/RadMon_config" - exit +radmon_config=${radmon_config:-${top_parm}/RadMon_config} +if [[ ! -e ${radmon_config} ]]; then + echo "Unable to source ${radmon_config}" + exit 2 fi -if [[ -s ${top_parm}/RadMon_user_settings ]]; then - . ${top_parm}/RadMon_user_settings -else - echo "ERROR: Unable to source ${top_parm}/RadMon_user_settings" - exit +. ${radmon_config} +if [[ $? -ne 0 ]]; then + echo "Error detected while sourcing ${radmon_config} file" + exit $? fi +radmon_user_settings=${radmon_user_settings:-${top_parm}/RadMon_user_settings} +if [[ ! -e ${radmon_user_settings} ]]; then + echo "Unable to locate ${radmon_user_settings} file" + exit 4 +fi + +. ${radmon_user_settings} +if [[ $? -ne 0 ]]; then + echo "Error detected while sourcing ${radmon_user_settings} file" + exit $? +fi + #-------------------------------------------------------------- # call the appropriate child script for glb or rgn # if [[ $RAD_AREA == "glb" ]]; then - ${RADMON_IMAGE_GEN}/html/install_glb.sh $RADMON_SUFFIX + ${RADMON_IMAGE_GEN}/html/install_glb.sh else - ${RADMON_IMAGE_GEN}/html/install_rgn.sh $RADMON_SUFFIX + ${RADMON_IMAGE_GEN}/html/install_rgn.sh fi diff --git a/util/Radiance_Monitor/image_gen/html/comp.html.glb b/util/Radiance_Monitor/image_gen/html/comp.html.glb deleted file mode 100755 index e1f9d07680..0000000000 --- a/util/Radiance_Monitor/image_gen/html/comp.html.glb +++ /dev/null @@ -1,90 +0,0 @@ - - - - - Radiance Monitoring - - - - - - - - -

Comparison plots

- -
- - - - - -
Select Platform:
-

-

Click button to...
-
-
- - -
- - - - - diff --git a/util/Radiance_Monitor/image_gen/html/comp.html.rgn b/util/Radiance_Monitor/image_gen/html/comp.html.rgn deleted file mode 100644 index e1f9d07680..0000000000 --- a/util/Radiance_Monitor/image_gen/html/comp.html.rgn +++ /dev/null @@ -1,90 +0,0 @@ - - - - - Radiance Monitoring - - - - - - - - -

Comparison plots

- -
- - - - - -
Select Platform:
-

-

Click button to...
-
-
- - -
- - - - - diff --git a/util/Radiance_Monitor/image_gen/html/install_glb.sh b/util/Radiance_Monitor/image_gen/html/install_glb.sh index 4c70951df9..5a54b5ff97 100755 --- a/util/Radiance_Monitor/image_gen/html/install_glb.sh +++ b/util/Radiance_Monitor/image_gen/html/install_glb.sh @@ -1,5 +1,4 @@ -#!/bin/sh -set -ax +#!/bin/bash #-------------------------------------------------------------------- #-------------------------------------------------------------------- @@ -11,40 +10,57 @@ set -ax #-------------------------------------------------------------------- #-------------------------------------------------------------------- -function usage { - echo "Usage: install_glb.sh suffix" - echo " Suffix is data source identifier that matches data in " - echo " the $TANKDIR/stats directory." -} - echo "BEGIN install_glb.sh" echo "" +echo "" + +do_cmp=0 +cmp_src="" + +#-------------------------------------------------------------- +# Allow user to enable comparison plots +# +echo "Do you wish to enable data plots to include comparison to" +echo " operational GDAS data, or another data source?" +echo "" +echo -n " Enter YES to enable comparison plots, any other input to disable. > " +read text +short=`echo $text | cut -c1` + +if [[ $short = "Y" || $short = "y" ]]; then + do_cmp=1 + cmp_src="GDAS" + + echo "Please specify the suffix of your comparison data source," + echo " or just hit the return key to use the operational GDAS as " + echo " the comparison source" + echo "" + echo -n " > " + read text + + if [[ ${#text} -gt 0 ]]; then + cmp_src=${text} + fi -nargs=$# -if [[ $nargs -lt 0 ]]; then - usage - exit 2 + echo "${cmp_src} will be used as the comparison source." fi -SUFFIX=$1 -echo SUFFIX = $SUFFIX + +SUFFIX=$RADMON_SUFFIX RAD_AREA="glb" this_file=`basename $0` this_dir=`dirname $0` - - new_webdir=${WEBDIR}/${SUFFIX} -echo RAD_AREA = $RAD_AREA -echo TANKverf = $TANKverf - #-------------------------------------------------------------- # Create a temporary working directory. # workdir=$STMP_USER/${SUFFIX}_html -rmdir $workdir +if [[ -e $workdir ]]; then + rm -rf $workdir +fi mkdir $workdir cd $workdir @@ -56,13 +72,10 @@ cd $workdir #----------------------------------------------------------- # Find the first date with data. Start at today and work -# backwards. Stop after 90 days and exit. +# backwards. If not found stop after 90 days and exit. # PDATE=`${IG_SCRIPTS}/nu_find_cycle.pl --dir ${TANKverf} --cyc 1` -echo PDATE= $PDATE - limit=`$NDATE -2160 $PDATE` # 90 days -echo limit, PDATE = $limit, $PDATE #----------------------------------------------------------- # Build test_list which will contain all data files for @@ -81,14 +94,22 @@ while [[ data_found -eq 0 && $PDATE -ge $limit ]]; do test_dir=${TANKverf}/${RUN}.${PDY} fi - echo "test_dir = ${test_dir}" - if [[ -d ${test_dir} ]]; then echo " test_dir is GO " - if [[ -e ${test_dir}/angle.tar ]]; then - test_list=`tar -tv ${test_dir}/angle.tar` + if [[ -e ${test_dir}/radmon_angle.tar || -e ${test_dir}/radmon_angle.tar.gz ]]; then + gzipped=0 + if [[ -e ${test_dir}/radmon_angle.tar.gz ]]; then + gunzip ${test_dir}/radmon_angle.tar.gz + gzipped=1 + fi + + test_list=`tar -tf ${test_dir}/radmon_angle.tar` data_found=1 + + if [[ $gzipped -eq 1 ]]; then + gzip ${test_dir}/radmon_angle.tar + fi else test=`ls ${test_dir}/angle.*${PDATE}*.ieee_d* | wc -l` if [[ $test -gt 0 ]]; then @@ -102,10 +123,10 @@ while [[ data_found -eq 0 && $PDATE -ge $limit ]]; do if [[ data_found -eq 0 ]]; then PDATE=`$NDATE -24 $PDATE` - echo PDATE = $PDATE fi done + if [[ $data_found -eq 0 ]]; then echo Unable to locate any data files in the past 90 days for $SUFFIX echo in $TANKverf/angle. @@ -115,25 +136,24 @@ fi #----------------------------------------------------------- # Go through test_list and identify all unique # sat_instrument combinations. The results are the -# SATYPE list for this source. +# SATYPE list for this source. Avoid _anl and ctl files +# since they will give us duplicate satypes. # - for test in ${test_list}; do this_file=`basename $test` + + test_anl=`echo $this_file | grep "_anl"` + test_ctl=`echo $this_file | grep "ctl"` + if [[ $test_anl != "" || $test_ctl != "" ]]; then + continue + fi + tmp=`echo "$this_file" | cut -d. -f1` if [[ $tmp == "angle" ]]; then tmp=`echo "$this_file" | cut -d. -f2` fi - - #---------------------------------------------------------- - # remove sat/instrument_anl names so we don't end up - # with both "airs_aqua" and "airs_aqua_anl" if analysis - # files are being generated for this source. - #---------------------------------------------------------- - test_anl=`echo $tmp | grep "_anl"` - if [[ $test_anl = "" ]]; then - SATYPE_LIST="$SATYPE_LIST $tmp" - fi + + SATYPE_LIST="$SATYPE_LIST $tmp" done export SATYPE=$SATYPE_LIST @@ -143,7 +163,6 @@ if [[ ${#SATYPE} -le 0 ]]; then exit fi -echo $SATYPE #-------------------------------------------------------------- # Use the SATYPE list to construct the platform table. @@ -171,6 +190,7 @@ for satype in $SATYPE; do sat="GOES-${sat_num}" elif [[ $char == "n" ]]; then sat="NOAA-${sat_num}" + else sat=`echo $sat | tr 'a-z' 'A-Z'` fi @@ -178,9 +198,8 @@ for satype in $SATYPE; do sat=`echo $sat | tr 'a-z' 'A-Z'` fi - #----------------------------------------------------------------- - # + # Certain instruments require specific formatting. # amsu_test=`expr match ${ins} "amsu"` hirs_test=`expr match ${ins} "hirs"` @@ -203,7 +222,7 @@ done #-------------------------------------------------------------- # Sort the list by Satellite # -`sort -d $UNSORTED_LIST > $SORTED_LIST` +`sort -d -u $UNSORTED_LIST > $SORTED_LIST` #-------------------------------------------------------------- # Read the sorted list and create the platform table @@ -234,7 +253,6 @@ done < "$SORTED_LIST" imgndir=`dirname ${IMGNDIR}` - #-------------------------------------------------------------- # Edit the html files to add the platform table to each. # @@ -261,39 +279,16 @@ for html_file in $mod_html_files; do # switch all 'INSERT_SUFFIX' tags to the actual suffix # and route output to $html_file and we're done. sed s/INSERT_SUFFIX/${SUFFIX}/g ${tmp_html} > ${html_file} + rm ${tmp_html} done #-------------------------------------------------------------- -# Optionally enable comparison plots to the operational -# GDAS data +# Enable comparison plots # -set +ax -comp_html_files="plot_summary.html plot_time.html" - -echo "Do you wish to enable data plots to include comparison to" -echo " operational GDAS data, or another data source?" -echo "" -echo -n " Enter YES to enable comparison plots, any other input to disable. > " -read text -short=`echo $text | cut -c1` +if [[ $do_cmp == 1 ]]; then -if [[ $short = "Y" || $short = "y" ]]; then - - cmp_src="GDAS" - - echo "Please specify the suffix of your comparison data source," - echo " or just hit the return key to use the operational GDAS as " - echo " the comparison source" - echo "" - echo -n " > " - read text - - if [[ ${#text} -gt 0 ]]; then - cmp_src=${text} - fi - - echo "Enabling ${cmp_src} as the comparison source." + comp_html_files="plot_summary.html plot_time.html" #------------------------------------------------------------------------- # If cmp_src == GDAS we only have to uncomment the comparison check box @@ -327,8 +322,6 @@ if [[ $short = "Y" || $short = "y" ]]; then done fi -set -ax - #-------------------------------------------------------------- # Generate the intro.html file. # @@ -336,21 +329,7 @@ $NCP ${RADMON_IMAGE_GEN}/html/mk_intro.sh . $NCP ${RADMON_IMAGE_GEN}/html/intro.html intro.html.stock ./mk_intro.sh - - -#-------------------------------------------------------------- -# Copy the menu.html file and change "Experimental" to -# "Operational" if the suffix is wopr or nrx (operational GDAS -# or NDAS. -# -$NCP ${RADMON_IMAGE_GEN}/html/menu.html.$RAD_AREA . - -if [[ $SUFFIX == "wopr" || $SUFFIX == "nrx" ]]; then - tmp_menu=./tmp_menu.html.${RAD_AREA} - sed s/Experimental/Operational/1 menu.html.${RAD_AREA} > ${tmp_menu} - mv -f ${tmp_menu} menu.html.${RAD_AREA} -fi - +rm mk_intro.sh #-------------------------------------------------------------- # Copy the index.html file and change INSERT_SUFFIX to actual suffix. @@ -360,7 +339,7 @@ new_index="index.html" $NCP ${RADMON_IMAGE_GEN}/html/${index_file} . sed s/INSERT_SUFFIX/${SUFFIX}/g $index_file > ${tmp_index} -if [[ $SUFFIX == "wopr" || $SUFFIX == "nrx" ]]; then +if [[ $SUFFIX == "GFS" || $SUFFIX == "nrx" ]]; then sed s/Experimental/Operational/1 ${tmp_index} > ${new_index} fi @@ -372,16 +351,13 @@ if [[ ! -s ${new_index} ]]; then fi fi -support_files="jsuri-1.1.1.js stats.js latest_cycle.php" +rm ./${index_file} #-------------------------------------------------------------- # Make starting directory in $imgndir and copy over html, # misc, and thumb images. # -subdirs="angle bcoef bcor comp horiz summary time" -subdirs="summary" - if [[ ! -d ${IMGNDIR} ]]; then mkdir -p ${IMGNDIR} fi @@ -390,7 +366,7 @@ imgndir=`dirname ${IMGNDIR}` #----------------------- # move html files to imgndir # -all_html_files="${mod_html_files} index.html menu.html intro.html" +all_html_files="${mod_html_files} index.html intro.html" for file in $all_html_files; do $NCP ${file} ${imgndir}/${file} done @@ -398,6 +374,7 @@ done #----------------------- # mk image dirs # +subdirs="angle bcoef summary time" for dir in $subdirs; do mkdir -p ${imgndir}/pngs/${dir} done @@ -405,6 +382,7 @@ done #----------------------- # js files # +support_files="jsuri-1.1.1.js stats.js latest_cycle.php" for file in $support_files; do $NCP ${RADMON_IMAGE_GEN}/html/${file} ${imgndir}/. done @@ -419,7 +397,7 @@ done #----------------------- # summary thumb images -# If any are missing dummy one in using a copy of sndrdr1_g15. +# If any are missing dummy one in using ssmis_f18. # thumbs="sum_thumbs.tar" $NCP ${RADMON_IMAGE_GEN}/html/${thumbs} ${imgndir}/pngs/summary/. @@ -429,7 +407,7 @@ rm -f ${thumbs} for satype in $SATYPE; do if [[ ! -e ${satype}.summary.png ]]; then - $NCP sndrd1_g15.summary.png ${satype}.summary.png + $NCP ssmis_f18.summary.png ${satype}.summary.png fi done @@ -447,19 +425,24 @@ done #--------------------------------------------------- # if on wcoss then cd $imgndir and do the rsync here # -if [[ $MY_MACHINE = "wcoss" ]]; then +if [[ $MY_MACHINE = "wcoss_d" || $MY_MACHINE = "wcoss2" ]]; then + if [[ ${imgndir} != "/" ]]; then # sanity check to avoid serious embarrassment /usr/bin/rsync -ave ssh --exclude *.ctl.${Z} ${imgndir}/ \ - ${WEB_USER}@${WEB_SVR}.ncep.noaa.gov:${WEBDIR}/ + ${WEB_USER}@${WEB_SVR}.ncep.noaa.gov:${WEBDIR}/${SUFFIX}/ fi + + ssh ${WEB_USER}@${WEB_SVR} mkdir ${WEBDIR}/${SUFFIX}/gdas + ssh ${WEB_USER}@${WEB_SVR} ln -s ${WEBDIR}/${SUFFIX}/pngs ${WEBDIR}/${SUFFIX}/gdas/pngs + fi #------------------------ # clean up $workdir # -#cd $workdir -#cd ../ -#rm -rf $workdir +cd $workdir +cd ../ +rm -rf $workdir echo "" echo "END install_glb.sh" diff --git a/util/Radiance_Monitor/image_gen/html/install_rgn.sh b/util/Radiance_Monitor/image_gen/html/install_rgn.sh index 6ceb1a4caf..4a1a1666e6 100755 --- a/util/Radiance_Monitor/image_gen/html/install_rgn.sh +++ b/util/Radiance_Monitor/image_gen/html/install_rgn.sh @@ -1,5 +1,4 @@ #!/bin/sh -set -ax #-------------------------------------------------------------------- #-------------------------------------------------------------------- @@ -10,60 +9,23 @@ set -ax # the server (rgn only, regional is handled by Install_html.sh). #-------------------------------------------------------------------- #-------------------------------------------------------------------- - -function usage { - echo "Usage: install_rgn.sh suffix" - echo " Suffix is data source identifier that matches data in " - echo " the $TANKDIR/stats directory." -} - echo "BEGIN install_rgn.sh" echo "" -nargs=$# -if [[ $nargs -lt 0 ]]; then - usage - exit 2 -fi - -SUFFIX=$1 -echo SUFFIX = $SUFFIX +SUFFIX=$RADMON_SUFFIX RAD_AREA="rgn" this_file=`basename $0` this_dir=`dirname $0` -top_parm=${this_dir}/../../parm - -if [[ -s ${top_parm}/RadMon_config ]]; then - . ${top_parm}/RadMon_config -else - echo "ERROR: Unable to source ${top_parm}/RadMon_config" - exit -fi - -if [[ -s ${top_parm}/RadMon_user_settings ]]; then - . ${top_parm}/RadMon_user_settings -else - echo "ERROR: Unable to source ${top_parm}/RadMon_user_settings" - exit -fi - - -#-------------------------------------------------------------- -# Get the area for this SUFFIX from the data_map file -# -new_webdir=${WEBDIR}/${SUFFIX} - -echo RAD_AREA = $RAD_AREA -echo TANKverf = $TANKverf - #-------------------------------------------------------------- # Create a temporary working directory. # workdir=$STMP_USER/${SUFFIX}_html -rmdir $workdir +if [[ -d ${workdir} ]]; then + rm -rf $workdir +fi mkdir $workdir cd $workdir @@ -75,13 +37,10 @@ cd $workdir #----------------------------------------------------------- # Find the first date with data. Start at today and work -# backwards. Stop after 90 days and exit. +# backwards. If not found, stop after 90 days and exit. # -PDATE=`${IG_SCRIPTS}/find_cycle.pl 1 ${TANKverf}` -echo PDATE= $PDATE - +PDATE=`${IG_SCRIPTS}/find_cycle.pl --cyc 1 --dir ${TANKverf}` limit=`$NDATE -2160 $PDATE` # 90 days -echo limit, PDATE = $limit, $PDATE #----------------------------------------------------------- # Build test_list which will contain all data files for @@ -106,8 +65,6 @@ while [[ data_found -eq 0 && $PDATE -ge $limit ]]; do test_dir=${TANKverf}/${MONITOR}.${PDY} fi - echo "test_dir = ${test_dir}" - if [[ -d ${test_dir} ]]; then echo " test_dir is GO " @@ -134,7 +91,6 @@ while [[ data_found -eq 0 && $PDATE -ge $limit ]]; do if [[ data_found -eq 0 ]]; then PDATE=`$NDATE -24 $PDATE` - echo PDATE = $PDATE fi done @@ -183,8 +139,6 @@ if [[ ${#SATYPE} -le 0 ]]; then exit fi -echo "SATYPE = $SATYPE" - #-------------------------------------------------------------- # Use the SATYPE list to construct the platform table. # @@ -245,9 +199,6 @@ done # `sort -d $UNSORTED_LIST > $SORTED_LIST` -echo -echo -cat $SORTED_LIST #-------------------------------------------------------------- # Read the sorted list and create the platform table @@ -285,7 +236,6 @@ imgndir=`dirname ${IMGNDIR}` mod_html_files="plot_summary.html.rgn plot_time.html.rgn plot_angle.html.rgn plot_bcoef.html.rgn" for html_file in $mod_html_files; do - echo "processing ${html_file}" $NCP ${RADMON_IMAGE_GEN}/html/${html_file} . tmp_html=./tmp_${html_file} @@ -304,7 +254,6 @@ for html_file in $mod_html_files; do # rm the .rgn extension on the html file names new_html_file=`echo $html_file | cut -d'.' -f-2` - echo "new_html_file = $new_html_file" # switch all 'INSERT_SUFFIX' tags to the actual suffix # and route output to $new_html_file and we're done. @@ -364,14 +313,6 @@ if [[ $SUFFIX == "wopr" || $SUFFIX == "nrx" ]]; then sed s/Experimental/Operational/1 ${tmp_index} > ${new_index} fi -#if [[ ! -s ${new_index} ]]; then -# if [[ -s ${tmp_index} ]]; then -# $NCP ${tmp_index} ${new_index} -# else -# $NCP ${index_file} ${new_index} -# fi -#fi - if [[ -s ${tmp_index} ]]; then rm ${tmp_index} fi @@ -440,7 +381,6 @@ done img_list=`ls *.png` # rm any images for sources not in $SATYPE for img in ${img_list}; do tmp=`echo "$img" | cut -d. -f1` - echo $tmp img_match=`echo $SATYPE | grep $tmp` if [[ ${#img_match} -le 0 ]]; then rm -f ${img} @@ -448,22 +388,24 @@ for img in ${img_list}; do done -#--------------------------------------------------- -# if on wcoss then cd $imgndir and do the rsync here +#------------------------------------------------------ +# if on wcoss_d then cd $imgndir and do the rsync here # -if [[ $MY_MACHINE = "wcoss" ]]; then +new_webdir=${WEBDIR}/regional/${SUFFIX} + +if [[ $MY_MACHINE = "wcoss_d" || $MY_MACHINE = "wcoss2" ]]; then if [[ ${imgndir} != "/" ]]; then # sanity check to avoid serious embarrassment /usr/bin/rsync -ave ssh --exclude *.ctl.${Z} ${imgndir}/ \ - ${WEB_USER}@${WEB_SVR}.ncep.noaa.gov:${WEBDIR}/ + ${WEB_USER}@${WEB_SVR}.ncep.noaa.gov:${new_webdir}/ fi fi #------------------------ # clean up $workdir # -#cd $workdir -#cd ../ -#rm -rf $workdir +cd $workdir +cd ../ +rm -rf $workdir echo "" echo "END install_rgn.sh" diff --git a/util/Radiance_Monitor/image_gen/html/menu.html.glb b/util/Radiance_Monitor/image_gen/html/menu.html.glb deleted file mode 100755 index 40a6a6fe45..0000000000 --- a/util/Radiance_Monitor/image_gen/html/menu.html.glb +++ /dev/null @@ -1,53 +0,0 @@ - - - - - Radiance Data Monitoring - - - - - - - -

Experimental GDAS

-

Monitoring Plots

- - - - - - - -
- GDAS home page
- Radiance monitoring home page
- - Time and Angle dependent plots
- Horizontal data plots
- Bias correction coefficients
- Bias correction terms
- Summary plots
- Comparison plots
- -
- - diff --git a/util/Radiance_Monitor/image_gen/html/menu.html.rgn b/util/Radiance_Monitor/image_gen/html/menu.html.rgn deleted file mode 100644 index 95ae4be968..0000000000 --- a/util/Radiance_Monitor/image_gen/html/menu.html.rgn +++ /dev/null @@ -1,53 +0,0 @@ - - - - - Radiance Data Monitoring - - - - - - - -

Experimental NDAS

-

Monitoring Plots

- - - - - - - -
- GDAS home page
- Radiance monitoring home page
- - Time and Angle dependent plots
- Horizontal data plots
- Bias correction coefficients
- Bias correction terms
- Summary plots
- Comparison plots
- -
- - diff --git a/util/Radiance_Monitor/image_gen/html/mk_intro.sh b/util/Radiance_Monitor/image_gen/html/mk_intro.sh index 0e2ae95301..9274395e6c 100755 --- a/util/Radiance_Monitor/image_gen/html/mk_intro.sh +++ b/util/Radiance_Monitor/image_gen/html/mk_intro.sh @@ -111,7 +111,6 @@ while read line; do echo ' ' >> $table # end table row if [[ $ctr != $num_sats ]]; then -# row_end=`expr $row_end + $num_rows` row_end=`expr $row_end + $cols` if [[ $extra < $remainder ]]; then row_end=`expr $row_end + 1` @@ -125,54 +124,20 @@ while read line; do done < "$SORTED_LIST" -#echo ' ' >> $table - -#---------------------------------------------------------------- -# process stock intro.html file and add in the customized table -# -table_start=0 -table_end=0 -table_dump=0 #-------------------------------------------------------------- -# Edit the html files to add the platform table to each. +# Edit infile to add the platform table. # -echo "processing ${infile}" - - # copy the $file from start to tag - sed -e '/
/,$d' ${infile} > ${outfile} - - # add the $PLATFORM_TBL (built above) - `cat $table >> ${outfile}` - # copy the $file from 'END_TABLE_INSERT' comment to end - sed -n '/<\/table>/,$p' ${infile} >> ${outfile} +# copy the $file from start to
tag +sed -e '/
/,$d' ${infile} > ${outfile} +# add the $PLATFORM_TBL (built above) +`cat $table >> ${outfile}` -#while read line; do -# echo $line -# test_start=`echo $line | grep "
"` -# test_end=`echo $line | grep "
"` -## echo "test_start len = ${#test_start} -## echo "test_end len = ${#test_end} -# -# if [[ $table_start -eq 0 && ${#test_start} -gt 0 ]]; then -# table_start=1 -# fi -# -# if [[ ${table_start} -eq 0 || ${table_end} -eq 1 ]]; then -# echo $line >> $outfile -# elif [[ ${table_start} -eq 1 && ${table_dump} -eq 0 ]]; then -# cat $table >> $outfile -# table_dump=1 -# fi -# -# if [[ $table_end -eq 0 && ${#test_end} -gt 0 ]]; then -# table_end=1 -# fi -# -#done < $infile - +# copy the $file from 'END_TABLE_INSERT' comment to end +sed -n '/<\/table>/,$p' ${infile} >> ${outfile} +rm ${infile} echo END mk_intro.sh echo diff --git a/util/Radiance_Monitor/image_gen/html/sum_thumbs.tar b/util/Radiance_Monitor/image_gen/html/sum_thumbs.tar new file mode 100644 index 0000000000000000000000000000000000000000..95d6498ba093c4883755062fe1ff789b0d4270b2 GIT binary patch literal 204800 zcmeFYWl&r}v@SXjk^n<+*Fa#9ATzi_a2Z^J!{Dw1A-KB_GB^Q(1Shx$cL)SY2oAv^ zcp%7~oLlwIId$Lrck8}-Rj;R}yLa1OtGoBwYkjLX%-qh@2Exzj;pyZAbAQk2>TCn} z-yOiHz|YGI0E4-C!Tb2HuJQ%{q&j$v)1pn_Q zU=o~Xs!Hv<5KkpKXhYfhVnD5~O_xw5SOSr66-G2_NgFPTep=cFuzD(&J-P=V7+=Nw#x}VWkop2`0^Y@73Bc zd+{gEs|(msP$YnTT&1-YOP;}gsp%7&41Ds^pH7xw=chqMYSx6cR8-fAK=um;hC=*a zp=N`s(fWpYT^*94)5}Y~yTB*tfc>_WBd-j3g+Ct?In{Mn{~_h}H) zp*l_bPLsg9KN^?>vR~I-k*cK{#8~tp4f^T&4ipwLa6xK(0zk(#0Gwb!ESs1Z2a(ZM z3O{k9eq53nt#>4k=vz-Jv0G}Q6ZgZH0!S`+N?|fkM!@RyFnF-(8|LHcst@o_H#Q0n z+Oo5=M=T3d*~ZQ?Oatt)fjcHqN(WeYjlIRSD=B6+EVB~`^_zIIEcdn{Hwn|n>l(!%?`hz96+vJNr{GVx2 zzs_&gU0*U-;8;MFSi_<=Tbv95;K7$0j{F{zbDr^0^9GZx_~xxf=#pk@!@LX@hbpfr zEMf}DBhUcFse?$#*S~|KB-g%Da0f3Od5I*v0^n>)zJk_|xVLt6cjr%78!x0;Q8;qv zPd!m_s7mENY&-IzYRS-4Vx_QnpBHI7N^TKHgZuo%UxQdv=|qe{iFF8evG#)858sh` z;_sMf0y{n$HU`=n8bz{=U+~J`_V)Hd7Qr1&rJbL1b2`Mb>sHUJFRUOJ{wMfgTNc~V zApPQzMuVp#bn$n5S`sfsAO8N?92@$*9Llx)0j#nRX@!QlZ~UHSJ8q=ud4$oYVg1-$ zwCx*iqQrz?yq+X9DWKE#zBRAc)MvTzU^-TG0NTro(8o%F3<)lZyqMblX z(jL=#jwp_;!C=ExTgcXXym48YzsS71c09-c&vgUjP|kGd1BZysX^q^zW7wU^Cp(gHv zm+g4-P-fxHaK7mruxyKVUncKI@?PSbxmNx?Jp2?xV>&d$tn=%x?k!=*6Pem)5LHIc6B=bwa^{GI?osl?6dROndMFNT&`}dcDu~Q9xi2zbZ=?TeDrlBBKZ6bbx zFUE8$?pru|&vt_;mokhkI=9g+J~A`JLVP|E|EO#xn`Khg4UR-S;-68U-RW3G*CMa; zJa7TRU(XMKv?{a;ek!wj7gA)8Ht^1k+mgSgh@_@3OZh1(^LJBDj%W4`X;W=aMemUVe|tNO`n^?yn0J ziIepMiznFHeLKf&StC~CX_j$8gMG(|iQV^H)1RB@R3ua}0}Y@TqNjuO1wI+)EKj#2 zcjF-w)iv|YP9Cjy6=vvwt2yiBk46;&_gAkO6j>YRCtHJJ&wzd zVG&?Ypb$QhMvUXPy(Qip0N>_ycZ^Rbte*blHjD zlN!HWeVbZ5Tv@JJ{8Si&77O9@tJ@l@TKdJ_Y$;E^15Wd8zIEPi|8zo?M7U3B&BLOs zK=|8-)>JzEYk-8pvjJSe#}E%eN$SvRGaO3_Vp;zC>OFCnfWI7b!aV8L<}ID!yia#} z$Q86$u#uBZn-OHbBU3XIoLME1~5&jurCe@B($qt5n=zl z$ps|Wq3Gh7bpEL_<@2PZD4{2L-G5G9hkmDk@9=s}%!ywxGst858l2nf^I~2~KO|M2 zdM80}zE`jpJna)dFD1fFa+K>LUQ@c4dV*>d(l~uBnEmx$%)q9l^Us5XW+QZy$rX2t zz8QcKCnhxf_;>oU^IZa<7=L(CFd+^y*}CS(`z&e~l8vzz$67Vy#yiOJk>TF*$+iFL z<-H~My33Bc%NLDeH|p#XtV`L|>%%D`^q|$i$4++Tll$-AGd1nomkE>SZx0U-udd`U z(Smez!bfE>-TeKH7d)4mX7_lvN>+QgUIuwQm%L~S&6aRR?QEGP`bLP z7hZmW!Ac%xLoaDbF@?hBZ#j~szTFUs2QBh$mM&KhPCmK$@{8?0pi$}7QvmJmYjc`e zRY^(mpWcQ7YyCG?{1=N&g(E%M9eH)jw%pKF*^K6k@+zZIb3eF?gOh%`mA;9INyLwH zCOzr6&1N5_FU!{N+V9~3t)ltWApw^PoROT|u4WB>(`?{SfeuyXc(f4cHl;>t>bIzUh&QD$T0Idz2XFinK68;qSK z?XO#!*egQu_|N$j%U-6F7Z*N|#KzxvkIwwrf(d7S?hKC^#kVY)QP`1vA6Q~O3<$=5 zGu4#lg8Rj?-qxaQS|I$~dC`l)< zt)FgLhuNhDXXGP)_x37b+$wJY+!<|birh})m0~}m4Nl#`vifDT5;@4%C6KL{ymH2m} z>@bjhN6lkI0!{w4tcLbNqelYm8r3RZR1R7qd57`Gz^BDCHe*1BGmM?s>!!cFh-h?K zmxUx#)2GsZst<&}HO(9B{IrYQg;TA~9-W4dYHbaLuioBMlrwiQpanNhk`u{l2s-5_ zQx2Wp-U=(G2*{!bO;7kbR^0)A;WoXP_%Y;VX~rVjt0;x9``)i%9z6(uL!+auy5f@o zXNKnJTcq5blD9VlfZ)=9wdqJ=LO1VAK=cWFu`f5#Lnm4yVAz$gWDmx zz{-O*yDXNN$*HN(+vcUCiHZ36+eo^Amd;M5!gnVt+jx!FPqiV#F9S~ zD+Y(;-X2xFCXMV%Ex&MN#WSoYZLrqIZs!ofLqvP`1D$)wvnQ-q^O_Q5`pg*mjmBo9 z9K_J@eJ5n)Hhl`Cx9-;IcI+17cRL&%#ns;UeV14vdRbq}nDhCa2wy;|gh|bn$SU<8 z;zh^mlVj@FnS!!babNqoZsn5M9GFvpJ(q&2f?R)-F$p>FJn7<1lvWbGT1SOM>Q^sH z#4b*GKgneeU=xgb$%GkUF=>~pyoGD9(($SB%Zo}NPwkkN!=m9Nz60eZ}gNEMl6mT>qr>O}KuF3!PIb9*OB571f^PzDURcl|8*C3kRgjMq?oebe8CVmgGmc1hEetH<_B+1_at7iEpJ%jpb5?=_D_EP3ZbL=JqF zJ5|W(sC>&JzVrUjxu4_C5{-vRjBIW~Ho}up>g{#DKwrTQ8A9R4$%dMHdrkeiE?sqW zrY*+}eK&WkljNQy3_7(@LY`tIS1!0bY;ls+{>g>iEkR!{3XB7F9m!Gepo_}M2d$u; z4^v1J|NM}YV*X;Ft9;A{eev2h7|4X&!VQ0t;C~P~q4s!QGyDB;?#Na2?6wQftwS^w zF(RsJKKVrmo)>aX5CmDPDEHHeI(M}HO=Q#1nIMbfEJpuaSy zraONoAWVK65ahrHCyGyi_6aJH*fsU?X^>VIX4bO^H0s`gr;0c8V?}#i zQkS+LVdRyKGTOwl$j*EHH(l8%nhPOrMy&w#Q|`>EYF{(NV@>l=>Rj=85k-%#T+;Rk=_jr(&< zqkvYrJm2KB)P$;9ZhsfZCi-&2qV;Oktt08H=MxsJ5W3UH^OMLEu31t4!i@=^*}7Zv zfcK^k^I7JfNO?wG5b0P1Fh<&=&#NaxZd*Cse8&N?*|%Q#?95PC04CzTnejc3r%dgn zQw8o|zc=k`gYn=fxr4yV?t^3hyvh6xEl(t=NN2}G^~qZudfTkmZ&#|Sa1gP8xH|T+ zEN4cai{L;LzWAz}dyhu)h}aXZ1Umr+bZ@mNdJJ*qjp(@ZTY681LZ4L^ABR1#RrA}C z?KU-Op*E)A{+rQJ_3#Y6qfvYKk8th!lh_7dd~4jSC``b$ZYAqT9xZ| z@9)w#@CJS_n{wiAS1%IVRJ8Imx%^d`b?UL*TDXCy%fZGyJD@VPT>fA2-g@M+sLCw0 zyMvFH2cYGC64e!Fv`cel?i{uY0;^MG`7wq2L;8hk3-LvpWxALObl6vl58A~&_l$(0 zjT39Mf7?v!(mn??)DCNBVgRMw#0orrRZ@9WNpyE>PNbBEX6T$QLfA%y-LZ-FeA!u(@n5~hMo`E&)`U5n&#v05UBWB==2UAu(Z zj6NrtTXO8geHmdu*KpTG0Gf^YseG41OlXE8=v9OxflXw|8Ai^!^o!qo-kNjIOtd&~ z*67-uGb{8?-!?r#dx`&Msn5|z6uU7gSYote=RYEeI$oU(pxu@5m5?~P<= zVRJ2IxY(aI>aOo4MlJHH=%-iAu{Cs?b`d(41t1irUK+*66%LbYEnNDEKM$&BFJjSD zdeLVVH&P!LHCjJ@o&(WMXW9`(fpT|#O4GikIy%TAoIZR`f(Vgd#E6W;GD_pWR!kwc z_|cR*VXZmZb<#KfXvlu2_RL2I6`lgt6U=OfP*9HRT-`xiTQypb5SjA{dXG;hXAC-^ zx8aWBjm7w@U)G-(&aA{ ze8mKp)q`8Y;8wyG?pCM-zzv4*azMZwU_ouLfH04MFrUzCFjyE2KFHM5_-`Fh(4n36 z`~T+-S^aBAs15){Syh>8DYKA&`TrvU_8<7K!2gT*FApCt3KRan@Lyi=f8f8o{Qp1x z`+s5s{IB>g1{45zqKWYYb@zB}xkFS!r4H4WeSCaGWq`1*C*W38^fw>` zFgy&?2|;zS*{D{$rA=pr0}_5IIlhyHszITA9BPT3l8_pRB-07%ju_pNS@D|1ake zT0^cVZgs