[1158] | 1 | ! $Id$ |
---|
| 2 | module regr_lat_time_climoz_m |
---|
| 3 | |
---|
| 4 | ! Author: Lionel GUEZ |
---|
| 5 | |
---|
| 6 | implicit none |
---|
| 7 | |
---|
| 8 | private |
---|
| 9 | public regr_lat_time_climoz |
---|
| 10 | |
---|
| 11 | contains |
---|
| 12 | |
---|
| 13 | subroutine regr_lat_time_climoz |
---|
| 14 | |
---|
| 15 | ! "regr_lat_time_climoz" stands for "regrid latitude time |
---|
| 16 | ! climatology ozone". |
---|
| 17 | |
---|
| 18 | ! This procedure reads a climatology of ozone from a NetCDF file, |
---|
| 19 | ! regrids it in latitude and time, and writes the regridded field |
---|
| 20 | ! to a new NetCDF file. |
---|
| 21 | |
---|
| 22 | ! The input field depends on time, pressure level and latitude. |
---|
| 23 | ! We assume that the input field is a step function of latitude |
---|
| 24 | ! and that the input latitude coordinate gives the centers of steps. |
---|
| 25 | ! Regridding in latitude is made by averaging, with a cosine of |
---|
| 26 | ! latitude factor. |
---|
| 27 | ! The target LMDZ latitude grid is the "scalar" grid: "rlatu". |
---|
| 28 | ! The values of "rlatu" are taken to be the centers of intervals. |
---|
| 29 | |
---|
| 30 | ! Regridding in time is by linear interpolation. |
---|
| 31 | ! Monthly values are processed to get daily values, on the basis |
---|
| 32 | ! of a 360-day calendar. |
---|
| 33 | ! The input file may contain either values for 12 months or values |
---|
| 34 | ! for 14 months. |
---|
| 35 | ! If there are 14 months then we assume that we have (in that order): |
---|
| 36 | ! December, January, February, ..., November, December, January |
---|
| 37 | ! We use the first December value to interpolate values between January |
---|
| 38 | ! 1st and mid-January. |
---|
| 39 | ! We use the last January value to interpolate values between |
---|
| 40 | ! mid-December and end of December. |
---|
| 41 | ! If there are only 12 months in the input file than we assume |
---|
| 42 | ! periodicity for interpolation at the beginning and at the end of the |
---|
| 43 | ! year. |
---|
| 44 | |
---|
| 45 | ! We assume that in the input file: |
---|
| 46 | ! -- the latitude is in degrees and strictly monotonic (as all |
---|
| 47 | ! NetCDF coordinate variables should be); |
---|
| 48 | ! -- time increases (even though we do not use values of the input |
---|
| 49 | ! time coordinate); |
---|
| 50 | ! -- pressure is in hPa and in strictly ascending order (even |
---|
| 51 | ! though we do not use pressure values here, we write the unit of |
---|
| 52 | ! pressure in the NetCDF header, and we will use the assumptions later, |
---|
| 53 | ! when we regrid in pressure). |
---|
| 54 | |
---|
| 55 | use regr1_step_av_m, only: regr1_step_av |
---|
| 56 | use regr3_lint_m, only: regr3_lint |
---|
| 57 | use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, handle_err, & |
---|
| 58 | nf95_put_var, nf95_gw_var |
---|
| 59 | use netcdf, only: nf90_nowrite, nf90_get_var |
---|
| 60 | |
---|
| 61 | ! Variables local to the procedure: |
---|
| 62 | |
---|
| 63 | include "dimensions.h" |
---|
| 64 | ! (for "jjm") |
---|
| 65 | include "paramet.h" |
---|
| 66 | ! (for the other included files) |
---|
| 67 | include "comgeom2.h" |
---|
| 68 | ! (for "rlatv") |
---|
| 69 | include "comconst.h" |
---|
| 70 | ! (for "pi") |
---|
| 71 | |
---|
| 72 | integer ncid_in, ncid_out ! NetCDF IDs for input and output files |
---|
| 73 | integer n_plev ! number of pressure levels in the input data |
---|
| 74 | integer n_lat! number of latitudes in the input data |
---|
| 75 | |
---|
| 76 | real, pointer:: latitude(:) |
---|
| 77 | ! (of input data, converted to rad, sorted in strictly ascending order) |
---|
| 78 | |
---|
| 79 | real, allocatable:: lat_in_edg(:) |
---|
| 80 | ! (edges of latitude intervals for input data, in rad, in strictly |
---|
| 81 | ! ascending order) |
---|
| 82 | |
---|
| 83 | real, pointer:: plev(:) ! pressure level of input data |
---|
| 84 | logical desc_lat ! latitude in descending order in the input file |
---|
| 85 | |
---|
| 86 | real, pointer:: o3_in(:, :, :) ! (n_lat, n_plev, 12 or 0:13) |
---|
| 87 | ! (ozone climatology from the input file) |
---|
| 88 | ! ("o3_in(j, l, month)" is at latitude "latitude(j)" and pressure |
---|
| 89 | ! level "plev(l)". "month" is between 1 and 12 or between 0 and 13) |
---|
| 90 | |
---|
| 91 | real, allocatable:: o3_regr_lat(:, :, :) ! (jjm + 1, n_plev, 0:13) |
---|
| 92 | ! (mean of "o3_in" over a latitude interval of LMDZ) |
---|
| 93 | ! (First dimension is latitude interval. |
---|
| 94 | ! The latitude interval for "o3_regr_lat(j,:, :)" contains "rlatu(j)". |
---|
| 95 | ! If "j" is between 2 and "jjm" then the interval is: |
---|
| 96 | ! [rlatv(j), rlatv(j-1)] |
---|
| 97 | ! If "j" is 1 or "jjm + 1" then the interval is: |
---|
| 98 | ! [rlatv(1), pi / 2] |
---|
| 99 | ! or: |
---|
| 100 | ! [- pi / 2, rlatv(jjm)] |
---|
| 101 | ! respectively. |
---|
| 102 | ! "o3_regr_lat(:, l, :)" is for pressure level "plev(l)". |
---|
| 103 | ! Last dimension is month number.) |
---|
| 104 | |
---|
| 105 | real, allocatable:: o3_out(:, :, :) ! (jjm + 1, n_plev, 360) |
---|
| 106 | ! (regridded ozone climatology) |
---|
| 107 | ! ("o3_out(j, l, day)" is at latitude "rlatu(j)", pressure |
---|
| 108 | ! level "plev(l)" and date "January 1st 0h" + "tmidday(day)", in a |
---|
| 109 | ! 360-day calendar.) |
---|
| 110 | |
---|
| 111 | integer j |
---|
| 112 | |
---|
| 113 | integer varid_in, varid_out, varid_plev, varid_time |
---|
| 114 | integer varid |
---|
| 115 | ! (for NetCDF) |
---|
| 116 | |
---|
| 117 | real, parameter:: tmidmonth(0:13) = (/(-15. + 30. * j, j = 0, 13)/) |
---|
| 118 | ! (time to middle of month, in days since January 1st 0h, in a |
---|
| 119 | ! 360-day calendar) |
---|
| 120 | ! (We add values -15 and 375 so that, for example, day 3 of the year is |
---|
| 121 | ! interpolated between the December and the January value.) |
---|
| 122 | |
---|
| 123 | real, parameter:: tmidday(360) = (/(j + 0.5, j = 0, 359)/) |
---|
| 124 | ! (time to middle of day, in days since January 1st 0h, in a |
---|
| 125 | ! 360-day calendar) |
---|
| 126 | |
---|
| 127 | !--------------------------------- |
---|
| 128 | |
---|
| 129 | print *, "Call sequence information: regr_lat_time_climoz" |
---|
| 130 | |
---|
| 131 | call nf95_open("climoz.nc", nf90_nowrite, ncid_in) |
---|
| 132 | |
---|
| 133 | ! Get coordinates from the input file: |
---|
| 134 | |
---|
| 135 | call nf95_inq_varid(ncid_in, "latitude", varid) |
---|
| 136 | call nf95_gw_var(ncid_in, varid, latitude) |
---|
| 137 | ! Convert from degrees to rad, because we will take the sine of latitude: |
---|
| 138 | latitude = latitude / 180. * pi |
---|
| 139 | n_lat = size(latitude) |
---|
| 140 | ! We need to supply the latitudes to "regr1_step_av" in |
---|
| 141 | ! ascending order, so invert order if necessary: |
---|
| 142 | desc_lat = latitude(1) > latitude(n_lat) |
---|
| 143 | if (desc_lat) latitude = latitude(n_lat:1:-1) |
---|
| 144 | |
---|
| 145 | ! Compute edges of latitude intervals: |
---|
| 146 | allocate(lat_in_edg(n_lat + 1)) |
---|
| 147 | lat_in_edg(1) = - pi / 2 |
---|
| 148 | forall (j = 2:n_lat) lat_in_edg(j) = (latitude(j - 1) + latitude(j)) / 2 |
---|
| 149 | lat_in_edg(n_lat + 1) = pi / 2 |
---|
| 150 | deallocate(latitude) ! pointer |
---|
| 151 | |
---|
| 152 | call nf95_inq_varid(ncid_in, "plev", varid) |
---|
| 153 | call nf95_gw_var(ncid_in, varid, plev) |
---|
| 154 | n_plev = size(plev) |
---|
| 155 | ! (We only need the pressure coordinate to copy it to the output file.) |
---|
| 156 | |
---|
| 157 | ! Create the output file and get the variable IDs: |
---|
| 158 | call prepare_out(ncid_in, n_plev, ncid_out, varid_out, varid_plev, & |
---|
| 159 | varid_time) |
---|
| 160 | |
---|
| 161 | ! Write remaining coordinate variables: |
---|
| 162 | call nf95_put_var(ncid_out, varid_plev, plev) |
---|
| 163 | call nf95_put_var(ncid_out, varid_time, tmidday) |
---|
| 164 | |
---|
| 165 | deallocate(plev) ! pointer |
---|
| 166 | |
---|
| 167 | call nf95_inq_varid(ncid_in, "tro3", varid_in) |
---|
| 168 | call nf95_gw_var(ncid_in, varid_in, o3_in) |
---|
| 169 | if (desc_lat) o3_in = o3_in(n_lat:1:-1, :, :) |
---|
| 170 | |
---|
| 171 | call nf95_close(ncid_in) |
---|
| 172 | |
---|
| 173 | allocate(o3_regr_lat(jjm + 1, n_plev, 0:13)) |
---|
| 174 | allocate(o3_out(jjm + 1, n_plev, 360)) |
---|
| 175 | |
---|
| 176 | ! Regrid in latitude: |
---|
| 177 | ! We average with respect to sine of latitude, which is |
---|
| 178 | ! equivalent to weighting by cosine of latitude: |
---|
| 179 | if (size(o3_in, 3) == 12) then |
---|
| 180 | print *, "Found 12 months in ozone climatology, assuming periodicity..." |
---|
| 181 | o3_regr_lat(jjm+1:1:-1, :, 1:12) = regr1_step_av(o3_in, & |
---|
| 182 | xs=sin(lat_in_edg), xt=sin((/- pi / 2, rlatv(jjm:1:-1), pi / 2/))) |
---|
| 183 | ! (invert order of indices in "o3_regr_lat" because "rlatu" is |
---|
| 184 | ! in descending order) |
---|
| 185 | |
---|
| 186 | ! Duplicate January and December values, in preparation of time |
---|
| 187 | ! interpolation: |
---|
| 188 | o3_regr_lat(:, :, 0) = o3_regr_lat(:, :, 12) |
---|
| 189 | o3_regr_lat(:, :, 13) = o3_regr_lat(:, :, 1) |
---|
| 190 | else |
---|
| 191 | print *, "Using 14 months in ozone climatology..." |
---|
| 192 | o3_regr_lat(jjm+1:1:-1, :, :) = regr1_step_av(o3_in, & |
---|
| 193 | xs=sin(lat_in_edg), xt=sin((/- pi / 2, rlatv(jjm:1:-1), pi / 2/))) |
---|
| 194 | ! (invert order of indices in "o3_regr_lat" because "rlatu" is |
---|
| 195 | ! in descending order) |
---|
| 196 | end if |
---|
| 197 | |
---|
| 198 | deallocate(o3_in) ! pointer |
---|
| 199 | |
---|
| 200 | ! Regrid in time by linear interpolation: |
---|
| 201 | o3_out = regr3_lint(o3_regr_lat, tmidmonth, tmidday) |
---|
| 202 | |
---|
| 203 | ! Write to file: |
---|
| 204 | call nf95_put_var(ncid_out, varid_out, o3_out(jjm+1:1:-1, :, :)) |
---|
| 205 | ! (The order of "rlatu" is inverted in the output file) |
---|
| 206 | |
---|
| 207 | call nf95_close(ncid_out) |
---|
| 208 | |
---|
| 209 | end subroutine regr_lat_time_climoz |
---|
| 210 | |
---|
| 211 | !******************************************** |
---|
| 212 | |
---|
| 213 | subroutine prepare_out(ncid_in, n_plev, ncid_out, varid_out, varid_plev, & |
---|
| 214 | varid_time) |
---|
| 215 | |
---|
| 216 | ! This subroutine creates the NetCDF output file, defines |
---|
| 217 | ! dimensions and variables, and writes one of the coordinate variables. |
---|
| 218 | |
---|
| 219 | use netcdf95, only: nf95_create, nf95_def_dim, nf95_def_var, & |
---|
| 220 | nf95_put_att, nf95_enddef, nf95_copy_att, nf95_put_var |
---|
| 221 | use netcdf, only: nf90_clobber, nf90_float, nf90_copy_att, nf90_global |
---|
| 222 | |
---|
| 223 | integer, intent(in):: ncid_in, n_plev |
---|
| 224 | integer, intent(out):: ncid_out, varid_out, varid_plev, varid_time |
---|
| 225 | |
---|
| 226 | ! Variables local to the procedure: |
---|
| 227 | |
---|
| 228 | include "dimensions.h" |
---|
| 229 | ! (for "jjm") |
---|
| 230 | include "paramet.h" |
---|
| 231 | ! (for the other included files) |
---|
| 232 | include "comgeom2.h" |
---|
| 233 | ! (for "rlatu") |
---|
| 234 | include "comconst.h" |
---|
| 235 | ! (for "pi") |
---|
| 236 | |
---|
| 237 | integer ncerr |
---|
| 238 | integer dimid_rlatu, dimid_plev, dimid_time |
---|
| 239 | integer varid_rlatu |
---|
| 240 | |
---|
| 241 | !--------------------------- |
---|
| 242 | |
---|
| 243 | print *, "Call sequence information: prepare_out" |
---|
| 244 | |
---|
| 245 | call nf95_create("climoz_LMDZ.nc", nf90_clobber, ncid_out) |
---|
| 246 | |
---|
| 247 | ! Dimensions: |
---|
| 248 | call nf95_def_dim(ncid_out, "time", 360, dimid_time) |
---|
| 249 | call nf95_def_dim(ncid_out, "plev", n_plev, dimid_plev) |
---|
| 250 | call nf95_def_dim(ncid_out, "rlatu", jjm + 1, dimid_rlatu) |
---|
| 251 | |
---|
| 252 | ! Define coordinate variables: |
---|
| 253 | |
---|
| 254 | call nf95_def_var(ncid_out, "time", nf90_float, dimid_time, varid_time) |
---|
| 255 | call nf95_put_att(ncid_out, varid_time, "units", "days since 2000-1-1") |
---|
| 256 | call nf95_put_att(ncid_out, varid_time, "calendar", "360_day") |
---|
| 257 | call nf95_put_att(ncid_out, varid_time, "standard_name", "time") |
---|
| 258 | |
---|
| 259 | call nf95_def_var(ncid_out, "plev", nf90_float, dimid_plev, varid_plev) |
---|
| 260 | call nf95_put_att(ncid_out, varid_plev, "units", "millibar") |
---|
| 261 | call nf95_put_att(ncid_out, varid_plev, "standard_name", "air_pressure") |
---|
| 262 | call nf95_put_att(ncid_out, varid_plev, "long_name", "air pressure") |
---|
| 263 | |
---|
| 264 | call nf95_def_var(ncid_out, "rlatu", nf90_float, dimid_rlatu, varid_rlatu) |
---|
| 265 | call nf95_put_att(ncid_out, varid_rlatu, "units", "degrees_north") |
---|
| 266 | call nf95_put_att(ncid_out, varid_rlatu, "standard_name", "latitude") |
---|
| 267 | |
---|
| 268 | ! Define the primary variable: |
---|
| 269 | |
---|
| 270 | call nf95_def_var(ncid_out, "tro3", nf90_float, & |
---|
| 271 | (/dimid_rlatu, dimid_plev, dimid_time/), varid_out) |
---|
| 272 | call nf95_put_att(ncid_out, varid_out, "long_name", "ozone mole fraction") |
---|
| 273 | call nf95_put_att(ncid_out, varid_out, "standard_name", & |
---|
| 274 | "mole_fraction_of_ozone_in_air") |
---|
| 275 | |
---|
| 276 | ! Global attributes: |
---|
| 277 | |
---|
| 278 | call nf95_put_att(ncid_out, nf90_global, "comment", "Regridded for LMDZ") |
---|
| 279 | |
---|
| 280 | ! The following commands, copying attributes, may fail. |
---|
| 281 | ! That is OK. |
---|
| 282 | ! It should just mean that the attribute is not defined in the input file. |
---|
| 283 | |
---|
| 284 | call nf95_copy_att(ncid_in, nf90_global, "Conventions", ncid_out, & |
---|
| 285 | nf90_global, ncerr) |
---|
| 286 | call handle_err_copy_att("Conventions") |
---|
| 287 | |
---|
| 288 | call nf95_copy_att(ncid_in, nf90_global, "title", ncid_out, nf90_global, & |
---|
| 289 | ncerr) |
---|
| 290 | call handle_err_copy_att("title") |
---|
| 291 | |
---|
| 292 | call nf95_copy_att(ncid_in, nf90_global, "source", ncid_out, nf90_global, & |
---|
| 293 | ncerr) |
---|
| 294 | call handle_err_copy_att("source") |
---|
| 295 | |
---|
| 296 | call nf95_enddef(ncid_out) |
---|
| 297 | |
---|
| 298 | ! Write one of the coordinate variables: |
---|
| 299 | call nf95_put_var(ncid_out, varid_rlatu, rlatu(jjm+1:1:-1) / pi * 180.) |
---|
| 300 | ! (convert from rad to degrees and sort in ascending order) |
---|
| 301 | |
---|
| 302 | contains |
---|
| 303 | |
---|
| 304 | subroutine handle_err_copy_att(att_name) |
---|
| 305 | |
---|
| 306 | use netcdf, only: nf90_noerr, nf90_strerror |
---|
| 307 | |
---|
| 308 | character(len=*), intent(in):: att_name |
---|
| 309 | |
---|
| 310 | !---------------------------------------- |
---|
| 311 | |
---|
| 312 | if (ncerr /= nf90_noerr) then |
---|
| 313 | print *, "regr_lat_time_climoz_m prepare_out nf90_copy_att " & |
---|
| 314 | // att_name // " -- " // trim(nf90_strerror(ncerr)) |
---|
| 315 | end if |
---|
| 316 | |
---|
| 317 | end subroutine handle_err_copy_att |
---|
| 318 | |
---|
| 319 | end subroutine prepare_out |
---|
| 320 | |
---|
| 321 | end module regr_lat_time_climoz_m |
---|