Page 1 Source Listing BUFR_TRANAMSUA 2012-11-20 14:00 tranamsua.f 1 PROGRAM BUFR_TRANAMSUA 2 C$$$ MAIN PROGRAM DOCUMENTATION BLOCK 3 C 4 C MAIN PROGRAM: BUFR_TRANAMSUA 5 C PRGMMR: KEYSER ORG: NP22 DATE: 2012-10-11 6 C 7 C ABSTRACT: Read raw AMSU-A 1B format file, decode, write selected 8 C Tb and Ta observations to unique output BUFR files. 9 C 10 C PROGRAM HISTORY LOG: 11 C 1998-06-15 Treadon -- Original author 12 C 2000-09-06 Woollen -- Added second output in BUFR 13 C 2000-11-20 Treadon -- Generalized to allow for different antenna 14 C to brightness temperature conversion constants for NOAA-15 vs. 15 C NOAA-16; Changed to properly relabel NOAA-16 satellite id from 16 C 2 to 16 to be consistent with the convention followed in the 17 C global and regional analysis systems; Added error handling when 18 C no output is created so that subsequent TRANJB's are skipped 19 C 2002-02-11 Woollen -- Modifications and corrections to output BUFR 20 C dataset: "SAID" (0-01-007) corrected to proper WMO Code Table 21 C value (was 14 for NOAA-14, etc.), "SIID" (0-02-019) repl. 22 C "SIDU" (0-02-021) which didn't seem to be correct, "HMSL" 23 C (0-07-002) corrected to proper units of meters (was being 24 C stored in km), "LSQL" (0-08-012) corrected to proper WMO Code 25 C Table value (0-land/1-sea) (was backwards), "TMBR" (0-12-163) 26 C corrected to proper units of K (was being stored as K + 27 C 273.15), channel 20 "TMBR" set to missing for HIRS-2 and HIRS-3 28 C types 29 C 2002-07-08 Keyser -- Accounts for NOAA-17 (converts NESDIS sat. 30 C no. from 6 to 17); Sets default values for Ta to Tb 31 C coefficients which are used when coefficient file is empty or 32 C missing (before program failed in this situation) 33 C 2002-10-23 Treadon -- Use lbyte instead of mbyte for unpacking 34 C counts 35 C 2004-01-23 Keyser -- Based on new namelist switch "compress", now 36 C has option to write compressed BUFR messages using WRITCP 37 C instead of WRITSB (removes the need for the downstream program 38 C BUFR_COMPRESS) 39 C 2005-04-29 Keyser -- Added processing of Ta reports to new, unique 40 C BUFR file in unit 53 41 C 2005-06-21 Keyser -- Modified to handle processing of NOAA-18 data 42 C 2006-04-21 Derber -- Modified to estimate solar and satellite 43 C azimuth (via new subroutine SATAZIMUTH), and to update time 44 C within a scan line 45 C 2007-02-09 Keyser -- Modified to encode the following new 46 C information into output BUFR file: estimated solar azimuth 47 C (SOLAZI) and estimated satellite azimuth (BEARAZ) for each 48 C subset (retrieval), and cold space temperature correction (CSTC) 49 C for each channel in subset (note: this is not added to output 50 C IEEE file); the "report" time in the BUFR and IEEE files now 51 C varies across a scan line as a result of 2006-04-21 change; 52 C modified subr. DATTIM to input real double precision second-of- 53 C day and output real double precision second-of-minute (both had 54 C been integer); array passed into subr. BUFR1B now contains spot 55 C hour-of-day, minute-of-hour and second-of-minute rather than 56 C only second-of-day (since the hour, minute and second are stored 57 C in BUFR, second now rounded to nearest whole second rather than Page 2 Source Listing BUFR_TRANAMSUA 2012-11-20 14:00 tranamsua.f 58 C truncated) (IEEE files still store second-of-day); increased 59 C limit for i/o filename length from 80 characters to 500 60 C characters; modified to write BUFR code table value 0-01-007, 61 C the BUFR value for satellite id, into word 1 of output "bdata" 62 C array, rather than the actual satellite number as before (this 63 C simplifies subroutine BUFR1B which then encodes this value 64 C directly into BUFR) (note: the actual satellite number is still 65 C written to the output IEEE file); now accounts for new METOP-1 66 C satellite 67 C 2007-04-12 Keyser -- Modified to correct METOP satellite number 68 C to METOP-2 with BUFR satellite id = 4 (was incorrectly set to 69 C METOP-1 with BUFR satellite id = 3) 70 C 2009-05-08 Krasowski -- Modified to handle processing of NOAA-19 71 C data 72 C 2012-10-11 Keyser -- Changes to run on WCOSS. Modified to handle 73 C processing of METOP-1 data. Removed IEEE output processing 74 C 75 C USAGE: 76 C 77 C INPUT FILES: 78 C UNIT 05 - Standard input (namelist "input") 79 C UNIT 11 - Binary file containing raw 1B AMSU-A data 80 C UNIT 12 - BUFR mnemonic table 81 C UNIT 21 - Binary file containing Ta to Tb coefficients 82 C UNIT 41 - Binary file containing topography information 83 C used by function LANSEA 84 C 85 C ***NOTE*** 86 C Function LANSEA assumes this information is in 87 C a file named 'lowtopog.dat' which is local to 88 C the working directory. 89 C 90 C OUTPUT FILES: 91 C UNIT 06 - Printout 92 C UNIT 52 - BUFR file containing AMSU-A Tb data 93 C UNIT 53 - BUFR file containing AMSU-A Ta data 94 C 95 C SUBPROGRAMS CALLED: 96 C UNIQUE: - AMSUA CHARS DATTIM ICHARS LANSEA LBIT 97 C MBYTE LBYTE XFLOAT SATAZIMUTH BUFR1B 98 C SYSTEM: - SYSTEM 99 C LIBRARY: 100 C W3NCO - W3TAGB W3TAGE ERREXIT W3FS26 W3MOVDAT W3DOXDAT 101 C BUFRLIB - OPENBF CLOSBF OPENMB WRITSB WRITCP UFBSEQ 102 C MESGBC 103 C 104 C 105 C EXIT STATES 106 C 0 = No errors detected 107 C 1 = Data type id decoded from header is not for AMSU-A 108 C 2 = Problem reading Ta to Tb coefficient file 109 C 3 = Problem reading header record of 1b AMSU-A file 110 C 6 = Unknown satellite id 111 C 7 = Unknown satellite instrument 112 C 113 C REMARKS: 114 C Switches read in Namelist INPUT: Page 3 Source Listing BUFR_TRANAMSUA 2012-11-20 14:00 tranamsua.f 115 C INFILE - Path to input 1B data file 116 C COMPRESS - BUFR compression switch (YES or NO) 117 C COEFILE - Path to input coefficient file 118 C PROCESS_Tb - Process brightness temps into BUFR files? 119 C PROCESS_Ta - Process antenna temps into BUFR files? 120 C 121 C#################################################################### 122 C NOTE: This program is designed to process BOTH Ta and Tb into 123 C separate BUFR files in a single run. However, BUFRLIB 124 C routine WRITCP (which writes COMPRESSED messages) cannot 125 C operate on two output files at the same time. If it is 126 C ever modified to do so (like WRITSB which writes uncompressed 127 C BUFR messages), then the code to ready to handle this. In 128 C the meantime, this code must be executed twice to process 129 C both Ta and Tb (once with process_Ta=YES and process_Tb=NO 130 C and again with process_Ta=NO and process_Tb=YES). 131 C#################################################################### 132 C 133 C ATTRIBUTES: 134 C LANGUAGE: FORTRAN 90 135 C MACHINE: NCEP WCOSS 136 C 137 C$$$ 138 139 C Declare namelist variables and namelist 140 C --------------------------------------- 141 142 integer stdout 143 character*500 infile,coefile 144 character*8 compress,process_Tb,process_Ta 145 namelist /input/ infile,coefile,compress,process_Tb,process_Ta 146 147 common/switches/compress,process_Tb,process_Ta 148 149 C Set I/O unit numbers 150 C -------------------- 151 152 data lunam, stdout / 5, 6 / 153 data lunin, luncof / 11, 21 / 154 155 call w3tagb('BUFR_TRANAMSUA',2012,0265,0068,'NP22') 156 157 print * 158 print *, 'WELCOME TO BUFR_TRANAMSUA - Version 10/11/2012' 159 print * 160 161 C Get Namelist input 162 C ------------------ 163 164 read(lunam,input) 165 write(stdout,*)'namelist input below' 166 write(stdout,input) 167 168 C Read/decode/output data records scan by scan 169 C -------------------------------------------- 170 171 call amsua(lunin,infile,luncof,coefile) Page 4 Source Listing BUFR_TRANAMSUA 2012-11-20 14:00 tranamsua.f 172 173 call w3tage('BUFR_TRANAMSUA') 174 175 stop 176 end ENTRY POINTS Name MAIN__ SYMBOL CROSS REFERENCE Name Object Declared Type Bytes Dimen Elements Attributes References AMSUA Subr 171 171 BUFR_TRANAMSUA Prog 1 COEFILE Local 143 CHAR 500 scalar 145,171 COMPRESS Scalar 144 CHAR 8 scalar COM 145 INFILE Local 143 CHAR 500 scalar 145,171 INPUT Local 145 scalar 164,166 LUNAM Local 152 I(4) 4 scalar 152,164 LUNCOF Local 153 I(4) 4 scalar 153,171 LUNIN Local 153 I(4) 4 scalar 153,171 PROCESS_TA Scalar 144 CHAR 8 scalar COM 145 PROCESS_TB Scalar 144 CHAR 8 scalar COM 145 STDOUT Local 142 I(4) 4 scalar 152,165,166 SWITCHES Common 147 24 W3TAGB Subr 155 155 W3TAGE Subr 173 173 Page 5 Source Listing BUFR_TRANAMSUA 2012-11-20 14:00 tranamsua.f 177 178 SUBROUTINE AMSUA(LUNIN,RAWAMSU,LUNCOF,COEFILE) 179 C$$$ SUBPROGRAM DOCUMENTATION BLOCK 180 C 181 C SUBPROGRAM: AMSUA 182 C PRGMMR: KEYSER ORG: NP22 DATE: 2012-10-11 183 C 184 C ABSTRACT: Read raw AMSU-A 1B format file, decode, write selected 185 C Tb and Ta observations to unique output BUFR files. 186 C 187 C PROGRAM HISTORY LOG: 188 C 1998-06-15 Treadon -- Original author 189 C 2002-07-08 Keyser -- Accounts for NOAA-17 (converts NESDIS sat. 190 C no. from 6 to 17) 191 C 2004-05-03 Keyser -- Added processing of Ta reports to new, unique 192 C BUFR file in UNIT 53 193 C 2005-06-21 Keyser -- Modified to handle processing of NOAA-18 data 194 C 2006-04-21 Derber -- Modified to estimate solar and satellite 195 C azimuth (via new subroutine SATAZIMUTH), and to update time 196 C within a scan line 197 C 2006-07-20 Keyser -- Modified to encode the following new 198 C information into output BUFR file: estimated solar azimuth 199 C (SOLAZI) and estimated satellite azimuth (BEARAZ) for each 200 C subset (retrieval), and cold space temperature correction (CSTC) 201 C for each channel in subset (note: this is not added to output 202 C IEEE file); the "report" time in the BUFR and IEEE files now 203 C varies across a scan line as a result of 2006-04-21 change; 204 C increased limit for i/o filename length from 80 characters to 205 C 500 characters 206 C 2009-05-08 Krasowski -- Modified to handle processing of NOAA-19 207 C data 208 C 2012-10-11 Keyser -- Removed IEEE output processing 209 C 210 C USAGE: CALL AMSUA(LUNIN,RAWAMSU,LUNCOF,COEFILE) 211 C INPUT ARGUMENT LIST: 212 C LUNIN - Unit connected to raw 1B AMSU-A data file 213 C RAWAMSU - Name of raw 1B AMSU-A data file 214 C LUNCOF - Unit connected to Ta to Tb conversion coefficients 215 C 216 C INPUT FILES: 217 C UNIT 12 - BUFR mnemonic table 218 C UNIT 41 - Binary file containing topography information 219 C used by function LANSEA 220 C UNIT LUNIN - Binary file containing raw 1B AMSU-A data 221 C UNIT LUNCOF - Binary file containing Ta to Tb coefficients 222 C 223 C ***NOTE*** 224 C Function LANSEA assumes this information is in 225 C a file named 'lowtopog.dat' which is local to 226 C the working directory. 227 C 228 C OUTPUT FILES: 229 C UNIT 06 - Printout 230 C UNIT 52 - BUFR file containing AMSU-A Tb data 231 C UNIT 53 - BUFR file containing AMSU-A Ta data 232 C 233 C REMARKS: Page 6 Source Listing AMSUA 2012-11-20 14:00 tranamsua.f 234 C 235 C ATTRIBUTES: 236 C LANGUAGE: FORTRAN 90 237 C MACHINE: NCEP WCOSS 238 C 239 C$$$ 240 241 C Include machine dependent parameters 242 C ------------------------------------ 243 244 include 'rfac.inc' 245 246 C Declare/set parameters: 247 C NBYTE1 = Total number of bytes (2560) in AMSU-A data record 248 C NBYTE4 = Number of 4-byte words in NBYTE1 bytes (2560/4=640) 249 C NSET = Number of topography datasets for function LANSEA3 250 C (not currently used) 251 C EPS = A "small" number 252 C MCH = Number of channels 253 C MPOS = Number of spots (positions) on a scan line 254 271 integer,parameter::real_32=selected_real_kind(6,37) 272 integer,parameter::real_64=selected_real_kind(15,307) 273 real(real_64) eps 274 parameter (nbyte1=2560,nbyte4=nbyte1/4) 275 parameter (nset=3) 276 parameter (eps=1.d-12) 277 parameter (mch=15) 278 parameter (mpos=30) 279 280 C Set parameters for structure of output data file 281 C ------------------------------------------------ 282 283 parameter (nreal=18,nperchan=2,ntot=nreal+nperchan*mch) 284 285 C Declare variables 286 C ----------------- 287 288 character*1 kbuf(nbyte1) 289 character*4 indat(nbyte4),jbuf(nbyte4) 290 character*8 compress,process_Tb,process_Ta 291 character*40 mapfile(nset) 292 character*500 rawamsu,coefile 293 294 integer stdout 295 integer(8) itime 296 integer idt(2),ndt(5),idat(8),jdat(8) 297 integer ichan(mch),lndsea(mpos),ikeepb(mpos),ikeepa(mpos) 298 integer imx(nset),jmx(nset),ibadc(mch) 299 300 real(real_64) p1,p2,term1,term2,term3,ta(mch,mpos),b,c,clight,ta0 301 real(real_64) scale,scale5,scale6,scale9,scale13,scale19 302 real(real_64) sctime,xsec,counts,rads,rad0,sathgt,rnorm,tb0 303 real(real_64) tshelf11a,tshelf12,tshelf2 304 real(real_64) soza(mpos),saza(mpos),rlocaz(mpos),aazimuth(mpos) 305 real(real_64) slat(mpos),slon(mpos),sazimuth(mpos) 306 real(real_64) cwave(mch),cnst1(mch),cnst2(mch) Page 7 Source Listing AMSUA 2012-11-20 14:00 tranamsua.f 307 real(real_64) cfrq0(mch) 308 real(real_64) rad(mch,mpos),tb(mch,mpos),sfchgt(mpos) 309 real(real_64) c0(mch),c1(mch),c2(mch) 310 real(real_32) bdata(ntot),adata(ntot),rinc(5),sctime1 311 real(real_64) f0(mpos,mch),f1(mpos,mch),f2(mpos,mch) 312 real(real_64) eta(mch),tref(mch),badr(mch),badtb(mch),dt(mch) 313 real(real_64) badta(mch) 314 315 C Declare equivalences 316 C -------------------- 317 318 equivalence (kbuf(1),jbuf(1)) 319 320 common/switches/compress,process_Tb,process_Ta 321 322 C Set information for different resolution map datasets 323 C (NOTE: Currently we do not use this information) 324 C ----------------------------------------------------- 325 326 data imx / 360, 720, 1440 / 327 data jmx / 181, 361, 721 / 328 data mapfile / 'mapdat_100.iee', 'mapdat_050.iee', 329 x 'mapdat_025.iee' / 330 331 C Lower/upper limits for gross temperature check on Tb 332 C ---------------------------------------------------- 333 334 data tlo,thi / 100., 400. / 335 336 C Constants for Planck equation 337 C ----------------------------- 338 339 data p1,p2 / 1.1910659d-5, 1.438833d0 / 340 341 C Speed of light (m/s) and center frequency of AMSU-A channels (GHz) 342 C ------------------------------------------------------------------ 343 344 data clight / 2.99793d8 / 345 data cfrq0 / 23.8d0, 31.4d0, 50.3d0, 52.8d0, 53.596d0, 54.50d0, 346 x 54.94d0, 55.50d0, 6*57.290344d0, 89.0d0 / 347 348 C Set constants for Ta to Tb conversion 349 C ------------------------------------- 350 351 data eta /.01d0,.08d0,.03d0,.04d0,.04d0,.03d0,.03d0,7*.04d0,.11d0/ 352 353 C Missing data flag 354 C ----------------- 355 356 data rmiss / -999. / 357 358 C Channel numbers 359 C --------------- 360 361 data ichan / 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 362 x 12, 13, 14, 15 / 363 Page 8 Source Listing AMSUA 2012-11-20 14:00 tranamsua.f 364 365 C Set I/O unit numbers (including standard output) 366 C ------------------------------------------------ 367 368 data stdout / 6/ 369 data lundx /12/ 370 data lubfrb /52/, lubfra /53/ 371 372 write(stdout,*)' ' 373 write(stdout,*)' BEGIN AMSU-A 1B DECODE' 374 375 C Initialize arrays 376 C ----------------- 377 378 badr = 0. 379 badtb = 0. 380 badta = 0. 381 nprint = 500 ! skip between data record diagnostic prints 382 383 C Write header record to standard output 384 C -------------------------------------- 385 386 write(stdout,*)' ' 387 write(stdout,*)'header information below' 388 write(stdout,*)'nreal,mch = ',nreal,mch 389 write(stdout,*)'ntot = ',ntot 390 write(stdout,*)'channel numbers below' 391 write(stdout,*) (ichan(i),i=1,mch) 392 write(stdout,*)' ' 393 ccccc write(99,*) nreal-4,mch,(ichan(i),i=1,mch) 394 395 C Open output BUFR files 396 C ---------------------- 397 398 ccccc call openbf(lubfrb,'OUT',lundx) 399 call openbf(lubfrb,'NODX',lundx) 400 ccccc call openbf(lubfra,'OUT',lundx) 401 call openbf(lubfra,'NODX',lundx) 402 403 C Load arrays containing coefficients to convert the antenna 404 C temperature (Ta) to brightness temperature (Tb) - finish loading 405 C after read, then scale coefficients 406 C ----------------------------------------------------------------- 407 408 open(luncof,file=coefile) 409 410 C Default values for when file is missing 411 C --------------------------------------- 412 413 do n = 1,mch 414 do i = 1,mpos 415 f0(i,n)=100.0 416 f1(i,n)=0.0 417 f2(i,n)=0.0 418 end do 419 end do 420 do i = 1,mpos Page 9 Source Listing AMSUA 2012-11-20 14:00 tranamsua.f 421 read(luncof,*,err=1800,end=100) (f0(i,n),f1(i,n),f2(i,n),n=1,5) 422 end do 423 do i = 1,mpos 424 read(luncof,*,err=1800,end=100)(f0(i,n),f1(i,n),f2(i,n),n=6,10) 425 end do 426 write(stdout,*)'Read Ta to Tb conversion coefficients from ', 427 & 'coefficient file' 428 429 C If present, read eta factor from coefficient file 430 C ------------------------------------------------- 431 432 read(luncof,*,end=101) 433 write(stdout,*)'Read eta factor from coefficient file' 434 read(luncof,*) (eta(n),n=1,9),eta(15) 435 do j=10,14 436 eta(j)=eta(9) 437 end do 438 go to 102 439 100 continue 440 write(stdout,*)'Ta to Tb conversion coefficient file missing, ', 441 & 'use default values --> no Ta to Tb conversion' 442 101 continue 443 write(stdout,*)'Eta factor not found in coefficient file, use ', 444 & 'default values --> no Ta to Tb conversion' 445 102 continue 446 close(luncof) 447 448 do i = 1,mpos 449 f0(i,15) = f0(i,10) 450 f1(i,15) = f1(i,10) 451 f2(i,15) = f2(i,10) 452 end do 453 do i = 1,mpos 454 do j = 10,14 455 f0(i,j) = f0(i,9) 456 f1(i,j) = f1(i,9) 457 f2(i,j) = f2(i,9) 458 end do 459 end do 460 do i = 1,mpos 461 do j = 1,mch 462 f0(i,j) = f0(i,j)/100. 463 f1(i,j) = f1(i,j)/100. 464 f2(i,j) = f2(i,j)/100. 465 end do 466 end do 467 468 C Open unit to raw 1B AMSU-A data file - read header record, see if 469 C valid data type - if not, exit routine 470 C ----------------------------------------------------------------- 471 472 open(lunin,file=rawamsu,recl=nbyte1/rfac, 473 & access='direct',status='old') 474 nri = 1 475 read (lunin,rec=nri,err=1900) (kbuf(i),i=1,nbyte1) 476 477 C Load header record into work array Page 10 Source Listing AMSUA 2012-11-20 14:00 tranamsua.f 478 C ---------------------------------- 479 480 do i = 1,nbyte4 481 indat(i) = jbuf(i) 482 end do 483 484 C Extract NOAA spacecraft identification code (72*8+1=577) 485 C and convert it into BUFR value (CODE TABLE 0-01-007) 486 C -------------------------------------------------------- 487 488 jsat = lbyte(577,16,indat) 489 if (jsat.eq.4) then ! NOAA-15 490 jsat0 = jsat 491 jsat = 206 492 write(stdout,*) '***WARNING: reset NOAA-15 satellite id ', 493 x 'from ',jsat0,' to ',jsat 494 elseif (jsat.eq.2) then ! NOAA-16 495 jsat0 = jsat 496 jsat = 207 497 write(stdout,*) '***WARNING: reset NOAA-16 satellite id ', 498 x 'from ',jsat0,' to ',jsat 499 elseif (jsat.eq.6) then ! NOAA-17 500 jsat0 = jsat 501 jsat = 208 502 write(stdout,*) '***WARNING: reset NOAA-17 satellite id ', 503 x 'from ',jsat0,' to ',jsat 504 elseif (jsat.eq.7) then ! NOAA-18 505 jsat0 = jsat 506 jsat = 209 507 write(stdout,*) '***WARNING: reset NOAA-18 satellite id ', 508 x 'from ',jsat0,' to ',jsat 509 elseif (jsat.eq.8) then ! NOAA-19 510 jsat0 = jsat 511 jsat = 223 512 write(stdout,*) '***WARNING: reset NOAA-19 satellite id ', 513 x 'from ',jsat0,' to ',jsat 514 else if(jsat.eq.12) then ! METOP-2 515 jsat0 = jsat 516 jsat = 4 517 write(stdout,*) '***WARNING: reset METOP-2 satellite id ', 518 x 'from ',jsat0,' to ',jsat 519 else if(jsat.eq.11) then ! METOP-1 520 jsat0 = jsat 521 jsat = 3 522 write(stdout,*) '***WARNING: reset METOP-1 satellite id ', 523 x 'from ',jsat0,' to ',jsat 524 else 525 write(stdout,*) '***ERROR*** unknown satellite id ',jsat 526 call w3tage('BUFR_TRANAMSUA') 527 call errexit(6) 528 endif 529 530 C Extract data type code (76*8+1=609) 531 C ----------------------------------- 532 533 jtype = lbyte(609,16,indat) 534 Page 11 Source Listing AMSUA 2012-11-20 14:00 tranamsua.f 535 if (jtype.ne.10) then 536 write(stdout,*)'***ERROR*** Input data file does not contain', 537 x ' AMSU-A data (type=10). data type = ',jtype 538 call w3tage('BUFR_TRANAMSUA') 539 call errexit(1) 540 endif 541 write(stdout,*) 'Data type and BUFR satellite id = ',jtype,jsat 542 543 C Extract number of data records in data set (144*8+1=1153) 544 C and number of scans (146*8+1=1169) 545 C --------------------------------------------------------- 546 547 nrecs = lbyte(1153,16,indat) 548 nscan = lbyte(1169,16,indat) 549 write(stdout,*)'nrecs,nscan=',nrecs,nscan 550 551 C Extract coefficients for radiance to temperature conversion 552 C ----------------------------------------------------------- 553 554 scale6 = 1.d-6 555 do j = 1,mch 556 jb = 689 + (j-1)*12 557 jb0 = jb 558 jb1 = jb0 + 4 559 jb2 = jb1 + 4 560 cwave(j) = xfloat(1,kbuf(jb0))*scale6 561 cnst1(j) = xfloat(1,kbuf(jb1))*scale6 562 cnst2(j) = xfloat(1,kbuf(jb2))*scale6 563 end do 564 write(stdout,*)'chn 1 cwave,cnst1,cnst2=', 565 x cwave(1),cnst1(1),cnst2(1) 566 567 C Extract instrument temperature and load into reference array 568 C ------------------------------------------------------------ 569 570 scale = 1.d-2 571 tshelf11a = mbyte(1745,16,indat)*scale ! (218*8+1=1745) 572 tshelf12 = mbyte(1793,16,indat)*scale ! (224*8+1=1793) 573 tshelf2 = mbyte(1841,16,indat)*scale ! (230*8+1=1841) 574 575 C A2 antenna handles channels 1,2 576 C ------------------------------- 577 578 tref(1) = tshelf2 579 tref(2) = tshelf2 580 581 C A1-1 antenna handles channels 6,7,9-15 582 C -------------------------------------- 583 584 tref(6) = tshelf11a 585 tref(7) = tshelf11a 586 do i = 9,15 587 tref(i) = tshelf11a 588 end do 589 590 C A1-2 antenna handles channels 3,4,5,8 591 C ------------------------------------- Page 12 Source Listing AMSUA 2012-11-20 14:00 tranamsua.f 592 593 tref(3) = tshelf12 594 tref(4) = tshelf12 595 tref(5) = tshelf12 596 tref(8) = tshelf12 597 598 C QC instrument temperatures - the lower and upper bounds are 599 C aribtrary - if unreasonable, replace with 290K (the instrument 600 C temperature s/b is near 290K) 601 C --------------------------------------------------------------- 602 603 do i = 1,mch 604 if (tref(i).lt.270 .or. tref(i).gt.310) tref(i) = 290. 605 end do 606 607 C Extract cold space fixed bias correction 608 C ---------------------------------------- 609 610 scale = 1.d-3 611 do i = 1,mch 612 jb = (270 + (i-1)*8)*8 + 1 613 dtcold = mbyte(jb,16,indat)*scale 614 dt(i) = dtcold 615 end do 616 617 write(stdout,*)'instrument T=',(tref(i),i=1,mch) 618 write(stdout,*)'cold space dt=',(dt(i),i=1,mch) 619 write(stdout,*)' ' 620 621 C Prepatory initializations prior to reading in satellite data 622 C ------------------------------------------------------------ 623 624 nopos = 0 625 nqcbad = 0 626 nqctim = 0 627 nqccal = 0 628 nqcloc = 0 629 nermin = 0 630 nermaj = 0 631 nbadc = 0 632 nbadl = 0 633 nbadtb = 0 634 nbadta = 0 635 nbadr = 0 636 nrecb = 0 637 nreca = 0 638 nrepb = 0 639 nrepa = 0 640 nskipc = 0 641 nskipq = 0 642 nskipm = 0 643 nskiptb = 0 644 nskipta = 0 645 nlandb = 0 646 nlanda = 0 647 nseab = 0 648 nseaa = 0 Page 13 Source Listing AMSUA 2012-11-20 14:00 tranamsua.f 649 nlo = 0 650 651 1200 continue 652 653 C********************************************************************** 654 C MAIN LOOP OVER NUMBER OF SCANS 655 C********************************************************************** 656 657 nri = nri + 1 ! Increment record counter 658 659 C Read data record and load into local work array 660 C ----------------------------------------------- 661 662 read(lunin,rec=nri,err=1600) (kbuf(i),i=1,nbyte1) 663 664 do i = 1,nbyte4 665 indat(i) = jbuf(i) 666 end do 667 668 nlo = nlo + 1 ! Increment line counter 669 line = nlo 670 671 C Extract scan line number, start date/time, position, and type 672 C ------------------------------------------------------------- 673 674 iline = lbyte(1,16,indat) 675 iyear = lbyte(17,16,indat) ! (2*8+1=17) 676 iddd = lbyte(33,16,indat) ! (4*8+1=33) 677 itime = lbyte(65,32,indat) ! (8*8+1=65) 678 idt(1) = iddd ! day of the year 679 idt(2) = iyear ! 4-digit year 680 sctime = 1.d-3*itime ! second of the day 681 682 C Convert scan start time from year, day-of-year and second-of-day 683 C to YYYY,MM,DD,HH,mm,ss 684 C ---------------------------------------------------------------- 685 686 call dattim(idt,sctime,ndt,xsec) 687 rinc = 0 688 idat = 0 689 idat(1:3) = ndt(1:3) 690 idat(5:6) = ndt(4:5) 691 idat(7) = xsec 692 idat(8) = mod(xsec*1000._8,1000._8) 693 694 C Extract quality bits - if all good (=0) continue, else skip this scan 695 C --------------------------------------------------------------------- 696 697 isum = 0 698 iqcbad = lbyte(193,1,indat) ! (8*24+1=193) 699 iqctim = lbyte(194,1,indat) ! (8*24+1+1=194) 700 iqccal = lbyte(196,1,indat) ! (8*24+1+3=196) 701 iqcloc = lbyte(197,1,indat) ! (8*24+1+4=197) 702 iermin = lbyte(222,1,indat) ! (8*24+1+29=222) 703 iermaj = lbyte(223,1,indat) ! (8*24+1+30=223) 704 isum = iqcbad + iqctim + iqccal + iqcloc 705 x + iermin + iermaj Page 14 Source Listing AMSUA 2012-11-20 14:00 tranamsua.f 706 if (iqcbad.ne.0) nqcbad = nqcbad + 1 707 if (iqctim.ne.0) nqctim = nqctim + 1 708 if (iqccal.ne.0) nqccal = nqccal + 1 709 if (iqcloc.ne.0) nqcloc = nqcloc + 1 710 if (iermin.ne.0) nermin = nermin + 1 711 if (iermaj.ne.0) nermaj = nermaj + 1 712 if (isum.ne.0) then 713 nskipq = nskipq + 1 714 write(stdout,1000) nlo,iline,iqcbad,iqctim, 715 x iqccal,iqcloc,iermin,iermaj,nskipq 716 1000 format('***FAIL QC : nlo=',i6,' iline=',i6, 717 x ' bad=',i2,' time=',i2,' cali=',i2,' loc=',i2, 718 x ' min=',i2,' maj=',i2,' nskipq=',i6) 719 goto 1200 720 endif 721 722 C Extract calibration quality bits for each channel 723 C ------------------------------------------------- 724 725 isum = 0 726 do j = 1,mch 727 ibadc(j) = 0 728 jb = (32 + (j-1)*2)*8+1 729 iqccal = lbyte(jb,16,indat) 730 ibadc(j) = iqccal 731 isum = isum + ibadc(j) 732 end do 733 if (isum.ne.0) then 734 nskipc = nskipc + 1 735 write(stdout,1005) nlo,iline,(ibadc(j),j=1,mch) 736 1005 format('***FAIL CAL : nlo=',i6,' iline=',i6, 737 x ' badcal=',20(i2,1x)) 738 endif 739 740 C Extract calibration coefficients 741 C -------------------------------- 742 743 scale19 = 1.d-19 744 scale13 = 1.d-13 745 scale9 = 1.d-9 746 do j = 1,mch 747 jb2 = 81 + (j-1)*12 748 c2(j) = xfloat(1,kbuf(jb2))*scale19 749 jb1 = jb2 + 4 750 c1(j) = xfloat(1,kbuf(jb1))*scale13 751 jb0 = jb1 + 4 752 c0(j) = xfloat(1,kbuf(jb0))*scale9 753 754 C Code to pull out secondary calibration coefficients 755 C --------------------------------------------------- 756 757 ccccc jb2 = 261 + (j-1)*12 758 ccccc c2j = xfloat(1,kbuf(jb2))*scale19 759 ccccc jb1 = jb2 + 4 760 ccccc c1j = xfloat(1,kbuf(jb1))*scale13 761 ccccc jb0 = jb1 + 4 762 ccccc c0j = xfloat(1,kbuf(jb0))*scale9 Page 15 Source Listing AMSUA 2012-11-20 14:00 tranamsua.f 763 end do 764 765 C EXTRACT NAVIGATION DATA 766 C ----------------------- 767 C ----------------------- 768 769 C Extract spacecraft altitude (km) 770 C -------------------------------- 771 772 scale = 1.d-1 773 jb = 470*8+1 774 sathgt = lbyte(jb,16,indat)*scale 775 776 C Extract angular relationships 777 C ----------------------------- 778 779 scale = 1.d-2 780 do i = 1,mpos 781 jb0 = 472*8+1 + (i-1)*48 782 jb1 = jb0 + 16 783 jb2 = jb1 + 16 784 soza(i) = mbyte(jb0,16,indat)*scale 785 saza(i) = mbyte(jb1,16,indat)*scale 786 rlocaz(i) = mbyte(jb2,16,indat)*scale 787 end do 788 789 C Extract earth location data 790 C --------------------------- 791 792 scale = 1.d-4 793 do i = 1,mpos 794 jb0 = 653 + (i-1)*8 795 slat(i) = xfloat(1,kbuf(jb0))*scale 796 jb1 = jb0 + 4 797 slon(i) = xfloat(1,kbuf(jb1))*scale 798 lndsea(i) = rmiss 799 sfchgt(i) = rmiss 800 ikeepb(i) = 0 801 ikeepa(i) = 0 802 if ( (abs(slat(i)).gt.90.) .or. 803 x (abs(slon(i)).gt.180.) ) then 804 ikeepb(i) = 0 805 ikeepa(i) = 0 806 nbadl = nbadl + 1 807 elseif ( (abs(slat(i)).le.eps) .and. 808 x (abs(slon(i)).le.eps) ) then 809 ikeepb(i) = 0 810 ikeepa(i) = 0 811 nopos = nopos + 1 812 else 813 ikeepb(i) = 1 814 ikeepa(i) = 1 815 816 C Set surface type information based on resolution option 817 C If iresol = 1, use 1.0 degree dataset 818 C iresol = 2, use 0.5 degree dataset 819 C iresol = 3, use 0.25 degree dataset Page 16 Source Listing AMSUA 2012-11-20 14:00 tranamsua.f 820 C ------------------------------------------------------- 821 822 ccccc call lansea3(xlat,xlon,imx(iresol),jmx(iresol), 823 ccccc x mapfile(iresol),rmask,water,elev,stdev) 824 ccccc lndsea(i) = rmask + 1.d-3 825 ccccc sfchgt(i) = elev 826 827 ils = lansea(slat(i),slon(i),ll) 828 if (ils.eq.2) then 829 lndsea(i) = 0 830 sfchgt(i) = 0.0 831 elseif (ils.eq.1) then 832 lndsea(i) = 1 833 sfchgt(i) = 1.*ll 834 else 835 lndsea(i) = rmiss 836 sfchgt(i) = rmiss 837 endif 838 endif 839 end do 840 841 C Is scan line for full scan mode - if so, keep, else skip this scan 842 C ------------------------------------------------------------------ 843 844 imode1 = lbyte(7207,1,indat) ! (8*900+1+6=7207) 845 imode2 = lbyte(17511,1,indat) ! (8*2188+1+6=17511) 846 if ( imode1.ne.1 .or. imode2.ne.1) then 847 nskipm = nskipm + 1 848 write(stdout,1010) nlo,iline,imode1,imode2,nskipm 849 1010 format('***FAIL MODE: nlo=',i6,' iline=',i6, 850 x ' imode1=',i2,' imode2=',i2,' nskipm=',i6) 851 goto 1200 852 endif 853 854 C Extract AMSU-A counts for channels 3-15 and 1-2, then convert counts 855 C to radiances 856 C -------------------------------------------------------------------- 857 858 do i = 1,mpos 859 860 C Channels 3-15 861 C ------------- 862 863 jb = (904 + (i-1)*34 + 4*2)*8+1 864 do j = 3,mch 865 jb0 = jb + (j-3)*16 866 counts = lbyte(jb0,16,indat) 867 rads = c0(j) + (c1(j)+c2(j)*counts)*counts 868 if (rads.lt.0.) then 869 nbadr = nbadr + 1 870 badr(j) = badr(j) + 1 871 rads = rmiss 872 endif 873 rad(j,i) = rads 874 end do 875 876 C Channels 1-2 Page 17 Source Listing AMSUA 2012-11-20 14:00 tranamsua.f 877 C ------------ 878 879 jb = (2192 + (i-1)*8 + 2*2)*8+1 880 do j = 1,2 881 jb0 = jb + (j-1)*16 882 counts = lbyte(jb0,16,indat) 883 rads = c0(j) + (c1(j)+c2(j)*counts)*counts 884 if (rads.lt.0.) then 885 nbadr = nbadr + 1 886 badr(j) = badr(j) + 1 887 rads = rmiss 888 endif 889 rad(j,i) = rads 890 end do 891 end do 892 893 C----------------------------------------------------------------------- 894 C Convert radiances to antenna temperature (Ta), then convert antenna 895 C temperature (Ta) to brightness temperature (Tb) 896 C QC all channels - if all channels are bad for a given spot, set 897 C flag to omit data in final write (for both Tb and Ta) 898 C----------------------------------------------------------------------- 899 900 do i = 1,mpos 901 ibadtb = 0 902 ibadta = 0 903 do j = 1,mch 904 rads = rad(j,i) 905 term1 = p2*cwave(j) 906 if (rads.gt.eps) then 907 term2 = 1. + p1*cwave(j)**3/rads 908 if (term2.le.0) term2 = eps 909 term3 = log(term2) 910 ta0 = term1/term3 911 b = cnst1(j) 912 c = cnst2(j) 913 ta(j,i) = (ta0-b)/c 914 rnorm = eta(j)*f1(i,j) + f2(i,j) + f0(i,j) 915 tb(j,i) = (rnorm*ta(j,i) - eta(j)*f1(i,j)*tref(j) - 916 x f2(i,j)*(2.73+dt(j)))/f0(i,j) 917 tb0 = tb(j,i) 918 else 919 ta(j,i) = rmiss 920 tb(j,i) = rmiss 921 tb0 = rmiss 922 endif 923 924 C Apply gross check to Tb using limits TLO and THI set in data stmt 925 C ----------------------------------------------------------------- 926 927 if ( (tb(j,i).lt.tlo) .or. 928 x (tb(j,i).gt.thi) ) then 929 nbadtb = nbadtb + 1 930 badtb(j) = badtb(j) + 1 931 tb(j,i) = rmiss 932 endif 933 Page 18 Source Listing AMSUA 2012-11-20 14:00 tranamsua.f 934 C Apply gross check to Ta - only limit is for Ta > 1000 935 C ----------------------------------------------------- 936 937 if (ta(j,i).gt.1000.) then 938 nbadta = nbadta + 1 939 badta(j) = badta(j) + 1 940 ta(j,i) = rmiss 941 endif 942 943 C If calibration quality flag for this channel is nonzero, we do not 944 C want to use this channel for Tb 945 C ------------------------------------------------------------------ 946 947 if (ibadc(j).ne.0) then 948 nbadc = nbadc + nbadc 949 tb(j,i) = rmiss 950 endif 951 952 C Count number of bad channels for current scan position 953 C ------------------------------------------------------ 954 955 if (tb(j,i).lt.0.) ibadtb = ibadtb + 1 956 if (ta(j,i).lt.0.) ibadta = ibadta + 1 957 958 end do 959 960 C If all Tb channels are bad for current scan position, set keep flag 961 C to zero (this tells the code below to not write Tb for this spot 962 C to the output files) 963 C ------------------------------------------------------------------- 964 965 if (ibadtb.eq.mch) then 966 nskiptb = nskiptb + 1 967 ikeepb(i) = 0 968 endif 969 970 C If all Ta channels are bad for current scan position, set keep flag 971 C to zero (this tells the code below to not write Ta for this spot 972 C to the output BUFR file) 973 C ------------------------------------------------------------------- 974 975 if (ibadta.eq.mch) then 976 nskipta = nskipta + 1 977 ikeepa(i) = 0 978 endif 979 980 end do 981 982 C ------------------------------------------------------------- 983 C WRITE AMSU-A DATA FOR EACH SPOT POSITION ON CURRENT SCAN LINE 984 C ------------------------------------------------------------- 985 986 do i = 1,mpos 987 988 if (min(ikeepb(i),ikeepa(i)).eq.1) then 989 990 C First, update scan start date time to reflect date time for spot Page 19 Source Listing AMSUA 2012-11-20 14:00 tranamsua.f 991 C ---------------------------------------------------------------- 992 993 rinc(4) = 0.00355+(i-1)*.202516+.0845! # of seconds to 994 ! add to scan start 995 ! time for this spot 996 call w3movdat(rinc,idat,jdat) ! update date time 997 998 C Calculate updated day-of-year & second-of-day for use by satazimuth 999 C ------------------------------------------------------------------- 1000 1001 call w3doxdat(jdat,jdow,jdoy,jday) ! calc. updated DOY 1002 sctime = jdat(5)*3600 +jdat(6)*60 +jdat(7) +jdat(8)/1000. 1003 sctime1 = sctime 1004 1005 C Estimate satellite azimuth angle and solar azimuth angle 1006 C -------------------------------------------------------- 1007 1008 call satazimuth(jdat(1),jdoy,sctime,slat(i),slon(i), 1009 x rlocaz(i),aazimuth(i),sazimuth(i)) 1010 ccccc write(stdout,*)i,sctime,slat(i),slon(i),aazimuth(i), 1011 cccccx sazimuth(i),soza(i) 1012 1013 C Store output information for this spot 1014 C -------------------------------------- 1015 1016 bdata(1) = jsat ! satellite id 1017 bdata(2) = jtype ! data type indicator 1018 bdata(3) = jdat(1) ! 4-digit year 1019 bdata(4) = jdat(2) ! month 1020 bdata(5) = jdat(3) ! day 1021 bdata(6) = jdat(5) ! hour 1022 bdata(7) = jdat(6) ! minute 1023 bdata(8) = jdat(7) + jdat(8)/1000. ! second 1024 bdata(9) = lndsea(i) ! land/sea tag 1025 bdata(10)= i ! F-O-V (spot) number 1026 bdata(11)= slat(i) ! latitude 1027 bdata(12)= slon(i) ! longitude 1028 ccccc bdata(13)= rlocaz(i) 1029 bdata(13)= saza(i) ! sat. zenith angle 1030 bdata(14)= soza(i) ! solar zenith angle 1031 bdata(15)= sfchgt(i) ! surface height 1032 bdata(16)= sathgt ! satellite height 1033 bdata(17)= sazimuth(i) ! solar azimuth angle 1034 bdata(18)= aazimuth(i) ! sat. azimuth angle 1035 adata(1:18) = bdata(1:18) 1036 1037 if(process_Tb.eq.'YES') then 1038 if (ikeepb(i).eq.1) then 1039 if (lndsea(i).lt.0.5) nseab = nseab + 1 1040 if (lndsea(i).gt.0.5) nlandb = nlandb + 1 1041 do j = 1,mch 1042 bdata(17+j*nperchan) = tb(j,i) ! brightness temp 1043 bdata(18+j*nperchan) = dt(j) ! cold space 1044 ! temp corr 1045 end do 1046 nrecb = nrecb + 1 1047 call bufr1b(lubfrb,'NC021023',nreal,mch,bdata,nrepb) Page 20 Source Listing AMSUA 2012-11-20 14:00 tranamsua.f 1048 ccccc write(99,*) (bdata(j),j=1,5),sctime1, 1049 ccccc+ (bdata(j),j=9,nreal-2), 1050 ccccc+ (bdata(j),j=nreal+1,ntot,nperchan) 1051 endif 1052 endif 1053 1054 if(process_Ta.eq.'YES') then 1055 if (ikeepa(i).eq.1) then 1056 if (lndsea(i).lt.0.5) nseaa = nseaa + 1 1057 if (lndsea(i).gt.0.5) nlanda = nlanda + 1 1058 do j = 1,mch 1059 adata(17+j*nperchan) = ta(j,i) ! brightness temp 1060 adata(18+j*nperchan) = dt(j) ! cold space 1061 ! temp corr 1062 end do 1063 nreca = nreca + 1 1064 call bufr1b(lubfra,'NC021123',nreal,mch,adata,nrepa) 1065 endif 1066 endif 1067 1068 endif 1069 end do 1070 1071 C Every NPRINT scan lines, print mpos-th record 1072 C --------------------------------------------- 1073 1074 if (mod(nlo,nprint).eq.0) then 1075 if(process_Tb.eq.'YES') then 1076 write(stdout,*)' ' 1077 write(stdout,*)' Tb data for line,rec=',nlo,nrecb 1078 write(stdout,*) (bdata(i),i=1,ntot) 1079 write(stdout,*)' ' 1080 endif 1081 if(process_Ta.eq.'YES') then 1082 write(stdout,*)' ' 1083 write(stdout,*)' Ta data for line,rec=',nlo,nreca 1084 write(stdout,*) (adata(i),i=1,ntot) 1085 write(stdout,*)' ' 1086 endif 1087 endif 1088 1089 C********************************************************************** 1090 C DONE WITH THIS SCAN LINE, READ NEXT SCAN LINE 1091 C********************************************************************** 1092 1093 goto 1200 1094 1095 1600 continue 1096 1097 C All scan lines have been read and processed, summarize 1098 C ------------------------------------------------------ 1099 1100 write(stdout,*)' ' 1101 write(stdout,*)'Done reading raw 1b file' 1102 write(stdout,*)' ' 1103 write(stdout,*)'AMSU-A INGEST STATS:' 1104 write(stdout,*)' no. scan lines = ',nlo,nrecs,nscan Page 21 Source Listing AMSUA 2012-11-20 14:00 tranamsua.f 1105 write(stdout,*)' no. fail good qc = ',nqcbad 1106 write(stdout,*)' no. fail time qc = ',nqctim 1107 write(stdout,*)' no. fail cali qc = ',nqccal 1108 write(stdout,*)' no. fail loc qc = ',nqcloc 1109 write(stdout,*)' no. frame errors = ',nermin,nermaj 1110 write(stdout,*)' no. scans bad qc = ',nskipq 1111 write(stdout,*)' no. scans bad calibr. = ',nskipc 1112 write(stdout,*)' no. scans bad mode = ',nskipm 1113 write(stdout,*)' no. bad (lat,lon) = ',nbadl 1114 write(stdout,*)' no. zero (lat,lon) = ',nopos 1115 write(stdout,*)' no. bad radiances = ',nbadr 1116 if(process_Tb.eq.'YES') then 1117 write(stdout,*)' no. bad calibration = ',nbadc 1118 write(stdout,*)' no. scans with bad Tb = ',nskiptb 1119 write(stdout,*)' no. bad Tb values = ',nbadtb 1120 write(stdout,*)' no. land/sea Tb obs = ',nlandb,nseab 1121 write(stdout,*)' no. Tb recs written = ',nrecb 1122 write(stdout,*)' no. Tb BUFR rpts written = ',nrepb 1123 endif 1124 if(process_Ta.eq.'YES') then 1125 write(stdout,*)' no. scans with bad Ta = ',nskipta 1126 write(stdout,*)' no. bad Ta values = ',nbadta 1127 write(stdout,*)' no. land/sea Ta obs = ',nlanda,nseaa 1128 write(stdout,*)' no. Ta recs written = ',nreca 1129 write(stdout,*)' no. Ta BUFR rpts written = ',nrepa 1130 endif 1131 1132 write(stdout,*)' ' 1133 write(stdout,*)'bad radiance,temperature counts per channel' 1134 write(stdout,1020) 1135 1020 format(t1,'channel',t10,'bad rad',t20,'bad Tb',t30,'bad Ta') 1136 sumr = 0. 1137 sumtb = 0. 1138 sumta = 0. 1139 do j = 1,mch 1140 write(stdout,1030) j,badr(j),badtb(j),badta(j) 1141 1030 format(t1,i2,t10,f8.1,t20,f8.1,t30,f8.1) 1142 sumr = sumr + badr(j) 1143 sumtb = sumtb + badtb(j) 1144 sumta = sumta + badta(j) 1145 end do 1146 write(stdout,*)'nbadr,nbadtb,nbadta=',sumr,sumtb,sumta 1147 1148 write(stdout,*)' ' 1149 write(stdout,*)' AMSU-A 1B DECODE COMPLETED' 1150 write(stdout,*)' ' 1151 1152 C Close UNITs 1153 C ----------- 1154 1155 close(lunin) 1156 call closbf(lubfrb) 1157 call closbf(lubfra) 1158 1159 call system('echo YES > Tb') 1160 if(process_Tb.eq.'YES') then 1161 if(nrecb.eq.0) then Page 22 Source Listing AMSUA 2012-11-20 14:00 tranamsua.f 1162 write(stdout,1003) 1163 1003 format(/' NO Tb RECORDS WRITTEN -- DISABLING ALL ', 1164 1 'SUBSEQUENT Tb PROCESSING.'/) 1165 call system('echo NO > Tb') 1166 else 1167 call mesgbc(lubfrb,msgt,icomp) 1168 if(icomp.eq.1) then 1169 print'(/"OUTPUT Tb BUFR FILE MESSAGES '// 1170 . 'C O M P R E S S E D"/"FIRST MESSAGE TYPE FOUND IS",I5/)', 1171 . msgt 1172 elseif(icomp.eq.0) then 1173 print'(/"OUTPUT Tb BUFR FILE MESSAGES '// 1174 . 'U N C O M P R E S S E D"/"FIRST MESSAGE TYPE FOUND IS",'// 1175 . 'I5/)', msgt 1176 elseif(icomp.eq.-1) then 1177 print'(//"ERROR READING OUTPUT Tb BUFR FILE - MESSAGE '// 1178 . 'COMPRESSION UNKNOWN"/)' 1179 elseif(icomp.eq.-3) then 1180 print'(/"OUTPUT Tb BUFR FILE DOES NOT EXIST"/)' 1181 elseif(icomp.eq.-2) then 1182 print'(/"OUTPUT Tb BUFR FILE HAS NO DATA MESSAGES"/'// 1183 . '"FIRST MESSAGE TYPE FOUND IS",I5/)', msgt 1184 endif 1185 endif 1186 endif 1187 1188 call system('echo YES > Ta') 1189 if(process_Ta.eq.'YES') then 1190 if(nreca.eq.0) then 1191 write(stdout,1004) 1192 1004 format(/' NO Ta RECORDS WRITTEN -- DISABLING ALL ', 1193 . 'SUBSEQUENT Ta PROCESSING.'/) 1194 call system('echo NO > Ta') 1195 else 1196 call mesgbc(lubfra,msgt,icomp) 1197 if(icomp.eq.1) then 1198 print'(/"OUTPUT Ta BUFR FILE MESSAGES '// 1199 . 'C O M P R E S S E D"/"FIRST MESSAGE TYPE FOUND IS",I5/)', 1200 . msgt 1201 elseif(icomp.eq.0) then 1202 print'(/"OUTPUT Ta BUFR FILE MESSAGES '// 1203 . 'U N C O M P R E S S E D"/"FIRST MESSAGE TYPE FOUND IS",'// 1204 . 'I5/)', msgt 1205 elseif(icomp.eq.-1) then 1206 print'(//"ERROR READING OUTPUT Ta BUFR FILE - MESSAGE '// 1207 . 'COMPRESSION UNKNOWN"/)' 1208 elseif(icomp.eq.-3) then 1209 print'(/"OUTPUT Ta BUFR FILE DOES NOT EXIST"/)' 1210 elseif(icomp.eq.-2) then 1211 print'(/"OUTPUT Ta BUFR FILE HAS NO DATA MESSAGES"/'// 1212 . '"FIRST MESSAGE TYPE FOUND IS",I5/)', msgt 1213 endif 1214 endif 1215 endif 1216 1217 close(lubfrb) 1218 close(lubfra) Page 23 Source Listing AMSUA 2012-11-20 14:00 tranamsua.f 1219 1220 return 1221 1222 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1223 C ERRORS 1224 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1225 1226 C Error reading Ta to Tb coefficient file 1227 C --------------------------------------- 1228 1229 1800 continue 1230 write(stdout,*)'*** error reading coefficient file ',coefile 1231 close(luncof) 1232 close(lunin) 1233 call closbf(lubfrb) 1234 call closbf(lubfra) 1235 call w3tage('BUFR_TRANAMSUA') 1236 call errexit(2) 1237 1238 C Error reading 1B file header record 1239 C ----------------------------------- 1240 1241 1900 continue 1242 write(stdout,*)' *** error reading hdr record of file ',rawamsu 1243 close(luncof) 1244 close(lunin) 1245 call closbf(lubfrb) 1246 call closbf(lubfra) 1247 call w3tage('BUFR_TRANAMSUA') 1248 call errexit(3) 1249 1250 end Page 24 Source Listing AMSUA 2012-11-20 14:00 Entry Points tranamsua.f ENTRY POINTS Name amsua_ SYMBOL CROSS REFERENCE Name Object Declared Type Bytes Dimen Elements Attributes References 100 Label 423 405,408 1000 Label 700 698 1003 Label 1147 1146 1004 Label 1176 1175 1005 Label 720 719 101 Label 426 416 1010 Label 833 832 102 Label 429 422 1020 Label 1119 1118 1030 Label 1125 1124 1200 Label 635 703,835,1077 1600 Label 1079 646 1800 Label 1213 405,408 1900 Label 1225 459 AAZIMUTH Local 288 R(8) 8 1 30 993,1018 ABS Func 786 scalar 786,787,791,792 ADATA Local 294 R(4) 4 1 48 1019,1043,1044,1048,1068 AMSUA Subr 178 B Local 284 R(8) 8 scalar 895,897 BADR Local 296 R(8) 8 1 15 362,854,870,1124,1126 BADTA Local 297 R(8) 8 1 15 364,923,1124,1128 BADTB Local 296 R(8) 8 1 15 363,914,1124,1127 BDATA Local 294 R(4) 4 1 48 1000,1001,1002,1003,1004,1005,1006 ,1007,1008,1009,1010,1011,1013,101 4,1015,1016,1017,1018,1019,1026,10 27,1031,1062 BUFR1B Subr 1031 1031,1048 C Local 284 R(8) 8 scalar 896,897 C0 Local 293 R(8) 8 1 15 736,851,867 C1 Local 293 R(8) 8 1 15 734,851,867 C2 Local 293 R(8) 8 1 15 732,851,867 CFRQ0 Local 291 R(8) 8 1 15 329 CLIGHT Local 284 R(8) 8 scalar 328 CLOSBF Subr 1140 1140,1141,1217,1218,1229,1230 CNST1 Local 290 R(8) 8 1 15 545,549,895 CNST2 Local 290 R(8) 8 1 15 546,549,896 COEFILE Dummy 178 CHAR 500 scalar ARG,INOUT 392,1214 COMPRESS Scalar 274 CHAR 8 scalar COM COUNTS Local 286 R(8) 8 scalar 850,851,866,867 CWAVE Local 290 R(8) 8 1 15 544,549,889,891 DATTIM Subr 670 670 DT Local 296 R(8) 8 1 15 598,602,900,1027,1044 DTCOLD Local 597 R(4) 4 scalar 597,598 EPS Param 257 R(8) 8 scalar 791,792,890,892 Page 25 Source Listing AMSUA 2012-11-20 14:00 Symbol Table tranamsua.f Name Object Declared Type Bytes Dimen Elements Attributes References ERREXIT Subr 511 511,523,1220,1232 ETA Local 296 R(8) 8 1 15 335,418,420,898,899 F0 Local 295 R(8) 8 2 450 399,405,408,433,439,446,898,900 F1 Local 295 R(8) 8 2 450 400,405,408,434,440,447,898,899 F2 Local 295 R(8) 8 2 450 401,405,408,435,441,448,898,900 I Local 375 I(4) 4 scalar 375,398,399,400,401,404,405,407,40 8,432,433,434,435,437,439,440,441, 444,446,447,448,459,464,465,570,57 1,587,588,595,596,598,601,602,646, 648,649,764,765,768,769,770,777,77 8,779,781,782,783,784,785,786,787, 788,789,791,792,793,794,797,798,81 1,813,814,816,817,819,820,842,847, 857,863,873,884,888,897,898,899,90 0,901,903,904,911,912,915,921,924, 933,939,940,951,961,970,972,977,99 2,993,1008,1009,1010,1011,1013,101 4,1015,1017,1018,1022,1023,1024,10 26,1039,1040,1041,1043,1062,1068 IBADC Local 282 I(4) 4 1 15 711,714,715,719,931 IBADTA Local 886 I(4) 4 scalar 886,940,959 IBADTB Local 885 I(4) 4 scalar 885,939,949 ICHAN Local 281 I(4) 4 1 15 345,375 ICOMP Local 1151 I(4) 4 scalar 1151,1152,1156,1160,1163,1165,1180 ,1181,1185,1189,1192,1194 IDAT Local 280 I(4) 4 1 8 672,673,674,675,676,980 IDDD Local 660 I(4) 4 scalar 660,662 IDT Local 280 I(4) 4 1 2 662,663,670 IERMAJ Local 687 I(4) 4 scalar 687,689,695,699 IERMIN Local 686 I(4) 4 scalar 686,689,694,699 IKEEPA Local 281 I(4) 4 1 30 785,789,794,798,961,972,1039 IKEEPB Local 281 I(4) 4 1 30 784,788,793,797,951,972,1022 ILINE Local 658 I(4) 4 scalar 658,698,719,832 ILS Local 811 I(4) 4 scalar 811,812,815 IMODE1 Local 828 I(4) 4 scalar 828,830,832 IMODE2 Local 829 I(4) 4 scalar 829,830,832 IMX Local 282 I(4) 4 1 3 310 INDAT Local 273 CHAR 4 1 640 465,472,517,531,532,555,556,557,59 7,649,658,659,660,661,682,683,684, 685,686,687,713,758,768,769,770,82 8,829,850,866 IQCBAD Local 682 I(4) 4 scalar 682,688,690,698 IQCCAL Local 684 I(4) 4 scalar 684,688,692,699,713,714 IQCLOC Local 685 I(4) 4 scalar 685,688,693,699 IQCTIM Local 683 I(4) 4 scalar 683,688,691,698 ISUM Local 681 I(4) 4 scalar 681,688,696,709,715,717 ITIME Local 279 I(8) 8 scalar 661,664 IYEAR Local 659 I(4) 4 scalar 659,663 J Local 419 I(4) 4 scalar 419,420,438,439,440,441,445,446,44 7,448,539,540,544,545,546,710,711, 712,714,715,719,730,731,732,734,73 6,848,849,851,854,857,864,865,867, 870,873,887,888,889,891,895,896,89 7,898,899,900,901,903,904,911,912, 914,915,921,923,924,931,933,939,94 Page 26 Source Listing AMSUA 2012-11-20 14:00 Symbol Table tranamsua.f Name Object Declared Type Bytes Dimen Elements Attributes References 0,1025,1026,1027,1042,1043,1044,11 23,1124,1126,1127,1128 JB Local 540 I(4) 4 scalar 540,541,596,597,712,713,757,758,84 7,849,863,865 JB0 Local 541 I(4) 4 scalar 541,542,544,735,736,765,766,768,77 8,779,780,849,850,865,866 JB1 Local 542 I(4) 4 scalar 542,543,545,733,734,735,766,767,76 9,780,781 JB2 Local 543 I(4) 4 scalar 543,546,731,732,733,767,770 JBUF Local 273 CHAR 4 1 640 465,649 JDAT Local 280 I(4) 4 1 8 980,985,986,992,1002,1003,1004,100 5,1006,1007 JDAY Local 985 I(4) 4 scalar 985 JDOW Local 985 I(4) 4 scalar 985 JDOY Local 985 I(4) 4 scalar 985,992 JMX Local 282 I(4) 4 1 3 311 JSAT Local 472 I(4) 4 scalar 472,473,474,475,477,478,479,480,48 2,483,484,485,487,488,489,490,492, 493,494,495,497,498,499,500,502,50 3,504,505,507,509,525,1000 JSAT0 Local 474 I(4) 4 scalar 474,477,479,482,484,487,489,492,49 4,497,499,502,504,507 JTYPE Local 517 I(4) 4 scalar 517,519,521,525,1001 KBUF Local 272 CHAR 1 1 2560 459,544,545,546,646,732,734,736,77 9,781 LANSEA Func 811 I(4) 4 scalar 811 LBYTE Func 472 I(4) 4 scalar 472,517,531,532,658,659,660,661,68 2,683,684,685,686,687,713,758,828, 829,850,866 LINE Local 653 I(4) 4 scalar 653 LL Local 811 I(4) 4 scalar 811,817 LNDSEA Local 281 I(4) 4 1 30 782,813,816,819,1008,1023,1024,104 0,1041 LOG Func 893 scalar 893 LUBFRA Local 354 I(4) 4 scalar 354,385,1048,1141,1180,1202,1218,1 230 LUBFRB Local 354 I(4) 4 scalar 354,383,1031,1140,1151,1201,1217,1 229 LUNCOF Dummy 178 I(4) 4 scalar ARG,INOUT 392,405,408,416,418,430,1215,1227 LUNDX Local 353 I(4) 4 scalar 353,383,385 LUNIN Dummy 178 I(4) 4 scalar ARG,INOUT 456,459,646,1139,1216,1228 MAPFILE Local 275 CHAR 40 1 3 312 MBYTE Func 555 I(4) 4 scalar 555,556,557,597,768,769,770 MCH Param 261 I(4) 4 scalar 267,281,282,284,290,291,292,293,29 5,296,297,372,375,397,445,539,587, 595,601,602,710,719,730,848,887,94 9,959,1025,1031,1042,1048,1123 MESGBC Subr 1151 1151,1180 MIN Func 972 scalar 972 MOD Func 676 scalar 676,1058 MPOS Param 262 I(4) 4 scalar 281,284,288,289,292,295,398,404,40 7,432,437,444,764,777,842,884,970 MSGT Local 1151 I(4) 4 scalar 1151,1155,1159,1167,1180,1184,1188 ,1196 N Local 397 I(4) 4 scalar 397,399,400,401,405,408,418 Page 27 Source Listing AMSUA 2012-11-20 14:00 Symbol Table tranamsua.f Name Object Declared Type Bytes Dimen Elements Attributes References NBADC Local 615 I(4) 4 scalar 615,932,1101 NBADL Local 616 I(4) 4 scalar 616,790,1097 NBADR Local 619 I(4) 4 scalar 619,853,869,1099 NBADTA Local 618 I(4) 4 scalar 618,922,1110 NBADTB Local 617 I(4) 4 scalar 617,913,1103 NBYTE1 Param 258 I(4) 4 scalar 258,272,456,459,646 NBYTE4 Param 258 I(4) 4 scalar 273,464,648 NDT Local 280 I(4) 4 1 5 670,673,674 NERMAJ Local 614 I(4) 4 scalar 614,695,1093 NERMIN Local 613 I(4) 4 scalar 613,694,1093 NLANDA Local 630 I(4) 4 scalar 630,1041,1111 NLANDB Local 629 I(4) 4 scalar 629,1024,1104 NLO Local 633 I(4) 4 scalar 633,652,653,698,719,832,1058,1061, 1067,1088 NOPOS Local 608 I(4) 4 scalar 608,795,1098 NPERCHAN Param 267 I(4) 4 scalar 267,1026,1027,1043,1044 NPRINT Local 365 I(4) 4 scalar 365,1058 NQCBAD Local 609 I(4) 4 scalar 609,690,1089 NQCCAL Local 611 I(4) 4 scalar 611,692,1091 NQCLOC Local 612 I(4) 4 scalar 612,693,1092 NQCTIM Local 610 I(4) 4 scalar 610,691,1090 NREAL Param 267 I(4) 4 scalar 267,372,1031,1048 NRECA Local 621 I(4) 4 scalar 621,1047,1067,1112,1174 NRECB Local 620 I(4) 4 scalar 620,1030,1061,1105,1145 NRECS Local 531 I(4) 4 scalar 531,533,1088 NREPA Local 623 I(4) 4 scalar 623,1048,1113 NREPB Local 622 I(4) 4 scalar 622,1031,1106 NRI Local 458 I(4) 4 scalar 458,459,641,646 NSCAN Local 532 I(4) 4 scalar 532,533,1088 NSEAA Local 632 I(4) 4 scalar 632,1040,1111 NSEAB Local 631 I(4) 4 scalar 631,1023,1104 NSET Param 259 I(4) 4 scalar 275,282 NSKIPC Local 624 I(4) 4 scalar 624,718,1095 NSKIPM Local 626 I(4) 4 scalar 626,831,832,1096 NSKIPQ Local 625 I(4) 4 scalar 625,697,699,1094 NSKIPTA Local 628 I(4) 4 scalar 628,960,1109 NSKIPTB Local 627 I(4) 4 scalar 627,950,1102 NTOT Param 267 I(4) 4 scalar 294,373,1062,1068 OPENBF Subr 383 383,385 P1 Local 284 R(8) 8 scalar 323,891 P2 Local 284 R(8) 8 scalar 323,889 PROCESS_TA Scalar 274 CHAR 8 scalar COM 1038,1065,1108,1173 PROCESS_TB Scalar 274 CHAR 8 scalar COM 1021,1059,1100,1144 RAD Local 292 R(8) 8 2 450 857,873,888 RAD0 Local 286 R(8) 8 scalar RADS Local 286 R(8) 8 scalar 851,852,855,857,867,868,871,873,88 8,890,891 RAWAMSU Dummy 178 CHAR 500 scalar ARG,INOUT 456,1226 REAL_32 Param 255 I(4) 4 scalar 294 REAL_64 Param 256 I(4) 4 scalar 257,284,285,286,287,288,289,290,29 1,292,293,295,296,297 RFAC Param 15 I(4) 4 scalar 456 RINC Local 294 R(4) 4 1 5 671,977,980 RLOCAZ Local 288 R(8) 8 1 30 770,993 RMISS Local 340 R(4) 4 scalar 340,782,783,819,820,855,871,903,90 Page 28 Source Listing AMSUA 2012-11-20 14:00 Symbol Table tranamsua.f Name Object Declared Type Bytes Dimen Elements Attributes References 4,905,915,924,933 RNORM Local 286 R(8) 8 scalar 898,899 SATAZIMUTH Subr 992 992 SATHGT Local 286 R(8) 8 scalar 758,1016 SAZA Local 288 R(8) 8 1 30 769,1013 SAZIMUTH Local 289 R(8) 8 1 30 993,1017 SCALE Local 285 R(8) 8 scalar 554,555,556,557,594,597,756,758,76 3,768,769,770,776,779,781 SCALE13 Local 285 R(8) 8 scalar 728,734 SCALE19 Local 285 R(8) 8 scalar 727,732 SCALE5 Local 285 R(8) 8 scalar SCALE6 Local 285 R(8) 8 scalar 538,544,545,546 SCALE9 Local 285 R(8) 8 scalar 729,736 SCTIME Local 286 R(8) 8 scalar 664,670,986,987,992 SCTIME1 Local 294 R(4) 4 scalar 987 SELECTED_REAL_KIND Func 255 scalar 255,256 SFCHGT Local 292 R(8) 8 1 30 783,814,817,820,1015 SLAT Local 289 R(8) 8 1 30 779,786,791,811,992,1010 SLON Local 289 R(8) 8 1 30 781,787,792,811,992,1011 SOZA Local 288 R(8) 8 1 30 768,1014 STDOUT Local 278 I(4) 4 scalar 352,356,357,370,371,372,373,374,37 5,376,410,417,424,427,476,481,486, 491,496,501,506,509,520,525,533,54 8,601,602,603,698,719,832,1060,106 1,1062,1063,1066,1067,1068,1069,10 84,1085,1086,1087,1088,1089,1090,1 091,1092,1093,1094,1095,1096,1097, 1098,1099,1101,1102,1103,1104,1105 ,1106,1109,1110,1111,1112,1113,111 6,1117,1118,1124,1130,1132,1133,11 34,1146,1175,1214,1226 SUMR Local 1120 R(4) 4 scalar 1120,1126,1130 SUMTA Local 1122 R(4) 4 scalar 1122,1128,1130 SUMTB Local 1121 R(4) 4 scalar 1121,1127,1130 SWITCHES Common 304 24 SYSTEM Subr 1143 1143,1149,1172,1178 TA Local 284 R(8) 8 2 450 897,899,903,921,924,940,1043 TA0 Local 284 R(8) 8 scalar 894,897 TB Local 292 R(8) 8 2 450 899,901,904,911,912,915,933,939,10 26 TB0 Local 286 R(8) 8 scalar 901,905 TERM1 Local 284 R(8) 8 scalar 889,894 TERM2 Local 284 R(8) 8 scalar 891,892,893 TERM3 Local 284 R(8) 8 scalar 893,894 THI Local 318 R(4) 4 scalar 318,912 TLO Local 318 R(4) 4 scalar 318,911 TREF Local 296 R(8) 8 1 15 562,563,568,569,571,577,578,579,58 0,588,601,899 TSHELF11A Local 287 R(8) 8 scalar 555,568,569,571 TSHELF12 Local 287 R(8) 8 scalar 556,577,578,579,580 TSHELF2 Local 287 R(8) 8 scalar 557,562,563 W3DOXDAT Subr 985 985 W3MOVDAT Subr 980 980 W3TAGE Subr 510 510,522,1219,1231 XFLOAT Func 544 R(4) 4 scalar 544,545,546,732,734,736,779,781 Page 29 Source Listing AMSUA 2012-11-20 14:00 Symbol Table tranamsua.f Name Object Declared Type Bytes Dimen Elements Attributes References XSEC Local 286 R(8) 8 scalar 670,675,676 Page 30 Source Listing AMSUA 2012-11-20 14:00 tranamsua.f 1251 1252 SUBROUTINE CHARS(IWORD,LEN,CWORD) 1253 C$$$ SUBPROGRAM DOCUMENTATION BLOCK 1254 C 1255 C SUBPROGRAM: CHARS 1256 C PRGMMR: BERT KATZ ORG: NP20 DATE: 1997-11-06 1257 C 1258 C ABSTRACT: Turns integer into character string of specified length, 1259 C starting at low-order byte of integer. 1260 C 1261 C PROGRAM HISTORY LOG: 1262 C 1997-11-06 Katz -- Original author 1263 C 1264 C USAGE: CALL CHARS(IWORD,LEN,CWORD) 1265 C INPUT ARGUMENT LIST: 1266 C IWORD - INTEGER argument 1267 C LEN - INTEGER argument holding number of low-order bytes 1268 C of first argument to convert into character 1269 C 1270 C OUTPUT ARGUMENT LIST: 1271 C CWORD - CHARACTER argument 1272 C 1273 C REMARKS: 1274 C 1275 C ATTRIBUTES: 1276 C LANGUAGE: FORTRAN 90 1277 C MACHINE: NCEP WCOSS 1278 C 1279 C$$$ 1280 1281 character*1 cword(len) 1282 integer iword 1283 do ic = len , 1 , -1 1284 ibeg = (len - ic) * 8 1285 ichr = ibits(iword,ibeg,8) 1286 cword(ic) = char(ichr) 1287 enddo 1288 return 1289 end Page 31 Source Listing CHARS 2012-11-20 14:00 Entry Points tranamsua.f ENTRY POINTS Name chars_ SYMBOL CROSS REFERENCE Name Object Declared Type Bytes Dimen Elements Attributes References CHAR Func 1270 scalar 1270 CHARS Subr 1236 CWORD Dummy 1236 CHAR 1 1 0 ARG,INOUT 1270 IBEG Local 1268 I(4) 4 scalar 1268,1269 IBITS Func 1269 scalar 1269 IC Local 1267 I(4) 4 scalar 1267,1268,1270 ICHR Local 1269 I(4) 4 scalar 1269,1270 IWORD Dummy 1236 I(4) 4 scalar ARG,INOUT 1269 LEN Dummy 1236 I(4) 4 scalar ARG,INOUT 1265,1267,1268 Page 32 Source Listing CHARS 2012-11-20 14:00 tranamsua.f 1290 1291 SUBROUTINE DATTIM(IDT,SCTIME,NDT,XSEC) 1292 C$$$ SUBPROGRAM DOCUMENTATION BLOCK 1293 C 1294 C SUBPROGRAM: DATTIM 1295 C PRGMMR: KEYSER ORG: NP22 DATE: 2006-07-20 1296 C 1297 C ABSTRACT: Converts year, day-of-year, and second-of-day into year, 1298 C month-of-year, day-of-month, hour-of-day, minute-of-hour, and 1299 C second-of-minute. 1300 C 1301 C PROGRAM HISTORY LOG: 1302 C 1997-11-06 Katz -- Original author 1303 C 2006-07-20 Keyser -- Modified to input real double precision 1304 C second-of-day and output real double precision second-of-minute 1305 C (both had been integer) 1306 C 1307 C USAGE: CALL DATTIM(IDT,SCTIME,NDT,XSEC) 1308 C INPUT ARGUMENT LIST: 1309 C IDT - INTEGER array argument containing two members: 1310 C day-of-year and year. 1311 C SCTIME - REAL*8 argument containing second-of-day. 1312 C 1313 C OUTPUT ARGUMENT LIST: 1314 C NDT - INTEGER array argument containing five members: 1315 C year, month-of-year, day-of-month, hour-of-day, 1316 C and minute-of-hour. 1317 C XSEC - REAL*8 argument containing second-of-minute. 1318 C 1319 C REMARKS: 1320 C NONE 1321 C 1322 C ATTRIBUTES: 1323 C LANGUAGE: FORTRAN 90 1324 C MACHINE: NCEP WCOSS 1325 C 1326 C$$$ 1327 1328 integer,parameter::real_64=selected_real_kind(15,307) 1329 1330 integer idt(2),ndt(5) 1331 integer iday,ihr,imin,imon,isec,iyr,jday,jsec 1332 real(real_64) sctime,xsec 1333 1334 external w3fs26 1335 1336 intrinsic mod 1337 1338 JULIAN(IYR,IDYR) = -31739 + 1461 * (IYR + 4799) / 4 1339 & -3 * ((IYR + 4899) / 100) / 4 + IDYR 1340 1341 jday = idt(1) 1342 iyr = idt(2) 1343 1344 C If year is two digits, convert to 4 digits 1345 C ------------------------------------------ 1346 Page 33 Source Listing DATTIM 2012-11-20 14:00 tranamsua.f 1347 if (iyr.ge.0.and.iyr.le.99) then 1348 if (iyr.lt.21) then 1349 kyr = iyr + 2000 1350 else 1351 kyr = iyr + 1900 1352 endif 1353 else 1354 kyr = iyr 1355 endif 1356 1357 C Compute julian day number as number days after 4713 bc 1358 C ------------------------------------------------------ 1359 1360 idy = jday 1361 jdn = julian(kyr,idy) 1362 1363 call w3fs26(jdn,iyear,jmo,jda,idaywk,idayyr) 1364 1365 imon = jmo 1366 iday = jda 1367 jsec = sctime 1368 ihr = sctime/3600 1369 imin = mod(sctime,3600._8)/60 1370 xsec = sctime - 3600*ihr - 60*imin 1371 ndt(1) = iyr 1372 ndt(2) = imon 1373 ndt(3) = iday 1374 ndt(4) = ihr 1375 ndt(5) = imin 1376 return 1377 end Page 34 Source Listing DATTIM 2012-11-20 14:00 Entry Points tranamsua.f ENTRY POINTS Name dattim_ SYMBOL CROSS REFERENCE Name Object Declared Type Bytes Dimen Elements Attributes References DATTIM Subr 1275 IDAY Local 1315 I(4) 4 scalar 1350,1357 IDAYWK Local 1347 I(4) 4 scalar 1347 IDAYYR Local 1347 I(4) 4 scalar 1347 IDT Dummy 1275 I(4) 4 1 2 ARG,INOUT 1325,1326 IDY Local 1344 I(4) 4 scalar 1344,1345 IHR Local 1315 I(4) 4 scalar 1352,1354,1358 IMIN Local 1315 I(4) 4 scalar 1353,1354,1359 IMON Local 1315 I(4) 4 scalar 1349,1356 ISEC Local 1315 I(4) 4 scalar IYEAR Local 1347 I(4) 4 scalar 1347 IYR Local 1315 I(4) 4 scalar 1326,1331,1332,1333,1335,1338,1355 JDA Local 1347 I(4) 4 scalar 1347,1350 JDAY Local 1315 I(4) 4 scalar 1325,1344 JDN Local 1345 I(4) 4 scalar 1345,1347 JMO Local 1347 I(4) 4 scalar 1347,1349 JSEC Local 1315 I(4) 4 scalar 1351 JULIAN Local 1322 I(4) 4 scalar 1345 KYR Local 1333 I(4) 4 scalar 1333,1335,1338,1345 MOD Func 1320 scalar 1353 NDT Dummy 1275 I(4) 4 1 5 ARG,INOUT 1355,1356,1357,1358,1359 REAL_64 Param 1312 I(4) 4 scalar 1316 SCTIME Dummy 1275 R(8) 8 scalar ARG,INOUT 1351,1352,1353,1354 SELECTED_REAL_KIND Func 1312 scalar 1312 W3FS26 Subr 1318 1347 XSEC Dummy 1275 R(8) 8 scalar ARG,INOUT 1354 Page 35 Source Listing DATTIM 2012-11-20 14:00 tranamsua.f 1378 1379 INTEGER FUNCTION ICHARS(CWORD,LEN) 1380 C$$$ SUBPROGRAM DOCUMENTATION BLOCK 1381 C 1382 C SUBPROGRAM: ICHARS 1383 C PRGMMR: BERT KATZ ORG: NP20 DATE: 1997-11-05 1384 C 1385 C ABSTRACT: Turns character string of specified length into integer. 1386 C 1387 C PROGRAM HISTORY LOG: 1388 C 1997-11-05 Katz -- Original author 1389 C 1390 C USAGE: ICHARS(CWORD,LEN) 1391 C INPUT ARGUMENT LIST: 1392 C CWORD - CHARACTER*1 array argument 1393 C LEN - INTEGER argument holding length of cword 1394 C 1395 C REMARKS: 1396 C NONE 1397 C 1398 C ATTRIBUTES: 1399 C LANGUAGE: FORTRAN 90 1400 C MACHINE: NCEP WCOSS 1401 C 1402 C$$$ 1403 1404 character*1 cword(len) 1405 lchars = 0 1406 do ic = len , 1 , -1 1407 ibeg = (len - ic) * 8 1408 icmove = mova2i(cword(ic)) 1409 call mvbits(icmove,0,8,lchars,ibeg) 1410 enddo 1411 ichars = lchars 1412 return 1413 end Page 36 Source Listing ICHARS 2012-11-20 14:00 Entry Points tranamsua.f ENTRY POINTS Name ichars_ SYMBOL CROSS REFERENCE Name Object Declared Type Bytes Dimen Elements Attributes References CWORD Dummy 1363 CHAR 1 1 0 ARG,INOUT 1392 IBEG Local 1391 I(4) 4 scalar 1391,1393 IC Local 1390 I(4) 4 scalar 1390,1391,1392 ICHARS Func 1363 I(4) 4 scalar 1395 ICHARS@0 Local 1363 I(4) 4 scalar ICMOVE Local 1392 I(4) 4 scalar 1392,1393 LCHARS Local 1389 I(4) 4 scalar 1389,1393,1395 LEN Dummy 1363 I(4) 4 scalar ARG,INOUT 1388,1390,1391 MOVA2I Func 1392 I(4) 4 scalar 1392 MVBITS Intrin 1393 1393 Page 37 Source Listing ICHARS 2012-11-20 14:00 tranamsua.f 1414 cfpp$ expand(ichars,lbit) 1415 1416 INTEGER FUNCTION LANSEA(RLAT,RLON,LEVEL) 1417 C$$$ SUBPROGRAM DOCUMENTATION BLOCK 1418 C 1419 C SUBPROGRAM: LANSEA 1420 C PRGMMR: BERT KATZ ORG: NP20 DATE: 1997-11-05 1421 C 1422 C ABSTRACT: Calculates topography, land/sea status from latitude and 1423 C longitude. 1424 C 1425 C PROGRAM HISTORY LOG: 1426 C 1997-11-05 Katz -- Original author 1427 C 1428 C USAGE: LANSEA(RLAT,RLON,LEVEL) 1429 C INPUT ARGUMENT LIST: 1430 C RLAT - INTEGER argument containing scaled latitude 1431 C RLON - INTEGER argument containing scaled longitude 1432 C 1433 C OUTPUT ARGUMENT LIST: 1434 C LEVEL - INTEGER argument containing scaled topography 1435 C 1436 C INPUT FILES: 1437 C UNIT 41 - Binary low-resolution topography file 1438 C 1439 C REMARKS: 1440 C NONE 1441 C 1442 C ATTRIBUTES: 1443 C LANGUAGE: FORTRAN 90 1444 C MACHINE: NCEP WCOSS 1445 C 1446 C$$$ 1447 1448 include 'rfac.inc' 1449 1466 integer,parameter::real_64=selected_real_kind(15,307) 1467 integer ilat,ilon,level 1468 real(real_64) rlat,rlon 1469 real slat,slon,xlon 1470 integer iopn,iu,lan,last,lat,lenr,lon 1471 character*12 name 1472 character*4 iflag(12),kelev(192) 1473 character*2 ielev(360) 1474 1475 integer lbit,ichars 1476 1477 external lbit,ichars 1478 1479 intrinsic float,max0 1480 1481 equivalence (iflag(1),kelev(1)), (ielev(1),kelev(13)) 1482 1483 data name/'lowtopog.dat'/,iu/41/,lenr/768/,last/0/,iopn/0/ 1484 1485 save iopn,kelev,last 1486 Page 38 Source Listing LANSEA 2012-11-20 14:00 tranamsua.f 1487 if (iopn.eq.0) then 1488 open (iu,recl=lenr/rfac, 1489 & file=name,access='direct',status='old') 1490 iopn = 1 1491 endif 1492 lan = 0 1493 level = 0 1494 slat = rlat 1495 slon = rlon 1496 lat = slat + 1. 1497 if (slat.lt.0.) lat = slat 1498 lat = max0(lat,-87) 1499 lat = 91 - lat 1500 if (lat.eq.last) go to 100 1501 read (iu,rec=lat) kelev 1502 last = lat 1503 100 xlon = slon 1504 if (xlon.lt.0.) xlon = xlon + 360. 1505 lon = xlon 1506 if (lon.eq.360) lon = 0 1507 lon = lon + 1 1508 lan = lbit(lon,iflag) 1509 if (lan.ne.0) then 1510 ltemp = ichars(ielev(lon),2) 1511 if (btest(ltemp,15)) then 1512 ltemp = ior(ltemp,-65536) 1513 endif 1514 level = ltemp 1515 endif 1516 lansea = 2 - lan 1517 return 1518 end Page 39 Source Listing LANSEA 2012-11-20 14:00 Entry Points tranamsua.f ENTRY POINTS Name lansea_ SYMBOL CROSS REFERENCE Name Object Declared Type Bytes Dimen Elements Attributes References 100 Label 1471 1468 BTEST Func 1479 scalar 1479 FLOAT Func 1447 scalar ICHARS Func 1443 I(4) 4 scalar 1478 IELEV Local 1441 CHAR 2 1 360 1478 IFLAG Local 1440 CHAR 4 1 12 1476 ILAT Local 1435 I(4) 4 scalar ILON Local 1435 I(4) 4 scalar IOPN Local 1438 I(4) 4 scalar 1451,1455,1458 IOR Func 1480 scalar 1480 IU Local 1438 I(4) 4 scalar 1451,1456,1469 KELEV Local 1440 CHAR 4 1 192 1469 LAN Local 1438 I(4) 4 scalar 1460,1476,1477,1484 LANSEA Func 1400 I(4) 4 scalar 1484 LANSEA@0 Local 1400 I(4) 4 scalar LAST Local 1438 I(4) 4 scalar 1451,1468,1470 LAT Local 1438 I(4) 4 scalar 1464,1465,1466,1467,1468,1469,1470 LBIT Func 1443 I(4) 4 scalar 1476 LENR Local 1438 I(4) 4 scalar 1451,1456 LEVEL Dummy 1400 I(4) 4 scalar ARG,INOUT 1461,1482 LON Local 1438 I(4) 4 scalar 1473,1474,1475,1476,1478 LTEMP Local 1478 I(4) 4 scalar 1478,1479,1480,1482 MAX0 Func 1447 scalar 1466 NAME Local 1439 CHAR 12 scalar 1451,1457 REAL_64 Param 1434 I(4) 4 scalar 1436 RFAC Param 15 I(4) 4 scalar 1456 RLAT Dummy 1400 R(8) 8 scalar ARG,INOUT 1462 RLON Dummy 1400 R(8) 8 scalar ARG,INOUT 1463 SELECTED_REAL_KIND Func 1434 scalar 1434 SLAT Local 1437 R(4) 4 scalar 1462,1464,1465 SLON Local 1437 R(4) 4 scalar 1463,1471 XLON Local 1437 R(4) 4 scalar 1471,1472,1473 Page 40 Source Listing LANSEA 2012-11-20 14:00 tranamsua.f 1519 cfpp$ expand(ichars) 1520 1521 INTEGER FUNCTION LBIT(J,ARRAY) 1522 C$$$ SUBPROGRAM DOCUMENTATION BLOCK 1523 C 1524 C SUBPROGRAM: LBIT 1525 C PRGMMR: BERT KATZ ORG: NP20 DATE: 1997-11-05 1526 C 1527 C ABSTRACT: Extracts j'th bit from array of CHARACTER*4. 1528 C 1529 C PROGRAM HISTORY LOG: 1530 C 1997-11-05 Katz -- Original author 1531 C 1532 C USAGE: LBIT(J,ARRAY) 1533 C INPUT ARGUMENT LIST: 1534 C J - INTEGER argument 1535 C ARRAY - CHARACTER*4 array argument 1536 C 1537 C REMARKS: 1538 C NONE 1539 C 1540 C ATTRIBUTES: 1541 C LANGUAGE: FORTRAN 90 1542 C MACHINE: NCEP WCOSS 1543 C 1544 C$$$ 1545 1546 integer j 1547 character*4 array(*) 1548 integer ibit,jout,jw,jword,nbit 1549 1550 integer ichars 1551 external ichars 1552 1553 intrinsic btest 1554 1555 jw = (j-1)/32 1556 nbit = j - jw*32 1557 jword = ichars(array(jw+1),4) 1558 ibit = 32 - nbit 1559 jout = 0 1560 if (btest(jword,ibit)) jout = 1 1561 lbit = jout 1562 return 1563 end Page 41 Source Listing LBIT 2012-11-20 14:00 Entry Points tranamsua.f ENTRY POINTS Name lbit_ SYMBOL CROSS REFERENCE Name Object Declared Type Bytes Dimen Elements Attributes References ARRAY Dummy 1489 CHAR 4 1 0 ARG,INOUT 1525 BTEST Func 1521 scalar 1528 IBIT Local 1516 I(4) 4 scalar 1526,1528 ICHARS Func 1518 I(4) 4 scalar 1525 J Dummy 1489 I(4) 4 scalar ARG,INOUT 1523,1524 JOUT Local 1516 I(4) 4 scalar 1527,1528,1529 JW Local 1516 I(4) 4 scalar 1523,1524,1525 JWORD Local 1516 I(4) 4 scalar 1525,1528 LBIT Func 1489 I(4) 4 scalar 1529 LBIT@0 Local 1489 I(4) 4 scalar NBIT Local 1516 I(4) 4 scalar 1524,1526 Page 42 Source Listing LBIT 2012-11-20 14:00 tranamsua.f 1564 cfpp$ expand(ichars) 1565 1566 INTEGER FUNCTION MBYTE(J,LENGTH,JARRAY) 1567 C$$$ SUBPROGRAM DOCUMENTATION BLOCK 1568 C 1569 C SUBPROGRAM: MBYTE 1570 C PRGMMR: BERT KATZ ORG: NP20 DATE: 1997-11-05 1571 C 1572 C ABSTRACT: Extracts bit string from array of CHARACTER*4 and 1573 C converts it to INTEGER. Entry point MBYTE propagates sign bit 1574 C in result; entry point LBYTE does not. 1575 C 1576 C PROGRAM HISTORY LOG: 1577 C 1997-11-05 Katz -- Original author 1578 C 1579 C USAGE: MBYTE(J,LENGTH,JARRAY) 1580 C INPUT ARGUMENT LIST: 1581 C J - INTEGER argument containing starting bit 1582 C LENGTH - integer argument containing number of bits 1583 C (maximum value 32) 1584 C JARRAY - CHARACTER*4 array argument 1585 C 1586 C OUTPUT FILES: 1587 C UNIT 06 - printout 1588 C 1589 C REMARKS: 1590 C NONE 1591 C 1592 C ATTRIBUTES: 1593 C LANGUAGE: FORTRAN 90 1594 C MACHINE: NCEP WCOSS 1595 C 1596 C$$$ 1597 1598 integer j,length 1599 character*4 jarray(*) 1600 integer inleft,jbit,kompl,mflag,n,nlj,nrj, 1601 + nsh,nword 1602 integer(8) item,jleft,jrite,mask 1603 integer(8) kounts(33) 1604 1605 integer ichars 1606 external ichars 1607 1608 intrinsic iand,ior,mod 1609 1610 integer lbyte 1611 1612 data kounts/1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192, 1613 + 16384,32768,65536,131072,262144,524288,1048576,2097152, 1614 + 4194304,8388608,16777216,33554432,67108864,134217728, 1615 + 268435456,536870912,1073741824,2147483648_8,0/ 1616 1617 C ENTRY MBYTE 1618 C ----------- 1619 1620 mflag = 1 Page 43 Source Listing MBYTE 2012-11-20 14:00 tranamsua.f 1621 110 nword = (j-1)/32 + 1 1622 if (length.lt.1 .or. length.gt.32) write (*,fmt=120) length 1623 120 format (' improper byte length in mbyte or lbyte',i10) 1624 nlj = mod(j-1,32) 1625 nrj = 32 - length - nlj 1626 if (nrj.lt.0) go to 150 1627 item = ichars(jarray(nword),4) 1628 kompl = 33 - length - nlj 1629 mask = -kounts(kompl) 1630 item = iand(item,mask) 1631 item = item/kounts(nrj+1) 1632 mask = kounts(length+1) - 1 1633 item = iand(item,mask) 1634 130 if (mflag.eq.0) go to 140 1635 1636 c ... means logical byte 1637 1638 mbyte = item 1639 jbit = iand(kounts(length),item) 1640 if (jbit.eq.0) return 1641 1642 c ... need sign extension 1643 1644 mask = -mask - 1 1645 item = ior(item,mask) 1646 mbyte = item 1647 return 1648 1649 C ENTRY LBYTE 1650 C ----------- 1651 1652 entry lbyte(j,length,jarray) 1653 1654 mflag = 0 1655 go to 110 1656 1657 140 lbyte = item 1658 return 1659 1660 c ... byte spans two words 1661 1662 150 inleft = length + nrj 1663 mask = kounts(inleft+1) - 1 1664 jleft = ichars(jarray(nword),4) 1665 jleft = iand(jleft,mask) 1666 nsh = 1 - nrj 1667 jleft = jleft*kounts(nsh) 1668 n = 1 - nrj 1669 jrite = ichars(jarray(nword+1),4) 1670 kompl = 33 + nrj 1671 mask = -kounts(kompl) 1672 jrite = iand(jrite,mask) 1673 jrite = jrite/kounts(nrj+33) 1674 mask = kounts(n) - 1 1675 jrite = iand(jrite,mask) 1676 item = ior(jrite,jleft) 1677 mask = kounts(length+1) - 1 Page 44 Source Listing MBYTE 2012-11-20 14:00 tranamsua.f 1678 go to 130 1679 end ENTRY POINTS Name Name lbyte_ mbyte_ SYMBOL CROSS REFERENCE Name Object Declared Type Bytes Dimen Elements Attributes References 110 Label 1589 1623 120 Label 1591 1590 130 Label 1602 1646 140 Label 1625 1602 150 Label 1630 1594 IAND Func 1576 scalar 1598,1601,1607,1633,1640,1643 ICHARS Func 1573 I(4) 4 scalar 1595,1632,1637 INLEFT Local 1568 I(4) 4 scalar 1630,1631 IOR Func 1576 scalar 1613,1644 ITEM Local 1570 I(8) 8 scalar 1595,1598,1599,1601,1606,1607,1613 ,1614,1625,1644 J Dummy 1534 I(4) 4 scalar ARG,INOUT 1589,1592 JARRAY Dummy 1534 CHAR 4 1 0 ARG,INOUT 1595,1632,1637 JBIT Local 1568 I(4) 4 scalar 1607,1608 JLEFT Local 1570 I(8) 8 scalar 1632,1633,1635,1644 JRITE Local 1570 I(8) 8 scalar 1637,1640,1641,1643,1644 KOMPL Local 1568 I(4) 4 scalar 1596,1597,1638,1639 KOUNTS Local 1571 I(8) 8 1 33 1580,1597,1599,1600,1607,1631,1635 ,1639,1641,1642,1645 LBYTE Func 1578 I(4) 4 scalar 1625 LENGTH Dummy 1534 I(4) 4 scalar ARG,INOUT 1590,1593,1596,1600,1607,1630,1645 MASK Local 1570 I(8) 8 scalar 1597,1598,1600,1601,1612,1613,1631 ,1633,1639,1640,1642,1643,1645 MBYTE Func 1534 I(4) 4 scalar 1606,1614 MBYTE@0 Local 1534 I(4) 4 scalar MFLAG Local 1568 I(4) 4 scalar 1588,1602,1622 MOD Func 1576 scalar 1592 N Local 1568 I(4) 4 scalar 1636,1642 NLJ Local 1568 I(4) 4 scalar 1592,1593,1596 NRJ Local 1568 I(4) 4 scalar 1593,1594,1599,1630,1634,1636,1638 ,1641 NSH Local 1569 I(4) 4 scalar 1634,1635 NWORD Local 1569 I(4) 4 scalar 1589,1595,1632,1637 Page 45 Source Listing MBYTE 2012-11-20 14:00 tranamsua.f 1680 cfpp$ expand(ichars) 1681 1682 REAL FUNCTION XFLOAT(JB,IARRAY) 1683 C$$$ SUBPROGRAM DOCUMENTATION BLOCK 1684 C 1685 C SUBPROGRAM: XFLOAT 1686 C PRGMMR: BERT KATZ ORG: NP20 DATE: 1997-11-05 1687 C 1688 C ABSTRACT: Takes two consecutive elements of CHARACTER*2 array and 1689 C forms a floating point number from them. 1690 C 1691 C PROGRAM HISTORY LOG: 1692 C 1997-11-05 Katz -- Original author 1693 C 1694 C USAGE: XFLOAT(JB,IARRAY) 1695 C INPUT ARGUMENT LIST: 1696 C JB - INTEGER argument containing array location 1697 C IARRAY - CHARACTER*2 array argument 1698 C 1699 C REMARKS: 1700 C NONE 1701 C 1702 C ATTRIBUTES: 1703 C LANGUAGE: FORTRAN 90 1704 C MACHINE: NCEP WCOSS 1705 C 1706 C$$$ 1707 1708 integer jb 1709 integer(8) mask 1710 data mask/x'ffffffff00000000'/ 1711 ccccc character*2 iarray(2) 1712 character*2 iarray(2000) 1713 character*4 conv 1714 real xf 1715 integer j 1716 integer(8) jj 1717 integer in(2) 1718 1719 integer ichars 1720 external ichars 1721 1722 intrinsic btest,ior 1723 1724 j = jb 1725 conv(1:2) = iarray(j) 1726 conv(3:4) = iarray(j+1) 1727 jj = ichars(conv,4) 1728 if (btest(jj,31_8)) then 1729 jj = ior(jj,mask) 1730 endif 1731 xf = jj 1732 xfloat = xf 1733 return 1734 end Page 46 Source Listing XFLOAT 2012-11-20 14:00 Entry Points tranamsua.f ENTRY POINTS Name xfloat_ SYMBOL CROSS REFERENCE Name Object Declared Type Bytes Dimen Elements Attributes References BTEST Func 1690 scalar 1696 CONV Local 1681 CHAR 4 scalar 1693,1694,1695 IARRAY Dummy 1650 CHAR 2 1 2000 ARG,INOUT 1693,1694 ICHARS Func 1687 I(4) 4 scalar 1695 IN Local 1685 I(4) 4 1 2 IOR Func 1690 scalar 1697 J Local 1683 I(4) 4 scalar 1692,1693,1694 JB Dummy 1650 I(4) 4 scalar ARG,INOUT 1692 JJ Local 1684 I(8) 8 scalar 1695,1696,1697,1699 MASK Local 1677 I(8) 8 scalar 1678,1697 XF Local 1682 R(4) 4 scalar 1699,1700 XFLOAT Func 1650 R(4) 4 scalar 1700 XFLOAT@0 Local 1650 R(4) 4 scalar Page 47 Source Listing XFLOAT 2012-11-20 14:00 tranamsua.f 1735 1736 SUBROUTINE SATAZIMUTH(IYEAR,IDAY,STIME,SLAT,SLON,RAZIMUTH, 1737 + AAZIMUTH,SAZIMUTH) 1738 C$$$ SUBPROGRAM DOCUMENTATION BLOCK 1739 C 1740 C SUBPROGRAM: SATAZIMUTH 1741 C PRGMMR: DERBER ORG: NP2 DATE: 2006-04-06 1742 C 1743 C ABSTRACT: Calculates satellite azimuth angle and solar azimuth angle 1744 C given the relative azimuth angle and the date/time. 1745 C 1746 C PROGRAM HISTORY LOG: 1747 C 2006-04-06 Derber - Original author 1748 C 1749 C USAGE: CALL SATAZIMUTH(IYEAR,IDAY,STIME,RAZIMUTH,SLAT,SLON, 1750 C AAZIMUTH,SAZIMUTH) 1751 C INPUT ARGUMENT LIST: 1752 C IYEAR - Year 1753 C IDAY - Day of year 1754 C STIME - Time of day (seconds) 1755 C SLAT - Latitude (degrees N) 1756 C SLON - Longitude (degrees, + E, - W) 1757 C RAZIMUTH - Relative azimuth angle (degrees true) 1758 C 1759 C OUTPUT ARGUMENT LIST: 1760 C AAZIMUTH - Satellite azimuth angle (degrees true, 0-360) 1761 C SAZIMUTH - Solar azimuth angle (degrees true, 0-360) 1762 C 1763 C REMARKS: 1764 C NONE 1765 C 1766 C ATTRIBUTES: 1767 C LANGUAGE: FORTRAN 90 1768 C MACHINE: NCEP WCOSS 1769 C 1770 C$$$ 1771 1772 integer,parameter::real_64=selected_real_kind(15,307) 1773 real(real_64) slon,slat,tcen,omega,gmeananomsun,sec,y 1774 real(real_64) razimuth,aazimuth,deg2rad,rad2deg 1775 real(real_64) sinm,sin2m,sin3m,suneqofcenter,geommeanlongsun 1776 real(real_64) suntruelon,sunapparentlong,meanobliquity 1777 real(real_64) obliquitycorrection,decl,ayr 1778 real(real_64) byr,sin210,sin410,cos210,eccearthorbit 1779 real(real_64) eqtime,tst,sazimuth,sozanew,stime 1780 integer iday 1781 1782 external w3fs26 1783 1784 ajday(byr,jday) = floor(365.25*(byr+4716.)) - 1785 + floor(byr/100.)+ 1786 + floor(float(floor(byr/100.))/4.)-2452639.5+float(jday) 1787 deg2rad = acos(-1.)/180. 1788 rad2deg = 1./deg2rad 1789 1790 ayr = float(iyear-1) 1791 Page 48 Source Listing SATAZIMUTH 2012-11-20 14:00 tranamsua.f 1792 C Calculate relative time to 2000.0 (in centuries) 1793 C ------------------------------------------------ 1794 1795 tcen = (ajday(ayr,iday)+stime/86400.)/36525.0 1796 omega=(125.04-1934.136*tcen)*deg2rad 1797 gmeananomsun=(357.52911+tcen*(35999.05029- 1798 + 0.0001537*tcen))*deg2rad 1799 sinm=sin(gmeananomsun) 1800 sin2m=sin(2.0*gmeananomsun) 1801 sin3m=sin(3.0*gmeananomsun) 1802 suneqofcenter=(sinm*(1.914602- 1803 + tcen*(0.004817+0.000014*tcen))+sin2m*(0.019993- 1804 + 0.000101*tcen)+sin3m*0.000289)*deg2rad 1805 geommeanlongsun=280.46646+tcen*(36000.76983+0.0003032*tcen) 1806 if(geommeanlongsun > 360.0) geommeanlongsun=geommeanlongsun-360. 1807 if(geommeanlongsun < 0.0) geommeanlongsun=geommeanlongsun+360. 1808 geommeanlongsun=deg2rad*geommeanlongsun 1809 suntruelon=geommeanlongsun+suneqofcenter 1810 sunapparentlong=suntruelon-(0.00569-0.00478*sin(omega))*deg2rad 1811 sec=21.448-tcen*(46.8150+tcen*(0.00059-tcen*0.001813)) 1812 meanobliquity=23.0+(26.0+(sec/60.))/60. 1813 obliquitycorrection= (meanobliquity + 0.00256*cos(omega))*deg2rad 1814 decl = asin(sin(obliquitycorrection)*sin(sunapparentlong)) 1815 ccccc solarrightascension = atan2(cos(obliquitycorrection)* 1816 ccccc+ sin(sunapparentlong),cos(sunapparentlong)) 1817 y = tan(obliquitycorrection/2.0) 1818 y = y*y 1819 sin210=sin(2.0*geommeanlongsun) 1820 sin410=sin(4.0*geommeanlongsun) 1821 cos210=cos(2.0*geommeanlongsun) 1822 eccearthorbit=0.016708634-tcen*(0.000042037+ 1823 + 0.0000001287*tcen) 1824 eqtime=(y*sin210 1825 + -2.0*eccearthorbit*sinm 1826 + +4.0*eccearthorbit*y*sinm*cos210 1827 + -0.5*y*y*sin410 1828 + -1.25*eccearthorbit*eccearthorbit*sin2m)*4.0*rad2deg 1829 1830 tst=eqtime+4.*slon+stime/60. 1831 ha=(tst/4.0)-180. 1832 if(ha > 180.)ha=ha-360. 1833 if(ha < -180.)ha=360.+ha 1834 sozanew=acos(sin(deg2rad*slat)*sin(decl)+cos(deg2rad*slat)* 1835 + cos(decl)*cos(deg2rad*ha)) 1836 sazimuth = acos((sin(decl)-sin(deg2rad*slat)*cos(sozanew))/ 1837 + (cos(deg2rad*slat)*sin(sozanew)))*rad2deg 1838 if(ha > 0.)sazimuth=-sazimuth 1839 ccccc if(sazimuth > 180.)sazimuth=sazimuth-360. 1840 ccccc if(sazimuth < -180.)sazimuth=360.+sazimuth 1841 1842 aazimuth= razimuth+sazimuth 1843 if(aazimuth > 180.)aazimuth=aazimuth-360. 1844 if(aazimuth < -180.)aazimuth=360.+aazimuth 1845 1846 C Correct aziumuth angles to be 0-360 degrees true (for BUFR) 1847 C ----------------------------------------------------------- 1848 Page 49 Source Listing SATAZIMUTH 2012-11-20 14:00 tranamsua.f 1849 if(aazimuth < 0.)aazimuth=360.+aazimuth 1850 if(sazimuth < 0.)sazimuth=360.+sazimuth 1851 1852 return 1853 end ENTRY POINTS Name satazimuth_ SYMBOL CROSS REFERENCE Name Object Declared Type Bytes Dimen Elements Attributes References AAZIMUTH Dummy 1705 R(8) 8 scalar ARG,INOUT 1810,1811,1812,1817 ACOS Func 1755 scalar 1755,1802,1804 AJDAY Local 1752 R(4) 4 scalar 1763 ASIN Func 1782 scalar 1782 AYR Local 1745 R(8) 8 scalar 1758,1763 BYR Local 1746 R(8) 8 scalar COS Func 1781 scalar 1781,1789,1802,1803,1804,1805 COS210 Local 1746 R(8) 8 scalar 1789,1794 DECL Local 1745 R(8) 8 scalar 1782,1802,1803,1804 DEG2RAD Local 1742 R(8) 8 scalar 1755,1756,1764,1766,1772,1776,1778 ,1781,1802,1803,1804,1805 ECCEARTHORBIT Local 1746 R(8) 8 scalar 1790,1793,1794,1796 EQTIME Local 1747 R(8) 8 scalar 1792,1798 FLOAT Func 1754 scalar 1754,1758 FLOOR Func 1752 scalar 1752,1753,1754 GEOMMEANLONGSUN Local 1743 R(8) 8 scalar 1773,1774,1775,1776,1777,1787,1788 ,1789 GMEANANOMSUN Local 1741 R(8) 8 scalar 1765,1767,1768,1769 HA Local 1799 R(4) 4 scalar 1799,1800,1801,1803,1806 IDAY Dummy 1704 I(4) 4 scalar ARG,INOUT 1763 IYEAR Dummy 1704 I(4) 4 scalar ARG,INOUT 1758 MEANOBLIQUITY Local 1744 R(8) 8 scalar 1780,1781 OBLIQUITYCORRECTION Local 1745 R(8) 8 scalar 1781,1782,1785 OMEGA Local 1741 R(8) 8 scalar 1764,1778,1781 RAD2DEG Local 1742 R(8) 8 scalar 1756,1796,1805 RAZIMUTH Dummy 1704 R(8) 8 scalar ARG,INOUT 1810 REAL_64 Param 1740 I(4) 4 scalar 1741,1742,1743,1744,1745,1746,1747 SATAZIMUTH Subr 1704 SAZIMUTH Dummy 1705 R(8) 8 scalar ARG,INOUT 1804,1806,1810,1818 SEC Local 1741 R(8) 8 scalar 1779,1780 SELECTED_REAL_KIND Func 1740 scalar 1740 SIN Func 1767 scalar 1767,1768,1769,1778,1782,1787,1788 ,1802,1804,1805 SIN210 Local 1746 R(8) 8 scalar 1787,1792 SIN2M Local 1743 R(8) 8 scalar 1768,1771,1796 SIN3M Local 1743 R(8) 8 scalar 1769,1772 SIN410 Local 1746 R(8) 8 scalar 1788,1795 SINM Local 1743 R(8) 8 scalar 1767,1770,1793,1794 SLAT Dummy 1704 R(8) 8 scalar ARG,INOUT 1802,1804,1805 Page 50 Source Listing SATAZIMUTH 2012-11-20 14:00 Symbol Table tranamsua.f Name Object Declared Type Bytes Dimen Elements Attributes References SLON Dummy 1704 R(8) 8 scalar ARG,INOUT 1798 SOZANEW Local 1747 R(8) 8 scalar 1802,1804,1805 STIME Dummy 1704 R(8) 8 scalar ARG,INOUT 1763,1798 SUNAPPARENTLONG Local 1744 R(8) 8 scalar 1778,1782 SUNEQOFCENTER Local 1743 R(8) 8 scalar 1770,1777 SUNTRUELON Local 1744 R(8) 8 scalar 1777,1778 TAN Func 1785 scalar 1785 TCEN Local 1741 R(8) 8 scalar 1763,1764,1765,1766,1771,1772,1773 ,1779,1790,1791 TST Local 1747 R(8) 8 scalar 1798,1799 W3FS26 Subr 1750 scalar Y Local 1741 R(8) 8 scalar 1785,1786,1792,1794,1795 Page 51 Source Listing SATAZIMUTH 2012-11-20 14:00 Subprograms/Common Blocks tranamsua.f SUBPROGRAMS/COMMON BLOCKS Name Object Declared Type Bytes Dimen Elements Attributes References AMSUA Subr 178 BUFR_TRANAMSUA Prog 1 CHARS Subr 1236 DATTIM Subr 1275 ICHARS Func 1363 I(4) 4 scalar 1395 LANSEA Func 1400 I(4) 4 scalar 1484 LBIT Func 1489 I(4) 4 scalar 1529 MBYTE Func 1534 I(4) 4 scalar 1606,1614 SATAZIMUTH Subr 1704 SWITCHES Common 147 24 SWITCHES Common 304 24 XFLOAT Func 1650 R(4) 4 scalar 1700 COMPILER OPTIONS BEING USED -align nocommons -align nodcommons -align noqcommons -align records -align nosequence -align norec1byte -align norec2byte -align norec4byte -align norec8byte -align norec16byte -altparam -assume accuracy_sensitive -assume nobscc -assume nobuffered_io -assume byterecl -assume nocc_omp -assume nocstring -assume nodummy_aliases -assume nofpe_summary -assume noieee_fpe_flags -assume nominus0 -assume noold_boz -assume old_unit_star -assume noold_ldout_format -assume noold_logical_ldio -assume old_maxminloc -assume old_xor -assume protect_constants -assume noprotect_parens -assume split_common -assume source_include -assume nostd_intent_in -assume nostd_mod_proc_name -assume norealloc_lhs -assume underscore -assume no2underscores no -auto -auto_scalar no -bintext -ccdefault default -check noargs -check noarg_temp_created -check nobounds -check noformat -check nooutput_conversion -check nooverflow -check nopointers -check power -check noshape -check nounderflow -check nouninitialized -coarray-num-procs 0 no -coarray-config-file -convert big_endian -cross_reference -D __INTEL_COMPILER=1210 -D __unix__ -D __unix -D __linux__ -D __linux -D __gnu_linux__ -D unix -D linux -D __ELF__ -D __x86_64 -D __x86_64__ -D _MT -D __INTEL_COMPILER_BUILD_DATE=20120612 -D __i686 -D __i686__ -D __pentiumpro -D __pentiumpro__ Page 52 Source Listing SATAZIMUTH 2012-11-20 14:00 tranamsua.f -D __pentium4 -D __pentium4__ -D __tune_pentium4__ -D __SSE2__ -D __SSE__ -D __MMX__ -double_size 64 no -d_lines no -Qdyncom -error_limit 30 no -f66 no -f77rtl no -fast -fpscomp nofilesfromcmd -fpscomp nogeneral -fpscomp noioformat -fpscomp noldio_spacing -fpscomp nologicals -fixed no -fpconstant -fpe3 -fprm nearest no -ftz -fp_model noprecise -fp_model nofast -fp_model strict -fp_model nosource -fp_model nodouble -fp_model noextended -fp_model novery_fast -fp_model noexcept -fp_model nono_except -fp_modbits nofp_contract -fp_modbits nono_fp_contract -fp_modbits nofenv_access -fp_modbits nono_fenv_access -fp_modbits nocx_limited_range -fp_modbits nono_cx_limited_range -fp_modbits noprec_div -fp_modbits nono_prec_div -fp_modbits noprec_sqrt -fp_modbits nono_prec_sqrt -fp_modbits noftz -fp_modbits no_ftz -fp_modbits nointrin_limited_range -fp_modbits nono_intrin_limited_range -fp_modbits notrunc_compares -fp_modbits nono_trunc_compares -fp_modbits noieee_nan_compares -fp_modbits nono_ieee_nan_compares -fp_modbits nohonor_f32_conversion -fp_modbits nono_honor_f32_conversion -fp_modbits nohonor_f64_conversion -fp_modbits nono_honor_f64_conversion -fp_modbits nono_x87_copy -fp_modbits nono_no_x87_copy -fp_modbits noexception_semantics -fp_modbits nono_exception_semantics -fp_modbits noprecise_libm_functions -fp_modbits nono_precise_libm_functions -heap_arrays 0 no -threadprivate_compat -g0 -iface nomixed_str_len_arg -iface nono_mixed_str_len_arg no -intconstant -integer_size 32 no -mixed_str_len_arg no -module -names lowercase no -noinclude -O2 no -pad_source -real_size 32 no -recursive -reentrancy none no -sharable_localsaves -vec=simd -show noinclude -show map -show options no -syntax_only no -threadcom no -U no -vms -w noall -w nonone -w alignments -w noargument_checking -w nodeclarations -w general -w noignore_bounds -w noignore_loc -w nointerfaces -w notruncated_source -w uncalled -w uninitialized -w nounused -w usage -includepath : /gpfs/tp2/usrx/local/intel/composer_xe_2011_sp1.11.339/compiler/include/,.f,./.f,/usrx/local/intel/composerxe/mkl/include/.f, /usrx/local/intel/composerxe/tbb/include/.f,/gpfs/tp2/usrx/local/intel/composer_xe_2011_sp1.11.339/compiler/include/intel64/.f, /gpfs/tp2/usrx/local/intel/composer_xe_2011_sp1.11.339/compiler/include/.f,/usr/local/include/.f,/usr/lib/gcc/x86_64-redhat-linux/4.4.6/include/.f, /usr/include/.f,/usr/include/.f -list filename : tranamsua.lst -o filename : none Page 53 Source Listing SATAZIMUTH 2012-11-20 14:00 tranamsua.f COMPILER: Intel(R) Fortran 12.1-2100