PROGRAM SHALOW * In this version, shallow4.f, initial and calculated values * of U, V, and P are written to a netCDF file * for later use in visualizing the results. The netCDF data * management library is freely available from * http://www.unidata.ucar.edu/software/netcdf * This code is still serial but has been brought up to modern * Fortran constructs and uses portable intrinsic Fortran 90 timing routines * This can be compiled on the IBM SP using: * xlf90 -qmaxmem=-1 -g -o shallow4 -qfixed=132 -qsclk=micro \ * -I/usr/local/include shallow4.f -L/usr/local/lib32/r4i4 -l netcdf * where the -L and -I point to local installation of netCDF include 'netcdf.inc' C C BENCHMARK WEATHER PREDICTION PROGRAM FOR COMPARING THE C PREFORMANCE OF CURRENT SUPERCOMPUTERS. THE MODEL IS C BASED OF THE PAPER - THE DYNAMICS OF FINITE-DIFFERENCE C MODELS OF THE SHALLOW-WATER EQUATIONS, BY ROBERT SADOURNY C J. ATM. SCIENCES, VOL 32, NO 4, APRIL 1975. C C CODE BY PAUL N. SWARZTRAUBER, NATIONAL CENTER FOR C ATMOSPHERIC RESEARCH, BOULDER, CO, OCTOBER 1984. C Modified by Juliana Rew, NCAR, January 2006 integer mlen,nlen parameter (m_len = 65) parameter (n_len = 65) parameter (m = 64) parameter (n = 64) DIMENSION U(m_len,n_len),V(m_len,n_len),P(m_len,n_len),UNEW(m_len,n_len),VNEW(m_len,n_len), 1 PNEW(m_len,n_len),UOLD(m_len,n_len),VOLD(m_len,n_len),POLD(m_len,n_len), 2 CU(m_len,n_len),CV(m_len,n_len),Z(m_len,n_len),H(m_len,n_len),PSI(m_len,n_len) * REAL MFS100,MFS200,MFS300 double precision MFS100,MFS200,MFS300 double precision T100,T200,T300,TSTART,CTIME,TCYC integer c1,c2,r,max integer ncid,p_id,u_id,v_id integer istart(3) integer icount(3) * prepare netCDF file to receive model output data call netcdf_setup(m_len,n_len,ncid,p_id,u_id,v_id,istart,icount) C C NOTE BELOW THAT TWO DELTA T (TDT) IS SET TO DT ON THE FIRST C CYCLE AFTER WHICH IT IS RESET TO DT+DT. C DT = 90. TDT = DT C DX = 1.E5 DY = 1.E5 A = 1.E6 ALPHA = .001 ITMAX = 400 MPRINT = 100 * M and N could be set as constants above to make it easier * to change the size of the problem above * M = 64 * N = 64 MP1 = M+1 NP1 = N+1 EL = FLOAT(N)*DX PI = 4.*ATAN(1.) TPI = PI+PI DI = TPI/FLOAT(M) DJ = TPI/FLOAT(N) PCF = PI*PI*A*A/(EL*EL) C C INITIAL VALUES OF THE STREAM FUNCTION AND P C DO J=1,NP1 DO I=1,MP1 PSI(I,J) = A*SIN((FLOAT(I)-.5)*DI)*SIN((FLOAT(J)-.5)*DJ) * PSI(I,J) = 0 P(I,J) = PCF*(COS(2.*FLOAT(I-1)*DI) 1 +COS(2.*FLOAT(J-1)*DJ))+50000. END DO 50 END DO C C INITIALIZE VELOCITIES C DO J=1,N DO I=1,M U(I+1,J) = -(PSI(I+1,J+1)-PSI(I+1,J))/DY V(I,J+1) = (PSI(I+1,J+1)-PSI(I,J+1))/DX END DO 60 END DO C C PERIODIC CONTINUATION C DO J=1,N U(1,J) = U(M+1,J) V(M+1,J+1) = V(1,J+1) 70 END DO DO I=1,M U(I+1,N+1) = U(I+1,1) V(I,1) = V(I,N+1) 75 END DO U(1,N+1) = U(M+1,1) V(M+1,1) = V(1,N+1) DO J=1,NP1 DO I=1,MP1 UOLD(I,J) = U(I,J) VOLD(I,J) = V(I,J) POLD(I,J) = P(I,J) END DO 86 END DO C C PRINT INITIAL VALUES C WRITE(6,390) N,M,DX,DY,DT,ALPHA 390 FORMAT("1NUMBER OF POINTS IN THE X DIRECTION",I8,/ 1 " NUMBER OF POINTS IN THE Y DIRECTION",I8,/ 2 " GRID SPACING IN THE X DIRECTION ",F8.0,/ 3 " GRID SPACING IN THE Y DIRECTION ",F8.0,/ 4 " TIME STEP ",F8.0,/ 5 " TIME FILTER PARAMETER ",F8.3) MNMIN = MIN0(M,N) WRITE(6,391) (P(I,I),I=1,MNMIN) 391 FORMAT(/' INITIAL DIAGONAL ELEMENTS OF P ' //,(8E15.6)) WRITE(6,392) (U(I,I),I=1,MNMIN) 392 FORMAT(/' INITIAL DIAGONAL ELEMENTS OF U ' //,(8E15.6)) WRITE(6,393) (V(I,I),I=1,MNMIN) 393 FORMAT(/' INITIAL DIAGONAL ELEMENTS OF V ' //,(8E15.6)) * write intial values of p, u, and v into a netCDF file call my_ncwrite(ncid,p_id,istart,icount,p,m_len,n_len) call my_ncwrite(ncid,u_id,istart,icount,u,m_len,n_len) call my_ncwrite(ncid,v_id,istart,icount,v,m_len,n_len) * TSTART = SECOND(DUM) call system_clock (count=c1, count_rate=r, count_max=max) TSTART = c1 T300 = 1. TIME = 0. NCYCLE = 0 * 90 NCYCLE = NCYCLE + 1 DO ncycle=1,itmax C C COMPUTE CAPITAL U, CAPITAL V, Z AND H C FSDX = 4./DX FSDY = 4./DY * T100 = SECOND(DUM) call system_clock(count=c1, count_rate=r,count_max=max) T100 = c1 DO J=1,N DO I=1,M CU(I+1,J) = .5*(P(I+1,J)+P(I,J))*U(I+1,J) CV(I,J+1) = .5*(P(I,J+1)+P(I,J))*V(I,J+1) Z(I+1,J+1) =(FSDX*(V(I+1,J+1)-V(I,J+1))-FSDY*(U(I+1,J+1) 1 -U(I+1,J)))/(P(I,J)+P(I+1,J)+P(I+1,J+1)+P(I,J+1)) H(I,J) = P(I,J)+.25*(U(I+1,J)*U(I+1,J)+U(I,J)*U(I,J) 1 +V(I,J+1)*V(I,J+1)+V(I,J)*V(I,J)) END DO 100 END DO * T100 = SECOND(DUM)-T100 call system_clock(count=c2,count_rate=r,count_max=max) T100 = dble(c2-T100)/dble(r) C C PERIODIC CONTINUATION C DO J=1,N CU(1,J) = CU(M+1,J) CV(M+1,J+1) = CV(1,J+1) Z(1,J+1) = Z(M+1,J+1) H(M+1,J) = H(1,J) 110 END DO DO I=1,M CU(I+1,N+1) = CU(I+1,1) CV(I,1) = CV(I,N+1) Z(I+1,1) = Z(I+1,N+1) H(I,N+1) = H(I,1) 115 END DO CU(1,N+1) = CU(M+1,1) CV(M+1,1) = CV(1,N+1) Z(1,1) = Z(M+1,N+1) H(M+1,N+1) = H(1,1) C C COMPUTE NEW VALUES U,V AND P C TDTS8 = TDT/8. TDTSDX = TDT/DX TDTSDY = TDT/DY * T200 = SECOND(DUM) call system_clock(count=c1, count_rate=r, count_max=max) T200 = c1 DO J=1,N DO I=1,M UNEW(I+1,J) = UOLD(I+1,J)+ 1 TDTS8*(Z(I+1,J+1)+Z(I+1,J))*(CV(I+1,J+1)+CV(I,J+1)+CV(I,J) 2 +CV(I+1,J))-TDTSDX*(H(I+1,J)-H(I,J)) VNEW(I,J+1) = VOLD(I,J+1)-TDTS8*(Z(I+1,J+1)+Z(I,J+1)) 1 *(CU(I+1,J+1)+CU(I,J+1)+CU(I,J)+CU(I+1,J)) 2 -TDTSDY*(H(I,J+1)-H(I,J)) PNEW(I,J) = POLD(I,J)-TDTSDX*(CU(I+1,J)-CU(I,J)) 1 -TDTSDY*(CV(I,J+1)-CV(I,J)) END DO 200 END DO * T200 = SECOND(DUM)-T200 call system_clock(count=c2, count_rate=r, count_max=max) T200 = dble(c2 -T200)/dble(r) C C PERIODIC CONTINUATION C DO J=1,N UNEW(1,J) = UNEW(M+1,J) VNEW(M+1,J+1) = VNEW(1,J+1) PNEW(M+1,J) = PNEW(1,J) 210 END DO DO I=1,M UNEW(I+1,N+1) = UNEW(I+1,1) VNEW(I,1) = VNEW(I,N+1) PNEW(I,N+1) = PNEW(I,1) 215 END DO UNEW(1,N+1) = UNEW(M+1,1) VNEW(M+1,1) = VNEW(1,N+1) PNEW(M+1,N+1) = PNEW(1,1) * jr IF(NCYCLE .GT. ITMAX) CALL EXIT TIME = TIME + DT * jr IF(MOD(NCYCLE,MPRINT) .NE. 0) GO TO 370 IF(MOD(NCYCLE,MPRINT) .EQ. 0) then * jr 370 IF(NCYCLE .LE. 1) GO TO 310 PTIME = TIME/3600. WRITE(6,350) NCYCLE,PTIME 350 FORMAT(//' CYCLE NUMBER',I5,' MODEL TIME IN HOURS', F6.2) WRITE(6,355) (PNEW(I,I),I=1,MNMIN) 355 FORMAT(/' DIAGONAL ELEMENTS OF P ' //(8E15.6)) WRITE(6,360) (UNEW(I,I),I=1,MNMIN) 360 FORMAT(/' DIAGONAL ELEMENTS OF U ' //(8E15.6)) WRITE(6,365) (VNEW(I,I),I=1,MNMIN) 365 FORMAT(/' DIAGONAL ELEMENTS OF V ' //(8E15.6)) * jr added MFS310--don't know what actual mult factor should be * jr changed divide by 1 million to divide by 100K since system_clock * jr resolution is millisec rather than cpu_time's 10 millisec MFS310 = 24.*M*N/T310/1.D5 * MFS100 = 24.*M*N/T100/1.E6 MFS100 = 24.*M*N/T100/1.D5 * MFS200 = 26.*M*N/T200/1.E6 MFS200 = 26.*M*N/T200/1.D5 * MFS300 = 15.*M*N/T300/1.E6 MFS300 = 15.*M*N/T300/1.D5 * CTIME = SECOND(DUM)-TSTART call system_clock(count=c2, count_rate=r,count_max=max) CTIME = dble(c2 - TSTART)/dble(r) TCYC = CTIME/FLOAT(NCYCLE) WRITE(6,375) NCYCLE,CTIME,TCYC,T310,MFS310,T200,MFS200,T300,MFS300 * 375 FORMAT(' CYCLE NUMBER',I5,' TOTAL COMPUTER TIME', E15.6, * 1 ' TIME PER CYCLE', E15.6, / * 2 ' TIME AND MEGAFLOPS FOR LOOP 310 ', E15.6,F6.1/ * 3 ' TIME AND MEGAFLOPS FOR LOOP 200 ', E15.6,F6.1/ * 4 ' TIME AND MEGAFLOPS FOR LOOP 300 ', E15.6,F6.1/ ) 375 FORMAT(' CYCLE NUMBER',I5,' TOTAL COMPUTER TIME', D15.6, 1 ' TIME PER CYCLE', D15.6, / 2 ' TIME AND MEGAFLOPS FOR LOOP 310 ', D15.6,2x,D6.1/ 3 ' TIME AND MEGAFLOPS FOR LOOP 200 ', D15.6,2x,D6.1/ 4 ' TIME AND MEGAFLOPS FOR LOOP 300 ', D15.6,2x,D6.1/ ) endif * append calculated values of p, u, and v to netCDF file istart(3) = ncycle + 1 * shape of record to be written (one ncycle at a time) call my_ncwrite(ncid,p_id,istart,icount,p,m_len,n_len) call my_ncwrite(ncid,u_id,istart,icount,u,m_len,n_len) call my_ncwrite(ncid,v_id,istart,icount,v,m_len,n_len) C C TIME SMOOTHING AND UPDATE FOR NEXT CYCLE C * 370 IF(NCYCLE .LE. 1) GO TO 310 370 IF(NCYCLE .GT. 1) then * T300 = SECOND(DUM) call system_clock(count=c1,count_rate=r,count_max=max) T300 = c1 DO J=1,N DO I=1,M UOLD(I,J) = U(I,J)+ALPHA*(UNEW(I,J)-2.*U(I,J)+UOLD(I,J)) VOLD(I,J) = V(I,J)+ALPHA*(VNEW(I,J)-2.*V(I,J)+VOLD(I,J)) POLD(I,J) = P(I,J)+ALPHA*(PNEW(I,J)-2.*P(I,J)+POLD(I,J)) U(I,J) = UNEW(I,J) V(I,J) = VNEW(I,J) P(I,J) = PNEW(I,J) END DO 300 END DO * T300 = SECOND(DUM)-T300 call system_clock(count=c2,count_rate=r, count_max=max) T300 = dble(c2 - T300)/dble(r) C C PERIODIC CONTINUATION C DO J=1,N UOLD(M+1,J) = UOLD(1,J) VOLD(M+1,J) = VOLD(1,J) POLD(M+1,J) = POLD(1,J) U(M+1,J) = U(1,J) V(M+1,J) = V(1,J) P(M+1,J) = P(1,J) 320 END DO DO I=1,M UOLD(I,N+1) = UOLD(I,1) VOLD(I,N+1) = VOLD(I,1) POLD(I,N+1) = POLD(I,1) U(I,N+1) = U(I,1) V(I,N+1) = V(I,1) P(I,N+1) = P(I,1) 325 END DO UOLD(M+1,N+1) = UOLD(1,1) VOLD(M+1,N+1) = VOLD(1,1) POLD(M+1,N+1) = POLD(1,1) U(M+1,N+1) = U(1,1) V(M+1,N+1) = V(1,1) P(M+1,N+1) = P(1,1) else 310 TDT = TDT+TDT call system_clock(count=c1, count_rate=r,count_max=max) T310 = c1 DO J=1,NP1 DO I=1,MP1 UOLD(I,J) = U(I,J) VOLD(I,J) = V(I,J) POLD(I,J) = P(I,J) U(I,J) = UNEW(I,J) V(I,J) = VNEW(I,J) P(I,J) = PNEW(I,J) END DO 400 END DO call system_clock(count=c2, count_rate=r, count_max=max) T310 = dble(c2 - T310)/dble(r) endif * end big loop 90 End do * close the netCDF file iret = nf_close(ncid) call check_err(iret) END subroutine check_err(iret) integer iret include 'netcdf.inc' if(iret .ne. NF_NOERR) then print *, nf_strerror(iret) stop endif end subroutine netcdf_setup(m_len,n_len,ncid,p_id,u_id,v_id,istart,icount) * Input args: m_len, n_len * Output args: ncid, p_id,u_id,v_id,istart,icount) integer m_len,n_len * declarations for netCDF library include 'netcdf.inc' * error status return integer iret * netCDF id integer ncid * dimension ids integer m_dim integer n_dim integer time_dim * variable ids integer p_id integer u_id integer v_id * rank (number of dimensions) for each variable integer p_rank, u_rank, v_rank parameter (p_rank = 3) parameter (u_rank = 3) parameter (v_rank = 3) * variable shapes integer p_dims(p_rank) integer u_dims(u_rank) integer v_dims(v_rank) integer istart(p_rank) integer icount(p_rank) * enter define mode iret = nf_create('shallowdat.nc', NF_CLOBBER, + ncid) call check_err(iret) * define dimensions iret = nf_def_dim(ncid, 'm', m_len, m_dim) call check_err(iret) iret = nf_def_dim(ncid, 'n', n_len, n_dim) call check_err(iret) * time is an unlimited dimension so that any number of * records can be added iret = nf_def_dim(ncid, 'time', NF_UNLIMITED, time_dim) call check_err(iret) * define variables p_dims(1) = m_dim p_dims(2) = n_dim p_dims(3) = time_dim iret = nf_def_var(ncid, 'p', NF_REAL, p_rank, p_dims, p_id) call check_err(iret) u_dims(1) = m_dim u_dims(2) = n_dim u_dims(3) = time_dim iret = nf_def_var(ncid, 'u', NF_REAL, u_rank, u_dims, u_id) call check_err(iret) v_dims(1) = m_dim v_dims(2) = n_dim v_dims(3) = time_dim iret = nf_def_var(ncid, 'v', NF_REAL, v_rank, v_dims, v_id) call check_err(iret) * start netCDF write at the first (1,1,1) position of the array istart(1) = 1 istart(2) = 1 istart(3) = 1 * shape of record to be written (one ncycle at a time) icount(1) = m_len icount(2) = n_len icount(3) = 1 * leave define mode iret = nf_enddef(ncid) call check_err(iret) * end of netCDF definitions end subroutine netcdf_setup subroutine my_ncwrite(id,varid,istart,icount,var,m_len,n_len) * Input args: id, varid,istart,icount,var,m_len_n_len * Write a whole array out to the netCDF file integer id,varid,iret integer icount(3) integer istart(3) real var (m_len,n_len) iret = nf_put_vara_real(id,varid,istart,icount,var) call check_err(iret) end subroutine my_ncwrite