Changeset 2804 for trunk/LMDZ.GENERIC


Ignore:
Timestamp:
Oct 19, 2022, 11:24:44 AM (2 years ago)
Author:
aslmd
Message:

This is the RGCS scheme

Now we can use this scheme to take into account radiative effect of a Generic Condensable Specie. One needs to add the keyword is_rgcs=1 to the ModernTraceur?.def to activate the scheme + specify the number of aerogeneric specie in callphys.def. Also, one needs to go into suaer_corrk.F90 to add the Mie scattering file into the code for the specie you want. Currently, MnS, Fe, Cr, Mg2SiO4 are available. This has been tested on a Hot Jupiter case, might have issue in different configuration. Also, the numerical stability of the scheme is not guaranted when using strong scatterers such as silicates. Physical oscillations in the evaporation/condensation can occur, and need fixing. This will be for another commit.

Location:
trunk/LMDZ.GENERIC/libf/phystd
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90

    r2803 r2804  
    3131!     dust removal, simplification by Robin Wordsworth (2009)
    3232!     Generic n-layer aerosol - J. Vatant d'Ollone (2020)
    33 !     Radiative Generic Condensable aerosols - Lucas Teinturier (2022)
     33!     Radiative Generic Condensable Species - Lucas Teinturier (2022)
    3434!
    3535!     Input
    … …  
    9494      integer igen_ice,igen_vap ! to store the index of generic tracer
    9595      logical dummy_bool ! dummy boolean just in case we need one
     96      ! integer i_rgcs_ice(aerogeneric)
    9697      ! identify tracers
    9798      IF (firstcall) THEN
    98 
     99        ia =0
    99100        write(*,*) "Tracers found in aeropacity:"
    100101        do iq=1,nq
    … …  
    662663          ! Initialization
    663664          aerosol(1:ngrid,1:nlayer,iaer) = 0.D0
    664           igen_ice = -1
    665           igen_vap = -1
    666           ! Get index of the ice tracer (for pq, mol, rho_q)
    667           call generic_tracer_index(nq,2*iaer,igen_vap,igen_ice,dummy_bool)
    668           if (igen_ice .eq. -1) then
    669             !if igen_ice not changed by generic_tracer_index, it means it's an ice tracer => igen_ice=2*iaer
    670             igen_ice = 2*iaer
    671           endif
     665          igen_ice = i_rgcs_ice(ia)
    672666          ! Let's loop on the horizontal and vertical grid
    673667          do ig=1,ngrid
    … …  
    694688      end do
    695689
    696       do ig=1,ngrid
    697          do l=1,nlayer
    698             do iaer = 1, naerkind
    699                if(aerosol(ig,l,iaer).gt.1.e3)then
    700                   print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
    701                   print*,'at ig=',ig,',  l=',l,', iaer=',iaer
    702                   print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
    703                   print*,'reffrad=',reffrad(ig,l,iaer)
    704                endif
    705             end do
    706          end do
    707       end do
    708 
    709       do ig=1,ngrid
    710          if(tau_col(ig).gt.1.e3)then
    711             print*,'WARNING: tau_col=',tau_col(ig)
    712             print*,'at ig=',ig
    713             print*,'aerosol=',aerosol(ig,:,:)
    714             print*,'QREFvis3d=',QREFvis3d(ig,:,:)
    715             print*,'reffrad=',reffrad(ig,:,:)
    716          endif
    717       end do
     690      ! do ig=1,ngrid
     691      !    do l=1,nlayer
     692      !       do iaer = 1, naerkind
     693      !          if(aerosol(ig,l,iaer).gt.1.e3)then
     694      !             print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
     695      !             print*,'at ig=',ig,',  l=',l,', iaer=',iaer
     696      !             print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
     697      !             print*,'reffrad=',reffrad(ig,l,iaer)
     698      !          endif
     699      !       end do
     700      !    end do
     701      ! end do
     702
     703      ! do ig=1,ngrid
     704      !    if(tau_col(ig).gt.1.e3)then
     705      !       print*,'WARNING: tau_col=',tau_col(ig)
     706      !       print*,'at ig=',ig
     707      !       print*,'aerosol=',aerosol(ig,:,:)
     708      !       print*,'QREFvis3d=',QREFvis3d(ig,:,:)
     709      !       print*,'reffrad=',reffrad(ig,:,:)
     710      !    endif
     711      ! end do
    718712      return
    719713    end subroutine aeropacity

  • trunk/LMDZ.GENERIC/libf/phystd/aerosol_mod.F90

    r2803 r2804  
    2424! Generic aerosols
    2525      integer, dimension(:),allocatable,save :: iaero_generic
    26 !$OMP THREADPRIVATE(iaero_co2,iaero_h2o,iaero_dust,iaero_h2so4,noaero,iaero_back2lay,iaero_nh3,iaero_nlay,iaero_aurora,iaero_generic)
     26      integer,dimension(:),allocatable,save :: i_rgcs_ice
     27!$OMP THREADPRIVATE(iaero_co2,iaero_h2o,iaero_dust,iaero_h2so4,noaero,iaero_back2lay,iaero_nh3,iaero_nlay,iaero_aurora,iaero_generic,i_rgcs_ice)
    2728     
    2829!==================================================================

  • trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90

    r2803 r2804  
    501501           reffrad,QREFvis3d,QREFir3d,                             &
    502502           tau_col,cloudfrac,totcloudfrac,clearsky)               
    503          
    504 
    505 
     503 
    506504!-----------------------------------------------------------------------   
    507505      do ig=1,ngrid ! Starting Big Loop over every GCM column

  • trunk/LMDZ.GENERIC/libf/phystd/generic_tracer_index_mod.F90

    r2722 r2804  
    2020
    2121        USE tracer_h, only: noms, is_condensable
    22         integer nq,iq, igcm_generic_vap, igcm_generic_ice
     22        integer nq,iq, igcm_generic_vap, igcm_generic_ice, i_shortname, ii
    2323        logical call_ice_vap_generic
    24 
     24        character(len=10) :: shortname
    2525        if((is_condensable(iq)==1) .and. (index(noms(iq),"vap") .ne. 0) &
    2626                                                    .and. (index(noms(iq),"h2o") .eq. 0) .and. (index(noms(iq),"co2") .eq. 0)) then
    … …  
    3030            igcm_generic_ice = -1
    3131
    32             ! We look for the corresponding ice traceur (before or after in the list of traceurs, maybe could be generalised to the whole list)
    33             if (iq .ne. nq) then
    34                 if ((noms(iq)(1:len(trim(noms(iq)))-4) .eq. noms(iq+1)(1:len(trim(noms(iq+1)))-4)) .and. (index(noms(iq+1),"ice") .ne. 0)) then
    35                     igcm_generic_ice = iq+1
    36                 end if
    37             end if
    38             if ((iq .gt. 1)) then
    39                 if ((noms(iq)(1:len(trim(noms(iq)))-4) .eq. noms(iq-1)(1:len(trim(noms(iq-1)))-4)) .and. (index(noms(iq-1),"ice") .ne. 0)) then
    40                         igcm_generic_ice = iq-1
    41                 end if
    42             end if
    43             if (igcm_generic_ice .eq. -1) then
    44                 write(*,*) "ERROR : You set a vap traceur but you forgot to set the corresponding ice traceur, &
    45                 or the pair vap/ice is not written one after another in traceur.def"
     32            ! ! ! We look for the corresponding ice traceur (before or after in the list of traceurs, maybe could be generalised to the whole list)
     33            ! ! if (iq .ne. nq) then
     34            ! !     if ((noms(iq)(1:len(trim(noms(iq)))-4) .eq. noms(iq+1)(1:len(trim(noms(iq+1)))-4)) .and. (index(noms(iq+1),"ice") .ne. 0)) then
     35            ! !         igcm_generic_ice = iq+1
     36            ! !     end if
     37            ! ! end if
     38            ! ! if ((iq .gt. 1)) then
     39            ! !     if ((noms(iq)(1:len(trim(noms(iq)))-4) .eq. noms(iq-1)(1:len(trim(noms(iq-1)))-4)) .and. (index(noms(iq-1),"ice") .ne. 0)) then
     40            ! !             igcm_generic_ice = iq-1
     41            ! !     end if
     42            ! ! end if
     43            ! ! if (igcm_generic_ice .eq. -1) then
     44            ! !     write(*,*) "ERROR : You set a vap traceur but you forgot to set the corresponding ice traceur, &
     45            ! !     or the pair vap/ice is not written one after another in traceur.def"
     46            ! ! endif
     47
     48            ! ! call_ice_vap_generic = .true.
     49            i_shortname = index(noms(iq),'_')-1
     50            shortname=noms(iq)(1:i_shortname)
     51            ! We look for the corresponding ice traceur, first if it's before the vap traceur in pq
     52            if (iq>1 ) then
     53                do ii=1,iq-1
     54                    if (index(noms(ii),adjustl(trim(shortname))) .ne. 0) then
     55                        igcm_generic_ice=ii
     56                        exit
     57                    endif
     58                ENDDO
    4659            endif
    47 
     60            ! if we didn't find it before, we look after the vap tracer in pq
     61            if (igcm_generic_ice .eq. -1 .and. (iq<nq)) then
     62                do ii=iq+1,nq
     63                    if (index(noms(ii),adjustl(trim(shortname))) .ne. 0) then
     64                        igcm_generic_ice=ii
     65                        exit
     66                    endif
     67                enddo
     68            endif
    4869            call_ice_vap_generic = .true.
    49        
    5070        else
    5171

  • trunk/LMDZ.GENERIC/libf/phystd/iniaerosol.F

    r2803 r2804  
    44      use radinc_h, only: naerkind
    55      use aerosol_mod
     6      use tracer_h, only: n_rgcs, nqtot, is_rgcs
    67      use callkeys_mod, only: aeroco2,aeroh2o,dusttau,aeroh2so4,
    78     &          aeroback2lay,aeronh3, nlayaero, aeronlay, aeroaurora,
    … …  
    2021c=======================================================================
    2122
    22       integer i, ia
     23      integer i, ia, iq
    2324
    2425      ! Special case, dyn. allocation for n-layer depending on callphys.def
    2526      IF(.NOT.ALLOCATED(iaero_nlay)) ALLOCATE(iaero_nlay(nlayaero))
    2627      iaero_nlay(:) = 0
    27       ! Do the same for iaero_generic
     28      ! Do the same for iaero_generic and i_rgcs_ice
    2829      IF (.not. allocated(iaero_generic))
    2930     &  allocate(iaero_generic(aerogeneric))
     31      if (.not. allocated(i_rgcs_ice))
     32     & allocate(i_rgcs_ice(aerogeneric))
     33      ! Init of i_rgcs_ice
     34      i_rgcs_ice(:) =0.0
     35      ia = 1
     36      do iq=1,nqtot
     37         if (is_rgcs(iq) .eq. 1) then
     38            i_rgcs_ice(ia)=iq
     39            ia = ia+1
     40         endif
     41      enddo
     42      ! We check that we have the right number
     43      ! for aerogeneric (coherent between
     44      ! callphys and traceur.def)
     45   !    if (aerogeneric .ne. n_rgcs) then
     46   !       print*,'iniaerosol: aerogeneric is not'
     47   !   & 'equal to the number of RGCS in traceur.def'
     48   !       print*,'iniaerosol: check your'
     49   !   &'callphys.def and your traceur.def'
     50   !       call abort_physics("inaerosol: Aborting")
     51   !    endif
    3052      iaero_generic(:)=0
    3153      ia=0
    … …  
    86108         enddo
    87109      endif
    88       write(*,*)'--- Generic Condensable Aerosols',iaero_generic
     110      write(*,*)'--- Radiative Generic Condensable Species = '
     111     &,iaero_generic
    89112      write(*,*) '=== Number of aerosols= ', ia
    90113

  • trunk/LMDZ.GENERIC/libf/phystd/initracer.F90

    r2785 r2804  
    120120       ENDIF
    121121       IF (.NOT. allocated(is_condensable)) allocate(is_condensable(nq)) !LT
     122       IF (.NOT. allocated(is_rgcs)) allocate(is_rgcs(nq)) !LT
    122123       IF (.NOT. allocated(constants_mass)) allocate(constants_mass(nq))
    123124       IF (.NOT. allocated(constants_delta_vapH)) allocate(constants_delta_vapH(nq))
    … …  
    148149
    149150       is_condensable(:)= 0
    150 
     151       is_rgcs(:) = 0
    151152       constants_mass(:)=0
    152153       constants_delta_vapH(:)=0
    … …  
    489490      ngt = sum(is_condensable)
    490491      write(*,*) 'Number of generic tracer is  ngt = ',ngt
     492
     493      ! Calculate number of radiative generic condensable species
     494      n_rgcs = sum(is_rgcs)
     495      write(*,*)'Number of Radiative Generic Condensable Species is n_rgcs = ',n_rgcs
     496      if (n_rgcs> ngt/2) then
     497        write(*,*) 'You have more Radiative Generic Condensable Species than Generic Condensable Species'
     498        write(*,*)'This is not possible: check your Modern traceur.def'
     499        call abort_physic("initracer, issue with # of RGCS and GCS")
     500      endif
     501
    491502!     Processing modern traceur options
    492503      if(moderntracdef) then
    … …  
    693704              ' rho=', rho_q(iq)       
    694705          endif
     706          !option is_rgcs
     707          if (index(tracline,'is_rgcs') .ne. 0) then
     708            read(tracline(index(tracline,'is_rgcs=') &
     709              +len('is_rgcs='):),*) is_rgcs(iq)
     710            write(*,*)'Parameter value (traceur.def) :'// &
     711              'is_rgcs=',is_rgcs(iq)
     712          else
     713            write(*,*)'Parameter value (default) : '// &
     714              'is_rgcs = ',is_rgcs(iq)
     715          endif
    695716      end subroutine get_tracdat
    696717

  • trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90

    r2803 r2804  
    3030!     Jeremy Leconte (2012)
    3131!     Extended to dust, CO2, NH3, 2-lay,Nlay,auroral aerosols by ??
    32 !     Added Radiative Generic Condensable Aerosols effective radii
     32!     Added Radiative Generic Condensable Species effective radii
    3333!     calculations  (Lucas Teinturier 2022)
    3434!
    … …  
    4242      use aerosol_mod, only: iaero_back2lay, iaero_co2, iaero_dust, &
    4343                             iaero_h2o, iaero_h2so4, iaero_nh3, iaero_nlay, iaero_aurora, &
    44                              iaero_generic
     44                             iaero_generic, i_rgcs_ice
    4545      use callkeys_mod, only: size_nh3_cloud, nlayaero, aeronlay_size, aeronlay_nueff,aerogeneric
    46 
     46      use tracer_h, only: radius, nqtot, is_rgcs
    4747      Implicit none
    4848
    … …  
    5555      logical, save :: firstcall=.true.
    5656!$OMP THREADPRIVATE(firstcall)
    57       integer :: iaer, ia  
    58       real :: temprad !temporary variable for radius when we have Radiative Generic Condensable aerosols
     57      integer :: iaer, ia , iq, i_rad 
     58
    5959      do iaer=1,naerkind
    6060!     these values will change once the microphysics gets to work
    … …  
    107107         endif
    108108
    109          do ia=1,aerogeneric     ! Radiative Generic Condensable aerosols
     109         do ia=1,aerogeneric     ! Radiative Generic Condensable Species
    110110            if (iaer .eq. iaero_generic(ia)) then
    111                call generic_get_traceur_radius(iaer,temprad)
    112                reffrad(1:ngrid,1:nlayer,iaer)=temprad
     111               i_rad = i_rgcs_ice(ia)
     112               reffrad(1:ngrid,1:nlayer,iaer)=radius(i_rad)
    113113               nueffrad(1:ngrid,1:nlayer,iaer) = 0.1
    114114            endif
    115115         enddo  ! generic radiative condensable aerosols
    116 
     116         
    117117         if(iaer.gt.5)then ! I think this could be commented out by this needs testing. It works with 5 aerosols (LT 2022)
    118118            print*,'Error in callcorrk, naerkind is too high (>5).'
    … …  
    122122         endif
    123123      enddo ! iaer=1,naerkind
    124 
     124     
    125125
    126126      if (radfixed) then
    … …  
    390390!==================================================================
    391391
    392 !==================================================================
    393    subroutine generic_get_traceur_radius(iaer,radii)
    394 !==================================================================
    395 !     Purpose
    396 !     -------
    397 !     Get the generic traceur radius inputted in modern_traceur.def
    398 !      for the specie "specie"
    399 !
    400 !     Authors
    401 !     -------
    402 !     Lucas Teinturier (2022)
    403 !
    404 !==================================================================
    405       Use tracer_h, only: radius, noms, nqtot
    406       Implicit none
    407      
    408       integer, intent(in) :: iaer !index of aerosol which radius we're looking for
    409       real, intent(out) :: radii !radii(ngrid,nlayer,naerkind) ! Particle radii (m)
    410       integer :: iq !just a loop index
    411 
    412       if (index(noms(1),'ice') .ne. 0) then !traceur.def is organized with _ice then _vap
    413          radii=radius(2*iaer-1)
    414       elseif (index(noms(1),'vap') .ne. 0) then !traceur.def is organized with _vap then _ice
    415          radii = radius(2*iaer)
    416       else
    417          print*,"radii_mod: tracer is neither ice or vap :",noms(1)
    418          print*,"Aborting"
    419          call abort
    420       endif
    421      
    422    end subroutine generic_get_traceur_radius
    423 
    424392end module radii_mod
    425393!==================================================================

  • trunk/LMDZ.GENERIC/libf/phystd/suaer_corrk.F90

    r2803 r2804  
    233233       endif
    234234! the following was added by LT
    235        do ia=1,aerogeneric ! Read Radiative Generic Condensable Aerosols data
     235       do ia=1,aerogeneric ! Read Radiative Generic Condensable Species data
    236236         if (iaer .eq. iaero_generic(ia)) then
    237             if (index(noms(2*iaer),'Fe') .ne. 0) then
     237            if (index(noms(i_rgcs_ice(ia)),'Fe') .ne. 0) then
    238238               print*,"Reading Fe file"
    239239               file_id(iaer,1)='Fe.ocst.txt'
    … …  
    241241               lamrefvis(iaer) = 1.0E-6 !random pick
    242242               lamrefir(iaer) = 1.0E-6 !dummy but random pick
    243             else if (index(noms(2*iaer),'Mn') .ne. 0) then
     243            else if (index(noms(i_rgcs_ice(ia)),'Mn') .ne. 0) then
    244244               print*,"Reading MnS file"
    245245               file_id(iaer,1)='MnS.ocst.txt'
    … …  
    247247               lamrefvis(iaer) = 1.0E-6 !random pick
    248248               lamrefir(iaer) = 1.0E-6 !dummy but random pick   
    249             else if (index(noms(2*iaer),'Mg') .ne. 0) then 
     249            else if (index(noms(i_rgcs_ice(ia)),'Mg') .ne. 0) then 
    250250               print*,"Reading Mg2SiO4 file"
    251251               file_id(iaer,1)='Mg2SiO4.ocst.txt'
    … …  
    253253               lamrefvis(iaer) = 1.0E-6 !random pick
    254254               lamrefir(iaer) = 1.0E-6 !dummy but random pick 
    255             else if (index(noms(2*iaer),'Cr') .ne. 0) then
     255            else if (index(noms(i_rgcs_ice(ia)),'Cr') .ne. 0) then
    256256               print*,"Reading Cr file"
    257257               file_id(iaer,1)='Cr.ocst.txt'
    … …  
    261261            else
    262262! If you want to add another specie, copy,paste & adapt the elseif block up here to your new specie (LT 2022)
    263                call abort_physic("suaer_corrk", "Unknown specie in generic condensable aerosols",1)
     263               call abort_physic("suaer_corrk", "Unknown specie in radiative generic condensable species",1)
    264264            endif
    265265         endif
    266266       enddo ! ia=1,aerogeneric
    267267      enddo
    268 
     268     
    269269!------------------------------------------------------------------
    270270

  • trunk/LMDZ.GENERIC/libf/phystd/tracer_h.F90

    r2785 r2804  
    77       integer, save :: nesp  ! number of species in the chemistry
    88       integer, save :: ngt   ! number of generic tracers
     9<<<<<<< HEAD
    910!$OMP THREADPRIVATE(nqtot,nesp,ngt)
     11=======
     12       integer, save :: n_rgcs ! number of Radiative Generic Condensable Species
     13!$OMP THREADPRIVATE(nqtot,nesp,ngt,n_rgcs)
     14>>>>>>> d67f7c1... This is the RGCS scheme
    1015
    1116       logical :: moderntracdef=.false. ! Standard or modern traceur.def
    … …  
    4146!$OMP THREADPRIVATE(is_recomb,is_recomb_qset,is_recomb_qotf)
    4247       integer, save, allocatable :: is_condensable(:)      ! 1 if tracer is generic, else 0 (added LT)
     48<<<<<<< HEAD
    4349!$OMP THREADPRIVATE(is_condensable)   !also added by LT
    4450
     51=======
     52       integer,save,allocatable :: is_rgcs(:)               ! 1 if tracer is a radiative generic condensable specie, else 0 (added LT 2022)
     53>>>>>>> d67f7c1... This is the RGCS scheme
    4554       ! Lists of constants for condensable tracers
    4655       integer, save, allocatable :: constants_mass(:)                    ! molecular mass of the specie (g/mol)
    … …  
    5564!$OMP THREADPRIVATE(constants_RLVTT_generic,constants_metallicity_coeff)
    5665
     66<<<<<<< HEAD
     67=======
     68!$OMP THREADPRIVATE(is_condensable,is_rgcs)   !also added by LT
     69>>>>>>> d67f7c1... This is the RGCS scheme
    5770! tracer indexes: these are initialized in initracer and should be 0 if the
    5871!                 corresponding tracer does not exist
Note: See TracChangeset for help on using the changeset viewer.