Changeset 2974 for trunk/LMDZ.VENUS


Ignore:
Timestamp:
May 30, 2023, 10:56:31 AM (18 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.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/LMDZ.VENUS/libf/phyvenus/photolysis_online.F

    r2925 r2974  
    22
    33      subroutine photolysis_online(nlayer, 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_h2,
    66     $           i_oh, i_ho2, i_h2o2, i_h2o,i_h,i_hcl,
     
    2525     $                       i_no2,i_no, i_n2, i_n2d, i_h2
    2626
    27       real, dimension(nlayer), intent(in) :: press, temp,  zmmean    ! pressure (hpa)/temperature (k)/mean molecular mass (g.mol-1)
     27      real, dimension(nlayer), intent(in) :: press, temp, mmean      ! pressure (hpa)/temperature (k)/mean molecular mass (g.mol-1)
    2828      real, dimension(nlayer), intent(in) :: alt                     ! altitude (km)
    2929      real, dimension(nlayer,nesp), intent(in) :: rm                 ! mixing ratios
     
    4747!     atmosphere
    4848
    49       real, dimension(nlayer) :: colinc                              ! air column increment (molecule.cm-2)
    50       real, dimension(nlayer,nw) :: dtrl                             ! rayleigh optical depth
    51       real, dimension(nlayer,nw) :: dtaer                            ! aerosol optical depth
    52       real, dimension(nlayer,nw) :: omaer                            ! aerosol single scattering albedo
    53       real, dimension(nlayer,nw) :: gaer                             ! aerosol asymmetry parameter
    54       real, dimension(nlayer,nw) :: dtcld                            ! cloud optical depth
    55       real, dimension(nlayer,nw) :: omcld                            ! cloud single scattering albedo
    56       real, dimension(nlayer,nw) :: gcld                             ! cloud asymmetry parameter
    57       real, dimension(nlayer,nw,nabs) :: dtgas                       ! optical depth for each gas
    58       real, dimension(nlayer,nw) :: dagas                            ! total gas optical depth
    59       real, dimension(nlayer) :: edir, edn, eup                      ! normalised irradiances
    60       real, dimension(nlayer) :: fdir, fdn, fup                      ! normalised actinic fluxes
    61       real, dimension(nlayer) :: saflux                              ! total actinic flux
    62 
    63       integer, dimension(0:nlayer) :: nid
    64       real, dimension(0:nlayer,nlayer) :: dsdh
     49      real, dimension(nlayer+1) :: zpress, zalt, ztemp, zmmean       ! pressure (hpa)/altitude (km)/temperature (k)/mean molecular mass (g.mol-1)
     50
     51      real, dimension(nlayer+1) :: colinc                            ! air column increment (molecule.cm-2)
     52      real, dimension(nlayer+1,nw) :: dtrl                           ! rayleigh optical depth
     53      real, dimension(nlayer+1,nw) :: dtaer                          ! aerosol optical depth
     54      real, dimension(nlayer+1,nw) :: omaer                          ! aerosol single scattering albedo
     55      real, dimension(nlayer+1,nw) :: gaer                           ! aerosol asymmetry parameter
     56      real, dimension(nlayer+1,nw) :: dtcld                          ! cloud optical depth
     57      real, dimension(nlayer+1,nw) :: omcld                          ! cloud single scattering albedo
     58      real, dimension(nlayer+1,nw) :: gcld                           ! cloud asymmetry parameter
     59      real, dimension(nlayer+1,nw,nabs) :: dtgas                     ! optical depth for each gas
     60      real, dimension(nlayer+1,nw) :: dagas                          ! total gas optical depth
     61      real, dimension(nlayer+1) :: edir, edn, eup                    ! normalised irradiances
     62      real, dimension(nlayer+1) :: fdir, fdn, fup                    ! normalised actinic fluxes
     63      real, dimension(nlayer+1) :: saflux                            ! total actinic flux
     64
     65      integer, dimension(0:nlayer+1) :: nid
     66      real, dimension(0:nlayer+1,nlayer+1) :: dsdh
    6567
    6668      integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d, j_o3_o1d, j_o3_o,
     
    7274     $           a_hocl, a_so2, a_so, a_so3, a_s2, a_clo, a_ocs,
    7375     $           a_cocl2, a_h2so4, a_no2, a_no, a_n2, a_h2
    74       integer :: i, ilay, iw, ialt
     76      integer :: nlev, i, ilay, ilev, iw, ialt
    7577      real    :: deltaj
    7678      logical :: deutchem
     
    131133!     j_hdo_d   =        ! hdo + hv    -> d + oh
    132134
     135!==== define working vertical grid for the uv radiative code
     136
     137      nlev = nlayer + 1
     138
     139      do ilev = 1,nlev-1
     140         zpress(ilev) = press(ilev)
     141         zalt(ilev)   = alt(ilev)
     142         ztemp(ilev)  = temp(ilev)
     143         zmmean(ilev) = mmean(ilev)
     144      end do
     145
     146      zpress(nlev) = 0.                  ! top of the atmosphere
     147      zalt(nlev)   = zalt(nlev-1) + (zalt(nlev-1) - zalt(nlev-2))
     148      ztemp(nlev)  = ztemp(nlev-1)
     149      zmmean(nlev) = zmmean(nlev-1)
     150
    133151!==== air column increments and rayleigh optical depth
    134152
    135       call setair(nlayer, nw, wl, wc, press, temp, zmmean, colinc, dtrl)
     153      call setair(nlev, nw, wl, wc, zpress, ztemp, zmmean, colinc, dtrl)
    136154
    137155!==== set temperature-dependent cross-sections and optical depths
     156
     157      dtgas(:,:,:) = 0.
    138158
    139159!     o2:
     
    141161      call seto2(nphot, nlayer, nw, wc, mopt, temp, xso2_150, xso2_200,
    142162     $           xso2_250, xso2_300, yieldo2, j_o2_o, j_o2_o1d,
    143      $           colinc, rm(:,i_o2), dtgas(:,:,a_o2), sj)
     163     $           colinc(1:nlayer), rm(:,i_o2),
     164     $           dtgas(1:nlayer,:,a_o2), sj)
    144165
    145166!     co2:
     
    147168      call setco2(nphot, nlayer, nw, wc, temp, xsco2_195, xsco2_295,
    148169     $            xsco2_370, yieldco2, j_co2_o, j_co2_o1d,
    149      $            colinc, rm(:,i_co2), dtgas(:,:,a_co2), sj)
     170     $            colinc(1:nlayer), rm(:,i_co2),
     171     $            dtgas(1:nlayer,:,a_co2), sj)
    150172     
    151173!     o3:
    152174
    153175      call seto3(nphot, nlayer, nw, wc, temp, xso3_218, xso3_298,
    154      $           j_o3_o, j_o3_o1d,
    155      $           colinc, rm(:,i_o3), dtgas(:,:,a_o3), sj)
     176     $           j_o3_o, j_o3_o1d, colinc(1:nlayer), rm(:,i_o3),
     177     $           dtgas(1:nlayer,:,a_o3), sj)
    156178
    157179!     h2o2:
    158180
    159181      call seth2o2(nphot, nlayer, nw, wc, temp, xsh2o2, j_h2o2,
    160      $             colinc, rm(:,i_h2o2), dtgas(:,:,a_h2o2), sj)
     182     $             colinc(1:nlayer), rm(:,i_h2o2),
     183     $             dtgas(1:nlayer,:,a_h2o2), sj)
    161184
    162185!     so2:
    163186
    164187      call setso2(nphot, nlayer, nw, wc, temp, xsso2_200,xsso2_298,
    165      $            xsso2_360, j_so2,colinc, rm(:,i_so2),
    166      $            dtgas(:,:,a_so2), sj)
     188     $            xsso2_360, j_so2,colinc(1:nlayer), rm(:,i_so2),
     189     $            dtgas(1:nlayer,:,a_so2), sj)
    167190
    168191!     no2:
     
    170193      call setno2(nphot, nlayer, nw, wc, temp, xsno2, xsno2_220,
    171194     $            xsno2_294, yldno2_248, yldno2_298, j_no2,
    172      $            colinc, rm(:,i_no2), dtgas(:,:,a_no2), sj)
     195     $            colinc(1:nlayer), rm(:,i_no2),
     196     $            dtgas(1:nlayer,:,a_no2), sj)
    173197
    174198!==== temperature independent optical depths and cross-sections
     
    236260
    237261!      if (deutchem) then
    238 !         do ilay=1,nlayer
    239 !            do iw=1,nw-1
    240 !               !Two chanels for HDO, with same cross section
     262!         do ilay = 1,nlayer
     263!            do iw = 1,nw-1
     264!               !Two chanels for hdo, with same cross section
    241265!               sj(ilay,iw,j_hdo_od) = 0.5*xshdo(iw) ! hdo -> od + h
    242266!               sj(ilay,iw,j_hdo_d)  = 0.5*xshdo(iw) ! hdo -> d + oh
    243 !            enddo
    244 !         enddo
    245 !      endif
     267!            end do
     268!         end do
     269!      end if
    246270
    247271!==== set aerosol properties and optical depth
     
    249273      tau = 0. ! no solid aerosols for the present time
    250274
    251       call setaer(nlayer,alt,tau,nw,dtaer,omaer,gaer)
     275      call setaer(nlev,zalt,tau,nw,dtaer,omaer,gaer)
    252276
    253277!==== set cloud properties and optical depth
    254278
    255       call setcld(nlayer,alt,nw,wl,dtcld,omcld,gcld)
     279      call setcld(nlev,zalt,nw,wl,dtcld,omcld,gcld)
    256280
    257281!==== slant path lengths in spherical geometry
    258282
    259       call sphers(nlayer,alt,sza,dsdh,nid)
     283      call sphers(nlev,zalt,sza,dsdh,nid)
    260284
    261285!==== solar flux at venus
     
    280304!     dir = direct beam, dn = down-welling diffuse, up = up-welling diffuse
    281305
    282          call rtlink(nlayer, nw, iw, albedo(iw), sza, dsdh, nid, dtrl,
     306         call rtlink(nlev, nw, iw, albedo(iw), sza, dsdh, nid, dtrl,
    283307     $               dagas, dtcld, omcld, gcld, dtaer, omaer, gaer,
    284308     $               edir, edn, eup, fdir, fdn, fup)
     
    348372         colinc(ilev) = avo*0.1*dp/(zmmean(ilev)*g)
    349373      end do
    350       colinc(nlev) = avo*0.1*press(nlev)*100./(zmmean(nlev)*g)
     374      colinc(nlev) = 0.
    351375
    352376      do iw = 1, nw - 1
     
    834858!==============================================================================
    835859
    836       subroutine setaer(nlayer,alt,tau,nw,dtaer,omaer,gaer)
     860      subroutine setaer(nlev,zalt,tau,nw,dtaer,omaer,gaer)
    837861
    838862!-----------------------------------------------------------------------------*
     
    846870!     input
    847871
    848       integer :: nlayer                              ! number of vertical layers
    849       integer :: nw                                  ! number of wavelength grid points
    850       real, dimension(nlayer) :: alt                 ! altitude (km)
    851       real :: tau                                    ! integrated aerosol optical depth at the surface
     872      integer :: nlev                              ! number of vertical layers
     873      integer :: nw                                ! number of wavelength grid points
     874      real, dimension(nlev) :: zalt                ! altitude (km)
     875      real :: tau                                  ! integrated aerosol optical depth at the surface
    852876
    853877!     output
    854878
    855       real, dimension(nlayer,nw) :: dtaer            ! aerosol optical depth
    856       real, dimension(nlayer,nw) :: omaer            ! aerosol single scattering albedo
    857       real, dimension(nlayer,nw) :: gaer             ! aerosol asymmetry parameter
     879      real, dimension(nlev,nw) :: dtaer            ! aerosol optical depth
     880      real, dimension(nlev,nw) :: omaer            ! aerosol single scattering albedo
     881      real, dimension(nlev,nw) :: gaer             ! aerosol asymmetry parameter
    858882
    859883!     local
    860884
    861       integer :: ilay, iw
    862       real, dimension(nlayer) :: aer                                 ! dust extinction
     885      integer :: ilev, iw
     886      real, dimension(nlev) :: aer                 ! dust extinction
    863887      real :: omega, g, scaleh, gamma
    864888      real :: dz, tautot, q0
     
    879903
    880904       tautot = 0.
    881        do ilay = 1, nlayer-1
    882           dz = alt(ilay+1) - alt(ilay)
    883           tautot = tautot + exp(gamma*(1. - exp(alt(ilay)/scaleh)))*dz
     905       do ilev = 1, nlev-1
     906          dz = zalt(ilev+1) - zalt(ilev)
     907          tautot = tautot + exp(gamma*(1. - exp(zalt(ilev)/scaleh)))*dz
    884908       end do
    885909     
    886910       q0 = tau/tautot
    887        do ilay = 1, nlayer-1
    888           dz = alt(ilay+1) - alt(ilay)
    889           dtaer(ilay,:) = q0*exp(gamma*(1. - exp(alt(ilay)/scaleh)))*dz
    890           omaer(ilay,:) = omega
    891           gaer(ilay,:)  = g
     911       do ilev = 1, nlev-1
     912          dz = zalt(ilev+1) - zalt(ilev)
     913          dtaer(ilev,:) = q0*exp(gamma*(1. - exp(zalt(ilev)/scaleh)))*dz
     914          omaer(ilev,:) = omega
     915          gaer(ilev,:)  = g
    892916       end do
    893917
Note: See TracChangeset for help on using the changeset viewer.