Ticket #480: und_dat2nc.f90

File und_dat2nc.f90, 5.0 KB (added by Ian Culverwell, 7 years ago)

und_dat2nc.f90

Line 
1!export LIBDIR=/project/ukmo/rhel6/netcdf4/ifort_composerxe
2!ifort -o und_dat2nc.exe und_dat2nc.f90 -I$LIBDIR/include -L$LIBDIR/lib -lnetcdf -lnetcdff
3
4PROGRAM und_dat2nc
5
6!
7! Convert ropp_pp/data/egm96.dat and ropp_pp/data/corrcoef.dat into a
8! single netcdf file, which need not be zipped.
9!
10! IDC, 23/12/2016
11!
12
13 USE netcdf
14
15 IMPLICIT NONE
16
17 INTEGER, PARAMETER :: wp=KIND(1.D0)
18
19 INTEGER, PARAMETER :: ll=65341 ! = 361 * 362 / 2, i.e. for a 1/4 deg grid
20 ! (cos(i*360*lon) ~ (1, 0, -1, 0) per degree of lon)
21
22 CHARACTER(nf90_max_name) :: ifile1='egm96.dat'
23 INTEGER, PARAMETER :: lun1=11
24
25 CHARACTER(nf90_max_name) :: ifile2='corrcoef.dat'
26 INTEGER, PARAMETER :: lun2=12
27
28 REAL(wp), DIMENSION(ll) :: hc, hs, cc, cs, hc2, cc2
29 CHARACTER(nf90_max_name) :: ofile='und_coeffs.nc'
30
31 INTEGER :: i, n, m, ierr
32 INTEGER :: ncid, ll_dimid
33 INTEGER :: hc_varid, hs_varid, cs_varid, cc_varid
34 REAL(wp) :: c, s
35
36! 1. Input undulation coefficients
37! --------------------------------
38
39 OPEN ( UNIT=lun1, STATUS="OLD", FILE=ifile1, ACTION="READ" )
40
41 hc = 0.0_wp ; hs = 0.0_wp
42
43 DO
44 READ ( UNIT=lun1, FMT=*, IOSTAT=ierr ) n, m, c, s
45 IF ( ierr /= 0 ) EXIT
46 i = ( n * ( n + 1 ) ) / 2 + m + 1
47 hc(i) = c
48 hs(i) = s
49 END DO
50
51 CLOSE ( UNIT=lun1 )
52
53 WRITE(*, '(/,A,/,(5E30.20))') 'From ' // TRIM(ADJUSTL(ifile1)) // ': hc(1:10) = ', hc(1:10)
54
55! 2. Input correction coefficients
56! --------------------------------
57
58 OPEN ( UNIT=lun2, STATUS="OLD", FILE=ifile2, ACTION="READ" )
59
60 cc = 0.0_wp ; cs = 0.0_wp
61
62 DO
63 READ ( UNIT=lun2, FMT=*, IOSTAT=ierr ) n, m, c, s
64 IF ( ierr /= 0 ) EXIT
65 i = ( n * ( n + 1 ) ) / 2 + m + 1
66 cc(i) = c
67 cs(i) = s
68 END DO
69
70 CLOSE ( UNIT=lun2 )
71
72 WRITE(*, '(/,A,/,(5E30.20))') 'From ' // TRIM(ADJUSTL(ifile2)) // ': cc(1:10) = ', cc(1:10)
73
74! 3. Output both
75! --------------
76
77! Create output file
78 CALL check(nf90_create(ofile, nf90_clobber, ncid))
79
80! Define the dimensions
81 CALL check(nf90_def_dim(ncid, 'n', ll, ll_dimid))
82
83! Define the variables
84 CALL check(nf90_def_var(ncid, 'hc', nf90_double, ll_dimid, hc_varid))
85 CALL check(nf90_def_var(ncid, 'hs', nf90_double, ll_dimid, hs_varid))
86 CALL check(nf90_def_var(ncid, 'cc', nf90_double, ll_dimid, cc_varid))
87 CALL check(nf90_def_var(ncid, 'cs', nf90_double, ll_dimid, cs_varid))
88
89! Define attributes
90 CALL check(nf90_put_att(ncid, hc_varid, 'long_name', 'First (cosine) undulation coefficient'))
91 CALL check(nf90_put_att(ncid, hc_varid, 'units', '1'))
92 CALL check(nf90_put_att(ncid, hs_varid, 'long_name', 'Second (sine) undulation coefficient'))
93 CALL check(nf90_put_att(ncid, hs_varid, 'units', '1'))
94 CALL check(nf90_put_att(ncid, cc_varid, 'long_name', 'First (cosine) correction coefficient'))
95 CALL check(nf90_put_att(ncid, cc_varid, 'units', '1'))
96 CALL check(nf90_put_att(ncid, cs_varid, 'long_name', 'Second (sine) correction coefficient'))
97 CALL check(nf90_put_att(ncid, cs_varid, 'units', '1'))
98 CALL check(nf90_put_att(ncid, nf90_global, 'Source', 'Derived from WGS 84 EGM96 15-Minute geoid height and correction coefficients'))
99 CALL check(nf90_put_att(ncid, nf90_global, 'Reference', 'http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html'))
100 CALL check(nf90_put_att(ncid, nf90_global, 'Note', 'In early versions of ROPP, these data were held in ropp_pp/data/egm96.dat and ropp_pp/data/corrcoef.dat'))
101 CALL check(nf90_put_att(ncid, nf90_global, 'Author', 'ROM SAF, 2017'))
102
103! Come out of definition mode
104 CALL check(nf90_enddef(ncid))
105
106! Write out variables
107 CALL check(nf90_put_var(ncid, hc_varid, hc))
108 CALL check(nf90_put_var(ncid, hs_varid, hs))
109 CALL check(nf90_put_var(ncid, cc_varid, cc))
110 CALL check(nf90_put_var(ncid, cs_varid, cs))
111
112! Close the output file
113 CALL check(nf90_close(ncid))
114
115! Check that we can read it
116 hc2 = 0.0_wp
117 CALL check(nf90_open(ofile, nf90_nowrite, ncid))
118 CALL check(nf90_inq_varid(ncid, 'hc', hc_varid))
119 CALL check(nf90_get_var(ncid, hc_varid, hc2))
120 CALL check(nf90_close(ncid))
121
122 WRITE(*, '(/,A,/,(5E30.20))') 'From ' // TRIM(ADJUSTL(ofile)) // ': hc(1:10) = ', hc2(1:10)
123
124 WRITE(*, '(/,A,/,(5E30.20))') TRIM(ADJUSTL(ifile1)) // ' - ' // TRIM(ADJUSTL(ofile)) // ': hc(1:10) = ', hc(1:10)-hc2(1:10)
125
126 cc2 = 0.0_wp
127 CALL check(nf90_open(ofile, nf90_nowrite, ncid))
128 CALL check(nf90_inq_varid(ncid, 'cc', cc_varid))
129 CALL check(nf90_get_var(ncid, cc_varid, cc2))
130 CALL check(nf90_close(ncid))
131
132 WRITE(*, '(/,A,/,(5E30.20))') 'From ' // TRIM(ADJUSTL(ofile)) // ': cc(1:10) = ', cc2(1:10)
133
134 WRITE(*, '(/,A,/,(5E30.20))') TRIM(ADJUSTL(ifile1)) // ' - ' // TRIM(ADJUSTL(ofile)) // ': cc(1:10) = ', cc(1:10)-cc2(1:10)
135
136CONTAINS
137
138
139!----------------------------------------------------------------------------
140
141 SUBROUTINE check(status)
142
143 INTEGER, INTENT (IN) :: status
144
145 IF (status /= nf90_noerr) THEN
146 PRINT*, TRIM(nf90_strerror(status))
147 STOP "STOPPED"
148 END IF
149
150 END SUBROUTINE check
151
152
153END PROGRAM und_dat2nc