module satthin !$$$ 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 ! 2007-05-01 wu - correct error in subroutine makegvals which incorrectly ! defines longitude range on regional grid when domain ! includes north pole. ! 2009-04-21 derber - add ithin to call to makegrids - account for negative ithin ! 2009-08-19 guo - added assertions of ntguessig and ntguessfc. ! 2009-09-14 guo - added an experimenting description of the usecase. ! - added an an array size assertion on istart_val(:). ! 2011-04-01 li - add getnst to read nst fields, add destroy_nst ! 2011-05-26 todling - add create_nst ! 2012-01-31 hchuang - add read_nemsnst in sub getnst ! 2012-03-05 akella - remove create_nst,getnst and destroy_nst; nst fields now handled by gsi_nstcoupler ! 2015-05-01 li - modify to use single precision for the variables read from sfc files ! 2016-08-18 li - tic591: when use_readin_anl_sfcmask is true, ! add read sili_anl from analysis grid/resolution sfc file (sfcf06_anl) ! modify to use isli_anl ! determine sno2 with interpolate, accordingly ! use the modified 2d interpolation (sfc_interpolate to intrp22) ! ! 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 ! ! Usecase destription: ! read_obs --> read_airs, etc ! []_makegvals - set up for superob weighting ! []_getsfc - create full horizontal fields of surface arrays ! []_makegrids - set up thinning grids ! []_map2tgrid - map observation to location on thinning grid ! []_checkob - intermediate ob checking to see if it should not be used ! []_finalcheck - the final criterion check for sat obs and increments counters ! combine_radobs - ! []_destroygrids - deallocate thinning grid arrays ! []_destroy_sfc - deallocate full horizontal fields of surface arrays ! ! 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 flag true - no obs in this thinning box previously ! false - previous ob in 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 isli_anl - snow/land/ice mask mask at analysis grid resolution ! def sno_anl - snow-ice mask at analysis grid resolution ! 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,r_quad,r_single use mpeu_util, only: die, perr implicit none ! set default to private private ! set subroutines to public public :: makegvals public :: makegrids public :: getsfc public :: map2tgrid public :: destroygrids public :: destroy_sfc public :: indexx ! set passed variables to public public :: rlat_min,rlon_min,dlat_grid,dlon_grid,superp,super_val1,super_val public :: veg_type_full,soil_type_full,sfc_rough_full,sno_full,sst_full public :: fact10_full,isli_full,soil_moi_full,veg_frac_full,soil_temp_full public :: isli_anl,sno_anl public :: checkob,score_crit,itxmax,finalcheck,zs_full_gfs,zs_full integer(i_kind) mlat,superp,maxthin,itxmax integer(i_kind), save:: itx_all integer(i_kind),dimension(0:51):: istart_val integer(i_kind),allocatable,dimension(:):: mlon logical,allocatable,dimension(:)::icount real(r_kind) rlat_min,rlat_max,rlon_min,rlon_max,dlat_grid,dlon_grid real(r_kind), allocatable, dimension(:) :: glat,score_crit real(r_kind), allocatable, dimension(:) :: super_val,super_val1 real(r_kind), allocatable, dimension(:,:) :: glon,hll real(r_kind), allocatable, dimension(:,:) :: zs_full ! declare the dummy variables of routine read_gfssfc real(r_single), allocatable, dimension(:,:,:) :: fact10_full,sst_full,sno_full real(r_single), allocatable, dimension(:,:) :: veg_type_full real(r_single), allocatable, dimension(:,:,:) :: veg_frac_full real(r_single), allocatable, dimension(:,:) :: soil_type_full real(r_single), allocatable, dimension(:,:,:) :: soil_temp_full,soil_moi_full integer(i_kind),allocatable, dimension(:,:) :: isli_full real(r_single), allocatable, dimension(:,:,:) :: sfc_rough_full real(r_single), allocatable, dimension(:,:) :: zs_full_gfs ! declare the dummy variables of routine read_gfssfc_anl integer(i_kind),allocatable, dimension(:,:) :: isli_anl ! declare local array sno_anl real(r_single),allocatable, dimension(:,:,:) :: sno_anl logical use_all contains subroutine makegvals !$$$ 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 ! 2007-05-01 wu - correct error which incorrectly defines longitude range ! on regional grid when domain includes north pole. ! 2008-05-23 safford - rm unused vars ! 2008-09-08 lueken - merged ed's changes into q1fy09 code ! 2012-10-11 eliu/wu - make sure dlon_e is in the range of 0 and 360 ! 2013-01-09 collard - simplify regional dlon_e range check ! ! 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,& rad2deg,r1000 use obsmod, only: dmesh,dthin,ndat use gridmod, only: regional,nlat,nlon,txy2ll use mpeu_util, only: die 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 integer(i_kind) i,ii,j integer(i_kind) mlonx,icnt,mlony,mlonj real(r_kind) delat,dgv,dx,dy real(r_kind) twopi,dlon_g,dlat_g,dlon_e,dlat_e real(r_kind) factor,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 if (dlon_e < zero) then dlon_e = MOD(dlon_e,-r360) + r360 else if (dlon_e > r360) then dlon_e = MOD(dlon_e,r360) endif 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 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 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,ithin) !$$$ 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 ! 2008-05-23 safford - rm unused vars ! 2008-09-08 lueken - merged ed's changes into q1fy09 code ! 2015-03-23 zaizhong ma - changed itxmax=1e9 for Himawari-8 ahi read in ! ! 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,deg2rad,zero,half,one,two,pi implicit none real(r_kind) ,intent(in ) :: rmesh integer(i_kind),intent(in ) :: ithin real(r_kind),parameter:: r360 = 360.0_r_kind integer(i_kind) i,j integer(i_kind) mlonx,mlonj real(r_kind) dgv,halfpi,dx,dy real(r_kind) twopi real(r_kind) factor,delon real(r_kind) rkm2dg,glatm real(r_quad) delat ! If there is to be no thinning, simply return to calling routine use_all=.false. itx_all=0 if(abs(rmesh) <= one .or. ithin <= 0)then use_all=.true. itxmax=1e9 allocate(icount(itxmax)) allocate(score_crit(itxmax)) do j=1,itxmax icount(j) = .true. score_crit(j) = 9.99e10_r_kind end do 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 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(score_crit(itxmax)) do j=1,itxmax icount(j) = .true. score_crit(j) = 9.99e10_r_kind end do return end subroutine makegrids subroutine getsfc(mype,mype_io,use_sfc,use_sfc_any) !$$$ 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 ! 2008-12-05 todling - add ESMF blocks; update bias_tskin ! 2009-01-28 todling - remove reference to original GMAO interface ! 2010-04-01 treadon - move strip to gridmod ! 2013-10-19 todling - metguess now holds background ! 2013-10-25 todling - reposition ltosi and others to commvars ! 2014-10-05 todling - revisit bkg bias-tskin; rename bkg-bias-related interface ! 2014-12-03 derber - modify reading of surface fields ! 2015-05-01 li - modify to handle the single precision sfc fields read from sfc file ! 2017-08-31 li - modify to read a combined sfc & nst file ! (1) move gsi_nstcoupler_init and gsi_nstcoupler_read from read_obs.F90 to getsfc here ! (2) use sfcnst_comb from name list ! (3) modify subroutine getsfc to read a sfc & nst combined file ! ! input argument list: ! mype - current processor ! mype_io - surface IO processor ! use_sfc - true if processor uses extra surface fields ! use_sfc_any - true if any processor uses extra surface fields ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm rs/6000 sp ! !$$$ use kinds, only: r_kind,r_single use gridmod, only: nlat,nlon,lat2,lon2,lat1,lon1,jstart,& iglobal,itotsub,ijn,displs_g,regional,istart, & rlats,rlons,nlat_sfc,nlon_sfc,rlats_sfc,rlons_sfc,strip,& sfcnst_comb,use_gfs_nemsio,use_readin_anl_sfcmask use general_commvars_mod, only: ltosi,ltosj use guess_grids, only: ntguessig,isli,sfct,sno,fact10, & nfldsfc,ntguessfc,soil_moi,soil_temp,veg_type,soil_type, & veg_frac,sfc_rough,nfldsig,isli2,sno2 use m_gsiBiases, only: bkg_bias_model,bias_hour use jfunc, only: bcoption use mpimod, only: mpi_comm_world,ierror,mpi_rtype,mpi_rtype4 use constants, only: zero,half,pi,two,one use ncepgfs_io, only: read_gfssfc,read_gfssfc_anl use ncepnems_io, only: read_nemssfc,intrp22,read_nemssfc_anl use sfcio_module, only: sfcio_realfill use obsmod, only: lobserver use gsi_nstcouplermod, only: nst_gsi,gsi_nstcoupler_init,gsi_nstcoupler_read use gsi_nstcouplermod, only: tref_full,dt_cool_full,z_c_full,dt_warm_full,z_w_full,& c_0_full,c_d_full,w_0_full,w_d_full use gsi_metguess_mod, only: gsi_metguess_bundle use gsi_bundlemod, only: gsi_bundlegetpointer implicit none integer(i_kind),intent(in ) :: mype,mype_io logical ,intent(in ) :: use_sfc,use_sfc_any ! 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):: bias real(r_kind),dimension(nlat):: ailoc real(r_kind),dimension(nlon):: ajloc real(r_kind),allocatable,dimension(:)::wlatx,slatx real(r_kind) :: dlon, missing real(r_single),allocatable,dimension(:,:)::dum,work integer(i_kind) mm1,i,j,k,it,il,jl,jmax,idrt,istatus character(24) filename real(r_kind),dimension(:,:),pointer:: ges_z =>NULL() if( (ntguessig<1.or.ntguessig>nfldsig) .or. & (ntguessfc<1.or.ntguessfc>nfldsfc) ) then call perr('satthin.getsfc','ntguessig = ',ntguessig) call perr('satthin.getsfc','ntguessfc = ',ntguessfc) call die('satthin.getsfc') endif mm1=mype+1 ! if(mype == 0)write(6,*)'GETSFC: enter with nlat_sfc,nlon_sfc=',nlat_sfc,nlon_sfc,& ! ' and nlat,nlon=',nlat,nlon if(regional)then nlat_sfc=nlat nlon_sfc=nlon end if ! if(mype == 0)write(6,*)'GETSFC: set nlat_sfc,nlon_sfc=',nlat_sfc,nlon_sfc allocate(rlats_sfc(nlat_sfc),rlons_sfc(nlon_sfc)) allocate(isli_full(nlat_sfc,nlon_sfc),fact10_full(nlat_sfc,nlon_sfc,nfldsfc)) allocate(sst_full(nlat_sfc,nlon_sfc,nfldsfc),sno_full(nlat_sfc,nlon_sfc,nfldsfc)) allocate(zs_full(nlat,nlon)) allocate(sfc_rough_full(nlat_sfc,nlon_sfc,nfldsfc)) allocate(isli_anl(nlat,nlon)) allocate(sno_anl(nlat,nlon,nfldsfc)) allocate(soil_moi_full(nlat_sfc,nlon_sfc,nfldsfc),soil_temp_full(nlat_sfc,nlon_sfc,nfldsfc)) allocate(veg_frac_full(nlat_sfc,nlon_sfc,nfldsfc),soil_type_full(nlat_sfc,nlon_sfc)) allocate(veg_type_full(nlat_sfc,nlon_sfc)) do j=1,lon1*lat1 zsm(j)=zero end do ! Create full horizontal nst arrays if (nst_gsi > 0) call gsi_nstcoupler_init() ! Global read #ifndef HAVE_ESMF if (.not. regional) then if (nlon == nlon_sfc .and. nlat == nlat_sfc) then rlats_sfc=rlats rlons_sfc=rlons else idrt=4 jmax=nlat_sfc-2 allocate(slatx(jmax),wlatx(jmax)) call splat(idrt,jmax,slatx,wlatx) dlon=two*pi/float(nlon_sfc) do i=1,nlon_sfc rlons_sfc(i)=float(i-1)*dlon end do do i=1,(nlat_sfc-1)/2 rlats_sfc(i+1)=-asin(slatx(i)) rlats_sfc(nlat_sfc-i)=asin(slatx(i)) end do rlats_sfc(1)=-half*pi rlats_sfc(nlat_sfc)=half*pi deallocate(slatx,wlatx) end if allocate(zs_full_gfs(nlat_sfc,nlon_sfc)) if ( use_gfs_nemsio ) then if ( sfcnst_comb .and. nst_gsi > 0 ) then call read_nemssfc(mype_io, & sst_full,soil_moi_full,sno_full,soil_temp_full, & veg_frac_full,fact10_full,sfc_rough_full, & veg_type_full,soil_type_full,zs_full_gfs,isli_full,use_sfc_any,& tref_full,dt_cool_full,z_c_full,dt_warm_full,z_w_full,c_0_full,c_d_full,w_0_full,w_d_full) else call read_nemssfc(mype_io, & sst_full,soil_moi_full,sno_full,soil_temp_full, & veg_frac_full,fact10_full,sfc_rough_full, & veg_type_full,soil_type_full,zs_full_gfs,isli_full,use_sfc_any) endif ! if ( nst_gsi > 0 ) then if ( use_readin_anl_sfcmask ) then call read_nemssfc_anl(mype_io,isli_anl) endif else call read_gfssfc(mype_io, & sst_full,soil_moi_full,sno_full,soil_temp_full, & veg_frac_full,fact10_full,sfc_rough_full, & veg_type_full,soil_type_full,zs_full_gfs,isli_full,use_sfc_any) if ( use_readin_anl_sfcmask ) then call read_gfssfc_anl(mype_io,isli_anl) endif endif ! ! read NSST variables while .not. sfcnst_comb (in sigio or nemsio) ! if (nst_gsi > 0 .and. .not. sfcnst_comb) then call gsi_nstcoupler_read(mype_io) ! Read NST fields (each proc needs full NST fields) endif if(.not. use_sfc)then deallocate(soil_moi_full,soil_temp_full) deallocate(veg_frac_full,soil_type_full) deallocate(veg_type_full) end if if (bcoption>0) then if (mype==0) write(6,*)'GETSFC: add bias correction to guess field ', trim(filename) ! Correct Tskin over the full grid if(nlon == nlon_sfc .and. nlat == nlat_sfc)then call bkg_bias_model(work2,'sst',bias_hour,ierror) if (ierror==0) then ! successful application of bkg bias model ... call strip(work2,zsm) 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) bias(i,j)=nint(work1(k)) end do do it=1,nfldsfc do j=1,nlon do i=1,nlat sst_full(i,j,it)=sst_full(i,j,it)+bias(i,j) end do end do end do end if ! check for successful application of bkg bias model else write(6,*)'satthin bias correction - incompatible surface resolution',& nlon,nlon_sfc,nlat,nlat_sfc call stop2(82) end if end if else ! for regional #endif /* HAVE_ESMF */ it=ntguessfc rlats_sfc=rlats rlons_sfc=rlons ! isli_full do j=1,lon2 do i=1,lat2 work2(i,j)=isli(i,j,it) end do end do call strip(work2,zsm) 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) isli_full(i,j)=nint(work1(k)) end do ! Fields with multiple time levels do it=1,nfldsfc ! sst_full call strip(sfct(:,:,it),zsm) 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,it)=work1(k) end do ! fact10_full call strip(fact10(:,:,it),zsm) 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,it)=work1(k) end do ! sfc_rough_full call strip(sfc_rough(:,:,it),zsm) 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) sfc_rough_full(i,j,it)=work1(k) end do ! sno_full call strip(sno(:,:,it),zsm) 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,it)=work1(k) end do ! veg_frac_full call strip(veg_frac(:,:,it),zsm) call mpi_allgatherv(zsm,ijn(mm1),mpi_rtype,& work1,ijn,displs_g,mpi_rtype,& mpi_comm_world,ierror) if(use_sfc)then do k=1,iglobal i=ltosi(k) ; j=ltosj(k) veg_frac_full(i,j,it)=work1(k) end do end if ! soil_temp_full call strip(soil_temp(:,:,it),zsm) call mpi_allgatherv(zsm,ijn(mm1),mpi_rtype,& work1,ijn,displs_g,mpi_rtype,& mpi_comm_world,ierror) if(use_sfc)then do k=1,iglobal i=ltosi(k) ; j=ltosj(k) soil_temp_full(i,j,it)=work1(k) end do end if ! soil_moi_full call strip(soil_moi(:,:,it),zsm) call mpi_allgatherv(zsm,ijn(mm1),mpi_rtype,& work1,ijn,displs_g,mpi_rtype,& mpi_comm_world,ierror) if(use_sfc)then do k=1,iglobal i=ltosi(k) ; j=ltosj(k) soil_moi_full(i,j,it)=work1(k) end do end if end do !end do ntguessfc ! Now single time level surface fields it=ntguessfc ! soil_type_full call strip(soil_type(:,:,it),zsm) call mpi_allgatherv(zsm,ijn(mm1),mpi_rtype,& work1,ijn,displs_g,mpi_rtype,& mpi_comm_world,ierror) if(use_sfc)then do k=1,iglobal i=ltosi(k) ; j=ltosj(k) soil_type_full(i,j)=work1(k) end do end if ! veg_type_full call strip(veg_type(:,:,it),zsm) call mpi_allgatherv(zsm,ijn(mm1),mpi_rtype,& work1,ijn,displs_g,mpi_rtype,& mpi_comm_world,ierror) if(use_sfc)then do k=1,iglobal i=ltosi(k) ; j=ltosj(k) veg_type_full(i,j)=work1(k) end do end if #ifndef HAVE_ESMF end if ! if (.not. regional) then #endif /* HAVE_ESMF */ ! Now stuff that isn't model dependent ! zs_full it=ntguessig call gsi_bundlegetpointer (gsi_metguess_bundle(it),'z',ges_z,istatus) if (istatus==0) then do j=1,lon1*lat1 zsm(j)=zero end do call strip(ges_z,zsm) 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 endif ! old gfs surface files do not have a terrain record. ! therefore, need to interpolate the terrain from the atm ! grid to the surface grid, which may not be the same. if (.not. regional) then missing = sfcio_realfill + one if (maxval(zs_full_gfs) < missing) then if (nlon == nlon_sfc .and. nlat == nlat_sfc) then zs_full_gfs = zs_full else allocate(dum(nlat_sfc,nlon_sfc)) allocate(work(nlat,nlon)) work = zs_full call intrp22(work,rlons,rlats,nlon,nlat, & dum, rlons_sfc,rlats_sfc,nlon_sfc,nlat_sfc) zs_full_gfs = dum deallocate(dum) deallocate(work) endif endif endif ! find subdomain for isli2 if (nlon == nlon_sfc .and. nlat == nlat_sfc) then do j=1,lon2 jl=j+jstart(mm1)-2 jl=min0(max0(1,jl),nlon) do i=1,lat2 il=i+istart(mm1)-2 il=min0(max0(1,il),nlat) isli2(i,j)=isli_full(il,jl) do k=1,nfldsfc sno2(i,j,k)=sno_full(il,jl,k) end do end do end do else if ( use_readin_anl_sfcmask ) then do k = 1, nfldsfc call intrp22(sno_full(:,:,k),rlons_sfc,rlats_sfc,nlon_sfc,nlat_sfc, & sno_anl (:,:,k),rlons,rlats,nlon,nlat) enddo do j=1,lon2 jl=j+jstart(mm1)-2 jl=min0(max0(1,jl),nlon) do i=1,lat2 il=i+istart(mm1)-2 il=min0(max0(1,il),nlat) isli2(i,j)=isli_anl(il,jl) do k=1,nfldsfc sno2(i,j,k)=sno_anl(il,jl,k) if ( isli2(i,j) == 0 ) then sno2(i,j,k) = zero endif end do end do end do else ailoc=rlats ajloc=rlons call grdcrd(ailoc,nlat,rlats_sfc,nlat_sfc,1) call grdcrd(ajloc,nlon,rlons_sfc,nlon_sfc,1) do j=1,lon2 jl=j+jstart(mm1)-2 jl=min0(max0(1,jl),nlon) jl=nint(ajloc(jl)) jl=min0(max0(1,jl),nlon_sfc) do i=1,lat2 il=i+istart(mm1)-2 il=min0(max0(1,il),nlat) il=nint(ailoc(il)) il=min0(max0(1,il),nlat_sfc) isli2(i,j)=isli_full(il,jl) do k=1,nfldsfc sno2(i,j,k) =sno_full(il,jl,k) end do end do end do end if end if if (.not.lobserver) then if(allocated(veg_frac)) deallocate(veg_frac) if(allocated(veg_type)) deallocate(veg_type) if(allocated(soil_type)) deallocate(soil_type) if(allocated(soil_moi)) deallocate(soil_moi) if(allocated(sfc_rough)) deallocate(sfc_rough) endif return end subroutine getsfc subroutine map2tgrid(dlat_earth,dlon_earth,dist1,crit1,itx,ithin,itt,iuse,sis) !$$$ 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. ! 2013-01-26 parrish - change from grdcrd to grdcrd1 (to allow successful debug compile on WCOSS) ! ! 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 ! sis - sensor/instrument/satellite ! ! output argument list: ! itx - combined (i,j) index of observation on thinning grid ! itt - superobs thinning counter ! iuse - .true. if observation should be used ! ! ! attributes: ! language: f90 ! machine: ibm rs/6000 sp ! !$$$ 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 character(20) ,intent(in ) :: sis 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 .or. ithin <= 0)then iuse=.true. itt=1 dist1=one if(itx_all < itxmax) then itx_all=itx_all+1 else iuse = .false. write(6,*)'MAP2TGRID: ndata > maxobs when reading data for ',sis,itxmax end if itx=itx_all 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 grdcrd1(dlat1,glat,mlat,1) iy=int(dlat1) dy=dlat1-iy iy=max(1,min(iy,mlat)) call grdcrd1(dlon1,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))iuse=.false. return end subroutine map2tgrid subroutine checkob(dist1,crit1,itx,iuse) !$$$ 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 ! !$$$ 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))return if(crit1*dist1 > score_crit(itx))iuse=.false. return end subroutine checkob subroutine finalcheck(dist1,crit1,itx,iuse) !$$$ 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 ! ! output argument list: ! iuse - .true. if observation should be used ! ! ! attributes: ! language: f90 ! machine: ibm rs/6000 sp ! !$$$ implicit none logical ,intent(inout) :: iuse integer(i_kind),intent(in ) :: itx real(r_kind) ,intent(in ) :: dist1,crit1 real(r_kind) crit ! If not using data, simply return to calling routine if(.not. iuse)return crit=crit1*dist1 ! Increment obs counter if(icount(itx) .or. crit < score_crit(itx))then score_crit(itx)= crit icount(itx)=.false. else iuse = .false. end if return end subroutine finalcheck subroutine destroygrids !$$$ 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 deallocate(icount) deallocate(score_crit) ! 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) 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 if(allocated(sst_full))deallocate(sst_full) if(allocated(sno_full))deallocate(sno_full) if(allocated(fact10_full))deallocate(fact10_full) if(allocated(isli_full))deallocate(isli_full) if(allocated(veg_type_full))deallocate(veg_type_full) if(allocated(soil_type_full))deallocate(soil_type_full) if(allocated(veg_frac_full))deallocate(veg_frac_full) if(allocated(soil_temp_full))deallocate(soil_temp_full) if(allocated(soil_moi_full))deallocate(soil_moi_full) if(allocated(zs_full))deallocate(zs_full) if(allocated(sfc_rough_full))deallocate(sfc_rough_full) if(allocated(zs_full_gfs)) deallocate(zs_full_gfs) return end subroutine destroy_sfc subroutine indexx(n,arr,indx) !$$$ 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_double use constants, only: one implicit none integer(i_kind) maxblock parameter (maxblock=1000) integer(i_kind),intent(in ) :: n real(r_kind) ,intent(in ) :: arr(n) integer(i_kind),intent( out) :: indx(n) #ifdef ibm_sp real(r_kind),dimension(maxblock)::work real(r_double) oned oned=1._r_double 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 loop0: do if(ir-l<m)then loop: do j=l+1,ir indxt=indx(j) a=arr(indxt) do i=j-1,l,-1 if(arr(indx(i))<=a)then indx(i+1)=indxt cycle loop end if indx(i+1)=indx(i) end do i=l-1 indx(i+1)=indxt end do loop if(jstack==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))>arr(indx(ir)))then itemp=indx(l) indx(l)=indx(ir) indx(ir)=itemp endif if(arr(indx(l+1))>arr(indx(ir)))then itemp=indx(l+1) indx(l+1)=indx(ir) indx(ir)=itemp endif if(arr(indx(l))>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) loop1: do i=i+1 if(arr(indx(i))<a)cycle loop1 loop2: do j=j-1 if(arr(indx(j))<=a)exit loop2 end do loop2 if(j<i)exit loop1 itemp=indx(i) indx(i)=indx(j) indx(j)=itemp end do loop1 indx(l+1)=indx(j) indx(j)=indxt jstack=jstack+2 if(jstack>nstack)then write(6,*)'INDEXX: nstack=',nstack,' too small in indexx' call stop2(32) endif if(ir-i+1>=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 end do loop0 #endif end subroutine indexx end module satthin