| [282] | 1 | program extract |
|---|
| 2 | |
|---|
| 3 | ! program to extract (ie: interpolates) pointwise values of an atmospheric |
|---|
| 4 | ! variable from a 'zrecast'ed diagfi file (works if altitude is geometrical |
|---|
| 5 | ! height or a pressure vertical coordinates) |
|---|
| 6 | ! user has to specify: |
|---|
| 7 | ! - name of input file |
|---|
| 8 | ! - date (in sols) offset wrt the input file (e.g. if the input file "begins" |
|---|
| [2571] | 9 | ! at Ls=0, then the offset is 0; if the input file begins at Ls=90, the |
|---|
| [282] | 10 | ! offset date corresponding to the first 3 months is 61+66+66=193 sols, etc.) |
|---|
| [2571] | 11 | ! NB: that means the input file must have its "Time" coordinate in sols |
|---|
| [282] | 12 | ! - the "extraction mode": |
|---|
| 13 | ! 1: extract individual values; user will specify values of |
|---|
| 14 | ! lon lat alt Ls LT (all on a same line) |
|---|
| 15 | ! on as many lines as there are sought values |
|---|
| 16 | ! 2: extract a profile: user will specify on a first line the values of |
|---|
| 17 | ! lon lat Ls LT (all on a same line) |
|---|
| 18 | ! and then only specify values of altitudes (m or Pa depending on the |
|---|
| 19 | ! coordinate in the input file), one per line, at which values are |
|---|
| 20 | ! sought |
|---|
| [293] | 21 | ! - output values are sent to (ASCII) output file 'infile_var_.dat', where |
|---|
| 22 | ! 'infile' is the input file name (without trailing '.nc') and |
|---|
| 23 | ! 'var' is the sought variable, for extraction mode 1 as |
|---|
| 24 | ! lines of "lon lat alt Ls LT value" and for a profile (extraction mode 2) |
|---|
| [282] | 25 | ! as lines of "alt value" |
|---|
| 26 | ! |
|---|
| 27 | ! NB: If there is no data to do an appropriate interpolation to extract |
|---|
| 28 | ! the sought value, then a "missing_value" (taken from the variable's |
|---|
| 29 | ! attribute in the input file, most likely -9.99E33) is returned. |
|---|
| 30 | ! |
|---|
| 31 | ! EM. Sept. 2011 |
|---|
| 32 | |
|---|
| 33 | use netcdf |
|---|
| 34 | implicit none |
|---|
| 35 | |
|---|
| 36 | ! Input file: |
|---|
| 37 | character(len=256) :: infile |
|---|
| [293] | 38 | character(len=256) :: outfile |
|---|
| [282] | 39 | |
|---|
| [312] | 40 | character (len=256) :: text ! to store some text |
|---|
| [282] | 41 | character (len=64) :: varname ! to store the name of the variable to retreive |
|---|
| 42 | |
|---|
| 43 | ! NetCDF stuff |
|---|
| 44 | integer :: status ! NetCDF return code |
|---|
| 45 | integer :: inid ! NetCDF file IDs |
|---|
| 46 | integer :: varid ! to store the ID of a variable |
|---|
| 47 | integer :: lat_dimid,lon_dimid,alt_dimid,time_dimid |
|---|
| 48 | integer :: datashape(4) |
|---|
| 49 | |
|---|
| 50 | real,dimension(:),allocatable :: longitude ! longitude |
|---|
| 51 | integer lonlen ! # of grid points along longitude |
|---|
| 52 | real,dimension(:),allocatable :: latitude ! latitude |
|---|
| 53 | integer latlen ! # of grid points along latitude |
|---|
| 54 | real,dimension(:),allocatable :: altitude ! can be geometric heights or pressure |
|---|
| 55 | integer altlen ! # of grid point along altitude |
|---|
| 56 | real,dimension(:),allocatable :: time ! time |
|---|
| 57 | integer timelen ! # of points along time |
|---|
| 58 | character :: alttype ! altitude coord. type:'z' (altitude, m) 'p' (pressure, Pa) |
|---|
| 59 | real,dimension(:,:,:,:),allocatable :: field |
|---|
| 60 | real :: missing_value ! value to denote non-existant data |
|---|
| 61 | real :: starttimeoffset ! offset (in sols) wrt Ls=0 of sol 0 in file |
|---|
| 62 | integer :: extract_mode ! 1: point-by-point extraction 2: extract a profile |
|---|
| 63 | |
|---|
| 64 | ! point at which data is sought: |
|---|
| 65 | real :: lon,lat,alt,Ls,LT,value |
|---|
| 66 | real :: sol ! sol GCM date corresponding to sought Ls and LT |
|---|
| 67 | |
|---|
| 68 | integer :: nb |
|---|
| 69 | |
|---|
| 70 | !=============================================================================== |
|---|
| 71 | ! 1.1 Input file |
|---|
| 72 | !=============================================================================== |
|---|
| 73 | |
|---|
| 74 | write(*,*) "" |
|---|
| 75 | write(*,*) " Program valid for diagfi.nc or concatnc.nc files" |
|---|
| 76 | write(*,*) " processed by zrecast " |
|---|
| 77 | write(*,*) " Enter input file name:" |
|---|
| 78 | |
|---|
| 79 | read(*,'(a)') infile |
|---|
| 80 | write(*,*) "" |
|---|
| 81 | |
|---|
| 82 | ! open input file |
|---|
| 83 | status=nf90_open(infile,NF90_NOWRITE,inid) |
|---|
| 84 | if (status.ne.nf90_noerr) then |
|---|
| 85 | write(*,*)"Failed to open datafile ",trim(infile) |
|---|
| 86 | write(*,*)trim(nf90_strerror(status)) |
|---|
| 87 | stop |
|---|
| 88 | endif |
|---|
| 89 | |
|---|
| 90 | write(*,*) " Beginning date of the file?" |
|---|
| 91 | write(*,*) " (i.e. number of sols since Ls=0 that the Time=0.0 in the input" |
|---|
| 92 | write(*,*) " file corresponds to)" |
|---|
| 93 | read(*,*) starttimeoffset |
|---|
| 94 | |
|---|
| 95 | write(*,*) " Extraction mode?" |
|---|
| 96 | write(*,*) " ( 1: pointwise extraction , 2: profile extraction)" |
|---|
| 97 | read(*,*,iostat=status) extract_mode |
|---|
| 98 | if ((status.ne.0).or.(extract_mode.lt.1) & |
|---|
| 99 | .or.(extract_mode.gt.2)) then |
|---|
| 100 | write(*,*) "Error: invalid extraction mode:",extract_mode |
|---|
| 101 | stop |
|---|
| 102 | endif |
|---|
| 103 | |
|---|
| 104 | !=============================================================================== |
|---|
| 105 | ! 1.2 Input variable to extract |
|---|
| 106 | !=============================================================================== |
|---|
| 107 | |
|---|
| 108 | write(*,*) "Enter variable to extract:" |
|---|
| 109 | read(*,*) varname |
|---|
| 110 | ! check that input file contains that variable |
|---|
| 111 | status=nf90_inq_varid(inid,trim(varname),varid) |
|---|
| 112 | if (status.ne.nf90_noerr) then |
|---|
| 113 | write(*,*) "Failed to find variable ",trim(varname)," in ",trim(infile) |
|---|
| 114 | write(*,*)trim(nf90_strerror(status)) |
|---|
| 115 | write(*,*) " Might as well stop here." |
|---|
| 116 | stop |
|---|
| 117 | endif |
|---|
| 118 | |
|---|
| 119 | !=============================================================================== |
|---|
| 120 | ! 1.3 Get grids for dimensions lon,lat,alt,time |
|---|
| 121 | !=============================================================================== |
|---|
| 122 | |
|---|
| 123 | ! latitude |
|---|
| 124 | status=nf90_inq_dimid(inid,"latitude",lat_dimid) |
|---|
| 125 | if (status.ne.nf90_noerr) then |
|---|
| 126 | write(*,*)"Failed to find latitude dimension" |
|---|
| 127 | write(*,*)trim(nf90_strerror(status)) |
|---|
| 128 | stop |
|---|
| 129 | endif |
|---|
| 130 | status=nf90_inquire_dimension(inid,lat_dimid,len=latlen) |
|---|
| 131 | if (status.ne.nf90_noerr) then |
|---|
| 132 | write(*,*)"Failed to find latitude length" |
|---|
| 133 | write(*,*)trim(nf90_strerror(status)) |
|---|
| 134 | endif |
|---|
| 135 | allocate(latitude(latlen)) |
|---|
| 136 | status=nf90_inq_varid(inid,"latitude",varid) |
|---|
| 137 | if (status.ne.nf90_noerr) then |
|---|
| 138 | write(*,*) "Failed to find latitude ID" |
|---|
| 139 | write(*,*)trim(nf90_strerror(status)) |
|---|
| 140 | stop |
|---|
| 141 | endif |
|---|
| 142 | ! read latitude |
|---|
| 143 | status=nf90_get_var(inid,varid,latitude) |
|---|
| 144 | if (status.ne.nf90_noerr) then |
|---|
| 145 | write(*,*) "Failed to load latitude" |
|---|
| 146 | write(*,*)trim(nf90_strerror(status)) |
|---|
| 147 | stop |
|---|
| 148 | endif |
|---|
| 149 | |
|---|
| 150 | !longitude |
|---|
| 151 | status=nf90_inq_dimid(inid,"longitude",lon_dimid) |
|---|
| 152 | if (status.ne.nf90_noerr) then |
|---|
| 153 | write(*,*)"Failed to find longitude dimension" |
|---|
| 154 | write(*,*)trim(nf90_strerror(status)) |
|---|
| 155 | stop |
|---|
| 156 | endif |
|---|
| 157 | status=nf90_inquire_dimension(inid,lon_dimid,len=lonlen) |
|---|
| 158 | if (status.ne.nf90_noerr) then |
|---|
| 159 | write(*,*)"Failed to find longitude length" |
|---|
| 160 | write(*,*)trim(nf90_strerror(status)) |
|---|
| 161 | endif |
|---|
| 162 | allocate(longitude(lonlen)) |
|---|
| 163 | status=nf90_inq_varid(inid,"longitude",varid) |
|---|
| 164 | if (status.ne.nf90_noerr) then |
|---|
| 165 | write(*,*) "Failed to find longitude ID" |
|---|
| 166 | write(*,*)trim(nf90_strerror(status)) |
|---|
| 167 | stop |
|---|
| 168 | endif |
|---|
| 169 | ! read longitude |
|---|
| 170 | status=nf90_get_var(inid,varid,longitude) |
|---|
| 171 | if (status.ne.nf90_noerr) then |
|---|
| 172 | write(*,*) "Failed to load longitude" |
|---|
| 173 | write(*,*)trim(nf90_strerror(status)) |
|---|
| 174 | stop |
|---|
| 175 | endif |
|---|
| 176 | |
|---|
| 177 | !time |
|---|
| 178 | status=nf90_inq_dimid(inid,"Time",time_dimid) |
|---|
| 179 | if (status.ne.nf90_noerr) then |
|---|
| 180 | write(*,*)"Failed to find Time dimension" |
|---|
| 181 | write(*,*)trim(nf90_strerror(status)) |
|---|
| 182 | stop |
|---|
| 183 | endif |
|---|
| 184 | status=nf90_inquire_dimension(inid,time_dimid,len=timelen) |
|---|
| 185 | if (status.ne.nf90_noerr) then |
|---|
| 186 | write(*,*)"Failed to find Time length" |
|---|
| 187 | write(*,*)trim(nf90_strerror(status)) |
|---|
| 188 | endif |
|---|
| 189 | allocate(time(timelen)) |
|---|
| 190 | status=nf90_inq_varid(inid,"Time",varid) |
|---|
| 191 | if (status.ne.nf90_noerr) then |
|---|
| 192 | write(*,*) "Failed to find Time ID" |
|---|
| 193 | write(*,*)trim(nf90_strerror(status)) |
|---|
| 194 | stop |
|---|
| 195 | endif |
|---|
| 196 | ! read Time |
|---|
| 197 | status=nf90_get_var(inid,varid,time) |
|---|
| 198 | if (status.ne.nf90_noerr) then |
|---|
| 199 | write(*,*) "Failed to load Time" |
|---|
| 200 | write(*,*)trim(nf90_strerror(status)) |
|---|
| 201 | stop |
|---|
| 202 | endif |
|---|
| 203 | ! add the offset to time(:) |
|---|
| 204 | time(:)=time(:)+starttimeoffset |
|---|
| 205 | |
|---|
| 206 | !altitude |
|---|
| 207 | status=nf90_inq_dimid(inid,"altitude",alt_dimid) |
|---|
| 208 | if (status.ne.nf90_noerr) then |
|---|
| 209 | write(*,*)"Failed to find altitude dimension" |
|---|
| 210 | write(*,*)trim(nf90_strerror(status)) |
|---|
| 211 | stop |
|---|
| 212 | endif |
|---|
| 213 | status=nf90_inquire_dimension(inid,alt_dimid,len=altlen) |
|---|
| 214 | if (status.ne.nf90_noerr) then |
|---|
| 215 | write(*,*)"Failed to find altitude length" |
|---|
| 216 | write(*,*)trim(nf90_strerror(status)) |
|---|
| 217 | endif |
|---|
| 218 | allocate(altitude(altlen)) |
|---|
| 219 | status=nf90_inq_varid(inid,"altitude",varid) |
|---|
| 220 | if (status.ne.nf90_noerr) then |
|---|
| 221 | write(*,*) "Failed to find altitude ID" |
|---|
| 222 | write(*,*)trim(nf90_strerror(status)) |
|---|
| 223 | stop |
|---|
| 224 | endif |
|---|
| 225 | ! read altitude |
|---|
| 226 | status=nf90_get_var(inid,varid,altitude) |
|---|
| 227 | if (status.ne.nf90_noerr) then |
|---|
| 228 | write(*,*) "Failed to load altitude" |
|---|
| 229 | write(*,*)trim(nf90_strerror(status)) |
|---|
| 230 | stop |
|---|
| 231 | endif |
|---|
| 232 | ! check altitude attribute "units" to find out altitude type |
|---|
| 233 | status=nf90_get_att(inid,varid,"units",text) |
|---|
| 234 | if (status.ne.nf90_noerr) then |
|---|
| 235 | write(*,*) "Failed to load altitude units attribute" |
|---|
| 236 | write(*,*)trim(nf90_strerror(status)) |
|---|
| 237 | stop |
|---|
| 238 | else |
|---|
| 239 | if (trim(text).eq."Pa") then |
|---|
| 240 | alttype="p" ! pressure coordinate |
|---|
| 241 | else if (trim(text).eq."m") then |
|---|
| 242 | alttype="z" ! altitude coordinate |
|---|
| 243 | else |
|---|
| 244 | write(*,*)" I do not understand this unit ",trim(text)," for altitude!" |
|---|
| 245 | stop |
|---|
| 246 | endif |
|---|
| 247 | endif |
|---|
| 248 | |
|---|
| 249 | !=============================================================================== |
|---|
| 250 | ! 1.3 Get input dataset |
|---|
| 251 | !=============================================================================== |
|---|
| 252 | status=nf90_inq_varid(inid,trim(varname),varid) |
|---|
| 253 | if (status.ne.nf90_noerr) then |
|---|
| 254 | write(*,*) "Failed to find variable ",trim(varname)," in ",trim(infile) |
|---|
| 255 | write(*,*)trim(nf90_strerror(status)) |
|---|
| 256 | write(*,*) " Might as well stop here." |
|---|
| 257 | stop |
|---|
| 258 | endif |
|---|
| 259 | |
|---|
| 260 | ! sanity checks on the variable |
|---|
| 261 | status=nf90_inquire_variable(inid,varid,ndims=nb,dimids=datashape) |
|---|
| 262 | if (status.ne.nf90_noerr) then |
|---|
| 263 | write(*,*) "Failed to obtain information on variable ",trim(varname) |
|---|
| 264 | write(*,*)trim(nf90_strerror(status)) |
|---|
| 265 | write(*,*) " Might as well stop here." |
|---|
| 266 | stop |
|---|
| 267 | else |
|---|
| 268 | ! check that it is a 4D variable |
|---|
| 269 | if (nb.ne.4) then |
|---|
| 270 | write(*,*) "Error, expected a 4D (lon-lat-alt-time) variable!" |
|---|
| 271 | stop |
|---|
| 272 | endif |
|---|
| 273 | ! check that its dimensions are indeed lon,lat,alt,time |
|---|
| 274 | if (datashape(1).ne.lon_dimid) then |
|---|
| 275 | write(*,*) "Error, expected first dimension to be longitude!" |
|---|
| 276 | stop |
|---|
| 277 | endif |
|---|
| 278 | if (datashape(2).ne.lat_dimid) then |
|---|
| 279 | write(*,*) "Error, expected second dimension to be latitude!" |
|---|
| 280 | stop |
|---|
| 281 | endif |
|---|
| 282 | if (datashape(3).ne.alt_dimid) then |
|---|
| 283 | write(*,*) "Error, expected third dimension to be altitude!" |
|---|
| 284 | stop |
|---|
| 285 | endif |
|---|
| 286 | if (datashape(4).ne.time_dimid) then |
|---|
| 287 | write(*,*) "Error, expected fourth dimension to be time!" |
|---|
| 288 | stop |
|---|
| 289 | endif |
|---|
| 290 | endif |
|---|
| 291 | |
|---|
| 292 | allocate(field(lonlen,latlen,altlen,timelen)) |
|---|
| 293 | |
|---|
| 294 | ! load dataset |
|---|
| 295 | status=nf90_get_var(inid,varid,field) |
|---|
| 296 | if (status.ne.nf90_noerr) then |
|---|
| 297 | write(*,*) "Failed to load ",trim(varname) |
|---|
| 298 | write(*,*)trim(nf90_strerror(status)) |
|---|
| 299 | stop |
|---|
| 300 | else |
|---|
| 301 | write(*,*) "Loaded ",trim(varname) |
|---|
| 302 | endif |
|---|
| 303 | ! get dataset's missing_value attribute |
|---|
| 304 | status=nf90_get_att(inid,varid,"missing_value",missing_value) |
|---|
| 305 | if (status.ne.nf90_noerr) then |
|---|
| 306 | write(*,*) "Failed to load missing_value attribute" |
|---|
| 307 | write(*,*)trim(nf90_strerror(status)) |
|---|
| 308 | stop |
|---|
| 309 | else |
|---|
| [2140] | 310 | write(*,'(" with missing_value attribute : ",1pe12.5)') missing_value |
|---|
| [282] | 311 | endif |
|---|
| 312 | |
|---|
| 313 | !=============================================================================== |
|---|
| 314 | ! 2. Process and interpolate |
|---|
| 315 | !=============================================================================== |
|---|
| [293] | 316 | |
|---|
| 317 | !=============================================================================== |
|---|
| 318 | ! 2.1 Create output file |
|---|
| 319 | !=============================================================================== |
|---|
| 320 | |
|---|
| 321 | outfile=trim(infile(1:index(infile,".nc",back=.true.)-1))//"_"//& |
|---|
| 322 | trim(varname)//".dat" |
|---|
| 323 | open(42,file=outfile,form="formatted") |
|---|
| 324 | write(*,*) "Output file is: ",trim(outfile) |
|---|
| 325 | |
|---|
| 326 | !=============================================================================== |
|---|
| 327 | ! 2.2 Extract values |
|---|
| 328 | !=============================================================================== |
|---|
| 329 | |
|---|
| [282] | 330 | if (extract_mode==1) then ! pointwise extraction |
|---|
| 331 | write(*,*) " Enter values of: lon lat alt Ls LT (all on the same line," |
|---|
| 332 | write(*,*) " using as many lines as sought data values)" |
|---|
| 333 | write(*,*) " (enter anything else, e.g. blank line, nonesense,... to quit)" |
|---|
| 334 | else ! extract_mode==2 , profile extraction |
|---|
| 335 | write(*,*) " Enter values of: lon lat Ls LT (all on the same line)" |
|---|
| 336 | read(*,*) lon,lat,Ls,LT |
|---|
| 337 | write(*,*) " Enter values of altitude (one per line)" |
|---|
| 338 | write(*,*) " (enter anything else, e.g. blank line, nonesense,... to quit)" |
|---|
| 339 | endif |
|---|
| 340 | do |
|---|
| 341 | ! 2.1 read coordinates and do some sanity checks |
|---|
| [312] | 342 | text="" !initialize text to empty string |
|---|
| 343 | read(*,'(a)') text ! store input as text (to spot empty input line) |
|---|
| 344 | if (len_trim(adjustl(text)).eq.0) exit |
|---|
| 345 | |
|---|
| [282] | 346 | if (extract_mode==1) then |
|---|
| [312] | 347 | read(text,*,iostat=status) lon,lat,alt,Ls,LT |
|---|
| [282] | 348 | ! Ls in degrees, LT (local true solar time) in hours, i.e. in [0:24] |
|---|
| 349 | else ! extract_mode==2 , read alt only |
|---|
| [312] | 350 | read(text,*,iostat=status) alt |
|---|
| [282] | 351 | endif |
|---|
| 352 | if (status.ne.0) exit |
|---|
| 353 | |
|---|
| 354 | if ((lon.lt.-360.).or.(lon.gt.360.)) then |
|---|
| 355 | write(*,*) "Unexpected value for lon: ",lon |
|---|
| 356 | stop |
|---|
| 357 | endif |
|---|
| 358 | ! we want lon in [-180:180] |
|---|
| 359 | if (lon.lt.-180.) lon=lon+360. |
|---|
| [954] | 360 | if (lon.gt.180.) lon=lon-360. |
|---|
| [282] | 361 | |
|---|
| 362 | if ((lat.lt.-90.).or.(lat.gt.90.)) then |
|---|
| 363 | write(*,*) "Unexpected value for lat: ",lat |
|---|
| 364 | stop |
|---|
| 365 | endif |
|---|
| 366 | |
|---|
| 367 | if ((Ls.lt.0.).or.(Ls.gt.360.)) then |
|---|
| 368 | write(*,*) "Unexpected value for Ls: ",Ls |
|---|
| 369 | stop |
|---|
| 370 | endif |
|---|
| 371 | |
|---|
| 372 | if ((LT.lt.0.).or.(LT.gt.24.)) then |
|---|
| 373 | write(*,*) "Unexpected value for LT: ",LT |
|---|
| 374 | stop |
|---|
| 375 | endif |
|---|
| 376 | |
|---|
| 377 | ! 2.2 compute GCM sol date corresponding to sought Ls and LT |
|---|
| 378 | call ls2sol(Ls,sol) |
|---|
| 379 | !shift 'sol' decimal part to ensure a compatible LT with the desired one |
|---|
| 380 | sol=floor(sol)+(LT-lon/15.)/24. |
|---|
| 381 | ! handle bordeline cases: |
|---|
| 382 | if (sol.gt.669) sol=sol-669 |
|---|
| 383 | if (sol.lt.0) sol=sol+669 |
|---|
| 384 | |
|---|
| 385 | ! write(*,*) " Ls=",Ls," LT=",LT," => sol=",sol |
|---|
| 386 | |
|---|
| 387 | ! 2.3 do the interpolation |
|---|
| 388 | call extraction(lon,lat,alt,sol,& |
|---|
| 389 | lonlen,latlen,altlen,timelen,& |
|---|
| 390 | longitude,latitude,altitude,time,& |
|---|
| 391 | field,missing_value,alttype,varname,value) |
|---|
| 392 | ! 2.4 Write value to output |
|---|
| 393 | if (extract_mode==1) then ! pointwise extraction |
|---|
| [293] | 394 | write(42,'(6(1x,1pe12.5))')lon,lat,alt,Ls,LT,value |
|---|
| [282] | 395 | else ! profile |
|---|
| [293] | 396 | write(42,'(6(1x,1pe12.5))')alt,value |
|---|
| [282] | 397 | endif |
|---|
| 398 | enddo ! of do while |
|---|
| 399 | |
|---|
| [293] | 400 | ! close output file |
|---|
| 401 | close(42) |
|---|
| [282] | 402 | |
|---|
| 403 | end program extract |
|---|
| 404 | |
|---|
| 405 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 406 | subroutine extraction(lon,lat,alt,sol,& |
|---|
| 407 | lonlen,latlen,altlen,timelen,& |
|---|
| 408 | longitude,latitude,altitude,time,& |
|---|
| 409 | field,missing_value,alttype,varname,value) |
|---|
| 410 | |
|---|
| 411 | implicit none |
|---|
| 412 | ! Arguments: |
|---|
| 413 | real,intent(in) :: lon ! sought longitude (deg, in [-180:180]) |
|---|
| 414 | real,intent(in) :: lat ! sought latitude (deg, in [-90:90]) |
|---|
| 415 | real,intent(in) :: alt ! sought altitude (m or Pa) |
|---|
| 416 | real,intent(in) :: sol ! sought date (sols) |
|---|
| 417 | integer,intent(in) :: lonlen |
|---|
| 418 | integer,intent(in) :: latlen |
|---|
| 419 | integer,intent(in) :: altlen |
|---|
| 420 | integer,intent(in) :: timelen |
|---|
| 421 | real,intent(in) :: longitude(lonlen) |
|---|
| 422 | real,intent(in) :: latitude(latlen) |
|---|
| 423 | real,intent(in) :: altitude(altlen) |
|---|
| 424 | real,intent(in) :: time(timelen) |
|---|
| 425 | real,intent(in) :: field(lonlen,latlen,altlen,timelen) |
|---|
| 426 | real,intent(in) :: missing_value ! default value in GCM file for "no data" |
|---|
| 427 | character,intent(in) :: alttype ! altitude coord. type:'z' (altitude, m) |
|---|
| 428 | ! 'p' (pressure, Pa) |
|---|
| 429 | character(len=*),intent(in) :: varname ! variable name (in GCM file) |
|---|
| 430 | real,intent(out) :: value |
|---|
| 431 | |
|---|
| 432 | ! Local variables: |
|---|
| 433 | real,save :: prev_lon=-99 ! previous value of 'lon' routine was called with |
|---|
| 434 | real,save :: prev_lat=-99 ! previous value of 'lat' routine was called with |
|---|
| 435 | real,save :: prev_alt=-9.e20 ! ! previous value of 'alt' |
|---|
| 436 | real,save :: prev_sol=-99 ! previous value of 'sol' routine was called with |
|---|
| 437 | |
|---|
| 438 | ! encompasing indexes: |
|---|
| 439 | integer,save :: ilon_inf=-1,ilon_sup=-1 ! along longitude |
|---|
| 440 | integer,save :: ilat_inf=-1,ilat_sup=-1 ! along latitude |
|---|
| 441 | integer,save :: ialt_inf=-1,ialt_sup=-1 ! along altitude |
|---|
| 442 | integer,save :: itim_inf=-1,itim_sup=-1 ! along time |
|---|
| 443 | |
|---|
| 444 | ! intermediate interpolated values |
|---|
| 445 | real :: t_interp(2,2,2) ! after time interpolation |
|---|
| 446 | real :: zt_interp(2,2) ! after altitude interpolation |
|---|
| 447 | real :: yzt_interp(2) ! after latitude interpolation |
|---|
| 448 | real :: coeff ! interpolation coefficient |
|---|
| 449 | |
|---|
| 450 | integer :: i |
|---|
| 451 | |
|---|
| 452 | ! By default, set value to missing_value |
|---|
| 453 | value=missing_value |
|---|
| 454 | |
|---|
| 455 | ! 1. Find encompassing indexes |
|---|
| 456 | if (lon.ne.prev_lon) then |
|---|
| 457 | do i=1,lonlen-1 |
|---|
| 458 | if (longitude(i).le.lon) then |
|---|
| 459 | ilon_inf=i |
|---|
| 460 | else |
|---|
| 461 | exit |
|---|
| 462 | endif |
|---|
| 463 | enddo |
|---|
| 464 | ilon_sup=ilon_inf+1 |
|---|
| 465 | endif ! of if (lon.ne.prev_lon) |
|---|
| 466 | !write(*,*) 'ilon_inf=',ilon_inf," longitude(ilon_inf)=",longitude(ilon_inf) |
|---|
| 467 | |
|---|
| 468 | if (lat.ne.prev_lat) then |
|---|
| 469 | ! recall that latitudes start from north pole |
|---|
| 470 | do i=1,latlen-1 |
|---|
| 471 | if (latitude(i).ge.lat) then |
|---|
| 472 | ilat_inf=i |
|---|
| 473 | else |
|---|
| 474 | exit |
|---|
| 475 | endif |
|---|
| 476 | enddo |
|---|
| 477 | ilat_sup=ilat_inf+1 |
|---|
| 478 | endif ! of if (lat.ne.prev_lat) |
|---|
| 479 | !write(*,*) 'ilat_inf=',ilat_inf," latitude(ilat_inf)=",latitude(ilat_inf) |
|---|
| 480 | |
|---|
| 481 | if (alt.ne.prev_alt) then |
|---|
| 482 | if (alttype.eq.'p') then ! pressures are ordered from max to min |
|---|
| 483 | !handle special case for alt not in the altitude(1:altlen) interval |
|---|
| 484 | if ((alt.gt.altitude(1)).or.(alt.lt.altitude(altlen))) then |
|---|
| 485 | ialt_inf=-1 |
|---|
| 486 | ialt_sup=-1 |
|---|
| 487 | ! return to main program (with value=missing_value) |
|---|
| 488 | return |
|---|
| 489 | else ! general case |
|---|
| 490 | do i=1,altlen-1 |
|---|
| 491 | if (altitude(i).ge.alt) then |
|---|
| 492 | ialt_inf=i |
|---|
| 493 | else |
|---|
| 494 | exit |
|---|
| 495 | endif |
|---|
| 496 | enddo |
|---|
| 497 | ialt_sup=ialt_inf+1 |
|---|
| 498 | endif ! of if ((alt.gt.altitude(1)).or.(alt.lt.altitude(altlen))) |
|---|
| 499 | else ! altitudes (m) are ordered from min to max |
|---|
| 500 | !handle special case for alt not in the altitude(1:altlen) interval |
|---|
| 501 | if ((alt.lt.altitude(1)).or.(alt.gt.altitude(altlen))) then |
|---|
| 502 | ialt_inf=-1 |
|---|
| 503 | ialt_sup=-1 |
|---|
| 504 | ! return to main program (with value=missing_value) |
|---|
| 505 | return |
|---|
| 506 | else ! general case |
|---|
| 507 | do i=1,altlen-1 |
|---|
| 508 | if (altitude(i).le.alt) then |
|---|
| 509 | ialt_inf=i |
|---|
| 510 | else |
|---|
| 511 | exit |
|---|
| 512 | endif |
|---|
| 513 | enddo |
|---|
| 514 | ialt_sup=ialt_inf+1 |
|---|
| 515 | endif ! of if ((alt.lt.altitude(1)).or.(alt.gt.altitude(altlen))) |
|---|
| 516 | endif ! of if (alttype.eq.'p') |
|---|
| 517 | endif ! of if (alt.ne.prev_alt) |
|---|
| 518 | !write(*,*) 'ialt_inf=',ialt_inf," altitude(ialt_inf)=",altitude(ialt_inf) |
|---|
| 519 | |
|---|
| 520 | if (sol.ne.prev_sol) then |
|---|
| 521 | !handle special case for sol not in the time(1):time(timenlen) interval |
|---|
| 522 | if ((sol.lt.time(1)).or.(sol.gt.time(timelen))) then |
|---|
| 523 | itim_inf=-1 |
|---|
| 524 | itim_sup=-1 |
|---|
| 525 | ! return to main program (with value=missing_value) |
|---|
| 526 | return |
|---|
| 527 | else ! general case |
|---|
| 528 | do i=1,timelen-1 |
|---|
| 529 | if (time(i).le.sol) then |
|---|
| 530 | itim_inf=i |
|---|
| 531 | else |
|---|
| 532 | exit |
|---|
| 533 | endif |
|---|
| 534 | enddo |
|---|
| 535 | itim_sup=itim_inf+1 |
|---|
| 536 | endif ! of if ((sol.lt.time(1)).or.(sol.gt.time(timenlen))) |
|---|
| 537 | endif ! of if (sol.ne.prev_sol) |
|---|
| 538 | !write(*,*) 'itim_inf=',itim_inf," time(itim_inf)=",time(itim_inf) |
|---|
| 539 | !write(*,*) 'itim_sup=',itim_sup," time(itim_sup)=",time(itim_sup) |
|---|
| 540 | |
|---|
| 541 | ! 2. Interpolate |
|---|
| 542 | ! check that there are no "missing_value" in the field() elements we need |
|---|
| 543 | ! otherwise return to main program (with value=missing_value) |
|---|
| 544 | if (field(ilon_inf,ilat_inf,ialt_inf,itim_inf).eq.missing_value) return |
|---|
| 545 | if (field(ilon_inf,ilat_inf,ialt_inf,itim_sup).eq.missing_value) return |
|---|
| 546 | if (field(ilon_inf,ilat_inf,ialt_sup,itim_inf).eq.missing_value) return |
|---|
| 547 | if (field(ilon_inf,ilat_inf,ialt_sup,itim_sup).eq.missing_value) return |
|---|
| 548 | if (field(ilon_inf,ilat_sup,ialt_inf,itim_inf).eq.missing_value) return |
|---|
| 549 | if (field(ilon_inf,ilat_sup,ialt_inf,itim_sup).eq.missing_value) return |
|---|
| 550 | if (field(ilon_inf,ilat_sup,ialt_sup,itim_inf).eq.missing_value) return |
|---|
| 551 | if (field(ilon_inf,ilat_sup,ialt_sup,itim_sup).eq.missing_value) return |
|---|
| 552 | if (field(ilon_sup,ilat_inf,ialt_inf,itim_inf).eq.missing_value) return |
|---|
| 553 | if (field(ilon_sup,ilat_inf,ialt_inf,itim_sup).eq.missing_value) return |
|---|
| 554 | if (field(ilon_sup,ilat_inf,ialt_sup,itim_inf).eq.missing_value) return |
|---|
| 555 | if (field(ilon_sup,ilat_inf,ialt_sup,itim_sup).eq.missing_value) return |
|---|
| 556 | if (field(ilon_sup,ilat_sup,ialt_inf,itim_inf).eq.missing_value) return |
|---|
| 557 | if (field(ilon_sup,ilat_sup,ialt_inf,itim_sup).eq.missing_value) return |
|---|
| 558 | if (field(ilon_sup,ilat_sup,ialt_sup,itim_inf).eq.missing_value) return |
|---|
| 559 | if (field(ilon_sup,ilat_sup,ialt_sup,itim_sup).eq.missing_value) return |
|---|
| 560 | |
|---|
| 561 | ! 2.1 Linear interpolation in time |
|---|
| 562 | coeff=(sol-time(itim_inf))/(time(itim_sup)-time(itim_inf)) |
|---|
| 563 | t_interp(1,1,1)=field(ilon_inf,ilat_inf,ialt_inf,itim_inf)+ & |
|---|
| 564 | coeff*(field(ilon_inf,ilat_inf,ialt_inf,itim_sup)- & |
|---|
| 565 | field(ilon_inf,ilat_inf,ialt_inf,itim_inf)) |
|---|
| 566 | t_interp(1,1,2)=field(ilon_inf,ilat_inf,ialt_sup,itim_inf)+ & |
|---|
| 567 | coeff*(field(ilon_inf,ilat_inf,ialt_sup,itim_sup)- & |
|---|
| 568 | field(ilon_inf,ilat_inf,ialt_sup,itim_inf)) |
|---|
| 569 | t_interp(1,2,1)=field(ilon_inf,ilat_sup,ialt_inf,itim_inf)+ & |
|---|
| 570 | coeff*(field(ilon_inf,ilat_sup,ialt_inf,itim_sup)- & |
|---|
| 571 | field(ilon_inf,ilat_sup,ialt_inf,itim_inf)) |
|---|
| 572 | t_interp(1,2,2)=field(ilon_inf,ilat_sup,ialt_sup,itim_inf)+ & |
|---|
| 573 | coeff*(field(ilon_inf,ilat_sup,ialt_sup,itim_sup)- & |
|---|
| 574 | field(ilon_inf,ilat_sup,ialt_sup,itim_inf)) |
|---|
| 575 | t_interp(2,1,1)=field(ilon_sup,ilat_inf,ialt_inf,itim_inf)+ & |
|---|
| 576 | coeff*(field(ilon_sup,ilat_inf,ialt_inf,itim_sup)- & |
|---|
| 577 | field(ilon_sup,ilat_inf,ialt_inf,itim_inf)) |
|---|
| 578 | t_interp(2,1,2)=field(ilon_sup,ilat_inf,ialt_sup,itim_inf)+ & |
|---|
| 579 | coeff*(field(ilon_sup,ilat_inf,ialt_sup,itim_sup)- & |
|---|
| 580 | field(ilon_sup,ilat_inf,ialt_sup,itim_inf)) |
|---|
| 581 | t_interp(2,2,1)=field(ilon_sup,ilat_sup,ialt_inf,itim_inf)+ & |
|---|
| 582 | coeff*(field(ilon_sup,ilat_sup,ialt_inf,itim_sup)- & |
|---|
| 583 | field(ilon_sup,ilat_sup,ialt_inf,itim_inf)) |
|---|
| 584 | t_interp(2,2,2)=field(ilon_sup,ilat_sup,ialt_sup,itim_inf)+ & |
|---|
| 585 | coeff*(field(ilon_sup,ilat_sup,ialt_sup,itim_sup)- & |
|---|
| 586 | field(ilon_sup,ilat_sup,ialt_sup,itim_inf)) |
|---|
| 587 | |
|---|
| 588 | ! 2.2 Vertical interpolation |
|---|
| 589 | if (((varname=='rho').or.(varname=='pressure')).and.(alttype=='z')) then |
|---|
| 590 | ! do the interpolation on the log of the quantity |
|---|
| 591 | coeff=(alt-altitude(ialt_inf))/(altitude(ialt_sup)-altitude(ialt_inf)) |
|---|
| 592 | zt_interp(1,1)=log(t_interp(1,1,1))+coeff* & |
|---|
| 593 | (log(t_interp(1,1,2))-log(t_interp(1,1,1))) |
|---|
| 594 | zt_interp(1,2)=log(t_interp(1,2,1))+coeff* & |
|---|
| 595 | (log(t_interp(1,2,2))-log(t_interp(1,2,1))) |
|---|
| 596 | zt_interp(2,1)=log(t_interp(2,1,1))+coeff* & |
|---|
| 597 | (log(t_interp(2,1,2))-log(t_interp(2,1,1))) |
|---|
| 598 | zt_interp(2,2)=log(t_interp(2,2,1))+coeff* & |
|---|
| 599 | (log(t_interp(2,2,2))-log(t_interp(2,2,1))) |
|---|
| 600 | zt_interp(1:2,1:2)=exp(zt_interp(1:2,1:2)) |
|---|
| 601 | else ! general case |
|---|
| 602 | coeff=(alt-altitude(ialt_inf))/(altitude(ialt_sup)-altitude(ialt_inf)) |
|---|
| 603 | zt_interp(1,1)=t_interp(1,1,1)+coeff*(t_interp(1,1,2)-t_interp(1,1,1)) |
|---|
| 604 | zt_interp(1,2)=t_interp(1,2,1)+coeff*(t_interp(1,2,2)-t_interp(1,2,1)) |
|---|
| 605 | zt_interp(2,1)=t_interp(2,1,1)+coeff*(t_interp(2,1,2)-t_interp(2,1,1)) |
|---|
| 606 | zt_interp(2,2)=t_interp(2,2,1)+coeff*(t_interp(2,2,2)-t_interp(2,2,1)) |
|---|
| 607 | endif |
|---|
| 608 | |
|---|
| 609 | ! 2.3 Latitudinal interpolation |
|---|
| 610 | coeff=(lat-latitude(ilat_inf))/(latitude(ilat_sup)-latitude(ilat_inf)) |
|---|
| 611 | yzt_interp(1)=zt_interp(1,1)+coeff*(zt_interp(1,2)-zt_interp(1,1)) |
|---|
| 612 | yzt_interp(2)=zt_interp(2,1)+coeff*(zt_interp(2,2)-zt_interp(2,1)) |
|---|
| 613 | |
|---|
| 614 | ! 2.4 longitudinal interpolation |
|---|
| 615 | coeff=(lon-longitude(ilon_inf))/(longitude(ilon_sup)-longitude(ilon_inf)) |
|---|
| 616 | value=yzt_interp(1)+coeff*(yzt_interp(2)-yzt_interp(1)) |
|---|
| 617 | |
|---|
| 618 | end subroutine extraction |
|---|
| [312] | 619 | |
|---|
| [282] | 620 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 621 | subroutine ls2sol(ls,sol) |
|---|
| 622 | |
|---|
| 623 | implicit none |
|---|
| 624 | ! Arguments: |
|---|
| 625 | real,intent(in) :: ls |
|---|
| 626 | real,intent(out) :: sol |
|---|
| 627 | |
|---|
| 628 | ! Local: |
|---|
| 629 | double precision xref,zx0,zteta,zz |
|---|
| 630 | !xref: mean anomaly, zteta: true anomaly, zx0: eccentric anomaly |
|---|
| 631 | double precision year_day |
|---|
| 632 | double precision peri_day,timeperi,e_elips |
|---|
| 633 | double precision pi,degrad |
|---|
| 634 | parameter (year_day=668.6d0) ! number of sols in a martian year |
|---|
| 635 | parameter (peri_day=485.35d0) ! date (in sols) of perihelion |
|---|
| 636 | ! timeperi: 2*pi*( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99 |
|---|
| 637 | parameter (timeperi=1.90258341759902d0) |
|---|
| 638 | parameter (e_elips=0.0934d0) ! eccentricity of orbit |
|---|
| 639 | parameter (pi=3.14159265358979d0) |
|---|
| 640 | parameter (degrad=57.2957795130823d0) |
|---|
| 641 | |
|---|
| 642 | if (abs(ls).lt.1.0e-5) then |
|---|
| 643 | if (ls.ge.0.0) then |
|---|
| 644 | sol = 0.0 |
|---|
| 645 | else |
|---|
| 646 | sol = real(year_day) |
|---|
| 647 | end if |
|---|
| 648 | return |
|---|
| 649 | end if |
|---|
| 650 | |
|---|
| 651 | zteta = ls/degrad + timeperi |
|---|
| 652 | zx0 = 2.0*datan(dtan(0.5*zteta)/dsqrt((1.+e_elips)/(1.-e_elips))) |
|---|
| 653 | xref = zx0-e_elips*dsin(zx0) |
|---|
| 654 | zz = xref/(2.*pi) |
|---|
| 655 | sol = real(zz*year_day + peri_day) |
|---|
| 656 | if (sol.lt.0.0) sol = sol + real(year_day) |
|---|
| 657 | if (sol.ge.year_day) sol = sol - real(year_day) |
|---|
| 658 | |
|---|
| 659 | |
|---|
| 660 | end subroutine ls2sol |
|---|
| [312] | 661 | |
|---|
