!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! YZhang 2013-12-24 !!! ! Sample Fortran Program: read_o2_ISCCP_GWS.f ! (rvsd from read_ISCCP_GLWS.f) !--read the global, 3-hrly, ISCCP-D1 based Global Weather State (GWS) numbers ! (GWS number of 1 to 12) of reordered/current version 'o2'. ! ! It uses five subroutines: ! (1) SUBROUTINE MDYTOJ (MONTH,IDAY,IYEAR,JULIAN) ! (2) SUBROUTINE FEBDAY (Iyear,NDSFEB) ! (3) SUBROUTINE JYTMDP (JULIAN,IYEAR,MONTH,IDAY,PRJDAT) ! (4) SUBROUTINE EQAR25 ( DLON,DLAT) ! (5) SUBROUTINE JRPE25 (IOPTN,ARYBOX,EQANGL) ! !-- Users can customize the program to their needs ! ! Import: ! Current (o2-version) Global Weather States (GWS) numbers (K=1->12) of global, ! 3-hrly, 280-km-equal-area maps (EQA) as specified below, ! ! /your_dir/o2_dw_GL_8309_43.YYYY ! ! where YYYY = year (1983--2009) ! ! The Global Weather States (GWS) values are specified as: ! ! Ocean GWS #: 1 -> 12 for 12 WS' ( 12 = clear sky); undefined = -999 ! Land GWS #: 101 -> 112 for 12 WS' (112 = clear sky); undefined = -999 ! ! The data files are in 'direct access'/'unformatted' binary; each record ! is for one Julian day on 280-km equal-area map for 8 UTC's, ! (UTC = 00h, 03h, 06h, .. 21h, every third hour), in Integer*2: Ind1ws(6596,8). ! Each data file is for a year (but July to December for 1983 only). ! ! The direct-access record lenghth is ! RecL = (6596*8*2) bytes/4=105536/4=26384 for 4-byte-based RecL, ! (if your OS uses 1-byte-based RecL, it is then RecL=105536) ! ! ! Export: ! global maps of GWS numbers = 1 --> 12, where #12 is for clear-sky state, ! for all land and ocean areas ! * Users must specify (or comment out what they don't need): !(1) if IPRINT=1, will printout WS # for a desiganated latitude & longitude, * prtlat/prtlon = latitude/longitude (prtlat=-90->90; prtlon=0->360 E) * of which the cell's WS # can be printed out * on screen for monitoring/checking *(2) if IWRITE = 1, will write binary output files for a global map ! (S->N; Eastwrad) for year/month/date/UTC = iyou/imou/idou/ihou ! in one of three maps as specified by IOPOUT: * IOPOUT = option for three different output maps: * 0 equal-area (EQA) map of interger*4(6596) in gwsisccp_YYMMDDHH * 1 square lat/lon (SQD) map of interger*4(144,72) in gwsisccp.YYMMDDHH ! (longitude starts from dateline) * 2 square lat/lon (SQG) map of interger*4(144,72) in gwsisccp-YYMMDDHH ! (longitude starts from Greenwich line) ! ! for a specified year, month, date and UTC (iyou/imou/idou/ihou), YYMMDDHH, ! e.g., YYMMDDHH=99071521 for UTC=21 of 7/15, 1999; ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Note: LrecL is the 4-byte based record length for direct-access ! if your machine uses 1-byte based length, multiply LrecL by 4 ! parameter (nclus=12, Iprint=0,Iwrite=1,iopout=1, LrecL=26384) parameter (iyr1=1983,iyr2=2009,nmon=12,ndays=31,ngmt=8, & maxlon=144,maxlat=72,maxbox=6596,iundef=-999, & prtlat=31.25, prtlon=271.25, & iyou=1983,imou= 7,idou=15,ihou=12) parameter (ibx1=1,ibx2=6596, & alat1=-90.,alat2=90.,alon1=0.,alon2=360., & nlat=72,lat1=1,lat2=72,lon1=1,lon2=144) Integer*2 Ind1ws(ibx1:ibx2,ngmt),MONDAY(12) INTEGER Ieq(6596),Iea(144,72),jeq(maxbox) logical Lexist character ufin*80,ufou*80,ft*8 INTEGER*2 BOXNMB,LONEQR,LATEQR,NCELLS common/eqreqa/ BOXNMB(144,72),LONEQR(6596),LATEQR(6596), & NCELLS(72),DLONTB(72),XLATB(72),XLATE(72) !.......................................................... !: Speifying input data file location data UFIN / &'/Users/yzhang11/WS/d1__ws/o2_dw_GL_8309_43.YYYY'/ ! 26 & ipin/26/ !: Speifying output data file location data ufou/ &'/Users/yzhang11/WS/temp/gwsisccp*YYMMDDHH'/, ! 24 & ipou/24/ ! data monday /31,28,31,30,31,30,31,31,30,31,30,31/ !------------------------------------------------------------- ! Get 280-km (2.5 X 2.5 degree on Equator) equal area map info: ! CALL EQAR25(2.5,2.5) ! ! do loop for all the years: ! DO 377 iyr = iyr1, iyr2 write(ufin(ipin+18:ipin+21),'(i4.4)')iyr inquire(file=ufin, exist=Lexist) IF(Lexist .NEQV. .true.)THEN PRINT *,'................................................' PRINT *,'---- NON-EXISTENT FILE: ',ufin stop 1111 END IF print *,'---> input 11: ',ufin open(unit=11,file=ufin, * status='old',access='direct',recl=LrecL, !'GL' & form='unformatted',iostat=ios) if(iyr .eq. 1983)then km1= 7 else km1=1 end if km2=12 ! krec=0 ! do loop for all the months: ! DO 366 km = km1, km2 !! imn=km if(iyr .ge. 2000)then iy=iyr - 2000 else iy=iyr - 1900 end if call MDYTOJ(Imn,1,IYR,JD1) !get serial # = Julian date for a WHOLE year mdays=monday(imn) if(imn.eq.2)call FEBDAY(iyr,mdays) call MDYTOJ(Imn,mdays,IYR,JD2) if(iyr.eq.1983)then call MDYTOJ(6,30,IYR,JD0) jd1=jd1-jd0 !get serial record # for 1983 that is for July to Dec jd2=jd2-jd0 end if ! do loop for all the Julian dates of a month for the year (iyr) do 344 jd = jd1, jd2 read(11,rec=jd)((Ind1ws(i,it),i=ibx1,ibx2),it=1,ngmt) krec=krec+1 jjdd=jd if(iyr.eq.1983)jjdd=jd+jd0 call JYTMDP (JJDD,IYR,MONTH,IDAY,IPRJDAT) do 311 it=1, ngmt ih=(it-1)*3 write(ft(1:8),'(4i2.2)')iy,month,iday,ih Ieq(:)=Iundef do ibx = 1, maxbox Ieq(ibx)=Ind1ws(ibx,it) ! print *,'Ieq(',ibx,')=',ieq(ibx) end do ! jeq(:)=iundef iea(:,:)=iundef do 309 ibx = 1, maxbox if(Ieq(ibx).ne.iundef)then ! Ocean GWS #: 1 -> 12 for 12 WS' ( 12 = clear sky); undefined = -999 ! Land GWS #: 101 -> 112 for 12 WS' (112 = clear sky); undefined = -999 ! Here gives merged (Land + Ocean) GWS number = 1-->12 if(Ieq(ibx) .ge. 100)Ieq(ibx)=Ieq(ibx)-100 jeq(ibx)=Ieq(ibx) end if 309 continue ! ! print out a cell's GWS value: ! IF(Iprint .EQ. 1)then lat=(prtlat+90.)/2.5 + 1 if(lat .gt. maxlat)lat=maxlat lon=prtlon/dlontb(lat)+1 if(lon .gt. ncells(lat))lon=ncells(lat) jbx=boxnmb(lon,lat) ! printout WS # of the cell for the designated longitude/latitude write(6,6221)jeq(jbx),prtlon,prtlat,jbx,iyr,month,iday,ih 6221 format('WS = ',I4,' .. for prtlon/prtlat = ',f6.2,'/',f6.2, & ' & eq-area box #=',i4,/, & 'for Year/month/date/UTC = ',i4,3('/',i2)) pause '6221' END IF ! ! output binary global WS map for designated yy/mm/dd/hh ! = iyou/imou/idou/ihou and in the file types for iopout = 0, or 1, or 2 ! IF(IWRITE .NE. 1)GO TO 311 if(iyou.eq.iyr .and. imou.eq.month .and. & idou.eq.iday .and. ihou.eq.ih)then print *,'---> output binary file for Year/month/date/UTC = ', & iyou,'/',imou,'/',idou,'/',ihou,' :' if(iopout .eq. 0)then ! output EQA ufou(ipou+9:ipou+17)='_'//ft(1:8) PRINT *,'++++ 21 ++++ UFOU = ',UFOU OPEN(UNIT=21,FILE=Ufou, * STATUS='unknown',ACCESS='sequential',FORM='UNFORMATTED', * IOSTAT=IOS) WRITE(21)jeq close (21) else if(iopout .eq. 1)then ! output SQD c: produce center longitude/latitude of each cell for SQD call jrpe25(1,jeq,iea) ufou(ipou+9:ipou+17)='.'//ft(1:8) PRINT *,'++++ 21 ++++ UFOU = ',UFOU OPEN(UNIT=21,FILE=Ufou, * STATUS='unknown',ACCESS='sequential',FORM='UNFORMATTED', * IOSTAT=IOS) WRITE(21)iea close (21) else if(iopout .eq. 2)then ! output SQG c: produce center longitude/latitude of each cell for SQG call jrpe25(0,jeq,iea) ufou(ipou+9:ipou+17)='-'//ft(1:8) PRINT *,'++++ 21 ++++ UFOU = ',UFOU OPEN(UNIT=21,FILE=Ufou, * STATUS='unknown',ACCESS='sequential',FORM='UNFORMATTED', * IOSTAT=IOS) WRITE(21)iea close (21) end if end if 311 continue ! 344 continue ! 366 continue ! print *,krec,' records have been read for year = ',iyr ! close (11) ! 377 continue ! stop end C===================================================================== ****************************************************************** * CONVERT MONTH/DAY/YEAR TO JULIAN (1 TO 365 OR 366) ****************************************************************** SUBROUTINE MDYTOJ (MONTH,IDAY,IYEAR,JULIAN) INTEGER PRJDAT DIMENSION MNTAB(12) DATA MNTAB /31,28,31,30,31,30,31,31,30,31,30,31/ C----------------------------------------------------------------- IF((MONTH.LT.1).OR.(MONTH.GT.12)) GO TO 500 MNTAB(2)=28 IF(MOD(IYEAR,4).EQ.0)THEN MNTAB(2)=29 IF(MOD(IYEAR,100).EQ.0)THEN IF(MOD(IYEAR,400).NE.0)MNTAB(2)=28 END IF END IF IF((IDAY.LT.1).OR.(IDAY.GT.MNTAB(MONTH))) GO TO 500 JULIAN=0 DO 100 IMON=1,12 IF(IMON.EQ.MONTH) GO TO 200 JULIAN=JULIAN+MNTAB(IMON) 100 CONTINUE 200 JULIAN=JULIAN+IDAY COR PRINT 301,MONTH,IDAY,IYEAR,JULIAN 301 FORMAT(/,2X,'MONTH DAY YEAR',3I6,2X,'JULIAN',I6,/) RETURN C--- C DATE ERROR C--- 500 PRINT 501,MONTH,IDAY,IYEAR 501 FORMAT(/,10X,'CNVJUL ERROR MONTH DAY YEAR',3I8,/) JULIAN=0 RETURN END C======================================================================= * Import Iyear = YYYY and get the dates of Februry -- 01-03-02 *................................................................... SUBROUTINE FEBDAY(Iyear,NDSFEB) NDSFEB=28 IF(MOD(IYEAR,4).EQ.0)THEN NDSFEB=29 IF(MOD(IYEAR,100).EQ.0)THEN IF(MOD(IYEAR,400).NE.0)NDSFEB=28 END IF END IF return end C===================================================================== C*********************************************************************** C--- DATE CONVERSION ROUTINES (JULCNV, CNVPRJ, & CNVJUL) C JULIAN DAY TO MONTH-DAY CONVERSION ROUTINE ! SUBROUTINE JYTMDP (JULIAN,IYEAR,MONTH,IDAY,PRJDAT) ! DIMENSION MNTAB(12) INTEGER PRJDAT DATA MNTAB /31,28,31,30,31,30,31,31,30,31,30,31/ C----------------------------------------------------------------------- C PRINT *,' JULIAN DATE OF',IYEAR,' ----> ',JULIAN MNTAB(2) = 28 IF(MOD(IYEAR,4).EQ.0)THEN MNTAB(2)=29 IF(MOD(IYEAR,100).EQ.0)THEN IF(MOD(IYEAR,400).NE.0)MNTAB(2)=28 END IF END IF MONTH = 0 IDAY = 0 IF(JULIAN.EQ.0) GO TO 800 ISUM = 0 DO 500 MNTH = 1,12 ISUM = ISUM + MNTAB(MNTH) IF(JULIAN.GT.ISUM) GO TO 500 MONTH = MNTH ISUM = ISUM - MNTAB(MNTH) IDAY = JULIAN - ISUM C----------------------------------------------------------------------- C COMPUTE PROJECT DATE C----------------------------------------------------------------------- PRJDAT = -365 - 181 DO 100 IYR = 1982,IYEAR-1 JUL = 365 IF(MOD(IYR,4).EQ.0) JUL = 366 PRJDAT = PRJDAT + JUL 100 CONTINUE PRJDAT = PRJDAT + JULIAN C PRINT *,'M/D/Y = ',MONTH,IDAY,IYEAR RETURN 500 CONTINUE C----------------------------------------------------------------------- C JULIAN DATE ERROR C----------------------------------------------------------------------- 800 PRINT 801,JULIAN,IYEAR 801 FORMAT(/,10X,'JULCNV ERROR JULIAN',I5,2X,'YEAR',I8,/) PRJDAT = -1 RETURN END C===================================================================== C DLON AND DLAT CAN BE ANY VALUES SPECIFIED IN MAIN PROGRAM. C--- NOTE: DIMENSION: 1:# SUBROUTINE EQAR25 ( DLON,DLAT) C C EQUAL-AREA COMPUTATIONS C ** ** REAL DLAT, DLON INTEGER LATMAX, LONMAX C PARAMETER ( DLON=0.5, DLAT=0.5 ) PARAMETER ( LONMIN=1, LATMIN=1 ) C PARAMETER ( LONMAX=360.0/DLON, LATMAX=180.0/DLAT ) ** ** INTEGER*2 BOXNMB,LONEQR,LATEQR,NCELLS common/eqreqa/ BOXNMB(144,72),LONEQR(6596),LATEQR(6596), & NCELLS(72),DLONTB(72),XLATB(72),XLATE(72) ** ** DATA RE /6371.2/ C C BEGIN LONMAX=360.0/DLON LATMAX=180.0/DLAT C PI=2.*ASIN(1.) TWOPI=2.*PI TWOPIR=TWOPI/360. RDLAT=DLAT*TWOPIR HEZON=RE*SIN(RDLAT) AEZON=TWOPI*RE*HEZON co AECELL=(AEZON*DLAT)/360. AECELL=(AEZON*DLON)/360. C PRINT 10,DLAT,LATMAX,AECELL 10 FORMAT(/,2X,'DLAT',F8.2,2X,'LATMAX',I5,2X,'AREA EQ CELL', 1 F8.1,/) C DO 100 LAT=LATMIN, LATMAX XLATB( LAT )=DLAT*(LAT-1)-90.0 XLATE( LAT )=XLATB( LAT )+DLAT IF(XLATE( LAT ).GT.90.) XLATE( LAT )=90. IF(XLATB( LAT ).GT.90.) XLATB( LAT )=90. RLATB=TWOPIR*XLATB( LAT ) RLATE=TWOPIR*XLATE( LAT ) HTB=RE*SIN(RLATB) HTE=RE*SIN(RLATE) HTZONE=ABS(HTE-HTB) AZONE=TWOPI*RE*HTZONE CELLS=AZONE/AECELL NCELLS( LAT )=CELLS+0.5 100 CONTINUE C LATLFT=LATMIN IF ( MOD( LATMAX-LATMIN+1, 2 ) .EQ. 1 ) THEN LATRIT=LATMAX-1 ELSE LATRIT=LATMAX ENDIF DO 200 LAT=LATMIN, LATMAX IF ( LATLFT .GE. LATRIT ) GOTO 201 NCELLS( LATRIT )=NCELLS( LATLFT ) XLATB( LATRIT )=-XLATE( LATLFT ) XLATE( LATRIT )=-XLATB( LATLFT ) LATLFT=LATLFT+1 LATRIT=LATRIT-1 200 CONTINUE 201 CONTINUE C C WRITE(31,3000) 3000 FORMAT(1X,'LAT',5X,'XLATB',5X,'XLATE',5X,'NCELLS',5X,'DLONTB') DO 300 LAT=LATMIN, LATMAX IF ( NCELLS( LAT ) .GT. 0 ) THEN DLONTB( LAT )=360./NCELLS( LAT ) ELSE DLONTB( LAT )=360.0 ENDIF IF(XLATE( LAT ).GT.90.) XLATE( LAT )=90. C WRITE(21,301) LAT,XLATB( LAT ),XLATE( LAT ), C & NCELLS( LAT ),DLONTB( LAT ) 301 FORMAT(1X,'LAT',I4,1X,'XLAT',2F8.1,2X,'CELLS', 1 I8,2X,'DLON',F8.1) 300 CONTINUE c:: Added from cveq25() DO 401 J=1,72 DO 401 I=1,144 BOXNMB(I,J)=-1000 401 CONTINUE KOUNT=0 DO 411 J=1,72 DO 411 I=1,NCELLS(J) KOUNT=KOUNT+1 BOXNMB(I,J)=KOUNT LONEQR(KOUNT)=I LATEQR(KOUNT)=J 411 CONTINUE CC IF(KOUNT .NE. 6596)STOP 4444 c:: RETURN END C===================================================================== *------------------------------------------------------------------- * REPLICATING INTEGER*4 ARRAY(6596) TO (144,72) STARTING * AT EITHER GREENWICH LINE (IOPTN=0) OR DATELINE (IOPTN=1)! *------------------------------------------------------------------- SUBROUTINE JRPE25(IOPTN,ARYBOX,EQANGL) INTEGER ARYBOX(6596),EQANGL(144,72) ! INTEGER*2 BOXNMB(144,72),LONEQR(6596),LATEQR(6596),NCELLS(72) ! REAL DLONTB(72),XLATB(72),XLATE(72) INTEGER*2 BOXNMB,LONEQR,LATEQR,NCELLS common/eqreqa/ BOXNMB(144,72),LONEQR(6596),LATEQR(6596), & NCELLS(72),DLONTB(72),XLATB(72),XLATE(72) C--- DO 101 J=1,72 DO 101 I=1,144 EQANGL(I,J)=-999 101 CONTINUE DO 111 J=1,72 LAT=J DO 111 I=1,144 IF(IOPTN .EQ.0)THEN IADJ=I ELSE IF(IOPTN .EQ. 1)THEN IADJ=I+72 IF(IADJ.GT.144)IADJ=IADJ-144 END IF LON=((FLOAT(I)-0.5)*2.5)/DLONTB(LAT) + 1 IF(LON .GT. NCELLS(LAT))LON=NCELLS(LAT) EQANGL(IADJ,J)=ARYBOX(BOXNMB(LON,J)) 111 CONTINUE RETURN END