PROGRAM testdods !---------------------------------------------------------------------- ! Purpose: test accessability of GSWP2 data served from COLA GDS (GrADs ! DODS Server) ! Zhichang Guo, April, 2003 !---------------------------------------------------------------------- USE netcdf IMPLICIT NONE INTEGER,PARAMETER :: lpnts = 15238 INTEGER,PARAMETER :: strlen_max = 128 INTEGER,PARAMETER :: adtt = 10800 INTEGER,PARAMETER :: varnums = 9 INTEGER :: year, month, day, hour, year0, month0, day0, hour0, months REAL :: time0, time1 INTEGER :: i, iret, varnum, varid INTEGER :: record0, record1, record2 INTEGER :: sec_met, sec_lag, julian INTEGER :: iyr, imn, irc, days, drc INTEGER :: drv_readnf2 INTEGER :: ncids(varnums), varids(varnums) REAL, ALLOCATABLE :: var2d(:,:) CHARACTER(LEN=strlen_max) :: baseurl, fnames(varnums), urls(varnums), fname CHARACTER(LEN=20) :: vnames(varnums) LOGICAL :: LEAPYR TYPE :: time_pointer INTEGER :: year ! INTEGER :: month ! INTEGER :: day ! INTEGER :: hour ! INTEGER :: minute ! INTEGER :: second ! END TYPE time_pointer TYPE(time_pointer) :: time_cur, time_met, time_tmp INTEGER, DIMENSION(12) :: monlen= (/ 31,28,31,30,31,30,31,31,30,31,30,31 /) DATA fnames/'tair_cru', 'wind_ncep' ,'swdown_srb','lwdown_srb','snowf_gswp','rainf_gswp' & ,'qair_ncep','rainf_c_gswp','psurf_ecor'/ DATA vnames/'tair', 'wind', 'swdown', 'lwdown' ,'snowf', 'rainf' & ,'qair', 'rainf_c', 'psurf'/ ! ------------------ parameters ------------- baseurl = 'http://www.monsoondata.org:9090/dods/gswp/vector/' varnum = 9 ! number of variables days = 2 ! how many days data will be read for each access transaction year0 = 1982 ! The first year to be read month0 = 7 ! The first month to be read day0 = 1 ! The first day to be read hour0 = 0 ! The first hour to be read months = 1 ! How many months will be read ! ------------------ open urls ------------- day = day0 hour = hour0 DO i = 1, varnum fname = baseurl(1:len_trim(baseurl))//fnames(i)(1:len_trim(fnames(i))) urls(i) = fname ! WRITE(*,'(A14,A66)') ' Open file: ',fname(1:66) iret = NF90_OPEN(fname,NF90_NOWRITE,ncids(i)) CALL check_err(iret) iret = NF90_INQ_VARID(ncids(i), vnames(i)(1:len_trim(vnames(i))), varids(i)) CALL check_err(iret) ! WRITE(*,'(A20,A12,A22,2I6)') ' Found variable ',vnames(i)(1:12),' with varid and ncid: ',varids(i),ncids(i) ENDDO time_met = time_pointer(1982,07,01,0,0,0) ! ------------------ read data ------------- year = year0 month = month0-1 drc = days * 86400 / adtt ALLOCATE(var2d(lpnts, drc)) DO imn = month0, month0+months-1 month = month + 1 IF(month > 12) THEN year = year + 1 month = month - 12 ENDIF monlen(2) = 28 IF(LEAPYR(iyr)) monlen(2) = 29 time_cur = time_pointer(year,month,day,hour,0,0) CALL drv_timestamp(time_met,time_cur,sec_met,julian,'S') record0 = sec_met/adtt iret = NF90_INQ_VARID(ncids(1), 'time', varid) CALL check_err(iret) iret = NF90_GET_VAR(ncids(1), varid, time0, start = (/1/)) CALL check_err(iret) iret = NF90_GET_VAR(ncids(1), varid, time1, start = (/record0+1/)) CALL check_err(iret) sec_lag = INT((time1-time0)*86400) CALL drv_timestamp(time_met,time_tmp,sec_lag,julian,'D') IF(time_tmp%year /= time_cur%year .OR. time_tmp%month /= time_cur%month .OR. & time_tmp%day /= time_cur%day .OR. time_tmp%hour /= time_cur%hour) THEN WRITE(*,*) time1-time0, record0, time_tmp%year, time_cur%year, time_tmp%month, time_cur%month WRITE(*,*) time_tmp%day, time_cur%day, time_tmp%hour, time_cur%hour STOP 'Met. forcing data date does not match' ENDIF WRITE(*,'(A47,4I6,A4,I10)') ' First met. record for the date yy/mm/dd/hh: ' & ,year0,month0,day0,hour0,' is ',record0+1 record1 = 1 record2 = monlen(month)*86400/adtt WRITE(*,'(A37,2I6)') ' Read data for the year and month: ',year,month DO irc = record0+record1, record0+record2, drc WRITE(*,'(A30,I8,A4,I8)') ' Reading data from record ',irc,' to ',irc+drc-1 DO i = 1, varnum WRITE(*,'(A24,A12,A11,A80)') ' Reading variable ',vnames(i)(1:12),' from URL: ', urls(i)(1:80) iret = drv_readnf2(ncids(i), varids(i), irc, drc, var2d, 1, lpnts) CALL check_err(iret) ENDDO ! variables loop ENDDO ! record loop ENDDO ! month loop DEALLOCATE(var2d) DO i = 1, varnum iret = NF90_CLOSE(ncids(i)) CALL check_err(iret) ENDDO STOP 'Program ended normally' END PROGRAM testdods !---------------------------------------------------------------------------------------- FUNCTION drv_readnf2(ncid, varid, record_beg, record_len, var, x_dim0, lpnts) ! ! This routine reads 2-dimensional varibale from a NetCDF file ! USE netcdf IMPLICIT NONE ! INTEGER :: drv_readnf2 ! Type of the function INTEGER :: ncid ! NCID of field INTEGER :: varid ! VARID of field INTEGER, INTENT(IN) :: record_beg ! The beginning record to read INTEGER, INTENT(IN) :: record_len ! The record length to read INTEGER, INTENT(IN) :: x_dim0, lpnts ! Dimension of var REAL, INTENT(OUT) :: var(lpnts,record_len) ! Variable read and returned ! ! LOCAL ! INTEGER :: numDims ! Number of dimensions for variable tmpVar INTEGER :: numADims ! Number of actual dimensions for variable tmpVar INTEGER, DIMENSION(4) :: dimIDs ! Dimension IDs for local variable tmpVar INTEGER, DIMENSION(4) :: lenDims ! Length for each dimension INTEGER, DIMENSION(4) :: lenADims ! Actual length for each dimension INTEGER :: i, n, iret, istart(3), icount(3) ! istart = 1 icount = 1 iret = NF90_INQUIRE_VARIABLE(ncid,varid,ndims = numDims, dimids = dimIDs) IF (iret .NE. NF90_NOERR ) THEN WRITE (*,*) 'readnf2: Could not inquire dimensions of NetCDF variable ', & & varid,' from file ',ncid STOP 'drv_readnf2' ENDIF numADims = 0 n = 0 DO i = 1, numDims iret = NF90_INQUIRE_DIMENSION(ncid, dimIDs(i), len = lenDims(i)) IF (iret .NE. NF90_NOERR) THEN WRITE (*,*) 'readnf2: Could not inquire dimension length of NetCDF variable ', & & varid,' from file ',ncid STOP 'drv_readnf2' ENDIF IF (lenDims(i) > 1) THEN n = n + 1 numADims = numADims + 1 lenADims(n) = lenDims(i) ENDIF ENDDO IF (numADims == 2) THEN istart(1) = x_dim0 icount(1) = lpnts istart(numDims) = record_beg icount(numDims) = record_len iret = NF90_GET_VAR(ncid,varid,var,istart,icount) IF (iret .NE. NF90_NOERR ) THEN WRITE (*,1000) istart,icount,varid,ncid 1000 FORMAT('readnf2: Could not read istart=',2i6,', icount=',2i6,' for NetCDF variable ',i4,' from file ',i4) STOP 'drv_readnf2 ' ENDIF ELSE WRITE (*,*) 'readnf2: this NetCDF variable is not a 2-D matrix ',ncid,varid,numADims STOP 'drv_readnf2' ENDIF drv_readnf2 = iret END FUNCTION drv_readnf2 !--------------------------------------------------------------------- SUBROUTINE check_err(iret) USE netcdf IMPLICIT NONE INTEGER :: iret IF(iret /= NF90_NOERR) THEN PRINT*, NF90_STRERROR(iret) STOP ENDIF END SUBROUTINE check_err !---------------------------------------------------------------------------------------- function LEAPYR(year) integer :: year logical :: LEAPYR LEAPYR = .false. if(mod(year,400) == 0) then LEAPYR = .true. else if(mod(year,4) == 0 .and. mod(year,100) /= 0) then LEAPYR = .true. endif endif end function LEAPYR !--------------------------------------------------------------------- SUBROUTINE drv_timestamp(time_origin,time_cur,seconds,julian,flag) !====================================================================== ! Description: Convert between Y/M/D H:M:S and seconds from origin ! flag = 's' convert date to seconds and julian ! flag = 'd' convert seconds to date and julian ! Note: This only works for 1901-2099! !====================================================================== IMPLICIT NONE TYPE :: time_pointer INTEGER :: year ! INTEGER :: month ! INTEGER :: day ! INTEGER :: hour ! INTEGER :: minute ! INTEGER :: second ! END TYPE time_pointer TYPE(time_pointer) :: time_origin, time_cur INTEGER, INTENT(INOUT) :: seconds ! seconds from origin time INTEGER, INTENT(OUT) :: julian ! julian day of year CHARACTER(LEN=1), INTENT(IN):: flag !** Other internal variables INTEGER :: d_iyr, d_mth, d_day, d_nhr, d_min, d_sec INTEGER :: leap, s_iyr, e_iyr, s_mth, e_mth, off, msum INTEGER :: i INTEGER, DIMENSION(12) :: monlen= (/ 31,28,31,30,31,30,31,31,30,31,30,31 /) INTEGER :: work !***** Given date, find seconds from origin IF (flag == 's' .or. flag == 'S') THEN seconds = 0 !*** Second, Minute, Hour, Day, Month d_sec = time_cur%second - time_origin%second d_min = (time_cur%minute - time_origin%minute)*60 d_nhr = (time_cur%hour - time_origin%hour)*3600 d_day = (time_cur%day - time_origin%day)*86400 d_mth = time_cur%month - time_origin%month ! msum = 0 IF (d_mth > 0) THEN DO i = 1, d_mth off = MOD(time_origin%month+i-2, 12) + 1 msum = msum + monlen(off) ENDDO ELSE IF (d_mth < 0) THEN DO i = d_mth, -1 off = MOD(time_origin%month+i+11, 12) + 1 msum = msum - monlen(off) ENDDO ENDIF d_mth = msum * 86400 !*** Year leap = 0 IF (time_origin%year == time_cur%year) THEN IF (MOD(time_origin%year, 4) == 0) THEN IF (time_origin%month < 3 .AND. time_cur%month > 2) leap = 1 IF (time_origin%month > 2 .AND. time_cur%month < 3) leap = -1 ENDIF ELSE s_iyr = time_origin%year e_iyr = time_cur%year s_mth = time_origin%month e_mth = time_cur%month IF (time_origin%year > time_cur%year) THEN s_iyr = time_cur%year e_iyr = time_origin%year s_mth = time_cur%month e_mth = time_origin%month ENDIF IF (MOD(s_iyr,4) == 0 .AND. s_mth < 3) leap = leap + 1 IF (MOD(e_iyr,4) == 0 .AND. e_mth > 2) leap = leap + 1 IF (e_iyr-s_iyr > 1) THEN DO i=s_iyr+1,e_iyr-1 IF (MOD(i,4) == 0) leap = leap + 1 ENDDO ENDIF IF (time_origin%year > time_cur%year) leap = -leap ENDIF ! d_iyr = ((time_cur%year - time_origin%year) * 365 + leap) * 86400 seconds = d_iyr + d_mth + d_day + d_nhr + d_min + d_sec !***** Given seconds from origin, find date ELSEIF (flag == 'd' .OR. flag == 'D') THEN !*** Second work = seconds d_sec = MOD(work,60) !*** Minute work = (work - d_sec)/60 d_min = MOD(work,60) !*** Hour work = (work - d_min)/60 d_nhr = MOD(work,24) !*** Day work = (work - d_nhr)/24 d_day = work !*** Second part 2 time_cur%second = time_origin%second + d_sec IF (time_cur%second < 0) THEN time_cur%second = time_cur%second + 60 d_min = d_min - 1 ENDIF IF (time_cur%second >= 60) THEN time_cur%second = time_cur%second - 60 d_min = d_min + 1 ENDIF !*** Minute part 2 time_cur%minute = time_origin%minute + d_min IF (time_cur%minute < 0) THEN time_cur%minute = time_cur%minute + 60 d_nhr = d_nhr - 1 ENDIF IF (time_cur%minute >= 60) THEN time_cur%minute = time_cur%minute - 60 d_nhr = d_nhr + 1 ENDIF !*** Hour part 2 time_cur%hour = time_origin%hour + d_nhr IF (time_cur%hour < 0) THEN time_cur%hour = time_cur%hour + 24 d_day = d_day - 1 ENDIF IF (time_cur%hour >= 24) THEN time_cur%hour = time_cur%hour - 24 d_day = d_day + 1 ENDIF !*** Day, Month, Year leap = 0 time_cur%day = time_origin%day time_cur%month = time_origin%month time_cur%year = time_origin%year IF (d_day > 0) THEN DO WHILE (d_day > 0) d_day = d_day - monlen(time_cur%month) IF (time_cur%month == 2 .AND. MOD(time_cur%year,4) == 0) THEN d_day = d_day - 1 leap = leap + 1 ENDIF time_cur%month = time_cur%month + 1 IF (time_cur%month > 12) THEN time_cur%month = 1 time_cur%year = time_cur%year + 1 ENDIF ENDDO time_cur%day = time_cur%day + d_day IF (time_cur%day < 1) THEN time_cur%month = time_cur%month - 1 IF (time_cur%month < 1) THEN time_cur%month = 12 time_cur%year = time_cur%year - 1 ENDIF time_cur%day = monlen(time_cur%month) + time_cur%day IF (time_cur%month == 2 .AND. MOD(time_cur%year,4) == 0) THEN time_cur%day = time_cur%day + 1 leap = leap - 1 ENDIF ELSE IF (time_cur%day > monlen(time_cur%month)) THEN time_cur%day = time_cur%day - monlen(time_cur%month) IF (time_cur%month == 2 .AND. MOD(time_cur%year,4) == 0) THEN time_cur%day = time_cur%day - 1 leap = leap + 1 IF (time_cur%day == 0) THEN ! special case when you hit 29 Feb on the nose time_cur%day = 29 time_cur%month = time_cur%month - 1 ENDIF ENDIF time_cur%month = time_cur%month + 1 IF (time_cur%month > 12) THEN !this should never be invoked, since Dec has max days (31) time_cur%month = 1 time_cur%year = time_cur%year + 1 ENDIF ENDIF ELSE DO WHILE (d_day < 0) time_cur%month = time_cur%month - 1 IF (time_cur%month == 2 .AND. MOD(time_cur%year,4) == 0) THEN d_day = d_day + 1 leap = leap - 1 ENDIF IF (time_cur%month < 1) THEN time_cur%month = 12 time_cur%year = time_cur%year - 1 ENDIF d_day = d_day + monlen(time_cur%month) ENDDO time_cur%day = time_cur%day + d_day IF (time_cur%month == 2 .AND. MOD(time_cur%year,4) == 0 .AND. time_cur%day == 29) THEN ELSE IF (time_cur%day > monlen(time_cur%month)) THEN time_cur%day = time_cur%day - monlen(time_cur%month) IF (time_cur%month == 2 .AND. MOD(time_cur%year,4) == 0) THEN time_cur%day = time_cur%day - 1 leap = leap + 1 ENDIF time_cur%month = time_cur%month + 1 IF (time_cur%month > 12) THEN time_cur%month = 1 time_cur%year = time_cur%year + 1 ENDIF ENDIF ENDIF ENDIF ELSE WRITE(6,*) ' SUBROUTINE timestamp : unknown flag = ',flag STOP 'timestamp' ENDIF !*** Do the julian day julian = 0 IF (time_cur%month > 1) THEN DO i= 1,time_cur%month-1 julian= julian + monlen(i) ENDDO ENDIF julian = julian + time_cur%day IF (time_cur%month > 2 .AND. MOD(time_cur%year,4) == 0) julian = julian + 1 ! WRITE(6,1000)iyr,mth,day,nhr,min,sec,julian,seconds,leap 1000 FORMAT('Timestamp: ',i5,3i3,2(':',i2.2),' - ',i4,' - ',f14.1,i3) END SUBROUTINE drv_timestamp !----------------------------------------------------------------------------------------