module aeropacity_mod implicit none contains Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pt, pq, & aerosol,reffrad,nueffrad, QREFvis3d,QREFir3d,tau_col, & cloudfrac,totcloudfrac,clearsky) use radinc_h, only : L_TAUMAX,naerkind use aerosol_mod, only: iaero_nlay, iaero_generic, & iaero_aurora, iaero_back2lay, iaero_n2, & iaero_dust, iaero_h2so4, & iaero_nh3, i_rgcs_ice, noaero USE tracer_h, only: noms,rho_n2,rho_ice,rho_q,mmol use comcstfi_mod, only: g, pi, mugaz, avocado use geometry_mod, only: latitude use callkeys_mod, only: aerofixn2,kastprof,cloudlvl, & dusttau, & pres_bottom_tropo,pres_top_tropo,obs_tau_col_tropo, & pres_bottom_strato,pres_top_strato,obs_tau_col_strato, & tau_nh3_cloud, pres_nh3_cloud, & nlayaero, aeronlay_tauref, aeronlay_choice, & aeronlay_pbot, aeronlay_ptop, aeronlay_sclhght, & aerogeneric use generic_tracer_index_mod, only: generic_tracer_index implicit none !================================================================== ! ! Purpose ! ------- ! Compute aerosol optical depth in each gridbox. ! ! Authors ! ------- ! F. Forget ! F. Montmessin (water ice scheme) ! update J.-B. Madeleine (2008) ! dust removal, simplification by Robin Wordsworth (2009) ! Generic n-layer aerosol - J. Vatant d'Ollone (2020) ! Radiative Generic Condensable Species - Lucas Teinturier (2022) ! ! Input ! ----- ! ngrid Number of horizontal gridpoints ! nlayer Number of layers ! nq Number of tracers ! pplev Pressure (Pa) at each layer boundary ! pq Aerosol mixing ratio ! reffrad(ngrid,nlayer,naerkind) Aerosol effective radius ! QREFvis3d(ngrid,nlayer,naerkind) \ 3d extinction coefficients ! QREFir3d(ngrid,nlayer,naerkind) / at reference wavelengths ! ! Output ! ------ ! aerosol Aerosol optical depth in layer l, grid point ig ! tau_col Total column optical depth at grid point ig ! !======================================================================= INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers INTEGER,INTENT(IN) :: nq ! number of tracers REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa) REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa) REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air) REAL,INTENT(IN) :: pt(ngrid,nlayer) ! mid-layer temperature (K) REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! aerosol optical depth REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind) ! aerosol effective radius REAL,INTENT(IN) :: nueffrad(ngrid,nlayer,naerkind) ! aerosol effective variance REAL,INTENT(IN) :: QREFvis3d(ngrid,nlayer,naerkind) ! extinction coefficient in the visible REAL,INTENT(IN) :: QREFir3d(ngrid,nlayer,naerkind) REAL,INTENT(OUT):: tau_col(ngrid) !column integrated visible optical depth ! BENJAMIN MODIFS real,intent(in) :: cloudfrac(ngrid,nlayer) ! cloud fraction real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction logical,intent(in) :: clearsky real aerosol0, obs_tau_col_aurora, pm real pcloud_deck, cloud_slope real dp_strato(ngrid) real dp_tropo(ngrid) real dp_layer(ngrid) INTEGER l,ig,iq,iaer,ia LOGICAL,SAVE :: firstcall=.true. !$OMP THREADPRIVATE(firstcall) REAL CBRT EXTERNAL CBRT INTEGER,SAVE :: i_n2ice=0 ! n2 ice !$OMP THREADPRIVATE(i_n2ice) CHARACTER(LEN=20) :: tracername ! to temporarily store text ! for fixed dust profiles real topdust, expfactor, zp REAL taudusttmp(ngrid) ! Temporary dust opacity used before scaling REAL tauh2so4tmp(ngrid) ! Temporary h2so4 opacity used before scaling real CLFtot integer igen_ice,igen_vap ! to store the index of generic tracer logical dummy_bool ! dummy boolean just in case we need one ! integer i_rgcs_ice(aerogeneric) ! for venus clouds real :: p_bot,p_top,h_bot,h_top,mode_dens,h_lay ! identify tracers IF (firstcall) THEN ia =0 write(*,*) "Tracers found in aeropacity:" do iq=1,nq tracername=noms(iq) if (tracername.eq."n2_ice") then i_n2ice=iq write(*,*) "i_n2ice=",i_n2ice endif enddo if (noaero) then print*, "No active aerosols found in aeropacity" else print*, "If you would like to use aerosols, make sure any old" print*, "start files are updated in newstart using the option" print*, "q=0" write(*,*) "Active aerosols found in aeropacity:" endif if ((iaero_n2.ne.0).and.(.not.noaero)) then print*, 'iaero_n2= ',iaero_n2 endif if (iaero_dust.ne.0) then print*,'iaero_dust= ',iaero_dust endif if (iaero_h2so4.ne.0) then print*,'iaero_h2so4= ',iaero_h2so4 endif if (iaero_back2lay.ne.0) then print*,'iaero_back2lay= ',iaero_back2lay endif if (iaero_nh3.ne.0) then print*,'iaero_nh3= ',iaero_nh3 endif if (iaero_nlay(1).ne.0) then print*,'iaero_nlay= ',iaero_nlay(:) endif if (iaero_aurora.ne.0) then print*,'iaero_aurora= ',iaero_aurora endif if (aerogeneric .ne. 0) then print*,"iaero_generic= ",iaero_generic(:) endif firstcall=.false. ENDIF ! of IF (firstcall) ! --------------------------------------------------------- !================================================================== ! N2 ice aerosols !================================================================== if (iaero_n2.ne.0) then iaer=iaero_n2 ! 1. Initialization aerosol(1:ngrid,1:nlayer,iaer)=0.0 ! 2. Opacity calculation if (noaero) then ! aerosol set to zero aerosol(1:ngrid,1:nlayer,iaer)=0.0 elseif (aerofixn2.or.(i_n2ice.eq.0)) then ! N2 ice cloud prescribed aerosol(1:ngrid,1:nlayer,iaer)=1.e-9 !aerosol(1:ngrid,12,iaer)=4.0 ! single cloud layer option else DO ig=1, ngrid DO l=1,nlayer-1 ! to stop the rad tran bug aerosol0 = & ( 0.75 * QREFvis3d(ig,l,iaer) / & ( rho_n2 * reffrad(ig,l,iaer) ) ) * & ( pq(ig,l,i_n2ice) + 1.E-9 ) * & ( pplev(ig,l) - pplev(ig,l+1) ) / g aerosol0 = max(aerosol0,1.e-9) aerosol0 = min(aerosol0,L_TAUMAX) aerosol(ig,l,iaer) = aerosol0 ! aerosol(ig,l,iaer) = 0.0 ! print*, aerosol(ig,l,iaer) ! using cloud fraction ! aerosol(ig,l,iaer) = -log(1 - CLF + CLF*exp(-aerosol0/CLF)) ! aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),L_TAUMAX) ENDDO ENDDO end if ! if fixed or varying end if ! if N2 aerosols !================================================================== ! Water ice / liquid (AF24: removed) !================================================================== !================================================================== ! Dust (AF24: removed) !================================================================== !================================================================== ! H2SO4 (AF24: removed) !================================================================== !================================================================== ! Two-layer aerosols (unknown composition) ! S. Guerlet (2013) - Modif by J. Vatant d'Ollone (2020) ! ! This scheme is deprecated and left for retrocompatibility ! You should use the n-layer scheme below ! ! !================================================================== !================================================================== ! Saturn/Jupiter ammonia cloud = thin cloud (scale height 0.2 hard coded...) ! S. Guerlet (2013) ! JVO 20 : You should now use the generic n-layer scheme below !================================================================== !========================================================================================================= ! Generic N-layers aerosols/clouds ! Author : J. Vatant d'Ollone (2020) ! ! Purpose: Replaces and extents the former buggy 2-layer scheme as well as hard-coded NH3 cloud ! ! + Each layer can have different optical properties, size of particle ... ! + Enables up to n=4 layers as we apparently cannot run with more scatterers (could be worth checking...) ! + You have different choices for vertical profile of the aerosol layers : ! * aeronlay_choice = 1 : Layer tau is spread between ptop and pbot following atm scale height. ! * aeronlay_choice = 2 : Layer tau follows its own scale height above cloud deck (pbot). ! In this case ptop is dummy and sclhght gives the ratio H_cl/H_atm. ! * aeronlay_choice = ... feel free to add more cases ! ! + Layers can overlap if needed (if you want a 'transition layer' as in the 2-scheme, just add it) ! !========================================================================================================= if (iaero_nlay(1) .ne.0) then DO ia=1,nlayaero iaer=iaero_nlay(ia) ! a. Initialization aerosol(1:ngrid,1:nlayer,iaer)=0.D0 ! b. Opacity calculation ! Case 1 : Follows atmospheric scale height between boundaries pressures ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF (aeronlay_choice(ia).eq.1) THEN dp_layer(:)=0.D0 DO ig=1,ngrid DO l=1,nlayer-1 !! i. Opacity follows scale height IF ( pplev(ig,l).le.aeronlay_pbot(ia) .AND. & pplev(ig,l).ge.aeronlay_ptop(ia) ) THEN aerosol(ig,l,iaer) = ( pplev(ig,l) - pplev(ig,l+1) ) dp_layer(ig) = dp_layer(ig) + aerosol(ig,l,iaer) !! ii. Outside aerosol layer boundaries: no aerosols ELSE aerosol(ig,l,iaer) = 0.D0 ENDIF ENDDO ENDDO ! iii. Re-normalize to required total opacity DO ig=1,ngrid DO l=1,nlayer-1 IF ( pplev(ig,l).le.aeronlay_pbot(ia) .AND. & pplev(ig,l).ge.aeronlay_ptop(ia) ) THEN aerosol(ig,l,iaer) = aerosol(ig,l,iaer) / dp_layer(ig) & * aeronlay_tauref(ia) ENDIF ENDDO ENDDO ! Case 2 : Follows input scale height ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ELSE IF (aeronlay_choice(ia).eq.2) THEN cloud_slope = 1.D0/aeronlay_sclhght(ia) pcloud_deck = aeronlay_pbot(ia) dp_layer(:) = 0.D0 DO ig=1,ngrid DO l=1,nlayer-1 !! i. Below cloud layer: no opacity IF (pplev(ig,l) .gt. pcloud_deck) THEN aerosol(ig,l,iaer) = 0.D0 !! ii. Follows scale height above cloud deck ELSEIF (pplev(ig,l) .le. pcloud_deck) THEN aerosol(ig,l,iaer) = ((pplev(ig,l)/pcloud_deck)**(cloud_slope)) dp_layer(ig) = dp_layer(ig) + aerosol(ig,l,iaer) ENDIF ENDDO ENDDO ! iii. Re-normalize to required total opacity DO ig=1,ngrid DO l=1,nlayer-1 IF (pplev(ig,l) .le. pcloud_deck) THEN aerosol(ig,l,iaer) = aerosol(ig,l,iaer) / dp_layer(ig) & * aeronlay_tauref(ia) ENDIF ENDDO ENDDO ENDIF ! aeronlay_choice ENDDO ! loop on n aerosol layers end if ! if N-layer aerosols !================================================================== ! Jovian auroral aerosols (unknown composition) NON-GENERIC: vertical and meridional profile tuned to observations ! S. Guerlet (2015) !================================================================== if (iaero_aurora .ne.0) then iaer=iaero_aurora ! 1. Initialization aerosol(1:ngrid,1:nlayer,iaer)=0.D0 pm = 2000. !!case study: maxi aerosols at 20 hPa ! 2. Opacity calculation DO ig=1,ngrid !! Test Jupiter (based on Zhang et al 2013 observations, but a bit different), decembre 2015 DO l=1,nlayer aerosol(ig,l,iaer) = (pplev(ig,l)/pm)**2 * exp(-(pplev(ig,l)/pm)**2) ENDDO ENDDO ! 3. Meridional distribution, and re-normalize to observed total column dp_layer(:)=0.D0 DO ig=1,ngrid !!Jupiter !!Hem sud: IF (latitude(ig)*180.D0/pi .lt. -45.D0 .and. latitude(ig)*180.D0/pi .gt. -70.) THEN obs_tau_col_aurora= 10.D0**(-0.06D0*latitude(ig)*180.D0/pi-3.4D0) ELSEIF (latitude(ig)*180.D0/pi .lt. -37.D0 .and. latitude(ig)*180.D0/pi .ge. -45.) THEN obs_tau_col_aurora= 10.D0**(-0.3D0*latitude(ig)*180.D0/pi-14.3D0) ELSEIF (latitude(ig)*180./pi .le. -70. ) THEN obs_tau_col_aurora= 10**(0.06*70.-3.4D0) !!Hem Nord: ELSEIF (latitude(ig)*180.D0/pi .gt. 30.D0 .and. latitude(ig)*180.D0/pi .lt. 70.) THEN obs_tau_col_aurora= 10.D0**(0.03D0*latitude(ig)*180.D0/pi-1.17D0) ELSEIF (latitude(ig)*180.D0/pi .gt. 22.D0 .and. latitude(ig)*180.D0/pi .le. 30.) THEN obs_tau_col_aurora= 10.D0**(0.3D0*latitude(ig)*180.D0/pi-9.4D0) ELSEIF (latitude(ig)*180.D0/pi .ge. 70.) THEN obs_tau_col_aurora= 10**(0.03*70.-1.17D0) ELSEIF (latitude(ig)*180.D0/pi .ge. -37. .and. latitude(ig)*180.D0/pi .le. 22.) THEN obs_tau_col_aurora = 0.001D0 !!Jupiter: mini pas a zero ENDIF DO l=1,nlayer dp_layer(ig) = dp_layer(ig) + aerosol(ig,l,iaer)/obs_tau_col_aurora ENDDO ENDDO DO ig=1,ngrid DO l=1,nlayer-1 aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/dp_layer(ig) ENDDO ENDDO end if ! if Auroral aerosols !=========================================================================== ! Radiative Generic Condensable aerosols scheme ! Only used when we give aerogeneric != 0 in callphys.def ! Computes the generic aerosols' opacity in the same fashion as water of ! dust, using the QREFvis3d of the concerned specie ! Lucas Teinturier (2022) !=========================================================================== if (aerogeneric .ne. 0) then ! we enter the scheme do ia=1,aerogeneric iaer = iaero_generic(ia) ! Initialization aerosol(1:ngrid,1:nlayer,iaer) = 0.D0 igen_ice = i_rgcs_ice(ia) ! Let's loop on the horizontal and vertical grid do ig=1,ngrid do l=1,nlayer aerosol(ig,l,iaer) = ( 0.75*QREFvis3d(ig,l,iaer) / & (rho_q(igen_ice) * reffrad(ig,l,iaer)) ) * & (pq(ig,l,igen_ice)+1E-9 ) * & (pplev(ig,l) - pplev(ig,l+1)) /g enddo !l=1,nlayer enddo !ig=1,ngrid enddo !ia=1,aerogeneric endif !aerogeneric .ne. 0 ! -------------------------------------------------------------------------- ! Column integrated visible optical depth in each point (used for diagnostic) tau_col(:)=0.0 do iaer = 1, naerkind do l=1,nlayer do ig=1,ngrid tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer) end do end do end do ! do ig=1,ngrid ! do l=1,nlayer ! do iaer = 1, naerkind ! if(aerosol(ig,l,iaer).gt.1.e3)then ! print*,'WARNING: aerosol=',aerosol(ig,l,iaer) ! print*,'at ig=',ig,', l=',l,', iaer=',iaer ! print*,'QREFvis3d=',QREFvis3d(ig,l,iaer) ! print*,'reffrad=',reffrad(ig,l,iaer) ! endif ! end do ! end do ! end do ! do ig=1,ngrid ! if(tau_col(ig).gt.1.e3)then ! print*,'WARNING: tau_col=',tau_col(ig) ! print*,'at ig=',ig ! print*,'aerosol=',aerosol(ig,:,:) ! print*,'QREFvis3d=',QREFvis3d(ig,:,:) ! print*,'reffrad=',reffrad(ig,:,:) ! endif ! end do end subroutine aeropacity end module aeropacity_mod