PROGRAM looper ! ! Author: Paul A. Dirmeyer ! COLA/IGES ! Calverton, MD, USA ! October 1998 ! Revised Feb 2004 for new data sources, OMP implementation ! Revised Jul 2008 back to single thread for Linux ! USE forcing_mod ! ! INTEGER :: x0,x1,y0,y1,ii,jj,tt,t0,t1,ix0,iy0,t2 INTEGER :: isink ! Total number of grid boxes INTEGER :: ofil ! Output file unit number INTEGER :: jt ! Output file record number INTEGER :: lwindo ! Number of days in launch interval INTEGER :: ntsps ! number of trajectory periods in the year INTEGER :: year ! 4-digit year of job INTEGER :: pentst ! First pentad of job (typically 1) INTEGER :: pentend ! Last pentad of job (typically 73) INTEGER :: pentad ! Current pentad (FORTRAN code runs once per pentad) INTEGER :: slack ! max number days outside pentad to track parcels INTEGER :: juls ! Julian start date of reanalysis data files INTEGER :: leap ! Is this pentad after a leap day in a leap year (0/1)? INTEGER :: leapday ! Does this pentad contain a leap day (0/1)? INTEGER :: mstat ! Script execution status INTEGER :: xarec ! Current record in CMAP pentad precip. data INTEGER :: cmorphst ! Start record in CMORPH diurnal precip. data INTEGER :: pstep ! Number of hybrid precip time steps in pentad INTEGER :: ftmax ! Number of time steps in the reanalysis sample INTEGER :: partot ! Number of parcels to launch in pentad INTEGER :: pnmax ! # parcels to launch (at least 2*time rean steps) INTEGER :: maxcpu ! # parallel processes to run at one time INTEGER :: tid ! current thread number ! INTEGER :: ivec(ixmax*iymax) ! Vector of i coords of sink points INTEGER :: jvec(ixmax*iymax) ! Vector of j coords of sink points INTEGER :: limnorth,limsouth ! North and south latitudes of diurnal precip hybrid algorithm ! REAL :: lmask(ixmax,iymax) ! Land/sea mask (1=land, 0=water) REAL :: xapp(pxmax,pymax) ! Stdqc daily US gauge precip rate (mm/d) REAL :: hybp(ptmax,pxmax,pymax) ! Hybrid precip during pentad (mm/d - same as CMAP "truth") REAL :: cmorp(pxmax,pymax,ptmax) ! CMORPH diurnal precip during pentad (mm/s) REAL :: dlyp(pxmax,pymax,6) ! Daily rean precip during pentad (mm/s) REAL :: pwork(pxmax,pymax) ! Working array for precip REAL :: pwork6(pxmax,pymax,6) ! Working array for precip ! REAL :: tstep ! Number of time steps per rean data interval (<0 back traj) REAL :: stday REAL :: kickoff REAL :: pqmin ! Minimum threshold for parcel water content (%) REAL :: dtmin ! Trajectory integration time step (minutes) REAL :: dtt ! Trajectory integration time step (hours) REAL :: decang ! Solar declination (for computing range for applying CMORPH diurnal correction) REAL :: latnorth,latsouth ! North and south latitudes of diurnal precip hybrid algorithm REAL :: diulat ! Poleward extent of diurnal hyb. prec. on the equinox real :: tmr1,tmr2 ! CHARACTER(LEN=25) :: filnam CHARACTER(LEN=4) :: ayear CHARACTER(LEN=3) :: ax0, ay0 diulat = 30.0 ! !>>>>> Open data files ! !>>>>> File 1 = input parameters (ASCII file) ! OPEN (1,FILE='wcr.in',STATUS='OLD') ! READ(1,*) year,pentst,pentend,pentad,slack,dtmin,juls,leap,leapday,pqmin,partot,maxcpu,mstat dtt = dtmin / 60.0 ! !>>>>> Teens = parsed reanalysis fields ! OPEN (10,FILE='prec.grd', & FORM='UNFORMATTED',STATUS='OLD',ACCESS='DIRECT', & RECL=ixmax*iymax) OPEN (11,FILE='landsfc.grd', & FORM='UNFORMATTED',STATUS='OLD',ACCESS='DIRECT', & RECL=ixmax*iymax) OPEN (12,FILE='evap.grd', & FORM='UNFORMATTED',STATUS='OLD',ACCESS='DIRECT', & RECL=ixmax*iymax) OPEN (13,FILE='ps.grd', & FORM='UNFORMATTED',STATUS='OLD',ACCESS='DIRECT', & RECL=ixmax*iymax) OPEN (14,FILE='q0.grd', & FORM='UNFORMATTED',STATUS='OLD',ACCESS='DIRECT', & RECL=ixmax*iymax) OPEN (15,FILE='t0.grd', & FORM='UNFORMATTED',STATUS='OLD',ACCESS='DIRECT', & RECL=ixmax*iymax) OPEN (16,FILE='u0.grd', & FORM='UNFORMATTED',STATUS='OLD',ACCESS='DIRECT', & RECL=ixmax*iymax) OPEN (17,FILE='v0.grd', & FORM='UNFORMATTED',STATUS='OLD',ACCESS='DIRECT', & RECL=ixmax*iymax) ! !>>>>> Twenties = observational precipitation data sets ! OPEN (20,FILE='stdqc/stdqc_trp.grd', & FORM='UNFORMATTED',STATUS='OLD',ACCESS='DIRECT', & RECL=pxmax*pymax) OPEN (21,FILE='cmorph/cmorph_diurnal.grd', & FORM='UNFORMATTED',STATUS='OLD',ACCESS='DIRECT', & RECL=pxmax*pymax) ! !>>>>> Thirties+ are for the output files ! ! !>>>> Set critical arguments ! lwindo = 5 + leapday ! Number of days in launch interval pnmax = MAX(partot,lwindo*irtpd*2) ! # parcels to launch (at least 2*time rean steps) tstep = -REAL(24/(irtpd*dtt)) ! optimum magnitude=precip_data_int/reanal_data_int !!! kickoff = 20.0 ! Offset from start of first month of data to end of 1st lwindo ftmax = (lwindo + slack) * irtpd + 1 ! Number of rean time steps in files ! !>>>>> Memory allocation for reanalysis data (to be in shared memory for all threads ! ALLOCATE(prec(ixmax,iymax,ftmax)) ALLOCATE(evap(ixmax,iymax,ftmax)) ALLOCATE(psfc(ixmax,iymax,ftmax)) ALLOCATE(qhum(ixmax,iymax,izmax,ftmax)) ALLOCATE(tmpr(ixmax,iymax,izmax,ftmax)) ALLOCATE(uwnd(ixmax,iymax,izmax,ftmax)) ALLOCATE(vwnd(ixmax,iymax,izmax,ftmax)) ! !>>>>> Time navigation of the precipitation datasets ! xarec = (pentad - 1) * 5 + 1 + leap print *,"xarec=",xarec xapp = 0.0 DO ii = xarec,xarec+4 READ(20,REC=ii) pwork ! For mm/day input; divide by 5 by daily average within the pentad xapp = xapp + pwork * 1.0 / 5.0 ! For in/day input; divide by 5 by daily average within the pentad ! xapp = xapp + pwork * 25.4 / 5.0 ENDDO WRITE(6,3000)(jj,xapp(jj,58),jj=39,45) 3000 FORMAT(i4,2x,F11.6) cmorphst = (juls + slack - leap + leapday - 1) * ipdiur + 1 print *,juls,slack,leap,leapday cmorphst = MOD(cmorphst-1,365*ipdiur)+1 print *,'cmorphst=',cmorphst pstep = (5 + leapday) * ipdiur t0 = cmorphst t1 = t0 + pstep - 1 ii = 0 DO tt = t0,t1 ii = ii + 1 READ(21,REC=tt)cmorp(:,:,ii) END DO t0 = 1 t1 = ftmax ii = 0 DO tt = t0,t1 ii = ii + 1 READ(10,REC=tt)prec(:,:,ii) END DO ! !>>>>> Create hybrid 3-hour precip data for this pentad ! decang = 23.45 * sin((juls+2+slack+284)*2*pi/365) latnorth = decang*0.5+diulat latsouth = decang*0.5-diulat print *,'latnorth=',latnorth,', latsouth=',latsouth DO jj = 1, iymax IF (latsouth <= gauss(jj)) THEN limsouth = jj EXIT ENDIF ENDDO DO jj = iymax, 1, -1 IF (latnorth >= gauss(jj)) THEN limnorth = jj EXIT ENDIF ENDDO print *,'limnorth=',limnorth,', limsouth=',limsouth !d print *,xapp(166,32) !d print * !d print *,prec(166,32,61:80) !d print * !d print *,cmorp(166,32,:) ! !>>>>> First fill the global grid with interpolated reanalysis precipitation ! DO tt = 1, 5 + leapday DO ii = 1,ipdiur t0 = ii + ipdiur * (tt-1) t1 = (t0+1) / 2 + slack * irtpd t2 = t1 + 1 hybp(t0,:,:) = prec(:,:,t1) * 0.5 ! 0.5 because each rean ! timestep is applied twice (6h vs 3h time step) END DO END DO t0 = t0 + 1 t1 = t1 + 1 hybp(t0,:,:) = prec(:,:,t1) * 0.5 DO t0 = 1, ipdiur * (5 + leapday) t1 = t0 + 1 hybp(t0,:,:) = (hybp(t0,:,:) + hybp(t1,:,:)) * 0.5 ENDDO ! !>>>>> Go back through the tropics and replace with diurnal (CMORPH) based hybrid ! pwork6 = 0.0 DO tt = 1, 5 + leapday DO ii = 1,ipdiur t0 = ii + ipdiur * (tt-1) pwork6(:,:,tt) = pwork6(:,:,tt) + cmorp(:,:,t0) ! Daily sum of CMORPH P END DO END DO DO tt = 1, 5 + leapday t1 = irtpd * tt + slack * irtpd dlyp(:,:,tt) = prec(:,:,t1) + prec(:,:,t1-1) + prec(:,:,t1-2) + prec(:,:,t1-3) ! Daily sum of rean P END DO WHERE (pwork6 > 0.0) dlyp = dlyp / pwork6 ! P(daily) / sum[P(diurnal)] for each day of pentad DO tt = 1, 5 + leapday DO ii = 1,ipdiur t0 = ii + ipdiur * (tt-1) hybp(t0,:,limsouth:limnorth) = cmorp(:,limsouth:limnorth,t0) * dlyp(:,limsouth:limnorth,tt) END DO END DO ! !>>>>> At this point, each grid box contains the right local distribution of P in time within the pentad. !>>>>> Ex.Trop: P(intermed) = P(rean) interp to 3-hr !>>>>> Tropics: P(intermed) = P(diurn)*P(rean_daily)/sum[P(diurn)] !>>>>> Now the precipitation must be scaled by the Xie-Arkin pentad estimates to give proper magnitudes. ! pwork = 0.0 DO tt = 1, 5 + leapday DO ii = 1,ipdiur t0 = ii + ipdiur * (tt-1) pwork = pwork + hybp(t0,:,:) ! pentad sum of intermediate precip estimate END DO END DO WHERE (pwork > 0.0) xapp = xapp / pwork ! P(pentad) / sum[P(intermed)] ! !>>>>> Ex.Trop: P(hybrid) = P(rean)*P(pentad)/sum[P(rean)] !>>>>> Tropics: P(hybrid) = P(diurn)*P(rean_daily)*P(pentad)/sum[P(diurn)]/sum[P(rean)] ! DO tt = 1, 5 + leapday DO ii = 1,ipdiur t0 = ii + ipdiur * (tt-1) !PAD09/2006: This should be *5 as xapp is mm/d not mm/pentad !!! hybp(t0,:,:) = hybp(t0,:,:) * xapp END DO END DO ! ! !y WRITE(ayear,FMT='(i4)')year !y OPEN (44,FILE='/ptmp/dirmeyer/'//ayear//'/hybp.grd', & !y FORM='UNFORMATTED',STATUS='UNKNOWN',ACCESS='DIRECT', & !y RECL=4*pxmax*pymax) !y DO tt = 1, 5 + leapday !y DO ii = 1,ipdiur !y t0 = ii + ipdiur * (tt-1) !y WRITE(44,REC=t0)hybp(t0,:,:) !y END DO !y END DO !y CLOSE(44) !d print * !d print *,hybp(:,166,32) !d print * pwork = SUM(hybp,DIM=1)-hybp(41,:,:) !d print *,pwork(166,32) !c print * !c print *,xapp(144,69) !c print *,dlyp(144,69,:) !c print *,hybp(:,144,69) !c print * !c print *,pwork(144,69) !c print *,pwork6(144,69,:) !c print * !c pwork = SUM(hybp,DIM=1) !c print *,pwork(144,69) ! !>>>>> Read met forcing data into memory ! DO tt = 1,ftmax !! print *,"Reading time step ",tt READ(12,REC=tt)evap(:,:,tt) READ(13,REC=tt)psfc(:,:,tt) DO ii = 1,izmax jj = (tt - 1) * izmax + ii READ(14,REC=jj)qhum(:,:,ii,tt) READ(15,REC=jj)tmpr(:,:,ii,tt) READ(16,REC=jj)uwnd(:,:,ii,tt) READ(17,REC=jj)vwnd(:,:,ii,tt) ENDDO ENDDO ! !>>>>> Land mask ! READ(11,REC=1) lmask ! jt = pentad-pentst+1 jt = pentad stday = float(juls) basin = 0.0 WRITE(6,1100) pentad,jt,stday 1100 FORMAT('Starting pentad #',i3,', output record=',i3,', STDAY=',f6.0) ! !>>>>> Data are ready, Start SINK grid point loop ! ! ivec/jvec stores the coordinates of the points to be used by onebox ! ii = 0 DO y0 = domainy0, domainyf DO x0 = domainx0, domainxf IF (lmask(x0,y0) == 1.0) THEN ii = ii + 1 ivec(ii) = x0 jvec(ii) = y0 ENDIF ENDDO ENDDO isink = ii print *,'Total number of grid points =',isink ofil = 30 ! !///////////////////////////////////////////////////////////////////////////////////////// ! DO ii = 1, isink tid = 1 ofil = 30 + tid ix0 = ivec(ii) iy0 = jvec(ii) WRITE(ayear,FMT='(i4)')year WRITE(ax0,FMT='(i3.3)')ix0 WRITE(ay0,FMT='(i3.3)')iy0 filnam = 'source.' // ayear // '.x' // ax0 // '.y' // ay0 // '.grd' WRITE(6,1110) tid,ii,ax0,ay0 1110 FORMAT(i2,'> Point #',i6,', at X=',a3,', Y=',a3) print *,'output/',filnam OPEN (ofil,FILE='output/'//filnam, & FORM='UNFORMATTED',STATUS='UNKNOWN',ACCESS='DIRECT', & RECL=ixmax*iymax) ! ! Usage of onebox in onebox.f: ! SUBROUTINE onebox(thybp,ofil,orec,lwindo,x000,y000,tstep,ftmax,pstep,slack,pnmax,pqmin) ! CALL onebox(hybp(:,ix0,iy0),ofil,jt,lwindo,ix0,iy0,tstep,ftmax,pstep,slack,pnmax,pqmin) CLOSE (ofil) END DO ! !\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ ! print *,'CPUs used=',maxcpu CALL CPU_TIME(tmr1) print *,'Total CPU time (sank) =',tmr1 !O CALL setrteopts ('cpu_time_type=usertime') CALL CPU_TIME(tmr2) print *,'Total User time (sank) =',tmr2 ! !>>>>> Memory deallocation for reanalysis data ! DEALLOCATE(evap) DEALLOCATE(psfc) DEALLOCATE(qhum) DEALLOCATE(tmpr) DEALLOCATE(uwnd) DEALLOCATE(vwnd) ! END PROGRAM looper