PROGRAM sstuvhflux !$$$ main program documentation block ! ! main program: sstuvhflux ! PRGMMR: RICHARD YABLONSKY, URI/GSO, 2014-12-02 ! LANG: FORTRAN 90 ! Abstract: Read 1.5-hourly SST from MPIPOM-TC output (NetCDF) and ! calculate average, minimum, and maximum SST within 60, 100, 150, and ! 200 km radius around the storm center, as well as the SST difference ! from both phase 1 and the end of phase 2 (start of coupled forecast). ! ! ! program history log: ! USE netcdf IMPLICIT NONE INTERFACE SUBROUTINE check(index) USE netcdf IMPLICIT NONE INTEGER, INTENT(IN) :: index END SUBROUTINE END INTERFACE INTEGER, PARAMETER :: im = 869, jm = 449 INTEGER :: i, j, k, fid, n, nf, nprinto2 INTEGER :: ilonid, ilatid, isstid INTEGER :: wrad, fhr1, fhr2 REAL :: d2r, rearth, tavr, counter REAL :: tavrd1, tavrd2 REAL :: avgsst, avgsstd1, avgsstd2 REAL :: minsst, minsstd1, minsstd2 REAL :: maxsst, maxsstd1, maxsstd2 REAL :: xcen, ycen, fhr, dx, dy, r REAL :: xcen1, ycen1, xcen2, ycen2 REAL, DIMENSION(4) :: dist REAL, DIMENSION(im,jm) :: alon, alat, sst, sst0ph1, sst0cpl CHARACTER(LEN=18) :: sstinfile CHARACTER(LEN=12) :: sstoutfile CHARACTER(LEN=29) :: junk1 CHARACTER(LEN=1) :: ns, ew, junk2, junk3 CHARACTER(LEN=16) :: junk4 CHARACTER(LEN=190) :: junk5 dist(1) = 60000. dist(2) = 100000. dist(3) = 150000. dist(4) = 200000. DO n = 1,4 ! Open file for writing table WRITE(sstoutfile,'(''sst'',i3.3,''km.dat'')') int(dist(n)/1000.) OPEN(9,FILE=sstoutfile) ! Write header information to the table WRITE(9,100) "%FHOUR LON LAT AVGSST MINSST MAXSST & AVGSSTD1 MINSSTD1 MAXSSTD1 AVGSSTD2 MINSSTD2 MAXSSTD2" 100 FORMAT(a96) ! Read SST data from phase 1 hour 0 MPIPOM-TC output file (NetCDF) WRITE(sstinfile,'(''../pom/output/OCEAN/PHASE1/sstuvhflux.'', & i4.4,''.nc'')') 0 CALL check(NF90_OPEN(sstinfile,NF90_NOWRITE,fid)) CALL check(NF90_INQ_VARID(fid,"sst",isstid)) CALL check(NF90_GET_VAR(fid,isstid,sst0ph1)) CALL check(NF90_CLOSE(fid)) ! Read SST data from coupled hour 0 MPIPOM-TC output file (NetCDF) WRITE(sstinfile,'(''sstuvhflux.'',i4.4,''.nc'')') 0 CALL check(NF90_OPEN(sstinfile,NF90_NOWRITE,fid)) CALL check(NF90_INQ_VARID(fid,"sst",isstid)) CALL check(NF90_GET_VAR(fid,isstid,sst0cpl)) CALL check(NF90_CLOSE(fid)) OPEN(10) DO nf = 1,85 ! Read SST data from 1.5-hourly MPIPOM-TC output file (NetCDF) nprinto2 = nf-1 WRITE(sstinfile,'(''sstuvhflux.'',i4.4,''.nc'')') nprinto2 CALL check(NF90_OPEN(sstinfile,NF90_NOWRITE,fid)) CALL check(NF90_INQ_VARID(fid,"sst",isstid)) CALL check(NF90_GET_VAR(fid,isstid,sst)) CALL check(NF90_INQ_VARID(fid,"east_e",ilonid)) CALL check(NF90_GET_VAR(fid,ilonid,alon)) CALL check(NF90_INQ_VARID(fid,"north_e",ilatid)) CALL check(NF90_GET_VAR(fid,ilatid,alat)) CALL check(NF90_CLOSE(fid)) ! Read storm center lat/lon from 3-hourly track file (fort.10) 101 FORMAT(a29,i4,a1,f4.0,a1,a1,f5.0,a1,a16,i4,a190) IF(nprinto2.eq.0) THEN wrad = 0 fhr1 = 1 DO WHILE(wrad.ne.34.or.MOD(fhr1,3).ne.0) READ(10,101) junk1,fhr1,junk2,ycen1,ns,junk3,xcen1,ew,junk4,wrad,junk5 xcen1 = xcen1/10. ycen1 = ycen1/10. IF(ew.eq.'W') xcen1 = -xcen1 IF(ns.eq.'S') ycen1 = -ycen1 END DO wrad = 0 fhr2 = 1 DO WHILE(wrad.ne.34.or.MOD(fhr2,3).ne.0) READ(10,101) junk1,fhr2,junk2,ycen2,ns,junk3,xcen2,ew,junk4,wrad,junk5 xcen2 = xcen2/10. ycen2 = ycen2/10. IF(ew.eq.'W') xcen2 = -xcen2 IF(ns.eq.'S') ycen2 = -ycen2 END DO fhr = fhr1 xcen = xcen1 ycen = ycen1 ELSE IF(MOD(nprinto2,2).eq.0) THEN fhr1 = fhr2 xcen1 = xcen2 ycen1 = ycen2 wrad = 0 fhr2 = 1 DO WHILE((wrad.ne.34.or.MOD(fhr2,3).ne.0).and.fhr1.lt.126) READ(10,101) junk1,fhr2,junk2,ycen2,ns,junk3,xcen2,ew,junk4,wrad,junk5 xcen2 = xcen2/10. ycen2 = ycen2/10. IF(ew.eq.'W') xcen2 = -xcen2 IF(ns.eq.'S') ycen2 = -ycen2 END DO fhr = fhr1 xcen = xcen1 ycen = ycen1 ELSE fhr = 0.5*(fhr1+fhr2) xcen = 0.5*(xcen1+xcen2) ycen = 0.5*(ycen1+ycen2) END IF ! Write data to the table rearth = 6371000. d2r = 3.141592654/180. tavr = 0. tavrd1 = 0. tavrd2 = 0. maxsst = 0. maxsstd1 = 0. maxsstd2 = 0. minsst = 100. minsstd1 = 0. minsstd2 = 0. counter = 0. DO j = 1,jm DO i = 1,im dx = rearth*(alon(i,j)-xcen)*d2r*cos(ycen*d2r) dy = rearth*(alat(i,j)-ycen)*d2r r = sqrt(dx**2.+dy**2.) IF(r.lt.dist(n).and.sst(i,j).ne.0.) THEN tavr = tavr+sst(i,j) tavrd1 = tavrd1+sst(i,j)-sst0ph1(i,j) tavrd2 = tavrd2+sst(i,j)-sst0cpl(i,j) maxsst = max(sst(i,j),maxsst) maxsstd1 = max(sst(i,j)-sst0ph1(i,j),maxsstd1) maxsstd2 = max(sst(i,j)-sst0cpl(i,j),maxsstd2) minsst = min(sst(i,j),minsst) minsstd1 = min(sst(i,j)-sst0ph1(i,j),minsstd1) minsstd2 = min(sst(i,j)-sst0cpl(i,j),minsstd2) counter = counter+1. END IF END DO END DO IF(counter.gt.0.) THEN avgsst = tavr/counter avgsstd1 = tavrd1/counter avgsstd2 = tavrd2/counter ELSE avgsst = 0. avgsstd1 = 0. avgsstd2 = 0. END IF WRITE(9,102) fhr,xcen,ycen,avgsst,minsst,maxsst, & avgsstd1,minsstd1,maxsstd1,avgsstd2,minsstd2,maxsstd2 102 FORMAT(f6.1,f8.2,f7.2,3f7.2,6f9.2) END DO CLOSE(10) END DO END PROGRAM SUBROUTINE check(index) USE netcdf IMPLICIT NONE INTEGER, INTENT(IN) :: index IF( index /= nf90_noerr) THEN PRINT*, nf90_strerror(index) STOP ENDIF END SUBROUTINE