Changeset 3530 for trunk


Ignore:
Timestamp:
Dec 2, 2024, 5:02:15 PM (39 hours ago)
Author:
flefevre
Message:

Implementation of temperature-dependent SO absorption cross-sections from Heays et al., 2023.

The SO cross-section file must be downloaded from the link below and put in the /cross_sections directory:

https://owncloud.latmos.ipsl.fr/index.php/s/3xLzwn7G587sTdB

Location:
trunk/LMDZ.VENUS/libf/phyvenus
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/LMDZ.VENUS/libf/phyvenus/photolysis_mod.F90

    r2925 r3530  
    4040  real, dimension(nw), save :: xshocl                                 ! hocl absorption cross-section (cm2)
    4141  real, dimension(nw), save :: xsso2_200, xsso2_298, xsso2_360        ! so2 absorption cross-section (cm2)
    42   real, dimension(nw), save :: xsso                                   ! so absorption cross-section (cm2)
     42  real, dimension(nw), save :: xsso_150, xsso_250                     ! so absorption cross-section (cm2)
    4343  real, dimension(nw), save :: xsso3                                  ! so3 absorption cross-section (cm2)
    4444  real, dimension(nw), save :: xsclo                                  ! clo absorption cross-section (cm2)
     
    5858  real, dimension(nw), save :: albedo                                 ! surface albedo
    5959
    60 !$OMP THREADPRIVATE(f,xsco2_195, xsco2_295, xsco2_370,yieldco2,xso2_150, xso2_200, xso2_250, xso2_300,yieldo2,xso3_218, xso3_298,xsh2o,xsh2o2,xsho2,xsh2,yieldh2,xsno2, xsno2_220, xsno2_294,yldno2_248,yldno2_298)
    61 !$OMP THREADPRIVATE(xsno,yieldno,xsn2,yieldn2,xshdo,albedo)
     60!$OMP THREADPRIVATE(f, xsco2_195, xsco2_295, xsco2_370, yieldco2, xso2_150, xso2_200, xso2_250, xso2_300, yieldo2, xso3_218, xso3_298, xsh2o, xsh2o2, xsho2, xsh2, yieldh2, xsno2, xsno2_220, xsno2_294, yldno2_248, yldno2_298)
     61!$OMP THREADPRIVATE(xsno, yieldno, xsn2, yieldn2, xshdo, albedo)
    6262
    6363contains
     
    123123! read and grid so cross-sections
    124124
    125   call rdxsso(nw,wl,xsso)
     125  call rdxsso(nw,wl,xsso_150,xsso_250)
    126126
    127127! read and grid so3 cross-sections
     
    23272327!==============================================================================
    23282328
    2329       subroutine rdxsso(nw, wl, yg)
     2329      subroutine rdxsso(nw,wl,xsso_150,xsso_250)
    23302330
    23312331!-----------------------------------------------------------------------------*
    23322332!=  PURPOSE:                                                                 =*
    2333 !=  Read SO cross-sections                                                  =*
    2334 !=  JPL 2006 recommendation                                                  =*
     2333!=  Read SO cross-sections                                                   =*
     2334!=  Phillips (1981) or Heays et al. (2023)                                   =*
    23352335!-----------------------------------------------------------------------------*
    23362336!=  PARAMETERS:                                                              =*
     
    23422342      USE mod_phys_lmdz_transfert_para, ONLY: bcast
    23432343
    2344       IMPLICIT NONE
     2344      implicit none
    23452345
    23462346!     input
     2347
    23472348      integer :: nw               ! number of wavelength grid points
    2348       real, dimension(nw) :: wl  ! lower wavelength for each interval
     2349      real, dimension(nw) :: wl   ! lower wavelength for each interval
    23492350
    23502351!     output
    23512352
    2352       real, dimension(nw) :: yg   ! SO cross-sections (cm2)
     2353      real, dimension(nw) :: xsso_150, xsso_250 ! so cross-sections (cm2)
    23532354
    23542355!     local
     2356
     2357      integer, parameter :: kdata = 1000
     2358      integer :: i, l, n, n1, n2, ierr, iopt
     2359      integer :: kin, kout ! input/output logical units
     2360      real, dimension(nw) :: yg, yg1, yg2
     2361      real, dimension(kdata) :: x1, y1, x2, y2
    23552362      real :: dummy
    23562363      real, parameter :: deltax = 1.e-4
    2357       integer, parameter :: kdata = 1000
    2358       real, dimension(kdata) :: x1, y1
    2359       integer :: i, n, ierr
    23602364      character*100 fil
    2361       integer :: kin, kout ! input/output logical units
    23622365
    23632366      kin = 10
    23642367
    2365 !*** cross sections from JPL [2006]
    2366 
    2367       fil = 'cross_sections/so_marcq.txt'
    2368       print*, 'section efficace SO: ', fil
    2369 
    2370       if(is_master) then
     2368!     iopt = 1 : Phillips (1981)                       
     2369!     iopt = 2 : Heays et al., Molecular Physics, 2023
     2370
     2371      iopt = 2
     2372
     2373      if (iopt == 1) then
     2374
     2375!        cross sections from Phillips (1981)
     2376
     2377         fil = 'cross_sections/so_marcq.txt'
     2378
     2379         if (is_master) then
     2380
     2381         print*, 'section efficace SO: ', fil
    23712382
    23722383         n = 851
    2373       OPEN(kin,FILE=fil,STATUS='OLD')
    2374       DO i = 1,n
    2375         READ(kin,*) x1(i), dummy, dummy, dummy, dummy, y1(i)
    2376         y1(i) = y1(i) * 0.88
    2377       ENDDO
    2378       CLOSE(kin)
    2379 
    2380       CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
    2381       CALL addpnt(x1,y1,kdata,n,          0.,0.)
    2382       CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
    2383       CALL addpnt(x1,y1,kdata,n,        1E38,0.)
    2384 
    2385       CALL inter2(nw,wl,yg,n,x1,y1,ierr)
     2384         OPEN(kin,FILE=fil,STATUS='OLD')
     2385         DO i = 1,n
     2386            READ(kin,*) x1(i), dummy, dummy, dummy, dummy, y1(i)
     2387            y1(i) = y1(i)*0.88  ! from Mills, PhD, 1998
     2388         ENDDO
     2389         CLOSE(kin)
     2390
     2391         CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
     2392         CALL addpnt(x1,y1,kdata,n,          0.,0.)
     2393         CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
     2394         CALL addpnt(x1,y1,kdata,n,        1E38,0.)
     2395
     2396         do i = 1,n
     2397            print*, x1(i), y1(i)
     2398         end do
     2399
     2400         CALL inter2(nw,wl,yg,n,x1,y1,ierr)
    23862401 
    2387       IF (ierr .NE. 0) THEN
    2388         WRITE(*,*) ierr, fil
    2389         STOP
    2390       ENDIF
    2391 
    2392       endif !is_master
     2402         IF (ierr .NE. 0) THEN
     2403           WRITE(*,*) ierr, fil
     2404           STOP
     2405         ENDIF
     2406
     2407         DO l = 1, nw-1
     2408            xsso_150(l) = yg(l)
     2409            xsso_250(l) = yg(l)
     2410         END DO
     2411
     2412         endif !is_master
     2413
     2414      else if (iopt == 2) then
     2415
     2416!        cross sections from Heays et al. (2023)
     2417!        values given at 150 K and 250 K
     2418
     2419         fil = 'cross_sections/SO_heays_2023.txt'
     2420
     2421         if (is_master) then
     2422
     2423         print*, 'section efficace SO: ', fil
     2424         
     2425         n1 = 525
     2426         n2 = n1
     2427         OPEN(kin,FILE=fil,STATUS='OLD')
     2428
     2429         DO i = 1, 3
     2430            READ(kin,*)
     2431         ENDDO
     2432
     2433         DO i = 1,n1
     2434            READ(kin,*) x1(i), dummy, dummy, dummy, dummy, dummy, &
     2435                        dummy, dummy, dummy, y1(i), y2(i)
     2436            x2(i) = x1(i)
     2437         ENDDO
     2438         CLOSE(kin)
     2439
     2440         CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
     2441         CALL addpnt(x1,y1,kdata,n1,               0.,0.)
     2442         CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
     2443         CALL addpnt(x1,y1,kdata,n1,           1.e+38,0.)
     2444
     2445         CALL inter2(nw,wl,yg1,n1,x1,y1,ierr)
     2446
     2447         IF (ierr .NE. 0) THEN
     2448            WRITE(*,*) ierr, fil
     2449            STOP
     2450         ENDIF
     2451
     2452         CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.)
     2453         CALL addpnt(x2,y2,kdata,n2,               0.,0.)
     2454         CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
     2455         CALL addpnt(x2,y2,kdata,n2,           1.e+38,0.)
     2456
     2457         CALL inter2(nw,wl,yg2,n2,x2,y2,ierr)
     2458
     2459         IF (ierr .NE. 0) THEN
     2460            WRITE(*,*) ierr, fil
     2461            STOP
     2462         ENDIF
     2463
     2464         DO l = 1, nw-1
     2465            xsso_150(l) = yg1(l)
     2466            xsso_250(l) = yg2(l)
     2467         END DO
     2468
     2469         endif !is_master
     2470
     2471      end if ! iopt
    23932472
    23942473      call bcast(yg)
  • trunk/LMDZ.VENUS/libf/phyvenus/photolysis_online.F

    r2974 r3530  
    183183     $             dtgas(1:nlayer,:,a_h2o2), sj)
    184184
     185!     so:
     186
     187      call setso(nphot, nlayer, nw, wc, temp, xsso_150, xsso_250,
     188     $           j_so, colinc(1:nlayer), rm(:,i_so),
     189     $           dtgas(1:nlayer,:,a_so), sj)
     190
    185191!     so2:
    186192
    187       call setso2(nphot, nlayer, nw, wc, temp, xsso2_200,xsso2_298,
    188      $            xsso2_360, j_so2,colinc(1:nlayer), rm(:,i_so2),
     193      call setso2(nphot, nlayer, nw, wc, temp, xsso2_200, xsso2_298,
     194     $            xsso2_360, j_so2, colinc(1:nlayer), rm(:,i_so2),
    189195     $            dtgas(1:nlayer,:,a_so2), sj)
    190196
     
    209215            dtgas(ilay,iw,a_hocl) = colinc(ilay)*rm(ilay,i_hocl)
    210216     $      *xshocl(iw)
    211             dtgas(ilay,iw,a_so) = colinc(ilay)*rm(ilay,i_so)*xsso(iw)
    212217            dtgas(ilay,iw,a_so3) = colinc(ilay)*rm(ilay,i_so3)*xsso3(iw)
    213218            dtgas(ilay,iw,a_s2) = colinc(ilay)*rm(ilay,i_s2)*xss2(iw)
     
    245250            sj(ilay,iw,j_cl2) = xscl2(iw)                ! cl2
    246251            sj(ilay,iw,j_hocl) = xshocl(iw)              ! hocl
    247             sj(ilay,iw,j_so) = xsso(iw)                  ! so
    248252            sj(ilay,iw,j_s2) = xss2(iw)                  ! s2
    249253            sj(ilay,iw,j_so3) = xsso3(iw)                ! so3
     
    647651
    648652      end subroutine seto3
     653!==============================================================================
     654
     655      subroutine setso(nd, nlayer, nw, wc, tlay, xsso_150, xsso_250,
     656     $                 j_so, colinc, rm, dt, sj)
     657
     658!-----------------------------------------------------------------------------*
     659!=  PURPOSE:                                                                 =*
     660!=  Set up the so temperature dependent cross-sections and optical depth     =*
     661!-----------------------------------------------------------------------------*
     662
     663      implicit none
     664
     665!     input:
     666
     667      integer :: nd                                               ! number of photolysis rates
     668      integer :: nlayer                                           ! number of vertical layers
     669      integer :: nw                                               ! number of wavelength grid points
     670      integer :: j_so                                             ! photolysis indexe
     671      real, dimension(nw)     :: wc                               ! central wavelength for each interval
     672      real, dimension(nw)     :: xsso_150, xsso_250               ! so cross-sections (cm2)
     673      real, dimension(nlayer) :: tlay                             ! temperature (k)
     674      real, dimension(nlayer) :: rm                               ! so mixing ratio
     675      real, dimension(nlayer) :: colinc                           ! air column increment (molecule.cm-2)
     676
     677!     output:
     678
     679      real, dimension(nlayer,nw)    :: dt                         !  optical depth
     680      real, dimension(nlayer,nw,nd) :: sj                         !  cross-section array (cm2)
     681
     682!     local:
     683!
     684      integer :: ilev, iw
     685      real :: temp, xsso
     686
     687!     calculate temperature dependance
     688
     689      do ilev = 1,nlayer
     690         temp = max(tlay(ilev),150.)
     691         temp = min(temp, 250.)
     692         do iw = 1, nw-1
     693            if (xsso_150(iw) /= 0. .and. xsso_250(iw) /= 0.) then
     694               xsso = log(xsso_150(iw))
     695     $              + (log(xsso_250(iw)) - log(xsso_150(iw)))
     696     $               /(250. - 150.)*(temp - 150.)
     697
     698               sj(ilev,iw,j_so) = exp(xsso)
     699            else
     700               sj(ilev,iw,j_so) = 0.
     701            end if
     702
     703!     optical depth
     704
     705            dt(ilev,iw) = colinc(ilev)*rm(ilev)*sj(ilev,iw,j_so)
     706
     707         end do
     708      end do
     709
     710      end subroutine setso
     711
    649712!==============================================================================
    650713
Note: See TracChangeset for help on using the changeset viewer.