Changeset 2028 for trunk


Ignore:
Timestamp:
Oct 26, 2018, 8:25:32 PM (6 years ago)
Author:
flefevre
Message:

Diverses corrections dans la photolyse on-line (encore non operationnelle)

Location:
trunk/LMDZ.MARS/libf/aeronomars
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90

    r2027 r2028  
    9393real :: dt_guess            ! first-guess timestep (s)
    9494real :: dt_corrected        ! corrected timestep (s)
    95 real :: j(nlayer,nd)        ! interpolated photolysis rates (s-1)
    9695real :: time                ! internal time (between 0 and ptimestep, in s)
    9796
     
    147146if (jonline) then
    148147   tau = tau*press(1)/6.1  ! temporary
    149    call photolysis_online(nlayer,lswitch, alt, press, temp,                  &
    150                           zmmean, i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,  &
    151                           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,                  &
    152                           i_n, i_n2d, i_no, i_no2, i_n2, nesp, rm,           &
     148   call photolysis_online(nlayer, alt, press, temp, zmmean,          &
     149                          i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,  &
     150                          i_h2, i_oh, i_ho2, i_h2o2, i_h2o,          &
     151                          i_n, i_n2d, i_no, i_no2, i_n2, nesp, rm,   &
    153152                          tau, sza, dist_sol, v_phot)
    154153else
     
    18251824!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    18261825
    1827       do l = 1,lswitch-1
     1826      do l = 1,nlayer
    18281827         rm(l,i_co2)  = zycol(l, igcm_co2)
    18291828         rm(l,i_co)   = zycol(l, igcm_co)
  • trunk/LMDZ.MARS/libf/aeronomars/photolysis_online.F

    r2027 r2028  
    11!==============================================================================
    22
    3       subroutine photolysis_online(nlayer,lswitch, alt, press, temp,
     3      subroutine photolysis_online(nlayer, alt, press, temp,
    44     $           zmmean, i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,
    55     $           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               
     
    1414
    1515      integer :: nesp                                                ! total number of chemical species
    16       integer :: nlayer, lswitch
     16      integer :: nlayer
    1717      integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,           
    1818     $           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,                 
     
    3434!     integer, parameter :: nw = 3789                                ! number of spectral intervals (high-res)
    3535      integer, parameter :: nw = 152                                 ! number of spectral intervals (low-res)
    36       real, dimension(nw) :: wl, wc, wu                              ! lower, center, upper wavelength for each interval
    37       save wl, wc, wu
     36      real, dimension(nw), save :: wl, wc, wu                        ! lower, center, upper wavelength for each interval
    3837
    3938!     solar flux
    4039
    41       real, dimension(nw) :: f                                       !  solar flux (w.m-2.nm-1) at 1 au
     40      real, dimension(nw), save :: f                                 !  solar flux (w.m-2.nm-1) at 1 au
    4241      real :: factor
    43       save f
    4442
    4543!     cross-sections and quantum yields
    4644
    47       integer, parameter  :: nabs = 10                               ! number of absorbing gases
    48       real, dimension(nlayer,nw,nd) :: sj                            ! general cross-section array (cm2)
    49       real, dimension(nw) :: xsco2_195, xsco2_295, xsco2_370         ! co2 absorption cross-section at 195-295-370 k (cm2)
    50       real, dimension(nw) :: yieldco2                                ! co2 photodissociation yield
    51       real, dimension(nw) :: xso2_150, xso2_200, xso2_250, xso2_300  ! o2 absorption cross-section at 150-200-250-300 k (cm2)
    52       real, dimension(nw) :: yieldo2                                 ! o2 photodissociation yield
    53       real, dimension(nw) :: xso3_218, xso3_298                      ! o3 absorption cross-section at 218-298 k (cm2)
    54       real, dimension(nw) :: xsh2o                                   ! h2o absorption cross-section (cm2)
    55       real, dimension(nw) :: xsh2o2                                  ! h2o2 absorption cross-section (cm2)
    56       real, dimension(nw) :: xsho2                                   ! ho2 absorption cross-section (cm2)
    57       real, dimension(nw) :: xsh2                                    ! h2 absorption cross-section (cm2)
    58       real, dimension(nw) :: yieldh2                                 ! h2 photodissociation yield
    59       real, dimension(nw) :: xsno2, xsno2_220, xsno2_294             ! no2 absorption cross-section at 220-294 k (cm2)
    60       real, dimension(nw) :: yldno2_248, yldno2_298                  ! no2 quantum yield at 248-298 k
    61       real, dimension(nw) :: xsno                                    ! no absorption cross-section (cm2)
    62       real, dimension(nw) :: yieldno                                 ! no photodissociation yield
    63       real, dimension(nw) :: xsn2                                    ! n2 absorption cross-section (cm2)
    64       real, dimension(nw) :: yieldn2                                 ! n2 photodissociation yield
    65       real, dimension(nw) :: albedo                                  ! surface albedo
     45      integer, parameter  :: nabs  = 10                              ! number of absorbing gases
     46      integer, parameter  :: nphot = 13                              ! number of photolysis
     47      real, dimension(nlayer,nw,nphot) :: sj                         ! general cross-section array (cm2)
     48      real, dimension(nw), save :: xsco2_195, xsco2_295, xsco2_370   ! co2 absorption cross-section at 195-295-370 k (cm2)
     49      real, dimension(nw), save :: yieldco2                          ! co2 photodissociation yield
     50      real, dimension(nw), save :: xso2_150, xso2_200,               ! o2 absorption cross-section at 150-200-250-300 k (cm2)
     51     $                             xso2_250, xso2_300
     52      real, dimension(nw), save :: yieldo2                           ! o2 photodissociation yield
     53      real, dimension(nw), save :: xso3_218, xso3_298                ! o3 absorption cross-section at 218-298 k (cm2)
     54      real, dimension(nw), save :: xsh2o                             ! h2o absorption cross-section (cm2)
     55      real, dimension(nw), save :: xsh2o2                            ! h2o2 absorption cross-section (cm2)
     56      real, dimension(nw), save :: xsho2                             ! ho2 absorption cross-section (cm2)
     57      real, dimension(nw), save :: xsh2                              ! h2 absorption cross-section (cm2)
     58      real, dimension(nw), save :: yieldh2                           ! h2 photodissociation yield
     59      real, dimension(nw), save :: xsno2, xsno2_220, xsno2_294       ! no2 absorption cross-section at 220-294 k (cm2)
     60      real, dimension(nw), save :: yldno2_248, yldno2_298            ! no2 quantum yield at 248-298 k
     61      real, dimension(nw), save :: xsno                              ! no absorption cross-section (cm2)
     62      real, dimension(nw), save :: yieldno                           ! no photodissociation yield
     63      real, dimension(nw), save :: xsn2                              ! n2 absorption cross-section (cm2)
     64      real, dimension(nw), save :: yieldn2                           ! n2 photodissociation yield
     65      real, dimension(nw), save :: albedo                            ! surface albedo
    6666
    6767!     atmosphere
     
    8989      integer :: a_o2, a_co2, a_o3, a_h2o, a_h2o2, a_ho2, a_h2, a_no,
    9090     $           a_no2, a_n2
    91       integer :: i, ilay, iw, ialt, mopt
     91      integer :: i, ilay, iw, ialt
     92      integer, save :: mopt
    9293      real    :: deltaj
    9394
    94       logical, save :: firstcall
    95       data firstcall / .true. /
    96 
    97 !     high-res/low-res switch
    98 !
    99 !     mopt = 1 high-resolution
    100 !     mopt = 2 low-resolution (recommended for gcm use)
    101 
    102       if (nw >= 3789) then
    103          mopt = 1
    104       else
    105          mopt = 2
    106       end if
     95      logical, save :: firstcall = .true.
    10796
    10897!     absorbing gas numbering
     
    119108      a_n2      = 10     ! no2
    120109
    121 !     photodissociation rates numbering
     110!     photodissociation rates numbering.
     111!     photodissociations must be ordered the same way in subroutine "indice"
    122112
    123113      j_o2_o    = 1      ! o2 + hv     -> o + o
     
    143133      if (firstcall) then
    144134
     135!        high-res/low-res switch
     136!
     137!        mopt = 1 high-resolution
     138!        mopt = 2 low-resolution (recommended for gcm use)
     139
     140         mopt = 2
     141
    145142!        set wavelength grid
    146143
     
    208205!     o2:
    209206
    210       call seto2(nd, nlayer, nw, wc, mopt, temp, xso2_150, xso2_200,
     207      call seto2(nphot, nlayer, nw, wc, mopt, temp, xso2_150, xso2_200,
    211208     $           xso2_250, xso2_300, yieldo2, j_o2_o, j_o2_o1d,
    212209     $           colinc, rm(:,i_o2), dtgas(:,:,a_o2), sj)
     
    214211!     co2:
    215212
    216       call setco2(nd, nlayer, nw, wc, temp, xsco2_195, xsco2_295,
     213      call setco2(nphot, nlayer, nw, wc, temp, xsco2_195, xsco2_295,
    217214     $            xsco2_370, yieldco2, j_co2_o, j_co2_o1d,
    218215     $            colinc, rm(:,i_co2), dtgas(:,:,a_co2), sj)
     
    220217!     o3:
    221218
    222       call seto3(nd, nlayer, nw, wc, temp, xso3_218, xso3_298,
     219      call seto3(nphot, nlayer, nw, wc, temp, xso3_218, xso3_298,
    223220     $           j_o3_o, j_o3_o1d,
    224221     $           colinc, rm(:,i_o3), dtgas(:,:,a_o3), sj)
     
    226223!     h2o2:
    227224
    228       call seth2o2(nd, nlayer, nw, wc, temp, xsh2o2, j_h2o2,
     225      call seth2o2(nphot, nlayer, nw, wc, temp, xsh2o2, j_h2o2,
    229226     $             colinc, rm(:,i_h2o2), dtgas(:,:,a_h2o2), sj)
    230227
    231228!     no2:
    232229
    233       call setno2(nd, nlayer, nw, wc, temp, xsno2, xsno2_220,
     230      call setno2(nphot, nlayer, nw, wc, temp, xsno2, xsno2_220,
    234231     $            xsno2_294, yldno2_248, yldno2_298, j_no2,
    235232     $            colinc, rm(:,i_no2), dtgas(:,:,a_no2), sj)
     
    301298!==== initialise photolysis rates
    302299
    303       v_phot(:,1:nd) = 0.
     300      v_phot(:,1:nphot) = 0.
    304301
    305302!==== start of wavelength lopp
     
    320317!     spherical actinic flux
    321318
    322          do ilay = 1,lswitch-1
     319         do ilay = 1,nlayer
    323320            saflux(ilay) = f(iw)*(fdir(ilay) + fdn(ilay) + fup(ilay))
    324321         end do
     
    326323!     photolysis rate integration
    327324
    328          do i = 1,nd
    329             do ilay = 1, lswitch-1
     325         do i = 1,nphot
     326            do ilay = 1,nlayer
    330327               deltaj = saflux(ilay)*sj(ilay,iw,i)
    331328               v_phot(ilay,i) = v_phot(ilay,i) + deltaj*(wu(iw)-wl(iw))
     
    335332!     eliminate small values
    336333
    337          where (v_phot(:,1:nd) < 1.e-30)
    338             v_phot(:,1:nd) = 0.
     334         where (v_phot(:,1:nphot) < 1.e-30)
     335            v_phot(:,1:nphot) = 0.
    339336         end where
    340337
     
    343340      else   ! night
    344341
    345          v_phot(:,1:nd) = 0.
     342         v_phot(:,1:nphot) = 0.
    346343
    347344      end if ! sza
     
    945942!     local
    946943
    947       integer, parameter :: kdata = 36000
     944      integer, parameter :: kdata = 42000
    948945      real, parameter :: deltax = 1.e-4
    949946      real, dimension(kdata) :: x1, y1, y2, y3, xion, ion
     
    964961ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    965962c
    966 c     195K: Chan et al. (1993) + Starck et al. (2006)
     963c     195K: huestis and berkowitz (2010) + Starck et al. (2006)
    967964c           + Yoshino et al. (1996) + Parkinson et al. (2003) + extrapolation
    968 c           11 lignes d'en-tete et n1 = 34537
    969 c           fichier: co2_euv_uv_195k.txt
    970965c
    971 c     295K: Chan et al. (1993) + Starck et al. (2006)
     966c     295K: huestis and berkowitz (2010) + Starck et al. (2006)
    972967c           + Yoshino et al. (1996) + Parkinson et al. (2003) + extrapolation
    973 c           11 lignes d'en-tete et n2 = 35360
    974 c           fichier: co2_euv_uv_295k.txt
    975968c
    976 c     370K: Chan et al. (1993) + Starck et al. (2006)
     969c     370K: huestis and berkowitz (2010) + Starck et al. (2006)
    977970c           + Lewis and Carver (1983) + extrapolation
    978 c           11 lignes d'en-tete et n3 = 3884
    979 c           fichier: co2_euv_uv_370k.txt
    980971c
    981972ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    982 c
    983       n1 = 34537
    984       n2 = 35360
    985       n3 = 3884
    986 c
    987 c     195K:
    988 c
    989       fil = trim(datadir)//'/cross_sections/co2_euv_uv_195k.txt'
     973
     974      n1 = 40769
     975      n2 = 41586
     976      n3 = 10110
     977
     978!     195K:
     979
     980      fil = trim(datadir)//'/cross_sections/co2_euv_uv_2018_195k.txt'
    990981      print*, 'section efficace CO2 195K: ', fil
    991982
     
    10131004         xsco2_195(l) = yg(l)
    10141005      END DO
    1015 c
    1016 c     295K:
    1017 c
    1018       fil = trim(datadir)//'/cross_sections/co2_euv_uv_295k.txt'
     1006
     1007!     295K:
     1008
     1009      fil = trim(datadir)//'/cross_sections/co2_euv_uv_2018_295k.txt'
    10191010      print*, 'section efficace CO2 295K: ', fil
    1020 c
     1011
    10211012      OPEN(UNIT=kin,FILE=fil,STATUS='old')
    10221013      DO i = 1,11
    10231014         read(kin,*)
    10241015      END DO
    1025 c
     1016
    10261017      DO i = 1, n2
    10271018         READ(kin,*) x1(i), y1(i)
     
    10381029         STOP
    10391030      ENDIF
    1040 c
     1031
    10411032      DO l = 1, nw-1
    10421033         xsco2_295(l) = yg(l)
    10431034      END DO
    1044 c
    1045 c     370K:
    1046 c
    1047       fil = trim(datadir)//'/cross_sections/co2_euv_uv_370k.txt'
     1035
     1036!     370K:
     1037
     1038      fil = trim(datadir)//'/cross_sections/co2_euv_uv_2018_370k.txt'
    10481039      print*, 'section efficace CO2 370K: ', fil
    1049 c
     1040
    10501041      OPEN(UNIT=kin,FILE=fil,STATUS='old')
    10511042      DO i = 1,11
    10521043         read(kin,*)
    10531044      END DO
    1054 c
     1045
    10551046      DO i = 1, n3
    10561047         READ(kin,*) x1(i), y1(i)
    10571048      END DO
    10581049      CLOSE (kin)
    1059 c
     1050
    10601051      CALL addpnt(x1,y1,kdata,n3,x1(1)*(1.-deltax),0.)
    10611052      CALL addpnt(x1,y1,kdata,n3,          0.,0.)
Note: See TracChangeset for help on using the changeset viewer.