Subversion Repositories lagranto.um

      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'
CCC   Include UM dimensions(nxmax,nymax,nzmax,ntmax)
      include 'um_dims.inc'
C     ===  definitions of routinespecific variables  =======================

CCC   maximum array dimensions (change them for your purpose)
      integer   nxi,nyi,nlevmax,nwet
      real      pmin,pmax,pdiff
      parameter(nlevmax=nzmax)
      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*200 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, pdiff39C     Define constants file name

      if (cfn.eq.'default') then
        cfn=trim(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
C NH: changed to modplevs from modlevs
 
        call modplevs(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 ',trim(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=trim(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.lt.10) then
          write(char1,'(i1)') dmins 
          ext2='0'//char1
        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=trim(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
cSue modify to get expected structure for new trajectory code
        cdf_file=cdfname(1:7)//date_str//'_'//ext1//ext2
        call crecdf(cdf_file,cdfid,varmin,varmax,
     &              3,trim(cfn),ierr)
        if (ierr.ne.0) goto 996
        write(*,*)
c       write(*,*)'*** NetCDF file ',trim(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
cSue - leave in hours for new traj code  
c         i=ifix(tstep)  ! round down to nearest hour
c         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,trim(varname),4,mdv,vardim,
     &              varmin,varmax,stag,ierr)
        if (ierr.ne.0) goto 996
        write(*,*)
        write(*,*)'*** variable ',trim(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 ',
     &  trim(varname),' at time ',
     &  tstep,' on level ',lev,' written on cdffile'
      call putdat(cdfid,trim(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*120
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
      write(*,*) INFILE
      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(nxmax*nymax),data(nxmax,nymax)
      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 (.not. NewDyn) 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 (.not. NewDyn) 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 modplevs(aklev,bklev,aklay,bklay,plevs,nlevmax,pm,pd)
C NH: changed to modplevs from modlevs
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