有偿找人帮忙Fortran求解微分方程

123456789

求visual studio2019版，Fortran2022版求解这种微分方程应该怎么编写代码啊，完全小白一点不懂，想得到s在不同的z(z可以只取总的高度下的几个有代表性的高度就行）下随时间t（t要分的细一点，想要得到的m值多一点）变化的值，我在网上拼拼凑凑瞎写了一下代码好像不能直接识别d和德尔塔，果然错，提示有很多问题，我也不知道咋改，有搜到要用四阶龙格库塔法，但我完全看不懂，球球大神救命！    program Composition
    implicit none
    ! Variables declaration
    real :: ca ! Specific heat capacity of dry air J /kg*k
    real :: cv ! specific heat capacity of water vapor J /kg*k
    real :: cg ! Specific heat capacity of dry grains J /kg*k
    real :: cw ! Specific heat capacity of liquid water J /kg*k
    real :: Tem ! Air temperature before drying ℃
    real :: H ! Absolute humidity of air before drying g/kg
    real :: b! Dry grain capacity kg/m³
    real :: G ! Grain weight t
    real :: t ! drying time s
    real :: z ! Height of grain pile m
    real :: r ! heat of vaporization of water  J/kg
    real :: S0 ! Initial moisture content of grain  kg/kg
    real :: S ! moisture content of grain  kg/kg
    real :: Se ! Grain Balance Moisture kg/kg
    real :: K ! Parameters of the thin-layer drying equation
    real :: a ! Parameters of the thin-layer drying equation
    !assigning values
    ca=1029
    cv=1880
    cg=3858
    cw=4200
    Tem=35
    H=12.9
    b=600
    G=100
    r=12200
    S0=0.03158
    z=3.75
    !Humidity control equation
    dH=-(b/G)*(dS/dt)*dz
    !Temperature control equation
    dTem=(((cv-cw)*Tem+r)*dS)/((cg++cw*S)+(cw-cv)*dS+(G/b)*(ca+cv*H)*(dt/dz))
    !Moisture control equation
    (S-Se)/(S0-Se)=exp(-K*(t**a))
    !Thin-layer drying equation
    dS/dt=-k*(S-Se)
    k=0.242261*exp(-2530.2/(Tem+273.15)
    ！Initial boundary condition
    Tem(0,t)=Tem
    H(0,t)=H
    S(z,0)=S0
    write(*,*）"the drying time is: ", t
    ! output
    print *, "End-of-cereal moisture content=", S
     end program Composition

123456789

我知道我根本没有用到什么微分方程求解的，因为根本看不懂