subroutine map_ctp (ib,jb,nx,ny,nn_obs,numsao,data_s,sat_ctp,sat_tem,w_frac) ! !$$$ subprogram documentation block ! . . . . ! subprogram: map_ctp map GOES cloud product to analysis grid ! ! PRGMMR: Ming Hu ORG: GSD/AMB DATE: 2006-03_10 ! ! ABSTRACT: ! This subroutine map GOES cloud product to analysis grid ! ! PROGRAM HISTORY LOG: ! 2009-01-20 Hu Add NCO document block ! ! ! input argument list: ! ib - begin i point of this domain ! jb - begin j point of this domain ! nx - no. of lons on subdomain (buffer points on ends) ! ny - no. of lats on subdomain (buffer points on ends) ! nn_obs - 1st dimension of observation arry data_s ! numsao - number of observation ! data_s - observation array for GOES cloud products ! ! output argument list: ! sat_ctp - GOES cloud top pressure in analysis grid ! sat_tem - GOES cloud top temperature in analysis grid ! w_frac - GOES cloud coverage in analysis grid ! ! USAGE: ! INPUT FILES: ! ! OUTPUT FILES: ! ! REMARKS: ! ! adapted according to RUC subroutine rd_cld ! * ! * This routine reads NESDIS (Madison, WI) cloud product produced ! * from GOES sounder data. The original product is reprocessed onto ! * MAPS40 grid boxes. There could be more than one cloud product ! * in a grid-box, so we use the nearest one that falls in the ! * grid. The routine combines GOES-8 and 10 products. ! ! ===== History ===== ! ! * Internal variables: ! CTP_E, CTP_W Soft-linked filename for ascii GOES Clouds ! ! * Working variables: ! ! * Working variables used for sorting max size of 10: ! Pxx, Txx, xdist,xxxdist (R4) ! Fxx, Nxx, index, jndex (I4) ! ioption (I4) = 1 if selection is nearest neighbor ! = 2 if selection is median of samples ! ! ! * Output variables on gridpoint (Nx,Ny): ! sat_ctp, sat_tem (R4) Cloud-top pressure and temperature ! w_frac (R4) Effective fractional cloud coverage, option=1 ! fractional coverage within RUC grid, option=2 ! w_eca (R4) Effective fractional cloud regardless option ! (effective cloud amount - eca) ! nlev_cld (I4) Number of cloud levels. TO BE USED LATER ! to incorporate multi-level cloud ! ! * Calling routines ! sorting ! sortmed ! ! * ! ! ATTRIBUTES: ! LANGUAGE: FORTRAN 90 ! MACHINE: Linux cluster (WJET) ! !$$$ ! !_____________________________________________________________________ ! use kinds, only: r_kind,r_single,i_kind use constants, only: zero,one_tenth,one,deg2rad implicit none ! input-file variables: INTEGER(i_kind),intent(in) :: Nx, Ny INTEGER(i_kind),intent(in) :: ib, jb INTEGER(i_kind),intent(in) :: numsao, nn_obs real(r_kind),dimension(nn_obs,numsao):: data_s ! Output real(r_single), intent(out) :: sat_ctp(Nx,Ny) real(r_single), intent(out) :: sat_tem(Nx,Ny) real(r_single), intent(out) :: w_frac(Nx,Ny) ! ! misc integer(i_kind) :: nfov parameter (nfov=60) character header*80 ! Working real(r_kind) :: Pxx(Nx,Ny,nfov),Txx(Nx,Ny,nfov) real(r_kind) :: xdist(Nx,Ny,nfov), xxxdist(nfov) real(r_kind) :: fr,sqrt, qc, type integer(i_kind) :: Nxx(Nx,Ny,nfov),index(Nx,Ny), jndex(nfov) integer(i_kind) :: ioption integer(i_kind) :: ipt,ixx,ii,jj,i,med_pt,igrid,jgrid & ,ncount,ncount1,ncount2,ii1,jj1,nobs,n real(r_kind) :: xc real(r_kind) :: yc real(r_single) :: w_eca(Nx,Ny) integer(i_kind) :: nlev_cld(Nx,Ny) integer(i_kind) :: ios ! ! * Initialize outputs since GOES sounder do not scan all MAPS domain ! do jj=1,Ny do ii=1,Nx w_eca (ii,jj) =-99999._r_kind index(ii,jj) = 0 enddo enddo ! -- set ios as failed unless valid data points are found below ios = 0 ! ----------------------------------------------------------- ! ----------------------------------------------------------- ! Map each FOV onto RR grid points ! ----------------------------------------------------------- ! ----------------------------------------------------------- do ipt=1,numsao xc=data_s(2,ipt) - ib + 1.0_r_kind yc=data_s(3,ipt) - jb + 1.0_r_kind if(data_s(8,ipt) > 50 ) cycle ! * XC,YC should be within subdomain boundary, i.e., XC,YC >0 if(XC >= 1._r_kind .and. XC < Nx .and. & YC >= 1._r_kind .and. YC < Ny) then ii1 = int(xc+0.5_r_kind) jj1 = int(yc+0.5_r_kind) do jj = max(1,jj1-2), min(ny,jj1+2) if (jj1-1 >= 1 .and. jj1+1 <= ny) then do ii = max(1,ii1-2), min(nx,ii1+2) if (ii1-1 >= 1 .and. ii1+1 <= nx) then ! * We check multiple data within gridbox if (index(ii,jj) < nfov) then index(ii,jj) = index(ii,jj) + 1 Pxx(ii,jj,index(ii,jj)) = data_s(4,ipt) Txx(ii,jj,index(ii,jj)) = data_s(6,ipt) !mhu Nxx(ii,jj,index(ii,jj)) = int(data_s(5,ipt)) !mhu no cloud amount available, assign to 100 Nxx(ii,jj,index(ii,jj)) = 100 nlev_cld(ii,jj) = 1 xdist(ii,jj,index(ii,jj)) = sqrt( & (XC+1-ii)**2 + (YC+1-jj)**2) end if endif enddo ! ii endif enddo ! jj endif ! observation is in the domain enddo ! ipt ! ! * ioption = 1 is nearest neighrhood ! * ioption = 2 is median of cloudy fov ioption = 2 ! do jj = 1,Ny do ii = 1,Nx if (index(ii,jj) < 3 ) then ! sat_ctp(ii,jj) = Pxx(ii,jj,1) ! sat_tem(ii,jj) = Txx(ii,jj,1) ! w_frac(ii,jj) = float(Nxx(ii,jj,1))/100. ! w_eca(ii,jj) = float(Nxx(ii,jj,1))/100. elseif(index(ii,jj) >= 3) then ! * We decided to use nearest neighborhood for ECA values, ! * a kind of convective signal from GOES platform... do i=1,index(ii,jj) jndex(i) = i xxxdist(i) = xdist(ii,jj,i) enddo call sorting(xxxdist,index(ii,jj),jndex) w_eca(ii,jj) = float(Nxx(ii,jj,jndex(1)))/100._r_kind ! * Sort to find closest distance if more than one sample if(ioption == 1) then !nearest neighborhood do i=1,index(ii,jj) jndex(i) = i xxxdist(i) = xdist(ii,jj,i) enddo call sorting(xxxdist,index(ii,jj),jndex) sat_ctp(ii,jj) = Pxx(ii,jj,jndex(1)) sat_tem(ii,jj) = Txx(ii,jj,jndex(1)) w_frac(ii,jj) = float(Nxx(ii,jj,jndex(1)))/100._r_kind endif ! * Sort to find median value if(ioption == 2) then !pick median do i=1,index(ii,jj) jndex(i) = i xxxdist(i) = Pxx(ii,jj,i) enddo call sortmed(xxxdist,index(ii,jj),jndex,fr) med_pt = index(ii,jj)/2 + 1 sat_ctp(ii,jj) = Pxx(ii,jj,jndex(med_pt)) sat_tem(ii,jj) = Txx(ii,jj,jndex(med_pt)) w_frac(ii,jj) = fr endif endif enddo !ii enddo !jj return end subroutine map_ctp subroutine sorting(d,n,is) use kinds, only: r_kind,i_kind implicit none integer(i_kind), intent(in) :: n real(r_kind) , intent(inout) :: d(n) integer(i_kind), intent(inout) :: is(n) ! integer(i_kind) :: nm1,ip1,iold,i,j real(r_kind) :: temp ! ! nm1 = n-1 do 10 i=1,nm1 ip1 = i+1 do 10 j=ip1,n if(d(i) <= d(j)) goto 10 temp = d(i) d(i) = d(j) d(j) = temp iold = is(i) is(i) = is(j) is(j) = iold 10 continue return end subroutine sorting subroutine sortmed(p,n,is,f) use kinds, only: r_kind,i_kind implicit none real(r_kind), intent(inout) :: p(n) integer(i_kind), intent(in) :: n integer(i_kind), intent(inout) :: is(n) ! * count cloudy fov real(r_kind), intent(out) :: f integer(i_kind) :: cfov ! integer(i_kind) :: i,j,nm1,ip1,iold real(r_kind) :: temp ! ! ! cfov = 0 do i=1,n if(p(i) < 999._r_kind) cfov = cfov + 1 enddo f = float(cfov)/(max(1,n)) ! cloud-top pressure is sorted high cld to clear nm1 = n-1 do 10 i=1,nm1 ip1 = i+1 do 10 j=ip1,n if(p(i)<=p(j)) goto 10 temp = p(i) p(i) = p(j) p(j) = temp iold = is(i) is(i) = is(j) is(j) = iold 10 continue return end subroutine sortmed