SUBROUTINE F(X, Y, YP, ERRMSG) * * This is an example on how to write SUBROUTINE F which is called by * SUBROUTINE DE. * * This subroutine sets up the ODEs for the quasi-static solar coronal loop * models. * IMPLICIT DOUBLE PRECISION (A-H, O-Z) REAL*8 KAPPA0 CHARACTER ERRMSG*78, INMSG*78 * DIMENSION Y(3), YP(3) * PARAMETER (NPTS=100) COMMON /CURRENT/ NZCURR COMMON /VARS/ FCOND(NPTS), T(NPTS), ELECN(NPTS), S(NPTS), Z(NPTS), 1 AREA(NPTS), GS(NPTS), E_INPUT(NPTS), RADLOSS(NPTS), NZONES * * Define constants. * KAPPA0 = 1.D-6 BK = 1.380662D-16 AMU = 1.6605654D-24 HMASS = 1.008D0*AMU * * ERRMSG on input contains information concerning the subroutine that called F. * INMSG = ERRMSG ERRMSG = ' ' * * Check for fatal errors: T (Y(1)) and ELECN (Y(3)) should not be less than 0. * IF (Y(1) .LE. 0.D0) ERRMSG = 1 'Temperature is less than or equal to 0!' IF (Y(3) .LE. 0.D0) ERRMSG = 1 'Electron density is less than or equal to 0!' * DO 1 I = 1, 3 1 YP(I) = 0.0D0 * IF (ERRMSG(1:1) .NE. ' ') THEN * WRITE(6, 5) X, Y(1), Y(2), Y(3), INMSG * 5 FORMAT('Error in F: s = ', 1PE11.4,', T(e) = ',E11.4, * 1 ', F(cond) = ',E11.4,',',/12X,'N(e) = ',E11.4,/'Called ', * 2 'from ',A78) RETURN ENDIF * * Get the radiative losses. * RLOSS = RADLOOKUP(Y(1)) * * Get energy input. * EINP = ENERGY_INP(Y(1), 0) * * The temperature gradient. * YP(1) = -Y(2) / KAPPA0 / Y(1)**2.5 / AREA(NZCURR) * * The conductive flux gradient. * YP(2) = AREA(NZCURR) * (EINP - RLOSS*Y(3)**2) * * The electron density gradient. * YP(3) = (Y(3)/Y(1)) * ((HMASS*GS(NZCURR)/2.D0/BK) - YP(1)) * RETURN END