Commit f60b5c7b authored by Burkhardt Rockel's avatar Burkhardt Rockel

Repository starts with cfu_v1.5

parents
#########
# Intel Fortran at DKRZ Mistral
#########
#
#... netCDF
NETCDFC_ROOT = /sw/rhel6-x64/netcdf/netcdf_c-4.3.2-gcc48
NETCDFF_ROOT = /sw/rhel6-x64/netcdf/netcdf_fortran-4.4.2-intel14
NC_INC = -I${NETCDFF_ROOT}/include
NC_LIB = -L${NETCDFF_ROOT}/lib -lnetcdff
NC_LIB += -L${NETCDFC_ROOT}/lib -lnetcdf
# standard binary
PROGRAM = bin/cfu
# compiler, options and libraries
F90 = ifort
#F90 = gfortran
COMFLG = -c ${NC_INC}
LDSEQ = ifort
#LDSEQ = gfortran
LDFLG =
LIB = ${NC_LIB}
LIB += -Wl,-rpath,$(NETCDFF_ROOT)/lib:$(NETCDFC_ROOT)/lib
#########
# OSX gfortran
#########
#
# standard binary
PROGRAM = bin/cfu
# compiler, options and libraries
F90 = gfortran
COMFLG = -c
LDSEQ = gfortran
LDFLG =
#########
# Intel Fortran at DKRZ Mistral
#########
#
#... netCDF
NETCDFC_ROOT = /sw/rhel6-x64/netcdf/netcdf_c-4.3.2-gcc48
NETCDFF_ROOT = /sw/rhel6-x64/netcdf/netcdf_fortran-4.4.2-intel14
NC_INC = -I${NETCDFF_ROOT}/include
NC_LIB = -L${NETCDFF_ROOT}/lib -lnetcdff
NC_LIB += -L${NETCDFC_ROOT}/lib -lnetcdf
# standard binary
PROGRAM = bin/cfu
# compiler, options and libraries
F90 = ifort
#F90 = gfortran
COMFLG = -c ${NC_INC}
LDSEQ = ifort
#LDSEQ = gfortran
LDFLG =
LIB = ${NC_LIB}
LIB += -Wl,-rpath,$(NETCDFF_ROOT)/lib:$(NETCDFC_ROOT)/lib
####################################################
#
# Makefile for the Fortran 90 program
#
####################################################
#
#
####################################################
#
# Declaration of Variables
#
####################################################
#
.SILENT:
#
STDROOT = $(PWD)
SRCDIR = $(STDROOT)/src
OBJDIR = $(STDROOT)/obj
#
#########
# include file containing system special configurations
#########
#
include Fopts
#
#
####################################################
#
# Declaration of the Object Files
#
####################################################
#
include ObjFiles
#
#
####################################################
#
# valid actions
#
####################################################
#
#
exe: $(OBJECTS)
echo "creating cfu executable"
$(LDSEQ) -o $(PROGRAM) $(OBJECTS) $(LIB)
#
clean:
echo cleaning up
rm -f $(PROGRAM)
rm -f $(OBJDIR)/*.o
rm -f $(OBJDIR)/*.mod
#
####################################################
#
# Dependencies of the Data Modules
#
####################################################
#
include ObjDependencies
#
####################################################
#
# Dependencies of the Data Modules
#
####################################################
#
$(OBJDIR)/add_hours.o: $(SRCDIR)/add_hours.f90
( cd $(OBJDIR) && echo Compile ${@F} && $(F90) $(COMFLG) $(SRCDIR)/add_hours.f90 )
$(OBJDIR)/add_lonlat.o: $(SRCDIR)/add_lonlat.f90 \
$(OBJDIR)/utilities.o
( cd $(OBJDIR) && echo Compile ${@F} && $(F90) $(COMFLG) $(SRCDIR)/add_lonlat.f90 )
$(OBJDIR)/add_months.o: $(SRCDIR)/add_months.f90
( cd $(OBJDIR) && echo Compile ${@F} && $(F90) $(COMFLG) $(SRCDIR)/add_months.f90 )
$(OBJDIR)/add_vertices.o: $(SRCDIR)/add_vertices.f90 \
$(OBJDIR)/utilities.o
( cd $(OBJDIR) && echo Compile ${@F} && $(F90) $(COMFLG) $(SRCDIR)/add_vertices.f90 )
$(OBJDIR)/check_files.o: $(SRCDIR)/check_files.f90 \
$(OBJDIR)/utilities.o
( cd $(OBJDIR) && echo Compile ${@F} && $(F90) $(COMFLG) $(SRCDIR)/check_files.f90 )
$(OBJDIR)/check_files_iso.o: $(SRCDIR)/check_files_iso.f90 \
$(OBJDIR)/utilities.o
( cd $(OBJDIR) && echo Compile ${@F} && $(F90) $(COMFLG) $(SRCDIR)/check_files_iso.f90 )
$(OBJDIR)/get_hours.o: $(SRCDIR)/get_hours.f90
( cd $(OBJDIR) && echo Compile ${@F} && $(F90) $(COMFLG) $(SRCDIR)/get_hours.f90 )
$(OBJDIR)/get_next_dates.o: $(SRCDIR)/get_next_dates.f90
( cd $(OBJDIR) && echo Compile ${@F} && $(F90) $(COMFLG) $(SRCDIR)/get_next_dates.f90 )
$(OBJDIR)/get_prev_date.o: $(SRCDIR)/get_prev_date.f90
( cd $(OBJDIR) && echo Compile ${@F} && $(F90) $(COMFLG) $(SRCDIR)/get_prev_date.f90 )
$(OBJDIR)/utilities.o: $(SRCDIR)/utilities.f90
( cd $(OBJDIR) && echo Compile ${@F} && $(F90) $(COMFLG) $(SRCDIR)/utilities.f90 )
$(OBJDIR)/shellpp.o: $(SRCDIR)/shellpp.f90
( cd $(OBJDIR) && echo Compile ${@F} && $(F90) $(COMFLG) $(SRCDIR)/shellpp.f90 )
$(OBJDIR)/cfu.o: $(SRCDIR)/cfu.f90 \
$(OBJDIR)/add_hours.o $(OBJDIR)/add_months.o \
$(OBJDIR)/check_files.o $(OBJDIR)/check_files_iso.o \
$(OBJDIR)/get_hours.o \
$(OBJDIR)/get_next_dates.o \
$(OBJDIR)/get_prev_date.o \
$(OBJDIR)/add_lonlat.o $(OBJDIR)/utilities.o \
$(OBJDIR)/shellpp.o
( cd $(OBJDIR) && echo Compile ${@F} && $(F90) $(COMFLG) $(SRCDIR)/cfu.f90 )
#
####################################################
#
# Declaration of the Object Files
#
####################################################
#
OBJECTS = $(OBJDIR)/add_hours.o \
$(OBJDIR)/add_lonlat.o \
$(OBJDIR)/add_months.o \
$(OBJDIR)/add_vertices.o \
$(OBJDIR)/check_files.o \
$(OBJDIR)/check_files_iso.o \
$(OBJDIR)/get_hours.o \
$(OBJDIR)/get_next_dates.o \
$(OBJDIR)/get_prev_date.o \
$(OBJDIR)/utilities.o \
$(OBJDIR)/shellpp.o \
$(OBJDIR)/cfu.o
SUBROUTINE add_hours (ystartdate, hadd, itype_calendar)
!-------------------------------------------------------------------------------
!
! Description:
! This program adds a number of hours to a date
! YYYYMMDDHH[MMSS](new) = YYYMMDDHH[MMSS](old) + hours
!
! Usage:
! add_hours <yyyymmddhh[mmss]> <hours> <itype_calendar>
!
! The result is written to standard output in the form YYYYMMDDHH
!
!-------------------------------------------------------------------------------
!
IMPLICIT NONE
! Input Parameter list:
! ---------------------
CHARACTER (LEN=14) , INTENT(IN) :: &
ystartdate ! start date of the forecast
REAL, INTENT(IN) :: &
hadd ! hours to be added to ystartdate
INTEGER, INTENT(IN) :: &
itype_calendar ! type of calendar
! 0 = standard calendar
! 1 = climatological calendar (360 days per year)
! 2 = 365 days per year
! Parameter written to standard output
! ------------------------------------
CHARACTER (LEN=14) :: &
yactdate1 ! actual date in the form yyyymmddhh
! Local variables:
! ------------------------------------
INTEGER :: &
month(12), monthsum(13), ileap, iweek, iy, m, &
idd, imm, iyy, ihh, iday, imonth, iyear, ihour, immhours, iyyhours
CHARACTER (LEN=3) :: yweek(7)
CHARACTER (LEN=20) :: chadd
!------------ End of header ----------------------------------------------------
! Begin subroutine add_hours
DATA month / 31 , 28 , 31 , 30 , 31 , 30 , &
31 , 31 , 30 , 31 , 30 , 31 /
DATA yweek /'MON', 'TUE', 'WED', 'THU', 'FRI', 'SAT', 'SUN' /
! Statementfunction: ileap(yy) = 0: no leap year,
! ileap(yy) = 1: leap year
ileap (iy) = IABS( MOD(iy, 4) - 4) / 4 & ! every 4 years is a leapyear
-IABS( MOD(iy,100) - 100) / 100 & ! every 100 years is no leapyear
+IABS( MOD(iy,400) - 400) / 400 ! but every 400 years is a leapyear
! Divide ystartdate in day, month, year and hour
! and calculate the sums of days from the beginning of the year to the
! end of the months
READ ( ystartdate, '(I4,3I2)' ) iyy, imm, idd, ihh
SELECT CASE (itype_calendar)
CASE (0) ! Standard year
month (2) = 28 + ileap (iyy)
monthsum(1) = 0
DO m = 2 , 13
monthsum(m) = monthsum(m-1) + month(m-1)
enddo
! Determine how many hours have passed in this year
iyyhours = (idd*24) + monthsum(imm)*24 + (ihh-24)
iyyhours = iyyhours + hadd
! Take turning of the year into account
IF (iyyhours < 0) THEN
iyear = iyy
DO WHILE (iyyhours >= (8760+ileap(iyear)*24))
iyyhours = iyyhours - (8760+ileap(iyear)*24)
iyear=iyear+1
month (2) = 28 + ileap (iyear)
monthsum(1) = 0
DO m = 2 , 13
monthsum(m) = monthsum(m-1) + month(m-1)
ENDDO
ENDDO
ELSE IF (iyyhours >= (8760+ileap(iyy)*24)) THEN
! Take also into account if the run lasts
! for several years
iyear = iyy
DO WHILE (iyyhours >= (8760+ileap(iyear)*24))
iyyhours = iyyhours - (8760+ileap(iyear)*24)
iyear=iyear+1
month (2) = 28 + ileap (iyear)
monthsum(1) = 0
DO m = 2 , 13
monthsum(m) = monthsum(m-1) + month(m-1)
ENDDO
ENDDO
ELSE
iyear = iyy
ENDIF
! Determine the actual date from iyyhours
m = 1
immhours = iyyhours
DO WHILE (immhours >= 0)
m = m+1
immhours = iyyhours - monthsum(m) * 24
ENDDO
imonth = m-1
immhours = iyyhours - monthsum(imonth)*24
iday = immhours/24 + 1
ihour = MOD(immhours,24)
iweek = MOD(monthsum(imonth) + iday + (iyear-1981)+(iyear-1981)/4+2 , 7) + 1
WRITE ( yactdate1(1:4) , '(I4.4)' ) iyear
WRITE ( yactdate1(5:6) , '(I2.2)' ) imonth
WRITE ( yactdate1(7:8) , '(I2.2)' ) iday
WRITE ( yactdate1(9:10), '(I2.2)' ) ihour
IF (LEN(TRIM(ystartdate)) == 14) THEN
yactdate1(11:14) = ystartdate(11:14)
write(*,'(A14)') yactdate1
ELSE
write(*,'(A10)') yactdate1(1:10)
ENDIF
!....
CASE (1) ! 360 days per year
monthsum(1) = 0
DO m = 2 , 13
monthsum(m) = monthsum(m-1) + 30
enddo
! Determine how many hours have passed in this year
iyyhours = (idd*24) + monthsum(imm)*24 + (ihh-24)
iyyhours = iyyhours + hadd
! Take turning of the year into account
IF (iyyhours < 0) THEN
iyear = iyy-1
iyyhours = 8640 + iyyhours
ELSE IF (iyyhours >= 8640) THEN
iyear = iyy
DO WHILE (iyyhours >= 8640)
iyear = iyear+1
iyyhours = iyyhours - 8640
ENDDO
ELSE
iyear = iyy
ENDIF
! Determine the actual date from iyyhours
m = 1
immhours = iyyhours
DO WHILE (immhours >= 0)
m = m+1
immhours = iyyhours - monthsum(m) * 24
ENDDO
imonth = m-1
immhours = iyyhours - monthsum(imonth)*24
iday = immhours/24 + 1
ihour = MOD(immhours,24)
iweek = MOD(monthsum(imonth) + iday + (iyear-1981) , 7) + 1
WRITE ( yactdate1(1:4) , '(I4.4)' ) iyear
WRITE ( yactdate1(5:6) , '(I2.2)' ) imonth
WRITE ( yactdate1(7:8) , '(I2.2)' ) iday
WRITE ( yactdate1(9:10), '(I2.2)' ) ihour
IF (LEN(TRIM(ystartdate)) == 14) THEN
yactdate1(11:14) = ystartdate(11:14)
write(*,'(A14)') yactdate1
ELSE
write(*,'(A10)') yactdate1(1:10)
ENDIF
!....
CASE (2) ! 365 days per year
monthsum(1) = 0
DO m = 2 , 13
monthsum(m) = monthsum(m-1) + month(m-1)
enddo
! Determine how many hours have passed in this year
iyyhours = (idd*24) + monthsum(imm)*24 + (ihh-24)
iyyhours = iyyhours + hadd
! Take turning of the year into account
IF (iyyhours < 0) THEN
iyear = iyy
DO WHILE (iyyhours >= 8760)
iyyhours = iyyhours - 8760
iyear=iyear+1
monthsum(1) = 0
DO m = 2 , 13
monthsum(m) = monthsum(m-1) + month(m-1)
ENDDO
ENDDO
ELSE IF (iyyhours >= 8760) THEN
! Take also into account if the run lasts
! for several years
iyear = iyy
DO WHILE (iyyhours >= 8760)
iyyhours = iyyhours - 8760
iyear=iyear+1
monthsum(1) = 0
DO m = 2 , 13
monthsum(m) = monthsum(m-1) + month(m-1)
ENDDO
ENDDO
ELSE
iyear = iyy
ENDIF
! Determine the actual date from iyyhours
m = 1
immhours = iyyhours
DO WHILE (immhours >= 0)
m = m+1
immhours = iyyhours - monthsum(m) * 24
ENDDO
imonth = m-1
immhours = iyyhours - monthsum(imonth)*24
iday = immhours/24 + 1
ihour = MOD(immhours,24)
iweek = MOD(monthsum(imonth) + iday + (iyear-1981)+(iyear-1981)/4+2 , 7) + 1
WRITE ( yactdate1(1:4) , '(I4.4)' ) iyear
WRITE ( yactdate1(5:6) , '(I2.2)' ) imonth
WRITE ( yactdate1(7:8) , '(I2.2)' ) iday
WRITE ( yactdate1(9:10), '(I2.2)' ) ihour
IF (LEN(TRIM(ystartdate)) == 14) THEN
yactdate1(11:14) = ystartdate(11:14)
write(*,'(A14)') yactdate1
ELSE
write(*,'(A10)') yactdate1(1:10)
ENDIF
END SELECT
END SUBROUTINE add_hours
SUBROUTINE add_lonlat (ncfile)
!-------------------------------------------------------------------------------
!
! Description:
! This program adds 2D geographical latitude and longitude to a netCDF file
! containing rotated coordinats
!
! Usage:
! add_lonlat <ncfile>
!
!-------------------------------------------------------------------------------
!
USE netcdf
USE utilities
IMPLICIT NONE
! Input Parameter list:
! ---------------------
CHARACTER (LEN=*) :: ncfile
! Local variables:
! ------------------------------------
INTEGER, PARAMETER :: &
FourByteReal = selected_real_kind(P = 6, R = 37), &
EightByteReal = selected_real_kind(P = 13, R = 307)
REAL :: pollon, pollat, polgam
CHARACTER (LEN=80) :: &
grid_name ! name of grid mapping
INTEGER :: i, j, istatus
INTEGER :: ncid
INTEGER :: varid_rlon, varid_rlat, dimid_rlon, dimid_rlat
INTEGER :: varid_lon, varid_lat, varid_rotated_pole
INTEGER :: nlon, nlat
REAL (kind=EightByteReal), ALLOCATABLE :: &
rlon(:), &
rlat(:)
REAL (kind=FourByteReal), ALLOCATABLE :: &
lon(:,:), &
lat(:,:)
!------------ End of header ----------------------------------------------------
! get file ID
istatus = nf90_open (ncfile, NF90_WRITE, ncid)
IF (istatus /= NF90_NOERR) THEN
PRINT *, 'Error in nf90_open'
PRINT *, TRIM(NF90_strerror(istatus))
STOP
ENDIF
istatus = nf90_inq_varid(ncid, 'rotated_pole', varid_rotated_pole)
IF (istatus == NF90_NOERR) THEN ! rotated latitude-longitude
istatus = nf90_get_att(ncid, varid_rotated_pole, 'grid_mapping_name', grid_name)
IF (istatus /= NF90_NOERR) THEN
PRINT *, 'Error in read_nc_gdefs / nf90_get_att'
PRINT *, 'Attribute "grid_mapping_name"'
PRINT *, TRIM(NF90_strerror(istatus))
ENDIF
IF (grid_name(1:26) /= 'rotated_latitude_longitude') THEN
PRINT *, 'Error in read_nc_gdefs'
PRINT *, 'Invalid value for attribute "grid_mapping_name"'
PRINT *, TRIM(grid_name)
ENDIF
istatus = nf90_get_att(ncid, varid_rotated_pole, 'grid_north_pole_latitude', pollat)
IF (istatus /= NF90_NOERR) THEN
PRINT *, 'Error in read_nc_gdefs / nf90_get_att'
PRINT *, 'Attribute "grid_north_pole_latitude"'
PRINT *, TRIM(NF90_strerror(istatus))
RETURN
ENDIF
istatus = nf90_get_att(ncid, varid_rotated_pole, 'grid_north_pole_longitude', pollon)
IF (istatus /= NF90_NOERR) THEN
PRINT *, 'Error in read_nc_gdefs / nf90_get_att'
PRINT *, 'Attribute "grid_north_pole_longitude"'
PRINT *, TRIM(NF90_strerror(istatus))
RETURN
ENDIF
istatus = nf90_get_att(ncid, varid_rotated_pole, 'north_pole_grid_longitude', polgam)
IF (istatus == NF90_ENOTATT) THEN
polgam = 0.
ELSE IF (istatus /= NF90_NOERR) THEN
PRINT *, 'Error in read_nc_gdefs / nf90_get_att'
PRINT *, 'Attribute "north_pole_grid_longitude"'
PRINT *, TRIM(NF90_strerror(istatus))
RETURN
ENDIF
ELSE IF (istatus == NF90_ENOTVAR) THEN ! true latitude-longitude
pollon = 0.
pollat = 90.
polgam = 0.
ELSE
PRINT *, 'Error in read_nc_gdefs / nf90_inq_varid'
PRINT *, 'Variable "rotated_pole"'
RETURN
ENDIF
istatus = nf90_inq_varid(ncid, 'rlon', varid_rlon)
IF (istatus /= NF90_NOERR) THEN
PRINT *, 'Error in nccd_get / nf90_inq_varid'
PRINT *, TRIM(NF90_strerror(istatus))
PRINT *, 'rlon'
STOP
ENDIF
istatus = nf90_inq_dimid(ncid, 'rlon', dimid_rlon)
IF (istatus /= NF90_NOERR) THEN
PRINT 'Error in rlon / nf90_inq_dimid -- '// &
TRIM(NF90_strerror(istatus))
STOP
ENDIF
istatus = nf90_inquire_dimension(ncid, dimid_rlon, len=nlon)
IF (istatus /= NF90_NOERR) THEN
PRINT *, 'Error in rlon / nf90_inquire_dimension -- '// &
TRIM(NF90_strerror(istatus))
STOP
ENDIF
istatus = nf90_inq_varid(ncid, 'rlat', varid_rlat)
IF (istatus /= NF90_NOERR) THEN
PRINT *, 'Error in nccd_get / nf90_inq_varid'
PRINT *, TRIM(NF90_strerror(istatus))
PRINT *, 'rlat'
STOP
ENDIF
istatus = nf90_inq_dimid(ncid, 'rlat', dimid_rlat)
IF (istatus /= NF90_NOERR) THEN
PRINT *, 'Error in rlat / nf90_inq_dimid -- '// &
TRIM(NF90_strerror(istatus))
STOP
ENDIF
istatus = nf90_inquire_dimension(ncid, dimid_rlat, len=nlat)
IF (istatus /= NF90_NOERR) THEN
PRINT *, 'Error in rlon / nf90_inquire_dimension -- '// &
TRIM(NF90_strerror(istatus))
STOP
ENDIF
ALLOCATE (rlon(nlon))
ALLOCATE (rlat(nlat))
ALLOCATE (lon(nlon, nlat))
ALLOCATE (lat(nlon, nlat))
! read the data values
istatus = nf90_get_var (ncid, varid_rlon, rlon)
IF (istatus /= NF90_NOERR) THEN
PRINT *, 'Error in rlon nf90_get_var'
PRINT *, TRIM(NF90_strerror(istatus))
STOP
ENDIF
istatus = nf90_get_var (ncid, varid_rlat, rlat)
IF (istatus /= NF90_NOERR) THEN
PRINT *, 'Error in rlon nf90_get_var'
PRINT *, TRIM(NF90_strerror(istatus))
STOP
ENDIF
istatus = nf90_redef(ncid)
! geographic longitude
istatus=nf90_def_var(ncid, "lon", NF90_FLOAT, (/ dimid_rlon, dimid_rlat /), varid_lon)
IF (istatus /= NF90_NOERR) THEN
PRINT *,TRIM(NF90_strerror(istatus))
STOP
ENDIF
istatus=nf90_put_att(ncid, varid_lon, "standard_name", "longitude")
IF (istatus /= NF90_NOERR) THEN
PRINT *,TRIM(NF90_strerror(istatus))
STOP
ENDIF
istatus=nf90_put_att(ncid, varid_lon, "long_name", "longitude")
IF (istatus /= NF90_NOERR) THEN
PRINT *,TRIM(NF90_strerror(istatus))
STOP
ENDIF
istatus=nf90_put_att(ncid, varid_lon, "units", "degrees_east")
IF (istatus /= NF90_NOERR) THEN
PRINT *,TRIM(NF90_strerror(istatus))
STOP
ENDIF
! geographic latitude
istatus=nf90_def_var(ncid, "lat", NF90_FLOAT, (/ dimid_rlon, dimid_rlat /), varid_lat)
IF (istatus /= NF90_NOERR) THEN
PRINT *,TRIM(NF90_strerror(istatus))
STOP
ENDIF
istatus=nf90_put_att(ncid, varid_lat, "standard_name", "latitude")
IF (istatus /= NF90_NOERR) THEN
PRINT *,TRIM(NF90_strerror(istatus))
STOP
ENDIF
istatus=nf90_put_att(ncid, varid_lat, "long_name", "latitude")
IF (istatus /= NF90_NOERR) THEN
PRINT *,TRIM(NF90_strerror(istatus))
STOP
ENDIF
istatus=nf90_put_att(ncid, varid_lat, "units", "degrees_north")
IF (istatus /= NF90_NOERR) THEN
PRINT *,TRIM(NF90_strerror(istatus))
STOP
ENDIF
istatus = nf90_enddef(ncid)
DO j = 1, nlat
DO i = 1, nlon
lon(i,j) = rlarot2rla (REAL(rlat(j)), REAL(rlon(i)), pollat, pollon, polgam)
lat(i,j) = phirot2phi (REAL(rlat(j)), REAL(rlon(i)), pollat, pollon, polgam)
ENDDO
ENDDO