***************************************************************** * PROGRAM: cmap_pen_rt_merging_pen_v0712.lnx.f * * LANGUAGE: Fortran77 * * MACHINE: Any LINUX BOX *----------------------------------------------------------------* * PURPOSE: to merge the IR, MW estimates, gauge * observation and model prediction of * monthly rainfall fields by Maximum * Likelihood Estimation *----------------------------------------------------------------* * USAGE: *----------------------------------------------------------------* * INPUT FILES: Unit 10 == input data availability * 20 == land / sea mask * 50 == individual input data * ch-01 --- GPI estimates * ch-02 --- GPI quality * ch-03 --- MW scattering estimates * ch-04 --- MW scattering quality * ch-05 --- OPI estimates (All Sources) * ch-06 --- OPI quality (All Sources) * ch-07 --- gauge observations * ch-08 --- gauge quality * * OUTPUT FILES: Unit 50 == merged analysis * ch-09 --- comb analysis with all sources * ch-10 --- merg analysis ---------------- * ch-11 --- error field ---------------- * 60 == yearly output file of final merged analysis *----------------------------------------------------------------* * SUBROUTINES USED: include 'CMAP_parameters.h' * * FUNCTIONS USED: (all Standard Fortran 77 functions) * *----------------------------------------------------------------* * INPUT VARIABLES: kyy == year of the target date * kmm == month of the target date * kdd == day of the target date * kjj == Julian Day * * LOCAL VARIABLES: kdata == data availability * rland == land / sea mask * rain == input precip * qual == quality of input precip * err == error of input precip * comb == combined satellites * blend == gauge-satellite merged * error == error of gauge-satellite merged * *----------------------------------------------------------------* * DATE 2007-12-09 *----------------------------------------------------------------* * MODIFICATIONS: ****************************************************************** c=====|==|=========|=========|=========|=========|=========|=========|== c program CMAP_merg_pen_2_5deg c----------------------------------------------------------------------- c nsate : number of satellite data sources c nx : x direction (east-west) dimension c ny : y direction (north-south) dimension c----------------------------------------------------------------------- include 'CMAP_parameters.h' character yyyymmdd*8 integer yyyy, mm, dd, days(12), leap, m, jjj, ddd, pen character*128 f_mask, f_work, f_cmap dimension rland(nx,ny), 1 rain (nx,ny,nsate), 2 qual (nx,ny,nsate), 3 err (nx,ny,nsate), 4 comb (nx,ny), 5 blend(nx,ny), 6 error(nx,ny) dimension err_atl(nsate), qua_atl(nsate) data days / 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 / c----------------------------------------------------------------------- c check parameter c----------------------------------------------------------------------- if( iargc().lt.1 ) then write(*, *) write(*, *) '#################################' write(*, *) '# #' write(*, *) '# parameter needed #' write(*, *) '# #' write(*, *) '#################################' write(*, *) stop endif c----------------------------------------------------------------------- c get target date c----------------------------------------------------------------------- call getarg(1, yyyymmdd) read(yyyymmdd, '(i4.4,2i2.2)') yyyy, mm, dd c----------------------------------------------------------------------- c adjust for leap-year c----------------------------------------------------------------------- if( mod(yyyy,4).eq.0 ) days(2)=29 if( mod(yyyy,100).eq.0 ) days(2)=28 if( mod(yyyy,400).eq.0 ) days(2)=29 if( days(2).eq.29 ) then leap=1 else leap=0 endif c----------------------------------------------------------------------- c calculate Julian day of target date c----------------------------------------------------------------------- jjj=dd do m=1, mm-1 jjj=jjj+days(m) enddo c----------------------------------------------------------------------- c determine the pentad c----------------------------------------------------------------------- if( jjj.le.59 ) then ddd=jjj else ddd=jjj-leap endif if( mod(ddd,5).ne.0 ) then write(*, *) '## NOT end of pentad: ', yyyy, mm, dd, jjj stop endif pen=ddd/5 write(*, *) yyyy, mm, dd, jjj, ' pentad: ', pen c----------------------------------------------------------------------- c construct file names c----------------------------------------------------------------------- c------------>123456789|123456789|123456789|123456789|123456789|1234567<-- f_mask='../library/mask_land_sea_edge.lnx' f_work='../output/CMAP_PEN_RT/yyyy/cmap_pen_rt_v0011.lnx.yyyypp' f_cmap='../output/CMAP_PEN_RT/yyyy/cmap_pen_rt_v0011_out.lnx.yyyy' write(f_work(23:26), '(i4.4)') yyyy write(f_work(50:55), '(i4.4,i2.2)') yyyy, pen write(f_cmap(23:26), '(i4.4)') yyyy write(f_cmap(54:57), '(i4.4)') yyyy c----------------------------------------------------------------------- c read in land mask c----------------------------------------------------------------------- c 0 ==> ocean c 1 ==> land c 2 ==> partial land c 3 ==> land boarding to ocean c 4 ==> ocean boarding to land c----------------------------------------------------------------------- open(10, file=f_mask, status='old', 1 access='direct', recl=nx*ny*4) read(10, rec=1) rland close(10) c----------------------------------------------------------------------- c set Arctic region to ocean c----------------------------------------------------------------------- do jj=70, 72 do ii=1, nx rland(ii,jj)=0.0 enddo enddo c----------------------------------------------------------------------- c set Hawaii to ocean c----------------------------------------------------------------------- rland(82,44)=0.0 c----------------------------------------------------------------------- c open pentad merging working file (subroutines read from unit 50) c----------------------------------------------------------------------- open(50, file=f_work, status='old', 1 access='direct', recl=nx*ny*4) write(*, *) ' <-- ', f_work c----------------------------------------------------------------------- c to define the error structure over tropical ocean c----------------------------------------------------------------------- c ndat : length of error array c ndel : +/- period used to calculate the random error c err_atl : error compared to atolls c----------------------------------------------------------------------- call def_err_atl(err_atl, qua_atl) c----------------------------------------------------------------------- c to combine the 3 satellite fields c----------------------------------------------------------------------- call com_sate(pen, rland, err_atl, qua_atl) c----------------------------------------------------------------------- c to define the bias structure of the combined analysis c over tropical Pacific c----------------------------------------------------------------------- c bias_atl : bias of the combined analysis in % c----------------------------------------------------------------------- call def_bias_atl(bias_atl) c----------------------------------------------------------------------- c to blend with gauge data c----------------------------------------------------------------------- ngauge0=1 call rmv_bias(rland, ngauge0, bias_atl) c----------------------------------------------------------------------- c to output the merged analysis (item 10) into the yearly file c----------------------------------------------------------------------- open(60, file=f_cmap, 1 access='direct', recl=nx*ny*4) read(50, rec=10) blend write(60, rec=pen) blend close(50) close(60) write(*, *) ' ==> ', f_cmap stop end c c=====|==|=========|=========|=========|=========|=========|=========|== c c c ************************************************************************* c *** *** c *** *** c *** ***** ***** ***** ***** ***** ***** * ***** * *** c *** * * * * * * * * * * * * * *** c *** * * ***** ***** ***** ***** ***** ***** * * *** c *** * * * * * * * * * * * * * *** c *** ***** ***** * *** ***** * * * * *** * * * ***** *** c *** *** c ************************************************************************* c subroutine def_err_atl ( err_atl, qua_atl ) c include 'CMAP_parameters.h' c dimension err_atl(nsate), qua_atl(nsate) c c err_atl (1) = 0.445 * 1.946 err_atl (2) = 0.533 * 1.946 err_atl (3) = 0.483 * 1.946 c qua_atl (1) = 1.0 qua_atl (2) = 1.0 qua_atl (3) = 1.0 c return end c c c ************************************************************************* c *** *** c *** *** c *** ***** ***** * * ***** * ***** ***** *** c *** * * * ** ** * * * * * *** c *** * * * * * * ***** ***** * ***** *** c *** * * * * * * * * * * *** c *** ***** ***** * * ***** ***** * * * ***** *** c *** *** c ************************************************************************* c c ************************************************************************* c *** *** c *** objective : *** c *** to combine the satellite estimates and *** c *** model precipitation fields by *** c *** the Maximum Likelihood Estimation *** c *** *** c ************************************************************************* c subroutine com_sate ( kpen, rland, err_atl, qua_atl ) c include 'CMAP_parameters.h' c dimension rland (nx,ny), # err_atl (nsate), # qua_atl (nsate) c dimension rain (nx,ny,nsate+1), # qua_sam (nx,ny,nsate+1), # error (nx,ny,nsate), # qua_alg (nx,ny,nsate), # comb_rn (nx,ny), # comb_er (nx,ny), # globe1 (nx,ny), # globe2 (nx,ny) c dimension error0(nsate), qual0(nsate) c c 1. to read in the individual data sources c call input_data ( rain, qua_sam ) c c 2. to define random error over land c call def_err_land ( rland, rain, error ) c c 3. to define the relative quality functions (algorithm) c kpp = kpen if (kpp.le.6) then kseason = 1 else if (kpp.le.12) then kseason = 2 else if (kpp.le.18) then kseason = 3 else if (kpp.le.24) then kseason = 4 else if (kpp.le.30) then kseason = 5 else if (kpp.le.36) then kseason = 6 else if (kpp.le.42) then kseason = 7 else if (kpp.le.49) then kseason = 8 else if (kpp.le.55) then kseason = 9 else if (kpp.le.61) then kseason = 10 else if (kpp.le.67) then kseason = 11 else kseason = 12 end if call def_quality ( kseason, rland, qua_alg ) c c 4. to define the random error over ocean c do 4001 kk = 1, nsate error0 (kk) = err_atl (kk) qual0 (kk) = qua_atl (kk) 4001 continue call def_err_ocn # ( rland, error0, qual0, qua_sam, qua_alg, error ) c c 5. to fill in the gap in the error fields c call fill_err_gap ( rain, error ) c c 6. to combine through the maximum like estimation (MLE) c call com_by_mle ( rain, error, comb_rn, comb_er ) kout1 = 9 kout2 = 11 write (50,rec=kout1) comb_rn write (50,rec=kout2) comb_er c c return end c c c ************************************************************************* c *** *** c *** *** c *** ***** ***** ***** ***** ***** * ***** * ***** * *** c *** * * * * * * * * * * * * * * *** c *** * * ***** ***** **** * ***** ***** ***** * * *** c *** * * * * * * * * * * * * * * *** c *** ***** ***** * ** ***** ***** * * ***** ** * * * ***** *** c *** *** c ************************************************************************* c c ************************************************************************* c *** *** c *** objective : *** c *** to calculate the bias of the combias analysis *** c *** over ocean *** c *** *** c ************************************************************************* c subroutine def_bias_atl ( bias_atl ) c c c 1. to define the bias based on over all comparison c bias_atl = -0.119 + 0.0144 c c return end c c c ************************************************************************* c *** *** c *** *** c *** ***** * * * * ***** ***** * ***** *** c *** * * ** ** * * * * * * * * *** c *** ***** * * * * * **** * ***** ***** *** c *** * * * * * * * * * * * * *** c *** * * * * * ***** ***** ***** * * ***** *** c *** *** c ************************************************************************* c c ************************************************************************* c *** *** c *** objective : *** c *** to remove the bias from the combined field *** c *** *** c ************************************************************************* c subroutine rmv_bias ( rland, ngauge0, bias_atl ) c include 'CMAP_parameters.h' c dimension rland(nx,ny) c dimension globe (nx,ny), # rain (nx,ny,nsate+1), # r_gauge (nx,ny), # q_gauge (nx,ny), # c_merge (nx,ny), # b_merge (nx,ny), # e_merge (nx,ny) c dimension kcheck(nx,ny),kbench(nx,ny), # eglobe(164,92),wgt(3),err_gag(11) c data wgt / 0.25, 1.00, 0.25 /, # err_gag / 0.60, 0.30, 0.25, 0.20, 0.15, 0.12, # 0.10, 0.08, 0.07, 0.06, 0.05/ c c c ************************************************************************* c *** *** c *** 1. to read in data *** c *** *** c ************************************************************************* c c ## to read in data ## c kinp1 = 7 kinp2 = 8 kinp3 = 9 kinp4 = 11 read (50,rec=kinp1) r_gauge read (50,rec=kinp2) q_gauge read (50,rec=kinp3) c_merge read (50,rec=kinp4) e_merge c c ## to define the data availability ## c do 1002 kk = 1, nsate+1 kinp = (kk-1)*2 + 1 read (50,rec=kinp) globe do 1001 jj = 1, ny do 1001 ii = 1, nx rain (ii,jj,kk) = globe (ii,jj) 1001 continue 1002 continue c nt = nsate + 1 do 1004 jj = 1, ny do 1004 ii = 1, nx kcheck(ii,jj) = 0 do 1003 kk = 1, nt if (rain(ii,jj,kk).ge.0.0) then kcheck(ii,jj) = 1 go to 1004 end if 1003 continue 1004 continue c c ************************************************************************* c *** *** c *** 2. to remove the bias over land area *** c *** *** c ************************************************************************* c call rmv_bias_land ( ngauge0, rland, kcheck, # r_gauge, q_gauge, c_merge, b_merge, kbench ) c c ##### redefine the merged field for the antaritic area ##### c do 2002 jj=1,12 do 2002 ii=1,nx if (rland(ii,jj).ge.1.0.and. # rland(ii,jj).le.3.0.and. # r_gauge(ii,jj).ge.0.0) then ng=nint(q_gauge(ii,jj)) if (ng.gt.10) then ng = 10 end if b_merge(ii,jj) = r_gauge(ii,jj) e_merge(ii,jj) = err_gag(ng+1) end if 2002 continue c c ##### redefine the merged error for benchmark grids ##### c do 2003 jj = 13, ny do 2003 ii = 1, nx if (rland(ii,jj).ge.1.0.and. # rland(ii,jj).le.3.0) then if (kbench(ii,jj).eq.1) then ng=nint(q_gauge(ii,jj)) if (ng.gt.10) then ng=10 end if e_merge(ii,jj)=err_gag(ng+1) end if end if 2003 continue c c ************************************************************************* c *** *** c *** 3. to remove the bias over ocean *** c *** *** c ************************************************************************* c bias0 = bias_atl call rmv_bias_ocean ( rland, bias0, c_merge, b_merge ) c c ************************************************************************* c *** *** c *** 4. to smooth the land/ocean boundaries *** c *** *** c ************************************************************************* c do 4002 jj = 1, ny do 4002 ii = 1, nx if (rland(ii,jj).eq.4.0.and. # b_merge(ii,jj).ge.0.0) then i1 = ii - 1 i2 = ii + 1 js = jj sw = 0.0 sr = 0.0 do 4001 i=i1,i2 w = wgt(i-i1+1) if (i.ge.1.and.i.le.nx) then is = i else if (i.lt.1) then is = i + nx else is = i - nx end if if (b_merge(is,js).ge.0.0) then sw = sw + w sr = sr + w*b_merge(is,js) end if 4001 continue if (sw.gt.0.0) then b_merge(ii,jj)=sr/sw else b_merge(ii,jj)=-999.0 end if end if 4002 continue c do 4004 jj = 1, ny do 4004 ii = 1, nx if (rland(ii,jj).eq.4.0.and. # b_merge(ii,jj).ge.0.0) then is = ii j1 = jj - 1 j2 = jj + 1 if (j1.lt.1) then j1 = 1 end if if (j2.gt.ny) then j2 = ny end if sw=0.0 sr=0.0 do 4003 js = j1, j2 w = wgt(js-j1+1) if (b_merge(is,js).ge.0.0) then sw=sw+w sr=sr+w*b_merge(is,js) end if 4003 continue if (sw.gt.0.0) then b_merge(ii,jj)=sr/sw else b_merge(ii,jj)=-999.0 end if end if 4004 continue c c ************************************************************************* c *** *** c *** 5. to re-define the error in mm/day *** c *** *** c ************************************************************************* c do 5001 jj = 1, ny do 5001 ii = 1, nx if (b_merge(ii,jj).ge.0.0) then c e_merge(ii,jj) = e_merge(ii,jj) * b_merge(ii,jj) e_merge(ii,jj) = e_merge(ii,jj) else e_merge(ii,jj) = -999.0 end if 5001 continue c c ************************************************************************* c *** *** c *** 6. to output *** c *** *** c ************************************************************************* c kout1 = 10 kout2 = 11 write (50,rec=kout1) b_merge write (50,rec=kout2) e_merge c c 9999 return end c c c ************************************************************************* c *** *** c *** *** c *** ***** * * ***** * * ***** ***** * ***** * *** c *** * ** * * * * * * * * * * * * * *** c *** * * * * ***** * * * * * ***** * ***** *** c *** * * ** * * * * * * * * * * * *** c *** ***** * * * ***** * ***** ***** * * * * * *** c *** *** c ************************************************************************* c c ************************************************************************* c *** *** c *** objective : *** c *** to input the individual data sources *** c *** *** c ************************************************************************* c subroutine input_data ( rain, qual ) c include 'CMAP_parameters.h' c dimension rain (nx,ny,nsate+1), # qual (nx,ny,nsate+1) c dimension globe1(nx,ny),globe2(nx,ny) c c 1. to read the rainfall and quality data c do 1002 kk = 1, nsate+1 kinp1 = (kk-1)*2+1 kinp2 = kk*2 read (50,rec=kinp1) globe1 read (50,rec=kinp2) globe2 do 1001 jj=1,ny do 1001 ii=1,nx rain(ii,jj,kk) = globe1(ii,jj) qual(ii,jj,kk) = globe2(ii,jj) 1001 continue 1002 continue c return end c c c ************************************************************************* c *** *** c *** *** c *** ***** ***** ***** ***** * * * ***** ***** ***** *** c *** * * * * * * * * ** * * * * * * *** c *** * * ***** ***** ***** ***** * * * ***** ***** ***** *** c *** * * * * * * * * * ** * * * * * *** c *** ***** ***** * *** * * * * * * *** ***** * * * * *** c *** *** c ************************************************************************* c c ************************************************************************* c *** *** c *** objective : *** c *** to define the random error *** c *** *** c ************************************************************************* c subroutine def_ran_err (nn, yy, xx, ran ) c dimension yy(2000), xx(2000) c c 1. to calculate the random error c nt = 0 sx = 0.0 sy = 0.0 dxy = 0.0 do 1001 kk=1,nn if (yy(kk).ge.0.0.and.xx(kk).ge.0.0) then nt = nt + 1 sx = sx + xx(kk) sy = sy + yy(kk) dxy = dxy + (yy(kk)-xx(kk))*(yy(kk)-xx(kk)) end if 1001 continue if (nt.gt.0) then ax = sx / nt ay = sy / nt bi = ax - ay dxy = dxy - nt*bi*bi dxy = sqrt (dxy/nt) if (ay.gt.0.0) then ran = dxy / ay else ran = 10.0 end if else ran = -999.0 end if c return end c c c ************************************************************************* c *** *** c *** *** c *** ***** ***** ***** ***** ***** ***** * * * * ***** *** c *** * * * * * * * * * * * * ** * * * *** c *** * * ***** ***** ***** ***** ***** * ***** * * * * * *** c *** * * * * * * * * * * * * * ** * * *** c *** ***** ***** * ** ***** * * * * ** ***** * * * * ***** *** c *** *** c ************************************************************************* c c ************************************************************************* c *** *** c *** objective : *** c *** to define error for individual sources *** c *** over land *** c *** *** c ************************************************************************* c subroutine def_err_land ( rland, rain, error ) c include 'CMAP_parameters.h' c dimension rland (nx,ny), # rain (nx,ny,nsate+1), # error (nx,ny,nsate) c dimension yy(2000), xx(2000) c c 1. to define the random error c do 1002 kk = 1, nsate do 1002 jj = 1, ny do 1002 ii = 1, nx if (rland(ii,jj).ge.1.0.and. # rland(ii,jj).le.3.0) then nn = 0 i0 = ii - 3 i1 = ii + 3 j0 = jj - 3 j1 = jj + 3 if (j0.lt.1) then j0 = 1 end if if (j1.gt.ny) then j1 = ny end if do 1001 j = j0, j1 do 1001 i = i0, i1 js = j if (i.ge.1.and.i.le.nx) then is = i else if (i.le.0) then is = i + nx else is = i - nx end if if (rland(is,js).ge.1.0.and. # rland(is,js).le.3.0.and. # rain(is,js,nsate+1).ge.0.0.and. # rain(is,js,kk).ge.0.0) then nn = nn + 1 yy (nn) = rain (is,js,nsate+1) xx (nn) = rain (is,js,kk) end if 1001 continue call def_ran_err (nn, yy, xx, ran) error (ii,jj,kk) = ran else error (ii,jj,kk) = -999.0 end if 1002 continue c return end c c c ************************************************************************* c *** *** c *** *** c *** ***** ***** ***** ***** * * * * ***** ***** * * *** c *** * * * * * * * * * * * * * * * *** c *** * * ***** ***** * * * * * ***** * * * * *** c *** * * * * * ** * * * * * * * * *** c *** ***** ***** * ***** ***** ***** * * ***** ***** * * *** c *** *** c ************************************************************************* c c ************************************************************************* c *** *** c *** objective : *** c *** to define relatively quality *** c *** for individual sources *** c *** (algorithm related) *** c *** *** c *** kseason === 1,2,...12 *** c *** Jan.,Feb.,... Dec. *** c *** *** c ************************************************************************* c subroutine def_quality ( kseason, rland, qua_alg ) c include 'CMAP_parameters.h' c dimension rland(nx,ny), # qua_alg(nx,ny,nsate), # qual(nx,ny) c c ************************************************************************* c *** *** c *** 1. to define the quality coefficients *** c *** for GPI estimates *** c *** *** c ************************************************************************* c call def_quality_gpi (kseason,rland,qual) do 1001 jj=1,ny do 1001 ii=1,nx qua_alg (ii,jj,1) = qual(ii,jj) 1001 continue c c ************************************************************************* c *** *** c *** 2. to define the quality coefficients *** c *** for scattering estimates *** c *** *** c ************************************************************************* c call def_quality_scatter (kseason,rland,qual) do 2001 jj=1,ny do 2001 ii=1,nx qua_alg(ii,jj,2) = qual(ii,jj) 2001 continue c c ************************************************************************* c *** *** c *** 3. to define the weighting coefficients *** c *** for OPI estimates *** c *** *** c ************************************************************************* c call def_quality_opi (kseason,rland,qual) do 3001 jj=1,ny do 3001 ii=1,nx qua_alg(ii,jj,3) = qual (ii,jj) 3001 continue c c ************************************************************************* c *** *** c *** 7. to check the quality coefficients *** c *** *** c ************************************************************************* c return end c c c ************************************************************************* c *** *** c *** *** c *** ***** ***** ***** ***** ***** ***** ***** ***** * * *** c *** * * * * * * * * * * * * ** * *** c *** * * ***** ***** ***** ***** ***** * * * * * * *** c *** * * * * * * * * * * * * * ** *** c *** ***** ***** * ** ***** * * * * ** ***** ***** * * *** c *** *** c ************************************************************************* c c ************************************************************************* c *** *** c *** objective : *** c *** to define error for individual sources *** c *** over ocean *** c *** *** c ************************************************************************* c subroutine def_err_ocn # ( rland, error0, qual0, # qua_sam, qua_alg, error ) c include 'CMAP_parameters.h' c dimension rland (nx,ny), # error0 (nsate), # qual0 (nsate), # qua_sam (nx,ny,nsate+1), # qua_alg (nx,ny,nsate), # error (nx,ny,nsate) c c 1. to define the random error over ocean c do 1001 kk = 1, nsate do 1001 jj = 1, ny do 1001 ii = 1, nx if (rland(ii,jj).eq.0.0.or. # rland(ii,jj).eq.4.0) then error(ii,jj,kk) = -999.0 err0 = error0 (kk) qalg = qua_alg (ii,jj,kk) if (kk.eq.1) then qsam = qua_sam (ii,jj,kk) qua0 = qual0 (kk) else qsam = 1.0 qua0 = 1.0 end if if (err0.ge.0.0.and. # qua0.gt.0.0.and. # qalg.gt.0.0.and. # qsam.gt.0.0) then error(ii,jj,kk) = err0*qua0/(qsam*qalg) end if end if 1001 continue c return end c c c ************************************************************************* c *** *** c *** *** c *** ***** ***** * * ***** ***** ***** ***** *** c *** * * * * * * * * * * *** c *** ***** * * * ***** ***** ***** * ** *** c *** * * * * * * * * * * * *** c *** * ***** ***** ***** **** ***** * * * * *** ***** *** c *** *** c ************************************************************************* c c ************************************************************************* c *** *** c *** objective : *** c *** to fill in the gap in the error fields *** c *** *** c ************************************************************************* c subroutine fill_err_gap ( rain, error ) c include 'CMAP_parameters.h' c dimension rain (nx,ny,nsate+1), # error (nx,ny,nsate) c c 1. to fill in the gap by interpolation c do 1004 kk = 1, nsate 1001 nn = 0 do 1003 jj = 1, ny do 1003 ii = 1, nx if (rain(ii,jj,kk).ge.0.0.and. # error(ii,jj,kk).lt.0.0) then j0 = jj - 3 j1 = jj + 3 i0 = ii - 3 i1 = ii + 3 if (j0.lt.1) then j0 = 1 end if if (j1.gt.ny) then j1 = ny end if nt = 0 se = 0.0 do 1002 j = j0, j1 do 1002 i = i0, i1 js = j if (i.ge.1.and.i.le.nx) then is = i else if (i.lt.1) then is = i + nx else is = i - nx end if if (error(is,js,kk).ge.0.0) then nt = nt + 1 se = se + error(is,js,kk) end if 1002 continue if (nt.ge.1) then error(ii,jj,kk) = se/nt else nn = nn + 1 end if end if 1003 continue if (nn.gt.0) then go to 1001 end if 1004 continue c do 1005 kk = 1, nsate do 1005 jj = 1, ny do 1005 ii = 1, nx if (rain(ii,jj,kk).lt.0.0) then error(ii,jj,kk) = -999.0 end if 1005 continue c return end c c c ************************************************************************* c *** *** c *** ***** ***** * * ***** * * * * * ***** *** c *** * * * ** ** * * * * ** ** * * *** c *** * * * * * * **** * * * * * ***** *** c *** * * * * * * * * * * * * *** c *** ***** ***** * * *** ***** * *** * * ***** ***** *** c *** *** c ************************************************************************* c c ************************************************************************* c *** *** c *** subroutine : *** c *** *** c *** com_by_mle *** c *** *** c *** objective : *** c *** to combine the satellite and model fields *** c *** by the maximum likelihood estimation *** c *** *** c ************************************************************************* c subroutine com_by_mle ( rain, error, comb_rn, comb_er ) c include 'CMAP_parameters.h' c dimension rain (nx,ny,nsate+1), # error (nx,ny,nsate), # comb_rn (nx,ny), # comb_er (nx,ny) c c c 1. to combine c nt = nsate do 1002 jj = 1, ny do 1002 ii = 1, nx sw = 0.0 sr = 0.0 do 1001 ll = 1, nt rn = rain(ii,jj,ll) er = error(ii,jj,ll) if (er.eq.0.0) then er = 0.001 end if if (rn.ge.0.0.and.er.gt.0.0) then w = 1.0 / (er*er) sw = sw + w sr = sr + w * rn end if 1001 continue if (sw.gt.0.0) then comb_rn(ii,jj) = sr/sw comb_er(ii,jj) = 1.0 / sqrt(sw) else comb_rn(ii,jj) = -999.0 comb_er(ii,jj) = -999.0 end if 1002 continue c return end c c c ************************************************************************* c *** *** c *** *** c *** ***** * * * * ***** ***** * ***** *** c *** * * ** ** * * * * * * * * *** c *** ***** * * * * * **** * ***** ***** *** c *** * * * * * * * * * * * * *** c *** * * * * * ***** ***** ***** * * ***** *** c *** *** c *** * * * * ***** *** c *** * * * ** * * * *** c *** * ***** * * * * * *** c *** * * * * ** * * *** c *** ***** ***** * * * * ***** *** c *** *** c *** *** c ************************************************************************* c c ************************************************************************* c *** *** c *** objective : *** c *** to remove the bias from the combined field *** c *** for land areas by blending *** c *** *** c ************************************************************************* c subroutine rmv_bias_land (n_gage_lim,rland,kcheck,gage, # gagen,satcom,blend,kbench) c include 'CMAP_parameters.h' c dimension rland(nx,ny),gage(nx,ny),gagen(nx,ny), # satcom(nx,ny),blend(nx,ny) c dimension kcheck(nx,ny),kbench(nx,ny) c c c ************************************************************************* c *** *** c *** 1. to check the data availability *** c *** *** c ************************************************************************* c ndata0 = 0 ndata1 = 0 ndata2 = 0 do 1001 jj = 1, ny do 1001 ii = 1, nx if (rland(ii,jj).ge.1.0.and. # rland(ii,jj).le.3.0) then ndata0 = ndata0 + 1 if (gage (ii,jj).ge.0.0) then ndata1 = ndata1 + 1 end if if (satcom(ii,jj).ge.0.0) then ndata2 = ndata2 + 1 end if end if 1001 continue ndata0 = ndata0 * 0.3 if (ndata1.le.ndata0.or. # ndata1.le.ndata0) then do 1002 jj = 1, ny do 1002 ii = 1, nx blend (ii,jj) = -999.0 1002 continue go to 9999 end if c c ************************************************************************* c *** *** c *** 1. to call the blending sub-routine version 9403 *** c *** *** c ************************************************************************* c call blending_v9410 (n_gage_lim,rland,kcheck,gage, # gagen,satcom,blend,kbench) do 2001 jj = 1, ny do 2001 ii = 1, nx if (rland(ii,jj).ge.1.0.and. # rland(ii,jj).le.3.0) then if (blend(ii,jj).lt.0.0.and. # blend(ii,jj).ne.-999.0) then blend(ii,jj)=0.0 end if end if 2001 continue c c ************************************************************************* c *** *** c *** 2. to check the results *** c *** *** c ************************************************************************* c 9999 return end c c c ************************************************************************* c *** *** c *** *** c *** ***** * * * * ***** ***** * ***** *** c *** * * ** ** * * * * * * * * *** c *** ***** * * * * * **** * ***** ***** *** c *** * * * * * * * * * * * * *** c *** * * * * * ***** ***** ***** * * ***** *** c *** *** c *** ***** ***** ***** * * * *** c *** * * * * * * ** * *** c *** * * * ***** ***** * * * *** c *** * * * * * * * ** *** c *** ***** ***** ***** ***** * * * * *** c *** *** c *** *** c ************************************************************************* c c ************************************************************************* c *** *** c *** objective : *** c *** to remove the bias from the combined field *** c *** for oceanic areas by comparing with atoll *** c *** *** c ************************************************************************* c subroutine rmv_bias_ocean # ( rland, bias_atl, c_merge, b_merge ) c include 'CMAP_parameters.h' c dimension rland(nx,ny), # bias(nx,ny), # c_merge(nx,ny), # b_merge(nx,ny) c c ************************************************************************* c *** *** c *** 1. to define the bias structure over ocean *** c *** *** c ************************************************************************* c tar=bias_atl do 1002 jj=1,ny if (jj.le.20) then wlat=0.0 else if (jj.le.28) then wlat=float(jj-20)/8.0 else if (jj.le.44) then wlat=1.0 else if (jj.le.52) then wlat=float(52-jj)/8.0 else wlat=0.0 end if do 1001 ii=1,nx if (rland(ii,jj).eq.0.0.or. # rland(ii,jj).eq.4.0) then bias(ii,jj)=wlat*tar end if 1001 continue 1002 continue c c ************************************************************************* c *** *** c *** 2. to remove the bias from the combined satellite/model *** c *** *** c ************************************************************************* c do 2001 jj=1,ny do 2001 ii=1,nx if (rland(ii,jj).eq.0.0.or. # rland(ii,jj).eq.4.0) then if (c_merge(ii,jj).ne.-999.0.and. # bias(ii,jj).ne.-999.0) then c b_merge(ii,jj)=c_merge(ii,jj)/(1.0+bias(ii,jj)) b_merge(ii,jj)=c_merge(ii,jj)*(1.0-bias(ii,jj)) else b_merge(ii,jj)=-999.0 end if if (b_merge(ii,jj).lt.0.0.and. # b_merge(ii,jj).ne.-999.0) then b_merge(ii,jj)=0.0 end if end if 2001 continue c c ************************************************************************* c *** *** c *** 4. to check the results *** c *** *** c ************************************************************************* c return end c c c ************************************************************************* c *** *** c *** *** c *** ***** ***** ***** ***** * * * * ***** ***** * * *** c *** * * * * * * * * * * * * * * * *** c *** * * ***** ***** * * * * * ***** * * * * *** c *** * * * * * ** * * * * * * * * *** c *** ***** ***** * ***** ***** ***** * * ***** ***** * * *** c *** *** c *** ***** ***** ***** *** c *** * * * * *** c *** * *** ***** * *** c *** * * * * * *** c *** ***** ***** * ***** *** c *** *** c ************************************************************************* c c ************************************************************************* c *** *** c *** subroutine : *** c *** *** c *** def_quality_gpi *** c *** *** c *** objective : *** c *** *** c *** To define the quality coefficients *** c *** for IR-based GPI estimates *** c *** *** c ************************************************************************* subroutine def_quality_gpi (kseason,rland,quality) c include 'CMAP_parameters.h' c dimension rland(nx,ny),quality(nx,ny) c dimension bsouth1(12),bsouth2(12),bnorth1(12),bnorth2(12) data bsouth1/-40.0,-40.0,-40.0,-40.0,-40.0,-40.0, # -40.0,-40.0,-40.0,-40.0,-40.0,-40.0/, # bsouth2/-25.0,-25.0,-25.0,-21.7,-18.3,-15.0, # -15.0,-15.0,-15.0,-18.3,-21.7,-25.0/, # bnorth1/ 15.0, 15.0, 15.0, 18.3, 21.7, 25.0, # 25.0, 25.0, 25.0, 21.7, 18.3, 15.0/, # bnorth2/ 40.0, 40.0, 40.0, 40.0, 40.0, 40.0, # 40.0, 40.0, 40.0, 40.0, 40.0, 40.0/ c c ************************************************************************* c *** *** c *** 1. to define the boundaries *** c *** *** c ************************************************************************* c bs1=bsouth1(kseason) bs2=bsouth2(kseason) bn1=bnorth1(kseason) bn2=bnorth2(kseason) c c ************************************************************************* c *** *** c *** 2. to define the weighting coefficients *** c *** *** c ************************************************************************* c do 2002 jj=1,72 rlat=(jj-1)*2.5-88.75 if (rlat.le.bs1) then wgt_lat=0.0 else if (rlat.le.bs2) then wgt_lat=(rlat-bs1)/(bs2-bs1) else if (rlat.le.bn1) then wgt_lat=1.0 else if (rlat.le.bn2) then wgt_lat=(bn2-rlat)/(bn2-bn1) else wgt_lat=0.0 end if do 2001 ii=1,144 quality(ii,jj)=wgt_lat if (rland(ii,jj).ge.1.0.and. # rland(ii,jj).le.3.0) then quality(ii,jj)=0.0 end if 2001 continue 2002 continue c return end c c c ************************************************************************* c *** *** c *** *** c *** ***** ***** ***** ***** * * * * ***** ***** * * *** c *** * * * * * * * * * * * * * * * *** c *** * * ***** ***** * * * * * ***** * * * * *** c *** * * * * * ** * * * * * * * * *** c *** ***** ***** * ***** ***** ***** * * ***** ***** * * *** c *** *** c *** *** c *** ***** ***** * ***** ***** ***** ***** *** c *** * * * * * * * * * *** c *** ***** * ***** * * ***** ***** *** c *** * * * * * * * * * *** c *** ***** ***** ***** * * * * ***** * * *** c *** *** c ************************************************************************* c c ************************************************************************* c *** *** c *** subroutine : *** c *** *** c *** def_quality_scatter *** c *** *** c *** objective : *** c *** *** c *** To define the quality coefficients *** c *** for MW scattering estimates *** c *** *** c ************************************************************************* c subroutine def_quality_scatter (kseason,rland,quality) c include 'CMAP_parameters.h' c dimension rland(nx,ny),quality(nx,ny) c dimension bsouth1(12),bsouth2(12),bnorth1(12),bnorth2(12), # vs1_oce(12),vs2_oce(12),vn1_oce(12),vn2_oce(12) c data bsouth1/-70.0,-70.0,-70.0,-70.0,-70.0,-70.0, # -70.0,-70.0,-70.0,-70.0,-70.0,-70.0/, # bsouth2/-25.0,-25.0,-25.0,-21.7,-18.3,-15.0, # -15.0,-15.0,-15.0,-18.3,-21.7,-25.0/, # bnorth1/ 15.0, 15.0, 15.0, 18.3, 21.7, 25.0, # 25.0, 25.0, 25.0, 21.7, 18.3, 15.0/, # bnorth2/ 70.0, 70.0, 70.0, 70.0, 70.0, 70.0, # 70.0, 70.0, 70.0, 70.0, 70.0, 70.0/, # vs1_oce/ 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, # 0.2, 0.2, 0.2, 0.2, 0.2, 0.2/, # vs2_oce/ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, # 1.0, 1.0, 1.0, 1.0, 1.0, 1.0/, # vn1_oce/ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, # 1.0, 1.0, 1.0, 1.0, 1.0, 1.0/, # vn2_oce/ 0.5, 0.5, 0.55, 0.60, 0.65, 0.7, # 0.7, 0.7, 0.7, 0.63, 0.57, 0.5/ c c ************************************************************************* c *** *** c *** 1. to define the boundaries *** c *** *** c ************************************************************************* c bs1=bsouth1(kseason) bs2=bsouth2(kseason) bn1=bnorth1(kseason) bn2=bnorth2(kseason) c c ************************************************************************* c *** *** c *** 2. to define the vaules at the boundaries *** c *** *** c ************************************************************************* c vs1=vs1_oce(kseason) vs2=vs2_oce(kseason) vn1=vn1_oce(kseason) vn2=vn2_oce(kseason) c c c ************************************************************************* c *** *** c *** 3. to define the quality coefficients *** c *** *** c ************************************************************************* c do 3001 jj=1,72 do 3001 ii=1,144 if (rland(ii,jj).eq.0.0.or. # rland(ii,jj).eq.4.0) then rlat=(jj-1)*2.5-88.75 if (rlat.le.bs1) then wgt_lat=vs1 else if (rlat.le.bs2) then wgt_lat=(rlat-bs1)/(bs2-bs1)*(vs2-vs1)+vs1 else if (rlat.le.bn1) then wgt_lat=vs2 else if (rlat.lt.bn2) then wgt_lat=(bn2-rlat)/(bn2-bn1)*(vn1-vn2)+vn2 else wgt_lat=vn2 end if quality(ii,jj)=wgt_lat else quality(ii,jj)=0.0 end if 3001 continue c return end c c c ************************************************************************* c *** *** c *** *** c *** ***** ***** ***** ***** * * * * ***** ***** * * *** c *** * * * * * * * * * * * * * * * *** c *** * * ***** ***** * * * * * ***** * * * * *** c *** * * * * * ** * * * * * * * * *** c *** ***** ***** * ***** ***** ***** * * ***** ***** * * *** c *** *** c *** ***** ***** ***** *** c *** * * * * * *** c *** * * ***** * *** c *** * * * * *** c *** ***** ***** * ***** *** c *** *** c ************************************************************************* c c ************************************************************************* c *** *** c *** subroutine : *** c *** *** c *** def_quality_opi *** c *** *** c *** objective : *** c *** *** c *** To define the quality coefficients *** c *** for OPI estimates *** c *** *** c ************************************************************************* c subroutine def_quality_opi (kseason,rland,quality) c include 'CMAP_parameters.h' c dimension rland(nx,ny),quality(nx,ny),weight(nx,ny) c c c ************************************************************************* c *** *** c *** 1. to define the quality coefficients *** c *** *** c ************************************************************************* c do 2001 jj=1,ny do 2001 ii=1,nx quality(ii,jj) = 1.0 2001 continue c return end c c c ************************************************************************* c *** *** c *** *** c *** ******* * ******* * * ******* *** c *** * * * * ** * * * *** c *** * * * * * * * * * *** c *** ***** * ******* * * * * * *** c *** * * * * * * * * * *** c *** * * * * * ** * * *** c *** ******* ******* ******* * * ******* *** c *** *** c *** *** c *** *** c *** * * ***** ***** ***** ***** ***** ***** *** c *** * * * * * * * * * * * *** c *** * * ***** ***** ***** ***** * * ***** *** c *** * * * * * ** * * * * * *** c *** * ***** * * ** ***** ***** ***** ***** *** c *** *** c *** *** c ************************************************************************* c c ************************************************************************* c *** *** c *** Program : *** c *** blending.f *** c *** *** c *** Objective : *** c *** *** c *** to blend the satellite estimates and gage observation *** c *** *** c *** of monthly rainfall by applying Laplacian interplotation *** c *** *** c *** Programmer : *** c *** *** c *** Pingping Xie *** c *** National Meteorological Center *** c *** Rm. 800, World Weather Building *** c *** 5200 Auth Road *** c *** Camp Springs, MD 20746 *** c *** *** c *** Tel. 301-763-8167 *** c *** Fax. 301-763-8125 *** c *** E-mail xping@sgi17.wwb.noaa.gov *** c *** *** c ************************************************************************* c subroutine blending_v9410 (n_gage_lim,rland,kcheck, # gage,gagen,satcom,blend,kbench) c dimension rland(144,72),kcheck(144,72),gage(144,72), # gagen(144,72),satcom(144,72),blend(144,72) c dimension first(144,72),kbench(144,72) c c ************************************************************************* c *** *** c *** 1. to define the blending parameters *** c *** *** c ************************************************************************* c call def_para (flp_lim,alpha,mod_smooth) c c ************************************************************************* c *** *** c *** 2. to creat the first guess of blend field *** c *** *** c ************************************************************************* c call def_first (rland,gage,gagen,satcom,first) c c ************************************************************************* c *** *** c *** 3. to define the bench mark grids *** c *** *** c ************************************************************************* c call def_bench (n_gage_lim,rland,kcheck,gagen,kbench) c c ************************************************************************* c *** *** c *** 4. to blend the satellite and gage data *** c *** *** c ************************************************************************* c call cal_blend (flp_lim,alpha,mod_smooth, # first,satcom,kbench,blend) c 9999 return end c c c ************************************************************************* c *** *** c *** *** c *** ***** ***** ***** ***** * ***** * *** c *** * * * * * * * * * * * * *** c *** * * ***** ***** ***** ***** ***** ***** *** c *** * * * * * * * * * * * *** c *** ***** ***** * ***** * * * * * * * *** c *** *** c *** *** c ************************************************************************* c c ************************************************************************* c *** *** c *** Subroutine : *** c *** def_para *** c *** *** c *** Objective : *** c *** to define the operation modes and calculation *** c *** *** c *** parameters to be used in the followed *** c *** *** c *** subroutines *** c *** *** c ************************************************************************* c subroutine def_para (flp_lim,alpha,mod_smooth) c c 1. first guess definition c c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c +++ +++ c +++ n_gage_lim == limitation of gage number. The gage +++ c +++ observation at grids with at least n_gage_lim +++ c +++ gage will be used as bench mark +++ c +++ +++ c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c n_gage_lim = 1 c c 2. blend procedures c c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c +++ +++ c +++ flp_lim == limitation of difference in Laplacian +++ c +++ between two consecssive calculation. +++ c +++ Calculation will be regarded as completed +++ c +++ when the differences for all grids are smaller +++ c +++ than the limitation. +++ c +++ +++ c +++ alpha == coeffifient used to revise the field values +++ c +++ ( 0.0 -- 0.5 ) +++ c +++ +++ c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c flp_lim = 0.1 alpha = 0.25 c c c 3. smooth process c c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c +++ +++ c +++ mod_smooth == smooth process mode +++ c +++ +++ c +++ 0 --- no smooth process +++ c +++ 1 --- bi-directional linear smoothing +++ c +++ 2 --- smoothing but keep benchmark unchange +++ c +++ 3 --- reversed for future extension +++ c +++ +++ c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c mod_smooth = 2 c return end c c c ************************************************************************* c *** *** c *** *** c *** ***** ***** ***** ***** ***** ***** ***** ***** *** c *** * * * * * * * * * * *** c *** * * ***** ***** ***** * ***** ***** * *** c *** * * * * * * * * * * *** c *** ***** ***** * ***** * ***** * * ***** * *** c *** *** c *** *** c ************************************************************************* c c ************************************************************************* c *** *** c *** Subroutine : *** c *** *** c *** def_first *** c *** *** c *** Objective : *** c *** *** c *** To define the first guess of the blend field *** c *** *** c ************************************************************************* c subroutine def_first (rland,gage,gagen,satcom,first) c dimension rland(144,72),gage(144,72),gagen(144,72), # satcom(144,72),first(144,72) c dimension egage(164,92),egagen(164,92),wgt(6) data wgt/0.1,1.0,2.0,3.0,5.0,10.0/ c c ************************************************************************* c *** *** c *** 1. to fill in the undefined gauge field resulted from *** c *** different land/ocean masks used by GPCC and us *** c *** *** c ************************************************************************* c call def_edge (gagen,egagen) c krad=0 1001 krad=krad+1 dis0=(krad+1.0)*(krad+1.0) call def_edge (gage,egage) do 1003 jj=1,72 do 1003 ii=1,144 if (rland(ii,jj).ge.1.0.and. # rland(ii,jj).le.3.0.and. # gage(ii,jj).lt.0.0) then i0=ii+10 j0=jj+10 i1=i0-krad i2=i0+krad j1=j0-krad j2=j0+krad sw=0.0 sr=0.0 do 1002 j=j1,j2 do 1002 i=i1,i2 if (i.ge.1.and. # i.le.164.and. # j.ge.1.and. # j.le.92.and. # egage(i,j).ge.0.0) then kk=nint(egagen(i,j))+1 if (kk.gt.6) then kk=6 end if w_gage=wgt(kk) dis=(j-j0)*(j-j0)+(i-i0)*(i-i0) w_dis=(dis0-dis)/(dis0+dis) if (w_dis.lt.0.0) then w_dis=0.0 end if w=w_dis*w_gage sw=sw+w sr=sr+w*egage(i,j) end if 1002 continue if (sw.gt.0.0) then gage(ii,jj)=sr/sw else go to 1001 end if end if 1003 continue c c ************************************************************************* c *** *** c *** 2. to define the first guess *** c *** == gauge analysis over land *** c *** == satellite over ocean *** c *** *** c ************************************************************************* c do 2001 jj=1,72 do 2001 ii=1,144 if (rland(ii,jj).ge.1.0.and. # rland(ii,jj).le.3.0) then first(ii,jj)=gage(ii,jj) else first(ii,jj)=satcom(ii,jj) end if 2001 continue c return end c c c ************************************************************************* c *** *** c *** *** c *** ***** ***** ***** ***** ***** * * ***** * * *** c *** * * * * * * * ** * * * * *** c *** * * ***** ***** **** ***** * * * * ***** *** c *** * * * * * * * * ** * * * *** c *** ***** ***** * ***** ***** ***** * * ***** * * *** c *** *** c *** *** c ************************************************************************* c c ************************************************************************* c *** *** c *** Subroutine : *** c *** def_bench *** c *** *** c *** Objective : *** c *** to define the bench mark grids *** c *** *** c ************************************************************************* c subroutine def_bench (n_gage_lim,rland,kcheck,gagen,kbench) c dimension rland(144,72),kcheck(144,72),gagen(144,72), # kbench(144,72) c c ************************************************************************* c *** *** c *** 1. to define the bench mark grids *** c *** *** c ************************************************************************* c do 1001 jj=1,72 do 1001 ii=1,144 ld=nint(rland(ii,jj)) if (ld.eq.0.or.ld.eq.4) then kbench(ii,jj)=1 else if (ld.ge.1.and.ld.le.2) then ng=nint(gagen(ii,jj)) if (ng.ge.n_gage_lim) then kbench(ii,jj)=1 else kbench(ii,jj)=0 end if if (kcheck(ii,jj).eq.0) then kbench(ii,jj)=1 end if else if (ld.eq.3) then kbench(ii,jj)=1 else kbench(ii,jj)=0 end if 1001 continue c return end c c c ************************************************************************* c *** *** c *** *** c *** ***** * * ***** * ***** * * ***** *** c *** * * * * * * * * ** * * * *** c *** * ***** * **** * ***** * * * * * *** c *** * * * * * * * * * ** * * *** c *** ***** * * ***** ***** ***** ***** ***** * * ***** *** c *** *** c *** *** c ************************************************************************* c c ************************************************************************* c *** *** c *** Subroutine : *** c *** *** c *** cal_blend *** c *** *** c *** Objective : *** c *** *** c *** To construct the blend analysis of gage *** c *** *** c *** observation and satellite estimate of *** c *** *** c *** monthly rainfall using Laplacian *** c *** *** c *** interplotation *** c *** *** c ************************************************************************* c subroutine cal_blend (flp_lim,alpha,mod_smooth, # first,satcom,kbench,blend) c dimension first(144,72),satcom(144,72),kbench(144,72), # blend(144,72) c dimension field(164,92),flap(144,72), # rain(144,72),erain(164,92), # kchk(144,72) c c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c +++ +++ c +++ flp_lim == limitation of difference in Laplacian +++ c +++ between two consecssive calculation. +++ c +++ Calculation will be regarded as completed +++ c +++ when the differences for all grids are smaller +++ c +++ than the limitation. +++ c +++ +++ c +++ alpha == coeffifient used to revise the field values +++ c +++ ( 0.0 -- 0.5 ) +++ c +++ +++ c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c c ************************************************************************* c *** *** c *** 1. to define the satellite estimate to define the field shape *** c *** *** c ************************************************************************* c 1001 call def_edge (satcom,field) c c ************************************************************************* c *** *** c *** 2. to define the field of Laplacian of satellite estimate *** c *** *** c ************************************************************************* c do 2001 jj=1,72 do 2001 ii=1,144 iic=ii+10 jjc=jj+10 fc=field(iic,jjc) fn=field(iic,jjc+1) fs=field(iic,jjc-1) fw=field(iic-1,jjc) fe=field(iic+1,jjc) if (fc.ne.-999.0.and.fn.ne.-999.0.and. # fs.ne.-999.0.and. # fw.ne.-999.0.and.fe.ne.-999.0) then flap(ii,jj)=fn+fs+fw+fe-4.0*fc else flap(ii,jj)=0.0 end if 2001 continue c c ************************************************************************* c *** *** c *** 3. to define the blend field *** c *** *** c ************************************************************************* c do 3001 jj=1,72 do 3001 ii=1,144 rain(ii,jj)=first(ii,jj) 3001 continue c 3002 n_undef=0 call def_edge (rain,erain) do 3003 jj=1,72 do 3003 ii=1,144 if (kbench(ii,jj).eq.0) then rc=erain(ii+10,jj+10) rn=erain(ii+10,jj+11) rs=erain(ii+10,jj+9) rw=erain(ii+9,jj+10) re=erain(ii+11,jj+10) rlap=rn+rs+rw+re-4.0*rc rdif=rlap-flap(ii,jj) if (abs(rdif).gt.flp_lim) then n_undef=n_undef+1 rain(ii,jj)=rain(ii,jj)+alpha*rdif end if end if 3003 continue if (n_undef.ne.0) then go to 3002 else do 3004 jj=1,72 do 3004 ii=1,144 if (rain(ii,jj).ge.0.0) then blend(ii,jj)=rain(ii,jj) else blend(ii,jj)=0.0 rain(ii,jj)=0.0 end if 3004 continue end if c c ************************************************************************* c *** *** c *** 4. to smooth the results *** c *** *** c ************************************************************************* c if (mod_blend.eq.2) then call smoothing (mod_smooth,rain,kbench,blend) end if c 9999 return end c c c ************************************************************************* c *** *** c *** *** c *** ***** * * ***** ***** ***** * * ***** * * ***** *** c *** * ** ** * * * * * * * * ** * * *** c *** ***** * * * * * * * * ***** * * * * * ** *** c *** * * * * * * * * * * * * ** * * *** c *** ***** * * ***** ***** * * * ***** * * ***** *** c *** *** c *** *** c ************************************************************************* c c ************************************************************************* c *** *** c *** Subroutine : *** c *** *** c *** smoothing *** c *** *** c *** Objective : *** c *** *** c *** To remove the unreal and unwanted small scale *** c *** *** c *** features in the original blended analysis *** c *** *** c *** by using a linear ninomial (1-2-1) smoothing *** c *** *** c *** filter *** c *** *** c ************************************************************************* c subroutine smoothing (mod_smooth,origin,kcheck,smooth) c dimension origin(144,72),kcheck(144,72),smooth(144,72), # rain(164,92),wei(3) data wei/1.0,2.0,1.0/ c c c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c +++ +++ c +++ mod_smooth == smooth process mode +++ c +++ +++ c +++ 0 --- no smooth process +++ c +++ 1 --- bi-directional linear smoothing +++ c +++ 2 --- smoothing but keep benchmark unchange +++ c +++ 3 --- reversed for future extension +++ c +++ +++ c +++ kcheck == bench mark check +++ c +++ 0 --- no bench mark +++ c +++ 1 --- bench mark +++ c +++ +++ c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c c 1. check the operation mode c if (mod_smooth.eq.1) then go to 2001 else if (mod_smooth.eq.2) then go to 2001 else write (6,*) 'Wrong smoothing mode !!!' Write (6,*) 'Please choose 1 or 2 !!!' stop end if c c 2. linear bi-nomial smoothing c 2001 call def_edge (origin,rain) do 2003 jj=11,82 do 2003 ii=11,154 if (mod_smooth.ne.2.or. # kcheck(ii-10,jj-10).ne.1) then sw=0.0 sr=0.0 do 2002 i=1,3 if (rain(ii-2+i,jj).ne.-999.0) then sw=sw+wei(i) sr=sr+wei(i)*rain(ii-2+i,jj) end if 2002 continue if (sw.ne.0.0) then rain(ii,jj)=sr/sw else rain(ii,jj)=-999.0 end if end if 2003 continue do 2004 jj=1,72 do 2004 ii=1,144 smooth(ii,jj)=rain(ii+10,jj+10) 2004 continue call def_edge (smooth,rain) do 2006 jj=11,82 do 2006 ii=11,154 if (mod_smooth.ne.2.or. # kcheck(ii-10,jj-10).ne.1) then sw=0.0 sr=0.0 do 2005 j=1,3 if (rain(ii,jj-2+j).ne.-999.0) then sw=sw+wei(j) sr=sr+wei(j)*rain(ii,jj-2+j) end if 2005 continue if (sw.ne.0.0) then rain(ii,jj)=sr/sw else rain(ii,jj)=-999.0 end if end if 2006 continue c do 2007 jj=1,72 do 2007 ii=1,144 smooth(ii,jj)=rain(ii+10,jj+10) 2007 continue c 9999 return end c c c c ************************************************************************* c *** *** c *** *** c *** ***** ***** ***** ***** ***** ***** ***** *** c *** * * * * * * * * * *** c *** * * ***** ***** ***** * * * ** ***** *** c *** * * * * * * * * * * *** c *** ***** ***** * ***** ***** ***** ***** ***** *** c *** *** c *** *** c ************************************************************************* c c ************************************************************************* c *** *** c *** subroutine : *** c *** *** c *** def_edge *** c *** *** c *** objective : *** c *** *** c *** To define the extra-boundary vaules *** c *** *** c *** by using mirror effect *** c *** *** c ************************************************************************* c subroutine def_edge (rinp,rout) c dimension rinp(144,72),rout(164,92) c c ************************************************************************* c *** *** c *** 1. to define the central portion *** c *** *** c ************************************************************************* c do 1001 jj=11,82 do 1001 ii=11,154 rout(ii,jj)=rinp(ii-10,jj-10) 1001 continue c c ************************************************************************* c *** *** c *** 2. to define the mirror edge *** c *** *** c ************************************************************************* c do 2003 jj=11,82 do 2001 ii=1,10 rout(ii,jj)=rout(ii+144,jj) 2001 continue do 2002 ii=155,164 rout(ii,jj)=rout(ii-144,jj) 2002 continue 2003 continue c do 2006 ii=1,164 if (ii.le.92) then iim=ii+72 else iim=ii-72 end if iim=ii do 2004 jj=1,10 rout(ii,jj)=rout(iim,22-jj) 2004 continue do 2005 jj=83,92 rout(ii,jj)=rout(iim,164-jj) 2005 continue 2006 continue c return end