WebSVN – lagranto.um – Diff – /tags/2.2/convert/pp2cdf_p20_Andy_Elvidge.f


      program pp2cdf
      implicit none
C
C     Purpose.
C     --------
C
C       Reads UKMO MM pressure data and writes the data to a NetCDF-file.
C       It works on 20 pressure levels, and uses stashcodes 90003 and 90004 
C       (u and v wrt model grid).
C
C     Interface.
C     ----------
C
C       File with MM data attached as unit 8.
C       File with NetCDF filename and missing data flag is unit 15.
C       Varfile (variable informations) attached as unit 12.
C
C     Reference.
C     ----------
C
C       None.
C
C     Author.
C     -------
C
C
C     Modifications.
C     --------------
C
C       None.
C
C     -----------------------------------------------------------------
C
 
CCC   Include your netcdf.inc
      include 'netcdf.inc'
 
C     ===  definitions of routinespecific variables  =======================
 
CCC   maximum array dimensions (change them for your purpose)
      integer   nxmax,nymax,nxi,nyi,nlevmax,nwet
      real      pmin,pmax,pdiff
      parameter(nxmax=270,nymax=270,nlevmax=20,nwet=20)
      parameter(nxi=nxmax-2,nyi=nymax-2)
 
      real      fld(nxmax,nymax)
      real	field(nxi*nyi)
      real	pstar(nxmax,nymax)
      real	data(nxi,nyi)
      real      varmin(3), varmax(3), stag(3)
      integer   vardim(4)
      real      mdv
      integer	datar(14),idate(5),stdate(5),cstid
      integer	lev,cdfid
      real	tstep
      real	aklev(nlevmax),bklev(nlevmax),aklay(nlevmax),bklay(nlevmax)
      real	latn,lats,lonw,lone,dx,dy,pole_lon,pole_lat
      real      p1,p2,pex1,pex2,exner,pres,plevs(nlevmax)
      integer	nx,ny,nlev,nvars,ntimes,tcode,qcode,pcode
      integer	i,j,ierr,iform,ieee,stash_code,iend,ilev
      integer	dyear,dmonth,dday,dhour,dmins
      logical	opnfile
      integer	strend,saveplevs,tind,irec
      integer	done_u,done_v,done_t,done_q,done_o,done_p
      character*80 cfn,wrcstflag,cdfname,varname
      logical    NewDyn
 
c     include 'ABmeso31.f'     ! define Hybrid coord. coefficients
 
C R IS GAS CONSTANT FOR DRY AIR
C CP IS SPECIFIC HEAT OF DRY AIR AT CONSTANT PRESSURE
C PREF IS REFERENCE SURFACE PRESSURE
      REAL R,CP,KAPPA,PREF  ! PREFTOKI
      PARAMETER(R=287.05,CP=1005.0,KAPPA=R/CP,PREF=100000.0)
!     PARAMETER(pmin=1000.0,pmax=200.0,pdiff=25.0,pcode=409)
 
      character*100 infile
      character*60 cdf_file,ctime*2,char1*1,ext1*2,ext2*2,date_str*2
 
      write(*,*)'===  start of program pp2cdf_p20 ===  (+++)'
      write(*,*)'===  This Program works on UKMM pressure level data'
 
C     Read parameters from input-file fort.15
 
      read(15,'(a)') infile		! name of PP Format input file
      read(15,'(a)') cdfname	! name of NetCDF file
      read(15,'(a)') wrcstflag	! write-cst-file flag
      read(15,'(a)') cfn           ! name of constants file
      read(15,*) mdv			! missing data value
      read(15,*) tcode             ! stash code of T/Theta field
      read(15,*) qcode             ! stash code of Q field
      read(15,*) pcode             ! stash code of PS field
      read(15,*) pmin              ! min value of PS field
      read(15,*) pmax              ! max. value of PS field
      read(15,*) pdiff             ! level difference of PS field
      read(15,*) NewDyn            ! New Dynamics (version 5 onwards)
C      write(*,*) NewDyn
C     write(*,*)pmin, pmax, pdiff 
C     Define constants file name
 
      if (cfn.eq.'default') then
        cfn=cdfname(1:strend(cdfname))//'_cst'
      endif
      opnfile=.false.
      nvars=0
      ntimes=1
 
C     ======================================================
C     Start of loop to write specified fields on NetCDF file
C     loop 30 --> goto 30
C     ======================================================
 
      iform=0
      ieee=1
      irec=0
      done_u=0
      done_v=0
      done_t=0
      done_q=0
      done_o=0
      done_p=0
      saveplevs=0
      nlev=pmax  !50.0
      call open_ppf(infile,iform,ieee)  ! open PP Format input file
 
   30 continue
 
CCC   call here a subroutine that reads 1 variable on 1 level from
CCC   a data file
CCC   call another subroutine or do everything in one, as you like,
CCC   that gets some general information about the data
CCC   anyhow, what you do here should define the following parameters:
      
C     dyear	integer		year of date valid for current field
C     dmonth	integer		month of date valid for current field
C     dday	integer		day of date valid for current field
C     dhour	integer		hour of date valid for current field
C     dmins	integer		minute of date valid for current field
      
C     lonw	real		coordinates of data domain (in deg):
C     lone	real		      latn
C     lats	real       	 lonw --|-- lone
C     late	real		      lats
C     dx	real		grid size in W-E direction (in deg)
C     dy	real		grid size in S-N direction (in deg)
C     nx	integer	number of data points in W-E direction
C     ny	integer	number of data points in S-N direction
C     nlev	integer	number of data levels
 
C     stash_code	integer	your field code
C     lev	integer	number of the level for the current field
C     fld 	real		2-dim array with dimensions (nx,ny),
C				contains the field values
 
CCC   ... call here your subroutine(s)...
      call read_field(iform,nx,ny,fld,stash_code,lev,
     &     dx,dy,dyear,dmonth,dday,dhour,dmins,iend,
     &     lonw,lone,latn,lats,pole_lon,pole_lat,pres,nxmax,nymax,
     &     NewDyn)
c      write(*,*)'stash and level',stash_code,lev,iend,pres
      if(iend.eq.1) goto 999
 
CCC   Determine the name of the variable from the field code 'num'
CCC   If a variable is not in the list, it is not written to the NetCDF file
      if (stash_code.eq.90003) then
        varname='U'			! zonal velocity
        if(lev.eq.nlev) done_u=1
c       if(saveplevs.lt.nlevmax+1) then
c         if(saveplevs.eq.0) tind=1
c         plevs(tind)=pres
c         tind=tind+1
c       endif
      else if (stash_code.eq.90004) then
        varname='V'			! meridional velocity
        if(lev.eq.nlev) done_v=1
      else if (stash_code.eq.90087) then
        varname='OMEGA'			! vertical velocity
        if(lev.eq.nlev) done_o=1
      else if (stash_code.eq.tcode) then
        varname='T'			! temperature
        if(lev.eq.nlev) done_t=1
c if Theta supplied then convert to Temp
        if(tcode.eq.90006) then
          do j = 1,ny
            do i = 1,nx
c     form Temp. at full levels from Theta
              exner=(pres/pref)**kappa
              fld(i,j)=fld(i,j)*exner
            enddo
          enddo
        endif
      else if (stash_code.eq.pcode ) then
        varname='PS'			! surface pressure
cwang   nlev=1
cwang   lev=1        ! reset lev to 1 from 9999
        done_p=1
        do j = 1,ny
          do i = 1,nx
            pstar(i,j)=fld(i,j)
            fld(i,j)=fld(i,j)/100.0    ! Pstar needed in hPa
          enddo
        enddo
      else if (stash_code.eq.qcode) then
        varname='Q'			! specific humidity
        if(lev.eq.nlev) done_q=1
      else
        goto 30				! don't put fld on NetCDF file
      endif
 
      irec=irec+1
c     N.B. U, V and Omega are all on P-Grid
        do i=1,nxi
          do j=1,nyi
            data(i,j)=fld(i+1,j+1)
          enddo
        enddo
 
C     Define staggering coefficients
CCC   change the ??? to your num-value for vertical velocity
 
      stag(1)=0.0
      stag(2)=0.0
      stag(3)=0.0
c     if (stash_code.eq.???) then	! no staggering for vertical velocity
c       stag(3)=0.0
c     else
c       stag(3)=-0.5		! vertical staggering for all other variables
c     endif
 
C     Define data region
 
      varmin(1)=lonw
      varmin(2)=lats
      varmax(1)=lone
      varmax(2)=latn
C     varmin/max(3) changed to take into account reverse order of data
      varmin(3)=pmin   ! Not used
      varmax(3)=pmax  !50.0
      write(*,*) 'lonw,lats,lone,latn',lonw,lats,lone,latn
      write(*,*) 'pmax, pmin',pmax, pmin 
 
C     Define dimensions
 
      vardim(1)=nxi
      vardim(2)=nyi
CCC   Distinction between 2dim and 3dim fields
CCC   Insert num-values for 2dim fields
      if (stash_code.eq.pcode) then
        vardim(3)=1
      else
        vardim(3)=nlevmax
      endif
      vardim(4)=1
 
C     Define the date array
 
      idate(1)=dyear
      idate(2)=dmonth
      idate(3)=dday
      idate(4)=dhour
      idate(5)=dmins
 
C     set stdate on first read of PP record from file
 
c     if (irec.eq.1) then
c     endif
      if (wrcstflag.eq.'yes') then
 
C       stdate is the simulation start time
        stdate(1)=dyear
        stdate(2)=dmonth
        stdate(3)=dday
        stdate(4)=dhour
        stdate(5)=dmins
C       Define the datar array (used to define constants file)
 
        datar(1)=nxi
        datar(2)=nyi
        datar(3)=nint(1000.0*latn) ! I agree this is not nice
        datar(4)=nint(1000.0*lonw) ! but for historical reasons
        datar(5)=nint(1000.0*lats) ! subroutine wricst contains
        datar(6)=nint(1000.0*lone) ! this factor 1000. things
        datar(7)=nint(1000.0*dx)
        datar(8)=nint(1000.0*dy)
        datar(9)=nlevmax  ! just in case Pstar first in file
        datar(10)=0         ! not used
        datar(11)=0         ! not used
        datar(12)=0         ! not used
        datar(13)=nint(1000.0*pole_lon)
        datar(14)=nint(1000.0*pole_lat)
 
C       Get the 4 levels/layers arrays
 
        call modlevs(aklev,bklev,aklay,bklay,plevs,nlevmax,pmin,pdiff)
 
C       Write the constants file
 
        call wricst(cfn,datar,aklev,bklev,aklay,bklay,stdate)
        write(*,*)
        write(*,*)'*** constants-file ',cfn(1:strend(cfn)),' created'
        wrcstflag='no'
      else           ! assume that the constants-file already exists
C       Inquire start time
        call cdfopn(cfn,cstid,ierr)              ! open constants file
        call getstart(cstid,stdate,ierr)
        call ncclos(cstid,ierr)                  ! close constants file
      endif
 
C     Create the NetCDF file (only once)
 
      if (.not.opnfile) then
c    Code to derive Index No. for filenames
c       if(ntimes.lt.10) then
c         write(char1,'(i1)') ntimes
c         ctime='0'//char1
c       else
c         write(ctime,'(i2)') ntimes
c       endif
c       cdf_file=cdfname(1:strend(cdfname))//ctime
c    Code to derive timestamp for filenames
c       dhour=dhour-stdate(4)  ! subtract start time so always atart at 0
        if(dhour.eq.0) then
          ext1='00'
        elseif(dhour.lt.10) then
          write(char1,'(i1)') dhour
          ext1='0'//char1
        else
          write(ext1,'(i2)') dhour
        endif
        if(dmins.eq.0) then
          ext2='00'
        else
          write(ext2,'(i2)') dmins
        endif
c       n.b if data runs into next day then can't use root of first
c       filename for rest, so form filename manually
c       cdf_file=cdfname(1:strend(cdfname))//ext1//ext2
c       below assumes year and month are constant
        if(dday.lt.10) then
          write(char1,'(i1)') dday
          date_str='0'//char1
        else
          write(date_str,'(i2)') dday
        endif
        cdf_file=cdfname(1:5)//date_str//'_'//ext1//ext2
        call crecdf(cdf_file,cdfid,varmin,varmax,
     &              3,cfn(1:strend(cfn)),ierr)
        if (ierr.ne.0) goto 996
        write(*,*)
c       write(*,*)'*** NetCDF file ',cdfname(1:strend(cdfname)),' create
        write(*,*)'*** NetCDF file ',cdf_file,' created'
        opnfile=.true.
      endif
 
C     Put 2dim array in 1dim array
 
      do i=1,nxi
        do j=1,nyi
          field(i+(j-1)*nxi)=data(i,j)
        enddo
      enddo
 
C     Calculate tstep for actual date (relative to first time of
c     simulation) and convert to minutes
 
      call timediff(idate,stdate,tstep)
c     tstep=tstep*100.0
c     tstep=float(nint(tstep*100.0))
      write(6,*)"idate,stdate,tstep:",idate,stdate,tstep
      i=ifix(tstep)  ! round down to nearest hour
      tstep=float(i*60+nint(100.0*tstep)-100*i)
cwang     ilev=1+((1000-lev)/50) this should be changed to /25     
cwang     change ilev=1+((1000-lev)/50) statement when pdiff read in
      write(6,*)"tstep:",tstep
      ilev=1+((pmin-lev)/pdiff)
      if(stash_code.eq.pcode) ilev=1
 
      if (ilev.eq.1) then
        call putdef(cdfid,varname(1:strend(varname)),4,mdv,vardim,
     &              varmin,varmax,stag,ierr)
        if (ierr.ne.0) goto 996
        write(*,*)
        write(*,*)'*** variable ',varname(1:strend(varname)),' created'
        nvars=nvars+1
      endif
 
C     Put data on file
 
c     Convert units of Temp to Degrees C
      if (stash_code.eq.90002.or.stash_code.eq.90006) then
        do i=1,nxi*nyi
          field(i)=field(i)-273.15
        enddo
      endif
 
      write(*,*)
      write(*,'(a,a,a,f7.1,a,i4,a)')'*** var ',
     &  varname(1:strend(varname)),' at time ',
     &  tstep,' on level ',lev,' written on cdffile'
      call putdat(cdfid,varname(1:strend(varname)),tstep,ilev,
     &  field,ierr)
      if (ierr.ne.0) goto 996
      write(*,*)'cdfid',cdfid
c     if(nvars.eq.6.and.lev.eq.nlev) then   ! required variables written
      write(*,*)done_u,done_v,done_t,done_q,done_o,done_p
      if(done_u+done_v+done_t+done_q+done_o+done_p.eq.6) then
        opnfile=.false.
        ntimes=ntimes+1
        nvars=0
C     Close NetCDF file
 
        call ncclos(cdfid,ierr)
        if (ierr.ne.0) stop 'error when closing NetCDF file'
c       if(ntimes.eq.2) stop
        done_u=0
        done_v=0
        done_t=0
        done_q=0
        done_o=0
        done_p=0
      endif
      goto 30
 
  996 stop 'error: could not create NetCDF file'
  999 continue
      write(*,*)lev,nlev,ierr
      if(lev.gt.nlev.or.ierr.ne.0) then
        write(*,*)'===  Error: NetCDF file: ',cdf_file,' incomplete'
        write(*,*)'See fort.12 for list of data processed'
        stop
      endif
 
      write(*,*)'===  end of program pp2cdf_p20: status normal  ==='
 
      stop
      end
 
      subroutine open_ppf(infile,iform,ieee)
C    I/O STREAMS ARE AS FOLLOWS
C     1) 6   STANDARD OUTPUT
C     2) 8   PPUNIT - READ UNIT FOR PP FILE
      IMPLICIT NONE
      INTEGER ierr,IFORM,IEEE
      CHARACTER INFILE*100
c
      IF(IFORM.EQ.1) THEN
        OPEN(UNIT=8,FILE=INFILE,FORM='FORMATTED')
      ELSE
        IF(IEEE.EQ.1) THEN
c     2 lines below are Cray specific calls to identify IEEE file
c     n.b 'newlocal' refers to assign environment, not file
c         call asnctl('newlocal',1,ierr)
c         write(*,*) 'ierr from asnctl ',ierr
c         call asnunit(8,'-F f77 -N ieee',ierr)
c         write(*,*) 'ierr from asnunit ',ierr
        ENDIF
        OPEN(UNIT=8,FILE=INFILE,FORM='UNFORMATTED')
      ENDIF
      return
      END
c********************************************************
      subroutine read_field(iform,nx,ny,data,stash_code,level,
     &     dxdeg,dydeg,dyear,dmonth,dday,dhour,dmins,iend,
     &     lonw,lone,latn,lats,pole_lon,pole_lat,pres,nxmax,nymax,
     &      NewDyn)
C    I/O STREAMS ARE AS FOLLOWS
C     1) 6   STANDARD OUTPUT
C     2) 8   PPUNIT - READ UNIT FOR PP FILE
C     3) 12  HEADER INFORMATION STREAM
      IMPLICIT NONE
c     
      REAL ROOK(19), FIELD(72900),data(270,270)
      INTEGER LOOK(45),I,J,K,N,NUM_VALUES,NX,NY,IFORM,nxmax,nymax,
     &     nf,stash_code,level,dyear,dmonth,dday,dhour,dmins,iend
      real LONG1,LAT1,dxdeg,dydeg,long_nw,lat_nw,
     &     fac,lonw,lone,latn,lats,pole_lon,pole_lat,pres
      logical NewDyn
      data nf/0/
      save nf
      nf=nf+1
c
      IF(IFORM.EQ.1) THEN
        READ(8,12,END=99) LOOK
        READ(8,13) ROOK
      ELSE
        READ(8,END=99) LOOK,ROOK
      ENDIF
      NX=LOOK(19)
      NY=LOOK(18)
      NUM_VALUES=NX*NY
C     WRITE OUT CONTENTS OF INFILE FOR INFORMATION
      WRITE(12,106) LOOK(33),LOOK(14),nf,
     &      NUM_VALUES,LOOK(42),LOOK(18),LOOK(19)
      IF(IFORM.EQ.1) THEN
        READ(8,13) (FIELD(I),I=1,NUM_VALUES)
      ELSE
      write(*,*) '1'
        READ(8) (FIELD(I),I=1,NUM_VALUES)
      ENDIF
      write(*,*) '2'
      LEVEL=LOOK(33)
      stash_code=LOOK(42)
c     Assign exact time of field from "Validity time"
      dyear=look(1)
      if(dyear.gt.1900) dyear=dyear-1900   ! wants only year as 9?
      dmonth=look(2)
      dday=look(3)
      dhour=look(4)
      dmins=look(5)
      LONG1=ROOK(61-45)
      if(LONG1 .lt. 0.0 ) then 
        LONG1=LONG1+360.0
      endif
      LAT1=ROOK(59-45)
      dxdeg=abs(ROOK(62-45))
      dydeg=abs(ROOK(60-45))
cwang      write(*,*)'dxdy',ROOK(62-45),ROOK(60-45)
cwang      write(*,*)'LONG1Lat1',ROOK(61-45),ROOK(59-45)
      fac=0.0
c     if(stash_code.eq.2.or.stash_code.eq.3.or.stash_code.eq.12201)
c    &      fac=0.5
C PAC Verified 
      long_nw=LONG1+dxdeg-fac*dxdeg
cwang      lat_nw=Lat1-dydeg+fac*dydeg
c     define limits of inner 2:nxmax-1 X 2:nymax-1 grid
      lonw=long_nw+dxdeg
      lone=lonw+float(nxmax-3)*dxdeg
cwang change statements to calculate latn and lats
cwang to next if loop so that it counts for dy for positive condition
      if (NewDyn .eq. .false.) then 
        lat_nw=Lat1-dydeg+fac*dydeg
        latn=lat_nw-dydeg
        lats=latn-float(nymax-3)*dydeg
      else
C PAC Verified
        lat_nw=Lat1+dydeg-fac*dydeg
        lats=lat_nw+dydeg
        latn=lats+float(nymax-3)*dydeg
      endif
      pole_lon=ROOK(57-45)
      pole_lat=ROOK(56-45)
      pres=100.0*rook(52-45)  ! in Pa
C-------------------------------------------------------------------
 12   FORMAT(12I10)
 13   FORMAT(10E12.5)
  106 FORMAT(' LEVEL=',I5,' FCST=',I6,' FIELD NO',I4,' NVALS='
     &    ,I5,' STASH CODE=',I5,' NO.ROWS=',I3,' NO.COLUMNS=',I3)
  100 FORMAT('  THE NUMBER OF FIELDS IN THE INPUT FILE IS   ',I8)
C-----------------------Process Data--------------------------------
C     Data WRITTEN STARTING AT N-W CORNER, FOLLOWING LOOPS
C     ASSIGN DATA FROM S-W CORNER AND THEN WRITE OUT
      if (NewDyn .eq. .false.) then 
        write(*,*)'data read starting at N-W'
        DO J=NY-1,0,-1 !LOOP SOUTH-NORTH (ROWS)
          I=J*NX       !FIELD RUNS OPPOSITE
          DO N=1,NX
            DATA(N,NY-J)=FIELD(N+I)
          END DO
        END DO
      else
C PAC Verified
         write(*,*)'data read dstarting at S-W'
        DO J=1,NY
          DO N=1,NX
            DATA(N,J)=FIELD((J-1)*NX+N)
          END DO
        END DO
      endif
C     CODE TO PREVENT UNDERFLOWS ON CONVERSION TO IEEE
      DO J = 1,NY
        DO I = 1,NX
           IF(ABS(DATA(I,J)).LT.1.0E-30.AND.DATA(I,J).NE.0.0) THEN
c            WRITE(6,*) 'UNDERFLOW AT: ',I,J,DATA(I,J)
             DATA(I,J) = 0.0
           ENDIF
        END DO
      END DO
      iend=0
      return
 99   continue
      iend=1
      write(*,*) 'End of file reached'
      return
      END
      subroutine modlevs(aklev,bklev,aklay,bklay,plevs,nlevmax,pm,pd)
C------------------------------------------------------------------------
C
C     Defines the level- and layer-arrays given the number of levels nlev.
C
C     nlev    int    input  number of levels/layers
C     aklev   real   output array contains ak values for levels
C     bklev     real    output  array contains bk values for levels
C     aklay     real    output  array contains ak values for layers
C     bklay     real    output  array contains bk values for layers
C------------------------------------------------------------------------
 
      integer   nlev,k
 
      real      aklev(nlevmax),bklev(nlevmax),    ! level coefficients
     &          aklay(nlevmax),bklay(nlevmax)     ! layer coefficients
     &         ,plevs(nlevmax)
 
      do k=1,nlevmax
        aklay(k)=0.0
        bklay(k)=0.0
c       aklev(k)=plevs(k)/100.0    ! P needed in hPa
        aklev(k)=(pm-(k-1)*pd)    ! P needed in hPa
        bklev(k)=0.0
      enddo
      return
      end

Rev 3	Rev 14
1	`program pp2cdf`	1	`program pp2cdf`
2	`implicit none`	2	`implicit none`
3	`C`	3	`C`
4	`C Purpose.`	4	`C Purpose.`
5	`C --------`	5	`C --------`
6	`C`	6	`C`
7	`C Reads UKMO MM pressure data and writes the data to a NetCDF-file.`	7	`C Reads UKMO MM pressure data and writes the data to a NetCDF-file.`
8	`C It works on 20 pressure levels, and uses stashcodes 90003 and 90004`	8	`C It works on 20 pressure levels, and uses stashcodes 90003 and 90004`
9	`C (u and v wrt model grid).`	9	`C (u and v wrt model grid).`
10	`C`	10	`C`
11	`C Interface.`	11	`C Interface.`
12	`C ----------`	12	`C ----------`
13	`C`	13	`C`
14	`C File with MM data attached as unit 8.`	14	`C File with MM data attached as unit 8.`
15	`C File with NetCDF filename and missing data flag is unit 15.`	15	`C File with NetCDF filename and missing data flag is unit 15.`
16	`C Varfile (variable informations) attached as unit 12.`	16	`C Varfile (variable informations) attached as unit 12.`
17	`C`	17	`C`
18	`C Reference.`	18	`C Reference.`
19	`C ----------`	19	`C ----------`
20	`C`	20	`C`
21	`C None.`	21	`C None.`
22	`C`	22	`C`
23	`C Author.`	23	`C Author.`
24	`C -------`	24	`C -------`
25	`C`	25	`C`
26	`C`	26	`C`
27	`C Modifications.`	27	`C Modifications.`
28	`C --------------`	28	`C --------------`
29	`C`	29	`C`
30	`C None.`	30	`C None.`
31	`C`	31	`C`
32	`C -----------------------------------------------------------------`	32	`C -----------------------------------------------------------------`
33	`C`	33	`C`
34		34
35	`CCC Include your netcdf.inc`	35	`CCC Include your netcdf.inc`
36	`include 'netcdf.inc'`	36	`include 'netcdf.inc'`
37		37
38	`C === definitions of routinespecific variables =======================`	38	`C === definitions of routinespecific variables =======================`
39		39
40	`CCC maximum array dimensions (change them for your purpose)`	40	`CCC maximum array dimensions (change them for your purpose)`
41	`integer nxmax,nymax,nxi,nyi,nlevmax,nwet`	41	`integer nxmax,nymax,nxi,nyi,nlevmax,nwet`
42	`real pmin,pmax,pdiff`	42	`real pmin,pmax,pdiff`
43	`parameter(nxmax=270,nymax=270,nlevmax=20,nwet=20)`	43	`parameter(nxmax=270,nymax=270,nlevmax=20,nwet=20)`
44	`parameter(nxi=nxmax-2,nyi=nymax-2)`	44	`parameter(nxi=nxmax-2,nyi=nymax-2)`
45		45
46	`real fld(nxmax,nymax)`	46	`real fld(nxmax,nymax)`
47	`real field(nxi*nyi)`	47	`real field(nxi*nyi)`
48	`real pstar(nxmax,nymax)`	48	`real pstar(nxmax,nymax)`
49	`real data(nxi,nyi)`	49	`real data(nxi,nyi)`
50	`real varmin(3), varmax(3), stag(3)`	50	`real varmin(3), varmax(3), stag(3)`
51	`integer vardim(4)`	51	`integer vardim(4)`
52	`real mdv`	52	`real mdv`
53	`integer datar(14),idate(5),stdate(5),cstid`	53	`integer datar(14),idate(5),stdate(5),cstid`
54	`integer lev,cdfid`	54	`integer lev,cdfid`
55	`real tstep`	55	`real tstep`
56	`real aklev(nlevmax),bklev(nlevmax),aklay(nlevmax),bklay(nlevmax)`	56	`real aklev(nlevmax),bklev(nlevmax),aklay(nlevmax),bklay(nlevmax)`
57	`real latn,lats,lonw,lone,dx,dy,pole_lon,pole_lat`	57	`real latn,lats,lonw,lone,dx,dy,pole_lon,pole_lat`
58	`real p1,p2,pex1,pex2,exner,pres,plevs(nlevmax)`	58	`real p1,p2,pex1,pex2,exner,pres,plevs(nlevmax)`
59	`integer nx,ny,nlev,nvars,ntimes,tcode,qcode,pcode`	59	`integer nx,ny,nlev,nvars,ntimes,tcode,qcode,pcode`
60	`integer i,j,ierr,iform,ieee,stash_code,iend,ilev`	60	`integer i,j,ierr,iform,ieee,stash_code,iend,ilev`
61	`integer dyear,dmonth,dday,dhour,dmins`	61	`integer dyear,dmonth,dday,dhour,dmins`
62	`logical opnfile`	62	`logical opnfile`
63	`integer strend,saveplevs,tind,irec`	63	`integer strend,saveplevs,tind,irec`
64	`integer done_u,done_v,done_t,done_q,done_o,done_p`	64	`integer done_u,done_v,done_t,done_q,done_o,done_p`
65	`character*80 cfn,wrcstflag,cdfname,varname`	65	`character*80 cfn,wrcstflag,cdfname,varname`
66	`logical NewDyn`	66	`logical NewDyn`
67		67
68	`c include 'ABmeso31.f' ! define Hybrid coord. coefficients`	68	`c include 'ABmeso31.f' ! define Hybrid coord. coefficients`
69		69
70	`C R IS GAS CONSTANT FOR DRY AIR`	70	`C R IS GAS CONSTANT FOR DRY AIR`
71	`C CP IS SPECIFIC HEAT OF DRY AIR AT CONSTANT PRESSURE`	71	`C CP IS SPECIFIC HEAT OF DRY AIR AT CONSTANT PRESSURE`
72	`C PREF IS REFERENCE SURFACE PRESSURE`	72	`C PREF IS REFERENCE SURFACE PRESSURE`
73	`REAL R,CP,KAPPA,PREF ! PREFTOKI`	73	`REAL R,CP,KAPPA,PREF ! PREFTOKI`
74	`PARAMETER(R=287.05,CP=1005.0,KAPPA=R/CP,PREF=100000.0)`	74	`PARAMETER(R=287.05,CP=1005.0,KAPPA=R/CP,PREF=100000.0)`
75	`! PARAMETER(pmin=1000.0,pmax=200.0,pdiff=25.0,pcode=409)`	75	`! PARAMETER(pmin=1000.0,pmax=200.0,pdiff=25.0,pcode=409)`
76		76
77	`character*100 infile`	77	`character*100 infile`
78	`character60 cdf_file,ctime2,char11,ext12,ext22,date_str2`	78	`character60 cdf_file,ctime2,char11,ext12,ext22,date_str2`
79		79
80	`write(,)'=== start of program pp2cdf_p20 === (+++)'`	80	`write(,)'=== start of program pp2cdf_p20 === (+++)'`
81	`write(,)'=== This Program works on UKMM pressure level data'`	81	`write(,)'=== This Program works on UKMM pressure level data'`
82		82
83	`C Read parameters from input-file fort.15`	83	`C Read parameters from input-file fort.15`
84		84
85	`read(15,'(a)') infile ! name of PP Format input file`	85	`read(15,'(a)') infile ! name of PP Format input file`
86	`read(15,'(a)') cdfname ! name of NetCDF file`	86	`read(15,'(a)') cdfname ! name of NetCDF file`
87	`read(15,'(a)') wrcstflag ! write-cst-file flag`	87	`read(15,'(a)') wrcstflag ! write-cst-file flag`
88	`read(15,'(a)') cfn ! name of constants file`	88	`read(15,'(a)') cfn ! name of constants file`
89	`read(15,*) mdv ! missing data value`	89	`read(15,*) mdv ! missing data value`
90	`read(15,*) tcode ! stash code of T/Theta field`	90	`read(15,*) tcode ! stash code of T/Theta field`
91	`read(15,*) qcode ! stash code of Q field`	91	`read(15,*) qcode ! stash code of Q field`
92	`read(15,*) pcode ! stash code of PS field`	92	`read(15,*) pcode ! stash code of PS field`
93	`read(15,*) pmin ! min value of PS field`	93	`read(15,*) pmin ! min value of PS field`
94	`read(15,*) pmax ! max. value of PS field`	94	`read(15,*) pmax ! max. value of PS field`
95	`read(15,*) pdiff ! level difference of PS field`	95	`read(15,*) pdiff ! level difference of PS field`
96	`read(15,*) NewDyn ! New Dynamics (version 5 onwards)`	96	`read(15,*) NewDyn ! New Dynamics (version 5 onwards)`
97	`C write(,) NewDyn`	97	`C write(,) NewDyn`
98	`C write(,)pmin, pmax, pdiff`	98	`C write(,)pmin, pmax, pdiff`
99	`C Define constants file name`	99	`C Define constants file name`
100		100
101	`if (cfn.eq.'default') then`	101	`if (cfn.eq.'default') then`
102	`cfn=cdfname(1:strend(cdfname))//'_cst'`	102	`cfn=cdfname(1:strend(cdfname))//'_cst'`
103	`endif`	103	`endif`
104	`opnfile=.false.`	104	`opnfile=.false.`
105	`nvars=0`	105	`nvars=0`
106	`ntimes=1`	106	`ntimes=1`
107		107
108	`C ======================================================`	108	`C ======================================================`
109	`C Start of loop to write specified fields on NetCDF file`	109	`C Start of loop to write specified fields on NetCDF file`
110	`C loop 30 --> goto 30`	110	`C loop 30 --> goto 30`
111	`C ======================================================`	111	`C ======================================================`
112		112
113	`iform=0`	113	`iform=0`
114	`ieee=1`	114	`ieee=1`
115	`irec=0`	115	`irec=0`
116	`done_u=0`	116	`done_u=0`
117	`done_v=0`	117	`done_v=0`
118	`done_t=0`	118	`done_t=0`
119	`done_q=0`	119	`done_q=0`
120	`done_o=0`	120	`done_o=0`
121	`done_p=0`	121	`done_p=0`
122	`saveplevs=0`	122	`saveplevs=0`
123	`nlev=pmax !50.0`	123	`nlev=pmax !50.0`
124	`call open_ppf(infile,iform,ieee) ! open PP Format input file`	124	`call open_ppf(infile,iform,ieee) ! open PP Format input file`
125		125
126	`30 continue`	126	`30 continue`
127		127
128	`CCC call here a subroutine that reads 1 variable on 1 level from`	128	`CCC call here a subroutine that reads 1 variable on 1 level from`
129	`CCC a data file`	129	`CCC a data file`
130	`CCC call another subroutine or do everything in one, as you like,`	130	`CCC call another subroutine or do everything in one, as you like,`
131	`CCC that gets some general information about the data`	131	`CCC that gets some general information about the data`
132	`CCC anyhow, what you do here should define the following parameters:`	132	`CCC anyhow, what you do here should define the following parameters:`
133		133
134	`C dyear integer year of date valid for current field`	134	`C dyear integer year of date valid for current field`
135	`C dmonth integer month of date valid for current field`	135	`C dmonth integer month of date valid for current field`
136	`C dday integer day of date valid for current field`	136	`C dday integer day of date valid for current field`
137	`C dhour integer hour of date valid for current field`	137	`C dhour integer hour of date valid for current field`
138	`C dmins integer minute of date valid for current field`	138	`C dmins integer minute of date valid for current field`
139		139
140	`C lonw real coordinates of data domain (in deg):`	140	`C lonw real coordinates of data domain (in deg):`
141	`C lone real latn`	141	`C lone real latn`
142	`C lats real lonw --\|-- lone`	142	`C lats real lonw --\|-- lone`
143	`C late real lats`	143	`C late real lats`
144	`C dx real grid size in W-E direction (in deg)`	144	`C dx real grid size in W-E direction (in deg)`
145	`C dy real grid size in S-N direction (in deg)`	145	`C dy real grid size in S-N direction (in deg)`
146	`C nx integer number of data points in W-E direction`	146	`C nx integer number of data points in W-E direction`
147	`C ny integer number of data points in S-N direction`	147	`C ny integer number of data points in S-N direction`
148	`C nlev integer number of data levels`	148	`C nlev integer number of data levels`
149		149
150	`C stash_code integer your field code`	150	`C stash_code integer your field code`
151	`C lev integer number of the level for the current field`	151	`C lev integer number of the level for the current field`
152	`C fld real 2-dim array with dimensions (nx,ny),`	152	`C fld real 2-dim array with dimensions (nx,ny),`
153	`C contains the field values`	153	`C contains the field values`
154		154
155	`CCC ... call here your subroutine(s)...`	155	`CCC ... call here your subroutine(s)...`
156	`call read_field(iform,nx,ny,fld,stash_code,lev,`	156	`call read_field(iform,nx,ny,fld,stash_code,lev,`
157	`& dx,dy,dyear,dmonth,dday,dhour,dmins,iend,`	157	`& dx,dy,dyear,dmonth,dday,dhour,dmins,iend,`
158	`& lonw,lone,latn,lats,pole_lon,pole_lat,pres,nxmax,nymax,`	158	`& lonw,lone,latn,lats,pole_lon,pole_lat,pres,nxmax,nymax,`
159	`& NewDyn)`	159	`& NewDyn)`
160	`c write(,)'stash and level',stash_code,lev,iend,pres`	160	`c write(,)'stash and level',stash_code,lev,iend,pres`
161	`if(iend.eq.1) goto 999`	161	`if(iend.eq.1) goto 999`
162		162
163	`CCC Determine the name of the variable from the field code 'num'`	163	`CCC Determine the name of the variable from the field code 'num'`
164	`CCC If a variable is not in the list, it is not written to the NetCDF file`	164	`CCC If a variable is not in the list, it is not written to the NetCDF file`
165	`if (stash_code.eq.90003) then`	165	`if (stash_code.eq.90003) then`
166	`varname='U' ! zonal velocity`	166	`varname='U' ! zonal velocity`
167	`if(lev.eq.nlev) done_u=1`	167	`if(lev.eq.nlev) done_u=1`
168	`c if(saveplevs.lt.nlevmax+1) then`	168	`c if(saveplevs.lt.nlevmax+1) then`
169	`c if(saveplevs.eq.0) tind=1`	169	`c if(saveplevs.eq.0) tind=1`
170	`c plevs(tind)=pres`	170	`c plevs(tind)=pres`
171	`c tind=tind+1`	171	`c tind=tind+1`
172	`c endif`	172	`c endif`
173	`else if (stash_code.eq.90004) then`	173	`else if (stash_code.eq.90004) then`
174	`varname='V' ! meridional velocity`	174	`varname='V' ! meridional velocity`
175	`if(lev.eq.nlev) done_v=1`	175	`if(lev.eq.nlev) done_v=1`
176	`else if (stash_code.eq.90087) then`	176	`else if (stash_code.eq.90087) then`
177	`varname='OMEGA' ! vertical velocity`	177	`varname='OMEGA' ! vertical velocity`
178	`if(lev.eq.nlev) done_o=1`	178	`if(lev.eq.nlev) done_o=1`
179	`else if (stash_code.eq.tcode) then`	179	`else if (stash_code.eq.tcode) then`
180	`varname='T' ! temperature`	180	`varname='T' ! temperature`
181	`if(lev.eq.nlev) done_t=1`	181	`if(lev.eq.nlev) done_t=1`
182	`c if Theta supplied then convert to Temp`	182	`c if Theta supplied then convert to Temp`
183	`if(tcode.eq.90006) then`	183	`if(tcode.eq.90006) then`
184	`do j = 1,ny`	184	`do j = 1,ny`
185	`do i = 1,nx`	185	`do i = 1,nx`
186	`c form Temp. at full levels from Theta`	186	`c form Temp. at full levels from Theta`
187	`exner=(pres/pref)**kappa`	187	`exner=(pres/pref)**kappa`
188	`fld(i,j)=fld(i,j)*exner`	188	`fld(i,j)=fld(i,j)*exner`
189	`enddo`	189	`enddo`
190	`enddo`	190	`enddo`
191	`endif`	191	`endif`
192	`else if (stash_code.eq.pcode ) then`	192	`else if (stash_code.eq.pcode ) then`
193	`varname='PS' ! surface pressure`	193	`varname='PS' ! surface pressure`
194	`cwang nlev=1`	194	`cwang nlev=1`
195	`cwang lev=1 ! reset lev to 1 from 9999`	195	`cwang lev=1 ! reset lev to 1 from 9999`
196	`done_p=1`	196	`done_p=1`
197	`do j = 1,ny`	197	`do j = 1,ny`
198	`do i = 1,nx`	198	`do i = 1,nx`
199	`pstar(i,j)=fld(i,j)`	199	`pstar(i,j)=fld(i,j)`
200	`fld(i,j)=fld(i,j)/100.0 ! Pstar needed in hPa`	200	`fld(i,j)=fld(i,j)/100.0 ! Pstar needed in hPa`
201	`enddo`	201	`enddo`
202	`enddo`	202	`enddo`
203	`else if (stash_code.eq.qcode) then`	203	`else if (stash_code.eq.qcode) then`
204	`varname='Q' ! specific humidity`	204	`varname='Q' ! specific humidity`
205	`if(lev.eq.nlev) done_q=1`	205	`if(lev.eq.nlev) done_q=1`
206	`else`	206	`else`
207	`goto 30 ! don't put fld on NetCDF file`	207	`goto 30 ! don't put fld on NetCDF file`
208	`endif`	208	`endif`
209		209
210	`irec=irec+1`	210	`irec=irec+1`
211	`c N.B. U, V and Omega are all on P-Grid`	211	`c N.B. U, V and Omega are all on P-Grid`
212	`do i=1,nxi`	212	`do i=1,nxi`
213	`do j=1,nyi`	213	`do j=1,nyi`
214	`data(i,j)=fld(i+1,j+1)`	214	`data(i,j)=fld(i+1,j+1)`
215	`enddo`	215	`enddo`
216	`enddo`	216	`enddo`
217		217
218	`C Define staggering coefficients`	218	`C Define staggering coefficients`
219	`CCC change the ??? to your num-value for vertical velocity`	219	`CCC change the ??? to your num-value for vertical velocity`
220		220
221	`stag(1)=0.0`	221	`stag(1)=0.0`
222	`stag(2)=0.0`	222	`stag(2)=0.0`
223	`stag(3)=0.0`	223	`stag(3)=0.0`
224	`c if (stash_code.eq.???) then ! no staggering for vertical velocity`	224	`c if (stash_code.eq.???) then ! no staggering for vertical velocity`
225	`c stag(3)=0.0`	225	`c stag(3)=0.0`
226	`c else`	226	`c else`
227	`c stag(3)=-0.5 ! vertical staggering for all other variables`	227	`c stag(3)=-0.5 ! vertical staggering for all other variables`
228	`c endif`	228	`c endif`
229		229
230	`C Define data region`	230	`C Define data region`
231		231
232	`varmin(1)=lonw`	232	`varmin(1)=lonw`
233	`varmin(2)=lats`	233	`varmin(2)=lats`
234	`varmax(1)=lone`	234	`varmax(1)=lone`
235	`varmax(2)=latn`	235	`varmax(2)=latn`
236	`C varmin/max(3) changed to take into account reverse order of data`	236	`C varmin/max(3) changed to take into account reverse order of data`
237	`varmin(3)=pmin ! Not used`	237	`varmin(3)=pmin ! Not used`
238	`varmax(3)=pmax !50.0`	238	`varmax(3)=pmax !50.0`
239	`write(,) 'lonw,lats,lone,latn',lonw,lats,lone,latn`	239	`write(,) 'lonw,lats,lone,latn',lonw,lats,lone,latn`
240	`write(,) 'pmax, pmin',pmax, pmin`	240	`write(,) 'pmax, pmin',pmax, pmin`
241		241
242	`C Define dimensions`	242	`C Define dimensions`
243		243
244	`vardim(1)=nxi`	244	`vardim(1)=nxi`
245	`vardim(2)=nyi`	245	`vardim(2)=nyi`
246	`CCC Distinction between 2dim and 3dim fields`	246	`CCC Distinction between 2dim and 3dim fields`
247	`CCC Insert num-values for 2dim fields`	247	`CCC Insert num-values for 2dim fields`
248	`if (stash_code.eq.pcode) then`	248	`if (stash_code.eq.pcode) then`
249	`vardim(3)=1`	249	`vardim(3)=1`
250	`else`	250	`else`
251	`vardim(3)=nlevmax`	251	`vardim(3)=nlevmax`
252	`endif`	252	`endif`
253	`vardim(4)=1`	253	`vardim(4)=1`
254		254
255	`C Define the date array`	255	`C Define the date array`
256		256
257	`idate(1)=dyear`	257	`idate(1)=dyear`
258	`idate(2)=dmonth`	258	`idate(2)=dmonth`
259	`idate(3)=dday`	259	`idate(3)=dday`
260	`idate(4)=dhour`	260	`idate(4)=dhour`
261	`idate(5)=dmins`	261	`idate(5)=dmins`
262		262
263	`C set stdate on first read of PP record from file`	263	`C set stdate on first read of PP record from file`
264		264
265	`c if (irec.eq.1) then`	265	`c if (irec.eq.1) then`
266	`c endif`	266	`c endif`
267	`if (wrcstflag.eq.'yes') then`	267	`if (wrcstflag.eq.'yes') then`
268		268
269	`C stdate is the simulation start time`	269	`C stdate is the simulation start time`
270	`stdate(1)=dyear`	270	`stdate(1)=dyear`
271	`stdate(2)=dmonth`	271	`stdate(2)=dmonth`
272	`stdate(3)=dday`	272	`stdate(3)=dday`
273	`stdate(4)=dhour`	273	`stdate(4)=dhour`
274	`stdate(5)=dmins`	274	`stdate(5)=dmins`
275	`C Define the datar array (used to define constants file)`	275	`C Define the datar array (used to define constants file)`
276		276
277	`datar(1)=nxi`	277	`datar(1)=nxi`
278	`datar(2)=nyi`	278	`datar(2)=nyi`
279	`datar(3)=nint(1000.0*latn) ! I agree this is not nice`	279	`datar(3)=nint(1000.0*latn) ! I agree this is not nice`
280	`datar(4)=nint(1000.0*lonw) ! but for historical reasons`	280	`datar(4)=nint(1000.0*lonw) ! but for historical reasons`
281	`datar(5)=nint(1000.0*lats) ! subroutine wricst contains`	281	`datar(5)=nint(1000.0*lats) ! subroutine wricst contains`
282	`datar(6)=nint(1000.0*lone) ! this factor 1000. things`	282	`datar(6)=nint(1000.0*lone) ! this factor 1000. things`
283	`datar(7)=nint(1000.0*dx)`	283	`datar(7)=nint(1000.0*dx)`
284	`datar(8)=nint(1000.0*dy)`	284	`datar(8)=nint(1000.0*dy)`
285	`datar(9)=nlevmax ! just in case Pstar first in file`	285	`datar(9)=nlevmax ! just in case Pstar first in file`
286	`datar(10)=0 ! not used`	286	`datar(10)=0 ! not used`
287	`datar(11)=0 ! not used`	287	`datar(11)=0 ! not used`
288	`datar(12)=0 ! not used`	288	`datar(12)=0 ! not used`
289	`datar(13)=nint(1000.0*pole_lon)`	289	`datar(13)=nint(1000.0*pole_lon)`
290	`datar(14)=nint(1000.0*pole_lat)`	290	`datar(14)=nint(1000.0*pole_lat)`
291		291
292	`C Get the 4 levels/layers arrays`	292	`C Get the 4 levels/layers arrays`
293		293
294	`call modlevs(aklev,bklev,aklay,bklay,plevs,nlevmax,pmin,pdiff)`	294	`call modlevs(aklev,bklev,aklay,bklay,plevs,nlevmax,pmin,pdiff)`
295		295
296	`C Write the constants file`	296	`C Write the constants file`
297		297
298	`call wricst(cfn,datar,aklev,bklev,aklay,bklay,stdate)`	298	`call wricst(cfn,datar,aklev,bklev,aklay,bklay,stdate)`
299	`write(,)`	299	`write(,)`
300	`write(,)'*** constants-file ',cfn(1:strend(cfn)),' created'`	300	`write(,)'*** constants-file ',cfn(1:strend(cfn)),' created'`
301	`wrcstflag='no'`	301	`wrcstflag='no'`
302	`else ! assume that the constants-file already exists`	302	`else ! assume that the constants-file already exists`
303	`C Inquire start time`	303	`C Inquire start time`
304	`call cdfopn(cfn,cstid,ierr) ! open constants file`	304	`call cdfopn(cfn,cstid,ierr) ! open constants file`
305	`call getstart(cstid,stdate,ierr)`	305	`call getstart(cstid,stdate,ierr)`
306	`call ncclos(cstid,ierr) ! close constants file`	306	`call ncclos(cstid,ierr) ! close constants file`
307	`endif`	307	`endif`
308		308
309	`C Create the NetCDF file (only once)`	309	`C Create the NetCDF file (only once)`
310		310
311	`if (.not.opnfile) then`	311	`if (.not.opnfile) then`
312	`c Code to derive Index No. for filenames`	312	`c Code to derive Index No. for filenames`
313	`c if(ntimes.lt.10) then`	313	`c if(ntimes.lt.10) then`
314	`c write(char1,'(i1)') ntimes`	314	`c write(char1,'(i1)') ntimes`
315	`c ctime='0'//char1`	315	`c ctime='0'//char1`
316	`c else`	316	`c else`
317	`c write(ctime,'(i2)') ntimes`	317	`c write(ctime,'(i2)') ntimes`
318	`c endif`	318	`c endif`
319	`c cdf_file=cdfname(1:strend(cdfname))//ctime`	319	`c cdf_file=cdfname(1:strend(cdfname))//ctime`
320	`c Code to derive timestamp for filenames`	320	`c Code to derive timestamp for filenames`
321	`c dhour=dhour-stdate(4) ! subtract start time so always atart at 0`	321	`c dhour=dhour-stdate(4) ! subtract start time so always atart at 0`
322	`if(dhour.eq.0) then`	322	`if(dhour.eq.0) then`
323	`ext1='00'`	323	`ext1='00'`
324	`elseif(dhour.lt.10) then`	324	`elseif(dhour.lt.10) then`
325	`write(char1,'(i1)') dhour`	325	`write(char1,'(i1)') dhour`
326	`ext1='0'//char1`	326	`ext1='0'//char1`
327	`else`	327	`else`
328	`write(ext1,'(i2)') dhour`	328	`write(ext1,'(i2)') dhour`
329	`endif`	329	`endif`
330	`if(dmins.eq.0) then`	330	`if(dmins.eq.0) then`
331	`ext2='00'`	331	`ext2='00'`
332	`else`	332	`else`
333	`write(ext2,'(i2)') dmins`	333	`write(ext2,'(i2)') dmins`
334	`endif`	334	`endif`
335	`c n.b if data runs into next day then can't use root of first`	335	`c n.b if data runs into next day then can't use root of first`
336	`c filename for rest, so form filename manually`	336	`c filename for rest, so form filename manually`
337	`c cdf_file=cdfname(1:strend(cdfname))//ext1//ext2`	337	`c cdf_file=cdfname(1:strend(cdfname))//ext1//ext2`
338	`c below assumes year and month are constant`	338	`c below assumes year and month are constant`
339	`if(dday.lt.10) then`	339	`if(dday.lt.10) then`
340	`write(char1,'(i1)') dday`	340	`write(char1,'(i1)') dday`
341	`date_str='0'//char1`	341	`date_str='0'//char1`
342	`else`	342	`else`
343	`write(date_str,'(i2)') dday`	343	`write(date_str,'(i2)') dday`
344	`endif`	344	`endif`
345	`cdf_file=cdfname(1:5)//date_str//'_'//ext1//ext2`	345	`cdf_file=cdfname(1:5)//date_str//'_'//ext1//ext2`
346	`call crecdf(cdf_file,cdfid,varmin,varmax,`	346	`call crecdf(cdf_file,cdfid,varmin,varmax,`
347	`& 3,cfn(1:strend(cfn)),ierr)`	347	`& 3,cfn(1:strend(cfn)),ierr)`
348	`if (ierr.ne.0) goto 996`	348	`if (ierr.ne.0) goto 996`
349	`write(,)`	349	`write(,)`
350	`c write(,)'*** NetCDF file ',cdfname(1:strend(cdfname)),' create`	350	`c write(,)'*** NetCDF file ',cdfname(1:strend(cdfname)),' create`
351	`write(,)'*** NetCDF file ',cdf_file,' created'`	351	`write(,)'*** NetCDF file ',cdf_file,' created'`
352	`opnfile=.true.`	352	`opnfile=.true.`
353	`endif`	353	`endif`
354		354
355	`C Put 2dim array in 1dim array`	355	`C Put 2dim array in 1dim array`
356		356
357	`do i=1,nxi`	357	`do i=1,nxi`
358	`do j=1,nyi`	358	`do j=1,nyi`
359	`field(i+(j-1)*nxi)=data(i,j)`	359	`field(i+(j-1)*nxi)=data(i,j)`
360	`enddo`	360	`enddo`
361	`enddo`	361	`enddo`
362		362
363	`C Calculate tstep for actual date (relative to first time of`	363	`C Calculate tstep for actual date (relative to first time of`
364	`c simulation) and convert to minutes`	364	`c simulation) and convert to minutes`
365		365
366	`call timediff(idate,stdate,tstep)`	366	`call timediff(idate,stdate,tstep)`
367	`c tstep=tstep*100.0`	367	`c tstep=tstep*100.0`
368	`c tstep=float(nint(tstep*100.0))`	368	`c tstep=float(nint(tstep*100.0))`
369	`write(6,*)"idate,stdate,tstep:",idate,stdate,tstep`	369	`write(6,*)"idate,stdate,tstep:",idate,stdate,tstep`
370	`i=ifix(tstep) ! round down to nearest hour`	370	`i=ifix(tstep) ! round down to nearest hour`
371	`tstep=float(i60+nint(100.0tstep)-100*i)`	371	`tstep=float(i60+nint(100.0tstep)-100*i)`
372	`cwang ilev=1+((1000-lev)/50) this should be changed to /25`	372	`cwang ilev=1+((1000-lev)/50) this should be changed to /25`
373	`cwang change ilev=1+((1000-lev)/50) statement when pdiff read in`	373	`cwang change ilev=1+((1000-lev)/50) statement when pdiff read in`
374	`write(6,*)"tstep:",tstep`	374	`write(6,*)"tstep:",tstep`
375	`ilev=1+((pmin-lev)/pdiff)`	375	`ilev=1+((pmin-lev)/pdiff)`
376	`if(stash_code.eq.pcode) ilev=1`	376	`if(stash_code.eq.pcode) ilev=1`
377		377
378	`if (ilev.eq.1) then`	378	`if (ilev.eq.1) then`
379	`call putdef(cdfid,varname(1:strend(varname)),4,mdv,vardim,`	379	`call putdef(cdfid,varname(1:strend(varname)),4,mdv,vardim,`
380	`& varmin,varmax,stag,ierr)`	380	`& varmin,varmax,stag,ierr)`
381	`if (ierr.ne.0) goto 996`	381	`if (ierr.ne.0) goto 996`
382	`write(,)`	382	`write(,)`
383	`write(,)'*** variable ',varname(1:strend(varname)),' created'`	383	`write(,)'*** variable ',varname(1:strend(varname)),' created'`
384	`nvars=nvars+1`	384	`nvars=nvars+1`
385	`endif`	385	`endif`
386		386
387	`C Put data on file`	387	`C Put data on file`
388		388
389	`c Convert units of Temp to Degrees C`	389	`c Convert units of Temp to Degrees C`
390	`if (stash_code.eq.90002.or.stash_code.eq.90006) then`	390	`if (stash_code.eq.90002.or.stash_code.eq.90006) then`
391	`do i=1,nxi*nyi`	391	`do i=1,nxi*nyi`
392	`field(i)=field(i)-273.15`	392	`field(i)=field(i)-273.15`
393	`enddo`	393	`enddo`
394	`endif`	394	`endif`
395		395
396	`write(,)`	396	`write(,)`
397	`write(,'(a,a,a,f7.1,a,i4,a)')'** var ',`	397	`write(,'(a,a,a,f7.1,a,i4,a)')'** var ',`
398	`& varname(1:strend(varname)),' at time ',`	398	`& varname(1:strend(varname)),' at time ',`
399	`& tstep,' on level ',lev,' written on cdffile'`	399	`& tstep,' on level ',lev,' written on cdffile'`
400	`call putdat(cdfid,varname(1:strend(varname)),tstep,ilev,`	400	`call putdat(cdfid,varname(1:strend(varname)),tstep,ilev,`
401	`& field,ierr)`	401	`& field,ierr)`
402	`if (ierr.ne.0) goto 996`	402	`if (ierr.ne.0) goto 996`
403	`write(,)'cdfid',cdfid`	403	`write(,)'cdfid',cdfid`
404	`c if(nvars.eq.6.and.lev.eq.nlev) then ! required variables written`	404	`c if(nvars.eq.6.and.lev.eq.nlev) then ! required variables written`
405	`write(,)done_u,done_v,done_t,done_q,done_o,done_p`	405	`write(,)done_u,done_v,done_t,done_q,done_o,done_p`
406	`if(done_u+done_v+done_t+done_q+done_o+done_p.eq.6) then`	406	`if(done_u+done_v+done_t+done_q+done_o+done_p.eq.6) then`
407	`opnfile=.false.`	407	`opnfile=.false.`
408	`ntimes=ntimes+1`	408	`ntimes=ntimes+1`
409	`nvars=0`	409	`nvars=0`
410	`C Close NetCDF file`	410	`C Close NetCDF file`
411		411
412	`call ncclos(cdfid,ierr)`	412	`call ncclos(cdfid,ierr)`
413	`if (ierr.ne.0) stop 'error when closing NetCDF file'`	413	`if (ierr.ne.0) stop 'error when closing NetCDF file'`
414	`c if(ntimes.eq.2) stop`	414	`c if(ntimes.eq.2) stop`
415	`done_u=0`	415	`done_u=0`
416	`done_v=0`	416	`done_v=0`
417	`done_t=0`	417	`done_t=0`
418	`done_q=0`	418	`done_q=0`
419	`done_o=0`	419	`done_o=0`
420	`done_p=0`	420	`done_p=0`
421	`endif`	421	`endif`
422	`goto 30`	422	`goto 30`
423		423
424	`996 stop 'error: could not create NetCDF file'`	424	`996 stop 'error: could not create NetCDF file'`
425	`999 continue`	425	`999 continue`
426	`write(,)lev,nlev,ierr`	426	`write(,)lev,nlev,ierr`
427	`if(lev.gt.nlev.or.ierr.ne.0) then`	427	`if(lev.gt.nlev.or.ierr.ne.0) then`
428	`write(,)'=== Error: NetCDF file: ',cdf_file,' incomplete'`	428	`write(,)'=== Error: NetCDF file: ',cdf_file,' incomplete'`
429	`write(,)'See fort.12 for list of data processed'`	429	`write(,)'See fort.12 for list of data processed'`
430	`stop`	430	`stop`
431	`endif`	431	`endif`
432		432
433	`write(,)'=== end of program pp2cdf_p20: status normal ==='`	433	`write(,)'=== end of program pp2cdf_p20: status normal ==='`
434		434
435	`stop`	435	`stop`
436	`end`	436	`end`
437		437
438	`subroutine open_ppf(infile,iform,ieee)`	438	`subroutine open_ppf(infile,iform,ieee)`
439	`C I/O STREAMS ARE AS FOLLOWS`	439	`C I/O STREAMS ARE AS FOLLOWS`
440	`C 1) 6 STANDARD OUTPUT`	440	`C 1) 6 STANDARD OUTPUT`
441	`C 2) 8 PPUNIT - READ UNIT FOR PP FILE`	441	`C 2) 8 PPUNIT - READ UNIT FOR PP FILE`
442	`IMPLICIT NONE`	442	`IMPLICIT NONE`
443	`INTEGER ierr,IFORM,IEEE`	443	`INTEGER ierr,IFORM,IEEE`
444	`CHARACTER INFILE*100`	444	`CHARACTER INFILE*100`
445	`c`	445	`c`
446	`IF(IFORM.EQ.1) THEN`	446	`IF(IFORM.EQ.1) THEN`
447	`OPEN(UNIT=8,FILE=INFILE,FORM='FORMATTED')`	447	`OPEN(UNIT=8,FILE=INFILE,FORM='FORMATTED')`
448	`ELSE`	448	`ELSE`
449	`IF(IEEE.EQ.1) THEN`	449	`IF(IEEE.EQ.1) THEN`
450	`c 2 lines below are Cray specific calls to identify IEEE file`	450	`c 2 lines below are Cray specific calls to identify IEEE file`
451	`c n.b 'newlocal' refers to assign environment, not file`	451	`c n.b 'newlocal' refers to assign environment, not file`
452	`c call asnctl('newlocal',1,ierr)`	452	`c call asnctl('newlocal',1,ierr)`
453	`c write(,) 'ierr from asnctl ',ierr`	453	`c write(,) 'ierr from asnctl ',ierr`
454	`c call asnunit(8,'-F f77 -N ieee',ierr)`	454	`c call asnunit(8,'-F f77 -N ieee',ierr)`
455	`c write(,) 'ierr from asnunit ',ierr`	455	`c write(,) 'ierr from asnunit ',ierr`
456	`ENDIF`	456	`ENDIF`
457	`OPEN(UNIT=8,FILE=INFILE,FORM='UNFORMATTED')`	457	`OPEN(UNIT=8,FILE=INFILE,FORM='UNFORMATTED')`
458	`ENDIF`	458	`ENDIF`
459	`return`	459	`return`
460	`END`	460	`END`
461	`c********************************************************`	461	`c********************************************************`
462	`subroutine read_field(iform,nx,ny,data,stash_code,level,`	462	`subroutine read_field(iform,nx,ny,data,stash_code,level,`
463	`& dxdeg,dydeg,dyear,dmonth,dday,dhour,dmins,iend,`	463	`& dxdeg,dydeg,dyear,dmonth,dday,dhour,dmins,iend,`
464	`& lonw,lone,latn,lats,pole_lon,pole_lat,pres,nxmax,nymax,`	464	`& lonw,lone,latn,lats,pole_lon,pole_lat,pres,nxmax,nymax,`
465	`& NewDyn)`	465	`& NewDyn)`
466	`C I/O STREAMS ARE AS FOLLOWS`	466	`C I/O STREAMS ARE AS FOLLOWS`
467	`C 1) 6 STANDARD OUTPUT`	467	`C 1) 6 STANDARD OUTPUT`
468	`C 2) 8 PPUNIT - READ UNIT FOR PP FILE`	468	`C 2) 8 PPUNIT - READ UNIT FOR PP FILE`
469	`C 3) 12 HEADER INFORMATION STREAM`	469	`C 3) 12 HEADER INFORMATION STREAM`
470	`IMPLICIT NONE`	470	`IMPLICIT NONE`
471	`c`	471	`c`
472	`REAL ROOK(19), FIELD(72900),data(270,270)`	472	`REAL ROOK(19), FIELD(72900),data(270,270)`
473	`INTEGER LOOK(45),I,J,K,N,NUM_VALUES,NX,NY,IFORM,nxmax,nymax,`	473	`INTEGER LOOK(45),I,J,K,N,NUM_VALUES,NX,NY,IFORM,nxmax,nymax,`
474	`& nf,stash_code,level,dyear,dmonth,dday,dhour,dmins,iend`	474	`& nf,stash_code,level,dyear,dmonth,dday,dhour,dmins,iend`
475	`real LONG1,LAT1,dxdeg,dydeg,long_nw,lat_nw,`	475	`real LONG1,LAT1,dxdeg,dydeg,long_nw,lat_nw,`
476	`& fac,lonw,lone,latn,lats,pole_lon,pole_lat,pres`	476	`& fac,lonw,lone,latn,lats,pole_lon,pole_lat,pres`
477	`logical NewDyn`	477	`logical NewDyn`
478	`data nf/0/`	478	`data nf/0/`
479	`save nf`	479	`save nf`
480	`nf=nf+1`	480	`nf=nf+1`
481	`c`	481	`c`
482	`IF(IFORM.EQ.1) THEN`	482	`IF(IFORM.EQ.1) THEN`
483	`READ(8,12,END=99) LOOK`	483	`READ(8,12,END=99) LOOK`
484	`READ(8,13) ROOK`	484	`READ(8,13) ROOK`
485	`ELSE`	485	`ELSE`
486	`READ(8,END=99) LOOK,ROOK`	486	`READ(8,END=99) LOOK,ROOK`
487	`ENDIF`	487	`ENDIF`
488	`NX=LOOK(19)`	488	`NX=LOOK(19)`
489	`NY=LOOK(18)`	489	`NY=LOOK(18)`
490	`NUM_VALUES=NX*NY`	490	`NUM_VALUES=NX*NY`
491	`C WRITE OUT CONTENTS OF INFILE FOR INFORMATION`	491	`C WRITE OUT CONTENTS OF INFILE FOR INFORMATION`
492	`WRITE(12,106) LOOK(33),LOOK(14),nf,`	492	`WRITE(12,106) LOOK(33),LOOK(14),nf,`
493	`& NUM_VALUES,LOOK(42),LOOK(18),LOOK(19)`	493	`& NUM_VALUES,LOOK(42),LOOK(18),LOOK(19)`
494	`IF(IFORM.EQ.1) THEN`	494	`IF(IFORM.EQ.1) THEN`
495	`READ(8,13) (FIELD(I),I=1,NUM_VALUES)`	495	`READ(8,13) (FIELD(I),I=1,NUM_VALUES)`
496	`ELSE`	496	`ELSE`
497	`write(,) '1'`	497	`write(,) '1'`
498	`READ(8) (FIELD(I),I=1,NUM_VALUES)`	498	`READ(8) (FIELD(I),I=1,NUM_VALUES)`
499	`ENDIF`	499	`ENDIF`
500	`write(,) '2'`	500	`write(,) '2'`

Subversion Repositories lagranto.um

(root)/tags/2.2/convert/pp2cdf_p20_Andy_Elvidge.f – Rev 3 → 14