!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE INSERT_SAO1 ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
SUBROUTINE insert_sao1(nx,ny,nz,dx,dy,xs,ys,zs,topo,t_sfc_k, & 1,24
cldcv,wtcldcv,t_modelfg,cf_modelfg, &
nobsng,stn,obstype,xsng,ysng,ista_snd, &
obstime,latsta,lonsta,elevsta, &
kcloud,store_amt,store_hgt, &
l_stn_clouds,n_cld_snd,cld_snd,wt_snd, &
i_snd,j_snd, &
istatus)
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Ingest the SAO cloud coverage observations.
!
!-----------------------------------------------------------------------
!
! AUTHOR: (Jian Zhang)
! 02/96. Based on the LAPS cloud analysis code by Steve Albers,
! 1995.
!
! MODIFICATION HISTORY:
!
! 02/06/96 J. Zhang
! Modified for ADAS grid. Added documents.
! 03/14/97 J. Zhang
! Clean up the code and implemented for the official
! arps4.2.4 version
! 04/15/97 J. Zhang
! Added error message output for the case when the
! actual # of cloud soundings (N_CLD_SND) exceeds the
! maximum allowed # (MAX_CLD_SND) defined in
! adascld24.inc.
! 08/06/97 J. Zhang
! Change adascld24.inc to adascld25.inc
! 09/09/97 J. Zhang
! Assign a weight to a SAO cloud sounding base on
! the cloud coverage amount.
! 09/11/97 J. Zhang
! Change adascld25.inc to adascld26.inc
! 02/17/98 Keith Brewster
! Changed logic to skip cloud layers with below zero
! height -- missing, in other words.
! 05/05/98 J. Zhang
! Abandoned the cloud grid, using the ARPS grid instead.
! 05/01/98 Keith Brewster
! Removed abort for missing data, instead data are
! checked for missing at every step.
!
!-----------------------------------------------------------------------
!
! INCLUDE: (from adas.inc)
!
! mx_sng ! max. possible # of surface obs.
! max_cld_snd ! max. possible # of stations
! with cloud coverage reports.
!
! INPUT:
!
! nx,ny,nz ! Number of ADAS grid pts.in 3 directions.
!
! l_stn_clouds ! Using SAO stations' cloud obs?
!
!c ARPS grid variables.
!
! xs (nx) ! The x-coord. of the physical and
! computational grid. Defined at p-point.
! ys (ny) ! The y-coord. of the physical and
! computational grid. Defined at p-point.
! zs (nx,ny,nz) ! The physical height coordinate defined at
! p-point of the staggered grid.
! topo (nx,ny) ! The height of the terrain (equivalent
!
! First guess fields
!
! cf_modelfg (nx,ny,nz) ! first guess cloud cover field
! t_modelfg (nx,ny,nz) ! first guess temperature field
! t_sfc_k (nx,ny) ! surface temperature field
!
! cldcv (nx,ny,nz) 3D gridded fractional cloud cover analysis.
! (when input, it's the first guess field)
! wtcldcv (nx,ny,nz) weights assigned to gridded fractional
! cloud cover analysis. (when input, it's
! the weights for first guess field)
!
! Single-level (e.g., surface) station variables
!
! stn (mx_sng) ! station name of single-lvel data
! obstype (mx_sng) ! names of sfc sources
! obstime (mx_sng) ! time for the observation.
! latsta (mx_sng) ! latitude of single-level data
! lonsta (mx_sng) ! longitude of single-level data
! elevsta (mx_sng) ! height of single-level data
! xsng (mx_sng) ! x location of single-level data
! ysng (mx_sng) ! y location of single-level data
!
! kcloud (mx_sng) ! number of obs. cloud layers.
! store_amt (mx_sng,5) ! cloud coverage (ea. layer,
! ea. station).
! store_hgt (mx_sng,5) ! height of obs. cloud layers.
!
! OUTPUT:
!
! istatus The flag indicating process status.
!
! n_cld_snd number of cloud soundings created
! cld_snd (max_cld_snd,nz) Cld snding obtained from SAO data.
! wt_snd (max_cld_snd,nz) weights assigned to cloud sounding
! obtained from SAO data.
! i_snd (max_cld_snd) i-location of each cloud sounding
! station in ADAS grid
! j_snd (max_cld_snd) j-location of each cloud sounding
! station in ADAS grid
!
! LOCAL:
!
! name_array (max_cld_snd) station names (first letter) for
! each cloud sounding
! ista_snd (max_cld_snd) index of cloud sounding stations
! cvr_snd (max_cld_snd) column integrated cloud cover
! cloud_top(nx,ny) Analyzed cloud top heights (m ASL).
! cloud_base(nx,ny) Analyzed cloud base heights (m ASL).
! cloud_ceiling(nx,ny) Analyzed cloud ceiling heights (m ASL).
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'adas.inc'
!
!-----------------------------------------------------------------------
!
! INCLUDE: (from adas.inc)
!
! integer mx_sng ! max. possible # of surface obs.
! integer max_cld_snd ! max. possible # of stations
! with cloud coverage reports.
!
! INPUT:
INTEGER :: nx,ny,nz
REAL :: dx,dy ! ADAS grid spacing
LOGICAL :: l_stn_clouds ! Using SAO stns' cloud obs? (or bkgrnd)
REAL :: xs (nx) ! The x-coord. of the physical and
! computational grid. Defined at p-point.
REAL :: ys (ny) ! The y-coord. of the physical and
! computational grid. Defined at p-point.
REAL :: zs (nx,ny,nz) ! The physical height coordinate defined
! at p-point of the staggered grid.
REAL :: topo (nx,ny) ! The height of the terrain (equivalent
!
! First guess fields
!
REAL :: cf_modelfg (nx,ny,nz) ! first guess cloud cover field
REAL :: t_modelfg (nx,ny,nz) ! first guess temperature field
REAL :: t_sfc_k (nx,ny) ! surface temperature field
!
REAL :: cldcv (nx,ny,nz) ! 3D gridded frac cld cv analysis.
! (it's also bkgrnd field when input)
REAL :: wtcldcv (nx,ny,nz) ! wts assgned to cld cvr analysis
! (it's also bkgrnd field when input)
!
! Surface cloud coverage reports
!
INTEGER :: nobsng
CHARACTER (LEN=5) :: stn (mx_sng) ! station name of single-lvel data
CHARACTER (LEN=8) :: obstype (mx_sng)! names of sfc sources
INTEGER :: obstime (mx_sng) ! time for the observation.
REAL :: xsng (mx_sng) ! x location of single-level data
REAL :: ysng (mx_sng) ! y location of single-level data
REAL :: latsta (mx_sng) ! latitude of single-level data
REAL :: lonsta (mx_sng) ! longitude of single-level data
REAL :: elevsta (mx_sng) ! height of single-level data
!
INTEGER :: kcloud (mx_sng) ! number of cloud layers.
CHARACTER (LEN=4) :: store_amt (mx_sng,5) ! cld coverage (ea.lyr, ea. stn).
REAL :: store_hgt (mx_sng,5) ! height of cloud layers.
!
! OUTPUT:
!
INTEGER :: istatus ! Flag for process status
!
INTEGER :: n_cld_snd ! # of cloud snds created
REAL :: cld_snd (max_cld_snd,nz) ! Cld snd obtained from SAO data.
REAL :: wt_snd (max_cld_snd,nz) ! wgt for SAO cld snd
INTEGER :: i_snd (max_cld_snd) ! i-lctn of cld snd stn in ADAS grid
INTEGER :: j_snd (max_cld_snd) ! j-lctn of cld snd stn in ADAS grid
!
! LOCAL:
!
CHARACTER (LEN=1) :: name_array (nx,ny) ! Cld snd stn name (1st letter)
INTEGER :: ista_snd (max_cld_snd) ! index of cld snd stns
REAL :: cvr_snd (max_cld_snd) ! column integrated cloud cover
!
!-----------------------------------------------------------------------
!
! Empirical factors (parameters) for cloud cover analysis
!
!-----------------------------------------------------------------------
!
!c Using SAO stations slightly outside the ARPS domain.
INTEGER :: ix_low,ix_high,iy_low,iy_high
!
!-----------------------------------------------------------------------
!
! Misc local variables
!
!-----------------------------------------------------------------------
!
REAL :: ht_defined ! the reachable height of the SAOs
REAL :: ht_base,ht_top ! base and top heights of a cloud layer
INTEGER :: n_analyzed ! # of stations analyzed in this routine
INTEGER :: i,k,l
REAL :: ri(mx_sng),rj(mx_sng)
! i-,j-lctns of each stn in ADAS grid
INTEGER :: iloc,jloc ! i- and j-index of each sta.
INTEGER :: igrd,jgrd ! i- and j-index of each sta. in grid
REAL :: cover ! cloud fractional cover values
REAL :: cld_thk
INTEGER :: k_ceil
LOGICAL :: l_out_of_bounds,l_dry,cldprt
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
istatus = 0
WRITE(6,*)
IF (nobsng < 1) THEN
PRINT*,'## No SAO data available. Returning from INSERT_SAO...'
istatus = 1
RETURN
END IF
WRITE(6,*) 'Inserting SAO data from ',nobsng,' stations.'
cldprt=(clddiag == 1)
!
ix_low = 1 - i_perimeter
ix_high = nx + i_perimeter
iy_low = 1 - i_perimeter
iy_high = ny + i_perimeter
!
!-----------------------------------------------------------------------
!
! Initialize the cloud sounding array.
!
!-----------------------------------------------------------------------
!
DO i=1,max_cld_snd
DO k=1,nz
cld_snd(i,k) = 0.0
wt_snd(i,k) = 0.0
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Initially assign a unreal number for the sensor's reachable
! height.
!
!-----------------------------------------------------------------------
!
ht_defined = 99999.
n_analyzed = 0
!
!-----------------------------------------------------------------------
!
! Loop through the stations
!
!-----------------------------------------------------------------------
!
DO i=1,nobsng
IF( obstype(i)(1:4) == 'meso') GO TO 125
IF( obstype(i)(1:4) == 'armmn') GO TO 125
IF( obstype(i)(1:4) == 'isws') GO TO 125
IF( obstype(i)(1:4) == 'coagmet') GO TO 125
IF( obstype(i)(1:4) == 'hplains') GO TO 125
IF( obstype(i)(1:4) == 'wpdn') GO TO 125
!
!-----------------------------------------------------------------------
!
! Place station at proper ADAS grid point.
!
!-----------------------------------------------------------------------
!
ri(i) = 1. + (xsng(i) - xs(1))/dx
rj(i) = 1. + (ysng(i) - ys(1))/dy
iloc = nint(ri(i))
jloc = nint(rj(i))
IF( iloc < ix_low .OR. iloc > ix_high &
.OR. jloc < iy_low .OR. jloc > iy_high) THEN
! write(6,*) 'Station is outside domain ',stn(i)
GO TO 125
END IF
IF( kcloud(i) == 0) THEN ! Kick out AMOS stations not
! reporting clouds this time
! write(6,*)' No cloud layers - probably AMOS -',
! : ' goto end of loop ',stn(i),obstype(i)
GO TO 125
END IF
n_analyzed = n_analyzed + 1
n_cld_snd = n_cld_snd + 1
IF(n_cld_snd > max_cld_snd) THEN
PRINT*,'##ERROR: actual # of cldsnd=',n_cld_snd, &
' exceeds max. # allowed =',max_cld_snd
PRINT*,'Please increase MAX_CLD_SND in the .inc file.'
PRINT*,'ABORTING from INSERT_SAO......'
STOP
END IF
ista_snd(n_cld_snd) = i
IF(obstype(i)(5:5) /= ' ' .AND. obstype(i)(4:7) /= 'AMOS') THEN
!
!-----------------------------------------------------------------------
!
! Automated Station (12000' limit)
!
!-----------------------------------------------------------------------
!
ht_defined = elevsta(i) + 12000./3.281
ELSE
ht_defined = 99999.
END IF
IF(cldprt) THEN
WRITE(6,1,ERR=110) stn(i),latsta(i),lonsta(i) &
,iloc,jloc,kcloud(i) &
,obstype(i),ht_defined
1 FORMAT(1X,a5,' lat=',f7.2,' lon=',f7.2,' i=',i3,' j=',i3 &
,' kld=',i1,' obsty=',a8,' h_def=',f8.0)
110 WRITE(6,2,ERR=3) &
(store_amt(i,l),store_hgt(i,l),l=1,kcloud(i))
2 FORMAT(1X,5(a4,f8.0))
3 CONTINUE
END IF
IF( iloc < 1 .OR. iloc > (nx-1) .OR. jloc < 1 .OR. jloc > (ny-1) ) THEN
l_out_of_bounds = .true.
igrd=MIN(MAX(iloc,1),(nx-1))
jgrd=MIN(MAX(jloc,1),(ny-1))
ELSE
l_out_of_bounds = .false.
igrd=iloc
jgrd=jloc
name_array(iloc,jloc )=stn(i)(1:1)
name_array(iloc,MIN(jloc+1,ny)) =stn(i)(1:1)
name_array(iloc,MAX(jloc-1,1 ))=stn(i)(1:1)
END IF
cvr_snd (n_cld_snd) = 0. ! column intg. cloud cover.
!
!-----------------------------------------------------------------------
!
! For each station, loop through each cloud layer observed
! by the SAO.
!
!-----------------------------------------------------------------------
!
DO l = 1,kcloud(i)
cover = 0.
IF (store_hgt(i,l) < 0.0 .AND. store_amt(i,l) /= ' CLR' ) THEN
! ht_base = cld_base(l)
! print*,' stn',i,' level',l,' cld.ht.=',store_hgt(i,l)
! : ,' adj. cld. ht.=',ht_base
! ht_base = elevsta(i) + ht_base
CYCLE
ELSE
ht_base = elevsta(i) + store_hgt(i,l)
END IF
IF( ht_base > ht_defined+1.) THEN
IF( store_amt(i,l) /= ' CLR') THEN ! Clouds
WRITE(6,*)' Error, inconsistent SAO data, cloud base is' &
//' reported to be too high for this sensor type'
WRITE(6,*) ht_base,ht_defined,obstype(i)
WRITE(6,*)' Please check cloud layer heights in the LSO' &
//' file to see that they are compatable with the'
WRITE(6,*)' types of sensors used.'
istatus = 0
RETURN
ELSE ! CLR
WRITE(6,*)' Note, CLR sky cloud base does not' &
//' reflect this sensor type'
WRITE(6,*) ht_base,ht_defined,obstype(i)
END IF ! Clouds
END IF !"ht_base > ht_defined+1"
!
!-----------------------------------------------------------------------
!
! CLOUDS ARE NOW IN MSL
! Fill in clear for entire column for SAO or up to ht_base for AWOS
!
!-----------------------------------------------------------------------
!
IF( store_amt(i,l) == ' CLR') THEN
cover=.01
DO k=1,nz
IF( zs(igrd,jgrd,k) <= ht_base &
.AND. zs(igrd,jgrd,k) <= ht_defined ) THEN
CALL spread2
(cld_snd,wt_snd,i_snd,j_snd,n_cld_snd &
,max_cld_snd,nz,iloc,jloc,k,cover,1.)
END IF
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Assign cloud cover values for entire column for SAO
! or up to ht_base for AWOS. First check if there is dry
! layer or inversion.
!
!-----------------------------------------------------------------------
!
IF( store_amt(i,l) == '-SCT' ) THEN
cover=.15 ! .2
ht_top=ht_base+cld_thk(ht_base-elevsta(i),cover)
DO k=2,nz-2
IF(zs(igrd,jgrd,k) >= ht_base .AND. zs(igrd,jgrd,k) <= ht_top) THEN
!
!-----------------------------------------------------------------------
!
! Check if inversion or dry layer exist.
!
!-----------------------------------------------------------------------
!
CALL modify_sounding(cf_modelfg,t_modelfg,topo &
,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs &
,ht_base,ht_top,0,l_dry)
IF(.NOT. l_dry) THEN
CALL spread2(cld_snd,wt_snd,i_snd,j_snd,n_cld_snd &
,max_cld_snd,nz,iloc,jloc,k,cover,1.)
END IF
ELSE IF(zs(igrd,jgrd,k) < ht_base) THEN
!
!-----------------------------------------------------------------------
!
! Initialize the modify sounding routine
!
!-----------------------------------------------------------------------
!
CALL modify_sounding(cf_modelfg,t_modelfg,topo &
,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs &
,ht_base,ht_top,1,l_dry)
END IF ! ht_base < zs < ht_top
END DO ! k = 1, nz
END IF
IF( store_amt(i,l) == ' SCT') THEN
cover=.25 ! .3
ht_top=ht_base+cld_thk(ht_base-elevsta(i),cover)
DO k=2,nz-2
IF(zs(igrd,jgrd,k) >= ht_base .AND. zs(igrd,jgrd,k) <= ht_top) THEN
!
!-----------------------------------------------------------------------
!
! Check if inversion or dry layer exist.
!
!-----------------------------------------------------------------------
!
CALL modify_sounding(cf_modelfg,t_modelfg,topo &
,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs &
,ht_base,ht_top,0,l_dry)
IF(.NOT. l_dry) THEN
CALL spread2(cld_snd,wt_snd,i_snd,j_snd,n_cld_snd &
,max_cld_snd,nz,iloc,jloc,k,cover,1.)
END IF
ELSE IF(zs(igrd,jgrd,k) < ht_base) THEN
!
!-----------------------------------------------------------------------
!
! Initialize the modify sounding routine
!
!-----------------------------------------------------------------------
!
CALL modify_sounding(cf_modelfg,t_modelfg,topo &
,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs &
,ht_base,ht_top,1,l_dry)
END IF
END DO
END IF
IF( store_amt(i,l) == '-BKN' ) THEN
cover=.4 ! .5
ht_top=ht_base+cld_thk(ht_base-elevsta(i),cover)
DO k=2,nz-2
IF(zs(igrd,jgrd,k) >= ht_base .AND. zs(igrd,jgrd,k) <= ht_top) THEN
!
!-----------------------------------------------------------------------
!
! Check if inversion or dry layer exist.
!
!-----------------------------------------------------------------------
!
CALL modify_sounding(cf_modelfg,t_modelfg,topo &
,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs &
,ht_base,ht_top,0,l_dry)
IF(.NOT. l_dry)THEN
CALL spread2(cld_snd,wt_snd,i_snd,j_snd,n_cld_snd &
,max_cld_snd,nz,iloc,jloc,k,cover,1.)
END IF
ELSE IF(zs(igrd,jgrd,k) < ht_base) THEN
!
!-----------------------------------------------------------------------
!
! Initialize the modify sounding routine
!
!-----------------------------------------------------------------------
!
CALL modify_sounding(cf_modelfg,t_modelfg,topo &
,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs &
,ht_base,ht_top,1,l_dry)
END IF
END DO
END IF
IF( store_amt(i,l) == ' BKN' ) THEN
cover=.7
ht_top = ht_base + cld_thk(ht_base-elevsta(i),cover) ! 1500.
DO k=2,nz-2
IF(zs(igrd,jgrd,k) >= ht_base .AND. zs(igrd,jgrd,k) <= ht_top) THEN
!
!-----------------------------------------------------------------------
!
! Check if inversion or dry layer exist.
!
!-----------------------------------------------------------------------
!
CALL modify_sounding(cf_modelfg,t_modelfg,topo &
,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs &
,ht_base,ht_top,0,l_dry)
IF(.NOT. l_dry) THEN
CALL spread2(cld_snd,wt_snd,i_snd,j_snd,n_cld_snd &
,max_cld_snd,nz,iloc,jloc,k,cover,1.)
END IF
ELSE IF(zs(igrd,jgrd,k) < ht_base) THEN
!
!-----------------------------------------------------------------------
!
! Initialize the modify sounding routine
!
!-----------------------------------------------------------------------
!
CALL modify_sounding(cf_modelfg,t_modelfg,topo &
,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs &
,ht_base,ht_top,1,l_dry)
END IF
END DO
END IF
IF( store_amt(i,l) == '-OVC') THEN
cover=.6 ! .9
ht_top=ht_base+cld_thk(ht_base-elevsta(i),cover)
DO k=1,nz
IF(zs(igrd,jgrd,k) >= ht_base .AND. zs(igrd,jgrd,k) <= ht_top) THEN
!
!-----------------------------------------------------------------------
!
! Check if inversion or dry layer exist.
!
!-----------------------------------------------------------------------
!
CALL modify_sounding(cf_modelfg,t_modelfg,topo &
,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs &
,ht_base,ht_top,0,l_dry)
IF(.NOT. l_dry) THEN
CALL spread2(cld_snd,wt_snd,i_snd,j_snd,n_cld_snd &
,max_cld_snd,nz,iloc,jloc,k,cover,1.)
END IF
ELSE IF(zs(igrd,jgrd,k) < ht_base) THEN
!
!-----------------------------------------------------------------------
!
! Initialize the modify sounding routine
!
!-----------------------------------------------------------------------
!
CALL modify_sounding(cf_modelfg,t_modelfg,topo &
,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs &
,ht_base,ht_top,1,l_dry)
END IF
END DO
END IF
IF( store_amt(i,l) == ' OVC') THEN
cover=1.01
ht_top = ht_base + cld_thk(ht_base-elevsta(i),cover) ! 1500.
DO k=2,nz-2
IF(zs(igrd,jgrd,k) >= ht_base .AND. zs(igrd,jgrd,k) <= ht_top) THEN
!
!-----------------------------------------------------------------------
!
! Check if inversion or dry layer exist.
!
!-----------------------------------------------------------------------
!
CALL modify_sounding(cf_modelfg,t_modelfg,topo &
,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs &
,ht_base,ht_top,0,l_dry)
IF(.NOT. l_dry)THEN
CALL spread2(cld_snd,wt_snd,i_snd,j_snd,n_cld_snd &
,max_cld_snd,nz,iloc,jloc,k,cover,1.)
END IF
ELSE IF(zs(igrd,jgrd,k) < ht_base) THEN
!
!-----------------------------------------------------------------------
!
! Initialize the modify sounding routine
!
!-----------------------------------------------------------------------
!
CALL modify_sounding(cf_modelfg,t_modelfg,topo &
,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs &
,ht_base,ht_top,1,l_dry)
END IF
END DO
END IF
IF( store_amt(i,l) == ' X ') THEN
cover=1.01
ht_top = ht_base + cld_thk(ht_base-elevsta(i),cover) ! 1500.
DO k=2,nz-2
IF(zs(igrd,jgrd,k) >= ht_base .AND. zs(igrd,jgrd,k) <= ht_top) THEN
!
!-----------------------------------------------------------------------
!
! Check if inversion or dry layer exist.
!
!-----------------------------------------------------------------------
!
CALL modify_sounding(cf_modelfg,t_modelfg,topo &
,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs &
,ht_base,ht_top,0,l_dry)
IF(.NOT. l_dry) THEN
CALL spread2(cld_snd,wt_snd,i_snd,j_snd,n_cld_snd &
,max_cld_snd,nz,iloc,jloc,k,cover,1.)
END IF
ELSE IF(zs(igrd,jgrd,k) < ht_base) THEN
!
!-----------------------------------------------------------------------
!
! Initialize the modify sounding routine
!
!-----------------------------------------------------------------------
!
CALL modify_sounding(cf_modelfg,t_modelfg,topo &
,t_sfc_k,igrd,jgrd,k,nx,ny,nz,zs &
,ht_base,ht_top,1,l_dry)
END IF
END DO
END IF
cvr_snd(n_cld_snd) = 1. - ((1. - cvr_snd(n_cld_snd)) * cover)
END DO
!
!-----------------------------------------------------------------------
!
! Locate the highest ceiling
!
!-----------------------------------------------------------------------
k_ceil = nz-2
IF( l_stn_clouds) THEN
DO k=nz-2,2,-1
IF( wt_snd(n_cld_snd,k) == 1.00 .AND. cld_snd(n_cld_snd,k) > 0.5) THEN
k_ceil = k
GO TO 1001
END IF
END DO
ELSE
DO k=nz-2,2,-1
IF( wtcldcv(igrd,jgrd,k) == 1.00 .AND. cldcv(igrd,jgrd,k) > 0.5 ) THEN
k_ceil = k
GO TO 1001
END IF
END DO
END IF
!
!-----------------------------------------------------------------------
!
! Fill in other clear layers outside of clouds, below the ceiling,
! and within defined height range of sensor.
!
!-----------------------------------------------------------------------
!
1001 cover = .01
IF( l_stn_clouds) THEN
DO k=2,k_ceil
IF( wt_snd(n_cld_snd,k) /= 1.00 &
.AND. zs(igrd,jgrd,k) <= ht_defined ) THEN
CALL spread2(cld_snd,wt_snd,i_snd,j_snd,n_cld_snd &
,max_cld_snd,nz,iloc,jloc,k,cover,1.)
END IF
END DO
ELSE
DO k=2,k_ceil
IF( wtcldcv(igrd,jgrd,k) /= 1.00 &
.AND. zs(igrd,jgrd,k) <= ht_defined ) THEN
CALL spread2(cld_snd,wt_snd,i_snd,j_snd,n_cld_snd &
,max_cld_snd,nz,iloc,jloc,k,cover,1.)
END IF
END DO
END IF
125 CONTINUE
END DO ! I=1,NOBSNG
!
!-----------------------------------------------------------------------
!
WRITE(6,*)' Num stations analyzed/cloud soundings = ' &
,n_analyzed,n_cld_snd
!
!-----------------------------------------------------------------------
!
istatus = 1
! 999 CONTINUE
RETURN
END SUBROUTINE insert_sao1
!
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE MODIFY_SOUNDING ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
SUBROUTINE modify_sounding (cf_modelfg,t_modelfg,topo, & 14
t_sfc_k,i_in,j_in,k,ni,nj,nk,zs, &
ht_base,ht_top,init,l_dry)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Check if there is dry layer or inversion (If there is one,
! ( the observed cloud layer will be cleared out).
!
!-----------------------------------------------------------------------
!
! AUTHOR: (Jian Zhang)
! 03/96 Based on the LAPS cloud analysis code, 1995
!
! MODIFICATION HISTORY:
!
! 02/06/96 J. Zhang
! Modified for ADAS grid. Added documents.
!
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
! INPUT:
INTEGER :: init ! flag value indicating if the layer is
! within the cloud.
INTEGER :: ni,nj ! horizontal grid domain size
INTEGER :: nk ! # of cloud analysis levels
INTEGER :: i_in,j_in ! i- and j- location in ADAS grid
INTEGER :: k ! the cloud height level
REAL :: cf_modelfg (ni,nj,nk) ! first guess cloud cover field
REAL :: t_modelfg (ni,nj,nk) ! first guess temperature field
REAL :: t_sfc_k (ni,nj) ! surface temperature field
REAL :: topo(ni,nj) ! height of the terrain
!
REAL :: zs(ni,nj,nk) ! cloud analysis heights
REAL :: ht_base,ht_top
!
!-----------------------------------------------------------------------
!
!c OUTPUT:
LOGICAL :: l_dry ! if there is inversion or dry layer?
!
!-----------------------------------------------------------------------
!
!c LOCAL:
INTEGER :: i,j
REAL :: cf_model_base,t_model_base,t_subcloud
REAL :: t_dry_adiabat,t_inversion_strength
LOGICAL :: l_wait_for_base,l_cf,l_inversion
SAVE l_wait_for_base,cf_model_base,t_model_base, &
l_inversion,t_subcloud
!-----------------------------------------------------------------------
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
l_dry = .false.
l_cf = .false.
!
!-----------------------------------------------------------------------
!
! Find ARPS grid point nearest the SAO if it is out of bounds
!
!-----------------------------------------------------------------------
!
i = MAX(MIN(i_in,ni-1),1)
j = MAX(MIN(j_in,nj-1),1)
IF (init == 1) THEN ! (below base )
! reset to wait for the beginning of the layer
l_wait_for_base = .true.
l_inversion = .false.
t_subcloud = t_modelfg(i,j,k)
RETURN
ELSE ! init = 0 (inside cloud layer)
! For k=1 and init=0, l_wait_base=.F., and there will be NO
! t_subcloud available.
IF(l_wait_for_base .OR. k <= 2) THEN
!
!-----------------------------------------------------------------------
!
! Set reference (just within cloud base)
!
!-----------------------------------------------------------------------
!
l_wait_for_base = .false.
cf_model_base = cf_modelfg(i,j,k)
t_model_base = t_modelfg(i,j,k)
! write(6,21) t_subcloud
!21 format( /' model T_subcloud = ',f7.2)
! write(6,1)
!1 format(1x,21x,' cf_fg',' t_fg ',' dt_inv'
! : ,' lcf',' linv',' i j k ',' cldht')
END IF
!
!-----------------------------------------------------------------------
!
! Determine if the layer is dry or it has inversion.
! (in either case, the cloud will be cleared out)
!
!-----------------------------------------------------------------------
!
IF(.true.) THEN ! Set inversion strength flag
!
t_dry_adiabat = t_sfc_k(i,j) &
-.0098 * (zs(i,j,k) - topo(i,j))
t_inversion_strength = t_modelfg(i,j,k) - t_dry_adiabat
IF( ( (t_modelfg(i,j,k) > t_model_base) &
.OR. &
(k >= 2 .AND. t_modelfg(i,j,k) > t_subcloud) ) &
.AND. &
(t_modelfg(i,j,k) > 283.15) & ! temp check
.AND. &
(t_inversion_strength > 4.) & ! delta theta chk
) THEN
!
l_inversion = .true. ! Inversion exists
! write(6,2) cf_modelfg(i,j,k),t_modelfg(i,j,k)
! : ,t_inversion_strength,l_cf,l_inversion
! : ,i,j,k,nint(zs(i,j,k))
!
!2 format(' Inversion detected = ',3f7.2,2l4,3i4,i6)
ELSE IF (cf_modelfg(i,j,k) < cf_model_base - 0.3 & ! cf search
.AND. zs(i,j,k) - ht_base >= 500.) THEN
!
l_cf = .true. ! Dry layer exists
! write(6,3) cf_modelfg(i,j,k),t_modelfg(i,j,k)
! : ,t_inversion_strength,l_cf,l_inversion
! : ,i,j,k,nint(zs(i,j,k))
!
!3 format(' Dry layer detected = ',3f7.2,2l4,3i4,i6)
ELSE ! neither dry nor inversion
!
!-----------------------------------------------------------------------
!
! If l_inversion and l_cf are not reset to .false. here, the WHOLE
! cloud layer above a SINGLE level will be cleared out IF that one
! level has inversion.
! This is effective through the "save" statement at the beginning
! of this subroutine.
!
!-----------------------------------------------------------------------
!
! write(6,4) cf_modelfg(i,j,k),t_modelfg(i,j,k)
! : ,t_inversion_strength,l_cf,l_inversion
! : ,i,j,k,nint(zs(i,j,k))
!
!4 format(' model RH/T in cloud =',3f7.2,2l4,3i4,i6)
END IF ! inversion check
IF( l_cf.OR.l_inversion ) THEN
l_dry = .true.
END IF
END IF ! .true. for dry-inversion check.
END IF ! INIT=1
!
!-----------------------------------------------------------------------
!
RETURN
END SUBROUTINE modify_sounding
!
!##################################################################
!##################################################################
!###### ######
!###### FUNCTION CLD_THK ######
!###### ######
!##################################################################
!##################################################################
!
FUNCTION cld_thk (ht_base,cover)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Obtain thickness of a cloud layer. The thickness is simply
! set to 1 km for cloud below 7 km and 1.5km for those above
! 7 km.
!
!-----------------------------------------------------------------------
!
! AUTHOR: (Jian Zhang)
! 03/96
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
! ht_base ! the height of cloud base
! cover ! the cloud amount of the layer
!
! OUTPUT:
! cld_thk ! the thickness of the cloud layer
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
REAL :: ht_base,cover,cld_thk
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
IF(ht_base > 7000.)THEN
cld_thk = 1500.
ELSE IF(ht_base > 1000.) THEN
cld_thk = 1000.0+(ht_base-1000.)/12.
ELSE
cld_thk = MIN(1000.0,ht_base)
END IF
IF(cover > 0.01 .AND. cover < 0.5)THEN
cld_thk = MIN(1000.,cld_thk)
END IF
RETURN
END FUNCTION cld_thk
!
!##################################################################
!##################################################################
!###### ######
!###### FUNCTION CLD_BASE ######
!###### ######
!##################################################################
!##################################################################
!
FUNCTION cld_base (k)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Obtain the base height of a cloud layer.
! The base height is currently an artificial function of
! the reported cloud layer index.
!
!-----------------------------------------------------------------------
!
! AUTHOR: (Jian Zhang)
! 04/17/96
!
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
! k ! the index of the reported cloud layer
!
! OUTPUT:
! cld_base ! the base height of the cloud layer
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
REAL :: cld_base
INTEGER :: k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
IF(k == 1) cld_base = 500.
IF(k == 2) cld_base = 2000.
IF(k == 3) cld_base = 6000.
IF(k == 4) cld_base = 8000.
IF(k == 5) cld_base = 10000.
RETURN
END FUNCTION cld_base
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE INSERT_RADAR ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE insert_radar(nx,ny,nz,clddiag,topo,zs,temp_3d,z_lcl & 1
,ref_base1,ref_base2,hgt_thresh_ref &
,grid_ref,cldcv &
,cloud_base,cloud_base_buf,l_unresolved &
,istatus)
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Insert radar data into cloud grid
!
!-----------------------------------------------------------------------
!
! AUTHOR: (Jian Zhang)
! 05/1996
!
! MODIFICATION HISTORY:
!
! 05/08/96 J. Zhang
! Modified for the ARPSDAS format. Added full
! documentation.
! 07/26/96 J. Zhang
! Added quality check to avoid out-bounds cloud.
! 03/27/97 J. Zhang
! Updated for the official version of arps4.2.4
! 03/28/97 J. Zhang
! Using the lifting condensation level as the cloud base
! when there is no SAO cloud base.
! 09/01/97 J. Zhang
! Document cleanup for ADASCLD version 25.
! Using the lifting condensation level as the cloud base
! when the radar echo top is below the analyzed SAO cloud
! base
! 09/10/97 J. Zhang
! Put lcl in the input argument list.
! 04/20/98 J. Zhang
! Based on the arps4.3.3 version, abandon the cloud grid.
! Using arps grid.
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
IMPLICIT NONE
!-----------------------------------------------------------------------
!
INTEGER :: nx,ny,nz
! INPUT/OUTPUT:
REAL :: cldcv(nx,ny,nz)
! INPUT:
INTEGER :: clddiag
REAL :: temp_3d(nx,ny,nz)
REAL :: zs(nx,ny,nz)
REAL :: topo(nx,ny)
REAL :: z_lcl(nx,ny) ! lifting condensation level
REAL :: grid_ref(nx,ny,nz)
REAL :: ref_base1 ! "significant" radar echo at lower levels
REAL :: ref_base2 ! "significant" radar echo at upper levels
REAL :: hgt_thresh_ref ! height criteria for "significant" radar
! echo thresholds
! OUTPUT:
INTEGER :: istatus
REAL :: cloud_base(nx,ny)
REAL :: cloud_base_buf(nx,ny) ! Lowest SAO/IR base within search radius
!
! LOCAL:
REAL :: radar_cover
PARAMETER(radar_cover=1.0)
REAL :: thresh_cvr ! lower radar echo threshold for cloud filling
PARAMETER (thresh_cvr = 0.2)
LOGICAL :: l_below_base
LOGICAL :: l_unresolved(nx,ny)
!
!-----------------------------------------------------------------------
!
! Misc local variables
!
!-----------------------------------------------------------------------
!
LOGICAL :: cldprt
REAL :: unlimited_ceiling,ref_thresh
REAL :: zs_1d(nz)
INTEGER :: i,j,k,kp1
INTEGER :: icount_below,isearch_base,insert_count_tot
INTEGER :: icount_radar_lvl,insert_count_lvl
INTEGER :: i_out_bnd,i_lowecho
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
WRITE(6,*) 'Inserting radar reflectivity '
cldprt=(clddiag == 1)
!
!-----------------------------------------------------------------------
!
! Calculate Cloud Bases
!
!-----------------------------------------------------------------------
!
unlimited_ceiling = 200000.
DO j = 1,ny-1
DO i = 1,nx-1
cloud_base(i,j) = unlimited_ceiling
!
DO k = nz-1,1,-1
IF(cldcv(i,j,k) < thresh_cvr .AND. cldcv(i,j,k+1) >= thresh_cvr) THEN
cloud_base(i,j) = 0.5 * (zs(i,j,k) + zs(i,j,k+1))
END IF
END DO ! k
!
IF (cloud_base(i,j) > z_lcl(i,j)) THEN
cloud_base(i,j) = z_lcl(i,j) ! using lcl
END IF
cloud_base_buf(i,j) = cloud_base(i,j)
l_unresolved(i,j) = .false.
END DO ! i
END DO ! j
icount_below = 0 ! tot # of rdr echo pts below cloud base
isearch_base = 0 ! # of neighb cloud bases being succes. found
insert_count_tot = 0 ! tot # of data pts inserted the radar_cover
!
!-----------------------------------------------------------------------
!
! Essentially, this go downward to detect radar tops in time
! to search for a new cloud base
!
!-----------------------------------------------------------------------
!
i_out_bnd=0
i_lowecho=0
DO i = 1,nx-1
DO j = 1,ny-1
DO k=1,nz
zs_1d(k) = zs(i,j,k)
END DO
icount_radar_lvl = 0 ! # of lvls with refl > ref_thresh
insert_count_lvl = 0 ! # of cloud lvls inserted radar_cover
DO k = nz-1,1,-1
kp1 = MIN(k+1,nz)
!
!-----------------------------------------------------------------------
!
! Define the "significant" reflectivity value based on height (to
! eliminate ground clutters and/or other non-weather radar echoes).
!
!-----------------------------------------------------------------------
!
IF ((zs(i,j,k)-topo(i,j)) <= hgt_thresh_ref) THEN
ref_thresh = ref_base1
ELSE
ref_thresh = ref_base2
END IF
IF(grid_ref(i,j,k) > ref_thresh) THEN
icount_radar_lvl = icount_radar_lvl + 1
l_below_base = .false.
!
!-----------------------------------------------------------------------
!
! Test if we are at echo top
!
!-----------------------------------------------------------------------
!
IF(k == nz-1 .OR. grid_ref(i,j,kp1) < ref_thresh) THEN
!
!-----------------------------------------------------------------------
!
! Test if we are below the cloud base.
!
!-----------------------------------------------------------------------
!
IF(zs(i,j,k) < cloud_base_buf(i,j)) THEN
!
!-----------------------------------------------------------------------
!
! Radar Echo Top is below the analyzed cloud base.
! Using the lifting condensation level as the modified cloud base
! if it is lower than the current buffer value.
!
!-----------------------------------------------------------------------
!
cloud_base_buf(i,j)=MIN(z_lcl(i,j),cloud_base_buf(i,j))
IF(cloud_base_buf(i,j) < zs(i,j,k)) THEN
isearch_base = isearch_base + 1
IF(isearch_base < 50) THEN ! limit log output
WRITE(6,71) i,j,k,zs(i,j,k),cloud_base(i,j) &
,cloud_base_buf(i,j)
71 FORMAT(' Rdr Top < Bse*gp=',3I3,' zs=',f7.0 &
& ,' cld_b=',f7.0,'lcl=',f7.0,' Resolved')
END IF
ELSE ! Probably Unresolved base
WRITE(6,72) i,j,k,zs(i,j,k),cloud_base(i,j) &
,cloud_base_buf(i,j)
72 FORMAT(1X,' **** Prob Unresolved ****'/ &
& 1X,'Rdr Top < Bse*gp=',3I3,' zs=',f7.0 &
& ,' cld_b=',f7.0,' lcl=',f7.0)
IF(cloud_base_buf(i,j) == unlimited_ceiling) THEN
l_unresolved(i,j) = .true.
PRINT*,'Error, no LCL found for grid,',i,j, &
' aborting from INSERTRAD...'
STOP
END IF ! cloud_base(i,j).eq.unlimited_ceiling
GO TO 600
END IF ! cloud_base_buf(i,j) .lt. zs(i,j,k)
END IF ! Blw Cld Base: zs.lt.cloud_base_buf(i,j)
END IF ! At Echo Top: ref(k)>thr & (k.eq.nz .or. ref(kp1)<thr
!
!-----------------------------------------------------------------------
!
! Loop through range of cloud grid levels for this model level
!
!-----------------------------------------------------------------------
!
IF(zs(i,j,k) > cloud_base_buf(i,j)) THEN
! Insert radar if we are above cloud base
cldcv(i,j,k) = radar_cover
insert_count_lvl = insert_count_lvl + 1
insert_count_tot = insert_count_tot + 1
ELSE ! Radar Echo below cloud base
l_below_base = .true.
END IF
IF(l_below_base) THEN
icount_below = icount_below + 1
IF(icount_below <= 50) THEN
WRITE(6,81) i,j,k,zs(i,j,k) &
,cloud_base_buf(i,j)
81 FORMAT(1X,'Rdr < Bse* g.p.:',3I3 &
,' zs(i,j,k)=',f7.0,' cld_base=',f7.0)
END IF
GO TO 600
END IF ! l_below_base =.true.
ELSE ! grid_ref(i,j,k) <= ref_thresh
IF(cldcv(i,j,k) > 0.4.AND.i_lowecho <= 25) THEN
i_lowecho=i_lowecho+1
IF(cldprt) WRITE(6,642) grid_ref(i,j,k),cldcv(i,j,k),i,j,k
642 FORMAT(1X,' Rdr echo', f8.1,'<< Cld cvr',f6.2,' at',3I3)
END IF
GO TO 600
END IF ! grid_ref(i,j,k) > ref_thresh ?
IF (cldprt .AND. MOD(insert_count_tot,100) == 0) THEN
WRITE(6,591) k,icount_radar_lvl &
,insert_count_lvl,insert_count_tot
591 FORMAT(1X,'Inserted radar** k=',i3/ &
' n_rdrl=',i5,' n_insl=',i5,' n_instot=',i5)
END IF
600 CONTINUE
END DO ! k
END DO ! j
END DO ! i
WRITE(6,*)' Total cloud grid points modified by radar = ' &
,insert_count_tot
istatus=1
RETURN
END SUBROUTINE insert_radar
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE INSERT_VIS ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE insert_vis (nx,ny,nz,zs,topo,clouds_3d & 1
,albedo,cld_frac_vis_a &
,vis_rad_thcvr,vis_rad_thdbz &
,istat_radar,radar_ref_3d,ref_base,dbz_max_2d &
,r_missing,sfc_sao_buffer,istatus)
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! To process satellite data for clouds and to modify the
! 3d cloud cover array.
! Currently Band 8 (11.2 micron) brightness temps are used with a
! dummy call for the CO2 slicing method.
!
!-----------------------------------------------------------------------
!
! AUTHOR: (Jian Zhang)
! 04/26/96 (Based on LAPS code of 09/95)
!
! MODIFICATION HISTORY:
!
! 04/26/96 (J. Zhang)
! Modified for the ARPSDAS format.
! Added full documentation.
! 09/02/97 (J. Zhang)
! Added if statements to avoid putting in the reflectivity
! threshold values in the data holes.
! 04/20/98 (J. Zhang)
! Abandoned the cloud analysis grid, using the ARPS grid
! instead.
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
INTEGER :: nx,ny,nz
! INPUT/OUTPUT:
REAL :: clouds_3d(nx,ny,nz) ! cloud cover analysis
!
! INPUT:
REAL :: r_missing
REAL :: zs(nx,ny,nz)
REAL :: topo(nx,ny)
!
! INPUT: parameters
REAL :: vis_rad_thcvr ! =0.2, defined in cloud_cv.f
REAL :: vis_rad_thdbz ! =10dbz, defined in cloud_cv.f
REAL :: sfc_sao_buffer ! =800m, defined in cloud_cv.f
!
! INPUT: satellite visible channel data
REAL :: albedo(nx,ny)
REAL :: cld_frac_vis_a(nx,ny)
!
! INPUT: radar data
INTEGER :: istat_radar
REAL :: dbz_max_2d(nx,ny)
REAL :: radar_ref_3d(nx,ny,nz)
REAL :: ref_base
!
! OUTPUT:
INTEGER :: istatus
!
! LOCAL:
INTEGER :: ih_alb(-10:20)
INTEGER :: ih_cf_sat(-10:20) ! Histogram for VIS cf array
INTEGER :: ih_cfin(-10:20) ! Histogram for input cf array
INTEGER :: ih_cfout(-10:20) ! Histogram for output cf array
INTEGER :: ih_cfin_out(-10:20,-10:20) ! Histogram for the comparison
! between input cloud cover
! array and modified cf array.
INTEGER :: ih_cfin_sat(-10:20,-10:20) ! Histogram for the comparison
! between input cloud cover
! array and sat. VIS cf array.
INTEGER :: ih_cmaxin_sat(-10:20,-10:20)
INTEGER :: ih_cmaxout_sat(-10:20,-10:20)
!
!-----------------------------------------------------------------------
!
! Misc local variables
!
!-----------------------------------------------------------------------
!
REAL :: viscorfrc
PARAMETER (viscorfrc=0.29)
INTEGER :: i,j,k,kl
INTEGER :: n_missing_albedo,n_vis_mod
INTEGER :: i_sat,i_in,i_out,i_cmaxin,i_cmaxout
REAL :: diffin,diffout,diffin_sum,diffout_sum &
,diffin_sumsq,diffout_sumsq
REAL :: cld_frac_vis,cld_frac_in,cld_frac_out
REAL :: cushion
REAL :: cmaxin,cmaxout
INTEGER :: iblank_radar,iset_vis
REAL :: r_present
LOGICAL :: l_prt
!
!-----------------------------------------------------------------------
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
istatus=0
!
!-----------------------------------------------------------------------
!
! Initialize all histogram arrays.
!
!-----------------------------------------------------------------------
!
DO i=-10,20
ih_alb(i) = 0
ih_cf_sat(i) = 0
ih_cfin(i) = 0
ih_cfout(i) = 0
DO j=-10,20
ih_cfin_out(i,j) = 0
ih_cfin_sat(i,j) = 0
ih_cmaxin_sat(i,j) = 0
ih_cmaxout_sat(i,j) = 0
END DO
END DO
!
n_missing_albedo = 0
n_vis_mod = 0 ! # of data points modified by VIS data
diffin_sum = 0. ! sum of diff. between column max. INPUT cld
! fraction and the satellite vis cld frac'n.
diffout_sum = 0. ! sum of diff. between column max. OUTPUT cld
! fraction and the satellite vis cld frac'n.
diffin_sumsq = 0.
diffout_sumsq = 0.
!
!-----------------------------------------------------------------------
!
! Horizontal array loop
!
!-----------------------------------------------------------------------
!
l_prt = .false.
DO i = 1,nx-1
DO j = 1,ny-1
!
kl = nint(albedo(i,j)*10.)
kl = MIN(MAX(kl,-10),20)
ih_alb(kl) = ih_alb(kl) + 1
IF(cld_frac_vis_a(i,j) /= r_missing) THEN
cld_frac_vis = cld_frac_vis_a(i,j)
i_sat = nint(cld_frac_vis*10.)
i_sat = MIN(MAX(i_sat,-10),20)
ih_cf_sat(i_sat) = ih_cf_sat(i_sat) + 1 ! # of data points
! in a VIS cf value catagory
!
!-----------------------------------------------------------------------
!
! Make sure satellite cloud fraction is between 0 and 1
!
!-----------------------------------------------------------------------
!
IF(cld_frac_vis < 0.0) cld_frac_vis = 0.0
IF(cld_frac_vis > 1.0) cld_frac_vis = 1.0
i_sat = nint(cld_frac_vis*10.)
cmaxin = 0.
cmaxout = 0.
iblank_radar = 0
iset_vis = 0
DO k = 1,nz-1
! l_out(1)=k.eq.19.and.i.le.17.and.i.ge.2
! : .and.j.le.43.and.j.ge.36
! l_out(3)=k.eq.26.and.i.le.10.and.i.ge.6
! : .and.j.le.26.and.j.ge.30
! l_prt=l_out(1).or.l_out(3)
cld_frac_in = MIN(clouds_3d(i,j,k),1.0) ! save input cf value
!
!-----------------------------------------------------------------------
!
! Modify the cloud field with the vis input - allow .3 vis err?
! The scheme:
! Above the sao buffer layer, the previoud cloud analysis is
! retained if it is within cld_frac_vis + 0.3.
! Below the layer, the cloud analysis must be smaller than
! cld_frac_vis. If not, reset it to cld_frac_vis.
!
!-----------------------------------------------------------------------
!
IF(zs(i,j,k) > topo(i,j)+sfc_sao_buffer) THEN
cushion = 0.3 ! Allow 0.3 error in VIS cld cvr.
ELSE
cushion = 0.0 ! VIS cld cover is the max. allowed.
END IF
IF(cld_frac_vis < viscorfrc .AND. &
((clouds_3d(i,j,k)-cld_frac_vis) > cushion)) THEN
cld_frac_out = cld_frac_vis ! VIS cf value
!
IF(l_prt) THEN
WRITE(6,621) i,j,k,clouds_3d(i,j,k) &
,radar_ref_3d(i,j,k),dbz_max_2d(i,j) &
,cld_frac_vis_a(i,j) &
,cld_frac_out
621 FORMAT(1X,3I3,' oldcld=',f5.2,' ref=',f8.1,' dbz_m=',f8.1 &
,' viscld=',f5.2,' newcld=',f5.2)
END IF
!
!-----------------------------------------------------------------------
!
! Determine if we need to reconcile VIS with radar
!
!-----------------------------------------------------------------------
!
IF(istat_radar == 1 .AND. dbz_max_2d(i,j) /= r_missing &
.AND. dbz_max_2d(i,j) > ref_base) THEN
! Valid radar echo
IF(cld_frac_out < vis_rad_thcvr) THEN
IF(dbz_max_2d(i,j) < vis_rad_thdbz) THEN
! Blank out Radar, Normal VIS Clearing
iblank_radar = 1
ELSE ! dbz_max_2d(i,j) >= vis_rad_thdbz
cld_frac_out = vis_rad_thcvr ! use radar cf
iset_vis = 1
END IF
END IF
END IF
IF(cld_frac_in - cld_frac_out > .01) THEN
n_vis_mod = n_vis_mod + 1
END IF
clouds_3d(i,j,k) = cld_frac_out ! Modify the output
ELSE ! clouds_3d - cld_frac_vis .le. cushion
cld_frac_out = cld_frac_in ! previous cloud analysis
END IF ! clouds_3d - cld_frac_vis .gt. cushion
!
!-----------------------------------------------------------------------
!
! Update Histograms
!
!-----------------------------------------------------------------------
!
i_in = nint(cld_frac_in*10.)
ih_cfin(i_in) = ih_cfin(i_in) + 1
i_out = nint(cld_frac_out*10.)
ih_cfout(i_out) = ih_cfout(i_out) + 1
ih_cfin_sat(i_in,i_sat) &
= ih_cfin_sat(i_in,i_sat) + 1
ih_cfin_out(i_in,i_out) &
= ih_cfin_out(i_in,i_out) + 1
cmaxin = MAX(cmaxin,cld_frac_in) !col.max of inp cldcvr
cmaxout = MAX(cmaxout,cld_frac_out) !col.max of outp cldcvr
END DO ! k
!
!-----------------------------------------------------------------------
!
! Reconcile VIS with radar
!
!-----------------------------------------------------------------------
!
IF(iblank_radar == 1) THEN ! NO VIS / WEAK ECHO
!
!-----------------------------------------------------------------------
!
! Blank out radar column for this grid point
!
!-----------------------------------------------------------------------
!
IF(cmaxout <= vis_rad_thcvr) THEN
WRITE(6,1) i,j,cmaxout,dbz_max_2d(i,j),cld_frac_vis
1 FORMAT(' VIS_RDR - Blank out radar: cvr/dbz/vis << ' &
,2I3,f8.2,f8.1,f8.2)
ELSE ! Some of the column remains above the VIS threshold
! We are "saved" by the VIS cushion
! May not show up in comparisons
WRITE(6,2) i,j,cmaxout,dbz_max_2d(i,j),cld_frac_vis
2 FORMAT(' VIS_RDR - Blank out radar: cvr/dbz/vis-s >> ' &
,2I3,f8.2,f8.1,f8.2)
END IF
IF (dbz_max_2d(i,j) > ref_base) dbz_max_2d(i,j) = ref_base
DO kl = 1,nz
IF (radar_ref_3d(i,j,kl) > ref_base) radar_ref_3d(i,j,kl) = ref_base
END DO ! kl
ELSE IF (iset_vis == 1) THEN ! NO VIS / STRONG ECHO
!
!-----------------------------------------------------------------------
!
! Cloud cvr has been reset to threshold value above VIS
!
!-----------------------------------------------------------------------
!
IF(cmaxout <= vis_rad_thcvr) THEN
WRITE(6,3) i,j,cmaxout,dbz_max_2d(i,j),cld_frac_vis &
,vis_rad_thcvr
3 FORMAT(' VIS_RDR - Reset vis: cvr/dbz/vis/thr << ' &
,2I4,f5.2,f8.1,2F5.2)
ELSE ! Some of the column remains above the VIS threshold
! We are "saved" by the VIS cushion
! May not show up in comparisons
! Is resetting the VIS perhaps not necessary?
WRITE(6,4)i,j,cmaxout,dbz_max_2d(i,j),cld_frac_vis &
,vis_rad_thcvr
4 FORMAT(' VIS_RDR - Reset vis: cvr/dbz/vis/thr-s >> ' &
,2I4,f8.2,f8.1,2F8.2)
END IF
END IF ! iblank_radar=1
IF(j == ny/2 .AND. cmaxin-cmaxout > 0.1) THEN
WRITE(6,12) i,j,cmaxin,cmaxout
12 FORMAT(1X,'Vismod',2I3,' colcfi=',f5.2,' colcfo=',f5.2)
END IF
i_cmaxin = nint(cmaxin*10.)
i_cmaxout = nint(cmaxout*10.)
ih_cmaxin_sat(i_cmaxin,i_sat) &
= ih_cmaxin_sat(i_cmaxin,i_sat) + 1
ih_cmaxout_sat(i_cmaxout,i_sat) &
= ih_cmaxout_sat(i_cmaxout,i_sat) + 1
diffin = cmaxin - cld_frac_vis
diffout = cmaxout - cld_frac_vis
diffin_sum = diffin_sum + diffin
diffout_sum = diffout_sum + diffout
diffin_sumsq = diffin_sumsq + diffin**2
diffout_sumsq = diffout_sumsq + diffout**2
ELSE ! cld_frac_vis_a(i,j) = r_missing
n_missing_albedo = n_missing_albedo + 1
END IF ! cld_frac_vis_a(i,j) .ne. r_missing
END DO ! i
END DO ! j
WRITE(6,*)
WRITE(6,*)' N_MISSING_ALBEDO = ',n_missing_albedo
WRITE(6,*)' N_VIS_MOD = ',n_vis_mod
WRITE(6,*)
WRITE(6,*)' HISTOGRAMS'
WRITE(6,*)' I Alb CFsat CFin CFout'
DO i = -5,15
WRITE(6,11)i,ih_alb(i),ih_cf_sat(i),ih_cfin(i),ih_cfout(i)
11 FORMAT(1X,i3,4I7)
END DO
WRITE(6,*)
WRITE(6,*)' Input vs. Satellite Cloud Fraction Histogram'
WRITE(6,*)
WRITE(6,*)' X-axis: input cld frac, Y-axis: Satellite cld frac'
WRITE(6,20) (i,i=0,10)
20 FORMAT(1X,3X,11I6)
DO j = 0,10
WRITE(6,21) j,(ih_cfin_sat(i,j),i=0,10)
21 FORMAT(1X,i3,11I6)
END DO
WRITE(6,*)
WRITE(6,*)' Input vs. Output Cloud Fraction Histogram'
WRITE(6,*)
WRITE(6,*)' X-axis: input cld frac, Y-axis: output cld fraction'
WRITE(6,20) (i,i=0,10)
DO j = 0,10
WRITE(6,21) j,(ih_cfin_out(i,j),i=0,10)
END DO
WRITE(6,*)
WRITE(6,*)' Column Max Input vs. Satellite Fraction Histogram'
WRITE(6,*)
WRITE(6,*)' X-axis: column max.input cf, Y-axis: col.max.sat.cf.'
WRITE(6,20) (i,i=0,10)
DO j = 0,10
WRITE(6,21) j,(ih_cmaxin_sat(i,j),i=0,10)
END DO
WRITE(6,*)
WRITE(6,*)' Column Max Output vs. Satellite Fraction Histogram'
WRITE(6,*)
WRITE(6,*)' X-axis: column max.output cf, Y-axis: col.max.sat.cf.'
WRITE(6,20) (i,i=0,10)
DO j = 0,10
WRITE(6,21) j,(ih_cmaxout_sat(i,j),i=0,10)
END DO
r_present = nx*ny - n_missing_albedo
IF(r_present > 0.) THEN ! write stats
WRITE(6,31)diffin_sum/r_present,SQRT(diffin_sumsq/r_present)
31 FORMAT(' VIS STATS: Mean/RMS input residual = ',2F8.3)
WRITE(6,41) diffout_sum/r_present,SQRT(diffout_sumsq/r_present)
41 FORMAT(' VIS STATS: Mean/RMS output residual = ',2F8.3)
END IF
WRITE(6,*)
istatus = 1
RETURN
END SUBROUTINE insert_vis
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE INSERT_IR ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE insert_ir (nx,ny,nz,rlat,rlon,rland_frac,cvr_snow & 1,5
,topo,zs,z_lcl,temp_3d,p_3d &
,t_sfc_k,t_gnd_k,cldcv_sao &
,solar_alt,solar_ha,solar_dec &
,tb8_k,tb8_cold_k,cldtop_m_tb8,cldtop_m &
,grid_spacing_m,sfc_sao_buffer &
,cldcv,istatus)
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! To process satellite data for clouds and to modify the
! 3d cloud cover array.
! Currently Band 8 (11.2 micron) brightness temps are used.
!
!-----------------------------------------------------------------------
!
! AUTHOR: (Jian Zhang)
! 04/1996 Based on the LAPS cloud analysis code (Steve Albers,
! 09/1995)
!
! MODIFICATION HISTORY:
!
! 04/26/96 J. Zhang
! Modified for the ARPSDAS format.
! Added full documentation.
! 07/19/96 J. Zhang
! Added the quality check part to deal with the tall
! and steep cloud towers.
! 07/24/96 J. Zhang
! Added the quality check for cases when the cloud top
! is higher than the model top boundary
! 11/01/96 J. Zhang
! Added the rawinsonde data for determing the cloud
! top height for the low clouds.
! 03/18/97 J. Zhang
! Cleanup the code and implemented for the official
! arps4.2.4 version
! 08/06/97 J. Zhang
! Change adascld24.inc to adascld25.inc
! 09/09/97 J. Zhang
! Fixed a bug when calling COMPARE_RAD. Further cleanup.
! 09/10/97 J. Zhang
! Lifting condensation level is now an input argument.
! Change adascld25.inc to adascld26.inc
! 04/20/98 J. Zhang
! Abandoned the cloud analysis grid, using the ARPS grid
! instead.
!
!wdt update
! 2002-06-12 G. Bassett
! Redid the way tb8_cold_k & ht_sao_base are computed so that no
! grid points in the area concerned are skipped.
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! Include files:
INCLUDE 'phycst.inc'
INCLUDE 'adas.inc'
!
!-----------------------------------------------------------------------
!
INTEGER :: nx,ny,nz ! ARPS grid size
! INCLUDE: ( from adas.inc)
! real r_missing
!
! INPUT/OUTPUT:
REAL :: cldcv(nx,ny,nz) ! 3D cldcv analx after useing sat.data
!
! INPUT: general
REAL :: grid_spacing_m
!
! INPUT: for model grid geometry
REAL :: topo(nx,ny)
REAL :: rlat(nx,ny),rlon(nx,ny) ! Lat. and lon. of the grid pt.
REAL :: zs(nx,ny,nz) ! Phys. ht (m) at model grid point.
REAL :: z_lcl(nx,ny) ! Lifting condensation level (MSL)
REAL :: temp_3d(nx,ny,nz) ! temp. analysis on model grid
REAL :: p_3d(nx,ny,nz) ! pres. analysis on model grid
!
! INPUT: for model grid sfc characteristic fields
REAL :: rland_frac(nx,ny) ! land coverage fraction at a g.p.
REAL :: cvr_snow(nx,ny) ! snow coverage at a grid point.
REAL :: solar_alt(nx,ny) ! solar altitude angle
REAL :: solar_ha(nx,ny) ! solar hour angle
REAL :: solar_dec
!
! INPUT: for model grid analysis fields
REAL :: t_sfc_k(nx,ny) ! sfc air temp.
!
! INPUT: for SAO
REAL :: sfc_sao_buffer ! No clearing cloud by sat. below
! this ht.(m AGL) (hence letting
! SAOs dominate)
REAL :: cldcv_sao(nx,ny,nz) ! 3D Cloud cover array
!
! OUTPUT:
INTEGER :: istatus
REAL :: t_gnd_k(nx,ny) ! ground skin temp.
REAL :: tb8_k(nx,ny) ! Obs. sate. band8 bright. temp.
REAL :: tb8_cold_k(nx,ny) ! Cold filtered band8 bright. temp.
REAL :: cldtop_m(nx,ny) ! Adj. Sat. cloud top height (m)
REAL :: cldtop_m_tb8(nx,ny) ! Sat. cld top ht (m) from band 8 data
!
! Local: factors and parameters
REAL :: sfc_ir_buffer !No adding cld by sat. below this ht.
! (m AGL) (hence letting SAOs dominate)
PARAMETER (sfc_ir_buffer = 3000.)
REAL :: thk_def !Def. cld thkness inserted by sat.data
PARAMETER (thk_def = 1500.)
REAL :: thr_sao_cvr ! cldcvr thres. in evaluating presence
! of SAO cld layers
PARAMETER (thr_sao_cvr = 0.1)
REAL :: thresh2 ! Thres.for IR cloud detection
PARAMETER (thresh2 = 3.5) ! (Tsfc - T_IR) Were 5K and 15K also
!
! LOCAL: for band 8 brightness temp.
REAL :: tb8_est ! estimated tb8 temp. by using
! cld cvr analysis
! LOCAL: for interpolations
REAL :: zs_1d(nz)
REAL :: t_1d(nz) ! 3D temperature analysis on model grid
REAL :: p_1d(nz) ! pres. analysis on model grid
!
!wdt update
! LOCAL: for Controling search box for SAO analyzed data
! INTEGER :: idelt(3)
! INTEGER :: jdelt(3)
!
! LOCAL: for satellite cloud presence and cloud top height
LOGICAL :: l_tb8,l_cloud_present
REAL :: cldtop_m_old
INTEGER :: nlyr ! for iterative adj. of cld cvr field
!
! FUNCTIONS:
REAL :: t_ground_k
REAL :: temp_to_rad
REAL :: rad_to_temp
REAL :: band8_cover
!
!-----------------------------------------------------------------------
!
! Misc local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k,i_sao,j_sao
INTEGER :: il,ih,jl,jh,ii,jj
INTEGER :: n_no_sao1,n_no_sao2,n_no_sao3
INTEGER :: i_reset_base
!wdt update
REAL :: htbase,htbase_init,ht_cloud
PARAMETER (ht_cloud = 2000.0)
INTEGER :: i_delt,idelt_max,idelt_index,jdelt_index
!wdt update
REAL :: sat_cover,temp_thresh,cldcv_above,cover,cover_new,buffer
REAL :: tmin
!wdt update
!INTEGER*4 nidelt,njdelt
!DATA nidelt/3/,njdelt/3/
!wdt update
REAL :: ht_sao_base(nx,ny) ! automatic array
REAL :: ht_sao_base_min
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
! Print some band 8 brightness temp. samples.
!
!-----------------------------------------------------------------------
!
! open (29,file='cloud_diag',form='unformatted') ! XXX wdt tmp
WRITE(6,*) 'Insert satellite IR data routines called.'
WRITE(6,*)
!
il=MAX(1,nx/2-5)
ih=MIN(nx,nx/2+5)
jl=MAX(1,ny/2-5)
jh=MIN(ny,ny/2+5)
WRITE(6,*) ' Band 8 brightness temp values in the mid-domain:'
WRITE(6,601) (i,i=il,ih)
601 FORMAT(/1X,3X,11(2X,i2,2X))
WRITE(6,603) (j,(tb8_k(i,j),i=il,ih),j=jh,jl,-1)
603 FORMAT(1X,i3,11F6.1)
!
!-----------------------------------------------------------------------
!
! Define a box for searching SAO cloud base
!
!-----------------------------------------------------------------------
!
idelt_max = nint(50000. / grid_spacing_m)
!wdt update
!idelt(1) = -idelt_max
!idelt(2) = 0
!idelt(3) = +idelt_max
!jdelt(1) = 0
!jdelt(2) = +idelt_max
!jdelt(3) = -idelt_max
!wdt update
! ht_sao_base flags: -1E25 => not computed yet, 1E25 => no sao base.
ht_sao_base = -1E25
!
!-----------------------------------------------------------------------
!
! First guess conversion from cloud height grid to ADAS grid
! This has to err slightly on the high side
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! QC check on band 8 (11.2 mm) brightness temps
!
!-----------------------------------------------------------------------
!
DO j = 1,ny
DO i = 1,nx
IF( tb8_k(i,j) < 173. .OR. tb8_k(i,j) > 350.) tb8_k(i,j)=-999.
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Calculate cold filtered temperatures
!
!-----------------------------------------------------------------------
!
! write (29) tb8_k !XXX wdt tmp
!wdt FIXME: this assumes dx=dy
!wdt update
!i_delt = MAX(1,nint(grid_spacing_m/10000.))
i_delt = MAX(1,nint(10000./grid_spacing_m))
DO j = 1,ny
jl = MAX(1 ,j-i_delt)
jh = MIN(ny,j+i_delt)
DO i = 1,nx
il = MAX(1 ,i-i_delt)
ih = MIN(nx,i+i_delt)
tb8_cold_k(i,j) = 1.0E15
!wdt update i_delt => 1
DO jj = jl,jh,1
DO ii = il,ih,1
IF (tb8_k(ii,jj) > 0.) &
tb8_cold_k(i,j) = MIN(tb8_cold_k(i,j),tb8_k(ii,jj))
END DO ! ii
END DO ! jj
IF(tb8_cold_k(i,j) > 350.) tb8_cold_k(i,j) = -999.
END DO ! i
END DO ! j
! write (29) tb8_cold_k !XXX wdt tmp
!
!-----------------------------------------------------------------------
!
! Calculate ground temperature (for now equate to sfc air temp)
!
!-----------------------------------------------------------------------
!
WRITE(6,*) 'Computing ground skin temperature'
DO j = 1,ny-1
DO i = 1,nx-1
IF ( rland_frac(i,j) >= 0.5) THEN
t_gnd_k(i,j) = t_ground_k(t_sfc_k(i,j),solar_alt(i,j) &
,solar_ha(i,j),solar_dec &
,rlat(i,j),cvr_snow(i,j) &
,r_missing,i,j,nx,ny)
ELSE ! water environment
t_gnd_k(i,j) = t_sfc_k(i,j)
END IF
END DO
END DO
! write (29) t_gnd_k !XXX wdt tmp
!
!-----------------------------------------------------------------------
!
! Check number of points thrown out as a function of offset to check
! for IR navigation errors. This is only a diagnostic process,
! does not change anything.
!
!-----------------------------------------------------------------------
!
CALL correlation(t_gnd_k,tb8_k,thresh2,nx,ny)
!
!-----------------------------------------------------------------------
!
! Initializes routine
!
!-----------------------------------------------------------------------
!
n_no_sao1 = 0 ! # of g.p with sat.cld but no SAO cld at the grid pt
n_no_sao2 = 0 ! # of g.p with sat.cld but no SAO cld in the box
n_no_sao3 = 0 ! # of pts with sat.cldtop lower than SAO cld base.
i_reset_base = 1 ! # of pts that sate. ceiling being resetted
!
!-----------------------------------------------------------------------
!
WRITE(6,*) 'Find cloud_top for each grid point.'
DO j=1,ny-1
DO i=1,nx-1
IF(tb8_k(i,j) > 0.) THEN
DO k=1,nz
zs_1d(k) = zs(i,j,k)
t_1d(k) = temp_3d(i,j,k)
p_1d(k) = p_3d(i,j,k)
END DO
tmin=temp_3d(i,j,2)
DO k=2,nz-1
tmin=MIN(tmin,temp_3d(i,j,k))
END DO
!
!-----------------------------------------------------------------------
!
! Calculate cloud top height from Band 8 method
!
!-----------------------------------------------------------------------
!
CALL cloud_top(nx,ny,nz, &
i,j,zs_1d,t_1d,p_1d,z_ref_lcl,z_lcl(i,j), &
topo(i,j),tmin,r_missing, &
tb8_k(i,j),t_gnd_k(i,j),thresh2, &
cldtop_m_tb8(i,j),cldtop_m(i,j), &
l_tb8,l_cloud_present,sat_cover)
!
!-----------------------------------------------------------------------
!
! Modify those levels where the satellite shows warmer than the
! calculated brightness temp from the analysis of SAO/Pireps
!
!-----------------------------------------------------------------------
!
DO k = nz-1,1,-1
IF (cldcv(i,j,k) > 0.04) THEN ! Efficiency test
!
!-----------------------------------------------------------------------
!
! Estimate the brightness temperature using the cloud cover
! field.
!
!-----------------------------------------------------------------------
!
IF (cldcv(i,j,k) /= r_missing) THEN
tb8_est = temp_to_rad(temp_3d(i,j,k)) &
* cldcv(i,j,k) + temp_to_rad(t_gnd_k(i,j)) &
* (1.-cldcv(i,j,k))
tb8_est = rad_to_temp(tb8_est) ! estimated tb8 temp.
! tb8_est_a = temp_3d(i,j,k) * cldcv(i,j,k)
! : + t_gnd_k(i,j) * (1.-cldcv(i,j,k))
ELSE
tb8_est = t_gnd_k(i,j)
END IF
!
!-----------------------------------------------------------------------
!
! Test if clouds detected by SAO/PIREP should have been
! detected by the satellite (if satellite is warmer than analysis)
!
!-----------------------------------------------------------------------
!
IF (tb8_k(i,j) > tb8_est) THEN !seems too much SAO cld
! Don't touch points within buffer of surface
IF (zs(i,j,k)-topo(i,j) > sfc_sao_buffer) THEN
! Does satellite still imply at least some cloud?
IF(t_gnd_k(i,j)-tb8_k(i,j) > thresh2) THEN ! Some cloud
IF(cldcv(i,j,k) > 0.9) THEN
cldcv(i,j,k)=.01 ! Lower top of solid cld
ELSE ! Cover < 0.9, correct it
cldcv(i,j,k) = band8_cover( tb8_k(i,j) &
,t_gnd_k(i,j),temp_3d(i,j,k))
cldcv(i,j,k) = MAX(0.0,MIN(1.0,cldcv(i,j,k)))
END IF
ELSE !tg-tb<8 Band 8 nearly matches grnd, clear it
!
!-----------------------------------------------------------------------
!
! Insure that "black (or grey) stratus" is not present
!
!-----------------------------------------------------------------------
!
temp_thresh = MIN(t_gnd_k(i,j),t_sfc_k(i,j)-10.0)
IF(temp_3d(i,j,k) < temp_thresh)THEN
cldcv(i,j,k)=.01 ! not in inversion, clear it out
END IF
END IF ! tg-tb>8 IR signature present
END IF ! cldht-topo.gt.sfc_sao_buffer is .true.
END IF ! tb8_k(i,j).gt.tb8_est
END IF ! Current Cloud Cover is significant (> .04)
END DO ! k (for clearing clouds)
!
!-----------------------------------------------------------------------
!
! Insert satellite clouds
!
!-----------------------------------------------------------------------
!
IF(l_cloud_present) THEN ! Insert satellite clouds
!
!-----------------------------------------------------------------------
!
! Set initial satellite cloud base.
!
!-----------------------------------------------------------------------
!
htbase_init=cldtop_m(i,j) - thk_def
htbase = htbase_init
!
!-----------------------------------------------------------------------
!
! Locate lowest SAO cloud base
!
!-----------------------------------------------------------------------
!
!wdt update: only compute ht_sao_base once for each grid point
IF (ht_sao_base(i,j) < -1E20) THEN ! haven't tried yet
ht_sao_base(i,j) = 1E25
cldcv_above = cldcv_sao(i,j,nz-1)
DO k=nz-2,1,-1
IF (cldcv_sao(i,j,k) <= thr_sao_cvr &
.AND. cldcv_above > thr_sao_cvr) THEN
ht_sao_base(i,j) = zs(i,j,k+1)
END IF
cldcv_above = cldcv_sao(i,j,k)
END DO ! k
ENDIF
i_sao = i
j_sao = j
ht_sao_base_min = ht_sao_base(i,j)
!jz IF (ht_sao_base.eq.1e30) THEN ! Srch for nearby SAO cld lyrs
IF (ht_sao_base(i,j) > 1E20) THEN ! Srch for nearby SAO cld lyrs
! because the sat. says cld
! and the SAO doesn't
n_no_sao1 = n_no_sao1 + 1 ! Sat.cld but no SAO cld at same pt.
ht_sao_base_min = -1E25
idelt_index = 1
DO ! search for neighboring ht_sao_base values
IF (idelt_index > idelt_max .OR. ht_sao_base_min > -1E20) EXIT
DO jj = j-idelt_index,j+idelt_index
DO ii = i-idelt_index,i+idelt_index,2*idelt_index
IF(ii >= 1.AND.ii <= nx-1 .AND. jj >= 1.AND.jj <= ny-1 ) THEN
IF (ht_sao_base(ii,jj) < -1E20) THEN
! ht_sao_base(ii,jj) not computed yet
ht_sao_base(ii,jj) = 1E25
cldcv_above = cldcv_sao(ii,jj,nz-1)
DO k = nz-2,1,-1
IF (cldcv_sao(ii,jj,k) <= thr_sao_cvr &
.AND. cldcv_above > thr_sao_cvr) THEN
ht_sao_base(ii,jj) = zs(i,j,k+1)
END IF
cldcv_above = cldcv_sao(ii,jj,k)
END DO ! k
END IF ! compute ht_sao_base(ii,jj)
IF(ht_sao_base(ii,jj) < 1E20 .AND. &
ht_sao_base(ii,jj) < ht_sao_base_min)THEN
i_sao = ii
j_sao = jj
ht_sao_base_min = ht_sao_base(ii,jj)
END IF
END IF ! In bounds
END DO
END DO
DO jj = j-idelt_index,j+idelt_index,2*idelt_index
DO ii = i-idelt_index+1,i+idelt_index-1
IF(ii >= 1.AND.ii <= nx-1 .AND. jj >= 1.AND.jj <= ny-1 ) THEN
IF (ht_sao_base(ii,jj) < -1E20) THEN
! ht_sao_base(ii,jj) not computed yet
ht_sao_base(ii,jj) = 1E25
cldcv_above = cldcv_sao(ii,jj,nz-1)
DO k = nz-2,1,-1
IF (cldcv_sao(ii,jj,k) <= thr_sao_cvr &
.AND. cldcv_above > thr_sao_cvr) THEN
ht_sao_base(ii,jj) = zs(i,j,k+1)
END IF
cldcv_above = cldcv_sao(ii,jj,k)
END DO ! k
END IF ! compute ht_sao_base(ii,jj)
IF(ht_sao_base(ii,jj) < 1E20 .AND. &
ht_sao_base(ii,jj) < ht_sao_base_min)THEN
i_sao = ii
j_sao = jj
ht_sao_base_min = ht_sao_base(ii,jj)
END IF
END IF ! In bounds
END DO
END DO
idelt_index = idelt_index + 1
END DO
IF (ht_sao_base_min < -1E20) ht_sao_base_min = ht_sao_base(i,j)
END IF ! ht_sao_base.eq.1e30
IF (ht_sao_base_min > 1E20) THEN ! Sat. cld but no SAO cld
n_no_sao2 = n_no_sao2 + 1 ! even in neighbor points.
cover=sat_cover
htbase_init = ht_sao_base_min
IF(tb8_k(i,j)-t_gnd_k(i,j) < -21.) THEN ! We more likely
! have a cloud
buffer = 2100.
ELSE
buffer = sfc_ir_buffer ! Weed out IR tops < ~5000m AGL
END IF
!
!-----------------------------------------------------------------------
!
! Calculate new cloud top and cover
!
!-----------------------------------------------------------------------
!
cldtop_m_old = cldtop_m(i,j)
!
CALL cloud_top(nx,ny,nz, &
i,j,zs_1d,t_1d,p_1d,z_ref_lcl,z_lcl(i,j), &
topo(i,j),tmin,r_missing, &
tb8_cold_k(i,j),t_gnd_k(i,j),thresh2, &
cldtop_m_tb8(i,j),cldtop_m(i,j), &
l_tb8,l_cloud_present,sat_cover)
!
!-----------------------------------------------------------------------
!
! Change to cover
!
!-----------------------------------------------------------------------
!
cover=band8_cover(tb8_k(i,j),t_gnd_k(i,j),tb8_cold_k(i,j))
htbase = MAX(topo(i,j)+buffer, cldtop_m(i,j)-thk_def )
i_reset_base = i_reset_base +1
IF (htbase > cldtop_m(i,j)) THEN
n_no_sao3 = n_no_sao3 + 1 ! Sat.cld, but no SAO cld,
! & sat.cld is too low.
ELSE IF(i_reset_base/50*50 == i_reset_base) THEN
WRITE(6,611) i,j
611 FORMAT(/1X,' Correction at point(',2I2 &
,') because there is low sat.cld,' &
,' but no SAO cld.')
WRITE(6,612,ERR=212) tb8_k(i,j),tb8_cold_k(i,j) &
,cldtop_m_old,cldtop_m(i,j),cover
612 FORMAT(1X,'tb=',f8.1,' tb_cold=',f8.1, &
' ctop=',f9.0,' ctop_cold=',f9.0,' cvr=',f5.2/)
212 END IF
ELSE IF (ht_sao_base_min > cldtop_m(i,j)) THEN ! Sat.cld, but no
! SAO cld at same pt. Do have
! SAO cld in nearby pt. AND
! Sat.top is below ceiling
cover=sat_cover
htbase_init = ht_sao_base_min ! using SAO cld base
htbase = htbase_init
cldtop_m_old = cldtop_m(i,j)
cldtop_m(i,j) = htbase_init + thk_def ! correct sat. cldtop
!
!-----------------------------------------------------------------------
!
! Find a thinner value for cloud cover consistent with the new
! higher cloud top and the known brightness temperature.
! Note that cover is not really used here as an input
!
!-----------------------------------------------------------------------
!
CALL correct_cover(cover,cover_new &
,cldtop_m_old,cldtop_m(i,j),i,j,nx,ny,nz &
,zs_1d,t_1d,tb8_k(i,j),t_gnd_k(i,j) &
,istatus)
!
IF (istatus /= 1) THEN
WRITE(6,*)' Error in correct_cover'
RETURN
END IF
cover = cover_new
ELSE ! Normal use of satellite data
! There is a ht_sao_base below cldtop_m in the area
! near the grid pt.
cover=sat_cover
!
!-----------------------------------------------------------------------
!
! Locate SAO cloud base below satellite cloud top, modify
! satellite cloud base. Highest SAO ceiling within default
! thickness range of satellite layer is used.
!
!-----------------------------------------------------------------------
!
DO k=nz-1,1,-1
IF (zs(i,j,k) >= htbase_init .AND. &
zs(i,j,k) <= cldtop_m(i,j)) THEN
IF (cldcv_sao(i_sao,j_sao,k) <= thr_sao_cvr .AND. &
cldcv_sao(i_sao,j_sao,k+1) > thr_sao_cvr) THEN
!c We have an SAO base
htbase = zs(i,j,k+1)
!c If SAO (hence satellite) base is above the satellite
!c cloud top, lower the satellite base by one grid level
IF (htbase > cldtop_m(i,j)) htbase = zs(i,j,k)
GO TO 301
END IF ! We have an SAO base
END IF ! in satellite layer
END DO ! k
END IF ! HT_SAO_BASE .EQ. 1E30
! No SAO cloud base in the area near the grid pt.
301 IF (htbase /= htbase_init .AND. &
i_reset_base/50*50 == i_reset_base) THEN
WRITE(6,*)' Satellite ceiling reset by SAO.'
WRITE(6,660) i,j,htbase,htbase_init,cldtop_m(i,j)
660 FORMAT(1X,'i/j=',2I3,' htbase_new/old=',2F7.0, &
' cldtop_ir=',f7.0)
END IF
!
!-----------------------------------------------------------------------
!
! Add satellite cloud to array
!
!-----------------------------------------------------------------------
!
DO k=nz,1,-1
IF (zs(i,j,k) >= htbase .AND. &
zs(i,j,k) <= cldtop_m(i,j)) & ! in satellite layer
cldcv(i,j,k)=cover
END DO
END IF ! l_cloud_present (Cloudy)
END IF ! non-missing IR data
END DO ! i=1,nx
END DO ! j=1,ny
!wdt tmp
! write (29) cldtop_m_tb8 !XXX
! write (29) cldtop_m !XXX
! close(29)
istatus = 1
!
!-----------------------------------------------------------------------
!
! Write stats on Band 8 (11.2mm) methods
!
!-----------------------------------------------------------------------
!
WRITE(6,*)' n_no_sao (1/2/3) = ',n_no_sao1,n_no_sao2,n_no_sao3
CALL compare_rad (nx,ny,nz,clddiag,r_missing,zs,temp_3d &
,cldcv,t_sfc_k,t_gnd_k,tb8_k &
,cvr_snow,nlyr)
!
!-----------------------------------------------------------------------
! 999 RETURN
RETURN
END SUBROUTINE insert_ir
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE CORRELATION ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE correlation(t,tb8_k,thresh,ni,nj) 1
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Checking the navigation of the IR satellite data
! It's only a diagnostic procedure.
!
!-----------------------------------------------------------------------
!
! AUTHOR: (Dan Birkenheuer?)
! 09/95.
!
! MODIFICATION HISTORY:
!
! 04/29/96 J. Zhang
! Modified for the ARPSDAS format.
! Added full documentation.
!
! 05/01/98 Keith Brewster
! Added check for missing data.
!
!-----------------------------------------------------------------------
!
! Variable Declaration
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
! INPUT:
INTEGER :: ni,nj ! Model grid size
REAL :: t(ni,nj) ! ground skin temperature
REAL :: tb8_k(ni,nj) ! satellite IR band 8 brightness temperature
REAL :: thresh ! threshold define the bias criteria
!
! LOCAL:
INTEGER :: ibox
PARAMETER (ibox = 3)
INTEGER*4 i_corr_array(-ibox:ibox,-ibox:ibox)
!
!-----------------------------------------------------------------------
!
! Misc local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: min_count,ioff,joff,i,j,ii,jj,icount &
,ioff_min,joff_min
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
WRITE(6,*)
WRITE(6,*)' Checking Navigation of IR satellite, thresh = ' &
,thresh
min_count = 999999
DO joff = -ibox,ibox
DO ioff = -ibox,ibox
icount = 0
DO j = 1,nj-1
jj = MAX(MIN(j+joff,nj-1),1)
DO i = 1,ni-1
ii = MAX(MIN(i+ioff,ni-1),1)
IF(tb8_k(ii,jj) > 0. .AND. (t(i,j) - tb8_k(ii,jj)) > thresh) &
icount = icount+1 ! Clouds Indicated
END DO ! i
END DO ! j
IF(icount < min_count)THEN
min_count = icount
ioff_min = ioff
joff_min = joff
END IF
i_corr_array(ioff,joff) = icount
END DO ! IOFF
END DO ! JOFF
WRITE(6,99)
99 FORMAT(//1X,' Number of cloudy points indicated' &
,' by offseting tb8:')
WRITE(6,100) (i,i=-ibox,ibox)
100 FORMAT(1X,'joff:ioff',7I5)
DO j = ibox,-ibox,-1
WRITE(6,101) j,(i_corr_array(i,j),i=-ibox,ibox)
101 FORMAT(1X,1X,i2,5X,7I5)
END DO
WRITE(6,*)
WRITE(6,201) min_count,ioff_min,joff_min
201 FORMAT(' Min. of',i5,' pts flagged with an offset of',2I4/)
RETURN
END SUBROUTINE correlation
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE CLOUD_TOP ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE cloud_top(nx,ny,nz, & 2
i,j,zs_1d,t_1d,p_1d,z_ref_lcl,z_lcl, &
topo,tmin,r_missing, &
tb8_k,t_gnd_k,thresh2, &
cldtop_m_tb8,cldtop_m, &
l_tb8,l_cloud_present,sat_cover)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! This routine computes the cloud top height given a band 8
! brightness temperature and 3D fields of temp and height.
!
!-----------------------------------------------------------------------
!
! AUTHOR: (Jian Zhang)
! 04/96 Based on the LAPS cloud analysis code (09/95 Dan Birkenheuer)
!
! MODIFICATION HISTORY:
!
! 04/26/96 (J. Zhang)
! Modified for the ARPSDAS format.
! Added full documentation.
!
! 11/01/96 (J. Zhang)
! Added a scheme for calculation of cloud top height of low clouds
! (e.g., persistent stratocumulus), where temperature inversion
! may exist and the simple matching scheme may fail.
! Reference:
! Macpherson et al., 1996: The impact of MOPS moisture data in
! the U.K. Meteorological Office mesoscale data assimilation
! scheme. MWR, vol. 124, 1746-1766
!
! 09/10/97 (J. Zhang)
! Lifting condensation level is input through calling
! argument.
!
!-----------------------------------------------------------------------
!
! Variable Declarition
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! For tb8 method
!
!-----------------------------------------------------------------------
!
! INPUT:
INTEGER :: nx,ny,nz ! The model grid size
!
! INPUT: Vertical 1D arrays
REAL :: zs_1d(nz) ! physical height (m) at scalar points
REAL :: t_1d(nz) ! 3D temperature analysis on model grid
REAL :: p_1d(nz) ! pres. analysis on model grid
!
! INPUT: At grid point (i,j)
INTEGER :: i,j ! grid point location indices
REAL :: z_lcl ! Lifting condensation level at i,j
REAL :: tb8_k ! band 8 beightness temperature(K)
REAL :: t_gnd_k ! gound temperature
REAL :: topo ! terrain height at the grid point
!
! INPUT: Constants
REAL :: tmin
REAL :: r_missing ! missing data flag value
REAL :: thresh2 ! Input from parent routine (=8K)
REAL :: z_ref_lcl ! ref. level for computing LCL
!
! OUTPUT:
LOGICAL :: l_cloud_present ! "cloudy" indicated by sate. IR data?
LOGICAL :: l_tb8 ! "cloudy" indicated by band 8 data?
REAL :: cldtop_m_tb8 ! cld top height(m) obtained from tb8_k
REAL :: cldtop_m ! cld top height(m) obtained from IR data
REAL :: sat_cover ! cloud fraction obtained from IR data
!
! LOCAL:
REAL :: gamma_s ! ref. moist adiabatic lapse rate
!
! LOCAL: For the conceptual model for the cloud top
REAL :: t_base,p_base_pa,zbase,ztop,dzinc,zht
PARAMETER(dzinc=100.0)
REAL :: p,press,temp,eso,dlvdt,des,dtz,rl,es
REAL :: tb8_loc
INTEGER :: nlevel,nlm1
!
! CONSTANTS:
REAL :: gamma_d ! dry adiabatic lapse rate (K/m)
PARAMETER (gamma_d = 9.8/1004.0)
!
!-----------------------------------------------------------------------
!
! Misc local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: kl,j1
REAL :: arg,frac_k,frac_z,z_ref,t_ref
!
!-----------------------------------------------------------------------
!
! Include files.
!
INCLUDE 'phycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
! Define constants
!
!-----------------------------------------------------------------------
!
dlvdt=-2.3693E+3 ! J/kg/K
eso=610.78 ! pa
!
!-----------------------------------------------------------------------
!
! This section finds the cloud top using Band 8 data
!
!-----------------------------------------------------------------------
!
cldtop_m_tb8 = r_missing ! or zeros
l_tb8 = .false.
!
IF (tb8_k - t_gnd_k < -thresh2) THEN ! probably clds
l_tb8 = .true.
IF (tb8_k < 253.15) GO TO 300 ! Using scan&match method
!
!-----------------------------------------------------------------------
!
! When Tb >= -20 deg C, find the cloud top height using a
! conceptual model presented in Macpherson et al. (1996).
!
!-----------------------------------------------------------------------
!
zbase = z_lcl
IF (zbase >= zs_1d(nz-1)) GO TO 300
! Using scanning & matching method
!
!-----------------------------------------------------------------------
!
! Find the model pressure and temperature at the lifting
! condensation level.
!
!-----------------------------------------------------------------------
!
DO kl = 2,nz-1
z_ref = z_ref_lcl + topo
IF (z_ref <= zs_1d(kl)) THEN
frac_z = (z_ref-zs_1d(kl-1))/(zs_1d(kl)-zs_1d(kl-1))
t_ref = t_1d(kl-1) + frac_z * (t_1d(kl) - t_1d(kl-1))
GO TO 50
END IF
END DO ! kl
PRINT*,'Error for Tbase,aborting....'
STOP
50 CONTINUE
t_base = t_ref - (z_lcl -z_ref)*gamma_d
!
DO kl = 2,nz-1
IF (zbase <= zs_1d(kl)) THEN
frac_z = (zbase-zs_1d(kl-1))/(zs_1d(kl)-zs_1d(kl-1))
arg = ALOG(p_1d(kl-1)) &
+ frac_z * (ALOG(p_1d(kl))-ALOG(p_1d(kl-1)))
p_base_pa = EXP(arg)
GO TO 100
END IF
END DO ! kl
PRINT*,'Error for Zbase,aborting....'
STOP
100 CONTINUE
IF (t_base <= tb8_k) GO TO 300
! Using scanning & matching method
!
!-----------------------------------------------------------------------
!
! Find the cloud top height, which is the level where the temp.
! reaches the satellite brightness temperature throught the moist
! adiabatic cooling process.
!
!-----------------------------------------------------------------------
!
temp = t_base ! temp = t_base_k
press = p_base_pa
gamma_s = 0.004 ! K/m
nlevel = (temp -tb8_k)/gamma_s/dzinc ! a guess for cloud top
IF(nlevel < 0) THEN
PRINT*,'Error in zlevel. aborting....'
PRINT*,' nlevel=',nlevel,' temp=',temp,' gam=',gamma_s, &
' dzinc=',dzinc,' tb8=',tb8_k
STOP
END IF
nlm1 = nlevel
IF(nlm1 < 1) nlm1=1
zht = zbase
!
DO j1=1,nlm1
rl = lathv+(273.15-temp)*dlvdt ! Lv as func of Temp.
arg = rl*(temp-273.15)/273.15/temp/rv
es = eso*EXP(arg) ! satur. vapor press.
arg = -g*dzinc/rd/temp
p = press*EXP(arg) ! hydrosta. press.decrease
!
!-----------------------------------------------------------------------
!
! Calculate saturated adiabatic lapse rate
!
!-----------------------------------------------------------------------
!
des = es*rl/temp/temp/rv
dtz = -g*((1.0+0.621*es*rl/(press*rd*temp))/ &
(cp+0.621*rl*des/press))
IF(ABS(dtz) < 1.0E-15) THEN
PRINT*,' Zero dt/dz:',dtz,' g=',g,' es=',es,' press=',press
PRINT*,'temp=',temp,' cp=',cp,' rl=',rl,' des=',des
STOP
END IF
zht = zht+dzinc ! moist adiabatic ascent every 100m
temp = temp+dtz*dzinc
IF(temp <= tb8_k) THEN ! Cloud top is found
ztop = zht + (tb8_k-temp)/dtz
GO TO 150
END IF
press = p
END DO ! j=1,nlm1
ztop = zht + (tb8_k-temp)/dtz
150 CONTINUE
!
!-----------------------------------------------------------------------
!
! Apply quality check to the cloud top height.
!
!-----------------------------------------------------------------------
!
IF(ztop <= zs_1d(1)) THEN
PRINT*,' Error, ztop for cloud top model is out of' &
,' bound at grid pt(',i,j,').'
PRINT*,' ztop=',ztop,' topo=',topo, &
' z_bottom=',zs_1d(1)
PRINT*,'aborting...'
STOP
END IF
!
IF(ztop > zs_1d(nz-1)) THEN
PRINT*,' Warning, ztop for cloud top model is out of' &
,' bound at grid pt(',i,j,').'
PRINT*,' ztop=',ztop,' topo=',topo &
,' z_top=',zs_1d(nz-1)
PRINT*,' ztop is set to the highest model level'
cldtop_m_tb8 = zs_1d(nz-1)
GO TO 900
END IF
cldtop_m_tb8 = ztop
GO TO 900
!
300 CONTINUE
!
!-----------------------------------------------------------------------
!
! When Tb < -20 deg C, use the scheme that finds the cloud top
! height by scanning the model temperature profile for a value
! matching tb8.
!
!-----------------------------------------------------------------------
!
tb8_loc=tb8_k
IF(tb8_k <= tmin) THEN
PRINT*,' Warning: cloud top colder than tmin in column'
WRITE(6,*) ' i j tb8_k t_top'
WRITE(6,600) i,j,tb8_k,t_1d(nz-1)
600 FORMAT(1X,2I3,2F8.2)
PRINT*,' Cloud top temperature is set to tmin'
tb8_loc=tmin
END IF
!
DO kl = nz-2,1,-1
!
!-----------------------------------------------------------------------
!
! Find the lowest temp. crossing point in the model temperature
! profile.
!
!-----------------------------------------------------------------------
!
IF( (t_1d(kl)-tb8_loc) * (t_1d(kl+1)-tb8_loc) <= 0.) THEN ! Crossing Pt
frac_k = (tb8_loc - t_1d(kl)) &
/ (t_1d(kl+1) - t_1d(kl))
arg = zs_1d(kl) + frac_k * (zs_1d(kl+1)-zs_1d(kl))
IF(arg >= topo) THEN
cldtop_m_tb8 = arg
ELSE
WRITE(6,*)' Warning: Cloud Top Below Ground - flagged'
WRITE(6,601)
601 FORMAT(1X,' i j kl tb8_loc t_abv t_blw frac' &
,' hgt_m topo ctop_8 ')
WRITE (6,602) i,j,kl,tb8_loc,t_1d(kl+1) &
,t_1d(kl),frac_k,arg,topo,cldtop_m_tb8
602 FORMAT(1X,3I3,3F7.2,f5.2,3F8.0)
END IF
END IF ! Crossing point found
END DO ! kl
ELSE ! No clouds according to SATELLITE (Band 8 - 11.2mm)
l_tb8 = .false.
END IF ! tb8_k-t_gnd_k .lt. -thresh2
!
!-----------------------------------------------------------------------
!
! Set variables
!
!-----------------------------------------------------------------------
!
900 CONTINUE
l_cloud_present = l_tb8
cldtop_m = cldtop_m_tb8
sat_cover = 1.0
RETURN
END SUBROUTINE cloud_top
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE CORRECT_COVER ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE correct_cover(cover_in,cover_new_f & 1,2
,cldtop_old,cldtop_new,i,j,nx,ny,nz &
,zs_1d,t_1d,tb8_k,t_gnd_k &
,istatus)
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Find a thinner value for cloud cover consistent with the new
! higher cloud top and the known brightness temperature.
!
!-----------------------------------------------------------------------
!
! AUTHOR: (Dan Birkenheuer?)
! 09/95.
!
! MODIFICATION HISTORY:
!
! 04/29/96 (J. Zhang)
! Modified for the ARPSDAS format.
! Added full documentation.
!
!-----------------------------------------------------------------------
!
! Variable Declaration
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
! INPUT:
INTEGER :: i,j ! grid point index
INTEGER :: nx,ny,nz ! grid size
REAL :: t_1d(nz) ! temp. analysis field
REAL :: zs_1d(nz) ! phycisal heights (m) for a column of scalar
! grid points
REAL :: tb8_k ! satellite band 8 brightness temp.
REAL :: t_gnd_k ! ground skin temp.
REAL :: cover_in ! cloud cover before correction
REAL :: cldtop_old ! cloud top height before correction
!
! OUTPUT:
INTEGER :: istatus
REAL :: cover_new_f ! cloud cover after correction
REAL :: cldtop_new ! cloud top height after correction
!
!-----------------------------------------------------------------------
!
! Misc local variables
!
!-----------------------------------------------------------------------
!
REAL :: z_temp,temp_old,temp_new,frac
REAL :: band8_cover ! a function
REAL :: cover_old,cover_new
INTEGER :: iz_temp
INTEGER :: iwrite
DATA iwrite /0/
SAVE iwrite
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
! Find Temperature of old cloud top
!
!-----------------------------------------------------------------------
!
CALL hgt_to_zk (cldtop_old,z_temp,nz,zs_1d,istatus)
IF(istatus /= 1)THEN
PRINT*,' Old cldtop is out of domain range:', z_temp
PRINT*,' Old cldtop=',cldtop_old,' zs range from',zs_1d(1) &
,' to ',zs_1d(nz-1)
END IF
z_temp = MAX(1.,MIN(z_temp,FLOAT(nz-1)-.001))
iz_temp = INT(z_temp)
frac = z_temp - iz_temp
temp_old = t_1d(iz_temp) * (1. - frac) &
+ t_1d(iz_temp+1) * frac
!
!-----------------------------------------------------------------------
!
! Find Temperature of new cloud top
!
!-----------------------------------------------------------------------
!
CALL hgt_to_zk (cldtop_new,z_temp,nz,zs_1d,istatus)
IF(istatus /= 1)THEN
PRINT*,' New cldtop is out of domain range:', z_temp
PRINT*,' New cldtop=',cldtop_new,' zs range from',zs_1d(1) &
,' to ',zs_1d(nz-1)
! return
END IF
z_temp = MAX(1.,MIN(z_temp,FLOAT(nz-1)-.001))
iz_temp = INT(z_temp)
frac = z_temp - iz_temp
temp_new = t_1d(iz_temp) * (1. - frac) &
+ t_1d(iz_temp+1) * frac
!
!-----------------------------------------------------------------------
!
! This one utilizes a linear approximation to
! the sigma T**4 relationship
!
!-----------------------------------------------------------------------
!
cover_old = MIN(cover_in,1.0)
cover_new = cover_old * (tb8_k-t_gnd_k)/(temp_new-t_gnd_k)
!
!-----------------------------------------------------------------------
!
! This one utilizes the sigma T**4 relationship
!
!-----------------------------------------------------------------------
!
cover_new_f = band8_cover(tb8_k,t_gnd_k,temp_new)
IF((j-1) == INT(j-1)/10*10) THEN
iwrite = iwrite + 1
IF(iwrite < 15) THEN
WRITE(6,601,ERR=2)
601 FORMAT(//1X,'**CORR_CVR** i j tg_k told tnew ' &
,' ctold ctnew cvnew cvrnewf')
WRITE(6,1,ERR=2) i,j,t_gnd_k,temp_old,temp_new,cldtop_old &
,cldtop_new,cover_new,cover_new_f
1 FORMAT(1X,12X,2I3,3F7.0,2F8.0,2F8.2)
2 END IF
END IF
istatus=1
RETURN
END SUBROUTINE correct_cover
!
!##################################################################
!##################################################################
!###### ######
!###### FUNCTION BAND8_COVER ######
!###### ######
!##################################################################
!##################################################################
!
FUNCTION band8_cover(tb8_k,t_gnd_k,t_cld)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Convert the radiance to the brightness temperature
!
!-----------------------------------------------------------------------
!
! AUTHOR:
! 07/95
!
! MODIFICATION HISTORY:
! 05/01/96 (Jian Zhang)
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
! INPUT:
REAL :: tb8_k ! satellite band 8 brightness temp.
REAL :: t_gnd_k ! ground skin temp.
REAL :: t_cld ! estimated brightness temp. with SAO cld
!
! OUTPUT:
REAL :: band8_cover ! modified cld cover using sat. band 8 data
!
! LOCAL:
REAL :: r_sfc,r_sat,r_cld
REAL :: temp_to_rad ! a function
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
r_sfc = temp_to_rad(t_gnd_k)
r_sat = temp_to_rad(tb8_k)
r_cld = temp_to_rad(t_cld)
band8_cover = (r_sat-r_sfc) / (r_cld-r_sfc)
RETURN
END FUNCTION band8_cover
!
!##################################################################
!##################################################################
!###### ######
!###### FUNCTION TEMP_TO_RAD ######
!###### ######
!##################################################################
!##################################################################
!
FUNCTION temp_to_rad(temp),1
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Convert the radiance to the brightness temperature
!
!-----------------------------------------------------------------------
!
! AUTHOR:
! 07/95
!
! MODIFICATION HISTORY:
! 05/01/96 (Jian Zhang)
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
! INPUT:
REAL :: temp
!
! OUTPUT:
REAL :: temp_to_rad
!
! LOCAL:
INTEGER :: init
DATA init /0/
SAVE init
INTEGER :: nsat
!
! FUNCTIONS: (in /vortex/lib/goeslib/invers.f)
REAL :: vplanc ! function that convert a brightness temperature
! to the radiance for a specific band
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
INCLUDE 'adas.inc'
IF(init == 0)THEN
init = 1
nsat = 3
PRINT*,' Calling PLNKIV. calibration file:',calib_fname
CALL plnkiv(nsat,calib_fname)
PRINT*,' End calling PLNKIV'
END IF
temp_to_rad = temp
temp_to_rad = vplanc(temp,8)
RETURN
END FUNCTION temp_to_rad
!
!##################################################################
!##################################################################
!###### ######
!###### FUNCTION RAD_TO_TEMP ######
!###### ######
!##################################################################
!##################################################################
!
FUNCTION rad_to_temp(rad),1
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! Convert the radiance to the brightness temperature
!
!-----------------------------------------------------------------------
!
! AUTHOR:
! 07/95
!
! MODIFICATION HISTORY:
! 05/01/96 (Jian Zhang)
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
! INPUT:
REAL :: rad
!
! OUTPUT:
REAL :: rad_to_temp
!
! LOCAL:
INTEGER :: init
DATA init /0/
SAVE init
INTEGER :: nsat
!
! FUNCTIONS: (in /vortex/lib/goeslib/invers.f)
REAL :: vbrite ! function that convert a radiance to
! a brightness temp for a specific band
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
INCLUDE 'adas.inc'
IF(init == 0) THEN
init = 1
nsat = 3
CALL plnkiv(nsat,calib_fname)
END IF
rad_to_temp = rad
rad_to_temp = vbrite(rad,8)
RETURN
END FUNCTION rad_to_temp
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE COMPARE_RAD ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE compare_rad (nx,ny,nz,clddiag,r_missing,zs,t_3d & 1,6
,cldcv,t_sfc_k,t_gnd_k,tb8_k &
,cvr_snow,nlyr)
!
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! This routine compares the analyzed clouds to the 11.2mm radiation
! and determines adjusted cloud fractions of the cloud layers to yield
! a better fit.
!
!-----------------------------------------------------------------------
!
! AUTHOR: (?)
! 09/95.
!
! MODIFICATION HISTORY:
!
! 05/01/96 J. Zhang
! Modified for the ARPSDAS format.
! Added full documentation.
! 05/01/99 K. Brewster
! Added check for missing IR data.
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! INPUT:
INTEGER :: nx,ny,nz ! grid size
!
! INPUT:
INTEGER :: clddiag ! print diagnostic fields?
REAL :: r_missing ! bad flag
REAL :: t_3d(nx,ny,nz) ! temp. analysis on the model grid
REAL :: zs(nx,ny,nz) ! physical heights (m) at model scalar pts.
REAL :: cvr_snow(nx,ny) ! sfc snow coverage fraction
REAL :: tb8_k(nx,ny) ! satellite band8 brightness temp.
REAL :: t_gnd_k(nx,ny) ! ground skin temp.
REAL :: t_sfc_k(nx,ny) ! surface air temperature
!
! OUTPUT:
REAL :: cldcv(nx,ny,nz) ! adjusted cloud cover analysis
!
! LOCAL:
REAL :: zs_1d(nz) ! phy. heights (m) in a column of grid
REAL :: t_1d(nz) ! temp in a column of grid
REAL :: cldcv_1d(nz) ! cloud cover in one culomn of grid.
REAL :: cldcv_1dref(nz) ! cloud cover in one culomn of grid.
REAL :: a(nz) ! max. cldcvr in each cloud layer
REAL :: f(nz)
REAL :: a_new(nz)
REAL :: f_new(nz)
INTEGER :: ilyr(nz) ! Dims needs to be changed to nz
INTEGER :: ilyr_new(nz) ! Dims needs to be changed to nz
INTEGER :: nlyr,nlyr_new
LOGICAL :: l_correct
!
!-----------------------------------------------------------------------
!
! Misc local variables
!
!-----------------------------------------------------------------------
!
REAL :: corr_thr,cover_step
REAL :: delta_cover,delta_cover_ref
INTEGER :: i_correct,iter,iter_max,n_iter
REAL :: t_effective,tdiff,tdiff_ref
INTEGER :: i,j,k,l,iwrite,imid,jmid,kmid,ibeg,iend
INTEGER :: n_clear,n_clear_ns,n_clear_sn,n_cloudy,n_total &
,n_correct
REAL :: tdiff_sumsq,tdiff_cld_sumsq,tdiff_corr_sumsq &
,tdiff_corr_sumsq2,tdiff_corr_sumsq3,tdiff_corr_cld_sumsq
REAL :: tb8_g_clr_sum,tb8_a_clr_sum,tb8_g_clr_sumsq &
,tb8_a_clr_sumsq,tb8_g_clr_ns_sum,tb8_a_clr_ns_sum &
,tb8_g_clr_ns_sumsq,tb8_a_clr_ns_sumsq,tb8_g_clr_sn_sum &
,tb8_a_clr_sn_sum,tb8_g_clr_sn_sumsq,tb8_a_clr_sn_sumsq
REAL :: frac,frac_clouds
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
imid=nx/2
jmid=ny/2
ibeg=MAX( 1,imid-5)
iend=MIN(nx,imid+5)
kmid=nz/2
IF(clddiag == 1) THEN
WRITE(6,'(a)')' Comparing radiation'
WRITE(6,611) nx,ny,nz
611 FORMAT(1X,' nx=',i3,' ny=',i3,' nz=',i3)
! 612 FORMAT(1X,10F7.0)
WRITE(6,'(a)')'grid point heights:'
WRITE(6,613) (zs(i,jmid,kmid),i=ibeg,iend)
613 FORMAT(10F8.0)
WRITE(6,'(a)')'grid point temp:'
WRITE(6,613) (t_3d(i,jmid,kmid),i=ibeg,iend)
WRITE(6,'(a)')
!
WRITE(6,'(a)') 'Surface temperatures: j=jmid, i=imid-5,+5'
WRITE(6,614) (t_sfc_k(i,jmid),i=ibeg,iend)
614 FORMAT(1X,10F7.0)
WRITE(6,'(a)') 'Ground skin temperatures: j=jmid, i=imid-5,+5'
WRITE(6,614) (t_gnd_k(i,jmid),i=ibeg,iend)
WRITE(6,'(a)') 'Band 8 temperatures: j=jmid, i=imid-5,+5'
WRITE(6,614) (tb8_k(i,jmid),i=ibeg,iend)
WRITE(6,'(a)')
END IF
corr_thr = 4. ! if diff of tb8-tb8e > 4K, correction
cover_step = 0.05 ! add/minus 0.05 cldcvr each correction
iter_max = 10 ! we'll do maximum of 10 iterns
iwrite = 0 ! output message flag
n_clear = 0 !
n_clear_ns = 0
n_clear_sn = 0
n_cloudy = 0
n_total = 0
n_correct = 0
n_iter = 0
tdiff_sumsq = 0. ! sum of sqr(tb8o-tb8e) for all pts
! before correction
tdiff_cld_sumsq = 0. ! sum of sqare(tb8o-tb8e) for cloudy pts
! before correction
tdiff_corr_sumsq3 = 0. ! sum of sqr(tb8o-tb8e) for cloudy pts
! to be corrected, before the correction
tdiff_corr_sumsq = 0. ! sum of sqr(tb8o-tb8e) for all pts
! after correction
tdiff_corr_cld_sumsq = 0. ! sum of sqr(tb8o-tb8e) for all cldy pts
! after correction
tdiff_corr_sumsq2 = 0. ! sum of sqr(tb8o-tb8e) for cloudy pts
! to be corrected, after the correction
tb8_g_clr_sum = 0. ! sum of tb8-tg for clear air pts
tb8_a_clr_sum = 0. ! sum of tb8-tsfc for clear air pts
tb8_g_clr_sumsq = 0. ! sum of sqare(tb8-tg) for clear air pts
tb8_a_clr_sumsq = 0. ! sum of sqare(tb8-tsfc) for clear air pts
tb8_g_clr_ns_sum = 0. ! for clear air with < 0.5 snow cover pts
tb8_a_clr_ns_sum = 0.
tb8_g_clr_ns_sumsq = 0.
tb8_a_clr_ns_sumsq = 0.
tb8_g_clr_sn_sum = 0. ! for clear air with > 0.5 snow cover pts
tb8_a_clr_sn_sum = 0.
tb8_g_clr_sn_sumsq = 0.
tb8_a_clr_sn_sumsq = 0.
!
!-----------------------------------------------------------------------
!
! Compare and correct (if necessary) the cloud cover analysis field
! for each grid point.
!
!-----------------------------------------------------------------------
!
DO j = 1,ny-1
DO i = 1,nx-1
IF(tb8_k(i,j) > 0.) THEN
!jz do k = 1,nz-1
DO k = 1,nz
zs_1d(k) = zs(i,j,k)
t_1d(k) = t_3d(i,j,k)
END DO
DO k = 1,nz
cldcv_1dref(k) = cldcv(i,j,k)
END DO
CALL cvr_to_tb8e (nz,zs_1d,t_1d &
,cldcv_1dref,t_gnd_k(i,j) &
,a,f,ilyr,t_effective,nlyr)
!
tdiff = tb8_k(i,j)-t_effective ! Band 8 - calculated
frac_clouds = 1.-f(nlyr) ! column total cloud cover
IF(nlyr == 1) THEN
tb8_g_clr_sum = tb8_g_clr_sum + tdiff
tb8_g_clr_sumsq = tb8_g_clr_sumsq + tdiff**2
tb8_a_clr_sum = tb8_a_clr_sum + tb8_k(i,j)-t_sfc_k(i,j)
tb8_a_clr_sumsq=tb8_a_clr_sumsq+(tb8_k(i,j)-t_sfc_k(i,j))**2
n_clear = n_clear + 1
IF(cvr_snow(i,j) /= r_missing) THEN
IF(cvr_snow(i,j) < 0.5) THEN
tb8_g_clr_ns_sum = tb8_g_clr_ns_sum + tdiff
tb8_g_clr_ns_sumsq = tb8_g_clr_ns_sumsq + tdiff**2
tb8_a_clr_ns_sum = tb8_a_clr_ns_sum &
+ tb8_k(i,j)-t_sfc_k(i,j)
tb8_a_clr_ns_sumsq = tb8_a_clr_ns_sumsq &
+ (tb8_k(i,j) - t_sfc_k(i,j))**2
n_clear_ns = n_clear_ns + 1
ELSE
tb8_g_clr_sn_sum = tb8_g_clr_sn_sum + tdiff
tb8_g_clr_sn_sumsq = tb8_g_clr_sn_sumsq + tdiff**2
tb8_a_clr_sn_sum = tb8_a_clr_sn_sum &
+ tb8_k(i,j)-t_sfc_k(i,j)
tb8_a_clr_sn_sumsq = tb8_a_clr_sn_sumsq &
+(tb8_k(i,j)-t_sfc_k(i,j))**2
n_clear_sn = n_clear_sn + 1
END IF
END IF
ELSE
n_cloudy = n_cloudy + 1
tdiff_cld_sumsq = tdiff_cld_sumsq + tdiff**2
END IF
n_total = n_total + 1
tdiff_sumsq = tdiff_sumsq + tdiff**2
!
!-----------------------------------------------------------------------
!
! Apply corrections?
!
! Only when the following four conditions are satisfied:
! 1. at least one cloud layer
! 2. the difference between the observed and estimated tb8
! is larger than the threshold (4K)
! 3. tb8_obs - t_gnd < -8K
! 4. column total cloud cover is larger than 40%
!
!-----------------------------------------------------------------------
!
IF (nlyr >= 2 .AND. ABS(tdiff) > corr_thr &
.AND. t_gnd_k(i,j) - tb8_k(i,j) > 8. &
.AND. frac_clouds > 0.4 ) THEN
l_correct = .true.
n_correct = n_correct + 1
tdiff_corr_sumsq3 = tdiff_corr_sumsq3 + tdiff**2
ELSE
l_correct = .false.
END IF
!
!-----------------------------------------------------------------------
!
! This corrective section is now turned on
!
!-----------------------------------------------------------------------
!
iter = 0
delta_cover = 0.
905 IF (l_correct) THEN
iter = iter + 1
tdiff_ref = tdiff ! save initial difference
delta_cover_ref = delta_cover ! save initial correction
IF(tdiff < 0.) THEN
delta_cover = delta_cover + cover_step
ELSE
delta_cover = delta_cover - cover_step
END IF
CALL apply_correction (cldcv_1dref,cldcv_1d &
,nz,ilyr,f,delta_cover)
CALL cvr_to_tb8e (nz,zs_1d,t_1d &
,cldcv_1d,t_gnd_k(i,j) &
,a_new,f_new,ilyr_new,t_effective,nlyr_new)
!
tdiff = tb8_k(i,j)-t_effective
frac_clouds = 1.-f_new(nlyr_new)
!
!-----------------------------------------------------------------------
!
! Continue to apply corrections?
!
!-----------------------------------------------------------------------
!
IF(nlyr_new >= 2 .AND. frac_clouds > 0.4) THEN
i_correct = 1
ELSE
i_correct = 0
END IF
IF(MOD(iwrite,50) == 0) THEN
iwrite = iwrite + 1
WRITE(6,*) 'iter=',iter ,' tdiff_ref=',tdiff_ref &
,'tg=',t_gnd_k(i,j)
WRITE(6,*) (a(l),l=nlyr-1,1,-1)
WRITE(6,641)
641 FORMAT(1X,' i j ncn tb8o tb8e tdiff cldcvm i_c' &
,6(' cvrln'))
WRITE(6,642,ERR=912)i,j,nlyr_new-1,tb8_k(i,j),t_effective &
,tdiff,frac_clouds,i_correct &
,(a_new(l),l=nlyr_new-1,1,-1)
642 FORMAT(1X,2I3,i4,2F7.0,f7.1,f7.3,i2,2X,10F6.2)
912 END IF
IF(i_correct == 1 .AND. iter < iter_max &
.AND. ABS(tdiff) < ABS(tdiff_ref) &
.AND. tdiff * tdiff_ref > 0.) GO TO 905
! Loop back & increment cover
n_iter = n_iter + iter
!
!-----------------------------------------------------------------------
!
! Final iteration
!
!-----------------------------------------------------------------------
!
IF(.false. .OR. tdiff*tdiff_ref >= 0.) THEN
! Select best of last two increments
IF(ABS(tdiff) >= ABS(tdiff_ref)) THEN
delta_cover = delta_cover_ref
CALL apply_correction (cldcv_1dref,cldcv_1d &
,nz,ilyr,f,delta_cover)
tdiff = tdiff_ref
END IF
ELSE ! Do one Newton iteration
frac = tdiff / (tdiff - tdiff_ref)
delta_cover=delta_cover+frac*(delta_cover_ref-delta_cover)
CALL apply_correction (cldcv_1dref,cldcv_1d &
,nz,ilyr,f,delta_cover)
CALL cvr_to_tb8e (nz,zs_1d,t_1d &
,cldcv_1d,t_gnd_k(i,j) &
,a_new,f_new,ilyr_new,t_effective,nlyr_new)
tdiff = tb8_k(i,j)-t_effective
n_iter = n_iter + 1
END IF
!
!-----------------------------------------------------------------------
!
! Put the corrected column cloud coverback to 3D cloud cover field
!
!-----------------------------------------------------------------------
!
DO k=1,nz
cldcv(i,j,k) = cldcv_1d(k)
END DO
END IF ! Corrections were made
tdiff_corr_sumsq = tdiff_corr_sumsq + tdiff**2
IF(nlyr > 1) & ! cloudy (at least before adjustment)
tdiff_corr_cld_sumsq = tdiff_corr_cld_sumsq + tdiff**2
IF(l_correct) THEN
tdiff_corr_sumsq2 = tdiff_corr_sumsq2 + tdiff**2
END IF
END IF
END DO ! I
END DO ! J
!
!-----------------------------------------------------------------------
!
! Write out statistics on consistency of Band 8 data and cloud cover
!
!-----------------------------------------------------------------------
!
IF(n_clear > 0) THEN
WRITE(6,951,ERR=960) n_clear, tb8_g_clr_sum/FLOAT(n_clear) &
,SQRT(tb8_g_clr_sumsq/FLOAT(n_clear))
951 FORMAT(/' Mean/RMS band 8/gnd temp residual in clear skies =' &
,i5,2F9.3)
960 WRITE(6,961,ERR=970) n_clear, tb8_a_clr_sum/FLOAT(n_clear) &
,SQRT(tb8_a_clr_sumsq/FLOAT(n_clear))
961 FORMAT(' Mean/RMS band 8/air temp residual in clear skies =' &
,i5,2F9.3)
970 END IF
IF(n_clear_ns > 0) THEN
WRITE(6,956,ERR=965) n_clear_ns &
,tb8_g_clr_ns_sum/FLOAT(n_clear_ns) &
,SQRT(tb8_g_clr_ns_sumsq/FLOAT(n_clear_ns))
956 FORMAT(/' Mean/RMS band 8/gnd temp resid in clear/nsnow skies =' &
,i5,2F9.3)
965 WRITE(6,966,ERR=975) n_clear_ns &
,tb8_a_clr_ns_sum/FLOAT(n_clear_ns) &
,SQRT(tb8_a_clr_ns_sumsq/FLOAT(n_clear_ns))
966 FORMAT(' Mean/RMS band 8/air temp resid in clear/nsnow skies =' &
,i5,2F9.3)
975 END IF
IF(n_clear_sn > 0) THEN
WRITE(6,976,ERR=985) n_clear_sn &
,tb8_g_clr_sn_sum/FLOAT(n_clear_sn) &
,SQRT(tb8_g_clr_sn_sumsq/FLOAT(n_clear_sn))
976 FORMAT(/' Mean/RMS band 8/gnd temp resid in clear/snow skies =' &
,i5,2F9.3)
985 WRITE(6,986,ERR=995) n_clear_sn &
,tb8_a_clr_sn_sum/FLOAT(n_clear_sn) &
,SQRT(tb8_a_clr_sn_sumsq/FLOAT(n_clear_sn))
986 FORMAT(' Mean/RMS band 8/air temp resid in clear/snow skies =' &
,i5,2F9.3)
995 END IF
IF(n_total > 0) THEN
WRITE(6,971,ERR=980)n_total,SQRT(tdiff_sumsq/FLOAT(n_total))
971 FORMAT(/' RMS band 8/teff residual (bfr corr) in all skies =' &
,i5, f9.3)
980 WRITE(6,981,ERR=990) n_total, &
SQRT(tdiff_corr_sumsq/FLOAT(n_total))
981 FORMAT(' RMS band 8/teff residual (aft corr) in all skies =' &
,i5, f9.3)
990 END IF
IF(n_cloudy > 0) THEN
WRITE(6,991,ERR=1000) n_cloudy, &
SQRT(tdiff_cld_sumsq/FLOAT(n_cloudy))
991 FORMAT(' RMS band 8/teff residual (bfr corr) in cldy skies =' &
,i5, f9.3)
1000 WRITE(6,1001,ERR=1010) n_cloudy, &
SQRT(tdiff_corr_cld_sumsq/FLOAT(n_cloudy))
1001 FORMAT(' RMS band 8/teff residual (aft corr) in cldy skies =' &
,i5, f9.3)
1010 END IF
IF(n_correct > 0) THEN
WRITE(6,1011,ERR=1020) &
n_correct,SQRT(tdiff_corr_sumsq3/FLOAT(n_correct))
1011 FORMAT(/' RMS band 8/teff resid (bfr corr - corrected pts) =' &
,i5, f9.3)
1020 WRITE(6,1021,ERR=1030) &
n_correct,SQRT(tdiff_corr_sumsq2/FLOAT(n_correct))
1021 FORMAT(' RMS band 8/teff resid (aft corr - corrected pts) =' &
,i5, f9.3)
WRITE(6,*)' Total/Average # of iterations = ',n_iter, &
FLOAT(n_iter)/FLOAT(n_correct)
WRITE(6,*)
1030 END IF
RETURN
END SUBROUTINE compare_rad
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE CVR_TO_TB8E ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE cvr_to_tb8e (nz,zs_1d,t_1d & 3
,cvr,t_gnd_k &
,a,f,ilyr,t_effective,nlyr)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! To calculate the estimated band8 brightness temperature using
! the cloud cover analysis.
!
!-----------------------------------------------------------------------
!
! AUTHOR: (Jian Zhang)
! 04/1996 Based on the LAPS cloud analysis code of 09/95
!
! MODIFICATION HISTORY:
!
! 04/26/96 J. Zhang
! Modified for the ARPSDAS format. Added full
! documentation.
! 07/27/96 J. Zhang
! Modified code so that the cloud layer with top at the
! upper boundary of the cloud grid is counted.
! 04/11/97 J. Zhang
! Include adascld24.inc for ncloud
! 08/06/97 J. Zhang
! Change adascld24.inc to adascld25.inc
! 09/11/97 J. Zhang
! Change adascld25.inc to adascld26.inc
! 04/20/98 J. Zhang
! Based on the arps4.3.3 version. Abandoned cloud grid,
! using the model grid.
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! INPUT
INTEGER :: nz
!
! INPUT:
REAL :: cvr(nz) ! Cloud cover from analysis
REAL :: zs_1d(nz)
REAL :: t_1d(nz)
REAL :: t_gnd_k
!
! OUTPUT:
REAL :: t_effective ! effective brightness temp. estimated
! from cloud cover analysis
INTEGER :: nlyr
INTEGER :: ilyr(nz) ! Layer index for each cloud lvl (needs nz)
REAL :: a(nz) ! Cloud fractions of layers
REAL :: f(nz) ! Apparnt "x-sectn" of cldlyrs seen from above
!
! LOCAL:
INTEGER :: ik(nz) ! Height level representative of cloud layers
REAL :: temp_lyr(nz) ! temp. at "ik" lavels
!
! FUNCTIONS:
REAL :: temp_to_rad
REAL :: rad_to_temp
!
!-----------------------------------------------------------------------
!
! Misc local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: k,n
REAL :: rsum,sumf
REAL :: cvr_thresh
PARAMETER (cvr_thresh = 0.1)
!
!-----------------------------------------------------------------------
!
! Convert from cld cvr to discreet cloud layer indices (cvr to a)
!
! Find the maximum cloud cover "a" in each cloud deck, and
! the level index "ik" associated with the maximum cloud cover.
!
!-----------------------------------------------------------------------
!
nlyr = 0
IF(cvr(nz) >= cvr_thresh) THEN ! cldtop at the upper boundary
nlyr = nlyr + 1
a(nlyr) = cvr(nz)
ik(nlyr) = nz
ilyr(nz) = nlyr
ELSE
ilyr(nz)=0
END IF
DO k = nz-1,1,-1
IF(cvr(k) >= cvr_thresh .AND. cvr(k+1) < cvr_thresh) THEN
nlyr = nlyr + 1
a(nlyr) = cvr(k)
ik(nlyr) = k
ELSE
IF(nlyr >= 1) THEN
IF(cvr(k) > a(nlyr)) THEN
a(nlyr) = cvr(k) ! Max cldcvr within a layer
ik(nlyr) = k ! the lvl index for the max. cldcv
END IF
END IF
END IF
IF(cvr(k) >= cvr_thresh) THEN ! Still within layer
ilyr(k) = nlyr
ELSE ! Below layer
ilyr(k) = 0
END IF
END DO ! k
!
!-----------------------------------------------------------------------
!
! Get temperatures of the maximum cldcvr level in each cld deck.
!
!-----------------------------------------------------------------------
!
DO n = 1,nlyr
k = ik(n)
temp_lyr(n) = t_1d(k)
END DO ! n
!
!-----------------------------------------------------------------------
!
! Add a layer for the ground
!
!-----------------------------------------------------------------------
!
nlyr = nlyr + 1
a(nlyr) = 1.0
temp_lyr(nlyr) = t_gnd_k
!
!-----------------------------------------------------------------------
!
! Convert cloud layer fractions to "cross-section" seen from
! satellite. This solves for the "f" array given the "a" array
!
!-----------------------------------------------------------------------
!
a(1) = MIN( f(1) = a(1)
sumf = f(1)
IF(nlyr >= 2) THEN
DO n = 2,nlyr
a(n) = MIN( f(n) = a(n) * (1.0 - sumf) ! fraction of radiation reaches
! top of atm from each cld layer
sumf = sumf + f(n) ! fraction of radiation blked
! /attened by all cld lyrs above
END DO ! n
END IF ! nlyr
!
!-----------------------------------------------------------------------
!
! Calculate total radiance from all cloud layers + ground
!
!-----------------------------------------------------------------------
!
rsum = 0
DO n = 1,nlyr
rsum = rsum + temp_to_rad(temp_lyr(n)) * f(n)
END DO ! n
!
!-----------------------------------------------------------------------
!
! Convert to effective temperature and compare to
! the observed brightness temp
!
!-----------------------------------------------------------------------
!
t_effective = rad_to_temp(rsum)
RETURN
END SUBROUTINE cvr_to_tb8e
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE APPLY_CORRECTION ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE apply_correction (cldcv_1dref,cldcv_1d & 3
,nz,ilyr,f,delcv)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! To apply the correction to the cloud cover field.
! The amount for correction is an input parameter.
!
!-----------------------------------------------------------------------
!
! AUTHOR: (Jian Zhang)
! 05/02/96
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! INPUT:
INTEGER :: nz ! number of cloud analysis levels
REAL :: cldcv_1dref(nz) ! the column cldcvr to be corr.
INTEGER :: ilyr(nz) ! cld deck(layer) indices for ea.cldy level
REAL :: f(nz) ! fracn of rad. blocked by each cloud deck
REAL :: delcv ! the cld cvr increment added to the input
!
! OUTPUT:
REAL :: cldcv_1d(nz) ! corrected cloud cover
!
! LOCAL:
INTEGER :: k
!
!-----------------------------------------------------------------------
!
! Apply correction to 1D cloud cover field
!
!-----------------------------------------------------------------------
!
DO k = 1,nz
!jz if( ilyr(k).gt.0 .and. f(ilyr(k)).gt.0.0) then
IF(ilyr(k) > 0) THEN
IF( cldcv_1d(k) = MIN(cldcv_1dref(k)+delcv, 1.0)
cldcv_1d(k) = MAX(cldcv_1d(k), 0.0)
ELSE
cldcv_1d(k) = MAX(MIN(cldcv_1dref(k), 1.0), 0.0)
END IF
ELSE
cldcv_1d(k) = MAX(MIN(cldcv_1dref(k), 1.0), 0.0)
END IF
END DO
RETURN
END SUBROUTINE apply_correction