Changeset 2968


Ignore:
Timestamp:
May 24, 2023, 8:26:35 AM (19 months ago)
Author:
flefevre
Message:

on-line photolysis: bug fix in the vertical grid used by the UV code at the top of the atmosphere. This correction leads to significant changes in the chemistry with low top level PCM configurations. No signicant change is observed with PCM configurations including the thermosphere.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/aeronomars/photolysis_online.F

    r2615 r2968  
    22
    33      subroutine photolysis_online(nlayer, deutchem, nb_phot_max,
    4      $           alt, press, temp, zmmean,
     4     $           alt, press, temp, mmean,
    55     $           i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,
    66     $           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               
    … …  
    2222     $                       i_n, i_n2d, i_no, i_no2, i_n2
    2323
    24       real, dimension(nlayer), intent(in) :: press, temp, zmmean     ! pressure (hpa)/temperature (k)/mean molecular mass (g.mol-1)
     24      real, dimension(nlayer), intent(in) :: press, temp, mmean      ! pressure (hpa)/temperature (k)/mean molecular mass (g.mol-1)
    2525      real, dimension(nlayer), intent(in) :: alt                     ! altitude (km)
    2626      real, dimension(nlayer,nesp), intent(in) :: rm                 ! mixing ratios
    … …  
    4444!     atmosphere
    4545
    46       real, dimension(nlayer) :: colinc                              ! air column increment (molecule.cm-2)
    47       real, dimension(nlayer,nw) :: dtrl                             ! rayleigh optical depth
    48       real, dimension(nlayer,nw) :: dtaer                            ! aerosol optical depth
    49       real, dimension(nlayer,nw) :: omaer                            ! aerosol single scattering albedo
    50       real, dimension(nlayer,nw) :: gaer                             ! aerosol asymmetry parameter
    51       real, dimension(nlayer,nw) :: dtcld                            ! cloud optical depth
    52       real, dimension(nlayer,nw) :: omcld                            ! cloud single scattering albedo
    53       real, dimension(nlayer,nw) :: gcld                             ! cloud asymmetry parameter
    54       real, dimension(nlayer,nw,nabs) :: dtgas                       ! optical depth for each gas
    55       real, dimension(nlayer,nw) :: dagas                            ! total gas optical depth
    56       real, dimension(nlayer) :: edir, edn, eup                      ! normalised irradiances
    57       real, dimension(nlayer) :: fdir, fdn, fup                      ! normalised actinic fluxes
    58       real, dimension(nlayer) :: saflux                              ! total actinic flux
    59 
    60       integer, dimension(0:nlayer) :: nid
    61       real, dimension(0:nlayer,nlayer) :: dsdh
     46      real, dimension(nlayer+1) :: zpress, zalt, ztemp, zmmean       ! pressure (hpa)/altitude (km)/temperature (k)/mean molecular mass (g.mol-1)
     47
     48      real, dimension(nlayer+1) :: colinc                            ! air column increment (molecule.cm-2)
     49      real, dimension(nlayer+1,nw) :: dtrl                           ! rayleigh optical depth
     50      real, dimension(nlayer+1,nw) :: dtaer                          ! aerosol optical depth
     51      real, dimension(nlayer+1,nw) :: omaer                          ! aerosol single scattering albedo
     52      real, dimension(nlayer+1,nw) :: gaer                           ! aerosol asymmetry parameter
     53      real, dimension(nlayer+1,nw) :: dtcld                          ! cloud optical depth
     54      real, dimension(nlayer+1,nw) :: omcld                          ! cloud single scattering albedo
     55      real, dimension(nlayer+1,nw) :: gcld                           ! cloud asymmetry parameter
     56      real, dimension(nlayer+1,nw,nabs) :: dtgas                     ! optical depth for each gas
     57      real, dimension(nlayer+1,nw) :: dagas                          ! total gas optical depth
     58      real, dimension(nlayer+1) :: edir, edn, eup                    ! normalised irradiances
     59      real, dimension(nlayer+1) :: fdir, fdn, fup                    ! normalised actinic fluxes
     60      real, dimension(nlayer+1) :: saflux                            ! total actinic flux
     61
     62      integer, dimension(0:nlayer+1) :: nid
     63      real, dimension(0:nlayer+1,nlayer+1) :: dsdh
    6264
    6365      integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d, j_o3_o1d, j_o3_o,
    … …  
    6769      integer :: a_o2, a_co2, a_o3, a_h2o, a_h2o2, a_ho2, a_h2, a_no,
    6870     $           a_no2, a_n2
    69       integer :: i, ilay, iw, ialt
     71      integer :: nlev, i, ilay, ilev, iw, ialt
    7072      real    :: deltaj
    7173
    … …  
    102104      j_hdo_d   = 15     ! hdo + hv    -> d + oh
    103105
     106!==== define working vertical grid for the uv radiative code
     107
     108      nlev = nlayer + 1
     109
     110      do ilev = 1,nlev-1
     111         zpress(ilev) = press(ilev)
     112         zalt(ilev)   = alt(ilev)
     113         ztemp(ilev)  = temp(ilev)
     114         zmmean(ilev) = mmean(ilev)
     115      end do
     116
     117      zpress(nlev) = 0.                  ! top of the atmosphere
     118      zalt(nlev)   = zalt(nlev-1) + (zalt(nlev-1) - zalt(nlev-2))
     119      ztemp(nlev)  = ztemp(nlev-1)
     120      zmmean(nlev) = zmmean(nlev-1)
     121     
    104122!==== air column increments and rayleigh optical depth
    105123
    106       call setair(nlayer, nw, wl, wc, press, temp, zmmean, colinc, dtrl)
     124      call setair(nlev, nw, wl, wc, zpress, ztemp, zmmean, colinc, dtrl)
    107125
    108126!==== set temperature-dependent cross-sections and optical depths
     127
     128      dtgas(:,:,:) = 0.
    109129
    110130!     o2:
    … …  
    112132      call seto2(nphot, nlayer, nw, wc, mopt, temp, xso2_150, xso2_200,
    113133     $           xso2_250, xso2_300, yieldo2, j_o2_o, j_o2_o1d,
    114      $           colinc, rm(:,i_o2), dtgas(:,:,a_o2), sj)
     134     $           colinc(1:nlayer), rm(:,i_o2),
     135     $           dtgas(1:nlayer,:,a_o2), sj)
    115136
    116137!     co2:
    … …  
    118139      call setco2(nphot, nlayer, nw, wc, temp, xsco2_195, xsco2_295,
    119140     $            xsco2_370, yieldco2, j_co2_o, j_co2_o1d,
    120      $            colinc, rm(:,i_co2), dtgas(:,:,a_co2), sj)
     141     $            colinc(1:nlayer), rm(:,i_co2),
     142     $            dtgas(1:nlayer,:,a_co2), sj)
    121143
    122144!     o3:
    … …  
    124146      call seto3(nphot, nlayer, nw, wc, temp, xso3_218, xso3_298,
    125147     $           j_o3_o, j_o3_o1d,
    126      $           colinc, rm(:,i_o3), dtgas(:,:,a_o3), sj)
     148     $           colinc(1:nlayer), rm(:,i_o3),
     149     $           dtgas(1:nlayer,:,a_o3), sj)
    127150
    128151!     h2o2:
    129152
    130153      call seth2o2(nphot, nlayer, nw, wc, temp, xsh2o2, j_h2o2,
    131      $             colinc, rm(:,i_h2o2), dtgas(:,:,a_h2o2), sj)
     154     $             colinc(1:nlayer), rm(:,i_h2o2),
     155     $             dtgas(1:nlayer,:,a_h2o2), sj)
    132156
    133157!     no2:
    … …  
    135159      call setno2(nphot, nlayer, nw, wc, temp, xsno2, xsno2_220,
    136160     $            xsno2_294, yldno2_248, yldno2_298, j_no2,
    137      $            colinc, rm(:,i_no2), dtgas(:,:,a_no2), sj)
    138 
    139 !==== temperature independent optical depths and cross-sections
     161     $            colinc(1:nlayer), rm(:,i_no2),
     162     $            dtgas(1:nlayer,:,a_no2), sj)
     163
     164!==== temperature-independent optical depths and cross-sections
    140165
    141166!     optical depths
    … …  
    175200      end do
    176201
    177       !HDO cross section
     202!     if deuterium chemistry: hdo cross-section
     203
    178204      if (deutchem) then
    179          do ilay=1,nlayer
    180             do iw=1,nw-1
    181                !Two chanels for HDO, with same cross section
     205         do ilay = 1,nlayer
     206            do iw = 1,nw-1
     207               ! two chanels for hdo, with same cross-section
    182208               sj(ilay,iw,j_hdo_od) = 0.5*xshdo(iw) ! hdo -> od + h
    183209               sj(ilay,iw,j_hdo_d)  = 0.5*xshdo(iw) ! hdo -> d + oh
    184             enddo
    185          enddo
    186       endif
     210            end do
     211         end do
     212      end if
    187213
    188214!==== set aerosol properties and optical depth
    189215
    190       call setaer(nlayer,alt,tau,nw,dtaer,omaer,gaer)
     216      call setaer(nlev,zalt,tau,nw,dtaer,omaer,gaer)
    191217
    192218!==== set cloud properties and optical depth
    193219
    194       call setcld(nlayer,nw,dtcld,omcld,gcld)
     220      call setcld(nlev,nw,dtcld,omcld,gcld)
    195221
    196222!==== slant path lengths in spherical geometry
    197223
    198       call sphers(nlayer,alt,sza,dsdh,nid)
     224      call sphers(nlev,zalt,sza,dsdh,nid)
    199225
    200226!==== solar flux at mars
    … …  
    214240
    215241!     monochromatic radiative transfer. outputs are:
    216 !     normalized irradiances     edir(nlayer), edn(nlayer), eup(nlayer)
    217 !     normalized actinic fluxes  fdir(nlayer), fdn(nlayer), fup(nlayer)
     242!     normalized irradiances     edir(nlev), edn(nlev), eup(nlev)
     243!     normalized actinic fluxes  fdir(nlev), fdn(nlev), fup(nlev)
    218244!     where
    219245!     dir = direct beam, dn = down-welling diffuse, up = up-welling diffuse
    220246
    221          call rtlink(nlayer, nw, iw, albedo(iw), sza, dsdh, nid, dtrl,
     247         call rtlink(nlev, nw, iw, albedo(iw), sza, dsdh, nid, dtrl,
    222248     $               dagas, dtcld, omcld, gcld, dtaer, omaer, gaer,
    223249     $               edir, edn, eup, fdir, fdn, fup)
    … …  
    287313         colinc(ilev) = avo*0.1*dp/(zmmean(ilev)*g)
    288314      end do
    289       colinc(nlev) = avo*0.1*press(nlev)*100./(zmmean(nlev)*g)
     315      colinc(nlev) = 0.
    290316
    291317      do iw = 1, nw - 1
    … …  
    711737!==============================================================================
    712738
    713       subroutine setaer(nlayer,alt,tau,nw,dtaer,omaer,gaer)
     739      subroutine setaer(nlev,zalt,tau,nw,dtaer,omaer,gaer)
    714740
    715741!-----------------------------------------------------------------------------*
    … …  
    723749!     input
    724750
    725       integer :: nlayer                              ! number of vertical layers
     751      integer :: nlev                                ! number of vertical levels
    726752      integer :: nw                                  ! number of wavelength grid points
    727       real, dimension(nlayer) :: alt                 ! altitude (km)
     753      real, dimension(nlev) :: zalt                  ! altitude (km)
    728754      real :: tau                                    ! integrated aerosol optical depth at the surface
    729755
    730756!     output
    731757
    732       real, dimension(nlayer,nw) :: dtaer            ! aerosol optical depth
    733       real, dimension(nlayer,nw) :: omaer            ! aerosol single scattering albedo
    734       real, dimension(nlayer,nw) :: gaer             ! aerosol asymmetry parameter
     758      real, dimension(nlev,nw) :: dtaer              ! aerosol optical depth
     759      real, dimension(nlev,nw) :: omaer              ! aerosol single scattering albedo
     760      real, dimension(nlev,nw) :: gaer               ! aerosol asymmetry parameter
    735761
    736762!     local
    737763
    738       integer :: ilay, iw
    739       real, dimension(nlayer) :: aer                                 ! dust extinction
     764      integer :: ilev, iw
     765      real, dimension(nlev) :: aer                   ! dust extinction
    740766      real :: omega, g, scaleh, gamma
    741767      real :: dz, tautot, q0
    … …  
    753779
    754780      tautot = 0.
    755       do ilay = 1, nlayer-1
    756          dz = alt(ilay+1) - alt(ilay)
    757          tautot = tautot + exp(gamma*(1. - exp(alt(ilay)/scaleh)))*dz
     781      do ilev = 1, nlev-1
     782         dz = zalt(ilev+1) - zalt(ilev)
     783         tautot = tautot + exp(gamma*(1. - exp(zalt(ilev)/scaleh)))*dz
    758784      end do
    759785     
    760786      q0 = tau/tautot
    761       do ilay = 1, nlayer-1
    762          dz = alt(ilay+1) - alt(ilay)
    763          dtaer(ilay,:) = q0*exp(gamma*(1. - exp(alt(ilay)/scaleh)))*dz
    764          omaer(ilay,:) = omega
    765          gaer(ilay,:)  = g
     787      do ilev = 1, nlev-1
     788         dz = zalt(ilev+1) - zalt(ilev)
     789         dtaer(ilev,:) = q0*exp(gamma*(1. - exp(zalt(ilev)/scaleh)))*dz
     790         omaer(ilev,:) = omega
     791         gaer(ilev,:)  = g
    766792      end do
    767793
    … …  
    770796!==============================================================================
    771797
    772       subroutine setcld(nlayer,nw,dtcld,omcld,gcld)
     798      subroutine setcld(nlev,nw,dtcld,omcld,gcld)
    773799
    774800!-----------------------------------------------------------------------------*
    … …  
    782808!     input
    783809
    784       integer :: nlayer                              ! number of vertical layers
     810      integer :: nlev                                ! number of vertical levels
    785811      integer :: nw                                  ! number of wavelength grid points
    786812
    787813!     output
    788814
    789       real, dimension(nlayer,nw) :: dtcld            ! cloud optical depth
    790       real, dimension(nlayer,nw) :: omcld            ! cloud single scattering albedo
    791       real, dimension(nlayer,nw) :: gcld             ! cloud asymmetry parameter
     815      real, dimension(nlev,nw) :: dtcld            ! cloud optical depth
     816      real, dimension(nlev,nw) :: omcld            ! cloud single scattering albedo
     817      real, dimension(nlev,nw) :: gcld             ! cloud asymmetry parameter
    792818
    793819!     local
    794820
    795       integer :: ilay, iw
     821      integer :: ilev, iw
    796822
    797823!     dtcld : optical depth
    … …  
    799825!     gcld  : asymmetry factor
    800826
    801       do ilay = 1, nlayer - 1
    802          do iw = 1, nw - 1
    803             dtcld(ilay,iw) = 0.       ! no clouds for the moment
    804             omcld(ilay,iw) = 0.99
    805             gcld(ilay,iw)  = 0.85
     827      dtcld(:,:) = 0.
     828
     829      do ilev = 1,nlev-1
     830         do iw = 1,nw-1
     831            dtcld(ilev,iw) = 0.       ! no clouds for the moment
     832            omcld(ilev,iw) = 0.99
     833            gcld(ilev,iw)  = 0.85
    806834         end do
    807835      end do
Note: See TracChangeset for help on using the changeset viewer.