Ignore:
Timestamp:
Apr 26, 2012, 3:22:19 PM (13 years ago)
Author:
emillour
Message:

Mars GCM: Update of the chemistry package, including:

  • 93 reactions are accounted for (instead of 22); tracking 28 species (instead of 11)
  • computation of photoabsorption using raytracing
  • improved time stepping in the photochemistry
  • updated parameters (cross-sections); with this new version input files

are in 'EUV/param_v5' of "datafile" directory.

  • transition between lower and upper atmosphere chemistry set to 0.1 Pa (calchim.F90)
  • Lots of code clean-up: removed obsolete files column.F, param_v3.h, flujo.F, phdisrate.F, ch.F, interpfast.F, paramfoto.F, getch.F Converted chemtermos.F -> chemthermos.F90 and euvheat.F -> euvheat.F90. Added paramfoto_compact.F , param_v4.h and iono.h
  • Upadted surfacearea.F
  • Cleaned initracer.F and callkeys.h (removed use of obsolete "nqchem" and "oldnames" case when initializing tracers).
  • Minor correction in "callsedim": compute "rdust" and/or "rice" only when it makes sense.

FGG+FL+EM

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/aeronomars/calchim.F90

    r618 r635  
    9898      integer :: ig,l,i,iq,iqmax
    9999      integer :: foundswitch, lswitch
     100      integer,save :: chemthermod
    100101
    101102      integer,save :: i_co2  = 0
    … …  
    113114      integer,save :: i_n2   = 0
    114115      integer,save :: i_h2o  = 0
     116      integer,save :: i_n    = 0
     117      integer,save :: i_no   = 0
     118      integer,save :: i_no2  = 0
     119      integer,save :: i_n2d  = 0
     120      integer,save :: i_co2plus=0
     121      integer,save :: i_oplus=0
     122      integer,save :: i_o2plus=0
     123      integer,save :: i_coplus=0
     124      integer,save :: i_cplus=0
     125      integer,save :: i_nplus=0
     126      integer,save :: i_noplus=0
     127      integer,save :: i_n2plus=0
     128      integer,save :: i_hplus=0
     129      integer,save :: i_hco2plus=0
     130      integer,save :: i_elec=0
    115131
    116132      integer :: ig_vl1
    … …  
    157173
    158174         ! find index of chemical tracers to use
     175         ! Listed here are all tracers that can go into photochemistry
    159176         nbq = 0 ! to count number of tracers
     177         ! Species ALWAYS present if photochem=.T. or thermochem=.T.
    160178         i_co2 = igcm_co2
    161179         if (i_co2 == 0) then
    … …  
    248266         i_ch4 = igcm_ch4
    249267         if (i_ch4 == 0) then
    250             write(*,*) "calchim: no CH4 tracer !!!"
     268            write(*,*) "calchim: Error; no CH4 tracer !!!"
    251269            write(*,*) "CH4 will be ignored in the chemistry"
    252270         else
    … …  
    270288            niq(nbq) = i_h2o
    271289         end if
     290         !Check tracers needed for thermospheric chemistry
     291         if(thermochem) then
     292            chemthermod=0  !Default: C/O/H chemistry
     293            !Nitrogen chemistry
     294            !NO is used to determine if N chemistry is wanted
     295            !chemthermod=2 -> N chemistry
     296            i_no = igcm_no
     297            if (i_no == 0) then
     298               write(*,*) "calchim: no NO tracer"
     299               write(*,*) "C/O/H themosp chemistry only "
     300            else
     301               nbq = nbq + 1
     302               niq(nbq) = i_no
     303               chemthermod=2
     304               write(*,*) "calchim: NO in traceur.def"
     305               write(*,*) "Nitrogen chemistry included"
     306            end if
     307            ! N
     308            i_n = igcm_n
     309            if(chemthermod == 2) then
     310               if (i_n == 0) then
     311                  write(*,*) "calchim: Error; no N tracer !!!"
     312                  write(*,*) "N is needed if NO is in traceur.def"
     313                  stop
     314               else
     315                  nbq = nbq + 1
     316                  niq(nbq) = i_n
     317               end if
     318            else
     319               if (i_n /= 0) then
     320                  write(*,*) "calchim: Error: N is present, but NO is not!!!"
     321                  write(*,*) "Both must be in traceur.def if N chemistry is wanted"
     322                  stop
     323               endif
     324            endif    !Of if(chemthermod == 2)
     325            ! NO2
     326            i_no2 = igcm_no2
     327            if(chemthermod == 2) then
     328               if (i_no2 == 0) then
     329                  write(*,*) "calchim: Error; no NO2 tracer !!!"
     330                  write(*,*) "NO2 is needed if NO is in traceur.def"
     331                  stop
     332               else
     333                  nbq = nbq + 1
     334                  niq(nbq) = i_no2
     335               end if
     336            else
     337               if (i_no2 /= 0) then
     338                  write(*,*) "calchim: Error: N is present, but NO is not!!!"
     339                  write(*,*) "Both must be in traceur.def if N chemistry is wanted"
     340                  stop
     341               endif
     342            endif     !Of if(chemthermod == 2)
     343            ! N(2D)
     344            if(chemthermod == 2) then
     345               i_n2d = igcm_n2d
     346               if (i_n2d == 0) then
     347                  write(*,*) "calchim: Error; no N2D !!!"
     348                  write(*,*) "N2D is needed if NO is in traceur.def"
     349                  stop
     350               else
     351                  nbq = nbq + 1
     352                  niq(nbq) = i_n2d
     353               end if
     354            else
     355               if (i_n2d /= 0) then
     356                  write(*,*) "calchim: Error: N2D is present, but NO is not!!!"
     357                  write(*,*) "Both must be in traceur.def if N chemistry wanted"
     358                  stop
     359               endif
     360            endif    !Of if(chemthermod == 2)
     361            ! Ions
     362            ! O2+ is used to determine if ion chemistry is needed
     363            ! chemthermod=3 -> ion chemistry
     364            i_o2plus = igcm_o2plus
     365            if(chemthermod == 2) then
     366               if (i_o2plus == 0) then
     367                  write(*,*) "calchim: no O2+ tracer; no ion chemistry"
     368               else
     369                  nbq = nbq + 1
     370                  niq(nbq) = i_o2plus
     371                  chemthermod = 3
     372                  write(*,*) "calchim: O2+ in traceur.def"
     373                  write(*,*) "Ion chemistry included"
     374               end if
     375            else
     376               if (i_o2plus /= 0) then
     377                  write(*,*) "calchim: O2+ is present, but NO is not!!!"
     378                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
     379                  stop
     380               endif
     381            endif   !Of if(chemthermod == 2)
     382            ! CO2+
     383            i_co2plus = igcm_co2plus
     384            if(chemthermod == 3) then
     385               if (i_co2plus == 0) then
     386                  write(*,*) "calchim: Error; no CO2+ tracer !!!"
     387                  write(*,*) "CO2+ is needed if O2+ is in traceur.def"
     388                  stop
     389               else
     390                  nbq = nbq + 1
     391                  niq(nbq) = i_co2plus
     392               end if
     393            else
     394               if (i_co2plus /= 0) then
     395                  write(*,*) "calchim: Error: CO2+ is present, but O2+ is not!!!"
     396                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
     397                  stop
     398               endif
     399            endif    !Of if(chemthermod == 3)
     400            ! O+
     401            i_oplus = igcm_oplus
     402            if(chemthermod == 3) then
     403               if (i_oplus == 0) then
     404                  write(*,*) "calchim: Error; no O+ tracer !!!"
     405                  write(*,*) "O+ is needed if O2+ is in traceur.def"
     406                  stop
     407               else
     408                  nbq = nbq + 1
     409                  niq(nbq) = i_oplus
     410               end if
     411            else
     412               if (i_oplus /= 0) then
     413                  write(*,*) "calchim: Error: O+ is present, but O2+ is not!!!"
     414                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
     415                  stop
     416               endif
     417            endif   !Of if (chemthermod == 3)
     418            ! CO+
     419            i_coplus = igcm_coplus
     420            if(chemthermod == 3) then
     421               if (i_coplus == 0) then
     422                  write(*,*) "calchim: Error; no CO+ tracer !!!"
     423                  write(*,*) "CO+ is needed if O2+ is in traceur.def"
     424                  stop
     425               else
     426                  nbq = nbq + 1
     427                  niq(nbq) = i_coplus
     428               end if
     429            else
     430               if (i_coplus /= 0) then
     431                  write(*,*) "calchim: Error: CO+ is present, but O2+ is not!!!"
     432                  write(*,*) " Both must be in traceur.def if ionosphere wanted"
     433                  stop
     434               endif
     435            endif   ! Of if (chemthermod == 3)
     436            ! C+
     437            i_cplus = igcm_cplus
     438            if(chemthermod == 3) then
     439               if (i_cplus == 0) then
     440                  write(*,*) "calchim: Error; no C+ tracer !!!"
     441                  write(*,*) "C+ is needed if O2+ is in traceur.def"
     442                  stop
     443               else
     444                  nbq = nbq + 1
     445                  niq(nbq) = i_cplus
     446               end if
     447            else
     448               if (i_cplus /= 0) then
     449                  write(*,*) "calchim: Error; C+ is present, but O2+ is not!!!"
     450                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
     451                  stop
     452               endif
     453            endif   ! Of if (chemthermod == 3)
     454            ! N+
     455            i_nplus = igcm_nplus
     456            if(chemthermod == 3) then
     457               if (i_nplus == 0) then
     458                  write(*,*) "calchim: Error; no N+ tracer !!!"
     459                  write(*,*) "N+ is needed if O2+ is in traceur.def"
     460                  stop
     461               else
     462                  nbq = nbq + 1
     463                  niq(nbq) = i_nplus
     464               end if
     465            else
     466               if (i_nplus /= 0) then
     467                  write(*,*) "calchim: Error: N+ is present, but O2+ is not!!!"
     468                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
     469                  stop
     470               endif
     471            endif   !Of if (chemthermod == 3)
     472            ! NO+
     473            i_noplus = igcm_noplus
     474            if(chemthermod == 3) then
     475               if (i_noplus == 0) then
     476                  write(*,*) "calchim: Error; no NO+ tracer !!!"
     477                  write(*,*) "NO+ is needed if O2+ is in traceur.def"
     478                  stop
     479               else
     480                  nbq = nbq + 1
     481                  niq(nbq) = i_noplus
     482               end if
     483            else
     484               if (i_noplus /= 0) then
     485                  write(*,*) "calchim: Error: NO+ is present, but O2+ is not!!!"
     486                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
     487                  stop
     488               endif
     489            endif   !Of if (chemthermod == 3)
     490            ! N2+
     491            i_n2plus = igcm_n2plus
     492            if (chemthermod == 3) then
     493               if (i_n2plus == 0) then
     494                  write(*,*) "calchim: Error; no N2+ tracer !!!"
     495                  write(*,*) "N2+ is needed if O2+ is in traceur.def"
     496                  stop
     497               else
     498                  nbq = nbq + 1
     499                  niq(nbq) = i_n2plus
     500               end if
     501            else
     502               if (i_n2plus /= 0) then
     503                  write(*,*) "calchim: Error: N2+ is present, but O2+ is not!!!"
     504                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
     505                  stop
     506               endif
     507            endif   !Of if (chemthermod == 3)
     508            !H+
     509            i_hplus = igcm_hplus
     510            if (chemthermod == 3) then
     511               if (i_hplus == 0) then
     512                  write(*,*) "calchim: Error; no H+ tracer !!!"
     513                  write(*,*) "H+ is needed if O2+ is in traceur.def"
     514                  stop
     515               else
     516                  nbq = nbq + 1
     517                  niq(nbq) = i_hplus
     518               end if
     519            else
     520               if (i_hplus /= 0) then
     521                  write(*,*) "calchim: Error: H+ is present, but O2+ is not!!!"
     522                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
     523                  stop
     524               endif
     525            endif   !Of if (chemthermod == 3)
     526            ! HCO2+
     527            i_hco2plus = igcm_hco2plus
     528            if(chemthermod == 3) then
     529               if (i_hco2plus == 0) then
     530                  write(*,*) "calchim: Error; no HCO2+ tracer !!!"
     531                  write(*,*) "HCO2+ is needed if O2+ is in traceur.def"
     532                  stop
     533               else
     534                  nbq = nbq + 1
     535                  niq(nbq) = i_hco2plus
     536               end if
     537            else
     538               if (i_hco2plus /= 0) then
     539                  write(*,*) "calchim: Error: HCO2+ is present, but O2+ is not!!!"
     540                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
     541                  stop
     542               endif
     543            endif    !Of if(chemthermod == 3)
     544            !e-
     545            i_elec = igcm_elec
     546            if(chemthermod == 3) then
     547               if (i_elec == 0) then
     548                  write(*,*) "calchim: Error; no e- tracer !!!"
     549                  write(*,*) "e- is needed if O2+ is in traceur.def"
     550                  stop
     551               else
     552                  nbq = nbq + 1
     553                  niq(nbq) = i_elec
     554               end if
     555            else
     556               if(i_elec /= 0) then
     557                  write(*,*) "calchim: Error: e- is present, but O2+ is not!!!"
     558                  write(*,*) "Both must be in traceur.def if ionosphere wanted"
     559                  stop
     560               endif
     561            endif   !Of if (chemthermod == 3)
     562         endif      !Of thermochem
    272563
    273564         write(*,*) 'calchim: found nbq    = ',nbq,' tracers'
    274          
     565              
    275566         firstcall = .false.
    276567      end if ! if (firstcall)
    … …  
    279570
    280571      zycol(:,:)    = 0.
    281       dqchim(:,:,:) = 0
    282       dqschim(:,:)  = 0
     572      dqchim(:,:,:) = 0.
     573      dqschim(:,:)  = 0.
    283574
    284575!     latvl1= 22.27
    … …  
    292583
    293584      do ig = 1,ngridmx
    294 
     585         
    295586         foundswitch = 0
    296587         do l = 1,nlayermx
    … …  
    316607
    317608            if (photochem .and. thermochem) then
    318                if (foundswitch == 0 .and. pplay(ig,l) < 1.e-3) then
     609               if (foundswitch == 0 .and. pplay(ig,l) < 1.e-1) then
    319610                  lswitch = l
    320611                  foundswitch = 1
    … …  
    357648
    358649!        chemistry in upper atmosphere
    359 
     650       
    360651         if (thermochem) then
    361             call chemthermos(ig,lswitch,zycol,ztemp,zdens,zpress,  &
    362                              zlocal,szacol,ptimestep,zday)
     652            call chemthermos(ig,lswitch,chemthermod,zycol,ztemp,zdens,  &
     653                             zpress,zlocal,szacol,ptimestep,zday)
    363654         end if
    364655
    … …  
    392683            end do
    393684         end do ! of do l = 1,nlayermx
     685!         if(ig.eq.800)write(*,*)'calchim/686',dqchim(ig,27,23)
    394686
    395687!=======================================================================
    … …  
    419711      return
    420712      end
     713
Note: See TracChangeset for help on using the changeset viewer.