subroutine getcoordinate(nx,ny,lat,lon) 1,4 ! ! !----------------------------------------------------------------------- ! ! get coordinate of arps grid ! ! Author: Ming Hu, CAPS. University of Oklahma. ! ! First written: 04/06/2006. ! !----------------------------------------------------------------------- ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- IMPLICIT NONE ! INCLUDE 'globcst.inc' INCLUDE 'grid.inc' INTEGER :: nx, ny REAL :: x(nx) ! The x-coord. of the u grid point (m) REAL :: y(ny) ! The y-coord. of the v grid point (m) REAL :: xh(nx,ny) ! X-coord. of the terrain (scalar) point. REAL :: yh(nx,ny) ! Y-coord. of the terrain (scalar) point. REAL :: xh1d(0:nx) REAL :: yh1d(0:ny) REAL :: lat(0:nx,0:ny) ! Latitude of each terrain point REAL :: lon(0:nx,0:ny) ! Longitude of each terrain point ! NAMELIST /grid/ dx,dy,ctrlat,ctrlon ! NAMELIST /projection/ mapproj,trulat1,trulat2,trulon,sclfct !----------------------------------------------------------------------- ! ! LOCAL VARIABLES: ! !----------------------------------------------------------------------- ! REAL :: trulat (2) ! Latitude at which the map projection ! is true (degees N) REAL :: latmax,latmin,lonmax,lonmin ! max/min lat/lon of the projected grid INTEGER :: i,j,ih,jh INTEGER :: imin,imax,jmin,jmax,kmin,kmax, istatus REAL :: swx,swy,ctrx,ctry ! !----------------------------------------------------------------------- ! ! TEMPORARY ARRAYS: ! !----------------------------------------------------------------------- ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! IF ( mapproj == 0 ) THEN trulat1 = ctrlat trulat2 = ctrlat trulon = ctrlon END IF !----------------------------------------------------------------------- ! IF( ctrlat > 90 .OR. ctrlat < -90. )then PRINT *,'ctrlat is too large or too small. ', & 'must be set between 0. AND +90. degree north.' STOP END IF ! !----------------------------------------------------------------------- ! ! test to see if trulat is within limits... ! !----------------------------------------------------------------------- ! IF ( trulat1 > 90 .OR. trulat1 < -90 ) then PRINT *,'Input parameter trulat1 is too large or too small.', & 'must be set between 0. AND +90. degree north.' STOP END IF IF ( trulat2 > 90 .OR. trulat2 < -90 ) then PRINT *,'Input parameter trulat2 is too large or too small.', & 'must be set between 0. AND +90. degree north.' STOP END IF ! !----------------------------------------------------------------------- ! ! Set up map projection ! !----------------------------------------------------------------------- ! trulat(1) = trulat1 trulat(2) = trulat2 CALL setmapr( mapproj, 1.0 , trulat , trulon) ! set up parameters for map projection CALL lltoxy( 1,1, ctrlat,ctrlon, ctrx, ctry ) swx = ctrx - (FLOAT(nx-3)/2.) * dx swy = ctry - (FLOAT(ny-3)/2.) * dy CALL setorig( 1, swx, swy) ! !----------------------------------------------------------------------- ! ! Set the horizontal coordinates of the model grid ! !----------------------------------------------------------------------- ! DO i=1,nx x(i) = dx * (i-2) END DO DO j=1,ny y(j) = dy * (j-2) END DO DO j=1,ny DO i=1,nx xh(i,j)=(i-1.5)*dx yh(i,j)=(j-1.5)*dy END DO END DO DO j=0,ny DO i=0,nx xh1d(i)=(i-1.5)*dx yh1d(j)=(j-1.5)*dy END DO END DO !----------------------------------------------------------------------- ! ! Final out the lat/lon of each scalar (terrain data) point. ! !----------------------------------------------------------------------- ! CALL xytoll(nx+1,ny+1,xh1d,yh1d,lat,lon) ! !----------------------------------------------------------------------- ! ! Final out the max/min lat/lon of each scalar (terrain data) point. ! !----------------------------------------------------------------------- write(*,*) ' END of get coordinate' return 980 CONTINUE WRITE(6,'(/1x,a,a)') & 'Error occured when reading namelist input file. Program stopped.' STOP 990 CONTINUE WRITE(6,'(/1x,a,a)') & 'End of file reached when reading namelist input file. ', & 'Program stopped.' STOP END subroutine getCoordinate subroutine getVertcoordinate(nx,ny,nz,zp) 1,4 ! ! !----------------------------------------------------------------------- ! ! get coordinate of arps grid ! ! Author: Ming Hu, CAPS. University of Oklahma. ! ! First written: 04/06/2006. ! !----------------------------------------------------------------------- ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- IMPLICIT NONE INCLUDE 'globcst.inc' INCLUDE 'grid.inc' INTEGER :: nx, ny, nz REAL :: z(nz) ! z-coord. of the computational grid. ! Defined at w-point on the staggered grid. REAL :: zp(nx,ny,nz) ! Physical height coordinate defined at ! w-point of the staggered grid. REAL :: hterain(nx,ny) ! Terrain height. REAL :: zp1d (nz) ! Temporary array REAL :: dzp1d(nz) ! Temporary array REAL :: tem1(nx,ny,nz) ! Temporary array !----------------------------------------------------------------------- ! ! LOCAL VARIABLES: ! !----------------------------------------------------------------------- ! INTEGER :: i,j,k REAL :: zflat1,z1,z2 REAL :: zpmax ! !------------------------------- ! ! ! DO k=1,nz z(k) = zrefsfc + (k-2) * dz END DO ! ! !----------------------------------------------------------------------- ! ! Specify the terrain ! !----------------------------------------------------------------------- ! IF( ternopt == 0 ) THEN ! No terrain, the ground is flat DO j=1,ny-1 DO i=1,nx-1 hterain(i,j) = zrefsfc END DO END DO ELSE IF( ternopt == 1 ) THEN ! Bell-shaped mountain write(*,*) ' Wrong choice of ternopt!' stop 123 ELSE IF( ternopt == 2 ) THEN ! Read from terrain data base ! !----------------------------------------------------------------------- ! ! Read in the terrain data. ! !----------------------------------------------------------------------- ! ! blocking inserted for ordering i/o for message passing ! DO i=0,nprocs-1,readstride ! IF(myproc >= i.AND.myproc <= i+readstride-1)THEN ! ! IF (mp_opt > 0 .AND. readsplit > 0) THEN ! ! CALL readsplittrn( nx,ny,dx,dy, terndta, & ! mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, & ! hterain ) !! ! ! ELSE CALL readtrn( nx,ny,dx,dy, terndta, & mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, & hterain ) ! END IF ! ! END IF ! ! IF (mp_opt > 0) CALL mpbarrier ! END DO END IF ! !----------------------------------------------------------------------- ! ! Set up a stretched vertical grid. ! ! For strhopt=1, function y = a+b*x**3 is used to specify dz as a ! function of k. ! For strhopt=2, function y = c + a*tanh(b*x) is used to specify dz ! as a function of k. ! !----------------------------------------------------------------------- ! IF ( strhopt == 0 ) THEN DO k=1,nz zp1d(k) = z(k) END DO ELSE IF ( strhopt == 1 .OR.strhopt == 2 ) THEN z1 = zrefsfc + MAX(0.0, MIN(dlayer1, z(nz-2)-zrefsfc )) z2 = z1 + MAX(0.0, MIN(dlayer2, z(nz-1)-z1 )) IF( dlayer1 >= (nz-3)*dzmin ) THEN WRITE(6,'(/1x,a,f13.3,/a,f13.3,a,a)') & 'Can not setup a vertical grid with uniform dz=',dzmin, & ' over the depth of ',dlayer1,' please specify a smaller ', & 'value of dlayer1. Program stopped INIGRD.' CALL arpsstop('arpsstop called from INIGRD with ther vertical grid ' & ,1) END IF CALL strhgrd_local(nz,strhopt,zrefsfc,z1,z2,z(nz-1), & dzmin,strhtune, zp1d,dzp1d) ELSE WRITE(6,'(1x,a,i3,a/)') & 'Invalid vertical grid stretching option, strhopt was ',strhopt, & '. Program stopped in INIGRD.' CALL arpsstop('arpsstop called from INIGRD with stretching ' ,1) END IF ! ! !----------------------------------------------------------------------- ! ! Physical height of computational grid defined as ! ! Zp=(z-zrefsfc)*(Zm-hterain)/(Zm-zrefsfc)+hterain for z=<Zm. ! ZP=z for z>Zm ! ! where Zm the height at which the grid levels becomes flat. ! Hm < Zm =< Ztop, hm is the height of mountain and Ztop the height ! of model top. ! !----------------------------------------------------------------------- ! DO k=nz-1,2,-1 IF(zp1d(k) <= zflat) THEN zflat1 = zp1d(k) EXIT END IF END DO zflat1=MAX(MIN(z(nz-1),zflat1),zrefsfc) DO k=2,nz-1 IF(zp1d(k) > zflat1) THEN DO j=1,ny-1 DO i=1,nx-1 zp(i,j,k)=zp1d(k) END DO END DO ELSE DO j=1,ny-1 DO i=1,nx-1 zp(i,j,k)=(zp1d(k)-zrefsfc)*(zflat1-hterain(i,j)) & /(zflat1-zrefsfc)+hterain(i,j) END DO END DO END IF END DO DO j=1,ny-1 DO i=1,nx-1 zp(i,j,2)=hterain(i,j) zp(i,j,1)=2.0*zp(i,j,2)-zp(i,j,3) zp(i,j,nz)=2.0*zp(i,j,nz-1)-zp(i,j,nz-2) END DO END DO END subroutine getVertCoordinate !################################################################## !################################################################## !###### ###### !###### SUBROUTINE STRHGRD ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE strhgrd_local(nz,strhopt,z0,z1,z2,ztop,dzmin,strhtune, z,dzk) 1,53 ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! To construct a vertically stretched grid. ! !----------------------------------------------------------------------- ! ! AUTHOR: Ming Xue ! 10/17/94. ! ! MODIFICATION HISTORY: ! ! 05/11/95 (Jinxing Zong and MX) ! ! A bug fix for the case of nonzero zrefsfc, the reference height of ! the surface. Results not affected for zrefsfc=0 (default value). ! !----------------------------------------------------------------------- ! ! ! INPUT: ! ! nz The vertical dimension of ARPS grid. ! ! ! OUTPUT: ! ! z Array containing the height of veritical grid levels. ! dzk Array containing the grid spacing between vertical levels ! !----------------------------------------------------------------------- IMPLICIT NONE INTEGER :: nz INTEGER :: strhopt REAL :: z0 REAL :: z1 REAL :: z2 REAL :: ztop REAL :: dzmin REAL :: strhtune REAL :: z (nz) REAL :: dzk (nz) REAL :: rnzh,dzm REAL :: a,b,c,hnew,zkmid,dzu INTEGER :: nzh,nzl,k REAL :: dz IF( (z1-z0) == (nz-3)*dzmin.AND.(z1-z0) == (ztop-z0) ) THEN dz = (ztop-z0)/(nz-3) DO k=1,nz-1 dzk(k)= dz END DO DO k=1,nz z(k)=z0 + (k-2) * dz END DO WRITE(6,'(/1x,a,f13.3,/a,f13.3)') & 'Layer 1 depth was as deep as the entire domain. i', & 'A uniform vertical grid is assumed with dz=',dz, & ' over the model depth of ',ztop-z0 RETURN END IF IF(z1 < z0) z1 = z0 IF(z2 > ztop) z2 = ztop nzl = (z1-z0)/dzmin IF( (z1-z0) >= (nz-3)*dzmin ) THEN WRITE(6,'(/1x,a,f13.3,/a,f13.3,a,a)') & 'Can not setup a vertical grid with uniform dz=',dzmin, & ' over the depth of ',z1-z0,' please specify a smaller', & ' value of dlayer1 ' CALL arpsstop('arpsstop called from STRHGRD with stretching ' ,1) END IF IF( z2 >= ztop ) THEN dzm = (ztop-z0-nzl*dzmin)/(nz-3-nzl) ! print*, nzl*dzmin + (nz-3-nzl)*dzm nzh = 0 dzu = 2*dzm - dzmin ELSE a = 2*(nz-3-nzl) b = 2*z0-ztop-z2-(nz-3-3*nzl)*dzmin c = dzmin*(z2-z0-nzl*dzmin) dzm = (-b + SQRT(b*b-4*a*c) )/(2*a) rnzh = (ztop-z2)/(2*dzm-dzmin) nzh = INT(rnzh) hnew = nzl*dzmin + nzh*(2*dzm-dzmin) + & (nz-3-nzl-nzh)*dzm + z0 IF( nzh /= 0 ) THEN dzu = (2*dzm-dzmin) + (ztop-hnew)/nzh ELSE IF( nz-3-nzl-nzh /= 0 ) THEN dzm = dzm + (ztop-hnew)/(nz-3-nzl-nzh) dzu = (2*dzm-dzmin) END IF END IF DO k=1,nzl+1 dzk(k)=dzmin END DO IF( strhopt == 1 ) THEN a = (dzm-dzmin) DO k=nzl+2,nz-2-nzh dzk(k)= dzm+a* & ((2.0*FLOAT(k-nzl-2)/FLOAT(nz-4-nzh-nzl)-1.0) )**3 END DO ELSE zkmid=0.5*FLOAT( nz-nzh+nzl) IF( nzl+2-zkmid == 0.0 ) THEN b = 0.0 ELSE b= strhtune* 2.0/(nzl+2-zkmid) END IF a=(dzmin-dzm)/TANH( strhtune* 2.0) DO k=nzl+2,nz-2-nzh dzk(k)=dzm + a*TANH(b*(FLOAT(k)-zkmid)) END DO END IF DO k=nz-2-nzh+1, nz-2 dzk(k)= dzu END DO dzk(nz-1) = dzk(nz-2) dzk(nz ) = dzk(nz-1) z(2) = z0 DO k=3,nz-1 z(k) = z(k-1)+dzk(k-1) END DO z(1) =z(2)-dzk(1) z(nz-1)=ztop z(nz)=z(nz-1)+dzk(nz-1) RETURN END SUBROUTINE strhgrd_local ! Subroutine OPEN_Mosaic(mosaicfile,NCID) 1 ! IMPLICIT NONE INCLUDE 'netcdf.inc' CHARACTER*120 :: mosaicfile integer :: NCID integer :: status ! STATUS=NF_OPEN(trim(mosaicfile),0,NCID) ! IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) end Subroutine OPEN_Mosaic ! Subroutine CLOSE_Mosaic(NCID) 1 ! IMPLICIT NONE INCLUDE 'netcdf.inc' integer :: NCID integer :: status ! STATUS=NF_CLOSE(NCID) ! IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) end Subroutine CLOSE_Mosaic Subroutine GET_Mosaic_sngl_Mosaic(NCID,mscNlon,mscNlat,mscNlev,mscValue) 1 ! ! Author: Ming Hu, CAPS. University of Oklahma. ! ! First written: 04/06/2006. ! ! IN: ! mscNlon ! mscNlan ! mscNlev ! NCID ! out: ! mscValue ! IMPLICIT NONE INCLUDE 'netcdf.inc' INTEGER :: mscNlon ! number of longitude of mosaic data INTEGER :: mscNlat ! number of latitude of mosaic data INTEGER :: mscNlev ! number of level of mosaic data INTEGER :: NCID, STATUS, MSID INTEGER :: NDIMS PARAMETER (NDIMS=4) ! number of dimensions INTEGER START(NDIMS), COUNT(NDIMS) REAL :: mscValue(mscNlon,mscNlat,1,1) INTEGER :: i,j START(1)=1 START(2)=1 START(3)=mscNlev START(4)=1 COUNT(1)=mscNlon COUNT(2)=mscNlat COUNT(3)=1 COUNT(4)=1 STATUS = NF_INQ_VARID (NCID, 'mrefl_mosaic', MSID) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_GET_VARA_REAL (NCID, MSID, START, COUNT, mscValue) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) end subroutine GET_Mosaic_sngl_Mosaic Subroutine Check_DIM_ATT_Mosaic8(NCID,mscNlon,mscNlat,mscNlev, & 1 lonMin,latMin,lonMax,latMax,dlon,dlat,mosaic_opt) ! ! Author: Ming Hu, CAPS. University of Oklahma. ! ! First written: 04/06/2006. ! ! History: ! ! 05/15/2007 Fanyou Kong ! Add 2D field capability ! ! IN: ! mscNlon ! mscNlan ! mscNlev ! lonMin,latMin,lonMax,latMax,dlon,dlat ! ! out: ! NCID ! IMPLICIT NONE INCLUDE 'netcdf.inc' INTEGER :: mscNlon ! number of longitude of mosaic data INTEGER :: mscNlat ! number of latitude of mosaic data INTEGER :: mscNlev ! number of level of mosaic data REAL :: lonMin,latMin,lonMax,latMax,dlon,dlat INTEGER :: mosaic_opt INTEGER :: NCID, STATUS INTEGER :: LONID, LATID, LEVID INTEGER :: LONLEN, LATLEN, LEVLEN REAL :: lonMinVal REAL :: latMinVal REAL :: lonMaxVal REAL :: latMaxVal REAL :: dlonVal REAL :: dlatVal STATUS = NF_INQ_DIMID(NCID, 'Lon', LONID) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_INQ_DIMID(NCID, 'Lat', LATID) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) IF(mosaic_opt == 1) THEN STATUS = NF_INQ_DIMID(NCID, 'Ht', LEVID) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) END IF STATUS = NF_INQ_DIMLEN(NCID, LONID, LONLEN) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_INQ_DIMLEN(NCID, LATID, LATLEN) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) IF(mosaic_opt == 1) THEN STATUS = NF_INQ_DIMLEN(NCID, LEVID, LEVLEN) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) ELSE LEVLEN=1 END IF if( mscNlon.ne.LONLEN .or. mscNlat.ne.LATLEN .or. mscNlev.ne.LEVLEN) then write(*,*) 'Dimension is inconsistent' write(*,*) 'INPUT:', mscNlon,mscNlat,mscNlev write(*,*) 'Decoding', LONLEN,LATLEN,LEVLEN stop 123 endif STATUS = NF_GET_ATT_REAL (NCID, NF_GLOBAL, 'Longitude', lonMinVal) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_GET_ATT_REAL (NCID, NF_GLOBAL, 'Latitude', latMaxVal) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_GET_ATT_REAL (NCID, NF_GLOBAL, 'LonGridSpacing', dlonVal) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_GET_ATT_REAL (NCID, NF_GLOBAL, 'LatGridSpacing', dlatVal) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) if( abs(lonMin-lonMinVal) > 0.00001 .or. & abs(latMax-latMaxVal) > 0.00001 .or. & abs(dlon-dlonVal) > 0.00001 .or. & abs(dlat-dlatVal) > 0.00001) then write(*,*) ' Attributes are inconsistent. Check' stop 123 endif END SUBROUTINE Check_DIM_ATT_Mosaic8 Subroutine Check_DIM_ATT_Mosaic(NCID,mscNlon,mscNlat,mscNlev, & 1 lonMin,latMin,lonMax,latMax,dlon,dlat) ! ! Author: Ming Hu, CAPS. University of Oklahma. ! ! First written: 04/06/2006. ! ! IN: ! mscNlon ! mscNlan ! mscNlev ! lonMin,latMin,lonMax,latMax,dlon,dlat ! ! out: ! NCID ! IMPLICIT NONE INCLUDE 'netcdf.inc' INTEGER :: mscNlon ! number of longitude of mosaic data INTEGER :: mscNlat ! number of latitude of mosaic data INTEGER :: mscNlev ! number of level of mosaic data REAL :: lonMin,latMin,lonMax,latMax,dlon,dlat INTEGER :: NCID, STATUS INTEGER :: LONID, LATID, LEVID INTEGER :: LONLEN, LATLEN, LEVLEN REAL :: lonMinVal REAL :: latMinVal REAL :: lonMaxVal REAL :: latMaxVal REAL :: dlonVal REAL :: dlatVal STATUS = NF_INQ_DIMID(NCID, 'x', LONID) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_INQ_DIMID(NCID, 'y', LATID) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_INQ_DIMID(NCID, 'mrefl_mosaic_levels', LEVID) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_INQ_DIMLEN(NCID, LONID, LONLEN) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_INQ_DIMLEN(NCID, LATID, LATLEN) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_INQ_DIMLEN(NCID, LEVID, LEVLEN) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) if( mscNlon.ne.LONLEN .or. mscNlat.ne.LATLEN .or. mscNlev.ne.LEVLEN) then write(*,*) 'Dimension is inconsistent' write(*,*) 'INPUT:', mscNlon,mscNlat,mscNlev write(*,*) 'Decoding', LONLEN,LATLEN,LEVLEN stop 123 endif STATUS = NF_GET_ATT_REAL (NCID, NF_GLOBAL, 'xMin', lonMinVal) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_GET_ATT_REAL (NCID, NF_GLOBAL, 'yMin', latMinVal) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_GET_ATT_REAL (NCID, NF_GLOBAL, 'xMax', lonMaxVal) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_GET_ATT_REAL (NCID, NF_GLOBAL, 'yMax', latMaxVal) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_GET_ATT_REAL (NCID, NF_GLOBAL, 'dx', dlonVal) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_GET_ATT_REAL (NCID, NF_GLOBAL, 'dy', dlatVal) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) if( abs(lonMin-lonMinVal) > 0.00001 .or. & abs(lonMax-lonMaxVal) > 0.00001 .or. & abs(latMin-latMinVal) > 0.00001 .or. & abs(latMax-latMaxVal) > 0.00001 .or. & abs(dlon-dlonVal) > 0.00001 .or. & abs(dlat-dlatVal) > 0.00001) then write(*,*) ' Attributes are inconsistent. Check' stop 123 endif END SUBROUTINE Check_DIM_ATT_Mosaic Subroutine GET_DIM_ATT_Mosaic8(mosaicsngle,mscNlon,mscNlat,mscNlev, & 1 lonMin,latMin,lonMax,latMax,dlon,dlat,mosaic_opt) ! ! Author: Ming Hu, CAPS. University of Oklahma. ! ! First written: 04/06/2006. ! ! New verison of Mosaic file that has 8 tiles ! ! History: ! ! 05/15/2007 Fanyou Kong ! Add 2D field capability ! ! IN: ! mosaicPath : path of mosaic file ! mosaicsngle : name of mosaic file ! OUT ! mscNlon ! mscNlan ! mscNlev ! lonMin,latMin,lonMax,latMax,dlon,dlat ! IMPLICIT NONE INCLUDE 'netcdf.inc' CHARACTER*120 mosaicsngle INTEGER :: mscNlon ! number of longitude of mosaic data INTEGER :: mscNlat ! number of latitude of mosaic data INTEGER :: mscNlev ! number of level of mosaic data REAL :: lonMin,latMin,lonMax,latMax,dlon,dlat INTEGER :: mosaic_opt INTEGER :: NCID, STATUS INTEGER :: LONID, LATID, LEVID INTEGER :: LONLEN, LATLEN, LEVLEN REAL :: lonMinVal REAL :: latMinVal REAL :: lonMaxVal REAL :: latMaxVal REAL :: dlonVal REAL :: dlatVal STATUS = NF_OPEN(trim(mosaicsngle), 0, NCID) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_INQ_DIMID(NCID, 'Lon', LONID) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_INQ_DIMID(NCID, 'Lat', LATID) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) IF(mosaic_opt == 1) THEN STATUS = NF_INQ_DIMID(NCID, 'Ht', LEVID) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) END IF STATUS = NF_INQ_DIMLEN(NCID, LONID, LONLEN) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_INQ_DIMLEN(NCID, LATID, LATLEN) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) mscNlon=LONLEN mscNlat=LATLEN IF(mosaic_opt == 1) THEN STATUS = NF_INQ_DIMLEN(NCID, LEVID, LEVLEN) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) mscNlev=LEVLEN ELSE mscNlev=1 END IF STATUS = NF_GET_ATT_REAL (NCID, NF_GLOBAL, 'Longitude', lonMinVal) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_GET_ATT_REAL (NCID, NF_GLOBAL, 'Latitude', latMaxVal) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_GET_ATT_REAL (NCID, NF_GLOBAL, 'LatGridSpacing', dlonVal) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_GET_ATT_REAL (NCID, NF_GLOBAL, 'LatGridSpacing', dlatVal) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) lonMin=lonMinVal latMax=latMaxVal dlon=dlonVal dlat=dlatVal lonMax=lonMinVal+dlon*mscNlon latMin=latMaxVal-dlat*mscNlat STATUS = NF_CLOSE(NCID) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) END SUBROUTINE GET_DIM_ATT_Mosaic8 Subroutine GET_DIM_ATT_Mosaic(mosaicsngle,mscNlon,mscNlat,mscNlev, & 1 lonMin,latMin,lonMax,latMax,dlon,dlat) ! ! Author: Ming Hu, CAPS. University of Oklahma. ! ! First written: 04/06/2006. ! ! IN: ! mosaicPath : path of mosaic file ! mosaicsngle : name of mosaic file ! OUT ! mscNlon ! mscNlan ! mscNlev ! lonMin,latMin,lonMax,latMax,dlon,dlat ! IMPLICIT NONE INCLUDE 'netcdf.inc' CHARACTER*120 mosaicsngle INTEGER :: mscNlon ! number of longitude of mosaic data INTEGER :: mscNlat ! number of latitude of mosaic data INTEGER :: mscNlev ! number of level of mosaic data REAL :: lonMin,latMin,lonMax,latMax,dlon,dlat INTEGER :: NCID, STATUS INTEGER :: LONID, LATID, LEVID INTEGER :: LONLEN, LATLEN, LEVLEN REAL :: lonMinVal REAL :: latMinVal REAL :: lonMaxVal REAL :: latMaxVal REAL :: dlonVal REAL :: dlatVal STATUS = NF_OPEN(trim(mosaicsngle), 0, NCID) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_INQ_DIMID(NCID, 'x', LONID) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_INQ_DIMID(NCID, 'y', LATID) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_INQ_DIMID(NCID, 'mrefl_mosaic_levels', LEVID) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_INQ_DIMLEN(NCID, LONID, LONLEN) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_INQ_DIMLEN(NCID, LATID, LATLEN) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_INQ_DIMLEN(NCID, LEVID, LEVLEN) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) mscNlon=LONLEN mscNlat=LATLEN mscNlev=LEVLEN STATUS = NF_GET_ATT_REAL (NCID, NF_GLOBAL, 'xMin', lonMinVal) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_GET_ATT_REAL (NCID, NF_GLOBAL, 'yMin', latMinVal) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_GET_ATT_REAL (NCID, NF_GLOBAL, 'xMax', lonMaxVal) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_GET_ATT_REAL (NCID, NF_GLOBAL, 'yMax', latMaxVal) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_GET_ATT_REAL (NCID, NF_GLOBAL, 'dx', dlonVal) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_GET_ATT_REAL (NCID, NF_GLOBAL, 'dy', dlatVal) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) lonMin=lonMinVal lonMax=lonMaxVal latMin=latMinVal latMax=latMaxVal dlon=dlonVal dlat=dlatVal STATUS = NF_CLOSE(NCID) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) END SUBROUTINE GET_DIM_ATT_Mosaic ! SUBROUTINE HANDLE_ERR_Mosaic(STATUS) 56 INCLUDE 'netcdf.inc' INTEGER STATUS IF (STATUS .NE. NF_NOERR) THEN PRINT *, NF_STRERROR(STATUS) STOP 'Stopped' ENDIF END SUBROUTINE HANDLE_ERR_Mosaic SUBROUTINE vert_interp_ref(nx,ny,nz,Nmsclvl,ref_mos_3d,ref_arps_3d,zp) 1 ! ! vertical interpolation NSSL mosica reflectiivty into arps grid ! ! by MIng Hu (01/26/2007) ! implicit none INTEGER :: nx,ny,nz INTEGER :: Nmsclvl real :: zp(nx,ny,nz) ! height in w grid real :: ref_arps_3d(nx,ny,nz) ! reflectivity in arps grid real :: ref_mos_3d(nx,ny,Nmsclvl) ! reflectivity in mosaic vertical grid real :: msclvlAll(31) real :: msclvl21(21),msclvl31(31) DATA msclvl21/1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 6, 7, & 8, 9, 10, 11, 12, 13, 14, 15, 16, 17/ DATA msclvl31/0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, & 3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, & 9, 10, 11, 12, 13, 14, 15, 16, 18/ REAL :: heightGSI,upref,downref,wght INTEGER :: ilvl,numref INTEGER :: istat_radar real :: zs(nz) ! height in maa grid INTEGER :: i,j, k2, k ref_arps_3d=-9999.0 numref=0 if (Nmsclvl == 31 ) then DO k=1,Nmsclvl msclvlAll(k)=msclvl31(k)*1000.0 ENDDO elseif( Nmsclvl == 21 ) then msclvlAll=0 DO k=1,Nmsclvl msclvlAll(k)=msclvl21(k)*1000.0 ENDDO else write(*,*) ' Wrong vertical radar mosaic levels' stop 123 endif DO j=1,ny DO i=1,nx DO k2=2,nz-1 zs(k2)=(zp(i,j,k2)+zp(i,j,k2+1))/2 enddo zs(1)=zp(i,j,1) zs(nz)=zp(i,j,nz) DO k2=1,nz heightGSI=zs(k2) ! heightGSI=zp(i,j,k2) if(heightGSI >= msclvlAll(1) .and. heightGSI < msclvlAll(Nmsclvl) ) then do k=1,Nmsclvl-1 if( heightGSI >=msclvlAll(k) .and. heightGSI < msclvlAll(k+1) ) ilvl=k enddo upref=ref_mos_3d(i,j,ilvl+1) downref=ref_mos_3d(i,j,ilvl) if(abs(upref) < 100.0 .and. abs(downref) <100.0 ) then wght=(heightGSI-msclvlAll(ilvl))/(msclvlAll(ilvl+1)-msclvlAll(ilvl)) ref_arps_3d(i,j,k2)=(1-wght)*downref + wght*upref numref=numref+1 else ref_arps_3d(i,j,k2)=-9999.0 endif else ref_arps_3d(i,j,k2)=-9999.0 endif ENDDO !nz ENDDO ENDDO if(numref>0) istat_radar=1 END SUBROUTINE vert_interp_ref ! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE WTRADCOL ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE wtradcol_mosaic(nx,ny,nz,rfname, & 1,2 iyr,imon,iday,ihr,imin,isec, & grdlatc,grdlonc,zpsc,gridref) ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Writes gridded radar data to a file as columns with ! individual lat,lons. ! !----------------------------------------------------------------------- ! ! AUTHOR: Keith Brewster ! 06/22/95 ! ! MODIFICATION HISTORY: ! ! ! 01/26/07 ming HU ! for Mosaic file !----------------------------------------------------------------------- ! ! INPUT: ! dmpfmt file format (1:binary, 2:hdf) ! iradfmt binary format ! hdf4cmpr hdf4 compression level ! rfname radar file name (character*80) ! radid radar id (character*4) ! latrad latitude of radar (degrees N) ! lonrad longitude of radar (degrees E) ! elvrad elevation of radar (m MSL) ! iyr year ! imon month ! iday day ! ihr hour ! imin min ! isec sec ! vcpnum VCP (scan type) number ! isource) source number ! 1= WSR-88D raw ! 2= WSR-88D NIDS ! ! OUTPUT: ! data are written to file ! !----------------------------------------------------------------------- ! IMPLICIT NONE ! INTEGER :: nx,ny,nz ! INTEGER :: dmpfmt INTEGER :: iradfmt INTEGER :: hdf4cmpr CHARACTER (LEN=*) :: rfname CHARACTER (LEN=4) :: radid INTEGER :: iyr,imon,iday,ihr,imin,isec INTEGER :: isource ! REAL :: zs(nz) REAL :: zpsc(nx,ny,nz) REAL :: gridref(nx,ny,nz) ! REAL :: outk(nz) REAL :: outhgt(nz) REAL :: outref(nz) ! !----------------------------------------------------------------------- ! ! Include files ! !----------------------------------------------------------------------- ! INCLUDE 'grid.inc' ! !----------------------------------------------------------------------- ! ! Radar output descriptors ! !----------------------------------------------------------------------- ! !----------------------------------------------------------------------- ! ! Radar output thresholds ! !----------------------------------------------------------------------- ! REAL :: refmin,refmax PARAMETER(refmin=-5.0, refmax=100) REAL :: misval PARAMETER(misval=-999.0) ! !----------------------------------------------------------------------- ! ! Radar output variables ! !----------------------------------------------------------------------- ! REAL :: grdlatc(nx,ny) REAL :: grdlonc(nx,ny) ! !----------------------------------------------------------------------- ! ! Misc local variables ! !----------------------------------------------------------------------- ! INTEGER :: iunit,myr,itime INTEGER :: i,j,k,klev,kk,kntcol,nn INTEGER :: idummy INTEGER :: istat,sd_id INTEGER :: irngmin,irngmax REAL :: gridlat,gridlon,elev,rdummy ! INTEGER(2), allocatable :: itmp(:,:,:) ! Temporary array ! REAL, allocatable :: hmax(:), hmin(:) ! Temporary array CHARACTER (LEN=256) :: filename CHARACTER (LEN=80) :: runname ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! print *, ' inside wtradcol_mosaic' idummy=-999 rdummy=-999. dmpfmt=1 iradfmt=1 hdf4cmpr=1 radid='MOSC' runname='Mosaic2arps' ! IF (dmpfmt > 1 .AND. hdf4cmpr > 3) THEN write(*,*) 'HDF is not available!' stop 123 END IF if( iyr < 100 ) then myr=1900+iyr else myr=iyr endif IF(myr < 1960) myr=myr+100 CALL ctim2abss(myr,imon,iday,ihr,imin,isec,itime) IF(dmpfmt == 1)THEN CALL getunit(iunit) ! !----------------------------------------------------------------------- ! ! Open file for output ! !----------------------------------------------------------------------- ! filename=TRIM(rfname)//'.MscARPS' ! write(*,*) filename OPEN(iunit,FILE=TRIM(filename),STATUS='UNKNOWN',FORM='UNFORMATTED') ! !----------------------------------------------------------------------- ! ! Write radar description variables ! !----------------------------------------------------------------------- ! WRITE(iunit) radid WRITE(iunit) idummy,itime,idummy,isource,idummy, & idummy,idummy,idummy,idummy,idummy ! !----------------------------------------------------------------------- ! ! Write grid description variables ! This should provide enough info to verify that the ! proper grid has been chosen. To recreate the grid, ! icluding elevation information, ! the reading program should get a grid-base-file ! named runname.grdbasfil ! !----------------------------------------------------------------------- ! idummy=0 rdummy=0. WRITE(iunit) runname WRITE(iunit) iradfmt,strhopt,mapproj,idummy,idummy, & idummy,idummy,idummy,idummy,idummy WRITE(iunit) dx,dy,dz,dzmin,ctrlat, & ctrlon,trulat1,trulat2,trulon,sclfct, & rdummy,rdummy,rdummy,rdummy,rdummy WRITE(iunit) idummy,idummy ELSE !HDF4 format ! write(*,*) 'HDF is not availble now!' stop 123 ENDIF ! !----------------------------------------------------------------------- ! ! For each horizontal grid point form a column of remapped ! data containing the non-missing grid points ! !----------------------------------------------------------------------- ! IF(dmpfmt==1)THEN kntcol=0 DO j=1,ny DO i=1,nx DO k=1,nz outk(k)=misval outhgt(k)=misval outref(k)=misval END DO klev=0 DO k=2,nz-1 zs(k)=(zpsc(i,j,k)+zpsc(i,j,k+1))/2 enddo zs(1)=zpsc(i,j,1) zs(nz)=zpsc(i,j,nz) DO k=1,nz-1 IF(gridref(i,j,k)>refmin .AND. gridref(i,j,k)<refmax) THEN klev=klev+1 outk(klev)=FLOAT(k) outhgt(klev)=zs(k) outref(klev)=gridref(i,j,k) END IF END DO ! !----------------------------------------------------------------------- ! ! If there are data in this column, write them to the file. ! !----------------------------------------------------------------------- ! IF(klev > 0) THEN kntcol=kntcol+1 elev=zpsc(i,j,2) WRITE(iunit) i,j,rdummy,rdummy, & grdlatc(i,j),grdlonc(i,j),elev,klev WRITE(iunit) (outk(k),k=1,klev) WRITE(iunit) (outhgt(k),k=1,klev) WRITE(iunit) (outref(k),k=1,klev) END IF END DO END DO ELSE !HDF4 format ! write(*,*) 'HDF is not availble now!' stop 123 ENDIF ! IF(dmpfmt==1)THEN CLOSE(iunit) CALL retunit(iunit) ELSE ! write(*,*) 'HDF is not availble now!' stop 123 ENDIF ! !----------------------------------------------------------------------- ! ! Report on what data were written ! !----------------------------------------------------------------------- ! WRITE(6,'(//a,i4.4,i2.2,i2.2,a1,i2.2,a1,i2.2)') & ' Output statistics for time ', & iyr,imon,iday,' ',ihr,':',imin WRITE(6,'(a,i6,a,/a,i6,a//)') & ' There were ',kntcol,' columns written ', & ' of a total ',(nx*ny),' possible.' ! RETURN END SUBROUTINE wtradcol_mosaic ! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE READADCOL ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE readadcol_Mosaic(nx,ny,nz,rfname,gridref) 2,8 ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! read in gridded radar data to a file as columns with ! individual lat,lons. ! !----------------------------------------------------------------------- ! ! 01/28/06 MIng HU ! Modified to fit the need to write NSSL mosaic reflectiivty ! !----------------------------------------------------------------------- ! ! INPUT: ! dmpfmt file format (1:binary, 2:hdf) ! iradfmt binary format ! hdf4cmpr hdf4 compression level ! rfname radar file name (character*80) ! iyr year ! imon month ! iday day ! ihr hour ! imin min ! isec sec ! isource) source number ! 1= WSR-88D raw ! 2= WSR-88D NIDS ! 3= NSSL mosaic reflectivity ! ! OUTPUT: ! data are written to file ! !----------------------------------------------------------------------- ! IMPLICIT NONE INCLUDE 'grid.inc' ! INTEGER :: nx,ny,nz ! INTEGER :: dmpfmt INTEGER :: iradfmt INTEGER :: hdf4cmpr CHARACTER (LEN=*) :: rfname INTEGER :: iyr,imon,iday,ihr,imin,isec INTEGER :: isource ! REAL :: zs(nz) REAL :: zpsc(nx,ny,nz) REAL :: gridref(nx,ny,nz) ! REAL :: readk(nz) REAL :: readhgt(nz) REAL :: readref(nz) ! !----------------------------------------------------------------------- ! ! Include files ! !----------------------------------------------------------------------- ! ! Radar output descriptors ! !----------------------------------------------------------------------- ! ! INTEGER :: mxradvr,nradvr ! PARAMETER(mxradvr=10,nradvr=6) ! INTEGER :: iradvr(mxradvr) ! DATA iradvr /1,2,3,4,5,6,0,0,0,0/ ! !----------------------------------------------------------------------- ! ! Radar output thresholds ! !----------------------------------------------------------------------- ! REAL :: refmin,refmax,velmin,velmax PARAMETER(refmin=-5.0, refmax=100., & velmin=-200.,velmax=200.) REAL :: misval PARAMETER(misval=-999.0) ! !----------------------------------------------------------------------- ! ! Radar output variables ! !----------------------------------------------------------------------- ! REAL :: grdlatc(nx,ny) REAL :: grdlonc(nx,ny) ! !----------------------------------------------------------------------- ! ! Misc local variables ! !----------------------------------------------------------------------- ! CHARACTER (LEN=4) :: radid CHARACTER (LEN=12) :: runname INTEGER :: iunit,myr,itime INTEGER :: i,j,k,klev,kk,kntcol,nn INTEGER :: idummy INTEGER :: istat,sd_id INTEGER :: ipt REAL :: gridlat,gridlon,elev,rdummy INTEGER(2), allocatable :: itmp(:,:,:) ! Temporary array REAL, allocatable :: hmax(:), hmin(:) ! Temporary array CHARACTER (LEN=256) :: filename ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! print *, ' inside readadcol_mosaic' idummy=-999 rdummy=-999. dmpfmt=1 iradfmt=1 hdf4cmpr=1 ! IF (dmpfmt > 1 .AND. hdf4cmpr > 3) THEN write(*,*) 'HDF is not available!' stop 123 END IF IF(dmpfmt == 1)THEN CALL getunit(iunit) ! !----------------------------------------------------------------------- ! ! Open file for output ! !----------------------------------------------------------------------- ! filename=TRIM(rfname)//'.MscARPS' write(*,*) filename OPEN(iunit,FILE=TRIM(filename),STATUS='UNKNOWN',FORM='UNFORMATTED') ! !----------------------------------------------------------------------- ! ! Write radar description variables ! !----------------------------------------------------------------------- ! read(iunit) radid read(iunit) idummy,itime,idummy,isource,idummy, & idummy,idummy,idummy,idummy,idummy ! CALL abss2ctim(itime,myr,imon,iday,ihr,imin,isec) ! write(*,*) 'The time of the data is: ',myr,imon,iday,ihr,imin,isec !----------------------------------------------------------------------- ! ! Write grid description variables ! This should provide enough info to verify that the ! proper grid has been chosen. To recreate the grid, ! icluding elevation information, ! the reading program should get a grid-base-file ! named runname.grdbasfil ! !----------------------------------------------------------------------- ! idummy=0 rdummy=0. read(iunit) runname read(iunit) iradfmt,strhopt,mapproj,idummy,idummy, & idummy,idummy,idummy,idummy,idummy read(iunit) dx,dy,dz,dzmin,ctrlat, & ctrlon,trulat1,trulat2,trulon,sclfct, & rdummy,rdummy,rdummy,rdummy,rdummy read(iunit) idummy,idummy ELSE !HDF4 format ! write(*,*) 'HDF is not availble now!' stop 123 ENDIF ! !----------------------------------------------------------------------- ! ! For each horizontal grid point form a column of remapped ! data containing the non-missing grid points ! !----------------------------------------------------------------------- ! IF(dmpfmt==1)THEN DO ipt=1,(nx*ny) read(iunit,END=51) i,j,rdummy,rdummy, & gridlat,gridlon,elev,klev read(iunit,END=51) (readk(k),k=1,klev) read(iunit,END=51) (readhgt(k),k=1,klev) read(iunit,END=51) (readref(k),k=1,klev) IF(i <= nx.AND.i >= 1 .AND. j <= ny.AND.j >= 1) THEN DO kk=1,klev k=nint(readk(kk)) IF(k <= nz.AND.k >= 1) THEN gridref(i,j,k)=readref(kk) END IF ! 1 < k < nz END DO ! kk = 1, klev END IF ! 1 < i < nx & 1 < j < ny END DO ! ipt = 1, nx*ny 51 continue ipt=ipt-1 WRITE(6,'(a,i6,a)') ' End of file reached after reading', & ipt,' columns' ELSE !HDF4 format ! write(*,*) 'HDF is not availble now!' stop 123 ENDIF ! RETURN END SUBROUTINE readadcol_Mosaic Subroutine GET_Mosaic_cref_Mosaic(NCID,mscNlon,mscNlat,mscValue) ! ! Read in 2D field CREF ! ! Author: Fanyou Kong, CAPS. University of Oklahma. ! First written: 05/15/2007. ! ! IN: ! mscNlon ! mscNlan ! NCID ! out: ! mscValue ! IMPLICIT NONE INCLUDE 'netcdf.inc' INTEGER :: mscNlon ! number of longitude of mosaic data INTEGER :: mscNlat ! number of latitude of mosaic data INTEGER :: NCID, STATUS, MSID INTEGER :: NDIMS PARAMETER (NDIMS=2) ! number of dimensions INTEGER START(NDIMS), COUNT(NDIMS) REAL :: mscValue(mscNlon,mscNlat) INTEGER :: i,j START(1)=1 START(2)=1 COUNT(1)=mscNlon COUNT(2)=mscNlat STATUS = NF_INQ_VARID (NCID, 'cref', MSID) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_GET_VARA_REAL (NCID, MSID, START, COUNT, mscValue) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) end subroutine GET_Mosaic_cref_Mosaic Subroutine GET_Mosaic_hsr_Mosaic(NCID,mscNlon,mscNlat,mscValue) 1 ! ! Read in 2D field 1h_rad_hsr ! ! Author: Fanyou Kong, CAPS. University of Oklahma. ! First written: 05/15/2007. ! ! IN: ! mscNlon ! mscNlan ! NCID ! out: ! mscValue ! IMPLICIT NONE INCLUDE 'netcdf.inc' INTEGER :: mscNlon ! number of longitude of mosaic data INTEGER :: mscNlat ! number of latitude of mosaic data INTEGER :: NCID, STATUS, MSID INTEGER :: NDIMS PARAMETER (NDIMS=2) ! number of dimensions INTEGER START(NDIMS), COUNT(NDIMS) REAL :: mscValue(mscNlon,mscNlat) INTEGER :: i,j START(1)=1 START(2)=1 COUNT(1)=mscNlon COUNT(2)=mscNlat STATUS = NF_INQ_VARID (NCID, '1h_rad_hsr', MSID) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) STATUS = NF_GET_VARA_REAL (NCID, MSID, START, COUNT, mscValue) IF (STATUS .NE. NF_NOERR) CALL HANDLE_ERR_Mosaic(STATUS) end subroutine GET_Mosaic_hsr_Mosaic