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