subroutine radiate !JS Okt. 2014 -- Translated old radiate.f routine to new .f90 !VH-1 code. ! ! calculate the energy loss in the zone due to radiative cooling ! use global use sweepsize use zone implicit none ! integer :: i, j, n, k, index, neloss, jtp, jbt, itp, ibt real, dimension(maxsweep) :: tmpold, tmpnew, tbase, & & scrch, dtemp, dis, drhox, drhoy, anrm, & & elimit, qloss real, dimension(imax,jmax,kmax) :: dei real :: const1, const2, hdt, lambda0, lambda1, tmp0, elim real, parameter :: xlamfac=1.0 !JS-NOTE: xlamfac indat Input in old code, may be changed ! Lambda(T), for log T = 10^4 to 10^8, in steps of 0.01 in log T: real, dimension(0:400) :: lamdat data lamdat/ & ! Raymond cooling curve & 2.92836E-24, 3.43695E-24, 4.04705E-24, 4.77552E-24, 5.64662E-24, & & 6.68847E-24, 7.93165E-24, 9.40783E-24, 1.11466E-23, 1.31700E-23, & & 1.54870E-23, 1.80856E-23, 2.09255E-23, 2.39336E-23, 2.70041E-23, & & 3.00054E-23, 3.27936E-23, 3.52359E-23, 3.72282E-23, 3.87122E-23, & & 3.95926E-23, 3.98181E-23, 3.97595E-23, 3.94672E-23, 3.89973E-23, & & 3.84209E-23, 3.77888E-23, 3.71552E-23, 3.65688E-23, 3.60691E-23, & & 3.56886E-23, 3.54524E-23, 3.53792E-23, 3.54833E-23, 3.57763E-23, & & 3.62654E-23, 3.69593E-23, 3.78658E-23, 3.89938E-23, 4.03668E-23, & & 4.19640E-23, 4.37999E-23, 4.58700E-23, 4.81624E-23, 5.06556E-23, & & 5.33229E-23, 5.61447E-23, 5.90883E-23, 6.21137E-23, 6.51956E-23, & & 6.83216E-23, 7.14935E-23, 7.47379E-23, 7.80639E-23, 8.15185E-23, & & 8.51414E-23, 8.89808E-23, 9.30795E-23, 9.74830E-23, 1.02239E-22, & & 1.07396E-22, 1.13001E-22, 1.19088E-22, 1.25704E-22, 1.32888E-22, & & 1.40671E-22, 1.49081E-22, 1.58142E-22, 1.67879E-22, 1.78315E-22, & & 1.89474E-22, 2.01386E-22, 2.14083E-22, 2.27602E-22, 2.41972E-22, & & 2.57223E-22, 2.73376E-22, 2.90472E-22, 3.08617E-22, 3.27484E-22, & & 3.47171E-22, 3.67693E-22, 3.88872E-22, 4.10530E-22, 4.32448E-22, & & 4.54371E-22, 4.76017E-22, 4.97090E-22, 5.17311E-22, 5.36440E-22, & & 5.54275E-22, 5.70672E-22, 5.85537E-22, 5.98803E-22, 6.10413E-22, & & 6.20295E-22, 6.28375E-22, 6.34534E-22, 6.38667E-22, 6.40659E-22, & & 6.40466E-22, 6.37853E-22, 6.33152E-22, 6.26635E-22, 6.18700E-22, & & 6.09823E-22, 6.00532E-22, 5.91277E-22, 5.82580E-22, 5.74755E-22, & & 5.68150E-22, 5.62873E-22, 5.59016E-22, 5.56631E-22, 5.55652E-22, & & 5.55965E-22, 5.57384E-22, 5.59731E-22, 5.62861E-22, 5.66504E-22, & & 5.70487E-22, 5.74636E-22, 5.78707E-22, 5.82621E-22, 5.86247E-22, & & 5.89538E-22, 5.92521E-22, 5.95179E-22, 5.97622E-22, 5.99854E-22, & & 6.01960E-22, 6.03899E-22, 6.05605E-22, 6.06949E-22, 6.07754E-22, & & 6.07698E-22, 6.06500E-22, 6.03715E-22, 5.98990E-22, 5.91838E-22, & & 5.81914E-22, 5.68882E-22, 5.52572E-22, 5.32994E-22, 5.10377E-22, & & 4.85170E-22, 4.58018E-22, 4.29652E-22, 4.00890E-22, 3.72490E-22, & & 3.45100E-22, 3.19301E-22, 2.95282E-22, 2.73322E-22, 2.53556E-22, & & 2.35844E-22, 2.20166E-22, 2.06461E-22, 1.94458E-22, 1.84128E-22, & & 1.75194E-22, 1.67608E-22, 1.61123E-22, 1.55676E-22, 1.51054E-22, & & 1.47192E-22, 1.43933E-22, 1.41135E-22, 1.38739E-22, 1.36593E-22, & & 1.34645E-22, 1.32786E-22, 1.30921E-22, 1.29028E-22, 1.27051E-22, & & 1.24970E-22, 1.22767E-22, 1.20510E-22, 1.18233E-22, 1.15978E-22, & & 1.13785E-22, 1.11743E-22, 1.09884E-22, 1.08239E-22, 1.06818E-22, & & 1.05658E-22, 1.04753E-22, 1.04093E-22, 1.03660E-22, 1.03435E-22, & & 1.03384E-22, 1.03503E-22, 1.03754E-22, 1.04107E-22, 1.04533E-22, & & 1.05004E-22, 1.05494E-22, 1.05969E-22, 1.06423E-22, 1.06836E-22, & & 1.07193E-22, 1.07489E-22, 1.07722E-22, 1.07899E-22, 1.08016E-22, & & 1.08073E-22, 1.08090E-22, 1.08072E-22, 1.08012E-22, 1.07923E-22, & & 1.07798E-22, 1.07638E-22, 1.07445E-22, 1.07210E-22, 1.06932E-22, & & 1.06594E-22, 1.06204E-22, 1.05741E-22, 1.05201E-22, 1.04558E-22, & & 1.03771E-22, 1.02889E-22, 1.01811E-22, 1.00607E-22, 9.92107E-23, & & 9.76254E-23, 9.58223E-23, 9.37777E-23, 9.15551E-23, 8.91429E-23, & & 8.60895E-23, 8.33823E-23, 8.05534E-23, 7.76966E-23, 7.48259E-23, & & 7.19784E-23, 6.91894E-23, 6.64788E-23, 6.38798E-23, 6.13985E-23, & & 5.90465E-23, 5.67682E-23, 5.46860E-23, 5.27256E-23, 5.09142E-23, & & 4.92163E-23, 4.76640E-23, 4.62352E-23, 4.49228E-23, 4.37225E-23, & & 4.26227E-23, 4.16270E-23, 4.07269E-23, 3.99151E-23, 3.91863E-23, & & 3.85367E-23, 3.79509E-23, 3.74488E-23, 3.70152E-23, 3.66467E-23, & & 3.63123E-23, 3.60697E-23, 3.58855E-23, 3.57569E-23, 3.56813E-23, & & 3.56554E-23, 3.56762E-23, 3.57403E-23, 3.58367E-23, 3.59898E-23, & & 3.61226E-23, 3.63162E-23, 3.65379E-23, 3.67721E-23, 3.70341E-23, & & 3.72856E-23, 3.75251E-23, 3.77640E-23, 3.79964E-23, 3.82086E-23, & & 3.83945E-23, 3.85050E-23, 3.86249E-23, 3.86801E-23, 3.86797E-23, & & 3.86372E-23, 3.85180E-23, 3.83437E-23, 3.81109E-23, 3.77901E-23, & & 3.74187E-23, 3.69709E-23, 3.64503E-23, 3.58748E-23, 3.52477E-23, & & 3.45769E-23, 3.38525E-23, 3.31142E-23, 3.23473E-23, 3.14172E-23, & & 3.06235E-23, 2.98223E-23, 2.90281E-23, 2.82402E-23, 2.74777E-23, & & 2.67383E-23, 2.60250E-23, 2.53450E-23, 2.46977E-23, 2.40845E-23, & & 2.35069E-23, 2.29648E-23, 2.24517E-23, 2.19748E-23, 2.15378E-23, & & 2.10915E-23, 2.07138E-23, 2.03370E-23, 2.00227E-23, 1.97349E-23, & & 1.94718E-23, 1.90757E-23, 1.88575E-23, 1.86603E-23, 1.84824E-23, & & 1.83226E-23, 1.81798E-23, 1.80402E-23, 1.78860E-23, 1.77401E-23, & & 1.76579E-23, 1.75881E-23, 1.75299E-23, 1.74806E-23, 1.74332E-23, & & 1.74064E-23, 1.73855E-23, 1.73733E-23, 1.73689E-23, 1.73635E-23, & & 1.73747E-23, 1.73927E-23, 1.74170E-23, 1.74455E-23, 1.74820E-23, & & 1.75233E-23, 1.75676E-23, 1.76165E-23, 1.76701E-23, 1.77281E-23, & & 1.77811E-23, 1.78476E-23, 1.79175E-23, 1.79910E-23, 1.80674E-23, & & 1.81467E-23, 1.82287E-23, 1.83134E-23, 1.84002E-23, 1.84891E-23, & & 1.85801E-23, 1.86460E-23, 1.87406E-23, 1.88364E-23, 1.89334E-23, & & 1.90310E-23, 1.91001E-23, 1.92019E-23, 1.91526E-23, 1.92560E-23, & & 1.93593E-23, 1.94807E-23, 1.95842E-23, 1.96853E-23, 1.97857E-23, & & 1.98847E-23, 1.99802E-23, 2.00765E-23, 2.01713E-23, 2.02648E-23, & & 2.03565E-23, 2.04439E-23, 2.05294E-23, 2.06132E-23, 2.06938E-23, & & 2.07723E-23, 2.08488E-23, 2.09225E-23, 2.09925E-23, 2.10622E-23, & & 2.11274E-23, 2.11898E-23, 2.12490E-23, 2.13049E-23, 2.12422E-23, & & 2.12611E-23, 2.12861E-23, 2.13316E-23, 2.13757E-23, 2.14163E-23, & & 2.14537E-23 / ! !------------------------------------------------------------------- const2 = - 1.0*gamm/avgmass/boltzman const1 = avgmass/boltzman hdt = 0.5*dt ! do k = 1, kmax do j = 1, jmax ! ! dTemp = .5 * anrm*zro*(Lambda(oldtemp) + Lambda(netmin)) ! iterate a few times, replacing netmin with tmpnew + dTemp ! do i = 1, imax tmpold(i) = const1 * zpr(i,j,k) / zro(i,j,k) tmpold(i) = max(tmpold(i),tempw) anrm (i) = const2 * zro(i,j,k) * xlamfac enddo do i = 1, imax index = nint (100.* (log10(tmpold(i)) - 4.)) index = max (0, min (400,index)) qloss(i) = anrm(i) * lamdat(index) tbase(i) = tmpold(i) + hdt * qloss(i) tbase(i) = max(tbase(i),tempw) tmpnew(i) = tbase(i) enddo do neloss = 1, 3 do i = 1, imax index = nint (100.* (log10 (tmpnew(i)) - 4.)) index = max (0, min (400,index)) qloss(i) = anrm(i) * lamdat(index) tmpnew(i) = tbase(i) + hdt * qloss(i) tmpnew(i) = max(tmpnew(i),tempw) enddo enddo ! do i = 1, imax dei(i,j,k) = (tmpold(i) - tmpnew(i))*zro(i,j,k) enddo ! enddo enddo ! ! replace deltaE with minimum of neighboring deltaEs at cloud interface ! and subtract dT from the temperature given by the hydro ! do k = 1, kmax do j = 1, jmax ! ! Locate cloud interfaces => dis = 1, otherwise, dis = -1 ! jtp = min(j+1,jmax) jbt = max(j-1,1) do i = 1+1, imax-1 scrch(i) = min(zro(i+1,j,k),zro(i-1,j,k)) drhox(i) = zro(i+1,j,k) - zro(i-1,j,k) drhox(i) = sign(drhox(i),(zpr(i-1,j,k)/zro(i-1,j,k) & & - zpr(i+1,j,k)/zro(i+1,j,k)))/scrch(i)-2 scrch(i) = min(zro(i,jtp,k),zro(i,jbt,k)) drhoy(i) = zro(i,jtp,k) - zro(i,jbt,k) drhoy(i) = sign(drhoy(i),(zpr(i,jbt,k)/zro(i,jbt,k) & & - zpr(i,jtp,k)/zro(i,jtp,k)))/scrch(i)-2 dis (i) = max(drhox(i),drhoy(i)) enddo dis(1) = -1. dis(imax) = -1. ! elim=small do i = 1, imax itp = min(i+1,imax) ibt = max(i-1,1) scrch (i) = min(dei(ibt,j,k),dei(itp,j,k), & & dei(i,jtp,k),dei(i,jbt,k)) ! dei(i,j,k) = cvmgm(dei(i,j,k),scrch(i),dis(i)) ! if (dis(i).gt.0.0) dei(i,j,k) = scrch(i) !NOTE- cvmgm function redundant in new version !Using merge from Runacres old version instead dei(i,j,k) = merge(dei(i,j,k),scrch(i),dis(i)<0) tmpnew(i) = const1 * zpr(i,j,k) / zro(i,j,k) dtemp(i) = dei(i,j,k)/zro(i,j,k) tmpnew(i) = max((tmpnew(i)-dtemp(i)),tempw) zpr(i,j,k) = tmpnew(i) * zro(i,j,k) / const1 elimit(i) = max(dtemp(i),small) elimit(i) = tmpnew(i) / elimit(i) elim = min(elim,elimit(i)) enddo ! ! elim is the smallest T/dT ! enddo enddo ! return end subroutine radiate