subroutine get_patmos(debug, hdf_file_id, cloud_probability, land_cover, * cloud_phase, avhrr_ctp) c----------------------------------------------------------------------------------- c Description: Acquire PATMOS-X level 2 cloud probability for each AVHRR (GAC) c pixel in the current orbit. c Also get "packed_land_cover", "cloud_phase" SDSs and unpack. c c Input parameters: c debug Debug write flag c hdf_file_id HDF file handle c c Output parameters: c cloud_probability PATMOS probability of cloud for all AVHRR (GAC) pixels c in current orbit. c land_cover PATMOS land cover index (same as in MODIS geolocation c files) c 0 = shallow ocean c 1 = land c 2 = coast c 3 = shallow inland water c 4 = ephemeral water c 5 = deep inland water c 6 = moderate or continental ocean c 7 = deep ocean c cloud_phase PATMOS cloud phase product c 0 = clear c 1 = water c 2 = supercooled water c 3 = mixed ice and water c 4 = ice c 5 = unknown c avhrr_ctp PATMOS cloud top pressure product c c Revision history: 12/10 R. Frey originial version c 03/11 R. Frey added land cover c 12/13 R. Frey added cloud phase c 03/14 R. Frey added cloud top pressure c c Non-HDF calls: none c c---------------------------------------------------------------------------------- IMPLICIT NONE save c---------------------------------------------------------------------------------- include 'hdf.inc' c---------------------------------------------------------------------------------- integer ncp parameter (ncp = 409 * 15000) c---------------------------------------------------------------------------------- c Calling arguments. real cloud_probability(409,15000) real avhrr_ctp(409,15000) integer land_cover(409,15000) integer cloud_phase(409,15000) integer debug, hdf_file_id c Local variables. character*160 errmsg character*48 SDS_name real scf, offset integer sfn2index, sfselect, sfginfo, sfrdata, sffattr, num_type, nattrs, rank integer sfrattr integer sds_indx, pid, irt, dims(3), start(3), stride(3), edge(3), atr_index integer k, j, level, status integer mask, lcv integer*2 ctp(409,15000) byte cp(409,15000) byte plc(409,15000) byte phs(409,15000) c---------------------------------------------------------------------------------------- data cp / ncp*0 / data plc / ncp*0 / c---------------------------------------------------------------------------------------- c Initializations. dims(1) = 0 dims(2) = 0 dims(3) = 0 start(1) = 0 start(2) = 0 start(3) = 0 stride(1) = 1 stride(2) = 1 stride(3) = 1 edge(1) = 0 edge(2) = 0 edge(3) = 0 c---------------------------------------------------------------------------------------- c Read cloud probability values. c Get SDS index and array identifier. sds_indx = sfn2index(hdf_file_id,'cloud_probability') pid = sfselect(hdf_file_id,sds_indx) c Get array dimensions. irt = sfginfo(pid,SDS_name,rank,dims,num_type,nattrs) edge(1) = dims(1) edge(2) = dims(2) edge(3) = dims(3) c Read array. irt = sfrdata(pid,start,stride,edge,cp(1,1)) if(irt .ne. 0) then write(*,'(/,''sds index, pid, dims'',5i10)') sds_indx, pid, dims level = 3 status = -1 write( errmsg,'('' Cannot read GAC cloud probabilities'')') call message( 'get_patmos', errmsg, status, level ) end if c---------------------------------------------------------------------------------------- c Unscale values. do j = 1, 15000 do k = 1, 409 if(cp(k,j) .ne. -128) then cloud_probability(k,j) = cp(k,j) * 0.004 + 0.5 if(cloud_probability(k,j) .gt. 1.0) cloud_probability(k,j) = 1.0 if(cloud_probability(k,j) .lt. 0.0) cloud_probability(k,j) = 0.0 else cloud_probability(k,j) = 9.99 end if enddo enddo c---------------------------------------------------------------------------------------- c Read packed land cover values. c Get SDS index and array identifier. sds_indx = sfn2index(hdf_file_id,'packed_land_cover') pid = sfselect(hdf_file_id,sds_indx) c Get array dimensions. irt = sfginfo(pid,SDS_name,rank,dims,num_type,nattrs) edge(1) = dims(1) edge(2) = dims(2) edge(3) = dims(3) c Read array. irt = sfrdata(pid,start,stride,edge,plc(1,1)) if(irt .ne. 0) then write(*,'(/,''sds index, pid, dims'',5i10)') sds_indx, pid, dims level = 3 status = -1 write( errmsg,'('' Cannot read GAC packed land cover SDS'')') call message( 'get_patmos', errmsg, status, level ) end if c---------------------------------------------------------------------------------------- c Unscale values. do j = 1, 15000 do k = 1, 409 if(plc(k,j) .lt. 0) then lcv = plc(k,j) + 256 else lcv = plc(k,j) end if mask = 224 land_cover(k,j) = ishft( iand(lcv, mask), -5) c write(*,'(''land_cover: '',4i10)') k,j,plc(k,j),land_cover(k,j) enddo enddo c---------------------------------------------------------------------------------------- c Read cloud phase values. c Get SDS index and array identifier. sds_indx = sfn2index(hdf_file_id,'cloud_phase') pid = sfselect(hdf_file_id,sds_indx) c Get array dimensions. irt = sfginfo(pid,SDS_name,rank,dims,num_type,nattrs) edge(1) = dims(1) edge(2) = dims(2) edge(3) = dims(3) c Read array. irt = sfrdata(pid,start,stride,edge,phs(1,1)) if(irt .ne. 0) then write(*,'(/,''sds index, pid, dims'',5i10)') sds_indx, pid, dims level = 3 status = -1 write( errmsg,'('' Cannot read GAC cloud phase SDS'')') call message( 'get_patmos', errmsg, status, level ) end if c---------------------------------------------------------------------------------------- c Transfer values to integer array. do j = 1, 15000 do k = 1, 409 if(phs(k,j) .lt. 0) then cloud_phase(k,j) = 9.99 else cloud_phase(k,j) = phs(k,j) end if enddo enddo c---------------------------------------------------------------------------------------- c Read cloud top pressure values ('cld_press_acha'). c Get SDS index and array identifier. sds_indx = sfn2index(hdf_file_id,'cld_press_acha') pid = sfselect(hdf_file_id,sds_indx) c Get array dimensions. irt = sfginfo(pid,SDS_name,rank,dims,num_type,nattrs) edge(1) = dims(1) edge(2) = dims(2) edge(3) = dims(3) c ... Extract the scale factor and offset for this parameter atr_index = sffattr(pid,'scale_factor') irt = sfrattr(pid,atr_index,scf) if (irt .lt. 0) then write(*,'( 2x,''Unable to extract scale factor for '', + ''SDS cld_press_acha'')') endif atr_index = sffattr(pid,'add_offset') irt = sfrattr(pid,atr_index,offset) if (irt .lt. 0) then write(*,'( 2x,''Unable to extract additive offset for '', + ''SDS cld_press_acha'')') endif c Read array. irt = sfrdata(pid,start,stride,edge,ctp(1,1)) if(irt .ne. 0) then write(*,'(/,''sds index, pid, dims'',5i10)') sds_indx, pid, dims level = 3 status = -1 write( errmsg,'('' Cannot read GAC cloud top pressure SDS'')') call message( 'get_patmos', errmsg, status, level ) end if c---------------------------------------------------------------------------------------- c Unscale the byte values. do j = 1, 15000 do k = 1, 409 if(ctp(k,j) .eq. -32768) then avhrr_ctp(k,j) = -99.99 else avhrr_ctp(k,j) = (ctp(k,j) * scf) + offset end if enddo enddo c---------------------------------------------------------------------------------------- c Transfer values to integer array. if(debug .gt. 0) then write(*,'(''PATMOS L2 file dimensions: '',3i10)') dims write(*,'(/,''Sample AVHRR (GAC) cloud probabilities:'')') do j = 660, 670 write(*,'(6i5,5x,5f6.3)') j,(cp(k,j),k=213,217),(cloud_probability(k,j),k=213,217) enddo write(*,'(/,''Sample AVHRR (GAC) land cover values:'')') do j = 660, 670 write(*,'(6i5,5x,5i5)') j,(plc(k,j),k=213,217),(land_cover(k,j),k=213,217) enddo write(*,'(/,''Sample AVHRR (GAC) cloud phase values:'')') do j = 660, 670 write(*,'(6i5,5x,5i5)') j,(phs(k,j),k=213,217),(cloud_phase(k,j),k=213,217) enddo write(*,'(/,''Sample AVHRR (GAC) cloud top pressure values:'')') do j = 660, 670 write(*,'(6i10,5x,5f10.2)') j,(ctp(k,j),k=213,217),(avhrr_ctp(k,j),k=213,217) enddo end if c---------------------------------------------------------------------------------------- return end