Index: trunk/LMDZ.PLUTO/libf/phypluto/aeropacity.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/aeropacity.F90 (revision 3957) +++ trunk/LMDZ.PLUTO/libf/phypluto/aeropacity.F90 (revision 3959) @@ -13,5 +13,5 @@ use comcstfi_mod, only: g, pi, mugaz, avocado use geometry_mod, only: latitude - use callkeys_mod, only: kastprof, callmufi + use callkeys_mod, only: kastprof, callmufi, callmuclouds use mp2m_diagnostics implicit none @@ -32,11 +32,11 @@ ! 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 +! 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 ! @@ -48,18 +48,18 @@ !======================================================================= - 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) :: zzlev(ngrid,nlayer) ! Altitude at the layer boundaries. - REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air) - REAL,INTENT(IN) :: pt(ngrid,nlayer) ! mid-layer temperature (K) - REAL,INTENT(OUT) :: dtau_aer(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 + 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) :: zzlev(ngrid,nlayer) ! Altitude at the layer boundaries. + REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! Tracers (.../kg_of_air) + REAL,INTENT(IN) :: pt(ngrid,nlayer) ! Mid-layer temperature (K) + REAL,INTENT(OUT) :: dtau_aer(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) ! Extinction coefficient in the infrared + REAL,INTENT(OUT) :: tau_col(ngrid) ! Column integrated visible optical depth real aerosol0, obs_tau_col_aurora, pm @@ -74,4 +74,5 @@ real m0as(ngrid,nlayer) real m0af(ngrid,nlayer) + real m0ccn(ngrid,nlayer) INTEGER l,ig,iq,iaer,ia @@ -135,6 +136,14 @@ endwhere - ! write(*,*) 'dtau_as :', MINVAL(dtau_aer(:,:,1)), '-', MAXVAL(dtau_aer(:,:,1)) - ! write(*,*) 'dtau_af :', MINVAL(dtau_aer(:,:,2)), '-', MAXVAL(dtau_aer(:,:,2)) + ! Cloud drops + if (callmuclouds) then + m0ccn(:,:) = pq(:,:,micro_indx(5)) * (pplev(:,1:nlayer) - pplev(:,2:nlayer+1)) / g + sig = 0.2 + where ((m0ccn(:,:) >= 1e-8).and.(mp2m_rc_cld(:,:) >= 5e-9)) + dtau_aer(:,:,3) = m0ccn(:,:) * QREFvis3d(:,:,3) * pi * mp2m_rc_cld(:,:)**2 * exp(2*sig**2) + elsewhere + dtau_aer(:,:,3) = 0d0 + endwhere + endif ! end callmuclouds else Index: trunk/LMDZ.PLUTO/libf/phypluto/callcorrk.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/callcorrk.F90 (revision 3957) +++ trunk/LMDZ.PLUTO/libf/phypluto/callcorrk.F90 (revision 3959) @@ -37,5 +37,5 @@ methane,carbox,cooling,nlte,strobel,& ch4fix,vmrch4_proffix,vmrch4fix,& - callmufi,triton + callmufi,callmuclouds,triton use optcv_mod, only: optcv use optci_mod, only: optci @@ -492,4 +492,14 @@ endwhere nueffrad(:,:,2) = exp(sig**2) - 1 + ! Cloud drops + if (callmuclouds) then + sig = 0.2 + where (mp2m_rc_cld(:,:) >= 5e-9) + reffrad(:,:,3) = mp2m_rc_cld(:,:) * exp(5.*sig**2 / 2.) + elsewhere + reffrad(:,:,3) = 0d0 + endwhere + nueffrad(:,:,3) = exp(sig**2) - 1 + endif ! end callmuclouds else @@ -593,6 +603,4 @@ !----------------------------------------------------------------------- - - ! AF24: for now only consider one aerosol (=haze) if (optichaze) then do iaer=1,naerkind @@ -676,5 +684,5 @@ end do ! naerkind - ! Test / Correct for freaky s. s. albedo values. + ! Test / Correct for freaky single scattering albedo values. do iaer=1,naerkind do k=1,L_LEVELS Index: trunk/LMDZ.PLUTO/libf/phypluto/callkeys_mod.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/callkeys_mod.F90 (revision 3957) +++ trunk/LMDZ.PLUTO/libf/phypluto/callkeys_mod.F90 (revision 3959) @@ -217,5 +217,8 @@ !! Variables for aerosol absorption real,save :: Fabs_aers_VI, Fabs_aerf_VI, Fabs_aers_IR, Fabs_aerf_IR -!$OMP THREADPRIVATE(Fabs_aers_VI, Fabs_aerf_VI, Fabs_aers_IR, Fabs_aerf_IR) +!$OMP THREADPRIVATE(Fabs_aers_VI,Fabs_aerf_VI,Fabs_aers_IR,Fabs_aerf_IR) +!! Variables for cloud drop absorption + real,save :: Fabs_cldd_VI, Fabs_cldd_IR +!$OMP THREADPRIVATE(Fabs_cldd_VI,Fabs_cldd_IR) integer,save :: iddist Index: trunk/LMDZ.PLUTO/libf/phypluto/datafile_mod.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/datafile_mod.F90 (revision 3957) +++ trunk/LMDZ.PLUTO/libf/phypluto/datafile_mod.F90 (revision 3959) @@ -22,4 +22,5 @@ character(len=300),save :: aersprop_file character(len=300),save :: aerfprop_file + character(len=300),save :: clddprop_file ! surfdir stores planetary topography, albedo, etc. (surface.nc files) Index: trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90 (revision 3957) +++ trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90 (revision 3959) @@ -14,5 +14,5 @@ use radii_mod, only: radfixed use datafile_mod, only: datadir,hazeprop_file,hazerad_file,hazemmr_file,hazedens_file, & - config_mufi, mugasflux_file, aersprop_file, aerfprop_file + config_mufi, mugasflux_file, aersprop_file, aerfprop_file, clddprop_file use comdiurn_h, only: sinlat, coslat, sinlon, coslon use comgeomfi_h, only: totarea, totarea_planet @@ -801,4 +801,9 @@ if (is_master) write(*,*) trim(rname)//" aerfprop_file = ",trim(aerfprop_file) + if (is_master) write(*,*) "Cloud drop optical properties datafile" + clddprop_file="optprop_rannou_r2-200nm_nu003.dat" ! default file + call getin_p("clddprop_file",clddprop_file) + if (is_master) write(*,*) trim(rname)//" clddprop_file = ",trim(clddprop_file) + if (is_master) write(*,*) "Use haze production from CH4 photolysis or production rate?" call_haze_prod_pCH4=.false. ! default value @@ -879,4 +884,14 @@ call getin_p("Fabs_aerf_IR",Fabs_aerf_IR) if (is_master) write(*,*)" Fabs_aerf_IR = ",Fabs_aerf_IR + + if (is_master) write(*,*) "Cloud drop absorption correction in VI?" + Fabs_cldd_VI=1. ! default value + call getin_p("Fabs_cldd_VI",Fabs_cldd_VI) + if (is_master) write(*,*)" Fabs_cldd_VI = ",Fabs_cldd_VI + + if (is_master) write(*,*) "Cloud drop absorption correction in IR?" + Fabs_cldd_IR=1. ! default value + call getin_p("Fabs_cldd_IR",Fabs_cldd_IR) + if (is_master) write(*,*)" Fabs_cldd_IR = ",Fabs_cldd_IR ! Pluto haze model @@ -1454,7 +1469,10 @@ call abort_physic(rname, 'if microphysics is on, naerkind must be > 1!', 1) endif - ! if ((callmufi).and.(callmuclouds).and..not.(naerkind.gt.2)) then - ! call abort_physic(rname, 'if microphysical clouds are on, naerkind must be > 2!', 1) - ! endif + if ((callmufi).and..not.(callmuclouds).and.(naerkind.gt.2)) then + call abort_physic(rname, 'Warning: here microphysical clouds are on, naerkind must be = 2!', 1) + endif + if ((callmufi).and.(callmuclouds).and..not.(naerkind.gt.2)) then + call abort_physic(rname, 'if microphysical clouds are on, naerkind must be > 2!', 1) + endif if (.not.(callmufi.or.haze).and.(optichaze)) then call abort_physic(rname, 'if microphysics and haze are off, optichaze must be deactivated!', 1) Index: trunk/LMDZ.PLUTO/libf/phypluto/optci.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/optci.F90 (revision 3957) +++ trunk/LMDZ.PLUTO/libf/phypluto/optci.F90 (revision 3959) @@ -15,5 +15,6 @@ igas_CH4, igas_N2 use comcstfi_mod, only: g, r, mugaz - use callkeys_mod, only: kastprof,continuum,graybody,callmufi,Fabs_aers_IR,Fabs_aerf_IR + use callkeys_mod, only: kastprof,continuum,graybody,callmufi,callmuclouds, & + Fabs_aers_IR,Fabs_aerf_IR,Fabs_cldd_IR use recombin_corrk_mod, only: corrk_recombin, gasi_recomb use tpindex_mod, only: tpindex @@ -143,10 +144,13 @@ lkcoef(:,:) = 0.0 - if(callmufi) then + if (callmufi) then Fabs_aer(1) = Fabs_aers_IR Fabs_aer(2) = Fabs_aerf_IR + if (callmuclouds) then + Fabs_aer(3) = Fabs_cldd_IR + endif else Fabs_aer(:) = 1.0 - endif + endif ! end callmufi do K=2,L_LEVELS @@ -186,5 +190,5 @@ do K=2,L_LEVELS TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXIAER(K,NW,IAER) - ! Aerosol absorption correction --> impact on opacity. + ! Aerosol/Cloud absorption correction --> impact on opacity. if (callmufi .and. (TAEROS(K,NW,IAER).gt.0.d0)) then TAEROS(K,NW,IAER) = TAEROS(K,NW,IAER) * & @@ -343,4 +347,5 @@ ! Aerosol absorption correction --> impact on single scattering albedo. if (callmufi) then + ! Spherical aerosols if ((TAEROS(K,NW,1).gt.0.d0) .and. TAEROS(K+1,NW,1).gt.0.d0) then btemp(L,NW) = TAUAEROLK(K,NW,1) / & @@ -353,4 +358,5 @@ btemp(L,NW) = TAUAEROLK(K,NW,1) + TAUAEROLK(K+1,NW,1) endif + ! Fractal aerosols if ((TAEROS(K,NW,2).gt.0.d0) .and. TAEROS(K+1,NW,2).gt.0.d0) then btemp(L,NW) = btemp(L,NW) + & @@ -364,4 +370,18 @@ btemp(L,NW) = btemp(L,NW) + TAUAEROLK(K,NW,2) + TAUAEROLK(K+1,NW,2) endif + ! Cloud absorption correction --> impact on single scattering albedo. + if (callmuclouds) then + if ((TAEROS(K,NW,3).gt.0.d0) .and. TAEROS(K+1,NW,3).gt.0.d0) then + btemp(L,NW) = btemp(L,NW) + & + (TAUAEROLK(K,NW,3) / & + (QSIAER(K,NW,3)/QXIAER(K,NW,3) + & + Fabs_aer(3)*(1.-QSIAER(K,NW,3)/QXIAER(K,NW,3)))) + & + (TAUAEROLK(K+1,NW,3) / & + (QSIAER(K+1,NW,3)/QXIAER(K+1,NW,3) + & + Fabs_aer(3)*(1.-QSIAER(K+1,NW,3)/QXIAER(K+1,NW,3)))) + else + btemp(L,NW) = btemp(L,NW) + TAUAEROLK(K,NW,3) + TAUAEROLK(K+1,NW,3) + endif + endif ! end callmuclouds else btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind)) @@ -375,4 +395,5 @@ ! Aerosol absorption correction --> impact on single scattering albedo. if (callmufi) then + ! Spherical aerosols if (TAEROS(K,NW,1).gt.0.d0) then btemp(L,NW) = TAUAEROLK(K,NW,1) / & @@ -382,4 +403,5 @@ btemp(L,NW) = TAUAEROLK(K,NW,1) endif + ! Fractal aerosols if (TAEROS(K,NW,2).gt.0.d0) then btemp(L,NW) = btemp(L,NW) + & @@ -390,4 +412,15 @@ btemp(L,NW) = btemp(L,NW) + TAUAEROLK(K,NW,2) endif + ! Cloud absorption correction --> impact on single scattering albedo. + if (callmuclouds) then + if (TAEROS(K,NW,3).gt.0.d0) then + btemp(L,NW) = btemp(L,NW) + & + (TAUAEROLK(K,NW,3) / & + (QSIAER(K,NW,3)/QXIAER(K,NW,3) + & + Fabs_aer(3)*(1.-QSIAER(K,NW,3)/QXIAER(K,NW,3)))) + else + btemp(L,NW) = btemp(L,NW) + TAUAEROLK(K,NW,3) + endif + endif ! end callmuclouds else btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) Index: trunk/LMDZ.PLUTO/libf/phypluto/optcv.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/optcv.F90 (revision 3957) +++ trunk/LMDZ.PLUTO/libf/phypluto/optcv.F90 (revision 3959) @@ -15,5 +15,6 @@ igas_CH4, igas_N2 use comcstfi_mod, only: g, r, mugaz - use callkeys_mod, only: kastprof,continuum,graybody,callgasvis,callmufi,Fabs_aers_VI,Fabs_aerf_VI + use callkeys_mod, only: kastprof,continuum,graybody,callgasvis,callmufi,callmuclouds, & + Fabs_aers_VI,Fabs_aerf_VI,Fabs_cldd_VI use recombin_corrk_mod, only: corrk_recombin, gasv_recomb use tpindex_mod, only: tpindex @@ -146,7 +147,10 @@ wbarv_aer(:,:,:) = 1.0 - if(callmufi) then + if (callmufi) then Fabs_aer(1) = Fabs_aers_VI Fabs_aer(2) = Fabs_aerf_VI + if (callmuclouds) then + Fabs_aer(3) = Fabs_cldd_VI + endif else Fabs_aer(:) = 1.0 @@ -185,5 +189,5 @@ do K=5,L_LEVELS TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXVAER(K,NW,IAER) - ! Aerosol absorption correction --> impact on opacity. + ! Aerosol/Cloud absorption correction --> impact on opacity. if (callmufi .and. (TAEROS(K,NW,IAER).gt.0.d0)) then TAEROS(K,NW,IAER) = TAEROS(K,NW,IAER) * & @@ -376,4 +380,19 @@ btemp(L,NW) = btemp(L,NW) + TAUAEROLK(K,NW,2) + TAUAEROLK(K+1,NW,2) endif + + ! Cloud absorption correction --> impact on single scattering albedo. + if (callmuclouds) then + if ((TAEROS(K,NW,3).gt.0.d0) .and. TAEROS(K+1,NW,3).gt.0.d0) then + btemp(L,NW) = btemp(L,NW) + & + (TAUAEROLK(K,NW,3) / & + (QSVAER(K,NW,3)/QXVAER(K,NW,3) + & + Fabs_aer(3)*(1.-QSVAER(K,NW,3)/QXVAER(K,NW,3)))) + & + (TAUAEROLK(K+1,NW,3) / & + (QSVAER(K+1,NW,3)/QXVAER(K+1,NW,3) + & + Fabs_aer(3)*(1.-QSVAER(K+1,NW,3)/QXVAER(K+1,NW,3)))) + else + btemp(L,NW) = btemp(L,NW) + TAUAEROLK(K,NW,3) + TAUAEROLK(K+1,NW,3) + endif + endif ! end callmuclouds else btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind)) @@ -409,4 +428,15 @@ btemp(L,NW) = btemp(L,NW) + TAUAEROLK(K,NW,2) endif + ! Cloud absorption correction --> impact on single scattering albedo. + if (callmuclouds) then + if (TAEROS(K,NW,3).gt.0.d0) then + btemp(L,NW) = btemp(L,NW) + & + (TAUAEROLK(K,NW,3) / & + (QSVAER(K,NW,3)/QXVAER(K,NW,3) + & + Fabs_aer(3)*(1.-QSVAER(K,NW,3)/QXVAER(K,NW,3)))) + else + btemp(L,NW) = btemp(L,NW) + TAUAEROLK(K,NW,3) + endif + endif ! end callmuclouds else btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) @@ -421,21 +451,18 @@ DO NW=1,L_NSPECTV DO L=1,L_NLAYRAD-1 - K = 2*L+1 DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) + DTAUKV(K+1,NW,NG) WBARV(L,nw,ng) = ctemp(L,NW) / DTAUV(L,nw,ng) - END DO ! L vertical loop - ! Last level - - L = L_NLAYRAD - K = 2*L+1 - DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) - - WBARV(L,NW,NG) = ctemp(L,NW) / DTAUV(L,NW,NG) - - END DO ! NW spectral loop - END DO ! NG Gauss loop + ! Last level + !----------- + L = L_NLAYRAD + K = 2*L+1 + + DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) + WBARV(L,NW,NG) = ctemp(L,NW) / DTAUV(L,NW,NG) + END DO ! NW spectral loop + END DO ! NG Gauss loop ! Aerosols extinction optical depths Index: trunk/LMDZ.PLUTO/libf/phypluto/suaer_corrk.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/suaer_corrk.F90 (revision 3957) +++ trunk/LMDZ.PLUTO/libf/phypluto/suaer_corrk.F90 (revision 3959) @@ -11,5 +11,5 @@ use radcommon_h, only: blamv,blami,lamrefir,lamrefvis use datafile_mod, only: datadir, aerdir, & - hazeprop_file, aersprop_file, aerfprop_file + hazeprop_file, aersprop_file, aerfprop_file, clddprop_file ! outputs @@ -18,5 +18,5 @@ use radcommon_h, only: qrefvis,qrefir,omegarefir !,omegarefvis use aerosol_mod, only: iaero_haze - use callkeys_mod, only: tplanet, callmufi + use callkeys_mod, only: tplanet, callmufi, callmuclouds use tracer_h, only: noms @@ -127,7 +127,6 @@ !-------------------------------------------------------------- - ! allocate file_id, as naerkind is a variable + ! allocate file_id, as naerkind is a variable (VIS & IR) allocate(file_id(naerkind,2)) - if (callmufi) then @@ -137,11 +136,21 @@ write(*,*)'Suaer fractal aerosols optical properties, using: ', & TRIM(aerfprop_file) + if (callmuclouds) then + write(*,*)'Suaer cloud drop optical properties, using: ', & + TRIM(clddprop_file) + endif endif ! Visible file_id(1,1)=TRIM(aersprop_file) file_id(2,1)=TRIM(aerfprop_file) + if (callmuclouds) then + file_id(3,1)=TRIM(clddprop_file) + endif ! Infrared file_id(1,2)=file_id(1,1) file_id(2,2)=file_id(2,1) + if (callmuclouds) then + file_id(3,2)=file_id(3,1) + endif do iaer=1,naerkind