Ignore:
Timestamp:
Nov 23, 2022, 4:41:34 PM (2 years ago)
Author:
emillour
Message:

Generic PCM:
Add the possibility to include Venus-like aerosols (triggered by option
aerovenus=.true. in callphys.def); baseline is to use 5 distinct scatterers
but each may be turned on/off (via aerovenus1, aerovenus2, aerovenus2p,
aerovenus3, aerovenusUV flags which may be specified in callphys.def).
GG

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r2810 r2831  
    1       Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pq, &
    2          aerosol,reffrad,QREFvis3d,QREFir3d,tau_col,cloudfrac,totcloudfrac,clearsky)
     1      Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pt, pq, &
     2         aerosol,reffrad,nueffrad, QREFvis3d,QREFir3d,tau_col, &
     3         cloudfrac,totcloudfrac,clearsky)
    34
    45       use radinc_h, only : L_TAUMAX,naerkind
    … …  
    67                              iaero_aurora, iaero_back2lay, iaero_co2, &
    78                              iaero_dust, iaero_h2o, iaero_h2so4, &
    8                               iaero_nh3, i_rgcs_ice, noaero
     9                              iaero_nh3, i_rgcs_ice, noaero, &
     10                              iaero_venus1, iaero_venus2, iaero_venus2p, &
     11                              iaero_venus3, iaero_venusUV
    912       USE tracer_h, only: noms,rho_co2,rho_ice,rho_q,mmol
    1013       use comcstfi_mod, only: g, pi, mugaz, avocado
    … …  
    6063      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
    6164      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
     65      REAL,INTENT(IN) :: pt(ngrid,nlayer) ! mid-layer temperature (K)
    6266      REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! aerosol optical depth
    6367      REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind) ! aerosol effective radius
     68      REAL,INTENT(IN) :: nueffrad(ngrid,nlayer,naerkind) ! aerosol effective variance
    6469      REAL,INTENT(IN) :: QREFvis3d(ngrid,nlayer,naerkind) ! extinction coefficient in the visible
    6570      REAL,INTENT(IN) :: QREFir3d(ngrid,nlayer,naerkind)
    … …  
    98103      logical dummy_bool ! dummy boolean just in case we need one
    99104      ! integer i_rgcs_ice(aerogeneric)
     105      !  for venus clouds
     106      real      :: p_bot,p_top,h_bot,h_top,mode_dens,h_lay
     107
    100108      ! identify tracers
    101109      IF (firstcall) THEN
    … …  
    148156          print*,'iaero_aurora= ',iaero_aurora
    149157        endif
     158
     159        if (iaero_venus1.ne.0) then
     160          print*,'iaero_venus1= ',iaero_venus1
     161        endif
     162        if (iaero_venus2.ne.0) then
     163          print*,'iaero_venus2= ',iaero_venus2
     164        endif
     165        if (iaero_venus2p.ne.0) then
     166          print*,'iaero_venus2p= ',iaero_venus2p
     167        endif
     168        if (iaero_venus3.ne.0) then
     169          print*,'iaero_venus3= ',iaero_venus3
     170        endif
     171        if (iaero_venusUV.ne.0) then
     172          print*,'iaero_venusUV= ',iaero_venusUV
     173        endif
     174
    150175        if (aerogeneric .ne. 0) then
    151176          print*,"iaero_generic= ",iaero_generic(:)
    … …  
    679704      endif !aerogeneric .ne. 0
    680705
     706!==================================================================
     707!         Venus clouds (4 modes)
     708!   S. Lebonnois (jan 2016)
     709!==================================================================
     710! distributions from Haus et al, 2013
     711! mode             1      2      2p     3
     712! r (microns)     0.30   1.05   1.40   3.65
     713! sigma           1.56   1.29   1.23   1.28
     714! reff (microns)  0.49   1.23   1.56   4.25
     715! nueff           0.21   0.067  0.044  0.063
     716! (nueff=exp(ln^2 sigma)-1)
     717!
     718! p_bot <=> zb ; p_top <=> zb+zc ; h_bot <=> Hlo ; h_top <=> Hup
     719! p<p_top: N=No*(p/p_top)**(h_lay/h_top)      h_lay=RT/g  (in m)
     720! p>p_bot: N=No*(p_bot/p)**(h_lay/h_bot)      R=8.31/mu (mu in kg/mol)
     721! N is in m-3
     722!
     723! dTau = Qext*[pi*reff**2*exp(-3*ln(1+nueff))]*N*h_lay*(-dp)/p
     724
     725! Mode 1
     726      if (iaero_venus1 .ne.0) then
     727          iaer=iaero_venus1
     728
     729!       1. Initialization
     730          aerosol(1:ngrid,1:nlayer,iaer)=0.0
     731          p_bot = 1.e5
     732          p_top = 1.e4
     733          h_bot = 1.0e3 ! m
     734          h_top = 5.0e3
     735         
     736!       2. Opacity calculation
     737
     738          DO ig=1,ngrid
     739           DO l=1,nlayer-1
     740
     741             h_lay=8.31*pt(ig,l)/(g*0.044)
     742
     743             !! 1. below 2e5 Pa: no aerosols
     744             IF (pplay(ig,l) .gt. 2.e5) THEN
     745               mode_dens = 0.
     746
     747             !! 2. cloud layer
     748             ELSEIF (pplay(ig,l) .le. 2.e5 .and. pplay(ig,l) .gt. p_bot) THEN
     749               mode_dens = 1.81e8*(p_bot/pplay(ig,l))**(h_lay/h_bot)
     750               
     751             ELSEIF (pplay(ig,l) .le. p_bot .and. pplay(ig,l) .gt. p_top) THEN
     752               mode_dens = 1.81e8  ! m-3
     753               
     754             ELSEIF (pplay(ig,l) .le. p_top .and. pplay(ig,l) .gt. 1.e2) THEN
     755               mode_dens = 1.81e8*(pplay(ig,l)/p_top)**(h_lay/h_top)
     756               
     757             !! 3. above 1.e2 Pa: no aerosols
     758             ELSEIF (pplay(ig,l) .le. 1.e2) THEN
     759               mode_dens = 0.
     760             ENDIF
     761
     762             aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
     763              pi*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* &
     764              mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
     765
     766           ENDDO
     767          ENDDO
     768
     769      end if ! mode 1
     770
     771! Mode 2
     772      if (iaero_venus2 .ne.0) then
     773          iaer=iaero_venus2
     774
     775!       1. Initialization
     776          aerosol(1:ngrid,1:nlayer,iaer)=0.0
     777          p_bot = 1.1e4
     778          p_top = 1.e4
     779          h_bot = 3.0e3
     780          h_top = 3.5e3
     781         
     782!       2. Opacity calculation
     783
     784          DO ig=1,ngrid
     785           DO l=1,nlayer-1
     786
     787             h_lay=8.31*pt(ig,l)/(g*0.044)
     788
     789             !! 1. below 2e5 Pa: no aerosols
     790             IF (pplay(ig,l) .gt. 2.e5) THEN
     791               mode_dens = 0.
     792
     793             !! 2. cloud layer
     794             ELSEIF (pplay(ig,l) .le. 2.e5 .and. pplay(ig,l) .gt. p_bot) THEN
     795               mode_dens = 1.00e8*(p_bot/pplay(ig,l))**(h_lay/h_bot)
     796               
     797             ELSEIF (pplay(ig,l) .le. p_bot .and. pplay(ig,l) .gt. p_top) THEN
     798               mode_dens = 1.00e8
     799               
     800             ELSEIF (pplay(ig,l) .le. p_top .and. pplay(ig,l) .gt. 1.e2) THEN
     801               mode_dens = 1.00e8*(pplay(ig,l)/p_top)**(h_lay/h_top)
     802               
     803             !! 3. above 1.e2 Pa: no aerosols
     804             ELSEIF (pplay(ig,l) .le. 1.e2) THEN
     805               mode_dens = 0.
     806             ENDIF
     807
     808             aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
     809              pi*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* &
     810              mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
     811
     812           ENDDO
     813          ENDDO
     814
     815      end if ! mode 2
     816
     817! Mode 2p
     818      if (iaero_venus2p .ne.0) then
     819          iaer=iaero_venus2p
     820
     821!       1. Initialization
     822          aerosol(1:ngrid,1:nlayer,iaer)=0.0
     823          p_bot = 1.e5
     824          p_top = 2.3e4
     825          h_bot = 0.1e3
     826          h_top = 1.0e3
     827         
     828!       2. Opacity calculation
     829
     830          DO ig=1,ngrid
     831           DO l=1,nlayer-1
     832
     833             h_lay=8.31*pt(ig,l)/(g*0.044)
     834
     835             !! 1. below 2e5 Pa: no aerosols
     836             IF (pplay(ig,l) .gt. 2.e5) THEN
     837               mode_dens = 0.
     838
     839             !! 2. cloud layer
     840             ELSEIF (pplay(ig,l) .le. 2.e5 .and. pplay(ig,l) .gt. p_bot) THEN
     841               mode_dens = 5.00e7*(p_bot/pplay(ig,l))**(h_lay/h_bot)
     842               
     843             ELSEIF (pplay(ig,l) .le. p_bot .and. pplay(ig,l) .gt. p_top) THEN
     844               mode_dens = 5.00e7
     845               
     846             ELSEIF (pplay(ig,l) .le. p_top .and. pplay(ig,l) .gt. 1.e2) THEN
     847               mode_dens = 5.00e7*(pplay(ig,l)/p_top)**(h_lay/h_top)
     848               
     849             !! 3. above 1.e2 Pa: no aerosols
     850             ELSEIF (pplay(ig,l) .le. 1.e2) THEN
     851               mode_dens = 0.
     852             ENDIF
     853
     854             aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
     855              pi*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* &
     856              mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
     857
     858           ENDDO
     859          ENDDO
     860
     861      end if ! mode 2p
     862
     863! Mode 3
     864      if (iaero_venus3 .ne.0) then
     865          iaer=iaero_venus3
     866
     867!       1. Initialization
     868          aerosol(1:ngrid,1:nlayer,iaer)=0.0
     869          p_bot = 1.e5
     870          p_top = 4.e4
     871          h_bot = 0.5e3
     872          h_top = 1.0e3
     873         
     874!       2. Opacity calculation
     875
     876          DO ig=1,ngrid
     877           DO l=1,nlayer-1
     878 
     879              h_lay=8.31*pt(ig,l)/(g*0.044)
     880
     881             !! 1. below 2e5 Pa: no aerosols
     882             IF (pplay(ig,l) .gt. 2.e5) THEN
     883               mode_dens = 0.
     884
     885             !! 2. cloud layer
     886             ELSEIF (pplay(ig,l) .le. 2.e5 .and. pplay(ig,l) .gt. p_bot) THEN
     887               mode_dens = 1.40e7*(p_bot/pplay(ig,l))**(h_lay/h_bot)
     888               
     889             ELSEIF (pplay(ig,l) .le. p_bot .and. pplay(ig,l) .gt. p_top) THEN
     890               mode_dens = 1.40e7
     891               
     892             ELSEIF (pplay(ig,l) .le. p_top .and. pplay(ig,l) .gt. 1.e2) THEN
     893               mode_dens = 1.40e7*(pplay(ig,l)/p_top)**(h_lay/h_top)
     894               
     895             !! 3. above 1.e2 Pa: no aerosols
     896             ELSEIF (pplay(ig,l) .le. 1.e2) THEN
     897               mode_dens = 0.
     898             ENDIF
     899
     900             aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*                       &
     901              pi*(reffrad(ig,l,iaer))**2*exp(-3.*log(1+nueffrad(ig,l,iaer)))* &
     902              mode_dens*h_lay*(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
     903
     904           ENDDO
     905          ENDDO
     906
     907      end if ! mode 3
     908
     909! UV absorber
     910      if (iaero_venusUV .ne.0) then
     911          iaer=iaero_venusUV
     912
     913!       1. Initialization
     914          aerosol(1:ngrid,1:nlayer,iaer)=0.0
     915          p_bot = 3.3e4  ! 58 km
     916          p_top = 3.7e3 ! 70 km
     917          h_bot = 1.0e3
     918          h_top = 1.0e3
     919         
     920!       2. Opacity calculation
     921
     922          DO ig=1,ngrid
     923           DO l=1,nlayer-1
     924
     925             h_lay=8.31*pt(ig,l)/(g*0.044)
     926
     927             !! 1. below 7e4 Pa: no aerosols
     928             IF (pplay(ig,l) .gt. 7.e4) THEN
     929               mode_dens = 0.
     930
     931             !! 2. cloud layer
     932             ELSEIF (pplay(ig,l) .le. 7.e4 .and. pplay(ig,l) .gt. p_bot) THEN
     933               mode_dens = 1.00e7*(p_bot/pplay(ig,l))**(h_lay/h_bot)
     934               
     935             ELSEIF (pplay(ig,l) .le. p_bot .and. pplay(ig,l) .gt. p_top) THEN
     936               mode_dens = 1.00e7
     937               
     938             ELSEIF (pplay(ig,l) .le. p_top .and. pplay(ig,l) .gt. 1.e3) THEN
     939               mode_dens = 1.00e7*(pplay(ig,l)/p_top)**(h_lay/h_top)
     940               
     941             !! 3. above 1.e3 Pa: no aerosols
     942             ELSEIF (pplay(ig,l) .le. 1.e3) THEN
     943               mode_dens = 0.
     944             ENDIF
     945
     946! normalized to 0.35 microns (peak of absorption)
     947             aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*mode_dens
     948
     949           ENDDO
     950          ENDDO
     951
     952!       3. Re-normalize to Haus et al 2015 total column optical depth
     953         tau_col(:)=0.0
     954         DO l=1,nlayer
     955          DO ig=1,ngrid
     956               tau_col(ig) = tau_col(ig) &
     957                     + aerosol(ig,l,iaer)
     958            ENDDO
     959         ENDDO
     960         DO ig=1,ngrid
     961           DO l=1,nlayer-1
     962                aerosol(ig,l,iaer)=aerosol(ig,l,iaer)*0.205/tau_col(ig)
     963           ENDDO
     964         ENDDO
     965
     966      end if ! UV absorber
     967
     968!==================================================================
     969!     ig=10
     970!      do l=1,nlayer
     971!          print*,8.31*pt(ig,l)/(g*0.044),pplay(ig,l),aerosol(ig,l,1),aerosol(ig,l,2),aerosol(ig,l,3),aerosol(ig,l,4)
     972!         print*,l,pplay(ig,l),aerosol(ig,l,5)
     973!      enddo
     974!      stop           
     975!==================================================================
     976
     977
    681978! --------------------------------------------------------------------------
    682979! Column integrated visible optical depth in each point (used for diagnostic)
    … …  
    7131010      !    endif
    7141011      ! end do
    715       return
     1012
    7161013    end subroutine aeropacity
    7171014     
Note: See TracChangeset for help on using the changeset viewer.