module crtm_interface !$$$ module documentation block ! . . . ! module: crtm_interface module for setuprad. Calculates profile and calls crtm ! prgmmr: ! ! abstract: crtm_interface module for setuprad. Initializes CRTM, Calculates profile and ! calls CRTM and destroys initialization ! ! program history log: ! 2010-08-17 Derber - initial creation from intrppx ! 2011-05-06 merkova/todling - add use of q-clear calculation for AIRS ! 2011-04-08 li - (1) Add nst_gsi, itref,idtw, idtc, itz_tr to apply NSST. ! - (2) Use Tz instead of Ts as water surface temperature when nst_gsi > 1 ! - (3) add tzbgr as one of the out dummy variable ! - (4) Include tz_tr in ts calculation over water ! - (5) Change minmum temperature of water surface from 270.0 to 271.0 ! 2011-07-04 todling - fixes to run either single or double precision ! 2011-09-20 hclin - modified for modis_aod ! (1) The jacobian of wrfchem/gocart p25 species (not calculated in CRTM) ! is derived from dust1 and dust2 ! (2) skip loading geometry and surface structures for modis_aod ! (3) separate jacobian calculation for modis_aod ! 2012-01-17 sienkiewicz - pass date to crtm for SSU cell pressure ! 2013-02-25 zhu - add cold_start option for regional applications ! 2013-10-19 todling - metguess now holds background ! 2013-11-16 todling - merge in latest DTC AOD development; ! revisit handling of green-house-gases ! 2014-01-01 li - change the protection of data_s(itz_tr) ! 2014-02-26 zhu - add non zero jacobian ! 2014-04-27 eliu - add call crtm_forward to calculate clear-sky Tb under all-sky condition ! 2015-09-10 zhu - generalize enabling all-sky and using aerosol (radiance_mod & radmod) in radiance ! assimilation. use n_clouds_jac,cloud_names_jac,n_aerosols_jac,aerosol_names_jac, ! n_clouds_fwd,cloud_names_fwd, etc for difference sensors and channels ! - add handling of mixed_use of channels in a sensor (some are clear-sky, others all-sky) ! 2016-06-03 collard - Added changes to allow for historical naming conventions ! 2017-02-24 zhu/todling - remove gmao cloud fraction treatment ! 2018-01-12 collard - Force all satellite and solar zenith angles to be >= 0. ! ! ! subroutines included: ! sub init_crtm ! sub call_crtm ! sub destroy_crtm ! ! attributes: ! language: f90 ! machine: ! !$$$ end documentation block use kinds,only: r_kind,i_kind,r_single use crtm_module, only: crtm_atmosphere_type,crtm_surface_type,crtm_geometry_type, & crtm_options_type,crtm_rtsolution_type,crtm_destroy,crtm_options_destroy, & crtm_options_create,crtm_options_associated,success,crtm_atmosphere_create, & crtm_surface_create,crtm_k_matrix,crtm_forward, & ssu_input_setvalue, & crtm_channelinfo_type, & crtm_surface_destroy, crtm_surface_associated, crtm_surface_zero, & crtm_atmosphere_associated, & crtm_atmosphere_destroy,crtm_atmosphere_zero, & crtm_rtsolution_type, crtm_rtsolution_create, & crtm_rtsolution_destroy, crtm_rtsolution_associated, & crtm_irlandcoeff_classification, & crtm_kind => fp, & crtm_microwave_sensor => microwave_sensor use gridmod, only: lat2,lon2,nsig,msig,nvege_type,regional,wrf_mass_regional,netcdf,use_gfs_ozone use mpeu_util, only: die use crtm_aod_module, only: crtm_aod_k use radiance_mod, only: n_actual_clouds,cloud_names,n_clouds_fwd,cloud_names_fwd, & n_clouds_jac,cloud_names_jac,n_actual_aerosols,aerosol_names,n_aerosols_fwd,aerosol_names_fwd, & n_aerosols_jac,aerosol_names_jac,rad_obs_type,cw_cv implicit none private public init_crtm ! Subroutine initializes crtm for specified instrument public call_crtm ! Subroutine creates profile for crtm, calls crtm, then adjoint of create public destroy_crtm ! Subroutine destroys initialization for crtm public sensorindex public surface public isatid ! = 1 index of satellite id public itime ! = 2 index of analysis relative obs time public ilon ! = 3 index of grid relative obs location (x) public ilat ! = 4 index of grid relative obs location (y) public ilzen_ang ! = 5 index of local (satellite) zenith angle (radians) public ilazi_ang ! = 6 index of local (satellite) azimuth angle (radians) public iscan_ang ! = 7 index of scan (look) angle (radians) public iscan_pos ! = 8 index of integer scan position public iszen_ang ! = 9 index of solar zenith angle (degrees) public isazi_ang ! = 10 index of solar azimuth angle (degrees) public ifrac_sea ! = 11 index of ocean percentage public ifrac_lnd ! = 12 index of land percentage public ifrac_ice ! = 13 index of ice percentage public ifrac_sno ! = 14 index of snow percentage public its_sea ! = 15 index of ocean temperature public its_lnd ! = 16 index of land temperature public its_ice ! = 17 index of ice temperature public its_sno ! = 18 index of snow temperature public itsavg ! = 19 index of average temperature public ivty ! = 20 index of vegetation type public ivfr ! = 21 index of vegetation fraction public isty ! = 22 index of soil type public istp ! = 23 index of soil temperature public ism ! = 24 index of soil moisture public isn ! = 25 index of snow depth public izz ! = 26 index of surface height public idomsfc ! = 27 index of dominate surface type public isfcr ! = 28 index of surface roughness public iff10 ! = 29 index of ten meter wind factor public ilone ! = 30 index of earth relative longitude (degrees) public ilate ! = 31 index of earth relative latitude (degrees) public iclr_sky ! = 7 index of clear sky amount (goes_img, seviri, abi) public isst_navy ! = 7 index of navy sst retrieval (K) (avhrr_navy) public idata_type ! = 32 index of data type (151=day, 152=night, avhrr_navy) public iclavr ! = 32 index of clavr cloud flag (avhrr) public isst_hires ! = 33 index of interpolated hires sst public itref ! = 34/36 index of Tr public idtw ! = 35/37 index of d(Tw) public idtc ! = 36/38 index of d(Tc) public itz_tr ! = 37/39 index of d(Tz)/d(Tr) ! Note other module variables are only used within this routine character(len=*), parameter :: myname='crtm_interface' ! Indices for the CRTM NPOESS EmisCoeff file integer(i_kind), parameter :: INVALID_LAND = 0 integer(i_kind), parameter :: COMPACTED_SOIL = 1 integer(i_kind), parameter :: TILLED_SOIL = 2 integer(i_kind), parameter :: IRRIGATED_LOW_VEGETATION = 5 integer(i_kind), parameter :: MEADOW_GRASS = 6 integer(i_kind), parameter :: SCRUB = 7 integer(i_kind), parameter :: BROADLEAF_FOREST = 8 integer(i_kind), parameter :: PINE_FOREST = 9 integer(i_kind), parameter :: TUNDRA = 10 integer(i_kind), parameter :: GRASS_SOIL = 11 integer(i_kind), parameter :: BROADLEAF_PINE_FOREST = 12 integer(i_kind), parameter :: GRASS_SCRUB = 13 integer(i_kind), parameter :: URBAN_CONCRETE = 15 integer(i_kind), parameter :: BROADLEAF_BRUSH = 17 integer(i_kind), parameter :: WET_SOIL = 18 integer(i_kind), parameter :: SCRUB_SOIL = 19 real(r_kind) , save ,allocatable,dimension(:,:) :: aero ! aerosol (guess) profiles at obs location real(r_kind) , save ,allocatable,dimension(:,:) :: aero_conc ! aerosol (guess) concentrations at obs location real(r_kind) , save ,allocatable,dimension(:) :: auxrh ! temporary array for rh profile as seen by CRTM character(len=20),save,allocatable,dimension(:) :: ghg_names ! names of green-house gases integer(i_kind), save ,allocatable,dimension(:) :: icloud ! cloud index for those considered here integer(i_kind), save ,allocatable,dimension(:) :: jcloud ! cloud index for those fed to CRTM real(r_kind) , save ,allocatable,dimension(:,:) :: cloud ! cloud considered here real(r_kind) , save ,allocatable,dimension(:,:) :: cloudefr ! effective radius of cloud type in CRTM real(r_kind) , save ,allocatable,dimension(:,:) :: cloud_cont ! cloud content fed into CRTM real(r_kind) , save ,allocatable,dimension(:,:) :: cloud_efr ! effective radius of cloud type in CRTM real(r_kind) , save ,allocatable,dimension(:,:,:,:) :: gesqsat ! qsat to calc rh for aero particle size estimate real(r_kind) , save ,allocatable,dimension(:) :: lcloud4crtm_wk ! cloud info usage index for each channel integer(i_kind),save, allocatable,dimension(:) :: map_to_crtm_ir integer(i_kind),save, allocatable,dimension(:) :: map_to_crtm_mwave integer(i_kind),save, allocatable,dimension(:) :: icw integer(i_kind),save, allocatable,dimension(:) :: iaero_jac integer(i_kind),save :: isatid,itime,ilon,ilat,ilzen_ang,ilazi_ang,iscan_ang integer(i_kind),save :: iscan_pos,iszen_ang,isazi_ang,ifrac_sea,ifrac_lnd,ifrac_ice integer(i_kind),save :: ifrac_sno,its_sea,its_lnd,its_ice,its_sno,itsavg integer(i_kind),save :: ivty,ivfr,isty,istp,ism,isn,izz,idomsfc,isfcr,iff10,ilone,ilate integer(i_kind),save :: iclr_sky,isst_navy,idata_type,isst_hires,iclavr integer(i_kind),save :: itref,idtw,idtc,itz_tr,istype integer(i_kind),save :: sensorindex integer(i_kind),save :: ico2,ico24crtm integer(i_kind),save :: n_actual_aerosols_wk ! number of aerosols considered integer(i_kind),save :: n_aerosols_fwd_wk ! number of aerosols considered integer(i_kind),save :: n_aerosols_jac_wk ! number of aerosols considered integer(i_kind),save :: n_actual_clouds_wk ! number of clouds considered integer(i_kind),save :: n_clouds_fwd_wk ! number of clouds considered integer(i_kind),save :: n_clouds_jac_wk ! number of clouds considered integer(i_kind),save :: n_ghg ! number of green-house gases integer(i_kind),save :: itv,iqv,ioz,ius,ivs,isst integer(i_kind),save :: indx_p25, indx_dust1, indx_dust2 logical ,save :: lwind logical ,save :: cld_sea_only_wk logical ,save :: mixed_use integer(i_kind), parameter :: min_n_absorbers = 2 type(crtm_atmosphere_type),save,dimension(1) :: atmosphere type(crtm_surface_type),save,dimension(1) :: surface type(crtm_geometry_type),save,dimension(1) :: geometryinfo type(crtm_options_type),save,dimension(1) :: options type(crtm_channelinfo_type),save,dimension(1) :: channelinfo type(crtm_atmosphere_type),save,allocatable,dimension(:,:):: atmosphere_k type(crtm_atmosphere_type),save,allocatable,dimension(:,:):: atmosphere_k_clr type(crtm_surface_type),save,allocatable,dimension(:,:):: surface_k type(crtm_surface_type),save,allocatable,dimension(:,:):: surface_k_clr type(crtm_rtsolution_type),save,allocatable,dimension(:,:):: rtsolution type(crtm_rtsolution_type),save,allocatable,dimension(:,:):: rtsolution0 type(crtm_rtsolution_type),save,allocatable,dimension(:,:):: rtsolution_clr type(crtm_rtsolution_type),save,allocatable,dimension(:,:):: rtsolution_k type(crtm_rtsolution_type),save,allocatable,dimension(:,:):: rtsolution_k_clr ! Mapping land surface type of GFS to CRTM ! Notes: index 0 is water, and index 13 is ice. The two indices are not ! used and just assigned to COMPACTED_SOIL. Also, since there ! is currently one relevant mapping for the global we apply ! 'crtm' in the naming convention. integer(i_kind), parameter, dimension(0:13) :: gfs_to_crtm=(/COMPACTED_SOIL, & BROADLEAF_FOREST, BROADLEAF_FOREST, BROADLEAF_PINE_FOREST, PINE_FOREST, & PINE_FOREST, BROADLEAF_BRUSH, SCRUB, SCRUB, SCRUB_SOIL, TUNDRA, & COMPACTED_SOIL, TILLED_SOIL, COMPACTED_SOIL/) ! Mapping surface classification to CRTM integer(i_kind), parameter :: USGS_N_TYPES = 24 integer(i_kind), parameter :: IGBP_N_TYPES = 20 integer(i_kind), parameter :: GFS_N_TYPES = 13 integer(i_kind), parameter :: SOIL_N_TYPES = 16 integer(i_kind), parameter :: GFS_SOIL_N_TYPES = 9 integer(i_kind), parameter :: GFS_VEGETATION_N_TYPES = 13 integer(i_kind), parameter, dimension(1:USGS_N_TYPES) :: usgs_to_npoess=(/URBAN_CONCRETE, & COMPACTED_SOIL, IRRIGATED_LOW_VEGETATION, GRASS_SOIL, MEADOW_GRASS, & MEADOW_GRASS, MEADOW_GRASS, SCRUB, GRASS_SCRUB, MEADOW_GRASS, & BROADLEAF_FOREST, PINE_FOREST, BROADLEAF_FOREST, PINE_FOREST, & BROADLEAF_PINE_FOREST, INVALID_LAND, WET_SOIL, WET_SOIL, & IRRIGATED_LOW_VEGETATION, TUNDRA, TUNDRA, TUNDRA, TUNDRA, & INVALID_LAND/) integer(i_kind), parameter, dimension(1:IGBP_N_TYPES) :: igbp_to_npoess=(/PINE_FOREST, & BROADLEAF_FOREST, PINE_FOREST, BROADLEAF_FOREST, BROADLEAF_PINE_FOREST, & SCRUB, SCRUB_SOIL, BROADLEAF_BRUSH, BROADLEAF_BRUSH, SCRUB, BROADLEAF_BRUSH, & TILLED_SOIL, URBAN_CONCRETE, TILLED_SOIL, INVALID_LAND, COMPACTED_SOIL, & INVALID_LAND, TUNDRA, TUNDRA, TUNDRA/) integer(i_kind), parameter, dimension(1:USGS_N_TYPES) :: usgs_to_usgs=(/1, & 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, & 20, 21, 22, 23, 24/) integer(i_kind), parameter, dimension(1:IGBP_N_TYPES) :: igbp_to_igbp=(/1, & 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, & 20/) integer(i_kind), parameter, dimension(1:IGBP_N_TYPES) :: igbp_to_gfs=(/4, & 1, 5, 2, 3, 8, 9, 6, 6, 7, 8, 12, 7, 12, 13, 11, 0, 10, 10, 11/) integer(i_kind), parameter, dimension(1:USGS_N_TYPES) :: usgs_to_gfs=(/7, & 12, 12, 12, 12, 12, 7, 9, 8, 6, 2, 5, 1, 4, 3, 0, 8, 8, 11, 10, 10, & 10, 11, 13/) ! Mapping soil types to CRTM ! The CRTM soil types for microwave calculations are based on the ! GFS use of the 9 category Zobler dataset. The regional soil types ! are based on a 16 category representation of FAO/STATSGO. integer(i_kind), parameter, dimension(1:SOIL_N_TYPES) :: map_soil_to_crtm=(/1, & 1, 4, 2, 2, 8, 7, 2, 6, 5, 2, 3, 8, 1, 6, 9/) contains subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmod) !$$$ subprogram documentation block ! . . . . ! subprogram: init_crtm initializes things for use with call to crtm from setuprad ! ! prgmmr: derber org: np2 date: 2010-08-17 ! ! abstract: initialize things for use with call to crtm from setuprad. ! ! program history log: ! 2010-08-17 derber ! 2011-02-16 todling - add calculation of rh when aerosols are available ! 2011-05-03 todling - merge with Min-Jeong's MW cloudy radiance; combine w/ metguess ! 2011-05-20 mccarty - add atms wmo_sat_id hack (currently commented out) ! 2011-07-20 zhu - modified codes for lcw4crtm ! 2012-03-12 yang - modify to use ch4,n2o,and co ! 2012-12-03 eliu - add logic for RH total ! 2014-01-31 mkim - add flexibility in the variable lcw4crtm for the case when ql and ! qi are separate control variables for all-sky MW radiance DA ! 2014-04-27 eliu - add capability to call CRTM forward model to calculate ! clear-sky Tb under all-sky condition ! 2015-09-20 zhu - use centralized radiance info from radiance_mod: rad_obs_type, ! n_clouds_jac,cloud_names_jac,n_aerosols_jac,aerosol_names_jac,etc ! 2015-09-04 J.Jung - Added mods for CrIS full spectral resolution (FSR) and ! CRTM subset code for CrIS. ! 2019-04-25 H.Liu - Add nreal into the namelist to allow flexible SST variable index assigned ! ! input argument list: ! init_pass - state of "setup" processing ! mype_diaghdr - processor to produce output from crtm ! mype - current processor ! nchanl - number of channels ! isis - instrument/sensor character string ! obstype - observation type ! nreal - number of descriptor information in data_s ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use gsi_bundlemod, only: gsi_bundlegetpointer use gsi_chemguess_mod, only: gsi_chemguess_bundle ! for now, a common block use gsi_chemguess_mod, only: gsi_chemguess_get use gsi_metguess_mod, only: gsi_metguess_bundle ! for now, a common block use gsi_metguess_mod, only: gsi_metguess_get use crtm_module, only: mass_mixing_ratio_units,co2_id,o3_id,crtm_init, & crtm_channelinfo_subset, crtm_channelinfo_n_channels, toa_pressure,max_n_layers, & volume_mixing_ratio_units,h2o_id,ch4_id,n2o_id,co_id use radinfo, only: crtm_coeffs_path use radinfo, only: radjacindxs,radjacnames,jpch_rad,nusis,nuchan use aeroinfo, only: aerojacindxs use guess_grids, only: ges_tsen,ges_prsl,nfldsig use mpeu_util, only: getindex use constants, only: zero,max_varname_length use gsi_io, only: verbose implicit none ! argument logical ,intent(in) :: init_pass integer(i_kind),intent(in) :: nchanl,mype_diaghdr,mype,nreal character(20) ,intent(in) :: isis character(10) ,intent(in) :: obstype type(rad_obs_type),intent(in) :: radmod ! local parameters character(len=*), parameter :: myname_=myname//'*init_crtm' ! local variables integer(i_kind) :: ier,ii,error_status,iderivative integer(i_kind) :: k, subset_start, subset_end logical :: ice,Load_AerosolCoeff,Load_CloudCoeff character(len=20),dimension(1) :: sensorlist integer(i_kind) :: indx,iii,icloud4crtm ! ...all "additional absorber" variables integer(i_kind) :: j,icount integer(i_kind) :: ig integer(i_kind) :: n_absorbers logical quiet logical print_verbose print_verbose=.false. if(verbose)print_verbose=.true. isst=-1 ivs=-1 ius=-1 ioz=-1 iqv=-1 itv=-1 ! Get indexes of variables composing the jacobian indx =getindex(radjacnames,'tv') if(indx>0) itv=radjacindxs(indx) indx =getindex(radjacnames,'q' ) if(indx>0) iqv=radjacindxs(indx) indx =getindex(radjacnames,'oz') if(indx>0) ioz=radjacindxs(indx) indx =getindex(radjacnames,'u') if(indx>0) ius=radjacindxs(indx) indx =getindex(radjacnames,'v') if(indx>0) ivs=radjacindxs(indx) lwind=ius>0.and.ivs>0 indx=getindex(radjacnames,'sst') if(indx>0) isst=radjacindxs(indx) ! Get indexes of variables for cloud jacobians if (n_clouds_jac>0) then allocate(icw(max(n_clouds_jac,1))) icw=-1 icount=0 do ii=1,n_clouds_jac indx=getindex(radjacnames,trim(cloud_names_jac(ii))) if (indx>0) then icount=icount+1 icw(icount)=radjacindxs(indx) end if end do end if ! Get indexes of variables composing the jacobian_aero if (n_actual_aerosols > 0) then indx_p25 = getindex(aerosol_names,'p25') indx_dust1 = getindex(aerosol_names,'dust1') indx_dust2 = getindex(aerosol_names,'dust2') if (n_aerosols_jac >0) then allocate(iaero_jac(n_aerosols_jac)) iaero_jac=-1 icount=0 do ii=1,n_actual_aerosols indx=getindex(aerosol_names_jac,trim(aerosol_names(ii))) if(indx>0) then icount=icount+1 iaero_jac(icount)=aerojacindxs(indx) endif end do endif endif ! When Cloud is available in MetGuess, defined Cloudy Radiance mixed_use=.false. if (radmod%lcloud_fwd) then allocate(lcloud4crtm_wk(radmod%nchannel)) lcloud4crtm_wk(:) = radmod%lcloud4crtm(:) do ii=1,radmod%nchannel if (lcloud4crtm_wk(ii)<0) then mixed_use=.true. exit end if end do allocate(cloud_cont(msig,n_clouds_fwd)) allocate(cloud_efr(msig,n_clouds_fwd)) allocate(jcloud(n_clouds_fwd)) allocate(cloud(nsig,n_clouds_fwd)) allocate(cloudefr(nsig,n_clouds_fwd)) allocate(icloud(n_actual_clouds)) cloud_cont=zero cloud_efr =zero cloud =zero cloudefr =zero call gsi_bundlegetpointer(gsi_metguess_bundle(1),cloud_names,icloud,ier) iii=0 do ii=1,n_actual_clouds call gsi_metguess_get ( 'i4crtm::'//trim(cloud_names(ii)), icloud4crtm, ier ) if (icloud4crtm>10) then iii=iii+1 jcloud(iii)=ii endif end do if(iii/=n_clouds_fwd) call die(myname_,'inconsistent cloud count',1) n_actual_clouds_wk = n_actual_clouds n_clouds_fwd_wk = n_clouds_fwd n_clouds_jac_wk = n_clouds_jac cld_sea_only_wk = radmod%cld_sea_only Load_CloudCoeff = .true. else n_actual_clouds_wk = 0 n_clouds_fwd_wk = 0 n_clouds_jac_wk = 0 cld_sea_only_wk = .false. Load_CloudCoeff = .false. endif ! Set up index for input satellite data array isatid = 1 ! index of satellite id itime = 2 ! index of analysis relative obs time ilon = 3 ! index of grid relative obs location (x) ilat = 4 ! index of grid relative obs location (y) ilzen_ang = 5 ! index of local (satellite) zenith angle (radians) ilazi_ang = 6 ! index of local (satellite) azimuth angle (radians) iscan_ang = 7 ! index of scan (look) angle (radians) iscan_pos = 8 ! index of integer scan position iszen_ang = 9 ! index of solar zenith angle (degrees) isazi_ang = 10 ! index of solar azimuth angle (degrees) ifrac_sea = 11 ! index of ocean percentage ifrac_lnd = 12 ! index of land percentage ifrac_ice = 13 ! index of ice percentage ifrac_sno = 14 ! index of snow percentage its_sea = 15 ! index of ocean temperature its_lnd = 16 ! index of land temperature its_ice = 17 ! index of ice temperature its_sno = 18 ! index of snow temperature itsavg = 19 ! index of average temperature ivty = 20 ! index of vegetation type ivfr = 21 ! index of vegetation fraction isty = 22 ! index of soil type istp = 23 ! index of soil temperature ism = 24 ! index of soil moisture isn = 25 ! index of snow depth izz = 26 ! index of surface height idomsfc = 27 ! index of dominate surface type isfcr = 28 ! index of surface roughness iff10 = 29 ! index of ten meter wind factor ilone = 30 ! index of earth relative longitude (degrees) ilate = 31 ! index of earth relative latitude (degrees) itref = nreal-3 ! index of foundation temperature: Tr idtw = nreal-2 ! index of diurnal warming: d(Tw) at depth zob idtc = nreal-1 ! index of sub-layer cooling: d(Tc) at depth zob itz_tr = nreal ! index of d(Tz)/d(Tr) if (obstype == 'goes_img' .or. obstype == 'abi') then iclr_sky = 7 ! index of clear sky amount elseif (obstype == 'avhrr_navy') then isst_navy = 7 ! index of navy sst (K) retrieval idata_type = 32 ! index of data type (151=day, 152=night) isst_hires = 33 ! index of interpolated hires sst (K) elseif (obstype == 'avhrr') then iclavr = 32 ! index CLAVR cloud flag with AVHRR data isst_hires = 33 ! index of interpolated hires sst (K) elseif (obstype == 'seviri') then iclr_sky = 7 ! index of clear sky amount endif ! get the number of trace gases present in the chemguess bundle n_ghg=0 if(size(gsi_chemguess_bundle)>0) then call gsi_chemguess_get('ghg',n_ghg,ier) if (n_ghg>0) then allocate(ghg_names(n_ghg)) call gsi_chemguess_get('ghg',ghg_names,ier) endif endif n_absorbers = min_n_absorbers + n_ghg ! Are there aerosols to affect CRTM? if (radmod%laerosol_fwd) then if(.not.allocated(aero)) allocate(aero(nsig,n_actual_aerosols)) if(.not.allocated(aero_conc)) allocate(aero_conc(msig,n_actual_aerosols),auxrh(msig)) n_actual_aerosols_wk=n_actual_aerosols n_aerosols_fwd_wk=n_aerosols_fwd n_aerosols_jac_wk=n_aerosols_jac Load_AerosolCoeff=.true. else n_actual_aerosols_wk=0 n_aerosols_fwd_wk=0 n_aerosols_jac_wk=0 Load_AerosolCoeff=.false. endif ! Initialize radiative transfer sensorlist(1)=isis quiet=.not. print_verbose if( crtm_coeffs_path /= "" ) then if(init_pass .and. mype==mype_diaghdr .and. print_verbose) & write(6,*)myname_,': crtm_init() on path "'//trim(crtm_coeffs_path)//'"' error_status = crtm_init(sensorlist,channelinfo,& Process_ID=mype,Output_Process_ID=mype_diaghdr, & Load_CloudCoeff=Load_CloudCoeff,Load_AerosolCoeff=Load_AerosolCoeff, & File_Path = crtm_coeffs_path,quiet=quiet ) else error_status = crtm_init(sensorlist,channelinfo,& Process_ID=mype,Output_Process_ID=mype_diaghdr, & Load_CloudCoeff=Load_CloudCoeff,Load_AerosolCoeff=Load_AerosolCoeff,& quiet=quiet) endif if (error_status /= success) then write(6,*)myname_,': ***ERROR*** crtm_init error_status=',error_status,& ' TERMINATE PROGRAM EXECUTION' call stop2(71) endif sensorindex = 0 if (channelinfo(1)%sensor_id == isis) then sensorindex = 1 if (isis(1:4) == 'iasi' .or. & trim(isis) == 'amsua_aqua' .or. & isis(1:4) == 'airs' .or. & isis(1:4) == 'cris' ) then subset_start = 0 subset_end = 0 do k=1, jpch_rad if (isis == nusis(k)) then if (subset_start == 0) subset_start = k subset_end = k endif end do error_status = crtm_channelinfo_subset(channelinfo(1), & channel_subset = nuchan(subset_start:subset_end)) endif ! This is to try to keep the CrIS naming conventions more flexible. ! The consistency of CRTM and BUFR files is checked in read_cris: else if (channelinfo(1)%sensor_id(1:8) == 'cris-fsr' .AND. isis(1:8) == 'cris-fsr') then sensorindex = 1 subset_start = 0 subset_end = 0 do k=1, jpch_rad if (isis == nusis(k)) then if (subset_start == 0) subset_start = k subset_end = k endif end do error_status = crtm_channelinfo_subset(channelinfo(1), & channel_subset = nuchan(subset_start:subset_end)) else if (channelinfo(1)%sensor_id(1:4) == 'cris' .AND. isis(1:4) == 'cris') then sensorindex = 1 subset_start = 0 subset_end = 0 do k=1, jpch_rad if (isis == nusis(k)) then if (subset_start == 0) subset_start = k subset_end = k endif end do error_status = crtm_channelinfo_subset(channelinfo(1), & channel_subset = nuchan(subset_start:subset_end)) else if (channelinfo(1)%sensor_id(1:4) == 'iasi' .AND. isis(1:4) == 'iasi') then sensorindex = 1 subset_start = 0 subset_end = 0 do k=1, jpch_rad if (isis == nusis(k)) then if (subset_start == 0) subset_start = k subset_end = k endif end do error_status = crtm_channelinfo_subset(channelinfo(1), & channel_subset = nuchan(subset_start:subset_end)) else if (channelinfo(1)%sensor_id(1:4) == 'airs' .AND. isis(1:4) == 'airs') then sensorindex = 1 subset_start = 0 subset_end = 0 do k=1, jpch_rad if (isis == nusis(k)) then if (subset_start == 0) subset_start = k subset_end = k endif end do error_status = crtm_channelinfo_subset(channelinfo(1), & channel_subset = nuchan(subset_start:subset_end)) endif if (sensorindex == 0 ) then write(6,*)myname_,': ***WARNING*** problem with sensorindex=',isis,& ' --> CAN NOT PROCESS isis=',isis,' TERMINATE PROGRAM EXECUTION found ',& channelinfo(1)%sensor_id call stop2(71) endif ! Check for consistency between user specified number of channels (nchanl) ! and those defined by CRTM channelinfo structure. Return to calling ! routine if there is a mismatch. if (nchanl /= crtm_channelinfo_n_channels(channelinfo(sensorindex))) then write(6,*)myname_,': ***WARNING*** mismatch between nchanl=',& nchanl,' and n_channels=',crtm_channelinfo_n_channels(channelinfo(sensorindex)),& ' --> CAN NOT PROCESS isis=',isis,' TERMINATE PROGRAM EXECUTION' call stop2(71) endif ! Allocate structures for radiative transfer if (radmod%lcloud_fwd .and. (.not. mixed_use) .and. (.not. allocated(rtsolution0)) ) & allocate(rtsolution0(channelinfo(sensorindex)%n_channels,1)) allocate(& rtsolution (channelinfo(sensorindex)%n_channels,1),& rtsolution_k(channelinfo(sensorindex)%n_channels,1),& atmosphere_k(channelinfo(sensorindex)%n_channels,1),& surface_k (channelinfo(sensorindex)%n_channels,1)) if (mixed_use) allocate(& rtsolution_clr (channelinfo(sensorindex)%n_channels,1),& rtsolution_k_clr(channelinfo(sensorindex)%n_channels,1),& atmosphere_k_clr(channelinfo(sensorindex)%n_channels,1),& surface_k_clr (channelinfo(sensorindex)%n_channels,1)) ! Check to ensure that number of levels requested does not exceed crtm max if(msig > max_n_layers)then write(6,*) myname_,': msig > max_n_layers - increase crtm max_n_layers ',& msig,max_n_layers call stop2(36) end if ! Create structures for radiative transfer call crtm_atmosphere_create(atmosphere(1),msig,n_absorbers,n_clouds_fwd_wk,n_aerosols_fwd_wk) !_RTod-NOTE if(r_kind==r_single .and. crtm_kind/=r_kind) then ! take care of case: GSI(single); CRTM(double) !_RTod-NOTE call crtm_surface_create(surface(1),channelinfo(sensorindex)%n_channels,tolerance=1.0e-5_crtm_kind) !_RTod-NOTE else !_RTod-NOTE: the following will work in single precision but issue lots of msg and remove more obs than needed if ( channelinfo(sensorindex)%sensor_type == crtm_microwave_sensor ) then call crtm_surface_create(surface(1),channelinfo(sensorindex)%n_channels) if (.NOT.(crtm_surface_associated(surface(1)))) then write(6,*)myname_,' ***ERROR** creating surface.' else surface(1)%sensordata%sensor_id = channelinfo(sensorindex)%sensor_id surface(1)%sensordata%wmo_sensor_id = channelinfo(sensorindex)%wmo_sensor_id surface(1)%sensordata%wmo_satellite_id = channelinfo(sensorindex)%wmo_satellite_id surface(1)%sensordata%sensor_channel = channelinfo(sensorindex)%sensor_channel end if end if !_RTod-NOTE endif if (radmod%lcloud_fwd .and. (.not. mixed_use)) & call crtm_rtsolution_create(rtsolution0,msig) call crtm_rtsolution_create(rtsolution,msig) call crtm_rtsolution_create(rtsolution_k,msig) call crtm_options_create(options,nchanl) if (mixed_use) then call crtm_rtsolution_create(rtsolution_clr,msig) call crtm_rtsolution_create(rtsolution_k_clr,msig) end if if (.NOT.(crtm_atmosphere_associated(atmosphere(1)))) & write(6,*)myname_,' ***ERROR** creating atmosphere.' if (.NOT.(ANY(crtm_rtsolution_associated(rtsolution)))) & write(6,*)myname_,' ***ERROR** creating rtsolution.' if (radmod%lcloud_fwd .and. (.not. mixed_use)) then if (.NOT.(ANY(crtm_rtsolution_associated(rtsolution0)))) & write(6,*)' ***ERROR** creating rtsolution0.' endif if (.NOT.(ANY(crtm_rtsolution_associated(rtsolution_k)))) & write(6,*)myname_,' ***ERROR** creating rtsolution_k.' if (.NOT.(ANY(crtm_options_associated(options)))) & write(6,*)myname_,' ***ERROR** creating options.' ! Turn off antenna correction options(1)%use_antenna_correction = .false. ! Load surface sensor data structure surface(1)%sensordata%n_channels = channelinfo(sensorindex)%n_channels !! REL-1.2 CRTM !! surface(1)%sensordata%select_wmo_sensor_id = channelinfo(1)%wmo_sensor_id !! RB-1.1.rev1855 CRTM atmosphere(1)%n_layers = msig atmosphere(1)%absorber_id(1) = H2O_ID atmosphere(1)%absorber_id(2) = O3_ID atmosphere(1)%absorber_units(1) = MASS_MIXING_RATIO_UNITS atmosphere(1)%absorber_units(2) = VOLUME_MIXING_RATIO_UNITS atmosphere(1)%level_pressure(0) = TOA_PRESSURE ! Currently all considered trace gases affect CRTM. Load trace gases into CRTM atmosphere ico2=-1 if (n_ghg>0) then do ig=1,n_ghg j = min_n_absorbers + ig select case(trim(ghg_names(ig))) case('co2'); atmosphere(1)%absorber_id(j) = CO2_ID case('ch4'); atmosphere(1)%absorber_id(j) = CH4_ID case('n2o'); atmosphere(1)%absorber_id(j) = N2O_ID case('co') ; atmosphere(1)%absorber_id(j) = CO_ID case default call die(myname_,': invalid absorber TERMINATE PROGRAM'//trim(ghg_names(ig)),71) end select atmosphere(1)%absorber_units(j) = VOLUME_MIXING_RATIO_UNITS if (trim(ghg_names(ig))=='co2') ico2=j enddo endif ico24crtm=-1 if (ico2>0) call gsi_chemguess_get ( 'i4crtm::co2', ico24crtm, ier ) ! Allocate structure for _k arrays (jacobians) do ii=1,nchanl atmosphere_k(ii,1) = atmosphere(1) surface_k(ii,1) = surface(1) end do if (mixed_use) then do ii=1,nchanl atmosphere_k_clr(ii,1) = atmosphere(1) surface_k_clr(ii,1) = surface(1) end do end if ! Mapping land surface type to CRTM surface fields if (regional .or. nvege_type==IGBP_N_TYPES) then allocate(map_to_crtm_ir(nvege_type)) allocate(map_to_crtm_mwave(nvege_type)) if(nvege_type==USGS_N_TYPES)then ! Assign mapping for CRTM microwave calculations map_to_crtm_mwave=usgs_to_gfs ! map usgs to CRTM select case ( TRIM(CRTM_IRlandCoeff_Classification()) ) case('NPOESS'); map_to_crtm_ir=usgs_to_npoess case('USGS') ; map_to_crtm_ir=usgs_to_usgs end select else if(nvege_type==IGBP_N_TYPES)then ! Assign mapping for CRTM microwave calculations map_to_crtm_mwave=igbp_to_gfs ! nmm igbp to CRTM select case ( TRIM(CRTM_IRlandCoeff_Classification()) ) case('NPOESS'); map_to_crtm_ir=igbp_to_npoess case('IGBP') ; map_to_crtm_ir=igbp_to_igbp end select else write(6,*)myname_,': ***ERROR*** invalid vegetation types' & //' for the CRTM IRland EmisCoeff file used.', & ' (only 20 and 24 are setup) nvege_type=',nvege_type, & ' ***STOP IN SETUPRAD***' call stop2(71) endif ! nvege_type endif ! regional or IGBP ! Calculate RH when aerosols are present and/or cloud-fraction used if (n_actual_aerosols_wk>0) then allocate(gesqsat(lat2,lon2,nsig,nfldsig)) ice=.true. iderivative=0 do ii=1,nfldsig call genqsat(gesqsat(1,1,1,ii),ges_tsen(1,1,1,ii),ges_prsl(1,1,1,ii),lat2,lon2,nsig,ice,iderivative) end do endif return end subroutine init_crtm subroutine destroy_crtm !$$$ subprogram documentation block ! . . . . ! subprogram: destroy_crtm deallocates crtm arrays ! prgmmr: parrish org: np22 date: 2005-01-22 ! ! abstract: deallocates crtm arrays ! ! program history log: ! 2010-08-17 derber ! ! input argument list: ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ implicit none character(len=*),parameter::myname_ = myname//'*destroy_crtm' integer(i_kind) error_status error_status = crtm_destroy(channelinfo) if (error_status /= success) & write(6,*)myname_,': ***ERROR*** error_status=',error_status if (n_actual_aerosols_wk>0) then deallocate(gesqsat) endif call crtm_atmosphere_destroy(atmosphere(1)) call crtm_surface_destroy(surface(1)) if (n_clouds_fwd_wk>0 .and. (.not. mixed_use)) & call crtm_rtsolution_destroy(rtsolution0) call crtm_rtsolution_destroy(rtsolution) call crtm_rtsolution_destroy(rtsolution_k) if (mixed_use) then call crtm_rtsolution_destroy(rtsolution_clr) call crtm_rtsolution_destroy(rtsolution_k_clr) end if call crtm_options_destroy(options) if (crtm_atmosphere_associated(atmosphere(1))) & write(6,*)myname_,' ***ERROR** destroying atmosphere.' if (crtm_surface_associated(surface(1))) & write(6,*)myname_,' ***ERROR** destroying surface.' if (ANY(crtm_rtsolution_associated(rtsolution))) & write(6,*)myname_,' ***ERROR** destroying rtsolution.' if (n_clouds_fwd_wk>0 .and. (.not. mixed_use)) then if (ANY(crtm_rtsolution_associated(rtsolution0))) & write(6,*)' ***ERROR** destroying rtsolution0.' endif if (ANY(crtm_rtsolution_associated(rtsolution_k))) & write(6,*)myname_,' ***ERROR** destroying rtsolution_k.' if (ANY(crtm_options_associated(options))) & write(6,*)myname_,' ***ERROR** destroying options.' deallocate(rtsolution,atmosphere_k,surface_k,rtsolution_k) if (mixed_use) deallocate(rtsolution_clr,atmosphere_k_clr, & surface_k_clr,rtsolution_k_clr) if (n_clouds_fwd_wk>0 .and. (.not. mixed_use)) & deallocate(rtsolution0) if(n_actual_aerosols_wk>0)then deallocate(aero,aero_conc,auxrh) if(n_aerosols_jac_wk>0) deallocate(iaero_jac) endif if (n_ghg>0) then deallocate(ghg_names) endif if(allocated(icloud)) deallocate(icloud) if(allocated(cloud)) deallocate(cloud) if(allocated(cloudefr)) deallocate(cloudefr) if(allocated(jcloud)) deallocate(jcloud) if(allocated(cloud_cont)) deallocate(cloud_cont) if(allocated(cloud_efr)) deallocate(cloud_efr) if(allocated(icw)) deallocate(icw) if(allocated(lcloud4crtm_wk)) deallocate(lcloud4crtm_wk) if(regional .or. nvege_type==IGBP_N_TYPES)deallocate(map_to_crtm_ir) if(regional .or. nvege_type==IGBP_N_TYPES)deallocate(map_to_crtm_mwave) return end subroutine destroy_crtm subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & h,q,clw_guess,prsl,prsi, & trop5,tzbgr,dtsavg,sfc_speed,& tsim,emissivity,ptau5,ts, & emissivity_k,temp,wmix,jacobian,error_status,tsim_clr, & layer_od,jacobian_aero) !$$$ subprogram documentation block ! . . . . ! subprogram: call_crtm creates vertical profile of t,q,oz,p,zs,etc., ! calls crtm, and does adjoint of creation (where necessary) for setuprad ! prgmmr: parrish org: np22 date: 1990-10-11 ! ! abstract: creates vertical profile of t,q,oz,p,zs,etc., ! calls crtm, and does adjoint of creation (where necessary) for setuprad ! ! program history log: ! 2010-08-17 derber - modify from intrppx and add threading ! 2011-02-23 todling/da silva - revisit interface to fill in aerosols ! 2011-03-25 yang - turn off the drop-off of co2 amount when using climatological CO2 ! 2011-05-03 todling - merge with Min-Jeong's MW cloudy radiance; combine w/ metguess ! (did not include tendencies since they were calc but not used) ! 2011-05-17 auligne/todling - add handling for hydrometeors ! 2011-06-29 todling - no explict reference to internal bundle arrays ! 2011-07-05 zhu - add cloud_efr & cloudefr; add cloud_efr & jcloud in the interface of Set_CRTM_Cloud ! 2011-07-05 zhu - rewrite cloud_cont & cwj for cloud control variables (lcw4crtm) ! 2012-03-12 veldelst-- add a internal interpolation function (option) ! 2012-04-25 yang - modify to use trace gas chem_bundle. Trace gas variables are ! invoked by the global_anavinfo.ghg.l64.txt ! 2013-02-25 zhu - add cold_start option for regional applications ! 2014-01-31 mkim-- remove 60.0degree boundary for icmask for all-sky MW radiance DA ! 2014-02-26 zhu - add non zero jacobian so jacobian will be produced for ! clear-sky background or background with small amount of cloud ! 2014-04-27 eliu - add option to calculate clear-sky Tb under all-sky condition ! 2015-02-27 eliu-- wind direction fix for using CRTM FASTEM model ! 2015-03-23 zaizhong ma - add Himawari-8 ahi ! 2015-09-10 zhu - generalize enabling all-sky and aerosol usage in radiance assimilation, ! use n_clouds_fwd_wk,n_aerosols_fwd_wk,cld_sea_only_wk, cld_sea_only_wk,cw_cv,etc ! ! input argument list: ! obstype - type of observations for which to get profile ! obstime - time of observations for which to get profile ! data_s - array containing input data information ! nchanl - number of channels ! nreal - number of descriptor information in data_s ! ich - channel number array ! ! output argument list: ! h - interpolated temperature ! q - interpolated specific humidity (max(qsmall,q)) ! prsl - interpolated layer pressure (nsig) ! prsi - interpolated level pressure (nsig+1) ! trop5 - interpolated tropopause pressure ! tzbgr - water surface temperature used in Tz retrieval ! dtsavg - delta average skin temperature over surface types ! uu5 - interpolated bottom sigma level zonal wind ! vv5 - interpolated bottom sigma level meridional wind ! tsim - simulated brightness temperatures ! emissivity - surface emissivities ! ptau5 - level transmittances ! ts - skin temperature sensitivities ! emissivity_k - surface emissivity sensitivities ! temp - temperature sensitivities ! wmix - humidity sensitivities ! jacobian - nsigradjac level jacobians for use in intrad and stprad ! error_status - error status from crtm ! layer_od - layer optical depth ! jacobian_aero- nsigaerojac level jacobians for use in intaod ! tsim_clr - option to output simulated brightness temperatures for clear sky ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ !-------- use kinds, only: r_kind,i_kind use mpimod, only: mype use radinfo, only: ifactq use radinfo, only: nsigradjac use gsi_nstcouplermod, only: nst_gsi use guess_grids, only: ges_tsen,& ges_prsl,ges_prsi,tropprs,dsfct,add_rtm_layers, & hrdifsig,nfldsig,hrdifsfc,nfldsfc,ntguessfc,isli2,sno2 use cloud_efr_mod, only: efr_ql,efr_qi,efr_qr,efr_qs,efr_qg,efr_qh use gsi_bundlemod, only: gsi_bundlegetpointer use gsi_chemguess_mod, only: gsi_chemguess_bundle ! for now, a common block use gsi_chemguess_mod, only: gsi_chemguess_get use gsi_metguess_mod, only: gsi_metguess_bundle ! for now, a common block use gsi_metguess_mod, only: gsi_metguess_get use gridmod, only: istart,jstart,nlon,nlat,lon1 use wrf_params_mod, only: cold_start use constants, only: zero,half,one,one_tenth,fv,r0_05,r10,r100,r1000,constoz,grav,rad2deg, & sqrt_tiny_r_kind,constoz,two, three, four,five,t0c use constants, only: max_varname_length,pi use set_crtm_aerosolmod, only: set_crtm_aerosol use set_crtm_cloudmod, only: set_crtm_cloud use crtm_module, only: limit_exp,o3_id use obsmod, only: iadate use aeroinfo, only: nsigaerojac implicit none ! Declare passed variables real(r_kind) ,intent(in ) :: obstime integer(i_kind) ,intent(in ) :: nchanl,nreal integer(i_kind),dimension(nchanl) ,intent(in ) :: ich real(r_kind) ,intent( out) :: trop5,tzbgr real(r_kind),dimension(nsig) ,intent( out) :: h,q,prsl real(r_kind),dimension(nsig+1) ,intent( out) :: prsi real(r_kind) ,intent( out) :: sfc_speed,dtsavg real(r_kind),dimension(nchanl+nreal) ,intent(in ) :: data_s real(r_kind),dimension(nchanl) ,intent( out) :: tsim,emissivity,ts,emissivity_k character(10) ,intent(in ) :: obstype integer(i_kind) ,intent( out) :: error_status real(r_kind),dimension(nsig,nchanl) ,intent( out) :: temp,ptau5,wmix real(r_kind),dimension(nsigradjac,nchanl),intent(out):: jacobian real(r_kind) ,intent( out) :: clw_guess real(r_kind),dimension(nchanl) ,intent( out), optional :: tsim_clr real(r_kind),dimension(nsigaerojac,nchanl),intent(out),optional :: jacobian_aero real(r_kind),dimension(nsig,nchanl) ,intent( out) ,optional :: layer_od ! Declare local parameters character(len=*),parameter::myname_=myname//'*call_crtm' real(r_kind),parameter:: minsnow=one_tenth real(r_kind),parameter:: qsmall = 1.e-6_r_kind real(r_kind),parameter:: ozsmall = 1.e-10_r_kind real(r_kind),parameter:: small_wind = 1.e-3_r_kind real(r_kind),parameter:: windscale = 999999.0_r_kind real(r_kind),parameter:: windlimit = 0.0001_r_kind real(r_kind),parameter:: quadcof (4, 2 ) = & reshape((/0.0_r_kind, 1.0_r_kind, 1.0_r_kind, 2.0_r_kind, 1.0_r_kind, & -1.0_r_kind, 1.0_r_kind, -1.0_r_kind/), (/4, 2/)) ! Declare local variables integer(i_kind):: iquadrant integer(i_kind):: ier,ii,kk,kk2,i,itype,leap_day,day_of_year integer(i_kind):: ig,istatus integer(i_kind):: j,k,m1,ix,ix1,ixp,iy,iy1,iyp,m,iii integer(i_kind):: itsig,itsigp,itsfc,itsfcp integer(i_kind):: istyp00,istyp01,istyp10,istyp11 integer(i_kind):: iqs,iozs integer(i_kind):: error_status_clr integer(i_kind),dimension(8)::obs_time,anal_time integer(i_kind),dimension(msig) :: klevel ! ****************************** ! Constrained indexing for lai ! CRTM 2.1 implementation change ! ****************************** integer(i_kind):: lai_type real(r_kind):: wind10,wind10_direction,windratio,windangle real(r_kind):: w00,w01,w10,w11,kgkg_kgm2,f10,panglr,dx,dy ! real(r_kind):: w_weights(4) real(r_kind):: delx,dely,delx1,dely1,dtsig,dtsigp,dtsfc,dtsfcp real(r_kind):: sst00,sst01,sst10,sst11,total_od,term,uu5,vv5, ps real(r_kind):: sno00,sno01,sno10,sno11,secant_term real(r_kind),dimension(0:3):: wgtavg real(r_kind),dimension(nsig,nchanl):: omix real(r_kind),dimension(nsig,nchanl,n_aerosols_jac):: jaero real(r_kind),dimension(nchanl) :: uwind_k,vwind_k real(r_kind),dimension(msig+1) :: prsi_rtm real(r_kind),dimension(msig) :: prsl_rtm real(r_kind),dimension(msig) :: auxq,auxdp real(r_kind),dimension(nsig) :: poz real(r_kind),dimension(nsig) :: rh,qs real(r_kind),dimension(5) :: tmp_time real(r_kind),dimension(0:3) :: dtskin real(r_kind),dimension(msig) :: c6 real(r_kind),dimension(nsig) :: c2,c3,c4,c5 real(r_kind),dimension(nsig) :: ugkg_kgm2,cwj real(r_kind),allocatable,dimension(:,:) :: tgas1d real(r_kind),pointer,dimension(:,: )::psges_itsig =>NULL() real(r_kind),pointer,dimension(:,: )::psges_itsigp=>NULL() real(r_kind),pointer,dimension(:,:,:)::uges_itsig =>NULL() real(r_kind),pointer,dimension(:,:,:)::uges_itsigp=>NULL() real(r_kind),pointer,dimension(:,:,:)::vges_itsig =>NULL() real(r_kind),pointer,dimension(:,:,:)::vges_itsigp=>NULL() real(r_kind),pointer,dimension(:,:,:)::qges_itsig =>NULL() real(r_kind),pointer,dimension(:,:,:)::qges_itsigp=>NULL() real(r_kind),pointer,dimension(:,:,:)::ozges_itsig =>NULL() real(r_kind),pointer,dimension(:,:,:)::ozges_itsigp=>NULL() real(r_kind),pointer,dimension(:,:,:)::tgasges_itsig =>NULL() real(r_kind),pointer,dimension(:,:,:)::tgasges_itsigp=>NULL() real(r_kind),pointer,dimension(:,:,:)::aeroges_itsig =>NULL() real(r_kind),pointer,dimension(:,:,:)::aeroges_itsigp=>NULL() logical :: sea,icmask integer(i_kind),parameter,dimension(12):: mday=(/0,31,59,90,& 120,151,181,212,243,273,304,334/) real(r_kind) :: lai m1=mype+1 dx = data_s(ilat) ! grid relative latitude dy = data_s(ilon) ! grid relative longitude ! Set spatial interpolation indices and weights ix1=dx ix1=max(1,min(ix1,nlat)) delx=dx-ix1 delx=max(zero,min(delx,one)) ix=ix1-istart(m1)+2 ixp=ix+1 if(ix1==nlat) then ixp=ix end if delx1=one-delx iy1=dy dely=dy-iy1 iy=iy1-jstart(m1)+2 if(iy<1) then iy1=iy1+nlon iy=iy1-jstart(m1)+2 end if if(iy>lon1+1) then iy1=iy1-nlon iy=iy1-jstart(m1)+2 end if iyp=iy+1 dely1=one-dely w00=delx1*dely1; w10=delx*dely1; w01=delx1*dely; w11=delx*dely ! w_weights = (/w00,w10,w01,w11/) ! Get time interpolation factors for sigma files if(obstime > hrdifsig(1) .and. obstime < hrdifsig(nfldsig))then do j=1,nfldsig-1 if(obstime > hrdifsig(j) .and. obstime <= hrdifsig(j+1))then itsig=j itsigp=j+1 dtsig=((hrdifsig(j+1)-obstime)/(hrdifsig(j+1)-hrdifsig(j))) end if end do else if(obstime <=hrdifsig(1))then itsig=1 itsigp=1 dtsig=one else itsig=nfldsig itsigp=nfldsig dtsig=one end if dtsigp=one-dtsig ! Get time interpolation factors for surface files if(obstime > hrdifsfc(1) .and. obstime < hrdifsfc(nfldsfc))then do j=1,nfldsfc-1 if(obstime > hrdifsfc(j) .and. obstime <= hrdifsfc(j+1))then itsfc=j itsfcp=j+1 dtsfc=((hrdifsfc(j+1)-obstime)/(hrdifsfc(j+1)-hrdifsfc(j))) end if end do else if(obstime <=hrdifsfc(1))then itsfc=1 itsfcp=1 dtsfc=one else itsfc=nfldsfc itsfcp=nfldsfc dtsfc=one end if dtsfcp=one-dtsfc ier=0 call gsi_bundlegetpointer(gsi_metguess_bundle(itsig ),'ps',psges_itsig ,istatus) ier=ier+istatus call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'ps',psges_itsigp,istatus) ier=ier+istatus call gsi_bundlegetpointer(gsi_metguess_bundle(itsig ),'u' ,uges_itsig ,istatus) ier=ier+istatus call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'u' ,uges_itsigp ,istatus) ier=ier+istatus call gsi_bundlegetpointer(gsi_metguess_bundle(itsig ),'v' ,vges_itsig ,istatus) ier=ier+istatus call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'v' ,vges_itsigp ,istatus) ier=ier+istatus call gsi_bundlegetpointer(gsi_metguess_bundle(itsig ),'oz',ozges_itsig ,iozs) iozs=istatus call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'oz',ozges_itsigp,iozs) iozs=iozs+istatus call gsi_bundlegetpointer(gsi_metguess_bundle(itsig ),'q',qges_itsig ,istatus) iqs=istatus call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'q',qges_itsigp,istatus) iqs=iqs+istatus ! Space-time interpolation of temperature (h) and q fields from sigma files !$omp parallel do schedule(dynamic,1) private(k,ii,iii) do k=1,nsig if(k == 1)then jacobian=zero ! Set surface type flag. (Same logic as in subroutine deter_sfc) istyp00 = isli2(ix ,iy ) istyp10 = isli2(ixp,iy ) istyp01 = isli2(ix ,iyp) istyp11 = isli2(ixp,iyp) sno00= sno2(ix ,iy ,itsfc)*dtsfc+sno2(ix ,iy ,itsfcp)*dtsfcp sno01= sno2(ix ,iyp,itsfc)*dtsfc+sno2(ix ,iyp,itsfcp)*dtsfcp sno10= sno2(ixp,iy ,itsfc)*dtsfc+sno2(ixp,iy ,itsfcp)*dtsfcp sno11= sno2(ixp,iyp,itsfc)*dtsfc+sno2(ixp,iyp,itsfcp)*dtsfcp if(istyp00 >= 1 .and. sno00 > minsnow)istyp00 = 3 if(istyp01 >= 1 .and. sno01 > minsnow)istyp01 = 3 if(istyp10 >= 1 .and. sno10 > minsnow)istyp10 = 3 if(istyp11 >= 1 .and. sno11 > minsnow)istyp11 = 3 ! Find delta Surface temperatures for all surface types sst00= dsfct(ix ,iy,ntguessfc) ; sst01= dsfct(ix ,iyp,ntguessfc) sst10= dsfct(ixp,iy,ntguessfc) ; sst11= dsfct(ixp,iyp,ntguessfc) dtsavg=sst00*w00+sst10*w10+sst01*w01+sst11*w11 dtskin(0:3)=zero wgtavg(0:3)=zero if(istyp00 == 1)then wgtavg(1) = wgtavg(1) + w00 dtskin(1)=dtskin(1)+w00*sst00 else if(istyp00 == 2)then wgtavg(2) = wgtavg(2) + w00 dtskin(2)=dtskin(2)+w00*sst00 else if(istyp00 == 3)then wgtavg(3) = wgtavg(3) + w00 dtskin(3)=dtskin(3)+w00*sst00 else wgtavg(0) = wgtavg(0) + w00 dtskin(0)=dtskin(0)+w00*sst00 end if if(istyp01 == 1)then wgtavg(1) = wgtavg(1) + w01 dtskin(1)=dtskin(1)+w01*sst01 else if(istyp01 == 2)then wgtavg(2) = wgtavg(2) + w01 dtskin(2)=dtskin(2)+w01*sst01 else if(istyp01 == 3)then wgtavg(3) = wgtavg(3) + w01 dtskin(3)=dtskin(3)+w01*sst01 else wgtavg(0) = wgtavg(0) + w01 dtskin(0)=dtskin(0)+w01*sst01 end if if(istyp10 == 1)then wgtavg(1) = wgtavg(1) + w10 dtskin(1)=dtskin(1)+w10*sst10 else if(istyp10 == 2)then wgtavg(2) = wgtavg(2) + w10 dtskin(2)=dtskin(2)+w10*sst10 else if(istyp10 == 3)then wgtavg(3) = wgtavg(3) + w10 dtskin(3)=dtskin(3)+w10*sst10 else wgtavg(0) = wgtavg(0) + w10 dtskin(0)=dtskin(0)+w10*sst10 end if if(istyp11 == 1)then wgtavg(1) = wgtavg(1) + w11 dtskin(1)=dtskin(1)+w11*sst11 else if(istyp11 == 2)then wgtavg(2) = wgtavg(2) + w11 dtskin(2)=dtskin(2)+w11*sst11 else if(istyp11 == 3)then wgtavg(3) = wgtavg(3) + w11 dtskin(3)=dtskin(3)+w11*sst11 else wgtavg(0) = wgtavg(0) + w11 dtskin(0)=dtskin(0)+w11*sst11 end if if(wgtavg(0) > zero)then dtskin(0) = dtskin(0)/wgtavg(0) else dtskin(0) = dtsavg end if if(wgtavg(1) > zero)then dtskin(1) = dtskin(1)/wgtavg(1) else dtskin(1) = dtsavg end if if(wgtavg(2) > zero)then dtskin(2) = dtskin(2)/wgtavg(2) else dtskin(2) = dtsavg end if if(wgtavg(3) > zero)then dtskin(3) = dtskin(3)/wgtavg(3) else dtskin(3) = dtsavg end if if (n_clouds_fwd_wk>0) then ps=(psges_itsig (ix,iy )*w00+psges_itsig (ixp,iy )*w10+ & psges_itsig (ix,iyp)*w01+psges_itsig (ixp,iyp)*w11)*dtsig + & (psges_itsigp(ix,iy )*w00+psges_itsigp(ixp,iy )*w10+ & psges_itsigp(ix,iyp)*w01+psges_itsigp(ixp,iyp)*w11)*dtsigp endif ! skip loading surface structure if obstype is modis_aod if (trim(obstype) /= 'modis_aod') then ! Load surface structure ! **NOTE: The model surface type --> CRTM surface type ! mapping below is specific to the versions NCEP ! GFS and NNM as of Summer 2016 itype = nint(data_s(ivty)) istype = nint(data_s(isty)) if (regional .or. nvege_type==IGBP_N_TYPES) then itype = min(max(1,itype),nvege_type) istype = min(max(1,istype),SOIL_N_TYPES) if (ChannelInfo(sensorindex)%sensor_type == crtm_microwave_sensor)then surface(1)%land_type = max(1,map_to_crtm_mwave(itype)) else surface(1)%land_type = max(1,map_to_crtm_ir(itype)) end if surface(1)%Vegetation_Type = max(1,map_to_crtm_mwave(itype)) surface(1)%Soil_Type = map_soil_to_crtm(istype) lai_type = map_to_crtm_mwave(itype) elseif (nvege_type==GFS_N_TYPES) then itype = min(max(0,itype),GFS_VEGETATION_N_TYPES) istype = min(max(1,istype),GFS_SOIL_N_TYPES) surface(1)%land_type = gfs_to_crtm(itype) surface(1)%Vegetation_Type = max(1,itype) surface(1)%Soil_Type = istype lai_type = itype else write(6,*)myname_,': ***ERROR*** invalid vegetation types' & //' the information does not match any currenctly.', & ' supported surface type maps to the CRTM,', & ' ***STOP IN SETUPRAD***' call stop2(71) end if if (lwind) then ! Interpolate lowest level winds to observation location uu5=(uges_itsig (ix,iy ,1)*w00+uges_itsig (ixp,iy ,1)*w10+ & uges_itsig (ix,iyp,1)*w01+uges_itsig (ixp,iyp,1)*w11)*dtsig + & (uges_itsigp(ix,iy ,1)*w00+uges_itsigp(ixp,iy ,1)*w10+ & uges_itsigp(ix,iyp,1)*w01+uges_itsigp(ixp,iyp,1)*w11)*dtsigp vv5=(vges_itsig (ix,iy ,1)*w00+vges_itsig (ixp,iy ,1)*w10+ & vges_itsig (ix,iyp,1)*w01+vges_itsig (ixp,iyp,1)*w11)*dtsig + & (vges_itsigp(ix,iy ,1)*w00+vges_itsigp(ixp,iy ,1)*w10+ & vges_itsigp(ix,iyp,1)*w01+vges_itsigp(ixp,iyp,1)*w11)*dtsigp f10=data_s(iff10) sfc_speed = f10*sqrt(uu5*uu5+vv5*vv5) wind10 = sfc_speed if (uu5*f10 >= 0.0_r_kind .and. vv5*f10 >= 0.0_r_kind) iquadrant = 1 if (uu5*f10 >= 0.0_r_kind .and. vv5*f10 < 0.0_r_kind) iquadrant = 2 if (uu5*f10 < 0.0_r_kind .and. vv5*f10 >= 0.0_r_kind) iquadrant = 4 if (uu5*f10 < 0.0_r_kind .and. vv5*f10 < 0.0_r_kind) iquadrant = 3 if (abs(vv5*f10) >= windlimit) then windratio = (uu5*f10) / (vv5*f10) else windratio = 0.0_r_kind if (abs(uu5*f10) > windlimit) then windratio = windscale * uu5*f10 endif endif windangle = atan(abs(windratio)) ! wind azimuth is in radians wind10_direction = quadcof(iquadrant, 1) * pi + windangle * quadcof(iquadrant, 2) surface(1)%wind_speed = sfc_speed surface(1)%wind_direction = rad2deg*wind10_direction else !RTodling: not sure the following option makes any sense surface(1)%wind_speed = zero surface(1)%wind_direction = zero endif ! CRTM will reject surface coverages if greater than one and it is possible for ! these values to be larger due to round off. surface(1)%water_coverage = min(max(zero,data_s(ifrac_sea)),one) surface(1)%land_coverage = min(max(zero,data_s(ifrac_lnd)),one) surface(1)%ice_coverage = min(max(zero,data_s(ifrac_ice)),one) surface(1)%snow_coverage = min(max(zero,data_s(ifrac_sno)),one) ! ! get vegetation lai from summer and winter values. ! surface(1)%Lai = zero if (surface(1)%land_coverage>zero) then if(lai_type>0)then call get_lai(data_s,nchanl,nreal,itime,ilate,lai_type,lai) surface(1)%Lai = lai ! LAI endif ! for Glacial land ice soil type and vegetation type if(surface(1)%Soil_Type == 9 .OR. surface(1)%Vegetation_Type == 13) then surface(1)%ice_coverage = min(surface(1)%ice_coverage + surface(1)%land_coverage, one) surface(1)%land_coverage = zero endif endif surface(1)%water_temperature = max(data_s(its_sea)+dtskin(0),270._r_kind) if(nst_gsi>1 .and. surface(1)%water_coverage>zero) then surface(1)%water_temperature = max(data_s(itref)+data_s(idtw)-data_s(idtc)+dtskin(0),271._r_kind) endif surface(1)%land_temperature = data_s(its_lnd)+dtskin(1) surface(1)%ice_temperature = min(data_s(its_ice)+dtskin(2),280._r_kind) surface(1)%snow_temperature = min(data_s(its_sno)+dtskin(3),280._r_kind) surface(1)%soil_moisture_content = data_s(ism) surface(1)%vegetation_fraction = data_s(ivfr) surface(1)%soil_temperature = data_s(istp) surface(1)%snow_depth = data_s(isn) sea = min(max(zero,data_s(ifrac_sea)),one) >= 0.99_r_kind icmask = (sea .and. cld_sea_only_wk) .or. (.not. cld_sea_only_wk) ! assign tzbgr for Tz retrieval when necessary tzbgr = surface(1)%water_temperature endif ! end of loading surface structure ! Load geometry structure ! skip loading geometry structure if obstype is modis_aod ! iscan_ang,ilzen_ang,ilazi_ang are not available in the modis aod bufr file ! also, geometryinfo is not needed in crtm aod calculation if ( trim(obstype) /= 'modis_aod' ) then panglr = data_s(iscan_ang) if(obstype == 'goes_img' .or. obstype == 'seviri' .or. obstype == 'abi')panglr = zero geometryinfo(1)%sensor_zenith_angle = abs(data_s(ilzen_ang)*rad2deg) ! local zenith angle geometryinfo(1)%source_zenith_angle = abs(data_s(iszen_ang)) ! solar zenith angle geometryinfo(1)%sensor_azimuth_angle = data_s(ilazi_ang) ! local azimuth angle geometryinfo(1)%source_azimuth_angle = data_s(isazi_ang) ! solar azimuth angle geometryinfo(1)%sensor_scan_angle = panglr*rad2deg ! scan angle geometryinfo(1)%ifov = nint(data_s(iscan_pos)) ! field of view position ! For some microwave instruments the solar and sensor azimuth angles can be ! missing (given a value of 10^11). Set these to zero to get past CRTM QC. if (geometryinfo(1)%source_azimuth_angle > 360.0_r_kind .OR. & geometryinfo(1)%source_azimuth_angle < zero ) & geometryinfo(1)%source_azimuth_angle = zero if (geometryinfo(1)%sensor_azimuth_angle > 360.0_r_kind .OR. & geometryinfo(1)%sensor_azimuth_angle < zero ) & geometryinfo(1)%sensor_azimuth_angle = zero endif ! end of loading geometry structure ! Special block for SSU cell pressure leakage correction. Need to compute ! observation time and load into Time component of geometryinfo structure. ! geometryinfo%time is only defined in CFSRR CRTM. if (obstype == 'ssu') then ! Compute absolute observation time anal_time=0 obs_time=0 tmp_time=zero tmp_time(2)=obstime anal_time(1)=iadate(1) anal_time(2)=iadate(2) anal_time(3)=iadate(3) anal_time(5)=iadate(4) !external-subroutine w3movdat() call w3movdat(tmp_time,anal_time,obs_time) ! Compute decimal year, for example 1/10/1983 ! d_year = 1983.0 + 10.0/365.0 leap_day = 0 if( mod(obs_time(1),4)==0 ) then if( (mod(obs_time(1),100)/=0).or.(mod(obs_time(1),400)==0) ) leap_day = 1 endif day_of_year = mday(obs_time(2)) + obs_time(3) if(obs_time(2) > 2) day_of_year = day_of_year + leap_day call ssu_input_setvalue( options%SSU, & Time=dble(obs_time(1)) + dble(day_of_year)/(365.0_r_kind+leap_day)) endif ! Load surface sensor data structure do i=1,nchanl ! Set-up to return Tb jacobians. rtsolution_k(i,1)%radiance = zero rtsolution_k(i,1)%brightness_temperature = one if (mixed_use) then rtsolution_k_clr(i,1)%radiance = zero rtsolution_k_clr(i,1)%brightness_temperature = one end if if (trim(obstype) /= 'modis_aod')then ! Pass CRTM array of tb for surface emissiviy calculations if ( channelinfo(1)%sensor_type == crtm_microwave_sensor .and. & crtm_surface_associated(surface(1)) ) & surface(1)%sensordata%tb(i) = data_s(nreal+i) ! set up to return layer_optical_depth jacobians rtsolution_k(i,1)%layer_optical_depth = one if (mixed_use) rtsolution_k_clr(i,1)%layer_optical_depth = one endif end do end if h(k) =(ges_tsen(ix ,iy ,k,itsig )*w00+ & ges_tsen(ixp,iy ,k,itsig )*w10+ & ges_tsen(ix ,iyp,k,itsig )*w01+ & ges_tsen(ixp,iyp,k,itsig )*w11)*dtsig + & (ges_tsen(ix ,iy ,k,itsigp)*w00+ & ges_tsen(ixp,iy ,k,itsigp)*w10+ & ges_tsen(ix ,iyp,k,itsigp)*w01+ & ges_tsen(ixp,iyp,k,itsigp)*w11)*dtsigp ! Interpolate layer pressure to observation point prsl(k)=(ges_prsl(ix ,iy ,k,itsig )*w00+ & ges_prsl(ixp,iy ,k,itsig )*w10+ & ges_prsl(ix ,iyp,k,itsig )*w01+ & ges_prsl(ixp,iyp,k,itsig )*w11)*dtsig + & (ges_prsl(ix ,iy ,k,itsigp)*w00+ & ges_prsl(ixp,iy ,k,itsigp)*w10+ & ges_prsl(ix ,iyp,k,itsigp)*w01+ & ges_prsl(ixp,iyp,k,itsigp)*w11)*dtsigp ! Interpolate level pressure to observation point prsi(k)=(ges_prsi(ix ,iy ,k,itsig )*w00+ & ges_prsi(ixp,iy ,k,itsig )*w10+ & ges_prsi(ix ,iyp,k,itsig )*w01+ & ges_prsi(ixp,iyp,k,itsig )*w11)*dtsig + & (ges_prsi(ix ,iy ,k,itsigp)*w00+ & ges_prsi(ixp,iy ,k,itsigp)*w10+ & ges_prsi(ix ,iyp,k,itsigp)*w01+ & ges_prsi(ixp,iyp,k,itsigp)*w11)*dtsigp if (iqs==0) then q(k) =(qges_itsig (ix ,iy ,k)*w00+ & qges_itsig (ixp,iy ,k)*w10+ & qges_itsig (ix ,iyp,k)*w01+ & qges_itsig (ixp,iyp,k)*w11)*dtsig + & (qges_itsigp(ix ,iy ,k)*w00+ & qges_itsigp(ixp,iy ,k)*w10+ & qges_itsigp(ix ,iyp,k)*w01+ & qges_itsigp(ixp,iyp,k)*w11)*dtsigp ! Ensure q is greater than or equal to qsmall q(k)=max(qsmall,q(k)) else q(k) = qsmall endif c2(k)=one/(one+fv*q(k)) c3(k)=one/(one-q(k)) c4(k)=fv*h(k)*c2(k) c5(k)=r1000*c3(k)*c3(k) ! Space-time interpolation of ozone(poz) if (iozs==0) then poz(k)=((ozges_itsig (ix ,iy ,k)*w00+ & ozges_itsig (ixp,iy ,k)*w10+ & ozges_itsig (ix ,iyp,k)*w01+ & ozges_itsig (ixp,iyp,k)*w11)*dtsig + & (ozges_itsigp(ix ,iy ,k)*w00+ & ozges_itsigp(ixp,iy ,k)*w10+ & ozges_itsigp(ix ,iyp,k)*w01+ & ozges_itsigp(ixp,iyp,k)*w11)*dtsigp)*constoz ! Ensure ozone is greater than ozsmall poz(k)=max(ozsmall,poz(k)) endif ! oz ! Quantities required for MW cloudy radiance calculations if (n_clouds_fwd_wk>0) then do ii=1,n_clouds_fwd_wk iii=jcloud(ii) cloud(k,ii) =(gsi_metguess_bundle(itsig )%r3(icloud(iii))%q(ix ,iy ,k)*w00+ & ! kg/kg gsi_metguess_bundle(itsig )%r3(icloud(iii))%q(ixp,iy ,k)*w10+ & gsi_metguess_bundle(itsig )%r3(icloud(iii))%q(ix ,iyp,k)*w01+ & gsi_metguess_bundle(itsig )%r3(icloud(iii))%q(ixp,iyp,k)*w11)*dtsig + & (gsi_metguess_bundle(itsigp)%r3(icloud(iii))%q(ix ,iy ,k)*w00+ & gsi_metguess_bundle(itsigp)%r3(icloud(iii))%q(ixp,iy ,k)*w10+ & gsi_metguess_bundle(itsigp)%r3(icloud(iii))%q(ix ,iyp,k)*w01+ & gsi_metguess_bundle(itsigp)%r3(icloud(iii))%q(ixp,iyp,k)*w11)*dtsigp cloud(k,ii)=max(cloud(k,ii),zero) if (regional .and. (.not. wrf_mass_regional)) then if (trim(cloud_names(iii))== 'ql' ) then cloudefr(k,ii)=(efr_ql(ix ,iy ,k,itsig)*w00+efr_ql(ixp,iy ,k,itsig)*w10+ & efr_ql(ix ,iyp,k,itsig)*w01+efr_ql(ixp,iyp,k,itsig)*w11)*dtsig + & (efr_ql(ix ,iy ,k,itsigp)*w00+efr_ql(ixp,iy ,k,itsigp)*w10+ & efr_ql(ix ,iyp,k,itsigp)*w01+efr_ql(ixp,iyp,k,itsigp)*w11)*dtsigp else if (trim(cloud_names(iii))== 'qi' ) then cloudefr(k,ii)=(efr_qi(ix ,iy ,k,itsig)*w00+efr_qi(ixp,iy ,k,itsig)*w10+ & efr_qi(ix ,iyp,k,itsig)*w01+efr_qi(ixp,iyp,k,itsig)*w11)*dtsig + & (efr_qi(ix ,iy ,k,itsigp)*w00+efr_qi(ixp,iy ,k,itsigp)*w10+ & efr_qi(ix ,iyp,k,itsigp)*w01+efr_qi(ixp,iyp,k,itsigp)*w11)*dtsigp else if (trim(cloud_names(iii))== 'qs' ) then cloudefr(k,ii)=(efr_qs(ix ,iy ,k,itsig)*w00+efr_qs(ixp,iy ,k,itsig)*w10+ & efr_qs(ix ,iyp,k,itsig)*w01+efr_qs(ixp,iyp,k,itsig)*w11)*dtsig + & (efr_qs(ix ,iy ,k,itsigp)*w00+efr_qs(ixp,iy ,k,itsigp)*w10+ & efr_qs(ix ,iyp,k,itsigp)*w01+efr_qs(ixp,iyp,k,itsigp)*w11)*dtsigp else if (trim(cloud_names(iii))== 'qg' ) then cloudefr(k,ii)=(efr_qg(ix ,iy ,k,itsig)*w00+efr_qg(ixp,iy ,k,itsig)*w10+ & efr_qg(ix ,iyp,k,itsig)*w01+efr_qg(ixp,iyp,k,itsig)*w11)*dtsig + & (efr_qg(ix ,iy ,k,itsigp)*w00+efr_qg(ixp,iy ,k,itsigp)*w10+ & efr_qg(ix ,iyp,k,itsigp)*w01+efr_qg(ixp,iyp,k,itsigp)*w11)*dtsigp else if (trim(cloud_names(iii))== 'qh' ) then cloudefr(k,ii)=(efr_qh(ix ,iy ,k,itsig)*w00+efr_qh(ixp,iy ,k,itsig)*w10+ & efr_qh(ix ,iyp,k,itsig)*w01+efr_qh(ixp,iyp,k,itsig)*w11)*dtsig + & (efr_qh(ix ,iy ,k,itsigp)*w00+efr_qh(ixp,iy ,k,itsigp)*w10+ & efr_qh(ix ,iyp,k,itsigp)*w01+efr_qh(ixp,iyp,k,itsigp)*w11)*dtsigp else if (trim(cloud_names(iii))== 'qr' ) then cloudefr(k,ii)=(efr_qr(ix ,iy ,k,itsig)*w00+efr_qr(ixp,iy ,k,itsig)*w10+ & efr_qr(ix ,iyp,k,itsig)*w01+efr_qr(ixp,iyp,k,itsig)*w11)*dtsig + & (efr_qr(ix ,iy ,k,itsigp)*w00+efr_qr(ixp,iy ,k,itsigp)*w10+ & efr_qr(ix ,iyp,k,itsigp)*w01+efr_qr(ixp,iyp,k,itsigp)*w11)*dtsigp end if end if end do endif ! end do ! Interpolate level pressure to observation point for top interface prsi(nsig+1)=(ges_prsi(ix ,iy ,nsig+1,itsig )*w00+ & ges_prsi(ixp,iy ,nsig+1,itsig )*w10+ & ges_prsi(ix ,iyp,nsig+1,itsig )*w01+ & ges_prsi(ixp,iyp,nsig+1,itsig )*w11)*dtsig + & (ges_prsi(ix ,iy ,nsig+1,itsigp)*w00+ & ges_prsi(ixp,iy ,nsig+1,itsigp)*w10+ & ges_prsi(ix ,iyp,nsig+1,itsigp)*w01+ & ges_prsi(ixp,iyp,nsig+1,itsigp)*w11)*dtsigp ! if(any(prsl0) then allocate (tgas1d(nsig,n_ghg)) do ig=1,n_ghg if(size(gsi_chemguess_bundle)==1) then call gsi_bundlegetpointer(gsi_chemguess_bundle(1), ghg_names(ig),tgasges_itsig ,ier) do k=1,nsig ! choice: use the internal interpolation function ! or just explicitly code, not sure which one is efficient ! tgas1d(k,ig) = crtm_interface_interp(tgasges_itsig(ix:ixp,iy:iyp,:),& ! w_weights, & ! 1.0_r_kind) tgas1d(k,ig) =(tgasges_itsig(ix ,iy ,k)*w00+ & tgasges_itsig(ixp,iy ,k)*w10+ & tgasges_itsig(ix ,iyp,k)*w01+ & tgasges_itsig(ixp,iyp,k)*w11) enddo else call gsi_bundlegetpointer(gsi_chemguess_bundle(itsig ),ghg_names(ig),tgasges_itsig ,ier) call gsi_bundlegetpointer(gsi_chemguess_bundle(itsigp),ghg_names(ig),tgasges_itsigp,ier) do k=1,nsig ! tgas1d(k,ig) = crtm_interface_interp(tgasges_itsig(ix:ixp,iy:iyp,k),& ! w_weights, & ! dtsig) + & ! crtm_interface_interp(tgasges_itsigp(ix:ixp,iy:iyp,k),& ! w_weights, & ! dtsigp) tgas1d(k,ig) =(tgasges_itsig (ix ,iy ,k)*w00+ & tgasges_itsig (ixp,iy ,k)*w10+ & tgasges_itsig (ix ,iyp,k)*w01+ & tgasges_itsig (ixp,iyp,k)*w11)*dtsig + & (tgasges_itsigp(ix ,iy ,k)*w00+ & tgasges_itsigp(ixp,iy ,k)*w10+ & tgasges_itsigp(ix ,iyp,k)*w01+ & tgasges_itsigp(ixp,iyp,k)*w11)*dtsigp enddo endif enddo endif ! Space-time interpolation of aerosol fields from sigma files if(n_actual_aerosols_wk>0)then if(size(gsi_chemguess_bundle)==1) then do ii=1,n_actual_aerosols_wk call gsi_bundlegetpointer(gsi_chemguess_bundle(1),aerosol_names(ii),aeroges_itsig ,ier) do k=1,nsig aero(k,ii) =(aeroges_itsig(ix ,iy ,k)*w00+ & aeroges_itsig(ixp,iy ,k)*w10+ & aeroges_itsig(ix ,iyp,k)*w01+ & aeroges_itsig(ixp,iyp,k)*w11) end do enddo else do ii=1,n_actual_aerosols_wk call gsi_bundlegetpointer(gsi_chemguess_bundle(itsig ),aerosol_names(ii),aeroges_itsig ,ier) call gsi_bundlegetpointer(gsi_chemguess_bundle(itsigp),aerosol_names(ii),aeroges_itsigp,ier) do k=1,nsig aero(k,ii) =(aeroges_itsig (ix ,iy ,k)*w00+ & aeroges_itsig (ixp,iy ,k)*w10+ & aeroges_itsig (ix ,iyp,k)*w01+ & aeroges_itsig (ixp,iyp,k)*w11)*dtsig + & (aeroges_itsigp(ix ,iy ,k)*w00+ & aeroges_itsigp(ixp,iy ,k)*w10+ & aeroges_itsigp(ix ,iyp,k)*w01+ & aeroges_itsigp(ixp,iyp,k)*w11)*dtsigp end do enddo endif do k=1,nsig rh(k) = q(k)/qs(k) end do endif ! Find tropopause height at observation trop5= one_tenth*(tropprs(ix,iy )*w00+tropprs(ixp,iy )*w10+ & tropprs(ix,iyp)*w01+tropprs(ixp,iyp)*w11) ! Zero atmosphere jacobian structures call crtm_atmosphere_zero(atmosphere_k(:,:)) call crtm_surface_zero(surface_k(:,:)) if (mixed_use) then call crtm_atmosphere_zero(atmosphere_k_clr(:,:)) call crtm_surface_zero(surface_k_clr(:,:)) end if clw_guess = zero if (n_actual_aerosols_wk>0) then do k = 1, nsig ! Convert mixing-ratio to concentration ugkg_kgm2(k)=1.0e-9_r_kind*(prsi(k)-prsi(k+1))*r1000/grav aero(k,:)=aero(k,:)*ugkg_kgm2(k) enddo endif sea = min(max(zero,data_s(ifrac_sea)),one) >= 0.99_r_kind icmask = (sea .and. cld_sea_only_wk) .or. (.not. cld_sea_only_wk) do k = 1,msig ! Load profiles into extended RTM model layers kk = msig - k + 1 atmosphere(1)%level_pressure(k) = r10*prsi_rtm(kk) atmosphere(1)%pressure(k) = r10*prsl_rtm(kk) kk2 = klevel(kk) atmosphere(1)%temperature(k) = h(kk2) atmosphere(1)%absorber(k,1) = r1000*q(kk2)*c3(kk2) if(iozs==0) then atmosphere(1)%absorber(k,2) = poz(kk2) else atmosphere(1)%absorber(k,2) = O3_ID endif if (n_ghg > 0) then do ig=1,n_ghg j=min_n_absorbers+ ig atmosphere(1)%absorber(k,j) = tgas1d(kk2,ig) enddo endif if (n_actual_aerosols_wk>0) then aero_conc(k,:)=aero(kk2,:) auxrh(k) =rh(kk2) endif ! Include cloud guess profiles in mw radiance computation 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) then if (icmask) then c6(k) = kgkg_kgm2 auxdp(k)=abs(prsi_rtm(kk+1)-prsi_rtm(kk))*r10 auxq (k)=q(kk2) if (regional .and. (.not. wrf_mass_regional) .and. (.not. cold_start)) then do ii=1,n_clouds_fwd_wk cloud_cont(k,ii)=cloud(kk2,ii)*c6(k) cloud_efr (k,ii)=cloudefr(kk2,ii) end do else do ii=1,n_clouds_fwd_wk cloud_cont(k,ii)=cloud(kk2,ii)*c6(k) end do end if clw_guess = clw_guess + cloud_cont(k,1) do ii=1,n_clouds_fwd_wk if (ii==1 .and. atmosphere(1)%temperature(k)-t0c>-20.0_r_kind) & cloud_cont(k,1)=max(1.001_r_kind*1.0E-6_r_kind, cloud_cont(k,1)) if (ii==2 .and. atmosphere(1)%temperature(k)-20.0_r_kind) & 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)0) then atmosphere(1)%n_clouds = n_clouds_fwd_wk call Set_CRTM_Cloud (msig,n_actual_clouds_wk,cloud_names,icmask,n_clouds_fwd_wk,cloud_cont,cloud_efr,jcloud,auxdp, & atmosphere(1)%temperature,atmosphere(1)%pressure,auxq,atmosphere(1)%cloud) endif ! Set aerosols for CRTM if(n_actual_aerosols_wk>0) then call Set_CRTM_Aerosol ( msig, n_actual_aerosols_wk, n_aerosols_fwd_wk, aerosol_names, aero_conc, auxrh, & atmosphere(1)%aerosol ) endif ! Call CRTM K Matrix model error_status = 0 if ( trim(obstype) /= 'modis_aod' ) then error_status = crtm_k_matrix(atmosphere,surface,rtsolution_k,& geometryinfo,channelinfo(sensorindex:sensorindex),atmosphere_k,& surface_k,rtsolution,options=options) if (mixed_use) then ! Zero out data array in cloud structure atmosphere(1)%n_clouds = 0 error_status_clr = crtm_k_matrix(atmosphere,surface,rtsolution_k_clr,& geometryinfo,channelinfo(sensorindex:sensorindex),atmosphere_k_clr,& surface_k_clr,rtsolution_clr,options=options) end if else error_status = crtm_aod_k(atmosphere,rtsolution_k,& channelinfo(sensorindex:sensorindex),rtsolution,atmosphere_k) end if ! If the CRTM returns an error flag, do not assimilate any channels for this ob ! and set the QC flag to 10 (done in setuprad). if (error_status /=0) then write(6,*)myname_,': ***ERROR*** during crtm_k_matrix call ',& error_status end if ! Calculate clear-sky Tb for AMSU-A over sea when allsky condition is on if (n_clouds_fwd_wk>0 .and. present(tsim_clr) .and. (.not. mixed_use)) then ! Zero out data array in cloud structure: water content, effective ! radius and variance atmosphere(1)%n_clouds = 0 ! call crtm_cloud_zero(atmosphere(1)%cloud) ! call crtm forward model for clear-sky calculation error_status = crtm_forward(atmosphere,surface,& geometryinfo,channelinfo(sensorindex:sensorindex),& rtsolution0,options=options) ! If the CRTM returns an error flag, do not assimilate any channels for this ob ! and set the QC flag to 10 (done in setuprad). if (error_status /=0) then write(6,*)'CRTM_FORWARD ***ERROR*** during crtm_forward call ',& error_status end if endif if (trim(obstype) /= 'modis_aod' ) then ! Secant of satellite zenith angle secant_term = one/cos(data_s(ilzen_ang)) if (mixed_use) then do i=1,nchanl if (lcloud4crtm_wk(i)<0) then rtsolution(i,1) = rtsolution_clr(i,1) rtsolution_k(i,1) = rtsolution_k_clr(i,1) atmosphere_k(i,1) = atmosphere_k_clr(i,1) surface_k(i,1) = surface_k_clr(i,1) end if end do end if !$omp parallel do schedule(dynamic,1) private(i) & !$omp private(total_od,k,kk,m,term,ii,cwj) do i=1,nchanl ! Zero jacobian and transmittance arrays do k=1,nsig omix(k,i)=zero temp(k,i)=zero ptau5(k,i)=zero wmix(k,i)=zero end do ! Simulated brightness temperatures tsim(i)=rtsolution(i,1)%brightness_temperature if (n_clouds_fwd_wk>0 .and. present(tsim_clr)) then if (mixed_use) then tsim_clr(i)=rtsolution_clr(i,1)%brightness_temperature else tsim_clr(i)=rtsolution0(i,1)%brightness_temperature end if end if ! Estimated emissivity emissivity(i) = rtsolution(i,1)%surface_emissivity ! Emissivity sensitivities emissivity_k(i) = rtsolution_k(i,1)%surface_emissivity ! Surface temperature sensitivity if(nst_gsi > 1 .and. (data_s(itz_tr) > half .and. data_s(itz_tr) <= one) ) then ts(i) = surface_k(i,1)%water_temperature*data_s(itz_tr) + & surface_k(i,1)%land_temperature + & surface_k(i,1)%ice_temperature + & surface_k(i,1)%snow_temperature else ts(i) = surface_k(i,1)%water_temperature + & surface_k(i,1)%land_temperature + & surface_k(i,1)%ice_temperature + & surface_k(i,1)%snow_temperature endif if (abs(ts(i))small_wind) then term = surface_k(i,1)%wind_speed * f10*f10 / surface(1)%wind_speed uwind_k(i) = term * uu5 vwind_k(i) = term * vv5 else uwind_k(i) = zero vwind_k(i) = zero endif total_od = zero ! Accumulate values from extended into model layers ! temp - temperature sensitivity ! wmix - moisture sensitivity ! omix - ozone sensitivity ! ptau5 - layer transmittance do k=1,msig kk = klevel(msig-k+1) temp(kk,i) = temp(kk,i) + atmosphere_k(i,1)%temperature(k) wmix(kk,i) = wmix(kk,i) + atmosphere_k(i,1)%absorber(k,1) omix(kk,i) = omix(kk,i) + atmosphere_k(i,1)%absorber(k,2) total_od = total_od + rtsolution(i,1)%layer_optical_depth(k) ptau5(kk,i) = exp(-min(limit_exp,total_od*secant_term)) end do ! Load jacobian array do k=1,nsig ! Small sensitivities for temp if (abs(temp(k,i)) ! Deflate moisture jacobian above the tropopause. if (itv>=0) then do k=1,nsig jacobian(itv+k,i)=temp(k,i)*c2(k) ! virtual temperature sensitivity end do ! endif if (iqv>=0) then m=ich(i) do k=1,nsig jacobian(iqv+k,i)=c5(k)*wmix(k,i)-c4(k)*temp(k,i) ! moisture sensitivity if (prsi(k) < trop5) then term = (prsi(k)-trop5)/(trop5-prsi(nsig)) jacobian(iqv+k,i) = exp(ifactq(m)*term)*jacobian(iqv+k,i) endif end do ! endif if (ioz>=0) then ! if (.not. regional .or. use_gfs_ozone)then do k=1,nsig jacobian(ioz+k,i)=omix(k,i)*constoz ! ozone sensitivity end do ! ! end if endif if (n_clouds_fwd_wk>0 .and. n_clouds_jac_wk>0) then if (lcloud4crtm_wk(i)<=0) then do ii=1,n_clouds_jac_wk do k=1,nsig jacobian(icw(ii)+k,i) = zero end do end do else if (icmask) then do ii=1,n_clouds_jac_wk do k=1,nsig cwj(k)=zero end do do k=1,msig kk = klevel(msig-k+1) cwj(kk) = cwj(kk) + atmosphere_k(i,1)%cloud(ii)%water_content(k)*c6(k) end do do k=1,nsig jacobian(icw(ii)+k,i) = cwj(k) end do ! end do else do ii=1,n_clouds_jac_wk do k=1,nsig jacobian(icw(ii)+k,i) = zero end do ! end do endif endif endif if (ius>=0) then jacobian(ius+1,i)=uwind_k(i) ! surface u wind sensitivity endif if (ivs>=0) then jacobian(ivs+1,i)=vwind_k(i) ! surface v wind sensitivity endif if (isst>=0) then jacobian(isst+1,i)=ts(i) ! surface skin temperature sensitivity endif end do else ! obstype == 'modis_aod' ! initialize intent(out) variables that are not available with modis_aod tzbgr = zero sfc_speed = zero tsim = zero emissivity = zero ts = zero emissivity_k = zero ptau5 = zero temp = zero wmix = zero jaero = zero if(present(layer_od)) layer_od = zero if(present(jacobian_aero)) jacobian_aero = zero do i=1,nchanl do k=1,msig kk = klevel(msig-k+1) if(present(layer_od)) then layer_od(kk,i) = layer_od(kk,i) + rtsolution(i,1)%layer_optical_depth(k) endif do ii=1,n_aerosols_jac_wk if ( n_aerosols_jac_wk > n_aerosols_fwd_wk .and. ii == indx_p25 ) then jaero(kk,i,ii) = jaero(kk,i,ii) + & (0.5_r_kind*(0.78_r_kind*atmosphere_k(i,1)%aerosol(indx_dust1)%concentration(k) + & 0.22_r_kind*atmosphere_k(i,1)%aerosol(indx_dust2)%concentration(k)) ) else jaero(kk,i,ii) = jaero(kk,i,ii) + atmosphere_k(i,1)%aerosol(ii)%concentration(k) endif enddo enddo if (present(jacobian_aero)) then do k=1,nsig do ii=1,n_aerosols_jac_wk jacobian_aero(iaero_jac(ii)+k,i) = jaero(k,i,ii)*ugkg_kgm2(k) end do enddo endif enddo endif if (n_ghg >0) deallocate (tgas1d) ! contains ! pure function crtm_interface_interp(a,w,dtsig) result(intresult) ! real(r_kind), intent(in) :: a(:,:) ! real(r_kind), intent(in) :: w(:,:) ! real(r_kind), intent(in) :: dtsig ! real(r_kind) :: intresult ! integer :: i, j, n ! n = size(a,dim=1) ! intresult = 0.0_r_kind ! do j = 1, n ! do i = 1, n ! intresult = intresult + a(i,j)*w(i,j) ! enddo ! enddo ! intresult = intresult * dtsig ! end function crtm_interface_interp end subroutine call_crtm subroutine get_lai(data_s,nchanl,nreal,itime,ilate,lai_type,lai) !$$$ subprogram documentation block ! . . . . ! subprogram: get_lai interpolate vegetation LAI data for call_crtm ! ! prgmmr: ! ! abstract: ! ! program history log: ! ! input argument list: ! data_s - array containing input data information ! nchanl - number of channels ! nreal - number of descriptor information in data_s ! itime - index of analysis relative obs time ! ilate - index of earth relative latitude (degrees) ! ! output argument list: ! lai - interpolated vegetation leaf-area-index for various types (13) ! ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ !-------- use kinds, only: r_kind,i_kind use constants, only: zero use obsmod, only: iadate implicit none ! Declare passed variables integer(i_kind) ,intent(in ) :: nchanl,nreal real(r_kind),dimension(nchanl+nreal) ,intent(in ) :: data_s integer(i_kind) ,intent(in ) :: itime, ilate,lai_type real(r_kind) ,intent( out) :: lai ! Declare local variables integer(i_kind),dimension(8)::obs_time,anal_time real(r_kind),dimension(5) :: tmp_time integer(i_kind) jdow, jdoy, jday real(r_kind) rjday real(r_kind),dimension(3):: dayhf data dayhf/15.5_r_kind, 196.5_r_kind, 380.5_r_kind/ real(r_kind),dimension(13):: lai_min, lai_max data lai_min/3.08_r_kind, 1.85_r_kind, 2.80_r_kind, 5.00_r_kind, 1.00_r_kind, & 0.50_r_kind, 0.52_r_kind, 0.60_r_kind, 0.50_r_kind, 0.60_r_kind, & 0.10_r_kind, 1.56_r_kind, 0.01_r_kind / data lai_max/6.48_r_kind, 3.31_r_kind, 5.50_r_kind, 6.40_r_kind, 5.16_r_kind, & 3.66_r_kind, 2.90_r_kind, 2.60_r_kind, 3.66_r_kind, 2.60_r_kind, & 0.75_r_kind, 5.68_r_kind, 0.01_r_kind / real(r_kind),dimension(2):: lai_season real(r_kind) wei1s, wei2s integer(i_kind) n1, n2, mm, mmm, mmp ! anal_time=0 obs_time=0 tmp_time=zero tmp_time(2)=data_s(itime) anal_time(1)=iadate(1) anal_time(2)=iadate(2) anal_time(3)=iadate(3) anal_time(5)=iadate(4) call w3movdat(tmp_time,anal_time,obs_time) jdow = 0 jdoy = 0 jday = 0 call w3doxdat(obs_time,jdow,jdoy,jday) rjday=jdoy+obs_time(5)/24.0_r_kind if(rjday.lt.dayhf(1)) rjday=rjday+365.0 DO MM=1,2 MMM=MM MMP=MM+1 IF(RJDAY.GE.DAYHF(MMM).AND.RJDAY.LT.DAYHF(MMP)) THEN N1=MMM N2=MMP EXIT ENDIF if(mm == 2)PRINT *,'WRONG RJDAY',RJDAY ENDDO WEI1S = (DAYHF(N2)-RJDAY)/(DAYHF(N2)-DAYHF(N1)) WEI2S = (RJDAY-DAYHF(N1))/(DAYHF(N2)-DAYHF(N1)) IF(N2.EQ.3) N2=1 lai_season(1) = lai_min(lai_type) lai_season(2) = lai_max(lai_type) if(data_s(ilate) < 0.0_r_kind) then lai = wei1s * lai_season(n2) + wei2s * lai_season(n1) else lai = wei1s * lai_season(n1) + wei2s * lai_season(n2) endif return end subroutine get_lai end module crtm_interface