C+ C NAME: C MkVMaptdN0n_inm C PURPOSE: C Change the velocity map from the line of sight (L.O.S.) values of velocity and C weights by the ratio of observed to model values for each line of sight. C The big step taken to make real and smooth the numbers of points that C constitude convergence was also included for modifications to this program. C This smoothing is now done using gridsphere and timesmooth. C C CATEGORY: C Data processing C CALLING SEQUENCE: C call MkVMaptdn0n_inm(LON,LAT,ITIM,RP,LON_in,LAT_in,ITIM_in,RP_in,XEV,VWTij,VWTij_in,FIX,NL,N_IN,NLOS,NLOSP1, C nLng,nLat,nT,VMAP,ANLIMIT, VLIM,VSSIG,SIGB,NS,NBADSV,SIGMAP, C ANMAP,AWTMAP,VMA,FIXM,WTSM,VMEANF,VFIXM2,WTP,WTP_in,Alng,ConstL,CONR,CONT,aNday,Ztmp) C C INPUTS: C LON (NLOSP1,NL) integer Projected Carr. rot. value of point on L.O.S. C LAT (NLOSP1,NL) integer Heliographic lat. point on L.O.S. C ITIM(NLOSP1,NL) integer Projected time values C RP (NLOSP1,NL) real Distance above Sun of point on L.O.S. C XEV(NL) real Source Elongation C VWTij(NLOS,NL) real Weights of each point along the L.O.S. C FIX(NL+I_IN) real Ratio observed/model V for total L.O.S. plus in-situ C LON_in (I_IN) integer Projected Carr. rot. value of in-situ point C LAT_in (I_IN) integer Heliographic lat. of in-situ point C ITIM_in (I_IN) integer Projected time values of in-situ point C RP_in (I_IN) real Distance above Sun of in-situ point C VWTij_in(I_IN) real Weights of each in-situ point C NL integer Number of remote-sensing velocity data points C N_IN integer Number of in-situ data points C NLOS integer Number of L.O.S. distance segments C NLOSP1 integer Number of L.O.S. distance segments + 1 C nLng integer Number of longitude points in VMAP C nLat integer Number of latitude points in VMAP C nT integer Number of time intervals C VMAP(nLng,nLat,nT) real Velocity map C SIGMAP real Velocity map convergence criteria C ANLIMIT integer Number of points constituting deconvolution C VLIM real Limit on velocity (maximum) C VSSIG(NL+N_IN) real Deviation of source from mean in sigma C VM(NL+N_IN) real Model value C SIGB real Velocity sigma cut-off C NS integer 1 = Activate sigma cut-off, 0 = don't C NBADSV(NL+N_IN) integer Bad source velocity < 1 - bad, 1 - good C OUTPUTS: C VMAP(nLng,nLat,nT) real Velocity map C NBADSV(NL+N_IN) integer Bad source velocity -1 = kinematic bad, -2 = mhd bad, 1 - good, set by NS, -NS = MHD C SIGMAP real Velocity map convergence criteria C AWTMAP(nLng,nLat,nT) real Combination Wt and LOS crossing number map C SCRATCH ARRAYS: C ANMAP(nLng,nLat,nT) integer C VMA(nLng,nLat,nT) real C FIXM(nLng,nLat,nT) real C WTSM(nLng,nLat,nT) real C VMEANF(nLng,nLat,nT) real C VFIXM2(nLng,nLat,nT) real C WTP(NLOS,NL) real C WTP_in(N_IN) real C Ztmp(nLng,nLat,nT) real C PROCEDURE: C MODIFICATION HISTORY: C NOV, 1995 B. Jackson (STEL,UCSD) C- subroutine MkVMaptdn0n_inm(LON,LAT,ITIM,RP,LON_in,LAT_in,ITIM_in,RP_in,XEV,VWTij,VWTij_in, & FIX,NL,N_IN,NLOS,NLOSP1, & nLng,nLat,nT,VMAP,ANLIMIT,VLIM,VSSIG,SIGB,NS,NBADSV,SIGMAP, & ANMAP,AWTMAP,VMA,FIXM,WTSM,VMEANF,VFIXM2,WTP,WTP_in, & Alng,ConstL,CONR,CONT,aNday,Ztmp) C integer LON(NLOSP1,NL), ! Projected Carr. rot. value of point on L.O.S. & LAT(NLOSP1,NL), ! Heliographic lat. (in deg.) point on L.O.S. & ITIM(NLOSP1,NL), ! Projected time values & LON_in(N_IN), ! Projected Carr. rot. value of in-situ point. & LAT_in(N_IN), ! Heliographic lat. (in deg.) of in-situ point & ITIM_in(N_IN), ! Projected time values of in-situ point & NBADSV(NL+N_IN) ! Bad source velocity < 1 - bad, 1 - good real RP (NLOSP1,NL), ! Distance above Sun of point on L.O.S. & VWTij(NLOS,NL), ! Weights of each point along the L.O.S. & RP_in (N_IN), ! Distance above Sun of point on L.O.S. & VWTij_in(N_IN), ! Weights of each point along the L.O.S. & VMAP(nLng,nLat,nT), ! Velocity map & FIX (NL+N_IN), ! Ratio observed/model G for total L.O.S. & XEV (NL), ! Source elongation & VSSIG(NL+N_IN) ! Deviation of source from mean in sigma C Scratch files real ANMAP(nLng,nLat,nT), ! Number of points per bin & AWTMAP(nLng,nLat,nT), & VMA(nLng,nLat,nT), & FIXM(nLng,nLat,nT), & WTSM(nLng,nLat,nT), & VMEANF(nLng,nLat,nT), & VFIXM2(nLng,nLat,nT), & WTP(NLOS,NL) , & WTP_in(N_IN) , & Ztmp(nLng,nLat,nT) character cStr*150 include 'MAPANGLES.H' C print *, 'In MkVMaptdn_inm',nLng,nLat,nT,nLOS,NL,N_IN,NS C if(NS.lt.0) 'This is an MHD source removal' nLngLat = nLng*nLat nLngLatnT = nLng*nLat*nT nLOSnL = NLOS*NL Bad = BadR4() nLat1 = nLat - 1 daLat2 = 90./nLat do I=1,nLng do J=1,nLat do k=1,nT VMA (I,J,K) = Bad AWTMAP (I,J,K) = Bad FIXM (I,J,K) = 0.0 WTSM (I,J,K) = 0.0 VMEANF (I,J,K) = 0.0 VFIXM2 (I,J,K) = 0.0 ANMAP (I,J,K) = Bad end do end do end do do I=1,NLOS do J=1,NL WTP(I,J) = 0.0 end do end do do I=1,N_IN WTP_in(I) = 0.0 end do C call ArrR4Bad (nLngLatnT,VMA) C call ArrR4Bad(nLngLatnT,AWTMAP) C call ArrR4Zero(nLngLatnT,FIXM) C call ArrR4Bad(nLngLatnT,WTSM) C call ArrR4Zero(nLngLatnT,VMEANF) C call ArrR4Zero(nLngLatnT,VFIXM2) C call ArrR4Bad(nLngLatnT,ANMAP) C call ArrR4Zero(nLOSnL,WTP) NSS = 0 do I=1,NL+N_IN if (NBADSV(I) .gt. 0) then if (abs(VSSIG(I)) .gt. SIGB) then ! source above limit NSS = NSS+1 if (NS .eq. 1) NBADSV(I) = -1 ! Mark source as bad kinematic if (NS .eq. -1) NBADSV(I) = -2 ! Mark source as bad MHD end if end if if (NBADSV(I) .gt. 0) then ! Use only valid sources LNlast = 0 LTlast = 0 ITlast = 0 if(I.le.NL) then sinE = sind(abs(XEV(I))) do J=1,NLOS LN = LON(J,I) ! map longitude LT = LAT(J,I) ! map latitude IT = ITIM(J,I) ! integer time c VSN = sinE/RP(J,I) WTSj = VWTij(J,I) Fws = WTSj*FIX(I) if(VMAP(LN,LT,IT) .ge. VLIM) Fws = WTSj C if(FIXM(LN,LT,IT).eq.Bad) FIXM(LN,LT,IT) = 0. FIXM(LN,LT,IT) = FIXM(LN,LT,IT) + Fws WTP(J,I) = WTSj if(WTSM(LN,LT,IT).eq.Bad) WTSM(LN,LT,IT) = 0. WTSM(LN,LT,IT) = WTSM(LN,LT,IT)+WTP(J,I) if(LN.ne.LNlast.or.LT.ne.LTlast.or.IT.ne.ITlast) then if(ANMAP(LN,LT,IT).eq.Bad) ANMAP(LN,LT,IT) = 0. C ANMAP(LN,LT,IT) = ANMAP(LN,LT,IT)+WTSj ! Tried 2016 ANMAP(LN,LT,IT) = ANMAP(LN,LT,IT)+1. ! The LOS crossings, even partial crossings, must be greater than 1.0 ! All WTSj need not be treated equally, however C if(LN.eq.62.and.LT.eq.11.and.IT.eq.46) ! Works on CR 2153 Nagoya Gen format C & write (*,'(A,4F10.3)') 'In Mkvmaptdn0n_inm 01 62, 11, 46',ANMAP(62,11,46),WTSM(62,11,46),FIXM(62,11,46),VMAP(62,11,46) end if LNlast = LN LTlast = LT ITlast = IT end do else ! Now do the in-situ values. LN = LON_in(I-NL) ! map longitude LT = LAT_in(I-NL) ! map latitude IT = ITIM_in(I-NL) ! integer time c VSN = sinE/RP(J,I) WTSj = VWTij_in(I-NL) Fws = WTSj*FIX(I) if (VMAP(LN,LT,IT) .ge. VLIM) Fws = WTSj C if(FIXM(LN,LT,IT).eq.Bad) FIXM(LN,LT,IT) = 0. FIXM(LN,LT,IT) = FIXM(LN,LT,IT) + Fws WTP_in(I-NL) = WTSj if(WTSM(LN,LT,IT).eq.Bad) WTSM(LN,LT,IT) = 0. WTSM(LN,LT,IT) = WTSM(LN,LT,IT)+WTP_in(I-NL) if(LN.ne.LNlast.or.LT.ne.LTlast.or.IT.ne.ITlast) then if(ANMAP(LN,LT,IT).eq.Bad) ANMAP(LN,LT,IT) = 0. C ANMAP(LN,LT,IT) = ANMAP(LN,LT,IT)+WTSj ! Tried 2016 ANMAP(LN,LT,IT) = ANMAP(LN,LT,IT)+1. ! The LOS crossings, even partial crossings, must be greater than 1.0 ! All WTSj need not be treated equally, however C if(LN.eq.62.and.LT.eq.11.and.IT.eq.46) ! Works on CR 2153 Nagoya Gen format C & write (*,'(A,4F10.3)') 'In Mkvmaptdn0n_inm 02 62, 11, 46',ANMAP(62,11,46),WTSM(62,11,46),FIXM(62,11,46),VMAP(62,11,46) end if C if((I-NL).ge.610.and.(I-NL).le.630) write(*,'(A,4I5,4F10.3)') ! Works on CR 2153 Nagoya Gen format C & 'In Mkvmaptdn0n_inm',I-NL,LN,LT,IT,ANMAP(LN,LT,IT),WTSM(LN,LT,IT),FIXM(LN,LT,IT),VMAP(LN,LT,IT) LNlast = LN LTlast = LT ITlast = IT end if end if end do C call Timesmooth(nLng,nLat,nT,FIXM, 1,1.0*CONT/aNday,0.,Ztmp) ! prior to 2016 BVJ commented out call Timesmooth(nLng,nLat,nT,FIXM, 0,1.0*CONT/aNday,0.,Ztmp) ! Worked 2016 BVJ 10n CR 2153 C call Timesmooth(nLng,nLat,nT,FIXM, 0,0.1*CONT/aNday,0.,Ztmp) ! Tried 2016 BVJ (smooth valid values) new for 20n C call Timesmooth(nLng,nLat,nT,WTSM, 1,1.0*CONT/aNday,0.,Ztmp) ! prior to 2016 BVJ commented out call Timesmooth(nLng,nLat,nT,WTSM, 0,1.0*CONT/aNday,0.,Ztmp) ! Worked 2016 BVJ 10n CR 2153 C call Timesmooth(nLng,nLat,nT,WTSM, 0,0.1*CONT/aNday,0.,Ztmp) ! Tried 2016 BVJ (smooth valid values) new for 20n C call Timesmooth(nLng,nLat,nT,ANMAP,1,1.0*CONT/aNday,0.,Ztmp) ! prior to 2016 BVJ do I=1,nLng ! This put in on 2016 and following do J=1,nLat WTSMcheck = 0.0 do N = 1, nT if(ANMAP(I,J,N).eq.0.0) ANMAP(I,J,N) = Bad if(WTSM (I,J,N).ne.Bad) WTSMcheck = WTSMcheck + WTSM (I,J,N) end do do N = 1, nT if(WTSMcheck.ge.1.0.and.WTSM(I,J,N).lt.1.0.and.WTSM(I,J,N).ne.Bad) then ANMAP(I,J,N) = Bad FIXM (I,J,N) = 0. ! BVJ 2016. This should force timesmooth to only use the in situ on this lat. WTSM (I,J,N) = 0. ! BVJ 2016. This should force timesmooth to only use the in situ on this lat. end if C if(I.eq.62.and.J.eq.11) then ! Works on CR 2153 Nagoya Gen format C if(VMAP(62,11,N).ne.Bad.and.ANMAP(62,11,N).eq.Bad) then C write(*,'(A,I5,A,3F10.3)') C & 'In Mkvmaptdn0n_inm pre 62, 11',N,' Bad ',FIXM(62,11,N),WTSM(62,11,N),VMAP(62,11,N) C end if C if(VMAP(62,11,N).eq.Bad.and.ANMAP(62,11,N).ne.Bad) then C write(*,'(A,I5,3F10.3,A)') C & 'In Mkvmaptdn0n_inm pre 62, 11',N,ANMAP(62,11,N),FIXM(62,11,N),WTSM(62,11,N),' Bad ' C end if C if(VMAP(62,11,N).ne.Bad.and.ANMAP(62,11,N).ne.Bad) then C if(WTSMcheck.ne.Bad) write(*,'(A,2F10.3)') C & 'In Mkvmaptdn0n_inm pre 62, 11 WTSMcheck.ne.Bad',WTSMcheck,FIXM(62,11,N) C write(*,'(A,I5,4F10.3)') C & 'In Mkvmaptdn0n_inm pre 62, 11',N,ANMAP(62,11,N),FIXM(62,11,N),WTSM(62,11,N),VMAP(62,11,N) C end if C end if end do end do end do do N = 1, nT C do I=1,nLng C do J=1,nLat C if(ANMAP(I,J,N).eq.0.0) ANMAP(I,J,N) = badR4() C end do C end do call ArrR4getminmax(nLngLat,FIXM(1,1,N),amin,amax) if(amax.ne.Bad) then ! Worked 2016 BVJ 10n CR 2153 C call GridSphere2D(ALng,nLng,nLat,1,FIXM(1,1,N),CONR,3,0.0,0.0) ! prior to 2016 BVJ (Commented out) call GridSphere2D(ALng,nLng,nLat,1,FIXM(1,1,N),CONR,0,0.0,0.0) ! this should force everything close in lat. to blend C call GridSphere2D(ALng,nLng,nLat,1,FIXM(1,1,N),CONR/2.0,0,0.0,0.0) ! Tried 2016 for 20n end if call ArrR4getminmax(nLngLat,WTSM(1,1,N),amin,amax) if(amax.ne.Bad) then ! Worked 2016 BVJ 10n CR 2153 C call GridSphere2D(ALng,nLng,nLat,1,WTSM(1,1,N),CONR,4,0.0,0.0) ! prior to 2016 BVJ (Commented out) call GridSphere2D(ALng,nLng,nLat,1,WTSM(1,1,N),CONR,0,0.0,0.0) ! this should force everything close in lat. to blend C call GridSphere2D(ALng,nLng,nLat,1,WTSM(1,1,N),CONR/2.0,0,0.0,0.0) ! Tried 2016 for 20n end if call ArrR4getminmax(nLngLat,ANMAP(1,1,N),amin,amax) C print *, 'In mkvmaptdN N, amax =', N, amax if(amax.ne.Bad) then C call GridSphere2D(ALng,nLng,nLat,1,ANMAP(1,1,N),CONR,3,0.0,0.0) ! prior to 2016 BVJ call GridSphere2D(ALng,nLng,nLat,1,ANMAP(1,1,N),CONR,0,0.0,0.0) C C Now modify the latitude map count to approximate the cos latitude filter effect C do J=1,nLat alat = abs(XLdeg(J)) if(alat.eq.90.) alat = alat - dalat2 sqonebycos = SQRT(1./cosd(alat)) do I=1,nLng if(ANMAP(I,J,N).ne.Bad) then ANMAP(I,J,N) = ANMAP(I,J,N)*sqonebycos end if end do end do C end if end do C write (*,'(A,4F10.3)') 'In Mkvmaptdn0n_inm 2 62, 11, 46',ANMAP(62,11,46),WTSM(62,11,46),FIXM(62,11,46),VMAP(62,11,46) C call Timesmooth(nLng,nLat,nT,FIXM, 0,1.0*CONT/aNday,0.,Ztmp) ! prior to 2016 BVJ commented out C call Timesmooth(nLng,nLat,nT,FIXM, 1,1.0*CONT/aNday,0.,Ztmp) ! tried (smooth and fill all values) see gridsphere above BVJ 2016 C call Timesmooth(nLng,nLat,nT,WTSM, 0,1.0*CONT/aNday,0.,Ztmp) ! prior to 2016 BVJ commented out C call Timesmooth(nLng,nLat,nT,WTSM, 1,1.0*CONT/aNday,0.,Ztmp) ! tried (smooth and fill all values) see gridsphere above BVJ 2016 call Timesmooth(nLng,nLat,nT,ANMAP,0,1.0*CONT/aNday,0.,Ztmp) ! prior to 2016 BVJ C call Timesmooth(nLng,nLat,nT,ANMAP,0,1.0*CONT,0.,Ztmp) ! prior to 2016 correct, but try this. (bad - aie3) C call Timesmooth(nLng,nLat,nT,ANMAP,2,1.0*CONT/aNday,0.,Ztmp) ! tried (smooth only invalid values) (age3) C call Timesmooth(nLng,nLat,nT,ANMAP,1,1.0*CONT/aNday,0.,Ztmp) ! tried (smooth and fill all ANMAP values) see gridsphere above BVJ 2016 C do N=1,nT C call ArrR4getminmax(nLngLat,WTSM(1,1,N),amin,amax) C a70 = WTSM(30,8,N) C print *, 'Final map wt min, max, +70 deg.', N, amin, amax, a70 C call ArrR4getminmax(nLngLat,ANMAP(1,1,N),amin,amax) C a70 = ANMAP(30,8,N) C print *, 'Final map # min, max, +70 deg.', N, amin, amax, a70 C end do do LN=1,nLng do LT=1,nLat do IT=1,nT if (FIXM(LN,LT,IT) .gt. 0. .and. WTSM(LN,LT,IT) .gt. 0.) then VMEANF(LN,LT,IT) = FIXM(LN,LT,IT)/WTSM(LN,LT,IT) ! FIXM is WTSM times FIX VMA(LN,LT,IT) = VMEANF(LN,LT,IT)*VMAP(LN,LT,IT) else ANMAP(LN,LT,IT) = 0. end if if(ANMAP(LN,LT,IT).eq.Bad) ANMAP(LN,LT,IT) = 0.0 ! After 2016 BVJ C if(LN.eq.62.and.LT.eq.11) then ! Works on CR 2153 Nagoya Gen format C if(VMA(62,11,IT).ne.Bad) then C write(*,'(A,I5,4F10.3)') C & 'In Mkvmaptdn0n_inm 3 62, 11',IT,ANMAP(62,11,IT),FIXM(62,11,IT),WTSM(62,11,IT),VMA(62,11,IT) C else C write(*,'(A,I5,3F10.3,A)') C & 'In Mkvmaptdn0n_inm 3 62, 11',IT,ANMAP(62,11,IT),FIXM(62,11,IT),WTSM(62,11,IT),' Bad' C end if C end if end do end do end do C write (*,'(A,4F10.3)') 'In Mkvmaptdn0n_inm 4 62, 11, 46',FIXM(62,11,46),VMEANF(62,11,46),VMAP(62,11,46),VMA(62,11,46) C calculation of the Chi square distribution for each latitude and longitude do I=1,NL+N_IN if (NBADSV(I) .gt. 0) then ! Use only valid sources if(I.le.NL) then do J=1,NLOS LN = LON(J,I) ! map longitude LT = LAT(J,I) ! map latitude IT = ITIM(J,I) ! map time if (VMA(LN,LT,IT) .ne. Bad) VFIXM2(LN,LT,IT) = VFIXM2(LN,LT,IT) & +WTP(J,I)*((FIX(I)-VMEANF(LN,LT,IT))**2) end do else LN = LON_in(I-NL) ! map longitude LT = LAT_in(I-NL) ! map latitude IT = ITIM_in(I-NL) ! map time if (VMA(LN,LT,IT) .ne. Bad) VFIXM2(LN,LT,IT) = VFIXM2(LN,LT,IT) & +WTP_in(I-NL)*((FIX(I)-VMEANF(LN,LT,IT))**2) end if end if end do C calculation of the "total" Chi square distribution for all latitudes and longitudes CONVER = 0 NCONVER = 0 do I=1,nLng do J=1,nLat do K=1,nT if (ANMAP(I,J,K) .lt. ANLIMIT) VMA(I,J,K) = Bad ! data not deconvolved if (ANMAP(I,J,K).ne.Bad.and.WTSM(I,J,K).ne.0.0.and.WTSM(I,J,K).ne.Bad) then AWTMAP(I,J,K) = WTSM(I,J,K) end if if (VMA(I,J,K) .ne. Bad) then CONVER = CONVER +VFIXM2(I,J,K) NCONVER = NCONVER+1 end if VMAP(I,J,K) = VMA(I,J,K) end do end do end do SIGOLD = SIGMAP SIGMAP = CONVER/NCONVER write (cStr,'(A,F4.1,A,F5.2,A)') 'velocity sources above',SIGB,' sigma:',100.*NSS/NL,'%' if (NS .eq. 1) write (cStr(itrim(cStr)+1:),'(A,I4,A)') '#you just removed',NSS,' kinematic model velocity sources' if (NS .eq. -1) write (cStr(itrim(cStr)+1:),'(A,I4,A)') '#you just removed',NSS,' MHD model velocity sources' write (cStr(itrim(cStr)+1:),'(A,F9.6,A,F8.3,A)') '#V map convergence =',SIGMAP, & ', from last =',abs(1-SIGMAP/SIGOLD)*100,'%' call Say('MkVMap','I','Info',cStr) return end