module oneobmod
!$$$ module documentation block
! . . . .
! module: oneobmod
! prgmmr: kleist org: np20 date: 2003-10-20
!
! abstract: module contains everything necessary for running single
! observation experiments
!
! program history log:
! 2003-10-20 kleist
! 2004-05-13 kleist, documentation
! 2005-04-04 todling, fixed little endian ouput of prepqc file
! 2009-04-28 sienkiewicz - add text output for ozone level obs testing
! 2012-07-14 todling - only do it once (in observer mode)
! 2014-05-29 thomas - add lsingleradob parameter for single radiance
! assimilation (originally of mccarty)
! 2016-08-12 lippi - added nml parameters for single radial wind
! experiment (anaz_rw,anel_rw,range_rw,sstn,lsingleradar,
! singleradar,learthrel_rw). Added use of radar look-up
! table. Added oneobmakebufr and invtllv subroutines to make
! the radial wind superobs. Added nml options for a
! "single radar" observation test (all obs from one
! radar).
! 2016-12-14 lippi - added nml option learthrel_rw to not rotate the winds
! from lat lon to xy.
!
! subroutines included:
! sub init_oneobmod
! sub oneobmakebufr
! sub oneobmakerwsupob
! sub oneobo3lv
! sub invtllv
!
! variable definitions:
! def maginnov - magnitude of innovation for one ob exp
! def magoberr - magnitude of observational error for one ob exp
! def oblat - observation latitude for one ob exp
! def oblon - observation longitude for one ob exp
! def obhourset - observation delta time from analysis time for
! one ob exp
! def obpres - observation pressure (hPa) or one ob exp
! def obdattim - observation date for one ob exp
! def oneob_type - observation type for one ob exp
! def oneobtest - single observation test flag (true=on)
! def pctswitch - if true, innovation and error expressed as percentage
! of background value
!
!$$$
use kinds, only: r_kind,i_kind,r_double
implicit none
! set default to private
private
! set subroutines to public
public :: init_oneobmod
public :: oneobmakebufr
public :: oneobmakerwsupob
public :: oneobo3lv
! set passed variables to public
public :: oneobtest,oneob_type,magoberr,pctswitch,maginnov
public :: oblat,oblon,obpres,obdattim,obhourset,anel_rw,anaz_rw,range_rw
public :: lsingleradob, obchan, sstn
public :: lsingleradar, singleradar
public :: learthrel_rw
public :: virtmp
real(r_kind) maginnov, magoberr, oblat, oblon,&
obhourset, obpres,anel_rw,anaz_rw,range_rw
integer(i_kind) obdattim
character(10) oneob_type
character(4) sstn, singleradar
logical oneobtest
logical pctswitch
logical lsingleradob, lsingleradar, learthrel_rw
logical virtmp
integer(i_kind) obchan
logical :: oneobmade
contains
subroutine init_oneobmod
!$$$ subprogram documentation block
! . . . .
! subprogram: init_oneobmod
! prgmmr: kleist org: np20 date: 2003-10-20
!
! abstract: initialize defaults for single ob experiment vars
!
! program history log:
! 2003-10-20 kleist
! 2004-05-13 kleist, documentation
!
! input argument list:
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!$$$
use constants, only: zero, one, r1000
implicit none
oneobtest=.false.
maginnov=one
magoberr=one
oneob_type=' '
oblat=zero
oblon=zero
obpres=r1000
obdattim=2000010100
obhourset=zero
pctswitch=.false.
lsingleradob=.false.
obchan=zero
sstn='KOUN'
anel_rw=zero
anaz_rw=zero
range_rw=r1000
virtmp=.false.
oneobmade=.false.
lsingleradar=.false.
learthrel_rw=.false.
singleradar='KOUN'
return
end subroutine init_oneobmod
subroutine oneobmakebufr
!$$$ subprogram documentation block
! . . . .
! subprogram: oneobmakebufr
! prgmmr: kleist org: np20 date: 2003-10-20
!
! abstract: create prepbufr file for single ob experiment
!
! program history log:
! 2003-10-20 kleist
! 2004-05-13 kleist documentation
! 2006-04-06 middlecoff - changed lumk from 52 to lendian_in so one-obs prepqc
! file can be read as little endian
! 2014-08-04 carley - modify for tcamt and howv obs
! 2014-08-18 carley - added td2m, mxtm, mitm, pmsl, and wsdp10m
! 2016-01-18 pondeca - added cldch
! 2016-06-17 Sienkiewicz- virtmp switch - if true write VIRTMP program code to indicate
! that the observation is virtual temperature
!
! input argument list:
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!$$$
use constants, only: zero, one, five, one_tenth, r100, r0_01
use gsi_io, only: lendian_in
use obsmod, only: offtime_data,iadate,bmiss
implicit none
real(r_kind),parameter:: r20=20.0_r_kind
integer(i_kind) ludx,nobs,nlev,idate
character(8) subset,sid(1)
real(r_kind),dimension(1):: typ
real(r_kind),dimension(1,1):: qob,tob,zob,uob,vob,cat
real(r_kind),dimension(1,1):: pqm,qqm,tqm,zqm,wqm
real(r_kind),dimension(1,1):: poe,qoe,toe,woe
real(r_kind),dimension(1):: xob,yob,dhr
real(r_kind),dimension(1,1):: pob
real(r_double) vtcd
integer(i_kind) n,k,iret
real(r_kind) hdr(10),obs(13,255),qms(10,255),err(10,255),cld2seq(2,1), &
cldseq(3,10),owave(1,255),maxtmint(2,255),cldceilh(1,255),&
pcd(10,255)
character(80):: hdrstr='SID XOB YOB DHR TYP'
character(80):: obsstr='POB QOB TOB ZOB UOB VOB CAT PWO MXGS HOVI PRSS TDO PMO'
character(80):: qmsstr='PQM QQM TQM ZQM WQM PWQ PMQ'
character(80):: errstr='POE QOE TOE WOE'
character(80):: pcdstr='PPC QPC TPC WPC'
character(80):: cld2seqstr='TOCC HBLCS' ! total cloud cover and height above surface of base of lowest cloud seen
character(80):: cldseqstr='VSSO CLAM HOCB' ! vertical significance, cloud amount and cloud base height
character(80):: owavestr='HOWV' ! significant wave height
character(80):: maxtmintstr='MXTM MITM' ! Max T and Min T
! character(80):: cldceilhstr='CEILING' ! cldch
if (oneobmade) return
if (oneob_type == 'o3lev') then
call oneobo3lv
return
end if
! set values from parameter list
xob=oblon
yob=oblat
dhr=obhourset
idate=iadate(1)*1000000+iadate(2)*10000+iadate(3)*100+iadate(4)
write(6,*)'OneObMake: ', idate
pob=obpres
! set default values for this routine
ludx=22
nobs=1
nlev=1
subset='ADPUPA'
sid='SID00001'
qob=r100
tob=r20
zob=zero
uob=five
vob=five
owave=bmiss
maxtmint=bmiss
cldceilh=bmiss
cld2seq=bmiss
cldseq=bmiss
pqm=one
qqm=one
tqm=one
zqm=one
wqm=one
offtime_data=.true.
if (oneob_type=='ps') then
typ(1)=87._r_kind
cat(1,1)=zero
else if (oneob_type=='tcamt') then
subset='ADPSFC'
typ(1)=87._r_kind
cat(1,1)=zero
cld2seq(1,1)=25._r_kind !TOCC - total cloud amount (%)
else if (oneob_type=='howv') then
subset='SFCSHP'
typ(1)=80._r_kind
cat(1,1)=zero
owave(1,1)=4._r_kind !Significant wave height (m - includes wind+swell waves)
else if (oneob_type=='td2m') then
subset='ADPSFC'
typ(1)=87._r_kind
cat(1,1)=zero
obs(12,1)=280.0_r_kind !Dewpoint in Kelvin
else if (oneob_type=='mxtm') then
subset='ADPSFC'
typ(1)=87._r_kind
cat(1,1)=zero
maxtmint(1,1)=280.0_r_kind !Max T in Kelvin
else if (oneob_type=='mitm') then
subset='ADPSFC'
typ(1)=87._r_kind
cat(1,1)=zero
maxtmint(2,1)=273.15_r_kind !Min T in Kelvin
else if (oneob_type=='pmsl') then
subset='ADPSFC'
typ(1)=87._r_kind
cat(1,1)=zero
obs(13,1)=1008.10_r_kind !PMSL in MB
qms(7,1)=one
else if (oneob_type=='wspd10m') then !10m AGL wind speed - need to store both u and v (already by this point done)
subset='ADPSFC'
typ(1)=87._r_kind
cat(1,1)=zero
else if (oneob_type=='cldch') then
subset='ADPSFC'
typ(1)=87._r_kind
cat(1,1)=zero
cldceilh(1,1)=1000._r_kind !clch in m
else
typ(1)=20._r_kind
cat(1,1)=one
endif
! keep errs small so the single ob passes the QC check
poe=r0_01
qoe=one_tenth
toe=one_tenth
woe=one_tenth
open(ludx,file='prepobs_prep.bufrtable',action='read')
#if defined(__osf__) || defined(__ia64__) && (__INTEL_COMPILER>799)
open(lendian_in,file='prepqc',action='write',form='unformatted',convert='little_endian')
#else
open(lendian_in,file='prepqc',action='write',form='unformatted')
#endif
call datelen(10)
call openbf(lendian_in,'OUT',ludx)
if (virtmp) then
call ufbqcd(lendian_in,'VIRTMP',vtcd)
end if
do n=1,nobs
hdr(1)=transfer(sid(n),hdr(1))
hdr(2)=xob(n)
hdr(3)=yob(n)
hdr(4)=dhr(n)
hdr(5)=r100+typ(n)
obs=bmiss
qms=bmiss
err=bmiss
pcd=bmiss
do k=1,nlev
obs(1,k)=pob(k,n)
obs(2,k)=qob(k,n)
obs(3,k)=tob(k,n)
obs(4,k)=zob(k,n)
if (virtmp) then
pcd(3,k)=vtcd
end if
obs(7,k)=cat(k,n)
qms(1,k)=pqm(k,n)
qms(2,k)=qqm(k,n)
qms(3,k)=tqm(k,n)
qms(4,k)=zqm(k,n)
err(1,k)=poe(k,n)
err(2,k)=qoe(k,n)
err(3,k)=toe(k,n)
enddo
call openmb(lendian_in,subset,idate)
call ufbint(lendian_in,hdr,10,1,iret,hdrstr)
call ufbint(lendian_in,obs,13,nlev,iret,obsstr)
call ufbint(lendian_in,qms,10,nlev,iret,qmsstr)
call ufbint(lendian_in,err,10,nlev,iret,errstr)
call ufbint(lendian_in,pcd,10,nlev,iret,pcdstr)
if (oneob_type=='tcamt') then
call ufbint(lendian_in,cldseq,3,10,iret,cldseqstr)
call ufbint(lendian_in,cld2seq,2,1,iret,cld2seqstr)
else if (oneob_type=='howv') then
call ufbint(lendian_in,owave,1,nlev,iret,owavestr)
else if ( oneob_type=='mxtm' .or. oneob_type=='mitm') then
call ufbint(lendian_in,maxtmint,2,nlev,iret,maxtmintstr)
end if
call writsb(lendian_in)
hdr(1)=transfer(sid(n),hdr(1))
hdr(2)=xob(n)
hdr(3)=yob(n)
hdr(4)=dhr(n)
hdr(5)=200_r_kind+typ(n)
obs=bmiss
qms=bmiss
err=bmiss
pcd=bmiss
do k=1,nlev
obs(1,k)=pob(k,n)
obs(5,k)=uob(k,n)
obs(6,k)=vob(k,n)
obs(7,k)=cat(k,n)
qms(1,k)=pqm(k,n)
qms(5,k)=wqm(k,n)
err(1,k)=poe(k,n)
err(4,k)=woe(k,n)
enddo
call openmb(lendian_in,subset,idate)
call ufbint(lendian_in,hdr,10,1,iret,hdrstr)
call ufbint(lendian_in,obs,13,nlev,iret,obsstr)
call ufbint(lendian_in,qms,10,nlev,iret,qmsstr)
call ufbint(lendian_in,err,10,nlev,iret,errstr)
call writsb(lendian_in)
enddo
call closbf(lendian_in)
oneobmade=.true.
return
end subroutine oneobmakebufr
subroutine oneobmakerwsupob
! parts borrowed from subroutine radar_bufr_read_all.
use kinds, only: r_kind,i_kind
use constants, only: zero,half,one,two,zero_quad,one_quad
implicit none
character(4) :: this_staid,isstn
real(r_kind),parameter:: four_thirds = 4.0_r_kind / 3.0_r_kind
real(r_kind),parameter:: r8 = 8.0_r_kind
real(r_kind),parameter:: r360 = 360.0_r_kind
real(r_kind),parameter:: r720 = 720.0_r_kind
real(r_kind),parameter:: r89_5 = 89.5_r_kind
real(r_kind) :: erad,aactual,a43,thistilt,pi,deg2rad,rad2deg
real(r_kind) :: reradius,rearth
real(r_kind) :: selev0,celev0,b,c,ha,epsh,h,celev,selev
real(r_kind) :: thisazimuth,deltiltmax,deltiltmin
real(r_kind) :: deldistmax,deldistmin
real(r_kind) :: this_stalat,this_stalon,this_stahgt,thistime,&
thislat,thislon,thishgt,thisvr,corrected_azimuth,&
thiserr,corrected_tilt,gamma,thisrange,thistiltr,&
this_stalatr
real(r_kind) :: thisazimuthr,rlonloc,rad_per_meter,rlatloc,rlon0
real(r_kind) :: clat0,slat0,rlonglob,rlatglob,clat1,caz0,saz0
real(r_kind) :: cdlon,sdlon,caz1,saz1,delazmmax
real(r_kind) :: ilat,ilon,ihgt
real(r_kind),parameter :: r85 = 85._r_kind
integer(i_kind) :: iret!,i,num_supob_boxes
write(6,*) '*******************************************'
reradius=one/6370.e03_r_kind
rearth=one/reradius
rad_per_meter=one/rearth
erad=rearth
pi = acos(-one)
deg2rad = pi/180.0_r_kind
rad2deg = one/deg2rad
thisazimuth=anaz_rw
thistilt=anel_rw
thisrange=range_rw
thisvr=5.0_r_kind
thiserr=magoberr
thistime=obhourset
! Get station lat/lon/hgt from namelist station ID and radar list file.
open(unit=333,file="radar_list",status="old",action="read")
do
read(333,*,iostat=iret) isstn,ilat,ilon,ihgt
if(iret/=0) exit
if(isstn==sstn) then
this_staid=isstn
this_stalat=ilat
this_stalon=ilon
this_stahgt=ihgt
end if
end do
close(333)
rlon0=deg2rad*this_stalon
this_stalatr=this_stalat*deg2rad
clat0=cos(this_stalatr) ; slat0=sin(this_stalatr)
aactual=erad+this_stahgt
a43=four_thirds*aactual
thistiltr=thistilt*deg2rad
selev0=sin(thistiltr)
celev0=cos(thistiltr)
b=thisrange*(thisrange+two*aactual*selev0)
c=sqrt(aactual*aactual+b)
ha=b/(aactual+c)
epsh=(thisrange*thisrange-ha*ha)/(r8*aactual)
h=ha-epsh
thishgt=this_stahgt+h
! Get corrected tilt angle
celev=celev0
selev=selev0
if(thisrange>=one) then
celev=a43*celev0/(a43+h)
selev=(thisrange*thisrange+h*h+two*a43*h)/(two*thisrange*(a43+h))
end if
corrected_tilt=atan2(selev,celev)*rad2deg
deltiltmax=max(corrected_tilt-thistilt,deltiltmax)
deltiltmin=min(corrected_tilt-thistilt,deltiltmin)
gamma=half*thisrange*(celev0+celev)
deldistmax=max(gamma-thisrange,deldistmax)
deldistmin=min(gamma-thisrange,deldistmin)
! Get earth lat lon of superob
thisazimuthr=thisazimuth*deg2rad
rlonloc=rad_per_meter*gamma*cos(thisazimuthr)
rlatloc=rad_per_meter*gamma*sin(thisazimuthr)
call invtllv(rlonloc,rlatloc,rlon0,clat0,slat0,rlonglob,rlatglob)
thislat=rlatglob*rad2deg
thislon=rlonglob*rad2deg
! Get corrected azimuth
clat1=cos(rlatglob)
caz0=cos(thisazimuthr)
saz0=sin(thisazimuthr)
cdlon=cos(rlonglob-rlon0)
sdlon=sin(rlonglob-rlon0)
caz1=clat0*caz0/clat1
saz1=saz0*cdlon-caz0*sdlon*slat0
corrected_azimuth=atan2(saz1,caz1)*rad2deg
delazmmax=max(min(abs(corrected_azimuth-thisazimuth-r720),&
abs(corrected_azimuth-thisazimuth-r360),&
abs(corrected_azimuth-thisazimuth ),&
abs(corrected_azimuth-thisazimuth+r360),&
abs(corrected_azimuth-thisazimuth+r720)),delazmmax)
! Ensure elevation angle is exactly 90 or 0 degress for these one ob tests.
corrected_azimuth=anaz_rw
corrected_tilt=anel_rw
thishgt=range_rw
write(6,*) 'Single radial wind observation.'
write(6,*) '*******************************************'
write(6,*) 'this_staid = ',this_staid
write(6,*) 'this_stalat = ',this_stalat
write(6,*) 'this_stalon = ',this_stalon
write(6,*) 'this_stahgt = ',this_stahgt
write(6,*) 'thistime = ',thistime
write(6,*) 'thislat = ',thislat
write(6,*) 'thislon = ',thislon
write(6,*) 'thishgt = ',thishgt
write(6,*) 'maginnov = ',maginnov
!write(6,*) 'thisvr = ',thisvr ! This quantity is not used for SOT.
write(6,*) 'corrected_az= ',corrected_azimuth
write(6,*) 'corrected_tl= ',corrected_tilt
write(6,*) 'thiserr = ',thiserr
open(10,file='radar_supobs_from_level2',form='unformatted',iostat=iret)
write(10) this_staid,this_stalat,this_stalon,this_stahgt, &
thistime,thislat,thislon,thishgt,thisvr,corrected_azimuth,&
thiserr,corrected_tilt
close(10)
write(6,*) '*******************************************'
write(6,*) '*******************************************'
return
end subroutine oneobmakerwsupob
subroutine oneobo3lv
!$$$ subprogram documentation block
! . . . .
! subprogram: oneobmls
!
! abstract: create ozone level text file for single ob experiment
!
! program history log:
! 2007-09-11 Sienkiewicz - extend to create MLS text file on request
!
! input argument list:
!
! output argument list:
!
! attributes:
! language: f90
! machine: ?
!
!$$$
use constants, only: zero, one
implicit none
integer(i_kind) lumk ! output unit
integer(i_kind) ilev, isnd
integer(i_kind) ildat(8),jldat(8) ! "local" date/time
real(r_kind) rsec,rlnc(5)
real(r_kind) ppmv
2 format(i5,4i3,f6.2,i7,i5,f10.4,f11.4,e16.7,i7,i5,g16.7,g15.7,f6.3)
lumk = 22
ilev = 1 ! ilev > 24 is passive
isnd = 1
ppmv = one ! dummy value
! obdattim=2000010100
rlnc = zero
rlnc(2) = obhourset
ildat(1) = obdattim / 1000000 ! year
ildat(2) = mod(obdattim,1000000)/10000 ! month
ildat(3) = mod(obdattim,10000)/100 ! day
ildat(4) = 0
ildat(5) = mod(obdattim,100) ! hour
ildat(6:8) = 0 ! (no minute/sec in obdattim)
call w3movdat(rlnc,ildat,jldat)
rsec = jldat(7)+jldat(8)*1.e-3_r_kind
! open data file for output. for oneobtype gsimain sets the dfile(1)
! to be prepqc
open(unit=lumk,file='prepqc',form='formatted')
write(lumk,2) jldat(1),jldat(2),jldat(3),jldat(5),jldat(6),rsec,isnd,&
ilev,oblat,oblon,ppmv,isnd,isnd,magoberr,obpres,magoberr
close(lumk)
return
end subroutine oneobo3lv
SUBROUTINE invtllv(ALM,APH,TLMO,CTPH0,STPH0,TLM,TPH)
! borrowed from read_radar
use kinds, only: r_kind
implicit none
real(r_kind),intent(in ) :: alm,aph,tlmo,ctph0,stph0
real(r_kind),intent( out) :: tlm,tph
real(r_kind):: relm,srlm,crlm,sph,cph,cc,anum,denom
RELM=ALM
SRLM=SIN(RELM)
CRLM=COS(RELM)
SPH=SIN(APH)
CPH=COS(APH)
CC=CPH*CRLM
ANUM=CPH*SRLM
DENOM=CTPH0*CC-STPH0*SPH
TLM=tlmo+ATAN2(ANUM,DENOM)
TPH=ASIN(CTPH0*SPH+STPH0*CC)
END SUBROUTINE invtllv
end module oneobmod