C Copyright(c) 2004, Space Science and Engineering Center, UW-Madison C Refer to "McIDAS Software Acquisition and Distribution Policies" C in the file mcidas/data/license.txt C *** $Id: imgconv.f,v 1.10 2019/01/15 19:12:47 daves Exp $ *** SUBROUTINE MAIN0 C ? IMGCONV - Image Conversion C ? IMGCONV sdataset1..sdatasetn ddataset "comment C ? Parameters: C ? sdataset1..n | source ADDE dataset name(s) and position(s); specify C ? as alias.position or group/descriptor.position (no def C ? for alias or group/descriptor, def=0 for position) C ? ddataset | destination ADDE dataset name and absolute position; C ? specify as alias.position or group/descriptor.position; C ? only positive integers are valid for the position number C ? "comment | description placed in the image file directory's memo field; C ? 32 characters maximum (def=no comment) C ? Keywords: C ? ********************** Image Positioning Keywords ********************** C ? LATLON=lat lon | latitude and longitude of sdataset image(s) to place C ? at ddataset image location specified in PLACE keyword C ? LINELE=line ele sys | line and element numbers and coordinate system C ? (F or I for file and image, respectively) of the C ? sdataset image to place at the ddataset image C ? location specified in the PLACE keyword C ? (line def=0, ele def=line, sys def=F) C ? PLACE= | location in the ddataset image to place the point specified in C ? the LATLON, STATION or LINELE keyword; use one of the following: C ? ULEFT upper-left corner of image (def when LINELE keyword used) C ? CENTER center of image (def when LATLON or STATION keyword used) C ? STATION= | station ID, e.g., MSN, YSSY, 72201 to place at the ddataset C ? image location specified in PLACE keyword C ? *********************** Image Selection Keywords *********************** C ? DAY= | selects the most recent sdataset image from the specified day C ? (def=most recent image) C ? RTIME=bmin emin | selects the most recent sdataset image in the C ? specified range of minutes in the current hour; C ? overrides the TIME keyword (def=most recent image) C ? TIME=btime etime | selects the most recent sdataset image in the C ? specified time range (def=most recent image) C ? *********************** Image Sampling Keywords *********************** C ? BAND=b1..bn | applies the function to the specified image bands C ? (def=existing band if image is single-band; no def if C ? image is multi-band) C ? MAG=lmag emag | line and element magnification of the data (def=1 lmag) C ? SIZE=line ele | applies the function to the specified number of sdataset C ? lines and elements (def=480 640) C ? =ALL | applies the function to the entire sdataset image; all Image C ? Positioning keywords are invalid with SIZE=ALL C ? TRACKing=YES | tracks image if allowed by the server; tracking allows C ? commands to receive data as it's being ingested by the C ? server; commands complete when their data request has C ? been fulfilled; tracking is only supported on some C ? servers (e.g., Area, GVAR and MTSAT) and it can be C ? disabled by the server administrator C ? =NO | does not track image so server only sends data that is C ? currently available C ? (def=value of MCTRACK environment variable; if MCTRACK is C ? not set, def=value set by server or server administrator) C ? UNIT=| calibration types for the data (def=BRIT) C ? **************************** Input Keywords **************************** C ? LLMT=l1..ln | lower limit thresholds for each input data value in the C ? equation C ? ULMT=u1..un | upper limit thresholds for each input data value in the C ? equation C ? SPOS=longitude | Nominal longitude of satellite subpoint, which is C ? used by Limb Darkening option C ? (default from image navigation, if available) C ? IDT=limb darkening coefficient (1-99) C ? default is satellite and band dependent C ? *********************** Conversion Format Keywords ******************** C ? CONv= NOR | GOES normalize visible brightness for changes in sun angle C ? ANOR | Other satellite normalize visible brightness C ? = BNOR | Blue channel normalization with Rayleigh scattering correction C ? HN | convert IR counts to height in feet/250 for navigated images C ? LD | corrects for limb darkening C ? ASCI | converts height images to ASCII file HTASCII. C ? OFF=o1..oN | offsets for each term in the equation (def=0.0) C ? SAMpl=s1 | sampling size (def=20.0) C ? **************************** Output Keywords *************************** C ? PROD= | product name for destination image, up to 4 chars (def=PROD) C ? SCALE=prodlo prodhi britlo brithi | minimum and maximum product and C ? brightness values for the calibration (def=0 255 0 255) C ? ZERO=DATA | assigns zero pixel values as data and includes them in the C ? analysis (def=MISS, meaning treat as missing data) C ? **************************** Remarks *********************************** C ? This is a program that is used to generate the brightness C ? normalization functions along with other experimental conversions. C ? The NOR option is used to do the brightness normalization on eight C ? bit visible images that have a square root digization, such as the C ? GOES data. The ANOR option is used for linear digization data C ? normalization that end up as a square root digization brightness C ? normalized image that can be merged with GOES visible normalized data. C ? The original brightness normalization would square the brightness C ? count, divide by the cosine of the solar zenith angle, and then square C ? root the image again. Previously this algorithm over brightened the C ? clouds near the eastern horizon in the morning and near the western C ? horizon in the evening because of forward scattering by the clouds. C ? Recently, a correction factor was put in that computes a C ? bi-directional correction for the divide by the cosine of the solar C ? zenith angle. The bi-directional correction was measured using C ? low level cumulus clouds in the tropics. C ? The ANOR is similar to NOR, C ? but it doesn't first square the brightness value, but C ? rather multiplies it by 255. This ANOR is used for the Meteosat-5 C ? and -7 images. The new MSG Met-8 double byte data appears to be put C ? into a square root digization for the single byte data so NOR is C ? used for Met-8 as well as for the GAC and HRPT NOAA data. C ? The BNOR option is a blue channel brightness normalization. It is C ? the same normalization as NOR except it also removes the blue C ? light scattered by the atmosphere. C ? The LD limb darkening correction for GOES-16 through GOES-19 is now C ? channel dependent. The LD correction for all other satellites C ? is only for the 6.7 micron water vapor channel. C ? The limb darkening coeficient idt depends on the path length through C ? the atmosphere and the absorbing gas in the band of interest. C ? The specific idt depends on if one wants correction for the surface C ? or if one wants correction for the clouds. The default idt C ? for water vapor bands is for cloud level. The default idt for C ? IR bands is for ground level. If one wants correction for levels C ? other than the default, the key word IDT= has been added for C ? user supplied idt coeficients. C ? C ? This is site specific code that is maintained in McIDAS -XRD as the C ? IMGNORM command. C ? ---------- implicit none include 'areaparm.inc' c --- constants integer MAXBND ! maximum band number parameter ( MAXBND = 32 ) integer MAXELE ! maximum size of line buffer parameter ( MAXELE = 22000 ) integer MAXOPN ! maximum number of open images parameter ( MAXOPN = 100 ) integer MAXSRT ! maximum number of sort clauses parameter ( MAXSRT = 20 ) integer NUMKEY ! number of keywords parameter ( NUMKEY = 23 ) integer DIRSIZ ! file directory size parameter ( DIRSIZ = 64 ) integer XDIRSIZ ! expanded file directory size parameter ( XDIRSIZ = 65 ) integer MAXCAL ! calibration array size parameter ( MAXCAL = 128 ) integer MAXCARD ! number of comment cards parameter ( MAXCARD = 500 ) integer MXCDSZ ! maximum navigation size parameter (MXCDSZ = 5*128) c --- external functions character*12 cfj character*12 cfu,cfg character*4 clit integer idnint integer ifc integer iftok integer index integer isdgch integer lbo integer len_trim integer lwtrunc integer lit integer mcacou integer mcadir integer mcadrd integer mcaget integer mcalin integer mcanav integer mcaout integer mcaput integer mcasort integer mccmd integer mccmddbl integer mccmdint integer mccmdkey integer mccmdnum integer mccmdquo integer mccmdstr integer nkwp integer nvprep integer nv1ini integer nv1sae integer nv1opt real cos real acos double precision dkwp c --- local variables character* 4 csat ! satellite type flag character* 4 def_unit ! default units character* 4 def_zero ! default zero character* 4 prd_name ! product name character*8 cday ! day character*8 ctime ! time character*6 chgt ! cloud height character*8 cnum ! number of points character*8 clat ! latitude *1000 character*9 clon ! longitude *1000 character*35 cout ! output line of data character*12 coutname ! output file name character*12 cal ! calibration units character*12 carea ! destination dataset character*12 cband ! character form of band number character*12 cconv ! conversion option character*12 cdarea ! character form of destination character*12 chms ! time character*12 cname ! dataset name character*12 cpos ! dataset position character*12 cunit(MAXOPN) ! units of data values character*12 cunitd ! default unit character*12 cval1 ! generic output char variable character*12 cval2 ! generic output char varible character*12 cyd ! date character*12 czero(MAXOPN) ! zeros of data values character*12 czr ! zero pixel values from cmd line character*12 keywords(NUMKEY) ! array of keywords character*12 size_val ! character value of SIZE keyword character*40 c1 ! source dataset name (without group) character*40 c2 ! temp variable character*40 dest ! destination dataset name character*40 def_name ! default destination dataset name character*40 dir_sorts(3) ! directory sort list character*40 name ! source dataset name character*40 sorts(MAXSRT) ! sort clauses for server character*160 cmd ! command line character*160 ctext ! memo text character*240 card ! comment card double precision dblo ! minimum scale brightness value double precision dbo(MAXOPN) ! offset array double precision ddbud ! upper limit default double precision ddbld ! lower limit default double precision ddbod ! offset default double precision ddsmp ! sampling size double precision dbhi ! maximum scale brightness double precision dbl(MAXOPN) ! lower limit array double precision dbu(MAXOPN) ! upper limit array double precision dnlo ! minimum product value double precision dnhi ! maximum product value double precision dsmp(MAXOPN) ! sampling size array double precision dval ! floating point element value double precision dvax ! modified element value double precision dsf(MAXOPN) ! scaling factor array double precision dsum(MAXDFELEMENTS) ! summation array double precision scale ! scale of source data double precision xzero ! floating point zero data value integer band_pos ! position of bandin sort list integer calsiz ! size of the calibration array integer cards(MAXCARD) ! source comment cards integer cmd_len ! length of the command line integer cnt ! ir count integer darea ! destination image number integer day ! image date integer def_element_size ! default element size integer def_line_size ! default line size integer dest_pos ! pos of destination image integer dlen ! length of destination name integer dot ! position of dot in name integer elem_status ! default element size status integer element_size ! sample element size integer emag ! element magnification integer handle(MAXOPN) ! array of source/destiantion handles integer hhmmss ! time integer i ! generic loop index integer ia(MAXOPN) ! image array integer ib(MAXOPN) ! band array integer iband ! band number integer icalb(MAXCAL) ! calibration array integer iconv ! conversion flag integer iday ! integer idd ! integer idir1(XDIRSIZ) ! source image file directory integer idir2(XDIRSIZ) ! destination image file directory integer idt ! limb darkening correction coefficient integer idt_status ! idt parameter call status integer iele ! element index integer iflag(MAXDFELEMENTS) ! data error flag array integer ihgt ! height in feet integer ilat ! latitude*1000 integer ilon ! longitude*1000 integer inum ! number flag integer iposit ! character position of "." in name integer ioff ! offset integer iout(MAXDFELEMENTS) ! output array integer irstat ! generic function return status integer isamp ! granularity of processing same sun angle integer iscale ! scaling factor integer iss ! satellite type integer izero(MAXOPN) ! zero data array integer ival ! generic output variable integer j ! generic loop index integer lat ! point latitude integer line ! line index integer line_size ! sample line size integer line_status ! default line size status integer lmag ! line magnification integer look(256) ! look-up table for cloud tops integer lthck ! cloud top thickness in meters integer mag_flag ! blowup magnification flag integer mareas ! number of entered source areas integer mbands ! number of entered source bands integer msorts ! copy of sort count integer name_len ! length of source dataset name integer nareas ! number of images to read/write integer navarr(MXCDSZ) ! navigation block integer navsiz ! size of the navigation block integer nb ! band index integer nband ! band identifer integer nbands ! number of bands integer ncard ! numebr of comment cards integer ncoefs ! number of coefficients integer ne ! element index integer nelems ! number of source elements integer nl ! line index integer nlines ! number of source lines integer nsorts ! number of sort clauses integer num ! generic number integer numt ! number of thunderstorm points integer pbot ! pressure at bottom of cloud integer phase ! integer pmb ! cloud top pressure in millibars integer pos_len ! position length integer satln ! satellite line integer start ! data starting point in a line integer status ! generic function status integer temp ! ir temperature integer time ! image time integer valcod ! validity code integer xcor ! image x-coordinate integer ycor ! image y-coordinate integer xres ! image element resolution integer yres ! image line resolution integer yyddd ! julian day integer z ! cloud top height (km * 10) integer ibd ! array pointer for bi-directional scattering real cosang (MAXDFELEMENTS) ! cosine zenith angle of sun real cosx real dec ! declination of sun real dlin ! image line real dele ! image element real dumb ! dummy variable real gha ! Greenwich hour angle of sun real pi ! real satang ! satellite zenith angle real cossat(MAXDFELEMENTS) ! cos of sat zenith angle real x ! real xlat ! point latitude real xlon ! point longitude real xin ! navigation input array real xout(2) ! naviation output array real sunang ! angle to sun real xsunang real xsatang real relang ! relative angle between sun and sat real xrelang ! absolute value of relang real splat ! subpoint latitude of satellite real splon ! subpoint longitude of satellite real corr(MAXDFELEMENTS) ! correction for sun angle real sang ! solar zenith angle real xrel ! relative angle in radians real bd(181) ! full bidirectional reflectance array real bdang ! bidirectional angle calculation real ssat(MAXDFELEMENTS) ! sine of satellite angle real ssun(MAXDFELEMENTS) ! sine of solar angle real scat ! scattering angle between sun and satellite real zlon ! -xlon common /imgconvcommon/dsum, & iflag, & iout, & cosang, & cossat, & corr CHARACTER*2 CR CHARACTER*2 LF INTEGER*2 ICR INTEGER*2 ILF EQUIVALENCE(ICR,CR) EQUIVALENCE(ILF,LF) common /cosx/ cosx data pi/3.14159265/ data keywords/'LAT.LON', 'LIN.ELE', 'PLA.CE', 'STA.TION', 'BAN.D', & 'PRO.D', 'SIZ.E', 'CON.V', & 'LLM.T', 'MCO.N', 'OFF', 'SCA.LE', 'SAM.PL', & 'ULM.T', 'UNI.T', 'ZER.O', 'DAY', 'TIM.E', 'RTI.ME', & 'MAG','TRACK.ING','SPOS','IDT'/ C CARRAGE RETURN AND LINE FEED CHARACTERS DATA ICR/Z'00d'/ DATA ILF/Z'00a'/ c fill bidirectional array do 15 j=1,100 15 bd(j)=1.0 do 16 j=101,140 16 bd(j)=bd(j-1)+0.03 do 17 j=141,165 17 bd(j)=bd(j-1)+0.06 do 18 j=166,180 18 bd(j)=bd(j-1)+.35 c --- assume fail condition and set return code to 1; return code c --- is set to 0 upon success of the program call mccodeset(1) inum=0 c --- validate the keywords if( mccmdkey(NUMKEY,keywords).ne.0 ) goto 1000 c --- check the number of positional parameters c --- this will indicate the destination file dest_pos = mccmdnum(' ') c --- get the command line status = mccmd( cmd ) cmd_len = len_trim( cmd ) call mcgetdaytime(yyddd,hhmmss) cyd = cfj(yyddd) chms = cfu(hhmmss) ncard = cmd_len / 67 + 1 card(1:80) = cyd(8:12)//' '//chms(1:6)//' '//cmd(1:67) card(81:160) = ' '//cmd(68:134) card(161:199) = ' '//cmd(135:160) c --- get name of destination file if( mccmdstr(' ', dest_pos, ' ', c1).lt.0 ) then call edest('Invalid Destination File ='//c1,0) goto 1000 elseif (c1 .eq. ' ') then call edest('A destination file is required', 0) call edest('Use DEST= to specify a destination file', 0) goto 1000 endif c --- check if name is old area number form dot = index(c1,'.') if( isdgch(c1).eq.1 ) then darea = iftok(c1) c1 = ' ' elseif( dot.eq.0 ) then call edest('destination image POSITION not specified',0) goto 1000 else darea = iftok( c1(dot+1:) ) c2 = c1(1:dot-1) c1=c2 endif dest = c1 carea = cfu(darea) c --- squeeze out the extranious blank characters call bsquez(dest) c --- name of the product (def=PROD) if( mccmdstr('PRO.D', 1, 'PROD', prd_name).lt.0 ) then call edest('Invalid Product Name ='//cname,0) goto 1000 endif CALL DDEST('OPTION IS '//CCONV(1:3),0) c --- conversion option if( mccmdstr('CON.V',1,'NONE',cconv).lt.0 ) goto 1000 call mcupcase( cconv ) if( cconv(1:3).eq.'NOR' ) then iconv = 1 elseif (cconv(1:2).eq.'HN' ) then iconv = 2 elseif (cconv(1:2).eq.'LD' ) then iconv = 3 elseif (cconv(1:4).eq.'ASCI') then iconv = 4 elseif (cconv(1:4).eq.'ANOR') then iconv = 5 elseif (cconv(1:4).eq.'BNOR') then iconv = 6 else call edest('Invalid Conversion Option ='//cconv,0) goto 1000 endif c **************** SOURCE IMAGES ****************** c --- number of source files will be the number of positional c --- parameters-1 (last name is always the destiantion name) nareas = dest_pos-1 mareas = nareas if (nareas .le. 0) then call edest('At least one input image is required', 0) goto 1000 endif c --- check for the number of spectral bands c --- If BAND is specified, then the number of bands could c --- change the number of "open" image files nbands = mccmdnum('BAN.D') mbands = nbands if( nbands.eq.0 ) then nbands = nareas else if( nareas.eq.1 ) then nareas = nbands elseif( nbands.eq.1 ) then nbands = nareas else if( nbands.ne.nareas ) then call edest('Invalid Command --- ',0) call edest('Number of source images must match '// & 'the number of spectral bands',0) goto 1000 endif endif endif c --- make sure the number of source images does not c --- exceed the maximum allowable. the output image c --- counts as one open image files (MAXOPN-1) if (nareas .gt. MAXOPN-1) then call edest('Too many input images specified: ', nareas) call edest('Maximum input images allowed is: ', MAXOPN) goto 1000 endif c --- collect all of the global keywords irstat = mcasort(nsorts, sorts, 1) if (irstat .lt. 0) then call mccodeset(2) goto 1000 endif c --- scan sort list for BAND= band_pos = 0 do 1 i=1,nsorts 1 if( sorts(i)(1:4).eq.'BAND' ) band_pos = i c --- check for MAG keyword call ddest('to mag',0) mag_flag = 0 def_line_size = 480 def_element_size = 640 if( mccmdnum( 'MAG' ).ne.0 ) then if( mccmdint('MAG',1,'Line Magnification',1,1,-1,lmag) & .lt.0 ) goto 1000 if( mccmdint('MAG',2,'Element Magnification',lmag,1,-1,emag) & .lt.0 ) goto 1000 c --- for line blowup: reduce the default line request size if( lmag.gt.1 ) then mag_flag = mag_flag+1 def_line_size = def_line_size / lmag endif c --- for element blowup: reduce the default element request size if( emag.gt.1 ) then mag_flag = mag_flag+1 def_element_size = def_element_size / emag endif endif call ddest('mag_flag=',mag_flag) c --- sample size if( mccmdstr('SIZ.E',1,' ',size_val).lt.0 ) then call edest('Invalid SIZE value='//size_val,0) goto 1000 else if(size_val(1:3).eq.'ALL') then line_size = 99999 element_size = 99999 c --- make sure this does not interfer with other keywords if( mccmdnum( 'LAT.LON' ).ne.0 .or. & mccmdnum( 'LIN.ELE' ).ne.0 .or. & mccmdnum( 'PLA.CE' ).ne.0 .or. & mccmdnum( 'STA.TION').ne.0 ) then call edest('Invalid keyword combination',0) call edest('SIZE=ALL can NOT be used with ',0) call edest('keywords LATLON, LINELE,'// & ' PLACE or STATION',0) goto 1000 endif else line_status=mccmdint('SIZ.E',1,'Sample Line Size', & def_line_size,1,99999,line_size) if( line_status.lt.0 ) goto 1000 elem_status=mccmdint('SIZ.E',2,'Sample Element Size', & def_element_size, 1,99999,element_size) if( elem_status.lt.0 ) goto 1000 endif cval1 = cfu(line_size) cval2 = cfu(element_size) nsorts = nsorts+1 sorts(nsorts) = 'SIZE '//cval1//cval2 call bsquez( sorts(nsorts) ) endif c --- setup default unit def_unit = 'BRIT' c --- initialize the server for each source area, placing c --- the returned handle in the array 'handle' def_name = ' ' do 5 i=1,nareas c --- re-set sort list msorts = nsorts c --- Name and position of data status = mccmdstr(' ', i, def_name, c1) if( status.lt.0 ) then call edest('Invalid Image file name='//c1,0) goto 1000 else def_name = c1 endif c --- check for presence of a dot in the name c --- if no dot in name, position is default of 0 dot = index(c1,'.') if( dot.eq.0 ) then name = c1 iposit = 0 else name = c1(1:dot-1) iposit = iftok(c1(dot+1:)) endif c --- If numeric, will be absolute number on default path if( isdgch(name).eq.1 ) then ia(i) = ifc(name) name = cfu(ia(i)) else msorts = msorts+1 sorts(msorts) = 'POS '//cfu(iposit) endif c --- squeeze out the extranious blank characters call bsquez(name) c --- for blowup we need to check if source image is a multiple c --- of the blowup factor if( mag_flag.ne.0 ) then cpos = cfu(iposit) pos_len = len_trim( cpos ) dir_sorts(1) = 'SUBSET '//cpos(1:pos_len)//' '//cpos dir_sorts(2) = 'AUX NO' dir_sorts(3) = 'BAND ALL' if( mcadir(name, 3, dir_sorts, 1).lt.0 ) goto 1000 status = mcadrd(idir1, cards) c --- check if blowup is possible if( lmag.gt.1 .and. MOD(idir1(13),lmag).ne.0 ) then name_len = len_trim( name ) call edest('Line Resolution of '//name(1:name_len)//'.'// & cpos(1:pos_len)//' not divisible by ',lmag) call edest('Line Resolution =',idir1(13)) goto 1000 endif if( emag.gt.1 .and. MOD(idir1(14),emag).ne.0 ) then name_len = len_trim( name ) call edest('Element Resolution of '//name(1:name_len)// & '.'//cpos(1:pos_len)//' not divisible by ',emag) call edest('Element Resolution =',idir1(14)) goto 1000 endif endif c --- calibration type if( mccmdstr('UNI.T', i, def_unit, cal).lt.0 ) then call edest('Invalid Calibration Type ='//cal,0) goto 1000 endif call mcupcase( cal ) def_unit = cal(1:4) c --- set band if( band_pos.ne.0 ) then ib(i) = 0 if( mareas.gt.1 .and. mbands.eq.1 ) then if(mccmdint('BAN.D',i,'Spectral Band',ib(1),1,MAXBND,ib(i)) & .lt.0 ) goto 1000 else if(mccmdint('BAN.D',i,'Spectral Band',ib(i),1,MAXBND,ib(i)) & .lt.0 ) goto 1000 endif cband = cfu( ib(i) ) sorts(band_pos) = 'BAND '//cband(1:2) endif do 3 j=1,msorts 3 call ddest(sorts(j)(1:40)//' index=',j) c --- get the image directory from the server call ddest('name = '//name,0) if( mcaget(name,msorts,sorts,cal,'I4',MAXDFELEMENTS*4,1,idir1, & handle(i) ).ne.0 ) then call edest('Unable to initialize data read for ' // NAME, 0) goto 1000 endif c --- save the scaling factor dsf(i) = idir1(59) call ddest('scale factor=',NINT(dsf(i))) if(idir1(59) .gt. 20000) then call edest('Unable to process VIIRS DNB radiances', 0) go to 1000 endif c c satellite type used to determine limb darkening coefficient iss=idir1(3) nband=idir1(19) num=1 iband=0 do 27 j=1,16 if(nband.eq.num) iband=j 27 num=num*2 idt_status=mccmdint('IDT',1,'LIMB DARKENING',0,0,99,idt) if(idt.ne.0) go to 28 c default coefficient idt=30 c for GOES-12 if(iss.eq.78)then idt=11 c for GOES-13&14 elseif(iss.eq.180.or.iss.eq.182)then idt=11 c for GOES 9 thru 11 elseif(iss.eq.74.or.iss.eq.76.or.iss.eq.72)then idt=20 c c GMS CORRECTTION elseif((iss.eq.12) .or. & (iss.eq.13) .or. & (iss.eq.82) .or. & (iss.eq.83) .or. & (iss.eq.84) .or. & (iss.eq.85)) then idt=18 c c METEOSAT CORRECTION elseif(iss.ge.51.and.iss.le.58)then idt=16 c c GOES-16 TO 19 COEFFICIENTS for each ir band elseif(iss.ge.186.and.iss.le.192) then if (iband.eq.1) idt=0 if (iband.eq.2) idt=0 if (iband.eq.3) idt=0 if (iband.eq.4) idt=0 if (iband.eq.5) idt=0 if (iband.eq.6) idt=0 if (iband.eq.7) idt=20 if (iband.eq.8) idt=35 if (iband.eq.9) idt=60 if (iband.eq.10)idt=70 if (iband.eq.11)idt=30 if (iband.eq.12)idt=60 if (iband.eq.13)idt=20 if (iband.eq.14)idt=20 if (iband.eq.15)idt=20 if (iband.eq.16)idt=40 endif 28 if(iconv.eq.3) call sdest('limb darkening correction=',idt) 5 continue c **************** SOURCE IMAGES ****************** c **************** DESTINATION IMAGE ************** c --- create the destination files directory entry call movw(DIRSIZ,idir1,idir2) idir2(1) = 0 idir2(11) = 1 c --- allow for blowups if( lmag.gt.1 ) idir2(12) = idir2(12)/lmag if( emag.gt.1 ) idir2(13) = idir2(13)/emag idir2(14) = 1 idir2(15) = 0 c creation date/time statements were removed. c m0makara, makara take care of these. idir2(19) = 1 cc add word 20 as zero for done flag idir2(20) = 0 idir2(33) = darea idir2(36) = 0 idir2(49) = 0 idir2(50) = 0 idir2(51) = 0 idir2(52) = lit('PRD ') idir2(53) = lit('BRIT') idir2(57) = 0 idir2(60) = 0 idir2(61) = 0 idir2(64) = ncard c --- get the "text field if( mccmdquo(ctext).lt.0 ) then call edest('Invalid Text Field',0) goto 1000 else call movcw(ctext(1:32),idir2(25)) endif c --- calculate an offset into the line of data c --- based on the word length of the offset + 1 c --- this is where the data for each line will start start = idir2(15) / 4 + 1 c --- store the valcode from the area directory valcod = idir2(36) c --- move the DATA (34) and CAL (63) offsets to make c --- room for the PROD navigation (35) navsiz = idir2(63) - idir2(35) if( idir2(63).eq.0 ) navsiz = idir2(34)-idir2(35) idir2(63) = idir2(35) + navsiz calsiz = MAXCAL*4 idir2(34) = idir2(63) + calsiz c --- get the navigation block from the 1st source image file if( mcanav(handle(1), navarr).lt.0 ) then call edest('Unable to obtain navigation for ' & // 'destination image',0) goto 1000 endif c --- number of coefficients ncoefs = nareas+1 c --- scale factor dsf(ncoefs) = 1 c --- minimum product value if( mccmddbl('SCA.LE',1,'Scale',0.d0,1.d0,-1.d0,dnlo) & .lt.0 ) goto 1000 dnlo = dnlo * dsf(ncoefs) c --- maximum product value if( mccmddbl('SCA.LE',2,'Scale',255.d0,1.d0,-1.d0,dnhi) & .lt.0 ) goto 1000 dnhi = dnhi * dsf(ncoefs) c --- minimum brightness value if( mccmddbl('SCA.LE',3,'Scale',0.d0,1.d0,-1.d0,dblo) & .lt.0 ) goto 1000 c --- maximum brightness value if( mccmddbl('SCA.LE',4,'Scale',255.d0,1.d0,-1.d0,dbhi) & .lt.0 ) goto 1000 c --- compute the scale of the data scale = (dbhi-dblo) / (dnhi-dnlo) c --- create the products calibration block call zerow(calsiz/4, icalb) icalb(1) = lit( prd_name ) icalb(2) = idnint(dnlo) icalb(3) = idnint(dnhi) icalb(4) = idnint(dblo) icalb(5) = idnint(dbhi) c--- intialize destination file with nav and cal sorts(1) = 'POS ' // cfu(darea) if( mcaput(dest, 1, sorts, idir2, navarr, icalb).ne.0 ) then call edest('Unable to initialize destination file: '//dest, 0) goto 1000 endif c --- set default values for the first group of equation parameters ddsmp = 20.0D0 ddbod = 0.0D0 ddbld = -1.0D+35 ddbud = 1.0D+35 cunitd = 'BRIT' def_zero = 'MISS' c --- specify equation parameters do 10 nb = 1,nareas c --- re-set default values for bands if( nb.gt.1 ) then ddbod = dbo(1) ddbud = dbu(1) cunitd = cunit(1) endif c --- offsets if( mccmddbl('OFF',nb,'Offset value',ddbod,1.d0,-1.d0,dbo(nb)) & .lt.0 ) goto 1000 c --- sampling size if( mccmddbl('SAM.PL',nb,'Sampling value',ddsmp,1.d0,20.d0, & dsmp(nb)).lt.0 ) goto 1000 isamp = idnint( dsmp(nb) ) call ddest('Sampling Size=', isamp) c --- lower limits if( mccmddbl('LLM.T',nb,'Lower Limit',ddbld,1.d0,-1.d0,dbl(nb)) & .lt.0 ) goto 1000 call ddest('Lower Limit=',NINT( dbl(nb) )) ddbld = dbl(nb) c --- upper limits if( mccmddbl('ULM.T',nb,'Upper Limit',ddbud,1.d0,-1.d0,dbu(nb)) & .lt.0 ) goto 1000 call ddest('Upper Limit=',NINT( dbu(nb) )) ddbud = dbu(nb) 10 continue c --- set zero pixel value option do 20 nb = 1,nareas izero(nb) = 0 if( mccmdstr('ZER.O',nb,def_zero,czero(nb)).lt.0 ) then call edest('Invalid Zero Pixel value ='//czero(nb),0) goto 1000 else call mcupcase( czero(nb) ) def_zero = czero(nb)(1:4) czr = czero(nb) if( czr(1:1).lt.'A' .or. czr(1:1).gt.'Z' ) then if( mccmddbl('ZER.O',nb,'Missing Value',0.D0,1.D0, & -1.D0,xzero).lt.0 ) goto 1000 izero(nb) = idnint( xzero * dsf(nb) ) czero(nb) = 'MISS' endif endif 20 continue c *************** READ/WRITE IMAGE FILES **************** c --- get area information day = idir1(4) time = idir1(5) ycor = idir1(6) xcor = idir1(7) yres = idir1(12) xres = idir1(13) c --- read the navigation status = mcanav( handle(1), navarr ) if( status.lt.0 ) then call edest('Failed to read source navigation',0) goto 1000 endif c ------ initialize the source navigation if( nvprep( 1, navarr ).ne.0 ) then call edest('Failed to initialize navigation',0) goto 1000 endif c --- initialize the navigation status = nv1ini(2, 'LL ') c c find subpoint of satellite splat=0.0 splon=0.0 status = nv1opt('SPOS',xin,xout) c --- limb darkening (LD) option requires subpoint if(iconv.eq.3 .and. status.ne.0 .and. nkwp('SPOS').eq.0) then call edest('Satellite subpoint not available',0) call edest('Enter using keyword SPOS',0) return endif splat=xout(1) splon=-xout(2) splon=-dkwp('SPOS', 1, -dble(splon)) c --- delete contents of file HTASCII if(iconv.eq.4) then coutname='HTASCII' status=lwtrunc(coutname) numt=0 endif c --- number of lines and elements to read nlines = idir2(9) nelems = idir2(10) c --- line loop do 100 nl = 1,nlines do 151 ne =1,nelems iflag(ne) = 0 dsum(ne) = 1.0D0 151 continue c --- image loop do 200 nb = 1,nareas c --- read a line of data if( mcalin(handle(nb), iout).ne.0 ) then call edest('Failed to read data at line: ', nl-1) call mccodeset(2) goto 1000 endif ccc if (iconv.eq.2 ) then if(MOD(nl,20).eq.1) then iscale = idnint(dbhi) ioff = idnint(dbo(nb)) satln = ycor + yres * nl + 20/2 c c --- convert the satellite line into approximate latitude c --- using center of line c dlin = FLOAT(satln) dele = xcor + nelems * xres/2 if(nv1sae(dlin,dele,dumb,xlat,xlon,dumb).eq.0) then lat = INT(xlat + .5) phase = 0 idd = MOD(day,1000) if(lat.LT.0) phase = 180 iday = idd + phase if(iday.gt.365) iday = iday - 365 x = iday * 2.0 * PI/365. cosx = COS(x) endif c c --- assume a constant correction for a block of lines. c --- generate a look-up table for the block c do 400 i = 1,256 cnt = i - 1 c c --- convert the ir count to temperature c if(cnt.lt.176) then temp = (660 - cnt) * 10/2 else temp = (418 - cnt) * 10 endif call ttoz(lat,idd,temp,z,pmb,lthck,pbot) c c --- convert to feet and then apply scaling factors c look(i) = z * 328/iscale + ioff if(look(i).lt.0) look(i) = 0 if(look(i).gt.255) look(i) = 255 400 continue endif endif line=ycor+nl*yres c --- element loop do 300 ne = 1,nelems c --- modify the source value if((iconv.eq.1).or.(iconv.eq.5).or.(iconv.eq.6)) then ioff = idnint(dbo(nb))+1 c --- get solar declination and hour angle if (nl.eq.1)call solarp(day,time,gha,dec,xlat,xlon) c --- get solar angle every sample interval ccc if (MOD(ne,isamp).eq.1) then ccc if (MOD(nl,isamp).eq.1) then ccc iele=xcor+(ne+isamp/2)*xres if (MOD(ne,isamp).eq.0) then if (MOD(nl,isamp).eq.0) then iele=xcor+ne*xres dlin=float(line) dele=float(iele) if(nv1sae(dlin,dele,0.0,xlat,xlon,dumb).eq.0)then zlon=-xlon call zsun(day,time,xlat,xlon,gha,dec,cosang(ne)) sang=acos(cosang(ne)) call angless(day,time,xlat,zlon,gha,dec,splat,splon, &satang,sunang,relang) xlon=-xlon xrelang=abs(relang) xsatang=satang*pi/180 ssat(ne)=sin(xsatang) ssun(ne)=sin(sang) ibd=xrelang*ssat(ne)*ssun(ne)+1 if(ibd.gt.180) ibd=180 if(ibd.lt.1) ibd=1 corr(ne)=1.0/(bd(ibd)) cc if(nv1sae(dlin,dele,0.0,xlat,xlon,dumb).eq.0)then cc call zsun(day,time,xlat,xlon,gha,dec,cosang(ne)) cc sang=acos(cosang(ne)) cc corr(ne)=1./(1.+sang/1.57) endif endif if (iconv.eq.6 ) then cossat(ne)=1. iele=xcor+(ne+isamp/2)*xres dlin=float(line) dele=float(iele) if(nv1sae(dlin,dele,0.0,xlat,xlon,dumb).eq.0) then xlon=-xlon call zangsat(splat,splon,xlat,xlon,satang) cossat(ne)=satang endif endif cc if(cosang(ne).ge.0.05*ioff) then if(cosang(ne).ge.0.015*ioff) then do 420 i=1,isamp if((iconv.eq.1).or.(iconv.eq.6)) then iout(ne+i-1)=sqrt(((float(iout(ne+i-1)))**2/cosang(ne))*corr(ne)) if(iconv.eq.6) then iout(ne+i-1)=iout(ne+i-1)-(23*(1+1.3*ssun(ne)*ssat(ne)))+.5 endif if(iout(ne+i-1).gt.255)iout(ne+i-1)=255 if(iout(ne+i-1).lt.0)iout(ne+i-1)=0 else iout(ne+i-1)=sqrt(255*corr(ne)*(float(iout(ne+i-1))/cosang(ne))) if(iout(ne+i-1).gt.255)iout(ne+i-1)=255 if(iout(ne+i-1).lt.0)iout(ne+i-1)=0 endif 420 continue else do 421 i=1,isamp iout(ne+i-1)=0 421 continue endif endif endif if (iconv.eq.2 ) then iscale = idnint(dbhi) ioff = idnint(dbo(nb)) iout(ne) = look(iout(ne)) endif if (iconv.eq.3 ) then line=ycor+nl*yres if (MOD(ne,isamp).eq.1) then if((nl/isamp)*isamp.eq.nl)then cossat(ne)=1. iele=xcor+(ne+isamp/2)*xres dlin=float(line) dele=float(iele) if(nv1sae(dlin,dele,0.0,xlat,xlon,dumb).eq.0) then xlon=-xlon call zangsat(splat,splon,xlat,xlon,satang) cossat(ne)=satang endif endif do 440 i=1,isamp iout(ne+i-1)=iout(ne+i-1)-(idt*(1-cossat(ne))+.5) if(iout(ne+i-1).gt.255)iout(ne+i-1)=255 if(iout(ne+i-1).lt.0)iout(ne+i-1)=0 440 continue endif endif if(iconv.eq.4) then if(iout(ne).gt.0) then numt=numt+1 if(numt.lt.4) then call sdest('data value=',iout(ne)) call sdest('line=',line) call sdest('ele =',iele) endif iele=xcor+ne*xres+xres/2 line=ycor+nl*yres+yres/2 dlin=float(line) dele=float(iele) status=nv1sae(dlin,dele,0.0,xlat,xlon,dumb) ilat=xlat*1000 ilon=xlon*1000 ihgt=iout(ne)*250 clat=cfu(ilat) clon=cfu(ilon) chgt=cfu(ihgt) cnum=cfu(numt) cout=lf//cr//cnum//clat//clon//chgt c write ascii information to file status=lbo(coutname,numt*35,35,cout) endif endif 300 continue 200 continue c --- pack the pixels to 1 byte values call mpixel(nelems, 4, 1, iout) c --- write this line to the destination file if( mcaout(iout).ne.0 ) then call edest('Failed to write destination file at line: ', nl) goto 1000 endif 100 continue if(iconv.eq.4) then cday=cfu(day) ctime=cfu(time) cout=cday//ctime//cnum//' points' status=lbo(coutname,0,35,cout) call sdest(cout,0) endif c --- write the comment card if( mcacou( card ).ne.0 ) then call edest('Failed to write comment card',0) goto 1000 endif c *************** READ/WRITE IMAGE FILES **************** dlen = len_trim(dest) cdarea = cfu(darea) call sdest('Results written to destination file: '//dest(1:dlen)// & '.'//cdarea,0) call mccodeset(0) c --- All error conditions end up here. 1000 continue call edest('Done',0) c --- end of program return end C $ SUBROUTINE ZANGSAT(SPLAT,SPLON,XLAT,XLON,SATANG) C $ SSANGLES - computes zenith angle of satellite C $ INPUT: C $ SPLAT= (R) latitude of subpoint of satellite C $ SPLON= (R) longitude of subpoint of satellite C $ XLAT = (R) latitude of point C $ XLON = (R) longitude of point (WEST LON NEGATIVE!!) C $ OUTPUT: C $ SATANG = (R) cos of zenith angle of satellite C $$ ANGLES = COMPUTATION, NAVIGATION C ANGLES MOSHER 1074 WINLIB ZENITH ANGLES TO SAT,SUN,AND REL AZIMUTH AN C SUBROUTINE ZANGSAT(SPLAT,SPLON,XLAT,XLON,SATANG) DATA PI/3.14159265/ DATA R/6371.221/ RDPDG=PI/180.0 C DETERMINE SATELLITE POSITION ASSUMING CIRCULAR GEOSTATIONARY ORBIT YSPLAT=RDPDG*SPLAT YSPLON=RDPDG*SPLON XSAT=42164.545*COS(YSPLON) YSAT=42164.545*SIN(YSPLON) ZSAT=0. YLAT=RDPDG*XLAT YLAT=GEOLAT(YLAT,1) YLON=RDPDG*XLON SLAT=SIN(YLAT) CLAT=COS(YLAT) SLON=SIN(YLON) CLON=COS(YLON) XSAM=R*CLAT*CLON YSAM=R*CLAT*SLON ZSAM=R*SLAT C C DETERMINE ZENITH ANGLE OF SATELLITE C XVEC=XSAT-XSAM YVEC=YSAT-YSAM ZVEC=ZSAT-ZSAM XFACT=SQRT(XVEC**2+YVEC**2+ZVEC**2) SATANG=(XVEC*XSAM+YVEC*YSAM+ZVEC*ZSAM)/(R*XFACT) RETURN END SUBROUTINE TTOZ(LAT,IDATE,ICTEMP,IZ,IPMB,LTHICK,IPBOT) C MVTTOZ MOSHER 1074 WINLIB TEMPERATURE TO HEIGHT USING STND ATMOS C C $ CONVERTS TEMPERATURE TO HEIGHT AND PRESSURE. C $ LAT = (I) INPUT LATITUDE C $ IDATE = (I) INPUT DAY C $ ICTEMP = (I) INPUT CLOUD TOP TEMPERATURE IN DEG K*10 C $ LTHICK = (I) INPUT CLOUD THICKNESS IN METERS C $ IPMB = (I) OUTPUT CLOUD TOP PRESSURE IN MILLIBARS C $ IPBOT = (I) OUTPUT PRESSURE AT THE BOTTOM OF THE CLOUD C $ IZ = (I) OUTPUT CLOUD TOP HEIGHT IN KILOMETERS*10 C $ MVTTOZ = MV C* C* FUNCTION C* USES THE MODEL ATMOSPHERES OF THE 1962 STANDARD FOR A TEMPERATURE C* TO HEIGHT AND PRESSURE CONVERSION. HEIGHT IS GIVEN IN KM*10 WHILE C* PRESSURE IS GIVEN IN MB. ALL PRESSURES GREATER THAN 999 ARE SET T C* A LINEAR INTERPOLATION IS USED FOR LATITUDE POSITION WHILE A COSIN C* INTERPOLATION IS USED FOR JULIAN DATE. C* INTEGER T15(4),T30JAN(4),T30JUL(4),T45JAN(4),T45JUL(4),T60JAN(4),T *60JUL(4),T75JAN(4),T75JUL(4) INTEGER Z15(4),Z30JAN(4),Z30JUL(4),Z45JAN(4),Z45JUL(4),Z60JAN(4),Z *60JUL(4),Z75JAN(4),Z75JUL(4) INTEGER P15(4),P30JAN(4),P30JUL(4),P45JAN(4),P45JUL(4),P60JAN(4),P *60JUL(4),P75JAN(4),P75JUL(4) DATA T15/2997,2871,2870,1932/ DATA T30JAN/2872,2812,2162,2032/ DATA T30JUL/3012,2937,2662,2032/ DATA T45JAN/2722,2617,2197,2082/ DATA T45JUL/2942,2852,2612,2157/ C DATA T60JAN/2593,2592,2512,2172/ C MODIFY TO WORK BETTER OVER THE OCEAN WATERS DATA T60JAN/2693,2592,2512,2172/ DATA T60JUL/2872,2602,2252,2251/ DATA T75JAN/2538,2537,2152,2137/ DATA T75JUL/2782,2717,2262,2261/ DATA Z15/0,23,25,165/ DATA Z30JAN/0,20,120,140/ DATA Z30JUL/0,10,60,150/ DATA Z45JAN/0,30,100,126/ DATA Z45JUL/0,20,60,130/ DATA Z60JAN/0,10,35,85/ DATA Z60JUL/0,50,100,101/ DATA Z75JAN/0,15,85,115/ DATA Z75JUL/0,25,95,96/ DATA P15/1013,780,757,103/ DATA P30JUL/1013,905,493,133/ DATA P30JAN/1013,804,203,145/ DATA P45JAN/1018,694,257,171/ DATA P45JUL/1013,802,487,179/ DATA P60JAN/1013,888,636,307/ DATA P60JUL/1010,541,268,269/ DATA P75JAN/1013,828,299,185/ DATA P75JUL/1012,742,284,280/ ILAT=IABS(LAT) IPHASE=0 IF(LAT.LT.0) IPHASE=180 IDAY=IDATE+IPHASE IF(IDAY.GT.365) IDAY=IDAY-365 C C DETERMINE WHICH STANDARD ATMOSPHERE TO USE C ITHICK=IROUND(FLOAT(LTHICK)/100.) DO 10 I=1,5 ISTAND=15*I IDELT=15-(ISTAND-ILAT) IF(ILAT.LT.ISTAND)GO TO(15,30,45,60,75),I 10 CONTINUE GO TO 75 15 CALL INTERP(IDAY ,ICTEMP,0,T15,T15,T15,T15,Z15,Z15,Z15,Z15,P15,P15 *,P15,P15,IZ,IPMB,ITHICK,IPBOT) RETURN 30 CALL INTERP(IDAY ,ICTEMP,IDELT,T15,T15,T30JAN,T30JUL,Z15,Z15,Z30JA *N,Z30JUL,P15,P15,P30JAN,P30JUL,IZ,IPMB,ITHICK,IPBOT) RETURN 45 CALL INTERP(IDAY ,ICTEMP,IDELT,T30JAN,T30JUL,T45JAN,T45JUL,Z30JAN, *Z30JUL,Z45JAN,Z45JUL,P30JAN,P30JUL,P45JAN,P45JUL,IZ,IPMB,ITHICK, *IPBOT) RETURN 60 CALL INTERP(IDAY ,ICTEMP,IDELT,T45JAN,T45JUL,T60JAN,T60JUL,Z45JAN, *Z45JUL,Z60JAN,Z60JUL,P45JAN,P45JUL,P60JAN,P60JUL,IZ,IPMB,ITHICK, *IPBOT) RETURN 75 IF(IDELT.GT.15)IDELT=15 CALL INTERP(IDAY ,ICTEMP,IDELT,T60JAN,T60JUL,T75JAN,T75JUL,Z60JAN, *Z60JUL,Z75JAN,Z75JUL,P60JAN,P60JUL,P75JAN,P75JUL,IZ,IPMB,ITHICK, *IPBOT) RETURN END C $ SUBROUTINE ZSUN(JDAY,JTIME,XLAT,XLON,GHA,DEC,COSSUN) C $ ZSUN - computes COS zenith angles of sun C $ INPUT: C $ JDAY = (I) picture day (YYDDD) C $ JTIME = (I) picture start time C $ XLAT = (R) latitude of point C $ XLON = (R) longitude of point C $ GHA = (R) Greenwich hour angle of sun C $ DEC = (R) declination of sun C $ OUTPUT: C $ COSSUN = (R) COS zenith angle of sun C $$ ANGLES = COMPUTATION, NAVIGATION C ANGLES MOSHER 1074 WINLIB ZENITH ANGLES TO SAT,SUN,AND REL AZIMUTH AN C SUBROUTINE ZSUN(JDAY,JTIME,XLAT,XLON,GHA,DEC,COSSUN) DATA PI/3.14159265/ DATA R/6371.221/ RDPDG=PI/180.0 ZLON=-XLON PICTIM=FTIME(JTIME) C YLAT=RDPDG*XLAT YLAT=GEOLAT(YLAT,1) YLON=RDPDG*ZLON SLAT=SIN(YLAT) CLAT=COS(YLAT) SLON=SIN(YLON) CLON=COS(YLON) XSAM=R*CLAT*CLON YSAM=R*CLAT*SLON ZSAM=R*SLAT C C DETERMINE ZENITH ANGLE OF SUN C SNLG=-PICTIM*PI/12.0-RDPDG*GHA SNDC=RDPDG*DEC COSDEC=COS(SNDC) US=COS(SNLG)*COSDEC VS=SIN(SNLG)*COSDEC WS=SIN(SNDC) COSSUN=(US*XSAM+VS*YSAM+WS*ZSAM)/R RETURN END SUBROUTINE INTERP(IDATE,ICTEMP,IDELT,T1JAN,T1JUL,T2JAN,T2JUL,Z1JAN *,Z1JUL,Z2JAN,Z2JUL,P1JAN,P1JUL,P2JAN,P2JUL,IZ,IPMB,ITHICK,IPBOT) C C $ STANDARD ATMOSPHERE INTERPOLATION C $ INPUT: C $ IDATE = (I) DAY C $ ICTEMP = (I) CLOUD TOP TEMPERATURE C $ IDELT = (I) DISTANCE IN DEGREES FROM THE LOWER LATITUDE C $ T1JAN = (I) JANUARY TEMPERATURE BREAKPOINTS FOR LOWER LATITUDE C $ T1JUL = (I) JULY TEMPERATURE BREAKPOINTS FOR LOWER LATITUDE C $ T2JAN = (I) JANUARY TEMPERATURE BREAKPOINTS FOR UPPER LATITUDE C $ T2JUL = (I) JULY TEMPERATURE BREAKPOINTS FOR UPPER LATITUDE C $ Z1JAN = (I) JANUARY HEIGHT BREAKPOINTS FOR LOWER LATITUDE C $ Z1JUL = (I) JULY HEIGHT BREAKPOINTS FOR LOWER LATITUDE C $ Z2JAN = (I) JANUARY HEIGHT BREAKPOINTS FOR UPPER LATITUDE C $ Z2JUL = (I) JULY HEIGHT BREAKPOINTS FOR UPPER LATITUDE C $ P1JAN = (I) JANUARY PRESSURE BREAKPOINTS FOR LOWER LATITUDE C $ P1JUL = (I) JULY PRESSURE BREAKPOINTS FOR LOWER LATITUDE C $ P2JAN = (I) JANUARY PRESSURE BREAKPOINTS FOR UPPER LATITUDE C $ P2JUL = (I) JULY PRESSURE BREAKPOINTS FOR UPPER LATITUDE C $ ITHICK = (I) THICKNESS OF CLOUD C $ C $ OUTPUT: C $ IZ = (I) HEIGHT C $ IPMB = (I) PRESSURE AT THE TOP OF THE CLOUD C $ IPBOT = (I) PRESSURE AT THE BOTTOM OF THE CLOUD C $ INTERP = COMPUTATION, OBSERVATION, MV C* INTEGER T1JAN(4),T1JUL(4),T2JAN(4),T2JUL(4) INTEGER Z1JAN(4),Z1JUL(4),Z2JAN(4),Z2JUL(4) INTEGER P1JAN(4),P1JUL(4),P2JAN(4),P2JUL(4) COMMON /COSX/COSX C C LINEAR INTERPOLATION BETWEEN BREAK POINTS IN STANDARD ATMOSPHERE C CALL LININT(ICTEMP,T1JAN,Z1JAN,P1JAN,IZ1JAN,IP1JAN,ITHICK,IB1JAN) CALL LININT(ICTEMP,T2JAN,Z2JAN,P2JAN,IZ2JAN,IP2JAN,ITHICK,IB2JAN) CALL LININT(ICTEMP,T1JUL,Z1JUL,P1JUL,IZ1JUL,IP1JUL,ITHICK,IB1JUL) CALL LININT(ICTEMP,T2JUL,Z2JUL,P2JUL,IZ2JUL,IP2JUL,ITHICK,IB2JUL) C C LINEAR INTERPOLATION BETWEEN LATITUDES C IZJAN=(IZ2JAN-IZ1JAN)*IDELT/15+IZ1JAN IZJUL=(IZ2JUL-IZ1JUL)*IDELT/15+IZ1JUL C C COSINE INTERPOLATION BETWEEN DATES C IY=IZJUL-IZJAN Y=(-COSX+1.0)/2.0 IZ=IY*Y+IZJAN RETURN END SUBROUTINE LININT(ITEMP,IT,IH,IP,IZ,IPMB,ITHICK,IBOT) C $ C $ INTERPOLATES BETWEEN BREAK POINTS IN STANDARD ATMOSPHERES C $ ITEMP = (I) INPUT CLOUD TOP TEMPERATURE*10 C $ IT = (I) INPUT TEMPERATURE BREAKPOINTS C $ IH = (I) INPUT HEIGHT BREAKPOINTS C $ IP = (I) INPUT PRESSURE BREAKPOINTS C $ ITHICK = (I) INPUT THICKNESS OF CLOUD C $ IZ = (I) OUTPUT HEIGHT C $ IPMB = (I) OUTPUT PRESSURE AT THE TOP OF THE CLOUD C $ IBOT = (I) OUTPUT PRESSURE AT THE BOTTOM OF THE CLOUD C $ LININT = COMPUTATION, MV C* DIMENSION IT(4),IH(4),IP(4) IF(ITEMP.LT.IT(1))GO TO 1 C C CLOUD TEMPERATURE IS WARMER THAN WARMEST STANDARD ATMOSPHERE TEMPERATU C IZ=0 RETURN C C DETERMINE PROFILE POSITION C 1 DO 2 I=2,4 IF(ITEMP.LT.IT(I))GO TO 2 I1=I I2=I-1 GO TO 3 2 CONTINUE C C CLOUD IS COLDER THAN COLDEST STANDARD TEMPERATURE C I1=4 I2=3 C C DETERMINE PRESSURE AT HEIGHT IZ USING THE HYPSOMETRIC EQUATION C 3 IZ=(IH(I2)-IH(I1))*(ITEMP-IT(I1))/(IT(I2)-IT(I1))+IH(I1) RETURN END C $ SUBROUTINE ANGLESS(JDAY,JTIME,XLAT,XLON,GHA,DEC,SPLAT,SPLON,SATANG,SUNANG,RELANG) C $ ANGLESS - computes zenith angles of sun and satellite and relative C $ azimuth angle (DAS) C $ INPUT: C $ JDAY = (I) picture day (YYDDD) C $ JTIME = (I) picture start time C $ XLAT = (R) latitude of point C $ XLON = (R) longitude of point C $ GHA = (R) Greenwich hour angle of sun C $ DEC = (R) declination of sun C $ SPLAT= (R) Satellite subpoint lat C $ SPLON= (R) Satellite subpoint lon C $ OUTPUT: C $ SATANG = (R) zenith angle of satellite C $ SUNANG = (R) zenith angle of sun C $ RELANG = (R) relative angle C $$ ANGLES = COMPUTATION, NAVIGATION C ANGLES MOSHER 1074 WINLIB ZENITH ANGLES TO SAT,SUN,AND REL AZIMUTH AN C SUBROUTINE ANGLESS(JDAY,JTIME,XLAT,XLON,GHA,DEC,SPLAT,SPLON, &SATANG,SUNANG,RELANG) REAL CLAT, CLON, COSDEC DATA IDAY/0/ DATA PI/3.14159265/ DATA R/6371.221/ RDPDG=PI/180.0 IF(IDAY.EQ.JDAY)GO TO 1 IDAY=JDAY INORB=0 1 PICTIM=FTIME(JTIME) C C DETERMINE SATELLITE POSITION C cc CALL SATPOS(INORB,JTIME,XSAT,YSAT,ZSAT) cc HEIGHT=SQRT(XSAT**2+YSAT**2+ZSAT**2) cc YLAT=RDPDG*XLAT cc YLAT=GEOLAT(YLAT,1) cc YLON=RDPDG*XLON C DETERMINE SATELLITE POSITION ASSUMING CIRCULAR GEOSTATIONARY ORBIT YSPLAT=RDPDG*SPLAT YSPLON=RDPDG*SPLON XSAT=42164.545*COS(YSPLON) YSAT=42164.545*SIN(YSPLON) ZSAT=0. HEIGHT=SQRT(XSAT**2+YSAT**2+ZSAT**2) YLAT=RDPDG*XLAT YLAT=GEOLAT(YLAT,1) YLON=RDPDG*XLON SLAT=SIN(YLAT) CLAT=COS(YLAT) SLON=SIN(YLON) CLON=COS(YLON) XSAM=R*CLAT*CLON YSAM=R*CLAT*SLON ZSAM=R*SLAT C C DETERMINE ZENITH ANGLE OF SUN C SNLG=-PICTIM*PI/12.0-RDPDG*GHA SNDC=RDPDG*DEC COSDEC=COS(SNDC) US=COS(SNLG)*COSDEC VS=SIN(SNLG)*COSDEC WS=SIN(SNDC) SUNANG=ACOS((US*XSAM+VS*YSAM+WS*ZSAM)/R)/RDPDG C C DETERMINE ZENITH ANGLE OF SATELLITE C XVEC=XSAT-XSAM YVEC=YSAT-YSAM ZVEC=ZSAT-ZSAM XFACT=SQRT(XVEC**2+YVEC**2+ZVEC**2) SATANG=ACOS((XVEC*XSAM+YVEC*YSAM+ZVEC*ZSAM)/(R*XFACT))/RDPDG C C DETERMINE RELATIVE ANGLE C X1=CLAT*CLON Y1=CLAT*SLON Z1=SLAT X2=SLON Y2=-CLON X3=-SLAT*CLON Y3=-SLAT*SLON Z3=CLAT XC1=US-X1 YC1=VS-Y1 ZC1=WS-Z1 XC2=XSAT/HEIGHT-X1 YC2=YSAT/HEIGHT-Y1 ZC2=ZSAT/HEIGHT-Z1 XAN1=XC1*X3+YC1*Y3+ZC1*Z3 XAN2=XC2*X3+YC2*Y3+ZC2*Z3 YAN1=XC1*X2+YC1*Y2 YAN2=XC2*X2+YC2*Y2 XAN3=XAN1*XAN2+YAN1*YAN2 YAN3=-YAN1*XAN2+XAN1*YAN2 RELANG=ATAN2(YAN3,XAN3)/RDPDG RELANG=ABS(RELANG) RETURN END