! subroutine nmczday(lu,mssfile,mxrec,iyd0,iyd1,iyd,z,glon,nlon, | nlat,npr,ier,iprint) implicit none c c (BF: modified 6/9 to include 3 pressure levels 30.,10.,5.) c Get NMC geopotential heights (cm) for year-day iyd (yyddd) c Return heights in z(nlon,nlat) at longitudes glon(nlon), c (at gcm latitude grid) c On input: c lu = fortran unit to use for necessary i/o ! mssfile = direct access mss file with nmc data ! mxrec = number of records in mssfile ! iyd0,iyd1 = first,last days (yyddd) on mssfile ! iyd = requested day (yyddd) c glon(nlon) = longitudes at which to get heights c nlat = 36 (for dimension of z) c On output: c z(nlon,nlat) is defined with the needed heights c ier = 0 if no error has occurred c Note: c mssfile is the mass store path to file containing NMC data c (see vishnu.hao:~foster/ccm2/nmclb). This is a fortran unformatted c direct access file containing coefficients for geopotential height c for the days 79001 through 93120, one record per day (5234 total c records). Each record (day) of data is 901 words, or 7208 bytes. c (See the read statement below) c 3/13/95: nmc data file updated to 94151 (5630 records) ! 6/10/99: pass in mssfile, mxrec, iyd0, and iyd1 from nmcbnd. c ! ! Args: integer,intent(in) :: lu,mxrec,iyd0,iyd1,iyd,nlon,nlat,npr,iprint real ,intent(in) :: glon(nlon) integer,intent(out) :: ier real ,intent(out) :: z(nlon,nlat,npr) character(len=32),intent(in) :: mssfile ! ! Local: integer,parameter :: nw=12,jmx=36 real,parameter :: pi=3.14159, dtr=pi/180. real :: zm(jmx,npr),zs(jmx,nw,npr),zc(jmx,nw,npr) real :: start(nlon), stop(nlon) character(len=24) :: diskfile = 'nmcz.dat ' ! character(len=24) :: diskfile = '98116-98120.dat ' character*128 errmsg integer :: ncalls=0 integer :: irecl,lmssfile,iier,istat,irec,iydrd,i,j,k,iw real :: zmmin,zmmax,zsmin,zsmax,zcmin,zcmax,alfa,angle,wt1,wt, | zmin,zmax integer,external :: iyddif c c Check nlat and lu: c ncalls = ncalls+1 ier = 0 if (nlat.ne.jmx) then write(6,"(' ')") write(6,"('>>> nmczday: bad nlat=',i4,' (should be ', + i2,')')") nlat,jmx ier = 1 return endif if (lu.lt.1.or.lu.gt.99.or.lu.eq.2.or.lu.eq.5.or.lu.eq.6) then write(6,"(' ')") write(6,"('>>> nmczday: bad lu=',i4)") lu if (lu.eq.2) write(6,"(' (lu=2 potentially conflicts ', + 'with ncar graphics unit)')") ier = 2 return endif c c Get and open direct access mss file: c (assume file needs to be acquired and opened only if this is first call) c irecl = 8*(1+nlat*npr+2*nlat*nw*npr) if (ncalls.eq.1) then ! get mss file lmssfile = len_trim(mssfile) if (iprint > 0) write(6,"('Nmczday: acquiring mss file ',a)") + mssfile(1:lmssfile) ! call msread(iier,diskfile,mssfile,' ',' ') if (iier.ne.3.and.iier.ne.0) then write(6,"(' ')") write(6,"('>>> nmczday: error from msread: mssfile=',a)") + mssfile(1:lmssfile) ! call mserror(errmsg) write(6,"(' mss error: ',a)") errmsg write(6,"(' ')") ier = 3 return endif open(lu,file=diskfile,access='DIRECT',status='OLD', + form='UNFORMATTED',recl=irecl,err=900,iostat=istat) endif c c Determine record number in file and do the read: c irec = iyddif(iyd0,iyd)+1 if (irec.le.0.or.irec.gt.mxrec) then write(6,"('>>> nmczday: bad irec for iyd=',i5,' irec=', + i5)") iyd,irec ier = 5 return endif ! ! dimension zm(jmx,npr),zs(jmx,nw,npr),zc(jmx,nw,npr) ! read(lu,rec=irec,err=901,iostat=istat) iydrd,zm,zs,zc if (iydrd.ne.iyd) then write(6,"(' ')") write(6,"('>>> nmczday: iydrd=',i5,' not equal to iyd=',i5, + ' (irec=',i5,')')") iydrd,iyd,irec ier = 7 return endif ! ! Report to stdout if requested: if (iprint > 0) then call fminmax(zm,nlat*npr,zmmin,zmmax) write(6,"('Nmczday: iyd=',i5,' zm(nlat,npr) min,max=',2e12.4)") | iyd,zmmin,zmmax do i=1,nw call fminmax(zs(:,i,:),nlat*npr,zsmin,zsmax) write(6,"('i=',i3,' nw=',i3,' zs(nlat,i,npr) min,max=', | 2e12.4)") i,nw,zsmin,zsmax call fminmax(zc(:,i,:),nlat*npr,zcmin,zcmax) write(6,"(14x,'zc(nlat,i,npr) min,max=',2e12.4)") zcmin,zcmax enddo endif c c Do the expansion to desired longitudes: C*************************************************** C C do j=1,jmx C do i=1,nlon C sum = 0. C do kk=1,nw C sum = sum + (zs(j,kk) * sin(float(kk)*glon(i)*dtr)+ C + zc(j,kk) * cos(float(kk)*glon(i)*dtr)) C enddo C z(i,j) = (zm(j)+sum) * 100. ! m to cm C enddo C enddo ALFA = 3. DO K=1,NPR DO J = 1,JMX DO I = 1,NLON Z(I,J,K) = ZM(J,K) ENDDO ENDDO ENDDO DO K=1,NPR DO J = 1,JMX ANGLE = ABS(-87.5 + (J-1)*5.)*2.*ALFA ANGLE = AMIN1(ANGLE,180.) WT1 = (1.-COS(ANGLE*DTR))/2. DO IW = 1,4 IF(IW.LE.2)THEN WT = 1. ELSE WT = WT1 ENDIF DO I = 1,NLON Z(I,J,K) = Z(I,J,K) + WT*(zs(j,iw,k) * sin(float(iw)* + glon(i)*dtr)+ zc(j,iw,k) * cos(float(iw)*glon(i)*dtr)) enddo ENDDO ENDDO ENDDO C **** C **** (1:2:1) latitudinal smoothing and x 100 for cm C **** DO K=1,NPR DO I = 1,NLON START(I) = 25.*(Z(I,1,K) + Z(1+MOD(I+35,72),1,K)) ENDDO DO J = 1,JMX-1 DO I = 1,NLON Z(I,J,K) = 25.*(Z(I,J,K) + Z(I,J+1,K)) ENDDO ENDDO DO I = 1,NLON STOP(I) = 25.*(Z(I,JMX,K) + Z(1+MOD(I+35,72),JMX,K)) ENDDO DO I = 1,NLON Z(I,JMX,K) = STOP(I) ENDDO DO J = JMX,2,-1 DO I = 1,NLON Z(I,J,K) = Z(I,J,K) + Z(I,J-1,K) ENDDO ENDDO DO I = 1,NLON Z(I,1,K) = Z(I,1,K) + START(I) ENDDO ENDDO if (iprint > 0) then call fminmax(z,nlon*nlat*npr,zmin,zmax) write(6,"('Nmczday output global z(nlon,nlat,npr) ', | 'min,max=',2e12.4)") zmin,zmax endif C call cpezct(z,nlon,jmx) C call cpezct(z(1,13),nlon,12) C*************************************************** return ! done c c i/o error traps: c 900 continue write(6,"(' ')") write(6,"('>>> nmczday: error opening file ',a,': ', + 'istat=',i5,' lu=',i2,' irecl=',i4)") diskfile,istat,lu,irecl ier = 4 return 901 continue write(6,"(' ')") write(6,"('>>> nmczday: error reading record ',i5, + ' from file ',a,': istat=',i4,' lu=',i2)") + irec,diskfile,istat,lu ier = 6 return end