***************************************************************** * PROGRAM: globe_opi_pen_def_opi_pen_v0709.lnx.f * * LANGUAGE: Fortran77 * * MACHINE: Any LINUX BOX *----------------------------------------------------------------* * PURPOSE: to define pentad OPI on 2.5deg lat/lon *----------------------------------------------------------------* * USAGE: *----------------------------------------------------------------* * INPUT FILES: Unit 10 == availability of pentad OLR data * 20 == yearly data file of pentad OLR * OUTPUT FILES: Unit 50 == GrADS formatted pentad OPI on 2.5deg lat/lon * *----------------------------------------------------------------* * SUBROUTINES USED: NONE * * 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 * * LOCAL VARIABLES: clim_olr == pentad OLR climatology * clim_pre == pentad precipitation climatology * aa1 == coefficients to convert OLR into precip * bb1 == coefficients to convert OLR into precip * olr == OLR on 144x72 grid system * olr0 == OLR on 144x73 grid system * opi0 == OPI estimates based on anomaly alone * olr1 == OPI estimates based on total OLR * olr == OPI estimates based on OPI0 and OPI1 (final product) * *----------------------------------------------------------------* * AUTHOR(S): Pingping Xie *----------------------------------------------------------------* * DATE 2007-09-12 *----------------------------------------------------------------* * MODIFICATIONS: ****************************************************************** c=====|==|=========|=========|=========|=========|=========|=========|== c program OPI_def_pen_2_5deg implicit none integer NX, NY, NY1, NC parameter( NX=144, NY=72, NY1=73, NC=NX*NY1*4) character yyyymmdd*8, f_olr*23, f_clm*38, f_pen*46 character cin*(NC), cout*(NC) integer yyyy, mm, dd, days(12), leap, m, jjj, ddd, pen, mon real olr0(NX,NY1) real clim_olr(NX,NY), clim_pre(NX,NY), aa1(NX,NY), bb1(NX,NY) real opi0(NX,NY), opi1(NX,NY), opi(NX,NY) real aa0(2), bb0(2), a0, b0 integer ii, jj, mdl, kk, k1, k2, k3, k4 data aa0 / -0.02774, -0.02794 / data bb0 / -0.02943, -0.03760 / 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 get Julian day of target date c----------------------------------------------------------------------- jjj=dd do m=1, mm-1 jjj=jjj+days(m) enddo c----------------------------------------------------------------------- c determine if target date is end of a 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 determine month c----------------------------------------------------------------------- mon=(pen-1)/6+1 if( mon.gt.12 ) mon=12 c----------------------------------------------------------------------- c construct file names c----------------------------------------------------------------------- c----------->123456789|123456789|123456789|123456789|123456<-- f_olr='../input/pent.olr.y????' f_clm='../library/cmap_opi_pen_coef_v0003.lnx' f_pen='../output/OPI/yyyy/cmap_opi_pen_v0003.lnx.yyyy' write(f_olr(20:23), '(i4.4)') yyyy write(f_pen(15:18), '(i4.4)') yyyy write(f_pen(43:46), '(i4.4)') yyyy open(10, file=f_olr, status='old', 1 access='direct', recl=NX*NY1*4) open(11, file=f_clm, status='old', 1 access='direct', recl=NX*NY*4) write(*, *) ' <-- ', f_olr open(20, file=f_pen, access='direct', recl=NX*NY*4) c----------------------------------------------------------------------- c repeat for two model options c----------------------------------------------------------------------- do mdl=1, 2 c----------------------------------------------------------------------- c initialize arrays c----------------------------------------------------------------------- do jj=1, NY do ii=1, NX opi0(ii,jj)=-999.0 opi1(ii,jj)=-999.0 opi(ii,jj) =-999.0 enddo enddo c----------------------------------------------------------------------- c swap endian and write back (only do this once) c----------------------------------------------------------------------- if( mdl.eq.1 ) then read(10, rec=pen) cin do jj=1, NY1 do ii=1, NX kk=(jj-1)*NX+ii k1=(kk-1)*4+1 k2=(kk-1)*4+2 k3=(kk-1)*4+3 k4=(kk-1)*4+4 cout(k1:k1)=cin(k4:k4) cout(k2:k2)=cin(k3:k3) cout(k3:k3)=cin(k2:k2) cout(k4:k4)=cin(k1:k1) enddo enddo write(10, rec=pen) cout endif c----------------------------------------------------------------------- c read raw OLR 144*73 data and coversion coef c----------------------------------------------------------------------- read(10, rec=pen) olr0 read(11, rec=(pen-1)*8+mdl+0) clim_olr read(11, rec=(pen-1)*8+mdl+2) clim_pre read(11, rec=(pen-1)*8+mdl+4) aa1 read(11, rec=(pen-1)*8+mdl+6) bb1 c----------------------------------------------------------------------- c define OPI based on all sources c----------------------------------------------------------------------- a0=aa0(mdl) b0=bb0(mdl) call get_opi(mon, olr0, clim_olr, clim_pre, a0, b0, aa1, bb1, 1 opi0, opi1, opi) c----------------------------------------------------------------------- c output c----------------------------------------------------------------------- write(20, rec=(pen-1)*6+mdl+0) opi0 write(20, rec=(pen-1)*6+mdl+2) opi1 write(20, rec=(pen-1)*6+mdl+4) opi write(*, '(i4.4,2i2.2,1x,i3.3,3(a6,i3))') yyyy, mm, dd, jjj, 1 ' pen: ', pen, ' mon: ', mon, ' mdl: ', mdl enddo close(10) close(11) close(20) write(*, *) ' ==> ', f_pen stop 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 * subroutine: * c * * c * get_opi * c * * c * objective : * c * * c * to define global monthly precipitation estimates * c * from OLR data observed from NOAA polar satellites * c * * c * input : * c * * c * kmon : calendar month of the data * c * 1 - 12 ==> (Jan. .. Dec.) * c * * c * olr0 : NOAA OLR data (W/m^2) at 144*73 grid * c * olr0 (i,j) * c * i = 1, 144 ==> 0.0, 2.5E, ... 2.5W * c * j = 1, 73 ==> 90.0N, 87.5N , .. 90.0S * c * * c * clim_olr: monthly climatology (W/m^2) of OLR * c * clim_olr (i,j) * c * i = 1, 144 ==> 1.25E, 3.75E, ... 1.25W * c * j = 1, 72 ==> 88.75S, 86.25N, ... 88.75N* c * * c * clim_pre: monthly climatology (mm/day) of precip * c * clim_pre (i,j) * c * i = 1, 144 ==> 1.25E, 3.75E, ... 1.25W * c * j = 1, 72 ==> 88.75S, 86.25N, ... 88.75N* c * * c * aa0/bb0 : coefficients used to define OPI0 * c * PC = aa0 * clim + bb0 * c * precip anom = PC * OLR anom * c * precip totl = precip clim + precip anom * c * * c * aa1/bb1 : coefficients used to define OPI1 * c * precip total = aa1 * OLR totl + bb1 * c * aa1 (144,72), bb1 (144,72) * c * * c * * c * * c * output : * c * * c * opi0 : precipitation estimates (mm/day) * c * based on the original OPI algorithm * c * (Xie and Arkin, 1998) * c * * c * opi1 : precipitation estimates (mm/day) * c * based on regression between total precip * c * and total OLR * c * * c * opi : combination of opi0 and opi1 * c * (better than either opi0 or opi1) * c * * c * opi0 (i,j), opi1 (i,j), opi (i,j) * c * i = 1, 144 ==> 1.25E, 3.75E, ... 1.25W * c * j = 1, 72 ==> 88.75S, 86.26S, .. 88.75N * c * * c * programmer: * c * * c * Pingping Xie * c * NOAA Climate Prediction Center * c * 5200 Auth Road, #805B * c * Camp Springs, MD 20746 * c * * c * Tel. (301)763-8000, ext.7572 * c * Fax (301)763-8125 * c * E-mail xping@sgi17.wwb.noaa.gov * c * * c * * c ************************************************************************ c subroutine get_opi ( kmon, olr0, clim_olr, clim_pre, # aa0 , bb0, aa1, bb1, # opi0, opi1, opi ) c dimension olr0 (144,73), clim_olr (144,72), clim_pre (144,72), # aa1 (144,72), bb1 (144,72), # opi0 (144,72), opi1 (144,72), opi (144,72) c dimension olr (144,72) c c 1. to convert the OLR data from 144*73 array to 144*72 array c call get_olr ( olr0, olr ) c write (6,*) ' finished convert OLR data !!' c c 2. to define OPI0 c call def_opi0 ( olr, clim_olr, clim_pre, aa0, bb0, opi0 ) c write (6,*) ' finished defining OPI0 !!' c c 3. to define OPI1 c call def_opi1 ( olr, aa1, bb1, opi1 ) c write (6,*) ' finished defining OPI1 !!' c c 4. to combine OPI0 and OPI1 c call com_opis ( kmon, opi0, opi1, opi ) c write (6,*) ' finished combining OPIs !!' c c return end c c c ************************************************************************ c * * c * * c * ******* ******* ******* ******* * ******* * c * * * * * * * * * * c * * * * * * * * * * c * * **** ******* * * * * ******* * c * * * * * * * * * * * * c * * * * * * * * * * * c * ******* ******* * ******* ******* ******* * * * c * * c * * c ************************************************************************ c c c ************************************************************************ c * * c * subroutine : * c * * c * get_olr * c * * c * objective : * c * * c * to convert the OLR data from 144*73 array * c * to 144*72 array * c * * c * input : * c * * c * olr0 == original OLR data from ftp site * c * olr0 (i,j) * c * i = 1, 144 ==> (0.0, 2.5E, ... 2.5W) * c * j = 1, 73 ==> (90.0N, 87.5N, ... 90.0S) * c * * c * output : * c * * c * olr == converted OLR * c * olr (i,j) * c * i = 1, 144 ==> (1.25E, 3.75E, ... 1.25W) * c * j = 1, 72 ==> (88.75S,86.25S,... 88.75N)* c * * c ************************************************************************ c subroutine get_olr ( olr0, olr ) c dimension olr0 (144,72), olr (144,72) c c 1. to convert the data arrays for OLR c do 1002 jj=1,72 do 1002 ii=1,144 j1 = jj j2 = jj + 1 i1 = ii i2 = ii + 1 nn = 0 ss = 0.0 do 1001 j = j1, j2 do 1001 i = i1, i2 js = j if (i.le.144) then is = i else is = i - 144 end if if (olr0(is,js).ge.0.0) then nn = nn + 1 ss = ss + olr0(is,js) end if 1001 continue if (nn.gt.0) then olr(ii,73-jj) = ss/nn else olr(ii,73-jj) = -999.0 end if 1002 continue c c return end c c c ************************************************************************ c * * c * * c * ******* ******* ******* ******* ******* ******* ****** * c * * * * * * * * * * * * * c * * * * * * * * * * * * * c * * * ******* ******* * * ******* * * * * c * * * * * * * * * * * * c * * * * * * * * * * * * c * ******* ******* * ******* ******* * ******* ****** * c * * c * * c ************************************************************************ c c ************************************************************************ c * * c * subroutine : * c * * c * def_opi0 * c * * c * objective : * c * * c * to define the preciptation estimates from NOAA * c * OLR using the original OPI algorithm * c * * c * input : * c * * c * olr == total olr data * c * * c * clim_olr == olr climatology * c * * c * clim_pre == precip climatology * c * * c * aa0/bb0 == regressional coefficients between * c * proportional constant PC and mean * c * precipitation (R) * c * PC = aa0 * R + bb0 * c * * c * output : * c * * c * opi0 == precip estimates based on the * c * * c * algorithm : * c * * c * 1) to determine the proportional constant PC * c * between the precip anomaly and OLR anomaly * c * PC = aa0 * R + bb0 * c * * c * 2) to calculate the precip anomaly from OLR * c * anomaly * c * precip anom = PC * OLR anom * c * * c * 3) to define the total precip estimates as * c * the summation of anomaly and climatology * c * precip totl = precip anom + precip clim * c * * c ************************************************************************ c subroutine def_opi0 ( olr, clim_olr, clim_pre, # aa0, bb0 , opi0 ) c dimension olr (144,72), clim_olr (144,72), clim_pre (144,72), # opi0 (144,72) c dimension anom_olr (144,72), anom_pre (144,72) c c c 1. to calculate the OLR anomaly c do 1001 jj = 1, 72 do 1001 ii = 1, 144 if (olr (ii,jj).gt.0.0.and. # clim_olr(ii,jj).gt.0.0) then anom_olr (ii,jj) = olr (ii,jj) # - clim_olr (ii,jj) else anom_olr (ii,jj) = -999.0 end if 1001 continue c c 2. to define the precipitation estimates from OLR using the c original OPI algorithm c do 2001 jj = 1, 72 do 2001 ii = 1, 144 if (anom_olr(ii,jj).ne.-999.0.and. # clim_pre(ii,jj).ne.-999.0.and. # aa0 .ne.-999.0.and. # bb0 .ne.-999.0) then pc = aa0 * clim_pre (ii,jj) + bb0 anom_pre (ii,jj) = pc * anom_olr (ii,jj) opi0 (ii,jj) = clim_pre (ii,jj) # + anom_pre (ii,jj) if (opi0(ii,jj).lt.0.0) then opi0 (ii,jj) = 0.0 end if else anom_pre (ii,jj) = -999.0 opi0 (ii,jj) = -999.0 end if 2001 continue c c return end c c c ************************************************************************ c * * c * * c * ******* ******* ******* ******* ******* ******* * * c * * * * * * * * * * * * c * * * * * * * * * * * * c * * * ******* ******* * * ******* * * * c * * * * * * * * * * * c * * * * * * * * * * * c * ******* ******* * ******* ******* * ******* * * c * * c * * c ************************************************************************ c c ************************************************************************ c * * c * subroutine : * c * * c * def_opi1 * c * * c * objective : * c * * c * to define the total precip estimates from total * c * olr * c * * c * input : * c * * c * olr == total OLR (W/m^2) from NOAA polar * c * satellites * c * * c * aa1/bb1 == totl precip = aa1 * totl olr + bb1 * c * regressional coefficients * c * * c * output : * c * * c * opi1 == total precip estimates (mm/day) * c * from total OLR * c * * c ************************************************************************ c subroutine def_opi1 ( olr, aa1, bb1, opi1 ) c dimension olr (144,72), # aa1 (144,72), bb1 (144,72), # opi1 (144,72) c c 1. to destimate the total precipitation from total OLR c do 1001 jj = 1, 72 do 1001 ii = 1, 144 if (olr(ii,jj).gt.0.0.and. # aa1(ii,jj).ne.-999.0.and. # bb1(ii,jj).ne.-999.0) then opi1 (ii,jj) = aa1 (ii,jj) * olr (ii,jj) # + bb1 (ii,jj) if (opi1(ii,jj).lt.0.0) then opi1 (ii,jj) = 0.0 end if else opi1 (ii,jj) = -999.0 end if 1001 continue c c return end c c c ************************************************************************ c * * c * * c * ******* ******* * * ******* ******* ******* ******* * c * * * * ** ** * * * * * * * c * * * * * * * * * * * * * * * c * * * * * * * * * ******* * ******* * c * * * * * * * * * * * * c * * * * * * * * * * * * c * ******* ******* * * ******* ******* * ******* ******* * c * * c * * c ************************************************************************ c c ************************************************************************ c * * c * subroutine : * c * * c * com_opis * c * * c * objective : * c * * c * to combine the OPI0 and OPI1 * c * * c * input : * c * * c * kmon == calendar month of the target data * c * ( will be used to determine the weight ) * c * * c * opi0 == precip estimates by the original OPI * c * algorithm * c * * c * opi1 == precip estimates by the regressional * c * between total precip and total OLR * c * * c * opi == OPI estimates combined from OPI0 and * c * OPI1 * c * * c ************************************************************************ c subroutine com_opis ( kmon, opi0, opi1, opi ) c dimension opi0(144,72), opi1(144,72), opi (144,72) c dimension cor(144,72), # wei_lat (72), wei_opi (144,72) c dimension bsouth1(12), bsouth2(12), bnorth1(12), bnorth2(12) data bsouth1/-40.0,-40.0,-40.0,-37.0,-33.0,-30.0, # -30.0,-30.0,-30.0,-33.0,-37.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/ 30.0, 30.0, 30.0, 33.0, 37.0, 40.0, # 40.0, 40.0, 40.0, 37.0, 33.0, 40.0/ c parameter ( a = -49.0/1125.0, b = 49.0/45.0 ) c c 1. to define the weighting coefficients as a function of latitude c bs1 = bsouth1 (kmon) bs2 = bsouth2 (kmon) bn1 = bnorth1 (kmon) bn2 = bnorth2 (kmon) do 1001 jj=1,72 rlat = jj*2.5 - 91.25 if (rlat.le.bs1) then wei_lat(jj) = 0.0 else if (rlat.le.bs2) then wei_lat(jj) = (rlat-bs1)/(bs2-bs1) else if (rlat.le.bn1) then wei_lat(jj) = 1.0 else if (rlat.le.bn2) then wei_lat(jj) = (bn2-rlat)/(bn2-bn1) else wei_lat(jj) = 0.0 end if 1001 continue c c 2. to calculate the local pattern correlation between OPI0 and OPI1 c do 2001 jj = 1, 72 do 2001 ii = 1, 144 cor (ii,jj) = -999.0 2001 continue do 2003 jj = 1, 72 do 2003 ii = 1, 144 i1 = ii - 1 i2 = ii + 1 j1 = jj - 3 j2 = jj + 3 if (j1.le.0) then j1 = 1 end if if (j2.ge.72) then j2 = 72 end if nn = 0 sx = 0.0 sy = 0.0 sxx = 0.0 sxy = 0.0 syy = 0.0 do 2002 j=j1,j2 do 2002 i=i1,i2 js = j if (i.ge.1.and.i.le.144) then is = i else if (i.le.0) then is = i + 144 else is = i - 144 end if yy = opi1 (is,js) xx = opi0 (is,js) if (xx.ge.0.0.and.yy.ge.0.0) then nn = nn + 1 sx = sx + xx sy = sy + yy sxx = sxx + xx*xx sxy = sxy + xx*yy syy = syy + yy*yy end if 2002 continue if (nn.ge.10) then ax = sx / nn ay = sy / nn axx = sxx - nn*ax*ax axy = sxy - nn*ax*ay ayy = syy - nn*ay*ay if (axx.gt.0.0.and.ayy.gt.0.0) then cor(ii,jj) = axy / sqrt(axx*ayy) end if end if 2003 continue c c 3. to combine the OPI0 and OPI1 c do 3001 jj=1,72 do 3001 ii=1,144 c c = cor (ii,jj) c if (c.eq.-999.0) then c wei_opi (ii,jj) = 0.0 c else if (c.ge.0.8) then c wei_opi (ii,jj) = 0.0 c else if (c.ge.0.3) then c r = (1.0-c)*(1.0-c) c wei_opi (ii,jj) = a / r + b c else c wei_opi (ii,jj) = 1.0 c end if c = cor (ii,jj) if (c.eq.-999.0) then wei_opi (ii,jj) = 0.0 else if (c.ge.0.7) then wei_opi (ii,jj) = 0.0 else if (c.ge.0.3) then r = (1.0-c)*(1.0-c) wei_opi (ii,jj) = a / r + b else wei_opi (ii,jj) = 1.0 end if 3001 continue do 3002 jj = 1, 72 do 3002 ii = 1, 144 w1 = wei_lat (jj) * wei_opi (ii,jj) if (w1.lt.0.0) then w1 = 0.0 else if (w1.gt.1.0) then w1 = 1.0 end if w0 = 1.0 - w1 r0 = opi0 (ii,jj) r1 = opi1 (ii,jj) if (r0.ge.0.0.and.r1.ge.0.0) then opi (ii,jj) = r0*w0 +r1*w1 else if (r0.ge.0.0) then opi (ii,jj) = r0 else opi (ii,jj) = r1 end if 3002 continue c return end