diff -r -x .git SWMF/CON/Coupler/src/CON_couple_points.f90 SWMF_yesterday/CON/Coupler/src/CON_couple_points.f90
267c267
<        if(DoDebug) iUnitTest = i_proc() + 100
---
>        if(DoDebug) iUnitTest = i_proc() + 10
diff -r -x .git SWMF/GM/BATSRUS/srcInterface/GM_couple_pc.f90 SWMF_yesterday/GM/BATSRUS/srcInterface/GM_couple_pc.f90
28,29c28,29
<     use ModMain, ONLY: Dt
<     use ModPhysics, ONLY: No2Si_V, UnitT_
---
>     use ModMain,            ONLY: Dt
>     use ModPhysics,         ONLY: No2Si_V, UnitT_
45c45
<     use ModPIC, ONLY:DoRestartPicStatus
---
>     use ModPIC,      ONLY:DoRestartPicStatus
68,70c68,70
<     use ModAdvance, ONLY: UseMultiSpecies, nSpecies
<     use ModPhysics, ONLY: No2Si_V, UnitX_, PePerPtotal, rPlanetSi
<     use ModPIC, ONLY: XyzMinPic_DI, nRegionPiC, &
---
>     use ModAdvance,    ONLY: UseMultiSpecies, nSpecies
>     use ModPhysics,    ONLY: No2Si_V, UnitX_, PePerPtotal, rPlanetSi
>     use ModPIC,        ONLY: XyzMinPic_DI, nRegionPiC, &
73c73
<     use BATL_lib, ONLY: x_, y_, z_, nDim
---
>     use BATL_lib,      ONLY: x_, y_, z_, nDim
272,273c272,274
<     use BATL_lib, ONLY: nDim, Xyz_DGB, nBlock, Unused_B, nI, nJ, nK
<     use ModPIC, ONLY: DxyzPic_DI, LenPic_DI, &
---
>     use BATL_lib,     ONLY: nDim, Xyz_DGB, nBlock, Unused_B, &
>          nI, nJ, nK
>     use ModPIC,       ONLY: DxyzPic_DI, LenPic_DI, &
276,279c277,280
<     use ModPhysics, ONLY: No2Si_V, UnitX_, Si2No_V, iUnitCons_V
<     use ModMain, ONLY: UseB0, UseHyperbolicDivB
<     use ModB0, ONLY: B0_DGB
<     use ModAdvance, ONLY: State_VGB, Bx_, Bz_, Hyp_, HypE_
---
>     use ModPhysics,   ONLY: No2Si_V, UnitX_, Si2No_V, iUnitCons_V
>     use ModMain,      ONLY: UseB0, UseHyperbolicDivB
>     use ModB0,        ONLY: B0_DGB
>     use ModAdvance,   ONLY: State_VGB, Bx_, Bz_, Hyp_, HypE_
330c331
<        do iBlock = 1, nBlock
---
>        do iBlock=1, nBlock
333c334
<           do k = 1, nK; do j = 1, nJ; do i = 1, nI
---
>           do k = 1,nK; do j = 1,nJ; do i=1,nI
402c403,404
<   subroutine GM_put_from_pc(NameVar, nVar, nPoint, Data_VI, iPoint_I, Pos_DI)
---
>   subroutine GM_put_from_pc( &
>        NameVar, nVar, nPoint, Data_VI, iPoint_I, Pos_DI)
404c406
<     use BATL_lib, ONLY: nDim
---
>     use BATL_lib,    ONLY: nDim
diff -r -x .git SWMF/GM/BATSRUS/srcInterface/GM_couple_ua.f90 SWMF_yesterday/GM/BATSRUS/srcInterface/GM_couple_ua.f90
20c20
<     ! This function will be called 3 times:
---
>     ! This function will be called 3 times :
26c26
<     ! 3) Receive Data_VI from GM and put them into State_VGB.
---
>     ! 3) Recieve Data_VI from GM and put them into State_VGB.
30,34c30,33
<     use BATL_lib, ONLY: Xyz_DGB, nBlock, Unused_B, nI, nJ, nK, MaxBlock, &
<          test_start ! , iTest, jTest, kTest, iBlockTest
<     use ModMain, ONLY: UaState_VCB
<     use ModGeometry, ONLY: r_GB
<     use ModPhysics, ONLY: No2Si_V, UnitX_, rBody
---
>     use BATL_lib,     ONLY: Xyz_DGB, nBlock, Unused_B, nI, nJ, nK, MaxBlock
>     use ModMain,      ONLY: UaState_VCB
>     use ModGeometry,  ONLY: r_GB
>     use ModPhysics,   ONLY: No2Si_V, UnitX_, rBody
47d45
<     logical:: DoTest ! DoTestCell
50,51d47
<     call test_start(NameSub, DoTest)
< 
55a52,54
>     ! AltMaxUA = Grid_C(UA_)%Coord1_I(2)
>     ! rMaxUA = (rPlanetSi + AltMaxUA)*Si2No_V(UnitX_)
> 
68,71d66
<           ! DoTestCell = DoTest .and. iBlock == iBlockTest .and. &
<           !     i == iTest .and. j == jTest .and. k == kTest
<           ! if(DoTestCell) write(*,*) NameSub,': iPoint=', iPoint
< 
75,77c70
<              UaState_VCB(:5,i,j,k,iBlock) = Data_VI(:5,iPoint_I(iPoint))
<              ! if(DoTestCell) write(*,*) NameSub, ': UaState=', &
<              !     UaState_VCB(:5,i,j,k,iBlock)
---
>              UaState_VCB(:5,i,j,k,iBlock) = Data_VI(:5,iPoint)
81d73
<              ! if(DoTestCell) write(*,*) NameSub, ': Pos_D=', Pos_DI(:,iPoint)
95c87
<     use ModMain, ONLY: IsNewUaState
---
>     use ModMain,      ONLY: IsNewUaState
diff -r -x .git SWMF/GM/BATSRUS/srcInterface/GM_wrapper.f90 SWMF_yesterday/GM/BATSRUS/srcInterface/GM_wrapper.f90
126c126
<     call CON_set_do_test(NameSub, DoTest, DoTestMe)
---
>     call CON_set_do_test(NameSub,DoTest,DoTestMe)
128,129c128,129
<     if(DoTest)write(*,*)NameSub,' called with TypeAction,iProc=', &
<          TypeAction, iProc
---
>     if(DoTest)write(*,*)NameSub,' called with TypeAction,iProc=',&
>          TypeAction,iProc
133c133,135
<        call put(CompInfo, Use=.true., NameVersion='BATSRUS')
---
>        call put(CompInfo,                              &
>             Use        =.true.,                        &
>             NameVersion='BATSRUS (Univ. of Michigan)')
222c224
<     call CON_set_do_test('test_grids', DoTest, DoTestMe)
---
>     call CON_set_do_test('test_grids',DoTest,DoTestMe)
281,282d282
< 
<     logical:: DoTest, DoTestMe
285d284
<     call CON_set_do_test(NameSub, DoTest, DoTestMe)
291,293d289
<     if(DoTestMe)write(*,*) NameSub,': nDimOut, iGridOut, iDecompOut=', &
<          nDimOut, iGridOut, iDecompOut
<     
337c333
<     call CON_set_do_test(NameSub, DoTest, DoTestMe)
---
>     call CON_set_do_test(NameSub,DoTest, DoTestMe)
418c414
<     call CON_set_do_test(NameSub, DoTest, DoTestMe)
---
>     call CON_set_do_test(NameSub,DoTest,DoTestMe)
diff -r -x .git SWMF/output/test12/GM.log SWMF_yesterday/output/test12/GM.log
3,12c3,12
<        1  0.00000E+00  0.00000E+000  2.82470E+001  1.96765E+000  6.70884E-001  2.17171E-003  1.08552E-001 -7.53007E-017  2.76758E-018 -6.84350E-019 -6.63380E-002  1.45484E-001 -4.71319E-002  3.37500E-011  8.75764E-003  1.59842E+000
<        2  0.00000E+00  0.00000E+000  2.82833E+001  2.00521E+000  6.70929E-001  2.19666E-003  1.08482E-001 -1.08670E-016  4.06396E-018 -1.00024E-018 -9.51797E-002  2.09699E-001 -6.78262E-002  4.28035E-011  8.61148E-003  1.59842E+000
<        3  0.00000E+00  0.00000E+000  2.83166E+001  2.03971E+000  6.70975E-001  2.22197E-003  1.08411E-001 -1.36728E-016  5.17993E-018 -1.27305E-018 -1.19686E-001  2.65348E-001 -8.56877E-002  5.05104E-011  8.47653E-003  1.59841E+000
<        4  0.00000E+00  0.00000E+000  2.83505E+001  2.07491E+000  6.71022E-001  2.24763E-003  1.08339E-001 -1.62388E-016  6.21585E-018 -1.52724E-018 -1.41471E-001  3.16276E-001 -1.01996E-001  5.76167E-011  8.35091E-003  1.59841E+000
<        5  0.00000E+00  0.00000E+000  2.83859E+001  2.11098E+000  6.71056E-001  2.29049E-003  1.08283E-001 -1.86311E-016  7.19121E-018 -1.76772E-018 -1.61266E-001  3.63686E-001 -1.17149E-001  6.42997E-011  8.23596E-003  1.59840E+000
<        6  0.00000E+00  0.00000E+000  2.84221E+001  2.14790E+000  6.71091E-001  2.33428E-003  1.08226E-001 -2.08866E-016  8.11537E-018 -1.99727E-018 -1.79520E-001  4.08258E-001 -1.31365E-001  7.06557E-011  8.12213E-003  1.59840E+000
<        7  0.00000E+00  0.00000E+000  2.84591E+001  2.18566E+000  6.71126E-001  2.37903E-003  1.08167E-001 -2.30297E-016  8.99497E-018 -2.21692E-018 -1.96534E-001  4.50466E-001 -1.44798E-001  7.67476E-011  8.01016E-003  1.59840E+000
<        8  0.00000E+00  0.00000E+000  2.84969E+001  2.22416E+000  6.71162E-001  2.42477E-003  1.08108E-001 -2.50778E-016  9.83706E-018 -2.42793E-018 -2.12526E-001  4.90647E-001 -1.57564E-001  8.26208E-011  7.90397E-003  1.59839E+000
<        9  0.00000E+00  0.00000E+000  2.85353E+001  2.26328E+000  6.71199E-001  2.47152E-003  1.08047E-001 -2.70441E-016  1.06437E-017 -2.63105E-018 -2.27651E-001  5.29048E-001 -1.69752E-001  8.83086E-011  7.80335E-003  1.59839E+000
<       10  0.00000E+00  0.00000E+000  2.85740E+001  2.30279E+000  6.71237E-001  2.51930E-003  1.07985E-001 -2.89382E-016  1.14187E-017 -2.82661E-018 -2.42022E-001  5.65870E-001 -1.81437E-001  9.38327E-011  7.70830E-003  1.59838E+000
---
>        1  0.00000E+00  0.00000E+000  2.79501E+001  1.96765E+000  6.64351E-001  2.15180E-003  1.06564E-001 -7.53145E-017  2.76751E-018 -6.85986E-019 -6.63380E-002  1.45484E-001 -4.71319E-002  3.37499E-011  7.12802E-003  1.54639E+000
>        2  0.00000E+00  0.00000E+000  2.77025E+001  2.00521E+000  6.58197E-001  2.15662E-003  1.04556E-001 -1.08698E-016  4.06391E-018 -1.00381E-018 -9.51797E-002  2.09699E-001 -6.78262E-002  4.28034E-011  5.92307E-003  1.49605E+000
>        3  0.00000E+00  0.00000E+000  2.74521E+001  2.03971E+000  6.52001E-001  2.16153E-003  1.02586E-001 -1.36770E-016  5.17975E-018 -1.27868E-018 -1.19686E-001  2.65348E-001 -8.56877E-002  5.05102E-011  4.93640E-003  1.44735E+000
>        4  0.00000E+00  0.00000E+000  2.72016E+001  2.07491E+000  6.45731E-001  2.16653E-003  1.00651E-001 -1.62444E-016  6.21537E-018 -1.53496E-018 -1.41471E-001  3.16276E-001 -1.01996E-001  5.76164E-011  4.17219E-003  1.40041E+000
>        5  0.00000E+00  0.00000E+000  2.69426E+001  2.11098E+000  6.39133E-001  2.17621E-003  9.87400E-002 -1.86379E-016  7.18983E-018 -1.77675E-018 -1.61266E-001  3.63686E-001 -1.17149E-001  6.42993E-011  3.68200E-003  1.35499E+000
>        6  0.00000E+00  0.00000E+000  2.66861E+001  2.14790E+000  6.32536E-001  2.18609E-003  9.68655E-002 -2.08946E-016  8.11306E-018 -2.00758E-018 -1.79520E-001  4.08258E-001 -1.31365E-001  7.06553E-011  3.37139E-003  1.33734E+000
>        7  0.00000E+00  0.00000E+000  2.64325E+001  2.18566E+000  6.25954E-001  2.19619E-003  9.50274E-002 -2.30388E-016  8.99173E-018 -2.22850E-018 -1.96534E-001  4.50466E-001 -1.44798E-001  7.67471E-011  2.94156E-003  1.33734E+000
>        8  0.00000E+00  0.00000E+000  2.61820E+001  2.22417E+000  6.19396E-001  2.20649E-003  9.32252E-002 -2.50881E-016  9.83290E-018 -2.44077E-018 -2.12526E-001  4.90647E-001 -1.57564E-001  8.26202E-011  2.64107E-003  1.33733E+000
>        9  0.00000E+00  0.00000E+000  2.59348E+001  2.26328E+000  6.12871E-001  2.21703E-003  9.14585E-002 -2.70555E-016  1.06386E-017 -2.64512E-018 -2.27652E-001  5.29048E-001 -1.69752E-001  8.83079E-011  2.43980E-003  1.33733E+000
>       10  0.00000E+00  0.00000E+000  2.56906E+001  2.30279E+000  6.06381E-001  2.22780E-003  8.97266E-002 -2.89507E-016  1.14127E-017 -2.84190E-018 -2.42022E-001  5.65870E-001 -1.81438E-001  9.38319E-011  2.30368E-003  1.33733E+000
diff -r -x .git SWMF/output/test12/UA.log SWMF_yesterday/output/test12/UA.log
19,28c19,28
<        2 2015  3  8 10  0  4 584  4.5842    169.81615    170.21887    169.99869    -13.02736     15.62500      1.44275    125.0    125.0
<        3 2015  3  8 10  0  8 987  4.4037    169.69296    170.29690    169.99016    -20.50833     21.48438      2.99014    125.0    125.0
<        4 2015  3  8 10  0 13 326  4.3388    162.61081    170.33136    168.74579    -23.31302     23.68164      4.37333    125.0    125.0
<        5 2015  3  8 10  0 17 640  4.3143    162.60040    170.34670    168.75554    -24.36720     24.50562      5.36325    125.0    125.0
<        6 2015  3  8 10  0 21 945  4.3046    162.58796    170.38226    168.75656    -24.76266     24.81461      6.13729    125.0    125.0
<        7 2015  3  8 10  0 26 245  4.3004    162.57311    170.45826    168.75326    -24.91093     24.93048      6.77878    125.0    125.0
<        8 2015  3  8 10  0 30 544  4.2980    162.55538    170.53247    168.74838    -24.96654     24.97393      7.33106    125.0    125.0
<        9 2015  3  8 10  0 34 840  4.2966    162.53527    170.60517    168.74317    -24.98743     24.99022      7.81676    125.0    125.0
<       10 2015  3  8 10  0 39 136  4.2955    162.51357    170.67641    168.73797    -24.99528     24.99633      8.24784    125.0    125.0
<       11 2015  3  8 10  0 43 430  4.2946    162.49168    170.74622    168.73277    -24.99823     24.99863      8.63183    125.0    125.0
---
>        2 2015  3  8 10  0  4 584  4.5842    169.81613    170.21872    169.99953    -13.02736     15.62500      1.44275    125.0    125.0
>        3 2015  3  8 10  0  8 987  4.4037    169.69294    170.29606    169.99180    -20.47994     21.48438      3.01861    125.0    125.0
>        4 2015  3  8 10  0 13 326  4.3388    162.61051    170.32995    168.74782    -23.30399     23.68164      4.40700    125.0    125.0
>        5 2015  3  8 10  0 17 640  4.3143    162.60008    170.34492    168.75837    -24.36324     24.50562      5.39797    125.0    125.0
>        6 2015  3  8 10  0 21 945  4.3046    162.58760    170.38226    168.76019    -24.76115     24.81461      6.17202    125.0    125.0
>        7 2015  3  8 10  0 26 245  4.3004    162.57270    170.45826    168.75770    -24.91043     24.93048      6.81314    125.0    125.0
>        8 2015  3  8 10  0 30 544  4.2980    162.55492    170.53247    168.75362    -24.96641     24.97393      7.36460    125.0    125.0
>        9 2015  3  8 10  0 34 840  4.2966    162.53475    170.60517    168.74922    -24.98740     24.99022      7.84941    125.0    125.0
>       10 2015  3  8 10  0 39 136  4.2955    162.51298    170.67641    168.74483    -24.99527     24.99633      8.27949    125.0    125.0
>       11 2015  3  8 10  0 43 430  4.2946    162.49102    170.74622    168.74045    -24.99822     24.99863      8.66246    125.0    125.0
Only in SWMF/UA/GITM: share
Only in SWMF/UA/GITM: util
diff -r -x .git SWMF/UA/MGITM/src/init_msis.Mars.f90 SWMF_yesterday/UA/MGITM/src/init_msis.Mars.f90
33,34c33,34
<   ! nO2 = InNDensityS(i,iO2_)
<   ! nO   = InNDensityS(i,iO_)
---
> ! nO2 = InNDensityS(i,iO2_)
> ! nO   = InNDensityS(i,iO_)
39,46c39,46
<   ! m = (nCO2 * mass(iCO2_) + &
<   !      nO2 * mass(iO2_) + &
<   !      nN2 * mass(iN2_) + &
<   !      nCO * mass(iCO_)) / (nCO2 + nO2 + nN2 + nCO)
<   ! m = (nCO2 * mass(iCO2_) + &
<   !      nO   * mass(iO_) + &
<   !      nN2 * mass(iN2_) + &
<   !      nCO * mass(iCO_)) / (nCO2 + nO  + nN2 + nCO)
---
> ! m = (nCO2 * mass(iCO2_) + &
> !      nO2 * mass(iO2_) + &
> !      nN2 * mass(iN2_) + &
> !      nCO * mass(iCO_)) / (nCO2 + nO2 + nN2 + nCO)
> ! m = (nCO2 * mass(iCO2_) + &
> !      nO   * mass(iO_) + &
> !      nN2 * mass(iN2_) + &
> !      nCO * mass(iCO_)) / (nCO2 + nO  + nN2 + nCO)
54c54
<   !  g = Gravitational_Constant * (RBody/r) ** 2
---
> !  g = Gravitational_Constant * (RBody/r) ** 2
140a141,142
>      write(*,*) '==> Now Initializing Mars Background Composition', iBlock
> 
301a304
>      write(*,*) '============> init_msis.Mars.f90 Major Diagnostics:  Begin'
302a306
>      !\
404a409,410
>      write(*,*) '==> Now Completing Mars Background Composition: END', iBlock
> 
410c416
< !==============================================================================
---
> 
412c418
<      F107A,F107,AP,LogNS, Temp, LogRho)
---
>              F107A,F107,AP,LogNS, Temp, LogRho)
420c426
< !==============================================================================
---
> 
438,441c444,447
<   call readDustHeader
<   do iBlock = 1, nBlocks
<      call setTau(iBlock)
<   enddo
---
>    call readDustHeader
>    do iBlock = 1, nBlocks
>       call setTau(iBlock)
>    enddo
443c449
<   call cleanDust
---
>    call cleanDust
447c453,454
< !==============================================================================
---
> 
> 
498c505,507
< !==============================================================================
---
> 
> 
> 
519c528
<   open(unit=iInputUnit_,file=cDustFile)
---
>  open(unit=iInputUnit_,file=cDustFile)
528c537
<   if (DustFileType .eq. "FullHorizontal") then
---
> if (DustFileType .eq. "FullHorizontal") then
530,541c539,549
<      !The flies are ordered by lats first then lons.  I.e.,
<      do while (iError .eq. 0)
<         read(iInputUnit_,*,iostat=iError) TimeArray(1:6),Temp
<         if (iproc .eq. 0) then
<            write(*,*) TimeArray(1:6)
<         endif
<         i = 1
<         do iLat = 1, ndustlats
<            do iLon =1 ,ndustlons
<               TempDust(iLat,iLon) = Temp(i)
<               i = i + 1
<            enddo
---
>   !The flies are ordered by lats first then lons.  I.e.,
>   do while (iError .eq. 0)
>      read(iInputUnit_,*,iostat=iError) TimeArray(1:6),Temp
>      if (iproc .eq. 0) then
>         write(*,*) TimeArray(1:6)
>      endif
>      i = 1
>      do iLat = 1, ndustlats
>         do iLon =1 ,ndustlons
>            TempDust(iLat,iLon) = Temp(i)
>            i = i + 1
543,567c551,553
<         TimeArray(7) = 0
<         call time_int_to_real(TimeArray,TimeDust(iLine))
< 
<         do iLat = 1, nLats
<            do iLon = 1, nLons
< 
<               latFind = Latitude(ilat,iBlock)*180/pi
<               lonFind = Longitude(ilon,iBlock)*180/pi
< 
<               templat = DustLatitude
<               templon = DustLongitude
< 
<               where(LatFind - tempLat .lt. -0.00001) tempLat = -1.0e9
<               ilatlow =  maxloc(tempLat)
< 
<               where(LonFind - tempLon < -0.00001) tempLon = -1.0e9
<               ilonlow = maxloc(tempLon)
< 
<               if (ilatlow(1) .eq. nDustLats) ilatlow = ilatlow - 1
<               if (ilonlow(1) .eq. nDustLons) ilonlow = ilonlow - 1
< 
<               Latlow = DustLatitude(ilatlow(1))
<               LatHigh = DustLatitude(ilatlow(1)+1)
<               Lonlow = DustLongitude(ilonlow(1))
<               LonHigh = DustLongitude(ilonlow(1)+1)
---
>      enddo
>      TimeArray(7) = 0
>      call time_int_to_real(TimeArray,TimeDust(iLine))
569,570c555,556
<               invLatdiff = 1/(Lathigh - Latlow)
<               invLondiff = 1/(Lonhigh - Lonlow)
---
>      do iLat = 1, nLats
>         do iLon = 1, nLons
572,577c558,559
<               if (LatFind .lt. DustLatitude(1) .or. LonFind .lt. DustLongitude(1) &
<                    .or. LatFind .gt. DustLatitude(nDustLats) .or. LonFind .ge. DustLongitude(nDustLons)) then
<                  write(*,*) 'Dust grid does not cover GITM grid'
<                  write(*,*) 'Stopping...'
<                  call stop_gitm('Stopping in init_msis.Mars')
<               endif
---
>            latFind = Latitude(ilat,iBlock)*180/pi
>            lonFind = Longitude(ilon,iBlock)*180/pi
578a561,562
>            templat = DustLatitude
>            templon = DustLongitude
580,598c564,607
<               V11 = TempDust(ilatlow(1),ilonlow(1))
<               V12 = TempDust(ilatlow(1)+1,ilonlow(1))
<               V21 = TempDust(ilatlow(1),ilonlow(1)+1)
<               V22 = TempDust(ilatlow(1)+1,ilonlow(1)+1)
< 
<               rlat1 = (Lathigh - LatFind)*V11*invLatDiff + &
<                    (LatFind-LatLow)*V21*invLatDiff
< 
<               rlat2 = (Lathigh - LatFind)*V12*invLatDiff + &
<                    (LatFind-LatLow)*V22*invLatDiff
< 
<               Dust = (Lonhigh - LonFind)*invLonDiff*rlat1 +  (LonFind - LonLow)*invLonDiff*rlat2
<               !           if (iproc .eq. 1 .and. ilat .eq. 1 .and. ilon .eq. 5 .and. iline .eq. 54) then
<               !              write(*,*) iproc,ilat,ilon,iline,dust,lonhigh,lonfind,lonlow,lathigh,latfind,latlow
<               !              write(*,*) V11, V12, V21,V22
<               !              write(*,*) "lats: ",ilatlow(1)+1, ilonlow(1),ndustlats,ndustlons
<               !              stop
<               !           endif
<               HorizontalDustProfile(iline,iLat,iLon,iblock) = Dust
---
>            where(LatFind - tempLat .lt. -0.00001) tempLat = -1.0e9
>            ilatlow =  maxloc(tempLat)
> 
>            where(LonFind - tempLon < -0.00001) tempLon = -1.0e9
>            ilonlow = maxloc(tempLon)
> 
>            if (ilatlow(1) .eq. nDustLats) ilatlow = ilatlow - 1
>            if (ilonlow(1) .eq. nDustLons) ilonlow = ilonlow - 1
> 
>            Latlow = DustLatitude(ilatlow(1))
>            LatHigh = DustLatitude(ilatlow(1)+1)
>            Lonlow = DustLongitude(ilonlow(1))
>            LonHigh = DustLongitude(ilonlow(1)+1)
> 
>            invLatdiff = 1/(Lathigh - Latlow)
>            invLondiff = 1/(Lonhigh - Lonlow)
> 
>            if (LatFind .lt. DustLatitude(1) .or. LonFind .lt. DustLongitude(1) &
>                 .or. LatFind .gt. DustLatitude(nDustLats) .or. LonFind .ge. DustLongitude(nDustLons)) then
>               write(*,*) 'Dust grid does not cover GITM grid'
>               write(*,*) 'Stopping...'
>               call stop_gitm('Stopping in init_msis.Mars')
>            endif
> 
> 
>            V11 = TempDust(ilatlow(1),ilonlow(1))
>            V12 = TempDust(ilatlow(1)+1,ilonlow(1))
>            V21 = TempDust(ilatlow(1),ilonlow(1)+1)
>            V22 = TempDust(ilatlow(1)+1,ilonlow(1)+1)
> 
>            rlat1 = (Lathigh - LatFind)*V11*invLatDiff + &
>                 (LatFind-LatLow)*V21*invLatDiff
> 
>            rlat2 = (Lathigh - LatFind)*V12*invLatDiff + &
>                 (LatFind-LatLow)*V22*invLatDiff
> 
>            Dust = (Lonhigh - LonFind)*invLonDiff*rlat1 +  (LonFind - LonLow)*invLonDiff*rlat2
> !           if (iproc .eq. 1 .and. ilat .eq. 1 .and. ilon .eq. 5 .and. iline .eq. 54) then
> !              write(*,*) iproc,ilat,ilon,iline,dust,lonhigh,lonfind,lonlow,lathigh,latfind,latlow
> !              write(*,*) V11, V12, V21,V22
> !              write(*,*) "lats: ",ilatlow(1)+1, ilonlow(1),ndustlats,ndustlons
> !              stop
> !           endif
>            HorizontalDustProfile(iline,iLat,iLon,iblock) = Dust
600d608
<            enddo
602,603d609
<         iline = iline + 1
< 
605c611,614
<      nDustTimes = iline-1
---
>      iline = iline + 1
> 
>   enddo
>   nDustTimes = iline-1
607,609c616,618
<      nConrathTimes = nDustTimes
<      TimeConrath = TimeDust
<      HorizontalConrathProfile = 0.03
---
>   nConrathTimes = nDustTimes
>   TimeConrath = TimeDust
>   HorizontalConrathProfile = 0.03
612,614c621,623
<      !Read in data
<      !Files should be yyyy mm dd hh mm ss Latitude pressureLevel(not used) pressure(Pa) CumulativeTau  MixingRatio
<      !and organized by time, lat, alt to fill the dust grid
---
>     !Read in data
>     !Files should be yyyy mm dd hh mm ss Latitude pressureLevel(not used) pressure(Pa) CumulativeTau  MixingRatio
>     !and organized by time, lat, alt to fill the dust grid
619,625c628,634
<               if (iError .ne. 0) then
<                  write(*,*) "Error reading dustfile"
<                  call stop_gitm('Stopping in init_msis.Mars')
<               endif
<               DustPressureLevel(iAlt) = MCSTemp(3)
<               CumulativeTau(iTime,iLat,iAlt) = MCSTemp(4)
<               DustMixingRatio(iTime,iLat,iAlt) = MCSTemp(5)
---
>              if (iError .ne. 0) then
>                 write(*,*) "Error reading dustfile"
>                 call stop_gitm('Stopping in init_msis.Mars')
>              endif
>              DustPressureLevel(iAlt) = MCSTemp(3)
>              CumulativeTau(iTime,iLat,iAlt) = MCSTemp(4)
>              DustMixingRatio(iTime,iLat,iAlt) = MCSTemp(5)
628c637
<            DustLatitude(iLat) = MCSTemp(1)
---
>             DustLatitude(iLat) = MCSTemp(1)
636,638c645,647
<      do iLat = 1, nLats
<         latFind = Latitude(ilat,iBlock)*180/pi
<         templat = DustLatitude
---
>     do iLat = 1, nLats
>       latFind = Latitude(ilat,iBlock)*180/pi
>       templat = DustLatitude
640,641c649,650
<         where(LatFind - tempLat .lt. -0.00001) tempLat = -1.0e9
<         ilatlow =  maxloc(tempLat)
---
>       where(LatFind - tempLat .lt. -0.00001) tempLat = -1.0e9
>       ilatlow =  maxloc(tempLat)
643c652
<         if (ilatlow(1) .eq. nDustLats) ilatlow = ilatlow - 1
---
>       if (ilatlow(1) .eq. nDustLats) ilatlow = ilatlow - 1
645,647c654,656
<         Latlow = DustLatitude(ilatlow(1))
<         LatHigh = DustLatitude(ilatlow(1)+1)
<         invLatdiff = 1/(Lathigh - Latlow)
---
>       Latlow = DustLatitude(ilatlow(1))
>       LatHigh = DustLatitude(ilatlow(1)+1)
>       invLatdiff = 1/(Lathigh - Latlow)
649,651c658,660
<         if (LatFind .lt. DustLatitude(1)) then
<            CumulativeTauProfile(1:nDustTimes,ilat,1:nDustAlts,iblock) = &
<                 CumulativeTau(1:nDustTimes,1,1:nDustAlts)
---
>       if (LatFind .lt. DustLatitude(1)) then
>           CumulativeTauProfile(1:nDustTimes,ilat,1:nDustAlts,iblock) = &
>                CumulativeTau(1:nDustTimes,1,1:nDustAlts)
653,654c662,663
<            DustMixingRatioProfile(1:nDustTimes,ilat,1:nDustAlts,iblock) = &
<                 DustMixingRatio(1:nDustTimes,1,1:nDustAlts)
---
>           DustMixingRatioProfile(1:nDustTimes,ilat,1:nDustAlts,iblock) = &
>                DustMixingRatio(1:nDustTimes,1,1:nDustAlts)
656c665
<         else if (LatFind .gt. DustLatitude(nDustLats)) then
---
>        else if (LatFind .gt. DustLatitude(nDustLats)) then
658,659c667,668
<            CumulativeTauProfile(1:nDustTimes,ilat,1:nDustAlts,iblock) = &
<                 CumulativeTau(1:nDustTimes,nDustLats,1:nDustAlts)
---
>           CumulativeTauProfile(1:nDustTimes,ilat,1:nDustAlts,iblock) = &
>                CumulativeTau(1:nDustTimes,nDustLats,1:nDustAlts)
661,662c670,671
<            DustMixingRatioProfile(1:nDustTimes,ilat,1:nDustAlts,iblock) = &
<                 DustMixingRatio(1:nDustTimes,nDustLats,1:nDustAlts)
---
>           DustMixingRatioProfile(1:nDustTimes,ilat,1:nDustAlts,iblock) = &
>                DustMixingRatio(1:nDustTimes,nDustLats,1:nDustAlts)
668,683c677,692
<            CumulativeTauProfile(1:nDustTimes,ilat,1:nDustAlts,iblock) = &
<                 CumulativeTau(1:nDustTimes,ilatlow(1)+1,1:nDustAlts) - &
<                 (Lathigh-LatFind)*0.999*invLatDiff*&
<                 (CumulativeTau(1:nDustTimes,ilatlow(1)+1,1:nDustAlts)-&
<                 CumulativeTau(1:nDustTimes,ilatlow(1),1:nDustAlts))
< 
<            DustMixingRatioProfile(1:nDustTimes,ilat,1:nDustAlts,iblock) = &
<                 DustMixingRatio(1:nDustTimes,ilatlow(1)+1,1:nDustAlts) - &
<                 (Lathigh-LatFind)*0.999*invLatDiff*&
<                 (DustMixingRatio(1:nDustTimes,ilatlow(1)+1,1:nDustAlts)-&
<                 DustMixingRatio(1:nDustTimes,ilatlow(1),1:nDustAlts))
<         endif
< 
<      enddo
<   endif
< 
---
>           CumulativeTauProfile(1:nDustTimes,ilat,1:nDustAlts,iblock) = &
>                CumulativeTau(1:nDustTimes,ilatlow(1)+1,1:nDustAlts) - &
>                (Lathigh-LatFind)*0.999*invLatDiff*&
>                (CumulativeTau(1:nDustTimes,ilatlow(1)+1,1:nDustAlts)-&
>                CumulativeTau(1:nDustTimes,ilatlow(1),1:nDustAlts))
> 
>           DustMixingRatioProfile(1:nDustTimes,ilat,1:nDustAlts,iblock) = &
>                DustMixingRatio(1:nDustTimes,ilatlow(1)+1,1:nDustAlts) - &
>                (Lathigh-LatFind)*0.999*invLatDiff*&
>                (DustMixingRatio(1:nDustTimes,ilatlow(1)+1,1:nDustAlts)-&
>                DustMixingRatio(1:nDustTimes,ilatlow(1),1:nDustAlts))
>        endif
> 
>     enddo
>  endif
>  
697,735c706,744
<   !        templat = DustLatitude
<   !        where(LatFind - tempLat .lt. -0.00001) tempLat = -1.0e9
<   !
<   !        ilatlow =  maxloc(tempLat)
<   !
<   !        if (ilatlow(1) .eq. nDustLats) ialtlow = ialtlow - 1
<   !
<   !        Latlow = DustLatitude(ilatlow(1))
<   !        LatHigh = DustLatitude(ilatlow(1)+1)
<   !
<   !        invLatdiff = 1/(Lathigh - Latlow)
<   !
<   !        if (LatFind .lt. DustLatitude(1)) then
<   !
<   !           rlat = (latlow-latFind)*(TempConrath(ilatlow(1) + 1)-TempConrath(ilatlow(1))) * &
<   !                invLatDiff
<   !           Conrath = TempConrath(ilat) - rlat
<   !
<   !        else
<   !
<   !           if (LatFind .ge. DustLatitude(nDustLats)) then
<   !
<   !              rlat = (LatFind-Lathigh)* &
<   !                   (TempConrath(ilatlow(1) + 1) - TempConrath(ilatlow(1))) * invLatDiff
<   !
<   !              Conrath = TempConrath(ilatlow(1) + 1) + rlat
<   !
<   !           else
<   !              rlat = (Lathigh - LatFind)*(TempConrath(ilatlow(1) + 1) - &
<   !                   TempConrath(ilatlow(1))) * invLatDiff
<   !              Conrath = TempConrath(ilatlow(1) + 1) - rlat
<   !
<   !           endif
<   !        endif
<   !        HorizontalConrathProfile(iline,iLat,iblock) = Conrath
<   !
<   !     enddo
<   !
<   !     iline = iline + 1
---
>  !        templat = DustLatitude
>  !        where(LatFind - tempLat .lt. -0.00001) tempLat = -1.0e9
>  !
>  !        ilatlow =  maxloc(tempLat)
>  !
>  !        if (ilatlow(1) .eq. nDustLats) ialtlow = ialtlow - 1
>  !
>  !        Latlow = DustLatitude(ilatlow(1))
>  !        LatHigh = DustLatitude(ilatlow(1)+1)
>  !
>  !        invLatdiff = 1/(Lathigh - Latlow)
>  !
>  !        if (LatFind .lt. DustLatitude(1)) then
>  !
>  !           rlat = (latlow-latFind)*(TempConrath(ilatlow(1) + 1)-TempConrath(ilatlow(1))) * &
>  !                invLatDiff
>  !           Conrath = TempConrath(ilat) - rlat
>  !
>  !        else
>  !
>  !           if (LatFind .ge. DustLatitude(nDustLats)) then
>  !
>  !              rlat = (LatFind-Lathigh)* &
>  !                   (TempConrath(ilatlow(1) + 1) - TempConrath(ilatlow(1))) * invLatDiff
>  !
>  !              Conrath = TempConrath(ilatlow(1) + 1) + rlat
>  !
>  !           else
>  !              rlat = (Lathigh - LatFind)*(TempConrath(ilatlow(1) + 1) - &
>  !                   TempConrath(ilatlow(1))) * invLatDiff
>  !              Conrath = TempConrath(ilatlow(1) + 1) - rlat
>  !
>  !           endif
>  !        endif
>  !        HorizontalConrathProfile(iline,iLat,iblock) = Conrath
>  !
>  !     enddo
>  !
>  !     iline = iline + 1
739,740d747
< !==============================================================================
< subroutine cleanDust
742d748
<   use GITM_planet, only : DustLatitude,DustLongitude
744,749c750,758
<   if (allocated(DustLatitude)) then
<      deallocate(DustLatitude)
<   endif
<   if (allocated(DustLongitude)) then
<      deallocate(DustLongitude)
<   endif
---
> subroutine cleanDust
> use GITM_planet, only : DustLatitude,DustLongitude
> 
> if (allocated(DustLatitude)) then
>    deallocate(DustLatitude)
> endif
> if (allocated(DustLongitude)) then
>    deallocate(DustLongitude)
> endif
752d760
< !==============================================================================
diff -r -x .git SWMF/UA/MGITM/src/ModMars.f90 SWMF_yesterday/UA/MGITM/src/ModMars.f90
7,11c7,11
<   ! Modified (01/18/07) : SWB :   Aij, s-exponents for mutual diffusion
<   ! Modified (06/12/08) : SWB :   ordering to species revised
<   ! Modified (06/12/08) : SWB :   nSpecies = 6; nSpeciesTotal = 11
<   ! Modified (04/21/17) : SWB :   nEmissions = 10; iENOUV = 1
<   ! Majors (6):  COntrol the Pressures Gradients and winds
---
> ! Modified (01/18/07) : SWB :   Aij, s-exponents for mutual diffusion
> ! Modified (06/12/08) : SWB :   ordering to species revised
> ! Modified (06/12/08) : SWB :   nSpecies = 6; nSpeciesTotal = 11
> ! Modified (04/21/17) : SWB :   nEmissions = 10; iENOUV = 1
> ! Majors (6):  COntrol the Pressures Gradients and winds
22c22
<   ! Minors (6) : Ride on background of mean winds derived from majors
---
> ! Minors (6) : Ride on background of mean winds derived from majors
29,31c29,31
<   ! Major Ions (5):  Most Important to MWACM code
<   ! Modified (05/21/08) : SWB :   Add N2+ to major ions list
<   ! DJP (7/21/2025): Added CO+ and C+
---
> ! Major Ions (5):  Most Important to MWACM code
> ! Modified (05/21/08) : SWB :   Add N2+ to major ions list
> ! DJP (7/21/2025): Added CO+ and C+
60,65c60,65
<   ! integer, parameter :: iE7320_ = 2
<   ! integer, parameter :: iE3726_ = 3
<   ! integer, parameter :: iE5200_ = 4
<   ! integer, parameter :: iE10400_ = 5
<   ! integer, parameter :: iE6300_ = 6
<   ! integer, parameter :: iE6364_ = 7
---
> ! integer, parameter :: iE7320_ = 2
> ! integer, parameter :: iE3726_ = 3
> ! integer, parameter :: iE5200_ = 4
> ! integer, parameter :: iE10400_ = 5
> ! integer, parameter :: iE6300_ = 6
> ! integer, parameter :: iE6364_ = 7
94,98c94,98
<   ! real, parameter :: SunOrbit_A = 1.52
<   ! real, parameter :: SunOrbit_B = 0.04
<   ! real, parameter :: SunOrbit_C = 0.15
<   ! real, parameter :: SunOrbit_D = 0.00
<   ! real, parameter :: SunOrbit_E = 0.00
---
> ! real, parameter :: SunOrbit_A = 1.52
> ! real, parameter :: SunOrbit_B = 0.04
> ! real, parameter :: SunOrbit_C = 0.15
> ! real, parameter :: SunOrbit_D = 0.00
> ! real, parameter :: SunOrbit_E = 0.00
100c100
<   real, parameter :: SunOrbit_A = 1.52369
---
>  real, parameter :: SunOrbit_A = 1.52369
131,132c131,132
<   ! Photochemistry pathways
<   ! Some species have multiple photodissociation / photoionization pathways
---
> ! Photochemistry pathways
> ! Some species have multiple photodissociation / photoionization pathways
135c135
<   ! CO2
---
> ! CO2
139c139
<   !CO
---
> !CO
142c142
<   !N2
---
> !N2
145c145
<   !O2
---
> !O2
148c148
<   !Total
---
> !Total
152c152
<   !CO2
---
> !CO2
160c160
<   !CO
---
> !CO
163c163
<   !O2 
---
> !O2 
166c166
<   !O
---
> !O
169c169
<   !N2
---
> !N2
172c172
<   !Total
---
> !Total
175c175
<   !  real :: KappaTemp0 = 2.22e-4
---
> !  real :: KappaTemp0 = 2.22e-4
177c177
< !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
---
>   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
180c180
< !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
---
>   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
182,214c182,214
<   ! real, parameter, dimension(6, 6) :: Diff0 = 1.0e17 * reshape( (/ &
<   !  ! These are Aij coefficients from B&K (1973) formulation: Aij*1.0E+17
<   !  ! Use Jared Bell's Titan GITM formulation for Mars GITM
<   !! integer, parameter :: iCO2_    = 1
<   !! integer, parameter :: iCO_     = 2
<   !! integer, parameter :: iO_      = 3
<   !! integer, parameter :: iN2_     = 4
<   !! integer, parameter :: iO2_     = 5
<   !! integer, parameter :: iAr_     = 6
<   !     !------------------------------------------------+
<   !     ! i=C02      CO      O        N2      O2      Ar
<   !     !------------------------------------------------+
<   !       0.0000, 0.7762, 0.2219,  0.6580,  0.5770, 1.1920,         &  ! CO2
<   !       0.7762, 0.0000, 0.9466,  0.9280,  0.8300, 0.6625,         &  ! CO
<   !       0.2219, 0.9466, 0.0000,  0.9690,  0.9690, 0.5510,         &  ! O
<   !       0.6580, 0.9280, 0.9690,  0.0000,  0.7150, 0.6640,         &  ! N2
<   !       0.5770, 0.8300, 0.9690,  0.7150,  0.0000, 0.7170,         &  ! O2
<   !       1.1920, 0.6625, 0.5510,  0.6640,  0.7170, 0.000 /), (/6,6/) )! Ar
<   !
<   !  ! These are s-exponents from B&K (1973) formulation: T**s
<   !
<   !   real, parameter, dimension(6, 6) :: DiffExp = reshape( (/ &
<   !     !------------------------------------------------+
<   !     ! i=C02      CO      O     N2     O2     Ar
<   !     !------------------------------------------------+
<   !       0.000,  0.750,  0.750, 0.752, 0.749, 0.750,  &           ! CO2
<   !       0.750,  0.000,  0.750, 0.710, 0.724, 0.750,  &           ! CO
<   !       0.750,  0.750,  0.000, 0.774, 0.774, 0.841,  &           ! O
<   !       0.752,  0.710,  0.774, 0.000, 0.750, 0.752,  &           ! N2
<   !       0.749,  0.724,  0.774, 0.750, 0.000, 0.736,  &           ! O2
<   !       0.750,  0.750,  0.841, 0.752, 0.736, 0.000 /), (/6,6/) ) ! AR
<   !
<   !!     Arrays filled in init_radcool in data statements (np = 68)
---
> ! real, parameter, dimension(6, 6) :: Diff0 = 1.0e17 * reshape( (/ &
> !  ! These are Aij coefficients from B&K (1973) formulation: Aij*1.0E+17
> !  ! Use Jared Bell's Titan GITM formulation for Mars GITM
> !! integer, parameter :: iCO2_    = 1
> !! integer, parameter :: iCO_     = 2
> !! integer, parameter :: iO_      = 3
> !! integer, parameter :: iN2_     = 4
> !! integer, parameter :: iO2_     = 5
> !! integer, parameter :: iAr_     = 6
> !     !------------------------------------------------+
> !     ! i=C02      CO      O        N2      O2      Ar
> !     !------------------------------------------------+
> !       0.0000, 0.7762, 0.2219,  0.6580,  0.5770, 1.1920,         &  ! CO2
> !       0.7762, 0.0000, 0.9466,  0.9280,  0.8300, 0.6625,         &  ! CO
> !       0.2219, 0.9466, 0.0000,  0.9690,  0.9690, 0.5510,         &  ! O
> !       0.6580, 0.9280, 0.9690,  0.0000,  0.7150, 0.6640,         &  ! N2
> !       0.5770, 0.8300, 0.9690,  0.7150,  0.0000, 0.7170,         &  ! O2
> !       1.1920, 0.6625, 0.5510,  0.6640,  0.7170, 0.000 /), (/6,6/) )! Ar
> !
> !  ! These are s-exponents from B&K (1973) formulation: T**s
> !
> !   real, parameter, dimension(6, 6) :: DiffExp = reshape( (/ &
> !     !------------------------------------------------+
> !     ! i=C02      CO      O     N2     O2     Ar
> !     !------------------------------------------------+
> !       0.000,  0.750,  0.750, 0.752, 0.749, 0.750,  &           ! CO2
> !       0.750,  0.000,  0.750, 0.710, 0.724, 0.750,  &           ! CO
> !       0.750,  0.750,  0.000, 0.774, 0.774, 0.841,  &           ! O
> !       0.752,  0.710,  0.774, 0.000, 0.750, 0.752,  &           ! N2
> !       0.749,  0.724,  0.774, 0.750, 0.000, 0.736,  &           ! O2
> !       0.750,  0.750,  0.841, 0.752, 0.736, 0.000 /), (/6,6/) ) ! AR
> !
> !!     Arrays filled in init_radcool in data statements (np = 68)
216c216
< !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
---
>   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
219,253c219,253
< !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
<   !real, parameter, dimension(nSpecies, nSpecies) :: Diff0 = 1.0e17 * reshape( (/ &
<   !     !----------------------------------------------------------+
<   !     ! i=C02      CO      O        N2      O2      Ar    He
<   !     !----------------------------------------------------------+
<   !       0.0000, 0.7762, 0.2219,  0.6580,  0.5770, 1.1920, 2.4292, &  ! CO2
<   !       0.7762, 0.0000, 0.9466,  0.9280,  0.8300, 0.6625,1159.55, &  ! CO
<   !       0.2219, 0.9466, 0.0000,  0.9690,  0.9690, 0.5510, 3.4346, &  ! O
<   !       0.6580, 0.9280, 0.9690,  0.0000,  0.7150, 0.6640,1159.55, &  ! N2
<   !       0.5770, 0.8300, 0.9690,  0.7150,  0.0000, 0.7170, 3.2070, &  ! O2
<   !       1.1920, 0.6625, 0.5510,  0.6640,  0.7170, 0.0000,1000.00, &  ! Ar
<   !       2.4292,1159.55, 3.4346,1159.550,  3.2070,1000.00, 0.0000 /),  &  ! He
<   !       (/nSpecies,nSpecies/) )
<   !
<   !
<   !  ! These are s-exponents from B&K (1973) formulation: T**s
<   ! real, parameter, dimension(nSpecies, nSpecies) :: DiffExp = reshape( (/ &
<   !     !----------------------------------------------------------+
<   !     ! i=C02      CO      O     N2     O2     Ar     He
<   !     !----------------------------------------------------------+
<   !       0.000,  0.750,  0.750, 0.752, 0.749, 0.750, 0.720, &     ! CO2
<   !       0.750,  0.000,  0.750, 0.710, 0.724, 0.750, 0.524, &     ! CO
<   !       0.750,  0.750,  0.000, 0.774, 0.774, 0.841, 0.749, &     ! O
<   !       0.752,  0.710,  0.774, 0.000, 0.750, 0.752, 0.524, &     ! N2
<   !       0.749,  0.724,  0.774, 0.750, 0.000, 0.736, 0.710, &     ! O2
<   !       0.750,  0.750,  0.841, 0.752, 0.736, 0.000, 0.524, &     ! Ar
<   !       0.720,  0.524,  0.749, 0.524, 0.710, 0.524, 0.000 /), &  ! He
<   !       (/nSpecies,nSpecies/) )
< 
< 
< !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
<   real, parameter, dimension(nSpecies, nSpecies) :: Diff0 = 1.0e17 * reshape( (/ &
<                                 !----------------------------------------------------------------+
<                                 ! i=C02      CO      O        N2      O2      Ar    He       N4S
<                                 !----------------------------------------------------------------+
---
>   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
> !real, parameter, dimension(nSpecies, nSpecies) :: Diff0 = 1.0e17 * reshape( (/ &
> !     !----------------------------------------------------------+
> !     ! i=C02      CO      O        N2      O2      Ar    He
> !     !----------------------------------------------------------+
> !       0.0000, 0.7762, 0.2219,  0.6580,  0.5770, 1.1920, 2.4292, &  ! CO2
> !       0.7762, 0.0000, 0.9466,  0.9280,  0.8300, 0.6625,1159.55, &  ! CO
> !       0.2219, 0.9466, 0.0000,  0.9690,  0.9690, 0.5510, 3.4346, &  ! O
> !       0.6580, 0.9280, 0.9690,  0.0000,  0.7150, 0.6640,1159.55, &  ! N2
> !       0.5770, 0.8300, 0.9690,  0.7150,  0.0000, 0.7170, 3.2070, &  ! O2
> !       1.1920, 0.6625, 0.5510,  0.6640,  0.7170, 0.0000,1000.00, &  ! Ar
> !       2.4292,1159.55, 3.4346,1159.550,  3.2070,1000.00, 0.0000 /),  &  ! He
> !       (/nSpecies,nSpecies/) )
> !
> !
> !  ! These are s-exponents from B&K (1973) formulation: T**s
> ! real, parameter, dimension(nSpecies, nSpecies) :: DiffExp = reshape( (/ &
> !     !----------------------------------------------------------+
> !     ! i=C02      CO      O     N2     O2     Ar     He
> !     !----------------------------------------------------------+
> !       0.000,  0.750,  0.750, 0.752, 0.749, 0.750, 0.720, &     ! CO2
> !       0.750,  0.000,  0.750, 0.710, 0.724, 0.750, 0.524, &     ! CO
> !       0.750,  0.750,  0.000, 0.774, 0.774, 0.841, 0.749, &     ! O
> !       0.752,  0.710,  0.774, 0.000, 0.750, 0.752, 0.524, &     ! N2
> !       0.749,  0.724,  0.774, 0.750, 0.000, 0.736, 0.710, &     ! O2
> !       0.750,  0.750,  0.841, 0.752, 0.736, 0.000, 0.524, &     ! Ar
> !       0.720,  0.524,  0.749, 0.524, 0.710, 0.524, 0.000 /), &  ! He
> !       (/nSpecies,nSpecies/) )
> 
> 
>   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
>  real, parameter, dimension(nSpecies, nSpecies) :: Diff0 = 1.0e17 * reshape( (/ &
>      !----------------------------------------------------------------+
>      ! i=C02      CO      O        N2      O2      Ar    He       N4S
>      !----------------------------------------------------------------+
266,269c266,269
<   real, parameter, dimension(nSpecies, nSpecies) :: DiffExp = reshape( (/ &
<                                 !----------------------------------------------------------+
<                                 ! i=C02      CO      O     N2     O2     Ar     He
<                                 !----------------------------------------------------------+
---
>    real, parameter, dimension(nSpecies, nSpecies) :: DiffExp = reshape( (/ &
>      !----------------------------------------------------------+
>      ! i=C02      CO      O     N2     O2     Ar     He
>      !----------------------------------------------------------+
293,294c293
<   integer, parameter:: &
<        nDustLinesMax = 400,nDustPointsMax = 200,nDustLatsMax = 50,nDustTimesMax = 20
---
>   integer, parameter:: nDustLinesMax = 400,nDustPointsMax = 200,nDustLatsMax = 50,nDustTimesMax = 20
302c301
<        CumulativeTauProfile,DustMixingRatioProfile
---
>     CumulativeTauProfile,DustMixingRatioProfile
317c316
< 
---
>   
325,326c324,325
<        (/ 6.6092890e-14, 3.2531944e-14, 2.2676610e-14, 4.6901680e-14, 0.0, &
<        4.1007410e-14, 0.0, 0.0 /)
---
>   (/ 6.6092890e-14, 3.2531944e-14, 2.2676610e-14, 4.6901680e-14, 0.0, &
>   4.1007410e-14, 0.0, 0.0 /)
329,339c328,338
<                                 !----------------------------------------------------------+      
<                                 ! i=CS      C0      C1     C2     C3     C4     C5
<                                 !----------------------------------------------------------+              
<        -193.89664, 140.54578, -39.939533, 5.5677812, -0.38114739, 0.010236020, & !CO2
<        -196.60813, 142.51595, -40.500339, 5.6459650, -0.38648628, 0.010378626, & !CO 
<        -197.56426, 143.27663, -40.739811, 5.6832256, -0.38935008, 0.010465609, & !O  
<        -194.11932, 140.71999, -39.997168, 5.5776025, -0.38198639, 0.010264152, & !N2 
<        0.0,          0.0,         0.0,          0.0,          0.0,        0.0, & !O2 
<        -201.19056, 145.89558, -41.478495, 5.7849232, -0.39618319, 0.010644507, & !AR 
<        0.0,          0.0,         0.0,          0.0,          0.0,        0.0,&  !He 
<        0.0,          0.0,         0.0,          0.0,          0.0,        0.0 /), & !N4S
---
>      !----------------------------------------------------------+                                        
>      ! i=CS      C0      C1     C2     C3     C4     C5                                                  
>      !----------------------------------------------------------+                                        
>         -193.89664, 140.54578, -39.939533, 5.5677812, -0.38114739, 0.010236020, & !CO2                   
>         -196.60813, 142.51595, -40.500339, 5.6459650, -0.38648628, 0.010378626, & !CO                    
>         -197.56426, 143.27663, -40.739811, 5.6832256, -0.38935008, 0.010465609, & !O                     
>       	-194.11932, 140.71999, -39.997168, 5.5776025, -0.38198639, 0.010264152, & !N2                    
>         0.0,          0.0,         0.0,          0.0,          0.0,        0.0, & !O2                    
>       	-201.19056, 145.89558, -41.478495, 5.7849232, -0.39618319, 0.010644507, & !AR                    
>         0.0,          0.0,         0.0,          0.0,          0.0,        0.0,&  !He                    
>       	0.0,          0.0,         0.0,          0.0,          0.0,        0.0 /), & !N4S                
341a341
> 
357,570c357,893
<   !################ Nelli, April 07 ##########################
<   !Setting up parameters needed by the correlated k lower
<   !atmosphere radiation code
< 
<   !     Number of atmospheric layers
<   integer, PARAMETER :: LL_LAYERS  = nAlts
< 
<   !     Number of atmospheric levels:   2 * LL_LAYERS + 3
<   integer, PARAMETER :: LL_LEVELS  = 2*LL_LAYERS+3
< 
<   !C     LL_NLAYRAD is the number of radiation code layers
<   !C     LL_NLEVRAD is the number of radiation code levels.  Level N is the
<   !C
<   integer, PARAMETER :: LL_NLAYRAD  = LL_LAYERS+1
<   integer, PARAMETER :: LL_NLEVRAD  = LL_LAYERS+2
< 
<   ! Bottom layer subsurface temperature
<   real, parameter :: CoreTemp = 175.0
< 
<   ! Surface and subsurface temperature constants
<   real, parameter :: Pa = 5.927E+7
<   real, parameter :: Pd = 88775.0
< 
<   ! Stefan-Boltzmann constant in SI
<   real, parameter :: SBconstant = 5.67E-8
< 
<   !      Ls variable
<   real :: ell_s = 0.0
<   !C======================================================================C
<   !C
<   !C     RADINC.H    RADiation INCludes
<   !C
<   !C     Includes for the radiation code; RADIATION LAYERS, LEVELS,
<   !C     number of spectral intervals. . .
<   !C
<   !C     GCM2.0  Feb 2003
<   !C
<   !C======================================================================C
< 
<   !C     RADIATION parameters
< 
<   !C     In radiation code, layer 1 corresponds to the stratosphere.  Level
<   !C     1 is the top of the stratosphere.  The dummy layer is at the same
<   !C     temperature as the (vertically isothermal) stratosphere, and
<   !C     any time it is explicitly needed, the appropriate quantities will
<   !C     be dealt with (aka "top". . .)
< 
< 
<   !C
<   !C     L_NSPECTI is the number of IR spectral intervals
<   !C     L_NSPECTV is the number of Visual(or Solar) spectral intervals
<   !C     L_NGAUSS  is the number of Gauss points for K-coefficients
<   !C               GAUSS POINT 9 (aka the last one) is the special case
<   !C     L_NWNGI   is L_NSPECTI*L_NGAUSS;  the total number of "intervals"
<   !C               in the IR
<   !C     L_NWNGV   is L_NSPECTV*L_NGAUSS;  the total number of "intervals"
<   !C               in the VISUAL
<   !C
<   !C     L_NPREF   is the number of reference pressures that the
<   !C               k-coefficients are calculated on
<   !C     L_PINT    is the number of Lagrange interpolated reference
<   !C               pressures for the CO2 k-coefficients.
<   !C     L_NTREF   is the number of refernce temperatures for the
<   !C               k-coefficients
<   !C     L_TAUMAX  is the largest optical depth - larger ones are set
<   !C               to this value.
<   !C
<   !C     L_REFH2O  The number of different water-mixing ratio values for
<   !C               the k-coefficients that are now CO2+H2O.
<   !C
<   !C     L_NREFV   The spectral interval number of the visible reference
<   !C               wavelength (i.e. the 0.67 micron band)
<   !C
<   !C
<   !C     L_NLTE    The number of different LTE/NLTE heating rate conversion
<   !C               factors
<   !C
<   !C----------------------------------------------------------------------C
< 
< 
<   integer, parameter :: L_NSPECTI =  5
<   integer, parameter :: L_NSPECTV =  7
<   integer, parameter :: L_NGAUSS  = 17
< 
<   integer, parameter :: L_NPREF   = 11
<   integer, parameter :: L_NTREF   =  7
<   integer, parameter :: L_TAUMAX  = 35
< 
<   integer, parameter :: L_PINT    = 51
< 
<   integer, parameter :: L_REFH2O  = 10
< 
<   integer, parameter :: L_NREFV   = 6
< 
<   integer, parameter :: L_NLTE    = 37
< 
<   !C----------------------------------------------------------------------C
<   !C
<   !C                             radcommon.h
<   !C                         FORTRAN PARAMETERS
<   !C                          GCM2.0  Feb 2003
<   !C
<   !C----------------------------------------------------------------------C
<   !C
<   !C  "Include" grid.h and radinc.h before this file in code that uses
<   !C  some or all of this common data set
<   !C
<   !C     WNOI       - Array of wavenumbers at the spectral interval
<   !C                  centers for the infrared.  Array is NSPECTI
<   !C                  elements long.
<   !C     DWNI       - Array of "delta wavenumber", i.e., the width,
<   !C                  in wavenumbers (cm^-1) of each IR spectral
<   !C                  interval.  NSPECTI elements long.
<   !C     WAVEI      - Array (NSPECTI elements long) of the wavelenght
<   !C                  (in microns) at the center of each IR spectral
<   !C                  interval.
<   !C     WNOV       - Array of wavenumbers at the spectral interval
<   !C                  center for the VISUAL.  Array is NSPECTV
<   !C                  elements long.
<   !C     DWNV       - Array of "delta wavenumber", i.e., the width,
<   !C                  in wavenumbers (cm^-1) of each VISUAL spectral
<   !C                  interval.  NSPECTV elements long.
<   !C     WAVEV      - Array (NSPECTV elements long) of the wavelenght
<   !C                  (in microns) at the center of each VISUAL spectral
<   !C                  interval.
<   !C     SOLARF     - Array (NSPECTV elements) of solar flux (W/M^2) in
<   !C                  each spectral interval.  Values are for 1 AU, and
<   !C                  are scaled to the Mars distance elsewhere.
<   !C     SOL        - Solar flux at Mars (W/M^2)
<   !C     TAURAY     - Array (NSPECTV elements) of the pressure-independent
<   !C                  part of Rayleigh scattering optical depth.
<   !C     PTOP       - Pressure at the top of the radiation code coordinate;
<   !C                  = 0.5*Ptrop
<   !C     FZEROI     - Fraction of zeros in the IR CO2 k-coefficients, for
<   !C                  each temperature, pressure, and spectral interval
<   !C     FZEROV     - Fraction of zeros in the VISUAL CO2 k-coefficients, for
<   !C                  each temperature, pressure, and spectral interval
<   !C
<   !C     AEROSOL RADIATIVE OPTICAL CONSTANTS
<   !C     Values are at the wavelenght interval center
<   !C
<   !C     MIE SCATTERING - Size distribution weighted
<   !C     Qextv    - Extinction efficiency - in the visual.
<   !C     QextREF  - Reference visual wavelength (.67 micron band)
<   !C     Qscatv   - Scattering efficiency - in the visual.
<   !C     WV       - Single scattering albedo - in the visual.
<   !C     GV       - Asymmetry parameter - in the visual.
<   !C
<   !C     Qexti    - Extinction efficiency - in the infrared.
<   !C     Qscati   - Scattering efficiency - in the infrared.
<   !C     WI       - Single scattering albedo - in the infrared.
<   !C     GI       - Asymmetry parameter - in the infrared.
<   !C
<   !C     VIS2IR   - VISIBLE (0.67 micron band) to IR (9 micron band) ratio.
<   !C
<   !C     XLTEFACTOR - correction factor for over-prediction of heating rates
<   !C                  by the LTE code. Occurs at pressures where NLTE processes
<   !C                  become important
<   !C
<   !C     XLTEPRESSURE - pressure regime at which each factor is to be applied
<   !C
<   !C  "Include" grid.h and radinc.h before this file in code that uses
<   !C  some or all of this common data set
< 
<   REAL :: WNOI(L_NSPECTI), DWNI(L_NSPECTI), WAVEI(L_NSPECTI)
<   REAL :: WNOV(L_NSPECTV), DWNV(L_NSPECTV), WAVEV(L_NSPECTV)
<   REAL :: SOLARF(L_NSPECTV), TAURAY(L_NSPECTV),SOL(L_NSPECTV)
< 
<   real :: CO2I(L_NTREF,L_PINT,L_REFH2O,L_NSPECTI,L_NGAUSS)
<   real :: CO2V(L_NTREF,L_PINT,L_REFH2O,L_NSPECTV,L_NGAUSS)
<   real :: FZEROI(L_NSPECTI)
<   real :: FZEROV(L_NSPECTV)
<   real :: PGASREF(L_NPREF), TGASREF(L_NTREF)
< 
<   real :: qextv(L_NSPECTV), qscatv(L_NSPECTV), wv(L_NSPECTV)
<   real :: gv(L_NSPECTV)
<   real :: QextREF, VIS2IR, Cmk, tlimit
< 
<   real :: qexti(L_NSPECTI), qscati(L_NSPECTI), wi(L_NSPECTI)
<   real :: gi(L_NSPECTI)
< 
<   real :: planckir(L_NSPECTI,8501)
< 
<   real :: PTOP, UBARI,  GWEIGHT(L_NGAUSS)
<   real :: PFGASREF(L_PINT)
< 
<   !C  H2O and CO2 k-coefficients mixed
< 
<   real :: WREFCO2(L_REFH2O), WREFH2O(L_REFH2O)
< 
<   !C  LTE/NLTE heating rate conversion factors and associated pressures
< 
<   real :: XLTEFACTOR(L_NLTE),XLTEPRESSURE(L_NLTE)
< 
<   !  These are for the Gauss-split 0.95 case
< 
<   DATA GWEIGHT  / 4.8083554740D-02, 1.0563099137D-01,&
<        1.4901065679D-01, 1.7227479710D-01,&
<        1.7227479710D-01, 1.4901065679D-01,&
<        1.0563099137D-01, 4.8083554740D-02,&
<        2.5307134073D-03, 5.5595258613D-03,&
<        7.8426661469D-03, 9.0670945845D-03,&
<        9.0670945845D-03, 7.8426661469D-03,&
<        5.5595258613D-03, 2.5307134073D-03,  0.0D0 /
< 
<   DATA UBARI    / 0.5  /
< 
<   !C  These are for the CO2+H2O k-coefficients
< 
<   DATA WREFCO2 / 9.999999D-1, 9.99999D-1, 9.9999D-1, 9.999D-1,&
<        9.99D-1, 9.9D-1, 9.0D-1, 8.0D-1, 7.0D-1, 6.0D-1 /
< 
<   DATA WREFH2O / 1.0D-7, 1.0D-6, 1.0D-5, 1.0D-4, 1.0D-3, 1.0D-2,&
<        1.0D-1, 2.0D-1, 3.0D-1, 4.0D-1                  /
---
> !################ Nelli, April 07 ##########################
> !Setting up parameters needed by the correlated k lower
> !atmosphere radiation code
> 
> !     Number of atmospheric layers
>       integer, PARAMETER :: LL_LAYERS  = nAlts
> 
> !     Number of atmospheric levels:   2 * LL_LAYERS + 3
>       integer, PARAMETER :: LL_LEVELS  = 2*LL_LAYERS+3
> 
> !C     LL_NLAYRAD is the number of radiation code layers
> !C     LL_NLEVRAD is the number of radiation code levels.  Level N is the
> !C
>        integer, PARAMETER :: LL_NLAYRAD  = LL_LAYERS+1
>        integer, PARAMETER :: LL_NLEVRAD  = LL_LAYERS+2
> 
> ! Bottom layer subsurface temperature
>        real, parameter :: CoreTemp = 175.0
> 
> ! Surface and subsurface temperature constants
>        real, parameter :: Pa = 5.927E+7
>        real, parameter :: Pd = 88775.0
> 
> ! Stefan-Boltzmann constant in SI
>        real, parameter :: SBconstant = 5.67E-8
> 
> !      Ls variable
>        real :: ell_s = 0.0
> !C======================================================================C
> !C
> !C     RADINC.H    RADiation INCludes
> !C
> !C     Includes for the radiation code; RADIATION LAYERS, LEVELS,
> !C     number of spectral intervals. . .
> !C
> !C     GCM2.0  Feb 2003
> !C
> !C======================================================================C
> 
> !C     RADIATION parameters
> 
> !C     In radiation code, layer 1 corresponds to the stratosphere.  Level
> !C     1 is the top of the stratosphere.  The dummy layer is at the same
> !C     temperature as the (vertically isothermal) stratosphere, and
> !C     any time it is explicitly needed, the appropriate quantities will
> !C     be dealt with (aka "top". . .)
> 
> 
> !C
> !C     L_NSPECTI is the number of IR spectral intervals
> !C     L_NSPECTV is the number of Visual(or Solar) spectral intervals
> !C     L_NGAUSS  is the number of Gauss points for K-coefficients
> !C               GAUSS POINT 9 (aka the last one) is the special case
> !C     L_NWNGI   is L_NSPECTI*L_NGAUSS;  the total number of "intervals"
> !C               in the IR
> !C     L_NWNGV   is L_NSPECTV*L_NGAUSS;  the total number of "intervals"
> !C               in the VISUAL
> !C
> !C     L_NPREF   is the number of reference pressures that the
> !C               k-coefficients are calculated on
> !C     L_PINT    is the number of Lagrange interpolated reference
> !C               pressures for the CO2 k-coefficients.
> !C     L_NTREF   is the number of refernce temperatures for the
> !C               k-coefficients
> !C     L_TAUMAX  is the largest optical depth - larger ones are set
> !C               to this value.
> !C
> !C     L_REFH2O  The number of different water-mixing ratio values for
> !C               the k-coefficients that are now CO2+H2O.
> !C
> !C     L_NREFV   The spectral interval number of the visible reference
> !C               wavelength (i.e. the 0.67 micron band)
> !C
> !C
> !C     L_NLTE    The number of different LTE/NLTE heating rate conversion
> !C               factors
> !C
> !C----------------------------------------------------------------------C
> 
> 
>       integer, parameter :: L_NSPECTI =  5
>       integer, parameter :: L_NSPECTV =  7
>       integer, parameter :: L_NGAUSS  = 17
> 
>       integer, parameter :: L_NPREF   = 11
>       integer, parameter :: L_NTREF   =  7
>       integer, parameter :: L_TAUMAX  = 35
> 
>       integer, parameter :: L_PINT    = 51
> 
>       integer, parameter :: L_REFH2O  = 10
> 
>       integer, parameter :: L_NREFV   = 6
> 
>       integer, parameter :: L_NLTE    = 37
> 
> !C----------------------------------------------------------------------C
> !C
> !C                             radcommon.h
> !C                         FORTRAN PARAMETERS
> !C                          GCM2.0  Feb 2003
> !C
> !C----------------------------------------------------------------------C
> !C
> !C  "Include" grid.h and radinc.h before this file in code that uses
> !C  some or all of this common data set
> !C
> !C     WNOI       - Array of wavenumbers at the spectral interval
> !C                  centers for the infrared.  Array is NSPECTI
> !C                  elements long.
> !C     DWNI       - Array of "delta wavenumber", i.e., the width,
> !C                  in wavenumbers (cm^-1) of each IR spectral
> !C                  interval.  NSPECTI elements long.
> !C     WAVEI      - Array (NSPECTI elements long) of the wavelenght
> !C                  (in microns) at the center of each IR spectral
> !C                  interval.
> !C     WNOV       - Array of wavenumbers at the spectral interval
> !C                  center for the VISUAL.  Array is NSPECTV
> !C                  elements long.
> !C     DWNV       - Array of "delta wavenumber", i.e., the width,
> !C                  in wavenumbers (cm^-1) of each VISUAL spectral
> !C                  interval.  NSPECTV elements long.
> !C     WAVEV      - Array (NSPECTV elements long) of the wavelenght
> !C                  (in microns) at the center of each VISUAL spectral
> !C                  interval.
> !C     SOLARF     - Array (NSPECTV elements) of solar flux (W/M^2) in
> !C                  each spectral interval.  Values are for 1 AU, and
> !C                  are scaled to the Mars distance elsewhere.
> !C     SOL        - Solar flux at Mars (W/M^2)
> !C     TAURAY     - Array (NSPECTV elements) of the pressure-independent
> !C                  part of Rayleigh scattering optical depth.
> !C     PTOP       - Pressure at the top of the radiation code coordinate;
> !C                  = 0.5*Ptrop
> !C     FZEROI     - Fraction of zeros in the IR CO2 k-coefficients, for
> !C                  each temperature, pressure, and spectral interval
> !C     FZEROV     - Fraction of zeros in the VISUAL CO2 k-coefficients, for
> !C                  each temperature, pressure, and spectral interval
> !C
> !C     AEROSOL RADIATIVE OPTICAL CONSTANTS
> !C     Values are at the wavelenght interval center
> !C
> !C     MIE SCATTERING - Size distribution weighted
> !C     Qextv    - Extinction efficiency - in the visual.
> !C     QextREF  - Reference visual wavelength (.67 micron band)
> !C     Qscatv   - Scattering efficiency - in the visual.
> !C     WV       - Single scattering albedo - in the visual.
> !C     GV       - Asymmetry parameter - in the visual.
> !C
> !C     Qexti    - Extinction efficiency - in the infrared.
> !C     Qscati   - Scattering efficiency - in the infrared.
> !C     WI       - Single scattering albedo - in the infrared.
> !C     GI       - Asymmetry parameter - in the infrared.
> !C
> !C     VIS2IR   - VISIBLE (0.67 micron band) to IR (9 micron band) ratio.
> !C
> !C     XLTEFACTOR - correction factor for over-prediction of heating rates
> !C                  by the LTE code.  Occurs at pressures where NLTE processes
> !C                  become important
> !C
> !C     XLTEPRESSURE - pressure regime at which each factor is to be applied
> !C
> !C  "Include" grid.h and radinc.h before this file in code that uses
> !C  some or all of this common data set
> 
>       REAL :: WNOI(L_NSPECTI), DWNI(L_NSPECTI), WAVEI(L_NSPECTI)
>       REAL :: WNOV(L_NSPECTV), DWNV(L_NSPECTV), WAVEV(L_NSPECTV)
>       REAL :: SOLARF(L_NSPECTV), TAURAY(L_NSPECTV),SOL(L_NSPECTV)
> 
>       real :: CO2I(L_NTREF,L_PINT,L_REFH2O,L_NSPECTI,L_NGAUSS)
>       real :: CO2V(L_NTREF,L_PINT,L_REFH2O,L_NSPECTV,L_NGAUSS)
>       real :: FZEROI(L_NSPECTI)
>       real :: FZEROV(L_NSPECTV)
>       real :: PGASREF(L_NPREF), TGASREF(L_NTREF)
> 
>       real :: qextv(L_NSPECTV), qscatv(L_NSPECTV), wv(L_NSPECTV)
>       real :: gv(L_NSPECTV)
>       real :: QextREF, VIS2IR, Cmk, tlimit
> 
>       real :: qexti(L_NSPECTI), qscati(L_NSPECTI), wi(L_NSPECTI)
>       real :: gi(L_NSPECTI)
> 
>       real :: planckir(L_NSPECTI,8501)
> 
>       real :: PTOP, UBARI,  GWEIGHT(L_NGAUSS)
>       real :: PFGASREF(L_PINT)
> 
> !C  H2O and CO2 k-coefficients mixed
> 
>       real :: WREFCO2(L_REFH2O), WREFH2O(L_REFH2O)
> 
> !C  LTE/NLTE heating rate conversion factors and associated pressures
> 
>       real :: XLTEFACTOR(L_NLTE),XLTEPRESSURE(L_NLTE)
> 
> !  These are for the Gauss-split 0.95 case
> 
>       DATA GWEIGHT  / 4.8083554740D-02, 1.0563099137D-01,&
>                       1.4901065679D-01, 1.7227479710D-01,&
>                       1.7227479710D-01, 1.4901065679D-01,&
>                       1.0563099137D-01, 4.8083554740D-02,&
>                       2.5307134073D-03, 5.5595258613D-03,&
>                       7.8426661469D-03, 9.0670945845D-03,&
>                       9.0670945845D-03, 7.8426661469D-03,&
>                       5.5595258613D-03, 2.5307134073D-03,  0.0D0 /
> 
>       DATA UBARI    / 0.5  /
> 
> !C  These are for the CO2+H2O k-coefficients
> 
>       DATA WREFCO2 / 9.999999D-1, 9.99999D-1, 9.9999D-1, 9.999D-1,&
>                      9.99D-1, 9.9D-1, 9.0D-1, 8.0D-1, 7.0D-1, 6.0D-1 /
> 
>       DATA WREFH2O / 1.0D-7, 1.0D-6, 1.0D-5, 1.0D-4, 1.0D-3, 1.0D-2,&
>                      1.0D-1, 2.0D-1, 3.0D-1, 4.0D-1                  /
> 
>       DATA  Cmk     / 3.51E+22 /
> 
>       DATA  TLIMIT  / 0.5      /
> 
> ! ---------------------------------------------------------------------------
> !  Lopez-Valverde et al. (1998) Table 1, pg. 16809.
> !
> !      DATA  XLTEFACTOR / 1.007,  1.007,  1.007,  1.007,  1.008,&
> !                        1.009,  1.013,  1.017,  1.021,  1.027,&
> !                        1.033,  1.040,  1.047,  1.054,  1.061,&
> !                        1.069,  1.078,  1.087,  1.099,  1.112,&
> !                        1.129,  1.149,  1.174,  1.209,  1.263,&
> !                        1.344,  1.459,  1.608,  1.796,  2.033,&
> !                        2.343,  2.777,  3.404,  4.366,  5.856,&
> !                        8.161, 11.908, 17.791, 27.030, 40.897,&
> !                       60.798, 82.660,111.855,146.365,185.922,&
> !                      237.435,284.487,346.883,422.071,513.978,&
> !                      635.594 /
> !
> !      DATA XLTEPRESSURE / 1.122E-1, 8.822E-2, 6.912E-2, 5.397E-2, 4.202E-2,&
> !                         3.261E-2, 2.523E-2, 1.946E-2, 1.497E-2, 1.149E-2,&
> !                         8.793E-3, 6.713E-3, 5.115E-3, 3.890E-3, 2.953E-3,&
> !                         2.239E-3, 1.696E-3, 1.284E-3, 9.719E-4, 7.355E-4,&
> !                         5.566E-4, 4.212E-4, 3.187E-4, 2.412E-4, 1.825E-4,&
> !                         1.381E-4, 1.045E-4, 7.908E-5, 5.984E-5, 4.530E-5,&
> !                         3.430E-5, 2.600E-5, 1.973E-5, 1.500E-5, 1.145E-5,&
> !                         8.761E-6, 6.731E-6, 5.192E-6, 4.024E-6, 3.141E-6,&
> !                         2.478E-6, 1.981E-6, 1.601E-6, 1.306E-6, 1.076E-6,&
> !                         8.939E-7, 7.462E-7, 6.242E-7, 5.234E-7, 4.397E-7,&
> !                         3.702E-7 /
> !
> ! ---------------------------------------------------------------------------
> !  Bougher (2011) Calculated wrt QNLTE from MTGCM case for SZA = 0.0
> !  -- Ls = 90;  F10.7 = 70; dust = weak
> 
>       DATA  XLTEFACTOR / 1.00,  1.000,  1.000,  1.000,  1.000,&
>                         1.000,  1.000,  1.000,  1.000,  1.000,&
>                         1.000,  1.000,  1.000,  1.000,  1.000,&
>                         1.000,  1.000,  1.020,  1.040,  1.081,&
>                        1.13225, 1.13704, 1.2334, 1.32968, 2.47989,&
>                        5.17858, 9.21052, 13.9471, 19.1293, 25.9336,&
>                        45.9666, 113.783, 146.843, 212.959, 500.000,&
>                        1000.00, 1000.00 /
> !  mbar scale
>       DATA XLTEPRESSURE / 1.122E-1, 8.822E-2, 6.912E-2, 5.397E-2, 4.202E-2,&
>                          3.261E-2, 2.523E-2, 1.946E-2, 1.497E-2, 1.149E-2,&
>                          8.793E-3, 6.713E-3, 5.115E-3, 3.890E-3, 2.953E-3,&
>                          2.239E-3, 1.696E-3, 1.284E-3, 9.719E-4, 7.982E-4,&
>                          4.841E-4, 2.936E-4, 1.781E-4, 1.080E-4, 6.550E-5,&
>                          3.970E-5, 2.410E-5, 1.460E-5, 8.900E-6, 5.400E-6,&
>                          3.300E-6, 2.000E-6, 1.200E-6, 7.000E-7, 5.234E-7,&
>                          4.397E-7, 3.702E-7 /
> ! ---------------------------------------------------------------------------
> 
> !C  PLANCK defined variable
> 
> !      planckir
> 
> !C  SETSPI and SETSPV defined variables
> 
> !      WNOI, DWNI, WAVEI, WNOV, DWNV, WAVEV
> !      SOLARF, TAURAY, PTOP, TAUREF, GWEIGHT, UBARI
> !      PFGASREF
> 
> !C  SETRAD defined variables
> 
> !      CO2I, CO2V, PGASREF, TGASREF
> !      QEXTV, QSCATV, WV, GV
> !      QEXTI, QSCATI, WI, GI, QextREF, VIS2IR
> 
> !      fzeroi, fzerov
> 
> !      WREFCO2, WREFH2O
> 
> real*4 :: dummyalbedo(24,36), dummyti(24,36)
> 
> !##########################################################################
> !C
> !C     Includes for the nlte_tcool code: nlte_paramdef.h, nlte_commons.h
> !C
> !C     November 2017; S. W. Bougher
> !C
> !##########################################################################
> !
> !       Merging of different parameter definitions for new NLTE 15um param
> !       jul 2012    fgg+malv
> !       Modifications for MGITM
> !       nov 2017    swb
> !****************************************************************************
> ! *** Old mz1d.par ***
> ! Grids parameters :
> 
> 	integer,parameter :: nztabul = 79  ! # points in tabulation of Tesc & VC (ISO)
> 
> ! NLTE parameters:
> 
> 	integer,parameter :: nltot = 20	! incluye el actual # alt in NLTE module
> 	                                ! y el # alturas del Tstar110
> 	integer,parameter ::  nl = 12	! actual # alt in NLTE module & C.Matrix
> 	integer,parameter ::  nl2= nl-2 ! = nl-2, needed for matrix inversion (mmh2)
> 	integer,parameter ::  nzy = (nl-1)*4 + 1  ! Fine grid for C.Matrix
> 
> 	integer,parameter :: nl_cts = 2 + nltot-nl ! actual # alt para Tstar110
> 	integer,parameter ::  nzy_cts = (nl_cts-1)*4 + 1 ! fine grid for transmit calc
> 
> !  Other NLTE parameters:
> 	integer,parameter ::	nisot = 4	! number of isotopes considered
> 	integer,parameter :: 	nb = 41		! number of bands included
> 	integer,parameter :: 	nhist = 36	! # of temps in histogr.
> 	                                        ! (get it from histograms!)
> 	integer,parameter ::   nbox_max = 4     ! max.# boxes in histogram
> 
> ! *** Old tcr_15um.h ***
> 
>         integer,parameter ::  itt_cza = 13     ! Selection of NLTE scheme
>                                                ! Upper and lower limits of
>                                                ! NLTE model
> 
> !       Ptop_atm and Pbot_atm to Mars:  Domain of calculation applicability?
>         real, parameter  ::  Ptop_atm = 3.e-10, Pbottom_atm = 2.e-5
> 
>         real*8,parameter ::  rf19 = 1.d0, rf20 = 1.d0, rf21a = 1.d0, &
>                              rf21b = 1.d0, rf21c = 1.d0, rf33bc = 1.d0
> 
> ! *** Old bloque_dlvr11.f ***
> 
> ! data
> !        data nu(1,1),nu(1,2) /667.3801, 1335.1317/
> !        data nu(2,1)/662.3734/
> !        data nu(3,1)/648.4784/
> !        data nu(4,1)/664.7289/
>         real :: nu(nisot,8) = reshape( [ &
>              667.3801, 662.3734, 648.4784, 664.7289, &
>              1335.1317, 0., 0., 0., &
>              0., 0., 0., 0., &
>              0., 0., 0., 0., &
>              0., 0., 0., 0., &
>              0., 0., 0., 0., &
>              0., 0., 0., 0., &
>              0., 0., 0., 0. ], [nisot,8] )
> 
>         real,parameter ::  nu12_0200 = 1285.4087, nu12_1000 = 1388.1847
> 
>         integer :: indexisot(nisot)
>         data indexisot/26,28,36,27/
> 
> 	! ctes en el sistema cgs
> !       real*8, parameter ::  vlight = 2.9979245e10, ee = 1.43876866,&
> !                  hplanck = 6.6260755e-27, gamma = 1.191043934e-5
>         real*8, parameter ::  vlight = 2.9979245e10, ee = 1.43876866,&
>                    hplanck = 6.6260755e-27, gamma1 = 1.191043934e-5
> 
> 	! datos de marte
>         real ::  imr(nisot)
>         data imr /0.987, 0.00408, 0.0112, 0.000742/
> 
> !****************************************************************************
> !
> !       Merging of different common blocks used in the new NLTE 15um param
> !       jan 2012    fgg+malv
> !       Modification for MGITM usage in ModMars.f90
> !       nov 2017     swb
> !****************************************************************************
> ! *** Old datitos.cmn ***
> !       common /spectralv11/ elow, deltanu
> 	real :: elow(nisot,nb), deltanu(nisot,nb)
> 
> !       common/nu_levs_bands_v11/ nu11, nu12, nu121,
> !         &  	nu21, nu31, nu41
> 	real*8 :: nu11, nu12, nu121,nu21,nu31,nu41
> 
> !       common /aeinstein1v11/ a1_010_000, a1_020_010
> !       common /aeinstein2v11/ a2_010_000
> !       common /aeinstein3v11/ a3_010_000
> !       common /aeinstein4v11/ a4_010_000
> 	real*8 :: a1_010_000, a1_020_010, a2_010_000, &
>                   a3_010_000, a4_010_000
> 
> ! *** Old tabulation.cmn ***
> !common/input_tab_v11/ lnpnbtab,
> !          &  	tstar11tab, tstar21tab, tstar31tab, tstar41tab,
> !          &  	vc210tab, vc310tab, vc410tab
> 	real*8 :: lnpnbtab(nztabul),vc210tab(nztabul),vc310tab(nztabul), &
> 	           vc410tab(nztabul)
> 	real*8 :: tstar11tab(nztabul),tstar21tab(nztabul), &
>                tstar31tab(nztabul),tstar41tab(nztabul)
> 
> 
> ! *** Old nlte_results.cmn ***
> !       common/input_avilable_from/ input_cza
> 	integer :: input_cza
> 
> ! temperatura vibracional de entrada:
> !       common/temp626/ v626t1
> !       common/temp628/ v628t1
> !       common/temp636/ v636t1
> !       common/temp627/ v627t1
> 	real*8 :: v626t1(nl),v628t1(nl), v636t1(nl), v627t1(nl)
> 
> ! output de cza.for
> !       common /tv15um/	vt11, vt12, vt21, vt31, vt41
> 	real*8 ::  vt11(nl), vt12(nl), vt21(nl), vt31(nl), vt41(nl)
> !       common /hr15um/	hr110,hr210,hr310,hr410,hr121
> 	real*8 :: hr110(nl),hr121(nl),hr210(nl),hr310(nl),hr410(nl)
> !       common/sf15um/ el11,el12, el21, el31, el41
>         real*8 :: el11(nl), el12(nl),el21(nl), el31(nl), el41(nl)
> !       common/sl15um/ sl110,sl121, sl210,sl310,sl410
>         real*8 :: sl110(nl),sl121(nl),sl210(nl),sl310(nl),sl410(nl)
> 
> ! *** Old matrices.cmn ***
> 
> 
> ! curtis matrix de cza:
> !       common/curtis_matrixes_15um/ c110,c121, c210,
> !        &  	c310,c410,
> !        &  	vc110,vc121,vc210,vc310,vc410
> 	real*8 :: c110(nl,nl),c121(nl,nl),c210(nl,nl),c310(nl,nl), &
> 	          c410(nl,nl),vc110(nl),vc121(nl),vc210(nl), &
> 	          vc310(nl),vc410(nl)
> 
> ! for the cool-to-space formulation:
> !
> !      common/taustar_15um/ taustar11, taustar21, taustar31,
> !        &         taustar41, taustar12, taustar11_cts
> 	real*8 :: taustar11(nl), taustar21(nl), taustar31(nl), &
> 	      taustar41(nl), taustar12(nl), taustar11_cts(nl_cts)
> 
> ! *** Old atmref.cmn ***
> 
> 
> ! *******************************************************************
> ! NLTE Subgrid
> ! Revision for conflicts with the rest of the code via ModMars.f90
> !
> !       common /atm_nl/ zl, t, pl, nt, co2, n2, co, o3p,
> !          &     co2vmr, n2vmr, covmr, o3pvmr,
> !          &     hrkday_factor
> !       real, dimension (nl) ::  zl,t,pl,nt,co2, n2, co, o3p, &
> !         co2vmr, n2vmr, covmr, o3pvmr,hrkday_factor
> !       real, dimension (nl) ::  zl,tl,pl,nt,co2, n2, co, o3p, &
> !         co2vmr, n2vmr, covmr, o3pvmr,hrkday_factor
>         real, dimension (nl) ::  zl,tl,pl,ntl,co2, n2, co, o3p, &
>           co2vmr, n2vmr, covmr, o3pvmr,hrkday_factor
> ! *******************************************************************
> 
> ! Subgrid Transmittances
> !
> !       common /atm_ny/ zy, ty, py, nty, co2y
>         real :: zy(nzy), ty(nzy), py(nzy), nty(nzy), co2y(nzy)
> 
> ! Grids and indexes
> !       common/deltazetas/ deltaz, deltazy, deltaz_cts, deltazy_cts,
> !         &    jlowerboundary, jtopboundary, jtopCTS
> 	real  ::    deltaz, deltazy, deltaz_cts, deltazy_cts
> 	integer :: jlowerboundary, jtopboundary, jtopCTS
> 
> ! NLTE-CTS Subgrid
> !       common /atm_nl_cts/ zl_cts, t_cts, pl_cts, nt_cts,
> !        &    co2_cts, n2_cts, co_cts, o3p_cts,
> !        &    co2vmr_cts, n2vmr_cts, covmr_cts, o3pvmr_cts,
> !        &    hrkday_factor_cts,mmean_cts,cpnew_cts
> 
>         real :: zl_cts(nl_cts), t_cts(nl_cts), pl_cts(nl_cts), &
>           nt_cts(nl_cts), co2_cts(nl_cts), &
>           n2_cts(nl_cts), co_cts(nl_cts), &
>           o3p_cts(nl_cts), co2vmr_cts(nl_cts), n2vmr_cts(nl_cts), &
>           covmr_cts(nl_cts), o3pvmr_cts(nl_cts), &
>           hrkday_factor_cts(nl_cts),mmean_cts(nl_cts), &
>           cpnew_cts(nl_cts)
> 
> ! CTS Subgrid Transmittances
> !       common /atm_ny_cts/ zy_cts, ty_cts, py_cts, nty_cts, co2y_cts
>         real :: zy_cts(nzy_cts), ty_cts(nzy_cts), py_cts(nzy_cts), &
>                 nty_cts(nzy_cts), co2y_cts(nzy_cts)
> 
> ! *** Old rates.cmn ***
> !common/rates_vt/
> !       &      k19ba(4),k19bb(4),k19bc(4), k19bap(4),k19bbp(4),k19bcp(4),
> !       &      k19ca(4),k19cb(4),k19cc(4), k19cap(4),k19cbp(4),k19ccp(4),
> !       &      k20b(4),k20c(4), k20bp(4),k20cp(4)
> 	real*8 :: k19ba(4),k19bb(4),k19bc(4), k19bap(4),k19bbp(4),k19bcp(4)
> 	real*8 :: k19ca(4),k19cb(4),k19cc(4), k19cap(4),k19cbp(4),k19ccp(4)
> 	real*8 :: k20b(4),k20c(4), k20bp(4),k20cp(4)
> 
> !      common/rates_vv/
> !       &      	k21b(4),k21c(4), k21bp(4),k21cp(4),
> !       &      	k33c, k33cp(2:4)
> 	real*8 :: k21b(4),k21c(4), k21bp(4),k21cp(4)
> 	real*8 :: k33c, k33cp(2:4)
> 
> !      common/rates_last/ k23k21c, k24k21c, k34k21c, &
> !              	k23k21cp, k24k21cp, k34k21cp
> 	real*8 :: k23k21c,k24k21c,k34k21c, k23k21cp,k24k21cp,k34k21cp
> 
> ! *** Old curtis.cmn ***
> !       common /ini_file/ ibcode1
> 	character :: ibcode1*1
> !       common/block1/ alsa,alda,ka,kr
> 	real*8 :: ka(nbox_max),alsa(nbox_max),alda(nbox_max)
> 	integer :: kr
> !       common/block2/ hisfile
> 	character :: hisfile*75
> !       common/block3/ pp,ta,w
> 	real*8 :: pp,ta(nbox_max),w
> !       common/block4/ no,sk1,xls1,xld1,thist,nbox
> 	integer :: nbox
> 	real*8	:: sk1(nhist,nbox_max), xls1(nhist,nbox_max), &
> 	           xld1(nhist,nbox_max),thist(nhist), &
> 	           no(nbox_max)
> !       common/block5/eqw, aa,  cc, dd, ddbox, ccbox, mr, mr_cts
> 	real*8 :: eqw, aa, cc, dd
> 	real*8 :: ddbox(nbox_max), ccbox(nbox_max)
> 	real*8  :: mr(nzy), mr_cts(nzy_cts)
> !       common/blockstore/no_stored, sk1_stored, xls1_stored,
> !        &          xld1_stored, thist_stored, nbox_stored,
> !        &          mm_stored
>          real*8 :: sk1_stored(nb,nhist,nbox_max)
> 	 real*8 :: xls1_stored(nb,nhist,nbox_max)
> 	 real*8 :: xld1_stored(nb,nhist,nbox_max)
> 	 real*8 :: thist_stored(nb,nhist)
> 	 real*8 :: no_stored(nb,nbox_max)
> 	 integer :: nbox_stored(nb), mm_stored(nb)
572,895c895
<   DATA  Cmk     / 3.51E+22 /
< 
<   DATA  TLIMIT  / 0.5      /
< 
<   ! ---------------------------------------------------------------------------
<   !  Lopez-Valverde et al. (1998) Table 1, pg. 16809.
<   !
<   !      DATA  XLTEFACTOR / 1.007,  1.007,  1.007,  1.007,  1.008,&
<   !                        1.009,  1.013,  1.017,  1.021,  1.027,&
<   !                        1.033,  1.040,  1.047,  1.054,  1.061,&
<   !                        1.069,  1.078,  1.087,  1.099,  1.112,&
<   !                        1.129,  1.149,  1.174,  1.209,  1.263,&
<   !                        1.344,  1.459,  1.608,  1.796,  2.033,&
<   !                        2.343,  2.777,  3.404,  4.366,  5.856,&
<   !                        8.161, 11.908, 17.791, 27.030, 40.897,&
<   !                       60.798, 82.660,111.855,146.365,185.922,&
<   !                      237.435,284.487,346.883,422.071,513.978,&
<   !                      635.594 /
<   !
<   !      DATA XLTEPRESSURE / 1.122E-1, 8.822E-2, 6.912E-2, 5.397E-2, 4.202E-2,&
<   !                         3.261E-2, 2.523E-2, 1.946E-2, 1.497E-2, 1.149E-2,&
<   !                         8.793E-3, 6.713E-3, 5.115E-3, 3.890E-3, 2.953E-3,&
<   !                         2.239E-3, 1.696E-3, 1.284E-3, 9.719E-4, 7.355E-4,&
<   !                         5.566E-4, 4.212E-4, 3.187E-4, 2.412E-4, 1.825E-4,&
<   !                         1.381E-4, 1.045E-4, 7.908E-5, 5.984E-5, 4.530E-5,&
<   !                         3.430E-5, 2.600E-5, 1.973E-5, 1.500E-5, 1.145E-5,&
<   !                         8.761E-6, 6.731E-6, 5.192E-6, 4.024E-6, 3.141E-6,&
<   !                         2.478E-6, 1.981E-6, 1.601E-6, 1.306E-6, 1.076E-6,&
<   !                         8.939E-7, 7.462E-7, 6.242E-7, 5.234E-7, 4.397E-7,&
<   !                         3.702E-7 /
<   !
<   ! ---------------------------------------------------------------------------
<   !  Bougher (2011) Calculated wrt QNLTE from MTGCM case for SZA = 0.0
<   !  -- Ls = 90;  F10.7 = 70; dust = weak
< 
<   DATA  XLTEFACTOR / 1.00,  1.000,  1.000,  1.000,  1.000,&
<        1.000,  1.000,  1.000,  1.000,  1.000,&
<        1.000,  1.000,  1.000,  1.000,  1.000,&
<        1.000,  1.000,  1.020,  1.040,  1.081,&
<        1.13225, 1.13704, 1.2334, 1.32968, 2.47989,&
<        5.17858, 9.21052, 13.9471, 19.1293, 25.9336,&
<        45.9666, 113.783, 146.843, 212.959, 500.000,&
<        1000.00, 1000.00 /
<   !  mbar scale
<   DATA XLTEPRESSURE / 1.122E-1, 8.822E-2, 6.912E-2, 5.397E-2, 4.202E-2,&
<        3.261E-2, 2.523E-2, 1.946E-2, 1.497E-2, 1.149E-2,&
<        8.793E-3, 6.713E-3, 5.115E-3, 3.890E-3, 2.953E-3,&
<        2.239E-3, 1.696E-3, 1.284E-3, 9.719E-4, 7.982E-4,&
<        4.841E-4, 2.936E-4, 1.781E-4, 1.080E-4, 6.550E-5,&
<        3.970E-5, 2.410E-5, 1.460E-5, 8.900E-6, 5.400E-6,&
<        3.300E-6, 2.000E-6, 1.200E-6, 7.000E-7, 5.234E-7,&
<        4.397E-7, 3.702E-7 /
<   ! ---------------------------------------------------------------------------
< 
<   !C  PLANCK defined variable
< 
<   !      planckir
< 
<   !C  SETSPI and SETSPV defined variables
< 
<   !      WNOI, DWNI, WAVEI, WNOV, DWNV, WAVEV
<   !      SOLARF, TAURAY, PTOP, TAUREF, GWEIGHT, UBARI
<   !      PFGASREF
< 
<   !C  SETRAD defined variables
< 
<   !      CO2I, CO2V, PGASREF, TGASREF
<   !      QEXTV, QSCATV, WV, GV
<   !      QEXTI, QSCATI, WI, GI, QextREF, VIS2IR
< 
<   !      fzeroi, fzerov
< 
<   !      WREFCO2, WREFH2O
< 
<   real*4 :: dummyalbedo(24,36), dummyti(24,36)
< 
<   !##########################################################################
<   !C
<   !C     Includes for the nlte_tcool code: nlte_paramdef.h, nlte_commons.h
<   !C
<   !C     November 2017; S. W. Bougher
<   !C
<   !##########################################################################
<   !
<   !       Merging of different parameter definitions for new NLTE 15um param
<   !       jul 2012    fgg+malv
<   !       Modifications for MGITM
<   !       nov 2017    swb
<   !****************************************************************************
<   ! *** Old mz1d.par ***
<   ! Grids parameters :
< 
<   integer,parameter :: nztabul = 79  ! # points in tabulation of Tesc & VC (ISO)
< 
<   ! NLTE parameters:
< 
<   integer,parameter :: nltot = 20	! incluye el actual # alt in NLTE module
<   ! y el # alturas del Tstar110
<   integer,parameter ::  nl = 12	! actual # alt in NLTE module & C.Matrix
<   integer,parameter ::  nl2= nl-2 ! = nl-2, needed for matrix inversion (mmh2)
<   integer,parameter ::  nzy = (nl-1)*4 + 1  ! Fine grid for C.Matrix
< 
<   integer,parameter :: nl_cts = 2 + nltot-nl ! actual # alt para Tstar110
<   integer,parameter ::  nzy_cts = (nl_cts-1)*4 + 1 ! fine grid for transmit calc
< 
<   !  Other NLTE parameters:
<   integer,parameter ::	nisot = 4	! number of isotopes considered
<   integer,parameter :: 	nb = 41		! number of bands included
<   integer,parameter :: 	nhist = 36	! # of temps in histogr.
<   ! (get it from histograms!)
<   integer,parameter ::   nbox_max = 4     ! max.# boxes in histogram
< 
<   ! *** Old tcr_15um.h ***
< 
<   integer,parameter ::  itt_cza = 13     ! Selection of NLTE scheme
<   ! Upper and lower limits of
<   ! NLTE model
< 
<   !       Ptop_atm and Pbot_atm to Mars:  Domain of calculation applicability?
<   real, parameter  ::  Ptop_atm = 3.e-10, Pbottom_atm = 2.e-5
< 
<   real*8,parameter ::  rf19 = 1.d0, rf20 = 1.d0, rf21a = 1.d0, &
<        rf21b = 1.d0, rf21c = 1.d0, rf33bc = 1.d0
< 
<   ! *** Old bloque_dlvr11.f ***
< 
<   ! data
<   !        data nu(1,1),nu(1,2) /667.3801, 1335.1317/
<   !        data nu(2,1)/662.3734/
<   !        data nu(3,1)/648.4784/
<   !        data nu(4,1)/664.7289/
<   real :: nu(nisot,8) = reshape( [ &
<        667.3801, 662.3734, 648.4784, 664.7289, &
<        1335.1317, 0., 0., 0., &
<        0., 0., 0., 0., &
<        0., 0., 0., 0., &
<        0., 0., 0., 0., &
<        0., 0., 0., 0., &
<        0., 0., 0., 0., &
<        0., 0., 0., 0. ], [nisot,8] )
< 
<   real,parameter ::  nu12_0200 = 1285.4087, nu12_1000 = 1388.1847
< 
<   integer :: indexisot(nisot)
<   data indexisot/26,28,36,27/
< 
<   ! ctes en el sistema cgs
<   !       real*8, parameter ::  vlight = 2.9979245e10, ee = 1.43876866,&
<   !                  hplanck = 6.6260755e-27, gamma = 1.191043934e-5
<   real*8, parameter ::  vlight = 2.9979245e10, ee = 1.43876866,&
<        hplanck = 6.6260755e-27, gamma1 = 1.191043934e-5
< 
<   ! datos de marte
<   real ::  imr(nisot)
<   data imr /0.987, 0.00408, 0.0112, 0.000742/
< 
<   !****************************************************************************
<   !
<   !       Merging of different common blocks used in the new NLTE 15um param
<   !       jan 2012    fgg+malv
<   !       Modification for MGITM usage in ModMars.f90
<   !       nov 2017     swb
<   !****************************************************************************
<   ! *** Old datitos.cmn ***
<   !       common /spectralv11/ elow, deltanu
<   real :: elow(nisot,nb), deltanu(nisot,nb)
< 
<   !       common/nu_levs_bands_v11/ nu11, nu12, nu121,
<   !         &  	nu21, nu31, nu41
<   real*8 :: nu11, nu12, nu121,nu21,nu31,nu41
< 
<   !       common /aeinstein1v11/ a1_010_000, a1_020_010
<   !       common /aeinstein2v11/ a2_010_000
<   !       common /aeinstein3v11/ a3_010_000
<   !       common /aeinstein4v11/ a4_010_000
<   real*8 :: a1_010_000, a1_020_010, a2_010_000, &
<        a3_010_000, a4_010_000
< 
<   ! *** Old tabulation.cmn ***
<   !common/input_tab_v11/ lnpnbtab,
<   !          &  	tstar11tab, tstar21tab, tstar31tab, tstar41tab,
<   !          &  	vc210tab, vc310tab, vc410tab
<   real*8 :: lnpnbtab(nztabul),vc210tab(nztabul),vc310tab(nztabul), &
<        vc410tab(nztabul)
<   real*8 :: tstar11tab(nztabul),tstar21tab(nztabul), &
<        tstar31tab(nztabul),tstar41tab(nztabul)
< 
< 
<   ! *** Old nlte_results.cmn ***
<   !       common/input_avilable_from/ input_cza
<   integer :: input_cza
< 
<   ! temperatura vibracional de entrada:
<   !       common/temp626/ v626t1
<   !       common/temp628/ v628t1
<   !       common/temp636/ v636t1
<   !       common/temp627/ v627t1
<   real*8 :: v626t1(nl),v628t1(nl), v636t1(nl), v627t1(nl)
< 
<   ! output de cza.for
<   !       common /tv15um/	vt11, vt12, vt21, vt31, vt41
<   real*8 ::  vt11(nl), vt12(nl), vt21(nl), vt31(nl), vt41(nl)
<   !       common /hr15um/	hr110,hr210,hr310,hr410,hr121
<   real*8 :: hr110(nl),hr121(nl),hr210(nl),hr310(nl),hr410(nl)
<   !       common/sf15um/ el11,el12, el21, el31, el41
<   real*8 :: el11(nl), el12(nl),el21(nl), el31(nl), el41(nl)
<   !       common/sl15um/ sl110,sl121, sl210,sl310,sl410
<   real*8 :: sl110(nl),sl121(nl),sl210(nl),sl310(nl),sl410(nl)
< 
<   ! *** Old matrices.cmn ***
< 
< 
<   ! curtis matrix de cza:
<   !       common/curtis_matrixes_15um/ c110,c121, c210,
<   !        &  	c310,c410,
<   !        &  	vc110,vc121,vc210,vc310,vc410
<   real*8 :: c110(nl,nl),c121(nl,nl),c210(nl,nl),c310(nl,nl), &
<        c410(nl,nl),vc110(nl),vc121(nl),vc210(nl), &
<        vc310(nl),vc410(nl)
< 
<   ! for the cool-to-space formulation:
<   !
<   !      common/taustar_15um/ taustar11, taustar21, taustar31,
<   !        &         taustar41, taustar12, taustar11_cts
<   real*8 :: taustar11(nl), taustar21(nl), taustar31(nl), &
<        taustar41(nl), taustar12(nl), taustar11_cts(nl_cts)
< 
<   ! *** Old atmref.cmn ***
< 
< 
<   ! *******************************************************************
<   ! NLTE Subgrid
<   ! Revision for conflicts with the rest of the code via ModMars.f90
<   !
<   !       common /atm_nl/ zl, t, pl, nt, co2, n2, co, o3p,
<   !          &     co2vmr, n2vmr, covmr, o3pvmr,
<   !          &     hrkday_factor
<   !       real, dimension (nl) ::  zl,t,pl,nt,co2, n2, co, o3p, &
<   !         co2vmr, n2vmr, covmr, o3pvmr,hrkday_factor
<   !       real, dimension (nl) ::  zl,tl,pl,nt,co2, n2, co, o3p, &
<   !         co2vmr, n2vmr, covmr, o3pvmr,hrkday_factor
<   real, dimension (nl) ::  zl,tl,pl,ntl,co2, n2, co, o3p, &
<        co2vmr, n2vmr, covmr, o3pvmr,hrkday_factor
<   ! *******************************************************************
< 
<   ! Subgrid Transmittances
<   !
<   !       common /atm_ny/ zy, ty, py, nty, co2y
<   real :: zy(nzy), ty(nzy), py(nzy), nty(nzy), co2y(nzy)
< 
<   ! Grids and indexes
<   !       common/deltazetas/ deltaz, deltazy, deltaz_cts, deltazy_cts,
<   !         &    jlowerboundary, jtopboundary, jtopCTS
<   real  ::    deltaz, deltazy, deltaz_cts, deltazy_cts
<   integer :: jlowerboundary, jtopboundary, jtopCTS
< 
<   ! NLTE-CTS Subgrid
<   !       common /atm_nl_cts/ zl_cts, t_cts, pl_cts, nt_cts,
<   !        &    co2_cts, n2_cts, co_cts, o3p_cts,
<   !        &    co2vmr_cts, n2vmr_cts, covmr_cts, o3pvmr_cts,
<   !        &    hrkday_factor_cts,mmean_cts,cpnew_cts
< 
<   real :: zl_cts(nl_cts), t_cts(nl_cts), pl_cts(nl_cts), &
<        nt_cts(nl_cts), co2_cts(nl_cts), &
<        n2_cts(nl_cts), co_cts(nl_cts), &
<        o3p_cts(nl_cts), co2vmr_cts(nl_cts), n2vmr_cts(nl_cts), &
<        covmr_cts(nl_cts), o3pvmr_cts(nl_cts), &
<        hrkday_factor_cts(nl_cts),mmean_cts(nl_cts), &
<        cpnew_cts(nl_cts)
< 
<   ! CTS Subgrid Transmittances
<   !       common /atm_ny_cts/ zy_cts, ty_cts, py_cts, nty_cts, co2y_cts
<   real :: zy_cts(nzy_cts), ty_cts(nzy_cts), py_cts(nzy_cts), &
<        nty_cts(nzy_cts), co2y_cts(nzy_cts)
< 
<   ! *** Old rates.cmn ***
<   !common/rates_vt/
<   !       &      k19ba(4),k19bb(4),k19bc(4), k19bap(4),k19bbp(4),k19bcp(4),
<   !       &      k19ca(4),k19cb(4),k19cc(4), k19cap(4),k19cbp(4),k19ccp(4),
<   !       &      k20b(4),k20c(4), k20bp(4),k20cp(4)
<   real*8 :: k19ba(4),k19bb(4),k19bc(4), k19bap(4),k19bbp(4),k19bcp(4)
<   real*8 :: k19ca(4),k19cb(4),k19cc(4), k19cap(4),k19cbp(4),k19ccp(4)
<   real*8 :: k20b(4),k20c(4), k20bp(4),k20cp(4)
< 
<   !      common/rates_vv/
<   !       &      	k21b(4),k21c(4), k21bp(4),k21cp(4),
<   !       &      	k33c, k33cp(2:4)
<   real*8 :: k21b(4),k21c(4), k21bp(4),k21cp(4)
<   real*8 :: k33c, k33cp(2:4)
< 
<   !      common/rates_last/ k23k21c, k24k21c, k34k21c, &
<   !              	k23k21cp, k24k21cp, k34k21cp
<   real*8 :: k23k21c,k24k21c,k34k21c, k23k21cp,k24k21cp,k34k21cp
< 
<   ! *** Old curtis.cmn ***
<   !       common /ini_file/ ibcode1
<   character :: ibcode1*1
<   !       common/block1/ alsa,alda,ka,kr
<   real*8 :: ka(nbox_max),alsa(nbox_max),alda(nbox_max)
<   integer :: kr
<   !       common/block2/ hisfile
<   character :: hisfile*75
<   !       common/block3/ pp,ta,w
<   real*8 :: pp,ta(nbox_max),w
<   !       common/block4/ no,sk1,xls1,xld1,thist,nbox
<   integer :: nbox
<   real*8	:: sk1(nhist,nbox_max), xls1(nhist,nbox_max), &
<        xld1(nhist,nbox_max),thist(nhist), &
<        no(nbox_max)
<   !       common/block5/eqw, aa,  cc, dd, ddbox, ccbox, mr, mr_cts
<   real*8 :: eqw, aa, cc, dd
<   real*8 :: ddbox(nbox_max), ccbox(nbox_max)
<   real*8  :: mr(nzy), mr_cts(nzy_cts)
<   !       common/blockstore/no_stored, sk1_stored, xls1_stored,
<   !        &          xld1_stored, thist_stored, nbox_stored,
<   !        &          mm_stored
<   real*8 :: sk1_stored(nb,nhist,nbox_max)
<   real*8 :: xls1_stored(nb,nhist,nbox_max)
<   real*8 :: xld1_stored(nb,nhist,nbox_max)
<   real*8 :: thist_stored(nb,nhist)
<   real*8 :: no_stored(nb,nbox_max)
<   integer :: nbox_stored(nb), mm_stored(nb)
< 
<   !****************************************************************************
---
> !****************************************************************************
921c921
< 
---
>  
1006c1006
<             
---
> 
1033,1303c1033,1305
<   !################ Nelli, April 07 ##########################
<   !Filling arrays needed by the correlated k lower
<   !atmosphere radiation code
< 
<   subroutine radsetup
< 
<     !C  GCM2.0  Feb 2003
<     !C
<     !C     PURPOSE:
<     !C        Bundle the new radiation code setup subroutines and call
<     !C     this one subroutine from main, where the three included files
<     !C     are also listed.  Quantities are passed between this driver
<     !C     and the radiation code via common (in radcommon.h).
<     !C
<     !C----------------------------------------------------------------------C
< 
<     implicit none
< 
< 
<     REAL :: FACTOR
<     INTEGER :: NW
< 
<     !C======================================================================C
< 
<     call setspv
<     call setspi
<     call setrad
< 
<     !C  Scale IR opacities (Qexti and Qscati) such that
<     !C  TAU(0.67 micron)/TAU(9 micron) = VIS2IR, which nominally is 2.
< 
<     QextREF = Qextv(L_NREFV)
<     VIS2IR  = 2.75D0
< 
<     factor  = Qextv(6)/(VIS2IR*Qexti(4))
< 
<     DO NW=1,L_NSPECTI
<        Qexti(NW)  = Qexti(NW)*factor
<        Qscati(NW) = Qscati(NW)*factor
<     END DO
< 
<     PTOP = 10.0**PFGASREF(1)
< 
<   end subroutine radsetup
< 
<   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
<   subroutine SETSPV
< 
<     !C  GCM2.0  Feb 2003
<     !C
<     !C     PURPOSE:
<     !C        Set up the spectral intervals in the Visual (solar).  Based
<     !C     on Chris McKay's SETSPV code.
<     !C
<     !C     AUTHOR
<     !C        Jim Schaeffer
<     !C
<     !C     UPDATES FOR
<     !C        Bob Haberle
<     !C
<     !C     ENVIRONMENT
<     !C        Sun ULTRA-2       SOLARIS 2.5.1      FORTRAN
<     !C
<     !C     REVISION HISTORY
<     !C        Original 9/17/1998
<     !C     VERSION 2.0  OCT 2001
<     !C
<     !C     INPUT PARAMETERS
<     !C     L_NSPECTV  - Number of spectral intervals in the Visual
<     !C
<     !C     OUTPUT PARAMETERS
<     !C     WNOV       - Array of wavenumbers at the spectral interval
<     !C                  center for the VISUAL.  Array is NSPECTV
<     !C                  elements long.
<     !C     DWNV       - Array of "delta wavenumber", i.e., the width,
<     !C                  in wavenumbers (cm^-1) of each VISUAL spectral
<     !C                  interval.  NSPECTV elements long.
<     !C     WAVEV      - Array (NSPECTV elements long) of the wavelenght
<     !C                  (in microns) at the center of each VISUAL spectral
<     !C                  interval.
<     !C     SOLARF     - Array (NSPECTV elements) of solar flux (W/M^2) in
<     !C                  each spectral interval.  Values are for 1 AU, and
<     !C                  are scaled to the Mars distance elsewhere.
<     !C     TAURAY     - Array (NSPECTV elements) of the wavelength dependent
<     !C                  part of Rayleigh Scattering.  The pressure dependent
<     !C                  part is computed elsewhere (OPTCV).
<     !C     CALLED BY
<     !C        RADIATION
<     !C
<     !C     SUBROUTINES CALLED
<     !C        NONE
<     !C
<     !C**********************************************************************C
< 
<     implicit none
< 
< 
<     !C     BWNV - Bin wavenumber of the edges of the VISUAL spectral bins
<     !C     units are inverse centimeters.  Dimension needs to be changed
<     !C     if the number of VISUAL bins changes.
< 
<     REAL :: BWNV(L_NSPECTV+1)
<     REAL :: SOLAR(L_NSPECTV)
< 
<     REAL ::  P0, GRAV, SCALEP, SUM, WL
<     INTEGER :: N, M
< 
<     !C     P0      - Rayleigh scattering reference pressure in pascals.
<     !C     GRAV    - Acceleration due to gravity (g) - MKS
<     !C     SCALEP  - multiply by 100 to convert pressure from millibars
<     !C               to pascals.
< 
<     DATA P0     / 9.423D+6 /
<     DATA GRAV   / 3.72     /
<     DATA SCALEP / 100.0    /
< 
<     !C     Bin wavenumber - wavenumber [cm^(-1)] at the edges of the VISUAL
<     !C     spectral bins.  Go from smaller to larger wavenumbers, the same as
<     !C     in the IR.
< 
<     DATA BWNV / 2222.22D0,  3087.37D0,  4030.63D0,  5370.57D0,&
<          7651.11D0, 12500.00D0, 25000.00D0, 41666.67D0 /
< 
<     !C     Solar flux within each spectral interval, at 1AU (W/M^2)
<     !C     Sum equals 1356 W/m^2 (values from Allen, 4th edition)
< 
<     DATA SOLAR /  17.0, 29.0, 52.0, 148.0, 348.0, 643.0, 118.0 /
< 
<     !C======================================================================C
< 
<     !C     Set up mean wavenumbers and wavenumber deltas.  Units of
<     !C     wavenumbers is cm^(-1); units of wavelengths is microns.
< 
<     do M=1,L_NSPECTV
<        WNOV(M)  = 0.5*(BWNV(M+1)+BWNV(M))
<        DWNV(M)  = BWNV(M+1)-BWNV(M)
<        WAVEV(M) = 1.0E+4/WNOV(M)
<     end do
< 
<     !C     Sum the solar flux, and write out the result.
< 
<     sum = 0.0
<     do N=1,L_NSPECTV
<        SOLARF(N) = SOLAR(N)
<        sum       = sum+SOLARF(N)
<     end do
<     write(6,'("Solar flux at 1AU = ",f7.2," W/M^2")') sum
< 
<     !C     Set up the wavelength dependent part of Rayleigh Scattering.
<     !C     The pressure dependent part will be computed elsewhere (OPTCV).
<     !C     WAVEV is in microns.  There is no Rayleigh scattering in the IR.
< 
<     do N=1,L_NSPECTV
<        WL        = WAVEV(N)
<        TAURAY(N) = (8.7/grav)*(1.527*(1.0+0.013/wl**2)/wl**4)*&
<             scalep/P0
<     end do
< 
< 
<   end subroutine SETSPV
< 
<   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
< 
<   subroutine setspi
< 
<     !C  GCM2.0  Feb 2003
<     !C
<     !C     PURPOSE:
<     !C        Set up the spectral intervals in the infrared.  Based on
<     !C     Chris McKay's SETSPI code.
<     !C
<     !C     AUTHOR
<     !C        Jim Schaeffer
<     !C
<     !C     UPDATES FOR
<     !C        Bob Haberle
<     !C
<     !C     ENVIRONMENT
<     !C        Sun ULTRA-2       SOLARIS 2.5.1      FORTRAN
<     !C
<     !C     REVISION HISTORY
<     !C        Original 9/17/1998
<     !C        VERSION 2  OCT 2001
<     !C
<     !C     INPUT PARAMETERS
<     !C     L_NSPECTI  - Number of spectral intervals in the INFRARED
<     !C
<     !C     OUTPUT PARAMETERS
<     !C     WNOI       - Array of wavenumbers at the spectral interval
<     !C                  centers for the infrared.  Array is NSPECTI
<     !C                  elements long.
<     !C     DWNI       - Array of "delta wavenumber", i.e., the width,
<     !C                  in wavenumbers (cm^-1) of each IR spectral
<     !C                  interval.  NSPECTI elements long.
<     !C     WAVEI      - Array (NSPECTI elements long) of the wavelenght
<     !C                  (in microns) at the center of each IR spectral
<     !C                  interval.
<     !C
<     !C     CALLED BY
<     !C        RADIATION
<     !C
<     !C     SUBROUTINES CALLED
<     !C        NONE
<     !C
<     !C**********************************************************************C
< 
<     implicit none
< 
< 
<     !C     BWNI - Bin wavenumber of the edges of the IR spectral bins
<     !C     units are inverse centimeters.  Dimension needs to be changed
<     !C     if the number of IR bins changes.
< 
<     REAL :: BWNI(L_NSPECTI+1)
< 
<     real :: a, b, x(12), w(12), ans, y, bpa, bma, T
<     real :: c1, c2, wn1, wn2
<     integer :: n, nw, nt, m
< 
<     !C  C1 and C2 values from Goody and Yung (2nd edition)  MKS units
<     !C  These values lead to a "sigma" (sigma*T^4) of 5.67032E-8 W m^-2 K^-4
< 
<     data c1 / 3.741832D-16 /     ! W m^-2
<     data c2 / 1.438786D-2  /     ! m K
< 
<     data x / -0.981560634246719D0,  -0.904117256370475D0,&
<          -0.769902674194305D0,  -0.587317954286617D0,&
<          -0.367831498998180D0,  -0.125233408511469D0,&
<          0.125233408511469D0,   0.367831498998180D0,&
<          0.587317954286617D0,   0.769902674194305D0,&
<          0.904117256370475D0,   0.981560634246719D0    /
< 
<     data w /  0.047175336386512D0,   0.106939325995318D0,&
<          0.160078328543346D0,   0.203167426723066D0,&
<          0.233492536538355D0,   0.249147045813403D0,&
<          0.249147045813403D0,   0.233492536538355D0,&
<          0.203167426723066D0,   0.160078328543346D0,&
<          0.106939325995318D0,   0.047175336386512D0   /
< 
<     !C======================================================================C
< 
<     !C     Bin wavenumber - wavenumber [cm^(-1)] at the edges of the IR
<     !C     spectral bins.
< 
<     BWNI( 1) =   10.000D0
<     BWNI( 2) =  166.667D0
<     BWNI( 3) =  416.667D0
<     BWNI( 4) =  833.333D0
<     BWNI( 5) = 1250.000D0
<     BWNI( 6) = 2500.000D0
< 
<     !C     Set up mean wavenumbers and wavenumber deltas.  Units of
<     !C     wavenumbers is cm^(-1); units of wavelengths is microns.
< 
<     do M=1,L_NSPECTI
<        WNOI(M)  = 0.5*(BWNI(M+1)+BWNI(M))
<        DWNI(M)  = BWNI(M+1)-BWNI(M)
<        WAVEI(M) = 1.0E+4/WNOI(M)
<     end do
< 
<     !C  For each IR wavelength interval, compute the integral of B(T), the
<     !C  Planck function, divided by the wavelength interval, in cm-1.  The
<     !C  integration is in MKS units, the final answer is the same as the
<     !C  original planck.f; W m^-2 wavenumber^-1, where wavenumber is in CM^-1.
< 
<     DO NW=1,L_NSPECTI
<        a = 1.0D-2/BWNI(NW+1)
<        b = 1.0D-2/BWNI(NW)
<        bpa = (b+a)/2.0
<        bma = (b-a)/2.0
<        do nt=500,9000
---
> !################ Nelli, April 07 ##########################
> !Filling arrays needed by the correlated k lower
> !atmosphere radiation code
> 
>       subroutine radsetup
> 
> !C  GCM2.0  Feb 2003
> !C
> !C     PURPOSE:
> !C        Bundle the new radiation code setup subroutines and call
> !C     this one subroutine from main, where the three included files
> !C     are also listed.  Quantities are passed between this driver
> !C     and the radiation code via common (in radcommon.h).
> !C
> !C----------------------------------------------------------------------C
> 
>       implicit none
> 
> 
>       REAL :: FACTOR
>       INTEGER :: NW
> 
> !C======================================================================C
> 
>       call setspv
>       call setspi
>       call setrad
> 
> !C  Scale IR opacities (Qexti and Qscati) such that
> !C  TAU(0.67 micron)/TAU(9 micron) = VIS2IR, which nominally is 2.
> 
>       QextREF = Qextv(L_NREFV)
>       VIS2IR  = 2.75D0
> 
>       factor  = Qextv(6)/(VIS2IR*Qexti(4))
> 
>       DO NW=1,L_NSPECTI
>         Qexti(NW)  = Qexti(NW)*factor
>         Qscati(NW) = Qscati(NW)*factor
>       END DO
> 
>       PTOP = 10.0**PFGASREF(1)
> 
>       print*,'vis2ir=',vis2ir
> 
>     end subroutine radsetup
> 
> !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>       SUBROUTINE SETSPV
> 
> !C  GCM2.0  Feb 2003
> !C
> !C     PURPOSE:
> !C        Set up the spectral intervals in the Visual (solar).  Based
> !C     on Chris McKay's SETSPV code.
> !C
> !C     AUTHOR
> !C        Jim Schaeffer
> !C
> !C     UPDATES FOR
> !C        Bob Haberle
> !C
> !C     ENVIRONMENT
> !C        Sun ULTRA-2       SOLARIS 2.5.1      FORTRAN
> !C
> !C     REVISION HISTORY
> !C        Original 9/17/1998
> !C     VERSION 2.0  OCT 2001
> !C
> !C     INPUT PARAMETERS
> !C     L_NSPECTV  - Number of spectral intervals in the Visual
> !C
> !C     OUTPUT PARAMETERS
> !C     WNOV       - Array of wavenumbers at the spectral interval
> !C                  center for the VISUAL.  Array is NSPECTV
> !C                  elements long.
> !C     DWNV       - Array of "delta wavenumber", i.e., the width,
> !C                  in wavenumbers (cm^-1) of each VISUAL spectral
> !C                  interval.  NSPECTV elements long.
> !C     WAVEV      - Array (NSPECTV elements long) of the wavelenght
> !C                  (in microns) at the center of each VISUAL spectral
> !C                  interval.
> !C     SOLARF     - Array (NSPECTV elements) of solar flux (W/M^2) in
> !C                  each spectral interval.  Values are for 1 AU, and
> !C                  are scaled to the Mars distance elsewhere.
> !C     TAURAY     - Array (NSPECTV elements) of the wavelength dependent
> !C                  part of Rayleigh Scattering.  The pressure dependent
> !C                  part is computed elsewhere (OPTCV).
> !C     CALLED BY
> !C        RADIATION
> !C
> !C     SUBROUTINES CALLED
> !C        NONE
> !C
> !C**********************************************************************C
> 
>       implicit none
> 
> 
> !C     BWNV - Bin wavenumber of the edges of the VISUAL spectral bins
> !C     units are inverse centimeters.  Dimension needs to be changed
> !C     if the number of VISUAL bins changes.
> 
>       REAL :: BWNV(L_NSPECTV+1)
>       REAL :: SOLAR(L_NSPECTV)
> 
>       REAL ::  P0, GRAV, SCALEP, SUM, WL
>       INTEGER :: N, M
> 
> !C     P0      - Rayleigh scattering reference pressure in pascals.
> !C     GRAV    - Acceleration due to gravity (g) - MKS
> !C     SCALEP  - multiply by 100 to convert pressure from millibars
> !C               to pascals.
> 
>       DATA P0     / 9.423D+6 /
>       DATA GRAV   / 3.72     /
>       DATA SCALEP / 100.0    /
> 
> !C     Bin wavenumber - wavenumber [cm^(-1)] at the edges of the VISUAL
> !C     spectral bins.  Go from smaller to larger wavenumbers, the same as
> !C     in the IR.
> 
>       DATA BWNV / 2222.22D0,  3087.37D0,  4030.63D0,  5370.57D0,&
>                   7651.11D0, 12500.00D0, 25000.00D0, 41666.67D0 /
> 
> !C     Solar flux within each spectral interval, at 1AU (W/M^2)
> !C     Sum equals 1356 W/m^2 (values from Allen, 4th edition)
> 
>       DATA SOLAR /  17.0, 29.0, 52.0, 148.0, 348.0, 643.0, 118.0 /
> 
> !C======================================================================C
> 
> !C     Set up mean wavenumbers and wavenumber deltas.  Units of
> !C     wavenumbers is cm^(-1); units of wavelengths is microns.
> 
>       do M=1,L_NSPECTV
>         WNOV(M)  = 0.5*(BWNV(M+1)+BWNV(M))
>         DWNV(M)  = BWNV(M+1)-BWNV(M)
>         WAVEV(M) = 1.0E+4/WNOV(M)
>       end do
> 
> !C     Sum the solar flux, and write out the result.
> 
>       sum = 0.0
>       do N=1,L_NSPECTV
>         SOLARF(N) = SOLAR(N)
>         sum       = sum+SOLARF(N)
>       end do
>       write(6,'("Solar flux at 1AU = ",f7.2," W/M^2")') sum
> 
> !C     Set up the wavelength dependent part of Rayleigh Scattering.
> !C     The pressure dependent part will be computed elsewhere (OPTCV).
> !C     WAVEV is in microns.  There is no Rayleigh scattering in the IR.
> 
>       do N=1,L_NSPECTV
>         WL        = WAVEV(N)
>         TAURAY(N) = (8.7/grav)*(1.527*(1.0+0.013/wl**2)/wl**4)*&
>                      scalep/P0
>       end do
> 
> 
>       END SUBROUTINE SETSPV
> 
> !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> 
>       subroutine setspi
> 
> !C  GCM2.0  Feb 2003
> !C
> !C     PURPOSE:
> !C        Set up the spectral intervals in the infrared.  Based on
> !C     Chris McKay's SETSPI code.
> !C
> !C     AUTHOR
> !C        Jim Schaeffer
> !C
> !C     UPDATES FOR
> !C        Bob Haberle
> !C
> !C     ENVIRONMENT
> !C        Sun ULTRA-2       SOLARIS 2.5.1      FORTRAN
> !C
> !C     REVISION HISTORY
> !C        Original 9/17/1998
> !C        VERSION 2  OCT 2001
> !C
> !C     INPUT PARAMETERS
> !C     L_NSPECTI  - Number of spectral intervals in the INFRARED
> !C
> !C     OUTPUT PARAMETERS
> !C     WNOI       - Array of wavenumbers at the spectral interval
> !C                  centers for the infrared.  Array is NSPECTI
> !C                  elements long.
> !C     DWNI       - Array of "delta wavenumber", i.e., the width,
> !C                  in wavenumbers (cm^-1) of each IR spectral
> !C                  interval.  NSPECTI elements long.
> !C     WAVEI      - Array (NSPECTI elements long) of the wavelenght
> !C                  (in microns) at the center of each IR spectral
> !C                  interval.
> !C
> !C     CALLED BY
> !C        RADIATION
> !C
> !C     SUBROUTINES CALLED
> !C        NONE
> !C
> !C**********************************************************************C
> 
>       implicit none
> 
> 
> !C     BWNI - Bin wavenumber of the edges of the IR spectral bins
> !C     units are inverse centimeters.  Dimension needs to be changed
> !C     if the number of IR bins changes.
> 
>       REAL :: BWNI(L_NSPECTI+1)
> 
>       real :: a, b, x(12), w(12), ans, y, bpa, bma, T
>       real :: c1, c2, wn1, wn2
>       integer :: n, nw, nt, m
> 
> !C  C1 and C2 values from Goody and Yung (2nd edition)  MKS units
> !C  These values lead to a "sigma" (sigma*T^4) of 5.67032E-8 W m^-2 K^-4
> 
>       data c1 / 3.741832D-16 /     ! W m^-2
>       data c2 / 1.438786D-2  /     ! m K
> 
>       data x / -0.981560634246719D0,  -0.904117256370475D0,&
>                -0.769902674194305D0,  -0.587317954286617D0,&
>                -0.367831498998180D0,  -0.125233408511469D0,&
>                 0.125233408511469D0,   0.367831498998180D0,&
>                 0.587317954286617D0,   0.769902674194305D0,&
>                 0.904117256370475D0,   0.981560634246719D0    /
> 
>       data w /  0.047175336386512D0,   0.106939325995318D0,&
>                 0.160078328543346D0,   0.203167426723066D0,&
>                 0.233492536538355D0,   0.249147045813403D0,&
>                 0.249147045813403D0,   0.233492536538355D0,&
>                 0.203167426723066D0,   0.160078328543346D0,&
>                 0.106939325995318D0,   0.047175336386512D0   /
> 
> !C======================================================================C
> 
> !C     Bin wavenumber - wavenumber [cm^(-1)] at the edges of the IR
> !C     spectral bins.
> 
>       BWNI( 1) =   10.000D0
>       BWNI( 2) =  166.667D0
>       BWNI( 3) =  416.667D0
>       BWNI( 4) =  833.333D0
>       BWNI( 5) = 1250.000D0
>       BWNI( 6) = 2500.000D0
> 
> !C     Set up mean wavenumbers and wavenumber deltas.  Units of
> !C     wavenumbers is cm^(-1); units of wavelengths is microns.
> 
>       do M=1,L_NSPECTI
>         WNOI(M)  = 0.5*(BWNI(M+1)+BWNI(M))
>         DWNI(M)  = BWNI(M+1)-BWNI(M)
>         WAVEI(M) = 1.0E+4/WNOI(M)
>       end do
> 
> !C  For each IR wavelength interval, compute the integral of B(T), the
> !C  Planck function, divided by the wavelength interval, in cm-1.  The
> !C  integration is in MKS units, the final answer is the same as the
> !C  original planck.f; W m^-2 wavenumber^-1, where wavenumber is in CM^-1.
> 
>       DO NW=1,L_NSPECTI
>         a = 1.0D-2/BWNI(NW+1)
>         b = 1.0D-2/BWNI(NW)
>         bpa = (b+a)/2.0
>         bma = (b-a)/2.0
>         do nt=500,9000
1307,1308c1309,1310
<              y    = bma*x(m)+bpa
<              ans  = ans + w(m)*c1/(y**5*(exp(c2/(y*T))-1.0D0))
---
>             y    = bma*x(m)+bpa
>             ans  = ans + w(m)*c1/(y**5*(exp(c2/(y*T))-1.0D0))
1311,1312c1313,1314
<        end do
<     END DO
---
>         end do
>       END DO
1315c1317
<   end subroutine setspi
---
>       end subroutine setspi
1317c1319
<   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
---
> !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1319,1464c1321,1469
<   subroutine setrad
< 
<     !C  GCM2.0  Feb 2003
<     !C
<     !C     PURPOSE:
<     !C        Set up values used by the radiation code, such as the CO2 gas
<     !C     absorption coefficients.  True constants are defined, and the
<     !C     time-independent quantities used by the radiation code are
<     !C     calculated.
<     !C
<     !C     AUTHOR
<     !C
<     !C
<     !C     UPDATES FOR
<     !C        Jim Pollack
<     !C
<     !C     ENVIRONMENT
<     !C        NAS Cray-YMP      UNICOS       FORTRAN
<     !C
<     !C     REVISION HISTORY
<     !C        See notes from J. Pollack dated 2/20/90 and 3/12/90.
<     !C        Modified Cmk  See note from Jim Pollack dated 1/11/91.
<     !C        Several - See note dated 1/23/91 from Jim Pollack.
<     !C     VERSION 2.0  OCT 2001
<     !C
<     !C     INPUT PARAMETERS
<     !C     DTAU(L,M)      - Dust optical depth of layer L, and for aerosol
<     !C                      species M.
<     !C     ptrop          - Pressure of the tropopause (mb)
<     !C     SCALEP         - Factor to convert pressures from millibars to
<     !C                      Pascals
<     !C
<     !C     OUTPUT PARAMETERS
<     !C     PTOP       - Pressure at the TOP of the stratosphere ( or
<     !C                  equivalently, at the bottom of the dummy layer.
<     !C                  PTOP = PTROP*SDUMMY
<     !C
<     !C     AEROSOL RADIATIVE OPTICAL CONSTANTS
<     !C     Values are at the wavelenght interval center
<     !C
<     !C     MIE SCATTERING - Size distribution weighted
<     !C     Qextv    - Extinction efficiency - in the visual.
<     !C     Qscatv   - Scattering efficiency - in the visual.
<     !C     WV       - Single scattering albedo - in the visual.
<     !C     GV       - Asymmetry parameter - in the visual.
<     !C
<     !C     Qexti    - Extinction efficiency - in the infrared.
<     !C     Qscati   - Scattering efficiency - in the infrared.
<     !C     WI       - Single scattering albedo - in the infrared.
<     !C     GI       - Asymmetry parameter - in the infrared.
<     !C
<     !C     CALLED BY
<     !C        RAD
<     !C  
<     !C     SUBROUTINES CALLED
<     !C        DMIESS, PLNK
<     !C
<     !C----------------------------------------------------------------------C
< 
<     implicit none
< 
<     integer :: N, NS
< 
<     real :: qev1(L_NSPECTV)
<     real :: qsv1(L_NSPECTV)
<     real :: gv1(L_NSPECTV)
< 
<     real :: qei1(L_NSPECTI)
<     real :: qsi1(L_NSPECTI)
<     real :: gi1(L_NSPECTI)
< 
<     integer :: nt, nw, ng
< 
<     !C----------------------------------------------------------------------C
< 
<     !C  Visible dust properties:  Ockert-Bell Planck-weighted values (T=6000K)
< 
<     !C     Qext - Ockert-Bell values (order is increasing waveNUMBER)
<     !C     VISULAL WAVELENGTHS.
< 
<     data qev1 / 2.529D0, 2.949D0, 3.209D0, 3.337D0, 3.207D0,&
<          2.938D0, 2.622D0                                   /
< 
<     !C     Qscat - Ockert-Bell values
<     !C     VISUAL wavelengths
< 
<     data qsv1 / 2.374D0, 2.637D0, 3.049D0, 3.201D0, 3.045D0,&
<          2.513D0, 1.623D0                                   /
< 
<     !C     G - Ockert-Bell values
<     !C     VISUAL wavelengths
< 
<     data gv1  / 0.635D0, 0.646D0, 0.630D0, 0.630D0, 0.634D0,&
<          0.700D0, 0.856D0                                   /
< 
< 
<     !C     And now the INFRARED
< 
<     !C     Qext for a modified-gamma distribution, ALPHA=2, GAMMA=0.5,
<     !C     Rm = 0.4 microns, using the Forget optical constants (Nr and Ni).
<     !C     Planck-weighted values (T=215K)
<     !C     INFRARED wavelengths.  (The order is increasing waveNUMBER.)
< 
<     data qei1 / 0.193D0, 0.867D0, 1.209D0, 2.173D0, 0.638D0  /
< 
<     !C     Qsca for a modified gamma-distribution, using the Forget
<     !C     optical constants (Nr and Ni).       INFRARED wavelengths
< 
<     data qsi1 / 0.027D0, 0.319D0, 0.558D0, 1.136D0, 0.237D0  /
< 
<     !C     g for a modified gamma-distribution, using the Forget
<     !C     optical constants (Nr and Ni).       INFRARED wavelengths
< 
<     data gi1  / 0.024D0, 0.127D0, 0.288D0, 0.423D0, 0.548D0  /
< 
<     !C=======================================================================
< 
<     !C     Set the reference pressure and temperature arrays.  These are
<     !C     the pressures and temperatures at which we have k-coefficients.
< 
<     pgasref( 1) = 1.0E-6
<     pgasref( 2) = 1.0E-5
<     pgasref( 3) = 1.0E-4
<     pgasref( 4) = 1.0E-3
<     pgasref( 5) = 1.0E-2
<     pgasref( 6) = 1.0E-1
<     pgasref( 7) = 1.0
<     pgasref( 8) = 1.0E+1
<     pgasref( 9) = 1.0E+2
<     pgasref(10) = 1.0E+3
<     pgasref(11) = 1.0E+4
< 
<     tgasref(1)  =  50.0
<     tgasref(2)  = 100.0
<     tgasref(3)  = 150.0
<     tgasref(4)  = 200.0
<     tgasref(5)  = 250.0
<     tgasref(6)  = 300.0
<     tgasref(7)  = 350.0
< 
<     !C     Fill the (VISUAL) arrays Qextv, Qscatv, WV, GV
< 
<     DO N=1,L_NSPECTV
<        Qextv(n)  = qev1(n)
<        Qscatv(n) = qsv1(n)
<        IF(Qscatv(n).GE.Qextv(n)) then
---
>       subroutine setrad
> 
> !C  GCM2.0  Feb 2003
> !C
> !C     PURPOSE:
> !C        Set up values used by the radiation code, such as the CO2 gas
> !C     absorption coefficients.  True constants are defined, and the
> !C     time-independent quantities used by the radiation code are
> !C     calculated.
> !C
> !C     AUTHOR
> !C
> !C
> !C     UPDATES FOR
> !C        Jim Pollack
> !C
> !C     ENVIRONMENT
> !C        NAS Cray-YMP      UNICOS       FORTRAN
> !C
> !C     REVISION HISTORY
> !C        See notes from J. Pollack dated 2/20/90 and 3/12/90.
> !C        Modified Cmk  See note from Jim Pollack dated 1/11/91.
> !C        Several - See note dated 1/23/91 from Jim Pollack.
> !C     VERSION 2.0  OCT 2001
> !C
> !C     INPUT PARAMETERS
> !C     DTAU(L,M)      - Dust optical depth of layer L, and for aerosol
> !C                      species M.
> !C     ptrop          - Pressure of the tropopause (mb)
> !C     SCALEP         - Factor to convert pressures from millibars to
> !C                      Pascals
> !C
> !C     OUTPUT PARAMETERS
> !C     PTOP       - Pressure at the TOP of the stratosphere ( or
> !C                  equivalently, at the bottom of the dummy layer.
> !C                  PTOP = PTROP*SDUMMY
> !C
> !C     AEROSOL RADIATIVE OPTICAL CONSTANTS
> !C     Values are at the wavelenght interval center
> !C
> !C     MIE SCATTERING - Size distribution weighted
> !C     Qextv    - Extinction efficiency - in the visual.
> !C     Qscatv   - Scattering efficiency - in the visual.
> !C     WV       - Single scattering albedo - in the visual.
> !C     GV       - Asymmetry parameter - in the visual.
> !C
> !C     Qexti    - Extinction efficiency - in the infrared.
> !C     Qscati   - Scattering efficiency - in the infrared.
> !C     WI       - Single scattering albedo - in the infrared.
> !C     GI       - Asymmetry parameter - in the infrared.
> !C
> !C     CALLED BY
> !C        RAD
> !C  
> !C     SUBROUTINES CALLED
> !C        DMIESS, PLNK
> !C
> !C----------------------------------------------------------------------C
> 
>       implicit none
> 
> 
>       integer :: N, NS
> 
> 
>       real :: qev1(L_NSPECTV)
>       real :: qsv1(L_NSPECTV)
>       real :: gv1(L_NSPECTV)
> 
>       real :: qei1(L_NSPECTI)
>       real :: qsi1(L_NSPECTI)
>       real :: gi1(L_NSPECTI)
> 
>       integer :: nt, nw, ng
> 
> 
> !C----------------------------------------------------------------------C
> 
> !C  Visible dust properties:  Ockert-Bell Planck-weighted values (T=6000K)
> 
> !C     Qext - Ockert-Bell values (order is increasing waveNUMBER)
> !C     VISULAL WAVELENGTHS.
> 
>       data qev1 / 2.529D0, 2.949D0, 3.209D0, 3.337D0, 3.207D0,&
>                   2.938D0, 2.622D0                                   /
> 
> !C     Qscat - Ockert-Bell values
> !C     VISUAL wavelengths
> 
>       data qsv1 / 2.374D0, 2.637D0, 3.049D0, 3.201D0, 3.045D0,&
>                   2.513D0, 1.623D0                                   /
> 
> !C     G - Ockert-Bell values
> !C     VISUAL wavelengths
> 
>       data gv1  / 0.635D0, 0.646D0, 0.630D0, 0.630D0, 0.634D0,&
>                   0.700D0, 0.856D0                                   /
> 
> 
> !C     And now the INFRARED
> 
> !C     Qext for a modified-gamma distribution, ALPHA=2, GAMMA=0.5,
> !C     Rm = 0.4 microns, using the Forget optical constants (Nr and Ni).
> !C     Planck-weighted values (T=215K)
> !C     INFRARED wavelengths.  (The order is increasing waveNUMBER.)
> 
>       data qei1 / 0.193D0, 0.867D0, 1.209D0, 2.173D0, 0.638D0  /
> 
> !C     Qsca for a modified gamma-distribution, using the Forget
> !C     optical constants (Nr and Ni).       INFRARED wavelengths
> 
>       data qsi1 / 0.027D0, 0.319D0, 0.558D0, 1.136D0, 0.237D0  /
> 
> !C     g for a modified gamma-distribution, using the Forget
> !C     optical constants (Nr and Ni).       INFRARED wavelengths
> 
>       data gi1  / 0.024D0, 0.127D0, 0.288D0, 0.423D0, 0.548D0  /
> 
> !C=======================================================================
> 
> !C     Set the reference pressure and temperature arrays.  These are
> !C     the pressures and temperatures at which we have k-coefficients.
> 
>       pgasref( 1) = 1.0E-6
>       pgasref( 2) = 1.0E-5
>       pgasref( 3) = 1.0E-4
>       pgasref( 4) = 1.0E-3
>       pgasref( 5) = 1.0E-2
>       pgasref( 6) = 1.0E-1
>       pgasref( 7) = 1.0
>       pgasref( 8) = 1.0E+1
>       pgasref( 9) = 1.0E+2
>       pgasref(10) = 1.0E+3
>       pgasref(11) = 1.0E+4
> 
>       tgasref(1)  =  50.0
>       tgasref(2)  = 100.0
>       tgasref(3)  = 150.0
>       tgasref(4)  = 200.0
>       tgasref(5)  = 250.0
>       tgasref(6)  = 300.0
>       tgasref(7)  = 350.0
> 
> !C     Fill the (VISUAL) arrays Qextv, Qscatv, WV, GV
> 
>       DO N=1,L_NSPECTV
>         Qextv(n)  = qev1(n)
>         Qscatv(n) = qsv1(n)
>         IF(Qscatv(n).GE.Qextv(n)) then
1466,1476c1471,1481
<        END IF
<        WV(n)     = Qscatv(n)/Qextv(n)
<        GV(n)     = gv1(n)
<     END DO
< 
<     !C     Fill the (INFRARED) arrays Qexti, Qscati, WI, GI
< 
<     DO N=1,L_NSPECTI
<        Qexti(n)  = qei1(n)
<        Qscati(n) = qsi1(n)
<        IF(Qscati(n).GE.Qexti(n)) then
---
>         END IF
>         WV(n)     = Qscatv(n)/Qextv(n)
>         GV(n)     = gv1(n)
>       END DO
> 
> !C     Fill the (INFRARED) arrays Qexti, Qscati, WI, GI
> 
>       DO N=1,L_NSPECTI
>         Qexti(n)  = qei1(n)
>         Qscati(n) = qsi1(n)
>         IF(Qscati(n).GE.Qexti(n)) then
1478,1481c1483,1486
<        END IF
<        WI(n)     = Qscati(n)/Qexti(n)
<        GI(n)     = gi1(n)
<     END DO
---
>         END IF
>         WI(n)     = Qscati(n)/Qexti(n)
>         GI(n)     = gi1(n)
>       END DO
1483,1484c1488,1489
<     !C     Get CO2 k coefficients, and interpolate them to the finer
<     !C     pressure grid.
---
> !C     Get CO2 k coefficients, and interpolate them to the finer
> !C     pressure grid.
1487c1492
<     call laginterp(pgasref,pfgasref)
---
>       call laginterp(pgasref,pfgasref)
1490c1495
<   end subroutine setrad
---
>       end subroutine setrad
1492c1497
<   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
---
> !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1494c1499
<   subroutine laginterp(pgref,pint)
---
>       subroutine laginterp(pgref,pint)
1496c1501
<     !C  GCM2.0  Feb 2003
---
> !C  GCM2.0  Feb 2003
1498,1500c1503,1505
<     !C  Lagrange interpolation (linear in log pressure) of the CO2
<     !C  k-coefficients in the pressure domain.  Subsequent use of these
<     !C  values will use a simple linear interpolation in pressure.
---
> !C  Lagrange interpolation (linear in log pressure) of the CO2
> !C  k-coefficients in the pressure domain.  Subsequent use of these
> !C  values will use a simple linear interpolation in pressure.
1502c1507
<     implicit none
---
>       implicit none
1505,1508c1510,1513
<     integer :: n, nt, np, nh, ng, nw, m, i, j, k, l
<     real :: co2i8(L_NTREF,L_NPREF,L_REFH2O,L_NSPECTI,L_NGAUSS)
<     real :: co2v8(L_NTREF,L_NPREF,L_REFH2O,L_NSPECTV,L_NGAUSS)
<     real :: pgref(L_NPREF)
---
>       integer :: n, nt, np, nh, ng, nw, m, i, j, k, l
>       real :: co2i8(L_NTREF,L_NPREF,L_REFH2O,L_NSPECTI,L_NGAUSS)
>       real :: co2v8(L_NTREF,L_NPREF,L_REFH2O,L_NSPECTV,L_NGAUSS)
>       real :: pgref(L_NPREF)
1511,1512c1516,1517
<     real :: x, xi(4), yi(4), ans
<     real :: pint(L_PINT), pin(L_PINT), pref(L_NPREF), p
---
>       real :: x, xi(4), yi(4), ans
>       real :: pint(L_PINT), pin(L_PINT), pref(L_NPREF), p
1514,1524c1519,1529
<     data pin / -6.0D0, -5.8D0, -5.6D0, -5.4D0, -5.2D0,&
<          -5.0D0, -4.8D0, -4.6D0, -4.4D0, -4.2D0,&
<          -4.0D0, -3.8D0, -3.6D0, -3.4D0, -3.2D0,&
<          -3.0D0, -2.8D0, -2.6D0, -2.4D0, -2.2D0,&
<          -2.0D0, -1.8D0, -1.6D0, -1.4D0, -1.2D0,&
<          -1.0D0, -0.8D0, -0.6D0, -0.4D0, -0.2D0,&
<          0.0D0,  0.2D0,  0.4D0,  0.6D0,  0.8D0,&
<          1.0D0,  1.2D0,  1.4D0,  1.6D0,  1.8D0,&
<          2.0D0,  2.2D0,  2.4D0,  2.6D0,  2.8D0,&
<          3.0D0,  3.2D0,  3.4D0,  3.6D0,  3.8D0,&
<          4.0D0                                    /
---
>       data pin / -6.0D0, -5.8D0, -5.6D0, -5.4D0, -5.2D0,&
>                  -5.0D0, -4.8D0, -4.6D0, -4.4D0, -4.2D0,&
>                  -4.0D0, -3.8D0, -3.6D0, -3.4D0, -3.2D0,&
>                  -3.0D0, -2.8D0, -2.6D0, -2.4D0, -2.2D0,&
>                  -2.0D0, -1.8D0, -1.6D0, -1.4D0, -1.2D0,&
>                  -1.0D0, -0.8D0, -0.6D0, -0.4D0, -0.2D0,&
>                   0.0D0,  0.2D0,  0.4D0,  0.6D0,  0.8D0,&
>                   1.0D0,  1.2D0,  1.4D0,  1.6D0,  1.8D0,&
>                   2.0D0,  2.2D0,  2.4D0,  2.6D0,  2.8D0,&
>                   3.0D0,  3.2D0,  3.4D0,  3.6D0,  3.8D0,&
>                   4.0D0                                    /
1526c1531
<     !C======================================================================!
---
> !C======================================================================!
1528c1533
<     !C  Fill pint for output from this subroutine
---
> !C  Fill pint for output from this subroutine
1530,1532c1535,1537
<     do n=1,L_PINT
<        PINT(n) = PIN(n)
<     end do
---
>       do n=1,L_PINT
>         PINT(n) = PIN(n)
>       end do
1534c1539
<     !C  Take log of the reference pressures
---
> !C  Take log of the reference pressures
1536,1538c1541,1547
<     do n=1,L_NPREF
<        pref(n) = LOG10(PGREF(n))
<     end do
---
>       do n=1,L_NPREF
>         pref(n) = LOG10(PGREF(n))
>       end do
> 
> !C     Get CO2 k coefficients
> 
> 	OPEN(Unit=25,file='UA/DataIn/CO2H2O_V_12_95_ASCII',ACTION = 'READ')
1540c1549,1559
<     !C     Get CO2 k coefficients
---
> 	DO i=1,L_NTREF
> 	   DO j=1,L_NPREF
> 	      DO k=1,L_REFH2O
> 		 DO l=1,L_NSPECTV
> 		    DO m=1,L_NGAUSS
> 		       read(25,*)co2v8(i,j,k,l,m)
> 		    ENDDO
> 		 ENDDO
>               ENDDO
> 	   ENDDO
> 	ENDDO
1542c1561,1564
<     OPEN(Unit=25,file='UA/DataIn/CO2H2O_V_12_95_ASCII',ACTION = 'READ')
---
> 	DO i=1,L_NSPECTV
> 	    read(25,*)fzerov(i)
> 	ENDDO
> 	CLOSE(25)
1544,1578c1566,1583
<     DO i=1,L_NTREF
<        DO j=1,L_NPREF
<           DO k=1,L_REFH2O
<              DO l=1,L_NSPECTV
<                 DO m=1,L_NGAUSS
<                    read(25,*)co2v8(i,j,k,l,m)
<                 ENDDO
<              ENDDO
<           ENDDO
<        ENDDO
<     ENDDO
< 
<     DO i=1,L_NSPECTV
<        read(25,*)fzerov(i)
<     ENDDO
<     CLOSE(25)
< 
<     OPEN(unit=35,file='UA/DataIn/CO2H2O_IR_12_95_ASCII',ACTION = 'READ')
< 
<     DO i=1,L_NTREF
<        DO j=1,L_NPREF
<           DO k=1,L_REFH2O
<              DO l=1,L_NSPECTI
<                 DO m=1,L_NGAUSS
<                    read(35,*)co2i8(i,j,k,l,m)
<                 ENDDO
<              ENDDO
<           ENDDO
<        ENDDO
<     ENDDO
< 
<     DO i=1,L_NSPECTI
<        read(35,*)fzeroi(i)
<     ENDDO
<     CLOSE(35)
---
> 	OPEN(unit=35,file='UA/DataIn/CO2H2O_IR_12_95_ASCII',ACTION = 'READ')
> 
> 	DO i=1,L_NTREF
> 	   DO j=1,L_NPREF
> 	      DO k=1,L_REFH2O
> 		 DO l=1,L_NSPECTI
> 		    DO m=1,L_NGAUSS
> 		       read(35,*)co2i8(i,j,k,l,m)
> 		    ENDDO
> 		 ENDDO
>               ENDDO
> 	   ENDDO
> 	ENDDO
> 
> 	DO i=1,L_NSPECTI
> 	   read(35,*)fzeroi(i)
> 	ENDDO
> 	CLOSE(35)
1587,1588c1592,1593
<     !C  Take Log10 of the values - we interpolate the log10 of the values,
<     !C  not the values themselves.   Smallest value is 1.0E-200.
---
> !C  Take Log10 of the values - we interpolate the log10 of the values,
> !C  not the values themselves.   Smallest value is 1.0E-200.
1590,1591c1595,1596
<     do nt=1,L_NTREF
<        do np=1,L_NPREF
---
>       do nt=1,L_NTREF
>         do np=1,L_NPREF
1593c1598
<              do ng = 1,L_NGAUSS
---
>             do ng = 1,L_NGAUSS
1595,1609c1600,1614
<                 do nw=1,L_NSPECTV
<                    if(co2v8(nt,np,nh,nw,ng).gt.1.0e-200) then
<                       co2v8(nt,np,nh,nw,ng) = log10(co2v8(nt,np,nh,nw,ng))
<                    else
<                       co2v8(nt,np,nh,nw,ng) = -200.0
<                    end if
<                 end do
< 
<                 do nw=1,L_NSPECTI
<                    if(co2i8(nt,np,nh,nw,ng).gt.1.0e-200) then
<                       co2i8(nt,np,nh,nw,ng) = log10(co2i8(nt,np,nh,nw,ng))
<                    else
<                       co2i8(nt,np,nh,nw,ng) = -200.0
<                    end if
<                 end do
---
>               do nw=1,L_NSPECTV
>                 if(co2v8(nt,np,nh,nw,ng).gt.1.0e-200) then
>                   co2v8(nt,np,nh,nw,ng) = log10(co2v8(nt,np,nh,nw,ng))
>                 else
>                   co2v8(nt,np,nh,nw,ng) = -200.0
>                 end if
>               end do
> 
>               do nw=1,L_NSPECTI
>                 if(co2i8(nt,np,nh,nw,ng).gt.1.0e-200) then
>                   co2i8(nt,np,nh,nw,ng) = log10(co2i8(nt,np,nh,nw,ng))
>                 else
>                   co2i8(nt,np,nh,nw,ng) = -200.0
>                 end if
>               end do
1611c1616
<              end do
---
>             end do
1613,1614c1618,1619
<        end do
<     end do
---
>         end do
>       end do
1616c1621
<     !C  Interpolate the values:  first the IR
---
> !C  Interpolate the values:  first the IR
1618,1674c1623,1679
<     do nt=1,L_NTREF
<        do nh=1,L_REFH2O
<           do nw=1,L_NSPECTI
<              do ng=1,L_NGAUSS
< 
<                 !C  First, the initial interval (P=1e-6 to 1e-5)
< 
<                 n = 1
<                 do m=1,5
<                    x     = pint(m)
<                    xi(1) = pref(n)
<                    xi(2) = pref(n+1)
<                    xi(3) = pref(n+2)
<                    xi(4) = pref(n+3)
<                    yi(1) = co2i8(nt,n,nh,nw,ng)
<                    yi(2) = co2i8(nt,n+1,nh,nw,ng)
<                    yi(3) = co2i8(nt,n+2,nh,nw,ng)
<                    yi(4) = co2i8(nt,n+3,nh,nw,ng)
<                    call lagrange(x,xi,yi,ans)
<                    co2i(nt,m,nh,nw,ng) = 10.0**ans
<                 end do
< 
<                 do n=2,L_NPREF-2
<                    do m=1,5
<                       i     = (n-1)*5+m
<                       x     = pint(i)
<                       xi(1) = pref(n-1)
<                       xi(2) = pref(n)
<                       xi(3) = pref(n+1)
<                       xi(4) = pref(n+2)
<                       yi(1) = co2i8(nt,n-1,nh,nw,ng)
<                       yi(2) = co2i8(nt,n,nh,nw,ng)
<                       yi(3) = co2i8(nt,n+1,nh,nw,ng)
<                       yi(4) = co2i8(nt,n+2,nh,nw,ng)
<                       call lagrange(x,xi,yi,ans)
<                       co2i(nt,i,nh,nw,ng) = 10.0**ans
<                    end do
<                 end do
< 
<                 !C  Now, get the last interval (P=1e+3 to 1e+4)
< 
<                 n = L_NPREF-1
< 
<                 do m=1,5
<                    i     = (n-1)*5+m
<                    x     = pint(i)
<                    xi(1) = pref(n-2)
<                    xi(2) = pref(n-1)
<                    xi(3) = pref(n)
<                    xi(4) = pref(n+1)
<                    yi(1) = co2i8(nt,n-2,nh,nw,ng)
<                    yi(2) = co2i8(nt,n-1,nh,nw,ng)
<                    yi(3) = co2i8(nt,n,nh,nw,ng)
<                    yi(4) = co2i8(nt,n+1,nh,nw,ng)
<                    call lagrange(x,xi,yi,ans)
<                    co2i(nt,i,nh,nw,ng) = 10.0**ans
<                 end do
---
>       do nt=1,L_NTREF
>         do nh=1,L_REFH2O
>         do nw=1,L_NSPECTI
>           do ng=1,L_NGAUSS
> 
> !C  First, the initial interval (P=1e-6 to 1e-5)
> 
>             n = 1
>             do m=1,5
>               x     = pint(m)
>               xi(1) = pref(n)
>               xi(2) = pref(n+1)
>               xi(3) = pref(n+2)
>               xi(4) = pref(n+3)
>               yi(1) = co2i8(nt,n,nh,nw,ng)
>               yi(2) = co2i8(nt,n+1,nh,nw,ng)
>               yi(3) = co2i8(nt,n+2,nh,nw,ng)
>               yi(4) = co2i8(nt,n+3,nh,nw,ng)
>               call lagrange(x,xi,yi,ans)
>               co2i(nt,m,nh,nw,ng) = 10.0**ans
>             end do
> 
>             do n=2,L_NPREF-2
>               do m=1,5
>                 i     = (n-1)*5+m
>                 x     = pint(i)
>                 xi(1) = pref(n-1)
>                 xi(2) = pref(n)
>                 xi(3) = pref(n+1)
>                 xi(4) = pref(n+2)
>                 yi(1) = co2i8(nt,n-1,nh,nw,ng)
>                 yi(2) = co2i8(nt,n,nh,nw,ng)
>                 yi(3) = co2i8(nt,n+1,nh,nw,ng)
>                 yi(4) = co2i8(nt,n+2,nh,nw,ng)
>                 call lagrange(x,xi,yi,ans)
>                 co2i(nt,i,nh,nw,ng) = 10.0**ans
>               end do
>             end do
> 
> !C  Now, get the last interval (P=1e+3 to 1e+4)
> 
>             n = L_NPREF-1
> 
>             do m=1,5
>               i     = (n-1)*5+m
>               x     = pint(i)
>               xi(1) = pref(n-2)
>               xi(2) = pref(n-1)
>               xi(3) = pref(n)
>               xi(4) = pref(n+1)
>               yi(1) = co2i8(nt,n-2,nh,nw,ng)
>               yi(2) = co2i8(nt,n-1,nh,nw,ng)
>               yi(3) = co2i8(nt,n,nh,nw,ng)
>               yi(4) = co2i8(nt,n+1,nh,nw,ng)
>               call lagrange(x,xi,yi,ans)
>               co2i(nt,i,nh,nw,ng) = 10.0**ans
>             end do
1676c1681
<                 !C  Fill the last pressure point
---
> !C  Fill the last pressure point
1678c1683
<                 co2i(nt,L_PINT,nh,nw,ng) = 10.0**co2i8(nt,L_NPREF,nh,nw,ng)
---
>             co2i(nt,L_PINT,nh,nw,ng) = 10.0**co2i8(nt,L_NPREF,nh,nw,ng)
1680d1684
<              end do
1682,1743c1686,1748
<        end do
<     end do
< 
<     !C  Interpolate the values:  now the Visual
< 
<     do nt=1,L_NTREF
<        do nh=1,L_REFH2O
<           do nw=1,L_NSPECTV
<              do ng=1,L_NGAUSS
< 
<                 !C  First, the initial interval (P=1e-6 to 1e-5)
< 
<                 n = 1
<                 do m=1,5
<                    x     = pint(m)
<                    xi(1) = pref(n)
<                    xi(2) = pref(n+1)
<                    xi(3) = pref(n+2)
<                    xi(4) = pref(n+3)
<                    yi(1) = co2v8(nt,n,nh,nw,ng)
<                    yi(2) = co2v8(nt,n+1,nh,nw,ng)
<                    yi(3) = co2v8(nt,n+2,nh,nw,ng)
<                    yi(4) = co2v8(nt,n+3,nh,nw,ng)
<                    call lagrange(x,xi,yi,ans)
<                    co2v(nt,m,nh,nw,ng) = 10.0**ans
<                 end do
< 
<                 do n=2,L_NPREF-2
<                    do m=1,5
<                       i     = (n-1)*5+m
<                       x     = pint(i)
<                       xi(1) = pref(n-1)
<                       xi(2) = pref(n)
<                       xi(3) = pref(n+1)
<                       xi(4) = pref(n+2)
<                       yi(1) = co2v8(nt,n-1,nh,nw,ng)
<                       yi(2) = co2v8(nt,n,nh,nw,ng)
<                       yi(3) = co2v8(nt,n+1,nh,nw,ng)
<                       yi(4) = co2v8(nt,n+2,nh,nw,ng)
<                       call lagrange(x,xi,yi,ans)
<                       co2v(nt,i,nh,nw,ng) = 10.0**ans
<                    end do
<                 end do
< 
<                 !C  Now, get the last interval (P=1e+3 to 1e+4)
< 
<                 n = L_NPREF-1
< 
<                 do m=1,5
<                    i     = (n-1)*5+m
<                    x     = pint(i)
<                    xi(1) = pref(n-2)
<                    xi(2) = pref(n-1)
<                    xi(3) = pref(n)
<                    xi(4) = pref(n+1)
<                    yi(1) = co2v8(nt,n-2,nh,nw,ng)
<                    yi(2) = co2v8(nt,n-1,nh,nw,ng)
<                    yi(3) = co2v8(nt,n,nh,nw,ng)
<                    yi(4) = co2v8(nt,n+1,nh,nw,ng)
<                    call lagrange(x,xi,yi,ans)
<                    co2v(nt,i,nh,nw,ng) = 10.0**ans
<                 end do
---
>         end do
>         end do
>       end do
> 
> !C  Interpolate the values:  now the Visual
> 
>       do nt=1,L_NTREF
>         do nh=1,L_REFH2O
>         do nw=1,L_NSPECTV
>           do ng=1,L_NGAUSS
> 
> !C  First, the initial interval (P=1e-6 to 1e-5)
> 
>             n = 1
>             do m=1,5
>               x     = pint(m)
>               xi(1) = pref(n)
>               xi(2) = pref(n+1)
>               xi(3) = pref(n+2)
>               xi(4) = pref(n+3)
>               yi(1) = co2v8(nt,n,nh,nw,ng)
>               yi(2) = co2v8(nt,n+1,nh,nw,ng)
>               yi(3) = co2v8(nt,n+2,nh,nw,ng)
>               yi(4) = co2v8(nt,n+3,nh,nw,ng)
>               call lagrange(x,xi,yi,ans)
>               co2v(nt,m,nh,nw,ng) = 10.0**ans
>             end do
> 
>             do n=2,L_NPREF-2
>               do m=1,5
>                 i     = (n-1)*5+m
>                 x     = pint(i)
>                 xi(1) = pref(n-1)
>                 xi(2) = pref(n)
>                 xi(3) = pref(n+1)
>                 xi(4) = pref(n+2)
>                 yi(1) = co2v8(nt,n-1,nh,nw,ng)
>                 yi(2) = co2v8(nt,n,nh,nw,ng)
>                 yi(3) = co2v8(nt,n+1,nh,nw,ng)
>                 yi(4) = co2v8(nt,n+2,nh,nw,ng)
>                 call lagrange(x,xi,yi,ans)
>                 co2v(nt,i,nh,nw,ng) = 10.0**ans
>               end do
>             end do
> 
> !C  Now, get the last interval (P=1e+3 to 1e+4)
> 
>             n = L_NPREF-1
> 
>             do m=1,5
>               i     = (n-1)*5+m
>               x     = pint(i)
>               xi(1) = pref(n-2)
>               xi(2) = pref(n-1)
>               xi(3) = pref(n)
>               xi(4) = pref(n+1)
>               yi(1) = co2v8(nt,n-2,nh,nw,ng)
>               yi(2) = co2v8(nt,n-1,nh,nw,ng)
>               yi(3) = co2v8(nt,n,nh,nw,ng)
>               yi(4) = co2v8(nt,n+1,nh,nw,ng)
>               call lagrange(x,xi,yi,ans)
>               co2v(nt,i,nh,nw,ng) = 10.0**ans
>             end do
1745c1750
<                 !C  Fill the last pressure point
---
> !C  Fill the last pressure point
1747c1752
<                 co2v(nt,L_PINT,nh,nw,ng) = 10.0**co2v8(nt,L_NPREF,nh,nw,ng)
---
>             co2v(nt,L_PINT,nh,nw,ng) = 10.0**co2v8(nt,L_NPREF,nh,nw,ng)
1749d1753
<              end do
1751,1752c1755,1757
<        end do
<     end do
---
>         end do
>         end do
>       end do
1755c1760
<   end subroutine laginterp
---
>       end subroutine laginterp
1757c1762
<   !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
---
> !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1759c1764
<   subroutine lagrange(x, xi, yi, ans)
---
>       subroutine lagrange(x, xi, yi, ans)
1761,1764c1766,1769
<     !C  GCM2.0  Feb 2003
<     !C
<     !C  Lagrange interpolation - Polynomial interpolation at point x
<     !C  xi(1) <= x <= xi(4).  Yi(n) is the functional value at XI(n).
---
> !C  GCM2.0  Feb 2003
> !C
> !C  Lagrange interpolation - Polynomial interpolation at point x
> !C  xi(1) <= x <= xi(4).  Yi(n) is the functional value at XI(n).
1766c1771
<     implicit none
---
>       implicit none
1768,1769c1773,1774
<     real :: x, xi(4), yi(4), ans
<     real :: fm1, fm2, fm3, fm4
---
>       real :: x, xi(4), yi(4), ans
>       real :: fm1, fm2, fm3, fm4
1771c1776
<     !C======================================================================!
---
> !C======================================================================!
1773,1776c1778,1781
<     fm1   = x - XI(1)
<     fm2   = x - XI(2)
<     fm3   = x - XI(3)
<     fm4   = x - XI(4)
---
>       fm1   = x - XI(1)
>       fm2   = x - XI(2)
>       fm3   = x - XI(3)
>       fm4   = x - XI(4)
1778c1783
<     !C  Get the "answer" at the requested X
---
> !C  Get the "answer" at the requested X
1780,1787c1785,1792
<     ans = fm2*fm3*fm4*YI(1)/&
<          ((XI(1)-XI(2))*(XI(1)-XI(3))*(XI(1)-XI(4)))  +&
<          fm1*fm3*fm4*YI(2)/&
<          ((XI(2)-XI(1))*(XI(2)-XI(3))*(XI(2)-XI(4)))  +&
<          fm1*fm2*fm4*YI(3)/&
<          ((XI(3)-XI(1))*(XI(3)-XI(2))*(XI(3)-XI(4)))  +&
<          fm1*fm2*fm3*YI(4)/&
<          ((XI(4)-XI(1))*(XI(4)-XI(2))*(XI(4)-XI(3)))
---
>       ans = fm2*fm3*fm4*YI(1)/&
>                       ((XI(1)-XI(2))*(XI(1)-XI(3))*(XI(1)-XI(4)))  +&
>             fm1*fm3*fm4*YI(2)/&
>                       ((XI(2)-XI(1))*(XI(2)-XI(3))*(XI(2)-XI(4)))  +&
>             fm1*fm2*fm4*YI(3)/&
>                       ((XI(3)-XI(1))*(XI(3)-XI(2))*(XI(3)-XI(4)))  +&
>             fm1*fm2*fm3*YI(4)/&
>                       ((XI(4)-XI(1))*(XI(4)-XI(2))*(XI(4)-XI(3)))
1790c1795
<   end subroutine lagrange
---
>       end subroutine lagrange
1793c1798
<   !##########################################################
---
> !##########################################################
1801c1806
<     return
---
>   return
1805c1810
<     return
---
>   return
Only in SWMF_yesterday/UA/MGITM/src: user.f90
diff -r -x .git SWMF/UA/MGITM/srcInterface/UA_wrapper.f90 SWMF_yesterday/UA/MGITM/srcInterface/UA_wrapper.f90
33,34d32
<   ! real:: XyzTest_D(3) = [401037.03, 4071804.0, 402979.53]
< 
39,40c37,38
<     use ModInputs, ONLY: cInputText
<     use ModReadParam, ONLY: read_text
---
>     use ModInputs, only: cInputText
>     use ModReadParam, only: read_text
43c41
<     use ModInputs, ONLY: iStartTime, IsFramework, iOutputUnit_, set_defaults, &
---
>     use ModInputs, only: iStartTime, IsFramework, iOutputUnit_, set_defaults, &
46c44
<     use CON_physics, ONLY: get_time
---
>     use CON_physics,    ONLY: get_time
49,50c47,48
<     use ModGITM, ONLY: iCommGITM, nProcs, iProcGITM => iProc
<     use GITM_planet, ONLY: init_planet
---
>     use ModGITM, only: iCommGITM, nProcs, iProcGITM => iProc
>     use GITM_planet, only: init_planet
59a58
>     write(*,*) "-->Starting UA_set_param..."
62c61,64
<        call put(CompInfo, Use=.true., NameVersion='MGITM')
---
>        call put(CompInfo,&
>             Use=.true.,                                      &
>             NameVersion='Global Iono-Thermo Model (Ridley)', &
>             Version=2.0)
105,107c107,109
<     use ModInputs, ONLY: LatStart, LatEnd, LonStart, LonEnd, AltMin, AltMax
<     use ModSizeGitm, ONLY: nLons, nLats, nAlts
<     use ModInputs, ONLY: nBlocksLat, nBlocksLon
---
>     use ModInputs,    ONLY: LatStart, LatEnd, LonStart, LonEnd, AltMin, AltMax
>     use ModSizeGitm,  ONLY: nLons, nLats, nAlts
>     use ModInputs,    ONLY: nBlocksLat, nBlocksLon
109c111
<     use ModGeometry, ONLY: TypeGeometry
---
>     use ModGeometry,  ONLY: TypeGeometry
132,133c134,135
<     use CON_physics, ONLY: get_time
<     use ModTime, ONLY: StartTime, iTimeArray, CurrentTime
---
>     use CON_physics,    ONLY: get_time
>     use ModTime, only : StartTime, iTimeArray, CurrentTime
166,167c168,169
<     use ModGITM, ONLY: iProc, Dt
<     use ModTime, ONLY: StartTime, CurrentTime, EndTime, iStep
---
>     use ModGITM,   ONLY: iProc, Dt
>     use ModTime,   ONLY: StartTime, CurrentTime, EndTime, iStep
244c246
<     logical:: DoTest, DoTestMe, DoTestPoint
---
>     logical:: DoTest, DoTestMe
256,259d257
<        ! if(sum(abs(Xyz_DI(:,iPoint) - XyzTest_D)) < 100.0) &
<        !     write(*,*) NameSub,': iPoint, Lon, Lat, iProc_I=', &
<        !     iPoint, Lon, Lat, iProc_I(iPoint)
< 
330c328
<     logical:: DoTest, DoTestMe ! , DoTestCell
---
>     logical:: DoTest, DoTestMe
352,355d349
< 
<           ! DoTestCell = sum(abs(Xyz_DI(:,iPoint) - XyzTest_D)) < 100.0
<           ! if(DoTestCell)write(*,*) NameSub,': iProc, iPoint=', iProc, iPoint
< 
357d350
< 
362,364d354
<           ! if(DoTestCell)write(*,*) NameSub,': rBody, Alt, Lon, Lat=', &
<           !     rBody, Alt, Lon, Lat
< 
377,379d366
< 
<              ! if(DoTestCell) write(*,*) NameSub, &
<              !     ': iBlock, iLon, iLat, nAlts=', iBlock, iLon, iLat, nAlts
394,396d380
< 
<              ! if(DoTestCell) write(*,*) NameSub, &
<              !     ': iBlock, iLon, iLat, iAlt=', iBlock, iLon, iLat, iAlt
401,403c385
<     ! Calculate neutral temperature in [K]
<     if(.not.allocated(Tneu_GB)) &
<          allocate(Tneu_GB(-1:nLons+2,-1:nLats+2,-1:nAlts+2,nBlocks))
---
>     if(.not.allocated(Tneu_GB)) allocate(Tneu_GB(-1:nLons+2,-1:nLats+2,-1:nAlts+2,nBlocks))
419,422d400
<        ! DoTestCell = sum(abs(Xyz_DI(:,iPoint) - XyzTest_D)) < 100.0
<        ! if(DoTestCell) write(*,*) NameSub, &
<        !    ' iProc, iPoint, iBlock, iCell_D=', iProc, iPoint, iBlock, iCell_D
< 
445,449d422
< 
<           ! if(DoTestCell)write(*,*) NameSub, &
<           !     ': NCO2, NCO2(iCell), dH, HscaleCO2=', Data_VI(2,iPoint), &
<           !     NDensityS(iCell_D(1), iCell_D(2), nAlts,iCO2_,iBlock), &
<           !     dH, HscaleCO2
