Changeset 3794


Ignore:
Timestamp:
Jun 5, 2025, 3:39:30 PM (2 days ago)
Author:
slebonnois
Message:

SL: VenusPCM - update of CIA in solar module (adapted from Generic PCM) + potential output of band-resolved net flux

Location:
trunk/LMDZ.VENUS/libf/phyvenus
Files:
1 added
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.VENUS/libf/phyvenus/optcv.F90

    r2560 r3794  
    1111  use radinc_h, only: L_NLAYRAD, L_NLEVRAD, L_LEVELS, L_NSPECTV, L_NGAUSS, L_REFVAR, NAERKIND
    1212  use radcommon_h, only: gasv, tlimit, tgasref, pfgasref,wnov,scalep,indv
    13   use gases_h, only: gfrac, ngasmx, igas_H2, igas_H2O, igas_He, igas_N2
     13  use gases_h, only: gfrac, ngasmx, igas_H2O, igas_CO2
     14  use interpolate_continuum_mod, only: interpolate_continuum
    1415
    1516  implicit none
    … …  
    7273  real*8 DRAYAER
    7374  double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
    74   double precision p_cross
    7575
    7676  ! variable species mixing ratio variables
    … …  
    166166!        if(continuum)then
    167167           ! include continua if necessary
    168            wn_cont = dble(wnov(nw))
    169168           T_cont  = dble(TMID(k))
    170169           do igas=1,ngasmx
    … …  
    178177              dtemp=0.0
    179178
    180 ! For Venus: only H2O, using CKD
    181               if(igas.eq.igas_H2O.and.T_cont.gt.200.0)then
    182 
    183                  p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
    184 !                 if(H2Ocont_simple)then
    185 !                    call interpolateH2Ocont_PPC(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
    186 !                 else
    187                     interm = indv(nw,igas,igas)
    188                     call interpolateH2Ocont_CKD(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm)
    189                     indv(nw,igas,igas) = interm
    190 !                 endif
    191 
    192               endif
     179! For Venus: only H2O and CO2
     180              do jgas=1,ngasmx
     181                if ( ((igas .eq. igas_H2O) .and. (jgas .eq. igas_H2O)) .or.  &
     182                     ((igas .eq. igas_H2O) .and. (jgas .eq. igas_CO2)) .or.   &
     183                     ((igas .eq. igas_CO2) .and. (jgas .eq. igas_CO2)) ) then
     184                  call interpolate_continuum('',igas,jgas,'VI',nw,T_cont,p_cont,p_cont,dtemp,.false.)
     185                endif
     186              enddo
    193187
    194188              DCONT = DCONT + dtemp

  • trunk/LMDZ.VENUS/libf/phyvenus/radcommon_h.F90

    r3775 r3794  
    11module radcommon_h
    2       use radinc_h, only: L_NSPECTV, NTstart, NTstop, &
     2      use radinc_h, only: L_NSPECTI, L_NSPECTV, NTstart, NTstop, &
    33                          naerkind, nsizemax
    44      implicit none
    … …  
    1313!  some or all of this common data set
    1414!
     15!     WNOI       - Array of wavenumbers at the spectral interval
     16!                  centers for the infrared.  Array is NSPECTI
     17!                  elements long.
     18!     DWNI       - Array of "delta wavenumber", i.e., the width,
     19!                  in wavenumbers (cm^-1) of each IR spectral
     20!                  interval.  NSPECTI elements long.
     21!     WAVEI      - Array (NSPECTI elements long) of the wavelenght
     22!                  (in microns) at the center of each IR spectral
     23!                  interval.
    1524!     WNOV       - Array of wavenumbers at the spectral interval
    1625!                  center for the VISUAL.  Array is NSPECTV
    … …  
    2938!     TAURAYVAR  - Array (NSPECTV elements) of the pressure-independent
    3039!                  part of Rayleigh scattering optical depth for the variable gas.
     40!     FZEROI     - Fraction of zeros in the IR CO2 k-coefficients, for
     41!                  each temperature, pressure, and spectral interval
    3142!     FZEROV     - Fraction of zeros in the VISUAL CO2 k-coefficients, for
    3243!                  each temperature, pressure, and spectral interval
    … …  
    4253! gvis       :  assymetry factor
    4354!
     55!   Longwave
     56!   ~~~~~~~~
     57!
     58! For the "naerkind" kind of aerosol radiative properties :
     59! QIRsQREF :  Qext / Qext("longrefir")
     60! omegaIR  :  mean single scattering albedo
     61! gIR      :  mean assymetry factor
     62!
     63!
     64! Note - QIRsQREF in the martian model was scaled to longrefvis,
     65! here it is scaled to longrefir, which is actually a dummy parameter,
     66! as we do not output scaled aerosol opacity ...
     67!
    4468
     69      REAL*8 BWNI(L_NSPECTI+1), WNOI(L_NSPECTI), DWNI(L_NSPECTI), WAVEI(L_NSPECTI) !BWNI read by master in setspi
    4570      REAL*8 BWNV(L_NSPECTV+1), WNOV(L_NSPECTV), DWNV(L_NSPECTV), WAVEV(L_NSPECTV) !BWNV read by master in setspv
    4671      REAL*8 STELLARF(L_NSPECTV), TAURAY(L_NSPECTV), TAURAYVAR(L_NSPECTV)
    47 !$OMP THREADPRIVATE(WNOV,DWNV,WAVEV,&
     72!$OMP THREADPRIVATE(WNOI,DWNI,WAVEI,&
     73        !$OMP WNOV,DWNV,WAVEV,&
    4874        !$OMP STELLARF,TAURAY,TAURAYVAR)
    4975
     76      REAL*8 blami(L_NSPECTI+1)
    5077      REAL*8 blamv(L_NSPECTV+1) ! these are needed by suaer.F90
    51 !$OMP THREADPRIVATE(blamv)
     78!$OMP THREADPRIVATE(blami,blamv)
    5279
    5380      !! AS: introduced to avoid doing same computations again for continuum
     81      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: indi
    5482      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: indv
    55 !$OMP THREADPRIVATE(indv)
     83!$OMP THREADPRIVATE(indi,indv)
    5684
    5785      !!! ALLOCATABLE STUFF SO THAT DIMENSIONS ARE READ in *.dat FILES -- AS 12/2011 
    58       REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasv
     86      REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasi, gasv
    5987      REAL*8, DIMENSION(:), ALLOCATABLE :: PGASREF, TGASREF, PFGASREF, GWEIGHT
     88      real*8 FZEROI(L_NSPECTI)
    6089      real*8 FZEROV(L_NSPECTV)
    6190      real*8 pgasmin, pgasmax
    6291      real*8 tgasmin, tgasmax
    63 !$OMP THREADPRIVATE(gasv,&  !pgasref,tgasref,pfgasref read by master in sugas_corrk
    64         !$OMP FZEROV)           !pgasmin,pgasmax,tgasmin,tgasmax read by master in sugas_corrk
     92!$OMP THREADPRIVATE(gasi,gasv,&  !pgasref,tgasref,pfgasref read by master in sugas_corrk
     93        !$OMP FZEROI,FZEROV)           !pgasmin,pgasmax,tgasmin,tgasmax read by master in sugas_corrk
    6594
    6695      real QVISsQREF(L_NSPECTV,naerkind,nsizemax)
    6796      real omegavis(L_NSPECTV,naerkind,nsizemax)
    6897      real gvis(L_NSPECTV,naerkind,nsizemax)
    69 !$OMP THREADPRIVATE(QVISsQREF,omegavis,gvis)
     98      real QIRsQREF(L_NSPECTV,naerkind,nsizemax)
     99      real omegair(L_NSPECTV,naerkind,nsizemax)
     100      real gir(L_NSPECTV,naerkind,nsizemax)
     101!$OMP THREADPRIVATE(QVISsQREF,omegavis,gvis,QIRsQREF,omegair,gir)
    70102
    71103
    … …  
    73105! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    74106
    75       REAL lamrefvis(naerkind)
     107      REAL lamrefir(naerkind),lamrefvis(naerkind)
    76108
    77109! Actual number of grain size classes in each domain for a
    … …  
    87119
    88120! Extinction coefficient at reference wavelengths;
    89 !   These wavelengths are defined in aeroptproperties, and called longrefvis
     121!   These wavelengths are defined in aeroptproperties, and called longrefvis and longrefir
    90122
    91123      REAL :: QREFvis(naerkind,nsizemax)
     124      REAL :: QREFir(naerkind,nsizemax)
     125      REAL :: omegaREFir(naerkind,nsizemax)
    92126
    93127      REAL,SAVE :: tstellar ! Stellar brightness temperature (SW)
    94 
     128      REAL*8,SAVE :: planckir(naerkind,nsizemax)
    95129      real*8,save :: PTOP
    96130
    97131      real*8,parameter :: UBARI = 0.5D0
    98132
    99 !$OMP THREADPRIVATE(QREFvis,&   ! gweight read by master in sugas_corrk
    100                 !$OMP tstellar,PTOP)
     133!$OMP THREADPRIVATE(QREFvis,QREFir,omegaREFir,&         ! gweight read by master in sugas_corrk
     134                !$OMP tstellar,planckir,PTOP)
    101135
    102136!     If the gas optical depth (top to the surface) is less than

  • trunk/LMDZ.VENUS/libf/phyvenus/radinc_h.F90

    r2560 r3794  
    7979
    8080      integer, parameter :: L_NSPECTV = NBvisible
     81      integer, parameter :: L_NSPECTI = NBinfrared
    8182
    8283      integer, parameter :: NAERKIND  = 5

  • trunk/LMDZ.VENUS/libf/phyvenus/sfluxv.F

    r2560 r3794  
    11      SUBROUTINE SFLUXV(DTAUV,TAUV,TAUCUMV,RSFV,DWNV,WBARV,COSBV,
    22     *                  UBAR0,STEL,NFLUXTOPV,FLUXTOPVDN,
    3      *                  NFLUXOUTV_nu,NFLUXGNDV_nu,
     3     *                  NFLUXV_nu,NFLUXGNDV_nu,
    44     *                  FMNETV,FLUXUPV,FLUXDNV,FZEROV,taugsurf)
    55
    … …  
    1919      real*8 FLUXUPV(L_NLAYRAD), FLUXDNV(L_NLAYRAD)
    2020      real*8 NFLUXTOPV, FLUXUP, FLUXDN,FLUXTOPVDN
    21       real*8 NFLUXOUTV_nu(L_NSPECTV)
     21      real*8 NFLUXV_nu(L_NLAYRAD,L_NSPECTV)
    2222      real*8 NFLUXGNDV_nu(L_NSPECTV)
    2323
    … …  
    4040
    4141      DO NW=1,L_NSPECTV
    42          NFLUXOUTV_nu(NW)=0.0
     42         NFLUXV_nu(:,NW)=0.0
    4343         NFLUXGNDV_nu(NW)=0.0
    4444      END DO
    … …  
    115115          END DO
    116116
    117 c     band-resolved flux leaving TOA (RDW)
    118           NFLUXOUTV_nu(NW) = NFLUXOUTV_nu(NW)
    119      *      +FLUXUP*GWEIGHT(NG)*(1.0-FZEROV(NW))
     117c     band-resolved net flux
     118          NFLUXV_nu(:,NW) = NFLUXV_nu(:,NW)
     119     *      +(FMUPV-FMDV)*GWEIGHT(NG)*(1.0-FZEROV(NW))
    120120
    121121c     band-resolved flux at ground (RDW)
    … …  
    173173        END DO
    174174
    175 c     band-resolved flux leaving TOA (RDW)
    176           NFLUXOUTV_nu(NW) = NFLUXOUTV_nu(NW)
    177      *      +FLUXUP*FZERO
     175c     band-resolved net flux
     176          NFLUXV_nu(:,NW) = NFLUXV_nu(:,NW)
     177     *      +(FMUPV-FMDV)*FZERO
    178178
    179179c     band-resolved flux at ground (RDW)

  • trunk/LMDZ.VENUS/libf/phyvenus/su_gases.F90

    r2560 r3794  
    3939   igas_OCS=5
    4040
     41  ! We don't care about the others
     42  ! but initialized in order to avoid problems in interpolate_continuum.F90
     43  igas_H2=99
     44  igas_He=99
     45  igas_N2=99
     46  igas_CH4=99
     47  igas_O2=99
     48
    4149  vgas=0
    4250  if(.not.allocated(gfrac)) allocate(gfrac(ngasmx))

  • trunk/LMDZ.VENUS/libf/phyvenus/sugas_corrk.F90

    r2560 r3794  
    2727      use radcommon_h, only : WNOV
    2828      use datafile_mod, only: datadir,corrkdir,banddir
    29       use gases_h, only: gnom, ngasmx, igas_H2O
     29      use gases_h, only: gnom, ngasmx, igas_H2O, igas_CO2
     30      use interpolate_continuum_mod, only: interpolate_continuum
    3031      implicit none
    3132#include "YOMCST.h"
    … …  
    124125                 'match that in gases.def (',trim(gnom(igas)),'). You should compare ', &
    125126                 'gases.def with Q.dat in your radiative transfer directory.'
     127            print*,'For Venus, no gases.def, but su_gases.F90 hardcoded... Should be identical ', &
     128                 'to Q.dat in your radiative transfer directory.'
    126129            message='exiting.'
    127130            call abort_physic(subname,message,1)
    … …  
    452455!=======================================================================
    453456!     Initialise the continuum absorption data
    454 !      if(continuum)then
    455       do igas=1,ngasmx
    456 
    457 ! For Venus: only H2O, using CKD
    458          if (igas .eq. igas_H2O) then
    459 
    460             ! H2O is special
    461 !            if(H2Ocont_simple)then
    462 !               call interpolateH2Ocont_PPC(990.D+0,296.D+0,683.2D+0*2,0.D+0,testcont,.true.)
    463 !            else
    464                dummy = -9999
    465                call interpolateH2Ocont_CKD(990.D+0,296.D+0,683.2D+0*2,0.D+0,testcont,.true.,dummy)
    466 !            endif
    467 
    468          endif
    469 
    470       enddo
    471 !      endif
     457!     if(continuum)then
     458       do igas=1,ngasmx ! we loop on all pairs of molecules that have data available
     459       ! data can be downloaded from https://web.lmd.jussieu.fr/~lmdz/planets/generic/datagcm/continuum_data/
     460
     461! For Venus: only H2O and CO2
     462        if (igas .eq. igas_H2O) then
     463         file_id='/continuum_data/' // 'H2O-H2O_2022.cia'
     464         ! H2O-H2O_2022.cia = H2O-H2O_continuum_100-2000K_2025.dat
     465         file_path=TRIM(datadir)//TRIM(file_id)
     466         call interpolate_continuum(file_path,igas_H2O,igas_H2O,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.)
     467         do jgas=1,ngasmx
     468          if (jgas .eq. igas_CO2) then
     469           file_id='/continuum_data/' // 'H2O-CO2_2022.cia'
     470           ! H2O-CO2_2022.cia = H2O-CO2_continuum_100-800K_2025.dat
     471           file_path=TRIM(datadir)//TRIM(file_id)
     472           call interpolate_continuum(file_path,igas_H2O,igas_CO2,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.)
     473          endif
     474         enddo   
     475        elseif (igas .eq. igas_CO2) then
     476         file_id='/continuum_data/' // 'CO2-CO2_2025.cia'
     477         ! CO2-CO2_2025.cia = CO2-CO2_continuum_100-800K_2025.dat
     478         file_path=TRIM(datadir)//TRIM(file_id)
     479         call interpolate_continuum(file_path,igas_CO2,igas_CO2,'IR',1,300.D+0,10000.D+0,20000.D+0,testcont,.true.)
     480        endif
     481       enddo ! igas=1,ngasmx
     482!     endif
    472483
    473484!      print*,'----------------------------------------------------'

  • trunk/LMDZ.VENUS/libf/phyvenus/sw_venus_corrk.F90

    r2560 r3794  
    99                             gweight, pfgasref, pgasmax, pgasmin, &
    1010                             pgasref, tgasmax, tgasmin, tgasref, scalep, &
    11                              stellarf, dwnv, tauray
     11                             stellarf, dwnv, tauray, wavev, wnov
    1212      use datafile_mod, only: datadir,banddir,corrkdir
    1313      use ioipsl_getin_p_mod, only: getin_p
    … …  
    9999      REAL*8 tauaero(L_LEVELS,naerkind)
    100100      REAL*8 nfluxtopv,fluxtopvdn
    101       REAL*8 nfluxoutv_nu(L_NSPECTV)                 ! Outgoing band-resolved VI flux at TOA (W/m2).
     101      REAL*8 nfluxv_nu(L_NLAYRAD,L_NSPECTV)  ! vertical profil, band-resolved VI net flux (W/m2)
     102      REAL*8 heatrate_nu(L_NLAYRAD,L_NSPECTV)  ! vertical profil, band-resolved VI heating rate (K/s/micron)
    102103      REAL*8 fmnetv(L_NLAYRAD)
    103104      REAL*8 fluxupv(L_NLAYRAD)
    … …  
    504505            call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv,  &
    505506                 acosz,stel_fract,                                 &
    506                  nfluxtopv,fluxtopvdn,nfluxoutv_nu,nfluxgndv_nu,   &
     507                 nfluxtopv,fluxtopvdn,nfluxv_nu,nfluxgndv_nu,      &
    507508                 fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
    508509
    … …  
    510511            nfluxtopv       = 0.0d0
    511512            fluxtopvdn      = 0.0d0
    512             nfluxoutv_nu(:) = 0.0d0
     513            nfluxv_nu(:,:) = 0.0d0
    513514            nfluxgndv_nu(:) = 0.0d0
    514515            do l=1,L_NLAYRAD
    … …  
    564565             *RG/(cpdet(tmid(1))*scalep*(plevrad(3)-plevrad(2)))
    565566
     567
    566568!-----------------------------------------------------------------------   
    567569!-----------------------------------------------------------------------   
    … …  
    570572!-----------------------------------------------------------------------
    571573
     574! Diagnostic in 1D:
     575!     output the vertical profil of band-resolved heating rates, in K/s/microns
     576      if ((1.eq.1).and.is_master.and.(ngrid.eq.1).and.firstcall) then
     577         open(unit=11,file="nfluxnu.dat")
     578         write(11,*) "lambda(microns)",wavev(:)
     579         write(11,*) "dlbda(microns)",1.e4*dwnv(:)/wnov(:)**2
     580         write(11,*) "pressure(Pa) dT_nu(K/s/micron)"
     581         DO l=2,L_NLAYRAD
     582            heatrate_nu(L_NLAYRAD+1-l,:)=(nfluxv_nu(l,:)-nfluxv_nu(l-1,:))  &
     583                *RG/(cpdet(tmid(l))*scalep*(plevrad(2*l+1)-plevrad(2*l-1)))  &
     584            ! dlbda(microns)=-dnu(cm-1)*1E4/nu(cm-1)^2
     585                / (1.e4*dwnv(:)/wnov(:)**2)
     586            write(11,*) pplay(1,L_NLAYRAD+1-l),heatrate_nu(L_NLAYRAD+1-l,:)
     587         END DO
     588         close(11)
     589      endif
     590         
    572591    end subroutine sw_venus_corrk
    573592
Note: See TracChangeset for help on using the changeset viewer.