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 |
|
---|
4 | PROGRAM 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 |
|
---|
136 | CONTAINS
|
---|
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 |
|
---|
153 | END PROGRAM und_dat2nc
|
---|