#ifdef NMMB_CLOUDANALYSIS
SUBROUTINE gsdcloudanalysis4NMMB(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-20 s.liu - use cloud analysis for NMMB
! 2014-11-24 s.liu - use radar-derived qr and qs in cloud analysis
! 2015-01-24 s.liu - move cloud analysis related parameters used for NMMB to this subroutine
! 2015-03-20 s.liu - tune the method of using radar-derived qr and qs
! 2016-02-26 s.liu - use namelist to control cloud analysis parameters
! 2016-02-29 s.liu - define l_use_hydroretrieval_all as public variable
! 2016-05-05 s.liu - consider possible bias of derived reflecitvity from
! lightnint density based on region mask
!
!
! 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,rd_over_cp,zero_single
use constants, only: rd_over_cp, h1000
use kinds, only: r_single,i_kind, r_kind
use mpimod, only: ierror,mpi_real4,mpi_real8,mpi_sum,mpi_comm_world
use gridmod, only: pt_ll,eta1_ll,aeta1_ll
use gridmod, only: regional,wrf_mass_regional,regional_time,region_lon,region_lat
use gridmod, only: nsig,lat2,lon2,istart,jstart,nlon,nlat
use gridmod, only: itotsub,lon1,lat1,nlon_regional,nlat_regional,ijn,displs_g
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
use guess_grids, only: soil_temp,isli2
use guess_grids, only: ges_tsen,geop_hgtl
use guess_grids, only: ges_prsl
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, &
l_use_hydroretrieval_all, &
iclean_hydro_withRef, &
iclean_hydro_withRef_allcol, &
ioption
use gsi_metguess_mod, only: GSI_MetGuess_Bundle
use gsi_metguess_mod, only: GSI_MetGuess_get
use gsi_bundlemod, only: gsi_bundlegetpointer
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)
real(r_kind),pointer,dimension(:,:) :: ges_z =>NULL()
real(r_kind),pointer,dimension(:,:) :: ges_ps=>NULL()
real(r_kind),pointer,dimension(:,:,:) :: ges_q =>NULL()
! real(r_kind),pointer,dimension(:,:,:) :: ges_tv=>NULL()
!
! 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(:)
!
! lightning observation: 2D field in RR grid
!
real(r_kind),allocatable :: lightning(:,:),lghtn_region_mask(:,:)
real(r_kind),allocatable :: lghtn_ref_bias(:,:)
!
! 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(:,:)
real(r_single), allocatable :: sat_ctp_nesdis(:,:)
real(r_single), allocatable :: sat_tem_nesdis(:,:)
real(r_single), allocatable :: w_frac_nesdis(:,:)
!
! 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 :: max_obs_qrqs_save(:,:)
real(r_single),allocatable :: max_bk_qrqs_save(:,:)
real(r_single),allocatable :: ratio_qrqs_save(:,:)
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
REAL(r_single) :: miss_obs_single
PARAMETER ( miss_obs_int = -99999 )
PARAMETER ( miss_obs_real = -99999.0_r_kind )
PARAMETER ( miss_obs_single = -99999.0_r_single )
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
REAL(r_single),allocatable :: sumq(:,:,:) ! total liquid water
!
! 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_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_ref(:,:,:)=>NULL() ! ref
REAL(r_kind), pointer :: ges_tten(:,:,:)=>NULL() ! t-sensible
REAL(r_kind), pointer :: dfi_tten(:,:,:)=>NULL() ! tten
!
! 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,imaxlvl_modref,i1,j1,i2,j2
real(r_kind) :: max_retrieved_qrqs,max_bk_qrqs,ratio_hyd_bk2obs
REAL(r_kind) :: qrlimit
character(10) :: obstype
integer(i_kind) :: lunin, is, ier, istatus, ivar
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
real(r_kind) :: t, qfi, qfs, qfr, qfw, t_inc
real(r_kind) :: iss, jss
!* output grads
real(r_single)::outputgrads(lat2,lon2,nsig)
real(r_single)::cref(lon2,lat2)
real(r_single),parameter:: ze_qr_const=3.6308*1.0e9_r_single
real(r_single),parameter:: ze_qs_const=3.268*1.0e9_r_single
!
!
itsig=1 ! _RT shouldn't this be ntguessig?
itsfc=1 ! _RT shouldn't this be ntguessig?
!
call gsi_metguess_get('var::qi',ivar,ier)
ier=0
! 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),'z' ,ges_z, istatus);ier=ier+istatus
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(itsig),'ps',ges_ps,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),'ql',ges_ql,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),'ref',ges_ref,istatus);ier=ier+istatus
call gsi_bundlegetpointer (GSI_MetGuess_Bundle(itsig),'tten',dfi_tten,istatus);ier=ier+istatus
allocate(ges_tten(lat2,lon2,nsig))
if(mype==0)then
write(6,*)'gsdcloudanalysis4NMMB start ::'
! write(6,*)'sliu output bundle gsd::',GSI_MetGuess_Bundle(itsig)%r3(8)%shortname
end if
if(ier/=0) write(6,*)'gsdcloudanalysis :: no cloud guess, nothing to do RETURN'
if(ier/=0) return ! no guess, nothing to do
if(mype==0) then
write(6,*) '========================================'
write(6,*) 'gsdcloudanalysis: Start generalized cloud analysis '
write(6,*) '========================================'
endif
!
!
!
krad_bot=7.0_r_single
r_radius=metar_impact_radius
r_radius_lowCloud=metar_impact_radius_lowCloud
dfi_lhtp=dfi_radar_latent_heat_time_period
if(mype==0) then
write(6,*)'metar_impact_radius::',metar_impact_radius
write(6,*)'metar_impact_radius_lowCloud::',metar_impact_radius_lowCloud
write(6,*)'l_cld_bld::,',l_cld_bld
write(6,*)'cld_bld_hgt::,',cld_bld_hgt
write(6,*)'build_cloud_frac_p::,',build_cloud_frac_p
write(6,*)'clear_cloud_frac_p::,',clear_cloud_frac_p
write(6,*)'l_use_hydroretrieval_all::',l_use_hydroretrieval_all
write(6,*)'iclean_hydro_withRef::',iclean_hydro_withRef
write(6,*)'iclean_hydro_withRef_allcol::',iclean_hydro_withRef_allcol
write(6,*)'dfi_radar_latent_heat_time_period::',dfi_lhtp
write(6,*)'l_conserve_thetaV::',l_conserve_thetaV
write(6,*)'i_conserve_thetaV_iternum::',i_conserve_thetaV_iternum
end if
!* change option to Ferrier
opt_hydrometeor_retri=2 ! 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
!
! initialize the observation flag
!
istat_Surface=0
istat_NESDIS=0
istat_radar=1
istat_lightning=0
istat_NASALaRC=0
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))
allocate(lghtn_region_mask(lon2,lat2))
allocate(lghtn_ref_bias(lon2,lat2))
lightning=-9999.0_r_kind
lghtn_region_mask=0.0_r_kind
lghtn_ref_bias=0.0_r_kind
allocate(sat_ctp_nesdis(lon2,lat2))
allocate(sat_tem_nesdis(lon2,lat2))
allocate(w_frac_nesdis(lon2,lat2))
sat_ctp_nesdis=miss_obs_real
sat_tem_nesdis=miss_obs_real
w_frac_nesdis=miss_obs_real
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
!
! 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))
allocate(kwatericemax(lon2,lat2))
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
elseif( dtype(is) == 'gos_ctp' ) then
!
! 1.2.4 read in NESDIS cloud products
!
write(6,*) 'gsdcloudanalysis: ', &
'start to read NESDIS cloud products'
call read_NESDIS(mype,lunin,nsat1(is),regional_time,istart(mype+1), &
jstart(mype+1),lon2,lat2,sat_ctp,sat_tem,w_frac,ioption)
sat_ctp_nesdis=sat_ctp
sat_tem_nesdis=sat_tem
w_frac_nesdis=w_frac
istat_NESDIS = 1
elseif( dtype(is) == 'rad_ref' ) then
!
! 1.2.6 read in reflectivity mosaic
!
! allocate( ref_mosaic31(lon2,lat2,31) )
! ref_mosaic31=-9999.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
elseif( dtype(is)=='lghtn' ) then
!
! 1.2.8 read in lightning
!
call read_Lightningbufr2cld(mype,lunin,regional_time,istart(mype+1),jstart(mype+1), &
lon2,lat2,nsat1(is),lightning)
istat_lightning = 1
elseif( dtype(is) =='larccld' ) then
!
! 1.2.9 read in NASA LaRC cloud products
!
allocate(nasalarc_cld(lon2,lat2,5))
nasalarc_cld=miss_obs_real
call read_NASALaRC(mype,lunin,nsat1(is),regional_time,istart(mype+1), &
jstart(mype+1),lon2,lat2,nasalarc_cld)
if(mype == 0) write(6,*) 'gsdcloudanalysis:', &
'NASA LaRC cloud products are read in successfully'
istat_NASALaRC = 1
else
!
! 1.2.12 all other observations
!
read(lunin) obstype,isis,nreal,nchanl
read(lunin)
endif ! dtype
endif
enddo ! is
close(lunin)
88 format('gsdcloud Lightning::',7i6)
!
! 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
!
! for NMMB, just merge NESDIS and NASA cloud product
if(istat_NASALaRC == 1 ) then
DO j=2,lat2-1
DO i=2,lon2-1
if(sat_ctp(i,j) < -99990.0_r_single) 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
! else
! write(6,*)'merge sat data::',i,j,sat_tem(i,j)
endif
ENDDO
ENDDO
endif
do k=1,nsig
do j=1,lat2
do i=1,lon2
ref_mos_3d(i,j,k)=ges_ref(j,i,k) !
if(abs(ref_mos_3d(i,j,k))<100.0_r_single) istat_radar=1
end do
end do
end do
!
! 1.6 collect the cloud information from whole domain
! and assign the background cloud to each METAR obs
!
allocate(sumqci(lon2,lat2,nsig))
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
! 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_regional
! jm=nlat_regional
! 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_single(all_loc(1,1),strp,1)
! 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),1),im)
! job = min(max(int(ojstation(ista)+0.5),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)
!
! 1.6 check if data available: if no data in this subdomain, return.
!
if( (istat_radar + istat_Surface + istat_NESDIS + istat_lightning ) == 0 ) then
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(sumq(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))
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
!
! 2.4 read in background fields
!
do j=1,lat2
do i=1,lon2
iss=jstart(mype+1)+i-1
jss=istart(mype+1)+j-1
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)
soiltbk(i,j)=soil_temp(j,i,itsfc) ! soil temperature
xlon(i,j) =region_lon(jss,iss)*rad2deg ! longitude back to degree
xlat(i,j) =region_lat(jss,iss)*rad2deg ! latitude back to degree
! xland(i,j) =isli_cld(j,i,itsfc) ! 0=water, 1=land, 2=ice
! 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
! soiltbk(i,j)=0.0
! xlon(i,j) =0.0
! xlat(i,j) =0.0
! write(6,*)'xlon, xlat::',iss,jss,region_lon(jss,iss)*rad2deg,region_lat(jss,iss)*rad2deg,xland(i,j)
ENDDO
ENDDO
do k=1,nsig
do j=1,lat2
do i=1,lon2
q_bk(i,j,k)=ges_q(j,i,k)/(1.0_r_kind-ges_q(j,i,k)) ! specific humidity to mixing ratio, ges_q is spfh
h_bk(i,j,k)=geop_hgtl(j,i,k,1) !
! t_bk(i,j,k)=ges_tv(j,i,k,itsig)/ &
! (1.0_r_single+0.61_r_single*q_bk(i,j,k)) ! virtual temp to temp
t_bk(i,j,k)=ges_tsen(j,i,k,1)*(100.0_r_kind/ges_prsl(j,i,k,1))**rd_over_cp !to potential temperature
p_bk(i,j,k)=ges_prsl(j,i,k,1)*10.0_r_kind
ENDDO
ENDDO
ENDDO
! 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)
!
! 2.5 calculate PBL height
!
call calc_pbl_height(mype,lat2,lon2,nsig,q_bk,t_bk,p_bk,pblh)
!
! 2.6 vertical interpolation of radar reflectivity
!
if(istat_radar == 1 ) then
! call vinterp_radar_ref(mype,lon2,lat2,nsig,Nmsclvl_radar, &
! ref_mos_3d,ref_mosaic31,h_bk,zh)
! deallocate( ref_mosaic31 )
ref_mos_3d_tten=ref_mos_3d
call build_missing_REFcone(mype,lon2,lat2,nsig,krad_bot,ref_mos_3d_tten,h_bk,pblh)
end if
!
! 2.8 convert lightning to reflectivity
!
do j=1,lat2
do i=1,lon2
cref(i,j)=-999.0_r_single
do k=1,nsig
if(abs(ref_mos_3d(i,j,k))<100.0_r_single)then
cref(i,j) = max(cref(i,j),ref_mos_3d(i,j,k))
end if
end do
end do
end do
do j=1,lat2
do i=1,lon2
do i1=-1,1
do j1=-1,1
i2=i+i1; j2=j+j1
if(i2>0.and.j2>0) then
if(abs(cref(i2,j2))<100.0_r_single) then
! lightning(i,j)=-999.0_r_kind
lghtn_region_mask(i,j)=1.0_r_kind !* with radar coverage
end if
end if
end do
end do
end do
end do
if(istat_lightning == 1 ) then
call convert_lghtn2ref(mype,lon2,lat2,nsig,ref_mos_3d_tten, &
lightning,lghtn_region_mask,lghtn_ref_bias,h_bk)
endif
!
!
!----------------------------------------------
! 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
!
!
if(istat_Surface == 1) then
wimaxstation=0.0_r_kind; kwatericemax=0
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)
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,zh,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)
if(mype == 0) write(6,*) 'gsdcloudanalysis:', &
' success in cloud cover analysis using NESDIS data'
endif
! for Rapid Refresh application, turn off the radar reflectivity impact
! if(istat_radar == 1 .or. istat_lightning == 1 ) then
! write(6,*)'gsdcloudanalysis sliu max reflectivity::', maxval(ref_mos_3d)
! call cloudCover_radar(mype,lat2,lon2,nsig,h_bk,zh,ref_mos_3d, &
! cld_cover_3d,cld_type_3d,wthr_type_2d)
! if(mype == 0) write(6,*) 'gsdcloudanalysis: ', &
! ' success in cloud cover analysis using radar data'
! endif
! for NMMB, turn on the radar reflectivity impact
call cloudCover_radar(mype,lat2,lon2,nsig,h_bk,zh,ref_mos_3d_tten, &
cld_cover_3d,cld_type_3d,wthr_type_2d)
!
!----------------------------------------------
! 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
!
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'
!
! 4.4 decide the cloud type
!
call cloudType(lat2,lon2,nsig,h_bk,t_bk,p_bk,ref_mos_3d, &
cld_cover_3d,cld_type_3d,wthr_type_2d,cloudlayers_i)
if(mype==0) write(6,*) 'gsdcloudanalysis: success in deciding cloud types'
!
! 4.6 calculate liquid water content
!
if(opt_cloudwaterice_retri == 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 '
elseif (opt_cloudwaterice_retri == 2) then
call cloudLWC_Cumulus(lat2,lon2,nsig,h_bk,t_bk,p_bk, &
cld_cover_3d,cld_type_3d,wthr_type_2d,cloudlayers_i, &
cldwater_3d,cldice_3d,cldtmp_3d)
if(mype==0) write(6,*) 'gsdcloudanalysis: ', &
' success in modifying hydrometeors from radar reflectivity'
else
write(6,*)'gsdcloudanalysis: ', &
' Invalid cloud water calculation option, check opt_cloudwaterice_retri'
call stop2(113)
endif
!
! 4.8 calculate hydrometeors
!
istat_radar=1
where(rain_3d<0.0_r_single) rain_3d=0.0_r_single
where(snow_3d<0.0_r_single) snow_3d=0.0_r_single
if(istat_radar == 1 .or. istat_lightning == 1) then
call PrecipType(lat2,lon2,nsig,t_bk,p_bk,q_bk,ref_mos_3d, &
wthr_type_2d,pcp_type_3d)
if(mype==0) write(6,*) 'gsdcloudanalysis: ', &
' success in deciding precipitation type'
call PrecipMxR_radar(mype,lat2,lon2,nsig, &
t_bk,p_bk,ref_mos_3d, q_bk, &
pcp_type_3d,rain_3d,nrain_3d,snow_3d,graupel_3d,opt_hydrometeor_retri)
!* note: for NMMB get qr, qs in kg/kg
if(mype==0) write(6,*) 'gsdcloudanalysis: ', &
' success in determining hydrometeor types from radar refl'
endif
!
! 4.10 radar temperature tendency for DFI
!
! dfi_lhtp=20 ! 20 min forward intergration
if (istat_radar==1) then
if (istat_NESDIS==1) then
call radar_ref2tten(mype,istat_radar,istat_lightning,lon2,lat2,nsig,ref_mos_3d_tten, &
cld_cover_3d,p_bk,t_bk,ges_tten,dfi_lhtp,krad_bot,pblh,sat_ctp)
if(mype==0) write(6,*) 'gsdcloudanalysis: ', &
' success in getting dfi_tten from radar refl, istat_NESDIS==1'
else
call radar_ref2tten_nosat(mype,istat_radar,istat_lightning,lon2,lat2,nsig, &
ref_mos_3d_tten,cld_cover_3d,p_bk,t_bk,ges_tten,dfi_lhtp,krad_bot,pblh)
if(mype==0) write(6,*) 'gsdcloudanalysis: ', &
' success in getting dfi_tten from radar refl, istat_NESDIS==0'
endif
endif
!
! 4.12 temperature adjustment
!
! call TempAdjust(mype,lat2,lon2,nsig,opt_cloudtemperature, t_bk, p_bk, w_bk, q_bk, &
! cldwater_3d,cldice_3d,cldtmp_3d)
!
!----------------------------------------------
! 5. the final analysis or update background
!----------------------------------------------
!
! the final analysis of cloud
DO k=1,nsig
DO j=2,lat2-1
DO i=2,lon2-1
sumq(i,j,k)= ges_ql(j,i,k) + ges_qi(j,i,k)
if( cld_cover_3d(i,j,k) > -0.001_r_single ) then
if( cld_cover_3d(i,j,k) > 0.6_r_single ) then
cldwater_3d(i,j,k) = max(0.001_r_single*cldwater_3d(i,j,k),ges_ql(j,i,k))
cldice_3d(i,j,k) = max(0.001_r_single*cldice_3d(i,j,k),ges_qi(j,i,k))
else if(cld_cover_3d(i,j,k)<0.1_r_single) then ! clean cloud
cldwater_3d(i,j,k) = zero_single
cldice_3d(i,j,k) = zero_single
else
cldwater_3d(i,j,k) = ges_ql(j,i,k)
cldice_3d(i,j,k) = ges_qi(j,i,k)
endif
else ! unknown, using background values
cldwater_3d(i,j,k) = ges_ql(j,i,k)
cldice_3d(i,j,k) = ges_qi(j,i,k)
endif
END DO
END DO
END DO
!
! the final analysis of precipitation
!
allocate(max_bk_qrqs_save(lon2,lat2))
allocate(max_obs_qrqs_save(lon2,lat2))
allocate(ratio_qrqs_save(lon2,lat2))
max_bk_qrqs_save=-999.0_r_single
max_obs_qrqs_save=-999.0_r_single
ratio_qrqs_save=-999.0_r_single
ref_mos_3d=ref_mos_3d_tten
if(l_use_hydroretrieval_all) then !* use retr_qr, retr_qs block
qrlimit=15.0_r_kind*0.001_r_kind
DO k=1,nsig-1
DO j=2,lat2-1
DO i=2,lon2-1
snowtemp=snow_3d(i,j,k)
raintemp=rain_3d(i,j,k)
! nraintemp=nrain_3d(i,j,k)
rain_3d(i,j,k) = ges_qr(j,i,k)
! nrain_3d(i,j,k)= ges_qnr(j,i,k)
snow_3d(i,j,k) = ges_qs(j,i,k)
! graupel_3d(i,j,k) = ges_qg(j,i,k)
if(ref_mos_3d(i,j,k) > zero ) then
! snow_3d(i,j,k) = MIN(max(max(snowtemp,zero)*0.001_r_kind,ges_qs(j,i,k)),qrlimit)
snow_3d(i,j,k) = MIN(max(max(snowtemp,zero),ges_qs(j,i,k)),qrlimit)
! rain_3d(i,j,k) = MIN(max(max(raintemp,zero)*0.001_r_kind,ges_qr(j,i,k)),qrlimit)
! raintemp = max(raintemp,zero)*0.001_r_kind
raintemp = max(raintemp,zero)
if(raintemp > ges_qr(j,i,k) ) then
if(raintemp <= qrlimit) then
rain_3d(i,j,k) = raintemp
! nrain_3d(i,j,k)= nraintemp
else
rain_3d(i,j,k) = qrlimit
! nrain_3d(i,j,k)= nraintemp*(qrlimit/raintemp)
endif
else
rain_3d(i,j,k) = MIN(ges_qr(j,i,k),qrlimit)
endif
elseif( ref_mos_3d(i,j,k) <= zero .and. &
ref_mos_3d(i,j,k) > -100.0_r_kind ) then
rain_3d(i,j,k) = zero
! nrain_3d(i,j,k) = zero
snow_3d(i,j,k) = zero
! graupel_3d(i,j,k) = zero
else
rain_3d(i,j,k) = ges_qr(j,i,k)
! nrain_3d(i,j,k)= ges_qnr(j,i,k)
snow_3d(i,j,k) = ges_qs(j,i,k)
! graupel_3d(i,j,k) = ges_qg(j,i,k)
endif
END DO
END DO
END DO
else ! hydrometeor anlysis for NMM forecast
qrlimit=3.0_r_kind*0.001_r_kind
DO j=2,lat2-1
DO i=2,lon2-1
refmax=-999.0_r_kind
imaxlvl_ref=0
DO k=1,nsig
if(ref_mos_3d(i,j,k) > refmax) then
imaxlvl_ref=k
refmax=ref_mos_3d(i,j,k)
endif
!! rain_3d(i,j,k)=max(rain_3d(i,j,k)*0.001_r_kind,zero)
!! snow_3d(i,j,k)=max(snow_3d(i,j,k)*0.001_r_kind,zero)
rain_3d(i,j,k)=max(rain_3d(i,j,k),zero)
snow_3d(i,j,k)=max(snow_3d(i,j,k),zero)
!! ges_qnr(i,j,k)=max(ges_qnr(i,j,k),zero)
ENDDO
!* if refmax > 0 and refmax above the model bottom and below model top
if( refmax > 0 .and. (imaxlvl_ref > 0 .and. imaxlvl_ref < nsig ) ) then ! use retrieval hybrometeors
! tsfc=t_bk(i,j,1)*(p_bk(i,j,1)/h1000)**rd_over_cp - 273.15_r_kind
! tsfc=ges_tsen(i,j,1,1)- 273.15_r_kind
!* compare ges snow and retr snow, using the larger one
!* this will add snow based on retr snow.
!* the net add snow will be removed from rain at the same place
!* the may not work for NMMB
! r_cleanSnow_WarmTs_threshold = 4.0 !to delete
! if(tsfc < r_cleanSnow_WarmTs_threshold) then ! add snow on cold sfc
! DO k=1,nsig
! snowtemp=snow_3d(i,j,k)
! rain_3d(i,j,k) = ges_qr(j,i,k)
! snow_3d(i,j,k) = ges_qs(j,i,k)
! if(ref_mos_3d(i,j,k) > zero ) then
! snowtemp = MIN(max(snowtemp,ges_qs(j,i,k)),qrlimit)
! snowadd = max(snowtemp - snow_3d(i,j,k),zero)
! snow_3d(i,j,k) = snowtemp
! raintemp=rain_3d(i,j,k) !+ graupel_3d(i,j,k)
! if(raintemp > snowadd ) then
! if(raintemp > 1.0e-6_r_kind) then
! ratio2=1.0_r_kind - snowadd/raintemp
! rain_3d(i,j,k) = rain_3d(i,j,k) * ratio2
!! graupel_3d(i,j,k) = graupel_3d(i,j,k) * ratio2
! endif
! else
! rain_3d(i,j,k) = 0.0_r_kind
!! graupel_3d(i,j,k) = 0.0_r_kind
! endif
! endif
! END DO
! else ! adjust hydrometeors based on maximum reflectivity level
!* get retr hydrometeors on the maximum ref level
!* get bk maximum hydrometeors in the column
max_retrieved_qrqs=snow_3d(i,j,imaxlvl_ref)+rain_3d(i,j,imaxlvl_ref)
max_bk_qrqs=-999.0_r_kind
imaxlvl_modref=-999
max_obs_qrqs_save(i,j)=max_retrieved_qrqs
DO k=1,nsig
if(ges_qr(j,i,k)+ges_qs(j,i,k) > max_bk_qrqs) then
max_bk_qrqs = ges_qr(j,i,k)+ges_qs(j,i,k)
imaxlvl_modref=k
endif
ENDDO
max_bk_qrqs_save(i,j)=max_bk_qrqs
!* if bk max_hydro is larger, get ratio, based on ratio reduce bk qr and qs
!* if no retr max_hydro, what to do???
if( max_bk_qrqs > max_retrieved_qrqs) then ! tune background hyhro
ratio_hyd_bk2obs=max(min(max_retrieved_qrqs/max_bk_qrqs,1.0_r_kind),0.80_r_kind)
ratio_qrqs_save(i,j)=ratio_hyd_bk2obs
do k=1,nsig
!! graupel_3d(i,j,k) = ges_qg(j,i,k)
! rain_3d(i,j,k) = ges_qr(j,i,k)
!! nrain_3d(i,j,k)= ges_qnr(j,i,k)
! snow_3d(i,j,k) = ges_qs(j,i,k)
if(ges_qr(j,i,k) > qrlimit/2.5_r_kind) then !* based on qr retr curve
! rain_3d(i,j,k) = ges_qr(j,i,k)*ratio_hyd_bk2obs
! rain_3d(i,j,k) = MAX(rain_3d(i,j,k),ges_qr(j,i,k))
!! nrain_3d(i,j,k)= ges_qnr(j,i,k)*ratio_hyd_bk2obs
rain_3d(i,j,k) = MAX(MIN(rain_3d(i,j,k),qrlimit),ges_qr(j,i,k)*ratio_hyd_bk2obs) ! do we need qrlimit?
else
rain_3d(i,j,k) = ges_qr(j,i,k)
endif
if(ges_qs(j,i,k) > qrlimit/3.0_r_kind) then !* based on qs retr curve
! snow_3d(i,j,k) = ges_qs(j,i,k)*ratio_hyd_bk2obs
! snow_3d(i,j,k) = MAX(snow_3d(i,j,k),ges_qs(j,i,k))
snow_3d(i,j,k) = MAX(MIN(snow_3d(i,j,k),qrlimit),ges_qs(j,i,k)*ratio_hyd_bk2obs)
else
snow_3d(i,j,k) = ges_qs(j,i,k)
end if
enddo
else ! use hydro in max refl level
!* if bk max_hydro is smaller, just use retr hydr at max ref levl
DO k=1,nsig
! graupel_3d(i,j,k) = ges_qg(j,i,k)
if(k==imaxlvl_ref) then
snow_3d(i,j,k) = MAX(MIN(snow_3d(i,j,k),qrlimit),ges_qs(j,i,k))
rain_3d(i,j,k) = MAX(MIN(rain_3d(i,j,k),qrlimit),ges_qr(j,i,k)) ! do we need qrlimit?
! snow_3d(i,j,k) = MIN(snow_3d(i,j,k),qrlimit)
! rain_3d(i,j,k) = MIN(rain_3d(i,j,k),qrlimit) ! do we need qrlimit?
! nrain_3d(i,j,k) = nrain_3d(i,j,k)
else
snow_3d(i,j,k) = MAX(MIN(snow_3d(i,j,k),qrlimit)*1.0_r_kind,ges_qs(j,i,k))
rain_3d(i,j,k) = MAX(MIN(rain_3d(i,j,k),qrlimit)*1.0_r_kind,ges_qr(j,i,k)) ! do we need qrlimit?
! rain_3d(i,j,k) = ges_qr(j,i,k)
! snow_3d(i,j,k) = ges_qs(j,i,k)
! nrain_3d(i,j,k) = ges_qnr(j,i,k)
endif
! else
! !use retr qr and qs
! snow_3d(i,j,k) = MIN(snow_3d(i,j,k),qrlimit)
! rain_3d(i,j,k) = MIN(rain_3d(i,j,k),qrlimit) ! do we need qrlimit?
! endif
END DO
end if
else ! clean if ref=0 or use background hydrometeors
!* !if( refmax > 0 .and. (imaxlvl_ref > 0 .and. imaxlvl_ref < nsig ) ) then ! use retrieval hybrometeors
DO k=1,nsig
rain_3d(i,j,k) = ges_qr(j,i,k)
!! nrain_3d(i,j,k)= ges_qnr(j,i,k)
snow_3d(i,j,k) = ges_qs(j,i,k)
!! graupel_3d(i,j,k) = ges_qg(j,i,k)
!* using obsref -99 value
if((iclean_hydro_withRef==1)) then
if( iclean_hydro_withRef_allcol==1 .and. &
(refmax <= zero .and. refmax >= -100_r_kind) .and. &
(sat_ctp(i,j) >=1010.0_r_kind .and. sat_ctp(i,j) <1050._r_kind)) then
rain_3d(i,j,k) = rain_3d(i,j,k)*0.80_r_single
! nrain_3d(i,j,k)= zero
snow_3d(i,j,k) = snow_3d(i,j,k)*0.80_r_single
! graupel_3d(i,j,k) = zero
else
if((ref_mos_3d(i,j,k) <= zero .and. &
ref_mos_3d(i,j,k) > -100.0_r_kind)) then
rain_3d(i,j,k) = rain_3d(i,j,k)*0.80_r_single
! nrain_3d(i,j,k)= zero
snow_3d(i,j,k) = snow_3d(i,j,k)*0.80_r_single
! graupel_3d(i,j,k) = zero
endif
endif
endif
END DO
endif
END DO
END DO
endif !* use retr_qr retr_qs block
!
! remove any negative hydrometeor mixing ratio values
!
where(rain_3d<0.0_r_single) rain_3d=0.0_r_single
where(snow_3d<0.0_r_single) snow_3d=0.0_r_single
DO k=1,nsig
DO j=2,lat2-1
DO i=2,lon2-1
cldwater_3d(i,j,k)= max(0.0_r_single,cldwater_3d(i,j,k))
cldice_3d(i,j,k) = max(0.0_r_single,cldice_3d(i,j,k))
rain_3d(i,j,k) = max(0.0_r_single,rain_3d(i,j,k))
snow_3d(i,j,k) = max(0.0_r_single,snow_3d(i,j,k))
graupel_3d(i,j,k) = max(0.0_r_single,graupel_3d(i,j,k))
END DO
END DO
END DO
! l_conserve_thetaV = .true.
! i_conserve_thetaV_iternum = 3
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)
989 format("sliu debug1::",4e15.4)
988 format("ratio_hyd::",4i8,3e15.4)
! call check_cloud(mype,lat2,lon2,nsig,q_bk,rain_3d,snow_3d,graupel_3d, &
! cldwater_3d,cldice_3d,t_bk,p_bk,h_bk, &
! numsao,NVARCLD_P,numsao,OI,OJ,OCLD,OWX,Oelvtn,cstation, &
! sat_ctp,cld_cover_3d)
!----------------------------------------------
! 6. save the analysis results
!----------------------------------------------
do k=1,nsig
do j=1,lat2
do i=1,lon2
ges_q(j,i,k)=q_bk(i,j,k)/(1_r_kind+q_bk(i,j,k)) ! Here q is mixing ratio kg/kg,
! need to convert to specific humidity
! t_inc=min(t_bk(i,j,k)/((100.0/ges_prsl(j,i,k,1))**rd_over_cp)-ges_tsen(j,i,k,1),0.3)
! ges_tsen(j,i,k,1)=ges_tsen(j,i,k,1)+t_inc
ges_tsen(j,i,k,1)= t_bk(i,j,k)/((100.0_r_kind/ges_prsl(j,i,k,1))**rd_over_cp)
dfi_tten(j,i,k)=0.0_r_kind
!* total heat=ges_tten*20min*60second, for each step: total_heat/(3600 filter window)*(23 time step)
!* for this case: ges_tten*2400/3600*23=ges_tten*(10 to 15)
if(ges_tten(j,i,k)<1.0_r_kind.and. ges_tten(j,i,k)>-0.0000001_r_kind) then
dfi_tten(j,i,k)=ges_tten(j,i,k) !*20.0_r_kind
if(dfi_tten(j,i,k)>0.0025_r_kind) dfi_tten(j,i,k)=0.0025_r_kind
end if
if(xland(i,j)==0.0_r_single)dfi_tten(j,i,k)=dfi_tten(j,i,k)*0.2_r_kind
if(dfi_tten(j,i,k)<0.0_r_kind)dfi_tten(j,i,k)=0.0_r_kind
ges_qr(j,i,k)=rain_3d(i,j,k)
ges_qs(j,i,k)=snow_3d(i,j,k)
ges_ql(j,i,k)=cldwater_3d(i,j,k)
ges_qi(j,i,k)=cldice_3d(i,j,k)
! t=ges_tsen(j,i,k,1)
qfi=ges_qi(j,i,k)
qfs=ges_qs(j,i,k)
qfr=ges_qr(j,i,k)
qfw=ges_ql(j,i,k)
ges_qi(j,i,k)=qfi+qfs
ges_qg(j,i,k)=qfi+qfs+qfr+qfw
ENDDO
ENDDO
ENDDO
!
!----------------------------------------------
! 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,sumq,sumqci,pblh)
deallocate(xlon,xlat,xland,soiltbk)
deallocate(cldwater_3d,cldice_3d,rain_3d,snow_3d,graupel_3d,cldtmp_3d)
if(istat_Surface == 1 ) then
deallocate(OI,OJ,OCLD,OWX,Oelvtn,Odist,cstation)
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==0) then
write(6,*) '========================================'
write(6,*) 'gsdcloudanalysis: generalized cloud analysis finished:',mype
write(6,*) '========================================'
endif
END SUBROUTINE gsdcloudanalysis4NMMB
#else /* Start no NNMB cloud analysis library block */
SUBROUTINE gsdcloudanalysis4NMMB(mype)
use kinds, only: i_kind
implicit none
integer(i_kind),intent(in) :: mype
if(mype==0) then
write(*,*) 'dummy subroutine gsdcloudanalysis4NMMB'
endif
END SUBROUTINE gsdcloudanalysis4NMMB
#endif /* End no NNMB cloud analysis library block */