A Discrete Model Atmosphere, UV-updated

I have updated my Discrete Model Atmosphere so that it now includes UV-forcing of the upper layers giving rise to a thermosphere. I have also included a "troposphere" where the heat absorption is uniform leading to a constant lapse rate in that region. I stress that these kind of toy models may very well become superfluous after a more complete simulation using the Navier-Stokes equations. Maybe Claes has some update on this. Anyway, regardless of its usefulness it gives rise to some questions of pure academic interest. Here is what it looks like now:

The thing I wanted to point out is the annoying jump discontinuity at the surface. Numerical experiments suggest that this jump can not be made smaller than F/(2k) regardless of the meshsize and other factors. I would very much like in input from some clever mathematician about the significance of this and how it could perhaps be circumvented in a more developed model. As a side remark, I think that Miskolczi discusses the topic of jump discontinuities at the surface in his paper. The problem is that I don't understand his paper (nor find it on the internet anymore). Input is very welcome.

Updated script:

import numpy as np
import matplotlib.pyplot as plt


interval = 6
trop = 1        ## height of the "troposphere"
meshsize = 0.1

N = int (interval/meshsize)

weight = np.zeros(N)

## Please note that the "weight" is not the actual weight but a positive
## function taking values between 0 and 1 which increases monotonically on the
## actual weight (mass), meant to quantify the "heat-absorption"

weight[0] = 1   ## The surface is given complete heat absoption


for idx in range(1,N):

    if idx*meshsize < trop:
        weight[idx] = 1*meshsize    ## Troposheric weight put to 1 (times meshsize)
    else:
        weight[idx] =  np.exp(-(idx*meshsize - trop))*meshsize
   


A = np.zeros([N,N])
   
for idx1 in range(N):
    for idx2 in range(N):
        if idx1 == idx2:
            if idx1 == 0:
                A[idx1,idx2] = -1
            else:
                A[idx1,idx2] = -2
           
        else:
            A[idx1,idx2] = weight[idx2]

            if idx2>idx1:

                for idx3 in range(idx1+1, idx2):
                    A[idx1,idx2] = A[idx1,idx2]*(1-weight[idx3])

            else:
                for idx3 in range(idx2+1, idx1):
                    A[idx1,idx2] = A[idx1,idx2]*(1-weight[idx3])    
           

k = 1           ## Conductivity
uv = 1          ## UV-factor
screen = 0.5    ## Screening of UV-light (not rigorous, just toy model)

F = np.zeros(N)

F[0] = 1.0    ## Solar radiation incident on surface


for idx in range(N):
    F[idx] = F[idx] + uv*screen**(N-idx)    ## adding UV-forcing

 

temp = np.linalg.solve(A,-F/k)  ## Forcing vector divided by the "conductivity"
                                ## or perhaps more accurately, the diffusion parameter


x = np.arange(0,interval,meshsize)



plt.plot(x, temp)
plt.ylabel('Temperature')
plt.xlabel('Position')

plt.show()