;**************************************************************************** ; Author : Josh Grant, COAPS ; ; Description : This program is used to calculate divergence at each valid ; point in the NSCAT dataset. The data is read in by calling ; read_NSCAT_winds. The array of latitudes, longitudes, winds ; directions, wind speeds, and dimensions of the arrays are ; returned. From here the divergence of every point in the ; dataset is calculated using the same methodology used by ; Sonya C. Kulkarni in her paper "Variablility of Surface Wind ; Convergence", also known as the "Doughnut method". ; ; **Warning** : It is highly adviced that this code by optimized if it is to ; be used in the future. The origninal code was written in ; FORTRAN and later written in IDL, but it was never optimized ; properly for IDL. ;**************************************************************************** @read_NSCAT_winds.pro ;************************************************************************** ; Author : Josh Grant, COAPS ; Procedure Name : MINMAX ; Description : Given an array of latitudes and longitudes, the minimum ; and maximum of both arrays is determined. This is used ; to determine starting and stopping points of each sector ; in the doughnut. ; ; amax = maximum latitude ; amin = minimum latitude ; omax = maximum longitude ; omin = minimum longitude ; ; amax_ ; *' | `* ; .' | `. ; / | \ ; / 4 | 1 \ ; * | * ; omin---------.--------omax ; * | * ; \ 3 | 2 / ; \ | / ; `. | .' ; `.__amin_' ; ; The numbers correspond to the different sectors. ; ; ***If someone else feels they could make a better circle out of ASCII ; text then go right ahead. I didn't want to spend too much time on ; this. ;************************************************************************** PRO MINMAX, alat, alon, np, amax, amin, omax, omin amax = alat( 0 ) amin = alat( 0 ) omax = alon( 0 ) omin = alon( 0 ) FOR loop = 0, np - 1 DO BEGIN IF ( alat( loop ) LT amin ) THEN amin = alat( loop ) IF ( alat( loop ) GT amax ) THEN amax = alat( loop ) IF ( alon( loop ) LT omin ) THEN omin = alon( loop ) IF ( alon( loop ) GT omax ) THEN omax = alon( loop ) ENDFOR END ;************************************************************************** ; Author : Josh Grant, COAPS ; Procedure Name : ORDER ; Description : Given an array of latitudes, longitudes, wind ; magnitudes, wind directions the angle from the ; line drawn out from points (clat, clon) and (at, al) is ; calculated. Then each point is ordered by that angle, ; with the smallest angle being first. ; ; n = total number of points in sector, this is thus the ; the total valid points in each array. ; la = all the latitudes in the sector. ; lo = all the longitudes in the sector. ; ma = all the wind magnitudes in the sector. ; com = all the wind directions in the sector. ; clat= the latitude value of the center of the doughnut ; clon= the longitude value of the center of the doughnut ; at = the latitude value of the sector starting point. ; al = the longitude value of the sector starting point. ; _____ ; *' | `* ; .' | `. ; / | \ ; / 4 | 1 \ ; * | * ; |---------.--------| ; * | * ; \ 3 | 2 / ; \ | / ; `. | .' ; `.___|___' ; ; Each sector is a quadrant in the circle. If sector ; 1 was being ordered then the arrays would contain all ; longitudes, latitudes, wind magnitudes, and wind ; directions within that sector. clat and clon would be ; the center of the doughnut (this is the same no matter ; which sector is being ordered). Most important though ; is that (at, al) would equal the point at the top of the ; doughnut (at = maximum latitude, al = center point ; longitude. ; ; ***If someone else feels they could make a better circle out of ASCII ; text then go right ahead. I didn't want to spend too much time on ; this. ;************************************************************************** PRO ORDER, n, la, lo, ma, com, clat, clon, at, al angle = FLTARR( 300 ) angle( * ) = 0.0 dumang = FLTARR( 300 ) dumang( * ) = 0.0 rad = 180.0 / ( ACOS( -1 ) ) aladum = la alodum = lo amadum = ma icomdum = com c = SQRT( ( clat - at ) ^ 2 + ( clon - al ) ^ 2 ) FOR m = 0, n - 1 DO BEGIN a = SQRT( ( la( m ) - at ) ^ 2 + ( lo( m ) - al ) ^ 2 ) b = SQRT( ( la( m ) - clat ) ^ 2 + ( lo( m ) - clon ) ^ 2 ) cosA = ( ( b * b ) + ( c * c ) - ( a * a ) ) / ( 2.0 * b * c ) IF ( cosA GT 1.0 ) THEN BEGIN cosA = 0.9999999 ENDIF IF ( cosA LT -1.0 ) THEN BEGIN cosA = -0.9999999 ENDIF angle( m ) = ACOS( cosA ) * rad ENDFOR dumang = angle FOR m = 0, n - 2 DO BEGIN FOR l = m + 1, n - 1 DO BEGIN IF ( angle( m ) GT angle( l ) ) THEN BEGIN temp = angle( l ) angle( l ) = angle( m ) angle( m ) = temp ENDIF ENDFOR ENDFOR FOR m = 0, n - 1 DO BEGIN FOR l = 0, n - 1 DO BEGIN IF ( angle( l ) EQ dumang( m ) ) THEN BEGIN la( l ) = aladum( m ) lo( l ) = alodum( m ) ma( l ) = amadum( m ) com( l ) = icomdum( m ) ENDIF ENDFOR ENDFOR END ;************************************************************************** ; Author : Josh Grant, COAPS ; Procedure Name : DTHETA ; Description : Given an array of latitudes, longitudes, wind ; magnitudes, and wind directions of an ordered doughnut. ; The area in between every value is computed. Each area ; can be thought of as a wedge in the doughnut. These ; wedges are added together to return a total area. ; ; n = total number of points in the doughnut, and is ; thus the total number of points in each array. ; la = all the latitudes in the doughnut ; lo = all the longitudes in the doughnut ; ma = all the wind magnitudes in the doughnut ; com = all the wind directions in the doughnut ; clat= the latitude value of the center of the doughnut ; clon= the longitude value of the center of the doughnut ; area= the total area of the doughnut ; al = the longitude value of the sector starting point. ; ; When the area of a wedge is calculated... ; ; c ; ._________________________________ ; -----_____ | ; -----_____ | a ; b -----_____ | ; --- ; a, b, and c are calculated by finding their distances. ; Next by using the laws of cosine the angle swept out by ; the wedge is calculated and then the area. At the end ; all the areas are added together. ; ;************************************************************************** PRO DTHETA, n, la, lo, ma, com, clat, clon, area, theta wedge = FLTARR( 300 ) theta[ * ] = 0.0 rad = 180.0 / ( ACOS( -1 ) ) check = 0.0 area = 0.0 FOR k = 0, n - 2 DO BEGIN a = SQRT((la( k ) - la( k + 1 )) ^ 2 + (lo( k ) - lo( k + 1 )) ^ 2) b = SQRT((clat - la( k + 1 )) ^ 2 + ( clon - lo( k + 1 )) ^ 2) c = SQRT((la( k ) - clat ) ^ 2 + ( lo( k ) - clon ) ^ 2) cosA = (( b * b ) + ( c * c ) - ( a * a )) / ( 2.0 * b * c ) IF ( cosA GT 1.0 ) THEN cosA = 0.9999999 IF ( cosA LT -1.0 ) THEN cosA = -0.999999 theta( k ) = ACOS( cosA ) * rad tmp = ( ( la( k ) - clat ) * ( lo( k + 1 ) - clon ) ) - $ ( ( lo( k ) - clon ) * ( la( k + 1 ) - clat ) ) IF ( tmp LT 0.0 ) THEN theta( k ) = 360.0 - theta( k ) check = check + theta( k ) wedge( k ) = 0.5 * ABS( b * c * sin( theta( k ) / rad ) ) area = area + wedge( k ) ENDFOR a = SQRT((la( n - 1 ) - la( 0 )) ^ 2 + (lo( n - 1 ) - lo( 0 )) ^ 2) b = SQRT((clat - la( 0 )) ^ 2 + ( clon - lo( 0 )) ^ 2) c = SQRT((la( n - 1 ) - clat ) ^ 2 + ( lo( n - 1 ) - clon ) ^ 2) cosA = (( b * b ) + ( c * c ) - ( a * a )) / ( 2.0 * b * c ) IF ( cosA GT 1.0 ) THEN cosA = 0.9999999 IF ( cosA LT -1.0 ) THEN cosA = -0.999999 theta( n - 1 ) = ACOS( cosA ) * rad tmp = ( ( la( n - 1 ) - clat ) * ( lo( 0 ) - clon ) ) - $ ( ( lo( n - 1 ) - clon ) * ( la( 0 ) - clat ) ) IF ( tmp LT 0.0 ) THEN theta( n - 1 ) = 360.0 - theta( n - 1 ) wedge( n - 1 ) = 0.5 * ABS( b * c * sin( theta( n - 1 ) / rad ) ) area = area + wedge( n - 1 ) END ;************************************************************************** ; Author : Josh Grant, COAPS ; Procedure Name : DIVER ; Description : Once the area of the doughnut has been found, it is time ; to calculate the divergence. The divergence in each ; wedge is calculated and the all the values are added ; together. ; ;************************************************************************** PRO DIVER, n, alat, alon, amag, idir, clat, clon, area, div, totalarea, $ count, dr2 earth = 6.37 rad = 180.0 / ( ACOS( -1 ) ) trlat = earth * sin( 1.0 / rad ) trlon = trlat * cos( clat / rad ) wedgediv = FLTARR( 300 ) wedgediv( * ) = 0.0 area = area * trlat * trlon div = 0.0 FOR k = 0, n - 2 DO BEGIN mag = ( amag( k ) + amag( k + 1 ) ) * 0.5 dir = ( idir( k ) + idir( k + 1 ) ) * 0.5 u = mag * sin( dir / rad ) v = mag * cos( dir / rad ) xnorm = ( alat( k ) - alat( k + 1 ) ) * trlat ynorm = ( alon( k + 1 ) - alon( k ) ) * trlon wedgediv( k ) = ( u * xnorm + v * ynorm ) div = div - wedgediv( k ) ENDFOR mag = ( amag( n - 1 ) + amag( 0 ) ) * 0.5 dir = ( idir( n - 1 ) + idir( 0 ) ) * 0.5 u = mag * sin( dir / rad ) v = mag * cos( dir / rad ) xnorm = ( alat( n - 1 ) - alat( 0 ) ) * trlat ynorm = ( alon( 0 ) - alon( n - 1 ) ) * trlon wedgediv( n - 1 ) = ( u * xnorm + v * ynorm ) div = div - wedgediv( n - 1 ) div = div / area END ;************************************************************************** ; Author : Josh Grant, COAPS ; Procedure Name : MAIN ; Description : For each point in the dataset not equal to missing the ; divergence is calculated. When a point is valid every ; point around the current point is checked for validity. ; If not missing and its distance from the current point ; also known as the center point, is greater than dr1 and ; less than dr2 then it is added to the list of points. ; Once all points have been added the maximum and minimum ; latitude and longitude values are determined. Next the ; points are ordered such that point3 in the arrays has ; a greater angle than point2 from the North line. Also ; in order to properly order the points they must be ; separated into sectors like below. ; North Line ; | ; | ; _____ ; *' | `* ; .' | `. ; / | \ ; / 4 | 1 \ ; * | * ; |---------.--------| ; * | * ; \ 3 | 2 / ; \ | / ; `. | .' ; `.___|___' ; ; Once all the sectors are ordered they are then combined ; so that all the points are in one array. Then the area ; of each doughnut is determined and then the ; divergence. ; ; nx = x dimension of arrays lat, lon, dir, spd ; ny = y dimension of arrays lat, lon, dir, spd ; lat = latitudes of each point ; lon = longitudes of each point ; dir = wind direction of each point ; spd = wind magnitude of each point ; time = time of each swath ; dr1 = inner radius of doughnut ; dr2 = outer radius of doughnut ; diverge = divergence of each point, has dimensions ; nx and ny ; ; ;************************************************************************** miss = 99999 dividend = 4 modby = 5 number = 175 + dividend limit = 180 FOR year = 1997, 1997 DO BEGIN FOR num = 179, limit DO BEGIN IF ( ( num MOD modby ) EQ dividend ) THEN BEGIN stringFormat = STRCOMPRESS( '(i' + STRING( 3 ) + '.' $ + STRING( 3 ) + ')', /REMOVE_ALL ) nscatDirectory = '/Net/NSCAT/nscat/data/nc_W25/NSCAT.' file_name = nscatDirectory + STRTRIM( year, 1 ) + $ STRING( number, format = stringFormat ) + '.nc' dataDirectory = '/Net/NSCAT/people/grant/divfiles/' diverFile = STRTRIM(year,1) + STRING(number,format = stringFormat) $ + '.div' number = number + modby time = LONG(miss) lat = miss lon = miss dir = miss spd = miss flag = read_NSCAT_winds( file_name, miss, time, lat, lon, dir, spd, x, $ y ) IF ( flag GE 0 ) THEN BEGIN ; print, 'Time: ', time(55) ; print, 'Lats: ' ; print, lat(*,55) ; print, 'Lons: ' ; print, lon(*,55) ; print, 'Speeds: ' ; print, spd(*,55) ; print, 'Dirs: ' ; print, dir(*,55) diverge = spd diverge( *, * ) = 0.0 nx = x ny = y ; print, nx ; print, ny count = 0 totalarea = 0 dr1 = 1.25 dr2 = 2.25 nd = 300 alat = FLTARR( 300 ) alon = FLTARR( nd ) amag = FLTARR( nd ) idir = INTARR( nd ) alatmax = 0 alatmin = 0 alonmin = 0 alonmax = 0 np = INTARR( 4 ) ala = FLTARR( 4, nd ) alo = FLTARR( 4, nd ) ama = FLTARR( 4, nd ) icom = FLTARR( 4, nd ) area = 0.0 theta = FLTARR( 300 ) error = 0 FOR j = 0, ny - 1 DO BEGIN FOR i = 0, nx - 1 DO BEGIN clat = lat( i, j ) clon = lon( i, j ) IF ( spd( i, j ) LE 0.0 ) THEN BEGIN error = 1 diverge( i, j ) = miss ; ENDIF ; IF ( error EQ 1 OR (spd(i,j) EQ miss AND dir(i,j) EQ miss) ) $ ; THEN BEGIN ; diverge( i, j ) = miss ENDIF ELSE BEGIN IF ( spd(i,j) EQ miss AND dir(i,j) EQ miss ) THEN BEGIN diverge( i, j ) = miss ENDIF ELSE BEGIN npts = 0 IF ( i LT 20 ) THEN xbegin = 0 ELSE xbegin = i - 20 IF ( i GT nx - 21 ) THEN xend = nx - 1 ELSE xend = i + 20 IF ( j LT 20 ) THEN ybegin = 0 ELSE ybegin = j - 20 IF ( j GT ny - 21 ) THEN yend = ny - 1 ELSE yend = j + 20 FOR x = xbegin, xend DO BEGIN FOR y = ybegin, yend DO BEGIN IF ( spd( x, y ) GT 0.0 AND ( spd(x,y) NE miss AND $ dir(x,y) NE miss ) ) THEN BEGIN rr = SQRT(( clat - lat( x, y )) ^ 2 + $ ( clon - lon( x, y )) ^ 2 ) IF ( (rr GE dr1) AND (rr LE dr2) ) THEN BEGIN alat( npts ) = lat( x, y ) alon( npts ) = lon( x, y ) amag( npts ) = spd( x, y ) idir( npts ) = FIX( dir( x, y ) ) npts = npts + 1 ENDIF ENDIF ENDFOR ENDFOR IF( npts LT 35 ) THEN BEGIN error = 2 diverge( i, j ) = miss ENDIF ELSE BEGIN ; print, npts MINMAX, alat, alon, npts, alatmax, alatmin, alonmax, alonmin np( * ) = 0 FOR k = 0, npts - 1 DO BEGIN IF ( ( alat( k ) GT clat ) AND $ ( alon( k ) GE clon ) ) THEN BEGIN ala( 0, np( 0 ) ) = alat( k ) alo( 0, np( 0 ) ) = alon( k ) ama( 0, np( 0 ) ) = amag( k ) icom( 0, np( 0 ) ) = idir( k ) np( 0 ) = np( 0 ) + 1 ENDIF ENDFOR temp1 = ala( 0, * ) temp2 = alo( 0, * ) temp3 = ama( 0, * ) temp4 = icom( 0, * ) ORDER, np( 0 ), temp1, temp2, temp3, $ temp4, clat, clon, alatmax, clon ala( 0, * ) = temp1 alo( 0, * ) = temp2 ama( 0, * ) = temp3 icom( 0, * ) = temp4 FOR k = 0, npts - 1 DO BEGIN IF ( ( alat( k ) LE clat ) AND $ ( alon( k ) GT clon ) ) THEN BEGIN ala( 1, np( 1 ) ) = alat( k ) alo( 1, np( 1 ) ) = alon( k ) ama( 1, np( 1 ) ) = amag( k ) icom( 1, np( 1 ) ) = idir( k ) np( 1 ) = np( 1 ) + 1 ENDIF ENDFOR temp1 = ala( 1, * ) temp2 = alo( 1, * ) temp3 = ama( 1, * ) temp4 = icom( 1, * ) ORDER, np( 1 ), temp1, temp2, temp3, $ temp4, clat, clon, clat, alonmax ala( 1, * ) = temp1 alo( 1, * ) = temp2 ama( 1, * ) = temp3 icom( 1, * ) = temp4 ; ORDER, np( 1 ), ala( 1, * ), alo( 1, * ), ama( 1, * ), $ ; icom( 1, * ), clat, clon, clat, alonmax FOR k = 0, npts - 1 DO BEGIN IF ( ( alat( k ) LT clat ) AND $ ( alon( k ) LE clon ) ) THEN BEGIN ala( 2, np( 2 ) ) = alat( k ) alo( 2, np( 2 ) ) = alon( k ) ama( 2, np( 2 ) ) = amag( k ) icom( 2, np( 2 ) ) = idir( k ) np( 2 ) = np( 2 ) + 1 ENDIF ENDFOR temp1 = ala( 2, * ) temp2 = alo( 2, * ) temp3 = ama( 2, * ) temp4 = icom( 2, * ) ORDER, np( 2 ), temp1, temp2, temp3, $ temp4, clat, clon, alatmin, clon ala( 2, * ) = temp1 alo( 2, * ) = temp2 ama( 2, * ) = temp3 icom( 2, * ) = temp4 ; ORDER, np( 2 ), ala( 2, * ), alo( 2, * ), ama( 2, * ), $ ; icom( 2, * ), clat, clon, alatmin, clon FOR k = 0, npts - 1 DO BEGIN IF ( ( alat( k ) GE clat ) AND $ ( alon( k ) LT clon ) ) THEN BEGIN ala( 3, np( 3 ) ) = alat( k ) alo( 3, np( 3 ) ) = alon( k ) ama( 3, np( 3 ) ) = amag( k ) icom( 3, np( 3 ) ) = idir( k ) np( 3 ) = np( 3 ) + 1 ENDIF ENDFOR temp1 = ala( 3, * ) temp2 = alo( 3, * ) temp3 = ama( 3, * ) temp4 = icom( 3, * ) ORDER, np( 3 ), temp1, temp2, temp3, $ temp4, clat, clon, clat, alonmin ala( 3, * ) = temp1 alo( 3, * ) = temp2 ama( 3, * ) = temp3 icom( 3, * ) = temp4 ; ORDER, np( 3 ), ala( 3, * ), alo( 3, * ), ama( 3, * ), $ ; icom( 3, * ), clat, clon, clat, alonmin counter = 0 IF ( (np( 0 ) + np( 1 ) + np( 2 ) + np( 3 )) $ NE npts ) THEN print, 'Error with points' ELSE BEGIN counter = 0 FOR sector = 0, 3 DO BEGIN FOR point = 0, np( sector ) - 1 DO BEGIN alat( counter ) = ala( sector, point ) alon( counter ) = alo( sector, point ) amag( counter ) = ama( sector, point ) idir( counter ) = icom( sector, point ) counter = counter + 1 ENDFOR ENDFOR ENDELSE IF ( npts NE counter ) THEN print, 'Error with counter' DTHETA, npts, alat, alon, amag, idir, clat, clon, area, theta error = 0 FOR k = 0, npts - 1 DO BEGIN ; print, theta(k ) IF ( theta( k ) GT 160 ) THEN BEGIN error = 3 k = npts ENDIF ENDFOR IF ( error EQ 3 ) THEN diverge( i, j ) = miss ELSE BEGIN DIVER, npts, alat, alon, amag, idir, clat, clon,area, $ div, totalarea, count, dr2 diverge( i, j ) = div ENDELSE ENDELSE ENDELSE ENDELSE ; IF ( diverge[i,j] EQ miss ) THEN BEGIN ; print, FORMAT = '(2F9.2,1F14.5,I9)', lat[i,j], lon[i,j], $ ; diverge[i,j],time[j] ; ENDIF ELSE BEGIN ; print, FORMAT = '(2F9.2,1F14.5,I9)', lat[i,j], lon[i,j], $ ; diverge[i,j] * 1000000, time[j] ; ENDELSE ; print, 'divergence = ' + STRTRIM( diverge( i, j ), 1 ) ; print, 'speed = ' + STRTRIM( spd( i, j ), 1 ) ; print, 'direction = ' + STRTRIM( dir( i, j ), 1 ) ; print, 'latitude = ' + STRTRIM( lat( i, j ), 1 ) ; print, 'longitude = ' + STRTRIM( lon( i, j ), 1 ) ; print, ' ' ; print, 'j = ' + STRTRIM(j, 1) + ' is done' ; print, ' ' ENDFOR ; IF ( (j MOD 100) EQ 0 ) THEN print, 'j = ' + STRTRIM(j,1) + ' is done' ; print, 'j = ' + STRTRIM(j,1) ENDFOR OPENW, data, dataDirectory + diverFile, /GET_LUN ; OPENW, data, dataFile, /GET_LUN ; FOR y = 0, ny - 1 DO BEGIN ; printf, diver, FORMAT = '(48F10.4)', diverge[ *, y ] ; ENDFOR FOR y = 0, ny - 1 DO BEGIN FOR x = 0, nx - 1 DO BEGIN printf, data, FORMAT = '(2F9.2,1F13.5,1I9)', lat[x,y], lon[x,y],$ diverge[x,y],time[y] ENDFOR ENDFOR FREE_LUN, data ; FREE_LUN, diver ; ; OPENW, output, dataDirectory + latFile, /GET_LUN ; FOR y = 0, ny - 1 DO BEGIN ; printf, output, FORMAT = '(48F9.2)', lat[ *, y ] ; ENDFOR ; FREE_LUN, output ; ; OPENW, output, dataDirectory + lonFile, /GET_LUN ; FOR y = 0, ny - 1 DO BEGIN ; printf, output, FORMAT = '(48F9.2)', lon[ *, y ] ; ENDFOR ; FREE_LUN, output ; ENDIF ENDIF ENDFOR ENDFOR END