module aero_setup implicit none private public:: setup interface setup; module procedure setupaod; end interface contains subroutine setupaod(obsLL,odiagLL,lunin,mype,nchanl,nreal,nobs,& obstype,isis,is,aero_diagsave,init_pass) !$$$ subprogram documentation block ! . . . . ! subprogram: setupaod compute rhs of oi equation for aod ! prgmmr: hclin org: ncar/mmm date: 2010-10-20 ! ! abstract: read in data, first guess, and obtain rhs of oi equation ! for aod. ! ! program history log: ! 2010-10-20 hclin - modified from setuprad for aod ! 2014-01-28 todling - write sensitivity slot indicator (ioff) to header of diagfile ! 2014-12-30 derber - Modify for possibility of not using obsdiag ! 2015-10-01 guo - full res obvsr: index to allow redistribution of obsdiags ! 2015-09-10 zhu - generalize enabling all-sky and aerosol usage in radiance ! assimilation. Use radiance_obstype_search & type extentions ! 2016-02-20 pagowski - added NASA nnr AOD ! 2016-05-18 guo - replaced ob_type with polymorphic obsNode through type casting ! 2016-06-24 guo - fixed the default value of obsdiags(:,:)%tail%luse to luse(i) ! . removed (%dlat,%dlon) debris. ! 2017-02-09 guo - Remove m_alloc, n_alloc. ! . Remove my_node with corrected typecast(). ! 2018-05-19 eliu - updated crtm interface ! 2019-03-20 martin - added VIIRS AOD and ncdiag (from S-W Wei and M. Pagowski) ! ! input argument list: ! lunin - unit from which to read radiance (brightness temperature, tb) obs ! mype - mpi task id ! nchanl - number of channels per obs ! nreal - number of pieces of non-tb information per obs ! nobs - number of tb observations to process ! obstype - type of tb observation ! isis - sensor/instrument/satellite id ex.amsua_n15 ! is - integer counter for number of observation types to process ! aero_diagsave - logical to switch on diagnostic output (.false.=no output) ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use radinfo, only: nsigradjac use aeroinfo, only: nsigaerojac use crtm_interface, only: init_crtm,call_crtm,destroy_crtm,sensorindex, & isatid,itime,ilon,ilat,iszen_ang,isazi_ang use mpeu_util, only: die,perr use kinds, only: r_kind,r_single,i_kind use crtm_spccoeff, only: sc use obsmod, only: ianldate,mype_diaghdr,nchan_total, & dplat,lobsdiagsave,lobsdiag_allocated,& dirname,time_offset,luse_obsdiag use nc_diag_write_mod, only: nc_diag_init, nc_diag_header, nc_diag_metadata, & nc_diag_write, nc_diag_data2d, nc_diag_chaninfo_dim_set, nc_diag_chaninfo use nc_diag_read_mod, only: nc_diag_read_init, nc_diag_read_get_dim, nc_diag_read_close use gsi_4dvar, only: nobs_bins,hr_obsbin use gridmod, only: nsig,get_ij use constants, only: tiny_r_kind,zero,one,three,r10,max_varname_length use jfunc, only: jiter,miter use m_dtime, only: dtime_setup, dtime_check use chemmod, only: laeroana_gocart, l_aoderr_table use aeroinfo, only: jpch_aero, nusis_aero, nuchan_aero, iuse_aero, & error_aero, gross_aero use m_obsdiagNode, only: obs_diag use m_obsdiagNode, only: obs_diags use m_obsdiagNode, only: obsdiagLList_nextNode use m_obsdiagNode, only: obsdiagNode_set use m_obsdiagNode, only: obsdiagNode_get use m_obsdiagNode, only: obsdiagNode_assert use m_obsNode, only: obsNode use m_aeroNode, only: aeroNode_appendto use m_obsLList, only: obsLList use m_aeroNode, only: aeroNode, aeroNode_typecast use m_obsLList, only: obsLList_appendNode use m_obsLlist, only: obsLList_tailNode use obsmod, only: rmiss_single, netcdf_diag, binary_diag use qcmod, only: ifail_crtm_qc use radiance_mod, only: rad_obs_type,radiance_obstype_search use radiance_mod, only: n_aerosols_fwd use guess_grids, only: ntguessig,nfldsig use gsi_chemguess_mod, only: gsi_chemguess_get use gsi_bundlemod, only : gsi_bundlegetpointer use gsi_metguess_mod, only : gsi_metguess_get,gsi_metguess_bundle implicit none ! Declare passed variables type(obsLList ),target,dimension(:),intent(in):: obsLL type(obs_diags),target,dimension(:),intent(in):: odiagLL logical ,intent(in ) :: aero_diagsave character(10) ,intent(in ) :: obstype character(20) ,intent(in ) :: isis integer(i_kind) ,intent(in ) :: lunin,mype,nchanl,nreal,nobs,is logical ,intent(in ) :: init_pass ! state of "setup" processing ! Declare external calls for code analysis external:: stop2 ! Declare local parameters integer(i_kind),parameter:: ipchan=4 integer(i_kind),parameter:: ireal=5 integer(i_kind),parameter:: iversion_aerodiag=1 real(r_kind),parameter:: r1e10=1.0e10_r_kind ! Declare local variables character(128) diag_aero_file,guess_aero_file integer(i_kind) :: nvars integer(i_kind) error_status integer(i_kind) m,jc integer(i_kind) icc integer(i_kind) j,k,ncnt,i integer(i_kind) mm1 integer(i_kind) n,ibin,ioff,ioff0,iii integer(i_kind) ii,jj,idiag real(r_single) freq4,pol4,wave4,varch4 real(r_kind) errinv,useflag real(r_kind) trop5,pangs real(r_kind) cenlon,cenlat,slats,slons,dtime real(r_kind) val_obs ! Declare local arrays real(r_single),dimension(ireal):: diagbuf real(r_single),allocatable,dimension(:,:):: diagbufchan real(r_kind),dimension(nchanl):: varinv,error0 real(r_kind),dimension(nchanl):: tnoise,errmax real(r_kind),dimension(nchanl):: var,ratio_aoderr,aodinv real(r_kind),dimension(nreal+nchanl,nobs)::data_s real(r_kind),dimension(nsig):: prsltmp real(r_kind),dimension(nsig):: qsat,rh real(r_kind),dimension(nsig):: qvp,tvp real(r_kind),dimension(nsig+1):: prsitmp real(r_kind) :: psfc real(r_kind) dtsavg integer(i_kind),dimension(nchanl):: ich,id_qc real(r_kind), dimension(:,:), allocatable :: aerosols character(len=max_varname_length), dimension(:), allocatable :: & &aerosol_names character(len=56), dimension(:), allocatable :: varnames logical toss,l_may_be_passive logical,dimension(nobs):: luse integer(i_kind),dimension(nobs):: ioid ! initial (pre-distribution) obs ID integer(i_kind):: nperobs ! No. of data points, in channels, levels, or components, per obs. logical:: in_curbin, in_anybin type(aeroNode),pointer:: my_head type(obs_diag),pointer:: my_diag, obsptr type(obs_diags),pointer:: my_diagLL type(rad_obs_type) :: radmod character(len=*),parameter:: myname="setupaod" real(r_kind), dimension(nchanl) :: total_aod, aod_obs, aod integer(i_kind), parameter :: n_viirs_550nm=4 integer(i_kind) :: istyp, idbcf, ilone, ilate integer(i_kind) :: iqcall, ismask, nestat, istat real(r_kind) :: qcall, smask real(r_kind) :: styp, dbcf real(r_kind),dimension(nchanl):: emissivity,ts,emissivity_k real(r_kind),dimension(nchanl):: tsim real(r_kind),dimension(nsig,nchanl):: wmix,temp,ptau5 real(r_kind),dimension(nsigradjac,nchanl):: jacobian real(r_kind),dimension(nsigaerojac,nchanl):: jacobian_aero real(r_kind),dimension(nsig,nchanl):: layer_od real(r_kind) :: clw_guess, tzbgr, sfc_speed,ciw_guess,rain_guess,snow_guess type(obsLList),pointer,dimension(:):: aerohead aerohead => obsLL(:) if ( .not. laeroana_gocart ) then return endif !************************************************************************************** ! Initialize variables and constants. mm1 = mype+1 ncnt = 0 icc = 0 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) ilone = 5 ! index of earth relative longitude (degrees) ilate = 6 ! index of earth relative latitude (degrees) iszen_ang = 8 ! index of solar zenith angle (degrees) isazi_ang = 9 ! index of solar azimuth angle (degrees) istyp = 10 ! index of surface type idbcf = 11 ! index of deep blue confidence flag if ( obstype == 'viirs_aod' .or. obstype == 'modis_aod' ) then iqcall = 7 ! index of overall quality flag for AOD ismask = 10 ! index of surface type mask else ! obstype /= 'modis_aod' or 'viirs_aod' write(6,*)'SETUP_AOD: *** WARNING: unknown aerosol input type, obstype=',obstype end if ! Determine cloud & aerosol usages in radiance assimilation call radiance_obstype_search(obstype,radmod) ! Initialize channel related information tnoise = r1e10 errmax = r1e10 l_may_be_passive = .false. toss = .true. jc=0 do j=1,jpch_aero if(isis == nusis_aero(j))then jc=jc+1 if(jc > nchanl)then write(6,*)'setupaod: ***ERROR*** in channel numbers, jc,nchanl=',jc,nchanl,& ' ***STOP IN setupaod***' call stop2(71) end if ! Load channel numbers into local array based on satellite type ich(jc)=j ! ! Set error instrument channels tnoise(jc)=error_aero(j) errmax(jc)=gross_aero(j) if (iuse_aero(j)< -1 .or. (iuse_aero(j) == -1 .and. & .not.aero_diagsave)) tnoise(jc)=r1e10 if (iuse_aero(j)>-1) l_may_be_passive=.true. if (tnoise(jc) < 1.e4_r_kind) toss = .false. end if end do if ( mype == 0 .and. .not.l_may_be_passive) write(6,*)mype,'setupaod: passive obs',is,isis if(nchanl > jc) write(6,*)'setupaod: channel number reduced for ', & obstype,nchanl,' --> ',jc if(jc == 0) then if(mype == 0) write(6,*)'setupaod: No channels found for ', & obstype,isis if(nobs > 0)read(lunin) return end if if (toss) then if(mype == 0)write(6,*)'setupaod: all obs var > 1e4. do not use ',& 'data from satellite is=',isis if(nobs >0)read(lunin) return endif ioff0=0 if (lobsdiagsave) then if (l_may_be_passive) then ioff0=4 else ioff0=5 endif endif ! Initialize radiative transfer call init_crtm(init_pass,mype_diaghdr(is),mype,nchanl,nreal,isis,obstype,radmod) ! If diagnostic file requested, allocate arrays and init output file if (aero_diagsave) then allocate(aerosols(nsig,n_aerosols_fwd),aerosol_names(n_aerosols_fwd)) nvars=5+n_aerosols_fwd allocate(varnames(nvars)) call gsi_chemguess_get('aerosols::3d',aerosol_names,istat) varnames(1:5) = (/'air_temperature ','humidity_mixing_ratio', & 'relative_humidity ','air_pressure ','air_pressure_levels '/) varnames(6:) = aerosol_names if (binary_diag) call init_binary_diag_ if (netcdf_diag) call init_netcdf_diag_ end if idiag=ipchan if (lobsdiagsave) idiag=idiag+4*miter+1 allocate(diagbufchan(idiag,nchanl)) ! Load data array for current satellite read(lunin) data_s,luse,ioid write(*,*) 'read in AOD data ',nobs ! Loop over data in this block call dtime_setup() do n = 1,nobs ! Extract analysis relative observation time. dtime = data_s(itime,n) call dtime_check(dtime, in_curbin, in_anybin) if(.not.in_anybin) cycle if(in_curbin) then id_qc = 0 ! Extract lon and lat. slons = data_s(ilon,n) ! grid relative longitude slats = data_s(ilat,n) ! grid relative latitude cenlon = data_s(ilone,n) ! earth relative longitude (degrees) cenlat = data_s(ilate,n) ! earth relative latitude (degrees) pangs = data_s(iszen_ang,n) if ( obstype == 'modis_aod' ) then styp = data_s(istyp,n) dbcf = data_s(idbcf,n) else if ( obstype == 'viirs_aod' ) then qcall = data_s(iqcall,n) smask = data_s(ismask,n) end if ! Set relative weight value val_obs=one ! Load channel data into work array. aod_obs = rmiss_single do i = 1, nchanl ! fix channel issue for VIIRS except channel 4 if (obstype == 'viirs_aod' .and. i /= n_viirs_550nm) cycle aod_obs(i) = data_s(i+nreal,n) end do if ( .not. l_aoderr_table ) then ! set observation error if ( obstype == 'modis_aod' ) then select case ( nint(styp) ) case ( 0 ) ! water tnoise = 0.03_r_kind+0.05_r_kind*aod_obs case ( 1, 2, 3 ) ! coast, desert, land tnoise = 0.05_r_kind+0.15_r_kind*aod_obs case ( 4 ) ! deep blue if ( nint(dbcf) >= 0 .and. nint(dbcf) <= 3 ) then tnoise = 0.05_r_kind+0.15_r_kind*aod_obs+0.01_r_kind*(three-dbcf) end if case ( 5 ) ! nnr ocean tnoise = 0.2_r_kind*(aod_obs+0.01_r_kind) case ( 6 ) ! nnr land tnoise = 0.2_r_kind*(aod_obs+0.01_r_kind) end select else if ( obstype == 'viirs_aod' ) then nestat = nint(qcall)+nint(smask)*10 select case (nestat) case( 2 ) ! over water surface, medium-quality tnoise = 0.0416146_r_kind+0.0808841_r_kind*aod_obs case( 3 ) ! over water surface, high quality tnoise = 0.00784394_r_kind+0.219923_r_kind*aod_obs case( 12 ) ! over dark land surface, medium-quality tnoise = 0.0374849_r_kind+0.266073_r_kind*aod_obs case( 13 ) ! over dark land surface, high quality tnoise = 0.111431_r_kind+0.128699_r_kind*aod_obs case( 22 ) ! over bright land surface, medium-quality tnoise = 0.0693246_r_kind+0.270070_r_kind*aod_obs case( 23 ) ! over bright land surface, high quality tnoise = 0.0550472_r_kind+ 0.299558_r_kind*aod_obs end select else if (mype == 0) then write(6,*),'unknown obstype = ',obstype call stop2(283) end if end if ! end if obstype end if ! end if not l_aoderr_table ! Interpolate model fields to observation location, call crtm and create jacobians call call_crtm(obstype,dtime,data_s(:,n),nchanl,nreal,ich, & tvp,qvp,clw_guess,ciw_guess,rain_guess,snow_guess,prsltmp,prsitmp, & trop5,tzbgr,dtsavg,sfc_speed, & tsim,emissivity,ptau5,ts,emissivity_k, & temp,wmix,jacobian,error_status,layer_od=layer_od,jacobian_aero=jacobian_aero) ! interpolate aerosols at observation locations for diag files here if (aero_diagsave) then call genqsat(qsat,tvp,prsltmp,1,1,nsig,.true.,0) rh = qvp/qsat call aero_guess_at_obs_locations(dtime,data_s(:,n),& nchanl,nreal,nsig, n_aerosols_fwd, aerosols, aerosol_names) end if ! If the CRTM returns an error flag, do not assimilate any channels for this ob ! and set the QC flag to ifail_crtm_qc. ! We currently go through the rest of the QC steps, ensuring that the diagnostic ! files are populated, but this could be changed if it causes problems. if (error_status /=0) then id_qc(1:nchanl) = ifail_crtm_qc varinv(1:nchanl) = zero endif total_aod = zero do i = 1, nchanl total_aod(i) =sum(layer_od(:,i)) enddo do i = 1, nchanl aod(i) = aod_obs(i) - total_aod(i) error0(i) = tnoise(i) if(aod_obs(i)>zero .and. tnoise(i) < 1.e4_r_kind .or. (iuse_aero(ich(i))==-1 & .and. aero_diagsave))then varinv(i) = val_obs/tnoise(i)**2 else if(id_qc(i) == 0)id_qc(i)=1 varinv(i) = zero endif end do icc = 0 do i = 1, nchanl ! Only process observations to be assimilated if (varinv(i) > tiny_r_kind ) then m = ich(i) ! Only "good" obs are included in J calculation. if (iuse_aero(m) >= 1)then icc = icc + 1 aodinv(icc) = aod(i) ! obs-ges innovation var(icc) = one/error0(i)**2 ! 1/(obs error)**2 (original uninflated error) ratio_aoderr(icc)=error0(i)**2*varinv(i) ! (original error)/(inflated error) endif endif end do endif ! (in_curbin) ! In principle, we want ALL obs in the diagnostics structure but for ! passive obs (monitoring), it is difficult to do if aero_diagsave ! is not on in the first outer loop. For now we use l_may_be_passive... if (l_may_be_passive) then ! Link observation to appropriate observation bin if (nobs_bins>1) then ibin = NINT( dtime/hr_obsbin ) + 1 else ibin = 1 endif if (ibin<1.OR.ibin>nobs_bins) write(6,*)mype,'Error nobs_bins,ibin= ',nobs_bins,ibin if(luse_obsdiag) my_diagLL => odiagLL(ibin) if (in_curbin) then ! Load data into output arrays if (icc > 0) then ncnt =ncnt+1 nchan_total=nchan_total+icc allocate(my_head) call aeroNode_appendto(my_head,aerohead(ibin)) my_head%idv = is my_head%iob = ioid(n) my_head%elat= data_s(ilate,n) my_head%elon= data_s(ilone,n) allocate(my_head%res(icc),my_head%err2(icc), & my_head%raterr2(icc), & my_head%daod_dvar(nsigaerojac,icc), & my_head%ich(icc),& my_head%icx(icc)) if(luse_obsdiag)allocate (my_head%diags(icc)) my_head%nlaero = icc ! profile observation count call get_ij(mm1,slats,slons,my_head%ij,my_head%wij) my_head%time=dtime my_head%luse=luse(n) my_head%ich(:)=-1 iii=0 do ii=1,nchanl m=ich(ii) if (varinv(ii)>tiny_r_kind .and. iuse_aero(m)>=1) then iii=iii+1 my_head%res(iii)=aodinv(iii) my_head%err2(iii)=var(iii) my_head%raterr2(iii)=ratio_aoderr(iii) my_head%icx(iii)=m do k = 1, nsigaerojac my_head%daod_dvar(k,iii)=jacobian_aero(k,ii) end do my_head%ich(iii)=ii end if end do my_head => null() end if ! icc endif ! (in_curbin) ! Link obs to diagnostics structure if(luse_obsdiag) then iii=0 obsptr => null() do ii=1,nchanl nperobs=-99999; if(ii==1) nperobs=nchanl my_diag => obsdiagLList_nextNode(my_diagLL ,& create = .not.lobsdiag_allocated ,& idv = is ,& iob = ioid(n) ,& ich = ii ,& elat = data_s(ilate,n) ,& elon = data_s(ilone,n) ,& luse = luse(n) ,& miter = miter ) if(.not.associated(my_diag)) call die(myname,'not associated(my_diag)') if (ii==1) obsptr => my_diag ! this is the lead node if (in_curbin.and.icc>0) then my_head => tailNode_typecast_(aerohead(ibin)) if(.not.associated(my_head)) & call die(myname,'unexpected, associated(my_head) =',associated(my_head)) call obsdiagNode_set(my_diag, wgtjo=varinv(ii), jiter=jiter, nldepart=aod(ii) ) ! Load data into output arrays m=ich(ii) if (varinv(ii)>tiny_r_kind .and. iuse_aero(m)>=1) then iii=iii+1 call obsdiagNode_assert(my_diag,my_head%idv,my_head%iob,my_head%ich(iii),myname,'my_diag:my_head error') call obsdiagNode_set(my_diag, jiter=jiter, muse=.true.) my_head%diags(iii)%ptr => my_diag endif my_head => null() endif ! (in_curbin) enddo ! do ii=1,nchanl if (in_curbin) then if( iii/=icc ) then write(6,*)'setupaod: error iii icc',iii,icc call stop2(279) endif endif ! (in_curbin) endif ! (luse_obsdiag) ! End of l_may_be_passive block endif if(in_curbin) then ! Write diagnostics to output file. if (aero_diagsave .and. luse(n)) then diagbuf(1) = cenlat ! observation latitude (degrees) diagbuf(2) = cenlon ! observation longitude (degrees) diagbuf(3) = dtime-time_offset ! observation time (hours relative to analysis time) diagbuf(4) = pangs ! solar zenith angle (degrees) diagbuf(5) = data_s(isazi_ang,n) ! solar azimuth angle (degrees) do i=1,nchanl diagbufchan(1,i)=aod_obs(i) ! observed brightness temperature (K) ! diagbufchan(2,i)=total_aod(i) ! observed - simulated Tb with no bias corrrection (K) - this should be innovation diagbufchan(2,i)=aod(i) ! innovation errinv = sqrt(varinv(i)) diagbufchan(3,i)=errinv ! inverse observation error useflag=one if (iuse_aero(ich(i)) < 1) useflag=-one diagbufchan(4,i)= id_qc(i)*useflag! quality control mark or event indicator end do if (lobsdiagsave) then if (l_may_be_passive) then do ii=1,nchanl if (.not.associated(obsptr)) then write(6,*)'setupaod: error obsptr' call stop2(280) end if ioff=ioff0 do jj=1,miter ioff=ioff+1 if (obsptr%muse(jj)) then diagbufchan(ioff,ii) = one else diagbufchan(ioff,ii) = -one endif enddo do jj=1,miter+1 ioff=ioff+1 diagbufchan(ioff,ii) = obsptr%nldepart(jj) enddo do jj=1,miter ioff=ioff+1 diagbufchan(ioff,ii) = obsptr%tldepart(jj) enddo do jj=1,miter ioff=ioff+1 diagbufchan(ioff,ii) = obsptr%obssen(jj) enddo obsptr => obsptr%next enddo else ioff=ioff0 diagbufchan(ioff+1:ioff+4*miter+1,1:nchanl) = zero endif endif psfc=prsitmp(1)*r10 ! convert to hPa write(4) psfc,diagbuf,diagbufchan if (binary_diag) call contents_binary_diag_ if (netcdf_diag) call contents_netcdf_diag_ end if endif ! (in_curbin) 100 continue ! End of n-loop over obs end do ! Jump here when there is no data to process for current satellite ! Deallocate arrays deallocate(diagbufchan) if (aero_diagsave) then close(4) if (binary_diag) call final_binary_diag_ if (netcdf_diag) call nc_diag_write endif call destroy_crtm ! End of routine return contains function tailNode_typecast_(oll) result(ptr_) !> Cast the tailNode of oll to an aeroNode, as in !> ptr_ => typecast_(tailNode_(oll)) use m_aeroNode, only: aeroNode, typecast_ => aeroNode_typecast use m_obsLList, only: obsLList, tailNode_ => obsLList_tailNode use m_obsNode , only: obsNode implicit none type(aeroNode),pointer:: ptr_ type(obsLList),target ,intent(in):: oll class(obsNode),pointer:: inode_ inode_ => tailNode_(oll) ptr_ => typecast_(inode_) end function tailNode_typecast_ subroutine init_binary_diag_ ! subroutine to initialize binary diag files ! original: pagowski ! modified: 2019-03-20 - martin - cleaned up to fit GSI coding norms implicit none character(10) :: filex character(12) :: string filex=obstype write(string,1976) jiter 1976 format('_',i2.2) diag_aero_file= trim(dirname) // trim(filex) // '_' // trim(dplat(is)) // trim(string) guess_aero_file= trim(dirname) // trim(filex) // '_' // trim(dplat(is)) // '_vars' // trim(string) if(init_pass) then open(4,file=trim(diag_aero_file),form='unformatted',status='unknown',position='rewind') open(41,file=trim(guess_aero_file),form='unformatted',status='unknown',position='rewind') else open(4,file=trim(diag_aero_file),form='unformatted',status='old',position='append') open(41,file=trim(guess_aero_file),form='unformatted',status='old',position='append') endif ! Initialize/write parameters for satellite diagnostic file on ! first outer iteration. if (init_pass .and. mype==mype_diaghdr(is)) then write(4)isis,dplat(is),obstype,jiter,nchanl,ianldate,ireal,ipchan,nsig,ioff0 write(41)nsig,nvars,n_aerosols_fwd,ianldate write(41)varnames write(6,*)'setupaod: write header record for ',& isis,ireal,' to file ',trim(diag_aero_file),' ',ianldate do i=1,nchanl n=ich(i) if( iuse_aero(n) < 0 ) cycle !if( n < 1 )cycle varch4=error_aero(n) freq4=sc(sensorindex)%frequency(i) pol4=sc(sensorindex)%polarization(i) wave4=sc(sensorindex)%wavenumber(i) write(4)freq4,pol4,wave4,varch4,iuse_aero(n),& nuchan_aero(n),ich(i) end do end if end subroutine init_binary_diag_ subroutine init_netcdf_diag_ ! subroutine to initialize netcdf diag files ! original: pagowski ! modified: 2019-03-21 - martin - cleaned up to fit GSI coding norms implicit none character(10) :: filex character(12) :: string filex=obstype write(string,1976) jiter 1976 format('_',i2.2) diag_aero_file= trim(dirname) // trim(filex) // '_' // trim(dplat(is)) //trim(string) // '.nc4' if (init_pass .and. nobs > 0) then call nc_diag_init(diag_aero_file) call nc_diag_chaninfo_dim_set(nchanl) end if if (init_pass) then call nc_diag_header("Satellite_Sensor", isis ) call nc_diag_header("Satellite", dplat(is) ) call nc_diag_header("Observation_type", "aod" ) call nc_diag_header("Number_of_channels", nchanl ) call nc_diag_header("date_time", ianldate ) do i=1,nchanl n=ich(i) if( iuse_aero(n) < 0 ) cycle call nc_diag_chaninfo("frequency",sngl(sc(sensorindex)%frequency(i))) call nc_diag_chaninfo("polarization",sc(sensorindex)%polarization(i)) call nc_diag_chaninfo("wavenumber",sngl(sc(sensorindex)%wavenumber(i))) call nc_diag_chaninfo("use_flag", iuse_aero(n)) call nc_diag_chaninfo("sensor_chan", nuchan_aero(n)) end do end if end subroutine init_netcdf_diag_ subroutine contents_binary_diag_ ! subroutine to write contents to binary diag files ! original: pagowski ! modified: 2019-03-21 - martin - cleaned up to fit GSI coding norms implicit none diagbuf(1) = cenlat ! observation latitude (degrees) diagbuf(2) = cenlon ! observation longitude (degrees) diagbuf(3) = dtime!-time_offset ! observation time (hours relative to analysis time) diagbuf(4) = pangs ! solar zenith angle (degrees) diagbuf(5) = data_s(isazi_ang,n) ! solar azimuth angle (degrees) do i=1,nchanl diagbufchan(1,i)=aod_obs(i) ! observed brightness temperature (K) diagbufchan(2,i)=aod(i) ! innovation errinv = sqrt(varinv(i)) diagbufchan(3,i)=errinv ! inverse observation error useflag=one if (iuse_aero(ich(i)) < 1) useflag=-one diagbufchan(4,i)= id_qc(i)*useflag! quality control mark or event indicator end do if (lobsdiagsave) then if (l_may_be_passive) then do ii=1,nchanl if (.not.associated(obsptr)) then write(6,*)'setupaod: error obsptr' call stop2(280) end if ioff=ioff0 do jj=1,miter ioff=ioff+1 if (obsptr%muse(jj)) then diagbufchan(ioff,ii) = one else diagbufchan(ioff,ii) = -one endif enddo do jj=1,miter+1 ioff=ioff+1 diagbufchan(ioff,ii) = obsptr%nldepart(jj) enddo do jj=1,miter ioff=ioff+1 diagbufchan(ioff,ii) = obsptr%tldepart(jj) enddo do jj=1,miter ioff=ioff+1 diagbufchan(ioff,ii) = obsptr%obssen(jj) enddo obsptr => obsptr%next enddo else ioff=ioff0 diagbufchan(ioff+1:ioff+4*miter+1,1:nchanl) = zero endif endif write(4) diagbuf,diagbufchan write(41)real(tvp,r_single),real(qvp/(one-qvp),r_single),& &real(rh,r_single),& &real(prsltmp,r_single),real(prsitmp,r_single) write(41)real(aerosols,r_single) end subroutine contents_binary_diag_ subroutine contents_netcdf_diag_ ! subroutine to write contents to netcdf diag files ! original: pagowski ! modified: 2019-03-21 - martin - cleaned up to fit GSI coding norms implicit none character(7),parameter :: obsclass = ' aod' character(128) :: fieldname integer(i_kind) :: iaero,k,l real(r_single), dimension(nsig+1) :: tmp real(r_single),parameter:: missing = -9.99e9_r_single do i=1,nchanl l=ich(i) if ( iuse_aero(l) < 0 ) cycle call nc_diag_metadata("Channel_Index", i) call nc_diag_metadata("Observation_Class", obsclass) call nc_diag_metadata("Latitude", sngl(cenlat)) ! observation latitude (degrees) call nc_diag_metadata("Longitude", sngl(cenlon)) ! observation longitude (degrees) call nc_diag_metadata("Obs_Time", sngl(dtime))!-time_offset)) ! observation time (hours relative to analysis time) call nc_diag_metadata("Sol_Zenith_Angle", sngl(pangs)) ! solar zenith angle (degrees) call nc_diag_metadata("Sol_Azimuth_Angle", sngl(data_s(isazi_ang,n))) ! solar azimuth angle (degrees) call nc_diag_metadata("Surface_type", nint(data_s(istyp,n))) call nc_diag_metadata("MODIS_deep_blue_flag", nint(dbcf) ) call nc_diag_metadata("Observation", sngl(diagbufchan(1,i)) ) ! observed aod call nc_diag_metadata("Obs_Minus_Forecast_adjusted",sngl(diagbufchan(2,i))) call nc_diag_metadata("Obs_Minus_Forecast_unadjusted",sngl(diagbufchan(2,i)))! obs - sim aod with no bias correction if (diagbufchan(3,i) > tiny_r_kind) then tmp(1)=one/diagbufchan(3,i) else tmp(1)=missing end if call nc_diag_metadata("Observation_Error",tmp(1)) call nc_diag_metadata("QC_Flag", sngl(diagbufchan(4,i))) !quality control mark or event indicator tmp(1)=get_zsfc() call nc_diag_metadata("sfc_height",tmp(1)) ! height in meters do k=1,nsig tmp(k)=tvp(nsig-k+1) end do call nc_diag_data2d("air_temperature", tmp(1:nsig)) ! K do k=1,nsig tmp(k)=qvp(nsig-k+1)/(1_r_kind-qvp(nsig-k+1)) end do call nc_diag_data2d("humidity_mixing_ratio", tmp(1:nsig)) ! kg/kg do k=1,nsig tmp(k)=rh(nsig-k+1) end do call nc_diag_data2d("relative_humidity", tmp(1:nsig)) ! 0-1 do k=1,nsig tmp(k)=1000_r_single*prsltmp(nsig-k+1) end do call nc_diag_data2d("air_pressure", tmp(1:nsig)) ! Pa do k=1,nsig+1 tmp(k)=1000_r_single*prsitmp(nsig-k+2) end do call nc_diag_data2d("air_pressure_levels", tmp(1:nsig+1)) ! Pa do iaero = 1, n_aerosols_fwd write (fieldname, "(A,I0.2)") aerosol_names(iaero) do k=1,nsig tmp(k)=aerosols(nsig-k+1,iaero) end do call nc_diag_data2d(trim(fieldname), tmp(1:nsig)) !mixing ratios in ug/kg end do end do end subroutine contents_netcdf_diag_ subroutine final_binary_diag_ ! subroutine to finalize binary diag files ! original: pagowski ! modified: 2019-03-21 - martin - cleaned up to fit GSI coding norms close(4) close(41) end subroutine final_binary_diag_ ! nc_diag_write is a generic routine that takes care of finalizing the netcdf diag file ! so no need for final_netcdf_diag_ subroutine function get_zsfc() RESULT(zsfc) ! function to get surface height from GSI bundle ! original: pagowski ! modified: 2019-03-21 - martin - cleaned up to fit GSI coding norms implicit none real(r_kind) :: zsfc real(r_kind),dimension(:,: ),pointer:: rank2 character(len=5) :: varname integer(i_kind) :: istatus,ifld real(r_kind),allocatable,dimension(:,:,: ) :: ges_z varname='z' call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname)& &,rank2,istatus) if (istatus==0) then if(allocated(ges_z)) then write(6,*) trim(myname), ': ', trim(varname), ' already& & incorrectly allocated ' call stop2(111) end if allocate(ges_z(size(rank2,1),size(rank2,2),nfldsig)) ges_z(:,:,1)=rank2 do ifld=2,nfldsig call gsi_bundlegetpointer(gsi_metguess_bundle(ifld)& &,trim(varname),rank2,istatus) ges_z(:,:,ifld)=rank2 end do call intrp2a11(ges_z(1,1,ntguessig),zsfc,slats,slons,mype) else write(6,*) trim(myname),': ', trim(varname), ' not found in& & met bundle, ier= ',istatus call stop2(112) end if end function get_zsfc end subroutine setupaod end module aero_setup