SUBROUTINE cloudCover_Surface(mype,nlat,nlon,nsig,r_radius,thunderRadius,& t_bk,p_bk,q,h_bk,zh, & mxst_p,NVARCLD_P,numsao,OI,OJ,OCLD,OWX,Oelvtn,Odist,& cld_cover_3d,cld_type_3d,wthr_type,pcp_type_3d, & watericemax, kwatericemax) ! !$$$ subprogram documentation block ! . . . . ! subprogram: cloudCover_Surface cloud cover analysis using surface observation ! ! PRGMMR: Ming Hu ORG: GSD/AMB DATE: 2006-10-30 ! ! ABSTRACT: ! This subroutine determines 3D cloud fractional cover using surface observations ! Code based on RUC assimilation code (hybfront/hybcloud.f) ! ! PROGRAM HISTORY LOG: ! 2009-01-20 Hu Add NCO document block ! ! ! input argument list: ! mype - processor ID ! nlon - no. of lons on subdomain (buffer points on ends) ! nlat - no. of lats on subdomain (buffer points on ends) ! nsig - no. of levels ! r_radius - influence radius of the cloud observation ! thunderRadius - ! ! t_bk - 3D background potentional temperature (K) ! p_bk - 3D background pressure (hPa) ! q - 3D moisture (water vapor mixing ratio) ! h_bk - 3D background height (m) ! zh - terrain (m) ! ! mxst_p - maximum observation number ! NVARCLD_P - first dimension of OLCD ! numsao - observation number ! OI - observation x location ! OJ - observation y location ! OLCD - cloud amount, cloud height, visibility ! OWX - weather observation ! Oelvtn - observation elevation ! Odist - distance from the nearest station ! ! output argument list: ! cld_cover_3d- 3D cloud cover ! cld_type_3d - 3D cloud type ! wthr_type - 3D weather type ! pcp_type_3d - 3D weather precipitation type ! ! USAGE: ! INPUT FILES: ! ! OUTPUT FILES: ! ! ! REMARKS: ! ! ATTRIBUTES: ! LANGUAGE: FORTRAN 90 ! MACHINE: Linux cluster (WJET) ! !$$$ ! !_____________________________________________________________________ ! use kinds, only: r_single,i_kind,r_kind implicit none integer(i_kind),intent(in) :: mype REAL(r_single), intent(in) :: r_radius integer(i_kind),intent(in) :: nlat,nlon,nsig real(r_single), intent(in) :: thunderRadius ! ! surface observation ! INTEGER(i_kind),intent(in) :: mxst_p,NVARCLD_P ! PARAMETER (LSTAID_P=9) INTEGER(i_kind),intent(in) :: numsao real(r_single), intent(in) :: OI(mxst_p) ! x location real(r_single), intent(in) :: OJ(mxst_p) ! y location INTEGER(i_kind),intent(in) :: OCLD(NVARCLD_P,mxst_p) ! cloud amount, cloud height, ! visibility CHARACTER*10, intent(in) :: OWX(mxst_p) ! weather real(r_single), intent(in) :: Oelvtn(mxst_p) ! elevation real(r_single), intent(in) :: Odist(mxst_p) ! distance from the nearest station ! ! background ! real(r_single),intent(in) :: t_bk(nlon,nlat,nsig) ! temperature real(r_single),intent(in) :: p_bk(nlon,nlat,nsig) ! pressure real(r_single),intent(in) :: zh(nlon,nlat) ! terrain real(r_single),intent(in) :: q(nlon,nlat,nsig) ! moisture, water vapor mixing ratio (kg/kg) real(r_single),intent(in) :: h_bk(nlon,nlat,nsig) ! height ! REAL(r_single),intent(in) :: watericemax(mxst_p) ! max of background total liquid water in station INTEGER(i_kind),intent(in):: kwatericemax(nlon,nlat) ! lowest level of background total liquid water in grid ! ! Variables for cloud analysis ! real (r_single),intent(inout) :: cld_cover_3d(nlon,nlat,nsig) integer(i_kind),intent(inout) :: cld_type_3d(nlon,nlat,nsig) integer(i_kind),intent(inout) :: wthr_type(nlon,nlat) integer(i_kind),intent(inout) :: pcp_type_3d(nlon,nlat,nsig) ! ! local ! real (r_single) :: vis2qc(nlon,nlat) real (r_single) :: cloud_zthick_p data cloud_zthick_p /300._r_kind/ ! REAL (r_kind) :: spval_p PARAMETER ( spval_p = 99999.0_r_kind ) INTEGER(i_kind) :: i,j,k,k1 INTEGER(i_kind) :: i1,j1,ic INTEGER(i_kind) :: nx_p, ny_p, nztn_p INTEGER(i_kind) :: ista INTEGER(i_kind) :: ich, iob,job REAL(r_kind) :: min_dist, dist REAL(r_kind) :: zdiff REAL(r_kind) :: zlev_clr,cloud_dz,cl_base_ista,betav ! ! ! real(r_single):: tbk_k(nlon,nlat,nsig) real(r_single):: cv_bk(nlon,nlat,nsig) real(r_single):: z_lcl(nlon,nlat) REAL(r_kind) :: cf_model_base,t_model_base, ht_base REAL(r_kind) :: t_dry_adiabat,t_inversion_strength LOGICAL :: l_cf,l_inversion LOGICAL :: if_cloud_exist integer(i_kind) :: firstcloud,cl_base_broken_k real(r_single) :: underlim integer(i_kind) :: npts_near_clr !==================================================================== ! Begin ! ! set constant names consistent with original RUC code ! nx_p=nlon ny_p=nlat nztn_p=nsig vis2qc=-9999.0_r_kind npts_near_clr=0 ! ! !***************************************************************** ! analysis of surface/METAR cloud observations ! ***************************************************************** DO ista=1,numsao i1 = int(oi(ista)+0.0001_r_kind) j1 = int(oj(ista)+0.0001_r_kind) min_dist = Odist(ista) !mh - grid point has the closest cloud station ! -- find out if any precip is present do ich=1,1 if ( owx(ista)(ich:ich+1)=='SH' ) wthr_type(i1,j1)=16 if ( owx(ista)(ich:ich+1)=='TH' .and. & min_dist < thunderRadius) wthr_type(i1,j1)=1 if ( owx(ista)(ich:ich+1)=='RA' ) wthr_type(i1,j1)=11 if ( owx(ista)(ich:ich+1)=='SN' ) wthr_type(i1,j1)=12 if ( owx(ista)(ich:ich+1)=='PL' ) wthr_type(i1,j1)=13 if ( owx(ista)(ich:ich+1)=='DZ' ) wthr_type(i1,j1)=14 if ( owx(ista)(ich:ich+1)=='UP' ) wthr_type(i1,j1)=15 if ( owx(ista)(ich:ich+1)=='BR' ) wthr_type(i1,j1)=21 if ( owx(ista)(ich:ich+1)=='FG' ) wthr_type(i1,j1)=22 enddo ! Consider clear condition case ! ----------------------------- if (ocld(1,ista)==0) then do ic=1,6 if(float(abs(ocld(6+ic,ista))) < 55555) then write(6,*) 'cloudCover_Surface: Observed cloud above the clear level !!!' write(6,*) 'cloudCover_Surface: some thing is wrong in surface cloud observation !' write(6,*) 'cloudCover_Surface: check the station no.', ista, 'at process ', mype write(6,*) ic,OI(ista),OJ(ista) write(6,*) (ocld(k,ista),k=1,12) call stop2(114) endif enddo ! clean the whole column up to ceilometer height (12 kft) if ob is CLR ! h_bk is AGL, not ASL (per Ming Hu's notes below ! ! zlev_clr = Oelvtn(ista)+3650. ! Upcoming mods commented out below for this commit - 4/3/2010 ! PH: added in column cleaning up to ceilometer height if ob is CLR zlev_clr = 3650. do k=1,nztn_p if (h_bk(i1,j1,k) < zlev_clr) then cld_cover_3d(i1,j1,k)=0.0_r_kind pcp_type_3d(i1,j1,k)=0 endif end do wthr_type(i1,j1)=0 ! -- Now consider non-clear obs ! -------------------------- else ! increase zthick by 1.5x factor for ceiling < 900 m (~3000 ft - MVFR) cloud_dz = cloud_zthick_p cl_base_broken_k = -9 ! ????? check with Stan O(h_p) if (Oelvtn(ista).lt.900.) cloud_dz = cloud_zthick_p * 2 do ic = 1,6 if (ocld(ic,ista)>0 .and. ocld(ic,ista)<50) then ! if ( csza(i,j)>=0.10 .and. sat_ctp(i1,j1)>1010.0 & ! .and. sat_ctp(i1,j1)<1050.) go to 1850 ! ! New tweak - 11/07/2009 ! If there was cloud in background over station but if there ! was partial cloudiness within volume and this is one of the ! clear columns within the polygonal area for this METAR, ! then leave it that way and skip. ! if (watericemax(iob,job).gt.0. .and. ! 1 kwatericemax(iob,job).gt.0 .and. ! 1 kwatericemax(iob,job).le.12) then ! npts_cld_match = npts_cld_match + 1 ! dzbase = cl_base_ista - g3(iob,job,kwatericemax(iob,job),h_p) ! sum_dzbase = sum_dzbase + dzbase ! sum_dzbase_abs = sum_dzbase_abs + abs(dzbase) ! end if ! if (watericemax(ista) > 0._r_single .and. kwatericemax(i1,j1)==-1) then ! npts_near_clr = npts_near_clr + 1 ! cycle ! skip cloud build at point (i,j) because background is clear ! end if if(ocld(ic,ista) == 4) then if(wthr_type(i1,j1) > 10 .and. wthr_type(i1,j1) < 20) cloud_dz = 1000._r_kind ! precipitation + highest level if(wthr_type(i1,j1) == 1) cloud_dz = 10000._r_kind ! thunderstorm endif ! --- calculate cloud ceiling level, not exactly, FEW SCT are also considered now ! iob = int(oi(ista)-idw+0.5) ! job = int(oj(ista)-ids+0.5) ! cl_base_ista = (float(ocld(6+ic,ista))+zh(iob,job)) ! cl_base_ista = (float(ocld(6+ic,ista))+Oelvtn(ista)) ! the h_bk is AGL. So observation cloud base should be AGL too, delete Oelvtn(ista) ! cover cloud base observation from AGL to ASL cl_base_ista = float(ocld(6+ic,ista)) + Oelvtn(ista) - zh(i1,j1) if(zh(i1,j1) < 1.0_r_kind .and. Oelvtn(ista) > 20.0_r_kind & .and. float(ocld(6+ic,ista)) < 250.0_r_kind) then cycle ! limit the use of METAR station over oceas for low cloud base endif firstcloud = 0 underlim = 10._r_kind ! do k=1,nztn_p zdiff = cl_base_ista - h_bk(i1,j1,k) ! Must be within cloud_dz meters (300 or 1000 currently) ! ------------------------------------------------------------------- ! -- Bring in the clouds if model level is within 10m under cloud level. if(k==1) underlim=(h_bk(i1,j1,k+1)-h_bk(i1,j1,k))*0.5_r_kind if(k==2) underlim=10.0_r_kind ! 100 feet if(k==3) underlim=20.0_r_kind ! 300 feet if(k==4) underlim=15.0_r_kind ! 500 feet if(k==5) underlim=33.0_r_kind ! 1000 feet if (k>=6 .and. k <= 7) underlim = (h_bk(i1,j1,k+1)-h_bk(i1,j1,k))*0.6_r_kind if(k==8) underlim=95.0_r_kind ! 3000 feet if(k>=9 .and. k= 1.0 .and. (firstcloud==0 .or. abs(zdiff) 10 .and. wthr_type(i1,j1) < 20) then ! cld_type_3d(i1,j1,k)=5 pcp_type_3d(i1,j1,k)=1 endif else write(6,*) 'cloudCover_Surface: wrong cloud coverage observation!' call stop2(114) endif firstcloud = firstcloud + 1 end if ! zdiff < cloud_dz else ! ---- Clear up to cloud base of first cloud level if (ic==1) cld_cover_3d(i1,j1,k)=0 if (ocld(ic,ista) == 1) pcp_type_3d(i1,j1,k)=0 if (ocld(ic,ista) == 3 .or. ocld(ic,ista) == 4) then if( (wthr_type(i1,j1) > 10 .and. wthr_type(i1,j1) < 20) & .or. wthr_type(i1,j1) == 1 ) then pcp_type_3d(i1,j1,k)=1 endif endif end if ! underlim end do ! end K loop ! ----clean cloud above stratusphere do k=1,nztn_p if( h_bk(i1,j1,k) > 18000 ) cld_cover_3d(i1,j1,k)=0 enddo ! end if ! end if ocld > 0 end do ! end IC loop ! ! clean up to broken (3) or if cloud cover less than 2, clean to cloud top ! if(cl_base_broken_k > 0 .and. cl_base_broken_k < nztn_p) then do k=1, cl_base_broken_k if( cld_cover_3d(i1,j1,k) < -0.001_r_kind ) cld_cover_3d(i1,j1,k)=0 enddo else if(ocld(1,ista) == 1 .or. ocld(1,ista) == 2 ) then do k=1, nztn_p if( cld_cover_3d(i1,j1,k) < -0.001_r_kind ) cld_cover_3d(i1,j1,k)=0 enddo endif endif end if ! end if cloudy ob ocld(1,ista) > 0 ! -- Use visibility for low-level cloud whether if (wthr_type(i1,j1) < 30 .and. wthr_type(i1,j1) > 20 .and. & ocld(13,ista) < 5000 .and. ocld(13,ista) > 1) then cld_type_3d(i1,j1,1) = 2 cld_type_3d(i1,j1,2) = 2 betav = 3.912_r_kind / (float(ocld(13,ista)) / 1000._r_kind) vis2qc(i1,j1) = ( (betav/144.7_r_kind) ** 1.14_r_kind) / 1000._r_kind endif ! cloud or clear ENDDO ! ista ! Determine if the layer is dry or it has inversion. ! (in either case, the cloud will be cleared out) ! IF(.false.) THEN ! Set inversion strength flag call BckgrndCC(nlon,nlat,nsig, & t_bk,p_bk,q,h_bk,zh, & cv_bk,tbk_k,z_lcl) ! out DO j = 1,nlat DO i = 1,nlon if_cloud_exist=.false. do k=nsig-1,2,-1 if(cld_cover_3d(i,j,k) > 0.01_r_kind) then cf_model_base = cv_bk(i,j,k) t_model_base = tbk_k(i,j,k) ht_base=h_bk(i,j,k) if_cloud_exist=.true. endif enddo ! ! note, do we need to consider cloud base from background if(if_cloud_exist) then do k=2, nsig-1 if(cld_cover_3d(i,j,k) > 0.01_r_kind) then l_cf=.false. l_inversion=.false. t_dry_adiabat = tbk_k(i,j,2) -.0098_r_kind * (h_bk(i,j,k) - h_bk(i,j,2)) t_inversion_strength = tbk_k(i,j,k) - t_dry_adiabat IF( (tbk_k(i,j,k) > t_model_base) .and. & (tbk_k(i,j,k) > 283.15_r_kind) .and. & ! temp check (t_inversion_strength > 4._r_kind) ) then ! delta theta chk l_inversion = .true. ! Inversion exists endif IF( (cv_bk(i,j,k) < cf_model_base - 0.3_r_kind) .and. & (h_bk(i,j,k) - ht_base >= 500._r_kind) ) THEN l_cf = .true. ! Dry layer exists ENDIF if(l_inversion) then cld_cover_3d(i,j,k) =0.0_r_kind endif endif ! in cloud enddo ! k endif ! if_cloud_exist = true ENDDO ! i ENDDO ! j END IF ! .true. for dry-inversion check. END SUBROUTINE cloudCover_Surface