module satthin 16,1 !$$$ subprogram documenation block ! . . . . ! subprogram: satthin ! prgmmr: treadon org: np20 date: 2002-10-17 ! ! abstract: This module contains code which may be used to selectively ! thin satellite data. ! ! program history log: ! 2002-10-17 treadon ! 2004-05-28 kleist, subroutine call update ! 2005-03-28 wu, remove unused variable mlath ! 2005-10-06 treadon - add routine destroy_sfc ! 2005-11-29 parrish - add routines getsfc_nmm_binary, getsfc_nmm_netcdf, ! getsfc_mass_binary, getsfc_mass_netcdf ! 2005-12-08 treadon - rename global getsfc to getsfc_global, add getsfc ! to branch to appropriate getsfc* routine ! 2006-02-01 parrish - remove all getsfc routines, and create new getsfc ! which makes full sfc fields from guess_grids. ! rename surface fields to sst_full,sno_full, etc. ! to eliminate conflict with some fields with same ! names in guess_grids. ! also add new sfc fields u10_full,v10_full ! 2006-05-23 parrish - add new sfc field zs_full, the model terrain height ! 2006-07-28 derber - use r1000 from constants ! ! Subroutines Included: ! sub makegvals - set up for superob weighting ! sub makegrids - set up thinning grids ! sub getsfc - create full horizontal fields of surface arrays ! sub map2tgrid - map observation to location on thinning grid ! sub destroygrids - deallocate thinning grid arrays ! sub destroy_sfc - deallocate full horizontal surface arrays ! sub indexx - sort array into ascending order ! ! Variable Definitions: ! def mlat - number of latitudes in thinning grid ! def mlon - number of longitudes in thinning grid ! def superp - maximum number of data types ! def maxthin - maximum number of obs to retain in thinning grid box ! def itxmax - total number of points in thinning grid ! def istart_val - starting location on thinning grid for superobs ! def icount - observation counter ! def ibest_obs - index of "best" quality obs in thinning grid box ! def isli_full - snow/land/ice mask ! def glat - latitudes on thinning grid ! def super_val - super obs factor across data types ! def super_val1 - super obs factors summed over all mpi tasks (data types) ! def glon - longitudes on thinning grid ! def hll - (i,j) index of point on thinning grid ! def sli_full - 0=sea/1=land/2=ice mask ! def sst_full - skin temperature ! def sno_full - snow-ice mask ! def u10_full - 10 meter zonal wind ! def v10_full - 10 meter meridional wind ! def zs_full - model terrain elevation ! def score_crit - "best" quality obs score in thinning grid box ! def use_all - parameter for turning satellite thinning algorithm off ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ end documentation block use kinds, only: r_kind,i_kind implicit none integer(i_kind) mlat,superp,maxthin,itxmax integer(i_kind),dimension(0:51):: istart_val integer(i_kind),allocatable,dimension(:):: mlon,icount,ibest_obs integer(i_kind),allocatable,dimension(:,:):: isli_full real(r_kind) rlat_min,rlat_max,rlon_min,rlon_max,dlat_grid,dlon_grid real(r_kind),allocatable,dimension(:):: glat real(r_kind),allocatable,dimension(:):: super_val,super_val1 real(r_kind),allocatable,dimension(:,:):: glon,hll,sli_full,sst_full,sno_full real(r_kind),allocatable,dimension(:,:):: u10_full,v10_full,zs_full real(r_kind),allocatable,dimension(:):: score_crit logical use_all contains subroutine makegvals 1,4 !$$$ subprogram documentation block ! . . . . ! subprogram: makegvals ! prgmmr: derber org: np23 date: 2002-10-17 ! ! abstract: This routine allocates and initializes arrays ! used for superobs weighting. ! ! program history log: ! 2002-10-17 derber ! 2004-06-22 treadon - update documentation ! 2004-12-09 treadon - allocate thinning grids consistent with analysis domain ! 2006-07-28 derber - use r1000 from constants ! ! input argument list: ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm rs/6000 sp ! !$$$ use constants, only: deg2rad,rearth_equator,zero,two,pi,half,one,quarter,& rad2deg,r1000 use obsmod, only: dmesh,dthin,ndat use gridmod, only: regional,region_lat,region_lon,nlat,nlon,txy2ll implicit none real(r_kind),parameter:: r90 = 90.0_r_kind real(r_kind),parameter:: r360 = 360.0_r_kind real(r_kind),parameter:: r999 = 999.0_r_kind logical odd integer(i_kind) i,ii,j,k integer(i_kind) mlonx,icnt,mlony,mlonj real(r_kind) delonx,delat,dgv,dx,dy real(r_kind) twopi,dlon_g,dlat_g,dlon_e,dlat_e real(r_kind) factor,factors,delon real(r_kind) rkm2dg,glatm,glatx ! Initialize variables, set constants maxthin=0 do i=1,ndat maxthin=max(maxthin,abs(dthin(i))) end do istart_val=0 twopi = two*pi rkm2dg = r360/(twopi*rearth_equator)*r1000 ! Set lat,lon limits for analysis grid. rlat_min = -r90 rlat_max = r90 rlon_min = zero rlon_max = r360 if (regional) then rlat_min = r999 rlat_max = -r999 rlon_min = r999 rlon_max = -r999 do j=1,nlon dlon_g=j do i=1,nlat dlat_g=i call txy2ll(dlon_g,dlat_g,dlon_e,dlat_e) dlat_e=dlat_e*rad2deg dlon_e=dlon_e*rad2deg rlat_min = min(rlat_min,dlat_e) rlat_max = max(rlat_max,dlat_e) rlon_min = min(rlon_min,dlon_e) rlon_max = max(rlon_max,dlon_e) end do end do if (rlon_min < zero) rlon_min = rlon_min + r360 if (rlon_max < zero) rlon_max = rlon_max + r360 endif dlat_grid = rlat_max - rlat_min dlon_grid = rlon_max - rlon_min do ii=1,maxthin ! Set up dimensions for thinning grids istart_val(ii+1)=istart_val(ii) if(abs(dmesh(ii)) > one)then dx = dmesh(ii)*rkm2dg dy = dx mlat = dlat_grid/dy + half mlonx = dlon_grid/dx + half delat = dlat_grid/mlat delonx= dlon_grid/mlonx dgv = delat*half mlat=max(2,mlat); mlonx=max(2,mlonx) icnt=0 do j = 1,mlat glatx = rlat_min + (j-1)*delat glatx = glatx*deg2rad glatm = glatx + dgv*deg2rad factor = abs(cos(abs(glatm))) if (dmesh(ii)>zero) then mlonj = nint(mlonx*factor) mlony = max(2,mlonj) else delon = factor*dmesh(ii) delon = min(delon,r360) mlony = dlon_grid/delon endif do i = 1,mlony icnt=icnt+1 istart_val(ii+1)=istart_val(ii+1)+1 enddo enddo end if end do superp=istart_val(maxthin+1) ! Allocate and initialize arrays for superobs weighthing allocate(super_val(0:superp),super_val1(0:superp)) do i=0,superp super_val(i)=zero end do return end subroutine makegvals subroutine makegrids(rmesh) 9,1 !$$$ subprogram documentation block ! . . . . ! subprogram: makegrids ! prgmmr: treadon org: np23 date: 2002-10-17 ! ! abstract: This routine sets up dimensions for and allocates ! thinning grids. ! ! program history log: ! 2002-10-17 treadon ! 2004-06-22 treadon - update documentation ! 2004-12-09 treadon - allocate thinning grids consistent with analysis domain ! ! input argument list: ! rmesh - mesh size (km) of thinning grid. If (rmesh <= one), ! then no thinning of the data will occur. Instead, ! all data will be used without thinning. ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm rs/6000 sp ! !$$$ use constants, only: rearth_equator,two,deg2rad,zero,half,one,quarter,pi implicit none real(r_kind),parameter:: r360 = 360.0_r_kind logical odd integer(i_kind) i,j,k integer(i_kind) mlonx,mlonj real(r_kind) delonx,delat,dgv,halfpi,dx,dy real(r_kind) twopi real(r_kind) factor,factors,delon,rmesh real(r_kind) rkm2dg,glatm ! If there is to be no thinning, simply return to calling routine use_all=.false. if(abs(rmesh) <= one)then use_all=.true. itxmax=2.e6 return end if ! Set constants halfpi = half*pi twopi = two*pi rkm2dg = r360/(twopi*rearth_equator)*1.e3_r_kind ! Set up dimensions and allocate arrays for thinning grids if (rmesh<zero) rkm2dg=one dx = rmesh*rkm2dg dy = dx mlat = dlat_grid/dy + half mlonx = dlon_grid/dx + half delat = dlat_grid/mlat delonx= dlon_grid/mlonx dgv = delat*half mlat=max(2,mlat); mlonx=max(2,mlonx) allocate(mlon(mlat),glat(mlat),glon(mlonx,mlat),hll(mlonx,mlat)) ! Set up thinning grid lon & lat. The lon & lat represent the location of the ! lower left corner of the thinning grid box. itxmax=0 do j = 1,mlat glat(j) = rlat_min + (j-1)*delat glat(j) = glat(j)*deg2rad glatm = glat(j) + dgv*deg2rad factor = abs(cos(abs(glatm))) if (rmesh>zero) then mlonj = nint(mlonx*factor) mlon(j) = max(2,mlonj) delon = dlon_grid/mlon(j) else delon = factor*rmesh delon = min(delon,r360) mlon(j) = dlon_grid/delon endif glat(j) = min(max(-halfpi,glat(j)),halfpi) do i = 1,mlon(j) itxmax=itxmax+1 hll(i,j)=itxmax glon(i,j) = rlon_min + (i-1)*delon glon(i,j) = glon(i,j)*deg2rad glon(i,j) = min(max(zero,glon(i,j)),twopi) enddo end do ! Allocate and initialize arrays allocate(icount(itxmax)) allocate(ibest_obs(itxmax)) allocate(score_crit(itxmax)) do j=1,itxmax icount(j)=0 end do do j=1,itxmax score_crit(j) =9.99e6_r_kind end do return end subroutine makegrids subroutine getsfc(mype) 1,12 !$$$ subprogram documentation block ! . . . . ! subprogram: getsfc ! prgmmr: parrish org: np23 date: 2006-02-02 ! ! abstract: This routine converts subdomain surface fields in ! guess_grids to full horizontal fields for use in ! reading of observations. ! ! program history log: ! 2002-10-17 parrish ! ! input argument list: ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm rs/6000 sp ! !$$$ use kinds, only: r_kind,i_kind use gridmod, only: nlat,nlon,lat2,lon2,lat1,lon1,& ltosi,ltosj,iglobal,itotsub,ijn,displs_g use guess_grids, only: ntguessig,isli,sfct,sno,fact10,ges_u,ges_v,ges_z use mpimod, only: mpi_comm_world,ierror,mpi_rtype,strip use constants, only: zero implicit none integer(i_kind) mype ! Local variables real(r_kind),dimension(lat1*lon1):: zsm real(r_kind),dimension(itotsub):: work1 real(r_kind),dimension(lat2,lon2):: work2 real(r_kind),dimension(nlat,nlon):: fact10_full integer(i_kind) mm1,i,j,k,it mm1=mype+1 it=ntguessig allocate(isli_full(nlat,nlon),sli_full(nlat,nlon)) allocate(sst_full(nlat,nlon),sno_full(nlat,nlon)) allocate(u10_full(nlat,nlon),v10_full(nlat,nlon)) allocate(zs_full(nlat,nlon)) ! isli_full,sli_full do j=1,lon1*lat1 zsm(j)=zero end do do j=1,lon2 do i=1,lat2 work2(i,j)=isli(i,j,it) end do end do call strip(work2,zsm,1) call mpi_allgatherv(zsm,ijn(mm1),mpi_rtype,& work1,ijn,displs_g,mpi_rtype,& mpi_comm_world,ierror) do k=1,iglobal i=ltosi(k) ; j=ltosj(k) sli_full(i,j)=work1(k) isli_full(i,j)=nint(work1(k)) end do ! sst_full do j=1,lon1*lat1 zsm(j)=zero end do call strip(sfct(1,1,it),zsm,1) call mpi_allgatherv(zsm,ijn(mm1),mpi_rtype,& work1,ijn,displs_g,mpi_rtype,& mpi_comm_world,ierror) do k=1,iglobal i=ltosi(k) ; j=ltosj(k) sst_full(i,j)=work1(k) end do ! sno_full do j=1,lon1*lat1 zsm(j)=zero end do call strip(sno(1,1,it),zsm,1) call mpi_allgatherv(zsm,ijn(mm1),mpi_rtype,& work1,ijn,displs_g,mpi_rtype,& mpi_comm_world,ierror) do k=1,iglobal i=ltosi(k) ; j=ltosj(k) sno_full(i,j)=work1(k) end do ! fact10_full do j=1,lon1*lat1 zsm(j)=zero end do call strip(fact10(1,1,it),zsm,1) call mpi_allgatherv(zsm,ijn(mm1),mpi_rtype,& work1,ijn,displs_g,mpi_rtype,& mpi_comm_world,ierror) do k=1,iglobal i=ltosi(k) ; j=ltosj(k) fact10_full(i,j)=work1(k) end do ! u10_full do j=1,lon1*lat1 zsm(j)=zero end do call strip(ges_u(1,1,1,it),zsm,1) call mpi_allgatherv(zsm,ijn(mm1),mpi_rtype,& work1,ijn,displs_g,mpi_rtype,& mpi_comm_world,ierror) do k=1,iglobal i=ltosi(k) ; j=ltosj(k) u10_full(i,j)=work1(k) end do ! v10_full do j=1,lon1*lat1 zsm(j)=zero end do call strip(ges_v(1,1,1,it),zsm,1) call mpi_allgatherv(zsm,ijn(mm1),mpi_rtype,& work1,ijn,displs_g,mpi_rtype,& mpi_comm_world,ierror) do k=1,iglobal i=ltosi(k) ; j=ltosj(k) v10_full(i,j)=work1(k) end do ! zs_full do j=1,lon1*lat1 zsm(j)=zero end do call strip(ges_z(1,1,it),zsm,1) call mpi_allgatherv(zsm,ijn(mm1),mpi_rtype,& work1,ijn,displs_g,mpi_rtype,& mpi_comm_world,ierror) do k=1,iglobal i=ltosi(k) ; j=ltosj(k) zs_full(i,j)=work1(k) end do end subroutine getsfc subroutine map2tgrid(dlat_earth,dlon_earth,dist1,crit1,itx,ithin,itt,iuse) 9,4 !$$$ subprogram documentation block ! . . . . ! subprogram: map2tgrid ! prgmmr: derber org: np2 date: 2006-05-03 ! ! abstract: This routine maps observations to the thinning grid. ! ! program history log: ! 2006-05-03 derber (created from map2grids) ! 2006-09-13 treadon - set itx=1 for the case use_all=.true. ! ! input argument list: ! dlat_earth - earth relative observation latitude (radians) ! dlon_earth - earth relative observation longitude (radians) ! crit1 - quality indicator for observation (smaller = better) ! ithin - number of obs to retain per thinning grid box ! ! output argument list: ! iobs - observation counter ! itx - combined (i,j) index of observation on thinning grid ! itt - superobs thinning counter ! iobsout- location for observation to be put ! iuse - .true. if observation should be used ! ! ! attributes: ! language: f90 ! machine: ibm rs/6000 sp ! !$$$ use kinds, only: r_kind,i_kind use constants, only: one, half implicit none logical,intent(out):: iuse integer(i_kind),intent(in):: ithin integer(i_kind),intent(out):: itt,itx real(r_kind),intent(in):: dlat_earth,dlon_earth,crit1 real(r_kind),intent(out):: dist1 integer(i_kind) ix,iy real(r_kind) dlat1,dlon1,dx,dy,dxx,dyy ! If using all data (no thinning), simply return to calling routine if(use_all)then iuse=.true. itt=0 itx=1 return end if ! Compute (i,j) indices of coarse mesh grid (grid number 1) which ! contains the current observation. dlat1=dlat_earth dlon1=dlon_earth call grdcrd(dlat1,1,glat,mlat,1) iy=int(dlat1) dy=dlat1-iy iy=max(1,min(iy,mlat)) call grdcrd(dlon1,1,glon(1,iy),mlon(iy),1) ix=int(dlon1) dx=dlon1-ix ix=max(1,min(ix,mlon(iy))) dxx=half-min(dx,one-dx) dyy=half-min(dy,one-dy) dist1=dxx*dxx+dyy*dyy+half itx=hll(ix,iy) itt=istart_val(ithin)+itx if(ithin == 0) itt=0 ! Increment obs counter on coarse mesh grid. Also accumulate observation ! score and distance functions ! ratio=1.e9 ! if ( dx > zero ) ratio=dy/dx ! dista=sin(two*atan(ratio)) ! distb=sin(pi*dx) !dista+distb is max at grid box center ! dist1=one - quarter*(dista + distb) !dist1 is min at grid box center and !ranges from 1 (at corners)to !.5 (at center of box) iuse=.true. if(dist1*crit1 > score_crit(itx) .and. icount(itx) == 0)iuse=.false. return end subroutine map2tgrid subroutine checkob(dist1,crit1,itx,iuse) 7,1 !$$$ subprogram documentation block ! . . . . ! subprogram: checkob ! prgmmr: checkob org: np23 date: 2002-10-17 ! ! abstract: intermediate ob checking routine to see if we should not use. ! ! program history log: ! 2006-05-03 derber ! ! input argument list: ! dist1 - quality indicator for distance (smaller = better) ! crit1 - quality indicator for observation (smaller = better) ! itx - combined (i,j) index of observation on thinning grid ! iuse - .true. if observation should be used ! ! output argument list: ! iuse - .true. if observation should be used ! ! ! attributes: ! language: f90 ! machine: ibm rs/6000 sp ! !$$$ use kinds, only: r_kind,i_kind implicit none logical,intent(inout):: iuse integer(i_kind),intent(in):: itx real(r_kind),intent(in):: dist1,crit1 ! If using all data (no thinning), simply return to calling routine if(use_all .or. .not. iuse .or. icount(itx)==0)return if(crit1*dist1 > score_crit(itx))iuse=.false. return end subroutine checkob subroutine finalcheck(dist1,crit1,iobs,itx,iobsout,iuse,sis) 9,1 !$$$ subprogram documentation block ! . . . . ! subprogram: finalcheck ! prgmmr: derber org: np23 date: 2002-10-17 ! ! abstract: This routine performs the final criterion check for sat ! obs and increments counters ! ! program history log: ! 2006-05-03 derber ! ! input argument list: ! dist1 - quality indicator for distance (smaller = better) ! crit1 - quality indicator for observation (smaller = better) ! itx - combined (i,j) index of observation on thinning grid ! iobs - observation counter ! iuse - .true. if observation should be used ! sis - sensor/instrument/satellite ! ! output argument list: ! iobs - observation counter ! iobsout- location for observation to be put ! iuse - .true. if observation should be used ! ! ! attributes: ! language: f90 ! machine: ibm rs/6000 sp ! !$$$ use kinds, only: r_kind,i_kind implicit none logical,intent(inout):: iuse integer(i_kind),intent(inout):: iobs,iobsout integer(i_kind),intent(in):: itx real(r_kind),intent(in):: dist1,crit1 character(20),intent(in):: sis real(r_kind) crit if(.not. iuse)return ! If using all data (no thinning), simply return to calling routine if(use_all)then if(iobs < itxmax)then iobs=iobs+1 iobsout=iobs else iuse = .false. write(6,*)' ndata > maxobs when reading data for ',sis,itxmax end if return end if crit=crit1*dist1 if(icount(itx) == 0)then ! Increment obs counter if(iobs < itxmax)then iobs=iobs+1 iobsout=iobs score_crit(itx)= crit ibest_obs(itx) = iobs icount(itx)=icount(itx)+1 else iuse = .false. write(6,*)' ndata > maxobs when reading data for ',sis,itxmax end if else if(crit < score_crit(itx))then score_crit(itx)= crit iobsout=ibest_obs(itx) icount(itx)=icount(itx)+1 else iuse = .false. end if return end subroutine finalcheck subroutine destroygrids 9 !$$$ subprogram documentation block ! . . . . ! subprogram: destroygrids ! prgmmr: treadon org: np23 date: 2002-10-17 ! ! abstract: This deallocate arrays used in thinning of satellite data. ! ! program history log: ! 2002-10-17 treadon ! 2004-06-22 treadon - update documentation ! ! input argument list: ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm rs/6000 sp ! !$$$ implicit none ! If using all data (no thinning), the arrays following the ! return are never allocated --> therefore nothing to deallocate if(use_all) return deallocate(mlon,glat,glon,hll) deallocate(icount) deallocate(ibest_obs) deallocate(score_crit) return end subroutine destroygrids subroutine destroy_sfc !$$$ subprogram documentation block ! . . . . ! subprogram: destroy_sfc ! prgmmr: treadon org: np23 date: 2005-10-06 ! ! abstract: This deallocate surface arrays ! ! program history log: ! 2005-10-06 treadon ! ! input argument list: ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm rs/6000 sp ! !$$$ implicit none deallocate(sst_full,sli_full,sno_full,isli_full) deallocate(u10_full,v10_full,zs_full) return end subroutine destroy_sfc subroutine indexx(n,arr,indx) 1,2 !$$$ subprogram documentation block ! . . . . ! subprogram: indexx ! prgmmr: treadon org: np23 date: 2002-10-17 ! ! abstract: This routine indexes a sort key array, arr, such that ! arr(indx(i),i=1,n) is the sort key in ascending order. ! ! program history log: ! 2002-10-17 treadon ! 2004-06-22 treadon - update documentation ! ! input argument list: ! n - size of sort key array ! arr - sort key array ! ! output argument list: ! indx - index array ! ! attributes: ! language: f90 ! machine: ibm rs/6000 sp ! !$$$ use kinds, only: r_kind,r_double,i_kind implicit none integer(i_kind) maxblock parameter (maxblock=1000) integer(i_kind) n,indx(n) real(r_kind) arr(n) #ifdef ibm_sp real(r_kind),dimension(maxblock)::work real(r_double) oned real(r_kind) one one=1.; oned=1. if (digits(one)<digits(oned)) then call ssorts(arr,1,n,indx,work,maxblock) else call dsorts(arr,1,n,indx,work,maxblock) endif #else integer(i_kind) m,nstack parameter (m=7,nstack=500) integer(i_kind) i,indxt,ir,itemp,j,jstack,k,l,istack(nstack) real(r_kind) a do j=1,n indx(j)=j end do jstack=0 l=1 ir=n 1 continue if(ir-l.lt.m)then do j=l+1,ir indxt=indx(j) a=arr(indxt) do i=j-1,l,-1 if(arr(indx(i)).le.a)goto 2 indx(i+1)=indx(i) end do i=l-1 2 continue indx(i+1)=indxt end do if(jstack.eq.0)return ir=istack(jstack) l=istack(jstack-1) jstack=jstack-2 else k=(l+ir)/2 itemp=indx(k) indx(k)=indx(l+1) indx(l+1)=itemp if(arr(indx(l)).gt.arr(indx(ir)))then itemp=indx(l) indx(l)=indx(ir) indx(ir)=itemp endif if(arr(indx(l+1)).gt.arr(indx(ir)))then itemp=indx(l+1) indx(l+1)=indx(ir) indx(ir)=itemp endif if(arr(indx(l)).gt.arr(indx(l+1)))then itemp=indx(l) indx(l)=indx(l+1) indx(l+1)=itemp endif i=l+1 j=ir indxt=indx(l+1) a=arr(indxt) 3 continue i=i+1 if(arr(indx(i)).lt.a)goto 3 4 continue j=j-1 if(arr(indx(j)).gt.a)goto 4 if(j.lt.i)goto 5 itemp=indx(i) indx(i)=indx(j) indx(j)=itemp goto 3 5 continue indx(l+1)=indx(j) indx(j)=indxt jstack=jstack+2 if(jstack.gt.nstack)then write(6,*)'INDEXX: nstack=',nstack,' too small in indexx' call stop2(32) endif if(ir-i+1.ge.j-l)then istack(jstack)=ir istack(jstack-1)=i ir=j-1 else istack(jstack)=j-1 istack(jstack-1)=l l=i endif endif goto 1 #endif end subroutine indexx end module satthin