C C+ C NAME: C Write3D_bbtm4d_3 C PURPOSE: C Write magnetic field output to bb3d* file C CATEGORY: C I/O C CALLING SEQUENCE: subroutine Write3D_bbtm4d_3(bMagexx,XCbeGG,iT3D,nT3D,nTmaxG,BFACTORRR,BFACTORRC,BFACTORTC,BFACTORNC,RADMS,RAD1,dRR, & TT3D,XC3D,VV3D,BRRR,BRCC,BTCC,BNCC,BRss,BRCs,BTCs,BNCs,nLLng,nLLat,nMap,nT,BR3D,BT3D,BN3D,BB3) C INPUTS: C bMagexx logical If .true. write out magnetic 3D files. C XCbeGG(2,nTim) real Somehow related to the tomography times C Only XCbeGG(1,1) and XCbeGG(2,1) are used C iT3D integer 1st time index in TT3D processed C TT3D(iT3D) must be a valid time. C (the last time index processed is determined C internally as the last valid time in TT3D) C nT3D integer size of time dimension in arrays C TT3D,XC3D,VV3D and BRss. C (essential to control access to array C XC3D because of the extra trailing dimension) C nTmaxG integer Used to dimension BB3 in time C BFACTORRR real A multiplicative factor used to provide a magnetic field enhancement over the CSSS C as suggested by T. Hoeksema (Fall AGU 2013) C BFACTORRC real A multiplicative factor used to provide a closed r field enhancement over the usual C C BFACTORTC real A multiplicative factor used to provide a closed t field enhancement over the usual C C BFACTORNC real A multiplicative factor used to provide a closed n field enhancement C RADMS real The Magnetic field source surface C RAD1 real The tomography source surface. This and the above surface can be different C C The following three arrays are set up by Write3D_info. C C TT3D(nT3D) real Carrington times C These most likely will be the times spaced at C regular intervals in UT converted to Carrington times C XC3D(nLng,nLat,nRad,nT3D,3) real Shift needed to go back to source surface C (in Carrington variable). Will be negative. C Only XC3D(*,*,*,iT3D:jT3D,*) is used C VV3D(nLng,nLat,nRad,nT3D) real Velocity array C Only VV3D(*,*,*,iT3D:jT3D) is used C C BRRR(nLng,nLat,jT3D-kT3D+1) real scratch array; source surface magnetic field C C BRCC(nLng,nLat,jT3D-kT3D+1) real scratch array; Closed radial source surface magnetic field C C BTCC(nLng,nLat,jT3D-kT3D+1) real scratch array; Closed tangential source surface magnetic field C C BNCC(nLng,nLat,jT3D-kT3D+1) real scratch array; Closed normal source surface magnetic field C C BR3D(nLng,nLat,nRad) real scratch array; radial component written to file C C BT3D(nLng,nLat,nRad) real scratch array; tangential component written to file C C BN3D(nLng,nLat,nRad) real scratch array; normal component written to file C OUTPUTS: C (to file) C C BRss(nLng,nLat,jT3D-kT3D+1) real output source surface magnetic field C C BRCs(nLng,nLat,jT3D-kT3D+1) real output closed radial source surface magnetic field C C BTCs(nLng,nLat,jT3D-kT3D+1) real output closed tangential source surface magnetic field C C BNCs(nLng,nLat,jT3D-kT3D+1) real output closed normal source surface magnetic field C C BB3 (nLng,nLat,nRad,nMAXTG*4,3) real Maps of magnetic field component values at all heights and times C INCLUDE: include 'sun.h' include 't3d_array_3.h' include 'b3d_param_3.h' C include 't3d_grid_fnc_M.h' C CALLS: C FLINT, T3D_write_bb, T3D_get, T3D_set, T3D_iget, T3D_iset, T3D_get_grid C ArrR4Step, ArrR4Constant, ArrR4Interpolate, BadR4, BField_Get, Say C ArrR4TimesConstant C PROCEDURE: C > TT3D may be only partially filled with valid times, e.g TT3D(kT3D:jT3D). C kT3D and jT3D are determined internally as the 1st and last valid time C (.ne. BadR4()) in TT3D. The input value iT3D must be greater/equal kT3D C and less/equal jT3D). Output files are produced only for times iT3D to jT3D. C kT3D usually is smaller than jT3D. Magnetic source surface data are collected C between time indices kT3D and jT3D. This ensures that heliospheric grid elements C at early times (just above iT3D) after traceback will still have a valid C time associated with them (above kT3D) at the source surface. C C > The radial dependence of Br follows from flux conservation: C Br = BB*(RR/r)^2 C The tangential component follows from geometrical considerations C by assuming that the magnetic field lies along the Parker spiral: C (Omega is angular velocity of Sun): C Bt/Br = r*cos(Lat)*Omega/V C Bt = BB*(Omega*RR*cos(Lat)/V)*(RR/r) C Output are 'normalized' magnetic field components (making the C values independent of radial distance): C r^2*Br = RR^2*BB C r *Bt = RR^2*BB *Omega*cos(Lat)/V C Expressing radial distances in units of r0 = 1 AU): C (r/r0)^2*Br = (RR/r0)^2*BB C (r/r0) *Bt = (RR/r0)^2*BB * r0*Omega*cos(Lat)/V C In terms of quantities defined in include file sun.h: C r0*Omega/V = Sun__AU*10^13*Sun__Omega*10^-6/(V[km/s]*10^5) = 100*Sun__AU*Sun__Omega/V[km/s] C C We are mostly interested in magnetic field at Earth, i.e. about r0=1 AU. C Insitu magnetic field are usually expressed in nT, so we use the same units here. C We write 'normalized magnetic fields' in units of nT to file: (RR/r0)^2*Br for the radial C and (RR/r0)*Bt for the tangential component. These quantities should be more or less C constant across the heliosphere, since the radial dropoff has been taken out. At 1 AU C (RR/r0 = 1) it should reproduce the insitu values. C C *TJD*'s unit troubles C X's output is in mT or 10^-2G. to bring to nT: C (15/215)^2 * (10^-6 T / mT) (10^9 nT/ T) C C MODIFICATION HISTORY: C SEP-1999, Paul Hick (UCSD/CASS) C ???-2002, Tamsen Dunn (UCSD/CASS; tdunn@ucsd.edu) Made time dep., fleshed out. C JUL-2002, Paul Hick (UCSD/CASS; pphick@ucsd.edu) C Valid magnetic fields are returned only in grid points where the C source surface Br and V are valid. Previously Br would be valid C if only the source surface Br was valid (even if V was bad). C AUG-2002, Tamsen Dunn (UCSD/CASS; tdunn@ucsd.edu) C Altered Bt component to depend entirely on V3D. If V3D is bad, use VSS. C This revealed bug below: C SEP-2002, Paul Hick (UCSD/CASS; pphick@ucsd.edu) C Fixed big error in use of XC3D (did not take into account that the time dimension C was different from the V3D array). Now three pieces of information are expected C that are calculated by Write3D_info: TT3D, times at which magnetic fields are C required; XC3D and V3D: shifts and velocity at times TT3D. XC3D and V3D now have C the same dimensions. C JUN-2003, Paul Hick (UCSD/CASS) C Added processing of multiple magnetic field types C OCT-2003, Paul Hick (UCSD/CASS; pphick@ucsd.edu) C Added extra arguments iT3D and jT3D C- real XCbeGG(*) ! (2,nTim) integer iT3D integer nT3D real TT3D(*) real XC3D(*) C real XC3D(nLLng,nLLat,nMap,nTmaxG*4,3) real VV3D(*) real BRss(*) real BRCs(*) real BTCs(*) real BNCs(*) real BRRR(nLLng,nLLat,4*nTmaxG) ! Scratch file real BRCC(nLLng,nLLat,4*nTmaxG) ! Scratch file real BTCC(nLLng,nLLat,4*nTmaxG) ! Scratch file real BNCC(nLLng,nLLat,4*nTmaxG) ! Scratch file real BR3D(*) real BT3D(*) real BN3D(*) character cSay*16 /'Write3D_bbtm4d_3'/ real eps /0.00002/ integer iBBFlag(4) integer Nss(3) ! Array dimensions at source surface real Pss(3) ! Location at source surface equivalence (Nss(1),nLng) equivalence (Nss(2),nLat) equivalence (Nss(3),nTim) real Scale(4) /0.0,1.0,0.0,1.0/ real PwrR(2) /2.0,1.0/ character cFile*256 logical BField_Get_3 logical bMagexx integer iB parameter (nMaxT = 600) real XCbegMAT(nMaxT) real XCendMAT(nMaxT) real BB3 (nLLng,nLLat,nMap+1,nTmaxG,3) ! Output volumetric magnetic fields C real BB3 (nLLng,nLLat,nMap,nT,3) ! Output volumetric magnetic fields include 't3d_grid_fnc_M.h' C include 't3d_grid_fnc.h' C external fncrrindx C external fncrrhght external rrindx external rrhght common RADMSS,RAD11,dRR1 !------- ! loc3D converts from 3D location (iLng,iLat,iRad,iTim,iPos) to linear ! offset relative to start of array. loc3D(iLng,iLat,iRad,iTim,iPos) = ((((iPos-1)*nTmaxG*4+iTim-1)*nRad+iRad-1)*nLat+iLat-1)*nLng+iLng print*,' ' write(*,'(A,3F10.5,I4)') 'Inside write3d_bbtm4D_M RADMS,RAD1,dRR,nTim',RADMS,RAD1,dRR,nTim RADMSS = RADMS RAD11 = RAD1 dRR1 = dRR print *, 'Test of internal functions RRhght(I) and RRindx(R)' do I=1,3 RRvalue = RRhght(I) RRin = RRindx(RRvalue+0.01) print *, 'I, RRvalue, RRindx', I, RRvalue, RRvalue+0.01, RRin end do print *, ' ' Bad = BadR4() VT = 100.0*Sun__Omega*Sun__AU if (TT3D(iT3D) .eq. Bad) then ! Make sure TT3D(iT3D) is valid call Say(cSay,'W','time','at index iT3D is bad. No magnetic field calculated') return end if kT3D = iT3D-1 ! Search for first good time preceding iT3D do while (kT3D .ge. 1 .and. TT3D(kT3D) .ne. Bad) kT3D = kT3D-1 end do kT3D = kT3D+1 ! Index of earliest good time, or 1 jT3D = iT3D+1 ! Search for last good time after iT3D do while (jT3D .le. nT3D .and. TT3D(jT3D) .ne. Bad) jT3D = jT3D+1 end do jT3D = jT3D-1 ! Index of last good time, or nT3D print *, cSay, kT3D,iT3D,jT3D,nT3D !------- ! Setting the iteration counter to one more than the total number of ! iterations will suppress adding the iteration number to the output file name. call T3D_iget(T3D__ITERATIONS, 0, L ) call T3D_iset(T3D__ITERATION , 0, L+1) Scale(2) = Scale(2)*BFACTORR Scale(4) = Scale(4)*BFACTORR call T3D_set(T3D__SCALE,4,Scale) ! No scale change, but set anyways call T3D_set(T3D__R_PWR,2,PwrR ) ! Radial normalization powers !------- ! Update T3D structure. Storing nT3D is essential because BField_Get will use this ! to determine the dimension of the TT3D array. call T3D_iset(T3D__NTIM , 0, jT3D-kT3D+1) ! Must stay before T3D_get_grid call !! !------- ! Extract grid information. Note that nTim = jT3D. TT, dTT and dTTi (not used) are zero. call T3D_get_grid_3(TT,dTT,RR,dRR,nLng,nLat,nRad,nTim, dTTi,nLng1,nLat1) if (nTim .gt. nMaxT) call Say(cSay,'E','nMaxT','parameter too small') call ArrR4Constant(nTim,XCbeGG(1),XCbegMAT) call ArrR4Constant(nTim,XCbeGG(2),XCendMAT) iXTimf = (TT - int(TT))*24.0 ! This is usually 3 UT for the IPS STEL tomography FIXBRCs = ((1.0/RADMS)**(2.0-1.34)) ! Fixes the components to comply with Parker falloff (BVJ 2/2/2016) FIXBTCs = ((1.0/RADMS)**(1.0-1.34)) ! Needs the height of the magnetic surface here and not RR FIXBNCs = -1.0 FIXBRRR = (RADMS/RAD1)**(2.0) FIXBRCC = (RADMS/RAD1)**(2.0) FIXBTCC = (RADMS/RAD1) FIXBNCC = (RADMS/RAD1)**(1.34) do iB=1,B3D__N ! there are four types of magnetic field data possible so far do iBB=1,4 iBBFlag(iBB) = 0 end do do iBB=1,4 ! there are four possible component maps for each magnetic field !------- ! The magnetic field is returned in nano Tesla (nT). C if (BField_Get(iB,TT3D(kT3D),XCbegMAT,XCendMAT,BRss)) then if (iBB .eq. 1) then if (BField_Get_3(iB,iBB,TT3D(kT3D),XCbegMAT,XCendMAT,BRss)) then iBBFlag(IBB) = 1 C print *, 'Out of BField_Get_3 iBB = ', IBB go to 22 ! rr field end if else if (iBB .eq. 2) then if (BField_Get_3(iB,iBB,TT3D(kT3D),XCbegMAT,XCendMAT,BRCs)) then iBBFlag(IBB) = 1 C print *, 'Out of BField_Get_3 iBB = ', IBB go to 22 ! rc field end if else if (iBB .eq. 3) then if (BField_Get_3(iB,iBB,TT3D(kT3D),XCbegMAT,XCendMAT,BTCs)) then iBBFlag(IBB) = 1 C print *, 'Out of BField_Get_3 iBB = ', IBB go to 22 ! rt field end if else if (iBB .eq. 4) then if (BField_Get_3(iB,iBB,TT3D(kT3D),XCbegMAT,XCendMAT,BNCs)) then iBBFlag(IBB) = 1 C print *, 'Out of BField_Get_3 iBB = ', IBB go to 22 ! rn field end if end if if((iBBFlag(1)+iBBFlag(2)+iBBFlag(3)+iBBFlag(4)).eq.0) go to 23 ! no rr, rc, rt, or rn field data 22 continue print *, 'Out of BField_Get_3, into 22 continue = ', iB, iBB, iBBflag !------- ! Distance normalization: multiply radial component BR3D by RR*RR. ! The resulting r^2*Br should be in the range 1-10 nT. C call ArrR4TimesConstant(-nLng*nLat*nTim,BRss,RR*RR*BFACTOR,BRss) if(iBB.eq.1) call ArrR4TimesConstant(-nLng*nLat*nTim,BRss,RR*RR*BFACTORRR,BRss) if(iBB.eq.2) call ArrR4TimesConstant(-nLng*nLat*nTim,BRCs,RR*RR*FIXBRCs*BFACTORRC,BRCs) !((RR/1.3)**(2.0-1.34)) if(iBB.eq.3) call ArrR4TimesConstant(-nLng*nLat*nTim,BTCs,RR*FIXBTCs*BFACTORTC,BTCs) if(iBB.eq.4) call ArrR4TimesConstant(-nLng*nLat*nTim,BNCs,RR*FIXBNCs*BFACTORNC,BNCs) if(iBB.eq.1) call ArrR4TimesConstant(-nLng*nLat*nTim,BRRR,Bad,BRRR) if(iBB.eq.2) call ArrR4TimesConstant(-nLng*nLat*nTim,BRCC,Bad,BRCC) if(iBB.eq.3) call ArrR4TimesConstant(-nLng*nLat*nTim,BTCC,Bad,BTCC) if(iBB.eq.4) call ArrR4TimesConstant(-nLng*nLat*nTim,BNCC,Bad,BNCC) end do iTTTTi = 0 iTTTTim = 0 NNNL = 0 do iTim=iT3D,jT3D iTTTTi = iTTTTi +1 TTi = TT3D(iTim) print *, iTim, iTTTTi XCbeg = XCbegMAT(iTim-kT3D+1) ! Needed for(nTmaxG-iT3D) statement fnc XCvar XCend = XCendMAT(iTim-kT3D+1) XCrange = XCend-XCbeg XCb = TTi-0.5 ! Define region of interest (ROI) XCe = TTi+0.5 if (XCb .lt. XCbeg) then ! Keep ROI inside [XCbeg,XCend] XCb = XCbeg XCe = XCb+1.0 end if if (XCe .gt. XCend) then XCe = XCend XCb = XCe-1.0 end if call T3D_set (T3D__TIME , 0, TTi ) ! Time call T3D_iset(T3D__ITIME , 0, iTim-kT3D+1) ! Time index call T3D_set (T3D__XCMAT , 0, XCbeg) ! Begin of matrix call T3D_set (T3D__XCMAT+1, 0, XCend) ! End of matrix call T3D_set (T3D__XCROI , 0, XCb ) ! Begin ROI call T3D_set (T3D__XCROI+1, 0, XCe ) ! End ROI do iBB=1,4 if(iBB.eq.1) then ! First through at any time, make all bad to start do iRad=1,nRad do iLat=1,nLat do iLng=1,nLng iPos = loc3D(iLng,iLat,iRad,1,1) BR3D(iPos) = Bad BT3D(iPos) = Bad BN3D(iPos) = Bad end do end do end do end if do iRad=1,nRad C print *, 'iRad, RRhght =', iRad, RRhght(iRad) do iLat=1,nLat XL = cosd(XLdeg(iLat)) do iLng=1,nLng iPos = loc3D(iLng,iLat,iRad,1,1) !------- ! Trace heliospheric location back to the source surface. ! There are no bad values in XC3D. Statement functions can be used ! for longitude and latitude but NOT for the time. Pss(1) = XCindx( XCvar(iLng) + XC3D( loc3D(iLng,iLat,iRad,iTim,1) ) ) !------- ! We assume that the latitude shift is stored in units of grid spacings. ! This may be wrong. Currently the latitude shifts are zero so it doesn't ! matter. Pss(2) = float(iLat) + XC3D( loc3D(iLng,iLat,iRad,iTim,2) ) !------- ! To get the time index at the source surface we use the ArrR4Interpolate ! function. DO NOT USE the statement functions TTindx and TTtime here! ! Pss(3) will be bad if the traced-back time is outside the range of TT3D. ! Setting it to -1.0 forces FLINT to return a bad BR value. ! The shift needs to be in terms of Carrington variable since write3d_infotd3dMM.f ! outputs shifts in terms of time - BVJ Dec 5 2015. Pss(3) = ArrR4Interpolate(nTim, TT3D(kT3D), TTi+(XC3D( loc3D(iLng,iLat,iRad,iTim,3) )/sngl(SUN__SYNODICP))) if (Pss(3) .eq. Bad) Pss(3) = -1.0 !------ ! BR is the normalized radial magnetic field at source surface location Pss ! (interpolate on source surface array BRss) C BR = FLINT(-3,Nss,BRss,Pss,eps) ! r^2*Br at source surface if(iBB.eq.1) then BR = FLINT(-3,Nss,BRss,Pss,eps) ! r^2*Br at source surface end if if(iBB.eq.2) then BRC = FLINT(-3,Nss,BRCs,Pss,eps) end if if(iBB.eq.3) BTC = FLINT(-3,Nss,BTCs,Pss,eps) ! r^1*Btc at source surface if(iBB.eq.4) BNC = FLINT(-3,Nss,BNCs,Pss,eps) ! r^1*Bnc at source surface C if (iRad .eq. 10 .and. iLat .eq. nLat/2 .and. 133 .le. iTim .and. iTim .le. 143 ) C & print *, iLng, iLat, iTim, Pss(1), Pss(2), Pss(3), C & XC3D(loc3D(iLng,iLat,iRad,iTim,1)), TTi,BR if(iBB .eq. 1) then if(BR .ne. Bad) then BR3D(iPos) = BR else BR3D(iPos) = Bad end if if(RAD1.eq.RRhght(iRad)) BRRR(iLng,iLat,iTTTTi) = BR3D(iPos) end if C iLngLatRd2 = loc3D(nLng/2,nLat/2,nRad/2,1,1) C if(iPos.eq.iLngLatRd2.and.IBB.ge.1) Print *, 'Test 1 BR3D(iPos) =', IBB, iTim, BR3D(iPos) if(iBB.eq.2) then if (BR3D(iPos).eq.Bad) BR3D(iPos) = 0.0 if (BRC .ne. Bad .and. BR3D(iPos) .ne. Bad) then BR3D(iPos) = BR3D(iPos) + BRC ! r^2*Br else BR3D(iPos) = Bad end if if (BRC .ne. Bad) then if(RAD1.eq.RRhght(iRad)) BRCC(iLng,iLat,iTTTTi) = BRC ! r^2*Br else if(RAD1.eq.RRhght(iRad)) BRCC(iLng,iLat,iTTTTi) = Bad end if end if C if(iPos.eq.iLngLatRd2.and.IBB.ge.2) Print *, 'Test 2 BR3D(iPos) =', iBB, iTim, BR3D(iPos), BRC if(iBB.eq.3) then if (BTC .ne. Bad) then BT3D(iPos) = BTC ! r*Btc else BT3D(iPos) = Bad end if if(RAD1.eq.RRhght(iRad)) BTCC(iLng,iLat,iTTTTi) = BT3D(iPos) end if C if(iPos.eq.iLngLatRd2.and.IBB.ge.3) Print *, 'Test 1 BT3D(iPos) =', iBB, iTim, BT3D(iPos) if(iBB .eq. 4) VV = VV3D(loc3D(iLng,iLat,iRad,iTim,1)) ! Velocity in heliosphere if (iBB .eq. 4) then if(BT3D(iPos) .eq. Bad .and. BR3D(iPos).ne. Bad) BT3D(iPos) = 0.0 if(BT3D(iPos) .ne. Bad .and. BR3D(iPos).eq. Bad) BR3D(iPos) = 0.0 if(BR3D(iPos) .ne. Bad .and. BT3D(iPos) .ne. Bad .and. VV .ne. Bad) then BT3D(iPos) = -BR3D(iPos)*(VT*XL/VV) + BT3D(iPos) ! r*Bt + r*Btc else BT3D(iPos) = Bad end if if(BNC .ne. Bad) then BN3D(iPos) = BNC ! r*Bn else BN3D(iPos) = Bad end if if(RAD1.eq.RRhght(iRad)) BNCC(iLng,iLat,iTTTTi) = BN3D(iPos) aTTi = (TT+float(iTim-1)*dTT)*24.0 iaTTi = int(aTTi/24.0)*24 iXTim = nint(aTTi - float(iaTTi)) if(iXTim.eq.iXTimf) then ! write out every iXTimf UT (usually 3 UT) if(iLng.eq.1.and.iLat.eq.1.and.iRad.eq.1) iTTTTim = iTTTTim +1 BB3(iLng,iLat,iRad,iTTTTim,1) = BR3D(iPos) BB3(iLng,iLat,iRad,iTTTTim,2) = BT3D(iPos) BB3(iLng,iLat,iRad,iTTTTim,3) = BN3D(iPos) end if end if C iLngLatRd2 = loc3D(nLng/2,nLat/2,nRad/2,1,1) C if(iPos.eq.iLngLatRd2.and.iTim.lt.89) print *, iTim, iBB, BR3D(iPos), BT3D(iPos), BN3D(iPos) C if(iPos.eq.iLngLatRd2.and.iBB.eq.4.and.iTim.lt.89) print *, ' ' end do end do end do if(bMagexx.and.iBB.eq.4) call T3D_write_bb_3(B3D__PREFIX(iB),t3d, BR3D, BT3D, BN3D, 0, ' ', cFile) iprint = 0 ! iprint = 1 prints out the following NNN=iTTTTim if(NNN.gt.NNNL.and.iprint.eq.1) then do i=1,nLng do j=1,nLat do K=1,nRad if(I.eq.5.and.J.eq.5) write(*,'(A,4I4,3E11.3))'), & 'I,J,K,NNN BB31 BB32 BB33',I,J,K,NNN,BB3(I,J,K,NNN,1),BB3(I,J,K,NNN,2),BB3(I,J,K,NNN,3) if(I.eq.5.and.J.eq.5.and.K.eq.nRad) print *, 'Inside write3d_bbtm4d_3.f' end do end do end do NNNL = NNN end if C call T3D_write_bb(B3D__PREFIX(iB),t3d, BR3D, BT3D, 0, ' ', cFile) end do end do call ArrR4TimesConstant(-nLng*nLat*nTim,BRRR,FIXBRRR,BRss) call ArrR4TimesConstant(-nLng*nLat*nTim,BRCC,FIXBRCC,BRCs) call ArrR4TimesConstant(-nLng*nLat*nTim,BTCC,FIXBTCC,BTCs) call ArrR4TimesConstant(-nLng*nLat*nTim,BNCC,FIXBNCC,BNCs) 23 continue end do return end ********************************************************************************************** function rrindx(R) common RADMSS,RAD11,dRR1 if(R.lt.RAD11) then dRRR = RAD11 - RADMSS rrindx = 1+(R-RADMSS)/dRRR else rrindx = 2+(R-RAD11)/dRR1 end if return end ********************************************************************************************** function rrhght(I) common RADMSS,RAD11,dRR1 if(I.eq.1) then rrhght = RADMSS else rrhght = RAD11+(I-2)*dRR1 end if return end