C ============================================================================= C ROUTINE FOR GENERATING [H](z) PROFILE FOR MSIS THERMOSPHERE. C - diffusion coefficients updated [13MAY99 JBishop] C ----------------------------------------------------------------------------- C BASIC INFORMATION (ALTITUDES IN KM) | C | C MOD_MSIS . . FLAG TO SPECIFY WHETHER OR NOT MSIS [H](z) IS "CORRECTED" | C 0 (OR NEGATIVE): MSIS [H](z) NOT "CORRECTED" | C 1 (OR > 1): MSIS [H](z) "CORRECTED" | C | C D_EXO . . . USER-SPECIFIED PARAMETER: EXOBASE DENSITY | C IF VALUE < 0, THEN MSIS [H](z) EXOBASE VALUE USED | C | C FLUX . . . USER-SPECIFIED PARAMETER: VERTICAL ATOMIC HYDROGEN FLUX | C IF VALUE < 0, THEN NOMINAL MSIS VALUE IS ADOPTED | C -- 2.2D8 cm^(-2) s^(-1) | C | C D_MAX . . . USER-SPECIFIED PARAMETER: MESOSPHERIC PEAK DENSITY | C IF VALUE < 0, THEN NOMINAL MSIS VALUE IS ADOPTED | C -- 1.5D8 cm^(-3) | C | C DEQUIL . . . DIFFUSIVE EQUILIBRIUM [H](z) PROFILE BASED ON D_EXO AND | C MSIS TEMPERATURE PROFILE | C | C DFLUX . . . DIFFUSIVE FLOW [H](z) PROFILE BASED ON D_EXO, FLUX & D_MAX | C WITH MSIS [N2](z), [O2](z), [O](z) AND T(z) PROFILES | C | C DMOCK . . . MOCK CHAPMAN [H](z) PROFILE FOR MESOSPHERIC REGION | C | C ALT_MAX, ALT_JNT are not considered to be user parameters | C C ----------------------------------------------------------------------------- SUBROUTINE H_DENSITY IMPLICIT REAL*8 (A-H,O-Z) REAL*4 SATT, SATD, SCALN, ADJT, BASE, TOP DIMENSION THERMO_MSIS(7,61) DIMENSION DEQUIL(61), DFLUX(61), DMOCK(61) COMMON/PARA_EXOS/ IGEO,SATT,SATD,SCALN,ADJT COMMON/H_PROFILE/ D_EXO,FLUX,D_MAX,ALT_JNT,ALT_MAX,T_EXO,MOD_MSIS COMMON/FINE_GRID/ THERMO_MSIS,BASE,TOP,ITHERM C ----------------------------------------------------------------------------- C SETUP C ----------------------------------------------------------------------------- IKNT = ITHERM + 1 ADJT = 0.0 T_EXO = THERMO_MSIS(3,ITHERM+1) ALT_MAX = 85.0d0 ALT_JNT = 112.5d0 C ----------------------------------------------------------------------------- C ORIGINAL MSIS [H](z) PROFILE KEPT C ----------------------------------------------------------------------------- IF ( MOD_MSIS .LE. 0 ) THEN SCALN = 1.0 D_EXO = -1.0d0 FLUX = -1.0d0 D_MAX = -1.0d0 ELSE C ----------------------------------------------------------------------------- C DIFFUSIVE FLOW [H](z) PROFILE GENERATED C ----------------------------------------------------------------------------- IF (D_EXO.LE.0.0d0) D_EXO = THERMO_MSIS(4,ITHERM+1) IF (FLUX .LE.0.0d0) FLUX = 2.2D8 ! nominal MSIS value IF (D_MAX.LE.0.0d0) D_MAX = 1.5D8 ! nominal MSIS value SCALN = SNGL(D_EXO / THERMO_MSIS(4,ITHERM+1)) C get diffusive flux profile and overwrite THERMO_MSIS(4,*) CALL VERTICAL_FLUX(DEQUIL,DFLUX) DO II = 1,IKNT THERMO_MSIS(4,II) = DFLUX(II) END DO C ----------------------------------------------------------------------------- C MOCK-CHAPMAN PROFILE APPLIED (IF CALLED FOR) C ----------------------------------------------------------------------------- C 1) for Lyman beta & the higher-lying Lyman lines, the "black level" C lies above the mesospheric peak. hence this peak is not within C the RT altitude grid and adjusments to ALT_MAX are not justified. C C 2) for some lines (e.g., Lyman beta), the black level lies beneath C the "joint", calling for evaluation of the mock-Chapman segment. C for other lines (e.g., Lyman gamma), the black level lies above C the "joint" and the mock-Chapman profile is not used (i.e., the C D_MAX parameter has no direct impact). C ----------------------------------------------------------------------------- IF ( THERMO_MSIS(1,1) .LT. ALT_JNT ) THEN IAJNT = 0 ATEST = 100.0d0 DO II = 1,IKNT ADUM = ABS(THERMO_MSIS(1,II) - ALT_JNT) IF ( ADUM .LT. ATEST ) THEN ATEST = ADUM IAJNT = II END IF END DO ZJNT = THERMO_MSIS(1,IAJNT) DJNT = THERMO_MSIS(4,IAJNT) C _____ silly little iterative relaxation for effective Chapman scale height SCAL_EFF = (ZJNT - ALT_MAX) / LOG(D_MAX / DJNT) DO KK = 1,10 XPNT = (ZJNT - ALT_MAX) / SCAL_EFF CHAPMAN = D_MAX * EXP(1.0d0 - XPNT - EXP(-XPNT)) SCAL_EFF = SCAL_EFF - LOG(CHAPMAN/DJNT) END DO C _____ mock-Chapman profile DO II = 1,IKNT XPNT = (THERMO_MSIS(1,II) - ALT_MAX) / SCAL_EFF DMOCK(II) = D_MAX * EXP(1.0d0 - XPNT - EXP(-XPNT)) END DO C _____ impose mock-Chapman profile below the "joint" DO II = 1,IAJNT THERMO_MSIS(4,II) = DMOCK(II) END DO END IF ! end of imposing mock-Chapman MLT profile END IF ! end of adjusting (correcting) MSIS profile C ----------------------------------------------------------------------------- RETURN END C ============================================================================= C ============================================================================= C procedure to evaluate diffusive flow density profiles for atomic hydrogen C using MSIS-generated thermospheric profiles for temperature, [N2], [O2] & [O] SUBROUTINE VERTICAL_FLUX(DEQUIL,DFLUX) IMPLICIT REAL*8 (A-H,O-Z) REAL*4 BASE, TOP DIMENSION THERMO_MSIS(7,61) DIMENSION DEQUIL(61), DFLUX(61) DIMENSION SCAL_D(61), XPNT(61), D_mean(61) COMMON/H_PROFILE/ D_EXO,FLUX,D_MAX,ALT_JNT,ALT_MAX,T_EXO,MOD_MSIS COMMON/FINE_GRID/ THERMO_MSIS,BASE,TOP,ITHERM C ----------------------------------------------------------------------------- C SETUP C ----------------------------------------------------------------------------- IKNT = ITHERM + 1 C thermal diffusion coefficient (value from Hodges 1993) TDC = -0.20d0 DO II = 1,IKNT C pressure scale heights GRAV = 3.9898d20 / THERMO_MSIS(2,II)**2 SCAL_P = 1.3806d-16*THERMO_MSIS(3,II) / (1.6735d-24*GRAV) C temperature derivative IF (II.EQ.1) THEN dTdz = (THERMO_MSIS(3,2)-THERMO_MSIS(3,1)) & / (THERMO_MSIS(2,2)-THERMO_MSIS(2,1)) ELSE IF (II.EQ.IKNT) THEN dTdz = (THERMO_MSIS(3,IKNT)-THERMO_MSIS(3,IKNT-1)) & / (THERMO_MSIS(2,IKNT)-THERMO_MSIS(2,IKNT-1)) ELSE dTdz = (THERMO_MSIS(3,II+1)-THERMO_MSIS(3,II-1)) & / (THERMO_MSIS(2,II+1)-THERMO_MSIS(2,II-1)) END IF C density scale heights SCAL_D(II) = (1.0d0/SCAL_P) & + (1.0d0+TDC) * (dTdz/THERMO_MSIS(3,II)) END DO C ----------------------------------------------------------------------------- C approximate evaluation of scale height integral C ----------------------------------------------------------------------------- XPNT(IKNT) = 0.0d0 DO II = IKNT-1,1,-1 DZ = THERMO_MSIS(2,II+1) - THERMO_MSIS(2,II) SCAL_AVE = (SCAL_D(II+1) + SCAL_D(II)) / 2.0d0 XPNT(II) = XPNT(II+1) + DZ * SCAL_AVE END DO C ----------------------------------------------------------------------------- C DIFFUSION COEFFICIENTS C ----------------------------------------------------------------------------- DO II = 1,IKNT C simple diffusion equilibrium T_TERM = (T_EXO/THERMO_MSIS(3,II))**(1.0d0+TDC) DEQUIL(II) = D_EXO * T_TERM * EXP(XPNT(II)) C mean diffusion coefficient C D_ao . . . atomic oxygen, from Hodges (1993) C D_o2 . . . molecular oxygen, from Banks & Kockarts (1973) C D_n2 . . . molecular nitrogen, from Stallcop et al (1992) D_ao = 4.43d17*(THERMO_MSIS(3,II)**0.750d0)/THERMO_MSIS(7,II) D_o2 = 4.75d17*(THERMO_MSIS(3,II)**0.711d0)/THERMO_MSIS(5,II) D_n2 = 4.12d17*(THERMO_MSIS(3,II)**0.750d0)/THERMO_MSIS(6,II) D_inv = (1.0d0/D_ao) + (1.0d0/D_o2) + (1.0d0/D_n2) D_mean(II) = 1.0d0 / D_inv END DO C ----------------------------------------------------------------------------- C diffusion equation propogator: C ----------------------------------------------------------------------------- C - note that integration is from upper boundary downward C - Runge-Kutta 4th order (even though vertical grid is nonuniform) C - ignoring eddy mixing C ----------------------------------------------------------------------------- DFLUX(IKNT) = D_EXO DO II = IKNT-1,1,-1 C _____ preliminary DZ = THERMO_MSIS(2,II+1) - THERMO_MSIS(2,II) DZHALF = DZ / 2.0d0 SCAL_top = SCAL_D(II+1) FLUX_top = FLUX/D_mean(II+1) SCAL_bot = SCAL_D(II) FLUX_bot = FLUX/D_mean(II) SCAL_mid = (SCAL_top + SCAL_bot) / 2.0d0 FLUX_mid = (FLUX_top + FLUX_bot) / 2.0d0 C _____ step 1 dndz_1 = FLUX_top + DFLUX(II+1) * SCAL_top dmid_1 = DFLUX(II+1) + dndz_1 * DZHALF C _____ step 2 dndz_2 = FLUX_mid + dmid_1 * SCAL_mid dmid_2 = DFLUX(II+1) + dndz_2 * DZHALF C _____ step 3 dndz_3 = FLUX_mid + dmid_2 * SCAL_mid dmid_3 = DFLUX(II+1) + dndz_3 * DZ C _____ step 4 dndz_4 = FLUX_bot + dmid_3 * SCAL_bot C _____ final rksum = dndz_1 + 2.0d0*dndz_2 + 2.0d0*dndz_3 + dndz_4 DFLUX(II) = DFLUX(II+1) + DZ * rksum / 6.0d0 C simple euler propagation C dndz = FLUX_mid + DFLUX(II+1)*SCAL_mid C DFLUX(II) = DFLUX(II+1) + dndz * DZ END DO C ----------------------------------------------------------------------------- RETURN END C =============================================================================