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