C+ C NAME: C Write3D_nv_UT C PURPOSE: C Writes output from tomographic reconstruction to ascii t3d* file at specified times. C The times are equally separated in UT. C Ouput are the 3D velocity and normalized density arrays, and, optionally, C the source surface maps, and los crossings maps. C C As a secondary function the input arrays G2S and G3D are converted C to normalized density. C CATEGORY: C Tomography: I/O C CALLING SEQUENCE: subroutine Write3D_nv_UT(NT,nCar,JDCar,XCbegMAT,XCendMAT,XCbegROI,XCendROI,VSS,G2S,XC3D,V3D,G3D,Vtmp,Gtmp) C INPUTS: C NT integer # times per day written to separate files C All times inside the time range covered by the reconstruction C rounded to a fraction 1.NT days is written, i.e. C for each day, times (I-1)*24.0/NT hours (I=1,NT) UT is used. C C nCar integer # entries in JDCar C JDCar(nCar) real*8 Start times of Carrington rotation in Julian days C C (except for the scratch arrays all input is read-only) C C XCbegMAT(nTim) real start Carrington variable of matrix (- NCoff) C XCendMAT(nTim) real end Carrington variable of matrix (- NCoff) C (range of matrix typically is three Carrington rotations) C XCbegROI(nTim) real start Carrington variable of region of interest C XCendROI(nTim) real end Carrington variable of region of interest C (the region of interest typically covers one whole C Carrington rotation C VSS(nLng,nLat,nTim) real velocity at source surface RR (with bad values) C G2S(nLng,nLat,nTim) real nr^2 (normalized density) or g^2 at source surface at RR C (with bad values) C XC3D(3,nLng,nLat,nRad,nTim) C real shift from source surface (in Carrington variable (- NCoff)) C V3D(nLng,nLat,nRad,nTim)real velocity ratio or velocity (no bad values) C G3D(nLng,nLat,nRad,nTim)real nr^2 / g^2 ratio or nr^2 / g^2 (no bad values) C C Vtmp(nLng,nLat,nRad) real scratch array C Gtmp(nLng,nLat,nRad) real scratch array C GStmp(nLng,nLat) real scratch array C OUTPUTS: C (to t3d file) C CALLS: C T3D_set, T3D_get, T3D_iset, T3D_iget, T3D_get_grid, Str2Str, FLINT, FLINT8, CvG2D C BadR4, WR2DARR, T3D_write_nv C RESTRICTIONS: C V3D and G3D are not allowed to have bad values (the source surface arrays are used C to determine what is bad). V3D and G3D usually will be output from the last call C to href=SW_Model_Kinematic=. C INCLUDE: include 't3d_param.h' include 't3d_array.h' include 't3d_index.h' C include 't3d_grid_fnc.h' C include 't3d_loc_fnc.h' C PROCEDURE: C > This function should be called prior to href=Write3D_nv= (which will modify its C input arrays V3D and G3D) C C > G2S and V2S contain the source surface maps, which may contain bad values. C V3D and G3D contain velocity, and g^2 or nr^2, information without bad values. C If TOM__RATIO is not set then the source surface data are used C only to flag array elements as bad. C C Several bits in MODE to control several aspects of the calculations. C C TOM__MOD Carrington variables at the source surface are shifted into C the range [XCbeg,XCend] by applying a mod(XCend-XCbeg) operation. C C TOM__G2 If set then input is in g^2 rather than normalized density nr^2 C If set: C G2S is g^2 at source surface (with bad values) C G3D is 3D g^2 matrix (if TOM__RATIO not set) or C g^2 ratio (if TOM__RATIO set) C If not set: C G2S is nr^2 at source surface (with bad values) C G3D is 3D nr^2 matrix (if TOM__RATIO not set) or C nr^2 ratio (if TOM__RATIO set) C C TOM__NORATIO If not set then V3D and G3D are ratios C If set: C V3D is the 3D velocity matrix C G3D is the g^2 matrix (if TOM__G2 set) or the C nr^2 matrix (if TOM__G2 not set) C MODIFICATION HISTORY: C SEP-2001, Paul Hick (UCSD/CASS; pphick@ucsd.edu), base on Write3D_nv C- integer NT integer nCar real*8 JDCar(nCar) real XCbegMAT(*) real XCendMAT(*) real XCbegROI(*) real XCendROI(*) real VSS (*) ! Velocity at source surface real G2S (*) ! G^2 or nr^2 at source surface real XC3D (*) real V3D (*) real G3D (*) real Vtmp (*) real Gtmp (*) real*8 JDb real*8 JDe real*8 dJD real*8 FLINT8 integer NDSS(3) integer ND3D(4) integer NDCC(5) real PPSS(3) real PPCC(5) equivalence (ND3D, NDCC(2)) real Scale(4) /0.0,1.0, 0.0,1.0/ real R_pwr(2) /0.0,2.0/ character cFile *120 character cFmt2D*7 /'(F10.4)'/ character cUser *80 character cBinX *80 /'Velocity/normalized density line of sight crossings at source surface'/ character cShift*80 /'Shift in Carrington variable (add to go back to source surface)'/ integer Str2Str logical bG2 logical bRatio logical bBinX logical bShift logical bMod logical bPointP logical bInclV logical bInclG real Tiny /0.0/ real*8 Tiny8 /0d0/ include 't3d_grid_fnc.h' include 't3d_loc_fnc.h' if (NT .le. 0) return call T3D_iget(T3D__MODE,0,MODE) !------- ! Affect calculation bMod = iand(MODE,TOM__MOD ) .ne. 0 ! Enforce mod(XCend-XCbeg) bG2 = iand(MODE,TOM__G2 ) .ne. 0 ! Input is in G2 (not normalized density) bRatio = iand(MODE,TOM__NORATIO) .eq. 0 ! V3D, G3D are ratios !------- ! Only affect output quantities bBinX = iand(MODE,TOM__BINX ) .ne. 0 ! Add los crossings array to t3d file bShift = iand(MODE,TOM__SHIFT ) .ne. 0 !------- ! Only used to build user string bPointP = iand(MODE,TOM__POINTP ) .ne. 0 ! Identifies initial source surface maps (point-P map or not) bInclV = iand(MODE,TOM__INCLV ) .ne. 0 bInclG = iand(MODE,TOM__INCLG ) .ne. 0 !------- ! Summary of various logical settings is translated into a line ! for the header of the output t3d* file I = 0 I = I+Str2Str('IPSG2',cUser(I+1:))+2 if ( bPointP) I = I+Str2Str('POINTP',cUser(I+1:))+4 if (.not. bPointP) I = I+Str2Str('FLAT' ,cUser(I+1:))+4 if (.not. bInclV ) I = I+Str2Str('NO',cUser(I+1:))+1 I = I+Str2Str('INCLV',cUser(I+1:))+4 if (.not. bInclG ) I = I+Str2Str('NO',cUser(I+1:))+1 I = I+Str2Str('INCLG',cUser(I+1:)) call T3D_set(T3D__SCALE,4,Scale) call T3D_set(T3D__R_PWR,2,R_pwr) call T3D_get_grid(TT,dTT,RR,dRR,nLng,nLat,nRad,nTim,dTTi,nLng1,nLat1) call T3D_get(T3D__PWN_V,0,PWN) NDSS(1) = nLng NDSS(2) = nLat NDSS(3) = nTim ! May be degenerate (=1) !------- ! This also fills ND3D (equivalenced to NDCC) NDCC(1) = PRJ__N ! Dimensions of XC3D array NDCC(2) = nLng ! The last 4 entries are dimensions of V3D, G3D NDCC(3) = nLat NDCC(4) = nRad NDCC(5) = nTim ! May be degenerate (=1) Bad = BadR4() dJD = 1d0/NT ! Time increment in days JDb = TTtime( 1) ! First and last time in Carrington variable JDe = TTtime(nTim) JDb = FLINT8(1,nCar,JDCar,JDb,Tiny8) ! First and last time in Julian dates JDe = FLINT8(1,nCar,JDCar,JDe,Tiny8) JDb = int(JDb)+( 1+int((JDb-int(JDb))/dJD) )*dJD! First time rounded to dJD (inisde [JDb,JDe]) JDe = int(JDe)+( int((JDe-int(JDe))/dJD) )*dJD! Last time rounded to dJD (inisde [JDb,JDe]) nTim0 = 1+nint((JDe-JDb)/dJD) ! # times in equal steps of UT call T3D_iset(T3D__MODE,0,ior(MODE,TOM__NOSURF))! Switch off writing of source surface arrays ! Old MODE needs to be restored before returning) do L=0,nTim0-1 !------- ! Fill in time-dependent parameters in t3d array call Julian(1,iYr,Doy,JDb+L*dJD,JDe) ! Time: Julian day --> iYr,Doy TL = XMAP_SC_POS(EARTH,iYr,Doy,nCar,JDCar) ! Time (Carrington time) call T3D_iset(T3D__ITIME ,0,0 ) call T3D_set (T3D__TIME ,0,TL) TL = TTindx(TL) ! Time index XCbegL = FLINT(1,nTim,XCbegMAT,TL,Tiny) ! Matrix limits at time TL XCendL = FLINT(1,nTim,XCendMAT,TL,Tiny) call T3D_set (T3D__XCMAT ,0,XCbegL) call T3D_set (T3D__XCMAT+1,0,XCendL) call T3D_set (T3D__XCROI ,0,FLINT(1,nTim,XCbegROI,TL,Tiny)) call T3D_set (T3D__XCROI+1,0,FLINT(1,nTim,XCendROI,TL,Tiny)) PPCC(5) = TL ! Time index do K=1,nRad PPCC(4) = float(K) ! Radial index do J=1,nLat PPCC(3) = float(J) ! Latitude index XL = cosd(XLdeg(J)) do I=1,nLng PPCC(2) = float(I) ! Longitude index !------ ! PPCC contains the position for interpolation on the 5-dim XC3D array. ! The first four are set to integer array indices (converted to float). ! Only the last (time) dimension will be interpolated. ! The 4-dim location PPCC(2:5) is used to interpolate on the 4-dim V3D and G3D ! arrays (remember that ND3D is equivalenced to NDCC(2:5). ! Note that we assume that there are no bad values in V3D, G3D and XC3D. V3D_ = FLINT(5,ND3D,V3D,PPCC(2),Tiny)! Velocity in grid point I,J,K at time TL G3D_ = FLINT(5,ND3D,G3D,PPCC(2),Tiny)! G2-level in grid point I,J,K at time TL PPCC(1) = PRJ__XC ! Longitude shift index XCshift = FLINT(5,NDCC,XC3D,PPCC,Tiny)! Shift to SS in XC from grid point I,J,K at time TL PPCC(1) = PRJ__TT ! Time shift index TTshift = FLINT(5,NDCC,XC3D,PPCC,Tiny)! Shift to SS in TT from grid point I,J,K at time TL !------- ! Find the Carrington variable of origin at the source surface of grid ! point I,J,K at time TL Note that XCbeg and XCrange are set to the ! values at time L. XCbeg = XCbegL ! Matrix limits at time TL: needed for statement fnc XCvar XCend = XCendL XCrange = XCend-XCbeg PPSS(1) = XCvar(I)+XCshift ! Carrington variable at source surface !------- ! Find the time of origin at the source surface of grid point I,J,K,L, PPSS(3) = TL+TTshift ! Time at source surface PPSS(3) = TTindx(PPSS(3)) ! Index for time in source surface array !------ ! To convert the Carrington longitude to an array index, the values of ! XCbeg and XCrange need to be set to the proper values at time PPSS(3) ! (the time of origin at the source surface). ! This is done by interpolation on the XCbegMAT and XCendMAT arrays. XCbeg = FLINT(1,nTim,XCbegMAT,PPSS(3),Tiny)! Matrix limits at time TL+TTshift: XCend = FLINT(1,nTim,XCendMAT,PPSS(3),Tiny)! .. needed for statement fncs XCindx, XCmod if (XCbeg .eq. Bad .or. XCend .eq. Bad) then V = Bad D = Bad else XCrange = XCend-XCbeg if (bMod) PPSS(1) = XCmod(PPSS(1)) ! Map to [XCbeg, XCend] if bMod is set PPSS(1) = XCindx(PPSS(1)) ! Index for longitude in source surface array PPSS(2) = float(J) ! Latitude index !------ ! Find V and g^2 at source surface location (longitude PPSS(1), latitude ! PPSS(2), time PPSS(3)). There may be bad values in the source surface ! VSS, G2S arrays. V = FLINT(-3,NDSS,VSS,PPSS,Tiny) ! Velocity at source surface D = FLINT(-3,NDSS,G2S,PPSS,Tiny) ! nr^2 or g^2 at source surface if (bRatio) then ! V3D, G3D contain ratios if (V .ne. Bad) V = V3D_*V if (D .ne. Bad) D = G3D_*D else ! V3D, G3D do not contain ratios if (V .ne. Bad) V = V3D_ if (D .ne. Bad) D = G3D_ end if end if iloc = locVOL(I,J,K,1) ! Location in nLng x nLat x nRad array Vtmp(iloc) = V ! Velocity Gtmp(iloc) = D ! nr^2 or g^2 end do ! END: longitude loop end do ! END: latitude loop end do ! END: radial loop !------- ! At this point Vtmp is the 3D velocity matrix. Gtmp is the nr^2 or g^2 matrix. if (bG2) call CvG2D(PWN,-nLng*nLat*nRad,Gtmp,Gtmp) ! Convert to nr^2 !------- ! Write velocity and nr^2. Note that we don't write any source surface arrays. call T3D_write_nv('nvut',t3d,t3d,t3d,Vtmp,Gtmp,1,cUser,cFile) end do call T3D_iset(T3D__MODE,0,MODE) ! Restore old mode return end