!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE CLOUD_CV ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
SUBROUTINE cloud_cv (nx,ny,nz,i4time,time,dirname,runname, & 1,24
xs,ys,zs,dx,dy,topo,latgr,longr, &
p_3d,t_3d,rh_3d, &
nobsng,stn,obstype,xsng,ysng, &
obstime,latsta,lonsta,elevsta, &
kcloud,store_amt,store_hgt, &
ref_mos_3d,istat_radar, &
clouds_3d,cloud_ceiling, &
istatus)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Perform 3D fractional cloud cover analysis using SAO
! cloud coverage observations, visible (VIS) and infrared (IR)
! satellite measurements, and radar reflectivity data.
!
!-----------------------------------------------------------------------
!
! AUTHOR: (Jian Zhang)
! 02/01/96
!
! MODIFICATION HISTORY:
!
! 07/12/96 J. Zhang
! Added remapping of the cloud cover fields onto the
! ARPS grid.
! 03/14/97 J. Zhang
! Cleaning up the code and implemented for the official
! arps4.2.4 version.
! 03/28/97 J. Zhang
! Using the lifiting condensation level as the cloud
! base when inserting radar reflectivity with no SAO
! cloud base.
! 04/11/97 J. Zhang
! Included dims.inc for nx,ny,nz.
! 04/15/97 J. Zhang
! Change the size of the arrays istn, jstn, hstn, and
! cstn from fixed (40) to adjustable (mx_sng)
! 08/06/97 J. Zhang
! Change adascld24.inc to adascld25.inc
! 09/04/97 J. Zhang
! Changed the radar echo thresholds for inserting clouds
! from radar reflectivities.
! 09/10/97 J. Zhang
! Calculate lifting condensation level in INICLDGRD,
! and used by other routines (INSERTIR and INSERTRAD).
! Change adascld25.inc to adascld26.inc
! 11/17/97 J. Zhang
! Change the length of "runname" from fixed number
! of letters (6) to flexible numbers.
! 11/18/97 J. Zhang
! Added flag 'clddiag' in adascld26.inc for whether to
! output the special cloud fields.
! 03/26/98 J. Zhang
! Some modifications based on the official adas 4.3.3.
! 04/27/98 J. Zhang
! Conform WRTCLDVAR to WRTVAR
! 05/05/98 J. Zhang
! Abandoned the cloud grid, using the ARPS grid instead.
! 10/25/98 K. Brewster
! Modified creation of file name for cloud observation
! diagnostic file.
! 08/21/01 K. Brewster
! Changed wrtvar calls to wrtvar1 to match new routine
! for writing 3-D arrays to be read by arpsplt as an
! arbitary array.
!
!-----------------------------------------------------------------------
!
! INCLUDE: (from adas.inc)
!
! mx_sng ! max. possible # of surface obs.
! max_cld_snd ! max. possible # of stations
! with cloud coverage reports.
! refthr1 ! "significant" radar echo at lower levels
! refthr2 ! "significant" radar echo at upper levels
! hgtrefthr ! height criteria for "significant" radar
! echo thresholds
! clddiag ! flag for whether to output the special cloud
! fields
!
! INPUT:
!
! nx,ny,nz ! ARPS grid size
!
! i4time ! analysis time in seconds from 00:00 UTC
! Jan. 1, 1960
! time ! analysis time in seconds from the model
! initial time
! dirname ! directory for data dump
!
! INPUT (for 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.
! temp_3d(nx,ny,nz)! the temperature field
!
! topo (nx,ny) ! The height of the terrain (equivalent
! dx, dy ! grid spacing
!
! INPUT for SAO obs.
!
! 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.
!
! LOCAL
!
! cf_modelfg (nx,ny,nz) ! first guess cloud cover field
! rh_modelfg (nx,ny,nz) ! first guess relative humidity field
! t_sfc_k (nx,ny) ! surface temperature field
!
! OUTPUT:
!
! istatus ! The flag indicating process status.
! clouds_3d (nx,ny,nz) ! final gridded cloud cover analysis
! cloud_ceiling (nx,ny) ! cloud ceiling (M AGL)
!
! LOCAL:
! 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.
! cvr_max (nx,ny) ! final column maximum cloud cover
! cvr_tot (nx,ny) ! final column total cloud cover
! cloud_top (nx,ny) ! cloud top (M ASL)
!
! LOCAL for SAO:
!
! ista_snd (max_cld_snd) index of cloud sounding stations
! i_snd (max_cld_snd) i-locn of each cld snd stn in ADAS grid
! j_snd (max_cld_snd) j-locn of each cld snd stn in ADAS grid
!
!
! LOCAL for gridded cloud cover analysis
!
! cldcv (nx,ny,nz) 3D gridded fractional cloud cover analysis.
! (when input, it's the first guess field)
! wtcldcv (nx,ny,nz) wts given to gridded fracn cld cvr anx.
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
! Include files:
!
INCLUDE 'adas.inc' ! ADAS parameters.
!
!-----------------------------------------------------------------------
!
! INPUT: general
INTEGER :: nx,ny,nz ! ARPS grid size
INTEGER :: i4time ! analysis time in seconds from
! 00:00 UTC Jan. 1, 1960
REAL :: time ! analysis time in seconds from the
! model initiali time
CHARACTER (LEN=*) :: dirname ! directory for data dump
CHARACTER (LEN=*) :: runname ! a string identify the run
!
! INPUT: model grid
!
REAL :: dx,dy ! ADAS grid spacing
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 :: t_3d(nx,ny,nz) ! temperature (K)
REAL :: p_3d(nx,ny,nz)
REAL :: rh_3d(nx,ny,nz)
REAL :: topo (nx,ny) ! The height of the terrain (equivalent
!
! INPUT for SAO obs.
!
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.
!
! INPUT: for radar
INTEGER :: istat_radar ! Valid radar reflectivity data?
REAL :: ref_mos_3d(nx,ny,nz) ! 3D radar reflectivity on model grid
! may be modified by sate VIS data
!
! OUTPUT:
INTEGER :: istatus ! Flag for process status
REAL :: clouds_3d(nx,ny,nz) ! Final cloud cover field
REAL :: cloud_ceiling (nx,ny) ! cloud ceiling (M AGL)
!
! LOCAL: for first guess fields
REAL :: cf_modelfg (nx,ny,nz) ! first guess cloud cover field
REAL :: rh_modelfg (nx,ny,nz) ! first guess relative humidity fld
REAL :: z_lcl (nx,ny) ! lifting condensation level(MSL)
!
! LOCAL: for surface arrays
REAL :: latgr (nx,ny) ! latitude at each grid point
REAL :: longr (nx,ny) ! longitude at each grid point
REAL :: pres_sfc_pa(nx,ny) ! sfc pressure (Pascal)
REAL :: t_gnd_k(nx,ny) ! ground skin temp.
REAL :: t_sfc_k (nx,ny) ! surface temperature field
REAL :: cvr_snow (nx,ny) ! surface snow coverage
REAL :: rland_frac (nx,ny) ! land coverage fraction
!
! LOCAL: for SAO cloud sounding
INTEGER :: ista_snd (max_cld_snd) ! seq. index of cld snd stn
INTEGER :: i_snd (max_cld_snd) ! i-loctn of cld snd
INTEGER :: j_snd (max_cld_snd) ! j-loctn of cld snd
INTEGER :: n_cld_snd ! # of cloud snds created
REAL :: cld_snd (max_cld_snd,nz) ! SAO cld snd
REAL :: wt_snd (max_cld_snd,nz) ! wgt for SAO cld snd
!
! LOCAL for gridded cloud cover analysis
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)
! LOCAL: for solar positions
REAL :: solar_alt(nx,ny) ! solar altitude angle
REAL :: solar_ha(nx,ny) ! solar hour angle
REAL :: solar_dec ! solar declination angle
!
! LOCAL: for satellite insertion
INTEGER :: io_vis ! Valid satellite VIS albedo data?
REAL :: albedo(nx,ny) ! vis. albedo derived from sat.vis.data
REAL :: cloud_frac_vis_a(nx,ny)! cldcvr derived from vis. albedo
INTEGER :: istat_vis ! status of inserting VIS data
INTEGER :: io_ir ! Valid satellite bright. temp. data?
REAL :: tb8_k(nx,ny) ! Obs. sate. band8 bright. temp.
REAL :: tb8_cold_k(nx,ny) ! Cold filtered band8 bright. temp.
REAL :: cldtop_m_tb8(nx,ny) ! Sat. cld top ht (m) from band 8 data
REAL :: cldtop_m (nx,ny) ! cldtop height
REAL :: cloud_top (nx,ny) ! cloud top (M ASL)
INTEGER :: istatus_ir ! status of inserting IR data
!
REAL :: sfc_sao_buffer ! No clearing cloud by sat. below
! this ht.(m AGL) (hence letting
! SAOs dominate)
PARAMETER (sfc_sao_buffer = 800.) ! meters, keep lower SAO cld from
! cleared out by sate. ir data
!
! LOCAL: for radar data insertion
REAL :: dbz_max_2d(nx,ny) ! column max. radar refl. (dbZ)
REAL :: cloud_base (nx,ny) ! cloud base heights from SAO+sate
REAL :: cloud_base_buf(nx,ny) ! lowest cloud base heights within
! a searching radius
LOGICAL :: l_unresolved(nx,ny) ! no SAO+sate cld, yet radar echo is
! above the threshold
REAL :: vis_rad_thdbz
PARAMETER (vis_rad_thdbz=10.)! threshold define cloudy refl.
REAL :: vis_rad_thcvr
PARAMETER (vis_rad_thcvr=0.2)! VIS cld cvr threshold,
! below this thresh. may cause
! radar echo being cleared out.
!
! LOCAL: intermidiate products
REAL :: cvr_max (nx,ny) ! column maximum cloud cover
REAL :: cvr_tot (nx,ny) ! final column total cloud cover
REAL :: cvr_1d (nz) ! cloud cover array for 1 grid col.
!
! LOCAL: for post-process
REAL :: thresh_cvr_base,thresh_cvr_top,thresh_cvr_ceiling
PARAMETER (thresh_cvr_base = 0.1)
PARAMETER (thresh_cvr_top = 0.1)
PARAMETER (thresh_cvr_ceiling = 0.65)
REAL :: default_top,default_base,default_ceiling
PARAMETER (default_base = 0.0)
PARAMETER (default_top = 0.0)
PARAMETER (default_ceiling = 0.0)
!
!-----------------------------------------------------------------------
!
! Misc local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: k1,l,j1,j2,ibeg
LOGICAL :: l_stn_clouds ! Using SAO stns' cloud obs?
PARAMETER (l_stn_clouds = .true.)
LOGICAL :: cldprt
!
REAL :: default_clear_cover
PARAMETER (default_clear_cover=0.01)
!
INTEGER :: istatus_sao,istatus_fg,istatus_barnes &
,istatus_rad,istatus_vis
INTEGER :: i,j,k,imid,jmid
REAL :: grid_spacing_m
REAL :: albmin,albmax,vcfmin,vcfmax
CHARACTER (LEN=6) :: varid
CHARACTER (LEN=20) :: varname
CHARACTER (LEN=20) :: varunits
INTEGER :: nstn,istn(mx_sng),jstn(mx_sng)
REAL :: hstn(mx_sng)
CHARACTER (LEN=5) :: cstn(mx_sng)
INTEGER :: istat_col
CHARACTER (LEN=6) :: satnam
CHARACTER (LEN=80) :: fname
REAL :: latsat,lonsat
INTEGER :: itime,isource,ifield,lrunnam,istatwrt
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
WRITE(6,*)' Welcome to the ADAS Cloud Analysis'
cldprt=(clddiag == 1)
ibeg=MAX(1,(nx-10))
istatus=0
imid=nx/2
jmid=ny/2
!
!-----------------------------------------------------------------------
!
! Initialize cloud cover fields on the analysis grid
!
!-----------------------------------------------------------------------
!
PRINT*,'Call background fields routine'
CALL inicldgrd
(nx,ny,nz,zs, &
bgqcopt,default_clear_cover, &
topo,t_3d,p_3d,rh_3d,cf_modelfg,t_sfc_k,pres_sfc_pa, &
cldcv,wtcldcv,z_ref_lcl,z_lcl,rh_modelfg, &
r_missing,istatus_fg)
!
IF (cldprt) THEN
WRITE(6,'(/a)') ' ==== CLOUD_CV cloud cover first guess===='
WRITE(6,632) (i,i=ibeg,nx)
WRITE(6,633) (k,(cf_modelfg(i,jmid,k),i=ibeg,nx) &
,k=nz,1,-1)
END IF
632 FORMAT(/1X,3X,11(1X,i4,1X))
633 FORMAT(1X,i3,11F6.1)
!
IF (istatus_fg /= 1)THEN
WRITE(6,*)' Error in background field: Aborting cloud analysis'
GO TO 999
END IF
!
IF (cloudopt == 2) THEN ! Using radar data only
PRINT*,'## NOTE ##: Using radar data only. LCL is used.'
GO TO 330
END IF
!
!-----------------------------------------------------------------------
!
! Create cloud sounding from the SAO cloud observations.
!
!-----------------------------------------------------------------------
!
WRITE(6,*)
WRITE(6,*)' Call Ingest/Insert SAO routines'
WRITE(6,*)
n_cld_snd = 0
IF (nobsng >= 1) THEN
!
CALL insert_sao1 (nx,ny,nz,dx,dy,xs,ys,zs,topo,t_sfc_k &
,cldcv,wtcldcv,t_3d,rh_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_sao)
!
IF(istatus_sao /= 1) THEN
WRITE(6,*)' Error inserting SAO data: Aborting cloud analysis'
GO TO 999
END IF
ELSE
WRITE(6,'(a,a,/a)') &
' !!!WARNING!!! No surface cloud obs. available.', &
' Skipping SAOs insertion'
END IF
!
!-----------------------------------------------------------------------
!
! Initialize clouds_3d array
!
!-----------------------------------------------------------------------
!
DO k=1,nz
DO j=1,ny
DO i=1,nx
clouds_3d(i,j,k)=0.0
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Do grid analysis on SAO data using Barnes interpolation.
!
!-----------------------------------------------------------------------
!
WRITE(6,*)
WRITE(6,*)' Analyzing SAO cloud soundings'
!
grid_spacing_m = dx
CALL barnes_r5 (nx,ny,nz,cldcv,clouds_3d,wtcldcv,cf_modelfg &
,l_stn_clouds,default_clear_cover,cld_snd,wt_snd &
,grid_spacing_m,i_snd,j_snd,n_cld_snd &
,istatus_barnes)
!
IF(cldprt) THEN
WRITE(6,'(/a)') ' ==== CLOUD_CV cloud cover barnes===='
WRITE(6,632) (i,i=ibeg,nx)
WRITE(6,633) (k,(clouds_3d(i,jmid,k),i=ibeg,nx) &
,k=nz,1,-1)
END IF
!632 format(/1x,3x,11(1x,i4,1x))
!633 format(1x,i3,11f6.1)
IF (istatus_barnes /= 1) THEN
WRITE(6,*)' Error in Barnes interp, Aborting cld analysis...'
GO TO 999
END IF
!
!-----------------------------------------------------------------------
!
! Find all stations within the model domain and print the cloud
! soundings at these station points after the Barnes interpolation.
!
!-----------------------------------------------------------------------
!
k=0
DO i=1,n_cld_snd
IF(i_snd(i) >= 1.AND.i_snd(i) <= nx &
.AND. j_snd(i) <= ny.AND.j_snd(i) >= 1) THEN
k=k+1
cstn(k)=stn(ista_snd(i))
istn(k)=i_snd(i)
jstn(k)=j_snd(i)
hstn(k)=0.001*elevsta(ista_snd(i))
END IF
END DO
nstn=k
!
IF (cld_files == 1) THEN
CALL gtlfnkey( runname, lrunnam )
WRITE(fname,'(a,a)') runname(1:lrunnam),'.cld_brns_snd'
OPEN (UNIT=14,FILE=trim(fname),STATUS='unknown')
WRITE(14,404) nobsng,nstn,n_cld_snd
404 FORMAT(1X,' Total/in-domain # of stations:',2I5/ &
' Cloud observations used for deriving cloud soundings:',i3/)
i = nstn/24
IF(i < 1) THEN
WRITE(14,405) (cstn(k),k=1,nstn)
405 FORMAT(1X,' hgt(m)',1X,24A5)
WRITE(14,407) (istn(k),k=1,nstn)
WRITE(14,407) (jstn(k),k=1,nstn)
407 FORMAT(1X,7X,24I5)
WRITE(14,408) (hstn(k),k=1,nstn)
408 FORMAT(1X,7X,24F5.2)
DO k1=1,nz
WRITE(14,406) zs(istn(1),jstn(1),k1) &
,(clouds_3d(istn(l),jstn(l),k1),l=1,nstn)
406 FORMAT(1X,f7.1,24F5.2)
END DO !k1
ELSE
DO k=1,i
j1=(k-1)*24+1
j2=k*24
WRITE(14,405) (cstn(l),l=j1,j2)
WRITE(14,407) (istn(l),l=j1,j2)
WRITE(14,407) (jstn(l),l=j1,j2)
WRITE(14,408) (hstn(l),l=j1,j2)
DO k1=1,nz
WRITE(14,'(1x,f7.1,24f5.2)') zs(istn(j1),jstn(j1),k1) &
,(clouds_3d(istn(l),jstn(l),k1),l=j1,j2)
END DO !k1
END DO !k
j=i*24+1
WRITE(14,405) (cstn(l),l=j,nstn)
WRITE(14,407) (istn(l),l=j,nstn)
WRITE(14,407) (jstn(l),l=j,nstn)
WRITE(14,408) (hstn(l),l=j,nstn)
DO k1=1,nz
WRITE(14,406) zs(istn(j),jstn(j),k1) &
,(clouds_3d(istn(l),jstn(l),k1),l=j,nstn)
END DO !k1
END IF
CLOSE(14)
END IF
!
IF (cld_files == 1) THEN
WRITE(6,*) 'Writing 3D cloud fraction after inserting SAOs.'
varid='cldsao'
varname='Clouds after SAO'
varunits='Fraction'
CALL wrtvar1(nx,ny,nz,clouds_3d,varid,varname,varunits, &
time,runname,dirname,istatwrt)
END IF
!
!-----------------------------------------------------------------------
!
! Read in the surface snow coverage data.
! Currently, the snow cover field is simply set to zero
!
!-----------------------------------------------------------------------
!
! varname='snwcvr'
! CALL READVAR(nx,ny,1, varname,time,cvr_snow, runname)
DO j = 1,ny
DO i = 1,nx
cvr_snow(i,j) = 0.0
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Read in the land (water) coverage fraction data.
! Currently, this is simply set to 1.
!
!-----------------------------------------------------------------------
!
! varname='lndfrc'
! CALL READVAR(nx,ny,1, varname,time,rland_frac, runname)
DO j = 1,ny
DO i = 1,nx
rland_frac(i,j) = 1.0
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Calculate solar altitude
!
!-----------------------------------------------------------------------
!
WRITE(6,*) 'Calculating solar position parameters.'
!
DO j = 1,ny
DO i = 1,nx
CALL solar_pos(latgr(i,j),longr(i,j) &
,i4time,solar_alt(i,j),solar_dec,solar_ha(i,j))
END DO
END DO
imid = nx/2
jmid = ny/2
PRINT*,' pt(',imid,',',jmid,'): Solar alt=',solar_alt(imid,jmid), &
' declin=',solar_dec,' hour_angle=',solar_ha(imid,jmid)
!
!-----------------------------------------------------------------------
!
! READ IN SATELLITE VISIBLE DATA
!
!-----------------------------------------------------------------------
!
ifield=1
isource=1
WRITE(6,*) 'Getting satellite visible data remapped on the', &
' model grid'
!
CALL rdsatfld(nx,ny,1,vis_fname,satnam,latsat,lonsat, &
itime,isource,varname,albedo,io_vis)
!
IF (cld_files == 1) THEN
varid='albedo'
varname='Albedo'
varunits='Fraction'
CALL wrtvar1(nx,ny,1 ,albedo, varid,varname,varunits, &
time,runname,dirname,istatwrt)
END IF
IF(io_vis == 0) THEN
WRITE(6,'(a)') ' Read successful, calling get_vis'
CALL get_vis (solar_alt,nx,ny &
,cloud_frac_vis_a,albedo,r_missing,istat_vis)
albmax = -999.0
albmin = 999.0
vcfmax = -999.0
vcfmin = 999.0
DO j=1,ny
DO i=1,nx
IF(albedo(i,j) >= 0.) THEN
albmax=MAX(albmax,albedo(i,j))
albmin=MIN(albmin,albedo(i,j))
END IF
IF(cloud_frac_vis_a(i,j) >= 0.) THEN
vcfmax=MAX(vcfmax,cloud_frac_vis_a(i,j))
vcfmin=MIN(vcfmin,cloud_frac_vis_a(i,j))
END IF
END DO
END DO
WRITE(6,'(a,a)') 'Visible Satellite Data: ',vis_fname
WRITE(6,'(a,f10.4,a,f10.4)') &
' min albedo:',albmin,' max albedo:',albmax
WRITE(6,'(a,f10.4,a,f10.4)') &
' min vis cf:',vcfmin,' max vis cf:',vcfmax
ELSE
WRITE(6,'(a,a,/a)') ' ###Problem reading ',vis_fname, &
' Skipping vis processing'
END IF
!
IF (cld_files == 1) THEN
varid ='cldalb'
varname='Cloud from Albedo'
varunits='Fraction'
CALL wrtvar1(nx,ny,1,cloud_frac_vis_a,varid,varname,varunits, &
time,runname,dirname,istatwrt)
END IF
!
!-----------------------------------------------------------------------
!
! Insert satellite IR data
!
!-----------------------------------------------------------------------
!
ifield=1
isource=1
WRITE(6,*) 'Inserting satellite IR data.'
WRITE(6,*) 'Getting satellite IR data remapped on the model', &
' grid'
!
CALL rdsatfld(nx,ny,1,ir_fname,satnam,latsat,lonsat, &
itime,isource,varname,tb8_k,io_ir)
!
IF (cld_files == 1) THEN
varid='cttemp'
varname='Cloud Top Temperature'
varunits='Degrees K'
CALL wrtvar1(nx,ny,1, tb8_k,varid,varname,varunits, &
time,runname,dirname,istatwrt)
END IF
!
IF(io_ir == 0) THEN
WRITE(6,'(a)') ' Read successful, calling insert_ir'
CALL insert_ir (nx,ny,nz,latgr,longr,rland_frac,cvr_snow, &
topo,zs,z_lcl,t_3d,p_3d, &
t_sfc_k,t_gnd_k,clouds_3d, &
solar_alt,solar_ha,solar_dec, &
tb8_k,tb8_cold_k,cldtop_m_tb8,cldtop_m, &
grid_spacing_m,sfc_sao_buffer, &
clouds_3d,istatus_ir)
ELSE
WRITE(6,'(a,a,/a)') &
' !!!WARNING!!! Problem using satellite IR data', &
' Skipping IR insertion'
END IF
!
IF (cld_files == 1) THEN
WRITE(6,*) 'Writing cloud fraction after inserting IR data.'
varid='cld_ir'
varname='Clouds after IR Data'
varunits='Fraction'
CALL wrtvar1(nx,ny,nz,clouds_3d,varid,varname,varunits, &
time,runname,dirname,istatwrt)
END IF
!
!-----------------------------------------------------------------------
!
! Getting the 2-D maximum radar reflectivity field
!
!-----------------------------------------------------------------------
!
DO i=1,nx
DO j=1,ny
dbz_max_2d(i,j) = r_missing
DO k=1,nz
dbz_max_2d(i,j) = MAX(dbz_max_2d(i,j),ref_mos_3d(i,j,k))
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Insert radar data
!
!-----------------------------------------------------------------------
!
330 CONTINUE
IF (istat_radar == 1) THEN
WRITE(6,*) 'Inserting radar reflectivity data.'
CALL insert_radar(nx,ny,nz,clddiag,topo,zs,t_3d,z_lcl &
,refthr1,refthr2,hgtrefthr &
,ref_mos_3d,clouds_3d &
,cloud_base,cloud_base_buf,l_unresolved &
,istatus_rad)
!
IF (istatus_rad /= 1) THEN
WRITE(6,*)' Error inserting radar data,' &
,' Aborting cld analysis'
GO TO 999
END IF
ELSE
WRITE(6,'(a,a,/a)') ' !!!WARNING!!! Problem using radar data', &
' Skipping RADAR insertion'
END IF
!
IF (cld_files == 1) THEN
WRITE(6,*) 'Writing cloud fraction after inserting radar data.'
varid='cldrad'
varname='Clouds after Radar'
varunits='Fraction'
CALL wrtvar1(nx,ny,nz,clouds_3d,varid,varname,varunits, &
time,runname,dirname,istatwrt)
END IF
!
IF (cloudopt == 2) GO TO 340 ! Using radar data only
!
!-----------------------------------------------------------------------
!
! Inserting the satellite visible data
!
!-----------------------------------------------------------------------
!
IF (io_vis == 0) THEN
WRITE(6,*) 'Inserting satellite VIS data.'
CALL insert_vis (nx,ny,nz,zs,topo,clouds_3d &
,albedo,cloud_frac_vis_a &
,vis_rad_thcvr,vis_rad_thdbz &
,istat_radar,ref_mos_3d,refthr2,dbz_max_2d &
,r_missing,sfc_sao_buffer,istatus_vis)
!
IF (istatus_vis /= 1) THEN
WRITE(6,*)' Error inserting VIS data, Aborting cld analysis.'
GO TO 999
END IF
ELSE
WRITE(6,'(a,a,/a)') &
' !!!WARNING!!! Problem using satellite VIS data', &
' Skipping VIS data insertion'
END IF
!
IF (cld_files == 1) THEN
WRITE(6,*) 'Writing cloud fraction after inserting VIS data.'
varid='cldvis'
varname='Clouds after Vis'
varunits='Fraction'
CALL wrtvar1(nx,ny,nz,clouds_3d,varid,varname,varunits, &
time,runname,dirname,istatwrt)
END IF
!
340 CONTINUE
!
!-----------------------------------------------------------------------
!
! Calculate the column maximum and column total cloud cover.
!
!-----------------------------------------------------------------------
!
DO i=1,nx
DO j=1,ny
DO k=1,nz
cvr_1d (k) = clouds_3d(i,j,k)
END DO
CALL col_max_tot (nz,cvr_1d,cvr_max(i,j),cvr_tot(i,j) &
,istat_col)
END DO
END DO
IF (cld_files == 1) THEN
WRITE(6,*) 'Writing column total cloud cover.'
varid='totcvr'
varname='Total Column Cloud Cover'
varunits='Fraction'
CALL wrtvar1(nx,ny,1,cvr_tot,varid,varname,varunits, &
time,runname,dirname,istatwrt)
END IF
!
!-----------------------------------------------------------------------
!
! Compare the final cloud analysis with radar reflectivity
!
!-----------------------------------------------------------------------
!
IF (istat_radar == 1) THEN
CALL compare_radar(nx,ny,nz,ref_mos_3d,dbz_max_2d &
,cvr_max,refthr2,cloud_frac_vis_a &
,vis_rad_thcvr,vis_rad_thdbz,r_missing)
END IF
!
!-----------------------------------------------------------------------
!
! Get Cloud Bases and Tops
!
!-----------------------------------------------------------------------
!
WRITE(6,*)' Calculating Cloud Ceiling, Base and Top'
!
DO j = 1,ny-1
DO i = 1,nx-1
cloud_base(i,j) = default_base
cloud_ceiling(i,j) = default_ceiling
cloud_top(i,j) = default_top
!
!-----------------------------------------------------------------------
!
! Test for Cloud Base (MSL)
!
!-----------------------------------------------------------------------
!
IF (cvr_max(i,j) >= thresh_cvr_base) THEN
DO k = nz-1,1,-1
IF(clouds_3d(i,j,k) < thresh_cvr_base .AND. &
clouds_3d(i,j,k+1) >= thresh_cvr_base) THEN
cloud_base(i,j) = 0.5 * (zs(i,j,k) + zs(i,j,k+1))
END IF
END DO ! k
END IF ! Clouds exist in this column
!
!-----------------------------------------------------------------------
!
! Test for Cloud Top (MSL)
!
!-----------------------------------------------------------------------
!
IF (cvr_max(i,j) >= thresh_cvr_top) THEN
DO k = 1,nz-1
IF(clouds_3d(i,j,k) > thresh_cvr_top .AND. &
clouds_3d(i,j,k+1) <= thresh_cvr_top) THEN
cloud_top(i,j) = 0.5 * (zs(i,j,k) + zs(i,j,k+1))
END IF
END DO ! k
END IF ! Clouds exist in this column
!
!-----------------------------------------------------------------------
!
! Test for Cloud Ceiling (AGL)
!
!-----------------------------------------------------------------------
!
IF (cvr_max(i,j) >= thresh_cvr_ceiling) THEN
DO k = nz-1,1,-1
IF(clouds_3d(i,j,k) < thresh_cvr_ceiling .AND. &
clouds_3d(i,j,k+1) >= thresh_cvr_ceiling) THEN
cloud_ceiling(i,j) = 0.5 * (zs(i,j,k)+zs(i,j,k+1)) &
- topo(i,j)
END IF
END DO ! k
END IF ! Clouds exist in this column
END DO ! i
END DO ! j
!
!-----------------------------------------------------------------------
!
! Set the boundary values
!
!-----------------------------------------------------------------------
!
DO i = 1,nx-1
cloud_base(i,ny) = cloud_base(i,ny-1)
cloud_ceiling(i,ny) = cloud_ceiling(i,ny-1)
cloud_top(i,ny) = cloud_top(i,ny-1)
END DO ! i
!
DO j = 1,ny
cloud_base(nx,j) = cloud_base(nx-1,j)
cloud_ceiling(nx,j) = cloud_ceiling(nx-1,j)
cloud_top(nx,j) = cloud_top(nx-1,j)
END DO ! i
!
IF (cld_files == 1) THEN
varid='cldbas'
varname='Cloudbase'
varunits='Meters'
CALL wrtvar1(nx,ny,1,cloud_base,varid,varname,varunits, &
time,runname,dirname,istatwrt)
varid='cldtop'
varname='Cloud Top Height'
varunits='Meters'
CALL wrtvar1(nx,ny,1,cloud_top,varid,varname,varunits, &
time,runname,dirname,istatwrt)
varname='ceilng'
varname='Cloud Ceiling'
varunits='Meters'
CALL wrtvar1(nx,ny,1,cloud_ceiling,varid,varname,varunits, &
time,runname,dirname,istatwrt)
END IF
!
!-----------------------------------------------------------------------
istatus=1
999 CONTINUE
!
!-----------------------------------------------------------------------
!
RETURN
END SUBROUTINE cloud_cv
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE COL_MAX_TOT ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE col_max_tot (nz,cvr,cvr_max,cvr_tot,istatus) 1
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! To find the number of cloud decks in one column of grid.
!
!-----------------------------------------------------------------------
!
! AUTHOR: (Jian Zhang)
! 06/01/96.
!
! MODIFICATION HISTORY:
! 04/11/97 J. Zhang
! Included adascld24.inc for ncloud
! 08/06/97 J. Zhang
! Change adascld24.inc to adascld26.inc
! 05/05/98 J. Zhang
! Abandon cloud grid. Using only ARPS grid.
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! INPUT:
INTEGER :: nz
!
! INPUT:
REAL :: cvr(nz) ! Cloud cover in one column of grid.
!
! OUTPUT:
INTEGER :: istatus
REAL :: cvr_max ! column maximum cldcvr
REAL :: cvr_tot ! column integrated cldcvr
!
! LOCAL:
INTEGER :: nlyr ! number of cloud decks (layers) in the column.
INTEGER :: ilyr(nz) ! Layer index for each cloud lvl (needs)
REAL :: a(nz) ! max. cloud cover of each cloud layer
REAL :: f(nz) ! Apparnt "x-sectn" of cldlyrs seen from above
INTEGER :: ik(nz) ! Height level representative of cloud layers
!
!-----------------------------------------------------------------------
!
! Misc local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: k,n
REAL :: cvr_thresh,sumf
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.
!
!-----------------------------------------------------------------------
!
istatus=0
nlyr = 0
IF(cvr(nz) > cvr_thresh) THEN ! cld_top at upper bndry
nlyr = nlyr + 1
a(nlyr) = cvr(nz)
ik(nlyr) = nz
END IF
DO k = nz-1,1,-1
IF(cvr(k) >= cvr_thresh .AND. cvr(k+1) < cvr_thresh) THEN
nlyr = nlyr + 1 ! at top of a new cld layer
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
IF(cvr(k) >= cvr_thresh) THEN ! Still within layer
ilyr(k) = nlyr
ELSE ! Below layer
ilyr(k) = 0
END IF
END DO ! k
!
!-----------------------------------------------------------------------
!
! Convert cloud layer fractions to "cross-section" seen from
! satellite. This solves for the "f" array given the "a" array
!
!-----------------------------------------------------------------------
!
IF(nlyr == 0) THEN
istatus=1
cvr_max=0.0
cvr_tot=0.0
RETURN
END IF
a(1) = MIN( f(1) = a(1)
sumf = f(1)
cvr_max=a(1)
cvr_tot=a(1)
IF(nlyr >= 2) THEN
DO n = 2,nlyr
a(n) = MIN( cvr_max = MAX( 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
!
!-----------------------------------------------------------------------
!
! Calculate column integrated cloud cover as if the cloud
! particles are uniformly (fusely) distributed in the atmosphere.
!
!-----------------------------------------------------------------------
!
cvr_tot =sumf
!
END IF ! nlyr
!
!-----------------------------------------------------------------------
!
istatus=1
RETURN
END SUBROUTINE col_max_tot
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE COMPARE_RADAR ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE compare_radar (nx,ny,nz,ref_mos_3d,dbz_max_2d & 1
,cvr_max,ref_base,cf_vis_a &
,vis_rad_thcvr,vis_rad_thdbz,r_missing)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
! This routine compares the cloud and radar fields and flags
! remaining differences that weren't caught in earlier processing
! This routine does not alter any arrays or pass anything back,
! it is diagnostic only.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Jian Zhang
! 07/15/96 Based on the LAPS cloud analysis code of 9/1995
!
! MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
! INPUT:
INTEGER :: nx,ny,nz ! Model grid size(dbz)
REAL :: ref_mos_3d(nx,ny,nz) ! radar reflectivity remapped on grid
REAL :: dbz_max_2d(nx,ny) ! Column max. radar reflectivity (dbz)
REAL :: cvr_max(nx,ny) ! column maximum cloud cover
REAL :: cf_vis_a(nx,ny) ! sate. visible albedo derived cldcvr
REAL :: ref_base ! thrshld define signif. radar echo
REAL :: vis_rad_thcvr ! thrshld define signif. vis cldcvr (0.2)
REAL :: vis_rad_thdbz ! thrshld define signif. cldy radar ehco (10dbz)
REAL :: r_missing ! flag for missing or bad data
!
! OUTPUT:
! None
!
!-----------------------------------------------------------------------
!
! Misc local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j
!
!-----------------------------------------------------------------------
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
WRITE(6,*) 'Comparing clouds and radar (_RDR)'
WRITE(6,*) 'vis_rad_thcvr = ',vis_rad_thcvr
WRITE(6,*) 'vis_rad_thdbz = ',vis_rad_thdbz
DO i = 1,nx
DO j = 1,ny
IF(cvr_max(i,j) < vis_rad_thcvr .AND. dbz_max_2d(i,j) > ref_base) THEN
!
!-----------------------------------------------------------------------
!
! We have a discrepancy between the VIS and radar
!
!-----------------------------------------------------------------------
!
IF(cvr_max(i,j) == cf_vis_a(i,j)) THEN ! CVR = VIS
IF(dbz_max_2d(i,j) < vis_rad_thdbz) THEN
WRITE(6,1) i,j,cvr_max(i,j),dbz_max_2d(i,j),cf_vis_a(i,j)
1 FORMAT(' VIS_RDR: cvr_mx/dbz/vis <',2I4,f5.2,f8.1,f5.2)
ELSE
WRITE(6,11) i,j,cvr_max(i,j),dbz_max_2d(i,j),cf_vis_a(i,j)
11 FORMAT(' VIS_RDR: cvr_mx/dbz/vis >',2I4,f5.2,f8.1,f5.2)
END IF ! dbz_max_2d < vis_rad_thdbz
ELSE IF (cvr_max(i,j) < cf_vis_a(i,j)) THEN
!
!-----------------------------------------------------------------------
!
! Don't blame the VIS CVR < VIS
!
!-----------------------------------------------------------------------
!
IF(dbz_max_2d(i,j) < vis_rad_thdbz) THEN
!
!-----------------------------------------------------------------------
!
! We can potentially block out the radar
!
!-----------------------------------------------------------------------
!
WRITE(6,2) i,j,cvr_max(i,j),dbz_max_2d(i,j),cf_vis_a(i,j)
2 FORMAT(' CLD_RDR: cvr_mx/dbz/vis <',2I4,f5.2,f8.1,f5.2)
ELSE ! Radar is too strong to block out
WRITE(6,3) i,j,cvr_max(i,j),dbz_max_2d(i,j),cf_vis_a(i,j)
3 FORMAT(' CLD_RDR: cvr_mx/dbz/vis >',2I4,f5.2,f8.1,f5.2)
END IF ! dbz_max_2d < vis_rad_thdbz
ELSE IF (cvr_max(i,j) > cf_vis_a(i,j)) THEN
!
!-----------------------------------------------------------------------
!
! Don't know if VIS lowered cloud cover CVR > VIS
! At least if difference is less than VIS "cushion"
!
!-----------------------------------------------------------------------
!
IF(dbz_max_2d(i,j) < vis_rad_thdbz) THEN
!
!-----------------------------------------------------------------------
!
! We can potentially block out the radar
!
!-----------------------------------------------------------------------
!
WRITE(6,4) i,j,cvr_max(i,j),dbz_max_2d(i,j),cf_vis_a(i,j)
4 FORMAT(' ???_RDR: cvr_mx/dbz/vis <',2I4,f5.2,f8.1,f5.2)
ELSE ! Radar is too strong to block out
WRITE(6,5) i,j,cvr_max(i,j),dbz_max_2d(i,j),cf_vis_a(i,j)
5 FORMAT(' ???_RDR: cvr_mx/dbz/vis >',2I4,f5.2,f8.1,f5.2)
END IF ! dbz_max_2d < vis_rad_thdbz
END IF ! cvr_max(i,j) = cf_vis_a(i,j)
END IF ! cvr_max < vis_rad_thcvr & dbz_max_2d > ref_base
END DO ! j
END DO ! i
RETURN
END SUBROUTINE compare_radar