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