fixed rhot print in central ray output (fort.4)

src/gray_params.f90

@@ -232,8 +232,8 @@ contains
     read(u,*) prfparam%iprof, prfparam%irho  ! irho=0,1,2 -> num profiles vs rhot,rhop,psin
     read(u,*) prfparam%filenm
 !     psbnd     value of psi ( > 1 ) of density boundary
-    read(u,*) prfparam%psnbnd !, prfparam%sspld
-    prfparam%sspld=0.001_wp_
+    read(u,*) prfparam%psnbnd, prfparam%sspld
+!    prfparam%sspld=0.001_wp_
 !     iscal    ne Te scaling 0: nustar=const,  1: n_greenw=const; 2 no rescaling
 !     factT factn   factor for Te&ne scaling 
     read(u,*) prfparam%factte, prfparam%factne, prfparam%iscal

src/graycore.f90

@@ -1959,7 +1959,7 @@ bb: do
     if(jk.eq.1) then
       phideg=atan2(yym,xxm)/degree
       if(psinv>=zero .and. psinv<=one) then
-        rhot=frhotor(psinv)
+        rhot=frhotor(sqrt(psinv))
       else
         rhot=1.0_wp_
       end if