SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls, & pq,ccn,tauref,tau,aerosol,reffrad,nueffrad, & QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d) ! to use 'getin' USE ioipsl_getincom IMPLICIT NONE c======================================================================= c subject: c -------- c Computing aerosol optical depth in each gridbox. c c author: F.Forget c ------ c update F. Montmessin (water ice scheme) c and S. Lebonnois (12/06/2003) compatibility dust/ice/chemistry c update J.-B. Madeleine 2008-2009: c - added 3D scattering by aerosols; c - dustopacity transferred from physiq.F to callradite.F, c and renamed into aeropacity.F; c c input: c ----- c ngrid Number of gridpoint of horizontal grid c nlayer Number of layer c nq Number of tracer c zday Date (time since Ls=0, in martian days) c ls Solar longitude (Ls) , radian c pplay,pplev pressure (Pa) in the middle and boundary of each layer c pq Dust mixing ratio (used if tracer =T and active=T). c reffrad(ngrid,nlayer,naerkind) Aerosol effective radius c QREFvis3d(ngridmx,nlayermx,naerkind) \ 3d extinction coefficients c QREFir3d(ngridmx,nlayermx,naerkind) / at reference wavelengths; c omegaREFvis3d(ngridmx,nlayermx,naerkind) \ 3d single scat. albedo c omegaREFir3d(ngridmx,nlayermx,naerkind) / at reference wavelengths; c c output: c ------- c tauref Prescribed mean column optical depth at 700 Pa c tau Column total visible dust optical depth at each point c aerosol aerosol(ig,l,1) is the dust optical c depth in layer l, grid point ig c c======================================================================= #include "dimensions.h" #include "dimphys.h" #include "callkeys.h" #include "comcstfi.h" #include "comgeomfi.h" #include "dimradmars.h" #include "yomaer.h" #include "tracer.h" #include "planete.h" #include "aerkind.h" c----------------------------------------------------------------------- c c Declarations : c -------------- c c Input/Output c ------------ INTEGER ngrid,nlayer,nq REAL ls,zday,expfactor REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer) REAL pq(ngrid,nlayer,nq) REAL tauref(ngrid), tau(ngrid,naerkind) REAL aerosol(ngrid,nlayer,naerkind) REAL dsodust(ngridmx,nlayermx) REAL reffrad(ngrid,nlayer,naerkind) REAL nueffrad(ngrid,nlayer,naerkind) REAL QREFvis3d(ngridmx,nlayermx,naerkind) REAL QREFir3d(ngridmx,nlayermx,naerkind) REAL omegaREFvis3d(ngridmx,nlayermx,naerkind) REAL omegaREFir3d(ngridmx,nlayermx,naerkind) c c Local variables : c ----------------- INTEGER l,ig,iq,i,j INTEGER iaer ! Aerosol index real topdust(ngridmx) real zlsconst, zp real taueq,tauS,tauN c Mean Qext(vis)/Qext(ir) profile real msolsir(nlayermx,naerkind) c Mean Qext(ir)/Qabs(ir) profile real mqextsqabs(nlayermx,naerkind) c Variables used when multiple particle sizes are used c for dust or water ice particles in the radiative transfer c (see callradite.F for more information). REAL taudusttmp(ngridmx)! Temporary dust opacity ! used before scaling REAL taudustvis(ngridmx) ! Dust opacity after scaling REAL taudusttes(ngridmx) ! Dust opacity at IR ref. wav. as ! "seen" by the GCM. REAL taucloudvis(ngridmx)! Cloud opacity at visible ! reference wavelength REAL taucloudtes(ngridmx)! Cloud opacity at infrared ! reference wavelength using ! Qabs instead of Qext ! (direct comparison with TES) REAL qdust(ngridmx,nlayermx) ! True dust mass mixing ratio REAL ccn(ngridmx,nlayermx) ! Cloud condensation nuclei ! (particules kg-1) REAL qtot(ngridmx) ! Dust column (kg m-2) c CCN reduction factor REAL, PARAMETER :: ccn_factor = 4.5 !! comme TESTS_JB // 1. avant c c local saved variables c --------------------- REAL topdust0(ngridmx) SAVE topdust0 c Level under which the dust mixing ratio is held constant c when computing the dust opacity in each layer c (this applies when doubleq and active are true) INTEGER, PARAMETER :: cstdustlevel = 7 LOGICAL firstcall DATA firstcall/.true./ SAVE firstcall ! indexes of water ice and dust tracers: INTEGER,SAVE :: nqdust(nqmx) ! to store the indexes of dust tracers INTEGER,SAVE :: i_ice=0 ! water ice CHARACTER(LEN=20) :: txt ! to temporarly store text CHARACTER(LEN=1) :: txt2 ! to temporarly store text ! indexes of dust scatterers: INTEGER,SAVE :: iaerdust(naerkind) INTEGER,SAVE :: naerdust ! number of dust scatterers tau(1:ngrid,1:naerkind)=0 ! identify tracers IF (firstcall) THEN ! identify scatterers that are dust naerdust=0 DO iaer=1,naerkind txt=name_iaer(iaer) IF (txt(1:4).eq."dust") THEN naerdust=naerdust+1 iaerdust(naerdust)=iaer ENDIF ENDDO ! identify tracers which are dust i=0 DO iq=1,nq txt=noms(iq) IF (txt(1:4).eq."dust") THEN i=i+1 nqdust(i)=iq ENDIF ENDDO IF (water.AND.activice) THEN i_ice=igcm_h2o_ice write(*,*) "aeropacity: i_ice=",i_ice ENDIF c altitude of the top of the aerosol layer (km) at Ls=2.76rad: c in the Viking year scenario DO ig=1,ngrid topdust0(ig)=60. -22.*SIN(lati(ig))**2 END DO c typical profile of solsir and (1-w)^(-1): msolsir(1:nlayer,1:naerkind)=0 mqextsqabs(1:nlayer,1:naerkind)=0 WRITE(*,*) "Typical profiles of solsir and Qext/Qabs(IR):" DO iaer = 1, naerkind ! Loop on aerosol kind WRITE(*,*) "Aerosol # ",iaer DO l=1,nlayer DO ig=1,ngridmx msolsir(l,iaer)=msolsir(l,iaer)+ & QREFvis3d(ig,l,iaer)/ & QREFir3d(ig,l,iaer) mqextsqabs(l,iaer)=mqextsqabs(l,iaer)+ & (1.E0-omegaREFir3d(ig,l,iaer))**(-1) ENDDO msolsir(l,iaer)=msolsir(l,iaer)/REAL(ngridmx) mqextsqabs(l,iaer)=mqextsqabs(l,iaer)/REAL(ngridmx) ENDDO WRITE(*,*) "solsir: ",msolsir(:,iaer) WRITE(*,*) "Qext/Qabs(IR): ",mqextsqabs(:,iaer) ENDDO ! load value of tauvis from callphys.def (if given there, ! otherwise default value read from starfi.nc file will be used) call getin("tauvis",tauvis) firstcall=.false. END IF c Vertical column optical depth at 700.Pa c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF(iaervar.eq.1) THEN do ig=1, ngridmx tauref(ig)=max(tauvis,1.e-9) ! tauvis=cste (set in callphys.def ! or read in starfi end do ELSE IF (iaervar.eq.2) THEN ! << "Viking" Scenario>> tauref(1) = 0.7+.3*cos(ls+80.*pi/180.) ! like seen by VL1 do ig=2,ngrid tauref(ig) = tauref(1) end do ELSE IF (iaervar.eq.3) THEN ! << "MGS" scenario >> taueq= 0.2 +(0.5-0.2) *(cos(0.5*(ls-4.363)))**14 tauS= 0.1 +(0.5-0.1) *(cos(0.5*(ls-4.363)))**14 tauN = 0.1 c if (peri_day.eq.150) then c tauS=0.1 c tauN=0.1 +(0.5-0.1) *(cos(0.5*(ls+pi-4.363)))**14 c taueq= 0.2 +(0.5-0.2) *(cos(0.5*(ls+pi-4.363)))**14 c endif do ig=1,ngrid/2 ! Northern hemisphere tauref(ig)= tauN + & (taueq-tauN)*0.5*(1+tanh((45-lati(ig)*180./pi)*6/60)) end do do ig=ngrid/2+1, ngridmx ! Southern hemisphere tauref(ig)= tauS + & (taueq-tauS)*0.5*(1+tanh((45+lati(ig)*180./pi)*6/60)) end do ELSE IF ((iaervar.eq.4).or. & ((iaervar.ge.24).and.(iaervar.le.26))) & THEN ! << "TES assimilated dust scenarios >> call readtesassim(ngrid,nlayer,zday,pplev,tauref) ELSE IF (iaervar.eq.5) THEN ! << Escalier Scenario>> c tauref(1) = 0.2 c if ((ls.ge.210.*pi/180.).and.(ls.le.330.*pi/180.)) c & tauref(1) = 2.5 tauref(1) = 2.5 if ((ls.ge.30.*pi/180.).and.(ls.le.150.*pi/180.)) & tauref(1) = .2 do ig=2,ngrid tauref(ig) = tauref(1) end do ELSE stop 'problem with iaervar in aeropacity.F' ENDIF c ----------------------------------------------------------------- c Computing the opacity in each layer c ----------------------------------------------------------------- DO iaer = 1, naerkind ! Loop on aerosol kind c -------------------------------------------- aerkind: SELECT CASE (name_iaer(iaer)) c================================================================== CASE("dust_conrath") aerkind ! Typical dust profile c================================================================== c Altitude of the top of the dust layer c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ zlsconst=SIN(ls-2.76) if (iddist.eq.1) then do ig=1,ngrid topdust(ig)=topdustref ! constant dust layer top end do else if (iddist.eq.2) then ! "Viking" scenario do ig=1,ngrid topdust(ig)=topdust0(ig)+18.*zlsconst end do else if(iddist.eq.3) then !"MGS" scenario do ig=1,ngrid topdust(ig)=60.+18.*zlsconst & -(32+18*zlsconst)*sin(lati(ig))**4 & - 8*zlsconst*(sin(lati(ig)))**5 end do endif c Optical depth in each layer : c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if(iddist.ge.1) then expfactor=0. DO l=1,nlayer DO ig=1,ngrid c Typical mixing ratio profile if(pplay(ig,l).gt.700. $ /(988.**(topdust(ig)/70.))) then zp=(700./pplay(ig,l))**(70./topdust(ig)) expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3) else expfactor=1.e-3 endif c Vertical scaling function aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) * & expfactor * & QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer) ENDDO ENDDO else if(iddist.eq.0) then c old dust vertical distribution function (pollack90) DO l=1,nlayer DO ig=1,ngrid zp=700./pplay(ig,l) aerosol(ig,l,1)= tauref(ig)/700. * s (pplev(ig,l)-pplev(ig,l+1)) s *max( exp(.03*(1.-max(zp,1.))) , 1.E-3 ) ENDDO ENDDO end if c================================================================== CASE("dust_doubleq") aerkind! Two-moment scheme for dust c (transport of mass and number mixing ratio) c================================================================== DO l=1,nlayer IF (l.LE.cstdustlevel) THEN c Opacity in the first levels is held constant to c avoid unrealistic values due to constant lifting: DO ig=1,ngrid aerosol(ig,l,iaer) = & ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) / & ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) * & pq(ig,cstdustlevel,igcm_dust_mass) * & ( pplev(ig,l) - pplev(ig,l+1) ) / g ENDDO ELSE DO ig=1,ngrid aerosol(ig,l,iaer) = & ( 0.75 * QREFvis3d(ig,l,iaer) / & ( rho_dust * reffrad(ig,l,iaer) ) ) * & pq(ig,l,igcm_dust_mass) * & ( pplev(ig,l) - pplev(ig,l+1) ) / g ENDDO ENDIF ENDDO c================================================================== CASE("dust_submicron") aerkind ! Small dust population c================================================================== DO l=1,nlayer IF (l.LE.cstdustlevel) THEN c Opacity in the first levels is held constant to c avoid unrealistic values due to constant lifting: DO ig=1,ngrid aerosol(ig,l,iaer) = & ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) / & ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) * & pq(ig,cstdustlevel,igcm_dust_submicron) * & ( pplev(ig,l) - pplev(ig,l+1) ) / g ENDDO ELSE DO ig=1,ngrid aerosol(ig,l,iaer) = & ( 0.75 * QREFvis3d(ig,l,iaer) / & ( rho_dust * reffrad(ig,l,iaer) ) ) * & pq(ig,l,igcm_dust_submicron) * & ( pplev(ig,l) - pplev(ig,l+1) ) / g ENDDO ENDIF ENDDO c================================================================== CASE("h2o_ice") aerkind ! Water ice crystals c================================================================== c 1. Initialization aerosol(1:ngrid,1:nlayer,iaer) = 0. taucloudvis(1:ngrid) = 0. taucloudtes(1:ngrid) = 0. c 2. Opacity calculation DO ig=1, ngrid DO l=1,nlayer aerosol(ig,l,iaer) = max(1E-20, & ( 0.75 * QREFvis3d(ig,l,iaer) / & ( rho_ice * reffrad(ig,l,iaer) ) ) * & pq(ig,l,i_ice) * & ( pplev(ig,l) - pplev(ig,l+1) ) / g & ) taucloudvis(ig) = taucloudvis(ig) + aerosol(ig,l,iaer) taucloudtes(ig) = taucloudtes(ig) + aerosol(ig,l,iaer)* & QREFir3d(ig,l,iaer) / QREFvis3d(ig,l,iaer) * & ( 1.E0 - omegaREFir3d(ig,l,iaer) ) ENDDO ENDDO c 3. Outputs IF (ngrid.NE.1) THEN CALL WRITEDIAGFI(ngridmx,'tauVIS','tauext VIS refwvl', & ' ',2,taucloudvis) CALL WRITEDIAGFI(ngridmx,'tauTES','tauabs IR refwvl', & ' ',2,taucloudtes) IF (callstats) THEN CALL wstats(ngridmx,'tauVIS','tauext VIS refwvl', & ' ',2,taucloudvis) CALL wstats(ngridmx,'tauTES','tauabs IR refwvl', & ' ',2,taucloudtes) ENDIF ELSE c CALL writeg1d(ngrid,1,taucloudtes,'tautes','NU') ENDIF c================================================================== END SELECT aerkind c ----------------------------------- ENDDO ! iaer (loop on aerosol kind) c ----------------------------------------------------------------- c Rescaling each layer to reproduce the choosen (or assimilated) c dust extinction opacity at visible reference wavelength, which c is originally scaled to an equivalent 700Pa pressure surface. c ----------------------------------------------------------------- taudusttmp(1:ngrid)=0. DO iaer=1,naerdust DO l=1,nlayer DO ig=1,ngrid c Scaling factor taudusttmp(ig) = taudusttmp(ig) + & aerosol(ig,l,iaerdust(iaer)) ENDDO ENDDO ENDDO DO iaer=1,naerdust DO l=1,nlayer DO ig=1,ngrid aerosol(ig,l,iaerdust(iaer)) = max(1E-20, & tauref(ig) * & pplev(ig,1) / 700.E0 * & aerosol(ig,l,iaerdust(iaer)) / & taudusttmp(ig) & ) ENDDO ENDDO ENDDO c ----------------------------------------------------------------- c Computing the number of condensation nuclei c ----------------------------------------------------------------- DO iaer = 1, naerkind ! Loop on aerosol kind c -------------------------------------------- aerkind2: SELECT CASE (name_iaer(iaer)) c================================================================== CASE("dust_conrath") aerkind2 ! Typical dust profile c================================================================== DO l=1,nlayer DO ig=1,ngrid ccn(ig,l) = max(aerosol(ig,l,iaer) / & pi / QREFvis3d(ig,l,iaer) * & (1.+nueffrad(ig,l,iaer))**3. / & reffrad(ig,l,iaer)**2. * g / & (pplev(ig,l)-pplev(ig,l+1)),1e-30) ENDDO ENDDO c================================================================== CASE("dust_doubleq") aerkind2!Two-moment scheme for dust c (transport of mass and number mixing ratio) c================================================================== qtot(1:ngridmx) = 0. DO l=1,nlayer DO ig=1,ngrid qdust(ig,l) = pq(ig,l,igcm_dust_mass) * tauref(ig) * & pplev(ig,1) / 700.E0 / taudusttmp(ig) qtot(ig) = qtot(ig) + qdust(ig,l) * & (pplev(ig,l)-pplev(ig,l+1)) / g ccn(ig,l) = max( ( ref_r0 / & reffrad(ig,l,iaer) )**3. * & r3n_q * qdust(ig,l) ,1e-30) ENDDO ENDDO c================================================================== END SELECT aerkind2 c ----------------------------------- ENDDO ! iaer (loop on aerosol kind) c ----------------------------------------------------------------- c ----------------------------------------------------------------- c Reduce number of nuclei ! TEMPORAIRE : r�duction du nombre de nuclei FF 04/200 ! reduction facteur 3 ! ccn(ig,l) = ccn(ig,l) / 27. ! reduction facteur 2 ! ccn(ig,l) = ccn(ig,l) / 8. c ----------------------------------------------------------------- write(*,*) "water_param CCN reduc. fac. ", ccn_factor DO l=1,nlayer DO ig=1,ngrid ccn(ig,l) = ccn(ig,l) / ccn_factor ENDDO ENDDO c ----------------------------------------------------------------- c ----------------------------------------------------------------- c ----------------------------------------------------------------- c Column integrated visible optical depth in each point c ----------------------------------------------------------------- DO iaer=1,naerkind do l=1,nlayer do ig=1,ngrid tau(ig,iaer) = tau(ig,iaer) + aerosol(ig,l,iaer) end do end do ENDDO c ----------------------------------------------------------------- c Density scaled opacity and column opacity output c ----------------------------------------------------------------- dsodust(1:ngrid,1:nlayer) = 0. DO iaer=1,naerdust DO l=1,nlayermx DO ig=1,ngrid dsodust(ig,l) = dsodust(ig,l) + & aerosol(ig,l,iaerdust(iaer)) * g / & (pplev(ig,l) - pplev(ig,l+1)) ENDDO ENDDO IF (ngrid.NE.1) THEN write(txt2,'(i1.1)') iaer call WRITEDIAGFI(ngridmx,'taudust'//txt2, & 'Dust col opacity', & ' ',2,tau(1,iaerdust(iaer))) IF (callstats) THEN CALL wstats(ngridmx,'taudust'//txt2, & 'Dust col opacity', & ' ',2,tau(1,iaerdust(iaer))) ENDIF ENDIF ENDDO IF (ngrid.NE.1) THEN c CALL WRITEDIAGFI(ngridmx,'dsodust','tau*g/dp', c & 'm2.kg-1',3,dsodust) IF (callstats) THEN CALL wstats(ngridmx,'dsodust', & 'tau*g/dp', & 'm2.kg-1',3,dsodust) ENDIF c CALL WRITEDIAGFI(ngridmx,'ccn','Cond. nuclei', c & 'part kg-1',3,ccn) ELSE CALL writeg1d(ngrid,nlayer,dsodust,'dsodust','m2.kg-1') ENDIF c ----------------------------------------------------------------- return end