subroutine gsdcloudanalysis4gfs(mype)
!
!$$$ subprogram documentation block
! . . . .
! subprogram: gsdcloudanalysis driver for generalized cloud/hydrometeor analysis
!
! PRGMMR: Ming Hu ORG: GSD/AMB DATE: 2006-10-27
!
! ABSTRACT:
! This subroutine serves as a driver for generalized cloud/hydrometeor analysis
!
! PROGRAM HISTORY LOG:
! 2008-12-20 Hu Add NCO document block
! 2010-04-30 Hu Clean the code to meet GSI standard
! 2011-05-29 Todling - extra cloud-guess from MetGuess-Bundle (see Remark 1)
! some fields now from wrf_mass_guess_mod
! 2013-10-19 todling - metguess now holds background
! 2013-10-24 todling - revisit strip interface
! 2014-10-22 Hu - Add analysis for rain number concentation
! reflectivity between 15-28dBZ
! 2014-12-22 Hu - Add light rain in precipiation analysis using radar
! reflectivity between 15-28dBZ
!
!
! input argument list:
! mype - processor ID that does this IO
!
! output argument list:
!
! USAGE:
! INPUT FILES:
!
! OUTPUT FILES:
!
! REMARKS:
! 1. Notice that now the fields point to the instance of the GUESS at ntguessig
! no longer wired to 1 (as originally) - however, make sure when running RUC
! ntguessig is set correctly (see questions around definition of itsig)
! 2. Notice that some WRF-variable and grid specific is now defined in wrf_guess_mod
!
! ATTRIBUTES:
! LANGUAGE: FORTRAN 90
! MACHINE: Linux cluster (WJET) at NOAA/ESRL - Boulder, CO
!
!$$$
!
!_____________________________________________________________________
!
!
use constants, only: zero,one,rad2deg,fv
use constants, only: rd_over_cp, h1000
use kinds, only: r_single,i_kind, r_kind
use gridmod, only: pt_ll,eta1_ll,aeta1_ll
use gridmod, only: regional,wrf_mass_regional,regional_time
use gridmod, only: nsig,lat2,lon2,istart,jstart,twodvar_regional
use gridmod, only: itotsub,lon1,lat1,nlon,nlat,ijn,displs_g,strip
use gridmod, only: rlats,rlons
use guess_grids, only: pbl_height, load_gsdpbl_hgt
use obsmod, only: obs_setup,nsat1,ndat,dtype
use guess_grids, only: ntguessig,ntguessfc
use wrf_mass_guess_mod, only: soil_temp_cld,isli_cld,ges_xlon,ges_xlat,ges_tten
use guess_grids, only: isli,soil_temp,isli2
use guess_grids, only: ges_tsen,ges_prsl
use jfunc, only: tsensible
use mpimod, only: mpi_comm_world,ierror,mpi_real4
use rapidrefresh_cldsurf_mod, only: dfi_radar_latent_heat_time_period, &
metar_impact_radius, &
metar_impact_radius_lowCloud, &
l_cleanSnow_WarmTs,l_conserve_thetaV,&
r_cleanSnow_WarmTs_threshold, &
i_conserve_thetaV_iternum, &
l_cld_bld, cld_bld_hgt, &
build_cloud_frac_p, clear_cloud_frac_p, &
nesdis_npts_rad, &
iclean_hydro_withRef, iclean_hydro_withRef_allcol, &
i_lightpcp,i_gsdcldanal_type,ioption
use gsi_metguess_mod, only: GSI_MetGuess_Bundle
use gsi_bundlemod, only: gsi_bundlegetpointer
#ifdef RR_CLOUDANALYSIS
implicit none
! Declare passed variables
integer(i_kind),intent(in):: mype
!
! background
!
real(r_single),allocatable:: t_bk(:,:,:)
real(r_single),allocatable:: h_bk(:,:,:)
real(r_single),allocatable:: p_bk(:,:,:)
real(r_single),allocatable:: ps_bk(:,:)
real(r_single),allocatable:: zh(:,:)
real(r_single),allocatable:: q_bk(:,:,:)
real(r_single),allocatable:: xlon(:,:) ! 2D longitude in each grid
real(r_single),allocatable:: xlat(:,:) ! 2D latitude in each grid
real(r_single),allocatable:: gsfc(:,:,:)
real(r_single), allocatable:: xland(:,:)
real(r_single),allocatable:: soiltbk(:,:)
! real(r_single),allocatable:: z_lcl(:,:) ! lifting condensation level
! real(r_single),allocatable:: pblh(:,:) ! PBL height (grid coordinate)
!
! surface observation
!
integer(i_kind) :: nvarcld_p
parameter (nvarcld_p=13)
integer(i_kind) :: numsao
real(r_single), allocatable :: oi(:)
real(r_single), allocatable :: oj(:)
integer(i_kind),allocatable :: ocld(:,:)
character*10, allocatable :: owx(:)
real(r_single), allocatable :: oelvtn(:)
real(r_single), allocatable :: odist(:)
character(8), allocatable :: cstation(:)
real(r_single), allocatable :: oistation(:)
real(r_single), allocatable :: ojstation(:)
real(r_single), allocatable :: wimaxstation(:)
!
integer(i_kind),allocatable :: osfc_station_map(:,:)
!
! lightning observation: 2D field in RR grid
!
real(r_single),allocatable :: lightning(:,:)
!
! GOES - NASA LaRc cloud products: several 2D fields in RR grid
!
real(r_single),allocatable :: nasalarc_cld(:,:,:)
!
! radar observation : 3D reflectvity in RR grid
!
real(r_kind),allocatable :: ref_mos_3d(:,:,:)
real(r_kind),allocatable :: ref_mos_3d_tten(:,:,:)
real(r_kind),allocatable :: ref_mosaic31(:,:,:)
integer(i_kind) :: nmsclvl_radar
!
! GOES - NESDIS cloud products : 2d fields
!
real(r_single), allocatable :: sat_ctp(:,:)
real(r_single), allocatable :: sat_tem(:,:)
real(r_single), allocatable :: w_frac(:,:)
integer(i_kind),allocatable :: nlev_cld(:,:)
!
! cloud/hydrometeor analysis variables
!
!=========================================================================
! cld_cover_3d in the Generalized Cloud/Hydrometeor Analysis code
! Definition: 3-d gridded observation-based information
! including 0-1 cloud-fraction (with zero value for clear)
! and negative values indicating "unknown" status.
! cld_cover_3d is initialized with negative values.
! cld_type_3d, pcp_type_3d, wthr_type_2d - similar to cld_cover_3d
!=========================================================================
real(r_single), allocatable :: cld_cover_3d(:,:,:) ! cloud cover
integer(i_kind),allocatable :: cld_type_3d(:,:,:) ! cloud type
integer(i_kind),allocatable :: pcp_type_3d(:,:,:) ! precipitation type
integer(i_kind),allocatable :: wthr_type_2d(:,:) ! weather type type
integer(i_kind),allocatable :: cloudlayers_i(:,:,:) ! 5 different layers
! 1= the number of layers
! 2,4,... bottom
! 3,5,... top
!
real(r_single),allocatable :: cldwater_3d(:,:,:) ! cloud water
real(r_single),allocatable :: cldice_3d(:,:,:) ! cloud ice
real(r_single),allocatable :: rain_3d(:,:,:) ! rain
real(r_single),allocatable :: nrain_3d(:,:,:) ! rain number concentration
real(r_single),allocatable :: snow_3d(:,:,:) ! snow
real(r_single),allocatable :: graupel_3d(:,:,:) ! graupel
real(r_single),allocatable :: cldtmp_3d(:,:,:) ! cloud temperature
real(r_single),allocatable :: rain_1d_save(:) ! rain
real(r_single),allocatable :: nrain_1d_save(:) ! rain number concentration
real(r_single),allocatable :: snow_1d_save(:) ! snow
real(r_single),allocatable :: vis2qc(:,:) ! fog
real(r_kind) :: thunderRadius=2.5_r_kind
real(r_single) :: r_radius ! influence radius of cloud based on METAR obs
real(r_single) :: r_radius_lowCloud ! influence radius of low cloud to cloud top pressure
integer(i_kind) :: miss_obs_int
real(r_kind) :: miss_obs_real
parameter ( miss_obs_int = -99999 )
parameter ( miss_obs_real = -99999.0_r_kind )
real(r_single) :: krad_bot ! radar bottom level
!
! collect cloud
real(r_kind) :: cloud_def_p
data cloud_def_p / 0.000001_r_kind/
real(r_kind),allocatable :: sumqci(:,:,:) ! total liquid water
real(r_kind),allocatable :: watericemax(:,:) ! max of total liquid water
integer(i_kind),allocatable :: kwatericemax(:,:) ! lowest level of total liquid water
real(r_single),allocatable::temp1(:,:),tempa(:)
real(r_single),allocatable::all_loc(:,:)
real(r_single),allocatable::strp(:)
integer(i_kind) :: im,jm
!
! option in namelist
!
integer(i_kind) :: opt_cloudwaterice_retri ! method for cloud water retrieval
integer(i_kind) :: opt_hydrometeor_retri ! method for precipitation retrieval
integer(i_kind) :: opt_cloudtemperature ! if open temperature adjustment scheme
integer(i_kind) :: istat_surface,istat_nesdis,istat_radar ! 1 has observation
integer(i_kind) :: istat_nasalarc,istat_lightning ! 0 no observation
integer(i_kind) :: imerge_nesdis_nasalarc ! =1 merge NASA LaRC with NESDIS
! =2 use NASA LaRC only
! = other, use NESDIS only
!
!
real(r_kind), pointer :: ges_z (:,: )=>NULL() ! geopotential height
real(r_kind), pointer :: ges_ps(:,: )=>NULL() ! surface pressure
real(r_kind), pointer :: ges_tv(:,:,:)=>NULL() ! virtual temperature
real(r_kind), pointer :: ges_q (:,:,:)=>NULL() ! specifici humidity
real(r_kind), pointer :: ges_ql(:,:,:)=>NULL() ! cloud water
real(r_kind), pointer :: ges_qi(:,:,:)=>NULL() ! could ice
real(r_kind), pointer :: ges_qr(:,:,:)=>NULL() ! rain
real(r_kind), pointer :: ges_qs(:,:,:)=>NULL() ! snow
real(r_kind), pointer :: ges_qg(:,:,:)=>NULL() ! graupel
real(r_kind), pointer :: ges_qnr(:,:,:)=>NULL() ! rain number concentration
real(r_kind), pointer,dimension(:,:,:)::ges_cwmr=>NULL()
!
! misc.
!
logical :: l_use_hydroretrieval_all
integer(i_kind) :: i,j,k,itsig,itsfc
integer(i_kind) :: iglobal,jglobal,ilocal,jlocal
logical :: ifindomain
integer(i_kind) :: imaxlvl_ref
real(r_kind) :: max_retrieved_qrqs,max_bk_qrqs,ratio_hyd_bk2obs
real(r_kind) :: qrlimit,qrlimit_lightpcp
character(10) :: obstype
integer(i_kind) :: lunin, is, ier, istatus
integer(i_kind) :: nreal,nchanl,ilat1s,ilon1s
character(20) :: isis
real(r_kind) :: refmax,snowtemp,raintemp,nraintemp,graupeltemp
real(r_kind) :: snowadd,ratio2
integer(i_kind) :: imax, jmax, ista, iob, job
real(r_kind) :: dfi_lhtp, qmixr, tsfc
integer(i_kind) :: ntsig,ntsfc
integer(i_kind) :: ii,jj
!
!
if(mype==0) then
write(6,*) '========================================'
write(6,*) 'gsdcloudanalysis4gfs: Start generalized cloud analysis:'
write(6,*) '========================================'
endif
!
if(i_gsdcldanal_type==30) then
ntsig = ntguessig
ntsfc = ntguessfc
ier=0
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(ntsig),'ps',ges_ps,istatus);ier=ier+istatus
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(ntsig),'z' ,ges_z ,istatus);ier=ier+istatus
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(ntsig),'tv',ges_tv,istatus);ier=ier+istatus
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(ntsig),'q', ges_q ,istatus);ier=ier+istatus
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(ntsig),'cw', ges_cwmr,istatus);ier=ier+istatus
! sfct(:,:,ntsig)=ges_tsen(:,:,1,ntsig)
soil_temp(:,:,ntsig)=ges_tsen(:,:,1,ntsig)
ges_ql=>ges_cwmr
else
itsig=1 ! _RT shouldn't this be ntguessig?
itsfc=1 ! _RT shouldn't this be ntguessig?
ier=0
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(itsfc),'ps',ges_ps,istatus);ier=ier+istatus
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(itsfc),'z' ,ges_z ,istatus);ier=ier+istatus
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(itsig),'tv',ges_tv,istatus);ier=ier+istatus
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(itsig),'q', ges_q ,istatus);ier=ier+istatus
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(itsig),'ql',ges_ql,istatus);ier=ier+istatus
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(itsig),'qi',ges_qi,istatus);ier=ier+istatus
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(itsig),'qr',ges_qr,istatus);ier=ier+istatus
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(itsig),'qs',ges_qs,istatus);ier=ier+istatus
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(itsig),'qg',ges_qg,istatus);ier=ier+istatus
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(itsig),'qnr',ges_qnr,istatus);ier=ier+istatus
endif
if(ier/=0) return ! no guess, nothing to do
!mhu if(mype<10) ges_cwmr=ges_cwmr+0.001
!
! if(mype==10) then
! write(*,*) lat2,lon2,nsig
! write(*,*) 'ps=',maxval(ges_ps(:,:)),minval(ges_ps(:,:))
! do k=1,nsig
! write(*,*) 't=',k,maxval(ges_tv(:,:,k)),minval(ges_tv(:,:,k))
! enddo
! do k=1,nsig
! write(*,*) 'q=',k,maxval(ges_q(:,:,k)),minval(ges_q(:,:,k))
! enddo
! do k=1,nsig
! write(*,*) 'cw=',k,maxval(ges_cwmr(:,:,k)),minval(ges_cwmr(:,:,k))
! enddo
! write(*,*) 'z=',maxval(ges_z(:,:)),minval(ges_z(:,:))
! do k=1,nsig
! write(*,*) 'tsen=',k,maxval(ges_tsen(:,:,k,ntsig)),minval(ges_tsen(:,:,k,ntsig))
! enddo
! do k=1,nsig
! write(*,*) 'prsl=',k,maxval(ges_prsl(:,:,k,ntsig)),minval(ges_prsl(:,:,k,ntsig))
! enddo
! write(*,*) 'sfct=',k,maxval(sfct(:,:,ntsig)),minval(sfct(:,:,ntsig))
! write(*,*) 'isli=',mype,maxval(isli2(:,:)),minval(isli2(:,:))
! write(*,*) 'soilT=',k,maxval(soil_temp(:,:,ntsig)),minval(soil_temp(:,:,ntsig))
! write(*,*) rlats
! write(*,*) rlons
! endif
!
!
l_use_hydroretrieval_all=.false.
krad_bot=7.0_r_single
r_radius=metar_impact_radius
r_radius_lowCloud=metar_impact_radius_lowCloud
opt_hydrometeor_retri=3 ! 1=Kessler 2=Lin 3=Thompson
opt_cloudtemperature=3 ! 3=latent heat, 4,5,6 = adiabat profile
opt_cloudwaterice_retri=1 ! 1 = RUC layer saturation and autoconvert
! 2 = convective
imerge_nesdis_nasalarc=2 ! =1 merge NASA LaRC with NESDIS
! =2 use NASA LaRC only
! =3 No Satellite cloud top used
! = other, use NESDIS only
!
! initialize the observation flag
!
istat_surface=0
istat_nesdis=0
istat_radar=0
istat_lightning=0
istat_nasalarc=0
!mhu need debug more: call load_gsdpbl_hgt(mype)
!
! check consistency of the options
!
! Now either stratiform or cumulus cloud is considered in the cloud
! water calculation. This leads to a limitation for the temperature
! adjustment when stratiform cloud is chosen because adiabat profile
! scheme based on the convection. This limitation may change when
! stratiform and cumulus cloud are both considered at the same time in the future.
if(opt_cloudwaterice_retri == 1 .and. opt_cloudtemperature >= 4) then
write(6,*) 'gsdcloudanalysis: ',&
'inconsistent option for opt_cloudwaterice_retri and opt_cloudtemperature'
write(6,*) 'gsdcloudanalysis: ',&
'opt_cloudtemperature must be set to 3 when opt_cloudwaterice_retri =1'
call stop2(113)
endif
!
!----------------------------------------------
! 2. read observations
!----------------------------------------------
!
! 1.1 allocate observation fields
!
allocate(ref_mos_3d(lon2,lat2,nsig))
allocate(ref_mos_3d_tten(lon2,lat2,nsig))
ref_mos_3d=miss_obs_real
ref_mos_3d_tten=miss_obs_real
allocate(lightning(lon2,lat2))
lightning=-9999.0_r_kind
allocate(sat_ctp(lon2,lat2))
allocate(sat_tem(lon2,lat2))
allocate(w_frac(lon2,lat2))
allocate(nlev_cld(lon2,lat2))
sat_ctp=miss_obs_real
sat_tem=miss_obs_real
w_frac=miss_obs_real
nlev_cld=miss_obs_int
allocate(osfc_station_map(lon2,lat2))
osfc_station_map=miss_obs_int
!
! 1.2 start to read observations
!
nmsclvl_radar = -999
lunin=55
open(lunin,file=obs_setup,form='unformatted')
rewind lunin
numsao=0
do is=1,ndat
if(dtype(is) /= ' ' .and. nsat1(is) > 0)then
!
! 1.2.2 read in surface observations
!
if( dtype(is) == 'mta_cld' ) then
numsao=nsat1(is)
allocate(oi(numsao))
allocate(oj(numsao))
allocate(ocld(nvarcld_p,numsao))
allocate(owx(numsao))
allocate(oelvtn(numsao))
allocate(odist(numsao))
allocate(cstation(numsao))
allocate(oistation(numsao))
allocate(ojstation(numsao))
allocate(wimaxstation(numsao))
call read_Surface(mype,lunin,regional_time,istart(mype+1),jstart(mype+1),lon2,lat2, &
numsao,nvarcld_p,oi,oj,ocld,owx,oelvtn,odist,cstation,oistation,ojstation)
if(mype == 0) write(6,*) 'gsdcloudanalysis: ', &
'Surface cloud observations are read in successfully'
istat_surface=1
!
! 1.2.4 read in NESDIS cloud products
!
elseif( dtype(is) == 'gos_ctp' ) then
call read_NESDIS(mype,lunin,nsat1(is),regional_time,istart(mype+1), &
jstart(mype+1),lon2,lat2,sat_ctp,sat_tem,w_frac,nesdis_npts_rad,ioption)
if(mype == 0) write(6,*) 'gsdcloudanalysis: ', &
'NESDIS cloud products are read in successfully'
istat_nesdis = 1
!
! 1.2.6 read in reflectivity mosaic
!
elseif( dtype(is) == 'rad_ref' ) then
allocate( ref_mosaic31(lon2,lat2,31) )
ref_mosaic31=-99999.0_r_kind
call read_radar_ref(mype,lunin,regional_time,istart(mype+1),jstart(mype+1), &
lon2,lat2,nmsclvl_radar,nsat1(is),ref_mosaic31)
if(mype == 0) write(6,*) 'gsdcloudanalysis: ', &
' radar reflectivity is read in successfully'
istat_radar=1
!
! 1.2.8 read in lightning
!
elseif( dtype(is)=='lghtn' ) then
call read_Lightning2cld(mype,lunin,regional_time,istart(mype+1),jstart(mype+1), &
lon2,lat2,nsat1(is),lightning)
if(mype == 0) write(6,*) 'gsdcloudanalysis: Lightning is read in successfully'
istat_lightning = 1
!
! 1.2.9 read in NASA LaRC cloud products
!
elseif( dtype(is) =='larcglb' ) then
allocate(nasalarc_cld(lon2,lat2,5))
nasalarc_cld=miss_obs_real
call read_map_nasalarc(mype,lunin,nsat1(is),regional_time,istart(mype+1), &
jstart(mype+1),lon2,lat2,nasalarc_cld,ioption)
if(mype == 0) write(6,*) 'gsdcloudanalysis:', &
'NASA LaRC cloud products are read in successfully', &
mype
istat_nasalarc = 1
!
! 1.2.12 all other observations
!
else
read(lunin) obstype,isis,nreal,nchanl
read(lunin)
endif ! dtype
endif
enddo ! is
close(lunin)
!
! 1.4 if there are NASA LaRC cloud products, use them to replace NESDIS ones.
! So we use NASA LaRC data in the same way as NESDIS ones
!
if(imerge_nesdis_nasalarc == 1 ) then
if(istat_nasalarc == 1 ) then
do j=2,lat2-1
do i=2,lon2-1
if(sat_ctp(i,j) < -99990.0_r_kind) then ! missing value is -999999.0
sat_ctp(i,j) = nasalarc_cld(i,j,1)
sat_tem(i,j) = nasalarc_cld(i,j,2)
w_frac(i,j) = nasalarc_cld(i,j,3)
nlev_cld(i,j)= int(nasalarc_cld(i,j,5))
istat_nesdis =istat_nasalarc
endif
enddo
enddo
endif
elseif ( imerge_nesdis_nasalarc == 2) then
if(istat_nasalarc == 1 ) then
sat_ctp(:,:) = nasalarc_cld(:,:,1)
sat_tem(:,:) = nasalarc_cld(:,:,2)
w_frac(:,:) = nasalarc_cld(:,:,3)
nlev_cld(:,:)= int(nasalarc_cld(:,:,5))
istat_nesdis =istat_nasalarc
endif
elseif ( imerge_nesdis_nasalarc == 3) then
istat_nesdis = 0
endif
!
!
! 1.6 collect the cloud information from whole domain
! and assign the background cloud to each METAR obs
!
allocate(sumqci(lon2,lat2,nsig))
if(i_gsdcldanal_type==30) then
do k=1,nsig
do j=1,lat2
do i=1,lon2
sumqci(i,j,k)= ges_ql(j,i,k)
enddo
enddo
enddo
else
do k=1,nsig
do j=1,lat2
do i=1,lon2
sumqci(i,j,k)= ges_ql(j,i,k) + ges_qi(j,i,k)
enddo
enddo
enddo
endif
allocate(watericemax(lon2,lat2))
allocate(kwatericemax(lon2,lat2))
watericemax=0._r_kind
wimaxstation=0.0_r_single
kwatericemax=-1
do j=1,lat2
do i=1,lon2
do k = 1,nsig
watericemax(i,j) = max(watericemax(i,j),sumqci(i,j,k))
end do
do k=1,nsig
if (sumqci(i,j,k) > cloud_def_p .and. kwatericemax(i,j) == -1) then
kwatericemax(i,j) = k
end if
end do
enddo
enddo
!
im=nlon
jm=nlat
allocate(all_loc(lat2,lon2))
allocate(strp(lat1*lon1))
allocate(tempa(itotsub))
allocate(temp1(im,jm))
do j=1,lat2
do i=1,lon2
all_loc(j,i) = watericemax(i,j)
enddo
enddo
call strip(all_loc,strp)
call mpi_allgatherv(strp,ijn(mype+1),mpi_real4, &
tempa,ijn,displs_g,mpi_real4,mpi_comm_world,ierror)
ierror=0
if(ierror /= 0 ) then
write(*,*) 'MPI error: cloud analysis=',mype
endif
temp1=0.0_r_single
call unfill_mass_grid2t(tempa,im,jm,temp1)
if(istat_surface==1) then
do ista=1,numsao
iob = min(max(int(oistation(ista)+0.5_r_kind),1),im)
job = min(max(int(ojstation(ista)+0.5_r_kind),1),jm)
wimaxstation(ista)=temp1(iob,job)
if(wimaxstation(ista) > 0._r_single) then
i=int(oi(ista))
j=int(oj(ista))
endif
enddo
endif
deallocate(all_loc,strp,tempa,temp1)
! make a surface station map in grid coordinate
if(istat_surface==1) then
do ista=1,numsao
iob = int(oistation(ista)-jstart(mype+1)+2)
job = int(ojstation(ista)-istart(mype+1)+2)
if(iob >=1 .and. iob<=lon2-1 .and. job >=1 .and. job<=lat2-1) then
osfc_station_map(iob,job)=1
osfc_station_map(iob+1,job)=1
osfc_station_map(iob,job+1)=1
osfc_station_map(iob+1,job+1)=1
endif
enddo
endif
!
! 1.8 check if data available: if no data in this subdomain, return.
!
!mhu if( (istat_radar + istat_surface + istat_nesdis + istat_lightning ) == 0 ) then
if( (istat_nesdis + istat_surface) == 0 ) then
!mhu ges_cwmr(:,:,nsig)=1003.0
write(6,*) ' No cloud observations available, return', mype
deallocate(ref_mos_3d,ref_mos_3d_tten,lightning,sat_ctp,sat_tem,w_frac,nlev_cld)
return
endif
!
!----------------------------------------------
! 2. allocated background arrays and read background
! further observation data process before cloud analysis
!----------------------------------------------
!
! 2.2 allocate background and analysis fields
!
allocate(t_bk(lon2,lat2,nsig))
allocate(h_bk(lon2,lat2,nsig))
allocate(p_bk(lon2,lat2,nsig))
allocate(ps_bk(lon2,lat2))
allocate(zh(lon2,lat2))
allocate(q_bk(lon2,lat2,nsig))
allocate(xlon(lon2,lat2))
allocate(xlat(lon2,lat2))
allocate(xland(lon2,lat2))
allocate(soiltbk(lon2,lat2))
! allocate(z_lcl(lon2,lat2))
! allocate(pblh(lon2,lat2))
t_bk=miss_obs_real
h_bk=miss_obs_real
p_bk=miss_obs_real
ps_bk=miss_obs_real
zh=miss_obs_real
q_bk=miss_obs_real
xlon=miss_obs_real
xlat=miss_obs_real
xland=miss_obs_real
soiltbk=miss_obs_real
! pblh=miss_obs_real
allocate(cldwater_3d(lon2,lat2,nsig))
allocate(cldice_3d(lon2,lat2,nsig))
allocate(rain_3d(lon2,lat2,nsig))
allocate(nrain_3d(lon2,lat2,nsig))
allocate(snow_3d(lon2,lat2,nsig))
allocate(graupel_3d(lon2,lat2,nsig))
allocate(cldtmp_3d(lon2,lat2,nsig))
cldwater_3d=miss_obs_real
cldice_3d=miss_obs_real
rain_3d=miss_obs_real
nrain_3d=miss_obs_real
snow_3d=miss_obs_real
graupel_3d=miss_obs_real
cldtmp_3d=miss_obs_real
allocate(rain_1d_save(nsig))
allocate(nrain_1d_save(nsig))
allocate(snow_1d_save(nsig))
rain_1d_save=miss_obs_real
nrain_1d_save=miss_obs_real
snow_1d_save=miss_obs_real
!
! 2.4 read in background fields
!
if(i_gsdcldanal_type==30) then
!
do j=1,lat2
do i=1,lon2
ii = i + jstart(mype+1) - 2 ! covert it to the global grid
jj = j + istart(mype+1) - 2 ! covert it to the global grid
ii=max(1,min(ii,nlon))
jj=max(1,min(jj,nlat))
xlon(i,j) = rlons(ii)*rad2deg ! longitude back to degree
xlat(i,j) = rlats(jj)*rad2deg ! latitude back to degree
enddo
enddo
do j=1,lat2
do i=1,lon2
zh(i,j) =ges_z(j,i) ! terrain in meter
ps_bk(i,j) =ges_ps(j,i)*10.0_r_single ! surace pressure in mb
xland(i,j) =isli2(j,i) ! 0=water, 1=land, 2=ice
soiltbk(i,j)=soil_temp(j,i,ntsig) ! soil temperature
enddo
enddo
do k=1,nsig
do j=1,lat2
do i=1,lon2
p_bk(i,j,k)=ges_prsl(j,i,k,ntsig)*10.0_r_single
enddo
enddo
enddo
else
do j=1,lat2
do i=1,lon2
zh(i,j) =ges_z(j,i) ! terrain in meter
ps_bk(i,j) =ges_ps(j,i)*10.0_r_single ! surace pressure in mb
xland(i,j) =isli_cld(j,i,itsfc) ! 0=water, 1=land, 2=ice
soiltbk(i,j)=soil_temp_cld(j,i,itsfc) ! soil temperature
xlon(i,j) =ges_xlon(j,i,itsfc)*rad2deg ! longitude back to degree
xlat(i,j) =ges_xlat(j,i,itsfc)*rad2deg ! latitude back to degree
enddo
enddo
endif
do k=1,nsig
do j=1,lat2
do i=1,lon2
q_bk(i,j,k)=ges_q(j,i,k) ! specific humidity
qmixr = q_bk(i,j,k)/(one - q_bk(i,j,k)) ! covert from specific
! humidity to mixing ratio
t_bk(i,j,k)=ges_tsen(j,i,k,ntsig)
enddo
enddo
enddo
if(i_gsdcldanal_type==30) then
call BackgroundCldgfs(mype,lon2,lat2,nsig,t_bk,p_bk,ps_bk,q_bk,h_bk)
else
call BackgroundCld(mype,lon2,lat2,nsig,t_bk,p_bk,ps_bk,q_bk,h_bk, &
zh,pt_ll,eta1_ll,aeta1_ll,regional,wrf_mass_regional)
endif
!if(mype==25) then
! write(*,*) lat2,lon2,nsig
! write(*,*) 'ps=',maxval(ps_bk),minval(ps_bk)
! write(*,*) 'zh=',maxval(zh),minval(zh)
! write(*,*) 'xland=',maxval(xland),minval(xland)
! write(*,*) 'soiltbk=',maxval(soiltbk),minval(soiltbk)
! write(*,*) 'xlon=',maxval(xlon),minval(xlon)
! write(*,*) 'xlat=',maxval(xlat),minval(xlat)
! do k=1,nsig
! write(*,*) 't=',k,maxval(t_bk(:,:,k)),minval(t_bk(:,:,k))
! enddo
! do k=1,nsig
! write(*,*) 'q=',k,maxval(q_bk(:,:,k)),minval(q_bk(:,:,k))
! enddo
! do k=1,nsig
! write(*,*) 'cw=',k,maxval(ges_cwmr(:,:,k)),minval(ges_cwmr(:,:,k))
! enddo
! do k=1,nsig
! write(*,*) 'p_bk=',k,maxval(p_bk(:,:,k)),minval(p_bk(:,:,k))
! enddo
! do k=1,nsig
! write(*,*) 'h_bk=',k,maxval(h_bk(:,:,k)),minval(h_bk(:,:,k))
! enddo
!endif
!
! 2.5 calculate PBL height
!
!mhu call calc_pbl_height(mype,lat2,lon2,nsig,q_bk,t_bk,p_bk,pblh)
!
!
!----------------------------------------------
! 3. Calculate 3-d cloud cover obs-information field (cld_cover_3d),
! cloud type, precipitation type
!----------------------------------------------
!
allocate(cld_cover_3d(lon2,lat2,nsig))
allocate(cld_type_3d(lon2,lat2,nsig))
allocate(wthr_type_2d(lon2,lat2))
allocate(pcp_type_3d(lon2,lat2,nsig))
allocate(cloudlayers_i(lon2,lat2,21))
cld_cover_3d=miss_obs_real
cld_type_3d =miss_obs_int
wthr_type_2d=miss_obs_int
pcp_type_3d =miss_obs_int
allocate(vis2qc(lon2,lat2))
vis2qc=miss_obs_real
!
!
if(istat_surface == 1) then
call cloudCover_surface(mype,lat2,lon2,nsig,r_radius,thunderRadius, &
t_bk,p_bk,q_bk,h_bk,zh, &
numsao,nvarcld_p,numsao,oi,oj,ocld,owx,oelvtn,odist, &
cld_cover_3d,cld_type_3d,wthr_type_2d,pcp_type_3d, &
wimaxstation, kwatericemax,vis2qc)
if(mype == 0) write(6,*) 'gsdcloudanalysis:', &
'success in cloud cover analysis using surface data'
endif
if(istat_nesdis == 1 ) then
call cloudCover_NESDIS(mype,regional_time,lat2,lon2,nsig, &
xlon,xlat,t_bk,p_bk,h_bk,xland, &
soiltbk,sat_ctp,sat_tem,w_frac, &
l_cld_bld,cld_bld_hgt, &
build_cloud_frac_p,clear_cloud_frac_p,nlev_cld, &
cld_cover_3d,cld_type_3d,wthr_type_2d,osfc_station_map)
if(mype == 0) write(6,*) 'gsdcloudanalysis:', &
' success in cloud cover analysis using NESDIS data'
endif
!----------------------------------------------
! 4. Calculate 3-d cloud water and ice
! Calculate 3-d hydrometeors
! Calculate radar temperature tendency
! Calculate in cloud temperature
! moisture adjustment (to do)
!----------------------------------------------
!
! 4.2 find the cloud layers
!
if(istat_nesdis == 1 ) then
call cloudLayers(lat2,lon2,nsig,h_bk,zh,cld_cover_3d, &
cld_type_3d,cloudlayers_i)
if(mype==0) write(6,*) 'gsdcloudanalysis: success in finding cloud layers'
endif
!
! 4.4 decide the cloud type
!
!
! 4.6 calculate liquid water content
!
if(istat_nesdis == 1 ) then
call cloudLWC_stratiform(mype,lat2,lon2,nsig,q_bk,t_bk,p_bk, &
cld_cover_3d,cld_type_3d,wthr_type_2d,cloudlayers_i, &
cldwater_3d,cldice_3d)
if(mype==0) write(6,*) 'gsdcloudanalysis: ', &
'success in modifying hydrometeors for stratiform clouds '
endif
!
!----------------------------------------------
! 5. the final analysis or update background
!----------------------------------------------
!
! the final analysis of cloud
!
if(i_gsdcldanal_type==30 ) then
!if(istat_nasalarc == 1) then
if((istat_nasalarc + istat_nesdis + istat_surface) > 1) then
do k=1,nsig
do j=2,lat2-1
do i=2,lon2-1
cldwater_3d(i,j,k)=cldwater_3d(i,j,k)+cldice_3d(i,j,k)
if(cldwater_3d(i,j,k) < 0.0_r_kind) cldwater_3d(i,j,k)=zero
if( cld_cover_3d(i,j,k) > -0.001_r_kind ) then
if( cld_cover_3d(i,j,k) > 0.6_r_kind ) then
cldwater_3d(i,j,k) = max(0.001_r_kind*cldwater_3d(i,j,k), &
ges_cwmr(j,i,k))
else ! clean cloud
cldwater_3d(i,j,k) = zero
endif
else ! unknown, using background values
cldwater_3d(i,j,k) = ges_cwmr(j,i,k)
endif
cldwater_3d(i,j,k)= max(0.0_r_single,cldwater_3d(i,j,k))
end do
end do
end do
! call cloud_saturation(mype,l_conserve_thetaV,i_conserve_thetaV_iternum, &
! lat2,lon2,nsig,q_bk,t_bk,p_bk, &
! cld_cover_3d,wthr_type_2d,cldwater_3d,cldice_3d,sumqci)
!----------------------------------------------
! 6. save the analysis results
!----------------------------------------------
!
!
! do j=2,lat2-1
! do i=2,lon2-1
! write(*,*) mype,i,j,sat_ctp(i,j)
! enddo
! enddo
!do k=1,nsig
! write(*,'(a,2I5,2e15.5)') 'update the ges_cwmr for',mype,k, &
! maxval(ges_cwmr(:,:,k)),minval(ges_cwmr(:,:,k))
!enddo
do k=1,nsig-1
do j=2,lat2-1
do i=2,lon2-1
ges_cwmr(j,i,k)=cldwater_3d(i,j,k)
if(ges_cwmr(j,i,k)<0.0_r_kind) ges_cwmr(j,i,k)=0.0_r_kind
enddo
enddo
enddo
!
! do j=1,lat2
! write(400+mype,'(2I10,40f8.1)') mype,j,(sat_ctp(i,j),i=1,lon2)
! enddo
!!
! ges_cwmr(:,:,nsig)=1003.0
! do j=2,lat2-1
!! do i=2,lon2-1
! if(sat_ctp(i,j) < 0.0 .or. sat_ctp(i,j) > 1015.0 ) sat_ctp(i,j)=1015.0
! ges_cwmr(j,i,k)=sat_ctp(i,j)
! enddo
! enddo
!
!do k=1,nsig
! write(*,'(a,2I5,2e15.5)') 'after update the ges_cwmr for',mype,k,&
! maxval(ges_cwmr(:,:,k)),minval(ges_cwmr(:,:,k))
!enddo
endif
else !regional case
endif ! regional or global
!
!----------------------------------------------
! 7. release space
!----------------------------------------------
!
deallocate(cld_cover_3d,cld_type_3d,wthr_type_2d, &
pcp_type_3d,cloudlayers_i)
deallocate(t_bk,h_bk,p_bk,ps_bk,zh,q_bk,sumqci)
deallocate(xlon,xlat,xland,soiltbk)
deallocate(cldwater_3d,cldice_3d,rain_3d,nrain_3d,snow_3d,graupel_3d,cldtmp_3d)
deallocate(vis2qc)
if(istat_surface == 1 ) then
deallocate(oi,oj,ocld,owx,oelvtn,odist,cstation,oistation,ojstation,wimaxstation)
deallocate(watericemax,kwatericemax)
endif
if(istat_nasalarc == 1 ) then
deallocate(nasalarc_cld)
endif
deallocate(sat_ctp,sat_tem,w_frac,nlev_cld)
deallocate(ref_mos_3d,ref_mos_3d_tten,lightning)
if(mype==3) then
write(6,*) '========================================'
write(6,*) 'gsdcloudanalysis: generalized cloud analysis finished:',mype
write(6,*) '========================================'
endif
#else /* Start no RR cloud analysis library block */
implicit none
! Declare passed variables
integer(i_kind),intent(in):: mype
!
if( mype == 0) write(6,*)'gsdcloudanalysis: dummy routine, does nothing!'
#endif /* End no RR cloud analysis library block */
end subroutine gsdcloudanalysis4gfs